Context Navigation

source: trunk/source/processes/electromagnetic/standard/test/test90Ne10CO2pai.cc @ 1326

Last change on this file since 1326 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: 30.6 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: test90Ne10CO2pai.cc,v 1.8 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	// 11.04.08, V. Grichine test of PAI predictions for 90% Ne + 10% CO2
39	// ALICE TPC gas mixture
40	// 8.12.09 V. Grichine update for T2K gas mixture
41
42	#include "G4ios.hh"
43	#include <fstream>
44	#include <cmath>
45
46
47	#include "globals.hh"
48	#include "Randomize.hh"
49
50	#include "G4Isotope.hh"
51	#include "G4Element.hh"
52	#include "G4Material.hh"
53	#include "G4MaterialTable.hh"
54	#include "G4SandiaTable.hh"
55
56	// #include "G4PAIonisation.hh"
57	#include "G4PAIxSection.hh"
58
59	// Peter:
60	// From the AliRoot code in $ALICE_ROOT/TPC/AliTPCv2.cxx:
61	// const Float_t kprim = 14.35; // number of primary collisions per 1 cm
62	// const Float_t kpoti = 20.77e-9; // first ionization potential for Ne/CO2
63	// const Float_t kwIon = 35.97e-9; // energy for the ion-electron pair creation
64	//
65	// kprim is the number of primary collisions per cm for a MIP!
66	// kpoti = I.
67	// kwIon = W.
68
69	G4double FitALICE(G4double bg)
70	{
71	//
72	// Bethe-Bloch energy loss formula from ALICE TPC TRD
73	//
74	const G4double kp1 = 0.76176e-1;
75	const G4double kp2 = 10.632;
76	const G4double kp3 = 0.13279e-4;
77	const G4double kp4 = 1.8631;
78	const G4double kp5 = 1.9479;
79	const G4double dn1 = 14.35;
80
81	G4double dbg = (G4double) bg;
82
83	G4double beta = dbg/std::sqrt(1.+dbg*dbg);
84
85	G4double aa = std::pow(beta,kp4);
86	G4double bb = std::pow(1./dbg,kp5);
87
88
89	bb = std::log(kp3 + bb);
90
91	G4double result = ( kp2 - aa - bb)*kp1/aa;
92
93	result *= dn1;
94
95	return result;
96	}
97
98
99	G4double FitBichsel(G4double bg)
100	{
101	//
102	// Primary ionisation from Hans Bichsel fit
103	//
104	const G4double kp1 = 0.686e-1;
105	const G4double kp2 = 11.714;
106	const G4double kp3 = 0.218e-4;
107	const G4double kp4 = 1.997;
108	const G4double kp5 = 2.133;
109	const G4double dn1 = 13.32;
110
111	G4double dbg = (G4double) bg;
112
113	G4double beta = dbg/std::sqrt(1.+dbg*dbg);
114
115	G4double aa = std::pow(beta,kp4);
116	G4double bb = std::pow(1./dbg,kp5);
117
118
119	bb=std::log(kp3 + bb);
120
121	G4double result = ( kp2 - aa - bb)*kp1/aa;
122
123	result *= dn1;
124
125	return result;
126
127	}
128
129	////////////////////////////////////////////////////////////
130
131
132	G4double GetIonisation(G4double transfer)
133	{
134	G4double W = 34.75*eV;
135	G4double I1 = 13.62*eV; // first ionisation potential in mixture
136	I1 *= 0.9;
137
138	G4double result = W;
139
140	// result /= 1.-I1/transfer;
141
142	return transfer/result;
143
144	}
145
146
147	/////////////////////////////////////////////////
148
149
150
151	int main()
152	{
153	// std::ofstream outFile("90Ne10CO2pai.dat", std::ios::out );
154	std::ofstream outFile("e5GeVt2kPhilippe.dat", std::ios::out );
155	outFile.setf( std::ios::scientific, std::ios::floatfield );
156
157	std::ofstream fileOut("PAICerPlasm90Ne10CO2.dat", std::ios::out );
158	fileOut.setf( std::ios::scientific, std::ios::floatfield );
159
160	// std::ifstream fileRead("exp.dat", std::ios::out );
161	// fileRead.setf( std::ios::scientific, std::ios::floatfield );
162
163	std::ofstream fileWrite("exp.dat", std::ios::out );
164	fileWrite.setf( std::ios::scientific, std::ios::floatfield );
165
166	std::ofstream fileWrite1("mprrpai.dat", std::ios::out );
167	fileWrite1.setf( std::ios::scientific, std::ios::floatfield );
168
169	// Create materials
170
171
172	G4int iz , n, nel, ncomponents;
173	G4double a, z, ez, density , temperature, pressure, fractionmass;
174	G4State state;
175	G4String name, symbol;
176
177
178	a = 1.01*g/mole;
179	G4Isotope* ih1 = new G4Isotope("Hydrogen",iz=1,n=1,a);
180
181	a = 2.01*g/mole;
182	G4Isotope* ih2 = new G4Isotope("Deuterium",iz=1,n=2,a);
183
184	G4Element* elH = new G4Element(name="Hydrogen",symbol="H",2);
185	elH->AddIsotope(ih1,.999);
186	elH->AddIsotope(ih2,.001);
187
188	a = 12.01*g/mole;
189	G4Element* elC = new G4Element(name="Carbon",symbol="C", ez=6., a);
190
191	a = 14.01*g/mole;
192	G4Element* elN = new G4Element(name="Nitrogen", symbol="N", ez=7., a);
193
194	a = 16.00*g/mole;
195	G4Element* elO = new G4Element(name="Oxygen",symbol="O", ez=8., a);
196
197
198	a = 19.00*g/mole;
199	G4Element* elF = new G4Element(name="Fluorine", symbol="F", z=9., a);
200
201	a = 39.948*g/mole;
202	G4Element* elAr = new G4Element(name="Argon", symbol="Ar", z=18., a);
203
204	// Neon as detector gas, STP
205
206	density = 0.900*mg/cm3;
207	a = 20.179*g/mole;
208	G4Material* Ne = new G4Material(name="Ne",z=10., a, density );
209
210	// Carbone dioxide, CO2 STP
211
212	density = 1.977*mg/cm3;
213	G4Material* CarbonDioxide = new G4Material(name="CO2", density, nel=2);
214	CarbonDioxide->AddElement(elC,1);
215	CarbonDioxide->AddElement(elO,2);
216
217
218
219	// 90% Ne + 10% CO2, STP
220
221	density = 1.0077*mg/cm3;
222	G4Material* Ne10CO2 = new G4Material(name="Ne10CO2" , density,
223	ncomponents=2);
224	Ne10CO2->AddMaterial( Ne, fractionmass = 0.8038 );
225	Ne10CO2->AddMaterial( CarbonDioxide, fractionmass = 0.1962 );
226
227	density *= 273./293.;
228
229	G4Material* Ne10CO2T293 = new G4Material(name="Ne10CO2T293" , density,
230	ncomponents=2);
231	Ne10CO2T293->AddMaterial( Ne, fractionmass = 0.8038 );
232	Ne10CO2T293->AddMaterial( CarbonDioxide, fractionmass = 0.1962 );
233
234
235	density = 1.25053*mg/cm3; // STP
236	G4Material* Nitrogen = new G4Material(name="N2" , density, ncomponents=1);
237	Nitrogen->AddElement(elN, 2);
238
239	// 85.7% Ne + 9.5% CO2 +4.8% N2, STP
240
241	density = 1.0191*mg/cm3;
242	G4Material* Ne857CO295N2 = new G4Material(name="Ne857CO295N2" , density,
243	ncomponents=3);
244	Ne857CO295N2->AddMaterial( Ne, fractionmass = 0.7568 );
245	Ne857CO295N2->AddMaterial( CarbonDioxide, fractionmass = 0.1843 );
246	Ne857CO295N2->AddMaterial( Nitrogen, fractionmass = 0.0589 );
247
248	density *= 273./292.;
249	density *= 0.966/1.01325;
250
251	// G4cout<<"density of Ne857CO295N2T292 = "<<density*cm3/mg<<" mg/cm3"<<G4endl;
252
253	G4Material* Ne857CO295N2T292 = new G4Material(name="Ne857CO295N2T292" , density,
254	ncomponents=3);
255	Ne857CO295N2T292->AddMaterial( Ne, fractionmass = 0.76065 );
256	Ne857CO295N2T292->AddMaterial( CarbonDioxide, fractionmass = 0.18140 );
257	Ne857CO295N2T292->AddMaterial( Nitrogen, fractionmass = 0.05795 );
258
259
260
261	// Ar as detector gas,STP
262
263	density = 1.7836*mg/cm3 ; // STP
264	G4Material* Argon = new G4Material(name="Argon" , density, ncomponents=1);
265	Argon->AddElement(elAr, 1);
266	/*
267	// iso-Butane (methylpropane), STP
268
269	density = 2.67*mg/cm3 ;
270	G4Material* isobutane = new G4Material(name="isoC4H10",density,nel=2) ;
271	isobutane->AddElement(elC,4) ;
272	isobutane->AddElement(elH,10) ;
273
274	// CF4 from ATLAS TRT estimation
275
276	G4double TRT_CF4_density = 3.9*mg/cm3;
277	G4Material* TRT_CF4 = new G4Material(name="TRT_CF4", TRT_CF4_density, nel=2,
278	kStateGas,293.15kelvin,1.atmosphere);
279	*/
280
281	// Philippe Gros T2K mixture version
282	// Argon
283	/*
284	density = 1.66*mg/cm3;
285	pressure = 1*atmosphere;
286	temperature = 288.15*kelvin;
287	G4Material* Argon = new G4Material(name="Ar", // z=18., a=39.948*g/mole,
288	density, ncomponents=1); // kStateGas,temperature,pressure);
289	Argon->AddElement(elAr, 1);
290	*/
291
292	// IsoButane
293
294	density = 2.51*mg/cm3;
295	G4Material* Isobu = new G4Material(name="isoC4H10", z=34.,a=58.123*g/mole,
296	density, kStateGas,temperature,pressure);
297
298	// Tetrafluoromethane
299
300	density = 3.72*mg/cm3;
301	G4Material* FlMet = new
302	G4Material(name="CF4",z=42.,a=88.01*g/mole,density,kStateGas,temperature,pressure);
303
304	// Argon + 3% tetrafluoromethane + 2% iso-butane
305
306	density = 1.748*mg/cm3;
307	G4Material* t2kGasMixture = new G4Material(name="t2kGasMixture", density, ncomponents=3);
308	t2kGasMixture->AddMaterial(Argon, fractionmass = 90.9*perCent);
309	t2kGasMixture->AddMaterial(FlMet, fractionmass = 6.3*perCent);
310	t2kGasMixture->AddMaterial(Isobu, fractionmass = 2.8*perCent);
311
312
313	G4int i, j, jMax, k, numOfMaterials, iSan, nbOfElements, sanIndex, row;
314	G4double maxEnergyTransfer, kineticEnergy;
315	G4double tau, gamma, bg2, bg, beta2, rateMass, Tmax, Tmin, Tkin;
316
317	const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
318
319	numOfMaterials = theMaterialTable->size();
320
321	G4cout<<"Available materials under test : "<< G4endl<<G4endl;
322	// outFile<<"Available materials under test : "<< G4endl<<G4endl;
323
324	for( k = 0; k < numOfMaterials; k++ )
325	{
326	G4cout <<k<<"\t"<< " Material : " <<(*theMaterialTable)[k]->GetName() << G4endl;
327	// outFile <<k<<"\t"<< " Material : " <<(*theMaterialTable)[k]->GetName() << G4endl;
328	}
329
330
331
332	// G4String testName = "N2";
333	// G4String testName = "Ne10CO2";
334	G4String testName = "t2kGasMixture";
335	// G4String testName = "Ar";
336	// G4String testName = "Argon";
337	// G4String testName = "Ne10CO2T293";
338	// G4String testName = "Ne857CO295N2T292";
339
340
341	// G4cout<<"Enter material name for test : "<<std::flush;
342	// G4cin>>testName;
343
344
345
346	for( k = 0; k < numOfMaterials; k++ )
347	{
348	if((*theMaterialTable)[k]->GetName() != testName) continue;
349
350	// outFile << "Material : " <<(*theMaterialTable)[k]->GetName() << G4endl;
351	G4cout << "Material : " <<(*theMaterialTable)[k]->GetName() << G4endl;
352
353	nbOfElements = (*theMaterialTable)[k]->GetNumberOfElements();
354
355	G4cout<<"Sandia cof according old PAI stuff"<<G4endl<<G4endl;
356	// outFile<<"Sandia cof according old PAI stuff"<<G4endl<<G4endl;
357
358	G4int* thisMaterialZ = new G4int[nbOfElements];
359
360	for( iSan = 0; iSan < nbOfElements; iSan++ )
361	{
362	thisMaterialZ[iSan] = (G4int)(*theMaterialTable)[k]->
363	GetElement(iSan)->GetZ();
364	}
365	G4SandiaTable sandia(k);
366	sanIndex = sandia.SandiaIntervals(thisMaterialZ,nbOfElements);
367	sanIndex = sandia.SandiaMixing( thisMaterialZ ,
368	(*theMaterialTable)[k]->GetFractionVector(),
369	nbOfElements,sanIndex);
370
371	for(row = 0; row < sanIndex-1; row++ )
372	{
373	G4cout<<row+1<<"\t"<<sandia.GetPhotoAbsorpCof(row+1,0)/keV;
374	// outFile<<row+1<<" "<<sandia.GetPhotoAbsorpCof(row+1,0)/keV;
375
376	for( iSan = 1; iSan < 5; iSan++ )
377	{
378	G4cout<<"\t"<<sandia.GetPhotoAbsorpCof(row+1,iSan);
379	// (theMaterialTable)[k]->GetDensity();
380
381	// outFile<<" "<<sandia.GetPhotoAbsorpCof(row+1,iSan);
382	// (theMaterialTable)[k]->GetDensity();
383	}
384	G4cout<<G4endl;
385	// outFile<<G4endl;
386	}
387	G4cout<<G4endl;
388	// outFile<<G4endl;
389
390
391	// outFile<<G4endl;
392	maxEnergyTransfer = 100*keV;
393	gamma = 4.0;
394	bg2 = gamma*gamma - 1;
395
396	G4PAIxSection testPAI( k, maxEnergyTransfer, bg2);
397
398	G4cout<<"Interval no."<<"\t"<<"Energy interval"<<G4endl<<G4endl;
399	// outFile<<"Interval no."<<"\t"<<"Energy interval"<<G4endl<<G4endl;
400
401	for( j = 1; j <= testPAI.GetIntervalNumber(); j++ )
402	{
403	G4cout<<j<<"\t\t"<<testPAI.GetEnergyInterval(j)/keV<<G4endl;
404	// outFile<<j<<"\t\t"<<testPAI.GetEnergyInterval(j)/keV<<G4endl;
405	}
406	G4cout<<G4endl;
407	// outFile<<G4endl;
408
409	G4cout << "Actual spline size = "<<testPAI.GetSplineSize()<<G4endl;
410	G4cout <<"Normalization Cof = "<<testPAI.GetNormalizationCof()<<G4endl;
411	G4cout << G4endl;
412
413	// outFile<<"Actual spline size = "<<testPAI.GetSplineSize()<<G4endl;
414	// outFile<<"Normalization Cof = "<<testPAI.GetNormalizationCof()<<G4endl;
415	// outFile<<G4endl;
416
417
418	Tmin = sandia.GetPhotoAbsorpCof(1,0); // 0.02*keV;
419	G4cout<<"Tmin = "<<Tmin/eV<<" eV"<<G4endl;
420
421	G4cout
422	// <<"Tkin, keV"<<"\t"
423	<< "bg"<<"\t\t"
424	// <<"Max E transfer, kev"<<"\t"
425	<< "<dN/dxC>, 1/cm"<<"\t"
426	<< "<dN/dxMM>, 1/cm"<<"\t"
427	<< "<dN/dxP>, 1/cm"<<"\t"
428	// << "<dN/dxC+dN/dxP>"<<"\t"
429	<<"<dN/dx>, 1/cm"<<G4endl<<G4endl;
430
431	/*
432	outFile
433	// <<"Tkin, keV"<<"\t"
434	<<"gamma"<<"\t\t"
435	// <<"Max E transfer, kev"<<"\t"
436	<<"<dN/dxC>, 1/cm"<<"\t"
437	<< "<dN/dxP>, 1/cm"<<"\t"
438	<<"<dN/dxC+dN/dxP>"<<"\t"
439	<<"<dN/dx>, 1/cm"<<G4endl<<G4endl;
440	*/
441	// G4PAIxSection testPAIproton(k,maxEnergyTransfer);
442
443
444
445
446
447	// kineticEnergy = 10.0keV; // 100.GeV; // 10.0keV; // 110MeV; // for proton
448
449	// kineticEnergy = 5*GeV; // for electrons
450
451	kineticEnergy = 5GeVproton_mass_c2/electron_mass_c2;
452
453	// kineticEnergy = 5*GeV;
454
455	// for(j=1;j<testPAIproton.GetNumberOfGammas();j++)
456
457	// jMax = 70; // 70;
458	jMax = 1; // 70;
459
460	// outFile<<jMax<<G4endl;
461
462	for( j = 0; j < jMax; j++ )
463	{
464	tau = kineticEnergy/proton_mass_c2;
465	// tau = kineticEnergy/electron_mass_c2;
466	gamma = tau +1.0;
467	bg2 = tau*(tau + 2.0);
468	bg = std::sqrt(bg2);
469	beta2 = bg2/(gamma*gamma);
470	G4cout<<"bg = "<<bg<<"; b2 = "<<beta2<<G4endl<<G4endl;
471	rateMass = electron_mass_c2/proton_mass_c2;
472
473	Tmax = 2.0electron_mass_c2bg2/(1.0+2.0gammarateMass+rateMass*rateMass);
474	// Tmax = 0.5*kineticEnergy;
475
476	Tkin = maxEnergyTransfer;
477
478	if ( maxEnergyTransfer > Tmax)
479	{
480	Tkin = Tmax;
481	}
482	if ( Tmax <= Tmin + 0.5*eV )
483	{
484	Tkin = Tmin + 0.5*eV;
485	}
486	G4PAIxSection testPAIproton(k,Tkin,bg2);
487	/*
488	G4cout
489	// << kineticEnergy/keV<<"\t\t"
490	// << gamma << "\t\t"
491	<< bg << "\t\t"
492	// << Tkin/keV<<"\t\t"
493	<< testPAIproton.GetIntegralCerenkov(1)*cm << "\t"
494	<< testPAIproton.GetIntegralMM(1)*cm << "\t"
495	<< testPAIproton.GetIntegralPlasmon(1)*cm << "\t"
496	<< testPAIproton.GetIntegralResonance(1)*cm << "\t"
497	// << testPAIproton.GetIntegralCerenkov(1)*cm +
498	// testPAIproton.GetIntegralPlasmon(1)*cm << "\t"
499	// << FitALICE(bg) << "\t"
500	// << FitBichsel(bg) << "\t"
501	<< testPAIproton.GetIntegralPAIxSection(1)*cm << "\t\t"
502	<< G4endl;
503
504
505
506	outFile
507	// << kineticEnergy/keV<<"\t"
508	// << gamma << "\t"
509	<< bg << "\t\t"
510	// << Tkin/keV<<"\t"
511	<< testPAIproton.GetIntegralCerenkov(1)*cm << "\t"
512	<< testPAIproton.GetIntegralMM(1)*cm << "\t"
513	<< testPAIproton.GetIntegralPlasmon(1)*cm << "\t"
514	<< testPAIproton.GetIntegralResonance(1)*cm << "\t"
515	// << testPAIproton.GetIntegralCerenkov(1)*cm +
516	// testPAIproton.GetIntegralPlasmon(1)*cm << "\t"
517	// << FitALICE(bg) << "\t"
518	// << FitBichsel(bg) << "\t"
519	<< testPAIproton.GetIntegralPAIxSection(1)*cm << "\t"
520	<< G4endl;
521
522	// outFile<<testPAIproton.GetLorentzFactor(j)<<"\t"
523	// <<maxEnergyTransfer/keV<<"\t\t"
524	// <<testPAIproton.GetPAItable(0,j)*cm/keV<<"\t\t"
525	// <<testPAIproton.GetPAItable(1,j)*cm<<"\t\t"<<G4endl;
526
527	*/
528
529	outFile<<testPAIproton.GetSplineSize()-1<<G4endl;
530
531	for( i = 1; i < testPAIproton.GetSplineSize(); i++)
532	{
533	outFile
534	<< testPAIproton.GetSplineEnergy(i)/keV << "\t"
535	// << testPAIproton.GetIntegralCerenkov(i)*cm << "\t"
536	// << testPAIproton.GetIntegralMM(i)*cm << "\t"
537	// << testPAIproton.GetIntegralPlasmon(i)*cm << "\t"
538	// << testPAIproton.GetIntegralResonance(i)*cm << "\t"
539	<< testPAIproton.GetIntegralPAIxSection(i)*cm << "\t"
540	<< G4endl;
541
542	G4cout
543	<< testPAIproton.GetSplineEnergy(i)/keV << "\t"
544	<< testPAIproton.GetIntegralCerenkov(i)*cm << "\t"
545	<< testPAIproton.GetIntegralMM(i)*cm << "\t"
546	<< testPAIproton.GetIntegralPlasmon(i)*cm << "\t"
547	<< testPAIproton.GetIntegralResonance(i)*cm << "\t"
548	<< testPAIproton.GetIntegralPAIxSection(i)*cm << "\t"
549	<< G4endl;
550
551	}
552
553
554
555
556	/*
557	G4double position, transfer, lambda, range, r2cer=0., r2res=0., r2ruth=0., r2tot=0.;
558	G4int nCer = 0, nRes = 0, nRuth = 0, nTot = 0;
559	G4double rBin[100], rDistr[100], rTemp, rTemp2, sumDistr = 0., rSum = 0;
560	G4double ionBin[100], ionDistr[100], ionMean, ionRand, F = 0.19, ionSum=0., ionSigma;
561
562
563	for( i = 0; i < 100; i++)
564	{
565	ionBin[i] = i*1.;
566	ionDistr[i] = 0.;
567	rBin[i] = i/200.;
568	rDistr[i] = 0.;
569	}
570	for( i = 0; i < 10000; i++)
571	{
572
573	position = testPAIproton.GetIntegralPAIxSection(1)*G4UniformRand();
574
575	if( position < testPAIproton.GetIntegralCerenkov(1) )
576	{
577	transfer = testPAIproton.GetCerenkovEnergyTransfer();
578	lambda = testPAIproton.GetPhotonRange(transfer);
579	range = testPAIproton.GetElectronRange(transfer);
580	r2cer += 0.67(lambda+range)(lambda+range);
581	r2tot += 0.67(lambda+range)(lambda+range);
582	rTemp2 = 0.67(lambda+range)(lambda+range);
583	nCer++;
584	}
585	else if( position < (testPAIproton.GetIntegralCerenkov(1)+
586	testPAIproton.GetIntegralResonance(1) ))
587	{
588	transfer = testPAIproton.GetResonanceEnergyTransfer();
589	range = testPAIproton.GetElectronRange(transfer);
590	r2res += 0.67rangerange;
591	r2tot += 0.67rangerange;
592	rTemp2 = 0.67rangerange;
593	nRes++;
594	}
595	else
596	{
597	transfer = testPAIproton.GetRutherfordEnergyTransfer();
598	range = testPAIproton.GetElectronRange(transfer);
599	r2ruth += range*range;
600	r2tot += range*range;
601	rTemp2 = range*range;
602	nRuth++;
603	}
604	nTot++;
605
606	rTemp = std::sqrt(rTemp2);
607	rSum += rTemp;
608
609	// rTemp = rTemp2;
610
611	for( j = 0; j < 100; j++ )
612	{
613	if( rTemp <= rBin[j] )
614	{
615	rDistr[j] += 1.;
616	break;
617	}
618	}
619
620
621	transfer = testPAIproton.GetEnergyTransfer();
622	ionMean = GetIonisation(transfer);
623	ionSigma = std::sqrt(F*ionMean);
624
625	// ionRand = G4RandGauss::shoot(ionMean, ionSigma);
626	ionRand = ionMean;
627
628	if( ionRand < 0.) ionRand =0.;
629
630	ionSum += ionRand;
631	nTot++;
632
633	for( j = 0; j < 100; j++ )
634	{
635	if( ionRand <= ionBin[j] )
636	{
637	ionDistr[j] += 1.;
638	break;
639	}
640	}
641	}
642
643	if(nCer >0) r2cer /= nCer;
644	if(nRes >0) r2res /= nRes;
645	if(nRuth >0) r2ruth /= nRuth;
646	r2tot /= nTot;
647	rSum /= nTot;
648	G4cout<<"nCer = "<<nCer<<"; nRes = "<<nRes<<"; nRuth = "<<nRuth<<G4endl;
649	G4cout<<"sum of n = "<<nCer+nRes+nRuth<<"; nTot = "<<nTot<<G4endl;
650	G4cout<<"rCer = "<<std::sqrt(r2cer)<<" mm; rRes = "<<std::sqrt(r2res)<<" mm"<<G4endl;
651	G4cout<<"rRuth = "<<std::sqrt(r2ruth)<<" mm; rTot = "<<std::sqrt(r2tot)<<" mm"<<G4endl;
652	G4cout<<"rSum = "<<rSum<<" mm; "<<G4endl;
653
654	ionSum /= nTot;
655	G4cout<<"ionSum = "<<ionSum<<" electrons"<<G4endl;
656
657	outFile<<100<<G4endl;
658
659	for( j = 0; j < 100; j++ )
660	{
661	// outFile<<rBin[j]<<"\t"<<rDistr[j]<<G4endl;
662	outFile<<ionBin[j]<<"\t"<<ionDistr[j]<<G4endl;
663	sumDistr += rDistr[j];
664	}
665	G4cout<<"sumDistr = "<<sumDistr<<G4endl;
666	*/
667
668
669
670	kineticEnergy *= 1.41; // was 1.4; 1.5;
671	}
672
673	G4cout<<G4endl;
674	// outFile<<G4endl;
675	}
676
677
678	return 1; // end of test
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694	G4String confirm;
695	G4cout<<"Enter 'y' , if you would like to get dE/dx-distribution : "
696	<<std::flush;
697
698	G4cin>>confirm;
699	if(confirm != "y" ) return 1;
700	G4cout<<G4endl;
701
702	for(k=0;k<numOfMaterials;k++)
703	{
704	G4cout <<k<< " Material : " <<(*theMaterialTable)[k]->GetName() << G4endl;
705	}
706	G4cout<<"Enter material name for dE/dx-distribution : "<<std::flush;
707	G4cin>>testName;
708	G4cout<<G4endl;
709
710	G4int iLoss, iStat, iStatMax, nGamma;
711	G4double energyLoss[50], Ebin, delta, delta1, delta2, delta3, step, y, pos;
712	G4double intProb[200], colDist, sum, fact, GF, lambda, aaa;
713
714	G4double alphaCrossTalk = -0.055, betaS = 0.20.4keV;
715	G4int spectrum[50];
716
717	G4cout << " Enter nGamma 1<nGamma<10 : " <<std::flush;
718	G4cin>>nGamma;
719	G4cout<<G4endl;
720
721	for(k=0;k<numOfMaterials;k++)
722	{
723	if((*theMaterialTable)[k]->GetName() != testName) continue;
724
725	G4cout << "Material : " <<(*theMaterialTable)[k]->GetName() << G4endl<<G4endl;
726
727
728	G4cout << " Enter Lorentz factor : " <<std::flush;
729	G4cin>>gamma;
730	G4cout<<G4endl;
731
732	G4cout << " Enter step in mm : " <<std::flush;
733	G4cin>>step;
734	G4cout<<G4endl;
735	step *= mm;
736
737	G4cout << " Enter energy bin in keV : " <<std::flush;
738	G4cin>>Ebin;
739	G4cout<<G4endl;
740	Ebin *= keV;
741
742	G4cout << " Enter number of events : " <<std::flush;
743	G4cin>>iStatMax;
744
745	G4cout<<G4endl<<"Start dE/dx distribution"<<G4endl<<G4endl;
746
747	maxEnergyTransfer = 100*keV;
748	bg2 = gamma*gamma - 1;
749	rateMass = electron_mass_c2/proton_mass_c2;
750
751	Tmax = 2.0electron_mass_c2bg2
752	/(1.0+2.0gammarateMass+rateMass*rateMass);
753
754	if ( maxEnergyTransfer > Tmax) maxEnergyTransfer = Tmax;
755
756	G4PAIxSection testPAIenergyLoss(k,maxEnergyTransfer,bg2);
757
758	for( iLoss = 0; iLoss < 50; iLoss++ )
759	{
760	energyLoss[iLoss] = Ebin*iLoss;
761	spectrum[iLoss] = 0;
762	}
763	for(iStat=0;iStat<iStatMax;iStat++)
764	{
765
766	// aaa = (G4double)nGamma;
767	// lambda = aaa/step;
768	// colDist = RandGamma::shoot(aaa,lambda);
769
770	// delta = testPAIenergyLoss.GetStepEnergyLoss(colDist);
771
772	// delta = testPAIenergyLoss.GetStepEnergyLoss(step);
773
774	delta1 = testPAIenergyLoss.GetStepEnergyLoss(step);
775
776	delta = G4RandGauss::shoot(delta1,0.3*delta1);
777	if( delta < 0.0 ) delta = 0.0;
778
779	// delta2 = testPAIenergyLoss.GetStepEnergyLoss(step);
780	// delta3 = testPAIenergyLoss.GetStepEnergyLoss(step);
781
782	// delta = alphaCrossTalk*delta1 +
783	// delta2 + alphaCrossTalk*delta3 - betaS;
784
785	for(iLoss=0;iLoss<50;iLoss++)
786	{
787	if(delta <= energyLoss[iLoss]) break;
788	}
789	spectrum[iLoss-1]++;
790	}
791	G4double meanLoss = 0.0;
792
793	outFile<<"E, keV"<<"\t\t"<<"Distribution"<<G4endl<<G4endl;
794	G4cout<<"E, keV"<<"\t\t"<<"Distribution"<<G4endl<<G4endl;
795	G4cout<<G4endl;
796	for(iLoss=0;iLoss<50;iLoss++) // with last bin
797	{
798	fileOut<<energyLoss[iLoss]/keV<<"\t\t"<<spectrum[iLoss]<<G4endl;
799	G4cout<<energyLoss[iLoss]/keV<<"\t\t"<<spectrum[iLoss]<<G4endl;
800	meanLoss +=energyLoss[iLoss]*spectrum[iLoss];
801	}
802	G4cout<<G4endl;
803	G4cout<<"Mean loss over spectrum = "<<meanLoss/keV/iStatMax<<" keV"<<G4endl;
804	}
805
806	G4int exit = 1;
807
808	while(exit)
809	{
810	G4cout<<"Enter 'y' , if you would like to compare with exp. data : "<<std::flush;
811	G4cin>>confirm;
812	if(confirm != "y" ) break;
813	G4cout<<G4endl;
814
815	// Read experimental data file
816
817	G4double delExp[200], distr[200], deltaBin, sumPAI, sumExp;
818	G4int numberOfExpPoints;
819
820	G4cout<<G4endl;
821	G4cout << " Enter number of experimental points : " <<std::flush;
822	G4cin>>numberOfExpPoints;
823	G4cout<<G4endl;
824	G4cout << " Enter energy bin in keV : " <<std::flush;
825	G4cin>>deltaBin;
826	G4cout<<G4endl;
827	deltaBin *= keV;
828
829	std::ifstream fileRead;
830	fileRead.open("input.dat");
831	for(i=0;i<numberOfExpPoints;i++)
832	{
833	fileRead>>delExp[i]>>distr[i];
834	delExp[i] *= keV;
835	G4cout<<i<<"\t"<<delExp[i]<<"\t"<<distr[i]<<G4endl;
836	}
837	fileRead.close();
838
839	// Adjust statistics of experiment to PAI simulation
840
841	sumExp = 0.0;
842	for(i=0;i<numberOfExpPoints;i++) sumExp +=distr[i];
843	sumExp *= deltaBin;
844
845	sumPAI = 0.0;
846	for(i=0;i<49;i++) sumPAI +=spectrum[i];
847	sumPAI *= Ebin;
848
849	for(i=0;i<numberOfExpPoints;i++) distr[i] *= sumPAI/sumExp;
850
851	for(i=0;i<numberOfExpPoints;i++)
852	{
853	fileWrite<<delExp[i]/keV<<"\t"<<distr[i]<<G4endl;
854	G4cout<<delExp[i]/keV<<"\t"<<distr[i]<<G4endl;
855	}
856	exit = 0;
857	}
858
859	G4cout<<"Enter 'y' , if you would like to get most probable delta : "<<std::flush;
860	G4cin>>confirm;
861	if(confirm != "y" ) return 1;
862	G4cout<<G4endl;
863
864	G4int kGamma, iMPLoss, maxSpectrum, iMax;
865	G4double mpDelta[50], meanDelta[50], rrMP[50], rrMean[50];
866	G4double mpLoss, tmRatio, mpSum, mpStat;
867
868	G4double aGamma[33] =
869	{
870	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
871	20., 40.0, 60.0, 80.0, 100.0, 200.0, 400.0, 600.0, 800.0, 1000.0, // 23
872	2000.0, 4000.0, 6000.0, 8000.0, 100000.0, 20000.0, // 29
873	40000.0, 60000.0, 80000.0, 100000.0 // 33
874	};
875
876	for(k=0;k<numOfMaterials;k++)
877	{
878	G4cout <<k<< " Material : " <<(*theMaterialTable)[k]->GetName() << G4endl;
879	}
880	G4cout<<"Enter material name for dE/dx-distribution : "<<std::flush;
881	G4cin>>testName;
882	G4cout<<G4endl;
883
884
885	for(k=0;k<numOfMaterials;k++)
886	{
887	if((*theMaterialTable)[k]->GetName() != testName) continue;
888
889	G4cout << "Material : " <<(*theMaterialTable)[k]->GetName() << G4endl<<G4endl;
890
891	G4cout << " Enter nGamma 1<nGamma<10 : " <<std::flush;
892	G4cin>>nGamma;
893	G4cout<<G4endl;
894
895
896	G4cout << " Enter step in mm : " <<std::flush;
897	G4cin>>step;
898	G4cout<<G4endl;
899	step *= mm;
900
901	G4cout << " Enter energy bin in keV : " <<std::flush;
902	G4cin>>Ebin;
903	G4cout<<G4endl;
904	Ebin *= keV;
905
906	G4cout << " Enter trancated mean ration <1.0 : " <<std::flush;
907	G4cin>>tmRatio;
908	G4cout<<G4endl;
909
910
911	G4cout << " Enter number of events : " <<std::flush;
912	G4cin>>iStatMax;
913	G4cout<<G4endl;
914
915	G4cout<<"no."<<"\t"<<"Gamma"<<"\t"<<"Rel. rise"<<"\t"<<"M.P. loss, keV"
916	<<"\t"<<"Mean loss, keV"<<G4endl<<G4endl;
917	// outFile<<"no."<<"\t"<<"Gamma"<<"\t"<<"M.P. loss, keV"
918	// <<"\t"<<"Mean loss, keV"<<G4endl<<G4endl;
919
920
921	// gamma = 1.1852;
922
923	for(kGamma=0;kGamma<33;kGamma++)
924	{
925	// G4cout<<G4endl<<"Start dE/dx distribution"<<G4endl<<G4endl;
926
927	gamma = aGamma[kGamma];
928	maxEnergyTransfer = 100*keV;
929	bg2 = gamma*gamma - 1;
930	rateMass = electron_mass_c2/proton_mass_c2;
931
932	Tmax = 2.0electron_mass_c2bg2
933	/(1.0+2.0gammarateMass+rateMass*rateMass);
934
935	if ( maxEnergyTransfer > Tmax) maxEnergyTransfer = Tmax;
936
937	G4PAIxSection testPAIenergyLoss(k,maxEnergyTransfer,bg2);
938
939	for( iLoss = 0; iLoss < 50; iLoss++ )
940	{
941	energyLoss[iLoss] = Ebin*iLoss;
942	spectrum[iLoss] = 0;
943	}
944	for(iStat=0;iStat<iStatMax;iStat++)
945	{
946
947	// aaa = (G4double)nGamma;
948	// lambda = aaa/step;
949	// colDist = RandGamma::shoot(aaa,lambda);
950
951	// delta = testPAIenergyLoss.GetStepEnergyLoss(colDist);
952
953	delta = testPAIenergyLoss.GetStepEnergyLoss(step);
954
955	// delta1 = testPAIenergyLoss.GetStepEnergyLoss(step);
956	// delta2 = testPAIenergyLoss.GetStepEnergyLoss(step);
957	// delta3 = testPAIenergyLoss.GetStepEnergyLoss(step);
958
959	// delta = alphaCrossTalk*delta1 +
960	// delta2 + alphaCrossTalk*delta3 - betaS;
961
962	for(iLoss=0;iLoss<50;iLoss++)
963	{
964	if(delta <= energyLoss[iLoss]) break;
965	}
966	spectrum[iLoss-1]++;
967	}
968	G4int sumStat = 0;
969	for(iLoss=0;iLoss<49;iLoss++) // without last bin
970	{
971	sumStat += spectrum[iLoss];
972	if( sumStat > tmRatio*iStatMax ) break;
973	}
974	if(iLoss == 50) iLoss--;
975	iMPLoss = iLoss;
976	G4double meanLoss = 0.0;
977	maxSpectrum = 0;
978
979	for(iLoss=0;iLoss<iMPLoss;iLoss++) // without last bin
980	{
981	// fileOut<<energyLoss[iLoss]/keV<<"\t\t"<<spectrum[iLoss]<<G4endl;
982	// G4cout<<energyLoss[iLoss]/keV<<"\t\t"<<spectrum[iLoss]<<G4endl;
983
984	meanLoss += energyLoss[iLoss]*spectrum[iLoss];
985
986	if( spectrum[iLoss] > maxSpectrum )
987	{
988	maxSpectrum = spectrum[iLoss] ;
989	mpLoss = energyLoss[iLoss];
990	iMax = iLoss;
991	}
992	}
993	mpSum = 0.;
994	mpStat = 0;
995	for(iLoss = iMax-5;iLoss<=iMax+5;iLoss++)
996	{
997	mpSum += energyLoss[iLoss]*spectrum[iLoss];
998	mpStat += spectrum[iLoss];
999	}
1000	mpLoss = mpSum/mpStat;
1001	mpLoss /= keV;
1002	meanLoss /= keV*sumStat;
1003	meanDelta[kGamma] = meanLoss;
1004	mpDelta[kGamma] = mpLoss;
1005
1006	if(kGamma > 0)
1007	{
1008	rrMP[kGamma] = mpLoss/mpDelta[0];
1009	G4cout<<kGamma<<"\t"<<gamma<<"\t"<<rrMP[kGamma]<<"\t"<<mpLoss<<G4endl;
1010	// outFile<<gamma<<"\t"<<rrMP[kGamma]<<G4endl;
1011	fileWrite1<<gamma<<"\t"<<rrMP[kGamma]<<G4endl;
1012	}
1013
1014	// gamma *= 1.5;
1015	}
1016	G4cout<<G4endl;
1017	outFile<<G4endl;
1018	}
1019
1020	return EXIT_SUCCESS;
1021
1022	}
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032

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

Download in other formats: