GENFIT  Rev:NoNumberAvailable
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules
MaterialEffects.cc
Go to the documentation of this file.
1 /* Copyright 2008-2014, Technische Universitaet Muenchen,
2  Authors: Christian Hoeppner & Sebastian Neubert
3 
4  This file is part of GENFIT.
5 
6  GENFIT is free software: you can redistribute it and/or modify
7  it under the terms of the GNU Lesser General Public License as published
8  by the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  GENFIT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU Lesser General Public License for more details.
15 
16  You should have received a copy of the GNU Lesser General Public License
17  along with GENFIT. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #include "MaterialEffects.h"
21 #include "Exception.h"
22 #include "IO.h"
23 
24 #include <stdexcept>
25 #include <string>
26 #include <stdlib.h>
27 #include <math.h>
28 #include <assert.h>
29 
30 #include <TDatabasePDG.h>
31 #include <TMath.h>
32 
33 #include <TH1D.h>
34 #include <TFile.h>
35 
36 
37 namespace genfit {
38 
39 MaterialEffects* MaterialEffects::instance_ = nullptr;
40 
41 
43  noEffects_(false),
44  energyLossBetheBloch_(true), noiseBetheBloch_(true),
45  noiseCoulomb_(true),
46  energyLossBrems_(true), noiseBrems_(true),
47  ignoreBoundariesBetweenEqualMaterials_(true),
48  me_(0.510998910E-3),
49  stepSize_(0),
50  dEdx_(0),
51  E_(0),
52  matDensity_(0),
53  matZ_(0),
54  matA_(0),
55  radiationLength_(0),
56  mEE_(0),
57  pdg_(0),
58  charge_(0),
59  mass_(0),
60  mscModelCode_(0),
61  materialInterface_(nullptr),
62  debugLvl_(0)
63 {
64 }
65 
67 {
68  if (materialInterface_ != nullptr) delete materialInterface_;
69 }
70 
72 {
73  if (instance_ == nullptr) instance_ = new MaterialEffects();
74  return instance_;
75 }
76 
78 {
79  if (instance_ != nullptr) {
80  delete instance_;
81  instance_ = nullptr;
82  }
83 }
84 
86 {
87  if (materialInterface_ != nullptr) {
88  std::string msg("MaterialEffects::initMaterialInterface(): Already initialized! ");
89  std::runtime_error err(msg);
90  }
91  materialInterface_ = matIfc;
92 }
93 
94 
95 
96 void MaterialEffects::setMscModel(const std::string& modelName)
97 {
98  if (modelName == std::string("GEANE")) {
99  mscModelCode_ = 0;
100  } else if (modelName == std::string("Highland")) {
101  mscModelCode_ = 1;
102  } else {// throw exception
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__);
105  exc.setFatal();
106  errorOut << exc.what();
107  throw exc;
108  }
109 }
110 
111 
112 double MaterialEffects::effects(const std::vector<RKStep>& steps,
113  int materialsFXStart,
114  int materialsFXStop,
115  const double& mom,
116  const int& pdg,
117  M7x7* noise)
118 {
119 
120  if (debugLvl_ > 0) {
121  debugOut << " MaterialEffects::effects \n";
122  }
123 
124  /*debugOut << "noEffects_ " << noEffects_ << "\n";
125  debugOut << "energyLossBetheBloch_ " << energyLossBetheBloch_ << "\n";
126  debugOut << "noiseBetheBloch_ " << noiseBetheBloch_ << "\n";
127  debugOut << "noiseCoulomb_ " << noiseCoulomb_ << "\n";
128  debugOut << "energyLossBrems_ " << energyLossBrems_ << "\n";
129  debugOut << "noiseBrems_ " << noiseBrems_ << "\n";*/
130 
131 
132  if (noEffects_) return 0.;
133 
134  if (materialInterface_ == nullptr) {
135  std::string msg("MaterialEffects hasn't been initialized with a correct AbsMaterialInterface pointer!");
136  std::runtime_error err(msg);
137  throw err;
138  }
139 
140  bool doNoise(noise != nullptr);
141 
142  pdg_ = pdg;
144 
145  double momLoss = 0.;
146 
147  for ( std::vector<RKStep>::const_iterator it = steps.begin() + materialsFXStart; it != steps.begin() + materialsFXStop; ++it) { // loop over steps
148 
149  double realPath = it->matStep_.stepSize_;
150  if (fabs(realPath) < 1.E-8) {
151  // do material effects only if distance is not too small
152  continue;
153  }
154 
155  if (debugLvl_ > 0) {
156  debugOut << " calculate matFX ";
157  if (doNoise)
158  debugOut << "and noise";
159  debugOut << " for ";
160  debugOut << "stepSize = " << it->matStep_.stepSize_ << "\t";
161  it->matStep_.materialProperties_.Print();
162  }
163 
164  double stepSign(1.);
165  if (realPath < 0)
166  stepSign = -1.;
167  realPath = fabs(realPath);
168  stepSize_ = realPath;
169 
170  it->matStep_.materialProperties_.getMaterialProperties(matDensity_, matZ_, matA_, radiationLength_, mEE_);
171 
172  if (matZ_ > 1.E-3) { // don't calculate energy loss for vacuum
173 
174  momLoss += momentumLoss(stepSign, mom - momLoss, false);
175 
176  if (doNoise){
177  // get values for the "effective" energy of the RK step E_
178  double p(0), gammaSquare(0), gamma(0), betaSquare(0);
179  this->getMomGammaBeta(E_, p, gammaSquare, gamma, betaSquare);
180  double pSquare = p*p;
181 
183  this->noiseBetheBloch(*noise, p, betaSquare, gamma, gammaSquare);
184 
185  if (noiseCoulomb_)
186  this->noiseCoulomb(*noise, *((M1x3*) &it->state7_[3]), pSquare, betaSquare);
187 
189  this->noiseBrems(*noise, pSquare, betaSquare);
190  } // end doNoise
191 
192  }
193 
194  } // end loop over steps
195 
196  if (momLoss >= mom) {
197  Exception exc("MaterialEffects::effects ==> momLoss >= momentum, aborting extrapolation!",__LINE__,__FILE__);
198  exc.setFatal();
199  throw exc;
200  }
201 
202  return momLoss;
203 }
204 
205 
207  M1x7& state7,
208  const double& mom, // momentum
209  double& relMomLoss, // relative momloss for the step will be added
210  const int& pdg,
211  MaterialProperties& currentMaterial,
212  StepLimits& limits,
213  bool varField)
214 {
215 
216  static const double maxRelMomLoss = .01; // maximum relative momentum loss allowed
217  static const double Pmin = 4.E-3; // minimum momentum for propagation [GeV]
218  static const double minStep = 1.E-4; // 1 µm
219 
220  // check momentum
221  if(mom < Pmin){
222  std::ostringstream sstream;
223  sstream << "MaterialEffects::stepper ==> momentum too low: " << mom*1000. << " MeV";
224  Exception exc(sstream.str(),__LINE__,__FILE__);
225  exc.setFatal();
226  throw exc;
227  }
228 
229  // Trivial cases
230  if (noEffects_)
231  return;
232 
233  if (materialInterface_ == nullptr) {
234  std::string msg("MaterialEffects hasn't been initialized with a correct AbsMaterialInterface pointer!");
235  std::runtime_error err(msg);
236  throw err;
237  }
238 
239  if (relMomLoss > maxRelMomLoss) {
240  limits.setLimit(stp_momLoss, 0);
241  return;
242  }
243 
244 
245  double sMax = limits.getLowestLimitSignedVal(); // signed
246 
247  if (fabs(sMax) < minStep)
248  return;
249 
250 
251 
252  pdg_ = pdg;
254 
255 
256  // make 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];
260 
261  materialInterface_->initTrack(state7[0], state7[1], state7[2],
262  limits.getStepSign() * state7[3], limits.getStepSign() * state7[4], limits.getStepSign() * state7[5]);
263 
266 
267 
268  if (debugLvl_ > 0) {
269  debugOut << " currentMaterial "; currentMaterial.Print();
270  }
271 
272  // limit due to momloss
273  double relMomLossPer_cm(0);
274  stepSize_ = 1.; // set stepsize for momLoss calculation
275 
276  if (matZ_ > 1.E-3) { // don't calculate energy loss for vacuum
277  relMomLossPer_cm = this->momentumLoss(limits.getStepSign(), mom, true) / mom;
278  }
279 
280  double maxStepMomLoss = fabs((maxRelMomLoss - fabs(relMomLoss)) / relMomLossPer_cm); // >= 0
281  limits.setLimit(stp_momLoss, maxStepMomLoss);
282 
283  if (debugLvl_ > 0) {
284  debugOut << " momLoss exceeded after a step of " << maxStepMomLoss
285  << "; relMomLoss up to now = " << relMomLoss << "\n";
286  }
287 
288 
289  // now look for boundaries
290  sMax = limits.getLowestLimitSignedVal();
291 
292  stepSize_ = limits.getStepSign() * minStep;
293  MaterialProperties materialAfter;
294  M1x3 SA;
295  double boundaryStep(sMax);
296 
297  for (unsigned int i=0; i<100; ++i) {
298  if (debugLvl_ > 0) {
299  debugOut << " find next boundary\n";
300  }
301  double step = materialInterface_->findNextBoundary(rep, state7, boundaryStep, varField);
302 
303  if (debugLvl_ > 0) {
304  if (step == 0) {
305  debugOut << " materialInterface_ returned a step of 0 \n";
306  }
307  }
308 
309  stepSize_ += step;
310  boundaryStep -= step;
311 
312  if (debugLvl_ > 0) {
313  debugOut << " made a step of " << step << "\n";
314  }
315 
317  break;
318 
319  if (fabs(stepSize_) >= fabs(sMax))
320  break;
321 
322  // propagate with found step to boundary
323  rep->RKPropagate(state7, NULL, SA, step, varField);
324 
325  // make minStep to cross boundary
326  state7[0] += limits.getStepSign() * minStep * state7[3];
327  state7[1] += limits.getStepSign() * minStep * state7[4];
328  state7[2] += limits.getStepSign() * minStep * state7[5];
329 
330  materialInterface_->initTrack(state7[0], state7[1], state7[2],
331  limits.getStepSign() * state7[3], limits.getStepSign() * state7[4], limits.getStepSign() * state7[5]);
332 
334 
335  if (debugLvl_ > 0) {
336  debugOut << " material after step: "; materialAfter.Print();
337  }
338 
339  if (materialAfter != currentMaterial)
340  break;
341  }
342 
344 
345 
346  relMomLoss += relMomLossPer_cm * limits.getLowestLimitVal();
347 }
348 
349 
351 {
352  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(pdg_);
353  charge_ = int(part->Charge() / 3.); // We only ever use the square
354  mass_ = part->Mass(); // GeV
355 }
356 
357 
359  double& mom, double& gammaSquare, double& gamma, double& betaSquare) const {
360 
361  if (Energy <= mass_) {
362  Exception exc("MaterialEffects::getMomGammaBeta - Energy <= mass",__LINE__,__FILE__);
363  exc.setFatal();
364  throw exc;
365  }
366  gamma = Energy/mass_;
367  gammaSquare = gamma*gamma;
368  betaSquare = 1.-1./gammaSquare;
369  mom = Energy*sqrt(betaSquare);
370 }
371 
372 
373 
374 //---- Energy-loss and Noise calculations -----------------------------------------
375 
376 double MaterialEffects::momentumLoss(double stepSign, double mom, bool linear)
377 {
378  double E0 = hypot(mom, mass_);
379  double step = stepSize_*stepSign; // signed
380 
381 
382  // calc dEdx_, also needed in noiseBetheBloch!
383  // using fourth order Runge Kutta
384  //k1 = f(t0,y0)
385  //k2 = f(t0 + h/2, y0 + h/2 * k1)
386  //k3 = f(t0 + h/2, y0 + h/2 * k2)
387  //k4 = f(t0 + h, y0 + h * k3)
388 
389  // This means in our case:
390  //dEdx1 = dEdx(x0, E0)
391  //dEdx2 = dEdx(x0 + h/2, E1); E1 = E0 + h/2 * dEdx1
392  //dEdx3 = dEdx(x0 + h/2, E2); E2 = E0 + h/2 * dEdx2
393  //dEdx4 = dEdx(x0 + h, E3); E3 = E0 + h * dEdx3
394 
395  double dEdx1 = dEdx(E0); // dEdx(x0,p0)
396 
397  if (linear) {
398  dEdx_ = dEdx1;
399  }
400  else { // RK4
401  double E1 = E0 - dEdx1*step/2.;
402  double dEdx2 = dEdx(E1); // dEdx(x0 + h/2, E0 + h/2 * dEdx1)
403 
404  double E2 = E0 - dEdx2*step/2.;
405  double dEdx3 = dEdx(E2); // dEdx(x0 + h/2, E0 + h/2 * dEdx2)
406 
407  double E3 = E0 - dEdx3*step;
408  double dEdx4 = dEdx(E3); // dEdx(x0 + h, E0 + h * dEdx3)
409 
410  dEdx_ = (dEdx1 + 2.*dEdx2 + 2.*dEdx3 + dEdx4)/6.;
411  }
412 
413  E_ = E0 - dEdx_*step*0.5;
414 
415  double dE = step*dEdx_; // positive for positive stepSign
416 
417  double momLoss(0);
418 
419  if (E0 - dE <= mass_) {
420  // Step would stop particle (E_kin <= 0).
421  return momLoss = mom;
422  }
423  else momLoss = mom - sqrt(pow(E0 - dE, 2) - mass_*mass_); // momLoss; positive for positive stepSign
424 
425  if (debugLvl_ > 0) {
426  debugOut << " MaterialEffects::momentumLoss: mom = " << mom << "; E0 = " << E0
427  << "; dEdx = " << dEdx_
428  << "; dE = " << dE << "; mass = " << mass_ << "\n";
429  }
430 
431  //assert(momLoss * stepSign >= 0);
432 
433  return momLoss;
434 }
435 
436 
437 double MaterialEffects::dEdx(double Energy) const {
438 
439  double mom(0), gammaSquare(0), gamma(0), betaSquare(0);
440  this->getMomGammaBeta(Energy, mom, gammaSquare, gamma, betaSquare);
441 
442  double result(0);
443 
445  result += dEdxBetheBloch(betaSquare, gamma, gammaSquare);
446 
447  if (energyLossBrems_)
448  result += dEdxBrems(mom);
449 
450  return result;
451 }
452 
453 
454 double MaterialEffects::dEdxBetheBloch(double betaSquare, double gamma, double gammaSquare) const
455 {
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__);
459  exc.setFatal();
460  throw exc;
461  }
462 
463  // calc dEdx_, also needed in noiseBetheBloch!
464  double result( 0.307075 * matZ_ / matA_ * matDensity_ / betaSquare * charge_ * charge_ );
465  double massRatio( me_ / mass_ );
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; // Bethe-Bloch [MeV/cm]
469  result *= 1.E-3; // in GeV/cm, hence 1.e-3
470  if (result < 0.) {
471  result = 0;
472  }
473 
474  return result;
475 }
476 
477 
478 void MaterialEffects::noiseBetheBloch(M7x7& noise, double mom, double betaSquare, double gamma, double gammaSquare) const
479 {
480  // Code ported from GEANT 3
481 
482  // ENERGY LOSS FLUCTUATIONS; calculate sigma^2(E);
483  double sigma2E ( 0. );
484  double zeta ( 153.4E3 * charge_ * charge_ / betaSquare * matZ_ / matA_ * matDensity_ * fabs(stepSize_) ); // eV
485  double Emax ( 2.E9 * me_ * betaSquare * gammaSquare / (1. + 2.*gamma * me_ / mass_ + (me_ / mass_) * (me_ / mass_)) ); // eV
486  double kappa ( zeta / Emax );
487 
488  if (kappa > 0.01) { // Vavilov-Gaussian regime
489  sigma2E += zeta * Emax * (1. - betaSquare / 2.); // eV^2
490  } else { // Urban/Landau approximation
491  // calculate number of collisions Nc
492  double I = 16. * pow(matZ_, 0.9); // eV
493  double f2 = 0.;
494  if (matZ_ > 2.) f2 = 2. / matZ_;
495  double f1 = 1. - f2;
496  double e2 = 10.*matZ_ * matZ_; // eV
497  double e1 = pow((I / pow(e2, f2)), 1. / f1); // eV
498 
499  double mbbgg2 = 2.E9 * mass_ * betaSquare * gammaSquare; // eV
500  double Sigma1 = dEdx_ * 1.0E9 * f1 / e1 * (log(mbbgg2 / e1) - betaSquare) / (log(mbbgg2 / I) - betaSquare) * 0.6; // 1/cm
501  double Sigma2 = dEdx_ * 1.0E9 * f2 / e2 * (log(mbbgg2 / e2) - betaSquare) / (log(mbbgg2 / I) - betaSquare) * 0.6; // 1/cm
502  double Sigma3 = dEdx_ * 1.0E9 * Emax / (I * (Emax + I) * log((Emax + I) / I)) * 0.4; // 1/cm
503 
504  double Nc = (Sigma1 + Sigma2 + Sigma3) * fabs(stepSize_);
505 
506  if (Nc > 50.) { // truncated Landau distribution
507  double sigmaalpha = 15.76;
508  // calculate sigmaalpha (see GEANT3 manual W5013)
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);
511  // from lambda max to sigmaalpha=sigma (empirical polynomial)
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; }
520  // alpha=54.6 corresponds to a 0.9996 maximum cut
521  if (sigmaalpha > 54.6) sigmaalpha = 54.6;
522  sigma2E += sigmaalpha * sigmaalpha * zeta * zeta; // eV^2
523  } else { // Urban model
524  static const double alpha = 0.996;
525  double Ealpha = I / (1. - (alpha * Emax / (Emax + I))); // eV
526  double meanE32 = I * (Emax + I) / Emax * (Ealpha - I); // eV^2
527  sigma2E += fabs(stepSize_) * (Sigma1 * e1 * e1 + Sigma2 * e2 * e2 + Sigma3 * meanE32); // eV^2
528  }
529  }
530 
531  sigma2E *= 1.E-18; // eV -> GeV
532 
533  // update noise matrix, using linear error propagation from E to q/p
534  noise[6 * 7 + 6] += charge_*charge_/betaSquare / pow(mom, 4) * sigma2E;
535 }
536 
537 
539  const M1x3& direction, double momSquare, double betaSquare) const
540 {
541 
542  // MULTIPLE SCATTERING; calculate sigma^2
543  double sigma2 = 0;
544  assert(mscModelCode_ == 0 || mscModelCode_ == 1);
545  const double step = fabs(stepSize_);
546  const double step2 = step * step;
547  if (mscModelCode_ == 0) {// PANDA report PV/01-07 eq(43); linear in step length
548  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)
549 
550  } else if (mscModelCode_ == 1) { //Highland not linear in step length formula taken from PDG book 2011 edition
551  double stepOverRadLength = step / radiationLength_;
552  double logCor = (1 + 0.038 * log(stepOverRadLength));
553  sigma2 = 0.0136 * 0.0136 * charge_ * charge_ / (betaSquare * momSquare) * stepOverRadLength * logCor * logCor;
554  }
555  //assert(sigma2 >= 0.0);
556  sigma2 = (sigma2 > 0.0 ? sigma2 : 0.0);
557  //XXX debugOut << "MaterialEffects::noiseCoulomb the MSC variance is " << sigma2 << std::endl;
558 
559  M7x7 noiseAfter; // will hold the new MSC noise to cause by the current stepSize_ length
560  std::fill(noiseAfter.begin(), noiseAfter.end(), 0);
561 
562  const M1x3& a = direction; // as an abbreviation
563  // This calculates the MSC angular spread in the 7D global
564  // coordinate system. See PDG 2010, Sec. 27.3 for formulae.
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]; // Cov(x,a_y) = Cov(y,a_x)
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]; // Cov(z,a_x) = Cov(x,a_z)
581  noiseAfter[4 * 7 + 2] = noiseAfter[5 * 7 + 1]; // Cov(y,a_z) = Cov(z,a_y)
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]);
601 // debugOut << "new noise\n";
602 // RKTools::printDim(noiseAfter, 7,7);
603  for (unsigned int i = 0; i < 7 * 7; ++i) {
604  noise[i] += noiseAfter[i];
605  }
606 }
607 
608 
609 double MaterialEffects::dEdxBrems(double mom) const
610 {
611 
612  // Code ported from GEANT 3
613 
614  if (abs(pdg_) != 11) return 0; // only for electrons and positrons
615 
616 #if !defined(BETHE)
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;
619 #endif
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;
623 #endif
624 
625  double BCUT = 10000.; // energy up to which soft bremsstrahlung energy loss is calculated
626 
627  static const double THIGH = 100., CHIGH = 50.;
628  double dedxBrems = 0.;
629 
630  if (BCUT > 0.) {
631  double T, kc;
632 
633  if (BCUT >= mom) BCUT = mom; // confine BCUT to mom_
634 
635  // T=mom_, confined to THIGH
636  // kc=BCUT, confined to CHIGH ??
637  if (mom >= THIGH) {
638  T = THIGH;
639  if (BCUT >= THIGH) kc = CHIGH;
640  else kc = BCUT;
641  } else {
642  T = mom;
643  kc = BCUT;
644  }
645 
646  double E = T + me_; // total electron energy
647  if (BCUT > T) kc = T;
648 
649  double X = log(T / me_);
650  double Y = log(kc / (E * vl));
651 
652  double XX;
653  int K;
654  double S = 0., YY = 1.;
655 
656  for (unsigned int I = 1; I <= 2; ++I) {
657  XX = 1.;
658  for (unsigned int J = 1; J <= 6; ++J) {
659  K = 6 * I + J - 6;
660  S = S + C[K] * XX * YY;
661  XX = XX * X;
662  }
663  YY = YY * Y;
664  }
665 
666  for (unsigned int I = 3; I <= 6; ++I) {
667  XX = 1.;
668  for (unsigned int J = 1; J <= 6; ++J) {
669  K = 6 * I + J - 6;
670  if (Y <= 0.) S = S + C[K] * XX * YY;
671  else S = S + C[K + 24] * XX * YY;
672  XX = XX * X;
673  }
674  YY = YY * Y;
675  }
676 
677  double SS = 0.;
678  YY = 1.;
679 
680  for (unsigned int I = 1; I <= 2; ++I) {
681  XX = 1.;
682  for (unsigned int J = 1; J <= 5; ++J) {
683  K = 5 * I + J + 55;
684  SS = SS + C[K] * XX * YY;
685  XX = XX * X;
686  }
687  YY = YY * Y;
688  }
689 
690  for (unsigned int I = 3; I <= 5; ++I) {
691  XX = 1.;
692  for (unsigned int J = 1; J <= 5; ++J) {
693  K = 5 * I + J + 55;
694  if (Y <= 0.) SS = SS + C[K] * XX * YY;
695  else SS = SS + C[K + 15] * XX * YY;
696  XX = XX * X;
697  }
698  YY = YY * Y;
699  }
700 
701  S = S + matZ_ * SS;
702 
703  if (S > 0.) {
704  double CORR = 1.;
705 #if !defined(BETHE)
706  CORR = 1. / (1. + 0.805485E-10 * matDensity_ * matZ_ * E * E / (matA_ * kc * kc)); // MIGDAL correction factor
707 #endif
708 
709  // We use exp(beta * log(...) here because pow(..., beta) is
710  // REALLY slow and we don't need ultimate numerical precision
711  // for this approximation.
712  double FAC = matZ_ * (matZ_ + xi) * E * E / (E + me_);
713  if (beta == 1.) { // That is the #ifdef BETHE case
714  FAC *= kc * CORR / T;
715  } else {
716  FAC *= exp(beta * log(kc * CORR / T));
717  }
718  if (FAC <= 0.) return 0.;
719  dedxBrems = FAC * S;
720 
721 
722  if (mom >= THIGH) {
723  double RAT;
724  if (BCUT < THIGH) {
725  RAT = BCUT / mom;
726  S = (1. - 0.5 * RAT + 2.*RAT * RAT / 9.);
727  RAT = BCUT / T;
728  S = S / (1. - 0.5 * RAT + 2.*RAT * RAT / 9.);
729  } else {
730  RAT = BCUT / mom;
731  S = BCUT * (1. - 0.5 * RAT + 2.*RAT * RAT / 9.);
732  RAT = kc / T;
733  S = S / (kc * (1. - 0.5 * RAT + 2.*RAT * RAT / 9.));
734  }
735  dedxBrems = dedxBrems * S; // GeV barn
736  }
737 
738  dedxBrems = 0.60221367 * matDensity_ * dedxBrems / matA_; // energy loss dE/dx [GeV/cm]
739  }
740  }
741 
742  if (dedxBrems < 0.) dedxBrems = 0;
743 
744  double factor = 1.; // positron correction factor
745 
746  if (pdg_ == -11) {
747  static const double AA = 7522100., A1 = 0.415, A3 = 0.0021, A5 = 0.00054;
748 
749  double ETA = 0.;
750  if (matZ_ > 0.) {
751  double X = log(AA * mom / (matZ_ * matZ_));
752  if (X > -8.) {
753  if (X >= +9.) ETA = 1.;
754  else {
755  double W = A1 * X + A3 * pow(X, 3.) + A5 * pow(X, 5.);
756  ETA = 0.5 + atan(W) / M_PI;
757  }
758  }
759  }
760 
761  if (ETA < 0.0001) factor = 1.E-10;
762  else if (ETA > 0.9999) factor = 1.;
763  else {
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;
768  }
769  }
770 
771  return factor * dedxBrems; //always positive
772 }
773 
774 
775 void MaterialEffects::noiseBrems(M7x7& noise, double momSquare, double betaSquare) const
776 {
777 
778  // Code ported from GEANT 3 and simplified
779  // this formula assumes p >> m and therefore p^2 + m^2 = p^2
780  // the factor 1.44 is not in the original Behte-Heitler model.
781  // It seems to be some empirical correction copied over from some other project.
782 
783  if (abs(pdg_) != 11) return; // only for electrons and positrons
784 
785  double minusXOverLn2 = -1.442695 * fabs(stepSize_) / radiationLength_;
786  double sigma2 = 1.44*(pow(3., minusXOverLn2) - pow(4., minusXOverLn2)) / momSquare;
787  //XXX debugOut << "breams sigma: " << sigma2E << std::endl;
788  //assert(sigma2 >= 0.0);
789  sigma2 = (sigma2 > 0.0 ? sigma2 : 0.0);
790  noise[6 * 7 + 6] += charge_*charge_/betaSquare / (momSquare*momSquare) * sigma2;
791 
792 }
793 
794 
795 void MaterialEffects::setDebugLvl(unsigned int lvl) {
796  debugLvl_ = lvl;
797  if (materialInterface_ and debugLvl_ > 1)
799 }
800 
801 
803  pdg_ = pdg;
804  this->getParticleParameters();
805 
806  stepSize_ = 1;
807 
808  materialInterface_->initTrack(0, 0, 0, 1, 1, 1);
810 
811 
812  double minMom = 0.00001;
813  double maxMom = 10000;
814  int nSteps(10000);
815  double logStepSize = (log10(maxMom) - log10(minMom)) / (nSteps-1);
816 
817  TH1D hdEdxBethe("dEdxBethe", "dEdxBethe; log10(mom)", nSteps, log10(minMom), log10(maxMom));
818  TH1D hdEdxBrems("dEdxBrems", "dEdxBrems; log10(mom)", nSteps, log10(minMom), log10(maxMom));
819 
820  for (int i=0; i<nSteps; ++i) {
821  double mom = pow(10., log10(minMom) + i*logStepSize);
822  double E = hypot(mom, mass_);
823 
824  energyLossBrems_ = false;
825  energyLossBetheBloch_ = true;
826 
827  try {
828  hdEdxBethe.Fill(log10(mom), dEdx(E));
829  }
830  catch (...) {
831 
832  }
833 
834 
835  //debugOut<< "E = " << E << "; dEdx = " << dEdx(E) <<"\n";
836 
837  energyLossBrems_ = true;
838  energyLossBetheBloch_ = false;
839  try {
840  hdEdxBrems.Fill(log10(mom), dEdx(E));
841  }
842  catch (...) {
843 
844  }
845  }
846 
847  energyLossBrems_ = true;
848  energyLossBetheBloch_ = true;
849 
850  std::string Result;//string which will contain the result
851  std::stringstream convert; // stringstream used for the conversion
852  convert << pdg;//add the value of Number to the characters in the stream
853  Result = convert.str();//set Result to the content of the stream
854 
855  TFile outfile("dEdx_" + TString(Result) + ".root", "recreate");
856  outfile.cd();
857  hdEdxBethe.Write();
858  hdEdxBrems.Write();
859  outfile.Close();
860 }
861 
862 } /* End of namespace genfit */
863 
864 
Helper to store different limits on the stepsize for the RKTRackRep.
Definition: StepLimits.h:54
void stepper(const RKTrackRep *rep, M1x7 &state7, const double &mom, double &relMomLoss, const int &pdg, MaterialProperties &currentMaterial, StepLimits &limits, bool varField=true)
Returns maximum length so that a specified momentum loss will not be exceeded.
static MaterialEffects * instance_
double * end()
Definition: RKTools.h:46
double getLowestLimitSignedVal(double margin=1.E-3) const
Get the numerical value of the lowest limit, signed with stepSign_.
Definition: StepLimits.h:80
double RKPropagate(M1x7 &state7, M7x7 *jacobian, M1x3 &SA, double S, bool varField=true, bool calcOnlyLastRowOfJ=false) const
The actual Runge Kutta propagation.
Definition: RKTrackRep.cc:1272
double dEdxBrems(double mom) const
Returns dEdx.
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.
std::ostream debugOut
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...
Definition: StepLimits.h:89
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.
double * begin()
Definition: RKTools.h:45
virtual const char * what() const
Standard error message handling for exceptions. use like "std::cerr << e.what();".
Definition: Exception.cc:52
void setFatal(bool b=true)
Set fatal flag.
Definition: Exception.h:61
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) ...
Definition: Exception.h:48
void setMscModel(const std::string &modelName)
Select the multiple scattering model that will be used during track fit.
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) ...
Definition: RKTrackRep.h:71
Material properties needed e.g. for material effects calculation.
void setMaterialProperties(const double &density, const double &Z, const double &A, const double &radiationLength, const double &mEE)
std::ostream errorOut
void drawdEdx(int pdg=11)
void Print(const Option_t *="") const
static MaterialEffects * getInstance()
virtual void setDebugLvl(unsigned int lvl=1)
char getStepSign() const
Definition: StepLimits.h:84
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.
Definition: StepLimits.cc:65
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.