From 9ef364069e0bf1468ab8bd0496182e95e7eadd3e Mon Sep 17 00:00:00 2001 From: Anton Riedel Date: Wed, 14 Jan 2026 16:11:13 +0100 Subject: [PATCH] Feat: add 3 body processing --- PWGCF/Femto/Core/closePairRejection.h | 19 +- PWGCF/Femto/Core/closeTripletRejection.h | 96 +++ PWGCF/Femto/Core/mcBuilder.h | 7 +- PWGCF/Femto/Core/pairBuilder.h | 52 +- PWGCF/Femto/Core/pairHistManager.h | 73 ++- PWGCF/Femto/Core/pairProcessHelpers.h | 172 ++--- PWGCF/Femto/Core/trackBuilder.h | 18 +- PWGCF/Femto/Core/trackHistManager.h | 2 + PWGCF/Femto/Core/tripletBuilder.h | 366 +++++++++++ PWGCF/Femto/Core/tripletCleaner.h | 64 ++ PWGCF/Femto/Core/tripletHistManager.h | 588 ++++++++++++++++++ PWGCF/Femto/Core/tripletProcessHelpers.h | 541 ++++++++++++++++ PWGCF/Femto/TableProducer/femtoProducer.cxx | 70 +-- PWGCF/Femto/Tasks/CMakeLists.txt | 5 + .../Tasks/femtoTripletTrackTrackTrack.cxx | 188 ++++++ 15 files changed, 1991 insertions(+), 270 deletions(-) create mode 100644 PWGCF/Femto/Core/closeTripletRejection.h create mode 100644 PWGCF/Femto/Core/tripletBuilder.h create mode 100644 PWGCF/Femto/Core/tripletCleaner.h create mode 100644 PWGCF/Femto/Core/tripletHistManager.h create mode 100644 PWGCF/Femto/Core/tripletProcessHelpers.h create mode 100644 PWGCF/Femto/Tasks/femtoTripletTrackTrackTrack.cxx diff --git a/PWGCF/Femto/Core/closePairRejection.h b/PWGCF/Femto/Core/closePairRejection.h index 9ac2ea0ad8b..1739fdefa5e 100644 --- a/PWGCF/Femto/Core/closePairRejection.h +++ b/PWGCF/Femto/Core/closePairRejection.h @@ -65,8 +65,8 @@ struct ConfCpr : o2::framework::ConfigurableGroup { o2::framework::Configurable dphistarMax{"dphistarMax", 0.01f, "Maximum dphistar"}; o2::framework::Configurable detaCenter{"detaCenter", 0.f, "Center of deta cut"}; o2::framework::Configurable dphistarCenter{"dphistarCenter", 0.f, "Center of dphistar cut"}; - o2::framework::Configurable kstarMin{"kstarMin", -1.f, "Minimum kstar of pair for plotting (Set to negative value to turn off the cut)"}; - o2::framework::Configurable kstarMax{"kstarMax", -1.f, "Maximum kstar of pair for plotting (Set to negative value to turn off the cut)"}; + o2::framework::Configurable kinematicMin{"kinematicMin", -1.f, "Minimum kstar/Q3 of pair/triplet for plotting (Set to negative value to turn off the cut)"}; + o2::framework::Configurable kinematicMax{"kinematicMax", -1.f, "Maximum kstar/Q3 of pair/triplet for plotting (Set to negative value to turn off the cut)"}; o2::framework::ConfigurableAxis binningDeta{"binningDeta", {{250, -0.5, 0.5}}, "deta"}; o2::framework::ConfigurableAxis binningDphistar{"binningDphistar", {{250, -0.5, 0.5}}, "dphi"}; }; @@ -79,6 +79,7 @@ constexpr const char PrefixCprV0DaughterV0DaughterPos[] = "CprV0DaughterV0Daught constexpr const char PrefixCprV0DaughterV0DaughterNeg[] = "CprV0DaughterV0DaughterNeg"; constexpr const char PrefixCprTrackCascadeBachelor[] = "CprTrackCascadeBachelor"; +// pairs using ConfCprTrackTrack = ConfCpr; using ConfCprTrackV0Daughter = ConfCpr; using ConfCprTrackResonanceDaughter = ConfCpr; @@ -168,8 +169,8 @@ class CloseTrackRejection mCutAverage = confCpr.cutAverage.value; mCutAnyRadius = confCpr.cutAnyRadius.value; - mKstarMin = confCpr.kstarMin.value; - mKstarMax = confCpr.kstarMax.value; + mKinematicMin = confCpr.kinematicMin.value; + mKinematicMax = confCpr.kinematicMax.value; mPlotAverage = confCpr.plotAverage.value; mPlotAllRadii = confCpr.plotAllRadii.value; @@ -229,17 +230,17 @@ class CloseTrackRejection } } - void fill(float kstar) + void fill(float kinematic) { if (!mIsActivated) { return; } - if (mKstarMin > 0.f && kstar < mKstarMin) { + if (mKinematicMin > 0.f && kinematic < mKinematicMin) { return; } - if (mKstarMax > 0.f && kstar > mKstarMax) { + if (mKinematicMax > 0.f && kinematic > mKinematicMax) { return; } @@ -312,8 +313,8 @@ class CloseTrackRejection bool mPlotAllRadii = false; bool mPlotAverage = false; - float mKstarMin = -1.f; - float mKstarMax = -1.f; + float mKinematicMin = -1.f; + float mKinematicMax = -1.f; bool mCutAverage = false; bool mCutAnyRadius = false; diff --git a/PWGCF/Femto/Core/closeTripletRejection.h b/PWGCF/Femto/Core/closeTripletRejection.h new file mode 100644 index 00000000000..e929f0c4078 --- /dev/null +++ b/PWGCF/Femto/Core/closeTripletRejection.h @@ -0,0 +1,96 @@ +// Copyright 2019-2022 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 closeTripletRejection.h +/// \brief Definition of CloseTripletRejection class +/// \author Anton Riedel, TU München, anton.riedel@tum.de + +#ifndef PWGCF_FEMTO_CORE_CLOSETRIPLETREJECTION_H_ +#define PWGCF_FEMTO_CORE_CLOSETRIPLETREJECTION_H_ + +#include "PWGCF/Femto/Core/closePairRejection.h" + +#include "Framework/HistogramRegistry.h" +#include "Framework/HistogramSpec.h" + +#include +#include +#include + +namespace o2::analysis::femto +{ +namespace closetripletrejection +{ + +constexpr const char PrefixCtrTrackTrackTrack[] = "CtrTrackTrackTrack"; +using ConfCtrTrackTrackTrack = closepairrejection::ConfCpr; + +// directory names +constexpr char PrefixTrack1Track2Se[] = "CPR_Track1Track2/SE/"; +constexpr char PrefixTrack2Track3Se[] = "CPR_Track2Track3/SE/"; +constexpr char PrefixTrack1Track3Se[] = "CPR_Track1Track3/SE/"; +constexpr char PrefixTrack1Track2Me[] = "CPR_Track1Track2/ME/"; +constexpr char PrefixTrack2Track3Me[] = "CPR_Track2Track3/ME/"; +constexpr char PrefixTrack1Track3Me[] = "CPR_Track1Track3/ME/"; + +template +class CloseTripletRejectionTrackTrackTrack +{ + public: + template + void init(o2::framework::HistogramRegistry* registry, + std::map> const& specs, + T const& confCpr, + int absChargeTrack1, + int absChargeTrack2, + int absChargeTrack3) + { + mCtr1.init(registry, specs, confCpr, absChargeTrack1, absChargeTrack2); + mCtr2.init(registry, specs, confCpr, absChargeTrack2, absChargeTrack3); + mCtr3.init(registry, specs, confCpr, absChargeTrack1, absChargeTrack3); + } + + void setMagField(float magField) + { + mCtr1.setMagField(magField); + mCtr2.setMagField(magField); + mCtr3.setMagField(magField); + } + template + void setTriplet(T1 const& track1, T2 const& track2, T3 const& track3, T4 const& /*tracks*/) + { + mCtr1.compute(track1, track2); + mCtr2.compute(track2, track3); + mCtr3.compute(track1, track3); + } + bool isCloseTriplet() const + { + return mCtr1.isClosePair() || mCtr2.isClosePair() || mCtr3.isClosePair(); + } + + void fill(float q3) + { + mCtr1.fill(q3); + mCtr2.fill(q3); + mCtr3.fill(q3); + } + + private: + closepairrejection::CloseTrackRejection mCtr1; + closepairrejection::CloseTrackRejection mCtr2; + closepairrejection::CloseTrackRejection mCtr3; +}; + +}; // namespace closetripletrejection +}; // namespace o2::analysis::femto +#endif // PWGCF_FEMTO_CORE_CLOSETRIPLETREJECTION_H_ diff --git a/PWGCF/Femto/Core/mcBuilder.h b/PWGCF/Femto/Core/mcBuilder.h index bf4f023f78d..67385354077 100644 --- a/PWGCF/Femto/Core/mcBuilder.h +++ b/PWGCF/Femto/Core/mcBuilder.h @@ -180,12 +180,17 @@ class McBuilder bool fillAnyTable() const { return mFillAnyTable; } - void reset() + template + void reset(T1 const& mcCollisions, T2 const& mcParticles) { mCollisionMap.clear(); + mCollisionMap.reserve(mcCollisions.size()); mMcParticleMap.clear(); + mMcParticleMap.reserve(mcParticles.size()); mMcMotherMap.clear(); + mMcMotherMap.reserve(mcParticles.size()); mMcPartonicMotherMap.clear(); + mMcPartonicMotherMap.reserve(mcParticles.size()); } private: diff --git a/PWGCF/Femto/Core/pairBuilder.h b/PWGCF/Femto/Core/pairBuilder.h index 41d8b8d894b..b9ac36cf4ff 100644 --- a/PWGCF/Femto/Core/pairBuilder.h +++ b/PWGCF/Femto/Core/pairBuilder.h @@ -45,12 +45,7 @@ namespace o2::analysis::femto namespace pairbuilder { -template +template class PairTrackTrackBuilder { public: @@ -128,6 +123,7 @@ class PairTrackTrackBuilder randomSeed = static_cast(confMixing.seed.value); } mRng = std::mt19937(randomSeed); + mDist = std::uniform_int_distribution(static_cast(pairprocesshelpers::kOrder12), static_cast(pairprocesshelpers::kOrder21)); } } @@ -142,7 +138,11 @@ class PairTrackTrackBuilder } mColHistManager.template fill(col); mCprSe.setMagField(col.magField()); - pairprocesshelpers::processSameEvent(trackSlice1, trackTable, col, mTrackHistManager1, mPairHistManagerSe, mCprSe, mPc, mRng, mMixIdenticalParticles); + pairprocesshelpers::PairOrder pairOrder = pairprocesshelpers::kOrder12; + if (mMixIdenticalParticles) { + pairOrder = static_cast(mDist(mRng)); + } + pairprocesshelpers::processSameEvent(trackSlice1, trackTable, col, mTrackHistManager1, mPairHistManagerSe, mCprSe, mPc, pairOrder); } else { auto trackSlice1 = partition1->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); auto trackSlice2 = partition2->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); @@ -166,7 +166,11 @@ class PairTrackTrackBuilder } mColHistManager.template fill(col, mcCols); mCprSe.setMagField(col.magField()); - pairprocesshelpers::processSameEvent(trackSlice1, trackTable, mcParticles, mcMothers, mcPartonicMothers, col, mcCols, mTrackHistManager1, mPairHistManagerSe, mCprSe, mPc, mRng, mMixIdenticalParticles); + pairprocesshelpers::PairOrder pairOrder = pairprocesshelpers::kOrder12; + if (mMixIdenticalParticles) { + pairOrder = static_cast(mDist(mRng)); + } + pairprocesshelpers::processSameEvent(trackSlice1, trackTable, mcParticles, mcMothers, mcPartonicMothers, col, mcCols, mTrackHistManager1, mPairHistManagerSe, mCprSe, mPc, pairOrder); } else { auto trackSlice1 = partition1->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); auto trackSlice2 = partition2->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); @@ -186,13 +190,13 @@ class PairTrackTrackBuilder if (mSameSpecies) { switch (mMixingPolicy) { case static_cast(pairhistmanager::kVtxMult): - pairprocesshelpers::processMixedEvent(cols, partition1, trackTable, cache, binsVtxMult, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + pairprocesshelpers::processMixedEvent(cols, partition1, partition1, trackTable, cache, binsVtxMult, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); break; case static_cast(pairhistmanager::kVtxCent): - pairprocesshelpers::processMixedEvent(cols, partition1, trackTable, cache, binsVtxCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + pairprocesshelpers::processMixedEvent(cols, partition1, partition1, trackTable, cache, binsVtxCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); break; case static_cast(pairhistmanager::kVtxMultCent): - pairprocesshelpers::processMixedEvent(cols, partition1, trackTable, cache, binsVtxMultCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + pairprocesshelpers::processMixedEvent(cols, partition1, partition1, trackTable, cache, binsVtxMultCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); break; default: LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; @@ -220,13 +224,13 @@ class PairTrackTrackBuilder if (mSameSpecies) { switch (mMixingPolicy) { case static_cast(pairhistmanager::kVtxMult): - pairprocesshelpers::processMixedEvent(cols, mcCols, partition1, trackTable, mcParticles, cache, binsVtxMult, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + pairprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition1, trackTable, mcParticles, cache, binsVtxMult, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); break; case static_cast(pairhistmanager::kVtxCent): - pairprocesshelpers::processMixedEvent(cols, mcCols, partition1, trackTable, mcParticles, cache, binsVtxCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + pairprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition1, trackTable, mcParticles, cache, binsVtxCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); break; case static_cast(pairhistmanager::kVtxMultCent): - pairprocesshelpers::processMixedEvent(cols, mcCols, partition1, trackTable, mcParticles, cache, binsVtxMultCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + pairprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition1, trackTable, mcParticles, cache, binsVtxMultCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); break; default: LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; @@ -257,11 +261,12 @@ class PairTrackTrackBuilder closepairrejection::ClosePairRejectionTrackTrack mCprSe; closepairrejection::ClosePairRejectionTrackTrack mCprMe; paircleaner::TrackTrackPairCleaner mPc; - std::mt19937 mRng; pairhistmanager::MixingPolicy mMixingPolicy = pairhistmanager::MixingPolicy::kVtxMult; - bool mSameSpecies = false; int mMixingDepth = 5; + bool mSameSpecies = false; bool mMixIdenticalParticles = false; + std::mt19937 mRng; + std::uniform_int_distribution<> mDist; }; template (col); mCprSe.setMagField(col.magField()); - pairprocesshelpers::processSameEvent(v0Slice1, trackTable, col, mV0HistManager1, mPairHistManagerSe, mCprSe, mPc, mRng, mMixIdenticalParticles); + pairprocesshelpers::PairOrder pairOrder = pairprocesshelpers::kOrder12; + if (mMixIdenticalParticles) { + pairOrder = static_cast(mDist(mRng)); + } + pairprocesshelpers::processSameEvent(v0Slice1, trackTable, col, mV0HistManager1, mPairHistManagerSe, mCprSe, mPc, pairOrder); } else { auto v0Slice1 = partition1->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); auto v0Slice2 = partition2->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); @@ -395,13 +404,13 @@ class PairV0V0Builder if (mSameSpecies) { switch (mMixingPolicy) { case static_cast(pairhistmanager::kVtxMult): - pairprocesshelpers::processMixedEvent(cols, partition1, trackTable, cache, binsVtxMult, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + pairprocesshelpers::processMixedEvent(cols, partition1, partition1, trackTable, cache, binsVtxMult, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); break; case static_cast(pairhistmanager::kVtxCent): - pairprocesshelpers::processMixedEvent(cols, partition1, trackTable, cache, binsVtxCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + pairprocesshelpers::processMixedEvent(cols, partition1, partition1, trackTable, cache, binsVtxCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); break; case static_cast(pairhistmanager::kVtxMultCent): - pairprocesshelpers::processMixedEvent(cols, partition1, trackTable, cache, binsVtxMultCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + pairprocesshelpers::processMixedEvent(cols, partition1, partition1, trackTable, cache, binsVtxMultCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); break; default: LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; @@ -432,11 +441,12 @@ class PairV0V0Builder closepairrejection::ClosePairRejectionV0V0 mCprSe; closepairrejection::ClosePairRejectionV0V0 mCprMe; paircleaner::V0V0PairCleaner mPc; - std::mt19937 mRng; pairhistmanager::MixingPolicy mMixingPolicy = pairhistmanager::MixingPolicy::kVtxMult; bool mSameSpecies = false; int mMixingDepth = 5; bool mMixIdenticalParticles = false; + std::mt19937 mRng; + std::uniform_int_distribution<> mDist; }; template , kPairHistogramLast> {kKstarVsMass1VsMass2VsMult, o2::framework::kTHnSparseF, "hKstarVsMass1VsMass2VsMult", "k* vs m_{1} vs m_{2} vs multiplicity; k* (GeV/#it{c}); m_{1} (GeV/#it{c}^{2}); m_{2} (GeV/#it{c}^{2}); Multiplicity"}, {kTrueKstarVsKstar, o2::framework::kTH2F, "hTrueKstarVsKstar", "k*_{True} vs k*; k*_{True} (GeV/#it{c}); k* (GeV/#it{c})"}, {kTrueKtVsKt, o2::framework::kTH2F, "hTrueKtVsKt", "k_{T,True} vs k_{T}; k_{T,True} (GeV/#it{c}); k_{T} (GeV/#it{c})"}, - {kTrueMtVsMt, o2::framework::kTH2F, "hTrueMtVsMt", "m_{T,True} vs m_{T}; m_{T,True} (GeV/#it{c}^{2}); m_{T,True} (GeV/#it{c}^{2})"}, + {kTrueMtVsMt, o2::framework::kTH2F, "hTrueMtVsMt", "m_{T,True} vs m_{T}; m_{T,True} (GeV/#it{c}^{2}); m_{T} (GeV/#it{c}^{2})"}, {kTrueMultVsMult, o2::framework::kTH2F, "hTrueMultVsMult", "Multiplicity_{True} vs Multiplicity; Multiplicity_{True} ; Multiplicity"}, {kTrueCentVsCent, o2::framework::kTH2F, "hTrueCentVsCent", "Centrality_{True} vs Centrality; Centrality_{True} (%); Centrality (%)"}, }}; @@ -319,25 +319,20 @@ class PairHistManager { // pt in track table is calculated from 1/signedPt from the original track table // in case of He with Z=2, we have to rescale the pt with the absolute charge - mParticle1 = ROOT::Math::PtEtaPhiMVector{mAbsCharge1 * particle1.pt(), particle1.eta(), particle1.phi(), mPdgMass1}; - mParticle2 = ROOT::Math::PtEtaPhiMVector{mAbsCharge2 * particle2.pt(), particle2.eta(), particle2.phi(), mPdgMass2}; - auto partSum = mParticle1 + mParticle2; + mParticle1 = ROOT::Math::PtEtaPhiMVector(mAbsCharge1 * particle1.pt(), particle1.eta(), particle1.phi(), mPdgMass1); + mParticle2 = ROOT::Math::PtEtaPhiMVector(mAbsCharge2 * particle2.pt(), particle2.eta(), particle2.phi(), mPdgMass2); // set kT - mKt = 0.5f * partSum.Pt(); + mKt = getKt(mParticle1, mParticle2); // set mT - mMt = computeMt(partSum); + mMt = getMt(mParticle1, mParticle2); - // Boost particle to the pair rest frame (Prf) and calculate k* (would be equivalent using particle 2) - // make a copy of particle 1 - auto particle1Prf = ROOT::Math::PtEtaPhiMVector(mParticle1); - // get lorentz boost into pair rest frame - ROOT::Math::Boost boostPrf(partSum.BoostToCM()); - // boost particle 1 into pair rest frame and calculate its momentum, which has the same value as k* - mKstar = boostPrf(particle1Prf).P(); + // set kstar + mKstar = getKstar(mParticle1, mParticle2); - // if one of the particles has a mass getter, we cache the value for the filling later + // if one of the particles has a mass getter (like lambda), we cache the value for the filling later + // for the compuation of kinematic variables like kstar we use the pdg values if constexpr (modes::hasMass(particleType1)) { mMass1 = particle1.mass(); } @@ -374,23 +369,16 @@ class PairHistManager auto mcParticle1 = particle1.template fMcParticle_as(); auto mcParticle2 = particle2.template fMcParticle_as(); - mParticle1 = ROOT::Math::PtEtaPhiMVector{mAbsCharge1 * mcParticle1.pt(), mcParticle1.eta(), mcParticle1.phi(), mPdgMass1}; - mParticle2 = ROOT::Math::PtEtaPhiMVector{mAbsCharge2 * mcParticle2.pt(), mcParticle2.eta(), mcParticle2.phi(), mPdgMass2}; - auto partSum = mParticle1 + mParticle2; + mTrueParticle1 = ROOT::Math::PtEtaPhiMVector(mAbsCharge1 * mcParticle1.pt(), mcParticle1.eta(), mcParticle1.phi(), mPdgMass1); + mTrueParticle2 = ROOT::Math::PtEtaPhiMVector(mAbsCharge2 * mcParticle2.pt(), mcParticle2.eta(), mcParticle2.phi(), mPdgMass2); // set kT - mTrueKt = partSum.Pt() / 2.f; + mTrueKt = getKt(mTrueParticle1, mTrueParticle2); // set mT - mTrueMt = computeMt(partSum); + mTrueMt = getMt(mTrueParticle1, mTrueParticle2); - // Boost particle to the pair rest frame (Prf) and calculate k* (would be equivalent using particle 2) - // make a copy of particle 1 - auto particle1Prf = ROOT::Math::PtEtaPhiMVector(mParticle1); - // get lorentz boost into pair rest frame - ROOT::Math::Boost boostPrf(partSum.BoostToCM()); - // boost particle 1 into pair rest frame and calculate its momentum, which has the same value as k* - mTrueKstar = boostPrf(particle1Prf).P(); + mTrueKstar = getKstar(mTrueParticle1, mTrueParticle2); } template @@ -606,26 +594,45 @@ class PairHistManager } } - float computeMt(ROOT::Math::PtEtaPhiMVector const& PairMomentum) + float getKt(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) + { + auto sum = part1 + part2; + return 0.5 * sum.Pt(); + } + + float getMt(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) { + auto sum = part1 + part2; float mt = 0; switch (mMtType) { case modes::TransverseMassType::kAveragePdgMass: - mt = std::hypot(0.5 * PairMomentum.Pt(), mAverageMass); + mt = std::hypot(0.5 * sum.Pt(), mAverageMass); break; case modes::TransverseMassType::kReducedPdgMass: - mt = std::hypot(0.5 * PairMomentum.Pt(), mReducedMass); + mt = std::hypot(0.5 * sum.Pt(), mReducedMass); break; case modes::TransverseMassType::kMt4Vector: - mt = PairMomentum.Mt() / 2.f; + mt = sum.Mt() / 2.f; break; default: - LOG(warn) << "Invalid transverse mass type, falling back to default..."; - mt = std::hypot(0.5 * PairMomentum.Pt(), mAverageMass); + LOG(fatal) << "Invalid transverse mass type, breaking..."; } return mt; } + float getKstar(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) + { + // compute pair momentum + auto sum = part1 + part2; + // Boost particle 1 to the pair rest frame (Prf) and calculate k* (would be equivalent using particle 2) + // make a copy of particle 1 + auto particle1Prf = ROOT::Math::PtEtaPhiMVector(mParticle1); + // get lorentz boost into pair rest frame + ROOT::Math::Boost boostPrf(sum.BoostToCM()); + // boost particle 1 into pair rest frame and calculate its momentum, which has the same value as k* + return boostPrf(particle1Prf).P(); + } + o2::framework::HistogramRegistry* mHistogramRegistry = nullptr; float mPdgMass1 = 0.f; float mPdgMass2 = 0.f; @@ -647,6 +654,8 @@ class PairHistManager float mCent = 0.f; // mc + ROOT::Math::PtEtaPhiMVector mTrueParticle1{}; + ROOT::Math::PtEtaPhiMVector mTrueParticle2{}; float mTrueKstar = 0.f; float mTrueKt = 0.f; float mTrueMt = 0.f; diff --git a/PWGCF/Femto/Core/pairProcessHelpers.h b/PWGCF/Femto/Core/pairProcessHelpers.h index 6380685bdc2..6190c29d10a 100644 --- a/PWGCF/Femto/Core/pairProcessHelpers.h +++ b/PWGCF/Femto/Core/pairProcessHelpers.h @@ -19,16 +19,18 @@ #include "PWGCF/Femto/Core/modes.h" #include "PWGCF/Femto/DataModel/FemtoTables.h" -#include "CommonConstants/MathConstants.h" #include "Framework/ASoAHelpers.h" -#include - namespace o2::analysis::femto { namespace pairprocesshelpers { +enum PairOrder : uint8_t { + kOrder12, + kOrder21 +}; + // process same event for identical particles template + typename T7> void processSameEvent(T1 const& SliceParticle, T2 const& TrackTable, T3 const& Collision, @@ -46,13 +47,11 @@ void processSameEvent(T1 const& SliceParticle, T5& PairHistManager, T6& CprManager, T7& PcManager, - T8& rng, - bool randomize) + PairOrder pairOrder) { for (auto const& part : SliceParticle) { ParticleHistManager.template fill(part, TrackTable); } - std::uniform_real_distribution dist(0.f, 1.f); for (auto const& [p1, p2] : o2::soa::combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(SliceParticle, SliceParticle))) { // check if pair is clean if (!PcManager.isCleanPair(p1, p2, TrackTable)) { @@ -64,12 +63,15 @@ void processSameEvent(T1 const& SliceParticle, continue; } // Randomize pair order if enabled - float threshold = 0.5f; - bool swapPair = randomize ? (dist(rng) > threshold) : false; - if (swapPair) { - PairHistManager.setPair(p2, p1, Collision); - } else { - PairHistManager.setPair(p1, p2, Collision); + switch (pairOrder) { + case kOrder12: + PairHistManager.setPair(p1, p2, Collision); + break; + case kOrder21: + PairHistManager.setPair(p2, p1, Collision); + break; + default: + PairHistManager.setPair(p1, p2, Collision); } // fill deta-dphi histograms with kstar cutoff CprManager.fill(PairHistManager.getKstar()); @@ -92,8 +94,7 @@ template + typename T11> void processSameEvent(T1 const& SliceParticle, T2 const& TrackTable, T3 const& mcParticles, @@ -105,13 +106,11 @@ void processSameEvent(T1 const& SliceParticle, T9& PairHistManager, T10& CprManager, T11& PcManager, - T12& rng, - bool randomize) + PairOrder pairOrder) { for (auto const& part : SliceParticle) { ParticleHistManager.template fill(part, TrackTable, mcParticles, mcMothers, mcPartonicMothers); } - std::uniform_real_distribution dist(0.f, 1.f); for (auto const& [p1, p2] : o2::soa::combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(SliceParticle, SliceParticle))) { // check if pair is clean if (!PcManager.isCleanPair(p1, p2, TrackTable, mcPartonicMothers)) { @@ -123,12 +122,15 @@ void processSameEvent(T1 const& SliceParticle, continue; } // Randomize pair order if enabled - float threshold = 0.5f; - bool swapPair = randomize ? (dist(rng) > threshold) : false; - if (swapPair) { - PairHistManager.setPairMc(p2, p1, mcParticles, Collision, mcCollisions); - } else { - PairHistManager.setPairMc(p1, p2, mcParticles, Collision, mcCollisions); + switch (pairOrder) { + case kOrder12: + PairHistManager.setPairMc(p1, p2, mcParticles, Collision, mcCollisions); + break; + case kOrder21: + PairHistManager.setPairMc(p2, p1, mcParticles, Collision, mcCollisions); + break; + default: + PairHistManager.setPairMc(p1, p2, mcParticles, Collision, mcCollisions); } // fill deta-dphi histograms with kstar cutoff CprManager.fill(PairHistManager.getKstar()); @@ -239,111 +241,7 @@ void processSameEvent(T1 const& SliceParticle1, } } -// process mixed event for identical particles -template -void processMixedEvent(T1 const& Collisions, - T2& Partition, - T3 const& TrackTable, - T4& cache, - T5 const& policy, - T6 const& depth, - T7& PairHistManager, - T8& CprManager, - T9& PcManager) -{ - for (auto const& [collision1, collision2] : o2::soa::selfCombinations(policy, depth, -1, Collisions, Collisions)) { - if (!(std::fabs(collision1.magField() - collision2.magField()) < o2::constants::math::Epsilon)) { - continue; - } - CprManager.setMagField(collision1.magField()); - auto sliceParticle1 = Partition->sliceByCached(o2::aod::femtobase::stored::fColId, collision1.globalIndex(), cache); - auto sliceParticle2 = Partition->sliceByCached(o2::aod::femtobase::stored::fColId, collision2.globalIndex(), cache); - if (sliceParticle1.size() == 0 || sliceParticle2.size() == 0) { - continue; - } - for (auto const& [p1, p2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(sliceParticle1, sliceParticle2))) { - // pair cleaning - if (!PcManager.isCleanPair(p1, p2, TrackTable)) { - continue; - } - // Close pair rejection - CprManager.setPair(p1, p2, TrackTable); - if (CprManager.isClosePair()) { - continue; - } - PairHistManager.setPair(p1, p2, collision1, collision2); - CprManager.fill(PairHistManager.getKstar()); - if (PairHistManager.checkPairCuts()) { - PairHistManager.template fill(); - } - } - } -} - -// process mixed event for identical particles with mc information -template -void processMixedEvent(T1 const& Collisions, - T2 const& mcCollisions, - T3& Partition, - T4 const& TrackTable, - T5 const& mcParticles, - T6& cache, - T7 const& policy, - T8 const& depth, - T9& PairHistManager, - T10& CprManager, - T11& PcManager) -{ - for (auto const& [collision1, collision2] : o2::soa::selfCombinations(policy, depth, -1, Collisions, Collisions)) { - if (!(std::fabs(collision1.magField() - collision2.magField()) < o2::constants::math::Epsilon)) { - continue; - } - CprManager.setMagField(collision1.magField()); - auto sliceParticle1 = Partition->sliceByCached(o2::aod::femtobase::stored::fColId, collision1.globalIndex(), cache); - auto sliceParticle2 = Partition->sliceByCached(o2::aod::femtobase::stored::fColId, collision2.globalIndex(), cache); - if (sliceParticle1.size() == 0 || sliceParticle2.size() == 0) { - continue; - } - for (auto const& [p1, p2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(sliceParticle1, sliceParticle2))) { - // pair cleaning - if (!PcManager.isCleanPair(p1, p2, TrackTable)) { - continue; - } - // Close pair rejection - CprManager.setPair(p1, p2, TrackTable); - if (CprManager.isClosePair()) { - continue; - } - PairHistManager.setPairMc(p1, p2, mcParticles, collision1, collision2, mcCollisions); - CprManager.fill(PairHistManager.getKstar()); - if (PairHistManager.checkPairCuts()) { - PairHistManager.template fill(); - } - } - } -} - -// process mixed event different particles +// process mixed event template + typename T12, + typename T11> void processMixedEvent(T1 const& Collisions, T2 const& mcCollisions, T3& Partition1, @@ -416,14 +314,14 @@ void processMixedEvent(T1 const& Collisions, T5 const& TrackTable, T6 const& mcParticles, T7& cache, - T8& policy, - T9& depth, + T8 const& policy, + T9 const& depth, T10& PairHistManager, T11& CprManager, T12& PcManager) { for (auto const& [collision1, collision2] : o2::soa::selfCombinations(policy, depth, -1, Collisions, Collisions)) { - if (!(std::fabs(collision1.magField() - collision2.magField()) < o2::constants::math::Epsilon)) { + if (collision1.magField() != collision2.magField()) { continue; } CprManager.setMagField(collision1.magField()); diff --git a/PWGCF/Femto/Core/trackBuilder.h b/PWGCF/Femto/Core/trackBuilder.h index 7e1c90317bc..b3356b8f1d1 100644 --- a/PWGCF/Femto/Core/trackBuilder.h +++ b/PWGCF/Femto/Core/trackBuilder.h @@ -97,8 +97,8 @@ struct ConfTrackBits : o2::framework::ConfigurableGroup { o2::framework::Configurable> itsProton{"itsProton", {}, "Maximum |nsigma| for Proton PID"}; o2::framework::Configurable> tpcProton{"tpcProton", {}, "Maximum |nsigma| for Proton PID"}; o2::framework::Configurable> tofProton{"tofProton", {}, "Maximum |nsigma| for Proton PID"}; - o2::framework::Configurable> tpcitsProton{"tpcitsProton", {3.f}, "Maximum |nsigma| for Proton PID"}; - o2::framework::Configurable> tpctofProton{"tpctofProton", {3.f}, "Maximum |nsigma| for Proton PID"}; + o2::framework::Configurable> tpcitsProton{"tpcitsProton", {}, "Maximum |nsigma| for Proton PID"}; + o2::framework::Configurable> tpctofProton{"tpctofProton", {}, "Maximum |nsigma| for Proton PID"}; // Deuteron PID cuts o2::framework::Configurable requirePidDeuteron{"requirePidDeuteron", false, "Make election PID optional"}; @@ -592,10 +592,6 @@ class TrackBuilder if (!mFillAnyTable) { return; } - - // clear index map before processing next batch of tracks - indexMap.clear(); - for (const auto& track : tracks) { if (!mTrackSelection.checkFilters(track)) { continue; @@ -673,8 +669,6 @@ class TrackBuilder if (!mFillAnyTable) { return; } - // clear index map before processing next batch of tracks - indexMap.clear(); for (const auto& trackWithItsPid : tracksWithItsPid) { if (!mTrackSelection.checkFilters(trackWithItsPid)) { continue; @@ -683,7 +677,6 @@ class TrackBuilder if (!mTrackSelection.passesAllRequiredSelections()) { continue; } - collisionBuilder.template fillMcCollision(collisionProducts, col, mcCols, mcProducts, mcBuilder); // get track from the track table so we can dereference mc particle properly auto track = tracks.iteratorAt(trackWithItsPid.index()); @@ -728,6 +721,13 @@ class TrackBuilder } } + template + void reset(T const& tracks) + { + indexMap.clear(); + indexMap.reserve(tracks.size()); + }; + private: TrackSelection mTrackSelection; bool mFillAnyTable = false; diff --git a/PWGCF/Femto/Core/trackHistManager.h b/PWGCF/Femto/Core/trackHistManager.h index 4977aa4bd30..be9830dc8f9 100644 --- a/PWGCF/Femto/Core/trackHistManager.h +++ b/PWGCF/Femto/Core/trackHistManager.h @@ -150,6 +150,7 @@ struct ConfTrackBinning : o2::framework::ConfigurableGroup { constexpr const char PrefixTrackBinning1[] = "TrackBinning1"; constexpr const char PrefixTrackBinning2[] = "TrackBinning2"; +constexpr const char PrefixTrackBinning3[] = "TrackBinning3"; constexpr const char PrefixResonancePosDauBinning[] = "ResonancePosDauBinning"; constexpr const char PrefixResonanceNegDauBinning[] = "ResonanceNegDauBinning"; constexpr const char PrefixV0PosDauBinning[] = "V0PosDauBinning"; @@ -161,6 +162,7 @@ constexpr const char PrefixKinkChaDauBinning[] = "KinkChaDauBinning"; using ConfTrackBinning1 = ConfTrackBinning; using ConfTrackBinning2 = ConfTrackBinning; +using ConfTrackBinning3 = ConfTrackBinning; using ConfResonancePosDauBinning = ConfTrackBinning; using ConfResonanceNegDauBinning = ConfTrackBinning; using ConfV0PosDauBinning = ConfTrackBinning; diff --git a/PWGCF/Femto/Core/tripletBuilder.h b/PWGCF/Femto/Core/tripletBuilder.h new file mode 100644 index 00000000000..1487d3766c1 --- /dev/null +++ b/PWGCF/Femto/Core/tripletBuilder.h @@ -0,0 +1,366 @@ +// Copyright 2019-2025 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 tripletBuilder.h +/// \brief histogram manager for pair tasks +/// \author anton.riedel@tum.de, TU München, anton.riedel@tum.de + +#ifndef PWGCF_FEMTO_CORE_TRIPLETBUILDER_H_ +#define PWGCF_FEMTO_CORE_TRIPLETBUILDER_H_ + +#include "PWGCF/Femto/Core/closeTripletRejection.h" +#include "PWGCF/Femto/Core/collisionHistManager.h" +#include "PWGCF/Femto/Core/modes.h" +#include "PWGCF/Femto/Core/trackHistManager.h" +#include "PWGCF/Femto/Core/tripletCleaner.h" +#include "PWGCF/Femto/Core/tripletHistManager.h" +#include "PWGCF/Femto/Core/tripletProcessHelpers.h" +#include "PWGCF/Femto/DataModel/FemtoTables.h" + +#include "Framework/HistogramRegistry.h" +#include "Framework/HistogramSpec.h" + +#include "fairlogger/Logger.h" + +#include +#include +#include +#include +#include + +namespace o2::analysis::femto +{ +namespace tripletbuilder +{ + +template +class TripletTrackTrackTrackBuilder +{ + public: + TripletTrackTrackTrackBuilder() = default; + ~TripletTrackTrackTrackBuilder() = default; + + template + void init(o2::framework::HistogramRegistry* registry, + T1 const& confTrackSelection1, + T2 const& confTrackSelection2, + T3 const& confTrackSelection3, + T4 const& confCtr, + T5 const& confMixing, + T6 const& confTripletBinning, + T7 const& confTripletCuts, + std::map> const& colHistSpec, + std::map> const& trackHistSpec1, + std::map> const& trackHistSpec2, + std::map> const& trackHistSpec3, + std::map> const& pairHistSpec, + std::map> const& cprHistSpec) + { + // check if correlate the same tracks or not + mTrack1Track2Track3AreSameSpecies = confMixing.particle123AreSameSpecies.value; + mTrack1Track2AreSameSpecies = confMixing.particle12AreSameSpecies.value; + + if (mTrack1Track2Track3AreSameSpecies && mTrack1Track2AreSameSpecies) { + LOG(fatal) << "Option Track 1&2 are identical and Option Track 1&2&3 are identical is activated. Breaking..."; + } + + mColHistManager.template init(registry, colHistSpec); + mTripletHistManagerSe.template init(registry, pairHistSpec, confTripletBinning, confTripletCuts); + mTripletHistManagerMe.template init(registry, pairHistSpec, confTripletBinning, confTripletCuts); + + mTc.template init(confTripletCuts); + + if (mTrack1Track2Track3AreSameSpecies) { + // Track1 & Track2 & Track3 are the same particle species + mTrackHistManager1.template init(registry, trackHistSpec1, confTrackSelection1); + + mTripletHistManagerSe.setMass(confTrackSelection1.pdgCodeAbs.value, confTrackSelection1.pdgCodeAbs.value, confTrackSelection1.pdgCodeAbs.value); + mTripletHistManagerSe.setCharge(confTrackSelection1.chargeAbs.value, confTrackSelection1.chargeAbs.value, confTrackSelection1.chargeAbs.value); + mCtrSe.init(registry, cprHistSpec, confCtr, confTrackSelection1.chargeAbs.value, confTrackSelection1.chargeAbs.value, confTrackSelection1.chargeAbs.value); + + mTripletHistManagerMe.setMass(confTrackSelection1.pdgCodeAbs.value, confTrackSelection1.pdgCodeAbs.value, confTrackSelection1.pdgCodeAbs.value); + mTripletHistManagerMe.setCharge(confTrackSelection1.chargeAbs.value, confTrackSelection1.chargeAbs.value, confTrackSelection1.chargeAbs.value); + mCtrMe.init(registry, cprHistSpec, confCtr, confTrackSelection1.chargeAbs.value, confTrackSelection1.chargeAbs.value, confTrackSelection1.chargeAbs.value); + } else if (mTrack1Track2AreSameSpecies) { + // Track1 & Track2 & are the same particle species and track 3 is something else + mTrackHistManager1.template init(registry, trackHistSpec1, confTrackSelection1); + mTrackHistManager3.template init(registry, trackHistSpec3, confTrackSelection2); + + mTripletHistManagerSe.setMass(confTrackSelection1.pdgCodeAbs.value, confTrackSelection1.pdgCodeAbs.value, confTrackSelection3.pdgCodeAbs.value); + mTripletHistManagerSe.setCharge(confTrackSelection1.chargeAbs.value, confTrackSelection1.chargeAbs.value, confTrackSelection3.chargeAbs.value); + mCtrSe.init(registry, cprHistSpec, confCtr, confTrackSelection1.chargeAbs.value, confTrackSelection1.chargeAbs.value, confTrackSelection3.chargeAbs.value); + + mTripletHistManagerMe.setMass(confTrackSelection1.pdgCodeAbs.value, confTrackSelection1.pdgCodeAbs.value, confTrackSelection3.pdgCodeAbs.value); + mTripletHistManagerMe.setCharge(confTrackSelection1.chargeAbs.value, confTrackSelection1.chargeAbs.value, confTrackSelection3.chargeAbs.value); + mCtrMe.init(registry, cprHistSpec, confCtr, confTrackSelection1.chargeAbs.value, confTrackSelection1.chargeAbs.value, confTrackSelection3.chargeAbs.value); + } else { + // all three tracks are different + mTrackHistManager1.template init(registry, trackHistSpec1, confTrackSelection1); + mTrackHistManager2.template init(registry, trackHistSpec2, confTrackSelection2); + mTrackHistManager3.template init(registry, trackHistSpec3, confTrackSelection2); + + mTripletHistManagerSe.setMass(confTrackSelection1.pdgCodeAbs.value, confTrackSelection2.pdgCodeAbs.value, confTrackSelection3.pdgCodeAbs.value); + mTripletHistManagerSe.setCharge(confTrackSelection1.chargeAbs.value, confTrackSelection2.chargeAbs.value, confTrackSelection3.chargeAbs.value); + mCtrSe.init(registry, cprHistSpec, confCtr, confTrackSelection1.chargeAbs.value, confTrackSelection2.chargeAbs.value, confTrackSelection3.chargeAbs.value); + + mTripletHistManagerMe.setMass(confTrackSelection1.pdgCodeAbs.value, confTrackSelection2.pdgCodeAbs.value, confTrackSelection3.pdgCodeAbs.value); + mTripletHistManagerMe.setCharge(confTrackSelection1.chargeAbs.value, confTrackSelection2.chargeAbs.value, confTrackSelection3.chargeAbs.value); + mCtrMe.init(registry, cprHistSpec, confCtr, confTrackSelection1.chargeAbs.value, confTrackSelection2.chargeAbs.value, confTrackSelection3.chargeAbs.value); + } + + // setup mixing + mMixingPolicy = static_cast(confMixing.policy.value); + mMixingDepth = confMixing.depth.value; + + // setup rng if necessary + if (confMixing.seed.value >= 0) { + uint64_t randomSeed = 0; + mMixIdenticalParticles = true; + if (confMixing.seed.value == 0) { + randomSeed = static_cast(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count()); + } else { + randomSeed = static_cast(confMixing.seed.value); + } + mRng = std::mt19937(randomSeed); + if (mTrack1Track2Track3AreSameSpecies) { + mDist = std::uniform_int_distribution<>(tripletprocesshelpers::kOrder123, tripletprocesshelpers::kOrder321); + } + if (mTrack1Track2AreSameSpecies) { + mDist = std::uniform_int_distribution<>(tripletprocesshelpers::kOrder123, tripletprocesshelpers::kOrder213); + } + } + } + + // data + template + void processSameEvent(T1 const& col, T2& trackTable, T3& partition1, T4& partition2, T5& partition3, T6& cache) + { + tripletprocesshelpers::TripletOrder tripletOrder = tripletprocesshelpers::kOrder123; + if (mTrack1Track2Track3AreSameSpecies) { + auto trackSlice1 = partition1->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + if (trackSlice1.size() == 0) { + return; + } + mColHistManager.template fill(col); + mCtrSe.setMagField(col.magField()); + if (mMixIdenticalParticles) { + tripletOrder = static_cast(mDist(mRng)); + } + tripletprocesshelpers::processSameEvent(trackSlice1, trackTable, col, mTrackHistManager1, mTripletHistManagerSe, mCtrSe, mTc, tripletOrder); + } else if (mTrack1Track2AreSameSpecies) { + auto trackSlice1 = partition1->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + auto trackSlice3 = partition3->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + if (trackSlice1.size() == 0 || trackSlice3.size() == 0) { + return; + } + mColHistManager.template fill(col); + mCtrSe.setMagField(col.magField()); + if (mMixIdenticalParticles) { + tripletOrder = static_cast(mDist(mRng)); + } + tripletprocesshelpers::processSameEvent(trackSlice1, trackSlice3, trackTable, col, mTrackHistManager1, mTrackHistManager3, mTripletHistManagerSe, mCtrSe, mTc, tripletOrder); + } else { + auto trackSlice1 = partition1->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + auto trackSlice2 = partition2->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + auto trackSlice3 = partition3->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + if (trackSlice1.size() == 0 || trackSlice2.size() == 0 || trackSlice3.size() == 0) { + return; + } + mColHistManager.template fill(col); + mCtrSe.setMagField(col.magField()); + tripletprocesshelpers::processSameEvent(trackSlice1, trackSlice2, trackSlice3, trackTable, col, mTrackHistManager1, mTrackHistManager2, mTrackHistManager3, mTripletHistManagerSe, mCtrSe, mTc); + } + } + + // mc + template + void processSameEvent(T1 const& col, T2 const& mcCols, T3& trackTable, T4& partition1, T5& partition2, T6& partition3, T7 const& mcParticles, T8 const& mcMothers, T9 const& mcPartonicMothers, T10& cache) + { + tripletprocesshelpers::TripletOrder tripletOrder = tripletprocesshelpers::kOrder123; + if (mTrack1Track2Track3AreSameSpecies) { + auto trackSlice1 = partition1->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + if (trackSlice1.size() == 0) { + return; + } + mColHistManager.template fill(col); + mCtrSe.setMagField(col.magField()); + if (mMixIdenticalParticles) { + tripletOrder = static_cast(mDist(mRng)); + } + tripletprocesshelpers::processSameEvent(trackSlice1, trackTable, mcParticles, mcMothers, mcPartonicMothers, col, mcCols, mTrackHistManager1, mTripletHistManagerSe, mCtrSe, mTc, tripletOrder); + } else if (mTrack1Track2AreSameSpecies) { + auto trackSlice1 = partition1->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + auto trackSlice3 = partition3->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + if (trackSlice1.size() == 0 || trackSlice3.size() == 0) { + return; + } + mColHistManager.template fill(col); + mCtrSe.setMagField(col.magField()); + if (mMixIdenticalParticles) { + tripletOrder = static_cast(mDist(mRng)); + } + tripletprocesshelpers::processSameEvent(trackSlice1, trackSlice3, trackTable, mcParticles, mcMothers, mcPartonicMothers, col, mcCols, mTrackHistManager1, mTrackHistManager3, mTripletHistManagerSe, mCtrSe, mTc, tripletOrder); + } else { + auto trackSlice1 = partition1->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + auto trackSlice2 = partition2->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + auto trackSlice3 = partition3->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + if (trackSlice1.size() == 0 || trackSlice2.size() == 0 || trackSlice3.size() == 0) { + return; + } + mColHistManager.template fill(col); + mCtrSe.setMagField(col.magField()); + tripletprocesshelpers::processSameEvent(trackSlice1, trackSlice2, trackSlice3, trackTable, mcParticles, mcMothers, mcPartonicMothers, col, mcCols, mTrackHistManager1, mTrackHistManager2, mTrackHistManager3, mTripletHistManagerSe, mCtrSe, mTc); + } + } + + template + void processMixedEvent(T1 const& cols, T2& trackTable, T3& partition1, T4& partition2, T5& partition3, T6& cache, T7& binsVtxMult, T8& binsVtxCent, T9& binsVtxMultCent) + { + if (mTrack1Track2Track3AreSameSpecies) { + switch (mMixingPolicy) { + case static_cast(pairhistmanager::kVtxMult): + tripletprocesshelpers::processMixedEvent(cols, partition1, partition1, partition1, trackTable, cache, binsVtxMult, mMixingDepth, mTripletHistManagerMe, mCtrMe, mTc); + break; + case static_cast(pairhistmanager::kVtxCent): + tripletprocesshelpers::processMixedEvent(cols, partition1, partition1, partition1, trackTable, cache, binsVtxCent, mMixingDepth, mTripletHistManagerMe, mCtrMe, mTc); + break; + case static_cast(pairhistmanager::kVtxMultCent): + tripletprocesshelpers::processMixedEvent(cols, partition1, partition1, partition1, trackTable, cache, binsVtxMultCent, mMixingDepth, mTripletHistManagerMe, mCtrMe, mTc); + break; + default: + LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; + } + } else if (mTrack1Track2AreSameSpecies) { + switch (mMixingPolicy) { + case static_cast(pairhistmanager::kVtxMult): + tripletprocesshelpers::processMixedEvent(cols, partition1, partition1, partition3, trackTable, cache, binsVtxMult, mMixingDepth, mTripletHistManagerMe, mCtrMe, mTc); + break; + case static_cast(pairhistmanager::kVtxCent): + tripletprocesshelpers::processMixedEvent(cols, partition1, partition1, partition3, trackTable, cache, binsVtxCent, mMixingDepth, mTripletHistManagerMe, mCtrMe, mTc); + break; + case static_cast(pairhistmanager::kVtxMultCent): + tripletprocesshelpers::processMixedEvent(cols, partition1, partition1, partition3, trackTable, cache, binsVtxMultCent, mMixingDepth, mTripletHistManagerMe, mCtrMe, mTc); + break; + default: + LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; + } + } else { + switch (mMixingPolicy) { + case static_cast(pairhistmanager::kVtxMult): + tripletprocesshelpers::processMixedEvent(cols, partition1, partition2, partition3, trackTable, cache, binsVtxMult, mMixingDepth, mTripletHistManagerMe, mCtrMe, mTc); + break; + case static_cast(pairhistmanager::kVtxCent): + tripletprocesshelpers::processMixedEvent(cols, partition1, partition2, partition3, trackTable, cache, binsVtxCent, mMixingDepth, mTripletHistManagerMe, mCtrMe, mTc); + break; + case static_cast(pairhistmanager::kVtxMultCent): + tripletprocesshelpers::processMixedEvent(cols, partition1, partition2, partition3, trackTable, cache, binsVtxMultCent, mMixingDepth, mTripletHistManagerMe, mCtrMe, mTc); + break; + default: + LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; + } + } + } + + template + void processMixedEvent(T1 const& cols, T2 const& mcCols, T3& trackTable, T4& partition1, T5& partition2, T6& partition3, T7 const& mcParticles, T8& cache, T9& binsVtxMult, T10& binsVtxCent, T11& binsVtxMultCent) + { + if (mTrack1Track2Track3AreSameSpecies) { + switch (mMixingPolicy) { + case static_cast(pairhistmanager::kVtxMult): + tripletprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition1, partition1, trackTable, mcParticles, cache, binsVtxMult, mMixingDepth, mTripletHistManagerMe, mCtrMe, mTc); + break; + case static_cast(pairhistmanager::kVtxCent): + tripletprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition1, partition1, trackTable, mcParticles, cache, binsVtxCent, mMixingDepth, mTripletHistManagerMe, mCtrMe, mTc); + break; + case static_cast(pairhistmanager::kVtxMultCent): + tripletprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition1, partition1, trackTable, mcParticles, cache, binsVtxMultCent, mMixingDepth, mTripletHistManagerMe, mCtrMe, mTc); + break; + default: + LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; + } + } else if (mTrack1Track2AreSameSpecies) { + switch (mMixingPolicy) { + case static_cast(pairhistmanager::kVtxMult): + tripletprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition1, partition3, trackTable, mcParticles, cache, binsVtxMult, mMixingDepth, mTripletHistManagerMe, mCtrMe, mTc); + break; + case static_cast(pairhistmanager::kVtxCent): + tripletprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition1, partition3, trackTable, mcParticles, cache, binsVtxCent, mMixingDepth, mTripletHistManagerMe, mCtrMe, mTc); + break; + case static_cast(pairhistmanager::kVtxMultCent): + tripletprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition1, partition3, trackTable, mcParticles, cache, binsVtxMultCent, mMixingDepth, mTripletHistManagerMe, mCtrMe, mTc); + break; + default: + LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; + } + } else { + switch (mMixingPolicy) { + case static_cast(pairhistmanager::kVtxMult): + tripletprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition2, partition3, trackTable, mcParticles, cache, binsVtxMult, mMixingDepth, mTripletHistManagerMe, mCtrMe, mTc); + break; + case static_cast(pairhistmanager::kVtxCent): + tripletprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition2, partition3, trackTable, mcParticles, cache, binsVtxCent, mMixingDepth, mTripletHistManagerMe, mCtrMe, mTc); + break; + case static_cast(pairhistmanager::kVtxMultCent): + tripletprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition2, partition3, trackTable, mcParticles, cache, binsVtxMultCent, mMixingDepth, mTripletHistManagerMe, mCtrMe, mTc); + break; + default: + LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; + } + } + } + + private: + colhistmanager::CollisionHistManager mColHistManager; + trackhistmanager::TrackHistManager mTrackHistManager1; + trackhistmanager::TrackHistManager mTrackHistManager2; + trackhistmanager::TrackHistManager mTrackHistManager3; + triplethistmanager::TripletHistManager mTripletHistManagerSe; + triplethistmanager::TripletHistManager mTripletHistManagerMe; + + closetripletrejection::CloseTripletRejectionTrackTrackTrack mCtrSe; + closetripletrejection::CloseTripletRejectionTrackTrackTrack mCtrMe; + tripletcleaner::TrackTrackTrackTripletCleaner mTc; + triplethistmanager::MixingPolicy mMixingPolicy = triplethistmanager::MixingPolicy::kVtxMult; + bool mTrack1Track2Track3AreSameSpecies = false; + bool mTrack1Track2AreSameSpecies = false; + int mMixingDepth = 5; + bool mMixIdenticalParticles = false; + std::mt19937 mRng; + std::uniform_int_distribution<> mDist; +}; + +} // namespace tripletbuilder +} // namespace o2::analysis::femto + +#endif // PWGCF_FEMTO_CORE_TRIPLETBUILDER_H_ diff --git a/PWGCF/Femto/Core/tripletCleaner.h b/PWGCF/Femto/Core/tripletCleaner.h new file mode 100644 index 00000000000..fd8ebc71aad --- /dev/null +++ b/PWGCF/Femto/Core/tripletCleaner.h @@ -0,0 +1,64 @@ +// Copyright 2019-2022 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 tripletCleaner.h +/// \brief triplet cleaner class +/// \author anton.riedel@tum.de, TU München, anton.riedel@tum.de + +#ifndef PWGCF_FEMTO_CORE_TRIPLETCLEANER_H_ +#define PWGCF_FEMTO_CORE_TRIPLETCLEANER_H_ + +#include "PWGCF/Femto/Core/pairBuilder.h" + +namespace o2::analysis::femto +{ +namespace tripletcleaner +{ + +class TrackTrackTrackTripletCleaner : public paircleaner::BasePairCleaner +{ + public: + TrackTrackTrackTripletCleaner() = default; + + template + bool isCleanTriplet(T1 const& track1, T2 const& track2, T3 const& track3, T4 const& /*trackTable*/) const + { + return this->isCleanTrackPair(track1, track2) && + this->isCleanTrackPair(track2, track3) && + this->isCleanTrackPair(track1, track3); + } + + template + bool isCleanTriplet(T1 const& track1, T2 const& track2, T3 const& track3, T4 const& trackTable, T5 const& partonicMothers) const + { + if (!this->isCleanTriplet(track1, track2, track3, trackTable)) { + return false; + } + // pair is clean + // no check if we require common or non-common ancestry + if (mMixPairsWithCommonAncestor) { + return this->pairHasCommonAncestor(track1, track2, partonicMothers) && + this->pairHasCommonAncestor(track2, track3, partonicMothers) && + this->pairHasCommonAncestor(track1, track3, partonicMothers); + } + if (mMixPairsWithNonCommonAncestor) { + return this->pairHasNonCommonAncestor(track1, track2, partonicMothers) && + this->pairHasNonCommonAncestor(track2, track3, partonicMothers) && + this->pairHasNonCommonAncestor(track1, track3, partonicMothers); + } + return true; + } +}; + +} // namespace tripletcleaner +} // namespace o2::analysis::femto + +#endif // PWGCF_FEMTO_CORE_TRIPLETCLEANER_H_ diff --git a/PWGCF/Femto/Core/tripletHistManager.h b/PWGCF/Femto/Core/tripletHistManager.h new file mode 100644 index 00000000000..3ad3062c294 --- /dev/null +++ b/PWGCF/Femto/Core/tripletHistManager.h @@ -0,0 +1,588 @@ +// Copyright 2019-2025 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 tripletHistManager.h +/// \brief histogram manager for triplet tasks +/// \author anton.riedel@tum.de, TU München, anton.riedel@tum.de + +#ifndef PWGCF_FEMTO_CORE_TRIPLETHISTMANAGER_H_ +#define PWGCF_FEMTO_CORE_TRIPLETHISTMANAGER_H_ + +#include "PWGCF/Femto/Core/femtoUtils.h" +#include "PWGCF/Femto/Core/histManager.h" +#include "PWGCF/Femto/Core/modes.h" + +#include "Framework/Configurable.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/HistogramSpec.h" + +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace o2::analysis::femto +{ +namespace triplethistmanager +{ +// enum for pair histograms +enum TripletHist { + // standard 1D + kQ3, + kMt, // averate mt of all pairs in the triplet + // standard 2D + kPt1VsQ3, + kPt2VsQ3, + kPt3VsQ3, + kQ3VsMt, + kQ3VsMult, + kQ3VsCent, + // 2D with mass + // kQ3VsMass1, + // kQ3VsMass2, + // kQ3VsMass3, + // higher dimensions + kPt1VsPt2VsPt3, + kQ3VsPt1VsPt2VsPt3, + kQ3VsMtVsMult, + kQ3VsMtVsMultVsCent, + kQ3VsMtVsPt1VsPt2VsPt3VsMult, + kQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent, + // higher dimensions with mass + // mc + kTrueQ3VsQ3, + kTrueMtVsMt, + kTrueMultVsMult, + kTrueCentVsCent, + + kTripletHistogramLast +}; + +enum MixingPolicy { + kVtxMult, + kVtxCent, + kVtxMultCent, + kMixingPolicyLast +}; + +// Mixing configurables +struct ConfMixing : o2::framework::ConfigurableGroup { + std::string prefix = std::string("Mixing"); + o2::framework::ConfigurableAxis multBins{"multBins", {o2::framework::VARIABLE_WIDTH, 0.0f, 4.0f, 8.0f, 12.0f, 16.0f, 20.0f, 24.0f, 28.0f, 32.0f, 36.0f, 40.0f, 44.0f, 48.0f, 52.0f, 56.0f, 60.0f, 64.0f, 68.0f, 72.0f, 76.0f, 80.0f, 84.0f, 88.0f, 92.0f, 96.0f, 100.0f, 200.0f}, "Mixing bins - multiplicity"}; + o2::framework::ConfigurableAxis centBins{"centBins", {o2::framework::VARIABLE_WIDTH, 0.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.0f}, "Mixing bins - centrality"}; + o2::framework::ConfigurableAxis vtxBins{"vtxBins", {o2::framework::VARIABLE_WIDTH, -10.0f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; + o2::framework::Configurable depth{"depth", 5, "Number of events for mixing"}; + o2::framework::Configurable policy{"policy", 0, "Binning policy for mixing (alywas in combination with z-vertex) -> 0: multiplicity, -> 1: centrality, -> 2: both"}; + o2::framework::Configurable seed{"seed", -1, "Seed to randomize particle 1/2/3 (if they are identical). Set to negative value to deactivate. Set to 0 to generate unique seed in time."}; + o2::framework::Configurable particle123AreSameSpecies{"particle123AreSameSpecies", false, "Particle 1,2 and 3 are of the same species"}; + o2::framework::Configurable particle12AreSameSpecies{"particle12AreSameSpecies", false, "Particle 1 and 2 are of the same species"}; +}; + +struct ConfTripletBinning : o2::framework::ConfigurableGroup { + std::string prefix = std::string("TripletBinning"); + o2::framework::Configurable plot1D{"plot1D", true, "Enable 1D histograms"}; + o2::framework::Configurable plot2D{"plot2D", true, "Enable 2D histograms"}; + o2::framework::Configurable plotPt1VsPt2VsPt3{"plotPt1VsPt2VsPt3", false, "Enable 3D histogram (Pt1 vs Pt2 vs Pt3)"}; + o2::framework::Configurable plotQ3VsPt1VsPt2VsPt3{"plotQ3VsPt1VsPt2VsPt3", false, "Enable 4D histogram (Q3 vs Pt1 vs Pt2 vs Pt3)"}; + o2::framework::Configurable plotQ3VsMtVsMult{"plotQ3VsMtVsMult", false, "Enable 3D histogram (Q3 vs Mt vs Mult)"}; + o2::framework::Configurable plotQ3VsMtVsMultVsCent{"plotQ3VsMtVsMultVsCent", false, "Enable 4D histogram (Q3 vs Mt vs Mult Vs Cent)"}; + o2::framework::Configurable plotQ3VsMtVsPt1VsPt2VsPt3VsMult{"plotQ3VsMtVsPt1VsPt2VsPt3VsMult", false, "Enable 6D histogram (Q3 vs Mt Vs Pt1 vs Pt2 vs Pt3 vs Mult)"}; + o2::framework::Configurable plotQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent{"plotQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent", false, "Enable 3D histogram (Q3 vs Mt Vs Pt1 vs Pt2 vs Pt3 vs Mult vs Cent)"}; + o2::framework::ConfigurableAxis q3{"q3", {{600, 0, 6}}, "q3"}; + o2::framework::ConfigurableAxis mt{"mt", {{500, 0.8, 5.8}}, "mt"}; + o2::framework::ConfigurableAxis multiplicity{"multiplicity", {{50, 0, 200}}, "multiplicity"}; + o2::framework::ConfigurableAxis centrality{"centrality", {{10, 0, 100}}, "centrality (mult. percentile)"}; + o2::framework::ConfigurableAxis pt1{"pt1", {{100, 0, 6}}, "Pt binning for particle 1"}; + o2::framework::ConfigurableAxis pt2{"pt2", {{100, 0, 6}}, "Pt binning for particle 2"}; + o2::framework::ConfigurableAxis pt3{"pt3", {{100, 0, 6}}, "Pt binning for particle 3"}; + o2::framework::ConfigurableAxis mass1{"mass1", {{100, 0, 2}}, "Mass binning for particle 1 (if particle has mass getter)"}; + o2::framework::ConfigurableAxis mass2{"mass2", {{100, 0, 2}}, "Mass binning for particle 2 (if particle has mass getter)"}; + o2::framework::ConfigurableAxis mass3{"mass3", {{100, 0, 2}}, "Mass binning for particle 3 (if particle has mass getter)"}; + o2::framework::Configurable transverseMassType{"transverseMassType", static_cast(modes::TransverseMassType::kAveragePdgMass), "Type of transverse mass (0-> Average Pdg Mass, 1-> Reduced Pdg Mass, 2-> Mt from combined 4 vector)"}; +}; + +struct ConfTripletCuts : o2::framework::ConfigurableGroup { + std::string prefix = std::string("TripletCuts"); + o2::framework::Configurable q3Max{"q3Max", -1, "Maximal kstar (set to -1 to deactivate)"}; + o2::framework::Configurable q3Min{"q3Min", -1, "Minimal kstar (set to -1 to deactivate)"}; + o2::framework::Configurable mtMax{"mtMax", -1, "Maximal mt (set to -1 to deactivate)"}; + o2::framework::Configurable mtMin{"mtMin", -1, "Minimal mt (set to -1 to deactivate)"}; + o2::framework::Configurable mixOnlyCommonAncestor{"mixOnlyCommonAncestor", false, "Require pair to have common anchestor (in the same event)"}; + o2::framework::Configurable mixOnlyNonCommonAncestor{"mixOnlyNonCommonAncestor", false, "Require pair to have non-common anchestor (in the same event)"}; +}; + +// the enum gives the correct index in the array +constexpr std::array, kTripletHistogramLast> + HistTable = { + { + // 1D + {kQ3, o2::framework::kTH1F, "hQ3", "Q_{3}; Q_{3} (GeV/#it{c}); Entries"}, + {kMt, o2::framework::kTH1F, "hMt", "transverse mass; m_{T} (GeV/#it{c}^{2}); Entries"}, + // 2D + {kPt1VsQ3, o2::framework::kTH2F, "hPt1VsQ3", " p_{T,1} vs Q_{3}; p_{T,1} (GeV/#it{c}); Q_{3} (GeV/#it{c})"}, + {kPt2VsQ3, o2::framework::kTH2F, "hPt2VsQ3", "p_{T,2} vs Q_{3}; p_{T,2} (GeV/#it{c}); Q_{3} (GeV/#it{c})"}, + {kPt3VsQ3, o2::framework::kTH2F, "hPt3VsQ3", "p_{T,3} vs Q_{3}; p_{T,2} (GeV/#it{c}); Q_{3} (GeV/#it{c})"}, + {kQ3VsMt, o2::framework::kTH2F, "hQ3VsMt", "Q_{3} vs m_{T}; Q_{3} (GeV/#it{c}); m_{T} (GeV/#it{c}^{2})"}, + {kQ3VsMult, o2::framework::kTH2F, "hQ3VsMult", "Q_{3} vs Multiplicity; Q_{3} (GeV/#it{c}); Multiplicity"}, + {kQ3VsCent, o2::framework::kTH2F, "hQ3VsCent", "Q_{3} vs Centrality; Q_{3} (GeV/#it{c}); Centrality"}, + // n-D + {kPt1VsPt2VsPt3, o2::framework::kTHnSparseF, "hPt1VsPt2VsPt3", "k* vs m_{T} vs multiplicity; k* (GeV/#it{c}); m_{T} (GeV/#it{c}^{2}); Multiplicity"}, + {kQ3VsPt1VsPt2VsPt3, o2::framework::kTHnSparseF, "hQ3VsPt1VsPt2VsPt3", "Q_{3} vs p_{T,1} vs p_{T,2} vs p_{T,3}; Q_{3} (GeV/#it{c}); p_{T,1} (GeV/#it{c}) ; p_{T,2} (GeV/#it{c}); p_{T,3} (GeV/#it{c})"}, + {kQ3VsMtVsMult, o2::framework::kTHnSparseF, "hQ3VsMtVsMult", "Q_{3} vs m_{T} vs multiplicity; Q_{3} (GeV/#it{c}); m_{T} (GeV/#it{c}^{2}); Multiplicity"}, + {kQ3VsMtVsMultVsCent, o2::framework::kTHnSparseF, "hQ3VsMtVsMultVsCent", "Q_{3} vs m_{T} vs multiplicity vs centrality; Q_{3} (GeV/#it{c}); m_{T} (GeV/#it{c}^{2}); Multiplicity; Centrality"}, + // n-D with mass + {kQ3VsMtVsPt1VsPt2VsPt3VsMult, o2::framework::kTHnSparseF, "hQ3VsMtVsPt1VsPt2VsPt3VsMult", "Q_{3} vs m_{T} vs p_{T,1} vs p_{T,2} vs p_{T,3} vs multiplicity; Q_{3} (GeV/#it{c}); m_{T} (GeV/#it{c}^{2}); p_{T,1} (GeV/#it{c}) ; p_{T,2} (GeV/#it{c}); p_{T,3} (GeV/#it{c}) ; Multiplicity"}, + {kQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent, o2::framework::kTHnSparseF, "hQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent", "Q_{3} vs m_{T} vs p_{T,1} vs p_{T,2} vs p_{T,3} vs multiplicity vs centrality; Q_{3} (GeV/#it{c}); m_{T} (GeV/#it{c}^{2}); p_{T,1} (GeV/#it{c}) ; p_{T,2} (GeV/#it{c}); p_{T,3} (GeV/#it{c}) vs multiplicity vs centrality"}, + {kTrueQ3VsQ3, o2::framework::kTH2F, "hTrueQ3VsQ3", "Q_{3,True} vs Q_{3}; Q_{3,true} (GeV/#it{c}); Q_{3} (GeV/#it{c})"}, + {kTrueMtVsMt, o2::framework::kTH2F, "hTrueMtVsMt", "m_{T,True} vs m_{T}; m_{T,True} (GeV/#it{c}^{2}); m_{T} (GeV/#it{c}^{2})"}, + {kTrueMultVsMult, o2::framework::kTH2F, "hTrueMultVsMult", "Multiplicity_{True} vs Multiplicity; Multiplicity_{True} ; Multiplicity"}, + {kTrueCentVsCent, o2::framework::kTH2F, "hTrueCentVsCent", "Centrality_{True} vs Centrality; Centrality_{True} (%); Centrality (%)"}, + }}; + +#define TRIPLET_HIST_ANALYSIS_MAP(conf) \ + {kQ3, {conf.q3}}, \ + {kMt, {conf.mt}}, \ + {kPt1VsQ3, {conf.pt1, conf.q3}}, \ + {kPt2VsQ3, {conf.pt2, conf.q3}}, \ + {kPt3VsQ3, {conf.pt3, conf.q3}}, \ + {kQ3VsMt, {conf.q3, conf.mt}}, \ + {kQ3VsMult, {conf.q3, conf.multiplicity}}, \ + {kQ3VsCent, {conf.q3, conf.centrality}}, \ + {kPt1VsPt2VsPt3, {conf.pt1, conf.pt2, conf.pt3}}, \ + {kQ3VsPt1VsPt2VsPt3, {conf.q3, conf.pt1, conf.pt2, conf.pt3}}, \ + {kQ3VsMtVsMult, {conf.q3, conf.mt, conf.multiplicity}}, \ + {kQ3VsMtVsMultVsCent, {conf.q3, conf.mt, conf.multiplicity, conf.centrality}}, \ + {kQ3VsMtVsPt1VsPt2VsPt3VsMult, {conf.q3, conf.mt, conf.pt1, conf.pt2, conf.pt3, conf.multiplicity}}, \ + {kQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent, {conf.q3, conf.mt, conf.pt1, conf.pt2, conf.pt3, conf.multiplicity, conf.centrality}}, + +#define TRIPLET_HIST_MC_MAP(conf) \ + {kTrueQ3VsQ3, {conf.q3, conf.q3}}, \ + {kTrueMtVsMt, {conf.mt, conf.mt}}, \ + {kTrueMultVsMult, {conf.multiplicity, conf.multiplicity}}, \ + {kTrueCentVsCent, {conf.centrality, conf.centrality}}, + +template +auto makeTripletHistSpecMap(const T& confPairBinning) +{ + return std::map>{ + TRIPLET_HIST_ANALYSIS_MAP(confPairBinning)}; +}; + +template +auto makeTripletMcHistSpecMap(const T& confPairBinning) +{ + return std::map>{ + TRIPLET_HIST_ANALYSIS_MAP(confPairBinning) + TRIPLET_HIST_MC_MAP(confPairBinning)}; +}; + +#undef TRIPLET_HIST_ANALYSIS_MAP +#undef TRIPLET_HIST_MC_MAP + +constexpr char PrefixTrackTrackTrackSe[] = "TrackTrackTrack/SE/"; +constexpr char PrefixTrackTrackTrackMe[] = "TrackTrackTrack/ME/"; + +constexpr std::string_view AnalysisDir = "Analysis/"; +constexpr std::string_view QaDir = "QA/"; +constexpr std::string_view McDir = "MC/"; + +template +class TripletHistManager +{ + public: + TripletHistManager() = default; + ~TripletHistManager() = default; + + template + void init(o2::framework::HistogramRegistry* registry, + std::map> const& Specs, + T1 const& ConfTripletBinning, + T2 const& ConfTripletCuts) + { + mHistogramRegistry = registry; + + // flags for histograms + mPlot1d = ConfTripletBinning.plot1D.value; + mPlot2d = ConfTripletBinning.plot2D.value; + + // transverse mass type + mMtType = static_cast(ConfTripletBinning.transverseMassType.value); + + // values for cuts + mQ3Min = ConfTripletCuts.q3Min.value; + mQ3Max = ConfTripletCuts.q3Max.value; + mMtMin = ConfTripletCuts.mtMin.value; + mMtMax = ConfTripletCuts.mtMax.value; + + mPlotPt1VsPt2VsPt3 = ConfTripletBinning.plotPt1VsPt2VsPt3.value; + mPlotQ3VsPt1VsPt2VsPt3 = ConfTripletBinning.plotQ3VsPt1VsPt2VsPt3.value; + mPlotQ3VsMtVsMult = ConfTripletBinning.plotQ3VsMtVsMult.value; + mPlotQ3VsMtVsMultVsCent = ConfTripletBinning.plotQ3VsMtVsMultVsCent.value; + mPlotQ3VsMtVsPt1VsPt2VsPt3VsMult = ConfTripletBinning.plotQ3VsMtVsPt1VsPt2VsPt3VsMult.value; + mPlotQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent = ConfTripletBinning.plotQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent.value; + + if constexpr (isFlagSet(mode, modes::Mode::kAnalysis)) { + initAnalysis(Specs); + } + + if constexpr (isFlagSet(mode, modes::Mode::kMc)) { + initMc(Specs); + } + } + + void setMass(int PdgParticle1, int PdgParticle2, int PdgParticle3) + { + mPdgMass1 = o2::analysis::femto::utils::getMass(PdgParticle1); + mPdgMass2 = o2::analysis::femto::utils::getMass(PdgParticle2); + mPdgMass3 = o2::analysis::femto::utils::getMass(PdgParticle3); + } + void setCharge(int chargeAbsParticle1, int chargeAbsParticle2, int chargeAbsParticle3) + { + // the pt stored is actually as pt/z for tracks, so in case of particles with z > 1, we have to rescale the pt (this is so far only for He3 the case) + mAbsCharge1 = std::abs(chargeAbsParticle1); + mAbsCharge2 = std::abs(chargeAbsParticle2); + mAbsCharge3 = std::abs(chargeAbsParticle3); + } + + template + void setTriplet(const T1& particle1, const T2& particle2, const T3& particle3) + { + // pt in track table is calculated from 1/signedPt from the original track table + // in case of He with Z=2, we have to rescale the pt with the absolute charge + mParticle1 = ROOT::Math::PtEtaPhiMVector(mAbsCharge1 * particle1.pt(), particle1.eta(), particle1.phi(), mPdgMass1); + mParticle2 = ROOT::Math::PtEtaPhiMVector(mAbsCharge2 * particle2.pt(), particle2.eta(), particle2.phi(), mPdgMass2); + mParticle3 = ROOT::Math::PtEtaPhiMVector(mAbsCharge3 * particle3.pt(), particle3.eta(), particle3.phi(), mPdgMass3); + + // set mT + mMt = getMt(mParticle1, mPdgMass1, mParticle2, mPdgMass2, mParticle3, mPdgMass3); + // set Q3 + mQ3 = getQ3(mParticle1, mPdgMass1, mParticle2, mPdgMass2, mParticle3, mPdgMass3); + + // if one of the particles has a mass getter, we cache the value for the filling later + if constexpr (modes::hasMass(particleType1)) { + mMass1 = particle1.mass(); + } + if constexpr (modes::hasMass(particleType2)) { + mMass2 = particle2.mass(); + } + if constexpr (modes::hasMass(particleType3)) { + mMass3 = particle2.mass(); + } + } + + template + void setTriplet(T1 const& particle1, T2 const& particle2, T3 const& particle3, T4 const& col) + { + setTriplet(particle1, particle2, particle3); + mMult = col.mult(); + mCent = col.cent(); + } + + template + void setTriplet(T1 const& particle1, T2 const& particle2, T3 const& particle3, T4 const& col1, T4 const& col2, T4 const& col3) + { + setTriplet(particle1, particle2, particle3); + mMult = (col1.mult() + col2.mult() + col3.mult()) / 3.f; // if mixing with multiplicity, should be in the same mixing bin + mCent = (col1.cent() + col2.cent() + col3.cent()) / 3.f; // if mixing with centrality, should be in the same mixing bin + } + + template + void setTripletMc(T1 const& particle1, T2 const& particle2, T3 const& particle3, const T4& /*mcParticles*/) + { + if (!particle1.has_fMcParticle() || !particle2.has_fMcParticle() || !particle3.has_fMcParticle()) { + mHasMcTriplet = false; + return; + } + mHasMcTriplet = true; + + auto mcParticle1 = particle1.template fMcParticle_as(); + auto mcParticle2 = particle2.template fMcParticle_as(); + auto mcParticle3 = particle3.template fMcParticle_as(); + + mTrueParticle1 = ROOT::Math::PtEtaPhiMVector(mAbsCharge1 * mcParticle1.pt(), mcParticle1.eta(), mcParticle1.phi(), mPdgMass1); + mTrueParticle2 = ROOT::Math::PtEtaPhiMVector(mAbsCharge2 * mcParticle2.pt(), mcParticle2.eta(), mcParticle2.phi(), mPdgMass2); + mTrueParticle3 = ROOT::Math::PtEtaPhiMVector(mAbsCharge3 * mcParticle3.pt(), mcParticle3.eta(), mcParticle3.phi(), mPdgMass3); + + // set true mT + mTrueMt = getMt(mTrueParticle1, mPdgMass1, mTrueParticle2, mPdgMass2, mTrueParticle3, mPdgMass3); + // set true Q3 + mTrueQ3 = getQ3(mTrueParticle1, mPdgMass1, mTrueParticle2, mPdgMass2, mTrueParticle3, mPdgMass3); + } + + template + void setTripletMc(T1 const& particle1, T2 const& particle2, T3 const& particle3, T4 const& mcParticles, T5 const& col, const T6& /*mcCols*/) + { + setTriplet(particle1, particle2, particle3, col); + setTripletMc(particle1, particle2, particle3, mcParticles); + if (!col.has_fMcCol()) { + mHasMcCol = false; + return; + } + mHasMcCol = true; + auto mcCol = col.template fMcCol_as(); + mTrueMult = mcCol.mult(); + mTrueCent = mcCol.cent(); + } + + template + void setTripletMc(T1 const& particle1, T2 const& particle2, T3 const& particle3, T4 const& mcParticles, T5 const& col1, T5 const& col2, T5 const& col3, const T6& /*mcCols*/) + { + setTriplet(particle1, particle2, particle3, col1, col2, col3); + setTripletMc(particle1, particle2, particle3, mcParticles); + if (!col1.has_fMcCol() || !col2.has_fMcCol() || !col2.has_fMcCol()) { + mHasMcCol = false; + return; + } + mHasMcCol = true; + auto mcCol1 = col1.template fMcCol_as(); + auto mcCol2 = col2.template fMcCol_as(); + auto mcCol3 = col2.template fMcCol_as(); + mTrueMult = (mcCol1.mult() + mcCol2.mult() + mcCol3.mult()) / 3.f; + mTrueCent = (mcCol1.cent() + mcCol2.cent() + mcCol3.cent()) / 3.f; + } + + bool checkTripletCuts() const + { + return (!(mQ3Min > 0.f) || mQ3 > mQ3Min) && + (!(mQ3Max > 0.f) || mQ3 < mQ3Max) && + (!(mMtMin > 0.f) || mMt > mMtMin) && + (!(mMtMax > 0.f) || mMt < mMtMax); + } + + template + void fill() + { + if constexpr (isFlagSet(mode, modes::Mode::kAnalysis)) { + fillAnalysis(); + } + if constexpr (isFlagSet(mode, modes::Mode::kMc)) { + fillMc(); + } + } + + float getQ3() const { return mQ3; } + + private: + template + T getqij(const T& vecparti, const T& vecpartj) + { + T trackSum = vecparti + vecpartj; + T trackDifference = vecparti - vecpartj; + float scaling = trackDifference.Dot(trackSum) / trackSum.Dot(trackSum); + return trackDifference - scaling * trackSum; + } + + template + float getQ3(const T& part1, const float mass1, const T& part2, const float mass2, const T& part3, const float mass3) + { + ROOT::Math::PtEtaPhiMVector vecpart1(part1.pt(), part1.eta(), part1.phi(), mass1); + ROOT::Math::PtEtaPhiMVector vecpart2(part2.pt(), part2.eta(), part2.phi(), mass2); + ROOT::Math::PtEtaPhiMVector vecpart3(part3.pt(), part3.eta(), part3.phi(), mass3); + auto q12 = getqij(vecpart1, vecpart2); + auto q23 = getqij(vecpart2, vecpart3); + auto q31 = getqij(vecpart3, vecpart1); + float q = q12.M2() + q23.M2() + q31.M2(); + return std::sqrt(-q); + } + + void initAnalysis(std::map> const& Specs) + { + std::string analysisDir = std::string(prefix) + std::string(AnalysisDir); + if (mPlot1d) { + mHistogramRegistry->add(analysisDir + getHistNameV2(kQ3, HistTable), getHistDesc(kQ3, HistTable), getHistType(kQ3, HistTable), {Specs.at(kQ3)}); + mHistogramRegistry->add(analysisDir + getHistNameV2(kMt, HistTable), getHistDesc(kMt, HistTable), getHistType(kMt, HistTable), {Specs.at(kMt)}); + } + if (mPlot2d) { + mHistogramRegistry->add(analysisDir + getHistNameV2(kPt1VsQ3, HistTable), getHistDesc(kPt1VsQ3, HistTable), getHistType(kPt1VsQ3, HistTable), {Specs.at(kPt1VsQ3)}); + mHistogramRegistry->add(analysisDir + getHistNameV2(kPt2VsQ3, HistTable), getHistDesc(kPt2VsQ3, HistTable), getHistType(kPt2VsQ3, HistTable), {Specs.at(kPt2VsQ3)}); + mHistogramRegistry->add(analysisDir + getHistNameV2(kPt3VsQ3, HistTable), getHistDesc(kPt3VsQ3, HistTable), getHistType(kPt3VsQ3, HistTable), {Specs.at(kPt3VsQ3)}); + mHistogramRegistry->add(analysisDir + getHistNameV2(kQ3VsMt, HistTable), getHistDesc(kQ3VsMt, HistTable), getHistType(kQ3VsMt, HistTable), {Specs.at(kQ3VsMt)}); + mHistogramRegistry->add(analysisDir + getHistNameV2(kQ3VsMult, HistTable), getHistDesc(kQ3VsMult, HistTable), getHistType(kQ3VsMult, HistTable), {Specs.at(kQ3VsMult)}); + mHistogramRegistry->add(analysisDir + getHistNameV2(kQ3VsCent, HistTable), getHistDesc(kQ3VsCent, HistTable), getHistType(kQ3VsCent, HistTable), {Specs.at(kQ3VsCent)}); + // TODO: implement histograms differential im mass of the particles + } + + if (mPlotPt1VsPt2VsPt3) { + mHistogramRegistry->add(analysisDir + getHistNameV2(kPt1VsPt2VsPt3, HistTable), getHistDesc(kPt1VsPt2VsPt3, HistTable), getHistType(kPt1VsPt2VsPt3, HistTable), {Specs.at(kPt1VsPt2VsPt3)}); + } + if (mPlotQ3VsPt1VsPt2VsPt3) { + mHistogramRegistry->add(analysisDir + getHistNameV2(kQ3VsPt1VsPt2VsPt3, HistTable), getHistDesc(kQ3VsPt1VsPt2VsPt3, HistTable), getHistType(kQ3VsPt1VsPt2VsPt3, HistTable), {Specs.at(kQ3VsPt1VsPt2VsPt3)}); + } + if (mPlotQ3VsMtVsMult) { + mHistogramRegistry->add(analysisDir + getHistNameV2(kQ3VsMtVsMult, HistTable), getHistDesc(kQ3VsMtVsMult, HistTable), getHistType(kQ3VsMtVsMult, HistTable), {Specs.at(kQ3VsMtVsMult)}); + } + if (mPlotQ3VsMtVsMultVsCent) { + mHistogramRegistry->add(analysisDir + getHistNameV2(kQ3VsMtVsMultVsCent, HistTable), getHistDesc(kQ3VsMtVsMultVsCent, HistTable), getHistType(kQ3VsMtVsMultVsCent, HistTable), {Specs.at(kQ3VsMtVsMultVsCent)}); + } + if (mPlotQ3VsMtVsPt1VsPt2VsPt3VsMult) { + mHistogramRegistry->add(analysisDir + getHistNameV2(kQ3VsMtVsPt1VsPt2VsPt3VsMult, HistTable), getHistDesc(kQ3VsMtVsPt1VsPt2VsPt3VsMult, HistTable), getHistType(kQ3VsMtVsPt1VsPt2VsPt3VsMult, HistTable), {Specs.at(kQ3VsMtVsPt1VsPt2VsPt3VsMult)}); + } + if (mPlotQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent) { + mHistogramRegistry->add(analysisDir + getHistNameV2(kQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent, HistTable), getHistDesc(kQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent, HistTable), getHistType(kQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent, HistTable), {Specs.at(kQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent)}); + } + } + + void initMc(std::map> const& Specs) + { + std::string mcDir = std::string(prefix) + std::string(McDir); + mHistogramRegistry->add(mcDir + getHistNameV2(kTrueQ3VsQ3, HistTable), getHistDesc(kTrueQ3VsQ3, HistTable), getHistType(kTrueQ3VsQ3, HistTable), {Specs.at(kTrueQ3VsQ3)}); + mHistogramRegistry->add(mcDir + getHistNameV2(kTrueMtVsMt, HistTable), getHistDesc(kTrueMtVsMt, HistTable), getHistType(kTrueMtVsMt, HistTable), {Specs.at(kTrueMtVsMt)}); + mHistogramRegistry->add(mcDir + getHistNameV2(kTrueMultVsMult, HistTable), getHistDesc(kTrueMultVsMult, HistTable), getHistType(kTrueMultVsMult, HistTable), {Specs.at(kTrueMultVsMult)}); + mHistogramRegistry->add(mcDir + getHistNameV2(kTrueCentVsCent, HistTable), getHistDesc(kTrueCentVsCent, HistTable), getHistType(kTrueCentVsCent, HistTable), {Specs.at(kTrueCentVsCent)}); + } + + void fillAnalysis() + { + if (mPlot1d) { + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQ3, HistTable)), mQ3); + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kMt, HistTable)), mMt); + } + if (mPlot2d) { + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kPt1VsQ3, HistTable)), mParticle1.Pt(), mQ3); + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kPt2VsQ3, HistTable)), mParticle2.Pt(), mQ3); + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kPt3VsQ3, HistTable)), mParticle3.Pt(), mQ3); + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQ3VsMt, HistTable)), mQ3, mMt); + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQ3VsMult, HistTable)), mQ3, mMult); + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQ3VsCent, HistTable)), mQ3, mCent); + } + + if (mPlotPt1VsPt2VsPt3) { + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kPt1VsPt2VsPt3, HistTable)), mParticle1.Pt(), mParticle2.Pt(), mParticle3.Pt()); + } + if (mPlotQ3VsPt1VsPt2VsPt3) { + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQ3VsPt1VsPt2VsPt3, HistTable)), mQ3, mParticle1.Pt(), mParticle2.Pt(), mParticle3.Pt()); + } + if (mPlotQ3VsMtVsMult) { + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQ3VsMtVsMult, HistTable)), mQ3, mMt, mMult); + } + if (mPlotQ3VsMtVsMultVsCent) { + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQ3VsMtVsMultVsCent, HistTable)), mQ3, mMt, mMult, mCent); + } + if (mPlotQ3VsMtVsPt1VsPt2VsPt3VsMult) { + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQ3VsMtVsPt1VsPt2VsPt3VsMult, HistTable)), mQ3, mMt, mParticle1.Pt(), mParticle2.Pt(), mParticle3.Pt(), mMult); + } + if (mPlotQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent) { + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent, HistTable)), mQ3, mMt, mParticle1.Pt(), mParticle2.Pt(), mParticle3.Pt(), mMult, mCent); + } + } + + void fillMc() + { + if (mHasMcTriplet) { + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTrueQ3VsQ3, HistTable)), mTrueQ3, mQ3); + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTrueMtVsMt, HistTable)), mTrueMt, mMt); + } + if (mHasMcCol) { + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTrueMultVsMult, HistTable)), mTrueMult, mMult); + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTrueCentVsCent, HistTable)), mTrueCent, mCent); + } + } + + float getMt(ROOT::Math::PtEtaPhiMVector const& part1, + float pdgMass1, + ROOT::Math::PtEtaPhiMVector const& part2, + float pdgMass2) + { + auto sum = part1 + part2; + float mt = 0; + float averageMass = 0; + float reducedMass = 0; + switch (mMtType) { + case modes::TransverseMassType::kAveragePdgMass: + averageMass = 0.5 * (pdgMass1 + pdgMass2); + mt = std::hypot(0.5 * sum.Pt(), averageMass); + break; + case modes::TransverseMassType::kReducedPdgMass: + reducedMass = 2.f * (pdgMass1 * pdgMass2 / (pdgMass1 + pdgMass2)); + mt = std::hypot(0.5 * sum.Pt(), reducedMass); + break; + case modes::TransverseMassType::kMt4Vector: + mt = 0.5 * sum.Mt(); + break; + default: + LOG(fatal) << "Invalid transverse mass type, breaking..."; + } + return mt; + } + + float getMt(ROOT::Math::PtEtaPhiMVector const& part1, + float mass1, + ROOT::Math::PtEtaPhiMVector const& part2, + float mass2, + ROOT::Math::PtEtaPhiMVector const& part3, + float mass3) + { + return (getMt(part1, mass1, part2, mass2) + getMt(part2, mass2, part3, mass3) + getMt(part1, mass1, part3, mass3)) / 3.; + } + + o2::framework::HistogramRegistry* mHistogramRegistry = nullptr; + float mPdgMass1 = 0.f; + float mPdgMass2 = 0.f; + float mPdgMass3 = 0.f; + + modes::TransverseMassType mMtType = modes::TransverseMassType::kAveragePdgMass; + + int mAbsCharge1 = 1; + int mAbsCharge2 = 1; + int mAbsCharge3 = 1; + ROOT::Math::PtEtaPhiMVector mParticle1{}; + ROOT::Math::PtEtaPhiMVector mParticle2{}; + ROOT::Math::PtEtaPhiMVector mParticle3{}; + float mMass1 = 0.f; + float mMass2 = 0.f; + float mMass3 = 0.f; + float mQ3 = 0.f; + float mMt = 0.f; + float mMult = 0.f; + float mCent = 0.f; + + // mc + ROOT::Math::PtEtaPhiMVector mTrueParticle1{}; + ROOT::Math::PtEtaPhiMVector mTrueParticle2{}; + ROOT::Math::PtEtaPhiMVector mTrueParticle3{}; + float mTrueQ3 = 0.f; + float mTrueMt = 0.f; + float mTrueMult = 0.f; + float mTrueCent = 0.f; + + // cuts + bool mHasMcTriplet = false; + bool mHasMcCol = false; + float mQ3Min = -1.f; + float mQ3Max = -1.f; + float mKtMin = -1.f; + float mKtMax = -1.f; + float mMtMin = -1.f; + float mMtMax = -1.f; + + // flags + bool mPlot1d = true; + bool mPlot2d = true; + + bool mPlotPt1VsPt2VsPt3 = false; + bool mPlotQ3VsPt1VsPt2VsPt3 = false; + bool mPlotQ3VsMtVsMult = false; + bool mPlotQ3VsMtVsMultVsCent = false; + bool mPlotQ3VsMtVsPt1VsPt2VsPt3VsMult = false; + bool mPlotQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent = false; +}; + +}; // namespace triplethistmanager +}; // namespace o2::analysis::femto +#endif // PWGCF_FEMTO_CORE_TRIPLETHISTMANAGER_H_ diff --git a/PWGCF/Femto/Core/tripletProcessHelpers.h b/PWGCF/Femto/Core/tripletProcessHelpers.h new file mode 100644 index 00000000000..fc52c62da0e --- /dev/null +++ b/PWGCF/Femto/Core/tripletProcessHelpers.h @@ -0,0 +1,541 @@ +// Copyright 2019-2025 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 tripletProcessHelpers.h +/// \brief process functions used in pair tasks +/// \author anton.riedel@tum.de, TU München, anton.riedel@tum.de + +#ifndef PWGCF_FEMTO_CORE_TRIPLETPROCESSHELPERS_H_ +#define PWGCF_FEMTO_CORE_TRIPLETPROCESSHELPERS_H_ + +#include "PWGCF/Femto/Core/modes.h" +#include "PWGCF/Femto/DataModel/FemtoTables.h" + +#include "Framework/ASoAHelpers.h" + +namespace o2::analysis::femto +{ +namespace tripletprocesshelpers +{ + +enum TripletOrder : uint8_t { + kOrder123, // no swap + kOrder213, // first two swap 1&2 so we can use them for the case that particle 1 & 2 are the same species, particle 3 is something else + kOrder132, // swap 2&3 + kOrder321, // swap 1&2&3 +}; + +// process same event for identical 3 particles +template +void processSameEvent(T1 const& SliceParticle, + T2 const& TrackTable, + T3 const& Collision, + T4& ParticleHistManager, + T5& TripletHistManager, + T6& CtrManager, + T7& TcManager, + TripletOrder tripletOrder) +{ + for (auto const& part : SliceParticle) { + ParticleHistManager.template fill(part, TrackTable); + } + + for (auto const& [p1, p2, p3] : o2::soa::combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(SliceParticle, SliceParticle, SliceParticle))) { + + // check if triplet is clean + if (!TcManager.isCleanTriplet(p1, p2, p3, TrackTable)) { + continue; + } + + // check if triplet is close + CtrManager.setTriplet(p1, p2, p3, TrackTable); + if (CtrManager.isCloseTriplet()) { + continue; + } + + // Randomize pair order if enabled + switch (tripletOrder) { + case kOrder123: + TripletHistManager.setTriplet(p1, p2, p3, Collision); + break; + case kOrder213: + TripletHistManager.setTriplet(p2, p1, p3, Collision); + break; + case kOrder132: + TripletHistManager.setTriplet(p1, p3, p2, Collision); + break; + case kOrder321: + TripletHistManager.setTriplet(p3, p2, p1, Collision); + break; + default: + TripletHistManager.setTriplet(p1, p2, p3, Collision); + break; + } + + // fill deta-dphi histograms with !3 cutoff + CtrManager.fill(TripletHistManager.getQ3()); + + // if triplet cuts are configured check them before filling + if (TripletHistManager.checkTripletCuts()) { + TripletHistManager.template fill(); + } + } +} + +// process same event for identical 2 particles and 1 other particle +template +void processSameEvent(T1 const& SliceParticle1, // 1&2 have same species + T2 const& SliceParticle3, + T3 const& TrackTable, + T4 const& Collision, + T5& ParticleHistManager1, + T6& ParticleHistManager3, + T7& TripletHistManager, + T8& CtrManager, + T9& TcManager, + TripletOrder triplerOrder) +{ + for (auto const& part : SliceParticle1) { + ParticleHistManager1.template fill(part, TrackTable); + } + for (auto const& part : SliceParticle3) { + ParticleHistManager3.template fill(part, TrackTable); + } + + for (auto const& [p1, p2, p3] : o2::soa::combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(SliceParticle1, SliceParticle1, SliceParticle3))) { + + // check if triplet is clean + if (!TcManager.isCleanTriplet(p1, p2, p3, TrackTable)) { + continue; + } + + // check if triplet is close + CtrManager.setTriplet(p1, p2, p3, TrackTable); + if (CtrManager.isCloseTriplet()) { + continue; + } + + // Randomize triplet order if enabled + switch (triplerOrder) { + case kOrder123: + TripletHistManager.setTriplet(p1, p2, p3, Collision); + break; + case kOrder213: + TripletHistManager.setTriplet(p1, p2, p3, Collision); + break; + default: + TripletHistManager.setTriplet(p1, p2, p3, Collision); + break; + } + + // fill deta-dphi histograms with !3 cutoff + CtrManager.fill(TripletHistManager.getQ3()); + + // if triplet cuts are configured check them before filling + if (TripletHistManager.checkTripletCuts()) { + TripletHistManager.template fill(); + } + } +} + +// process same event for 3 different particles +template +void processSameEvent(T1 const& SliceParticle1, + T2 const& SliceParticle2, + T3 const& SliceParticle3, + T4 const& TrackTable, + T5 const& Collision, + T6& ParticleHistManager1, + T7& ParticleHistManager2, + T8& ParticleHistManager3, + T9& TripletHistManager, + T10& CtrManager, + T11& TcManager) +{ + for (auto const& part : SliceParticle1) { + ParticleHistManager1.template fill(part, TrackTable); + } + for (auto const& part : SliceParticle2) { + ParticleHistManager2.template fill(part, TrackTable); + } + for (auto const& part : SliceParticle3) { + ParticleHistManager3.template fill(part, TrackTable); + } + + for (auto const& [p1, p2, p3] : o2::soa::combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(SliceParticle1, SliceParticle2, SliceParticle3))) { + + // check if triplet is clean + if (!TcManager.isCleanTriplet(p1, p2, p3, TrackTable)) { + continue; + } + + // check if triplet is close + CtrManager.setTriplet(p1, p2, p3, TrackTable); + if (CtrManager.isCloseTriplet()) { + continue; + } + + TripletHistManager.setTriplet(p1, p2, p3, Collision); + + // fill deta-dphi histograms with !3 cutoff + CtrManager.fill(TripletHistManager.getQ3()); + + // if triplet cuts are configured check them before filling + if (TripletHistManager.checkTripletCuts()) { + TripletHistManager.template fill(); + } + } +} + +// // process same event for 3 identical particles with mc information +template +void processSameEvent(T1 const& SliceParticle, + T2 const& TrackTable, + T3 const& mcParticles, + T4 const& mcMothers, + T5 const& mcPartonicMothers, + T6 const& Collision, + T7 const& mcCollisions, + T8& ParticleHistManager, + T9& TripletHistManager, + T10& CtrManager, + T11& TcManager, + int swapTriplet) +{ + for (auto const& part : SliceParticle) { + ParticleHistManager.template fill(part, TrackTable, mcParticles, mcMothers, mcPartonicMothers); + } + for (auto const& [p1, p2, p3] : o2::soa::combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(SliceParticle, SliceParticle, SliceParticle))) { + // check if pair is clean + if (!TcManager.isCleanTriplet(p1, p2, p3, TrackTable, mcPartonicMothers)) { + continue; + } + // check if pair is close + CtrManager.setTriplet(p1, p2, p3, TrackTable); + if (CtrManager.isCloseTriplet()) { + continue; + } + // Randomize pair order if enabled + // + switch (swapTriplet) { + case 3: + TripletHistManager.setTripletMc(p1, p3, p2, mcParticles, Collision, mcCollisions); + break; + case 2: + TripletHistManager.setTripletMc(p2, p1, p3, mcParticles, Collision, mcCollisions); + break; + case 1: + TripletHistManager.setTripletMc(p1, p2, p3, mcParticles, Collision, mcCollisions); + break; + default: + TripletHistManager.setTripletMc(p1, p2, p3, mcParticles, Collision, mcCollisions); + break; + } + // float threshold = 0.5f; + // fill deta-dphi histograms with q3 cutoff + CtrManager.fill(TripletHistManager.getQ3()); + // if pair cuts are configured check them before filling + if (TripletHistManager.checkTripletCuts()) { + TripletHistManager.template fill(); + } + } +} + +// process same event for 2 identical particles and one other with mc information +template +void processSameEvent(T1 const& SliceParticle1, + T2 const& SliceParticle3, + T3 const& TrackTable, + T4 const& mcParticles, + T5 const& mcMothers, + T6 const& mcPartonicMothers, + T7 const& Collision, + T8 const& mcCollisions, + T9& ParticleHistManager1, + T10& ParticleHistManager3, + T11& TripletHistManager, + T12& CtrManager, + T13& TcManager, + int swapTriplet) +{ + for (auto const& part : SliceParticle1) { + ParticleHistManager1.template fill(part, TrackTable, mcParticles, mcMothers, mcPartonicMothers); + } + for (auto const& part : SliceParticle3) { + ParticleHistManager3.template fill(part, TrackTable, mcParticles, mcMothers, mcPartonicMothers); + } + for (auto const& [p1, p2, p3] : o2::soa::combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(SliceParticle1, SliceParticle1, SliceParticle3))) { + // check if pair is clean + if (!TcManager.isCleanTriplet(p1, p2, p3, TrackTable, mcPartonicMothers)) { + continue; + } + // check if pair is close + CtrManager.setTriplet(p1, p2, p3, TrackTable); + if (CtrManager.isCloseTriplet()) { + continue; + } + // Randomize pair order if enabled + // + switch (swapTriplet) { + case 2: + TripletHistManager.setTripletMc(p2, p1, p3, mcParticles, Collision, mcCollisions); + break; + case 1: + TripletHistManager.setTripletMc(p1, p2, p3, mcParticles, Collision, mcCollisions); + break; + default: + TripletHistManager.setTripletMc(p1, p2, p3, mcParticles, Collision, mcCollisions); + break; + } + // fill deta-dphi histograms with q3 cutoff + CtrManager.fill(TripletHistManager.getQ3()); + // if pair cuts are configured check them before filling + if (TripletHistManager.checkTripletCuts()) { + TripletHistManager.template fill(); + } + } +} + +// process same event for 3 different particles +template +void processSameEvent(T1 const& SliceParticle1, + T2 const& SliceParticle2, + T3 const& SliceParticle3, + T4 const& TrackTable, + T5 const& mcParticles, + T6 const& mcMothers, + T7 const& mcPartonicMothers, + T8 const& Collision, + T9 const& mcCollisions, + T10& ParticleHistManager1, + T11& ParticleHistManager2, + T12& ParticleHistManager3, + T13& TripletHistManager, + T14& CtrManager, + T15& TcManager) +{ + for (auto const& part : SliceParticle1) { + ParticleHistManager1.template fill(part, TrackTable, mcParticles, mcMothers, mcPartonicMothers); + } + for (auto const& part : SliceParticle2) { + ParticleHistManager2.template fill(part, TrackTable, mcParticles, mcMothers, mcPartonicMothers); + } + for (auto const& part : SliceParticle3) { + ParticleHistManager3.template fill(part, TrackTable, mcParticles, mcMothers, mcPartonicMothers); + } + for (auto const& [p1, p2, p3] : o2::soa::combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(SliceParticle1, SliceParticle2, SliceParticle3))) { + // check if pair is clean + if (!TcManager.isCleanTriplet(p1, p2, p3, TrackTable, mcPartonicMothers)) { + continue; + } + // check if pair is close + CtrManager.setTriplet(p1, p2, p3, TrackTable); + if (CtrManager.isCloseTriplet()) { + continue; + } + TripletHistManager.setTripletMc(p1, p2, p3, mcParticles, Collision, mcCollisions); + // fill deta-dphi histograms with q3 cutoff + CtrManager.fill(TripletHistManager.getQ3()); + // if pair cuts are configured check them before filling + if (TripletHistManager.checkTripletCuts()) { + TripletHistManager.template fill(); + } + } +} + +// process mixed event +template +void processMixedEvent(T1 const& Collisions, + T2& Partition1, + T3& Partition2, + T4& Partition3, + T5 const& TrackTable, + T6& cache, + T7 const& policy, + T8 const& depth, + T9& TripletHistManager, + T10& CtrManager, + T11& TcManager) +{ + for (auto const& [collision1, collision2, collision3] : o2::soa::selfCombinations(policy, depth, -1, Collisions, Collisions, Collisions)) { + if (collision1.magField() != collision2.magField() || + collision2.magField() != collision3.magField() || + collision1.magField() != collision3.magField()) { + continue; + } + CtrManager.setMagField(collision1.magField()); + auto sliceParticle1 = Partition1->sliceByCached(o2::aod::femtobase::stored::fColId, collision1.globalIndex(), cache); + auto sliceParticle2 = Partition2->sliceByCached(o2::aod::femtobase::stored::fColId, collision2.globalIndex(), cache); + auto sliceParticle3 = Partition3->sliceByCached(o2::aod::femtobase::stored::fColId, collision3.globalIndex(), cache); + if (sliceParticle1.size() == 0 || sliceParticle2.size() == 0 || sliceParticle3.size() == 0) { + continue; + } + for (auto const& [p1, p2, p3] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(sliceParticle1, sliceParticle2, sliceParticle3))) { + // pair cleaning + if (!TcManager.isCleanTriplet(p1, p2, p3, TrackTable)) { + continue; + } + // Close pair rejection + CtrManager.setTriplet(p1, p2, p3, TrackTable); + if (CtrManager.isCloseTriplet()) { + continue; + } + TripletHistManager.setTriplet(p1, p2, p3, collision1, collision2, collision3); + CtrManager.fill(TripletHistManager.getQ3()); + if (TripletHistManager.checkTripletCuts()) { + TripletHistManager.template fill(); + } + } + } +} + +// process mixed event in mc +template +void processMixedEvent(T1 const& Collisions, + T2 const& mcCollisions, + T3& Partition1, + T4& Partition2, + T5& Partition3, + T6 const& TrackTable, + T7 const& mcParticles, + T8& cache, + T9 const& policy, + T10 const& depth, + T11& TripletHistManager, + T12& CtrManager, + T13& TcManager) +{ + for (auto const& [collision1, collision2, collision3] : o2::soa::selfCombinations(policy, depth, -1, Collisions, Collisions, Collisions)) { + if (collision1.magField() != collision2.magField() || + collision2.magField() != collision3.magField() || + collision1.magField() != collision3.magField()) { + continue; + } + CtrManager.setMagField(collision1.magField()); + auto sliceParticle1 = Partition1->sliceByCached(o2::aod::femtobase::stored::fColId, collision1.globalIndex(), cache); + auto sliceParticle2 = Partition2->sliceByCached(o2::aod::femtobase::stored::fColId, collision2.globalIndex(), cache); + auto sliceParticle3 = Partition3->sliceByCached(o2::aod::femtobase::stored::fColId, collision3.globalIndex(), cache); + if (sliceParticle1.size() == 0 || sliceParticle2.size() == 0 || sliceParticle3.size() == 0) { + continue; + } + for (auto const& [p1, p2, p3] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(sliceParticle1, sliceParticle2, sliceParticle3))) { + // pair cleaning + if (!TcManager.isCleanTriplet(p1, p2, p3, TrackTable)) { + continue; + } + // Close pair rejection + CtrManager.setTriplet(p1, p2, p3, TrackTable); + if (CtrManager.isCloseTriplet()) { + continue; + } + TripletHistManager.setTripletMc(p1, p2, p3, mcParticles, collision1, collision2, collision3, mcCollisions); + CtrManager.fill(TripletHistManager.getQ3()); + if (TripletHistManager.checkTripletCuts()) { + TripletHistManager.template fill(); + } + } + } +} + +} // namespace tripletprocesshelpers +} // namespace o2::analysis::femto + +#endif // PWGCF_FEMTO_CORE_TRIPLETPROCESSHELPERS_H_ diff --git a/PWGCF/Femto/TableProducer/femtoProducer.cxx b/PWGCF/Femto/TableProducer/femtoProducer.cxx index 36897c02afb..a86ce731cb9 100644 --- a/PWGCF/Femto/TableProducer/femtoProducer.cxx +++ b/PWGCF/Femto/TableProducer/femtoProducer.cxx @@ -19,7 +19,6 @@ #include "PWGCF/Femto/Core/mcBuilder.h" #include "PWGCF/Femto/Core/modes.h" #include "PWGCF/Femto/Core/trackBuilder.h" -#include "PWGCF/Femto/Core/twoTrackResonanceBuilder.h" #include "PWGCF/Femto/Core/v0Builder.h" #include "PWGLF/DataModel/LFKinkDecayTables.h" #include "PWGLF/DataModel/LFStrangenessTables.h" @@ -38,7 +37,6 @@ #include "Framework/AnalysisHelpers.h" #include "Framework/AnalysisTask.h" #include "Framework/Configurable.h" -#include "Framework/Expressions.h" #include "Framework/HistogramRegistry.h" #include "Framework/InitContext.h" #include "Framework/OutputObjHeader.h" @@ -132,36 +130,6 @@ struct FemtoProducer { kinkbuilder::ConfSigmaPlusBits confSigmaPlusBits; kinkbuilder::KinkBuilder sigmaPlusBuilder; - // resonance daughter filters and partitions - twotrackresonancebuilder::ConfTwoTrackResonanceDaughterFilters confResonanceDaughterFilters; - // caching and preslicing - SliceCache cache; - Preslice perColTracks = track::collisionId; - Partition partitionPositiveDaughters = - (track::signed1Pt > 0.f) && - (track::pt > confResonanceDaughterFilters.ptMin && track::pt < confResonanceDaughterFilters.ptMax) && - (track::eta > confResonanceDaughterFilters.etaMin && track::eta < confResonanceDaughterFilters.etaMax) && - (track::phi > confResonanceDaughterFilters.phiMin && track::phi < confResonanceDaughterFilters.phiMax); - Partition partitionNegativeDaughters = - (track::signed1Pt < 0.f) && - (track::pt > confResonanceDaughterFilters.ptMin && track::pt < confResonanceDaughterFilters.ptMax) && - (track::eta > confResonanceDaughterFilters.etaMin && track::eta < confResonanceDaughterFilters.etaMax) && - (track::phi > confResonanceDaughterFilters.phiMin && track::phi < confResonanceDaughterFilters.phiMax); - - // resonance builders - twotrackresonancebuilder::TwoTrackResonanceBuilderProducts twoTrackResonanceBuilderProducts; - twotrackresonancebuilder::ConfTwoTrackResonanceTables confTwoTrackResonanceTables; - twotrackresonancebuilder::ConfRhoFilters confRhoFilters; - twotrackresonancebuilder::ConfRho0Bits confRho0Bits; - twotrackresonancebuilder::TwoTrackResonanceBuilder rho0Builder; - twotrackresonancebuilder::ConfPhiFilters confPhiFilters; - twotrackresonancebuilder::ConfPhiBits confPhiBits; - twotrackresonancebuilder::TwoTrackResonanceBuilder phiBuilder; - twotrackresonancebuilder::ConfKstarFilters confKstarFilters; - twotrackresonancebuilder::ConfKstar0Bits confKstar0Bits; - twotrackresonancebuilder::TwoTrackResonanceBuilder kstar0Builder; - twotrackresonancebuilder::TwoTrackResonanceBuilder kstar0barBuilder; - // mc builder mcbuilder::ConfMc confMc; mcbuilder::ConfMcTables confMcTables; @@ -216,12 +184,6 @@ struct FemtoProducer { xiBuilder.init(&hRegistry, confXiBits, confCascadeFilters, confCascadeTables, context); omegaBuilder.init(&hRegistry, confOmegaBits, confCascadeFilters, confCascadeTables, context); - // configure resonance selections - rho0Builder.init(&hRegistry, confRho0Bits, confRhoFilters, confResonanceDaughterFilters, confTwoTrackResonanceTables, context); - phiBuilder.init(&hRegistry, confPhiBits, confPhiFilters, confResonanceDaughterFilters, confTwoTrackResonanceTables, context); - kstar0Builder.init(&hRegistry, confKstar0Bits, confKstarFilters, confResonanceDaughterFilters, confTwoTrackResonanceTables, context); - kstar0barBuilder.init(&hRegistry, confKstar0Bits, confKstarFilters, confResonanceDaughterFilters, confTwoTrackResonanceTables, context); - // configure mcBuilder mcBuilder.init(confMc, confMcTables, context); @@ -243,11 +205,11 @@ struct FemtoProducer { return true; } - template - bool processMcCollisions(T1 const& col, T2 const& mcCols, T3 const& /* bcs*/, T4 const& tracks) + template + bool processMcCollisions(T1 const& col, T2 const& mcCols, T3 const& /* bcs*/, T4 const& tracks, T5 const& mcParticles) { - mcBuilder.reset(); collisionBuilder.reset(); + mcBuilder.reset(mcCols, mcParticles); // we call this function always first so we reset the mcBuilder here auto bc = col.template bc_as(); collisionBuilder.initCollision(bc, col, tracks, ccdb, hRegistry); if (!collisionBuilder.checkCollision(col, mcCols)) { @@ -261,26 +223,17 @@ struct FemtoProducer { template void processTracks(T1 const& col, T2 const& tracksWithItsPid) { + trackBuilder.reset(tracksWithItsPid); trackBuilder.fillTracks(col, collisionBuilder, collisionBuilderProducts, tracksWithItsPid, trackBuilderProducts); } template void processMcTracks(T1 const& col, T2 const& mcCols, T3 const& tracks, T4 const& tracksWithItsPid, T5 const& mcParticles) { + trackBuilder.reset(tracksWithItsPid); trackBuilder.fillMcTracks(col, collisionBuilder, collisionBuilderProducts, mcCols, tracks, tracksWithItsPid, trackBuilderProducts, mcParticles, mcBuilder, mcProducts); } - template - void processResonances(T1 const& col, T2 const& tracks, T3 const& tracksWithItsPid) - { - auto groupPositiveTracks = partitionPositiveDaughters->sliceByCached(track::collisionId, col.globalIndex(), cache); - auto groupNegativeTracks = partitionNegativeDaughters->sliceByCached(track::collisionId, col.globalIndex(), cache); - rho0Builder.fillResonances(col, collisionBuilder, collisionBuilderProducts, trackBuilderProducts, twoTrackResonanceBuilderProducts, groupPositiveTracks, groupNegativeTracks, tracks, tracksWithItsPid, trackBuilder); - phiBuilder.fillResonances(col, collisionBuilder, collisionBuilderProducts, trackBuilderProducts, twoTrackResonanceBuilderProducts, groupPositiveTracks, groupNegativeTracks, tracks, tracksWithItsPid, trackBuilder); - kstar0Builder.fillResonances(col, collisionBuilder, collisionBuilderProducts, trackBuilderProducts, twoTrackResonanceBuilderProducts, groupPositiveTracks, groupNegativeTracks, tracks, tracksWithItsPid, trackBuilder); - kstar0barBuilder.fillResonances(col, collisionBuilder, collisionBuilderProducts, trackBuilderProducts, twoTrackResonanceBuilderProducts, groupPositiveTracks, groupNegativeTracks, tracks, tracksWithItsPid, trackBuilder); - } - // add v0s template void processV0s(T1 const& col, T2 const& tracks, T3 const& tracksWithItsPid, T4 const& v0s) @@ -334,7 +287,6 @@ struct FemtoProducer { auto tracksWithItsPid = o2::soa::Attach(tracks); processTracks(col, tracksWithItsPid); - processResonances(col, tracks, tracksWithItsPid); } PROCESS_SWITCH(FemtoProducer, processTracksRun3pp, "Process tracks", true); @@ -350,7 +302,6 @@ struct FemtoProducer { auto tracksWithItsPid = o2::soa::Attach(tracks); processTracks(col, tracksWithItsPid); - processResonances(col, tracks, tracksWithItsPid); processV0s(col, tracks, tracksWithItsPid, v0s); }; PROCESS_SWITCH(FemtoProducer, processTracksV0sRun3pp, "Process tracks and v0s", false); @@ -367,7 +318,6 @@ struct FemtoProducer { auto tracksWithItsPid = o2::soa::Attach(tracks); processTracks(col, tracksWithItsPid); - processResonances(col, tracks, tracksWithItsPid); processKinks(col, tracks, tracksWithItsPid, kinks); } PROCESS_SWITCH(FemtoProducer, processTracksKinksRun3pp, "Process tracks and kinks", false); @@ -385,7 +335,6 @@ struct FemtoProducer { auto tracksWithItsPid = o2::soa::Attach(tracks); processTracks(col, tracksWithItsPid); - processResonances(col, tracks, tracksWithItsPid); processV0s(col, tracks, tracksWithItsPid, v0s); processCascades(col, tracks, tracksWithItsPid, cascades); } @@ -405,7 +354,6 @@ struct FemtoProducer { auto tracksWithItsPid = o2::soa::Attach(tracks); processTracks(col, tracksWithItsPid); - processResonances(col, tracks, tracksWithItsPid); processV0s(col, tracks, tracksWithItsPid, v0s); processKinks(col, tracks, tracksWithItsPid, kinks); processCascades(col, tracks, tracksWithItsPid, cascades); @@ -419,7 +367,7 @@ struct FemtoProducer { consumeddata::Run3McRecoTracks const& tracks, consumeddata::Run3McGenParticles const& mcParticles) { - if (!processMcCollisions(col, mcCols, bcs, tracks)) { + if (!processMcCollisions(col, mcCols, bcs, tracks, mcParticles)) { return; } auto tracksWithItsPid = o2::soa::Attach(col, mcCols, bcs, tracks)) { + if (!processMcCollisions(col, mcCols, bcs, tracks, mcParticles)) { return; } auto tracksWithItsPid = o2::soa::Attach(col, mcCols, bcs, tracks)) { + if (!processMcCollisions(col, mcCols, bcs, tracks, mcParticles)) { return; } auto tracksWithItsPid = o2::soa::Attach(col, mcCols, bcs, tracks)) { + if (!processMcCollisions(col, mcCols, bcs, tracks, mcParticles)) { return; } auto tracksWithItsPid = o2::soa::Attach +#include + +using namespace o2; +using namespace o2::aod; +using namespace o2::soa; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::analysis::femto; + +struct FemtoTripletTrackTrackTrack { + + // setup tables + using FemtoCollisions = Join; + using FilteredFemtoCollisions = o2::soa::Filtered; + using FilteredFemtoCollision = FilteredFemtoCollisions::iterator; + + using FemtoCollisionsWithLabel = o2::soa::Join; + using FilteredFemtoCollisionsWithLabel = o2::soa::Filtered; + using FilteredFemtoCollisionWithLabel = FilteredFemtoCollisionsWithLabel::iterator; + + using FemtoTracks = o2::soa::Join; + + using FemtoTracksWithLabel = o2::soa::Join; + + SliceCache cache; + + // setup collisions + collisionbuilder::ConfCollisionSelection collisionSelection; + Filter collisionFilter = MAKE_COLLISION_FILTER(collisionSelection); + colhistmanager::ConfCollisionBinning confCollisionBinning; + + // setup tracks + trackbuilder::ConfTrackSelection1 confTrackSelections1; + trackhistmanager::ConfTrackBinning1 confTrackBinning1; + trackbuilder::ConfTrackSelection2 confTrackSelections2; + trackhistmanager::ConfTrackBinning2 confTrackBinning2; + trackbuilder::ConfTrackSelection3 confTrackSelections3; + trackhistmanager::ConfTrackBinning3 confTrackBinning3; + + Partition trackPartition1 = MAKE_TRACK_PARTITION(confTrackSelections1); + Partition trackPartition2 = MAKE_TRACK_PARTITION(confTrackSelections2); + Partition trackPartition3 = MAKE_TRACK_PARTITION(confTrackSelections3); + Preslice perColtracks = aod::femtobase::stored::fColId; + + Partition trackWithLabelPartition1 = MAKE_TRACK_PARTITION(confTrackSelections1); + Partition trackWithLabelPartition2 = MAKE_TRACK_PARTITION(confTrackSelections2); + Partition trackWithLabelPartition3 = MAKE_TRACK_PARTITION(confTrackSelections3); + Preslice perColtracksWithLabel = aod::femtobase::stored::fColId; + + // setup triplets + triplethistmanager::ConfTripletBinning confTripletBinning; + triplethistmanager::ConfTripletCuts confTripletCuts; + + closetripletrejection::ConfCtrTrackTrackTrack confCtr; + + tripletbuilder::TripletTrackTrackTrackBuilder< + trackhistmanager::PrefixTrack1, + trackhistmanager::PrefixTrack2, + trackhistmanager::PrefixTrack3, + triplethistmanager::PrefixTrackTrackTrackSe, + triplethistmanager::PrefixTrackTrackTrackMe, + closetripletrejection::PrefixTrack1Track2Se, + closetripletrejection::PrefixTrack2Track3Se, + closetripletrejection::PrefixTrack1Track3Se, + closetripletrejection::PrefixTrack1Track2Me, + closetripletrejection::PrefixTrack2Track3Me, + closetripletrejection::PrefixTrack1Track3Me> + tripletTrackTrackTrackBuilder; + + // setup mixing + std::vector defaultVtxBins{10, -10, 10}; + std::vector defaultMultBins{50, 0, 200}; + std::vector defaultCentBins{10, 0, 100}; + ColumnBinningPolicy mixBinsVtxMult{{defaultVtxBins, defaultMultBins}, true}; + ColumnBinningPolicy mixBinsVtxCent{{defaultVtxBins, defaultCentBins}, true}; + ColumnBinningPolicy mixBinsVtxMultCent{{defaultVtxBins, defaultMultBins, defaultCentBins}, true}; + triplethistmanager::ConfMixing confMixing; + + HistogramRegistry hRegistry{"FemtoTrackTrackTrack", {}, OutputObjHandlingPolicy::AnalysisObject}; + + void init(InitContext&) + { + if ((doprocessSameEvent + doprocessSameEventMc) > 1 || (doprocessMixedEvent + doprocessMixedEventMc) > 1) { + LOG(fatal) << "More than 1 same or mixed event process function is activated. Breaking..."; + } + bool processData = doprocessSameEvent || doprocessMixedEvent; + bool processMc = doprocessSameEventMc || doprocessMixedEventMc; + if (processData && processMc) { + LOG(fatal) << "Both data and mc processing is activated. Breaking..."; + } + + // setup columnpolicy for binning + // default values are used during instantiation, so we need to explicity update them here + mixBinsVtxMult = {{confMixing.vtxBins, confMixing.multBins.value}, true}; + mixBinsVtxCent = {{confMixing.vtxBins.value, confMixing.centBins.value}, true}; + mixBinsVtxMultCent = {{confMixing.vtxBins.value, confMixing.multBins.value, confMixing.centBins.value}, true}; + + // setup histogram specs + auto ctrHistSpec = closepairrejection::makeCprHistSpecMap(confCtr); + + if (processData) { + auto colHistSpec = colhistmanager::makeColHistSpecMap(confCollisionBinning); + auto trackHistSpec1 = trackhistmanager::makeTrackHistSpecMap(confTrackBinning1); + auto trackHistSpec2 = trackhistmanager::makeTrackHistSpecMap(confTrackBinning2); + auto trackHistSpec3 = trackhistmanager::makeTrackHistSpecMap(confTrackBinning3); + auto tripletHistSpec = triplethistmanager::makeTripletHistSpecMap(confTripletBinning); + tripletTrackTrackTrackBuilder.init(&hRegistry, confTrackSelections1, confTrackSelections2, confTrackSelections3, confCtr, confMixing, confTripletBinning, confTripletCuts, colHistSpec, trackHistSpec1, trackHistSpec2, trackHistSpec3, tripletHistSpec, ctrHistSpec); + } else { + auto colHistSpec = colhistmanager::makeColMcHistSpecMap(confCollisionBinning); + auto trackHistSpec1 = trackhistmanager::makeTrackMcHistSpecMap(confTrackBinning1); + auto trackHistSpec2 = trackhistmanager::makeTrackMcHistSpecMap(confTrackBinning2); + auto trackHistSpec3 = trackhistmanager::makeTrackMcHistSpecMap(confTrackBinning3); + auto tripletHistSpec = triplethistmanager::makeTripletMcHistSpecMap(confTripletBinning); + tripletTrackTrackTrackBuilder.init(&hRegistry, confTrackSelections1, confTrackSelections2, confTrackSelections3, confCtr, confMixing, confTripletBinning, confTripletCuts, colHistSpec, trackHistSpec1, trackHistSpec2, trackHistSpec3, tripletHistSpec, ctrHistSpec); + } + hRegistry.print(); + }; + + void processSameEvent(FilteredFemtoCollision const& col, FemtoTracks const& tracks) + { + tripletTrackTrackTrackBuilder.processSameEvent(col, tracks, trackPartition1, trackPartition2, trackPartition3, cache); + } + PROCESS_SWITCH(FemtoTripletTrackTrackTrack, processSameEvent, "Enable processing same event processing", true); + + void processSameEventMc(FilteredFemtoCollisionWithLabel const& col, FMcCols const& mcCols, FemtoTracksWithLabel const& tracks, FMcParticles const& mcParticles, FMcMothers const& mcMothers, FMcPartMoths const& mcPartonicMothers) + { + tripletTrackTrackTrackBuilder.processSameEvent(col, mcCols, tracks, trackWithLabelPartition1, trackWithLabelPartition2, trackWithLabelPartition3, mcParticles, mcMothers, mcPartonicMothers, cache); + } + PROCESS_SWITCH(FemtoTripletTrackTrackTrack, processSameEventMc, "Enable processing same event processing", false); + + void processMixedEvent(FilteredFemtoCollisions const& cols, FemtoTracks const& tracks) + { + tripletTrackTrackTrackBuilder.processMixedEvent(cols, tracks, trackPartition1, trackPartition2, trackPartition3, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent); + } + PROCESS_SWITCH(FemtoTripletTrackTrackTrack, processMixedEvent, "Enable processing mixed event processing", true); + + void processMixedEventMc(FilteredFemtoCollisionsWithLabel const& cols, FMcCols const& mcCols, FemtoTracksWithLabel const& tracks, FMcParticles const& mcParticles) + { + tripletTrackTrackTrackBuilder.processMixedEvent(cols, mcCols, tracks, trackWithLabelPartition1, trackWithLabelPartition2, trackWithLabelPartition3, mcParticles, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent); + } + PROCESS_SWITCH(FemtoTripletTrackTrackTrack, processMixedEventMc, "Enable processing mixed event processing", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{ + adaptAnalysisTask(cfgc), + }; + return workflow; +}