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;
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");