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

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

cvs update phys-lists-V09-03-03
File size: 29.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//
49// --------------------------------------------------------------
50
51
52#include "globals.hh"
53#include "G4ProcessManager.hh"
54#include "G4ProcessVector.hh"
55
56#include "G4ParticleDefinition.hh"
57#include "G4ParticleWithCuts.hh"
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
238#include "G4LowEnergyRayleigh.hh"
239#include "G4LowEnergyPhotoElectric.hh"
240#include "G4LowEnergyCompton.hh" 
241#include "G4LowEnergyGammaConversion.hh"
242
243
244// e-
245#include "G4LowEnergyIonisation.hh"
246#include "G4LowEnergyBremsstrahlung.hh"
247
248// e+
249#include "G4eIonisation.hh"
250#include "G4eBremsstrahlung.hh"
251#include "G4eplusAnnihilation.hh"
252
253
254// alpha and GenericIon and deuterons, triton, He3:
255#include "G4hLowEnergyIonisation.hh"
256#include "G4EnergyLossTables.hh"
257// hLowEnergyIonisation uses Ziegler 1988 as the default
258
259
260//muon:
261#include "G4MuIonisation.hh"
262#include "G4MuBremsstrahlung.hh"
263#include "G4MuPairProduction.hh"
264#include "G4MuonMinusCaptureAtRest.hh"
265
266//OTHERS:
267//#include "G4hIonisation.hh" // standard hadron ionisation
268
269 template<class T> void TLBE<T>::ConstructEM() {
270
271// processes:
272
273  G4LowEnergyPhotoElectric* lowePhot = new G4LowEnergyPhotoElectric();
274  G4LowEnergyIonisation* loweIon  = new G4LowEnergyIonisation();
275  G4LowEnergyBremsstrahlung* loweBrem = new G4LowEnergyBremsstrahlung();
276
277  // note LowEIon uses proton as basis for its data-base, therefore
278  // cannot specify different LowEnergyIonisation models for different
279  // particles, but can change model globally for Ion, Alpha and Proton.
280
281
282  //fluorescence apply specific cut for fluorescence from photons, electrons
283  //and bremsstrahlung photons:
284  G4double fluorcut = 250*eV;
285  lowePhot->SetCutForLowEnSecPhotons(fluorcut);
286  loweIon->SetCutForLowEnSecPhotons(fluorcut);
287  loweBrem->SetCutForLowEnSecPhotons(fluorcut);
288 
289  // setting tables explicitly for electronic stopping power
290  //  ahadronLowEIon->SetElectronicStoppingPowerModel
291  //  (G4GenericIon::GenericIonDefinition(), "ICRU_R49p") ;
292  //  ahadronLowEIon->SetElectronicStoppingPowerModel
293  //  (G4Proton::ProtonDefinition(), "ICRU_R49p") ;
294
295  // Switch off the Barkas and Bloch corrections
296  //  ahadronLowEIon->SetBarkasOff();
297
298
299  this->theParticleIterator->reset();
300  while( (*(this->theParticleIterator))() ){
301    G4ParticleDefinition* particle = this->theParticleIterator->value();
302    G4ProcessManager* pmanager = particle->GetProcessManager();
303    G4String particleName = particle->GetParticleName();
304    G4String particleType = particle->GetParticleType();
305    G4double charge = particle->GetPDGCharge();
306   
307    if (particleName == "gamma")
308      {
309        //gamma
310        pmanager->AddDiscreteProcess(new G4LowEnergyRayleigh());
311        pmanager->AddDiscreteProcess(lowePhot);
312        pmanager->AddDiscreteProcess(new G4LowEnergyCompton());
313        pmanager->AddDiscreteProcess(new G4LowEnergyGammaConversion());
314      }
315    else if (particleName == "e-")
316      {
317        //electron
318        // process ordering: AddProcess(name, at rest, along step, post step)
319        // -1 = not implemented, then ordering
320        G4eMultipleScattering* aMultipleScattering = new G4eMultipleScattering();
321        pmanager->AddProcess(aMultipleScattering,     -1, 1, 1);
322        pmanager->AddProcess(loweIon,                 -1, 2, 2);
323        pmanager->AddProcess(loweBrem,                -1,-1, 3);
324      }
325    else if (particleName == "e+")
326      {
327        //positron
328        G4eMultipleScattering* aMultipleScattering = new G4eMultipleScattering();
329        pmanager->AddProcess(aMultipleScattering,     -1, 1, 1);
330        pmanager->AddProcess(new G4eIonisation(),     -1, 2, 2);
331        pmanager->AddProcess(new G4eBremsstrahlung(), -1,-1, 3);
332        pmanager->AddProcess(new G4eplusAnnihilation(),0,-1, 4);     
333      }
334    else if( particleName == "mu+" ||
335             particleName == "mu-"    )
336      {
337        //muon 
338        G4MuMultipleScattering* aMultipleScattering = new G4MuMultipleScattering();
339        pmanager->AddProcess(aMultipleScattering,           -1, 1, 1);
340        pmanager->AddProcess(new G4MuIonisation(),          -1, 2, 2);
341        pmanager->AddProcess(new G4MuBremsstrahlung(),      -1,-1, 3);
342        pmanager->AddProcess(new G4MuPairProduction(),      -1,-1, 4);
343        if( particleName == "mu-" )
344          pmanager->AddProcess(new G4MuonMinusCaptureAtRest(), 0,-1,-1);
345      }
346    else if (particleName == "proton"     ||
347             particleName == "alpha"      ||
348             particleName == "deuteron"   ||
349             particleName == "triton"     ||
350             particleName == "He3"        ||
351             particleName == "GenericIon" ||
352             (particleType == "nucleus" && charge != 0))
353      {
354        // OBJECT may be dynamically created as either a GenericIon or nucleus
355        // G4Nucleus exists and therefore has particle type nucleus
356        // genericIon:
357        G4hMultipleScattering* aMultipleScattering = new G4hMultipleScattering();
358        G4hLowEnergyIonisation* ahadronLowEIon = new G4hLowEnergyIonisation();
359        pmanager->AddProcess(aMultipleScattering,-1,1,1);
360        pmanager->AddProcess(ahadronLowEIon,-1,2,2);
361        // ahadronLowEIon->SetNuclearStoppingOff() ;
362        //        ahadronLowEIon->SetNuclearStoppingPowerModel("ICRU_R49") ;
363        //        ahadronLowEIon->SetNuclearStoppingOn() ;
364 
365        //fluorescence switch off for hadrons (for now) PIXE:
366        ahadronLowEIon->SetFluorescence(false);
367      }
368    else if ((!particle->IsShortLived()) &&
369             (charge != 0.0) &&
370             (particle->GetParticleName() != "chargedgeantino"))
371      {
372        //all others charged particles except geantino
373        G4hMultipleScattering* aMultipleScattering = new G4hMultipleScattering();
374        G4hLowEnergyIonisation* ahadronLowEIon = new G4hLowEnergyIonisation();
375        pmanager->AddProcess(aMultipleScattering,-1,1,1);
376        pmanager->AddProcess(ahadronLowEIon,       -1,2,2);     
377        //      pmanager->AddProcess(new G4hIonisation(),       -1,2,2);     
378      }
379   
380  }
381}
382
383
384// Optical Processes ////////////////////////////////////////////////////////
385#include "G4Scintillation.hh"
386#include "G4OpAbsorption.hh"
387//#include "G4OpRayleigh.hh"
388#include "G4OpBoundaryProcess.hh"
389
390 template<class T> void TLBE<T>::ConstructOp()
391{
392  // default scintillation process
393  G4Scintillation* theScintProcessDef = new G4Scintillation("Scintillation");
394  // theScintProcessDef->DumpPhysicsTable();
395  theScintProcessDef->SetTrackSecondariesFirst(true);
396  theScintProcessDef->SetScintillationYieldFactor(1.0); //
397  theScintProcessDef->SetScintillationExcitationRatio(0.0); //
398  theScintProcessDef->SetVerboseLevel(OpVerbLevel);
399
400  // scintillation process for alpha:
401  G4Scintillation* theScintProcessAlpha = new G4Scintillation("Scintillation");
402  // theScintProcessNuc->DumpPhysicsTable();
403  theScintProcessAlpha->SetTrackSecondariesFirst(true);
404  theScintProcessAlpha->SetScintillationYieldFactor(1.1);
405  theScintProcessAlpha->SetScintillationExcitationRatio(1.0);
406  theScintProcessAlpha->SetVerboseLevel(OpVerbLevel);
407
408  // scintillation process for heavy nuclei
409  G4Scintillation* theScintProcessNuc = new G4Scintillation("Scintillation");
410  // theScintProcessNuc->DumpPhysicsTable();
411  theScintProcessNuc->SetTrackSecondariesFirst(true);
412  theScintProcessNuc->SetScintillationYieldFactor(0.2);
413  theScintProcessNuc->SetScintillationExcitationRatio(1.0);
414  theScintProcessNuc->SetVerboseLevel(OpVerbLevel);
415
416  // optical processes
417  G4OpAbsorption* theAbsorptionProcess = new G4OpAbsorption();
418  //  G4OpRayleigh* theRayleighScatteringProcess = new G4OpRayleigh();
419  G4OpBoundaryProcess* theBoundaryProcess = new G4OpBoundaryProcess();
420  //  theAbsorptionProcess->DumpPhysicsTable();
421  //  theRayleighScatteringProcess->DumpPhysicsTable();
422  theAbsorptionProcess->SetVerboseLevel(OpVerbLevel);
423  // theRayleighScatteringProcess->SetVerboseLevel(OpVerbLevel);
424  theBoundaryProcess->SetVerboseLevel(OpVerbLevel);
425  G4OpticalSurfaceModel themodel = unified;
426  theBoundaryProcess->SetModel(themodel);
427
428  this->theParticleIterator->reset();
429  while( (*(this->theParticleIterator))() )
430    {
431      G4ParticleDefinition* particle = this->theParticleIterator->value();
432      G4ProcessManager* pmanager = particle->GetProcessManager();
433      G4String particleName = particle->GetParticleName();
434      if (theScintProcessDef->IsApplicable(*particle)) {
435        //      if(particle->GetPDGMass() > 5.0*GeV)
436        if(particle->GetParticleName() == "GenericIon") {
437          pmanager->AddProcess(theScintProcessNuc); // AtRestDiscrete
438          pmanager->SetProcessOrderingToLast(theScintProcessNuc,idxAtRest);
439          pmanager->SetProcessOrderingToLast(theScintProcessNuc,idxPostStep);
440        }         
441        else if(particle->GetParticleName() == "alpha") {
442          pmanager->AddProcess(theScintProcessAlpha);
443          pmanager->SetProcessOrderingToLast(theScintProcessAlpha,idxAtRest);
444          pmanager->SetProcessOrderingToLast(theScintProcessAlpha,idxPostStep);
445        }
446        else {
447          pmanager->AddProcess(theScintProcessDef);
448          pmanager->SetProcessOrderingToLast(theScintProcessDef,idxAtRest);
449          pmanager->SetProcessOrderingToLast(theScintProcessDef,idxPostStep);
450        }         
451      }
452     
453      if (particleName == "opticalphoton") {
454        pmanager->AddDiscreteProcess(theAbsorptionProcess);
455        //      pmanager->AddDiscreteProcess(theRayleighScatteringProcess);
456        pmanager->AddDiscreteProcess(theBoundaryProcess);
457      }
458    }
459}
460
461
462// Hadronic processes ////////////////////////////////////////////////////////
463
464// Elastic processes:
465#include "G4HadronElasticProcess.hh"
466
467// Inelastic processes:
468#include "G4PionPlusInelasticProcess.hh"
469#include "G4PionMinusInelasticProcess.hh"
470#include "G4KaonPlusInelasticProcess.hh"
471#include "G4KaonZeroSInelasticProcess.hh"
472#include "G4KaonZeroLInelasticProcess.hh"
473#include "G4KaonMinusInelasticProcess.hh"
474#include "G4ProtonInelasticProcess.hh"
475#include "G4AntiProtonInelasticProcess.hh"
476#include "G4NeutronInelasticProcess.hh"
477#include "G4AntiNeutronInelasticProcess.hh"
478#include "G4DeuteronInelasticProcess.hh"
479#include "G4TritonInelasticProcess.hh"
480#include "G4AlphaInelasticProcess.hh"
481
482// Low-energy Models: < 20GeV
483#include "G4LElastic.hh"
484#include "G4LEPionPlusInelastic.hh"
485#include "G4LEPionMinusInelastic.hh"
486#include "G4LEKaonPlusInelastic.hh"
487#include "G4LEKaonZeroSInelastic.hh"
488#include "G4LEKaonZeroLInelastic.hh"
489#include "G4LEKaonMinusInelastic.hh"
490#include "G4LEProtonInelastic.hh"
491#include "G4LEAntiProtonInelastic.hh"
492#include "G4LENeutronInelastic.hh"
493#include "G4LEAntiNeutronInelastic.hh"
494#include "G4LEDeuteronInelastic.hh"
495#include "G4LETritonInelastic.hh"
496#include "G4LEAlphaInelastic.hh"
497
498// High-energy Models: >20 GeV
499#include "G4HEPionPlusInelastic.hh"
500#include "G4HEPionMinusInelastic.hh"
501#include "G4HEKaonPlusInelastic.hh"
502#include "G4HEKaonZeroInelastic.hh"
503#include "G4HEKaonZeroInelastic.hh"
504#include "G4HEKaonMinusInelastic.hh"
505#include "G4HEProtonInelastic.hh"
506#include "G4HEAntiProtonInelastic.hh"
507#include "G4HENeutronInelastic.hh"
508#include "G4HEAntiNeutronInelastic.hh"
509
510// Neutron high-precision models: <20 MeV
511#include "G4NeutronHPElastic.hh"
512#include "G4NeutronHPElasticData.hh"
513#include "G4NeutronHPCapture.hh"
514#include "G4NeutronHPCaptureData.hh"
515#include "G4NeutronHPInelastic.hh"
516#include "G4NeutronHPInelasticData.hh"
517#include "G4LCapture.hh"
518
519// Stopping processes
520#include "G4PiMinusAbsorptionAtRest.hh"
521#include "G4KaonMinusAbsorptionAtRest.hh"
522#include "G4AntiProtonAnnihilationAtRest.hh"
523#include "G4AntiNeutronAnnihilationAtRest.hh"
524
525
526// ConstructHad()
527// Makes discrete physics processes for the hadrons, at present limited
528// to those particles with GHEISHA interactions (INTRC > 0).
529// The processes are: Elastic scattering and Inelastic scattering.
530// F.W.Jones  09-JUL-1998
531 template<class T> void TLBE<T>::ConstructHad()
532{
533  G4HadronElasticProcess* theElasticProcess = new G4HadronElasticProcess;
534  G4LElastic* theElasticModel = new G4LElastic;
535  theElasticProcess->RegisterMe(theElasticModel);
536 
537  this->theParticleIterator->reset();
538  while ((*(this->theParticleIterator))())
539    {
540      G4ParticleDefinition* particle = this->theParticleIterator->value();
541      G4ProcessManager* pmanager = particle->GetProcessManager();
542      G4String particleName = particle->GetParticleName();
543
544      if (particleName == "pi+")
545        {
546          pmanager->AddDiscreteProcess(theElasticProcess);
547          G4PionPlusInelasticProcess* theInelasticProcess =
548            new G4PionPlusInelasticProcess("inelastic");
549          G4LEPionPlusInelastic* theLEInelasticModel =
550            new G4LEPionPlusInelastic;
551          theInelasticProcess->RegisterMe(theLEInelasticModel);
552          G4HEPionPlusInelastic* theHEInelasticModel =
553            new G4HEPionPlusInelastic;
554          theInelasticProcess->RegisterMe(theHEInelasticModel);
555          pmanager->AddDiscreteProcess(theInelasticProcess);
556        }
557
558      else if (particleName == "pi-")
559        {
560          pmanager->AddDiscreteProcess(theElasticProcess);
561          G4PionMinusInelasticProcess* theInelasticProcess =
562            new G4PionMinusInelasticProcess("inelastic");
563          G4LEPionMinusInelastic* theLEInelasticModel =
564            new G4LEPionMinusInelastic;
565          theInelasticProcess->RegisterMe(theLEInelasticModel);
566          G4HEPionMinusInelastic* theHEInelasticModel =
567            new G4HEPionMinusInelastic;
568          theInelasticProcess->RegisterMe(theHEInelasticModel);
569          pmanager->AddDiscreteProcess(theInelasticProcess);
570          G4String prcNam;
571          pmanager->AddRestProcess(new G4PiMinusAbsorptionAtRest, ordDefault);
572        }
573     
574      else if (particleName == "kaon+")
575        {
576          pmanager->AddDiscreteProcess(theElasticProcess);
577          G4KaonPlusInelasticProcess* theInelasticProcess =
578            new G4KaonPlusInelasticProcess("inelastic");
579          G4LEKaonPlusInelastic* theLEInelasticModel =
580            new G4LEKaonPlusInelastic;
581          theInelasticProcess->RegisterMe(theLEInelasticModel);
582          G4HEKaonPlusInelastic* theHEInelasticModel =
583            new G4HEKaonPlusInelastic;
584          theInelasticProcess->RegisterMe(theHEInelasticModel);
585          pmanager->AddDiscreteProcess(theInelasticProcess);
586        }
587     
588      else if (particleName == "kaon0S")
589        {
590          pmanager->AddDiscreteProcess(theElasticProcess);
591          G4KaonZeroSInelasticProcess* theInelasticProcess =
592            new G4KaonZeroSInelasticProcess("inelastic");
593          G4LEKaonZeroSInelastic* theLEInelasticModel =
594            new G4LEKaonZeroSInelastic;
595          theInelasticProcess->RegisterMe(theLEInelasticModel);
596          G4HEKaonZeroInelastic* theHEInelasticModel =
597            new G4HEKaonZeroInelastic;
598          theInelasticProcess->RegisterMe(theHEInelasticModel);
599          pmanager->AddDiscreteProcess(theInelasticProcess);
600        }
601
602      else if (particleName == "kaon0L")
603        {
604          pmanager->AddDiscreteProcess(theElasticProcess);
605          G4KaonZeroLInelasticProcess* theInelasticProcess =
606            new G4KaonZeroLInelasticProcess("inelastic");
607          G4LEKaonZeroLInelastic* theLEInelasticModel =
608            new G4LEKaonZeroLInelastic;
609          theInelasticProcess->RegisterMe(theLEInelasticModel);
610          G4HEKaonZeroInelastic* theHEInelasticModel =
611            new G4HEKaonZeroInelastic;
612          theInelasticProcess->RegisterMe(theHEInelasticModel);
613          pmanager->AddDiscreteProcess(theInelasticProcess);
614        }
615
616      else if (particleName == "kaon-")
617        {
618          pmanager->AddDiscreteProcess(theElasticProcess);
619          G4KaonMinusInelasticProcess* theInelasticProcess =
620            new G4KaonMinusInelasticProcess("inelastic");
621          G4LEKaonMinusInelastic* theLEInelasticModel =
622            new G4LEKaonMinusInelastic;
623          theInelasticProcess->RegisterMe(theLEInelasticModel);
624          G4HEKaonMinusInelastic* theHEInelasticModel =
625            new G4HEKaonMinusInelastic;
626          theInelasticProcess->RegisterMe(theHEInelasticModel);
627          pmanager->AddDiscreteProcess(theInelasticProcess);
628          pmanager->AddRestProcess(new G4KaonMinusAbsorptionAtRest, ordDefault);
629        }
630
631      else if (particleName == "proton")
632        {
633          pmanager->AddDiscreteProcess(theElasticProcess);
634          G4ProtonInelasticProcess* theInelasticProcess =
635            new G4ProtonInelasticProcess("inelastic");
636          G4LEProtonInelastic* theLEInelasticModel = new G4LEProtonInelastic;
637          theInelasticProcess->RegisterMe(theLEInelasticModel);
638          G4HEProtonInelastic* theHEInelasticModel = new G4HEProtonInelastic;
639          theInelasticProcess->RegisterMe(theHEInelasticModel);
640          pmanager->AddDiscreteProcess(theInelasticProcess);
641        }
642
643      else if (particleName == "anti_proton")
644        {
645          pmanager->AddDiscreteProcess(theElasticProcess);
646          G4AntiProtonInelasticProcess* theInelasticProcess =
647            new G4AntiProtonInelasticProcess("inelastic");
648          G4LEAntiProtonInelastic* theLEInelasticModel =
649            new G4LEAntiProtonInelastic;
650          theInelasticProcess->RegisterMe(theLEInelasticModel);
651          G4HEAntiProtonInelastic* theHEInelasticModel =
652            new G4HEAntiProtonInelastic;
653          theInelasticProcess->RegisterMe(theHEInelasticModel);
654          pmanager->AddDiscreteProcess(theInelasticProcess);
655        }
656
657      else if (particleName == "neutron") {
658        // elastic scattering
659        G4HadronElasticProcess* theNeutronElasticProcess =
660          new G4HadronElasticProcess;
661        G4LElastic* theElasticModel1 = new G4LElastic;
662        G4NeutronHPElastic * theElasticNeutron = new G4NeutronHPElastic;
663        theNeutronElasticProcess->RegisterMe(theElasticModel1);
664        theElasticModel1->SetMinEnergy(19*MeV);
665        theNeutronElasticProcess->RegisterMe(theElasticNeutron);
666        G4NeutronHPElasticData * theNeutronData = new G4NeutronHPElasticData;
667        theNeutronElasticProcess->AddDataSet(theNeutronData);
668        pmanager->AddDiscreteProcess(theNeutronElasticProcess);
669        // inelastic scattering
670        G4NeutronInelasticProcess* theInelasticProcess =
671          new G4NeutronInelasticProcess("inelastic");
672        G4LENeutronInelastic* theInelasticModel = new G4LENeutronInelastic;
673        theInelasticModel->SetMinEnergy(19*MeV);
674        theInelasticProcess->RegisterMe(theInelasticModel);
675        G4NeutronHPInelastic * theLENeutronInelasticModel =
676          new G4NeutronHPInelastic;
677        theInelasticProcess->RegisterMe(theLENeutronInelasticModel);
678        G4NeutronHPInelasticData * theNeutronData1 =
679          new G4NeutronHPInelasticData;
680        theInelasticProcess->AddDataSet(theNeutronData1);
681        pmanager->AddDiscreteProcess(theInelasticProcess);
682        // capture
683        G4HadronCaptureProcess* theCaptureProcess =
684          new G4HadronCaptureProcess;
685        G4LCapture* theCaptureModel = new G4LCapture;
686        theCaptureModel->SetMinEnergy(19*MeV);
687        theCaptureProcess->RegisterMe(theCaptureModel);
688        G4NeutronHPCapture * theLENeutronCaptureModel = new G4NeutronHPCapture;
689        theCaptureProcess->RegisterMe(theLENeutronCaptureModel);
690        G4NeutronHPCaptureData * theNeutronData3 = new G4NeutronHPCaptureData;
691        theCaptureProcess->AddDataSet(theNeutronData3);
692        pmanager->AddDiscreteProcess(theCaptureProcess);
693        //  G4ProcessManager* pmanager = G4Neutron::Neutron->GetProcessManager();
694        //  pmanager->AddProcess(new G4UserSpecialCuts(),-1,-1,1);
695      }
696      else if (particleName == "anti_neutron")
697        {
698          pmanager->AddDiscreteProcess(theElasticProcess);
699          G4AntiNeutronInelasticProcess* theInelasticProcess =
700            new G4AntiNeutronInelasticProcess("inelastic");
701          G4LEAntiNeutronInelastic* theLEInelasticModel =
702            new G4LEAntiNeutronInelastic;
703          theInelasticProcess->RegisterMe(theLEInelasticModel);
704          G4HEAntiNeutronInelastic* theHEInelasticModel =
705            new G4HEAntiNeutronInelastic;
706          theInelasticProcess->RegisterMe(theHEInelasticModel);
707          pmanager->AddDiscreteProcess(theInelasticProcess);
708        }
709
710      else if (particleName == "deuteron")
711        {
712          pmanager->AddDiscreteProcess(theElasticProcess);
713          G4DeuteronInelasticProcess* theInelasticProcess =
714            new G4DeuteronInelasticProcess("inelastic");
715          G4LEDeuteronInelastic* theLEInelasticModel =
716            new G4LEDeuteronInelastic;
717          theInelasticProcess->RegisterMe(theLEInelasticModel);
718          pmanager->AddDiscreteProcess(theInelasticProcess);
719        }
720     
721      else if (particleName == "triton")
722        {
723          pmanager->AddDiscreteProcess(theElasticProcess);
724          G4TritonInelasticProcess* theInelasticProcess =
725            new G4TritonInelasticProcess("inelastic");
726          G4LETritonInelastic* theLEInelasticModel =
727            new G4LETritonInelastic;
728          theInelasticProcess->RegisterMe(theLEInelasticModel);
729          pmanager->AddDiscreteProcess(theInelasticProcess);
730        }
731
732      else if (particleName == "alpha")
733        {
734          pmanager->AddDiscreteProcess(theElasticProcess);
735          G4AlphaInelasticProcess* theInelasticProcess =
736            new G4AlphaInelasticProcess("inelastic");
737          G4LEAlphaInelastic* theLEInelasticModel =
738            new G4LEAlphaInelastic;
739          theInelasticProcess->RegisterMe(theLEInelasticModel);
740          pmanager->AddDiscreteProcess(theInelasticProcess);
741        }
742
743    }
744}
745
746
747// Decays ///////////////////////////////////////////////////////////////////
748#include "G4Decay.hh"
749#include "G4RadioactiveDecay.hh"
750#include "G4IonTable.hh"
751#include "G4Ions.hh"
752
753 template<class T> void TLBE<T>::ConstructGeneral() {
754
755  // Add Decay Process
756  G4Decay* theDecayProcess = new G4Decay();
757  this->theParticleIterator->reset();
758  while( (*(this->theParticleIterator))() )
759    {
760      G4ParticleDefinition* particle = this->theParticleIterator->value();
761      G4ProcessManager* pmanager = particle->GetProcessManager();
762     
763      if (theDecayProcess->IsApplicable(*particle) && !particle->IsShortLived())
764        {
765          pmanager ->AddProcess(theDecayProcess);
766          // set ordering for PostStepDoIt and AtRestDoIt
767          pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
768          pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
769        }
770    }
771
772  // Declare radioactive decay to the GenericIon in the IonTable.
773  const G4IonTable *theIonTable =
774    G4ParticleTable::GetParticleTable()->GetIonTable();
775  G4RadioactiveDecay *theRadioactiveDecay = new G4RadioactiveDecay();
776
777  for (G4int i=0; i<theIonTable->Entries(); i++)
778    {
779      G4String particleName = theIonTable->GetParticle(i)->GetParticleName();
780      G4String particleType = theIonTable->GetParticle(i)->GetParticleType();
781     
782      if (particleName == "GenericIon")
783        {
784          G4ProcessManager* pmanager =
785            theIonTable->GetParticle(i)->GetProcessManager();
786          pmanager->SetVerboseLevel(VerboseLevel);
787          pmanager ->AddProcess(theRadioactiveDecay);
788          pmanager ->SetProcessOrdering(theRadioactiveDecay, idxPostStep);
789          pmanager ->SetProcessOrdering(theRadioactiveDecay, idxAtRest);
790        }
791    }
792}
793
794// Cuts /////////////////////////////////////////////////////////////////////
795template<class T> void TLBE<T>::SetCuts()
796{
797 
798  if (this->verboseLevel >1)
799    G4cout << "LBE::SetCuts:";
800 
801  if (this->verboseLevel>0){
802    G4cout << "LBE::SetCuts:";
803    G4cout << "CutLength : "
804           << G4BestUnit(this->defaultCutValue,"Length") << G4endl;
805  }
806
807  //special for low energy physics
808  G4double lowlimit=250*eV; 
809  G4ProductionCutsTable * aPCTable = G4ProductionCutsTable::GetProductionCutsTable();
810  aPCTable->SetEnergyRange(lowlimit,100*GeV);
811   
812  // set cut values for gamma at first and for e- second and next for e+,
813  // because some processes for e+/e- need cut values for gamma
814  this->SetCutValue(cutForGamma, "gamma");
815  this->SetCutValue(cutForElectron, "e-");
816  this->SetCutValue(cutForPositron, "e+");
817 
818  //  this->SetCutValue(cutForProton, "proton");
819  //  this->SetCutValue(cutForProton, "anti_proton");
820  //  this->SetCutValue(cutForAlpha,  "alpha");
821  //  this->SetCutValue(cutForGenericIon,  "GenericIon");
822 
823  //  this->SetCutValueForOthers(this->defaultCutValue);
824 
825  if (this->verboseLevel>0) this->DumpCutValuesTable();
826}
827
Note: See TracBrowser for help on using the repository browser.