[Common] single-device TPC PID service for testing

Common/Core/PID/TPCPIDResponse.h

@@ -1,4 +1,4 @@
 // Copyright 2019-2020 CERN and copyright holders of ALICE O2.
 // See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
 // All rights not expressly granted are reserved.
 //
@@ -12,7 +12,7 @@
 /// \file   TPCPIDResponse.h
 /// \author Annalena Kalteyer <annalena.sophie.kalteyer@cern.ch>
 /// \author Christian Sonnabend
 /// \brief
 
 #ifndef COMMON_CORE_PID_TPCPIDRESPONSE_H_
 #define COMMON_CORE_PID_TPCPIDRESPONSE_H_
@@ -39,13 +39,13 @@
   ~Response() = default;
 
   /// Setter and Getter for the private parameters
   void SetBetheBlochParams(const std::array<float, 5>& betheBlochParams) { mBetheBlochParams = betheBlochParams; }
   void SetResolutionParamsDefault(const std::array<float, 2>& resolutionParamsDefault) { mResolutionParamsDefault = resolutionParamsDefault; }
   void SetResolutionParams(const std::vector<double>& resolutionParams) { mResolutionParams = resolutionParams; }
   void SetMIP(const float mip) { mMIP = mip; }
   void SetChargeFactor(const float chargeFactor) { mChargeFactor = chargeFactor; }
   void SetMultiplicityNormalization(const float multNormalization) { mMultNormalization = multNormalization; }
   void SetNClNormalization(const float nclnorm) { nClNorm = nclnorm; }
   void SetUseDefaultResolutionParam(const bool useDefault) { mUseDefaultResolutionParam = useDefault; }
   void SetParameters(const Response* response)
   {
@@ -74,12 +74,18 @@
   /// Gets the expected resolution of the track
   template <typename CollisionType, typename TrackType>
   float GetExpectedSigma(const CollisionType& collision, const TrackType& trk, const o2::track::PID::ID id) const;
+  /// Gets the expected resolution of the track with multTPC explicitly provided
+  template <typename TrackType>
+  float GetExpectedSigmaAtMultiplicity(const long multTPC, const TrackType& trk, const o2::track::PID::ID id) const;
   /// Gets the number of sigmas with respect the expected value
   template <typename CollisionType, typename TrackType>
   float GetNumberOfSigma(const CollisionType& collision, const TrackType& trk, const o2::track::PID::ID id) const;
   // Number of sigmas with respect to expected for MC, defining a tune-on-data signal value
   template <typename CollisionType, typename TrackType>
   float GetNumberOfSigmaMCTuned(const CollisionType& collision, const TrackType& trk, const o2::track::PID::ID id, float mcTunedTPCSignal) const;
+  // Number of sigmas with respect to expected for MC, defining a tune-on-data signal value, explicit multTPC
+  template <typename TrackType>
+  float GetNumberOfSigmaMCTunedAtMultiplicity(const long multTPC, const TrackType& trk, const o2::track::PID::ID id, float mcTunedTPCSignal) const;
   /// Gets the deviation to the expected signal
   template <typename TrackType>
   float GetSignalDelta(const TrackType& trk, const o2::track::PID::ID id) const;
@@ -116,6 +122,14 @@
 /// Gets the expected resolution of the measurement
 template <typename CollisionType, typename TrackType>
 inline float Response::GetExpectedSigma(const CollisionType& collision, const TrackType& track, const o2::track::PID::ID id) const
+{
+  // use multTPC (legacy behaviour) if multTPC not provided
+  return Response::GetExpectedSigmaAtMultiplicity(collision.multTPC(), track, id);
+}
+
+/// Gets the expected resolution of the measurement
+template <typename TrackType>
+inline float Response::GetExpectedSigmaAtMultiplicity(const long multTPC, const TrackType& track, const o2::track::PID::ID id) const
 {
   if (!track.hasTPC()) {
     return -999.f;
@@ -133,9 +147,9 @@
     const double dEdx = o2::tpc::BetheBlochAleph(static_cast<float>(bg), mBetheBlochParams[0], mBetheBlochParams[1], mBetheBlochParams[2], mBetheBlochParams[3], mBetheBlochParams[4]) * std::pow(static_cast<float>(o2::track::pid_constants::sCharges[id]), mChargeFactor);
     const double relReso = GetRelativeResolutiondEdx(p, mass, o2::track::pid_constants::sCharges[id], mResolutionParams[3]);
 
-    const std::vector<double> values{1.f / dEdx, track.tgl(), std::sqrt(ncl), relReso, track.signed1Pt(), collision.multTPC() / mMultNormalization};
+    const std::vector<double> values{1.f / dEdx, track.tgl(), std::sqrt(ncl), relReso, track.signed1Pt(), multTPC / mMultNormalization};
 
     const float reso = sqrt(pow(mResolutionParams[0], 2) * values[0] + pow(mResolutionParams[1], 2) * (values[2] * mResolutionParams[5]) * pow(values[0] / sqrt(1 + pow(values[1], 2)), mResolutionParams[2]) + values[2] * pow(values[3], 2) + pow(mResolutionParams[4] * values[4], 2) + pow(values[5] * mResolutionParams[6], 2) + pow(values[5] * (values[0] / sqrt(1 + pow(values[1], 2))) * mResolutionParams[7], 2)) * dEdx * mMIP;
     reso >= 0.f ? resolution = reso : resolution = -999.f;
   }
   return resolution;
@@ -160,7 +174,13 @@
 template <typename CollisionType, typename TrackType>
 inline float Response::GetNumberOfSigmaMCTuned(const CollisionType& collision, const TrackType& trk, const o2::track::PID::ID id, float mcTunedTPCSignal) const
 {
-  if (GetExpectedSigma(collision, trk, id) < 0.) {
+  return Response::GetNumberOfSigmaMCTunedAtMultiplicity(collision.multTPC(), trk, id, mcTunedTPCSignal);
+}
+
+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.) {
     return -999.f;
   }
   if (GetExpectedSignal(trk, id) < 0.) {
@@ -169,7 +189,7 @@
   if (!trk.hasTPC()) {
     return -999.f;
   }
-  return ((mcTunedTPCSignal - GetExpectedSignal(trk, id)) / GetExpectedSigma(collision, trk, id));
+  return ((mcTunedTPCSignal - GetExpectedSignal(trk, id)) / GetExpectedSigmaAtMultiplicity(multTPC, trk, id));
 }
 
 /// Gets the deviation between the actual signal and the expected signal

Common/TableProducer/PID/CMakeLists.txt

@@ -43,6 +43,11 @@ o2physics_add_dpl_workflow(pid-tof-full
 
 # TPC
 
+o2physics_add_dpl_workflow(pid-tpc-service
+                    SOURCES pidTPCService.cxx
+                    PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::MLCore O2Physics::AnalysisCCDB
+                    COMPONENT_NAME Analysis)
+
 o2physics_add_dpl_workflow(pid-tpc-base
                     SOURCES pidTPCBase.cxx
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::AnalysisCCDB

Common/TableProducer/PID/pidTPC.cxx

@@ -435,6 +435,11 @@ struct tpcPid {
     }
 
     float nSigma = -999.f;
+    int multTPC = 0;
+    if (trk.has_collision()) {
+      auto collision = collisions.rawIteratorAt(trk.collisionId());
+      multTPC = collision.multTPC();
+    }
     float bg = trk.tpcInnerParam() / o2::track::pid_constants::sMasses[pid]; // estimated beta-gamma for network cutoff
     if (useNetworkCorrection && speciesNetworkFlags[pid] && trk.has_collision() && bg > networkBetaGammaCutoff) {
 
@@ -457,7 +462,7 @@ struct tpcPid {
         LOGF(fatal, "Network output-dimensions incompatible!");
       }
     } else {
-      nSigma = response->GetNumberOfSigmaMCTuned(collisions.iteratorAt(trk.collisionId()), trk, pid, tpcSignal);
+      nSigma = response->GetNumberOfSigmaMCTunedAtMultiplicity(multTPC, trk, pid, tpcSignal);
     }
     if (flagFull)
       tableFull(expSigma, nSigma);