From 3912c43599863b9373dde7e872a0d371c231fb21 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Fri, 1 Aug 2025 19:51:48 +0000 Subject: [PATCH] Please consider the following formatting changes --- PWGJE/Tasks/fullJetSpectra.cxx | 3651 ++++++++++++++++---------------- 1 file changed, 1833 insertions(+), 1818 deletions(-) diff --git a/PWGJE/Tasks/fullJetSpectra.cxx b/PWGJE/Tasks/fullJetSpectra.cxx index 045b303d11d..16ecbc5c892 100644 --- a/PWGJE/Tasks/fullJetSpectra.cxx +++ b/PWGJE/Tasks/fullJetSpectra.cxx @@ -22,10 +22,10 @@ #include "PWGJE/DataModel/JetReducedData.h" #include "Common/CCDB/TriggerAliases.h" -#include "CCDB/BasicCCDBManager.h" #include "EventFiltering/Zorro.h" #include "EventFiltering/ZorroSummary.h" +#include "CCDB/BasicCCDBManager.h" #include "Framework/ASoA.h" #include "Framework/AnalysisTask.h" #include "Framework/HistogramRegistry.h" @@ -68,7 +68,7 @@ struct FullJetSpectra { Configurable doEMCALEventWorkaround{"doEMCALEventWorkaround", false, "apply the workaround to read the EMC trigger bit by requiring a cell content in the EMCAL"}; Configurable doMBGapTrigger{"doMBGapTrigger", true, "set to true only when using MB-Gap Trigger JJ MC to reject MB events at the collision and track level"}; - //Software Trigger configurables + // Software Trigger configurables Configurable doSoftwareTriggerSelection{"doSoftwareTriggerSelection", false, "set to true when using triggered datasets"}; Configurable triggerMasks{"triggerMasks", "fJetFullHighPt", "possible JE Trigger masks: fJetChLowPt,fJetChHighPt,fEMCALReadout,fJetFullHighPt,fJetFullLowPt"}; @@ -83,7 +83,7 @@ struct FullJetSpectra { Configurable jetPhiMax{"jetPhiMax", 2.86, "maximum jet phi"}; // phi_jet_min for R = 0.4 is 2.86 Configurable jetAreaFractionMin{"jetAreaFractionMin", -99.0, "used to make a cut on the jet areas"}; - //Leading track and cluster pT configurables + // Leading track and cluster pT configurables Configurable leadingTrackPtMin{"leadingTrackPtMin", -99.0, "minimum pT selection on jet tracks"}; Configurable leadingTrackPtMax{"leadingTrackPtMax", 999.0, "maximum pT selection on jet tracks"}; Configurable leadingClusterPtMin{"leadingClusterPtMin", -99.0, "minimum pT selection on jet clusters"}; @@ -136,23 +136,23 @@ struct FullJetSpectra { Zorro zorro; OutputObj zorroSummary{"zorroSummary"}; - //Multiplicity Utilities - // struct CentClass { - // const char* name; - // float min; - // float max; - // }; - // // Define multiplicity classes here (example: MB(0-100), HM(0-1), 1-10, 10-20, 20-40, 40-60, 60-100) - // static constexpr int nCentClasses = 4; - // CentClass centClasses[nCentClasses] = { - // {"MB", 0.0, 100.0}, - // {"HM", 0.0, 1.0}, - // {"1_10", 1.0, 10.0}, - // {"10_20", 10.0, 20.0}, - // {"20_40", 20.0, 40.0}, - // {"40_60", 40.0, 60.0}, - // {"60_100", 60.0, 100.0} - // }; + // Multiplicity Utilities + // struct CentClass { + // const char* name; + // float min; + // float max; + // }; + // // Define multiplicity classes here (example: MB(0-100), HM(0-1), 1-10, 10-20, 20-40, 40-60, 60-100) + // static constexpr int nCentClasses = 4; + // CentClass centClasses[nCentClasses] = { + // {"MB", 0.0, 100.0}, + // {"HM", 0.0, 1.0}, + // {"1_10", 1.0, 10.0}, + // {"10_20", 10.0, 20.0}, + // {"20_40", 20.0, 40.0}, + // {"40_60", 40.0, 60.0}, + // {"60_100", 60.0, 100.0} + // }; // Random splitter instance /* TRandom3 randGen; @@ -211,7 +211,7 @@ struct FullJetSpectra { hDetcollisionCounter->GetXaxis()->SetBinLabel(8, "EMCAcceptedDetColl"); } - if(doprocessJetsTriggeredData) { + if (doprocessJetsTriggeredData) { auto hDetTrigcollisionCounter = registry.get(HIST("hDetTrigcollisionCounter")); hDetTrigcollisionCounter->GetXaxis()->SetBinLabel(1, "allDetTrigColl"); hDetTrigcollisionCounter->GetXaxis()->SetBinLabel(2, "DetTrigCollAfterZorroSelection"); @@ -312,1688 +312,1685 @@ struct FullJetSpectra { hSpliteventSelector->GetXaxis()->SetBinLabel(3, "MatchedforRM"); } */ -void init(o2::framework::InitContext&) -{ + void init(o2::framework::InitContext&) + { - trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast(trackSelections)); - eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventSelections)); - triggerMaskBits = jetderiveddatautilities::initialiseTriggerMaskBits(triggerMasks); - // triggerMaskBits = jetderiveddatautilities::initialiseTriggerMaskBits(jetderiveddatautilities::JTriggerMasks); - particleSelection = static_cast(particleSelections); - jetRadiiValues = (std::vector)jetRadii; - - /* if (doMcClosure) { - // randGen.SetSeed(mcSplitSeed); - // randGen.SetSeed(static_cast(std::time(nullptr))); - // int seed = mcSplitSeed >= 0 ? mcSplitSeed : static_cast(std::time(nullptr)); - // randGen.SetSeed(seed); - // LOGF(info, "MC closure seed = %d", seed); - - int seed = mcSplitSeed >= 0 ? mcSplitSeed : static_cast(std::time(nullptr)); - randGen.SetSeed(seed); - LOGF(info, "MC closure splitting enabled with seed = %d, split fraction = %.2f", seed, static_cast(mcClosureSplitFrac)); - - registry.add("hSpliteventSelector", "Random MC Split Selector;Split Type;Entries",{HistType::kTH1F, {{3, 0.0, 3.0}}}); // 0=MCD, 1=MCP, 2=RM - - //individual processes' event counters for sanity checks - registry.add("h_MCD_splitevent_counter", "Events into MCD split", {HistType::kTH1F, {{1, 0.0, 1.0}}}); - registry.add("h_MCP_splitevent_counter", "Events into MCP split", {HistType::kTH1F, {{1, 0.0, 1.0}}}); - registry.add("h_Matched_splitevent_counter", "Events into Matched split", {HistType::kTH1F, {{1, 0.0, 1.0}}}); - registry.add("hRandomValueDebug", "Random values for debugging;Random Value;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}); - - // DEBUG: Add counters for total events processed (before splitting) - registry.add("h_MCD_total_events", "Total MCD events processed", {HistType::kTH1F, {{1, 0.0, 1.0}}}); - registry.add("h_MCP_total_events", "Total MCP events processed", {HistType::kTH1F, {{1, 0.0, 1.0}}}); - registry.add("h_Matched_total_events", "Total Matched events processed", {HistType::kTH1F, {{1, 0.0, 1.0}}}); - registry.add("hMCCollisionIdDebug_MCP", "MC Collision Ids being processed", {HistType::kTH1F, {{100000, 0.0, 100000.0}}}); + trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast(trackSelections)); + eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventSelections)); + triggerMaskBits = jetderiveddatautilities::initialiseTriggerMaskBits(triggerMasks); + // triggerMaskBits = jetderiveddatautilities::initialiseTriggerMaskBits(jetderiveddatautilities::JTriggerMasks); + particleSelection = static_cast(particleSelections); + jetRadiiValues = (std::vector)jetRadii; -} -*/ -for (std::size_t iJetRadius = 0; iJetRadius < jetRadiiValues.size(); iJetRadius++) { - filledJetR.push_back(0.0); -} -auto jetRadiiBins = (std::vector)jetRadii; -if (jetRadiiBins.size() > 1) { - jetRadiiBins.push_back(jetRadiiBins[jetRadiiBins.size() - 1] + (std::abs(jetRadiiBins[jetRadiiBins.size() - 1] - jetRadiiBins[jetRadiiBins.size() - 2]))); -} else { - jetRadiiBins.push_back(jetRadiiBins[jetRadiiBins.size() - 1] + 0.1); -} + /* if (doMcClosure) { + // randGen.SetSeed(mcSplitSeed); + // randGen.SetSeed(static_cast(std::time(nullptr))); + // int seed = mcSplitSeed >= 0 ? mcSplitSeed : static_cast(std::time(nullptr)); + // randGen.SetSeed(seed); + // LOGF(info, "MC closure seed = %d", seed); -// Track QA histograms -if (doprocessDataTracks || doprocessMCTracks || doprocessTracksWeighted) { - registry.add("hCollisionsUnweighted", "event status; event status;entries", {HistType::kTH1F, {{12, 0., 12.0}}}); + int seed = mcSplitSeed >= 0 ? mcSplitSeed : static_cast(std::time(nullptr)); + randGen.SetSeed(seed); + LOGF(info, "MC closure splitting enabled with seed = %d, split fraction = %.2f", seed, static_cast(mcClosureSplitFrac)); - registry.add("h_track_pt", "track pT;#it{p}_{T,track} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h_track_eta", "track #eta;#eta_{track};entries", {HistType::kTH1F, {{100, -1., 1.}}}); - registry.add("h_track_phi", "track #varphi;#varphi_{track};entries", {HistType::kTH1F, {{160, 0., 7.}}}); - registry.add("h_track_energy", "track energy;Energy of tracks;entries", {HistType::kTH1F, {{400, 0., 400.}}}); - registry.add("h_track_energysum", "track energy sum;Sum of track energy per event;entries", {HistType::kTH1F, {{400, 0., 400.}}}); + registry.add("hSpliteventSelector", "Random MC Split Selector;Split Type;Entries",{HistType::kTH1F, {{3, 0.0, 3.0}}}); // 0=MCD, 1=MCP, 2=RM - // Cluster QA histograms - registry.add("h_clusterTime", "Time of cluster", HistType::kTH1F, {{500, -250, 250, "#it{t}_{cls} (ns)"}}); - registry.add("h_cluster_pt", "cluster pT;#it{p}_{T_cluster} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}); - registry.add("h_cluster_eta", "cluster #eta;#eta_{cluster};entries", {HistType::kTH1F, {{100, -1., 1.}}}); - registry.add("h_cluster_phi", "cluster #varphi;#varphi_{cluster};entries", {HistType::kTH1F, {{160, 0., 7.}}}); - registry.add("h_cluster_energy", "cluster energy;Energy of cluster;entries", {HistType::kTH1F, {{400, 0., 400.}}}); - registry.add("h_cluster_energysum", "cluster energy sum;Sum of cluster energy per event;entries", {HistType::kTH1F, {{400, 0., 400.}}}); + //individual processes' event counters for sanity checks + registry.add("h_MCD_splitevent_counter", "Events into MCD split", {HistType::kTH1F, {{1, 0.0, 1.0}}}); + registry.add("h_MCP_splitevent_counter", "Events into MCP split", {HistType::kTH1F, {{1, 0.0, 1.0}}}); + registry.add("h_Matched_splitevent_counter", "Events into Matched split", {HistType::kTH1F, {{1, 0.0, 1.0}}}); + registry.add("hRandomValueDebug", "Random values for debugging;Random Value;Entries", {HistType::kTH1F, {{100, 0.0, 1.0}}}); + + // DEBUG: Add counters for total events processed (before splitting) + registry.add("h_MCD_total_events", "Total MCD events processed", {HistType::kTH1F, {{1, 0.0, 1.0}}}); + registry.add("h_MCP_total_events", "Total MCP events processed", {HistType::kTH1F, {{1, 0.0, 1.0}}}); + registry.add("h_Matched_total_events", "Total Matched events processed", {HistType::kTH1F, {{1, 0.0, 1.0}}}); + registry.add("hMCCollisionIdDebug_MCP", "MC Collision Ids being processed", {HistType::kTH1F, {{100000, 0.0, 100000.0}}}); - if (doprocessTracksWeighted) { - registry.add("hCollisionsWeighted", "event status;event status;entries", {HistType::kTH1F, {{12, 0.0, 12.0}}}); } -} + */ + for (std::size_t iJetRadius = 0; iJetRadius < jetRadiiValues.size(); iJetRadius++) { + filledJetR.push_back(0.0); + } + auto jetRadiiBins = (std::vector)jetRadii; + if (jetRadiiBins.size() > 1) { + jetRadiiBins.push_back(jetRadiiBins[jetRadiiBins.size() - 1] + (std::abs(jetRadiiBins[jetRadiiBins.size() - 1] - jetRadiiBins[jetRadiiBins.size() - 2]))); + } else { + jetRadiiBins.push_back(jetRadiiBins[jetRadiiBins.size() - 1] + 0.1); + } -// Jet QA histograms -if (doprocessJetsData || doprocessJetsMCD || doprocessJetsMCDWeighted) { - - registry.add("hDetcollisionCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.}}}); - - registry.add("h_full_jet_pt", "#it{p}_{T,jet};#it{p}_{T_jet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h_full_jet_pt_pTHatcut", "#it{p}_{T,jet};#it{p}_{T_jet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h_full_jet_eta", "jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1., 1.}}}); - registry.add("h_full_jet_phi", "jet #varphi;#varphi_{jet};entries", {HistType::kTH1F, {{160, 0., 7.}}}); - registry.add("h_full_jet_clusterTime", "Time of cluster", HistType::kTH1F, {{500, -250, 250, "#it{t}_{cls} (ns)"}}); - registry.add("h2_full_jet_nef", "#it{p}_{T,jet} vs nef at Det Level; #it{p}_{T,jet} (GeV/#it{c});nef", {HistType::kTH2F, {{350, 0., 350.}, {105, 0., 1.05}}}); - registry.add("h2_full_jet_nef_rejected", "#it{p}_{T,jet} vs nef at Det Level for rejected events; #it{p}_{T,jet} (GeV/#it{c});nef", {HistType::kTH2F, {{350, 0., 350.}, {105, 0., 1.05}}}); - - registry.add("h_Detjet_ntracks", "#it{p}_{T,track};#it{p}_{T,track} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h2_full_jet_chargedconstituents", "Number of charged constituents at Det Level;#it{p}_{T,jet} (GeV/#it{c});N_{ch}", {HistType::kTH2F, {{350, 0., 350.}, {100, 0., 100.}}}); - registry.add("h2_full_jet_neutralconstituents", "Number of neutral constituents at Det Level;#it{p}_{T,jet} (GeV/#it{c});N_{ne}", {HistType::kTH2F, {{350, 0., 350.}, {100, 0., 100.}}}); - registry.add("h_full_jet_chargedconstituents_pt", "track pT;#it{p}^{T,jet}_{track} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h_full_jet_chargedconstituents_eta", "track #eta;#eta^{jet}_{track};entries", {HistType::kTH1F, {{100, -1., 1.}}}); - registry.add("h_full_jet_chargedconstituents_phi", "track #varphi;#varphi^{jet}_{track};entries", {HistType::kTH1F, {{160, 0., 7.}}}); - registry.add("h_full_jet_chargedconstituents_energy", "track energy;Energy of tracks;entries", {HistType::kTH1F, {{400, 0., 400.}}}); - registry.add("h_full_jet_chargedconstituents_energysum", "track energy sum;Sum of track energy per event;entries", {HistType::kTH1F, {{400, 0., 400.}}}); - registry.add("h_full_jet_neutralconstituents_pt", "cluster pT;#it{p}^{T,jet}_{cluster} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}); - registry.add("h_full_jet_neutralconstituents_eta", "cluster #eta;#eta^{jet}_{cluster};entries", {HistType::kTH1F, {{100, -1., 1.}}}); - registry.add("h_full_jet_neutralconstituents_phi", "cluster #varphi;#varphi^{jet}_{cluster};entries", {HistType::kTH1F, {{160, 0., 7.}}}); - registry.add("h_full_jet_neutralconstituents_energy", "cluster energy;Energy of cluster;entries", {HistType::kTH1F, {{400, 0., 400.}}}); - registry.add("h_full_jet_neutralconstituents_energysum", "cluster energy sum;Sum of cluster energy per event;entries", {HistType::kTH1F, {{400, 0., 400.}}}); - registry.add("h2_full_jettrack_pt", "#it{p}_{T,jet} vs #it{p}_{T,track}; #it{p}_{T,jet} (GeV/#it{c});#it{p}_{T,track} (GeV/#it{c})", {HistType::kTH2F, {{350, 0., 350.}, {200, 0., 200.}}}); - registry.add("h2_full_jettrack_eta", "jet #eta vs jet_track #eta; #eta_{jet};#eta_{track}", {HistType::kTH2F, {{100, -1., 1.}, {500, -5., 5.}}}); - registry.add("h2_full_jettrack_phi", "jet #varphi vs jet_track #varphi; #varphi_{jet}; #varphi_{track}", {HistType::kTH2F, {{160, 0., 7.}, {160, -1., 7.}}}); - - registry.add("h2_track_etaphi", "jet_track #eta vs jet_track #varphi; #eta_{track};#varphi_{track}", {HistType::kTH2F, {{500, -5., 5.}, {160, -1., 7.}}}); - registry.add("h2_jet_etaphi", "jet #eta vs jet #varphi; #eta_{jet};#varphi_{jet}", {HistType::kTH2F, {{100, -1., 1.}, {160, -1., 7.}}}); -} -if (doprocessJetsTriggeredData) { - registry.add("hDetTrigcollisionCounter", "event status;;entries", {HistType::kTH1F, {{14, 0.0, 14.}}}); -} -if (doprocessJetsMCP || doprocessJetsMCPWeighted) { - registry.add("hPartcollisionCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}}); - - registry.add("h_full_mcpjet_tablesize", "", {HistType::kTH1F, {{4, 0., 5.}}}); - registry.add("h_full_mcpjet_ntracks", "", {HistType::kTH1F, {{200, -0.5, 200.}}}); - registry.add("h_full_jet_pt_part", "jet pT;#it{p}_{T_jet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h_full_jet_eta_part", "jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1., 1.}}}); - registry.add("h_full_jet_phi_part", "jet #varphi;#varphi_{jet};entries", {HistType::kTH1F, {{160, 0., 7.}}}); - registry.add("h2_full_jet_nef_part", "#it{p}_{T,jet} vs nef at Part Level;#it{p}_{T,jet} (GeV/#it{c});nef", {HistType::kTH2F, {{350, 0., 350.}, {105, 0., 1.05}}}); - - registry.add("h_Partjet_ntracks", "#it{p}_{T,constituent};#it{p}_{T_constituent} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h2_full_jet_chargedconstituents_part", "Number of charged constituents at Part Level;#it{p}_{T,jet} (GeV/#it{c});N_{ch}", {HistType::kTH2F, {{350, 0., 350.}, {100, 0., 100.}}}); - registry.add("h2_full_jet_neutralconstituents_part", "Number of neutral constituents at Part Level;#it{p}_{T,jet} (GeV/#it{c});N_{ne}", {HistType::kTH2F, {{350, 0., 350.}, {100, 0., 100.}}}); - registry.add("h_full_jet_neutralconstituents_pt_part", "#it{p}_{T} of neutral constituents at Part Level;#it{p}_{T,ne} (GeV/#it{c}); entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h_full_jet_neutralconstituents_eta_part", "#eta of neutral constituents at Part Level;#eta_{ne};entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h_full_jet_neutralconstituents_phi_part", "#varphi of neutral constituents at Part Level;#varphi_{ne};entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h_full_jet_neutralconstituents_energy_part", "neutral constituents' energy;Energy of neutral constituents;entries", {HistType::kTH1F, {{400, 0., 400.}}}); - registry.add("h_full_jet_neutralconstituents_energysum_part", "neutral constituents' energy sum;Sum of neutral constituents' energy per event;entries", {HistType::kTH1F, {{400, 0., 400.}}}); - - registry.add("h2_jettrack_pt_part", "#it{p}_{T,jet} vs #it{p}_{T_track}; #it{p}_{T_jet} (GeV/#it{c});#it{p}_{T_track} (GeV/#it{c})", {HistType::kTH2F, {{350, 0., 350.}, {200, 0., 200.}}}); - registry.add("h2_jettrack_eta_part", "jet #eta vs jet_track #eta; #eta_{jet};#eta_{track}", {HistType::kTH2F, {{100, -1., 1.}, {500, -5., 5.}}}); - registry.add("h2_jettrack_phi_part", "jet #varphi vs jet_track #varphi; #varphi_{jet}; #varphi_{track}", {HistType::kTH2F, {{160, 0., 7.}, {160, -1., 7.}}}); - - registry.add("h2_track_etaphi_part", "jet_track #eta vs jet_track #varphi; #eta_{track};#varphi_{track}", {HistType::kTH2F, {{500, -5., 5.}, {160, -1., 7.}}}); - registry.add("h2_jet_etaphi_part", "jet #eta vs jet #varphi; #eta_{jet};#varphi_{jet}", {HistType::kTH2F, {{100, -1., 1.}, {160, -1., 7.}}}); - - // registry.add("h_NOmcpemcalcollisions", "event status;entries", {HistType::kTH1F, {{100, 0., 100.}}}); - // registry.add("h_mcpemcalcollisions", "event status;entries", {HistType::kTH1F, {{100, 0., 100.}}}); - registry.add("h2_full_mcpjetOutsideFiducial_pt", "MCP jet outside EMC Fiducial Acceptance #it{p}_{T,part};#it{p}_{T,part} (GeV/c); Ncounts", {HistType::kTH2F, {{350, 0., 350.}, {10000, 0., 10000.}}}); - registry.add("h_full_mcpjetOutside_eta_part", "MCP jet #eta outside EMC Fiducial Acceptance;#eta_{jet};entries", {HistType::kTH1F, {{100, -1., 1.}}}); - registry.add("h_full_mcpjetOutside_phi_part", "MCP jet #varphi outside EMC Fiducial Acceptance;#varphi_{jet};entries", {HistType::kTH1F, {{160, 0., 7.}}}); - registry.add("h2_full_mcpjetInsideFiducial_pt", "MCP jet #it{p}_{T,part} inside EMC Fiducial Acceptance;#it{p}_{T,part} (GeV/c); Ncounts", {HistType::kTH2F, {{350, 0., 350.}, {10000, 0., 10000.}}}); - registry.add("h_full_mcpjetInside_eta_part", "MCP jet #eta inside EMC Fiducial Acceptance;#eta_{jet};entries", {HistType::kTH1F, {{100, -1., 1.}}}); - registry.add("h_full_mcpjetInside_phi_part", "MCP jet #varphi inside EMC Fiducial Acceptance;#varphi_{jet};entries", {HistType::kTH1F, {{160, 0., 7.}}}); -} + // Track QA histograms + if (doprocessDataTracks || doprocessMCTracks || doprocessTracksWeighted) { + registry.add("hCollisionsUnweighted", "event status; event status;entries", {HistType::kTH1F, {{12, 0., 12.0}}}); -if (doprocessJetsMCPMCDMatched || doprocessJetsMCPMCDMatchedWeighted) { - registry.add("hMatchedcollisionCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}}); - - registry.add("h_full_matchedmcdjet_tablesize", "", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h_full_matchedmcpjet_tablesize", "", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h_full_matchedmcdjet_ntracks", "", {HistType::kTH1F, {{200, -0.5, 200.}}}); - registry.add("h_full_matchedmcpjet_ntracks", "", {HistType::kTH1F, {{200, -0.5, 200.}}}); - registry.add("h_full_matchedmcdjet_eta", "Matched MCD jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1., 1.}}}); - registry.add("h_full_matchedmcdjet_phi", "Matched MCD jet #varphi;#varphi_{jet};entries", {HistType::kTH1F, {{160, 0., 7.}}}); - registry.add("h_full_matchedmcpjet_eta", "Matched MCP jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1., 1.}}}); - registry.add("h_full_matchedmcpjet_phi", "Matched MCP jet #varphi;#varphi_{jet};entries", {HistType::kTH1F, {{160, 0., 7.}}}); - registry.add("h_full_jet_deltaR", "Distance between matched Det Jet and Part Jet; #Delta R; entries", {HistType::kTH1F, {{100, 0., 1.}}}); - - registry.add("h2_full_jet_energyscaleDet", "Jet Energy Scale (det); p_{T,det} (GeV/c); (p_{T,det} - p_{T,part})/p_{T,part}", {HistType::kTH2F, {{400, 0., 400.}, {200, -1., 1.}}}); - - registry.add("h2_matchedjet_etaphiDet", "Det jet #eta vs jet #varphi; #eta_{jet};#varphi_{jet}", {HistType::kTH2F, {{100, -1., 1.}, {160, -1., 7.}}}); - registry.add("h2_matchedjet_etaphiPart", "Part jet #eta vs jet #varphi; #eta_{jet};#varphi_{jet}", {HistType::kTH2F, {{100, -1., 1.}, {160, -1., 7.}}}); - registry.add("h2_matchedjet_deltaEtaCorr", "Correlation between Det Eta and Part Eta; #eta_{jet,det}; #eta_{jet,part}", {HistType::kTH2F, {{100, -1., 1.}, {100, -1., 1.}}}); - registry.add("h2_matchedjet_deltaPhiCorr", "Correlation between Det Phi and Part Phi; #varphi_{jet,det}; #varphi_{jet,part}", {HistType::kTH2F, {{160, 0., 7.}, {160, 0., 7.}}}); - - registry.add("h2_full_jet_energyscalePart", "Jet Energy Scale (part); p_{T,part} (GeV/c); (p_{T,det} - p_{T,part})/p_{T,part}", {HistType::kTH2F, {{400, 0., 400.}, {200, -1., 1.}}}); - registry.add("h3_full_jet_energyscalePart", "R dependence of Jet Energy Scale (Part); #it{R}_{jet};p_{T,det} (GeV/c); (p_{T,det} - p_{T,part})/p_{T,part}", {HistType::kTH3F, {{jetRadiiBins, ""}, {400, 0., 400.}, {200, -1., 1.}}}); - registry.add("h2_full_jet_etaresolutionPart", ";p_{T,part} (GeV/c); (#eta_{jet,det} - #eta_{jet,part})/#eta_{jet,part}", {HistType::kTH2F, {{400, 0., 400.}, {100, -1., 1.}}}); - registry.add("h2_full_jet_phiresolutionPart", ";p_{T,part} (GeV/c); (#varphi_{jet,det} - #varphi_{jet,part})/#varphi_{jet,part}", {HistType::kTH2F, {{400, 0., 400.}, {160, -1., 7.}}}); - registry.add("h2_full_jet_energyscaleChargedPart", "Jet Energy Scale (charged part); p_{T,part} (GeV/c); (p_{T,det} - p_{T,part})/p_{T,part}", {HistType::kTH2F, {{400, 0., 400.}, {200, -1., 1.}}}); - registry.add("h2_full_jet_energyscaleNeutralPart", "Jet Energy Scale (neutral part); p_{T,part} (GeV/c); (p_{T,det} - p_{T,part})/p_{T,part}", {HistType::kTH2F, {{400, 0., 400.}, {200, -1., 1.}}}); - registry.add("h2_full_jet_energyscaleChargedVsFullPart", "Jet Energy Scale (charged part, vs. full jet pt); p_{T,part} (GeV/c); (p_{T,det} - p_{T,part})/p_{T,part}", {HistType::kTH2F, {{400, 0., 400.}, {200, -1., 1.}}}); - registry.add("h2_full_jet_energyscaleNeutralVsFullPart", "Jet Energy Scale (neutral part, vs. full jet pt); p_{T,part} (GeV/c); (p_{T,det} - p_{T,part})/p_{T,part}", {HistType::kTH2F, {{400, 0., 400.}, {200, -1., 1.}}}); - registry.add("h2_full_fakemcdjets", "Fake MCD Jets; p_{T,det} (GeV/c); NCounts", {HistType::kTH2F, {{350, 0., 350.}, {100, 0., 100.}}}); - registry.add("h2FullfakeMcpJets", "Fake MCP Jets; p_{T,part} (GeV/c); NCounts", {HistType::kTH2F, {{350, 0., 350.}, {100, 0., 100.}}}); - registry.add("h2_full_matchedmcpjet_pt", "Matched MCP jet in EMC Fiducial Acceptance #it{p}_{T,part};#it{p}_{T,part} (GeV/c); Ncounts", {HistType::kTH2F, {{350, 0., 350.}, {10000, 0., 10000.}}}); - - // Response Matrix - registry.add("h_full_jet_ResponseMatrix", "Full Jets Response Matrix; p_{T,det} (GeV/c); p_{T,part} (GeV/c)", {HistType::kTH2F, {{500, 0., 500.}, {500, 0., 500.}}}); -} + registry.add("h_track_pt", "track pT;#it{p}_{T,track} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_track_eta", "track #eta;#eta_{track};entries", {HistType::kTH1F, {{100, -1., 1.}}}); + registry.add("h_track_phi", "track #varphi;#varphi_{track};entries", {HistType::kTH1F, {{160, 0., 7.}}}); + registry.add("h_track_energy", "track energy;Energy of tracks;entries", {HistType::kTH1F, {{400, 0., 400.}}}); + registry.add("h_track_energysum", "track energy sum;Sum of track energy per event;entries", {HistType::kTH1F, {{400, 0., 400.}}}); -if (doprocessMBCollisionsDATAWithMultiplicity || doprocessMBMCDCollisionsWithMultiplicity || doprocessMCDCollisionsWeightedWithMultiplicity) { - registry.add("hEventmultiplicityCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}}); - registry.add("h_FT0Mults_occupancy", "", {HistType::kTH1F, {{3500, 0., 3500.}}}); - - registry.add("h_all_fulljet_Njets", "Full Jet Multiplicity (per Event)", {HistType::kTH1F, {{20, 0., 20.}}}); - registry.add("h_Leading_full_jet_pt", "#it{p}_{T,leading jet};#it{p}_{T_leading jet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h2_full_jet_leadingJetPt_vs_counts", ";#it{p}_{T_leading jet} (GeV/#it{c}); Counts", {HistType::kTH2F, {{350, 0., 350.}, {20, 0., 20.}}}); - registry.add("h_SubLeading_full_jet_pt", "#it{p}_{T,leading jet};#it{p}_{T_leading jet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h2_full_jet_subLeadingJetPt_vs_counts", ";#it{p}_{T_leading jet} (GeV/#it{c}); Counts", {HistType::kTH2F, {{350, 0., 350.}, {20, 0., 20.}}}); - //Inside Jet Loop: - //CASE 1: - registry.add("h_all_fulljet_pt", "#it{p}_{T,fulljet};#it{p}_{T_fulljet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h_all_fulljet_Nch", ";N_{ch};", {HistType::kTH1F, {{50, 0., 50.}}}); - registry.add("h_all_fulljet_NEF", ";NEF;", {HistType::kTH1F, {{105, 0., 1.05}}}); - registry.add("h2_all_fulljet_jetpTDet_vs_FT0Mults", "; p_{T,det} (GeV/c); FT0M Multiplicity", {HistType::kTH2F, {{350, 0., 350.}, {3500, 0., 3500.}}}); - registry.add("h2_all_fulljet_jetpTDet_vs_Nch", ";#it{p}_{T_fulljet} (GeV/#it{c}); N_{ch}", {HistType::kTH2F, {{350, 0., 350.}, {50, 0., 50.}}}); - registry.add("h3_full_jet_jetpTDet_FT0Mults_nef", "; p_{T,det} (GeV/c); FT0M Multiplicity, nef", {HistType::kTH3F, {{350, 0., 350.}, {50, 0., 50.}, {105, 0.0, 1.05}}}); - //CASE 2: - registry.add("h_leading_fulljet_pt", "#it{p}_{T,Leading fulljet};#it{p}_{T_Leadingfulljet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h_leading_fulljet_Nch", ";N_{ch};", {HistType::kTH1F, {{50, 0., 50.}}}); - registry.add("h_leading_fulljet_NEF", ";NEF;", {HistType::kTH1F, {{105, 0., 1.05}}}); - registry.add("h2_leading_fulljet_jetpTDet_vs_FT0Mults", ";Leading p_{T,det} (GeV/c); FT0M Multiplicity", {HistType::kTH2F, {{350, 0., 350.}, {3500, 0., 3500.}}}); - registry.add("h2_leading_fulljet_jetpTDet_vs_Nch", ";#it{p}_{T_Leadingfulljet} (GeV/#it{c}); N_{ch}", {HistType::kTH2F, {{350, 0., 350.}, {50, 0., 50.}}}); - registry.add("h3_leading_fulljet_jetpTDet_FT0Mults_nef", "; Leading p_{T,det} (GeV/c); FT0M Multiplicity, nef", {HistType::kTH3F, {{350, 0., 350.}, {50, 0., 50.}, {105, 0.0, 1.05}}}); - //CASE 3: - registry.add("h_subleading_fulljet_pt", "#it{p}_{T,SubLeading fulljet};#it{p}_{T_SubLeadingfulljet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h_subleading_fulljet_Nch", ";N_{ch};", {HistType::kTH1F, {{50, 0., 50.}}}); - registry.add("h_subleading_fulljet_NEF", ";NEF;", {HistType::kTH1F, {{105, 0., 1.05}}}); - registry.add("h2_subleading_fulljet_jetpTDet_vs_FT0Mults", ";SubLeading p_{T,det} (GeV/c); FT0M Multiplicity", {HistType::kTH2F, {{350, 0., 350.}, {3500, 0., 3500.}}}); - registry.add("h2_subleading_fulljet_jetpTDet_vs_Nch", ";#it{p}_{T_SubLeadingfulljet} (GeV/#it{c}); N_{ch}", {HistType::kTH2F, {{350, 0., 350.}, {50, 0., 50.}}}); - registry.add("h3_subleading_fulljet_jetpTDet_FT0Mults_nef", "; SubLeading p_{T,det} (GeV/c); FT0M Multiplicity, nef", {HistType::kTH3F, {{350, 0., 350.}, {50, 0., 50.}, {105, 0.0, 1.05}}}); + // Cluster QA histograms + registry.add("h_clusterTime", "Time of cluster", HistType::kTH1F, {{500, -250, 250, "#it{t}_{cls} (ns)"}}); + registry.add("h_cluster_pt", "cluster pT;#it{p}_{T_cluster} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}); + registry.add("h_cluster_eta", "cluster #eta;#eta_{cluster};entries", {HistType::kTH1F, {{100, -1., 1.}}}); + registry.add("h_cluster_phi", "cluster #varphi;#varphi_{cluster};entries", {HistType::kTH1F, {{160, 0., 7.}}}); + registry.add("h_cluster_energy", "cluster energy;Energy of cluster;entries", {HistType::kTH1F, {{400, 0., 400.}}}); + registry.add("h_cluster_energysum", "cluster energy sum;Sum of cluster energy per event;entries", {HistType::kTH1F, {{400, 0., 400.}}}); -} + if (doprocessTracksWeighted) { + registry.add("hCollisionsWeighted", "event status;event status;entries", {HistType::kTH1F, {{12, 0.0, 12.0}}}); + } + } -if (doprocessMBMCPCollisionsWithMultiplicity || doprocessMBMCPCollisionsWeightedWithMultiplicity) { - registry.add("hPartEventmultiplicityCounter", "event status;event status;entries", {HistType::kTH1F, {{11, 0.0, 11.0}}}); - registry.add("hRecoMatchesPerMcCollision", "split vertices QA;;entries", {HistType::kTH1F, {{5, 0.0, 5.0}}}); - registry.add("hMCCollMatchedFT0Mult", "", {HistType::kTH1F, {{3500, 0., 3500.}}}); - registry.add("hMCCollMatchedFT0Cent", "", {HistType::kTH1F, {{105, 0., 105.}}}); - - registry.add("h_all_fulljet_Njets_part", "Full Jet Multiplicity (per Event)", {HistType::kTH1F, {{20, 0., 20.}}}); - registry.add("h_Leading_full_jet_pt_part", "#it{p}_{T,leading jet};#it{p}_{T_leading jet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h2_full_jet_leadingJetPt_vs_counts_part", ";#it{p}_{T_leading jet} (GeV/#it{c}); Counts", {HistType::kTH2F, {{350, 0., 350.}, {20, 0., 20.}}}); - registry.add("h_SubLeading_full_jet_pt_part", "#it{p}_{T,leading jet};#it{p}_{T_leading jet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h2_full_jet_subLeadingJetPt_vs_counts_part", ";#it{p}_{T_leading jet} (GeV/#it{c}); Counts", {HistType::kTH2F, {{350, 0., 350.}, {20, 0., 20.}}}); - - //Inside Jet Loop: - //CASE 1: - registry.add("h_all_fulljet_pt_part", "#it{p}_{T,fulljet};#it{p}_{T_fulljet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h_all_fulljet_Nch_part", ";N_{ch};", {HistType::kTH1F, {{50, 0., 50.}}}); - registry.add("h_all_fulljet_Nne_part", ";N_{ne};", {HistType::kTH1F, {{100, 0., 100.}}}); - registry.add("h_all_fulljet_NEF_part", ";NEF;", {HistType::kTH1F, {{105, 0., 1.05}}}); - registry.add("h2_all_fulljet_jetpT_vs_FT0Mults_part", "; p_{T,part} (GeV/c); FT0M Multiplicity", {HistType::kTH2F, {{350, 0., 350.}, {3500, 0., 3500.}}}); - registry.add("h2_all_fulljet_jetpT_vs_Nch_part", ";#it{p}_{T_fulljet} (GeV/#it{c}); N_{ch}", {HistType::kTH2F, {{350, 0., 350.}, {50, 0., 50.}}}); - registry.add("h3_full_jet_jetpT_FT0Mults_nef_part", "; p_{T,part} (GeV/c); FT0M Multiplicity, nef", {HistType::kTH3F, {{350, 0., 350.}, {50, 0., 50.}, {105, 0.0, 1.05}}}); - //CASE 2: - registry.add("h_leading_fulljet_pt_part", "#it{p}_{T,Leading fulljet};#it{p}_{T_Leadingfulljet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h_leading_fulljet_Nch_part", ";N_{ch};", {HistType::kTH1F, {{50, 0., 50.}}}); - registry.add("h_leading_fulljet_Nne_part", ";N_{ne};", {HistType::kTH1F, {{100, 0., 100.}}}); - registry.add("h_leading_fulljet_NEF_part", ";NEF;", {HistType::kTH1F, {{105, 0., 1.05}}}); - registry.add("h2_leading_fulljet_jetpT_vs_FT0Mults_part", ";Leading p_{T,part} (GeV/c); FT0M Multiplicity", {HistType::kTH2F, {{350, 0., 350.}, {3500, 0., 3500.}}}); - registry.add("h2_leading_fulljet_jetpT_vs_Nch_part", ";#it{p}_{T_Leadingfulljet} (GeV/#it{c}); N_{ch}", {HistType::kTH2F, {{350, 0., 350.}, {50, 0., 50.}}}); - registry.add("h3_leading_fulljet_jetpT_FT0Mults_nef_part", "; Leading p_{T,part} (GeV/c); FT0M Multiplicity, nef", {HistType::kTH3F, {{350, 0., 350.}, {50, 0., 50.}, {105, 0.0, 1.05}}}); - //CASE 3: - registry.add("h_subleading_fulljet_pt_part", "#it{p}_{T,SubLeading fulljet};#it{p}_{T_SubLeadingfulljet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); - registry.add("h_subleading_fulljet_Nch_part", ";N_{ch};", {HistType::kTH1F, {{50, 0., 50.}}}); - registry.add("h_subleading_fulljet_Nne_part", ";N_{ne};", {HistType::kTH1F, {{100, 0., 100.}}}); - registry.add("h_subleading_fulljet_NEF_part", ";NEF;", {HistType::kTH1F, {{105, 0., 1.05}}}); - registry.add("h2_subleading_fulljet_jetpT_vs_FT0Mults_part", ";SubLeading p_{T,part} (GeV/c); FT0M Multiplicity", {HistType::kTH2F, {{350, 0., 350.}, {3500, 0., 3500.}}}); - registry.add("h2_subleading_fulljet_jetpT_vs_Nch_part", ";#it{p}_{T_SubLeadingfulljet} (GeV/#it{c}); N_{ch}", {HistType::kTH2F, {{350, 0., 350.}, {50, 0., 50.}}}); - registry.add("h3_subleading_fulljet_jetpT_FT0Mults_nef_part", "; SubLeading p_{T,part} (GeV/c); FT0M Multiplicity, nef", {HistType::kTH3F, {{350, 0., 350.}, {50, 0., 50.}, {105, 0.0, 1.05}}}); + // Jet QA histograms + if (doprocessJetsData || doprocessJetsMCD || doprocessJetsMCDWeighted) { -} + registry.add("hDetcollisionCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.}}}); + + registry.add("h_full_jet_pt", "#it{p}_{T,jet};#it{p}_{T_jet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_full_jet_pt_pTHatcut", "#it{p}_{T,jet};#it{p}_{T_jet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_full_jet_eta", "jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1., 1.}}}); + registry.add("h_full_jet_phi", "jet #varphi;#varphi_{jet};entries", {HistType::kTH1F, {{160, 0., 7.}}}); + registry.add("h_full_jet_clusterTime", "Time of cluster", HistType::kTH1F, {{500, -250, 250, "#it{t}_{cls} (ns)"}}); + registry.add("h2_full_jet_nef", "#it{p}_{T,jet} vs nef at Det Level; #it{p}_{T,jet} (GeV/#it{c});nef", {HistType::kTH2F, {{350, 0., 350.}, {105, 0., 1.05}}}); + registry.add("h2_full_jet_nef_rejected", "#it{p}_{T,jet} vs nef at Det Level for rejected events; #it{p}_{T,jet} (GeV/#it{c});nef", {HistType::kTH2F, {{350, 0., 350.}, {105, 0., 1.05}}}); + + registry.add("h_Detjet_ntracks", "#it{p}_{T,track};#it{p}_{T,track} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h2_full_jet_chargedconstituents", "Number of charged constituents at Det Level;#it{p}_{T,jet} (GeV/#it{c});N_{ch}", {HistType::kTH2F, {{350, 0., 350.}, {100, 0., 100.}}}); + registry.add("h2_full_jet_neutralconstituents", "Number of neutral constituents at Det Level;#it{p}_{T,jet} (GeV/#it{c});N_{ne}", {HistType::kTH2F, {{350, 0., 350.}, {100, 0., 100.}}}); + registry.add("h_full_jet_chargedconstituents_pt", "track pT;#it{p}^{T,jet}_{track} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_full_jet_chargedconstituents_eta", "track #eta;#eta^{jet}_{track};entries", {HistType::kTH1F, {{100, -1., 1.}}}); + registry.add("h_full_jet_chargedconstituents_phi", "track #varphi;#varphi^{jet}_{track};entries", {HistType::kTH1F, {{160, 0., 7.}}}); + registry.add("h_full_jet_chargedconstituents_energy", "track energy;Energy of tracks;entries", {HistType::kTH1F, {{400, 0., 400.}}}); + registry.add("h_full_jet_chargedconstituents_energysum", "track energy sum;Sum of track energy per event;entries", {HistType::kTH1F, {{400, 0., 400.}}}); + registry.add("h_full_jet_neutralconstituents_pt", "cluster pT;#it{p}^{T,jet}_{cluster} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}); + registry.add("h_full_jet_neutralconstituents_eta", "cluster #eta;#eta^{jet}_{cluster};entries", {HistType::kTH1F, {{100, -1., 1.}}}); + registry.add("h_full_jet_neutralconstituents_phi", "cluster #varphi;#varphi^{jet}_{cluster};entries", {HistType::kTH1F, {{160, 0., 7.}}}); + registry.add("h_full_jet_neutralconstituents_energy", "cluster energy;Energy of cluster;entries", {HistType::kTH1F, {{400, 0., 400.}}}); + registry.add("h_full_jet_neutralconstituents_energysum", "cluster energy sum;Sum of cluster energy per event;entries", {HistType::kTH1F, {{400, 0., 400.}}}); + registry.add("h2_full_jettrack_pt", "#it{p}_{T,jet} vs #it{p}_{T,track}; #it{p}_{T,jet} (GeV/#it{c});#it{p}_{T,track} (GeV/#it{c})", {HistType::kTH2F, {{350, 0., 350.}, {200, 0., 200.}}}); + registry.add("h2_full_jettrack_eta", "jet #eta vs jet_track #eta; #eta_{jet};#eta_{track}", {HistType::kTH2F, {{100, -1., 1.}, {500, -5., 5.}}}); + registry.add("h2_full_jettrack_phi", "jet #varphi vs jet_track #varphi; #varphi_{jet}; #varphi_{track}", {HistType::kTH2F, {{160, 0., 7.}, {160, -1., 7.}}}); + + registry.add("h2_track_etaphi", "jet_track #eta vs jet_track #varphi; #eta_{track};#varphi_{track}", {HistType::kTH2F, {{500, -5., 5.}, {160, -1., 7.}}}); + registry.add("h2_jet_etaphi", "jet #eta vs jet #varphi; #eta_{jet};#varphi_{jet}", {HistType::kTH2F, {{100, -1., 1.}, {160, -1., 7.}}}); + } + if (doprocessJetsTriggeredData) { + registry.add("hDetTrigcollisionCounter", "event status;;entries", {HistType::kTH1F, {{14, 0.0, 14.}}}); + } + if (doprocessJetsMCP || doprocessJetsMCPWeighted) { + registry.add("hPartcollisionCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}}); + + registry.add("h_full_mcpjet_tablesize", "", {HistType::kTH1F, {{4, 0., 5.}}}); + registry.add("h_full_mcpjet_ntracks", "", {HistType::kTH1F, {{200, -0.5, 200.}}}); + registry.add("h_full_jet_pt_part", "jet pT;#it{p}_{T_jet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_full_jet_eta_part", "jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1., 1.}}}); + registry.add("h_full_jet_phi_part", "jet #varphi;#varphi_{jet};entries", {HistType::kTH1F, {{160, 0., 7.}}}); + registry.add("h2_full_jet_nef_part", "#it{p}_{T,jet} vs nef at Part Level;#it{p}_{T,jet} (GeV/#it{c});nef", {HistType::kTH2F, {{350, 0., 350.}, {105, 0., 1.05}}}); + + registry.add("h_Partjet_ntracks", "#it{p}_{T,constituent};#it{p}_{T_constituent} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h2_full_jet_chargedconstituents_part", "Number of charged constituents at Part Level;#it{p}_{T,jet} (GeV/#it{c});N_{ch}", {HistType::kTH2F, {{350, 0., 350.}, {100, 0., 100.}}}); + registry.add("h2_full_jet_neutralconstituents_part", "Number of neutral constituents at Part Level;#it{p}_{T,jet} (GeV/#it{c});N_{ne}", {HistType::kTH2F, {{350, 0., 350.}, {100, 0., 100.}}}); + registry.add("h_full_jet_neutralconstituents_pt_part", "#it{p}_{T} of neutral constituents at Part Level;#it{p}_{T,ne} (GeV/#it{c}); entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_full_jet_neutralconstituents_eta_part", "#eta of neutral constituents at Part Level;#eta_{ne};entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_full_jet_neutralconstituents_phi_part", "#varphi of neutral constituents at Part Level;#varphi_{ne};entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_full_jet_neutralconstituents_energy_part", "neutral constituents' energy;Energy of neutral constituents;entries", {HistType::kTH1F, {{400, 0., 400.}}}); + registry.add("h_full_jet_neutralconstituents_energysum_part", "neutral constituents' energy sum;Sum of neutral constituents' energy per event;entries", {HistType::kTH1F, {{400, 0., 400.}}}); + + registry.add("h2_jettrack_pt_part", "#it{p}_{T,jet} vs #it{p}_{T_track}; #it{p}_{T_jet} (GeV/#it{c});#it{p}_{T_track} (GeV/#it{c})", {HistType::kTH2F, {{350, 0., 350.}, {200, 0., 200.}}}); + registry.add("h2_jettrack_eta_part", "jet #eta vs jet_track #eta; #eta_{jet};#eta_{track}", {HistType::kTH2F, {{100, -1., 1.}, {500, -5., 5.}}}); + registry.add("h2_jettrack_phi_part", "jet #varphi vs jet_track #varphi; #varphi_{jet}; #varphi_{track}", {HistType::kTH2F, {{160, 0., 7.}, {160, -1., 7.}}}); + + registry.add("h2_track_etaphi_part", "jet_track #eta vs jet_track #varphi; #eta_{track};#varphi_{track}", {HistType::kTH2F, {{500, -5., 5.}, {160, -1., 7.}}}); + registry.add("h2_jet_etaphi_part", "jet #eta vs jet #varphi; #eta_{jet};#varphi_{jet}", {HistType::kTH2F, {{100, -1., 1.}, {160, -1., 7.}}}); + + // registry.add("h_NOmcpemcalcollisions", "event status;entries", {HistType::kTH1F, {{100, 0., 100.}}}); + // registry.add("h_mcpemcalcollisions", "event status;entries", {HistType::kTH1F, {{100, 0., 100.}}}); + registry.add("h2_full_mcpjetOutsideFiducial_pt", "MCP jet outside EMC Fiducial Acceptance #it{p}_{T,part};#it{p}_{T,part} (GeV/c); Ncounts", {HistType::kTH2F, {{350, 0., 350.}, {10000, 0., 10000.}}}); + registry.add("h_full_mcpjetOutside_eta_part", "MCP jet #eta outside EMC Fiducial Acceptance;#eta_{jet};entries", {HistType::kTH1F, {{100, -1., 1.}}}); + registry.add("h_full_mcpjetOutside_phi_part", "MCP jet #varphi outside EMC Fiducial Acceptance;#varphi_{jet};entries", {HistType::kTH1F, {{160, 0., 7.}}}); + registry.add("h2_full_mcpjetInsideFiducial_pt", "MCP jet #it{p}_{T,part} inside EMC Fiducial Acceptance;#it{p}_{T,part} (GeV/c); Ncounts", {HistType::kTH2F, {{350, 0., 350.}, {10000, 0., 10000.}}}); + registry.add("h_full_mcpjetInside_eta_part", "MCP jet #eta inside EMC Fiducial Acceptance;#eta_{jet};entries", {HistType::kTH1F, {{100, -1., 1.}}}); + registry.add("h_full_mcpjetInside_phi_part", "MCP jet #varphi inside EMC Fiducial Acceptance;#varphi_{jet};entries", {HistType::kTH1F, {{160, 0., 7.}}}); + } -// Label the histograms -labelCollisionHistograms(registry); -// labelMCSplitHistogram(registry); -} // init + if (doprocessJetsMCPMCDMatched || doprocessJetsMCPMCDMatchedWeighted) { + registry.add("hMatchedcollisionCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}}); + + registry.add("h_full_matchedmcdjet_tablesize", "", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_full_matchedmcpjet_tablesize", "", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_full_matchedmcdjet_ntracks", "", {HistType::kTH1F, {{200, -0.5, 200.}}}); + registry.add("h_full_matchedmcpjet_ntracks", "", {HistType::kTH1F, {{200, -0.5, 200.}}}); + registry.add("h_full_matchedmcdjet_eta", "Matched MCD jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1., 1.}}}); + registry.add("h_full_matchedmcdjet_phi", "Matched MCD jet #varphi;#varphi_{jet};entries", {HistType::kTH1F, {{160, 0., 7.}}}); + registry.add("h_full_matchedmcpjet_eta", "Matched MCP jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1., 1.}}}); + registry.add("h_full_matchedmcpjet_phi", "Matched MCP jet #varphi;#varphi_{jet};entries", {HistType::kTH1F, {{160, 0., 7.}}}); + registry.add("h_full_jet_deltaR", "Distance between matched Det Jet and Part Jet; #Delta R; entries", {HistType::kTH1F, {{100, 0., 1.}}}); + + registry.add("h2_full_jet_energyscaleDet", "Jet Energy Scale (det); p_{T,det} (GeV/c); (p_{T,det} - p_{T,part})/p_{T,part}", {HistType::kTH2F, {{400, 0., 400.}, {200, -1., 1.}}}); + + registry.add("h2_matchedjet_etaphiDet", "Det jet #eta vs jet #varphi; #eta_{jet};#varphi_{jet}", {HistType::kTH2F, {{100, -1., 1.}, {160, -1., 7.}}}); + registry.add("h2_matchedjet_etaphiPart", "Part jet #eta vs jet #varphi; #eta_{jet};#varphi_{jet}", {HistType::kTH2F, {{100, -1., 1.}, {160, -1., 7.}}}); + registry.add("h2_matchedjet_deltaEtaCorr", "Correlation between Det Eta and Part Eta; #eta_{jet,det}; #eta_{jet,part}", {HistType::kTH2F, {{100, -1., 1.}, {100, -1., 1.}}}); + registry.add("h2_matchedjet_deltaPhiCorr", "Correlation between Det Phi and Part Phi; #varphi_{jet,det}; #varphi_{jet,part}", {HistType::kTH2F, {{160, 0., 7.}, {160, 0., 7.}}}); + + registry.add("h2_full_jet_energyscalePart", "Jet Energy Scale (part); p_{T,part} (GeV/c); (p_{T,det} - p_{T,part})/p_{T,part}", {HistType::kTH2F, {{400, 0., 400.}, {200, -1., 1.}}}); + registry.add("h3_full_jet_energyscalePart", "R dependence of Jet Energy Scale (Part); #it{R}_{jet};p_{T,det} (GeV/c); (p_{T,det} - p_{T,part})/p_{T,part}", {HistType::kTH3F, {{jetRadiiBins, ""}, {400, 0., 400.}, {200, -1., 1.}}}); + registry.add("h2_full_jet_etaresolutionPart", ";p_{T,part} (GeV/c); (#eta_{jet,det} - #eta_{jet,part})/#eta_{jet,part}", {HistType::kTH2F, {{400, 0., 400.}, {100, -1., 1.}}}); + registry.add("h2_full_jet_phiresolutionPart", ";p_{T,part} (GeV/c); (#varphi_{jet,det} - #varphi_{jet,part})/#varphi_{jet,part}", {HistType::kTH2F, {{400, 0., 400.}, {160, -1., 7.}}}); + registry.add("h2_full_jet_energyscaleChargedPart", "Jet Energy Scale (charged part); p_{T,part} (GeV/c); (p_{T,det} - p_{T,part})/p_{T,part}", {HistType::kTH2F, {{400, 0., 400.}, {200, -1., 1.}}}); + registry.add("h2_full_jet_energyscaleNeutralPart", "Jet Energy Scale (neutral part); p_{T,part} (GeV/c); (p_{T,det} - p_{T,part})/p_{T,part}", {HistType::kTH2F, {{400, 0., 400.}, {200, -1., 1.}}}); + registry.add("h2_full_jet_energyscaleChargedVsFullPart", "Jet Energy Scale (charged part, vs. full jet pt); p_{T,part} (GeV/c); (p_{T,det} - p_{T,part})/p_{T,part}", {HistType::kTH2F, {{400, 0., 400.}, {200, -1., 1.}}}); + registry.add("h2_full_jet_energyscaleNeutralVsFullPart", "Jet Energy Scale (neutral part, vs. full jet pt); p_{T,part} (GeV/c); (p_{T,det} - p_{T,part})/p_{T,part}", {HistType::kTH2F, {{400, 0., 400.}, {200, -1., 1.}}}); + registry.add("h2_full_fakemcdjets", "Fake MCD Jets; p_{T,det} (GeV/c); NCounts", {HistType::kTH2F, {{350, 0., 350.}, {100, 0., 100.}}}); + registry.add("h2FullfakeMcpJets", "Fake MCP Jets; p_{T,part} (GeV/c); NCounts", {HistType::kTH2F, {{350, 0., 350.}, {100, 0., 100.}}}); + registry.add("h2_full_matchedmcpjet_pt", "Matched MCP jet in EMC Fiducial Acceptance #it{p}_{T,part};#it{p}_{T,part} (GeV/c); Ncounts", {HistType::kTH2F, {{350, 0., 350.}, {10000, 0., 10000.}}}); -// Initialize CCDB access and histogram registry for Zorro processing -template -void initCCDB(const BCType& bc) -{ - if (doSoftwareTriggerSelection) { - zorro.initCCDB(ccdb.service, bc.runNumber(), bc.timestamp(), triggerMasks.value); - zorro.populateHistRegistry(registry, bc.runNumber()); + // Response Matrix + registry.add("h_full_jet_ResponseMatrix", "Full Jets Response Matrix; p_{T,det} (GeV/c); p_{T,part} (GeV/c)", {HistType::kTH2F, {{500, 0., 500.}, {500, 0., 500.}}}); + } + + if (doprocessMBCollisionsDATAWithMultiplicity || doprocessMBMCDCollisionsWithMultiplicity || doprocessMCDCollisionsWeightedWithMultiplicity) { + registry.add("hEventmultiplicityCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}}); + registry.add("h_FT0Mults_occupancy", "", {HistType::kTH1F, {{3500, 0., 3500.}}}); + + registry.add("h_all_fulljet_Njets", "Full Jet Multiplicity (per Event)", {HistType::kTH1F, {{20, 0., 20.}}}); + registry.add("h_Leading_full_jet_pt", "#it{p}_{T,leading jet};#it{p}_{T_leading jet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h2_full_jet_leadingJetPt_vs_counts", ";#it{p}_{T_leading jet} (GeV/#it{c}); Counts", {HistType::kTH2F, {{350, 0., 350.}, {20, 0., 20.}}}); + registry.add("h_SubLeading_full_jet_pt", "#it{p}_{T,leading jet};#it{p}_{T_leading jet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h2_full_jet_subLeadingJetPt_vs_counts", ";#it{p}_{T_leading jet} (GeV/#it{c}); Counts", {HistType::kTH2F, {{350, 0., 350.}, {20, 0., 20.}}}); + // Inside Jet Loop: + // CASE 1: + registry.add("h_all_fulljet_pt", "#it{p}_{T,fulljet};#it{p}_{T_fulljet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_all_fulljet_Nch", ";N_{ch};", {HistType::kTH1F, {{50, 0., 50.}}}); + registry.add("h_all_fulljet_NEF", ";NEF;", {HistType::kTH1F, {{105, 0., 1.05}}}); + registry.add("h2_all_fulljet_jetpTDet_vs_FT0Mults", "; p_{T,det} (GeV/c); FT0M Multiplicity", {HistType::kTH2F, {{350, 0., 350.}, {3500, 0., 3500.}}}); + registry.add("h2_all_fulljet_jetpTDet_vs_Nch", ";#it{p}_{T_fulljet} (GeV/#it{c}); N_{ch}", {HistType::kTH2F, {{350, 0., 350.}, {50, 0., 50.}}}); + registry.add("h3_full_jet_jetpTDet_FT0Mults_nef", "; p_{T,det} (GeV/c); FT0M Multiplicity, nef", {HistType::kTH3F, {{350, 0., 350.}, {50, 0., 50.}, {105, 0.0, 1.05}}}); + // CASE 2: + registry.add("h_leading_fulljet_pt", "#it{p}_{T,Leading fulljet};#it{p}_{T_Leadingfulljet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_leading_fulljet_Nch", ";N_{ch};", {HistType::kTH1F, {{50, 0., 50.}}}); + registry.add("h_leading_fulljet_NEF", ";NEF;", {HistType::kTH1F, {{105, 0., 1.05}}}); + registry.add("h2_leading_fulljet_jetpTDet_vs_FT0Mults", ";Leading p_{T,det} (GeV/c); FT0M Multiplicity", {HistType::kTH2F, {{350, 0., 350.}, {3500, 0., 3500.}}}); + registry.add("h2_leading_fulljet_jetpTDet_vs_Nch", ";#it{p}_{T_Leadingfulljet} (GeV/#it{c}); N_{ch}", {HistType::kTH2F, {{350, 0., 350.}, {50, 0., 50.}}}); + registry.add("h3_leading_fulljet_jetpTDet_FT0Mults_nef", "; Leading p_{T,det} (GeV/c); FT0M Multiplicity, nef", {HistType::kTH3F, {{350, 0., 350.}, {50, 0., 50.}, {105, 0.0, 1.05}}}); + // CASE 3: + registry.add("h_subleading_fulljet_pt", "#it{p}_{T,SubLeading fulljet};#it{p}_{T_SubLeadingfulljet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_subleading_fulljet_Nch", ";N_{ch};", {HistType::kTH1F, {{50, 0., 50.}}}); + registry.add("h_subleading_fulljet_NEF", ";NEF;", {HistType::kTH1F, {{105, 0., 1.05}}}); + registry.add("h2_subleading_fulljet_jetpTDet_vs_FT0Mults", ";SubLeading p_{T,det} (GeV/c); FT0M Multiplicity", {HistType::kTH2F, {{350, 0., 350.}, {3500, 0., 3500.}}}); + registry.add("h2_subleading_fulljet_jetpTDet_vs_Nch", ";#it{p}_{T_SubLeadingfulljet} (GeV/#it{c}); N_{ch}", {HistType::kTH2F, {{350, 0., 350.}, {50, 0., 50.}}}); + registry.add("h3_subleading_fulljet_jetpTDet_FT0Mults_nef", "; SubLeading p_{T,det} (GeV/c); FT0M Multiplicity, nef", {HistType::kTH3F, {{350, 0., 350.}, {50, 0., 50.}, {105, 0.0, 1.05}}}); + } + + if (doprocessMBMCPCollisionsWithMultiplicity || doprocessMBMCPCollisionsWeightedWithMultiplicity) { + registry.add("hPartEventmultiplicityCounter", "event status;event status;entries", {HistType::kTH1F, {{11, 0.0, 11.0}}}); + registry.add("hRecoMatchesPerMcCollision", "split vertices QA;;entries", {HistType::kTH1F, {{5, 0.0, 5.0}}}); + registry.add("hMCCollMatchedFT0Mult", "", {HistType::kTH1F, {{3500, 0., 3500.}}}); + registry.add("hMCCollMatchedFT0Cent", "", {HistType::kTH1F, {{105, 0., 105.}}}); + + registry.add("h_all_fulljet_Njets_part", "Full Jet Multiplicity (per Event)", {HistType::kTH1F, {{20, 0., 20.}}}); + registry.add("h_Leading_full_jet_pt_part", "#it{p}_{T,leading jet};#it{p}_{T_leading jet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h2_full_jet_leadingJetPt_vs_counts_part", ";#it{p}_{T_leading jet} (GeV/#it{c}); Counts", {HistType::kTH2F, {{350, 0., 350.}, {20, 0., 20.}}}); + registry.add("h_SubLeading_full_jet_pt_part", "#it{p}_{T,leading jet};#it{p}_{T_leading jet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h2_full_jet_subLeadingJetPt_vs_counts_part", ";#it{p}_{T_leading jet} (GeV/#it{c}); Counts", {HistType::kTH2F, {{350, 0., 350.}, {20, 0., 20.}}}); + + // Inside Jet Loop: + // CASE 1: + registry.add("h_all_fulljet_pt_part", "#it{p}_{T,fulljet};#it{p}_{T_fulljet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_all_fulljet_Nch_part", ";N_{ch};", {HistType::kTH1F, {{50, 0., 50.}}}); + registry.add("h_all_fulljet_Nne_part", ";N_{ne};", {HistType::kTH1F, {{100, 0., 100.}}}); + registry.add("h_all_fulljet_NEF_part", ";NEF;", {HistType::kTH1F, {{105, 0., 1.05}}}); + registry.add("h2_all_fulljet_jetpT_vs_FT0Mults_part", "; p_{T,part} (GeV/c); FT0M Multiplicity", {HistType::kTH2F, {{350, 0., 350.}, {3500, 0., 3500.}}}); + registry.add("h2_all_fulljet_jetpT_vs_Nch_part", ";#it{p}_{T_fulljet} (GeV/#it{c}); N_{ch}", {HistType::kTH2F, {{350, 0., 350.}, {50, 0., 50.}}}); + registry.add("h3_full_jet_jetpT_FT0Mults_nef_part", "; p_{T,part} (GeV/c); FT0M Multiplicity, nef", {HistType::kTH3F, {{350, 0., 350.}, {50, 0., 50.}, {105, 0.0, 1.05}}}); + // CASE 2: + registry.add("h_leading_fulljet_pt_part", "#it{p}_{T,Leading fulljet};#it{p}_{T_Leadingfulljet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_leading_fulljet_Nch_part", ";N_{ch};", {HistType::kTH1F, {{50, 0., 50.}}}); + registry.add("h_leading_fulljet_Nne_part", ";N_{ne};", {HistType::kTH1F, {{100, 0., 100.}}}); + registry.add("h_leading_fulljet_NEF_part", ";NEF;", {HistType::kTH1F, {{105, 0., 1.05}}}); + registry.add("h2_leading_fulljet_jetpT_vs_FT0Mults_part", ";Leading p_{T,part} (GeV/c); FT0M Multiplicity", {HistType::kTH2F, {{350, 0., 350.}, {3500, 0., 3500.}}}); + registry.add("h2_leading_fulljet_jetpT_vs_Nch_part", ";#it{p}_{T_Leadingfulljet} (GeV/#it{c}); N_{ch}", {HistType::kTH2F, {{350, 0., 350.}, {50, 0., 50.}}}); + registry.add("h3_leading_fulljet_jetpT_FT0Mults_nef_part", "; Leading p_{T,part} (GeV/c); FT0M Multiplicity, nef", {HistType::kTH3F, {{350, 0., 350.}, {50, 0., 50.}, {105, 0.0, 1.05}}}); + // CASE 3: + registry.add("h_subleading_fulljet_pt_part", "#it{p}_{T,SubLeading fulljet};#it{p}_{T_SubLeadingfulljet} (GeV/#it{c});entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_subleading_fulljet_Nch_part", ";N_{ch};", {HistType::kTH1F, {{50, 0., 50.}}}); + registry.add("h_subleading_fulljet_Nne_part", ";N_{ne};", {HistType::kTH1F, {{100, 0., 100.}}}); + registry.add("h_subleading_fulljet_NEF_part", ";NEF;", {HistType::kTH1F, {{105, 0., 1.05}}}); + registry.add("h2_subleading_fulljet_jetpT_vs_FT0Mults_part", ";SubLeading p_{T,part} (GeV/c); FT0M Multiplicity", {HistType::kTH2F, {{350, 0., 350.}, {3500, 0., 3500.}}}); + registry.add("h2_subleading_fulljet_jetpT_vs_Nch_part", ";#it{p}_{T_SubLeadingfulljet} (GeV/#it{c}); N_{ch}", {HistType::kTH2F, {{350, 0., 350.}, {50, 0., 50.}}}); + registry.add("h3_subleading_fulljet_jetpT_FT0Mults_nef_part", "; SubLeading p_{T,part} (GeV/c); FT0M Multiplicity, nef", {HistType::kTH3F, {{350, 0., 350.}, {50, 0., 50.}, {105, 0.0, 1.05}}}); + } + + // Label the histograms + labelCollisionHistograms(registry); + // labelMCSplitHistogram(registry); + } // init + + // Initialize CCDB access and histogram registry for Zorro processing + template + void initCCDB(const BCType& bc) + { + if (doSoftwareTriggerSelection) { + zorro.initCCDB(ccdb.service, bc.runNumber(), bc.timestamp(), triggerMasks.value); + zorro.populateHistRegistry(registry, bc.runNumber()); + } } -} -// Get or generate random value for a specific MC collision -/* float getMCCollisionRandomValue(int64_t mcCollisionId) { -if (!doMcClosure) return 0.0f; + // Get or generate random value for a specific MC collision + /* float getMCCollisionRandomValue(int64_t mcCollisionId) { + if (!doMcClosure) return 0.0f; -// Check if I already have a random value for this MC collision -auto it = mcCollisionRandomValues.find(mcCollisionId); -if (it != mcCollisionRandomValues.end()) { -LOGF(debug, "Using cached random value %.4f for MC collision %lld", it->second, mcCollisionId); -return it->second; -} + // Check if I already have a random value for this MC collision + auto it = mcCollisionRandomValues.find(mcCollisionId); + if (it != mcCollisionRandomValues.end()) { + LOGF(debug, "Using cached random value %.4f for MC collision %lld", it->second, mcCollisionId); + return it->second; + } -// Generate new random value for this MC collision -float randomVal = randGen.Uniform(0.0, 1.0); -mcCollisionRandomValues[mcCollisionId] = randomVal; + // Generate new random value for this MC collision + float randomVal = randGen.Uniform(0.0, 1.0); + mcCollisionRandomValues[mcCollisionId] = randomVal; -// Debug histogram -registry.fill(HIST("hRandomValueDebug"), randomVal); + // Debug histogram + registry.fill(HIST("hRandomValueDebug"), randomVal); -LOGF(info, "Generated NEW random value %.4f for MC collision %lld", randomVal, mcCollisionId); -return randomVal; -} -*/ -using EMCCollisionsData = o2::soa::Join; // JetCollisions with EMCAL Collision Labels -using EMCCollisionsTriggeredData = o2::soa::Join; + LOGF(info, "Generated NEW random value %.4f for MC collision %lld", randomVal, mcCollisionId); + return randomVal; + } + */ + using EMCCollisionsData = o2::soa::Join; // JetCollisions with EMCAL Collision Labels + using EMCCollisionsTriggeredData = o2::soa::Join; -using EMCCollisionsMCD = o2::soa::Join; // where, JetCollisionsMCD = JetCollisions+JMcCollisionLbs + using EMCCollisionsMCD = o2::soa::Join; // where, JetCollisionsMCD = JetCollisions+JMcCollisionLbs -using FullJetTableDataJoined = soa::Join; -using JetTableMCDJoined = soa::Join; -using JetTableMCDWeightedJoined = soa::Join; -using JetTableMCPJoined = soa::Join; -using JetTableMCPWeightedJoined = soa::Join; + using FullJetTableDataJoined = soa::Join; + using JetTableMCDJoined = soa::Join; + using JetTableMCDWeightedJoined = soa::Join; + using JetTableMCPJoined = soa::Join; + using JetTableMCPWeightedJoined = soa::Join; -using JetTableMCDMatchedJoined = soa::Join; -using jetMcpPerMcCollision = soa::Join; + using JetTableMCDMatchedJoined = soa::Join; + using jetMcpPerMcCollision = soa::Join; -using JetTableMCDMatchedWeightedJoined = soa::Join; -using JetTableMCPMatchedWeightedJoined = soa::Join; + using JetTableMCDMatchedWeightedJoined = soa::Join; + using JetTableMCPMatchedWeightedJoined = soa::Join; -// Applying some cuts(filters) on collisions, tracks, clusters + // Applying some cuts(filters) on collisions, tracks, clusters -Filter eventCuts = (nabs(aod::jcollision::posZ) < vertexZCut && aod::jcollision::centFT0M >= centralityMin && aod::jcollision::centFT0M < centralityMax); -// Filter EMCeventCuts = (nabs(aod::collision::posZ) < vertexZCut && aod::collision::centrality >= centralityMin && aod::collision::centrality < centralityMax); -Filter trackCuts = (aod::jtrack::pt >= trackpTMin && aod::jtrack::pt < trackpTMax && aod::jtrack::eta > trackEtaMin && aod::jtrack::eta < trackEtaMax && aod::jtrack::phi >= trackPhiMin && aod::jtrack::phi <= trackPhiMax); -aod::EMCALClusterDefinition clusterDefinition = aod::emcalcluster::getClusterDefinitionFromString(clusterDefinitionS.value); -Filter clusterFilter = (aod::jcluster::definition == static_cast(clusterDefinition) && aod::jcluster::eta > clusterEtaMin && aod::jcluster::eta < clusterEtaMax && aod::jcluster::phi >= clusterPhiMin && aod::jcluster::phi <= clusterPhiMax && aod::jcluster::energy >= clusterEnergyMin && aod::jcluster::time > clusterTimeMin && aod::jcluster::time < clusterTimeMax && (clusterRejectExotics && aod::jcluster::isExotic != true)); -Preslice JetMCPPerMcCollision = aod::jet::mcCollisionId; -PresliceUnsorted> CollisionsPerMCPCollision = aod::jmccollisionlb::mcCollisionId; + Filter eventCuts = (nabs(aod::jcollision::posZ) < vertexZCut && aod::jcollision::centFT0M >= centralityMin && aod::jcollision::centFT0M < centralityMax); + // Filter EMCeventCuts = (nabs(aod::collision::posZ) < vertexZCut && aod::collision::centrality >= centralityMin && aod::collision::centrality < centralityMax); + Filter trackCuts = (aod::jtrack::pt >= trackpTMin && aod::jtrack::pt < trackpTMax && aod::jtrack::eta > trackEtaMin && aod::jtrack::eta < trackEtaMax && aod::jtrack::phi >= trackPhiMin && aod::jtrack::phi <= trackPhiMax); + aod::EMCALClusterDefinition clusterDefinition = aod::emcalcluster::getClusterDefinitionFromString(clusterDefinitionS.value); + Filter clusterFilter = (aod::jcluster::definition == static_cast(clusterDefinition) && aod::jcluster::eta > clusterEtaMin && aod::jcluster::eta < clusterEtaMax && aod::jcluster::phi >= clusterPhiMin && aod::jcluster::phi <= clusterPhiMax && aod::jcluster::energy >= clusterEnergyMin && aod::jcluster::time > clusterTimeMin && aod::jcluster::time < clusterTimeMax && (clusterRejectExotics && aod::jcluster::isExotic != true)); + Preslice JetMCPPerMcCollision = aod::jet::mcCollisionId; + PresliceUnsorted> CollisionsPerMCPCollision = aod::jmccollisionlb::mcCollisionId; -template -bool isAcceptedRecoJet(U const& jet) -{ - // if (jetAreaFractionMin > kJetAreaFractionMinThreshold) { - // if (jet.area() < jetAreaFractionMin * o2::constants::math::PI * (jet.r() / 100.0) * (jet.r() / 100.0)) { - // return false; - // } - // } - /* if (leadingTrackPtMin > kLeadingTrackPtMinThreshold) { - bool isMinleadingTrack = false; + template + bool isAcceptedRecoJet(U const& jet) + { + // if (jetAreaFractionMin > kJetAreaFractionMinThreshold) { + // if (jet.area() < jetAreaFractionMin * o2::constants::math::PI * (jet.r() / 100.0) * (jet.r() / 100.0)) { + // return false; + // } + // } + /* if (leadingTrackPtMin > kLeadingTrackPtMinThreshold) { + bool isMinleadingTrack = false; + for (const auto& constituent : jet.template tracks_as()) { + if (constituent.pt() >= leadingTrackPtMin) { + isMinleadingTrack = true; + break; + } + } + + if (!isMinleadingTrack) { + return false; + } + } + if (leadingTrackPtMax < kLeadingTrackPtMaxThreshold) { + bool isMaxleadingTrack = false; for (const auto& constituent : jet.template tracks_as()) { - if (constituent.pt() >= leadingTrackPtMin) { - isMinleadingTrack = true; + if (constituent.pt() <= leadingTrackPtMax) { + isMaxleadingTrack = true; break; -} -} - -if (!isMinleadingTrack) { -return false; -} -} -if (leadingTrackPtMax < kLeadingTrackPtMaxThreshold) { -bool isMaxleadingTrack = false; -for (const auto& constituent : jet.template tracks_as()) { -if (constituent.pt() <= leadingTrackPtMax) { -isMaxleadingTrack = true; -break; -} -} - -if (!isMaxleadingTrack) { -return false; -} -} - -if (leadingClusterPtMin > kLeadingClusterPtMinThreshold) { -bool isMinleadingCluster = false; -for (const auto& cluster : jet.template clusters_as()) { -double clusterpt = cluster.energy() / std::cosh(cluster.eta()); -if (clusterpt >= leadingClusterPtMin) { -isMinleadingCluster = true; -break; -} -} + } + } -if (!isMinleadingCluster) { -return false; -} -} -if (leadingClusterPtMax < kLeadingClusterPtMaxThreshold) { -bool isMaxleadingCluster = false; -for (const auto& cluster : jet.template clusters_as()) { -double clusterpt = cluster.energy() / std::cosh(cluster.eta()); -if (clusterpt <= leadingClusterPtMax) { -isMaxleadingCluster = true; -break; -} -} + if (!isMaxleadingTrack) { + return false; + } + } -if (!isMaxleadingCluster) { -return false; -} -} -return true; -*/ -// --- Track cuts: ALL tracks must satisfy 0.15 <= pT <= 200 or 150 GeV/c--- -if (leadingTrackPtMin > kLeadingTrackPtMinThreshold || leadingTrackPtMax < kLeadingTrackPtMaxThreshold) { - bool hasValidTrack = false; - for (const auto& constituent : jet.template tracks_as()) { - const float pt = constituent.pt(); - // Reject if ANY track fails min or max cut - if ((leadingTrackPtMin > kLeadingTrackPtMinThreshold && pt < leadingTrackPtMin) || - (leadingTrackPtMax < kLeadingTrackPtMaxThreshold && pt > leadingTrackPtMax)) { - return false; - } - hasValidTrack = true; // At least one track exists (if needed) + if (leadingClusterPtMin > kLeadingClusterPtMinThreshold) { + bool isMinleadingCluster = false; + for (const auto& cluster : jet.template clusters_as()) { + double clusterpt = cluster.energy() / std::cosh(cluster.eta()); + if (clusterpt >= leadingClusterPtMin) { + isMinleadingCluster = true; + break; } - // Reject if no tracks exist (edge case) - if (leadingTrackPtMin > kLeadingTrackPtMinThreshold && !hasValidTrack) { - return false; } -} -// --- Cluster cuts: ALL clusters must satisfy min <= pT <= max == 0.3 <= pT <= 250 -if (leadingClusterPtMin > kLeadingClusterPtMinThreshold || leadingClusterPtMax < kLeadingClusterPtMaxThreshold) { - bool hasValidCluster = false; + if (!isMinleadingCluster) { + return false; + } + } + if (leadingClusterPtMax < kLeadingClusterPtMaxThreshold) { + bool isMaxleadingCluster = false; for (const auto& cluster : jet.template clusters_as()) { - const double pt = cluster.energy() / std::cosh(cluster.eta()); - // Reject if ANY cluster fails min or max cut - if ((leadingClusterPtMin > kLeadingClusterPtMinThreshold && pt < leadingClusterPtMin) || - (leadingClusterPtMax < kLeadingClusterPtMaxThreshold && pt > leadingClusterPtMax)) { - return false; - } - hasValidCluster = true; // At least one cluster exists + double clusterpt = cluster.energy() / std::cosh(cluster.eta()); + if (clusterpt <= leadingClusterPtMax) { + isMaxleadingCluster = true; + break; } - // Reject if no clusters exist (edge case) - if (leadingClusterPtMin > kLeadingClusterPtMinThreshold && !hasValidCluster) { - return false; } -} -return true; -} //isAcceptedRecoJet ends -template -bool isAcceptedPartJet(U const& jet) -{ - // if (jetAreaFractionMin > kJetAreaFractionMinThreshold) { - // if (jet.area() < jetAreaFractionMin * o2::constants::math::PI * (jet.r() / 100.0) * (jet.r() / 100.0)) { - // return false; - // } - // } - //track pt Min cut at the Part level: 0 < pT <= 200 or 150 GeV/c - if (leadingTrackPtMin > kLeadingTrackPtMinThreshold || leadingTrackPtMax < kLeadingTrackPtMaxThreshold) { - bool hasValidParticle = false; - for (const auto& constituent : jet.template tracks_as()) { - const float pt = constituent.pt(); - // Reject if ANY particle fails min or max cut - if ((leadingTrackPtMin > kLeadingTrackPtMinThreshold && pt < leadingTrackPtMin) || - (leadingTrackPtMax < kLeadingTrackPtMaxThreshold && pt > leadingTrackPtMax)) { + if (!isMaxleadingCluster) { + return false; + } + } + return true; + */ + // --- Track cuts: ALL tracks must satisfy 0.15 <= pT <= 200 or 150 GeV/c--- + if (leadingTrackPtMin > kLeadingTrackPtMinThreshold || leadingTrackPtMax < kLeadingTrackPtMaxThreshold) { + bool hasValidTrack = false; + for (const auto& constituent : jet.template tracks_as()) { + const float pt = constituent.pt(); + // Reject if ANY track fails min or max cut + if ((leadingTrackPtMin > kLeadingTrackPtMinThreshold && pt < leadingTrackPtMin) || + (leadingTrackPtMax < kLeadingTrackPtMaxThreshold && pt > leadingTrackPtMax)) { + return false; + } + hasValidTrack = true; // At least one track exists (if needed) + } + // Reject if no tracks exist (edge case) + if (leadingTrackPtMin > kLeadingTrackPtMinThreshold && !hasValidTrack) { + return false; + } + } + + // --- Cluster cuts: ALL clusters must satisfy min <= pT <= max == 0.3 <= pT <= 250 + if (leadingClusterPtMin > kLeadingClusterPtMinThreshold || leadingClusterPtMax < kLeadingClusterPtMaxThreshold) { + bool hasValidCluster = false; + for (const auto& cluster : jet.template clusters_as()) { + const double pt = cluster.energy() / std::cosh(cluster.eta()); + // Reject if ANY cluster fails min or max cut + if ((leadingClusterPtMin > kLeadingClusterPtMinThreshold && pt < leadingClusterPtMin) || + (leadingClusterPtMax < kLeadingClusterPtMaxThreshold && pt > leadingClusterPtMax)) { + return false; + } + hasValidCluster = true; // At least one cluster exists + } + // Reject if no clusters exist (edge case) + if (leadingClusterPtMin > kLeadingClusterPtMinThreshold && !hasValidCluster) { return false; } - hasValidParticle = true; // At least one track exists (if needed) } - // Reject if no particle exist (edge case) - if (leadingTrackPtMin > kLeadingTrackPtMinThreshold && !hasValidParticle) { - return false; + return true; + } // isAcceptedRecoJet ends + + template + bool isAcceptedPartJet(U const& jet) + { + // if (jetAreaFractionMin > kJetAreaFractionMinThreshold) { + // if (jet.area() < jetAreaFractionMin * o2::constants::math::PI * (jet.r() / 100.0) * (jet.r() / 100.0)) { + // return false; + // } + // } + // track pt Min cut at the Part level: 0 < pT <= 200 or 150 GeV/c + if (leadingTrackPtMin > kLeadingTrackPtMinThreshold || leadingTrackPtMax < kLeadingTrackPtMaxThreshold) { + bool hasValidParticle = false; + for (const auto& constituent : jet.template tracks_as()) { + const float pt = constituent.pt(); + // Reject if ANY particle fails min or max cut + if ((leadingTrackPtMin > kLeadingTrackPtMinThreshold && pt < leadingTrackPtMin) || + (leadingTrackPtMax < kLeadingTrackPtMaxThreshold && pt > leadingTrackPtMax)) { + return false; + } + hasValidParticle = true; // At least one track exists (if needed) + } + // Reject if no particle exist (edge case) + if (leadingTrackPtMin > kLeadingTrackPtMinThreshold && !hasValidParticle) { + return false; + } } + return true; + } + /* + if (leadingTrackPtMin > kLeadingTrackPtMinThreshold) { + bool isMinleadingTrack = false; + for (const auto& constituent : jet.template tracks_as()) { + auto pdgParticle = pdgDatabase->GetParticle(constituent.pdgCode()); + if (pdgParticle->Charge() != 0 && constituent.pt() >= leadingTrackPtMin) { + isMinleadingTrack = true; + break; + } } - return true; -} -/* -if (leadingTrackPtMin > kLeadingTrackPtMinThreshold) { -bool isMinleadingTrack = false; -for (const auto& constituent : jet.template tracks_as()) { -auto pdgParticle = pdgDatabase->GetParticle(constituent.pdgCode()); -if (pdgParticle->Charge() != 0 && constituent.pt() >= leadingTrackPtMin) { -isMinleadingTrack = true; -break; -} -} - -if (!isMinleadingTrack) { -return false; -} -} -//Leading track pt Max cut -if (leadingTrackPtMax < kLeadingTrackPtMaxThreshold) { -bool isMaxleadingTrack = false; -for (const auto& constituent : jet.template tracks_as()) { -auto pdgParticle = pdgDatabase->GetParticle(constituent.pdgCode()); -if (pdgParticle->Charge() != 0 && constituent.pt() <= leadingTrackPtMax) { -isMaxleadingTrack = true; -break; -} -} -if (!isMaxleadingTrack) { -return false; -} -} -//Leading cluster pt Min cut -if (leadingClusterPtMin > kLeadingClusterPtMinThreshold) { -bool isMinleadingCluster = false; -for (const auto& constituent : jet.template clusters_as()) { -auto pdgParticle = pdgDatabase->GetParticle(constituent.pdgCode()); -if (pdgParticle->Charge() == 0) { -double clusterPt = constituent.e() / std::cosh(constituent.eta()); -if (clusterPt >= leadingClusterPtMin) { -isMinleadingCluster = true; -break; -} -} -} -if (!isMinleadingCluster) { -return false; -} -} -//Leading cluster pt Max cut -if (leadingClusterPtMax < kLeadingClusterPtMaxThreshold) { -bool isMaxleadingCluster = false; -for (const auto& constituent : jet.template clusters_as()) { -auto pdgParticle = pdgDatabase->GetParticle(constituent.pdgCode()); -if (pdgParticle->Charge() == 0) { -double clusterPt = constituent.e() / std::cosh(constituent.eta()); -if (clusterPt <= leadingClusterPtMax) { -isMaxleadingCluster = true; -break; -} -} -} -if (!isMaxleadingCluster) { -return false; -} -} -return true; -} -*/ + if (!isMinleadingTrack) { + return false; + } + } + //Leading track pt Max cut + if (leadingTrackPtMax < kLeadingTrackPtMaxThreshold) { + bool isMaxleadingTrack = false; + for (const auto& constituent : jet.template tracks_as()) { + auto pdgParticle = pdgDatabase->GetParticle(constituent.pdgCode()); + if (pdgParticle->Charge() != 0 && constituent.pt() <= leadingTrackPtMax) { + isMaxleadingTrack = true; + break; + } + } -template -void fillJetHistograms(T const& jet, float weight = 1.0) -{ - float neutralEnergy = 0.0; - double sumtrackE = 0.0; - if (jet.r() == round(selectedJetsRadius * 100.0f)) { - registry.fill(HIST("h_full_jet_pt"), jet.pt(), weight); - registry.fill(HIST("h_full_jet_eta"), jet.eta(), weight); - registry.fill(HIST("h_full_jet_phi"), jet.phi(), weight); - registry.fill(HIST("h2_jet_etaphi"), jet.eta(), jet.phi(), weight); - - for (const auto& cluster : jet.template clusters_as()) { - registry.fill(HIST("h2_full_jet_neutralconstituents"), jet.pt(), jet.clustersIds().size(), weight); - - neutralEnergy += cluster.energy(); - double clusterpt = cluster.energy() / std::cosh(cluster.eta()); - registry.fill(HIST("h_full_jet_clusterTime"), cluster.time(), weight); - registry.fill(HIST("h_full_jet_neutralconstituents_pt"), clusterpt, weight); - registry.fill(HIST("h_full_jet_neutralconstituents_eta"), cluster.eta(), weight); - registry.fill(HIST("h_full_jet_neutralconstituents_phi"), cluster.phi(), weight); - registry.fill(HIST("h_full_jet_neutralconstituents_energy"), cluster.energy(), weight); - registry.fill(HIST("h_full_jet_neutralconstituents_energysum"), neutralEnergy, weight); - } - auto nef = neutralEnergy / jet.energy(); - registry.fill(HIST("h2_full_jet_nef"), jet.pt(), nef, weight); - - for (const auto& jettrack : jet.template tracks_as()) { - sumtrackE += jettrack.energy(); - - registry.fill(HIST("h_Detjet_ntracks"), jettrack.pt(), weight); - registry.fill(HIST("h2_full_jet_chargedconstituents"), jet.pt(), jet.tracksIds().size(), weight); - registry.fill(HIST("h2_full_jettrack_pt"), jet.pt(), jettrack.pt(), weight); - registry.fill(HIST("h2_full_jettrack_eta"), jet.eta(), jettrack.eta(), weight); - registry.fill(HIST("h2_full_jettrack_phi"), jet.phi(), jettrack.phi(), weight); - - registry.fill(HIST("h2_track_etaphi"), jettrack.eta(), jettrack.phi(), weight); - registry.fill(HIST("h_full_jet_chargedconstituents_pt"), jettrack.pt(), weight); - registry.fill(HIST("h_full_jet_chargedconstituents_eta"), jettrack.eta(), weight); - registry.fill(HIST("h_full_jet_chargedconstituents_phi"), jettrack.phi(), weight); - registry.fill(HIST("h_full_jet_chargedconstituents_energy"), jettrack.energy(), weight); - registry.fill(HIST("h_full_jet_chargedconstituents_energysum"), sumtrackE, weight); - } - } // jet.r() -} + if (!isMaxleadingTrack) { + return false; + } + } + //Leading cluster pt Min cut + if (leadingClusterPtMin > kLeadingClusterPtMinThreshold) { + bool isMinleadingCluster = false; + for (const auto& constituent : jet.template clusters_as()) { + auto pdgParticle = pdgDatabase->GetParticle(constituent.pdgCode()); + if (pdgParticle->Charge() == 0) { + double clusterPt = constituent.e() / std::cosh(constituent.eta()); + if (clusterPt >= leadingClusterPtMin) { + isMinleadingCluster = true; + break; + } + } + } + if (!isMinleadingCluster) { + return false; + } + } + //Leading cluster pt Max cut + if (leadingClusterPtMax < kLeadingClusterPtMaxThreshold) { + bool isMaxleadingCluster = false; + for (const auto& constituent : jet.template clusters_as()) { + auto pdgParticle = pdgDatabase->GetParticle(constituent.pdgCode()); + if (pdgParticle->Charge() == 0) { + double clusterPt = constituent.e() / std::cosh(constituent.eta()); + if (clusterPt <= leadingClusterPtMax) { + isMaxleadingCluster = true; + break; + } + } + } + if (!isMaxleadingCluster) { + return false; + } + } + return true; + } + */ -// check for nef distribution for rejected events -template -void fillRejectedJetHistograms(T const& jet, float weight = 1.0) -{ - float neutralEnergy = 0.0; - if (jet.r() == round(selectedJetsRadius * 100.0f)) { - for (const auto& cluster : jet.template clusters_as()) { - neutralEnergy += cluster.energy(); - } - auto nef = neutralEnergy / jet.energy(); - registry.fill(HIST("h2_full_jet_nef_rejected"), jet.pt(), nef, weight); - } // jet.r() -} + template + void fillJetHistograms(T const& jet, float weight = 1.0) + { + float neutralEnergy = 0.0; + double sumtrackE = 0.0; + if (jet.r() == round(selectedJetsRadius * 100.0f)) { + registry.fill(HIST("h_full_jet_pt"), jet.pt(), weight); + registry.fill(HIST("h_full_jet_eta"), jet.eta(), weight); + registry.fill(HIST("h_full_jet_phi"), jet.phi(), weight); + registry.fill(HIST("h2_jet_etaphi"), jet.eta(), jet.phi(), weight); + + for (const auto& cluster : jet.template clusters_as()) { + registry.fill(HIST("h2_full_jet_neutralconstituents"), jet.pt(), jet.clustersIds().size(), weight); + + neutralEnergy += cluster.energy(); + double clusterpt = cluster.energy() / std::cosh(cluster.eta()); + registry.fill(HIST("h_full_jet_clusterTime"), cluster.time(), weight); + registry.fill(HIST("h_full_jet_neutralconstituents_pt"), clusterpt, weight); + registry.fill(HIST("h_full_jet_neutralconstituents_eta"), cluster.eta(), weight); + registry.fill(HIST("h_full_jet_neutralconstituents_phi"), cluster.phi(), weight); + registry.fill(HIST("h_full_jet_neutralconstituents_energy"), cluster.energy(), weight); + registry.fill(HIST("h_full_jet_neutralconstituents_energysum"), neutralEnergy, weight); + } + auto nef = neutralEnergy / jet.energy(); + registry.fill(HIST("h2_full_jet_nef"), jet.pt(), nef, weight); + + for (const auto& jettrack : jet.template tracks_as()) { + sumtrackE += jettrack.energy(); + + registry.fill(HIST("h_Detjet_ntracks"), jettrack.pt(), weight); + registry.fill(HIST("h2_full_jet_chargedconstituents"), jet.pt(), jet.tracksIds().size(), weight); + registry.fill(HIST("h2_full_jettrack_pt"), jet.pt(), jettrack.pt(), weight); + registry.fill(HIST("h2_full_jettrack_eta"), jet.eta(), jettrack.eta(), weight); + registry.fill(HIST("h2_full_jettrack_phi"), jet.phi(), jettrack.phi(), weight); + + registry.fill(HIST("h2_track_etaphi"), jettrack.eta(), jettrack.phi(), weight); + registry.fill(HIST("h_full_jet_chargedconstituents_pt"), jettrack.pt(), weight); + registry.fill(HIST("h_full_jet_chargedconstituents_eta"), jettrack.eta(), weight); + registry.fill(HIST("h_full_jet_chargedconstituents_phi"), jettrack.phi(), weight); + registry.fill(HIST("h_full_jet_chargedconstituents_energy"), jettrack.energy(), weight); + registry.fill(HIST("h_full_jet_chargedconstituents_energysum"), sumtrackE, weight); + } + } // jet.r() + } -template -void fillMCPHistograms(T const& jet, float weight = 1.0) -{ - float neutralEnergy = 0.0; - int neutralconsts = 0; - int chargedconsts = 0; - int mcpjetOutsideFid = 0; - int mcpjetInsideFid = 0; - - auto isInFiducial = [&](auto const& jet) { - return jet.eta() >= jetEtaMin && jet.eta() <= jetEtaMax && - jet.phi() >= jetPhiMin && jet.phi() <= jetPhiMax; - }; - - if (jet.r() == round(selectedJetsRadius * 100.0f)) { - registry.fill(HIST("h_full_mcpjet_tablesize"), jet.size(), weight); - registry.fill(HIST("h_full_mcpjet_ntracks"), jet.tracksIds().size(), weight); - registry.fill(HIST("h_full_jet_pt_part"), jet.pt(), weight); - registry.fill(HIST("h_full_jet_eta_part"), jet.eta(), weight); - registry.fill(HIST("h_full_jet_phi_part"), jet.phi(), weight); - registry.fill(HIST("h2_jet_etaphi_part"), jet.eta(), jet.phi(), weight); - - if (!isInFiducial(jet)) { - // jet is outside - mcpjetOutsideFid++; - registry.fill(HIST("h2_full_mcpjetOutsideFiducial_pt"), jet.pt(), mcpjetOutsideFid, weight); - registry.fill(HIST("h_full_mcpjetOutside_eta_part"), jet.eta(), weight); - registry.fill(HIST("h_full_mcpjetOutside_phi_part"), jet.phi(), weight); - } else { - // jet is inside - mcpjetInsideFid++; - registry.fill(HIST("h2_full_mcpjetInsideFiducial_pt"), jet.pt(), mcpjetInsideFid, weight); - registry.fill(HIST("h_full_mcpjetInside_eta_part"), jet.eta(), weight); - registry.fill(HIST("h_full_mcpjetInside_phi_part"), jet.phi(), weight); - } - - for (const auto& constituent : jet.template tracks_as()) { - auto pdgParticle = pdgDatabase->GetParticle(constituent.pdgCode()); - if (pdgParticle->Charge() == 0) { - neutralconsts++; - neutralEnergy += constituent.e(); // neutral jet constituents at particle level - double clusterpt = constituent.e() / std::cosh(constituent.eta()); - registry.fill(HIST("h2_full_jet_neutralconstituents_part"), jet.pt(), neutralconsts, weight); - registry.fill(HIST("h_full_jet_neutralconstituents_pt_part"), clusterpt, weight); - registry.fill(HIST("h_full_jet_neutralconstituents_eta_part"), constituent.eta(), weight); - registry.fill(HIST("h_full_jet_neutralconstituents_phi_part"), constituent.phi(), weight); - registry.fill(HIST("h_full_jet_neutralconstituents_energy_part"), constituent.e(), weight); - registry.fill(HIST("h_full_jet_neutralconstituents_energysum_part"), neutralEnergy, weight); + // check for nef distribution for rejected events + template + void fillRejectedJetHistograms(T const& jet, float weight = 1.0) + { + float neutralEnergy = 0.0; + if (jet.r() == round(selectedJetsRadius * 100.0f)) { + for (const auto& cluster : jet.template clusters_as()) { + neutralEnergy += cluster.energy(); + } + auto nef = neutralEnergy / jet.energy(); + registry.fill(HIST("h2_full_jet_nef_rejected"), jet.pt(), nef, weight); + } // jet.r() + } + template + void fillMCPHistograms(T const& jet, float weight = 1.0) + { + float neutralEnergy = 0.0; + int neutralconsts = 0; + int chargedconsts = 0; + int mcpjetOutsideFid = 0; + int mcpjetInsideFid = 0; + + auto isInFiducial = [&](auto const& jet) { + return jet.eta() >= jetEtaMin && jet.eta() <= jetEtaMax && + jet.phi() >= jetPhiMin && jet.phi() <= jetPhiMax; + }; + + if (jet.r() == round(selectedJetsRadius * 100.0f)) { + registry.fill(HIST("h_full_mcpjet_tablesize"), jet.size(), weight); + registry.fill(HIST("h_full_mcpjet_ntracks"), jet.tracksIds().size(), weight); + registry.fill(HIST("h_full_jet_pt_part"), jet.pt(), weight); + registry.fill(HIST("h_full_jet_eta_part"), jet.eta(), weight); + registry.fill(HIST("h_full_jet_phi_part"), jet.phi(), weight); + registry.fill(HIST("h2_jet_etaphi_part"), jet.eta(), jet.phi(), weight); + + if (!isInFiducial(jet)) { + // jet is outside + mcpjetOutsideFid++; + registry.fill(HIST("h2_full_mcpjetOutsideFiducial_pt"), jet.pt(), mcpjetOutsideFid, weight); + registry.fill(HIST("h_full_mcpjetOutside_eta_part"), jet.eta(), weight); + registry.fill(HIST("h_full_mcpjetOutside_phi_part"), jet.phi(), weight); } else { - chargedconsts++; - registry.fill(HIST("h2_full_jet_chargedconstituents_part"), jet.pt(), chargedconsts, weight); // charged jet constituents at particle level - registry.fill(HIST("h2_jettrack_pt_part"), jet.pt(), constituent.pt(), weight); - registry.fill(HIST("h2_jettrack_eta_part"), jet.eta(), constituent.eta(), weight); - registry.fill(HIST("h2_jettrack_phi_part"), jet.phi(), constituent.phi(), weight); - registry.fill(HIST("h2_track_etaphi_part"), constituent.eta(), constituent.phi(), weight); - } - } // constituent loop - auto nef = neutralEnergy / jet.energy(); - registry.fill(HIST("h2_full_jet_nef_part"), jet.pt(), nef, weight); - } // jet.r() -} - -template -void fillTrackHistograms(T const& tracks, U const& clusters, float weight = 1.0) -{ - double sumtrackE = 0.0; - - for (auto const& track : tracks) { - if (!jetderiveddatautilities::selectTrack(track, trackSelection)) { - continue; - } - sumtrackE += track.energy(); - registry.fill(HIST("h_track_pt"), track.pt(), weight); - registry.fill(HIST("h_track_eta"), track.eta(), weight); - registry.fill(HIST("h_track_phi"), track.phi(), weight); - registry.fill(HIST("h_track_energysum"), sumtrackE, weight); - } - double sumclusterE = 0.0; - for (auto const& cluster : clusters) { - double clusterpt = cluster.energy() / std::cosh(cluster.eta()); - sumclusterE += cluster.energy(); - - registry.fill(HIST("h_clusterTime"), cluster.time(), weight); - registry.fill(HIST("h_cluster_pt"), clusterpt, weight); - registry.fill(HIST("h_cluster_eta"), cluster.eta(), weight); - registry.fill(HIST("h_cluster_phi"), cluster.phi(), weight); - registry.fill(HIST("h_cluster_energy"), cluster.energy(), weight); - registry.fill(HIST("h_cluster_energysum"), sumclusterE, weight); - } -} + // jet is inside + mcpjetInsideFid++; + registry.fill(HIST("h2_full_mcpjetInsideFiducial_pt"), jet.pt(), mcpjetInsideFid, weight); + registry.fill(HIST("h_full_mcpjetInside_eta_part"), jet.eta(), weight); + registry.fill(HIST("h_full_mcpjetInside_phi_part"), jet.phi(), weight); + } -template -void fillMatchedHistograms(T const& jetBase, float weight = 1.0) -{ - if (jetBase.has_matchedJetGeo()) { // geometrical jet matching only needed for pp - here,matching Base(Det.level) with Tag (Part. level) jets - registry.fill(HIST("h_full_matchedmcdjet_tablesize"), jetBase.size(), weight); - registry.fill(HIST("h_full_matchedmcdjet_ntracks"), jetBase.tracksIds().size(), weight); - registry.fill(HIST("h2_matchedjet_etaphiDet"), jetBase.eta(), jetBase.phi(), weight); - - for (const auto& jetTag : jetBase.template matchedJetGeo_as>()) { - auto deltaEta = jetBase.eta() - jetTag.eta(); - auto deltaPhi = jetBase.phi() - jetTag.phi(); - auto deltaR = jetutilities::deltaR(jetBase, jetTag); - - registry.fill(HIST("h_full_jet_deltaR"), deltaR, weight); - registry.fill(HIST("h_full_matchedmcpjet_tablesize"), jetTag.size(), weight); - registry.fill(HIST("h_full_matchedmcpjet_ntracks"), jetTag.tracksIds().size(), weight); - registry.fill(HIST("h2_matchedjet_etaphiPart"), jetTag.eta(), jetTag.phi(), weight); - registry.fill(HIST("h2_matchedjet_deltaEtaCorr"), jetBase.eta(), jetTag.eta(), weight); - registry.fill(HIST("h2_matchedjet_deltaPhiCorr"), jetBase.phi(), jetTag.phi(), weight); - - // JES for fulljets - registry.fill(HIST("h2_full_jet_energyscaleDet"), jetBase.pt(), (jetBase.pt() - jetTag.pt()) / jetTag.pt(), weight); - registry.fill(HIST("h2_full_jet_energyscalePart"), jetTag.pt(), (jetBase.pt() - jetTag.pt()) / jetTag.pt(), weight); - registry.fill(HIST("h3_full_jet_energyscalePart"), jetBase.r() / 100.0, jetTag.pt(), (jetBase.pt() - jetTag.pt()) / jetTag.pt(), weight); - registry.fill(HIST("h2_full_jet_etaresolutionPart"), jetTag.pt(), deltaEta / jetTag.eta(), weight); - registry.fill(HIST("h2_full_jet_phiresolutionPart"), jetTag.pt(), deltaPhi / jetTag.phi(), weight); + for (const auto& constituent : jet.template tracks_as()) { + auto pdgParticle = pdgDatabase->GetParticle(constituent.pdgCode()); + if (pdgParticle->Charge() == 0) { + neutralconsts++; + neutralEnergy += constituent.e(); // neutral jet constituents at particle level + double clusterpt = constituent.e() / std::cosh(constituent.eta()); + registry.fill(HIST("h2_full_jet_neutralconstituents_part"), jet.pt(), neutralconsts, weight); + registry.fill(HIST("h_full_jet_neutralconstituents_pt_part"), clusterpt, weight); + registry.fill(HIST("h_full_jet_neutralconstituents_eta_part"), constituent.eta(), weight); + registry.fill(HIST("h_full_jet_neutralconstituents_phi_part"), constituent.phi(), weight); + registry.fill(HIST("h_full_jet_neutralconstituents_energy_part"), constituent.e(), weight); + registry.fill(HIST("h_full_jet_neutralconstituents_energysum_part"), neutralEnergy, weight); - // Response Matrix - registry.fill(HIST("h_full_jet_ResponseMatrix"), jetBase.pt(), jetTag.pt(), weight); // MCD vs MCP jet pT - } // jetTag - } // jetBase -} + } else { + chargedconsts++; + registry.fill(HIST("h2_full_jet_chargedconstituents_part"), jet.pt(), chargedconsts, weight); // charged jet constituents at particle level + registry.fill(HIST("h2_jettrack_pt_part"), jet.pt(), constituent.pt(), weight); + registry.fill(HIST("h2_jettrack_eta_part"), jet.eta(), constituent.eta(), weight); + registry.fill(HIST("h2_jettrack_phi_part"), jet.phi(), constituent.phi(), weight); + registry.fill(HIST("h2_track_etaphi_part"), constituent.eta(), constituent.phi(), weight); + } + } // constituent loop + auto nef = neutralEnergy / jet.energy(); + registry.fill(HIST("h2_full_jet_nef_part"), jet.pt(), nef, weight); + } // jet.r() + } -void processDummy(aod::JetCollisions const&) -{ -} -PROCESS_SWITCH(FullJetSpectra, processDummy, "dummy task", true); + template + void fillTrackHistograms(T const& tracks, U const& clusters, float weight = 1.0) + { + double sumtrackE = 0.0; -void processJetsData(soa::Filtered::iterator const& collision, FullJetTableDataJoined const& jets, aod::JetTracks const&, aod::JetClusters const&) -{ - bool eventAccepted = false; + for (auto const& track : tracks) { + if (!jetderiveddatautilities::selectTrack(track, trackSelection)) { + continue; + } + sumtrackE += track.energy(); + registry.fill(HIST("h_track_pt"), track.pt(), weight); + registry.fill(HIST("h_track_eta"), track.eta(), weight); + registry.fill(HIST("h_track_phi"), track.phi(), weight); + registry.fill(HIST("h_track_energysum"), sumtrackE, weight); + } + double sumclusterE = 0.0; + for (auto const& cluster : clusters) { + double clusterpt = cluster.energy() / std::cosh(cluster.eta()); + sumclusterE += cluster.energy(); - registry.fill(HIST("hDetcollisionCounter"), 0.5); // allDetColl - if (std::fabs(collision.posZ()) > vertexZCut) { - return; + registry.fill(HIST("h_clusterTime"), cluster.time(), weight); + registry.fill(HIST("h_cluster_pt"), clusterpt, weight); + registry.fill(HIST("h_cluster_eta"), cluster.eta(), weight); + registry.fill(HIST("h_cluster_phi"), cluster.phi(), weight); + registry.fill(HIST("h_cluster_energy"), cluster.energy(), weight); + registry.fill(HIST("h_cluster_energysum"), sumclusterE, weight); + } } - registry.fill(HIST("hDetcollisionCounter"), 1.5); // DetCollWithVertexZ - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { - registry.fill(HIST("hDetcollisionCounter"), 2.5); // EventsNotSatisfyingEventSelection - return; + template + void fillMatchedHistograms(T const& jetBase, float weight = 1.0) + { + if (jetBase.has_matchedJetGeo()) { // geometrical jet matching only needed for pp - here,matching Base(Det.level) with Tag (Part. level) jets + registry.fill(HIST("h_full_matchedmcdjet_tablesize"), jetBase.size(), weight); + registry.fill(HIST("h_full_matchedmcdjet_ntracks"), jetBase.tracksIds().size(), weight); + registry.fill(HIST("h2_matchedjet_etaphiDet"), jetBase.eta(), jetBase.phi(), weight); + + for (const auto& jetTag : jetBase.template matchedJetGeo_as>()) { + auto deltaEta = jetBase.eta() - jetTag.eta(); + auto deltaPhi = jetBase.phi() - jetTag.phi(); + auto deltaR = jetutilities::deltaR(jetBase, jetTag); + + registry.fill(HIST("h_full_jet_deltaR"), deltaR, weight); + registry.fill(HIST("h_full_matchedmcpjet_tablesize"), jetTag.size(), weight); + registry.fill(HIST("h_full_matchedmcpjet_ntracks"), jetTag.tracksIds().size(), weight); + registry.fill(HIST("h2_matchedjet_etaphiPart"), jetTag.eta(), jetTag.phi(), weight); + registry.fill(HIST("h2_matchedjet_deltaEtaCorr"), jetBase.eta(), jetTag.eta(), weight); + registry.fill(HIST("h2_matchedjet_deltaPhiCorr"), jetBase.phi(), jetTag.phi(), weight); + + // JES for fulljets + registry.fill(HIST("h2_full_jet_energyscaleDet"), jetBase.pt(), (jetBase.pt() - jetTag.pt()) / jetTag.pt(), weight); + registry.fill(HIST("h2_full_jet_energyscalePart"), jetTag.pt(), (jetBase.pt() - jetTag.pt()) / jetTag.pt(), weight); + registry.fill(HIST("h3_full_jet_energyscalePart"), jetBase.r() / 100.0, jetTag.pt(), (jetBase.pt() - jetTag.pt()) / jetTag.pt(), weight); + registry.fill(HIST("h2_full_jet_etaresolutionPart"), jetTag.pt(), deltaEta / jetTag.eta(), weight); + registry.fill(HIST("h2_full_jet_phiresolutionPart"), jetTag.pt(), deltaPhi / jetTag.phi(), weight); + + // Response Matrix + registry.fill(HIST("h_full_jet_ResponseMatrix"), jetBase.pt(), jetTag.pt(), weight); // MCD vs MCP jet pT + } // jetTag + } // jetBase + } + + void processDummy(aod::JetCollisions const&) + { } - if (doEMCALEventWorkaround) { - if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content - eventAccepted = true; - if (collision.alias_bit(kTVXinEMC)) { + PROCESS_SWITCH(FullJetSpectra, processDummy, "dummy task", true); + + void processJetsData(soa::Filtered::iterator const& collision, FullJetTableDataJoined const& jets, aod::JetTracks const&, aod::JetClusters const&) + { + bool eventAccepted = false; + + registry.fill(HIST("hDetcollisionCounter"), 0.5); // allDetColl + if (std::fabs(collision.posZ()) > vertexZCut) { + return; + } + registry.fill(HIST("hDetcollisionCounter"), 1.5); // DetCollWithVertexZ + + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { + registry.fill(HIST("hDetcollisionCounter"), 2.5); // EventsNotSatisfyingEventSelection + return; + } + if (doEMCALEventWorkaround) { + if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content + eventAccepted = true; + if (collision.alias_bit(kTVXinEMC)) { + registry.fill(HIST("hDetcollisionCounter"), 3.5); // EMCreadoutDetEventsWithkTVXinEMC + } + } + } else { + if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { registry.fill(HIST("hDetcollisionCounter"), 3.5); // EMCreadoutDetEventsWithkTVXinEMC + eventAccepted = true; } } - } else { - if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { - registry.fill(HIST("hDetcollisionCounter"), 3.5); // EMCreadoutDetEventsWithkTVXinEMC - eventAccepted = true; + + if (!eventAccepted) { + for (auto const& jet : jets) { + if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax) || !isAcceptedRecoJet(jet)) { + fillRejectedJetHistograms(jet, 1.0); + } + } + registry.fill(HIST("hDetcollisionCounter"), 4.5); // AllRejectedDetEventsAfterEMCEventSelection + return; } - } + registry.fill(HIST("hDetcollisionCounter"), 5.5); // EMCAcceptedDetColl - if (!eventAccepted) { for (auto const& jet : jets) { - if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax) || !isAcceptedRecoJet(jet)) { - fillRejectedJetHistograms(jet, 1.0); + if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { + continue; + } + if (jet.phi() < jetPhiMin || jet.phi() > jetPhiMax) { + continue; + } + if (!isAcceptedRecoJet(jet)) { + continue; } + fillJetHistograms(jet); } - registry.fill(HIST("hDetcollisionCounter"), 4.5); // AllRejectedDetEventsAfterEMCEventSelection - return; } - registry.fill(HIST("hDetcollisionCounter"), 5.5); // EMCAcceptedDetColl + PROCESS_SWITCH(FullJetSpectra, processJetsData, "Full Jets Data", false); - for (auto const& jet : jets) { - if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { - continue; - } - if (jet.phi() < jetPhiMin || jet.phi() > jetPhiMax) { - continue; - } - if (!isAcceptedRecoJet(jet)) { - continue; - } - fillJetHistograms(jet); - } -} -PROCESS_SWITCH(FullJetSpectra, processJetsData, "Full Jets Data", false); + void processJetsTriggeredData(soa::Filtered::iterator const& collision, FullJetTableDataJoined const& /*jets*/, + aod::JetTracks const&, aod::JetClusters const&, aod::JBCs const&) + { + bool eventAccepted = false; -void processJetsTriggeredData(soa::Filtered::iterator const& collision, FullJetTableDataJoined const& /*jets*/, - aod::JetTracks const&, aod::JetClusters const&, aod::JBCs const&) -{ - bool eventAccepted = false; + registry.fill(HIST("hDetTrigcollisionCounter"), 0.5); // allDetTrigColl - registry.fill(HIST("hDetTrigcollisionCounter"), 0.5); // allDetTrigColl + // Get BC info associated with the collision before applying any event selections + auto bc = collision.bc_as(); + // Initialize CCDB objects using the BC info + initCCDB(bc); + // If SoftwareTriggerSelection (i.e. skimming) is enabled, skip this event unless it passes Zorro selection + if (doSoftwareTriggerSelection && !zorro.isSelected(bc.globalBC())) { + return; + } + registry.fill(HIST("hDetTrigcollisionCounter"), 1.5); // DetTrigCollAfterZorroSelection - //Get BC info associated with the collision before applying any event selections - auto bc = collision.bc_as(); - // Initialize CCDB objects using the BC info - initCCDB(bc); - // If SoftwareTriggerSelection (i.e. skimming) is enabled, skip this event unless it passes Zorro selection - if (doSoftwareTriggerSelection && !zorro.isSelected(bc.globalBC())) { - return; - } - registry.fill(HIST("hDetTrigcollisionCounter"), 1.5); // DetTrigCollAfterZorroSelection + if (std::fabs(collision.posZ()) > vertexZCut) { + return; + } + registry.fill(HIST("hDetTrigcollisionCounter"), 2.5); // DetTrigCollWithVertexZ - if (std::fabs(collision.posZ()) > vertexZCut) { - return; + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || !jetderiveddatautilities::selectTrigger(collision, triggerMaskBits)) { + registry.fill(HIST("hDetTrigcollisionCounter"), 3.5); // EventsNotSatisfyingEvent+TriggerSelection + return; + } + //- should this kTVX HW trigger be still in place?? + if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { + eventAccepted = true; + registry.fill(HIST("hDetTrigcollisionCounter"), 4.5); // EMCreadoutDetTrigEventsWithkTVXinEMC + } + // split event selections based on selected triggers - + // make sure there're no trigger overlaps: when analysing JetFullHighPt-> check no JetFullLowPt and kTVXinEMC are fired + // when analysing JetFullLowPt, check kTVXinEMC isn't fired! + // apply exclusive trigger bit selections for event selection + // use a flag and fill QA for trigger overlap - + // - how often you reject a higher pT trig because lower trigs were fired : 5 cases -> 2D hist as a funtn of jet pT + // - check how often the ChJet Trigs were fired for every fullJetTrig fired.(don't reject these events but only for QA) + + /*if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::EMCALReadout)) { + registry.fill(HIST("hDetTrigcollisionCounter"), 4.5); // EMCreadoutDetTrigEventsWithMBTrigs + eventAccepted = true; } - registry.fill(HIST("hDetTrigcollisionCounter"), 2.5); // DetTrigCollWithVertexZ - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || !jetderiveddatautilities::selectTrigger(collision, triggerMaskBits)) { - registry.fill(HIST("hDetTrigcollisionCounter"), 3.5); // EventsNotSatisfyingEvent+TriggerSelection - return; + if (jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::JetFullLowPt)) { + registry.fill(HIST("hDetTrigcollisionCounter"), 6.5); // EMCAcceptedDetTrigCollWithLowFullJetTriggers + eventAccepted = true; } - //- should this kTVX HW trigger be still in place?? - if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { - eventAccepted = true; - registry.fill(HIST("hDetTrigcollisionCounter"), 4.5); // EMCreadoutDetTrigEventsWithkTVXinEMC - } - //split event selections based on selected triggers - - // make sure there're no trigger overlaps: when analysing JetFullHighPt-> check no JetFullLowPt and kTVXinEMC are fired - // when analysing JetFullLowPt, check kTVXinEMC isn't fired! - //apply exclusive trigger bit selections for event selection - // use a flag and fill QA for trigger overlap - - // - how often you reject a higher pT trig because lower trigs were fired : 5 cases -> 2D hist as a funtn of jet pT - // - check how often the ChJet Trigs were fired for every fullJetTrig fired.(don't reject these events but only for QA) - - /*if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::EMCALReadout)) { - registry.fill(HIST("hDetTrigcollisionCounter"), 4.5); // EMCreadoutDetTrigEventsWithMBTrigs + if (jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::JetFullHighPt)) { + registry.fill(HIST("hDetTrigcollisionCounter"), 7.5); // EMCAcceptedDetTrigCollWithHighFullJetTriggers eventAccepted = true; -} + } + */ + // Get trigger status + bool hasFullJetHighPt = jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::JetFullHighPt); + bool hasFullJetLowPt = jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::JetFullLowPt); + bool hasMB = jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::EMCALReadout); -if (jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::JetFullLowPt)) { -registry.fill(HIST("hDetTrigcollisionCounter"), 6.5); // EMCAcceptedDetTrigCollWithLowFullJetTriggers -eventAccepted = true; -} -if (jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::JetFullHighPt)) { -registry.fill(HIST("hDetTrigcollisionCounter"), 7.5); // EMCAcceptedDetTrigCollWithHighFullJetTriggers -eventAccepted = true; -} -*/ -// Get trigger status -bool hasFullJetHighPt = jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::JetFullHighPt); -bool hasFullJetLowPt = jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::JetFullLowPt); -bool hasMB = jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::EMCALReadout); + //*****Step 1: Check for pure triggers (NO overlaps)***** -//*****Step 1: Check for pure triggers (NO overlaps)***** + // Case 1: hasFullJetHighPt && !hasFullJetLowPt && !hasMB : Pure FullJetHighPt + // i.e. for every JetFullHighPt trig that was fired, check the low triggers weren't fired + if (hasFullJetHighPt && !hasFullJetLowPt && !hasMB) { + registry.fill(HIST("hDetTrigcollisionCounter"), 5.5); // OnlyHighPt+NoLowPt+NoMB + } + // Case 2: hasFullJetLowPt && !hasMB : Pure FullJetLowPt + // i.e. for every hasFullJetLowPt trig that was fired, check the MB trig wasn't fired + if (hasFullJetLowPt && !hasMB) { + registry.fill(HIST("hDetTrigcollisionCounter"), 6.5); // OnlyLowPt+NoMB + } + // Case 3: hasMB && !hasFullJetLowPt && !hasFullJetHighPt : Pure MB + // i.e. for every MB trig that was fired, check the higher trigs weren't fired + if (hasMB && !hasFullJetLowPt && !hasFullJetHighPt) { + registry.fill(HIST("hDetTrigcollisionCounter"), 7.5); // OnlyMB + } -//Case 1: hasFullJetHighPt && !hasFullJetLowPt && !hasMB : Pure FullJetHighPt -// i.e. for every JetFullHighPt trig that was fired, check the low triggers weren't fired -if (hasFullJetHighPt && !hasFullJetLowPt && !hasMB) { - registry.fill(HIST("hDetTrigcollisionCounter"), 5.5); //OnlyHighPt+NoLowPt+NoMB -} -//Case 2: hasFullJetLowPt && !hasMB : Pure FullJetLowPt -// i.e. for every hasFullJetLowPt trig that was fired, check the MB trig wasn't fired -if (hasFullJetLowPt && !hasMB) { - registry.fill(HIST("hDetTrigcollisionCounter"), 6.5); //OnlyLowPt+NoMB -} -//Case 3: hasMB && !hasFullJetLowPt && !hasFullJetHighPt : Pure MB -// i.e. for every MB trig that was fired, check the higher trigs weren't fired -if (hasMB && !hasFullJetLowPt && !hasFullJetHighPt) { - registry.fill(HIST("hDetTrigcollisionCounter"), 7.5); //OnlyMB -} + //*****Step 2: Check for trigger overlap cases (for QA):***** -//*****Step 2: Check for trigger overlap cases (for QA):***** + if (hasFullJetHighPt && hasFullJetLowPt) { + registry.fill(HIST("hDetTrigcollisionCounter"), 8.5); // FullJetHighPt+FullJetLowPt + } + if (hasFullJetHighPt && hasMB) { + registry.fill(HIST("hDetTrigcollisionCounter"), 9.5); // FullJetHighPt+MB + } + if (hasFullJetLowPt && hasMB) { + registry.fill(HIST("hDetTrigcollisionCounter"), 10.5); // FullJetLowPt+MB + } -if (hasFullJetHighPt && hasFullJetLowPt) { - registry.fill(HIST("hDetTrigcollisionCounter"), 8.5); //FullJetHighPt+FullJetLowPt -} -if (hasFullJetHighPt && hasMB) { - registry.fill(HIST("hDetTrigcollisionCounter"), 9.5); //FullJetHighPt+MB -} -if (hasFullJetLowPt && hasMB) { - registry.fill(HIST("hDetTrigcollisionCounter"), 10.5); //FullJetLowPt+MB -} + //*****Step 3: Reject ALL overlapping events by applying EXCLUSIVE Trigger Selections ***** + // Skip further processing if ANY overlaps exist + if ((hasFullJetHighPt && (hasFullJetLowPt || hasMB)) || (hasFullJetLowPt && hasMB)) { + registry.fill(HIST("hDetTrigcollisionCounter"), 11.5); // AllRejectedTrigOverlaps + return; + } + registry.fill(HIST("hDetTrigcollisionCounter"), 12.5); // EMCAcceptedDetTrigColl + // if (!eventAccepted) { + // // for (auto const& jet : jets) { + // // if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax) || !isAcceptedRecoJet(jet)) { + // // fillRejectedJetHistograms(jet, 1.0); + // // } + // // } + // registry.fill(HIST("hDetTrigcollisionCounter"), 4.5); // AllRejectedDetTrigEventsAfterEMCEventSelection + // return; + // } + // registry.fill(HIST("hDetTrigcollisionCounter"), 5.5); // EMCAcceptedDetTrigCollWithkTVXinEMC + // -//*****Step 3: Reject ALL overlapping events by applying EXCLUSIVE Trigger Selections ***** -// Skip further processing if ANY overlaps exist -if ((hasFullJetHighPt && (hasFullJetLowPt || hasMB)) || (hasFullJetLowPt && hasMB)) { - registry.fill(HIST("hDetTrigcollisionCounter"), 11.5); //AllRejectedTrigOverlaps - return; -} -registry.fill(HIST("hDetTrigcollisionCounter"), 12.5); //EMCAcceptedDetTrigColl -// if (!eventAccepted) { -// // for (auto const& jet : jets) { -// // if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax) || !isAcceptedRecoJet(jet)) { -// // fillRejectedJetHistograms(jet, 1.0); -// // } -// // } -// registry.fill(HIST("hDetTrigcollisionCounter"), 4.5); // AllRejectedDetTrigEventsAfterEMCEventSelection -// return; -// } -// registry.fill(HIST("hDetTrigcollisionCounter"), 5.5); // EMCAcceptedDetTrigCollWithkTVXinEMC -// + // if (jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::JetChLowPt)) { + // registry.fill(HIST("hDetTrigcollisionCounter"), 8.5); // EMCAcceptedDetTrigCollWithLowChargedJetTriggers + // eventAccepted = true; + // } + // if (jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::JetChHighPt)) { + // registry.fill(HIST("hDetTrigcollisionCounter"), 9.5); // EMCAcceptedDetTrigCollWithHighChargedJetTriggers + // eventAccepted = true; + // } -// if (jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::JetChLowPt)) { -// registry.fill(HIST("hDetTrigcollisionCounter"), 8.5); // EMCAcceptedDetTrigCollWithLowChargedJetTriggers -// eventAccepted = true; -// } -// if (jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::JetChHighPt)) { -// registry.fill(HIST("hDetTrigcollisionCounter"), 9.5); // EMCAcceptedDetTrigCollWithHighChargedJetTriggers -// eventAccepted = true; -// } - -// if (jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::JetFullLowPt) && jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::JetFullHighPt)) { -// registry.fill(HIST("hDetTrigcollisionCounter"), 8.5); // EMCAcceptedDetTrigCollWithLow+HighFullJetTriggers -// } -// for (auto const& jet : jets) { -// if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { -// continue; -// } -// if (jet.phi() < jetPhiMin || jet.phi() > jetPhiMax) { -// continue; -// } -// if (!isAcceptedRecoJet(jet)) { -// continue; -// } -// fillJetHistograms(jet); -// } -} -PROCESS_SWITCH(FullJetSpectra, processJetsTriggeredData, "Full Jets Triggered Data", false); + // if (jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::JetFullLowPt) && jetderiveddatautilities::selectTrigger(collision, jetderiveddatautilities::JTrigSel::JetFullHighPt)) { + // registry.fill(HIST("hDetTrigcollisionCounter"), 8.5); // EMCAcceptedDetTrigCollWithLow+HighFullJetTriggers + // } + // for (auto const& jet : jets) { + // if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { + // continue; + // } + // if (jet.phi() < jetPhiMin || jet.phi() > jetPhiMax) { + // continue; + // } + // if (!isAcceptedRecoJet(jet)) { + // continue; + // } + // fillJetHistograms(jet); + // } + } + PROCESS_SWITCH(FullJetSpectra, processJetsTriggeredData, "Full Jets Triggered Data", false); + void processJetsMCD(soa::Filtered::iterator const& collision, JetTableMCDJoined const& jets, aod::JetTracks const&, aod::JetClusters const&) + { + bool eventAccepted = false; + double weight = 1.0; + float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent)); -void processJetsMCD(soa::Filtered::iterator const& collision, JetTableMCDJoined const& jets, aod::JetTracks const&, aod::JetClusters const&) -{ - bool eventAccepted = false; - double weight = 1.0; - float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent)); + /* if (doMcClosure) { + // Count total events processed (before splitting decision) + registry.fill(HIST("h_MCD_total_events"), 0.5); - /* if (doMcClosure) { - // Count total events processed (before splitting decision) - registry.fill(HIST("h_MCD_total_events"), 0.5); + // DEBUG: Let's verify what collision IDs we're actually seeing + LOGF(info, "[MCD DEBUG] Processing MC collision ID: %lld", collision.mcCollisionId()); - // DEBUG: Let's verify what collision IDs we're actually seeing - LOGF(info, "[MCD DEBUG] Processing MC collision ID: %lld", collision.mcCollisionId()); + // Get random value for this MC collision + float eventRandomValue = getMCCollisionRandomValue(collision.mcCollisionId()); - // Get random value for this MC collision - float eventRandomValue = getMCCollisionRandomValue(collision.mcCollisionId()); + // MCD gets events with random value <= split fraction (20%) + if (eventRandomValue > mcClosureSplitFrac) { + LOGF(debug, "[MCD] Event REJECTED: rand = %.4f > split = %.2f (MC collision %d)", + eventRandomValue, static_cast(mcClosureSplitFrac), collision.mcCollisionId()); + return; // This event goes to MCP & Matched processes + } - // MCD gets events with random value <= split fraction (20%) - if (eventRandomValue > mcClosureSplitFrac) { - LOGF(debug, "[MCD] Event REJECTED: rand = %.4f > split = %.2f (MC collision %d)", + LOGF(info, "[MCD] Event ACCEPTED: rand = %.4f <= split = %.2f (MC collision %d)", eventRandomValue, static_cast(mcClosureSplitFrac), collision.mcCollisionId()); - return; // This event goes to MCP & Matched processes -} -LOGF(info, "[MCD] Event ACCEPTED: rand = %.4f <= split = %.2f (MC collision %d)", -eventRandomValue, static_cast(mcClosureSplitFrac), collision.mcCollisionId()); + registry.fill(HIST("hSpliteventSelector"), 0.5); // 20% Closure input for the measured spectra (reco) + registry.fill(HIST("h_MCD_splitevent_counter"), 0.5); + } + */ + registry.fill(HIST("hDetcollisionCounter"), 0.5); // allDetColl + if (std::fabs(collision.posZ()) > vertexZCut) { + return; + } + registry.fill(HIST("hDetcollisionCounter"), 1.5); // DetCollWithVertexZ -registry.fill(HIST("hSpliteventSelector"), 0.5); // 20% Closure input for the measured spectra (reco) -registry.fill(HIST("h_MCD_splitevent_counter"), 0.5); -} -*/ -registry.fill(HIST("hDetcollisionCounter"), 0.5); // allDetColl -if (std::fabs(collision.posZ()) > vertexZCut) { - return; -} -registry.fill(HIST("hDetcollisionCounter"), 1.5); // DetCollWithVertexZ + // outlier check: for every outlier jet, reject the whole event + for (auto const& jet : jets) { + if (jet.pt() > pTHatMaxMCD * pTHat || pTHat < pTHatAbsoluteMin) { // for MCD jets only to remove outliers; setting pTHatMaxMCD = 1 improves purity + registry.fill(HIST("hDetcollisionCounter"), 2.5); // RejectedDetCollWithOutliers + return; + } + // this cut only to be used for calculating Jet Purity and not for Response Matrix + // this is mainly applied to remove all high weight jets causing big fluctuations + if (jet.pt() > 1 * pTHat) { + registry.fill(HIST("h_full_jet_pt_pTHatcut"), jet.pt(), weight); + } + } + if (doMBGapTrigger && collision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { + registry.fill(HIST("hDetcollisionCounter"), 3.5); // MBRejectedDetEvents + return; + } + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { + registry.fill(HIST("hDetcollisionCounter"), 4.5); // EventsNotSatisfyingEventSelection + return; + } + if (doEMCALEventWorkaround) { + if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content + if (collision.alias_bit(kTVXinEMC)) { + eventAccepted = true; + registry.fill(HIST("hDetcollisionCounter"), 5.5); // EMCreadoutDetEventsWithkTVXinEMC + } + } + } else { + if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { + eventAccepted = true; + registry.fill(HIST("hDetcollisionCounter"), 5.5); // EMCreadoutDetEventsWithkTVXinEMC + } + } -// outlier check: for every outlier jet, reject the whole event -for (auto const& jet : jets) { - if (jet.pt() > pTHatMaxMCD * pTHat || pTHat < pTHatAbsoluteMin) { // for MCD jets only to remove outliers; setting pTHatMaxMCD = 1 improves purity - registry.fill(HIST("hDetcollisionCounter"), 2.5); // RejectedDetCollWithOutliers - return; - } - // this cut only to be used for calculating Jet Purity and not for Response Matrix - // this is mainly applied to remove all high weight jets causing big fluctuations - if (jet.pt() > 1 * pTHat) { - registry.fill(HIST("h_full_jet_pt_pTHatcut"), jet.pt(), weight); - } -} -if (doMBGapTrigger && collision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { - registry.fill(HIST("hDetcollisionCounter"), 3.5); // MBRejectedDetEvents - return; -} -if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { - registry.fill(HIST("hDetcollisionCounter"), 4.5); // EventsNotSatisfyingEventSelection - return; -} -if (doEMCALEventWorkaround) { - if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content - if (collision.alias_bit(kTVXinEMC)) { - eventAccepted = true; - registry.fill(HIST("hDetcollisionCounter"), 5.5); // EMCreadoutDetEventsWithkTVXinEMC + if (!eventAccepted) { + for (auto const& jet : jets) { + if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax) || !isAcceptedRecoJet(jet)) { + fillRejectedJetHistograms(jet, 1.0); + } + } + registry.fill(HIST("hDetcollisionCounter"), 6.5); // AllRejectedDetEventsAfterEMCEventSelection + return; } - } -} else { - if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { - eventAccepted = true; - registry.fill(HIST("hDetcollisionCounter"), 5.5); // EMCreadoutDetEventsWithkTVXinEMC - } -} + registry.fill(HIST("hDetcollisionCounter"), 7.5); // EMCAcceptedDetColl -if (!eventAccepted) { - for (auto const& jet : jets) { - if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax) || !isAcceptedRecoJet(jet)) { - fillRejectedJetHistograms(jet, 1.0); + for (auto const& jet : jets) { + if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { + continue; + } + if (jet.phi() < jetPhiMin || jet.phi() > jetPhiMax) { + continue; + } + if (!isAcceptedRecoJet(jet)) { + continue; + } + fillJetHistograms(jet); } } - registry.fill(HIST("hDetcollisionCounter"), 6.5); // AllRejectedDetEventsAfterEMCEventSelection - return; -} -registry.fill(HIST("hDetcollisionCounter"), 7.5); // EMCAcceptedDetColl + PROCESS_SWITCH(FullJetSpectra, processJetsMCD, "Full Jets at Detector Level", false); -for (auto const& jet : jets) { - if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { - continue; - } - if (jet.phi() < jetPhiMin || jet.phi() > jetPhiMax) { - continue; - } - if (!isAcceptedRecoJet(jet)) { - continue; - } - fillJetHistograms(jet); -} -} -PROCESS_SWITCH(FullJetSpectra, processJetsMCD, "Full Jets at Detector Level", false); + void processJetsMCDWeighted(soa::Filtered::iterator const& collision, JetTableMCDWeightedJoined const& jets, aod::JMcCollisions const&, aod::JetTracks const&, aod::JetClusters const&) + { + bool eventAccepted = false; + double pTHat = 10. / (std::pow(collision.mcCollision().weight(), 1.0 / pTHatExponent)); -void processJetsMCDWeighted(soa::Filtered::iterator const& collision, JetTableMCDWeightedJoined const& jets, aod::JMcCollisions const&, aod::JetTracks const&, aod::JetClusters const&) -{ - bool eventAccepted = false; - double pTHat = 10. / (std::pow(collision.mcCollision().weight(), 1.0 / pTHatExponent)); + /* if (doMcClosure) { + // Count total events processed (before splitting decision) + registry.fill(HIST("h_MCD_total_events"), 0.5); - /* if (doMcClosure) { - // Count total events processed (before splitting decision) - registry.fill(HIST("h_MCD_total_events"), 0.5); + // DEBUG: Let's verify what collision IDs we're actually seeing + LOGF(info, "[MCD DEBUG] Processing MC collision ID: %lld", collision.mcCollisionId()); - // DEBUG: Let's verify what collision IDs we're actually seeing - LOGF(info, "[MCD DEBUG] Processing MC collision ID: %lld", collision.mcCollisionId()); + // Get random value for this MC collision + float eventRandomValue = getMCCollisionRandomValue(collision.mcCollisionId()); - // Get random value for this MC collision - float eventRandomValue = getMCCollisionRandomValue(collision.mcCollisionId()); + // MCD gets events with random value <= split fraction (20%) + if (eventRandomValue > mcClosureSplitFrac) { + LOGF(debug, "[MCD] Event REJECTED: rand = %.4f > split = %.2f (MC collision %d)", + eventRandomValue, static_cast(mcClosureSplitFrac), collision.mcCollisionId()); + return; // This event goes to MCP & Matched processes + } - // MCD gets events with random value <= split fraction (20%) - if (eventRandomValue > mcClosureSplitFrac) { - LOGF(debug, "[MCD] Event REJECTED: rand = %.4f > split = %.2f (MC collision %d)", + LOGF(info, "[MCD] Event ACCEPTED: rand = %.4f <= split = %.2f (MC collision %d)", eventRandomValue, static_cast(mcClosureSplitFrac), collision.mcCollisionId()); - return; // This event goes to MCP & Matched processes -} - -LOGF(info, "[MCD] Event ACCEPTED: rand = %.4f <= split = %.2f (MC collision %d)", -eventRandomValue, static_cast(mcClosureSplitFrac), collision.mcCollisionId()); -registry.fill(HIST("hSpliteventSelector"), 0.5); // 20% Closure input for the measured spectra (reco) -registry.fill(HIST("h_MCD_splitevent_counter"), 0.5); -} -*/ -registry.fill(HIST("hDetcollisionCounter"), 0.5, collision.mcCollision().weight()); // allDetColl -if (std::fabs(collision.posZ()) > vertexZCut) { - return; -} -registry.fill(HIST("hDetcollisionCounter"), 1.5, collision.mcCollision().weight()); // DetCollWithVertexZ -// outlier check: for every outlier jet, reject the whole event -for (auto const& jet : jets) { - if (jet.pt() > pTHatMaxMCD * pTHat || pTHat < pTHatAbsoluteMin) { // for MCD jets only to remove outliers; setting pTHatMaxMCD = 1 improves purity - registry.fill(HIST("hDetcollisionCounter"), 2.5, collision.mcCollision().weight()); // RejectedDetCollWithOutliers - return; - } - // this cut only to be used for calculating Jet Purity and not for Response Matrix - // this is mainly applied to remove all high weight jets causing big fluctuations - if (jet.pt() > 1 * pTHat) { - registry.fill(HIST("h_full_jet_pt_pTHatcut"), jet.pt(), collision.mcCollision().weight()); + registry.fill(HIST("hSpliteventSelector"), 0.5); // 20% Closure input for the measured spectra (reco) + registry.fill(HIST("h_MCD_splitevent_counter"), 0.5); } -} -if (doMBGapTrigger && collision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { - registry.fill(HIST("hDetcollisionCounter"), 3.5, collision.mcCollision().weight()); // MBRejectedDetEvents - return; -} -if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { - registry.fill(HIST("hDetcollisionCounter"), 4.5, collision.mcCollision().weight()); // EventsNotSatisfyingEventSelection - return; -} -if (doEMCALEventWorkaround) { - if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content - if (collision.alias_bit(kTVXinEMC)) { - eventAccepted = true; - registry.fill(HIST("hDetcollisionCounter"), 5.5, collision.mcCollision().weight()); // EMCreadoutDetEventsWithkTVXinEMC + */ + registry.fill(HIST("hDetcollisionCounter"), 0.5, collision.mcCollision().weight()); // allDetColl + if (std::fabs(collision.posZ()) > vertexZCut) { + return; + } + registry.fill(HIST("hDetcollisionCounter"), 1.5, collision.mcCollision().weight()); // DetCollWithVertexZ + // outlier check: for every outlier jet, reject the whole event + for (auto const& jet : jets) { + if (jet.pt() > pTHatMaxMCD * pTHat || pTHat < pTHatAbsoluteMin) { // for MCD jets only to remove outliers; setting pTHatMaxMCD = 1 improves purity + registry.fill(HIST("hDetcollisionCounter"), 2.5, collision.mcCollision().weight()); // RejectedDetCollWithOutliers + return; + } + // this cut only to be used for calculating Jet Purity and not for Response Matrix + // this is mainly applied to remove all high weight jets causing big fluctuations + if (jet.pt() > 1 * pTHat) { + registry.fill(HIST("h_full_jet_pt_pTHatcut"), jet.pt(), collision.mcCollision().weight()); + } + } + if (doMBGapTrigger && collision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { + registry.fill(HIST("hDetcollisionCounter"), 3.5, collision.mcCollision().weight()); // MBRejectedDetEvents + return; + } + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { + registry.fill(HIST("hDetcollisionCounter"), 4.5, collision.mcCollision().weight()); // EventsNotSatisfyingEventSelection + return; + } + if (doEMCALEventWorkaround) { + if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content + if (collision.alias_bit(kTVXinEMC)) { + eventAccepted = true; + registry.fill(HIST("hDetcollisionCounter"), 5.5, collision.mcCollision().weight()); // EMCreadoutDetEventsWithkTVXinEMC + } + } + } else { + if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { + eventAccepted = true; + registry.fill(HIST("hDetcollisionCounter"), 5.5, collision.mcCollision().weight()); // EMCreadoutDetEventsWithkTVXinEMC + } } - } -} else { - if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { - eventAccepted = true; - registry.fill(HIST("hDetcollisionCounter"), 5.5, collision.mcCollision().weight()); // EMCreadoutDetEventsWithkTVXinEMC - } -} -if (!eventAccepted) { - for (auto const& jet : jets) { - if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax) || !isAcceptedRecoJet(jet)) { - fillRejectedJetHistograms(jet, collision.mcCollision().weight()); + if (!eventAccepted) { + for (auto const& jet : jets) { + if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax) || !isAcceptedRecoJet(jet)) { + fillRejectedJetHistograms(jet, collision.mcCollision().weight()); + } + } + registry.fill(HIST("hDetcollisionCounter"), 6.5, collision.mcCollision().weight()); // AllRejectedDetEventsAfterEMCEventSelection + return; } - } - registry.fill(HIST("hDetcollisionCounter"), 6.5, collision.mcCollision().weight()); // AllRejectedDetEventsAfterEMCEventSelection - return; -} -registry.fill(HIST("hDetcollisionCounter"), 7.5, collision.mcCollision().weight()); // EMCAcceptedDetColl + registry.fill(HIST("hDetcollisionCounter"), 7.5, collision.mcCollision().weight()); // EMCAcceptedDetColl -for (auto const& jet : jets) { - if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { - continue; - } - if (jet.phi() < jetPhiMin || jet.phi() > jetPhiMax) { - continue; - } - if (!isAcceptedRecoJet(jet)) { - continue; + for (auto const& jet : jets) { + if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { + continue; + } + if (jet.phi() < jetPhiMin || jet.phi() > jetPhiMax) { + continue; + } + if (!isAcceptedRecoJet(jet)) { + continue; + } + fillJetHistograms(jet, collision.mcCollision().weight()); + } } - fillJetHistograms(jet, collision.mcCollision().weight()); -} -} -PROCESS_SWITCH(FullJetSpectra, processJetsMCDWeighted, "Full Jets at Detector Level on weighted events", false); - -void processJetsMCP(aod::JetMcCollision const& mccollision, JetTableMCPJoined const& jets, aod::JetParticles const&, soa::SmallGroups const& collisions) -{ - bool eventAccepted = false; - double weight = 1.0; - float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent)); + PROCESS_SWITCH(FullJetSpectra, processJetsMCDWeighted, "Full Jets at Detector Level on weighted events", false); - /* if (doMcClosure) { - // Count total events processed (before splitting decision) - registry.fill(HIST("h_MCP_total_events"), 0.5); + void processJetsMCP(aod::JetMcCollision const& mccollision, JetTableMCPJoined const& jets, aod::JetParticles const&, soa::SmallGroups const& collisions) + { + bool eventAccepted = false; + double weight = 1.0; + float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent)); - // DEBUG: Let's verify what collision IDs we're actually seeing - LOGF(info, "[MCP DEBUG] Processing MC collision ID: %lld", mccollision.globalIndex()); + /* if (doMcClosure) { + // Count total events processed (before splitting decision) + registry.fill(HIST("h_MCP_total_events"), 0.5); - // Get random value for this MC collision - float eventRandomValue = getMCCollisionRandomValue(mccollision.globalIndex()); + // DEBUG: Let's verify what collision IDs we're actually seeing + LOGF(info, "[MCP DEBUG] Processing MC collision ID: %lld", mccollision.globalIndex()); - // DEBUG: Track which MC collisions we're processing - registry.fill(HIST("hMCCollisionIdDebug_MCP"), static_cast(mccollision.globalIndex() % 100000)); + // Get random value for this MC collision + float eventRandomValue = getMCCollisionRandomValue(mccollision.globalIndex()); - // MCP gets events with random value > split fraction (80%) - if (eventRandomValue <= mcClosureSplitFrac) { - LOGF(debug, "[MCP] Event REJECTED: rand = %.4f <= split = %.2f (MC collision %lld)", - eventRandomValue, static_cast(mcClosureSplitFrac), mccollision.globalIndex()); - return; // This event goes to MCD only -} + // DEBUG: Track which MC collisions we're processing + registry.fill(HIST("hMCCollisionIdDebug_MCP"), static_cast(mccollision.globalIndex() % 100000)); -LOGF(info, "[MCP] Event ACCEPTED: rand = %.4f > split = %.2f (MC collision %lld)", -eventRandomValue, static_cast(mcClosureSplitFrac), mccollision.globalIndex()); + // MCP gets events with random value > split fraction (80%) + if (eventRandomValue <= mcClosureSplitFrac) { + LOGF(debug, "[MCP] Event REJECTED: rand = %.4f <= split = %.2f (MC collision %lld)", + eventRandomValue, static_cast(mcClosureSplitFrac), mccollision.globalIndex()); + return; // This event goes to MCD only + } -registry.fill(HIST("hSpliteventSelector"), 1.5); // remaining 80% input for MCP -registry.fill(HIST("h_MCP_splitevent_counter"), 0.5); -} -*/ -registry.fill(HIST("hPartcollisionCounter"), 0.5); // allMcColl -if (std::fabs(mccollision.posZ()) > vertexZCut) { - return; -} -registry.fill(HIST("hPartcollisionCounter"), 1.5); // McCollWithVertexZ -// if (collisions.size() < 1) { -// return; -// } -// registry.fill(HIST("hPartcollisionCounter"), 2.5); // PartCollWithSize>1 -// -// if (collisions.size() == 0) { -// registry.fill(HIST("hPartcollisionCounter"), 3.5); // RejectedPartCollForDetCollWithSize0 -// return; -// } + LOGF(info, "[MCP] Event ACCEPTED: rand = %.4f > split = %.2f (MC collision %lld)", + eventRandomValue, static_cast(mcClosureSplitFrac), mccollision.globalIndex()); -// outlier check: for every outlier jet, reject the whole event -for (auto const& jet : jets) { - if (jet.pt() > pTHatMaxMCP * pTHat || pTHat < pTHatAbsoluteMin) { - registry.fill(HIST("hPartcollisionCounter"), 2.5); // RejectedPartCollWithOutliers - return; + registry.fill(HIST("hSpliteventSelector"), 1.5); // remaining 80% input for MCP + registry.fill(HIST("h_MCP_splitevent_counter"), 0.5); } -} + */ + registry.fill(HIST("hPartcollisionCounter"), 0.5); // allMcColl + if (std::fabs(mccollision.posZ()) > vertexZCut) { + return; + } + registry.fill(HIST("hPartcollisionCounter"), 1.5); // McCollWithVertexZ + // if (collisions.size() < 1) { + // return; + // } + // registry.fill(HIST("hPartcollisionCounter"), 2.5); // PartCollWithSize>1 + // + // if (collisions.size() == 0) { + // registry.fill(HIST("hPartcollisionCounter"), 3.5); // RejectedPartCollForDetCollWithSize0 + // return; + // } -if (doMBGapTrigger && mccollision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { - // Fill rejected MB events; - registry.fill(HIST("hPartcollisionCounter"), 3.5); // MBRejectedPartEvents - return; -} + // outlier check: for every outlier jet, reject the whole event + for (auto const& jet : jets) { + if (jet.pt() > pTHatMaxMCP * pTHat || pTHat < pTHatAbsoluteMin) { + registry.fill(HIST("hPartcollisionCounter"), 2.5); // RejectedPartCollWithOutliers + return; + } + } -auto collisionspermcpjet = collisions.sliceBy(CollisionsPerMCPCollision, mccollision.globalIndex()); -registry.fill(HIST("hRecoMatchesPerMcCollision"), collisionspermcpjet.size()); // for split vertices QA + if (doMBGapTrigger && mccollision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { + // Fill rejected MB events; + registry.fill(HIST("hPartcollisionCounter"), 3.5); // MBRejectedPartEvents + return; + } -if (collisionspermcpjet.size() == 0 || collisionspermcpjet.size() < 1) { - registry.fill(HIST("hPartcollisionCounter"), 4.5); // RejectedPartCollForDetCollWithSize0or<1 - return; -} -registry.fill(HIST("hPartcollisionCounter"), 5.5); // AcceptedPartCollWithSize>=1 + auto collisionspermcpjet = collisions.sliceBy(CollisionsPerMCPCollision, mccollision.globalIndex()); + registry.fill(HIST("hRecoMatchesPerMcCollision"), collisionspermcpjet.size()); // for split vertices QA -for (auto const& collision : collisionspermcpjet) { - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { - continue; - } - if (doEMCALEventWorkaround) { - if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content - if (collision.alias_bit(kTVXinEMC)) { - eventAccepted = true; - registry.fill(HIST("hPartcollisionCounter"), 6.5); // EMCreadoutDetEventsWithkTVXinEMC + if (collisionspermcpjet.size() == 0 || collisionspermcpjet.size() < 1) { + registry.fill(HIST("hPartcollisionCounter"), 4.5); // RejectedPartCollForDetCollWithSize0or<1 + return; + } + registry.fill(HIST("hPartcollisionCounter"), 5.5); // AcceptedPartCollWithSize>=1 + + for (auto const& collision : collisionspermcpjet) { + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { + continue; + } + if (doEMCALEventWorkaround) { + if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content + if (collision.alias_bit(kTVXinEMC)) { + eventAccepted = true; + registry.fill(HIST("hPartcollisionCounter"), 6.5); // EMCreadoutDetEventsWithkTVXinEMC + } + } + } else { + if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { + eventAccepted = true; + registry.fill(HIST("hPartcollisionCounter"), 6.5); // EMCreadoutDetEventsWithkTVXinEMC + } } } - } else { - if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { - eventAccepted = true; - registry.fill(HIST("hPartcollisionCounter"), 6.5); // EMCreadoutDetEventsWithkTVXinEMC + if (!eventAccepted) { + registry.fill(HIST("hPartcollisionCounter"), 7.5); // AllRejectedPartEventsAfterEMCEventSelection + return; } - } -} -if (!eventAccepted) { - registry.fill(HIST("hPartcollisionCounter"), 7.5); // AllRejectedPartEventsAfterEMCEventSelection - return; -} -registry.fill(HIST("hPartcollisionCounter"), 8.5); // EMCAcceptedPartColl + registry.fill(HIST("hPartcollisionCounter"), 8.5); // EMCAcceptedPartColl -for (auto const& jet : jets) { - if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { - continue; - } - if (!isAcceptedPartJet(jet)) { - continue; + for (auto const& jet : jets) { + if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { + continue; + } + if (!isAcceptedPartJet(jet)) { + continue; + } + fillMCPHistograms(jet); + } } - fillMCPHistograms(jet); -} -} -PROCESS_SWITCH(FullJetSpectra, processJetsMCP, "Full Jets at Particle Level", false); - -void processJetsMCPWeighted(aod::JetMcCollision const& mccollision, JetTableMCPWeightedJoined const& jets, aod::JetParticles const&, soa::SmallGroups const& collisions) -{ - bool eventAccepted = false; - float pTHat = 10. / (std::pow(mccollision.weight(), 1.0 / pTHatExponent)); + PROCESS_SWITCH(FullJetSpectra, processJetsMCP, "Full Jets at Particle Level", false); - /* if (doMcClosure) { - // Count total events processed (before splitting decision) - registry.fill(HIST("h_MCP_total_events"), 0.5); + void processJetsMCPWeighted(aod::JetMcCollision const& mccollision, JetTableMCPWeightedJoined const& jets, aod::JetParticles const&, soa::SmallGroups const& collisions) + { + bool eventAccepted = false; + float pTHat = 10. / (std::pow(mccollision.weight(), 1.0 / pTHatExponent)); - // DEBUG: Let's verify what collision IDs we're actually seeing - LOGF(info, "[MCP DEBUG] Processing MC collision ID: %lld", mccollision.globalIndex()); + /* if (doMcClosure) { + // Count total events processed (before splitting decision) + registry.fill(HIST("h_MCP_total_events"), 0.5); - // Get random value for this MC collision - float eventRandomValue = getMCCollisionRandomValue(mccollision.globalIndex()); + // DEBUG: Let's verify what collision IDs we're actually seeing + LOGF(info, "[MCP DEBUG] Processing MC collision ID: %lld", mccollision.globalIndex()); - // DEBUG: Track which MC collisions we're processing - registry.fill(HIST("hMCCollisionIdDebug_MCP"), static_cast(mccollision.globalIndex() % 100000)); + // Get random value for this MC collision + float eventRandomValue = getMCCollisionRandomValue(mccollision.globalIndex()); - // MCP gets events with random value > split fraction (80%) - if (eventRandomValue <= mcClosureSplitFrac) { - LOGF(debug, "[MCP] Event REJECTED: rand = %.4f <= split = %.2f (MC collision %lld)", - eventRandomValue, static_cast(mcClosureSplitFrac), mccollision.globalIndex()); - return; // This event goes to MCD only -} + // DEBUG: Track which MC collisions we're processing + registry.fill(HIST("hMCCollisionIdDebug_MCP"), static_cast(mccollision.globalIndex() % 100000)); -LOGF(info, "[MCP] Event ACCEPTED: rand = %.4f > split = %.2f (MC collision %lld)", -eventRandomValue, static_cast(mcClosureSplitFrac), mccollision.globalIndex()); + // MCP gets events with random value > split fraction (80%) + if (eventRandomValue <= mcClosureSplitFrac) { + LOGF(debug, "[MCP] Event REJECTED: rand = %.4f <= split = %.2f (MC collision %lld)", + eventRandomValue, static_cast(mcClosureSplitFrac), mccollision.globalIndex()); + return; // This event goes to MCD only + } -registry.fill(HIST("hSpliteventSelector"), 1.5); // remaining 80% input for MCP -registry.fill(HIST("h_MCP_splitevent_counter"), 0.5); -} -*/ -registry.fill(HIST("hPartcollisionCounter"), 0.5, mccollision.weight()); // allMcColl -if (std::fabs(mccollision.posZ()) > vertexZCut) { - return; -} -registry.fill(HIST("hPartcollisionCounter"), 1.5, mccollision.weight()); // McCollWithVertexZ -// if (collisions.size() < 1) { -// return; -// } -// registry.fill(HIST("hPartcollisionCounter"), 2.5, mccollision.weight()); // PartCollWithSize>1 -// -// if (collisions.size() == 0) { -// registry.fill(HIST("hPartcollisionCounter"), 3.5, mccollision.weight()); // RejectedPartCollForDetCollWithSize0 -// return; -// } + LOGF(info, "[MCP] Event ACCEPTED: rand = %.4f > split = %.2f (MC collision %lld)", + eventRandomValue, static_cast(mcClosureSplitFrac), mccollision.globalIndex()); -// outlier check: for every outlier jet, reject the whole event -for (auto const& jet : jets) { - if (jet.pt() > pTHatMaxMCP * pTHat || pTHat < pTHatAbsoluteMin) { - registry.fill(HIST("hPartcollisionCounter"), 2.5, mccollision.weight()); // RejectedPartCollWithOutliers - return; + registry.fill(HIST("hSpliteventSelector"), 1.5); // remaining 80% input for MCP + registry.fill(HIST("h_MCP_splitevent_counter"), 0.5); } -} + */ + registry.fill(HIST("hPartcollisionCounter"), 0.5, mccollision.weight()); // allMcColl + if (std::fabs(mccollision.posZ()) > vertexZCut) { + return; + } + registry.fill(HIST("hPartcollisionCounter"), 1.5, mccollision.weight()); // McCollWithVertexZ + // if (collisions.size() < 1) { + // return; + // } + // registry.fill(HIST("hPartcollisionCounter"), 2.5, mccollision.weight()); // PartCollWithSize>1 + // + // if (collisions.size() == 0) { + // registry.fill(HIST("hPartcollisionCounter"), 3.5, mccollision.weight()); // RejectedPartCollForDetCollWithSize0 + // return; + // } -if (doMBGapTrigger && mccollision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { - // Fill rejected MB events - registry.fill(HIST("hPartcollisionCounter"), 3.5, mccollision.weight()); // MBRejectedPartEvents - return; -} + // outlier check: for every outlier jet, reject the whole event + for (auto const& jet : jets) { + if (jet.pt() > pTHatMaxMCP * pTHat || pTHat < pTHatAbsoluteMin) { + registry.fill(HIST("hPartcollisionCounter"), 2.5, mccollision.weight()); // RejectedPartCollWithOutliers + return; + } + } -auto collisionspermcpjet = collisions.sliceBy(CollisionsPerMCPCollision, mccollision.globalIndex()); -registry.fill(HIST("hRecoMatchesPerMcCollision"), collisionspermcpjet.size(), mccollision.weight()); // for split vertices QA + if (doMBGapTrigger && mccollision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { + // Fill rejected MB events + registry.fill(HIST("hPartcollisionCounter"), 3.5, mccollision.weight()); // MBRejectedPartEvents + return; + } -if (collisionspermcpjet.size() == 0 || collisionspermcpjet.size() < 1) { - registry.fill(HIST("hPartcollisionCounter"), 4.5, mccollision.weight()); // RejectedPartCollForDetCollWithSize0or<1 - return; -} -registry.fill(HIST("hPartcollisionCounter"), 5.5, mccollision.weight()); // AcceptedPartCollWithSize>=1 + auto collisionspermcpjet = collisions.sliceBy(CollisionsPerMCPCollision, mccollision.globalIndex()); + registry.fill(HIST("hRecoMatchesPerMcCollision"), collisionspermcpjet.size(), mccollision.weight()); // for split vertices QA -for (auto const& collision : collisionspermcpjet) { - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { - continue; - } - if (doEMCALEventWorkaround) { - if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content - if (collision.alias_bit(kTVXinEMC)) { - eventAccepted = true; - registry.fill(HIST("hPartcollisionCounter"), 6.5, mccollision.weight()); // EMCreadoutDetEventsWithkTVXinEMC + if (collisionspermcpjet.size() == 0 || collisionspermcpjet.size() < 1) { + registry.fill(HIST("hPartcollisionCounter"), 4.5, mccollision.weight()); // RejectedPartCollForDetCollWithSize0or<1 + return; + } + registry.fill(HIST("hPartcollisionCounter"), 5.5, mccollision.weight()); // AcceptedPartCollWithSize>=1 + + for (auto const& collision : collisionspermcpjet) { + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { + continue; + } + if (doEMCALEventWorkaround) { + if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content + if (collision.alias_bit(kTVXinEMC)) { + eventAccepted = true; + registry.fill(HIST("hPartcollisionCounter"), 6.5, mccollision.weight()); // EMCreadoutDetEventsWithkTVXinEMC + } + } + } else { + if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { + eventAccepted = true; + registry.fill(HIST("hPartcollisionCounter"), 6.5, mccollision.weight()); // EMCreadoutDetEventsWithkTVXinEMC + } } } - } else { - if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { - eventAccepted = true; - registry.fill(HIST("hPartcollisionCounter"), 6.5, mccollision.weight()); // EMCreadoutDetEventsWithkTVXinEMC + if (!eventAccepted) { + registry.fill(HIST("hPartcollisionCounter"), 7.5, mccollision.weight()); // AllRejectedPartEventsAfterEMCEventSelection + return; } - } -} -if (!eventAccepted) { - registry.fill(HIST("hPartcollisionCounter"), 7.5, mccollision.weight()); // AllRejectedPartEventsAfterEMCEventSelection - return; -} -// Fill EMCAL JJ Part events -registry.fill(HIST("hPartcollisionCounter"), 8.5, mccollision.weight()); // EMCAcceptedWeightedPartColl + // Fill EMCAL JJ Part events + registry.fill(HIST("hPartcollisionCounter"), 8.5, mccollision.weight()); // EMCAcceptedWeightedPartColl -for (auto const& jet : jets) { - if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { - return; - } - if (!isAcceptedPartJet(jet)) { - return; - } - if (doMBGapTrigger && jet.eventWeight() == 1) { - return; + for (auto const& jet : jets) { + if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { + return; + } + if (!isAcceptedPartJet(jet)) { + return; + } + if (doMBGapTrigger && jet.eventWeight() == 1) { + return; + } + fillMCPHistograms(jet, jet.eventWeight()); + } } - fillMCPHistograms(jet, jet.eventWeight()); -} -} -PROCESS_SWITCH(FullJetSpectra, processJetsMCPWeighted, "Full Jets at Particle Level on weighted events", false); + PROCESS_SWITCH(FullJetSpectra, processJetsMCPWeighted, "Full Jets at Particle Level on weighted events", false); -void processJetsMCPMCDMatched(soa::Filtered::iterator const& collision, JetTableMCDMatchedJoined const& mcdjets, jetMcpPerMcCollision const& mcpjets, aod::JMcCollisions const&, aod::JetTracks const&, aod::JetClusters const&, aod::JetParticles const&) -{ - bool eventAccepted = false; - int fakeMcdJet = 0; - int fakeMcpJet = 0; - double weight = 1.0; - float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent)); - const auto mcpJetsPerMcCollision = mcpjets.sliceBy(JetMCPPerMcCollision, collision.mcCollisionId()); - - /* if (doMcClosure) { - // Count total events processed (before splitting decision) - registry.fill(HIST("h_Matched_total_events"), 0.5); - - // Use consistent MC collision ID - same as MCD - int64_t mcCollisionId = collision.mcCollisionId(); - - // DEBUG: Let's verify what collision IDs we're actually seeing - LOGF(info, "[Matched DEBUG] Processing MC collision ID: %lld", mcCollisionId); - float eventRandomValue = getMCCollisionRandomValue(mcCollisionId); - - // Matched gets events with random value > split fraction (80%) - same as MCP - if (eventRandomValue <= mcClosureSplitFrac) { - LOGF(debug, "[Matched] Event REJECTED: rand = %.4f <= split = %.2f (MC collision %lld)", - eventRandomValue, static_cast(mcClosureSplitFrac), mcCollisionId); - return; // This event goes to MCD only -} + void processJetsMCPMCDMatched(soa::Filtered::iterator const& collision, JetTableMCDMatchedJoined const& mcdjets, jetMcpPerMcCollision const& mcpjets, aod::JMcCollisions const&, aod::JetTracks const&, aod::JetClusters const&, aod::JetParticles const&) + { + bool eventAccepted = false; + int fakeMcdJet = 0; + int fakeMcpJet = 0; + double weight = 1.0; + float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent)); + const auto mcpJetsPerMcCollision = mcpjets.sliceBy(JetMCPPerMcCollision, collision.mcCollisionId()); -LOGF(info, "[Matched] Event ACCEPTED: rand = %.4f > split = %.2f (MC collision %lld)", -eventRandomValue, static_cast(mcClosureSplitFrac), mcCollisionId); + /* if (doMcClosure) { + // Count total events processed (before splitting decision) + registry.fill(HIST("h_Matched_total_events"), 0.5); -registry.fill(HIST("hSpliteventSelector"), 2.5); // Bin for Response Matrix -registry.fill(HIST("h_Matched_splitevent_counter"), 0.5); -} -*/ -registry.fill(HIST("hMatchedcollisionCounter"), 0.5); // allDetColl + // Use consistent MC collision ID - same as MCD + int64_t mcCollisionId = collision.mcCollisionId(); -if (std::fabs(collision.posZ()) > vertexZCut) { // making double sure this condition is satisfied - return; -} -registry.fill(HIST("hMatchedcollisionCounter"), 1.5); // DetCollWithVertexZ + // DEBUG: Let's verify what collision IDs we're actually seeing + LOGF(info, "[Matched DEBUG] Processing MC collision ID: %lld", mcCollisionId); + float eventRandomValue = getMCCollisionRandomValue(mcCollisionId); -// outlier check: for every outlier jet, reject the whole event -for (auto const& mcdjet : mcdjets) { - if (mcdjet.pt() > pTHatMaxMCD * pTHat || pTHat < pTHatAbsoluteMin) { - registry.fill(HIST("hMatchedcollisionCounter"), 2.5); // RejectedDetCollWithOutliers - return; + // Matched gets events with random value > split fraction (80%) - same as MCP + if (eventRandomValue <= mcClosureSplitFrac) { + LOGF(debug, "[Matched] Event REJECTED: rand = %.4f <= split = %.2f (MC collision %lld)", + eventRandomValue, static_cast(mcClosureSplitFrac), mcCollisionId); + return; // This event goes to MCD only } -} -// //outlier check for Part collisions: commenting out this for now otherwise this rejects all Det Colls -// for (auto const& mcpjet : mcpjets) { -// if (mcpjet.pt() > pTHatMaxMCP * pTHat || pTHat < pTHatAbsoluteMin) { -// registry.fill(HIST("hMatchedcollisionCounter"),3.5); //RejectedPartCollWithOutliers -// return; -// } -// } - -if (doMBGapTrigger && collision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { - registry.fill(HIST("hMatchedcollisionCounter"), 4.5); // EMCMBRejectedDetColl - return; -} -if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { - registry.fill(HIST("hMatchedcollisionCounter"), 5.5); // EventsNotSatisfyingEventSelection - return; -} -if (doEMCALEventWorkaround) { - if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content - if (collision.alias_bit(kTVXinEMC)) { - eventAccepted = true; - registry.fill(HIST("hMatchedcollisionCounter"), 6.5); // EMCreadoutDetEventsWithkTVXinEMC - } - } -} else { - if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { - eventAccepted = true; - registry.fill(HIST("hMatchedcollisionCounter"), 6.5); // EMCreadoutDetEventsWithkTVXinEMC - } -} -if (!eventAccepted) { - registry.fill(HIST("hMatchedcollisionCounter"), 7.5); // AllRejectedDetEventsAfterEMCEventSelection - return; -} -registry.fill(HIST("hMatchedcollisionCounter"), 8.5); // EMCAcceptedDetColl + LOGF(info, "[Matched] Event ACCEPTED: rand = %.4f > split = %.2f (MC collision %lld)", + eventRandomValue, static_cast(mcClosureSplitFrac), mcCollisionId); -for (const auto& mcdjet : mcdjets) { - if (!isAcceptedRecoJet(mcdjet)) { - continue; - } - if (!jetfindingutilities::isInEtaAcceptance(mcdjet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { - continue; - } - // Check if MCD jet is within the EMCAL fiducial region; if not then flag it as a fake jet - if (mcdjet.phi() < jetPhiMin || mcdjet.phi() > jetPhiMax || mcdjet.eta() < jetEtaMin || mcdjet.eta() > jetEtaMax) { - fakeMcdJet++; - registry.fill(HIST("h2_full_fakemcdjets"), mcdjet.pt(), fakeMcdJet, 1.0); - continue; + registry.fill(HIST("hSpliteventSelector"), 2.5); // Bin for Response Matrix + registry.fill(HIST("h_Matched_splitevent_counter"), 0.5); } + */ + registry.fill(HIST("hMatchedcollisionCounter"), 0.5); // allDetColl - for (const auto& mcpjet : mcdjet.template matchedJetGeo_as()) { - if (!isAcceptedPartJet(mcpjet)) { + if (std::fabs(collision.posZ()) > vertexZCut) { // making double sure this condition is satisfied return; } - // apply emcal fiducial cuts to the matched particle level jets - if (mcpjet.eta() > jetEtaMax || mcpjet.eta() < jetEtaMin || mcpjet.phi() > jetPhiMax || mcpjet.phi() < jetPhiMin) { - fakeMcpJet++; - registry.fill(HIST("h2_full_fakemcpjets"), mcpjet.pt(), fakeMcpJet, 1.0); - continue; + registry.fill(HIST("hMatchedcollisionCounter"), 1.5); // DetCollWithVertexZ + + // outlier check: for every outlier jet, reject the whole event + for (auto const& mcdjet : mcdjets) { + if (mcdjet.pt() > pTHatMaxMCD * pTHat || pTHat < pTHatAbsoluteMin) { + registry.fill(HIST("hMatchedcollisionCounter"), 2.5); // RejectedDetCollWithOutliers + return; + } } - } // mcpjet loop - fillMatchedHistograms(mcdjet); -} // mcdjet loop -} -PROCESS_SWITCH(FullJetSpectra, processJetsMCPMCDMatched, "Full Jet finder MCP matched to MCD", false); + // //outlier check for Part collisions: commenting out this for now otherwise this rejects all Det Colls + // for (auto const& mcpjet : mcpjets) { + // if (mcpjet.pt() > pTHatMaxMCP * pTHat || pTHat < pTHatAbsoluteMin) { + // registry.fill(HIST("hMatchedcollisionCounter"),3.5); //RejectedPartCollWithOutliers + // return; + // } + // } -void processJetsMCPMCDMatchedWeighted(soa::Filtered::iterator const& collision, JetTableMCDMatchedWeightedJoined const& mcdjets, JetTableMCPMatchedWeightedJoined const& mcpjets, aod::JMcCollisions const&, aod::JetTracks const&, aod::JetClusters const&, aod::JetParticles const&) -{ - bool eventAccepted = false; - int fakeMcdJet = 0; - int fakeMcpJet = 0; - int NPartJetFid = 0; - float eventWeight = collision.mcCollision().weight(); - float pTHat = 10. / (std::pow(eventWeight, 1.0 / pTHatExponent)); - const auto mcpJetsPerMcCollision = mcpjets.sliceBy(JetMCPPerMcCollision, collision.mcCollisionId()); - - /* if (doMcClosure) { - // Count total events processed (before splitting decision) - registry.fill(HIST("h_Matched_total_events"), 0.5); - - // Use consistent MC collision ID - same as MCD - int64_t mcCollisionId = collision.mcCollisionId(); - - // DEBUG: Let's verify what collision IDs we're actually seeing - LOGF(info, "[Matched DEBUG] Processing MC collision ID: %lld", mcCollisionId); - float eventRandomValue = getMCCollisionRandomValue(mcCollisionId); - - // Matched gets events with random value > split fraction (80%) - same as MCP - if (eventRandomValue <= mcClosureSplitFrac) { - LOGF(debug, "[Matched] Event REJECTED: rand = %.4f <= split = %.2f (MC collision %lld)", - eventRandomValue, static_cast(mcClosureSplitFrac), mcCollisionId); - return; // This event goes to MCD only -} + if (doMBGapTrigger && collision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { + registry.fill(HIST("hMatchedcollisionCounter"), 4.5); // EMCMBRejectedDetColl + return; + } -LOGF(info, "[Matched] Event ACCEPTED: rand = %.4f > split = %.2f (MC collision %lld)", -eventRandomValue, static_cast(mcClosureSplitFrac), mcCollisionId); + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { + registry.fill(HIST("hMatchedcollisionCounter"), 5.5); // EventsNotSatisfyingEventSelection + return; + } + if (doEMCALEventWorkaround) { + if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content + if (collision.alias_bit(kTVXinEMC)) { + eventAccepted = true; + registry.fill(HIST("hMatchedcollisionCounter"), 6.5); // EMCreadoutDetEventsWithkTVXinEMC + } + } + } else { + if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { + eventAccepted = true; + registry.fill(HIST("hMatchedcollisionCounter"), 6.5); // EMCreadoutDetEventsWithkTVXinEMC + } + } + if (!eventAccepted) { + registry.fill(HIST("hMatchedcollisionCounter"), 7.5); // AllRejectedDetEventsAfterEMCEventSelection + return; + } + registry.fill(HIST("hMatchedcollisionCounter"), 8.5); // EMCAcceptedDetColl -registry.fill(HIST("hSpliteventSelector"), 2.5); // Bin for Response Matrix -registry.fill(HIST("h_Matched_splitevent_counter"), 0.5); -} -*/ -registry.fill(HIST("hMatchedcollisionCounter"), 0.5, eventWeight); // allDetColl -if (std::fabs(collision.posZ()) > vertexZCut) { // making double sure this condition is satisfied - return; -} -registry.fill(HIST("hMatchedcollisionCounter"), 1.5, eventWeight); // DetCollWithVertexZ + for (const auto& mcdjet : mcdjets) { + if (!isAcceptedRecoJet(mcdjet)) { + continue; + } + if (!jetfindingutilities::isInEtaAcceptance(mcdjet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { + continue; + } + // Check if MCD jet is within the EMCAL fiducial region; if not then flag it as a fake jet + if (mcdjet.phi() < jetPhiMin || mcdjet.phi() > jetPhiMax || mcdjet.eta() < jetEtaMin || mcdjet.eta() > jetEtaMax) { + fakeMcdJet++; + registry.fill(HIST("h2_full_fakemcdjets"), mcdjet.pt(), fakeMcdJet, 1.0); + continue; + } -// outlier check: for every outlier jet, reject the whole event -for (auto const& mcdjet : mcdjets) { - if (mcdjet.pt() > pTHatMaxMCD * pTHat || pTHat < pTHatAbsoluteMin) { - registry.fill(HIST("hMatchedcollisionCounter"), 2.5, eventWeight); // RejectedDetCollWithOutliers - return; + for (const auto& mcpjet : mcdjet.template matchedJetGeo_as()) { + if (!isAcceptedPartJet(mcpjet)) { + return; + } + // apply emcal fiducial cuts to the matched particle level jets + if (mcpjet.eta() > jetEtaMax || mcpjet.eta() < jetEtaMin || mcpjet.phi() > jetPhiMax || mcpjet.phi() < jetPhiMin) { + fakeMcpJet++; + registry.fill(HIST("h2_full_fakemcpjets"), mcpjet.pt(), fakeMcpJet, 1.0); + continue; + } + } // mcpjet loop + fillMatchedHistograms(mcdjet); + } // mcdjet loop } -} -// outlier check for Part collisions: commenting out this for now otherwise this rejects all Det Colls -// for (auto const& mcpjet : mcpjets) { -// if (mcpjet.pt() > pTHatMaxMCP * pTHat || pTHat < pTHatAbsoluteMin) { -// registry.fill(HIST("hMatchedcollisionCounter"),3.5, eventWeight); //RejectedPartCollWithOutliers -// return; -// } -// } - -if (doMBGapTrigger && collision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { - registry.fill(HIST("hMatchedcollisionCounter"), 4.5, eventWeight); // EMCMBRejectedDetColl - return; -} + PROCESS_SWITCH(FullJetSpectra, processJetsMCPMCDMatched, "Full Jet finder MCP matched to MCD", false); -if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { - registry.fill(HIST("hMatchedcollisionCounter"), 5.5, eventWeight); // EventsNotSatisfyingEventSelection - return; -} + void processJetsMCPMCDMatchedWeighted(soa::Filtered::iterator const& collision, JetTableMCDMatchedWeightedJoined const& mcdjets, JetTableMCPMatchedWeightedJoined const& mcpjets, aod::JMcCollisions const&, aod::JetTracks const&, aod::JetClusters const&, aod::JetParticles const&) + { + bool eventAccepted = false; + int fakeMcdJet = 0; + int fakeMcpJet = 0; + int NPartJetFid = 0; + float eventWeight = collision.mcCollision().weight(); + float pTHat = 10. / (std::pow(eventWeight, 1.0 / pTHatExponent)); + const auto mcpJetsPerMcCollision = mcpjets.sliceBy(JetMCPPerMcCollision, collision.mcCollisionId()); -for (auto const& mcpjet : mcpJetsPerMcCollision) { - if (mcpjet.pt() > pTHatMaxMCP * pTHat) { // outlier rejection for MCP: Should I remove this cut as I'm already doing MC outlier rejection @L1071? - return; -} -} -if (doEMCALEventWorkaround) { - if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content - if (collision.alias_bit(kTVXinEMC)) { - eventAccepted = true; - registry.fill(HIST("hMatchedcollisionCounter"), 6.5, eventWeight); // EMCreadoutDetJJEventsWithkTVXinEMC - } - } -} else { - if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { - eventAccepted = true; - registry.fill(HIST("hMatchedcollisionCounter"), 6.5, eventWeight); // EMCreadoutDetJJEventsWithkTVXinEMC - } -} -if (!eventAccepted) { - registry.fill(HIST("hMatchedcollisionCounter"), 7.5, eventWeight); // AllRejectedDetEventsAfterEMCEventSelection - return; -} -registry.fill(HIST("hMatchedcollisionCounter"), 8.5, eventWeight); // EMCAcceptedDetColl + /* if (doMcClosure) { + // Count total events processed (before splitting decision) + registry.fill(HIST("h_Matched_total_events"), 0.5); -for (const auto& mcdjet : mcdjets) { - if (!isAcceptedRecoJet(mcdjet)) { - continue; - } - // Check if MCD jet is within the EMCAL fiducial region; if not then flag it as a fake jet - if (mcdjet.phi() < jetPhiMin || mcdjet.phi() > jetPhiMax || mcdjet.eta() < jetEtaMin || mcdjet.eta() > jetEtaMax) { - fakeMcdJet++; - registry.fill(HIST("h2_full_fakemcdjets"), mcdjet.pt(), fakeMcdJet, eventWeight); - continue; + // Use consistent MC collision ID - same as MCD + int64_t mcCollisionId = collision.mcCollisionId(); + + // DEBUG: Let's verify what collision IDs we're actually seeing + LOGF(info, "[Matched DEBUG] Processing MC collision ID: %lld", mcCollisionId); + float eventRandomValue = getMCCollisionRandomValue(mcCollisionId); + + // Matched gets events with random value > split fraction (80%) - same as MCP + if (eventRandomValue <= mcClosureSplitFrac) { + LOGF(debug, "[Matched] Event REJECTED: rand = %.4f <= split = %.2f (MC collision %lld)", + eventRandomValue, static_cast(mcClosureSplitFrac), mcCollisionId); + return; // This event goes to MCD only } - if (!jetfindingutilities::isInEtaAcceptance(mcdjet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { - continue; + + LOGF(info, "[Matched] Event ACCEPTED: rand = %.4f > split = %.2f (MC collision %lld)", + eventRandomValue, static_cast(mcClosureSplitFrac), mcCollisionId); + + registry.fill(HIST("hSpliteventSelector"), 2.5); // Bin for Response Matrix + registry.fill(HIST("h_Matched_splitevent_counter"), 0.5); } + */ + registry.fill(HIST("hMatchedcollisionCounter"), 0.5, eventWeight); // allDetColl + if (std::fabs(collision.posZ()) > vertexZCut) { // making double sure this condition is satisfied + return; + } + registry.fill(HIST("hMatchedcollisionCounter"), 1.5, eventWeight); // DetCollWithVertexZ - for (const auto& mcpjet : mcdjet.template matchedJetGeo_as()) { - if (!isAcceptedPartJet(mcpjet)) { + // outlier check: for every outlier jet, reject the whole event + for (auto const& mcdjet : mcdjets) { + if (mcdjet.pt() > pTHatMaxMCD * pTHat || pTHat < pTHatAbsoluteMin) { + registry.fill(HIST("hMatchedcollisionCounter"), 2.5, eventWeight); // RejectedDetCollWithOutliers + return; + } + } + // outlier check for Part collisions: commenting out this for now otherwise this rejects all Det Colls + // for (auto const& mcpjet : mcpjets) { + // if (mcpjet.pt() > pTHatMaxMCP * pTHat || pTHat < pTHatAbsoluteMin) { + // registry.fill(HIST("hMatchedcollisionCounter"),3.5, eventWeight); //RejectedPartCollWithOutliers + // return; + // } + // } + + if (doMBGapTrigger && collision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { + registry.fill(HIST("hMatchedcollisionCounter"), 4.5, eventWeight); // EMCMBRejectedDetColl return; } - // apply emcal fiducial cuts to the matched particle level jets - if the matched mcp jet lies outside of the EMCAL fiducial, flag it as a fake jet - if (mcpjet.eta() > jetEtaMax || mcpjet.eta() < jetEtaMin || mcpjet.phi() > jetPhiMax || mcpjet.phi() < jetPhiMin) { - fakeMcpJet++; - registry.fill(HIST("h2FullfakeMcpJets"), mcpjet.pt(), fakeMcpJet, eventWeight); - continue; + + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { + registry.fill(HIST("hMatchedcollisionCounter"), 5.5, eventWeight); // EventsNotSatisfyingEventSelection + return; + } + + for (auto const& mcpjet : mcpJetsPerMcCollision) { + if (mcpjet.pt() > pTHatMaxMCP * pTHat) { // outlier rejection for MCP: Should I remove this cut as I'm already doing MC outlier rejection @L1071? + return; + } + } + if (doEMCALEventWorkaround) { + if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content + if (collision.alias_bit(kTVXinEMC)) { + eventAccepted = true; + registry.fill(HIST("hMatchedcollisionCounter"), 6.5, eventWeight); // EMCreadoutDetJJEventsWithkTVXinEMC + } + } } else { - NPartJetFid++; - // // If both MCD-MCP matched jet pairs are within the EMCAL fiducial region, fill these histos - registry.fill(HIST("h2_full_matchedmcpjet_pt"), mcpjet.pt(), NPartJetFid, eventWeight); - registry.fill(HIST("h_full_matchedmcpjet_eta"), mcpjet.eta(), eventWeight); - registry.fill(HIST("h_full_matchedmcpjet_phi"), mcpjet.phi(), eventWeight); - } - } // mcpjet - fillMatchedHistograms(mcdjet, eventWeight); -} // mcdjet -} -PROCESS_SWITCH(FullJetSpectra, processJetsMCPMCDMatchedWeighted, "Full Jet finder MCP matched to MCD on weighted events", false); + if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { + eventAccepted = true; + registry.fill(HIST("hMatchedcollisionCounter"), 6.5, eventWeight); // EMCreadoutDetJJEventsWithkTVXinEMC + } + } + if (!eventAccepted) { + registry.fill(HIST("hMatchedcollisionCounter"), 7.5, eventWeight); // AllRejectedDetEventsAfterEMCEventSelection + return; + } + registry.fill(HIST("hMatchedcollisionCounter"), 8.5, eventWeight); // EMCAcceptedDetColl -// Periodic cleanup to prevent unbounded memory growth -/*void processCleanup(aod::Collision const&) { -static int callCount = 0; -callCount++; + for (const auto& mcdjet : mcdjets) { + if (!isAcceptedRecoJet(mcdjet)) { + continue; + } + // Check if MCD jet is within the EMCAL fiducial region; if not then flag it as a fake jet + if (mcdjet.phi() < jetPhiMin || mcdjet.phi() > jetPhiMax || mcdjet.eta() < jetEtaMin || mcdjet.eta() > jetEtaMax) { + fakeMcdJet++; + registry.fill(HIST("h2_full_fakemcdjets"), mcdjet.pt(), fakeMcdJet, eventWeight); + continue; + } + if (!jetfindingutilities::isInEtaAcceptance(mcdjet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) { + continue; + } -// Clean up cache every 50000 calls to prevent memory issues -if (doMcClosure && callCount % 50000 == 0 && mcCollisionRandomValues.size() > 20000) { -LOGF(info, "Cleaning up MC collision random values cache (size: %zu)", mcCollisionRandomValues.size()); -mcCollisionRandomValues.clear(); + for (const auto& mcpjet : mcdjet.template matchedJetGeo_as()) { + if (!isAcceptedPartJet(mcpjet)) { + return; + } + // apply emcal fiducial cuts to the matched particle level jets - if the matched mcp jet lies outside of the EMCAL fiducial, flag it as a fake jet + if (mcpjet.eta() > jetEtaMax || mcpjet.eta() < jetEtaMin || mcpjet.phi() > jetPhiMax || mcpjet.phi() < jetPhiMin) { + fakeMcpJet++; + registry.fill(HIST("h2FullfakeMcpJets"), mcpjet.pt(), fakeMcpJet, eventWeight); + continue; + } else { + NPartJetFid++; + // // If both MCD-MCP matched jet pairs are within the EMCAL fiducial region, fill these histos + registry.fill(HIST("h2_full_matchedmcpjet_pt"), mcpjet.pt(), NPartJetFid, eventWeight); + registry.fill(HIST("h_full_matchedmcpjet_eta"), mcpjet.eta(), eventWeight); + registry.fill(HIST("h_full_matchedmcpjet_phi"), mcpjet.phi(), eventWeight); + } + } // mcpjet + fillMatchedHistograms(mcdjet, eventWeight); + } // mcdjet + } + PROCESS_SWITCH(FullJetSpectra, processJetsMCPMCDMatchedWeighted, "Full Jet finder MCP matched to MCD on weighted events", false); -// IMPROVEMENT: Add logging to verify our split ratios -float mcdCount = registry.get(HIST("h_MCD_splitevent_counter"))->GetBinContent(1); -float mcpCount = registry.get(HIST("h_MCP_splitevent_counter"))->GetBinContent(1); -float matchedCount = registry.get(HIST("h_Matched_splitevent_counter"))->GetBinContent(1); + // Periodic cleanup to prevent unbounded memory growth + /*void processCleanup(aod::Collision const&) { + static int callCount = 0; + callCount++; -float totalEvents = mcdCount + mcpCount; // MCP and Matched should be the same, so don't double count -float actualSplitFrac = totalEvents > 0 ? mcdCount / totalEvents : 0.0f; + // Clean up cache every 50000 calls to prevent memory issues + if (doMcClosure && callCount % 50000 == 0 && mcCollisionRandomValues.size() > 20000) { + LOGF(info, "Cleaning up MC collision random values cache (size: %zu)", mcCollisionRandomValues.size()); + mcCollisionRandomValues.clear(); -LOGF(info, "Current split statistics: MCD=%.1f, MCP=%.1f, Matched=%.1f", mcdCount, mcpCount, matchedCount); -LOGF(info, "Actual split fraction: %.3f (target: %.3f)", actualSplitFrac, static_cast(mcClosureSplitFrac)); -} -} -PROCESS_SWITCH(FullJetSpectra, processCleanup, "Periodic cleanup", true); -*/ -void processDataTracks(soa::Filtered::iterator const& collision, soa::Filtered const& tracks, soa::Filtered const& clusters) -{ - bool eventAccepted = false; + // IMPROVEMENT: Add logging to verify our split ratios + float mcdCount = registry.get(HIST("h_MCD_splitevent_counter"))->GetBinContent(1); + float mcpCount = registry.get(HIST("h_MCP_splitevent_counter"))->GetBinContent(1); + float matchedCount = registry.get(HIST("h_Matched_splitevent_counter"))->GetBinContent(1); - registry.fill(HIST("hCollisionsUnweighted"), 0.5); // allDetColl - if (std::fabs(collision.posZ()) > vertexZCut) { - return; - } - registry.fill(HIST("hCollisionsUnweighted"), 1.5); // DetCollWithVertexZ + float totalEvents = mcdCount + mcpCount; // MCP and Matched should be the same, so don't double count + float actualSplitFrac = totalEvents > 0 ? mcdCount / totalEvents : 0.0f; - if (doMBGapTrigger && collision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { - registry.fill(HIST("hCollisionsUnweighted"), 2.5); // MBRejectedDetEvents - return; + LOGF(info, "Current split statistics: MCD=%.1f, MCP=%.1f, Matched=%.1f", mcdCount, mcpCount, matchedCount); + LOGF(info, "Actual split fraction: %.3f (target: %.3f)", actualSplitFrac, static_cast(mcClosureSplitFrac)); } - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { - registry.fill(HIST("hCollisionsUnweighted"), 3.5); // EventsNotSatisfyingEventSelection - return; } - // needed for the workaround to access EMCAL trigger bits. - This is needed for the MC productions in which the EMC trigger bits are missing. (MB MC LHC24f3, for ex.) - // It first requires for atleast a cell in EMCAL to have energy content. - // Once it finds a cell content, - // it then checks if the collision is not an ambiguous collision (i.e. it has to be a unique collision = no bunch pile up) - // If all of these conditions are satisfied then it checks for the required trigger bit in EMCAL. - // For LHC22o, since the EMCAL didn't have hardware triggers, one would only require MB trigger (kTVXinEMC) in the EMCAL. + PROCESS_SWITCH(FullJetSpectra, processCleanup, "Periodic cleanup", true); + */ + void processDataTracks(soa::Filtered::iterator const& collision, soa::Filtered const& tracks, soa::Filtered const& clusters) + { + bool eventAccepted = false; + + registry.fill(HIST("hCollisionsUnweighted"), 0.5); // allDetColl + if (std::fabs(collision.posZ()) > vertexZCut) { + return; + } + registry.fill(HIST("hCollisionsUnweighted"), 1.5); // DetCollWithVertexZ - if (doEMCALEventWorkaround) { - if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content - if (collision.alias_bit(kTVXinEMC)) { + if (doMBGapTrigger && collision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { + registry.fill(HIST("hCollisionsUnweighted"), 2.5); // MBRejectedDetEvents + return; + } + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { + registry.fill(HIST("hCollisionsUnweighted"), 3.5); // EventsNotSatisfyingEventSelection + return; + } + // needed for the workaround to access EMCAL trigger bits. - This is needed for the MC productions in which the EMC trigger bits are missing. (MB MC LHC24f3, for ex.) + // It first requires for atleast a cell in EMCAL to have energy content. + // Once it finds a cell content, + // it then checks if the collision is not an ambiguous collision (i.e. it has to be a unique collision = no bunch pile up) + // If all of these conditions are satisfied then it checks for the required trigger bit in EMCAL. + // For LHC22o, since the EMCAL didn't have hardware triggers, one would only require MB trigger (kTVXinEMC) in the EMCAL. + + if (doEMCALEventWorkaround) { + if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content + if (collision.alias_bit(kTVXinEMC)) { + eventAccepted = true; + registry.fill(HIST("hCollisionsUnweighted"), 4.5); // EMCreadoutDetEventsWithkTVXinEMC + } + } + } else { + // Check if EMCAL was readout with the MB trigger(kTVXinEMC) fired. If not then reject the event and exit the function. + // This is the default check for the simulations with proper trigger flags not requiring the above workaround. + if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { eventAccepted = true; registry.fill(HIST("hCollisionsUnweighted"), 4.5); // EMCreadoutDetEventsWithkTVXinEMC } } - } else { - // Check if EMCAL was readout with the MB trigger(kTVXinEMC) fired. If not then reject the event and exit the function. - // This is the default check for the simulations with proper trigger flags not requiring the above workaround. - if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { - eventAccepted = true; - registry.fill(HIST("hCollisionsUnweighted"), 4.5); // EMCreadoutDetEventsWithkTVXinEMC - } - } - if (!eventAccepted) { - registry.fill(HIST("hCollisionsUnweighted"), 5.5); // AllRejectedDetEventsAfterEMCEventSelection - return; - } - registry.fill(HIST("hCollisionsUnweighted"), 6.5); // EMCAcceptedDetColl + if (!eventAccepted) { + registry.fill(HIST("hCollisionsUnweighted"), 5.5); // AllRejectedDetEventsAfterEMCEventSelection + return; + } + registry.fill(HIST("hCollisionsUnweighted"), 6.5); // EMCAcceptedDetColl - for (auto const& track : tracks) { - if (!jetderiveddatautilities::selectTrack(track, trackSelection)) { - continue; + for (auto const& track : tracks) { + if (!jetderiveddatautilities::selectTrack(track, trackSelection)) { + continue; + } + // Fill Accepted events histos + fillTrackHistograms(tracks, clusters, 1.0); } - // Fill Accepted events histos - fillTrackHistograms(tracks, clusters, 1.0); + registry.fill(HIST("hCollisionsUnweighted"), 7.5); // EMCAcceptedCollAfterTrackSel } - registry.fill(HIST("hCollisionsUnweighted"), 7.5); // EMCAcceptedCollAfterTrackSel -} -PROCESS_SWITCH(FullJetSpectra, processDataTracks, "Full Jet tracks for Data", false); + PROCESS_SWITCH(FullJetSpectra, processDataTracks, "Full Jet tracks for Data", false); -void processMCTracks(soa::Filtered::iterator const& collision, soa::Filtered const& tracks, soa::Filtered const& clusters) -{ - bool eventAccepted = false; - double weight = 1.0; - float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent)); + void processMCTracks(soa::Filtered::iterator const& collision, soa::Filtered const& tracks, soa::Filtered const& clusters) + { + bool eventAccepted = false; + double weight = 1.0; + float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent)); - registry.fill(HIST("hCollisionsUnweighted"), 0.5); // allDetColl - if (std::fabs(collision.posZ()) > vertexZCut) { - return; - } - registry.fill(HIST("hCollisionsUnweighted"), 1.5); // DetCollWithVertexZ + registry.fill(HIST("hCollisionsUnweighted"), 0.5); // allDetColl + if (std::fabs(collision.posZ()) > vertexZCut) { + return; + } + registry.fill(HIST("hCollisionsUnweighted"), 1.5); // DetCollWithVertexZ - // for (auto const& track : tracks) { - if (pTHat < pTHatAbsoluteMin) { // Track outlier rejection: should this be for every track iteration or for every collision? - return; - } - // } - if (doMBGapTrigger && collision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { - registry.fill(HIST("hCollisionsUnweighted"), 2.5); // MBRejectedDetEvents - return; - } - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { - registry.fill(HIST("hCollisionsUnweighted"), 3.5); // EventsNotSatisfyingEventSelection - return; - } + // for (auto const& track : tracks) { + if (pTHat < pTHatAbsoluteMin) { // Track outlier rejection: should this be for every track iteration or for every collision? + return; + } + // } + if (doMBGapTrigger && collision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { + registry.fill(HIST("hCollisionsUnweighted"), 2.5); // MBRejectedDetEvents + return; + } + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { + registry.fill(HIST("hCollisionsUnweighted"), 3.5); // EventsNotSatisfyingEventSelection + return; + } - if (doEMCALEventWorkaround) { - if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content - if (collision.alias_bit(kTVXinEMC)) { + if (doEMCALEventWorkaround) { + if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content + if (collision.alias_bit(kTVXinEMC)) { + eventAccepted = true; + registry.fill(HIST("hCollisionsUnweighted"), 4.5); // EMCreadoutDetEventsWithkTVXinEMC + } + } + } else { + if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { eventAccepted = true; registry.fill(HIST("hCollisionsUnweighted"), 4.5); // EMCreadoutDetEventsWithkTVXinEMC } } - } else { - if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { - eventAccepted = true; - registry.fill(HIST("hCollisionsUnweighted"), 4.5); // EMCreadoutDetEventsWithkTVXinEMC - } - } - if (!eventAccepted) { - registry.fill(HIST("hCollisionsUnweighted"), 5.5); // AllRejectedDetEventsAfterEMCEventSelection - return; - } - registry.fill(HIST("hCollisionsUnweighted"), 6.5); // EMCAcceptedDetColl + if (!eventAccepted) { + registry.fill(HIST("hCollisionsUnweighted"), 5.5); // AllRejectedDetEventsAfterEMCEventSelection + return; + } + registry.fill(HIST("hCollisionsUnweighted"), 6.5); // EMCAcceptedDetColl - for (auto const& track : tracks) { - if (!jetderiveddatautilities::selectTrack(track, trackSelection)) { - continue; + for (auto const& track : tracks) { + if (!jetderiveddatautilities::selectTrack(track, trackSelection)) { + continue; + } + // Fill Accepted events histos + fillTrackHistograms(tracks, clusters, 1.0); } - // Fill Accepted events histos - fillTrackHistograms(tracks, clusters, 1.0); + registry.fill(HIST("hCollisionsUnweighted"), 7.5); // EMCAcceptedCollAfterTrackSel } - registry.fill(HIST("hCollisionsUnweighted"), 7.5); // EMCAcceptedCollAfterTrackSel -} -PROCESS_SWITCH(FullJetSpectra, processMCTracks, "Full Jet tracks for MC", false); + PROCESS_SWITCH(FullJetSpectra, processMCTracks, "Full Jet tracks for MC", false); -void processTracksWeighted(soa::Filtered::iterator const& collision, - aod::JMcCollisions const&, - soa::Filtered const& tracks, - soa::Filtered const& clusters) + void processTracksWeighted(soa::Filtered::iterator const& collision, + aod::JMcCollisions const&, + soa::Filtered const& tracks, + soa::Filtered const& clusters) { bool eventAccepted = false; float eventWeight = collision.mcCollision().weight(); @@ -2103,13 +2100,16 @@ void processTracksWeighted(soa::Filtered::iterator const& coll for (auto const& jet : jets) { if (jet.collisionId() != collision.globalIndex()) { std::cout << "WARNING: Jet with pT " << jet.pt() << " belongs to collision " << jet.collisionId() - << " but processing collision " << collision.globalIndex() << std::endl; + << " but processing collision " << collision.globalIndex() << std::endl; continue; } - if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) continue; - if (jet.phi() < jetPhiMin || jet.phi() > jetPhiMax) continue; - if (!isAcceptedRecoJet(jet)) continue; + if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) + continue; + if (jet.phi() < jetPhiMin || jet.phi() > jetPhiMax) + continue; + if (!isAcceptedRecoJet(jet)) + continue; selectedJets.push_back(jet); nJetsThisEvent++; @@ -2117,7 +2117,7 @@ void processTracksWeighted(soa::Filtered::iterator const& coll } // 1. Sort selected jets by pT before processing std::sort(selectedJets.begin(), selectedJets.end(), - [](auto const& a, auto const& b) { return (*a).pt() > (*b).pt(); }); + [](auto const& a, auto const& b) { return (*a).pt() > (*b).pt(); }); // std::cout << "Number of selected jets: " << selectedJets.size() << std::endl; @@ -2139,12 +2139,12 @@ void processTracksWeighted(soa::Filtered::iterator const& coll } if (selectedJets.size() > 0) { - //Jet multiplicity per event + // Jet multiplicity per event registry.fill(HIST("h_all_fulljet_Njets"), selectedJets.size(), 1.0); - //Select Leading Jet for N_ch calculation (for every leading jet that is found). There's always one leading jet per event! + // Select Leading Jet for N_ch calculation (for every leading jet that is found). There's always one leading jet per event! auto const& leadingJet = *selectedJets[0]; - auto const& leadingJetPt = leadingJet.pt(); //jet pT distribution of the leading jet + auto const& leadingJetPt = leadingJet.pt(); // jet pT distribution of the leading jet // std::cout << "Leading Jet pT: " << leadingJetPt << std::endl; registry.fill(HIST("h_Leading_full_jet_pt"), leadingJetPt, 1.0); registry.fill(HIST("h2_full_jet_leadingJetPt_vs_counts"), leadingJetPt, nJetsThisEvent, 1.0); @@ -2152,7 +2152,7 @@ void processTracksWeighted(soa::Filtered::iterator const& coll if (selectedJets.size() > 1) { auto const& subLeadingJet = *selectedJets[1]; - auto const& subLeadingJetPt = subLeadingJet.pt(); //jet pT distribution of the subleading jet i.e. 2nd leading jet + auto const& subLeadingJetPt = subLeadingJet.pt(); // jet pT distribution of the subleading jet i.e. 2nd leading jet registry.fill(HIST("h_SubLeading_full_jet_pt"), subLeadingJetPt, 1.0); registry.fill(HIST("h2_full_jet_subLeadingJetPt_vs_counts"), subLeadingJetPt, nJetsThisEvent, 1.0); } @@ -2161,14 +2161,15 @@ void processTracksWeighted(soa::Filtered::iterator const& coll auto const& jet = *selectedJets[i]; float jetPt = jet.pt(); bool isLeading = (i == 0); - bool isSubLeading = (i == 1 && selectedJets.size() > 1); //first sub-leading jet + bool isSubLeading = (i == 1 && selectedJets.size() > 1); // first sub-leading jet // Count charged particles(NCh) for this jet int numberOfChargedParticles = 0; for (const auto& jettrack : jet.tracks_as()) { if (jetderiveddatautilities::selectTrack(jettrack, trackSelection)) { numberOfChargedParticles++; - } else continue; + } else + continue; } // Calculate neutral energy fraction for this jet @@ -2196,7 +2197,7 @@ void processTracksWeighted(soa::Filtered::iterator const& coll registry.fill(HIST("h3_leading_fulljet_jetpTDet_FT0Mults_nef"), jetPt, collision.multFT0M(), nef, 1.0); } - //CASE 3: Additional first sub-leading jet processing + // CASE 3: Additional first sub-leading jet processing if (isSubLeading) { registry.fill(HIST("h_subleading_fulljet_pt"), jetPt, 1.0); registry.fill(HIST("h_subleading_fulljet_Nch"), numberOfChargedParticles, 1.0); @@ -2258,13 +2259,16 @@ void processTracksWeighted(soa::Filtered::iterator const& coll for (auto const& mcdjet : mcdjets) { if (mcdjet.collisionId() != collision.globalIndex()) { std::cout << "WARNING: Jet with pT " << mcdjet.pt() << " belongs to collision " << mcdjet.collisionId() - << " but processing collision " << collision.globalIndex() << std::endl; + << " but processing collision " << collision.globalIndex() << std::endl; continue; } - if (!jetfindingutilities::isInEtaAcceptance(mcdjet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) continue; - if (mcdjet.phi() < jetPhiMin || mcdjet.phi() > jetPhiMax) continue; - if (!isAcceptedRecoJet(mcdjet)) continue; + if (!jetfindingutilities::isInEtaAcceptance(mcdjet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) + continue; + if (mcdjet.phi() < jetPhiMin || mcdjet.phi() > jetPhiMax) + continue; + if (!isAcceptedRecoJet(mcdjet)) + continue; selectedJets.push_back(mcdjet); nJetsThisEvent++; @@ -2272,19 +2276,19 @@ void processTracksWeighted(soa::Filtered::iterator const& coll } // 1. Sort selected jets by pT before processing std::sort(selectedJets.begin(), selectedJets.end(), - [](auto const& a, auto const& b) { return (*a).pt() > (*b).pt(); }); + [](auto const& a, auto const& b) { return (*a).pt() > (*b).pt(); }); if (selectedJets.size() == 0) { // no jets = no leading jet return; } if (selectedJets.size() > 0) { - //Jet multiplicity per event + // Jet multiplicity per event registry.fill(HIST("h_all_fulljet_Njets"), selectedJets.size(), 1.0); - //Select Leading Jet for N_ch calculation (for every leading jet that is found). There's always one leading jet per event! + // Select Leading Jet for N_ch calculation (for every leading jet that is found). There's always one leading jet per event! auto const& leadingJet = *selectedJets[0]; - auto const& leadingJetPt = leadingJet.pt(); //jet pT distribution of the leading jet + auto const& leadingJetPt = leadingJet.pt(); // jet pT distribution of the leading jet // std::cout << "Leading Jet pT: " << leadingJetPt << std::endl; registry.fill(HIST("h_Leading_full_jet_pt"), leadingJetPt, 1.0); registry.fill(HIST("h2_full_jet_leadingJetPt_vs_counts"), leadingJetPt, nJetsThisEvent, 1.0); @@ -2292,7 +2296,7 @@ void processTracksWeighted(soa::Filtered::iterator const& coll if (selectedJets.size() > 1) { auto const& subLeadingJet = *selectedJets[1]; - auto const& subLeadingJetPt = subLeadingJet.pt(); //jet pT distribution of the subleading jet i.e. 2nd leading jet + auto const& subLeadingJetPt = subLeadingJet.pt(); // jet pT distribution of the subleading jet i.e. 2nd leading jet registry.fill(HIST("h_SubLeading_full_jet_pt"), subLeadingJetPt, 1.0); registry.fill(HIST("h2_full_jet_subLeadingJetPt_vs_counts"), subLeadingJetPt, nJetsThisEvent, 1.0); } @@ -2301,14 +2305,15 @@ void processTracksWeighted(soa::Filtered::iterator const& coll auto const& jet = *selectedJets[i]; float jetPt = jet.pt(); bool isLeading = (i == 0); - bool isSubLeading = (i == 1 && selectedJets.size() > 1); //first sub-leading jet + bool isSubLeading = (i == 1 && selectedJets.size() > 1); // first sub-leading jet // Count charged particles(NCh) for this jet int numberOfChargedParticles = 0; for (const auto& jettrack : jet.tracks_as()) { if (jetderiveddatautilities::selectTrack(jettrack, trackSelection)) { numberOfChargedParticles++; - } else continue; + } else + continue; } // Calculate neutral energy fraction for this jet @@ -2336,7 +2341,7 @@ void processTracksWeighted(soa::Filtered::iterator const& coll registry.fill(HIST("h3_leading_fulljet_jetpTDet_FT0Mults_nef"), jetPt, collision.multFT0M(), nef, 1.0); } - //CASE 3: Additional first sub-leading jet processing + // CASE 3: Additional first sub-leading jet processing if (isSubLeading) { registry.fill(HIST("h_subleading_fulljet_pt"), jetPt, 1.0); registry.fill(HIST("h_subleading_fulljet_Nch"), numberOfChargedParticles, 1.0); @@ -2392,7 +2397,7 @@ void processTracksWeighted(soa::Filtered::iterator const& coll registry.fill(HIST("hEventmultiplicityCounter"), 5.5, eventWeight); // AllRejectedWeightedDetEventsAfterEMCEventSelection return; } - registry.fill(HIST("hEventmultiplicityCounter"), 6.5, eventWeight); // EMCAcceptedWeightedDetColl + registry.fill(HIST("hEventmultiplicityCounter"), 6.5, eventWeight); // EMCAcceptedWeightedDetColl for (auto const& track : tracks) { if (!jetderiveddatautilities::selectTrack(track, trackSelection)) { continue; @@ -2410,7 +2415,7 @@ void processTracksWeighted(soa::Filtered::iterator const& coll if (mcdjet.collisionId() != collision.globalIndex()) { std::cout << "WARNING: Jet with pT " << mcdjet.pt() << " belongs to collision " << mcdjet.collisionId() - << " but processing collision " << collision.globalIndex() << std::endl; + << " but processing collision " << collision.globalIndex() << std::endl; continue; } @@ -2423,7 +2428,7 @@ void processTracksWeighted(soa::Filtered::iterator const& coll if (mcdjet.phi() < jetPhiMin || mcdjet.phi() > jetPhiMax) { continue; } - if (!isAcceptedRecoJet(mcdjet)) { + if (!isAcceptedRecoJet(mcdjet)) { continue; } selectedJets.push_back(mcdjet); @@ -2431,18 +2436,18 @@ void processTracksWeighted(soa::Filtered::iterator const& coll } // 1. Sort selected jets by pT before processing std::sort(selectedJets.begin(), selectedJets.end(), - [](auto const& a, auto const& b) { return (*a).pt() > (*b).pt(); }); + [](auto const& a, auto const& b) { return (*a).pt() > (*b).pt(); }); if (selectedJets.size() == 0) { // no jets = no leading jet return; } if (selectedJets.size() > 0) { - //Jet multiplicity per event + // Jet multiplicity per event registry.fill(HIST("h_all_fulljet_Njets"), selectedJets.size(), eventWeight); - //Select Leading Jet for N_ch calculation (for every leading jet that is found). There's always one leading jet per event! + // Select Leading Jet for N_ch calculation (for every leading jet that is found). There's always one leading jet per event! auto const& leadingJet = *selectedJets[0]; - auto const& leadingJetPt = leadingJet.pt(); //jet pT distribution of the leading jet + auto const& leadingJetPt = leadingJet.pt(); // jet pT distribution of the leading jet // std::cout << "Leading Jet pT: " << leadingJetPt << std::endl; registry.fill(HIST("h_Leading_full_jet_pt"), leadingJetPt, eventWeight); registry.fill(HIST("h2_full_jet_leadingJetPt_vs_counts"), leadingJetPt, nJetsThisEvent, eventWeight); @@ -2450,7 +2455,7 @@ void processTracksWeighted(soa::Filtered::iterator const& coll if (selectedJets.size() > 1) { auto const& subLeadingJet = *selectedJets[1]; - auto const& subLeadingJetPt = subLeadingJet.pt(); //jet pT distribution of the subleading jet i.e. 2nd leading jet + auto const& subLeadingJetPt = subLeadingJet.pt(); // jet pT distribution of the subleading jet i.e. 2nd leading jet registry.fill(HIST("h_SubLeading_full_jet_pt"), subLeadingJetPt, eventWeight); registry.fill(HIST("h2_full_jet_subLeadingJetPt_vs_counts"), subLeadingJetPt, nJetsThisEvent, eventWeight); } @@ -2459,14 +2464,15 @@ void processTracksWeighted(soa::Filtered::iterator const& coll auto const& jet = *selectedJets[i]; float jetPt = jet.pt(); bool isLeading = (i == 0); - bool isSubLeading = (i == 1 && selectedJets.size() > 1); //first sub-leading jet + bool isSubLeading = (i == 1 && selectedJets.size() > 1); // first sub-leading jet // Count charged particles(NCh) for this jet int numberOfChargedParticles = 0; for (const auto& jettrack : jet.tracks_as()) { if (jetderiveddatautilities::selectTrack(jettrack, trackSelection)) { numberOfChargedParticles++; - } else continue; + } else + continue; } // Calculate neutral energy fraction for this jet @@ -2494,7 +2500,7 @@ void processTracksWeighted(soa::Filtered::iterator const& coll registry.fill(HIST("h3_leading_fulljet_jetpTDet_FT0Mults_nef"), jetPt, collision.multFT0M(), nef, eventWeight); } - //CASE 3: Additional first sub-leading jet processing + // CASE 3: Additional first sub-leading jet processing if (isSubLeading) { registry.fill(HIST("h_subleading_fulljet_pt"), jetPt, eventWeight); registry.fill(HIST("h_subleading_fulljet_Nch"), numberOfChargedParticles, eventWeight); @@ -2508,331 +2514,340 @@ void processTracksWeighted(soa::Filtered::iterator const& coll PROCESS_SWITCH(FullJetSpectra, processMCDCollisionsWeightedWithMultiplicity, "Weighted MCD Collisions for Full Jets Multiplicity Studies", false); void processMBMCPCollisionsWithMultiplicity(aod::JetMcCollision const& mccollision, - JetTableMCPJoined const& jets, aod::JetParticles const& /*particles*/, - soa::SmallGroups const& collisions) - { - bool eventAccepted = false; - double weight = 1.0; - float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent)); - float mcpMult = 0.0f; float mcpCent= 0.0f; - - registry.fill(HIST("hPartEventmultiplicityCounter"), 0.5); // allMcColl - if (std::fabs(mccollision.posZ()) > vertexZCut) { - return; - } - registry.fill(HIST("hPartEventmultiplicityCounter"), 1.5); // McCollWithVertexZ + JetTableMCPJoined const& jets, aod::JetParticles const& /*particles*/, + soa::SmallGroups const& collisions) + { + bool eventAccepted = false; + double weight = 1.0; + float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent)); + float mcpMult = 0.0f; + float mcpCent = 0.0f; - // outlier check: for every outlier jet, reject the whole event - for (auto const& jet : jets) { - if (jet.pt() > pTHatMaxMCP * pTHat || pTHat < pTHatAbsoluteMin) { - registry.fill(HIST("hPartEventmultiplicityCounter"), 2.5); // RejectedPartCollWithOutliers - return; - } - } + registry.fill(HIST("hPartEventmultiplicityCounter"), 0.5); // allMcColl + if (std::fabs(mccollision.posZ()) > vertexZCut) { + return; + } + registry.fill(HIST("hPartEventmultiplicityCounter"), 1.5); // McCollWithVertexZ - if (doMBGapTrigger && mccollision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { - // Fill rejected MB events; - registry.fill(HIST("hPartEventmultiplicityCounter"), 3.5); // MBRejectedPartEvents + // outlier check: for every outlier jet, reject the whole event + for (auto const& jet : jets) { + if (jet.pt() > pTHatMaxMCP * pTHat || pTHat < pTHatAbsoluteMin) { + registry.fill(HIST("hPartEventmultiplicityCounter"), 2.5); // RejectedPartCollWithOutliers return; } + } - //Perform MC Collision matching, i.e. match the current MC collision to its associated reco (MCD) collision - // to get the corresponding FT0M component at the particle level - auto collisionspermcpjet = collisions.sliceBy(CollisionsPerMCPCollision, mccollision.globalIndex()); - registry.fill(HIST("hRecoMatchesPerMcCollision"), collisionspermcpjet.size()); // for split vertices QA + if (doMBGapTrigger && mccollision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { + // Fill rejected MB events; + registry.fill(HIST("hPartEventmultiplicityCounter"), 3.5); // MBRejectedPartEvents + return; + } - if (collisionspermcpjet.size() == 0 || collisionspermcpjet.size() < 1) { - registry.fill(HIST("hPartEventmultiplicityCounter"), 4.5); // RejectedPartCollForDetCollWithSize0or<1 - return; - } - registry.fill(HIST("hPartEventmultiplicityCounter"), 5.5); // AcceptedPartCollWithSize>1 + // Perform MC Collision matching, i.e. match the current MC collision to its associated reco (MCD) collision + // to get the corresponding FT0M component at the particle level + auto collisionspermcpjet = collisions.sliceBy(CollisionsPerMCPCollision, mccollision.globalIndex()); + registry.fill(HIST("hRecoMatchesPerMcCollision"), collisionspermcpjet.size()); // for split vertices QA - for (auto const& collision : collisionspermcpjet) { - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { - continue; - } - if (doEMCALEventWorkaround) { - if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content - if (collision.alias_bit(kTVXinEMC)) { - eventAccepted = true; - registry.fill(HIST("hPartEventmultiplicityCounter"), 6.5); // EMCreadoutDetEventsWithkTVXinEMC - } - } - } else { - if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { + if (collisionspermcpjet.size() == 0 || collisionspermcpjet.size() < 1) { + registry.fill(HIST("hPartEventmultiplicityCounter"), 4.5); // RejectedPartCollForDetCollWithSize0or<1 + return; + } + registry.fill(HIST("hPartEventmultiplicityCounter"), 5.5); // AcceptedPartCollWithSize>1 + + for (auto const& collision : collisionspermcpjet) { + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { + continue; + } + if (doEMCALEventWorkaround) { + if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content + if (collision.alias_bit(kTVXinEMC)) { eventAccepted = true; registry.fill(HIST("hPartEventmultiplicityCounter"), 6.5); // EMCreadoutDetEventsWithkTVXinEMC } } - mcpMult += collision.multFT0M(); - mcpCent += collision.centFT0M(); - }//collision loop ends - - if (!eventAccepted) { - registry.fill(HIST("hPartEventmultiplicityCounter"), 7.5); // AllRejectedPartEventsAfterEMCEventSelection - return; + } else { + if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { + eventAccepted = true; + registry.fill(HIST("hPartEventmultiplicityCounter"), 6.5); // EMCreadoutDetEventsWithkTVXinEMC + } } - registry.fill(HIST("hPartEventmultiplicityCounter"), 8.5); // EMCAcceptedPartColl - registry.fill(HIST("hMCCollMatchedFT0Mult"), mcpMult); - registry.fill(HIST("hMCCollMatchedFT0Cent"), mcpCent); + mcpMult += collision.multFT0M(); + mcpCent += collision.centFT0M(); + } // collision loop ends - std::vector::iterator> selectedJets; - int nJetsThisEvent = 0; + if (!eventAccepted) { + registry.fill(HIST("hPartEventmultiplicityCounter"), 7.5); // AllRejectedPartEventsAfterEMCEventSelection + return; + } + registry.fill(HIST("hPartEventmultiplicityCounter"), 8.5); // EMCAcceptedPartColl + registry.fill(HIST("hMCCollMatchedFT0Mult"), mcpMult); + registry.fill(HIST("hMCCollMatchedFT0Cent"), mcpCent); - for (auto const& jet : jets) { - if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) continue; - if (jet.phi() < jetPhiMin || jet.phi() > jetPhiMax) continue; - if (!isAcceptedPartJet(jet)) continue; + std::vector::iterator> selectedJets; + int nJetsThisEvent = 0; - selectedJets.push_back(jet); - nJetsThisEvent++; - } - // 1. Sort selected jets by pT before processing - std::sort(selectedJets.begin(), selectedJets.end(), - [](auto const& a, auto const& b) { return (*a).pt() > (*b).pt(); }); + for (auto const& jet : jets) { + if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) + continue; + if (jet.phi() < jetPhiMin || jet.phi() > jetPhiMax) + continue; + if (!isAcceptedPartJet(jet)) + continue; - if (selectedJets.size() == 0) { // no jets = no leading jet - return; - } + selectedJets.push_back(jet); + nJetsThisEvent++; + } + // 1. Sort selected jets by pT before processing + std::sort(selectedJets.begin(), selectedJets.end(), + [](auto const& a, auto const& b) { return (*a).pt() > (*b).pt(); }); - if (selectedJets.size() > 0) { - //Jet multiplicity per event - registry.fill(HIST("h_all_fulljet_Njets_part"), selectedJets.size(), 1.0); - - //Select Leading Jet for N_ch calculation (for every leading jet that is found). There's always one leading jet per event! - auto const& leadingJet = *selectedJets[0]; - auto const& leadingJetPt = leadingJet.pt(); //jet pT distribution of the leading jet - // std::cout << "Leading Jet pT: " << leadingJetPt << std::endl; - registry.fill(HIST("h_Leading_full_jet_pt_part"), leadingJetPt, 1.0); - registry.fill(HIST("h2_full_jet_leadingJetPt_vs_counts_part"), leadingJetPt, nJetsThisEvent, 1.0); - } - - if (selectedJets.size() > 1) { - auto const& subLeadingJet = *selectedJets[1]; - auto const& subLeadingJetPt = subLeadingJet.pt(); //jet pT distribution of the subleading jet i.e. 2nd leading jet - registry.fill(HIST("h_SubLeading_full_jet_pt_part"), subLeadingJetPt, 1.0); - registry.fill(HIST("h2_full_jet_subLeadingJetPt_vs_counts_part"), subLeadingJetPt, nJetsThisEvent, 1.0); - } - // Process ALL selected jets (not just leading) - for (size_t i = 0; i < selectedJets.size(); i++) { - auto const& jet = *selectedJets[i]; - float jetPt = jet.pt(); - bool isLeading = (i == 0); - bool isSubLeading = (i == 1 && selectedJets.size() > 1); //first sub-leading jet - - // Count charged particles(NCh) and neutral particles (NNe) for this jet - int numberOfChargedParticles = 0; int numberOfNeutralParticles = 0; - float neutralEnergy = 0.0f; - for (const auto& constituent : jet.tracks_as()) { - auto pdgParticle = pdgDatabase->GetParticle(constituent.pdgCode()); - if (pdgParticle->Charge() == 0) { - numberOfNeutralParticles++; - neutralEnergy += constituent.e(); - } else { - numberOfChargedParticles++; - } - } - float nef = neutralEnergy / jet.energy(); - - // CASE 1: Fill histograms for ALL selected jets - registry.fill(HIST("h_all_fulljet_pt_part"), jetPt, 1.0); - registry.fill(HIST("h_all_fulljet_Nch_part"), numberOfChargedParticles, 1.0); - registry.fill(HIST("h_all_fulljet_Nne_part"), numberOfNeutralParticles, 1.0); - registry.fill(HIST("h_all_fulljet_NEF_part"), nef, 1.0); - registry.fill(HIST("h2_all_fulljet_jetpT_vs_FT0Mults_part"), jetPt, mcpMult, 1.0); - registry.fill(HIST("h2_all_fulljet_jetpT_vs_Nch_part"), jetPt, numberOfChargedParticles, 1.0); - registry.fill(HIST("h3_full_jet_jetpT_FT0Mults_nef_part"), jetPt, mcpMult, nef, 1.0); - - // CASE 2: Additional leading jet processing - if (isLeading) { - registry.fill(HIST("h_leading_fulljet_pt_part"), jetPt, 1.0); - registry.fill(HIST("h_leading_fulljet_Nch_part"), numberOfChargedParticles, 1.0); - registry.fill(HIST("h_leading_fulljet_Nne_part"), numberOfNeutralParticles, 1.0); - registry.fill(HIST("h_leading_fulljet_NEF_part"), nef, 1.0); - registry.fill(HIST("h2_leading_fulljet_jetpT_vs_FT0Mults_part"), jetPt, mcpMult, 1.0); - registry.fill(HIST("h2_leading_fulljet_jetpT_vs_Nch_part"), jetPt, numberOfChargedParticles, 1.0); - registry.fill(HIST("h3_leading_fulljet_jetpT_FT0Mults_nef_part"), jetPt, mcpMult, nef, 1.0); - } + if (selectedJets.size() == 0) { // no jets = no leading jet + return; + } + + if (selectedJets.size() > 0) { + // Jet multiplicity per event + registry.fill(HIST("h_all_fulljet_Njets_part"), selectedJets.size(), 1.0); + + // Select Leading Jet for N_ch calculation (for every leading jet that is found). There's always one leading jet per event! + auto const& leadingJet = *selectedJets[0]; + auto const& leadingJetPt = leadingJet.pt(); // jet pT distribution of the leading jet + // std::cout << "Leading Jet pT: " << leadingJetPt << std::endl; + registry.fill(HIST("h_Leading_full_jet_pt_part"), leadingJetPt, 1.0); + registry.fill(HIST("h2_full_jet_leadingJetPt_vs_counts_part"), leadingJetPt, nJetsThisEvent, 1.0); + } + + if (selectedJets.size() > 1) { + auto const& subLeadingJet = *selectedJets[1]; + auto const& subLeadingJetPt = subLeadingJet.pt(); // jet pT distribution of the subleading jet i.e. 2nd leading jet + registry.fill(HIST("h_SubLeading_full_jet_pt_part"), subLeadingJetPt, 1.0); + registry.fill(HIST("h2_full_jet_subLeadingJetPt_vs_counts_part"), subLeadingJetPt, nJetsThisEvent, 1.0); + } + // Process ALL selected jets (not just leading) + for (size_t i = 0; i < selectedJets.size(); i++) { + auto const& jet = *selectedJets[i]; + float jetPt = jet.pt(); + bool isLeading = (i == 0); + bool isSubLeading = (i == 1 && selectedJets.size() > 1); // first sub-leading jet - //CASE 3: Additional first sub-leading jet processing - if (isSubLeading) { - registry.fill(HIST("h_subleading_fulljet_pt_part"), jetPt, 1.0); - registry.fill(HIST("h_subleading_fulljet_Nch_part"), numberOfChargedParticles, 1.0); - registry.fill(HIST("h_subleading_fulljet_Nne_part"), numberOfNeutralParticles, 1.0); - registry.fill(HIST("h_subleading_fulljet_NEF_part"), nef, 1.0); - registry.fill(HIST("h2_subleading_fulljet_jetpT_vs_FT0Mults_part"), jetPt, mcpMult, 1.0); - registry.fill(HIST("h2_subleading_fulljet_jetpT_vs_Nch_part"), jetPt, numberOfChargedParticles, 1.0); - registry.fill(HIST("h3_subleading_fulljet_jetpT_FT0Mults_nef_part"), jetPt, mcpMult, nef, 1.0); + // Count charged particles(NCh) and neutral particles (NNe) for this jet + int numberOfChargedParticles = 0; + int numberOfNeutralParticles = 0; + float neutralEnergy = 0.0f; + for (const auto& constituent : jet.tracks_as()) { + auto pdgParticle = pdgDatabase->GetParticle(constituent.pdgCode()); + if (pdgParticle->Charge() == 0) { + numberOfNeutralParticles++; + neutralEnergy += constituent.e(); + } else { + numberOfChargedParticles++; } } + float nef = neutralEnergy / jet.energy(); + + // CASE 1: Fill histograms for ALL selected jets + registry.fill(HIST("h_all_fulljet_pt_part"), jetPt, 1.0); + registry.fill(HIST("h_all_fulljet_Nch_part"), numberOfChargedParticles, 1.0); + registry.fill(HIST("h_all_fulljet_Nne_part"), numberOfNeutralParticles, 1.0); + registry.fill(HIST("h_all_fulljet_NEF_part"), nef, 1.0); + registry.fill(HIST("h2_all_fulljet_jetpT_vs_FT0Mults_part"), jetPt, mcpMult, 1.0); + registry.fill(HIST("h2_all_fulljet_jetpT_vs_Nch_part"), jetPt, numberOfChargedParticles, 1.0); + registry.fill(HIST("h3_full_jet_jetpT_FT0Mults_nef_part"), jetPt, mcpMult, nef, 1.0); + + // CASE 2: Additional leading jet processing + if (isLeading) { + registry.fill(HIST("h_leading_fulljet_pt_part"), jetPt, 1.0); + registry.fill(HIST("h_leading_fulljet_Nch_part"), numberOfChargedParticles, 1.0); + registry.fill(HIST("h_leading_fulljet_Nne_part"), numberOfNeutralParticles, 1.0); + registry.fill(HIST("h_leading_fulljet_NEF_part"), nef, 1.0); + registry.fill(HIST("h2_leading_fulljet_jetpT_vs_FT0Mults_part"), jetPt, mcpMult, 1.0); + registry.fill(HIST("h2_leading_fulljet_jetpT_vs_Nch_part"), jetPt, numberOfChargedParticles, 1.0); + registry.fill(HIST("h3_leading_fulljet_jetpT_FT0Mults_nef_part"), jetPt, mcpMult, nef, 1.0); + } + + // CASE 3: Additional first sub-leading jet processing + if (isSubLeading) { + registry.fill(HIST("h_subleading_fulljet_pt_part"), jetPt, 1.0); + registry.fill(HIST("h_subleading_fulljet_Nch_part"), numberOfChargedParticles, 1.0); + registry.fill(HIST("h_subleading_fulljet_Nne_part"), numberOfNeutralParticles, 1.0); + registry.fill(HIST("h_subleading_fulljet_NEF_part"), nef, 1.0); + registry.fill(HIST("h2_subleading_fulljet_jetpT_vs_FT0Mults_part"), jetPt, mcpMult, 1.0); + registry.fill(HIST("h2_subleading_fulljet_jetpT_vs_Nch_part"), jetPt, numberOfChargedParticles, 1.0); + registry.fill(HIST("h3_subleading_fulljet_jetpT_FT0Mults_nef_part"), jetPt, mcpMult, nef, 1.0); + } } - PROCESS_SWITCH(FullJetSpectra, processMBMCPCollisionsWithMultiplicity, "MB MCP Collisions for Full Jets Multiplicity Studies", false); + } + PROCESS_SWITCH(FullJetSpectra, processMBMCPCollisionsWithMultiplicity, "MB MCP Collisions for Full Jets Multiplicity Studies", false); - void processMBMCPCollisionsWeightedWithMultiplicity(aod::JetMcCollision const& mccollision, - JetTableMCPWeightedJoined const& jets, aod::JetParticles const& /*particles*/, - soa::SmallGroups const& collisions) - { - bool eventAccepted = false; - float pTHat = 10. / (std::pow(mccollision.weight(), 1.0 / pTHatExponent)); - float mcpMult = 0.0f; - float mcpCent= 0.0f; + void processMBMCPCollisionsWeightedWithMultiplicity(aod::JetMcCollision const& mccollision, + JetTableMCPWeightedJoined const& jets, aod::JetParticles const& /*particles*/, + soa::SmallGroups const& collisions) + { + bool eventAccepted = false; + float pTHat = 10. / (std::pow(mccollision.weight(), 1.0 / pTHatExponent)); + float mcpMult = 0.0f; + float mcpCent = 0.0f; - registry.fill(HIST("hPartEventmultiplicityCounter"), 0.5, mccollision.weight()); // allWeightedMcColl - if (std::fabs(mccollision.posZ()) > vertexZCut) { - return; - } - registry.fill(HIST("hPartEventmultiplicityCounter"), 1.5, mccollision.weight()); // WeightedMcCollWithVertexZ + registry.fill(HIST("hPartEventmultiplicityCounter"), 0.5, mccollision.weight()); // allWeightedMcColl + if (std::fabs(mccollision.posZ()) > vertexZCut) { + return; + } + registry.fill(HIST("hPartEventmultiplicityCounter"), 1.5, mccollision.weight()); // WeightedMcCollWithVertexZ - // outlier check: for every outlier jet, reject the whole event - for (auto const& jet : jets) { - if (jet.pt() > pTHatMaxMCP * pTHat || pTHat < pTHatAbsoluteMin) { - registry.fill(HIST("hPartEventmultiplicityCounter"), 2.5, mccollision.weight()); // RejectedWeightedPartCollWithOutliers - return; - } - } + // outlier check: for every outlier jet, reject the whole event + for (auto const& jet : jets) { + if (jet.pt() > pTHatMaxMCP * pTHat || pTHat < pTHatAbsoluteMin) { + registry.fill(HIST("hPartEventmultiplicityCounter"), 2.5, mccollision.weight()); // RejectedWeightedPartCollWithOutliers + return; + } + } - if (doMBGapTrigger && mccollision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { - // Fill rejected MB events; - registry.fill(HIST("hPartEventmultiplicityCounter"), 3.5, mccollision.weight()); // MBRejectedPartEvents - return; - } + if (doMBGapTrigger && mccollision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { + // Fill rejected MB events; + registry.fill(HIST("hPartEventmultiplicityCounter"), 3.5, mccollision.weight()); // MBRejectedPartEvents + return; + } - //Perform MC Collision matching, i.e. match the current MC collision to its associated reco (MCD) collision - // to get the corresponding FT0M component at the particle level - auto collisionspermcpjet = collisions.sliceBy(CollisionsPerMCPCollision, mccollision.globalIndex()); - registry.fill(HIST("hRecoMatchesPerMcCollision"), collisionspermcpjet.size(), mccollision.weight()); // for split vertices QA + // Perform MC Collision matching, i.e. match the current MC collision to its associated reco (MCD) collision + // to get the corresponding FT0M component at the particle level + auto collisionspermcpjet = collisions.sliceBy(CollisionsPerMCPCollision, mccollision.globalIndex()); + registry.fill(HIST("hRecoMatchesPerMcCollision"), collisionspermcpjet.size(), mccollision.weight()); // for split vertices QA - if (collisionspermcpjet.size() == 0 || collisionspermcpjet.size() < 1) { - registry.fill(HIST("hPartEventmultiplicityCounter"), 4.5, mccollision.weight()); // RejectedWeightedPartCollForDetCollWithSize0or<1 - return; - } - registry.fill(HIST("hPartEventmultiplicityCounter"), 5.5, mccollision.weight()); // AcceptedWeightedPartCollWithSize>1 + if (collisionspermcpjet.size() == 0 || collisionspermcpjet.size() < 1) { + registry.fill(HIST("hPartEventmultiplicityCounter"), 4.5, mccollision.weight()); // RejectedWeightedPartCollForDetCollWithSize0or<1 + return; + } + registry.fill(HIST("hPartEventmultiplicityCounter"), 5.5, mccollision.weight()); // AcceptedWeightedPartCollWithSize>1 - for (auto const& collision : collisionspermcpjet) { - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { - continue; - } - if (doEMCALEventWorkaround) { - if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content - if (collision.alias_bit(kTVXinEMC)) { - eventAccepted = true; - registry.fill(HIST("hPartEventmultiplicityCounter"), 6.5, mccollision.weight()); // EMCreadoutDetEventsWithkTVXinEMC - } - } - } else { - if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { - eventAccepted = true; - registry.fill(HIST("hPartEventmultiplicityCounter"), 6.5, mccollision.weight()); // EMCreadoutDetEventsWithkTVXinEMC - } + for (auto const& collision : collisionspermcpjet) { + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, doMBGapTrigger)) { + continue; + } + if (doEMCALEventWorkaround) { + if (collision.isEmcalReadout() && !collision.isAmbiguous()) { // i.e. EMCAL has a cell content + if (collision.alias_bit(kTVXinEMC)) { + eventAccepted = true; + registry.fill(HIST("hPartEventmultiplicityCounter"), 6.5, mccollision.weight()); // EMCreadoutDetEventsWithkTVXinEMC } - mcpMult += collision.multFT0M(); - mcpCent += collision.centFT0M(); - }//collision loop ends - - if (!eventAccepted) { - registry.fill(HIST("hPartEventmultiplicityCounter"), 7.5, mccollision.weight()); // AllRejectedWeightedPartEventsAfterEMCEventSelection - return; } - registry.fill(HIST("hPartEventmultiplicityCounter"), 8.5, mccollision.weight()); // EMCAcceptedWeightedPartColl - registry.fill(HIST("hMCCollMatchedFT0Mult"), mcpMult, mccollision.weight()); - registry.fill(HIST("hMCCollMatchedFT0Cent"), mcpCent, mccollision.weight()); + } else { + if (!collision.isAmbiguous() && jetderiveddatautilities::eventEMCAL(collision) && collision.alias_bit(kTVXinEMC)) { + eventAccepted = true; + registry.fill(HIST("hPartEventmultiplicityCounter"), 6.5, mccollision.weight()); // EMCreadoutDetEventsWithkTVXinEMC + } + } + mcpMult += collision.multFT0M(); + mcpCent += collision.centFT0M(); + } // collision loop ends - std::vector::iterator> selectedJets; - int nJetsThisEvent = 0; + if (!eventAccepted) { + registry.fill(HIST("hPartEventmultiplicityCounter"), 7.5, mccollision.weight()); // AllRejectedWeightedPartEventsAfterEMCEventSelection + return; + } + registry.fill(HIST("hPartEventmultiplicityCounter"), 8.5, mccollision.weight()); // EMCAcceptedWeightedPartColl + registry.fill(HIST("hMCCollMatchedFT0Mult"), mcpMult, mccollision.weight()); + registry.fill(HIST("hMCCollMatchedFT0Cent"), mcpCent, mccollision.weight()); - for (auto const& jet : jets) { - if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) continue; - if (jet.phi() < jetPhiMin || jet.phi() > jetPhiMax) continue; - if (!isAcceptedPartJet(jet)) continue; + std::vector::iterator> selectedJets; + int nJetsThisEvent = 0; - selectedJets.push_back(jet); - nJetsThisEvent++; - } - // 1. Sort selected jets by pT before processing - std::sort(selectedJets.begin(), selectedJets.end(), - [](auto const& a, auto const& b) { return (*a).pt() > (*b).pt(); }); + for (auto const& jet : jets) { + if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) + continue; + if (jet.phi() < jetPhiMin || jet.phi() > jetPhiMax) + continue; + if (!isAcceptedPartJet(jet)) + continue; - if (selectedJets.size() == 0) { // no jets = no leading jet - return; - } + selectedJets.push_back(jet); + nJetsThisEvent++; + } + // 1. Sort selected jets by pT before processing + std::sort(selectedJets.begin(), selectedJets.end(), + [](auto const& a, auto const& b) { return (*a).pt() > (*b).pt(); }); - if (selectedJets.size() > 0) { - //Jet multiplicity per event - registry.fill(HIST("h_all_fulljet_Njets_part"), selectedJets.size(), mccollision.weight()); + if (selectedJets.size() == 0) { // no jets = no leading jet + return; + } - //Select Leading Jet for N_ch calculation (for every leading jet that is found). There's always one leading jet per event! - auto const& leadingJet = *selectedJets[0]; - auto const& leadingJetPt = leadingJet.pt(); //jet pT distribution of the leading jet - // std::cout << "Leading Jet pT: " << leadingJetPt << std::endl; - registry.fill(HIST("h_Leading_full_jet_pt_part"), leadingJetPt, mccollision.weight()); - registry.fill(HIST("h2_full_jet_leadingJetPt_vs_counts_part"), leadingJetPt, nJetsThisEvent, mccollision.weight()); - } + if (selectedJets.size() > 0) { + // Jet multiplicity per event + registry.fill(HIST("h_all_fulljet_Njets_part"), selectedJets.size(), mccollision.weight()); - if (selectedJets.size() > 1) { - auto const& subLeadingJet = *selectedJets[1]; - auto const& subLeadingJetPt = subLeadingJet.pt(); //jet pT distribution of the subleading jet i.e. 2nd leading jet - registry.fill(HIST("h_SubLeading_full_jet_pt_part"), subLeadingJetPt, mccollision.weight()); - registry.fill(HIST("h2_full_jet_subLeadingJetPt_vs_counts_part"), subLeadingJetPt, nJetsThisEvent, mccollision.weight()); - } - // Process ALL selected jets (not just leading) - for (size_t i = 0; i < selectedJets.size(); i++) { - auto const& jet = *selectedJets[i]; - float jetPt = jet.pt(); - bool isLeading = (i == 0); - bool isSubLeading = (i == 1 && selectedJets.size() > 1); //first sub-leading jet - - // Count charged particles(NCh) and neutral particles (NNe) for this jet - int numberOfChargedParticles = 0; int numberOfNeutralParticles = 0; - float neutralEnergy = 0.0f; - for (const auto& constituent : jet.tracks_as()) { - auto pdgParticle = pdgDatabase->GetParticle(constituent.pdgCode()); - if (pdgParticle->Charge() == 0) { - numberOfNeutralParticles++; - neutralEnergy += constituent.e(); - } else { - numberOfChargedParticles++; - } - } - float nef = neutralEnergy / jet.energy(); - - // CASE 1: Fill histograms for ALL selected jets - registry.fill(HIST("h_all_fulljet_pt_part"), jetPt, mccollision.weight()); - registry.fill(HIST("h_all_fulljet_Nch_part"), numberOfChargedParticles, mccollision.weight()); - registry.fill(HIST("h_all_fulljet_Nne_part"), numberOfNeutralParticles, mccollision.weight()); - registry.fill(HIST("h_all_fulljet_NEF_part"), nef, mccollision.weight()); - registry.fill(HIST("h2_all_fulljet_jetpT_vs_FT0Mults_part"), jetPt, mcpMult, mccollision.weight()); - registry.fill(HIST("h2_all_fulljet_jetpT_vs_Nch_part"), jetPt, numberOfChargedParticles, mccollision.weight()); - registry.fill(HIST("h3_full_jet_jetpT_FT0Mults_nef_part"), jetPt, mcpMult, nef, mccollision.weight()); - - // CASE 2: Additional leading jet processing - if (isLeading) { - registry.fill(HIST("h_leading_fulljet_pt_part"), jetPt, mccollision.weight()); - registry.fill(HIST("h_leading_fulljet_Nch_part"), numberOfChargedParticles, mccollision.weight()); - registry.fill(HIST("h_leading_fulljet_Nne_part"), numberOfNeutralParticles, mccollision.weight()); - registry.fill(HIST("h_leading_fulljet_NEF_part"), nef, mccollision.weight()); - registry.fill(HIST("h2_leading_fulljet_jetpT_vs_FT0Mults_part"), jetPt, mcpMult, mccollision.weight()); - registry.fill(HIST("h2_leading_fulljet_jetpT_vs_Nch_part"), jetPt, numberOfChargedParticles, mccollision.weight()); - registry.fill(HIST("h3_leading_fulljet_jetpT_FT0Mults_nef_part"), jetPt, mcpMult, nef, mccollision.weight()); - } + // Select Leading Jet for N_ch calculation (for every leading jet that is found). There's always one leading jet per event! + auto const& leadingJet = *selectedJets[0]; + auto const& leadingJetPt = leadingJet.pt(); // jet pT distribution of the leading jet + // std::cout << "Leading Jet pT: " << leadingJetPt << std::endl; + registry.fill(HIST("h_Leading_full_jet_pt_part"), leadingJetPt, mccollision.weight()); + registry.fill(HIST("h2_full_jet_leadingJetPt_vs_counts_part"), leadingJetPt, nJetsThisEvent, mccollision.weight()); + } - //CASE 3: Additional first sub-leading jet processing - if (isSubLeading) { - registry.fill(HIST("h_subleading_fulljet_pt_part"), jetPt, mccollision.weight()); - registry.fill(HIST("h_subleading_fulljet_Nch_part"), numberOfChargedParticles, mccollision.weight()); - registry.fill(HIST("h_subleading_fulljet_Nne_part"), numberOfNeutralParticles, mccollision.weight()); - registry.fill(HIST("h_subleading_fulljet_NEF_part"), nef, mccollision.weight()); - registry.fill(HIST("h2_subleading_fulljet_jetpT_vs_FT0Mults_part"), jetPt, mcpMult, mccollision.weight()); - registry.fill(HIST("h2_subleading_fulljet_jetpT_vs_Nch_part"), jetPt, numberOfChargedParticles, mccollision.weight()); - registry.fill(HIST("h3_subleading_fulljet_jetpT_FT0Mults_nef_part"), jetPt, mcpMult, nef, mccollision.weight()); - } + if (selectedJets.size() > 1) { + auto const& subLeadingJet = *selectedJets[1]; + auto const& subLeadingJetPt = subLeadingJet.pt(); // jet pT distribution of the subleading jet i.e. 2nd leading jet + registry.fill(HIST("h_SubLeading_full_jet_pt_part"), subLeadingJetPt, mccollision.weight()); + registry.fill(HIST("h2_full_jet_subLeadingJetPt_vs_counts_part"), subLeadingJetPt, nJetsThisEvent, mccollision.weight()); + } + // Process ALL selected jets (not just leading) + for (size_t i = 0; i < selectedJets.size(); i++) { + auto const& jet = *selectedJets[i]; + float jetPt = jet.pt(); + bool isLeading = (i == 0); + bool isSubLeading = (i == 1 && selectedJets.size() > 1); // first sub-leading jet + + // Count charged particles(NCh) and neutral particles (NNe) for this jet + int numberOfChargedParticles = 0; + int numberOfNeutralParticles = 0; + float neutralEnergy = 0.0f; + for (const auto& constituent : jet.tracks_as()) { + auto pdgParticle = pdgDatabase->GetParticle(constituent.pdgCode()); + if (pdgParticle->Charge() == 0) { + numberOfNeutralParticles++; + neutralEnergy += constituent.e(); + } else { + numberOfChargedParticles++; } } - PROCESS_SWITCH(FullJetSpectra, processMBMCPCollisionsWeightedWithMultiplicity, "MB MCP Weighted Collisions for Full Jets Multiplicity Studies", false); + float nef = neutralEnergy / jet.energy(); - }; // struct + // CASE 1: Fill histograms for ALL selected jets + registry.fill(HIST("h_all_fulljet_pt_part"), jetPt, mccollision.weight()); + registry.fill(HIST("h_all_fulljet_Nch_part"), numberOfChargedParticles, mccollision.weight()); + registry.fill(HIST("h_all_fulljet_Nne_part"), numberOfNeutralParticles, mccollision.weight()); + registry.fill(HIST("h_all_fulljet_NEF_part"), nef, mccollision.weight()); + registry.fill(HIST("h2_all_fulljet_jetpT_vs_FT0Mults_part"), jetPt, mcpMult, mccollision.weight()); + registry.fill(HIST("h2_all_fulljet_jetpT_vs_Nch_part"), jetPt, numberOfChargedParticles, mccollision.weight()); + registry.fill(HIST("h3_full_jet_jetpT_FT0Mults_nef_part"), jetPt, mcpMult, nef, mccollision.weight()); + + // CASE 2: Additional leading jet processing + if (isLeading) { + registry.fill(HIST("h_leading_fulljet_pt_part"), jetPt, mccollision.weight()); + registry.fill(HIST("h_leading_fulljet_Nch_part"), numberOfChargedParticles, mccollision.weight()); + registry.fill(HIST("h_leading_fulljet_Nne_part"), numberOfNeutralParticles, mccollision.weight()); + registry.fill(HIST("h_leading_fulljet_NEF_part"), nef, mccollision.weight()); + registry.fill(HIST("h2_leading_fulljet_jetpT_vs_FT0Mults_part"), jetPt, mcpMult, mccollision.weight()); + registry.fill(HIST("h2_leading_fulljet_jetpT_vs_Nch_part"), jetPt, numberOfChargedParticles, mccollision.weight()); + registry.fill(HIST("h3_leading_fulljet_jetpT_FT0Mults_nef_part"), jetPt, mcpMult, nef, mccollision.weight()); + } - WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) - { - return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; + // CASE 3: Additional first sub-leading jet processing + if (isSubLeading) { + registry.fill(HIST("h_subleading_fulljet_pt_part"), jetPt, mccollision.weight()); + registry.fill(HIST("h_subleading_fulljet_Nch_part"), numberOfChargedParticles, mccollision.weight()); + registry.fill(HIST("h_subleading_fulljet_Nne_part"), numberOfNeutralParticles, mccollision.weight()); + registry.fill(HIST("h_subleading_fulljet_NEF_part"), nef, mccollision.weight()); + registry.fill(HIST("h2_subleading_fulljet_jetpT_vs_FT0Mults_part"), jetPt, mcpMult, mccollision.weight()); + registry.fill(HIST("h2_subleading_fulljet_jetpT_vs_Nch_part"), jetPt, numberOfChargedParticles, mccollision.weight()); + registry.fill(HIST("h3_subleading_fulljet_jetpT_FT0Mults_nef_part"), jetPt, mcpMult, nef, mccollision.weight()); } + } + } + PROCESS_SWITCH(FullJetSpectra, processMBMCPCollisionsWeightedWithMultiplicity, "MB MCP Weighted Collisions for Full Jets Multiplicity Studies", false); + +}; // struct + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +}