diff --git a/Common/Core/PID/TPCPIDResponse.h b/Common/Core/PID/TPCPIDResponse.h index 7831039f4d8..5120950fa5f 100644 --- a/Common/Core/PID/TPCPIDResponse.h +++ b/Common/Core/PID/TPCPIDResponse.h @@ -90,12 +90,17 @@ class Response /// Gets the deviation to the expected signal template 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 + float sigmaFromSignal(float expectedSignal, const long multTPC, const TrackType& track, const o2::track::PID::ID id) const; + std::array mBetheBlochParams = {0.03209809958934784, 19.9768009185791, 2.5266601063857674e-16, 2.7212300300598145, 6.080920219421387}; std::array mResolutionParamsDefault = {0.07, 0.0}; std::vector 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 +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(track.tpcNClsFound()) > 0 ? std::sqrt(1. + mResolutionParamsDefault[1] / static_cast(track.tpcNClsFound())) : 1.f); + const float reso = expectedSignal * mResolutionParamsDefault[0] * (static_cast(track.tpcNClsFound()) > 0 ? std::sqrt(1. + mResolutionParamsDefault[1] / static_cast(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 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 @@ -181,29 +190,26 @@ inline float Response::GetNumberOfSigmaMCTuned(const CollisionType& collision, c template 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 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