140 if (mode == 0)
return e3 * xi1 / (xi1 + xi2);
142 static TDatabasePDG* ptrTDatabasePDG = TDatabasePDG::Instance();
143 const double mass = ptrTDatabasePDG->GetParticle(pdg)->Mass();
144 const double chrg = fabs(ptrTDatabasePDG->GetParticle(pdg)->Charge() / 3);
149 return e3 * me1 / me0;
152 if (mode == 1)
return e3 * me1 / (me1 + me2);
157 const double mom2 = mom * mom;
158 const double beta2 = mom2 / (mass * mass + mom2);
159 const double buf = 0.1535e-3 * dx * chrg * chrg / beta2;
162 const double xi1i = 1 / xi1;
163 const double xi2i = 1 / xi2;
164 const double lm1min = -me1 * xi1i;
165 const double lm1max = (e3 - me1) * xi1i;
166 double lm1 = (e3 * me1 / (me1 + me2) - me1) * xi1i;
167 lm1 = min(lm1, lm1max);
168 lm1 = max(lm1, lm1min);
173 const double lm2b = (e3 - me1 - me2) * xi2i;
174 const double r = xi1 * xi2i;
177 while (itry <= 2 || fabs(lm1 - lm1old) > 0.001) {
185 lm1 = max(lm1min, lm1);
186 lm1 = min(lm1max, lm1);
187 double lm2 = lm2b - r * lm1;
188 double emlm1 = exp(-lm1);
189 double emlm2 = exp(-lm2);
190 double dldlm1 = 1 - emlm1 - r + r * emlm2;
191 double dl2dlm12 = emlm1 + r * r * emlm2;
192 deltal = dldlm1 / dl2dlm12;
196 if (mode == 3)
return xi1 * lm1 + me1;
198 double lm2 = lm2b - r * lm1;
199 double lmin = 0.5 * (lm1 + lm2 + exp(-lm1) + exp(-lm2));
200 const double ymax = exp(-lmin);
204 const double lm1w = lm1max - lm1min;
213 lm1 = (e3 * me1 / (me1 + me2) - me1) * xi1i;
216 gRandom->RndmArray(2, urndm);
217 lm1 = lm1min + lm1w * urndm[0];
218 lm2 = lm2b - r * lm1;
219 double l = 0.5 * (lm1 + lm2 + exp(-lm1) + exp(-lm2));
224 return xi1 * lm1 + me1;
228 double A,
double I,
double rho)
const
243 const double psq = p * p;
244 const double masq = mass * mass;
245 const double tote = psq + masq;
246 const double beta2 = psq / tote;
247 const double bgam2 = psq / masq;
248 const double gamma =
sqrt(tote) / mass;
249 const double alft = 0.1535e6 * rho * dx * Z / A;
250 const double k = log(
m_massE * alft / I / I) + 0.891;
251 const double k1 = k - log(dx);
252 const double alngam = log(gamma);
253 const double alggam = log10(gamma);
258 if (alggam > 3.6) irgn = 1;
262 }
else if (k1 < 7.) {
263 if (alggam > 3.6) irgn = 1;
267 }
else if (k1 < 8.) {
268 if (alggam > 3.6) irgn = 1;
272 }
else if (k1 < 9.) {
273 if (alggam > 3.6) irgn = 1;
277 }
else if (k1 < 10.) {
278 if (alggam > 3.6) irgn = 1;
282 }
else if (k1 < 11.) {
283 if (alggam > 4.0) irgn = 1;
288 if (alggam > 4.0) irgn = 1;
294 double delta = aa + alggam * (bb + cc * alggam);
295 delta = max(delta, 0.);
296 delta = min(delta, fabs(log(bgam2)));
298 double val = zi * zi * alft * (k + 2 * alngam - beta2 - delta) / beta2;
double getMostProbabEDep(double p, double mass, double zi, double dx, double z, double a, double i, double rho) const
Return the energy deosite in the material.