Geant4  10.00.p03
LBE.icc
Go to the documentation of this file.
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 //
29 // For information related to this code contact: Alex Howard
30 // e-mail: alexander.howard@cern.ch
31 // --------------------------------------------------------------
32 // Comments
33 //
34 // Underground Advanced
35 //
36 // This physics list is taken from the underground_physics example with small
37 // modifications. It is an example of a "flat" physics list with no dependence
38 // on builders. The physics covered would be suitable for a low background
39 // experiment including the neutron_hp package
40 //
41 //
42 //
43 // PhysicsList program
44 //
45 // Modified:
46 //
47 // 14-02-03 Fix bugs in msc and hIon instanciation + cut per region
48 // 16-08-10 Remove inclusion of obsolete class of G4ParticleWithCuts
49 // 20-10-10 Migrate LowEnergy process to Livermore models, LP
50 // 28-03-13 Replace LEP/HEP with FTFP+BERT (A.R.)
51 // --------------------------------------------------------------
52 
53 #include <iomanip>
54 
55 #include "globals.hh"
56 #include "G4ios.hh"
57 #include "G4ProcessManager.hh"
58 #include "G4ProcessVector.hh"
59 
60 #include "G4ParticleTypes.hh"
61 #include "G4ParticleTable.hh"
62 #include "G4ProductionCutsTable.hh"
63 
64 #include "G4UserLimits.hh"
65 #include "G4DataQuestionaire.hh"
66 #include "G4WarnPLStatus.hh"
67 
68 // Builder for all stopping processes
69 #include "G4StoppingPhysics.hh"
70 
71 // Constructor /////////////////////////////////////////////////////////////
72 template<class T> TLBE<T>::TLBE() :T()
73 {
74 
75  G4DataQuestionaire it(photon, lowenergy, neutron, radioactive);
76  G4cout << "You are using the simulation engine: LBE 5.3"<<G4endl;
77  G4cout <<G4endl<<G4endl;
78  this->defaultCutValue = 1.0*CLHEP::micrometer; //
79  cutForGamma = this->defaultCutValue;
80  cutForElectron = 1.0*CLHEP::nanometer;
81  cutForPositron = this->defaultCutValue;
82  cutForProton = this->defaultCutValue;
83  cutForAlpha = 1.0*CLHEP::nanometer;
84  cutForGenericIon = 1.0*CLHEP::nanometer;
85 
86  stoppingPhysics = new G4StoppingPhysics;
87 
88  VerboseLevel = 1;
89  OpVerbLevel = 0;
90 
91  this->SetVerboseLevel(VerboseLevel);
92 }
93 
94 
95 // Destructor //////////////////////////////////////////////////////////////
96 template<class T> TLBE<T>::~TLBE()
97 {
98  delete stoppingPhysics;
99 }
100 
101 
102 // Construct Particles /////////////////////////////////////////////////////
103  template<class T> void TLBE<T>::ConstructParticle()
104 {
105 
106  // In this method, static member functions should be called
107  // for all particles which you want to use.
108  // This ensures that objects of these particle types will be
109  // created in the program.
110 
111  ConstructMyBosons();
112  ConstructMyLeptons();
113  ConstructMyMesons();
114  ConstructMyBaryons();
115  ConstructMyIons();
116  ConstructMyShortLiveds();
117  stoppingPhysics->ConstructParticle(); // Anything not included above
118 }
119 
120 
121 // construct Bosons://///////////////////////////////////////////////////
122  template<class T> void TLBE<T>::ConstructMyBosons()
123 {
124  // pseudo-particles
125  G4Geantino::GeantinoDefinition();
126  G4ChargedGeantino::ChargedGeantinoDefinition();
127 
128  // gamma
129  G4Gamma::GammaDefinition();
130 
131  //OpticalPhotons
132  G4OpticalPhoton::OpticalPhotonDefinition();
133 }
134 
135 
136 // construct Leptons://///////////////////////////////////////////////////
137  template<class T> void TLBE<T>::ConstructMyLeptons()
138 {
139  // leptons
140  G4Electron::ElectronDefinition();
141  G4Positron::PositronDefinition();
142  G4MuonPlus::MuonPlusDefinition();
143  G4MuonMinus::MuonMinusDefinition();
144 
145  G4NeutrinoE::NeutrinoEDefinition();
146  G4AntiNeutrinoE::AntiNeutrinoEDefinition();
147  G4NeutrinoMu::NeutrinoMuDefinition();
148  G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
149 }
150 
151 #include "G4MesonConstructor.hh"
152 #include "G4BaryonConstructor.hh"
153 #include "G4IonConstructor.hh"
154 
155 
156 // construct Mesons://///////////////////////////////////////////////////
157  template<class T> void TLBE<T>::ConstructMyMesons()
158 {
159  // mesons
160  G4MesonConstructor mConstructor;
161  mConstructor.ConstructParticle();
162 
163 }
164 
165 
166 // construct Baryons://///////////////////////////////////////////////////
167  template<class T> void TLBE<T>::ConstructMyBaryons()
168 {
169  // baryons
170  G4BaryonConstructor bConstructor;
171  bConstructor.ConstructParticle();
172 
173 }
174 
175 
176 // construct Ions://///////////////////////////////////////////////////
177  template<class T> void TLBE<T>::ConstructMyIons()
178 {
179  // ions
180  G4IonConstructor iConstructor;
181  iConstructor.ConstructParticle();
182 
183 }
184 
185 // construct Shortliveds://///////////////////////////////////////////////////
186  template<class T> void TLBE<T>::ConstructMyShortLiveds()
187 {
188  // ShortLiveds
189  ;
190 }
191 
192 
193 
194 
195 // Construct Processes //////////////////////////////////////////////////////
196  template<class T> void TLBE<T>::ConstructProcess()
197 {
198  AddTransportation();
199  ConstructEM();
200  ConstructOp();
201  ConstructHad();
202  ConstructGeneral();
203 }
204 
205 
206 // Transportation ///////////////////////////////////////////////////////////
207 #include "G4MaxTimeCuts.hh"
208 #include "G4MinEkineCuts.hh"
209 
210  template<class T> void TLBE<T>::AddTransportation() {
211 
212  G4VUserPhysicsList::AddTransportation();
213 
214  theParticleIterator->reset();
215  while( (*(theParticleIterator))() ){
216  G4ParticleDefinition* particle = theParticleIterator->value();
217  G4ProcessManager* pmanager = particle->GetProcessManager();
218  G4String particleName = particle->GetParticleName();
219  // time cuts for ONLY neutrons:
220  if(particleName == "neutron")
221  pmanager->AddDiscreteProcess(new G4MaxTimeCuts());
222  // Energy cuts to kill charged (embedded in method) particles:
223  pmanager->AddDiscreteProcess(new G4MinEkineCuts());
224  }
225 }
226 
227 
228 // Electromagnetic Processes ////////////////////////////////////////////////
229 // all charged particles
230 
231 #include "G4eMultipleScattering.hh"
232 #include "G4MuMultipleScattering.hh"
233 #include "G4hMultipleScattering.hh"
234 
235 // gamma. Use Livermore models
236 #include "G4PhotoElectricEffect.hh"
237 #include "G4LivermorePhotoElectricModel.hh"
238 #include "G4ComptonScattering.hh"
239 #include "G4LivermoreComptonModel.hh"
240 #include "G4GammaConversion.hh"
241 #include "G4LivermoreGammaConversionModel.hh"
242 #include "G4RayleighScattering.hh"
243 #include "G4LivermoreRayleighModel.hh"
244 
245 
246 // e-
247 #include "G4eMultipleScattering.hh"
248 #include "G4UniversalFluctuation.hh"
249 #include "G4UrbanMscModel.hh"
250 
251 #include "G4eIonisation.hh"
252 #include "G4LivermoreIonisationModel.hh"
253 
254 #include "G4eBremsstrahlung.hh"
255 #include "G4LivermoreBremsstrahlungModel.hh"
256 
257 // e+
258 #include "G4eplusAnnihilation.hh"
259 
260 
261 // alpha and GenericIon and deuterons, triton, He3:
262 #include "G4ionIonisation.hh"
263 #include "G4hIonisation.hh"
264 #include "G4hBremsstrahlung.hh"
265 //
266 #include "G4IonParametrisedLossModel.hh"
267 #include "G4NuclearStopping.hh"
268 #include "G4EnergyLossTables.hh"
269 
270 //muon:
271 #include "G4MuIonisation.hh"
272 #include "G4MuBremsstrahlung.hh"
273 #include "G4MuPairProduction.hh"
274 #include "G4MuonMinusCaptureAtRest.hh"
275 
276 //OTHERS:
277 //#include "G4hIonisation.hh" // standard hadron ionisation
278 
279  template<class T> void TLBE<T>::ConstructEM() {
280 
281  // models & processes:
282  // Use Livermore models up to 20 MeV, and standard
283  // models for higher energy
284  G4double LivermoreHighEnergyLimit = 20*CLHEP::MeV;
285  //
286  theParticleIterator->reset();
287  while( (*(theParticleIterator))() ){
288  G4ParticleDefinition* particle = theParticleIterator->value();
289  G4ProcessManager* pmanager = particle->GetProcessManager();
290  G4String particleName = particle->GetParticleName();
291  G4String particleType = particle->GetParticleType();
292  G4double charge = particle->GetPDGCharge();
293 
294  if (particleName == "gamma")
295  {
296  G4PhotoElectricEffect* thePhotoElectricEffect = new G4PhotoElectricEffect();
297  G4LivermorePhotoElectricModel* theLivermorePhotoElectricModel =
298  new G4LivermorePhotoElectricModel();
299  theLivermorePhotoElectricModel->SetHighEnergyLimit(LivermoreHighEnergyLimit);
300  thePhotoElectricEffect->AddEmModel(0, theLivermorePhotoElectricModel);
301  pmanager->AddDiscreteProcess(thePhotoElectricEffect);
302 
303  G4ComptonScattering* theComptonScattering = new G4ComptonScattering();
304  G4LivermoreComptonModel* theLivermoreComptonModel =
305  new G4LivermoreComptonModel();
306  theLivermoreComptonModel->SetHighEnergyLimit(LivermoreHighEnergyLimit);
307  theComptonScattering->AddEmModel(0, theLivermoreComptonModel);
308  pmanager->AddDiscreteProcess(theComptonScattering);
309 
310  G4GammaConversion* theGammaConversion = new G4GammaConversion();
311  G4LivermoreGammaConversionModel* theLivermoreGammaConversionModel =
312  new G4LivermoreGammaConversionModel();
313  theLivermoreGammaConversionModel->SetHighEnergyLimit(LivermoreHighEnergyLimit);
314  theGammaConversion->AddEmModel(0, theLivermoreGammaConversionModel);
315  pmanager->AddDiscreteProcess(theGammaConversion);
316 
317  G4RayleighScattering* theRayleigh = new G4RayleighScattering();
318  G4LivermoreRayleighModel* theRayleighModel = new G4LivermoreRayleighModel();
319  theRayleighModel->SetHighEnergyLimit(LivermoreHighEnergyLimit);
320  theRayleigh->AddEmModel(0, theRayleighModel);
321  pmanager->AddDiscreteProcess(theRayleigh);
322 
323  }
324  else if (particleName == "e-")
325  {
326  //electron
327  // process ordering: AddProcess(name, at rest, along step, post step)
328  // -1 = not implemented, then ordering
329  G4eMultipleScattering* msc = new G4eMultipleScattering();
330  //msc->AddEmModel(0, new G4UrbanMscModel());
331  msc->SetStepLimitType(fUseDistanceToBoundary);
332  pmanager->AddProcess(msc, -1, 1, 1);
333 
334  // Ionisation
335  G4eIonisation* eIoni = new G4eIonisation();
336  G4LivermoreIonisationModel* theIoniLivermore = new
337  G4LivermoreIonisationModel();
338  theIoniLivermore->SetHighEnergyLimit(1*CLHEP::MeV);
339  eIoni->AddEmModel(0, theIoniLivermore, new G4UniversalFluctuation() );
340  eIoni->SetStepFunction(0.2, 100*CLHEP::um); //
341  pmanager->AddProcess(eIoni, -1, 2, 2);
342 
343  // Bremsstrahlung
344  G4eBremsstrahlung* eBrem = new G4eBremsstrahlung();
345  G4LivermoreBremsstrahlungModel* theBremLivermore = new
346  G4LivermoreBremsstrahlungModel();
347  theBremLivermore->SetHighEnergyLimit(LivermoreHighEnergyLimit);
348  eBrem->AddEmModel(0, theBremLivermore);
349  pmanager->AddProcess(eBrem, -1,-3, 3);
350  }
351  else if (particleName == "e+")
352  {
353  //positron
354  G4eMultipleScattering* msc = new G4eMultipleScattering();
355  //msc->AddEmModel(0, new G4UrbanMscModel());
356  msc->SetStepLimitType(fUseDistanceToBoundary);
357  pmanager->AddProcess(msc, -1, 1, 1);
358  G4eIonisation* eIoni = new G4eIonisation();
359  eIoni->SetStepFunction(0.2, 100*CLHEP::um);
360  pmanager->AddProcess(eIoni, -1, 2, 2);
361  pmanager->AddProcess(new G4eBremsstrahlung, -1,-3, 3);
362  pmanager->AddProcess(new G4eplusAnnihilation,0,-1, 4);
363  }
364  else if( particleName == "mu+" ||
365  particleName == "mu-" )
366  {
367  //muon
368  G4MuMultipleScattering* aMultipleScattering = new G4MuMultipleScattering();
369  pmanager->AddProcess(aMultipleScattering, -1, 1, 1);
370  pmanager->AddProcess(new G4MuIonisation(), -1, 2, 2);
371  pmanager->AddProcess(new G4MuBremsstrahlung(), -1,-1, 3);
372  pmanager->AddProcess(new G4MuPairProduction(), -1,-1, 4);
373  if( particleName == "mu-" )
374  pmanager->AddProcess(new G4MuonMinusCaptureAtRest(), 0,-1,-1);
375  }
376  else if (particleName == "GenericIon")
377  {
378  pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
379  G4ionIonisation* ionIoni = new G4ionIonisation();
380  ionIoni->SetEmModel(new G4IonParametrisedLossModel());
381  ionIoni->SetStepFunction(0.1, 10*CLHEP::um);
382  pmanager->AddProcess(ionIoni, -1, 2, 2);
383  pmanager->AddProcess(new G4NuclearStopping(), -1, 3,-1);
384  }
385  else if (particleName == "alpha" || particleName == "He3")
386  {
387  //MSC, ion-Ionisation, Nuclear Stopping
388  pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
389 
390  G4ionIonisation* ionIoni = new G4ionIonisation();
391  ionIoni->SetStepFunction(0.1, 20*CLHEP::um);
392  pmanager->AddProcess(ionIoni, -1, 2, 2);
393  pmanager->AddProcess(new G4NuclearStopping(), -1, 3,-1);
394  }
395  else if (particleName == "proton" ||
396  particleName == "deuteron" ||
397  particleName == "triton" ||
398  particleName == "pi+" ||
399  particleName == "pi-" ||
400  particleName == "kaon+" ||
401  particleName == "kaon-")
402  {
403  //MSC, h-ionisation, bremsstrahlung
404  pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
405  G4hIonisation* hIoni = new G4hIonisation();
406  hIoni->SetStepFunction(0.2, 50*CLHEP::um);
407  pmanager->AddProcess(hIoni, -1, 2, 2);
408  pmanager->AddProcess(new G4hBremsstrahlung, -1,-3, 3);
409  }
410  else if ((!particle->IsShortLived()) &&
411  (charge != 0.0) &&
412  (particle->GetParticleName() != "chargedgeantino"))
413  {
414  //all others charged particles except geantino
415  pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
416  pmanager->AddProcess(new G4hIonisation, -1, 2, 2);
417  }
418 
419  }
420 }
421 
422 
423 // Optical Processes ////////////////////////////////////////////////////////
424 #include "G4Scintillation.hh"
425 #include "G4OpAbsorption.hh"
426 //#include "G4OpRayleigh.hh"
427 #include "G4OpBoundaryProcess.hh"
428 
429  template<class T> void TLBE<T>::ConstructOp()
430 {
431  // default scintillation process
432  //Coverity report: check that the process is actually used, if not must delete
433  G4bool theScintProcessDefNeverUsed = true;
434  G4Scintillation* theScintProcessDef = new G4Scintillation("Scintillation");
435  // theScintProcessDef->DumpPhysicsTable();
436  theScintProcessDef->SetTrackSecondariesFirst(true);
437  theScintProcessDef->SetScintillationYieldFactor(1.0); //
438  theScintProcessDef->SetScintillationExcitationRatio(0.0); //
439  theScintProcessDef->SetVerboseLevel(OpVerbLevel);
440 
441  // scintillation process for alpha:
442  G4bool theScintProcessAlphaNeverUsed = true;
443  G4Scintillation* theScintProcessAlpha = new G4Scintillation("Scintillation");
444  // theScintProcessNuc->DumpPhysicsTable();
445  theScintProcessAlpha->SetTrackSecondariesFirst(true);
446  theScintProcessAlpha->SetScintillationYieldFactor(1.1);
447  theScintProcessAlpha->SetScintillationExcitationRatio(1.0);
448  theScintProcessAlpha->SetVerboseLevel(OpVerbLevel);
449 
450  // scintillation process for heavy nuclei
451  G4bool theScintProcessNucNeverUsed = true;
452  G4Scintillation* theScintProcessNuc = new G4Scintillation("Scintillation");
453  // theScintProcessNuc->DumpPhysicsTable();
454  theScintProcessNuc->SetTrackSecondariesFirst(true);
455  theScintProcessNuc->SetScintillationYieldFactor(0.2);
456  theScintProcessNuc->SetScintillationExcitationRatio(1.0);
457  theScintProcessNuc->SetVerboseLevel(OpVerbLevel);
458 
459  // optical processes
460  G4bool theAbsorptionProcessNeverUsed = true;
461  G4OpAbsorption* theAbsorptionProcess = new G4OpAbsorption();
462  // G4OpRayleigh* theRayleighScatteringProcess = new G4OpRayleigh();
463  G4bool theBoundaryProcessNeverUsed = true;
464  G4OpBoundaryProcess* theBoundaryProcess = new G4OpBoundaryProcess();
465  // theAbsorptionProcess->DumpPhysicsTable();
466  // theRayleighScatteringProcess->DumpPhysicsTable();
467  theAbsorptionProcess->SetVerboseLevel(OpVerbLevel);
468  // theRayleighScatteringProcess->SetVerboseLevel(OpVerbLevel);
469  theBoundaryProcess->SetVerboseLevel(OpVerbLevel);
470 
471  theParticleIterator->reset();
472  while( (*(theParticleIterator))() )
473  {
474  G4ParticleDefinition* particle = theParticleIterator->value();
475  G4ProcessManager* pmanager = particle->GetProcessManager();
476  G4String particleName = particle->GetParticleName();
477  if (theScintProcessDef->IsApplicable(*particle)) {
478  // if(particle->GetPDGMass() > 5.0*CLHEP::GeV)
479  if(particle->GetParticleName() == "GenericIon") {
480  pmanager->AddProcess(theScintProcessNuc); // AtRestDiscrete
481  pmanager->SetProcessOrderingToLast(theScintProcessNuc,idxAtRest);
482  pmanager->SetProcessOrderingToLast(theScintProcessNuc,idxPostStep);
483  theScintProcessNucNeverUsed = false;
484  }
485  else if(particle->GetParticleName() == "alpha") {
486  pmanager->AddProcess(theScintProcessAlpha);
487  pmanager->SetProcessOrderingToLast(theScintProcessAlpha,idxAtRest);
488  pmanager->SetProcessOrderingToLast(theScintProcessAlpha,idxPostStep);
489  theScintProcessAlphaNeverUsed = false;
490  }
491  else {
492  pmanager->AddProcess(theScintProcessDef);
493  pmanager->SetProcessOrderingToLast(theScintProcessDef,idxAtRest);
494  pmanager->SetProcessOrderingToLast(theScintProcessDef,idxPostStep);
495  theScintProcessDefNeverUsed = false;
496  }
497  }
498 
499  if (particleName == "opticalphoton") {
500  pmanager->AddDiscreteProcess(theAbsorptionProcess);
501  theAbsorptionProcessNeverUsed = false;
502  // pmanager->AddDiscreteProcess(theRayleighScatteringProcess);
503  theBoundaryProcessNeverUsed = false;
504  pmanager->AddDiscreteProcess(theBoundaryProcess);
505  }
506  }
507  if ( theScintProcessDefNeverUsed ) delete theScintProcessDef;
508  if ( theScintProcessAlphaNeverUsed ) delete theScintProcessAlpha;
509  if ( theScintProcessNucNeverUsed ) delete theScintProcessNuc;
510  if ( theBoundaryProcessNeverUsed ) delete theBoundaryProcess;
511  if ( theAbsorptionProcessNeverUsed ) delete theAbsorptionProcess;
512 }
513 
514 
515 // Hadronic processes ////////////////////////////////////////////////////////
516 
517 // Elastic processes:
518 #include "G4HadronElasticProcess.hh"
519 #include "G4HadronCaptureProcess.hh"
520 #include "G4HadronElastic.hh"
521 #include "G4ChipsElasticModel.hh"
522 #include "G4ElasticHadrNucleusHE.hh"
523 #include "G4AntiNuclElastic.hh"
524 #include "G4BGGNucleonElasticXS.hh"
525 #include "G4BGGPionElasticXS.hh"
526 #include "G4NeutronElasticXS.hh"
527 #include "G4CrossSectionDataSetRegistry.hh"
528 #include "G4ChipsProtonElasticXS.hh"
529 #include "G4ChipsNeutronElasticXS.hh"
530 #include "G4ComponentAntiNuclNuclearXS.hh"
531 #include "G4CrossSectionElastic.hh"
532 
533 // Inelastic processes:
534 #include "G4PionPlusInelasticProcess.hh"
535 #include "G4PionMinusInelasticProcess.hh"
536 #include "G4KaonPlusInelasticProcess.hh"
537 #include "G4KaonZeroSInelasticProcess.hh"
538 #include "G4KaonZeroLInelasticProcess.hh"
539 #include "G4KaonMinusInelasticProcess.hh"
540 #include "G4ProtonInelasticProcess.hh"
541 #include "G4AntiProtonInelasticProcess.hh"
542 #include "G4NeutronInelasticProcess.hh"
543 #include "G4AntiNeutronInelasticProcess.hh"
544 #include "G4DeuteronInelasticProcess.hh"
545 #include "G4TritonInelasticProcess.hh"
546 #include "G4AlphaInelasticProcess.hh"
547 
548 // FTFP + BERT model
549 #include "G4TheoFSGenerator.hh"
550 #include "G4ExcitationHandler.hh"
551 #include "G4PreCompoundModel.hh"
552 #include "G4GeneratorPrecompoundInterface.hh"
553 #include "G4FTFModel.hh"
554 #include "G4LundStringFragmentation.hh"
555 #include "G4ExcitedStringDecay.hh"
556 #include "G4CascadeInterface.hh"
557 #include "G4CrossSectionInelastic.hh"
558 #include "G4PiNuclearCrossSection.hh"
559 #include "G4CrossSectionPairGG.hh"
560 #include "G4ChipsKaonMinusInelasticXS.hh"
561 #include "G4ChipsKaonPlusInelasticXS.hh"
562 #include "G4ChipsKaonZeroInelasticXS.hh"
563 #include "G4CrossSectionDataSetRegistry.hh"
564 #include "G4BGGNucleonInelasticXS.hh"
565 #include "G4ComponentAntiNuclNuclearXS.hh"
566 #include "G4GGNuclNuclCrossSection.hh"
567 
568 // Neutron high-precision models: <20 MeV
569 #include "G4NeutronHPElastic.hh"
570 #include "G4NeutronHPElasticData.hh"
571 #include "G4NeutronHPCapture.hh"
572 #include "G4NeutronHPCaptureData.hh"
573 #include "G4NeutronHPInelastic.hh"
574 #include "G4NeutronHPInelasticData.hh"
575 #include "G4NeutronCaptureXS.hh"
576 #include "G4NeutronRadCapture.hh"
577 
578 // ConstructHad()
579 // Makes discrete physics processes for the hadrons, at present limited
580 // to those particles with GHEISHA interactions (INTRC > 0).
581 // The processes are: Elastic scattering and Inelastic scattering.
582 // F.W.Jones 09-JUL-1998
583  template<class T> void TLBE<T>::ConstructHad()
584 {
585  // Elastic scattering
586  const G4double elastic_elimitPi = 1.0*CLHEP::GeV;
587 
588  G4HadronElastic* elastic_lhep0 = new G4HadronElastic();
589  G4HadronElastic* elastic_lhep1 = new G4HadronElastic();
590  elastic_lhep1->SetMaxEnergy( elastic_elimitPi );
591 
592  G4ChipsElasticModel* elastic_chip = new G4ChipsElasticModel();
593 
594  G4ElasticHadrNucleusHE* elastic_he = new G4ElasticHadrNucleusHE();
595  elastic_he->SetMinEnergy( elastic_elimitPi );
596 
597  // Inelastic scattering
598  const G4double theFTFMin0 = 0.0*CLHEP::GeV;
599  const G4double theFTFMin1 = 4.0*CLHEP::GeV;
600  const G4double theFTFMax = 100.0*CLHEP::TeV;
601  const G4double theBERTMin0 = 0.0*CLHEP::GeV;
602  const G4double theBERTMin1 = 19.0*CLHEP::MeV;
603  const G4double theBERTMax = 5.0*CLHEP::GeV;
604  const G4double theHPMin = 0.0*CLHEP::GeV;
605  const G4double theHPMax = 20.0*CLHEP::MeV;
606 
607  G4FTFModel * theStringModel = new G4FTFModel;
608  G4ExcitedStringDecay * theStringDecay = new G4ExcitedStringDecay( new G4LundStringFragmentation );
609  theStringModel->SetFragmentationModel( theStringDecay );
610  G4PreCompoundModel * thePreEquilib = new G4PreCompoundModel( new G4ExcitationHandler );
611  G4GeneratorPrecompoundInterface * theCascade = new G4GeneratorPrecompoundInterface( thePreEquilib );
612 
613  G4TheoFSGenerator * theFTFModel0 = new G4TheoFSGenerator( "FTFP" );
614  theFTFModel0->SetHighEnergyGenerator( theStringModel );
615  theFTFModel0->SetTransport( theCascade );
616  theFTFModel0->SetMinEnergy( theFTFMin0 );
617  theFTFModel0->SetMaxEnergy( theFTFMax );
618 
619  G4TheoFSGenerator * theFTFModel1 = new G4TheoFSGenerator( "FTFP" );
620  theFTFModel1->SetHighEnergyGenerator( theStringModel );
621  theFTFModel1->SetTransport( theCascade );
622  theFTFModel1->SetMinEnergy( theFTFMin1 );
623  theFTFModel1->SetMaxEnergy( theFTFMax );
624 
625  G4CascadeInterface * theBERTModel0 = new G4CascadeInterface;
626  theBERTModel0->SetMinEnergy( theBERTMin0 );
627  theBERTModel0->SetMaxEnergy( theBERTMax );
628 
629  G4CascadeInterface * theBERTModel1 = new G4CascadeInterface;
630  theBERTModel1->SetMinEnergy( theBERTMin1 );
631  theBERTModel1->SetMaxEnergy( theBERTMax );
632 
633  G4VCrossSectionDataSet * thePiData = new G4CrossSectionPairGG( new G4PiNuclearCrossSection, 91*GeV );
634  G4VCrossSectionDataSet * theAntiNucleonData = new G4CrossSectionInelastic( new G4ComponentAntiNuclNuclearXS );
635  G4VCrossSectionDataSet * theGGNuclNuclData = G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4GGNuclNuclCrossSection::Default_Name());
636 
637  theParticleIterator->reset();
638  while ((*(theParticleIterator))())
639  {
640  G4ParticleDefinition* particle = theParticleIterator->value();
641  G4ProcessManager* pmanager = particle->GetProcessManager();
642  G4String particleName = particle->GetParticleName();
643 
644  if (particleName == "pi+")
645  {
646  // Elastic scattering
647  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
648  theElasticProcess->AddDataSet( new G4BGGPionElasticXS( particle ) );
649  theElasticProcess->RegisterMe( elastic_lhep1 );
650  theElasticProcess->RegisterMe( elastic_he );
651  pmanager->AddDiscreteProcess( theElasticProcess );
652  // Inelastic scattering
653  G4PionPlusInelasticProcess* theInelasticProcess = new G4PionPlusInelasticProcess("inelastic");
654  theInelasticProcess->AddDataSet( thePiData );
655  theInelasticProcess->RegisterMe( theFTFModel1 );
656  theInelasticProcess->RegisterMe( theBERTModel0 );
657  pmanager->AddDiscreteProcess( theInelasticProcess );
658  }
659 
660  else if (particleName == "pi-")
661  {
662  // Elastic scattering
663  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
664  theElasticProcess->AddDataSet( new G4BGGPionElasticXS( particle ) );
665  theElasticProcess->RegisterMe( elastic_lhep1 );
666  theElasticProcess->RegisterMe( elastic_he );
667  pmanager->AddDiscreteProcess( theElasticProcess );
668  // Inelastic scattering
669  G4PionMinusInelasticProcess* theInelasticProcess = new G4PionMinusInelasticProcess("inelastic");
670  theInelasticProcess->AddDataSet( thePiData );
671  theInelasticProcess->RegisterMe( theFTFModel1 );
672  theInelasticProcess->RegisterMe( theBERTModel0 );
673  pmanager->AddDiscreteProcess( theInelasticProcess );
674  }
675 
676  else if (particleName == "kaon+")
677  {
678  // Elastic scattering
679  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
680  theElasticProcess->RegisterMe( elastic_lhep0 );
681  pmanager->AddDiscreteProcess( theElasticProcess );
682  // Inelastic scattering
683  G4KaonPlusInelasticProcess* theInelasticProcess = new G4KaonPlusInelasticProcess("inelastic");
684  theInelasticProcess->AddDataSet( G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsKaonPlusInelasticXS::Default_Name()));
685  theInelasticProcess->RegisterMe( theFTFModel1 );
686  theInelasticProcess->RegisterMe( theBERTModel0 );
687  pmanager->AddDiscreteProcess( theInelasticProcess );
688  }
689 
690  else if (particleName == "kaon0S")
691  {
692  // Elastic scattering
693  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
694  theElasticProcess->RegisterMe( elastic_lhep0 );
695  pmanager->AddDiscreteProcess( theElasticProcess );
696  // Inelastic scattering
697  G4KaonZeroSInelasticProcess* theInelasticProcess = new G4KaonZeroSInelasticProcess("inelastic");
698  theInelasticProcess->AddDataSet( G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsKaonZeroInelasticXS::Default_Name()));
699  theInelasticProcess->RegisterMe( theFTFModel1 );
700  theInelasticProcess->RegisterMe( theBERTModel0 );
701  pmanager->AddDiscreteProcess( theInelasticProcess );
702  }
703 
704  else if (particleName == "kaon0L")
705  {
706  // Elastic scattering
707  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
708  theElasticProcess->RegisterMe( elastic_lhep0 );
709  pmanager->AddDiscreteProcess( theElasticProcess );
710  // Inelastic scattering
711  G4KaonZeroLInelasticProcess* theInelasticProcess = new G4KaonZeroLInelasticProcess("inelastic");
712  theInelasticProcess->AddDataSet( G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsKaonZeroInelasticXS::Default_Name()));
713  theInelasticProcess->RegisterMe( theFTFModel1 );
714  theInelasticProcess->RegisterMe( theBERTModel0 );
715  pmanager->AddDiscreteProcess( theInelasticProcess );
716  }
717 
718  else if (particleName == "kaon-")
719  {
720  // Elastic scattering
721  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
722  theElasticProcess->RegisterMe( elastic_lhep0 );
723  pmanager->AddDiscreteProcess( theElasticProcess );
724  // Inelastic scattering
725  G4KaonMinusInelasticProcess* theInelasticProcess = new G4KaonMinusInelasticProcess("inelastic");
726  theInelasticProcess->AddDataSet( G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsKaonMinusInelasticXS::Default_Name()));
727  theInelasticProcess->RegisterMe( theFTFModel1 );
728  theInelasticProcess->RegisterMe( theBERTModel0 );
729  pmanager->AddDiscreteProcess( theInelasticProcess );
730  }
731 
732  else if (particleName == "proton")
733  {
734  // Elastic scattering
735  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
736  theElasticProcess->AddDataSet(G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsProtonElasticXS::Default_Name()));
737  theElasticProcess->RegisterMe( elastic_chip );
738  pmanager->AddDiscreteProcess( theElasticProcess );
739  // Inelastic scattering
740  G4ProtonInelasticProcess* theInelasticProcess = new G4ProtonInelasticProcess("inelastic");
741  theInelasticProcess->AddDataSet( new G4BGGNucleonInelasticXS( G4Proton::Proton() ) );
742  theInelasticProcess->RegisterMe( theFTFModel1 );
743  theInelasticProcess->RegisterMe( theBERTModel0 );
744  pmanager->AddDiscreteProcess( theInelasticProcess );
745  }
746 
747  else if (particleName == "anti_proton")
748  {
749  // Elastic scattering
750  const G4double elastic_elimitAntiNuc = 100.0*CLHEP::MeV;
751  G4AntiNuclElastic* elastic_anuc = new G4AntiNuclElastic();
752  elastic_anuc->SetMinEnergy( elastic_elimitAntiNuc );
753  G4CrossSectionElastic* elastic_anucxs = new G4CrossSectionElastic( elastic_anuc->GetComponentCrossSection() );
754  G4HadronElastic* elastic_lhep2 = new G4HadronElastic();
755  elastic_lhep2->SetMaxEnergy( elastic_elimitAntiNuc );
756  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
757  theElasticProcess->AddDataSet( elastic_anucxs );
758  theElasticProcess->RegisterMe( elastic_lhep2 );
759  theElasticProcess->RegisterMe( elastic_anuc );
760  pmanager->AddDiscreteProcess( theElasticProcess );
761  // Inelastic scattering
762  G4AntiProtonInelasticProcess* theInelasticProcess = new G4AntiProtonInelasticProcess("inelastic");
763  theInelasticProcess->AddDataSet( theAntiNucleonData );
764  theInelasticProcess->RegisterMe( theFTFModel0 );
765  pmanager->AddDiscreteProcess( theInelasticProcess );
766  }
767 
768  else if (particleName == "neutron") {
769  // elastic scattering
770  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
771  theElasticProcess->AddDataSet(G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsNeutronElasticXS::Default_Name()));
772  G4HadronElastic* elastic_neutronChipsModel = new G4ChipsElasticModel();
773  elastic_neutronChipsModel->SetMinEnergy( 19.0*CLHEP::MeV );
774  theElasticProcess->RegisterMe( elastic_neutronChipsModel );
775  G4NeutronHPElastic * theElasticNeutronHP = new G4NeutronHPElastic;
776  theElasticNeutronHP->SetMinEnergy( theHPMin );
777  theElasticNeutronHP->SetMaxEnergy( theHPMax );
778  theElasticProcess->RegisterMe( theElasticNeutronHP );
779  theElasticProcess->AddDataSet( new G4NeutronHPElasticData );
780  pmanager->AddDiscreteProcess( theElasticProcess );
781  // inelastic scattering
782  G4NeutronInelasticProcess* theInelasticProcess = new G4NeutronInelasticProcess("inelastic");
783  theInelasticProcess->AddDataSet( new G4BGGNucleonInelasticXS( G4Neutron::Neutron() ) );
784  theInelasticProcess->RegisterMe( theFTFModel1 );
785  theInelasticProcess->RegisterMe( theBERTModel1 );
786  G4NeutronHPInelastic * theNeutronInelasticHPModel = new G4NeutronHPInelastic;
787  theNeutronInelasticHPModel->SetMinEnergy( theHPMin );
788  theNeutronInelasticHPModel->SetMaxEnergy( theHPMax );
789  theInelasticProcess->RegisterMe( theNeutronInelasticHPModel );
790  theInelasticProcess->AddDataSet( new G4NeutronHPInelasticData );
791  pmanager->AddDiscreteProcess(theInelasticProcess);
792  // capture
793  G4HadronCaptureProcess* theCaptureProcess = new G4HadronCaptureProcess;
794  G4NeutronHPCapture * theNeutronCaptureHPModel = new G4NeutronHPCapture;
795  theNeutronCaptureHPModel->SetMinEnergy( theHPMin );
796  theNeutronCaptureHPModel->SetMaxEnergy( theHPMax );
797  G4NeutronRadCapture* theNeutronRadCapture = new G4NeutronRadCapture();
798  theNeutronRadCapture->SetMinEnergy(theHPMax*0.99);
799  theCaptureProcess->RegisterMe( theNeutronCaptureHPModel );
800  theCaptureProcess->RegisterMe( theNeutronRadCapture);
801  theCaptureProcess->AddDataSet( new G4NeutronHPCaptureData );
802  theCaptureProcess->AddDataSet( new G4NeutronCaptureXS);
803  pmanager->AddDiscreteProcess(theCaptureProcess);
804  }
805  else if (particleName == "anti_neutron")
806  {
807  // Elastic scattering
808  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
809  theElasticProcess->RegisterMe( elastic_lhep0 );
810  pmanager->AddDiscreteProcess( theElasticProcess );
811  // Inelastic scattering
812  G4AntiNeutronInelasticProcess* theInelasticProcess = new G4AntiNeutronInelasticProcess("inelastic");
813  theInelasticProcess->AddDataSet( theAntiNucleonData );
814  theInelasticProcess->RegisterMe( theFTFModel0 );
815  pmanager->AddDiscreteProcess( theInelasticProcess );
816  }
817 
818  else if (particleName == "deuteron")
819  {
820  // Elastic scattering
821  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
822  theElasticProcess->RegisterMe( elastic_lhep0 );
823  pmanager->AddDiscreteProcess( theElasticProcess );
824  // Inelastic scattering
825  G4DeuteronInelasticProcess* theInelasticProcess = new G4DeuteronInelasticProcess("inelastic");
826  theInelasticProcess->AddDataSet( theGGNuclNuclData );
827  theInelasticProcess->RegisterMe( theFTFModel1 );
828  theInelasticProcess->RegisterMe( theBERTModel0 );
829  pmanager->AddDiscreteProcess( theInelasticProcess );
830  }
831 
832  else if (particleName == "triton")
833  {
834  // Elastic scattering
835  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
836  theElasticProcess->RegisterMe( elastic_lhep0 );
837  pmanager->AddDiscreteProcess( theElasticProcess );
838  // Inelastic scattering
839  G4TritonInelasticProcess* theInelasticProcess = new G4TritonInelasticProcess("inelastic");
840  theInelasticProcess->AddDataSet( theGGNuclNuclData );
841  theInelasticProcess->RegisterMe( theFTFModel1 );
842  theInelasticProcess->RegisterMe( theBERTModel0 );
843  pmanager->AddDiscreteProcess( theInelasticProcess );
844  }
845 
846  else if (particleName == "alpha")
847  {
848  // Elastic scattering
849  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
850  theElasticProcess->RegisterMe( elastic_lhep0 );
851  pmanager->AddDiscreteProcess( theElasticProcess );
852  // Inelastic scattering
853  G4AlphaInelasticProcess* theInelasticProcess = new G4AlphaInelasticProcess("inelastic");
854  theInelasticProcess->AddDataSet( theGGNuclNuclData );
855  theInelasticProcess->RegisterMe( theFTFModel1 );
856  theInelasticProcess->RegisterMe( theBERTModel0 );
857  pmanager->AddDiscreteProcess( theInelasticProcess );
858  }
859  } // while ((*(theParticleIterator))())
860 
861  // Add stopping processes with builder
862  stoppingPhysics->ConstructProcess();
863 }
864 
865 
866 // Decays ///////////////////////////////////////////////////////////////////
867 #include "G4Decay.hh"
868 #include "G4RadioactiveDecay.hh"
869 #include "G4IonTable.hh"
870 #include "G4Ions.hh"
871 
872  template<class T> void TLBE<T>::ConstructGeneral() {
873 
874  // Add Decay Process
875  G4Decay* theDecayProcess = new G4Decay();
876  G4bool theDecayProcessNeverUsed = true; //Check if theDecayProcess will be used
877  theParticleIterator->reset();
878  while( (*(theParticleIterator))() )
879  {
880  G4ParticleDefinition* particle = theParticleIterator->value();
881  G4ProcessManager* pmanager = particle->GetProcessManager();
882 
883  if (theDecayProcess->IsApplicable(*particle) && !particle->IsShortLived())
884  {
885  theDecayProcessNeverUsed = false;
886  pmanager ->AddProcess(theDecayProcess);
887  // set ordering for PostStepDoIt and AtRestDoIt
888  pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
889  pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
890  }
891  }
892 
893  // Declare radioactive decay to the GenericIon in the IonTable.
894  const G4IonTable *theIonTable =
895  G4ParticleTable::GetParticleTable()->GetIonTable();
896  G4RadioactiveDecay *theRadioactiveDecay = new G4RadioactiveDecay();
897 
898  for (G4int i=0; i<theIonTable->Entries(); i++)
899  {
900  G4String particleName = theIonTable->GetParticle(i)->GetParticleName();
901  G4String particleType = theIonTable->GetParticle(i)->GetParticleType();
902 
903  if (particleName == "GenericIon")
904  {
905  G4ProcessManager* pmanager =
906  theIonTable->GetParticle(i)->GetProcessManager();
907  pmanager->SetVerboseLevel(VerboseLevel);
908  pmanager ->AddProcess(theRadioactiveDecay);
909  pmanager ->SetProcessOrdering(theRadioactiveDecay, idxPostStep);
910  pmanager ->SetProcessOrdering(theRadioactiveDecay, idxAtRest);
911  }
912  }
913  //If we actually never used the process, delete it
914  //From Coverity report
915  if ( theDecayProcessNeverUsed ) delete theDecayProcess;
916 }
917 
918 // Cuts /////////////////////////////////////////////////////////////////////
919 template<class T> void TLBE<T>::SetCuts()
920 {
921 
922  if (this->verboseLevel >1)
923  G4cout << "LBE::SetCuts:";
924 
925  if (this->verboseLevel>0){
926  G4cout << "LBE::SetCuts:";
927  G4cout << "CutLength : "
928  << G4BestUnit(this->defaultCutValue,"Length") << G4endl;
929  }
930 
931  //special for low energy physics
932  G4double lowlimit=250*CLHEP::eV;
933  G4ProductionCutsTable * aPCTable = G4ProductionCutsTable::GetProductionCutsTable();
934  aPCTable->SetEnergyRange(lowlimit,100*CLHEP::GeV);
935 
936  // set cut values for gamma at first and for e- second and next for e+,
937  // because some processes for e+/e- need cut values for gamma
938  this->SetCutValue(cutForGamma, "gamma");
939  this->SetCutValue(cutForElectron, "e-");
940  this->SetCutValue(cutForPositron, "e+");
941 
942  // this->SetCutValue(cutForProton, "proton");
943  // this->SetCutValue(cutForProton, "anti_proton");
944  // this->SetCutValue(cutForAlpha, "alpha");
945  // this->SetCutValue(cutForGenericIon, "GenericIon");
946 
947  // this->SetCutValueForOthers(this->defaultCutValue);
948 
949  if (this->verboseLevel>0) this->DumpCutValuesTable();
950 }
951