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

Last change on this file since 836 was 819, checked in by garnier, 16 years ago

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