Geant4  10.02.p01
Em10PhysicsList.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
28 //
29 //
30 // $Id: Em10PhysicsList.cc 90463 2015-06-01 09:33:38Z gcosmo $
31 //
32 
33 #include "Em10PhysicsList.hh"
36 
37 #include "G4ParticleDefinition.hh"
38 #include "G4ProcessManager.hh"
39 #include "G4ParticleTable.hh"
40 #include "G4ParticleTypes.hh"
41 #include "G4Material.hh"
42 
43 #include "G4UnitsTable.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "G4ios.hh"
46 #include <iomanip>
47 
48 #include "G4Region.hh"
49 #include "G4RegionStore.hh"
50 
51 #include "G4ProductionCuts.hh"
52 #include "G4EmProcessOptions.hh"
53 
54 #include "G4ComptonScattering.hh"
55 #include "G4GammaConversion.hh"
56 #include "G4PhotoElectricEffect.hh"
57 
58 #include "G4eMultipleScattering.hh"
60 #include "G4hMultipleScattering.hh"
61 
62 #include "G4eIonisation.hh"
63 #include "G4eBremsstrahlung.hh"
64 #include "G4eplusAnnihilation.hh"
65 #include "G4PAIModel.hh"
66 #include "G4PAIPhotModel.hh"
67 
69 
70 #include "G4MuIonisation.hh"
71 #include "G4MuBremsstrahlung.hh"
72 #include "G4MuPairProduction.hh"
73 
74 #include "G4hIonisation.hh"
75 
76 #include "G4Decay.hh"
77 
78 #include "G4VXTRenergyLoss.hh"
79 #include "G4RegularXTRadiator.hh"
81 #include "G4GammaXTRadiator.hh"
82 #include "G4StrawTubeXTRadiator.hh"
83 
84 #include "G4XTRGammaRadModel.hh"
85 #include "G4XTRRegularRadModel.hh"
88 
89 #include "Em10StepCut.hh"
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
93 
96  MaxChargedStep(DBL_MAX),
97  theeminusStepCut(0), theeplusStepCut(0),
98  fRadiatorCuts(0),fDetectorCuts(0),fXTRModel("transpM")
99 {
100  pDet = p;
101 
102  // world cuts
103 
104  defaultCutValue = 1.*mm;
108 
109  // Region cuts
110 
114 
115  SetVerboseLevel(1);
117 }
118 
120 {
121  delete physicsListMessenger;
122 }
123 
124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
125 
127 {
128  // In this method, static member functions should be called
129  // for all particles which you want to use.
130  // This ensures that objects of these particle types will be
131  // created in the program.
132 
133  ConstructBosons();
135  ConstructMesons();
137 }
138 
139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
140 
142 {
143  // gamma
145 }
146 
147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
148 
150 {
151  // leptons
152 
157 
162 }
163 
164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
165 
167 {
168  // mesons
169 
175 }
176 
177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
178 
180 {
181 // barions
182 
185 }
186 
187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
188 
190 {
192  ConstructEM();
194 }
195 
196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
197 
199 {
200 
201  // G4cout<<"fMinElectronEnergy = "<<fMinElectronEnergy/keV<<" keV"<<G4endl;
202  // G4cout<<"fMinGammaEnergy = "<<fMinGammaEnergy/keV<<" keV"<<G4endl;
203  G4cout<<"XTR model = "<<fXTRModel<<G4endl;
204 
205  const G4RegionStore* theRegionStore = G4RegionStore::GetInstance();
206  G4Region* gas = theRegionStore->GetRegion("XTRdEdxDetector");
207 
208  G4VXTRenergyLoss* processXTR = 0;
209 
210  if(fXTRModel == "gammaR" )
211  {
212  // G4GammaXTRadiator*
213  processXTR = new G4GammaXTRadiator(pDet->GetLogicalRadiator(),
214  100.,
215  100.,
217  pDet->GetGasMaterial(),
218  pDet->GetFoilThick(),
219  pDet->GetGasThick(),
220  pDet->GetFoilNumber(),
221  "GammaXTRadiator");
222  }
223  else if(fXTRModel == "gammaM" )
224  {
225  // G4XTRGammaRadModel*
226  processXTR = new G4XTRGammaRadModel(pDet->GetLogicalRadiator(),
227  100.,
228  100.,
230  pDet->GetGasMaterial(),
231  pDet->GetFoilThick(),
232  pDet->GetGasThick(),
233  pDet->GetFoilNumber(),
234  "GammaXTRadiator");
235  }
236  else if(fXTRModel == "strawR" )
237  {
238  // G4StrawTubeXTRadiator*
239  processXTR = new G4StrawTubeXTRadiator(pDet->GetLogicalRadiator(),
241  pDet->GetGasMaterial(),
242  0.53, // pDet->GetFoilThick(),
243  3.14159, // pDet->GetGasThick(),
245  true,
246  "strawXTRadiator");
247  }
248  else if(fXTRModel == "regR" )
249  {
250  // G4RegularXTRadiator*
251  processXTR = new G4RegularXTRadiator(pDet->GetLogicalRadiator(),
253  pDet->GetGasMaterial(),
254  pDet->GetFoilThick(),
255  pDet->GetGasThick(),
256  pDet->GetFoilNumber(),
257  "RegularXTRadiator");
258  }
259  else if(fXTRModel == "transpR" )
260  {
261  // G4TransparentRegXTRadiator*
264  pDet->GetGasMaterial(),
265  pDet->GetFoilThick(),
266  pDet->GetGasThick(),
267  pDet->GetFoilNumber(),
268  "RegularXTRadiator");
269  }
270  else if(fXTRModel == "regM" )
271  {
272  // G4XTRRegularRadModel*
273  processXTR = new G4XTRRegularRadModel(pDet->GetLogicalRadiator(),
275  pDet->GetGasMaterial(),
276  pDet->GetFoilThick(),
277  pDet->GetGasThick(),
278  pDet->GetFoilNumber(),
279  "RegularXTRadiator");
280 
281  }
282  else if(fXTRModel == "transpM" )
283  {
284  // G4XTRTransparentRegRadModel*
285  // processXTR = new G4XTRTransparentRegRadModel(pDet->GetLogicalRadiator(),
288  pDet->GetGasMaterial(),
289  pDet->GetFoilThick(),
290  pDet->GetGasThick(),
291  pDet->GetFoilNumber(),
292  "RegularXTRadiator");
293  }
294  else
295  {
296  G4Exception("Invalid XTR model name", "InvalidSetup",
297  FatalException, "XTR model name is out of the name list");
298  }
299  // processXTR->SetCompton(true);
300  processXTR->SetVerboseLevel(1);
301 
302  theParticleIterator->reset();
303 
304  while( (*theParticleIterator)() )
305  {
306  G4ParticleDefinition* particle = theParticleIterator->value();
307  G4ProcessManager* pmanager = particle->GetProcessManager();
308  G4String particleName = particle->GetParticleName();
309 
310  if (particleName == "gamma")
311  {
312  // Construct processes for gamma
313 
316  pmanager->AddDiscreteProcess(new G4GammaConversion);
317 
318  }
319  else if (particleName == "e-")
320  {
321  // Construct processes for electron
324  G4eIonisation* eioni = new G4eIonisation();
325  G4PAIModel* pai = new G4PAIModel(particle,"PAIModel");
326  eioni->AddEmModel(0,pai,pai,gas);
327 
328  pmanager->AddProcess(new G4eMultipleScattering,-1,1,1);
329  pmanager->AddProcess(eioni,-1,2,2);
330  pmanager->AddProcess(new G4eBremsstrahlung,-1,3,3);
331  pmanager->AddDiscreteProcess(processXTR);
334 
335  }
336  else if (particleName == "e+")
337  {
338  // Construct processes for positron
339 
342  G4eIonisation* eioni = new G4eIonisation();
343  G4PAIModel* pai = new G4PAIModel(particle,"PAIModel");
344  eioni->AddEmModel(0,pai,pai,gas);
345 
346  pmanager->AddProcess(new G4eMultipleScattering,-1,1,1);
347  pmanager->AddProcess(eioni,-1,2,2);
348  pmanager->AddProcess(new G4eBremsstrahlung,-1,3,3);
349  pmanager->AddProcess(new G4eplusAnnihilation,0,-1,4);
350  pmanager->AddDiscreteProcess(processXTR);
353 
354  }
355  else if( particleName == "mu+" ||
356  particleName == "mu-" )
357  {
358  // Construct processes for muon+
359 
360  Em10StepCut* muonStepCut = new Em10StepCut();
361  muonStepCut->SetMaxStep(MaxChargedStep) ;
362 
363  G4MuIonisation* muioni = new G4MuIonisation() ;
364 
365  G4PAIModel* pai = new G4PAIModel(particle,"PAIModel");
366  muioni->AddEmModel(0,pai,pai,gas);
367 
368  pmanager->AddProcess(new G4MuMultipleScattering(),-1,1,1);
369  pmanager->AddProcess(muioni,-1,2,2);
370  pmanager->AddProcess(new G4MuBremsstrahlung(),-1,3,3);
371  pmanager->AddProcess(new G4MuPairProduction(),-1,4,4);
372  pmanager->AddProcess( muonStepCut,-1,-1,5);
373 
374  }
375  else if (
376  particleName == "proton"
377  || particleName == "antiproton"
378  || particleName == "pi+"
379  || particleName == "pi-"
380  || particleName == "kaon+"
381  || particleName == "kaon-"
382  )
383  {
384  Em10StepCut* thehadronStepCut = new Em10StepCut();
385  thehadronStepCut->SetMaxStep(MaxChargedStep) ;
386 
387  G4hIonisation* thehIonisation = new G4hIonisation();
388  G4PAIModel* pai = new G4PAIModel(particle,"PAIModel");
389  thehIonisation->AddEmModel(0,pai,pai,gas);
390 
391  pmanager->AddProcess(new G4hMultipleScattering,-1,1,1);
392  pmanager->AddProcess(thehIonisation,-1,2,2);
393  pmanager->AddProcess( thehadronStepCut,-1,-1,3);
394 
395  }
396  }
397  G4EmProcessOptions opt;
398  opt.SetApplyCuts(true);
399 }
400 
401 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
402 
404 {
405  // Add Decay Process
406 
407  G4Decay* theDecayProcess = new G4Decay();
408  theParticleIterator->reset();
409 
410  while( (*theParticleIterator)() )
411  {
412  G4ParticleDefinition* particle = theParticleIterator->value();
413  G4ProcessManager* pmanager = particle->GetProcessManager();
414 
415  if (theDecayProcess->IsApplicable(*particle))
416  {
417  pmanager ->AddProcess(theDecayProcess);
418 
419  // set ordering for PostStepDoIt and AtRestDoIt
420 
421  pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
422  pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
423  }
424  }
425 }
426 
427 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
428 
430 {
431  // set cut values for gamma at first and for e- second and next for e+,
432  // because some processes for e+/e- need cut values for gamma
433 
434  SetCutValue(cutForGamma, "gamma", "DefaultRegionForTheWorld");
435  SetCutValue(cutForElectron, "e-", "DefaultRegionForTheWorld");
436  SetCutValue(cutForPositron, "e+", "DefaultRegionForTheWorld");
437 
438  if (verboseLevel > 0)
439  {
440  G4cout << "Em10PhysicsList::SetCuts:";
441  G4cout << "CutLength for e-, e+ and gamma is: "
442  << G4BestUnit(defaultCutValue,"Length") << G4endl;
443  }
444 
446 
447  G4Region* region = (G4RegionStore::GetInstance())->GetRegion("XTRradiator");
449  G4cout << "Radiator cuts are set" << G4endl;
450 
452  region = (G4RegionStore::GetInstance())->GetRegion("XTRdEdxDetector");
454  G4cout << "Detector cuts are set" << G4endl;
455 
456  if (verboseLevel > 1) DumpCutValuesTable();
457 }
458 
459 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
460 
462 {
463  cutForGamma = val;
464 }
465 
466 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
467 
469 {
470  cutForElectron = val;
471 }
472 
473 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
474 
476 {
477  MaxChargedStep = step ;
478  G4cout << " MaxChargedStep=" << MaxChargedStep << G4endl;
479  G4cout << G4endl;
480 }
481 
482 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
483 
485 {
487 
491 
492  G4cout<<"Radiator gamma cut = "<<fGammaCut/mm<<" mm"<<G4endl;
493  G4cout<<"Radiator electron cut = "<<fElectronCut/mm<<" mm"<<G4endl;
494  G4cout<<"Radiator positron cut = "<<fPositronCut/mm<<" mm"<<G4endl;
495 }
496 
497 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
498 
500 {
502 
506 
507  G4cout<<"Detector gamma cut = "<<fGammaCut/mm<<" mm"<<G4endl;
508  G4cout<<"Detector electron cut = "<<fElectronCut/mm<<" mm"<<G4endl;
509  G4cout<<"Detector positron cut = "<<fPositronCut/mm<<" mm"<<G4endl;
510 
511 }
512 
513 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
G4ProductionCuts * fDetectorCuts
static G4KaonPlus * KaonPlusDefinition()
Definition: G4KaonPlus.cc:108
G4Region * GetRegion(const G4String &name, G4bool verbose=true) const
static G4MuonPlus * MuonPlusDefinition()
Definition: G4MuonPlus.cc:94
void SetCutValue(G4double aCut, const G4String &pname)
void SetMaxStep(G4double)
Definition: Em10StepCut.cc:56
virtual G4bool IsApplicable(const G4ParticleDefinition &)
Definition: G4Decay.cc:89
static G4KaonMinus * KaonMinusDefinition()
Definition: G4KaonMinus.cc:108
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
G4int AddDiscreteProcess(G4VProcess *aProcess, G4int ord=ordDefault)
G4double MaxChargedStep
Definition of the Em10StepCut class.
G4double cutForElectron
void SetProductionCut(G4double cut, G4int index=-1)
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
void SetMaxStep(G4double)
static G4AntiProton * AntiProtonDefinition()
Definition: G4AntiProton.cc:88
G4ProcessManager * GetProcessManager() const
static G4PionZero * PionZeroDefinition()
Definition: G4PionZero.cc:103
const G4String & GetParticleName() const
G4ProductionCuts * fRadiatorCuts
static G4RegionStore * GetInstance()
void DumpCutValuesTable(G4int flag=1)
G4GLOB_DLL std::ostream G4cout
static G4PionMinus * PionMinusDefinition()
Definition: G4PionMinus.cc:93
G4int AddProcess(G4VProcess *aProcess, G4int ordAtRestDoIt=ordInActive, G4int ordAlongSteptDoIt=ordInActive, G4int ordPostStepDoIt=ordInActive)
static G4NeutrinoE * NeutrinoEDefinition()
Definition: G4NeutrinoE.cc:80
static G4AntiNeutrinoMu * AntiNeutrinoMuDefinition()
static G4PionPlus * PionPlusDefinition()
Definition: G4PionPlus.cc:93
void SetGammaCut(G4double)
void AddEmModel(G4int, G4VEmModel *, G4VEmFluctuationModel *fluc=0, const G4Region *region=0)
void SetProcessOrdering(G4VProcess *aProcess, G4ProcessVectorDoItIndex idDoIt, G4int ordDoIt=ordDefault)
void SetVerboseLevel(G4int value)
Em10DetectorConstruction * pDet
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4Positron * PositronDefinition()
Definition: G4Positron.cc:89
Definition of the Em10XTRTransparentRegRadModel class.
Em10StepCut * theeminusStepCut
Definition of the Em10PhysicsListMessenger class.
Em10StepCut * theeplusStepCut
Definition of the Em10DetectorConstruction class.
void SetElectronCut(G4double)
static G4MuonMinus * MuonMinusDefinition()
Definition: G4MuonMinus.cc:95
Em10PhysicsListMessenger * physicsListMessenger
Em10PhysicsList(Em10DetectorConstruction *)
#define G4endl
Definition: G4ios.hh:61
void SetProductionCuts(G4ProductionCuts *cut)
static G4AntiNeutrinoE * AntiNeutrinoEDefinition()
Definition of the Em10PhysicsList class.
G4double cutForPositron
double G4double
Definition: G4Types.hh:76
static G4NeutrinoMu * NeutrinoMuDefinition()
Definition: G4NeutrinoMu.cc:80
#define DBL_MAX
Definition: templates.hh:83
static const double mm
Definition: G4SIunits.hh:114
G4LogicalVolume * GetLogicalRadiator()
#define theParticleIterator
void SetVerboseLevel(G4int value)
Definition: G4VProcess.hh:437
void SetApplyCuts(G4bool val)
static G4Gamma * GammaDefinition()
Definition: G4Gamma.cc:81