StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StiElossCalculator.cxx
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <cmath>
4 #include <assert.h>
5 #include "Sti/StiElossCalculator.h"
6 #include "StiUtilities/StiDebug.h"
7 #include "Stiostream.h"
8 static double gsigma2(double ZoverA,double DENS,double CHARGE2
9  ,double AMASS ,double BET2,double STEP );
10 static double gdrelx (double A ,double Z ,double DENS ,double T,double HMASS);
11 
12 //______________________________________________________________________________
13 StiElossCalculator::StiElossCalculator(double zOverA, double ionization, double A, double Z, double Dens)
14  : _zOverA(zOverA), _ionization2(ionization*ionization), _A(A), _Z(Z), _Dens(Dens)
15 {
16 static int nCall=0; nCall++;
17  assert(_A<=0 || _Z >0);
18  mId=nCall;
19 }
20 //______________________________________________________________________________
21 StiElossCalculator::StiElossCalculator()
22  : _zOverA(0), _ionization2(0), _A(0), _Z(0), _Dens(0)
23 {
24 static int nCall=0; nCall++;
25  mId=-nCall;
26 }
27 
28 //______________________________________________________________________________
29 void StiElossCalculator::set(double zOverA, double ionization, double A, double Z, double Dens)
30 {
31  _zOverA = zOverA; _ionization2 = ionization*ionization; _A = A; _Z = Z; _Dens = Dens;
32 }
33 //______________________________________________________________________________
34 
38 //______________________________________________________________________________
39 StiElossCalculator::~StiElossCalculator(){}
40 
55 //______________________________________________________________________________
56 double StiElossCalculator::calculate(double z2, double m, double beta2) const
57 {
58 static int nCall=0; nCall++;
59 static int noEloss = StiDebug::iFlag("NOELOSS");
60 if (noEloss) return 0;
61 
62  if (_A <=0.) return 0.;
63  double beta21 = 1 - beta2; if (beta21 < 1.e-10) beta21 = 1.e-10;
64  double T = m*(1./::sqrt(beta21) - 1);
65  double dedx = gdrelx(_A,_Z,_Dens,T,m)*z2*_Dens;
66 assert(!std::isnan(dedx));
67 assert(dedx>=0 && dedx<1e3);
68  return dedx;
69 }
70 //______________________________________________________________________________
71 double StiElossCalculator::calcError(double z2, double m, double beta2) const
72 {
73 static int noEloss = StiDebug::iFlag("NOELOSS");
74 if (noEloss) return 0;
75 
76  if (_A<=0.) return 0.;
77  double err2=gsigma2(_Z/_A,1.,z2,m ,beta2,1.);
78 assert(!std::isnan(err2));
79 assert(err2>=0 && err2<1e3);
80  return err2;
81 }
82 
83 //______________________________________________________________________________
84 ostream& operator<<(ostream& os, const StiElossCalculator& m) {
85  os << "zOverA "<< m.getzOverA()
86  << "\tionization2 "<< m.getionization2()
87  << "\tA "<< m.getA()
88  << "\tZ "<< m.getZ()
89  << "\tDens"<< m.getDens() << endl;
90  return os;
91 }
92 //______________________________________________________________________________
93 //* Revision 1.1.1.1 1995/10/24 10:21:24 cernlib
94 //* Geant
95 //*
96 //*
97 //#include "geant321/pilot.h"
98 //*CMZ : 3.21/03 06/10/94 18.22.43 by S.Giani
99 //*-- Author :
100 //______________________________________________________________________________
101 // SUBROUTINE GDRELX(A,Z,DENS,T,HMASS,dedx)
102 double gdrelx(double A,double Z,double DENS,double T,double HMASS)
103 {
104 //
105 // ******************************************************************
106 // * *
107 // * Calculates the mean 1/DENS*dE/dx of a particle with kinetic *
108 // * energy T in an element of atomic number Z, atomic weight A *
109 // * and density DENS ( the density is just used for the *
110 // * calculation of the density effect in the case of high T). *
111 // * The routine reproduces the experimental and/or tabulated *
112 // * energy losses rather well down to T -> 0. *
113 // * Simple parametrization is used for T .le. T2L=2 MeV (see *
114 // * H.H.Andersen,J.F.Ziegler:Hydrogen stopping powers and *
115 // * ranges in all element,Pergamon Press,1977.). *
116 // * For T .gt. T2L=2 MeV the corrected Bethe-Bloch stopping *
117 // * power / restricted energy loss formula is used. *
118 // * *
119 // * *
120 // * ==>Called by : GDRELA *
121 // * Author L.Urban ********* *
122 // * *
123 // ******************************************************************
124 //
125 //#include "geant321/gconsp.inc"
126 //#include "geant321/gccuts.inc"
127 //#include "geant321/gcunit.inc"
128 static const double AMUKEV=931494.32,AMUKEV50=pow(AMUKEV,0.50),AMUKEV45=pow(AMUKEV,0.45);
129 static const double D=0.000153537,T1L=0.00001,T2L=0.002;
130 static const double AVO=0.60221367,EMPROT=0.9382723,EMASS=0.0005109990615;
131 static const double DCUTM=9999.;
132 //DIMENSION B(6,92),C(6,92),CECOF[6]
133 //*
134 static const double B[93][6]={
135  {0.},
136  {1.262 ,1.44 ,242.6 ,12000. ,0.1159 ,18.8},
137  {1.229 ,1.397 ,484.5 ,5873. ,0.05225 ,41.7},
138  {1.411 ,1.6 ,725.6 ,3013. ,0.04578 ,47.6},
139  {2.248 ,2.59 ,966. ,153.8 ,0.03475 ,62.7},
140  {2.474 ,2.815 ,1206. ,1060. ,0.02855 ,76.0},
141  {2.631 ,2.989 ,1445. ,957.2 ,0.02819 ,77.3},
142  {2.954 ,3.35 ,1683. ,1900. ,0.02513 ,86.7},
143  {2.652 ,3. ,1920. ,2000. ,0.0223 ,97.7},
144  {2.085 ,2.352 ,2157. ,2634. ,0.01816 ,120.},
145  {1.951 ,2.199 ,2393. ,2699. ,0.01568 ,139.},
146  {2.542 ,2.869 ,2628. ,1854. ,0.01472 ,148.},
147  {3.792 ,4.293 ,2862. ,1009. ,0.01397 ,156.},
148  {4.154 ,4.739 ,2766. ,164.5 ,0.02023 ,162.},
149  {4.15 ,4.7 ,3329. ,550. ,0.01321 ,165.},
150  {3.232 ,3.647 ,3561. ,1560. ,0.01267 ,172.},
151  {3.447 ,3.891 ,3792. ,1219. ,0.01211 ,180.},
152  {5.047 ,5.714 ,4023. ,878.6 ,0.01178 ,185.},
153  {5.731 ,6.5 ,4253. ,530. ,0.01123 ,194.},
154  {5.151 ,5.833 ,4482. ,545.7 ,0.01129 ,193.},
155  {5.521 ,6.252 ,4710. ,553.3 ,0.01112 ,196.},
156  {5.201 ,5.884 ,4938. ,560.9 ,0.009995 ,218.},
157  {4.862 ,5.496 ,5165. ,568.5 ,0.009474 ,230.},
158  {4.48 ,5.055 ,5391. ,952.3 ,0.009117 ,239.},
159  {3.983 ,4.489 ,5616. ,1336. ,0.008413 ,259.},
160  {3.469 ,3.907 ,5725. ,1461. ,0.008829 ,270.},
161  {3.519 ,3.963 ,6065. ,1243. ,0.007782 ,280.},
162  {3.14 ,3.535 ,6288. ,1372. ,0.007361 ,296.},
163  {3.553 ,4.004 ,6205. ,555.1 ,0.008763 ,310.},
164  {3.696 ,4.175 ,4673. ,387.8 ,0.02188 ,322.},
165  {4.21 ,4.75 ,6953. ,295.2 ,0.006809 ,320.},
166  {5.041 ,5.697 ,7173. ,202.6 ,0.006725 ,324.},
167  {5.554 ,6.3 ,6496. ,110. ,0.009689 ,330.},
168  {5.323 ,6.012 ,7611. ,292.5 ,0.006447 ,338.},
169  {5.874 ,6.656 ,7395. ,117.5 ,0.007684 ,340.},
170  {5.611 ,6.335 ,8046. ,365.2 ,0.006244 ,349.},
171  {6.411 ,7.25 ,8262. ,220. ,0.006087 ,358.},
172  {5.694 ,6.429 ,8478. ,292.9 ,0.006087 ,358.},
173  {6.339 ,7.159 ,8693. ,330.3 ,0.006003 ,363.},
174  {6.407 ,7.234 ,8907. ,367.8 ,0.005889 ,370.},
175  {6.734 ,7.603 ,9120. ,405.2 ,0.005765 ,378.},
176  {6.902 ,7.791 ,9333. ,442.7 ,0.005587 ,390.},
177  {6.425 ,7.248 ,9545. ,480.2 ,0.005367 ,406.},
178  {6.799 ,7.671 ,9756. ,517.6 ,0.005315 ,410.},
179  {6.108 ,6.887 ,9966. ,555.1 ,0.005151 ,423.},
180  {5.924 ,6.677 ,10180. ,592.5 ,0.004919 ,443.},
181  {5.238 ,5.9 ,10380. ,630. ,0.004758 ,458.},
182  {5.623 ,6.354 ,7160. ,337.6 ,0.01394 ,466.},
183  {5.814 ,6.554 ,10800. ,355.5 ,0.004626 ,471.},
184  {6.23 ,7.024 ,11010. ,370.9 ,0.00454 ,480.},
185  {6.41 ,7.227 ,11210. ,386.4 ,0.004474 ,487.},
186  {7.5 ,8.48 ,8608. ,348. ,0.009074 ,494.},
187  {6.979 ,7.871 ,11620. ,392.4 ,0.004402 ,495.},
188  {7.725 ,8.716 ,11830. ,394.8 ,0.004376 ,498.},
189  {8.231 ,9.289 ,12030. ,397.3 ,0.004384 ,497.},
190  {7.287 ,8.218 ,12230. ,399.7 ,0.004447 ,490.},
191  {7.899 ,8.911 ,12430. ,402.1 ,0.004511 ,483.},
192  {8.041 ,9.071 ,12630. ,404.5 ,0.00454 ,480.},
193  {7.489 ,8.444 ,12830. ,406.9 ,0.00442 ,493.},
194  {7.291 ,8.219 ,13030. ,409.3 ,0.004298 ,507.},
195  {7.098 ,8. ,13230. ,411.8 ,0.004182 ,521.},
196  {6.91 ,7.786 ,13430. ,414.2 ,0.004058 ,537.},
197  {6.728 ,7.58 ,13620. ,416.6 ,0.003976 ,548.},
198  {6.551 ,7.38 ,13820. ,419. ,0.003877 ,562.},
199  {6.739 ,7.592 ,14020. ,421.4 ,0.003863 ,564.},
200  {6.212 ,6.996 ,14210. ,423.9 ,0.003725 ,585.},
201  {5.517 ,6.21 ,14400. ,426.3 ,0.003632 ,600.},
202  {5.219 ,5.874 ,14600. ,428.7 ,0.003498 ,623.},
203  {5.071 ,5.706 ,14790. ,433. ,0.003405 ,640.},
204  {4.926 ,5.542 ,14980. ,433.5 ,0.003342 ,652.},
205  {4.787 ,5.386 ,15170. ,435.9 ,0.003292 ,662.},
206  {4.893 ,5.505 ,15360. ,438.4 ,0.003243 ,672.},
207  {5.028 ,5.657 ,15550. ,440.8 ,0.003195 ,682.},
208  {4.738 ,5.329 ,15740. ,443.2 ,0.003186 ,684.},
209  {4.574 ,5.144 ,15930. ,442.4 ,0.003144 ,693.},
210  {5.2 ,5.851 ,16120. ,441.6 ,0.003122 ,698.},
211  {5.07 ,5.704 ,16300. ,440.9 ,0.003082 ,707.},
212  {4.945 ,5.563 ,16490. ,440.1 ,0.002965 ,735.},
213  {4.476 ,5.034 ,16670. ,439.3 ,0.002871 ,759.},
214  {4.856 ,5.46 ,18320. ,438.5 ,0.002542 ,755.},
215  {4.308 ,4.843 ,17040. ,487.8 ,0.002882 ,756.},
216  {4.723 ,5.311 ,17220. ,537. ,0.002913 ,748.},
217  {5.319 ,5.982 ,17400. ,586.3 ,0.002871 ,759.},
218  {5.956 ,6.7 ,17800. ,677. ,0.00266 ,765.},
219  {6.158 ,6.928 ,17770. ,586.3 ,0.002812 ,775.},
220  {6.204 ,6.979 ,17950. ,586.3 ,0.002776 ,785.},
221  {6.181 ,6.954 ,18120. ,586.3 ,0.002748 ,793.},
222  {6.949 ,7.82 ,18300. ,586.3 ,0.002737 ,796.},
223  {7.506 ,8.448 ,18480. ,586.3 ,0.002727 ,799.},
224  {7.649 ,8.609 ,18660. ,586.3 ,0.002697 ,808.},
225  {7.71 ,8.679 ,18830. ,586.3 ,0.002641 ,825.},
226  {7.407 ,8.336 ,19010. ,586.3 ,0.002603 ,837.},
227  {7.29 ,8.204 ,19180. ,586.3 ,0.002573 ,847.}};
228 
229 static const double CECOF[7]={0.,0.42237,0.0304,-0.00038,3.858,-0.1668,0.00158};
230 double C[6]={0};
231 
232 double poti,p,e,beta,bet2,tau,sl,sh,eta,eta2,b2g2,tmax,cc,x0,x1,xa,xm,delta;
233 double f1,f2,f3,f4,f5,tupp,ce,st,sbb,dedx;
234 //* ------------------------------------------------------------------
235 //* in the case of non-integer Z the low energy parameters
236 //* and the ionization potential are taken at INT(Z) !
237 //*
238  int iz=(int)(Z+1e-8); if (iz == 92) iz = 91;
239  double wt1=Z-iz,wt0 = 1-wt1;
240  assert((iz>0)&&(iz<92));
241 //*
242 //* Calculate coefficients C(I,J) if it has not been done already
243 //*
244  double fac=AVO/A;
245  C[0]=fac*AMUKEV50*(B[iz][0]*wt0+B[iz+1][0]*wt1);
246  C[1]=fac*AMUKEV45*(B[iz][1]*wt0+B[iz+1][1]*wt1);
247  C[2]=fac* (B[iz][2]*wt0+B[iz+1][2]*wt1)/AMUKEV;
248  C[3]= (B[iz][3]*wt0+B[iz+1][3]*wt1)/AMUKEV;
249  C[4]=AMUKEV* (B[iz][4]*wt0+B[iz+1][4]*wt1);
250 //* poti=16.E-9*Z**0.9
251  C[5]= (B[iz][5]*wt0+B[iz+1][5]*wt1)*1.E-9;
252 //*
253 //* ----------------------------------------------------------------
254  double hmass2 = HMASS*HMASS;
255  double T1LIM=HMASS*T1L/EMPROT;
256  double T2LIM=HMASS*T2L/EMPROT;
257 //*
258 //* Calculate dE/dx
259 //* ---> for T .le. T1LIM (very low energy)
260 //*
261  if(T<=T1LIM) {
262  tau=T/HMASS;
263  dedx=C[0]*pow(tau,0.5);
264  } else {
265 //*
266 //* ---> for T1LIM .lt. T and T .le. T2LIM (low energy)
267 //*
268  if(T<=T2LIM) {
269  tau=T/HMASS;
270  sl=C[1]*pow(tau,0.45);
271  sh=C[2]*log(1.+C[3]/tau+C[4]*tau)/tau;
272  dedx=sl*sh/(sl+sh);
273 //*
274 //* ---> for T .gt. T2LIM ( "high " energy , Bethe-Bloch formula)
275 //*
276  } else {
277  p=sqrt(T*(T+2.*HMASS));
278  e=T+HMASS;
279  beta=p/e;
280  bet2=beta*beta;
281  eta=p/HMASS;
282  eta2=eta*eta;
283 //*+++ new line follows.....
284  b2g2=eta*eta;
285 //*+++ end of correction
286  tmax=2.*EMASS*T*(T+2.*HMASS);
287 //*+++ correction of the next line
288 //* tmax=tmax/(HMASS**2+EMASS**2+EMASS*(T+HMASS));
289  tmax=tmax/(hmass2+EMASS*(EMASS+2.*e));
290 //*+++ end of correction
291 //*
292 //* density correction
293 //*
294  poti=C[5];
295  cc=1.+2.*log(poti/(28.8E-9*sqrt(DENS*Z/A)));
296 //* condensed material ? ( dens .gt. 0.05 ? )
297  if(DENS>0.05) {
298  if(poti < 1.E-7) {
299  if(cc < 3.681) {
300  x0=0.2;
301  } else {
302  x0=0.326*cc-1.;
303  }
304  x1=2.;
305  } else {
306  if(cc < 5.215) {
307  x0=0.2;
308  } else {
309  x0=0.326*cc-1.5;
310  }
311  x1=3.;
312  }
313 //* gas ? ( dens . le . 0.05 ? )
314  } else {
315  if(cc<=12.25) {
316  int ip=(int)((cc-10.)/0.5)+1;
317  if(ip < 0) ip=0;
318  if(ip>4) ip=4;
319  x0=1.6+0.1*ip;
320  x1=4.;
321  } else {
322  if(cc<=13.804) {
323  x0=2.;
324  x1=5.;
325  } else {
326  x0=0.326*cc-2.5;
327  x1=5.;
328  }
329  }
330  }
331 //*
332  xa=cc/4.606;
333  xm=3.;
334  double aa=4.606*(xa-x0)/pow(x1-x0,xm);
335 //*
336  double x=log10(eta);
337  delta=0.;
338  if(x>x0) {
339  delta=4.606*x-cc;
340  if(x < x1) delta=delta+aa*pow(x1-x,xm);
341  }
342 //*
343 //* calculate shell correction
344 //*
345  double potsq=poti*poti;
346  if(eta>0.13) {
347  f1=1./eta2;
348  f2=f1*f1;
349  f3=f1*f2;
350  f4=(f1*CECOF[1]+f2*CECOF[2]+f3*CECOF[3])*1.E+12;
351  f5=(f1*CECOF[4]+f2*CECOF[5]+f3*CECOF[6])*1.E+18;
352  ce=f4*potsq+f5*potsq*poti;
353  } else {
354  eta2=0.0169;
355  f1=1./eta2;
356  f2=f1*f1;
357  f3=f1*f2;
358  f4=(f1*CECOF[1]+f2*CECOF[2]+f3*CECOF[3])*1.E+12;
359  f5=(f1*CECOF[4]+f2*CECOF[5]+f3*CECOF[6])*1.E+18;
360  ce=f4*potsq+f5*potsq*poti;
361  ce=ce*log(T/T2LIM)/log(0.0079/T2LIM);
362  }
363 //*
364  f1=D*Z/(A*bet2);
365 //*
366 //* stopping power or restricted dE/dx ?
367 //*
368 //*+++ correction of the next few lines
369 //* if(DCUTM.GE.tmax) {
370 //* f2=2.*(log(tmax/poti)-bet2)
371 //* } else {
372 //* f2=log(tmax*DCUTM/potsq)-bet2*(1.+DCUTM/tmax)
373 //* }
374  tupp=DCUTM;
375  if(tmax < DCUTM) tupp=tmax;
376  f2=log(2.*EMASS*b2g2/poti)+log(tupp/poti)-bet2*(1.+tupp/tmax);
377 //*+++ end of correction
378  dedx=f1*(f2-delta-2.*ce/Z);
379 //*
380 //*
381  tau=T2LIM/HMASS;
382  sl=C[1]*pow(tau,0.45);
383  sh=C[2]*log(1.+C[3]/tau+C[4]*tau)/tau;
384  st=sl*sh/(sl+sh);
385 //*
386  tmax=2.*EMASS*T2LIM*(T2LIM+2.*HMASS);
387 //*+++ correction of the next line
388 //* tmax=tmax/(HMASS**2+EMASS**2+EMASS*(T2LIM+HMASS));
389  tmax=tmax/(HMASS*HMASS+EMASS*EMASS+2.*EMASS*(T2LIM+HMASS));
390 //*+++ end of correction
391  bet2=T2LIM*(T2LIM+2.*HMASS)/pow(T2LIM+HMASS,2);
392  sbb=2.*(log(tmax/poti)-bet2);
393  sbb=D*Z*sbb/(A*bet2);
394  double corbb=(st/sbb-1.)*T2LIM;
395 //*
396  dedx=dedx*(1.+corbb/T);
397 //*
398  }
399  }
400  return dedx;
401 }
402 //*CMZ : 3.21/02 29/03/94 15.41.21 by S.Giani
403 //-- Author :
404 //______________________________________________________________________________
405 double gsigma2(double ZoverA,double DENS,double CHARGE2
406  ,double AMASS ,double BET2,double STEP )
407 {
408 // SUBROUTINE GFLUCT(DEMEAN,DE)
409 
410 static const double DGEV=0.153536E-3,EMASS=0.0005109990615;
411 //
412  double gamm2 = 1./(1.-BET2);
413  double gamma = sqrt(gamm2);
414 //
415 // *** low energy transfer
416  double xi = DGEV*CHARGE2*STEP*DENS*ZoverA/(BET2);
417 //
418 // Energy straggling using Gaussian
419 // STEP = current step-length (cm)
420 // Author : G.N. Patrick
421 //
422 //
423 // Maximum energy transfer to atomic electron (GeV).
424  double eta2 = BET2*gamm2;
425  double ratio = EMASS/AMASS;
426  double emax =(2*EMASS*eta2)/(1+2*ratio*gamma+ratio*ratio);
427 //
428 // +-----------------------------------+
429 // I Sample from Gaussian distribution I
430 // +-----------------------------------+
431  double sigma2 = xi*(1.-0.5*BET2)*emax;
432  return sigma2;
433 }
StiElossCalculator::_ionization2
double _ionization2
square of the ionization potential.
Definition: StiElossCalculator.h:49
StiElossCalculator::calculate
double calculate(double charge2, double m, double beta2) const
Definition: StiElossCalculator.cxx:56
StiElossCalculator::_zOverA
double _zOverA
Ratio of Z to A of the scattering material.
Definition: StiElossCalculator.h:47
Generated by   doxygen 1.8.5