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