CbmRoot
CbmLitMaterialEffectsImp.cxx
Go to the documentation of this file.
1 
6 
10 
11 #include "TDatabasePDG.h"
12 #include "TParticlePDG.h"
13 
14 #include <cassert>
15 #include <cmath>
16 #include <cstdlib>
17 #include <iostream>
18 
20  : fMass(0.105), fDownstream(true), fIsElectron(false), fIsMuon(true) {}
21 
23 
25  const CbmLitMaterialInfo* mat,
26  int pdg,
27  bool downstream) {
28  if (mat->GetLength() * mat->GetRho() < 1e-10) { return kLITSUCCESS; }
29 
30  fDownstream = downstream;
31  TDatabasePDG* db = TDatabasePDG::Instance();
32  TParticlePDG* particle = db->GetParticle(pdg);
33  assert(particle != NULL);
34  fMass = particle->Mass();
35  fIsElectron = (std::abs(pdg) == 11) ? true : false;
36  fIsMuon = (std::abs(pdg) == 13) ? true : false;
37 
38  AddEnergyLoss(par, mat);
39 
40  // AddThinScatter(par, mat);
41  AddThickScatter(par, mat);
42 
43  return kLITSUCCESS;
44 }
45 
47  CbmLitTrackParam* par,
48  const CbmLitMaterialInfo* mat) const {
49  if (fIsElectron) {
50  litfloat radLength = mat->GetLength() / mat->GetRL();
51  litfloat t;
52 
53  if (!fDownstream) {
54  t = radLength;
55  } else {
56  t = -radLength;
57  }
58 
59  litfloat qp = par->GetQp();
60  qp *= std::exp(-t);
61  par->SetQp(qp);
62 
63  litfloat cov = par->GetCovariance(18);
64  cov += CalcSigmaSqQpElectron(par, mat);
65  par->SetCovariance(18, cov);
66  } else {
67  litfloat Eloss = EnergyLoss(par, mat);
68  par->SetQp(CalcQpAfterEloss(par->GetQp(), Eloss));
69 
70  litfloat cov = par->GetCovariance(18);
71  cov += CalcSigmaSqQp(par, mat);
72  par->SetCovariance(18, cov);
73  }
74 }
75 
77  CbmLitTrackParam* par,
78  const CbmLitMaterialInfo* mat) const {
79  if (mat->GetLength() < 1e-10) { return; }
80 
81  litfloat tx = par->GetTx();
82  litfloat ty = par->GetTy();
83  litfloat thickness = mat->GetLength(); //cm
84  litfloat thetaSq = CalcThetaSq(par, mat);
85 
86  litfloat t = 1 + tx * tx + ty * ty;
87 
88  litfloat Q33 = (1 + tx * tx) * t * thetaSq;
89  litfloat Q44 = (1 + ty * ty) * t * thetaSq;
90  litfloat Q34 = tx * ty * t * thetaSq;
91 
92  litfloat T23 = (thickness * thickness) / 3.0;
93  litfloat T2 = thickness / 2.0;
94 
95  litfloat D = (fDownstream) ? 1. : -1.;
96 
97  std::vector<litfloat> C = par->GetCovMatrix();
98 
99  C[0] += Q33 * T23;
100  C[1] += Q34 * T23;
101  C[2] += Q33 * D * T2;
102  C[3] += Q34 * D * T2;
103 
104  C[6] += Q44 * T23;
105  C[7] += Q34 * D * T2;
106  C[8] += Q44 * D * T2;
107 
108  C[11] += Q33;
109  C[12] += Q34;
110 
111  C[15] += Q44;
112 
113  par->SetCovMatrix(C);
114 }
115 
117  CbmLitTrackParam* par,
118  const CbmLitMaterialInfo* mat) const {
119  if (mat->GetLength() < 1e-10) { return; }
120  litfloat tx = par->GetTx();
121  litfloat ty = par->GetTy();
122  litfloat thetaSq = CalcThetaSq(par, mat);
123 
124  litfloat t = 1 + tx * tx + ty * ty;
125 
126  litfloat Q33 = (1 + tx * tx) * t * thetaSq;
127  litfloat Q44 = (1 + ty * ty) * t * thetaSq;
128  litfloat Q34 = tx * ty * t * thetaSq;
129 
130  std::vector<litfloat> C = par->GetCovMatrix();
131  C[11] += Q33;
132  C[15] += Q44;
133  C[12] += Q34;
134  par->SetCovMatrix(C);
135 }
136 
137 litfloat
139  const CbmLitMaterialInfo* mat) const {
140  litfloat p = std::abs(1. / par->GetQp()); //GeV
141  litfloat E = std::sqrt(fMass * fMass + p * p);
142  litfloat beta = p / E;
143  litfloat x = mat->GetLength(); //cm
144  litfloat X0 = mat->GetRL(); //cm
145  litfloat bcp = beta * p;
146  litfloat z = 1.;
147 
148  litfloat theta = 0.0136 * (1. / bcp) * z * std::sqrt(x / X0)
149  * (1. + 0.038 * std::log(x / X0));
150  return theta * theta;
151 }
152 
153 litfloat
155  const CbmLitMaterialInfo* mat) const {
156  litfloat length = mat->GetRho() * mat->GetLength();
157  return dEdx(par, mat) * length;
158  //return MPVEnergyLoss(par, mat);
159 }
160 
162  const CbmLitMaterialInfo* mat) const {
163  litfloat dedx = BetheBloch(par, mat);
164  // dedx += BetheHeitler(par, mat);
165  // if (fIsMuon) dedx += PairProduction(par, mat);
166  return dedx;
167 }
168 
169 litfloat
171  const CbmLitMaterialInfo* mat) const {
172  litfloat K = 0.000307075; // GeV * g^-1 * cm^2
173  litfloat z = (par->GetQp() > 0.) ? 1 : -1.;
174  litfloat Z = mat->GetZ();
175  litfloat A = mat->GetA();
176 
177  litfloat M = fMass;
178  litfloat p = std::abs(1. / par->GetQp()); //GeV
179  litfloat E = std::sqrt(M * M + p * p);
180  litfloat beta = p / E;
181  litfloat betaSq = beta * beta;
182  litfloat gamma = E / M;
183  litfloat gammaSq = gamma * gamma;
184 
185  litfloat I = CalcI(Z) * 1e-9; // GeV
186 
187  litfloat me = 0.000511; // GeV
188  litfloat ratio = me / M;
189  litfloat Tmax =
190  (2 * me * betaSq * gammaSq) / (1 + 2 * gamma * ratio + ratio * ratio);
191 
192  // density correction
193  litfloat dc = 0.;
194  if (p > 0.5) { // for particles above 1 Gev
195  litfloat rho = mat->GetRho();
196  litfloat hwp = 28.816 * std::sqrt(rho * Z / A) * 1e-9; // GeV
197  dc = std::log(hwp / I) + std::log(beta * gamma) - 0.5;
198  }
199 
200  return K * z * z * (Z / A) * (1. / betaSq)
201  * (0.5 * std::log(2 * me * betaSq * gammaSq * Tmax / (I * I)) - betaSq
202  - dc);
203 }
204 
206  const CbmLitTrackParam* par,
207  const CbmLitMaterialInfo* mat) const {
208  litfloat K = 0.000307075; // GeV * g^-1 * cm^2
209  //myf z = (par->GetQp() > 0.) ? 1 : -1.;
210  litfloat Z = mat->GetZ();
211  litfloat A = mat->GetA();
212 
213  litfloat me = 0.000511; // GeV;
214  litfloat p = std::abs(1. / par->GetQp()); //GeV
215  litfloat E = std::sqrt(me * me + p * p);
216  litfloat gamma = E / me;
217 
218  litfloat I = CalcI(Z) * 1e-9; // GeV
219 
220  if (par->GetQp() > 0) { // electrons
221  return K * (Z / A) * (std::log(2 * me / I) + 1.5 * std::log(gamma) - 0.975);
222  } else { //positrons
223  return K * (Z / A) * (std::log(2 * me / I) + 2. * std::log(gamma) - 1.);
224  }
225 }
226 
228  litfloat eloss) const {
229  litfloat massSq = fMass * fMass;
230  litfloat p = std::abs(1. / qp);
231  litfloat E = std::sqrt(p * p + massSq);
232  litfloat q = (qp > 0) ? 1. : -1.;
233  if (!fDownstream) { eloss *= -1.0; } // TODO check this
234  litfloat Enew = E - eloss;
235  litfloat pnew = (Enew > fMass) ? std::sqrt(Enew * Enew - massSq) : 0.;
236  if (pnew != 0) {
237  return q / pnew;
238  } else {
239  return 1e5;
240  }
241 
242  //if (!fDownstream) eloss *= -1.0;
243  //if (p > 0.) p -= eloss;
244  //else p += eloss;
245  //return 1./p;
246 }
247 
248 litfloat
250  const CbmLitMaterialInfo* mat) const {
251  litfloat P = std::abs(1. / par->GetQp()); // GeV
252  litfloat XMASS = fMass; // GeV
253  litfloat E = std::sqrt(P * P + XMASS * XMASS);
254  litfloat Z = mat->GetZ();
255  litfloat A = mat->GetA();
256  litfloat RHO = mat->GetRho();
257  litfloat STEP = mat->GetLength();
258  litfloat EMASS = 0.511 * 1e-3; // GeV
259 
260  litfloat BETA = P / E;
261  litfloat GAMMA = E / XMASS;
262 
263  // Calculate xi factor (KeV).
264  litfloat XI = (153.5 * Z * STEP * RHO) / (A * BETA * BETA);
265 
266  // Maximum energy transfer to atomic electron (KeV).
267  litfloat ETA = BETA * GAMMA;
268  litfloat ETASQ = ETA * ETA;
269  litfloat RATIO = EMASS / XMASS;
270  litfloat F1 = 2. * EMASS * ETASQ;
271  litfloat F2 = 1. + 2. * RATIO * GAMMA + RATIO * RATIO;
272  litfloat EMAX = 1e6 * F1 / F2;
273 
274  litfloat DEDX2 = XI * EMAX * (1. - (BETA * BETA / 2.)) * 1e-12;
275 
276  litfloat SDEDX = (E * E * DEDX2) / std::pow(P, 6);
277 
278  return std::abs(SDEDX);
279 }
280 
282  const CbmLitTrackParam* par,
283  const CbmLitMaterialInfo* mat) const {
284  litfloat x = mat->GetLength(); //cm
285  litfloat X0 = mat->GetRL(); //cm
286  return par->GetQp() * par->GetQp()
287  * (std::exp(-x / X0 * std::log(3.0) / std::log(2.0))
288  - std::exp(-2.0 * x / X0));
289 }
290 
292  // mean excitation energy in eV
293  if (Z > 16.) {
294  return 10 * Z;
295  } else {
296  return 16 * std::pow(Z, 0.9);
297  }
298 }
299 
300 litfloat
302  const CbmLitMaterialInfo* mat) const {
303  litfloat M = fMass; //GeV
304  litfloat p = std::abs(1. / par->GetQp()); // GeV
305  litfloat rho = mat->GetRho();
306  litfloat X0 = mat->GetRL();
307  litfloat me = 0.000511; // GeV
308  litfloat E = std::sqrt(M * M + p * p);
309  litfloat ratio = me / M;
310 
311  return (E * ratio * ratio) / (X0 * rho);
312 }
313 
314 litfloat
316  const CbmLitMaterialInfo* mat) const {
317  litfloat p = std::abs(1. / par->GetQp()); // GeV
318  litfloat M = fMass; //GeV
319  litfloat rho = mat->GetRho();
320  litfloat X0 = mat->GetRL();
321  litfloat E = std::sqrt(M * M + p * p);
322 
323  return 7e-5 * E / (X0 * rho);
324 }
325 
327  const CbmLitMaterialInfo* mat) const {
329  / mat->GetA();
330 }
331 
332 litfloat
334  const CbmLitMaterialInfo* mat) const {
335  litfloat M = fMass * 1e3; //MeV
336  litfloat p = std::abs(1. / par->GetQp()) * 1e3; // MeV
337 
338  // myf rho = mat->GetRho();
339  litfloat Z = mat->GetZ();
340  litfloat A = mat->GetA();
341  litfloat x = mat->GetRho() * mat->GetLength();
342 
343  litfloat I = CalcI(Z) * 1e-6; // MeV
344 
345  litfloat K = 0.307075; // MeV g^-1 cm^2
346  litfloat j = 0.200;
347 
348  litfloat E = std::sqrt(M * M + p * p);
349  litfloat beta = p / E;
350  litfloat betaSq = beta * beta;
351  litfloat gamma = E / M;
352  litfloat gammaSq = gamma * gamma;
353 
354  litfloat ksi = (K / 2.) * (Z / A) * (x / betaSq); // MeV
355 
356  // myf hwp = 28.816 * std::sqrt(rho*Z/A) * 1e-6 ; // MeV
357  // myf dc = std::log(hwp/I) + std::log(beta*gamma) - 0.5;
358  // dc *= 2;
359  litfloat dc = 0.;
360 
361  litfloat eloss = ksi
362  * (std::log(2 * M * betaSq * gammaSq / I) + std::log(ksi / I)
363  + j - betaSq - dc);
364 
365  return eloss * 1e-3; //GeV
366 }
exp
friend F32vec4 exp(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:134
CbmLitTrackParam.h
Data class for track parameters.
litfloat
double litfloat
Definition: CbmLitFloat.h:15
CbmLitMaterialEffectsImp::AddThinScatter
void AddThinScatter(CbmLitTrackParam *par, const CbmLitMaterialInfo *mat) const
Definition: CbmLitMaterialEffectsImp.cxx:116
CbmLitMaterialEffectsImp::fIsElectron
bool fIsElectron
Definition: CbmLitMaterialEffectsImp.h:95
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmLitTrackParam
Data class for track parameters.
Definition: CbmLitTrackParam.h:29
CbmLitMaterialEffectsImp::CalcI
litfloat CalcI(litfloat Z) const
Definition: CbmLitMaterialEffectsImp.cxx:291
CbmLitMaterialInfo::GetRL
litfloat GetRL() const
Definition: CbmLitMaterialInfo.h:57
CbmLitMaterialEffectsImp::~CbmLitMaterialEffectsImp
virtual ~CbmLitMaterialEffectsImp()
Destructor.
Definition: CbmLitMaterialEffectsImp.cxx:22
lit::CbmLitDefaultSettings::ENERGY_LOSS_CONST
static const litfloat ENERGY_LOSS_CONST
Definition: CbmLitDefaultSettings.h:20
CbmLitDefaultSettings.h
CbmLitTrackParam::GetQp
litfloat GetQp() const
Definition: CbmLitTrackParam.h:58
CbmLitTrackParam::GetCovariance
litfloat GetCovariance(int index) const
Definition: CbmLitTrackParam.h:60
CbmLitTrackParam::SetQp
void SetQp(litfloat qp)
Definition: CbmLitTrackParam.h:69
CbmLitMaterialEffectsImp::BetheBlochElectron
litfloat BetheBlochElectron(const CbmLitTrackParam *par, const CbmLitMaterialInfo *mat) const
Definition: CbmLitMaterialEffectsImp.cxx:205
kLITSUCCESS
@ kLITSUCCESS
Definition: CbmLitEnums.h:24
CbmLitMaterialInfo
Definition: CbmLitMaterialInfo.h:19
CbmLitMaterialEffectsImp::PairProduction
litfloat PairProduction(const CbmLitTrackParam *par, const CbmLitMaterialInfo *mat) const
Definition: CbmLitMaterialEffectsImp.cxx:315
CbmLitMaterialEffectsImp::fMass
litfloat fMass
Definition: CbmLitMaterialEffectsImp.h:94
CbmLitMaterialEffectsImp::MPVEnergyLoss
litfloat MPVEnergyLoss(const CbmLitTrackParam *par, const CbmLitMaterialInfo *mat) const
Definition: CbmLitMaterialEffectsImp.cxx:333
CbmLitMaterialEffectsImp::fIsMuon
bool fIsMuon
Definition: CbmLitMaterialEffectsImp.h:96
CbmLitMaterialEffectsImp::CbmLitMaterialEffectsImp
CbmLitMaterialEffectsImp()
Constructor.
Definition: CbmLitMaterialEffectsImp.cxx:19
CbmLitMaterialEffectsImp::dEdx
litfloat dEdx(const CbmLitTrackParam *par, const CbmLitMaterialInfo *mat) const
Definition: CbmLitMaterialEffectsImp.cxx:161
log
friend F32vec4 log(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:135
CbmLitMaterialEffectsImp::AddThickScatter
void AddThickScatter(CbmLitTrackParam *par, const CbmLitMaterialInfo *mat) const
Definition: CbmLitMaterialEffectsImp.cxx:76
CbmLitMaterialEffectsImp::Update
LitStatus Update(CbmLitTrackParam *par, const CbmLitMaterialInfo *mat, int pdg, bool downstream)
Inherited from CbmLitMaterialEffects.
Definition: CbmLitMaterialEffectsImp.cxx:24
CbmLitMaterialEffectsImp::EnergyLoss
litfloat EnergyLoss(const CbmLitTrackParam *par, const CbmLitMaterialInfo *mat) const
Definition: CbmLitMaterialEffectsImp.cxx:154
CbmLitMaterialEffectsImp::fDownstream
bool fDownstream
Definition: CbmLitMaterialEffectsImp.h:93
CbmLitMaterialEffectsImp::BetheBlochSimple
litfloat BetheBlochSimple(const CbmLitMaterialInfo *mat) const
Definition: CbmLitMaterialEffectsImp.cxx:326
CbmLitMaterialEffectsImp::CalcSigmaSqQp
litfloat CalcSigmaSqQp(const CbmLitTrackParam *par, const CbmLitMaterialInfo *mat) const
Definition: CbmLitMaterialEffectsImp.cxx:249
CbmLitMaterialInfo::GetLength
litfloat GetLength() const
Definition: CbmLitMaterialInfo.h:54
CbmLitMaterialEffectsImp::CalcQpAfterEloss
litfloat CalcQpAfterEloss(litfloat qp, litfloat eloss) const
Definition: CbmLitMaterialEffectsImp.cxx:227
CbmLitMaterialInfo::GetZ
litfloat GetZ() const
Definition: CbmLitMaterialInfo.h:63
CbmLitTrackParam::SetCovMatrix
void SetCovMatrix(const vector< litfloat > &C)
Definition: CbmLitTrackParam.h:71
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmLitMaterialInfo::GetA
litfloat GetA() const
Definition: CbmLitMaterialInfo.h:66
CbmLitTrackParam::GetCovMatrix
const vector< litfloat > & GetCovMatrix() const
Definition: CbmLitTrackParam.h:61
CbmLitTrackParam::GetTy
litfloat GetTy() const
Definition: CbmLitTrackParam.h:57
CbmLitMaterialEffectsImp::BetheBloch
litfloat BetheBloch(const CbmLitTrackParam *par, const CbmLitMaterialInfo *mat) const
Definition: CbmLitMaterialEffectsImp.cxx:170
CbmLitMaterialInfo.h
LitStatus
LitStatus
Definition: CbmLitEnums.h:23
CbmLitTrackParam::SetCovariance
void SetCovariance(int index, litfloat cov)
Definition: CbmLitTrackParam.h:74
CbmLitMaterialInfo::GetRho
litfloat GetRho() const
Definition: CbmLitMaterialInfo.h:60
CbmLitMaterialEffectsImp::BetheHeitler
litfloat BetheHeitler(const CbmLitTrackParam *par, const CbmLitMaterialInfo *mat) const
Definition: CbmLitMaterialEffectsImp.cxx:301
CbmLitMaterialEffectsImp::CalcSigmaSqQpElectron
litfloat CalcSigmaSqQpElectron(const CbmLitTrackParam *par, const CbmLitMaterialInfo *mat) const
Definition: CbmLitMaterialEffectsImp.cxx:281
CbmLitTrackParam::GetTx
litfloat GetTx() const
Definition: CbmLitTrackParam.h:56
CbmLitMaterialEffectsImp::AddEnergyLoss
void AddEnergyLoss(CbmLitTrackParam *par, const CbmLitMaterialInfo *mat) const
Definition: CbmLitMaterialEffectsImp.cxx:46
CbmLitMaterialEffectsImp::CalcThetaSq
litfloat CalcThetaSq(const CbmLitTrackParam *par, const CbmLitMaterialInfo *mat) const
Definition: CbmLitMaterialEffectsImp.cxx:138
CbmLitMaterialEffectsImp.h
Calculation of multiple scattering and energy loss.