Lit: Remove Bremsstrahlung energy loss for electrons and replace with existing...

Lit: Remove Bremsstrahlung energy loss for electrons and replace with existing Bethe Bloch for electrons #3019

reco/littrack/std/propagation/CbmLitMaterialEffectsImp.cxx

-/* 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);
   }