diff --git a/PWGEM/PhotonMeson/Core/DalitzEECut.cxx b/PWGEM/PhotonMeson/Core/DalitzEECut.cxx index df6418474a1..f129a1ba289 100644 --- a/PWGEM/PhotonMeson/Core/DalitzEECut.cxx +++ b/PWGEM/PhotonMeson/Core/DalitzEECut.cxx @@ -9,14 +9,18 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -// -// Class for dilepton Cut -// +/// \file EMCPhotonCut.cxx +/// \brief header of class for dalitz ee cuts. +/// \author D. Sekihata, daiki.sekihata@cern.ch #include "PWGEM/PhotonMeson/Core/DalitzEECut.h" #include "Framework/Logger.h" +#include + +#include +#include #include #include diff --git a/PWGEM/PhotonMeson/Core/DalitzEECut.h b/PWGEM/PhotonMeson/Core/DalitzEECut.h index ddc1da18340..c238b5bead3 100644 --- a/PWGEM/PhotonMeson/Core/DalitzEECut.h +++ b/PWGEM/PhotonMeson/Core/DalitzEECut.h @@ -9,9 +9,9 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -// -// Class for dalitz ee selection -// +/// \file EMCPhotonCut.cxx +/// \brief header of class for dalitz ee cuts. +/// \author D. Sekihata, daiki.sekihata@cern.ch #ifndef PWGEM_PHOTONMESON_CORE_DALITZEECUT_H_ #define PWGEM_PHOTONMESON_CORE_DALITZEECUT_H_ @@ -19,21 +19,19 @@ #include "PWGEM/Dilepton/Utils/EMTrackUtilities.h" #include "PWGEM/Dilepton/Utils/PairUtilities.h" -#include "Tools/ML/MlResponse.h" -#include "Tools/ML/model.h" +#include -#include "CommonConstants/PhysicsConstants.h" -#include "Framework/DataTypes.h" -#include "Framework/Logger.h" +#include // IWYU pragma: keep +#include +#include -#include "Math/Vector4D.h" -#include "TNamed.h" +#include #include +#include +#include #include -#include #include -#include using namespace o2::aod::pwgem::dilepton::utils::emtrackutil; @@ -74,13 +72,9 @@ class DalitzEECut : public TNamed kTPConly = 1, }; - template - bool IsSelected(TPair const& pair) const + template + bool IsSelected(TTrack1 const& t1, TTrack2 const& t2, float bz) const { - auto t1 = std::get<0>(pair); - auto t2 = std::get<1>(pair); - float bz = std::get<2>(pair); - if (!IsSelectedTrack(t1) || !IsSelectedTrack(t2)) { return false; } @@ -361,7 +355,7 @@ class DalitzEECut : public TNamed float mMinTPCNsigmaPi{0}, mMaxTPCNsigmaPi{0}; float mMinTOFNsigmaEl{-1e+10}, mMaxTOFNsigmaEl{+1e+10}; - ClassDef(DalitzEECut, 1); + ClassDef(DalitzEECut, 2); }; #endif // PWGEM_PHOTONMESON_CORE_DALITZEECUT_H_ diff --git a/PWGEM/PhotonMeson/Core/DiphotonHadronMPC.h b/PWGEM/PhotonMeson/Core/DiphotonHadronMPC.h index 49d405e0dd1..54381491765 100644 --- a/PWGEM/PhotonMeson/Core/DiphotonHadronMPC.h +++ b/PWGEM/PhotonMeson/Core/DiphotonHadronMPC.h @@ -8,12 +8,9 @@ // 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 DiphotonHadronMPC.h /// \brief This code is to analyze diphoton-hadron correlation. Keep in mind that cumulant method does not require event mixing. -/// /// \author D. Sekihata, daiki.sekihata@cern.ch #ifndef PWGEM_PHOTONMESON_CORE_DIPHOTONHADRONMPC_H_ @@ -21,37 +18,46 @@ #include "PWGEM/Dilepton/Core/EMTrackCut.h" #include "PWGEM/Dilepton/Utils/EMTrack.h" -#include "PWGEM/Dilepton/Utils/EMTrackUtilities.h" #include "PWGEM/Dilepton/Utils/EventMixingHandler.h" #include "PWGEM/PhotonMeson/Core/DalitzEECut.h" #include "PWGEM/PhotonMeson/Core/EMPhotonEventCut.h" #include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" #include "PWGEM/PhotonMeson/Utils/EventHistograms.h" -#include "PWGEM/PhotonMeson/Utils/NMHistograms.h" #include "PWGEM/PhotonMeson/Utils/PairUtilities.h" -#include "Common/CCDB/RCTSelectionFlags.h" -#include "Common/Core/RecoDecay.h" - -#include "CCDB/BasicCCDBManager.h" -#include "DataFormatsParameters/GRPMagField.h" -#include "DataFormatsParameters/GRPObject.h" -#include "DetectorsBase/GeometryManager.h" -#include "Framework/ASoAHelpers.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" - -#include "Math/Vector4D.h" -#include "TString.h" +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include // IWYU pragma: keep +#include #include -#include -#include -#include +#include +#include #include #include +#include #include #include #include @@ -353,11 +359,11 @@ struct DiphotonHadronMPC { // diphoton-hadron info const AxisSpec axis_deta{ConfDEtaBins, deta_axis_title}; - const AxisSpec axis_dphi{cfgNbinsDPhi, -M_PI / 2, +3 * M_PI / 2, dphi_axis_title}; + const AxisSpec axis_dphi{cfgNbinsDPhi, -o2::constants::math::PIHalf, +3 * o2::constants::math::PIHalf, dphi_axis_title}; const AxisSpec axis_pt_hadron{ConfPtHadronBins, "p_{T,h} (GeV/c)"}; const AxisSpec axis_eta_hadron{40, -2, +2, "#eta_{h}"}; - const AxisSpec axis_phi_hadron{36, 0, 2 * M_PI, "#varphi_{h} (rad.)"}; + const AxisSpec axis_phi_hadron{36, 0, o2::constants::math::TwoPI, "#varphi_{h} (rad.)"}; fRegistry.add("Hadron/hs", "hadron", kTHnSparseD, {axis_pt_hadron, axis_eta_hadron, axis_phi_hadron}, false); fRegistry.add("Hadron/hTrackBit", "track bit", kTH1D, {{65536, -0.5, 65535.5}}, false); @@ -370,7 +376,7 @@ struct DiphotonHadronMPC { // hadron-hadron const AxisSpec axis_deta_hh{60, -3, +3, "#Delta#eta = #eta_{h}^{ref1} - #eta_{h}^{ref2}"}; - const AxisSpec axis_dphi_hh{90, -M_PI / 2, +3 * M_PI / 2, "#Delta#varphi = #varphi_{h}^{ref1} - #varphi_{h}^{ref2} (rad.)"}; + const AxisSpec axis_dphi_hh{90, -o2::constants::math::PIHalf, +3 * o2::constants::math::PIHalf, "#Delta#varphi = #varphi_{h}^{ref1} - #varphi_{h}^{ref2} (rad.)"}; // const AxisSpec axis_cosndphi_hh{cfgNbinsCosNDPhi, -1, +1, std::format("cos({0:d}(#varphi_{{h}}^{{ref1}} - #varphi_{{h}}^{{ref2}}))", cfgNmod.value)}; fRegistry.add("HadronHadron/same/hDEtaDPhi", "hadron-hadron 2PC", kTH2D, {axis_dphi_hh, axis_deta_hh}, true); fRegistry.addClone("HadronHadron/same/", "HadronHadron/mix/"); @@ -564,7 +570,7 @@ struct DiphotonHadronMPC { auto photons2_per_collision = photons2.sliceBy(perCollision2, collision.globalIndex()); for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photons1_per_collision, photons2_per_collision))) { - if (!cut1.template IsSelected(g1) || !cut2.template IsSelected(g2)) { + if (!cut1.template IsSelected(g1) || !cut2.template IsSelected(g2)) { continue; } @@ -594,7 +600,7 @@ struct DiphotonHadronMPC { float deta = v12.Eta() - v3.Eta(); float dphi = v12.Phi() - v3.Phi(); // o2::math_utils::bringTo02Pi(dphi); - dphi = RecoDecay::constrainAngle(dphi, -M_PI / 2, 1U); + dphi = RecoDecay::constrainAngle(dphi, -o2::constants::math::PIHalf, 1U); fRegistry.fill(HIST("DiphotonHadron/same/hs"), v12.M(), v12.Pt(), deta, dphi); npair++; } @@ -624,7 +630,7 @@ struct DiphotonHadronMPC { auto electrons_per_collision = electrons->sliceByCached(o2::aod::emprimaryelectron::emeventId, collision.globalIndex(), cache); for (const auto& g1 : photons1_per_collision) { - if (!cut1.template IsSelected(g1)) { + if (!cut1.template IsSelected(g1)) { continue; } auto pos1 = g1.template posTrack_as(); @@ -670,7 +676,7 @@ struct DiphotonHadronMPC { float deta = veeg.Eta() - v3.Eta(); float dphi = veeg.Phi() - v3.Phi(); // o2::math_utils::bringTo02Pi(dphi); - dphi = RecoDecay::constrainAngle(dphi, -M_PI / 2, 1U); + dphi = RecoDecay::constrainAngle(dphi, -o2::constants::math::PIHalf, 1U); fRegistry.fill(HIST("DiphotonHadron/same/hs"), veeg.M(), veeg.Pt(), deta, dphi); npair++; @@ -720,7 +726,7 @@ struct DiphotonHadronMPC { float deta = ref1.eta() - ref2.eta(); float dphi = ref1.phi() - ref2.phi(); // o2::math_utils::bringTo02Pi(dphi); - dphi = RecoDecay::constrainAngle(dphi, -M_PI / 2, 1U); + dphi = RecoDecay::constrainAngle(dphi, -o2::constants::math::PIHalf, 1U); fRegistry.fill(HIST("HadronHadron/same/hDEtaDPhi"), dphi, deta); } } @@ -794,7 +800,7 @@ struct DiphotonHadronMPC { float deta = trg.eta() - ref.eta(); float dphi = trg.phi() - ref.phi(); // o2::math_utils::bringTo02Pi(dphi); - dphi = RecoDecay::constrainAngle(dphi, -M_PI / 2, 1U); + dphi = RecoDecay::constrainAngle(dphi, -o2::constants::math::PIHalf, 1U); fRegistry.fill(HIST("DiphotonHadron/mix/hs"), trg.mass(), trg.pt(), deta, dphi); } } @@ -889,7 +895,7 @@ struct DiphotonHadronMPC { float deta = trg.eta() - ref.eta(); float dphi = trg.phi() - ref.phi(); // o2::math_utils::bringTo02Pi(dphi); - dphi = RecoDecay::constrainAngle(dphi, -M_PI / 2, 1U); + dphi = RecoDecay::constrainAngle(dphi, -o2::constants::math::PIHalf, 1U); fRegistry.fill(HIST("DiphotonHadron/mix/hs"), trg.mass(), trg.pt(), deta, dphi); } } @@ -917,7 +923,7 @@ struct DiphotonHadronMPC { float deta = ref1.eta() - ref2.eta(); float dphi = ref1.phi() - ref2.phi(); // o2::math_utils::bringTo02Pi(dphi); - dphi = RecoDecay::constrainAngle(dphi, -M_PI / 2, 1U); + dphi = RecoDecay::constrainAngle(dphi, -o2::constants::math::PIHalf, 1U); fRegistry.fill(HIST("HadronHadron/mix/hDEtaDPhi"), dphi, deta); } } diff --git a/PWGEM/PhotonMeson/Core/EMCPhotonCut.h b/PWGEM/PhotonMeson/Core/EMCPhotonCut.h index e7becec48b5..7847253a4f0 100644 --- a/PWGEM/PhotonMeson/Core/EMCPhotonCut.h +++ b/PWGEM/PhotonMeson/Core/EMCPhotonCut.h @@ -16,12 +16,56 @@ #ifndef PWGEM_PHOTONMESON_CORE_EMCPHOTONCUT_H_ #define PWGEM_PHOTONMESON_CORE_EMCPHOTONCUT_H_ +#include "PWGEM/PhotonMeson/DataModel/gammaTables.h" + +#include + #include #include #include +#include +#include #include +#include + +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) { + // requires that the following are valid calls: + { cluster.deltaEta() } -> std::convertible_to>; + { cluster.deltaPhi() } -> std::convertible_to>; + { cluster.trackpt() } -> std::convertible_to>; + { cluster.trackp() } -> std::convertible_to>; +}; + +template +concept hasSecondaryMatching = requires(Cluster cluster) { + // requires that the following are valid calls: + { cluster.deltaEtaSec() } -> std::convertible_to>; + { cluster.deltaPhiSec() } -> std::convertible_to>; + { cluster.trackptSec() } -> std::convertible_to>; + { 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}; @@ -52,9 +96,11 @@ class EMCPhotonCut : public TNamed /// \brief check if given cluster survives all cuts /// \param cluster cluster to check + /// \param matchedTracks subtable of the matched primary tracks (optional) + /// \param matchedSecondaries subtable of the matched secondary tracks (optional) /// \return true if cluster survives all cuts else false - template - bool IsSelected(Cluster const& cluster) const + template + bool IsSelected(Cluster const& cluster, TMatchedTracks const& emcmatchedtracks = nullptr, TMatchedSecondaries const& secondaries = nullptr) const { if (!IsSelectedEMCal(EMCPhotonCuts::kDefinition, cluster)) { return false; @@ -71,10 +117,10 @@ class EMCPhotonCut : public TNamed if (!IsSelectedEMCal(EMCPhotonCuts::kTiming, cluster)) { return false; } - if (mUseTM && (!IsSelectedEMCal(EMCPhotonCuts::kTM, cluster))) { + if (mUseTM && (!IsSelectedEMCal(EMCPhotonCuts::kTM, cluster, emcmatchedtracks))) { return false; } - if (mUseSecondaryTM && (!IsSelectedEMCal(EMCPhotonCuts::kSecondaryTM, cluster))) { + if (mUseSecondaryTM && (!IsSelectedEMCal(EMCPhotonCuts::kSecondaryTM, cluster, secondaries))) { return false; } if (!IsSelectedEMCal(EMCPhotonCuts::kExotic, cluster)) { @@ -86,9 +132,11 @@ class EMCPhotonCut : public TNamed /// \brief check if given cluster survives a given cut /// \param cut enum of the cluster cut to check /// \param cluster cluster to check + /// \param matchedTracks subtable of the matched primary tracks (optional) + /// \param matchedSecondaries subtable of the matched secondary tracks (optional) /// \return true if cluster survives cut else false - template - bool IsSelectedEMCal(const EMCPhotonCuts& cut, Cluster const& cluster) const + template + bool IsSelectedEMCal(const EMCPhotonCuts& cut, Cluster const& cluster, TMatchedTracks const& emcmatchedtracks = nullptr) const { switch (cut) { case EMCPhotonCuts::kDefinition: @@ -107,34 +155,64 @@ class EMCPhotonCut : public TNamed return mMinTime <= cluster.time() && cluster.time() <= mMaxTime; case EMCPhotonCuts::kTM: { - auto dEtas = cluster.deltaEta(); // std:vector - auto dPhis = cluster.deltaPhi(); // std:vector - auto trackspt = cluster.trackpt(); // std:vector - auto tracksp = cluster.trackp(); // std:vector - int ntrack = tracksp.size(); - for (int itr = 0; itr < ntrack; itr++) { - float dEta = std::fabs(dEtas[itr]); - float dPhi = std::fabs(dPhis[itr]); - bool result = (dEta > GetTrackMatchingEta(trackspt[itr])) || (dPhi > GetTrackMatchingPhi(trackspt[itr])) || (cluster.e() / tracksp[itr] >= mMinEoverP); - if (!result) { - return false; + 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); + if (!result) { + return false; + } } + } else if constexpr (hasTrackMatching) { + auto dEtas = cluster.deltaEta(); // std:vector + auto dPhis = cluster.deltaPhi(); // std:vector + auto trackspt = cluster.trackpt(); // std:vector + auto tracksp = cluster.trackp(); // std:vector + int ntrack = tracksp.size(); + for (int itr = 0; itr < ntrack; itr++) { + float dEta = std::fabs(dEtas[itr]); + float dPhi = std::fabs(dPhis[itr]); + bool result = (dEta > GetTrackMatchingEta(trackspt[itr])) || (dPhi > GetTrackMatchingPhi(trackspt[itr])) || (cluster.e() / tracksp[itr] >= mMinEoverP); + if (!result) { + return false; + } + } + } else { + return true; } return true; // when we don't have any tracks the cluster should always survive the TM cut! } case EMCPhotonCuts::kSecondaryTM: { - auto dEtas = cluster.deltaEtaSec(); // std:vector - auto dPhis = cluster.deltaPhiSec(); // std:vector - auto trackspt = cluster.trackptSec(); // std:vector - auto tracksp = cluster.trackpSec(); // std:vector - int ntrack = tracksp.size(); - for (int itr = 0; itr < ntrack; itr++) { - float dEta = std::fabs(dEtas[itr]); - float dPhi = std::fabs(dPhis[itr]); - bool result = (dEta > GetSecTrackMatchingEta(trackspt[itr])) || (dPhi > GetSecTrackMatchingPhi(trackspt[itr])); - if (!result) { - return false; + 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); + if (!result) { + return false; + } + } + } else if constexpr (hasSecondaryMatching) { + auto dEtas = cluster.deltaEtaSec(); // std:vector + auto dPhis = cluster.deltaPhiSec(); // std:vector + auto trackspt = cluster.trackptSec(); // std:vector + auto tracksp = cluster.trackpSec(); // std:vector + int ntrack = tracksp.size(); + for (int itr = 0; itr < ntrack; itr++) { + float dEta = std::fabs(dEtas[itr]); + float dPhi = std::fabs(dPhis[itr]); + bool result = (dEta > GetSecTrackMatchingEta(trackspt[itr])) || (dPhi > GetSecTrackMatchingPhi(trackspt[itr])); + if (!result) { + return false; + } } + } else { + return true; } return true; // when we don't have any secondary tracks the cluster should always survive the secondary TM cut! } @@ -272,7 +350,7 @@ class EMCPhotonCut : public TNamed TrackMatchingParams mSecTrackMatchingEtaParams = {-1.f, 0.f, 0.f}; TrackMatchingParams mSecTrackMatchingPhiParams = {-1.f, 0.f, 0.f}; - ClassDef(EMCPhotonCut, 2); + ClassDef(EMCPhotonCut, 3); }; #endif // PWGEM_PHOTONMESON_CORE_EMCPHOTONCUT_H_ diff --git a/PWGEM/PhotonMeson/Core/PHOSPhotonCut.cxx b/PWGEM/PhotonMeson/Core/PHOSPhotonCut.cxx index 96db9e331b3..f2a2046b120 100644 --- a/PWGEM/PhotonMeson/Core/PHOSPhotonCut.cxx +++ b/PWGEM/PhotonMeson/Core/PHOSPhotonCut.cxx @@ -9,13 +9,16 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -// -// Class for track selection -// +/// \file PHOSPhotonCut.cxx +/// \brief Source of class for phos photon selection. +/// \author D. Sekihata, daiki.sekihata@cern.ch -#include "Framework/Logger.h" #include "PWGEM/PhotonMeson/Core/PHOSPhotonCut.h" +#include + +#include + ClassImp(PHOSPhotonCut); const char* PHOSPhotonCut::mCutNames[static_cast(PHOSPhotonCut::PHOSPhotonCuts::kNCuts)] = {"Energy", "Dispersion", "CPV"}; diff --git a/PWGEM/PhotonMeson/Core/PHOSPhotonCut.h b/PWGEM/PhotonMeson/Core/PHOSPhotonCut.h index c60944dd52b..ea7fd3440ab 100644 --- a/PWGEM/PhotonMeson/Core/PHOSPhotonCut.h +++ b/PWGEM/PhotonMeson/Core/PHOSPhotonCut.h @@ -9,21 +9,18 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -// -// Class for v0 photon selection -// +/// \file PHOSPhotonCut.h +/// \brief Header of class for phos photon selection. +/// \author D. Sekihata, daiki.sekihata@cern.ch #ifndef PWGEM_PHOTONMESON_CORE_PHOSPHOTONCUT_H_ #define PWGEM_PHOTONMESON_CORE_PHOSPHOTONCUT_H_ -#include -#include -#include -#include -#include "Framework/Logger.h" -#include "Framework/DataTypes.h" -#include "Rtypes.h" -#include "TNamed.h" +#include + +#include + +#include class PHOSPhotonCut : public TNamed { @@ -41,7 +38,7 @@ class PHOSPhotonCut : public TNamed static const char* mCutNames[static_cast(PHOSPhotonCuts::kNCuts)]; // Temporary function to check if track passes selection criteria. To be replaced by framework filters. - template + template bool IsSelected(Cluster const& cluster) const { // auto track = cluster.template MatchedTrack_as(); //please implement a column to point matched track index (DECLARE_SOA_ARRAY_INDEX_COLUMN) in SkimPHOSClusters table. @@ -74,7 +71,7 @@ class PHOSPhotonCut : public TNamed } // Temporary function to check if track passes a given selection criteria. To be replaced by framework filters. - template + template bool IsSelectedCluster(Cluster const& cls, const PHOSPhotonCuts& cut) const { switch (cut) { @@ -101,7 +98,7 @@ class PHOSPhotonCut : public TNamed private: float mMinEnergy{0.1f}, mMaxEnergy{1e+10f}; - ClassDef(PHOSPhotonCut, 1); + ClassDef(PHOSPhotonCut, 2); }; #endif // PWGEM_PHOTONMESON_CORE_PHOSPHOTONCUT_H_ diff --git a/PWGEM/PhotonMeson/Core/PhotonHBT.h b/PWGEM/PhotonMeson/Core/PhotonHBT.h index 2950c0c4ec6..08738ce8e07 100644 --- a/PWGEM/PhotonMeson/Core/PhotonHBT.h +++ b/PWGEM/PhotonMeson/Core/PhotonHBT.h @@ -8,45 +8,50 @@ // 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. -// -// ======================== -// -// This code loops over v0 photons and makes pairs for photon HBT analysis. -// Please write to: daiki.sekihata@cern.ch + +/// \file PhotonHBT.h +/// \brief This code loops over v0 photons and makes pairs for photon HBT analysis. +/// \author Daiki Sekihata, daiki.sekihata@cern.ch #ifndef PWGEM_PHOTONMESON_CORE_PHOTONHBT_H_ #define PWGEM_PHOTONMESON_CORE_PHOTONHBT_H_ #include "PWGEM/Dilepton/Utils/EMTrack.h" -#include "PWGEM/Dilepton/Utils/EMTrackUtilities.h" #include "PWGEM/Dilepton/Utils/EventMixingHandler.h" #include "PWGEM/PhotonMeson/Core/EMPhotonEventCut.h" #include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" #include "PWGEM/PhotonMeson/Utils/EventHistograms.h" - -// #include "CCDB/BasicCCDBManager.h" -// #include "DataFormatsParameters/GRPMagField.h" -// #include "DataFormatsParameters/GRPObject.h" -// #include "DetectorsBase/GeometryManager.h" - -#include "Framework/ASoAHelpers.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" -#include "MathUtils/Utils.h" - -#include "Math/GenVector/Boost.h" -#include "Math/Vector3D.h" -#include "Math/Vector4D.h" -#include "TString.h" +// +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include // IWYU pragma: keep +#include +#include // IWYU pragma: keep +#include +#include #include -#include -#include +#include +#include #include #include #include +#include #include #include #include @@ -63,7 +68,6 @@ using namespace o2::aod; using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::soa; -using namespace o2::aod::pwgem::dilepton::utils::emtrackutil; using namespace o2::aod::pwgem::dilepton::utils; using namespace o2::aod::pwgem::photon::core::photonhbt; @@ -95,7 +99,7 @@ struct PhotonHBT { Configurable ndiff_bc_mix{"ndiff_bc_mix", 594, "difference in global BC required in mixed events"}; ConfigurableAxis ConfVtxBins{"ConfVtxBins", {VARIABLE_WIDTH, -10.0f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; ConfigurableAxis ConfCentBins{"ConfCentBins", {VARIABLE_WIDTH, 0.0f, 5.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.f, 999.f}, "Mixing bins - centrality"}; - ConfigurableAxis ConfEPBins{"ConfEPBins", {16, -M_PI / 2, +M_PI / 2}, "Mixing bins - event plane angle"}; + ConfigurableAxis ConfEPBins{"ConfEPBins", {16, -o2::constants::math::PIHalf, +o2::constants::math::PIHalf}, "Mixing bins - event plane angle"}; ConfigurableAxis ConfOccupancyBins{"ConfOccupancyBins", {VARIABLE_WIDTH, -1, 1e+10}, "Mixing bins - occupancy"}; Configurable cfgUseLCMS{"cfgUseLCMS", true, "measure relative momentum in LCMS for 1D"}; // always in LCMS for 3D @@ -329,8 +333,8 @@ struct PhotonHBT { { // o2::aod::pwgem::dilepton::utils::eventhistogram::addEventHistograms<-1>(&fRegistry); static constexpr std::string_view qvec_det_names[6] = {"FT0M", "FT0A", "FT0C", "BTot", "BPos", "BNeg"}; - fRegistry.add("Event/before/hEP2_CentFT0C_forMix", Form("2nd harmonics event plane for mix;centrality FT0C (%%);#Psi_{2}^{%s} (rad.)", qvec_det_names[cfgEP2Estimator_for_Mix].data()), kTH2F, {{110, 0, 110}, {180, -M_PI_2, +M_PI_2}}, false); - fRegistry.add("Event/after/hEP2_CentFT0C_forMix", Form("2nd harmonics event plane for mix;centrality FT0C (%%);#Psi_{2}^{%s} (rad.)", qvec_det_names[cfgEP2Estimator_for_Mix].data()), kTH2F, {{110, 0, 110}, {180, -M_PI_2, +M_PI_2}}, false); + fRegistry.add("Event/before/hEP2_CentFT0C_forMix", Form("2nd harmonics event plane for mix;centrality FT0C (%%);#Psi_{2}^{%s} (rad.)", qvec_det_names[cfgEP2Estimator_for_Mix].data()), kTH2F, {{110, 0, 110}, {180, -o2::constants::math::PIHalf, +o2::constants::math::PIHalf}}, false); + fRegistry.add("Event/after/hEP2_CentFT0C_forMix", Form("2nd harmonics event plane for mix;centrality FT0C (%%);#Psi_{2}^{%s} (rad.)", qvec_det_names[cfgEP2Estimator_for_Mix].data()), kTH2F, {{110, 0, 110}, {180, -o2::constants::math::PIHalf, +o2::constants::math::PIHalf}}, false); // pair info const AxisSpec axis_kt{ConfKtBins, "k_{T} (GeV/c)"}; @@ -543,7 +547,7 @@ struct PhotonHBT { auto photons1_coll = photons1.sliceBy(perCollision1, collision.globalIndex()); auto photons2_coll = photons2.sliceBy(perCollision2, collision.globalIndex()); for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photons1_coll, photons2_coll))) { - if (!cut1.template IsSelected(g1) || !cut2.template IsSelected(g2)) { + if (!cut1.template IsSelected(g1) || !cut2.template IsSelected(g2)) { continue; } @@ -563,8 +567,8 @@ struct PhotonHBT { ROOT::Math::XYZVector cp2(g2.vx(), g2.vy(), g2.vz()); float opa = std::acos(cp1.Dot(cp2) / (std::sqrt(cp1.Mag2()) * std::sqrt(cp2.Mag2()))); // opening angle between 2 conversion points o2::math_utils::bringTo02Pi(opa); - if (opa > M_PI) { - opa -= M_PI; + if (opa > o2::constants::math::PI) { + opa -= o2::constants::math::PI; } float cosOA = std::cos(opa / 2.f); if (dr / cosOA < ggpaircuts.cfgMinDR_CosOA) { @@ -634,8 +638,8 @@ struct PhotonHBT { ROOT::Math::XYZVector cp2(g2.vx(), g2.vy(), g2.vz()); float opa = std::acos(cp1.Dot(cp2) / (std::sqrt(cp1.Mag2()) * std::sqrt(cp2.Mag2()))); // opening angle between 2 conversion points o2::math_utils::bringTo02Pi(opa); - if (opa > M_PI) { - opa -= M_PI; + if (opa > o2::constants::math::PI) { + opa -= o2::constants::math::PI; } float cosOA = std::cos(opa / 2.f); if (dr / cosOA < ggpaircuts.cfgMinDR_CosOA) { diff --git a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h index 7b845f5399c..b13ff0b8695 100644 --- a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h +++ b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h @@ -59,6 +59,7 @@ #include #include +#include #include #include #include @@ -67,7 +68,7 @@ #include #include -template +template struct Pi0EtaToGammaGamma { o2::framework::Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; o2::framework::Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; @@ -192,6 +193,8 @@ struct Pi0EtaToGammaGamma { o2::framework::Configurable EMC_Eoverp{"EMC_Eoverp", 1.75, "Minimum cluster energy over track momentum for EMCal track matching"}; o2::framework::Configurable EMC_UseExoticCut{"EMC_UseExoticCut", true, "FLag to use the EMCal exotic cluster cut"}; o2::framework::Configurable cfgDistanceToEdge{"cfgDistanceToEdge", 1, "Distance to edge in cells required for rotated cluster to be accepted"}; + o2::framework::Configurable emcUseTM{"emcUseTM", true, "flag to use EMCal track matching cut or not."}; + o2::framework::Configurable emcUseSecondaryTM{"emcUseSecondaryTM", false, "flag to use EMCal secondary track matching cut or not. Required for PCM-EMC analyses."}; } emccuts; PHOSPhotonCut fPHOSCut; @@ -204,6 +207,28 @@ struct Pi0EtaToGammaGamma { // static constexpr std::string_view event_types[2] = {"before/", "after/"}; // static constexpr std::string_view event_pair_types[2] = {"same/", "mix/"}; + //--------------------------------------------------------------------------- + // Preslices and partitions + o2::framework::SliceCache cache; + o2::framework::PresliceOptional>> perCollision_pcm = o2::aod::v0photonkf::emeventId; + o2::framework::PresliceOptional> perCollision_emc = o2::aod::emccluster::emeventId; + o2::framework::PresliceOptional> perCollision_phos = o2::aod::phoscluster::emeventId; + o2::framework::PresliceOptional>> perCollision_electron = o2::aod::emprimaryelectron::emeventId; + + o2::framework::PresliceOptional perEMCClusterMT = o2::aod::trackmatching::emEmcClusterId; + o2::framework::PresliceOptional perEMCClusterMS = o2::aod::trackmatching::emEmcClusterId; + + o2::framework::Partition>> positrons = o2::aod::emprimaryelectron::sign > int8_t(0) && dileptoncuts.cfg_min_pt_track < o2::aod::track::pt&& nabs(o2::aod::track::eta) < dileptoncuts.cfg_max_eta_track; + o2::framework::Partition>> electrons = o2::aod::emprimaryelectron::sign < int8_t(0) && dileptoncuts.cfg_min_pt_track < o2::aod::track::pt && nabs(o2::aod::track::eta) < dileptoncuts.cfg_max_eta_track; + + o2::aod::pwgem::dilepton::utils::EventMixingHandler, std::pair, o2::aod::pwgem::dilepton::utils::EMTrack>* emh1 = nullptr; + o2::aod::pwgem::dilepton::utils::EventMixingHandler, std::pair, o2::aod::pwgem::dilepton::utils::EMTrack>* emh2 = nullptr; + //--------------------------------------------------------------------------- + + std::vector used_photonIds_per_col; // + std::vector> used_dileptonIds_per_col; // + std::map, uint64_t> map_mixed_eventId_to_globalBC; + std::vector zvtx_bin_edges; std::vector cent_bin_edges; std::vector ep_bin_edges; @@ -215,6 +240,106 @@ struct Pi0EtaToGammaGamma { float d_bz; o2::emcal::Geometry* emcalGeom; + //--------------------------------------------------------------------------- + // In the following are tags defined which help to select the correct preslice and cuts + struct PCMTag { + + static auto& perCollision() + { + static auto slice{o2::aod::v0photonkf::emeventId}; + return slice; + } + + template + static auto& cut(Self& s) + { + return s.fV0PhotonCut; + } + + template + static bool applyCut(Self& s, Photon const& g) + { + return s.fV0PhotonCut.template IsSelected(g); + } + }; + + struct EMCTag { + + static auto& perCollision() + { + static auto slice{o2::aod::emccluster::emeventId}; + return slice; + } + + static auto& perClusterMT() + { + static auto slice{o2::aod::trackmatching::emEmcClusterId}; + return slice; + } + + static auto& perClusterMS() + { + static auto slice{o2::aod::trackmatching::emEmcClusterId}; + return slice; + } + + template + static auto& cut(Self& s) + { + return s.fEMCCut; + } + + // EMCal version has optional tables for matched tracks (global and secondaries) + template + static bool applyCut(Self& s, Cluster const& c, TMatchedTracks const& emcmatchedtracks = nullptr, TMatchedSecondaries const& secondaries = nullptr) + { + return s.fEMCCut.IsSelected(c, emcmatchedtracks, secondaries); + } + }; + + struct PHOSTag { + + static auto& perCollision() + { + static auto slice{o2::aod::phoscluster::emeventId}; + return slice; + } + + template + static auto& cut(Self& s) + { + return s.fPHOSCut; + } + + template + static bool applyCut(Self& s, Cluster const& c) + { + return s.fPHOSCut.IsSelected(c); + } + }; + + struct DalitzEETag { + + static auto& perCollision() + { + static auto slice{o2::aod::emprimaryelectron::emeventId}; + return slice; + } + + template + static auto& cut(Self& s) + { + return s.fDileptonCut; + } + + // Dalitz version has two tracks as argument + B_z + template + static bool applyCut(Self& s, Track1 const& track1, Track2 const& track2, float bz) + { + return s.fDileptonCut.IsSelected(track1, track2, bz); + } + }; + void init(o2::framework::InitContext&) { zvtx_bin_edges = std::vector(ConfVtxBins.value.begin(), ConfVtxBins.value.end()); @@ -407,6 +532,8 @@ struct Pi0EtaToGammaGamma { fEMCCut.SetMinEoverP(emccuts.EMC_Eoverp); fEMCCut.SetUseExoticCut(emccuts.EMC_UseExoticCut); + fEMCCut.SetUseTM(emccuts.emcUseTM.value); // disables or enables TM + fEMCCut.SetUseSecondaryTM(emccuts.emcUseSecondaryTM.value); // disables or enables secondary TM } void DefinePHOSCut() @@ -492,7 +619,7 @@ struct Pi0EtaToGammaGamma { // only combine rotated photons with other photons continue; } - if (!(fEMCCut.IsSelected::iterator>(photon))) { + if (!(fEMCCut.IsSelected(photon))) { continue; } @@ -513,27 +640,25 @@ struct Pi0EtaToGammaGamma { return; } - o2::framework::SliceCache cache; - o2::framework::Preslice>> perCollision_pcm = o2::aod::v0photonkf::emeventId; - o2::framework::Preslice> perCollision_emc = o2::aod::emccluster::emeventId; - o2::framework::Preslice> perCollision_phos = o2::aod::phoscluster::emeventId; - - o2::framework::Preslice>> perCollision_electron = o2::aod::emprimaryelectron::emeventId; - o2::framework::Partition>> positrons = o2::aod::emprimaryelectron::sign > int8_t(0) && dileptoncuts.cfg_min_pt_track < o2::aod::track::pt&& nabs(o2::aod::track::eta) < dileptoncuts.cfg_max_eta_track; - o2::framework::Partition>> electrons = o2::aod::emprimaryelectron::sign < int8_t(0) && dileptoncuts.cfg_min_pt_track < o2::aod::track::pt && nabs(o2::aod::track::eta) < dileptoncuts.cfg_max_eta_track; - - o2::aod::pwgem::dilepton::utils::EventMixingHandler, std::pair, o2::aod::pwgem::dilepton::utils::EMTrack>* emh1 = nullptr; - o2::aod::pwgem::dilepton::utils::EventMixingHandler, std::pair, o2::aod::pwgem::dilepton::utils::EMTrack>* emh2 = nullptr; - std::vector used_photonIds_per_col; // - std::vector> used_dileptonIds_per_col; // - std::map, uint64_t> map_mixed_eventId_to_globalBC; - - template + /// \brief function to run the photon pairing + /// \tparam TDetectorTag1 tag for TPhotons1 type to select the proper cut function and arguments + /// \tparam TDetectorTag2 tag for TPhotons2 type to select the proper cut function and arguments + /// \tparam TCombinationPolicy StrictlyUpperIndex or FullIndex depending if we combine from the same detector or different detectors + /// \tparam TCollisions collision table type + /// \tparam TPhotons1 first photon table type + /// \tparam TPhotons2 second photon table type + /// \tparam TLegs V0 leg table type used in TCut1 and/or TCut2 (only for PCM/DalitzEE) + /// \tparam TMatchedTracks matched global track table type for EMCal (only for EMC) + /// \tparam TMatchedSecondaries matched secondary track table type for EMCal (only for EMC) + /// \param collisions table of collisions (TCollisions) + /// \param photons1 table of photons (TPhotons1) + /// \param photons2 table of photons (TPhotons2) + /// \param legs placeholder argument used only for template deduction (PCM/DalitzEE) + /// \param matchedtracks table of matched global tracks to EMCal clusters (optional) + /// \param matchedsecondaries table of matched secondary tracks to EMCal clusters (optional) + template class TCombinationPolicy = o2::soa::CombinationsStrictlyUpperIndexPolicy, o2::soa::is_table TCollisions, o2::soa::is_table TPhotons1, o2::soa::is_table TPhotons2, typename TLegs = std::nullptr_t, typename TMatchedTracks = std::nullptr_t, typename TMatchedSecondaries = std::nullptr_t> void runPairing(TCollisions const& collisions, - TPhotons1 const& photons1, TPhotons2 const& photons2, - TSubInfos1 const& /*subinfos1*/, TSubInfos2 const& /*subinfos2*/, - TPreslice1 const& perCollision1, TPreslice2 const& perCollision2, - TCut1 const& cut1, TCut2 const& cut2) + TPhotons1 const& photons1, TPhotons2 const& photons2, TLegs const& /*legs*/ = nullptr, TMatchedTracks const& matchedTracks = nullptr, TMatchedSecondaries const& matchedSecondaries = nullptr) { for (const auto& collision : collisions) { initCCDB(collision); @@ -607,59 +732,26 @@ struct Pi0EtaToGammaGamma { std::tuple key_bin = std::make_tuple(zbin, centbin, epbin, occbin); std::pair key_df_collision = std::make_pair(ndf, collision.globalIndex()); - if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMPCM || pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPHOSPHOS || pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kEMCEMC) { // same kinds pairing - auto photons1_per_collision = photons1.sliceBy(perCollision1, collision.globalIndex()); - auto photons2_per_collision = photons2.sliceBy(perCollision2, collision.globalIndex()); - - for (const auto& [g1, g2] : o2::soa::combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(photons1_per_collision, photons2_per_collision))) { - if (!cut1.template IsSelected(g1) || !cut2.template IsSelected(g2)) { - continue; - } - - ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); - ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); - ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; - if (std::fabs(v12.Rapidity()) > maxY) { - continue; - } - - if (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kEMCEMC) { - float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P())); - if (openingAngle < emccuts.minOpenAngle) { - continue; - } - } - - fRegistry.fill(HIST("Pair/same/hs"), v12.M(), v12.Pt(), weight); - - if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kEMCEMC) { - RotationBackground>(v12, v1, v2, photons2_per_collision, g1.globalIndex(), g2.globalIndex(), weight); - } - - if (std::find(used_photonIds_per_col.begin(), used_photonIds_per_col.end(), g1.globalIndex()) == used_photonIds_per_col.end()) { - emh1->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::dilepton::utils::EMTrack(g1.pt(), g1.eta(), g1.phi(), 0)); - used_photonIds_per_col.emplace_back(g1.globalIndex()); - } - if (std::find(used_photonIds_per_col.begin(), used_photonIds_per_col.end(), g2.globalIndex()) == used_photonIds_per_col.end()) { - emh1->AddTrackToEventPool(key_df_collision, o2::aod::pwgem::dilepton::utils::EMTrack(g2.pt(), g2.eta(), g2.phi(), 0)); - used_photonIds_per_col.emplace_back(g2.globalIndex()); - } - ndiphoton++; - } // end of pairing loop - } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMDalitzEE) { - auto photons1_per_collision = photons1.sliceBy(perCollision1, collision.globalIndex()); + if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMDalitzEE) { + auto photons1_per_collision = photons1.sliceByCached(TDetectorTag1::perCollision(), collision.globalIndex(), cache); auto positrons_per_collision = positrons->sliceByCached(o2::aod::emprimaryelectron::emeventId, collision.globalIndex(), cache); auto electrons_per_collision = electrons->sliceByCached(o2::aod::emprimaryelectron::emeventId, collision.globalIndex(), cache); for (const auto& g1 : photons1_per_collision) { - if (!cut1.template IsSelected(g1)) { - continue; + if constexpr (std::is_same_v) { + if (!TDetectorTag1::applyCut(*this, g1)) { + continue; + } + } else { + if (!TDetectorTag2::applyCut(*this, g1)) { + continue; + } } - auto pos1 = g1.template posTrack_as(); - auto ele1 = g1.template negTrack_as(); + auto pos1 = g1.template posTrack_as(); + auto ele1 = g1.template negTrack_as(); ROOT::Math::PtEtaPhiMVector v_gamma(g1.pt(), g1.eta(), g1.phi(), 0.); - for (const auto& [pos2, ele2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(positrons_per_collision, electrons_per_collision))) { + for (const auto& [pos2, ele2] : o2::soa::combinations(TCombinationPolicy(positrons_per_collision, electrons_per_collision))) { if (pos2.trackId() == ele2.trackId()) { // this is protection against pairing identical 2 tracks. continue; @@ -668,13 +760,22 @@ struct Pi0EtaToGammaGamma { continue; } - if (!cut2.template IsSelectedTrack(pos2, collision) || !cut2.template IsSelectedTrack(ele2, collision)) { - continue; + if constexpr (std::is_same_v) { + if (!TDetectorTag2::applyCut(*this, pos2, ele2, d_bz)) { + continue; + } + } else { + if (!TDetectorTag1::applyCut(*this, pos2, ele2, d_bz)) { + continue; + } } + // if (!cut2.template IsSelectedTrack(pos2, collision) || !cut2.template IsSelectedTrack(ele2, collision)) { + // continue; + // } - if (!cut2.IsSelectedPair(pos2, ele2, d_bz)) { - continue; - } + // if (!cut2.IsSelectedPair(pos2, ele2, d_bz)) { + // continue; + // } ROOT::Math::PtEtaPhiMVector v_pos(pos2.pt(), pos2.eta(), pos2.phi(), o2::constants::physics::MassElectron); ROOT::Math::PtEtaPhiMVector v_ele(ele2.pt(), ele2.eta(), ele2.phi(), o2::constants::physics::MassElectron); @@ -698,14 +799,37 @@ struct Pi0EtaToGammaGamma { ndiphoton++; } // end of dielectron loop } // end of g1 loop - } else { // PCM-EMC, PCM-PHOS. Nightmare. don't run these pairs. - auto photons1_per_collision = photons1.sliceBy(perCollision1, collision.globalIndex()); - auto photons2_per_collision = photons2.sliceBy(perCollision2, collision.globalIndex()); + } else { // PCM-PCM, EMC-EMC, PHOS-PHOS, PCM-EMC and PCM-PHOS. + auto photons1_per_collision = photons1.sliceByCached(TDetectorTag1::perCollision(), collision.globalIndex(), cache); + auto photons2_per_collision = photons2.sliceByCached(TDetectorTag2::perCollision(), collision.globalIndex(), cache); - for (const auto& [g1, g2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(photons1_per_collision, photons2_per_collision))) { - if (!cut1.template IsSelected(g1) || !cut2.template IsSelected(g2)) { - continue; + for (const auto& [g1, g2] : o2::soa::combinations(TCombinationPolicy(photons1_per_collision, photons2_per_collision))) { + if constexpr (std::is_same_v) { + // For the EMCal case we need to get the primary and secondary matched tracks + auto matchedTracks1 = matchedTracks.sliceByCached(TDetectorTag1::perClusterMT(), g1.globalIndex(), cache); + auto matchedSecondaries1 = matchedSecondaries.sliceByCached(TDetectorTag1::perClusterMS(), g1.globalIndex(), cache); + + if (!TDetectorTag1::applyCut(*this, g1, matchedTracks1, matchedSecondaries1)) { + continue; + } + } else { + if (!TDetectorTag1::applyCut(*this, g1) || !TDetectorTag2::applyCut(*this, g2)) { + continue; + } } + + if constexpr (std::is_same_v) { + auto matchedTracks2 = matchedTracks.sliceByCached(TDetectorTag2::perClusterMT(), g2.globalIndex(), cache); + auto matchedSecondaries2 = matchedSecondaries.sliceByCached(TDetectorTag2::perClusterMS(), g2.globalIndex(), cache); + if (!TDetectorTag2::applyCut(*this, g2, matchedTracks2, matchedSecondaries2)) { + continue; + } + } else { + if (!TDetectorTag2::applyCut(*this, g2)) { + continue; + } + } + ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; @@ -868,34 +992,34 @@ struct Pi0EtaToGammaGamma { { // LOGF(info, "ndf = %d", ndf); if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMPCM) { - auto v0photons = std::get<0>(std::tie(args...)); - auto v0legs = std::get<1>(std::tie(args...)); - runPairing(collisions, v0photons, v0photons, v0legs, v0legs, perCollision_pcm, perCollision_pcm, fV0PhotonCut, fV0PhotonCut); + auto&& [v0photons, v0legs] = std::forward_as_tuple(args...); + // auto v0photons = std::get<0>(std::tie(args...)); + // auto v0legs = std::get<1>(std::tie(args...)); + runPairing(collisions, v0photons, v0photons, v0legs); } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMDalitzEE) { - auto v0photons = std::get<0>(std::tie(args...)); - auto v0legs = std::get<1>(std::tie(args...)); - auto emprimaryelectrons = std::get<2>(std::tie(args...)); + auto&& [v0photons, v0legs, emprimaryelectrons] = std::forward_as_tuple(args...); + // auto v0photons = std::get<0>(std::tie(args...)); + // auto v0legs = std::get<1>(std::tie(args...)); + // auto emprimaryelectrons = std::get<2>(std::tie(args...)); // LOGF(info, "electrons.size() = %d, positrons.size() = %d", electrons.size(), positrons.size()); - runPairing(collisions, v0photons, emprimaryelectrons, v0legs, emprimaryelectrons, perCollision_pcm, perCollision_electron, fV0PhotonCut, fDileptonCut); + runPairing(collisions, v0photons, emprimaryelectrons, v0legs); } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kEMCEMC) { - auto emcclusters = std::get<0>(std::tie(args...)); - runPairing(collisions, emcclusters, emcclusters, nullptr, nullptr, perCollision_emc, perCollision_emc, fEMCCut, fEMCCut); + auto&& [emcClusters, emcMatchedTracks, emcMatchedSecondaries] = std::forward_as_tuple(args...); + // auto emcClusters = std::get<0>(std::tie(args...)); + // auto emcMatchedTracks = std::get<1>(std::tie(args...)); + // auto emcMatchedSecondaries = std::get<2>(std::tie(args...)); + runPairing(collisions, emcClusters, emcClusters, nullptr, emcMatchedTracks, emcMatchedSecondaries); } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPHOSPHOS) { - auto phosclusters = std::get<0>(std::tie(args...)); - runPairing(collisions, phosclusters, phosclusters, nullptr, nullptr, perCollision_phos, perCollision_phos, fPHOSCut, fPHOSCut); + auto&& [phosclusters] = std::forward_as_tuple(args...); + // auto phosclusters = std::get<0>(std::tie(args...)); + runPairing(collisions, phosclusters, phosclusters); + } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMEMC) { + auto&& [v0photons, emcclusters, v0legs, emcMatchedTracks, emcMatchedSecondaries] = std::forward_as_tuple(args...); + runPairing(collisions, v0photons, emcclusters, v0legs, emcMatchedTracks, emcMatchedSecondaries); + } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMPHOS) { + auto&& [v0photons, v0legs, phosclusters] = std::forward_as_tuple(args...); + runPairing(collisions, v0photons, phosclusters, v0legs); } - // else if constexpr (pairtype == PairType::kPCMEMC) { - // auto v0photons = std::get<0>(std::tie(args...)); - // auto v0legs = std::get<1>(std::tie(args...)); - // auto emcclusters = std::get<2>(std::tie(args...)); - // auto emcmatchedtracks = std::get<3>(std::tie(args...)); - // runPairing(collisions, v0photons, emcclusters, v0legs, nullptr, perCollision_pcm, perCollision_emc, fV0PhotonCut, fEMCCut, emcmatchedtracks, nullptr); - // } else if constexpr (pairtype == PairType::kPCMPHOS) { - // auto v0photons = std::get<0>(std::tie(args...)); - // auto v0legs = std::get<1>(std::tie(args...)); - // auto phosclusters = std::get<2>(std::tie(args...)); - // runPairing(collisions, v0photons, phosclusters, v0legs, nullptr, perCollision_pcm, perCollision_phos, fV0PhotonCut, fPHOSCut, nullptr, nullptr); - // } ndf++; } PROCESS_SWITCH(Pi0EtaToGammaGamma, processAnalysis, "process pair analysis", true); @@ -905,34 +1029,25 @@ struct Pi0EtaToGammaGamma { { // LOGF(info, "ndf = %d", ndf); if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMPCM) { - auto v0photons = std::get<0>(std::tie(args...)); - auto v0legs = std::get<1>(std::tie(args...)); - runPairing(collisions, v0photons, v0photons, v0legs, v0legs, perCollision_pcm, perCollision_pcm, fV0PhotonCut, fV0PhotonCut); + auto&& [v0photons, v0legs] = std::forward_as_tuple(args...); + runPairing(collisions, v0photons, v0photons, v0legs); } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMDalitzEE) { - auto v0photons = std::get<0>(std::tie(args...)); - auto v0legs = std::get<1>(std::tie(args...)); - auto emprimaryelectrons = std::get<2>(std::tie(args...)); - // LOGF(info, "electrons.size() = %d, positrons.size() = %d", electrons.size(), positrons.size()); - runPairing(collisions, v0photons, emprimaryelectrons, v0legs, emprimaryelectrons, perCollision_pcm, perCollision_electron, fV0PhotonCut, fDileptonCut); + auto&& [v0photons, v0legs, emprimaryelectrons] = std::forward_as_tuple(args...); + runPairing(collisions, v0photons, emprimaryelectrons, v0legs); } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kEMCEMC) { - auto emcclusters = std::get<0>(std::tie(args...)); - runPairing(collisions, emcclusters, emcclusters, nullptr, nullptr, perCollision_emc, perCollision_emc, fEMCCut, fEMCCut); + auto&& [emcClusters, emcMatchedTracks, emcMatchedSecondaries] = std::forward_as_tuple(args...); + runPairing(collisions, emcClusters, emcClusters, nullptr, emcMatchedTracks, emcMatchedSecondaries); } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPHOSPHOS) { - auto phosclusters = std::get<0>(std::tie(args...)); - runPairing(collisions, phosclusters, phosclusters, nullptr, nullptr, perCollision_phos, perCollision_phos, fPHOSCut, fPHOSCut); + auto&& [phosclusters] = std::forward_as_tuple(args...); + // auto phosclusters = std::get<0>(std::tie(args...)); + runPairing(collisions, phosclusters, phosclusters); + } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMEMC) { + auto&& [v0photons, emcclusters, v0legs, emcMatchedTracks, emcMatchedSecondaries] = std::forward_as_tuple(args...); + runPairing(collisions, v0photons, emcclusters, v0legs, emcMatchedTracks, emcMatchedSecondaries); + } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMPHOS) { + auto&& [v0photons, phosclusters, v0legs] = std::forward_as_tuple(args...); + runPairing(collisions, v0photons, phosclusters, v0legs); } - // else if constexpr (pairtype == PairType::kPCMEMC) { - // auto v0photons = std::get<0>(std::tie(args...)); - // auto v0legs = std::get<1>(std::tie(args...)); - // auto emcclusters = std::get<2>(std::tie(args...)); - // auto emcmatchedtracks = std::get<3>(std::tie(args...)); - // runPairing(collisions, v0photons, emcclusters, v0legs, nullptr, perCollision_pcm, perCollision_emc, fV0PhotonCut, fEMCCut, emcmatchedtracks, nullptr); - // } else if constexpr (pairtype == PairType::kPCMPHOS) { - // auto v0photons = std::get<0>(std::tie(args...)); - // auto v0legs = std::get<1>(std::tie(args...)); - // auto phosclusters = std::get<2>(std::tie(args...)); - // runPairing(collisions, v0photons, phosclusters, v0legs, nullptr, perCollision_pcm, perCollision_phos, fV0PhotonCut, fPHOSCut, nullptr, nullptr); - // } ndf++; } PROCESS_SWITCH(Pi0EtaToGammaGamma, processAnalysisJJMC, "process pair analysis", false); diff --git a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h index df416caf9d4..586ea25c9e6 100644 --- a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h +++ b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h @@ -56,13 +56,14 @@ #include #include +#include #include #include #include #include #include -template +template struct Pi0EtaToGammaGammaMC { o2::framework::Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; o2::framework::Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; @@ -180,6 +181,8 @@ struct Pi0EtaToGammaGammaMC { o2::framework::Configurable> EMC_TM_Phi{"EMC_TM_Phi", {0.015f, 3.65f, -2.f}, "|phi| <= [0]+(pT+[1])^[2] for EMCal track matching"}; o2::framework::Configurable EMC_Eoverp{"EMC_Eoverp", 1.75, "Minimum cluster energy over track momentum for EMCal track matching"}; o2::framework::Configurable EMC_UseExoticCut{"EMC_UseExoticCut", true, "FLag to use the EMCal exotic cluster cut"}; + o2::framework::Configurable emcUseTM{"emcUseTM", true, "flag to use EMCal track matching cut or not."}; + o2::framework::Configurable emcUseSecondaryTM{"emcUseSecondaryTM", false, "flag to use EMCal secondary track matching cut or not. Required for PCM-EMC analyses."}; } emccuts; o2::framework::Configurable maxY_gen{"maxY_gen", 0.9, "maximum rapidity for generated particles"}; // for PCM and dielectron @@ -369,6 +372,8 @@ struct Pi0EtaToGammaGammaMC { fEMCCut.SetMinEoverP(emccuts.EMC_Eoverp); fEMCCut.SetUseExoticCut(emccuts.EMC_UseExoticCut); + fEMCCut.SetUseTM(emccuts.emcUseTM.value); // disables or enables TM + fEMCCut.SetUseSecondaryTM(emccuts.emcUseSecondaryTM.value); // disables or enables secondary TM } void DefinePHOSCut() @@ -377,21 +382,141 @@ struct Pi0EtaToGammaGammaMC { } o2::framework::SliceCache cache; - o2::framework::Preslice>> perCollision_pcm = o2::aod::v0photonkf::emeventId; - o2::framework::Preslice> perCollision_emc = o2::aod::emccluster::emeventId; - o2::framework::Preslice> perCollision_phos = o2::aod::phoscluster::emeventId; + o2::framework::PresliceOptional>> perCollision_pcm = o2::aod::v0photonkf::emeventId; + o2::framework::PresliceOptional> perCollision_emc = o2::aod::emccluster::emeventId; + o2::framework::PresliceOptional> perCollision_phos = o2::aod::phoscluster::emeventId; + o2::framework::PresliceOptional>> perCollision_electron = o2::aod::emprimaryelectron::emeventId; + + o2::framework::PresliceOptional perEMCClusterMT = o2::aod::trackmatching::emEmcClusterId; + o2::framework::PresliceOptional perEMCClusterMS = o2::aod::trackmatching::emEmcClusterId; - o2::framework::Preslice>> perCollision_electron = o2::aod::emprimaryelectron::emeventId; o2::framework::Partition>> positrons = o2::aod::emprimaryelectron::sign > int8_t(0) && dileptoncuts.cfg_min_pt_track < o2::aod::track::pt&& nabs(o2::aod::track::eta) < dileptoncuts.cfg_max_eta_track; o2::framework::Partition>> electrons = o2::aod::emprimaryelectron::sign < int8_t(0) && dileptoncuts.cfg_min_pt_track < o2::aod::track::pt && nabs(o2::aod::track::eta) < dileptoncuts.cfg_max_eta_track; - template + //--------------------------------------------------------------------------- + // In the following are tags defined which help to select the correct preslice and cuts + struct PCMTag { + + static auto& perCollision() + { + static auto slice{o2::aod::v0photonkf::emeventId}; + return slice; + } + + template + static auto& cut(Self& s) + { + return s.fV0PhotonCut; + } + + template + static bool applyCut(Self& s, Photon const& g) + { + return s.fV0PhotonCut.template IsSelected>(g); + } + }; + + struct EMCTag { + + static auto& perCollision() + { + static auto slice{o2::aod::emccluster::emeventId}; + return slice; + } + + static auto& perClusterMT() + { + static auto slice{o2::aod::trackmatching::emEmcClusterId}; + return slice; + } + + static auto& perClusterMS() + { + static auto slice{o2::aod::trackmatching::emEmcClusterId}; + return slice; + } + + template + static auto& cut(Self& s) + { + return s.fEMCCut; + } + + // EMCal version has optional tables for matched tracks (global and secondaries) + template + static bool applyCut(Self& s, Cluster const& c, TMatchedTracks const& emcmatchedtracks = nullptr, TMatchedSecondaries const& secondaries = nullptr) + { + return s.fEMCCut.IsSelected(c, emcmatchedtracks, secondaries); + } + }; + + struct PHOSTag { + + static auto& perCollision() + { + static auto slice{o2::aod::phoscluster::emeventId}; + return slice; + } + + template + static auto& cut(Self& s) + { + return s.fPHOSCut; + } + + template + static bool applyCut(Self& s, Cluster const& c) + { + return s.fPHOSCut.IsSelected(c); + } + }; + + struct DalitzEETag { + static auto& perCollision() + { + static auto slice{o2::aod::emprimaryelectron::emeventId}; + return slice; + } + + template + static auto& cut(Self& s) + { + return s.fDileptonCut; + } + + // Dalitz version has two tracks as argument + B_z + template + static bool applyCut(Self& s, Track1 const& track1, Track2 const& track2, float bz) + { + return s.fDileptonCut.IsSelected(track1, track2, bz); + } + }; + + /// \brief function to run the photon pairing + /// \tparam TDetectorTag1 tag for TPhotons1 type to select the proper cut function and arguments + /// \tparam TDetectorTag2 tag for TPhotons2 type to select the proper cut function and arguments + /// \tparam TCombinationPolicy StrictlyUpperIndex or FullIndex depending if we combine from the same detector or different detectors + /// \tparam TCollisions collision table type + /// \tparam TPhotons1 first photon table type + /// \tparam TPhotons2 second photon table type + /// \tparam TMCCollisions mc collision table type + /// \tparam TMCParticles mc particles table type + /// \tparam TLeg V0 leg table type used in TCut1 and/or TCut2 (for PCM/DalitzEE) + /// \tparam TMatched matched global track table type for EMCal (for EMC) + /// \tparam TSecondaries matched secondary track table type for EMCal (for EMC) + /// \param collisions table of collisions (TCollisions) + /// \param photons1 table of photons (TPhotons1) + /// \param photons2 table of photons (TPhotons2) + /// \param mccollisions table of mc collisions (TMCCollisions) + /// \param mcparticles table of mc particles (TMCParticles) + /// \param legs placeholder argument used only for template deduction (for PCM/DalitzEE) + /// \param matchedtracks table of matched global tracks to EMCal clusters (for EMC) + /// \param matchedsecondaries table of matched secondary tracks to EMCal clusters (for EMC) + template class TCombinationPolicy = o2::soa::CombinationsStrictlyUpperIndexPolicy, o2::soa::is_table TCollisions, o2::soa::is_table TPhotons1, o2::soa::is_table TPhotons2, o2::soa::is_table TMCCollisions, o2::soa::is_table TMCParticles, typename TLegs = std::nullptr_t, typename TMatchedTracks = std::nullptr_t, typename TMatchedSecondaries = std::nullptr_t> void runTruePairing(TCollisions const& collisions, TPhotons1 const& photons1, TPhotons2 const& photons2, - TSubInfos1 const& /*subinfos1*/, TSubInfos2 const& /*subinfos2*/, - TPreslice1 const& perCollision1, TPreslice2 const& perCollision2, - TCut1 const& cut1, TCut2 const& cut2, - TMCCollisions const& mccollisions, TMCParticles const& mcparticles) + TMCCollisions const& mccollisions, TMCParticles const& mcparticles, + TLegs const& /*legs*/ = nullptr, TMatchedTracks const& matchedTracks = nullptr, TMatchedSecondaries const& matchedSecondaries = nullptr) { for (auto& collision : collisions) { initCCDB(collision); @@ -424,19 +549,42 @@ struct Pi0EtaToGammaGammaMC { int photonid1 = -1, photonid2 = -1, pi0id = -1, etaid = -1; if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMPCM || pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPHOSPHOS || pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kEMCEMC) { // same kinds pairing - auto photons1_per_collision = photons1.sliceBy(perCollision1, collision.globalIndex()); - auto photons2_per_collision = photons2.sliceBy(perCollision2, collision.globalIndex()); - for (auto& [g1, g2] : o2::soa::combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(photons1_per_collision, photons2_per_collision))) { - if (!cut1.template IsSelected(g1) || !cut2.template IsSelected(g2)) { - continue; + auto photons1_per_collision = photons1.sliceByCached(TDetectorTag1::perCollision(), collision.globalIndex(), cache); + auto photons2_per_collision = photons2.sliceByCached(TDetectorTag2::perCollision(), collision.globalIndex(), cache); + + for (const auto& [g1, g2] : o2::soa::combinations(TCombinationPolicy(photons1_per_collision, photons2_per_collision))) { + if constexpr (std::is_same_v) { + // For the EMCal case we need to get the primary and secondary matched tracks + auto matchedTracks1 = matchedTracks.sliceByCached(TDetectorTag1::perClusterMT(), g1.globalIndex(), cache); + auto matchedSecondaries1 = matchedSecondaries.sliceByCached(TDetectorTag1::perClusterMS(), g1.globalIndex(), cache); + + if (!TDetectorTag1::applyCut(*this, g1, matchedTracks1, matchedSecondaries1)) { + continue; + } + } else { + if (!TDetectorTag1::applyCut(*this, g1) || !TDetectorTag2::applyCut(*this, g2)) { + continue; + } + } + + if constexpr (std::is_same_v) { + auto matchedTracks2 = matchedTracks.sliceByCached(TDetectorTag2::perClusterMT(), g2.globalIndex(), cache); + auto matchedSecondaries2 = matchedSecondaries.sliceByCached(TDetectorTag2::perClusterMS(), g2.globalIndex(), cache); + if (!TDetectorTag2::applyCut(*this, g2, matchedTracks2, matchedSecondaries2)) { + continue; + } + } else { + if (!TDetectorTag2::applyCut(*this, g2)) { + continue; + } } if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMPCM) { // check 2 legs - auto pos1 = g1.template posTrack_as(); - auto ele1 = g1.template negTrack_as(); - auto pos2 = g2.template posTrack_as(); - auto ele2 = g2.template negTrack_as(); + auto pos1 = g1.template posTrack_as(); + auto ele1 = g1.template negTrack_as(); + auto pos2 = g2.template posTrack_as(); + auto ele2 = g2.template negTrack_as(); auto pos1mc = pos1.template emmcparticle_as(); auto ele1mc = ele1.template emmcparticle_as(); @@ -512,16 +660,24 @@ struct Pi0EtaToGammaGammaMC { } } // end of pairing loop } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMDalitzEE) { - auto photons1_per_collision = photons1.sliceBy(perCollision1, collision.globalIndex()); + + auto photons1_per_collision = photons1.sliceByCached(TDetectorTag1::perCollision(), collision.globalIndex(), cache); auto positrons_per_collision = positrons->sliceByCached(o2::aod::emprimaryelectron::emeventId, collision.globalIndex(), cache); auto electrons_per_collision = electrons->sliceByCached(o2::aod::emprimaryelectron::emeventId, collision.globalIndex(), cache); - for (auto& g1 : photons1_per_collision) { - if (!cut1.template IsSelected(g1)) { - continue; + for (const auto& g1 : photons1_per_collision) { + if constexpr (std::is_same_v) { + if (!TDetectorTag1::applyCut(*this, g1)) { + continue; + } + } else { + if (!TDetectorTag2::applyCut(*this, g1)) { + continue; + } } - auto pos1 = g1.template posTrack_as(); - auto ele1 = g1.template negTrack_as(); + + auto pos1 = g1.template posTrack_as(); + auto ele1 = g1.template negTrack_as(); auto pos1mc = pos1.template emmcparticle_as(); auto ele1mc = ele1.template emmcparticle_as(); photonid1 = o2::aod::pwgem::dilepton::utils::mcutil::FindCommonMotherFrom2Prongs(pos1mc, ele1mc, -11, 11, 22, mcparticles); @@ -534,20 +690,22 @@ struct Pi0EtaToGammaGammaMC { } ROOT::Math::PtEtaPhiMVector v_gamma(g1.pt(), g1.eta(), g1.phi(), 0.f); - for (auto& [pos2, ele2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(positrons_per_collision, electrons_per_collision))) { // ULS - if (pos2.trackId() == ele2.trackId()) { // this is protection against pairing identical 2 tracks. + for (const auto& [pos2, ele2] : o2::soa::combinations(TCombinationPolicy(positrons_per_collision, electrons_per_collision))) { // ULS + if (pos2.trackId() == ele2.trackId()) { // this is protection against pairing identical 2 tracks. continue; } if (pos1.trackId() == pos2.trackId() || ele1.trackId() == ele2.trackId()) { continue; } - if (!cut2.template IsSelectedTrack(pos2, collision) || !cut2.template IsSelectedTrack(ele2, collision)) { - continue; - } - - if (!cut2.IsSelectedPair(pos2, ele2, d_bz)) { - continue; + if constexpr (std::is_same_v) { + if (!TDetectorTag2::applyCut(*this, pos2, ele2, d_bz)) { + continue; + } + } else { + if (!TDetectorTag1::applyCut(*this, pos2, ele2, d_bz)) { + continue; + } } auto pos2mc = mcparticles.iteratorAt(pos2.emmcparticleId()); @@ -578,14 +736,38 @@ struct Pi0EtaToGammaGammaMC { } } // end of dielectron loop } // end of pcm loop - } else { // PCM-EMC, PCM-PHOS. Nightmare. don't run these pairs. - auto photons1_per_collision = photons1.sliceBy(perCollision1, collision.globalIndex()); - auto photons2_per_collision = photons2.sliceBy(perCollision2, collision.globalIndex()); + } else { // PCM-EMC, PCM-PHOS. + // TODO: implement proper functionality if we ever want to run this in Pb-Pb + auto photons1_per_collision = photons1.sliceByCached(TDetectorTag1::perCollision(), collision.globalIndex(), cache); + auto photons2_per_collision = photons2.sliceByCached(TDetectorTag2::perCollision(), collision.globalIndex(), cache); + + for (const auto& [g1, g2] : o2::soa::combinations(TCombinationPolicy(photons1_per_collision, photons2_per_collision))) { + if constexpr (std::is_same_v) { + // For the EMCal case we need to get the primary and secondary matched tracks + auto matchedTracks1 = matchedTracks.sliceByCached(TDetectorTag1::perClusterMT(), g1.globalIndex(), cache); + auto matchedSecondaries1 = matchedSecondaries.sliceByCached(TDetectorTag1::perClusterMS(), g1.globalIndex(), cache); + + if (!TDetectorTag1::applyCut(*this, g1, matchedTracks1, matchedSecondaries1)) { + continue; + } + } else { + if (!TDetectorTag1::applyCut(*this, g1) || !TDetectorTag2::applyCut(*this, g2)) { + continue; + } + } - for (auto& [g1, g2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(photons1_per_collision, photons2_per_collision))) { - if (!cut1.template IsSelected(g1) || !cut2.template IsSelected(g2)) { - continue; + if constexpr (std::is_same_v) { + auto matchedTracks2 = matchedTracks.sliceByCached(TDetectorTag2::perClusterMT(), g2.globalIndex(), cache); + auto matchedSecondaries2 = matchedSecondaries.sliceByCached(TDetectorTag2::perClusterMS(), g2.globalIndex(), cache); + if (!TDetectorTag2::applyCut(*this, g2, matchedTracks2, matchedSecondaries2)) { + continue; + } + } else { + if (!TDetectorTag2::applyCut(*this, g2)) { + continue; + } } + ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; @@ -684,20 +866,16 @@ struct Pi0EtaToGammaGammaMC { void processAnalysis(o2::soa::Filtered> const& collisions, o2::soa::Join const& mccollisions, o2::aod::EMMCParticles const& mcparticles, Types const&... args) { if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMPCM) { - auto v0photons = std::get<0>(std::tie(args...)); - auto v0legs = std::get<1>(std::tie(args...)); - runTruePairing(collisions, v0photons, v0photons, v0legs, v0legs, perCollision_pcm, perCollision_pcm, fV0PhotonCut, fV0PhotonCut, mccollisions, mcparticles); + auto&& [v0photons, v0legs] = std::forward_as_tuple(args...); + runTruePairing(collisions, v0photons, v0photons, mccollisions, mcparticles, v0legs); runGenInfo(collisions, mccollisions, mcparticles); } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMDalitzEE) { - auto v0photons = std::get<0>(std::tie(args...)); - auto v0legs = std::get<1>(std::tie(args...)); - auto emprimaryelectrons = std::get<2>(std::tie(args...)); - // LOGF(info, "electrons.size() = %d, positrons.size() = %d", electrons.size(), positrons.size()); - runTruePairing(collisions, v0photons, emprimaryelectrons, v0legs, emprimaryelectrons, perCollision_pcm, perCollision_electron, fV0PhotonCut, fDileptonCut, mccollisions, mcparticles); + auto&& [v0photons, v0legs, emprimaryelectrons] = std::forward_as_tuple(args...); + runTruePairing(collisions, v0photons, emprimaryelectrons, mccollisions, mcparticles, v0legs); runGenInfo(collisions, mccollisions, mcparticles); } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kEMCEMC) { - auto emcclusters = std::get<0>(std::tie(args...)); - runTruePairing(collisions, emcclusters, emcclusters, nullptr, nullptr, perCollision_emc, perCollision_emc, fEMCCut, fEMCCut, mccollisions, mcparticles); + auto&& [emcclusters, emcMatchedTracks, emcMatchedSecondaries] = std::forward_as_tuple(args...); + runTruePairing(collisions, emcclusters, emcclusters, mccollisions, mcparticles, nullptr, emcMatchedTracks, emcMatchedSecondaries); runGenInfo(collisions, mccollisions, mcparticles); } @@ -724,20 +902,26 @@ struct Pi0EtaToGammaGammaMC { void processAnalysisJJMC(o2::soa::Filtered, o2::aod::EMEventsWeight>> const& collisions, o2::soa::Join const& mccollisions, o2::aod::EMMCParticles const& mcparticles, Types const&... args) { if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMPCM) { - auto v0photons = std::get<0>(std::tie(args...)); - auto v0legs = std::get<1>(std::tie(args...)); - runTruePairing(collisions, v0photons, v0photons, v0legs, v0legs, perCollision_pcm, perCollision_pcm, fV0PhotonCut, fV0PhotonCut, mccollisions, mcparticles); + auto&& [v0photons, v0legs] = std::forward_as_tuple(args...); + // auto v0photons = std::get<0>(std::tie(args...)); + // auto v0legs = std::get<1>(std::tie(args...)); + // runTruePairing(collisions, v0photons, v0photons, v0legs, v0legs, perCollision_pcm, perCollision_pcm, fV0PhotonCut, fV0PhotonCut, mccollisions, mcparticles); + runTruePairing(collisions, v0photons, v0photons, mccollisions, mcparticles, v0legs); runGenInfo(collisions, mccollisions, mcparticles); } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMDalitzEE) { - auto v0photons = std::get<0>(std::tie(args...)); - auto v0legs = std::get<1>(std::tie(args...)); - auto emprimaryelectrons = std::get<2>(std::tie(args...)); + auto&& [v0photons, v0legs, emprimaryelectrons] = std::forward_as_tuple(args...); + // auto v0photons = std::get<0>(std::tie(args...)); + // auto v0legs = std::get<1>(std::tie(args...)); + // auto emprimaryelectrons = std::get<2>(std::tie(args...)); // LOGF(info, "electrons.size() = %d, positrons.size() = %d", electrons.size(), positrons.size()); - runTruePairing(collisions, v0photons, emprimaryelectrons, v0legs, emprimaryelectrons, perCollision_pcm, perCollision_electron, fV0PhotonCut, fDileptonCut, mccollisions, mcparticles); + // runTruePairing(collisions, v0photons, emprimaryelectrons, v0legs, emprimaryelectrons, perCollision_pcm, perCollision_electron, fV0PhotonCut, fDileptonCut, mccollisions, mcparticles); + runTruePairing(collisions, v0photons, emprimaryelectrons, mccollisions, mcparticles, v0legs); runGenInfo(collisions, mccollisions, mcparticles); } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kEMCEMC) { - auto emcclusters = std::get<0>(std::tie(args...)); - runTruePairing(collisions, emcclusters, emcclusters, nullptr, nullptr, perCollision_emc, perCollision_emc, fEMCCut, fEMCCut, mccollisions, mcparticles); + auto&& [emcclusters, emcMatchedTracks, emcMatchedSecondaries] = std::forward_as_tuple(args...); + // auto emcclusters = std::get<0>(std::tie(args...)); + // runTruePairing(collisions, emcclusters, emcclusters, nullptr, nullptr, perCollision_emc, perCollision_emc, fEMCCut, fEMCCut, mccollisions, mcparticles); + runTruePairing(collisions, emcclusters, emcclusters, mccollisions, mcparticles, nullptr, emcMatchedTracks, emcMatchedSecondaries); runGenInfo(collisions, mccollisions, mcparticles); } diff --git a/PWGEM/PhotonMeson/Core/TaggingPi0.h b/PWGEM/PhotonMeson/Core/TaggingPi0.h index e374a21b26b..681b9d9938d 100644 --- a/PWGEM/PhotonMeson/Core/TaggingPi0.h +++ b/PWGEM/PhotonMeson/Core/TaggingPi0.h @@ -539,7 +539,7 @@ struct TaggingPi0 { auto electrons_per_collision = electrons->sliceByCached(o2::aod::emprimaryelectron::emeventId, collision.globalIndex(), cache); // electrons for (const auto& g1 : photons1_per_collision) { - if (!cut1.template IsSelected(g1)) { + if (!cut1.template IsSelected(g1)) { continue; } auto pos1 = g1.template posTrack_as(); @@ -588,7 +588,7 @@ struct TaggingPi0 { auto photons2_per_collision = photons2.sliceBy(perCollision2, collision.globalIndex()); // EMC or PHOS for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(photons1_per_collision, photons2_per_collision))) { - if (!cut1.template IsSelected(g1) || !cut2.template IsSelected(g2)) { + if (!cut1.template IsSelected(g1) || !cut2.template IsSelected(g2)) { continue; } ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); diff --git a/PWGEM/PhotonMeson/Core/TaggingPi0MC.h b/PWGEM/PhotonMeson/Core/TaggingPi0MC.h index c13b531bfdc..fb9ba22d87f 100644 --- a/PWGEM/PhotonMeson/Core/TaggingPi0MC.h +++ b/PWGEM/PhotonMeson/Core/TaggingPi0MC.h @@ -465,7 +465,7 @@ struct TaggingPi0MC { auto electrons_per_collision = electrons->sliceByCached(o2::aod::emprimaryelectron::emeventId, collision.globalIndex(), cache); for (const auto& g1 : photons1_per_collision) { - if (!cut1.template IsSelected(g1)) { + if (!cut1.template IsSelected(g1)) { continue; } @@ -562,7 +562,7 @@ struct TaggingPi0MC { auto photons2_per_collision = photons2.sliceBy(perCollision2, collision.globalIndex()); for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(photons1_per_collision, photons2_per_collision))) { - if (!cut1.template IsSelected(g1) || !cut2.template IsSelected(g2)) { + if (!cut1.template IsSelected(g1) || !cut2.template IsSelected(g2)) { continue; } // ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); diff --git a/PWGEM/PhotonMeson/Core/V0PhotonCut.cxx b/PWGEM/PhotonMeson/Core/V0PhotonCut.cxx index 4c149b7fd20..bebd8df962c 100644 --- a/PWGEM/PhotonMeson/Core/V0PhotonCut.cxx +++ b/PWGEM/PhotonMeson/Core/V0PhotonCut.cxx @@ -9,16 +9,21 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -// -// Class for v0 photon selection -// +/// \file V0PhotonCut.cxx +/// \brief Source of class for V0 photon selection. +/// \author D. Sekihata, daiki.sekihata@cern.ch + +#include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" +#include + +#include + +#include +#include #include #include -#include "Framework/Logger.h" -#include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" - ClassImp(V0PhotonCut); const std::pair> V0PhotonCut::its_ib_Requirement = {0, {0, 1, 2}}; // no hit on 3 ITS ib layers. diff --git a/PWGEM/PhotonMeson/Core/V0PhotonCut.h b/PWGEM/PhotonMeson/Core/V0PhotonCut.h index dcd9cd2c5ce..4c3033ba632 100644 --- a/PWGEM/PhotonMeson/Core/V0PhotonCut.h +++ b/PWGEM/PhotonMeson/Core/V0PhotonCut.h @@ -9,25 +9,29 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -// -// Class for v0 photon selection -// +/// \file V0PhotonCut.h +/// \brief Header of class for V0 photon selection. +/// \author D. Sekihata, daiki.sekihata@cern.ch #ifndef PWGEM_PHOTONMESON_CORE_V0PHOTONCUT_H_ #define PWGEM_PHOTONMESON_CORE_V0PHOTONCUT_H_ -#include "Rtypes.h" - #include "PWGEM/PhotonMeson/Utils/TrackSelection.h" -#include "TMath.h" -#include "TNamed.h" +#include + +#include +#include + +#include #include +#include +#include +#include #include -#include #include -#include + using namespace o2::pwgem::photonmeson; class V0PhotonCut : public TNamed @@ -74,7 +78,7 @@ class V0PhotonCut : public TNamed kNCuts }; - template + template bool IsSelected(TV0 const& v0) const { if (!IsSelectedV0(v0, V0PhotonCuts::kV0PtRange)) { @@ -241,7 +245,7 @@ class V0PhotonCut : public TNamed return true; } - template + template bool IsSelectedV0(T const& v0, const V0PhotonCuts& cut) const { switch (cut) { @@ -494,7 +498,7 @@ class V0PhotonCut : public TNamed bool mDisableITSonly{false}; bool mDisableTPConly{false}; - ClassDef(V0PhotonCut, 1); + ClassDef(V0PhotonCut, 2); }; #endif // PWGEM_PHOTONMESON_CORE_V0PHOTONCUT_H_ diff --git a/PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h b/PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h new file mode 100644 index 00000000000..15e66f8bd96 --- /dev/null +++ b/PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h @@ -0,0 +1,133 @@ +// 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 GammaTablesRedux.h +/// \brief This header provides the table definitions to store very lightweight EMCal clusters per collision +/// \author Marvin Hemmer (marvin.hemmer@cern.ch) - Goethe University Frankfurt +/// + +#ifndef PWGEM_PHOTONMESON_DATAMODEL_GAMMATABLESREDUX_H_ +#define PWGEM_PHOTONMESON_DATAMODEL_GAMMATABLESREDUX_H_ + +#include "PWGEM/PhotonMeson/DataModel/gammaTables.h" + +#include +#include + +#include + +#include +#include +#include + +namespace o2::aod +{ + +namespace emcdownscaling +{ +enum Observable { + kDefinition, + kEnergy, + kEta, + kPhi, + kNCellsExo, + kM02, + kTime, + kDeltaEta, + kDeltaPhi, + nObservables +}; + +// Values in tables are stored in downscaled format to save disk space +const float downscalingFactors[nObservables]{ + 1E0, // Cluster definition + 1E3, // Cluster energy + 1E4, // Cluster eta + 1E4, // Cluster phi + 1E0, // Number of cells + exotic + 1E4, // M02 + 1E2, // Cluster time + 1E5, // diff between cluster and track in eta + 1E5, // diff between cluster and track in phi +}; + +/// \brief convert values for storage into correspoding smaller types +/// \param valueIn input value that should be stored +/// \param observable type of observable that should be stored to determine downscaling +template +OutputType convertForStorage(InputType const& valueIn, Observable observable) +{ + InputType valueToBeChecked = std::lround(valueIn * downscalingFactors[observable]); + if (valueToBeChecked < std::numeric_limits::lowest()) { + LOG(warning) << "Value " << valueToBeChecked << " of observable " << observable << " below lowest possible value of " << typeid(OutputType).name() << ": " << static_cast(std::numeric_limits::lowest()); + valueToBeChecked = static_cast(std::numeric_limits::lowest()); + } + if (valueToBeChecked > std::numeric_limits::max()) { + LOG(warning) << "Value " << valueToBeChecked << " of observable " << observable << " obove highest possible value of " << typeid(OutputType).name() << ": " << static_cast(std::numeric_limits::max()); + valueToBeChecked = static_cast(std::numeric_limits::max()); + } + return static_cast(valueToBeChecked); +} +} // namespace emcdownscaling + +namespace mincluster +{ +DECLARE_SOA_COLUMN(StoredDefinition, storedDefinition, int8_t); //! cluster definition, see EMCALClusterDefinition.h +DECLARE_SOA_COLUMN(StoredE, storedE, uint16_t); //! cluster energy (1 MeV -> Maximum cluster energy of ~65 GeV) +DECLARE_SOA_COLUMN(StoredP, storedP, uint16_t); //! track momentum (1 MeV -> Maximum track momentum of ~65 GeV) +DECLARE_SOA_COLUMN(StoredPt, storedPt, uint16_t); //! track pT (1 MeV -> Maximum track pT of ~65 GeV) +DECLARE_SOA_COLUMN(StoredEta, storedEta, int16_t); //! cluster pseudorapidity (x10,000) +DECLARE_SOA_COLUMN(StoredPhi, storedPhi, uint16_t); //! cluster azimuthal angle (x10 000) from 0 to 2pi +DECLARE_SOA_COLUMN(StoredNCells, storedNCells, uint8_t); //! number of cells in cluster last 7 digits, first digit for exotic flag +DECLARE_SOA_COLUMN(StoredM02, storedM02, int16_t); //! shower shape long axis (x10 000) +DECLARE_SOA_COLUMN(StoredTime, storedTime, int16_t); //! cluster time (10 ps resolution) + +DECLARE_SOA_COLUMN(StoredDeltaEta, storedDeltaEta, int16_t); //! cluster track difference pseudorapidity (x100,000) +DECLARE_SOA_COLUMN(StoredDeltaPhi, storedDeltaPhi, int16_t); //! cluster track difference azimuthal angle (x100,000) + +DECLARE_SOA_DYNAMIC_COLUMN(Definition, definition, [](int8_t storedDefinition) -> int8_t { return storedDefinition; }); //! cluster definition, see EMCALClusterDefinition.h +DECLARE_SOA_DYNAMIC_COLUMN(E, e, [](uint16_t storedE) -> float { return storedE / emcdownscaling::downscalingFactors[emcdownscaling::kEnergy]; }); //! cluster energy (GeV) +DECLARE_SOA_DYNAMIC_COLUMN(Eta, eta, [](int16_t storedEta) -> float { return storedEta / emcdownscaling::downscalingFactors[emcdownscaling::kEta]; }); //! cluster pseudorapidity +DECLARE_SOA_DYNAMIC_COLUMN(Phi, phi, [](uint16_t storedPhi) -> float { return storedPhi / emcdownscaling::downscalingFactors[emcdownscaling::kPhi]; }); //! cluster azimuthal angle (0 to 2pi) +DECLARE_SOA_DYNAMIC_COLUMN(NCells, nCells, [](uint8_t storedNCells) -> uint8_t { return storedNCells & 0x7F; }); //! number of cells in cluster +DECLARE_SOA_DYNAMIC_COLUMN(M02, m02, [](int16_t storedM02) -> float { return storedM02 / emcdownscaling::downscalingFactors[emcdownscaling::kM02]; }); //! shower shape long axis +DECLARE_SOA_DYNAMIC_COLUMN(Time, time, [](int16_t storedTime) -> float { return storedTime / emcdownscaling::downscalingFactors[emcdownscaling::kTime]; }); //! cluster time (ns) +DECLARE_SOA_DYNAMIC_COLUMN(IsExotic, isExotic, [](int16_t storedNCells) -> bool { return static_cast(storedNCells & 0x80); }); //! flag to mark cluster as exotic +DECLARE_SOA_DYNAMIC_COLUMN(TrackP, trackP, [](uint16_t storedP) -> float { return storedP / emcdownscaling::downscalingFactors[emcdownscaling::kEnergy]; }); //! track momentum (GeV) +DECLARE_SOA_DYNAMIC_COLUMN(TrackPt, trackPt, [](uint16_t storedPt) -> float { return storedPt / emcdownscaling::downscalingFactors[emcdownscaling::kEnergy]; }); //! track pT (GeV) + +DECLARE_SOA_DYNAMIC_COLUMN(DeltaEta, deltaEta, [](int16_t storedDeltaEta) -> float { return storedDeltaEta / emcdownscaling::downscalingFactors[emcdownscaling::kDeltaEta]; }); //! cluster track difference pseudorapidity +DECLARE_SOA_DYNAMIC_COLUMN(DeltaPhi, deltaPhi, [](int16_t storedDeltaPhi) -> float { return storedDeltaPhi / emcdownscaling::downscalingFactors[emcdownscaling::kDeltaPhi]; }); //! cluster track difference azimuthal angle + +DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, [](float storedE, float storedEta) -> float { return storedE / emcdownscaling::downscalingFactors[emcdownscaling::kEnergy] / std::cosh(storedEta / emcdownscaling::downscalingFactors[emcdownscaling::kEta]); }); //! cluster pt, assuming m=0 (photons) +} // namespace mincluster + +DECLARE_SOA_TABLE(MinClusters, "AOD", "MINCLUSTER", //! table of skimmed EMCal clusters + o2::soa::Index<>, skimmedcluster::CollisionId, mincluster::StoredDefinition, mincluster::StoredE, mincluster::StoredEta, mincluster::StoredPhi, mincluster::StoredNCells, mincluster::StoredM02, mincluster::StoredTime, + mincluster::Definition, mincluster::E, mincluster::Eta, mincluster::Phi, mincluster::NCells, mincluster::M02, mincluster::Time, mincluster::IsExotic, + mincluster::Pt); + +namespace mintm +{ +DECLARE_SOA_INDEX_COLUMN(MinCluster, minCluster); //! +} // namespace mintm + +DECLARE_SOA_TABLE(MinMTracks, "AOD", "MINMTRACK", //! + mintm::MinClusterId, mincluster::StoredDeltaPhi, mincluster::StoredDeltaEta, mincluster::StoredP, mincluster::StoredPt, + mincluster::DeltaPhi, mincluster::DeltaEta, mincluster::TrackP, mincluster::TrackPt); + +DECLARE_SOA_TABLE(MinMSTracks, "AOD", "MINMSTRACK", //! + mintm::MinClusterId, mincluster::StoredDeltaPhi, mincluster::StoredDeltaEta, mincluster::StoredP, mincluster::StoredPt, + mincluster::DeltaPhi, mincluster::DeltaEta, mincluster::TrackP, mincluster::TrackPt); + +} // namespace o2::aod + +#endif // PWGEM_PHOTONMESON_DATAMODEL_GAMMATABLESREDUX_H_ diff --git a/PWGEM/PhotonMeson/DataModel/gammaTables.h b/PWGEM/PhotonMeson/DataModel/gammaTables.h index a151e3762c6..5fee7952cfd 100644 --- a/PWGEM/PhotonMeson/DataModel/gammaTables.h +++ b/PWGEM/PhotonMeson/DataModel/gammaTables.h @@ -9,6 +9,10 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +/// \file gammaTables.h +/// \brief This header provides the table definitions for the soft photon and neutral meson tables +/// \author D. Sekihata, daiki.sekihata@cern.ch; Marvin Hemmer (marvin.hemmer@cern.ch) - Goethe University Frankfurt + #include "PWGEM/Dilepton/DataModel/dileptonTables.h" #include "Common/Core/RecoDecay.h" @@ -17,6 +21,7 @@ #include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" +#include #include #include @@ -128,7 +133,7 @@ DECLARE_SOA_DYNAMIC_COLUMN(P, p, [](float px, float py, float pz) -> float { ret DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, [](float px, float py) -> float { return RecoDecay::sqrtSumOfSquares(px, py); }); DECLARE_SOA_DYNAMIC_COLUMN(Eta, eta, [](float px, float py, float pz) -> float { return RecoDecay::eta(std::array{px, py, pz}); }); DECLARE_SOA_DYNAMIC_COLUMN(Phi, phi, [](float px, float py) -> float { return RecoDecay::phi(px, py); }); -DECLARE_SOA_DYNAMIC_COLUMN(Tgl, tgl, [](float px, float py, float pz) -> float { return std::tan(M_PI_2 - 2 * std::atan(std::exp(-RecoDecay::eta(std::array{px, py, pz})))); }); +DECLARE_SOA_DYNAMIC_COLUMN(Tgl, tgl, [](float px, float py, float pz) -> float { return std::tan(o2::constants::math::PIHalf - 2 * std::atan(std::exp(-RecoDecay::eta(std::array{px, py, pz})))); }); DECLARE_SOA_DYNAMIC_COLUMN(MeanClusterSizeITS, meanClusterSizeITS, [](uint32_t itsClusterSizes) -> float { int total_cluster_size = 0, nl = 0; for (unsigned int layer = 0; layer < 7; layer++) { @@ -566,22 +571,30 @@ DECLARE_SOA_COLUMN(NLM, nlm, int); //! numbe namespace emccluster { -DECLARE_SOA_INDEX_COLUMN(EMEvent, emevent); //! -DECLARE_SOA_COLUMN(CoreEnergy, coreEnergy, float); //! cluster core energy (GeV) -DECLARE_SOA_COLUMN(Time, time, float); //! cluster time (ns) -DECLARE_SOA_COLUMN(IsExotic, isExotic, bool); //! flag to mark cluster as exotic -DECLARE_SOA_COLUMN(Definition, definition, int); //! cluster definition, see EMCALClusterDefinition.h -DECLARE_SOA_ARRAY_INDEX_COLUMN(Track, track); //! TrackIds -DECLARE_SOA_COLUMN(DeltaPhi, deltaPhi, std::vector); //! phi values of the matched tracks -DECLARE_SOA_COLUMN(DeltaEta, deltaEta, std::vector); //! eta values of the matched tracks -DECLARE_SOA_COLUMN(TrackP, trackp, std::vector); //! momentum values of the matched tracks -DECLARE_SOA_COLUMN(TrackPt, trackpt, std::vector); //! pt values of the matched tracks -DECLARE_SOA_COLUMN(DeltaPhiSec, deltaPhiSec, std::vector); //! phi values of the matched secondary tracks -DECLARE_SOA_COLUMN(DeltaEtaSec, deltaEtaSec, std::vector); //! eta values of the matched secondary tracks -DECLARE_SOA_COLUMN(TrackPSec, trackpSec, std::vector); //! momentum values of the matched secondary tracks -DECLARE_SOA_COLUMN(TrackPtSec, trackptSec, std::vector); //! pt values of the matched secondary tracks -DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, [](float e, float eta, float m = 0) -> float { return sqrt(e * e - m * m) / cosh(eta); }); //! cluster pt, mass to be given as argument when getter is called! +DECLARE_SOA_INDEX_COLUMN(EMEvent, emevent); //! +DECLARE_SOA_COLUMN(CoreEnergy, coreEnergy, float); //! cluster core energy (GeV) +DECLARE_SOA_COLUMN(Time, time, float); //! cluster time (ns) +DECLARE_SOA_COLUMN(IsExotic, isExotic, bool); //! flag to mark cluster as exotic +DECLARE_SOA_COLUMN(Definition, definition, int); //! cluster definition, see EMCALClusterDefinition.h +DECLARE_SOA_ARRAY_INDEX_COLUMN(Track, track); //! TrackIds +DECLARE_SOA_COLUMN(DeltaPhi, deltaPhi, std::vector); //! phi values of the matched tracks +DECLARE_SOA_COLUMN(DeltaEta, deltaEta, std::vector); //! eta values of the matched tracks +DECLARE_SOA_COLUMN(TrackP, trackp, std::vector); //! momentum values of the matched tracks +DECLARE_SOA_COLUMN(TrackPt, trackpt, std::vector); //! pt values of the matched tracks +DECLARE_SOA_COLUMN(DeltaPhiSec, deltaPhiSec, std::vector); //! phi values of the matched secondary tracks +DECLARE_SOA_COLUMN(DeltaEtaSec, deltaEtaSec, std::vector); //! eta values of the matched secondary tracks +DECLARE_SOA_COLUMN(TrackPSec, trackpSec, std::vector); //! momentum values of the matched secondary tracks +DECLARE_SOA_COLUMN(TrackPtSec, trackptSec, std::vector); //! pt values of the matched secondary tracks +DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, [](float e, float eta, float m = 0) -> float { return std::sqrt(e * e - m * m) / std::cosh(eta); }); //! cluster pt, mass to be given as argument when getter is called! } // namespace emccluster + +namespace emctm +{ +DECLARE_SOA_COLUMN(DeltaPhi, deltaPhi, float); //! phi values of the matched tracks +DECLARE_SOA_COLUMN(DeltaEta, deltaEta, float); //! eta values of the matched tracks +DECLARE_SOA_COLUMN(TrackP, trackp, float); //! momentum values of the matched tracks +DECLARE_SOA_COLUMN(TrackPt, trackpt, float); //! pt values of the matched tracks +} // namespace emctm DECLARE_SOA_TABLE_VERSIONED(SkimEMCClusters_000, "AOD", "SKIMEMCCLUSTER", 0, //! table of skimmed EMCal clusters o2::soa::Index<>, skimmedcluster::CollisionId, emccluster::Definition, skimmedcluster::E, skimmedcluster::Eta, skimmedcluster::Phi, skimmedcluster::M02, skimmedcluster::NCells, skimmedcluster::Time, emccluster::IsExotic, emccluster::DeltaPhi, @@ -597,23 +610,23 @@ DECLARE_SOA_TABLE_VERSIONED(SkimEMCClusters_001, "AOD", "SKIMEMCCLUSTER", 1, //! using SkimEMCClusters = SkimEMCClusters_001; using SkimEMCCluster = SkimEMCClusters_001::iterator; -// DECLARE_SOA_TABLE_VERSIONED(EmEmcClusters_000, "AOD", "EMEMCCLUSTER", 0, //! table of skimmed EMCal clusters -// o2::soa::Index<>, skimmedcluster::CollisionId, emccluster::Definition, skimmedcluster::E, skimmedcluster::Eta, skimmedcluster::Phi, -// skimmedcluster::M02, skimmedcluster::NCells, skimmedcluster::Time, emccluster::IsExotic, emccluster::Pt); +DECLARE_SOA_TABLE_VERSIONED(EmEmcClusters_000, "AOD", "EMEMCCLUSTER", 0, //! table of skimmed EMCal clusters + o2::soa::Index<>, skimmedcluster::CollisionId, emccluster::Definition, skimmedcluster::E, skimmedcluster::Eta, skimmedcluster::Phi, + skimmedcluster::M02, skimmedcluster::NCells, skimmedcluster::Time, emccluster::IsExotic, emccluster::Pt); -// using EmEmcClusters = EmEmcClusters_000; -// using EmEmcCluster = EmEmcClusters_000::iterator; +using EmEmcClusters = EmEmcClusters_000; +using EmEmcCluster = EmEmcClusters_000::iterator; -// namespace trackmatching -// { -// DECLARE_SOA_INDEX_COLUMN(EmEmcCluster, emEmcCluster); //! -// } // namespace trackmatching +namespace trackmatching +{ +DECLARE_SOA_INDEX_COLUMN(EmEmcCluster, emEmcCluster); //! +} // namespace trackmatching -// DECLARE_SOA_TABLE(EmEmcMTracks, "AOD", "EMEMCMTRACK", //! -// o2::soa::Index<>, trackmatching::EmEmcClusterId, emccluster::DeltaPhi, emccluster::DeltaEta, emccluster::TrackP, emccluster::TrackPt); +DECLARE_SOA_TABLE(EmEmcMTracks, "AOD", "EMEMCMTRACK", //! + trackmatching::EmEmcClusterId, emctm::DeltaPhi, emctm::DeltaEta, emctm::TrackP, emctm::TrackPt); -// DECLARE_SOA_TABLE(EmEmcMSTracks, "AOD", "EMEMCMSTRACK", //! -// o2::soa::Index<>, trackmatching::EmEmcClusterId, emccluster::DeltaPhiSec, emccluster::DeltaEtaSec, emccluster::TrackPSec, emccluster::TrackPtSec); +DECLARE_SOA_TABLE(EmEmcMSTracks, "AOD", "EMEMCMSTRACK", //! + trackmatching::EmEmcClusterId, emctm::DeltaPhi, emctm::DeltaEta, emctm::TrackP, emctm::TrackPt); DECLARE_SOA_TABLE(EMCEMEventIds, "AOD", "EMCEMEVENTID", emccluster::EMEventId); // To be joined with SkimEMCClusters table at analysis level. // iterators @@ -632,10 +645,10 @@ DECLARE_SOA_COLUMN(CellZ, cellz, int); // DECLARE_SOA_COLUMN(TrackPhi, trackphi, float); //! phi of the matched track // DECLARE_SOA_COLUMN(TrackP, trackp, float); //! momentum of the matched track // DECLARE_SOA_COLUMN(TrackPt, trackpt, float); //! pt of the matched track -DECLARE_SOA_DYNAMIC_COLUMN(Px, px, [](float e, float x, float y, float z, float m = 0) -> float { return x / RecoDecay::sqrtSumOfSquares(x, y, z) * sqrt(e * e - m * m); }); -DECLARE_SOA_DYNAMIC_COLUMN(Py, py, [](float e, float x, float y, float z, float m = 0) -> float { return y / RecoDecay::sqrtSumOfSquares(x, y, z) * sqrt(e * e - m * m); }); -DECLARE_SOA_DYNAMIC_COLUMN(Pz, pz, [](float e, float x, float y, float z, float m = 0) -> float { return z / RecoDecay::sqrtSumOfSquares(x, y, z) * sqrt(e * e - m * m); }); -DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, [](float e, float x, float y, float z, float m = 0) -> float { return RecoDecay::sqrtSumOfSquares(x, y) / RecoDecay::sqrtSumOfSquares(x, y, z) * sqrt(e * e - m * m); }); +DECLARE_SOA_DYNAMIC_COLUMN(Px, px, [](float e, float x, float y, float z, float m = 0) -> float { return x / RecoDecay::sqrtSumOfSquares(x, y, z) * std::sqrt(e * e - m * m); }); +DECLARE_SOA_DYNAMIC_COLUMN(Py, py, [](float e, float x, float y, float z, float m = 0) -> float { return y / RecoDecay::sqrtSumOfSquares(x, y, z) * std::sqrt(e * e - m * m); }); +DECLARE_SOA_DYNAMIC_COLUMN(Pz, pz, [](float e, float x, float y, float z, float m = 0) -> float { return z / RecoDecay::sqrtSumOfSquares(x, y, z) * std::sqrt(e * e - m * m); }); +DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, [](float e, float x, float y, float z, float m = 0) -> float { return RecoDecay::sqrtSumOfSquares(x, y) / RecoDecay::sqrtSumOfSquares(x, y, z) * std::sqrt(e * e - m * m); }); DECLARE_SOA_DYNAMIC_COLUMN(Eta, eta, [](float x, float y, float z) -> float { return RecoDecay::eta(std::array{x, y, z}); }); DECLARE_SOA_DYNAMIC_COLUMN(Phi, phi, [](float x, float y) -> float { return RecoDecay::phi(x, y); }); } // namespace phoscluster diff --git a/PWGEM/PhotonMeson/TableProducer/createPCM.cxx b/PWGEM/PhotonMeson/TableProducer/createPCM.cxx index 4ca74633e1b..138a190f7cd 100644 --- a/PWGEM/PhotonMeson/TableProducer/createPCM.cxx +++ b/PWGEM/PhotonMeson/TableProducer/createPCM.cxx @@ -8,32 +8,39 @@ // 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. -// -// ======================== -// -// This code produces photon data tables. -// Please write to: daiki.sekihata@cern.ch -#include -#include -#include -#include "Framework/runDataProcessing.h" -#include "Framework/AnalysisTask.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/ASoAHelpers.h" -#include "Common/Core/trackUtilities.h" -#include "Common/Core/RecoDecay.h" -#include "Common/DataModel/CollisionAssociationTables.h" -#include "DCAFitter/DCAFitterN.h" -#include "DetectorsBase/Propagator.h" -#include "DetectorsBase/GeometryManager.h" -#include "DataFormatsParameters/GRPObject.h" -#include "DataFormatsParameters/GRPMagField.h" -#include "CCDB/BasicCCDBManager.h" +/// \file createPCM.cxx +/// \brief This code produces photon data tables. +/// \author Daiki Sekihata , Tokyo + #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" -#include "PWGLF/DataModel/LFStrangenessTables.h" #include "PWGEM/PhotonMeson/Utils/PCMUtilities.h" #include "PWGEM/PhotonMeson/Utils/TrackSelection.h" +// +#include "PWGLF/DataModel/LFStrangenessTables.h" + +#include "Common/Core/RecoDecay.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/CollisionAssociationTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include using namespace o2; using namespace o2::soa; @@ -159,7 +166,7 @@ struct createPCM { d_bz = d_bz_input; fitter.setBz(d_bz); o2::parameters::GRPMagField grpmag; - if (fabs(d_bz) > 1e-5) { + if (std::fabs(d_bz) > 1e-5) { grpmag.setL3Current(30000.f / (d_bz / 5.0f)); } o2::base::Propagator::initFieldFromGRP(&grpmag); @@ -242,7 +249,7 @@ struct createPCM { float xyz[3] = {0.f, 0.f, 0.f}; Vtx_recalculation(o2::base::Propagator::Instance(), pos, ele, xyz, matCorr); - float recalculatedVtxR = std::sqrt(pow(xyz[0], 2) + pow(xyz[1], 2)); + float recalculatedVtxR = std::sqrt(std::pow(xyz[0], 2) + std::pow(xyz[1], 2)); // LOGF(info, "recalculated vtx : x = %f , y = %f , z = %f", xyz[0], xyz[1], xyz[2]); if (recalculatedVtxR > std::min(pos.x(), ele.x()) + margin_r && (pos.x() > 1.f && ele.x() > 1.f)) { return false; @@ -329,10 +336,10 @@ struct createPCM { template bool isSelected(TTrack const& track) { - if (track.pt() < minpt || abs(track.eta()) > maxeta) { + if (track.pt() < minpt || std::abs(track.eta()) > maxeta) { return false; } - if (abs(track.dcaXY()) < dcamin || dcamax < abs(track.dcaXY())) { + if (std::abs(track.dcaXY()) < dcamin || dcamax < std::abs(track.dcaXY())) { return false; } if (!track.hasITS() && !track.hasTPC()) { @@ -357,7 +364,7 @@ struct createPCM { return false; } - if (abs(track.z() / track.x() - track.tgl()) > 0.5) { + if (std::abs(track.z() / track.x() - track.tgl()) > 0.5) { return false; } @@ -402,8 +409,8 @@ struct createPCM { initCCDB(bc); // registry.fill(HIST("hEventCounter"), 1); - int32_t min_sw = std::max(int64_t(0), collision.globalIndex()); - int32_t max_sw = std::min(int64_t(min_sw + nsw), int64_t(collisions.size())); + int32_t min_sw = std::max(static_cast(0), collision.globalIndex()); + int32_t max_sw = std::min(static_cast(min_sw + nsw), static_cast(collisions.size())); // LOGF(info, "orphan_posTracks.size() = %d, orphan_negTracks.size() = %d", orphan_posTracks.size(), orphan_negTracks.size()); negTracks_sw.reserve(max_sw - min_sw); @@ -419,9 +426,9 @@ struct createPCM { } // LOGF(info, "min_sw = %d , max_sw = %d , collision.globalIndex() = %d , n posTracks_sw = %d , n negTracks_sw = %d", min_sw, max_sw, collision.globalIndex(), npos, nneg); - for (auto& negTracks_coll : negTracks_sw) { - for (auto& posTracks_coll : posTracks_sw) { - for (auto& [ele, pos] : combinations(CombinationsFullIndexPolicy(negTracks_coll, posTracks_coll))) { + for (const auto& negTracks_coll : negTracks_sw) { + for (const auto& posTracks_coll : posTracks_sw) { + for (const auto& [ele, pos] : combinations(CombinationsFullIndexPolicy(negTracks_coll, posTracks_coll))) { if (!isSelected(ele) || !isSelected(pos)) { continue; } @@ -443,9 +450,9 @@ struct createPCM { // collision_in_sw.globalIndex(), ele.collisionId(), pos.collisionId(), ele.globalIndex(), pos.globalIndex()); fillV0Table(collision_in_sw, ele, pos, false); } // end of searching window loop - } // end of pairing loop - } // end of pos track loop in sw - } // end of pos track loop in sw + } // end of pairing loop + } // end of pos track loop in sw + } // end of pos track loop in sw // LOGF(info, "possible number of V0 = %d", cospa_map.size()); std::map, bool> used_pair_map; @@ -522,7 +529,7 @@ struct createPCM { Preslice trackIndicesPerCollision = aod::track_association::collisionId; void processTrkCollAsso(aod::TrackAssoc const& trackIndices, FullTracksExtIU const&, aod::Collisions const& collisions, aod::BCsWithTimestamps const&) { - for (auto& collision : collisions) { + for (const auto& collision : collisions) { registry.fill(HIST("hEventCounter"), 1); auto bc = collision.bc_as(); @@ -530,7 +537,7 @@ struct createPCM { auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, collision.globalIndex()); // LOGF(info,"%d tracks in collision %d", trackIdsThisCollision.size(), collision.globalIndex()); - for (auto& [eleId, posId] : combinations(CombinationsStrictlyUpperIndexPolicy(trackIdsThisCollision, trackIdsThisCollision))) { + for (const auto& [eleId, posId] : combinations(CombinationsStrictlyUpperIndexPolicy(trackIdsThisCollision, trackIdsThisCollision))) { auto ele = eleId.track_as(); auto pos = posId.track_as(); // LOGF(info,"eleId = %d , posId = %d", ele.globalIndex(), pos.globalIndex()); @@ -550,7 +557,7 @@ struct createPCM { } } } // end of collision loop - } // end of process + } // end of process PROCESS_SWITCH(createPCM, processTrkCollAsso, "create V0s with track-to-collision associator", false); }; diff --git a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx index f84f260bc46..5240f826855 100644 --- a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx +++ b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx @@ -56,7 +56,6 @@ #include #include #include -#include #include #include #include diff --git a/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx b/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx index b7cefc4045e..e4f11543583 100644 --- a/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx +++ b/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx @@ -14,6 +14,7 @@ /// \author marvin.hemmer@cern.ch /// dependencies: emcal-correction-task +#include "PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" #include "PWGEM/PhotonMeson/Utils/emcalHistoDefinitions.h" #include "PWGJE/DataModel/EMCALClusters.h" @@ -37,11 +38,13 @@ #include #include #include +#include #include using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; +using namespace o2::aod::emcdownscaling; struct SkimmerGammaCalo { @@ -53,12 +56,21 @@ struct SkimmerGammaCalo { Produces tableEMCClusterMCLabels; Produces tableCellEMCReco; + Produces tableEmEmcClusters; + Produces tableEmEmcMTracks; + Produces tableEmEmcMSTracks; + + Produces tableMinClusters; + Produces tableMinMTracks; + Produces tableMinMSTracks; + // Configurable for filter/cuts Configurable minTime{"minTime", -200., "Minimum cluster time for time cut"}; Configurable maxTime{"maxTime", +200., "Maximum cluster time for time cut"}; Configurable minM02{"minM02", 0.0, "Minimum M02 for M02 cut"}; Configurable maxM02{"maxM02", 1.0, "Maximum M02 for M02 cut"}; Configurable minE{"minE", 0.5, "Minimum energy for energy cut"}; + Configurable maxE{"maxE", std::numeric_limits::max(), "Maximum energy for energy cut"}; Configurable removeExotic{"removeExotic", false, "Flag to enable the removal of exotic clusters."}; Configurable> clusterDefinitions{"clusterDefinitions", {0, 1, 2, 10, 11, 12, 13, 20, 21, 22, 30, 40, 41, 42, 43, 44, 45}, "Cluster definitions to be accepted (e.g. 13 for kV3MostSplitLowSeed)"}; Configurable maxdEta{"maxdEta", 0.1, "Set a maximum difference in eta for tracks and cluster to still count as matched"}; @@ -163,7 +175,7 @@ struct SkimmerGammaCalo { continue; } // Energy cut - if (emccluster.energy() < minE) { + if (emccluster.energy() < minE || emccluster.energy() > maxE) { historeg.fill(HIST("hCaloClusterFilter"), 2); continue; } @@ -234,6 +246,39 @@ struct SkimmerGammaCalo { tableGammaEMCReco(emccluster.collisionId(), emccluster.definition(), emccluster.energy(), emccluster.eta(), emccluster.phi(), emccluster.m02(), emccluster.nCells(), emccluster.time(), emccluster.isExotic(), vPhi, vEta, vP, vPt, vPhiSecondaries, vEtaSecondaries, vPSecondaries, vPtSecondaries); + + tableEmEmcClusters(emccluster.collisionId(), emccluster.definition(), emccluster.energy(), emccluster.eta(), emccluster.phi(), emccluster.m02(), + emccluster.nCells(), emccluster.time(), emccluster.isExotic()); + + tableMinClusters(emccluster.collisionId(), convertForStorage(emccluster.definition(), Observable::kDefinition), + convertForStorage(emccluster.energy(), Observable::kEnergy), + convertForStorage(emccluster.eta(), Observable::kEta), + convertForStorage(emccluster.phi(), Observable::kPhi), + convertForStorage((emccluster.nCells() & 0x7F) | (emccluster.isExotic() << 7), Observable::kNCellsExo), + convertForStorage(emccluster.m02(), Observable::kM02), + convertForStorage(emccluster.time(), Observable::kTime)); + + if (vEta.size() > 0) { + for (size_t iPart = 0; iPart < vEta.size(); ++iPart) { + tableEmEmcMTracks(tableEmEmcClusters.lastIndex(), vEta[iPart], vPhi[iPart], vP[iPart], vPt[iPart]); + tableMinMTracks(tableMinClusters.lastIndex(), + convertForStorage(vPhi[iPart], Observable::kDeltaPhi), + convertForStorage(vEta[iPart], Observable::kDeltaEta), + convertForStorage(vP[iPart], Observable::kEnergy), + convertForStorage(vPt[iPart], Observable::kEnergy)); + } + } + if (vEtaSecondaries.size() > 0) { + for (size_t iPart = 0; iPart < vEtaSecondaries.size(); ++iPart) { + tableEmEmcMSTracks(tableEmEmcClusters.lastIndex(), vEtaSecondaries[iPart], vPhiSecondaries[iPart], vPSecondaries[iPart], vPtSecondaries[iPart]); + tableMinMSTracks(tableMinClusters.lastIndex(), + convertForStorage(vPhiSecondaries[iPart], Observable::kDeltaPhi), + convertForStorage(vEtaSecondaries[iPart], Observable::kDeltaEta), + convertForStorage(vPSecondaries[iPart], Observable::kEnergy), + convertForStorage(vPtSecondaries[iPart], Observable::kEnergy)); + } + } + vEta.clear(); vPhi.clear(); vP.clear(); diff --git a/PWGEM/PhotonMeson/Tasks/CMakeLists.txt b/PWGEM/PhotonMeson/Tasks/CMakeLists.txt index 0e727db9e0c..ef6371b5e94 100644 --- a/PWGEM/PhotonMeson/Tasks/CMakeLists.txt +++ b/PWGEM/PhotonMeson/Tasks/CMakeLists.txt @@ -101,6 +101,11 @@ o2physics_add_dpl_workflow(pi0eta-to-gammagamma-emcemc PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(pi0eta-to-gammagamma-pcmemc + SOURCES Pi0EtaToGammaGammaPCMEMC.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(pi0eta-to-gammagamma-mc-pcmpcm SOURCES Pi0EtaToGammaGammaMCPCMPCM.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore O2Physics::MLCore @@ -181,3 +186,8 @@ o2physics_add_dpl_workflow(calib-task-emc SOURCES calibTaskEmc.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::EMCALBase O2::EMCALCalib O2Physics::PWGEMPhotonMesonCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(test-task-emc + SOURCES testTaskEmc.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::EMCALBase O2::EMCALCalib O2Physics::PWGEMPhotonMesonCore + COMPONENT_NAME Analysis) diff --git a/PWGEM/PhotonMeson/Tasks/MaterialBudget.cxx b/PWGEM/PhotonMeson/Tasks/MaterialBudget.cxx index e648dfe8fca..298f05a1198 100644 --- a/PWGEM/PhotonMeson/Tasks/MaterialBudget.cxx +++ b/PWGEM/PhotonMeson/Tasks/MaterialBudget.cxx @@ -14,46 +14,48 @@ /// \author S. Mrozinski, smrozins@cern.ch #include "PWGEM/Dilepton/Utils/MCUtilities.h" -#include "PWGEM/PhotonMeson/Core/CutsLibrary.h" +#include "PWGEM/PhotonMeson/Core/DalitzEECut.h" #include "PWGEM/PhotonMeson/Core/EMPhotonEventCut.h" -#include "PWGEM/PhotonMeson/Core/HistogramsLibrary.h" -#include "PWGEM/PhotonMeson/Core/PairCut.h" #include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" -#include "PWGEM/PhotonMeson/Utils/MCUtilities.h" -#include "PWGEM/PhotonMeson/Utils/PairUtilities.h" -#include "Common/Core/RecoDecay.h" -#include "Common/Core/TPCVDriftManager.h" -#include "Common/Core/TrackSelection.h" -#include "Common/Core/trackUtilities.h" +#include "Common/CCDB/EventSelectionParams.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/McCollisionExtra.h" -#include "Common/DataModel/TrackSelectionTables.h" #include +#include #include #include #include #include +#include #include #include +#include +#include #include +#include // IWYU pragma: keep +#include +#include +#include +#include #include #include -#include +#include +#include #include +#include +#include +#include +#include #include -#include #include -#include #include +#include #include -#include -#include #include #include @@ -62,7 +64,6 @@ using namespace o2::aod; using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::soa; -using namespace o2::aod::pwgem::photonmeson::photonpair; using namespace o2::aod::pwgem::photon; using namespace o2::aod::pwgem::dilepton::utils::emtrackutil; using namespace o2::aod::pwgem::dilepton::utils::mcutil; @@ -773,7 +774,7 @@ struct MaterialBudget { { for (const auto& prev : pool) { for (const auto& g1 : current) { - if (!fV0PhotonCut.IsSelected(g1)) { + if (!fV0PhotonCut.IsSelected(g1)) { continue; } for (const auto& g2 : prev) { @@ -796,7 +797,7 @@ struct MaterialBudget { { for (const auto& g1 : v0s) { - if (!fV0PhotonCut.IsSelected(g1)) { + if (!fV0PhotonCut.IsSelected(g1)) { continue; } ROOT::Math::PtEtaPhiMVector vGamma(g1.pt(), g1.eta(), g1.phi(), 0.); @@ -862,7 +863,7 @@ struct MaterialBudget { for (const auto& v0 : v0sThisCollision) { - if (!fV0PhotonCut.IsSelected(v0)) { + if (!fV0PhotonCut.IsSelected(v0)) { continue; } @@ -963,7 +964,7 @@ struct MaterialBudget { std::vector eventCopy; eventCopy.reserve(v0sThisCollision.size()); for (const auto& v0 : v0sThisCollision) { - if (!fV0PhotonCut.IsSelected(v0)) { + if (!fV0PhotonCut.IsSelected(v0)) { continue; } eventCopy.push_back({v0.pt(), v0.eta(), v0.phi(), v0.vz(), v0.v0radius()}); @@ -1045,7 +1046,7 @@ struct MaterialBudget { auto posmc = pos.template emmcparticle_as(); auto elemc = ele.template emmcparticle_as(); - if (!fV0PhotonCut.IsSelected(v0)) { + if (!fV0PhotonCut.IsSelected(v0)) { continue; } @@ -1133,7 +1134,7 @@ struct MaterialBudget { for (auto const& v0 : v0s) { - if (!fV0PhotonCut.IsSelected(v0)) { + if (!fV0PhotonCut.IsSelected(v0)) { continue; } diff --git a/PWGEM/PhotonMeson/Tasks/MaterialBudgetMC.cxx b/PWGEM/PhotonMeson/Tasks/MaterialBudgetMC.cxx index 2c0a33d32ec..a2cea8bd2b3 100644 --- a/PWGEM/PhotonMeson/Tasks/MaterialBudgetMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/MaterialBudgetMC.cxx @@ -8,39 +8,43 @@ // 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. -// -// ======================== -// -// This code loops over v0 photons for studying material budget. -// Please write to: daiki.sekihata@cern.ch -#include "PWGEM/Dilepton/Utils/MCUtilities.h" +/// \file MaterialBudgetMC.cxx +/// \brief Task to analyse and calculate the material budget weights in MC +/// \author D. Sekihata: daiki.sekihata@cern, S. Mrozinski: smrozins@cern.ch + #include "PWGEM/PhotonMeson/Core/CutsLibrary.h" +#include "PWGEM/PhotonMeson/Core/EMPhotonEventCut.h" #include "PWGEM/PhotonMeson/Core/HistogramsLibrary.h" #include "PWGEM/PhotonMeson/Core/PairCut.h" #include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" #include "PWGEM/PhotonMeson/Utils/MCUtilities.h" #include "PWGEM/PhotonMeson/Utils/PairUtilities.h" +// +#include "PWGEM/Dilepton/Utils/MCUtilities.h" -#include "Common/Core/RecoDecay.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/TrackSelectionTables.h" - -#include "Framework/ASoAHelpers.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" -#include "ReconstructionDataFormats/Track.h" - -#include "Math/Vector4D.h" -#include "TString.h" -#include -#include +#include +#include +#include +#include +#include +#include + +#include // IWYU pragma: keep +#include +#include +#include +#include + +#include +#include #include #include +#include #include using namespace o2; @@ -112,8 +116,8 @@ struct MaterialBudgetMC { template void add_pair_histograms(THashList* list_pair, const std::string pairname, TCuts1 const& tagcuts, TCuts2 const& probecuts, TCuts3 const& cuts3) { - for (auto& tagcut : tagcuts) { - for (auto& probecut : probecuts) { + for (const auto& tagcut : tagcuts) { + for (const auto& probecut : probecuts) { std::string cutname1 = tagcut.GetName(); std::string cutname2 = probecut.GetName(); @@ -126,7 +130,7 @@ struct MaterialBudgetMC { o2::aod::pwgem::photon::histogram::AddHistClass(list_pair_subsys, photon_cut_name.data()); THashList* list_pair_subsys_photoncut = reinterpret_cast(list_pair_subsys->FindObject(photon_cut_name.data())); - for (auto& cut3 : cuts3) { + for (const auto& cut3 : cuts3) { std::string pair_cut_name = cut3.GetName(); o2::aod::pwgem::photon::histogram::AddHistClass(list_pair_subsys_photoncut, pair_cut_name.data()); THashList* list_pair_subsys_paircut = reinterpret_cast(list_pair_subsys_photoncut->FindObject(pair_cut_name.data())); @@ -159,7 +163,7 @@ struct MaterialBudgetMC { o2::aod::pwgem::photon::histogram::DefineHistograms(list_v0_cut, "material_budget_study", "V0"); } - for (auto& pairname : fPairNames) { + for (const auto& pairname : fPairNames) { LOGF(info, "Enabled pairs = %s", pairname.data()); THashList* list_ev_pair = reinterpret_cast(o2::aod::pwgem::photon::histogram::AddHistClass(list_ev, pairname.data())); @@ -238,7 +242,7 @@ struct MaterialBudgetMC { THashList* list_ev_pair_after = static_cast(fMainList->FindObject("Event")->FindObject(pairnames[pairtype].data())->FindObject(event_types[1].data())); THashList* list_v0 = static_cast(fMainList->FindObject("V0")); double value[4] = {0.f}; - for (auto& collision : collisions) { + for (const auto& collision : collisions) { float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; if (centralities[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) { continue; @@ -253,10 +257,10 @@ struct MaterialBudgetMC { reinterpret_cast(list_ev_pair_after->FindObject("hCollisionCounter"))->Fill("accepted", 1.f); auto photons_coll = photons.sliceBy(perCollision, collision.globalIndex()); - for (auto& cut : cuts) { - for (auto& photon : photons_coll) { + for (const auto& cut : cuts) { + for (const auto& photon : photons_coll) { - if (!cut.template IsSelected(photon)) { + if (!cut.template IsSelected(photon)) { continue; } @@ -276,11 +280,11 @@ struct MaterialBudgetMC { continue; } - float phi_cp = atan2(photon.vy(), photon.vx()); - float eta_cp = std::atanh(photon.vz() / sqrt(pow(photon.vx(), 2) + pow(photon.vy(), 2) + pow(photon.vz(), 2))); + float phi_cp = std::atan2(photon.vy(), photon.vx()); + float eta_cp = std::atanh(photon.vz() / std::sqrt(std::pow(photon.vx(), 2) + std::pow(photon.vy(), 2) + std::pow(photon.vz(), 2))); value[0] = photon.pt(); value[1] = photon.v0radius(); - value[2] = phi_cp > 0 ? phi_cp : phi_cp + TMath::TwoPi(); + value[2] = RecoDecay::constrainAngle(phi_cp); value[3] = eta_cp; reinterpret_cast(list_v0->FindObject(cut.GetName())->FindObject("hs_conv_point"))->Fill(value); @@ -295,7 +299,7 @@ struct MaterialBudgetMC { { THashList* list_pair_ss = static_cast(fMainList->FindObject("Pair")->FindObject(pairnames[pairtype].data())); - for (auto& collision : collisions) { + for (const auto& collision : collisions) { float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; if (centralities[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) { continue; @@ -309,10 +313,10 @@ struct MaterialBudgetMC { double value[6] = {0.f}; float phi_cp2 = 0.f, eta_cp2 = 0.f; - for (auto& tagcut : tagcuts) { - for (auto& probecut : probecuts) { - for (auto& g1 : photons1_coll) { - for (auto& g2 : photons2_coll) { + for (const auto& tagcut : tagcuts) { + for (const auto& probecut : probecuts) { + for (const auto& g1 : photons1_coll) { + for (const auto& g2 : photons2_coll) { if (g1.globalIndex() == g2.globalIndex()) { continue; } @@ -320,7 +324,7 @@ struct MaterialBudgetMC { continue; } - for (auto& paircut : paircuts) { + for (const auto& paircut : paircuts) { if (!paircut.IsSelected(g1, g2)) { continue; } @@ -367,13 +371,13 @@ struct MaterialBudgetMC { ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); // tag ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); // probe ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; - phi_cp2 = atan2(g2.vy(), g2.vx()); - eta_cp2 = std::atanh(g2.vz() / sqrt(pow(g2.vx(), 2) + pow(g2.vy(), 2) + pow(g2.vz(), 2))); + phi_cp2 = std::atan2(g2.vy(), g2.vx()); + eta_cp2 = std::atanh(g2.vz() / std::sqrt(std::pow(g2.vx(), 2) + std::pow(g2.vy(), 2) + std::pow(g2.vz(), 2))); value[0] = v12.M(); value[1] = g1.pt(); value[2] = g2.pt(); value[3] = g2.v0radius(); - value[4] = phi_cp2 > 0.f ? phi_cp2 : phi_cp2 + TMath::TwoPi(); + value[4] = RecoDecay::constrainAngle(phi_cp2); value[5] = eta_cp2; reinterpret_cast(list_pair_ss->FindObject(Form("%s_%s", tagcut.GetName(), probecut.GetName()))->FindObject(paircut.GetName())->FindObject("hs_conv_point_same"))->Fill(value); } // end of pair cut loop @@ -385,7 +389,7 @@ struct MaterialBudgetMC { } Partition grouped_collisions = cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < cfgCentMax; // this goes to same event. - Filter collisionFilter_common = nabs(o2::aod::collision::posZ) < 10.f && o2::aod::collision::numContrib > (uint16_t)0 && o2::aod::evsel::sel8 == true; + Filter collisionFilter_common = nabs(o2::aod::collision::posZ) < 10.f && o2::aod::collision::numContrib > static_cast(0) && o2::aod::evsel::sel8 == true; Filter collisionFilter_centrality = (cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0A && o2::aod::cent::centFT0A < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < cfgCentMax); using MyFilteredCollisions = soa::Filtered; // this goes to same event pairing. @@ -401,7 +405,7 @@ struct MaterialBudgetMC { // loop over mc stack and fill histograms for pure MC truth signals // all MC tracks which belong to the MC event corresponding to the current reconstructed event - for (auto& collision : grouped_collisions) { + for (const auto& collision : grouped_collisions) { float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; if (centralities[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) { continue; @@ -421,7 +425,7 @@ struct MaterialBudgetMC { } reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(3.0); - if (abs(collision.posZ()) > 10.0) { + if (std::abs(collision.posZ()) > 10.0) { continue; } reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(4.0); diff --git a/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaEMCEMC.cxx b/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaEMCEMC.cxx index b1b657f8a1f..1e0567d6ee6 100644 --- a/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaEMCEMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaEMCEMC.cxx @@ -26,11 +26,11 @@ using namespace o2::aod; using namespace o2::framework; using namespace o2::aod::pwgem::photonmeson::photonpair; -using MyEMCClusters = soa::Join; +using MyEMCClusters = soa::Join; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask>(cfgc, TaskName{"pi0eta-to-gammagamma-emcemc"}), + adaptAnalysisTask>(cfgc, TaskName{"pi0eta-to-gammagamma-emcemc"}), }; } diff --git a/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMCEMCEMC.cxx b/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMCEMCEMC.cxx index 982eba314e1..36fa6c724b9 100644 --- a/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMCEMCEMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMCEMCEMC.cxx @@ -26,11 +26,11 @@ using namespace o2::aod; using namespace o2::framework; using namespace o2::aod::pwgem::photonmeson::photonpair; -using MyEMCClusters = soa::Join; +using MyEMCClusters = soa::Join; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask>(cfgc, TaskName{"pi0eta-to-gammagamma-mc-emcemc"}), + adaptAnalysisTask>(cfgc, TaskName{"pi0eta-to-gammagamma-mc-emcemc"}), }; } diff --git a/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaPCMDalitzEE.cxx b/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaPCMDalitzEE.cxx index b2c39f2093a..27d1c1cca7d 100644 --- a/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaPCMDalitzEE.cxx +++ b/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaPCMDalitzEE.cxx @@ -18,16 +18,18 @@ #include "PWGEM/PhotonMeson/Utils/PairUtilities.h" #include +#include #include #include using namespace o2; +using namespace o2::soa; using namespace o2::aod; using namespace o2::framework; using namespace o2::aod::pwgem::photonmeson::photonpair; -using MyV0Photons = o2::soa::Filtered>; -using MyPrimaryElectrons = o2::soa::Filtered>; +using MyV0Photons = Filtered>; +using MyPrimaryElectrons = Filtered>; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { diff --git a/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaPCMEMC.cxx b/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaPCMEMC.cxx new file mode 100644 index 00000000000..678d0f21f47 --- /dev/null +++ b/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaPCMEMC.cxx @@ -0,0 +1,38 @@ +// 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 Pi0EtaToGammaGammaEMCEMC.cxx +/// \brief This code loops over photons and makes pairs for neutral mesons analyses for EMC-EMC. +/// \author M. Hemmer, marvin.hemmer@cern.ch + +#include "PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h" +#include "PWGEM/PhotonMeson/DataModel/gammaTables.h" +#include "PWGEM/PhotonMeson/Utils/PairUtilities.h" + +#include +#include +#include + +using namespace o2; +using namespace o2::soa; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::aod::pwgem::photonmeson::photonpair; + +using MyV0Photons = Filtered>; +using MyEMCClusters = soa::Join; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask>(cfgc, TaskName{"pi0eta-to-gammagamma-pcmemc"}), + }; +} diff --git a/PWGEM/PhotonMeson/Tasks/SinglePhoton.cxx b/PWGEM/PhotonMeson/Tasks/SinglePhoton.cxx index cbb87e23b37..5d5972e3685 100644 --- a/PWGEM/PhotonMeson/Tasks/SinglePhoton.cxx +++ b/PWGEM/PhotonMeson/Tasks/SinglePhoton.cxx @@ -250,7 +250,7 @@ struct SinglePhoton { { bool is_selected = false; if constexpr (photontype == EMDetType::kPCM) { - is_selected = cut1.template IsSelected(g1); + is_selected = cut1.template IsSelected(g1); } else if constexpr (photontype == EMDetType::kPHOS) { is_selected = cut1.template IsSelected(g1); // dummy, because track matching is not ready. //} else if constexpr (photontype == EMDetType::kEMC) { diff --git a/PWGEM/PhotonMeson/Tasks/SinglePhotonMC.cxx b/PWGEM/PhotonMeson/Tasks/SinglePhotonMC.cxx index 7a5b3e42120..b9aeedf9ffa 100644 --- a/PWGEM/PhotonMeson/Tasks/SinglePhotonMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/SinglePhotonMC.cxx @@ -8,11 +8,10 @@ // 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. -// -// ======================== -// -// This code loops over photon candidate and fill histograms -// Please write to: daiki.sekihata@cern.ch + +/// \file SinglePhotonMC.cxx +/// \brief HeadThis code loops over photon candidate and fill histograms +/// \author D. Sekihata, daiki.sekihata@cern.ch #include "EMPhotonEventCut.h" @@ -252,7 +251,7 @@ struct SinglePhotonMC { { bool is_selected = false; if constexpr (photontype == EMDetType::kPCM) { - is_selected = cut1.template IsSelected(g1); + is_selected = cut1.template IsSelected(g1); } else if constexpr (photontype == EMDetType::kPHOS) { is_selected = cut1.template IsSelected(g1); // dummy, because track matching is not ready. //} else if constexpr (photontype == EMDetType::kEMC) { diff --git a/PWGEM/PhotonMeson/Tasks/TagAndProbe.cxx b/PWGEM/PhotonMeson/Tasks/TagAndProbe.cxx index 2def85ae4d2..a2d09adcc0f 100644 --- a/PWGEM/PhotonMeson/Tasks/TagAndProbe.cxx +++ b/PWGEM/PhotonMeson/Tasks/TagAndProbe.cxx @@ -29,6 +29,7 @@ #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" +#include #include #include #include @@ -285,15 +286,15 @@ struct TagAndProbe { for (auto& g1 : photons1_coll) { if constexpr (pairtype == PairType::kPCMPCM) { - if (!tagcut.template IsSelected(g1)) { + if (!tagcut.template IsSelected(g1)) { continue; } } else if constexpr (pairtype == PairType::kPHOSPHOS) { - if (!tagcut.template IsSelected(g1)) { + if (!tagcut.template IsSelected(g1)) { continue; } } else if constexpr (pairtype == PairType::kEMCEMC) { - if (!tagcut.template IsSelected(g1)) { + if (!tagcut.template IsSelected(g1)) { continue; } } @@ -318,15 +319,15 @@ struct TagAndProbe { reinterpret_cast(list_pair_ss->FindObject(Form("%s_%s", tagcut.GetName(), probecut.GetName()))->FindObject(paircut.GetName())->FindObject("hMggPt_Probe_Same"))->Fill(v12.M(), v2.Pt()); if constexpr (pairtype == PairType::kPCMPCM) { - if (!probecut.template IsSelected(g2)) { + if (!probecut.template IsSelected(g2)) { continue; } } else if constexpr (pairtype == PairType::kPHOSPHOS) { - if (!probecut.template IsSelected(g2)) { + if (!probecut.template IsSelected(g2)) { continue; } } else if constexpr (pairtype == PairType::kEMCEMC) { - if (!probecut.template IsSelected(g2)) { + if (!probecut.template IsSelected(g2)) { continue; } } @@ -379,15 +380,15 @@ struct TagAndProbe { for (auto& g1 : photons_coll1) { if constexpr (pairtype == PairType::kPCMPCM) { - if (!tagcut.template IsSelected(g1)) { + if (!tagcut.template IsSelected(g1)) { continue; } } else if constexpr (pairtype == PairType::kPHOSPHOS) { - if (!tagcut.template IsSelected(g1)) { + if (!tagcut.template IsSelected(g1)) { continue; } } else if constexpr (pairtype == PairType::kEMCEMC) { - if (!tagcut.template IsSelected(g1)) { + if (!tagcut.template IsSelected(g1)) { continue; } } @@ -409,15 +410,15 @@ struct TagAndProbe { reinterpret_cast(list_pair_ss->FindObject(Form("%s_%s", tagcut.GetName(), probecut.GetName()))->FindObject(paircut.GetName())->FindObject("hMggPt_Probe_Mixed"))->Fill(v12.M(), v2.Pt()); if constexpr (pairtype == PairType::kPCMPCM) { - if (!probecut.template IsSelected(g2)) { + if (!probecut.template IsSelected(g2)) { continue; } } else if constexpr (pairtype == PairType::kPHOSPHOS) { - if (!probecut.template IsSelected(g2)) { + if (!probecut.template IsSelected(g2)) { continue; } } else if constexpr (pairtype == PairType::kEMCEMC) { - if (!probecut.template IsSelected(g2)) { + if (!probecut.template IsSelected(g2)) { continue; } } @@ -439,7 +440,7 @@ struct TagAndProbe { if (photons_coll.size() < 3) { return; } - const double rotationAngle = M_PI / 2.0; // rotaion angle 90° + const double rotationAngle = o2::constants::math::PIHalf; // rotaion angle 90° // ROOT::Math::XYZVector meson3D(meson.Px(), meson.Py(), meson.Pz()); ROOT::Math::AxisAngle rotationAxis(meson.Vect(), rotationAngle); // LOG(info) << "rotationAxis.Angle() = " << rotationAxis.Angle(); @@ -461,7 +462,7 @@ struct TagAndProbe { // only combine rotated photons with other photons continue; } - if (!cut.template IsSelected(photon)) { + if (!cut.template IsSelected(photon)) { continue; } diff --git a/PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx b/PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx index 607a562601d..d9d75bf8ade 100644 --- a/PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx +++ b/PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx @@ -666,7 +666,7 @@ struct CalibTaskEmc { if (photon.globalIndex() == ig2) { continue; } - if (!(fPCMPhotonCut.IsSelected(photon))) { + if (!(fPCMPhotonCut.IsSelected(photon))) { continue; } ROOT::Math::PtEtaPhiMVector photon3(photon.pt(), photon.eta(), photon.phi(), 0.); @@ -839,7 +839,7 @@ struct CalibTaskEmc { } } for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(photonsEMCPerCollision, photonsPCMPerCollision))) { - if (!(fEMCCut.IsSelected(g1)) || !(fPCMPhotonCut.IsSelected(g2))) { + if (!(fEMCCut.IsSelected(g1)) || !(fPCMPhotonCut.IsSelected(g2))) { continue; } @@ -1005,7 +1005,7 @@ struct CalibTaskEmc { runBefore = runNow; } for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(photonEMC, photonPCM))) { - if (!(fEMCCut.IsSelected(g1)) || !(fPCMPhotonCut.IsSelected(g2))) { + if (!(fEMCCut.IsSelected(g1)) || !(fPCMPhotonCut.IsSelected(g2))) { continue; } // Cut edge clusters away, similar to rotation method to ensure same acceptance is used diff --git a/PWGEM/PhotonMeson/Tasks/emcalQC.cxx b/PWGEM/PhotonMeson/Tasks/emcalQC.cxx index 6c3ef79b152..b7f632512c2 100644 --- a/PWGEM/PhotonMeson/Tasks/emcalQC.cxx +++ b/PWGEM/PhotonMeson/Tasks/emcalQC.cxx @@ -199,8 +199,8 @@ struct EmcalQC { fRegistry.fill(HIST("Cluster/hClusterQualityCuts"), 0., cluster.e(), collision.weight()); // Define two boleans to see, whether the cluster "survives" the EMC cluster cuts to later check, whether the cuts in this task align with the ones in EMCPhotonCut.h: - bool survivesIsSelectedEMCalCuts = true; // Survives "manual" cuts listed in this task - bool survivesIsSelectedCuts = fEMCCut.IsSelected(cluster); // Survives the cutlist defines in EMCPhotonCut.h, which is also used in the Pi0Eta task + bool survivesIsSelectedEMCalCuts = true; // Survives "manual" cuts listed in this task + bool survivesIsSelectedCuts = fEMCCut.IsSelected(cluster); // Survives the cutlist defines in EMCPhotonCut.h, which is also used in the Pi0Eta task for (int icut = 0; icut < static_cast(EMCPhotonCut::EMCPhotonCuts::kNCuts); icut++) { // Loop through different cut observables EMCPhotonCut::EMCPhotonCuts specificcut = static_cast(icut); diff --git a/PWGEM/PhotonMeson/Tasks/gammaConversions.cxx b/PWGEM/PhotonMeson/Tasks/gammaConversions.cxx index 37e4bf4dde1..b9268ab3907 100644 --- a/PWGEM/PhotonMeson/Tasks/gammaConversions.cxx +++ b/PWGEM/PhotonMeson/Tasks/gammaConversions.cxx @@ -9,27 +9,38 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +/// \file gammaConversions.cxx /// \brief perform photon conversion analysis on V0 candidates from aod::StoredV0Datas -/// dependencies: o2-analysis-lf-lambdakzerobuilder /// \author stephan.friedrich.stiefelmaier@cern.ch +/// dependencies: o2-analysis-lf-lambdakzerobuilder #include "PWGEM/PhotonMeson/Tasks/gammaConversions.h" -#include -#include -#include -#include - #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" #include "PWGEM/PhotonMeson/Utils/gammaConvDefinitions.h" -#include "Framework/AnalysisTask.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/runDataProcessing.h" #include "Common/Core/RecoDecay.h" +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include #include -#include // for ATan2, Cos, Sin, Sqrt + +#include +#include +#include +#include +#include +#include using namespace o2; using namespace o2::framework; @@ -70,7 +81,7 @@ struct GammaConversions { Configurable fV0QtPtMultiplicator{"fV0QtPtMultiplicator", 0.125, "Multiply pt of V0s by this value to get the 2nd denominator in the armenteros cut. The products maximum value is fV0QtMax."}; Configurable fV0QtMax{"fV0QtMax", 0.050, "the maximum value of the product, that is the maximum qt"}; Configurable fLineCutZ0{"fLineCutZ0", 12.0, "The offset for the linecute used in the Z vs R plot"}; - Configurable fLineCutZRSlope{"fLineCutZRSlope", static_cast(TMath::Tan(2 * TMath::ATan(TMath::Exp(-fTruePhotonEtaMax)))), "The slope for the line cut"}; + Configurable fLineCutZRSlope{"fLineCutZRSlope", std::tan(2.f * std::atan(std::exp(-fTruePhotonEtaMax.value))), "The slope for the line cut"}; Configurable fPhysicalPrimaryOnly{"fPhysicalPrimaryOnly", true, "fPhysicalPrimaryOnly"}; std::map fPhotonCutLabels{ @@ -188,7 +199,7 @@ struct GammaConversions { auto lHisto = theContainer.find(theHistoName); if (lHisto != theContainer.end()) { TAxis* lXaxis = std::get>(lHisto->second)->GetXaxis(); - for (auto& lPairIt : theLables) { + for (const auto& lPairIt : theLables) { lXaxis->SetBinLabel(static_cast(lPairIt.first) + 1, lPairIt.second.data()); } } @@ -198,7 +209,7 @@ struct GammaConversions { auto lHisto = theContainer.find(theHistoName); if (lHisto != theContainer.end()) { TAxis* lXaxis = std::get>(lHisto->second)->GetXaxis(); - for (auto& lPairIt : theLables) { + for (const auto& lPairIt : theLables) { lXaxis->SetBinLabel(static_cast(lPairIt.first) + 1, lPairIt.second.data()); } } @@ -316,7 +327,7 @@ struct GammaConversions { fillReconstructedInfoHistogramsI(kAfterRecCuts); - return kTRUE; + return true; } template @@ -384,7 +395,7 @@ struct GammaConversions { return false; } - if (theMcPhoton.pdgCode() != 22) { + if (theMcPhoton.pdgCode() != PDG_t::kGamma) { fillV0McValidationHisto(eV0McValidation::kNoPhoton); // add here a fillRejectedV0HistosI() call if interested in properties of non-gamma V0s return false; @@ -482,27 +493,27 @@ struct GammaConversions { if (PDGCode[0] == 0) { fillTH2(theContainer, "hDecays", static_cast(eV0Decays::nomcparticle), 0); - } else if (sameMother && ((PDGCode[0] == 11 && PDGCode[1] == -11) || (PDGCode[0] == -11 && PDGCode[1] == 11))) { + } else if (sameMother && ((PDGCode[0] == PDG_t::kElectron && PDGCode[1] == PDG_t::kPositron) || (PDGCode[0] == PDG_t::kPositron && PDGCode[1] == PDG_t::kElectron))) { fillTH2(theContainer, "hDecays", static_cast(eV0Decays::ee1), MCV0p); - } else if (!sameMother && ((PDGCode[0] == 11 && PDGCode[1] == -11) || (PDGCode[0] == -11 && PDGCode[1] == 11))) { + } else if (!sameMother && ((PDGCode[0] == PDG_t::kElectron && PDGCode[1] == PDG_t::kPositron) || (PDGCode[0] == PDG_t::kPositron && PDGCode[1] == PDG_t::kElectron))) { fillTH2(theContainer, "hDecays", static_cast(eV0Decays::ee2), MCV0p); - } else if ((PDGCode[0] == 11 && PDGCode[1] == 211) || (PDGCode[0] == -11 && PDGCode[1] == -211)) { + } else if ((PDGCode[0] == PDG_t::kElectron && PDGCode[1] == PDG_t::kPiPlus) || (PDGCode[0] == PDG_t::kPositron && PDGCode[1] == PDG_t::kPiMinus)) { fillTH2(theContainer, "hDecays", static_cast(eV0Decays::epi), MCV0p); - } else if ((PDGCode[0] == 11 && PDGCode[1] == 321) || (PDGCode[0] == -11 && PDGCode[1] == -321)) { + } else if ((PDGCode[0] == PDG_t::kElectron && PDGCode[1] == PDG_t::kKPlus) || (PDGCode[0] == PDG_t::kPositron && PDGCode[1] == PDG_t::kKMinus)) { fillTH2(theContainer, "hDecays", static_cast(eV0Decays::ek), MCV0p); - } else if ((PDGCode[0] == 11 && PDGCode[1] == 2212) || (PDGCode[0] == -11 && PDGCode[1] == -2212)) { + } else if ((PDGCode[0] == PDG_t::kElectron && PDGCode[1] == PDG_t::kProton) || (PDGCode[0] == PDG_t::kPositron && PDGCode[1] == PDG_t::kProtonBar)) { fillTH2(theContainer, "hDecays", static_cast(eV0Decays::ep), MCV0p); - } else if ((PDGCode[0] == 11 && PDGCode[1] == -13) || (PDGCode[0] == -11 && PDGCode[1] == 13)) { + } else if ((PDGCode[0] == PDG_t::kElectron && PDGCode[1] == PDG_t::kMuonPlus) || (PDGCode[0] == PDG_t::kPositron && PDGCode[1] == PDG_t::kMuonMinus)) { fillTH2(theContainer, "hDecays", static_cast(eV0Decays::emu), MCV0p); - } else if ((PDGCode[0] == 211 && PDGCode[1] == -211) || (PDGCode[0] == -211 && PDGCode[1] == 211)) { + } else if ((PDGCode[0] == PDG_t::kPiPlus && PDGCode[1] == PDG_t::kPiMinus) || (PDGCode[0] == PDG_t::kPiMinus && PDGCode[1] == PDG_t::kPiPlus)) { fillTH2(theContainer, "hDecays", static_cast(eV0Decays::pipi), MCV0p); - } else if ((PDGCode[0] == 211 && PDGCode[1] == -321) || (PDGCode[0] == -211 && PDGCode[1] == 321)) { + } else if ((PDGCode[0] == PDG_t::kPiPlus && PDGCode[1] == PDG_t::kKMinus) || (PDGCode[0] == PDG_t::kPiMinus && PDGCode[1] == PDG_t::kKPlus)) { fillTH2(theContainer, "hDecays", static_cast(eV0Decays::pik), MCV0p); - } else if ((PDGCode[0] == 211 && PDGCode[1] == -2212) || (PDGCode[0] == -211 && PDGCode[1] == 2212)) { + } else if ((PDGCode[0] == PDG_t::kPiPlus && PDGCode[1] == PDG_t::kProtonBar) || (PDGCode[0] == PDG_t::kPiMinus && PDGCode[1] == PDG_t::kProton)) { fillTH2(theContainer, "hDecays", static_cast(eV0Decays::pip), MCV0p); - } else if ((PDGCode[0] == 211 && PDGCode[1] == 13) || (PDGCode[0] == -211 && PDGCode[1] == -13)) { + } else if ((PDGCode[0] == PDG_t::kPiPlus && PDGCode[1] == PDG_t::kMuonMinus) || (PDGCode[0] == PDG_t::kPiMinus && PDGCode[1] == PDG_t::kMuonPlus)) { fillTH2(theContainer, "hDecays", static_cast(eV0Decays::pimu), MCV0p); - } else if ((PDGCode[0] == 2212 && PDGCode[1] == -321) || (PDGCode[0] == -2212 && PDGCode[1] == 321) || (PDGCode[0] == 2212 && PDGCode[1] == 13) || (PDGCode[0] == -2212 && PDGCode[1] == -13)) { + } else if ((PDGCode[0] == PDG_t::kProton && PDGCode[1] == PDG_t::kKMinus) || (PDGCode[0] == PDG_t::kProtonBar && PDGCode[1] == PDG_t::kKPlus) || (PDGCode[0] == PDG_t::kProton && PDGCode[1] == PDG_t::kMuonMinus) || (PDGCode[0] == PDG_t::kProtonBar && PDGCode[1] == PDG_t::kMuonPlus)) { fillTH2(theContainer, "hDecays", static_cast(eV0Decays::pKmu), MCV0p); } else { fillTH2(theContainer, "hDecays", static_cast(eV0Decays::other), MCV0p); @@ -618,7 +629,7 @@ struct GammaConversions { theCollision.posZ()); auto theV0s_per_coll = theV0s.sliceBy(perCollision, theCollision.globalIndex()); - for (auto& lV0 : theV0s_per_coll) { + for (const auto& lV0 : theV0s_per_coll) { float lV0CosinePA = lV0.cospa(); @@ -645,7 +656,7 @@ struct GammaConversions { theCollision.posZ()); auto theV0s_per_coll = theV0s.sliceBy(perCollision, theCollision.globalIndex()); - for (auto& lV0 : theV0s) { + for (const auto& lV0 : theV0s) { float lV0CosinePA = lV0.cospa(); auto pos = lV0.posTrack_as(); @@ -776,31 +787,31 @@ struct GammaConversions { bool trackPassesCuts(const T& theTrack) { // single track eta cut - if (TMath::Abs(theTrack.eta()) > fTrackEtaMax) { + if (std::abs(theTrack.eta()) > fTrackEtaMax) { fillV0SelectionHisto(ePhotonCuts::kTrackEta); - return kFALSE; + return false; } // single track pt cut if (theTrack.pt() < fTrackPtMin) { fillV0SelectionHisto(ePhotonCuts::kTrackPt); - return kFALSE; + return false; } if (!(selectionPIDTPC_track(theTrack))) { - return kFALSE; + return false; } if (theTrack.tpcFoundOverFindableCls() < fMinTPCFoundOverFindableCls) { fillV0SelectionHisto(ePhotonCuts::kTPCFoundOverFindableCls); - return kFALSE; + return false; } if (theTrack.tpcCrossedRowsOverFindableCls() < fMinTPCCrossedRowsOverFindableCls) { fillV0SelectionHisto(ePhotonCuts::kTPCCrossedRowsOverFindableCls); - return kFALSE; + return false; } - return kTRUE; + return true; } template @@ -820,30 +831,30 @@ struct GammaConversions { { if (theV0.v0radius() < fV0RMin || theV0.v0radius() > fV0RMax) { fillV0SelectionHisto(ePhotonCuts::kV0Radius); - return kFALSE; + return false; } if (fV0PhotonAsymmetryMax > 0. && !ArmenterosQtCut(theV0.alpha(), theV0.qtarm(), theV0.pt())) { fillV0SelectionHisto(ePhotonCuts::kArmenteros); - return kFALSE; + return false; } // if (fV0PsiPairMax > 0. && TMath::Abs(theV0.psipair()) > fV0PsiPairMax) { - if (fV0PsiPairMax > 0. && TMath::Abs(0.f) > fV0PsiPairMax) { + if (fV0PsiPairMax > 0. && std::abs(0.f) > fV0PsiPairMax) { fillV0SelectionHisto(ePhotonCuts::kPsiPair); - return kFALSE; + return false; } if (fV0CosPAngleMin > 0. && theV0CosinePA < fV0CosPAngleMin) { fillV0SelectionHisto(ePhotonCuts::kCosinePA); - return kFALSE; + return false; } - if (TMath::Abs(theV0.vz()) > fLineCutZ0 + theV0.v0radius() * fLineCutZRSlope) { // as long as z recalculation is not fixed use this + if (std::abs(theV0.vz()) > fLineCutZ0 + theV0.v0radius() * fLineCutZRSlope) { // as long as z recalculation is not fixed use this fillV0SelectionHisto(ePhotonCuts::kRZLine); - return kFALSE; + return false; } - return kTRUE; + return true; } template @@ -869,17 +880,17 @@ struct GammaConversions { theV0); } - Bool_t ArmenterosQtCut(Double_t theAlpha, Double_t theQt, Double_t thePt) + bool ArmenterosQtCut(float theAlpha, float theQt, float thePt) { // in AliPhysics this is the cut for if fDo2DQt && fDoQtGammaSelection == 2 - Float_t lQtMaxPtDep = fV0QtPtMultiplicator * thePt; + float lQtMaxPtDep = fV0QtPtMultiplicator * thePt; if (lQtMaxPtDep > fV0QtMax) { lQtMaxPtDep = fV0QtMax; } - if ((TMath::Power(theAlpha / fV0PhotonAsymmetryMax, 2) + TMath::Power(theQt / lQtMaxPtDep, 2)) >= 1) { - return kFALSE; + if ((std::pow(theAlpha / fV0PhotonAsymmetryMax, 2) + std::pow(theQt / lQtMaxPtDep, 2)) >= 1) { + return false; } - return kTRUE; + return true; } template @@ -888,7 +899,7 @@ struct GammaConversions { // TPC Electron Line if (fPIDnSigmaElectronMin && (theTrack.tpcNSigmaEl() < fPIDnSigmaElectronMin || theTrack.tpcNSigmaEl() > fPIDnSigmaElectronMax)) { fillV0SelectionHisto(ePhotonCuts::kElectronPID); - return kFALSE; + return false; } // TPC Pion Line @@ -897,17 +908,17 @@ struct GammaConversions { if (theTrack.p() < fPIDPionRejectionPBoarder) { if (theTrack.tpcNSigmaEl() > fPIDnSigmaElectronMin && theTrack.tpcNSigmaEl() < fPIDnSigmaElectronMax && theTrack.tpcNSigmaPi() < fPIDnSigmaAbovePionLineLowPMin) { fillV0SelectionHisto(ePhotonCuts::kPionRejLowMom); - return kFALSE; + return false; } // High Pt Pion rej } else { if (theTrack.tpcNSigmaEl() > fPIDnSigmaElectronMin && theTrack.tpcNSigmaEl() < fPIDnSigmaElectronMax && theTrack.tpcNSigmaPi() < fPIDnSigmaAbovePionLineHighPMin) { fillV0SelectionHisto(ePhotonCuts::kPionRejHighMom); - return kFALSE; + return false; } } } - return kTRUE; + return true; } }; diff --git a/PWGEM/PhotonMeson/Tasks/gammaConversions.h b/PWGEM/PhotonMeson/Tasks/gammaConversions.h index 9e1c9ebe7db..7eb85f26efe 100644 --- a/PWGEM/PhotonMeson/Tasks/gammaConversions.h +++ b/PWGEM/PhotonMeson/Tasks/gammaConversions.h @@ -9,22 +9,35 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +/// \file gammaConversions.h /// \brief contains type definitions for gammaConversions.cxx /// \author stephan.friedrich.stiefelmaier@cern.ch -#include "Framework/AnalysisTask.h" +#ifndef PWGEM_PHOTONMESON_TASKS_GAMMACONVERSIONS_H_ +#define PWGEM_PHOTONMESON_TASKS_GAMMACONVERSIONS_H_ -using namespace o2::framework; +#include +#include +#include -typedef std::map mapStringHistPtr; +#include +#include + +#include +#include +#include +#include +#include + +typedef std::map mapStringHistPtr; // define this in order to have a constructor of the HistogramSpec which copies the name into the title and to have dataOnly flag struct MyHistogramSpec { - MyHistogramSpec(char const* const name_, char const* const title_, HistogramConfigSpec config_, bool callSumw2_ = false, bool dataOnly_ = false) + MyHistogramSpec(char const* const name_, char const* const title_, o2::framework::HistogramConfigSpec config_, bool callSumw2_ = false, bool dataOnly_ = false) : m{name_, title_, config_, callSumw2_}, fDataOnly(dataOnly_) {} - MyHistogramSpec(char const* const name_, HistogramConfigSpec config_, bool callSumw2_ = false, bool dataOnly_ = false) + MyHistogramSpec(char const* const name_, o2::framework::HistogramConfigSpec config_, bool callSumw2_ = false, bool dataOnly_ = false) : m{name_, name_, config_, callSumw2_}, fDataOnly(dataOnly_) {} - HistogramSpec m{}; + o2::framework::HistogramSpec m{}; bool fDataOnly{false}; }; @@ -82,11 +95,11 @@ enum eV0HistoFlavor { kRec, kRes }; struct tV0Kind { - tV0Kind(std::string const& thePath) : mPath{thePath} {} + explicit tV0Kind(std::string const& thePath) : mPath{thePath} {} // todo: remove template - void appendSuffixToTitleI(HistPtr& theHistPtr, std::string const* theSuffix) + void appendSuffixToTitleI(o2::framework::HistPtr& theHistPtr, std::string const* theSuffix) { auto lHisto = std::get>(theHistPtr); if (lHisto) { @@ -98,7 +111,7 @@ struct tV0Kind { } // todo: bind tV0Kind to gammaConverions such that theOfficialRegistry and theCheckDataOnly are already known - void addHistosToOfficalRegistry(HistogramRegistry& theOfficialRegistry, + void addHistosToOfficalRegistry(o2::framework::HistogramRegistry& theOfficialRegistry, std::vector const& theHistoDefinitions, std::string const* theSuffix = nullptr, bool theCheckDataOnly = false) @@ -109,14 +122,14 @@ struct tV0Kind { } std::string lFullName(mPath + tHisto.m.name + (theSuffix ? *theSuffix : std::string(""))); LOGF(info, "adding %s %d", lFullName, tHisto.fDataOnly); - HistPtr lHistPtr = theOfficialRegistry.add(lFullName.data(), tHisto.m.title.data(), tHisto.m.config); + o2::framework::HistPtr lHistPtr = theOfficialRegistry.add(lFullName.data(), tHisto.m.title.data(), tHisto.m.config); mContainer.insert(std::pair{tHisto.m.name, lHistPtr}); // todo ugly: remove if (theSuffix) { - if (tHisto.m.config.type == kTH1F) { + if (tHisto.m.config.type == o2::framework::kTH1F) { appendSuffixToTitleI(lHistPtr, theSuffix); - } else if (tHisto.m.config.type == kTH2F) { + } else if (tHisto.m.config.type == o2::framework::kTH2F) { appendSuffixToTitleI(lHistPtr, theSuffix); } } @@ -128,38 +141,40 @@ struct tV0Kind { }; struct tMotherDirV0Kinds { - tMotherDirV0Kinds(std::string const& thePath) : mPath{thePath} {} + explicit tMotherDirV0Kinds(std::string const& thePath) : mPath{thePath} {} std::string mPath{""}; - tV0Kind mV0Kind[4]{mPath + "Rec/", mPath + "MCTrue/", mPath + "MCVal/", mPath + "Res/"}; + tV0Kind mV0Kind[4]{tV0Kind{mPath + "Rec/"}, tV0Kind{mPath + "MCTrue/"}, tV0Kind{mPath + "MCVal/"}, tV0Kind{mPath + "Res/"}}; }; struct tMotherDirRejMc { - tMotherDirRejMc(std::string const& thePath) : mPath{thePath} {} + explicit tMotherDirRejMc(std::string const& thePath) : mPath{thePath} {} std::string mPath{""}; - tMotherDirV0Kinds mBeforeAfterRecCuts[2]{mPath + "beforeRecCuts/", mPath + "afterRecCuts/"}; + tMotherDirV0Kinds mBeforeAfterRecCuts[2]{tMotherDirV0Kinds{mPath + "beforeRecCuts/"}, tMotherDirV0Kinds{mPath + "afterRecCuts/"}}; }; // for collision track v0 struct tHistoFolderCTV { - tHistoFolderCTV(std::string const& thePath) : mPath{thePath} {} + explicit tHistoFolderCTV(std::string const& thePath) : mPath{thePath} {} std::string mPath{""}; tV0Kind mSpecialHistos{mPath}; - tMotherDirV0Kinds mBeforeAfterRecCuts[2]{mPath + "beforeRecCuts/", - mPath + "afterRecCuts/"}; + tMotherDirV0Kinds mBeforeAfterRecCuts[2]{tMotherDirV0Kinds{mPath + "beforeRecCuts/"}, + tMotherDirV0Kinds{mPath + "afterRecCuts/"}}; - tMotherDirRejMc mRejectedByMc[2]{mPath + "rejectedByMc/kNoPhysicalPrimary/", - mPath + "rejectedByMc/kOutsideMCEtaAcc/"}; + tMotherDirRejMc mRejectedByMc[2]{tMotherDirRejMc{mPath + "rejectedByMc/kNoPhysicalPrimary/"}, + tMotherDirRejMc{mPath + "rejectedByMc/kOutsideMCEtaAcc/"}}; }; enum class eMcRejectedSaved { kNoPhysicalPrimary, kOutsideMCEtaAcc }; struct tHistoRegistry { - tHistoRegistry(std::string const& thePath) : mPath{thePath} {} + explicit tHistoRegistry(std::string const& thePath) : mPath{thePath} {} std::string mPath{""}; tHistoFolderCTV mCollision{mPath + "Collision/"}; tHistoFolderCTV mTrack{mPath + "Track/"}; tHistoFolderCTV mV0{mPath + "V0/"}; -}; \ No newline at end of file +}; + +#endif // PWGEM_PHOTONMESON_TASKS_GAMMACONVERSIONS_H_ diff --git a/PWGEM/PhotonMeson/Tasks/gammaConversionsTruthOnlyMc.cxx b/PWGEM/PhotonMeson/Tasks/gammaConversionsTruthOnlyMc.cxx index 259a53cd037..ae2ad0c548f 100644 --- a/PWGEM/PhotonMeson/Tasks/gammaConversionsTruthOnlyMc.cxx +++ b/PWGEM/PhotonMeson/Tasks/gammaConversionsTruthOnlyMc.cxx @@ -9,19 +9,22 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -// SFS todo: update description +/// \file gammaConversionsTruthOnlyMc.cxx /// \brief extract relevant mc truth information that allows to compute efficiency, purity and more quantities for the photon conversion analysis. -/// dependencies: none /// \author stephan.friedrich.stiefelmaier@cern.ch +/// dependencies: none +// SFS todo: update description #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" #include "PWGEM/PhotonMeson/Utils/gammaConvDefinitions.h" -#include "TVector3.h" +#include +#include +#include + +#include -#include "Framework/AnalysisTask.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/runDataProcessing.h" +#include using namespace o2; using namespace o2::framework; @@ -34,7 +37,7 @@ struct gammaConversionsTruthOnlyMc { Configurable fV0RMin{"fV0RMin", 0., "minimum conversion radius of the V0s"}; Configurable fV0RMax{"fV0RMax", 180., "maximum conversion radius of the V0s"}; Configurable LineCutZ0{"fLineCutZ0", 7.0, "The offset for the linecute used in the Z vs R plot"}; - Configurable LineCutZRSlope{"LineCutZRSlope", (float)TMath::Tan(2 * TMath::ATan(TMath::Exp(-fEtaMax))), "The slope for the line cut"}; + Configurable LineCutZRSlope{"LineCutZRSlope", std::tan(2.f * std::atan(std::exp(-fEtaMax.value))), "The slope for the line cut"}; HistogramRegistry registry{ "registry", @@ -111,7 +114,7 @@ struct gammaConversionsTruthOnlyMc { return false; } - if (TMath::Abs(theMcGamma.conversionZ()) > LineCutZ0 + theMcGamma.v0Radius() * LineCutZRSlope) { + if (std::abs(theMcGamma.conversionZ()) > LineCutZ0 + theMcGamma.v0Radius() * LineCutZRSlope) { return false; } return true; @@ -130,11 +133,11 @@ struct gammaConversionsTruthOnlyMc { } registry.fill(HIST("hCollisionZ_MCTrue"), theMcCollision.posZ()); - for (auto& lCollision : theCollisions) { + for (const auto& lCollision : theCollisions) { registry.fill(HIST("hCollisionZ_MCRec"), lCollision.posZ()); } - for (auto& lMcGamma : theMcGammas) { + for (const auto& lMcGamma : theMcGammas) { if (!photonPassesCuts(lMcGamma)) { continue; @@ -163,4 +166,4 @@ struct gammaConversionsTruthOnlyMc { WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc)}; -} \ No newline at end of file +} diff --git a/PWGEM/PhotonMeson/Tasks/pcmQC.cxx b/PWGEM/PhotonMeson/Tasks/pcmQC.cxx index a2a06a832e9..3f8130abe0d 100644 --- a/PWGEM/PhotonMeson/Tasks/pcmQC.cxx +++ b/PWGEM/PhotonMeson/Tasks/pcmQC.cxx @@ -8,23 +8,35 @@ // 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. -// -// ======================== -// -// This code runs loop over v0 photons for PCM QC. -// Please write to: daiki.sekihata@cern.ch + +/// \file MaterialBudget.cxx +/// \brief This code runs loop over v0 photons for PCM QC. +/// \author Daiki Sekihata, daiki.sekihata@cern.ch #include "PWGEM/PhotonMeson/Core/EMPhotonEventCut.h" #include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" -#include "PWGEM/PhotonMeson/Utils/PCMUtilities.h" -#include "Framework/ASoAHelpers.h" +#include "Common/CCDB/EventSelectionParams.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" + #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" #include "Framework/runDataProcessing.h" +#include +#include +#include +#include +#include +#include +#include + +#include +#include #include +#include #include using namespace o2; @@ -134,7 +146,7 @@ struct PCMQC { // v0 info fRegistry.add("V0/hPt", "pT;p_{T,#gamma} (GeV/c)", kTH1F, {{2000, 0.0f, 20}}, false); - fRegistry.add("V0/hEtaPhi", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{90, 0, 2 * M_PI}, {200, -1.0f, 1.0f}}, false); + fRegistry.add("V0/hEtaPhi", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{90, 0, o2::constants::math::TwoPI}, {200, -1.0f, 1.0f}}, false); fRegistry.add("V0/hXY", "conversion point in XY;V_{x} (cm);V_{y} (cm)", kTH2F, {{400, -100.0f, 100.0f}, {400, -100.0f, 100.0f}}, false); fRegistry.add("V0/hRZ", "conversion point in RZ;Z (cm);R_{xy} (cm)", kTH2F, {{200, -100, 100}, {200, 0.0f, 100.0f}}, false); fRegistry.add("V0/hCosPA", "V0CosPA;cosine pointing angle in 3D", kTH1F, {{100, 0.99f, 1.0f}}, false); @@ -150,13 +162,13 @@ struct PCMQC { fRegistry.add("V0/hKFChi2vsX", "KF chi2 vs. conversion point in X;X (cm);KF chi2/NDF", kTH2F, {{200, -100.0f, 100.0f}, {100, 0.f, 100.0f}}, false); fRegistry.add("V0/hKFChi2vsY", "KF chi2 vs. conversion point in Y;Y (cm);KF chi2/NDF", kTH2F, {{200, -100.0f, 100.0f}, {100, 0.f, 100.0f}}, false); fRegistry.add("V0/hKFChi2vsZ", "KF chi2 vs. conversion point in Z;Z (cm);KF chi2/NDF", kTH2F, {{200, -100.0f, 100.0f}, {100, 0.f, 100.0f}}, false); - fRegistry.add("V0/hsConvPoint", "photon conversion point;r_{xy} (cm);#varphi (rad.);#eta;", kTHnSparseF, {{100, 0.0f, 100}, {90, 0, 2 * M_PI}, {80, -2, +2}}, false); + fRegistry.add("V0/hsConvPoint", "photon conversion point;r_{xy} (cm);#varphi (rad.);#eta;", kTHnSparseF, {{100, 0.0f, 100}, {90, 0, o2::constants::math::TwoPI}, {80, -2, +2}}, false); fRegistry.add("V0/hNgamma", "Number of #gamma candidates per collision", kTH1F, {{101, -0.5f, 100.5f}}); // v0leg info fRegistry.add("V0Leg/hPt", "pT;p_{T,e} (GeV/c)", kTH1F, {{1000, 0.0f, 10}}, false); fRegistry.add("V0Leg/hQoverPt", "q/pT;q/p_{T} (GeV/c)^{-1}", kTH1F, {{1000, -50, 50}}, false); - fRegistry.add("V0Leg/hEtaPhi", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{90, 0, 2 * M_PI}, {200, -1.0f, 1.0f}}, false); + fRegistry.add("V0Leg/hEtaPhi", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{90, 0, o2::constants::math::TwoPI}, {200, -1.0f, 1.0f}}, false); fRegistry.add("V0Leg/hDCAxyz", "DCA xy vs. z;DCA_{xy} (cm);DCA_{z} (cm)", kTH2F, {{200, -50.0f, 50.0f}, {200, -50.0f, 50.0f}}, false); fRegistry.add("V0Leg/hNclsTPC", "number of TPC clusters", kTH1F, {{161, -0.5, 160.5}}, false); fRegistry.add("V0Leg/hNcrTPC", "number of TPC crossed rows", kTH1F, {{161, -0.5, 160.5}}, false); @@ -347,7 +359,7 @@ struct PCMQC { auto pos = v0.posTrack_as(); auto ele = v0.negTrack_as(); - if (!fV0PhotonCut.IsSelected(v0)) { + if (!fV0PhotonCut.IsSelected(v0)) { continue; } fillV0Info(v0); diff --git a/PWGEM/PhotonMeson/Tasks/pcmQCMC.cxx b/PWGEM/PhotonMeson/Tasks/pcmQCMC.cxx index 359f7811795..f858ec60743 100644 --- a/PWGEM/PhotonMeson/Tasks/pcmQCMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/pcmQCMC.cxx @@ -8,25 +8,39 @@ // 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. -// -// ======================== -// -// This code runs loop over v0 photons for PCM QC. -// Please write to: daiki.sekihata@cern.ch + +/// \file MaterialBudget.cxx +/// \brief This code runs loop over v0 photons for PCM QC in MC. +/// \author Daiki Sekihata, daiki.sekihata@cern.ch #include "PWGEM/Dilepton/Utils/MCUtilities.h" #include "PWGEM/PhotonMeson/Core/EMPhotonEventCut.h" #include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" #include "PWGEM/PhotonMeson/Utils/MCUtilities.h" -#include "PWGEM/PhotonMeson/Utils/PCMUtilities.h" -#include "Framework/ASoAHelpers.h" +#include "Common/CCDB/EventSelectionParams.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" + #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" #include "Framework/runDataProcessing.h" - +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include #include +#include #include using namespace o2; @@ -145,11 +159,11 @@ struct PCMQCMC { fRegistry.add("Generated/hPtY", "Generated info", kTH2F, {axis_pt, axis_rapidity}, true); fRegistry.add("Generated/hPt_ConversionPhoton", "converted photon pT;p_{T} (GeV/c)", kTH1F, {axis_pt}, true); fRegistry.add("Generated/hY_ConversionPhoton", "converted photon y;rapidity y", kTH1F, {{40, -2.0f, 2.0f}}, true); - fRegistry.add("Generated/hPhi_ConversionPhoton", "converted photon #varphi;#varphi (rad.)", kTH1F, {{180, 0, 2 * M_PI}}, true); + fRegistry.add("Generated/hPhi_ConversionPhoton", "converted photon #varphi;#varphi (rad.)", kTH1F, {{180, 0, o2::constants::math::TwoPI}}, true); fRegistry.add("Generated/hXY", "conversion point in XY MC;V_{x} (cm);V_{y} (cm)", kTH2F, {{800, -100.0f, 100.0f}, {800, -100.0f, 100.0f}}, true); fRegistry.add("Generated/hRZ", "conversion point in RZ MC;V_{z} (cm);R_{xy} (cm)", kTH2F, {{400, -100.0f, 100.0f}, {400, 0.f, 100.0f}}, true); - fRegistry.add("Generated/hRPhi", "conversion point of #varphi vs. R_{xy} MC;#varphi (rad.);R_{xy} (cm);N_{e}", kTH2F, {{360, 0.0f, 2 * M_PI}, {400, 0, 100}}, true); - fRegistry.add("Generated/hsConvPoint", "photon conversion point;r_{xy} (cm);#varphi (rad.);#eta;", kTHnSparseF, {{100, 0.0f, 100}, {90, 0, 2 * M_PI}, {80, -2, +2}}, true); + fRegistry.add("Generated/hRPhi", "conversion point of #varphi vs. R_{xy} MC;#varphi (rad.);R_{xy} (cm);N_{e}", kTH2F, {{360, 0.0f, o2::constants::math::TwoPI}, {400, 0, 100}}, true); + fRegistry.add("Generated/hsConvPoint", "photon conversion point;r_{xy} (cm);#varphi (rad.);#eta;", kTHnSparseF, {{100, 0.0f, 100}, {90, 0, o2::constants::math::TwoPI}, {80, -2, +2}}, true); } // event info @@ -178,7 +192,7 @@ struct PCMQCMC { // v0 info fRegistry.add("V0/primary/hPt", "pT;p_{T,#gamma} (GeV/c)", kTH1F, {{2000, 0.0f, 20}}, false); - fRegistry.add("V0/primary/hEtaPhi", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{90, 0, 2 * M_PI}, {200, -1.0f, 1.0f}}, false); + fRegistry.add("V0/primary/hEtaPhi", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{90, 0, o2::constants::math::TwoPI}, {200, -1.0f, 1.0f}}, false); fRegistry.add("V0/primary/hXY", "conversion point in XY;V_{x} (cm);V_{y} (cm)", kTH2F, {{400, -100.0f, 100.0f}, {400, -100.0f, 100.0f}}, false); fRegistry.add("V0/primary/hRZ", "conversion point in RZ;Z (cm);R_{xy} (cm)", kTH2F, {{200, -100, 100}, {200, 0.0f, 100.0f}}, false); fRegistry.add("V0/primary/hCosPA", "V0CosPA;cosine pointing angle in 3D", kTH1F, {{100, 0.99f, 1.0f}}, false); @@ -207,7 +221,7 @@ struct PCMQCMC { fRegistry.add("V0/primary/hRxyGen_DeltaR", "photon #varphi resolution;R_{xy}^{gen} (cm);#varphi^{rec} - #varphi^{gen} (rad.)", kTH2F, {{100, 0, 100}, {100, 0, 100}}, true); fRegistry.add("V0/primary/hXY_MC", "X vs. Y of true photon conversion point.;X (cm);Y (cm)", kTH2F, {{400, -100.0f, +100}, {400, -100, +100}}, true); fRegistry.add("V0/primary/hRZ_MC", "R vs. Z of true photon conversion point;Z (cm);R_{xy} (cm)", kTH2F, {{200, -100.0f, +100}, {200, 0, 100}}, true); - fRegistry.add("V0/primary/hsConvPoint", "photon conversion point;r_{xy} (cm);#varphi (rad.);#eta;", kTHnSparseF, {{100, 0.0f, 100}, {90, 0, 2 * M_PI}, {80, -2, +2}}, false); + fRegistry.add("V0/primary/hsConvPoint", "photon conversion point;r_{xy} (cm);#varphi (rad.);#eta;", kTHnSparseF, {{100, 0.0f, 100}, {90, 0, o2::constants::math::TwoPI}, {80, -2, +2}}, false); fRegistry.addClone("V0/primary/", "V0/fromWD/"); // from weak decay fRegistry.addClone("V0/primary/", "V0/fromHS/"); // from hadronic shower in detector materials fRegistry.addClone("V0/primary/", "V0/fromPi0Dalitz/"); // misidentified dielectron from pi0 dalitz decay @@ -218,7 +232,7 @@ struct PCMQCMC { // v0leg info fRegistry.add("V0Leg/primary/hPt", "pT;p_{T,e} (GeV/c)", kTH1F, {{1000, 0.0f, 10}}, false); fRegistry.add("V0Leg/primary/hQoverPt", "q/pT;q/p_{T} (GeV/c)^{-1}", kTH1F, {{1000, -50, 50}}, false); - fRegistry.add("V0Leg/primary/hEtaPhi", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{90, 0, 2 * M_PI}, {200, -1.0f, 1.0f}}, false); + fRegistry.add("V0Leg/primary/hEtaPhi", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{90, 0, o2::constants::math::TwoPI}, {200, -1.0f, 1.0f}}, false); fRegistry.add("V0Leg/primary/hDCAxyz", "DCA xy vs. z;DCA_{xy} (cm);DCA_{z} (cm)", kTH2F, {{200, -50.0f, 50.0f}, {200, -50.0f, 50.0f}}, false); fRegistry.add("V0Leg/primary/hNclsTPC", "number of TPC clusters", kTH1F, {{161, -0.5, 160.5}}, false); fRegistry.add("V0Leg/primary/hNcrTPC", "number of TPC crossed rows", kTH1F, {{161, -0.5, 160.5}}, false); @@ -448,7 +462,7 @@ struct PCMQCMC { // LOGF(info, "posmc.isPhysicalPrimary() = %d, posmc.producedByGenerator() = %d, elemc.isPhysicalPrimary() = %d, elemc.producedByGenerator() = %d", posmc.isPhysicalPrimary(), posmc.producedByGenerator(), elemc.isPhysicalPrimary(), elemc.producedByGenerator()); - if (!fV0PhotonCut.IsSelected(v0)) { + if (!fV0PhotonCut.IsSelected(v0)) { continue; } diff --git a/PWGEM/PhotonMeson/Tasks/phosQC.cxx b/PWGEM/PhotonMeson/Tasks/phosQC.cxx index 5f42145883e..4f5ef2ec075 100644 --- a/PWGEM/PhotonMeson/Tasks/phosQC.cxx +++ b/PWGEM/PhotonMeson/Tasks/phosQC.cxx @@ -19,25 +19,22 @@ #include "PWGEM/PhotonMeson/Core/PHOSPhotonCut.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" -#include "Common/Core/RecoDecay.h" -#include "Common/Core/TrackSelection.h" -#include "Common/Core/trackUtilities.h" -#include "Common/DataModel/Centrality.h" -#include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/TrackSelectionTables.h" - -#include "Framework/ASoAHelpers.h" +#include "Common/CCDB/TriggerAliases.h" + #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" #include "Framework/runDataProcessing.h" -#include "ReconstructionDataFormats/Track.h" +#include +#include +#include -#include "THashList.h" -#include "TString.h" +#include +#include -#include +#include #include #include +#include #include using namespace o2; @@ -46,7 +43,6 @@ using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::soa; using namespace o2::aod::pwgem::photon; -using std::array; using MyCollisions = soa::Join; using MyCollision = MyCollisions::iterator; @@ -147,7 +143,7 @@ struct phosQC { int ng = 0; for (auto& cluster : clusters_per_coll) { - if (cut.IsSelected(cluster)) { + if (cut.IsSelected(cluster)) { o2::aod::pwgem::photon::histogram::FillHistClass(list_cluster_cut, "", cluster); ng++; } diff --git a/PWGEM/PhotonMeson/Tasks/prefilterPhoton.cxx b/PWGEM/PhotonMeson/Tasks/prefilterPhoton.cxx index 972c805f340..cfaa5371cac 100644 --- a/PWGEM/PhotonMeson/Tasks/prefilterPhoton.cxx +++ b/PWGEM/PhotonMeson/Tasks/prefilterPhoton.cxx @@ -8,25 +8,11 @@ // 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. -// -// ======================== -// -// This code produces information on prefilter for photon. -// Please write to: daiki.sekihata@cern.ch - -#include -#include -#include -#include -#include -// #include "TString.h" -#include "Framework/ASoAHelpers.h" -#include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" +/// \file prefilterPhoton.cxx +/// \brief This code produces information on prefilter for . +/// \author D. Sekihata, daiki.sekihata@cern.ch -#include "Math/Vector4D.h" -// #include "Common/Core/RecoDecay.h" #include "PWGEM/Dilepton/Utils/PairUtilities.h" #include "PWGEM/PhotonMeson/Core/DalitzEECut.h" #include "PWGEM/PhotonMeson/Core/EMPhotonEventCut.h" @@ -34,13 +20,37 @@ #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" #include "PWGEM/PhotonMeson/Utils/PairUtilities.h" -#include "Common/Core/trackUtilities.h" - -#include "CCDB/BasicCCDBManager.h" -#include "DataFormatsParameters/GRPMagField.h" -#include "DataFormatsParameters/GRPObject.h" -#include "DetectorsBase/GeometryManager.h" -#include "DetectorsBase/Propagator.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include // IWYU pragma: keep +#include + +#include +#include +#include +#include +#include using namespace o2; using namespace o2::aod; @@ -231,12 +241,12 @@ struct prefilterPhoton { { const AxisSpec axis_mass{200, 0, 0.8, "m_{#gamma#gamma} (GeV/c^{2})"}; const AxisSpec axis_pair_pt{100, 0, 10, "p_{T,#gamma#gamma} (GeV/c)"}; - const AxisSpec axis_phiv{180, 0, M_PI, "#varphi_{V} (rad.)"}; + const AxisSpec axis_phiv{180, 0, o2::constants::math::PI, "#varphi_{V} (rad.)"}; // for pair fRegistry.add("Pair/PCMPCM/before/hMvsPt", "m_{#gamma#gamma} vs. p_{T,#gamma#gamma}", kTH2D, {axis_mass, axis_pair_pt}, true); fRegistry.add("Pair/PCMDalitzEE/before/hMvsPt", "m_{ee#gamma} vs. p_{T,ee#gamma}", kTH2D, {axis_mass, axis_pair_pt}, true); - fRegistry.add("Pair/PCMDalitzEE/before/hMvsPhiV", "m_{ee} vs. #varphi_{V}", kTH2D, {{180, 0, M_PI}, {100, 0, 0.1}}, true); + fRegistry.add("Pair/PCMDalitzEE/before/hMvsPhiV", "m_{ee} vs. #varphi_{V}", kTH2D, {{180, 0, o2::constants::math::PI}, {100, 0, 0.1}}, true); fRegistry.addClone("Pair/PCMPCM/before/", "Pair/PCMPCM/after/"); fRegistry.addClone("Pair/PCMDalitzEE/before/", "Pair/PCMDalitzEE/after/"); } @@ -353,7 +363,7 @@ struct prefilterPhoton { continue; } for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photons1_per_collision, photons2_per_collision))) { - if (!cut1.template IsSelected(g1) || !cut2.template IsSelected(g2)) { + if (!cut1.template IsSelected(g1) || !cut2.template IsSelected(g2)) { continue; } // don't apply pair cut when you produce prefilter bit. @@ -396,7 +406,7 @@ struct prefilterPhoton { } for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photons1_per_collision, photons1_per_collision))) { // PCM-PCM // cut, and subinfo is different from kPCMPCM - if (!cut1.template IsSelected(g1) || !cut1.template IsSelected(g2)) { + if (!cut1.template IsSelected(g1) || !cut1.template IsSelected(g2)) { continue; } // don't apply pair cut when you produce prefilter bit. @@ -413,7 +423,7 @@ struct prefilterPhoton { } // end of 2photon pairing loop for (const auto& g1 : photons1_per_collision) { // PCM-DalitzEE - if (!cut1.template IsSelected(g1)) { + if (!cut1.template IsSelected(g1)) { continue; } auto pos1 = g1.template posTrack_as(); @@ -501,7 +511,7 @@ struct prefilterPhoton { auto photons2_per_collision = photons2.sliceBy(perCollision2, collision.globalIndex()); for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photons1_per_collision, photons2_per_collision))) { - if (!cut1.template IsSelected(g1) || !cut2.template IsSelected(g2)) { + if (!cut1.template IsSelected(g1) || !cut2.template IsSelected(g2)) { continue; } if (map_pfb_v0[g1.globalIndex()] != 0 || map_pfb_v0[g2.globalIndex()] != 0) { @@ -531,7 +541,7 @@ struct prefilterPhoton { auto electrons_per_collision = negTracks->sliceByCached(o2::aod::emprimaryelectron::emeventId, collision.globalIndex(), cache); for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photons1_per_collision, photons1_per_collision))) { - if (!cut1.template IsSelected(g1) || !cut1.template IsSelected(g2)) { + if (!cut1.template IsSelected(g1) || !cut1.template IsSelected(g2)) { continue; } if (map_pfb_v0[g1.globalIndex()] != 0 || map_pfb_v0[g2.globalIndex()] != 0) { @@ -545,7 +555,7 @@ struct prefilterPhoton { } for (const auto& g1 : photons1_per_collision) { - if (!cut1.template IsSelected(g1)) { + if (!cut1.template IsSelected(g1)) { continue; } auto pos1 = g1.template posTrack_as(); diff --git a/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx b/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx index 729759af7b2..9431f833041 100644 --- a/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx @@ -772,7 +772,7 @@ struct TaskPi0FlowEMC { // only combine rotated photons with other photons continue; } - if (!(fEMCCut.IsSelected(photon))) { + if (!(fEMCCut.IsSelected(photon))) { continue; } if (checkEtaPhi1D(photon.eta(), RecoDecay::constrainAngle(photon.phi())) >= cfgEMCalMapLevelBackground.value) { @@ -869,7 +869,7 @@ struct TaskPi0FlowEMC { // only combine rotated photons with other photons continue; } - if (!(fEMCCut.IsSelected(photon))) { + if (!(fEMCCut.IsSelected(photon))) { continue; } if (checkEtaPhi1D(photon.eta(), RecoDecay::constrainAngle(photon.phi())) >= cfgEMCalMapLevelBackground.value) { @@ -993,7 +993,7 @@ struct TaskPi0FlowEMC { for (const auto& photon : photonsPerCollision) { registry.fill(HIST("clusterQA/hEClusterBefore"), photon.e()); // before cuts registry.fill(HIST("clusterQA/hClusterEtaPhiBefore"), photon.phi(), photon.eta()); // before cuts - if (!(fEMCCut.IsSelected(photon))) { + if (!(fEMCCut.IsSelected(photon))) { continue; } if (cfgDistanceToEdge.value && (checkEtaPhi1D(photon.eta(), RecoDecay::constrainAngle(photon.phi())) >= cfgEMCalMapLevelSameEvent.value)) { @@ -1004,7 +1004,7 @@ struct TaskPi0FlowEMC { } } for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photonsPerCollision, photonsPerCollision))) { - if (!(fEMCCut.IsSelected(g1)) || !(fEMCCut.IsSelected(g2))) { + if (!(fEMCCut.IsSelected(g1)) || !(fEMCCut.IsSelected(g2))) { continue; } @@ -1113,7 +1113,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 (!(fEMCCut.IsSelected(g1)) || !(fEMCCut.IsSelected(g2))) { continue; } // Cut edge clusters away, similar to rotation method to ensure same acceptance is used @@ -1352,14 +1352,14 @@ struct TaskPi0FlowEMC { if (emccuts.cfgEnableQA.value) { for (const auto& photon : photonsPerCollision) { 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))) { + if (!(fEMCCut.IsSelected(g1)) || !(fEMCCut.IsSelected(g2))) { continue; } @@ -1467,7 +1467,7 @@ struct TaskPi0FlowEMC { registry.fill(HIST("clusterQA/hEClusterBefore"), photon.e()); // before cuts registry.fill(HIST("clusterQA/hClusterEtaPhiBefore"), photon.phi(), photon.eta()); // before cuts } - if (!(fEMCCut.IsSelected(photon))) { + if (!(fEMCCut.IsSelected(photon))) { continue; } if (cfgDistanceToEdge.value && (checkEtaPhi1D(photon.eta(), RecoDecay::constrainAngle(photon.phi())) >= cfgEMCalMapLevelSameEvent.value)) { diff --git a/PWGEM/PhotonMeson/Tasks/testTaskEmc.cxx b/PWGEM/PhotonMeson/Tasks/testTaskEmc.cxx new file mode 100644 index 00000000000..c10b30fc3a6 --- /dev/null +++ b/PWGEM/PhotonMeson/Tasks/testTaskEmc.cxx @@ -0,0 +1,178 @@ +// 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 testTaskEmc.cxx +/// \brief Task to test values for EMCal +/// \author M. Hemmer, marvin.hemmer@cern.ch + +#include "PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h" +#include "PWGEM/PhotonMeson/DataModel/gammaTables.h" +#include "PWGEM/PhotonMeson/Utils/emcalHistoDefinitions.h" +// + +#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; + +struct TestTaskEmc { + + HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; + + void init(InitContext&) + { + const AxisSpec energyAxis{10000, 0., 10., "#it{E}_{clus} (GeV)"}; + const AxisSpec energyAxis2{10000, 0.699, 0.701, "#it{E}_{clus} (GeV)"}; + const AxisSpec m02axis{200, 0., +2, "#it{M}_{02} (a.u.)"}; + const AxisSpec phiAxis{200, 0., o2::constants::math::TwoPI, "#it{#varphi} (rad)"}; + const AxisSpec etaAxis{140, -0.7, +0.7, "#it{#eta}"}; + const AxisSpec nCellsAxis{140, -0.7, +0.7, "#it{N}_{cell}"}; + const AxisSpec isExoticAxis{2, -0.5, +1.5, "isExotic"}; + const AxisSpec timeAxis{2000, -1000., +1000., "#it{N}_{cell}"}; + + const AxisSpec deltaEnergyAxis{1000, -0.01, 0.01, "#Delta#it{E}_{clus} (GeV)"}; + const AxisSpec deltaM02axis{1000, -0.01, 0.01, "#Delta#it{M}_{02} (a.u.)"}; + const AxisSpec deltaNCellsAxis{1000, -0.01, 0.01, "#Delta#it{N}_{cell}"}; + const AxisSpec deltaTimeAxis{1000, -0.01, 0.01, "#Delta#it{N}_{cell}"}; + + const AxisSpec deltaPhiAxis{400, -0.05, +0.05, "#Delta#it{#varphi} (rad)"}; + const AxisSpec deltaEtaAxis{400, -0.05, +0.05, "#Delta#it{#eta}"}; + + registry.add("Old/hEnergy", "hEnergy", HistType::kTH1D, {energyAxis}); + registry.add("Old/hEnergy2", "hEnergy2", HistType::kTH1D, {energyAxis2}); + registry.add("Old/hM02", "hM02", HistType::kTH1D, {m02axis}); + registry.add("Old/hEta", "hEta", HistType::kTH1D, {etaAxis}); + registry.add("Old/hPhi", "hPhi", HistType::kTH1D, {phiAxis}); + registry.add("Old/hNCells", "hNCells", HistType::kTH1D, {nCellsAxis}); + registry.add("Old/hExotic", "hExotic", HistType::kTH1D, {isExoticAxis}); + registry.add("Old/hTime", "hTime", HistType::kTH1D, {timeAxis}); + + registry.add("Old/Prim/hDeltaEta", "hDeltaEta", HistType::kTH1D, {deltaEtaAxis}); + registry.add("Old/Prim/hDeltaPhi", "hDeltaPhi", HistType::kTH1D, {deltaPhiAxis}); + registry.add("Old/Prim/hP", "hP", HistType::kTH1D, {energyAxis}); + registry.add("Old/Prim/hPt", "hPt", HistType::kTH1D, {energyAxis}); + + registry.add("Old/Sec/hDeltaEta", "hDeltaEta", HistType::kTH1D, {deltaEtaAxis}); + registry.add("Old/Sec/hDeltaPhi", "hDeltaPhi", HistType::kTH1D, {deltaPhiAxis}); + registry.add("Old/Sec/hP", "hP", HistType::kTH1D, {energyAxis}); + registry.add("Old/Sec/hPt", "hPt", HistType::kTH1D, {energyAxis}); + + registry.addClone("Old/", "Redux/"); + + registry.add("hDeltaEnergy", "hDeltaEnergy", HistType::kTH1D, {deltaEnergyAxis}); + registry.add("hDeltaM02", "hDeltaM02", HistType::kTH1D, {deltaM02axis}); + registry.add("hDeltaEta", "hDeltaEta", HistType::kTH1D, {deltaEtaAxis}); + registry.add("hDeltaPhi", "hDeltaPhi", HistType::kTH1D, {deltaPhiAxis}); + registry.add("hDeltaNCells", "hDeltaNCells", HistType::kTH1D, {deltaNCellsAxis}); + registry.add("hDeltaTime", "hDeltaTime", HistType::kTH1D, {deltaTimeAxis}); + + }; // end init + + // EMCal calibration same event + void processEMCalCalib(aod::SkimEMCClusters_001 const& oldClusters, aod::MinClusters const& reduxClusters, aod::MinMTracks const& reduxMatchedTracks, aod::MinMSTracks const& reduxMatchedSecondaries) + { + auto oldCluster = oldClusters.begin(); + auto reduxCluster = reduxClusters.begin(); + + while (oldCluster != oldClusters.end() && reduxCluster != reduxClusters.end()) { + registry.fill(HIST("hDeltaEnergy"), oldCluster.e() - reduxCluster.e()); + registry.fill(HIST("hDeltaM02"), oldCluster.m02() - reduxCluster.m02()); + registry.fill(HIST("hDeltaEta"), oldCluster.eta() - reduxCluster.eta()); + registry.fill(HIST("hDeltaPhi"), oldCluster.phi() - reduxCluster.phi()); + registry.fill(HIST("hDeltaNCells"), oldCluster.nCells() - reduxCluster.nCells()); + registry.fill(HIST("hDeltaTime"), oldCluster.time() - reduxCluster.time()); + ++oldCluster; + ++reduxCluster; + } + for (const auto& cluster : oldClusters) { + registry.fill(HIST("Old/hEnergy"), cluster.e()); + registry.fill(HIST("Old/hEnergy2"), cluster.e()); + registry.fill(HIST("Old/hM02"), cluster.m02()); + registry.fill(HIST("Old/hEta"), cluster.eta()); + registry.fill(HIST("Old/hPhi"), cluster.phi()); + registry.fill(HIST("Old/hNCells"), cluster.nCells()); + registry.fill(HIST("Old/hExotic"), cluster.isExotic()); + registry.fill(HIST("Old/hTime"), cluster.time()); + + for (const auto& var : cluster.deltaEta()) { + registry.fill(HIST("Old/Prim/hDeltaEta"), var); + } + for (const auto& var : cluster.deltaPhi()) { + registry.fill(HIST("Old/Prim/hDeltaPhi"), var); + } + for (const auto& dPhi : cluster.trackp()) { + registry.fill(HIST("Old/Prim/hP"), dPhi); + } + for (const auto& var : cluster.trackpt()) { + registry.fill(HIST("Old/Prim/hPt"), var); + } + for (const auto& var : cluster.deltaEtaSec()) { + registry.fill(HIST("Old/Sec/hDeltaEta"), var); + } + for (const auto& var : cluster.deltaPhiSec()) { + registry.fill(HIST("Old/Sec/hDeltaPhi"), var); + } + for (const auto& var : cluster.trackpSec()) { + registry.fill(HIST("Old/Sec/hP"), var); + } + for (const auto& var : cluster.trackptSec()) { + registry.fill(HIST("Old/Sec/hPt"), var); + } + } + + for (const auto& cluster : reduxClusters) { + registry.fill(HIST("Redux/hEnergy"), cluster.e()); + registry.fill(HIST("Redux/hEnergy2"), cluster.e()); + registry.fill(HIST("Redux/hM02"), cluster.m02()); + registry.fill(HIST("Redux/hEta"), cluster.eta()); + registry.fill(HIST("Redux/hPhi"), cluster.phi()); + registry.fill(HIST("Redux/hNCells"), cluster.nCells()); + registry.fill(HIST("Redux/hExotic"), cluster.isExotic()); + registry.fill(HIST("Redux/hTime"), cluster.time()); + } + + for (const auto& prim : reduxMatchedTracks) { + registry.fill(HIST("Redux/Prim/hDeltaEta"), prim.deltaEta()); + registry.fill(HIST("Redux/Prim/hDeltaPhi"), prim.deltaPhi()); + registry.fill(HIST("Redux/Prim/hP"), prim.trackP()); + registry.fill(HIST("Redux/Prim/hPt"), prim.trackPt()); + } + + for (const auto& sec : reduxMatchedSecondaries) { + registry.fill(HIST("Redux/Sec/hDeltaEta"), sec.deltaEta()); + registry.fill(HIST("Redux/Sec/hDeltaPhi"), sec.deltaPhi()); + registry.fill(HIST("Redux/Sec/hP"), sec.trackP()); + registry.fill(HIST("Redux/Sec/hPt"), sec.trackPt()); + } + } + + PROCESS_SWITCH(TestTaskEmc, processEMCalCalib, "Process EMCal calibration same event", true); +}; // End struct TestTaskEmc + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} diff --git a/PWGEM/PhotonMeson/Utils/PairUtilities.h b/PWGEM/PhotonMeson/Utils/PairUtilities.h index 02b6b548dca..416d8d2e3b3 100644 --- a/PWGEM/PhotonMeson/Utils/PairUtilities.h +++ b/PWGEM/PhotonMeson/Utils/PairUtilities.h @@ -52,8 +52,8 @@ bool IsSelectedPair(TG1 const& g1, TG2 const& g2, TCut1 const& cut1, TCut2 const { bool is_g1_selected = false; bool is_g2_selected = false; - is_g1_selected = cut1.template IsSelected(g1); - is_g2_selected = cut2.template IsSelected(g2); + is_g1_selected = cut1.template IsSelected(g1); + is_g2_selected = cut2.template IsSelected(g2); return (is_g1_selected && is_g2_selected); }