Context Navigation

test90Ne10CO2pai.cc @ 1350

Last change on this file since 1350 was 1347, checked in by garnier, 14 years ago
geant4 tag 9.4
File size: 30.7 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.9 2010/11/23 15:31:10 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	// 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	// sandia.SetVerbose(1);
367	sanIndex = sandia.SandiaIntervals(thisMaterialZ,nbOfElements);
368	sanIndex = sandia.SandiaMixing( thisMaterialZ ,
369	(*theMaterialTable)[k]->GetFractionVector(),
370	nbOfElements,sanIndex);
371
372	for(row = 0; row < sanIndex-1; row++ )
373	{
374	G4cout<<row+1<<"\t"<<sandia.GetPhotoAbsorpCof(row+1,0)/keV;
375	// outFile<<row+1<<" "<<sandia.GetPhotoAbsorpCof(row+1,0)/keV;
376
377	for( iSan = 1; iSan < 5; iSan++ )
378	{
379	G4cout<<"\t"<<sandia.GetPhotoAbsorpCof(row+1,iSan);
380	// (theMaterialTable)[k]->GetDensity();
381
382	// outFile<<" "<<sandia.GetPhotoAbsorpCof(row+1,iSan);
383	// (theMaterialTable)[k]->GetDensity();
384	}
385	G4cout<<G4endl;
386	// outFile<<G4endl;
387	}
388	G4cout<<G4endl;
389	// outFile<<G4endl;
390
391
392	// outFile<<G4endl;
393	maxEnergyTransfer = 100*keV;
394	gamma = 4.0;
395	bg2 = gamma*gamma - 1;
396
397	G4PAIxSection testPAI( k, maxEnergyTransfer, bg2);
398
399	G4cout<<"Interval no."<<"\t"<<"Energy interval"<<G4endl<<G4endl;
400	// outFile<<"Interval no."<<"\t"<<"Energy interval"<<G4endl<<G4endl;
401
402	for( j = 1; j <= testPAI.GetIntervalNumber(); j++ )
403	{
404	G4cout<<j<<"\t\t"<<testPAI.GetEnergyInterval(j)/keV<<G4endl;
405	// outFile<<j<<"\t\t"<<testPAI.GetEnergyInterval(j)/keV<<G4endl;
406	}
407	G4cout<<G4endl;
408	// outFile<<G4endl;
409
410	G4cout << "Actual spline size = "<<testPAI.GetSplineSize()<<G4endl;
411	G4cout <<"Normalization Cof = "<<testPAI.GetNormalizationCof()<<G4endl;
412	G4cout << G4endl;
413
414	// outFile<<"Actual spline size = "<<testPAI.GetSplineSize()<<G4endl;
415	// outFile<<"Normalization Cof = "<<testPAI.GetNormalizationCof()<<G4endl;
416	// outFile<<G4endl;
417
418
419	Tmin = sandia.GetPhotoAbsorpCof(1,0); // 0.02*keV;
420	G4cout<<"Tmin = "<<Tmin/eV<<" eV"<<G4endl;
421
422	G4cout
423	// <<"Tkin, keV"<<"\t"
424	<< "bg"<<"\t\t"
425	// <<"Max E transfer, kev"<<"\t"
426	<< "<dN/dxC>, 1/cm"<<"\t"
427	<< "<dN/dxMM>, 1/cm"<<"\t"
428	<< "<dN/dxP>, 1/cm"<<"\t"
429	// << "<dN/dxC+dN/dxP>"<<"\t"
430	<<"<dN/dx>, 1/cm"<<G4endl<<G4endl;
431
432	/*
433	outFile
434	// <<"Tkin, keV"<<"\t"
435	<<"gamma"<<"\t\t"
436	// <<"Max E transfer, kev"<<"\t"
437	<<"<dN/dxC>, 1/cm"<<"\t"
438	<< "<dN/dxP>, 1/cm"<<"\t"
439	<<"<dN/dxC+dN/dxP>"<<"\t"
440	<<"<dN/dx>, 1/cm"<<G4endl<<G4endl;
441	*/
442	// G4PAIxSection testPAIproton(k,maxEnergyTransfer);
443
444
445
446
447
448	// kineticEnergy = 10.0keV; // 100.GeV; // 10.0keV; // 110MeV; // for proton
449
450	// kineticEnergy = 5*GeV; // for electrons
451
452	// kineticEnergy = 5GeVproton_mass_c2/electron_mass_c2;
453
454	kineticEnergy = 3.*GeV;
455
456	// for(j=1;j<testPAIproton.GetNumberOfGammas();j++)
457
458	// jMax = 70; // 70;
459	jMax = 1; // 70;
460
461	// outFile<<jMax<<G4endl;
462
463	for( j = 0; j < jMax; j++ )
464	{
465	tau = kineticEnergy/proton_mass_c2;
466	// tau = kineticEnergy/electron_mass_c2;
467	gamma = tau +1.0;
468	bg2 = tau*(tau + 2.0);
469	bg = std::sqrt(bg2);
470	beta2 = bg2/(gamma*gamma);
471	G4cout<<"bg = "<<bg<<"; b2 = "<<beta2<<G4endl<<G4endl;
472	rateMass = electron_mass_c2/proton_mass_c2;
473
474	Tmax = 2.0electron_mass_c2bg2/(1.0+2.0gammarateMass+rateMass*rateMass);
475	// Tmax = 0.5*kineticEnergy;
476
477	Tkin = maxEnergyTransfer;
478
479	if ( maxEnergyTransfer > Tmax)
480	{
481	Tkin = Tmax;
482	}
483	if ( Tmax <= Tmin + 0.5*eV )
484	{
485	Tkin = Tmin + 0.5*eV;
486	}
487	G4PAIxSection testPAIproton(k,Tkin,bg2);
488	/*
489	G4cout
490	// << kineticEnergy/keV<<"\t\t"
491	// << gamma << "\t\t"
492	<< bg << "\t\t"
493	// << Tkin/keV<<"\t\t"
494	<< testPAIproton.GetIntegralCerenkov(1)*cm << "\t"
495	<< testPAIproton.GetIntegralMM(1)*cm << "\t"
496	<< testPAIproton.GetIntegralPlasmon(1)*cm << "\t"
497	<< testPAIproton.GetIntegralResonance(1)*cm << "\t"
498	// << testPAIproton.GetIntegralCerenkov(1)*cm +
499	// testPAIproton.GetIntegralPlasmon(1)*cm << "\t"
500	// << FitALICE(bg) << "\t"
501	// << FitBichsel(bg) << "\t"
502	<< testPAIproton.GetIntegralPAIxSection(1)*cm << "\t\t"
503	<< G4endl;
504
505
506
507	outFile
508	// << kineticEnergy/keV<<"\t"
509	// << gamma << "\t"
510	<< bg << "\t\t"
511	// << Tkin/keV<<"\t"
512	<< testPAIproton.GetIntegralCerenkov(1)*cm << "\t"
513	<< testPAIproton.GetIntegralMM(1)*cm << "\t"
514	<< testPAIproton.GetIntegralPlasmon(1)*cm << "\t"
515	<< testPAIproton.GetIntegralResonance(1)*cm << "\t"
516	// << testPAIproton.GetIntegralCerenkov(1)*cm +
517	// testPAIproton.GetIntegralPlasmon(1)*cm << "\t"
518	// << FitALICE(bg) << "\t"
519	// << FitBichsel(bg) << "\t"
520	<< testPAIproton.GetIntegralPAIxSection(1)*cm << "\t"
521	<< G4endl;
522
523	// outFile<<testPAIproton.GetLorentzFactor(j)<<"\t"
524	// <<maxEnergyTransfer/keV<<"\t\t"
525	// <<testPAIproton.GetPAItable(0,j)*cm/keV<<"\t\t"
526	// <<testPAIproton.GetPAItable(1,j)*cm<<"\t\t"<<G4endl;
527
528	*/
529
530	outFile<<testPAIproton.GetSplineSize()-1<<G4endl;
531
532	// for( i = 1; i < testPAIproton.GetSplineSize(); i++)
533	for( i = 1; i < 10; i++)
534	{
535	outFile
536	<< testPAIproton.GetSplineEnergy(i)/keV << "\t"
537	// << testPAIproton.GetIntegralCerenkov(i)*cm << "\t"
538	// << testPAIproton.GetIntegralMM(i)*cm << "\t"
539	// << testPAIproton.GetIntegralPlasmon(i)*cm << "\t"
540	// << testPAIproton.GetIntegralResonance(i)*cm << "\t"
541	<< testPAIproton.GetIntegralPAIxSection(i)*cm << "\t"
542	<< G4endl;
543
544	G4cout
545	<< testPAIproton.GetSplineEnergy(i)/keV << "\t"
546	<< testPAIproton.GetIntegralCerenkov(i)*cm << "\t"
547	<< testPAIproton.GetIntegralMM(i)*cm << "\t"
548	<< testPAIproton.GetIntegralPlasmon(i)*cm << "\t"
549	<< testPAIproton.GetIntegralResonance(i)*cm << "\t"
550	<< testPAIproton.GetIntegralPAIxSection(i)*cm << "\t"
551	<< G4endl;
552
553	}
554
555
556
557
558	/*
559	G4double position, transfer, lambda, range, r2cer=0., r2res=0., r2ruth=0., r2tot=0.;
560	G4int nCer = 0, nRes = 0, nRuth = 0, nTot = 0;
561	G4double rBin[100], rDistr[100], rTemp, rTemp2, sumDistr = 0., rSum = 0;
562	G4double ionBin[100], ionDistr[100], ionMean, ionRand, F = 0.19, ionSum=0., ionSigma;
563
564
565	for( i = 0; i < 100; i++)
566	{
567	ionBin[i] = i*1.;
568	ionDistr[i] = 0.;
569	rBin[i] = i/200.;
570	rDistr[i] = 0.;
571	}
572	for( i = 0; i < 10000; i++)
573	{
574
575	position = testPAIproton.GetIntegralPAIxSection(1)*G4UniformRand();
576
577	if( position < testPAIproton.GetIntegralCerenkov(1) )
578	{
579	transfer = testPAIproton.GetCerenkovEnergyTransfer();
580	lambda = testPAIproton.GetPhotonRange(transfer);
581	range = testPAIproton.GetElectronRange(transfer);
582	r2cer += 0.67(lambda+range)(lambda+range);
583	r2tot += 0.67(lambda+range)(lambda+range);
584	rTemp2 = 0.67(lambda+range)(lambda+range);
585	nCer++;
586	}
587	else if( position < (testPAIproton.GetIntegralCerenkov(1)+
588	testPAIproton.GetIntegralResonance(1) ))
589	{
590	transfer = testPAIproton.GetResonanceEnergyTransfer();
591	range = testPAIproton.GetElectronRange(transfer);
592	r2res += 0.67rangerange;
593	r2tot += 0.67rangerange;
594	rTemp2 = 0.67rangerange;
595	nRes++;
596	}
597	else
598	{
599	transfer = testPAIproton.GetRutherfordEnergyTransfer();
600	range = testPAIproton.GetElectronRange(transfer);
601	r2ruth += range*range;
602	r2tot += range*range;
603	rTemp2 = range*range;
604	nRuth++;
605	}
606	nTot++;
607
608	rTemp = std::sqrt(rTemp2);
609	rSum += rTemp;
610
611	// rTemp = rTemp2;
612
613	for( j = 0; j < 100; j++ )
614	{
615	if( rTemp <= rBin[j] )
616	{
617	rDistr[j] += 1.;
618	break;
619	}
620	}
621
622
623	transfer = testPAIproton.GetEnergyTransfer();
624	ionMean = GetIonisation(transfer);
625	ionSigma = std::sqrt(F*ionMean);
626
627	// ionRand = G4RandGauss::shoot(ionMean, ionSigma);
628	ionRand = ionMean;
629
630	if( ionRand < 0.) ionRand =0.;
631
632	ionSum += ionRand;
633	nTot++;
634
635	for( j = 0; j < 100; j++ )
636	{
637	if( ionRand <= ionBin[j] )
638	{
639	ionDistr[j] += 1.;
640	break;
641	}
642	}
643	}
644
645	if(nCer >0) r2cer /= nCer;
646	if(nRes >0) r2res /= nRes;
647	if(nRuth >0) r2ruth /= nRuth;
648	r2tot /= nTot;
649	rSum /= nTot;
650	G4cout<<"nCer = "<<nCer<<"; nRes = "<<nRes<<"; nRuth = "<<nRuth<<G4endl;
651	G4cout<<"sum of n = "<<nCer+nRes+nRuth<<"; nTot = "<<nTot<<G4endl;
652	G4cout<<"rCer = "<<std::sqrt(r2cer)<<" mm; rRes = "<<std::sqrt(r2res)<<" mm"<<G4endl;
653	G4cout<<"rRuth = "<<std::sqrt(r2ruth)<<" mm; rTot = "<<std::sqrt(r2tot)<<" mm"<<G4endl;
654	G4cout<<"rSum = "<<rSum<<" mm; "<<G4endl;
655
656	ionSum /= nTot;
657	G4cout<<"ionSum = "<<ionSum<<" electrons"<<G4endl;
658
659	outFile<<100<<G4endl;
660
661	for( j = 0; j < 100; j++ )
662	{
663	// outFile<<rBin[j]<<"\t"<<rDistr[j]<<G4endl;
664	outFile<<ionBin[j]<<"\t"<<ionDistr[j]<<G4endl;
665	sumDistr += rDistr[j];
666	}
667	G4cout<<"sumDistr = "<<sumDistr<<G4endl;
668	*/
669
670
671
672	kineticEnergy *= 1.41; // was 1.4; 1.5;
673	}
674
675	G4cout<<G4endl;
676	// outFile<<G4endl;
677	}
678
679
680	return 1; // end of test
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696	G4String confirm;
697	G4cout<<"Enter 'y' , if you would like to get dE/dx-distribution : "
698	<<std::flush;
699
700	G4cin>>confirm;
701	if(confirm != "y" ) return 1;
702	G4cout<<G4endl;
703
704	for(k=0;k<numOfMaterials;k++)
705	{
706	G4cout <<k<< " Material : " <<(*theMaterialTable)[k]->GetName() << G4endl;
707	}
708	G4cout<<"Enter material name for dE/dx-distribution : "<<std::flush;
709	G4cin>>testName;
710	G4cout<<G4endl;
711
712	G4int iLoss, iStat, iStatMax, nGamma;
713	G4double energyLoss[50], Ebin, delta, delta1, delta2, delta3, step, y, pos;
714	G4double intProb[200], colDist, sum, fact, GF, lambda, aaa;
715
716	G4double alphaCrossTalk = -0.055, betaS = 0.20.4keV;
717	G4int spectrum[50];
718
719	G4cout << " Enter nGamma 1<nGamma<10 : " <<std::flush;
720	G4cin>>nGamma;
721	G4cout<<G4endl;
722
723	for(k=0;k<numOfMaterials;k++)
724	{
725	if((*theMaterialTable)[k]->GetName() != testName) continue;
726
727	G4cout << "Material : " <<(*theMaterialTable)[k]->GetName() << G4endl<<G4endl;
728
729
730	G4cout << " Enter Lorentz factor : " <<std::flush;
731	G4cin>>gamma;
732	G4cout<<G4endl;
733
734	G4cout << " Enter step in mm : " <<std::flush;
735	G4cin>>step;
736	G4cout<<G4endl;
737	step *= mm;
738
739	G4cout << " Enter energy bin in keV : " <<std::flush;
740	G4cin>>Ebin;
741	G4cout<<G4endl;
742	Ebin *= keV;
743
744	G4cout << " Enter number of events : " <<std::flush;
745	G4cin>>iStatMax;
746
747	G4cout<<G4endl<<"Start dE/dx distribution"<<G4endl<<G4endl;
748
749	maxEnergyTransfer = 100*keV;
750	bg2 = gamma*gamma - 1;
751	rateMass = electron_mass_c2/proton_mass_c2;
752
753	Tmax = 2.0electron_mass_c2bg2
754	/(1.0+2.0gammarateMass+rateMass*rateMass);
755
756	if ( maxEnergyTransfer > Tmax) maxEnergyTransfer = Tmax;
757
758	G4PAIxSection testPAIenergyLoss(k,maxEnergyTransfer,bg2);
759
760	for( iLoss = 0; iLoss < 50; iLoss++ )
761	{
762	energyLoss[iLoss] = Ebin*iLoss;
763	spectrum[iLoss] = 0;
764	}
765	for(iStat=0;iStat<iStatMax;iStat++)
766	{
767
768	// aaa = (G4double)nGamma;
769	// lambda = aaa/step;
770	// colDist = RandGamma::shoot(aaa,lambda);
771
772	// delta = testPAIenergyLoss.GetStepEnergyLoss(colDist);
773
774	// delta = testPAIenergyLoss.GetStepEnergyLoss(step);
775
776	delta1 = testPAIenergyLoss.GetStepEnergyLoss(step);
777
778	delta = G4RandGauss::shoot(delta1,0.3*delta1);
779	if( delta < 0.0 ) delta = 0.0;
780
781	// delta2 = testPAIenergyLoss.GetStepEnergyLoss(step);
782	// delta3 = testPAIenergyLoss.GetStepEnergyLoss(step);
783
784	// delta = alphaCrossTalk*delta1 +
785	// delta2 + alphaCrossTalk*delta3 - betaS;
786
787	for(iLoss=0;iLoss<50;iLoss++)
788	{
789	if(delta <= energyLoss[iLoss]) break;
790	}
791	spectrum[iLoss-1]++;
792	}
793	G4double meanLoss = 0.0;
794
795	outFile<<"E, keV"<<"\t\t"<<"Distribution"<<G4endl<<G4endl;
796	G4cout<<"E, keV"<<"\t\t"<<"Distribution"<<G4endl<<G4endl;
797	G4cout<<G4endl;
798	for(iLoss=0;iLoss<50;iLoss++) // with last bin
799	{
800	fileOut<<energyLoss[iLoss]/keV<<"\t\t"<<spectrum[iLoss]<<G4endl;
801	G4cout<<energyLoss[iLoss]/keV<<"\t\t"<<spectrum[iLoss]<<G4endl;
802	meanLoss +=energyLoss[iLoss]*spectrum[iLoss];
803	}
804	G4cout<<G4endl;
805	G4cout<<"Mean loss over spectrum = "<<meanLoss/keV/iStatMax<<" keV"<<G4endl;
806	}
807
808	G4int exit = 1;
809
810	while(exit)
811	{
812	G4cout<<"Enter 'y' , if you would like to compare with exp. data : "<<std::flush;
813	G4cin>>confirm;
814	if(confirm != "y" ) break;
815	G4cout<<G4endl;
816
817	// Read experimental data file
818
819	G4double delExp[200], distr[200], deltaBin, sumPAI, sumExp;
820	G4int numberOfExpPoints;
821
822	G4cout<<G4endl;
823	G4cout << " Enter number of experimental points : " <<std::flush;
824	G4cin>>numberOfExpPoints;
825	G4cout<<G4endl;
826	G4cout << " Enter energy bin in keV : " <<std::flush;
827	G4cin>>deltaBin;
828	G4cout<<G4endl;
829	deltaBin *= keV;
830
831	std::ifstream fileRead;
832	fileRead.open("input.dat");
833	for(i=0;i<numberOfExpPoints;i++)
834	{
835	fileRead>>delExp[i]>>distr[i];
836	delExp[i] *= keV;
837	G4cout<<i<<"\t"<<delExp[i]<<"\t"<<distr[i]<<G4endl;
838	}
839	fileRead.close();
840
841	// Adjust statistics of experiment to PAI simulation
842
843	sumExp = 0.0;
844	for(i=0;i<numberOfExpPoints;i++) sumExp +=distr[i];
845	sumExp *= deltaBin;
846
847	sumPAI = 0.0;
848	for(i=0;i<49;i++) sumPAI +=spectrum[i];
849	sumPAI *= Ebin;
850
851	for(i=0;i<numberOfExpPoints;i++) distr[i] *= sumPAI/sumExp;
852
853	for(i=0;i<numberOfExpPoints;i++)
854	{
855	fileWrite<<delExp[i]/keV<<"\t"<<distr[i]<<G4endl;
856	G4cout<<delExp[i]/keV<<"\t"<<distr[i]<<G4endl;
857	}
858	exit = 0;
859	}
860
861	G4cout<<"Enter 'y' , if you would like to get most probable delta : "<<std::flush;
862	G4cin>>confirm;
863	if(confirm != "y" ) return 1;
864	G4cout<<G4endl;
865
866	G4int kGamma, iMPLoss, maxSpectrum, iMax;
867	G4double mpDelta[50], meanDelta[50], rrMP[50], rrMean[50];
868	G4double mpLoss, tmRatio, mpSum, mpStat;
869
870	G4double aGamma[33] =
871	{
872	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
873	20., 40.0, 60.0, 80.0, 100.0, 200.0, 400.0, 600.0, 800.0, 1000.0, // 23
874	2000.0, 4000.0, 6000.0, 8000.0, 100000.0, 20000.0, // 29
875	40000.0, 60000.0, 80000.0, 100000.0 // 33
876	};
877
878	for(k=0;k<numOfMaterials;k++)
879	{
880	G4cout <<k<< " Material : " <<(*theMaterialTable)[k]->GetName() << G4endl;
881	}
882	G4cout<<"Enter material name for dE/dx-distribution : "<<std::flush;
883	G4cin>>testName;
884	G4cout<<G4endl;
885
886
887	for(k=0;k<numOfMaterials;k++)
888	{
889	if((*theMaterialTable)[k]->GetName() != testName) continue;
890
891	G4cout << "Material : " <<(*theMaterialTable)[k]->GetName() << G4endl<<G4endl;
892
893	G4cout << " Enter nGamma 1<nGamma<10 : " <<std::flush;
894	G4cin>>nGamma;
895	G4cout<<G4endl;
896
897
898	G4cout << " Enter step in mm : " <<std::flush;
899	G4cin>>step;
900	G4cout<<G4endl;
901	step *= mm;
902
903	G4cout << " Enter energy bin in keV : " <<std::flush;
904	G4cin>>Ebin;
905	G4cout<<G4endl;
906	Ebin *= keV;
907
908	G4cout << " Enter trancated mean ration <1.0 : " <<std::flush;
909	G4cin>>tmRatio;
910	G4cout<<G4endl;
911
912
913	G4cout << " Enter number of events : " <<std::flush;
914	G4cin>>iStatMax;
915	G4cout<<G4endl;
916
917	G4cout<<"no."<<"\t"<<"Gamma"<<"\t"<<"Rel. rise"<<"\t"<<"M.P. loss, keV"
918	<<"\t"<<"Mean loss, keV"<<G4endl<<G4endl;
919	// outFile<<"no."<<"\t"<<"Gamma"<<"\t"<<"M.P. loss, keV"
920	// <<"\t"<<"Mean loss, keV"<<G4endl<<G4endl;
921
922
923	// gamma = 1.1852;
924
925	for(kGamma=0;kGamma<33;kGamma++)
926	{
927	// G4cout<<G4endl<<"Start dE/dx distribution"<<G4endl<<G4endl;
928
929	gamma = aGamma[kGamma];
930	maxEnergyTransfer = 100*keV;
931	bg2 = gamma*gamma - 1;
932	rateMass = electron_mass_c2/proton_mass_c2;
933
934	Tmax = 2.0electron_mass_c2bg2
935	/(1.0+2.0gammarateMass+rateMass*rateMass);
936
937	if ( maxEnergyTransfer > Tmax) maxEnergyTransfer = Tmax;
938
939	G4PAIxSection testPAIenergyLoss(k,maxEnergyTransfer,bg2);
940
941	for( iLoss = 0; iLoss < 50; iLoss++ )
942	{
943	energyLoss[iLoss] = Ebin*iLoss;
944	spectrum[iLoss] = 0;
945	}
946	for(iStat=0;iStat<iStatMax;iStat++)
947	{
948
949	// aaa = (G4double)nGamma;
950	// lambda = aaa/step;
951	// colDist = RandGamma::shoot(aaa,lambda);
952
953	// delta = testPAIenergyLoss.GetStepEnergyLoss(colDist);
954
955	delta = testPAIenergyLoss.GetStepEnergyLoss(step);
956
957	// delta1 = testPAIenergyLoss.GetStepEnergyLoss(step);
958	// delta2 = testPAIenergyLoss.GetStepEnergyLoss(step);
959	// delta3 = testPAIenergyLoss.GetStepEnergyLoss(step);
960
961	// delta = alphaCrossTalk*delta1 +
962	// delta2 + alphaCrossTalk*delta3 - betaS;
963
964	for(iLoss=0;iLoss<50;iLoss++)
965	{
966	if(delta <= energyLoss[iLoss]) break;
967	}
968	spectrum[iLoss-1]++;
969	}
970	G4int sumStat = 0;
971	for(iLoss=0;iLoss<49;iLoss++) // without last bin
972	{
973	sumStat += spectrum[iLoss];
974	if( sumStat > tmRatio*iStatMax ) break;
975	}
976	if(iLoss == 50) iLoss--;
977	iMPLoss = iLoss;
978	G4double meanLoss = 0.0;
979	maxSpectrum = 0;
980
981	for(iLoss=0;iLoss<iMPLoss;iLoss++) // without last bin
982	{
983	// fileOut<<energyLoss[iLoss]/keV<<"\t\t"<<spectrum[iLoss]<<G4endl;
984	// G4cout<<energyLoss[iLoss]/keV<<"\t\t"<<spectrum[iLoss]<<G4endl;
985
986	meanLoss += energyLoss[iLoss]*spectrum[iLoss];
987
988	if( spectrum[iLoss] > maxSpectrum )
989	{
990	maxSpectrum = spectrum[iLoss] ;
991	mpLoss = energyLoss[iLoss];
992	iMax = iLoss;
993	}
994	}
995	mpSum = 0.;
996	mpStat = 0;
997	for(iLoss = iMax-5;iLoss<=iMax+5;iLoss++)
998	{
999	mpSum += energyLoss[iLoss]*spectrum[iLoss];
1000	mpStat += spectrum[iLoss];
1001	}
1002	mpLoss = mpSum/mpStat;
1003	mpLoss /= keV;
1004	meanLoss /= keV*sumStat;
1005	meanDelta[kGamma] = meanLoss;
1006	mpDelta[kGamma] = mpLoss;
1007
1008	if(kGamma > 0)
1009	{
1010	rrMP[kGamma] = mpLoss/mpDelta[0];
1011	G4cout<<kGamma<<"\t"<<gamma<<"\t"<<rrMP[kGamma]<<"\t"<<mpLoss<<G4endl;
1012	// outFile<<gamma<<"\t"<<rrMP[kGamma]<<G4endl;
1013	fileWrite1<<gamma<<"\t"<<rrMP[kGamma]<<G4endl;
1014	}
1015
1016	// gamma *= 1.5;
1017	}
1018	G4cout<<G4endl;
1019	outFile<<G4endl;
1020	}
1021
1022	return EXIT_SUCCESS;
1023
1024	}
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034

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

Context Navigation

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

Download in other formats: