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