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

Last change on this file since 850 was 825, checked in by garnier, 16 years ago

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