Context Navigation

muEnergyLossTest.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago
geant4 tag 9.4
File size: 39.4 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: muEnergyLossTest.cc,v 1.7 2006/06/29 19:49:56 gunter Exp $
28	// GEANT4 tag $Name: geant4-09-04-ref-00 $
29	//
30	//-----------------------------------------------------------------
31	// new testprogram for testing the mu processes:
32	// G4MuEnergyLoss
33	// G4MuIonisation
34	// G4MuBremstrahlung
35	// G4MuPairProduction
36	// G4MuNuclearInteraction
37	//-----------------------------------------------------------------
38	// created by L. Urban , May 1998
39	//---------------------------------------------------------------
40
41	#include "G4ios.hh"
42	#include <fstream>
43	#include <iomanip>
44	#include "globals.hh"
45	#include "G4Timer.hh"
46	#include "G4MuEnergyLoss.hh"
47	#include "G4MuIonisation.hh"
48	#include "G4MuBremsstrahlung.hh"
49	#include "G4MuPairProduction.hh"
50	#include "G4MuNuclearInteraction.hh"
51	#include "G4DynamicParticle.hh"
52	#include "G4Element.hh"
53	#include "G4Material.hh"
54	#include "G4PVPlacement.hh"
55	#include "G4LogicalVolume.hh"
56	#include "G4GRSVolume.hh"
57	#include "G4Box.hh"
58	#include "G4ProcessManager.hh"
59	#include "G4Step.hh"
60	#include "G4StepPoint.hh"
61	#include "G4Track.hh"
62	#include "G4Gamma.hh"
63	#include "G4Electron.hh"
64	#include "G4Positron.hh"
65	#include "G4MuonPlus.hh"
66	#include "G4MuonMinus.hh"
67	#include "G4GPILSelection.hh"
68
69
70	G4VPhysicalVolume* BuildVolume(G4Material* matworld)
71	// it builds a simple box filled with material matword .......
72	{
73	G4Box myWorldBox= new G4Box ("WBox",10000.cm,10000.cm,10000.cm);
74
75	G4LogicalVolume *myWorldLog = new G4LogicalVolume(myWorldBox,matworld,
76	"WLog",0,0,0) ;
77
78	G4PVPlacement *myWorldPhys = new G4PVPlacement(0,G4ThreeVector(),
79	"WPhys",
80	myWorldLog,
81	0,false,0) ;
82
83
84	return myWorldPhys ;
85
86	}
87
88
89	int main()
90	{
91	//-------- set output format-------
92	G4cout.setf( std::ios::scientific, std::ios::floatfield );
93
94	G4int nrandom;
95	G4double ran ;
96
97	G4cout << "Give the number of random numbers you want to Generate at start!" << G4endl;
98	G4cin >> nrandom ;
99	if( nrandom>0)
100	{
101	for (G4int ir=0; ir<nrandom; ir++)
102	{
103	ran=G4UniformRand() ;
104	}
105	}
106	//--------- Material definition ---------
107	G4double a, z, ez, density ,temperature,pressure;
108	G4State state ;
109	G4String name, symbol;
110	G4int nel;
111
112	a = 9.012*g/mole;
113	density = 1.848*g/cm3;
114	G4Material* Be = new G4Material(name="Beryllium", z=4. , a, density);
115
116	a = 26.98*g/mole;
117	density = 2.7*g/cm3;
118	G4Material* Al = new G4Material(name="Aluminium", z=13., a, density);
119
120	a = 28.09*g/mole;
121	density = 2.33*g/cm3;
122	G4Material* Si = new G4Material(name="Silicon", z=14., a, density);
123
124	G4Element* elH = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
125
126	a = 14.01*g/mole;
127	G4Element* elN = new G4Element(name="Nitrogen", symbol="N", ez=7., a);
128	a = 16.00*g/mole;
129	G4Element* elO = new G4Element(name="Oxigen", symbol="O", ez=8., a);
130	density = 1.29e-03*g/cm3;
131	state = kStateGas ;
132	temperature = 273.*kelvin ;
133	pressure = 1.*atmosphere ;
134	G4Material* Air = new G4Material(name="Air", density, nel=2 ,
135	state ,temperature , pressure ) ;
136	Air->AddElement(elN, .7);
137	Air->AddElement(elO, .3);
138
139	G4Material* H2O = new G4Material ("Water" , 1.*g/cm3, 2);
140	H2O->AddElement(elH,2);
141	H2O->AddElement(elO,1);
142
143	a = 55.85*g/mole;
144	density = 7.87*g/cm3;
145	G4Material* Fe = new G4Material(name="Iron", z=26., a, density);
146
147	a = 196.97*g/mole;
148	density = 19.32*g/cm3;
149	G4Material* Au = new G4Material(name="Gold", z=79., a, density);
150
151	a = 207.19*g/mole;
152	density = 11.35*g/cm3;
153	G4Material* Pb = new G4Material(name="Lead", z=82., a, density);
154
155	a = 238.03*g/mole;
156	density = 18.95*g/cm3;
157	G4Material* U = new G4Material(name="Uranium", z=92., a, density);
158
159	//VacuumOnEarth (air with a very low density , like in the beam pipe of an accelerator)
160	density = 1.0e-10*g/cm3;
161	G4Material* VacuumOnEarth = new G4Material(name="VacuumOnEarth", density, nel=2);
162	VacuumOnEarth->AddElement(elN, .7);
163	VacuumOnEarth->AddElement(elO, .3);
164
165	//VacuumInSpace (H , density is extremely small)
166	a = 1.01*g/mole;
167	density = 1.0e-50*g/cm3;
168	G4Material* VacuumInSpace = new G4Material(name="VacuumInSpace",
169	z=1., a, density) ;
170
171	const G4MaterialTable* theMaterialTable ;
172	G4Material* apttoMaterial ;
173	G4String MaterialName ;
174	G4Timer theTimer ;
175	G4double mloss,sloss,dEdxdelta,dEdxbrems ;
176	//--------- Particle definition ---------
177
178	G4ParticleDefinition* theGamma = G4Gamma::GammaDefinition();
179	G4ParticleDefinition* theElectron = G4Electron::ElectronDefinition();
180	G4ParticleDefinition* thePositron = G4Positron::PositronDefinition();
181	G4ParticleDefinition* theMuonPlus = G4MuonPlus::MuonPlusDefinition();
182	G4ParticleDefinition* theMuonMinus = G4MuonMinus::MuonMinusDefinition();
183	G4ParticleDefinition* thePionZero = G4PionZero::PionZeroDefinition();
184
185	G4double* GammaKineticEnergyCuts ;
186	G4double* ElectronKineticEnergyCuts ;
187	G4double* PositronKineticEnergyCuts ;
188	G4double* ParticleKineticEnergyCuts ;
189
190	theMaterialTable = G4Material::GetMaterialTable() ;
191
192	G4double cutinrange,CutInRangeele,CutInRangepos ;
193
194	G4ParticleDefinition* theParticle ;
195	G4GPILSelection selection;
196
197	G4double energy, momentum, mass;
198	G4ProcessVector* palongget ;
199	G4ProcessVector* palongdo ;
200	G4ProcessVector* ppostget ;
201	G4ProcessVector* ppostdo ;
202
203	G4String confirm ;
204
205	G4cout << " Do you want the mu+ as particle (yes/no)? " << std::flush;
206	G4cin >> confirm ;
207	if(confirm == "yes")
208	{
209	mass=theMuonPlus->GetPDGMass();
210	theParticle = theMuonPlus;
211	}
212	else
213	{
214	G4cout << " Do you want the mu- as particle (yes/no)? " << std::flush;
215	G4cin >> confirm ;
216	if(confirm == "yes")
217	{
218	mass=theMuonMinus->GetPDGMass();
219	theParticle = theMuonMinus;
220	}
221	}
222
223
224	energy = 1.*GeV + mass ;
225	momentum=std::sqrt(energyenergy-massmass) ;
226
227	G4ParticleMomentum theMomentum(momentum,0.,0.);
228
229	G4double pModule = theMomentum.mag();
230
231	G4DynamicParticle aParticle(theParticle,energy,theMomentum);
232
233	aParticle.SetKineticEnergy(energy-mass);
234
235
236	G4MuIonisation theParticleIonisation ;
237	G4ProcessManager* theParticleProcessManager = theParticle->GetProcessManager();
238	theParticleProcessManager->AddProcess(&theParticleIonisation,-1,0,0) ;
239	G4MuBremsstrahlung theParticleBremsstrahlung ;
240	theParticleProcessManager->AddProcess(&theParticleBremsstrahlung,-1,-1,1) ;
241	G4MuPairProduction theParticlePairProduction ;
242	theParticleProcessManager->AddProcess(&theParticlePairProduction,-1,-1,2) ;
243	G4MuNuclearInteraction theParticleNuclearInteraction ;
244	theParticleProcessManager->AddProcess(&theParticleNuclearInteraction,-1,-1,3) ;
245
246	G4ForceCondition cond ;
247	G4ForceCondition* condition = &cond ;
248
249	G4double currentSafety ;
250	G4double& refsafety=currentSafety;
251
252	theTimer.Start() ;
253	G4cout << "cut for GAMMA in mm =" ;
254	G4cin >> cutinrange ; cutinrange *= mm;
255	theGamma->SetCuts(cutinrange) ;
256	G4cout << "gamma,cut in range(mm)=" << theGamma->GetCuts()/mm << G4endl ;
257
258	GammaKineticEnergyCuts = theGamma->GetCutsInEnergy() ;
259	for (G4int icut=0; icut<theMaterialTable->length(); icut++)
260	{
261	G4cout << "material index=" << icut << " kin.energy cut(MeV)=" <<
262	GammaKineticEnergyCuts[icut]/MeV << G4endl ;
263	}
264	theTimer.Stop() ;
265	G4cout << " time = " << theTimer.GetUserElapsed() << G4endl;
266
267	theTimer.Start() ;
268	G4cout << "cut for ELECTRON in mm =" ;
269	G4cin >> cutinrange ; cutinrange *= mm;
270	CutInRangeele = cutinrange ;
271	theElectron->SetCuts(cutinrange) ;
272	G4cout << "electron,cut in range(mm)=" << theElectron->GetCuts()/mm << G4endl ;
273
274	ElectronKineticEnergyCuts = theElectron->GetCutsInEnergy() ;
275	for ( icut=0; icut<theMaterialTable->length(); icut++)
276	{
277	G4cout << "material index=" << icut << " kin.energy cut(MeV)=" <<
278	ElectronKineticEnergyCuts[icut]/MeV << G4endl ;
279	}
280	theTimer.Stop() ;
281	G4cout << " time = " << theTimer.GetUserElapsed() << G4endl;
282
283	theTimer.Start() ;
284	G4cout << "cut for POSITRON in mm =" ;
285	G4cin >> cutinrange ; cutinrange *= mm;
286	CutInRangepos = cutinrange ;
287	thePositron->SetCuts(cutinrange) ;
288	G4cout << "positron,cut in range(mm)=" << thePositron->GetCuts()/mm << G4endl ;
289
290	PositronKineticEnergyCuts = thePositron->GetCutsInEnergy() ;
291	for ( icut=0; icut<theMaterialTable->length(); icut++)
292	{
293	G4cout << "material index=" << icut << " kin.energy cut(MeV)=" <<
294	PositronKineticEnergyCuts[icut]/MeV << G4endl ;
295	}
296	theTimer.Stop() ;
297	G4cout << " time = " << theTimer.GetUserElapsed() << G4endl;
298
299	theTimer.Start() ;
300	G4cout << "cut for muons in mm =" ;
301	G4cin >> cutinrange ; cutinrange *= mm;
302	theParticle->SetCuts(cutinrange) ;
303
304	G4cout << "cut in range(mm)=" << theParticle->GetLengthCuts()/mm << G4endl ;
305
306	ParticleKineticEnergyCuts = theParticle->GetEnergyCuts() ;
307
308	for ( icut=0; icut<theMaterialTable->length(); icut++)
309	{
310	G4cout << "material index=" << icut << " kin.energy cut(MeV)=" <<
311	ParticleKineticEnergyCuts[icut]/MeV << G4endl ;
312	}
313	theTimer.Stop() ;
314	G4cout << " time = " << theTimer.GetUserElapsed() << G4endl;
315
316
317	G4cout << " ------ ----- " << G4endl ;
318
319	palongget = aParticle.GetDefinition()->GetProcessManager()
320	->GetAlongStepProcessVector(typeGPIL);
321	ppostget = aParticle.GetDefinition()->GetProcessManager()
322	->GetPostStepProcessVector(typeGPIL);
323	palongdo = aParticle.GetDefinition()->GetProcessManager()
324	->GetAlongStepProcessVector(typeDoIt);
325	ppostdo = aParticle.GetDefinition()->GetProcessManager()
326	->GetPostStepProcessVector(typeDoIt);
327
328	//---------------------------------- Physics --------------------------------
329
330	G4int itry=1, Ntry=1, Nstart, ir;
331	G4double r ;
332
333	//**************************************************************************
334	const G4int Nbin=137 ;
335	G4double TkinMeV[Nbin] =
336	{0.00001,0.000015,0.00002,0.00003,0.00004,0.00005,0.00006,0.00008,
337	0.0001,0.00015,0.0002,0.0003,0.0004,0.0005,0.0006,0.0008,
338	0.001,0.0015,0.002,0.003,0.004,0.005,0.006,0.008,
339	0.01,0.015,0.02,0.03,0.04,0.05,0.06,0.08,
340	0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.8,
341	1.,1.5,2.,3.,4.,5.,6.,8.,
342	10.,15.,20.,30.,40.,50.,60.,80.,
343	100.,150.,200.,300.,400.,500.,600.,800.,
344	1.0e3,1.5e3,2.0e3,3.0e3,4.0e3,5.0e3,6.0e3,8.0e3,
345	1.0e4,1.5e4,2.0e4,3.0e4,4.0e4,5.0e4,6.0e4,8.0e4,
346	1.0e5,1.5e5,2.0e5,3.0e5,4.0e5,5.0e5,6.0e5,8.0e5,
347	1.0e6,1.5e6,2.0e6,3.0e6,4.0e6,5.0e6,6.0e6,8.0e6,
348	1.0e7,1.5e7,2.0e7,3.0e7,4.0e7,5.0e7,6.0e7,8.0e7,
349	1.0e8,1.5e8,2.0e8,3.0e8,4.0e8,5.0e8,6.0e8,8.0e8,
350	1.0e9,1.5e9,2.0e9,3.0e9,4.0e9,5.0e9,6.0e9,8.0e9,
351	1.e10,1.5e10,2.e10,3.e10,4.e10,5.e10,6.e10,8.e10,
352	1.e11,1.5e11,2.e11,3.e11,4.e11,5.e11,6.e11,8.e11,
353	1.e12} ;
354	for (G4int k=0; k<Nbin; k++) TkinMeV[k] *= MeV;
355
356	G4int J=-1 ;
357
358	G4double lambda,trueStep,geomStep,stepLimit,
359	previousStepSize,currentMinimumStep ;
360	G4ParticleChange* aParticleChange ;
361
362	G4double T,dEdx,range ;
363
364	NEXTMATERIAL: ;
365	J = J+1 ;
366	if ( J >= theMaterialTable->length() )
367	{ G4cout << "that was the last material in the table --> STOP" << G4endl;
368	return EXIT_FAILURE ; }
369
370	apttoMaterial = (*theMaterialTable)[ J ] ;
371	MaterialName = apttoMaterial->GetName() ;
372	G4cout << "material=" << MaterialName << G4endl ;
373	G4cout << "Do you want the Energyloss test 1. for this material?" << G4endl ;
374	G4cout << "type a positive number if the answer is YES" << G4endl ;
375	G4cout << "type a negative number if the answer is NO " << G4endl ;
376	G4int icont ;
377	G4cin >> icont ;
378	if ( icont < 0 )
379	goto NEXTMATERIAL ;
380
381	//---------- Volume definition ---------------------
382
383	G4VPhysicalVolume* myVolume ;
384
385	myVolume = BuildVolume(apttoMaterial) ;
386
387	//--------- track and Step definition (for this test ONLY!)------------
388	G4ThreeVector aPosition(0.,0.,0.);
389	const G4ThreeVector aDirection(0.,0.,1.) ;
390	const G4ThreeVector transl(0.,0.,0.) ;
391
392	G4double aTime = 0. ;
393
394	G4Track* tracke = new G4Track(&aParticle,aTime,aPosition) ;
395	G4Track& trackele = (*tracke) ;
396	//(*tracke).SetVolume(myVolume) ;
397	G4GRSVolume* touche = new G4GRSVolume(myVolume,NULL,transl);
398	(*tracke).SetTouchable(touche);
399
400	(*tracke).SetMomentumDirection(aDirection) ;
401
402
403	G4Step* Step = new G4Step() ;
404	G4Step& Step = (*Step) ;
405	tracke->SetStep(Step);
406
407	G4StepPoint* aPoint = new G4StepPoint();
408	(*aPoint).SetPosition(aPosition) ;
409	G4double safety = 10000.*cm ;
410	(*aPoint).SetSafety(safety) ;
411
412	(*Step).SetPostStepPoint(aPoint) ;
413
414	//**************************************************************************
415
416	G4cout << G4endl;
417	G4cout <<" " << MaterialName << " Energyloss test 1." << G4endl ;
418	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
419	G4cout << G4endl ;
420	G4cout << "kin.en.(MeV) dE/dx(MeV/mm) range(mm) Step(mm)"
421	<< " dEdx(GeVcm2/g)" << G4endl ;
422	G4cout << G4endl ;
423
424
425	for ( G4int i=0 ; i<Nbin ; i++)
426	{
427	trueStep = cutinrange ;
428	previousStepSize = cutinrange ;
429	currentMinimumStep = trueStep ;
430	(*tracke).SetKineticEnergy(TkinMeV[i]) ;
431	stepLimit = (*palongget)(0)->AlongStepGetPhysicalInteractionLength(
432	trackele,
433	previousStepSize,
434	currentMinimumStep,
435	refsafety,
436	&selection) ;
437
438	dEdx = theParticleIonisation.GetdEdx() ;
439	range = theParticleIonisation.GetRangeNow() ;
440	T = TkinMeV[i] ;
441
442	G4cout <<" " << T/MeV << " " << dEdx/(MeV/mm) << " " ;
443	G4cout << range/mm << " " << stepLimit/mm << " "
444	<< dEdx/(GeV/cm)/(apttoMaterial->GetDensity()/(g/cm3)) << G4endl ;
445
446	}
447
448	G4cout << G4endl;
449
450	ENERGYLOSS2: ;
451
452	G4cout << "material=" << MaterialName << G4endl ;
453	G4cout << "Do you want the Energyloss test 2. for this material?" << G4endl ;
454	G4cout << "type a positive number if the answer is YES" << G4endl ;
455	G4cout << "type a negative number if the answer is NO " << G4endl ;
456	G4cin >> icont ;
457	if ( icont < 0 )
458	goto ENERGYLOSS3 ;
459
460	G4double TMeV,stepmm,stepmx,meanloss,lossnow ;
461
462
463	G4cout << "give an energy value in MeV " ;
464	G4cin >> TMeV ; TMeV *= MeV;
465
466	trueStep = cutinrange ;
467	previousStepSize = cutinrange ;
468	currentMinimumStep = trueStep ;
469	(*tracke).SetKineticEnergy(TMeV) ;
470	stepmx = (*palongget)(0)->AlongStepGetPhysicalInteractionLength(
471	trackele,
472	previousStepSize,
473	currentMinimumStep,
474	refsafety,
475	&selection);
476
477	G4cout << " give a steplength in mm , the max. meaningful Step is " << stepmx/mm << " mm" <<G4endl;
478	G4cout << "Step:" ;
479	G4cin >> stepmm ; stepmm *= mm;
480
481	(*Step).SetTrack(tracke) ;
482	(*Step).SetStepLength(stepmm);
483
484
485	aParticleChange = (G4ParticleChange*)
486	((*palongdo)(0)->AlongStepDoIt(trackele,Step));
487	meanloss = theParticleIonisation.GetMeanLoss() ;
488	lossnow = TMeV-(*aParticleChange).GetEnergyChange();
489
490	G4cout << G4endl;
491	G4cout <<" " << MaterialName << " Energyloss test 2." << G4endl ;
492	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
493	G4cout << G4endl ;
494	G4cout << "kin.en.(MeV) Step(mm) meanloss(MeV) act.loss(MeV)" << G4endl ;
495	G4cout << TMeV/MeV << " " << stepmm/mm << " " << meanloss/MeV << " " << lossnow/MeV << G4endl ;
496	G4cout << " status change:" << (*aParticleChange).GetStatusChange() << G4endl ;
497	G4cout << G4endl ;
498
499	goto ENERGYLOSS2 ;
500
501	ENERGYLOSS3: ;
502
503	G4cout << "material=" << MaterialName << G4endl ;
504	G4cout << "Do you want the Energyloss test 3. for this material?" << G4endl ;
505	G4cout << "type a positive number if the answer is YES" << G4endl ;
506	G4cout << "type a negative number if the answer is NO " << G4endl ;
507	G4cin >> icont ;
508	if ( icont < 0 )
509	goto DELTARAY1 ;
510
511	G4cout << "give an energy value in MeV " ;
512	G4cin >> TMeV ; TMeV *= MeV;
513
514	trueStep = cutinrange ;
515	previousStepSize = cutinrange ;
516	currentMinimumStep = trueStep ;
517	(*tracke).SetKineticEnergy(TMeV) ;
518	stepmx = (*palongget)(0)->AlongStepGetPhysicalInteractionLength(
519	trackele,
520	previousStepSize,
521	currentMinimumStep,
522	refsafety,
523	&selection);
524
525	G4cout << " give a steplength in mm , the max. meaningful Step is " << stepmx << " mm" <<G4endl;
526	G4cout << "Step:" ;
527	G4cin >> stepmm ; stepmm *= mm;
528
529	(*Step).SetTrack(tracke) ;
530	(*Step).SetStepLength(stepmm);
531
532
533	G4cout << " give number of events you want " ;
534	G4int nbev,ibev ;
535	G4cin >> nbev ;
536
537	meanloss=0.;
538	theTimer.Start();
539
540	for ( ibev=0; ibev<nbev; ibev++)
541	{
542	aParticleChange = (G4ParticleChange*)
543	((*palongdo)(0)->AlongStepDoIt(trackele,Step));
544	lossnow = TMeV-(*aParticleChange).GetEnergyChange();
545
546	meanloss += lossnow ;
547	}
548
549	theTimer.Stop();
550	meanloss /= nbev ;
551	G4cout << G4endl;
552	G4cout <<" " << MaterialName << " Energyloss test 3." << G4endl ;
553	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
554	G4cout << G4endl ;
555	G4cout << "kin.en.(MeV) Step(mm) meanloss(MeV) time/event(sec) " << G4endl ;
556	G4cout << TMeV/MeV << " " << stepmm/mm << " " << meanloss/MeV << " " <<
557	theTimer.GetUserElapsed()/nbev << G4endl ;
558	G4cout << G4endl ;
559
560	goto ENERGYLOSS3 ;
561
562	DELTARAY1: ;
563	G4cout << "material=" << MaterialName << G4endl ;
564	G4cout << "Do you want the delta ray test 1. for this material?" << G4endl ;
565	G4cout << "type a positive number if the answer is YES" << G4endl ;
566	G4cout << "type a negative number if the answer is NO " << G4endl ;
567	G4cin >> icont ;
568	if ( icont < 0 )
569	goto DELTARAY2 ;
570
571
572	G4cout << G4endl;
573	G4cout <<" " << MaterialName << " delta ray test 1." << G4endl ;
574	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
575	G4cout << G4endl ;
576	G4cout << "kin.en.(MeV) mean free path(mm)" << G4endl ;
577	G4cout << G4endl ;
578
579	for ( i=0 ; i<Nbin ; i++)
580	{
581
582	previousStepSize = cutinrange ;
583	(*tracke).SetKineticEnergy(TkinMeV[i]) ;
584	stepLimit = theParticleIonisation.GetMeanFreePath(
585	trackele,
586	previousStepSize,
587	condition) ;
588
589	T = TkinMeV[i] ;
590
591	G4cout <<" " << T/MeV << " " << stepLimit/mm << G4endl ;
592
593	}
594
595	G4cout << G4endl;
596
597	//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
598	DELTARAY2: ;
599
600	G4cout << "material=" << MaterialName << G4endl ;
601	G4cout << "Do you want the deltaray test 2. for this material?" << G4endl ;
602	G4cout << "type a positive number if the answer is YES" << G4endl ;
603	G4cout << "type a negative number if the answer is NO " << G4endl ;
604	G4cin >> icont ;
605
606	G4double newenergy,dx,dy,dz,Tdelta,ddx,ddy,ddz ;
607	G4int nd ;
608	const G4ThreeVector* momdir ;
609	G4ParticleMomentum ddir ;
610
611	if ( icont < 0 )
612	goto BREMS1 ;
613
614	G4cout << "give an energy value in MeV " ;
615	G4cin >> TMeV ; TMeV *= MeV;
616
617	stepmm = 1.*mm ;
618
619	(*Step).SetTrack(tracke) ;
620	(*Step).SetStepLength(stepmm);
621
622
623	(*tracke).SetKineticEnergy(TMeV) ;
624	aParticleChange = (G4ParticleChange*)
625	((*ppostdo)(0)->PostStepDoIt(trackele,Step));
626
627	newenergy=(*aParticleChange).GetEnergyChange() ;
628	momdir=(*aParticleChange).GetMomentumChange();
629	dx = (*momdir).x();
630	dy = (*momdir).y();
631	dz = (*momdir).z();
632	nd=aParticleChange->GetNumberOfSecondaries();
633
634	if(nd>0)
635	{
636	Tdelta=aParticleChange->GetSecondary(0)->GetKineticEnergy();
637	ddir=aParticleChange->GetSecondary(0)->
638	GetMomentumDirection();
639	ddx = (ddir).x();
640	ddy = (ddir).y();
641	ddz = (ddir).z();
642	}
643	(*aParticleChange).Clear();
644
645
646	G4cout << G4endl;
647	G4cout <<" " << MaterialName << " delta ray test 2." << G4endl ;
648	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
649	G4cout << G4endl ;
650	G4cout << "T=" << TMeV/MeV << " newT=" << newenergy/MeV << " (MeV)" << G4endl ;
651	G4cout << " status change:" << (*aParticleChange).GetStatusChange() << G4endl ;
652	if(nd>0)
653	G4cout << "Tdelta=" << Tdelta/MeV << G4endl ;
654	G4cout << "new direction:" << dx << " " << dy << " " << dz << G4endl;
655	if(nd>0)
656	G4cout << "delta direction:" << ddx << " " << ddy << " " << ddz << G4endl ;
657	//............................................
658	nbev=50000 ;
659	mloss = 0. ;
660	sloss = 0. ;
661	theTimer.Start();
662
663	for (ibev=0 ; ibev<nbev ; ibev++ )
664	{
665	aParticleChange = (G4ParticleChange*)
666	((*ppostdo)(0)->PostStepDoIt(trackele,Step));
667	nd=aParticleChange->GetNumberOfSecondaries();
668	Tdelta=0. ;
669	if(nd>0)
670	{
671	Tdelta=aParticleChange->GetSecondary(0)->GetKineticEnergy();
672	}
673	mloss += Tdelta ;
674	sloss += Tdelta*Tdelta ;
675	(*aParticleChange).Clear();
676	}
677	theTimer.Stop();
678	mloss /= nbev ;
679	sloss /= nbev ;
680	sloss = (sloss-mloss*mloss)/nbev ;
681	if(sloss>0.)
682	sloss = std::sqrt(sloss) ;
683	else
684	sloss = 0. ;
685
686	previousStepSize = cutinrange ;
687	stepLimit = theParticleIonisation.GetMeanFreePath(
688	trackele,
689	previousStepSize,
690	condition) ;
691	dEdxdelta=mloss/stepLimit ;
692
693	G4cout << " mean energy loss due to delta production (in MeV)=" <<
694	mloss/MeV << " +- " << sloss/MeV << G4endl ;
695	G4cout << " dE/dx due to delta production (in MeV/mm,from " <<
696	nbev << " events )=" << dEdxdelta/(MeV/mm) << G4endl ;
697	G4cout << "in GeVcm2/g ="
698	<< dEdxdelta/(GeV/cm)/(apttoMaterial->GetDensity()/(g/cm3)) << G4endl ;
699	G4cout << " time/delta =" << theTimer.GetUserElapsed()/nbev << G4endl ;
700	G4cout << G4endl ;
701
702	//..................................
703
704	goto DELTARAY2 ;
705
706	BREMS1: ;
707	G4cout << "material=" << MaterialName << G4endl ;
708	G4cout << "Do you want the brems test 1. for this material?" << G4endl ;
709	G4cout << "type a positive number if the answer is YES" << G4endl ;
710	G4cout << "type a negative number if the answer is NO " << G4endl ;
711	G4cin >> icont ;
712	if ( icont < 0 )
713	goto BREMS2 ;
714
715	G4cout << G4endl;
716	G4cout <<" " << MaterialName << " brems test 1." << G4endl ;
717	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
718	G4cout << G4endl ;
719	G4cout << "kin.en.(MeV) mean free path(mm)" << G4endl ;
720	G4cout << G4endl ;
721
722	for ( i=0 ; i<Nbin ; i++)
723	{
724
725	previousStepSize = CutInRangeele ;
726	(*tracke).SetKineticEnergy(TkinMeV[i]) ;
727	stepLimit = theParticleBremsstrahlung.GetMeanFreePath(
728	trackele,
729	previousStepSize,
730	condition) ;
731
732	T = TkinMeV[i] ;
733
734	G4cout <<" " << T/MeV << " " << stepLimit/mm << G4endl ;
735
736	}
737
738	G4cout << G4endl;
739
740	//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
741	BREMS2: ;
742
743	G4cout << "material=" << MaterialName << G4endl ;
744	G4cout << "Do you want the brems test 2. for this material?" << G4endl ;
745	G4cout << "type a positive number if the answer is YES" << G4endl ;
746	G4cout << "type a negative number if the answer is NO " << G4endl ;
747	G4cin >> icont ;
748	if ( icont < 0 )
749	goto PAIR1 ;
750
751	G4cout << "give an energy value in MeV " ;
752	G4cin >> TMeV ; TMeV *= MeV;
753
754	stepmm = 1. ;
755
756	(*Step).SetTrack(tracke) ;
757	(*Step).SetStepLength(stepmm);
758	(*tracke).SetKineticEnergy(TMeV) ;
759	aParticleChange = (G4ParticleChange*)
760	((*ppostdo)(1)->PostStepDoIt(trackele,Step));
761
762	newenergy=(*aParticleChange).GetEnergyChange() ;
763	momdir=(*aParticleChange).GetMomentumChange();
764	dx = (*momdir).x();
765	dy = (*momdir).y();
766	dz = (*momdir).z();
767	nd=aParticleChange->GetNumberOfSecondaries();
768
769	if(nd>0)
770	{
771	Tdelta=aParticleChange->GetSecondary(0)->GetKineticEnergy();
772	ddir=aParticleChange->GetSecondary(0)->
773	GetMomentumDirection();
774	ddx = (ddir).x();
775	ddy = (ddir).y();
776	ddz = (ddir).z();
777	}
778	(*aParticleChange).Clear();
779
780
781	G4cout << G4endl;
782	G4cout <<" " << MaterialName << " brems test 2." << G4endl ;
783	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
784	G4cout << G4endl ;
785	G4cout << "T=" << TMeV/MeV << " newT=" << newenergy/MeV << " (MeV)" << G4endl ;
786	G4cout << " status change:" << (*aParticleChange).GetStatusChange() << G4endl ;
787	if(nd>0)
788	G4cout << "Tgamma=" << Tdelta/MeV << G4endl ;
789	G4cout << "new direction:" << dx << " " << dy << " " << dz << G4endl;
790	if(nd>0)
791	G4cout << "gamma direction:" << ddx << " " << ddy << " " << ddz << G4endl ;
792
793	nbev=50000 ;
794	mloss = 0. ;
795	sloss = 0. ;
796
797	theTimer.Start();
798
799	for (ibev=0 ; ibev<nbev ; ibev++ )
800	{
801	aParticleChange = (G4ParticleChange*)
802	((*ppostdo)(1)->PostStepDoIt(trackele,Step));
803	nd=aParticleChange->GetNumberOfSecondaries();
804	Tdelta=0. ;
805	if(nd>0)
806	{
807	Tdelta=aParticleChange->GetSecondary(0)->GetKineticEnergy();
808	}
809	mloss += Tdelta ;
810	sloss += Tdelta*Tdelta ;
811	(*aParticleChange).Clear();
812	}
813	theTimer.Stop();
814	mloss /= nbev ;
815	sloss /= nbev ;
816	sloss = std::sqrt((sloss-mloss*mloss)/nbev) ;
817
818	previousStepSize = cutinrange ;
819	stepLimit = theParticleBremsstrahlung.GetMeanFreePath(
820	trackele,
821	previousStepSize,
822	condition) ;
823	dEdxbrems=mloss/stepLimit ;
824
825	G4cout << " mean energy loss due to bremsstrahlung (in MeV)=" <<
826	mloss/MeV << " +- " << sloss/MeV << G4endl ;
827	G4cout << " dE/dx due to bremsstrahlung (in MeV/mm,from " <<
828	nbev << " events )=" << dEdxbrems/(MeV/mm) << G4endl;
829	G4cout << "in GeVcm2/g ="
830	<< dEdxbrems/(GeV/cm)/(apttoMaterial->GetDensity()/(g/cm3)) << G4endl ;
831	G4cout << " time/brems =" << theTimer.GetUserElapsed()/nbev << G4endl ;
832	G4cout << G4endl ;
833
834	goto BREMS2 ;
835
836	PAIR1: ;
837	G4cout << "material=" << MaterialName << G4endl ;
838	G4cout << "Do you want the pair test 1. for this material?" << G4endl ;
839	G4cout << "type a positive number if the answer is YES" << G4endl ;
840	G4cout << "type a negative number if the answer is NO " << G4endl ;
841	G4cin >> icont ;
842	if ( icont < 0 )
843	goto PAIR2 ;
844
845	G4cout << G4endl;
846	G4cout <<" " << MaterialName << " pair test 1." << G4endl ;
847	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
848	G4cout << G4endl ;
849	G4cout << "kin.en.(MeV) mean free path(mm)" << G4endl ;
850	G4cout << G4endl ;
851
852	for ( i=0 ; i<Nbin ; i++)
853	{
854
855	previousStepSize = CutInRangeele ;
856	(*tracke).SetKineticEnergy(TkinMeV[i]) ;
857	stepLimit = theParticlePairProduction.GetMeanFreePath(
858	trackele,
859	previousStepSize,
860	condition) ;
861
862	T = TkinMeV[i] ;
863
864	G4cout <<" " << T/MeV << " " << stepLimit/mm << G4endl ;
865
866	}
867
868	G4cout << G4endl;
869
870	//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
871	PAIR2: ;
872
873	G4double Tdelta1,ddx1,ddy1,ddz1,Tdelta2,ddx2,ddy2,ddz2 ;
874	G4ParticleMomentum ddir1,ddir2 ;
875	G4double dEdxpair ;
876
877	G4cout << "material=" << MaterialName << G4endl ;
878	G4cout << "Do you want the pair test 2. for this material?" << G4endl ;
879	G4cout << "type a positive number if the answer is YES" << G4endl ;
880	G4cout << "type a negative number if the answer is NO " << G4endl ;
881	G4cin >> icont ;
882	if ( icont < 0 )
883	goto NUCL1 ;
884
885	G4cout << "give an energy value in MeV " ;
886	G4cin >> TMeV ; TMeV *= MeV;
887
888	stepmm = 1. ;
889
890	(*Step).SetTrack(tracke) ;
891	(*Step).SetStepLength(stepmm);
892	(*tracke).SetKineticEnergy(TMeV) ;
893	aParticleChange = (G4ParticleChange*)
894	((*ppostdo)(2)->PostStepDoIt(trackele,Step));
895
896	newenergy=(*aParticleChange).GetEnergyChange() ;
897	momdir=(*aParticleChange).GetMomentumChange();
898	dx = (*momdir).x();
899	dy = (*momdir).y();
900	dz = (*momdir).z();
901	nd=aParticleChange->GetNumberOfSecondaries();
902
903	if(nd>0)
904	{
905	Tdelta1=aParticleChange->GetSecondary(0)->GetKineticEnergy();
906	ddir1=aParticleChange->GetSecondary(0)->
907	GetMomentumDirection();
908	ddx1 = (ddir1).x();
909	ddy1 = (ddir1).y();
910	ddz1 = (ddir1).z();
911	}
912	if(nd>1)
913	{
914	Tdelta2=aParticleChange->GetSecondary(1)->GetKineticEnergy();
915	ddir2=aParticleChange->GetSecondary(1)->
916	GetMomentumDirection();
917	ddx2 = (ddir2).x();
918	ddy2 = (ddir2).y();
919	ddz2 = (ddir2).z();
920	}
921
922	(*aParticleChange).Clear();
923
924
925	G4cout << G4endl;
926	G4cout <<" " << MaterialName << " pair test 2." << G4endl ;
927	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
928	G4cout << G4endl ;
929	G4cout << "T=" << TMeV/MeV << " newT=" << newenergy/MeV << " (MeV)" << G4endl ;
930	G4cout << " status change:" << (*aParticleChange).GetStatusChange() << G4endl ;
931	if(nd>0)
932	G4cout << "T1=" << Tdelta1/MeV << G4endl ;
933	if(nd>1)
934	G4cout << "T2=" << Tdelta2/MeV << G4endl ;
935
936	G4cout << "new direction:" << dx << " " << dy << " " << dz << G4endl;
937	if(nd>0)
938	G4cout << "direction1:" << ddx1 << " " << ddy1 << " " << ddz1 << G4endl ;
939	if(nd>1)
940	G4cout << "direction2:" << ddx2 << " " << ddy2 << " " << ddz2 << G4endl ;
941
942	nbev=50000 ;
943	mloss = 0. ;
944	sloss = 0. ;
945
946	theTimer.Start();
947
948	for (ibev=0 ; ibev<nbev ; ibev++ )
949	{
950	aParticleChange = (G4ParticleChange*)
951	((*ppostdo)(2)->PostStepDoIt(trackele,Step));
952	nd=aParticleChange->GetNumberOfSecondaries();
953	Tdelta=0. ;
954	if(nd>0)
955	{
956	Tdelta1=aParticleChange->GetSecondary(0)->GetKineticEnergy();
957	Tdelta=Tdelta1;
958	}
959	if(nd>1)
960	{
961	Tdelta2=aParticleChange->GetSecondary(1)->GetKineticEnergy();
962	Tdelta+=Tdelta2;
963	}
964
965	mloss += Tdelta ;
966	sloss += Tdelta*Tdelta ;
967	(*aParticleChange).Clear();
968	}
969	theTimer.Stop();
970	mloss /= nbev ;
971	sloss /= nbev ;
972	sloss = std::sqrt((sloss-mloss*mloss)/nbev) ;
973
974	previousStepSize = cutinrange ;
975	stepLimit = theParticlePairProduction.GetMeanFreePath(
976	trackele,
977	previousStepSize,
978	condition) ;
979	dEdxpair=mloss/stepLimit ;
980
981	G4cout << " mean energy loss due to pair production (in MeV)=" <<
982	mloss/MeV << " +- " << sloss/MeV << G4endl ;
983	G4cout << " dE/dx due to pair production (in MeV/mm,from " <<
984	nbev << " events )=" << dEdxpair/(MeV/mm) << G4endl;
985	G4cout << "in GeVcm2/g ="
986	<< dEdxpair/(GeV/cm)/(apttoMaterial->GetDensity()/(g/cm3)) << G4endl ;
987	G4cout << " time/pair =" << theTimer.GetUserElapsed()/nbev << G4endl ;
988	G4cout << G4endl ;
989
990	goto PAIR2 ;
991
992	NUCL1: ;
993	G4cout << "material=" << MaterialName << G4endl ;
994	G4cout << "Do you want the nucl.int. test 1. for this material?" << G4endl ;
995	G4cout << "type a positive number if the answer is YES" << G4endl ;
996	G4cout << "type a negative number if the answer is NO " << G4endl ;
997	G4cin >> icont ;
998	if ( icont < 0 )
999	goto NUCL2 ;
1000
1001	G4cout << G4endl;
1002	G4cout <<" " << MaterialName << " nucl.int. test 1." << G4endl ;
1003	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
1004	G4cout << G4endl ;
1005	G4cout << "kin.en.(MeV) mean free path(mm)" << G4endl ;
1006	G4cout << G4endl ;
1007
1008	for ( i=0 ; i<Nbin ; i++)
1009	{
1010
1011	previousStepSize = CutInRangeele ;
1012	(*tracke).SetKineticEnergy(TkinMeV[i]) ;
1013	stepLimit = theParticleNuclearInteraction.GetMeanFreePath(
1014	trackele,
1015	previousStepSize,
1016	condition) ;
1017
1018	T = TkinMeV[i] ;
1019
1020	G4cout <<" " << T/MeV << " " << stepLimit/mm << G4endl ;
1021
1022	}
1023
1024	G4cout << G4endl;
1025
1026	//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1027	NUCL2: ;
1028
1029	G4cout << "material=" << MaterialName << G4endl ;
1030	G4cout << "Do you want the nucl.int. test 2. for this material?" << G4endl ;
1031	G4cout << "type a positive number if the answer is YES" << G4endl ;
1032	G4cout << "type a negative number if the answer is NO " << G4endl ;
1033	G4cin >> icont ;
1034	if ( icont < 0 )
1035	goto NEXTMATERIAL ;
1036
1037	G4cout << "give an energy value in MeV " ;
1038	G4cin >> TMeV ; TMeV *= MeV;
1039
1040	stepmm = 1. ;
1041
1042	(*Step).SetTrack(tracke) ;
1043	(*Step).SetStepLength(stepmm);
1044	(*tracke).SetKineticEnergy(TMeV) ;
1045	aParticleChange = (G4ParticleChange*)
1046	((*ppostdo)(3)->PostStepDoIt(trackele,Step));
1047
1048	newenergy=(*aParticleChange).GetEnergyChange() ;
1049	momdir=(*aParticleChange).GetMomentumChange();
1050	dx = (*momdir).x();
1051	dy = (*momdir).y();
1052	dz = (*momdir).z();
1053	nd=aParticleChange->GetNumberOfSecondaries();
1054
1055	if(nd>0)
1056	{
1057	Tdelta1=aParticleChange->GetSecondary(0)->GetKineticEnergy();
1058	ddir1=aParticleChange->GetSecondary(0)->
1059	GetMomentumDirection();
1060	ddx1 = (ddir1).x();
1061	ddy1 = (ddir1).y();
1062	ddz1 = (ddir1).z();
1063	}
1064
1065	(*aParticleChange).Clear();
1066
1067
1068	G4cout << G4endl;
1069	G4cout <<" " << MaterialName << " nucl.int. test 2." << G4endl ;
1070	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
1071	G4cout << G4endl ;
1072	G4cout << "T=" << TMeV/MeV << " newT=" << newenergy/MeV << " (MeV)" << G4endl ;
1073	G4cout << " status change:" << (*aParticleChange).GetStatusChange() << G4endl ;
1074	if(nd>0)
1075	G4cout << "T1=" << Tdelta1/MeV << G4endl ;
1076	if(nd>1)
1077	G4cout << "T2=" << Tdelta2/MeV << G4endl ;
1078
1079	G4cout << "new direction:" << dx << " " << dy << " " << dz << G4endl;
1080	if(nd>0)
1081	G4cout << "direction1:" << ddx1 << " " << ddy1 << " " << ddz1 << G4endl ;
1082	if(nd>1)
1083	G4cout << "direction2:" << ddx2 << " " << ddy2 << " " << ddz2 << G4endl ;
1084
1085	nbev=50000 ;
1086	mloss = 0. ;
1087	sloss = 0. ;
1088
1089	theTimer.Start();
1090
1091	for (ibev=0 ; ibev<nbev ; ibev++ )
1092	{
1093	aParticleChange = (G4ParticleChange*)
1094	((*ppostdo)(3)->PostStepDoIt(trackele,Step));
1095	nd=aParticleChange->GetNumberOfSecondaries();
1096	Tdelta=0. ;
1097	// secondary is not generated ....................................
1098	Tdelta1 = TMeV - (*aParticleChange).GetEnergyChange() ;
1099	Tdelta=Tdelta1;
1100	if(nd>0)
1101	{
1102	// secondary is not generated ....................................
1103	//Tdelta1=aParticleChange->GetSecondary(0)->GetKineticEnergy();
1104	Tdelta=Tdelta1;
1105	}
1106
1107	mloss += Tdelta ;
1108	sloss += Tdelta*Tdelta ;
1109	(*aParticleChange).Clear();
1110	}
1111	theTimer.Stop();
1112	mloss /= nbev ;
1113	sloss /= nbev ;
1114	sloss = std::sqrt((sloss-mloss*mloss)/nbev) ;
1115
1116	previousStepSize = cutinrange ;
1117	stepLimit = theParticleNuclearInteraction.GetMeanFreePath(
1118	trackele,
1119	previousStepSize,
1120	condition) ;
1121	dEdxpair=mloss/stepLimit ;
1122
1123	G4cout << " mean energy loss due to nucl.int. production (in MeV)=" <<
1124	mloss/MeV << " +- " << sloss/MeV << G4endl ;
1125	G4cout << " dE/dx due to nucl.int. production (in MeV/mm,from " <<
1126	nbev << " events )=" << dEdxpair/(MeV/mm) << G4endl;
1127	G4cout << "in GeVcm2/g ="
1128	<< dEdxpair/(GeV/cm)/(apttoMaterial->GetDensity()/(g/cm3)) << G4endl ;
1129	G4cout << " time/nucl =" << theTimer.GetUserElapsed()/nbev << G4endl ;
1130	G4cout << G4endl ;
1131
1132	goto NUCL2 ;
1133	if( J < theMaterialTable->length()-1 )
1134	goto NEXTMATERIAL ;
1135
1136	return EXIT_SUCCESS;
1137	}

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

Context Navigation

source: trunk/source/processes/electromagnetic/muons/test/muEnergyLossTest.cc @ 1347

Download in other formats: