source: trunk/source/processes/electromagnetic/standard/test/G4PAIdNdxTest.cc@ 1330

Last change on this file since 1330 was 1315, checked in by garnier, 15 years ago

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