In This Package:

DecayRates.cc

Go to the documentation of this file.

#include "GenDecay/DecayRates.h"
#include "GenDecay/NucState.h"
#include "GenDecay/NucDecay.h"
#include "GenDecay/NucUtil.h"

#include "GaudiKernel/GaudiException.h"

#include <algorithm>

using namespace GenDecay;
using namespace std;

using namespace more;
using namespace more::phys;

DecayRates::DecayRates(const map<phys::nucleus,double>& abundance,
                       bool secularEquilibrium, double corrTime, double epsTime,
                       IRndmGenSvc *rgs)
    : NucVisitor()
    , m_abundance(abundance)
    , 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() 
{
}

// Return a randomly chosen decay.  
GenDecay::NucDecay* DecayRates::decay() 
{
    // 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);
}

// Visit each possible decay in the chain.  save uncorrelated mothers
// mean decay time of each decay of which the mother has some
// abundance.
bool GenDecay::DecayRates::visit(GenDecay::NucState* mother, GenDecay::NucState* daughter, 
                                 GenDecay::NucDecay* decay) 
{
    // 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;
}

/*
  For the given decay, pick a time for it to occur.  

  If the daughter is unstable with a correlated lifetime decay her and
  recuse.
*/

GenDecay::DecayRates::DecayValues_t GenDecay::DecayRates::decayTimes(NucDecay* decay) 
{ 
    DecayValues_t decayTime;
    get_decay_times(decay,decayTime,0.0);
    return decayTime;
}

void DecayRates::get_decay_times(NucDecay* decay, DecayRates::DecayValues_t& decayTime, double now) 
{
    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
}



NucDecay* DecayRates::decay_state(NucState* nuc)
{
    // 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;
}

---