Skip to content

Commit 41cc1e1

Browse files
shahor02chiarazampolli
shahor02
authored and
chiarazampolli
committed
Use iterative e.loss correction if kin.E loss fraction exceeds 2.5%
To account for the variation of the dedx when passing provided xrho material.
1 parent 868e77d commit 41cc1e1

6 files changed

Lines changed: 248 additions & 85 deletions

File tree

‎DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h‎

Lines changed: 28 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -115,8 +115,12 @@ GPUconstexpr() int CovarMap[kNParams][kNParams] = {{0, 1, 3, 6, 10},
115115
GPUconstexpr() int DiagMap[kNParams] = {0, 2, 5, 9, 14};
116116

117117
constexpr float HugeF = o2::constants::math::VeryBig;
118-
constexpr float MaxPT = 100000.; // do not allow pTs exceeding this value (to avoid NANs)
119-
constexpr float MinPTInv = 1. / MaxPT; // do not allow q/pTs less this value (to avoid NANs)
118+
constexpr float MaxPT = 100000.; // do not allow pTs exceeding this value (to avoid NANs)
119+
constexpr float MinPTInv = 1. / MaxPT; // do not allow q/pTs less this value (to avoid NANs)
120+
constexpr float ELoss2EKinThreshInv = 1. / 0.025; // do not allow E.Loss correction step with dE/Ekin above the inverse of this value
121+
constexpr int MaxELossIter = 50; // max number of iteration for the ELoss to account for BB dependence on beta*gamma
122+
// uncomment this to enable correction for BB dependence on beta*gamma via BB derivative
123+
// #define _BB_NONCONST_CORR_
120124

121125
template <typename value_T = float>
122126
class TrackParametrization
@@ -193,6 +197,11 @@ class TrackParametrization
193197
GPUd() value_t getPInv() const;
194198
GPUd() value_t getP() const;
195199
GPUd() value_t getPt() const;
200+
GPUd() value_t getE2() const;
201+
GPUd() value_t getE() const;
202+
GPUd() static inline value_t getdEdxBB(value_t betagamma) { return BetheBlochSolid(betagamma); }
203+
GPUd() static inline value_t getdEdxBBOpt(value_t betagamma) { return BetheBlochSolidOpt(betagamma); }
204+
GPUd() static inline value_t getBetheBlochSolidDerivativeApprox(value_T dedx, value_T bg) { return BetheBlochSolidDerivative(dedx, bg); }
196205

197206
GPUd() value_t getTheta() const;
198207
GPUd() value_t getEta() const;
@@ -210,7 +219,7 @@ class TrackParametrization
210219
GPUd() math_utils::Point3D<value_t> getXYZGloAt(value_t xk, value_t b, bool& ok) const;
211220

212221
// parameters manipulation
213-
GPUd() bool correctForELoss(value_t xrho, bool anglecorr = false, value_t dedx = kCalcdEdxAuto);
222+
GPUd() bool correctForELoss(value_t xrho, bool anglecorr = false);
214223
GPUd() bool rotateParam(value_t alpha);
215224
GPUd() bool propagateParamTo(value_t xk, value_t b);
216225
GPUd() bool propagateParamTo(value_t xk, const dim3_t& b);
@@ -585,6 +594,22 @@ GPUdi() auto TrackParametrization<value_T>::getP() const -> value_t
585594
return 1.f / getPInv(); // getPInv is already protected against being 0
586595
}
587596

597+
//____________________________________________________________
598+
template <typename value_T>
599+
GPUdi() auto TrackParametrization<value_T>::getE2() const -> value_t
600+
{
601+
// return the track energy^2
602+
return getP2() + getPID().getMass2();
603+
}
604+
605+
//____________________________________________________________
606+
template <typename value_T>
607+
GPUdi() auto TrackParametrization<value_T>::getE() const -> value_t
608+
{
609+
// return the track energy
610+
return gpu::CAMath::Sqrt(getE2());
611+
}
612+
588613
//____________________________________________________________
589614
template <typename value_T>
590615
GPUdi() auto TrackParametrization<value_T>::getPt() const -> value_t

‎DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrizationWithError.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -111,7 +111,7 @@ class TrackParametrizationWithError : public TrackParametrization<value_T>
111111
template <typename T>
112112
GPUd() bool update(const BaseCluster<T>& p);
113113

114-
GPUd() bool correctForMaterial(value_t x2x0, value_t xrho, bool anglecorr = false, value_t dedx = kCalcdEdxAuto);
114+
GPUd() bool correctForMaterial(value_t x2x0, value_t xrho, bool anglecorr = false);
115115

116116
GPUd() void resetCovariance(value_t s2 = 0);
117117
GPUd() void checkCovariance();

‎DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackUtils.h‎

Lines changed: 77 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -33,8 +33,11 @@ namespace track
3333
{
3434
// helper function
3535
template <typename value_T = float>
36-
GPUd() value_T BetheBlochSolid(value_T bg, value_T rho = 2.33, value_T kp1 = 0.20, value_T kp2 = 3.00, value_T meanI = 173e-9,
37-
value_T meanZA = 0.49848);
36+
GPUd() value_T BetheBlochSolid(value_T bg, value_T rho = 2.33, value_T kp1 = 0.20, value_T kp2 = 3.00, value_T meanI = 173e-9, value_T meanZA = 0.49848);
37+
38+
template <typename value_T = float>
39+
GPUd() value_T BetheBlochSolidOpt(value_T bg);
40+
3841
template <typename value_T = float>
3942
GPUd() void g3helx3(value_T qfield, value_T step, gpu::gpustd::array<value_T, 7>& vect);
4043

@@ -121,11 +124,11 @@ GPUd() value_T BetheBlochSolid(value_T bg, value_T rho, value_T kp1, value_T kp2
121124
static_assert(std::is_floating_point_v<value_T>);
122125
#endif
123126

124-
constexpr value_T mK = 0.307075e-3f; // [GeV*cm^2/g]
125-
constexpr value_T me = 0.511e-3f; // [GeV/c^2]
127+
constexpr value_T mK = 0.307075e-3; // [GeV*cm^2/g]
128+
constexpr value_T me = 0.511e-3; // [GeV/c^2]
126129
kp1 *= 2.303f;
127130
kp2 *= 2.303f;
128-
value_T bg2 = bg * bg;
131+
value_T bg2 = bg * bg, beta2 = bg2 / (1 + bg2);
129132
value_T maxT = 2.f * me * bg2; // neglecting the electron mass
130133

131134
//*** Density effect
@@ -138,7 +141,75 @@ GPUd() value_T BetheBlochSolid(value_T bg, value_T rho, value_T kp1, value_T kp2
138141
double r = (kp2 - x) / (kp2 - kp1);
139142
d2 = lhwI + x - 0.5f + (0.5f - lhwI - kp1) * r * r * r;
140143
}
141-
return mK * meanZA * (1 + bg2) / bg2 * (0.5f * gpu::CAMath::Log(2 * me * bg2 * maxT / (meanI * meanI)) - bg2 / (1 + bg2) - d2);
144+
auto dedx = mK * meanZA / beta2 * (0.5f * gpu::CAMath::Log(2 * me * bg2 * maxT / (meanI * meanI)) - beta2 - d2);
145+
return dedx > 0. ? dedx : 0.;
146+
}
147+
148+
//____________________________________________________
149+
template <typename value_T>
150+
GPUd() value_T BetheBlochSolidOpt(value_T bg)
151+
{
152+
//
153+
// This is the parameterization of the Bethe-Bloch formula inspired by Geant with hardcoded constants and better optimization
154+
//
155+
// bg - beta*gamma
156+
// rho - density [g/cm^3]
157+
// kp1 - density effect first junction point
158+
// kp2 - density effect second junction point
159+
// meanI - mean excitation energy [GeV]
160+
// meanZA - mean Z/A
161+
//
162+
// The default values for the kp* parameters are for silicon.
163+
// The returned value is in [GeV/(g/cm^2)].
164+
//
165+
constexpr value_T mK = 0.307075e-3; // [GeV*cm^2/g]
166+
constexpr value_T me = 0.511e-3; // [GeV/c^2]
167+
constexpr value_T rho = 2.33;
168+
constexpr value_T kp1 = 0.20 * 2.303;
169+
constexpr value_T kp2 = 3.00 * 2.303;
170+
constexpr value_T meanI = 173e-9;
171+
constexpr value_T meanZA = 0.49848;
172+
constexpr value_T lhwI = -1.7175226; // gpu::CAMath::Log(28.816 * 1e-9 * gpu::CAMath::Sqrt(rho * meanZA) / meanI);
173+
constexpr value_T log2muTomeanI = 8.6839805; // gpu::CAMath::Log( 2. * me / meanI);
174+
175+
value_T bg2 = bg * bg, beta2 = bg2 / (1 + bg2);
176+
value_T maxT = 2.f * me * bg2; // neglecting the electron mass
177+
178+
//*** Density effect
179+
value_T d2 = 0.;
180+
const value_T x = gpu::CAMath::Log(bg);
181+
if (x > kp2) {
182+
d2 = lhwI + x - 0.5f;
183+
} else if (x > kp1) {
184+
value_T r = (kp2 - x) / (kp2 - kp1);
185+
d2 = lhwI + x - 0.5 + (0.5 - lhwI - kp1) * r * r * r;
186+
}
187+
auto dedx = mK * meanZA / beta2 * (log2muTomeanI + x + x - beta2 - d2);
188+
return dedx > 0. ? dedx : 0.;
189+
}
190+
191+
//____________________________________________________
192+
template <typename value_T>
193+
GPUd() value_T inline BetheBlochSolidDerivative(value_T dedx, value_T bg)
194+
{
195+
//
196+
// This is approximate derivative of the BB over betagamm, NO check for the consistency of the provided dedx and bg is done
197+
// Charge 1 particle is assumed for the provied dedx. For charge > 1 particles dedx/q^2 should be provided and obtained value must be scaled by q^2
198+
// The call should be usually done as
199+
// auto dedx = BetheBlochSolidOpt(bg);
200+
// // if derivative needed
201+
// auto ddedx = BetheBlochSolidDerivative(dedx, bg, bg*bg)
202+
//
203+
// dedx - precalculate dedx for bg
204+
// bg - beta*gamma
205+
//
206+
constexpr value_T mK = 0.307075e-3; // [GeV*cm^2/g]
207+
constexpr value_T meanZA = 0.49848;
208+
auto bg2 = bg * bg;
209+
auto t1 = 1 + bg2;
210+
// auto derH = (mK * meanZA * (t1+bg2) - dedx*bg2)/(bg*t1);
211+
auto derH = (mK * meanZA * (t1 + 1. / bg2) - dedx) / (bg * t1);
212+
return derH + derH;
142213
}
143214

144215
} // namespace track

‎DataFormats/Reconstruction/src/ReconstructionDataFormatsLinkDef.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@
3030
#pragma link C++ class o2::track::TrackParFwd + ;
3131
#pragma link C++ class o2::track::PID + ;
3232
#pragma link C++ class o2::track::TrackLTIntegral + ;
33+
#pragma link C++ class std::vector < o2::track::TrackLTIntegral> + ;
3334

3435
#pragma link C++ class o2::track::TrackParCovFwd + ;
3536
#pragma link C++ class std::vector < o2::track::TrackParCovFwd> + ;

‎DataFormats/Reconstruction/src/TrackParametrization.cxx‎

Lines changed: 62 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -757,7 +757,7 @@ GPUd() bool TrackParametrization<value_T>::getXatLabR(value_t r, value_t& x, val
757757

758758
//______________________________________________
759759
template <typename value_T>
760-
GPUd() bool TrackParametrization<value_T>::correctForELoss(value_t xrho, bool anglecorr, value_t dedx)
760+
GPUd() bool TrackParametrization<value_T>::correctForELoss(value_t xrho, bool anglecorr)
761761
{
762762
//------------------------------------------------------------------
763763
// This function corrects the track parameters for the energy loss in crossed material.
@@ -767,41 +767,74 @@ GPUd() bool TrackParametrization<value_T>::correctForELoss(value_t xrho, bool an
767767
// "dedx" - mean enery loss (GeV/(g/cm^2), if <=kCalcdEdxAuto : calculate on the fly
768768
// "anglecorr" - switch for the angular correction
769769
//------------------------------------------------------------------
770-
constexpr value_t kMaxELossFrac = 0.3f; // max allowed fractional eloss
771770
constexpr value_t kMinP = 0.01f; // kill below this momentum
772771

773-
// Apply angle correction, if requested
774-
if (anglecorr) {
775-
value_t csp2 = (1.f - getSnp()) * (1.f + getSnp()); // cos(phi)^2
776-
value_t cst2I = (1.f + getTgl() * getTgl()); // 1/cos(lambda)^2
777-
value_t angle = gpu::CAMath::Sqrt(cst2I / (csp2));
778-
xrho *= angle;
779-
}
780-
value_t p = getP();
781-
value_t p2 = p * p;
782-
value_t e2 = p2 + getPID().getMass2();
783-
value_t beta2 = p2 / e2;
784-
785-
// Calculating the energy loss corrections************************
786-
if ((xrho != 0.f) && (beta2 < 1.f)) {
787-
if (dedx < kCalcdEdxAuto + constants::math::Almost1) { // request to calculate dedx on the fly
788-
dedx = BetheBlochSolid(p / getPID().getMass());
789-
if (mAbsCharge != 1) {
790-
dedx *= mAbsCharge * mAbsCharge;
791-
}
772+
auto m = getPID().getMass();
773+
if (m > 0 && xrho != 0.f) {
774+
// Apply angle correction, if requested
775+
if (anglecorr) {
776+
value_t csp2 = (1.f - getSnp()) * (1.f + getSnp()); // cos(phi)^2
777+
value_t cst2I = (1.f + getTgl() * getTgl()); // 1/cos(lambda)^2
778+
value_t angle = gpu::CAMath::Sqrt(cst2I / (csp2));
779+
xrho *= angle;
780+
}
781+
int charge2 = getAbsCharge() * getAbsCharge();
782+
value_t p = getP(), p0 = p, p2 = p * p, e2 = p2 + getPID().getMass2(), massInv = 1. / m, bg = p * massInv;
783+
value_t e = gpu::CAMath::Sqrt(e2), ekin = e - m, dedx1 = getdEdxBBOpt(bg), dedx = dedx1, dedxDer = 0.;
784+
if (charge2 != 1) {
785+
dedx *= charge2;
792786
}
793-
794787
value_t dE = dedx * xrho;
795-
value_t e = gpu::CAMath::Sqrt(e2);
796-
if (gpu::CAMath::Abs(dE) > kMaxELossFrac * e) {
797-
return false; // 30% energy loss is too much!
788+
int na = 1 + int(gpu::CAMath::Abs(dE) / ekin * ELoss2EKinThreshInv);
789+
if (na > MaxELossIter) {
790+
na = MaxELossIter;
791+
}
792+
if (na > 1) {
793+
dE /= na;
794+
xrho /= na;
795+
#ifdef _BB_NONCONST_CORR_
796+
dedxDer = getBetheBlochSolidDerivativeApprox(dedx1, bg); // require correction for non-constantness of dedx vs betagamma
797+
if (charge2 != 1) {
798+
dedxDer *= charge2;
799+
}
800+
#endif
801+
}
802+
while (na--) {
803+
#ifdef _BB_NONCONST_CORR_
804+
if (dedxDer != 0.) { // correction for non-constantness of dedx vs beta*gamma (in linear approximation): for a single step dE -> dE * [(exp(dedxDer) - 1)/dedxDer]
805+
if (xrho < 0) {
806+
dedxDer = -dedxDer; // E.loss ( -> positive derivative)
807+
}
808+
auto corrC = (gpu::CAMath::Exp(dedxDer) - 1.) / dedxDer;
809+
dE *= corrC;
810+
}
811+
#endif
812+
e += dE;
813+
if (e > m) { // stopped
814+
p = gpu::CAMath::Sqrt(e * e - getPID().getMass2());
815+
} else {
816+
return false;
817+
}
818+
if (na) {
819+
bg = p * massInv;
820+
dedx = getdEdxBBOpt(bg);
821+
#ifdef _BB_NONCONST_CORR_
822+
dedxDer = getBetheBlochSolidDerivativeApprox(dedx, bg);
823+
#endif
824+
if (charge2 != 1) {
825+
dedx *= charge2;
826+
#ifdef _BB_NONCONST_CORR_
827+
dedxDer *= charge2;
828+
#endif
829+
}
830+
dE = dedx * xrho;
831+
}
798832
}
799-
value_t eupd = e + dE;
800-
value_t pupd2 = eupd * eupd - getPID().getMass2();
801-
if (pupd2 < kMinP * kMinP) {
833+
834+
if (p < kMinP) {
802835
return false;
803836
}
804-
setQ2Pt(getQ2Pt() * p / gpu::CAMath::Sqrt(pupd2));
837+
setQ2Pt(getQ2Pt() * p0 / p);
805838
}
806839

807840
return true;

0 commit comments

Comments
 (0)