From b40ae22b2e428cf693e5fd1f4bb74b1f2c646f0b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tuba=20G=C3=BCndem?= <48834043+tubagundem@users.noreply.github.com> Date: Mon, 7 Aug 2023 10:22:56 +0200 Subject: [PATCH] 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 | 3 +- .../include/DataFormatsTPC/BetheBlochAleph.h | 13 +- .../TPC/include/DataFormatsTPC/PIDResponse.h | 127 ++++++++++++++++++ .../Detectors/TPC/src/DataFormatsTPCLinkDef.h | 1 + GPU/GPUTracking/Definitions/GPUSettingsList.h | 4 + GPU/GPUTracking/Merger/GPUTPCGMO2Output.cxx | 10 ++ 6 files changed, 151 insertions(+), 7 deletions(-) create mode 100644 DataFormats/Detectors/TPC/include/DataFormatsTPC/PIDResponse.h diff --git a/DataFormats/Detectors/TPC/CMakeLists.txt b/DataFormats/Detectors/TPC/CMakeLists.txt index 28c0bc4672561..b2f9eb9e53e85 100644 --- a/DataFormats/Detectors/TPC/CMakeLists.txt +++ b/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 diff --git a/DataFormats/Detectors/TPC/include/DataFormatsTPC/BetheBlochAleph.h b/DataFormats/Detectors/TPC/include/DataFormatsTPC/BetheBlochAleph.h index 3fff89b44ba0f..e8fe7457f3091 100644 --- a/DataFormats/Detectors/TPC/include/DataFormatsTPC/BetheBlochAleph.h +++ b/DataFormats/Detectors/TPC/include/DataFormatsTPC/BetheBlochAleph.h @@ -12,7 +12,8 @@ #ifndef AliceO2_TPC_BETHEBLOCH_H_ #define AliceO2_TPC_BETHEBLOCH_H_ -#include +#include "GPUCommonDef.h" +#include "GPUCommonMath.h" namespace o2 { @@ -20,13 +21,13 @@ namespace tpc { template -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(1.) + bg * bg); + T beta = bg / o2::gpu::GPUCommonMath::Sqrt(static_cast(1.) + bg * bg); - T aa = std::pow(beta, kp4); - T bb = std::pow(static_cast(1.) / bg, kp5); - bb = std::log(kp3 + bb); + T aa = o2::gpu::GPUCommonMath::Pow(beta, kp4); + T bb = o2::gpu::GPUCommonMath::Pow(static_cast(1.) / bg, kp5); + bb = o2::gpu::GPUCommonMath::Log(kp3 + bb); return (kp2 - aa - bb) * kp1 / aa; } diff --git a/DataFormats/Detectors/TPC/include/DataFormatsTPC/PIDResponse.h b/DataFormats/Detectors/TPC/include/DataFormatsTPC/PIDResponse.h new file mode 100644 index 0000000000000..da1df159da923 --- /dev/null +++ b/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, tuba.gundem@cern.ch +/// + +#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(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(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 diff --git a/DataFormats/Detectors/TPC/src/DataFormatsTPCLinkDef.h b/DataFormats/Detectors/TPC/src/DataFormatsTPCLinkDef.h index 1b9329d3e8280..023892f01b92a 100644 --- a/DataFormats/Detectors/TPC/src/DataFormatsTPCLinkDef.h +++ b/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 diff --git a/GPU/GPUTracking/Definitions/GPUSettingsList.h b/GPU/GPUTracking/Definitions/GPUSettingsList.h index 9c85fc3585f7b..046cfe9c46f4b 100644 --- a/GPU/GPUTracking/Definitions/GPUSettingsList.h +++ b/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") diff --git a/GPU/GPUTracking/Merger/GPUTPCGMO2Output.cxx b/GPU/GPUTracking/Merger/GPUTPCGMO2Output.cxx index bc6b5baae8a4d..7fbbeafed3cf6 100644 --- a/GPU/GPUTracking/Merger/GPUTPCGMO2Output.cxx +++ b/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(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);