Geant4  10.02.p02
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 "G4ComponentGGNuclNuclXsc.hh"
567 
568 // Neutron high-precision models: <20 MeV
569 #include "G4ParticleHPElastic.hh"
570 #include "G4ParticleHPElasticData.hh"
571 #include "G4ParticleHPCapture.hh"
572 #include "G4ParticleHPCaptureData.hh"
573 #include "G4ParticleHPInelastic.hh"
574 #include "G4ParticleHPInelasticData.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  G4ComponentGGNuclNuclXsc * ggNuclNuclXsec = new G4ComponentGGNuclNuclXsc();
647  G4VCrossSectionDataSet * theGGNuclNuclData = new G4CrossSectionInelastic(ggNuclNuclXsec);
648 
649  theParticleIterator->reset();
650  while ((*(theParticleIterator))())
651  {
652  G4ParticleDefinition* particle = theParticleIterator->value();
653  G4ProcessManager* pmanager = particle->GetProcessManager();
654  G4String particleName = particle->GetParticleName();
655 
656  if (particleName == "pi+")
657  {
658  // Elastic scattering
659  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
660  theElasticProcess->AddDataSet( new G4BGGPionElasticXS( particle ) );
661  theElasticProcess->RegisterMe( elastic_lhep1 );
662  theElasticProcess->RegisterMe( elastic_he );
663  pmanager->AddDiscreteProcess( theElasticProcess );
664  // Inelastic scattering
665  G4PionPlusInelasticProcess* theInelasticProcess = new G4PionPlusInelasticProcess("inelastic");
666  theInelasticProcess->AddDataSet( thePiData );
667  theInelasticProcess->RegisterMe( theFTFModel1 );
668  theInelasticProcess->RegisterMe( theBERTModel0 );
669  pmanager->AddDiscreteProcess( theInelasticProcess );
670  }
671 
672  else if (particleName == "pi-")
673  {
674  // Elastic scattering
675  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
676  theElasticProcess->AddDataSet( new G4BGGPionElasticXS( particle ) );
677  theElasticProcess->RegisterMe( elastic_lhep1 );
678  theElasticProcess->RegisterMe( elastic_he );
679  pmanager->AddDiscreteProcess( theElasticProcess );
680  // Inelastic scattering
681  G4PionMinusInelasticProcess* theInelasticProcess = new G4PionMinusInelasticProcess("inelastic");
682  theInelasticProcess->AddDataSet( thePiData );
683  theInelasticProcess->RegisterMe( theFTFModel1 );
684  theInelasticProcess->RegisterMe( theBERTModel0 );
685  pmanager->AddDiscreteProcess( theInelasticProcess );
686  }
687 
688  else if (particleName == "kaon+")
689  {
690  // Elastic scattering
691  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
692  theElasticProcess->RegisterMe( elastic_lhep0 );
693  pmanager->AddDiscreteProcess( theElasticProcess );
694  // Inelastic scattering
695  G4KaonPlusInelasticProcess* theInelasticProcess = new G4KaonPlusInelasticProcess("inelastic");
696  theInelasticProcess->AddDataSet( G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsKaonPlusInelasticXS::Default_Name()));
697  theInelasticProcess->RegisterMe( theFTFModel1 );
698  theInelasticProcess->RegisterMe( theBERTModel0 );
699  pmanager->AddDiscreteProcess( theInelasticProcess );
700  }
701 
702  else if (particleName == "kaon0S")
703  {
704  // Elastic scattering
705  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
706  theElasticProcess->RegisterMe( elastic_lhep0 );
707  pmanager->AddDiscreteProcess( theElasticProcess );
708  // Inelastic scattering
709  G4KaonZeroSInelasticProcess* theInelasticProcess = new G4KaonZeroSInelasticProcess("inelastic");
710  theInelasticProcess->AddDataSet( G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsKaonZeroInelasticXS::Default_Name()));
711  theInelasticProcess->RegisterMe( theFTFModel1 );
712  theInelasticProcess->RegisterMe( theBERTModel0 );
713  pmanager->AddDiscreteProcess( theInelasticProcess );
714  }
715 
716  else if (particleName == "kaon0L")
717  {
718  // Elastic scattering
719  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
720  theElasticProcess->RegisterMe( elastic_lhep0 );
721  pmanager->AddDiscreteProcess( theElasticProcess );
722  // Inelastic scattering
723  G4KaonZeroLInelasticProcess* theInelasticProcess = new G4KaonZeroLInelasticProcess("inelastic");
724  theInelasticProcess->AddDataSet( G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsKaonZeroInelasticXS::Default_Name()));
725  theInelasticProcess->RegisterMe( theFTFModel1 );
726  theInelasticProcess->RegisterMe( theBERTModel0 );
727  pmanager->AddDiscreteProcess( theInelasticProcess );
728  }
729 
730  else if (particleName == "kaon-")
731  {
732  // Elastic scattering
733  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
734  theElasticProcess->RegisterMe( elastic_lhep0 );
735  pmanager->AddDiscreteProcess( theElasticProcess );
736  // Inelastic scattering
737  G4KaonMinusInelasticProcess* theInelasticProcess = new G4KaonMinusInelasticProcess("inelastic");
738  theInelasticProcess->AddDataSet( G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsKaonMinusInelasticXS::Default_Name()));
739  theInelasticProcess->RegisterMe( theFTFModel1 );
740  theInelasticProcess->RegisterMe( theBERTModel0 );
741  pmanager->AddDiscreteProcess( theInelasticProcess );
742  }
743 
744  else if (particleName == "proton")
745  {
746  // Elastic scattering
747  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
748  theElasticProcess->AddDataSet(G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsProtonElasticXS::Default_Name()));
749  theElasticProcess->RegisterMe( elastic_chip );
750  pmanager->AddDiscreteProcess( theElasticProcess );
751  // Inelastic scattering
752  G4ProtonInelasticProcess* theInelasticProcess = new G4ProtonInelasticProcess("inelastic");
753  theInelasticProcess->AddDataSet( new G4BGGNucleonInelasticXS( G4Proton::Proton() ) );
754  theInelasticProcess->RegisterMe( theFTFModel1 );
755  theInelasticProcess->RegisterMe( theBERTModel0 );
756  pmanager->AddDiscreteProcess( theInelasticProcess );
757  }
758 
759  else if (particleName == "anti_proton")
760  {
761  // Elastic scattering
762  const G4double elastic_elimitAntiNuc = 100.0*CLHEP::MeV;
763  G4AntiNuclElastic* elastic_anuc = new G4AntiNuclElastic();
764  elastic_anuc->SetMinEnergy( elastic_elimitAntiNuc );
765  G4CrossSectionElastic* elastic_anucxs = new G4CrossSectionElastic( elastic_anuc->GetComponentCrossSection() );
766  G4HadronElastic* elastic_lhep2 = new G4HadronElastic();
767  elastic_lhep2->SetMaxEnergy( elastic_elimitAntiNuc );
768  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
769  theElasticProcess->AddDataSet( elastic_anucxs );
770  theElasticProcess->RegisterMe( elastic_lhep2 );
771  theElasticProcess->RegisterMe( elastic_anuc );
772  pmanager->AddDiscreteProcess( theElasticProcess );
773  // Inelastic scattering
774  G4AntiProtonInelasticProcess* theInelasticProcess = new G4AntiProtonInelasticProcess("inelastic");
775  theInelasticProcess->AddDataSet( theAntiNucleonData );
776  theInelasticProcess->RegisterMe( theFTFModel0 );
777  pmanager->AddDiscreteProcess( theInelasticProcess );
778  }
779 
780  else if (particleName == "neutron") {
781  // elastic scattering
782  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
783  theElasticProcess->AddDataSet(G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsNeutronElasticXS::Default_Name()));
784  G4HadronElastic* elastic_neutronChipsModel = new G4ChipsElasticModel();
785  elastic_neutronChipsModel->SetMinEnergy( 19.0*CLHEP::MeV );
786  theElasticProcess->RegisterMe( elastic_neutronChipsModel );
787  G4ParticleHPElastic * theElasticNeutronHP = new G4ParticleHPElastic;
788  theElasticNeutronHP->SetMinEnergy( theHPMin );
789  theElasticNeutronHP->SetMaxEnergy( theHPMax );
790  theElasticProcess->RegisterMe( theElasticNeutronHP );
791  theElasticProcess->AddDataSet( new G4ParticleHPElasticData );
792  pmanager->AddDiscreteProcess( theElasticProcess );
793  // inelastic scattering
794  G4NeutronInelasticProcess* theInelasticProcess = new G4NeutronInelasticProcess("inelastic");
795  theInelasticProcess->AddDataSet( new G4BGGNucleonInelasticXS( G4Neutron::Neutron() ) );
796  theInelasticProcess->RegisterMe( theFTFModel1 );
797  theInelasticProcess->RegisterMe( theBERTModel1 );
798  G4ParticleHPInelastic * theNeutronInelasticHPModel = new G4ParticleHPInelastic;
799  theNeutronInelasticHPModel->SetMinEnergy( theHPMin );
800  theNeutronInelasticHPModel->SetMaxEnergy( theHPMax );
801  theInelasticProcess->RegisterMe( theNeutronInelasticHPModel );
802  theInelasticProcess->AddDataSet( new G4ParticleHPInelasticData );
803  pmanager->AddDiscreteProcess(theInelasticProcess);
804  // capture
805  G4HadronCaptureProcess* theCaptureProcess = new G4HadronCaptureProcess;
806  G4ParticleHPCapture * theNeutronCaptureHPModel = new G4ParticleHPCapture;
807  theNeutronCaptureHPModel->SetMinEnergy( theHPMin );
808  theNeutronCaptureHPModel->SetMaxEnergy( theHPMax );
809  G4NeutronRadCapture* theNeutronRadCapture = new G4NeutronRadCapture();
810  theNeutronRadCapture->SetMinEnergy(theHPMax*0.99);
811  theCaptureProcess->RegisterMe( theNeutronCaptureHPModel );
812  theCaptureProcess->RegisterMe( theNeutronRadCapture);
813  theCaptureProcess->AddDataSet( new G4ParticleHPCaptureData );
814  theCaptureProcess->AddDataSet((G4NeutronCaptureXS*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4NeutronCaptureXS::Default_Name()));
815  pmanager->AddDiscreteProcess(theCaptureProcess);
816  }
817  else if (particleName == "anti_neutron")
818  {
819  // Elastic scattering
820  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
821  theElasticProcess->RegisterMe( elastic_lhep0 );
822  pmanager->AddDiscreteProcess( theElasticProcess );
823  // Inelastic scattering
824  G4AntiNeutronInelasticProcess* theInelasticProcess = new G4AntiNeutronInelasticProcess("inelastic");
825  theInelasticProcess->AddDataSet( theAntiNucleonData );
826  theInelasticProcess->RegisterMe( theFTFModel0 );
827  pmanager->AddDiscreteProcess( theInelasticProcess );
828  }
829 
830  else if (particleName == "deuteron")
831  {
832  // Elastic scattering
833  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
834  theElasticProcess->RegisterMe( elastic_lhep0 );
835  pmanager->AddDiscreteProcess( theElasticProcess );
836  // Inelastic scattering
837  G4DeuteronInelasticProcess* theInelasticProcess = new G4DeuteronInelasticProcess("inelastic");
838  theInelasticProcess->AddDataSet( theGGNuclNuclData );
839  theInelasticProcess->RegisterMe( theFTFModel1 );
840  theInelasticProcess->RegisterMe( theIonBC );
841  pmanager->AddDiscreteProcess( theInelasticProcess );
842  }
843 
844  else if (particleName == "triton")
845  {
846  // Elastic scattering
847  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
848  theElasticProcess->RegisterMe( elastic_lhep0 );
849  pmanager->AddDiscreteProcess( theElasticProcess );
850  // Inelastic scattering
851  G4TritonInelasticProcess* theInelasticProcess = new G4TritonInelasticProcess("inelastic");
852  theInelasticProcess->AddDataSet( theGGNuclNuclData );
853  theInelasticProcess->RegisterMe( theFTFModel1 );
854  theInelasticProcess->RegisterMe( theIonBC );
855  pmanager->AddDiscreteProcess( theInelasticProcess );
856  }
857 
858  else if (particleName == "alpha")
859  {
860  // Elastic scattering
861  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
862  theElasticProcess->RegisterMe( elastic_lhep0 );
863  pmanager->AddDiscreteProcess( theElasticProcess );
864  // Inelastic scattering
865  G4AlphaInelasticProcess* theInelasticProcess = new G4AlphaInelasticProcess("inelastic");
866  theInelasticProcess->AddDataSet( theGGNuclNuclData );
867  theInelasticProcess->RegisterMe( theFTFModel1 );
868  theInelasticProcess->RegisterMe( theIonBC );
869  pmanager->AddDiscreteProcess( theInelasticProcess );
870  }
871  } // while ((*(theParticleIterator))())
872 
873  // Add stopping processes with builder
874  stoppingPhysics->ConstructProcess();
875 }
876 
877 
878 // Decays ///////////////////////////////////////////////////////////////////
879 #include "G4Decay.hh"
880 #include "G4RadioactiveDecay.hh"
881 #include "G4IonTable.hh"
882 #include "G4Ions.hh"
883 
884  template<class T> void TLBE<T>::ConstructGeneral() {
885 
886  // Add Decay Process
887  G4Decay* theDecayProcess = new G4Decay();
888  G4bool theDecayProcessNeverUsed = true; //Check if theDecayProcess will be used
889  theParticleIterator->reset();
890  while( (*(theParticleIterator))() )
891  {
892  G4ParticleDefinition* particle = theParticleIterator->value();
893  G4ProcessManager* pmanager = particle->GetProcessManager();
894 
895  if (theDecayProcess->IsApplicable(*particle) && !particle->IsShortLived())
896  {
897  theDecayProcessNeverUsed = false;
898  pmanager ->AddProcess(theDecayProcess);
899  // set ordering for PostStepDoIt and AtRestDoIt
900  pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
901  pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
902  }
903  }
904 
905  // Declare radioactive decay to the GenericIon in the IonTable.
906  const G4IonTable *theIonTable =
907  G4ParticleTable::GetParticleTable()->GetIonTable();
908  G4RadioactiveDecay *theRadioactiveDecay = new G4RadioactiveDecay();
909 
910  for (G4int i=0; i<theIonTable->Entries(); i++)
911  {
912  G4String particleName = theIonTable->GetParticle(i)->GetParticleName();
913  G4String particleType = theIonTable->GetParticle(i)->GetParticleType();
914 
915  if (particleName == "GenericIon")
916  {
917  G4ProcessManager* pmanager =
918  theIonTable->GetParticle(i)->GetProcessManager();
919  pmanager->SetVerboseLevel(VerboseLevel);
920  pmanager ->AddProcess(theRadioactiveDecay);
921  pmanager ->SetProcessOrdering(theRadioactiveDecay, idxPostStep);
922  pmanager ->SetProcessOrdering(theRadioactiveDecay, idxAtRest);
923  }
924  }
925  //If we actually never used the process, delete it
926  //From Coverity report
927  if ( theDecayProcessNeverUsed ) delete theDecayProcess;
928 }
929 
930 // Cuts /////////////////////////////////////////////////////////////////////
931 template<class T> void TLBE<T>::SetCuts()
932 {
933 
934  if (this->verboseLevel >1)
935  G4cout << "LBE::SetCuts:";
936 
937  if (this->verboseLevel>0){
938  G4cout << "LBE::SetCuts:";
939  G4cout << "CutLength : "
940  << G4BestUnit(this->defaultCutValue,"Length") << G4endl;
941  }
942 
943  //special for low energy physics
944  G4double lowlimit=250*CLHEP::eV;
945  G4ProductionCutsTable * aPCTable = G4ProductionCutsTable::GetProductionCutsTable();
946  aPCTable->SetEnergyRange(lowlimit,100*CLHEP::GeV);
947 
948  // set cut values for gamma at first and for e- second and next for e+,
949  // because some processes for e+/e- need cut values for gamma
950  this->SetCutValue(cutForGamma, "gamma");
951  this->SetCutValue(cutForElectron, "e-");
952  this->SetCutValue(cutForPositron, "e+");
953 
954  // this->SetCutValue(cutForProton, "proton");
955  // this->SetCutValue(cutForProton, "anti_proton");
956  // this->SetCutValue(cutForAlpha, "alpha");
957  // this->SetCutValue(cutForGenericIon, "GenericIon");
958 
959  // this->SetCutValueForOthers(this->defaultCutValue);
960 
961  if (this->verboseLevel>0) this->DumpCutValuesTable();
962 }
963