Context Navigation

PAIenergyLossTest.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago
geant4 tag 9.4
File size: 29.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: PAIenergyLossTest.cc,v 1.8 2006/06/29 19:54:17 gunter Exp $
28	// GEANT4 tag $Name: geant4-09-04-ref-00 $
29	//
30	//
31	//-----------------------------------------------------------------
32	#include "G4ios.hh"
33	#include <fstream>
34	#include <iomanip>
35	#include "g4templates.hh"
36	#include "globals.hh"
37	#include "G4Timer.hh"
38
39	#include "G4PhysicsLogVector.hh"
40	#include "G4PhysicsFreeVector.hh"
41	#include "G4PhysicsVector.hh"
42
43	#include "G4eEnergyLoss.hh"
44	#include "G4eIonisation.hh"
45	#include "G4hEnergyLoss.hh"
46	// #include "G4PAIenergyLoss.hh"
47	#include "G4hIonisation.hh"
48	// #include "G4PAIonisation.hh"
49
50	#include "G4DynamicParticle.hh"
51	#include "G4Element.hh"
52	#include "G4Material.hh"
53	#include "G4PVPlacement.hh"
54	#include "G4LogicalVolume.hh"
55	#include "G4GRSVolume.hh"
56	#include "G4Box.hh"
57	#include "G4ProcessManager.hh"
58	#include "G4Step.hh"
59	#include "G4StepPoint.hh"
60	#include "G4Track.hh"
61
62	#include "G4Electron.hh"
63	#include "G4Positron.hh"
64	#include "G4Proton.hh"
65	#include "G4AntiProton.hh"
66	#include "G4PionPlus.hh"
67	#include "G4PionMinus.hh"
68	#include "G4KaonPlus.hh"
69	#include "G4KaonMinus.hh"
70
71	// It tests the G4PAIenergyLoss,G4PAIonisation processes -----------
72	// created by L.Urban on 06/06/97 and modified by V.Grichine on 30.11.97
73	//
74	// Modifications:
75	// 23-09-97, geometry adapted for the touchable
76	// 14.11.97 PAIonisation class is responsible for <dE/dx> calculation
77	//
78
79	G4VPhysicalVolume* BuildVolume(G4Material* matworld)
80	// it builds a simple box filled with material matword .......
81	{
82	G4Box myWorldBox= new G4Box ("WBox",10000.cm,10000.cm,10000.cm);
83
84	G4LogicalVolume *myWorldLog = new G4LogicalVolume(myWorldBox,matworld,
85	"WLog",0,0,0) ;
86
87	G4PVPlacement *myWorldPhys = new G4PVPlacement(0,G4ThreeVector(),
88	"WPhys",
89	myWorldLog,
90	0,false,0) ;
91
92
93	return myWorldPhys ;
94
95	}
96
97
98	int main()
99	{
100	//-------- set output format-------
101	G4cout.setf( std::ios::scientific, std::ios::floatfield );
102	//---write results to the file hloss -----
103	std::ofstream outFile("PAIdEdx.cc", std::ios::out ) ;
104	outFile.setf( std::ios::scientific, std::ios::floatfield );
105
106	//--------- Material definition ---------
107
108	G4double a, z, ez, density ,temperature,pressure;
109	G4State state ;
110	G4String name, symbol;
111	G4int nel;
112
113	a = 14.01*g/mole;
114	G4Element* elN = new G4Element(name="Nitrogen", symbol="N", ez=7., a);
115
116	a = 16.00*g/mole;
117	G4Element* elO = new G4Element(name="Oxigen", symbol="O", ez=8., a);
118
119	density = 1.29e-03*g/cm3;
120	state = kStateGas ;
121	temperature = 273.*kelvin ;
122	pressure = 1.*atmosphere ;
123	G4Material* Air = new G4Material(name="Air", density, nel=2 ,
124	state ,temperature , pressure ) ;
125	Air->AddElement(elN, .7);
126	Air->AddElement(elO, .3);
127
128	density = 5.85e-3*g/cm3 ;
129	a = 131.29*g/mole ;
130	G4Material* Xe = new G4Material(name="Xenon", ez=54., a, density);
131
132	density = 1.782e-3*g/cm3 ;
133	a = 39.948*g/mole ;
134	G4Material* Ar = new G4Material(name="Argon", ez=18., a, density);
135
136	a = 9.012*g/mole;
137	density = 1.848*g/cm3;
138	G4Material* Be = new G4Material(name="Beryllium", z=4. , a, density);
139
140	a = 26.98*g/mole;
141	density = 2.7*g/cm3;
142	G4Material* Al = new G4Material(name="Aluminium", z=13., a, density);
143
144	a = 28.09*g/mole;
145	density = 2.33*g/cm3;
146	G4Material* Si = new G4Material(name="Silicon", z=14., a, density);
147
148	G4Element* elH = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
149
150	G4Material* H2O = new G4Material ("Water" , 1.*g/cm3, 2);
151	H2O->AddElement(elH,2);
152	H2O->AddElement(elO,1);
153
154	a = 55.85*g/mole;
155	density = 7.87*g/cm3;
156	G4Material* Fe = new G4Material(name="Iron", z=26., a, density);
157
158	a = 196.97*g/mole;
159	density = 19.32*g/cm3;
160	G4Material* Au = new G4Material(name="Gold", z=79., a, density);
161
162	a = 207.19*g/mole;
163	density = 11.35*g/cm3;
164	G4Material* Pb = new G4Material(name="Lead", z=82., a, density);
165
166	a = 0. ;
167	density = 0. ;
168	G4Material* Vac= new G4Material(name="Vacuum",z=0., a, density,kVacuum);
169
170	const G4MaterialTable* theMaterialTable ;
171	G4Material* apttoMaterial ;
172	G4String MaterialName ;
173	G4Timer theTimer ;
174
175	//--------- Particle definition ---------
176
177	G4ParticleDefinition* theGamma = G4Gamma::GammaDefinition();
178	G4ParticleDefinition* theElectron = G4Electron::ElectronDefinition();
179	G4ParticleDefinition* thePositron = G4Positron::PositronDefinition();
180	G4ParticleDefinition* theProton = G4Proton::ProtonDefinition();
181	G4ParticleDefinition* theAntiProton = G4AntiProton::AntiProtonDefinition();
182	G4ParticleDefinition* thePionPlus = G4PionPlus::PionPlusDefinition();
183	G4ParticleDefinition* thePionMinus = G4PionMinus::PionMinusDefinition();
184	G4ParticleDefinition* theKaonPlus = G4KaonPlus::KaonPlusDefinition();
185	G4ParticleDefinition* theKaonMinus = G4KaonMinus::KaonMinusDefinition();
186
187	G4double* GammaKineticEnergyCuts ;
188	G4double* ElectronKineticEnergyCuts ;
189	G4double* PositronKineticEnergyCuts ;
190	G4double* ParticleKineticEnergyCuts ;
191
192	theMaterialTable = G4Material::GetMaterialTable() ;
193
194	G4double cutinrange,CutInRangeele,CutInRangepos ;
195
196	G4eIonisation theElectronIonisation , thePositronIonisation ;
197	G4ProcessManager* theElectronProcessManager = theElectron->GetProcessManager();
198	theElectronProcessManager->AddProcess(&theElectronIonisation,-1,0,0) ;
199	G4ProcessManager* thePositronProcessManager = thePositron->GetProcessManager();
200	thePositronProcessManager->AddProcess(&thePositronIonisation,-1,0,0) ;
201
202	G4ParticleWithCuts* theParticle ;
203
204	G4double energy, momentum, mass;
205	G4ProcessVector* palongget ;
206	G4ProcessVector* palongdo ;
207	G4ProcessVector* ppostget ;
208	G4ProcessVector* ppostdo ;
209
210	G4String confirm ;
211
212	G4cout << " Do you want the proton as particle (yes/no)? " << std::flush;
213	G4cin >> confirm ;
214	if(confirm == "yes")
215	{
216	mass=theProton->GetPDGMass();
217	theParticle = theProton;
218	}
219	else
220	{
221	G4cout << " Do you want the antiproton as particle (yes/no)? " << std::flush;
222	G4cin >> confirm ;
223	if(confirm == "yes")
224	{
225	mass=theAntiProton->GetPDGMass();
226	theParticle = theAntiProton;
227	}
228	else
229	{
230	G4cout << " Do you want the pi+ as particle (yes/no)? " << std::flush;
231	G4cin >> confirm ;
232	if(confirm == "yes")
233	{
234	mass=thePionPlus->GetPDGMass();
235	theParticle = thePionPlus;
236	}
237	else
238	{
239	G4cout << " Do you want the pi- as particle (yes/no)? " << std::flush;
240	G4cin >> confirm ;
241	if(confirm == "yes")
242	{
243	mass=thePionMinus->GetPDGMass();
244	theParticle = thePionMinus;
245	}
246	else
247	{
248	G4cout << " Do you want the K+ as particle (yes/no)? " << std::flush;
249	G4cin >> confirm ;
250	if(confirm == "yes")
251	{
252	mass=theKaonPlus->GetPDGMass();
253	theParticle = theKaonPlus;
254	}
255	else
256	{
257	G4cout << " Do you want the K- as particle (yes/no)? " << std::flush;
258	G4cin >> confirm ;
259	if(confirm == "yes")
260	{
261	mass=theKaonMinus->GetPDGMass();
262	theParticle = theKaonMinus;
263	}
264	else
265	{
266	G4cout << " There is no other particle in the test." << G4endl;
267	return EXIT_FAILURE;
268	}
269	}
270	}
271	}
272	}
273	}
274
275
276	energy = 1000.0GeV + mass ; // was 1.0GeV now 1.0*TeV
277	momentum=std::sqrt(energyenergy-massmass) ;
278
279	G4ParticleMomentum theMomentum(momentum,0.,0.);
280
281	G4double pModule = theMomentum.mag();
282
283	G4DynamicParticle aParticle(theParticle,energy,theMomentum);
284
285	aParticle.SetKineticEnergy(energy-mass);
286
287
288	// G4hIonisation theParticleIonisation ;
289
290	G4PAIonisation theParticleIonisation ;
291
292	G4ProcessManager* theParticleProcessManager = theParticle->GetProcessManager();
293	theParticleProcessManager->AddProcess(&theParticleIonisation,-1,0,0) ;
294
295	G4ForceCondition cond ;
296	G4ForceCondition* condition = &cond ;
297
298	G4double currentSafety ;
299	G4double& refsafety=currentSafety;
300
301
302	G4cout << "cut for GAMMA in mm =" ;
303	G4cin >> cutinrange ;
304	theGamma->SetCuts(cutinrange) ;
305	G4cout << "gamma,cut in range(mm)=" << theGamma->GetCuts() << G4endl ;
306	outFile << " ---------------------------------------" << G4endl ;
307	outFile << " gamma,cut in range(mm)=" << theGamma->GetCuts() << G4endl ;
308
309	GammaKineticEnergyCuts = theGamma->GetCutsInEnergy() ;
310	for (G4int icut=0; icut<theMaterialTable->length(); icut++)
311	{
312	G4cout << "material index=" << icut << " kin.energy cut(MeV)=" <<
313	GammaKineticEnergyCuts[icut] << G4endl ;
314	outFile << " material index=" << icut << " kin.energy cut(MeV)=" <<
315	GammaKineticEnergyCuts[icut] << G4endl ;
316	}
317
318
319	G4cout << "cut for ELECTRON in mm =" ;
320	G4cin >> cutinrange ;
321	CutInRangeele = cutinrange ;
322	theElectron->SetCuts(cutinrange) ;
323	G4cout << "electron,cut in range(mm)=" << theElectron->GetCuts() << G4endl ;
324	outFile << " ---------------------------------------" << G4endl ;
325	outFile << " electron,cut in range(mm)=" << theElectron->GetCuts() << G4endl ;
326
327	ElectronKineticEnergyCuts = theElectron->GetCutsInEnergy() ;
328	for ( icut=0; icut<theMaterialTable->length(); icut++)
329	{
330	G4cout << "material index=" << icut << " kin.energy cut(MeV)=" <<
331	ElectronKineticEnergyCuts[icut] << G4endl ;
332	outFile << " material index=" << icut << " kin.energy cut(MeV)=" <<
333	ElectronKineticEnergyCuts[icut] << G4endl ;
334	}
335
336	G4cout << "cut for POSITRON in mm =" ;
337	G4cin >> cutinrange ;
338	CutInRangepos = cutinrange ;
339	thePositron->SetCuts(cutinrange) ;
340	G4cout << "positron,cut in range(mm)=" << thePositron->GetCuts() << G4endl ;
341	outFile << " ---------------------------------------" << G4endl ;
342	outFile << " positron,cut in range(mm)=" << thePositron->GetCuts() << G4endl ;
343
344	PositronKineticEnergyCuts = thePositron->GetCutsInEnergy() ;
345	for ( icut=0; icut<theMaterialTable->length(); icut++)
346	{
347	G4cout << "material index=" << icut << " kin.energy cut(MeV)=" <<
348	PositronKineticEnergyCuts[icut] << G4endl ;
349	outFile << " material index=" << icut << " kin.energy cut(MeV)=" <<
350	PositronKineticEnergyCuts[icut] << G4endl ;
351	}
352
353	G4cout << "cut for hadrons in mm =" ;
354	G4cin >> cutinrange ;
355	theParticle->SetCuts(cutinrange) ;
356	G4cout << "after particle setcuts " << G4endl;
357
358
359	G4cout << "cut in range(mm)=" << theParticle->GetLengthCuts() << G4endl ;
360	outFile << " ---------------------------------------" << G4endl ;
361	outFile << " cut in range(mm)=" << theParticle->GetLengthCuts() << G4endl ;
362
363	ParticleKineticEnergyCuts = theParticle->GetEnergyCuts() ;
364
365	for ( icut=0; icut<theMaterialTable->length(); icut++)
366	{
367	G4cout << "material index=" << icut << " kin.energy cut(MeV)=" <<
368	ParticleKineticEnergyCuts[icut] << G4endl ;
369	outFile << " material index=" << icut << " kin.energy cut(MeV)=" <<
370	ParticleKineticEnergyCuts[icut] << G4endl ;
371	}
372
373	G4cout << " ------ ----- " << G4endl ;
374	outFile << " " << G4endl;
375
376	/* *********************************************************************
377	// Output to file of fPAItransferBank data
378
379	G4bool isOutRange ;
380	G4PhysicsLogVector*
381	aLogVector = new G4PhysicsLogVector(G4PAIonisation::GetMinKineticEnergy(),
382	G4PAIonisation::GetMaxKineticEnergy(),
383	G4PAIonisation::GetBinNumber() ) ;
384
385	for(G4int materialIndex=0 ;
386	materialIndex < theMaterialTable->length() ; materialIndex++)
387	{
388	apttoMaterial = (*theMaterialTable)[materialIndex] ;
389	MaterialName = apttoMaterial->GetName() ;
390	outFile<<"PAI transfer data for the material = "
391	<<MaterialName<<G4endl<<G4endl ;
392
393	outFile<<"Particle Tkin"<<"\t\t"
394	<<"Transfer energy, keV"<<"\t"
395	<<"Primary ionisation, 1/mm"<<G4endl<<G4endl ;
396
397	// G4cout<<"Read of TotBin ? TotBin = "<<G4PAIonisation::GetBinNumber()<<G4endl;
398
399	for(G4int iTkin=48;iTkin<G4PAIonisation::GetBinNumber();iTkin++)
400	{
401
402	outFile<<aLogVector->GetLowEdgeEnergy(iTkin)<<G4endl<<G4endl ;
403
404	G4PhysicsVector* trVec = (*G4PAIenergyLoss::GetPAItransferBank())
405	(materialIndex*G4PAIonisation::GetBinNumber() + iTkin) ;
406	G4cout<<materialIndex*G4PAIonisation::GetBinNumber() + iTkin
407	<<"\t\t"<<trVec->GetVectorLength()<<G4endl ;
408	for(G4int iTr=0;iTr<trVec->GetVectorLength();iTr++)
409	{
410	outFile<<"\t\t"<<trVec->GetLowEdgeEnergy(iTr)*1000.
411	<<"\t\t"<<(*trVec)(iTr)<<G4endl ;
412	}
413
414	}
415	outFile<<G4endl ;
416	}
417	**************************************************************** */
418
419
420	G4cout<<"Start dE/dx distribution"<<G4endl ;
421	G4int iLoss, iStat, iStatMax = 10000;
422	G4double energyLoss[200], delta ;
423	G4int spectrum[200] ;
424	const G4DynamicParticle* particlePtr = &aParticle ;
425
426	for(iLoss=0;iLoss<200;iLoss++)
427	{
428	energyLoss[iLoss] = 0.5iLosskeV ;
429	spectrum[iLoss] = 0 ;
430	}
431	for(iStat=0;iStat<iStatMax;iStat++)
432	{
433	delta = G4PAIenergyLoss::GetLossWithFluct(10.0*mm,particlePtr,Xe) ;
434	for(iLoss=0;iLoss<200;iLoss++)
435	{
436	if(delta <= energyLoss[iLoss]) break ;
437	}
438	spectrum[iLoss-1]++ ;
439	}
440	G4double meanLoss = 0.0 ;
441	outFile<<"Energy loss, keV"<<"\t\t"<<"Distribution"<<G4endl<<G4endl ;
442	for(iLoss=0;iLoss<200;iLoss++)
443	{
444	outFile<<energyLoss[iLoss]*1000.0<<"\t\t"<<spectrum[iLoss]<<G4endl ;
445	meanLoss +=energyLoss[iLoss]*spectrum[iLoss] ;
446	}
447	G4cout<<"Mean loss over spectrum = "<<meanLoss*1000./iStatMax<<" keV"<<G4endl ;
448
449	G4int message ;
450	G4cout<<"Continue ?(enter int>0) Break ? (enter int<=0)"<<G4endl ;
451	G4cin>>message ;
452	if(message<= 0) // break program at this point
453	{
454	return 1 ;
455	}
456
457
458
459	outFile << " ionisation test **************************************" << G4endl ;
460	outFile << " " << G4endl;
461	outFile << " particle = " <<
462	aParticle.GetDefinition()->GetParticleName() << G4endl ;
463	outFile << G4endl;
464
465	palongget = aParticle.GetDefinition()->GetProcessManager()
466	->GetAlongStepProcessVector(0);
467	ppostget = aParticle.GetDefinition()->GetProcessManager()
468	->GetPostStepProcessVector(0);
469	palongdo = aParticle.GetDefinition()->GetProcessManager()
470	->GetAlongStepProcessVector(1);
471	ppostdo = aParticle.GetDefinition()->GetProcessManager()
472	->GetPostStepProcessVector(1);
473
474	//---------------------------------- Physics --------------------------------
475
476	G4int itry=1, Ntry=1, Nstart, ir;
477	G4double r ;
478
479	//**************************************************************************
480	const G4int Nbin=97 ;
481	G4double TkinMeV[Nbin] =
482	{0.001,0.0015,0.002,0.003,0.004,0.005,0.006,0.008,
483	0.01,0.015,0.02,0.03,0.04,0.05,0.06,0.08,
484	0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.8,
485	1.,1.5,2.,3.,4.,5.,6.,8.,
486	10.,15.,20.,30.,40.,50.,60.,80.,
487	100.,150.,200.,300.,400.,500.,600.,800.,
488	1.0e3,1.5e3,2.0e3,3.0e3,4.0e3,5.0e3,6.0e3,8.0e3,
489	1.0e4,1.5e4,2.0e4,3.0e4,4.0e4,5.0e4,6.0e4,8.0e4,
490	1.0e5,1.5e5,2.0e5,3.0e5,4.0e5,5.0e5,6.0e5,8.0e5,
491	1.0e6,1.5e6,2.0e6,3.0e6,4.0e6,5.0e6,6.0e6,8.0e6,
492	1.0e7,1.5e7,2.0e7,3.0e7,4.0e7,5.0e7,6.0e7,8.0e7,
493	1.0e8,1.5e8,2.0e8,3.0e8,4.0e8,5.0e8,6.0e8,8.0e8,
494	1.0e9} ;
495
496	G4int J=-1 ;
497
498	G4double lambda,trueStep,geomStep,stepLimit,
499	previousStepSize,currentMinimumStep ;
500	G4ParticleChange* aParticleChange ;
501
502	G4double T,dEdx,range ;
503
504	NEXTMATERIAL: ;
505	J = J+1 ;
506	if ( J >= theMaterialTable->length() )
507	{ G4cout << "that was the last material in the table --> STOP" << G4endl;
508	return EXIT_FAILURE ; }
509
510	apttoMaterial = (*theMaterialTable)[ J ] ;
511	MaterialName = apttoMaterial->GetName() ;
512	G4cout << "material=" << MaterialName << G4endl ;
513	G4cout << "Do you want the Energyloss test 1. for this material?" << G4endl ;
514	G4cout << "type a positive number if the answer is YES" << G4endl ;
515	G4cout << "type a negative number if the answer is NO " << G4endl ;
516	G4int icont ;
517	G4cin >> icont ;
518	if ( icont < 0 )
519	goto NEXTMATERIAL ;
520
521	//---------- Volume definition ---------------------
522
523	G4VPhysicalVolume* myVolume ;
524
525	myVolume = BuildVolume(apttoMaterial) ;
526
527	//--------- track and Step definition (for this test ONLY!)------------
528	G4ThreeVector aPosition(0.,0.,0.);
529	const G4ThreeVector aDirection(0.,0.,1.) ;
530	G4double aTime = 0. ;
531
532	G4Track* tracke = new G4Track(&aParticle,aTime,aPosition) ;
533	G4Track& trackele = (*tracke) ;
534	//(*tracke).SetVolume(myVolume) ;
535	G4GRSVolume* touche = new G4GRSVolume(myVolume, NULL, aPosition);
536	(*tracke).SetTouchable(touche);
537	(*tracke).SetMomentumDirection(aDirection) ;
538
539
540	G4Step* Step = new G4Step() ;
541	G4Step& Step = (*Step) ;
542
543	G4StepPoint* aPoint = new G4StepPoint();
544	(*aPoint).SetPosition(aPosition) ;
545	G4double safety = 10000.*cm ;
546	(*aPoint).SetSafety(safety) ;
547
548	(*Step).SetPostStepPoint(aPoint) ;
549
550	//**************************************************************************
551
552	G4cout << G4endl;
553	G4cout <<" " << MaterialName << " Energyloss test 1." << G4endl ;
554	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
555	G4cout << G4endl ;
556	G4cout << "kin.en.(MeV) dE/dx(MeV/mm) range(mm) Step(mm)" << G4endl ;
557	G4cout << G4endl ;
558
559	outFile << G4endl;
560	outFile <<" " << MaterialName << " Energyloss test 1." << G4endl ;
561	outFile << " +++++++++++++++++++++++++++++++++++++++++" << G4endl ;
562	outFile << G4endl ;
563	outFile << "kin.en.(MeV) dE/dx(MeV/mm) range(mm) Step(mm)" << G4endl ;
564	outFile << G4endl ;
565
566
567	for ( G4int i=0 ; i<Nbin ; i++)
568	{
569	trueStep = cutinrange ;
570	previousStepSize = cutinrange ;
571	currentMinimumStep = trueStep ;
572	(*tracke).SetKineticEnergy(TkinMeV[i]) ;
573	stepLimit = (*palongget)(0)->AlongStepGetPhysicalInteractionLength(
574	trackele,
575	previousStepSize,
576	currentMinimumStep,
577	refsafety) ;
578
579	dEdx = theParticleIonisation.GetdEdx() ;
580	range = theParticleIonisation.GetRangeNow() ;
581	T = TkinMeV[i] ;
582
583	G4cout <<" " << T << " " << dEdx/mm << " " ;
584	G4cout << range/mm << " " << stepLimit/mm << G4endl ;
585
586	outFile <<" " << T << " " << dEdx/mm << " " ;
587	outFile << range/mm << " " << stepLimit/mm << G4endl ;
588
589	}
590
591	G4cout << G4endl;
592	outFile << G4endl;
593
594	ENERGYLOSS2: ;
595
596	G4cout << "material=" << MaterialName << G4endl ;
597	G4cout << "Do you want the Energyloss test 2. for this material?" << G4endl ;
598	G4cout << "type a positive number if the answer is YES" << G4endl ;
599	G4cout << "type a negative number if the answer is NO " << G4endl ;
600	G4cin >> icont ;
601	if ( icont < 0 )
602	goto ENERGYLOSS3 ;
603
604	G4double TMeV,stepmm,stepmx,meanloss,lossnow ;
605
606
607	G4cout << "give an energy value in MeV " ;
608	G4cin >> TMeV ;
609
610	trueStep = cutinrange ;
611	previousStepSize = cutinrange ;
612	currentMinimumStep = trueStep ;
613	(*tracke).SetKineticEnergy(TMeV) ;
614	stepmx = (*palongget)(0)->AlongStepGetPhysicalInteractionLength(
615	trackele,
616	previousStepSize,
617	currentMinimumStep,
618	refsafety);
619
620	G4cout << " give a steplength in mm , the max. meaningful Step is " << stepmx << " mm" <<G4endl;
621	G4cout << "Step:" ;
622	G4cin >> stepmm ;
623
624	(*Step).SetTrack(tracke) ;
625	(*Step).SetStepLength(stepmm);
626
627
628	aParticleChange = (*palongdo)(0)->AlongStepDoIt(trackele,Step);
629	meanloss = theParticleIonisation.GetMeanLoss() ;
630	lossnow = TMeV-(*aParticleChange).GetEnergyChange();
631
632	G4cout << G4endl;
633	G4cout <<" " << MaterialName << " Energyloss test 2." << G4endl ;
634	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
635	G4cout << G4endl ;
636	G4cout << "kin.en.(MeV) Step(mm) meanloss(MeV) act.loss(MeV)" << G4endl ;
637	G4cout << TMeV << " " << stepmm << " " << meanloss << " " << lossnow << G4endl ;
638	G4cout << " status change:" << (*aParticleChange).GetStatusChange() << G4endl ;
639	G4cout << G4endl ;
640
641	outFile << G4endl;
642	outFile <<" " << MaterialName << " Energyloss test 2." << G4endl ;
643	outFile << " +++++++++++++++++++++++++++++++++++++++++" << G4endl ;
644	outFile << G4endl ;
645	outFile << "kin.en.(MeV) Step(mm) meanloss(MeV) act.loss(MeV)" << G4endl ;
646	outFile << TMeV << " " << stepmm << " " << meanloss << " " << lossnow << G4endl ;
647	outFile << " status change:" << (*aParticleChange).GetStatusChange() << G4endl ;
648	outFile << G4endl ;
649
650	goto ENERGYLOSS2 ;
651
652	ENERGYLOSS3: ;
653
654	G4cout << "material=" << MaterialName << G4endl ;
655	G4cout << "Do you want the Energyloss test 3. for this material?" << G4endl ;
656	G4cout << "type a positive number if the answer is YES" << G4endl ;
657	G4cout << "type a negative number if the answer is NO " << G4endl ;
658	G4cin >> icont ;
659	if ( icont < 0 )
660	goto DELTARAY1 ;
661
662	G4cout << "give an energy value in MeV " ;
663	G4cin >> TMeV ;
664
665	trueStep = cutinrange ;
666	previousStepSize = cutinrange ;
667	currentMinimumStep = trueStep ;
668	(*tracke).SetKineticEnergy(TMeV) ;
669	stepmx = (*palongget)(0)->AlongStepGetPhysicalInteractionLength(
670	trackele,
671	previousStepSize,
672	currentMinimumStep,
673	refsafety);
674
675	G4cout << " give a steplength in mm , the max. meaningful Step is " << stepmx << " mm" <<G4endl;
676	G4cout << "Step:" ;
677	G4cin >> stepmm ;
678
679	(*Step).SetTrack(tracke) ;
680	(*Step).SetStepLength(stepmm);
681
682
683	G4cout << " give number of events you want " ;
684	G4int nbev,ibev ;
685	G4cin >> nbev ;
686
687	meanloss=0.;
688	theTimer.Start();
689
690	for ( ibev=0; ibev<nbev; ibev++)
691	{
692	aParticleChange = (*palongdo)(0)->AlongStepDoIt(trackele,Step);
693	lossnow = TMeV-(*aParticleChange).GetEnergyChange();
694
695	meanloss += lossnow ;
696	}
697
698	theTimer.Stop();
699	meanloss /= nbev ;
700	G4cout << G4endl;
701	G4cout <<" " << MaterialName << " Energyloss test 3." << G4endl ;
702	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
703	G4cout << G4endl ;
704	G4cout << "kin.en.(MeV) Step(mm) meanloss(MeV) time/event(sec) " << G4endl ;
705	G4cout << TMeV << " " << stepmm << " " << meanloss << " " <<
706	theTimer.GetUserElapsed()/nbev << G4endl ;
707	G4cout << G4endl ;
708
709	outFile << G4endl;
710	outFile <<" " << MaterialName << " Energyloss test 3." << G4endl ;
711	outFile << " +++++++++++++++++++++++++++++++++++++++++" << G4endl ;
712	outFile << G4endl ;
713	outFile << "kin.en.(MeV) Step(mm) meanloss(MeV) time/event(sec) " << G4endl ;
714	outFile << TMeV << " " << stepmm << " " << meanloss << " " <<
715	theTimer.GetUserElapsed()/nbev << G4endl ;
716	outFile << G4endl ;
717
718	goto ENERGYLOSS3 ;
719
720	DELTARAY1: ;
721	G4cout << "material=" << MaterialName << G4endl ;
722	G4cout << "Do you want the delta ray test 1. for this material?" << G4endl ;
723	G4cout << "type a positive number if the answer is YES" << G4endl ;
724	G4cout << "type a negative number if the answer is NO " << G4endl ;
725	G4cin >> icont ;
726	if ( icont < 0 )
727	goto DELTARAY2 ;
728
729
730	G4cout << G4endl;
731	G4cout <<" " << MaterialName << " delta ray test 1." << G4endl ;
732	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
733	G4cout << G4endl ;
734	G4cout << "kin.en.(MeV) mean free path(mm)" << G4endl ;
735	G4cout << G4endl ;
736
737	outFile << G4endl;
738	outFile <<" " << MaterialName << " delta ray test 1." << G4endl ;
739	outFile << " +++++++++++++++++++++++++++++++++++++++++" << G4endl ;
740	outFile << G4endl ;
741	outFile << "kin.en.(MeV) mean free path(mm)" << G4endl ;
742	outFile << G4endl ;
743
744	for ( i=0 ; i<Nbin ; i++)
745	{
746
747	previousStepSize = cutinrange ;
748	(*tracke).SetKineticEnergy(TkinMeV[i]) ;
749	stepLimit = theParticleIonisation.GetMeanFreePath(
750	trackele,
751	previousStepSize,
752	condition) ;
753
754	T = TkinMeV[i] ;
755
756	G4cout <<" " << T << " " << stepLimit/mm << G4endl ;
757
758	outFile <<" " << T << " " << stepLimit/mm << G4endl ;
759
760	}
761
762	G4cout << G4endl;
763	outFile << G4endl;
764
765	//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
766	DELTARAY2: ;
767
768	G4cout << "material=" << MaterialName << G4endl ;
769	G4cout << "Do you want the deltaray test 2. for this material?" << G4endl ;
770	G4cout << "type a positive number if the answer is YES" << G4endl ;
771	G4cout << "type a negative number if the answer is NO " << G4endl ;
772	G4cin >> icont ;
773
774	G4double newenergy,dx,dy,dz,Tdelta,ddx,ddy,ddz ;
775	G4int nd ;
776	const G4ThreeVector* momdir ;
777	G4ParticleMomentum ddir ;
778
779	if ( icont < 0 )
780	goto BREMS1 ;
781
782	G4cout << "give an energy value in MeV " ;
783	G4cin >> TMeV ;
784
785	stepmm = 1. ;
786
787	(*Step).SetTrack(tracke) ;
788	(*Step).SetStepLength(stepmm);
789
790
791	(*tracke).SetKineticEnergy(TMeV) ;
792	aParticleChange = (*ppostdo)(0)->PostStepDoIt(trackele,Step);
793
794	newenergy=(*aParticleChange).GetEnergyChange() ;
795	momdir=(*aParticleChange).GetMomentumChange();
796	dx = (*momdir).x();
797	dy = (*momdir).y();
798	dz = (*momdir).z();
799	nd=aParticleChange->GetNumberOfSecondaries();
800
801	if(nd>0)
802	{
803	Tdelta=aParticleChange->GetSecondary(0)->GetKineticEnergy();
804	ddir=aParticleChange->GetSecondary(0)->
805	GetMomentumDirection();
806	ddx = (ddir).x();
807	ddy = (ddir).y();
808	ddz = (ddir).z();
809	}
810
811	G4cout << G4endl;
812	G4cout <<" " << MaterialName << " delta ray test 2." << G4endl ;
813	G4cout << " ++++++++++++++++++++++++++++++++++++++++++++" << G4endl ;
814	G4cout << G4endl ;
815	G4cout << "T=" << TMeV << " newT=" << newenergy << " (MeV)" << G4endl ;
816	G4cout << " status change:" << (*aParticleChange).GetStatusChange() << G4endl ;
817	if(nd>0)
818	G4cout << "Tdelta=" << Tdelta << G4endl ;
819	G4cout << "new direction:" << dx << " " << dy << " " << dz << G4endl;
820	if(nd>0)
821	G4cout << "delta direction:" << ddx << " " << ddy << " " << ddz << G4endl ;
822	G4cout << G4endl ;
823
824	outFile << G4endl;
825	outFile <<" " << MaterialName << " delta ray test 2." << G4endl ;
826	outFile << " +++++++++++++++++++++++++++++++++++++++++" << G4endl ;
827	outFile << G4endl ;
828	outFile << "T=" << TMeV << " newT=" << newenergy << " (MeV)" << G4endl;
829	outFile << " status change:" << (*aParticleChange).GetStatusChange() << G4endl ;
830	if(nd>0)
831	outFile << "Tdelta=" << Tdelta << G4endl ;
832	outFile << "new direction:" << dx << " " << dy << " " << dz << G4endl;
833	if(nd>0)
834	outFile << "delta direction:" << ddx << " " << ddy << " " << ddz << G4endl ;
835	outFile << G4endl ;
836
837	(*aParticleChange).Clear();
838
839	goto DELTARAY2 ;
840
841	BREMS1: ;
842
843	if( J < theMaterialTable->length()-1 )
844	goto NEXTMATERIAL ;
845
846	return EXIT_SUCCESS;
847	}

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

Context Navigation

source: trunk/source/processes/electromagnetic/standard/test/PAIenergyLossTest.cc @ 1347

Download in other formats: