Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: