diff --git a/PWGEM/PhotonMeson/Core/CMakeLists.txt b/PWGEM/PhotonMeson/Core/CMakeLists.txt index 72a17aa1497..7e4eac119f4 100644 --- a/PWGEM/PhotonMeson/Core/CMakeLists.txt +++ b/PWGEM/PhotonMeson/Core/CMakeLists.txt @@ -18,6 +18,7 @@ o2physics_add_library(PWGEMPhotonMesonCore EMPhotonEventCut.cxx CutsLibrary.cxx HistogramsLibrary.cxx + EMBitFlags.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::MLCore O2Physics::PWGEMDileptonCore) o2physics_target_root_dictionary(PWGEMPhotonMesonCore @@ -29,4 +30,5 @@ o2physics_target_root_dictionary(PWGEMPhotonMesonCore EMPhotonEventCut.h CutsLibrary.h HistogramsLibrary.h + EMBitFlags.h LINKDEF PWGEMPhotonMesonCoreLinkDef.h) diff --git a/PWGEM/PhotonMeson/Core/EMBitFlags.cxx b/PWGEM/PhotonMeson/Core/EMBitFlags.cxx new file mode 100644 index 00000000000..d182eb03f40 --- /dev/null +++ b/PWGEM/PhotonMeson/Core/EMBitFlags.cxx @@ -0,0 +1,72 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file EMBitFlags.cxx +/// \brief Source of bit flag class for particle selection. +/// \author M. Hemmer, marvin.hemmer@cern.ch + +#include "EMBitFlags.h" + +#include // fill +#include // assert +#include // size_t + +EMBitFlags::EMBitFlags(std::size_t nBits) + : mBits((nBits + 63) / 64, ~0ULL), + mSize(nBits) +{ +} + +std::size_t EMBitFlags::size() const +{ + return mSize; +} + +bool EMBitFlags::test(std::size_t i) const +{ + assert(i < mSize); + return mBits[word(i)] & mask(i); +} + +void EMBitFlags::set(std::size_t i) +{ + ensureSize(i + 1); + mBits[word(i)] &= ~mask(i); +} + +void EMBitFlags::reset(std::size_t i) +{ + assert(i < mSize); + mBits[word(i)] |= mask(i); +} + +void EMBitFlags::clear() +{ + std::fill(mBits.begin(), mBits.end(), ~0ULL); +} + +void EMBitFlags::reserve(std::size_t nBits) +{ + mBits.reserve((nBits + 63) / 64); +} + +void EMBitFlags::resize(std::size_t nBits) +{ + mBits.resize((nBits + 63) / 64, ~0ULL); + mSize = nBits; +} + +void EMBitFlags::ensureSize(std::size_t nBits) +{ + if (nBits > mSize) { + resize(nBits); + } +} diff --git a/PWGEM/PhotonMeson/Core/EMBitFlags.h b/PWGEM/PhotonMeson/Core/EMBitFlags.h new file mode 100644 index 00000000000..d8a47390cc1 --- /dev/null +++ b/PWGEM/PhotonMeson/Core/EMBitFlags.h @@ -0,0 +1,73 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file EMBitFlags.h +/// \brief Header of bit flag class for particle selection. +/// \author M. Hemmer, marvin.hemmer@cern.ch + +#ifndef PWGEM_PHOTONMESON_CORE_EMBITFLAGS_H_ +#define PWGEM_PHOTONMESON_CORE_EMBITFLAGS_H_ + +#include // size_t +#include // uint64_t +#include // vector + +/// \class EMBitFlags +/// \brief Dynamically-sized bit container with bit-level storage. +/// +/// Bits can be set beyond the current size, in which case the container +/// grows automatically. Access via test() and reset() requires the index +/// to be within the current size. Bits are all on by default and will be +/// switched off if particle fails a cut +class EMBitFlags +{ + public: + explicit EMBitFlags(std::size_t nBits = 0); + + /// \brief get number of stored bits + std::size_t size() const; + + /// \brief check bit i + /// \param i index of bit that should be checked + bool test(std::size_t i) const; + + /// \brief set bit i + /// \param i index of bit which value should be set + void set(std::size_t i); + + /// \brief reset bit i + /// \param i index of bit which value should be reset + void reset(std::size_t i); + + /// \brief resetting all flags to false + void clear(); + + /// \brief reserve space in the underlying storage for nBits bits + /// \param nBits number of bits to reserve capacity for + void reserve(std::size_t nBits); + + /// \brief resize the container to hold nBits bits + /// \param nBits new number of bits (new bits are initialized to false) + void resize(std::size_t nBits); + + private: + /// \brief ensure that the container can hold at least nBits bits + /// \param nBits required number of bits + void ensureSize(std::size_t nBits); + + static constexpr std::size_t word(std::size_t i) { return i >> 6; } + static constexpr std::uint64_t mask(std::size_t i) { return 1ULL << (i & 63); } + + std::vector mBits; + std::size_t mSize = 0; +}; + +#endif // PWGEM_PHOTONMESON_CORE_EMBITFLAGS_H_ diff --git a/PWGEM/PhotonMeson/Core/EMCPhotonCut.h b/PWGEM/PhotonMeson/Core/EMCPhotonCut.h index 7847253a4f0..3105ece21cd 100644 --- a/PWGEM/PhotonMeson/Core/EMCPhotonCut.h +++ b/PWGEM/PhotonMeson/Core/EMCPhotonCut.h @@ -16,6 +16,7 @@ #ifndef PWGEM_PHOTONMESON_CORE_EMCPHOTONCUT_H_ #define PWGEM_PHOTONMESON_CORE_EMCPHOTONCUT_H_ +#include "PWGEM/PhotonMeson/Core/EMBitFlags.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" #include @@ -31,16 +32,25 @@ #include template -concept isTrackContainer = o2::soa::is_table && requires(T t) { +concept IsTrackIterator = o2::soa::is_iterator && requires(T t) { // Check that the *elements* of the container have the required methods: - { t.begin.deltaEta() } -> std::same_as; - { t.begin.deltaPhi() } -> std::same_as; - { t.begin.trackPt() } -> std::same_as; - { t.begin.trackP() } -> std::same_as; + { t.deltaEta() } -> std::same_as; + { t.deltaPhi() } -> std::same_as; + { t.trackPt() } -> std::same_as; + { t.trackP() } -> std::same_as; +}; + +template +concept IsTrackContainer = o2::soa::is_table && requires(T t) { + // Check that the *elements* of the container have the required methods: + { t.begin().deltaEta() } -> std::same_as; + { t.begin().deltaPhi() } -> std::same_as; + { t.begin().trackPt() } -> std::same_as; + { t.begin().trackP() } -> std::same_as; }; template -concept hasTrackMatching = requires(Cluster cluster) { +concept HasTrackMatching = requires(Cluster cluster) { // requires that the following are valid calls: { cluster.deltaEta() } -> std::convertible_to>; { cluster.deltaPhi() } -> std::convertible_to>; @@ -49,7 +59,7 @@ concept hasTrackMatching = requires(Cluster cluster) { }; template -concept hasSecondaryMatching = requires(Cluster cluster) { +concept HasSecondaryMatching = requires(Cluster cluster) { // requires that the following are valid calls: { cluster.deltaEtaSec() } -> std::convertible_to>; { cluster.deltaPhiSec() } -> std::convertible_to>; @@ -57,16 +67,6 @@ concept hasSecondaryMatching = requires(Cluster cluster) { { cluster.trackpSec() } -> std::convertible_to>; }; -template -concept isSkimEMCClusterLike = requires(Cluster c) { - { c.e() } -> std::convertible_to; - { c.nCells() } -> std::convertible_to; - { c.m02() } -> std::convertible_to; - { c.time() } -> std::convertible_to; - { c.definition() } -> std::convertible_to; -} && (!hasTrackMatching) // important: skim clusters don’t have TM - &&(!hasSecondaryMatching); - struct TrackMatchingParams { float a{0.01f}; float b{4.07f}; @@ -94,6 +94,168 @@ class EMCPhotonCut : public TNamed static const char* mCutNames[static_cast(EMCPhotonCuts::kNCuts)]; + constexpr auto getClusterId(o2::soa::is_iterator auto const& t) const + { + if constexpr (requires { t.emEmcClusterId(); }) { + return t.emEmcClusterId(); + } else if constexpr (requires { t.minClusterId(); }) { + return t.minClusterId(); + } else { + return -1; + } + } + + /// \brief performs check if track is matched with given cluster + /// \param cluster cluster to be checked + /// \param emcmatchedtrack matched track iterator + /// \param emcmatchedtrackEnd matched track end iterator + /// \param GetEtaCut lambda to get the eta cut value + /// \param GetPhiCut lambda to get the phi cut value + /// \param applyEoverP bool to check if E/p should be checked (for secondaries we do not check this!) + bool checkTrackMatching(o2::soa::is_iterator auto const& cluster, IsTrackIterator auto& emcmatchedtrack, o2::soa::RowViewSentinel const emcmatchedtrackEnd, + bool applyEoverP, auto GetEtaCut, auto GetPhiCut) const + { + // advance to cluster + while (emcmatchedtrack != emcmatchedtrackEnd && getClusterId(emcmatchedtrack) < cluster.globalIndex()) { + ++emcmatchedtrack; + } + // all matched tracks have been checked + if (emcmatchedtrack == emcmatchedtrackEnd) { + return true; + } + // if all remaining tracks are beyond this cluster, it survives + if (getClusterId(emcmatchedtrack) > cluster.globalIndex()) { + return true; + } + // iterate over tracks belonging to this cluster + while (emcmatchedtrack != emcmatchedtrackEnd && getClusterId(emcmatchedtrack) == cluster.globalIndex()) { + auto dEta = std::fabs(emcmatchedtrack.deltaEta()); + auto dPhi = std::fabs(emcmatchedtrack.deltaPhi()); + auto trackpt = emcmatchedtrack.trackPt(); + auto trackp = emcmatchedtrack.trackP(); + bool fail = (dEta > GetEtaCut(trackpt)) || + (dPhi > GetPhiCut(trackpt)) || + (applyEoverP && cluster.e() / trackp >= mMinEoverP); + if (!fail) { + return false; // cluster got a track matche to it + } + ++emcmatchedtrack; + } + return true; // all tracks checked, cluster survives + } + + /// \brief check if given clusters survives all cuts + /// \param flags EMBitFlags where results will be stored + /// \param cluster cluster table to check + /// \param matchedTracks matched primary tracks table + /// \param matchedSecondaries matched secondary tracks table + void AreSelectedRunning(EMBitFlags& flags, o2::soa::is_table auto const& clusters, IsTrackContainer auto const& emcmatchedtracks, IsTrackContainer auto const& secondaries) const + { + auto emcmatchedtrackIter = emcmatchedtracks.begin(); + auto emcmatchedtrackEnd = emcmatchedtracks.end(); + auto secondaryIter = secondaries.begin(); + auto secondaryEnd = secondaries.end(); + size_t iCluster = 0; + for (const auto& cluster : clusters) { + if (!IsSelectedRunning(cluster, emcmatchedtrackIter, emcmatchedtrackEnd, secondaryIter, secondaryEnd)) { + flags.set(iCluster); + } + ++iCluster; + } + } + + /// \brief check if given cluster survives all cuts + /// \param cluster cluster to check + /// \param emcmatchedtrackIter current iterator of matched primary tracks + /// \param emcmatchedtrackEnd end iterator of matched primary tracks + /// \param secondaryIter current iterator of matched secondary tracks + /// \param secondaryEnd end iterator of matched secondary tracks + /// \return true if cluster survives all cuts else false + bool IsSelectedRunning(o2::soa::is_iterator auto const& cluster, IsTrackIterator auto& emcmatchedtrackIter, o2::soa::RowViewSentinel const emcmatchedtrackEnd, IsTrackIterator auto& secondaryIter, o2::soa::RowViewSentinel const secondaryEnd) const + { + if (!IsSelectedEMCalRunning(EMCPhotonCuts::kDefinition, cluster)) { + return false; + } + if (!IsSelectedEMCalRunning(EMCPhotonCuts::kEnergy, cluster)) { + return false; + } + if (!IsSelectedEMCalRunning(EMCPhotonCuts::kNCell, cluster)) { + return false; + } + if (!IsSelectedEMCalRunning(EMCPhotonCuts::kM02, cluster)) { + return false; + } + if (!IsSelectedEMCalRunning(EMCPhotonCuts::kTiming, cluster)) { + return false; + } + if (mUseTM && (!IsSelectedEMCalRunning(EMCPhotonCuts::kTM, cluster, emcmatchedtrackIter, emcmatchedtrackEnd))) { + return false; + } + if (mUseSecondaryTM && (!IsSelectedEMCalRunning(EMCPhotonCuts::kSecondaryTM, cluster, secondaryIter, secondaryEnd))) { + return false; + } + if (!IsSelectedEMCalRunning(EMCPhotonCuts::kExotic, cluster)) { + return false; + } + return true; + } + + /// \brief check if given cluster survives a given cut + /// \param cut enum of the cluster cut to check + /// \param cluster cluster to check + /// \return true if cluster survives cut else false + bool IsSelectedEMCalRunning(const EMCPhotonCuts& cut, o2::soa::is_iterator auto const& cluster) const + { + switch (cut) { + case EMCPhotonCuts::kDefinition: + return cluster.definition() == mDefinition; + + case EMCPhotonCuts::kEnergy: + return cluster.e() > mMinE; + + case EMCPhotonCuts::kNCell: + return cluster.nCells() >= mMinNCell; + + case EMCPhotonCuts::kM02: + return (cluster.nCells() == 1 || (mMinM02 <= cluster.m02() && cluster.m02() <= mMaxM02)); + + case EMCPhotonCuts::kTiming: + return mMinTime <= cluster.time() && cluster.time() <= mMaxTime; + + case EMCPhotonCuts::kTM: + return false; + + case EMCPhotonCuts::kSecondaryTM: + return false; + + case EMCPhotonCuts::kExotic: + return mUseExoticCut ? !cluster.isExotic() : true; + + default: + return true; + } + } + + /// \brief check if given cluster survives a given cut + /// \param cut enum of the cluster cut to check + /// \param cluster cluster to check + /// \param matchedTrackIter current iterator of matched primary or secondary tracks + /// \param matchedTrackEnd end iterator of matched primary or secondary tracks + /// \return true if cluster survives cut else false + bool IsSelectedEMCalRunning(const EMCPhotonCuts& cut, o2::soa::is_iterator auto const& cluster, IsTrackIterator auto& matchedTrackIter, o2::soa::RowViewSentinel const matchedTrackEnd) const + { + switch (cut) { + case EMCPhotonCuts::kTM: + return checkTrackMatching(cluster, matchedTrackIter, matchedTrackEnd, true, [this](float pt) { return GetTrackMatchingEta(pt); }, [this](float pt) { return GetTrackMatchingPhi(pt); }); + + case EMCPhotonCuts::kSecondaryTM: + return checkTrackMatching(cluster, matchedTrackIter, matchedTrackEnd, false, [this](float pt) { return GetSecTrackMatchingEta(pt); }, [this](float pt) { return GetSecTrackMatchingPhi(pt); }); + + default: + return true; + } + } + /// \brief check if given cluster survives all cuts /// \param cluster cluster to check /// \param matchedTracks subtable of the matched primary tracks (optional) @@ -155,7 +317,7 @@ class EMCPhotonCut : public TNamed return mMinTime <= cluster.time() && cluster.time() <= mMaxTime; case EMCPhotonCuts::kTM: { - if constexpr (isTrackContainer) { + if constexpr (IsTrackContainer) { for (const auto& emcmatchedtrack : emcmatchedtracks) { auto dEta = std::fabs(emcmatchedtrack.deltaEta()); auto dPhi = std::fabs(emcmatchedtrack.deltaPhi()); @@ -166,7 +328,7 @@ class EMCPhotonCut : public TNamed return false; } } - } else if constexpr (hasTrackMatching) { + } else if constexpr (HasTrackMatching) { auto dEtas = cluster.deltaEta(); // std:vector auto dPhis = cluster.deltaPhi(); // std:vector auto trackspt = cluster.trackpt(); // std:vector @@ -186,18 +348,18 @@ class EMCPhotonCut : public TNamed return true; // when we don't have any tracks the cluster should always survive the TM cut! } case EMCPhotonCuts::kSecondaryTM: { - if constexpr (isTrackContainer) { + if constexpr (IsTrackContainer) { for (const auto& emcmatchedtrack : emcmatchedtracks) { auto dEta = std::fabs(emcmatchedtrack.deltaEta()); auto dPhi = std::fabs(emcmatchedtrack.deltaPhi()); auto trackpt = emcmatchedtrack.trackPt(); auto trackp = emcmatchedtrack.trackP(); - bool result = (dEta > GetTrackMatchingEta(trackpt)) || (dPhi > GetTrackMatchingPhi(trackpt)) || (cluster.e() / trackp >= mMinEoverP); + bool result = (dEta > GetSecTrackMatchingEta(trackpt)) || (dPhi > GetSecTrackMatchingPhi(trackpt)) || (cluster.e() / trackp >= mMinEoverP); if (!result) { return false; } } - } else if constexpr (hasSecondaryMatching) { + } else if constexpr (HasSecondaryMatching) { auto dEtas = cluster.deltaEtaSec(); // std:vector auto dPhis = cluster.deltaPhiSec(); // std:vector auto trackspt = cluster.trackptSec(); // std:vector diff --git a/PWGEM/PhotonMeson/Core/PWGEMPhotonMesonCoreLinkDef.h b/PWGEM/PhotonMeson/Core/PWGEMPhotonMesonCoreLinkDef.h index ff03b41a75e..2a1a0ff7226 100644 --- a/PWGEM/PhotonMeson/Core/PWGEMPhotonMesonCoreLinkDef.h +++ b/PWGEM/PhotonMeson/Core/PWGEMPhotonMesonCoreLinkDef.h @@ -22,5 +22,6 @@ #pragma link C++ class EMCPhotonCut + ; #pragma link C++ class EMPhotonEventCut + ; #pragma link C++ class PairCut + ; +#pragma link C++ class EMBitFlags + ; #endif // PWGEM_PHOTONMESON_CORE_PWGEMPHOTONMESONCORELINKDEF_H_ diff --git a/PWGEM/PhotonMeson/Tasks/CMakeLists.txt b/PWGEM/PhotonMeson/Tasks/CMakeLists.txt index ef6371b5e94..29cbae1336b 100644 --- a/PWGEM/PhotonMeson/Tasks/CMakeLists.txt +++ b/PWGEM/PhotonMeson/Tasks/CMakeLists.txt @@ -151,6 +151,11 @@ o2physics_add_dpl_workflow(pi0-flow-emc PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(task-flow-reso + SOURCES taskFlowReso.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(prefilter-photon SOURCES prefilterPhoton.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore diff --git a/PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx b/PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx index d9d75bf8ade..91469523448 100644 --- a/PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx +++ b/PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx @@ -13,9 +13,11 @@ /// \brief Task to produce calibration values for EMCal /// \author M. Hemmer, marvin.hemmer@cern.ch +#include "PWGEM/PhotonMeson/Core/EMBitFlags.h" #include "PWGEM/PhotonMeson/Core/EMCPhotonCut.h" #include "PWGEM/PhotonMeson/Core/EMPhotonEventCut.h" #include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" +#include "PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" #include "PWGEM/PhotonMeson/Utils/EventHistograms.h" #include "PWGEM/PhotonMeson/Utils/emcalHistoDefinitions.h" @@ -156,33 +158,38 @@ struct CalibTaskEmc { } emccuts; V0PhotonCut fPCMPhotonCut; - struct : ConfigurableGroup { + struct : o2::framework::ConfigurableGroup { std::string prefix = "pcmcuts"; - Configurable cfgRequireV0WithITSTPC{"cfgRequireV0WithITSTPC", false, "flag to select V0s with ITS-TPC matched tracks"}; - Configurable cfgRequireV0WithITSonly{"cfgRequireV0WithITSonly", false, "flag to select V0s with ITSonly tracks"}; - Configurable cfgRequireV0WithTPConly{"cfgRequireV0WithTPConly", false, "flag to select V0s with TPConly tracks"}; - Configurable cfgMinPtV0{"cfgMinPtV0", 0.1, "min pT for v0 photons at PV"}; - Configurable cfgMaxPtV0{"cfgMaxPtV0", 1e+10, "max pT for v0 photons at PV"}; - Configurable cfgMinEtaV0{"cfgMinEtaV0", -0.8, "min eta for v0 photons at PV"}; - Configurable cfgMaxEtaV0{"cfgMaxEtaV0", 0.8, "max eta for v0 photons at PV"}; - Configurable cfgMinV0Radius{"cfgMinV0Radius", 4.0, "min v0 radius"}; - Configurable cfgMaxV0Radius{"cfgMaxV0Radius", 90.0, "max v0 radius"}; - Configurable cfgMaxAlphaAP{"cfgMaxAlphaAP", 0.95, "max alpha for AP cut"}; - Configurable cfgMaxQtAP{"cfgMaxQtAP", 0.01, "max qT for AP cut"}; - Configurable cfgMinCosPA{"cfgMinCosPA", 0.999, "min V0 CosPA"}; - Configurable cfgMaxPCA{"cfgMaxPCA", 1.5, "max distance btween 2 legs"}; - Configurable cfgMaxChi2KF{"cfgMaxChi2KF", 1e+10, "max chi2/ndf with KF"}; - Configurable cfgRejectV0OnITSib{"cfgRejectV0OnITSib", true, "flag to reject V0s on ITSib"}; - - Configurable cfgMinNClusterTPC{"cfgMinNClusterTPC", 0, "min ncluster tpc"}; - Configurable cfgMinNCrossedRows{"cfgMinNCrossedRows", 40, "min ncrossed rows"}; - Configurable cfgMaxFracSharedClusterTPC{"cfgMaxFracSharedClusterTPC", 999.f, "max fraction of shared clusters in TPC"}; - Configurable cfgMaxChi2TPC{"cfgMaxChi2TPC", 4.0, "max chi2/NclsTPC"}; - Configurable cfgMaxChi2ITS{"cfgMaxChi2ITS", 36.0, "max chi2/NclsITS"}; - Configurable cfgMinTPCNSigmaEl{"cfgMinTPCNSigmaEl", -3.0, "min. TPC n sigma for electron"}; - Configurable cfgMaxTPCNSigmaEl{"cfgMaxTPCNSigmaEl", +3.0, "max. TPC n sigma for electron"}; - Configurable cfgDisableITSOnly{"cfgDisableITSOnly", false, "flag to disable ITSonly tracks"}; - Configurable cfgDisableTPCOnly{"cfgDisableTPCOnly", false, "flag to disable TPConly tracks"}; + o2::framework::Configurable requireV0WithITSTPC{"requireV0WithITSTPC", false, "flag to enforce V0s have ITS and TPC"}; + o2::framework::Configurable requireV0WithITSOnly{"requireV0WithITSOnly", false, "flag to select V0s with ITSonly tracks"}; + o2::framework::Configurable requireV0WithTPCOnly{"requireV0WithTPCOnly", false, "flag to select V0s with TPConly tracks"}; + o2::framework::Configurable minPtV0{"minPtV0", 0.1, "min pT for v0 photons at PV"}; + o2::framework::Configurable maxPtV0{"maxPtV0", 1e+10, "max pT for v0 photons at PV"}; + o2::framework::Configurable minEtaV0{"minEtaV0", -0.8, "min eta for v0 photons at PV"}; + o2::framework::Configurable maxEtaV0{"maxEtaV0", 0.8, "max eta for v0 photons at PV"}; + o2::framework::Configurable minRV0{"minRV0", 4.0, "min v0 radius"}; + o2::framework::Configurable maxRV0{"maxRV0", 90.0, "max v0 radius"}; + o2::framework::Configurable maxAlphaAP{"maxAlphaAP", 0.95, "max alpha for AP cut"}; + o2::framework::Configurable maxQtAP{"maxQtAP", 0.01, "max qT for AP cut"}; + o2::framework::Configurable minCosPA{"minCosPA", 0.999, "min V0 CosPA"}; + o2::framework::Configurable maxPCA{"maxPCA", 1.5, "max distance btween 2 legs"}; + o2::framework::Configurable maxChi2KF{"maxChi2KF", 1.e+10f, "max chi2/ndf with KF"}; + o2::framework::Configurable rejectV0onITSib{"rejectV0onITSib", true, "flag to reject V0s on ITSib"}; + o2::framework::Configurable applyPrefilter{"applyPrefilter", false, "flag to apply prefilter to V0"}; + + o2::framework::Configurable minNClusterTPC{"minNClusterTPC", 0, "min NCluster TPC"}; + o2::framework::Configurable minNCrossedRowsTPC{"minNCrossedRowsTPC", 40, "min ncrossed rows in TPC"}; + o2::framework::Configurable minNCrossedRowsOverFindableClustersTPC{"minNCrossedRowsOverFindableClustersTPC", 0.8f, "min fraction of crossed rows over findable clusters in TPC"}; + o2::framework::Configurable maxFracSharedClustersTPC{"maxFracSharedClustersTPC", 999.f, "max fraction of shared clusters in TPC"}; + o2::framework::Configurable minNClusterITS{"minNClusterITS", 0, "min NCluster ITS"}; + o2::framework::Configurable minMeanClusterSizeITSob{"minMeanClusterSizeITSob", 0.f, "min ITSob"}; + o2::framework::Configurable maxMeanClusterSizeITSob{"maxMeanClusterSizeITSob", 16.f, "max ITSob"}; + o2::framework::Configurable macChi2TPC{"macChi2TPC", 4.f, "max chi2/NclsTPC"}; + o2::framework::Configurable macChi2ITS{"macChi2ITS", 36.f, "max chi2/NclsITS"}; + o2::framework::Configurable minTPCNSigmaEl{"minTPCNSigmaEl", -3.0f, "min. TPC n sigma for electron"}; + o2::framework::Configurable maxTPCNSigmaEl{"maxTPCNSigmaEl", +3.0f, "max. TPC n sigma for electron"}; + o2::framework::Configurable disableITSOnly{"disableITSOnly", false, "flag to disable ITSonly tracks"}; + o2::framework::Configurable disableTPCOnly{"disableTPCOnly", false, "flag to disable TPConly tracks"}; } pcmcuts; struct : ConfigurableGroup { @@ -206,10 +213,7 @@ struct CalibTaskEmc { int runNow = 0; int runBefore = -1; - // Filter clusterFilter = aod::skimmedcluster::time >= emccuts.cfgEMCminTime && aod::skimmedcluster::time <= emccuts.cfgEMCmaxTime && aod::skimmedcluster::m02 >= emccuts.cfgEMCminM02 && aod::skimmedcluster::m02 <= emccuts.cfgEMCmaxM02 && skimmedcluster::e >= emccuts.cfgEMCminE; - // Filter collisionFilter = (nabs(aod::collision::posZ) <= eventcuts.cfgZvtxMax) && (aod::evsel::ft0cOccupancyInTimeRange <= eventcuts.cfgFT0COccupancyMax) && (aod::evsel::ft0cOccupancyInTimeRange >= eventcuts.cfgFT0COccupancyMin); - // using FilteredEMCalPhotons = soa::Filtered>; - using EMCalPhotons = soa::Join; + using EMCalPhotons = soa::Join; using PCMPhotons = soa::Join; using Colls = soa::Join; @@ -218,8 +222,8 @@ struct CalibTaskEmc { MyEMH* emh1 = nullptr; MyEMH* emh2 = nullptr; - Preslice perCollisionEMC = aod::emccluster::emeventId; - Preslice perCollisionPCM = aod::v0photonkf::emeventId; + PresliceOptional perCollisionEMC = aod::emccluster::emeventId; + PresliceOptional perCollisionPCM = aod::v0photonkf::emeventId; using BinningType = ColumnBinningPolicy; BinningType binningOnPositions{{mixingConfig.cfgVtxBins, mixingConfig.cfgCentBins}, true}; @@ -301,35 +305,35 @@ struct CalibTaskEmc { fEMCCut.SetUseSecondaryTM(emccuts.emcUseSecondaryTM.value); // disables or enables secondary TM } - void DefinePCMCut() + void definePCMCut() { fPCMPhotonCut = V0PhotonCut("fPCMPhotonCut", "fPCMPhotonCut"); // for v0 - fPCMPhotonCut.SetV0PtRange(pcmcuts.cfgMinPtV0, pcmcuts.cfgMaxPtV0); - fPCMPhotonCut.SetV0EtaRange(pcmcuts.cfgMinEtaV0, pcmcuts.cfgMaxEtaV0); - fPCMPhotonCut.SetMinCosPA(pcmcuts.cfgMinCosPA); - fPCMPhotonCut.SetMaxPCA(pcmcuts.cfgMaxPCA); - fPCMPhotonCut.SetMaxChi2KF(pcmcuts.cfgMaxChi2KF); - fPCMPhotonCut.SetRxyRange(pcmcuts.cfgMinV0Radius, pcmcuts.cfgMaxV0Radius); - fPCMPhotonCut.SetAPRange(pcmcuts.cfgMaxAlphaAP, pcmcuts.cfgMaxQtAP); - fPCMPhotonCut.RejectITSib(pcmcuts.cfgRejectV0OnITSib); + fPCMPhotonCut.SetV0PtRange(pcmcuts.minPtV0, pcmcuts.maxPtV0); + fPCMPhotonCut.SetV0EtaRange(pcmcuts.minEtaV0, pcmcuts.maxEtaV0); + fPCMPhotonCut.SetMinCosPA(pcmcuts.minCosPA); + fPCMPhotonCut.SetMaxPCA(pcmcuts.maxPCA); + fPCMPhotonCut.SetMaxChi2KF(pcmcuts.maxChi2KF); + fPCMPhotonCut.SetRxyRange(pcmcuts.minRV0, pcmcuts.maxRV0); + fPCMPhotonCut.SetAPRange(pcmcuts.maxAlphaAP, pcmcuts.maxQtAP); + fPCMPhotonCut.RejectITSib(pcmcuts.rejectV0onITSib); // for track - fPCMPhotonCut.SetMinNClustersTPC(pcmcuts.cfgMinNClusterTPC); - fPCMPhotonCut.SetMinNCrossedRowsTPC(pcmcuts.cfgMinNCrossedRows); - fPCMPhotonCut.SetMinNCrossedRowsOverFindableClustersTPC(0.8); - fPCMPhotonCut.SetMaxFracSharedClustersTPC(pcmcuts.cfgMaxFracSharedClusterTPC); - fPCMPhotonCut.SetChi2PerClusterTPC(0.0, pcmcuts.cfgMaxChi2TPC); - fPCMPhotonCut.SetTPCNsigmaElRange(pcmcuts.cfgMinTPCNSigmaEl, pcmcuts.cfgMaxTPCNSigmaEl); - fPCMPhotonCut.SetChi2PerClusterITS(-1e+10, pcmcuts.cfgMaxChi2ITS); - fPCMPhotonCut.SetNClustersITS(0, 7); - fPCMPhotonCut.SetMeanClusterSizeITSob(0.0, 16.0); - fPCMPhotonCut.SetDisableITSonly(pcmcuts.cfgDisableITSOnly); - fPCMPhotonCut.SetDisableTPConly(pcmcuts.cfgDisableTPCOnly); - fPCMPhotonCut.SetRequireITSTPC(pcmcuts.cfgRequireV0WithITSTPC); - fPCMPhotonCut.SetRequireITSonly(pcmcuts.cfgRequireV0WithITSonly); - fPCMPhotonCut.SetRequireTPConly(pcmcuts.cfgRequireV0WithTPConly); + fPCMPhotonCut.SetMinNClustersTPC(pcmcuts.minNClusterTPC); + fPCMPhotonCut.SetMinNCrossedRowsTPC(pcmcuts.minNCrossedRowsTPC); + fPCMPhotonCut.SetMinNCrossedRowsOverFindableClustersTPC(pcmcuts.minNCrossedRowsOverFindableClustersTPC); + fPCMPhotonCut.SetMaxFracSharedClustersTPC(pcmcuts.maxFracSharedClustersTPC); + fPCMPhotonCut.SetChi2PerClusterTPC(0.f, pcmcuts.macChi2TPC); + fPCMPhotonCut.SetTPCNsigmaElRange(pcmcuts.minTPCNSigmaEl, pcmcuts.maxTPCNSigmaEl); + fPCMPhotonCut.SetChi2PerClusterITS(0.f, pcmcuts.macChi2ITS); + fPCMPhotonCut.SetNClustersITS(pcmcuts.minNClusterITS, 7); + fPCMPhotonCut.SetMeanClusterSizeITSob(pcmcuts.minMeanClusterSizeITSob, pcmcuts.maxMeanClusterSizeITSob); + fPCMPhotonCut.SetDisableITSonly(pcmcuts.disableITSOnly); + fPCMPhotonCut.SetDisableTPConly(pcmcuts.disableTPCOnly); + fPCMPhotonCut.SetRequireITSTPC(pcmcuts.requireV0WithITSTPC); + fPCMPhotonCut.SetRequireITSonly(pcmcuts.requireV0WithITSOnly); + fPCMPhotonCut.SetRequireTPConly(pcmcuts.requireV0WithTPCOnly); } void init(InitContext&) @@ -692,10 +696,14 @@ struct CalibTaskEmc { } // EMCal calibration same event - void processEMCalCalib(Colls const& collisions, EMCalPhotons const& clusters, PCMPhotons const&) + void processEMCalCalib(Colls const& collisions, EMCalPhotons const& clusters, PCMPhotons const&, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds) { float energyCorrectionFactor = 1.f; int nColl = 1; + + EMBitFlags emcFlags(clusters.size()); + fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds); + for (const auto& collision : collisions) { auto photonsPerCollision = clusters.sliceBy(perCollisionEMC, collision.globalIndex()); o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<0>(®istry, collision); @@ -731,7 +739,7 @@ struct CalibTaskEmc { } } for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photonsPerCollision, photonsPerCollision))) { - if (!(fEMCCut.IsSelected(g1)) || !(fEMCCut.IsSelected(g2))) { + if (!(emcFlags.test(g1.globalIndex())) || !(emcFlags.test(g2.globalIndex()))) { continue; } @@ -798,10 +806,14 @@ struct CalibTaskEmc { PROCESS_SWITCH(CalibTaskEmc, processEMCalCalib, "Process EMCal calibration same event", true); // EMCal calibration - void processEMCalPCMCalib(Colls const& collisions, EMCalPhotons const& clusters, PCMPhotons const& photons, aod::V0Legs const&) + void processEMCalPCMCalib(Colls const& collisions, EMCalPhotons const& clusters, PCMPhotons const& photons, aod::V0Legs const&, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds) { float energyCorrectionFactor = 1.f; int nColl = 1; + + EMBitFlags emcFlags(clusters.size()); + fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds); + for (const auto& collision : collisions) { o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<0>(®istry, collision); if (!(fEMEventCut.IsSelected(collision))) { @@ -839,7 +851,7 @@ struct CalibTaskEmc { } } for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(photonsEMCPerCollision, photonsPCMPerCollision))) { - if (!(fEMCCut.IsSelected(g1)) || !(fPCMPhotonCut.IsSelected(g2))) { + if (!(emcFlags.test(g1.globalIndex())) || !(fPCMPhotonCut.IsSelected(g2))) { continue; } @@ -901,9 +913,11 @@ struct CalibTaskEmc { PROCESS_SWITCH(CalibTaskEmc, processEMCalPCMCalib, "Process EMCal calibration using PCM-EMC same event", true); // EMCal calibration mixed event - void processEMCalCalibMixed(Colls const&, EMCalPhotons const&, PCMPhotons const&) + void processEMCalCalibMixed(Colls const&, EMCalPhotons const& clusters, PCMPhotons const&, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds) { float energyCorrectionFactor = 1.f; + EMBitFlags emcFlags(clusters.size()); + fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds); SameKindPair pair{binningOnPositions, mixingConfig.cfgMixingDepth, -1, &cache}; // indicates that 5 events should be mixed and under/overflow (-1) to be ignored @@ -926,7 +940,7 @@ struct CalibTaskEmc { runBefore = runNow; } for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(clusters1, clusters2))) { - if (!(fEMCCut.IsSelected(g1)) || !(fEMCCut.IsSelected(g2))) { + if (!(emcFlags.test(g1.globalIndex())) || !(emcFlags.test(g2.globalIndex()))) { continue; } // Cut edge clusters away, similar to rotation method to ensure same acceptance is used @@ -980,11 +994,11 @@ struct CalibTaskEmc { PROCESS_SWITCH(CalibTaskEmc, processEMCalCalibMixed, "Process EMCal calibration mixed event", false); // EMCal calibration - void processEMCalPCMCalibMixed(Colls const&, EMCalPhotons const&, PCMPhotons const&, aod::V0Legs const&) + void processEMCalPCMCalibMixed(Colls const&, EMCalPhotons const& clusters, PCMPhotons const&, aod::V0Legs const&, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds) { float energyCorrectionFactor = 1.f; - - LOG(info) << "Beginning of processEMCalPCMCalibMixed"; + EMBitFlags emcFlags(clusters.size()); + fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds); for (const auto& [c1, photonEMC, c2, photonPCM] : pairPCMEMC) { if (!(fEMEventCut.IsSelected(c1)) || !(fEMEventCut.IsSelected(c2))) { @@ -1005,7 +1019,7 @@ struct CalibTaskEmc { runBefore = runNow; } for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(photonEMC, photonPCM))) { - if (!(fEMCCut.IsSelected(g1)) || !(fPCMPhotonCut.IsSelected(g2))) { + if (!(emcFlags.test(g1.globalIndex())) || !(fPCMPhotonCut.IsSelected(g2))) { continue; } // Cut edge clusters away, similar to rotation method to ensure same acceptance is used diff --git a/PWGEM/PhotonMeson/Tasks/taskFlowReso.cxx b/PWGEM/PhotonMeson/Tasks/taskFlowReso.cxx new file mode 100644 index 00000000000..10731a2de42 --- /dev/null +++ b/PWGEM/PhotonMeson/Tasks/taskFlowReso.cxx @@ -0,0 +1,522 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file taskFlowReso.cxx +/// \brief Analysis task for flow resolution +/// \author M. Hemmer, marvin.hemmer@cern.ch + +#include "PWGEM/PhotonMeson/Core/EMPhotonEventCut.h" +#include "PWGEM/PhotonMeson/DataModel/gammaTables.h" +#include "PWGEM/PhotonMeson/Utils/EventHistograms.h" + +#include "Common/Core/EventPlaneHelper.h" +#include "Common/Core/RecoDecay.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::soa; +using namespace o2::aod::pwgem::photon; + +enum QvecEstimator { + FT0M = 0, + FT0A = 1, + FT0C, + TPCPos, + TPCNeg, + TPCTot, + FV0A +}; + +enum CentralityEstimator { + None = 0, + CFT0A = 1, + CFT0C, + CFT0M, + NCentralityEstimators +}; + +enum Harmonics { + kNone = 0, + kDirect = 1, + kElliptic = 2, + kTriangluar = 3, + kQuadrangular = 4, + kPentagonal = 5, + kHexagonal = 6, + kHeptagonal = 7, + kOctagonal = 8 +}; + +struct TaskFlowReso { + + // configurable for flow + Configurable harmonic{"harmonic", 2, "harmonic number"}; + Configurable qvecDetector{"qvecDetector", 0, "Detector for Q vector estimation (FT0M: 0, FT0A: 1, FT0C: 2, TPC Pos: 3, TPC Neg: 4, TPC Tot: 5, FV0A: 6)"}; + Configurable qvecSubADetector{"qvecSubADetector", 3, "Sub A Detector for Q vector estimation for resolution (FT0M: 0, FT0A: 1, FT0C: 2, TPC Pos: 3, TPC Neg: 4, TPC Tot: 5, FV0A: 6)"}; + Configurable qvecSubBDetector{"qvecSubBDetector", 4, "Sub B Detector for Q vector estimation for resolution (FT0M: 0, FT0A: 1, FT0C: 2, TPC Pos: 3, TPC Neg: 4, TPC Tot: 5, FV0A: 6)"}; + Configurable centEstimator{"centEstimator", 2, "Centrality estimation (FT0A: 1, FT0C: 2, FT0M: 3)"}; + Configurable saveEpResoHisto{"saveEpResoHisto", true, "Flag to save event plane resolution histogram"}; + Configurable saveSPResoHist{"saveSPResoHist", true, "Flag to save scalar product resolution histogram"}; + Configurable cfgMaxQVector{"cfgMaxQVector", 999999999.f, "Maximum allowed absolute QVector value."}; + + // configurable axis + ConfigurableAxis thnConfigAxisCent{"thnConfigAxisCent", {20, 0., 100.}, "centrality axis for the current event"}; + ConfigurableAxis thnConfigAxisCosNPhi{"thnConfigAxisCosNPhi", {100, -1., 1.}, "cos(n*phi) axis for the current event"}; + + EMPhotonEventCut fEMEventCut; + struct : ConfigurableGroup { + std::string prefix = "eventcuts"; + Configurable cfgZvtxMax{"cfgZvtxMax", 10.f, "max. Zvtx"}; + Configurable cfgRequireSel8{"cfgRequireSel8", true, "require sel8 in event cut"}; + Configurable cfgRequireFT0AND{"cfgRequireFT0AND", true, "require FT0AND in event cut"}; + Configurable cfgRequireNoTFB{"cfgRequireNoTFB", false, "require No time frame border in event cut"}; + Configurable cfgRequireNoITSROFB{"cfgRequireNoITSROFB", false, "require no ITS readout frame border in event cut"}; + Configurable cfgRequireNoSameBunchPileup{"cfgRequireNoSameBunchPileup", false, "require no same bunch pileup in event cut"}; + Configurable cfgRequireVertexITSTPC{"cfgRequireVertexITSTPC", false, "require Vertex ITSTPC in event cut"}; // ITS-TPC matched track contributes PV. + Configurable cfgRequireGoodZvtxFT0vsPV{"cfgRequireGoodZvtxFT0vsPV", false, "require good Zvtx between FT0 vs. PV in event cut"}; + Configurable cfgRequireEMCReadoutInMB{"cfgRequireEMCReadoutInMB", true, "require the EMC to be read out in an MB collision (kTVXinEMC)"}; + Configurable cfgRequireEMCHardwareTriggered{"cfgRequireEMCHardwareTriggered", false, "require the EMC to be hardware triggered (kEMC7 or kDMC7)"}; + Configurable cfgTrackOccupancyMin{"cfgTrackOccupancyMin", -1, "min. track occupancy"}; + Configurable cfgTrackOccupancyMax{"cfgTrackOccupancyMax", 1000000000, "max. track occupancy"}; + Configurable cfgFT0COccupancyMin{"cfgFT0COccupancyMin", -1, "min. FT0C occupancy"}; + Configurable cfgFT0COccupancyMax{"cfgFT0COccupancyMax", 1000000000, "max. FT0C occupancy"}; + Configurable cfgMinCent{"cfgMinCent", 0, "min. centrality (%)"}; + Configurable cfgMaxCent{"cfgMaxCent", 90, "max. centrality (%)"}; + Configurable onlyKeepWeightedEvents{"onlyKeepWeightedEvents", false, "flag to keep only weighted events (for JJ MCs) and remove all MB events (with weight = 1)"}; + } eventcuts; + + SliceCache cache; + EventPlaneHelper epHelper; + + using CollsWithQvecs = soa::Join; + using CollWithQvec = CollsWithQvecs::iterator; + + static constexpr std::size_t NQVecEntries = 6; + + HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; + + void defineEMEventCut() + { + fEMEventCut = EMPhotonEventCut("fEMEventCut", "fEMEventCut"); + fEMEventCut.SetRequireSel8(eventcuts.cfgRequireSel8); + fEMEventCut.SetRequireFT0AND(eventcuts.cfgRequireFT0AND); + fEMEventCut.SetZvtxRange(-eventcuts.cfgZvtxMax, +eventcuts.cfgZvtxMax); + fEMEventCut.SetRequireNoTFB(eventcuts.cfgRequireNoTFB); + fEMEventCut.SetRequireNoITSROFB(eventcuts.cfgRequireNoITSROFB); + fEMEventCut.SetRequireNoSameBunchPileup(eventcuts.cfgRequireNoSameBunchPileup); + fEMEventCut.SetRequireVertexITSTPC(eventcuts.cfgRequireVertexITSTPC); + fEMEventCut.SetRequireGoodZvtxFT0vsPV(eventcuts.cfgRequireGoodZvtxFT0vsPV); + fEMEventCut.SetRequireEMCReadoutInMB(eventcuts.cfgRequireEMCReadoutInMB); + fEMEventCut.SetRequireEMCHardwareTriggered(eventcuts.cfgRequireEMCHardwareTriggered); + } + + void init(InitContext&) + { + if (harmonic != kElliptic && harmonic != kTriangluar) { + LOG(info) << "Harmonic was set to " << harmonic << " but can only be 2 or 3!"; + } + + defineEMEventCut(); + o2::aod::pwgem::photonmeson::utils::eventhistogram::addEventHistograms(®istry); + + const AxisSpec thnAxisCent{thnConfigAxisCent, "Centrality (%)"}; + const AxisSpec thnAxisCosNPhi{thnConfigAxisCosNPhi, Form("cos(%d#varphi)", harmonic.value)}; + const AxisSpec thAxisPsi{360 / harmonic.value, -(1. / static_cast(harmonic.value)) * std::numbers::pi_v, (1. / static_cast(harmonic.value)) * std::numbers::pi_v, Form("#Psi_{%d}", harmonic.value)}; + const AxisSpec thAxisCN{8, 0.5, 8.5, "#it{c}_{n}"}; + const AxisSpec thAxisSN{8, 0.5, 8.5, "#it{s}_{n}"}; + const AxisSpec thAxisAzimuth{360, -o2::constants::math::PI, o2::constants::math::PI, "#it{#varphi} (rad)"}; + + if (saveSPResoHist.value) { + registry.add("spReso/hSpResoFT0cFT0a", "hSpResoFT0cFT0a; centrality; Q_{FT0c} #bullet Q_{FT0a}", HistType::kTProfile, {thnAxisCent}); + registry.add("spReso/hSpResoFT0cTPCpos", "hSpResoFT0cTPCpos; centrality; Q_{FT0c} #bullet Q_{TPCpos}", HistType::kTProfile, {thnAxisCent}); + registry.add("spReso/hSpResoFT0cTPCneg", "hSpResoFT0cTPCneg; centrality; Q_{FT0c} #bullet Q_{TPCneg}", HistType::kTProfile, {thnAxisCent}); + registry.add("spReso/hSpResoFT0cTPCtot", "hSpResoFT0cTPCtot; centrality; Q_{FT0c} #bullet Q_{TPCtot}", HistType::kTProfile, {thnAxisCent}); + registry.add("spReso/hSpResoFT0aTPCpos", "hSpResoFT0aTPCpos; centrality; Q_{FT0a} #bullet Q_{TPCpos}", HistType::kTProfile, {thnAxisCent}); + registry.add("spReso/hSpResoFT0aTPCneg", "hSpResoFT0aTPCneg; centrality; Q_{FT0a} #bullet Q_{TPCneg}", HistType::kTProfile, {thnAxisCent}); + registry.add("spReso/hSpResoFT0aTPCtot", "hSpResoFT0aTPCtot; centrality; Q_{FT0m} #bullet Q_{TPCtot}", HistType::kTProfile, {thnAxisCent}); + registry.add("spReso/hSpResoFT0mTPCpos", "hSpResoFT0mTPCpos; centrality; Q_{FT0m} #bullet Q_{TPCpos}", HistType::kTProfile, {thnAxisCent}); + registry.add("spReso/hSpResoFT0mTPCneg", "hSpResoFT0mTPCneg; centrality; Q_{FT0m} #bullet Q_{TPCneg}", HistType::kTProfile, {thnAxisCent}); + registry.add("spReso/hSpResoFT0mTPCtot", "hSpResoFT0mTPCtot; centrality; Q_{FT0m} #bullet Q_{TPCtot}", HistType::kTProfile, {thnAxisCent}); + registry.add("spReso/hSpResoTPCposTPCneg", "hSpResoTPCposTPCneg; centrality; Q_{TPCpos} #bullet Q_{TPCneg}", HistType::kTProfile, {thnAxisCent}); + registry.add("spReso/hSpResoFV0aFT0c", "hSpResoFV0aFT0c; centrality; Q_{FV0a} #bullet Q_{FT0c}", HistType::kTProfile, {thnAxisCent}); + registry.add("spReso/hSpResoFV0aTPCpos", "hSpResoFV0aTPCpos; centrality; Q_{FV0a} #bullet Q_{TPCpos}", HistType::kTProfile, {thnAxisCent}); + registry.add("spReso/hSpResoFV0aTPCneg", "hSpResoFV0aTPCneg; centrality; Q_{FV0a} #bullet Q_{TPCneg}", HistType::kTProfile, {thnAxisCent}); + registry.add("spReso/hSpResoFV0aTPCtot", "hSpResoFV0aTPCtot; centrality; Q_{FV0a} #bullet Q_{TPCtot}", HistType::kTProfile, {thnAxisCent}); + } + + if (saveEpResoHisto.value) { + registry.add("hEventPlaneAngleFT0C", "hEventPlaneAngleFT0C", HistType::kTH2D, {thnAxisCent, thAxisPsi}); + registry.add("hEventPlaneAngleFT0A", "hEventPlaneAngleFT0A", HistType::kTH2D, {thnAxisCent, thAxisPsi}); + registry.add("hEventPlaneAngleFT0M", "hEventPlaneAngleFT0M", HistType::kTH2D, {thnAxisCent, thAxisPsi}); + registry.add("hEventPlaneAngleFV0A", "hEventPlaneAngleFV0A", HistType::kTH2D, {thnAxisCent, thAxisPsi}); + registry.add("hEventPlaneAngleTPCpos", "hEventPlaneAngleTPCpos", HistType::kTH2D, {thnAxisCent, thAxisPsi}); + registry.add("hEventPlaneAngleTPCneg", "hEventPlaneAngleTPCneg", HistType::kTH2D, {thnAxisCent, thAxisPsi}); + registry.add("epReso/hEpResoFT0cFT0a", "hEpResoFT0cFT0a; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); + registry.add("epReso/hEpResoFT0cTPCpos", "hEpResoFT0cTPCpos; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); + registry.add("epReso/hEpResoFT0cTPCneg", "hEpResoFT0cTPCneg; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); + registry.add("epReso/hEpResoFT0cTPCtot", "hEpResoFT0cTPCtot; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); + registry.add("epReso/hEpResoFT0aTPCpos", "hEpResoFT0aTPCpos; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); + registry.add("epReso/hEpResoFT0aTPCneg", "hEpResoFT0aTPCneg; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); + registry.add("epReso/hEpResoFT0aTPCtot", "hEpResoFT0aTPCtot; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); + registry.add("epReso/hEpResoFT0mTPCpos", "hEpResoFT0mTPCpos; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); + registry.add("epReso/hEpResoFT0mTPCneg", "hEpResoFT0mTPCneg; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); + registry.add("epReso/hEpResoFT0mTPCtot", "hEpResoFT0mTPCtot; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); + registry.add("epReso/hEpResoTPCposTPCneg", "hEpResoTPCposTPCneg; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); + registry.add("epReso/hEpResoFV0aFT0c", "hEpResoFV0aFT0c; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); + registry.add("epReso/hEpResoFV0aTPCpos", "hEpResoFV0aTPCpos; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); + registry.add("epReso/hEpResoFV0aTPCneg", "hEpResoFV0aTPCneg; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); + registry.add("epReso/hEpResoFV0aTPCtot", "hEpResoFV0aTPCtot; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); + registry.add("epReso/hEpCosCoefficientsFT0c", "hEpCosCoefficientsFT0c; centrality; c_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisCN}); + registry.add("epReso/hEpSinCoefficientsFT0c", "hEpSinCoefficientsFT0c; centrality; s_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisSN}); + registry.add("epReso/hEpCosCoefficientsFT0a", "hEpCosCoefficientsFT0a; centrality; c_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisCN}); + registry.add("epReso/hEpSinCoefficientsFT0a", "hEpSinCoefficientsFT0a; centrality; s_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisSN}); + registry.add("epReso/hEpCosCoefficientsFV0a", "hEpCosCoefficientsFV0a; centrality; c_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisCN}); + registry.add("epReso/hEpSinCoefficientsFV0a", "hEpSinCoefficientsFV0a; centrality; s_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisSN}); + registry.add("epReso/hEpCosCoefficientsFT0m", "hEpCosCoefficientsFT0m; centrality; c_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisCN}); + registry.add("epReso/hEpSinCoefficientsFT0m", "hEpSinCoefficientsFT0m; centrality; s_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisSN}); + registry.add("epReso/hEpCosCoefficientsTPCpos", "hEpCosCoefficientsTPCpos; centrality; c_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisCN}); + registry.add("epReso/hEpSinCoefficientsTPCpos", "hEpSinCoefficientsTPCpos; centrality; s_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisSN}); + registry.add("epReso/hEpCosCoefficientsTPCneg", "hEpCosCoefficientsTPCneg; centrality; c_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisCN}); + registry.add("epReso/hEpSinCoefficientsTPCneg", "hEpSinCoefficientsTPCneg; centrality; s_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisSN}); + registry.add("epReso/hEpCosCoefficientsTPCTots", "hEpCosCoefficientsTPCTots; centrality; c_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisCN}); + registry.add("epReso/hEpSinCoefficientsTPCTots", "hEpSinCoefficientsTPCTots; centrality; s_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisSN}); + registry.add("QVector/hQVecMeanRVsPhiFT0a", "hQVecMeanRVsPhiFT0a; centrality; #it{#varphi} (rad), <#it{r}> (a.u.)", HistType::kTProfile2D, {thnAxisCent, thAxisAzimuth}); + registry.add("QVector/hQVecMeanRVsPhiFT0c", "hQVecMeanRVsPhiFT0c; centrality; #it{#varphi} (rad), <#it{r}> (a.u.)", HistType::kTProfile2D, {thnAxisCent, thAxisAzimuth}); + registry.add("QVector/hQVecMeanRVsPhiFT0m", "hQVecMeanRVsPhiFT0m; centrality; #it{#varphi} (rad), <#it{r}> (a.u.)", HistType::kTProfile2D, {thnAxisCent, thAxisAzimuth}); + registry.add("QVector/hQVecMeanRVsPhiFV0a", "hQVecMeanRVsPhiFV0a; centrality; #it{#varphi} (rad), <#it{r}> (a.u.)", HistType::kTProfile2D, {thnAxisCent, thAxisAzimuth}); + registry.add("QVector/hQVecMeanRVsPhiTPCpos", "hQVecMeanRVsPhiTPCpos; centrality; #it{#varphi} (rad), <#it{r}> (a.u.)", HistType::kTProfile2D, {thnAxisCent, thAxisAzimuth}); + registry.add("QVector/hQVecMeanRVsPhiTPCneg", "hQVecMeanRVsPhiTPCneg; centrality; #it{#varphi} (rad), <#it{r}> (a.u.)", HistType::kTProfile2D, {thnAxisCent, thAxisAzimuth}); + registry.add("QVector/hQVecMeanRVsPhiTPCTot", "hQVecMeanRVsPhiTPCTot; centrality; #it{#varphi} (rad), <#it{r}> (a.u.)", HistType::kTProfile2D, {thnAxisCent, thAxisAzimuth}); + } + }; // end init + + /// Compute the delta psi in the range [0, pi/harmonic] + /// \param psi1 is the first angle + /// \param psi2 is the second angle + float getDeltaPsiInRange(float psi1, float psi2) + { + float deltaPsi = psi1 - psi2; + return RecoDecay::constrainAngle(deltaPsi, 0.f, harmonic); + } + + /// Get the centrality + /// \param collision is the collision with the centrality information + template + float getCentrality(TCollision const& collision) + { + float cent = -999.; + switch (centEstimator) { + case CentralityEstimator::CFT0M: + cent = collision.centFT0M(); + break; + case CentralityEstimator::CFT0A: + cent = collision.centFT0A(); + break; + case CentralityEstimator::CFT0C: + cent = collision.centFT0C(); + break; + default: + LOG(warning) << "Centrality estimator not valid. Possible values are T0M, T0A, T0C. Fallback to T0C"; + cent = collision.centFT0C(); + break; + } + return cent; + } + + /// Get all used Q vector + /// \param collision is the collision with the Q vector information + template + std::array getAllQvec(TCollision const& collision) + { + // Retrieve the Q vectors using the helper function for each detector + auto [xQVecMain, yQVecMain] = getQvec(collision, qvecDetector); + auto [xQVecSubA, yQVecSubA] = getQvec(collision, qvecSubADetector); + auto [xQVecSubB, yQVecSubB] = getQvec(collision, qvecSubBDetector); + + return {xQVecMain, yQVecMain, xQVecSubA, yQVecSubA, xQVecSubB, yQVecSubB}; + } + + /// Get the Q vector + /// \param collision is the collision with the Q vector information + template + std::pair getQvec(TCollision const& collision, int detector) + { + float xQVec = -999.f; + float yQVec = -999.f; + + switch (detector) { + case QvecEstimator::FT0M: + if (harmonic == kElliptic) { + xQVec = collision.q2xft0m(); + yQVec = collision.q2yft0m(); + } else if (harmonic == kTriangluar) { + xQVec = collision.q3xft0m(); + yQVec = collision.q3yft0m(); + } + break; + case QvecEstimator::FT0A: + if (harmonic == kElliptic) { + xQVec = collision.q2xft0a(); + yQVec = collision.q2yft0a(); + } else if (harmonic == kTriangluar) { + xQVec = collision.q3xft0a(); + yQVec = collision.q3yft0a(); + } + break; + case QvecEstimator::FT0C: + if (harmonic == kElliptic) { + xQVec = collision.q2xft0c(); + yQVec = collision.q2yft0c(); + } else if (harmonic == kTriangluar) { + xQVec = collision.q3xft0c(); + yQVec = collision.q3yft0c(); + } + break; + case QvecEstimator::TPCPos: + if (harmonic == kElliptic) { + xQVec = collision.q2xbpos(); + yQVec = collision.q2ybpos(); + } else if (harmonic == kTriangluar) { + xQVec = collision.q3xbpos(); + yQVec = collision.q3ybpos(); + } + break; + case QvecEstimator::TPCNeg: + if (harmonic == kElliptic) { + xQVec = collision.q2xbneg(); + yQVec = collision.q2ybneg(); + } else if (harmonic == kTriangluar) { + xQVec = collision.q3xbneg(); + yQVec = collision.q3ybneg(); + } + break; + case QvecEstimator::TPCTot: + if (harmonic == kElliptic) { + xQVec = collision.q2xbtot(); + yQVec = collision.q2ybtot(); + } else if (harmonic == kTriangluar) { + xQVec = collision.q3xbtot(); + yQVec = collision.q3ybtot(); + } + break; + case QvecEstimator::FV0A: + if (harmonic == kElliptic) { + xQVec = collision.q2xfv0a(); + yQVec = collision.q2yfv0a(); + } else if (harmonic == kTriangluar) { + xQVec = collision.q3xfv0a(); + yQVec = collision.q3yfv0a(); + } + break; + default: + LOG(warning) << "Q vector estimator not valid. Falling back to FT0M"; + if (harmonic == kElliptic) { + xQVec = collision.q2xft0m(); + yQVec = collision.q2yft0m(); + } else if (harmonic == kTriangluar) { + xQVec = collision.q3xft0m(); + yQVec = collision.q3yft0m(); + } + break; + } + return {xQVec, yQVec}; + } + + /// Check if the QVector values are within reasonable range + /// \param collision is the collision with the Q vector information + bool isQvecGood(std::array const& QVecs) + { + bool isgood = true; + for (const auto& QVec : QVecs) { + if (std::fabs(QVec) > cfgMaxQVector) { + isgood = false; + break; + } + } + return isgood; + } + + // Resolution + void processResolution(CollWithQvec const& collision) + { + o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<0>(®istry, collision); + if (!(fEMEventCut.IsSelected(collision))) { + // no selection on the centrality is applied on purpose to allow for the resolution study in post-processing + return; + } + if (!(eventcuts.cfgFT0COccupancyMin <= collision.ft0cOccupancyInTimeRange() && collision.ft0cOccupancyInTimeRange() < eventcuts.cfgFT0COccupancyMax)) { + return; + } + float cent = getCentrality(collision); + if (cent < eventcuts.cfgMinCent || cent > eventcuts.cfgMaxCent) { + // event selection + return; + } + if (!isQvecGood(getAllQvec(collision))) { + // selection based on QVector + return; + } + o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<1>(®istry, collision); + registry.fill(HIST("Event/before/hCollisionCounter"), 12.0); // accepted + registry.fill(HIST("Event/after/hCollisionCounter"), 12.0); // accepted + + float centrality = getCentrality(collision); // centrality not updated in the rejection mask function + float xQVecFT0a = -999.f; + float yQVecFT0a = -999.f; + float xQVecFT0c = -999.f; + float yQVecFT0c = -999.f; + float xQVecFT0m = -999.f; + float yQVecFT0m = -999.f; + float xQVecBPos = -999.f; + float yQVecBPos = -999.f; + float xQVecBNeg = -999.f; + float yQVecBNeg = -999.f; + float xQVecBTot = -999.f; + float yQVecBTot = -999.f; + float xQVecFV0a = -999.f; + float yQVecFV0a = -999.f; + if (harmonic == kElliptic) { + xQVecFT0a = collision.q2xft0a(); + yQVecFT0a = collision.q2yft0a(); + xQVecFT0c = collision.q2xft0c(); + yQVecFT0c = collision.q2yft0c(); + xQVecFT0m = collision.q2xft0m(); + yQVecFT0m = collision.q2yft0m(); + xQVecBPos = collision.q2xbpos(); + yQVecBPos = collision.q2ybpos(); + xQVecBNeg = collision.q2xbneg(); + yQVecBNeg = collision.q2ybneg(); + xQVecBTot = collision.q2xbtot(); + yQVecBTot = collision.q2ybtot(); + xQVecFV0a = collision.q2xfv0a(); + yQVecFV0a = collision.q2yfv0a(); + } else if (harmonic == kTriangluar) { + xQVecFT0a = collision.q3xft0a(); + yQVecFT0a = collision.q3yft0a(); + xQVecFT0c = collision.q3xft0c(); + yQVecFT0c = collision.q3yft0c(); + xQVecFT0m = collision.q3xft0m(); + yQVecFT0m = collision.q3yft0m(); + xQVecBPos = collision.q3xbpos(); + yQVecBPos = collision.q3ybpos(); + xQVecBNeg = collision.q3xbneg(); + yQVecBNeg = collision.q3ybneg(); + xQVecBTot = collision.q3xbtot(); + yQVecBTot = collision.q3ybtot(); + xQVecFV0a = collision.q3xfv0a(); + yQVecFV0a = collision.q3yfv0a(); + } + + if (saveSPResoHist) { + registry.fill(HIST("spReso/hSpResoFT0cFT0a"), centrality, xQVecFT0c * xQVecFT0a + yQVecFT0c * yQVecFT0a); + registry.fill(HIST("spReso/hSpResoFT0cTPCpos"), centrality, xQVecFT0c * xQVecBPos + yQVecFT0c * yQVecBPos); + registry.fill(HIST("spReso/hSpResoFT0cTPCneg"), centrality, xQVecFT0c * xQVecBNeg + yQVecFT0c * yQVecBNeg); + registry.fill(HIST("spReso/hSpResoFT0cTPCtot"), centrality, xQVecFT0c * xQVecBTot + yQVecFT0c * yQVecBTot); + registry.fill(HIST("spReso/hSpResoFT0aTPCpos"), centrality, xQVecFT0a * xQVecBPos + yQVecFT0a * yQVecBPos); + registry.fill(HIST("spReso/hSpResoFT0aTPCneg"), centrality, xQVecFT0a * xQVecBNeg + yQVecFT0a * yQVecBNeg); + registry.fill(HIST("spReso/hSpResoFT0aTPCtot"), centrality, xQVecFT0a * xQVecBTot + yQVecFT0a * yQVecBTot); + registry.fill(HIST("spReso/hSpResoFT0mTPCpos"), centrality, xQVecFT0m * xQVecBPos + yQVecFT0m * yQVecBPos); + registry.fill(HIST("spReso/hSpResoFT0mTPCneg"), centrality, xQVecFT0m * xQVecBNeg + yQVecFT0m * yQVecBNeg); + registry.fill(HIST("spReso/hSpResoFT0mTPCtot"), centrality, xQVecFT0m * xQVecBTot + yQVecFT0m * yQVecBTot); + registry.fill(HIST("spReso/hSpResoTPCposTPCneg"), centrality, xQVecBPos * xQVecBNeg + yQVecBPos * yQVecBNeg); + registry.fill(HIST("spReso/hSpResoFV0aFT0c"), centrality, xQVecFV0a * xQVecFT0c + yQVecFV0a * yQVecFT0c); + registry.fill(HIST("spReso/hSpResoFV0aTPCpos"), centrality, xQVecFV0a * xQVecBPos + yQVecFV0a * yQVecBPos); + registry.fill(HIST("spReso/hSpResoFV0aTPCneg"), centrality, xQVecFV0a * xQVecBNeg + yQVecFV0a * yQVecBNeg); + registry.fill(HIST("spReso/hSpResoFV0aTPCtot"), centrality, xQVecFV0a * xQVecBTot + yQVecFV0a * yQVecBTot); + } + + if (saveEpResoHisto) { + float epFT0a = epHelper.GetEventPlane(xQVecFT0a, yQVecFT0a, harmonic); + float epFT0c = epHelper.GetEventPlane(xQVecFT0c, yQVecFT0c, harmonic); + float epFT0m = epHelper.GetEventPlane(xQVecFT0m, yQVecFT0m, harmonic); + float epBPoss = epHelper.GetEventPlane(xQVecBPos, yQVecBPos, harmonic); + float epBNegs = epHelper.GetEventPlane(xQVecBNeg, yQVecBNeg, harmonic); + float epBTots = epHelper.GetEventPlane(xQVecBTot, yQVecBTot, harmonic); + float epFV0a = epHelper.GetEventPlane(xQVecFV0a, yQVecFV0a, harmonic); + + registry.fill(HIST("hEventPlaneAngleFT0M"), centrality, epFT0m); + registry.fill(HIST("hEventPlaneAngleFT0A"), centrality, epFT0c); + registry.fill(HIST("hEventPlaneAngleFT0C"), centrality, epFT0a); + registry.fill(HIST("hEventPlaneAngleTPCpos"), centrality, epBPoss); + registry.fill(HIST("hEventPlaneAngleTPCneg"), centrality, epBNegs); + registry.fill(HIST("hEventPlaneAngleFV0A"), centrality, epFV0a); + + registry.fill(HIST("QVector/hQVecMeanRVsPhiFT0a"), centrality, std::atan2(yQVecFT0a, xQVecFT0a), std::hypot(xQVecFT0a, yQVecFT0a)); + registry.fill(HIST("QVector/hQVecMeanRVsPhiFT0c"), centrality, std::atan2(yQVecFT0c, xQVecFT0c), std::hypot(xQVecFT0c, yQVecFT0c)); + registry.fill(HIST("QVector/hQVecMeanRVsPhiFT0m"), centrality, std::atan2(yQVecFT0m, xQVecFT0m), std::hypot(xQVecFT0m, yQVecFT0m)); + registry.fill(HIST("QVector/hQVecMeanRVsPhiTPCpos"), centrality, std::atan2(yQVecBPos, xQVecBPos), std::hypot(xQVecBPos, yQVecBPos)); + registry.fill(HIST("QVector/hQVecMeanRVsPhiTPCneg"), centrality, std::atan2(yQVecBNeg, xQVecBNeg), std::hypot(xQVecBNeg, yQVecBNeg)); + registry.fill(HIST("QVector/hQVecMeanRVsPhiTPCTot"), centrality, std::atan2(yQVecBTot, xQVecBTot), std::hypot(xQVecBTot, yQVecBTot)); + registry.fill(HIST("QVector/hQVecMeanRVsPhiFV0a"), centrality, std::atan2(yQVecFV0a, xQVecFV0a), std::hypot(xQVecFV0a, yQVecFV0a)); + + registry.fill(HIST("epReso/hEpResoFT0cFT0a"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0c, epFT0a))); + registry.fill(HIST("epReso/hEpResoFT0cTPCpos"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0c, epBPoss))); + registry.fill(HIST("epReso/hEpResoFT0cTPCneg"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0c, epBNegs))); + registry.fill(HIST("epReso/hEpResoFT0cTPCtot"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0c, epBTots))); + registry.fill(HIST("epReso/hEpResoFT0aTPCpos"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0a, epBPoss))); + registry.fill(HIST("epReso/hEpResoFT0aTPCneg"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0a, epBNegs))); + registry.fill(HIST("epReso/hEpResoFT0aTPCtot"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0a, epBTots))); + registry.fill(HIST("epReso/hEpResoFT0mTPCpos"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0m, epBPoss))); + registry.fill(HIST("epReso/hEpResoFT0mTPCneg"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0m, epBNegs))); + registry.fill(HIST("epReso/hEpResoFT0mTPCtot"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0m, epBTots))); + registry.fill(HIST("epReso/hEpResoTPCposTPCneg"), centrality, std::cos(harmonic * getDeltaPsiInRange(epBPoss, epBNegs))); + registry.fill(HIST("epReso/hEpResoFV0aFT0c"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFV0a, epFT0c))); + registry.fill(HIST("epReso/hEpResoFV0aTPCpos"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFV0a, epBPoss))); + registry.fill(HIST("epReso/hEpResoFV0aTPCneg"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFV0a, epBNegs))); + registry.fill(HIST("epReso/hEpResoFV0aTPCtot"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFV0a, epBTots))); + for (int n = 1; n <= kOctagonal; n++) { + registry.fill(HIST("epReso/hEpCosCoefficientsFT0c"), centrality, n, std::cos(n * harmonic * epFT0c)); + registry.fill(HIST("epReso/hEpSinCoefficientsFT0c"), centrality, n, std::sin(n * harmonic * epFT0c)); + registry.fill(HIST("epReso/hEpCosCoefficientsFT0a"), centrality, n, std::cos(n * harmonic * epFT0a)); + registry.fill(HIST("epReso/hEpSinCoefficientsFT0a"), centrality, n, std::sin(n * harmonic * epFT0a)); + registry.fill(HIST("epReso/hEpCosCoefficientsFT0m"), centrality, n, std::cos(n * harmonic * epFT0m)); + registry.fill(HIST("epReso/hEpSinCoefficientsFT0m"), centrality, n, std::sin(n * harmonic * epFT0m)); + registry.fill(HIST("epReso/hEpCosCoefficientsTPCpos"), centrality, n, std::cos(n * harmonic * epBPoss)); + registry.fill(HIST("epReso/hEpSinCoefficientsTPCpos"), centrality, n, std::sin(n * harmonic * epBPoss)); + registry.fill(HIST("epReso/hEpCosCoefficientsTPCneg"), centrality, n, std::cos(n * harmonic * epBNegs)); + registry.fill(HIST("epReso/hEpSinCoefficientsTPCneg"), centrality, n, std::sin(n * harmonic * epBNegs)); + registry.fill(HIST("epReso/hEpCosCoefficientsTPCTots"), centrality, n, std::cos(n * harmonic * epBTots)); + registry.fill(HIST("epReso/hEpSinCoefficientsTPCTots"), centrality, n, std::sin(n * harmonic * epBTots)); + registry.fill(HIST("epReso/hEpCosCoefficientsFV0a"), centrality, n, std::cos(n * harmonic * epFV0a)); + registry.fill(HIST("epReso/hEpSinCoefficientsFV0a"), centrality, n, std::sin(n * harmonic * epFV0a)); + } + } + } + PROCESS_SWITCH(TaskFlowReso, processResolution, "Process resolution", true); + +}; // End struct TaskFlowReso + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} diff --git a/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx b/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx index 9431f833041..517c514987e 100644 --- a/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx @@ -13,8 +13,11 @@ /// \brief Analysis task for neutral pion flow with EMCal /// \author M. Hemmer, marvin.hemmer@cern.ch +#include "PWGEM/PhotonMeson/Core/EMBitFlags.h" #include "PWGEM/PhotonMeson/Core/EMCPhotonCut.h" #include "PWGEM/PhotonMeson/Core/EMPhotonEventCut.h" +#include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" +#include "PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" #include "PWGEM/PhotonMeson/Utils/EventHistograms.h" @@ -54,9 +57,9 @@ #include #include +#include #include #include -#include #include #include #include @@ -116,13 +119,11 @@ struct TaskPi0FlowEMC { Configurable qvecSubADetector{"qvecSubADetector", 3, "Sub A Detector for Q vector estimation for resolution (FT0M: 0, FT0A: 1, FT0C: 2, TPC Pos: 3, TPC Neg: 4, TPC Tot: 5, FV0A: 6)"}; Configurable qvecSubBDetector{"qvecSubBDetector", 4, "Sub B Detector for Q vector estimation for resolution (FT0M: 0, FT0A: 1, FT0C: 2, TPC Pos: 3, TPC Neg: 4, TPC Tot: 5, FV0A: 6)"}; Configurable centEstimator{"centEstimator", 2, "Centrality estimation (FT0A: 1, FT0C: 2, FT0M: 3)"}; - Configurable saveEpResoHisto{"saveEpResoHisto", false, "Flag to save event plane resolution histogram"}; - Configurable saveSPResoHist{"saveSPResoHist", false, "Flag to save scalar product resolution histogram"}; Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - Configurable cfgDoRotation{"cfgDoRotation", true, "Flag to enable rotation background method"}; + Configurable cfgDoRotation{"cfgDoRotation", false, "Flag to enable rotation background method"}; Configurable cfgDownsampling{"cfgDownsampling", 1, "Calculate rotation background only for every collision"}; - Configurable cfgEMCalMapLevelBackground{"cfgEMCalMapLevelBackground", 2, "Different levels of correction for the background, the smaller number includes the level of the higher number (4: none, 3: only inside EMCal, 2: exclude bad channels, 1: remove edges)"}; - Configurable cfgEMCalMapLevelSameEvent{"cfgEMCalMapLevelSameEvent", 1, "Different levels of correction for the same event, the smaller number includes the level of the higher number (4: none, 3: only inside EMCal, 2: exclude bad channels, 1: remove edges)"}; + Configurable cfgEMCalMapLevelBackground{"cfgEMCalMapLevelBackground", 4, "Different levels of correction for the background, the smaller number includes the level of the higher number (4: none, 3: only inside EMCal, 2: exclude bad channels, 1: remove edges)"}; + Configurable cfgEMCalMapLevelSameEvent{"cfgEMCalMapLevelSameEvent", 4, "Different levels of correction for the same event, the smaller number includes the level of the higher number (4: none, 3: only inside EMCal, 2: exclude bad channels, 1: remove edges)"}; Configurable cfgRotAngle{"cfgRotAngle", std::move(const_cast(o2::constants::math::PIHalf)), "Angle used for the rotation method"}; Configurable cfgDistanceToEdge{"cfgDistanceToEdge", 1, "Distance to edge in cells required for rotated cluster to be accepted"}; Configurable cfgDoM02{"cfgDoM02", false, "Flag to enable flow vs M02 for single photons"}; @@ -153,8 +154,6 @@ struct TaskPi0FlowEMC { Configurable cfgRequireGoodZvtxFT0vsPV{"cfgRequireGoodZvtxFT0vsPV", false, "require good Zvtx between FT0 vs. PV in event cut"}; Configurable cfgRequireEMCReadoutInMB{"cfgRequireEMCReadoutInMB", true, "require the EMC to be read out in an MB collision (kTVXinEMC)"}; Configurable cfgRequireEMCHardwareTriggered{"cfgRequireEMCHardwareTriggered", false, "require the EMC to be hardware triggered (kEMC7 or kDMC7)"}; - Configurable cfgTrackOccupancyMin{"cfgTrackOccupancyMin", -1, "min. track occupancy"}; - Configurable cfgTrackOccupancyMax{"cfgTrackOccupancyMax", 1000000000, "max. track occupancy"}; Configurable cfgFT0COccupancyMin{"cfgFT0COccupancyMin", -1, "min. FT0C occupancy"}; Configurable cfgFT0COccupancyMax{"cfgFT0COccupancyMax", 1000000000, "max. FT0C occupancy"}; Configurable cfgMinCent{"cfgMinCent", 0, "min. centrality (%)"}; @@ -165,7 +164,7 @@ struct TaskPi0FlowEMC { EMCPhotonCut fEMCCut; struct : ConfigurableGroup { std::string prefix = "emccuts"; - Configurable clusterDefinition{"clusterDefinition", "kV3Default", "Clusterizer to be selected, e.g. V3Default"}; + Configurable clusterDefinition{"clusterDefinition", "kV3MostSplitSmallestTimeDiff", "Clusterizer to be selected, e.g. V3Default"}; Configurable cfgEMCminTime{"cfgEMCminTime", -25., "Minimum cluster time for EMCal time cut"}; Configurable cfgEMCmaxTime{"cfgEMCmaxTime", +30., "Maximum cluster time for EMCal time cut"}; Configurable cfgEMCminM02{"cfgEMCminM02", 0.1, "Minimum M02 for EMCal M02 cut"}; @@ -183,6 +182,41 @@ struct TaskPi0FlowEMC { Configurable cfgEnableQA{"cfgEnableQA", false, "flag to turn QA plots on/off"}; } emccuts; + V0PhotonCut fV0PhotonCut; + struct : o2::framework::ConfigurableGroup { + std::string prefix = "PCMcuts"; + o2::framework::Configurable requireV0WithITSTPC{"requireV0WithITSTPC", false, "flag to enforce V0s have ITS and TPC"}; + o2::framework::Configurable requireV0WithITSOnly{"requireV0WithITSOnly", false, "flag to select V0s with ITSonly tracks"}; + o2::framework::Configurable requireV0WithTPCOnly{"requireV0WithTPCOnly", false, "flag to select V0s with TPConly tracks"}; + o2::framework::Configurable minPtV0{"minPtV0", 0.1, "min pT for v0 photons at PV"}; + o2::framework::Configurable maxPtV0{"maxPtV0", 1e+10, "max pT for v0 photons at PV"}; + o2::framework::Configurable minEtaV0{"minEtaV0", -0.8, "min eta for v0 photons at PV"}; + o2::framework::Configurable maxEtaV0{"maxEtaV0", 0.8, "max eta for v0 photons at PV"}; + o2::framework::Configurable minRV0{"minRV0", 4.0, "min v0 radius"}; + o2::framework::Configurable maxRV0{"maxRV0", 90.0, "max v0 radius"}; + o2::framework::Configurable maxAlphaAP{"maxAlphaAP", 0.95, "max alpha for AP cut"}; + o2::framework::Configurable maxQtAP{"maxQtAP", 0.01, "max qT for AP cut"}; + o2::framework::Configurable minCosPA{"minCosPA", 0.999, "min V0 CosPA"}; + o2::framework::Configurable maxPCA{"maxPCA", 1.5, "max distance btween 2 legs"}; + o2::framework::Configurable maxChi2KF{"maxChi2KF", 1.e+10f, "max chi2/ndf with KF"}; + o2::framework::Configurable rejectV0onITSib{"rejectV0onITSib", true, "flag to reject V0s on ITSib"}; + o2::framework::Configurable applyPrefilter{"applyPrefilter", false, "flag to apply prefilter to V0"}; + + o2::framework::Configurable minNClusterTPC{"minNClusterTPC", 0, "min NCluster TPC"}; + o2::framework::Configurable minNCrossedRowsTPC{"minNCrossedRowsTPC", 40, "min ncrossed rows in TPC"}; + o2::framework::Configurable minNCrossedRowsOverFindableClustersTPC{"minNCrossedRowsOverFindableClustersTPC", 0.8f, "min fraction of crossed rows over findable clusters in TPC"}; + o2::framework::Configurable maxFracSharedClustersTPC{"maxFracSharedClustersTPC", 999.f, "max fraction of shared clusters in TPC"}; + o2::framework::Configurable minNClusterITS{"minNClusterITS", 0, "min NCluster ITS"}; + o2::framework::Configurable minMeanClusterSizeITSob{"minMeanClusterSizeITSob", 0.f, "min ITSob"}; + o2::framework::Configurable maxMeanClusterSizeITSob{"maxMeanClusterSizeITSob", 16.f, "max ITSob"}; + o2::framework::Configurable macChi2TPC{"macChi2TPC", 4.f, "max chi2/NclsTPC"}; + o2::framework::Configurable macChi2ITS{"macChi2ITS", 36.f, "max chi2/NclsITS"}; + o2::framework::Configurable minTPCNSigmaEl{"minTPCNSigmaEl", -3.0f, "min. TPC n sigma for electron"}; + o2::framework::Configurable maxTPCNSigmaEl{"maxTPCNSigmaEl", +3.0f, "max. TPC n sigma for electron"}; + o2::framework::Configurable disableITSOnly{"disableITSOnly", false, "flag to disable ITSonly tracks"}; + o2::framework::Configurable disableTPCOnly{"disableTPCOnly", false, "flag to disable TPConly tracks"}; + } pcmcuts; + struct : ConfigurableGroup { std::string prefix = "mesonConfig"; Configurable minOpenAngle{"minOpenAngle", 0.0202, "apply min opening angle. Default value one EMCal cell"}; @@ -216,15 +250,21 @@ struct TaskPi0FlowEMC { int runNow = 0; int runBefore = -1; - Filter clusterFilter = aod::skimmedcluster::time >= emccuts.cfgEMCminTime && aod::skimmedcluster::time <= emccuts.cfgEMCmaxTime && aod::skimmedcluster::m02 >= emccuts.cfgEMCminM02 && aod::skimmedcluster::m02 <= emccuts.cfgEMCmaxM02 && aod::skimmedcluster::e >= emccuts.cfgEMCminE; + // Filter clusterFilter = aod::skimmedcluster::time >= emccuts.cfgEMCminTime && aod::skimmedcluster::time <= emccuts.cfgEMCmaxTime && aod::skimmedcluster::m02 >= emccuts.cfgEMCminM02 && aod::skimmedcluster::m02 <= emccuts.cfgEMCmaxM02 && aod::skimmedcluster::e >= emccuts.cfgEMCminE; Filter collisionFilter = (nabs(aod::collision::posZ) <= eventcuts.cfgZvtxMax) && (aod::evsel::ft0cOccupancyInTimeRange <= eventcuts.cfgFT0COccupancyMax) && (aod::evsel::ft0cOccupancyInTimeRange >= eventcuts.cfgFT0COccupancyMin); - using FilteredEMCalPhotons = soa::Filtered>; - using EMCalPhotons = soa::Join; + // using FilteredEMCalPhotons = soa::Filtered>; + using EMCalPhotons = soa::Join; + using PCMPhotons = soa::Join; using FilteredCollsWithQvecs = soa::Filtered>; using CollsWithQvecs = soa::Join; using Colls = soa::Join; - Preslice perCollisionEMC = aod::emccluster::emeventId; + static constexpr std::size_t NQVecEntries = 6; + + PresliceOptional perCollisionEMC = o2::aod::emccluster::emeventId; + PresliceOptional perCollisionPCM = aod::v0photonkf::emeventId; + PresliceOptional perEMCClusterMT = o2::aod::mintm::minClusterId; + PresliceOptional perEMCClusterMS = o2::aod::mintm::minClusterId; HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; @@ -246,6 +286,8 @@ struct TaskPi0FlowEMC { // Usage when cfgEnableNonLin is enabled std::unique_ptr fEMCalCorrectionFactor; // ("fEMCalCorrectionFactor","(1 + [0]/x + [1]/x^2) / (1 + [2]/x)", 0.3, 100.); + float energyCorrectionFactor = 1.f; + uint8_t nColl = 1; // To access the 1D array inline int getIndex(int iEta, int iPhi) @@ -303,6 +345,37 @@ struct TaskPi0FlowEMC { fEMCCut.SetUseSecondaryTM(emccuts.emcUseSecondaryTM.value); // disables or enables secondary TM } + void definePCMCut() + { + fV0PhotonCut = V0PhotonCut("fV0PhotonCut", "fV0PhotonCut"); + + // for v0 + fV0PhotonCut.SetV0PtRange(pcmcuts.minPtV0, pcmcuts.maxPtV0); + fV0PhotonCut.SetV0EtaRange(pcmcuts.minEtaV0, pcmcuts.maxEtaV0); + fV0PhotonCut.SetMinCosPA(pcmcuts.minCosPA); + fV0PhotonCut.SetMaxPCA(pcmcuts.maxPCA); + fV0PhotonCut.SetMaxChi2KF(pcmcuts.maxChi2KF); + fV0PhotonCut.SetRxyRange(pcmcuts.minRV0, pcmcuts.maxRV0); + fV0PhotonCut.SetAPRange(pcmcuts.maxAlphaAP, pcmcuts.maxQtAP); + fV0PhotonCut.RejectITSib(pcmcuts.rejectV0onITSib); + + // for track + fV0PhotonCut.SetMinNClustersTPC(pcmcuts.minNClusterTPC); + fV0PhotonCut.SetMinNCrossedRowsTPC(pcmcuts.minNCrossedRowsTPC); + fV0PhotonCut.SetMinNCrossedRowsOverFindableClustersTPC(pcmcuts.minNCrossedRowsOverFindableClustersTPC); + fV0PhotonCut.SetMaxFracSharedClustersTPC(pcmcuts.maxFracSharedClustersTPC); + fV0PhotonCut.SetChi2PerClusterTPC(0.f, pcmcuts.macChi2TPC); + fV0PhotonCut.SetTPCNsigmaElRange(pcmcuts.minTPCNSigmaEl, pcmcuts.maxTPCNSigmaEl); + fV0PhotonCut.SetChi2PerClusterITS(0.f, pcmcuts.macChi2ITS); + fV0PhotonCut.SetNClustersITS(pcmcuts.minNClusterITS, 7); + fV0PhotonCut.SetMeanClusterSizeITSob(pcmcuts.minMeanClusterSizeITSob, pcmcuts.maxMeanClusterSizeITSob); + fV0PhotonCut.SetDisableITSonly(pcmcuts.disableITSOnly); + fV0PhotonCut.SetDisableTPConly(pcmcuts.disableTPCOnly); + fV0PhotonCut.SetRequireITSTPC(pcmcuts.requireV0WithITSTPC); + fV0PhotonCut.SetRequireITSonly(pcmcuts.requireV0WithITSOnly); + fV0PhotonCut.SetRequireTPConly(pcmcuts.requireV0WithTPCOnly); + } + void init(InitContext&) { if (harmonic != kElliptic && harmonic != kTriangluar) { @@ -311,29 +384,21 @@ struct TaskPi0FlowEMC { defineEMEventCut(); defineEMCCut(); + definePCMCut(); o2::aod::pwgem::photonmeson::utils::eventhistogram::addEventHistograms(®istry); const AxisSpec thnAxisInvMass{thnConfigAxisInvMass, "#it{M}_{#gamma#gamma} (GeV/#it{c}^{2})"}; const AxisSpec thnAxisPt{thnConfigAxisPt, "#it{p}_{T} (GeV/#it{c})"}; const AxisSpec thnAxisCent{thnConfigAxisCent, "Centrality (%)"}; - const AxisSpec thnAxisCosNPhi{thnConfigAxisCosNPhi, Form("cos(%d#varphi)", harmonic.value)}; const AxisSpec thnAxisCosDeltaPhi{thnConfigAxisCosDeltaPhi, Form("cos(%d(#varphi - #Psi_{sub}))", harmonic.value)}; const AxisSpec thnAxisM02{thnConfigAxisM02, "M_{02}"}; const AxisSpec thAxisTanThetaPhi{mesonConfig.thConfigAxisTanThetaPhi, "atan(#Delta#theta/#Delta#varphi)"}; const AxisSpec thAxisClusterEnergy{thnConfigAxisPt, "#it{E} (GeV)"}; const AxisSpec thAxisEnergyCalib{thnConfigAxisEnergyCalib, "#it{E}_{clus} (GeV)"}; const AxisSpec thAxisAlpha{100, -1., +1, "#alpha"}; - const AxisSpec thAxisMult{1000, 0., +1000, "#it{N}_{ch}"}; const AxisSpec thAxisEnergy{1000, 0., 100., "#it{E}_{clus} (GeV)"}; - const AxisSpec thAxisTime{1500, -600, 900, "#it{t}_{cl} (ns)"}; const AxisSpec thAxisEta{320, -0.8, 0.8, "#eta"}; const AxisSpec thAxisPhi{500, 0, 2 * 3.14159, "phi"}; - const AxisSpec thAxisNCell{17664, 0.5, +17664.5, "#it{N}_{cell}"}; - const AxisSpec thAxisPsi{360 / harmonic.value, -(1. / static_cast(harmonic.value)) * std::numbers::pi_v, (1. / static_cast(harmonic.value)) * std::numbers::pi_v, Form("#Psi_{%d}", harmonic.value)}; - const AxisSpec thAxisCN{8, 0.5, 8.5, "#it{c}_{n}"}; - const AxisSpec thAxisSN{8, 0.5, 8.5, "#it{s}_{n}"}; - const AxisSpec thAxisCPUTime{1000, 0, 10000, "#it{t} (#mus)"}; - const AxisSpec thAxisAzimuth{360, -o2::constants::math::PI, o2::constants::math::PI, "#it{#varphi} (rad)"}; const AxisSpec thnAxisMixingVtx{mixingConfig.cfgVtxBins, "#it{z} (cm)"}; const AxisSpec thnAxisMixingCent{mixingConfig.cfgCentBins, "Centrality (%)"}; @@ -354,84 +419,21 @@ struct TaskPi0FlowEMC { if (cfgDoPlaneQA.value) { registry.add("hSparsePi0FlowPlane", "THn for SP", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisCent, thnAxisCosDeltaPhi}); } - auto hClusterCuts = registry.add("hClusterCuts", "hClusterCuts;;Counts", kTH1D, {{6, 0.5, 6.5}}, false); - hClusterCuts->GetXaxis()->SetBinLabel(1, "in"); - hClusterCuts->GetXaxis()->SetBinLabel(2, "opening angle"); - hClusterCuts->GetXaxis()->SetBinLabel(3, "#it{M}_{#gamma#gamma}"); - hClusterCuts->GetXaxis()->SetBinLabel(4, "#it{p}_{T}"); - hClusterCuts->GetXaxis()->SetBinLabel(5, "conversion cut"); - hClusterCuts->GetXaxis()->SetBinLabel(6, "out"); - - auto hClusterCutsMixed = registry.add("hClusterCutsMixed", "hClusterCutsMixed;;Counts", kTH1D, {{6, 0.5, 6.5}}, false); - hClusterCutsMixed->GetXaxis()->SetBinLabel(1, "in"); - hClusterCutsMixed->GetXaxis()->SetBinLabel(2, "opening angle"); - hClusterCutsMixed->GetXaxis()->SetBinLabel(3, "#it{M}_{#gamma#gamma}"); - hClusterCutsMixed->GetXaxis()->SetBinLabel(4, "#it{p}_{T}"); - hClusterCutsMixed->GetXaxis()->SetBinLabel(5, "conversion cut"); - hClusterCutsMixed->GetXaxis()->SetBinLabel(6, "out"); - - if (saveSPResoHist.value) { - registry.add("spReso/hSpResoFT0cFT0a", "hSpResoFT0cFT0a; centrality; Q_{FT0c} #bullet Q_{FT0a}", HistType::kTProfile, {thnAxisCent}); - registry.add("spReso/hSpResoFT0cTPCpos", "hSpResoFT0cTPCpos; centrality; Q_{FT0c} #bullet Q_{TPCpos}", HistType::kTProfile, {thnAxisCent}); - registry.add("spReso/hSpResoFT0cTPCneg", "hSpResoFT0cTPCneg; centrality; Q_{FT0c} #bullet Q_{TPCneg}", HistType::kTProfile, {thnAxisCent}); - registry.add("spReso/hSpResoFT0cTPCtot", "hSpResoFT0cTPCtot; centrality; Q_{FT0c} #bullet Q_{TPCtot}", HistType::kTProfile, {thnAxisCent}); - registry.add("spReso/hSpResoFT0aTPCpos", "hSpResoFT0aTPCpos; centrality; Q_{FT0a} #bullet Q_{TPCpos}", HistType::kTProfile, {thnAxisCent}); - registry.add("spReso/hSpResoFT0aTPCneg", "hSpResoFT0aTPCneg; centrality; Q_{FT0a} #bullet Q_{TPCneg}", HistType::kTProfile, {thnAxisCent}); - registry.add("spReso/hSpResoFT0aTPCtot", "hSpResoFT0aTPCtot; centrality; Q_{FT0m} #bullet Q_{TPCtot}", HistType::kTProfile, {thnAxisCent}); - registry.add("spReso/hSpResoFT0mTPCpos", "hSpResoFT0mTPCpos; centrality; Q_{FT0m} #bullet Q_{TPCpos}", HistType::kTProfile, {thnAxisCent}); - registry.add("spReso/hSpResoFT0mTPCneg", "hSpResoFT0mTPCneg; centrality; Q_{FT0m} #bullet Q_{TPCneg}", HistType::kTProfile, {thnAxisCent}); - registry.add("spReso/hSpResoFT0mTPCtot", "hSpResoFT0mTPCtot; centrality; Q_{FT0m} #bullet Q_{TPCtot}", HistType::kTProfile, {thnAxisCent}); - registry.add("spReso/hSpResoTPCposTPCneg", "hSpResoTPCposTPCneg; centrality; Q_{TPCpos} #bullet Q_{TPCneg}", HistType::kTProfile, {thnAxisCent}); - registry.add("spReso/hSpResoFV0aFT0c", "hSpResoFV0aFT0c; centrality; Q_{FV0a} #bullet Q_{FT0c}", HistType::kTProfile, {thnAxisCent}); - registry.add("spReso/hSpResoFV0aTPCpos", "hSpResoFV0aTPCpos; centrality; Q_{FV0a} #bullet Q_{TPCpos}", HistType::kTProfile, {thnAxisCent}); - registry.add("spReso/hSpResoFV0aTPCneg", "hSpResoFV0aTPCneg; centrality; Q_{FV0a} #bullet Q_{TPCneg}", HistType::kTProfile, {thnAxisCent}); - registry.add("spReso/hSpResoFV0aTPCtot", "hSpResoFV0aTPCtot; centrality; Q_{FV0a} #bullet Q_{TPCtot}", HistType::kTProfile, {thnAxisCent}); - } - - if (saveEpResoHisto.value) { - registry.add("hEventPlaneAngleFT0C", "hEventPlaneAngleFT0C", HistType::kTH2D, {thnAxisCent, thAxisPsi}); - registry.add("hEventPlaneAngleFT0A", "hEventPlaneAngleFT0A", HistType::kTH2D, {thnAxisCent, thAxisPsi}); - registry.add("hEventPlaneAngleFT0M", "hEventPlaneAngleFT0M", HistType::kTH2D, {thnAxisCent, thAxisPsi}); - registry.add("hEventPlaneAngleFV0A", "hEventPlaneAngleFV0A", HistType::kTH2D, {thnAxisCent, thAxisPsi}); - registry.add("hEventPlaneAngleTPCpos", "hEventPlaneAngleTPCpos", HistType::kTH2D, {thnAxisCent, thAxisPsi}); - registry.add("hEventPlaneAngleTPCneg", "hEventPlaneAngleTPCneg", HistType::kTH2D, {thnAxisCent, thAxisPsi}); - registry.add("epReso/hEpResoFT0cFT0a", "hEpResoFT0cFT0a; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); - registry.add("epReso/hEpResoFT0cTPCpos", "hEpResoFT0cTPCpos; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); - registry.add("epReso/hEpResoFT0cTPCneg", "hEpResoFT0cTPCneg; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); - registry.add("epReso/hEpResoFT0cTPCtot", "hEpResoFT0cTPCtot; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); - registry.add("epReso/hEpResoFT0aTPCpos", "hEpResoFT0aTPCpos; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); - registry.add("epReso/hEpResoFT0aTPCneg", "hEpResoFT0aTPCneg; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); - registry.add("epReso/hEpResoFT0aTPCtot", "hEpResoFT0aTPCtot; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); - registry.add("epReso/hEpResoFT0mTPCpos", "hEpResoFT0mTPCpos; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); - registry.add("epReso/hEpResoFT0mTPCneg", "hEpResoFT0mTPCneg; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); - registry.add("epReso/hEpResoFT0mTPCtot", "hEpResoFT0mTPCtot; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); - registry.add("epReso/hEpResoTPCposTPCneg", "hEpResoTPCposTPCneg; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); - registry.add("epReso/hEpResoFV0aFT0c", "hEpResoFV0aFT0c; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); - registry.add("epReso/hEpResoFV0aTPCpos", "hEpResoFV0aTPCpos; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); - registry.add("epReso/hEpResoFV0aTPCneg", "hEpResoFV0aTPCneg; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); - registry.add("epReso/hEpResoFV0aTPCtot", "hEpResoFV0aTPCtot; centrality; #Delta#Psi_{sub}", HistType::kTH2D, {thnAxisCent, thnAxisCosNPhi}); - registry.add("epReso/hEpCosCoefficientsFT0c", "hEpCosCoefficientsFT0c; centrality; c_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisCN}); - registry.add("epReso/hEpSinCoefficientsFT0c", "hEpSinCoefficientsFT0c; centrality; s_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisSN}); - registry.add("epReso/hEpCosCoefficientsFT0a", "hEpCosCoefficientsFT0a; centrality; c_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisCN}); - registry.add("epReso/hEpSinCoefficientsFT0a", "hEpSinCoefficientsFT0a; centrality; s_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisSN}); - registry.add("epReso/hEpCosCoefficientsFV0a", "hEpCosCoefficientsFV0a; centrality; c_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisCN}); - registry.add("epReso/hEpSinCoefficientsFV0a", "hEpSinCoefficientsFV0a; centrality; s_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisSN}); - registry.add("epReso/hEpCosCoefficientsFT0m", "hEpCosCoefficientsFT0m; centrality; c_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisCN}); - registry.add("epReso/hEpSinCoefficientsFT0m", "hEpSinCoefficientsFT0m; centrality; s_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisSN}); - registry.add("epReso/hEpCosCoefficientsTPCpos", "hEpCosCoefficientsTPCpos; centrality; c_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisCN}); - registry.add("epReso/hEpSinCoefficientsTPCpos", "hEpSinCoefficientsTPCpos; centrality; s_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisSN}); - registry.add("epReso/hEpCosCoefficientsTPCneg", "hEpCosCoefficientsTPCneg; centrality; c_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisCN}); - registry.add("epReso/hEpSinCoefficientsTPCneg", "hEpSinCoefficientsTPCneg; centrality; s_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisSN}); - registry.add("epReso/hEpCosCoefficientsTPCTots", "hEpCosCoefficientsTPCTots; centrality; c_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisCN}); - registry.add("epReso/hEpSinCoefficientsTPCTots", "hEpSinCoefficientsTPCTots; centrality; s_{n}", HistType::kTProfile2D, {thnAxisCent, thAxisSN}); - registry.add("QVector/hQVecMeanRVsPhiFT0a", "hQVecMeanRVsPhiFT0a; centrality; #it{#varphi} (rad), <#it{r}> (a.u.)", HistType::kTProfile2D, {thnAxisCent, thAxisAzimuth}); - registry.add("QVector/hQVecMeanRVsPhiFT0c", "hQVecMeanRVsPhiFT0c; centrality; #it{#varphi} (rad), <#it{r}> (a.u.)", HistType::kTProfile2D, {thnAxisCent, thAxisAzimuth}); - registry.add("QVector/hQVecMeanRVsPhiFT0m", "hQVecMeanRVsPhiFT0m; centrality; #it{#varphi} (rad), <#it{r}> (a.u.)", HistType::kTProfile2D, {thnAxisCent, thAxisAzimuth}); - registry.add("QVector/hQVecMeanRVsPhiFV0a", "hQVecMeanRVsPhiFV0a; centrality; #it{#varphi} (rad), <#it{r}> (a.u.)", HistType::kTProfile2D, {thnAxisCent, thAxisAzimuth}); - registry.add("QVector/hQVecMeanRVsPhiTPCpos", "hQVecMeanRVsPhiTPCpos; centrality; #it{#varphi} (rad), <#it{r}> (a.u.)", HistType::kTProfile2D, {thnAxisCent, thAxisAzimuth}); - registry.add("QVector/hQVecMeanRVsPhiTPCneg", "hQVecMeanRVsPhiTPCneg; centrality; #it{#varphi} (rad), <#it{r}> (a.u.)", HistType::kTProfile2D, {thnAxisCent, thAxisAzimuth}); - registry.add("QVector/hQVecMeanRVsPhiTPCTot", "hQVecMeanRVsPhiTPCTot; centrality; #it{#varphi} (rad), <#it{r}> (a.u.)", HistType::kTProfile2D, {thnAxisCent, thAxisAzimuth}); - } + auto hMesonCuts = registry.add("hMesonCuts", "hMesonCuts;;Counts", kTH1D, {{6, 0.5, 6.5}}, false); + hMesonCuts->GetXaxis()->SetBinLabel(1, "in"); + hMesonCuts->GetXaxis()->SetBinLabel(2, "opening angle"); + hMesonCuts->GetXaxis()->SetBinLabel(3, "#it{M}_{#gamma#gamma}"); + hMesonCuts->GetXaxis()->SetBinLabel(4, "#it{p}_{T}"); + hMesonCuts->GetXaxis()->SetBinLabel(5, "conversion cut"); + hMesonCuts->GetXaxis()->SetBinLabel(6, "out"); + + auto hMesonCutsMixed = registry.add("hMesonCutsMixed", "hMesonCutsMixed;;Counts", kTH1D, {{6, 0.5, 6.5}}, false); + hMesonCutsMixed->GetXaxis()->SetBinLabel(1, "in"); + hMesonCutsMixed->GetXaxis()->SetBinLabel(2, "opening angle"); + hMesonCutsMixed->GetXaxis()->SetBinLabel(3, "#it{M}_{#gamma#gamma}"); + hMesonCutsMixed->GetXaxis()->SetBinLabel(4, "#it{p}_{T}"); + hMesonCutsMixed->GetXaxis()->SetBinLabel(5, "conversion cut"); + hMesonCutsMixed->GetXaxis()->SetBinLabel(6, "out"); if (emccuts.cfgEnableQA.value) { registry.add("clusterQA/hEClusterBefore", "Histo for cluster energy before cuts", HistType::kTH1D, {thAxisClusterEnergy}); @@ -478,9 +480,9 @@ struct TaskPi0FlowEMC { /// Change radians to degree /// \param angle in radians /// \return angle in degree - float getAngleDegree(float angle) + constexpr float getAngleDegree(float angle) { - return angle * 180.f * std::numbers::inv_pi_v; + return angle * o2::constants::math::Rad2Deg; } /// Compute the delta psi in the range [0, pi/harmonic] @@ -508,7 +510,7 @@ struct TaskPi0FlowEMC { /// Get the centrality /// \param collision is the collision with the centrality information - template + template float getCentrality(TCollision const& collision) { float cent = -999.; @@ -532,8 +534,8 @@ struct TaskPi0FlowEMC { /// Get all used Q vector /// \param collision is the collision with the Q vector information - template - std::vector getAllQvec(TCollision const& collision) + template + std::array getAllQvec(TCollision const& collision) { // Retrieve the Q vectors using the helper function for each detector auto [xQVecMain, yQVecMain] = getQvec(collision, qvecDetector); @@ -545,7 +547,7 @@ struct TaskPi0FlowEMC { /// Get the Q vector /// \param collision is the collision with the Q vector information - template + template std::pair getQvec(TCollision const& collision, int detector) { float xQVec = -999.f; @@ -631,7 +633,7 @@ struct TaskPi0FlowEMC { /// Check if the QVector values are within reasonable range /// \param collision is the collision with the Q vector information - bool isQvecGood(std::vector const& QVecs) + bool isQvecGood(std::array const& QVecs) { bool isgood = true; for (const auto& QVec : QVecs) { @@ -688,7 +690,7 @@ struct TaskPi0FlowEMC { return masked; } - template + template void initCCDB(TCollision const& collision) { // Load EMCal geometry @@ -731,7 +733,7 @@ struct TaskPi0FlowEMC { } /// \brief Calculate background using rotation background method - template + template void rotationBackground(const ROOT::Math::PtEtaPhiMVector& meson, ROOT::Math::PtEtaPhiMVector photon1, ROOT::Math::PtEtaPhiMVector photon2, TPhotons const& photons_coll, unsigned int ig1, unsigned int ig2, TCollision const& collision) { // if less than 3 clusters are present skip event since we need at least 3 clusters @@ -838,7 +840,7 @@ struct TaskPi0FlowEMC { } /// \brief Calculate background using rotation background method for calib - template + template void rotationBackgroundCalib(const ROOT::Math::PtEtaPhiMVector& meson, ROOT::Math::PtEtaPhiMVector photon1, ROOT::Math::PtEtaPhiMVector photon2, TPhotons const& photons_coll, unsigned int ig1, unsigned int ig2, TCollision const& collision) { // if less than 3 clusters are present skip event since we need at least 3 clusters @@ -923,7 +925,7 @@ struct TaskPi0FlowEMC { /// Compute the scalar product /// \param collision is the collision with the Q vector information and event plane /// \param meson are the selected candidates - template + template void runFlowAnalysis(TCollision const& collision, ROOT::Math::PtEtaPhiMVector const& meson) { auto [xQVec, yQVec] = getQvec(collision, qvecDetector); @@ -951,33 +953,114 @@ struct TaskPi0FlowEMC { return; } - // Pi0 from EMCal - void processEMCal(CollsWithQvecs const& collisions, EMCalPhotons const& clusters) + /// \brief check if standard event cuts + FT0 occupancy + centrality + QVec good is + /// \param collision collision that will be checked + /// \return true if collision survives all checks, otherwise false + template + bool isFullEventSelected(TCollision const& collision, bool fillHisto = false) { - int nColl = 1; - float energyCorrectionFactor = 1.f; - if (cfgDoReverseScaling.value) { - energyCorrectionFactor = 1.0505f; + if (fillHisto) { + o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<0>(®istry, collision); } - for (const auto& collision : collisions) { - auto photonsPerCollision = clusters.sliceBy(perCollisionEMC, collision.globalIndex()); + if (!(fEMEventCut.IsSelected(collision))) { + // general event selection + return false; + } + if (!(eventcuts.cfgFT0COccupancyMin <= collision.ft0cOccupancyInTimeRange() && collision.ft0cOccupancyInTimeRange() < eventcuts.cfgFT0COccupancyMax)) { + // occupancy selection + return false; + } + float cent = getCentrality(collision); + if (cent < eventcuts.cfgMinCent || cent > eventcuts.cfgMaxCent) { + // event selection + return false; + } + if (!isQvecGood(getAllQvec(collision))) { + // selection based on QVector + return false; + } + if (fillHisto) { + o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<1>(®istry, collision); + registry.fill(HIST("Event/before/hCollisionCounter"), 12.0); // accepted + registry.fill(HIST("Event/after/hCollisionCounter"), 12.0); // accepted + } + return true; + } - o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<0>(®istry, collision); - if (!(fEMEventCut.IsSelected(collision))) { - // general event selection + template + void runPairingLoop(TCollision const& collision, TPhotons1 const& photons1, TPhotons2 const& photons2, EMBitFlags const& flags1, EMBitFlags const& flags2) + { + for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photons1, photons2))) { + if (!(flags1.test(g1.globalIndex())) || !(flags2.test(g2.globalIndex()))) { continue; } - if (!(eventcuts.cfgFT0COccupancyMin <= collision.ft0cOccupancyInTimeRange() && collision.ft0cOccupancyInTimeRange() < eventcuts.cfgFT0COccupancyMax)) { - // occupancy selection + + // Cut edge clusters away, similar to rotation method to ensure same acceptance is used + if (cfgDistanceToEdge.value) { + if (checkEtaPhi1D(g1.eta(), RecoDecay::constrainAngle(g1.phi())) >= cfgEMCalMapLevelSameEvent.value) { + continue; + } + if (checkEtaPhi1D(g2.eta(), RecoDecay::constrainAngle(g2.phi())) >= cfgEMCalMapLevelSameEvent.value) { + continue; + } + } + if (correctionConfig.cfgEnableNonLin.value) { + energyCorrectionFactor = fEMCalCorrectionFactor->Eval(g1.e() > MinEnergy ? g1.e() : MinEnergy); + } + if (correctionConfig.cfgEnableNonLin.value) { + energyCorrectionFactor = fEMCalCorrectionFactor->Eval(g2.e() > MinEnergy ? g2.e() : MinEnergy); + } + + ROOT::Math::PtEtaPhiMVector v1(energyCorrectionFactor * g1.pt(), g1.eta(), g1.phi(), 0.); + ROOT::Math::PtEtaPhiMVector v2(energyCorrectionFactor * g2.pt(), g2.eta(), g2.phi(), 0.); + ROOT::Math::PtEtaPhiMVector vMeson = v1 + v2; + + float dTheta = v1.Theta() - v2.Theta(); + float dPhi = v1.Phi() - v2.Phi(); + float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P())); + + registry.fill(HIST("hMesonCuts"), 1); + if (openingAngle <= mesonConfig.minOpenAngle) { + registry.fill(HIST("hMesonCuts"), 2); continue; } - float cent = getCentrality(collision); - if (cent < eventcuts.cfgMinCent || cent > eventcuts.cfgMaxCent) { - // event selection + if (cfgDoRotation && nColl % cfgDownsampling.value == 0) { + rotationBackground(vMeson, v1, v2, photons1, g1.globalIndex(), g2.globalIndex(), collision); + } + if (thnConfigAxisInvMass.value[1] > vMeson.M() || thnConfigAxisInvMass.value.back() < vMeson.M()) { + registry.fill(HIST("hMesonCuts"), 3); continue; } - if (!isQvecGood(getAllQvec(collision))) { - // selection based on QVector + if (thnConfigAxisPt.value[1] > vMeson.Pt() || thnConfigAxisPt.value.back() < vMeson.Pt()) { + registry.fill(HIST("hMesonCuts"), 4); + continue; + } + if (mesonConfig.cfgEnableQA.value) { + registry.fill(HIST("mesonQA/hInvMassPt"), vMeson.M(), vMeson.Pt()); + registry.fill(HIST("mesonQA/hTanThetaPhi"), vMeson.M(), getAngleDegree(std::atan(dTheta / dPhi))); + registry.fill(HIST("mesonQA/hAlphaPt"), (v1.E() - v2.E()) / (v1.E() + v2.E()), vMeson.Pt()); + } + if (mesonConfig.enableTanThetadPhi && mesonConfig.minTanThetadPhi > std::fabs(getAngleDegree(std::atan(dTheta / dPhi)))) { + registry.fill(HIST("hMesonCuts"), 5); + continue; + } + registry.fill(HIST("hMesonCuts"), 6); + runFlowAnalysis<0>(collision, vMeson); + } + } + + // Pi0 from EMCal + void processEMCal(CollsWithQvecs const& collisions, EMCalPhotons const& clusters, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds) + { + EMBitFlags flags(clusters.size()); + fEMCCut.AreSelectedRunning(flags, clusters, matchedPrims, matchedSeconds); + + if (cfgDoReverseScaling.value) { + energyCorrectionFactor = 1.0505f; + } + for (const auto& collision : collisions) { + + if (!isFullEventSelected(collision)) { continue; } runNow = collision.runNumber(); @@ -985,9 +1068,8 @@ struct TaskPi0FlowEMC { initCCDB(collision); runBefore = runNow; } - o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<1>(®istry, collision); - registry.fill(HIST("Event/before/hCollisionCounter"), 12.0); // accepted - registry.fill(HIST("Event/after/hCollisionCounter"), 12.0); // accepted + + auto photonsPerCollision = clusters.sliceBy(perCollisionEMC, collision.globalIndex()); if (emccuts.cfgEnableQA.value) { for (const auto& photon : photonsPerCollision) { @@ -1003,62 +1085,7 @@ struct TaskPi0FlowEMC { registry.fill(HIST("clusterQA/hClusterEtaPhiAfter"), photon.phi(), photon.eta()); // after cuts } } - for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photonsPerCollision, photonsPerCollision))) { - if (!(fEMCCut.IsSelected(g1)) || !(fEMCCut.IsSelected(g2))) { - continue; - } - - // Cut edge clusters away, similar to rotation method to ensure same acceptance is used - if (cfgDistanceToEdge.value) { - if (checkEtaPhi1D(g1.eta(), RecoDecay::constrainAngle(g1.phi())) >= cfgEMCalMapLevelSameEvent.value) { - continue; - } - if (checkEtaPhi1D(g2.eta(), RecoDecay::constrainAngle(g2.phi())) >= cfgEMCalMapLevelSameEvent.value) { - continue; - } - } - if (correctionConfig.cfgEnableNonLin.value) { - energyCorrectionFactor = fEMCalCorrectionFactor->Eval(g1.e() > MinEnergy ? g1.e() : MinEnergy); - } - ROOT::Math::PtEtaPhiMVector v1(energyCorrectionFactor * g1.pt(), g1.eta(), g1.phi(), 0.); - if (correctionConfig.cfgEnableNonLin.value) { - energyCorrectionFactor = fEMCalCorrectionFactor->Eval(g2.e() > MinEnergy ? g2.e() : MinEnergy); - } - ROOT::Math::PtEtaPhiMVector v2(energyCorrectionFactor * g2.pt(), g2.eta(), g2.phi(), 0.); - ROOT::Math::PtEtaPhiMVector vMeson = v1 + v2; - float dTheta = v1.Theta() - v2.Theta(); - float dPhi = v1.Phi() - v2.Phi(); - float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P())); - registry.fill(HIST("hClusterCuts"), 1); - if (openingAngle <= mesonConfig.minOpenAngle) { - registry.fill(HIST("hClusterCuts"), 2); - continue; - } - if (cfgDoRotation) { - if (nColl % cfgDownsampling.value == 0) { - rotationBackground(vMeson, v1, v2, photonsPerCollision, g1.globalIndex(), g2.globalIndex(), collision); - } - } - if (thnConfigAxisInvMass.value[1] > vMeson.M() || thnConfigAxisInvMass.value.back() < vMeson.M()) { - registry.fill(HIST("hClusterCuts"), 3); - continue; - } - if (thnConfigAxisPt.value[1] > vMeson.Pt() || thnConfigAxisPt.value.back() < vMeson.Pt()) { - registry.fill(HIST("hClusterCuts"), 4); - continue; - } - if (mesonConfig.cfgEnableQA.value) { - registry.fill(HIST("mesonQA/hInvMassPt"), vMeson.M(), vMeson.Pt()); - registry.fill(HIST("mesonQA/hTanThetaPhi"), vMeson.M(), getAngleDegree(std::atan(dTheta / dPhi))); - registry.fill(HIST("mesonQA/hAlphaPt"), (v1.E() - v2.E()) / (v1.E() + v2.E()), vMeson.Pt()); - } - if (mesonConfig.enableTanThetadPhi && mesonConfig.minTanThetadPhi > std::fabs(getAngleDegree(std::atan(dTheta / dPhi)))) { - registry.fill(HIST("hClusterCuts"), 5); - continue; - } - registry.fill(HIST("hClusterCuts"), 6); - runFlowAnalysis<0>(collision, vMeson); - } + runPairingLoop(collision, photonsPerCollision, photonsPerCollision, flags, flags); if (cfgDoRotation) { if (nColl % cfgDownsampling.value == 0) { nColl = 1; // reset counter @@ -1071,23 +1098,20 @@ struct TaskPi0FlowEMC { PROCESS_SWITCH(TaskPi0FlowEMC, processEMCal, "Process EMCal Pi0 candidates", true); // Pi0 from EMCal - void processEMCalMixed(FilteredCollsWithQvecs const& collisions, FilteredEMCalPhotons const& clusters) + void processEMCalMixed(FilteredCollsWithQvecs const& collisions, EMCalPhotons const& clusters, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds) { + EMBitFlags flags(clusters.size()); + fEMCCut.AreSelectedRunning(flags, clusters, matchedPrims, matchedSeconds); float energyCorrectionFactor = 1.f; if (cfgDoReverseScaling.value) { energyCorrectionFactor = 1.0505f; } - auto getClustersSize = - [&clusters, this](FilteredCollsWithQvecs::iterator const& col) { - auto associatedClusters = clusters.sliceByCached(emccluster::emeventId, col.globalIndex(), this->cache); // it's cached, so slicing/grouping happens only once - return associatedClusters.size(); - }; - using BinningType = FlexibleBinningPolicy, aod::collision::PosZ, aod::cent::CentFT0C, emevent::EP2FT0M>; - BinningType binningWithLambda{{getClustersSize}, {mixingConfig.cfgVtxBins, mixingConfig.cfgCentBins, mixingConfig.cfgEPBins}, true}; + using BinningType = ColumnBinningPolicy>; + BinningType binningMixedEvent{{mixingConfig.cfgVtxBins, mixingConfig.cfgCentBins, mixingConfig.cfgEPBins}, true}; auto clustersTuple = std::make_tuple(clusters); - SameKindPair pair{binningWithLambda, mixingConfig.cfgMixingDepth, -1, collisions, clustersTuple, &cache}; // indicates that 5 events should be mixed and under/overflow (-1) to be ignored + SameKindPair pair{binningMixedEvent, mixingConfig.cfgMixingDepth, -1, collisions, clustersTuple, &cache}; // indicates that 5 events should be mixed and under/overflow (-1) to be ignored for (const auto& [c1, clusters1, c2, clusters2] : pair) { if (!(fEMEventCut.IsSelected(c1)) || !(fEMEventCut.IsSelected(c2))) { @@ -1113,7 +1137,7 @@ struct TaskPi0FlowEMC { } registry.fill(HIST("h3DMixingCount"), c1.posZ(), getCentrality(c1), c1.ep2ft0m()); for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(clusters1, clusters2))) { - if (!(fEMCCut.IsSelected(g1)) || !(fEMCCut.IsSelected(g2))) { + if (!(flags.test(g1.globalIndex())) || !(flags.test(g2.globalIndex()))) { continue; } // Cut edge clusters away, similar to rotation method to ensure same acceptance is used @@ -1139,17 +1163,17 @@ struct TaskPi0FlowEMC { float dPhi = v1.Phi() - v2.Phi(); float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P())); - registry.fill(HIST("hClusterCutsMixed"), 1); + registry.fill(HIST("hMesonCutsMixed"), 1); if (openingAngle <= mesonConfig.minOpenAngle) { - registry.fill(HIST("hClusterCutsMixed"), 2); + registry.fill(HIST("hMesonCutsMixed"), 2); continue; } if (thnConfigAxisInvMass.value[1] > vMeson.M() || thnConfigAxisInvMass.value.back() < vMeson.M()) { - registry.fill(HIST("hClusterCutsMixed"), 3); + registry.fill(HIST("hMesonCutsMixed"), 3); continue; } if (thnConfigAxisPt.value[1] > vMeson.Pt() || thnConfigAxisPt.value.back() < vMeson.Pt()) { - registry.fill(HIST("hClusterCutsMixed"), 4); + registry.fill(HIST("hMesonCutsMixed"), 4); continue; } if (mesonConfig.cfgEnableQA.value) { @@ -1158,186 +1182,28 @@ struct TaskPi0FlowEMC { registry.fill(HIST("mesonQA/hAlphaPtMixed"), (v1.E() - v2.E()) / (v1.E() + v2.E()), vMeson.Pt()); } if (mesonConfig.enableTanThetadPhi && mesonConfig.minTanThetadPhi > std::fabs(getAngleDegree(std::atan(dTheta / dPhi)))) { - registry.fill(HIST("hClusterCutsMixed"), 5); + registry.fill(HIST("hMesonCutsMixed"), 5); continue; } - registry.fill(HIST("hClusterCutsMixed"), 6); + registry.fill(HIST("hMesonCutsMixed"), 6); runFlowAnalysis<2>(c1, vMeson); } } } PROCESS_SWITCH(TaskPi0FlowEMC, processEMCalMixed, "Process EMCal Pi0 mixed event candidates", false); - // Resolution - void processResolution(CollsWithQvecs::iterator const& collision) + // PCM-EMCal same event + void processEMCalPCMC(CollsWithQvecs const& collisions, EMCalPhotons const& clusters, PCMPhotons const& photons, aod::V0Legs const&, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds) { - o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<0>(®istry, collision); - if (!(fEMEventCut.IsSelected(collision))) { - // no selection on the centrality is applied on purpose to allow for the resolution study in post-processing - return; - } - if (!(eventcuts.cfgFT0COccupancyMin <= collision.ft0cOccupancyInTimeRange() && collision.ft0cOccupancyInTimeRange() < eventcuts.cfgFT0COccupancyMax)) { - return; - } - float cent = getCentrality(collision); - if (cent < eventcuts.cfgMinCent || cent > eventcuts.cfgMaxCent) { - // event selection - return; - } - if (!isQvecGood(getAllQvec(collision))) { - // selection based on QVector - return; - } - o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<1>(®istry, collision); - registry.fill(HIST("Event/before/hCollisionCounter"), 12.0); // accepted - registry.fill(HIST("Event/after/hCollisionCounter"), 12.0); // accepted - - float centrality = getCentrality(collision); // centrality not updated in the rejection mask function - float xQVecFT0a = -999.f; - float yQVecFT0a = -999.f; - float xQVecFT0c = -999.f; - float yQVecFT0c = -999.f; - float xQVecFT0m = -999.f; - float yQVecFT0m = -999.f; - float xQVecBPos = -999.f; - float yQVecBPos = -999.f; - float xQVecBNeg = -999.f; - float yQVecBNeg = -999.f; - float xQVecBTot = -999.f; - float yQVecBTot = -999.f; - float xQVecFV0a = -999.f; - float yQVecFV0a = -999.f; - if (harmonic == kElliptic) { - xQVecFT0a = collision.q2xft0a(); - yQVecFT0a = collision.q2yft0a(); - xQVecFT0c = collision.q2xft0c(); - yQVecFT0c = collision.q2yft0c(); - xQVecFT0m = collision.q2xft0m(); - yQVecFT0m = collision.q2yft0m(); - xQVecBPos = collision.q2xbpos(); - yQVecBPos = collision.q2ybpos(); - xQVecBNeg = collision.q2xbneg(); - yQVecBNeg = collision.q2ybneg(); - xQVecBTot = collision.q2xbtot(); - yQVecBTot = collision.q2ybtot(); - xQVecFV0a = collision.q2xfv0a(); - yQVecFV0a = collision.q2yfv0a(); - } else if (harmonic == kTriangluar) { - xQVecFT0a = collision.q3xft0a(); - yQVecFT0a = collision.q3yft0a(); - xQVecFT0c = collision.q3xft0c(); - yQVecFT0c = collision.q3yft0c(); - xQVecFT0m = collision.q3xft0m(); - yQVecFT0m = collision.q3yft0m(); - xQVecBPos = collision.q3xbpos(); - yQVecBPos = collision.q3ybpos(); - xQVecBNeg = collision.q3xbneg(); - yQVecBNeg = collision.q3ybneg(); - xQVecBTot = collision.q3xbtot(); - yQVecBTot = collision.q3ybtot(); - xQVecFV0a = collision.q3xfv0a(); - yQVecFV0a = collision.q3yfv0a(); - } - - if (saveSPResoHist) { - registry.fill(HIST("spReso/hSpResoFT0cFT0a"), centrality, xQVecFT0c * xQVecFT0a + yQVecFT0c * yQVecFT0a); - registry.fill(HIST("spReso/hSpResoFT0cTPCpos"), centrality, xQVecFT0c * xQVecBPos + yQVecFT0c * yQVecBPos); - registry.fill(HIST("spReso/hSpResoFT0cTPCneg"), centrality, xQVecFT0c * xQVecBNeg + yQVecFT0c * yQVecBNeg); - registry.fill(HIST("spReso/hSpResoFT0cTPCtot"), centrality, xQVecFT0c * xQVecBTot + yQVecFT0c * yQVecBTot); - registry.fill(HIST("spReso/hSpResoFT0aTPCpos"), centrality, xQVecFT0a * xQVecBPos + yQVecFT0a * yQVecBPos); - registry.fill(HIST("spReso/hSpResoFT0aTPCneg"), centrality, xQVecFT0a * xQVecBNeg + yQVecFT0a * yQVecBNeg); - registry.fill(HIST("spReso/hSpResoFT0aTPCtot"), centrality, xQVecFT0a * xQVecBTot + yQVecFT0a * yQVecBTot); - registry.fill(HIST("spReso/hSpResoFT0mTPCpos"), centrality, xQVecFT0m * xQVecBPos + yQVecFT0m * yQVecBPos); - registry.fill(HIST("spReso/hSpResoFT0mTPCneg"), centrality, xQVecFT0m * xQVecBNeg + yQVecFT0m * yQVecBNeg); - registry.fill(HIST("spReso/hSpResoFT0mTPCtot"), centrality, xQVecFT0m * xQVecBTot + yQVecFT0m * yQVecBTot); - registry.fill(HIST("spReso/hSpResoTPCposTPCneg"), centrality, xQVecBPos * xQVecBNeg + yQVecBPos * yQVecBNeg); - registry.fill(HIST("spReso/hSpResoFV0aFT0c"), centrality, xQVecFV0a * xQVecFT0c + yQVecFV0a * yQVecFT0c); - registry.fill(HIST("spReso/hSpResoFV0aTPCpos"), centrality, xQVecFV0a * xQVecBPos + yQVecFV0a * yQVecBPos); - registry.fill(HIST("spReso/hSpResoFV0aTPCneg"), centrality, xQVecFV0a * xQVecBNeg + yQVecFV0a * yQVecBNeg); - registry.fill(HIST("spReso/hSpResoFV0aTPCtot"), centrality, xQVecFV0a * xQVecBTot + yQVecFV0a * yQVecBTot); - } - - if (saveEpResoHisto) { - float epFT0a = epHelper.GetEventPlane(xQVecFT0a, yQVecFT0a, harmonic); - float epFT0c = epHelper.GetEventPlane(xQVecFT0c, yQVecFT0c, harmonic); - float epFT0m = epHelper.GetEventPlane(xQVecFT0m, yQVecFT0m, harmonic); - float epBPoss = epHelper.GetEventPlane(xQVecBPos, yQVecBPos, harmonic); - float epBNegs = epHelper.GetEventPlane(xQVecBNeg, yQVecBNeg, harmonic); - float epBTots = epHelper.GetEventPlane(xQVecBTot, yQVecBTot, harmonic); - float epFV0a = epHelper.GetEventPlane(xQVecFV0a, yQVecFV0a, harmonic); - - registry.fill(HIST("hEventPlaneAngleFT0M"), centrality, epFT0m); - registry.fill(HIST("hEventPlaneAngleFT0A"), centrality, epFT0c); - registry.fill(HIST("hEventPlaneAngleFT0C"), centrality, epFT0a); - registry.fill(HIST("hEventPlaneAngleTPCpos"), centrality, epBPoss); - registry.fill(HIST("hEventPlaneAngleTPCneg"), centrality, epBNegs); - registry.fill(HIST("hEventPlaneAngleFV0A"), centrality, epFV0a); - - registry.fill(HIST("QVector/hQVecMeanRVsPhiFT0a"), centrality, std::atan2(yQVecFT0a, xQVecFT0a), std::hypot(xQVecFT0a, yQVecFT0a)); - registry.fill(HIST("QVector/hQVecMeanRVsPhiFT0c"), centrality, std::atan2(yQVecFT0c, xQVecFT0c), std::hypot(xQVecFT0c, yQVecFT0c)); - registry.fill(HIST("QVector/hQVecMeanRVsPhiFT0m"), centrality, std::atan2(yQVecFT0m, xQVecFT0m), std::hypot(xQVecFT0m, yQVecFT0m)); - registry.fill(HIST("QVector/hQVecMeanRVsPhiTPCpos"), centrality, std::atan2(yQVecBPos, xQVecBPos), std::hypot(xQVecBPos, yQVecBPos)); - registry.fill(HIST("QVector/hQVecMeanRVsPhiTPCneg"), centrality, std::atan2(yQVecBNeg, xQVecBNeg), std::hypot(xQVecBNeg, yQVecBNeg)); - registry.fill(HIST("QVector/hQVecMeanRVsPhiTPCTot"), centrality, std::atan2(yQVecBTot, xQVecBTot), std::hypot(xQVecBTot, yQVecBTot)); - registry.fill(HIST("QVector/hQVecMeanRVsPhiFV0a"), centrality, std::atan2(yQVecFV0a, xQVecFV0a), std::hypot(xQVecFV0a, yQVecFV0a)); - - registry.fill(HIST("epReso/hEpResoFT0cFT0a"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0c, epFT0a))); - registry.fill(HIST("epReso/hEpResoFT0cTPCpos"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0c, epBPoss))); - registry.fill(HIST("epReso/hEpResoFT0cTPCneg"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0c, epBNegs))); - registry.fill(HIST("epReso/hEpResoFT0cTPCtot"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0c, epBTots))); - registry.fill(HIST("epReso/hEpResoFT0aTPCpos"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0a, epBPoss))); - registry.fill(HIST("epReso/hEpResoFT0aTPCneg"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0a, epBNegs))); - registry.fill(HIST("epReso/hEpResoFT0aTPCtot"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0a, epBTots))); - registry.fill(HIST("epReso/hEpResoFT0mTPCpos"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0m, epBPoss))); - registry.fill(HIST("epReso/hEpResoFT0mTPCneg"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0m, epBNegs))); - registry.fill(HIST("epReso/hEpResoFT0mTPCtot"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0m, epBTots))); - registry.fill(HIST("epReso/hEpResoTPCposTPCneg"), centrality, std::cos(harmonic * getDeltaPsiInRange(epBPoss, epBNegs))); - registry.fill(HIST("epReso/hEpResoFV0aFT0c"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFV0a, epFT0c))); - registry.fill(HIST("epReso/hEpResoFV0aTPCpos"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFV0a, epBPoss))); - registry.fill(HIST("epReso/hEpResoFV0aTPCneg"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFV0a, epBNegs))); - registry.fill(HIST("epReso/hEpResoFV0aTPCtot"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFV0a, epBTots))); - for (int n = 1; n <= kOctagonal; n++) { - registry.fill(HIST("epReso/hEpCosCoefficientsFT0c"), centrality, n, std::cos(n * harmonic * epFT0c)); - registry.fill(HIST("epReso/hEpSinCoefficientsFT0c"), centrality, n, std::sin(n * harmonic * epFT0c)); - registry.fill(HIST("epReso/hEpCosCoefficientsFT0a"), centrality, n, std::cos(n * harmonic * epFT0a)); - registry.fill(HIST("epReso/hEpSinCoefficientsFT0a"), centrality, n, std::sin(n * harmonic * epFT0a)); - registry.fill(HIST("epReso/hEpCosCoefficientsFT0m"), centrality, n, std::cos(n * harmonic * epFT0m)); - registry.fill(HIST("epReso/hEpSinCoefficientsFT0m"), centrality, n, std::sin(n * harmonic * epFT0m)); - registry.fill(HIST("epReso/hEpCosCoefficientsTPCpos"), centrality, n, std::cos(n * harmonic * epBPoss)); - registry.fill(HIST("epReso/hEpSinCoefficientsTPCpos"), centrality, n, std::sin(n * harmonic * epBPoss)); - registry.fill(HIST("epReso/hEpCosCoefficientsTPCneg"), centrality, n, std::cos(n * harmonic * epBNegs)); - registry.fill(HIST("epReso/hEpSinCoefficientsTPCneg"), centrality, n, std::sin(n * harmonic * epBNegs)); - registry.fill(HIST("epReso/hEpCosCoefficientsTPCTots"), centrality, n, std::cos(n * harmonic * epBTots)); - registry.fill(HIST("epReso/hEpSinCoefficientsTPCTots"), centrality, n, std::sin(n * harmonic * epBTots)); - registry.fill(HIST("epReso/hEpCosCoefficientsFV0a"), centrality, n, std::cos(n * harmonic * epFV0a)); - registry.fill(HIST("epReso/hEpSinCoefficientsFV0a"), centrality, n, std::sin(n * harmonic * epFV0a)); - } - } - } - PROCESS_SWITCH(TaskPi0FlowEMC, processResolution, "Process resolution", false); + EMBitFlags emcFlags(clusters.size()); + fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds); - // EMCal calibration - void processEMCalCalib(Colls const& collisions, EMCalPhotons const& clusters) - { - float energyCorrectionFactor = 1.f; - if (!correctionConfig.doEMCalCalib) { - return; + if (cfgDoReverseScaling.value) { + energyCorrectionFactor = 1.0505f; } - int nColl = 1; for (const auto& collision : collisions) { - auto photonsPerCollision = clusters.sliceBy(perCollisionEMC, collision.globalIndex()); - o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<0>(®istry, collision); - if (!(fEMEventCut.IsSelected(collision))) { - // general event selection - continue; - } - if (!(eventcuts.cfgFT0COccupancyMin <= collision.ft0cOccupancyInTimeRange() && collision.ft0cOccupancyInTimeRange() < eventcuts.cfgFT0COccupancyMax)) { - // occupancy selection - continue; - } - float cent = getCentrality(collision); - if (cent < eventcuts.cfgMinCent || cent > eventcuts.cfgMaxCent) { - // event selection + + if (!isFullEventSelected(collision)) { continue; } runNow = collision.runNumber(); @@ -1345,21 +1211,21 @@ struct TaskPi0FlowEMC { initCCDB(collision); runBefore = runNow; } - o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<1>(®istry, collision); - registry.fill(HIST("Event/before/hCollisionCounter"), 12.0); // accepted - registry.fill(HIST("Event/after/hCollisionCounter"), 12.0); // accepted - if (emccuts.cfgEnableQA.value) { - for (const auto& photon : photonsPerCollision) { + auto photonsEMCPerCollision = clusters.sliceBy(perCollisionEMC, collision.globalIndex()); + auto photonsPCMPerCollision = photons.sliceBy(perCollisionPCM, collision.globalIndex()); + + if (emccuts.cfgEnableQA) { + for (const auto& photon : photonsEMCPerCollision) { registry.fill(HIST("clusterQA/hEClusterBefore"), photon.e()); // before cuts - if (!(fEMCCut.IsSelected(photon))) { + if (!(fEMCCut.IsSelected(photon))) { continue; } registry.fill(HIST("clusterQA/hEClusterAfter"), photon.e()); // accepted after cuts } } - for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photonsPerCollision, photonsPerCollision))) { - if (!(fEMCCut.IsSelected(g1)) || !(fEMCCut.IsSelected(g2))) { + for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(photonsEMCPerCollision, photonsPCMPerCollision))) { + if (!(emcFlags.test(g1.globalIndex())) || !(fV0PhotonCut.IsSelected(g2))) { continue; } @@ -1368,73 +1234,129 @@ struct TaskPi0FlowEMC { if (checkEtaPhi1D(g1.eta(), RecoDecay::constrainAngle(g1.phi())) >= cfgEMCalMapLevelSameEvent.value) { continue; } - if (checkEtaPhi1D(g2.eta(), RecoDecay::constrainAngle(g2.phi())) >= cfgEMCalMapLevelSameEvent.value) { - continue; - } } - if (correctionConfig.cfgEnableNonLin.value) { - energyCorrectionFactor = fEMCalCorrectionFactor->Eval(g1.e() > MinEnergy ? g1.e() : MinEnergy); - } + // EMCal photon v1 + energyCorrectionFactor = fEMCalCorrectionFactor->Eval(g1.e() > MinEnergy ? g1.e() : MinEnergy); ROOT::Math::PtEtaPhiMVector v1(energyCorrectionFactor * g1.pt(), g1.eta(), g1.phi(), 0.); - if (correctionConfig.cfgEnableNonLin.value) { - energyCorrectionFactor = fEMCalCorrectionFactor->Eval(g2.e() > MinEnergy ? g2.e() : MinEnergy); - } - ROOT::Math::PtEtaPhiMVector v2(energyCorrectionFactor * g2.pt(), g2.eta(), g2.phi(), 0.); + // PCM photon v2s + ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); ROOT::Math::PtEtaPhiMVector vMeson = v1 + v2; float dTheta = v1.Theta() - v2.Theta(); float dPhi = v1.Phi() - v2.Phi(); float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P())); - registry.fill(HIST("hClusterCuts"), 1); + registry.fill(HIST("hMesonCuts"), 1); if (openingAngle <= mesonConfig.minOpenAngle) { - registry.fill(HIST("hClusterCuts"), 2); + registry.fill(HIST("hMesonCuts"), 2); continue; } - if (cfgDoRotation) { - if (nColl % cfgDownsampling.value == 0) { - rotationBackgroundCalib(vMeson, v1, v2, photonsPerCollision, g1.globalIndex(), g2.globalIndex(), collision); - } - } if (thnConfigAxisInvMass.value[1] > vMeson.M() || thnConfigAxisInvMass.value.back() < vMeson.M()) { - registry.fill(HIST("hClusterCuts"), 3); + registry.fill(HIST("hMesonCuts"), 3); continue; } if (thnConfigAxisPt.value[1] > vMeson.Pt() || thnConfigAxisPt.value.back() < vMeson.Pt()) { - registry.fill(HIST("hClusterCuts"), 4); + registry.fill(HIST("hMesonCuts"), 4); continue; } - if (mesonConfig.cfgEnableQA.value) { + if (mesonConfig.cfgEnableQA) { registry.fill(HIST("mesonQA/hInvMassPt"), vMeson.M(), vMeson.Pt()); registry.fill(HIST("mesonQA/hTanThetaPhi"), vMeson.M(), getAngleDegree(std::atan(dTheta / dPhi))); registry.fill(HIST("mesonQA/hAlphaPt"), (v1.E() - v2.E()) / (v1.E() + v2.E()), vMeson.Pt()); } if (mesonConfig.enableTanThetadPhi && mesonConfig.minTanThetadPhi > std::fabs(getAngleDegree(std::atan(dTheta / dPhi)))) { - registry.fill(HIST("hClusterCuts"), 5); + registry.fill(HIST("hMesonCuts"), 5); continue; } - if (std::fabs((v1.E() - v2.E()) / (v1.E() + v2.E()) < cfgMaxAsymmetry)) { // only use symmetric decays - registry.fill(HIST("hClusterCuts"), 6); - registry.fill(HIST("hSparseCalibSE"), vMeson.M(), vMeson.E() / 2., getCentrality(collision)); - } + runFlowAnalysis<0>(collision, vMeson); } - if (cfgDoRotation) { - if (nColl % cfgDownsampling.value == 0) { - nColl = 1; // reset counter - } else { - nColl++; + } + } + PROCESS_SWITCH(TaskPi0FlowEMC, processEMCalPCMC, "Process neutral meson flow using PCM-EMC same event", false); + + // PCM-EMCal mixed event + void processEMCalPCMMixed(CollsWithQvecs const& collisions, EMCalPhotons const& clusters, PCMPhotons const& pcmPhotons, aod::V0Legs const&, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds) + { + using BinningTypeMixed = ColumnBinningPolicy>; + BinningTypeMixed binningOnPositions{{mixingConfig.cfgVtxBins, mixingConfig.cfgCentBins, mixingConfig.cfgEPBins}, true}; + + auto associatedTables = std::make_tuple(clusters, pcmPhotons); + + Pair pairPCMEMC{binningOnPositions, mixingConfig.cfgMixingDepth, -1, collisions, associatedTables, &cache}; // indicates that mixingConfig.cfgMixingDepth events should be mixed and under/overflow (-1) to be ignored + + float energyCorrectionFactor = 1.f; + EMBitFlags emcFlags(clusters.size()); + fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds); + for (const auto& [c1, photonEMC, c2, photonPCM] : pairPCMEMC) { + if (!(fEMEventCut.IsSelected(c1)) || !(fEMEventCut.IsSelected(c2))) { + // general event selection + continue; + } + if (!(eventcuts.cfgFT0COccupancyMin <= c1.ft0cOccupancyInTimeRange() && c1.ft0cOccupancyInTimeRange() < eventcuts.cfgFT0COccupancyMax) || !(eventcuts.cfgFT0COccupancyMin <= c2.ft0cOccupancyInTimeRange() && c2.ft0cOccupancyInTimeRange() < eventcuts.cfgFT0COccupancyMax)) { + // occupancy selection + continue; + } + if (getCentrality(c1) < eventcuts.cfgMinCent || getCentrality(c1) > eventcuts.cfgMaxCent || getCentrality(c2) < eventcuts.cfgMinCent || getCentrality(c2) > eventcuts.cfgMaxCent) { + // event selection + continue; + } + runNow = c1.runNumber(); + if (runNow != runBefore) { + initCCDB(c1); + runBefore = runNow; + } + for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(photonEMC, photonPCM))) { + if (!(emcFlags.test(g1.globalIndex())) || !(fV0PhotonCut.IsSelected(g2))) { + continue; + } + // Cut edge clusters away, similar to rotation method to ensure same acceptance is used + if (cfgDistanceToEdge.value) { + if (checkEtaPhi1D(g1.eta(), RecoDecay::constrainAngle(g1.phi())) >= cfgEMCalMapLevelBackground.value) { + continue; + } + } + energyCorrectionFactor = fEMCalCorrectionFactor->Eval(g1.e() > MinEnergy ? g1.e() : MinEnergy); + ROOT::Math::PtEtaPhiMVector v1(energyCorrectionFactor * g1.pt(), g1.eta(), g1.phi(), 0.); + ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); + ROOT::Math::PtEtaPhiMVector vMeson = v1 + v2; + + float dTheta = v1.Theta() - v2.Theta(); + float dPhi = v1.Phi() - v2.Phi(); + float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P())); + + registry.fill(HIST("hMesonCutsMixed"), 1); + if (openingAngle <= mesonConfig.minOpenAngle) { + registry.fill(HIST("hMesonCutsMixed"), 2); + continue; + } + if (thnConfigAxisInvMass.value[1] > vMeson.M() || thnConfigAxisInvMass.value.back() < vMeson.M()) { + registry.fill(HIST("hMesonCutsMixed"), 3); + continue; + } + if (thnConfigAxisPt.value[1] > vMeson.Pt() || thnConfigAxisPt.value.back() < vMeson.Pt()) { + registry.fill(HIST("hMesonCutsMixed"), 4); + continue; + } + if (mesonConfig.cfgEnableQA) { + registry.fill(HIST("mesonQA/hInvMassPtMixed"), vMeson.M(), vMeson.Pt()); + registry.fill(HIST("mesonQA/hTanThetaPhiMixed"), vMeson.M(), getAngleDegree(std::atan(dTheta / dPhi))); + registry.fill(HIST("mesonQA/hAlphaPtMixed"), (v1.E() - v2.E()) / (v1.E() + v2.E()), vMeson.Pt()); } + if (mesonConfig.enableTanThetadPhi && mesonConfig.minTanThetadPhi > std::fabs(getAngleDegree(std::atan(dTheta / dPhi)))) { + registry.fill(HIST("hMesonCutsMixed"), 5); + continue; + } + registry.fill(HIST("hMesonCutsMixed"), 6); + runFlowAnalysis<2>(c1, vMeson); } } } - PROCESS_SWITCH(TaskPi0FlowEMC, processEMCalCalib, "Process EMCal calibration mixed event", false); + PROCESS_SWITCH(TaskPi0FlowEMC, processEMCalPCMMixed, "Process neutral meson flow using PCM-EMC mixed event", false); // Pi0 from EMCal - void processM02(CollsWithQvecs const& collisions, EMCalPhotons const& clusters) + void processM02(CollsWithQvecs const& collisions, EMCalPhotons const& clusters, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds) { for (const auto& collision : collisions) { - auto photonsPerCollision = clusters.sliceBy(perCollisionEMC, collision.globalIndex()); - o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<0>(®istry, collision); if (!(fEMEventCut.IsSelected(collision))) { // general event selection @@ -1462,12 +1384,16 @@ struct TaskPi0FlowEMC { registry.fill(HIST("Event/before/hCollisionCounter"), 12.0); // accepted registry.fill(HIST("Event/after/hCollisionCounter"), 12.0); // accepted + auto photonsPerCollision = clusters.sliceBy(perCollisionEMC, collision.globalIndex()); + for (const auto& photon : photonsPerCollision) { if (emccuts.cfgEnableQA.value) { registry.fill(HIST("clusterQA/hEClusterBefore"), photon.e()); // before cuts registry.fill(HIST("clusterQA/hClusterEtaPhiBefore"), photon.phi(), photon.eta()); // before cuts } - if (!(fEMCCut.IsSelected(photon))) { + auto matchedPrimsPerCluster = matchedPrims.sliceBy(perEMCClusterMT, photon.globalIndex()); + auto matchedSecondsPerCluster = matchedSeconds.sliceBy(perEMCClusterMS, photon.globalIndex()); + if (!(fEMCCut.IsSelected(photon, matchedPrimsPerCluster, matchedSecondsPerCluster))) { continue; } if (cfgDistanceToEdge.value && (checkEtaPhi1D(photon.eta(), RecoDecay::constrainAngle(photon.phi())) >= cfgEMCalMapLevelSameEvent.value)) {