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

Last change on this file since 1351 was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 30.7 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: test90Ne10CO2pai.cc,v 1.9 2010/11/23 15:31:10 grichine Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30//
31//
32//
33//
34// Test routine for G4PAIxSection class code
35//
36// History:
37//
38// 11.04.08, V. Grichine test of PAI predictions for 90% Ne + 10% CO2
39// ALICE TPC gas mixture
40// 8.12.09 V. Grichine update for T2K gas mixture
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{
153 // std::ofstream outFile("90Ne10CO2pai.dat", std::ios::out );
154 std::ofstream outFile("e5GeVt2kPhilippe.dat", std::ios::out );
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
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
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
197
198 a = 19.00*g/mole;
199 G4Element* elF = new G4Element(name="Fluorine", symbol="F", z=9., a);
200
201 a = 39.948*g/mole;
202 G4Element* elAr = new G4Element(name="Argon", symbol="Ar", z=18., a);
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
251 // G4cout<<"density of Ne857CO295N2T292 = "<<density*cm3/mg<<" mg/cm3"<<G4endl;
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
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
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";
333 // G4String testName = "Ne10CO2";
334 // G4String testName = "t2kGasMixture";
335 // G4String testName = "Ar";
336 G4String testName = "Argon";
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 // sandia.SetVerbose(1);
367 sanIndex = sandia.SandiaIntervals(thisMaterialZ,nbOfElements);
368 sanIndex = sandia.SandiaMixing( thisMaterialZ ,
369 (*theMaterialTable)[k]->GetFractionVector(),
370 nbOfElements,sanIndex);
371
372 for(row = 0; row < sanIndex-1; row++ )
373 {
374 G4cout<<row+1<<"\t"<<sandia.GetPhotoAbsorpCof(row+1,0)/keV;
375 // outFile<<row+1<<" "<<sandia.GetPhotoAbsorpCof(row+1,0)/keV;
376
377 for( iSan = 1; iSan < 5; iSan++ )
378 {
379 G4cout<<"\t"<<sandia.GetPhotoAbsorpCof(row+1,iSan);
380 // *(*theMaterialTable)[k]->GetDensity();
381
382 // outFile<<" "<<sandia.GetPhotoAbsorpCof(row+1,iSan);
383 // *(*theMaterialTable)[k]->GetDensity();
384 }
385 G4cout<<G4endl;
386 // outFile<<G4endl;
387 }
388 G4cout<<G4endl;
389 // outFile<<G4endl;
390
391
392 // outFile<<G4endl;
393 maxEnergyTransfer = 100*keV;
394 gamma = 4.0;
395 bg2 = gamma*gamma - 1;
396
397 G4PAIxSection testPAI( k, maxEnergyTransfer, bg2);
398
399 G4cout<<"Interval no."<<"\t"<<"Energy interval"<<G4endl<<G4endl;
400 // outFile<<"Interval no."<<"\t"<<"Energy interval"<<G4endl<<G4endl;
401
402 for( j = 1; j <= testPAI.GetIntervalNumber(); j++ )
403 {
404 G4cout<<j<<"\t\t"<<testPAI.GetEnergyInterval(j)/keV<<G4endl;
405 // outFile<<j<<"\t\t"<<testPAI.GetEnergyInterval(j)/keV<<G4endl;
406 }
407 G4cout<<G4endl;
408 // outFile<<G4endl;
409
410 G4cout << "Actual spline size = "<<testPAI.GetSplineSize()<<G4endl;
411 G4cout <<"Normalization Cof = "<<testPAI.GetNormalizationCof()<<G4endl;
412 G4cout << G4endl;
413
414 // outFile<<"Actual spline size = "<<testPAI.GetSplineSize()<<G4endl;
415 // outFile<<"Normalization Cof = "<<testPAI.GetNormalizationCof()<<G4endl;
416 // outFile<<G4endl;
417
418
419 Tmin = sandia.GetPhotoAbsorpCof(1,0); // 0.02*keV;
420 G4cout<<"Tmin = "<<Tmin/eV<<" eV"<<G4endl;
421
422 G4cout
423 // <<"Tkin, keV"<<"\t"
424 << "bg"<<"\t\t"
425 // <<"Max E transfer, kev"<<"\t"
426 << "<dN/dxC>, 1/cm"<<"\t"
427 << "<dN/dxMM>, 1/cm"<<"\t"
428 << "<dN/dxP>, 1/cm"<<"\t"
429 // << "<dN/dxC+dN/dxP>"<<"\t"
430 <<"<dN/dx>, 1/cm"<<G4endl<<G4endl;
431
432 /*
433 outFile
434 // <<"Tkin, keV"<<"\t"
435 <<"gamma"<<"\t\t"
436 // <<"Max E transfer, kev"<<"\t"
437 <<"<dN/dxC>, 1/cm"<<"\t"
438 << "<dN/dxP>, 1/cm"<<"\t"
439 <<"<dN/dxC+dN/dxP>"<<"\t"
440 <<"<dN/dx>, 1/cm"<<G4endl<<G4endl;
441 */
442 // G4PAIxSection testPAIproton(k,maxEnergyTransfer);
443
444
445
446
447
448 // kineticEnergy = 10.0*keV; // 100.*GeV; // 10.0*keV; // 110*MeV; // for proton
449
450 // kineticEnergy = 5*GeV; // for electrons
451
452 // kineticEnergy = 5*GeV*proton_mass_c2/electron_mass_c2;
453
454 kineticEnergy = 3.*GeV;
455
456 // for(j=1;j<testPAIproton.GetNumberOfGammas();j++)
457
458 // jMax = 70; // 70;
459 jMax = 1; // 70;
460
461 // outFile<<jMax<<G4endl;
462
463 for( j = 0; j < jMax; j++ )
464 {
465 tau = kineticEnergy/proton_mass_c2;
466 // tau = kineticEnergy/electron_mass_c2;
467 gamma = tau +1.0;
468 bg2 = tau*(tau + 2.0);
469 bg = std::sqrt(bg2);
470 beta2 = bg2/(gamma*gamma);
471 G4cout<<"bg = "<<bg<<"; b2 = "<<beta2<<G4endl<<G4endl;
472 rateMass = electron_mass_c2/proton_mass_c2;
473
474 Tmax = 2.0*electron_mass_c2*bg2/(1.0+2.0*gamma*rateMass+rateMass*rateMass);
475 // Tmax = 0.5*kineticEnergy;
476
477 Tkin = maxEnergyTransfer;
478
479 if ( maxEnergyTransfer > Tmax)
480 {
481 Tkin = Tmax;
482 }
483 if ( Tmax <= Tmin + 0.5*eV )
484 {
485 Tkin = Tmin + 0.5*eV;
486 }
487 G4PAIxSection testPAIproton(k,Tkin,bg2);
488 /*
489 G4cout
490 // << kineticEnergy/keV<<"\t\t"
491 // << gamma << "\t\t"
492 << bg << "\t\t"
493 // << Tkin/keV<<"\t\t"
494 << testPAIproton.GetIntegralCerenkov(1)*cm << "\t"
495 << testPAIproton.GetIntegralMM(1)*cm << "\t"
496 << testPAIproton.GetIntegralPlasmon(1)*cm << "\t"
497 << testPAIproton.GetIntegralResonance(1)*cm << "\t"
498 // << testPAIproton.GetIntegralCerenkov(1)*cm +
499 // testPAIproton.GetIntegralPlasmon(1)*cm << "\t"
500 // << FitALICE(bg) << "\t"
501 // << FitBichsel(bg) << "\t"
502 << testPAIproton.GetIntegralPAIxSection(1)*cm << "\t\t"
503 << G4endl;
504
505
506
507 outFile
508 // << kineticEnergy/keV<<"\t"
509 // << gamma << "\t"
510 << bg << "\t\t"
511 // << Tkin/keV<<"\t"
512 << testPAIproton.GetIntegralCerenkov(1)*cm << "\t"
513 << testPAIproton.GetIntegralMM(1)*cm << "\t"
514 << testPAIproton.GetIntegralPlasmon(1)*cm << "\t"
515 << testPAIproton.GetIntegralResonance(1)*cm << "\t"
516 // << testPAIproton.GetIntegralCerenkov(1)*cm +
517 // testPAIproton.GetIntegralPlasmon(1)*cm << "\t"
518 // << FitALICE(bg) << "\t"
519 // << FitBichsel(bg) << "\t"
520 << testPAIproton.GetIntegralPAIxSection(1)*cm << "\t"
521 << G4endl;
522
523 // outFile<<testPAIproton.GetLorentzFactor(j)<<"\t"
524 // <<maxEnergyTransfer/keV<<"\t\t"
525 // <<testPAIproton.GetPAItable(0,j)*cm/keV<<"\t\t"
526 // <<testPAIproton.GetPAItable(1,j)*cm<<"\t\t"<<G4endl;
527
528 */
529
530 outFile<<testPAIproton.GetSplineSize()-1<<G4endl;
531
532 // for( i = 1; i < testPAIproton.GetSplineSize(); i++)
533 for( i = 1; i < 10; i++)
534 {
535 outFile
536 << testPAIproton.GetSplineEnergy(i)/keV << "\t"
537 // << testPAIproton.GetIntegralCerenkov(i)*cm << "\t"
538 // << testPAIproton.GetIntegralMM(i)*cm << "\t"
539 // << testPAIproton.GetIntegralPlasmon(i)*cm << "\t"
540 // << testPAIproton.GetIntegralResonance(i)*cm << "\t"
541 << testPAIproton.GetIntegralPAIxSection(i)*cm << "\t"
542 << G4endl;
543
544 G4cout
545 << testPAIproton.GetSplineEnergy(i)/keV << "\t"
546 << testPAIproton.GetIntegralCerenkov(i)*cm << "\t"
547 << testPAIproton.GetIntegralMM(i)*cm << "\t"
548 << testPAIproton.GetIntegralPlasmon(i)*cm << "\t"
549 << testPAIproton.GetIntegralResonance(i)*cm << "\t"
550 << testPAIproton.GetIntegralPAIxSection(i)*cm << "\t"
551 << G4endl;
552
553 }
554
555
556
557
558 /*
559 G4double position, transfer, lambda, range, r2cer=0., r2res=0., r2ruth=0., r2tot=0.;
560 G4int nCer = 0, nRes = 0, nRuth = 0, nTot = 0;
561 G4double rBin[100], rDistr[100], rTemp, rTemp2, sumDistr = 0., rSum = 0;
562 G4double ionBin[100], ionDistr[100], ionMean, ionRand, F = 0.19, ionSum=0., ionSigma;
563
564
565 for( i = 0; i < 100; i++)
566 {
567 ionBin[i] = i*1.;
568 ionDistr[i] = 0.;
569 rBin[i] = i/200.;
570 rDistr[i] = 0.;
571 }
572 for( i = 0; i < 10000; i++)
573 {
574
575 position = testPAIproton.GetIntegralPAIxSection(1)*G4UniformRand();
576
577 if( position < testPAIproton.GetIntegralCerenkov(1) )
578 {
579 transfer = testPAIproton.GetCerenkovEnergyTransfer();
580 lambda = testPAIproton.GetPhotonRange(transfer);
581 range = testPAIproton.GetElectronRange(transfer);
582 r2cer += 0.67*(lambda+range)*(lambda+range);
583 r2tot += 0.67*(lambda+range)*(lambda+range);
584 rTemp2 = 0.67*(lambda+range)*(lambda+range);
585 nCer++;
586 }
587 else if( position < (testPAIproton.GetIntegralCerenkov(1)+
588 testPAIproton.GetIntegralResonance(1) ))
589 {
590 transfer = testPAIproton.GetResonanceEnergyTransfer();
591 range = testPAIproton.GetElectronRange(transfer);
592 r2res += 0.67*range*range;
593 r2tot += 0.67*range*range;
594 rTemp2 = 0.67*range*range;
595 nRes++;
596 }
597 else
598 {
599 transfer = testPAIproton.GetRutherfordEnergyTransfer();
600 range = testPAIproton.GetElectronRange(transfer);
601 r2ruth += range*range;
602 r2tot += range*range;
603 rTemp2 = range*range;
604 nRuth++;
605 }
606 nTot++;
607
608 rTemp = std::sqrt(rTemp2);
609 rSum += rTemp;
610
611 // rTemp = rTemp2;
612
613 for( j = 0; j < 100; j++ )
614 {
615 if( rTemp <= rBin[j] )
616 {
617 rDistr[j] += 1.;
618 break;
619 }
620 }
621
622
623 transfer = testPAIproton.GetEnergyTransfer();
624 ionMean = GetIonisation(transfer);
625 ionSigma = std::sqrt(F*ionMean);
626
627 // ionRand = G4RandGauss::shoot(ionMean, ionSigma);
628 ionRand = ionMean;
629
630 if( ionRand < 0.) ionRand =0.;
631
632 ionSum += ionRand;
633 nTot++;
634
635 for( j = 0; j < 100; j++ )
636 {
637 if( ionRand <= ionBin[j] )
638 {
639 ionDistr[j] += 1.;
640 break;
641 }
642 }
643 }
644
645 if(nCer >0) r2cer /= nCer;
646 if(nRes >0) r2res /= nRes;
647 if(nRuth >0) r2ruth /= nRuth;
648 r2tot /= nTot;
649 rSum /= nTot;
650 G4cout<<"nCer = "<<nCer<<"; nRes = "<<nRes<<"; nRuth = "<<nRuth<<G4endl;
651 G4cout<<"sum of n = "<<nCer+nRes+nRuth<<"; nTot = "<<nTot<<G4endl;
652 G4cout<<"rCer = "<<std::sqrt(r2cer)<<" mm; rRes = "<<std::sqrt(r2res)<<" mm"<<G4endl;
653 G4cout<<"rRuth = "<<std::sqrt(r2ruth)<<" mm; rTot = "<<std::sqrt(r2tot)<<" mm"<<G4endl;
654 G4cout<<"rSum = "<<rSum<<" mm; "<<G4endl;
655
656 ionSum /= nTot;
657 G4cout<<"ionSum = "<<ionSum<<" electrons"<<G4endl;
658
659 outFile<<100<<G4endl;
660
661 for( j = 0; j < 100; j++ )
662 {
663 // outFile<<rBin[j]<<"\t"<<rDistr[j]<<G4endl;
664 outFile<<ionBin[j]<<"\t"<<ionDistr[j]<<G4endl;
665 sumDistr += rDistr[j];
666 }
667 G4cout<<"sumDistr = "<<sumDistr<<G4endl;
668 */
669
670
671
672 kineticEnergy *= 1.41; // was 1.4; 1.5;
673 }
674
675 G4cout<<G4endl;
676 // outFile<<G4endl;
677 }
678
679
680 return 1; // end of test
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696 G4String confirm;
697 G4cout<<"Enter 'y' , if you would like to get dE/dx-distribution : "
698 <<std::flush;
699
700 G4cin>>confirm;
701 if(confirm != "y" ) return 1;
702 G4cout<<G4endl;
703
704 for(k=0;k<numOfMaterials;k++)
705 {
706 G4cout <<k<< " Material : " <<(*theMaterialTable)[k]->GetName() << G4endl;
707 }
708 G4cout<<"Enter material name for dE/dx-distribution : "<<std::flush;
709 G4cin>>testName;
710 G4cout<<G4endl;
711
712 G4int iLoss, iStat, iStatMax, nGamma;
713 G4double energyLoss[50], Ebin, delta, delta1, delta2, delta3, step, y, pos;
714 G4double intProb[200], colDist, sum, fact, GF, lambda, aaa;
715
716 G4double alphaCrossTalk = -0.055, betaS = 0.2*0.4*keV;
717 G4int spectrum[50];
718
719 G4cout << " Enter nGamma 1<nGamma<10 : " <<std::flush;
720 G4cin>>nGamma;
721 G4cout<<G4endl;
722
723 for(k=0;k<numOfMaterials;k++)
724 {
725 if((*theMaterialTable)[k]->GetName() != testName) continue;
726
727 G4cout << "Material : " <<(*theMaterialTable)[k]->GetName() << G4endl<<G4endl;
728
729
730 G4cout << " Enter Lorentz factor : " <<std::flush;
731 G4cin>>gamma;
732 G4cout<<G4endl;
733
734 G4cout << " Enter step in mm : " <<std::flush;
735 G4cin>>step;
736 G4cout<<G4endl;
737 step *= mm;
738
739 G4cout << " Enter energy bin in keV : " <<std::flush;
740 G4cin>>Ebin;
741 G4cout<<G4endl;
742 Ebin *= keV;
743
744 G4cout << " Enter number of events : " <<std::flush;
745 G4cin>>iStatMax;
746
747 G4cout<<G4endl<<"Start dE/dx distribution"<<G4endl<<G4endl;
748
749 maxEnergyTransfer = 100*keV;
750 bg2 = gamma*gamma - 1;
751 rateMass = electron_mass_c2/proton_mass_c2;
752
753 Tmax = 2.0*electron_mass_c2*bg2
754 /(1.0+2.0*gamma*rateMass+rateMass*rateMass);
755
756 if ( maxEnergyTransfer > Tmax) maxEnergyTransfer = Tmax;
757
758 G4PAIxSection testPAIenergyLoss(k,maxEnergyTransfer,bg2);
759
760 for( iLoss = 0; iLoss < 50; iLoss++ )
761 {
762 energyLoss[iLoss] = Ebin*iLoss;
763 spectrum[iLoss] = 0;
764 }
765 for(iStat=0;iStat<iStatMax;iStat++)
766 {
767
768 // aaa = (G4double)nGamma;
769 // lambda = aaa/step;
770 // colDist = RandGamma::shoot(aaa,lambda);
771
772 // delta = testPAIenergyLoss.GetStepEnergyLoss(colDist);
773
774 // delta = testPAIenergyLoss.GetStepEnergyLoss(step);
775
776 delta1 = testPAIenergyLoss.GetStepEnergyLoss(step);
777
778 delta = G4RandGauss::shoot(delta1,0.3*delta1);
779 if( delta < 0.0 ) delta = 0.0;
780
781 // delta2 = testPAIenergyLoss.GetStepEnergyLoss(step);
782 // delta3 = testPAIenergyLoss.GetStepEnergyLoss(step);
783
784 // delta = alphaCrossTalk*delta1 +
785 // delta2 + alphaCrossTalk*delta3 - betaS;
786
787 for(iLoss=0;iLoss<50;iLoss++)
788 {
789 if(delta <= energyLoss[iLoss]) break;
790 }
791 spectrum[iLoss-1]++;
792 }
793 G4double meanLoss = 0.0;
794
795 outFile<<"E, keV"<<"\t\t"<<"Distribution"<<G4endl<<G4endl;
796 G4cout<<"E, keV"<<"\t\t"<<"Distribution"<<G4endl<<G4endl;
797 G4cout<<G4endl;
798 for(iLoss=0;iLoss<50;iLoss++) // with last bin
799 {
800 fileOut<<energyLoss[iLoss]/keV<<"\t\t"<<spectrum[iLoss]<<G4endl;
801 G4cout<<energyLoss[iLoss]/keV<<"\t\t"<<spectrum[iLoss]<<G4endl;
802 meanLoss +=energyLoss[iLoss]*spectrum[iLoss];
803 }
804 G4cout<<G4endl;
805 G4cout<<"Mean loss over spectrum = "<<meanLoss/keV/iStatMax<<" keV"<<G4endl;
806 }
807
808 G4int exit = 1;
809
810 while(exit)
811 {
812 G4cout<<"Enter 'y' , if you would like to compare with exp. data : "<<std::flush;
813 G4cin>>confirm;
814 if(confirm != "y" ) break;
815 G4cout<<G4endl;
816
817 // Read experimental data file
818
819 G4double delExp[200], distr[200], deltaBin, sumPAI, sumExp;
820 G4int numberOfExpPoints;
821
822 G4cout<<G4endl;
823 G4cout << " Enter number of experimental points : " <<std::flush;
824 G4cin>>numberOfExpPoints;
825 G4cout<<G4endl;
826 G4cout << " Enter energy bin in keV : " <<std::flush;
827 G4cin>>deltaBin;
828 G4cout<<G4endl;
829 deltaBin *= keV;
830
831 std::ifstream fileRead;
832 fileRead.open("input.dat");
833 for(i=0;i<numberOfExpPoints;i++)
834 {
835 fileRead>>delExp[i]>>distr[i];
836 delExp[i] *= keV;
837 G4cout<<i<<"\t"<<delExp[i]<<"\t"<<distr[i]<<G4endl;
838 }
839 fileRead.close();
840
841 // Adjust statistics of experiment to PAI simulation
842
843 sumExp = 0.0;
844 for(i=0;i<numberOfExpPoints;i++) sumExp +=distr[i];
845 sumExp *= deltaBin;
846
847 sumPAI = 0.0;
848 for(i=0;i<49;i++) sumPAI +=spectrum[i];
849 sumPAI *= Ebin;
850
851 for(i=0;i<numberOfExpPoints;i++) distr[i] *= sumPAI/sumExp;
852
853 for(i=0;i<numberOfExpPoints;i++)
854 {
855 fileWrite<<delExp[i]/keV<<"\t"<<distr[i]<<G4endl;
856 G4cout<<delExp[i]/keV<<"\t"<<distr[i]<<G4endl;
857 }
858 exit = 0;
859 }
860
861 G4cout<<"Enter 'y' , if you would like to get most probable delta : "<<std::flush;
862 G4cin>>confirm;
863 if(confirm != "y" ) return 1;
864 G4cout<<G4endl;
865
866 G4int kGamma, iMPLoss, maxSpectrum, iMax;
867 G4double mpDelta[50], meanDelta[50], rrMP[50], rrMean[50];
868 G4double mpLoss, tmRatio, mpSum, mpStat;
869
870 G4double aGamma[33] =
871 {
872 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
873 20., 40.0, 60.0, 80.0, 100.0, 200.0, 400.0, 600.0, 800.0, 1000.0, // 23
874 2000.0, 4000.0, 6000.0, 8000.0, 100000.0, 20000.0, // 29
875 40000.0, 60000.0, 80000.0, 100000.0 // 33
876 };
877
878 for(k=0;k<numOfMaterials;k++)
879 {
880 G4cout <<k<< " Material : " <<(*theMaterialTable)[k]->GetName() << G4endl;
881 }
882 G4cout<<"Enter material name for dE/dx-distribution : "<<std::flush;
883 G4cin>>testName;
884 G4cout<<G4endl;
885
886
887 for(k=0;k<numOfMaterials;k++)
888 {
889 if((*theMaterialTable)[k]->GetName() != testName) continue;
890
891 G4cout << "Material : " <<(*theMaterialTable)[k]->GetName() << G4endl<<G4endl;
892
893 G4cout << " Enter nGamma 1<nGamma<10 : " <<std::flush;
894 G4cin>>nGamma;
895 G4cout<<G4endl;
896
897
898 G4cout << " Enter step in mm : " <<std::flush;
899 G4cin>>step;
900 G4cout<<G4endl;
901 step *= mm;
902
903 G4cout << " Enter energy bin in keV : " <<std::flush;
904 G4cin>>Ebin;
905 G4cout<<G4endl;
906 Ebin *= keV;
907
908 G4cout << " Enter trancated mean ration <1.0 : " <<std::flush;
909 G4cin>>tmRatio;
910 G4cout<<G4endl;
911
912
913 G4cout << " Enter number of events : " <<std::flush;
914 G4cin>>iStatMax;
915 G4cout<<G4endl;
916
917 G4cout<<"no."<<"\t"<<"Gamma"<<"\t"<<"Rel. rise"<<"\t"<<"M.P. loss, keV"
918 <<"\t"<<"Mean loss, keV"<<G4endl<<G4endl;
919 // outFile<<"no."<<"\t"<<"Gamma"<<"\t"<<"M.P. loss, keV"
920 // <<"\t"<<"Mean loss, keV"<<G4endl<<G4endl;
921
922
923 // gamma = 1.1852;
924
925 for(kGamma=0;kGamma<33;kGamma++)
926 {
927 // G4cout<<G4endl<<"Start dE/dx distribution"<<G4endl<<G4endl;
928
929 gamma = aGamma[kGamma];
930 maxEnergyTransfer = 100*keV;
931 bg2 = gamma*gamma - 1;
932 rateMass = electron_mass_c2/proton_mass_c2;
933
934 Tmax = 2.0*electron_mass_c2*bg2
935 /(1.0+2.0*gamma*rateMass+rateMass*rateMass);
936
937 if ( maxEnergyTransfer > Tmax) maxEnergyTransfer = Tmax;
938
939 G4PAIxSection testPAIenergyLoss(k,maxEnergyTransfer,bg2);
940
941 for( iLoss = 0; iLoss < 50; iLoss++ )
942 {
943 energyLoss[iLoss] = Ebin*iLoss;
944 spectrum[iLoss] = 0;
945 }
946 for(iStat=0;iStat<iStatMax;iStat++)
947 {
948
949 // aaa = (G4double)nGamma;
950 // lambda = aaa/step;
951 // colDist = RandGamma::shoot(aaa,lambda);
952
953 // delta = testPAIenergyLoss.GetStepEnergyLoss(colDist);
954
955 delta = testPAIenergyLoss.GetStepEnergyLoss(step);
956
957 // delta1 = testPAIenergyLoss.GetStepEnergyLoss(step);
958 // delta2 = testPAIenergyLoss.GetStepEnergyLoss(step);
959 // delta3 = testPAIenergyLoss.GetStepEnergyLoss(step);
960
961 // delta = alphaCrossTalk*delta1 +
962 // delta2 + alphaCrossTalk*delta3 - betaS;
963
964 for(iLoss=0;iLoss<50;iLoss++)
965 {
966 if(delta <= energyLoss[iLoss]) break;
967 }
968 spectrum[iLoss-1]++;
969 }
970 G4int sumStat = 0;
971 for(iLoss=0;iLoss<49;iLoss++) // without last bin
972 {
973 sumStat += spectrum[iLoss];
974 if( sumStat > tmRatio*iStatMax ) break;
975 }
976 if(iLoss == 50) iLoss--;
977 iMPLoss = iLoss;
978 G4double meanLoss = 0.0;
979 maxSpectrum = 0;
980
981 for(iLoss=0;iLoss<iMPLoss;iLoss++) // without last bin
982 {
983 // fileOut<<energyLoss[iLoss]/keV<<"\t\t"<<spectrum[iLoss]<<G4endl;
984 // G4cout<<energyLoss[iLoss]/keV<<"\t\t"<<spectrum[iLoss]<<G4endl;
985
986 meanLoss += energyLoss[iLoss]*spectrum[iLoss];
987
988 if( spectrum[iLoss] > maxSpectrum )
989 {
990 maxSpectrum = spectrum[iLoss] ;
991 mpLoss = energyLoss[iLoss];
992 iMax = iLoss;
993 }
994 }
995 mpSum = 0.;
996 mpStat = 0;
997 for(iLoss = iMax-5;iLoss<=iMax+5;iLoss++)
998 {
999 mpSum += energyLoss[iLoss]*spectrum[iLoss];
1000 mpStat += spectrum[iLoss];
1001 }
1002 mpLoss = mpSum/mpStat;
1003 mpLoss /= keV;
1004 meanLoss /= keV*sumStat;
1005 meanDelta[kGamma] = meanLoss;
1006 mpDelta[kGamma] = mpLoss;
1007
1008 if(kGamma > 0)
1009 {
1010 rrMP[kGamma] = mpLoss/mpDelta[0];
1011 G4cout<<kGamma<<"\t"<<gamma<<"\t"<<rrMP[kGamma]<<"\t"<<mpLoss<<G4endl;
1012 // outFile<<gamma<<"\t"<<rrMP[kGamma]<<G4endl;
1013 fileWrite1<<gamma<<"\t"<<rrMP[kGamma]<<G4endl;
1014 }
1015
1016 // gamma *= 1.5;
1017 }
1018 G4cout<<G4endl;
1019 outFile<<G4endl;
1020 }
1021
1022 return EXIT_SUCCESS;
1023
1024}
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
Note: See TracBrowser for help on using the repository browser.