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

Last change on this file since 1245 was 962, checked in by garnier, 15 years ago

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