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

In This Package:

GtDecayerator.cc
Go to the documentation of this file.
00001 
00002 #include "GtDecayerator.h"
00003 #include "GenDecay/NucUtil.h"
00004 #include "GenDecay/NucState.h"
00005 #include "GenDecay/NucDecay.h"
00006 #include "GenDecay/NucVisitor.h"
00007 #include "GenDecay/Radiation.h"
00008 #include "GenDecay/DecayRates.h"
00009 
00010 #include "GaudiKernel/IRndmGenSvc.h"
00011 
00012 #include "HepMC/GenVertex.h"
00013 #include "HepMC/GenEvent.h"
00014 #include "CLHEP/Vector/LorentzVector.h"
00015 #include "CLHEP/Units/SystemOfUnits.h"
00016 
00017 #include <sstream>
00018 
00019 using namespace std;
00020 using namespace GenDecay;
00021 using namespace more;
00022 using namespace more::phys;
00023 
00024 GtDecayerator::GtDecayerator(const std::string& type,
00025                              const std::string& name, 
00026                              const IInterface* parent)
00027     : GaudiTool(type,name,parent)
00028     , m_rates(0)
00029 {
00030     declareInterface<IHepMCEventMutator>(this);
00031 
00032     declareProperty("ParentNuclide", m_parentNuclide="", 
00033                     "Name of the parent nuclide, eg U-238");
00034     declareProperty("ParentAbundance",m_parentAbundance=0,
00035                     "The parent nuclide abundance");
00036     declareProperty("AbundanceMap",m_abundanceMap,
00037                     "Map nuclide name to its abundance");
00038     declareProperty("SecularEquilibrium",m_secularEquilibrium=true,
00039                     "Set uncorrelated daughter nuclides in secular equilibrium with the parent");
00040     declareProperty("CorrelationTime",m_correlationTime=0,
00041                     "Nuclide with shorter lifetime will be corelated parent");
00042 }
00043 
00044 GtDecayerator::~GtDecayerator()
00045 {
00046 }
00047 
00048 
00049 HepMC::GenParticle* GtDecayerator::make_particle(GenDecay::NucDecay& dk)
00050 {
00051     const Radiation* rad = decay_radiation(dk);
00052     if (!rad) return 0;
00053 
00054     double kine = rad->kineticEnergy();
00055     double mass = rad->mass();
00056     double mom = sqrt((kine+mass)*(kine+mass) - mass*mass);
00057 
00058     debug() << "make_particle: type " << rad->type() << " pid " << rad->pid() << " KE(MeV) " << kine/CLHEP::MeV << " mass(MeV) " << mass/CLHEP::MeV << endreq;
00059     if ( rad->type() == EleCapture ) debug() << "make_particle: this type " << rad->type() << " is electron capture " << endreq;
00060 
00061     // Pick random direction
00062 
00063     // cos(angle) from mean direction
00064     double costh = m_rates->uni()*2.0 - 1.0;
00065     double sinth = sqrt(1-costh*costh);
00066 
00067     // azimuth angle around mean direction
00068     double phi = m_rates->uni()*360.0*CLHEP::degree;
00069     double cosphi = cos(phi);
00070     double sinphi = sin(phi);
00071         
00072     CLHEP::Hep3Vector dir(sinth*cosphi,sinth*sinphi,costh);
00073 
00074     CLHEP::HepLorentzVector lorvec(mom*dir,kine+mass);
00075     HepMC::FourVector fourmom(lorvec);
00076 
00077     // for electron capture, assign HEPEVT status code 3 signifying "a documentation line, defined separately from the event 
00078     // history. This could include the two incoming reacting particles, etc."
00079     // HEPEVT status code 1 is the standard code for particles to be propagated further
00080     // HEPEVT standard: http://cepa.fnal.gov/psm/simulation/mcgen/lund/pythia_manual/pythia6.3/pythia6301/node39.html
00081     int status = 1;
00082     if ( rad->type() == EleCapture && rad->pid() == 0 ) status = 3;
00083     debug() << "make_particle: assigning status " << status << " for this particle " << endreq;
00084 
00085     return new HepMC::GenParticle(fourmom,rad->pid(),status);
00086 }
00087 
00088 HepMC::GenParticle* GtDecayerator::make_particle(GenDecay::NucState& state, int status)
00089 {
00090   // PDG code for Ions (Geant4 version)
00091   // Nuclear codes are given as 10-digit numbers +-100ZZZAAAI.
00092   //For a nucleus consisting of np protons and nn neutrons
00093   // A = np + nn and Z = np.
00094   // I gives the isomer level, with I = 0 corresponding
00095   // to the ground state and I >0 to excitations
00097  
00098     more::phys::nucleus mo_nuc = state.nuc();
00099     int mo_pid = 1000000000 + mo_nuc.n_prot()*10000 + mo_nuc.n_part()*10;  
00100     return new HepMC::GenParticle(HepMC::FourVector(0,0,0,0),mo_pid,status);
00101 }
00102 
00103 StatusCode GtDecayerator::mutate(HepMC::GenEvent& event)
00104 {
00105     NucDecay* decay = 0;
00106     const Radiation* radiation = 0;
00107     do {
00108         decay = m_rates->decay();
00109         if (!decay) {
00110             error() << "Failed to retrieve a decay!" << endreq;
00111             return StatusCode::FAILURE;
00112         }
00113 
00114         // Check that primary radiation is known.  
00115         radiation = decay_radiation(*decay);
00116         if (!radiation) {
00117             warning() << "failed to get radiation from decay of " << *(decay->mother) 
00118                       << " via " << *decay << " skipping" << endreq;
00119         }
00120 
00121     } while (!radiation);
00122 
00123     DecayRates::DecayValues_t decayTimes = m_rates->decayTimes(decay);
00124     HepMC::GenVertex* last_vertex = 0;
00125     NucState* last_daughter = 0;
00126     for (size_t ind = 0; ind < decayTimes.size(); ++ind) {
00127 
00128         NucDecay* dk = decayTimes[ind].first;
00129 
00130         HepMC::GenParticle* particle = make_particle(*dk);
00131         if (!particle) {
00132             warning() << "failed to get radiation from decay of " << *(dk->mother) 
00133                       << " via " << *dk << " with fraction=" << dk->fraction
00134                       << " skipping" << endreq;
00135             continue;
00136         }
00137 
00138         double time = more_to_clhep_time(decayTimes[ind].second);
00139         HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(0,0,0,time));
00140         last_vertex = vertex;
00141 
00142         debug() << "Decay: " << *(dk->mother) 
00143                 << " to " << *(dk->daughter)
00144                 << " via " << *dk
00145                 << " @ " << time/CLHEP::second << " seconds"
00146                 << " (lifetime:" << dk->mother->lifetime() << " seconds)"
00147                 << endreq;
00148 
00149         event.add_vertex(vertex);
00150         if (decay == dk) {      // first in chain
00151             event.set_signal_process_vertex(vertex);
00152         }
00153         
00154         vertex->add_particle_out(particle);
00155 
00156         // Save mother particle as info particle (default status=2)
00157         particle = make_particle(*(dk->mother));
00158         if (particle) {
00159             vertex->add_particle_in(particle);
00160         }
00161 
00162         last_daughter = dk->daughter;
00163     }
00164 
00165     // Save final daughter as outgoing
00166     if (last_vertex && last_daughter) {
00167         HepMC::GenParticle* dp = make_particle(*last_daughter,1);
00168         if (dp) {
00169             last_vertex->add_particle_out(dp);
00170         }
00171     }
00172 
00173     debug() << "Decay of " << *(decay->mother) << " via " << *decay 
00174             << " produced " << event.particles_size() 
00175             << " particles in " << event.vertices_size() << " vertices" 
00176             << endreq;
00177     //    event.print() ;  // djaffe
00178 
00179 
00180 
00181     return StatusCode::SUCCESS;
00182 }
00183 
00184 
00185 StatusCode GtDecayerator::initialize()
00186 {
00187     StatusCode sc;
00188 
00189     IRndmGenSvc *rgs = 0;
00190     if (service("RndmGenSvc",rgs,true).isFailure()) {
00191         fatal() << "Failed to get random service" << endreq;
00192         return StatusCode::FAILURE;        
00193     }
00194 
00195     phys::nucleus nucl;
00196     istringstream iss(m_parentNuclide);
00197     iss >> nucl;
00198     debug() << "Staring with parent nuclide: \"" << nucl << "\"" << endreq;
00199     
00200     NucState* head = get_ground(nucl);
00201     if (!head) {
00202         fatal() << "Failed to get ground state for \"" << nucl << "\"" << endreq;
00203         return StatusCode::FAILURE;
00204     }
00205     chain(head,-1);
00206     std::map<NucDecay*,int> chain_map = getChainMap();
00207 
00208     map<phys::nucleus,double> abundance;
00209     map<string,double>::iterator sdit, sdend = m_abundanceMap.end();
00210     for (sdit=m_abundanceMap.begin(); sdit != sdend; ++sdit) {
00211         istringstream iss(sdit->first);
00212         phys::nucleus nuc;
00213         iss >> nuc;
00214         abundance[nuc] = sdit->second;
00215     }
00216 
00217     if (m_parentAbundance <= 0.0 ) {
00218         m_parentAbundance = abundance[nucl];
00219     }
00220     else {
00221         abundance[nucl] = m_parentAbundance;
00222     }
00223 
00224     m_rates = new DecayRates(abundance, m_secularEquilibrium, m_correlationTime/CLHEP::second, 1e-9/*1 ns*/, rgs);
00225     m_rates -> setMap(chain_map);
00226     try {
00227         m_rates->descend(head);
00228     }
00229     catch (const GaudiException& err) {
00230         error() << "Failed to calculate decay rates of chain" << endreq;
00231         error() << err.message() << endreq;
00232         delete m_rates;
00233         m_rates = 0;
00234         return StatusCode::FAILURE;
00235     }
00236 
00237     const vector<NucState*>& mothers = m_rates->mothers();
00238     info() << "I ("
00239            << this->name() << ") got " 
00240            << mothers.size() << " uncorrelated mothers:\n";
00241 
00242     for (size_t ind=0; ind<mothers.size(); ++ind) {
00243         NucState* mom = mothers[ind];
00244         info() << "\t" << *mom
00245                << " hl=" << mom->halflife() << "\n";
00246     }
00247     info () << "total rate = " << m_rates->totalRate() << endreq;
00248 
00249     if (!mothers.size()) return StatusCode::FAILURE;
00250     return StatusCode::SUCCESS;
00251 }
00252 
00253 StatusCode GtDecayerator::finalize()
00254 {
00255     // Clear out the cache at end of job.  This seems to work around bug #484
00256     NucDecayMap_t& dks = const_cast<NucDecayMap_t&>(get_decays());
00257     dks.clear();
00258     
00259     return StatusCode::SUCCESS;
00260 }
| 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