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

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

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