Context Navigation

testG4ForwardXrayTR.cc @ 1347

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

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

Context Navigation

source: trunk/source/processes/electromagnetic/xrays/test/testG4ForwardXrayTR.cc @ 1347

Download in other formats: