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

Last change on this file since 1345 was 1340, checked in by garnier, 15 years ago

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