Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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.hh //
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 "G4ProtonDecay.hh"
54 #include "G4NeutronDecay.hh"
55 #include "G4VDecayChannel.hh"
56 #include "G4NuclearDecay.hh"
58 #include "G4Fragment.hh"
59 #include "G4Ions.hh"
60 #include "G4IonTable.hh"
61 #include "G4BetaDecayType.hh"
62 #include "Randomize.hh"
63 #include "G4LogicalVolumeStore.hh"
64 #include "G4NuclearLevelData.hh"
65 #include "G4DeexPrecoParameters.hh"
66 #include "G4LevelManager.hh"
67 #include "G4ThreeVector.hh"
68 #include "G4Electron.hh"
69 #include "G4Positron.hh"
70 #include "G4Neutron.hh"
71 #include "G4Gamma.hh"
72 #include "G4Alpha.hh"
73 #include "G4Proton.hh"
74 
75 #include "G4HadronicProcessType.hh"
77 #include "G4HadronicException.hh"
78 #include "G4LossTableManager.hh"
79 #include "G4VAtomDeexcitation.hh"
80 #include "G4UAtomicDeexcitation.hh"
81 #include "G4PhotonEvaporation.hh"
82 
83 #include <vector>
84 #include <sstream>
85 #include <algorithm>
86 #include <fstream>
87 
88 using namespace CLHEP;
89 
92 
93 #ifdef G4MULTITHREADED
94 #include "G4AutoLock.hh"
95 G4Mutex G4RadioactiveDecayBase::radioactiveDecayMutex = G4MUTEX_INITIALIZER;
96 DecayTableMap* G4RadioactiveDecayBase::master_dkmap = 0;
97 #endif
98 
100  : G4VRestDiscreteProcess(processName, fDecay), isInitialised(false),
101  forceDecayDirection(0.,0.,0.), forceDecayHalfAngle(0.*deg), dirPath(""),
102  verboseLevel(0)
103 {
104 #ifdef G4VERBOSE
105  if (GetVerboseLevel() > 1) {
106  G4cout << "G4RadioactiveDecayBase constructor: processName = " << processName
107  << G4endl;
108  }
109 #endif
110 
112 
115 
116  // Set up photon evaporation for use in G4ITDecay
119  photonEvaporation->SetICM(true);
120 
122  deex->SetCorrelatedGamma(true);
123 
124  // Check data directory
125  char* path_var = getenv("G4RADIOACTIVEDATA");
126  if (!path_var) {
127  G4Exception("G4RadioactiveDecay()","HAD_RDM_200",FatalException,
128  "Environment variable G4RADIOACTIVEDATA is not set");
129  } else {
130  dirPath = path_var; // convert to string
131  std::ostringstream os;
132  os << dirPath << "/z1.a3"; // used as a dummy
133  std::ifstream testFile;
134  testFile.open(os.str() );
135  if (!testFile.is_open() )
136  G4Exception("G4RadioactiveDecay()","HAD_RDM_201",FatalException,
137  "Environment variable G4RADIOACTIVEDATA is set, but does not point to correct directory");
138  }
139 
140  // Reset the list of user defined data files
142 
143  // Instantiate the map of decay tables
144 #ifdef G4MULTITHREADED
145  G4AutoLock lk(&G4RadioactiveDecayBase::radioactiveDecayMutex);
146  if(!master_dkmap) master_dkmap = new DecayTableMap;
147 #endif
148  dkmap = new DecayTableMap;
149 
150  // Apply default values
151  applyARM = true;
152  applyICM = true; // Always on; keep only for backward compatibility
153 
154  // RDM applies to all logical volumes by default
155  isAllVolumesMode = true;
158 }
159 
160 void G4RadioactiveDecayBase::ProcessDescription(std::ostream& outFile) const
161 {
162  outFile << "The radioactive decay process (G4RadioactiveDecay) handles the\n"
163  << "alpha, beta+, beta-, electron capture and isomeric transition\n"
164  << "decays of nuclei (G4GenericIon) with masses A > 4.\n"
165  << "The required half-lives and decay schemes are retrieved from\n"
166  << "the RadioactiveDecay database which was derived from ENSDF.\n";
167 }
168 
169 
171 {
173  delete photonEvaporation;
174  for (DecayTableMap::iterator i = dkmap->begin(); i != dkmap->end(); i++) {
175  delete i->second;
176  }
177  dkmap->clear();
178  delete dkmap;
179 }
180 
181 
183 {
184  // All particles other than G4Ions, are rejected by default
185  if (((const G4Ions*)(&aParticle))->GetExcitationEnergy() > 0.) {return true;}
186  if (aParticle.GetParticleName() == "GenericIon") {
187  return true;
188  } else if (!(aParticle.GetParticleType() == "nucleus")
189  || aParticle.GetPDGLifeTime() < 0. ) {
190  return false;
191  }
192 
193  // Determine whether the nuclide falls into the correct A and Z range
194  G4int A = ((const G4Ions*) (&aParticle))->GetAtomicMass();
195  G4int Z = ((const G4Ions*) (&aParticle))->GetAtomicNumber();
196 
198  {return false;}
199  else if (Z > theNucleusLimits.GetZMax() || Z < theNucleusLimits.GetZMin())
200  {return false;}
201  return true;
202 }
203 
205 {
206  G4String key = aNucleus->GetParticleName();
207  DecayTableMap::iterator table_ptr = dkmap->find(key);
208 
209  G4DecayTable* theDecayTable = 0;
210  if (table_ptr == dkmap->end() ) { // If table not there,
211  theDecayTable = LoadDecayTable(*aNucleus); // load from file and
212  if(theDecayTable) (*dkmap)[key] = theDecayTable; // store in library
213  } else {
214  theDecayTable = table_ptr->second;
215  }
216 
217  return theDecayTable;
218 }
219 
220 
222 {
223  G4LogicalVolumeStore* theLogicalVolumes;
224  G4LogicalVolume* volume;
225  theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
226  for (size_t i = 0; i < theLogicalVolumes->size(); i++) {
227  volume=(*theLogicalVolumes)[i];
228  if (volume->GetName() == aVolume) {
229  ValidVolumes.push_back(aVolume);
230  std::sort(ValidVolumes.begin(), ValidVolumes.end());
231  // sort need for performing binary_search
232 #ifdef G4VERBOSE
233  if (GetVerboseLevel()>0)
234  G4cout << " RDM Applies to : " << aVolume << G4endl;
235 #endif
236  } else if(i == theLogicalVolumes->size()) {
237  G4cerr << "SelectAVolume: "<< aVolume
238  << " is not a valid logical volume name" << G4endl;
239  }
240  }
241 }
242 
243 
245 {
246  G4LogicalVolumeStore* theLogicalVolumes;
247  G4LogicalVolume* volume;
248  theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
249  for (size_t i = 0; i < theLogicalVolumes->size(); i++){
250  volume=(*theLogicalVolumes)[i];
251  if (volume->GetName() == aVolume) {
252  std::vector<G4String>::iterator location;
253  location = std::find(ValidVolumes.begin(),ValidVolumes.end(),aVolume);
254  if (location != ValidVolumes.end()) {
255  ValidVolumes.erase(location);
256  std::sort(ValidVolumes.begin(), ValidVolumes.end());
257  isAllVolumesMode =false;
258  } else {
259  G4cerr << " DeselectVolume:" << aVolume << " is not in the list "
260  << G4endl;
261  }
262 #ifdef G4VERBOSE
263  if (GetVerboseLevel() > 0)
264  G4cout << " DeselectVolume: " << aVolume << " is removed from list "
265  << G4endl;
266 #endif
267  } else if (i == theLogicalVolumes->size()) {
268  G4cerr << " DeselectVolume:" << aVolume
269  << "is not a valid logical volume name" << G4endl;
270  }
271  }
272 }
273 
274 
276 {
277  G4LogicalVolumeStore* theLogicalVolumes;
278  G4LogicalVolume* volume;
279  theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
280  ValidVolumes.clear();
281 #ifdef G4VERBOSE
282  if (GetVerboseLevel()>0)
283  G4cout << " RDM Applies to all Volumes" << G4endl;
284 #endif
285  for (size_t i = 0; i < theLogicalVolumes->size(); i++){
286  volume = (*theLogicalVolumes)[i];
287  ValidVolumes.push_back(volume->GetName());
288 #ifdef G4VERBOSE
289  if (GetVerboseLevel()>0)
290  G4cout << " RDM Applies to Volume " << volume->GetName() << G4endl;
291 #endif
292  }
293  std::sort(ValidVolumes.begin(), ValidVolumes.end());
294  // sort needed in order to allow binary_search
295  isAllVolumesMode=true;
296 }
297 
298 
300 {
301  ValidVolumes.clear();
302  isAllVolumesMode=false;
303 #ifdef G4VERBOSE
304  if (GetVerboseLevel() > 0) G4cout << "RDM removed from all volumes" << G4endl;
305 #endif
306 }
307 
308 
310 // //
311 // GetMeanLifeTime (required by the base class) //
312 // //
314 
317 {
318  // For variance reduction the time is set to 0 so as to force the particle
319  // to decay immediately.
320  // In analogueMC mode it returns the particle's mean-life.
321 
322  G4double meanlife = 0.;
323 // if (AnalogueMC) {
324  const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
325  const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
326  G4double theLife = theParticleDef->GetPDGLifeTime();
327 #ifdef G4VERBOSE
328  if (GetVerboseLevel() > 2) {
329  G4cout << "G4RadioactiveDecay::GetMeanLifeTime() " << G4endl;
330  G4cout << "KineticEnergy: " << theParticle->GetKineticEnergy()/GeV
331  << " GeV, Mass: " << theParticle->GetMass()/GeV
332  << " GeV, Life time: " << theLife/ns << " ns " << G4endl;
333  }
334 #endif
335  if (theParticleDef->GetPDGStable()) {meanlife = DBL_MAX;}
336  else if (theLife < 0.0) {meanlife = DBL_MAX;}
337  else {meanlife = theLife;}
338  // Set meanlife to zero for excited istopes which are not in the
339  // RDM database
340  if (((const G4Ions*)(theParticleDef))->GetExcitationEnergy() > 0. &&
341  meanlife == DBL_MAX) {meanlife = 0.;}
342 // }
343 #ifdef G4VERBOSE
344  if (GetVerboseLevel() > 1)
345  G4cout << " mean life time: " << meanlife/s << " s " << G4endl;
346 #endif
347 
348  return meanlife;
349 }
350 
352 // //
353 // GetMeanFreePath for decay in flight //
354 // //
356 
359 {
360  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
361  const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
362  G4double tau = aParticleDef->GetPDGLifeTime();
363  G4double aMass = aParticle->GetMass();
364 
365 #ifdef G4VERBOSE
366  if (GetVerboseLevel() > 2) {
367  G4cout << "G4RadioactiveDecay::GetMeanFreePath() " << G4endl;
368  G4cout << " KineticEnergy: " << aParticle->GetKineticEnergy()/GeV
369  << " GeV, Mass: " << aMass/GeV << " GeV, tau: " << tau << " ns "
370  << G4endl;
371  }
372 #endif
373  G4double pathlength = DBL_MAX;
374  if (tau != -1) {
375  // Ion can decay
376 
377  if (tau < -1000.0) {
378  pathlength = DBL_MIN; // nuclide had very short lifetime or wasn't in table
379 
380  } else if (tau < 0.0) {
381  G4cout << aParticleDef->GetParticleName() << " has lifetime " << tau << G4endl;
383  ed << "Ion has negative lifetime " << tau
384  << " but is not stable. Setting mean free path to DBL_MAX" << G4endl;
385  G4Exception("G4RadioactiveDecay::GetMeanFreePath()", "HAD_RDM_011",
386  JustWarning, ed);
387  pathlength = DBL_MAX;
388 
389  } else {
390  // Calculate mean free path
391  G4double betaGamma = aParticle->GetTotalMomentum()/aMass;
392  pathlength = c_light*tau*betaGamma;
393 
394  if (pathlength < DBL_MIN) {
395  pathlength = DBL_MIN;
396 #ifdef G4VERBOSE
397  if (GetVerboseLevel() > 2) {
398  G4cout << "G4Decay::GetMeanFreePath: "
399  << aParticleDef->GetParticleName()
400  << " stops, kinetic energy = "
401  << aParticle->GetKineticEnergy()/keV <<" keV " << G4endl;
402  }
403 #endif
404  }
405  }
406  }
407 
408 #ifdef G4VERBOSE
409  if (GetVerboseLevel() > 1) {
410  G4cout << "mean free path: "<< pathlength/m << " m" << G4endl;
411  }
412 #endif
413  return pathlength;
414 }
415 
417 // //
418 // BuildPhysicsTable - initialization of atomic de-excitation //
419 // //
421 
423 {
424  if (!isInitialised) {
425  isInitialised = true;
426  }
427 
430 }
431 
433 // //
434 // LoadDecayTable loads the decay scheme from the RadioactiveDecay database //
435 // for the parent nucleus. //
436 // //
438 
441 {
442  // Generate input data file name using Z and A of the parent nucleus
443  // file containing radioactive decay data.
444  G4int A = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
445  G4int Z = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
446 
447  G4double levelEnergy = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
448  G4Ions::G4FloatLevelBase floatingLevel =
449  ((const G4Ions*)(&theParentNucleus))->GetFloatLevelBase();
450 
451 #ifdef G4MULTITHREADED
452  G4AutoLock lk(&G4RadioactiveDecayBase::radioactiveDecayMutex);
453 
454  G4String key = theParentNucleus.GetParticleName();
455  DecayTableMap::iterator master_table_ptr = master_dkmap->find(key);
456 
457  if (master_table_ptr != master_dkmap->end() ) { // If table is there
458  return master_table_ptr->second;
459  }
460 #endif
461 
462  //Check if data have been provided by the user
464 
465  if (file == "") {
466 /*
467  if (!getenv("G4RADIOACTIVEDATA") ) {
468  G4cout << "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files."
469  << G4endl;
470  throw G4HadronicException(__FILE__, __LINE__, " Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files.");
471  }
472  G4String dirName = getenv("G4RADIOACTIVEDATA");
473 */
474  std::ostringstream os;
475  os << dirPath << "/z" << Z << ".a" << A << '\0';
476  file = os.str();
477  }
478 
479  G4DecayTable* theDecayTable = new G4DecayTable();
480  G4bool found(false); // True if energy level matches one in table
481 
482  std::ifstream DecaySchemeFile;
483  DecaySchemeFile.open(file);
484 
485  if (DecaySchemeFile.good()) {
486  // Initialize variables used for reading in radioactive decay data
487  G4bool floatMatch(false);
488  const G4int nMode = 9;
489  G4double modeTotalBR[nMode] = {0.0};
490  G4double modeSumBR[nMode];
491  for (G4int i = 0; i < nMode; i++) {
492  modeSumBR[i] = 0.0;
493  }
494 
495  char inputChars[120]={' '};
496  G4String inputLine;
497  G4String recordType("");
498  G4String floatingFlag("");
499  G4String daughterFloatFlag("");
500  G4Ions::G4FloatLevelBase daughterFloatLevel;
501  G4RadioactiveDecayMode theDecayMode;
502  G4double decayModeTotal(0.0);
503  G4double parentExcitation(0.0);
504  G4double a(0.0);
505  G4double b(0.0);
506  G4double c(0.0);
507  G4double dummy(0.0);
508  G4BetaDecayType betaType(allowed);
509 
510  // Loop through each data file record until you identify the decay
511  // data relating to the nuclide of concern.
512 
513  G4bool complete(false); // bool insures only one set of values read for any
514  // given parent energy level
515  G4int loop = 0;
517  ed << " While count exceeded " << G4endl;
518 
519  while (!complete && !DecaySchemeFile.getline(inputChars, 120).eof()) { /* Loop checking, 01.09.2015, D.Wright */
520  loop++;
521  if (loop > 100000) {
522  G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_100", JustWarning, ed);
523  break;
524  }
525 
526  inputLine = inputChars;
527  inputLine = inputLine.strip(1);
528  if (inputChars[0] != '#' && inputLine.length() != 0) {
529  std::istringstream tmpStream(inputLine);
530 
531  if (inputChars[0] == 'P') {
532  // Nucleus is a parent type. Check excitation level to see if it
533  // matches that of theParentNucleus
534  tmpStream >> recordType >> parentExcitation >> floatingFlag >> dummy;
535  // "dummy" takes the place of half-life
536  // Now read in from ENSDFSTATE in particle category
537 
538  if (found) {
539  complete = true;
540  } else {
541  // Take first level which matches excitation energy regardless of floating level
542  found = (std::abs(parentExcitation*keV - levelEnergy) < levelTolerance);
543  if (floatingLevel != noFloat) {
544  // If floating level specificed, require match of both energy and floating level
545  floatMatch = (floatingLevel == G4Ions::FloatLevelBase(floatingFlag.back()) );
546  if (!floatMatch) found = false;
547  }
548  }
549 
550  } else if (found) {
551  // The right part of the radioactive decay data file has been found. Search
552  // through it to determine the mode of decay of the subsequent records.
553 
554  // Store for later the total decay probability for each decay mode
555  if (inputLine.length() < 72) {
556  tmpStream >> theDecayMode >> dummy >> decayModeTotal;
557  switch (theDecayMode) {
558  case IT:
559  {
560  G4ITDecay* anITChannel = new G4ITDecay(&theParentNucleus, decayModeTotal,
561  0.0, 0.0, photonEvaporation);
562 // anITChannel->SetHLThreshold(halflifethreshold);
563  anITChannel->SetARM(applyARM);
564  theDecayTable->Insert(anITChannel);
565 // anITChannel->DumpNuclearInfo();
566  }
567  break;
568  case BetaMinus:
569  modeTotalBR[1] = decayModeTotal; break;
570  case BetaPlus:
571  modeTotalBR[2] = decayModeTotal; break;
572  case KshellEC:
573  modeTotalBR[3] = decayModeTotal; break;
574  case LshellEC:
575  modeTotalBR[4] = decayModeTotal; break;
576  case MshellEC:
577  modeTotalBR[5] = decayModeTotal; break;
578  case Alpha:
579  modeTotalBR[6] = decayModeTotal; break;
580  case Proton:
581  modeTotalBR[7] = decayModeTotal; break;
582  case Neutron:
583  modeTotalBR[8] = decayModeTotal; break;
584  case BDProton:
585  break;
586  case BDNeutron:
587  break;
588  case Beta2Minus:
589  break;
590  case Beta2Plus:
591  break;
592  case Proton2:
593  break;
594  case Neutron2:
595  break;
596  case SpFission:
597  break;
598  case RDM_ERROR:
599 
600  default:
601  G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_000",
602  FatalException, "Selected decay mode does not exist");
603  } // switch
604 
605  } else {
606  if (inputLine.length() < 84) {
607  tmpStream >> theDecayMode >> a >> daughterFloatFlag >> b >> c;
608  betaType = allowed;
609  } else {
610  tmpStream >> theDecayMode >> a >> daughterFloatFlag >> b >> c >> betaType;
611  }
612 
613  // Allowed transitions are the default. Forbidden transitions are
614  // indicated in the last column.
615  a /= 1000.;
616  c /= 1000.;
617  daughterFloatLevel = G4Ions::FloatLevelBase(daughterFloatFlag.back());
618 
619  switch (theDecayMode) {
620  case BetaMinus:
621  {
622  G4BetaMinusDecay* aBetaMinusChannel =
623  new G4BetaMinusDecay(&theParentNucleus, b, c*MeV, a*MeV,
624  daughterFloatLevel, betaType);
625 // aBetaMinusChannel->DumpNuclearInfo();
626 // aBetaMinusChannel->SetHLThreshold(halflifethreshold);
627  theDecayTable->Insert(aBetaMinusChannel);
628  modeSumBR[1] += b;
629  }
630  break;
631 
632  case BetaPlus:
633  {
634  G4BetaPlusDecay* aBetaPlusChannel =
635  new G4BetaPlusDecay(&theParentNucleus, b, c*MeV, a*MeV,
636  daughterFloatLevel, betaType);
637 // aBetaPlusChannel->DumpNuclearInfo();
638 // aBetaPlusChannel->SetHLThreshold(halflifethreshold);
639  theDecayTable->Insert(aBetaPlusChannel);
640  modeSumBR[2] += b;
641  }
642  break;
643 
644  case KshellEC: // K-shell electron capture
645  {
646  G4ECDecay* aKECChannel =
647  new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
648  daughterFloatLevel, KshellEC);
649 // aKECChannel->DumpNuclearInfo();
650 // aKECChannel->SetHLThreshold(halflifethreshold);
651  aKECChannel->SetARM(applyARM);
652  theDecayTable->Insert(aKECChannel);
653  modeSumBR[3] += b;
654  }
655  break;
656 
657  case LshellEC: // L-shell electron capture
658  {
659  G4ECDecay* aLECChannel =
660  new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
661  daughterFloatLevel, LshellEC);
662 // aLECChannel->DumpNuclearInfo();
663 // aLECChannel->SetHLThreshold(halflifethreshold);
664  aLECChannel->SetARM(applyARM);
665  theDecayTable->Insert(aLECChannel);
666  modeSumBR[4] += b;
667  }
668  break;
669 
670  case MshellEC: // M-shell electron capture
671  {
672  G4ECDecay* aMECChannel =
673  new G4ECDecay(&theParentNucleus, b, c*MeV, a*MeV,
674  daughterFloatLevel, MshellEC);
675 // aMECChannel->DumpNuclearInfo();
676 // aMECChannel->SetHLThreshold(halflifethreshold);
677  aMECChannel->SetARM(applyARM);
678  theDecayTable->Insert(aMECChannel);
679  modeSumBR[5] += b;
680  }
681  break;
682 
683  case Alpha:
684  {
685  G4AlphaDecay* anAlphaChannel =
686  new G4AlphaDecay(&theParentNucleus, b, c*MeV, a*MeV,
687  daughterFloatLevel);
688 // anAlphaChannel->DumpNuclearInfo();
689 // anAlphaChannel->SetHLThreshold(halflifethreshold);
690  theDecayTable->Insert(anAlphaChannel);
691  modeSumBR[6] += b;
692  }
693  break;
694 
695  case Proton:
696  {
697  G4ProtonDecay* aProtonChannel =
698  new G4ProtonDecay(&theParentNucleus, b, c*MeV, a*MeV,
699  daughterFloatLevel);
700 // aProtonChannel->DumpNuclearInfo();
701 // aProtonChannel->SetHLThreshold(halflifethreshold);
702  theDecayTable->Insert(aProtonChannel);
703  modeSumBR[7] += b;
704  }
705  break;
706 
707  case Neutron:
708  {
709  G4NeutronDecay* aNeutronChannel =
710  new G4NeutronDecay(&theParentNucleus, b, c*MeV, a*MeV,
711  daughterFloatLevel);
712 // aNeutronChannel->DumpNuclearInfo();
713 // aNeutronChannel->SetHLThreshold(halflifethreshold);
714  theDecayTable->Insert(aNeutronChannel);
715  modeSumBR[8] += b;
716  }
717  break;
718 
719  case BDProton:
720  // Not yet implemented
721  // G4cout << " beta-delayed proton decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
722  break;
723  case BDNeutron:
724  // Not yet implemented
725  // G4cout << " beta-delayed neutron decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
726  break;
727  case Beta2Minus:
728  // Not yet implemented
729  // G4cout << " Double beta- decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
730  break;
731  case Beta2Plus:
732  // Not yet implemented
733  // G4cout << " Double beta+ decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
734  break;
735  case Proton2:
736  // Not yet implemented
737  // G4cout << " Double proton decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
738  break;
739  case Neutron2:
740  // Not yet implemented
741  // G4cout << " Double beta- decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
742  break;
743  case SpFission:
744  // Not yet implemented
745  //G4cout<<"Sp fission channel"<<a<<'\t'<<b<<'\t'<<c<<std::endl;
746  break;
747  case RDM_ERROR:
748 
749  default:
750  G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_000",
751  FatalException, "Selected decay mode does not exist");
752  } // switch
753  } // line < 72
754  } // if char == P
755  } // if char != #
756  } // While
757 
758  // Go through the decay table and make sure that the branching ratios are
759  // correctly normalised.
760 
761  G4VDecayChannel* theChannel = 0;
762  G4NuclearDecay* theNuclearDecayChannel = 0;
763  G4String mode = "";
764 
765  G4double theBR = 0.0;
766  for (G4int i = 0; i < theDecayTable->entries(); i++) {
767  theChannel = theDecayTable->GetDecayChannel(i);
768  theNuclearDecayChannel = static_cast<G4NuclearDecay*>(theChannel);
769  theDecayMode = theNuclearDecayChannel->GetDecayMode();
770 
771  if (theDecayMode != IT) {
772  theBR = theChannel->GetBR();
773  theChannel->SetBR(theBR*modeTotalBR[theDecayMode]/modeSumBR[theDecayMode]);
774  }
775  }
776  } // decay file exists
777 
778  DecaySchemeFile.close();
779 
780  if (!found && levelEnergy > 0) {
781  // Case where IT cascade for excited isotopes has no entries in RDM database
782  // Decay mode is isomeric transition.
783  G4ITDecay* anITChannel = new G4ITDecay(&theParentNucleus, 1.0, 0.0, 0.0,
785 // anITChannel->SetHLThreshold(halflifethreshold);
786  anITChannel->SetARM(applyARM);
787  theDecayTable->Insert(anITChannel);
788  }
789 
790  if (theDecayTable && GetVerboseLevel() > 1) {
791  theDecayTable->DumpInfo();
792  }
793 
794 #ifdef G4MULTITHREADED
795  //(*master_dkmap)[key] = theDecayTable; // store in master library
796 #endif
797  return theDecayTable;
798 }
799 
800 void
802 {
803  if (Z < 1 || A < 2) G4cout << "Z and A not valid!" << G4endl;
804 
805  std::ifstream DecaySchemeFile(filename);
806  if (DecaySchemeFile) {
807  G4int ID_ion = A*1000 + Z;
808  theUserRadioactiveDataFiles[ID_ion] = filename;
809  } else {
810  G4cout << "The file " << filename << " does not exist!" << G4endl;
811  }
812 }
813 
815 // //
816 // DecayIt //
817 // //
819 
822 {
823  // Initialize G4ParticleChange object, get particle details and decay table
824 
827  const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
828  const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
829 
830  // First check whether RDM applies to the current logical volume
831  if (!isAllVolumesMode) {
832  if (!std::binary_search(ValidVolumes.begin(), ValidVolumes.end(),
833  theTrack.GetVolume()->GetLogicalVolume()->GetName())) {
834 #ifdef G4VERBOSE
835  if (GetVerboseLevel()>0) {
836  G4cout <<"G4RadioactiveDecay::DecayIt : "
837  << theTrack.GetVolume()->GetLogicalVolume()->GetName()
838  << " is not selected for the RDM"<< G4endl;
839  G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
840  G4cout << " The Valid volumes are " << G4endl;
841  for (size_t i = 0; i< ValidVolumes.size(); i++)
842  G4cout << ValidVolumes[i] << G4endl;
843  }
844 #endif
846 
847  // Kill the parent particle.
852  }
853  }
854 
855  // Now check if particle is valid for RDM
856  if (!(IsApplicable(*theParticleDef) ) ) {
857  // Particle is not an ion or is outside the nucleuslimits for decay
858 
859 #ifdef G4VERBOSE
860  if (GetVerboseLevel()>0) {
861  G4cerr << "G4RadioactiveDecay::DecayIt : "
862  << theParticleDef->GetParticleName()
863  << " is not a valid nucleus for the RDM"<< G4endl;
864  }
865 #endif
867 
868  // Kill the parent particle
873  }
874  G4DecayTable* theDecayTable = GetDecayTable(theParticleDef);
875 
876  if (theDecayTable == 0 || theDecayTable->entries() == 0) {
877  // No data in the decay table. Set particle change parameters
878  // to indicate this.
879 #ifdef G4VERBOSE
880  if (GetVerboseLevel() > 0) {
881  G4cerr <<"G4RadioactiveDecay::DecayIt : decay table not defined for ";
882  G4cerr <<theParticleDef->GetParticleName() <<G4endl;
883  }
884 #endif
886 
887  // Kill the parent particle.
892 
893  } else {
894  // Data found. Try to decay nucleus
895  G4double energyDeposit = 0.0;
896  G4double finalGlobalTime = theTrack.GetGlobalTime();
897  G4double finalLocalTime = theTrack.GetLocalTime();
898  G4int index;
899  G4ThreeVector currentPosition;
900  currentPosition = theTrack.GetPosition();
901 
902  G4DecayProducts* products = DoDecay(*theParticleDef);
903 
904  // If the product is the same as the input kill the track if
905  // necessary to prevent infinite loop (11/05/10, F.Lei)
906  if (products->entries() == 1) {
912  }
913 
914  // Get parent particle information and boost the decay products to the
915  // laboratory frame based on this information.
916 
917  // The Parent Energy used for the boost should be the total energy of
918  // the nucleus of the parent ion without the energy of the shell electrons
919  // (correction for bug 1359 by L. Desorgher)
920  G4double ParentEnergy = theParticle->GetKineticEnergy()
921  + theParticle->GetParticleDefinition()->GetPDGMass();
922  G4ThreeVector ParentDirection(theParticle->GetMomentumDirection());
923 
924  if (theTrack.GetTrackStatus() == fStopButAlive) {
925  // This condition seems to be always True, further investigation is needed
926  // (L.Desorgher)
927  // The particle is decayed at rest.
928  // since the time is still for rest particle in G4 we need to add the
929  // additional time lapsed between the particle come to rest and the
930  // actual decay. This time is simply sampled with the mean-life of
931  // the particle. But we need to protect the case PDGTime < 0.
932  // (F.Lei 11/05/10)
933  G4double temptime = -std::log( G4UniformRand())
934  *theParticleDef->GetPDGLifeTime();
935  if (temptime < 0.) temptime = 0.;
936  finalGlobalTime += temptime;
937  finalLocalTime += temptime;
938  energyDeposit += theParticle->GetKineticEnergy();
939  }
940  products->Boost(ParentEnergy, ParentDirection);
941 
942  // Add products in theParticleChangeForRadDecay.
943  G4int numberOfSecondaries = products->entries();
945 
946 #ifdef G4VERBOSE
947  if (GetVerboseLevel()>1) {
948  G4cout <<"G4RadioactiveDecay::DecayIt : Decay vertex :";
949  G4cout <<" Time: " <<finalGlobalTime/ns <<"[ns]";
950  G4cout <<" X:" <<(theTrack.GetPosition()).x() /cm <<"[cm]";
951  G4cout <<" Y:" <<(theTrack.GetPosition()).y() /cm <<"[cm]";
952  G4cout <<" Z:" <<(theTrack.GetPosition()).z() /cm <<"[cm]";
953  G4cout << G4endl;
954  G4cout <<"G4Decay::DecayIt : decay products in Lab. Frame" <<G4endl;
955  products->DumpInfo();
956  products->IsChecked();
957  }
958 #endif
959  for (index=0; index < numberOfSecondaries; index++) {
960  G4Track* secondary = new G4Track(products->PopProducts(),
961  finalGlobalTime, currentPosition);
962  secondary->SetGoodForTrackingFlag();
963  secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
965  }
966  delete products;
967 
968  // Kill the parent particle
972  // Reset NumberOfInteractionLengthLeft.
974 
976  }
977 }
978 
979 
982 {
983  G4DecayProducts* products = 0;
984  G4DecayTable* theDecayTable = GetDecayTable(&theParticleDef);
985 
986  // Choose a decay channel.
987 #ifdef G4VERBOSE
988  if (GetVerboseLevel() > 0) G4cout << "Select a channel..." << G4endl;
989 #endif
990 
991  // G4DecayTable::SelectADecayChannel checks to see if sum of daughter masses
992  // exceeds parent mass. Pass it the parent mass + maximum Q value to account
993  // for difference in mass defect.
994  G4double parentPlusQ = theParticleDef.GetPDGMass() + 30.*MeV;
995  G4VDecayChannel* theDecayChannel = theDecayTable->SelectADecayChannel(parentPlusQ);
996 
997  if (theDecayChannel == 0) {
998  // Decay channel not found.
1000  ed << " Cannot determine decay channel for " << theParticleDef.GetParticleName() << G4endl;
1001  G4Exception("G4RadioactiveDecay::DoDecay", "HAD_RDM_013",
1002  FatalException, ed);
1003  } else {
1004  // A decay channel has been identified, so execute the DecayIt.
1005 #ifdef G4VERBOSE
1006  if (GetVerboseLevel() > 1) {
1007  G4cerr << "G4RadioactiveDecay::DoIt : selected decay channel addr:";
1008  G4cerr << theDecayChannel << G4endl;
1009  }
1010 #endif
1011  products = theDecayChannel->DecayIt(theParticleDef.GetPDGMass() );
1012 
1013  // Apply directional bias if requested by user
1014  CollimateDecay(products);
1015  }
1016 
1017  return products;
1018 }
1019 
1020 
1021 // Apply directional bias for "visible" daughters (e+-, gamma, n, p, alpha)
1022 
1024  if (origin == forceDecayDirection) return; // No collimation requested
1025  if (180.*deg == forceDecayHalfAngle) return;
1026  if (0 == products || 0 == products->entries()) return;
1027 
1028 #ifdef G4VERBOSE
1029  if (GetVerboseLevel() > 0) G4cout << "Begin of CollimateDecay..." << G4endl;
1030 #endif
1031 
1032  // Particles suitable for directional biasing (for if-blocks below)
1036  static const G4ParticleDefinition* gamma = G4Gamma::Definition();
1039 
1040  G4ThreeVector newDirection; // Re-use to avoid memory churn
1041  for (G4int i=0; i<products->entries(); i++) {
1042  G4DynamicParticle* daughter = (*products)[i];
1043  const G4ParticleDefinition* daughterType =
1044  daughter->GetParticleDefinition();
1045  if (daughterType == electron || daughterType == positron ||
1046  daughterType == neutron || daughterType == gamma ||
1047  daughterType == alpha || daughterType == proton) CollimateDecayProduct(daughter);
1048  }
1049 }
1050 
1052 #ifdef G4VERBOSE
1053  if (GetVerboseLevel() > 1) {
1054  G4cout << "CollimateDecayProduct for daughter "
1055  << daughter->GetParticleDefinition()->GetParticleName() << G4endl;
1056  }
1057 #endif
1058 
1060  if (origin != collimate) daughter->SetMomentumDirection(collimate);
1061 }
1062 
1063 
1064 // Choose random direction within collimation cone
1065 
1067  if (origin == forceDecayDirection) return origin; // Don't do collimation
1068  if (forceDecayHalfAngle == 180.*deg) return origin;
1069 
1071 
1072  // Return direction offset by random throw
1073  if (forceDecayHalfAngle > 0.) {
1074  // Generate uniform direction around central axis
1075  G4double phi = 2.*pi*G4UniformRand();
1076  G4double cosMin = std::cos(forceDecayHalfAngle);
1077  G4double cosTheta = (1.-cosMin)*G4UniformRand() + cosMin; // [cosMin,1.)
1078 
1079  dir.setPhi(dir.phi()+phi);
1080  dir.setTheta(dir.theta()+std::acos(cosTheta));
1081  }
1082 
1083 #ifdef G4VERBOSE
1084  if (GetVerboseLevel()>1)
1085  G4cout << " ChooseCollimationDirection returns " << dir << G4endl;
1086 #endif
1087 
1088  return dir;
1089 }
1090 
Float_t x
Definition: compare.C:6
virtual void ProcessDescription(std::ostream &outFile) const
static G4GenericIon * GenericIon()
Definition: G4GenericIon.cc:93
G4double GetLocalTime() const
G4DecayTable * GetDecayTable(const G4ParticleDefinition *)
G4double GetBR() const
static G4Gamma * Definition()
Definition: G4Gamma.cc:49
static G4Proton * Definition()
Definition: G4Proton.cc:49
G4int GetAMax() const
void BuildPhysicsTable(const G4ParticleDefinition &)
virtual void SetICM(G4bool)
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:45
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4LogicalVolume * GetLogicalVolume() const
static G4Alpha * Definition()
Definition: G4Alpha.cc:49
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double GetPDGLifeTime() const
static G4Neutron * Definition()
Definition: G4Neutron.cc:54
void RegisterExtraProcess(G4VProcess *)
G4DeexPrecoParameters * GetParameters()
static constexpr double keV
Definition: G4SIunits.hh:216
void ProposeWeight(G4double finalWeight)
G4ThreeVector ChooseCollimationDirection() const
void AddSecondary(G4Track *aSecondary)
#define G4endl
Definition: G4ios.hh:61
Float_t y
Definition: compare.C:6
const G4ThreeVector & GetMomentumDirection() const
static G4NuclearLevelData * GetInstance()
G4RadioactiveDecayBaseMessenger * theRadioactiveDecayBaseMessenger
Double_t z
G4int GetAMin() const
const G4String & GetParticleName() const
const G4String & GetParticleType() const
const G4TouchableHandle & GetTouchableHandle() const
G4int GetZMax() const
G4PhotonEvaporation * photonEvaporation
G4VParticleChange * DecayIt(const G4Track &theTrack, const G4Step &theStep)
void setTheta(double)
void Insert(G4VDecayChannel *aChannel)
Definition: G4DecayTable.cc:62
G4bool IsChecked() const
G4double GetMeanLifeTime(const G4Track &theTrack, G4ForceCondition *condition)
G4DynamicParticle * PopProducts()
G4double GetPDGMass() const
void Boost(G4double totalEnergy, const G4ThreeVector &momentumDirection)
double theta() const
G4FloatLevelBase
Definition: G4Ions.hh:95
void SetMomentumDirection(const G4ThreeVector &aDirection)
G4double GetWeight() const
std::map< G4int, G4String > theUserRadioactiveDataFiles
G4int entries() const
Float_t Z
void SelectAVolume(const G4String aVolume)
static constexpr double m
Definition: G4SIunits.hh:129
const XML_Char * s
Definition: expat.h:262
void SetBR(G4double value)
static const G4ThreeVector origin
Definition: G4Ions.hh:51
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:88
G4ParticleDefinition * GetDefinition() const
G4VParticleChange * pParticleChange
Definition: G4VProcess.hh:283
static constexpr double deg
Definition: G4SIunits.hh:152
void CollimateDecay(G4DecayProducts *products)
void SetARM(G4bool onoff)
Definition: G4ITDecay.hh:59
void DumpInfo() const
G4VDecayChannel * GetDecayChannel(G4int index) const
std::vector< G4String > ValidVolumes
G4DecayProducts * DoDecay(const G4ParticleDefinition &theParticleDef)
static const G4double alpha
void DumpInfo() const
#define G4UniformRand()
Definition: Randomize.hh:53
G4double GetGlobalTime() const
G4String strip(G4int strip_Type=trailing, char c=' ')
double A(double temperature)
void SetTouchableHandle(const G4TouchableHandle &apValue)
const G4ThreeVector & GetPosition() const
static constexpr double eV
G4int GetZMin() const
Definition: G4Step.hh:76
G4ParticleChangeForRadDecay fParticleChangeForRadDecay
G4RadioactiveDecayMode
#define noFloat
Definition: G4Ions.hh:118
G4GLOB_DLL std::ostream G4cerr
void setPhi(double)
G4bool IsApplicable(const G4ParticleDefinition &)
G4bool GetPDGStable() const
static G4Electron * Definition()
Definition: G4Electron.cc:49
void DeselectAVolume(const G4String aVolume)
static G4HadronicProcessStore * Instance()
void SetGoodForTrackingFlag(G4bool value=true)
static const G4double levelTolerance
void SetARM(G4bool onoff)
Definition: G4ECDecay.hh:56
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
int G4int
Definition: G4Types.hh:78
double phi() const
static constexpr double c_light
TDirectory * dir
#define DBL_MIN
Definition: templates.hh:75
void AddUserDecayDataFile(G4int Z, G4int A, G4String filename)
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:435
G4BetaDecayType
G4ForceCondition
TFile * file
G4double GetKineticEnergy() const
static constexpr double cm
Definition: G4SIunits.hh:119
G4GLOB_DLL std::ostream G4cout
void CollimateDecayProduct(G4DynamicParticle *product)
static constexpr double deg
void ClearNumberOfInteractionLengthLeft()
Definition: G4VProcess.hh:450
virtual G4DecayProducts * DecayIt(G4double parentMass=-1.0)=0
virtual void RDMForced(G4bool)
std::map< G4String, G4DecayTable * > DecayTableMap
G4VDecayChannel * SelectADecayChannel(G4double parentMass=-1.)
Definition: G4DecayTable.cc:85
G4double GetMass() const
G4double GetTotalMomentum() const
static constexpr double pi
Definition: G4SIunits.hh:75
const G4ParticleDefinition * GetParticleDefinition() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4TrackStatus GetTrackStatus() const
G4DecayTable * LoadDecayTable(const G4ParticleDefinition &theParentNucleus)
static G4Positron * Definition()
Definition: G4Positron.cc:49
static G4Ions::G4FloatLevelBase FloatLevelBase(char flbChar)
Definition: G4Ions.hh:189
G4RadioactiveDecayMode GetDecayMode()
virtual void Initialize(const G4Track &)
#define DBL_MAX
Definition: templates.hh:83
G4RadioactiveDecayBase(const G4String &processName="RadioactiveDecayBase")
void ProposeTrackStatus(G4TrackStatus status)
static constexpr double GeV
Definition: G4SIunits.hh:217
#define ns
Definition: xmlparse.cc:614
const G4DynamicParticle * GetDynamicParticle() const
static G4LogicalVolumeStore * GetInstance()
void RegisterParticleForExtraProcess(G4VProcess *, const G4ParticleDefinition *)
G4double GetMeanFreePath(const G4Track &theTrack, G4double previousStepSize, G4ForceCondition *condition)
const G4String & GetName() const
G4VPhysicalVolume * GetVolume() const
std::mutex G4Mutex
Definition: G4Threading.hh:84
void SetNumberOfSecondaries(G4int totSecondaries)
G4int entries() const