From 0ad0ac484e1f9c393af4412b1a90e84180385439 Mon Sep 17 00:00:00 2001 From: Luca Barioglio Date: Mon, 19 Feb 2024 16:12:40 +0100 Subject: [PATCH 1/8] Fix last occurrence of tpcInnerParam (#4764) --- PWGLF/TableProducer/nucleiSpectra.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGLF/TableProducer/nucleiSpectra.cxx b/PWGLF/TableProducer/nucleiSpectra.cxx index 78fb208e2db..e1088ee41fd 100644 --- a/PWGLF/TableProducer/nucleiSpectra.cxx +++ b/PWGLF/TableProducer/nucleiSpectra.cxx @@ -420,7 +420,7 @@ struct nucleiSpectra { bool heliumPID = track.pidForTracking() == o2::track::PID::Helium3 || track.pidForTracking() == o2::track::PID::Alpha; float correctedTpcInnerParam = (heliumPID && cfgCompensatePIDinTracking) ? track.tpcInnerParam() / 2 : track.tpcInnerParam(); - spectra.fill(HIST("hTpcSignalData"), track.tpcInnerParam() * track.sign(), track.tpcSignal()); + spectra.fill(HIST("hTpcSignalData"), correctedTpcInnerParam * track.sign(), track.tpcSignal()); float nSigma[2][5]{ {-10., -10., -10., -10., -10.}, {0.f, 0.f, 0.f, 0.f, 0.f}}; /// then we will calibrate the TOF mass for the He3 and Alpha From 27eea4ccaca46c049db35563bf413e40872ac63a Mon Sep 17 00:00:00 2001 From: Nicolas Strangmann <77485327+nstrangm@users.noreply.github.com> Date: Mon, 19 Feb 2024 16:56:58 +0100 Subject: [PATCH 2/8] PWGEM/PhotomMeson: Add generated function for EMC in GammaGammaMC task (#4760) --- PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMC.cxx | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMC.cxx b/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMC.cxx index a44c838fe29..1c060c42173 100644 --- a/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMC.cxx @@ -610,7 +610,10 @@ struct Pi0EtaToGammaGammaMC { runGenInfo(collisions, mccollisions, mcparticles); } void processPHOSPHOS(MyCollisions const& collisions) {} - void processEMCEMC(MyCollisions const& collisions) {} + void processEMCEMC(MyCollisions const& collisions, aod::EMReducedMCEvents const& mccollisions, aod::EMMCParticles const& mcparticles) + { + runGenInfo(collisions, mccollisions, mcparticles); + } void processPCMPHOS(MyCollisions const& collisions) {} void processPCMEMC(MyCollisions const& collisions) {} From e82db61b95deb94f78ba840a39bd397746f549aa Mon Sep 17 00:00:00 2001 From: Junlee Kim Date: Tue, 20 Feb 2024 03:00:50 +0900 Subject: [PATCH 3/8] fix bug (#4758) Co-authored-by: junleekim --- Common/TableProducer/qVectorsTable.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Common/TableProducer/qVectorsTable.cxx b/Common/TableProducer/qVectorsTable.cxx index 9ad1640b78c..3e898b8bf38 100644 --- a/Common/TableProducer/qVectorsTable.cxx +++ b/Common/TableProducer/qVectorsTable.cxx @@ -356,7 +356,7 @@ struct qVectorsTable { // to ensure a proper channel number in FT0 as a whole. float ampl = ft0.amplitudeC()[iChC]; - histosQA.fill(HIST("FITAmp"), ampl, iChC); + histosQA.fill(HIST("FITAmp"), ampl, iChC + 96); helperEP.SumQvectors(0, iChC + 96, ampl, cfgnMod, QvecDet, sumAmplFT0C); helperEP.SumQvectors(0, iChC + 96, ampl, cfgnMod, QvecFT0M, sumAmplFT0M); From 3bc2fc41c782fac27fede56a874d88b90888d233 Mon Sep 17 00:00:00 2001 From: Ionut Cristian Arsene Date: Mon, 19 Feb 2024 19:20:53 +0100 Subject: [PATCH 4/8] [PWGDQ] collision association QA and additional process function in table-maker with associations (#4765) * Adding collision association QA in the table reader workflow; adding more process functions in table-maker * clang fixes --------- Co-authored-by: Ionut-Cristian Arsene --- PWGDQ/Core/HistogramsLibrary.cxx | 23 +- PWGDQ/Core/VarManager.cxx | 20 + PWGDQ/Core/VarManager.h | 54 +- PWGDQ/TableProducer/tableMaker_withAssoc.cxx | 150 ++++-- PWGDQ/Tasks/tableReader.cxx | 2 +- PWGDQ/Tasks/tableReader_withAssoc.cxx | 531 +++++++++++++------ 6 files changed, 553 insertions(+), 227 deletions(-) diff --git a/PWGDQ/Core/HistogramsLibrary.cxx b/PWGDQ/Core/HistogramsLibrary.cxx index b4f66474454..5ab4e8a1cf9 100644 --- a/PWGDQ/Core/HistogramsLibrary.cxx +++ b/PWGDQ/Core/HistogramsLibrary.cxx @@ -50,9 +50,9 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h } } if (subGroupStr.Contains("vtx")) { - hm->AddHistogram(histClass, "VtxX", "Vtx X", false, 100, -0.5, 0.5, VarManager::kVtxX); - hm->AddHistogram(histClass, "VtxY", "Vtx Y", false, 100, -0.5, 0.5, VarManager::kVtxY); - hm->AddHistogram(histClass, "VtxYVtxX", "Vtx Y vs Vtx X", false, 50, -0.5, 0.5, VarManager::kVtxX, 50, -0.5, 0.5, VarManager::kVtxY); + hm->AddHistogram(histClass, "VtxX", "Vtx X", false, 200, -0.1, 0.1, VarManager::kVtxX); + hm->AddHistogram(histClass, "VtxY", "Vtx Y", false, 200, -0.1, 0.1, VarManager::kVtxY); + hm->AddHistogram(histClass, "VtxYVtxX", "Vtx Y vs Vtx X", false, 100, -0.1, 0.1, VarManager::kVtxX, 100, -0.1, 0.1, VarManager::kVtxY); } if (subGroupStr.Contains("vtxpp")) { hm->AddHistogram(histClass, "VtxNContrib", "Vtx n contributors", false, 100, 0.0, 100.0, VarManager::kVtxNcontrib); @@ -173,6 +173,17 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h hm->AddHistogram(histClass, "IsSingleGapA", "Is single gap on side A", false, 2, -0.5, 1.5, VarManager::kIsSingleGapA); hm->AddHistogram(histClass, "IsSingleGapC", "Is single gap on side C", false, 2, -0.5, 1.5, VarManager::kIsSingleGapC); } + } // end "event" + + if (groupStr.CompareTo("two-collisions") == 0) { + hm->AddHistogram(histClass, "DeltaZ", "z_{1} - z_{2}", false, 400, -20., 20., VarManager::kTwoEvDeltaZ); + hm->AddHistogram(histClass, "DeltaZ_Z1", "z_{1} - z_{2} vs z_{1}", false, 24, -12., 12., VarManager::kTwoEvPosZ1, 300, -15., 15., VarManager::kTwoEvDeltaZ); + hm->AddHistogram(histClass, "DeltaR", "r_{1} - r_{2}", false, 200, -0.1, 0.1, VarManager::kTwoEvDeltaR); + hm->AddHistogram(histClass, "DeltaX", "x_{1} - x_{2}", false, 200, -0.1, 0.1, VarManager::kTwoEvDeltaX); + hm->AddHistogram(histClass, "DeltaY", "y_{1} - y_{2}", false, 200, -0.1, 0.1, VarManager::kTwoEvDeltaY); + hm->AddHistogram(histClass, "Z1vsZ2", "z_{1} vs z_{2}", false, 120, -12.0, 12.0, VarManager::kTwoEvPosZ1, 120, -12.0, 12.0, VarManager::kTwoEvPosZ2); + hm->AddHistogram(histClass, "R1vsR2", "r_{1} vs r_{2}", false, 100, -0.1, 0.1, VarManager::kTwoEvPosR1, 100, -0.1, 0.1, VarManager::kTwoEvPosR2); + hm->AddHistogram(histClass, "NContrib1vs2", "n.contrib 1 vs 2", false, 100, 0.0, 100.0, VarManager::kTwoEvPVcontrib1, 100, 0.0, 100.0, VarManager::kTwoEvPVcontrib2); } if (groupStr.Contains("track")) { @@ -503,6 +514,12 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h hm->AddHistogram(histClass, "Mass", "", false, 500, 0.0, 5.0, VarManager::kMass); hm->AddHistogram(histClass, "Rapidity", "", false, 400, -4.0, 4.0, VarManager::kRap); } + if (subGroupStr.Contains("ambiguity")) { + hm->AddHistogram(histClass, "AmbiguityInBunch", "in bunch collision ambiguity", false, 10, 0., 10., VarManager::kBarrelNAssocsInBunch); + hm->AddHistogram(histClass, "AmbiguityInBunch_pt", "in bunch collision ambiguity vs p_{T}", false, 50, 0.0, 10.0, VarManager::kPt, 10, 0., 10., VarManager::kBarrelNAssocsInBunch); + hm->AddHistogram(histClass, "AmbiguityOutOfBunch", "out of bunch collision ambiguity", false, 10, 0., 10., VarManager::kBarrelNAssocsOutOfBunch); + hm->AddHistogram(histClass, "AmbiguityOutOfBunch_pt", "out of bunch collision ambiguity vs p_{T}", false, 50, 0.0, 10.0, VarManager::kPt, 10, 0., 10., VarManager::kBarrelNAssocsOutOfBunch); + } } if (groupStr.Contains("mctruth_pair")) { hm->AddHistogram(histClass, "Mass_Pt", "", false, 500, 0.0, 5.0, VarManager::kMass, 200, 0.0, 20.0, VarManager::kPt); diff --git a/PWGDQ/Core/VarManager.cxx b/PWGDQ/Core/VarManager.cxx index 9bd573cb408..2f6c3b654e5 100644 --- a/PWGDQ/Core/VarManager.cxx +++ b/PWGDQ/Core/VarManager.cxx @@ -325,6 +325,26 @@ void VarManager::SetDefaultVarNames() fgVariableUnits[kMCEventTime] = ""; // TODO: add proper unit fgVariableUnits[kMCEventWeight] = ""; fgVariableUnits[kMCEventImpParam] = "b"; + fgVariableNames[kTwoEvPosZ1] = "vtx-z_{1}"; + fgVariableUnits[kTwoEvPosZ1] = "cm"; + fgVariableNames[kTwoEvPosZ2] = "vtx-z_{2}"; + fgVariableUnits[kTwoEvPosZ2] = "cm"; + fgVariableNames[kTwoEvPosR1] = "vtx-R_{1}"; + fgVariableUnits[kTwoEvPosR1] = "cm"; + fgVariableNames[kTwoEvPosR2] = "vtx-R_{2}"; + fgVariableUnits[kTwoEvPosR2] = "cm"; + fgVariableNames[kTwoEvDeltaZ] = "#Delta_{z}"; + fgVariableUnits[kTwoEvDeltaZ] = "cm"; + fgVariableNames[kTwoEvDeltaR] = "#Delta_{R}"; + fgVariableUnits[kTwoEvDeltaR] = "cm"; + fgVariableNames[kTwoEvDeltaX] = "#Delta_{x}"; + fgVariableUnits[kTwoEvDeltaX] = "cm"; + fgVariableNames[kTwoEvDeltaY] = "#Delta_{y}"; + fgVariableUnits[kTwoEvDeltaY] = "cm"; + fgVariableNames[kTwoEvPVcontrib1] = "n.contrib 1"; + fgVariableUnits[kTwoEvPVcontrib1] = ""; + fgVariableNames[kTwoEvPVcontrib2] = "n.contrib 2"; + fgVariableUnits[kTwoEvPVcontrib2] = ""; fgVariableNames[kPt] = "p_{T}"; fgVariableUnits[kPt] = "GeV/c"; fgVariableNames[kInvPt] = "1/p_{T}"; diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index fe7f068b700..f78a560337b 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -60,6 +60,9 @@ #include "KFParticleBase.h" #include "KFVertex.h" +using std::cout; +using std::endl; + using SMatrix55 = ROOT::Math::SMatrix>; using SMatrix5 = ROOT::Math::SVector; using Vec3D = ROOT::Math::SVector; @@ -215,6 +218,16 @@ class VarManager : public TObject kIsSingleGapA, // Rapidity gap on side A kIsSingleGapC, // Rapidity gap on side C kIsSingleGap, // Rapidity gap on either side + kTwoEvPosZ1, // vtx-z for collision 1 in two events correlations + kTwoEvPosZ2, // vtx-z for collision 2 in two events correlations + kTwoEvPosR1, // vtx-R for collision 1 in two events correlations + kTwoEvPosR2, + kTwoEvPVcontrib1, // n-contributors for collision 1 in two events correlations + kTwoEvPVcontrib2, + kTwoEvDeltaZ, // distance in z between collisions + kTwoEvDeltaX, // distance in x between collisions + kTwoEvDeltaY, // distance in y between collisions + kTwoEvDeltaR, // distance in (x,y) plane between collisions kNEventWiseVariables, // Basic track/muon/pair wise variables @@ -326,7 +339,9 @@ class VarManager : public TObject kIsLegFromAntiLambda, kIsLegFromOmega, kIsProtonFromLambdaAndAntiLambda, - kIsDalitzLeg, // Up to 8 dalitz selections + kIsDalitzLeg, // Up to 8 dalitz selections + kBarrelNAssocsInBunch, // number of in bunch collision associations + kBarrelNAssocsOutOfBunch, // number of out of bunch collision associations kNBarrelTrackVariables = kIsDalitzLeg + 8, // Muon track variables @@ -642,6 +657,8 @@ class VarManager : public TObject static void FillPropagateMuon(const T& muon, const C& collision, float* values = nullptr); template static void FillEvent(T const& event, float* values = nullptr); + template + static void FillTwoEvents(T const& event1, T const& event2, float* values = nullptr); template static void FillTrack(T const& track, float* values = nullptr); template @@ -852,14 +869,14 @@ o2::dataformats::GlobalFwdTrack VarManager::PropagateMuon(const T& muon, const C propmuon.setParameters(proptrack.getParameters()); propmuon.setZ(proptrack.getZ()); propmuon.setCovariances(proptrack.getCovariances()); - } else if (static_cast(muon.trackType()) < 2) { - double centerMFT[3] = {0, 0, -61.4}; + /*double centerMFT[3] = {0, 0, -61.4}; o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); auto Bz = field->getBz(centerMFT); // Get field at centre of MFT + //auto Bz = fgMagField; auto geoMan = o2::base::GeometryManager::meanMaterialBudget(muon.x(), muon.y(), muon.z(), collision.posX(), collision.posY(), collision.posZ()); auto x2x0 = static_cast(geoMan.meanX2X0); - fwdtrack.propagateToVtxhelixWithMCS(collision.posZ(), {collision.posX(), collision.posY()}, {collision.covXX(), collision.covYY()}, Bz, x2x0); + fwdtrack.propagateToVtxhelixWithMCS(collision.posZ(), {collision.posX(), collision.posY()}, {collision.covXX(), collision.covYY()}, Bz, x2x0);*/ propmuon.setParameters(fwdtrack.getParameters()); propmuon.setZ(fwdtrack.getZ()); propmuon.setCovariances(fwdtrack.getCovariances()); @@ -873,7 +890,6 @@ void VarManager::FillPropagateMuon(const T& muon, const C& collision, float* val if (!values) { values = fgValues; } - if constexpr ((fillMap & MuonCov) > 0 || (fillMap & ReducedMuonCov) > 0) { o2::dataformats::GlobalFwdTrack propmuon = PropagateMuon(muon, collision); @@ -1133,6 +1149,34 @@ void VarManager::FillEvent(T const& event, float* values) FillEventDerived(values); } +template +void VarManager::FillTwoEvents(T const& ev1, T const& ev2, float* values) +{ + if (!values) { + values = fgValues; + } + + values[kTwoEvPosZ1] = ev1.posZ(); + values[kTwoEvPosZ2] = ev2.posZ(); + values[kTwoEvPosR1] = std::sqrt(ev1.posX() * ev1.posX() + ev1.posY() * ev1.posY()); + values[kTwoEvPosR2] = std::sqrt(ev2.posX() * ev2.posX() + ev2.posY() * ev2.posY()); + values[kTwoEvPVcontrib1] = ev1.numContrib(); + values[kTwoEvPVcontrib2] = ev2.numContrib(); + if (ev1.numContrib() < ev2.numContrib()) { + values[kTwoEvPosZ1] = ev2.posZ(); + values[kTwoEvPosZ2] = ev1.posZ(); + values[kTwoEvPVcontrib1] = ev2.numContrib(); + values[kTwoEvPVcontrib2] = ev1.numContrib(); + values[kTwoEvPosR1] = std::sqrt(ev2.posX() * ev2.posX() + ev2.posY() * ev2.posY()); + ; + values[kTwoEvPosR2] = std::sqrt(ev1.posX() * ev1.posX() + ev1.posY() * ev1.posY()); + } + values[kTwoEvDeltaZ] = ev1.posZ() - ev2.posZ(); + values[kTwoEvDeltaX] = ev1.posX() - ev2.posX(); + values[kTwoEvDeltaY] = ev1.posY() - ev2.posY(); + values[kTwoEvDeltaR] = std::sqrt(values[kTwoEvDeltaX] * values[kTwoEvDeltaX] + values[kTwoEvDeltaY] * values[kTwoEvDeltaY]); +} + template void VarManager::FillTrack(T const& track, float* values) { diff --git a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx index 56e32544e7b..e0146b0695b 100644 --- a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx @@ -114,7 +114,7 @@ DECLARE_SOA_TABLE(AmbiguousTracksFwd, "AOD", "AMBIGUOUSFWDTR", //! Table for Fwd // constexpr static uint32_t gkEventFillMapWithMult = VarManager::ObjTypes::BC | VarManager::ObjTypes::Collision | VarManager::ObjTypes::CollisionMult; constexpr static uint32_t gkEventFillMapWithMultsAndEventFilter = VarManager::ObjTypes::BC | VarManager::ObjTypes::Collision | VarManager::ObjTypes::CollisionMult | VarManager::ObjTypes::EventFilter; // constexpr static uint32_t gkEventFillMapWithCent = VarManager::ObjTypes::BC | VarManager::ObjTypes::Collision | VarManager::ObjTypes::CollisionCent; -// constexpr static uint32_t gkEventFillMapWithCentAndMults = VarManager::ObjTypes::BC | VarManager::ObjTypes::Collision | VarManager::ObjTypes::CollisionCent | VarManager::CollisionMult; +constexpr static uint32_t gkEventFillMapWithCentAndMults = VarManager::ObjTypes::BC | VarManager::ObjTypes::Collision | VarManager::ObjTypes::CollisionCent | VarManager::CollisionMult; // constexpr static uint32_t gkEventFillMapWithCentRun2 = VarManager::ObjTypes::BC | VarManager::ObjTypes::Collision | VarManager::ObjTypes::CollisionCentRun2; // Unused variable // constexpr static uint32_t gkTrackFillMap = VarManager::ObjTypes::Track | VarManager::ObjTypes::TrackExtra | VarManager::ObjTypes::TrackDCA | VarManager::ObjTypes::TrackPID | VarManager::ObjTypes::TrackPIDExtra; constexpr static uint32_t gkTrackFillMapWithCov = VarManager::ObjTypes::Track | VarManager::ObjTypes::TrackExtra | VarManager::ObjTypes::TrackDCA | VarManager::ObjTypes::TrackCov | VarManager::ObjTypes::TrackPID | VarManager::ObjTypes::TrackPIDExtra; @@ -257,29 +257,10 @@ struct TableMaker { histClasses += "Event_AfterCuts;"; } - /* - bool enableBarrelHistos = (context.mOptions.get("processFull") || context.mOptions.get("processFullWithCov") || - context.mOptions.get("processFullWithCent") || context.mOptions.get("processFullWithCovAndEventFilter") || - context.mOptions.get("processFullWithCovMultsAndEventFilter") || - context.mOptions.get("processBarrelOnly") || context.mOptions.get("processBarrelOnlyWithCent") || context.mOptions.get("processBarrelOnlyWithCovWithCent") || - context.mOptions.get("processBarrelOnlyWithMults") || context.mOptions.get("processBarrelOnlyWithCentAndMults") || context.mOptions.get("processBarrelOnlyWithCovWithCentAndMults") || - context.mOptions.get("processBarrelOnlyWithCov") || context.mOptions.get("processBarrelOnlyWithEventFilter") || - context.mOptions.get("processBarrelOnlyWithMultsAndEventFilter") || context.mOptions.get("processBarrelOnlyWithCovAndEventFilter") || - context.mOptions.get("processBarrelOnlyWithDalitzBits") || context.mOptions.get("processBarrelOnlyWithV0Bits") || - context.mOptions.get("processBarrelOnlyWithV0BitsAndMaps") || context.mOptions.get("processAmbiguousBarrelOnly")); - bool enableMuonHistos = (context.mOptions.get("processFull") || context.mOptions.get("processFullWithCov") || - context.mOptions.get("processFullWithCent") || context.mOptions.get("processFullWithCovAndEventFilter") || - context.mOptions.get("processFullWithCovMultsAndEventFilter") || - context.mOptions.get("processMuonOnly") || context.mOptions.get("processMuonOnlyWithCent") || - context.mOptions.get("processMuonOnlyWithMults") || context.mOptions.get("processMuonOnlyWithCentAndMults") || - context.mOptions.get("processMuonOnlyWithCovAndCent") || - context.mOptions.get("processMuonOnlyWithCov") || context.mOptions.get("processMuonOnlyWithCovAndEventFilter") || context.mOptions.get("processMuonOnlyWithEventFilter") || - context.mOptions.get("processMuonOnlyWithMultsAndEventFilter") || context.mOptions.get("processAmbiguousMuonOnlyWithCov") || context.mOptions.get("processAmbiguousMuonOnly") || - context.mOptions.get("processMuonsAndMFT") || context.mOptions.get("processMuonsAndMFTWithFilter") || context.mOptions.get("processMuonMLOnly")); - */ - - bool enableBarrelHistos = (context.mOptions.get("processFullWithFilter")); - bool enableMuonHistos = (context.mOptions.get("processFullWithFilter")); + bool enableBarrelHistos = (context.mOptions.get("processPPWithFilter") || context.mOptions.get("processPPWithFilterBarrelOnly") || + context.mOptions.get("processPbPb") || context.mOptions.get("processPbPbBarrelOnly")); + bool enableMuonHistos = (context.mOptions.get("processPPWithFilter") || context.mOptions.get("processPPWithFilterMuonOnly") || context.mOptions.get("processPPWithFilterMuonMFT") || + context.mOptions.get("processPbPb") || context.mOptions.get("processPbPbMuonOnly") || context.mOptions.get("processPbPbMuonMFT")); if (enableBarrelHistos) { if (fDoDetailedQA) { @@ -412,15 +393,15 @@ struct TableMaker { fStatsList.setObject(new TList()); fStatsList->SetOwner(kTRUE); std::vector eventLabels{"BCs", "Collisions before filtering", "Before cuts", "After cuts"}; - TH2I* histEvents = new TH2I("EventStats", "Event statistics", eventLabels.size(), -0.5, eventLabels.size() - 0.5, kNaliases + 1, -0.5, kNaliases + 0.5); + TH2I* histEvents = new TH2I("EventStats", "Event statistics", eventLabels.size(), -0.5, eventLabels.size() - 0.5, o2::aod::evsel::kNsel + 1, -0.5, double(o2::aod::evsel::kNsel) + 0.5); int ib = 1; for (auto label = eventLabels.begin(); label != eventLabels.end(); label++, ib++) { histEvents->GetXaxis()->SetBinLabel(ib, (*label).Data()); } - for (int ib = 1; ib <= kNaliases; ib++) { - histEvents->GetYaxis()->SetBinLabel(ib, aliasLabels[ib - 1].data()); + for (int ib = 1; ib <= o2::aod::evsel::kNsel; ib++) { + histEvents->GetYaxis()->SetBinLabel(ib, o2::aod::evsel::selectionLabels[ib - 1]); } - histEvents->GetYaxis()->SetBinLabel(kNaliases + 1, "Total"); + histEvents->GetYaxis()->SetBinLabel(o2::aod::evsel::kNsel + 1, "Total"); fStatsList->Add(histEvents); // Track statistics: one bin for each track selection and 5 bins for V0 tags (gamma, K0s, Lambda, anti-Lambda, Omega) @@ -465,12 +446,12 @@ struct TableMaker { for (const auto& collision : collisions) { - for (int i = 0; i < kNaliases; i++) { - if (collision.alias_bit(i) > 0) { + for (int i = 0; i < o2::aod::evsel::kNsel; i++) { + if (collision.selection_bit(i)) { (reinterpret_cast(fStatsList->At(0)))->Fill(1.0, static_cast(i)); } } - (reinterpret_cast(fStatsList->At(0)))->Fill(1.0, static_cast(kNaliases)); + (reinterpret_cast(fStatsList->At(0)))->Fill(1.0, static_cast(o2::aod::evsel::kNsel)); // apply the event filter computed by filter-PP if constexpr ((TEventFillMap & VarManager::ObjTypes::EventFilter) > 0) { @@ -505,26 +486,25 @@ struct TableMaker { fHistMan->FillHistClass("Event_BeforeCuts", VarManager::fgValues); } - uint32_t triggerAliases = collision.alias_raw(); // fill stats information, before selections - for (int i = 0; i < kNaliases; i++) { - if (triggerAliases & (uint32_t(1) << i)) { + for (int i = 0; i < o2::aod::evsel::kNsel; i++) { + if (collision.selection_bit(i)) { (reinterpret_cast(fStatsList->At(0)))->Fill(2.0, static_cast(i)); } } - (reinterpret_cast(fStatsList->At(0)))->Fill(2.0, static_cast(kNaliases)); + (reinterpret_cast(fStatsList->At(0)))->Fill(2.0, static_cast(o2::aod::evsel::kNsel)); if (!fEventCut->IsSelected(VarManager::fgValues)) { continue; } // fill stats information, after selections - for (int i = 0; i < kNaliases; i++) { - if (triggerAliases & (uint32_t(1) << i)) { + for (int i = 0; i < o2::aod::evsel::kNsel; i++) { + if (collision.selection_bit(i)) { (reinterpret_cast(fStatsList->At(0)))->Fill(3.0, static_cast(i)); } } - (reinterpret_cast(fStatsList->At(0)))->Fill(3.0, static_cast(kNaliases)); + (reinterpret_cast(fStatsList->At(0)))->Fill(3.0, static_cast(o2::aod::evsel::kNsel)); fHistMan->FillHistClass("Event_AfterCuts", VarManager::fgValues); @@ -793,10 +773,11 @@ struct TableMaker { } // end skimMuons // Produce standard barrel + muon tables with event filter (typically for pp and p-Pb) ------------------------------------------------------ - template + template void fullSkimming(TEvents const& collisions, BCsWithTimestamps const& bcs, - TTracks const& tracksBarrel, TMuons const& muons, MFTTracks const& mftTracks, - TrackAssoc const& trackAssocs, FwdTrackAssoc const& fwdTrackAssocs, MFTTrackAssoc const& mftAssocs) + TTracks const& tracksBarrel, TMuons const& muons, TMFTTracks const& mftTracks, + TTrackAssoc const& trackAssocs, TFwdTrackAssoc const& fwdTrackAssocs, TMFTTrackAssoc const& mftAssocs) { if (fCurrentRun != bcs.begin().runNumber()) { @@ -879,28 +860,95 @@ struct TableMaker { } // end loop over skimmed collisions } - // - void processFullWithFilter(MyEventsWithMultsAndFilter const& collisions, BCsWithTimestamps const& bcs, - MyBarrelTracksWithCov const& tracksBarrel, - MyMuonsWithCov const& muons, MFTTracks const& mftTracks, - TrackAssoc const& trackAssocs, FwdTrackAssoc const& fwdTrackAssocs, - MFTTrackAssoc const& mftAssocs) + // produce the full DQ skimmed data model typically for pp/p-Pb or UPC Pb-Pb (no centrality), subscribe to the DQ event filter (filter-pp or filter-PbPb) + void processPPWithFilter(MyEventsWithMultsAndFilter const& collisions, BCsWithTimestamps const& bcs, + MyBarrelTracksWithCov const& tracksBarrel, + MyMuonsWithCov const& muons, MFTTracks const& mftTracks, + TrackAssoc const& trackAssocs, FwdTrackAssoc const& fwdTrackAssocs, + MFTTrackAssoc const& mftAssocs) { fullSkimming(collisions, bcs, tracksBarrel, muons, mftTracks, trackAssocs, fwdTrackAssocs, mftAssocs); } + // produce the barrel-only DQ skimmed data model typically for pp/p-Pb or UPC Pb-Pb (no centrality), subscribe to the DQ event filter (filter-pp or filter-PbPb) + void processPPWithFilterBarrelOnly(MyEventsWithMultsAndFilter const& collisions, BCsWithTimestamps const& bcs, + MyBarrelTracksWithCov const& tracksBarrel, + TrackAssoc const& trackAssocs) + { + /*const int& a = 0; + MFTTracks const& mftTracks = 0; + FwdTrackAssoc const& fwdTrackAssocs = 0; + MFTTrackAssoc const& mftAssocs = 0;*/ + fullSkimming(collisions, bcs, tracksBarrel, nullptr, nullptr, trackAssocs, nullptr, nullptr); + } + + // produce the muon-only DQ skimmed data model typically for pp/p-Pb or UPC Pb-Pb (no centrality), subscribe to the DQ event filter (filter-pp or filter-PbPb) + void processPPWithFilterMuonOnly(MyEventsWithMultsAndFilter const& collisions, BCsWithTimestamps const& bcs, + MyMuonsWithCov const& muons, FwdTrackAssoc const& fwdTrackAssocs) + { + fullSkimming(collisions, bcs, nullptr, muons, nullptr, nullptr, fwdTrackAssocs, nullptr); + } + + // produce the muon+mft DQ skimmed data model typically for pp/p-Pb or UPC Pb-Pb (no centrality), subscribe to the DQ event filter (filter-pp or filter-PbPb) + void processPPWithFilterMuonMFT(MyEventsWithMultsAndFilter const& collisions, BCsWithTimestamps const& bcs, + MyMuonsWithCov const& muons, MFTTracks const& mftTracks, + FwdTrackAssoc const& fwdTrackAssocs, MFTTrackAssoc const& mftAssocs) + { + fullSkimming(collisions, bcs, nullptr, muons, mftTracks, nullptr, fwdTrackAssocs, mftAssocs); + } + + // produce the full DQ skimmed data model typically for Pb-Pb (with centrality), no subscribtion to the DQ event filter + void processPbPb(MyEventsWithCentAndMults const& collisions, BCsWithTimestamps const& bcs, + MyBarrelTracksWithCov const& tracksBarrel, + MyMuonsWithCov const& muons, MFTTracks const& mftTracks, + TrackAssoc const& trackAssocs, FwdTrackAssoc const& fwdTrackAssocs, + MFTTrackAssoc const& mftAssocs) + { + fullSkimming(collisions, bcs, tracksBarrel, muons, mftTracks, trackAssocs, fwdTrackAssocs, mftAssocs); + } + + // produce the barrel only DQ skimmed data model typically for Pb-Pb (with centrality), no subscribtion to the DQ event filter + void processPbPbBarrelOnly(MyEventsWithCentAndMults const& collisions, BCsWithTimestamps const& bcs, + MyBarrelTracksWithCov const& tracksBarrel, + TrackAssoc const& trackAssocs) + { + fullSkimming(collisions, bcs, tracksBarrel, nullptr, nullptr, trackAssocs, nullptr, nullptr); + } + + // produce the muon only DQ skimmed data model typically for Pb-Pb (with centrality), no subscribtion to the DQ event filter + void processPbPbMuonOnly(MyEventsWithCentAndMults const& collisions, BCsWithTimestamps const& bcs, + MyMuonsWithCov const& muons, FwdTrackAssoc const& fwdTrackAssocs) + { + fullSkimming(collisions, bcs, nullptr, muons, nullptr, nullptr, fwdTrackAssocs, nullptr); + } + + // produce the muon+mft DQ skimmed data model typically for Pb-Pb (with centrality), no subscribtion to the DQ event filter + void processPbPbMuonMFT(MyEventsWithCentAndMults const& collisions, BCsWithTimestamps const& bcs, + MyMuonsWithCov const& muons, MFTTracks const& mftTracks, + FwdTrackAssoc const& fwdTrackAssocs, MFTTrackAssoc const& mftAssocs) + { + fullSkimming(collisions, bcs, nullptr, muons, mftTracks, nullptr, fwdTrackAssocs, mftAssocs); + } + // Process the BCs and store stats for luminosity retrieval ----------------------------------------------------------------------------------- void processOnlyBCs(soa::Join::iterator const& bc) { - for (int i = 0; i < kNaliases; i++) { - if (bc.alias_bit(i) > 0) { + for (int i = 0; i < o2::aod::evsel::kNsel; i++) { + if (bc.selection_bit(i) > 0) { (reinterpret_cast(fStatsList->At(0)))->Fill(0.0, static_cast(i)); } } - (reinterpret_cast(fStatsList->At(0)))->Fill(0.0, static_cast(kNaliases)); + (reinterpret_cast(fStatsList->At(0)))->Fill(0.0, static_cast(o2::aod::evsel::kNsel)); } - PROCESS_SWITCH(TableMaker, processFullWithFilter, "Build full DQ skimmed data model, w/ event filtering", false); + PROCESS_SWITCH(TableMaker, processPPWithFilter, "Build full DQ skimmed data model typically for pp/p-Pb and UPC Pb-Pb, w/ event filtering", false); + PROCESS_SWITCH(TableMaker, processPPWithFilterBarrelOnly, "Build barrel only DQ skimmed data model typically for pp/p-Pb and UPC Pb-Pb, w/ event filtering", false); + PROCESS_SWITCH(TableMaker, processPPWithFilterMuonOnly, "Build muon only DQ skimmed data model typically for pp/p-Pb and UPC Pb-Pb, w/ event filtering", false); + PROCESS_SWITCH(TableMaker, processPPWithFilterMuonMFT, "Build muon + mft DQ skimmed data model typically for pp/p-Pb and UPC Pb-Pb, w/ event filtering", false); + PROCESS_SWITCH(TableMaker, processPbPb, "Build full DQ skimmed data model typically for Pb-Pb, w/o event filtering", false); + PROCESS_SWITCH(TableMaker, processPbPbBarrelOnly, "Build barrel only DQ skimmed data model typically for Pb-Pb, w/o event filtering", false); + PROCESS_SWITCH(TableMaker, processPbPbMuonOnly, "Build muon only DQ skimmed data model typically for Pb-Pb, w/o event filtering", false); + PROCESS_SWITCH(TableMaker, processPbPbMuonMFT, "Build muon + mft DQ skimmed data model typically for Pb-Pb, w/o event filtering", false); PROCESS_SWITCH(TableMaker, processOnlyBCs, "Analyze the BCs to store sampled lumi", false); }; diff --git a/PWGDQ/Tasks/tableReader.cxx b/PWGDQ/Tasks/tableReader.cxx index c8efd3d31fe..eda688d8fba 100644 --- a/PWGDQ/Tasks/tableReader.cxx +++ b/PWGDQ/Tasks/tableReader.cxx @@ -1359,7 +1359,7 @@ struct AnalysisDileptonHadron { Produces BmesonsTable; Filter eventFilter = aod::dqanalysisflags::isEventSelected == 1; - Filter dileptonFilter = aod::reducedpair::pt > fConfigDileptonpTCut.value&& aod::reducedpair::mass > fConfigDileptonLowMass.value&& aod::reducedpair::mass < fConfigDileptonHighMass.value&& aod::reducedpair::sign == 0; + Filter dileptonFilter = aod::reducedpair::pt > fConfigDileptonpTCut&& aod::reducedpair::mass > fConfigDileptonLowMass&& aod::reducedpair::mass < fConfigDileptonHighMass&& aod::reducedpair::sign == 0; Filter filterBarrelTrackSelected = aod::dqanalysisflags::isBarrelSelected > 0; constexpr static uint32_t fgDileptonFillMap = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::Pair; // fill map diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index 0d0521d640c..7901ad55643 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -57,8 +57,10 @@ namespace o2::aod namespace dqanalysisflags { DECLARE_SOA_COLUMN(MixingHash, mixingHash, int); //! Hash used in event mixing -DECLARE_SOA_COLUMN(IsEventSelected, isEventSelected, int); //! Event decision +DECLARE_SOA_BITMAP_COLUMN(IsEventSelected, isEventSelected, 32); //! Event decision DECLARE_SOA_BITMAP_COLUMN(IsBarrelSelected, isBarrelSelected, 32); //! Barrel track decisions +DECLARE_SOA_COLUMN(BarrelAmbiguityInBunch, barrelAmbiguityInBunch, int8_t); //! Barrel track in-bunch ambiguity +DECLARE_SOA_COLUMN(BarrelAmbiguityOutOfBunch, barrelAmbiguityOutOfBunch, int8_t); //! Barrel track out of bunch ambiguity DECLARE_SOA_BITMAP_COLUMN(IsMuonSelected, isMuonSelected, 32); //! Muon track decisions (joinable to ReducedMuonsAssoc) DECLARE_SOA_BITMAP_COLUMN(IsBarrelSelectedPrefilter, isBarrelSelectedPrefilter, 32); //! Barrel prefilter decisions (joinable to ReducedTracksAssoc) // Bcandidate columns for ML analysis of B->Jpsi+K @@ -76,6 +78,7 @@ DECLARE_SOA_COLUMN(Chi2Bcandidate, chi2Bcandidate, float); DECLARE_SOA_TABLE(EventCuts, "AOD", "DQANAEVCUTS", dqanalysisflags::IsEventSelected); //! joinable to ReducedEvents DECLARE_SOA_TABLE(MixingHashes, "AOD", "DQANAMIXHASH", dqanalysisflags::MixingHash); //! joinable to ReducedEvents DECLARE_SOA_TABLE(BarrelTrackCuts, "AOD", "DQANATRKCUTS", dqanalysisflags::IsBarrelSelected); //! joinable to ReducedTracksAssoc +DECLARE_SOA_TABLE(BarrelAmbiguities, "AOD", "DQBARRELAMB", dqanalysisflags::BarrelAmbiguityInBunch, dqanalysisflags::BarrelAmbiguityOutOfBunch); //! joinable to ReducedBarrelTracks DECLARE_SOA_TABLE(MuonTrackCuts, "AOD", "DQANAMUONCUTS", dqanalysisflags::IsMuonSelected); //! joinable to ReducedMuonsAssoc DECLARE_SOA_TABLE(Prefilter, "AOD", "DQPREFILTER", dqanalysisflags::IsBarrelSelectedPrefilter); //! joinable to ReducedTracksAssoc DECLARE_SOA_TABLE(BmesonCandidates, "AOD", "DQBMESONS", dqanalysisflags::massBcandidate, dqanalysisflags::pTBcandidate, dqanalysisflags::LxyBcandidate, dqanalysisflags::LxyzBcandidate, dqanalysisflags::LzBcandidate, dqanalysisflags::TauxyBcandidate, dqanalysisflags::TauzBcandidate, dqanalysisflags::CosPBcandidate, dqanalysisflags::Chi2Bcandidate); @@ -92,7 +95,9 @@ using MyEventsQvector = soa::Join; using MyBarrelTracks = soa::Join; +using MyBarrelTracksWithAmbiguities = soa::Join; using MyBarrelTracksWithCov = soa::Join; +using MyBarrelTracksWithCovWithAmbiguities = soa::Join; using MyDielectronCandidates = soa::Join; using MyDimuonCandidates = soa::Join; using MyMuonTracks = soa::Join; @@ -142,6 +147,9 @@ struct AnalysisEventSelection { MixingHandler* fMixHandler = nullptr; AnalysisCompositeCut* fEventCut; + std::map fSelMap; // key: reduced event global index, value: event selection decision + std::map> fBCCollMap; // key: global BC, value: vector of reduced event global indices + void init(o2::framework::InitContext&) { fEventCut = new AnalysisCompositeCut(true); @@ -154,7 +162,7 @@ struct AnalysisEventSelection { fHistMan = new HistogramManager("analysisHistos", "", VarManager::kNVars); fHistMan->SetUseDefaultVariableNames(kTRUE); fHistMan->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); - DefineHistograms(fHistMan, "Event_BeforeCuts;Event_AfterCuts;", fConfigAddEventHistogram.value.data()); // define all histograms + DefineHistograms(fHistMan, "Event_BeforeCuts;Event_AfterCuts;SameBunchCorrelations", fConfigAddEventHistogram.value.data()); // define all histograms VarManager::SetUseVars(fHistMan->GetUsedVars()); // provide the list of required variables so that VarManager knows what to fill fOutputList.setObject(fHistMan->GetMainHistogramList()); } @@ -170,35 +178,99 @@ struct AnalysisEventSelection { } } - template - void runEventSelection(TEvent const& event) + template + void runEventSelection(TEvents const& events) { - // Reset the fValues array and fill event observables - VarManager::ResetValues(0, VarManager::kNEventWiseVariables); - VarManager::FillEvent(event); + fSelMap.clear(); + fBCCollMap.clear(); - // if QA is requested fill histograms before event selections - if (fConfigQA) { - fHistMan->FillHistClass("Event_BeforeCuts", VarManager::fgValues); // automatically fill all the histograms in the class Event - } - if (fEventCut->IsSelected(VarManager::fgValues)) { + for (auto& event : events) { + // Reset the fValues array and fill event observables + VarManager::ResetValues(0, VarManager::kNEventWiseVariables); + VarManager::FillEvent(event); + + bool decision = false; + // if QA is requested fill histograms before event selections if (fConfigQA) { - fHistMan->FillHistClass("Event_AfterCuts", VarManager::fgValues); + fHistMan->FillHistClass("Event_BeforeCuts", VarManager::fgValues); // automatically fill all the histograms in the class Event + } + if (fEventCut->IsSelected(VarManager::fgValues)) { + if (fConfigQA) { + fHistMan->FillHistClass("Event_AfterCuts", VarManager::fgValues); + } + decision = true; + } + fSelMap[event.globalIndex()] = decision; + if (fBCCollMap.find(event.globalBC()) == fBCCollMap.end()) { + std::vector evIndices = {event.globalIndex()}; + fBCCollMap[event.globalBC()] = evIndices; + } else { + auto& evIndices = fBCCollMap[event.globalBC()]; + evIndices.push_back(event.globalIndex()); + } + + if (fMixHandler != nullptr) { + int hh = fMixHandler->FindEventCategory(VarManager::fgValues); + hash(hh); } - eventSel(1); - } else { - eventSel(0); } + } + + template + void publishSelections(TEvents const& events) + { - if (fMixHandler != nullptr) { - int hh = fMixHandler->FindEventCategory(VarManager::fgValues); - hash(hh); + std::map collisionSplittingMap; // key: event global index, value: whether pileup event is a possible splitting + + // Reset the fValues array and fill event observables + VarManager::ResetValues(0, VarManager::kNEventWiseVariables); + // loop over the BC map, find BCs with more than one collision and compute 2-event correlation quantities + for (auto& [bc, evIndices] : fBCCollMap) { + if (evIndices.size() < 2) { + continue; + } + for (auto ev1Idx = evIndices.begin(); ev1Idx != evIndices.end(); ++ev1Idx) { + if (!fSelMap[*ev1Idx]) { + continue; + } + auto ev1 = events.rawIteratorAt(*ev1Idx); + for (auto ev2Idx = std::next(ev1Idx); ev2Idx != evIndices.end(); ++ev2Idx) { + if (!fSelMap[*ev2Idx]) { + continue; + } + auto ev2 = events.rawIteratorAt(*ev2Idx); + VarManager::FillTwoEvents(ev1, ev2); + if (TMath::Abs(VarManager::fgValues[VarManager::kTwoEvDeltaZ]) < 1.0) { // this is a possible collision split + collisionSplittingMap[*ev1Idx] = true; + collisionSplittingMap[*ev2Idx] = true; + } + fHistMan->FillHistClass("SameBunchCorrelations", VarManager::fgValues); + } + } + } + + // publish the table + uint32_t evSel = 0; + for (auto& event : events) { + evSel = 0; + if (fSelMap[event.globalIndex()]) { // event passed the user cuts + evSel |= (uint32_t(1) << 0); + } + std::vector sameBunchEvents = fBCCollMap[event.globalBC()]; + if (sameBunchEvents.size() > 1) { // event with in-bunch pileup + evSel |= (uint32_t(1) << 1); + } + if (collisionSplittingMap.find(event.globalIndex()) != collisionSplittingMap.end()) { // event with possible fake in-bunch pileup (collision splitting) + evSel |= (uint32_t(1) << 2); + } + eventSel(evSel); } } - void processSkimmed(MyEvents::iterator const& event) + void processSkimmed(MyEvents const& events) { - runEventSelection(event); + runEventSelection(events); + publishSelections(events); } void processDummy(MyEvents&) { @@ -213,6 +285,7 @@ struct AnalysisEventSelection { // Here one should add all the track cuts needed through the workflow (e.g. cuts for same-even pairing, electron prefiltering, track for dilepton-track correlations) struct AnalysisTrackSelection { Produces trackSel; + Produces trackAmbiguities; OutputObj fOutputList{"output"}; Configurable fConfigCuts{"cfgTrackCuts", "jpsiO2MCdebugCuts2", "Comma separated list of barrel track cuts"}; @@ -234,6 +307,9 @@ struct AnalysisTrackSelection { int fCurrentRun; // current run (needed to detect run changes for loading CCDB parameters) + std::map> fNAssocsInBunch; // key: track global index, value: vector of global index for events associated in-bunch (events that have in-bunch pileup or splitting) + std::map> fNAssocsOutOfBunch; // key: track global index, value: vector of global index for events associated out-of-bunch (events that have no in-bunch pileup) + void init(o2::framework::InitContext&) { fCurrentRun = 0; @@ -259,6 +335,7 @@ struct AnalysisTrackSelection { for (auto& cut : fTrackCuts) { histDirNames += Form("TrackBarrel_%s;", cut.GetName()); } + histDirNames += "TrackBarrel_AmbiguityInBunch;TrackBarrel_AmbiguityOutOfBunch;"; VarManager::SetRunlist((TString)fConfigRunPeriods); DefineHistograms(fHistMan, histDirNames.Data(), fConfigAddTrackHistogram.value.data()); // define all histograms @@ -281,6 +358,9 @@ struct AnalysisTrackSelection { template void runTrackSelection(ReducedTracksAssoc const& assocs, TEvents const& events, TTracks const& tracks) { + fNAssocsInBunch.clear(); + fNAssocsOutOfBunch.clear(); + if (fConfigComputeTPCpostCalib && events.size() > 0 && fCurrentRun != events.begin().runNumber()) { auto calibList = fCCDB->getForTimeStamp(fConfigCcdbPathTPC.value, events.begin().timestamp()); VarManager::SetCalibrationObject(VarManager::kTPCElectronMean, calibList->FindObject("mean_map_electron")); @@ -301,12 +381,13 @@ struct AnalysisTrackSelection { } trackSel.reserve(assocs.size()); + trackAmbiguities.reserve(tracks.size()); uint32_t filterMap = 0; int iCut = 0; for (auto& assoc : assocs) { auto event = assoc.template reducedevent_as(); - if (!event.isEventSelected()) { + if (!event.isEventSelected_bit(0)) { trackSel(0); continue; } @@ -333,7 +414,65 @@ struct AnalysisTrackSelection { } // end loop over cuts trackSel(filterMap); + + // count the number of associations per track + if (filterMap > 0) { + if (event.isEventSelected_bit(1)) { + if (fNAssocsInBunch.find(track.globalIndex()) == fNAssocsInBunch.end()) { + std::vector evVector = {event.globalIndex()}; + fNAssocsInBunch[track.globalIndex()] = evVector; + } else { + auto& evVector = fNAssocsInBunch[track.globalIndex()]; + evVector.push_back(event.globalIndex()); + } + } else { + if (fNAssocsOutOfBunch.find(track.globalIndex()) == fNAssocsOutOfBunch.end()) { + std::vector evVector = {event.globalIndex()}; + fNAssocsOutOfBunch[track.globalIndex()] = evVector; + } else { + auto& evVector = fNAssocsOutOfBunch[track.globalIndex()]; + evVector.push_back(event.globalIndex()); + } + } + } } // end loop over associations + + // QA the collision-track associations + for (auto& [trackIdx, evIndices] : fNAssocsInBunch) { + if (evIndices.size() == 1) { + continue; + } + auto track = tracks.rawIteratorAt(trackIdx); + VarManager::ResetValues(0, VarManager::kNBarrelTrackVariables); + VarManager::FillTrack(track); + VarManager::fgValues[VarManager::kBarrelNAssocsInBunch] = float(evIndices.size()); + fHistMan->FillHistClass("TrackBarrel_AmbiguityInBunch", VarManager::fgValues); + } // end loop over in-bunch ambiguous tracks + + for (auto& [trackIdx, evIndices] : fNAssocsOutOfBunch) { + if (evIndices.size() == 1) { + continue; + } + auto track = tracks.rawIteratorAt(trackIdx); + VarManager::ResetValues(0, VarManager::kNBarrelTrackVariables); + VarManager::FillTrack(track); + VarManager::fgValues[VarManager::kBarrelNAssocsOutOfBunch] = float(evIndices.size()); + fHistMan->FillHistClass("TrackBarrel_AmbiguityOutOfBunch", VarManager::fgValues); + } // end loop over out-of-bunch ambiguous tracks + + // publish the ambiguity table + for (auto& track : tracks) { + int8_t nInBunch = 0; + if (fNAssocsInBunch.find(track.globalIndex()) != fNAssocsInBunch.end()) { + nInBunch = fNAssocsInBunch[track.globalIndex()].size(); + } + int8_t nOutOfBunch = 0; + if (fNAssocsOutOfBunch.find(track.globalIndex()) != fNAssocsOutOfBunch.end()) { + nOutOfBunch = fNAssocsOutOfBunch[track.globalIndex()].size(); + } + trackAmbiguities(nInBunch, nOutOfBunch); + } + } // end runTrackSelection() void processSkimmed(ReducedTracksAssoc const& assocs, MyEventsSelected const& events, MyBarrelTracks const& tracks) @@ -428,7 +567,7 @@ struct AnalysisMuonSelection { for (auto& assoc : assocs) { auto event = assoc.template reducedevent_as(); - if (!event.isEventSelected()) { + if (!event.isEventSelected_bit(0)) { muonSel(0); continue; } @@ -663,7 +802,7 @@ struct AnalysisSameEventPairing { Service fCCDB; - Filter filterEventSelected = aod::dqanalysisflags::isEventSelected == 1; + Filter filterEventSelected = aod::dqanalysisflags::isEventSelected == uint32_t(1); HistogramManager* fHistMan; @@ -679,6 +818,11 @@ struct AnalysisSameEventPairing { int fNCutsMuon; int fNPairCuts; + bool fEnableBarrelMixingHistos; + bool fEnableBarrelHistos; + bool fEnableMuonHistos; + bool fEnableMuonMixingHistos; + NoBinningPolicy hashBin; Preslice> trackAssocsPerCollision = aod::reducedtrack_association::reducedeventId; @@ -692,13 +836,13 @@ struct AnalysisSameEventPairing { context.mOptions.get("processDecayToEEPrefilterSkimmed") || context.mOptions.get("processDecayToEEPrefilterSkimmedNoTwoProngFitter") || context.mOptions.get("processDecayToEESkimmedWithColl") || context.mOptions.get("processDecayToEESkimmedWithCollNoTwoProngFitter") || context.mOptions.get("processDecayToPiPiSkimmed") || context.mOptions.get("processAllSkimmed"); */ - bool enableBarrelHistos = context.mOptions.get("processAllSkimmed"); - bool enableBarrelMixingHistos = context.mOptions.get("processMixingAllSkimmed"); + fEnableBarrelHistos = context.mOptions.get("processAllSkimmed"); + fEnableBarrelMixingHistos = context.mOptions.get("processMixingAllSkimmed"); /*bool enableMuonHistos = context.mOptions.get("processDecayToMuMuSkimmed") || context.mOptions.get("processDecayToMuMuVertexingSkimmed") || context.mOptions.get("processDecayToMuMuSkimmedWithColl") || context.mOptions.get("processVnDecayToMuMuSkimmed") || context.mOptions.get("processAllSkimmed");*/ - bool enableMuonHistos = context.mOptions.get("processAllSkimmed"); - bool enableMuonMixingHistos = context.mOptions.get("processMixingAllSkimmed"); + fEnableMuonHistos = context.mOptions.get("processAllSkimmed"); + fEnableMuonMixingHistos = context.mOptions.get("processMixingAllSkimmed"); // Keep track of all the histogram class names to avoid composing strings in the pairing loop TString histNames = ""; @@ -737,18 +881,26 @@ struct AnalysisSameEventPairing { if (objArrayTrackCuts->FindObject(tempStr.Data()) != nullptr) { fTrackFilterMask |= (uint32_t(1) << icut); - if (enableBarrelHistos) { + if (fEnableBarrelHistos) { names = { Form("PairsBarrelSEPM_%s", objArray->At(icut)->GetName()), Form("PairsBarrelSEPP_%s", objArray->At(icut)->GetName()), Form("PairsBarrelSEMM_%s", objArray->At(icut)->GetName())}; histNames += Form("%s;%s;%s;", names[0].Data(), names[1].Data(), names[2].Data()); - if (enableBarrelMixingHistos) { + if (fEnableBarrelMixingHistos) { names.push_back(Form("PairsBarrelMEPM_%s", objArray->At(icut)->GetName())); names.push_back(Form("PairsBarrelMEPP_%s", objArray->At(icut)->GetName())); names.push_back(Form("PairsBarrelMEMM_%s", objArray->At(icut)->GetName())); histNames += Form("%s;%s;%s;", names[3].Data(), names[4].Data(), names[5].Data()); } + names.push_back(Form("PairsBarrelSEPM_ambiguousInBunch_%s", objArray->At(icut)->GetName())); + names.push_back(Form("PairsBarrelSEPP_ambiguousInBunch_%s", objArray->At(icut)->GetName())); + names.push_back(Form("PairsBarrelSEMM_ambiguousInBunch_%s", objArray->At(icut)->GetName())); + names.push_back(Form("PairsBarrelSEPM_ambiguousOutOfBunch_%s", objArray->At(icut)->GetName())); + names.push_back(Form("PairsBarrelSEPP_ambiguousOutOfBunch_%s", objArray->At(icut)->GetName())); + names.push_back(Form("PairsBarrelSEMM_ambiguousOutOfBunch_%s", objArray->At(icut)->GetName())); + histNames += Form("%s;%s;%s;", names[(fEnableBarrelMixingHistos ? 6 : 3)].Data(), names[(fEnableBarrelMixingHistos ? 7 : 4)].Data(), names[(fEnableBarrelMixingHistos ? 8 : 5)].Data()); + histNames += Form("%s;%s;%s;", names[(fEnableBarrelMixingHistos ? 9 : 6)].Data(), names[(fEnableBarrelMixingHistos ? 10 : 7)].Data(), names[(fEnableBarrelMixingHistos ? 11 : 8)].Data()); fTrackHistNames[icut] = names; TString cutNamesStr = fConfigPairCuts.value; @@ -779,14 +931,14 @@ struct AnalysisSameEventPairing { if (objArrayMuonCuts->FindObject(tempStr.Data()) != nullptr) { fMuonFilterMask |= (uint32_t(1) << icut); - if (enableMuonHistos) { + if (fEnableMuonHistos) { // no pair cuts names = { Form("PairsMuonSEPM_%s", objArray->At(icut)->GetName()), Form("PairsMuonSEPP_%s", objArray->At(icut)->GetName()), Form("PairsMuonSEMM_%s", objArray->At(icut)->GetName())}; histNames += Form("%s;%s;%s;", names[0].Data(), names[1].Data(), names[2].Data()); - if (enableMuonMixingHistos) { + if (fEnableMuonMixingHistos) { names.push_back(Form("PairsMuonMEPM_%s", objArray->At(icut)->GetName())); names.push_back(Form("PairsMuonMEPP_%s", objArray->At(icut)->GetName())); names.push_back(Form("PairsMuonMEMM_%s", objArray->At(icut)->GetName())); @@ -910,12 +1062,12 @@ struct AnalysisSameEventPairing { } // Template function to run same event pairing (barrel-barrel, muon-muon, barrel-muon) - template - void runSameEventPairing(TEvent const& event, TTrackAssocs const& assocs, TTracks const& tracks) + template + void runSameEventPairing(TEvents const& events, Preslice& preslice, TTrackAssocs const& assocs, TTracks const& tracks) { - if (fCurrentRun != event.runNumber()) { - initParamsFromCCDB(event.timestamp(), TTwoProngFitter); - fCurrentRun = event.runNumber(); + if (fCurrentRun != events.begin().runNumber()) { + initParamsFromCCDB(events.begin().timestamp(), TTwoProngFitter); + fCurrentRun = events.begin().runNumber(); } TString cutNames = fConfigTrackCuts.value; @@ -926,6 +1078,12 @@ struct AnalysisSameEventPairing { histNames = fMuonHistNames; ncuts = fNCutsMuon; } + int histIdxOffset = 0; + if constexpr (TPairType == pairTypeEE) { + if (fEnableBarrelMixingHistos) { + histIdxOffset = 3; + } + } /*if constexpr (TPairType == pairTypeEMu) { cutNames = fConfigMuonCuts.value; histNames = fTrackMuonHistNames; @@ -945,136 +1103,184 @@ struct AnalysisSameEventPairing { constexpr bool eventHasQvector = ((TEventFillMap & VarManager::ObjTypes::ReducedEventQvector) > 0); constexpr bool trackHasCov = ((TTrackFillMap & VarManager::ObjTypes::TrackCov) > 0 || (TTrackFillMap & VarManager::ObjTypes::ReducedTrackBarrelCov) > 0); - for (auto& [a1, a2] : o2::soa::combinations(assocs, assocs)) { + for (auto& event : events) { + if (!event.isEventSelected_bit(0)) { + continue; + } + // Reset the fValues array + VarManager::ResetValues(0, VarManager::kNVars); + VarManager::FillEvent(event, VarManager::fgValues); - if constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi) { - twoTrackFilter = a1.isBarrelSelected_raw() & a2.isBarrelSelected_raw() & a1.isBarrelSelectedPrefilter_raw() & a2.isBarrelSelectedPrefilter_raw() & fTrackFilterMask; + auto groupedAssocs = assocs.sliceBy(preslice, event.globalIndex()); + if (groupedAssocs.size() == 0) { + continue; + } - if (!twoTrackFilter) { // the tracks must have at least one filter bit in common to continue - continue; - } + for (auto& [a1, a2] : o2::soa::combinations(groupedAssocs, groupedAssocs)) { - auto t1 = a1.template reducedtrack_as(); - auto t2 = a2.template reducedtrack_as(); - sign1 = t1.sign(); - sign2 = t2.sign(); + if constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi) { + twoTrackFilter = a1.isBarrelSelected_raw() & a2.isBarrelSelected_raw() & a1.isBarrelSelectedPrefilter_raw() & a2.isBarrelSelectedPrefilter_raw() & fTrackFilterMask; - VarManager::FillPair(t1, t2); - if constexpr (TTwoProngFitter) { - VarManager::FillPairVertexing(event, t1, t2, fConfigPropToPCA); - } - if constexpr (eventHasQvector) { - VarManager::FillPairVn(t1, t2); - } + if (!twoTrackFilter) { // the tracks must have at least one filter bit in common to continue + continue; + } - dielectronList(event, VarManager::fgValues[VarManager::kMass], - VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], - t1.sign() + t2.sign(), twoTrackFilter, 0); + auto t1 = a1.template reducedtrack_as(); + auto t2 = a2.template reducedtrack_as(); + sign1 = t1.sign(); + sign2 = t2.sign(); + // store the ambiguity number of the two dilepton legs in the last 4 digits of the two-track filter + if (t1.barrelAmbiguityInBunch() > 1) { + twoTrackFilter |= (uint32_t(1) << 28); + } + if (t2.barrelAmbiguityInBunch() > 1) { + twoTrackFilter |= (uint32_t(1) << 29); + } + if (t1.barrelAmbiguityOutOfBunch() > 1) { + twoTrackFilter |= (uint32_t(1) << 30); + } + if (t2.barrelAmbiguityOutOfBunch() > 1) { + twoTrackFilter |= (uint32_t(1) << 31); + } - if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrackCollInfo) > 0) { - dileptonInfoList(t1.collisionId(), event.posX(), event.posY(), event.posZ()); - } - if constexpr (trackHasCov && TTwoProngFitter) { - dielectronsExtraList(t1.globalIndex(), t2.globalIndex(), VarManager::fgValues[VarManager::kVertexingTauz], VarManager::fgValues[VarManager::kVertexingLz], VarManager::fgValues[VarManager::kVertexingLxy]); - } - } + VarManager::FillPair(t1, t2); + if constexpr (TTwoProngFitter) { + VarManager::FillPairVertexing(event, t1, t2, fConfigPropToPCA); + } + if constexpr (eventHasQvector) { + VarManager::FillPairVn(t1, t2); + } - if constexpr (TPairType == VarManager::kDecayToMuMu) { + dielectronList(event.globalIndex(), VarManager::fgValues[VarManager::kMass], + VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], + t1.sign() + t2.sign(), twoTrackFilter, 0); - twoTrackFilter = a1.isMuonSelected_raw() & a2.isMuonSelected_raw() & fMuonFilterMask; - if (!twoTrackFilter) { // the tracks must have at least one filter bit in common to continue - continue; + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrackCollInfo) > 0) { + dileptonInfoList(t1.collisionId(), event.posX(), event.posY(), event.posZ()); + } + if constexpr (trackHasCov && TTwoProngFitter) { + dielectronsExtraList(t1.globalIndex(), t2.globalIndex(), VarManager::fgValues[VarManager::kVertexingTauz], VarManager::fgValues[VarManager::kVertexingLz], VarManager::fgValues[VarManager::kVertexingLxy]); + } } - auto t1 = a1.template reducedmuon_as(); - auto t2 = a2.template reducedmuon_as(); - sign1 = t1.sign(); - sign2 = t2.sign(); - - VarManager::FillPair(t1, t2); - if constexpr (TTwoProngFitter) { - VarManager::FillPairVertexing(event, t1, t2, fConfigPropToPCA); - } - if constexpr (eventHasQvector) { - VarManager::FillPairVn(t1, t2); - } + if constexpr (TPairType == VarManager::kDecayToMuMu) { + twoTrackFilter = a1.isMuonSelected_raw() & a2.isMuonSelected_raw() & fMuonFilterMask; + if (!twoTrackFilter) { // the tracks must have at least one filter bit in common to continue + continue; + } - dimuonList(event, VarManager::fgValues[VarManager::kMass], - VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], - t1.sign() + t2.sign(), twoTrackFilter, 0); - if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedMuonCollInfo) > 0) { - dileptonInfoList(t1.collisionId(), event.posX(), event.posY(), event.posZ()); - } + auto t1 = a1.template reducedmuon_as(); + auto t2 = a2.template reducedmuon_as(); + sign1 = t1.sign(); + sign2 = t2.sign(); - if constexpr (TTwoProngFitter) { - dimuonsExtraList(t1.globalIndex(), t2.globalIndex(), VarManager::fgValues[VarManager::kVertexingTauz], VarManager::fgValues[VarManager::kVertexingLz], VarManager::fgValues[VarManager::kVertexingLxy]); - if (fConfigFlatTables.value) { - dimuonAllList(event.posX(), event.posY(), event.posZ(), event.numContrib(), - -999., -999., -999., - VarManager::fgValues[VarManager::kMass], - false, - VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], t1.sign() + t2.sign(), VarManager::fgValues[VarManager::kVertexingChi2PCA], - VarManager::fgValues[VarManager::kVertexingTauz], VarManager::fgValues[VarManager::kVertexingTauzErr], - VarManager::fgValues[VarManager::kVertexingTauxy], VarManager::fgValues[VarManager::kVertexingTauxyErr], - VarManager::fgValues[VarManager::kCosPointingAngle], - VarManager::fgValues[VarManager::kPt1], VarManager::fgValues[VarManager::kEta1], VarManager::fgValues[VarManager::kPhi1], t1.sign(), - VarManager::fgValues[VarManager::kPt2], VarManager::fgValues[VarManager::kEta2], VarManager::fgValues[VarManager::kPhi2], t2.sign(), - t1.fwdDcaX(), t1.fwdDcaY(), t2.fwdDcaX(), t2.fwdDcaY(), - 0., 0., - t1.chi2MatchMCHMID(), t2.chi2MatchMCHMID(), - t1.chi2MatchMCHMFT(), t2.chi2MatchMCHMFT(), - t1.chi2(), t2.chi2(), - -999., -999., -999., -999., - -999., -999., -999., -999., - -999., -999., -999., -999., - -999., -999., -999., -999., - t1.isAmbiguous(), t2.isAmbiguous(), - VarManager::fgValues[VarManager::kU2Q2], VarManager::fgValues[VarManager::kU3Q3], - VarManager::fgValues[VarManager::kR2EP], VarManager::fgValues[VarManager::kR2SP], VarManager::fgValues[VarManager::kCentFT0C], - VarManager::fgValues[VarManager::kCos2DeltaPhi], VarManager::fgValues[VarManager::kCos3DeltaPhi], - VarManager::fgValues[VarManager::kVertexingPz], - VarManager::fgValues[VarManager::kVertexingSV]); + VarManager::FillPair(t1, t2); + if constexpr (TTwoProngFitter) { + VarManager::FillPairVertexing(event, t1, t2, fConfigPropToPCA); + } + if constexpr (eventHasQvector) { + VarManager::FillPairVn(t1, t2); } - } - } - // TODO: the model for the electron-muon combination has to be thought through - /*if constexpr (TPairType == VarManager::kElectronMuon) { - twoTrackFilter = a1.isBarrelSelected_raw() & a1.isBarrelSelectedPrefilter_raw() & a2.isMuonSelected_raw() & fTwoTrackFilterMask; - }*/ - if constexpr (eventHasQvector) { - dileptonFlowList(VarManager::fgValues[VarManager::kU2Q2], VarManager::fgValues[VarManager::kU3Q3], VarManager::fgValues[VarManager::kCos2DeltaPhi], VarManager::fgValues[VarManager::kCos3DeltaPhi]); - } + dimuonList(event.globalIndex(), VarManager::fgValues[VarManager::kMass], + VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], + t1.sign() + t2.sign(), twoTrackFilter, 0); + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedMuonCollInfo) > 0) { + dileptonInfoList(t1.collisionId(), event.posX(), event.posY(), event.posZ()); + } - // Fill histograms - for (int icut = 0; icut < ncuts; icut++) { - if (twoTrackFilter & (uint32_t(1) << icut)) { - if (sign1 * sign2 < 0) { - fHistMan->FillHistClass(histNames[icut][0].Data(), VarManager::fgValues); - } else { - if (sign1 > 0) { - fHistMan->FillHistClass(histNames[icut][1].Data(), VarManager::fgValues); - } else { - fHistMan->FillHistClass(histNames[icut][2].Data(), VarManager::fgValues); + if constexpr (TTwoProngFitter) { + dimuonsExtraList(t1.globalIndex(), t2.globalIndex(), VarManager::fgValues[VarManager::kVertexingTauz], VarManager::fgValues[VarManager::kVertexingLz], VarManager::fgValues[VarManager::kVertexingLxy]); + if (fConfigFlatTables.value) { + dimuonAllList(event.posX(), event.posY(), event.posZ(), event.numContrib(), + -999., -999., -999., + VarManager::fgValues[VarManager::kMass], + false, + VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], t1.sign() + t2.sign(), VarManager::fgValues[VarManager::kVertexingChi2PCA], + VarManager::fgValues[VarManager::kVertexingTauz], VarManager::fgValues[VarManager::kVertexingTauzErr], + VarManager::fgValues[VarManager::kVertexingTauxy], VarManager::fgValues[VarManager::kVertexingTauxyErr], + VarManager::fgValues[VarManager::kCosPointingAngle], + VarManager::fgValues[VarManager::kPt1], VarManager::fgValues[VarManager::kEta1], VarManager::fgValues[VarManager::kPhi1], t1.sign(), + VarManager::fgValues[VarManager::kPt2], VarManager::fgValues[VarManager::kEta2], VarManager::fgValues[VarManager::kPhi2], t2.sign(), + t1.fwdDcaX(), t1.fwdDcaY(), t2.fwdDcaX(), t2.fwdDcaY(), + 0., 0., + t1.chi2MatchMCHMID(), t2.chi2MatchMCHMID(), + t1.chi2MatchMCHMFT(), t2.chi2MatchMCHMFT(), + t1.chi2(), t2.chi2(), + -999., -999., -999., -999., + -999., -999., -999., -999., + -999., -999., -999., -999., + -999., -999., -999., -999., + t1.isAmbiguous(), t2.isAmbiguous(), + VarManager::fgValues[VarManager::kU2Q2], VarManager::fgValues[VarManager::kU3Q3], + VarManager::fgValues[VarManager::kR2EP], VarManager::fgValues[VarManager::kR2SP], VarManager::fgValues[VarManager::kCentFT0C], + VarManager::fgValues[VarManager::kCos2DeltaPhi], VarManager::fgValues[VarManager::kCos3DeltaPhi], + VarManager::fgValues[VarManager::kVertexingPz], + VarManager::fgValues[VarManager::kVertexingSV]); } } - for (unsigned int iPairCut = 0; iPairCut < fPairCuts.size(); iPairCut++) { - AnalysisCompositeCut cut = fPairCuts.at(iPairCut); - if (!(cut.IsSelected(VarManager::fgValues))) // apply pair cuts - continue; + } + // TODO: the model for the electron-muon combination has to be thought through + /*if constexpr (TPairType == VarManager::kElectronMuon) { + twoTrackFilter = a1.isBarrelSelected_raw() & a1.isBarrelSelectedPrefilter_raw() & a2.isMuonSelected_raw() & fTwoTrackFilterMask; + }*/ + + if constexpr (eventHasQvector) { + dileptonFlowList(VarManager::fgValues[VarManager::kU2Q2], VarManager::fgValues[VarManager::kU3Q3], VarManager::fgValues[VarManager::kCos2DeltaPhi], VarManager::fgValues[VarManager::kCos3DeltaPhi]); + } + + // Fill histograms + bool isAmbiInBunch = false; + bool isAmbiOutOfBunch = false; + for (int icut = 0; icut < ncuts; icut++) { + if (twoTrackFilter & (uint32_t(1) << icut)) { + isAmbiInBunch = (twoTrackFilter & (uint32_t(1) << 28)) || (twoTrackFilter & (uint32_t(1) << 29)); + isAmbiOutOfBunch = (twoTrackFilter & (uint32_t(1) << 30)) || (twoTrackFilter & (uint32_t(1) << 31)); if (sign1 * sign2 < 0) { - fHistMan->FillHistClass(histNames[ncuts + icut * ncuts + iPairCut][0].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNames[icut][0].Data(), VarManager::fgValues); + if (isAmbiInBunch) { + fHistMan->FillHistClass(histNames[icut][3 + histIdxOffset].Data(), VarManager::fgValues); + } + if (isAmbiOutOfBunch) { + fHistMan->FillHistClass(histNames[icut][3 + histIdxOffset + 3].Data(), VarManager::fgValues); + } } else { if (sign1 > 0) { - fHistMan->FillHistClass(histNames[ncuts + icut * ncuts + iPairCut][1].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNames[icut][1].Data(), VarManager::fgValues); + if (isAmbiInBunch) { + fHistMan->FillHistClass(histNames[icut][4 + histIdxOffset].Data(), VarManager::fgValues); + } + if (isAmbiOutOfBunch) { + fHistMan->FillHistClass(histNames[icut][4 + histIdxOffset + 3].Data(), VarManager::fgValues); + } } else { - fHistMan->FillHistClass(histNames[ncuts + icut * ncuts + iPairCut][2].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNames[icut][2].Data(), VarManager::fgValues); + if (isAmbiInBunch) { + fHistMan->FillHistClass(histNames[icut][5 + histIdxOffset].Data(), VarManager::fgValues); + } + if (isAmbiOutOfBunch) { + fHistMan->FillHistClass(histNames[icut][5 + histIdxOffset + 3].Data(), VarManager::fgValues); + } } } - } // end loop (pair cuts) - } - } // end loop (cuts) - } // end loop over pairs of track associations + for (unsigned int iPairCut = 0; iPairCut < fPairCuts.size(); iPairCut++) { + AnalysisCompositeCut cut = fPairCuts.at(iPairCut); + if (!(cut.IsSelected(VarManager::fgValues))) // apply pair cuts + continue; + if (sign1 * sign2 < 0) { + fHistMan->FillHistClass(histNames[ncuts + icut * ncuts + iPairCut][0].Data(), VarManager::fgValues); + } else { + if (sign1 > 0) { + fHistMan->FillHistClass(histNames[ncuts + icut * ncuts + iPairCut][1].Data(), VarManager::fgValues); + } else { + fHistMan->FillHistClass(histNames[ncuts + icut * ncuts + iPairCut][2].Data(), VarManager::fgValues); + } + } + } // end loop (pair cuts) + } + } // end loop (cuts) + } // end loop over pairs of track associations + } // end loop over events } template @@ -1118,7 +1324,7 @@ struct AnalysisSameEventPairing { twoTrackFilter = a1.isBarrelSelected_raw() & a1.isBarrelSelectedPrefilter_raw() & a2.isMuonSelected_raw() & fTrackFilterMask; }*/ - for (unsigned int icut = 0; icut < ncuts; icut++) { + for (int icut = 0; icut < ncuts; icut++) { if (!(twoTrackFilter & (uint32_t(1) << icut))) { continue; // cut not passed } @@ -1157,28 +1363,12 @@ struct AnalysisSameEventPairing { } void processAllSkimmed(MyEventsVtxCovSelected const& events, - soa::Join const& barrelAssocs, MyBarrelTracksWithCov const& barrelTracks, + soa::Join const& barrelAssocs, MyBarrelTracksWithCovWithAmbiguities const& barrelTracks, soa::Join const& muonAssocs, MyMuonTracksWithCov const& muons) { - for (auto& event : events) { - if (event.isEventSelected() == 0) { - continue; - } - // Reset the fValues array - VarManager::ResetValues(0, VarManager::kNVars); - VarManager::FillEvent(event, VarManager::fgValues); - - auto groupedBarrelAssocs = barrelAssocs.sliceBy(trackAssocsPerCollision, event.globalIndex()); - auto groupedMuonAssocs = muonAssocs.sliceBy(muonAssocsPerCollision, event.globalIndex()); - - if (groupedBarrelAssocs.size() > 0) { - runSameEventPairing(event, groupedBarrelAssocs, barrelTracks); - } - if (groupedMuonAssocs.size() > 0) { - runSameEventPairing(event, groupedMuonAssocs, muons); - } - // runSameEventPairing(event, tracks, muons); - } + runSameEventPairing(events, trackAssocsPerCollision, barrelAssocs, barrelTracks); + runSameEventPairing(events, muonAssocsPerCollision, muonAssocs, muons); + // runSameEventPairing(event, tracks, muons); } void processMixingAllSkimmed(soa::Filtered& events, @@ -1373,8 +1563,8 @@ struct AnalysisDileptonTrack { Service fCCDB; // TODO: The filter expressions seem to always use the default value of configurables, not the values from the actual configuration file - Filter eventFilter = aod::dqanalysisflags::isEventSelected == 1; - Filter dileptonFilter = aod::reducedpair::pt > fConfigDileptonpTCut.value&& aod::reducedpair::mass > fConfigDileptonLowMass.value&& aod::reducedpair::mass < fConfigDileptonHighMass.value&& aod::reducedpair::sign == 0; + Filter eventFilter = aod::dqanalysisflags::isEventSelected == uint32_t(1); + Filter dileptonFilter = aod::reducedpair::pt > fConfigDileptonpTCut&& aod::reducedpair::mass > fConfigDileptonLowMass&& aod::reducedpair::mass < fConfigDileptonHighMass&& aod::reducedpair::sign == 0; Filter filterBarrel = aod::dqanalysisflags::isBarrelSelected > uint32_t(0); Filter filterMuon = aod::dqanalysisflags::isMuonSelected > uint32_t(0); @@ -1711,6 +1901,10 @@ void DefineHistograms(HistogramManager* histMan, TString histClasses, const char dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "event", histName); } + if (classStr.Contains("SameBunchCorrelations")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "two-collisions", histName); + } + if (classStr.Contains("Track") && !classStr.Contains("Pairs")) { if (classStr.Contains("Barrel")) { dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "track", histName); @@ -1723,6 +1917,9 @@ void DefineHistograms(HistogramManager* histMan, TString histClasses, const char if (classStr.Contains("PIDCalibProton")) { dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "track", "postcalib_proton"); } + if (classStr.Contains("Ambiguity")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "track", Form("%s,ambiguity", histName.Data())); + } } if (classStr.Contains("Muon")) { dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "track", histName); From 4ae8639b1d3020c6c2ae1a97abcd342df20bb377 Mon Sep 17 00:00:00 2001 From: ldellost <47105254+DelloStritto@users.noreply.github.com> Date: Mon, 19 Feb 2024 22:02:20 +0100 Subject: [PATCH 5/8] PWGHF: New features in LcpKpi task and tree creator (#4769) Co-authored-by: Luigi Dello Stritto --- PWGHF/TableProducer/candidateSelectorLc.cxx | 20 +++++----- PWGHF/TableProducer/treeCreatorLcToPKPi.cxx | 42 ++++++++++++--------- 2 files changed, 35 insertions(+), 27 deletions(-) diff --git a/PWGHF/TableProducer/candidateSelectorLc.cxx b/PWGHF/TableProducer/candidateSelectorLc.cxx index ecb7156e6cf..e09534e5e35 100644 --- a/PWGHF/TableProducer/candidateSelectorLc.cxx +++ b/PWGHF/TableProducer/candidateSelectorLc.cxx @@ -365,17 +365,12 @@ struct HfCandidateSelectorLc { continue; } - if ((pidLcToPKPi == -1 || pidLcToPKPi == 1) && (pidBayesLcToPKPi == -1 || pidBayesLcToPKPi == 1) && topolLcToPKPi) { - statusLcToPKPi = 1; // identified as LcToPKPi - } - if ((pidLcToPiKP == -1 || pidLcToPiKP == 1) && (pidBayesLcToPiKP == -1 || pidBayesLcToPiKP == 1) && topolLcToPiKP) { - statusLcToPiKP = 1; // identified as LcToPiKP - } - + bool isSelectedMlLcToPKPi = true; + bool isSelectedMlLcToPiKP = true; if (applyMl) { // ML selections - bool isSelectedMlLcToPKPi = false; - bool isSelectedMlLcToPiKP = false; + isSelectedMlLcToPKPi = false; + isSelectedMlLcToPiKP = false; if ((pidLcToPKPi == -1 || pidLcToPKPi == 1) && (pidBayesLcToPKPi == -1 || pidBayesLcToPKPi == 1) && topolLcToPKPi) { std::vector inputFeaturesLcToPKPi = hfMlResponse.getInputFeatures(candidate, trackPos1, trackNeg, trackPos2); @@ -398,6 +393,13 @@ struct HfCandidateSelectorLc { } } + if ((pidLcToPKPi == -1 || pidLcToPKPi == 1) && (pidBayesLcToPKPi == -1 || pidBayesLcToPKPi == 1) && isSelectedMlLcToPKPi && topolLcToPKPi) { + statusLcToPKPi = 1; // identified as LcToPKPi + } + if ((pidLcToPiKP == -1 || pidLcToPiKP == 1) && (pidBayesLcToPiKP == -1 || pidBayesLcToPiKP == 1) && isSelectedMlLcToPiKP && topolLcToPiKP) { + statusLcToPiKP = 1; // identified as LcToPiKP + } + hfSelLcCandidate(statusLcToPKPi, statusLcToPiKP); } } diff --git a/PWGHF/TableProducer/treeCreatorLcToPKPi.cxx b/PWGHF/TableProducer/treeCreatorLcToPKPi.cxx index 53b3f8b7380..d0b20e1422d 100644 --- a/PWGHF/TableProducer/treeCreatorLcToPKPi.cxx +++ b/PWGHF/TableProducer/treeCreatorLcToPKPi.cxx @@ -25,6 +25,7 @@ #include "PWGHF/DataModel/CandidateReconstructionTables.h" #include "PWGHF/DataModel/CandidateSelectionTables.h" #include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/Centrality.h" using namespace o2; using namespace o2::framework; @@ -84,10 +85,11 @@ DECLARE_SOA_INDEX_COLUMN(McParticle, mcParticle); DECLARE_SOA_INDEX_COLUMN(McCollision, mcCollision); DECLARE_SOA_COLUMN(IsEventReject, isEventReject, int); DECLARE_SOA_COLUMN(RunNumber, runNumber, int); -DECLARE_SOA_COLUMN(MultZeqFT0A, multZeqFT0A, float); -DECLARE_SOA_COLUMN(MultZeqFT0C, multZeqFT0C, float); -DECLARE_SOA_COLUMN(MultFT0M, multFT0M, float); -DECLARE_SOA_COLUMN(MultZeqFV0A, multZeqFV0A, float); +DECLARE_SOA_COLUMN(CentFT0A, centFT0A, float); +DECLARE_SOA_COLUMN(CentFT0C, centFT0C, float); +DECLARE_SOA_COLUMN(CentFT0M, centFT0M, float); +DECLARE_SOA_COLUMN(CentFV0A, centFV0A, float); +DECLARE_SOA_COLUMN(CentFDDM, centFDDM, float); DECLARE_SOA_COLUMN(MultZeqNTracksPV, multZeqNTracksPV, float); } // namespace full @@ -227,10 +229,11 @@ DECLARE_SOA_TABLE(HfCandLcFullEvs, "AOD", "HFCANDLCFULLEV", collision::PosZ, full::IsEventReject, full::RunNumber, - full::MultZeqFT0A, - full::MultZeqFT0C, - full::MultFT0M, - full::MultZeqFV0A, + full::CentFT0A, + full::CentFT0C, + full::CentFT0M, + full::CentFV0A, + full::CentFDDM, full::MultZeqNTracksPV); DECLARE_SOA_TABLE(HfCandLcFullPs, "AOD", "HFCANDLCFULLP", @@ -263,12 +266,13 @@ struct HfTreeCreatorLcToPKPi { HfHelper hfHelper; using TracksWPid = soa::Join; + using Cents = soa::Join; void init(InitContext const&) { } - void processMc(soa::Join const& collisions, + void processMc(soa::Join const& collisions, aod::McCollisions const& mcCollisions, soa::Join const& candidates, soa::Join const& particles, @@ -287,10 +291,11 @@ struct HfTreeCreatorLcToPKPi { collision.posZ(), 0, collision.bc().runNumber(), - collision.multZeqFT0A(), - collision.multZeqFT0C(), - collision.multFT0M(), - collision.multZeqFV0A(), + collision.centFT0A(), + collision.centFT0C(), + collision.centFT0M(), + collision.centFV0A(), + collision.centFDDM(), collision.multZeqNTracksPV()); } @@ -474,7 +479,7 @@ struct HfTreeCreatorLcToPKPi { } PROCESS_SWITCH(HfTreeCreatorLcToPKPi, processMc, "Process MC tree writer", true); - void processData(soa::Join const& collisions, + void processData(soa::Join const& collisions, soa::Join const& candidates, TracksWPid const& tracks, aod::BCs const&) { @@ -491,10 +496,11 @@ struct HfTreeCreatorLcToPKPi { collision.posZ(), 0, collision.bc().runNumber(), - collision.multZeqFT0A(), - collision.multZeqFT0C(), - collision.multFT0M(), - collision.multZeqFV0A(), + collision.centFT0A(), + collision.centFT0C(), + collision.centFT0M(), + collision.centFV0A(), + collision.centFDDM(), collision.multZeqNTracksPV()); } From 74a1282197cc6e339fed0579b5ef3fd133ef5333 Mon Sep 17 00:00:00 2001 From: rbailhac Date: Mon, 19 Feb 2024 23:03:21 +0100 Subject: [PATCH 6/8] PbPb TOF required QA (#4771) --- PWGDQ/Core/CutsLibrary.cxx | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/PWGDQ/Core/CutsLibrary.cxx b/PWGDQ/Core/CutsLibrary.cxx index be99391fbba..c1c872ed803 100644 --- a/PWGDQ/Core/CutsLibrary.cxx +++ b/PWGDQ/Core/CutsLibrary.cxx @@ -1360,6 +1360,43 @@ AnalysisCompositeCut* o2::aod::dqcuts::GetCompositeCut(const char* cutName) return cut; } + // 4 cuts for QC + if (!nameStr.compare(Form("lmee_posTOFreqRun3_posEta%s_strongNSigEPbPb_rejBadTOF", vecPIDcase.at(icase).Data()))) { + cut->AddCut(GetAnalysisCut("posTrack")); + cut->AddCut(GetAnalysisCut("posEtaSel")); + cut->AddCut(GetAnalysisCut("TightGlobalTrackRun3")); + cut->AddCut(GetAnalysisCut("standardPrimaryTrackDCAz")); + cut->AddCut(GetAnalysisCut(Form("electronPID_TOFreq%s_strongNSigEPbPb_rejBadTOF", vecPIDcase.at(icase).Data()))); + return cut; + } + + if (!nameStr.compare(Form("lmee_posTOFreqRun3_negEta%s_strongNSigEPbPb_rejBadTOF", vecPIDcase.at(icase).Data()))) { + cut->AddCut(GetAnalysisCut("posTrack")); + cut->AddCut(GetAnalysisCut("negEtaSel")); + cut->AddCut(GetAnalysisCut("TightGlobalTrackRun3")); + cut->AddCut(GetAnalysisCut("standardPrimaryTrackDCAz")); + cut->AddCut(GetAnalysisCut(Form("electronPID_TOFreq%s_strongNSigEPbPb_rejBadTOF", vecPIDcase.at(icase).Data()))); + return cut; + } + + if (!nameStr.compare(Form("lmee_negTOFreqRun3_posEta%s_strongNSigEPbPb_rejBadTOF", vecPIDcase.at(icase).Data()))) { + cut->AddCut(GetAnalysisCut("negTrack")); + cut->AddCut(GetAnalysisCut("posEtaSel")); + cut->AddCut(GetAnalysisCut("TightGlobalTrackRun3")); + cut->AddCut(GetAnalysisCut("standardPrimaryTrackDCAz")); + cut->AddCut(GetAnalysisCut(Form("electronPID_TOFreq%s_strongNSigEPbPb_rejBadTOF", vecPIDcase.at(icase).Data()))); + return cut; + } + + if (!nameStr.compare(Form("lmee_negTOFreqRun3_negEta%s_strongNSigEPbPb_rejBadTOF", vecPIDcase.at(icase).Data()))) { + cut->AddCut(GetAnalysisCut("negTrack")); + cut->AddCut(GetAnalysisCut("negEtaSel")); + cut->AddCut(GetAnalysisCut("TightGlobalTrackRun3")); + cut->AddCut(GetAnalysisCut("standardPrimaryTrackDCAz")); + cut->AddCut(GetAnalysisCut(Form("electronPID_TOFreq%s_strongNSigEPbPb_rejBadTOF", vecPIDcase.at(icase).Data()))); + return cut; + } + if (!nameStr.compare(Form("lmee_eNSigmaRun3%s_TPC_PID", vecPIDcase.at(icase).Data()))) { cut->AddCut(GetAnalysisCut("lmeeStandardKine")); cut->AddCut(GetAnalysisCut("TightGlobalTrackRun3")); From 7c51975cad466aadf74bc25e4a54e7b8c29b2423 Mon Sep 17 00:00:00 2001 From: Nima Zardoshti Date: Mon, 19 Feb 2024 15:43:01 +0100 Subject: [PATCH 7/8] PWGJE: fixing bug in Substructure Output matching --- PWGJE/DataModel/JetReducedData.h | 2 +- PWGJE/DataModel/JetSubstructure.h | 86 ++++----- .../TableProducer/jetderiveddataproducer.cxx | 2 +- PWGJE/TableProducer/jetderiveddatawriter.cxx | 15 +- PWGJE/Tasks/jetsubstructurehfoutput.cxx | 178 ++++++++++-------- PWGJE/Tasks/jetsubstructureoutput.cxx | 136 +++++++------ 6 files changed, 231 insertions(+), 188 deletions(-) diff --git a/PWGJE/DataModel/JetReducedData.h b/PWGJE/DataModel/JetReducedData.h index 073dd9fba93..a17accca23e 100644 --- a/PWGJE/DataModel/JetReducedData.h +++ b/PWGJE/DataModel/JetReducedData.h @@ -175,7 +175,7 @@ DECLARE_SOA_COLUMN(Pt, pt, float); DECLARE_SOA_COLUMN(Eta, eta, float); DECLARE_SOA_COLUMN(Phi, phi, float); DECLARE_SOA_COLUMN(Energy, energy, float); -DECLARE_SOA_COLUMN(Sign, sign, float); +DECLARE_SOA_COLUMN(Sign, sign, int8_t); DECLARE_SOA_COLUMN(TrackSel, trackSel, uint8_t); DECLARE_SOA_DYNAMIC_COLUMN(Px, px, [](float pt, float phi) -> float { return pt * std::cos(phi); }); diff --git a/PWGJE/DataModel/JetSubstructure.h b/PWGJE/DataModel/JetSubstructure.h index 49be4461c6c..31a8872bd50 100644 --- a/PWGJE/DataModel/JetSubstructure.h +++ b/PWGJE/DataModel/JetSubstructure.h @@ -55,48 +55,50 @@ DECLARE_SOA_COLUMN(JetNConstituents, jetNConstituents, int); //! } // namespace jetoutput // Defines the jet substrcuture table definition -#define JETSUBSTRUCTURE_TABLE_DEF(_collision_type_, _table_type_, _jet_type_, _matched_jet_type_, _cand_type_, _name_, _description_) \ - \ - namespace _name_##collisionoutput \ - { \ - DECLARE_SOA_DYNAMIC_COLUMN(Dummy##_collision_type_, dummy##_collision_type_, []() -> int { return 0; }); \ - } \ - \ - DECLARE_SOA_TABLE(_table_type_##COs, "AOD", _description_ "CO", jetcollision::PosZ, jetcollision::Centrality, jetcollision::EventSel, _name_##collisionoutput::Dummy##_collision_type_<>); \ - using _table_type_##CO = _table_type_##COs::iterator; \ - \ - namespace _name_##geomatched \ - { \ - DECLARE_SOA_ARRAY_INDEX_COLUMN_FULL(_matched_jet_type_, matchedJetGeo, int32_t, _matched_jet_type_##s, "_geo"); \ - } \ - \ - namespace _name_##ptmatched \ - { \ - DECLARE_SOA_ARRAY_INDEX_COLUMN_FULL(_matched_jet_type_, matchedJetPt, int32_t, _matched_jet_type_##s, "_pt"); \ - } \ - \ - namespace _name_##candmatched \ - { \ - DECLARE_SOA_ARRAY_INDEX_COLUMN_FULL(_matched_jet_type_, matchedJetCand, int32_t, _matched_jet_type_##s, "_hf"); \ - } \ - namespace _name_##jetoutput \ - { \ - DECLARE_SOA_INDEX_COLUMN(_table_type_##CO, collision); \ - DECLARE_SOA_INDEX_COLUMN_FULL(Candidate, candidate, int, _cand_type_, "_0"); \ - } \ - DECLARE_SOA_TABLE(_table_type_##Os, "AOD", _description_ "O", _name_##jetoutput::_table_type_##COId, _name_##jetoutput::CandidateId, _name_##geomatched::_matched_jet_type_##Ids, _name_##ptmatched::_matched_jet_type_##Ids, _name_##candmatched::_matched_jet_type_##Ids, jetoutput::JetPt, jetoutput::JetPhi, jetoutput::JetEta, jetoutput::JetR, jetoutput::JetNConstituents); \ - using _table_type_##O = _table_type_##Os::iterator; \ - namespace _name_##substructure \ - { \ - DECLARE_SOA_INDEX_COLUMN(_table_type_##O, outputTable); \ - DECLARE_SOA_DYNAMIC_COLUMN(Dummy##_jet_type_, dummy##_jet_type_, []() -> int { return 0; }); \ - } \ - \ - DECLARE_SOA_TABLE(_table_type_##SSs, "AOD", _description_ "SS", jetsubstructure::EnergyMother, jetsubstructure::PtLeading, jetsubstructure::PtSubLeading, jetsubstructure::Theta, _name_##substructure::Dummy##_jet_type_<>); \ - DECLARE_SOA_TABLE(_table_type_##SSOs, "AOD", _description_ "SSO", _name_##substructure::_table_type_##OId, jetsubstructure::EnergyMother, jetsubstructure::PtLeading, jetsubstructure::PtSubLeading, jetsubstructure::Theta); \ - \ - using _table_type_##O = _table_type_##Os::iterator; \ - using _table_type_##SSO = _table_type_##SSOs::iterator; +#define JETSUBSTRUCTURE_TABLE_DEF(_collision_type_, _table_type_, _jet_type_, _matched_jet_type_, _cand_type_, _name_, _description_) \ + \ + namespace _name_##collisionoutput \ + { \ + DECLARE_SOA_DYNAMIC_COLUMN(Dummy##_collision_type_, dummy##_collision_type_, []() -> int { return 0; }); \ + } \ + \ + DECLARE_SOA_TABLE(_table_type_##COs, "AOD", _description_ "CO", jetcollision::PosZ, jetcollision::Centrality, jetcollision::EventSel, _name_##collisionoutput::Dummy##_collision_type_<>); \ + using _table_type_##CO = _table_type_##COs::iterator; \ + \ + namespace _name_##geomatched \ + { \ + DECLARE_SOA_ARRAY_INDEX_COLUMN_FULL(_matched_jet_type_, matchedJetGeo, int32_t, _matched_jet_type_##s, "_geo"); \ + } \ + \ + namespace _name_##ptmatched \ + { \ + DECLARE_SOA_ARRAY_INDEX_COLUMN_FULL(_matched_jet_type_, matchedJetPt, int32_t, _matched_jet_type_##s, "_pt"); \ + } \ + \ + namespace _name_##candmatched \ + { \ + DECLARE_SOA_ARRAY_INDEX_COLUMN_FULL(_matched_jet_type_, matchedJetCand, int32_t, _matched_jet_type_##s, "_hf"); \ + } \ + namespace _name_##jetoutput \ + { \ + DECLARE_SOA_INDEX_COLUMN(_table_type_##CO, collision); \ + DECLARE_SOA_INDEX_COLUMN_FULL(Candidate, candidate, int, _cand_type_, "_0"); \ + } \ + DECLARE_SOA_TABLE(_table_type_##Os, "AOD", _description_ "O", _name_##jetoutput::_table_type_##COId, _name_##jetoutput::CandidateId, jetoutput::JetPt, jetoutput::JetPhi, jetoutput::JetEta, jetoutput::JetR, jetoutput::JetNConstituents); \ + using _table_type_##O = _table_type_##Os::iterator; \ + namespace _name_##substructure \ + { \ + DECLARE_SOA_INDEX_COLUMN(_table_type_##O, outputTable); \ + DECLARE_SOA_DYNAMIC_COLUMN(Dummy##_jet_type_, dummy##_jet_type_, []() -> int { return 0; }); \ + } \ + \ + DECLARE_SOA_TABLE(_table_type_##SSs, "AOD", _description_ "SS", jetsubstructure::EnergyMother, jetsubstructure::PtLeading, jetsubstructure::PtSubLeading, jetsubstructure::Theta, _name_##substructure::Dummy##_jet_type_<>); \ + DECLARE_SOA_TABLE(_table_type_##SSOs, "AOD", _description_ "SSO", _name_##substructure::_table_type_##OId, jetsubstructure::EnergyMother, jetsubstructure::PtLeading, jetsubstructure::PtSubLeading, jetsubstructure::Theta); \ + DECLARE_SOA_TABLE(_table_type_##MOs, "AOD", _description_ "MO", _name_##substructure::_table_type_##OId, _name_##geomatched::_matched_jet_type_##Ids, _name_##ptmatched::_matched_jet_type_##Ids, _name_##candmatched::_matched_jet_type_##Ids); \ + \ + using _table_type_##O = _table_type_##Os::iterator; \ + using _table_type_##SSO = _table_type_##SSOs::iterator; \ + using _table_type_##MO = _table_type_##MOs::iterator; #define JETSUBSTRUCTURE_TABLES_DEF(table_type_, _jet_type_, _cand_type_, _hfparticle_type_, _description_) \ JETSUBSTRUCTURE_TABLE_DEF(JCollision, table_type_##Jet, _jet_type_##Jet, _jet_type_##EventWiseSubtractedJet, _cand_type_, _jet_type_##jet, _description_ "JET") \ diff --git a/PWGJE/TableProducer/jetderiveddataproducer.cxx b/PWGJE/TableProducer/jetderiveddataproducer.cxx index df878408be4..9d500b6b97c 100644 --- a/PWGJE/TableProducer/jetderiveddataproducer.cxx +++ b/PWGJE/TableProducer/jetderiveddataproducer.cxx @@ -112,7 +112,7 @@ struct JetDerivedDataProducerTask { void processTracks(soa::Join::iterator const& track) { - jTracksTable(track.collisionId(), track.pt(), track.eta(), track.phi(), jetderiveddatautilities::trackEnergy(track), track.sign(), jetderiveddatautilities::setTrackSelectionBit(track)); + jTracksTable(track.collisionId(), track.pt(), track.eta(), track.phi(), jetderiveddatautilities::trackEnergy(track), static_cast(track.sign()), jetderiveddatautilities::setTrackSelectionBit(track)); jTracksParentIndexTable(track.globalIndex()); } PROCESS_SWITCH(JetDerivedDataProducerTask, processTracks, "produces derived track table", true); diff --git a/PWGJE/TableProducer/jetderiveddatawriter.cxx b/PWGJE/TableProducer/jetderiveddatawriter.cxx index d44cbb10510..557b9ae8b90 100644 --- a/PWGJE/TableProducer/jetderiveddatawriter.cxx +++ b/PWGJE/TableProducer/jetderiveddatawriter.cxx @@ -16,6 +16,7 @@ /// \author Jochen Klein /// \author Nima Zardoshti +#include #include #include "Framework/AnalysisTask.h" @@ -181,7 +182,10 @@ struct JetDerivedDataWriter { storedJFullTriggerSelsTable(collision.fullTriggerSel()); for (const auto& track : tracks) { - storedJTracksTable(storedJCollisionsTable.lastIndex(), track.pt(), track.eta(), track.phi(), track.energy(), track.sign(), track.trackSel()); + if (track.trackSel() == 0) { // skips tracks that pass no selections. This might cause a problem with tracks matched with clusters. We should generate a track selection purely for cluster matched tracks so that they are kept + continue; + } + storedJTracksTable(storedJCollisionsTable.lastIndex(), o2::math_utils::detail::truncateFloatFraction(track.pt()), o2::math_utils::detail::truncateFloatFraction(track.eta()), o2::math_utils::detail::truncateFloatFraction(track.phi()), o2::math_utils::detail::truncateFloatFraction(track.energy()), track.sign(), track.trackSel()); storedJTracksParentIndexTable(track.trackId()); trackMapping.insert(std::make_pair(track.globalIndex(), storedJTracksTable.lastIndex())); } @@ -288,7 +292,7 @@ struct JetDerivedDataWriter { i++; } } - storedJMcParticlesTable(storedJMcCollisionsTable.lastIndex(), particle.pt(), particle.eta(), particle.phi(), particle.y(), particle.e(), particle.pdgCode(), particle.getGenStatusCode(), particle.getHepMCStatusCode(), particle.isPhysicalPrimary(), mothersId, daughtersId); + storedJMcParticlesTable(storedJMcCollisionsTable.lastIndex(), o2::math_utils::detail::truncateFloatFraction(particle.pt()), o2::math_utils::detail::truncateFloatFraction(particle.eta()), o2::math_utils::detail::truncateFloatFraction(particle.phi()), o2::math_utils::detail::truncateFloatFraction(particle.y()), o2::math_utils::detail::truncateFloatFraction(particle.e()), particle.pdgCode(), particle.getGenStatusCode(), particle.getHepMCStatusCode(), particle.isPhysicalPrimary(), mothersId, daughtersId); storedJParticlesParentIndexTable(particle.mcParticleId()); } @@ -350,7 +354,10 @@ struct JetDerivedDataWriter { const auto tracksPerCollision = tracks.sliceBy(TracksPerCollision, collision.globalIndex()); for (const auto& track : tracksPerCollision) { - storedJTracksTable(storedJCollisionsTable.lastIndex(), track.pt(), track.eta(), track.phi(), track.energy(), track.sign(), track.trackSel()); + if (track.trackSel() == 0) { // skips tracks that pass no selections. This might cause a problem with tracks matched with clusters. We should generate a track selection purely for cluster matched tracks so that they are kept + continue; + } + storedJTracksTable(storedJCollisionsTable.lastIndex(), o2::math_utils::detail::truncateFloatFraction(track.pt()), o2::math_utils::detail::truncateFloatFraction(track.eta()), o2::math_utils::detail::truncateFloatFraction(track.phi()), o2::math_utils::detail::truncateFloatFraction(track.energy()), track.sign(), track.trackSel()); storedJTracksParentIndexTable(track.trackId()); if (track.has_mcParticle()) { @@ -461,7 +468,7 @@ struct JetDerivedDataWriter { i++; } } - storedJMcParticlesTable(storedJMcCollisionsTable.lastIndex(), particle.pt(), particle.eta(), particle.phi(), particle.y(), particle.e(), particle.pdgCode(), particle.getGenStatusCode(), particle.getHepMCStatusCode(), particle.isPhysicalPrimary(), mothersId, daughtersId); + storedJMcParticlesTable(storedJMcCollisionsTable.lastIndex(), o2::math_utils::detail::truncateFloatFraction(particle.pt()), o2::math_utils::detail::truncateFloatFraction(particle.eta()), o2::math_utils::detail::truncateFloatFraction(particle.phi()), o2::math_utils::detail::truncateFloatFraction(particle.y()), o2::math_utils::detail::truncateFloatFraction(particle.e()), particle.pdgCode(), particle.getGenStatusCode(), particle.getHepMCStatusCode(), particle.isPhysicalPrimary(), mothersId, daughtersId); storedJParticlesParentIndexTable(particle.mcParticleId()); } diff --git a/PWGJE/Tasks/jetsubstructurehfoutput.cxx b/PWGJE/Tasks/jetsubstructurehfoutput.cxx index 690f16bffc0..40aa2809a44 100644 --- a/PWGJE/Tasks/jetsubstructurehfoutput.cxx +++ b/PWGJE/Tasks/jetsubstructurehfoutput.cxx @@ -40,20 +40,24 @@ using namespace o2::framework::expressions; // NB: runDataProcessing.h must be included after customize! #include "Framework/runDataProcessing.h" -template +template struct JetSubstructureHFOutputTask { Produces collisionOutputTableData; Produces jetOutputTableData; Produces jetSubstructureOutputTableData; + Produces jetMatchingOutputTableData; Produces collisionOutputTableDataSub; Produces jetOutputTableDataSub; Produces jetSubstructureOutputTableDataSub; + Produces jetMatchingOutputTableDataSub; Produces collisionOutputTableMCD; Produces jetOutputTableMCD; Produces jetSubstructureOutputTableMCD; + Produces jetMatchingOutputTableMCD; Produces collisionOutputTableMCP; Produces jetOutputTableMCP; Produces jetSubstructureOutputTableMCP; + Produces jetMatchingOutputTableMCP; Produces hfCollisionsTable; Produces candidateTable; Produces candidateParsTable; @@ -73,6 +77,11 @@ struct JetSubstructureHFOutputTask { std::map candidateMapping; std::map candidateCollisionMapping; + std::map jetMappingData; + std::map jetMappingDataSub; + std::map jetMappingMCD; + std::map jetMappingMCP; + std::vector jetRadiiValues; void init(InitContext const&) @@ -82,8 +91,8 @@ struct JetSubstructureHFOutputTask { Filter jetSelection = aod::jet::pt >= jetPtMin; - template - void fillTables(T const& jet, U const& cand, int32_t collisionIndex, int32_t candidateIndex, V& collisionOutputTable, M& jetOutputTable, N& jetSubstructureOutputTable, std::vector geoMatching, std::vector ptMatching, std::vector candMatching) + template + void fillTables(T const& jet, U const& cand, int32_t collisionIndex, int32_t candidateIndex, V& jetOutputTable, M& jetSubstructureOutputTable, std::map& jetMapping) { std::vector energyMotherVec; std::vector ptLeadingVec; @@ -97,12 +106,13 @@ struct JetSubstructureHFOutputTask { std::copy(ptLeadingSpan.begin(), ptLeadingSpan.end(), std::back_inserter(ptLeadingVec)); std::copy(ptSubLeadingSpan.begin(), ptSubLeadingSpan.end(), std::back_inserter(ptSubLeadingVec)); std::copy(thetaSpan.begin(), thetaSpan.end(), std::back_inserter(thetaVec)); - jetOutputTable(collisionIndex, candidateIndex, geoMatching, ptMatching, candMatching, jet.pt(), jet.phi(), jet.eta(), jet.r(), jet.tracksIds().size() + jet.hfcandidatesIds().size()); // here we take the decision to keep the collision index consistent with the JE framework in case it is later needed to join to other tables. The candidate Index however can be linked to the HF tables + jetOutputTable(collisionIndex, candidateIndex, jet.pt(), jet.phi(), jet.eta(), jet.r(), jet.tracksIds().size() + jet.hfcandidatesIds().size()); // here we take the decision to keep the collision index consistent with the JE framework in case it is later needed to join to other tables. The candidate Index however can be linked to the HF tables jetSubstructureOutputTable(jetOutputTable.lastIndex(), energyMotherVec, ptLeadingVec, ptSubLeadingVec, thetaVec); + jetMapping.insert(std::make_pair(jet.globalIndex(), jetOutputTable.lastIndex())); } - template - void analyseCharged(T const& collision, U const& jets, V const& jetsTag, M const& tracks, N const& candidates, O& collisionOutputTable, P& jetOutputTable, S& jetSubstructureOutputTable) + template + void analyseCharged(T const& collision, U const& jets, V const& candidates, M& collisionOutputTable, N& jetOutputTable, O& jetSubstructureOutputTable, std::map& jetMapping) { int nJetInCollision = 0; @@ -113,39 +123,20 @@ struct JetSubstructureHFOutputTask { } for (const auto& jetRadiiValue : jetRadiiValues) { if (jet.r() == round(jetRadiiValue * 100.0f)) { - auto candidate = jet.template hfcandidates_first_as(); - std::vector geoMatching; - std::vector ptMatching; - std::vector hfMatching; - if constexpr (isMatched) { - if (jet.has_matchedJetGeo()) { - for (auto& jetTag : jet.template matchedJetGeo_as()) { - geoMatching.push_back(jetTag.globalIndex()); - } - } - if (jet.has_matchedJetPt()) { - for (auto& jetTag : jet.template matchedJetPt_as()) { - ptMatching.push_back(jetTag.globalIndex()); - } - } - if (jet.has_matchedJetCand()) { - for (auto& jetTag : jet.template matchedJetCand_as()) { - hfMatching.push_back(jetTag.globalIndex()); - } - } - } + auto candidate = jet.template hfcandidates_first_as(); int32_t candidateIndex = -1; auto candidateTableIndex = candidateMapping.find(candidate.globalIndex()); if (candidateTableIndex != candidateMapping.end()) { candidateIndex = candidateTableIndex->second; } - - if (nJetInCollision == 0) { - collisionOutputTable(collision.posZ(), collision.centrality(), collision.eventSel()); - collisionIndex = collisionOutputTable.lastIndex(); + if constexpr (!isMc) { + if (nJetInCollision == 0) { + collisionOutputTable(collision.posZ(), collision.centrality(), collision.eventSel()); + collisionIndex = collisionOutputTable.lastIndex(); + } + nJetInCollision++; } - nJetInCollision++; - fillTables(jet, candidate, collisionIndex, candidateIndex, collisionOutputTable, jetOutputTable, jetSubstructureOutputTable, geoMatching, ptMatching, hfMatching); + fillTables(jet, candidate, collisionIndex, candidateIndex, jetOutputTable, jetSubstructureOutputTable, jetMapping); } } } @@ -190,6 +181,53 @@ struct JetSubstructureHFOutputTask { } } + template + void analyseMatched(T const& jets, U const& jetsTag, std::map& jetMapping, std::map& jetTagMapping, V& matchingOutputTable) + { + for (const auto& jet : jets) { + if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { + continue; + } + for (const auto& jetRadiiValue : jetRadiiValues) { + if (jet.r() == round(jetRadiiValue * 100.0f)) { + std::vector geoMatching; + std::vector ptMatching; + std::vector candMatching; + if (jet.has_matchedJetGeo()) { + for (auto& jetTagId : jet.matchedJetGeoIds()) { + auto jetTagIndex = jetTagMapping.find(jetTagId); + if (jetTagIndex != jetTagMapping.end()) { + geoMatching.push_back(jetTagIndex->second); + } + } + } + if (jet.has_matchedJetPt()) { + for (auto& jetTagId : jet.matchedJetPtIds()) { + auto jetTagIndex = jetTagMapping.find(jetTagId); + if (jetTagIndex != jetTagMapping.end()) { + ptMatching.push_back(jetTagIndex->second); + } + } + } + if (jet.has_matchedJetCand()) { + for (auto& jetTagId : jet.matchedJetCandIds()) { + auto jetTagIndex = jetTagMapping.find(jetTagId); + if (jetTagIndex != jetTagMapping.end()) { + candMatching.push_back(jetTagIndex->second); + } + } + } + int storedJetIndex = -1; + auto jetIndex = jetMapping.find(jet.globalIndex()); + if (jetIndex != jetMapping.end()) { + storedJetIndex = jetIndex->second; + } + matchingOutputTable(storedJetIndex, geoMatching, ptMatching, candMatching); + } + } + } + } + void processDummy(JetCollisions const& collisions) {} PROCESS_SWITCH(JetSubstructureHFOutputTask, processDummy, "Dummy process function turned on by default", true); @@ -225,81 +263,55 @@ struct JetSubstructureHFOutputTask { void processOutputData(JetCollision const& collision, JetTableData const& jets, - JetTracks const& tracks, - CandidateCollisionTable const& canidateCollisions, CandidateTable const& candidates) { - analyseCharged(collision, jets, jets, tracks, candidates, collisionOutputTableData, jetOutputTableData, jetSubstructureOutputTableData); + analyseCharged(collision, jets, candidates, collisionOutputTableData, jetOutputTableData, jetSubstructureOutputTableData, jetMappingData); } PROCESS_SWITCH(JetSubstructureHFOutputTask, processOutputData, "hf jet substructure output Data", false); void processOutputDataSub(JetCollision const& collision, - JetMatchedTableData const& jets, - JetTableDataSub const& jetsSub, - TracksSub const& tracks, - CandidateCollisionTable const& canidateCollisions, + JetTableDataSub const& jets, CandidateTable const& candidates) { - analyseCharged(collision, jets, jetsSub, tracks, candidates, collisionOutputTableData, jetOutputTableData, jetSubstructureOutputTableData); - analyseCharged(collision, jetsSub, jets, tracks, candidates, collisionOutputTableDataSub, jetOutputTableDataSub, jetSubstructureOutputTableDataSub); + analyseCharged(collision, jets, candidates, collisionOutputTableDataSub, jetOutputTableDataSub, jetSubstructureOutputTableDataSub, jetMappingDataSub); } PROCESS_SWITCH(JetSubstructureHFOutputTask, processOutputDataSub, "hf jet substructure output event-wise subtracted Data", false); + void processOutputMatchingData(JetMatchedTableData const& jets, + JetTableDataSub const& jetsSub) + { + analyseMatched(jets, jetsSub, jetMappingData, jetMappingDataSub, jetMatchingOutputTableData); + analyseMatched(jetsSub, jets, jetMappingDataSub, jetMappingData, jetMatchingOutputTableDataSub); + } + PROCESS_SWITCH(JetSubstructureHFOutputTask, processOutputMatchingData, "jet matching output Data", false); + void processOutputMCD(JetCollision const& collision, JetTableMCD const& jets, - JetTableMCP const& jetsTag, - JetTracks const& tracks, - CandidateCollisionTable const& canidateCollisions, CandidateTableMCD const& candidates) { - analyseCharged(collision, jets, jetsTag, tracks, candidates, collisionOutputTableMCD, jetOutputTableMCD, jetSubstructureOutputTableMCD); + analyseCharged(collision, jets, candidates, collisionOutputTableMCD, jetOutputTableMCD, jetSubstructureOutputTableMCD, jetMappingMCD); } PROCESS_SWITCH(JetSubstructureHFOutputTask, processOutputMCD, "hf jet substructure output MCD", false); void processOutputMCP(JetMcCollision const& collision, JetTableMCP const& jets, - JetTableMCD const& jetsTag, - JetParticles const& particles, CandidateTableMCP const& candidates) { - for (const auto& jet : jets) { - if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { - continue; - } - for (const auto& jetRadiiValue : jetRadiiValues) { - if (jet.r() == round(jetRadiiValue * 100.0f)) { - auto candidate = jet.template hfcandidates_first_as(); - std::vector geoMatching; - std::vector ptMatching; - std::vector hfMatching; - if (jet.has_matchedJetGeo()) { - for (auto& jetTag : jet.template matchedJetGeo_as()) { - geoMatching.push_back(jetTag.globalIndex()); - } - } - if (jet.has_matchedJetPt()) { - for (auto& jetTag : jet.template matchedJetPt_as()) { - ptMatching.push_back(jetTag.globalIndex()); - } - } - if (jet.has_matchedJetCand()) { - for (auto& jetTag : jet.template matchedJetCand_as()) { - hfMatching.push_back(jetTag.globalIndex()); - } - } - int32_t candidateIndex = -1; - jethfutilities::fillCandidateMcTable(candidate, hfParticlesTable, candidateIndex); - - fillTables(jet, candidate, -1, candidateIndex, collisionOutputTableMCP, jetOutputTableMCP, jetSubstructureOutputTableMCP, geoMatching, ptMatching, hfMatching); - } - } - } + analyseCharged(collision, jets, candidates, collisionOutputTableMCP, jetOutputTableMCP, jetSubstructureOutputTableMCP, jetMappingMCP); } PROCESS_SWITCH(JetSubstructureHFOutputTask, processOutputMCP, "hf jet substructure output MCP", false); + + void processOutputMatchingMC(JetTableMCD const& jetsMCD, + JetTableMCP const& jetsMCP) + { + analyseMatched(jetsMCD, jetsMCP, jetMappingMCD, jetMappingMCP, jetMatchingOutputTableMCD); + analyseMatched(jetsMCP, jetsMCD, jetMappingMCP, jetMappingMCD, jetMatchingOutputTableMCP); + } + PROCESS_SWITCH(JetSubstructureHFOutputTask, processOutputMatchingMC, "jet matching output MC", false); }; -using JetSubstructureOutputD0 = JetSubstructureHFOutputTask>, soa::Filtered>, aod::D0CJetCOs, aod::D0CJetOs, aod::D0CJetSSOs, soa::Filtered>, aod::D0CMCDJetCOs, aod::D0CMCDJetOs, aod::D0CMCDJetSSOs, soa::Filtered>, aod::D0CMCPJetCOs, aod::D0CMCPJetOs, aod::D0CMCPJetSSOs, soa::Filtered>, aod::D0CEWSJetCOs, aod::D0CEWSJetOs, aod::D0CEWSJetSSOs, aod::StoredHfD0CollBase, aod::StoredHfD0Bases, aod::StoredHfD0Pars, aod::StoredHfD0ParEs, aod::StoredHfD0Sels, aod::StoredHfD0Mcs, aod::StoredHfD0PBases>; -// using JetSubstructureOutputLc = JetSubstructureHFOutputTask>, soa::Filtered>, aod::LcCJetCOs, aod::LcCJetOs, aod::LcCJetSSOs, soa::Filtered>, aod::LcCMCDJetCOs, aod::LcMCDJetOs, aod::LcCMCDJetSSOs, soa::Filtered>, aod::LcCMCPJetCOs, aod::LcCMCPJetOs, aod::LcCMCPJetSSOs, soa::Filtered>, aod::LcCEWSJetCOs, aod::LcCEWSJetOs, aod::LcCEWSJetSSOs, aod::StoredHfLcCollBase, aod::StoredHfLcBases, aod::StoredHfLcPars, aod::StoredHfLcParEs, aod::StoredHfLcSels, aod::StoredHfLcMcs, aod::StoredHfLcPBases>; -// using JetSubstructureOutputBplus = JetSubstructureHFOutputTask>, soa::Filtered>, aod::BplusCJetCOs, aod::BplusCJetOs, aod::BplusCJetSSOs, soa::Filtered>, aod::BplusCMCDJetCOs, aod::BplusCMCDJetOs, aod::BplusCMCDJetSSOs, soa::Filtered>, aod::BplusCMCPJetCOs, aod::BplusCMCPJetOs, aod::BplusCMCPJetSSOs, soa::Filtered>, aod::BplusCEWSJetCOs, aod::BplusCEWSJetOs, aod::BplusCEWSJetSSOs, aod::StoredHfBplusCollBase, aod::StoredHfBplusBases, aod::StoredHfBplusPars, aod::StoredHfBplusParEs, aod::StoredHfBplusSels, aod::StoredHfBplusMcs, aod::StoredHfBplusPBases>; +using JetSubstructureOutputD0 = JetSubstructureHFOutputTask>, soa::Filtered>, aod::D0CJetCOs, aod::D0CJetOs, aod::D0CJetSSOs, aod::D0CJetMOs, soa::Filtered>, aod::D0CMCDJetCOs, aod::D0CMCDJetOs, aod::D0CMCDJetSSOs, aod::D0CMCDJetMOs, soa::Filtered>, aod::D0CMCPJetCOs, aod::D0CMCPJetOs, aod::D0CMCPJetSSOs, aod::D0CMCPJetMOs, soa::Filtered>, aod::D0CEWSJetCOs, aod::D0CEWSJetOs, aod::D0CEWSJetSSOs, aod::D0CEWSJetMOs, aod::StoredHfD0CollBase, aod::StoredHfD0Bases, aod::StoredHfD0Pars, aod::StoredHfD0ParEs, aod::StoredHfD0Sels, aod::StoredHfD0Mcs, aod::StoredHfD0PBases>; +// using JetSubstructureOutputLc = JetSubstructureHFOutputTask>, soa::Filtered>, aod::LcCJetCOs, aod::LcCJetOs, aod::LcCJetSSOs, aod::LcCJetMOs, soa::Filtered>, aod::LcCMCDJetCOs, aod::LcMCDJetOs, aod::LcCMCDJetSSOs, aod::LcCMCDJetMOs, soa::Filtered>, aod::LcCMCPJetCOs, aod::LcCMCPJetOs, aod::LcCMCPJetSSOs, aod::LcCMCPJetMOs, soa::Filtered>, aod::LcCEWSJetCOs, aod::LcCEWSJetOs, aod::LcCEWSJetSSOs, aod::LcCEWSJetMOs, aod::StoredHfLcCollBase, aod::StoredHfLcBases, aod::StoredHfLcPars, aod::StoredHfLcParEs, aod::StoredHfLcSels, aod::StoredHfLcMcs, aod::StoredHfLcPBases>; +// using JetSubstructureOutputBplus = JetSubstructureHFOutputTask>, soa::Filtered>, aod::BplusCJetCOs, aod::BplusCJetOs, aod::BplusCJetSSOs, aod::BplusCJetMOs, soa::Filtered>, aod::BplusCMCDJetCOs, aod::BplusCMCDJetOs, aod::BplusCMCDJetSSOs, aod::BplusCMCDJetMOs, soa::Filtered>, aod::BplusCMCPJetCOs, aod::BplusCMCPJetOs, aod::BplusCMCPJetSSOs, aod::BplusCMCPJetMOs, soa::Filtered>, aod::BplusCEWSJetCOs, aod::BplusCEWSJetOs, aod::BplusCEWSJetSSOs, aod::BplusCEWSJetMOs, aod::StoredHfBplusCollBase, aod::StoredHfBplusBases, aod::StoredHfBplusPars, aod::StoredHfBplusParEs, aod::StoredHfBplusSels, aod::StoredHfBplusMcs, aod::StoredHfBplusPBases>; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { diff --git a/PWGJE/Tasks/jetsubstructureoutput.cxx b/PWGJE/Tasks/jetsubstructureoutput.cxx index 09abdda8ff4..8d424c7db5b 100644 --- a/PWGJE/Tasks/jetsubstructureoutput.cxx +++ b/PWGJE/Tasks/jetsubstructureoutput.cxx @@ -42,15 +42,19 @@ struct JetSubstructureOutputTask { Produces collisionOutputTableData; Produces jetOutputTableData; Produces jetSubstructureOutputTableData; + Produces jetMatchingOutputTableData; Produces collisionOutputTableDataSub; Produces jetOutputTableDataSub; Produces jetSubstructureOutputTableDataSub; + Produces jetMatchingOutputTableDataSub; Produces collisionOutputTableMCD; Produces jetOutputTableMCD; Produces jetSubstructureOutputTableMCD; + Produces jetMatchingOutputTableMCD; Produces collisionOutputTableMCP; Produces jetOutputTableMCP; Produces jetSubstructureOutputTableMCP; + Produces jetMatchingOutputTableMCP; Configurable jetPtMin{"jetPtMin", 0.0, "minimum jet pT cut"}; Configurable> jetRadii{"jetRadii", std::vector{0.4}, "jet resolution parameters"}; @@ -60,6 +64,11 @@ struct JetSubstructureOutputTask { Configurable trackEtaMin{"trackEtaMin", -0.9, "minimum track pseudorapidity"}; Configurable trackEtaMax{"trackEtaMax", 0.9, "maximum track pseudorapidity"}; + std::map jetMappingData; + std::map jetMappingDataSub; + std::map jetMappingMCD; + std::map jetMappingMCP; + std::vector jetRadiiValues; void init(InitContext const&) @@ -69,8 +78,8 @@ struct JetSubstructureOutputTask { Filter jetSelection = aod::jet::pt >= jetPtMin; - template - void fillTables(T const& jet, int32_t collisionIndex, U& collisionOutputTable, V& jetOutputTable, M& jetSubstructureOutputTable, std::vector geoMatching, std::vector ptMatching, std::vector candMatching) + template + void fillTables(T const& jet, int32_t collisionIndex, U& jetOutputTable, V& jetSubstructureOutputTable, std::map& jetMapping) { std::vector energyMotherVec; std::vector ptLeadingVec; @@ -84,17 +93,40 @@ struct JetSubstructureOutputTask { std::copy(ptLeadingSpan.begin(), ptLeadingSpan.end(), std::back_inserter(ptLeadingVec)); std::copy(ptSubLeadingSpan.begin(), ptSubLeadingSpan.end(), std::back_inserter(ptSubLeadingVec)); std::copy(thetaSpan.begin(), thetaSpan.end(), std::back_inserter(thetaVec)); - jetOutputTable(collisionIndex, -1, geoMatching, ptMatching, candMatching, jet.pt(), jet.phi(), jet.eta(), jet.r(), jet.tracksIds().size()); + jetOutputTable(collisionIndex, -1, jet.pt(), jet.phi(), jet.eta(), jet.r(), jet.tracksIds().size()); jetSubstructureOutputTable(jetOutputTable.lastIndex(), energyMotherVec, ptLeadingVec, ptSubLeadingVec, thetaVec); + jetMapping.insert(std::make_pair(jet.globalIndex(), jetOutputTable.lastIndex())); } - template - void analyseCharged(T const& collision, U const& jets, V const& jetsTag, M& collisionOutputTable, N& jetOutputTable, O& jetSubstructureOutputTable) + template + void analyseCharged(T const& collision, U const& jets, V& collisionOutputTable, M& jetOutputTable, N& jetSubstructureOutputTable, std::map& jetMapping) { - std::vector candMatching{-1}; int nJetInCollision = 0; int32_t collisionIndex = -1; + for (const auto& jet : jets) { + if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { + continue; + } + for (const auto& jetRadiiValue : jetRadiiValues) { + if (jet.r() == round(jetRadiiValue * 100.0f)) { + if constexpr (!isMc) { + if (nJetInCollision == 0) { + collisionOutputTable(collision.posZ(), collision.centrality(), collision.eventSel()); + collisionIndex = collisionOutputTable.lastIndex(); + } + nJetInCollision++; + } + fillTables(jet, collisionIndex, jetOutputTable, jetSubstructureOutputTable, jetMapping); + } + } + } + } + + template + void analyseMatched(T const& jets, U const& jetsTag, std::map& jetMapping, std::map& jetTagMapping, V& matchingOutputTable) + { + std::vector candMatching; for (const auto& jet : jets) { if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { continue; @@ -103,24 +135,28 @@ struct JetSubstructureOutputTask { if (jet.r() == round(jetRadiiValue * 100.0f)) { std::vector geoMatching; std::vector ptMatching; - if constexpr (hasMatching) { - if (jet.has_matchedJetGeo()) { - for (auto& jetTag : jet.template matchedJetGeo_as()) { - geoMatching.push_back(jetTag.globalIndex()); + if (jet.has_matchedJetGeo()) { + for (auto& jetTagId : jet.matchedJetGeoIds()) { + auto jetTagIndex = jetTagMapping.find(jetTagId); + if (jetTagIndex != jetTagMapping.end()) { + geoMatching.push_back(jetTagIndex->second); } } + } if (jet.has_matchedJetPt()) { - for (auto& jetTag : jet.template matchedJetPt_as()) { - ptMatching.push_back(jetTag.globalIndex()); + for (auto& jetTagId : jet.matchedJetPtIds()) { + auto jetTagIndex = jetTagMapping.find(jetTagId); + if (jetTagIndex != jetTagMapping.end()) { + ptMatching.push_back(jetTagIndex->second); + } } } - } - if (nJetInCollision == 0) { - collisionOutputTable(collision.posZ(), collision.centrality(), collision.eventSel()); - collisionIndex = collisionOutputTable.lastIndex(); - } - nJetInCollision++; - fillTables(jet, collisionIndex, collisionOutputTable, jetOutputTable, jetSubstructureOutputTable, geoMatching, ptMatching, candMatching); + int storedJetIndex = -1; + auto jetIndex = jetMapping.find(jet.globalIndex()); + if (jetIndex != jetMapping.end()) { + storedJetIndex = jetIndex->second; + } + matchingOutputTable(storedJetIndex, geoMatching, ptMatching, candMatching); } } } @@ -130,62 +166,48 @@ struct JetSubstructureOutputTask { PROCESS_SWITCH(JetSubstructureOutputTask, processDummy, "Dummy process function turned on by default", true); void processOutputData(JetCollision const& collision, - soa::Filtered> const& jets, - JetTracks const& tracks) + soa::Filtered> const& jets) { - analyseCharged(collision, jets, jets, collisionOutputTableData, jetOutputTableData, jetSubstructureOutputTableData); + analyseCharged(collision, jets, collisionOutputTableData, jetOutputTableData, jetSubstructureOutputTableData, jetMappingData); } PROCESS_SWITCH(JetSubstructureOutputTask, processOutputData, "jet substructure output Data", false); void processOutputDataSub(JetCollision const& collision, - soa::Filtered> const& jets, - soa::Filtered> const& jetsSub, - JetTracks const& tracks) + soa::Filtered> const& jets) { - analyseCharged(collision, jets, jetsSub, collisionOutputTableData, jetOutputTableData, jetSubstructureOutputTableData); - analyseCharged(collision, jetsSub, jets, collisionOutputTableDataSub, jetOutputTableDataSub, jetSubstructureOutputTableDataSub); + analyseCharged(collision, jets, collisionOutputTableDataSub, jetOutputTableDataSub, jetSubstructureOutputTableDataSub, jetMappingDataSub); } PROCESS_SWITCH(JetSubstructureOutputTask, processOutputDataSub, "jet substructure output event-wise subtracted Data", false); + void processOutputMatchingData(soa::Filtered> const& jets, + soa::Filtered> const& jetsSub) + { + analyseMatched(jets, jetsSub, jetMappingData, jetMappingDataSub, jetMatchingOutputTableData); + analyseMatched(jetsSub, jets, jetMappingDataSub, jetMappingData, jetMatchingOutputTableDataSub); + } + PROCESS_SWITCH(JetSubstructureOutputTask, processOutputMatchingData, "jet matching output Data", false); + void processOutputMCD(JetCollision const& collision, - soa::Filtered> const& jets, - aod::ChargedMCParticleLevelJets const& jetsTag, - JetTracks const& tracks) + soa::Filtered> const& jets) { - analyseCharged(collision, jets, jetsTag, collisionOutputTableMCD, jetOutputTableMCD, jetSubstructureOutputTableMCD); + analyseCharged(collision, jets, collisionOutputTableMCD, jetOutputTableMCD, jetSubstructureOutputTableMCD, jetMappingMCD); } PROCESS_SWITCH(JetSubstructureOutputTask, processOutputMCD, "jet substructure output MCD", false); void processOutputMCP(JetMcCollision const& collision, - soa::Filtered> const& jets, - aod::ChargedMCDetectorLevelJets const& jetsTag, - JetParticles const& particles) + soa::Filtered> const& jets) { - std::vector candMatching{-1}; - for (const auto& jet : jets) { - if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { - continue; - } - for (const auto& jetRadiiValue : jetRadiiValues) { - if (jet.r() == round(jetRadiiValue * 100.0f)) { - std::vector geoMatching; - std::vector ptMatching; - if (jet.has_matchedJetGeo()) { - for (auto& jetTag : jet.template matchedJetGeo_as()) { - geoMatching.push_back(jetTag.globalIndex()); - } - } - if (jet.has_matchedJetPt()) { - for (auto& jetTag : jet.template matchedJetPt_as()) { - ptMatching.push_back(jetTag.globalIndex()); - } - } - fillTables(jet, -1, collisionOutputTableMCP, jetOutputTableMCP, jetSubstructureOutputTableMCP, geoMatching, ptMatching, candMatching); - } - } - } + analyseCharged(collision, jets, collisionOutputTableMCP, jetOutputTableMCP, jetSubstructureOutputTableMCP, jetMappingMCP); } PROCESS_SWITCH(JetSubstructureOutputTask, processOutputMCP, "jet substructure output MCP", false); + + void processOutputMatchingMC(soa::Filtered> const& jetsMCD, + soa::Filtered> const& jetsMCP) + { + analyseMatched(jetsMCD, jetsMCP, jetMappingMCD, jetMappingMCP, jetMatchingOutputTableMCD); + analyseMatched(jetsMCP, jetsMCD, jetMappingMCP, jetMappingMCD, jetMatchingOutputTableMCP); + } + PROCESS_SWITCH(JetSubstructureOutputTask, processOutputMatchingMC, "jet matching output MC", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From 987289b8f6242bede2163296b031e80dce05bf97 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Mon, 19 Feb 2024 23:51:29 +0000 Subject: [PATCH 8/8] MegaLinter fixes --- PWGJE/Tasks/jetsubstructureoutput.cxx | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/PWGJE/Tasks/jetsubstructureoutput.cxx b/PWGJE/Tasks/jetsubstructureoutput.cxx index 8d424c7db5b..257d8b4a447 100644 --- a/PWGJE/Tasks/jetsubstructureoutput.cxx +++ b/PWGJE/Tasks/jetsubstructureoutput.cxx @@ -143,20 +143,20 @@ struct JetSubstructureOutputTask { } } } - if (jet.has_matchedJetPt()) { - for (auto& jetTagId : jet.matchedJetPtIds()) { - auto jetTagIndex = jetTagMapping.find(jetTagId); - if (jetTagIndex != jetTagMapping.end()) { - ptMatching.push_back(jetTagIndex->second); - } + if (jet.has_matchedJetPt()) { + for (auto& jetTagId : jet.matchedJetPtIds()) { + auto jetTagIndex = jetTagMapping.find(jetTagId); + if (jetTagIndex != jetTagMapping.end()) { + ptMatching.push_back(jetTagIndex->second); } } - int storedJetIndex = -1; - auto jetIndex = jetMapping.find(jet.globalIndex()); - if (jetIndex != jetMapping.end()) { - storedJetIndex = jetIndex->second; - } - matchingOutputTable(storedJetIndex, geoMatching, ptMatching, candMatching); + } + int storedJetIndex = -1; + auto jetIndex = jetMapping.find(jet.globalIndex()); + if (jetIndex != jetMapping.end()) { + storedJetIndex = jetIndex->second; + } + matchingOutputTable(storedJetIndex, geoMatching, ptMatching, candMatching); } } }