In This Package:

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:

[legend]Collaboration diagram for GenDecay::DecayRates:[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.
virtual	~DecayRates ()
NucDecay *	decay ()
NucDecay *	decay_state (NucState *nuc)
bool	visit (NucState *mother, NucState *, NucDecay *decay)
	Subclass overrides.
const std::vector< NucState * > &	mothers () const
double	totalRate () const
DecayValues_t	decayTimes (NucDecay *decay)
	Access resulting decay times. Times are in seconds.
double	uni ()
void	descend (NucState *state)
Private Member Functions
void	get_decay_times (NucDecay *decay, DecayValues_t &decayTime, double now)
Private Attributes
AbundanceMap_t	m_abundance
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_motherBranching
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.

typedef std::pair<NucDecay*,double> GenDecay::DecayRates::DecayValuePair_t

Definition at line 35 of file DecayRates.h.

typedef std::vector<DecayValuePair_t> GenDecay::DecayRates::DecayValues_t

Definition at line 36 of file DecayRates.h.

typedef std::map<NucState*,double> GenDecay::DecayRates::StateValueMap_t

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

(

)

[virtual]

Definition at line 33 of file DecayRates.cc.

{
}

---

Member Function Documentation

GenDecay::NucDecay * DecayRates::decay

(

)

Definition at line 38 of file DecayRates.cc.

{
    // Initially the m_meanDecayTimes for mothers are assuming they
    // are decayed in isolation.  Since they all produces decays
    // together with total rate equal to the sum of their individual
    // ones we must use their individual rates to choose a decay and
    // replace "their" mean decay time with the total one so that
    // their chosen times are set correctly.  W/out this one gets N
    // times too much time covered where N is ~number of uncorrelated
    // branches (=N for full SE).  More in trac ticket #315.
    if (!m_motherBranching.size()) {
        for (size_t ind=0; ind<m_mothers.size(); ++ind) {
            NucState* mom = m_mothers[ind];
            double rate = 1.0/m_meanDecayTimes[mom];
            m_motherBranching.push_back(rate);
            m_meanDecayTimes[mom] = 1.0/m_totalDecayRate;
        }
    }

    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_motherBranching[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 212 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;
}

bool DecayRates::visit	(	GenDecay::NucState *	mother,
		GenDecay::NucState *	daughter,
		GenDecay::NucDecay *	decay
	)			[virtual]

Subclass overrides.

Reimplemented from GenDecay::NucVisitor.

Definition at line 81 of file DecayRates.cc.

{
    // don't bother visiting a decay that has a mother already seen
    if (find(m_mothers.begin(),m_mothers.end(),mother) != m_mothers.end()) 
        return true;

    //cerr << "DecayRates::visit: visiting " 
    //     << mother->nuc() << " -> " << daughter->nuc() << " "
    //     << *decay << endl;

    // First time through we have mother of whole chain so can take
    // her abundance as the one to use for secular equilibrium
    if (!m_mothers.size()) {  
        // mother of chain is of course an uncorrelated mother.
        m_mothers.push_back(mother);
        m_totalDecayRate = 0.0;
        double abundance = m_abundance[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);
        }
        double rate = abundance / mother->lifetime();
        if (m_secEq) m_secEqRate = rate;
        m_meanDecayTimes[mother] = 1.0/rate;
        m_totalDecayRate += rate;
        return true;
    }

    // For decays that are considered instantaneous, set their decay
    // time to exactly zero.  (assumes epsilon < correlation time)
    if (mother->lifetime() < m_epsilonTime) {
        m_meanDecayTimes[mother] = 0;
        return true;
    }

    // What is left is either a "fast" correlated decay or a "slow"
    // uncorrelated one.  In either case we find a decay rate based on
    // an abundance.

    double meanDecayTime = 0;
    double abundance = m_abundance[mother->nuc()];

    if (abundance > 0.0) {
        meanDecayTime = mother->lifetime() / abundance;
    }
    else if (m_secEq) {
        meanDecayTime = 1/m_secEqRate;
    }
    else {
        cerr << "Warning: No abundance for " << mother->nuc()
             << " via " << *decay
             << ", and no seqular equlibrium, this decay won't be allowed to happen\n";
        return true;
    }

    // Check if mother is correlated or not
    if (mother->lifetime() > m_correlationTime) {
        if (find(m_mothers.begin(),m_mothers.end(),mother) == m_mothers.end()) {
            m_mothers.push_back(mother);
            m_totalDecayRate += 1.0/meanDecayTime;
        }
    }

    // All decays from this mother have a mean decay time which is the
    // mother's lifetime per her abundance.
    m_meanDecayTimes[mother] = meanDecayTime;

    return true;
}

const std::vector<NucState*>& GenDecay::DecayRates::mothers

(

)

const [inline]

Definition at line 59 of file DecayRates.h.

{ return m_mothers; }

double GenDecay::DecayRates::totalRate

(

)

const [inline]

Definition at line 60 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 161 of file DecayRates.cc.

{ 
    DecayValues_t decayTime;
    get_decay_times(decay,decayTime,0.0);
    return decayTime;
}

double GenDecay::DecayRates::uni

(

)

[inline]

Definition at line 65 of file DecayRates.h.

{ return m_uni(); }

void DecayRates::get_decay_times	(	NucDecay *	decay,
		DecayValues_t &	decayTime,
		double	now
	)			[private]

Definition at line 168 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

)

[inherited]

Definition at line 16 of file NucVisitor.cc.

{
    if (m_countMap[state]) return;
    m_countMap[state] += 1;

    vector<NucDecay*> &decays = state->decays();
    for (size_t ind=0; ind < decays.size(); ++ind) {

        if (decays[ind]->fraction < m_minBranchFraction) continue;

        NucDecay* decay = decays[ind];
        NucState* daughter = decay->daughter;
        if (this->visit(state,daughter,decay)) {
            descend(daughter);
        }
    }
}

---

Member Data Documentation

AbundanceMap_t GenDecay::DecayRates::m_abundance [private]

Definition at line 70 of file DecayRates.h.

bool GenDecay::DecayRates::m_secEq [private]

Definition at line 71 of file DecayRates.h.

double GenDecay::DecayRates::m_correlationTime [private]

Definition at line 72 of file DecayRates.h.

double GenDecay::DecayRates::m_epsilonTime [private]

Definition at line 73 of file DecayRates.h.

double GenDecay::DecayRates::m_secEqRate [private]

Definition at line 76 of file DecayRates.h.

StateValueMap_t GenDecay::DecayRates::m_meanDecayTimes [private]

Definition at line 79 of file DecayRates.h.

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

Definition at line 82 of file DecayRates.h.

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

Definition at line 83 of file DecayRates.h.

double GenDecay::DecayRates::m_totalDecayRate [private]

Definition at line 84 of file DecayRates.h.

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

Definition at line 86 of file DecayRates.h.

---

The documentation for this class was generated from the following files:

• /home/dayabay/installs/trunk_2011_04_11_opt/NuWa-trunk/dybgaudi/Simulation/GenDecay/GenDecay/DecayRates.h
• /home/dayabay/installs/trunk_2011_04_11_opt/NuWa-trunk/dybgaudi/Simulation/GenDecay/src/lib/DecayRates.cc

---