|
| 1 | +// Copyright 2019-2022 CERN and copyright holders of ALICE O2. |
| 2 | +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. |
| 3 | +// All rights not expressly granted are reserved. |
| 4 | +// |
| 5 | +// This software is distributed under the terms of the GNU General Public |
| 6 | +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". |
| 7 | +// |
| 8 | +// In applying this license CERN does not waive the privileges and immunities |
| 9 | +// granted to it by virtue of its status as an Intergovernmental Organization |
| 10 | +// or submit itself to any jurisdiction. |
| 11 | + |
| 12 | +/// \file closePairRejection.h |
| 13 | +/// \brief Definition of ClosePairRejection class |
| 14 | +/// \author Anton Riedel, TU München, anton.riedel@tum.de |
| 15 | + |
| 16 | +#ifndef PWGCF_FEMTO_CORE_CLOSEPAIRREJECTION_H_ |
| 17 | +#define PWGCF_FEMTO_CORE_CLOSEPAIRREJECTION_H_ |
| 18 | + |
| 19 | +#include "PWGCF/Femto/Core/dataTypes.h" |
| 20 | +#include "PWGCF/Femto/Core/femtoUtils.h" |
| 21 | +#include "PWGCF/Femto/Core/histManager.h" |
| 22 | +#include "PWGCF/Femto/Core/modes.h" |
| 23 | +#include "PWGCF/Femto/DataModel/FemtoTables.h" |
| 24 | + |
| 25 | +#include "Framework/HistogramRegistry.h" |
| 26 | + |
| 27 | +#include <array> |
| 28 | +#include <map> |
| 29 | +#include <memory> |
| 30 | +#include <numeric> |
| 31 | +#include <string> |
| 32 | +#include <type_traits> |
| 33 | +#include <unordered_map> |
| 34 | +#include <vector> |
| 35 | + |
| 36 | +namespace o2::analysis::femto |
| 37 | +{ |
| 38 | +namespace closepairrejection |
| 39 | +{ |
| 40 | +// enum for track histograms |
| 41 | +enum CprHist { |
| 42 | + // kinemtics |
| 43 | + kAverage, |
| 44 | + kRadius0, |
| 45 | + kRadius1, |
| 46 | + kRadius2, |
| 47 | + kRadius3, |
| 48 | + kRadius4, |
| 49 | + kRadius5, |
| 50 | + kRadius6, |
| 51 | + kRadius7, |
| 52 | + kRadius8, |
| 53 | + kCprHistogramLast |
| 54 | +}; |
| 55 | + |
| 56 | +struct ConfCpr : o2::framework::ConfigurableGroup { |
| 57 | + std::string prefix = std::string("ClosePairRejection"); |
| 58 | + o2::framework::Configurable<bool> on{"on", true, "Turn on CPR"}; |
| 59 | + o2::framework::Configurable<float> detaMax{"detaMax", 0.01f, "Maximium deta"}; |
| 60 | + o2::framework::Configurable<float> dphistarMax{"dphistarMax", 0.01f, "Maximum dphistar"}; |
| 61 | + o2::framework::ConfigurableAxis binningDeta{"binningDeta", {{200, -0.2, 0.2}}, "deta"}; |
| 62 | + o2::framework::ConfigurableAxis binningDphistar{"binningDphistar", {{200, -0.2, 0.2}}, "dphi"}; |
| 63 | +}; |
| 64 | + |
| 65 | +// tpc radii for computing phistar |
| 66 | +constexpr int kNradii = 9; |
| 67 | +constexpr std::array<float, kNradii> kTpcRadius = {85., 105., 125., 145., 165., 185., 205., 225., 245.}; |
| 68 | + |
| 69 | +// directory names |
| 70 | +constexpr char PrefixTrackTrackSe[] = "CPR_TrackTrack/SE/"; |
| 71 | +constexpr char PrefixTrackTrackMe[] = "CPR_TrackTrack/ME/"; |
| 72 | +constexpr char PrefixTrackPosDauSe[] = "CPR_TrackPosDau/SE/"; |
| 73 | +constexpr char PrefixTrackNegDauSe[] = "CPR_TrackNegDau/SE/"; |
| 74 | +constexpr char PrefixTrackPosDauMe[] = "CPR_TrackPosDau/ME/"; |
| 75 | +constexpr char PrefixTrackNegDauMe[] = "CPR_TrackNegDau/ME/"; |
| 76 | + |
| 77 | +// must be in sync with enum TrackVariables |
| 78 | +// the enum gives the correct index in the array |
| 79 | +constexpr std::array<histmanager::HistInfo<CprHist>, kCprHistogramLast> HistTable = { |
| 80 | + {{kAverage, o2::framework::kTH2F, "hAverage", "#Delta #eta vs #Delta #phi* (averaged over all radii); #Delta #eta; #Delta #phi*"}, |
| 81 | + {kRadius0, o2::framework::kTH2F, "hRadius0", "Radius 0: #Delta #eta vs #Delta #phi*; #Delta #eta; #Delta #phi*"}, |
| 82 | + {kRadius1, o2::framework::kTH2F, "hRadius1", "Radius 1: #Delta #eta vs #Delta #phi*; #Delta #eta; #Delta #phi*"}, |
| 83 | + {kRadius2, o2::framework::kTH2F, "hRadius2", "Radius 2: #Delta #eta vs #Delta #phi*; #Delta #eta; #Delta #phi*"}, |
| 84 | + {kRadius3, o2::framework::kTH2F, "hRadius3", "Radius 3: #Delta #eta vs #Delta #phi*; #Delta #eta; #Delta #phi*"}, |
| 85 | + {kRadius4, o2::framework::kTH2F, "hRadius4", "Radius 4: #Delta #eta vs #Delta #phi*; #Delta #eta; #Delta #phi*"}, |
| 86 | + {kRadius5, o2::framework::kTH2F, "hRadius5", "Radius 5: #Delta #eta vs #Delta #phi*; #Delta #eta; #Delta #phi*"}, |
| 87 | + {kRadius6, o2::framework::kTH2F, "hRadius6", "Radius 6: #Delta #eta vs #Delta #phi*; #Delta #eta; #Delta #phi*"}, |
| 88 | + {kRadius7, o2::framework::kTH2F, "hRadius7", "Radius 7: #Delta #eta vs #Delta #phi*; #Delta #eta; #Delta #phi*"}, |
| 89 | + {kRadius8, o2::framework::kTH2F, "hRadius8", "Radius 8: #Delta #eta vs #Delta #phi*; #Delta #eta; #Delta #phi*"}}}; |
| 90 | + |
| 91 | +template <typename T> |
| 92 | +auto makeCprHistSpecMap(const T& confCpr) |
| 93 | +{ |
| 94 | + std::map<CprHist, std::vector<framework::AxisSpec>> specs; |
| 95 | + for (int i = 0; i < kCprHistogramLast; ++i) { |
| 96 | + specs[static_cast<CprHist>(i)] = {confCpr.binningDeta, confCpr.binningDphistar}; |
| 97 | + } |
| 98 | + return specs; |
| 99 | +}; |
| 100 | + |
| 101 | +template <const char* prefix> |
| 102 | +class CloseTrackRejection |
| 103 | +{ |
| 104 | + public: |
| 105 | + CloseTrackRejection() {} |
| 106 | + |
| 107 | + void init(o2::framework::HistogramRegistry* registry, std::map<CprHist, std::vector<o2::framework::AxisSpec>>& specs, float detaMax, float dphistarMax) |
| 108 | + { |
| 109 | + mHistogramRegistry = registry; |
| 110 | + mDetaMax = detaMax; |
| 111 | + mDphistarMax = dphistarMax; |
| 112 | + |
| 113 | + for (int i = 0; i < kCprHistogramLast; ++i) { |
| 114 | + auto hist = static_cast<CprHist>(i); |
| 115 | + mHistogramRegistry->add(std::string(prefix) + GetHistNamev2(hist, HistTable), GetHistDesc(hist, HistTable), GetHistType(hist, HistTable), {specs.at(hist)}); |
| 116 | + } |
| 117 | + } |
| 118 | + |
| 119 | + void setMagField(float magField) { mMagField = magField; } |
| 120 | + |
| 121 | + void reset() |
| 122 | + { |
| 123 | + mSameCharge = false; |
| 124 | + mAverageDphistar = 0.f; |
| 125 | + mDeta = 0.f; |
| 126 | + mDphistar.fill(0.f); |
| 127 | + } |
| 128 | + |
| 129 | + template <typename T1, typename T2> |
| 130 | + void compute(const T1& track1, const T2& track2) |
| 131 | + { |
| 132 | + reset(); |
| 133 | + if (track1.sign() != track2.sign()) { |
| 134 | + return; |
| 135 | + } |
| 136 | + mSameCharge = true; |
| 137 | + mDeta = track1.eta() - track2.eta(); |
| 138 | + for (size_t i = 0; i < kTpcRadius.size(); ++i) { |
| 139 | + auto phistar1 = utils::dphistar(mMagField, kTpcRadius[i], track1.sign(), track1.pt(), track1.phi()); |
| 140 | + auto phistar2 = utils::dphistar(mMagField, kTpcRadius[i], track2.sign(), track2.pt(), track2.phi()); |
| 141 | + |
| 142 | + if (phistar1 && phistar2) { |
| 143 | + // if the calculation for one phistar fails, keep the default value, which is 0 |
| 144 | + // this makes it more likelier for the pair to be rejected sind the averave will be biased towards lower values |
| 145 | + mDphistar[i] = phistar1.value() - phistar2.value(); |
| 146 | + } |
| 147 | + } |
| 148 | + mAverageDphistar = std::accumulate(mDphistar.begin(), mDphistar.end(), 0.f) / mDphistar.size(); |
| 149 | + } |
| 150 | + |
| 151 | + void fill() |
| 152 | + { |
| 153 | + // fill average hist |
| 154 | + mHistogramRegistry->fill(HIST(prefix) + HIST(GetHistName(kAverage, HistTable)), mDeta, mAverageDphistar); |
| 155 | + |
| 156 | + // fill radii hists |
| 157 | + for (int i = 0; i < kNradii; ++i) { |
| 158 | + mHistogramRegistry->fill(HIST(prefix) + HIST(GetHistName(static_cast<CprHist>(kRadius0 + 1), HistTable)), mDeta, mDphistar.at(i)); |
| 159 | + } |
| 160 | + } |
| 161 | + |
| 162 | + bool isClosePair() const |
| 163 | + { |
| 164 | + if (!mSameCharge) { |
| 165 | + return false; |
| 166 | + } else { |
| 167 | + return ((mAverageDphistar * mAverageDphistar) / (mDphistarMax * mDphistarMax) + (mDeta * mDeta) / (mDetaMax * mDetaMax)) < 1.f; |
| 168 | + } |
| 169 | + } |
| 170 | + |
| 171 | + private: |
| 172 | + float mMagField = 0.f; |
| 173 | + float mAverageDphistar = 0.f; |
| 174 | + bool mSameCharge = false; |
| 175 | + float mDeta = 0.f; |
| 176 | + float mDetaMax = 0.f; |
| 177 | + float mDphistarMax = 0.f; |
| 178 | + std::array<float, kNradii> mDphistar = {}; |
| 179 | + |
| 180 | + o2::framework::HistogramRegistry* mHistogramRegistry = nullptr; |
| 181 | +}; |
| 182 | + |
| 183 | +template <const char* prefix> |
| 184 | +class ClosePairRejectionTrackTrack |
| 185 | +{ |
| 186 | + public: |
| 187 | + void init(o2::framework::HistogramRegistry* registry, std::map<CprHist, std::vector<o2::framework::AxisSpec>>& specs, float detaMax, float dphistarMax, bool isActivated) |
| 188 | + { |
| 189 | + mIsActivated = isActivated; |
| 190 | + mCtr.init(registry, specs, detaMax, dphistarMax); |
| 191 | + } |
| 192 | + |
| 193 | + void setMagField(float magField) { mCtr.setMagField(magField); } |
| 194 | + template <typename T1, typename T2> |
| 195 | + void setPair(const T1& track1, const T2& track2) |
| 196 | + { |
| 197 | + mCtr.compute(track1, track2); |
| 198 | + } |
| 199 | + bool isClosePair() const { return mCtr.isClosePair(); } |
| 200 | + void fill() { mCtr.fill(); } |
| 201 | + bool isActivated() const { return mIsActivated; } |
| 202 | + |
| 203 | + private: |
| 204 | + CloseTrackRejection<prefix> mCtr; |
| 205 | + bool mIsActivated = true; |
| 206 | +}; |
| 207 | + |
| 208 | +template <const char* prefixPosDau, const char* prefixNegDau> |
| 209 | +class ClosePairRejectionTrackV0 // can also be used for any particle type that has pos/neg daughters, like resonances |
| 210 | +{ |
| 211 | + public: |
| 212 | + void init(o2::framework::HistogramRegistry* registry, std::map<CprHist, std::vector<o2::framework::AxisSpec>>& specs, float detaMax, float dphistarMax, bool isActivated) |
| 213 | + { |
| 214 | + mIsActivated = isActivated; |
| 215 | + mCtrPosDau.init(registry, specs, detaMax, dphistarMax); |
| 216 | + mCtrNegDau.init(registry, specs, detaMax, dphistarMax); |
| 217 | + } |
| 218 | + |
| 219 | + void setMagField(float magField) |
| 220 | + { |
| 221 | + mCtrPosDau.setMagField(magField); |
| 222 | + mCtrNegDau.setMagField(magField); |
| 223 | + } |
| 224 | + template <typename T1, typename T2, typename T3> |
| 225 | + void setPair(const T1& track, const T2& v0, const T3 /*trackTable*/) |
| 226 | + { |
| 227 | + auto posDaughter = v0.template posDau_as<T3>(); |
| 228 | + auto negDaughter = v0.template negDau_as<T3>(); |
| 229 | + mCtrPosDau.compute(track, posDaughter); |
| 230 | + mCtrNegDau.compute(track, negDaughter); |
| 231 | + } |
| 232 | + bool isClosePair() const { return mCtrPosDau.isClosePair() || mCtrNegDau.isClosePair(); } |
| 233 | + void fill() |
| 234 | + { |
| 235 | + mCtrPosDau.fill(); |
| 236 | + mCtrNegDau.fill(); |
| 237 | + } |
| 238 | + bool isActivated() const { return mIsActivated; } |
| 239 | + |
| 240 | + private: |
| 241 | + CloseTrackRejection<prefixPosDau> mCtrPosDau; |
| 242 | + CloseTrackRejection<prefixNegDau> mCtrNegDau; |
| 243 | + bool mIsActivated = true; |
| 244 | +}; |
| 245 | + |
| 246 | +}; // namespace closepairrejection |
| 247 | +}; // namespace o2::analysis::femto |
| 248 | +#endif // PWGCF_FEMTO_CORE_CLOSEPAIRREJECTION_H_ |
0 commit comments