Avoid recalculating BB over and over

Author: ktf

Common/Core/PID/TPCPIDResponse.h

@@ -90,12 +90,17 @@ class Response
   /// Gets the deviation to the expected signal
   template <typename TrackType>
   float GetSignalDelta(const TrackType& trk, const o2::track::PID::ID id) const;
+
   /// Gets relative dEdx resolution contribution due to relative pt resolution
   float GetRelativeResolutiondEdx(const float p, const float mass, const float charge, const float resol) const;
 
   void PrintAll() const;
 
  private:
+  /// Compute expected sigma given a pre-computed expected signal, avoiding a redundant Bethe-Bloch call.
+  template <typename TrackType>
+  float sigmaFromSignal(float expectedSignal, const long multTPC, const TrackType& track, const o2::track::PID::ID id) const;
+
   std::array<float, 5> mBetheBlochParams = {0.03209809958934784, 19.9768009185791, 2.5266601063857674e-16, 2.7212300300598145, 6.080920219421387};
   std::array<float, 2> mResolutionParamsDefault = {0.07, 0.0};
   std::vector<double> mResolutionParams = {5.43799e-7, 0.053044, 0.667584, 0.0142667, 0.00235175, 1.22482, 2.3501e-7, 0.031585};
@@ -135,13 +140,18 @@ inline float Response::GetExpectedSigmaAtMultiplicity(const long multTPC, const
   if (!track.hasTPC()) {
     return -999.f;
   }
-  float resolution = 0.;
+  return sigmaFromSignal(GetExpectedSignal(track, id), multTPC, track, id);
+}
+
+template <typename TrackType>
+inline float Response::sigmaFromSignal(float expectedSignal, const long multTPC, const TrackType& track, const o2::track::PID::ID id) const
+{
+  float resolution = 0.f;
   if (mUseDefaultResolutionParam) {
-    const float reso = GetExpectedSignal(track, id) * mResolutionParamsDefault[0] * (static_cast<float>(track.tpcNClsFound()) > 0 ? std::sqrt(1. + mResolutionParamsDefault[1] / static_cast<float>(track.tpcNClsFound())) : 1.f);
+    const float reso = expectedSignal * mResolutionParamsDefault[0] * (static_cast<float>(track.tpcNClsFound()) > 0 ? std::sqrt(1. + mResolutionParamsDefault[1] / static_cast<float>(track.tpcNClsFound())) : 1.f);
     reso >= 0.f ? resolution = reso : resolution = -999.f;
   } else {
-
-    const double ncl = nClNorm / track.tpcNClsFound(); //
+    const double ncl = nClNorm / track.tpcNClsFound();
     const double p = track.tpcInnerParam();
     const double mass = o2::track::pid_constants::sMasses[id];
     const double bg = p / mass;
@@ -160,16 +170,15 @@ inline float Response::GetExpectedSigmaAtMultiplicity(const long multTPC, const
 template <typename CollisionType, typename TrackType>
 inline float Response::GetNumberOfSigma(const CollisionType& collision, const TrackType& trk, const o2::track::PID::ID id) const
 {
-  if (GetExpectedSigma(collision, trk, id) < 0.) {
-    return -999.f;
-  }
-  if (GetExpectedSignal(trk, id) < 0.) {
+  const float signal = GetExpectedSignal(trk, id);
+  if (signal < 0.f) {
     return -999.f;
   }
-  if (!trk.hasTPC()) {
+  const float sigma = sigmaFromSignal(signal, collision.multTPC(), trk, id);
+  if (sigma < 0.f) {
     return -999.f;
   }
-  return ((trk.tpcSignal() - GetExpectedSignal(trk, id)) / GetExpectedSigma(collision, trk, id));
+  return (trk.tpcSignal() - signal) / sigma;
 }
 
 template <typename CollisionType, typename TrackType>
@@ -181,29 +190,26 @@ inline float Response::GetNumberOfSigmaMCTuned(const CollisionType& collision, c
 template <typename TrackType>
 inline float Response::GetNumberOfSigmaMCTunedAtMultiplicity(const long multTPC, const TrackType& trk, const o2::track::PID::ID id, float mcTunedTPCSignal) const
 {
-  if (GetExpectedSigmaAtMultiplicity(multTPC, trk, id) < 0.) {
+  const float signal = GetExpectedSignal(trk, id);
+  if (signal < 0.f) {
     return -999.f;
   }
-  if (GetExpectedSignal(trk, id) < 0.) {
+  const float sigma = sigmaFromSignal(signal, multTPC, trk, id);
+  if (sigma < 0.f) {
     return -999.f;
   }
-  if (!trk.hasTPC()) {
-    return -999.f;
-  }
-  return ((mcTunedTPCSignal - GetExpectedSignal(trk, id)) / GetExpectedSigmaAtMultiplicity(multTPC, trk, id));
+  return (mcTunedTPCSignal - signal) / sigma;
 }
 
 /// Gets the deviation between the actual signal and the expected signal
 template <typename TrackType>
 inline float Response::GetSignalDelta(const TrackType& trk, const o2::track::PID::ID id) const
 {
-  if (GetExpectedSignal(trk, id) < 0.) {
-    return -999.f;
-  }
-  if (!trk.hasTPC()) {
+  const float signal = GetExpectedSignal(trk, id);
+  if (signal < 0.f) {
     return -999.f;
   }
-  return (trk.tpcSignal() - GetExpectedSignal(trk, id));
+  return trk.tpcSignal() - signal;
 }
 
 //// Gets relative dEdx resolution contribution due relative pt resolution