GENFIT  Rev:NoNumberAvailable
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 
genfit::MaterialEffects::energyLossBrems_
bool energyLossBrems_
Definition: MaterialEffects.h:175
genfit::RKTrackRep::RKPropagate
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
genfit::MaterialEffects::mass_
double mass_
Definition: MaterialEffects.h:195
genfit::MaterialEffects::noiseCoulomb
void noiseCoulomb(M7x7 &noise, const M1x3 &direction, double momSquare, double betaSquare) const
calculation of multiple scattering
Definition: MaterialEffects.cc:538
genfit::MaterialEffects::getParticleParameters
void getParticleParameters()
sets charge_, mass_
Definition: MaterialEffects.cc:350
genfit::MaterialEffects::dEdxBrems
double dEdxBrems(double mom) const
Returns dEdx.
Definition: MaterialEffects.cc:609
genfit::MaterialEffects::noiseBrems_
bool noiseBrems_
Definition: MaterialEffects.h:176
genfit::RKMatrix< 7, 7 >
genfit::MaterialProperties::setMaterialProperties
void setMaterialProperties(const double &density, const double &Z, const double &A, const double &radiationLength, const double &mEE)
Definition: MaterialProperties.cc:57
genfit::MaterialEffects::materialInterface_
AbsMaterialInterface * materialInterface_
depending on this number a specific msc model is chosen in the noiseCoulomb function.
Definition: MaterialEffects.h:199
genfit::MaterialEffects::destruct
static void destruct()
Definition: MaterialEffects.cc:77
genfit::AbsMaterialInterface
Abstract base class for geometry interfacing.
Definition: AbsMaterialInterface.h:39
genfit::StepLimits::getLowestLimitVal
double getLowestLimitVal(double margin=1.E-3) const
Get the unsigned numerical value of the lowest limit.
Definition: StepLimits.cc:65
genfit::MaterialEffects::noEffects_
bool noEffects_
Definition: MaterialEffects.h:170
genfit::MaterialEffects::dEdx_
double dEdx_
Definition: MaterialEffects.h:185
genfit::MaterialEffects::stepper
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.
Definition: MaterialEffects.cc:206
genfit::Exception
Exception class for error handling in GENFIT (provides storage for diagnostic information)
Definition: Exception.h:48
genfit
Defines for I/O streams used for error and debug printing.
Definition: AbsFinitePlane.cc:22
genfit::MaterialEffects
Stepper and energy loss/noise matrix calculation.
Definition: MaterialEffects.h:50
genfit::MaterialEffects::charge_
int charge_
Definition: MaterialEffects.h:194
genfit::MaterialEffects::setDebugLvl
void setDebugLvl(unsigned int lvl=1)
Definition: MaterialEffects.cc:795
genfit::RKTrackRep
AbsTrackRep with 5D track parameterization in plane coordinates: (q/p, u', v', u, v)
Definition: RKTrackRep.h:71
genfit::MaterialEffects::dEdxBetheBloch
double dEdxBetheBloch(double betaSquare, double gamma, double gammasquare) const
Uses Bethe Bloch formula to calculate dEdx.
Definition: MaterialEffects.cc:454
genfit::MaterialEffects::~MaterialEffects
virtual ~MaterialEffects()
Definition: MaterialEffects.cc:66
genfit::MaterialEffects::MaterialEffects
MaterialEffects()
Definition: MaterialEffects.cc:42
genfit::MaterialEffects::noiseBetheBloch
void noiseBetheBloch(M7x7 &noise, double mom, double betaSquare, double gamma, double gammaSquare) const
calculation of energy loss straggeling
Definition: MaterialEffects.cc:478
genfit::RKMatrix::end
double * end()
Definition: RKTools.h:46
genfit::MaterialEffects::noiseBetheBloch_
bool noiseBetheBloch_
Definition: MaterialEffects.h:173
genfit::Exception::setFatal
void setFatal(bool b=true)
Set fatal flag.
Definition: Exception.h:61
genfit::MaterialEffects::pdg_
int pdg_
Definition: MaterialEffects.h:193
genfit::Exception::what
virtual const char * what() const
Standard error message handling for exceptions. use like "std::cerr << e.what();".
Definition: Exception.cc:52
MaterialEffects.h
genfit::MaterialEffects::stepSize_
double stepSize_
Definition: MaterialEffects.h:182
genfit::MaterialEffects::dEdx
double dEdx(double Energy) const
Calculate dEdx for a given energy.
Definition: MaterialEffects.cc:437
genfit::MaterialEffects::ignoreBoundariesBetweenEqualMaterials_
bool ignoreBoundariesBetweenEqualMaterials_
Definition: MaterialEffects.h:178
genfit::AbsMaterialInterface::setDebugLvl
virtual void setDebugLvl(unsigned int lvl=1)
Definition: AbsMaterialInterface.h:72
genfit::MaterialEffects::radiationLength_
double radiationLength_
Definition: MaterialEffects.h:190
genfit::stp_momLoss
Definition: StepLimits.h:39
IO.h
genfit::MaterialEffects::instance_
static MaterialEffects * instance_
Definition: MaterialEffects.h:57
genfit::MaterialEffects::effects
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.
Definition: MaterialEffects.cc:112
genfit::MaterialEffects::matDensity_
double matDensity_
Definition: MaterialEffects.h:187
genfit::MaterialProperties::Print
void Print(const Option_t *="") const
Definition: MaterialProperties.cc:70
genfit::StepLimits::getStepSign
char getStepSign() const
Definition: StepLimits.h:84
genfit::StepLimits::setLimit
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
genfit::AbsMaterialInterface::getMaterialParameters
virtual void getMaterialParameters(double &density, double &Z, double &A, double &radiationLength, double &mEE)=0
Get material parameters in current material.
genfit::MaterialEffects::momentumLoss
double momentumLoss(double stepSign, double mom, bool linear)
Returns momentum loss.
Definition: MaterialEffects.cc:376
genfit::stp_boundary
Definition: StepLimits.h:44
genfit::MaterialEffects::debugLvl_
unsigned int debugLvl_
Definition: MaterialEffects.h:201
genfit::AbsMaterialInterface::initTrack
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.
genfit::MaterialEffects::matZ_
double matZ_
Definition: MaterialEffects.h:188
genfit::MaterialEffects::setMscModel
void setMscModel(const std::string &modelName)
Select the multiple scattering model that will be used during track fit.
Definition: MaterialEffects.cc:96
genfit::MaterialEffects::E_
double E_
Definition: MaterialEffects.h:186
genfit::MaterialEffects::getMomGammaBeta
void getMomGammaBeta(double Energy, double &mom, double &gammaSquare, double &gamma, double &betaSquare) const
Definition: MaterialEffects.cc:358
genfit::StepLimits
Helper to store different limits on the stepsize for the RKTRackRep.
Definition: StepLimits.h:54
genfit::errorOut
std::ostream errorOut
genfit::MaterialEffects::noiseBrems
void noiseBrems(M7x7 &noise, double momSquare, double betaSquare) const
calculation of energy loss straggeling
Definition: MaterialEffects.cc:775
genfit::MaterialEffects::me_
const double me_
Definition: MaterialEffects.h:180
genfit::MaterialEffects::noiseCoulomb_
bool noiseCoulomb_
Definition: MaterialEffects.h:174
genfit::MaterialEffects::getInstance
static MaterialEffects * getInstance()
Definition: MaterialEffects.cc:71
genfit::RKMatrix::begin
double * begin()
Definition: RKTools.h:45
genfit::MaterialEffects::matA_
double matA_
Definition: MaterialEffects.h:189
Exception.h
genfit::MaterialEffects::mEE_
double mEE_
Definition: MaterialEffects.h:191
genfit::MaterialEffects::init
void init(AbsMaterialInterface *matIfc)
set the material interface here. Material interface classes must be derived from AbsMaterialInterface...
Definition: MaterialEffects.cc:85
genfit::AbsMaterialInterface::findNextBoundary
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.
genfit::MaterialEffects::mscModelCode_
int mscModelCode_
Definition: MaterialEffects.h:197
genfit::StepLimits::getLowestLimitSignedVal
double getLowestLimitSignedVal(double margin=1.E-3) const
Get the numerical value of the lowest limit, signed with stepSign_.
Definition: StepLimits.h:80
genfit::MaterialEffects::drawdEdx
void drawdEdx(int pdg=11)
Definition: MaterialEffects.cc:802
genfit::MaterialProperties
Material properties needed e.g. for material effects calculation.
Definition: MaterialProperties.h:35
genfit::debugOut
std::ostream debugOut
genfit::MaterialEffects::energyLossBetheBloch_
bool energyLossBetheBloch_
Definition: MaterialEffects.h:172