CbmRoot
Loading...
Searching...
No Matches
LitAddMaterial.h
Go to the documentation of this file.
1/* Copyright (C) 2009-2013 GSI/JINR-LIT, Darmstadt/Dubna
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Andrey Lebedev [committer] */
4
21#ifndef LITADDMATERIAL_H_
22#define LITADDMATERIAL_H_
23
24#include "LitMath.h"
25#include "LitTrackParam.h"
26
27namespace lit
28{
29 namespace parallel
30 {
31
39 template<class T>
40 inline void LitAddMaterial(LitTrackParam<T>& par, T siliconThickness)
41 {
42 // Silicon properties
43 static const T SILICON_DENSITY = 2.33; // g*cm^-3
44 // static const T SILICON_A = 28.08855; // silicon atomic weight
45 // static const T SILICON_Z = 14.0; // silicon atomic number
46 static const T SILICON_Z_OVER_A = 0.4984; //2.006325; // Z/A ratio for silicon
47 static const T SILICON_RAD_LENGTH = 9.365; // cm
48 static const T SILICON_I = 173 * 1e-9; // mean excitation energy [GeV]
49 static const T SILICON_I_SQ = SILICON_I * SILICON_I; // squared mean excitation energy
50
51 static const T ZERO = 0.0, ONE = 1., TWO = 2.;
52 static const T mass = 0.105658369; // muon mass [GeV/c]
53 static const T massSq = 0.105658369 * 0.105658369; // muon mass squared
54 static const T C1_2 = 0.5, C1_3 = 1. / 3.;
55 static const T me = 0.000511; // Electron mass [GeV/c]
56 static const T ratio = me / mass;
57
58 T p = sgn(par.Qp) / par.Qp; // Momentum [GeV/c]
59 T E = sqrt(massSq + p * p);
60 T beta = p / E;
61 T betaSq = beta * beta;
62 T gamma = E / mass;
63 T gammaSq = gamma * gamma;
64
65 // Scale material thickness
66 T norm = sqrt(ONE + par.Tx * par.Tx + par.Ty * par.Ty);
67 T thickness = norm * siliconThickness;
68 T radThick = thickness / SILICON_RAD_LENGTH;
69 T sqrtRadThick = sqrt(radThick);
70 T logRadThick = log(radThick);
71
72 /*
73 * Energy loss corrections
74 */
75
76 // Bethe-Block
77 static const T K = 0.000307075; // GeV * g^-1 * cm^2
78 T Tmax = (2 * me * betaSq * gammaSq) / (ONE + TWO * gamma * ratio + ratio * ratio);
79
80 // density correction
81 T dc = ZERO;
82 if (p > 0.5) { // for particles above 1 Gev
83 static const T c7 = 28.816;
84 static const T c8 = 1e-9;
85 T hwp = c7 * sqrt(SILICON_DENSITY * SILICON_Z_OVER_A) * c8; // GeV
86 dc = log(hwp / SILICON_I) + log(beta * gamma) - C1_2;
87 }
88
89 T bbLoss = K * SILICON_Z_OVER_A * rcp(betaSq)
90 * (C1_2 * log(TWO * me * betaSq * gammaSq * Tmax / SILICON_I_SQ) - betaSq - dc);
91
92 // static const T bbc = 0.00354;
93 // T bbLoss = bbc * SILICON_Z_OVER_A;
94
95 // Bethe-Heitler
96 // T bhLoss = (E * ratio * ratio)/(mat.X0 * mat.Rho);
97 T bhLoss = ZERO;
98
99 // Pair production approximation
100 // static const T c3 = 7e-5;
101 // T ppLoss = c3 * E / (mat.X0 * mat.Rho);
102 T ppLoss = ZERO;
103
104 // Integrated value of the energy loss
105 T energyLoss = (bbLoss + bhLoss + ppLoss) * SILICON_DENSITY * thickness;
106
107 // Correct Q/p value due to energy loss
108 T Enew = E - energyLoss;
109 T pnew = sqrt(Enew * Enew - massSq);
110 par.Qp = sgn(par.Qp) * rcp(pnew);
111
112
113 // Calculate Q/p correction in the covariance matrix
114 T betanew = pnew / Enew;
115 T betaSqnew = betanew * betanew;
116 T gammanew = Enew / mass;
117 // T gammaSqnew = gammanew * gammanew;
118
119 // Calculate xi factor (KeV).
120 static const T c4 = 153.5;
121 T XI = (c4 * SILICON_Z_OVER_A * thickness * SILICON_DENSITY) / betaSqnew;
122
123 // Maximum energy transfer to atomic electron (KeV).
124 T etanew = betanew * gammanew;
125 T etaSqnew = etanew * etanew;
126 T F1 = TWO * me * etaSqnew;
127 T F2 = ONE + TWO * ratio * gammanew + ratio * ratio;
128 static const T c5 = 1e6;
129 T emax = c5 * F1 / F2;
130
131 static const T c6 = 1e-12;
132 T dedxSq = XI * emax * (ONE - C1_2 * betaSqnew) * c6;
133
134 T p2 = pnew * pnew;
135 T p6 = p2 * p2 * p2;
136 T qpCorr = (Enew * Enew * dedxSq) / p6;
137 par.C14 += qpCorr;
138 // end calculate Q/p correction in the covariance matrix
139
140 /*
141 * End energy loss corrections
142 */
143
144 /*
145 * Multiple scattering corrections
146 */
147
148 T tx = par.Tx;
149 T ty = par.Ty;
150 T bcp = betanew * pnew;
151 //T bcp = beta * p;
152 static const T c1 = 0.0136, c2 = 0.038;
153 T theta = c1 * rcp(bcp) * sqrtRadThick * (ONE + c2 * logRadThick);
154 T thetaSq = theta * theta;
155
156 T t = ONE + tx * tx + ty * ty;
157
158 T Q33 = (1 + tx * tx) * t * thetaSq;
159 T Q44 = (1 + ty * ty) * t * thetaSq;
160 T Q34 = tx * ty * t * thetaSq;
161
162 T T23 = thickness * thickness * C1_3;
163 T T2 = thickness * C1_2;
164
165 par.C0 += Q33 * T23;
166 par.C1 += Q34 * T23;
167 par.C2 += Q33 * T2;
168 par.C3 += Q34 * T2;
169
170 par.C5 += Q44 * T23;
171 par.C6 += Q34 * T2;
172 par.C7 += Q44 * T2;
173
174 par.C9 += Q33;
175 par.C10 += Q34;
176
177 par.C12 += Q44;
178
179 /*
180 * End multiple scattering corrections
181 */
182 }
183
184
192 template<class T>
193 inline void LitAddMaterialElectron(LitTrackParam<T>& par, T siliconThickness)
194 {
195 // Silicon properties
196 static const T SILICON_RAD_LENGTH = 9.365; // cm
197
198 static const T ZERO = 0.0, ONE = 1., TWO = 2., THREE = 3., INF = 1.e5;
199 static const T C1_2 = 0.5, C1_3 = 1. / 3.;
200
201 //scale material thickness
202 static const T C_LOG = log(THREE) / log(TWO);
203 T norm = sqrt(ONE + par.Tx * par.Tx + par.Ty * par.Ty);
204 T thickness = norm * siliconThickness;
205 T radThick = thickness / SILICON_RAD_LENGTH;
206 T sqrtRadThick = sqrt(radThick);
207 T logRadThick = log(radThick);
208
209 // no material thickness scaling
210 // T thickness = mat.Thickness;
211 // T radThick = mat.RadThick;
212 // T sqrtRadThick = mat.SqrtRadThick;
213 // T logRadThick = mat.LogRadThick;
214
215 /*
216 * Energy loss corrections
217 */
218 par.Qp *= exp(radThick);
219 // par.C14 += par.Qp * par.Qp * mat.ElLoss; // no thickness scaling
220 par.C14 += (exp(radThick * C_LOG) - exp(-TWO * radThick)); // thickness scaling
221 /*
222 * End of energy loss corrections
223 */
224
225 /*
226 * Multiple scattering corrections
227 */
228 T pnew = sgn(par.Qp) / par.Qp; // Momentum [GeV/c]
229 T betanew = ONE;
230 T tx = par.Tx;
231 T ty = par.Ty;
232 T bcp = betanew * pnew;
233 static const T c1 = 0.0136, c2 = 0.038;
234 T theta = c1 * rcp(bcp) * sqrtRadThick * (ONE + c2 * logRadThick);
235 T thetaSq = theta * theta;
236
237 T t = ONE + tx * tx + ty * ty;
238
239 T Q33 = (1 + tx * tx) * t * thetaSq;
240 T Q44 = (1 + ty * ty) * t * thetaSq;
241 T Q34 = tx * ty * t * thetaSq;
242
243 T T23 = thickness * thickness * C1_3;
244 T T2 = thickness * C1_2;
245
246 par.C0 += Q33 * T23;
247 par.C1 += Q34 * T23;
248 par.C2 += Q33 * T2;
249 par.C3 += Q34 * T2;
250
251 par.C5 += Q44 * T23;
252 par.C6 += Q34 * T2;
253 par.C7 += Q44 * T2;
254
255 par.C9 += Q33;
256 par.C10 += Q34;
257
258 par.C12 += Q44;
259
260 /*
261 * End multiple scattering corrections
262 */
263 }
264
265 } // namespace parallel
266} // namespace lit
267#endif /* LITADDMATERIAL_H_ */
sqrt
friend fvec sqrt(const fvec &a)
Definition KfSimdPseudo.h:118
exp
friend fvec exp(const fvec &a)
Definition KfSimdPseudo.h:121
log
friend fvec log(const fvec &a)
Definition KfSimdPseudo.h:122
LitMath.h
Useful math functions.
LitTrackParam.h
Track parameters data class.
lit::parallel::LitTrackParam
Track parameters data class.
Definition LitTrackParam.h:40
lit::parallel::rcp
fscal rcp(const fscal &a)
Returns reciprocal.
Definition LitMath.h:32
lit::parallel::sgn
fscal sgn(const fscal &a)
Returns sign of the input number.
Definition LitMath.h:44
lit::parallel::LitAddMaterial
void LitAddMaterial(LitTrackParam< T > &par, T siliconThickness)
Definition LitAddMaterial.h:40
lit::parallel::LitAddMaterialElectron
void LitAddMaterialElectron(LitTrackParam< T > &par, T siliconThickness)
Definition LitAddMaterial.h:193