In This Package:

GtDecayerator.cc

Go to the documentation of this file.

#include "GtDecayerator.h"
#include "GenDecay/NucUtil.h"
#include "GenDecay/NucState.h"
#include "GenDecay/NucDecay.h"
#include "GenDecay/NucVisitor.h"
#include "GenDecay/Radiation.h"
#include "GenDecay/DecayRates.h"

#include "GaudiKernel/IRndmGenSvc.h"

#include "HepMC/GenVertex.h"
#include "HepMC/GenEvent.h"
#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Units/SystemOfUnits.h"

#include <sstream>

using namespace std;
using namespace GenDecay;
using namespace more;
using namespace more::phys;

GtDecayerator::GtDecayerator(const std::string& type,
                             const std::string& name, 
                             const IInterface* parent)
    : GaudiTool(type,name,parent)
    , m_rates(0)
{
    declareInterface<IHepMCEventMutator>(this);

    declareProperty("ParentNuclide", m_parentNuclide="", 
                    "Name of the parent nuclide, eg U-238");
    declareProperty("ParentAbundance",m_parentAbundance=0,
                    "The parent nuclide abundance");
    declareProperty("AbundanceMap",m_abundanceMap,
                    "Map nuclide name to its abundance");
    declareProperty("SecularEquilibrium",m_secularEquilibrium=true,
                    "Set uncorrelated daughter nuclides in secular equilibrium with the parent");
    declareProperty("CorrelationTime",m_correlationTime=0,
                    "Nuclide with shorter lifetime will be corelated parent");
}

GtDecayerator::~GtDecayerator()
{
}


HepMC::GenParticle* GtDecayerator::make_particle(GenDecay::NucDecay& dk)
{
    const Radiation* rad = decay_radiation(dk);
    if (!rad) return 0;

    double kine = rad->kineticEnergy();
    double mass = rad->mass();
    double mom = sqrt((kine+mass)*(kine+mass) - mass*mass);

    debug() << "make_particle: type " << rad->type() << " pid " << rad->pid() << " KE(MeV) " << kine/CLHEP::MeV << " mass(MeV) " << mass/CLHEP::MeV << endreq;
    if ( rad->type() == EleCapture ) debug() << "make_particle: this type " << rad->type() << " is electron capture " << endreq;

    // Pick random direction

    // cos(angle) from mean direction
    double costh = m_rates->uni()*2.0 - 1.0;
    double sinth = sqrt(1-costh*costh);

    // azimuth angle around mean direction
    double phi = m_rates->uni()*360.0*CLHEP::degree;
    double cosphi = cos(phi);
    double sinphi = sin(phi);
        
    CLHEP::Hep3Vector dir(sinth*cosphi,sinth*sinphi,costh);

    CLHEP::HepLorentzVector lorvec(mom*dir,kine+mass);
    HepMC::FourVector fourmom(lorvec);

    // for electron capture, assign HEPEVT status code 3 signifying "a documentation line, defined separately from the event 
    // history. This could include the two incoming reacting particles, etc."
    // HEPEVT status code 1 is the standard code for particles to be propagated further
    // HEPEVT standard: http://cepa.fnal.gov/psm/simulation/mcgen/lund/pythia_manual/pythia6.3/pythia6301/node39.html
    int status = 1;
    if ( rad->type() == EleCapture && rad->pid() == 0 ) status = 3;
    debug() << "make_particle: assigning status " << status << " for this particle " << endreq;

    return new HepMC::GenParticle(fourmom,rad->pid(),status);
}

HepMC::GenParticle* GtDecayerator::make_particle(GenDecay::NucState& state, int status)
{
  // PDG code for Ions (Geant4 version)
  // Nuclear codes are given as 10-digit numbers +-100ZZZAAAI.
  //For a nucleus consisting of np protons and nn neutrons
  // A = np + nn and Z = np.
  // I gives the isomer level, with I = 0 corresponding
  // to the ground state and I >0 to excitations
 
    more::phys::nucleus mo_nuc = state.nuc();
    int mo_pid = 1000000000 + mo_nuc.n_prot()*10000 + mo_nuc.n_part()*10;  
    return new HepMC::GenParticle(HepMC::FourVector(0,0,0,0),mo_pid,status);
}

StatusCode GtDecayerator::mutate(HepMC::GenEvent& event)
{
    NucDecay* decay = 0;
    const Radiation* radiation = 0;
    do {
        decay = m_rates->decay();
        if (!decay) {
            error() << "Failed to retrieve a decay!" << endreq;
            return StatusCode::FAILURE;
        }

        // Check that primary radiation is known.  
        radiation = decay_radiation(*decay);
        if (!radiation) {
            warning() << "failed to get radiation from decay of " << *(decay->mother) 
                      << " via " << *decay << " skipping" << endreq;
        }

    } while (!radiation);

    DecayRates::DecayValues_t decayTimes = m_rates->decayTimes(decay);
    HepMC::GenVertex* last_vertex = 0;
    NucState* last_daughter = 0;
    for (size_t ind = 0; ind < decayTimes.size(); ++ind) {

        NucDecay* dk = decayTimes[ind].first;

        HepMC::GenParticle* particle = make_particle(*dk);
        if (!particle) {
            warning() << "failed to get radiation from decay of " << *(dk->mother) 
                      << " via " << *dk << " with fraction=" << dk->fraction
                      << " skipping" << endreq;
            continue;
        }

        double time = more_to_clhep_time(decayTimes[ind].second);
        HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(0,0,0,time));
        last_vertex = vertex;

        debug() << "Decay: " << *(dk->mother) 
                << " to " << *(dk->daughter)
                << " via " << *dk
                << " @ " << time/CLHEP::second << " seconds"
                << " (lifetime:" << dk->mother->lifetime() << " seconds)"
                << endreq;

        event.add_vertex(vertex);
        if (decay == dk) {      // first in chain
            event.set_signal_process_vertex(vertex);
        }
        
        vertex->add_particle_out(particle);

        // Save mother particle as info particle (default status=2)
        particle = make_particle(*(dk->mother));
        if (particle) {
            vertex->add_particle_in(particle);
        }

        last_daughter = dk->daughter;
    }

    // Save final daughter as outgoing
    if (last_vertex && last_daughter) {
        HepMC::GenParticle* dp = make_particle(*last_daughter,1);
        if (dp) {
            last_vertex->add_particle_out(dp);
        }
    }

    debug() << "Decay of " << *(decay->mother) << " via " << *decay 
            << " produced " << event.particles_size() 
            << " particles in " << event.vertices_size() << " vertices" 
            << endreq;
    //    event.print() ;  // djaffe



    return StatusCode::SUCCESS;
}


StatusCode GtDecayerator::initialize()
{
    StatusCode sc;

    IRndmGenSvc *rgs = 0;
    if (service("RndmGenSvc",rgs,true).isFailure()) {
        fatal() << "Failed to get random service" << endreq;
        return StatusCode::FAILURE;        
    }

    phys::nucleus nucl;
    istringstream iss(m_parentNuclide);
    iss >> nucl;
    debug() << "Staring with parent nuclide: \"" << nucl << "\"" << endreq;
    
    NucState* head = get_ground(nucl);
    if (!head) {
        fatal() << "Failed to get ground state for \"" << nucl << "\"" << endreq;
        return StatusCode::FAILURE;
    }
    chain(head,-1);

    map<phys::nucleus,double> abundance;
    map<string,double>::iterator sdit, sdend = m_abundanceMap.end();
    for (sdit=m_abundanceMap.begin(); sdit != sdend; ++sdit) {
        istringstream iss(sdit->first);
        phys::nucleus nuc;
        iss >> nuc;
        abundance[nuc] = sdit->second;
    }

    if (m_parentAbundance <= 0.0 ) {
        m_parentAbundance = abundance[nucl];
    }
    else {
        abundance[nucl] = m_parentAbundance;
    }

    m_rates = new DecayRates(abundance, m_secularEquilibrium, m_correlationTime/CLHEP::second, 1e-9/*1 ns*/, rgs);
    try {
        m_rates->descend(head);
    }
    catch (const GaudiException& err) {
        error() << "Failed to calculate decay rates of chain" << endreq;
        error() << err.message() << endreq;
        delete m_rates;
        m_rates = 0;
        return StatusCode::FAILURE;
    }

    const vector<NucState*>& mothers = m_rates->mothers();
    info() << "Got " << mothers.size() << " uncorrelated mothers:\n";
    double totalRate = 0.0;
    for (size_t ind=0; ind<mothers.size(); ++ind) {
        NucState* mom = mothers[ind];
        info() << "\t" << *mom
               << " hl=" << mom->halflife() << "\n";
    }
    info () << "total rate = " << m_rates->totalRate() << endreq;

    if (!mothers.size()) return StatusCode::FAILURE;
    return StatusCode::SUCCESS;
}

StatusCode GtDecayerator::finalize()
{
    // Clear out the cache at end of job.  This seems to work around bug #484
    NucDecayMap_t& dks = const_cast<NucDecayMap_t&>(get_decays());
    dks.clear();
    
    return StatusCode::SUCCESS;
}

---