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