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

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

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