1 package org.lcsim.recon.tracking.trfeloss;
2 /**
3 *Class for evaluating the energy loss
4 * according to the Bethe-Bloch equation.
5 *
6 *@author Norman A. Graf
7 *@version 1.0
8 *
9 */
10
11 public class DeDxBethe extends DeDx
12 {
13
14
15 private double _density;
16 private double _particleMass;
17
18 // use notations from Sternheimer ...
19 private double _a;
20 private double _ionizationPot;
21 private double _cA;
22 private double _X0;
23 private double _X1;
24 private double _m;
25 private double _C;
26
27 /**
28 *Construct an instance given a material density.
29 *
30 * @param density The density of the material.
31 */
32 public DeDxBethe(double density)
33 {
34 _density = density;
35 // assume that the ionizing particle is a pion and
36 // all parameters of the absorber (except density)
37 // are the average of Si, Be, Al and polyvinyl Cl.
38 _particleMass = 0.13957;
39 _a = 0.074315;
40 _ionizationPot = 127.726;
41 _cA = 0.3123;
42 _X0 = 0.0966;
43 _X1 = 2.45;
44 _C = -3.878;
45 _m = 2.878;
46 }
47
48 /**
49 *Return energy loss (dE/dx) for a given energy.
50 *
51 * @param energy The energy of the particle.
52 * @return The average energy lost by a particle of this energy.
53 */
54 public double dEdX(double energy)
55 {
56 double de_dx;
57 de_dx = _a*_density;
58 // need a protection ...
59 double gamma = energy/_particleMass;
60 double beta = Math.sqrt(gamma*gamma-1)/gamma;
61 double weakDep = Math.log(1.044e12*Math.pow(beta,4)*Math.pow(gamma,4)
62 /Math.pow(_ionizationPot,2));
63 weakDep -= 2*beta*beta;
64
65 // include density effect correction
66 double X = Math.log(beta*gamma)/2.3026; // recall log10(x) = log(x)/log(10)
67 double delta;
68 if(X<_X1 && X>_X0)
69 {delta = 4.6052*X + _cA*Math.pow((_X1-X),_m) + _C;}
70 else if(X>_X1)
71 {delta = 4.6052*X + _C;}
72 else
73 {delta = 0;}
74 weakDep -= delta;
75
76 de_dx *= weakDep; // MeV/cm
77 de_dx /= beta*beta;
78
79 de_dx /= 1000.0; // GeV/cm
80 return de_dx;
81 }
82
83 /**
84 *Return the uncertainty in the energy lost sigma(E).
85 *
86 * @param energy The energy of the particle.
87 * @param x The amount of material.
88 * @return The uncertainty in the energy lost sigma(E).
89 */
90 public double sigmaEnergy(double energy, double x)
91 {
92 double factor = 1.022;
93 double sigma_e = factor*_density*_a*x;
94
95 // add relativistic correction
96 double gamma = energy/_particleMass;
97 double beta = Math.sqrt(gamma*gamma - 1)/gamma;
98 sigma_e *= (1-0.5*beta*beta)/(1-beta*beta); // MeV^2
99 sigma_e = Math.sqrt(sigma_e)/1000; // GeV
100 return sigma_e;
101 }
102
103 /**
104 *Return new energy for a given path distance.
105 * Energy increases if x < 0.
106 *
107 * @param energy The energy of the particle.
108 * @param x The amount of material.
109 * @return New energy for a given path distance.
110 */
111 public double loseEnergy(double energy, double x)
112 {
113 double deltaE = dEdX(energy)*x;
114 return energy-deltaE;
115 }
116
117
118 /**
119 *output stream
120 *
121 * @return A String representation of this instance.
122 */
123 public String toString()
124 {
125 return "DeDxBethe";
126 }
127
128 }