14 #include "GaudiKernel/MsgStream.h"
15 #include "GaudiKernel/SystemOfUnits.h"
26 , m_detailedEloss(true)
27 , m_optimalRadiation(true)
29 declareInterface<Trk::IEnergyLossUpdator>(
this);
47 double E = std::sqrt(
p *
p +
m *
m);
60 dEdX += -0.5345 /
mat.x0() + 6.803e-5 *
E /
mat.x0() +
61 2.278e-11 *
E *
E /
mat.x0() -
62 9.899e-18 *
E *
E *
E /
mat.x0();
64 dEdX += -2.986 /
mat.x0() + 9.253e-5 *
E /
mat.x0();
74 double pathcorrection,
81 "undefined ParticleHypothesis, energy loss calculation cancelled");
93 double pathLength = pathcorrection *
mat.thicknessInX0() *
mat.x0();
103 meanIoni =
sign * pathLength * meanIoni;
104 meanRad =
sign * pathLength * meanRad;
105 sigIoni = pathLength * sigIoni;
106 sigRad = pathLength * sigRad;
114 deltaE = meanIoni + meanRad;
116 double sigmaDeltaE = std::sqrt(sigIoni * sigIoni + sigRad * sigRad);
118 << deltaE <<
" meanIoni " << meanIoni <<
" meanRad " << meanRad
119 <<
" sigIoni " << sigIoni <<
" sigRad " << sigRad <<
" sign "
120 <<
sign <<
" pathLength " << pathLength);
123 sigmaDeltaE, meanIoni, sigIoni,
124 meanRad, sigRad, pathLength));
132 double caloEnergyError,
134 double momentumError,
135 int& elossFlag)
const
166 double deltaE = eLoss.
deltaE();
169 double deltaE_ioni = eLoss.
meanIoni();
170 double sigmaDeltaE_ioni = 0.45 * eLoss.
sigmaIoni();
171 double deltaE_rad = eLoss.
meanRad();
172 double sigmaDeltaE_rad =
178 if (eLoss.
meanRad() > 100000.) {
179 deltaE_rad = 100000.;
186 double MOP = deltaE_ioni - isign * 3.59524 * sigmaDeltaE_ioni;
192 isign * 50 * 10000. / pCaloEntry +
193 isign * 0.75 * std::sqrt(sigmaDeltaE_ioni * sigmaDeltaE_rad);
194 double MOPshiftNoRad = isign * 50 * 10000. / pCaloEntry;
198 double fracErad = sigmaDeltaE_rad +
199 std::abs(deltaE_rad) * pCaloEntry / (800000. + pCaloEntry);
200 double sigmaL = sigmaDeltaE_ioni + 0.8 * fracErad / 3.59524;
202 double sigmaMinus = 1.02 * sigmaDeltaE_ioni + 0.08 * sigmaDeltaE_rad;
203 double sigmaPlus = 4.65 * sigmaDeltaE_ioni + 1.16 * sigmaDeltaE_rad;
206 double xc = momentumError / (sigmaL > 0. ? sigmaL : 1.);
208 (0.3849 * xc * xc + 7.76672e-03 * xc * xc * xc) /
209 (1 + 2.8631 * xc + 0.3849 * xc * xc + 7.76672
e-03 * xc * xc * xc);
216 if (caloEnergyError <= 0) {
220 double MOPreso = isign * 3.59524 * sigmaL *
correction;
222 deltaE = MOP + MOPshift + MOPreso;
223 sigmaMinusDeltaE = sigmaMinus;
224 sigmaPlusDeltaE = sigmaPlus;
225 sigmaDeltaE = std::sqrt(0.5 * sigmaMinusDeltaE * sigmaMinusDeltaE +
226 0.5 * sigmaPlusDeltaE * sigmaPlusDeltaE);
228 if (m_optimalRadiation && std::abs(deltaE) < caloEnergy &&
229 pCaloEntry > 100000) {
238 sigmaL = sigmaDeltaE_ioni + 0.3 * fracErad / 3.59524;
239 xc = momentumError / (sigmaL > 0. ? sigmaL : 1.);
241 (0.3849 * xc * xc + 7.76672e-03 * xc * xc * xc) /
242 (1 + 2.8631 * xc + 0.3849 * xc * xc + 7.76672
e-03 * xc * xc * xc);
244 MOPreso = isign * 3.59524 * sigmaL *
correction;
245 deltaE = MOP + MOPshift + MOPreso;
246 sigmaMinusDeltaE = sigmaMinus;
247 sigmaPlusDeltaE = sigmaPlus;
248 sigmaDeltaE = std::sqrt(0.5 * sigmaMinusDeltaE * sigmaMinusDeltaE +
249 0.5 * sigmaPlusDeltaE * sigmaPlusDeltaE);
252 double sigmaPlusTot =
253 std::sqrt(sigmaPlus * sigmaPlus + caloEnergyError * caloEnergyError);
254 if (m_optimalRadiation) {
256 std::sqrt(4.65 * sigmaDeltaE_ioni * 4.65 * sigmaDeltaE_ioni +
257 caloEnergyError * caloEnergyError);
259 double MOPtot = std::abs(MOP + MOPshift);
260 if (m_optimalRadiation) {
261 MOPtot = std::abs(MOP + MOPshiftNoRad);
264 if (caloEnergy > MOPtot + 2 * sigmaPlusTot) {
271 double MOPreso = isign * 3.59524 * sigmaL *
correction;
272 deltaE = isign * caloEnergy + MOPreso;
273 sigmaMinusDeltaE = caloEnergyError + 0.08 * sigmaDeltaE_rad;
274 sigmaPlusDeltaE = caloEnergyError + 1.16 * sigmaDeltaE_rad;
275 sigmaDeltaE = std::sqrt(0.5 * sigmaMinusDeltaE * sigmaMinusDeltaE +
276 0.5 * sigmaPlusDeltaE * sigmaPlusDeltaE);
284 sigmaL = sigmaDeltaE_ioni + 0.3 * fracErad / 3.59524;
285 xc = momentumError / (sigmaL > 0. ? sigmaL : 1);
287 (0.3849 * xc * xc + 7.76672e-03 * xc * xc * xc) /
288 (1 + 2.8631 * xc + 0.3849 * xc * xc + 7.76672
e-03 * xc * xc * xc);
289 double MOPreso = isign * 3.59524 * sigmaL *
correction;
294 deltaE = MOP + MOPshift + MOPreso;
295 sigmaMinusDeltaE = sigmaMinus;
296 sigmaPlusDeltaE = sigmaPlus;
297 sigmaDeltaE = std::sqrt(0.5 * sigmaMinusDeltaE * sigmaMinusDeltaE +
298 0.5 * sigmaPlusDeltaE * sigmaPlusDeltaE);
319 double& ElossScale)
const
331 double X0CaloGirder[4] = {
332 -1.02877e-01, -2.74322e-02, 8.12989e-02, 9.73551e-01
336 constexpr
double X0CaloScale[60] = {
337 1.01685, 1.02092, 1.01875, 1.01812, 1.01791, 1.01345, 1.01354,
338 1.02145, 1.01645, 1.01585, 1.0172, 1.02262, 1.01464, 0.990931,
339 0.971953, 0.99845, 1.01433, 0.982143, 0.974015, 0.978742, 0.960029,
340 0.966766, 0.980199, 0.989586, 0.997144, 1.00169, 0.994166, 0.966332,
341 0.93671, 0.935656, 0.921994, 0.901489, 0.897799, 0.89638, 0.905629,
342 0.903374, 0.925922, 0.941203, 0.956273, 0.968618, 0.976883, 0.988349,
343 0.99855, 1.00212, 1.01456, 1.01541, 1.02532, 1.03238, 1.03688,
344 1.03783, 1.02078, 1.01529, 1.0156, 1.02212, 1.02226, 1.02406,
345 1.02188, 1.00661, 1.00661, 1.00661
349 constexpr
double ElossCaloScale[30] = {
350 1.06921, 1.06828, 1.06734, 1.06092, 1.06638, 1.06335, 1.07421, 1.05885,
351 1.07351, 1.07435, 1.06902, 1.07704, 1.08782, 1.09844, 1.115, 1.07609,
352 1.08233, 1.08764, 1.08209, 1.08255, 1.08008, 1.07573, 1.077, 1.07271,
353 1.07343, 1.07769, 1.07794, 1.08377, 1.08377, 1.08377
357 constexpr
double X0MuonScale[60] = {
358 -0.0320612, -0.0320612, -0.0320612, -0.0320612, -0.0693796, -0.0389677,
359 -0.0860891, -0.124606, -0.0882329, -0.100014, -0.0790912, -0.0745538,
360 -0.099088, -0.0933711, -0.0618782, -0.0619762, -0.0658361, -0.109704,
361 -0.129547, -0.143364, -0.0774768, -0.0739859, -0.0417835, -0.022119,
362 0.00308797, 0.0197657, -0.0137871, -0.036848, -0.0643794, -0.0514949,
363 -0.0317105, 0.016539, 0.0308435, -0.00056883, -0.00756813, -0.00760612,
364 -0.0234571, -0.0980915, -0.101175, -0.102354, -0.0920337, -0.100337,
365 -0.0887628, 0.0660931, 0.228999, 0.260675, 0.266301, 0.267907,
366 0.281668, 0.194433, 0.132954, 0.20707, 0.220466, 0.20936,
367 0.191441, 0.191441, 0.191441, 0.191441, 0.191441, 0.191441
370 int i60 = std::abs(eta) * 20.;
384 phi + 3.1416 - 3.1416 / 32. *
int((3.1416 +
phi) / (3.1416 / 32.));
386 if (
x >
M_PI / 64.) {
392 X0CaloGirder[0] * (1 -
x / 0.005) + X0CaloGirder[1] + X0CaloGirder[3];
393 }
else if (
x < 0.017) {
394 scale = X0CaloGirder[1] + X0CaloGirder[3];
395 }
else if (
x < 0.028) {
396 scale = X0CaloGirder[2] + X0CaloGirder[3];
398 scale = X0CaloGirder[3];
401 if (std::abs(eta) > 1.3) {
407 scale *= X0CaloScale[i60];
412 int i30 = std::abs(eta) * 10.;
420 double nfactor = 0.987363 / 1.05471;
422 ElossScale = nfactor * ElossCaloScale[i30] *
scale;
429 double scale = 1. + X0MuonScale[i60];
436 ElossScale = 0.93 *
scale;
443 double pathcorrection,
449 double pathLength = pathcorrection *
mat.thicknessInX0() *
mat.x0();
458 return (!m_detailedEloss