Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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$
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 "G4PAIPhotonModel.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 
92 //
93 //
94 
97  MaxChargedStep(DBL_MAX),
98  theeminusStepCut(0), theeplusStepCut(0),
99  fRadiatorCuts(0),fDetectorCuts(0),fXTRModel("transpM")
100 {
101  pDet = p;
102 
103  // world cuts
104 
105  defaultCutValue = 1.*mm;
106  cutForGamma = defaultCutValue;
107  cutForElectron = defaultCutValue;
108  cutForPositron = defaultCutValue;
109 
110  // Region cuts
111 
112  fGammaCut = defaultCutValue;
113  fElectronCut = defaultCutValue;
114  fPositronCut = defaultCutValue;
115 
116  SetVerboseLevel(1);
117  physicsListMessenger = new Em10PhysicsListMessenger(this);
118 }
119 
121 {
122  delete physicsListMessenger;
123 }
124 
126 //
127 //
128 
130 {
131  // In this method, static member functions should be called
132  // for all particles which you want to use.
133  // This ensures that objects of these particle types will be
134  // created in the program.
135 
136  ConstructBosons();
137  ConstructLeptons();
138  ConstructMesons();
139  ConstructBarions();
140 }
141 
143 //
144 //
145 
146 void Em10PhysicsList::ConstructBosons()
147 {
148  // gamma
150 }
151 
152 void Em10PhysicsList::ConstructLeptons()
153 {
154  // leptons
155 
160 
165 }
166 
167 void Em10PhysicsList::ConstructMesons()
168 {
169  // mesons
170 
176 }
177 
178 
179 void Em10PhysicsList::ConstructBarions()
180 {
181 // barions
182 
185 }
186 
187 
189 //
190 //
191 
193 {
195  ConstructEM();
196  ConstructGeneral();
197 }
198 
200 //
201 //
202 
203 void Em10PhysicsList::ConstructEM()
204 {
205 
206  // G4cout<<"fMinElectronEnergy = "<<fMinElectronEnergy/keV<<" keV"<<G4endl;
207  // G4cout<<"fMinGammaEnergy = "<<fMinGammaEnergy/keV<<" keV"<<G4endl;
208  G4cout<<"XTR model = "<<fXTRModel<<G4endl;
209 
210  const G4RegionStore* theRegionStore = G4RegionStore::GetInstance();
211  G4Region* gas = theRegionStore->GetRegion("XTRdEdxDetector");
212 
213  G4VXTRenergyLoss* processXTR = 0;
214 
215  if(fXTRModel == "gammaR" )
216  {
217  // G4GammaXTRadiator*
218  processXTR = new G4GammaXTRadiator(pDet->GetLogicalRadiator(),
219  100.,
220  100.,
221  pDet->GetFoilMaterial(),
222  pDet->GetGasMaterial(),
223  pDet->GetFoilThick(),
224  pDet->GetGasThick(),
225  pDet->GetFoilNumber(),
226  "GammaXTRadiator");
227  }
228  else if(fXTRModel == "gammaM" )
229  {
230  // G4XTRGammaRadModel*
231  processXTR = new G4XTRGammaRadModel(pDet->GetLogicalRadiator(),
232  100.,
233  100.,
234  pDet->GetFoilMaterial(),
235  pDet->GetGasMaterial(),
236  pDet->GetFoilThick(),
237  pDet->GetGasThick(),
238  pDet->GetFoilNumber(),
239  "GammaXTRadiator");
240  }
241  else if(fXTRModel == "strawR" )
242  {
243 
244  // G4StrawTubeXTRadiator*
245  processXTR = new G4StrawTubeXTRadiator(pDet->GetLogicalRadiator(),
246  pDet->GetFoilMaterial(),
247  pDet->GetGasMaterial(),
248  0.53, // pDet->GetFoilThick(),
249  3.14159, // pDet->GetGasThick(),
250  pDet->GetAbsorberMaterial(),
251  true,
252  "strawXTRadiator");
253  }
254  else if(fXTRModel == "regR" )
255  {
256  // G4RegularXTRadiator*
257  processXTR = new G4RegularXTRadiator(pDet->GetLogicalRadiator(),
258  pDet->GetFoilMaterial(),
259  pDet->GetGasMaterial(),
260  pDet->GetFoilThick(),
261  pDet->GetGasThick(),
262  pDet->GetFoilNumber(),
263  "RegularXTRadiator");
264  }
265  else if(fXTRModel == "transpR" )
266  {
267  // G4TransparentRegXTRadiator*
268  processXTR = new G4TransparentRegXTRadiator(pDet->GetLogicalRadiator(),
269  pDet->GetFoilMaterial(),
270  pDet->GetGasMaterial(),
271  pDet->GetFoilThick(),
272  pDet->GetGasThick(),
273  pDet->GetFoilNumber(),
274  "RegularXTRadiator");
275  }
276  else if(fXTRModel == "regM" )
277  {
278  // G4XTRRegularRadModel*
279  processXTR = new G4XTRRegularRadModel(pDet->GetLogicalRadiator(),
280  pDet->GetFoilMaterial(),
281  pDet->GetGasMaterial(),
282  pDet->GetFoilThick(),
283  pDet->GetGasThick(),
284  pDet->GetFoilNumber(),
285  "RegularXTRadiator");
286 
287  }
288  else if(fXTRModel == "transpM" )
289  {
290  // G4XTRTransparentRegRadModel*
291  // processXTR = new G4XTRTransparentRegRadModel(pDet->GetLogicalRadiator(),
292  processXTR = new Em10XTRTransparentRegRadModel(pDet->GetLogicalRadiator(),
293  pDet->GetFoilMaterial(),
294  pDet->GetGasMaterial(),
295  pDet->GetFoilThick(),
296  pDet->GetGasThick(),
297  pDet->GetFoilNumber(),
298  "RegularXTRadiator");
299  }
300  else
301  {
302  G4Exception("Invalid XTR model name", "InvalidSetup",
303  FatalException, "XTR model name is out of the name list");
304  }
305  // processXTR->SetCompton(true);
306  processXTR->SetVerboseLevel(1);
307 
309 
310  while( (*theParticleIterator)() )
311  {
313  G4ProcessManager* pmanager = particle->GetProcessManager();
314  G4String particleName = particle->GetParticleName();
315 
316  if (particleName == "gamma")
317  {
318  // Construct processes for gamma
319 
322  pmanager->AddDiscreteProcess(new G4GammaConversion);
323 
324  }
325  else if (particleName == "e-")
326  {
327  // Construct processes for electron
328  theeminusStepCut = new Em10StepCut();
329  theeminusStepCut->SetMaxStep(MaxChargedStep) ;
330  G4eIonisation* eioni = new G4eIonisation();
331  G4PAIModel* pai = new G4PAIModel(particle,"PAIModel");
332  eioni->AddEmModel(0,pai,pai,gas);
333 
334  pmanager->AddProcess(new G4eMultipleScattering,-1,1,1);
335  pmanager->AddProcess(eioni,-1,2,2);
336  pmanager->AddProcess(new G4eBremsstrahlung,-1,3,3);
337  pmanager->AddDiscreteProcess(processXTR);
339  pmanager->AddDiscreteProcess(theeminusStepCut);
340 
341  }
342  else if (particleName == "e+")
343  {
344  // Construct processes for positron
345 
346  theeplusStepCut = new Em10StepCut();
347  theeplusStepCut->SetMaxStep(MaxChargedStep) ;
348  G4eIonisation* eioni = new G4eIonisation();
349  G4PAIModel* pai = new G4PAIModel(particle,"PAIModel");
350  eioni->AddEmModel(0,pai,pai,gas);
351 
352  pmanager->AddProcess(new G4eMultipleScattering,-1,1,1);
353  pmanager->AddProcess(eioni,-1,2,2);
354  pmanager->AddProcess(new G4eBremsstrahlung,-1,3,3);
355  pmanager->AddProcess(new G4eplusAnnihilation,0,-1,4);
356  pmanager->AddDiscreteProcess(processXTR);
358  pmanager->AddDiscreteProcess(theeplusStepCut);
359 
360  }
361  else if( particleName == "mu+" ||
362  particleName == "mu-" )
363  {
364  // Construct processes for muon+
365 
366  Em10StepCut* muonStepCut = new Em10StepCut();
367  muonStepCut->SetMaxStep(MaxChargedStep) ;
368 
369  G4MuIonisation* muioni = new G4MuIonisation() ;
370 
371  G4PAIModel* pai = new G4PAIModel(particle,"PAIModel");
372  muioni->AddEmModel(0,pai,pai,gas);
373 
374  pmanager->AddProcess(new G4MuMultipleScattering(),-1,1,1);
375  pmanager->AddProcess(muioni,-1,2,2);
376  pmanager->AddProcess(new G4MuBremsstrahlung(),-1,3,3);
377  pmanager->AddProcess(new G4MuPairProduction(),-1,4,4);
378  pmanager->AddProcess( muonStepCut,-1,-1,5);
379 
380  }
381  else if (
382  particleName == "proton"
383  || particleName == "antiproton"
384  || particleName == "pi+"
385  || particleName == "pi-"
386  || particleName == "kaon+"
387  || particleName == "kaon-"
388  )
389  {
390  Em10StepCut* thehadronStepCut = new Em10StepCut();
391  thehadronStepCut->SetMaxStep(MaxChargedStep) ;
392 
393  G4hIonisation* thehIonisation = new G4hIonisation();
394  G4PAIModel* pai = new G4PAIModel(particle,"PAIModel");
395  thehIonisation->AddEmModel(0,pai,pai,gas);
396 
397  pmanager->AddProcess(new G4hMultipleScattering,-1,1,1);
398  pmanager->AddProcess(thehIonisation,-1,2,2);
399  pmanager->AddProcess( thehadronStepCut,-1,-1,3);
400 
401  }
402  }
403  G4EmProcessOptions opt;
404  opt.SetApplyCuts(true);
405 }
406 
407 void Em10PhysicsList::ConstructGeneral()
408 {
409  // Add Decay Process
410 
411  G4Decay* theDecayProcess = new G4Decay();
413 
414  while( (*theParticleIterator)() )
415  {
417  G4ProcessManager* pmanager = particle->GetProcessManager();
418 
419  if (theDecayProcess->IsApplicable(*particle))
420  {
421  pmanager ->AddProcess(theDecayProcess);
422 
423  // set ordering for PostStepDoIt and AtRestDoIt
424 
425  pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
426  pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
427  }
428  }
429 }
430 
432 
434 {
435  // set cut values for gamma at first and for e- second and next for e+,
436  // because some processes for e+/e- need cut values for gamma
437 
438  SetCutValue(cutForGamma, "gamma", "DefaultRegionForTheWorld");
439  SetCutValue(cutForElectron, "e-", "DefaultRegionForTheWorld");
440  SetCutValue(cutForPositron, "e+", "DefaultRegionForTheWorld");
441 
442  if (verboseLevel > 0)
443  {
444  G4cout << "Em10PhysicsList::SetCuts:";
445  G4cout << "CutLength for e-, e+ and gamma is: "
446  << G4BestUnit(defaultCutValue,"Length") << G4endl;
447  }
448 
449  if( !fRadiatorCuts ) SetRadiatorCuts();
450 
451  G4Region* region = (G4RegionStore::GetInstance())->GetRegion("XTRradiator");
452  region->SetProductionCuts(fRadiatorCuts);
453  G4cout << "Radiator cuts are set" << G4endl;
454 
455  if( !fDetectorCuts ) SetDetectorCuts();
456  region = (G4RegionStore::GetInstance())->GetRegion("XTRdEdxDetector");
457  region->SetProductionCuts(fDetectorCuts);
458  G4cout << "Detector cuts are set" << G4endl;
459 
460  if (verboseLevel > 1) DumpCutValuesTable();
461 }
462 
464 
466 {
467  cutForGamma = val;
468 }
469 
471 
473 {
474  cutForElectron = val;
475 }
476 
478 
480 {
481  MaxChargedStep = step ;
482  G4cout << " MaxChargedStep=" << MaxChargedStep << G4endl;
483  G4cout << G4endl;
484 }
485 
487 
489 {
490  if( !fRadiatorCuts ) fRadiatorCuts = new G4ProductionCuts();
491 
492  fRadiatorCuts->SetProductionCut(fGammaCut, idxG4GammaCut);
493  fRadiatorCuts->SetProductionCut(fElectronCut, idxG4ElectronCut);
494  fRadiatorCuts->SetProductionCut(fPositronCut, idxG4PositronCut);
495 
496  G4cout<<"Radiator gamma cut = "<<fGammaCut/mm<<" mm"<<G4endl;
497  G4cout<<"Radiator electron cut = "<<fElectronCut/mm<<" mm"<<G4endl;
498  G4cout<<"Radiator positron cut = "<<fPositronCut/mm<<" mm"<<G4endl;
499 }
500 
502 
504 {
505  if( !fDetectorCuts ) fDetectorCuts = new G4ProductionCuts();
506 
507  fDetectorCuts->SetProductionCut(fGammaCut, idxG4GammaCut);
508  fDetectorCuts->SetProductionCut(fElectronCut, idxG4ElectronCut);
509  fDetectorCuts->SetProductionCut(fPositronCut, idxG4PositronCut);
510 
511  G4cout<<"Detector gamma cut = "<<fGammaCut/mm<<" mm"<<G4endl;
512  G4cout<<"Detector electron cut = "<<fElectronCut/mm<<" mm"<<G4endl;
513  G4cout<<"Detector positron cut = "<<fPositronCut/mm<<" mm"<<G4endl;
514 
515 }