ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RadioactiveDecayBase.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RadioactiveDecayBase.cc
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 //
27 // //
28 // File: G4RadioactiveDecayBase.cc //
29 // Author: D.H. Wright (SLAC) //
30 // Date: 9 August 2017 //
31 // Description: version the G4RadioactiveDecay process by F. Lei and //
32 // P.R. Truscott with biasing and activation calculations //
33 // removed to a derived class. It performs alpha, beta, //
34 // electron capture and isomeric transition decays of //
35 // radioactive nuclei. //
36 // //
38 
41 
42 #include "G4SystemOfUnits.hh"
43 #include "G4DynamicParticle.hh"
44 #include "G4DecayProducts.hh"
45 #include "G4DecayTable.hh"
47 #include "G4ITDecay.hh"
48 #include "G4BetaDecayType.hh"
49 #include "G4BetaMinusDecay.hh"
50 #include "G4BetaPlusDecay.hh"
51 #include "G4ECDecay.hh"
52 #include "G4AlphaDecay.hh"
53 #include "G4TritonDecay.hh"
54 #include "G4ProtonDecay.hh"
55 #include "G4NeutronDecay.hh"
56 #include "G4SFDecay.hh"
57 #include "G4VDecayChannel.hh"
58 #include "G4NuclearDecay.hh"
60 #include "G4Fragment.hh"
61 #include "G4Ions.hh"
62 #include "G4IonTable.hh"
63 #include "G4BetaDecayType.hh"
64 #include "Randomize.hh"
65 #include "G4LogicalVolumeStore.hh"
66 #include "G4NuclearLevelData.hh"
67 #include "G4DeexPrecoParameters.hh"
68 #include "G4LevelManager.hh"
69 #include "G4ThreeVector.hh"
70 #include "G4Electron.hh"
71 #include "G4Positron.hh"
72 #include "G4Neutron.hh"
73 #include "G4Gamma.hh"
74 #include "G4Alpha.hh"
75 #include "G4Triton.hh"
76 #include "G4Proton.hh"
77 
78 #include "G4HadronicProcessType.hh"
80 #include "G4HadronicException.hh"
81 #include "G4LossTableManager.hh"
82 #include "G4VAtomDeexcitation.hh"
83 #include "G4UAtomicDeexcitation.hh"
84 #include "G4PhotonEvaporation.hh"
85 
86 #include <vector>
87 #include <sstream>
88 #include <algorithm>
89 #include <fstream>
90 
91 using namespace CLHEP;
92 
95 
96 #ifdef G4MULTITHREADED
97 #include "G4AutoLock.hh"
98 G4Mutex G4RadioactiveDecayBase::radioactiveDecayMutex = G4MUTEX_INITIALIZER;
99 DecayTableMap* G4RadioactiveDecayBase::master_dkmap = 0;
100 
101 G4int& G4RadioactiveDecayBase::NumberOfInstances()
102 {
103  static G4int numberOfInstances = 0;
104  return numberOfInstances;
105 }
106 #endif
107 
109  : G4VRestDiscreteProcess(processName, fDecay), isInitialised(false),
110  forceDecayDirection(0.,0.,0.), forceDecayHalfAngle(0.*deg), dirPath(""),
111  verboseLevel(1)
112 {
113 #ifdef G4VERBOSE
114  if (GetVerboseLevel() > 1) {
115  G4cout << "G4RadioactiveDecayBase constructor: processName = " << processName
116  << G4endl;
117  }
118 #endif
119 
121 
124 
125  // Set up photon evaporation for use in G4ITDecay
128  photonEvaporation->SetICM(true);
129 
130  // DHW G4DeexPrecoParameters* deex = G4NuclearLevelData::GetInstance()->GetParameters();
131  // DHW deex->SetCorrelatedGamma(true);
132 
133  // Check data directory
134  char* path_var = std::getenv("G4RADIOACTIVEDATA");
135  if (!path_var) {
136  G4Exception("G4RadioactiveDecay()","HAD_RDM_200",FatalException,
137  "Environment variable G4RADIOACTIVEDATA is not set");
138  } else {
139  dirPath = path_var; // convert to string
140  std::ostringstream os;
141  os << dirPath << "/z1.a3"; // used as a dummy
142  std::ifstream testFile;
143  testFile.open(os.str() );
144  if (!testFile.is_open() )
145  G4Exception("G4RadioactiveDecay()","HAD_RDM_201",FatalException,
146  "Environment variable G4RADIOACTIVEDATA is set, but does not point to correct directory");
147  }
148 
149  // Reset the list of user defined data files
151 
152  // Instantiate the map of decay tables
153 #ifdef G4MULTITHREADED
154  G4AutoLock lk(&G4RadioactiveDecayBase::radioactiveDecayMutex);
155  NumberOfInstances()++;
156  if(!master_dkmap) master_dkmap = new DecayTableMap;
157 #endif
158  dkmap = new DecayTableMap;
159 
160  // Apply default values
161  applyARM = true;
162  applyICM = true; // Always on; keep only for backward compatibility
163 
164  // RDM applies to all logical volumes by default
165  isAllVolumesMode = true;
168 }
169 
170 void G4RadioactiveDecayBase::ProcessDescription(std::ostream& outFile) const
171 {
172  outFile << "The radioactive decay process (G4RadioactiveDecay) handles the\n"
173  << "alpha, beta+, beta-, electron capture and isomeric transition\n"
174  << "decays of nuclei (G4GenericIon) with masses A > 4.\n"
175  << "The required half-lives and decay schemes are retrieved from\n"
176  << "the RadioactiveDecay database which was derived from ENSDF.\n";
177 }
178 
179 
181 {
183  delete photonEvaporation;
184  for (DecayTableMap::iterator i = dkmap->begin(); i != dkmap->end(); i++) {
185  delete i->second;
186  }
187  dkmap->clear();
188  delete dkmap;
189 #ifdef G4MULTITHREADED
190  G4AutoLock lk(&G4RadioactiveDecayBase::radioactiveDecayMutex);
191  --NumberOfInstances();
192  if(NumberOfInstances()==0)
193  {
194  for (DecayTableMap::iterator i = master_dkmap->begin(); i != master_dkmap->end(); i++) {
195  delete i->second;
196  }
197  master_dkmap->clear();
198  delete master_dkmap;
199  }
200 #endif
201 }
202 
203 
205 {
206  // All particles other than G4Ions, are rejected by default
207  if (((const G4Ions*)(&aParticle))->GetExcitationEnergy() > 0.) {return true;}
208  if (aParticle.GetParticleName() == "GenericIon") {
209  return true;
210  } else if (!(aParticle.GetParticleType() == "nucleus")
211  || aParticle.GetPDGLifeTime() < 0. ) {
212  return false;
213  }
214 
215  // Determine whether the nuclide falls into the correct A and Z range
216  G4int A = ((const G4Ions*) (&aParticle))->GetAtomicMass();
217  G4int Z = ((const G4Ions*) (&aParticle))->GetAtomicNumber();
218 
220  {return false;}
221  else if (Z > theNucleusLimits.GetZMax() || Z < theNucleusLimits.GetZMin())
222  {return false;}
223  return true;
224 }
225 
227 {
228  G4String key = aNucleus->GetParticleName();
229  DecayTableMap::iterator table_ptr = dkmap->find(key);
230 
231  G4DecayTable* theDecayTable = 0;
232  if (table_ptr == dkmap->end() ) { // If table not there,
233  theDecayTable = LoadDecayTable(*aNucleus); // load from file and
234  if(theDecayTable) (*dkmap)[key] = theDecayTable; // store in library
235  } else {
236  theDecayTable = table_ptr->second;
237  }
238  return theDecayTable;
239 }
240 
241 
243 {
244  G4LogicalVolumeStore* theLogicalVolumes;
246  theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
247  for (size_t i = 0; i < theLogicalVolumes->size(); i++) {
248  volume = (*theLogicalVolumes)[i];
249  if (volume->GetName() == aVolume) {
250  ValidVolumes.push_back(aVolume);
251  std::sort(ValidVolumes.begin(), ValidVolumes.end());
252  // sort need for performing binary_search
253 
254  if (GetVerboseLevel() > 0)
255  G4cout << " Radioactive decay applied to " << aVolume << G4endl;
256 
257  } else if (i == theLogicalVolumes->size() ) {
259  ed << aVolume << " is not a valid logical volume name. Decay not activated for it."
260  << G4endl;
261  G4Exception("G4RadioactiveDecayBase::SelectAVolume()", "HAD_RDM_300",
262  JustWarning, ed);
263  }
264  }
265 }
266 
267 
269 {
270  G4LogicalVolumeStore* theLogicalVolumes;
272  theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
273  for (size_t i = 0; i < theLogicalVolumes->size(); i++) {
274  volume = (*theLogicalVolumes)[i];
275  if (volume->GetName() == aVolume) {
276  std::vector<G4String>::iterator location;
277  location = std::find(ValidVolumes.begin(),ValidVolumes.end(),aVolume);
278  if (location != ValidVolumes.end() ) {
279  ValidVolumes.erase(location);
280  std::sort(ValidVolumes.begin(), ValidVolumes.end());
281  isAllVolumesMode = false;
282  if (GetVerboseLevel() > 0)
283  G4cout << " G4RadioactiveDecay::DeselectAVolume: " << aVolume
284  << " is removed from list " << G4endl;
285  } else {
287  ed << aVolume << " is not in the list. No action taken." << G4endl;
288  G4Exception("G4RadioactiveDecayBase::DeselectAVolume()", "HAD_RDM_300",
289  JustWarning, ed);
290  }
291  } else if (i == theLogicalVolumes->size()) {
293  ed << aVolume << " is not a valid logical volume name. No action taken."
294  << G4endl;
295  G4Exception("G4RadioactiveDecayBase::DeselectAVolume()", "HAD_RDM_300",
296  JustWarning, ed);
297  }
298  }
299 }
300 
301 
303 {
304  G4LogicalVolumeStore* theLogicalVolumes;
306  theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
307  ValidVolumes.clear();
308 #ifdef G4VERBOSE
309  if (GetVerboseLevel()>1)
310  G4cout << " RDM Applies to all Volumes" << G4endl;
311 #endif
312  for (size_t i = 0; i < theLogicalVolumes->size(); i++){
313  volume = (*theLogicalVolumes)[i];
314  ValidVolumes.push_back(volume->GetName());
315 #ifdef G4VERBOSE
316  if (GetVerboseLevel()>1)
317  G4cout << " RDM Applies to Volume " << volume->GetName() << G4endl;
318 #endif
319  }
320  std::sort(ValidVolumes.begin(), ValidVolumes.end());
321  // sort needed in order to allow binary_search
322  isAllVolumesMode=true;
323 }
324 
325 
327 {
328  ValidVolumes.clear();
329  isAllVolumesMode=false;
330 #ifdef G4VERBOSE
331  if (GetVerboseLevel() > 1) G4cout << "RDM removed from all volumes" << G4endl;
332 #endif
333 }
334 
335 
337 // //
338 // GetMeanLifeTime (required by the base class) //
339 // //
341 
344 {
345  G4double meanlife = 0.;
346  const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
347  const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
348  G4double theLife = theParticleDef->GetPDGLifeTime();
349 #ifdef G4VERBOSE
350  if (GetVerboseLevel() > 2) {
351  G4cout << "G4RadioactiveDecay::GetMeanLifeTime() " << G4endl;
352  G4cout << "KineticEnergy: " << theParticle->GetKineticEnergy()/GeV
353  << " GeV, Mass: " << theParticle->GetMass()/GeV
354  << " GeV, Life time: " << theLife/ns << " ns " << G4endl;
355  }
356 #endif
357  if (theParticleDef->GetPDGStable()) {meanlife = DBL_MAX;}
358  else if (theLife < 0.0) {meanlife = DBL_MAX;}
359  else {meanlife = theLife;}
360  // Set meanlife to zero for excited istopes which are not in the
361  // RDM database
362  if (((const G4Ions*)(theParticleDef))->GetExcitationEnergy() > 0. &&
363  meanlife == DBL_MAX) {meanlife = 0.;}
364 #ifdef G4VERBOSE
365  if (GetVerboseLevel() > 2)
366  G4cout << " mean life time: " << meanlife/s << " s " << G4endl;
367 #endif
368 
369  return meanlife;
370 }
371 
373 // //
374 // GetMeanFreePath for decay in flight //
375 // //
377 
380 {
381  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
382  const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
383  G4double tau = aParticleDef->GetPDGLifeTime();
384  G4double aMass = aParticle->GetMass();
385 
386 #ifdef G4VERBOSE
387  if (GetVerboseLevel() > 2) {
388  G4cout << "G4RadioactiveDecay::GetMeanFreePath() " << G4endl;
389  G4cout << " KineticEnergy: " << aParticle->GetKineticEnergy()/GeV
390  << " GeV, Mass: " << aMass/GeV << " GeV, tau: " << tau << " ns "
391  << G4endl;
392  }
393 #endif
394  G4double pathlength = DBL_MAX;
395  if (tau != -1) {
396  // Ion can decay
397 
398  if (tau < -1000.0) {
399  pathlength = DBL_MIN; // nuclide had very short lifetime or wasn't in table
400 
401  } else if (tau < 0.0) {
402  G4cout << aParticleDef->GetParticleName() << " has lifetime " << tau << G4endl;
404  ed << "Ion has negative lifetime " << tau
405  << " but is not stable. Setting mean free path to DBL_MAX" << G4endl;
406  G4Exception("G4RadioactiveDecay::GetMeanFreePath()", "HAD_RDM_011",
407  JustWarning, ed);
408  pathlength = DBL_MAX;
409 
410  } else {
411  // Calculate mean free path
412  G4double betaGamma = aParticle->GetTotalMomentum()/aMass;
413  pathlength = c_light*tau*betaGamma;
414 
415  if (pathlength < DBL_MIN) {
416  pathlength = DBL_MIN;
417 #ifdef G4VERBOSE
418  if (GetVerboseLevel() > 2) {
419  G4cout << "G4Decay::GetMeanFreePath: "
420  << aParticleDef->GetParticleName()
421  << " stops, kinetic energy = "
422  << aParticle->GetKineticEnergy()/keV <<" keV " << G4endl;
423  }
424 #endif
425  }
426  }
427  }
428 
429 #ifdef G4VERBOSE
430  if (GetVerboseLevel() > 2) {
431  G4cout << "mean free path: "<< pathlength/m << " m" << G4endl;
432  }
433 #endif
434  return pathlength;
435 }
436 
438 // //
439 // BuildPhysicsTable - initialization of atomic de-excitation //
440 // //
442 
444 {
445  if (!isInitialised) {
446  isInitialised = true;
447 #ifdef G4VERBOSE
449 #endif
450  }
453 }
454 
456 // //
457 // StreamInfo - stream out parameters //
458 // //
460 
461 void
462 G4RadioactiveDecayBase::StreamInfo(std::ostream& os, const G4String& endline)
463 {
464  G4DeexPrecoParameters* deex =
467 
468  G4int prec = os.precision(5);
469  os << "======================================================================"
470  << endline;
471  os << "====== Radioactive Decay Physics Parameters ======="
472  << endline;
473  os << "======================================================================"
474  << endline;
475  os << "Max life time "
476  << deex->GetMaxLifeTime()/CLHEP::ps << " ps" << endline;
477  os << "Internal e- conversion flag "
478  << deex->GetInternalConversionFlag() << endline;
479  os << "Stored internal conversion coefficients "
480  << deex->StoreICLevelData() << endline;
481  os << "Enable correlated gamma emission "
482  << deex->CorrelatedGamma() << endline;
483  os << "Max 2J for sampling of angular correlations "
484  << deex->GetTwoJMAX() << endline;
485  os << "Atomic de-excitation enabled "
486  << emparam->Fluo() << endline;
487  os << "Auger electron emission enabled "
488  << emparam->Auger() << endline;
489  os << "Auger cascade enabled "
490  << emparam->AugerCascade() << endline;
491  os << "Check EM cuts disabled for atomic de-excitation "
492  << emparam->DeexcitationIgnoreCut() << endline;
493  os << "Use Bearden atomic level energies "
494  << emparam->BeardenFluoDir() << endline;
495  os << "======================================================================"
496  << endline;
497  os.precision(prec);
498 }
499 
500 
502 // //
503 // LoadDecayTable loads the decay scheme from the RadioactiveDecay database //
504 // for the parent nucleus. //
505 // //
507 
510 {
511  // Generate input data file name using Z and A of the parent nucleus
512  // file containing radioactive decay data.
513  G4int A = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
514  G4int Z = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
515 
516  G4double levelEnergy = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
517  G4Ions::G4FloatLevelBase floatingLevel =
518  ((const G4Ions*)(&theParentNucleus))->GetFloatLevelBase();
519 
520 #ifdef G4MULTITHREADED
521  G4AutoLock lk(&G4RadioactiveDecayBase::radioactiveDecayMutex);
522 
523  G4String key = theParentNucleus.GetParticleName();
524  DecayTableMap::iterator master_table_ptr = master_dkmap->find(key);
525 
526  if (master_table_ptr != master_dkmap->end() ) { // If table is there
527  return master_table_ptr->second;
528  }
529 #endif
530 
531  //Check if data have been provided by the user
533 
534  if (file == "") {
535  std::ostringstream os;
536  os << dirPath << "/z" << Z << ".a" << A << '\0';
537  file = os.str();
538  }
539 
540  G4DecayTable* theDecayTable = new G4DecayTable();
541  G4bool found(false); // True if energy level matches one in table
542 
543  std::ifstream DecaySchemeFile;
544  DecaySchemeFile.open(file);
545 
546  if (DecaySchemeFile.good()) {
547  // Initialize variables used for reading in radioactive decay data
548  G4bool floatMatch(false);
549  const G4int nMode = 12;
550  G4double modeTotalBR[nMode] = {0.0};
551  G4double modeSumBR[nMode];
552  for (G4int i = 0; i < nMode; i++) {
553  modeSumBR[i] = 0.0;
554  }
555 
556  char inputChars[120]={' '};
557  G4String inputLine;
558  G4String recordType("");
559  G4String floatingFlag("");
560  G4String daughterFloatFlag("");
561  G4Ions::G4FloatLevelBase daughterFloatLevel;
562  G4RadioactiveDecayMode theDecayMode;
563  G4double decayModeTotal(0.0);
564  G4double parentExcitation(0.0);
565  G4double a(0.0);
566  G4double b(0.0);
567  G4double c(0.0);
568  G4double dummy(0.0);
569  G4BetaDecayType betaType(allowed);
570 
571  // Loop through each data file record until you identify the decay
572  // data relating to the nuclide of concern.
573 
574  G4bool complete(false); // bool insures only one set of values read for any
575  // given parent energy level
576  G4int loop = 0;
577  while (!complete && !DecaySchemeFile.getline(inputChars, 120).eof()) { /* Loop checking, 01.09.2015, D.Wright */
578  loop++;
579  if (loop > 100000) {
580  G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_100",
581  JustWarning, "While loop count exceeded");
582  break;
583  }
584 
585  inputLine = inputChars;
586  inputLine = inputLine.strip(1);
587  if (inputChars[0] != '#' && inputLine.length() != 0) {
588  std::istringstream tmpStream(inputLine);
589 
590  if (inputChars[0] == 'P') {
591  // Nucleus is a parent type. Check excitation level to see if it
592  // matches that of theParentNucleus
593  tmpStream >> recordType >> parentExcitation >> floatingFlag >> dummy;
594  // "dummy" takes the place of half-life
595  // Now read in from ENSDFSTATE in particle category
596 
597  if (found) {
598  complete = true;
599  } else {
600  // Take first level which matches excitation energy regardless of floating level
601  found = (std::abs(parentExcitation*keV - levelEnergy) < levelTolerance);
602  if (floatingLevel != noFloat) {
603  // If floating level specificed, require match of both energy and floating level
604  floatMatch = (floatingLevel == G4Ions::FloatLevelBase(floatingFlag.back()) );
605  if (!floatMatch) found = false;
606  }
607  }
608 
609  } else if (found) {
610  // The right part of the radioactive decay data file has been found. Search
611  // through it to determine the mode of decay of the subsequent records.
612 
613  // Store for later the total decay probability for each decay mode
614  if (inputLine.length() < 72) {
615  tmpStream >> theDecayMode >> dummy >> decayModeTotal;
616  switch (theDecayMode) {
617  case IT:
618  {
619  G4ITDecay* anITChannel = new G4ITDecay(&theParentNucleus, decayModeTotal,
620  0.0, 0.0, photonEvaporation);
621 // anITChannel->SetHLThreshold(halflifethreshold);
622  anITChannel->SetARM(applyARM);
623  theDecayTable->Insert(anITChannel);
624 // anITChannel->DumpNuclearInfo();
625  }
626  break;
627  case BetaMinus:
628  modeTotalBR[1] = decayModeTotal; break;
629  case BetaPlus:
630  modeTotalBR[2] = decayModeTotal; break;
631  case KshellEC:
632  modeTotalBR[3] = decayModeTotal; break;
633  case LshellEC:
634  modeTotalBR[4] = decayModeTotal; break;
635  case MshellEC:
636  modeTotalBR[5] = decayModeTotal; break;
637  case NshellEC:
638  modeTotalBR[6] = decayModeTotal; break;
639  case Alpha:
640  modeTotalBR[7] = decayModeTotal; break;
641  case Proton:
642  modeTotalBR[8] = decayModeTotal; break;
643  case Neutron:
644  modeTotalBR[9] = decayModeTotal; break;
645  case BDProton:
646  break;
647  case BDNeutron:
648  break;
649  case Beta2Minus:
650  break;
651  case Beta2Plus:
652  break;
653  case Proton2:
654  break;
655  case Neutron2:
656  break;
657  case SpFission:
658  modeTotalBR[10] = decayModeTotal; break;
659  case Triton:
660  modeTotalBR[11] = decayModeTotal; break;
661  case RDM_ERROR:
662 
663  default:
664  G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_000",
665  FatalException, "Selected decay mode does not exist");
666  } // switch
667 
668  } else {
669  if (inputLine.length() < 84) {
670  tmpStream >> theDecayMode >> a >> daughterFloatFlag >> b >> c;
671  betaType = allowed;
672  } else {
673  tmpStream >> theDecayMode >> a >> daughterFloatFlag >> b >> c >> betaType;
674  }
675 
676  // Allowed transitions are the default. Forbidden transitions are
677  // indicated in the last column.
678  a /= 1000.;
679  c /= 1000.;
680  b /= 100.;
681  daughterFloatLevel = G4Ions::FloatLevelBase(daughterFloatFlag.back());
682 
683  switch (theDecayMode) {
684  case BetaMinus:
685  {
686  G4BetaMinusDecay* aBetaMinusChannel =
687  new G4BetaMinusDecay(&theParentNucleus, b, c*MeV, a*MeV,
688  daughterFloatLevel, betaType);
689 // aBetaMinusChannel->DumpNuclearInfo();
690 // aBetaMinusChannel->SetHLThreshold(halflifethreshold);
691  theDecayTable->Insert(aBetaMinusChannel);
692  modeSumBR[1] += b;
693  }
694  break;
695 
696  case BetaPlus:
697  {
698  G4BetaPlusDecay* aBetaPlusChannel =
699  new G4BetaPlusDecay(&theParentNucleus, b, c*MeV, a*MeV,
700  daughterFloatLevel, betaType);
701 // aBetaPlusChannel->DumpNuclearInfo();
702 // aBetaPlusChannel->SetHLThreshold(halflifethreshold);
703  theDecayTable->Insert(aBetaPlusChannel);
704  modeSumBR[2] += b;
705  }
706  break;
707 
708  case KshellEC: // K-shell electron capture
709  {
710  G4ECDecay* aKECChannel =
711  new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
712  daughterFloatLevel, KshellEC);
713 // aKECChannel->DumpNuclearInfo();
714 // aKECChannel->SetHLThreshold(halflifethreshold);
715  aKECChannel->SetARM(applyARM);
716  theDecayTable->Insert(aKECChannel);
717  modeSumBR[3] += b;
718  }
719  break;
720 
721  case LshellEC: // L-shell electron capture
722  {
723  G4ECDecay* aLECChannel =
724  new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
725  daughterFloatLevel, LshellEC);
726 // aLECChannel->DumpNuclearInfo();
727 // aLECChannel->SetHLThreshold(halflifethreshold);
728  aLECChannel->SetARM(applyARM);
729  theDecayTable->Insert(aLECChannel);
730  modeSumBR[4] += b;
731  }
732  break;
733 
734  case MshellEC: // M-shell electron capture
735  {
736  G4ECDecay* aMECChannel =
737  new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
738  daughterFloatLevel, MshellEC);
739 // aMECChannel->DumpNuclearInfo();
740 // aMECChannel->SetHLThreshold(halflifethreshold);
741  aMECChannel->SetARM(applyARM);
742  theDecayTable->Insert(aMECChannel);
743  modeSumBR[5] += b;
744  }
745  break;
746 
747  case NshellEC: // N-shell electron capture
748  {
749  G4ECDecay* aNECChannel =
750  new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
751  daughterFloatLevel, NshellEC);
752 // aNECChannel->DumpNuclearInfo();
753 // aNECChannel->SetHLThreshold(halflifethreshold);
754  aNECChannel->SetARM(applyARM);
755  theDecayTable->Insert(aNECChannel);
756  modeSumBR[6] += b;
757  }
758  break;
759 
760  case Alpha:
761  {
762  G4AlphaDecay* anAlphaChannel =
763  new G4AlphaDecay(&theParentNucleus, b, c*MeV, a*MeV,
764  daughterFloatLevel);
765 // anAlphaChannel->DumpNuclearInfo();
766 // anAlphaChannel->SetHLThreshold(halflifethreshold);
767  theDecayTable->Insert(anAlphaChannel);
768  modeSumBR[7] += b;
769  }
770  break;
771 
772  case Proton:
773  {
774  G4ProtonDecay* aProtonChannel =
775  new G4ProtonDecay(&theParentNucleus, b, c*MeV, a*MeV,
776  daughterFloatLevel);
777 // aProtonChannel->DumpNuclearInfo();
778 // aProtonChannel->SetHLThreshold(halflifethreshold);
779  theDecayTable->Insert(aProtonChannel);
780  modeSumBR[8] += b;
781  }
782  break;
783 
784  case Neutron:
785  {
786  G4NeutronDecay* aNeutronChannel =
787  new G4NeutronDecay(&theParentNucleus, b, c*MeV, a*MeV,
788  daughterFloatLevel);
789 // aNeutronChannel->DumpNuclearInfo();
790 // aNeutronChannel->SetHLThreshold(halflifethreshold);
791  theDecayTable->Insert(aNeutronChannel);
792  modeSumBR[9] += b;
793  }
794  break;
795 
796  case BDProton:
797  // Not yet implemented
798  // G4cout << " beta-delayed proton decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
799  break;
800  case BDNeutron:
801  // Not yet implemented
802  // G4cout << " beta-delayed neutron decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
803  break;
804  case Beta2Minus:
805  // Not yet implemented
806  // G4cout << " Double beta- decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
807  break;
808  case Beta2Plus:
809  // Not yet implemented
810  // G4cout << " Double beta+ decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
811  break;
812  case Proton2:
813  // Not yet implemented
814  // G4cout << " Double proton decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
815  break;
816  case Neutron2:
817  // Not yet implemented
818  // G4cout << " Double beta- decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
819  break;
820  case SpFission:
821  {
822  G4SFDecay* aSpontFissChannel =
823 // new G4SFDecay(&theParentNucleus, decayModeTotal, 0.0, 0.0);
824  new G4SFDecay(&theParentNucleus, b, c*MeV, a*MeV,
825  daughterFloatLevel);
826  theDecayTable->Insert(aSpontFissChannel);
827  modeSumBR[10] += b;
828  }
829  break;
830  case Triton:
831  {
832  G4TritonDecay* aTritonChannel =
833  new G4TritonDecay(&theParentNucleus, b, c*MeV, a*MeV,
834  daughterFloatLevel);
835  // anAlphaChannel->DumpNuclearInfo();
836  // anAlphaChannel->SetHLThreshold(halflifethreshold);
837  theDecayTable->Insert(aTritonChannel);
838  modeSumBR[11] += b;
839  }
840  break;
841 
842  case RDM_ERROR:
843 
844  default:
845  G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_000",
846  FatalException, "Selected decay mode does not exist");
847  } // switch
848  } // line < 72
849  } // if char == P
850  } // if char != #
851  } // While
852 
853  // Go through the decay table and make sure that the branching ratios are
854  // correctly normalised.
855 
856  G4VDecayChannel* theChannel = 0;
857  G4NuclearDecay* theNuclearDecayChannel = 0;
858  G4String mode = "";
859 
860  G4double theBR = 0.0;
861  for (G4int i = 0; i < theDecayTable->entries(); i++) {
862  theChannel = theDecayTable->GetDecayChannel(i);
863  theNuclearDecayChannel = static_cast<G4NuclearDecay*>(theChannel);
864  theDecayMode = theNuclearDecayChannel->GetDecayMode();
865 
866  if (theDecayMode != IT) {
867  theBR = theChannel->GetBR();
868  theChannel->SetBR(theBR*modeTotalBR[theDecayMode]/modeSumBR[theDecayMode]);
869  }
870  }
871  } // decay file exists
872 
873  DecaySchemeFile.close();
874 
875  if (!found && levelEnergy > 0) {
876  // Case where IT cascade for excited isotopes has no entries in RDM database
877  // Decay mode is isomeric transition.
878  G4ITDecay* anITChannel = new G4ITDecay(&theParentNucleus, 1.0, 0.0, 0.0,
880 // anITChannel->SetHLThreshold(halflifethreshold);
881  anITChannel->SetARM(applyARM);
882  theDecayTable->Insert(anITChannel);
883  }
884 
885  if (theDecayTable && GetVerboseLevel() > 1) {
886  theDecayTable->DumpInfo();
887  }
888 
889 #ifdef G4MULTITHREADED
890  //(*master_dkmap)[key] = theDecayTable; // store in master library
891 #endif
892  return theDecayTable;
893 }
894 
895 void
897 {
898  if (Z < 1 || A < 2) G4cout << "Z and A not valid!" << G4endl;
899 
900  std::ifstream DecaySchemeFile(filename);
901  if (DecaySchemeFile) {
902  G4int ID_ion = A*1000 + Z;
904  } else {
906  ed << filename << " does not exist! " << G4endl;
907  G4Exception("G4RadioactiveDecayBase::AddUserDecayDataFile()", "HAD_RDM_001",
908  FatalException, ed);
909  }
910 }
911 
913 // //
914 // DecayIt //
915 // //
917 
920 {
921  // Initialize G4ParticleChange object, get particle details and decay table
924  const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
925  const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
926 
927  // First check whether RDM applies to the current logical volume
928  if (!isAllVolumesMode) {
929  if (!std::binary_search(ValidVolumes.begin(), ValidVolumes.end(),
930  theTrack.GetVolume()->GetLogicalVolume()->GetName())) {
931 #ifdef G4VERBOSE
932  if (GetVerboseLevel()>1) {
933  G4cout <<"G4RadioactiveDecay::DecayIt : "
934  << theTrack.GetVolume()->GetLogicalVolume()->GetName()
935  << " is not selected for the RDM"<< G4endl;
936  G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
937  G4cout << " The Valid volumes are " << G4endl;
938  for (size_t i = 0; i< ValidVolumes.size(); i++)
939  G4cout << ValidVolumes[i] << G4endl;
940  }
941 #endif
943 
944  // Kill the parent particle.
949  }
950  }
951 
952  // Now check if particle is valid for RDM
953  if (!(IsApplicable(*theParticleDef) ) ) {
954  // Particle is not an ion or is outside the nucleuslimits for decay
955 #ifdef G4VERBOSE
956  if (GetVerboseLevel() > 1) {
957  G4cout << "G4RadioactiveDecay::DecayIt : "
958  << theParticleDef->GetParticleName()
959  << " is not an ion or is outside (Z,A) limits set for the decay. "
960  << " Set particle change accordingly. "
961  << G4endl;
962  }
963 #endif
965 
966  // Kill the parent particle
971  }
972 
973  G4DecayTable* theDecayTable = GetDecayTable(theParticleDef);
974 
975  if (theDecayTable == 0 || theDecayTable->entries() == 0) {
976  // No data in the decay table. Set particle change parameters
977  // to indicate this.
978 #ifdef G4VERBOSE
979  if (GetVerboseLevel() > 1) {
980  G4cout << "G4RadioactiveDecay::DecayIt : "
981  << "decay table not defined for "
982  << theParticleDef->GetParticleName()
983  << ". Set particle change accordingly. "
984  << G4endl;
985  }
986 #endif
988 
989  // Kill the parent particle.
994 
995  } else {
996  // Data found. Try to decay nucleus
997 
998 /*
999  G4double energyDeposit = 0.0;
1000  G4double finalGlobalTime = theTrack.GetGlobalTime();
1001  G4double finalLocalTime = theTrack.GetLocalTime();
1002  G4int index;
1003  G4ThreeVector currentPosition;
1004  currentPosition = theTrack.GetPosition();
1005 
1006  G4DecayProducts* products = DoDecay(*theParticleDef);
1007 
1008  // If the product is the same as the input kill the track if
1009  // necessary to prevent infinite loop (11/05/10, F.Lei)
1010  if (products->entries() == 1) {
1011  fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1012  fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill);
1013  fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1014  ClearNumberOfInteractionLengthLeft();
1015  return &fParticleChangeForRadDecay;
1016  }
1017 
1018  // Get parent particle information and boost the decay products to the
1019  // laboratory frame based on this information.
1020 
1021  // The Parent Energy used for the boost should be the total energy of
1022  // the nucleus of the parent ion without the energy of the shell electrons
1023  // (correction for bug 1359 by L. Desorgher)
1024  G4double ParentEnergy = theParticle->GetKineticEnergy()
1025  + theParticle->GetParticleDefinition()->GetPDGMass();
1026  G4ThreeVector ParentDirection(theParticle->GetMomentumDirection());
1027 
1028  if (theTrack.GetTrackStatus() == fStopButAlive) {
1029  // This condition seems to be always True, further investigation is needed
1030  // (L.Desorgher)
1031  // The particle is decayed at rest.
1032  // since the time is still for rest particle in G4 we need to add the
1033  // additional time lapsed between the particle come to rest and the
1034  // actual decay. This time is simply sampled with the mean-life of
1035  // the particle. But we need to protect the case PDGTime < 0.
1036  // (F.Lei 11/05/10)
1037  G4double temptime = -std::log( G4UniformRand())
1038  *theParticleDef->GetPDGLifeTime();
1039  if (temptime < 0.) temptime = 0.;
1040  finalGlobalTime += temptime;
1041  finalLocalTime += temptime;
1042  energyDeposit += theParticle->GetKineticEnergy();
1043  }
1044  products->Boost(ParentEnergy, ParentDirection);
1045 
1046  // Add products in theParticleChangeForRadDecay.
1047  G4int numberOfSecondaries = products->entries();
1048  fParticleChangeForRadDecay.SetNumberOfSecondaries(numberOfSecondaries);
1049 
1050 #ifdef G4VERBOSE
1051  if (GetVerboseLevel()>1) {
1052  G4cout <<"G4RadioactiveDecay::DecayIt : Decay vertex :";
1053  G4cout <<" Time: " <<finalGlobalTime/ns <<"[ns]";
1054  G4cout <<" X:" <<(theTrack.GetPosition()).x() /cm <<"[cm]";
1055  G4cout <<" Y:" <<(theTrack.GetPosition()).y() /cm <<"[cm]";
1056  G4cout <<" Z:" <<(theTrack.GetPosition()).z() /cm <<"[cm]";
1057  G4cout << G4endl;
1058  G4cout <<"G4Decay::DecayIt : decay products in Lab. Frame" <<G4endl;
1059  products->DumpInfo();
1060  products->IsChecked();
1061  }
1062 #endif
1063  for (index=0; index < numberOfSecondaries; index++) {
1064  G4Track* secondary = new G4Track(products->PopProducts(),
1065  finalGlobalTime, currentPosition);
1066  secondary->SetGoodForTrackingFlag();
1067  secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
1068  fParticleChangeForRadDecay.AddSecondary(secondary);
1069  }
1070  delete products;
1071 
1072  // Kill the parent particle
1073  fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
1074  fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(energyDeposit);
1075  fParticleChangeForRadDecay.ProposeLocalTime(finalLocalTime);
1076  // Reset NumberOfInteractionLengthLeft.
1077  ClearNumberOfInteractionLengthLeft();
1078 */
1079  // Decay without variance reduction
1080  DecayAnalog(theTrack);
1081  return &fParticleChangeForRadDecay ;
1082  }
1083 }
1084 
1085 
1087 {
1088  const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
1089  const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
1090  G4DecayProducts* products = DoDecay(*theParticleDef);
1091 
1092  // Check if the product is the same as input and kill the track if
1093  // necessary to prevent infinite loop (11/05/10, F.Lei)
1094  if (products->entries() == 1) {
1099  return;
1100  }
1101 
1102  G4double energyDeposit = 0.0;
1103  G4double finalGlobalTime = theTrack.GetGlobalTime();
1104  G4double finalLocalTime = theTrack.GetLocalTime();
1105 
1106  // Get parent particle information and boost the decay products to the
1107  // laboratory frame
1108 
1109  // ParentEnergy used for the boost should be the total energy of the nucleus
1110  // of the parent ion without the energy of the shell electrons
1111  // (correction for bug 1359 by L. Desorgher)
1112  G4double ParentEnergy = theParticle->GetKineticEnergy()
1113  + theParticle->GetParticleDefinition()->GetPDGMass();
1114  G4ThreeVector ParentDirection(theParticle->GetMomentumDirection());
1115 
1116  if (theTrack.GetTrackStatus() == fStopButAlive) {
1117  // this condition seems to be always True, further investigation is needed (L.Desorgher)
1118 
1119  // The particle is decayed at rest
1120  // Since the time is for the particle at rest, need to add additional time
1121  // lapsed between particle coming to rest and the actual decay. This time
1122  // is sampled with the mean-life of the particle. Need to protect the case
1123  // PDGTime < 0. (F.Lei 11/05/10)
1124  G4double temptime = -std::log(G4UniformRand() ) *
1125  theParticleDef->GetPDGLifeTime();
1126  if (temptime < 0.) temptime = 0.;
1127  finalGlobalTime += temptime;
1128  finalLocalTime += temptime;
1129  energyDeposit += theParticle->GetKineticEnergy();
1130  }
1131  products->Boost(ParentEnergy, ParentDirection);
1132 
1133  // Add products in theParticleChangeForRadDecay.
1134  G4int numberOfSecondaries = products->entries();
1136 
1137  if (GetVerboseLevel() > 1) {
1138  G4cout << "G4RadioactiveDecay::DecayAnalog: Decay vertex :";
1139  G4cout << " Time: " << finalGlobalTime/ns << "[ns]";
1140  G4cout << " X:" << (theTrack.GetPosition()).x() /cm << "[cm]";
1141  G4cout << " Y:" << (theTrack.GetPosition()).y() /cm << "[cm]";
1142  G4cout << " Z:" << (theTrack.GetPosition()).z() /cm << "[cm]";
1143  G4cout << G4endl;
1144  G4cout << "G4Decay::DecayIt : decay products in Lab. Frame" << G4endl;
1145  products->DumpInfo();
1146  products->IsChecked();
1147  }
1148 
1149  for (G4int index = 0; index < numberOfSecondaries; index++) {
1150  G4Track* secondary = new G4Track(products->PopProducts(), finalGlobalTime,
1151  theTrack.GetPosition() );
1153  //Change for atomics relaxation
1154  if (theRadDecayMode == IT && index>0){
1155  if (index == numberOfSecondaries-1) secondary->SetCreatorModelIndex(IT);
1156  else secondary->SetCreatorModelIndex(30);
1157  }
1159  && index <numberOfSecondaries-1){
1160  secondary->SetCreatorModelIndex(30);
1161  }
1162  secondary->SetGoodForTrackingFlag();
1163  secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
1165  }
1166 
1167  delete products;
1168 
1169  // Kill the parent particle
1173 
1174  // Reset NumberOfInteractionLengthLeft.
1176 }
1177 
1178 
1181 {
1182  G4DecayProducts* products = 0;
1183  G4DecayTable* theDecayTable = GetDecayTable(&theParticleDef);
1184 
1185  // Choose a decay channel.
1186  // G4DecayTable::SelectADecayChannel checks to see if sum of daughter masses
1187  // exceeds parent mass. Pass it the parent mass + maximum Q value to account
1188  // for difference in mass defect.
1189  G4double parentPlusQ = theParticleDef.GetPDGMass() + 30.*MeV;
1190  G4VDecayChannel* theDecayChannel = theDecayTable->SelectADecayChannel(parentPlusQ);
1191 
1192  if (theDecayChannel == 0) {
1193  // Decay channel not found.
1195  ed << " Cannot determine decay channel for " << theParticleDef.GetParticleName() << G4endl;
1196  G4Exception("G4RadioactiveDecay::DoDecay", "HAD_RDM_013",
1197  FatalException, ed);
1198  } else {
1199  // A decay channel has been identified, so execute the DecayIt.
1200 #ifdef G4VERBOSE
1201  if (GetVerboseLevel() > 1) {
1202  G4cout << "G4RadioactiveDecay::DoIt : selected decay channel addr: "
1203  << theDecayChannel << G4endl;
1204  }
1205 #endif
1206  theRadDecayMode = (static_cast<G4NuclearDecay*>(theDecayChannel))->GetDecayMode();
1207  products = theDecayChannel->DecayIt(theParticleDef.GetPDGMass() );
1208 
1209  // Apply directional bias if requested by user
1210  CollimateDecay(products);
1211  }
1212 
1213  return products;
1214 }
1215 
1216 
1217 // Apply directional bias for "visible" daughters (e+-, gamma, n, p, alpha)
1218 
1220  if (origin == forceDecayDirection) return; // No collimation requested
1221  if (180.*deg == forceDecayHalfAngle) return;
1222  if (0 == products || 0 == products->entries()) return;
1223 
1224 #ifdef G4VERBOSE
1225  if (GetVerboseLevel() > 1) G4cout << "Begin of CollimateDecay..." << G4endl;
1226 #endif
1227 
1228  // Particles suitable for directional biasing (for if-blocks below)
1232  static const G4ParticleDefinition* gamma = G4Gamma::Definition();
1236 
1237  G4ThreeVector newDirection; // Re-use to avoid memory churn
1238  for (G4int i=0; i<products->entries(); i++) {
1239  G4DynamicParticle* daughter = (*products)[i];
1240  const G4ParticleDefinition* daughterType =
1241  daughter->GetParticleDefinition();
1242  if (daughterType == electron || daughterType == positron ||
1243  daughterType == neutron || daughterType == gamma ||
1244  daughterType == alpha || daughterType == triton || daughterType == proton) CollimateDecayProduct(daughter);
1245  }
1246 }
1247 
1249 #ifdef G4VERBOSE
1250  if (GetVerboseLevel() > 1) {
1251  G4cout << "CollimateDecayProduct for daughter "
1252  << daughter->GetParticleDefinition()->GetParticleName() << G4endl;
1253  }
1254 #endif
1255 
1257  if (origin != collimate) daughter->SetMomentumDirection(collimate);
1258 }
1259 
1260 
1261 // Choose random direction within collimation cone
1262 
1264  if (origin == forceDecayDirection) return origin; // Don't do collimation
1265  if (forceDecayHalfAngle == 180.*deg) return origin;
1266 
1268 
1269  // Return direction offset by random throw
1270  if (forceDecayHalfAngle > 0.) {
1271  // Generate uniform direction around central axis
1272  G4double phi = 2.*pi*G4UniformRand();
1273  G4double cosMin = std::cos(forceDecayHalfAngle);
1274  G4double cosTheta = (1.-cosMin)*G4UniformRand() + cosMin; // [cosMin,1.)
1275 
1276  dir.setPhi(dir.phi()+phi);
1277  dir.setTheta(dir.theta()+std::acos(cosTheta));
1278  }
1279 
1280 #ifdef G4VERBOSE
1281  if (GetVerboseLevel()>1)
1282  G4cout << " ChooseCollimationDirection returns " << dir << G4endl;
1283 #endif
1284 
1285  return dir;
1286 }
1287