source: trunk/examples/extended/electromagnetic/TestEm10/src/Em10PhysicsList.cc @ 1339

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

tag geant4.9.4 beta 1 + modifs locales
File size: 15.1 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// $Id: Em10PhysicsList.cc,v 1.25 2009/11/21 16:12:01 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-01 $
29//
30
31#include "Em10PhysicsList.hh"
32#include "Em10DetectorConstruction.hh"
33#include "Em10PhysicsListMessenger.hh"
34
35#include "G4ParticleDefinition.hh"
36#include "G4ProcessManager.hh"
37#include "G4ParticleTable.hh"
38#include "G4ParticleTypes.hh"
39#include "G4Material.hh"
40
41#include "G4UnitsTable.hh"
42#include "G4ios.hh"
43#include <iomanip>
44
45#include "G4Region.hh"
46#include "G4RegionStore.hh"
47
48#include "G4ProductionCuts.hh"
49#include "G4EmProcessOptions.hh"
50
51#include "G4ComptonScattering.hh"
52#include "G4GammaConversion.hh"
53#include "G4PhotoElectricEffect.hh"
54
55#include "G4eMultipleScattering.hh"
56#include "G4MuMultipleScattering.hh"
57#include "G4hMultipleScattering.hh"
58
59#include "G4eIonisation.hh"
60#include "G4eBremsstrahlung.hh"
61#include "G4eplusAnnihilation.hh"
62#include "G4PAIModel.hh"
63#include "G4PAIPhotonModel.hh"
64
65#include "G4SynchrotronRadiation.hh"
66
67#include "G4MuIonisation.hh"
68#include "G4MuBremsstrahlung.hh"
69#include "G4MuPairProduction.hh"
70
71#include "G4hIonisation.hh"
72
73#include "G4Decay.hh"
74
75#include "G4VXTRenergyLoss.hh"
76#include "G4RegularXTRadiator.hh"
77#include "G4TransparentRegXTRadiator.hh"
78#include "G4GammaXTRadiator.hh"
79#include "G4StrawTubeXTRadiator.hh"
80
81#include "G4XTRGammaRadModel.hh"
82#include "G4XTRRegularRadModel.hh"
83#include "G4XTRTransparentRegRadModel.hh"
84#include "Em10XTRTransparentRegRadModel.hh"
85
86#include "Em10StepCut.hh"
87
88/////////////////////////////////////////////////////////////
89//
90//
91
92Em10PhysicsList::Em10PhysicsList(Em10DetectorConstruction* p)
93  :  G4VUserPhysicsList(),
94     MaxChargedStep(DBL_MAX),
95     theeminusStepCut(0),            theeplusStepCut(0),
96     fRadiatorCuts(0),fDetectorCuts(0),fXTRModel("transpM")
97{
98  pDet = p;
99
100  // world cuts
101
102  defaultCutValue = 1.*mm;
103  cutForGamma     = defaultCutValue;
104  cutForElectron  = defaultCutValue;
105  cutForPositron  = defaultCutValue;
106
107  // Region cuts
108
109  fGammaCut    = defaultCutValue;
110  fElectronCut = defaultCutValue;
111  fPositronCut = defaultCutValue;
112
113  SetVerboseLevel(1);
114  physicsListMessenger = new Em10PhysicsListMessenger(this);
115}
116
117Em10PhysicsList::~Em10PhysicsList()
118{
119  delete physicsListMessenger; 
120}
121
122///////////////////////////////////////////////////////////////////////////
123//
124//
125
126void Em10PhysicsList::ConstructParticle()
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();
134  ConstructLeptons();
135  ConstructMesons();
136  ConstructBarions();
137}
138
139////////////////////////////////////////////////////////////////////////////
140//
141//
142
143void Em10PhysicsList::ConstructBosons()
144{
145  // gamma
146  G4Gamma::GammaDefinition();
147}
148
149void Em10PhysicsList::ConstructLeptons()
150{
151  // leptons
152
153  G4Electron::ElectronDefinition();
154  G4Positron::PositronDefinition();
155  G4MuonPlus::MuonPlusDefinition();
156  G4MuonMinus::MuonMinusDefinition();
157
158  G4NeutrinoE::NeutrinoEDefinition();
159  G4AntiNeutrinoE::AntiNeutrinoEDefinition();
160  G4NeutrinoMu::NeutrinoMuDefinition();
161  G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
162}
163
164void Em10PhysicsList::ConstructMesons()
165{
166 //  mesons
167
168  G4PionPlus::PionPlusDefinition();
169  G4PionMinus::PionMinusDefinition();
170  G4PionZero::PionZeroDefinition();
171  G4KaonPlus::KaonPlusDefinition();
172  G4KaonMinus::KaonMinusDefinition();
173}
174
175
176void Em10PhysicsList::ConstructBarions()
177{
178//  barions
179
180  G4Proton::ProtonDefinition();
181  G4AntiProton::AntiProtonDefinition();
182}
183
184
185///////////////////////////////////////////////////////////////////////
186//
187//
188
189void Em10PhysicsList::ConstructProcess()
190{
191  AddTransportation();
192  ConstructEM();
193  ConstructGeneral();
194}
195
196/////////////////////////////////////////////////////////////////////////////
197//
198//
199
200void Em10PhysicsList::ConstructEM()
201{
202 
203  // G4cout<<"fMinElectronEnergy = "<<fMinElectronEnergy/keV<<" keV"<<G4endl;
204  // G4cout<<"fMinGammaEnergy = "<<fMinGammaEnergy/keV<<" keV"<<G4endl;
205  G4cout<<"XTR model = "<<fXTRModel<<G4endl;
206
207  const G4RegionStore* theRegionStore = G4RegionStore::GetInstance();
208  G4Region* gas = theRegionStore->GetRegion("XTRdEdxDetector");
209
210  G4VXTRenergyLoss* processXTR = 0;
211
212  if(fXTRModel == "gammaR" )         
213  {     
214    // G4GammaXTRadiator*
215    processXTR = new G4GammaXTRadiator(pDet->GetLogicalRadiator(),
216                                       100.,
217                                       100.,
218                                       pDet->GetFoilMaterial(),
219                                       pDet->GetGasMaterial(),
220                                       pDet->GetFoilThick(),
221                                       pDet->GetGasThick(),
222                                       pDet->GetFoilNumber(),
223                                       "GammaXTRadiator");
224  }
225  else if(fXTRModel == "gammaM" ) 
226  {
227    // G4XTRGammaRadModel*
228    processXTR = new G4XTRGammaRadModel(pDet->GetLogicalRadiator(),
229                                       100.,
230                                       100.,
231                                       pDet->GetFoilMaterial(),
232                                       pDet->GetGasMaterial(),
233                                       pDet->GetFoilThick(),
234                                       pDet->GetGasThick(),
235                                       pDet->GetFoilNumber(),
236                                       "GammaXTRadiator");
237  }
238  else if(fXTRModel == "strawR" ) 
239  {
240
241    // G4StrawTubeXTRadiator*
242    processXTR = new G4StrawTubeXTRadiator(pDet->GetLogicalRadiator(),
243                                         pDet->GetFoilMaterial(),
244                                         pDet->GetGasMaterial(),
245                                0.53,      // pDet->GetFoilThick(),
246                                3.14159,           // pDet->GetGasThick(),
247                                         pDet->GetAbsorberMaterial(),
248                                         true,
249                                         "strawXTRadiator");
250  }
251  else if(fXTRModel == "regR" ) 
252  {     
253    // G4RegularXTRadiator*
254    processXTR = new G4RegularXTRadiator(pDet->GetLogicalRadiator(),
255                                         pDet->GetFoilMaterial(),
256                                         pDet->GetGasMaterial(),
257                                         pDet->GetFoilThick(),
258                                         pDet->GetGasThick(),
259                                         pDet->GetFoilNumber(),
260                                         "RegularXTRadiator");     
261  }
262  else if(fXTRModel == "transpR" ) 
263  {
264    // G4TransparentRegXTRadiator*
265    processXTR = new G4TransparentRegXTRadiator(pDet->GetLogicalRadiator(),
266                                         pDet->GetFoilMaterial(),
267                                         pDet->GetGasMaterial(),
268                                         pDet->GetFoilThick(),
269                                         pDet->GetGasThick(),
270                                         pDet->GetFoilNumber(),
271                                         "RegularXTRadiator");
272  }
273  else if(fXTRModel == "regM" ) 
274  {
275    // G4XTRRegularRadModel*
276    processXTR = new G4XTRRegularRadModel(pDet->GetLogicalRadiator(),
277                                         pDet->GetFoilMaterial(),
278                                         pDet->GetGasMaterial(),
279                                         pDet->GetFoilThick(),
280                                         pDet->GetGasThick(),
281                                         pDet->GetFoilNumber(),
282                                         "RegularXTRadiator");
283       
284  }
285  else if(fXTRModel == "transpM" ) 
286  { 
287    // G4XTRTransparentRegRadModel*
288    // processXTR = new G4XTRTransparentRegRadModel(pDet->GetLogicalRadiator(),
289    processXTR = new Em10XTRTransparentRegRadModel(pDet->GetLogicalRadiator(),
290                                         pDet->GetFoilMaterial(),
291                                         pDet->GetGasMaterial(),
292                                         pDet->GetFoilThick(),
293                                         pDet->GetGasThick(),
294                                         pDet->GetFoilNumber(),
295                                         "RegularXTRadiator");
296  }     
297  else
298  {
299    G4Exception("Invalid XTR model name", "InvalidSetup",
300                 FatalException, "XTR model name is out of the name list");
301  }     
302  //  processXTR->SetCompton(true);
303  processXTR->SetVerboseLevel(1);
304
305  theParticleIterator->reset();
306
307  while( (*theParticleIterator)() )
308  {
309    G4ParticleDefinition* particle = theParticleIterator->value();
310    G4ProcessManager* pmanager = particle->GetProcessManager();
311    G4String particleName = particle->GetParticleName();
312
313    if (particleName == "gamma")
314    {
315      // Construct processes for gamma
316
317      pmanager->AddDiscreteProcess(new G4PhotoElectricEffect);
318      pmanager->AddDiscreteProcess(new G4ComptonScattering);
319      pmanager->AddDiscreteProcess(new G4GammaConversion);
320
321    }
322    else if (particleName == "e-")
323    {
324      // Construct processes for electron
325      theeminusStepCut = new Em10StepCut();
326      theeminusStepCut->SetMaxStep(MaxChargedStep) ;
327      G4eIonisation* eioni = new G4eIonisation();
328      G4PAIModel*     pai = new G4PAIModel(particle,"PAIModel");
329      eioni->AddEmModel(0,pai,pai,gas);
330
331      pmanager->AddProcess(new G4eMultipleScattering,-1,1,1);
332      pmanager->AddProcess(eioni,-1,2,2);
333      pmanager->AddProcess(new G4eBremsstrahlung,-1,3,3);
334      pmanager->AddDiscreteProcess(processXTR);
335      pmanager->AddDiscreteProcess(new G4SynchrotronRadiation);
336      pmanager->AddDiscreteProcess(theeminusStepCut);
337
338    }
339    else if (particleName == "e+")
340    {
341      // Construct processes for positron
342
343      theeplusStepCut = new Em10StepCut();
344      theeplusStepCut->SetMaxStep(MaxChargedStep) ;
345      G4eIonisation* eioni = new G4eIonisation();
346      G4PAIModel*     pai = new G4PAIModel(particle,"PAIModel");
347      eioni->AddEmModel(0,pai,pai,gas);
348
349      pmanager->AddProcess(new G4eMultipleScattering,-1,1,1);
350      pmanager->AddProcess(eioni,-1,2,2);
351      pmanager->AddProcess(new G4eBremsstrahlung,-1,3,3);
352      pmanager->AddProcess(new G4eplusAnnihilation,0,-1,4);
353      pmanager->AddDiscreteProcess(processXTR);
354      pmanager->AddDiscreteProcess(new G4SynchrotronRadiation);
355      pmanager->AddDiscreteProcess(theeplusStepCut);
356
357    }
358    else if( particleName == "mu+" ||
359             particleName == "mu-"    )
360    {
361     // Construct processes for muon+
362
363      Em10StepCut* muonStepCut = new Em10StepCut();
364      muonStepCut->SetMaxStep(MaxChargedStep) ;
365
366      G4MuIonisation* muioni = new G4MuIonisation() ;
367
368      G4PAIModel*     pai = new G4PAIModel(particle,"PAIModel");
369      muioni->AddEmModel(0,pai,pai,gas);
370
371      pmanager->AddProcess(new G4MuMultipleScattering(),-1,1,1);
372      pmanager->AddProcess(muioni,-1,2,2);
373      pmanager->AddProcess(new G4MuBremsstrahlung(),-1,3,3);
374      pmanager->AddProcess(new G4MuPairProduction(),-1,4,4);
375      pmanager->AddProcess( muonStepCut,-1,-1,5);
376
377    }
378    else if (
379                particleName == "proton"
380               || particleName == "antiproton"
381               || particleName == "pi+"
382               || particleName == "pi-"
383               || particleName == "kaon+"
384               || particleName == "kaon-"
385              )
386    {
387      Em10StepCut* thehadronStepCut = new Em10StepCut();
388      thehadronStepCut->SetMaxStep(MaxChargedStep) ;
389
390      G4hIonisation* thehIonisation = new G4hIonisation();
391      G4PAIModel*     pai = new G4PAIModel(particle,"PAIModel");
392      thehIonisation->AddEmModel(0,pai,pai,gas);
393
394      pmanager->AddProcess(new G4hMultipleScattering,-1,1,1);
395      pmanager->AddProcess(thehIonisation,-1,2,2);
396      pmanager->AddProcess( thehadronStepCut,-1,-1,3);
397
398    }
399  }
400  G4EmProcessOptions opt;
401  opt.SetApplyCuts(true);
402}
403
404void Em10PhysicsList::ConstructGeneral()
405{
406  // Add Decay Process
407
408  G4Decay* theDecayProcess = new G4Decay();
409  theParticleIterator->reset();
410
411  while( (*theParticleIterator)() )
412  {
413    G4ParticleDefinition* particle = theParticleIterator->value();
414    G4ProcessManager* pmanager = particle->GetProcessManager();
415
416    if (theDecayProcess->IsApplicable(*particle)) 
417    { 
418      pmanager ->AddProcess(theDecayProcess);
419
420      // set ordering for PostStepDoIt and AtRestDoIt
421
422      pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
423      pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
424    }
425  }
426}
427
428/////////////////////////////////////////////////////////////////////////////
429
430void Em10PhysicsList::SetCuts()
431{
432  // set cut values for gamma at first and for e- second and next for e+,
433  // because some processes for e+/e- need cut values for gamma
434 
435  SetCutValue(cutForGamma, "gamma", "DefaultRegionForTheWorld");
436  SetCutValue(cutForElectron, "e-", "DefaultRegionForTheWorld");
437  SetCutValue(cutForPositron, "e+", "DefaultRegionForTheWorld");
438
439  if (verboseLevel > 0)
440  {
441    G4cout << "Em10PhysicsList::SetCuts:";
442    G4cout << "CutLength for e-, e+ and gamma is: " 
443           << G4BestUnit(defaultCutValue,"Length") << G4endl;
444  }
445 
446  if( !fRadiatorCuts ) SetRadiatorCuts();
447
448  G4Region* region = (G4RegionStore::GetInstance())->GetRegion("XTRradiator");
449  region->SetProductionCuts(fRadiatorCuts);
450  G4cout << "Radiator cuts are set" << G4endl;
451
452  if( !fDetectorCuts ) SetDetectorCuts();
453  region = (G4RegionStore::GetInstance())->GetRegion("XTRdEdxDetector");
454  region->SetProductionCuts(fDetectorCuts);
455  G4cout << "Detector cuts are set" << G4endl;
456
457  if (verboseLevel > 1)     DumpCutValuesTable();
458}
459
460///////////////////////////////////////////////////////////////////////////
461
462void Em10PhysicsList::SetGammaCut(G4double val)
463{
464  cutForGamma = val;
465}
466
467///////////////////////////////////////////////////////////////////////////
468
469void Em10PhysicsList::SetElectronCut(G4double val)
470{
471  cutForElectron = val;
472}
473
474////////////////////////////////////////////////////////////////////////////
475
476void Em10PhysicsList::SetMaxStep(G4double step)
477{
478  MaxChargedStep = step ;
479  G4cout << " MaxChargedStep=" << MaxChargedStep << G4endl;
480  G4cout << G4endl;
481}
482
483/////////////////////////////////////////////////////
484
485void Em10PhysicsList::SetRadiatorCuts()
486{
487  if( !fRadiatorCuts ) fRadiatorCuts = new G4ProductionCuts();
488
489  fRadiatorCuts->SetProductionCut(fGammaCut, idxG4GammaCut);
490  fRadiatorCuts->SetProductionCut(fElectronCut, idxG4ElectronCut);
491  fRadiatorCuts->SetProductionCut(fPositronCut, idxG4PositronCut);
492
493  G4cout<<"Radiator gamma cut    = "<<fGammaCut/mm<<" mm"<<G4endl;
494  G4cout<<"Radiator electron cut = "<<fElectronCut/mm<<" mm"<<G4endl;
495  G4cout<<"Radiator positron cut = "<<fPositronCut/mm<<" mm"<<G4endl;
496}
497
498/////////////////////////////////////////////////////////////
499
500void Em10PhysicsList::SetDetectorCuts()
501{
502  if( !fDetectorCuts ) fDetectorCuts = new G4ProductionCuts();
503
504  fDetectorCuts->SetProductionCut(fGammaCut, idxG4GammaCut);
505  fDetectorCuts->SetProductionCut(fElectronCut, idxG4ElectronCut);
506  fDetectorCuts->SetProductionCut(fPositronCut, idxG4PositronCut);
507
508  G4cout<<"Detector gamma cut    = "<<fGammaCut/mm<<" mm"<<G4endl;
509  G4cout<<"Detector electron cut = "<<fElectronCut/mm<<" mm"<<G4endl;
510  G4cout<<"Detector positron cut = "<<fPositronCut/mm<<" mm"<<G4endl;
511
512}
Note: See TracBrowser for help on using the repository browser.