AliceO2Group · ktf · Aug 15, 2021 · Apr 1, 2021 · Apr 20, 2021 · May 1, 2021
@@ -0,0 +1,44 @@
+// 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   MID.h
+/// \author Vít Kučera
+/// \note   Based on RICH.h
+/// \brief  Set of tables for the ALICE3 MID information
+///
+
+#ifndef O2_ANALYSIS_ALICE3_MID_H_
+#define O2_ANALYSIS_ALICE3_MID_H_
+
+// O2 includes
+#include "Framework/AnalysisDataModel.h"
+
+namespace o2::aod
+{
+namespace alice3mid
+{
+DECLARE_SOA_INDEX_COLUMN(Collision, collision);    //!
+DECLARE_SOA_INDEX_COLUMN(Track, track);            //!
+DECLARE_SOA_COLUMN(MIDIsMuon, midIsMuon, uint8_t); //! FIXME: To be changed to bool once bool columns are groupable.
+} // namespace alice3mid
+
+DECLARE_SOA_TABLE(MIDs, "AOD", "MID", //!
+                  o2::soa::Index<>,
+                  alice3mid::CollisionId,
+                  alice3mid::TrackId,
+                  alice3mid::MIDIsMuon);
+
+using MID = MIDs::iterator;
+
+} // namespace o2::aod
+
+#endif // O2_ANALYSIS_ALICE3_MID_H_
@@ -29,7 +29,8 @@ enum Code {
   kDPlus = 411,
   kLambdaCPlus = 4122,
   kXiCPlus = 4232,
-  kJpsi = 443
+  kJpsi = 443,
+  kX3872 = 9920443
 };
 } // namespace pdg
 
@@ -99,7 +100,7 @@ auto pTBinsVec = std::vector<double>{pTBins, pTBins + npTBins + 1};
 
 // default values for the cuts
 constexpr double cuts[npTBins][nCutVars] = {{1.65, 2.15, 0.5, 100.},  /* 1 < pt < 5 */
-                                            {1.65, 2.15, 0.5, 100.}}; /* 5 < pt > 1000 */
+                                            {1.65, 2.15, 0.5, 100.}}; /* 5 < pt < 1000 */
 
 // row labels
 static const std::vector<std::string> pTBinLabels{};
@@ -383,7 +384,7 @@ static const std::vector<std::string> cutVarLabels = {"m", "pT p", "pT K", "pT P
 namespace hf_cuts_jpsi_toee
 {
 static constexpr int npTBins = 9;
-static constexpr int nCutVars = 4;
+static constexpr int nCutVars = 5;
 // default values for the pT bin edges (can be used to configure histogram axis)
 // offset by 1 from the bin numbers in cuts array
 constexpr double pTBins[npTBins + 1] = {
@@ -401,15 +402,15 @@ constexpr double pTBins[npTBins + 1] = {
 auto pTBins_v = std::vector<double>{pTBins, pTBins + npTBins + 1};
 
 // default values for the cuts
-constexpr double cuts[npTBins][nCutVars] = {{0.5, 0.2, 0.4, 1},  /* 0   < pT < 0.5 */
-                                            {0.5, 0.2, 0.4, 1},  /* 0.5 < pT < 1   */
-                                            {0.5, 0.2, 0.4, 1},  /* 1   < pT < 2   */
-                                            {0.5, 0.2, 0.4, 1},  /* 2   < pT < 3   */
-                                            {0.5, 0.2, 0.4, 1},  /* 3   < pT < 4   */
-                                            {0.5, 0.2, 0.4, 1},  /* 4   < pT < 5   */
-                                            {0.5, 0.2, 0.4, 1},  /* 5   < pT < 7   */
-                                            {0.5, 0.2, 0.4, 1},  /* 7   < pT < 10  */
-                                            {0.5, 0.2, 0.4, 1}}; /* 10  < pT < 15  */
+constexpr double cuts[npTBins][nCutVars] = {{0.5, 0.2, 0.4, 1, 1.},  /* 0   < pT < 0.5 */
+                                            {0.5, 0.2, 0.4, 1, 1.},  /* 0.5 < pT < 1   */
+                                            {0.5, 0.2, 0.4, 1, 1.},  /* 1   < pT < 2   */
+                                            {0.5, 0.2, 0.4, 1, 1.},  /* 2   < pT < 3   */
+                                            {0.5, 0.2, 0.4, 1, 1.},  /* 3   < pT < 4   */
+                                            {0.5, 0.2, 0.4, 1, 1.},  /* 4   < pT < 5   */
+                                            {0.5, 0.2, 0.4, 1, 1.},  /* 5   < pT < 7   */
+                                            {0.5, 0.2, 0.4, 1, 1.},  /* 7   < pT < 10  */
+                                            {0.5, 0.2, 0.4, 1, 1.}}; /* 10  < pT < 15  */
 
 // row labels
 static const std::vector<std::string> pTBinLabels = {
@@ -424,8 +425,54 @@ static const std::vector<std::string> pTBinLabels = {
   "pT bin 8"};
 
 // column labels
-static const std::vector<std::string> cutVarLabels = {"m", "DCA_xy", "DCA_z", "pT El"};
+static const std::vector<std::string> cutVarLabels = {"m", "DCA_xy", "DCA_z", "pT El", "chi2PCA"};
 } // namespace hf_cuts_jpsi_toee
+
+namespace hf_cuts_x_tojpsipipi
+{
+static constexpr int npTBins = 9;
+static constexpr int nCutVars = 7;
+// default values for the pT bin edges (can be used to configure histogram axis)
+// offset by 1 from the bin numbers in cuts array
+constexpr double pTBins[npTBins + 1] = {
+  0,
+  0.5,
+  1.0,
+  2.0,
+  3.0,
+  4.0,
+  5.0,
+  7.0,
+  10.0,
+  15.0,
+};
+auto pTBins_v = std::vector<double>{pTBins, pTBins + npTBins + 1};
+
+// default values for the cuts
+//                                            m   CPA  d0Jpsi  d0Pi pTJpsi pTPi chi2PCA
+constexpr double cuts[npTBins][nCutVars] = {{0.5, 0.80, 0.001, 0.001, 3.0, 0.15, 1.},  /* 0<pt<0.5 */
+                                            {0.5, 0.80, 0.001, 0.001, 3.0, 0.15, 1.},  /* 0.5<pt<1 */
+                                            {0.5, 0.80, 0.001, 0.001, 3.0, 0.15, 1.},  /* 1<pt<2   */
+                                            {0.5, 0.80, 0.001, 0.001, 3.0, 0.15, 1.},  /* 2<pt<3   */
+                                            {0.5, 0.80, 0.001, 0.001, 3.0, 0.15, 1.},  /* 3<pt<4   */
+                                            {0.5, 0.80, 0.001, 0.001, 3.0, 0.15, 1.},  /* 4<pt<5   */
+                                            {0.5, 0.80, 0.001, 0.001, 3.0, 0.15, 1.},  /* 5<pt<7   */
+                                            {0.5, 0.80, 0.001, 0.001, 3.0, 0.15, 1.},  /* 7<pt<10  */
+                                            {0.5, 0.80, 0.001, 0.001, 3.0, 0.15, 1.}}; /* 10<pt<15 */
+// row labels
+static const std::vector<std::string> pTBinLabels = {
+  "pT bin 0",
+  "pT bin 1",
+  "pT bin 2",
+  "pT bin 3",
+  "pT bin 4",
+  "pT bin 5",
+  "pT bin 6",
+  "pT bin 7",
+  "pT bin 8"};
+// column labels
+static const std::vector<std::string> cutVarLabels = {"m", "CPA", "d0 Jpsi", "d0 Pi", "pT Jpsi", "pT Pi", "chi2PCA"};
+} // namespace hf_cuts_x_tojpsipipi
 } // namespace o2::analysis
 
 #endif // HF_SELECTOR_CUTS_H_
@@ -252,9 +252,159 @@ class TrackSelectorPID
     }
   }
 
-  /// Returns status of combined PID selection for a given track.
+  // RICH
+
+  /// Set pT range where RICH PID is applicable.
+  void setRangePtRICH(float ptMin, float ptMax)
+  {
+    mPtRICHMin = ptMin;
+    mPtRICHMax = ptMax;
+  }
+
+  /// Set RICH nσ range in which a track should be accepted.
+  void setRangeNSigmaRICH(float nsMin, float nsMax)
+  {
+    mNSigmaRICHMin = nsMin;
+    mNSigmaRICHMax = nsMax;
+  }
+
+  /// Set RICH nσ range in which a track should be conditionally accepted if combined with TOF.
+  void setRangeNSigmaRICHCondTOF(float nsMin, float nsMax)
+  {
+    mNSigmaRICHMinCondTOF = nsMin;
+    mNSigmaRICHMaxCondTOF = nsMax;
+  }
+
+  /// Checks if track is OK for RICH PID.
+  /// \param track  track
+  /// \return true if track is OK for RICH PID
+  template <typename T>
+  bool isValidTrackPIDRICH(const T& track)
+  {
+    if (track.richId() < 0) {
+      return false;
+    }
+    auto pt = track.pt();
+    return mPtRICHMin <= pt && pt <= mPtRICHMax;
+  }
+
+  /// Checks if track is compatible with given particle species hypothesis within given RICH nσ range.
+  /// \param track  track
+  /// \param conditionalTOF  variable to store the result of selection with looser cuts for conditional accepting of track if combined with TOF
+  /// \return true if track satisfies RICH PID hypothesis for given RICH nσ range
+  template <typename T>
+  bool isSelectedTrackPIDRICH(const T& track, bool& conditionalTOF)
+  {
+    // Accept if selection is disabled via large values.
+    if (mNSigmaRICHMin < -999. && mNSigmaRICHMax > 999.) {
+      return true;
+    }
+
+    // Get nσ for a given particle hypothesis.
+    double nSigma = 100.;
+    switch (mPdg) {
+      case kElectron: {
+        nSigma = track.rich().richNsigmaEl();
+        break;
+      }
+      case kMuonMinus: {
+        nSigma = track.rich().richNsigmaMu();
+        break;
+      }
+      case kPiPlus: {
+        nSigma = track.rich().richNsigmaPi();
+        break;
+      }
+      case kKPlus: {
+        nSigma = track.rich().richNsigmaKa();
+        break;
+      }
+      case kProton: {
+        nSigma = track.rich().richNsigmaPr();
+        break;
+      }
+      default: {
+        LOGF(error, "ERROR: RICH PID not implemented for PDG %d", mPdg);
+        assert(false);
+      }
+    }
+
+    if (mNSigmaRICHMinCondTOF < -999. && mNSigmaRICHMaxCondTOF > 999.) {
+      conditionalTOF = true;
+    } else {
+      conditionalTOF = mNSigmaRICHMinCondTOF <= nSigma && nSigma <= mNSigmaRICHMaxCondTOF;
+    }
+    return mNSigmaRICHMin <= nSigma && nSigma <= mNSigmaRICHMax;
+  }
+
+  /// Returns status of RICH PID selection for a given track.
+  /// \param track  track
+  /// \return RICH selection status (see TrackSelectorPID::Status)
+  template <typename T>
+  int getStatusTrackPIDRICH(const T& track)
+  {
+    if (isValidTrackPIDRICH(track)) {
+      bool condTOF = false;
+      if (isSelectedTrackPIDRICH(track, condTOF)) {
+        return Status::PIDAccepted; // accepted
+      } else if (condTOF) {
+        return Status::PIDConditional; // potential to be accepted if combined with TOF
+      } else {
+        return Status::PIDRejected; // rejected
+      }
+    } else {
+      return Status::PIDNotApplicable; // PID not applicable
+    }
+  }
+
+  // MID
+
+  /// Checks if track is OK for MID PID.
+  /// \param track  track
+  /// \return true if track is OK for MID PID
+  template <typename T>
+  bool isValidTrackPIDMID(const T& track)
+  {
+    return track.midId() > -1;
+  }
+
+  /// Checks if track is compatible with muon hypothesis in the MID detector.
+  /// \param track  track
+  /// \return true if track has been identified as muon by the MID detector
+  template <typename T>
+  bool isSelectedTrackPIDMID(const T& track)
+  {
+    if (mPdg != kMuonMinus) {
+      return false;
+    }
+    return track.mid().midIsMuon() == 1; // FIXME: change to return track.midIsMuon() once the column is bool.
+  }
+
+  /// Returns status of MID PID selection for a given track.
+  /// \param track  track
+  /// \return MID selection status (see TrackSelectorPID::Status)
+  template <typename T>
+  int getStatusTrackPIDMID(const T& track)
+  {
+    if (mPdg != kMuonMinus) {
+      return Status::PIDRejected;
+    }
+    if (isValidTrackPIDMID(track)) {
+      if (isSelectedTrackPIDMID(track)) {
+        return Status::PIDAccepted; // accepted
+      } else {
+        return Status::PIDRejected; // rejected
+      }
+    } else {
+      return Status::PIDNotApplicable; // PID not applicable
+    }
+  }
+
+  // Combined selection (TPC + TOF)
+
+  /// Returns status of combined PID (TPC + TOF) selection for a given track.
   /// \param track  track
-  /// \return combined-selection status (see TrackSelectorPID::Status)
+  /// \return status of combined PID (TPC + TOF) (see TrackSelectorPID::Status)
   template <typename T>
   int getStatusTrackPIDAll(const T& track)
   {
@@ -273,6 +423,60 @@ class TrackSelectorPID
     return Status::PIDNotApplicable; // (NotApplicable for one detector) and (NotApplicable or Conditional for the other)
   }
 
+  /// Checks whether a track is identified as electron and rejected as pion by TOF or RICH.
+  /// \param track  track
+  /// \param useTOF  switch to use TOF
+  /// \param useRICH  switch to use RICH
+  /// \return true if track is selected by TOF or RICH
+  /// \note Ported from https://github.com/feisenhu/ALICE3-LoI-LMee/blob/main/efficiency/macros/anaEEstudy.cxx
+  template <typename T>
+  bool isElectronAndNotPion(const T& track, bool useTOF = true, bool useRICH = true)
+  {
+    bool isSelTOF = false;
+    bool isSelRICH = false;
+    bool hasRICH = track.richId() > -1;
+    bool hasTOF = isValidTrackPIDTOF(track);
+    auto nSigmaTOFEl = track.tofNSigmaEl();
+    auto nSigmaTOFPi = track.tofNSigmaPi();
+    auto nSigmaRICHEl = hasRICH ? track.rich().richNsigmaEl() : -1000.;
+    auto nSigmaRICHPi = hasRICH ? track.rich().richNsigmaPi() : -1000.;
+    auto p = track.p();
+
+    // TOF
+    if (useTOF && hasTOF && (p < 0.6)) {
+      if (p > 0.4 && hasRICH) {
+        if ((std::abs(nSigmaTOFEl) < mNSigmaTOFMax) && (std::abs(nSigmaRICHEl) < mNSigmaRICHMax)) {
+          isSelTOF = true; // is selected as electron by TOF and RICH
+        }
+      } else if (p <= 0.4) {
+        if (std::abs(nSigmaTOFEl) < mNSigmaTOFMax) {
+          isSelTOF = true; // is selected as electron by TOF
+        }
+      } else {
+        isSelTOF = false; // This is rejecting all the heavier particles which do not have a RICH signal in the p area of 0.4-0.6 GeV/c
+      }
+      if (std::abs(nSigmaTOFPi) < mNSigmaTOFMax) {
+        isSelTOF = false; // is selected as pion by TOF
+      }
+    } else {
+      isSelTOF = false;
+    }
+
+    // RICH
+    if (useRICH && hasRICH) {
+      if (std::abs(nSigmaRICHEl) < mNSigmaRICHMax) {
+        isSelRICH = true; // is selected as electron by RICH
+      }
+      if ((std::abs(nSigmaRICHPi) < mNSigmaRICHMax) && (p > 1.0) && (p < 2.0)) {
+        isSelRICH = false; // is selected as pion by RICH
+      }
+    } else {
+      isSelRICH = false;
+    }
+
+    return isSelRICH || isSelTOF;
+  }
+
  private:
   uint mPdg = kPiPlus; ///< PDG code of the expected particle
 
@@ -291,6 +495,14 @@ class TrackSelectorPID
   float mNSigmaTOFMax = 3.;            ///< maximum number of TOF σ
   float mNSigmaTOFMinCondTPC = -1000.; ///< minimum number of TOF σ if combined with TPC
   float mNSigmaTOFMaxCondTPC = 1000.;  ///< maximum number of TOF σ if combined with TPC
+
+  // RICH
+  float mPtRICHMin = 0.;                ///< minimum pT for RICH PID [GeV/c]
+  float mPtRICHMax = 100.;              ///< maximum pT for RICH PID [GeV/c]
+  float mNSigmaRICHMin = -3.;           ///< minimum number of RICH σ
+  float mNSigmaRICHMax = 3.;            ///< maximum number of RICH σ
+  float mNSigmaRICHMinCondTOF = -1000.; ///< minimum number of RICH σ if combined with TOF
+  float mNSigmaRICHMaxCondTOF = 1000.;  ///< maximum number of RICH σ if combined with TOF
 };
 
 #endif // O2_ANALYSIS_TRACKSELECTORPID_H_