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

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

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 53.4 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 = 1e-6*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                    theDecayTable->Insert(aBetaMinusChannel);
775                    modeSumBR[1] += b;
776                    delete[] pdf;
777                    delete aBetaFermiFunction;
778                  }
779                }
780                break;
781              case BetaPlus:
782                //
783                //
784                // Decay mode is beta+.
785                //
786                if (modeFirstRecord[2])
787                  {modeFirstRecord[2] = false; modeTotalBR[2] = b;}
788                else {
789                  //  e0 = c*MeV/0.511;
790                  // bug fix (#662) (flei, 22/09/2004)
791                  // need to test e0 as there are some data files (e.g. z67.a162) which have entries for beta+
792                  // with Q < 2Me
793                  //
794                  e0 = c*MeV/0.511 -2.;
795                  if (e0 > 0.) {
796                    G4BetaFermiFunction* aBetaFermiFunction = new G4BetaFermiFunction (A, -(Z-1));
797                               
798                    n = aBetaFermiFunction->GetFFN(e0);
799                               
800                    // now to work out the histogram and initialise the random generator
801                    G4int npti = 100;                           
802                    G4double* pdf = new G4double[npti];
803                    G4int ptn;
804                    G4double g,e,ee,f;
805                    ee = e0+1.;
806                    for (ptn=0; ptn<npti; ptn++) {
807                      // e =e0*(ptn+1.)/102.;
808                      // bug fix (#662) (flei, 22/09/2004)
809                      e =e0*(ptn+0.5)/100.;
810                      g = e+1.;
811                      f = std::sqrt(g*g-1)*(ee-g)*(ee-g)*g;
812                      pdf[ptn] = f*aBetaFermiFunction->GetFF(e);
813                    }
814                    RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti); 
815                    G4BetaPlusDecayChannel *aBetaPlusChannel = new 
816                      G4BetaPlusDecayChannel (GetVerboseLevel(), &theParentNucleus,
817                                              b, c*MeV, a*MeV, n, FBeta, aRandomEnergy);
818                    theDecayTable->Insert(aBetaPlusChannel);
819                    modeSumBR[2] += b;
820                               
821                    delete[] pdf;
822                    delete aBetaFermiFunction;       
823                  }
824                }
825                break;
826              case KshellEC:
827                //
828                //
829                // Decay mode is K-electron capture.
830                //
831                if (modeFirstRecord[3])
832                  {modeFirstRecord[3] = false; modeTotalBR[3] = b;}
833                else {
834                  G4KshellECDecayChannel *aKECChannel = new
835                    G4KshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
836                                            b, c*MeV, a*MeV);
837                  aKECChannel->SetICM(applyICM);
838                  aKECChannel->SetARM(applyARM);
839                  aKECChannel->SetHLThreshold(halflifethreshold);
840                  theDecayTable->Insert(aKECChannel);
841                  //delete aKECChannel;
842                  modeSumBR[3] += b;
843                }
844                break;
845              case LshellEC:
846                //
847                //
848                // Decay mode is L-electron capture.
849                //
850                if (modeFirstRecord[4])
851                  {modeFirstRecord[4] = false; modeTotalBR[4] = b;}
852                else {
853                  G4LshellECDecayChannel *aLECChannel = new
854                    G4LshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
855                                            b, c*MeV, a*MeV);
856                  aLECChannel->SetICM(applyICM);
857                  aLECChannel->SetARM(applyARM);
858                  aLECChannel->SetHLThreshold(halflifethreshold);
859                  theDecayTable->Insert(aLECChannel);
860                  //delete aLECChannel;
861                  modeSumBR[4] += b;
862                }
863                break;
864              case MshellEC:
865                //
866                //
867                // Decay mode is M-electron capture. In this implementation it is added to L-shell case
868                //
869                if (modeFirstRecord[5])
870                  {modeFirstRecord[5] = false; modeTotalBR[5] = b;}
871                else {
872                  G4MshellECDecayChannel *aMECChannel = new
873                    G4MshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
874                                            b, c*MeV, a*MeV);
875                  aMECChannel->SetICM(applyICM);
876                  aMECChannel->SetARM(applyARM);
877                  aMECChannel->SetHLThreshold(halflifethreshold);
878                  theDecayTable->Insert(aMECChannel);
879                  modeSumBR[5] += b;
880                }
881                break;
882              case Alpha:
883                //
884                //
885                // Decay mode is alpha.
886                //
887                if (modeFirstRecord[6])
888                  {modeFirstRecord[6] = false; modeTotalBR[6] = b;}
889                else {
890                  G4AlphaDecayChannel *anAlphaChannel = new
891                    G4AlphaDecayChannel (GetVerboseLevel(), &theParentNucleus,
892                                         b, c*MeV, a*MeV);
893                  theDecayTable->Insert(anAlphaChannel);
894                  //          delete anAlphaChannel;
895                  modeSumBR[6] += b;
896                }
897                break;
898              case ERROR:
899              default:
900                G4Exception("G4RadioactiveDecay::LoadDecayTable()", "601",
901                            FatalException, "Error in decay mode selection");
902               
903              }
904            }
905          }
906        }
907      }
908    }   
909    //
910    //
911    // Go through the decay table and make sure that the branching ratios are
912    // correctly normalised.
913    //
914    G4VDecayChannel       *theChannel             = 0;
915    G4NuclearDecayChannel *theNuclearDecayChannel = 0;
916    G4String mode                     = "";
917    G4int j                           = 0;
918    G4double theBR                    = 0.0;
919    for (i=0; i<theDecayTable->entries(); i++) {
920      theChannel             = theDecayTable->GetDecayChannel(i);
921      theNuclearDecayChannel = static_cast<G4NuclearDecayChannel *>(theChannel);
922      theDecayMode           = theNuclearDecayChannel->GetDecayMode();
923      j          = 0;
924      if (theDecayMode != IT) {
925        theBR = theChannel->GetBR();
926        theChannel->SetBR(theBR*modeTotalBR[theDecayMode]/modeSumBR[theDecayMode]);
927      }
928    } 
929  }             
930  DecaySchemeFile.close();             
931  if (!found && E > 0.) {
932    // cases where IT cascade for exited isotopes without entry in RDM database
933    // Decay mode is isomeric transition.
934    //
935    G4ITDecayChannel *anITChannel = new G4ITDecayChannel
936      (GetVerboseLevel(), (const G4Ions*) &theParentNucleus, 1);
937    anITChannel->SetICM(applyICM);
938    anITChannel->SetARM(applyARM);
939    anITChannel->SetHLThreshold(halflifethreshold);
940    theDecayTable->Insert(anITChannel);
941  } 
942  if (!theDecayTable) {
943    //
944    // There is no radioactive decay data for this nucleus.  Return a null
945    // decay table.
946    //
947    G4cerr <<"G4RadoactiveDecay::LoadDecayTable() : cannot find ion radioactive decay file " <<G4endl;
948    theDecayTable = 0;
949    return theDecayTable;
950  }     
951  if (theDecayTable && GetVerboseLevel()>1)
952    {
953      G4cout <<"G4RadioactiveDecay::LoadDecayTable()" << G4endl;
954      G4cout << "  No. of  entries: "<< theDecayTable->entries() <<G4endl;
955      theDecayTable ->DumpInfo();
956    }
957
958  return theDecayTable;
959}
960
961////////////////////////////////////////////////////////////////////////
962//
963//
964void G4RadioactiveDecay::SetDecayRate(G4int theZ, G4int theA, G4double theE, 
965                                       G4int theG, std::vector<G4double> theRates, 
966                                       std::vector<G4double> theTaos)
967{ 
968  //fill the decay rate vector
969  theDecayRate.SetZ(theZ);
970  theDecayRate.SetA(theA);
971  theDecayRate.SetE(theE);
972  theDecayRate.SetGeneration(theG);
973  theDecayRate.SetDecayRateC(theRates);
974  theDecayRate.SetTaos(theTaos);
975}
976//////////////////////////////////////////////////////////////////////////
977//
978//
979void G4RadioactiveDecay::AddDecayRateTable(const G4ParticleDefinition &theParentNucleus)
980{
981  // 1) To calculate all the coefficiecies required to derive the radioactivities for all
982  // progeny of theParentNucleus
983  //
984  // 2) Add the coefficiencies to the decay rate table vector
985  //
986 
987  //
988  // Create and initialise variables used in the method.
989  //
990
991  theDecayRateVector.clear();
992 
993  G4int nGeneration = 0;
994  std::vector<G4double> rates;
995  std::vector<G4double> taos;
996 
997  // start rate is -1.
998  rates.push_back(-1.);
999  //
1000  //
1001  G4int A = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
1002  G4int Z = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
1003  G4double E = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
1004  G4double tao = theParentNucleus.GetPDGLifeTime();
1005  taos.push_back(tao);
1006  G4int nEntry = 0;
1007 
1008  //fill the decay rate with the intial isotope data
1009  SetDecayRate(Z,A,E,nGeneration,rates,taos);
1010
1011  // store the decay rate in decay rate vector
1012  theDecayRateVector.push_back(theDecayRate);
1013  nEntry++;
1014
1015  // now start treating the sencondary generations..
1016
1017  G4bool stable = false;
1018  G4int i;
1019  G4int j;
1020  G4VDecayChannel       *theChannel             = 0;
1021  G4NuclearDecayChannel *theNuclearDecayChannel = 0;
1022  G4ITDecayChannel *theITChannel = 0;
1023  G4BetaMinusDecayChannel *theBetaMinusChannel = 0;
1024  G4BetaPlusDecayChannel *theBetaPlusChannel = 0;
1025  G4AlphaDecayChannel *theAlphaChannel = 0;
1026  G4RadioactiveDecayMode theDecayMode;
1027  //  G4NuclearLevelManager levelManager;
1028  //const G4NuclearLevel* level;
1029  G4double theBR = 0.0;
1030  G4int AP = 0;
1031  G4int ZP = 0;
1032  G4int AD = 0;
1033  G4int ZD = 0;
1034  G4double EP = 0.;
1035  std::vector<G4double> TP;
1036  std::vector<G4double> RP;
1037  G4ParticleDefinition *theDaughterNucleus;
1038  G4double daughterExcitation;
1039  G4ParticleDefinition *aParentNucleus;
1040  G4IonTable* theIonTable;
1041  G4DecayTable *aTempDecayTable;
1042  G4double theRate;
1043  G4double TaoPlus;
1044  G4int nS = 0;
1045  G4int nT = nEntry;
1046  G4double brs[7];
1047  //
1048  theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
1049 
1050  while (!stable) {
1051    nGeneration++;
1052    for (j = nS; j< nT; j++) {
1053      ZP = theDecayRateVector[j].GetZ();
1054      AP = theDecayRateVector[j].GetA();
1055      EP = theDecayRateVector[j].GetE();
1056      RP = theDecayRateVector[j].GetDecayRateC();
1057      TP = theDecayRateVector[j].GetTaos();     
1058      if (GetVerboseLevel()>0){
1059        G4cout <<"G4RadioactiveDecay::AddDecayRateTable : "
1060               << " daughters of ("<< ZP <<", "<<AP<<", "
1061               << EP <<") "
1062               << " are being calculated. "       
1063               <<" generation = "
1064               << nGeneration << G4endl;
1065      }
1066      aParentNucleus = theIonTable->GetIon(ZP,AP,EP);
1067      if (!IsLoaded(*aParentNucleus)){
1068        aParentNucleus->SetDecayTable(LoadDecayTable(*aParentNucleus));
1069      }
1070      aTempDecayTable = aParentNucleus->GetDecayTable();
1071      //
1072      //
1073      // Go through the decay table and to combine the same decay channels
1074      //
1075      for (i=0; i< 7; i++) brs[i] = 0.0;
1076     
1077      G4DecayTable *theDecayTable = new G4DecayTable();
1078     
1079      for (i=0; i<aTempDecayTable->entries(); i++) {
1080        theChannel             = aTempDecayTable->GetDecayChannel(i);
1081        theNuclearDecayChannel = static_cast<G4NuclearDecayChannel *>(theChannel);
1082        theDecayMode           = theNuclearDecayChannel->GetDecayMode();
1083        daughterExcitation = theNuclearDecayChannel->GetDaughterExcitation ();
1084        theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus () ;
1085        AD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
1086        ZD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber(); 
1087        G4NuclearLevelManager * levelManager = G4NuclearLevelStore::GetInstance()->GetManager (ZD, AD);
1088        if ( levelManager->NumberOfLevels() ) {
1089          const G4NuclearLevel* level = levelManager->NearestLevel (daughterExcitation);
1090
1091          if (std::abs(daughterExcitation - level->Energy()) < levelTolerance) {
1092           
1093            // Level hafe life is in ns and the gate as 1 micros by default
1094            if ( theDecayMode == 0 && level->HalfLife()*ns >= halflifethreshold ){   
1095              // some further though may needed here
1096              theDecayTable->Insert(theChannel);
1097            } 
1098            else{
1099              brs[theDecayMode] += theChannel->GetBR();
1100            }
1101          }
1102          else {
1103            brs[theDecayMode] += theChannel->GetBR();
1104          }
1105        }
1106        else{
1107          brs[theDecayMode] += theChannel->GetBR();
1108        }
1109      }     
1110      brs[2] = brs[2]+brs[3]+brs[4]+brs[5];
1111      brs[3] = brs[4] =brs[5] =  0.0;
1112      for (i= 0; i<7; i++){
1113        if (brs[i] > 0) {
1114          switch ( i ) {
1115          case 0:
1116            //
1117            //
1118            // Decay mode is isomeric transition.
1119            //
1120           
1121            theITChannel =  new G4ITDecayChannel
1122              (0, (const G4Ions*) aParentNucleus, brs[0]);
1123           
1124            theDecayTable->Insert(theITChannel);
1125            break;
1126           
1127          case 1:
1128            //
1129            //
1130            // Decay mode is beta-.
1131            //
1132            theBetaMinusChannel = new G4BetaMinusDecayChannel (0, aParentNucleus,
1133                                                               brs[1], 0.*MeV, 0.*MeV, 1, false, 0);
1134            theDecayTable->Insert(theBetaMinusChannel);
1135           
1136            break;
1137           
1138          case 2:
1139            //
1140            //
1141            // Decay mode is beta+ + EC.
1142            //
1143            theBetaPlusChannel = new G4BetaPlusDecayChannel (GetVerboseLevel(), aParentNucleus,
1144                                                             brs[2], 0.*MeV, 0.*MeV, 1, false, 0);
1145            theDecayTable->Insert(theBetaPlusChannel);
1146            break;                   
1147           
1148          case 6:
1149            //
1150            //
1151            // Decay mode is alpha.
1152            //
1153            theAlphaChannel = new G4AlphaDecayChannel (GetVerboseLevel(), aParentNucleus,
1154                                                       brs[6], 0.*MeV, 0.*MeV);
1155            theDecayTable->Insert(theAlphaChannel);
1156            break;
1157           
1158          default:
1159            break;
1160          }
1161        }
1162      }
1163       
1164      //
1165      // loop over all braches in theDecayTable
1166      //
1167      for ( i=0; i<theDecayTable->entries(); i++){
1168        theChannel             = theDecayTable->GetDecayChannel(i);
1169        theNuclearDecayChannel = static_cast<G4NuclearDecayChannel *>(theChannel);
1170        theBR = theChannel->GetBR();
1171        theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus();
1172       
1173        //
1174        // now test if the daughterNucleus is a valid one
1175        //
1176        if (IsApplicable(*theDaughterNucleus) && theBR
1177            && aParentNucleus != theDaughterNucleus ) {
1178          // need to make sure daugher has decaytable
1179          if (!IsLoaded(*theDaughterNucleus)){
1180            theDaughterNucleus->SetDecayTable(LoadDecayTable(*theDaughterNucleus));
1181          }
1182          if (theDaughterNucleus->GetDecayTable()->entries() ) {
1183            //
1184            A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
1185            Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
1186            E = ((const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy();
1187         
1188            TaoPlus = theDaughterNucleus->GetPDGLifeTime();
1189            //          cout << TaoPlus <<G4endl;
1190            if (TaoPlus > 0.) {
1191              // first set the taos, one simply need to add to the parent ones
1192              taos.clear();
1193              taos = TP;
1194              taos.push_back(TaoPlus);
1195              // now calculate the coefficiencies
1196              //
1197              // they are in two parts, first the les than n ones
1198              rates.clear();
1199              size_t k;
1200              for (k = 0; k < RP.size(); k++){
1201                theRate = TP[k]/(TP[k]-TaoPlus) * theBR * RP[k];
1202                rates.push_back(theRate);
1203              }
1204              //
1205              // the sencond part: the n:n coefficiency
1206              theRate = 0.;
1207              for (k = 0; k < RP.size(); k++){
1208                theRate -=TaoPlus/(TP[k]-TaoPlus) * theBR * RP[k];
1209              }
1210              rates.push_back(theRate);               
1211              SetDecayRate (Z,A,E,nGeneration,rates,taos);
1212              theDecayRateVector.push_back(theDecayRate);
1213              nEntry++;
1214            } 
1215          }
1216        } 
1217        // end of testing daughter nucleus
1218      }
1219      // end of i loop( the branches)
1220    }
1221    //end of for j loop
1222    nS = nT;
1223    nT = nEntry;
1224    if (nS == nT) stable = true;
1225  }
1226 
1227  //end of while loop
1228  // the calculation completed here
1229 
1230 
1231  // fill the first part of the decay rate table
1232  // which is the name of the original particle (isotope)
1233  //
1234  theDecayRateTable.SetIonName(theParentNucleus.GetParticleName()); 
1235  //
1236  //
1237  // now fill the decay table with the newly completed decay rate vector
1238  //
1239 
1240  theDecayRateTable.SetItsRates(theDecayRateVector);
1241 
1242  //
1243  // finally add the decayratetable to the tablevector
1244  //
1245  theDecayRateTableVector.push_back(theDecayRateTable);
1246}
1247 
1248////////////////////////////////////////////////////////////////////////////////
1249//
1250//
1251// SetSourceTimeProfile
1252//
1253//  read in the source time profile function (histogram)
1254//
1255
1256void G4RadioactiveDecay::SetSourceTimeProfile(G4String filename)
1257{
1258  std::ifstream infile ( filename, std::ios::in );
1259  if ( !infile ) G4Exception(__FILE__, G4inttostring(__LINE__), FatalException,  "Unable to open source data file" );
1260 
1261  float bin, flux;
1262  NSourceBin = -1;
1263  while (infile >> bin >> flux ) {
1264    NSourceBin++;
1265    if (NSourceBin > 99)  G4Exception(__FILE__, G4inttostring(__LINE__), FatalException,  "input source data file too big (>100 rows)" );
1266    SBin[NSourceBin] = bin * s;
1267    SProfile[NSourceBin] = flux;
1268  }
1269  SetAnalogueMonteCarlo(0);
1270#ifdef G4VERBOSE
1271  if (GetVerboseLevel()>1)
1272    {G4cout <<" Source Timeprofile Nbin = " << NSourceBin <<G4endl;}
1273#endif
1274}
1275
1276////////////////////////////////////////////////////////////////////////////////
1277//
1278//
1279// SetDecayBiasProfile
1280//
1281// read in the decay bias scheme function (histogram)
1282//
1283void G4RadioactiveDecay::SetDecayBias(G4String filename)
1284{
1285  std::ifstream infile ( filename, std::ios::in);
1286  if ( !infile ) G4Exception(__FILE__, G4inttostring(__LINE__), FatalException,  "Unable to open bias data file" );
1287 
1288  float bin, flux;
1289  NDecayBin = -1;
1290  while (infile >> bin >> flux ) {
1291    NDecayBin++;
1292    if (NDecayBin > 99)  G4Exception(__FILE__, G4inttostring(__LINE__), FatalException,  "input bias data file too big (>100 rows)" );
1293    DBin[NDecayBin] = bin * s;
1294    DProfile[NDecayBin] = flux;
1295  }
1296  G4int i ;
1297  for ( i = 1; i<= NDecayBin; i++) DProfile[i] += DProfile[i-1];
1298  for ( i = 0; i<= NDecayBin; i++) DProfile[i] /= DProfile[NDecayBin];
1299  SetAnalogueMonteCarlo(0);
1300#ifdef G4VERBOSE
1301  if (GetVerboseLevel()>1)
1302    {G4cout <<" Decay Bias Profile  Nbin = " << NDecayBin <<G4endl;}
1303#endif
1304}
1305
1306////////////////////////////////////////////////////////////////////////////////
1307//
1308//
1309// DecayIt
1310//
1311G4VParticleChange* G4RadioactiveDecay::DecayIt(const G4Track& theTrack, const G4Step& )
1312{
1313  //
1314  // Initialize the G4ParticleChange object. Get the particle details and the
1315  // decay table.
1316  //
1317  fParticleChangeForRadDecay.Initialize(theTrack);
1318  const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
1319  G4ParticleDefinition *theParticleDef = theParticle->GetDefinition();
1320
1321  // First check whether RDM applies to the current logical volume
1322  //
1323  if(!std::binary_search(ValidVolumes.begin(),
1324                    ValidVolumes.end(), 
1325                    theTrack.GetVolume()->GetLogicalVolume()->GetName()))
1326    {
1327#ifdef G4VERBOSE
1328      if (GetVerboseLevel()>0)
1329        {
1330          G4cout <<"G4RadioactiveDecay::DecayIt : "
1331                 << theTrack.GetVolume()->GetLogicalVolume()->GetName()
1332                 << " is not selected for the RDM"<< G4endl;
1333          G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
1334          G4cout << " The Valid volumes are " << G4endl;
1335          for (size_t i = 0; i< ValidVolumes.size(); i++)
1336            G4cout << ValidVolumes[i] << G4endl;
1337        }
1338#endif
1339      fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1340      //
1341      //
1342      // Kill the parent particle.
1343      //
1344      fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
1345      fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1346      ClearNumberOfInteractionLengthLeft();
1347      return &fParticleChangeForRadDecay;
1348    }
1349   
1350  // now check is the particle is valid for RDM
1351  //
1352  if (!(IsApplicable(*theParticleDef)))
1353    { 
1354      //
1355      // The particle is not a Ion or outside the nucleuslimits for decay
1356      //
1357#ifdef G4VERBOSE
1358      if (GetVerboseLevel()>0)
1359        {
1360          G4cerr <<"G4RadioactiveDecay::DecayIt : "
1361                 <<theParticleDef->GetParticleName() 
1362                 << " is not a valid nucleus for the RDM"<< G4endl;
1363        }
1364#endif
1365      fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1366      //
1367      //
1368      // Kill the parent particle.
1369      //
1370      fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
1371      fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1372      ClearNumberOfInteractionLengthLeft();
1373      return &fParticleChangeForRadDecay;
1374    }
1375 
1376  if (!IsLoaded(*theParticleDef))
1377    {
1378      theParticleDef->SetDecayTable(LoadDecayTable(*theParticleDef));
1379    }
1380  G4DecayTable *theDecayTable = theParticleDef->GetDecayTable();
1381 
1382  if  (theDecayTable == 0 || theDecayTable->entries() == 0 )
1383    {
1384      //
1385      //
1386      // There are no data in the decay table.  Set the particle change parameters
1387      // to indicate this.
1388      //
1389#ifdef G4VERBOSE
1390      if (GetVerboseLevel()>0)
1391        {
1392          G4cerr <<"G4RadioactiveDecay::DecayIt : decay table not defined  for ";
1393          G4cerr <<theParticleDef->GetParticleName() <<G4endl;
1394        }
1395#endif
1396      fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1397      //
1398      //
1399      // Kill the parent particle.
1400      //
1401      fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
1402      fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1403      ClearNumberOfInteractionLengthLeft();
1404      return &fParticleChangeForRadDecay;
1405    }
1406  else 
1407    { 
1408      //
1409      // now try to  decay it
1410      //
1411      G4double energyDeposit = 0.0;
1412      G4double finalGlobalTime = theTrack.GetGlobalTime();
1413      G4int index;
1414      G4ThreeVector currentPosition;
1415      currentPosition = theTrack.GetPosition();
1416     
1417      // check whether use Analogue or VR implementation
1418      //
1419      if (AnalogueMC){
1420        //
1421        // Aanalogue MC
1422#ifdef G4VERBOSE
1423        if (GetVerboseLevel()>0)
1424          {
1425            G4cout <<"DecayIt:  Analogue MC version " << G4endl;
1426          }
1427#endif
1428        //
1429        G4DecayProducts *products = DoDecay(*theParticleDef);
1430        //
1431        // check if the product is the same as input and kill the track if necessary to prevent infinite loop
1432        // (11/05/10, F.Lei)
1433        //
1434        if ( products->entries() == 1) {
1435          fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1436          fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
1437          fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1438          ClearNumberOfInteractionLengthLeft();
1439          return &fParticleChangeForRadDecay;
1440        }       
1441        //
1442        // Get parent particle information and boost the decay products to the
1443        // laboratory frame based on this information.
1444        //
1445        G4double ParentEnergy = theParticle->GetTotalEnergy();
1446        G4ThreeVector ParentDirection(theParticle->GetMomentumDirection());
1447       
1448        if (theTrack.GetTrackStatus() == fStopButAlive )
1449          {
1450            //
1451            //
1452            // The particle is decayed at rest.
1453            //
1454            // since the time is still for rest particle in G4 we need to add the additional
1455            // time lapsed between the particle come to rest and the actual decay. This time
1456            // is simply sampled with the mean-life of the particle.
1457            // but we need to protect the case PDGTime < 0. (F.Lei 11/05/10)
1458            G4double temptime =  -std::log( G4UniformRand()) * theParticleDef->GetPDGLifeTime();
1459            if (temptime <0.) temptime =0.; 
1460            finalGlobalTime += temptime ;
1461            energyDeposit += theParticle->GetKineticEnergy();
1462          }
1463        else
1464          {
1465            //
1466            //
1467            // The particle is decayed in flight (PostStep case).
1468            //
1469            products->Boost( ParentEnergy, ParentDirection);
1470          }
1471        //
1472        //
1473        // Add products in theParticleChangeForRadDecay.
1474        //
1475        G4int numberOfSecondaries = products->entries();
1476        fParticleChangeForRadDecay.SetNumberOfSecondaries(numberOfSecondaries);
1477#ifdef G4VERBOSE
1478        if (GetVerboseLevel()>1) {
1479          G4cout <<"G4RadioactiveDecay::DecayIt : Decay vertex :";
1480          G4cout <<" Time: " <<finalGlobalTime/ns <<"[ns]";
1481          G4cout <<" X:" <<(theTrack.GetPosition()).x() /cm <<"[cm]";
1482          G4cout <<" Y:" <<(theTrack.GetPosition()).y() /cm <<"[cm]";
1483          G4cout <<" Z:" <<(theTrack.GetPosition()).z() /cm <<"[cm]";
1484          G4cout <<G4endl;
1485          G4cout <<"G4Decay::DecayIt  : decay products in Lab. Frame" <<G4endl;
1486          products->DumpInfo();
1487        }
1488#endif
1489        for (index=0; index < numberOfSecondaries; index++) 
1490          {
1491            G4Track* secondary = new G4Track
1492              (products->PopProducts(), finalGlobalTime, currentPosition);
1493            secondary->SetGoodForTrackingFlag();
1494                        secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
1495            fParticleChangeForRadDecay.AddSecondary(secondary);
1496          }
1497        delete products;
1498        //
1499        // end of analogue MC algarithm
1500        //
1501      }
1502      else {
1503        //
1504        // Varaice Reduction Method
1505        //
1506#ifdef G4VERBOSE
1507        if (GetVerboseLevel()>0)
1508          {
1509            G4cout <<"DecayIt:  Variance Reduction version " << G4endl;
1510          }
1511#endif
1512        if (!IsRateTableReady(*theParticleDef)) {
1513          // if the decayrates are not ready, calculate them and
1514          // add to the rate table vector
1515          AddDecayRateTable(*theParticleDef);
1516        }
1517        //retrieve the rates
1518        GetDecayRateTable(*theParticleDef);
1519        //
1520        // declare some of the variables required in the implementation
1521        //
1522        G4ParticleDefinition* parentNucleus;
1523        G4IonTable *theIonTable;
1524        G4int PZ;
1525        G4int PA;
1526        G4double PE;
1527        std::vector<G4double> PT;
1528        std::vector<G4double> PR;
1529        G4double taotime;
1530        G4double decayRate;
1531       
1532        size_t i;
1533        size_t j;
1534        G4int numberOfSecondaries;
1535        G4int totalNumberOfSecondaries = 0;
1536        G4double currentTime = 0.;
1537        G4int ndecaych;
1538        G4DynamicParticle* asecondaryparticle;
1539        //      G4DecayProducts* products = 0;
1540        std::vector<G4DynamicParticle*> secondaryparticles;
1541        std::vector<G4double> pw;
1542        std::vector<G4double> ptime;
1543        pw.clear();
1544        ptime.clear();
1545        //now apply the nucleus splitting
1546        //
1547        //
1548        for (G4int n = 0; n < NSplit; n++)
1549          {
1550            /*
1551            //
1552            // Get the decay time following the decay probability function
1553            // suppllied by user
1554            //
1555            G4double theDecayTime = GetDecayTime();
1556           
1557            G4int nbin = GetDecayTimeBin(theDecayTime);
1558           
1559            // claculate the first part of the weight function
1560           
1561            G4double weight1 =1./DProfile[nbin-1]
1562              *(DBin[nbin]-DBin[nbin-1])
1563              /NSplit;
1564            if (nbin > 1) {
1565               weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
1566                 *(DBin[nbin]-DBin[nbin-1])
1567                 /NSplit;}
1568            // it should be calculated in seconds
1569            weight1 /= s ;
1570            */
1571            //
1572            // loop over all the possible secondaries of the nucleus
1573            // the first one is itself.
1574            //
1575            for ( i = 0; i<theDecayRateVector.size(); i++){
1576              PZ = theDecayRateVector[i].GetZ();
1577              PA = theDecayRateVector[i].GetA();
1578              PE = theDecayRateVector[i].GetE();
1579              PT = theDecayRateVector[i].GetTaos();
1580              PR = theDecayRateVector[i].GetDecayRateC();
1581
1582              //
1583              // Get the decay time following the decay probability function
1584              // suppllied by user
1585              //
1586              G4double theDecayTime = GetDecayTime();
1587             
1588              G4int nbin = GetDecayTimeBin(theDecayTime);
1589             
1590              // claculate the first part of the weight function
1591             
1592              G4double weight1 =1./DProfile[nbin-1] 
1593                *(DBin[nbin]-DBin[nbin-1])
1594                /NSplit;
1595              if (nbin > 1) {
1596                weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
1597                  *(DBin[nbin]-DBin[nbin-1])
1598                  /NSplit;}
1599              // it should be calculated in seconds
1600              weight1 /= s ;
1601             
1602              // a temprary products buffer and its contents is transfered to
1603              // the products at the end of the loop
1604              //
1605              G4DecayProducts *tempprods = 0;
1606             
1607              // calculate the decay rate of the isotope
1608              // one need to fold the the source bias function with the decaytime
1609              //
1610              decayRate = 0.;
1611              for ( j = 0; j < PT.size(); j++){
1612                taotime = GetTaoTime(theDecayTime,PT[j]);
1613                decayRate -= PR[j] * taotime;
1614              }
1615             
1616              // decayRatehe radioactivity of isotope (PZ,PA,PE) at the
1617              // time 'theDecayTime'
1618              // it will be used to calculate the statistical weight of the
1619              // decay products of this isotope
1620             
1621             
1622              //
1623              // now calculate the statistical weight
1624              //
1625             
1626              G4double weight = weight1*decayRate; 
1627              // decay the isotope
1628              theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
1629              parentNucleus = theIonTable->GetIon(PZ,PA,PE);
1630             
1631              // decide whther to apply branching ratio bias or not
1632              //
1633              if (BRBias){
1634                G4DecayTable *theDecayTable = parentNucleus->GetDecayTable();
1635                ndecaych = G4int(theDecayTable->entries()*G4UniformRand());
1636                G4VDecayChannel *theDecayChannel = theDecayTable->GetDecayChannel(ndecaych);
1637                if (theDecayChannel == 0)
1638                  {
1639                    // Decay channel not found.
1640#ifdef G4VERBOSE
1641                    if (GetVerboseLevel()>0)
1642                      {
1643                        G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel";
1644                        G4cerr <<G4endl;
1645                        theDecayTable ->DumpInfo();
1646                      }
1647#endif
1648                  }
1649                else
1650                  {
1651                    // A decay channel has been identified, so execute the DecayIt.
1652                    G4double tempmass = parentNucleus->GetPDGMass();     
1653                    tempprods = theDecayChannel->DecayIt(tempmass);
1654                    weight *= (theDecayChannel->GetBR())*(theDecayTable->entries());
1655                  }
1656              }
1657              else {
1658                tempprods = DoDecay(*parentNucleus);
1659              }
1660              //
1661              // save the secondaries for buffers
1662              //
1663              numberOfSecondaries = tempprods->entries();
1664              currentTime = finalGlobalTime + theDecayTime;
1665              for (index=0; index < numberOfSecondaries; index++) 
1666                {
1667                  asecondaryparticle = tempprods->PopProducts();
1668                  if (asecondaryparticle->GetDefinition()->GetBaryonNumber() < 5){
1669                    pw.push_back(weight);
1670                    ptime.push_back(currentTime);
1671                    secondaryparticles.push_back(asecondaryparticle);
1672                  }
1673                }
1674              //
1675              delete tempprods;
1676             
1677              //end of i loop
1678            }
1679           
1680            // end of n loop
1681          } 
1682        // now deal with the secondaries in the two stl containers
1683        // and submmit them back to the tracking manager
1684        //
1685        totalNumberOfSecondaries = pw.size();
1686        fParticleChangeForRadDecay.SetNumberOfSecondaries(totalNumberOfSecondaries);
1687        for (index=0; index < totalNumberOfSecondaries; index++) 
1688          { 
1689            G4Track* secondary = new G4Track(
1690                                             secondaryparticles[index], ptime[index], currentPosition);
1691            secondary->SetGoodForTrackingFlag();           
1692                        secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
1693            secondary->SetWeight(pw[index]);       
1694            fParticleChangeForRadDecay.AddSecondary(secondary); 
1695          }
1696        //
1697        // make sure the original track is set to stop and its kinematic energy collected
1698        //
1699        //theTrack.SetTrackStatus(fStopButAlive);
1700        //energyDeposit += theParticle->GetKineticEnergy();
1701       
1702      }
1703   
1704      //
1705      // Kill the parent particle.
1706      //
1707      fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
1708      fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(energyDeposit);
1709      //
1710      fParticleChangeForRadDecay.ProposeGlobalTime( finalGlobalTime );
1711      //
1712      // Reset NumberOfInteractionLengthLeft.
1713      //
1714      ClearNumberOfInteractionLengthLeft();
1715     
1716      return &fParticleChangeForRadDecay ;
1717    }
1718} 
1719
1720////////////////////////////////////////////////////////////////////////////////
1721//
1722//
1723// DoDecay
1724//
1725G4DecayProducts* G4RadioactiveDecay::DoDecay(  G4ParticleDefinition& theParticleDef )
1726{
1727  G4DecayProducts *products = 0;
1728  //
1729  //
1730  // follow the decaytable and generate the secondaries...
1731  //
1732  //
1733#ifdef G4VERBOSE
1734  if (GetVerboseLevel()>0)
1735    {
1736      G4cout<<"Begin of DoDecay..."<<G4endl;
1737    }
1738#endif
1739  G4DecayTable *theDecayTable = theParticleDef.GetDecayTable();
1740  //
1741  // Choose a decay channel.
1742  //
1743#ifdef G4VERBOSE
1744  if (GetVerboseLevel()>0)
1745    {
1746      G4cout <<"Selecte a channel..."<<G4endl;
1747    }
1748#endif
1749  G4VDecayChannel *theDecayChannel = theDecayTable->SelectADecayChannel();
1750  if (theDecayChannel == 0)
1751    {
1752      // Decay channel not found.
1753      //
1754      G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel";
1755      G4cerr <<G4endl;
1756      theDecayTable ->DumpInfo();
1757    }
1758      else
1759    {
1760      //
1761      // A decay channel has been identified, so execute the DecayIt.
1762      //
1763#ifdef G4VERBOSE
1764      if (GetVerboseLevel()>1)
1765        {
1766          G4cerr <<"G4RadioactiveDecay::DoIt : selected decay channel  addr:";
1767          G4cerr <<theDecayChannel <<G4endl;
1768        }
1769#endif
1770     
1771      G4double tempmass = theParticleDef.GetPDGMass();
1772      //
1773     
1774      products = theDecayChannel->DecayIt(tempmass);
1775     
1776    }
1777  return products;
1778
1779}
1780
1781
1782
1783
1784
1785
1786
1787
1788
Note: See TracBrowser for help on using the repository browser.