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

Last change on this file since 1228 was 1199, checked in by garnier, 16 years ago

nvx fichiers dans CVS
File size: 28.4 KB
RevLine 
[1199]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.15 2006/06/29 19:54:15 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// 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
58int 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.15*kelvin,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.15*kelvin,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.15*kelvin,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.15*kelvin,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.15*kelvin,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 ;
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 = 10.0*keV ; // 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.0*electron_mass_c2*bg2
558 /(1.0+2.0*gamma*rateMass+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.2*0.4*keV ;
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.0*electron_mass_c2*bg2
656 /(1.0+2.0*gamma*rateMass+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.0*electron_mass_c2*bg2
837 /(1.0+2.0*gamma*rateMass+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.