30 #include <TDatabasePDG.h>
44 energyLossBetheBloch_(true), noiseBetheBloch_(true),
46 energyLossBrems_(true), noiseBrems_(true),
47 ignoreBoundariesBetweenEqualMaterials_(true),
61 materialInterface_(nullptr),
88 std::string msg(
"MaterialEffects::initMaterialInterface(): Already initialized! ");
89 std::runtime_error err(msg);
98 if (modelName == std::string(
"GEANE")) {
100 }
else if (modelName == std::string(
"Highland")) {
103 std::string errorMsg = std::string(
"There is no MSC model called \"") + modelName +
"\". Maybe it is not implemented or you misspelled the model name";
104 Exception exc(errorMsg, __LINE__, __FILE__);
113 int materialsFXStart,
121 debugOut <<
" MaterialEffects::effects \n";
135 std::string msg(
"MaterialEffects hasn't been initialized with a correct AbsMaterialInterface pointer!");
136 std::runtime_error err(msg);
140 bool doNoise(noise !=
nullptr);
147 for ( std::vector<RKStep>::const_iterator it = steps.begin() + materialsFXStart; it != steps.begin() + materialsFXStop; ++it) {
149 double realPath = it->matStep_.stepSize_;
150 if (fabs(realPath) < 1.E-8) {
160 debugOut <<
"stepSize = " << it->matStep_.stepSize_ <<
"\t";
161 it->matStep_.materialProperties_.Print();
167 realPath = fabs(realPath);
174 momLoss +=
momentumLoss(stepSign, mom - momLoss,
false);
178 double p(0), gammaSquare(0), gamma(0), betaSquare(0);
180 double pSquare = p*p;
186 this->
noiseCoulomb(*noise, *((
M1x3*) &it->state7_[3]), pSquare, betaSquare);
189 this->
noiseBrems(*noise, pSquare, betaSquare);
196 if (momLoss >= mom) {
197 Exception exc(
"MaterialEffects::effects ==> momLoss >= momentum, aborting extrapolation!",__LINE__,__FILE__);
216 static const double maxRelMomLoss = .01;
217 static const double Pmin = 4.E-3;
218 static const double minStep = 1.E-4;
222 std::ostringstream sstream;
223 sstream <<
"MaterialEffects::stepper ==> momentum too low: " << mom*1000. <<
" MeV";
224 Exception exc(sstream.str(),__LINE__,__FILE__);
234 std::string msg(
"MaterialEffects hasn't been initialized with a correct AbsMaterialInterface pointer!");
235 std::runtime_error err(msg);
239 if (relMomLoss > maxRelMomLoss) {
247 if (fabs(sMax) < minStep)
257 state7[0] += limits.
getStepSign() * minStep * state7[3];
258 state7[1] += limits.
getStepSign() * minStep * state7[4];
259 state7[2] += limits.
getStepSign() * minStep * state7[5];
273 double relMomLossPer_cm(0);
280 double maxStepMomLoss = fabs((maxRelMomLoss - fabs(relMomLoss)) / relMomLossPer_cm);
284 debugOut <<
" momLoss exceeded after a step of " << maxStepMomLoss
285 <<
"; relMomLoss up to now = " << relMomLoss <<
"\n";
295 double boundaryStep(sMax);
297 for (
unsigned int i=0; i<100; ++i) {
299 debugOut <<
" find next boundary\n";
305 debugOut <<
" materialInterface_ returned a step of 0 \n";
310 boundaryStep -= step;
313 debugOut <<
" made a step of " << step <<
"\n";
323 rep->
RKPropagate(state7, NULL, SA, step, varField);
326 state7[0] += limits.
getStepSign() * minStep * state7[3];
327 state7[1] += limits.
getStepSign() * minStep * state7[4];
328 state7[2] += limits.
getStepSign() * minStep * state7[5];
339 if (materialAfter != currentMaterial)
352 TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(
pdg_);
353 charge_ = int(part->Charge() / 3.);
354 mass_ = part->Mass();
359 double& mom,
double& gammaSquare,
double& gamma,
double& betaSquare)
const {
361 if (Energy <=
mass_) {
362 Exception exc(
"MaterialEffects::getMomGammaBeta - Energy <= mass",__LINE__,__FILE__);
366 gamma = Energy/
mass_;
367 gammaSquare = gamma*gamma;
368 betaSquare = 1.-1./gammaSquare;
369 mom = Energy*sqrt(betaSquare);
378 double E0 = hypot(mom,
mass_);
395 double dEdx1 =
dEdx(E0);
401 double E1 = E0 - dEdx1*step/2.;
402 double dEdx2 =
dEdx(E1);
404 double E2 = E0 - dEdx2*step/2.;
405 double dEdx3 =
dEdx(E2);
407 double E3 = E0 - dEdx3*step;
408 double dEdx4 =
dEdx(E3);
410 dEdx_ = (dEdx1 + 2.*dEdx2 + 2.*dEdx3 + dEdx4)/6.;
415 double dE = step*
dEdx_;
419 if (E0 - dE <=
mass_) {
421 return momLoss = mom;
423 else momLoss = mom - sqrt(pow(E0 - dE, 2) -
mass_*
mass_);
426 debugOut <<
" MaterialEffects::momentumLoss: mom = " << mom <<
"; E0 = " << E0
427 <<
"; dEdx = " << dEdx_
428 <<
"; dE = " << dE <<
"; mass = " << mass_ <<
"\n";
439 double mom(0), gammaSquare(0), gamma(0), betaSquare(0);
456 static const double betaGammaMin(0.05);
457 if (betaSquare*gammaSquare < betaGammaMin*betaGammaMin) {
458 Exception exc(
"MaterialEffects::dEdxBetheBloch ==> beta*gamma < 0.05, Bethe-Bloch implementation not valid anymore!",__LINE__,__FILE__);
466 double argument( gammaSquare * betaSquare *
me_ * 1.E3 * 2. / ((1.E-6 *
mEE_) *
467 sqrt(1. + 2. * gamma * massRatio + massRatio * massRatio)) );
468 result *= log(argument) - betaSquare;
483 double sigma2E ( 0. );
486 double kappa ( zeta / Emax );
489 sigma2E += zeta * Emax * (1. - betaSquare / 2.);
492 double I = 16. * pow(
matZ_, 0.9);
497 double e1 = pow((I / pow(e2, f2)), 1. / f1);
499 double mbbgg2 = 2.E9 *
mass_ * betaSquare * gammaSquare;
500 double Sigma1 =
dEdx_ * 1.0E9 * f1 / e1 * (log(mbbgg2 / e1) - betaSquare) / (log(mbbgg2 / I) - betaSquare) * 0.6;
501 double Sigma2 =
dEdx_ * 1.0E9 * f2 / e2 * (log(mbbgg2 / e2) - betaSquare) / (log(mbbgg2 / I) - betaSquare) * 0.6;
502 double Sigma3 =
dEdx_ * 1.0E9 * Emax / (I * (Emax + I) * log((Emax + I) / I)) * 0.4;
504 double Nc = (Sigma1 + Sigma2 + Sigma3) * fabs(
stepSize_);
507 double sigmaalpha = 15.76;
509 double RLAMED = -0.422784 - betaSquare - log(zeta / Emax);
510 double RLAMAX = 0.60715 + 1.1934 * RLAMED + (0.67794 + 0.052382 * RLAMED) * exp(0.94753 + 0.74442 * RLAMED);
512 if (RLAMAX <= 1010.) {
513 sigmaalpha = 1.975560
514 + 9.898841e-02 * RLAMAX
515 - 2.828670e-04 * RLAMAX * RLAMAX
516 + 5.345406e-07 * pow(RLAMAX, 3.)
517 - 4.942035e-10 * pow(RLAMAX, 4.)
518 + 1.729807e-13 * pow(RLAMAX, 5.);
519 }
else { sigmaalpha = 1.871887E+01 + 1.296254E-02 * RLAMAX; }
521 if (sigmaalpha > 54.6) sigmaalpha = 54.6;
522 sigma2E += sigmaalpha * sigmaalpha * zeta * zeta;
524 static const double alpha = 0.996;
525 double Ealpha = I / (1. - (alpha * Emax / (Emax + I)));
526 double meanE32 = I * (Emax + I) / Emax * (Ealpha - I);
527 sigma2E += fabs(
stepSize_) * (Sigma1 * e1 * e1 + Sigma2 * e2 * e2 + Sigma3 * meanE32);
534 noise[6 * 7 + 6] +=
charge_*
charge_/betaSquare / pow(mom, 4) * sigma2E;
539 const M1x3& direction,
double momSquare,
double betaSquare)
const
546 const double step2 = step * step;
552 double logCor = (1 + 0.038 * log(stepOverRadLength));
553 sigma2 = 0.0136 * 0.0136 *
charge_ *
charge_ / (betaSquare * momSquare) * stepOverRadLength * logCor * logCor;
556 sigma2 = (sigma2 > 0.0 ? sigma2 : 0.0);
560 std::fill(noiseAfter.
begin(), noiseAfter.
end(), 0);
562 const M1x3& a = direction;
565 noiseAfter[0 * 7 + 0] = sigma2 * step2 / 3.0 * (1 - a[0]*a[0]);
566 noiseAfter[1 * 7 + 0] = -sigma2 * step2 / 3.0 * a[0]*a[1];
567 noiseAfter[2 * 7 + 0] = -sigma2 * step2 / 3.0 * a[0]*a[2];
568 noiseAfter[3 * 7 + 0] = sigma2 * step * 0.5 * (1 - a[0]*a[0]);
569 noiseAfter[4 * 7 + 0] = -sigma2 * step * 0.5 * a[0]*a[1];
570 noiseAfter[5 * 7 + 0] = -sigma2 * step * 0.5 * a[0]*a[1];
571 noiseAfter[0 * 7 + 1] = noiseAfter[1 * 7 + 0];
572 noiseAfter[1 * 7 + 1] = sigma2 * step2 / 3.0 * (1 - a[1]*a[1]);
573 noiseAfter[2 * 7 + 1] = -sigma2 * step2 / 3.0 * a[1]*a[2];
574 noiseAfter[3 * 7 + 1] = noiseAfter[4 * 7 + 0];
575 noiseAfter[4 * 7 + 1] = sigma2 * step * 0.5 * (1 - a[1] * a[1]);
576 noiseAfter[5 * 7 + 1] = -sigma2 * step * 0.5 * a[1]*a[2];
577 noiseAfter[0 * 7 + 2] = noiseAfter[2 * 7 + 0];
578 noiseAfter[1 * 7 + 2] = noiseAfter[2 * 7 + 1];
579 noiseAfter[2 * 7 + 2] = sigma2 * step2 / 3.0 * (1 - a[2]*a[2]);
580 noiseAfter[3 * 7 + 2] = noiseAfter[5 * 7 + 0];
581 noiseAfter[4 * 7 + 2] = noiseAfter[5 * 7 + 1];
582 noiseAfter[5 * 7 + 2] = sigma2 * step * 0.5 * (1 - a[2]*a[2]);
583 noiseAfter[0 * 7 + 3] = noiseAfter[3 * 7 + 0];
584 noiseAfter[1 * 7 + 3] = noiseAfter[3 * 7 + 1];
585 noiseAfter[2 * 7 + 3] = noiseAfter[3 * 7 + 2];
586 noiseAfter[3 * 7 + 3] = sigma2 * (1 - a[0]*a[0]);
587 noiseAfter[4 * 7 + 3] = -sigma2 * a[0]*a[1];
588 noiseAfter[5 * 7 + 3] = -sigma2 * a[0]*a[2];
589 noiseAfter[0 * 7 + 4] = noiseAfter[4 * 7 + 0];
590 noiseAfter[1 * 7 + 4] = noiseAfter[4 * 7 + 1];
591 noiseAfter[2 * 7 + 4] = noiseAfter[4 * 7 + 2];
592 noiseAfter[3 * 7 + 4] = noiseAfter[4 * 7 + 3];
593 noiseAfter[4 * 7 + 4] = sigma2 * (1 - a[1]*a[1]);
594 noiseAfter[5 * 7 + 4] = -sigma2 * a[1]*a[2];
595 noiseAfter[0 * 7 + 5] = noiseAfter[5 * 7 + 0];
596 noiseAfter[1 * 7 + 5] = noiseAfter[5 * 7 + 1];
597 noiseAfter[2 * 7 + 5] = noiseAfter[5 * 7 + 2];
598 noiseAfter[3 * 7 + 5] = noiseAfter[5 * 7 + 3];
599 noiseAfter[4 * 7 + 5] = noiseAfter[5 * 7 + 4];
600 noiseAfter[5 * 7 + 5] = sigma2 * (1 - a[2]*a[2]);
603 for (
unsigned int i = 0; i < 7 * 7; ++i) {
604 noise[i] += noiseAfter[i];
614 if (abs(
pdg_) != 11)
return 0;
617 static const double C[101] = { 0.0, -0.960613E-01, 0.631029E-01, -0.142819E-01, 0.150437E-02, -0.733286E-04, 0.131404E-05, 0.859343E-01, -0.529023E-01, 0.131899E-01, -0.159201E-02, 0.926958E-04, -0.208439E-05, -0.684096E+01, 0.370364E+01, -0.786752E+00, 0.822670E-01, -0.424710E-02, 0.867980E-04, -0.200856E+01, 0.129573E+01, -0.306533E+00, 0.343682E-01, -0.185931E-02, 0.392432E-04, 0.127538E+01, -0.515705E+00, 0.820644E-01, -0.641997E-02, 0.245913E-03, -0.365789E-05, 0.115792E+00, -0.463143E-01, 0.725442E-02, -0.556266E-03, 0.208049E-04, -0.300895E-06, -0.271082E-01, 0.173949E-01, -0.452531E-02, 0.569405E-03, -0.344856E-04, 0.803964E-06, 0.419855E-02, -0.277188E-02, 0.737658E-03, -0.939463E-04, 0.569748E-05, -0.131737E-06, -0.318752E-03, 0.215144E-03, -0.579787E-04, 0.737972E-05, -0.441485E-06, 0.994726E-08, 0.938233E-05, -0.651642E-05, 0.177303E-05, -0.224680E-06, 0.132080E-07, -0.288593E-09, -0.245667E-03, 0.833406E-04, -0.129217E-04, 0.915099E-06, -0.247179E-07, 0.147696E-03, -0.498793E-04, 0.402375E-05, 0.989281E-07, -0.133378E-07, -0.737702E-02, 0.333057E-02, -0.553141E-03, 0.402464E-04, -0.107977E-05, -0.641533E-02, 0.290113E-02, -0.477641E-03, 0.342008E-04, -0.900582E-06, 0.574303E-05, 0.908521E-04, -0.256900E-04, 0.239921E-05, -0.741271E-07, -0.341260E-04, 0.971711E-05, -0.172031E-06, -0.119455E-06, 0.704166E-08, 0.341740E-05, -0.775867E-06, -0.653231E-07, 0.225605E-07, -0.114860E-08, -0.119391E-06, 0.194885E-07, 0.588959E-08, -0.127589E-08, 0.608247E-10};
618 static const double xi = 2.51, beta = 0.99, vl = 0.00004;
620 #if defined(BETHE) // no MIGDAL corrections
621 static const double C[101] = { 0.0, 0.834459E-02, 0.443979E-02, -0.101420E-02, 0.963240E-04, -0.409769E-05, 0.642589E-07, 0.464473E-02, -0.290378E-02, 0.547457E-03, -0.426949E-04, 0.137760E-05, -0.131050E-07, -0.547866E-02, 0.156218E-02, -0.167352E-03, 0.101026E-04, -0.427518E-06, 0.949555E-08, -0.406862E-02, 0.208317E-02, -0.374766E-03, 0.317610E-04, -0.130533E-05, 0.211051E-07, 0.158941E-02, -0.385362E-03, 0.315564E-04, -0.734968E-06, -0.230387E-07, 0.971174E-09, 0.467219E-03, -0.154047E-03, 0.202400E-04, -0.132438E-05, 0.431474E-07, -0.559750E-09, -0.220958E-02, 0.100698E-02, -0.596464E-04, -0.124653E-04, 0.142999E-05, -0.394378E-07, 0.477447E-03, -0.184952E-03, -0.152614E-04, 0.848418E-05, -0.736136E-06, 0.190192E-07, -0.552930E-04, 0.209858E-04, 0.290001E-05, -0.133254E-05, 0.116971E-06, -0.309716E-08, 0.212117E-05, -0.103884E-05, -0.110912E-06, 0.655143E-07, -0.613013E-08, 0.169207E-09, 0.301125E-04, -0.461920E-04, 0.871485E-05, -0.622331E-06, 0.151800E-07, -0.478023E-04, 0.247530E-04, -0.381763E-05, 0.232819E-06, -0.494487E-08, -0.336230E-04, 0.223822E-04, -0.384583E-05, 0.252867E-06, -0.572599E-08, 0.105335E-04, -0.567074E-06, -0.216564E-06, 0.237268E-07, -0.658131E-09, 0.282025E-05, -0.671965E-06, 0.565858E-07, -0.193843E-08, 0.211839E-10, 0.157544E-04, -0.304104E-05, -0.624410E-06, 0.120124E-06, -0.457445E-08, -0.188222E-05, -0.407118E-06, 0.375106E-06, -0.466881E-07, 0.158312E-08, 0.945037E-07, 0.564718E-07, -0.319231E-07, 0.371926E-08, -0.123111E-09};
622 static const double xi = 2.10, beta = 1.00, vl = 0.001;
625 double BCUT = 10000.;
627 static const double THIGH = 100., CHIGH = 50.;
628 double dedxBrems = 0.;
633 if (BCUT >= mom) BCUT = mom;
639 if (BCUT >= THIGH) kc = CHIGH;
647 if (BCUT > T) kc = T;
649 double X = log(T / me_);
650 double Y = log(kc / (E * vl));
654 double S = 0., YY = 1.;
656 for (
unsigned int I = 1; I <= 2; ++I) {
658 for (
unsigned int J = 1; J <= 6; ++J) {
660 S = S + C[K] * XX * YY;
666 for (
unsigned int I = 3; I <= 6; ++I) {
668 for (
unsigned int J = 1; J <= 6; ++J) {
670 if (Y <= 0.) S = S + C[K] * XX * YY;
671 else S = S + C[K + 24] * XX * YY;
680 for (
unsigned int I = 1; I <= 2; ++I) {
682 for (
unsigned int J = 1; J <= 5; ++J) {
684 SS = SS + C[K] * XX * YY;
690 for (
unsigned int I = 3; I <= 5; ++I) {
692 for (
unsigned int J = 1; J <= 5; ++J) {
694 if (Y <= 0.) SS = SS + C[K] * XX * YY;
695 else SS = SS + C[K + 15] * XX * YY;
712 double FAC =
matZ_ * (
matZ_ + xi) * E * E / (E + me_);
714 FAC *= kc * CORR / T;
716 FAC *= exp(beta * log(kc * CORR / T));
718 if (FAC <= 0.)
return 0.;
726 S = (1. - 0.5 * RAT + 2.*RAT * RAT / 9.);
728 S = S / (1. - 0.5 * RAT + 2.*RAT * RAT / 9.);
731 S = BCUT * (1. - 0.5 * RAT + 2.*RAT * RAT / 9.);
733 S = S / (kc * (1. - 0.5 * RAT + 2.*RAT * RAT / 9.));
735 dedxBrems = dedxBrems * S;
742 if (dedxBrems < 0.) dedxBrems = 0;
747 static const double AA = 7522100., A1 = 0.415, A3 = 0.0021, A5 = 0.00054;
753 if (X >= +9.) ETA = 1.;
755 double W = A1 * X + A3 * pow(X, 3.) + A5 * pow(X, 5.);
756 ETA = 0.5 + atan(W) / M_PI;
761 if (ETA < 0.0001) factor = 1.E-10;
762 else if (ETA > 0.9999) factor = 1.;
764 double E0 = BCUT / mom;
765 if (E0 > 1.) E0 = 1.;
766 if (E0 < 1.E-8) factor = 1.;
767 else factor = ETA * (1. - pow(1. - E0, 1. / ETA)) / E0;
771 return factor * dedxBrems;
783 if (abs(
pdg_) != 11)
return;
786 double sigma2 = 1.44*(pow(3., minusXOverLn2) - pow(4., minusXOverLn2)) / momSquare;
789 sigma2 = (sigma2 > 0.0 ? sigma2 : 0.0);
790 noise[6 * 7 + 6] +=
charge_*
charge_/betaSquare / (momSquare*momSquare) * sigma2;
812 double minMom = 0.00001;
813 double maxMom = 10000;
815 double logStepSize = (log10(maxMom) - log10(minMom)) / (nSteps-1);
817 TH1D hdEdxBethe(
"dEdxBethe",
"dEdxBethe; log10(mom)", nSteps, log10(minMom), log10(maxMom));
818 TH1D hdEdxBrems(
"dEdxBrems",
"dEdxBrems; log10(mom)", nSteps, log10(minMom), log10(maxMom));
820 for (
int i=0; i<nSteps; ++i) {
821 double mom = pow(10., log10(minMom) + i*logStepSize);
822 double E = hypot(mom,
mass_);
828 hdEdxBethe.Fill(log10(mom),
dEdx(E));
840 hdEdxBrems.Fill(log10(mom),
dEdx(E));
851 std::stringstream convert;
853 Result = convert.str();
855 TFile outfile(
"dEdx_" + TString(Result) +
".root",
"recreate");
Helper to store different limits on the stepsize for the RKTRackRep.
void stepper(const RKTrackRep *rep, M1x7 &state7, const double &mom, double &relMomLoss, const int &pdg, MaterialProperties ¤tMaterial, StepLimits &limits, bool varField=true)
Returns maximum length so that a specified momentum loss will not be exceeded.
static MaterialEffects * instance_
double getLowestLimitSignedVal(double margin=1.E-3) const
Get the numerical value of the lowest limit, signed with stepSign_.
double RKPropagate(M1x7 &state7, M7x7 *jacobian, M1x3 &SA, double S, bool varField=true, bool calcOnlyLastRowOfJ=false) const
The actual Runge Kutta propagation.
double dEdxBrems(double mom) const
Returns dEdx.
bool energyLossBetheBloch_
void setDebugLvl(unsigned int lvl=1)
Abstract base class for geometry interfacing.
Stepper and energy loss/noise matrix calculation.
void getParticleParameters()
sets charge_, mass_
double momentumLoss(double stepSign, double mom, bool linear)
Returns momentum loss.
virtual double findNextBoundary(const RKTrackRep *rep, const M1x7 &state7, double sMax, bool varField=true)=0
Make a step until maxStep or the next boundary is reached.
void setLimit(StepLimitType type, double value)
absolute of value will be taken! If limit is already lower, it will be set to value anyway...
void noiseBrems(M7x7 &noise, double momSquare, double betaSquare) const
calculation of energy loss straggeling
double dEdx(double Energy) const
Calculate dEdx for a given energy.
virtual const char * what() const
Standard error message handling for exceptions. use like "std::cerr << e.what();".
void setFatal(bool b=true)
Set fatal flag.
void getMomGammaBeta(double Energy, double &mom, double &gammaSquare, double &gamma, double &betaSquare) const
virtual void getMaterialParameters(double &density, double &Z, double &A, double &radiationLength, double &mEE)=0
Get material parameters in current material.
AbsMaterialInterface * materialInterface_
depending on this number a specific msc model is chosen in the noiseCoulomb function.
void noiseCoulomb(M7x7 &noise, const M1x3 &direction, double momSquare, double betaSquare) const
calculation of multiple scattering
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
void setMscModel(const std::string &modelName)
Select the multiple scattering model that will be used during track fit.
bool ignoreBoundariesBetweenEqualMaterials_
void init(AbsMaterialInterface *matIfc)
set the material interface here. Material interface classes must be derived from AbsMaterialInterface...
AbsTrackRep with 5D track parameterization in plane coordinates: (q/p, u', v', u, v) ...
Material properties needed e.g. for material effects calculation.
virtual ~MaterialEffects()
void setMaterialProperties(const double &density, const double &Z, const double &A, const double &radiationLength, const double &mEE)
void drawdEdx(int pdg=11)
void Print(const Option_t *="") const
static MaterialEffects * getInstance()
virtual void setDebugLvl(unsigned int lvl=1)
virtual bool initTrack(double posX, double posY, double posZ, double dirX, double dirY, double dirZ)=0
Initialize the navigator at given position and with given direction. Return true if volume changed...
void noiseBetheBloch(M7x7 &noise, double mom, double betaSquare, double gamma, double gammaSquare) const
calculation of energy loss straggeling
double getLowestLimitVal(double margin=1.E-3) const
Get the unsigned numerical value of the lowest limit.
double dEdxBetheBloch(double betaSquare, double gamma, double gammasquare) const
Uses Bethe Bloch formula to calculate dEdx.
double effects(const std::vector< RKStep > &steps, int materialsFXStart, int materialsFXStop, const double &mom, const int &pdg, M7x7 *noise=nullptr)
Calculates energy loss in the traveled path, optional calculation of noise matrix.
Defines for I/O streams used for error and debug printing.