| File: | genfit2/code2/trackReps/src/MaterialEffects.cc |
| Warning: | line 437, column 12 Although the value stored to 'momLoss' is used in the enclosing expression, the value is never actually read from 'momLoss' |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
| 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 "MonopoleConstants.h" |
| 32 | #include <TMath.h> |
| 33 | |
| 34 | #include <TH1D.h> |
| 35 | #include <TFile.h> |
| 36 | |
| 37 | |
| 38 | namespace genfit { |
| 39 | |
| 40 | MaterialEffects* MaterialEffects::instance_ = nullptr; |
| 41 | |
| 42 | |
| 43 | MaterialEffects::MaterialEffects(): |
| 44 | noEffects_(false), |
| 45 | energyLossBetheBloch_(true), noiseBetheBloch_(true), |
| 46 | noiseCoulomb_(true), |
| 47 | energyLossBrems_(true), noiseBrems_(true), |
| 48 | ignoreBoundariesBetweenEqualMaterials_(true), |
| 49 | me_(0.510998910E-3), |
| 50 | stepSize_(0), |
| 51 | dEdx_(0), |
| 52 | E_(0), |
| 53 | matDensity_(0), |
| 54 | matZ_(0), |
| 55 | matA_(0), |
| 56 | radiationLength_(0), |
| 57 | mEE_(0), |
| 58 | pdg_(0), |
| 59 | charge_(0), |
| 60 | mag_charge_(0), |
| 61 | mass_(0), |
| 62 | mscModelCode_(0), |
| 63 | materialInterface_(nullptr), |
| 64 | debugLvl_(0) |
| 65 | { |
| 66 | } |
| 67 | |
| 68 | MaterialEffects::~MaterialEffects() |
| 69 | { |
| 70 | if (materialInterface_ != nullptr) delete materialInterface_; |
| 71 | } |
| 72 | |
| 73 | MaterialEffects* MaterialEffects::getInstance() |
| 74 | { |
| 75 | if (instance_ == nullptr) instance_ = new MaterialEffects(); |
| 76 | return instance_; |
| 77 | } |
| 78 | |
| 79 | void MaterialEffects::destruct() |
| 80 | { |
| 81 | if (instance_ != nullptr) { |
| 82 | delete instance_; |
| 83 | instance_ = nullptr; |
| 84 | } |
| 85 | } |
| 86 | |
| 87 | void MaterialEffects::init(AbsMaterialInterface* matIfc) |
| 88 | { |
| 89 | if (materialInterface_ != nullptr) { |
| 90 | std::string msg("MaterialEffects::initMaterialInterface(): Already initialized! "); |
| 91 | std::runtime_error err(msg); |
| 92 | } |
| 93 | materialInterface_ = matIfc; |
| 94 | } |
| 95 | |
| 96 | |
| 97 | |
| 98 | void MaterialEffects::setMscModel(const std::string& modelName) |
| 99 | { |
| 100 | if (modelName == std::string("GEANE")) { |
| 101 | mscModelCode_ = 0; |
| 102 | } else if (modelName == std::string("Highland")) { |
| 103 | mscModelCode_ = 1; |
| 104 | } else {// throw exception |
| 105 | std::string errorMsg = std::string("There is no MSC model called \"") + modelName + "\". Maybe it is not implemented or you misspelled the model name"; |
| 106 | Exception exc(errorMsg, __LINE__106, __FILE__"genfit2/code2/trackReps/src/MaterialEffects.cc"); |
| 107 | exc.setFatal(); |
| 108 | errorOut << exc.what(); |
| 109 | throw exc; |
| 110 | } |
| 111 | } |
| 112 | |
| 113 | |
| 114 | double MaterialEffects::effects(const std::vector<RKStep>& steps, |
| 115 | int materialsFXStart, |
| 116 | int materialsFXStop, |
| 117 | const double& mom, |
| 118 | const int& pdg, |
| 119 | M7x7* noise) |
| 120 | { |
| 121 | |
| 122 | if (debugLvl_ > 0) { |
| 123 | debugOut << " MaterialEffects::effects \n"; |
| 124 | } |
| 125 | |
| 126 | /*debugOut << "noEffects_ " << noEffects_ << "\n"; |
| 127 | debugOut << "energyLossBetheBloch_ " << energyLossBetheBloch_ << "\n"; |
| 128 | debugOut << "noiseBetheBloch_ " << noiseBetheBloch_ << "\n"; |
| 129 | debugOut << "noiseCoulomb_ " << noiseCoulomb_ << "\n"; |
| 130 | debugOut << "energyLossBrems_ " << energyLossBrems_ << "\n"; |
| 131 | debugOut << "noiseBrems_ " << noiseBrems_ << "\n";*/ |
| 132 | |
| 133 | |
| 134 | if (noEffects_) return 0.; |
| 135 | |
| 136 | if (materialInterface_ == nullptr) { |
| 137 | std::string msg("MaterialEffects hasn't been initialized with a correct AbsMaterialInterface pointer!"); |
| 138 | std::runtime_error err(msg); |
| 139 | throw err; |
| 140 | } |
| 141 | |
| 142 | bool doNoise(noise != nullptr); |
| 143 | |
| 144 | pdg_ = pdg; |
| 145 | getParticleParameters(); |
| 146 | |
| 147 | double momLoss = 0.; |
| 148 | |
| 149 | for ( std::vector<RKStep>::const_iterator it = steps.begin() + materialsFXStart; it != steps.begin() + materialsFXStop; ++it) { // loop over steps |
| 150 | |
| 151 | double realPath = it->matStep_.stepSize_; |
| 152 | if (fabs(realPath) < 1.E-8) { |
| 153 | // do material effects only if distance is not too small |
| 154 | continue; |
| 155 | } |
| 156 | |
| 157 | if (debugLvl_ > 0) { |
| 158 | debugOut << " calculate matFX "; |
| 159 | if (doNoise) |
| 160 | debugOut << "and noise"; |
| 161 | debugOut << " for "; |
| 162 | debugOut << "stepSize = " << it->matStep_.stepSize_ << "\t"; |
| 163 | it->matStep_.material_.Print(); |
| 164 | } |
| 165 | |
| 166 | double stepSign(1.); |
| 167 | if (realPath < 0) |
| 168 | stepSign = -1.; |
| 169 | realPath = fabs(realPath); |
| 170 | stepSize_ = realPath; |
| 171 | |
| 172 | |
| 173 | const Material& currentMaterial = it->matStep_.material_; |
| 174 | matDensity_ = currentMaterial.density; |
| 175 | matZ_ = currentMaterial.Z; |
| 176 | matA_ = currentMaterial.A; |
| 177 | radiationLength_ = currentMaterial.radiationLength; |
| 178 | mEE_ = currentMaterial.mEE; |
| 179 | |
| 180 | |
| 181 | if (matZ_ > 1.E-3) { // don't calculate energy loss for vacuum |
| 182 | |
| 183 | momLoss += momentumLoss(stepSign, mom - momLoss, false); |
| 184 | |
| 185 | if (doNoise){ |
| 186 | // get values for the "effective" energy of the RK step E_ |
| 187 | double p(0), gammaSquare(0), gamma(0), betaSquare(0); |
| 188 | this->getMomGammaBeta(E_, p, gammaSquare, gamma, betaSquare); |
| 189 | double pSquare = p*p; |
| 190 | |
| 191 | if (pdg_ == c_monopolePDGCode) { |
| 192 | charge_ = mag_charge_ * mom / hypot(mom, mass_); //effective charge for monopoles |
| 193 | } |
| 194 | |
| 195 | if (energyLossBetheBloch_ && noiseBetheBloch_) |
| 196 | this->noiseBetheBloch(*noise, p, betaSquare, gamma, gammaSquare); |
| 197 | |
| 198 | if (noiseCoulomb_) |
| 199 | this->noiseCoulomb(*noise, *((M1x3*) &it->state7_[3]), pSquare, betaSquare); |
| 200 | |
| 201 | if (energyLossBrems_ && noiseBrems_) |
| 202 | this->noiseBrems(*noise, pSquare, betaSquare); |
| 203 | } // end doNoise |
| 204 | |
| 205 | } |
| 206 | |
| 207 | } // end loop over steps |
| 208 | |
| 209 | if (momLoss >= mom) { |
| 210 | Exception exc("MaterialEffects::effects ==> momLoss >= momentum, aborting extrapolation!",__LINE__210,__FILE__"genfit2/code2/trackReps/src/MaterialEffects.cc"); |
| 211 | exc.setFatal(); |
| 212 | throw exc; |
| 213 | } |
| 214 | |
| 215 | return momLoss; |
| 216 | } |
| 217 | |
| 218 | |
| 219 | void MaterialEffects::stepper(const RKTrackRep* rep, |
| 220 | M1x7& state7, |
| 221 | const double& mom, // momentum |
| 222 | double& relMomLoss, // relative momloss for the step will be added |
| 223 | const int& pdg, |
| 224 | Material& currentMaterial, |
| 225 | StepLimits& limits, |
| 226 | bool varField) |
| 227 | { |
| 228 | |
| 229 | static const double maxRelMomLoss = .01; // maximum relative momentum loss allowed |
| 230 | static const double Pmin = 4.E-3; // minimum momentum for propagation [GeV] |
| 231 | static const double minStep = 1.E-4; // 1 µm |
| 232 | |
| 233 | // check momentum |
| 234 | if(mom < Pmin){ |
| 235 | std::ostringstream sstream; |
| 236 | sstream << "MaterialEffects::stepper ==> momentum too low: " << mom*1000. << " MeV"; |
| 237 | Exception exc(sstream.str(),__LINE__237,__FILE__"genfit2/code2/trackReps/src/MaterialEffects.cc"); |
| 238 | exc.setFatal(); |
| 239 | throw exc; |
| 240 | } |
| 241 | |
| 242 | // Trivial cases |
| 243 | if (noEffects_) |
| 244 | return; |
| 245 | |
| 246 | if (materialInterface_ == nullptr) { |
| 247 | std::string msg("MaterialEffects hasn't been initialized with a correct AbsMaterialInterface pointer!"); |
| 248 | std::runtime_error err(msg); |
| 249 | throw err; |
| 250 | } |
| 251 | |
| 252 | if (relMomLoss > maxRelMomLoss) { |
| 253 | limits.setLimit(stp_momLoss, 0); |
| 254 | return; |
| 255 | } |
| 256 | |
| 257 | |
| 258 | double sMax = limits.getLowestLimitSignedVal(); // signed |
| 259 | |
| 260 | if (fabs(sMax) < minStep) |
| 261 | return; |
| 262 | |
| 263 | |
| 264 | |
| 265 | pdg_ = pdg; |
| 266 | getParticleParameters(); |
| 267 | |
| 268 | |
| 269 | // make minStep |
| 270 | state7[0] += limits.getStepSign() * minStep * state7[3]; |
| 271 | state7[1] += limits.getStepSign() * minStep * state7[4]; |
| 272 | state7[2] += limits.getStepSign() * minStep * state7[5]; |
| 273 | |
| 274 | materialInterface_->initTrack(state7[0], state7[1], state7[2], |
| 275 | limits.getStepSign() * state7[3], limits.getStepSign() * state7[4], limits.getStepSign() * state7[5]); |
| 276 | |
| 277 | |
| 278 | currentMaterial = materialInterface_->getMaterialParameters(); |
| 279 | matDensity_ = currentMaterial.density; |
| 280 | matZ_ = currentMaterial.Z; |
| 281 | matA_ = currentMaterial.A; |
| 282 | radiationLength_ = currentMaterial.radiationLength; |
| 283 | mEE_ = currentMaterial.mEE; |
| 284 | |
| 285 | if (debugLvl_ > 0) { |
| 286 | debugOut << " currentMaterial "; currentMaterial.Print(); |
| 287 | } |
| 288 | |
| 289 | // limit due to momloss |
| 290 | double relMomLossPer_cm(0); |
| 291 | stepSize_ = 1.; // set stepsize for momLoss calculation |
| 292 | |
| 293 | if (matZ_ > 1.E-3) { // don't calculate energy loss for vacuum |
| 294 | relMomLossPer_cm = this->momentumLoss(limits.getStepSign(), mom, true) / mom; |
| 295 | } |
| 296 | |
| 297 | double maxStepMomLoss = fabs((maxRelMomLoss - fabs(relMomLoss)) / relMomLossPer_cm); // >= 0 |
| 298 | limits.setLimit(stp_momLoss, maxStepMomLoss); |
| 299 | |
| 300 | if (debugLvl_ > 0) { |
| 301 | debugOut << " momLoss exceeded after a step of " << maxStepMomLoss |
| 302 | << "; relMomLoss up to now = " << relMomLoss << "\n"; |
| 303 | } |
| 304 | |
| 305 | |
| 306 | // now look for boundaries |
| 307 | sMax = limits.getLowestLimitSignedVal(); |
| 308 | |
| 309 | stepSize_ = limits.getStepSign() * minStep; |
| 310 | M1x3 SA; |
| 311 | double boundaryStep(sMax); |
| 312 | |
| 313 | for (unsigned int i=0; i<100; ++i) { |
| 314 | if (debugLvl_ > 0) { |
| 315 | debugOut << " find next boundary\n"; |
| 316 | } |
| 317 | double step = materialInterface_->findNextBoundary(rep, state7, boundaryStep, varField); |
| 318 | |
| 319 | if (debugLvl_ > 0) { |
| 320 | if (step == 0) { |
| 321 | debugOut << " materialInterface_ returned a step of 0 \n"; |
| 322 | } |
| 323 | } |
| 324 | |
| 325 | stepSize_ += step; |
| 326 | boundaryStep -= step; |
| 327 | |
| 328 | if (debugLvl_ > 0) { |
| 329 | debugOut << " made a step of " << step << "\n"; |
| 330 | } |
| 331 | |
| 332 | if (! ignoreBoundariesBetweenEqualMaterials_) |
| 333 | break; |
| 334 | |
| 335 | if (fabs(stepSize_) >= fabs(sMax)) |
| 336 | break; |
| 337 | |
| 338 | // propagate with found step to boundary |
| 339 | rep->RKPropagate(state7, nullptr, SA, step, varField); |
| 340 | |
| 341 | // make minStep to cross boundary |
| 342 | state7[0] += limits.getStepSign() * minStep * state7[3]; |
| 343 | state7[1] += limits.getStepSign() * minStep * state7[4]; |
| 344 | state7[2] += limits.getStepSign() * minStep * state7[5]; |
| 345 | |
| 346 | materialInterface_->initTrack(state7[0], state7[1], state7[2], |
| 347 | limits.getStepSign() * state7[3], limits.getStepSign() * state7[4], limits.getStepSign() * state7[5]); |
| 348 | |
| 349 | Material materialAfter = materialInterface_->getMaterialParameters(); |
| 350 | |
| 351 | if (debugLvl_ > 0) { |
| 352 | debugOut << " material after step: "; materialAfter.Print(); |
| 353 | } |
| 354 | |
| 355 | if (materialAfter != currentMaterial) |
| 356 | break; |
| 357 | } |
| 358 | |
| 359 | limits.setLimit(stp_boundary, stepSize_); |
| 360 | |
| 361 | |
| 362 | relMomLoss += relMomLossPer_cm * limits.getLowestLimitVal(); |
| 363 | } |
| 364 | |
| 365 | |
| 366 | void MaterialEffects::getParticleParameters() |
| 367 | { |
| 368 | TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(pdg_); |
| 369 | charge_ = part->Charge() / 3.; // We only ever use the square |
| 370 | mass_ = part->Mass(); // GeV |
| 371 | } |
| 372 | |
| 373 | |
| 374 | void MaterialEffects::getMomGammaBeta(double Energy, |
| 375 | double& mom, double& gammaSquare, double& gamma, double& betaSquare) const { |
| 376 | |
| 377 | if (Energy <= mass_) { |
| 378 | Exception exc("MaterialEffects::getMomGammaBeta - Energy <= mass",__LINE__378,__FILE__"genfit2/code2/trackReps/src/MaterialEffects.cc"); |
| 379 | exc.setFatal(); |
| 380 | throw exc; |
| 381 | } |
| 382 | gamma = Energy/mass_; |
| 383 | gammaSquare = gamma*gamma; |
| 384 | betaSquare = 1.-1./gammaSquare; |
| 385 | mom = Energy*sqrt(betaSquare); |
| 386 | } |
| 387 | |
| 388 | |
| 389 | |
| 390 | //---- Energy-loss and Noise calculations ----------------------------------------- |
| 391 | |
| 392 | double MaterialEffects::momentumLoss(double stepSign, double mom, bool linear) |
| 393 | { |
| 394 | double E0 = hypot(mom, mass_); |
| 395 | double step = stepSize_*stepSign; // signed |
| 396 | |
| 397 | |
| 398 | // calc dEdx_, also needed in noiseBetheBloch! |
| 399 | // using fourth order Runge Kutta |
| 400 | //k1 = f(t0,y0) |
| 401 | //k2 = f(t0 + h/2, y0 + h/2 * k1) |
| 402 | //k3 = f(t0 + h/2, y0 + h/2 * k2) |
| 403 | //k4 = f(t0 + h, y0 + h * k3) |
| 404 | |
| 405 | // This means in our case: |
| 406 | //dEdx1 = dEdx(x0, E0) |
| 407 | //dEdx2 = dEdx(x0 + h/2, E1); E1 = E0 + h/2 * dEdx1 |
| 408 | //dEdx3 = dEdx(x0 + h/2, E2); E2 = E0 + h/2 * dEdx2 |
| 409 | //dEdx4 = dEdx(x0 + h, E3); E3 = E0 + h * dEdx3 |
| 410 | |
| 411 | double dEdx1 = dEdx(E0); // dEdx(x0,p0) |
| 412 | |
| 413 | if (linear) { |
| 414 | dEdx_ = dEdx1; |
| 415 | } |
| 416 | else { // RK4 |
| 417 | double E1 = E0 - dEdx1*step/2.; |
| 418 | double dEdx2 = dEdx(E1); // dEdx(x0 + h/2, E0 + h/2 * dEdx1) |
| 419 | |
| 420 | double E2 = E0 - dEdx2*step/2.; |
| 421 | double dEdx3 = dEdx(E2); // dEdx(x0 + h/2, E0 + h/2 * dEdx2) |
| 422 | |
| 423 | double E3 = E0 - dEdx3*step; |
| 424 | double dEdx4 = dEdx(E3); // dEdx(x0 + h, E0 + h * dEdx3) |
| 425 | |
| 426 | dEdx_ = (dEdx1 + 2.*dEdx2 + 2.*dEdx3 + dEdx4)/6.; |
| 427 | } |
| 428 | |
| 429 | E_ = E0 - dEdx_*step*0.5; |
| 430 | |
| 431 | double dE = step*dEdx_; // positive for positive stepSign |
| 432 | |
| 433 | double momLoss(0); |
| 434 | |
| 435 | if (E0 - dE <= mass_) { |
| 436 | // Step would stop particle (E_kin <= 0). |
| 437 | return momLoss = mom; |
Although the value stored to 'momLoss' is used in the enclosing expression, the value is never actually read from 'momLoss' | |
| 438 | } |
| 439 | else momLoss = mom - sqrt(pow(E0 - dE, 2) - mass_*mass_); // momLoss; positive for positive stepSign |
| 440 | |
| 441 | if (debugLvl_ > 0) { |
| 442 | debugOut << " MaterialEffects::momentumLoss: mom = " << mom << "; E0 = " << E0 |
| 443 | << "; dEdx = " << dEdx_ |
| 444 | << "; dE = " << dE << "; mass = " << mass_ << "\n"; |
| 445 | } |
| 446 | |
| 447 | //assert(momLoss * stepSign >= 0); |
| 448 | |
| 449 | return momLoss; |
| 450 | } |
| 451 | |
| 452 | |
| 453 | double MaterialEffects::dEdx(double Energy) { |
| 454 | |
| 455 | double mom(0), gammaSquare(0), gamma(0), betaSquare(0); |
| 456 | this->getMomGammaBeta(Energy, mom, gammaSquare, gamma, betaSquare); |
| 457 | if (pdg_ == c_monopolePDGCode) { // if TParticlePDG also had magnetic charge, life would have been easier. |
| 458 | charge_ = mag_charge_ * sqrt(betaSquare); //effective charge for monopoles |
| 459 | } |
| 460 | |
| 461 | double result(0); |
| 462 | |
| 463 | if (energyLossBetheBloch_) |
| 464 | result += dEdxBetheBloch(betaSquare, gamma, gammaSquare); |
| 465 | |
| 466 | if (energyLossBrems_) |
| 467 | result += dEdxBrems(mom); |
| 468 | |
| 469 | return result; |
| 470 | } |
| 471 | |
| 472 | |
| 473 | double MaterialEffects::dEdxBetheBloch(double betaSquare, double gamma, double gammaSquare) const |
| 474 | { |
| 475 | static const double betaGammaMin(0.05); |
| 476 | if (betaSquare*gammaSquare < betaGammaMin*betaGammaMin) { |
| 477 | Exception exc("MaterialEffects::dEdxBetheBloch ==> beta*gamma < 0.05, Bethe-Bloch implementation not valid anymore!",__LINE__477,__FILE__"genfit2/code2/trackReps/src/MaterialEffects.cc"); |
| 478 | exc.setFatal(); |
| 479 | throw exc; |
| 480 | } |
| 481 | |
| 482 | // calc dEdx_, also needed in noiseBetheBloch! |
| 483 | double result( 0.307075 * matZ_ / matA_ * matDensity_ / betaSquare * charge_ * charge_ ); |
| 484 | double massRatio( me_ / mass_ ); |
| 485 | double argument( gammaSquare * betaSquare * me_ * 1.E3 * 2. / ((1.E-6 * mEE_) * |
| 486 | sqrt(1. + 2. * gamma * massRatio + massRatio * massRatio)) ); |
| 487 | result *= log(argument) - betaSquare; // Bethe-Bloch [MeV/cm] |
| 488 | result *= 1.E-3; // in GeV/cm, hence 1.e-3 |
| 489 | if (result < 0.) { |
| 490 | result = 0; |
| 491 | } |
| 492 | |
| 493 | return result; |
| 494 | } |
| 495 | |
| 496 | |
| 497 | void MaterialEffects::noiseBetheBloch(M7x7& noise, double mom, double betaSquare, double gamma, double gammaSquare) const |
| 498 | { |
| 499 | // Code ported from GEANT 3 (erland.F) |
| 500 | |
| 501 | // ENERGY LOSS FLUCTUATIONS; calculate sigma^2(E); |
| 502 | double sigma2E ( 0. ); |
| 503 | double zeta ( 153.4E3 * charge_ * charge_ / betaSquare * matZ_ / matA_ * matDensity_ * fabs(stepSize_) ); // eV |
| 504 | double Emax ( 2.E9 * me_ * betaSquare * gammaSquare / (1. + 2.*gamma * me_ / mass_ + (me_ / mass_) * (me_ / mass_)) ); // eV |
| 505 | double kappa ( zeta / Emax ); |
| 506 | |
| 507 | if (kappa > 0.01) { // Vavilov-Gaussian regime |
| 508 | sigma2E += zeta * Emax * (1. - betaSquare / 2.); // eV^2 |
| 509 | } else { // Urban/Landau approximation |
| 510 | // calculate number of collisions Nc |
| 511 | double I = 16. * pow(matZ_, 0.9); // eV |
| 512 | double f2 = 0.; |
| 513 | if (matZ_ > 2.) f2 = 2. / matZ_; |
| 514 | double f1 = 1. - f2; |
| 515 | double e2 = 10.*matZ_ * matZ_; // eV |
| 516 | double e1 = pow((I / pow(e2, f2)), 1. / f1); // eV |
| 517 | |
| 518 | double mbbgg2 = 2.E9 * mass_ * betaSquare * gammaSquare; // eV |
| 519 | double Sigma1 = dEdx_ * 1.0E9 * f1 / e1 * (log(mbbgg2 / e1) - betaSquare) / (log(mbbgg2 / I) - betaSquare) * 0.6; // 1/cm |
| 520 | double Sigma2 = dEdx_ * 1.0E9 * f2 / e2 * (log(mbbgg2 / e2) - betaSquare) / (log(mbbgg2 / I) - betaSquare) * 0.6; // 1/cm |
| 521 | double Sigma3 = dEdx_ * 1.0E9 * Emax / (I * (Emax + I) * log((Emax + I) / I)) * 0.4; // 1/cm |
| 522 | |
| 523 | double Nc = (Sigma1 + Sigma2 + Sigma3) * fabs(stepSize_); |
| 524 | |
| 525 | if (Nc > 50.) { // truncated Landau distribution |
| 526 | double sigmaalpha = 15.76; |
| 527 | // calculate sigmaalpha (see GEANT3 manual W5013) |
| 528 | double RLAMED = -0.422784 - betaSquare - log(zeta / Emax); |
| 529 | double RLAMAX = 0.60715 + 1.1934 * RLAMED + (0.67794 + 0.052382 * RLAMED) * exp(0.94753 + 0.74442 * RLAMED); |
| 530 | // from lambda max to sigmaalpha=sigma (empirical polynomial) |
| 531 | if (RLAMAX <= 1010.) { |
| 532 | sigmaalpha = 1.975560 |
| 533 | + 9.898841e-02 * RLAMAX |
| 534 | - 2.828670e-04 * RLAMAX * RLAMAX |
| 535 | + 5.345406e-07 * pow(RLAMAX, 3.) |
| 536 | - 4.942035e-10 * pow(RLAMAX, 4.) |
| 537 | + 1.729807e-13 * pow(RLAMAX, 5.); |
| 538 | } else { sigmaalpha = 1.871887E+01 + 1.296254E-02 * RLAMAX; } |
| 539 | // alpha=54.6 corresponds to a 0.9996 maximum cut |
| 540 | if (sigmaalpha > 54.6) sigmaalpha = 54.6; |
| 541 | sigma2E += sigmaalpha * sigmaalpha * zeta * zeta; // eV^2 |
| 542 | } else { // Urban model |
| 543 | static const double alpha = 0.996; |
| 544 | double Ealpha = I / (1. - (alpha * Emax / (Emax + I))); // eV |
| 545 | double meanE32 = I * (Emax + I) / Emax * (Ealpha - I); // eV^2 |
| 546 | sigma2E += fabs(stepSize_) * (Sigma1 * e1 * e1 + Sigma2 * e2 * e2 + Sigma3 * meanE32); // eV^2 |
| 547 | } |
| 548 | } |
| 549 | |
| 550 | sigma2E *= 1.E-18; // eV -> GeV |
| 551 | |
| 552 | // update noise matrix, using linear error propagation from E to q/p |
| 553 | noise[6 * 7 + 6] += charge_*charge_/betaSquare / pow(mom, 4) * sigma2E; |
| 554 | } |
| 555 | |
| 556 | |
| 557 | void MaterialEffects::noiseCoulomb(M7x7& noise, |
| 558 | const M1x3& direction, double momSquare, double betaSquare) const |
| 559 | { |
| 560 | |
| 561 | // MULTIPLE SCATTERING; calculate sigma^2 |
| 562 | double sigma2 = 0; |
| 563 | assert(mscModelCode_ == 0 || mscModelCode_ == 1)(static_cast <bool> (mscModelCode_ == 0 || mscModelCode_ == 1) ? void (0) : __assert_fail ("mscModelCode_ == 0 || mscModelCode_ == 1" , "genfit2/code2/trackReps/src/MaterialEffects.cc", 563, __extension__ __PRETTY_FUNCTION__)); |
| 564 | const double step = fabs(stepSize_); |
| 565 | const double step2 = step * step; |
| 566 | if (mscModelCode_ == 0) {// PANDA report PV/01-07 eq(43); linear in step length |
| 567 | 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) |
| 568 | |
| 569 | } else if (mscModelCode_ == 1) { //Highland not linear in step length formula taken from PDG book 2011 edition |
| 570 | double stepOverRadLength = step / radiationLength_; |
| 571 | double logCor = (1 + 0.038 * log(stepOverRadLength)); |
| 572 | sigma2 = 0.0136 * 0.0136 * charge_ * charge_ / (betaSquare * momSquare) * stepOverRadLength * logCor * logCor; |
| 573 | } |
| 574 | //assert(sigma2 >= 0.0); |
| 575 | sigma2 = (sigma2 > 0.0 ? sigma2 : 0.0); |
| 576 | //XXX debugOut << "MaterialEffects::noiseCoulomb the MSC variance is " << sigma2 << std::endl; |
| 577 | |
| 578 | M7x7 noiseAfter; // will hold the new MSC noise to cause by the current stepSize_ length |
| 579 | std::fill(noiseAfter.begin(), noiseAfter.end(), 0); |
| 580 | |
| 581 | const M1x3& a = direction; // as an abbreviation |
| 582 | // This calculates the MSC angular spread in the 7D global |
| 583 | // coordinate system. See PDG 2010, Sec. 27.3 for formulae. |
| 584 | noiseAfter[0 * 7 + 0] = sigma2 * step2 / 3.0 * (1 - a[0]*a[0]); |
| 585 | noiseAfter[1 * 7 + 0] = -sigma2 * step2 / 3.0 * a[0]*a[1]; |
| 586 | noiseAfter[2 * 7 + 0] = -sigma2 * step2 / 3.0 * a[0]*a[2]; |
| 587 | noiseAfter[3 * 7 + 0] = sigma2 * step * 0.5 * (1 - a[0]*a[0]); |
| 588 | noiseAfter[4 * 7 + 0] = -sigma2 * step * 0.5 * a[0]*a[1]; |
| 589 | noiseAfter[5 * 7 + 0] = -sigma2 * step * 0.5 * a[0]*a[1]; |
| 590 | noiseAfter[0 * 7 + 1] = noiseAfter[1 * 7 + 0]; |
| 591 | noiseAfter[1 * 7 + 1] = sigma2 * step2 / 3.0 * (1 - a[1]*a[1]); |
| 592 | noiseAfter[2 * 7 + 1] = -sigma2 * step2 / 3.0 * a[1]*a[2]; |
| 593 | noiseAfter[3 * 7 + 1] = noiseAfter[4 * 7 + 0]; // Cov(x,a_y) = Cov(y,a_x) |
| 594 | noiseAfter[4 * 7 + 1] = sigma2 * step * 0.5 * (1 - a[1] * a[1]); |
| 595 | noiseAfter[5 * 7 + 1] = -sigma2 * step * 0.5 * a[1]*a[2]; |
| 596 | noiseAfter[0 * 7 + 2] = noiseAfter[2 * 7 + 0]; |
| 597 | noiseAfter[1 * 7 + 2] = noiseAfter[2 * 7 + 1]; |
| 598 | noiseAfter[2 * 7 + 2] = sigma2 * step2 / 3.0 * (1 - a[2]*a[2]); |
| 599 | noiseAfter[3 * 7 + 2] = noiseAfter[5 * 7 + 0]; // Cov(z,a_x) = Cov(x,a_z) |
| 600 | noiseAfter[4 * 7 + 2] = noiseAfter[5 * 7 + 1]; // Cov(y,a_z) = Cov(z,a_y) |
| 601 | noiseAfter[5 * 7 + 2] = sigma2 * step * 0.5 * (1 - a[2]*a[2]); |
| 602 | noiseAfter[0 * 7 + 3] = noiseAfter[3 * 7 + 0]; |
| 603 | noiseAfter[1 * 7 + 3] = noiseAfter[3 * 7 + 1]; |
| 604 | noiseAfter[2 * 7 + 3] = noiseAfter[3 * 7 + 2]; |
| 605 | noiseAfter[3 * 7 + 3] = sigma2 * (1 - a[0]*a[0]); |
| 606 | noiseAfter[4 * 7 + 3] = -sigma2 * a[0]*a[1]; |
| 607 | noiseAfter[5 * 7 + 3] = -sigma2 * a[0]*a[2]; |
| 608 | noiseAfter[0 * 7 + 4] = noiseAfter[4 * 7 + 0]; |
| 609 | noiseAfter[1 * 7 + 4] = noiseAfter[4 * 7 + 1]; |
| 610 | noiseAfter[2 * 7 + 4] = noiseAfter[4 * 7 + 2]; |
| 611 | noiseAfter[3 * 7 + 4] = noiseAfter[4 * 7 + 3]; |
| 612 | noiseAfter[4 * 7 + 4] = sigma2 * (1 - a[1]*a[1]); |
| 613 | noiseAfter[5 * 7 + 4] = -sigma2 * a[1]*a[2]; |
| 614 | noiseAfter[0 * 7 + 5] = noiseAfter[5 * 7 + 0]; |
| 615 | noiseAfter[1 * 7 + 5] = noiseAfter[5 * 7 + 1]; |
| 616 | noiseAfter[2 * 7 + 5] = noiseAfter[5 * 7 + 2]; |
| 617 | noiseAfter[3 * 7 + 5] = noiseAfter[5 * 7 + 3]; |
| 618 | noiseAfter[4 * 7 + 5] = noiseAfter[5 * 7 + 4]; |
| 619 | noiseAfter[5 * 7 + 5] = sigma2 * (1 - a[2]*a[2]); |
| 620 | // debugOut << "new noise\n"; |
| 621 | // RKTools::printDim(noiseAfter, 7,7); |
| 622 | for (unsigned int i = 0; i < 7 * 7; ++i) { |
| 623 | noise[i] += noiseAfter[i]; |
| 624 | } |
| 625 | } |
| 626 | |
| 627 | |
| 628 | double MaterialEffects::dEdxBrems(double mom) const |
| 629 | { |
| 630 | |
| 631 | // Code ported from GEANT 3 (gbrele.F) |
| 632 | |
| 633 | if (abs(pdg_) != 11) return 0; // only for electrons and positrons |
| 634 | |
| 635 | #if !defined(BETHE) |
| 636 | 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}; |
| 637 | static const double xi = 2.51, beta = 0.99, vl = 0.00004; |
| 638 | #endif |
| 639 | #if defined(BETHE) // no MIGDAL corrections |
| 640 | 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}; |
| 641 | static const double xi = 2.10, beta = 1.00, vl = 0.001; |
| 642 | #endif |
| 643 | |
| 644 | double BCUT = 10000.; // energy up to which soft bremsstrahlung energy loss is calculated |
| 645 | |
| 646 | static const double THIGH = 100., CHIGH = 50.; |
| 647 | double dedxBrems = 0.; |
| 648 | |
| 649 | if (BCUT > 0.) { |
| 650 | double T, kc; |
| 651 | |
| 652 | if (BCUT > mom) BCUT = mom; // confine BCUT to mom_ |
| 653 | |
| 654 | // T=mom_, confined to THIGH |
| 655 | // kc=BCUT, confined to CHIGH ?? |
| 656 | if (mom > THIGH) { |
| 657 | T = THIGH; |
| 658 | if (BCUT >= THIGH) |
| 659 | kc = CHIGH; |
| 660 | else |
| 661 | kc = BCUT; |
| 662 | } else { |
| 663 | T = mom; |
| 664 | kc = BCUT; |
| 665 | } |
| 666 | |
| 667 | double E = T + me_; // total electron energy |
| 668 | if (BCUT > T) |
| 669 | kc = T; |
| 670 | |
| 671 | double X = log(T / me_); |
| 672 | double Y = log(kc / (E * vl)); |
| 673 | |
| 674 | double XX; |
| 675 | int K; |
| 676 | double S = 0., YY = 1.; |
| 677 | |
| 678 | for (unsigned int I = 1; I <= 2; ++I) { |
| 679 | XX = 1.; |
| 680 | for (unsigned int J = 1; J <= 6; ++J) { |
| 681 | K = 6 * I + J - 6; |
| 682 | S += C[K] * XX * YY; |
| 683 | XX *= X; |
| 684 | } |
| 685 | YY *= Y; |
| 686 | } |
| 687 | |
| 688 | for (unsigned int I = 3; I <= 6; ++I) { |
| 689 | XX = 1.; |
| 690 | for (unsigned int J = 1; J <= 6; ++J) { |
| 691 | K = 6 * I + J - 6; |
| 692 | if (Y > 0.) |
| 693 | K += 24; |
| 694 | S += C[K] * XX * YY; |
| 695 | XX *= X; |
| 696 | } |
| 697 | YY *= Y; |
| 698 | } |
| 699 | |
| 700 | double SS = 0.; |
| 701 | YY = 1.; |
| 702 | |
| 703 | for (unsigned int I = 1; I <= 2; ++I) { |
| 704 | XX = 1.; |
| 705 | for (unsigned int J = 1; J <= 5; ++J) { |
| 706 | K = 5 * I + J + 55; |
| 707 | SS += C[K] * XX * YY; |
| 708 | XX *= X; |
| 709 | } |
| 710 | YY *= Y; |
| 711 | } |
| 712 | |
| 713 | for (unsigned int I = 3; I <= 5; ++I) { |
| 714 | XX = 1.; |
| 715 | for (unsigned int J = 1; J <= 5; ++J) { |
| 716 | K = 5 * I + J + 55; |
| 717 | if (Y > 0.) |
| 718 | K += 15; |
| 719 | SS += C[K] * XX * YY; |
| 720 | XX *= X; |
| 721 | } |
| 722 | YY *= Y; |
| 723 | } |
| 724 | |
| 725 | S += matZ_ * SS; |
| 726 | |
| 727 | if (S > 0.) { |
| 728 | double CORR = 1.; |
| 729 | #if !defined(BETHE) |
| 730 | CORR = 1. / (1. + 0.805485E-10 * matDensity_ * matZ_ * E * E / (matA_ * kc * kc)); // MIGDAL correction factor |
| 731 | #endif |
| 732 | |
| 733 | // We use exp(beta * log(...) here because pow(..., beta) is |
| 734 | // REALLY slow and we don't need ultimate numerical precision |
| 735 | // for this approximation. |
| 736 | double FAC = matZ_ * (matZ_ + xi) * E * E / (E + me_); |
| 737 | if (beta == 1.) // That is the #ifdef BETHE case |
| 738 | FAC *= kc * CORR / T; |
| 739 | else |
| 740 | FAC *= exp(beta * log(kc * CORR / T)); |
| 741 | if (FAC <= 0.) |
| 742 | return 0.; |
| 743 | dedxBrems = FAC * S; |
| 744 | |
| 745 | |
| 746 | if (mom >= THIGH) { |
| 747 | double RAT; |
| 748 | if (BCUT < THIGH) { |
| 749 | RAT = BCUT / mom; |
| 750 | S = (1. - 0.5 * RAT + 2.*RAT * RAT / 9.); |
| 751 | RAT = BCUT / T; |
| 752 | S /= 1. - 0.5 * RAT + 2.*RAT * RAT / 9.; |
| 753 | } else { |
| 754 | RAT = BCUT / mom; |
| 755 | S = BCUT * (1. - 0.5 * RAT + 2.*RAT * RAT / 9.); |
| 756 | RAT = kc / T; |
| 757 | S /= kc * (1. - 0.5 * RAT + 2.*RAT * RAT / 9.); |
| 758 | } |
| 759 | dedxBrems *= S; // GeV barn |
| 760 | } |
| 761 | |
| 762 | dedxBrems *= 0.60221367 * matDensity_ / matA_; // energy loss dE/dx [GeV/cm] |
| 763 | } |
| 764 | } |
| 765 | |
| 766 | if (dedxBrems < 0.) |
| 767 | dedxBrems = 0; |
| 768 | |
| 769 | double factor = 1.; // positron correction factor |
| 770 | |
| 771 | if (pdg_ == -11) { |
| 772 | static const double AA = 7522100., A1 = 0.415, A3 = 0.0021, A5 = 0.00054; |
| 773 | |
| 774 | double ETA = 0.; |
| 775 | if (matZ_ > 0.) { |
| 776 | double X = log(AA * mom / (matZ_ * matZ_)); |
| 777 | if (X > -8.) { |
| 778 | if (X >= +9.) ETA = 1.; |
| 779 | else { |
| 780 | double W = A1 * X + A3 * pow(X, 3.) + A5 * pow(X, 5.); |
| 781 | ETA = 0.5 + atan(W) / M_PI3.14159265358979323846; |
| 782 | } |
| 783 | } |
| 784 | } |
| 785 | |
| 786 | if (ETA < 0.0001) |
| 787 | factor = 1.E-10; |
| 788 | else if (ETA > 0.9999) |
| 789 | factor = 1.; |
| 790 | else { |
| 791 | double E0 = BCUT / mom; |
| 792 | if (E0 > 1.) |
| 793 | E0 = 1.; |
| 794 | if (E0 < 1.E-8) |
| 795 | factor = 1.; |
| 796 | else |
| 797 | factor = ETA * (1. - pow(1. - E0, 1. / ETA)) / E0; |
| 798 | } |
| 799 | } |
| 800 | |
| 801 | return factor * dedxBrems; //always positive |
| 802 | } |
| 803 | |
| 804 | |
| 805 | void MaterialEffects::noiseBrems(M7x7& noise, double momSquare, double betaSquare) const |
| 806 | { |
| 807 | // Code ported from GEANT 3 (erland.F) and simplified |
| 808 | // E \approx p is assumed. |
| 809 | // the factor 1.44 is not in the original Bethe-Heitler model. |
| 810 | // It seems to be some empirical correction copied over from some other project. |
| 811 | |
| 812 | if (abs(pdg_) != 11) return; // only for electrons and positrons |
| 813 | |
| 814 | double minusXOverLn2 = -1.442695 * fabs(stepSize_) / radiationLength_; |
| 815 | double sigma2E = 1.44*(pow(3., minusXOverLn2) - pow(4., minusXOverLn2)) * momSquare; |
| 816 | sigma2E = std::max(sigma2E, 0.0); // must be positive |
| 817 | |
| 818 | // update noise matrix, using linear error propagation from E to q/p |
| 819 | noise[6 * 7 + 6] += charge_*charge_/betaSquare / pow(momSquare, 2) * sigma2E; |
| 820 | } |
| 821 | |
| 822 | |
| 823 | void MaterialEffects::setDebugLvl(unsigned int lvl) { |
| 824 | debugLvl_ = lvl; |
| 825 | if (materialInterface_ and debugLvl_ > 1) |
| 826 | materialInterface_->setDebugLvl(debugLvl_-1); |
| 827 | } |
| 828 | |
| 829 | |
| 830 | void MaterialEffects::drawdEdx(int pdg) { |
| 831 | pdg_ = pdg; |
| 832 | this->getParticleParameters(); |
| 833 | |
| 834 | stepSize_ = 1; |
| 835 | |
| 836 | materialInterface_->initTrack(0, 0, 0, 1, 1, 1); |
| 837 | auto currentMaterial = materialInterface_->getMaterialParameters(); |
| 838 | matDensity_ = currentMaterial.density; |
| 839 | matZ_ = currentMaterial.Z; |
| 840 | matA_ = currentMaterial.A; |
| 841 | radiationLength_ = currentMaterial.radiationLength; |
| 842 | mEE_ = currentMaterial.mEE; |
| 843 | |
| 844 | double minMom = 0.00001; |
| 845 | double maxMom = 10000; |
| 846 | int nSteps(10000); |
| 847 | double logStepSize = (log10(maxMom) - log10(minMom)) / (nSteps-1); |
| 848 | |
| 849 | TH1D hdEdxBethe("dEdxBethe", "dEdxBethe; log10(mom)", nSteps, log10(minMom), log10(maxMom)); |
| 850 | TH1D hdEdxBrems("dEdxBrems", "dEdxBrems; log10(mom)", nSteps, log10(minMom), log10(maxMom)); |
| 851 | |
| 852 | for (int i=0; i<nSteps; ++i) { |
| 853 | double mom = pow(10., log10(minMom) + i*logStepSize); |
| 854 | double E = hypot(mom, mass_); |
| 855 | if (pdg_ == c_monopolePDGCode) { |
| 856 | charge_ = mag_charge_ * mom / E; //effective charge for monopoles |
| 857 | } |
| 858 | |
| 859 | energyLossBrems_ = false; |
| 860 | energyLossBetheBloch_ = true; |
| 861 | |
| 862 | try { |
| 863 | hdEdxBethe.Fill(log10(mom), dEdx(E)); |
| 864 | } |
| 865 | catch (...) { |
| 866 | |
| 867 | } |
| 868 | |
| 869 | |
| 870 | //debugOut<< "E = " << E << "; dEdx = " << dEdx(E) <<"\n"; |
| 871 | |
| 872 | energyLossBrems_ = true; |
| 873 | energyLossBetheBloch_ = false; |
| 874 | try { |
| 875 | hdEdxBrems.Fill(log10(mom), dEdx(E)); |
| 876 | } |
| 877 | catch (...) { |
| 878 | |
| 879 | } |
| 880 | } |
| 881 | |
| 882 | energyLossBrems_ = true; |
| 883 | energyLossBetheBloch_ = true; |
| 884 | |
| 885 | std::string Result;//string which will contain the result |
| 886 | std::stringstream convert; // stringstream used for the conversion |
| 887 | convert << pdg;//add the value of Number to the characters in the stream |
| 888 | Result = convert.str();//set Result to the content of the stream |
| 889 | |
| 890 | TFile outfile("dEdx_" + TString(Result) + ".root", "recreate"); |
| 891 | outfile.cd(); |
| 892 | hdEdxBethe.Write(); |
| 893 | hdEdxBrems.Write(); |
| 894 | outfile.Close(); |
| 895 | } |
| 896 | |
| 897 | } /* End of namespace genfit */ |
| 898 | |
| 899 |