Context Navigation

source: trunk/source/processes/electromagnetic/standard/test/G4PAIxSectionTest.cc @ 1315

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago
update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 28.4 KB

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

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

Download in other formats: