| Geant4 Cross Reference |
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 // GEANT 4 - Underground Dark Matter Detector Advanced Example
29 //
30 // For information related to this code contact: Alex Howard
31 // e-mail: alexander.howard@cern.ch
32 // --------------------------------------------------------------
33 // Comments
34 //
35 // Underground Advanced
36 // by A. Howard and H. Araujo
37 // (27th November 2001)
38 //
39 // PhysicsList program
40 //
41 // Modified:
42 //
43 // 14-02-03 Fix bugs in msc and hIon instanciation + cut per region
44 //
45 // 05-02-05 AH - changes to G4Decay - added is not short lived protection
46 // and redefined particles to allow non-static creation
47 // i.e. changed construction to G4MesonConstructor, G4BaryonConstructor
48 //
49 // 23-10-09 LP - migrated EM physics from the LowEnergy processes (not supported) to
50 // the new G4Livermore model implementation. Results unchanged.
51 //
52 // --------------------------------------------------------------
53
54 #include <iomanip>
55
56 #include "DMXPhysicsList.hh"
57
58 #include "globals.hh"
59 #include "G4SystemOfUnits.hh"
60 #include "G4ProcessManager.hh"
61 #include "G4ProcessVector.hh"
62
63 #include "G4ParticleDefinition.hh"
64 #include "G4ParticleWithCuts.hh"
65 #include "G4ParticleTypes.hh"
66 #include "G4ParticleTable.hh"
67
68 #include "G4ios.hh"
69 #include "G4UserLimits.hh"
70
71 #include "G4MesonConstructor.hh"
72 #include "G4BaryonConstructor.hh"
73 #include "G4IonConstructor.hh"
74 #include "G4ShortLivedConstructor.hh"
75
76 #include "DMXMaxTimeCuts.hh"
77 #include "DMXMinEkineCuts.hh"
78 #include "G4StepLimiter.hh"
79
80 // gamma
81 #include "G4PhotoElectricEffect.hh"
82 #include "G4LivermorePhotoElectricModel.hh"
83
84 #include "G4ComptonScattering.hh"
85 #include "G4LivermoreComptonModel.hh"
86
87 #include "G4GammaConversion.hh"
88 #include "G4BetheHeitler5DModel.hh"
89
90 #include "G4RayleighScattering.hh"
91 #include "G4LivermoreRayleighModel.hh"
92
93 // e-
94 #include "G4eMultipleScattering.hh"
95
96 #include "G4eIonisation.hh"
97 #include "G4LivermoreIonisationModel.hh"
98
99 #include "G4eBremsstrahlung.hh"
100 #include "G4UniversalFluctuation.hh"
101
102 // e+
103 #include "G4eIonisation.hh"
104 #include "G4eBremsstrahlung.hh"
105 #include "G4eplusAnnihilation.hh"
106
107 // alpha and GenericIon and deuterons, triton, He3:
108 //muon:
109 #include "G4MuIonisation.hh"
110 #include "G4MuBremsstrahlung.hh"
111 #include "G4MuPairProduction.hh"
112 #include "G4MuMultipleScattering.hh"
113 #include "G4MuonMinusCapture.hh"
114
115 //OTHERS:
116 #include "G4hIonisation.hh"
117 #include "G4hMultipleScattering.hh"
118 #include "G4hBremsstrahlung.hh"
119 #include "G4ionIonisation.hh"
120 #include "G4IonParametrisedLossModel.hh"
121 #include "G4NuclearStopping.hh"
122
123 //em process options to allow msc step-limitation to be switched off
124 #include "G4EmParameters.hh"
125 #include "G4VAtomDeexcitation.hh"
126 #include "G4UAtomicDeexcitation.hh"
127 #include "G4LossTableManager.hh"
128
129 #include "G4Scintillation.hh"
130 #include "G4OpAbsorption.hh"
131 #include "G4OpBoundaryProcess.hh"
132 #include "G4OpticalParameters.hh"
133
134 // Elastic processes:
135 #include "G4HadronElasticProcess.hh"
136 #include "G4ChipsElasticModel.hh"
137 #include "G4ElasticHadrNucleusHE.hh"
138
139 // Inelastic processes:
140 #include "G4HadronInelasticProcess.hh"
141
142 // High energy FTFP model and Bertini cascade
143 #include "G4FTFModel.hh"
144 #include "G4LundStringFragmentation.hh"
145 #include "G4ExcitedStringDecay.hh"
146 #include "G4PreCompoundModel.hh"
147 #include "G4GeneratorPrecompoundInterface.hh"
148 #include "G4TheoFSGenerator.hh"
149 #include "G4CascadeInterface.hh"
150
151 // Cross sections
152 #include "G4VCrossSectionDataSet.hh"
153 #include "G4CrossSectionDataSetRegistry.hh"
154
155 #include "G4CrossSectionElastic.hh"
156 #include "G4CrossSectionInelastic.hh"
157 #include "G4BGGPionElasticXS.hh"
158 #include "G4BGGPionInelasticXS.hh"
159 #include "G4AntiNuclElastic.hh"
160
161 #include "G4CrossSectionInelastic.hh"
162 #include "G4BGGNucleonInelasticXS.hh"
163 #include "G4BGGNucleonElasticXS.hh"
164 #include "G4NeutronInelasticXS.hh"
165 #include "G4NeutronElasticXS.hh"
166 #include "G4ComponentAntiNuclNuclearXS.hh"
167 #include "G4ComponentGGNuclNuclXsc.hh"
168 #include "G4ComponentGGHadronNucleusXsc.hh"
169
170 #include "G4HadronElastic.hh"
171 #include "G4NeutronCaptureProcess.hh"
172
173 // Neutron high-precision models: <20 MeV
174 #include "G4ParticleHPElastic.hh"
175 #include "G4ParticleHPElasticData.hh"
176 #include "G4ParticleHPCapture.hh"
177 #include "G4ParticleHPCaptureData.hh"
178 #include "G4ParticleHPInelastic.hh"
179 #include "G4ParticleHPInelasticData.hh"
180
181 // Stopping processes
182 #include "G4HadronStoppingProcess.hh"
183 #include "G4HadronicAbsorptionBertini.hh"
184 #include "G4HadronicAbsorptionFritiof.hh"
185
186 #include "G4HadronicParameters.hh"
187
188 #include "G4Decay.hh"
189 #include "G4RadioactiveDecay.hh"
190 #include "G4PhysicsListHelper.hh"
191 #include "G4NuclideTable.hh"
192 #include "G4NuclearLevelData.hh"
193
194 // Constructor /////////////////////////////////////////////////////////////
195 DMXPhysicsList::DMXPhysicsList() : G4VUserPhysicsList()
196 {
197
198 defaultCutValue = 1.0*micrometer; //
199 cutForGamma = defaultCutValue;
200 cutForElectron = 1.0*nanometer;
201 cutForPositron = defaultCutValue;
202
203 VerboseLevel = 1;
204 OpVerbLevel = 0;
205
206 //set a finer grid of the physic tables in order to improve precision
207 //former LowEnergy models have 200 bins up to 100 GeV
208 G4EmParameters* param = G4EmParameters::Instance();
209 param->SetMaxEnergy(100*GeV);
210 param->SetNumberOfBinsPerDecade(20);
211 param->SetMscStepLimitType(fMinimal);
212 param->SetFluo(true);
213 param->SetPixe(true);
214 param->SetAuger(true);
215
216 G4EmParameters::Instance()->AddPhysics("World","G4RadioactiveDecay");
217 G4DeexPrecoParameters* deex = G4NuclearLevelData::GetInstance()->GetParameters();
218 deex->SetStoreICLevelData(true);
219 deex->SetMaxLifeTime(G4NuclideTable::GetInstance()->GetThresholdOfHalfLife()
220 /std::log(2.));
221 SetVerboseLevel(VerboseLevel);
222 }
223
224 // Destructor //////////////////////////////////////////////////////////////
225 DMXPhysicsList::~DMXPhysicsList()
226 {;}
227
228 // Construct Particles /////////////////////////////////////////////////////
229 void DMXPhysicsList::ConstructParticle()
230 {
231 // In this method, static member functions should be called
232 // for all particles which you want to use.
233 // This ensures that objects of these particle types will be
234 // created in the program.
235
236 ConstructMyBosons();
237 ConstructMyLeptons();
238 ConstructMyHadrons();
239 ConstructMyShortLiveds();
240 }
241
242 // construct Bosons://///////////////////////////////////////////////////
243 void DMXPhysicsList::ConstructMyBosons()
244 {
245 // pseudo-particles
246 G4Geantino::GeantinoDefinition();
247 G4ChargedGeantino::ChargedGeantinoDefinition();
248
249 // gamma
250 G4Gamma::GammaDefinition();
251
252 //OpticalPhotons
253 G4OpticalPhoton::OpticalPhotonDefinition();
254
255 }
256
257 // construct Leptons://///////////////////////////////////////////////////
258 void DMXPhysicsList::ConstructMyLeptons()
259 {
260 // leptons
261 G4Electron::ElectronDefinition();
262 G4Positron::PositronDefinition();
263 G4MuonPlus::MuonPlusDefinition();
264 G4MuonMinus::MuonMinusDefinition();
265
266 G4NeutrinoE::NeutrinoEDefinition();
267 G4AntiNeutrinoE::AntiNeutrinoEDefinition();
268 G4NeutrinoMu::NeutrinoMuDefinition();
269 G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
270 }
271
272 // construct Hadrons://///////////////////////////////////////////////////
273 void DMXPhysicsList::ConstructMyHadrons()
274 {
275 // mesons
276 G4MesonConstructor mConstructor;
277 mConstructor.ConstructParticle();
278
279 // baryons
280 G4BaryonConstructor bConstructor;
281 bConstructor.ConstructParticle();
282
283 // ions
284 G4IonConstructor iConstructor;
285 iConstructor.ConstructParticle();
286 }
287
288 // construct Shortliveds://///////////////////////////////////////////////////
289 void DMXPhysicsList::ConstructMyShortLiveds()
290 {
291 G4ShortLivedConstructor slConstructor;
292 slConstructor.ConstructParticle();
293 }
294
295 // Construct Processes //////////////////////////////////////////////////////
296 void DMXPhysicsList::ConstructProcess()
297 {
298 AddTransportation();
299
300 ConstructEM();
301
302 ConstructOp();
303
304 ConstructHad();
305
306 ConstructGeneral();
307 }
308
309 // Transportation ///////////////////////////////////////////////////////////
310 void DMXPhysicsList::AddTransportation() {
311
312 G4VUserPhysicsList::AddTransportation();
313
314 auto particleIterator=GetParticleIterator();
315 particleIterator->reset();
316 while( (*particleIterator)() ){
317 G4ParticleDefinition* particle = particleIterator->value();
318 G4ProcessManager* pmanager = particle->GetProcessManager();
319 G4String particleName = particle->GetParticleName();
320 // time cuts for ONLY neutrons:
321 if(particleName == "neutron")
322 pmanager->AddDiscreteProcess(new DMXMaxTimeCuts());
323 // Energy cuts to kill charged (embedded in method) particles:
324 pmanager->AddDiscreteProcess(new DMXMinEkineCuts());
325
326 // Step limit applied to all particles:
327 pmanager->AddDiscreteProcess(new G4StepLimiter);
328 }
329 }
330
331 // Electromagnetic Processes ////////////////////////////////////////////////
332 // all charged particles
333 void DMXPhysicsList::ConstructEM() {
334
335 G4LossTableManager* man = G4LossTableManager::Instance();
336 man->SetAtomDeexcitation(new G4UAtomicDeexcitation());
337
338 G4EmParameters* em_params = G4EmParameters::Instance();
339
340 auto particleIterator=GetParticleIterator();
341 particleIterator->reset();
342 while( (*particleIterator)() ){
343 G4ParticleDefinition* particle = particleIterator->value();
344 G4ProcessManager* pmanager = particle->GetProcessManager();
345 G4String particleName = particle->GetParticleName();
346 G4String particleType = particle->GetParticleType();
347 G4double charge = particle->GetPDGCharge();
348
349 if (particleName == "gamma")
350 {
351 //gamma
352 G4RayleighScattering* theRayleigh = new G4RayleighScattering();
353 pmanager->AddDiscreteProcess(theRayleigh);
354
355 G4PhotoElectricEffect* thePhotoElectricEffect = new G4PhotoElectricEffect();
356 thePhotoElectricEffect->SetEmModel(new G4LivermorePhotoElectricModel());
357 pmanager->AddDiscreteProcess(thePhotoElectricEffect);
358
359 G4ComptonScattering* theComptonScattering = new G4ComptonScattering();
360 theComptonScattering->SetEmModel(new G4LivermoreComptonModel());
361 pmanager->AddDiscreteProcess(theComptonScattering);
362
363 G4GammaConversion* theGammaConversion = new G4GammaConversion();
364 theGammaConversion->SetEmModel(new G4BetheHeitler5DModel());
365 pmanager->AddDiscreteProcess(theGammaConversion);
366
367 }
368 else if (particleName == "e-")
369 {
370 //electron
371 // process ordering: AddProcess(name, at rest, along step, post step)
372 // Multiple scattering
373 G4eMultipleScattering* msc = new G4eMultipleScattering();
374 em_params->SetMscStepLimitType(fUseDistanceToBoundary);
375 pmanager->AddProcess(msc,-1, 1, -1);
376
377 // Ionisation
378 G4eIonisation* eIonisation = new G4eIonisation();
379 G4VEmModel* theIoniLiv = new G4LivermoreIonisationModel();
380 theIoniLiv->SetHighEnergyLimit(0.1*MeV);
381 eIonisation->AddEmModel(0, theIoniLiv, new G4UniversalFluctuation() );
382 em_params->SetStepFunction(0.2, 100*um); //improved precision in tracking
383 pmanager->AddProcess(eIonisation,-1, 2, 1);
384
385 // Bremsstrahlung
386 G4eBremsstrahlung* eBremsstrahlung = new G4eBremsstrahlung();
387 pmanager->AddProcess(eBremsstrahlung, -1,-3, 2);
388 }
389 else if (particleName == "e+")
390 {
391 //positron
392 G4eMultipleScattering* msc = new G4eMultipleScattering();
393 msc->SetStepLimitType(fUseDistanceToBoundary);
394 pmanager->AddProcess(msc,-1, 1, -1);
395
396 // Ionisation
397 G4eIonisation* eIonisation = new G4eIonisation();
398 // eIonisation->SetStepFunction(0.2, 100*um); //
399 pmanager->AddProcess(eIonisation, -1, 2, 1);
400
401 //Bremsstrahlung (use default, no low-energy available)
402 pmanager->AddProcess(new G4eBremsstrahlung(), -1,-1, 2);
403
404 //Annihilation
405 pmanager->AddProcess(new G4eplusAnnihilation(),0,-1, 3);
406 }
407 else if( particleName == "mu+" ||
408 particleName == "mu-" )
409 {
410 //muon
411 pmanager->AddProcess(new G4MuMultipleScattering, -1, 1,-1);
412 pmanager->AddProcess(new G4MuIonisation(), -1, 2, 1);
413 pmanager->AddProcess(new G4MuBremsstrahlung(), -1,-1, 2);
414 pmanager->AddProcess(new G4MuPairProduction(), -1,-1, 3);
415 if( particleName == "mu-" )
416 pmanager->AddProcess(new G4MuonMinusCapture(), 0,-1,-1);
417 }
418 else if (particleName == "proton" ||
419 particleName == "pi+" ||
420 particleName == "pi-")
421 {
422 //multiple scattering
423 pmanager->AddProcess(new G4hMultipleScattering, -1, 1, -1);
424
425 //ionisation
426 G4hIonisation* hIonisation = new G4hIonisation();
427 em_params->SetStepFunctionMuHad(0.2, 50*um);
428 pmanager->AddProcess(hIonisation, -1, 2, 1);
429
430 //bremmstrahlung
431 pmanager->AddProcess(new G4hBremsstrahlung, -1,-3, 2);
432 }
433 else if(particleName == "alpha" ||
434 particleName == "deuteron" ||
435 particleName == "triton" ||
436 particleName == "He3")
437 {
438 //multiple scattering
439 pmanager->AddProcess(new G4hMultipleScattering,-1,1,-1);
440
441 //ionisation
442 G4ionIonisation* ionIoni = new G4ionIonisation();
443 em_params->SetStepFunctionLightIons(0.1, 1*CLHEP::um);
444 pmanager->AddProcess(ionIoni, -1, 2, 1);
445 pmanager->AddProcess(new G4NuclearStopping(), -1, 3,-1);
446 }
447 else if (particleName == "GenericIon")
448 {
449 // OBJECT may be dynamically created as either a GenericIon or nucleus
450 // G4Nucleus exists and therefore has particle type nucleus
451 // genericIon:
452
453 //multiple scattering
454 pmanager->AddProcess(new G4hMultipleScattering,-1,1,-1);
455
456 //ionisation
457 G4ionIonisation* ionIoni = new G4ionIonisation();
458 em_params->SetStepFunctionIons(0.1, 1*CLHEP::um);
459 pmanager->AddProcess(ionIoni, -1, 2, 1);
460 pmanager->AddProcess(new G4NuclearStopping(), -1, 3,-1);
461 }
462
463 else if ((!particle->IsShortLived()) &&
464 (charge != 0.0) &&
465 (particle->GetParticleName() != "chargedgeantino"))
466 {
467 //all others charged particles except geantino
468 G4hMultipleScattering* aMultipleScattering = new G4hMultipleScattering();
469 G4hIonisation* ahadronIon = new G4hIonisation();
470
471 //multiple scattering
472 pmanager->AddProcess(aMultipleScattering,-1,1,-1);
473
474 //ionisation
475 pmanager->AddProcess(ahadronIon, -1,2,1);
476 }
477 }
478 }
479
480 // Optical Processes ////////////////////////////////////////////////////////
481 void DMXPhysicsList::ConstructOp()
482 {
483 G4OpticalParameters* opParams = G4OpticalParameters::Instance();
484 G4Scintillation* theScintProcessDef = new G4Scintillation("Scintillation");
485 opParams->SetScintTrackSecondariesFirst(true);
486 opParams->SetScintByParticleType(true);
487
488 // optical processes
489 G4OpAbsorption* theAbsorptionProcess = new G4OpAbsorption();
490 G4OpBoundaryProcess* theBoundaryProcess = new G4OpBoundaryProcess();
491
492 auto particleIterator=GetParticleIterator();
493 particleIterator->reset();
494 while( (*particleIterator)() )
495 {
496 G4ParticleDefinition* particle = particleIterator->value();
497 G4ProcessManager* pmanager = particle->GetProcessManager();
498 G4String particleName = particle->GetParticleName();
499 if (theScintProcessDef->IsApplicable(*particle)) {
500 pmanager->AddProcess(theScintProcessDef);
501 pmanager->SetProcessOrderingToLast(theScintProcessDef,idxAtRest);
502 pmanager->SetProcessOrderingToLast(theScintProcessDef,idxPostStep);
503 }
504
505 if (particleName == "opticalphoton") {
506 pmanager->AddDiscreteProcess(theAbsorptionProcess);
507 pmanager->AddDiscreteProcess(theBoundaryProcess);
508 }
509 }
510 }
511
512 // Hadronic processes ////////////////////////////////////////////////////////
513
514 void DMXPhysicsList::ConstructHad()
515 {
516 //Elastic models
517 G4HadronElastic* elastic_lhep0 = new G4HadronElastic();
518 G4ChipsElasticModel* elastic_chip = new G4ChipsElasticModel();
519 G4ElasticHadrNucleusHE* elastic_he = new G4ElasticHadrNucleusHE();
520
521 // Inelastic scattering
522 const G4double theFTFMin0 = 0.0*GeV;
523 const G4double theFTFMin1 = 3.0*GeV;
524 const G4double theFTFMax = G4HadronicParameters::Instance()->GetMaxEnergy();
525 const G4double theBERTMin0 = 0.0*GeV;
526 const G4double theBERTMin1 = 19.0*MeV;
527 const G4double theBERTMax = 6.0*GeV;
528 const G4double theHPMin = 0.0*GeV;
529 const G4double theHPMax = 20.0*MeV;
530
531 G4FTFModel * theStringModel = new G4FTFModel;
532 G4ExcitedStringDecay * theStringDecay = new G4ExcitedStringDecay( new G4LundStringFragmentation );
533 theStringModel->SetFragmentationModel( theStringDecay );
534 G4PreCompoundModel * thePreEquilib = new G4PreCompoundModel( new G4ExcitationHandler );
535 G4GeneratorPrecompoundInterface * theCascade = new G4GeneratorPrecompoundInterface( thePreEquilib );
536
537 G4TheoFSGenerator * theFTFModel0 = new G4TheoFSGenerator( "FTFP" );
538 theFTFModel0->SetHighEnergyGenerator( theStringModel );
539 theFTFModel0->SetTransport( theCascade );
540 theFTFModel0->SetMinEnergy( theFTFMin0 );
541 theFTFModel0->SetMaxEnergy( theFTFMax );
542
543 G4TheoFSGenerator * theFTFModel1 = new G4TheoFSGenerator( "FTFP" );
544 theFTFModel1->SetHighEnergyGenerator( theStringModel );
545 theFTFModel1->SetTransport( theCascade );
546 theFTFModel1->SetMinEnergy( theFTFMin1 );
547 theFTFModel1->SetMaxEnergy( theFTFMax );
548
549 G4CascadeInterface * theBERTModel0 = new G4CascadeInterface;
550 theBERTModel0->SetMinEnergy( theBERTMin0 );
551 theBERTModel0->SetMaxEnergy( theBERTMax );
552
553 G4CascadeInterface * theBERTModel1 = new G4CascadeInterface;
554 theBERTModel1->SetMinEnergy( theBERTMin1 );
555 theBERTModel1->SetMaxEnergy( theBERTMax );
556
557 G4VCrossSectionDataSet * theAntiNucleonData = new G4CrossSectionInelastic( new G4ComponentAntiNuclNuclearXS );
558 G4ComponentGGNuclNuclXsc * ggNuclNuclXsec = new G4ComponentGGNuclNuclXsc();
559 G4VCrossSectionDataSet * theGGNuclNuclData = new G4CrossSectionInelastic(ggNuclNuclXsec);
560 G4VCrossSectionDataSet * theGGNNEl = new G4CrossSectionElastic(ggNuclNuclXsec);
561 G4ComponentGGHadronNucleusXsc * ggHNXsec = new G4ComponentGGHadronNucleusXsc();
562 G4VCrossSectionDataSet * theGGHNEl = new G4CrossSectionElastic(ggHNXsec);
563 G4VCrossSectionDataSet * theGGHNInel = new G4CrossSectionInelastic(ggHNXsec);
564
565 auto particleIterator=GetParticleIterator();
566 particleIterator->reset();
567 while ((*particleIterator)())
568 {
569 G4ParticleDefinition* particle = particleIterator->value();
570 G4ProcessManager* pmanager = particle->GetProcessManager();
571 G4String particleName = particle->GetParticleName();
572
573 if (particleName == "pi+")
574 {
575 // Elastic scattering
576 G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
577 theElasticProcess->AddDataSet( new G4BGGPionElasticXS( particle ) );
578 theElasticProcess->RegisterMe( elastic_he );
579 pmanager->AddDiscreteProcess( theElasticProcess );
580 //Inelastic scattering
581 G4HadronInelasticProcess* theInelasticProcess =
582 new G4HadronInelasticProcess( "inelastic", G4PionPlus::Definition() );
583 theInelasticProcess->AddDataSet( new G4BGGPionElasticXS( particle ) );
584 theInelasticProcess->RegisterMe( theFTFModel1 );
585 theInelasticProcess->RegisterMe( theBERTModel0 );
586 pmanager->AddDiscreteProcess( theInelasticProcess );
587 }
588
589 else if (particleName == "pi-")
590 {
591 // Elastic scattering
592 G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
593 theElasticProcess->AddDataSet( new G4BGGPionElasticXS( particle ) );
594 theElasticProcess->RegisterMe( elastic_he );
595 pmanager->AddDiscreteProcess( theElasticProcess );
596 //Inelastic scattering
597 G4HadronInelasticProcess* theInelasticProcess =
598 new G4HadronInelasticProcess( "inelastic", G4PionMinus::Definition() );
599 theInelasticProcess->AddDataSet( new G4BGGPionInelasticXS( particle ) );
600 theInelasticProcess->RegisterMe( theFTFModel1 );
601 theInelasticProcess->RegisterMe( theBERTModel0 );
602 pmanager->AddDiscreteProcess( theInelasticProcess );
603 //Absorption
604 pmanager->AddRestProcess(new G4HadronicAbsorptionBertini(G4PionMinus::Definition()), ordDefault);
605 }
606 else if (particleName == "kaon+")
607 {
608 // Elastic scattering
609 G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
610 theElasticProcess->AddDataSet( theGGHNEl );
611 theElasticProcess->RegisterMe( elastic_lhep0 );
612 pmanager->AddDiscreteProcess( theElasticProcess );
613 // Inelastic scattering
614 G4HadronInelasticProcess* theInelasticProcess =
615 new G4HadronInelasticProcess( "inelastic", G4KaonPlus::Definition() );
616 theInelasticProcess->AddDataSet( theGGHNInel );
617 theInelasticProcess->RegisterMe( theFTFModel1 );
618 theInelasticProcess->RegisterMe( theBERTModel0 );
619 pmanager->AddDiscreteProcess( theInelasticProcess );
620 }
621 else if (particleName == "kaon0S")
622 {
623 // Elastic scattering
624 G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
625 theElasticProcess->AddDataSet( theGGHNEl );
626 theElasticProcess->RegisterMe( elastic_lhep0 );
627 pmanager->AddDiscreteProcess( theElasticProcess );
628 // Inelastic scattering
629 G4HadronInelasticProcess* theInelasticProcess =
630 new G4HadronInelasticProcess( "inelastic", G4KaonZeroShort::Definition() );
631 theInelasticProcess->AddDataSet( theGGHNInel );
632 theInelasticProcess->RegisterMe( theFTFModel1 );
633 theInelasticProcess->RegisterMe( theBERTModel0 );
634 pmanager->AddDiscreteProcess( theInelasticProcess );
635 }
636
637 else if (particleName == "kaon0L")
638 {
639 // Elastic scattering
640 G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
641 theElasticProcess->AddDataSet( theGGHNEl );
642 theElasticProcess->RegisterMe( elastic_lhep0 );
643 pmanager->AddDiscreteProcess( theElasticProcess );
644 // Inelastic scattering
645 G4HadronInelasticProcess* theInelasticProcess =
646 new G4HadronInelasticProcess( "inelastic", G4KaonZeroLong::Definition() );
647 theInelasticProcess->AddDataSet( theGGHNInel );
648 theInelasticProcess->RegisterMe( theFTFModel1 );
649 theInelasticProcess->RegisterMe( theBERTModel0 );
650 pmanager->AddDiscreteProcess( theInelasticProcess );
651 }
652
653 else if (particleName == "kaon-")
654 {
655 // Elastic scattering
656 G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
657 theElasticProcess->AddDataSet( theGGHNEl );
658 theElasticProcess->RegisterMe( elastic_lhep0 );
659 pmanager->AddDiscreteProcess( theElasticProcess );
660 // Inelastic scattering
661 G4HadronInelasticProcess* theInelasticProcess =
662 new G4HadronInelasticProcess( "inelastic", G4KaonMinus::Definition() );
663 theInelasticProcess->AddDataSet( theGGHNInel );
664 theInelasticProcess->RegisterMe( theFTFModel1 );
665 theInelasticProcess->RegisterMe( theBERTModel0 );
666 pmanager->AddDiscreteProcess( theInelasticProcess );
667 pmanager->AddRestProcess(new G4HadronicAbsorptionBertini(G4KaonMinus::Definition()), ordDefault);
668 }
669
670 else if (particleName == "proton")
671 {
672 // Elastic scattering
673 G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
674 theElasticProcess->AddDataSet( new G4BGGNucleonElasticXS( G4Proton::Proton() ) );
675 theElasticProcess->RegisterMe( elastic_chip );
676 pmanager->AddDiscreteProcess( theElasticProcess );
677 // Inelastic scattering
678 G4HadronInelasticProcess* theInelasticProcess =
679 new G4HadronInelasticProcess( "inelastic", G4Proton::Definition() );
680 theInelasticProcess->AddDataSet( new G4BGGNucleonInelasticXS( G4Proton::Proton() ) );
681 theInelasticProcess->RegisterMe( theFTFModel1 );
682 theInelasticProcess->RegisterMe( theBERTModel0 );
683 pmanager->AddDiscreteProcess( theInelasticProcess );
684 }
685 else if (particleName == "anti_proton")
686 {
687 // Elastic scattering
688 const G4double elastic_elimitAntiNuc = 100.0*MeV;
689 G4AntiNuclElastic* elastic_anuc = new G4AntiNuclElastic();
690 elastic_anuc->SetMinEnergy( elastic_elimitAntiNuc );
691 G4CrossSectionElastic* elastic_anucxs = new G4CrossSectionElastic( elastic_anuc->GetComponentCrossSection() );
692 G4HadronElastic* elastic_lhep2 = new G4HadronElastic();
693 elastic_lhep2->SetMaxEnergy( elastic_elimitAntiNuc );
694 G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
695 theElasticProcess->AddDataSet( elastic_anucxs );
696 theElasticProcess->RegisterMe( elastic_lhep2 );
697 theElasticProcess->RegisterMe( elastic_anuc );
698 pmanager->AddDiscreteProcess( theElasticProcess );
699 // Inelastic scattering
700 G4HadronInelasticProcess* theInelasticProcess =
701 new G4HadronInelasticProcess( "inelastic", G4AntiProton::Definition() );
702 theInelasticProcess->AddDataSet( theAntiNucleonData );
703 theInelasticProcess->RegisterMe( theFTFModel0 );
704 pmanager->AddDiscreteProcess( theInelasticProcess );
705 // Absorption
706 pmanager->AddRestProcess(new G4HadronicAbsorptionFritiof(G4AntiProton::Definition()), ordDefault);
707 }
708 else if (particleName == "neutron") {
709 // elastic scattering
710 G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
711 theElasticProcess->AddDataSet(new G4NeutronElasticXS());
712 G4HadronElastic* elastic_neutronChipsModel = new G4ChipsElasticModel();
713 elastic_neutronChipsModel->SetMinEnergy( 19.0*MeV );
714 theElasticProcess->RegisterMe( elastic_neutronChipsModel );
715 G4ParticleHPElastic * theElasticNeutronHP = new G4ParticleHPElastic;
716 theElasticNeutronHP->SetMinEnergy( theHPMin );
717 theElasticNeutronHP->SetMaxEnergy( theHPMax );
718 theElasticProcess->RegisterMe( theElasticNeutronHP );
719 theElasticProcess->AddDataSet( new G4ParticleHPElasticData );
720 pmanager->AddDiscreteProcess( theElasticProcess );
721 // inelastic scattering
722 G4HadronInelasticProcess* theInelasticProcess =
723 new G4HadronInelasticProcess( "inelastic", G4Neutron::Definition() );
724 theInelasticProcess->AddDataSet( new G4NeutronInelasticXS() );
725 theInelasticProcess->RegisterMe( theFTFModel1 );
726 theInelasticProcess->RegisterMe( theBERTModel1 );
727 G4ParticleHPInelastic * theNeutronInelasticHPModel = new G4ParticleHPInelastic;
728 theNeutronInelasticHPModel->SetMinEnergy( theHPMin );
729 theNeutronInelasticHPModel->SetMaxEnergy( theHPMax );
730 theInelasticProcess->RegisterMe( theNeutronInelasticHPModel );
731 theInelasticProcess->AddDataSet( new G4ParticleHPInelasticData );
732 pmanager->AddDiscreteProcess(theInelasticProcess);
733 // capture
734 G4NeutronCaptureProcess* theCaptureProcess =
735 new G4NeutronCaptureProcess;
736 G4ParticleHPCapture * theLENeutronCaptureModel = new G4ParticleHPCapture;
737 theLENeutronCaptureModel->SetMinEnergy(theHPMin);
738 theLENeutronCaptureModel->SetMaxEnergy(theHPMax);
739 theCaptureProcess->RegisterMe(theLENeutronCaptureModel);
740 theCaptureProcess->AddDataSet( new G4ParticleHPCaptureData);
741 pmanager->AddDiscreteProcess(theCaptureProcess);
742 }
743 else if (particleName == "anti_neutron")
744 {
745 // Elastic scattering
746 G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
747 theElasticProcess->AddDataSet( theGGHNEl );
748 theElasticProcess->RegisterMe( elastic_lhep0 );
749 pmanager->AddDiscreteProcess( theElasticProcess );
750 // Inelastic scattering (include annihilation on-fly)
751 G4HadronInelasticProcess* theInelasticProcess =
752 new G4HadronInelasticProcess( "inelastic", G4AntiNeutron::Definition() );
753 theInelasticProcess->AddDataSet( theAntiNucleonData );
754 theInelasticProcess->RegisterMe( theFTFModel0 );
755 pmanager->AddDiscreteProcess( theInelasticProcess );
756 }
757 else if (particleName == "deuteron")
758 {
759 // Elastic scattering
760 G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
761 theElasticProcess->AddDataSet( theGGNNEl );
762 theElasticProcess->RegisterMe( elastic_lhep0 );
763 pmanager->AddDiscreteProcess( theElasticProcess );
764 // Inelastic scattering
765 G4HadronInelasticProcess* theInelasticProcess =
766 new G4HadronInelasticProcess( "inelastic", G4Deuteron::Definition() );
767 theInelasticProcess->AddDataSet( theGGNuclNuclData );
768 theInelasticProcess->RegisterMe( theFTFModel1 );
769 theInelasticProcess->RegisterMe( theBERTModel0 );
770 pmanager->AddDiscreteProcess( theInelasticProcess );
771 }
772 else if (particleName == "triton")
773 {
774 // Elastic scattering
775 G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
776 theElasticProcess->AddDataSet( theGGNNEl );
777 theElasticProcess->RegisterMe( elastic_lhep0 );
778 pmanager->AddDiscreteProcess( theElasticProcess );
779 // Inelastic scattering
780 G4HadronInelasticProcess* theInelasticProcess =
781 new G4HadronInelasticProcess( "inelastic", G4Triton::Definition() );
782 theInelasticProcess->AddDataSet( theGGNuclNuclData );
783 theInelasticProcess->RegisterMe( theFTFModel1 );
784 theInelasticProcess->RegisterMe( theBERTModel0 );
785 pmanager->AddDiscreteProcess( theInelasticProcess );
786 }
787 else if (particleName == "alpha")
788 {
789 // Elastic scattering
790 G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
791 theElasticProcess->AddDataSet( theGGNNEl );
792 theElasticProcess->RegisterMe( elastic_lhep0 );
793 pmanager->AddDiscreteProcess( theElasticProcess );
794 // Inelastic scattering
795 G4HadronInelasticProcess* theInelasticProcess =
796 new G4HadronInelasticProcess( "inelastic", G4Alpha::Definition() );
797 theInelasticProcess->AddDataSet( theGGNuclNuclData );
798 theInelasticProcess->RegisterMe( theFTFModel1 );
799 theInelasticProcess->RegisterMe( theBERTModel0 );
800 pmanager->AddDiscreteProcess( theInelasticProcess );
801 }
802 }
803 }
804
805 // Decays ///////////////////////////////////////////////////////////////////
806 void DMXPhysicsList::ConstructGeneral() {
807
808 // Add Decay Process
809 G4Decay* theDecayProcess = new G4Decay();
810 auto particleIterator=GetParticleIterator();
811 particleIterator->reset();
812 while( (*particleIterator)() )
813 {
814 G4ParticleDefinition* particle = particleIterator->value();
815 G4ProcessManager* pmanager = particle->GetProcessManager();
816
817 if (theDecayProcess->IsApplicable(*particle) && !particle->IsShortLived())
818 {
819 pmanager ->AddProcess(theDecayProcess);
820 // set ordering for PostStepDoIt and AtRestDoIt
821 pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
822 pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
823 }
824 }
825
826 // Declare radioactive decay to the GenericIon in the IonTable.
827 G4LossTableManager* man = G4LossTableManager::Instance();
828 G4VAtomDeexcitation* ad = man->AtomDeexcitation();
829 if(!ad) {
830 G4EmParameters::Instance()->SetAugerCascade(true);
831 ad = new G4UAtomicDeexcitation();
832 man->SetAtomDeexcitation(ad);
833 ad->InitialiseAtomicDeexcitation();
834 }
835
836 G4PhysicsListHelper::GetPhysicsListHelper()->
837 RegisterProcess(new G4RadioactiveDecay(), G4GenericIon::GenericIon());
838 }
839
840 // Cuts /////////////////////////////////////////////////////////////////////
841 void DMXPhysicsList::SetCuts()
842 {
843
844 if (verboseLevel >1)
845 G4cout << "DMXPhysicsList::SetCuts:";
846
847 if (verboseLevel>0){
848 G4cout << "DMXPhysicsList::SetCuts:";
849 G4cout << "CutLength : "
850 << G4BestUnit(defaultCutValue,"Length") << G4endl;
851 }
852
853 //special for low energy physics
854 G4double lowlimit=250*eV;
855 G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(lowlimit,100.*GeV);
856
857 // set cut values for gamma at first and for e- second and next for e+,
858 // because some processes for e+/e- need cut values for gamma
859 SetCutValue(cutForGamma, "gamma");
860 SetCutValue(cutForElectron, "e-");
861 SetCutValue(cutForPositron, "e+");
862
863 if (verboseLevel>0) DumpCutValuesTable();
864 }
865
866