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

Last change on this file since 1337 was 1315, checked in by garnier, 15 years ago

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