Dataformats/TPC: Add PID Response class (#11646)

* Dataformats/TPC: Add PID Response class

* Dataformats/TPC: PID Response class copyright headers corrected

* Dataformats/TPC:  Implement requested changes to PID response class
GPUTracking/Merger: Add PID assignment

* Dataformats/TPC: Add GPU definitions to PID response class

* Dataformats/TPC: Move PID response class functions to header for GPU

* Dataformats/TPC: Add condition to check dEdx for PID

* Dataformats/TPC: Add band crossing regions to PIDResponse

DataFormats/Detectors/TPC/CMakeLists.txt

@@ -62,7 +62,8 @@ o2_target_root_dictionary(
           include/DataFormatsTPC/LtrCalibData.h
           include/DataFormatsTPC/VDriftCorrFact.h
           include/DataFormatsTPC/CalibdEdxCorrection.h
-          include/DataFormatsTPC/BetheBlochAleph.h)
+          include/DataFormatsTPC/BetheBlochAleph.h
+          include/DataFormatsTPC/PIDResponse.h)
 
 o2_add_test(
   ClusterNative

DataFormats/Detectors/TPC/include/DataFormatsTPC/BetheBlochAleph.h

@@ -12,21 +12,22 @@
 #ifndef AliceO2_TPC_BETHEBLOCH_H_
 #define AliceO2_TPC_BETHEBLOCH_H_
 
-#include <cmath>
+#include "GPUCommonDef.h"
+#include "GPUCommonMath.h"
 
 namespace o2
 {
 namespace tpc
 {
 
 template <typename T>
-inline T BetheBlochAleph(T bg, T kp1, T kp2, T kp3, T kp4, T kp5)
+GPUdi() T BetheBlochAleph(T bg, T kp1, T kp2, T kp3, T kp4, T kp5)
 {
-  T beta = bg / std::sqrt(static_cast<T>(1.) + bg * bg);
+  T beta = bg / o2::gpu::GPUCommonMath::Sqrt(static_cast<T>(1.) + bg * bg);
 
-  T aa = std::pow(beta, kp4);
-  T bb = std::pow(static_cast<T>(1.) / bg, kp5);
-  bb = std::log(kp3 + bb);
+  T aa = o2::gpu::GPUCommonMath::Pow(beta, kp4);
+  T bb = o2::gpu::GPUCommonMath::Pow(static_cast<T>(1.) / bg, kp5);
+  bb = o2::gpu::GPUCommonMath::Log(kp3 + bb);
 
   return (kp2 - aa - bb) * kp1 / aa;
 }

DataFormats/Detectors/TPC/include/DataFormatsTPC/PIDResponse.h

@@ -0,0 +1,127 @@
+// 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.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+///
+/// @file   PIDResponse.h
+/// @author Tuba Gündem, [email protected]
+///
+
+#ifndef AliceO2_TPC_PIDResponse_H
+#define AliceO2_TPC_PIDResponse_H
+
+// o2 includes
+#include "GPUCommonDef.h"
+#include "GPUCommonRtypes.h"
+#include "GPUCommonMath.h"
+#include "ReconstructionDataFormats/PID.h"
+#include "DataFormatsTPC/PIDResponse.h"
+#include "DataFormatsTPC/BetheBlochAleph.h"
+#include "DataFormatsTPC/TrackTPC.h"
+
+namespace o2::tpc
+{
+class TrackTPC;
+
+/// \brief PID response class
+///
+/// This class is used to handle the TPC PID response.
+///
+
+class PIDResponse
+{
+ public:
+  /// default constructor
+  PIDResponse() CON_DEFAULT;
+
+  /// default destructor
+  ~PIDResponse() CON_DEFAULT;
+
+  /// setters
+  GPUd() void setBetheBlochParams(const float betheBlochParams[5]);
+  GPUd() void setMIP(float mip) { mMIP = mip; }
+  GPUd() void setChargeFactor(float chargeFactor) { mChargeFactor = chargeFactor; }
+
+  /// getters
+  GPUd() const float* getBetheBlochParams() const { return mBetheBlochParams; }
+  GPUd() float getMIP() const { return mMIP; }
+  GPUd() float getChargeFactor() const { return mChargeFactor; }
+
+  /// get expected signal of the track
+  GPUd() float getExpectedSignal(const TrackTPC& track, const o2::track::PID::ID id) const;
+
+  /// get most probable PID of the track
+  GPUd() o2::track::PID::ID getMostProbablePID(const TrackTPC& track, float PID_EKrangeMin, float PID_EKrangeMax, float PID_EPrangeMin, float PID_EPrangeMax) const;
+
+ private:
+  float mBetheBlochParams[5] = {0.19310481, 4.26696118, 0.00522579, 2.38124907, 0.98055396}; // BBAleph average fit parameters
+  float mMIP = 50.f;
+  float mChargeFactor = 2.299999952316284f;
+
+#ifndef GPUCA_ALIROOT_LIB
+  ClassDefNV(PIDResponse, 1);
+#endif
+};
+
+GPUd() void PIDResponse::setBetheBlochParams(const float betheBlochParams[5])
+{
+  for (int i = 0; i < 5; i++) {
+    mBetheBlochParams[i] = betheBlochParams[i];
+  }
+}
+
+GPUd() float PIDResponse::getExpectedSignal(const TrackTPC& track, const o2::track::PID::ID id) const
+{
+  const float bg = static_cast<float>(track.getP() / o2::track::pid_constants::sMasses[id]);
+  if (bg < 0.05) {
+    return -999.;
+  }
+  const float bethe = mMIP * o2::tpc::BetheBlochAleph(bg, mBetheBlochParams[0], mBetheBlochParams[1], mBetheBlochParams[2], mBetheBlochParams[3], mBetheBlochParams[4]) * o2::gpu::GPUCommonMath::Pow(static_cast<float>(o2::track::pid_constants::sCharges[id]), mChargeFactor);
+  return bethe >= 0. ? bethe : -999.;
+}
+
+// get most probable PID
+GPUd() o2::track::PID::ID PIDResponse::getMostProbablePID(const TrackTPC& track, float PID_EKrangeMin, float PID_EKrangeMax, float PID_EPrangeMin, float PID_EPrangeMax) const
+{
+  const float dEdx = track.getdEdx().dEdxTotTPC;
+
+  if (dEdx < 10.) {
+    return o2::track::PID::Pion;
+  }
+
+  auto id = o2::track::PID::Electron;
+  float distanceMin = o2::gpu::GPUCommonMath::Abs(dEdx - getExpectedSignal(track, id));
+
+  // calculate the distance to the expected dEdx signals
+  // start from Pion to exlude Muons
+  for (o2::track::PID::ID i = o2::track::PID::Pion; i < o2::track::PID::NIDs; i++) {
+    const float distance = o2::gpu::GPUCommonMath::Abs(dEdx - getExpectedSignal(track, i));
+    if (distance < distanceMin) {
+      id = i;
+      distanceMin = distance;
+    }
+  }
+
+  // override the electrons to the closest alternative PID in the bands crossing regions
+  if (id == o2::track::PID::Electron) {
+    const float p = track.getP();
+    if ((p > PID_EKrangeMin) && (p < PID_EKrangeMax)) {
+      id = o2::track::PID::Kaon;
+    } else if ((p > PID_EPrangeMin) && (p < PID_EPrangeMax)) {
+      id = o2::track::PID::Proton;
+    }
+  }
+
+  return id;
+}
+
+} // namespace o2::tpc
+
+#endif

DataFormats/Detectors/TPC/src/DataFormatsTPCLinkDef.h

@@ -69,5 +69,6 @@
 #pragma link C++ class o2::tpc::dcs::DataPoint < o2::tpc::dcs::HV::StackState> + ;
 #pragma link C++ class o2::tpc::dcs::DataPointVector < o2::tpc::dcs::HV::StackState> + ;
 #pragma link C++ class o2::tpc::dcs::Gas + ;
+#pragma link C++ class o2::tpc::PIDResponse + ;
 
 #endif

GPU/GPUTracking/Definitions/GPUSettingsList.h

@@ -56,6 +56,10 @@ AddOptionRTC(tubeMaxSize2, float, 2.5f * 2.5f, "", 0, "Square of max tube size (
 AddOptionRTC(clustersShiftTimebins, float, 0, "", 0, "Shift of TPC clusters (applied during CTF cluster decoding)")
 AddOptionRTC(clustersShiftTimebinsClusterizer, float, 0, "", 0, "Shift of TPC clusters (applied during CTF clusterization)")
 AddOptionRTC(defaultZOffsetOverR, float, 0.5210953f, "", 0, "Shift of TPC clusters (applied during CTF cluster decoding)")
+AddOptionRTC(PID_EKrangeMin, float, 0.38, "", 0, "min P of electron/K BB bands crossing")
+AddOptionRTC(PID_EKrangeMax, float, 0.48, "", 0, "max P of electron/K BB bands crossing")
+AddOptionRTC(PID_EPrangeMin, float, 0.75, "", 0, "min P of electron/p BB bands crossing")
+AddOptionRTC(PID_EPrangeMax, float, 0.85, "", 0, "max P of electron/p BB bands crossing")
 AddOptionRTC(maxTimeBinAboveThresholdIn1000Bin, unsigned short, 500, "", 0, "Except pad from cluster finding if total number of charges in a fragment is above this baseline (disable = 0)")
 AddOptionRTC(maxConsecTimeBinAboveThreshold, unsigned short, 200, "", 0, "Except pad from cluster finding if number of consecutive charges in a fragment is above this baseline (disable = 0)")
 AddOptionRTC(noisyPadSaturationThreshold, unsigned short, 700, "", 0, "Threshold where a timebin is considered saturated, disabling the noisy pad check for that pad")

GPU/GPUTracking/Merger/GPUTPCGMO2Output.cxx

@@ -17,6 +17,7 @@
 #include "GPUCommonAlgorithm.h"
 #include "DataFormatsTPC/TrackTPC.h"
 #include "DataFormatsTPC/Constants.h"
+#include "DataFormatsTPC/PIDResponse.h"
 #include "TPCFastTransform.h"
 #include "CorrectionMapsHelper.h"
 
@@ -134,12 +135,21 @@ GPUdii() void GPUTPCGMO2Output::Thread<GPUTPCGMO2Output::output>(int nBlocks, in
     if (merger.Param().par.dodEdx) {
       oTrack.setdEdx(tracksdEdx[i]);
     }
+
     oTrack.setOuterParam(o2::track::TrackParCov(
       outerPar.X, outerPar.alpha,
       {outerPar.P[0], outerPar.P[1], outerPar.P[2], outerPar.P[3], outerPar.P[4]},
       {outerPar.C[0], outerPar.C[1], outerPar.C[2], outerPar.C[3], outerPar.C[4], outerPar.C[5],
        outerPar.C[6], outerPar.C[7], outerPar.C[8], outerPar.C[9], outerPar.C[10], outerPar.C[11],
        outerPar.C[12], outerPar.C[13], outerPar.C[14]}));
+
+    if (merger.Param().par.dodEdx) {
+      PIDResponse pidResponse{};
+      const auto pid = pidResponse.getMostProbablePID(oTrack, merger.Param().rec.tpc.PID_EKrangeMin, merger.Param().rec.tpc.PID_EKrangeMax, merger.Param().rec.tpc.PID_EPrangeMin, merger.Param().rec.tpc.PID_EPrangeMax);
+      oTrack.setPID(pid);
+      oTrack.getParamOut().setPID(pid);
+    }
+
     unsigned int nOutCl = tmpData[i].x;
     unsigned int clBuff = tmpData[i].y;
     oTrack.setClusterRef(clBuff, nOutCl);