source: trunk/source/physics_lists/lists/include/LBE.icc

Last change on this file was 1340, checked in by garnier, 14 years ago

update ti head
File size: 31.0 KB
Line 
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//
51// --------------------------------------------------------------
52
53
54#include "globals.hh"
55#include "G4ProcessManager.hh"
56#include "G4ProcessVector.hh"
57
58#include "G4ParticleTypes.hh"
59#include "G4ParticleTable.hh"
60#include "G4ProductionCutsTable.hh"
61
62#include "G4ios.hh"
63#include <iomanip> 
64
65#include "G4UserLimits.hh"
66#include "G4DataQuestionaire.hh"
67#include "G4WarnPLStatus.hh"
68
69
70// Constructor /////////////////////////////////////////////////////////////
71template<class T> TLBE<T>::TLBE() : G4VUserPhysicsList()
72{
73
74  G4DataQuestionaire it(photon, lowenergy, neutron, radioactive);
75  G4cout << "You are using the simulation engine: LBE 5.3"<<G4endl;
76  G4cout <<G4endl<<G4endl;
77  this->defaultCutValue     = 1.0*micrometer; //
78  cutForGamma         = this->defaultCutValue;
79  cutForElectron      = 1.0*nanometer;
80  cutForPositron      = this->defaultCutValue;
81  cutForProton        = this->defaultCutValue;
82  cutForAlpha         = 1.0*nanometer;
83  cutForGenericIon    = 1.0*nanometer;
84
85  VerboseLevel = 1;
86  OpVerbLevel = 0;
87
88  this->SetVerboseLevel(VerboseLevel);
89}
90
91
92// Destructor //////////////////////////////////////////////////////////////
93template<class T> TLBE<T>::~TLBE()
94{;}
95
96
97// Construct Particles /////////////////////////////////////////////////////
98 template<class T> void TLBE<T>::ConstructParticle()
99{
100
101  // In this method, static member functions should be called
102  // for all particles which you want to use.
103  // This ensures that objects of these particle types will be
104  // created in the program.
105
106  ConstructMyBosons();
107  ConstructMyLeptons();
108  ConstructMyMesons();
109  ConstructMyBaryons();
110  ConstructMyIons();
111  ConstructMyShortLiveds();
112
113}
114
115
116// construct Bosons://///////////////////////////////////////////////////
117 template<class T> void TLBE<T>::ConstructMyBosons()
118{
119  // pseudo-particles
120  G4Geantino::GeantinoDefinition();
121  G4ChargedGeantino::ChargedGeantinoDefinition();
122 
123  // gamma
124  G4Gamma::GammaDefinition();
125
126  //OpticalPhotons
127  G4OpticalPhoton::OpticalPhotonDefinition();
128
129}
130
131
132// construct Leptons://///////////////////////////////////////////////////
133 template<class T> void TLBE<T>::ConstructMyLeptons()
134{
135  // leptons
136  G4Electron::ElectronDefinition();
137  G4Positron::PositronDefinition();
138  G4MuonPlus::MuonPlusDefinition();
139  G4MuonMinus::MuonMinusDefinition();
140
141  G4NeutrinoE::NeutrinoEDefinition();
142  G4AntiNeutrinoE::AntiNeutrinoEDefinition();
143  G4NeutrinoMu::NeutrinoMuDefinition();
144  G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
145}
146
147#include "G4MesonConstructor.hh"
148#include "G4BaryonConstructor.hh"
149#include "G4IonConstructor.hh"
150
151
152// construct Mesons://///////////////////////////////////////////////////
153 template<class T> void TLBE<T>::ConstructMyMesons()
154{
155 //  mesons
156  G4MesonConstructor mConstructor;
157  mConstructor.ConstructParticle();
158
159}
160
161
162// construct Baryons://///////////////////////////////////////////////////
163 template<class T> void TLBE<T>::ConstructMyBaryons()
164{
165 //  baryons
166  G4BaryonConstructor bConstructor;
167  bConstructor.ConstructParticle();
168
169}
170
171
172// construct Ions://///////////////////////////////////////////////////
173 template<class T> void TLBE<T>::ConstructMyIons()
174{
175 //  ions
176  G4IonConstructor iConstructor;
177  iConstructor.ConstructParticle();
178
179}
180
181// construct Shortliveds://///////////////////////////////////////////////////
182 template<class T> void TLBE<T>::ConstructMyShortLiveds()
183{
184  // ShortLiveds
185  ;
186}
187
188
189
190
191// Construct Processes //////////////////////////////////////////////////////
192 template<class T> void TLBE<T>::ConstructProcess()
193{
194
195  AddTransportation();
196
197  ConstructEM();
198
199  ConstructOp();
200
201  ConstructHad();
202
203  ConstructGeneral();
204
205}
206
207
208// Transportation ///////////////////////////////////////////////////////////
209#include "MaxTimeCuts.hh"
210#include "MinEkineCuts.hh"
211
212 template<class T> void TLBE<T>::AddTransportation() {
213 
214  G4VUserPhysicsList::AddTransportation();
215 
216  this->theParticleIterator->reset();
217  while( (*(this->theParticleIterator))() ){
218    G4ParticleDefinition* particle = this->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 MaxTimeCuts());
224    // Energy cuts to kill charged (embedded in method) particles:
225    pmanager->AddDiscreteProcess(new MinEkineCuts());
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 "G4UrbanMscModel93.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
273// msc models
274#include "G4UrbanMscModel93.hh"
275
276
277//muon:
278#include "G4MuIonisation.hh"
279#include "G4MuBremsstrahlung.hh"
280#include "G4MuPairProduction.hh"
281#include "G4MuonMinusCaptureAtRest.hh"
282
283//OTHERS:
284//#include "G4hIonisation.hh" // standard hadron ionisation
285
286 template<class T> void TLBE<T>::ConstructEM() {
287
288 // models & processes:
289 // Use Livermore models up to 20 MeV, and standard
290 // models for higher energy
291 G4double LivermoreHighEnergyLimit = 20*MeV;
292 //
293  this->theParticleIterator->reset();
294  while( (*(this->theParticleIterator))() ){
295    G4ParticleDefinition* particle = this->theParticleIterator->value();
296    G4ProcessManager* pmanager = particle->GetProcessManager();
297    G4String particleName = particle->GetParticleName();
298    G4String particleType = particle->GetParticleType();
299    G4double charge = particle->GetPDGCharge();
300   
301    if (particleName == "gamma")
302      {
303      G4PhotoElectricEffect* thePhotoElectricEffect = new G4PhotoElectricEffect();
304      G4LivermorePhotoElectricModel* theLivermorePhotoElectricModel =
305        new G4LivermorePhotoElectricModel();
306      theLivermorePhotoElectricModel->SetHighEnergyLimit(LivermoreHighEnergyLimit);
307      thePhotoElectricEffect->AddEmModel(0, theLivermorePhotoElectricModel);
308      pmanager->AddDiscreteProcess(thePhotoElectricEffect);
309
310      G4ComptonScattering* theComptonScattering = new G4ComptonScattering();
311      G4LivermoreComptonModel* theLivermoreComptonModel =
312        new G4LivermoreComptonModel();
313      theLivermoreComptonModel->SetHighEnergyLimit(LivermoreHighEnergyLimit);
314      theComptonScattering->AddEmModel(0, theLivermoreComptonModel);
315      pmanager->AddDiscreteProcess(theComptonScattering);
316
317      G4GammaConversion* theGammaConversion = new G4GammaConversion();
318      G4LivermoreGammaConversionModel* theLivermoreGammaConversionModel =
319        new G4LivermoreGammaConversionModel();
320      theLivermoreGammaConversionModel->SetHighEnergyLimit(LivermoreHighEnergyLimit);
321      theGammaConversion->AddEmModel(0, theLivermoreGammaConversionModel);
322      pmanager->AddDiscreteProcess(theGammaConversion);
323
324      G4RayleighScattering* theRayleigh = new G4RayleighScattering();
325      G4LivermoreRayleighModel* theRayleighModel = new G4LivermoreRayleighModel();
326      theRayleighModel->SetHighEnergyLimit(LivermoreHighEnergyLimit);
327      theRayleigh->AddEmModel(0, theRayleighModel);
328      pmanager->AddDiscreteProcess(theRayleigh);
329
330      }
331    else if (particleName == "e-")
332      {
333       //electron
334       // process ordering: AddProcess(name, at rest, along step, post step)
335       // -1 = not implemented, then ordering
336        G4eMultipleScattering* msc = new G4eMultipleScattering();     
337        msc->AddEmModel(0, new G4UrbanMscModel93());
338        msc->SetStepLimitType(fUseDistanceToBoundary);
339        pmanager->AddProcess(msc,                   -1, 1, 1);
340     
341       // Ionisation
342       G4eIonisation* eIoni = new G4eIonisation();
343       G4LivermoreIonisationModel* theIoniLivermore = new
344        G4LivermoreIonisationModel();
345       theIoniLivermore->SetHighEnergyLimit(1*MeV);
346       eIoni->AddEmModel(0, theIoniLivermore, new G4UniversalFluctuation() );
347       eIoni->SetStepFunction(0.2, 100*um); //     
348       pmanager->AddProcess(eIoni,                 -1, 2, 2);
349     
350       // Bremsstrahlung
351       G4eBremsstrahlung* eBrem = new G4eBremsstrahlung();
352       G4LivermoreBremsstrahlungModel* theBremLivermore = new
353         G4LivermoreBremsstrahlungModel();
354       theBremLivermore->SetHighEnergyLimit(LivermoreHighEnergyLimit);
355       eBrem->AddEmModel(0, theBremLivermore);
356       pmanager->AddProcess(eBrem, -1,-3, 3);   
357      }
358    else if (particleName == "e+")
359      {
360        //positron
361      G4eMultipleScattering* msc = new G4eMultipleScattering();
362      msc->AddEmModel(0, new G4UrbanMscModel93());     
363      msc->SetStepLimitType(fUseDistanceToBoundary);
364      pmanager->AddProcess(msc,                   -1, 1, 1);
365      G4eIonisation* eIoni = new G4eIonisation();
366      eIoni->SetStepFunction(0.2, 100*um);     
367      pmanager->AddProcess(eIoni,                 -1, 2, 2);
368      pmanager->AddProcess(new G4eBremsstrahlung, -1,-3, 3);     
369      pmanager->AddProcess(new G4eplusAnnihilation,0,-1, 4);
370      }
371    else if( particleName == "mu+" ||
372             particleName == "mu-"    )
373      {
374        //muon 
375        G4MuMultipleScattering* aMultipleScattering = new G4MuMultipleScattering();
376        pmanager->AddProcess(aMultipleScattering,           -1, 1, 1);
377        pmanager->AddProcess(new G4MuIonisation(),          -1, 2, 2);
378        pmanager->AddProcess(new G4MuBremsstrahlung(),      -1,-1, 3);
379        pmanager->AddProcess(new G4MuPairProduction(),      -1,-1, 4);
380        if( particleName == "mu-" )
381          pmanager->AddProcess(new G4MuonMinusCaptureAtRest(), 0,-1,-1);
382      }
383    else if (particleName == "GenericIon")
384    {
385      pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
386      G4ionIonisation* ionIoni = new G4ionIonisation();
387      ionIoni->SetEmModel(new G4IonParametrisedLossModel());
388      ionIoni->SetStepFunction(0.1, 10*um);
389      pmanager->AddProcess(ionIoni,                   -1, 2, 2);
390      pmanager->AddProcess(new G4NuclearStopping(),   -1, 3,-1);       
391    }
392    else if (particleName == "alpha" || particleName == "He3")
393    {
394      //MSC, ion-Ionisation, Nuclear Stopping
395      pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
396
397      G4ionIonisation* ionIoni = new G4ionIonisation();
398      ionIoni->SetStepFunction(0.1, 20*um);
399      pmanager->AddProcess(ionIoni,                   -1, 2, 2);
400      pmanager->AddProcess(new G4NuclearStopping(),   -1, 3,-1);
401    }
402    else if (particleName == "proton"     ||     
403             particleName == "deuteron"   ||
404             particleName == "triton"     ||
405             particleName == "pi+" ||
406             particleName == "pi-" ||
407             particleName == "kaon+" ||
408             particleName == "kaon-")
409      {
410       //MSC, h-ionisation, bremsstrahlung
411       pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);     
412       G4hIonisation* hIoni = new G4hIonisation();
413       hIoni->SetStepFunction(0.2, 50*um);
414       pmanager->AddProcess(hIoni,                     -1, 2, 2);     
415       pmanager->AddProcess(new G4hBremsstrahlung,     -1,-3, 3);   
416      }
417    else if ((!particle->IsShortLived()) &&
418             (charge != 0.0) &&
419             (particle->GetParticleName() != "chargedgeantino"))
420      {
421        //all others charged particles except geantino
422        pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
423        pmanager->AddProcess(new G4hIonisation,         -1, 2, 2);
424      }
425   
426  }
427}
428
429
430// Optical Processes ////////////////////////////////////////////////////////
431#include "G4Scintillation.hh"
432#include "G4OpAbsorption.hh"
433//#include "G4OpRayleigh.hh"
434#include "G4OpBoundaryProcess.hh"
435
436 template<class T> void TLBE<T>::ConstructOp()
437{
438  // default scintillation process
439  G4Scintillation* theScintProcessDef = new G4Scintillation("Scintillation");
440  // theScintProcessDef->DumpPhysicsTable();
441  theScintProcessDef->SetTrackSecondariesFirst(true);
442  theScintProcessDef->SetScintillationYieldFactor(1.0); //
443  theScintProcessDef->SetScintillationExcitationRatio(0.0); //
444  theScintProcessDef->SetVerboseLevel(OpVerbLevel);
445
446  // scintillation process for alpha:
447  G4Scintillation* theScintProcessAlpha = new G4Scintillation("Scintillation");
448  // theScintProcessNuc->DumpPhysicsTable();
449  theScintProcessAlpha->SetTrackSecondariesFirst(true);
450  theScintProcessAlpha->SetScintillationYieldFactor(1.1);
451  theScintProcessAlpha->SetScintillationExcitationRatio(1.0);
452  theScintProcessAlpha->SetVerboseLevel(OpVerbLevel);
453
454  // scintillation process for heavy nuclei
455  G4Scintillation* theScintProcessNuc = new G4Scintillation("Scintillation");
456  // theScintProcessNuc->DumpPhysicsTable();
457  theScintProcessNuc->SetTrackSecondariesFirst(true);
458  theScintProcessNuc->SetScintillationYieldFactor(0.2);
459  theScintProcessNuc->SetScintillationExcitationRatio(1.0);
460  theScintProcessNuc->SetVerboseLevel(OpVerbLevel);
461
462  // optical processes
463  G4OpAbsorption* theAbsorptionProcess = new G4OpAbsorption();
464  //  G4OpRayleigh* theRayleighScatteringProcess = new G4OpRayleigh();
465  G4OpBoundaryProcess* theBoundaryProcess = new G4OpBoundaryProcess();
466  //  theAbsorptionProcess->DumpPhysicsTable();
467  //  theRayleighScatteringProcess->DumpPhysicsTable();
468  theAbsorptionProcess->SetVerboseLevel(OpVerbLevel);
469  // theRayleighScatteringProcess->SetVerboseLevel(OpVerbLevel);
470  theBoundaryProcess->SetVerboseLevel(OpVerbLevel);
471  G4OpticalSurfaceModel themodel = unified;
472  theBoundaryProcess->SetModel(themodel);
473
474  this->theParticleIterator->reset();
475  while( (*(this->theParticleIterator))() )
476    {
477      G4ParticleDefinition* particle = this->theParticleIterator->value();
478      G4ProcessManager* pmanager = particle->GetProcessManager();
479      G4String particleName = particle->GetParticleName();
480      if (theScintProcessDef->IsApplicable(*particle)) {
481        //      if(particle->GetPDGMass() > 5.0*GeV)
482        if(particle->GetParticleName() == "GenericIon") {
483          pmanager->AddProcess(theScintProcessNuc); // AtRestDiscrete
484          pmanager->SetProcessOrderingToLast(theScintProcessNuc,idxAtRest);
485          pmanager->SetProcessOrderingToLast(theScintProcessNuc,idxPostStep);
486        }         
487        else if(particle->GetParticleName() == "alpha") {
488          pmanager->AddProcess(theScintProcessAlpha);
489          pmanager->SetProcessOrderingToLast(theScintProcessAlpha,idxAtRest);
490          pmanager->SetProcessOrderingToLast(theScintProcessAlpha,idxPostStep);
491        }
492        else {
493          pmanager->AddProcess(theScintProcessDef);
494          pmanager->SetProcessOrderingToLast(theScintProcessDef,idxAtRest);
495          pmanager->SetProcessOrderingToLast(theScintProcessDef,idxPostStep);
496        }         
497      }
498     
499      if (particleName == "opticalphoton") {
500        pmanager->AddDiscreteProcess(theAbsorptionProcess);
501        //      pmanager->AddDiscreteProcess(theRayleighScatteringProcess);
502        pmanager->AddDiscreteProcess(theBoundaryProcess);
503      }
504    }
505}
506
507
508// Hadronic processes ////////////////////////////////////////////////////////
509
510// Elastic processes:
511#include "G4HadronElasticProcess.hh"
512
513// Inelastic processes:
514#include "G4PionPlusInelasticProcess.hh"
515#include "G4PionMinusInelasticProcess.hh"
516#include "G4KaonPlusInelasticProcess.hh"
517#include "G4KaonZeroSInelasticProcess.hh"
518#include "G4KaonZeroLInelasticProcess.hh"
519#include "G4KaonMinusInelasticProcess.hh"
520#include "G4ProtonInelasticProcess.hh"
521#include "G4AntiProtonInelasticProcess.hh"
522#include "G4NeutronInelasticProcess.hh"
523#include "G4AntiNeutronInelasticProcess.hh"
524#include "G4DeuteronInelasticProcess.hh"
525#include "G4TritonInelasticProcess.hh"
526#include "G4AlphaInelasticProcess.hh"
527
528// Low-energy Models: < 20GeV
529#include "G4LElastic.hh"
530#include "G4LEPionPlusInelastic.hh"
531#include "G4LEPionMinusInelastic.hh"
532#include "G4LEKaonPlusInelastic.hh"
533#include "G4LEKaonZeroSInelastic.hh"
534#include "G4LEKaonZeroLInelastic.hh"
535#include "G4LEKaonMinusInelastic.hh"
536#include "G4LEProtonInelastic.hh"
537#include "G4LEAntiProtonInelastic.hh"
538#include "G4LENeutronInelastic.hh"
539#include "G4LEAntiNeutronInelastic.hh"
540#include "G4LEDeuteronInelastic.hh"
541#include "G4LETritonInelastic.hh"
542#include "G4LEAlphaInelastic.hh"
543
544// High-energy Models: >20 GeV
545#include "G4HEPionPlusInelastic.hh"
546#include "G4HEPionMinusInelastic.hh"
547#include "G4HEKaonPlusInelastic.hh"
548#include "G4HEKaonZeroInelastic.hh"
549#include "G4HEKaonZeroInelastic.hh"
550#include "G4HEKaonMinusInelastic.hh"
551#include "G4HEProtonInelastic.hh"
552#include "G4HEAntiProtonInelastic.hh"
553#include "G4HENeutronInelastic.hh"
554#include "G4HEAntiNeutronInelastic.hh"
555
556// Neutron high-precision models: <20 MeV
557#include "G4NeutronHPElastic.hh"
558#include "G4NeutronHPElasticData.hh"
559#include "G4NeutronHPCapture.hh"
560#include "G4NeutronHPCaptureData.hh"
561#include "G4NeutronHPInelastic.hh"
562#include "G4NeutronHPInelasticData.hh"
563#include "G4LCapture.hh"
564
565// Stopping processes
566#include "G4PiMinusAbsorptionAtRest.hh"
567#include "G4KaonMinusAbsorptionAtRest.hh"
568#include "G4AntiProtonAnnihilationAtRest.hh"
569#include "G4AntiNeutronAnnihilationAtRest.hh"
570
571
572// ConstructHad()
573// Makes discrete physics processes for the hadrons, at present limited
574// to those particles with GHEISHA interactions (INTRC > 0).
575// The processes are: Elastic scattering and Inelastic scattering.
576// F.W.Jones  09-JUL-1998
577 template<class T> void TLBE<T>::ConstructHad()
578{
579  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
580  G4LElastic* theElasticModel = new G4LElastic;
581  theElasticProcess->RegisterMe(theElasticModel);
582 
583  this->theParticleIterator->reset();
584  while ((*(this->theParticleIterator))())
585    {
586      G4ParticleDefinition* particle = this->theParticleIterator->value();
587      G4ProcessManager* pmanager = particle->GetProcessManager();
588      G4String particleName = particle->GetParticleName();
589
590      if (particleName == "pi+")
591        {
592          pmanager->AddDiscreteProcess(theElasticProcess);
593          G4PionPlusInelasticProcess* theInelasticProcess =
594            new G4PionPlusInelasticProcess("inelastic");
595          G4LEPionPlusInelastic* theLEInelasticModel =
596            new G4LEPionPlusInelastic;
597          theInelasticProcess->RegisterMe(theLEInelasticModel);
598          G4HEPionPlusInelastic* theHEInelasticModel =
599            new G4HEPionPlusInelastic;
600          theInelasticProcess->RegisterMe(theHEInelasticModel);
601          pmanager->AddDiscreteProcess(theInelasticProcess);
602        }
603
604      else if (particleName == "pi-")
605        {
606          pmanager->AddDiscreteProcess(theElasticProcess);
607          G4PionMinusInelasticProcess* theInelasticProcess =
608            new G4PionMinusInelasticProcess("inelastic");
609          G4LEPionMinusInelastic* theLEInelasticModel =
610            new G4LEPionMinusInelastic;
611          theInelasticProcess->RegisterMe(theLEInelasticModel);
612          G4HEPionMinusInelastic* theHEInelasticModel =
613            new G4HEPionMinusInelastic;
614          theInelasticProcess->RegisterMe(theHEInelasticModel);
615          pmanager->AddDiscreteProcess(theInelasticProcess);
616          G4String prcNam;
617          pmanager->AddRestProcess(new G4PiMinusAbsorptionAtRest, ordDefault);
618        }
619     
620      else if (particleName == "kaon+")
621        {
622          pmanager->AddDiscreteProcess(theElasticProcess);
623          G4KaonPlusInelasticProcess* theInelasticProcess =
624            new G4KaonPlusInelasticProcess("inelastic");
625          G4LEKaonPlusInelastic* theLEInelasticModel =
626            new G4LEKaonPlusInelastic;
627          theInelasticProcess->RegisterMe(theLEInelasticModel);
628          G4HEKaonPlusInelastic* theHEInelasticModel =
629            new G4HEKaonPlusInelastic;
630          theInelasticProcess->RegisterMe(theHEInelasticModel);
631          pmanager->AddDiscreteProcess(theInelasticProcess);
632        }
633     
634      else if (particleName == "kaon0S")
635        {
636          pmanager->AddDiscreteProcess(theElasticProcess);
637          G4KaonZeroSInelasticProcess* theInelasticProcess =
638            new G4KaonZeroSInelasticProcess("inelastic");
639          G4LEKaonZeroSInelastic* theLEInelasticModel =
640            new G4LEKaonZeroSInelastic;
641          theInelasticProcess->RegisterMe(theLEInelasticModel);
642          G4HEKaonZeroInelastic* theHEInelasticModel =
643            new G4HEKaonZeroInelastic;
644          theInelasticProcess->RegisterMe(theHEInelasticModel);
645          pmanager->AddDiscreteProcess(theInelasticProcess);
646        }
647
648      else if (particleName == "kaon0L")
649        {
650          pmanager->AddDiscreteProcess(theElasticProcess);
651          G4KaonZeroLInelasticProcess* theInelasticProcess =
652            new G4KaonZeroLInelasticProcess("inelastic");
653          G4LEKaonZeroLInelastic* theLEInelasticModel =
654            new G4LEKaonZeroLInelastic;
655          theInelasticProcess->RegisterMe(theLEInelasticModel);
656          G4HEKaonZeroInelastic* theHEInelasticModel =
657            new G4HEKaonZeroInelastic;
658          theInelasticProcess->RegisterMe(theHEInelasticModel);
659          pmanager->AddDiscreteProcess(theInelasticProcess);
660        }
661
662      else if (particleName == "kaon-")
663        {
664          pmanager->AddDiscreteProcess(theElasticProcess);
665          G4KaonMinusInelasticProcess* theInelasticProcess =
666            new G4KaonMinusInelasticProcess("inelastic");
667          G4LEKaonMinusInelastic* theLEInelasticModel =
668            new G4LEKaonMinusInelastic;
669          theInelasticProcess->RegisterMe(theLEInelasticModel);
670          G4HEKaonMinusInelastic* theHEInelasticModel =
671            new G4HEKaonMinusInelastic;
672          theInelasticProcess->RegisterMe(theHEInelasticModel);
673          pmanager->AddDiscreteProcess(theInelasticProcess);
674          pmanager->AddRestProcess(new G4KaonMinusAbsorptionAtRest, ordDefault);
675        }
676
677      else if (particleName == "proton")
678        {
679          pmanager->AddDiscreteProcess(theElasticProcess);
680          G4ProtonInelasticProcess* theInelasticProcess =
681            new G4ProtonInelasticProcess("inelastic");
682          G4LEProtonInelastic* theLEInelasticModel = new G4LEProtonInelastic;
683          theInelasticProcess->RegisterMe(theLEInelasticModel);
684          G4HEProtonInelastic* theHEInelasticModel = new G4HEProtonInelastic;
685          theInelasticProcess->RegisterMe(theHEInelasticModel);
686          pmanager->AddDiscreteProcess(theInelasticProcess);
687        }
688
689      else if (particleName == "anti_proton")
690        {
691          pmanager->AddDiscreteProcess(theElasticProcess);
692          G4AntiProtonInelasticProcess* theInelasticProcess =
693            new G4AntiProtonInelasticProcess("inelastic");
694          G4LEAntiProtonInelastic* theLEInelasticModel =
695            new G4LEAntiProtonInelastic;
696          theInelasticProcess->RegisterMe(theLEInelasticModel);
697          G4HEAntiProtonInelastic* theHEInelasticModel =
698            new G4HEAntiProtonInelastic;
699          theInelasticProcess->RegisterMe(theHEInelasticModel);
700          pmanager->AddDiscreteProcess(theInelasticProcess);
701        }
702
703      else if (particleName == "neutron") {
704        // elastic scattering
705        G4HadronElasticProcess* theNeutronElasticProcess =
706          new G4HadronElasticProcess;
707        G4LElastic* theElasticModel1 = new G4LElastic;
708        G4NeutronHPElastic * theElasticNeutron = new G4NeutronHPElastic;
709        theNeutronElasticProcess->RegisterMe(theElasticModel1);
710        theElasticModel1->SetMinEnergy(19*MeV);
711        theNeutronElasticProcess->RegisterMe(theElasticNeutron);
712        G4NeutronHPElasticData * theNeutronData = new G4NeutronHPElasticData;
713        theNeutronElasticProcess->AddDataSet(theNeutronData);
714        pmanager->AddDiscreteProcess(theNeutronElasticProcess);
715        // inelastic scattering
716        G4NeutronInelasticProcess* theInelasticProcess =
717          new G4NeutronInelasticProcess("inelastic");
718        G4LENeutronInelastic* theInelasticModel = new G4LENeutronInelastic;
719        theInelasticModel->SetMinEnergy(19*MeV);
720        theInelasticProcess->RegisterMe(theInelasticModel);
721        G4NeutronHPInelastic * theLENeutronInelasticModel =
722          new G4NeutronHPInelastic;
723        theInelasticProcess->RegisterMe(theLENeutronInelasticModel);
724        G4NeutronHPInelasticData * theNeutronData1 =
725          new G4NeutronHPInelasticData;
726        theInelasticProcess->AddDataSet(theNeutronData1);
727        pmanager->AddDiscreteProcess(theInelasticProcess);
728        // capture
729        G4HadronCaptureProcess* theCaptureProcess =
730          new G4HadronCaptureProcess;
731        G4LCapture* theCaptureModel = new G4LCapture;
732        theCaptureModel->SetMinEnergy(19*MeV);
733        theCaptureProcess->RegisterMe(theCaptureModel);
734        G4NeutronHPCapture * theLENeutronCaptureModel = new G4NeutronHPCapture;
735        theCaptureProcess->RegisterMe(theLENeutronCaptureModel);
736        G4NeutronHPCaptureData * theNeutronData3 = new G4NeutronHPCaptureData;
737        theCaptureProcess->AddDataSet(theNeutronData3);
738        pmanager->AddDiscreteProcess(theCaptureProcess);
739        //  G4ProcessManager* pmanager = G4Neutron::Neutron->GetProcessManager();
740        //  pmanager->AddProcess(new G4UserSpecialCuts(),-1,-1,1);
741      }
742      else if (particleName == "anti_neutron")
743        {
744          pmanager->AddDiscreteProcess(theElasticProcess);
745          G4AntiNeutronInelasticProcess* theInelasticProcess =
746            new G4AntiNeutronInelasticProcess("inelastic");
747          G4LEAntiNeutronInelastic* theLEInelasticModel =
748            new G4LEAntiNeutronInelastic;
749          theInelasticProcess->RegisterMe(theLEInelasticModel);
750          G4HEAntiNeutronInelastic* theHEInelasticModel =
751            new G4HEAntiNeutronInelastic;
752          theInelasticProcess->RegisterMe(theHEInelasticModel);
753          pmanager->AddDiscreteProcess(theInelasticProcess);
754        }
755
756      else if (particleName == "deuteron")
757        {
758          pmanager->AddDiscreteProcess(theElasticProcess);
759          G4DeuteronInelasticProcess* theInelasticProcess =
760            new G4DeuteronInelasticProcess("inelastic");
761          G4LEDeuteronInelastic* theLEInelasticModel =
762            new G4LEDeuteronInelastic;
763          theInelasticProcess->RegisterMe(theLEInelasticModel);
764          pmanager->AddDiscreteProcess(theInelasticProcess);
765        }
766     
767      else if (particleName == "triton")
768        {
769          pmanager->AddDiscreteProcess(theElasticProcess);
770          G4TritonInelasticProcess* theInelasticProcess =
771            new G4TritonInelasticProcess("inelastic");
772          G4LETritonInelastic* theLEInelasticModel =
773            new G4LETritonInelastic;
774          theInelasticProcess->RegisterMe(theLEInelasticModel);
775          pmanager->AddDiscreteProcess(theInelasticProcess);
776        }
777
778      else if (particleName == "alpha")
779        {
780          pmanager->AddDiscreteProcess(theElasticProcess);
781          G4AlphaInelasticProcess* theInelasticProcess =
782            new G4AlphaInelasticProcess("inelastic");
783          G4LEAlphaInelastic* theLEInelasticModel =
784            new G4LEAlphaInelastic;
785          theInelasticProcess->RegisterMe(theLEInelasticModel);
786          pmanager->AddDiscreteProcess(theInelasticProcess);
787        }
788
789    }
790}
791
792
793// Decays ///////////////////////////////////////////////////////////////////
794#include "G4Decay.hh"
795#include "G4RadioactiveDecay.hh"
796#include "G4IonTable.hh"
797#include "G4Ions.hh"
798
799 template<class T> void TLBE<T>::ConstructGeneral() {
800
801  // Add Decay Process
802  G4Decay* theDecayProcess = new G4Decay();
803  this->theParticleIterator->reset();
804  while( (*(this->theParticleIterator))() )
805    {
806      G4ParticleDefinition* particle = this->theParticleIterator->value();
807      G4ProcessManager* pmanager = particle->GetProcessManager();
808     
809      if (theDecayProcess->IsApplicable(*particle) && !particle->IsShortLived())
810        {
811          pmanager ->AddProcess(theDecayProcess);
812          // set ordering for PostStepDoIt and AtRestDoIt
813          pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
814          pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
815        }
816    }
817
818  // Declare radioactive decay to the GenericIon in the IonTable.
819  const G4IonTable *theIonTable =
820    G4ParticleTable::GetParticleTable()->GetIonTable();
821  G4RadioactiveDecay *theRadioactiveDecay = new G4RadioactiveDecay();
822
823  for (G4int i=0; i<theIonTable->Entries(); i++)
824    {
825      G4String particleName = theIonTable->GetParticle(i)->GetParticleName();
826      G4String particleType = theIonTable->GetParticle(i)->GetParticleType();
827     
828      if (particleName == "GenericIon")
829        {
830          G4ProcessManager* pmanager =
831            theIonTable->GetParticle(i)->GetProcessManager();
832          pmanager->SetVerboseLevel(VerboseLevel);
833          pmanager ->AddProcess(theRadioactiveDecay);
834          pmanager ->SetProcessOrdering(theRadioactiveDecay, idxPostStep);
835          pmanager ->SetProcessOrdering(theRadioactiveDecay, idxAtRest);
836        }
837    }
838}
839
840// Cuts /////////////////////////////////////////////////////////////////////
841template<class T> void TLBE<T>::SetCuts()
842{
843 
844  if (this->verboseLevel >1)
845    G4cout << "LBE::SetCuts:";
846 
847  if (this->verboseLevel>0){
848    G4cout << "LBE::SetCuts:";
849    G4cout << "CutLength : "
850           << G4BestUnit(this->defaultCutValue,"Length") << G4endl;
851  }
852
853  //special for low energy physics
854  G4double lowlimit=250*eV; 
855  G4ProductionCutsTable * aPCTable = G4ProductionCutsTable::GetProductionCutsTable();
856  aPCTable->SetEnergyRange(lowlimit,100*GeV);
857   
858  // set cut values for gamma at first and for e- second and next for e+,
859  // because some processes for e+/e- need cut values for gamma
860  this->SetCutValue(cutForGamma, "gamma");
861  this->SetCutValue(cutForElectron, "e-");
862  this->SetCutValue(cutForPositron, "e+");
863 
864  //  this->SetCutValue(cutForProton, "proton");
865  //  this->SetCutValue(cutForProton, "anti_proton");
866  //  this->SetCutValue(cutForAlpha,  "alpha");
867  //  this->SetCutValue(cutForGenericIon,  "GenericIon");
868 
869  //  this->SetCutValueForOthers(this->defaultCutValue);
870 
871  if (this->verboseLevel>0) this->DumpCutValuesTable();
872}
873
Note: See TracBrowser for help on using the repository browser.