diff --git a/reco/littrack/std/propagation/CbmLitMaterialEffectsImp.cxx b/reco/littrack/std/propagation/CbmLitMaterialEffectsImp.cxx index 66fa626e751967e7e1d4da13d5f9b77704ec46d5..419e1983a2d7444f6f2b63006353f3a080d9ec41 100644 --- a/reco/littrack/std/propagation/CbmLitMaterialEffectsImp.cxx +++ b/reco/littrack/std/propagation/CbmLitMaterialEffectsImp.cxx @@ -1,26 +1,26 @@ -/* Copyright (C) 2008-2017 GSI/JINR-LIT, Darmstadt/Dubna +/* Copyright (C) 2008-2023 GSI/JINR-LIT, Darmstadt/Dubna SPDX-License-Identifier: GPL-3.0-only - Authors: Andrey Lebedev [committer] */ + Authors: Andrey Lebedev [committer], Martin Beyer */ /** CbmLitMaterialEffectsImp.cxx *@author A.Lebedev <andrey.lebedev@gsi.de> *@since 2008 **/ + #include "propagation/CbmLitMaterialEffectsImp.h" #include "base/CbmLitDefaultSettings.h" #include "data/CbmLitTrackParam.h" #include "propagation/CbmLitMaterialInfo.h" -#include "TDatabasePDG.h" -#include "TParticlePDG.h" +#include <TDatabasePDG.h> +#include <TParticlePDG.h> #include <cassert> +#include <cmath> #include <cstdlib> #include <iostream> -#include <cmath> - CbmLitMaterialEffectsImp::CbmLitMaterialEffectsImp() : fMass(0.105) , fDownstream(true) @@ -34,19 +34,21 @@ CbmLitMaterialEffectsImp::~CbmLitMaterialEffectsImp() {} LitStatus CbmLitMaterialEffectsImp::Update(CbmLitTrackParam* par, const CbmLitMaterialInfo* mat, int pdg, bool downstream) { - if (mat->GetLength() * mat->GetRho() < 1e-10) { return kLITSUCCESS; } + if (mat->GetLength() * mat->GetRho() < 1e-10) { + return kLITSUCCESS; + } fDownstream = downstream; TDatabasePDG* db = TDatabasePDG::Instance(); TParticlePDG* particle = db->GetParticle(pdg); - assert(particle != NULL); + assert(particle != nullptr); fMass = particle->Mass(); fIsElectron = (std::abs(pdg) == 11) ? true : false; fIsMuon = (std::abs(pdg) == 13) ? true : false; AddEnergyLoss(par, mat); - // AddThinScatter(par, mat); + // AddThinScatter(par, mat); AddThickScatter(par, mat); return kLITSUCCESS; @@ -54,37 +56,19 @@ LitStatus CbmLitMaterialEffectsImp::Update(CbmLitTrackParam* par, const CbmLitMa void CbmLitMaterialEffectsImp::AddEnergyLoss(CbmLitTrackParam* par, const CbmLitMaterialInfo* mat) const { - if (fIsElectron) { - litfloat radLength = mat->GetLength() / mat->GetRL(); - litfloat t; - - if (!fDownstream) { t = radLength; } - else { - t = -radLength; - } - - litfloat qp = par->GetQp(); - qp *= std::exp(-t); - par->SetQp(qp); - - litfloat cov = par->GetCovariance(18); - cov += CalcSigmaSqQpElectron(par, mat); - par->SetCovariance(18, cov); - } - else { - litfloat Eloss = EnergyLoss(par, mat); - par->SetQp(CalcQpAfterEloss(par->GetQp(), Eloss)); + litfloat Eloss = EnergyLoss(par, mat); + par->SetQp(CalcQpAfterEloss(par->GetQp(), Eloss)); - litfloat cov = par->GetCovariance(18); - cov += CalcSigmaSqQp(par, mat); - par->SetCovariance(18, cov); - } + litfloat cov = par->GetCovariance(18); + cov += CalcSigmaSqQp(par, mat); // FIXME: Full bethe bloch is used even if energy loss is BetheBlochElectron + par->SetCovariance(18, cov); } void CbmLitMaterialEffectsImp::AddThickScatter(CbmLitTrackParam* par, const CbmLitMaterialInfo* mat) const { - if (mat->GetLength() < 1e-10) { return; } - + if (mat->GetLength() < 1e-10) { + return; + } litfloat tx = par->GetTx(); litfloat ty = par->GetTy(); litfloat thickness = mat->GetLength(); //cm @@ -122,7 +106,9 @@ void CbmLitMaterialEffectsImp::AddThickScatter(CbmLitTrackParam* par, const CbmL void CbmLitMaterialEffectsImp::AddThinScatter(CbmLitTrackParam* par, const CbmLitMaterialInfo* mat) const { - if (mat->GetLength() < 1e-10) { return; } + if (mat->GetLength() < 1e-10) { + return; + } litfloat tx = par->GetTx(); litfloat ty = par->GetTy(); litfloat thetaSq = CalcThetaSq(par, mat); @@ -163,7 +149,7 @@ litfloat CbmLitMaterialEffectsImp::EnergyLoss(const CbmLitTrackParam* par, const litfloat CbmLitMaterialEffectsImp::dEdx(const CbmLitTrackParam* par, const CbmLitMaterialInfo* mat) const { - litfloat dedx = BetheBloch(par, mat); + litfloat dedx = fIsElectron ? BetheBlochElectron(par, mat) : BetheBloch(par, mat); // dedx += BetheHeitler(par, mat); // if (fIsMuon) dedx += PairProduction(par, mat); return dedx; @@ -171,10 +157,10 @@ litfloat CbmLitMaterialEffectsImp::dEdx(const CbmLitTrackParam* par, const CbmLi litfloat CbmLitMaterialEffectsImp::BetheBloch(const CbmLitTrackParam* par, const CbmLitMaterialInfo* mat) const { - litfloat K = 0.000307075; // GeV * g^-1 * cm^2 + litfloat K = 0.000307075; // GeV * mol^-1 * cm^2 litfloat z = (par->GetQp() > 0.) ? 1 : -1.; - litfloat Z = mat->GetZ(); - litfloat A = mat->GetA(); + litfloat Z = mat->GetZ(); // no unit + litfloat A = mat->GetA(); // g * mol^-1 litfloat M = fMass; litfloat p = std::abs(1. / par->GetQp()); //GeV @@ -204,10 +190,9 @@ litfloat CbmLitMaterialEffectsImp::BetheBloch(const CbmLitTrackParam* par, const litfloat CbmLitMaterialEffectsImp::BetheBlochElectron(const CbmLitTrackParam* par, const CbmLitMaterialInfo* mat) const { - litfloat K = 0.000307075; // GeV * g^-1 * cm^2 - //myf z = (par->GetQp() > 0.) ? 1 : -1.; - litfloat Z = mat->GetZ(); - litfloat A = mat->GetA(); + litfloat K = 0.000307075; // GeV * mol^-1 * cm^2 + litfloat Z = mat->GetZ(); // no unit + litfloat A = mat->GetA(); // g * mol^-1 litfloat me = 0.000511; // GeV; litfloat p = std::abs(1. / par->GetQp()); //GeV @@ -216,10 +201,13 @@ litfloat CbmLitMaterialEffectsImp::BetheBlochElectron(const CbmLitTrackParam* pa litfloat I = CalcI(Z) * 1e-9; // GeV - if (par->GetQp() > 0) { // electrons + // TODO: check this, since such a asymmetry is not seen in geant3 or geant4 + // different energy loss should only be relevant for low momentum < 100 MeV + // For now only use the 2. equation, since it matches geant3 and geant4 better + if (false) { // electrons return K * (Z / A) * (std::log(2 * me / I) + 1.5 * std::log(gamma) - 0.975); } - else { //positrons + else { // positrons return K * (Z / A) * (std::log(2 * me / I) + 2. * std::log(gamma) - 1.); } } @@ -230,10 +218,14 @@ litfloat CbmLitMaterialEffectsImp::CalcQpAfterEloss(litfloat qp, litfloat eloss) litfloat p = std::abs(1. / qp); litfloat E = std::sqrt(p * p + massSq); litfloat q = (qp > 0) ? 1. : -1.; - if (!fDownstream) { eloss *= -1.0; } // TODO check this + if (!fDownstream) { + eloss *= -1.0; + } // TODO check this litfloat Enew = E - eloss; litfloat pnew = (Enew > fMass) ? std::sqrt(Enew * Enew - massSq) : 0.; - if (pnew != 0) { return q / pnew; } + if (pnew != 0) { + return q / pnew; + } else { return 1e5; } @@ -287,7 +279,9 @@ litfloat CbmLitMaterialEffectsImp::CalcSigmaSqQpElectron(const CbmLitTrackParam* litfloat CbmLitMaterialEffectsImp::CalcI(litfloat Z) const { // mean excitation energy in eV - if (Z > 16.) { return 10 * Z; } + if (Z > 16.) { + return 10 * Z; + } else { return 16 * std::pow(Z, 0.9); }