14 #include "GaudiKernel/MsgStream.h"
15 #include "GaudiKernel/SystemOfUnits.h"
27 declareInterface<Trk::IEnergyLossUpdator>(
this);
42 double E = std::sqrt(
p *
p +
m *
m);
55 dEdX += -0.5345 /
mat.x0() + 6.803e-5 *
E /
mat.x0() +
56 2.278e-11 *
E *
E /
mat.x0() -
57 9.899e-18 *
E *
E *
E /
mat.x0();
59 dEdX += -2.986 /
mat.x0() + 9.253e-5 *
E /
mat.x0();
69 double pathcorrection,
76 "undefined ParticleHypothesis, energy loss calculation cancelled");
88 double pathLength = pathcorrection *
mat.thicknessInX0() *
mat.x0();
98 meanIoni =
sign * pathLength * meanIoni;
99 meanRad =
sign * pathLength * meanRad;
100 sigIoni = pathLength * sigIoni;
101 sigRad = pathLength * sigRad;
109 deltaE = meanIoni + meanRad;
111 double sigmaDeltaE = std::sqrt(sigIoni * sigIoni + sigRad * sigRad);
113 << deltaE <<
" meanIoni " << meanIoni <<
" meanRad " << meanRad
114 <<
" sigIoni " << sigIoni <<
" sigRad " << sigRad <<
" sign "
115 <<
sign <<
" pathLength " << pathLength);
118 sigmaDeltaE, meanIoni, sigIoni,
119 meanRad, sigRad, pathLength));
127 double caloEnergyError,
129 double momentumError,
130 int& elossFlag)
const
161 double deltaE = eLoss.
deltaE();
164 double deltaE_ioni = eLoss.
meanIoni();
165 double sigmaDeltaE_ioni = 0.45 * eLoss.
sigmaIoni();
166 double deltaE_rad = eLoss.
meanRad();
167 double sigmaDeltaE_rad =
173 if (eLoss.
meanRad() > 100000.) {
174 deltaE_rad = 100000.;
181 double MOP = deltaE_ioni - isign * 3.59524 * sigmaDeltaE_ioni;
187 isign * 50 * 10000. / pCaloEntry +
188 isign * 0.75 * std::sqrt(sigmaDeltaE_ioni * sigmaDeltaE_rad);
189 double MOPshiftNoRad = isign * 50 * 10000. / pCaloEntry;
193 double fracErad = sigmaDeltaE_rad +
194 std::abs(deltaE_rad) * pCaloEntry / (800000. + pCaloEntry);
195 double sigmaL = sigmaDeltaE_ioni + 0.8 * fracErad / 3.59524;
197 double sigmaMinus = 1.02 * sigmaDeltaE_ioni + 0.08 * sigmaDeltaE_rad;
198 double sigmaPlus = 4.65 * sigmaDeltaE_ioni + 1.16 * sigmaDeltaE_rad;
201 double xc = momentumError / (sigmaL > 0. ? sigmaL : 1.);
203 (0.3849 * xc * xc + 7.76672e-03 * xc * xc * xc) /
204 (1 + 2.8631 * xc + 0.3849 * xc * xc + 7.76672
e-03 * xc * xc * xc);
211 if (caloEnergyError <= 0) {
215 double MOPreso = isign * 3.59524 * sigmaL *
correction;
217 deltaE = MOP + MOPshift + MOPreso;
218 sigmaMinusDeltaE = sigmaMinus;
219 sigmaPlusDeltaE = sigmaPlus;
220 sigmaDeltaE = std::sqrt(0.5 * sigmaMinusDeltaE * sigmaMinusDeltaE +
221 0.5 * sigmaPlusDeltaE * sigmaPlusDeltaE);
223 if (m_optimalRadiation && std::abs(deltaE) < caloEnergy &&
224 pCaloEntry > 100000) {
233 sigmaL = sigmaDeltaE_ioni + 0.3 * fracErad / 3.59524;
234 xc = momentumError / (sigmaL > 0. ? sigmaL : 1.);
236 (0.3849 * xc * xc + 7.76672e-03 * xc * xc * xc) /
237 (1 + 2.8631 * xc + 0.3849 * xc * xc + 7.76672
e-03 * xc * xc * xc);
239 MOPreso = isign * 3.59524 * sigmaL *
correction;
240 deltaE = MOP + MOPshift + MOPreso;
241 sigmaMinusDeltaE = sigmaMinus;
242 sigmaPlusDeltaE = sigmaPlus;
243 sigmaDeltaE = std::sqrt(0.5 * sigmaMinusDeltaE * sigmaMinusDeltaE +
244 0.5 * sigmaPlusDeltaE * sigmaPlusDeltaE);
247 double sigmaPlusTot =
248 std::sqrt(sigmaPlus * sigmaPlus + caloEnergyError * caloEnergyError);
249 if (m_optimalRadiation) {
251 std::sqrt(4.65 * sigmaDeltaE_ioni * 4.65 * sigmaDeltaE_ioni +
252 caloEnergyError * caloEnergyError);
254 double MOPtot = std::abs(MOP + MOPshift);
255 if (m_optimalRadiation) {
256 MOPtot = std::abs(MOP + MOPshiftNoRad);
259 if (caloEnergy > MOPtot + 2 * sigmaPlusTot) {
266 double MOPreso = isign * 3.59524 * sigmaL *
correction;
267 deltaE = isign * caloEnergy + MOPreso;
268 sigmaMinusDeltaE = caloEnergyError + 0.08 * sigmaDeltaE_rad;
269 sigmaPlusDeltaE = caloEnergyError + 1.16 * sigmaDeltaE_rad;
270 sigmaDeltaE = std::sqrt(0.5 * sigmaMinusDeltaE * sigmaMinusDeltaE +
271 0.5 * sigmaPlusDeltaE * sigmaPlusDeltaE);
279 sigmaL = sigmaDeltaE_ioni + 0.3 * fracErad / 3.59524;
280 xc = momentumError / (sigmaL > 0. ? sigmaL : 1);
282 (0.3849 * xc * xc + 7.76672e-03 * xc * xc * xc) /
283 (1 + 2.8631 * xc + 0.3849 * xc * xc + 7.76672
e-03 * xc * xc * xc);
284 double MOPreso = isign * 3.59524 * sigmaL *
correction;
289 deltaE = MOP + MOPshift + MOPreso;
290 sigmaMinusDeltaE = sigmaMinus;
291 sigmaPlusDeltaE = sigmaPlus;
292 sigmaDeltaE = std::sqrt(0.5 * sigmaMinusDeltaE * sigmaMinusDeltaE +
293 0.5 * sigmaPlusDeltaE * sigmaPlusDeltaE);
314 double& ElossScale)
const
326 double X0CaloGirder[4] = {
327 -1.02877e-01, -2.74322e-02, 8.12989e-02, 9.73551e-01
331 constexpr
double X0CaloScale[60] = {
332 1.01685, 1.02092, 1.01875, 1.01812, 1.01791, 1.01345, 1.01354,
333 1.02145, 1.01645, 1.01585, 1.0172, 1.02262, 1.01464, 0.990931,
334 0.971953, 0.99845, 1.01433, 0.982143, 0.974015, 0.978742, 0.960029,
335 0.966766, 0.980199, 0.989586, 0.997144, 1.00169, 0.994166, 0.966332,
336 0.93671, 0.935656, 0.921994, 0.901489, 0.897799, 0.89638, 0.905629,
337 0.903374, 0.925922, 0.941203, 0.956273, 0.968618, 0.976883, 0.988349,
338 0.99855, 1.00212, 1.01456, 1.01541, 1.02532, 1.03238, 1.03688,
339 1.03783, 1.02078, 1.01529, 1.0156, 1.02212, 1.02226, 1.02406,
340 1.02188, 1.00661, 1.00661, 1.00661
344 constexpr
double ElossCaloScale[30] = {
345 1.06921, 1.06828, 1.06734, 1.06092, 1.06638, 1.06335, 1.07421, 1.05885,
346 1.07351, 1.07435, 1.06902, 1.07704, 1.08782, 1.09844, 1.115, 1.07609,
347 1.08233, 1.08764, 1.08209, 1.08255, 1.08008, 1.07573, 1.077, 1.07271,
348 1.07343, 1.07769, 1.07794, 1.08377, 1.08377, 1.08377
352 constexpr
double X0MuonScale[60] = {
353 -0.0320612, -0.0320612, -0.0320612, -0.0320612, -0.0693796, -0.0389677,
354 -0.0860891, -0.124606, -0.0882329, -0.100014, -0.0790912, -0.0745538,
355 -0.099088, -0.0933711, -0.0618782, -0.0619762, -0.0658361, -0.109704,
356 -0.129547, -0.143364, -0.0774768, -0.0739859, -0.0417835, -0.022119,
357 0.00308797, 0.0197657, -0.0137871, -0.036848, -0.0643794, -0.0514949,
358 -0.0317105, 0.016539, 0.0308435, -0.00056883, -0.00756813, -0.00760612,
359 -0.0234571, -0.0980915, -0.101175, -0.102354, -0.0920337, -0.100337,
360 -0.0887628, 0.0660931, 0.228999, 0.260675, 0.266301, 0.267907,
361 0.281668, 0.194433, 0.132954, 0.20707, 0.220466, 0.20936,
362 0.191441, 0.191441, 0.191441, 0.191441, 0.191441, 0.191441
365 int i60 = std::abs(eta) * 20.;
379 phi + 3.1416 - 3.1416 / 32. *
int((3.1416 +
phi) / (3.1416 / 32.));
381 if (
x >
M_PI / 64.) {
387 X0CaloGirder[0] * (1 -
x / 0.005) + X0CaloGirder[1] + X0CaloGirder[3];
388 }
else if (
x < 0.017) {
389 scale = X0CaloGirder[1] + X0CaloGirder[3];
390 }
else if (
x < 0.028) {
391 scale = X0CaloGirder[2] + X0CaloGirder[3];
393 scale = X0CaloGirder[3];
396 if (std::abs(eta) > 1.3) {
402 scale *= X0CaloScale[i60];
407 int i30 = std::abs(eta) * 10.;
415 double nfactor = 0.987363 / 1.05471;
417 ElossScale = nfactor * ElossCaloScale[i30] *
scale;
424 double scale = 1. + X0MuonScale[i60];
431 ElossScale = 0.93 *
scale;
438 double pathcorrection,
444 double pathLength = pathcorrection *
mat.thicknessInX0() *
mat.x0();
453 return (!m_detailedEloss