Context Navigation

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

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

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

Download in other formats: