/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

Public Types | Public Member Functions | Protected Member Functions | Private Member Functions | Private Attributes
GenDecay::DecayRates Class Reference

Visit the chain, calculate mean and total rates given abundances and collect correlated decays. More...

#include <DecayRates.h>

Inheritance diagram for GenDecay::DecayRates:
Inheritance graph
[legend]
Collaboration diagram for GenDecay::DecayRates:
Collaboration graph
[legend]

List of all members.

Public Types

typedef std::map
< more::phys::nucleus, double > 
AbundanceMap_t
typedef std::pair< NucDecay
*, double > 
DecayValuePair_t
typedef std::vector
< DecayValuePair_t
DecayValues_t
typedef std::map< NucState
*, double > 
StateValueMap_t

Public Member Functions

 DecayRates (const AbundanceMap_t &abundance, bool secularEquilibrium, double corrTime, double epsTime, IRndmGenSvc *rgs)
 Create a decay rate with given abundance, whether to use secular equilibrium or not, correlation time (seconds, compared to lifetime not half life) and pointer to random number generator.
 DecayRates (bool secularEquilibrium, double corrTime, double epsTime, IRndmGenSvc *rgs)
 As above, but no abundance map specified.
virtual ~DecayRates ()
void addAbundance (more::phys::nucleus nuc, double abundance)
NucDecaydecay ()
NucDecaydecay_state (NucState *nuc)
const std::vector< NucState * > & mothers () const
double totalRate () const
DecayValues_t decayTimes (NucDecay *decay)
 Access resulting decay times. Times are in seconds.
double uni ()
const std::vector< double > & motherRate ()
 Mostly exposed for debugging purposes.
virtual void descend (NucState *state)
 Start the visitation with the given state.
virtual void preDescend ()
 Hooks called before and after actual descent through the hiearchy.
void setMap (std::map< NucDecay *, int > &map)

Protected Member Functions

bool visit (NucState *mother, NucState *, NucDecay *decay)
 Subclass must override.
virtual void postDescend ()

Private Member Functions

void get_decay_times (NucDecay *decay, DecayValues_t &decayTime, double now)

Private Attributes

AbundanceMap_t m_forcedAbundance
std::map< NucState *, double > m_bookkeeping
bool m_secEq
double m_correlationTime
double m_epsilonTime
double m_secEqRate
StateValueMap_t m_meanDecayTimes
std::vector< NucState * > m_mothers
std::vector< double > m_motherRate
double m_totalDecayRate
Rndm::Numbers m_uni

Detailed Description

Visit the chain, calculate mean and total rates given abundances and collect correlated decays.

Definition at line 31 of file DecayRates.h.


Member Typedef Documentation

typedef std::map<more::phys::nucleus,double> GenDecay::DecayRates::AbundanceMap_t

Definition at line 34 of file DecayRates.h.

Definition at line 35 of file DecayRates.h.

Definition at line 36 of file DecayRates.h.

Definition at line 37 of file DecayRates.h.


Constructor & Destructor Documentation

GenDecay::DecayRates::DecayRates ( const AbundanceMap_t abundance,
bool  secularEquilibrium,
double  corrTime,
double  epsTime,
IRndmGenSvc *  rgs 
)

Create a decay rate with given abundance, whether to use secular equilibrium or not, correlation time (seconds, compared to lifetime not half life) and pointer to random number generator.

epsTime is a small epsilon such that any lifetime less is considered zero.

DecayRates::DecayRates ( bool  secularEquilibrium,
double  corrTime,
double  epsTime,
IRndmGenSvc *  rgs 
)

As above, but no abundance map specified.

Use addAbundance() before calling anything else.

Definition at line 33 of file DecayRates.cc.

    : NucVisitor()
    , m_secEq(secularEquilibrium)
    , m_correlationTime(corrTime)
    , m_epsilonTime(epsTime)
    , m_secEqRate(0.0)
{
    StatusCode sc = m_uni.initialize(rgs, Rndm::Flat(0,1));
    if (sc.isFailure()) {
        throw GaudiException("Failed to initialize uniform random numbers",
                             "GtDecayerator::DecayRates",StatusCode::FAILURE);
    }
}
DecayRates::~DecayRates ( ) [virtual]

Definition at line 49 of file DecayRates.cc.

{
}

Member Function Documentation

void DecayRates::addAbundance ( more::phys::nucleus  nuc,
double  abundance 
)

Definition at line 53 of file DecayRates.cc.

{
    m_forcedAbundance[nuc] = abundance;
}
GenDecay::NucDecay * DecayRates::decay ( )

Definition at line 87 of file DecayRates.cc.

{
    double target = m_uni() * m_totalDecayRate;
    double arrow = 0.0;
    NucState* mother = 0;
    for (size_t ind=0; ind<m_mothers.size(); ++ind) {
        NucState* mom = m_mothers[ind];
        arrow += m_motherRate[ind];
        if (arrow > target) {
            mother = mom;
            break;
        }
    }
    if (!mother) {
        cerr << "DecayRates::decay() failed to find mother to decay out of " << m_mothers.size() 
             << " with target = " << target << " and total = " << m_totalDecayRate 
             << endl;
        throw GaudiException("Failed to find mother to decay, this should not happen",
                             "GtDecayerator::DecayRates",StatusCode::FAILURE);
    }
    return decay_state(mother);
}
NucDecay * DecayRates::decay_state ( NucState nuc)

Definition at line 286 of file DecayRates.cc.

{
    // Now, pick which daughter to follow
    vector<NucDecay*>& decays = nuc->decays();
    if (!decays.size()) return 0;

    double total_fraction = 0, target = m_uni();

    for (size_t ind=0; ind < decays.size(); ++ind) {
        NucDecay* dk = decays[ind];

        // Pick daughter
        total_fraction += dk->fraction;
        if (total_fraction > target) {
            //cerr << "Choose decay " << *dk << " at total_fraction = " << total_fraction
            //     << " and target fraction at " << target << endl;
            return dk;
        }
    }

    cerr << "Failed to choose a decay for " << *nuc
         << " target fraction = " << target << ", total fraction at end = " << total_fraction
         << " Precision error?\n";
    return 0;
}
const std::vector<NucState*>& GenDecay::DecayRates::mothers ( ) const [inline]

Definition at line 66 of file DecayRates.h.

{ return m_mothers; }
double GenDecay::DecayRates::totalRate ( ) const [inline]

Definition at line 67 of file DecayRates.h.

{ return m_totalDecayRate; }
GenDecay::DecayRates::DecayValues_t DecayRates::decayTimes ( NucDecay decay)

Access resulting decay times. Times are in seconds.

Definition at line 235 of file DecayRates.cc.

{ 
    DecayValues_t decayTime;
    get_decay_times(decay,decayTime,0.0);
    return decayTime;
}
double GenDecay::DecayRates::uni ( ) [inline]

Definition at line 72 of file DecayRates.h.

{ return m_uni(); }
const std::vector<double>& GenDecay::DecayRates::motherRate ( ) [inline]

Mostly exposed for debugging purposes.

Definition at line 75 of file DecayRates.h.

{ return m_motherRate; }
bool DecayRates::visit ( NucState ,
NucState ,
NucDecay  
) [protected, virtual]

Subclass must override.

Implements GenDecay::NucVisitor.

Definition at line 112 of file DecayRates.cc.

{
    //cerr << "DecayRates::visit: visiting " 
    //     << mother->nuc() << " -> " << daughter->nuc() << " "
    //     << *decay << endl;
    
    if (decay->fraction <= 0.0 || decay->fraction > 1.0) {
        stringstream err;
        err << "DecayRates::visit: illegal decay fraction: " << decay->fraction
            << " mother: " << mother->nuc()
            << " daughter: " << daughter->nuc();
        throw GaudiException(err.str(),"GtDecayerator::DecayRates",StatusCode::FAILURE);
    }

    // debugging
    if (true) {
        double totfrac = 0.0;
        std::vector<NucDecay*>& dks = mother->decays();
        for (size_t ind=0; ind < dks.size(); ++ind) {
            totfrac += dks[ind]->fraction;
        }
                 
        double diff = fabs(totfrac-1.0);
        if (diff > 0.000001) {
            stringstream err;
            err << "DecayRates::visit: illegal total decay fraction, differs from 1.0 by: " << diff
                << " mother: " << mother->nuc();
            throw GaudiException(err.str(),"GtDecayerator::DecayRates",StatusCode::FAILURE);
        }
    }

    // no moms yet, implies we are at the ultimate mother
    if (!m_mothers.size()) {  

        m_mothers.push_back(mother);
        double abundance = m_forcedAbundance[mother->nuc()];
        if (abundance == 0.0) { 
            stringstream err;
            err << "DecayRates::visit: no abundance set for chain's mother nuclide: " 
                << mother->nuc();
            throw GaudiException(err.str(),"GtDecayerator::DecayRates",StatusCode::FAILURE);
        }
        m_bookkeeping[mother] = abundance/mother->lifetime();
        cerr << "Ultimate mother rate: " << m_bookkeeping[mother] << " for " << *mother << endl;
    }
    else {
        if (mother->lifetime() > m_correlationTime) {
            if (find(m_mothers.begin(),m_mothers.end(),mother) == m_mothers.end()) {
                m_mothers.push_back(mother);
            }
        }
        // bookkeeping handled as a daughter
    }


    // Now the daughter

    // The user can force an abundance for ground states, check if
    // this has been done.
    if (daughter->ground_state()) {
        double abundance = m_forcedAbundance[daughter->nuc()];
        if (abundance > 0.0) {
            m_bookkeeping[daughter] = abundance/daughter->lifetime();
            return true;
        }
    }

    // if this decay is uncorrelated and the user hasn't provided an
    // abundance yet doesn't rely on secular equilibrium then we don't
    // know what to do. 
    if (daughter->lifetime() > m_correlationTime) {
        if (!m_secEq) {
            stringstream err;
            err << "DecayRates::visit: no explicit abundance and no "
                << "secular equilibrium for uncorrelated decay. mother: " 
                << mother->nuc()
                << " daughter: " << daughter->nuc();
            throw GaudiException(err.str(),"GtDecayerator::DecayRates",StatusCode::FAILURE);
        }
    }
    
    double rate = m_bookkeeping[mother];
    if (0.0 == rate) {
        stringstream err;
        err << "DecayRates::visit: mother, why u no have rate?"
            << " mother: " << mother->nuc()
            << " daughter: " << daughter->nuc();
        throw GaudiException(err.str(),"GtDecayerator::DecayRates",StatusCode::FAILURE);
    }

    double more_rate = rate * decay->fraction;
    m_bookkeeping[daughter] += more_rate;
    cerr << "Add " << more_rate << " to " << daughter->nuc() << endl;



    std::cout <<  (mother -> nuc()).name() << std::endl << "\n";

    std::cout << "Decay:" << std::endl;
    std::cout << *mother << "\t" << "Mother Rate:" << rate << "\n";
    std::cout << *daughter << "\n" << "Decay fraction:" << decay -> fraction << "\t" << *decay << "\n";
    std::cout << "Contribute to daughter rate:" << rate * decay -> fraction << std::endl;
    std::cout << "Rate:" << m_bookkeeping[daughter] << std::endl;

    std::cout << "There are " << daughter -> n_origin_decays() << " nucs decay into this daughter." <<std::endl;
    std::cout << "There are " << daughter -> ndecays() << " nucs this daughter can decay into." <<std::endl;
    if ( daughter->ndecays() == 0)
        std::cout << "This daughter is one of the final terminal nuc." << std::endl;
    std::cout << "----------------------------------" << std::endl;    



    return true;
}
void DecayRates::postDescend ( ) [protected, virtual]

Reimplemented from GenDecay::NucVisitor.

Definition at line 58 of file DecayRates.cc.

{
    // we should only be called once, guard against multiple calls.
    // If for some reason this must change, replace this return with a
    // wipe.
    if (m_motherRate.size()) return;

    m_totalDecayRate = 0.0;
    //cerr << "Calculating rates for " << m_mothers.size() << " mothers" << endl;
    for (size_t ind=0; ind<m_mothers.size(); ++ind) {
        NucState* mom = m_mothers[ind];
        double rate = m_bookkeeping[mom];
        
        if (rate == 0.0) {
            stringstream err;
            err << "DecayRates::decay: no rate set for nuclide: " 
                << mom->nuc();
            throw GaudiException(err.str(),
                                 "GtDecayerator::DecayRates",StatusCode::FAILURE);
        }
            
        m_totalDecayRate += rate;
        m_motherRate.push_back(rate);
        m_meanDecayTimes[mom] = 1.0/rate;
        //cerr << ind << ": rate " << rate << " for mother " << *mom << endl;
    }
}
void DecayRates::get_decay_times ( NucDecay decay,
DecayRates::DecayValues_t decayTime,
double  now 
) [private]

Definition at line 242 of file DecayRates.cc.

{
    double meanDecayTime = 0.0;

    if (decayTime.size()) {     // not the mother of the chain

        // check if this decay is considered correlated
        if (decay->mother->lifetime() > m_correlationTime) {
            return;    // uncorelated decay, stop descent
        }

        // got a specific decay, use specific branch lifetime
        meanDecayTime = decay->mother->lifetime() / decay->fraction;
        
        //cerr << "Daughter ";
    }
    else {                      // this is the mother of the chain
        // got a general decay, use abundance based lifetime
        meanDecayTime = m_meanDecayTimes[decay->mother];

        //cerr << "Mother:  ";
    }

    double dt = 0.0;
    if (meanDecayTime > 0.0) {
        dt = (-1.0 * log(m_uni())) * meanDecayTime;
    }
    now += dt;

    // cerr << "now="<<now << " meanDecayTime=" << meanDecayTime 
    //      << " decay=" << *decay << " mother=" << *(decay->mother) 
    //      << " daughter=" << *(decay->daughter)
    //      << endl;

    decayTime.push_back(DecayValuePair_t(decay, now));

    NucDecay* the_decay = decay_state(decay->daughter);    
    if (!the_decay) { return; }

    get_decay_times(the_decay,decayTime,now); // recurse on chosen daughter decay
}
void NucVisitor::descend ( NucState state) [virtual, inherited]

Start the visitation with the given state.

Definition at line 16 of file NucVisitor.cc.

{
    this->preDescend();
    this->_descend(state);
    this->addResidual();
    this->postDescend();
}
virtual void GenDecay::NucVisitor::preDescend ( ) [inline, virtual, inherited]

Hooks called before and after actual descent through the hiearchy.

Subclass may optionally implement.

Definition at line 43 of file NucVisitor.h.

{}
void GenDecay::NucVisitor::setMap ( std::map< NucDecay *, int > &  map) [inline, inherited]

Definition at line 45 of file NucVisitor.h.

{ m_countMap = map;}

Member Data Documentation

Definition at line 91 of file DecayRates.h.

std::map<NucState*,double> GenDecay::DecayRates::m_bookkeeping [private]

Definition at line 92 of file DecayRates.h.

Definition at line 94 of file DecayRates.h.

Definition at line 95 of file DecayRates.h.

Definition at line 96 of file DecayRates.h.

Definition at line 99 of file DecayRates.h.

Definition at line 102 of file DecayRates.h.

std::vector<NucState*> GenDecay::DecayRates::m_mothers [private]

Definition at line 105 of file DecayRates.h.

std::vector<double> GenDecay::DecayRates::m_motherRate [private]

Definition at line 106 of file DecayRates.h.

Definition at line 107 of file DecayRates.h.

Rndm::Numbers GenDecay::DecayRates::m_uni [private]

Definition at line 109 of file DecayRates.h.


The documentation for this class was generated from the following files:
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:21:06 for GenDecay by doxygen 1.7.4