Context Navigation

source: trunk/source/processes/electromagnetic/standard/test/G4InitXscPAItest.cc @ 1199

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

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

Download in other formats: