trackReps/src/MaterialEffects.cc

0001 /* Copyright 2008-2014, Technische Universitaet Muenchen,
0002    Authors: Christian Hoeppner & Sebastian Neubert
0003
0004    This file is part of GENFIT.
0005
0006    GENFIT is free software: you can redistribute it and/or modify
0007    it under the terms of the GNU Lesser General Public License as published
0008    by the Free Software Foundation, either version 3 of the License, or
0009    (at your option) any later version.
0010
0011    GENFIT is distributed in the hope that it will be useful,
0012    but WITHOUT ANY WARRANTY; without even the implied warranty of
0013    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0014    GNU Lesser General Public License for more details.
0015
0016    You should have received a copy of the GNU Lesser General Public License
0017    along with GENFIT.  If not, see <http://www.gnu.org/licenses/>.
0018 */
0019
0020 #include "MaterialEffects.h"
0021 #include "Exception.h"
0022 #include "IO.h"
0023
0024 #include <stdexcept>
0025 #include <string>
0026 #include <stdlib.h>
0027 #include <math.h>
0028 #include <assert.h>
0029
0030 #include <TDatabasePDG.h>
0031 #include "MonopoleConstants.h"
0032 #include <TMath.h>
0033
0034 #include <TH1D.h>
0035 #include <TFile.h>
0036
0037
0038 namespace genfit {
0039
0040 MaterialEffects* MaterialEffects::instance_ = nullptr;
0041
0042
0043 MaterialEffects::MaterialEffects():
0044   noEffects_(false),
0045   energyLossBetheBloch_(true), noiseBetheBloch_(true),
0046   noiseCoulomb_(true),
0047   energyLossBrems_(true), noiseBrems_(true),
0048   ignoreBoundariesBetweenEqualMaterials_(true),
0049   me_(0.510998910E-3),
0050   stepSize_(0),
0051   dEdx_(0),
0052   E_(0),
0053   matDensity_(0),
0054   matZ_(0),
0055   matA_(0),
0056   radiationLength_(0),
0057   mEE_(0),
0058   pdg_(0),
0059   charge_(0),
0060   mag_charge_(0),
0061   mass_(0),
0062   mscModelCode_(0),
0063   materialInterface_(nullptr),
0064   debugLvl_(0)
0065 {
0066 }
0067
0068 MaterialEffects::~MaterialEffects()
0069 {
0070   if (materialInterface_ != nullptr) delete materialInterface_;
0071 }
0072
0073 MaterialEffects* MaterialEffects::getInstance()
0074 {
0075   if (instance_ == nullptr) instance_ = new MaterialEffects();
0076   return instance_;
0077 }
0078
0079 void MaterialEffects::destruct()
0080 {
0081   if (instance_ != nullptr) {
0082     delete instance_;
0083     instance_ = nullptr;
0084   }
0085 }
0086
0087 void MaterialEffects::init(AbsMaterialInterface* matIfc)
0088 {
0089   if (materialInterface_ != nullptr) {
0090     std::string msg("MaterialEffects::initMaterialInterface(): Already initialized! ");
0091     std::runtime_error err(msg);
0092   }
0093   materialInterface_ = matIfc;
0094 }
0095
0096
0097
0098 void MaterialEffects::setMscModel(const std::string& modelName)
0099 {
0100   if (modelName == std::string("GEANE")) {
0101     mscModelCode_ = 0;
0102   } else if (modelName == std::string("Highland")) {
0103     mscModelCode_ = 1;
0104   } else {// throw exception
0105     std::string errorMsg = std::string("There is no MSC model called \"") + modelName + "\". Maybe it is not implemented or you misspelled the model name";
0106     Exception exc(errorMsg, __LINE__, __FILE__);
0107     exc.setFatal();
0108     errorOut << exc.what();
0109     throw exc;
0110   }
0111 }
0112
0113
0114 double MaterialEffects::effects(const std::vector<RKStep>& steps,
0115                                 int materialsFXStart,
0116                                 int materialsFXStop,
0117                                 const double& mom,
0118                                 const int& pdg,
0119                                 M7x7* noise)
0120 {
0121
0122   if (debugLvl_ > 0) {
0123     debugOut << "     MaterialEffects::effects \n";
0124   }
0125
0126   /*debugOut << "noEffects_ " << noEffects_ << "\n";
0127   debugOut << "energyLossBetheBloch_ " << energyLossBetheBloch_ << "\n";
0128   debugOut << "noiseBetheBloch_ " << noiseBetheBloch_ << "\n";
0129   debugOut << "noiseCoulomb_ " << noiseCoulomb_ << "\n";
0130   debugOut << "energyLossBrems_ " << energyLossBrems_ << "\n";
0131   debugOut << "noiseBrems_ " << noiseBrems_ << "\n";*/
0132
0133
0134   if (noEffects_) return 0.;
0135
0136   if (materialInterface_ == nullptr) {
0137     std::string msg("MaterialEffects hasn't been initialized with a correct AbsMaterialInterface pointer!");
0138     std::runtime_error err(msg);
0139     throw err;
0140   }
0141
0142   bool doNoise(noise != nullptr);
0143
0144   pdg_ = pdg;
0145   getParticleParameters();
0146
0147   double momLoss = 0.;
0148
0149   for ( std::vector<RKStep>::const_iterator it = steps.begin() + materialsFXStart; it !=  steps.begin() + materialsFXStop; ++it) { // loop over steps
0150
0151     double realPath = it->matStep_.stepSize_;
0152     if (fabs(realPath) < 1.E-8) {
0153       // do material effects only if distance is not too small
0154       continue;
0155     }
0156
0157     if (debugLvl_ > 0) {
0158       debugOut << "     calculate matFX ";
0159       if (doNoise)
0160         debugOut << "and noise";
0161       debugOut << " for ";
0162       debugOut << "stepSize = " << it->matStep_.stepSize_ << "\t";
0163       it->matStep_.material_.Print();
0164     }
0165
0166     double stepSign(1.);
0167     if (realPath < 0)
0168       stepSign = -1.;
0169     realPath = fabs(realPath);
0170     stepSize_ = realPath;
0171
0172
0173     const Material& currentMaterial = it->matStep_.material_;
0174     matDensity_ = currentMaterial.density;
0175     matZ_ = currentMaterial.Z;
0176     matA_ = currentMaterial.A;
0177     radiationLength_ = currentMaterial.radiationLength;
0178     mEE_ = currentMaterial.mEE;
0179
0180
0181     if (matZ_ > 1.E-3) { // don't calculate energy loss for vacuum
0182
0183       momLoss += momentumLoss(stepSign, mom - momLoss, false);
0184
0185       if (doNoise){
0186         // get values for the "effective" energy of the RK step E_
0187         double p(0), gammaSquare(0), gamma(0), betaSquare(0);
0188         this->getMomGammaBeta(E_, p, gammaSquare, gamma, betaSquare);
0189         double pSquare = p*p;
0190
0191         if (pdg_ == c_monopolePDGCode) {
0192           charge_ = mag_charge_ * mom / hypot(mom, mass_); //effective charge for monopoles
0193         }
0194
0195         if (energyLossBetheBloch_ && noiseBetheBloch_)
0196           this->noiseBetheBloch(*noise, p, betaSquare, gamma, gammaSquare);
0197
0198         if (noiseCoulomb_)
0199           this->noiseCoulomb(*noise, *((M1x3*) &it->state7_[3]), pSquare, betaSquare);
0200
0201         if (energyLossBrems_ && noiseBrems_)
0202           this->noiseBrems(*noise, pSquare, betaSquare);
0203       } // end doNoise
0204
0205     }
0206
0207   } // end loop over steps
0208
0209   if (momLoss >= mom) {
0210     Exception exc("MaterialEffects::effects ==> momLoss >= momentum, aborting extrapolation!",__LINE__,__FILE__);
0211     exc.setFatal();
0212     throw exc;
0213   }
0214
0215   return momLoss;
0216 }
0217
0218
0219 void MaterialEffects::stepper(const RKTrackRep* rep,
0220                               M1x7& state7,
0221                               const double& mom, // momentum
0222                               double& relMomLoss, // relative momloss for the step will be added
0223                               const int& pdg,
0224                               Material& currentMaterial,
0225                               StepLimits& limits,
0226                               bool varField)
0227 {
0228
0229   static const double maxRelMomLoss = .01; // maximum relative momentum loss allowed
0230   static const double Pmin   = 4.E-3;           // minimum momentum for propagation [GeV]
0231   static const double minStep = 1.E-4; // 1 µm
0232
0233   // check momentum
0234   if(mom < Pmin){
0235     std::ostringstream sstream;
0236     sstream << "MaterialEffects::stepper ==> momentum too low: " << mom*1000. << " MeV";
0237     Exception exc(sstream.str(),__LINE__,__FILE__);
0238     exc.setFatal();
0239     throw exc;
0240   }
0241
0242   // Trivial cases
0243   if (noEffects_)
0244     return;
0245
0246   if (materialInterface_ == nullptr) {
0247     std::string msg("MaterialEffects hasn't been initialized with a correct AbsMaterialInterface pointer!");
0248     std::runtime_error err(msg);
0249     throw err;
0250   }
0251
0252   if (relMomLoss > maxRelMomLoss) {
0253     limits.setLimit(stp_momLoss, 0);
0254     return;
0255   }
0256
0257
0258   double sMax = limits.getLowestLimitSignedVal(); // signed
0259
0260   if (fabs(sMax) < minStep)
0261     return;
0262
0263
0264
0265   pdg_ = pdg;
0266   getParticleParameters();
0267
0268
0269   // make minStep
0270   state7[0] += limits.getStepSign() * minStep * state7[3];
0271   state7[1] += limits.getStepSign() * minStep * state7[4];
0272   state7[2] += limits.getStepSign() * minStep * state7[5];
0273
0274   materialInterface_->initTrack(state7[0], state7[1], state7[2],
0275                                 limits.getStepSign() * state7[3], limits.getStepSign() * state7[4], limits.getStepSign() * state7[5]);
0276
0277
0278   currentMaterial = materialInterface_->getMaterialParameters();
0279   matDensity_ = currentMaterial.density;
0280   matZ_ = currentMaterial.Z;
0281   matA_ = currentMaterial.A;
0282   radiationLength_ = currentMaterial.radiationLength;
0283   mEE_ = currentMaterial.mEE;
0284
0285   if (debugLvl_ > 0) {
0286     debugOut << "     currentMaterial "; currentMaterial.Print();
0287   }
0288
0289   // limit due to momloss
0290   double relMomLossPer_cm(0);
0291   stepSize_ = 1.; // set stepsize for momLoss calculation
0292
0293   if (matZ_ > 1.E-3) { // don't calculate energy loss for vacuum
0294     relMomLossPer_cm = this->momentumLoss(limits.getStepSign(), mom, true) / mom;
0295   }
0296
0297   double maxStepMomLoss = fabs((maxRelMomLoss - fabs(relMomLoss)) / relMomLossPer_cm); // >= 0
0298   limits.setLimit(stp_momLoss, maxStepMomLoss);
0299
0300   if (debugLvl_ > 0) {
0301     debugOut << "     momLoss exceeded after a step of " <<  maxStepMomLoss
0302         << "; relMomLoss up to now = " << relMomLoss << "\n";
0303   }
0304
0305
0306   // now look for boundaries
0307   sMax = limits.getLowestLimitSignedVal();
0308
0309   stepSize_ = limits.getStepSign() * minStep;
0310   M1x3 SA;
0311   double boundaryStep(sMax);
0312
0313   for (unsigned int i=0; i<100; ++i) {
0314     if (debugLvl_ > 0) {
0315       debugOut << "     find next boundary\n";
0316     }
0317     double step =  materialInterface_->findNextBoundary(rep, state7, boundaryStep, varField);
0318
0319     if (debugLvl_ > 0) {
0320       if (step == 0) {
0321         debugOut << "     materialInterface_ returned a step of 0 \n";
0322       }
0323     }
0324
0325     stepSize_ += step;
0326     boundaryStep -= step;
0327
0328     if (debugLvl_ > 0) {
0329       debugOut << "     made a step of " << step << "\n";
0330     }
0331
0332     if (! ignoreBoundariesBetweenEqualMaterials_)
0333       break;
0334
0335     if (fabs(stepSize_) >= fabs(sMax))
0336       break;
0337
0338     // propagate with found step to boundary
0339     rep->RKPropagate(state7, nullptr, SA, step, varField);
0340
0341     // make minStep to cross boundary
0342     state7[0] += limits.getStepSign() * minStep * state7[3];
0343     state7[1] += limits.getStepSign() * minStep * state7[4];
0344     state7[2] += limits.getStepSign() * minStep * state7[5];
0345
0346     materialInterface_->initTrack(state7[0], state7[1], state7[2],
0347                                   limits.getStepSign() * state7[3], limits.getStepSign() * state7[4], limits.getStepSign() * state7[5]);
0348
0349     Material materialAfter = materialInterface_->getMaterialParameters();
0350
0351     if (debugLvl_ > 0) {
0352       debugOut << "     material after step: "; materialAfter.Print();
0353     }
0354
0355     if (materialAfter != currentMaterial)
0356       break;
0357   }
0358
0359   limits.setLimit(stp_boundary, stepSize_);
0360
0361
0362   relMomLoss += relMomLossPer_cm * limits.getLowestLimitVal();
0363 }
0364
0365
0366 void MaterialEffects::getParticleParameters()
0367 {
0368   TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(pdg_);
0369   charge_ = part->Charge() / 3.;  // We only ever use the square
0370   mass_ = part->Mass(); // GeV
0371 }
0372
0373
0374 void MaterialEffects::getMomGammaBeta(double Energy,
0375                      double& mom, double& gammaSquare, double& gamma, double& betaSquare) const {
0376
0377   if (Energy <= mass_) {
0378     Exception exc("MaterialEffects::getMomGammaBeta - Energy <= mass",__LINE__,__FILE__);
0379     exc.setFatal();
0380     throw exc;
0381   }
0382   gamma = Energy/mass_;
0383   gammaSquare = gamma*gamma;
0384   betaSquare = 1.-1./gammaSquare;
0385   mom = Energy*sqrt(betaSquare);
0386 }
0387
0388
0389
0390 //---- Energy-loss and Noise calculations -----------------------------------------
0391
0392 double MaterialEffects::momentumLoss(double stepSign, double mom, bool linear)
0393 {
0394   double E0 = hypot(mom, mass_);
0395   double step = stepSize_*stepSign; // signed
0396
0397
0398   // calc dEdx_, also needed in noiseBetheBloch!
0399   // using fourth order Runge Kutta
0400   //k1 = f(t0,y0)
0401   //k2 = f(t0 + h/2, y0 + h/2 * k1)
0402   //k3 = f(t0 + h/2, y0 + h/2 * k2)
0403   //k4 = f(t0 + h,   y0 + h   * k3)
0404
0405   // This means in our case:
0406   //dEdx1 = dEdx(x0,       E0)
0407   //dEdx2 = dEdx(x0 + h/2, E1); E1 = E0 + h/2 * dEdx1
0408   //dEdx3 = dEdx(x0 + h/2, E2); E2 = E0 + h/2 * dEdx2
0409   //dEdx4 = dEdx(x0 + h,   E3); E3 = E0 + h   * dEdx3
0410
0411   double dEdx1 = dEdx(E0); // dEdx(x0,p0)
0412
0413   if (linear) {
0414     dEdx_ = dEdx1;
0415   }
0416   else { // RK4
0417     double E1 = E0 - dEdx1*step/2.;
0418     double dEdx2 = dEdx(E1); // dEdx(x0 + h/2, E0 + h/2 * dEdx1)
0419
0420     double E2 = E0 - dEdx2*step/2.;
0421     double dEdx3 = dEdx(E2); // dEdx(x0 + h/2, E0 + h/2 * dEdx2)
0422
0423     double E3 = E0 - dEdx3*step;
0424     double dEdx4 = dEdx(E3); // dEdx(x0 + h, E0 + h * dEdx3)
0425
0426     dEdx_ = (dEdx1 + 2.*dEdx2 + 2.*dEdx3 + dEdx4)/6.;
0427   }
0428
0429   E_ = E0 - dEdx_*step*0.5;
0430
0431   double dE = step*dEdx_; // positive for positive stepSign
0432
0433   double momLoss(0);
0434
0435   if (E0 - dE <= mass_) {
0436     // Step would stop particle (E_kin <= 0).
0437     return momLoss = mom;
0438   }
0439   else momLoss = mom - sqrt(pow(E0 - dE, 2) - mass_*mass_); // momLoss; positive for positive stepSign
0440
0441   if (debugLvl_ > 0) {
0442     debugOut << "      MaterialEffects::momentumLoss: mom = " << mom << "; E0 = " << E0
0443         << "; dEdx = " << dEdx_
0444         << "; dE = " << dE << "; mass = " << mass_ << "\n";
0445   }
0446
0447   //assert(momLoss * stepSign >= 0);
0448
0449   return momLoss;
0450 }
0451
0452
0453 double MaterialEffects::dEdx(double Energy) {
0454
0455   double mom(0), gammaSquare(0), gamma(0), betaSquare(0);
0456   this->getMomGammaBeta(Energy, mom, gammaSquare, gamma, betaSquare);
0457   if (pdg_ == c_monopolePDGCode) { // if TParticlePDG also had magnetic charge, life would have been easier.
0458     charge_ = mag_charge_ * sqrt(betaSquare); //effective charge for monopoles
0459   }
0460
0461   double result(0);
0462
0463   if (energyLossBetheBloch_)
0464     result += dEdxBetheBloch(betaSquare, gamma, gammaSquare);
0465
0466   if (energyLossBrems_)
0467     result += dEdxBrems(mom);
0468
0469   return result;
0470 }
0471
0472
0473 double MaterialEffects::dEdxBetheBloch(double betaSquare, double gamma, double gammaSquare) const
0474 {
0475   static const double betaGammaMin(0.05);
0476   if (betaSquare*gammaSquare < betaGammaMin*betaGammaMin) {
0477     Exception exc("MaterialEffects::dEdxBetheBloch ==> beta*gamma < 0.05, Bethe-Bloch implementation not valid anymore!",__LINE__,__FILE__);
0478     exc.setFatal();
0479     throw exc;
0480   }
0481
0482   // calc dEdx_, also needed in noiseBetheBloch!
0483   double result( 0.307075 * matZ_ / matA_ * matDensity_ / betaSquare * charge_ * charge_ );
0484   double massRatio( me_ / mass_ );
0485   double argument( gammaSquare * betaSquare * me_ * 1.E3 * 2. / ((1.E-6 * mEE_) *
0486       sqrt(1. + 2. * gamma * massRatio + massRatio * massRatio)) );
0487   result *= log(argument) - betaSquare; // Bethe-Bloch [MeV/cm]
0488   result *= 1.E-3;  // in GeV/cm, hence 1.e-3
0489   if (result < 0.) {
0490     result = 0;
0491   }
0492
0493   return result;
0494 }
0495
0496
0497 void MaterialEffects::noiseBetheBloch(M7x7& noise, double mom, double betaSquare, double gamma, double gammaSquare) const
0498 {
0499   // Code ported from GEANT 3 (erland.F)
0500
0501   // ENERGY LOSS FLUCTUATIONS; calculate sigma^2(E);
0502   double sigma2E ( 0. );
0503   double zeta  ( 153.4E3 * charge_ * charge_ / betaSquare * matZ_ / matA_ * matDensity_ * fabs(stepSize_) ); // eV
0504   double Emax  ( 2.E9 * me_ * betaSquare * gammaSquare / (1. + 2.*gamma * me_ / mass_ + (me_ / mass_) * (me_ / mass_)) ); // eV
0505   double kappa ( zeta / Emax );
0506
0507   if (kappa > 0.01) { // Vavilov-Gaussian regime
0508     sigma2E += zeta * Emax * (1. - betaSquare / 2.); // eV^2
0509   } else { // Urban/Landau approximation
0510     // calculate number of collisions Nc
0511     double I = 16. * pow(matZ_, 0.9); // eV
0512     double f2 = 0.;
0513     if (matZ_ > 2.) f2 = 2. / matZ_;
0514     double f1 = 1. - f2;
0515     double e2 = 10.*matZ_ * matZ_; // eV
0516     double e1 = pow((I / pow(e2, f2)), 1. / f1); // eV
0517
0518     double mbbgg2 = 2.E9 * mass_ * betaSquare * gammaSquare; // eV
0519     double Sigma1 = dEdx_ * 1.0E9 * f1 / e1 * (log(mbbgg2 / e1) - betaSquare) / (log(mbbgg2 / I) - betaSquare) * 0.6; // 1/cm
0520     double Sigma2 = dEdx_ * 1.0E9 * f2 / e2 * (log(mbbgg2 / e2) - betaSquare) / (log(mbbgg2 / I) - betaSquare) * 0.6; // 1/cm
0521     double Sigma3 = dEdx_ * 1.0E9 * Emax / (I * (Emax + I) * log((Emax + I) / I)) * 0.4; // 1/cm
0522
0523     double Nc = (Sigma1 + Sigma2 + Sigma3) * fabs(stepSize_);
0524
0525     if (Nc > 50.) { // truncated Landau distribution
0526       double sigmaalpha = 15.76;
0527       // calculate sigmaalpha  (see GEANT3 manual W5013)
0528       double RLAMED = -0.422784 - betaSquare - log(zeta / Emax);
0529       double RLAMAX =  0.60715 + 1.1934 * RLAMED + (0.67794 + 0.052382 * RLAMED) * exp(0.94753 + 0.74442 * RLAMED);
0530       // from lambda max to sigmaalpha=sigma (empirical polynomial)
0531       if (RLAMAX <= 1010.) {
0532         sigmaalpha =  1.975560
0533                       + 9.898841e-02 * RLAMAX
0534                       - 2.828670e-04 * RLAMAX * RLAMAX
0535                       + 5.345406e-07 * pow(RLAMAX, 3.)
0536                       - 4.942035e-10 * pow(RLAMAX, 4.)
0537                       + 1.729807e-13 * pow(RLAMAX, 5.);
0538       } else { sigmaalpha = 1.871887E+01 + 1.296254E-02 * RLAMAX; }
0539       // alpha=54.6  corresponds to a 0.9996 maximum cut
0540       if (sigmaalpha > 54.6) sigmaalpha = 54.6;
0541       sigma2E += sigmaalpha * sigmaalpha * zeta * zeta; // eV^2
0542     } else { // Urban model
0543       static const double alpha = 0.996;
0544       double Ealpha  = I / (1. - (alpha * Emax / (Emax + I))); // eV
0545       double meanE32 = I * (Emax + I) / Emax * (Ealpha - I); // eV^2
0546       sigma2E += fabs(stepSize_) * (Sigma1 * e1 * e1 + Sigma2 * e2 * e2 + Sigma3 * meanE32); // eV^2
0547     }
0548   }
0549
0550   sigma2E *= 1.E-18; // eV -> GeV
0551
0552   // update noise matrix, using linear error propagation from E to q/p
0553   noise[6 * 7 + 6] += charge_*charge_/betaSquare / pow(mom, 4) * sigma2E;
0554 }
0555
0556
0557 void MaterialEffects::noiseCoulomb(M7x7& noise,
0558                                    const M1x3& direction, double momSquare, double betaSquare) const
0559 {
0560
0561   // MULTIPLE SCATTERING; calculate sigma^2
0562   double sigma2 = 0;
0563   assert(mscModelCode_ == 0 || mscModelCode_ == 1);
0564   const double step = fabs(stepSize_);
0565   const double step2 = step * step;
0566   if (mscModelCode_ == 0) {// PANDA report PV/01-07 eq(43); linear in step length
0567     sigma2 = 225.E-6 * charge_ * charge_ / (betaSquare * momSquare) * step / radiationLength_ * matZ_ / (matZ_ + 1) * log(159.*pow(matZ_, -1. / 3.)) / log(287.*pow(matZ_, -0.5)); // sigma^2 = 225E-6*z^2/mom^2 * XX0/beta_^2 * Z/(Z+1) * ln(159*Z^(-1/3))/ln(287*Z^(-1/2)
0568
0569   } else if (mscModelCode_ == 1) { //Highland not linear in step length formula taken from PDG book 2011 edition
0570     double stepOverRadLength = step / radiationLength_;
0571     double logCor = (1 + 0.038 * log(stepOverRadLength));
0572     sigma2 = 0.0136 * 0.0136 * charge_ * charge_ / (betaSquare * momSquare) * stepOverRadLength * logCor * logCor;
0573   }
0574   //assert(sigma2 >= 0.0);
0575   sigma2 = (sigma2 > 0.0 ? sigma2 : 0.0);
0576   //XXX debugOut << "MaterialEffects::noiseCoulomb the MSC variance is " << sigma2 << std::endl;
0577
0578   M7x7 noiseAfter; // will hold the new MSC noise to cause by the current stepSize_ length
0579   std::fill(noiseAfter.begin(), noiseAfter.end(), 0);
0580
0581   const M1x3& a = direction; // as an abbreviation
0582   // This calculates the MSC angular spread in the 7D global
0583   // coordinate system.  See PDG 2010, Sec. 27.3 for formulae.
0584   noiseAfter[0 * 7 + 0] =  sigma2 * step2 / 3.0 * (1 - a[0]*a[0]);
0585   noiseAfter[1 * 7 + 0] = -sigma2 * step2 / 3.0 * a[0]*a[1];
0586   noiseAfter[2 * 7 + 0] = -sigma2 * step2 / 3.0 * a[0]*a[2];
0587   noiseAfter[3 * 7 + 0] =  sigma2 * step * 0.5 * (1 - a[0]*a[0]);
0588   noiseAfter[4 * 7 + 0] = -sigma2 * step * 0.5 * a[0]*a[1];
0589   noiseAfter[5 * 7 + 0] = -sigma2 * step * 0.5 * a[0]*a[1];
0590   noiseAfter[0 * 7 + 1] = noiseAfter[1 * 7 + 0];
0591   noiseAfter[1 * 7 + 1] =  sigma2 * step2 / 3.0 * (1 - a[1]*a[1]);
0592   noiseAfter[2 * 7 + 1] = -sigma2 * step2 / 3.0 * a[1]*a[2];
0593   noiseAfter[3 * 7 + 1] = noiseAfter[4 * 7 + 0]; // Cov(x,a_y) = Cov(y,a_x)
0594   noiseAfter[4 * 7 + 1] =  sigma2 * step * 0.5 * (1 - a[1] * a[1]);
0595   noiseAfter[5 * 7 + 1] = -sigma2 * step * 0.5 * a[1]*a[2];
0596   noiseAfter[0 * 7 + 2] = noiseAfter[2 * 7 + 0];
0597   noiseAfter[1 * 7 + 2] = noiseAfter[2 * 7 + 1];
0598   noiseAfter[2 * 7 + 2] =  sigma2 * step2 / 3.0 * (1 - a[2]*a[2]);
0599   noiseAfter[3 * 7 + 2] = noiseAfter[5 * 7 + 0]; // Cov(z,a_x) = Cov(x,a_z)
0600   noiseAfter[4 * 7 + 2] = noiseAfter[5 * 7 + 1]; // Cov(y,a_z) = Cov(z,a_y)
0601   noiseAfter[5 * 7 + 2] =  sigma2 * step * 0.5 * (1 - a[2]*a[2]);
0602   noiseAfter[0 * 7 + 3] = noiseAfter[3 * 7 + 0];
0603   noiseAfter[1 * 7 + 3] = noiseAfter[3 * 7 + 1];
0604   noiseAfter[2 * 7 + 3] = noiseAfter[3 * 7 + 2];
0605   noiseAfter[3 * 7 + 3] =  sigma2 * (1 - a[0]*a[0]);
0606   noiseAfter[4 * 7 + 3] = -sigma2 * a[0]*a[1];
0607   noiseAfter[5 * 7 + 3] = -sigma2 * a[0]*a[2];
0608   noiseAfter[0 * 7 + 4] = noiseAfter[4 * 7 + 0];
0609   noiseAfter[1 * 7 + 4] = noiseAfter[4 * 7 + 1];
0610   noiseAfter[2 * 7 + 4] = noiseAfter[4 * 7 + 2];
0611   noiseAfter[3 * 7 + 4] = noiseAfter[4 * 7 + 3];
0612   noiseAfter[4 * 7 + 4] =  sigma2 * (1 - a[1]*a[1]);
0613   noiseAfter[5 * 7 + 4] = -sigma2 * a[1]*a[2];
0614   noiseAfter[0 * 7 + 5] = noiseAfter[5 * 7 + 0];
0615   noiseAfter[1 * 7 + 5] = noiseAfter[5 * 7 + 1];
0616   noiseAfter[2 * 7 + 5] = noiseAfter[5 * 7 + 2];
0617   noiseAfter[3 * 7 + 5] = noiseAfter[5 * 7 + 3];
0618   noiseAfter[4 * 7 + 5] = noiseAfter[5 * 7 + 4];
0619   noiseAfter[5 * 7 + 5] = sigma2 * (1 - a[2]*a[2]);
0620 //    debugOut << "new noise\n";
0621 //    RKTools::printDim(noiseAfter, 7,7);
0622   for (unsigned int i = 0; i < 7 * 7; ++i) {
0623     noise[i] += noiseAfter[i];
0624   }
0625 }
0626
0627
0628 double MaterialEffects::dEdxBrems(double mom) const
0629 {
0630
0631   // Code ported from GEANT 3 (gbrele.F)
0632
0633   if (abs(pdg_) != 11) return 0; // only for electrons and positrons
0634
0635 #if !defined(BETHE)
0636   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};
0637   static const double xi = 2.51, beta = 0.99, vl = 0.00004;
0638 #endif
0639 #if defined(BETHE) // no MIGDAL corrections
0640   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};
0641   static const double xi = 2.10, beta = 1.00, vl = 0.001;
0642 #endif
0643
0644   double BCUT = 10000.; // energy up to which soft bremsstrahlung energy loss is calculated
0645
0646   static const double THIGH = 100., CHIGH = 50.;
0647   double dedxBrems = 0.;
0648
0649   if (BCUT > 0.) {
0650     double T, kc;
0651
0652     if (BCUT > mom) BCUT = mom; // confine BCUT to mom_
0653
0654     // T=mom_,  confined to THIGH
0655     // kc=BCUT, confined to CHIGH ??
0656     if (mom > THIGH) {
0657       T = THIGH;
0658       if (BCUT >= THIGH)
0659         kc = CHIGH;
0660       else
0661         kc = BCUT;
0662     } else {
0663       T = mom;
0664       kc = BCUT;
0665     }
0666
0667     double E = T + me_; // total electron energy
0668     if (BCUT > T)
0669       kc = T;
0670
0671     double X = log(T / me_);
0672     double Y = log(kc / (E * vl));
0673
0674     double XX;
0675     int    K;
0676     double S = 0., YY = 1.;
0677
0678     for (unsigned int I = 1; I <= 2; ++I) {
0679       XX = 1.;
0680       for (unsigned int J = 1; J <= 6; ++J) {
0681         K = 6 * I + J - 6;
0682         S += C[K] * XX * YY;
0683         XX *= X;
0684       }
0685       YY *= Y;
0686     }
0687
0688     for (unsigned int I = 3; I <= 6; ++I) {
0689       XX = 1.;
0690       for (unsigned int J = 1; J <= 6; ++J) {
0691         K = 6 * I + J - 6;
0692         if (Y > 0.)
0693           K += 24;
0694         S += C[K] * XX * YY;
0695         XX *= X;
0696       }
0697       YY *= Y;
0698     }
0699
0700     double SS = 0.;
0701     YY = 1.;
0702
0703     for (unsigned int I = 1; I <= 2; ++I) {
0704       XX = 1.;
0705       for (unsigned int J = 1; J <= 5; ++J) {
0706         K = 5 * I + J + 55;
0707         SS += C[K] * XX * YY;
0708         XX *= X;
0709       }
0710       YY *= Y;
0711     }
0712
0713     for (unsigned int I = 3; I <= 5; ++I) {
0714       XX = 1.;
0715       for (unsigned int J = 1; J <= 5; ++J) {
0716         K = 5 * I + J + 55;
0717         if (Y > 0.)
0718           K += 15;
0719         SS += C[K] * XX * YY;
0720         XX *= X;
0721       }
0722       YY *= Y;
0723     }
0724
0725     S += matZ_ * SS;
0726
0727     if (S > 0.) {
0728       double CORR = 1.;
0729 #if !defined(BETHE)
0730       CORR = 1. / (1. + 0.805485E-10 * matDensity_ * matZ_ * E * E / (matA_ * kc * kc)); // MIGDAL correction factor
0731 #endif
0732
0733       // We use exp(beta * log(...) here because pow(..., beta) is
0734       // REALLY slow and we don't need ultimate numerical precision
0735       // for this approximation.
0736       double FAC = matZ_ * (matZ_ + xi) * E * E / (E + me_);
0737       if (beta == 1.)  // That is the #ifdef BETHE case
0738         FAC *= kc * CORR / T;
0739       else
0740         FAC *= exp(beta * log(kc * CORR / T));
0741       if (FAC <= 0.)
0742         return 0.;
0743       dedxBrems = FAC * S;
0744
0745
0746       if (mom >= THIGH) {
0747         double RAT;
0748         if (BCUT < THIGH) {
0749           RAT = BCUT / mom;
0750           S = (1. - 0.5 * RAT + 2.*RAT * RAT / 9.);
0751           RAT = BCUT / T;
0752           S /= 1. - 0.5 * RAT + 2.*RAT * RAT / 9.;
0753         } else {
0754           RAT = BCUT / mom;
0755           S = BCUT * (1. - 0.5 * RAT + 2.*RAT * RAT / 9.);
0756           RAT = kc / T;
0757           S /= kc * (1. - 0.5 * RAT + 2.*RAT * RAT / 9.);
0758         }
0759         dedxBrems *= S; // GeV barn
0760       }
0761
0762       dedxBrems *= 0.60221367 * matDensity_ / matA_; // energy loss dE/dx [GeV/cm]
0763     }
0764   }
0765
0766   if (dedxBrems < 0.)
0767     dedxBrems = 0;
0768
0769   double factor = 1.; // positron correction factor
0770
0771   if (pdg_ == -11) {
0772     static const double AA = 7522100., A1 = 0.415, A3 = 0.0021, A5 = 0.00054;
0773
0774     double ETA = 0.;
0775     if (matZ_ > 0.) {
0776       double X = log(AA * mom / (matZ_ * matZ_));
0777       if (X > -8.) {
0778         if (X >= +9.) ETA = 1.;
0779         else {
0780           double W = A1 * X + A3 * pow(X, 3.) + A5 * pow(X, 5.);
0781           ETA = 0.5 + atan(W) / M_PI;
0782         }
0783       }
0784     }
0785
0786     if (ETA < 0.0001)
0787       factor = 1.E-10;
0788     else if (ETA > 0.9999)
0789       factor = 1.;
0790     else {
0791       double E0 = BCUT / mom;
0792       if (E0 > 1.)
0793         E0 = 1.;
0794       if (E0 < 1.E-8)
0795         factor = 1.;
0796       else
0797         factor = ETA * (1. - pow(1. - E0, 1. / ETA)) / E0;
0798     }
0799   }
0800
0801   return factor * dedxBrems; //always positive
0802 }
0803
0804
0805 void MaterialEffects::noiseBrems(M7x7& noise, double momSquare, double betaSquare) const
0806 {
0807   // Code ported from GEANT 3 (erland.F) and simplified
0808   // E \approx p is assumed.
0809   // the factor  1.44 is not in the original Bethe-Heitler model.
0810   // It seems to be some empirical correction copied over from some other project.
0811
0812   if (abs(pdg_) != 11) return; // only for electrons and positrons
0813
0814   double minusXOverLn2  = -1.442695 * fabs(stepSize_) / radiationLength_;
0815   double sigma2E = 1.44*(pow(3., minusXOverLn2) - pow(4., minusXOverLn2)) * momSquare;
0816   sigma2E = std::max(sigma2E, 0.0); // must be positive
0817
0818   // update noise matrix, using linear error propagation from E to q/p
0819   noise[6 * 7 + 6] += charge_*charge_/betaSquare / pow(momSquare, 2) * sigma2E;
0820 }
0821
0822
0823 void MaterialEffects::setDebugLvl(unsigned int lvl) {
0824   debugLvl_ = lvl;
0825   if (materialInterface_ and debugLvl_ > 1)
0826     materialInterface_->setDebugLvl(debugLvl_-1);
0827 }
0828
0829
0830 void MaterialEffects::drawdEdx(int pdg) {
0831   pdg_ = pdg;
0832   this->getParticleParameters();
0833
0834   stepSize_ = 1;
0835
0836   materialInterface_->initTrack(0, 0, 0, 1, 1, 1);
0837   auto currentMaterial = materialInterface_->getMaterialParameters();
0838   matDensity_ = currentMaterial.density;
0839   matZ_ = currentMaterial.Z;
0840   matA_ = currentMaterial.A;
0841   radiationLength_ = currentMaterial.radiationLength;
0842   mEE_ = currentMaterial.mEE;
0843
0844   double minMom = 0.00001;
0845   double maxMom = 10000;
0846   int nSteps(10000);
0847   double logStepSize = (log10(maxMom) - log10(minMom)) / (nSteps-1);
0848
0849   TH1D hdEdxBethe("dEdxBethe", "dEdxBethe; log10(mom)", nSteps, log10(minMom), log10(maxMom));
0850   TH1D hdEdxBrems("dEdxBrems", "dEdxBrems; log10(mom)", nSteps, log10(minMom), log10(maxMom));
0851
0852   for (int i=0; i<nSteps; ++i) {
0853     double mom = pow(10., log10(minMom) + i*logStepSize);
0854     double E = hypot(mom, mass_);
0855     if (pdg_ == c_monopolePDGCode) {
0856       charge_ = mag_charge_ * mom / E; //effective charge for monopoles
0857     }
0858
0859     energyLossBrems_ = false;
0860     energyLossBetheBloch_ = true;
0861
0862     try {
0863       hdEdxBethe.Fill(log10(mom), dEdx(E));
0864     }
0865     catch (...) {
0866
0867     }
0868
0869
0870     //debugOut<< "E = " << E << "; dEdx = " << dEdx(E) <<"\n";
0871
0872     energyLossBrems_ = true;
0873     energyLossBetheBloch_ = false;
0874     try {
0875       hdEdxBrems.Fill(log10(mom), dEdx(E));
0876     }
0877     catch (...) {
0878
0879     }
0880   }
0881
0882   energyLossBrems_ = true;
0883   energyLossBetheBloch_ = true;
0884
0885   std::string Result;//string which will contain the result
0886   std::stringstream convert; // stringstream used for the conversion
0887   convert << pdg;//add the value of Number to the characters in the stream
0888   Result = convert.str();//set Result to the content of the stream
0889
0890   TFile outfile("dEdx_" + TString(Result) + ".root", "recreate");
0891   outfile.cd();
0892   hdEdxBethe.Write();
0893   hdEdxBrems.Write();
0894   outfile.Close();
0895 }
0896
0897 } /* End of namespace genfit */
0898
0899