Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

testG4EnergyLossTables.cc

Last change on this file was 1347, checked in by garnier, 14 years ago
geant4 tag 9.4
File size: 23.9 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: testG4EnergyLossTables.cc,v 1.8 2010/08/17 17:44:18 vnivanch Exp $
28	// GEANT4 tag $Name: geant4-09-04-ref-00 $
29	//
30	//-------------------------------------------------------------------
31	#include "G4ios.hh"
32	#include <fstream>
33	#include <iomanip>
34	#include "globals.hh"
35	#include "G4Timer.hh"
36	#include "G4MultipleScattering.hh"
37	#include "G4DynamicParticle.hh"
38	#include "G4Element.hh"
39	#include "G4Material.hh"
40	#include "G4GRSVolume.hh"
41	#include "G4PVPlacement.hh"
42	#include "G4LogicalVolume.hh"
43	#include "G4Box.hh"
44	#include "G4ProcessManager.hh"
45	#include "G4Step.hh"
46	#include "G4StepPoint.hh"
47	#include "G4Track.hh"
48	#include "G4Proton.hh"
49	#include "G4AntiProton.hh"
50	#include "G4PionPlus.hh"
51	#include "G4PionMinus.hh"
52	#include "G4KaonPlus.hh"
53	#include "G4KaonMinus.hh"
54	#include "G4MuonPlus.hh"
55	#include "G4MuonMinus.hh"
56	#include "G4Electron.hh"
57	#include "G4Positron.hh"
58	#include "G4Gamma.hh"
59	#include "G4eEnergyLoss.hh"
60	#include "G4eIonisation.hh"
61	#include "G4eBremsstrahlung.hh"
62	#include "G4MuEnergyLoss.hh"
63	#include "G4MuIonisation.hh"
64	#include "G4MuBremsstrahlung.hh"
65	#include "G4MuPairProduction.hh"
66	#include "G4hEnergyLoss.hh"
67	#include "G4hIonisation.hh"
68	#include "G4GPILSelection.hh"
69
70	#include "G4EnergyLossTables.hh"
71
72	// test for G4EnergyLossTables, based on tunemsc
73	// P. Urban
74
75	G4VPhysicalVolume* BuildVolume(G4Material* matworld)
76	// it builds a simple box filled with material matword .......
77	{
78	G4Box myWorldBox= new G4Box ("WBox",10000.cm,10000.cm,10000.cm);
79
80	G4LogicalVolume *myWorldLog = new G4LogicalVolume(myWorldBox,matworld,
81	"WLog",0,0,0) ;
82
83	G4PVPlacement *myWorldPhys = new G4PVPlacement(0,G4ThreeVector(),
84	"WPhys",
85	myWorldLog,
86	0,false,0) ;
87
88
89	return myWorldPhys ;
90
91	}
92
93
94	int main()
95	{
96	//-------- set output format-------
97	G4cout.setf( std::ios::scientific, std::ios::floatfield );
98	//---write results to the file msc.out-----
99	std::ofstream outFile("msc.out", std::ios::out ) ;
100	outFile.setf( std::ios::scientific, std::ios::floatfield );
101
102	//--------- Material definition ---------
103	G4Timer theTimer ;
104	G4double a, z, ez, density ,temperature,pressure;
105	G4State state ;
106	G4String name, symbol;
107	G4int nel;
108
109	a = 9.012*g/mole;
110	density = 1.848*g/cm3;
111	G4Material* Be = new G4Material(name="Beryllium", z=4. , a, density);
112
113	a = 12.011*g/mole;
114	density = 2.220*g/cm3;
115	G4Material* C = new G4Material(name="Carbon", z=6. , a, density);
116
117	a = 26.98*g/mole;
118	density = 2.7*g/cm3;
119	G4Material* Al = new G4Material(name="Aluminium", z=13., a, density);
120
121	a = 28.09*g/mole;
122	density = 2.33*g/cm3;
123	G4Material* Si = new G4Material(name="Silicon", z=14., a, density);
124
125	a = 55.85*g/mole;
126	density = 7.87*g/cm3;
127	G4Material* Fe = new G4Material(name="Iron", z=26., a, density);
128
129	a = 63.540*g/mole;
130	density = 8.960*g/cm3;
131	G4Material* Cu = new G4Material(name="Copper", z=29., a, density);
132
133	a = 196.97*g/mole;
134	density = 19.32*g/cm3;
135	G4Material* Au = new G4Material(name="Gold", z=79., a, density);
136
137	a = 207.19*g/mole;
138	density = 11.35*g/cm3;
139	G4Material* Pb = new G4Material(name="Lead", z=82., a, density);
140
141	a = 238.03*g/mole;
142	density = 18.7*g/cm3;
143	G4Material* U = new G4Material(name="Uranium", z=92., a, density);
144
145	a = 14.01*g/mole;
146	G4Element* elN = new G4Element(name="Nitrogen", symbol="N", ez=7., a);
147	a = 16.00*g/mole;
148	G4Element* elO = new G4Element(name="Oxigen", symbol="O", ez=8., a);
149	density = 1.29e-03*g/cm3;
150	state = kStateGas ;
151	temperature = 273.*kelvin ;
152	pressure = 1.*atmosphere ;
153	G4Material* Air = new G4Material(name="Air", density, nel=2 ,
154	state ,temperature , pressure ) ;
155	Air->AddElement(elN, .7);
156	Air->AddElement(elO, .3);
157
158	a = 0. ;
159	density = 0. ;
160	G4Material* Vac= new G4Material(name="Vacuum",z=0., a, density,kVacuum);
161
162	static const G4MaterialTable* theMaterialTable ;
163	G4Material* apttoMaterial ;
164	G4String MaterialName ;
165
166	//--------- Particle definition ---------
167
168	G4ParticleDefinition* theGamma = G4Gamma::GammaDefinition();
169
170	G4ParticleDefinition* theElectron = G4Electron::ElectronDefinition();
171	G4ParticleDefinition* thePositron = G4Positron::PositronDefinition();
172
173	G4ParticleDefinition* theMuonPlus = G4MuonPlus::MuonPlusDefinition();
174	G4ParticleDefinition* theMuonMinus = G4MuonMinus::MuonMinusDefinition();
175
176	G4ParticleDefinition* theProton = G4Proton::ProtonDefinition();
177	G4ParticleDefinition* theAntiProton = G4AntiProton::AntiProtonDefinition();
178	G4ParticleDefinition* thePionPlus = G4PionPlus::PionPlusDefinition();
179	G4ParticleDefinition* thePionMinus = G4PionMinus::PionMinusDefinition();
180	G4ParticleDefinition* theKaonPlus = G4KaonPlus::KaonPlusDefinition();
181	G4ParticleDefinition* theKaonMinus = G4KaonMinus::KaonMinusDefinition();
182
183	//--------- Process definition ---------
184	G4eIonisation theElectronIonisation,thePositronIonisation;
185	G4eBremsstrahlung theElectronBremsstrahlung,thePositronBremsstrahlung;
186	//..........e-/e+....................................................
187	G4MultipleScattering theElectronMultipleScattering,thePositronMultipleScattering ;
188
189	G4ProcessManager* theElectronProcessManager = theElectron->GetProcessManager();
190	G4ProcessManager* thePositronProcessManager = thePositron->GetProcessManager();
191
192	theElectronProcessManager->AddProcess(&theElectronMultipleScattering,-1,0,0) ;
193	thePositronProcessManager->AddProcess(&thePositronMultipleScattering,-1,0,0) ;
194	theElectronProcessManager->AddProcess(&theElectronIonisation,-1,1,1) ;
195	theElectronProcessManager->AddProcess(&theElectronBremsstrahlung,-1,-1,2) ;
196	thePositronProcessManager->AddProcess(&thePositronIonisation,-1,1,1) ;
197	thePositronProcessManager->AddProcess(&thePositronBremsstrahlung,-1,-1,2) ;
198
199	//..........mu+/mu-................................................
200	G4MuIonisation theMuonPlusIonisation,theMuonMinusIonisation ;
201	G4MuBremsstrahlung theMuonPlusBremsstrahlung,theMuonMinusBremsstrahlung ;
202	G4MuPairProduction theMuonPlusPairProduction,theMuonMinusPairProduction ;
203	G4MultipleScattering theMuonPlusMultipleScattering,theMuonMinusMultipleScattering;
204
205	G4ProcessManager* theMuonPlusProcessManager = theMuonPlus->GetProcessManager();
206	theMuonPlusProcessManager->AddProcess(&theMuonPlusMultipleScattering,-1,0,0) ;
207	theMuonPlusProcessManager->AddProcess(&theMuonPlusIonisation,-1,1,1);
208	theMuonPlusProcessManager->AddProcess(&theMuonPlusBremsstrahlung,-1,-1,2);
209	theMuonPlusProcessManager->AddProcess(&theMuonPlusPairProduction,-1,-1,3);
210
211	G4ProcessManager* theMuonMinusProcessManager = theMuonMinus->GetProcessManager();
212	theMuonMinusProcessManager->AddProcess(&theMuonMinusMultipleScattering,-1,0,0) ;
213	theMuonMinusProcessManager->AddProcess(&theMuonMinusIonisation,-1,1,1);
214	theMuonMinusProcessManager->AddProcess(&theMuonMinusBremsstrahlung,-1,-1,2);
215	theMuonMinusProcessManager->AddProcess(&theMuonMinusPairProduction,-1,-1,3);
216
217	//------multiple instantiation of G4MultipleScattering for hadrons-----------------
218	G4MultipleScattering theProtonMultipleScattering,theAntiProtonMultipleScattering;
219	G4hIonisation theProtonIonisation,theAntiProtonIonisation ;
220	G4MultipleScattering thePionPlusMultipleScattering,thePionMinusMultipleScattering;
221	G4hIonisation thePionPlusIonisation,thePionMinusIonisation;
222	G4MultipleScattering theKaonPlusMultipleScattering,theKaonMinusMultipleScattering;
223	G4hIonisation theKaonPlusIonisation,theKaonMinusIonisation;
224
225	G4ProcessManager* theProtonProcessManager = theProton->GetProcessManager();
226	theProtonProcessManager->AddProcess(&theProtonMultipleScattering,-1,0,0) ;
227	theProtonProcessManager->AddProcess(&theProtonIonisation,-1,1,1);
228
229	G4ProcessManager* thePionPlusProcessManager = thePionPlus->GetProcessManager();
230	thePionPlusProcessManager->AddProcess(&thePionPlusMultipleScattering,-1,0,0) ;
231	thePionPlusProcessManager->AddProcess(&thePionPlusIonisation,-1,1,1);
232
233	G4ProcessManager* theKaonPlusProcessManager = theKaonPlus->GetProcessManager();
234	theKaonPlusProcessManager->AddProcess(&theKaonPlusMultipleScattering,-1,0,0) ;
235	theKaonPlusProcessManager->AddProcess(&theKaonPlusIonisation,-1,1,1);
236
237	G4ProcessManager* theAntiProtonProcessManager = theAntiProton->GetProcessManager();
238	theAntiProtonProcessManager->AddProcess(&theAntiProtonMultipleScattering,-1,0,0) ;
239	theAntiProtonProcessManager->AddProcess(&theAntiProtonIonisation,-1,1,1);
240
241	G4ProcessManager* thePionMinusProcessManager = thePionMinus->GetProcessManager();
242	thePionMinusProcessManager->AddProcess(&thePionMinusMultipleScattering,-1,0,0) ;
243	thePionMinusProcessManager->AddProcess(&thePionMinusIonisation,-1,1,1);
244
245	G4ProcessManager* theKaonMinusProcessManager = theKaonMinus->GetProcessManager();
246	theKaonMinusProcessManager->AddProcess(&theKaonMinusMultipleScattering,-1,0,0) ;
247	theKaonMinusProcessManager->AddProcess(&theKaonMinusIonisation,-1,1,1);
248
249	G4GPILSelection selection;
250	G4ParticleWithCuts* theParticle ;
251	G4MultipleScattering* theParticleMultipleScattering;
252	G4ProcessManager* theParticleProcessManager ;
253	G4String confirm ;
254	G4int i1e ;
255	G4double theta1e,th1,th2 ;
256	G4double theta1edata,errth,lambdadata,lambdadatamin,lambdadatamax ;
257
258	// loop until no particle is selected
259	while (1) {
260
261	G4cout << "Do you want the electron as particle (yes/no)?" << std::flush;
262	G4cin >> confirm ;
263	if(confirm == "yes")
264	{
265	theParticle = theElectron ;
266	theParticleMultipleScattering=&theElectronMultipleScattering;
267	theParticleProcessManager=theElectronProcessManager;
268	outFile << " ----------particle = electron -------------" << G4endl;
269	}
270	else
271	{
272	G4cout << "Do you want the positron as particle (yes/no)?" << std::flush;
273	G4cin >> confirm ;
274	if(confirm == "yes")
275	{
276	theParticle = thePositron ;
277	theParticleMultipleScattering=&thePositronMultipleScattering;
278	theParticleProcessManager=thePositronProcessManager;
279	outFile << " ----------particle = positron -------------" << G4endl;
280	}
281	else
282	{
283	G4cout << "Do you want the mu+ as particle (yes/no)?" << std::flush;
284	G4cin >> confirm ;
285	if(confirm == "yes")
286	{
287	theParticle = theMuonPlus ;
288	theParticleMultipleScattering=&theMuonPlusMultipleScattering;
289	theParticleProcessManager=theMuonPlusProcessManager;
290	outFile << " --------particle = mu+ -------------" << G4endl;
291	}
292	else
293	{
294	G4cout << "Do you want the mu- as particle (yes/no)?" << std::flush;
295	G4cin >> confirm ;
296	if(confirm == "yes")
297	{
298	theParticle = theMuonMinus ;
299	theParticleMultipleScattering=&theMuonMinusMultipleScattering;
300	theParticleProcessManager=theMuonMinusProcessManager;
301	outFile << " --------particle = mu- -------------" << G4endl;
302	}
303	else
304	{
305	G4cout << " Do you want the proton as particle (yes/no)? " << std::flush;
306	G4cin >> confirm ;
307	if(confirm == "yes")
308	{
309	theParticle = theProton;
310	theParticleMultipleScattering=&theProtonMultipleScattering;
311	theParticleProcessManager=theProtonProcessManager;
312	outFile << " ---------- particle = proton ----------------" << G4endl;
313	}
314	else
315	{
316	G4cout << " Do you want the antiproton as particle (yes/no)? " << std::flush;
317	G4cin >> confirm ;
318	if(confirm == "yes")
319	{
320	theParticle = theAntiProton;
321	theParticleMultipleScattering=&theAntiProtonMultipleScattering;
322	theParticleProcessManager=theAntiProtonProcessManager;
323	outFile << " ---------- particle = antiproton ----------------" << G4endl;
324	}
325	else
326	{
327	G4cout << " Do you want the pi+ as particle (yes/no)? " << std::flush;
328	G4cin >> confirm ;
329	if(confirm == "yes")
330	{
331	theParticle = thePionPlus;
332	theParticleMultipleScattering=&thePionPlusMultipleScattering;
333	theParticleProcessManager=thePionPlusProcessManager;
334	outFile << " ---------- particle = pi+ ----------------" << G4endl;
335	}
336	else
337	{
338	G4cout << " Do you want the pi- as particle (yes/no)? " << std::flush;
339	G4cin >> confirm ;
340	if(confirm == "yes")
341	{
342	theParticle = thePionMinus;
343	theParticleMultipleScattering=&thePionMinusMultipleScattering;
344	theParticleProcessManager=thePionMinusProcessManager;
345	outFile << " ---------- particle = pi- ----------------" << G4endl;
346	}
347	else
348	{
349	G4cout << " Do you want the K+ as particle (yes/no)? " << std::flush;
350	G4cin >> confirm ;
351	if(confirm == "yes")
352	{
353	theParticle = theKaonPlus;
354	theParticleMultipleScattering=&theKaonPlusMultipleScattering;
355	theParticleProcessManager=theKaonPlusProcessManager;
356	outFile << " ---------- particle = K+ ----------------" << G4endl;
357	}
358	else
359	{
360	G4cout << " Do you want the K- as particle (yes/no)? " << std::flush;
361	G4cin >> confirm ;
362	if(confirm == "yes")
363	{
364	theParticle = theKaonMinus;
365	theParticleMultipleScattering=&theKaonMinusMultipleScattering;
366	theParticleProcessManager=theKaonMinusProcessManager;
367	outFile << " ---------- particle = K- ----------------" << G4endl;
368	}
369	else
370	{
371	G4cout << " There is no other particle in the test." << G4endl;
372	return EXIT_SUCCESS;
373	}
374	}
375	}
376	}
377	}
378	}
379	}
380	}
381	}
382	}
383
384	G4ForceCondition cond;
385	G4ForceCondition* condition=&cond ;
386
387	G4double* ElectronKineticEnergyCuts ;
388	G4double* PositronKineticEnergyCuts ;
389	G4double* GammaKineticEnergyCuts ;
390	G4double* ParticleKineticEnergyCuts ;
391
392	theMaterialTable = G4Material::GetMaterialTable() ;
393
394	G4double cutinrange ;
395
396
397	G4cout << "give cuts in range" << G4endl ;
398
399	G4cout << "cut for GAMMA in mm =" ;
400	G4cin >> cutinrange ;
401	theGamma->SetCuts(cutinrange) ;
402
403	GammaKineticEnergyCuts = theGamma->GetCutsInEnergy() ;
404
405	G4cout << "cut for ELECTRON in mm =" ;
406	G4cin >> cutinrange ;
407	theElectron->SetCuts(cutinrange) ;
408
409	ElectronKineticEnergyCuts = theElectron->GetCutsInEnergy() ;
410
411	G4cout << "cut for POSITRON in mm =" ;
412	G4cin >> cutinrange ;
413	thePositron->SetCuts(cutinrange) ;
414
415	PositronKineticEnergyCuts = thePositron->GetCutsInEnergy() ;
416
417	//****************************************************
418	// setcut for the selected particle (if it is not e-/e+)
419	if((theParticle != theElectron) && (theParticle != thePositron))
420	{
421	G4cout << "cut for the selected particle in mm =" ;
422	G4cin >> cutinrange ;
423	theParticle->SetCuts(cutinrange) ;
424
425	ParticleKineticEnergyCuts = theParticle->GetEnergyCuts() ;
426	}
427
428	G4double energy, momentum, mass;
429	G4ProcessVector* palongget ;
430	G4ProcessVector* palongdo ;
431	G4ProcessVector* ppostget ;
432	G4ProcessVector* ppostdo ;
433
434	mass = theParticle->GetPDGMass() ;
435	energy = 1.*GeV + mass ;
436	momentum=std::sqrt(energyenergy-massmass) ;
437	G4ParticleMomentum theMomentum(momentum,0.,0.);
438	G4double pModule = theMomentum.mag();
439	G4DynamicParticle aParticle(theParticle,energy,theMomentum);
440	aParticle.SetKineticEnergy(energy-mass);
441
442	outFile << " " << G4endl;
443	outFile << " M S C test **********************************************" << G4endl ;
444	outFile << " " << G4endl;
445	palongget = aParticle.GetParticleDefinition()->GetProcessManager()
446	->GetAlongStepProcessVector(0);
447	ppostget = aParticle.GetParticleDefinition()->GetProcessManager()
448	->GetPostStepProcessVector(0);
449	palongdo = aParticle.GetParticleDefinition()->GetProcessManager()
450	->GetAlongStepProcessVector(1);
451	ppostdo = aParticle.GetParticleDefinition()->GetProcessManager()
452	->GetPostStepProcessVector(1);
453
454	//---------------------------------- Physics --------------------------------
455
456	G4int itry=1, Ntry=1, Nstart, ir;
457	G4double r ;
458
459	//**************************************************************************
460	const G4int Nbin=100 ;
461	G4double TkinMeV[Nbin] =
462	{0.0002, 0.0005, 0.001,0.0015,0.002,0.003,0.004,0.005,0.006,0.008,
463	0.01,0.015,0.02,0.03,0.04,0.05,0.06,0.08,
464	0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.8,
465	1.,1.5,2.,3.,4.,5.,6.,8.,
466	10.,15.,20.,30.,40.,50.,60.,80.,
467	100.,150.,200.,300.,400.,500.,600.,800.,
468	1.0e3,1.5e3,2.0e3,3.0e3,4.0e3,5.0e3,6.0e3,8.0e3,
469	1.0e4,1.5e4,2.0e4,3.0e4,4.0e4,5.0e4,6.0e4,8.0e4,
470	1.0e5,1.5e5,2.0e5,3.0e5,4.0e5,5.0e5,6.0e5,8.0e5,
471	1.0e6,1.5e6,2.0e6,3.0e6,4.0e6,5.0e6,6.0e6,8.0e6,
472	1.0e7,1.5e7,2.0e7,3.0e7,4.0e7,5.0e7,6.0e7,8.0e7,
473	1.0e8,1.5e8,2.0e8,3.0e8,4.0e8,5.0e8,6.0e8,8.0e8,
474	1.0e9, 5.0e9} ;
475
476	G4int J=-1 ;
477
478	G4double lambda,trueStep,geomStep,stepLimit,
479	previousStepSize,currentMinimumStep ;
480	G4ParticleChange* aParticleChange ;
481
482	// material loop
483	for (J=0; J<theMaterialTable->length(); J++) {
484
485	apttoMaterial = (*theMaterialTable)[ J ] ;
486	MaterialName = apttoMaterial->GetName() ;
487	G4cout << G4endl << "material=" << MaterialName << G4endl << G4endl;
488
489	//---------- Volume definition ---------------------
490
491	G4VPhysicalVolume* myVolume ;
492
493	myVolume = BuildVolume(apttoMaterial) ;
494
495	//--------- track and Step definition (for this test ONLY!)------------
496	G4ThreeVector aPosition(0.,0.,0.);
497	const G4ThreeVector aDirection(0.,0.,1.) ;
498	const G4ThreeVector transl(0.,0.,0.) ;
499
500	G4double aTime = 0. ;
501
502	G4Track* tracke = new G4Track(&aParticle,aTime,aPosition) ;
503	G4Track& trackele = (*tracke) ;
504	G4GRSVolume* touche = new G4GRSVolume(myVolume,NULL,transl);
505	(*tracke).SetTouchable(touche);
506
507	(*tracke).SetMomentumDirection(aDirection) ;
508
509	G4Step* Step = new G4Step() ;
510	//******************************************
511	tracke->SetStep(Step);
512	(*Step).InitializeStep(tracke) ;
513
514	const G4Step& Step = (*Step) ;
515
516	G4double CurrentSafety=1000.;
517	G4double& currentSafety = CurrentSafety ;
518
519	G4StepPoint* aPoint = new G4StepPoint();
520	(*aPoint).SetPosition(aPosition) ;
521	G4double safety = 10000.*cm ;
522	(*aPoint).SetSafety(safety) ;
523
524	(*Step).SetPostStepPoint(aPoint) ;
525
526	G4int ib,ev ;
527	G4double xc,yc,zc,later ;
528	//**************************************************************************
529
530	//---------------------------------------------------------------------
531	G4double distr[100]=
532	{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
533	0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
534	0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
535	0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
536	0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
537	0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
538	0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
539	0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
540	0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
541	0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
542	G4double distrx[100]=
543	{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
544	0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
545	0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
546	0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
547	0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
548	0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
549	0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
550	0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
551	0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
552	0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
553
554	G4int nev[100]=
555	{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
556	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
557	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
558	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
559	G4int nevx[100]=
560	{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
561	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
562	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
563	0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
564
565	G4double thetamean,thetamean2,dthetamean,errdistr ;
566
567	outFile << G4endl ;
568	outFile << " " << MaterialName << " PostStepDoIt (scattering) test " << G4endl ;
569	outFile << " +++++++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
570	outFile << G4endl ;
571
572	G4double wg[100] ;
573
574	G4double over,overx ;
575	G4double fwg, costheta,theta ;
576	G4int ibin,iw ;
577	G4double TMeV ;
578	G4double dthetadeg,dtheta ;
579	G4int flagdeg ;
580
581	//*************************************************************************************
582	// example: how to use the GetRange() fuctions............
583
584	for (G4int i=0; i<Nbin; i++) {
585	G4double TMeV= TkinMeV[i]*MeV;
586
587	G4double rangestart;
588	if( theParticle == theElectron )
589	rangestart = theElectronIonisation.OldGetRange(theElectron,TMeV,apttoMaterial) ;
590	if( theParticle == thePositron )
591	rangestart = thePositronIonisation.OldGetRange(thePositron,TMeV,apttoMaterial) ;
592	if( theParticle == theMuonPlus )
593	rangestart = theMuonPlusIonisation.OldGetRange(theMuonPlus,TMeV,apttoMaterial) ;
594	if( theParticle == theMuonMinus )
595	rangestart = theMuonMinusIonisation.OldGetRange(theMuonMinus,TMeV,apttoMaterial) ;
596	if( theParticle == theProton )
597	rangestart = theProtonIonisation.OldGetRange(theProton,TMeV,apttoMaterial) ;
598	if( theParticle == theAntiProton )
599	rangestart = theAntiProtonIonisation.OldGetRange(theAntiProton,TMeV,apttoMaterial) ;
600	if( theParticle == thePionPlus )
601	rangestart = thePionPlusIonisation.OldGetRange(thePionPlus,TMeV,apttoMaterial) ;
602	if( theParticle == thePionMinus )
603	rangestart = thePionMinusIonisation.OldGetRange(thePionMinus,TMeV,apttoMaterial) ;
604	if( theParticle == theKaonPlus )
605	rangestart = theKaonPlusIonisation.OldGetRange(theKaonPlus,TMeV,apttoMaterial) ;
606	if( theParticle == theKaonMinus )
607	rangestart = theKaonMinusIonisation.OldGetRange(theKaonMinus,TMeV,apttoMaterial) ;
608	G4cout << " energy (MeV) = " << TMeV << G4endl;
609	G4cout << " range (OldGetRange) = " << rangestart << G4endl ;
610
611	G4double range2= G4EnergyLossTables::GetRange(theParticle, TMeV, apttoMaterial);
612
613	G4cout << " range (G4EnergyLossTables::GetRange) = " << range2 << G4endl ;
614	}
615
616	} // of material loop
617
618	} // of largest loop
619
620	return EXIT_SUCCESS;
621	}

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

Context Navigation

source: trunk/source/processes/electromagnetic/utils/test/testG4EnergyLossTables.cc

Download in other formats: