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

Last change on this file since 1036 was 962, checked in by garnier, 17 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.