source: trunk/source/processes/hadronic/models/radioactive_decay/src/G4RadioactiveDecay.cc @ 1340

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago

update ti head
File size: 53.8 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
27//
28// MODULE:              G4RadioactiveDecay.cc
29//
30// Author:              F Lei & P R Truscott
31// Organisation:        DERA UK
32// Customer:            ESA/ESTEC, NOORDWIJK
33// Contract:            12115/96/JG/NL Work Order No. 3
34//
35// Documentation avaialable at http://www.space.dera.gov.uk/space_env/rdm.html
36//   These include:
37//       User Requirement Document (URD)
38//       Software Specification Documents (SSD)
39//       Software User Manual (SUM)
40//       Technical Note (TN) on the physics and algorithms
41//
42//    The test and example programs are not included in the public release of
43//    G4 but they can be downloaded from
44//      http://www.space.qinetiq.com/space_env/rdm.html
45//
46// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
47//
48// CHANGE HISTORY
49// --------------
50// 16 February 2006, V.Ivanchenko fix problem in IsApplicable connected with
51//            8.0 particle design
52// 18 October 2002, F. Lei
53//            in the case of beta decay, added a check of the end-energy
54//            to ensure it is > 0.
55//            ENSDF occationally have beta decay entries with zero energies
56//
57// 27 Sepetember 2001, F. Lei
58//            verboselevel(0) used in constructor
59//
60// 01 November 2000, F.Lei
61//            added " ee = e0 +1. ;" as line 763
62//            tagged as "radiative_decay-V02-00-02"             
63// 28 October 2000, F Lei
64//            added fast beta decay mode. Many files have been changed.
65//            tagged as "radiative_decay-V02-00-01"
66//
67// 25 October 2000, F Lei, DERA UK
68//            1) line 1185 added 'const' to work with tag "Track-V02-00-00"
69//            tagged as "radiative_decay-V02-00-00"
70// 14 April 2000, F Lei, DERA UK
71// 0.b.4 release. Changes are:
72//            1) Use PhotonEvaporation instead of DiscreteGammaDeexcitation
73//            2) VR: Significant efficiency improvement
74//
75// 29 February 2000, P R Truscott, DERA UK
76// 0.b.3 release.
77//
78// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
79///////////////////////////////////////////////////////////////////////////////
80//
81#include "G4RadioactiveDecay.hh"
82#include "G4RadioactiveDecaymessenger.hh"
83
84#include "G4DynamicParticle.hh"
85#include "G4DecayProducts.hh"
86#include "G4DecayTable.hh"
87#include "G4PhysicsLogVector.hh"
88#include "G4ParticleChangeForRadDecay.hh"
89#include "G4ITDecayChannel.hh"
90#include "G4BetaMinusDecayChannel.hh"
91#include "G4BetaPlusDecayChannel.hh"
92#include "G4KshellECDecayChannel.hh"
93#include "G4LshellECDecayChannel.hh"
94#include "G4MshellECDecayChannel.hh"
95#include "G4AlphaDecayChannel.hh"
96#include "G4VDecayChannel.hh"
97#include "G4RadioactiveDecayMode.hh"
98#include "G4Ions.hh"
99#include "G4IonTable.hh"
100#include "G4RIsotopeTable.hh"
101#include "G4BetaFermiFunction.hh"
102#include "Randomize.hh"
103#include "G4LogicalVolumeStore.hh"
104#include "G4NuclearLevelManager.hh"
105#include "G4NuclearLevelStore.hh"
106
107#include "G4HadTmpUtil.hh"
108
109#include <vector>
110#include <sstream>
111#include <algorithm>
112#include <fstream>
113
114using namespace CLHEP;
115
116const G4double   G4RadioactiveDecay::levelTolerance =2.0*keV;
117
118///////////////////////////////////////////////////////////////////////////////
119//
120//
121// Constructor
122//
123G4RadioactiveDecay::G4RadioactiveDecay
124  (const G4String& processName)
125  :G4VRestDiscreteProcess(processName, fDecay), HighestBinValue(10.0),
126   LowestBinValue(1.0e-3), TotBin(200), verboseLevel(0)
127{
128#ifdef G4VERBOSE
129  if (GetVerboseLevel()>1) {
130    G4cout <<"G4RadioactiveDecay constructor    Name: ";
131    G4cout <<processName << G4endl;   }
132#endif
133
134  theRadioactiveDecaymessenger = new G4RadioactiveDecaymessenger(this);
135  theIsotopeTable              = new G4RIsotopeTable();
136  aPhysicsTable                = 0;
137  pParticleChange              = &fParticleChangeForRadDecay;
138 
139  //
140  // Now register the Isotopetable with G4IonTable.
141  //
142  G4IonTable *theIonTable =
143    (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
144  G4VIsotopeTable *aVirtualTable = theIsotopeTable;
145  theIonTable->RegisterIsotopeTable(aVirtualTable);
146  //
147  //
148  // Reset the contents of the list of nuclei for which decay scheme data
149  // have been loaded.
150  //
151  LoadedNuclei.clear();
152  //
153  //
154  // Apply default values.
155  //
156  NSourceBin  = 1;
157  SBin[0]     = 0.* s;
158  SBin[1]     = 1e10 * s;
159  SProfile[0] = 1.;
160  SProfile[1] = 1.;
161  NDecayBin   = 1;
162  DBin[0]     = 9.9e9 * s ;
163  DBin[1]     = 1e10 * s;
164  DProfile[0] = 1.;
165  DProfile[1] = 0.;
166  NSplit      = 1;
167  AnalogueMC  = true ;
168  FBeta       = false ;
169  BRBias      = true ;
170  applyICM    = true ;
171  applyARM    = true ;
172  halflifethreshold = 0.*second;
173  //
174  // RDM applies to xall logical volumes as default
175  SelectAllVolumes();
176}
177////////////////////////////////////////////////////////////////////////////////
178//
179//
180// Destructor
181//
182G4RadioactiveDecay::~G4RadioactiveDecay()
183{
184  if (aPhysicsTable != 0)
185    {
186      aPhysicsTable->clearAndDestroy();
187      delete aPhysicsTable;
188    }
189  //  delete theIsotopeTable;
190  delete theRadioactiveDecaymessenger;
191}
192
193///////////////////////////////////////////////////////////////////////////////
194//
195//
196// IsApplicable
197//
198G4bool G4RadioactiveDecay::IsApplicable( const G4ParticleDefinition &
199  aParticle)
200{
201  //
202  //
203  // All particles, other than G4Ions, are rejected by default.
204  //
205  if (((const G4Ions*)(&aParticle))->GetExcitationEnergy() > 0.) {return true;}
206  if (aParticle.GetParticleName() == "GenericIon")      {return true;}
207  else if (!(aParticle.GetParticleType() == "nucleus") || aParticle.GetPDGLifeTime() < 0. ) {return false;}
208  //
209  //
210  // Determine whether the nuclide falls into the correct A and Z range.
211  //
212  G4int A = ((const G4Ions*) (&aParticle))->GetAtomicMass();
213  G4int Z = ((const G4Ions*) (&aParticle))->GetAtomicNumber();
214  if( A> theNucleusLimits.GetAMax() || A< theNucleusLimits.GetAMin())
215    {return false;}
216  else if( Z> theNucleusLimits.GetZMax() || Z< theNucleusLimits.GetZMin())
217    {return false;}
218  return true;
219}
220///////////////////////////////////////////////////////////////////////////////
221//
222//
223// IsLoaded
224//
225G4bool G4RadioactiveDecay::IsLoaded(const G4ParticleDefinition &
226  aParticle)
227{
228  //
229  //
230  // Check whether the radioactive decay data on the ion have already been
231  // loaded.
232  //
233  return std::binary_search(LoadedNuclei.begin(),
234                       LoadedNuclei.end(),
235                       aParticle.GetParticleName());
236}
237////////////////////////////////////////////////////////////////////////////////
238//
239//
240// SelectAVolume
241//
242void G4RadioactiveDecay::SelectAVolume(const G4String aVolume)
243{
244 
245  G4LogicalVolumeStore *theLogicalVolumes;
246  G4LogicalVolume *volume;
247  theLogicalVolumes=G4LogicalVolumeStore::GetInstance();
248  for (size_t i = 0; i < theLogicalVolumes->size(); i++){
249    volume=(*theLogicalVolumes)[i];
250    if (volume->GetName() == aVolume) {
251      ValidVolumes.push_back(aVolume);
252      std::sort(ValidVolumes.begin(), ValidVolumes.end());
253      // sort need for performing binary_search
254#ifdef G4VERBOSE
255      if (GetVerboseLevel()>0)
256        G4cout << " RDM Applies to : " << aVolume << G4endl; 
257#endif
258    }else if(i ==  theLogicalVolumes->size())
259      {
260        G4cerr << "SelectAVolume: "<<aVolume << " is not a valid logical volume name"<< G4endl; 
261      }
262  }
263}
264////////////////////////////////////////////////////////////////////////////////
265//
266//
267// DeSelectAVolume
268//
269void G4RadioactiveDecay::DeselectAVolume(const G4String aVolume)
270{
271  G4LogicalVolumeStore *theLogicalVolumes;
272  G4LogicalVolume *volume;
273  theLogicalVolumes=G4LogicalVolumeStore::GetInstance();
274  for (size_t i = 0; i < theLogicalVolumes->size(); i++){
275    volume=(*theLogicalVolumes)[i];
276    if (volume->GetName() == aVolume) {
277      std::vector<G4String>::iterator location;
278      location = std::find(ValidVolumes.begin(),ValidVolumes.end(),aVolume);
279      if (location != ValidVolumes.end()) {
280        ValidVolumes.erase(location);
281        std::sort(ValidVolumes.begin(), ValidVolumes.end());
282      }else{
283        G4cerr << " DeselectVolume:" << aVolume << " is not in the list"<< G4endl; 
284      }   
285#ifdef G4VERBOSE
286      if (GetVerboseLevel()>0)
287        G4cout << " DeselectVolume: " << aVolume << " is removed from list"<<G4endl; 
288#endif
289    }else if(i ==  theLogicalVolumes->size())
290      {
291        G4cerr << " DeselectVolume:" << aVolume << "is not a valid logical volume name"<< G4endl; 
292      }
293  }
294}
295////////////////////////////////////////////////////////////////////////////////
296//
297//
298// SelectAllVolumes
299//
300void G4RadioactiveDecay::SelectAllVolumes() 
301{
302  G4LogicalVolumeStore *theLogicalVolumes;
303  G4LogicalVolume *volume;
304  theLogicalVolumes=G4LogicalVolumeStore::GetInstance();
305  ValidVolumes.clear();
306#ifdef G4VERBOSE
307  if (GetVerboseLevel()>0)
308    G4cout << " RDM Applies to all Volumes"  << G4endl;
309#endif
310  for (size_t i = 0; i < theLogicalVolumes->size(); i++){
311    volume=(*theLogicalVolumes)[i];
312    ValidVolumes.push_back(volume->GetName());   
313#ifdef G4VERBOSE
314    if (GetVerboseLevel()>0)
315      G4cout << "         RDM Applies to Volume "  << volume->GetName() << G4endl;
316#endif
317  }
318  std::sort(ValidVolumes.begin(), ValidVolumes.end());
319  // sort needed in order to allow binary_search
320}
321////////////////////////////////////////////////////////////////////////////////
322//
323//
324// DeSelectAllVolumes
325//
326void G4RadioactiveDecay::DeselectAllVolumes() 
327{
328  ValidVolumes.clear();
329#ifdef G4VERBOSE
330  if (GetVerboseLevel()>0)
331    G4cout << " RDM removed from all volumes" << G4endl; 
332#endif
333 
334}
335
336///////////////////////////////////////////////////////////////////////////////
337//
338//
339// IsRateTableReady
340//
341G4bool G4RadioactiveDecay::IsRateTableReady(const G4ParticleDefinition &
342  aParticle)
343{
344  //
345  //
346  // Check whether the radioactive decay rates table on the ion has already
347  // been calculated.
348  //
349  G4String aParticleName = aParticle.GetParticleName();
350  for (size_t i = 0; i < theDecayRateTableVector.size(); i++)
351    {
352      if (theDecayRateTableVector[i].GetIonName() == aParticleName)
353        return true ;
354    }
355  return false;
356}
357////////////////////////////////////////////////////////////////////////////////
358//
359//
360// GetDecayRateTable
361//
362// retrieve the decayratetable for the specified aParticle
363//
364void G4RadioactiveDecay::GetDecayRateTable(const G4ParticleDefinition &
365  aParticle)
366{
367
368  G4String aParticleName = aParticle.GetParticleName();
369
370  for (size_t i = 0; i < theDecayRateTableVector.size(); i++)
371    {
372      if (theDecayRateTableVector[i].GetIonName() == aParticleName)
373        {
374          theDecayRateVector = theDecayRateTableVector[i].GetItsRates() ;
375        }
376    }
377#ifdef G4VERBOSE
378        if (GetVerboseLevel()>0)
379          {
380            G4cout <<"The DecayRate Table for "
381                   << aParticleName << " is selected." <<  G4endl;
382          }
383#endif
384}
385////////////////////////////////////////////////////////////////////////////////
386//
387//
388// GetTaoTime
389//
390// to perform the convilution of the sourcetimeprofile function with the
391// decay constants in the decay chain.
392//
393G4double G4RadioactiveDecay::GetTaoTime(G4double t, G4double tao)
394{
395  G4double taotime =0.;
396  G4int nbin;
397  if ( t > SBin[NSourceBin]) {
398    nbin  = NSourceBin;}
399  else {
400    nbin = 0;
401    while (t > SBin[nbin]) nbin++;
402    nbin--;}
403  if (nbin > 0) { 
404    for (G4int i = 0; i < nbin; i++) 
405      {
406        taotime += SProfile[i] * (std::exp(-(t-SBin[i+1])/tao)-std::exp(-(t-SBin[i])/tao));
407      }
408  }
409  taotime +=  SProfile[nbin] * (1-std::exp(-(t-SBin[nbin])/tao));
410#ifdef G4VERBOSE
411  if (GetVerboseLevel()>1)
412    {G4cout <<" Tao time: " <<taotime <<G4endl;}
413#endif
414  return  taotime ;
415}
416 
417////////////////////////////////////////////////////////////////////////////////
418//
419//
420// GetDecayTime
421//
422// Randomly select a decay time for the decay process, following the supplied
423// decay time bias scheme.
424//
425G4double G4RadioactiveDecay::GetDecayTime()
426{
427  G4double decaytime = 0.;
428  G4double rand = G4UniformRand();
429  G4int i = 0;
430  while ( DProfile[i] < rand) i++;
431  rand = G4UniformRand();
432  decaytime = DBin[i] + rand*(DBin[i+1]-DBin[i]);
433#ifdef G4VERBOSE
434  if (GetVerboseLevel()>1)
435    {G4cout <<" Decay time: " <<decaytime <<"[ns]" <<G4endl;}
436#endif
437  return  decaytime;       
438}
439
440////////////////////////////////////////////////////////////////////////////////
441//
442//
443// GetDecayTimeBin
444//
445//
446//
447G4int G4RadioactiveDecay::GetDecayTimeBin(const G4double aDecayTime)
448{
449  for (G4int i = 0; i < NDecayBin; i++) 
450    {
451      if ( aDecayTime > DBin[i]) return i+1;     
452    }
453  return  1;
454}
455////////////////////////////////////////////////////////////////////////////////
456//
457//
458// GetMeanLifeTime
459//
460// this is required by the base class
461//
462G4double G4RadioactiveDecay::GetMeanLifeTime(const G4Track& theTrack,
463                                             G4ForceCondition* )
464{
465  // For varience reduction implementation the time is set to 0 so as to
466  // force the particle to decay immediately.
467  // in analogueMC mode it return the particles meanlife.
468  //
469  G4double meanlife = 0.;
470  if (AnalogueMC) {
471    const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
472    G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
473    G4double theLife = theParticleDef->GetPDGLifeTime();
474     
475#ifdef G4VERBOSE
476    if (GetVerboseLevel()>2)
477      {
478        G4cout <<"G4RadioactiveDecay::GetMeanLifeTime() " <<G4endl;
479        G4cout <<"KineticEnergy:" <<theParticle->GetKineticEnergy()/GeV <<"[GeV]";
480        G4cout <<"Mass:" <<theParticle->GetMass()/GeV <<"[GeV]"; 
481        G4cout <<"Life time: " <<theLife/ns <<"[ns]" << G4endl;
482      }
483#endif
484    if (theParticleDef->GetPDGStable()) {meanlife = DBL_MAX;}
485    else if (theLife < 0.0) {meanlife = DBL_MAX;}
486    else {meanlife = theLife;}
487  // set the meanlife to zero for excited istopes which is not in the RDM database
488    if (((const G4Ions*)(theParticleDef))->GetExcitationEnergy() > 0. && meanlife == DBL_MAX) {meanlife = 0.;}
489  }
490#ifdef G4VERBOSE
491  if (GetVerboseLevel()>1)
492    {G4cout <<"mean life time: " <<meanlife <<"[ns]" <<G4endl;}
493#endif
494 
495  return  meanlife;
496}
497////////////////////////////////////////////////////////////////////////////////
498//
499//
500// GetMeanFreePath
501//
502// it is of similar functions to GetMeanFreeTime
503//
504G4double G4RadioactiveDecay::GetMeanFreePath (const G4Track& aTrack,
505                                              G4double, G4ForceCondition*)
506{
507  // constants
508  G4bool isOutRange ;
509 
510  // get particle
511  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
512 
513  // returns the mean free path in GEANT4 internal units
514  G4double pathlength;
515  G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
516  G4double aCtau = c_light * aParticleDef->GetPDGLifeTime();
517  G4double aMass = aParticle->GetMass();
518 
519#ifdef G4VERBOSE
520  if (GetVerboseLevel()>2) {
521    G4cout << "G4RadioactiveDecay::GetMeanFreePath() "<< G4endl;
522    G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV <<"[GeV]";
523    G4cout << "Mass:" << aMass/GeV <<"[GeV]";
524    G4cout << "c*Tau:" << aCtau/m <<"[m]" <<G4endl;
525  }
526#endif
527 
528  // check if the particle is stable?
529  if (aParticleDef->GetPDGStable()) {
530    pathlength = DBL_MAX;
531   
532  } else if (aCtau < 0.0) {
533    pathlength =  DBL_MAX;
534   
535    //check if the particle has very short life time ?
536  } else if (aCtau < DBL_MIN) {
537    pathlength =  DBL_MIN;
538   
539    //check if zero mass
540  } else if (aMass <  DBL_MIN)  {
541    pathlength =  DBL_MAX;
542#ifdef G4VERBOSE
543    if (GetVerboseLevel()>1) {
544      G4cerr << " Zero Mass particle " << G4endl;
545    }
546#endif
547   } else {
548     //calculate the mean free path
549     // by using normalized kinetic energy (= Ekin/mass)
550     G4double   rKineticEnergy = aParticle->GetKineticEnergy()/aMass;
551     if ( rKineticEnergy > HighestBinValue) {
552       // beta >> 1
553       pathlength = ( rKineticEnergy + 1.0)* aCtau;
554     } else if ( rKineticEnergy > LowestBinValue) {
555       // check if aPhysicsTable exists
556       if (aPhysicsTable == 0) BuildPhysicsTable(*aParticleDef);
557       // beta is in the range valid for PhysicsTable
558       pathlength = aCtau *
559         ((*aPhysicsTable)(0))-> GetValue(rKineticEnergy,isOutRange);
560     } else if ( rKineticEnergy < DBL_MIN ) {
561       // too slow particle
562#ifdef G4VERBOSE
563       if (GetVerboseLevel()>2) {
564         G4cout << "G4Decay::GetMeanFreePath()   !!particle stops!!";
565         G4cout << aParticleDef->GetParticleName() << G4endl;
566         G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV <<"[GeV]";
567       }
568#endif
569       pathlength = DBL_MIN;
570     } else {
571       // beta << 1
572       pathlength = (aParticle->GetTotalMomentum())/aMass*aCtau ;
573     }
574   }
575#ifdef G4VERBOSE
576   if (GetVerboseLevel()>1) {
577     G4cout << "mean free path: "<< pathlength/m << "[m]" << G4endl;
578   }
579#endif
580   return  pathlength;
581}
582////////////////////////////////////////////////////////////////////////////////
583//
584//
585//
586//
587void G4RadioactiveDecay::BuildPhysicsTable(const G4ParticleDefinition&)
588{
589  // if aPhysicsTableis has already been created, do nothing
590  if (aPhysicsTable != 0) return;
591
592  // create  aPhysicsTable
593  if (GetVerboseLevel()>1) G4cerr <<" G4Decay::BuildPhysicsTable() "<< G4endl;
594  aPhysicsTable = new G4PhysicsTable(1);
595
596  //create physics vector
597  G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
598                                                       LowestBinValue,
599                                                       HighestBinValue,
600                                                       TotBin);
601
602  G4double beta, gammainv;
603  // fill physics Vector
604  G4int i;
605  for ( i = 0 ; i < TotBin ; i++ ) {
606      gammainv = 1.0/(aVector->GetLowEdgeEnergy(i) + 1.0);
607      beta  = std::sqrt((1.0 - gammainv)*(1.0 +gammainv));
608      aVector->PutValue(i, beta/gammainv);
609  }
610  aPhysicsTable->insert(aVector);
611}
612///////////////////////////////////////////////////////////////////////////////
613//
614// LoadDecayTable
615//
616//     To load the decay scheme from the Radioactivity database for
617//     theParentNucleus.
618//
619G4DecayTable *G4RadioactiveDecay::LoadDecayTable (G4ParticleDefinition &theParentNucleus)
620{
621  //
622  //
623  // Create and initialise variables used in the method.
624  //
625  G4DecayTable *theDecayTable = new G4DecayTable();
626  //
627  //
628  // Determine the filename of the file containing radioactive decay data.  Open
629  // it.
630  //
631  G4int A    = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
632  G4int Z    = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
633  G4double E = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
634
635  if ( !getenv("G4RADIOACTIVEDATA") ) {
636    G4cout << "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files." << G4endl;
637    throw G4HadronicException(__FILE__, __LINE__, 
638              "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files.");
639  }
640  G4String dirName = getenv("G4RADIOACTIVEDATA");
641  LoadedNuclei.push_back(theParentNucleus.GetParticleName());
642  std::sort( LoadedNuclei.begin(), LoadedNuclei.end() );
643  // sort needed to allow binary_search
644
645  std::ostringstream os;
646  os <<dirName <<"/z" <<Z <<".a" <<A <<'\0';
647  G4String file = os.str();
648
649 
650  std::ifstream DecaySchemeFile(file);
651  G4bool found(false);
652  if (DecaySchemeFile) { 
653    //
654    //
655    // Initialise variables used for reading in radioactive decay data.
656    //
657    G4int    nMode = 7;
658    G4bool   modeFirstRecord[7];
659    G4double modeTotalBR[7];
660    G4double modeSumBR[7];
661    G4int i;
662    for (i=0; i<nMode; i++) {
663      modeFirstRecord[i] = true;
664      modeSumBR[i]       = 0.0;
665    }
666   
667    G4bool complete(false);
668    char inputChars[80]={' '};
669    G4String inputLine;
670    G4String recordType("");
671    G4RadioactiveDecayMode theDecayMode;
672    G4double a(0.0);
673    G4double b(0.0);
674    G4double c(0.0);
675    G4double n(1.0);
676    G4double e0;
677    //
678    //
679    // Go through each record in the data file until you identify the decay
680    // data relating to the nuclide of concern.
681    //
682    //  while (!complete && -DecaySchemeFile.getline(inputChars, 80).eof() != EOF)
683    while (!complete && !DecaySchemeFile.getline(inputChars, 80).eof()) {
684      inputLine = inputChars;
685      //    G4String::stripType stripend(1);
686      //    inputLine = inputLine.strip(stripend);
687      inputLine = inputLine.strip(1);
688      if (inputChars[0] != '#' && inputLine.length() != 0) {
689        std::istringstream tmpStream(inputLine);
690        if (inputChars[0] == 'P') {
691          //
692          //
693          // Nucleus is a parent type.  Check the excitation level to see if it matches
694          // that of theParentNucleus
695          //
696          tmpStream >>recordType >>a >>b;
697          if (found) {complete = true;}
698          else {found = (std::abs(a*keV - E)<levelTolerance);}
699        } else if (found) {
700          //
701          //
702          // The right part of the radioactive decay data file has been found.  Search
703          // through it to determine the mode of decay of the subsequent records.
704          //
705          if (inputChars[0] == 'W') {
706#ifdef G4VERBOSE
707            if (GetVerboseLevel()>0) {
708              // a comment line identified and print out the message
709              //
710              G4cout << " Warning in G4RadioactiveDecay::LoadDecayTable " << G4endl;
711              G4cout << "   In data file " << file << G4endl;
712              G4cout << "   " << inputLine << G4endl;
713            }
714#endif
715          }     
716          else {
717            tmpStream >>theDecayMode >>a >>b >>c;
718            a/=1000.;
719            c/=1000.;
720                 
721            //  cout<< "The decay mode is [LoadTable] "<<theDecayMode<<G4endl;
722                 
723            switch (theDecayMode) {
724            case IT:
725              {
726                //
727                //
728                // Decay mode is isomeric transition.
729                //
730                G4ITDecayChannel *anITChannel = new G4ITDecayChannel
731                  (GetVerboseLevel(), (const G4Ions*) &theParentNucleus, b);
732                anITChannel->SetICM(applyICM);
733                anITChannel->SetARM(applyARM);
734                anITChannel->SetHLThreshold(halflifethreshold);
735                theDecayTable->Insert(anITChannel);
736                break;
737              }
738            case BetaMinus:
739              {
740                //
741                //
742                // Decay mode is beta-.
743                //
744                if (modeFirstRecord[1])
745                  {modeFirstRecord[1] = false; modeTotalBR[1] = b;}
746                else {
747                  if (c > 0.) {
748                    // to work out the Fermi function normalization factor first
749                    G4BetaFermiFunction* aBetaFermiFunction = new G4BetaFermiFunction (A, (Z+1));
750                    e0 = c*MeV/0.511;
751                    n = aBetaFermiFunction->GetFFN(e0);
752                               
753                    // now to work out the histogram and initialise the random generator
754                    G4int npti = 100;                           
755                    G4double* pdf = new G4double[npti];
756                    G4int ptn;
757                    G4double g,e,ee,f;
758                    ee = e0+1.;
759                    for (ptn=0; ptn<npti; ptn++) {
760                      // e =e0*(ptn+1.)/102.;
761                      // bug fix (#662) (flei, 22/09/2004)
762                      e =e0*(ptn+0.5)/100.;
763                      g = e+1.;
764                      f = std::sqrt(g*g-1)*(ee-g)*(ee-g)*g;
765                      // Special treatment for K-40 (problem #1068) (flei,06/05/2010)
766                      if (theParentNucleus.GetParticleName() == "K40[0.0]") f *=
767                        (std::pow((g*g-1),3)+std::pow((ee-g),6)+7*(g*g-1)*std::pow((ee-g),2)*(g*g-1+std::pow((ee-g),2)));
768                      pdf[ptn] = f*aBetaFermiFunction->GetFF(e);
769                    }             
770                    RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti); 
771                    G4BetaMinusDecayChannel *aBetaMinusChannel = new
772                      G4BetaMinusDecayChannel (GetVerboseLevel(), &theParentNucleus,
773                                               b, c*MeV, a*MeV, n, FBeta, aRandomEnergy);
774                    aBetaMinusChannel->SetICM(applyICM);
775                    aBetaMinusChannel->SetARM(applyARM);
776                    aBetaMinusChannel->SetHLThreshold(halflifethreshold);
777                    theDecayTable->Insert(aBetaMinusChannel);
778                    modeSumBR[1] += b;
779                    delete[] pdf;
780                    delete aBetaFermiFunction;
781                  }
782                }
783                break;
784              case BetaPlus:
785                //
786                //
787                // Decay mode is beta+.
788                //
789                if (modeFirstRecord[2])
790                  {modeFirstRecord[2] = false; modeTotalBR[2] = b;}
791                else {
792                  //  e0 = c*MeV/0.511;
793                  // bug fix (#662) (flei, 22/09/2004)
794                  // need to test e0 as there are some data files (e.g. z67.a162) which have entries for beta+
795                  // with Q < 2Me
796                  //
797                  e0 = c*MeV/0.511 -2.;
798                  if (e0 > 0.) {
799                    G4BetaFermiFunction* aBetaFermiFunction = new G4BetaFermiFunction (A, -(Z-1));
800                               
801                    n = aBetaFermiFunction->GetFFN(e0);
802                               
803                    // now to work out the histogram and initialise the random generator
804                    G4int npti = 100;                           
805                    G4double* pdf = new G4double[npti];
806                    G4int ptn;
807                    G4double g,e,ee,f;
808                    ee = e0+1.;
809                    for (ptn=0; ptn<npti; ptn++) {
810                      // e =e0*(ptn+1.)/102.;
811                      // bug fix (#662) (flei, 22/09/2004)
812                      e =e0*(ptn+0.5)/100.;
813                      g = e+1.;
814                      f = std::sqrt(g*g-1)*(ee-g)*(ee-g)*g;
815                      pdf[ptn] = f*aBetaFermiFunction->GetFF(e);
816                    }
817                    RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti); 
818                    G4BetaPlusDecayChannel *aBetaPlusChannel = new 
819                      G4BetaPlusDecayChannel (GetVerboseLevel(), &theParentNucleus,
820                                              b, c*MeV, a*MeV, n, FBeta, aRandomEnergy);
821                    aBetaPlusChannel->SetICM(applyICM);
822                    aBetaPlusChannel->SetARM(applyARM);
823                    aBetaPlusChannel->SetHLThreshold(halflifethreshold);
824                    theDecayTable->Insert(aBetaPlusChannel);
825                    modeSumBR[2] += b;
826                               
827                    delete[] pdf;
828                    delete aBetaFermiFunction;       
829                  }
830                }
831                break;
832              case KshellEC:
833                //
834                //
835                // Decay mode is K-electron capture.
836                //
837                if (modeFirstRecord[3])
838                  {modeFirstRecord[3] = false; modeTotalBR[3] = b;}
839                else {
840                  G4KshellECDecayChannel *aKECChannel = new
841                    G4KshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
842                                            b, c*MeV, a*MeV);
843                  aKECChannel->SetICM(applyICM);
844                  aKECChannel->SetARM(applyARM);
845                  aKECChannel->SetHLThreshold(halflifethreshold);
846                  theDecayTable->Insert(aKECChannel);
847                  //delete aKECChannel;
848                  modeSumBR[3] += b;
849                }
850                break;
851              case LshellEC:
852                //
853                //
854                // Decay mode is L-electron capture.
855                //
856                if (modeFirstRecord[4])
857                  {modeFirstRecord[4] = false; modeTotalBR[4] = b;}
858                else {
859                  G4LshellECDecayChannel *aLECChannel = new
860                    G4LshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
861                                            b, c*MeV, a*MeV);
862                  aLECChannel->SetICM(applyICM);
863                  aLECChannel->SetARM(applyARM);
864                  aLECChannel->SetHLThreshold(halflifethreshold);
865                  theDecayTable->Insert(aLECChannel);
866                  //delete aLECChannel;
867                  modeSumBR[4] += b;
868                }
869                break;
870              case MshellEC:
871                //
872                //
873                // Decay mode is M-electron capture. In this implementation it is added to L-shell case
874                //
875                if (modeFirstRecord[5])
876                  {modeFirstRecord[5] = false; modeTotalBR[5] = b;}
877                else {
878                  G4MshellECDecayChannel *aMECChannel = new
879                    G4MshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
880                                            b, c*MeV, a*MeV);
881                  aMECChannel->SetICM(applyICM);
882                  aMECChannel->SetARM(applyARM);
883                  aMECChannel->SetHLThreshold(halflifethreshold);
884                  theDecayTable->Insert(aMECChannel);
885                  modeSumBR[5] += b;
886                }
887                break;
888              case Alpha:
889                //
890                //
891                // Decay mode is alpha.
892                //
893                if (modeFirstRecord[6])
894                  {modeFirstRecord[6] = false; modeTotalBR[6] = b;}
895                else {
896                  G4AlphaDecayChannel *anAlphaChannel = new
897                    G4AlphaDecayChannel (GetVerboseLevel(), &theParentNucleus,
898                                         b, c*MeV, a*MeV);
899                  anAlphaChannel->SetICM(applyICM);
900                  anAlphaChannel->SetARM(applyARM);
901                  anAlphaChannel->SetHLThreshold(halflifethreshold);
902                  theDecayTable->Insert(anAlphaChannel);
903                  //          delete anAlphaChannel;
904                  modeSumBR[6] += b;
905                }
906                break;
907              case ERROR:
908              default:
909                G4Exception("G4RadioactiveDecay::LoadDecayTable()", "601",
910                            FatalException, "Error in decay mode selection");
911               
912              }
913            }
914          }
915        }
916      }
917    }   
918    //
919    //
920    // Go through the decay table and make sure that the branching ratios are
921    // correctly normalised.
922    //
923    G4VDecayChannel       *theChannel             = 0;
924    G4NuclearDecayChannel *theNuclearDecayChannel = 0;
925    G4String mode                     = "";
926    G4int j                           = 0;
927    G4double theBR                    = 0.0;
928    for (i=0; i<theDecayTable->entries(); i++) {
929      theChannel             = theDecayTable->GetDecayChannel(i);
930      theNuclearDecayChannel = static_cast<G4NuclearDecayChannel *>(theChannel);
931      theDecayMode           = theNuclearDecayChannel->GetDecayMode();
932      j          = 0;
933      if (theDecayMode != IT) {
934        theBR = theChannel->GetBR();
935        theChannel->SetBR(theBR*modeTotalBR[theDecayMode]/modeSumBR[theDecayMode]);
936      }
937    } 
938  }             
939  DecaySchemeFile.close();             
940  if (!found && E > 0.) {
941    // cases where IT cascade for exited isotopes without entry in RDM database
942    // Decay mode is isomeric transition.
943    //
944    G4ITDecayChannel *anITChannel = new G4ITDecayChannel
945      (GetVerboseLevel(), (const G4Ions*) &theParentNucleus, 1);
946    anITChannel->SetICM(applyICM);
947    anITChannel->SetARM(applyARM);
948    anITChannel->SetHLThreshold(halflifethreshold);
949    theDecayTable->Insert(anITChannel);
950  } 
951  if (!theDecayTable) {
952    //
953    // There is no radioactive decay data for this nucleus.  Return a null
954    // decay table.
955    //
956    G4cerr <<"G4RadoactiveDecay::LoadDecayTable() : cannot find ion radioactive decay file " <<G4endl;
957    theDecayTable = 0;
958    return theDecayTable;
959  }     
960  if (theDecayTable && GetVerboseLevel()>1)
961    {
962      G4cout <<"G4RadioactiveDecay::LoadDecayTable()" << G4endl;
963      G4cout << "  No. of  entries: "<< theDecayTable->entries() <<G4endl;
964      theDecayTable ->DumpInfo();
965    }
966
967  return theDecayTable;
968}
969
970////////////////////////////////////////////////////////////////////////
971//
972//
973void G4RadioactiveDecay::SetDecayRate(G4int theZ, G4int theA, G4double theE, 
974                                       G4int theG, std::vector<G4double> theRates, 
975                                       std::vector<G4double> theTaos)
976{ 
977  //fill the decay rate vector
978  theDecayRate.SetZ(theZ);
979  theDecayRate.SetA(theA);
980  theDecayRate.SetE(theE);
981  theDecayRate.SetGeneration(theG);
982  theDecayRate.SetDecayRateC(theRates);
983  theDecayRate.SetTaos(theTaos);
984}
985//////////////////////////////////////////////////////////////////////////
986//
987//
988void G4RadioactiveDecay::AddDecayRateTable(const G4ParticleDefinition &theParentNucleus)
989{
990  // 1) To calculate all the coefficiecies required to derive the radioactivities for all
991  // progeny of theParentNucleus
992  //
993  // 2) Add the coefficiencies to the decay rate table vector
994  //
995 
996  //
997  // Create and initialise variables used in the method.
998  //
999
1000  theDecayRateVector.clear();
1001 
1002  G4int nGeneration = 0;
1003  std::vector<G4double> rates;
1004  std::vector<G4double> taos;
1005 
1006  // start rate is -1.
1007  rates.push_back(-1.);
1008  //
1009  //
1010  G4int A = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
1011  G4int Z = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
1012  G4double E = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
1013  G4double tao = theParentNucleus.GetPDGLifeTime();
1014  taos.push_back(tao);
1015  G4int nEntry = 0;
1016 
1017  //fill the decay rate with the intial isotope data
1018  SetDecayRate(Z,A,E,nGeneration,rates,taos);
1019
1020  // store the decay rate in decay rate vector
1021  theDecayRateVector.push_back(theDecayRate);
1022  nEntry++;
1023
1024  // now start treating the sencondary generations..
1025
1026  G4bool stable = false;
1027  G4int i;
1028  G4int j;
1029  G4VDecayChannel       *theChannel             = 0;
1030  G4NuclearDecayChannel *theNuclearDecayChannel = 0;
1031  G4ITDecayChannel *theITChannel = 0;
1032  G4BetaMinusDecayChannel *theBetaMinusChannel = 0;
1033  G4BetaPlusDecayChannel *theBetaPlusChannel = 0;
1034  G4AlphaDecayChannel *theAlphaChannel = 0;
1035  G4RadioactiveDecayMode theDecayMode;
1036  //  G4NuclearLevelManager levelManager;
1037  //const G4NuclearLevel* level;
1038  G4double theBR = 0.0;
1039  G4int AP = 0;
1040  G4int ZP = 0;
1041  G4int AD = 0;
1042  G4int ZD = 0;
1043  G4double EP = 0.;
1044  std::vector<G4double> TP;
1045  std::vector<G4double> RP;
1046  G4ParticleDefinition *theDaughterNucleus;
1047  G4double daughterExcitation;
1048  G4ParticleDefinition *aParentNucleus;
1049  G4IonTable* theIonTable;
1050  G4DecayTable *aTempDecayTable;
1051  G4double theRate;
1052  G4double TaoPlus;
1053  G4int nS = 0;
1054  G4int nT = nEntry;
1055  G4double brs[7];
1056  //
1057  theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
1058 
1059  while (!stable) {
1060    nGeneration++;
1061    for (j = nS; j< nT; j++) {
1062      ZP = theDecayRateVector[j].GetZ();
1063      AP = theDecayRateVector[j].GetA();
1064      EP = theDecayRateVector[j].GetE();
1065      RP = theDecayRateVector[j].GetDecayRateC();
1066      TP = theDecayRateVector[j].GetTaos();     
1067      if (GetVerboseLevel()>0){
1068        G4cout <<"G4RadioactiveDecay::AddDecayRateTable : "
1069               << " daughters of ("<< ZP <<", "<<AP<<", "
1070               << EP <<") "
1071               << " are being calculated. "       
1072               <<" generation = "
1073               << nGeneration << G4endl;
1074      }
1075      aParentNucleus = theIonTable->GetIon(ZP,AP,EP);
1076      if (!IsLoaded(*aParentNucleus)){
1077        aParentNucleus->SetDecayTable(LoadDecayTable(*aParentNucleus));
1078      }
1079      aTempDecayTable = aParentNucleus->GetDecayTable();
1080      //
1081      //
1082      // Go through the decay table and to combine the same decay channels
1083      //
1084      for (i=0; i< 7; i++) brs[i] = 0.0;
1085     
1086      G4DecayTable *theDecayTable = new G4DecayTable();
1087     
1088      for (i=0; i<aTempDecayTable->entries(); i++) {
1089        theChannel             = aTempDecayTable->GetDecayChannel(i);
1090        theNuclearDecayChannel = static_cast<G4NuclearDecayChannel *>(theChannel);
1091        theDecayMode           = theNuclearDecayChannel->GetDecayMode();
1092        daughterExcitation = theNuclearDecayChannel->GetDaughterExcitation ();
1093        theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus () ;
1094        AD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
1095        ZD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber(); 
1096        G4NuclearLevelManager * levelManager = G4NuclearLevelStore::GetInstance()->GetManager (ZD, AD);
1097        if ( levelManager->NumberOfLevels() ) {
1098          const G4NuclearLevel* level = levelManager->NearestLevel (daughterExcitation);
1099
1100          if (std::abs(daughterExcitation - level->Energy()) < levelTolerance) {
1101           
1102            // Level hafe life is in ns and the gate as 1 micros by default
1103            if ( theDecayMode == 0 && level->HalfLife()*ns >= halflifethreshold ){   
1104              // some further though may needed here
1105              theDecayTable->Insert(theChannel);
1106            } 
1107            else{
1108              brs[theDecayMode] += theChannel->GetBR();
1109            }
1110          }
1111          else {
1112            brs[theDecayMode] += theChannel->GetBR();
1113          }
1114        }
1115        else{
1116          brs[theDecayMode] += theChannel->GetBR();
1117        }
1118      }     
1119      brs[2] = brs[2]+brs[3]+brs[4]+brs[5];
1120      brs[3] = brs[4] =brs[5] =  0.0;
1121      for (i= 0; i<7; i++){
1122        if (brs[i] > 0) {
1123          switch ( i ) {
1124          case 0:
1125            //
1126            //
1127            // Decay mode is isomeric transition.
1128            //
1129           
1130            theITChannel =  new G4ITDecayChannel
1131              (0, (const G4Ions*) aParentNucleus, brs[0]);
1132           
1133            theDecayTable->Insert(theITChannel);
1134            break;
1135           
1136          case 1:
1137            //
1138            //
1139            // Decay mode is beta-.
1140            //
1141            theBetaMinusChannel = new G4BetaMinusDecayChannel (0, aParentNucleus,
1142                                                               brs[1], 0.*MeV, 0.*MeV, 1, false, 0);
1143            theDecayTable->Insert(theBetaMinusChannel);
1144           
1145            break;
1146           
1147          case 2:
1148            //
1149            //
1150            // Decay mode is beta+ + EC.
1151            //
1152            theBetaPlusChannel = new G4BetaPlusDecayChannel (GetVerboseLevel(), aParentNucleus,
1153                                                             brs[2], 0.*MeV, 0.*MeV, 1, false, 0);
1154            theDecayTable->Insert(theBetaPlusChannel);
1155            break;                   
1156           
1157          case 6:
1158            //
1159            //
1160            // Decay mode is alpha.
1161            //
1162            theAlphaChannel = new G4AlphaDecayChannel (GetVerboseLevel(), aParentNucleus,
1163                                                       brs[6], 0.*MeV, 0.*MeV);
1164            theDecayTable->Insert(theAlphaChannel);
1165            break;
1166           
1167          default:
1168            break;
1169          }
1170        }
1171      }
1172       
1173      //
1174      // loop over all braches in theDecayTable
1175      //
1176      for ( i=0; i<theDecayTable->entries(); i++){
1177        theChannel             = theDecayTable->GetDecayChannel(i);
1178        theNuclearDecayChannel = static_cast<G4NuclearDecayChannel *>(theChannel);
1179        theBR = theChannel->GetBR();
1180        theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus();
1181       
1182        //
1183        // now test if the daughterNucleus is a valid one
1184        //
1185        if (IsApplicable(*theDaughterNucleus) && theBR
1186            && aParentNucleus != theDaughterNucleus ) {
1187          // need to make sure daugher has decaytable
1188          if (!IsLoaded(*theDaughterNucleus)){
1189            theDaughterNucleus->SetDecayTable(LoadDecayTable(*theDaughterNucleus));
1190          }
1191          if (theDaughterNucleus->GetDecayTable()->entries() ) {
1192            //
1193            A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
1194            Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
1195            E = ((const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy();
1196         
1197            TaoPlus = theDaughterNucleus->GetPDGLifeTime();
1198            //          cout << TaoPlus <<G4endl;
1199            if (TaoPlus > 0.) {
1200              // first set the taos, one simply need to add to the parent ones
1201              taos.clear();
1202              taos = TP;
1203              taos.push_back(TaoPlus);
1204              // now calculate the coefficiencies
1205              //
1206              // they are in two parts, first the les than n ones
1207              rates.clear();
1208              size_t k;
1209              for (k = 0; k < RP.size(); k++){
1210                theRate = TP[k]/(TP[k]-TaoPlus) * theBR * RP[k];
1211                rates.push_back(theRate);
1212              }
1213              //
1214              // the sencond part: the n:n coefficiency
1215              theRate = 0.;
1216              for (k = 0; k < RP.size(); k++){
1217                theRate -=TaoPlus/(TP[k]-TaoPlus) * theBR * RP[k];
1218              }
1219              rates.push_back(theRate);               
1220              SetDecayRate (Z,A,E,nGeneration,rates,taos);
1221              theDecayRateVector.push_back(theDecayRate);
1222              nEntry++;
1223            } 
1224          }
1225        } 
1226        // end of testing daughter nucleus
1227      }
1228      // end of i loop( the branches)
1229    }
1230    //end of for j loop
1231    nS = nT;
1232    nT = nEntry;
1233    if (nS == nT) stable = true;
1234  }
1235 
1236  //end of while loop
1237  // the calculation completed here
1238 
1239 
1240  // fill the first part of the decay rate table
1241  // which is the name of the original particle (isotope)
1242  //
1243  theDecayRateTable.SetIonName(theParentNucleus.GetParticleName()); 
1244  //
1245  //
1246  // now fill the decay table with the newly completed decay rate vector
1247  //
1248 
1249  theDecayRateTable.SetItsRates(theDecayRateVector);
1250 
1251  //
1252  // finally add the decayratetable to the tablevector
1253  //
1254  theDecayRateTableVector.push_back(theDecayRateTable);
1255}
1256 
1257////////////////////////////////////////////////////////////////////////////////
1258//
1259//
1260// SetSourceTimeProfile
1261//
1262//  read in the source time profile function (histogram)
1263//
1264
1265void G4RadioactiveDecay::SetSourceTimeProfile(G4String filename)
1266{
1267  std::ifstream infile ( filename, std::ios::in );
1268  if ( !infile ) G4Exception(__FILE__, G4inttostring(__LINE__), FatalException,  "Unable to open source data file" );
1269 
1270  float bin, flux;
1271  NSourceBin = -1;
1272  while (infile >> bin >> flux ) {
1273    NSourceBin++;
1274    if (NSourceBin > 99)  G4Exception(__FILE__, G4inttostring(__LINE__), FatalException,  "input source data file too big (>100 rows)" );
1275    SBin[NSourceBin] = bin * s;
1276    SProfile[NSourceBin] = flux;
1277  }
1278  SetAnalogueMonteCarlo(0);
1279#ifdef G4VERBOSE
1280  if (GetVerboseLevel()>1)
1281    {G4cout <<" Source Timeprofile Nbin = " << NSourceBin <<G4endl;}
1282#endif
1283}
1284
1285////////////////////////////////////////////////////////////////////////////////
1286//
1287//
1288// SetDecayBiasProfile
1289//
1290// read in the decay bias scheme function (histogram)
1291//
1292void G4RadioactiveDecay::SetDecayBias(G4String filename)
1293{
1294  std::ifstream infile ( filename, std::ios::in);
1295  if ( !infile ) G4Exception(__FILE__, G4inttostring(__LINE__), FatalException,  "Unable to open bias data file" );
1296 
1297  float bin, flux;
1298  NDecayBin = -1;
1299  while (infile >> bin >> flux ) {
1300    NDecayBin++;
1301    if (NDecayBin > 99)  G4Exception(__FILE__, G4inttostring(__LINE__), FatalException,  "input bias data file too big (>100 rows)" );
1302    DBin[NDecayBin] = bin * s;
1303    DProfile[NDecayBin] = flux;
1304  }
1305  G4int i ;
1306  for ( i = 1; i<= NDecayBin; i++) DProfile[i] += DProfile[i-1];
1307  for ( i = 0; i<= NDecayBin; i++) DProfile[i] /= DProfile[NDecayBin];
1308  SetAnalogueMonteCarlo(0);
1309#ifdef G4VERBOSE
1310  if (GetVerboseLevel()>1)
1311    {G4cout <<" Decay Bias Profile  Nbin = " << NDecayBin <<G4endl;}
1312#endif
1313}
1314
1315////////////////////////////////////////////////////////////////////////////////
1316//
1317//
1318// DecayIt
1319//
1320G4VParticleChange* G4RadioactiveDecay::DecayIt(const G4Track& theTrack, const G4Step& )
1321{
1322  //
1323  // Initialize the G4ParticleChange object. Get the particle details and the
1324  // decay table.
1325  //
1326  fParticleChangeForRadDecay.Initialize(theTrack);
1327  const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
1328  G4ParticleDefinition *theParticleDef = theParticle->GetDefinition();
1329
1330  // First check whether RDM applies to the current logical volume
1331  //
1332  if(!std::binary_search(ValidVolumes.begin(),
1333                    ValidVolumes.end(), 
1334                    theTrack.GetVolume()->GetLogicalVolume()->GetName()))
1335    {
1336#ifdef G4VERBOSE
1337      if (GetVerboseLevel()>0)
1338        {
1339          G4cout <<"G4RadioactiveDecay::DecayIt : "
1340                 << theTrack.GetVolume()->GetLogicalVolume()->GetName()
1341                 << " is not selected for the RDM"<< G4endl;
1342          G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
1343          G4cout << " The Valid volumes are " << G4endl;
1344          for (size_t i = 0; i< ValidVolumes.size(); i++)
1345            G4cout << ValidVolumes[i] << G4endl;
1346        }
1347#endif
1348      fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1349      //
1350      //
1351      // Kill the parent particle.
1352      //
1353      fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
1354      fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1355      ClearNumberOfInteractionLengthLeft();
1356      return &fParticleChangeForRadDecay;
1357    }
1358   
1359  // now check is the particle is valid for RDM
1360  //
1361  if (!(IsApplicable(*theParticleDef)))
1362    { 
1363      //
1364      // The particle is not a Ion or outside the nucleuslimits for decay
1365      //
1366#ifdef G4VERBOSE
1367      if (GetVerboseLevel()>0)
1368        {
1369          G4cerr <<"G4RadioactiveDecay::DecayIt : "
1370                 <<theParticleDef->GetParticleName() 
1371                 << " is not a valid nucleus for the RDM"<< G4endl;
1372        }
1373#endif
1374      fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1375      //
1376      //
1377      // Kill the parent particle.
1378      //
1379      fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
1380      fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1381      ClearNumberOfInteractionLengthLeft();
1382      return &fParticleChangeForRadDecay;
1383    }
1384 
1385  if (!IsLoaded(*theParticleDef))
1386    {
1387      theParticleDef->SetDecayTable(LoadDecayTable(*theParticleDef));
1388    }
1389  G4DecayTable *theDecayTable = theParticleDef->GetDecayTable();
1390 
1391  if  (theDecayTable == 0 || theDecayTable->entries() == 0 )
1392    {
1393      //
1394      //
1395      // There are no data in the decay table.  Set the particle change parameters
1396      // to indicate this.
1397      //
1398#ifdef G4VERBOSE
1399      if (GetVerboseLevel()>0)
1400        {
1401          G4cerr <<"G4RadioactiveDecay::DecayIt : decay table not defined  for ";
1402          G4cerr <<theParticleDef->GetParticleName() <<G4endl;
1403        }
1404#endif
1405      fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1406      //
1407      //
1408      // Kill the parent particle.
1409      //
1410      fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
1411      fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1412      ClearNumberOfInteractionLengthLeft();
1413      return &fParticleChangeForRadDecay;
1414    }
1415  else 
1416    { 
1417      //
1418      // now try to  decay it
1419      //
1420      G4double energyDeposit = 0.0;
1421      G4double finalGlobalTime = theTrack.GetGlobalTime();
1422      G4int index;
1423      G4ThreeVector currentPosition;
1424      currentPosition = theTrack.GetPosition();
1425     
1426      // check whether use Analogue or VR implementation
1427      //
1428      if (AnalogueMC){
1429        //
1430        // Aanalogue MC
1431#ifdef G4VERBOSE
1432        if (GetVerboseLevel()>0)
1433          {
1434            G4cout <<"DecayIt:  Analogue MC version " << G4endl;
1435          }
1436#endif
1437        //
1438        G4DecayProducts *products = DoDecay(*theParticleDef);
1439        //
1440        // check if the product is the same as input and kill the track if necessary to prevent infinite loop
1441        // (11/05/10, F.Lei)
1442        //
1443        if ( products->entries() == 1) {
1444          fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1445          fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
1446          fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1447          ClearNumberOfInteractionLengthLeft();
1448          return &fParticleChangeForRadDecay;
1449        }       
1450        //
1451        // Get parent particle information and boost the decay products to the
1452        // laboratory frame based on this information.
1453        //
1454        G4double ParentEnergy = theParticle->GetTotalEnergy();
1455        G4ThreeVector ParentDirection(theParticle->GetMomentumDirection());
1456       
1457        if (theTrack.GetTrackStatus() == fStopButAlive )
1458          {
1459            //
1460            //
1461            // The particle is decayed at rest.
1462            //
1463            // since the time is still for rest particle in G4 we need to add the additional
1464            // time lapsed between the particle come to rest and the actual decay. This time
1465            // is simply sampled with the mean-life of the particle.
1466            // but we need to protect the case PDGTime < 0. (F.Lei 11/05/10)
1467            G4double temptime =  -std::log( G4UniformRand()) * theParticleDef->GetPDGLifeTime();
1468            if (temptime <0.) temptime =0.; 
1469            finalGlobalTime += temptime ;
1470            energyDeposit += theParticle->GetKineticEnergy();
1471          }
1472        else
1473          {
1474            //
1475            //
1476            // The particle is decayed in flight (PostStep case).
1477            //
1478            products->Boost( ParentEnergy, ParentDirection);
1479          }
1480        //
1481        //
1482        // Add products in theParticleChangeForRadDecay.
1483        //
1484        G4int numberOfSecondaries = products->entries();
1485        fParticleChangeForRadDecay.SetNumberOfSecondaries(numberOfSecondaries);
1486#ifdef G4VERBOSE
1487        if (GetVerboseLevel()>1) {
1488          G4cout <<"G4RadioactiveDecay::DecayIt : Decay vertex :";
1489          G4cout <<" Time: " <<finalGlobalTime/ns <<"[ns]";
1490          G4cout <<" X:" <<(theTrack.GetPosition()).x() /cm <<"[cm]";
1491          G4cout <<" Y:" <<(theTrack.GetPosition()).y() /cm <<"[cm]";
1492          G4cout <<" Z:" <<(theTrack.GetPosition()).z() /cm <<"[cm]";
1493          G4cout <<G4endl;
1494          G4cout <<"G4Decay::DecayIt  : decay products in Lab. Frame" <<G4endl;
1495          products->DumpInfo();
1496        }
1497#endif
1498        for (index=0; index < numberOfSecondaries; index++) 
1499          {
1500            G4Track* secondary = new G4Track
1501              (products->PopProducts(), finalGlobalTime, currentPosition);
1502            secondary->SetGoodForTrackingFlag();
1503                        secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
1504            fParticleChangeForRadDecay.AddSecondary(secondary);
1505          }
1506        delete products;
1507        //
1508        // end of analogue MC algarithm
1509        //
1510      }
1511      else {
1512        //
1513        // Varaice Reduction Method
1514        //
1515#ifdef G4VERBOSE
1516        if (GetVerboseLevel()>0)
1517          {
1518            G4cout <<"DecayIt:  Variance Reduction version " << G4endl;
1519          }
1520#endif
1521        if (!IsRateTableReady(*theParticleDef)) {
1522          // if the decayrates are not ready, calculate them and
1523          // add to the rate table vector
1524          AddDecayRateTable(*theParticleDef);
1525        }
1526        //retrieve the rates
1527        GetDecayRateTable(*theParticleDef);
1528        //
1529        // declare some of the variables required in the implementation
1530        //
1531        G4ParticleDefinition* parentNucleus;
1532        G4IonTable *theIonTable;
1533        G4int PZ;
1534        G4int PA;
1535        G4double PE;
1536        std::vector<G4double> PT;
1537        std::vector<G4double> PR;
1538        G4double taotime;
1539        G4double decayRate;
1540       
1541        size_t i;
1542        size_t j;
1543        G4int numberOfSecondaries;
1544        G4int totalNumberOfSecondaries = 0;
1545        G4double currentTime = 0.;
1546        G4int ndecaych;
1547        G4DynamicParticle* asecondaryparticle;
1548        //      G4DecayProducts* products = 0;
1549        std::vector<G4DynamicParticle*> secondaryparticles;
1550        std::vector<G4double> pw;
1551        std::vector<G4double> ptime;
1552        pw.clear();
1553        ptime.clear();
1554        //now apply the nucleus splitting
1555        //
1556        //
1557        for (G4int n = 0; n < NSplit; n++)
1558          {
1559            /*
1560            //
1561            // Get the decay time following the decay probability function
1562            // suppllied by user
1563            //
1564            G4double theDecayTime = GetDecayTime();
1565           
1566            G4int nbin = GetDecayTimeBin(theDecayTime);
1567           
1568            // claculate the first part of the weight function
1569           
1570            G4double weight1 =1./DProfile[nbin-1]
1571              *(DBin[nbin]-DBin[nbin-1])
1572              /NSplit;
1573            if (nbin > 1) {
1574               weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
1575                 *(DBin[nbin]-DBin[nbin-1])
1576                 /NSplit;}
1577            // it should be calculated in seconds
1578            weight1 /= s ;
1579            */
1580            //
1581            // loop over all the possible secondaries of the nucleus
1582            // the first one is itself.
1583            //
1584            for ( i = 0; i<theDecayRateVector.size(); i++){
1585              PZ = theDecayRateVector[i].GetZ();
1586              PA = theDecayRateVector[i].GetA();
1587              PE = theDecayRateVector[i].GetE();
1588              PT = theDecayRateVector[i].GetTaos();
1589              PR = theDecayRateVector[i].GetDecayRateC();
1590
1591              //
1592              // Get the decay time following the decay probability function
1593              // suppllied by user
1594              //
1595              G4double theDecayTime = GetDecayTime();
1596             
1597              G4int nbin = GetDecayTimeBin(theDecayTime);
1598             
1599              // claculate the first part of the weight function
1600             
1601              G4double weight1 =1./DProfile[nbin-1] 
1602                *(DBin[nbin]-DBin[nbin-1])
1603                /NSplit;
1604              if (nbin > 1) {
1605                weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
1606                  *(DBin[nbin]-DBin[nbin-1])
1607                  /NSplit;}
1608              // it should be calculated in seconds
1609              weight1 /= s ;
1610             
1611              // a temprary products buffer and its contents is transfered to
1612              // the products at the end of the loop
1613              //
1614              G4DecayProducts *tempprods = 0;
1615             
1616              // calculate the decay rate of the isotope
1617              // one need to fold the the source bias function with the decaytime
1618              //
1619              decayRate = 0.;
1620              for ( j = 0; j < PT.size(); j++){
1621                taotime = GetTaoTime(theDecayTime,PT[j]);
1622                decayRate -= PR[j] * taotime;
1623              }
1624             
1625              // decayRatehe radioactivity of isotope (PZ,PA,PE) at the
1626              // time 'theDecayTime'
1627              // it will be used to calculate the statistical weight of the
1628              // decay products of this isotope
1629             
1630             
1631              //
1632              // now calculate the statistical weight
1633              //
1634             
1635              G4double weight = weight1*decayRate; 
1636              // decay the isotope
1637              theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
1638              parentNucleus = theIonTable->GetIon(PZ,PA,PE);
1639             
1640              // decide whther to apply branching ratio bias or not
1641              //
1642              if (BRBias){
1643                G4DecayTable *theDecayTable = parentNucleus->GetDecayTable();
1644                ndecaych = G4int(theDecayTable->entries()*G4UniformRand());
1645                G4VDecayChannel *theDecayChannel = theDecayTable->GetDecayChannel(ndecaych);
1646                if (theDecayChannel == 0)
1647                  {
1648                    // Decay channel not found.
1649#ifdef G4VERBOSE
1650                    if (GetVerboseLevel()>0)
1651                      {
1652                        G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel";
1653                        G4cerr <<G4endl;
1654                        theDecayTable ->DumpInfo();
1655                      }
1656#endif
1657                  }
1658                else
1659                  {
1660                    // A decay channel has been identified, so execute the DecayIt.
1661                    G4double tempmass = parentNucleus->GetPDGMass();     
1662                    tempprods = theDecayChannel->DecayIt(tempmass);
1663                    weight *= (theDecayChannel->GetBR())*(theDecayTable->entries());
1664                  }
1665              }
1666              else {
1667                tempprods = DoDecay(*parentNucleus);
1668              }
1669              //
1670              // save the secondaries for buffers
1671              //
1672              numberOfSecondaries = tempprods->entries();
1673              currentTime = finalGlobalTime + theDecayTime;
1674              for (index=0; index < numberOfSecondaries; index++) 
1675                {
1676                  asecondaryparticle = tempprods->PopProducts();
1677                  if (asecondaryparticle->GetDefinition()->GetBaryonNumber() < 5){
1678                    pw.push_back(weight);
1679                    ptime.push_back(currentTime);
1680                    secondaryparticles.push_back(asecondaryparticle);
1681                  }
1682                }
1683              //
1684              delete tempprods;
1685             
1686              //end of i loop
1687            }
1688           
1689            // end of n loop
1690          } 
1691        // now deal with the secondaries in the two stl containers
1692        // and submmit them back to the tracking manager
1693        //
1694        totalNumberOfSecondaries = pw.size();
1695        fParticleChangeForRadDecay.SetNumberOfSecondaries(totalNumberOfSecondaries);
1696        for (index=0; index < totalNumberOfSecondaries; index++) 
1697          { 
1698            G4Track* secondary = new G4Track(
1699                                             secondaryparticles[index], ptime[index], currentPosition);
1700            secondary->SetGoodForTrackingFlag();           
1701                        secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
1702            secondary->SetWeight(pw[index]);       
1703            fParticleChangeForRadDecay.AddSecondary(secondary); 
1704          }
1705        //
1706        // make sure the original track is set to stop and its kinematic energy collected
1707        //
1708        //theTrack.SetTrackStatus(fStopButAlive);
1709        //energyDeposit += theParticle->GetKineticEnergy();
1710       
1711      }
1712   
1713      //
1714      // Kill the parent particle.
1715      //
1716      fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
1717      fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(energyDeposit);
1718      //
1719      fParticleChangeForRadDecay.ProposeGlobalTime( finalGlobalTime );
1720      //
1721      // Reset NumberOfInteractionLengthLeft.
1722      //
1723      ClearNumberOfInteractionLengthLeft();
1724     
1725      return &fParticleChangeForRadDecay ;
1726    }
1727} 
1728
1729////////////////////////////////////////////////////////////////////////////////
1730//
1731//
1732// DoDecay
1733//
1734G4DecayProducts* G4RadioactiveDecay::DoDecay(  G4ParticleDefinition& theParticleDef )
1735{
1736  G4DecayProducts *products = 0;
1737  //
1738  //
1739  // follow the decaytable and generate the secondaries...
1740  //
1741  //
1742#ifdef G4VERBOSE
1743  if (GetVerboseLevel()>0)
1744    {
1745      G4cout<<"Begin of DoDecay..."<<G4endl;
1746    }
1747#endif
1748  G4DecayTable *theDecayTable = theParticleDef.GetDecayTable();
1749  //
1750  // Choose a decay channel.
1751  //
1752#ifdef G4VERBOSE
1753  if (GetVerboseLevel()>0)
1754    {
1755      G4cout <<"Selecte a channel..."<<G4endl;
1756    }
1757#endif
1758  G4VDecayChannel *theDecayChannel = theDecayTable->SelectADecayChannel();
1759  if (theDecayChannel == 0)
1760    {
1761      // Decay channel not found.
1762      //
1763      G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel";
1764      G4cerr <<G4endl;
1765      theDecayTable ->DumpInfo();
1766    }
1767      else
1768    {
1769      //
1770      // A decay channel has been identified, so execute the DecayIt.
1771      //
1772#ifdef G4VERBOSE
1773      if (GetVerboseLevel()>1)
1774        {
1775          G4cerr <<"G4RadioactiveDecay::DoIt : selected decay channel  addr:";
1776          G4cerr <<theDecayChannel <<G4endl;
1777        }
1778#endif
1779     
1780      G4double tempmass = theParticleDef.GetPDGMass();
1781      //
1782     
1783      products = theDecayChannel->DecayIt(tempmass);
1784     
1785    }
1786  return products;
1787
1788}
1789
1790
1791
1792
1793
1794
1795
1796
1797
Note: See TracBrowser for help on using the repository browser.