PIDResponse.h

Go to the documentation of this file.

// 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.
 
 
#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;
 
 
class PIDResponse
{
 public:
  PIDResponse() = default;
 
  ~PIDResponse() = default;
 
  GPUd() void setBetheBlochParams(const float betheBlochParams[5]);
  GPUd() void setMIP(float mip) { mMIP = mip; }
  GPUd() void setChargeFactor(float chargeFactor) { mChargeFactor = chargeFactor; }
 
  GPUd() const float* getBetheBlochParams() const { return mBetheBlochParams; }
  GPUd() float getMIP() const { return mMIP; }
  GPUd() float getChargeFactor() const { return mChargeFactor; }
 
  GPUd() float getExpectedSignal(const TrackTPC& track, const o2::track::PID::ID id) const;
 
  GPUd() o2::track::PID::ID getMostProbablePID(const TrackTPC& track, float PID_EKrangeMin, float PID_EKrangeMax, float PID_EPrangeMin, float PID_EPrangeMax, float PID_EDrangeMin, float PID_EDrangeMax, float PID_ETrangeMin, float PID_ETrangeMax, char PID_useNsigma, float PID_sigma) 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;
 
  ClassDefNV(PIDResponse, 1);
};
 
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, float PID_EDrangeMin, float PID_EDrangeMax, float PID_ETrangeMin, float PID_ETrangeMax, char PID_useNsigma, float PID_sigma) const
{
  const float dEdx = track.getdEdx().dEdxTotTPC;
 
  if (dEdx < 10.) {
    return o2::track::PID::Pion;
  }
 
  auto id = o2::track::PID::Electron;
  float distanceMin = 0.;
  float dEdxExpected = getExpectedSignal(track, id);
  if (PID_useNsigma) {
    // using nSigma
    distanceMin = o2::gpu::GPUCommonMath::Abs((dEdx - dEdxExpected) / (PID_sigma * dEdxExpected));
  } else {
    // using absolute distance
    distanceMin = o2::gpu::GPUCommonMath::Abs(dEdx - dEdxExpected);
  }
 
  // 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++) {
    float distance = 0.;
    dEdxExpected = getExpectedSignal(track, i);
    if (PID_useNsigma) {
      // using nSigma
      distance = o2::gpu::GPUCommonMath::Abs((dEdx - dEdxExpected) / (PID_sigma * dEdxExpected));
    } else {
      // using absolute distance
      distance = o2::gpu::GPUCommonMath::Abs(dEdx - dEdxExpected);
    }
    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;
    } else if ((p > PID_EDrangeMin) && (p < PID_EDrangeMax)) {
      id = o2::track::PID::Deuteron;
    } else if ((p > PID_ETrangeMin) && (p < PID_ETrangeMax)) {
      id = o2::track::PID::Triton;
    }
  }
 
  return id;
}
 
} // namespace o2::tpc
 
#endif