Context Navigation

G4PAIdNdxTest.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.0 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: G4PAIdNdxTest.cc,v 1.8 2008/04/14 14:53:39 grichine Exp $
28	// GEANT4 tag $Name: geant4-09-04-ref-00 $
29	//
30	//
31	//
32	//
33	//
34	// Test routine for G4PAIxSection class code
35	//
36	// History:
37	//
38	// 07.02.02, V. Grichine implementation based on G4PAIxSectionTest
39
40	#include "G4ios.hh"
41	#include <fstream>
42	#include <cmath>
43	#include "globals.hh"
44	#include "Randomize.hh"
45
46	#include "G4Isotope.hh"
47	#include "G4Element.hh"
48	#include "G4Material.hh"
49	#include "G4MaterialTable.hh"
50	#include "G4SandiaTable.hh"
51
52	// #include "G4PAIonisation.hh"
53	#include "G4PAIxSection.hh"
54
55	int main()
56	{
57	std::ofstream outFile("PAIdNdx.out", std::ios::out ) ;
58	outFile.setf( std::ios::scientific, std::ios::floatfield );
59
60	std::ofstream fileOut("PAICrenkovPlasmon.out", std::ios::out ) ;
61	fileOut.setf( std::ios::scientific, std::ios::floatfield );
62
63	// std::ifstream fileRead("exp.dat", std::ios::out ) ;
64	// fileRead.setf( std::ios::scientific, std::ios::floatfield );
65
66	std::ofstream fileWrite("exp.dat", std::ios::out ) ;
67	fileWrite.setf( std::ios::scientific, std::ios::floatfield );
68
69	std::ofstream fileWrite1("mprrpai.dat", std::ios::out ) ;
70	fileWrite1.setf( std::ios::scientific, std::ios::floatfield );
71
72	// Create materials
73
74
75	G4int iz , n, nel, ncomponents ;
76	G4double a, z, ez, density , temperature, pressure, fractionmass ;
77	G4State state ;
78	G4String name, symbol ;
79
80	// G4Element* elH = new G4Element ("Hydrogen", "H", 1. , 1.01*g/mole);
81
82	a = 14.01*g/mole;
83	G4Element* elN = new G4Element(name="Nitrogen", symbol="N", ez=7., a);
84
85	a = 16.00*g/mole;
86	// G4Element* elO = new G4Element(name="Oxigen", symbol="O", ez=8., a);
87
88	a = 12.01*g/mole;
89	G4Element* elC = new G4Element(name="Carbon",symbol="C", ez=6., a);
90
91	a = 55.85*g/mole;
92	G4Element* elFe = new G4Element(name="Iron",symbol="Fe", ez=26., a);
93
94	a = 16.00*g/mole;
95	G4Element* elO = new G4Element(name="Oxygen",symbol="O", ez=8., a);
96
97	a = 1.01*g/mole;
98	G4Isotope* ih1 = new G4Isotope("Hydrogen",iz=1,n=1,a);
99
100	a = 2.01*g/mole;
101	G4Isotope* ih2 = new G4Isotope("Deuterium",iz=1,n=2,a);
102
103	G4Element* elH = new G4Element(name="Hydrogen",symbol="H",2);
104	elH->AddIsotope(ih1,.999);
105	elH->AddIsotope(ih2,.001);
106
107	a = 39.948*g/mole;
108	G4Element* elAr = new G4Element(name="Argon", symbol="Ar", z=18., a);
109
110	a = 131.29*g/mole;
111	G4Element* elXe = new G4Element(name="Xenon", symbol="Xe", z=54., a);
112
113	a = 19.00*g/mole;
114	G4Element* elF = new G4Element(name="Fluorine", symbol="F", z=9., a);
115
116	a = 69.723*g/mole;
117	G4Element* elGa = new G4Element(name="Ga", symbol="Ga", z=31., a);
118
119	a = 74.9216*g/mole;
120	G4Element* elAs = new G4Element(name="As", symbol="As", z=33., a);
121
122
123	// G4Isotope::DumpInfo();
124	// G4Element::DumpInfo();
125	// G4Material::DumpInfo();
126
127	/* ***************************************************************
128
129	a = 9.012*g/mole;
130	density = 1.848*g/cm3;
131	G4Material* Be = new G4Material(name="Beryllium", z=4. , a, density);
132
133	density = 1.390*g/cm3;
134	a = 39.95*g/mole;
135	G4Material* lAr = new G4Material(name="liquidArgon", z=18., a, density);
136
137	density = 19.32*g/cm3;
138	a =196.97*g/mole;
139	G4Material* Au = new G4Material(name="Gold" , z=79., a, density);
140
141	// Carbon dioxide
142
143	density = 1.977*mg/cm3;
144	G4Material* CO2 = new G4Material(name="CO2", density, nel=2,
145	kStateGas,273.15kelvin,1.atmosphere);
146	CO2->AddElement(elC,1);
147	CO2->AddElement(elO,2);
148
149	density = 1.290*mg/cm3; // old air from elements
150	G4Material* air = new G4Material(name="air" , density, ncomponents=2);
151	Air->AddElement(elN, fractionmass=0.7);
152	Air->AddElement(elO, fractionmass=0.3);
153
154
155	density = 1.25053*mg/cm3 ; // STP
156	a = 14.01*g/mole ; // get atomic weight !!!
157	// a = 28.016*g/mole;
158	G4Material* newN2 = new G4Material(name="newN2", z= 7.,a,density) ;
159
160	density = 1.25053*mg/cm3 ; // STP
161	G4Material* anotherN2 = new G4Material(name="anotherN2", density,ncomponents=2);
162	anotherN2->AddElement(elN, 1);
163	anotherN2->AddElement(elN, 1);
164
165	density = 1.000*g/cm3;
166	G4Material* H2O = new G4Material(name="Water", density, ncomponents=2);
167	H2O->AddElement(elH, natoms=2);
168	H2O->AddElement(elO, natoms=1);
169
170
171
172	a = 26.98159*g/mole;
173	density = 2.7*g/cm3;
174	G4Material* Al = new G4Material(name="Aluminium", z=13., a, density);
175
176	// Silicon as detector material
177
178	density = 2.330*g/cm3;
179	a = 28.0855*g/mole;
180	G4Material* Si = new G4Material(name="Silicon", z=14., a, density);
181
182
183	density = 7.870*g/cm3;
184	a = 55.85*g/mole;
185	G4Material* Fe = new G4Material(name="Iron" , z=26., a, density);
186
187	density = 8.960*g/cm3;
188	a = 63.55*g/mole;
189	G4Material* Cu = new G4Material(name="Copper" , z=29., a, density);
190
191	density = 11.35*g/cm3;
192	a = 207.19*g/mole;
193	G4Material* Pb = new G4Material(name="Lead" , z=82., a, density);
194
195	// Polypropelene
196
197	G4Material* CH2 = new G4Material ("Polypropelene" , 0.91*g/cm3, 2);
198	CH2->AddElement(elH,2);
199	CH2->AddElement(elC,1);
200
201	// Kapton (polyimide)
202
203	density = 1.39*g/cm3;
204	G4Material* Kapton = new G4Material(name="Kapton", density, nel=3);
205	Kapton->AddElement(elO,2);
206	Kapton->AddElement(elC,5);
207	Kapton->AddElement(elH,4);
208
209
210	// Germanium as detector material
211
212	density = 5.323*g/cm3;
213	a = 72.59*g/mole;
214	G4Material* Ge = new G4Material(name="Ge", z=32., a, density);
215
216	// GaAs detectors
217
218	density = 5.32*g/cm3;
219	G4Material* GaAs = new G4Material(name="GaAs",density, nel=2);
220	GaAs->AddElement(elGa,1);
221	GaAs->AddElement(elAs,1);
222
223	// Diamond detectors
224
225	density = 3.5*g/cm3;
226	G4Material* Diamond = new G4Material(name="Diamond",density, nel=1);
227	Diamond->AddElement(elC,1);
228
229	G4double TRT_Xe_density = 5.485*mg/cm3;
230	G4Material* TRT_Xe = new G4Material(name="TRT_Xe", TRT_Xe_density, nel=1,
231	kStateGas,293.15kelvin,1.atmosphere);
232	TRT_Xe->AddElement(elXe,1);
233
234	G4double TRT_CO2_density = 1.842*mg/cm3;
235	G4Material* TRT_CO2 = new G4Material(name="TRT_CO2", TRT_CO2_density, nel=2,
236	kStateGas,293.15kelvin,1.atmosphere);
237	TRT_CO2->AddElement(elC,1);
238	TRT_CO2->AddElement(elO,2);
239
240	G4double TRT_CF4_density = 3.9*mg/cm3;
241	G4Material* TRT_CF4 = new G4Material(name="TRT_CF4", TRT_CF4_density, nel=2,
242	kStateGas,293.15kelvin,1.atmosphere);
243	TRT_CF4->AddElement(elC,1);
244	TRT_CF4->AddElement(elF,4);
245
246	// ATLAS TRT straw tube gas mixture (20 C, 1 atm)
247
248	G4double XeCO2CF4_density = 4.76*mg/cm3;
249	G4Material* XeCO2CF4 = new G4Material(name="XeCO2CF4", XeCO2CF4_density,
250	ncomponents=3,
251	kStateGas,293.15kelvin,1.atmosphere);
252	XeCO2CF4->AddMaterial(TRT_Xe,0.807);
253	XeCO2CF4->AddMaterial(TRT_CO2,0.039);
254	XeCO2CF4->AddMaterial(TRT_CF4,0.154);
255
256	// TRT_CH2
257
258	density = 0.935*g/cm3;
259	G4Material* TRT_CH2 = new G4Material(name="TRT_CH2",density, nel=2);
260	TRT_CH2->AddElement(elC,1);
261	TRT_CH2->AddElement(elH,2);
262
263	// Radiator
264
265	density = 0.059*g/cm3;
266	G4Material* Radiator = new G4Material(name="Radiator",density, nel=2);
267	Radiator->AddElement(elC,1);
268	Radiator->AddElement(elH,2);
269
270	// Carbon Fiber
271
272	density = 0.145*g/cm3;
273	G4Material* CarbonFiber = new G4Material(name="CarbonFiber",density, nel=1);
274	CarbonFiber->AddElement(elC,1);
275
276
277	// Dry air (average composition)
278
279
280	density = 1.25053*mg/cm3 ; // STP
281	G4Material* Nitrogen = new G4Material(name="N2" , density, ncomponents=1);
282	Nitrogen->AddElement(elN, 2);
283
284	density = 1.4289*mg/cm3 ; // STP
285	G4Material* Oxygen = new G4Material(name="O2" , density, ncomponents=1);
286	Oxygen->AddElement(elO, 2);
287
288	density = 1.7836*mg/cm3 ; // STP
289	G4Material* Argon = new G4Material(name="Argon" , density, ncomponents=1);
290	Argon->AddElement(elAr, 1);
291
292	density = 1.2928*mg/cm3 ; // STP
293	G4Material* Air = new G4Material(name="Air" , density, ncomponents=3);
294	Air->AddMaterial( Nitrogen, fractionmass = 0.7557 ) ;
295	Air->AddMaterial( Oxygen, fractionmass = 0.2315 ) ;
296	Air->AddMaterial( Argon, fractionmass = 0.0128 ) ;
297
298	// Xenon as detector gas, STP
299
300	density = 5.858*mg/cm3 ;
301	a = 131.29*g/mole ;
302	G4Material* Xe = new G4Material(name="Xenon",z=54., a, density );
303
304	// Helium as detector gas, STP
305
306	density = 0.178*mg/cm3 ;
307	a = 4.0026*g/mole ;
308	G4Material* He = new G4Material(name="He",z=2., a, density );
309
310
311	// Krypton as detector gas, STP
312
313	density = 3.700*mg/cm3 ;
314	a = 83.80*g/mole ;
315	G4Material* Kr = new G4Material(name="Kr",z=36., a, density );
316
317	****************************************************** */
318
319	// Neon as detector gas, STP
320
321	density = 0.900*mg/cm3 ;
322	a = 20.179*g/mole ;
323	G4Material* Ne = new G4Material(name="Ne",z=10., a, density );
324
325	// Carbone dioxide, CO2 STP
326
327	density = 1.977*mg/cm3 ;
328	G4Material* CarbonDioxide = new G4Material(name="CO2", density, nel=2) ;
329	CarbonDioxide->AddElement(elC,1) ;
330	CarbonDioxide->AddElement(elO,2) ;
331
332	/* *****************************************************
333
334	// Metane, STP
335
336	density = 0.7174*mg/cm3 ;
337	G4Material* metane = new G4Material(name="CH4",density,nel=2) ;
338	metane->AddElement(elC,1) ;
339	metane->AddElement(elH,4) ;
340
341	// Propane, STP
342
343	density = 2.005*mg/cm3 ;
344	G4Material* propane = new G4Material(name="C3H8",density,nel=2) ;
345	propane->AddElement(elC,3) ;
346	propane->AddElement(elH,8) ;
347
348	// iso-Butane (methylpropane), STP
349
350	density = 2.67*mg/cm3 ;
351	G4Material* isobutane = new G4Material(name="isoC4H10",density,nel=2) ;
352	isobutane->AddElement(elC,4) ;
353	isobutane->AddElement(elH,10) ;
354
355	// 87.5% Xe + 7.5% CH4 + 5% C3H8, 20 C, 1 atm
356
357	density = 4.9196*mg/cm3 ;
358
359	G4Material* XeCH4C3H8 = new G4Material(name="XeCH4C3H8" , density,
360	ncomponents=3);
361	XeCH4C3H8->AddMaterial( Xe, fractionmass = 0.971 ) ;
362	XeCH4C3H8->AddMaterial( metane, fractionmass = 0.010 ) ;
363	XeCH4C3H8->AddMaterial( propane, fractionmass = 0.019 ) ;
364
365	// Propane in MWPC, 2 atm, 20 C
366
367	// density = 3.758*mg/cm3 ;
368	density = 3.736*mg/cm3 ;
369	G4Material* propaneDet = new G4Material(name="detC3H8",density,nel=2) ;
370	propaneDet->AddElement(elC,3) ;
371	propaneDet->AddElement(elH,8) ;
372
373	************************************************** */
374
375	// 90% Ne + 10% CO2, STP
376
377	density = 1.0077*mg/cm3 ;
378	G4Material* Ne10CO2 = new G4Material(name="Ne10CO2" , density,
379
380	ncomponents=2);
381	Ne10CO2->AddMaterial( Ne, fractionmass = 0.8038 ) ;
382	Ne10CO2->AddMaterial( CarbonDioxide, fractionmass = 0.1962 ) ;
383
384	/* *****************************************************
385
386	// 80% Ar + 20% CO2, STP
387
388	density = 1.8223*mg/cm3 ;
389	G4Material* Ar20CO2 = new G4Material(name="Ar20CO2" , density,
390	ncomponents=2);
391	Ar20CO2->AddMaterial( Argon, fractionmass = 0.783 ) ;
392	Ar20CO2->AddMaterial( CarbonDioxide, fractionmass = 0.217 ) ;
393
394	// 93% Ar + 7% CH4, STP
395
396	density = 1.709*mg/cm3 ;
397	G4Material* Ar7CH4 = new G4Material(name="Ar7CH4" , density,
398	ncomponents=2);
399	Ar7CH4->AddMaterial( Argon, fractionmass = 0.971 ) ;
400	Ar7CH4->AddMaterial( metane, fractionmass = 0.029 ) ;
401
402	// 80% Xe + 20% CO2, STP
403
404	density = 5.0818*mg/cm3 ;
405	G4Material* Xe20CO2 = new G4Material(name="Xe20CO2" , density,
406	ncomponents=2);
407	Xe20CO2->AddMaterial( Xe, fractionmass = 0.922 ) ;
408	Xe20CO2->AddMaterial( CarbonDioxide, fractionmass = 0.078 ) ;
409
410	// 80% Kr + 20% CO2, STP
411
412	density = 3.601*mg/cm3 ;
413	G4Material* Kr20CO2 = new G4Material(name="Kr20CO2" , density,
414	ncomponents=2);
415	Kr20CO2->AddMaterial( Kr, fractionmass = 0.89 ) ;
416	Kr20CO2->AddMaterial( CarbonDioxide, fractionmass = 0.11 ) ;
417
418	// 80% He + 20% CO2, STP
419
420	density = 0.5378*mg/cm3 ;
421	G4Material* He20CO2 = new G4Material(name="He20CO2" , density,
422	ncomponents=2);
423	He20CO2->AddMaterial( He, fractionmass = 0.265 ) ;
424	He20CO2->AddMaterial( CarbonDioxide, fractionmass = 0.735 ) ;
425
426	*/ //////////////////////
427
428
429	// G4cout << *(G4Material::GetMaterialTable()) << G4endl;
430
431	//
432	// Create Sandia/PAI tables for given material
433	//
434
435	G4int i, j, k, numOfMaterials, iSan, nbOfElements, sanIndex, row ;
436	G4double maxEnergyTransfer, kineticEnergy ;
437	G4double tau, gamma, bg2, beta2, rateMass, Tmax, Tmin, Tkin ;
438
439	const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable() ;
440
441	numOfMaterials = theMaterialTable->size();
442
443	G4cout<<"Available materials under test : "<< G4endl<<G4endl ;
444	outFile<<"Available materials under test : "<< G4endl<<G4endl ;
445
446	for(k=0;k<numOfMaterials;k++)
447	{
448	G4cout <<k<<"\t"<< " Material : " <<(*theMaterialTable)[k]->GetName() << G4endl ;
449	outFile <<k<<"\t"<< " Material : " <<(*theMaterialTable)[k]->GetName() << G4endl ;
450	}
451	G4String testName ;
452	G4cout<<"Enter material name for test : "<<std::flush ;
453	// G4cin>>testName ;
454
455	for(k=0;k<numOfMaterials;k++)
456	{
457	// if((*theMaterialTable)[k]->GetName() != testName) continue ;
458
459	outFile << "Material : " <<(*theMaterialTable)[k]->GetName() << G4endl ;
460	G4cout << "Material : " <<(*theMaterialTable)[k]->GetName() << G4endl ;
461
462	nbOfElements = (*theMaterialTable)[k]->GetNumberOfElements() ;
463
464	G4cout<<"Sandia cof according old PAI stuff"<<G4endl<<G4endl ;
465	outFile<<"Sandia cof according old PAI stuff"<<G4endl<<G4endl ;
466
467	G4int* thisMaterialZ = new G4int[nbOfElements] ;
468	for(iSan=0;iSan<nbOfElements;iSan++)
469	{
470	thisMaterialZ[iSan] = (G4int)(*theMaterialTable)[k]->
471	GetElement(iSan)->GetZ() ;
472	}
473	G4SandiaTable sandia(k) ;
474	sanIndex = sandia.SandiaIntervals(thisMaterialZ,nbOfElements) ;
475	sanIndex = sandia.SandiaMixing( thisMaterialZ ,
476	(*theMaterialTable)[k]->GetFractionVector() ,
477	nbOfElements,sanIndex) ;
478
479	for(row=0;row<sanIndex-1;row++)
480	{
481	G4cout<<row+1<<"\t"<<sandia.GetPhotoAbsorpCof(row+1,0)/keV ;
482	outFile<<row+1<<" "<<sandia.GetPhotoAbsorpCof(row+1,0)/keV ;
483
484	for(iSan=1;iSan<5;iSan++)
485	{
486	G4cout<<"\t"<<sandia.GetPhotoAbsorpCof(row+1,iSan) ;
487	// (theMaterialTable)[k]->GetDensity() ;
488
489	outFile<<" "<<sandia.GetPhotoAbsorpCof(row+1,iSan) ;
490	// (theMaterialTable)[k]->GetDensity() ;
491	}
492	G4cout<<G4endl ;
493	outFile<<G4endl ;
494	}
495	G4cout<<G4endl ;
496	outFile<<G4endl ;
497
498
499	outFile<<G4endl ;
500	maxEnergyTransfer = 100*keV ;
501	gamma = 4.0 ;
502	bg2 = gamma*gamma - 1 ;
503
504	G4PAIxSection testPAI(k,maxEnergyTransfer,bg2) ;
505
506	G4cout<<"Interval no."<<"\t"<<"Energy interval"<<G4endl<<G4endl ;
507	outFile<<"Interval no."<<"\t"<<"Energy interval"<<G4endl<<G4endl ;
508
509	for(j=1;j<=testPAI.GetIntervalNumber();j++)
510	{
511	G4cout<<j<<"\t\t"<<testPAI.GetEnergyInterval(j)/keV<<G4endl ;
512	outFile<<j<<"\t\t"<<testPAI.GetEnergyInterval(j)/keV<<G4endl ;
513	}
514	G4cout<<G4endl ;
515	outFile<<G4endl ;
516
517	outFile<<"Actual spline size = "<<testPAI.GetSplineSize()<<G4endl ;
518	outFile<<"Normalization Cof = "<<testPAI.GetNormalizationCof()<<G4endl ;
519	outFile<<G4endl ;
520
521	G4cout << "Actual spline size = "<<testPAI.GetSplineSize()<<G4endl ;
522	G4cout <<"Normalization Cof = "<<testPAI.GetNormalizationCof()<<G4endl ;
523	G4cout << G4endl ;
524
525	Tmin = sandia.GetPhotoAbsorpCof(1,0) ; // 0.02*keV ;
526	G4cout<<"Tmin = "<<Tmin/eV<<" eV"<<G4endl;
527
528	outFile
529	// <<"Tkin, keV"<<"\t"
530	<<"gamma"<<"\t\t"
531	// <<"Max E transfer, kev"<<"\t"
532	<<"<dN/dxC>, 1/cm"<<"\t"
533	<< "<dN/dxP>, 1/cm"<<"\t"
534	<<"<dN/dxC+dN/dxP>"<<"\t"
535	<<"<dN/dx>, 1/cm"<<G4endl<<G4endl ;
536
537	G4cout
538	// <<"Tkin, keV"<<"\t"
539	<< "gamma"<<"\t\t"
540	// <<"Max E transfer, kev"<<"\t"
541	<< "<dN/dxC>, 1/cm"<<"\t"
542	<< "<dN/dxP>, 1/cm"<<"\t"
543	<< "<dN/dxC+dN/dxP>"<<"\t"
544	<<"<dN/dx>, 1/cm"<<G4endl<<G4endl ;
545
546
547	// G4PAIxSection testPAIproton(k,maxEnergyTransfer) ;
548
549	kineticEnergy = 10.0keV ; // 110MeV ;
550
551	// for(j=1;j<testPAIproton.GetNumberOfGammas();j++)
552
553	for(j=1;j<70;j++)
554	{
555	tau = kineticEnergy/proton_mass_c2 ;
556	gamma = tau +1.0 ;
557	bg2 = tau*(tau + 2.0) ;
558	beta2 = bg2/(gamma*gamma) ;
559	rateMass = electron_mass_c2/proton_mass_c2 ;
560
561	Tmax = 2.0electron_mass_c2bg2
562	/(1.0+2.0gammarateMass+rateMass*rateMass) ;
563
564
565	Tkin = maxEnergyTransfer ;
566
567	if ( maxEnergyTransfer > Tmax)
568	{
569	Tkin = Tmax ;
570	}
571	if ( Tmax <= Tmin + 0.5*eV )
572	{
573	Tkin = Tmin + 0.5*eV ;
574	}
575	G4PAIxSection testPAIproton(k,Tkin,bg2) ;
576
577	outFile
578	// << kineticEnergy/keV<<"\t"
579	<< gamma << "\t"
580	// << Tkin/keV<<"\t"
581	<< testPAIproton.GetIntegralCerenkov(1)*cm << "\t"
582	<< testPAIproton.GetIntegralPlasmon(1)*cm << "\t"
583	<< testPAIproton.GetIntegralCerenkov(1)*cm +
584	testPAIproton.GetIntegralPlasmon(1)*cm << "\t"
585	<< testPAIproton.GetIntegralPAIxSection(1)*cm << "\t" << G4endl ;
586	G4cout
587	// << kineticEnergy/keV<<"\t\t"
588	<< gamma << "\t\t"
589	// << Tkin/keV<<"\t\t"
590	<< testPAIproton.GetIntegralCerenkov(1)*cm << "\t"
591	<< testPAIproton.GetIntegralPlasmon(1)*cm << "\t"
592	<< testPAIproton.GetIntegralCerenkov(1)*cm +
593	testPAIproton.GetIntegralPlasmon(1)*cm << "\t"
594	<< testPAIproton.GetIntegralPAIxSection(1)*cm << "\t\t" << G4endl ;
595
596	// outFile<<testPAIproton.GetLorentzFactor(j)<<"\t"
597	// <<maxEnergyTransfer/keV<<"\t\t"
598	// <<testPAIproton.GetPAItable(0,j)*cm/keV<<"\t\t"
599	// <<testPAIproton.GetPAItable(1,j)*cm<<"\t\t"<<G4endl ;
600
601	kineticEnergy *= 1.41 ; // was 1.4 ; 1.5 ;
602	}
603	G4cout<<G4endl ;
604	outFile<<G4endl ;
605	}
606	return 1 ;
607
608	G4String confirm ;
609	G4cout<<"Enter 'y' , if you would like to get dE/dx-distribution : "
610	<<std::flush ;
611
612	G4cin>>confirm ;
613	if(confirm != "y" ) return 1 ;
614	G4cout<<G4endl ;
615
616	for(k=0;k<numOfMaterials;k++)
617	{
618	G4cout <<k<< " Material : " <<(*theMaterialTable)[k]->GetName() << G4endl ;
619	}
620	G4cout<<"Enter material name for dE/dx-distribution : "<<std::flush ;
621	G4cin>>testName ;
622	G4cout<<G4endl ;
623
624	G4int iLoss, iStat, iStatMax, nGamma ;
625	G4double energyLoss[50], Ebin, delta, delta1, delta2, delta3, step, y, pos ;
626	G4double intProb[200], colDist, sum, fact, GF, lambda, aaa ;
627
628	G4double alphaCrossTalk = -0.055, betaS = 0.20.4keV ;
629	G4int spectrum[50] ;
630
631	G4cout << " Enter nGamma 1<nGamma<10 : " <<std::flush ;
632	G4cin>>nGamma ;
633	G4cout<<G4endl ;
634
635	for(k=0;k<numOfMaterials;k++)
636	{
637	if((*theMaterialTable)[k]->GetName() != testName) continue ;
638
639	G4cout << "Material : " <<(*theMaterialTable)[k]->GetName() << G4endl<<G4endl ;
640
641
642	G4cout << " Enter Lorentz factor : " <<std::flush ;
643	G4cin>>gamma ;
644	G4cout<<G4endl ;
645
646	G4cout << " Enter step in mm : " <<std::flush ;
647	G4cin>>step ;
648	G4cout<<G4endl ;
649	step *= mm ;
650
651	G4cout << " Enter energy bin in keV : " <<std::flush ;
652	G4cin>>Ebin ;
653	G4cout<<G4endl ;
654	Ebin *= keV ;
655
656	G4cout << " Enter number of events : " <<std::flush ;
657	G4cin>>iStatMax ;
658
659	G4cout<<G4endl<<"Start dE/dx distribution"<<G4endl<<G4endl ;
660
661	maxEnergyTransfer = 100*keV ;
662	bg2 = gamma*gamma - 1 ;
663	rateMass = electron_mass_c2/proton_mass_c2 ;
664
665	Tmax = 2.0electron_mass_c2bg2
666	/(1.0+2.0gammarateMass+rateMass*rateMass) ;
667
668	if ( maxEnergyTransfer > Tmax) maxEnergyTransfer = Tmax ;
669
670	G4PAIxSection testPAIenergyLoss(k,maxEnergyTransfer,bg2) ;
671
672	for( iLoss = 0 ; iLoss < 50 ; iLoss++ )
673	{
674	energyLoss[iLoss] = Ebin*iLoss ;
675	spectrum[iLoss] = 0 ;
676	}
677	for(iStat=0;iStat<iStatMax;iStat++)
678	{
679
680	// aaa = (G4double)nGamma ;
681	// lambda = aaa/step ;
682	// colDist = RandGamma::shoot(aaa,lambda) ;
683
684	// delta = testPAIenergyLoss.GetStepEnergyLoss(colDist) ;
685
686	// delta = testPAIenergyLoss.GetStepEnergyLoss(step) ;
687
688	delta1 = testPAIenergyLoss.GetStepEnergyLoss(step) ;
689
690	delta = G4RandGauss::shoot(delta1,0.3*delta1) ;
691	if( delta < 0.0 ) delta = 0.0 ;
692
693	// delta2 = testPAIenergyLoss.GetStepEnergyLoss(step) ;
694	// delta3 = testPAIenergyLoss.GetStepEnergyLoss(step) ;
695
696	// delta = alphaCrossTalk*delta1 +
697	// delta2 + alphaCrossTalk*delta3 - betaS ;
698
699	for(iLoss=0;iLoss<50;iLoss++)
700	{
701	if(delta <= energyLoss[iLoss]) break ;
702	}
703	spectrum[iLoss-1]++ ;
704	}
705	G4double meanLoss = 0.0 ;
706
707	outFile<<"E, keV"<<"\t\t"<<"Distribution"<<G4endl<<G4endl ;
708	G4cout<<"E, keV"<<"\t\t"<<"Distribution"<<G4endl<<G4endl ;
709	G4cout<<G4endl ;
710	for(iLoss=0;iLoss<50;iLoss++) // with last bin
711	{
712	fileOut<<energyLoss[iLoss]/keV<<"\t\t"<<spectrum[iLoss]<<G4endl ;
713	G4cout<<energyLoss[iLoss]/keV<<"\t\t"<<spectrum[iLoss]<<G4endl ;
714	meanLoss +=energyLoss[iLoss]*spectrum[iLoss] ;
715	}
716	G4cout<<G4endl ;
717	G4cout<<"Mean loss over spectrum = "<<meanLoss/keV/iStatMax<<" keV"<<G4endl ;
718	}
719
720	G4int exit = 1 ;
721
722	while(exit)
723	{
724	G4cout<<"Enter 'y' , if you would like to compare with exp. data : "<<std::flush ;
725	G4cin>>confirm ;
726	if(confirm != "y" ) break ;
727	G4cout<<G4endl ;
728
729	// Read experimental data file
730
731	G4double delExp[200], distr[200], deltaBin, sumPAI, sumExp ;
732	G4int numberOfExpPoints ;
733
734	G4cout<<G4endl ;
735	G4cout << " Enter number of experimental points : " <<std::flush ;
736	G4cin>>numberOfExpPoints ;
737	G4cout<<G4endl ;
738	G4cout << " Enter energy bin in keV : " <<std::flush ;
739	G4cin>>deltaBin ;
740	G4cout<<G4endl ;
741	deltaBin *= keV ;
742
743	std::ifstream fileRead ;
744	fileRead.open("input.dat") ;
745	for(i=0;i<numberOfExpPoints;i++)
746	{
747	fileRead>>delExp[i]>>distr[i] ;
748	delExp[i] *= keV ;
749	G4cout<<i<<"\t"<<delExp[i]<<"\t"<<distr[i]<<G4endl ;
750	}
751	fileRead.close() ;
752
753	// Adjust statistics of experiment to PAI simulation
754
755	sumExp = 0.0 ;
756	for(i=0;i<numberOfExpPoints;i++) sumExp +=distr[i] ;
757	sumExp *= deltaBin ;
758
759	sumPAI = 0.0 ;
760	for(i=0;i<49;i++) sumPAI +=spectrum[i] ;
761	sumPAI *= Ebin ;
762
763	for(i=0;i<numberOfExpPoints;i++) distr[i] *= sumPAI/sumExp ;
764
765	for(i=0;i<numberOfExpPoints;i++)
766	{
767	fileWrite<<delExp[i]/keV<<"\t"<<distr[i]<<G4endl ;
768	G4cout<<delExp[i]/keV<<"\t"<<distr[i]<<G4endl ;
769	}
770	exit = 0 ;
771	}
772
773	G4cout<<"Enter 'y' , if you would like to get most probable delta : "<<std::flush ;
774	G4cin>>confirm ;
775	if(confirm != "y" ) return 1 ;
776	G4cout<<G4endl ;
777
778	G4int kGamma, iMPLoss, maxSpectrum, iMax ;
779	G4double mpDelta[50], meanDelta[50], rrMP[50], rrMean[50] ;
780	G4double mpLoss, tmRatio, mpSum, mpStat ;
781
782	G4double aGamma[33] =
783	{
784	4.0, 1.5, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 8.0, 10.0, // 13
785	20., 40.0, 60.0, 80.0, 100.0, 200.0, 400.0, 600.0, 800.0, 1000.0, // 23
786	2000.0, 4000.0, 6000.0, 8000.0, 100000.0, 20000.0, // 29
787	40000.0, 60000.0, 80000.0, 100000.0 // 33
788	} ;
789
790	for(k=0;k<numOfMaterials;k++)
791	{
792	G4cout <<k<< " Material : " <<(*theMaterialTable)[k]->GetName() << G4endl ;
793	}
794	G4cout<<"Enter material name for dE/dx-distribution : "<<std::flush ;
795	G4cin>>testName ;
796	G4cout<<G4endl ;
797
798
799	for(k=0;k<numOfMaterials;k++)
800	{
801	if((*theMaterialTable)[k]->GetName() != testName) continue ;
802
803	G4cout << "Material : " <<(*theMaterialTable)[k]->GetName() << G4endl<<G4endl ;
804
805	G4cout << " Enter nGamma 1<nGamma<10 : " <<std::flush ;
806	G4cin>>nGamma ;
807	G4cout<<G4endl ;
808
809
810	G4cout << " Enter step in mm : " <<std::flush ;
811	G4cin>>step ;
812	G4cout<<G4endl ;
813	step *= mm ;
814
815	G4cout << " Enter energy bin in keV : " <<std::flush ;
816	G4cin>>Ebin ;
817	G4cout<<G4endl ;
818	Ebin *= keV ;
819
820	G4cout << " Enter trancated mean ration <1.0 : " <<std::flush ;
821	G4cin>>tmRatio ;
822	G4cout<<G4endl ;
823
824
825	G4cout << " Enter number of events : " <<std::flush ;
826	G4cin>>iStatMax ;
827	G4cout<<G4endl ;
828
829	G4cout<<"no."<<"\t"<<"Gamma"<<"\t"<<"Rel. rise"<<"\t"<<"M.P. loss, keV"
830	<<"\t"<<"Mean loss, keV"<<G4endl<<G4endl ;
831	// outFile<<"no."<<"\t"<<"Gamma"<<"\t"<<"M.P. loss, keV"
832	// <<"\t"<<"Mean loss, keV"<<G4endl<<G4endl ;
833
834
835	// gamma = 1.1852 ;
836
837	for(kGamma=0;kGamma<33;kGamma++)
838	{
839	// G4cout<<G4endl<<"Start dE/dx distribution"<<G4endl<<G4endl ;
840
841	gamma = aGamma[kGamma] ;
842	maxEnergyTransfer = 100*keV ;
843	bg2 = gamma*gamma - 1 ;
844	rateMass = electron_mass_c2/proton_mass_c2 ;
845
846	Tmax = 2.0electron_mass_c2bg2
847	/(1.0+2.0gammarateMass+rateMass*rateMass) ;
848
849	if ( maxEnergyTransfer > Tmax) maxEnergyTransfer = Tmax ;
850
851	G4PAIxSection testPAIenergyLoss(k,maxEnergyTransfer,bg2) ;
852
853	for( iLoss = 0 ; iLoss < 50 ; iLoss++ )
854	{
855	energyLoss[iLoss] = Ebin*iLoss ;
856	spectrum[iLoss] = 0 ;
857	}
858	for(iStat=0;iStat<iStatMax;iStat++)
859	{
860
861	// aaa = (G4double)nGamma ;
862	// lambda = aaa/step ;
863	// colDist = RandGamma::shoot(aaa,lambda) ;
864
865	// delta = testPAIenergyLoss.GetStepEnergyLoss(colDist) ;
866
867	delta = testPAIenergyLoss.GetStepEnergyLoss(step) ;
868
869	// delta1 = testPAIenergyLoss.GetStepEnergyLoss(step) ;
870	// delta2 = testPAIenergyLoss.GetStepEnergyLoss(step) ;
871	// delta3 = testPAIenergyLoss.GetStepEnergyLoss(step) ;
872
873	// delta = alphaCrossTalk*delta1 +
874	// delta2 + alphaCrossTalk*delta3 - betaS ;
875
876	for(iLoss=0;iLoss<50;iLoss++)
877	{
878	if(delta <= energyLoss[iLoss]) break ;
879	}
880	spectrum[iLoss-1]++ ;
881	}
882	G4int sumStat = 0 ;
883	for(iLoss=0;iLoss<49;iLoss++) // without last bin
884	{
885	sumStat += spectrum[iLoss] ;
886	if( sumStat > tmRatio*iStatMax ) break ;
887	}
888	if(iLoss == 50) iLoss-- ;
889	iMPLoss = iLoss ;
890	G4double meanLoss = 0.0 ;
891	maxSpectrum = 0 ;
892
893	for(iLoss=0;iLoss<iMPLoss;iLoss++) // without last bin
894	{
895	// fileOut<<energyLoss[iLoss]/keV<<"\t\t"<<spectrum[iLoss]<<G4endl ;
896	// G4cout<<energyLoss[iLoss]/keV<<"\t\t"<<spectrum[iLoss]<<G4endl ;
897
898	meanLoss += energyLoss[iLoss]*spectrum[iLoss] ;
899
900	if( spectrum[iLoss] > maxSpectrum )
901	{
902	maxSpectrum = spectrum[iLoss] ;
903	mpLoss = energyLoss[iLoss] ;
904	iMax = iLoss ;
905	}
906	}
907	mpSum = 0. ;
908	mpStat = 0 ;
909	for(iLoss = iMax-5;iLoss<=iMax+5;iLoss++)
910	{
911	mpSum += energyLoss[iLoss]*spectrum[iLoss] ;
912	mpStat += spectrum[iLoss] ;
913	}
914	mpLoss = mpSum/mpStat ;
915	mpLoss /= keV ;
916	meanLoss /= keV*sumStat ;
917	meanDelta[kGamma] = meanLoss ;
918	mpDelta[kGamma] = mpLoss ;
919
920	if(kGamma > 0)
921	{
922	rrMP[kGamma] = mpLoss/mpDelta[0] ;
923	G4cout<<kGamma<<"\t"<<gamma<<"\t"<<rrMP[kGamma]<<"\t"<<mpLoss<<G4endl ;
924	// outFile<<gamma<<"\t"<<rrMP[kGamma]<<G4endl ;
925	fileWrite1<<gamma<<"\t"<<rrMP[kGamma]<<G4endl ;
926	}
927
928	// gamma *= 1.5 ;
929	}
930	G4cout<<G4endl ;
931	outFile<<G4endl ;
932	}
933
934	return EXIT_SUCCESS;
935
936	}
937
938
939
940
941
942
943
944
945
946

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

Context Navigation

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

Download in other formats: