diff --git a/PWGJE/Tasks/recoilJets.cxx b/PWGJE/Tasks/recoilJets.cxx index d81ec304d59..871aaef7d64 100644 --- a/PWGJE/Tasks/recoilJets.cxx +++ b/PWGJE/Tasks/recoilJets.cxx @@ -86,13 +86,13 @@ using TrackMCLbsTbl = soa::Filtered; using PartTbl = soa::Filtered; // --- Jets (with constituents) -using JetsDataTbl = soa::Filtered>; -using JetsDetTbl = soa::Filtered>; -using JetsPartTbl = soa::Filtered>; +using ChJetsDataTbl = soa::Filtered>; +using ChJetsDetTbl = soa::Filtered>; +using ChJetsPartTbl = soa::Filtered>; // --- Matched jets (det <-> part) -using MatchedJetsDetToPartTbl = soa::Filtered>; -using MatchedJetsPartToDetTbl = soa::Filtered>; +using MatchedChJetsDetToPartTbl = soa::Filtered>; +using MatchedChJetsPartToDetTbl = soa::Filtered>; // --- O2 collisions event selection (not JCollisions) using CollEvSelExtendedIt = soa::Filtered>::iterator; @@ -105,11 +105,11 @@ struct RecoilJets { // ---------- Event selection ---------- struct EvCfg : ConfigurableGroup { std::string prefix = "event"; - Configurable sel{"sel", "sel8", "Choose event selection"}, + Configurable sel{"sel", "sel8+IsGoodZvtxFT0vsPV+NoSameBunchPileup", "Choose event selection"}, triggerMasks{"triggerMasks", "", "Relevant trigger masks: fTrackLowPt,fTrackHighPt"}; Configurable vertexZCut{"vertexZCut", 10., "Accepted z-vertex range"}; - Configurable skipMBGapEvents{"skipMBGapEvents", false, + Configurable skipMBGapEvents{"skipMBGapEvents", true, "Flag to choose to reject min. bias gap events; jet-level rejection " "applied at the jet finder level, here rejection is applied for " "collision and track process functions"}; @@ -118,10 +118,10 @@ struct RecoilJets { // ---------- RCT / flag-based selections ---------- struct Rct : ConfigurableGroup { std::string prefix = "rct"; - Configurable label{"label", "CBT_hadronPID", "Apply rct flag"}; // CBT + Configurable label{"label", "CBT_hadronPID", "Apply rct flag"}; Configurable enable{"enable", true, "Apply RCT selections"}, - requireZDC{"requireZDC", false, "Require ZDC flag"}, + requireZDC{"requireZDC", true, "Require ZDC flag"}, rejectLimitedAcceptance{"rejectLimitedAcceptance", false, "Reject LimitedAcceptance flag"}; } rct; @@ -136,11 +136,11 @@ struct RecoilJets { } trk; // ---------- Jet reconstruction ---------- - struct Jet : ConfigurableGroup { + struct Jets : ConfigurableGroup { std::string prefix = "jet"; Configurable constituentPtMax{"constituentPtMax", 100., "Remove jets with constituent above this pT cut"}, radius{"radius", 0.4, "Jet cone radius"}; - } jet; + } jets; // ---------- Background tools ---------- struct Bkgd : ConfigurableGroup { @@ -155,15 +155,15 @@ struct RecoilJets { // ---------- Normalization FT0 by means ---------- struct FT0A : ConfigurableGroup { std::string prefix = "ft0a"; - Configurable mean{"mean", -1., "Mean FT0A signal"}, - meanPartLevel{"meanPartLevel", -1., "Mean Nch (part level) within FT0A acceptance"}, + Configurable mean{"mean", 3499., "Mean FT0A signal"}, + meanPartLevel{"meanPartLevel", 40.7, "Mean Nch (part level) within FT0A acceptance"}, meanZeq{"meanZeq", -1., "Mean equalized FT0A signal"}; } ft0a; struct FT0C : ConfigurableGroup { std::string prefix = "ft0c"; - Configurable mean{"mean", -1., "Mean FT0C signal"}, - meanPartLevel{"meanPartLevel", -1., "Mean Nch (part level) within FT0C acceptance"}, + Configurable mean{"mean", 844.6, "Mean FT0C signal"}, + meanPartLevel{"meanPartLevel", 43.11, "Mean Nch (part level) within FT0C acceptance"}, meanZeq{"meanZeq", -1., "Mean equalized FT0C signal"}; } ft0c; @@ -173,8 +173,8 @@ struct RecoilJets { Configurable fracSig{"fracSig", 0.9, "Fraction of events used for signal TT"}; Configurable recoilRegion{"recoilRegion", 0.6, "Width of recoil acceptance"}; - Configurable> refPtRange{"refPtRange", {5., 7}, "Reference TT pT range [min,max] (GeV/c)"}, - sigPtRange{"sigPtRange", {20., 50}, "Signal TT pT range [min,max] (GeV/c)"}; + Configurable> refPtRange{"refPtRange", {5., 7.}, "Reference TT pT range [min,max] (GeV/c)"}, + sigPtRange{"sigPtRange", {10., 18.}, "Signal TT pT range [min,max] (GeV/c)"}; Configurable> phiRestr{"phiRestr", {0., 6.3}, "Phi restriction [min,max] (rad) for TT search"}; } tt; @@ -195,13 +195,13 @@ struct RecoilJets { multNBins{"multNBins", 600, "Number of bins for scaled FT0M multiplicity"}, zdcTimeNBins{"zdcTimeNBins", 240, "Number of bins for ZDC timing histograms"}; - ConfigurableAxis multFT0CThresh{"multFT0CThresh", {VARIABLE_WIDTH, 0.0, 0.133, 0.233, 0.367, 0.567, 0.767, 1.067, 1.4, 1.867, 2.5, 3.9, 5.4, 6.9, 20.}, "Percentiles of scaled FT0C: 100%, 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 1%, 0.1%, 0.01%"}; // default values for raw data - ConfigurableAxis multFT0CThreshPartLevel{"multFT0CThreshPartLevel", {VARIABLE_WIDTH, 0.0, 0.133, 0.233, 0.367, 0.567, 0.767, 1.067, 1.4, 1.867, 2.5, 3.9, 5.4, 6.9, 20.}, "Percentiles of scaled FT0C: 100%, 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 1%, 0.1%, 0.01%"}; + ConfigurableAxis multFT0CThresh{"multFT0CThresh", {VARIABLE_WIDTH, 0.0, 0.133333, 0.233333, 0.366667, 0.533333, 0.733333, 1.0, 1.33333, 1.76667, 2.36667, 3.63333, 20.}, "Percentiles of scaled FT0C: 100%, 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 1%, 0%"}; + ConfigurableAxis multFT0CThreshPartLevel{"multFT0CThreshPartLevel", {VARIABLE_WIDTH, 0.0, 0.06, 0.14, 0.3, 0.46, 0.7, 1.0, 1.36, 1.82, 2.42, 3.64, 20.}, "Percentiles of scaled FT0C: 100%, 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 1%, 0%"}; - ConfigurableAxis multFT0MThresh{"multFT0MThresh", {VARIABLE_WIDTH, 0.0, 0.167, 0.267, 0.4, 0.567, 0.8, 1.067, 1.4, 1.833, 2.433, 3.667, 5.1, 6.433, 20.}, "Percentiles of scaled FT0M: 100%, 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 1%, 0.1%, 0.01%"}; // default values for raw data - ConfigurableAxis multFT0MThreshPartLevel{"multFT0MThreshPartLevel", {VARIABLE_WIDTH, 0.0, 0.167, 0.267, 0.4, 0.567, 0.8, 1.067, 1.4, 1.833, 2.433, 3.667, 5.1, 6.433, 20.}, "Percentiles of scaled FT0M: 100%, 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 1%, 0.1%, 0.01%"}; + ConfigurableAxis multFT0MThresh{"multFT0MThresh", {VARIABLE_WIDTH, 0.0, 0.133333, 0.266667, 0.366667, 0.533333, 0.766667, 1.0, 1.33333, 1.76667, 2.33333, 3.4, 20.}, "Percentiles of scaled FT0M: 100%, 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 1%, 0%"}; + ConfigurableAxis multFT0MThreshPartLevel{"multFT0MThreshPartLevel", {VARIABLE_WIDTH, 0.0, 0.08, 0.18, 0.32, 0.5, 0.7, 1.0, 1.36, 1.82, 2.4, 3.5, 20.}, "Percentiles of scaled FT0M: 100%, 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 1%, 0%"}; - ConfigurableAxis axisPtTrackEff{"axisPtTrackEff", {VARIABLE_WIDTH, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.2, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 70.0, 100.0}, "#it{p}_{T} (GeV/#it{c})"}; + ConfigurableAxis axisPtTrackEff{"axisPtTrackEff", {VARIABLE_WIDTH, 0.0, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.2, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 70.0, 100.0}, "#it{p}_{T} (GeV/#it{c})"}; ConfigurableAxis axisCentrality{"axisCentrality", {VARIABLE_WIDTH, -5.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0, 75.0, 80.0, 85.0, 90.0, 95.0, 100.0, 105.0}, "Centrality (%)"}; } hist; @@ -219,8 +219,8 @@ struct RecoilJets { Filter partFilter = nabs(aod::jmcparticle::eta) < trk.etaCut.node(); // Declare filter on jets - Filter jetRadiusFilter = aod::jet::r == nround(jet.radius.node() * 100.); - Filter jetEtaFilter = nabs(aod::jet::eta) < trk.etaCut.node() - jet.radius.node(); // 0.5 in our analysis + Filter jetRadiusFilter = aod::jet::r == nround(jets.radius.node() * 100.); + Filter jetEtaFilter = nabs(aod::jet::eta) < trk.etaCut.node() - jets.radius.node(); // 0.5 in our analysis HistogramRegistry spectra; @@ -228,8 +228,8 @@ struct RecoilJets { int trackSelection = -1; std::vector triggerMaskBits; - Service pdg; - Preslice partJetsByMcCollision = aod::jet::mcCollisionId; + Service pdg{}; + Preslice partJetsByMcCollision = aod::jet::mcCollisionId; Preslice particlesByMcCollision = aod::jet::mcCollisionId; template @@ -237,9 +237,9 @@ struct RecoilJets { AxisDesc(const char* label_, const AxisObject& axis_, const std::string& axisName_ = "") : label(label_), axis(axis_), axisName(axisName_) {} - const char* label; // "FT0C" / "FT0M" - const AxisObject& axis; // AxisSpec or ConfigurableAxis - std::string axisName; // Empty for AxisSpec + const char* label; // "FT0C" / "FT0M" + AxisObject axis; // AxisSpec or ConfigurableAxis + std::string axisName; // Empty for AxisSpec }; struct FT0Metrics { @@ -251,11 +251,55 @@ struct RecoilJets { float scaledFT0M = 0.f; }; + // Indices of the FT0 percentile boundaries in the configurable axis + enum FT0PercentileEdge { + kEA100 = 0, + kEA90, + kEA80, + kEA70, + kEA60, + kEA50, + kEA40, + kEA30, + kEA20, + kEA10, + kEA1, + kEA0 + }; + + // EA intervals used for the rho-shift correction + enum EAInterval { + kMB, + kEA0_20, + kEA0_10, + kEA20_40, + kEA60_80, + kEA50_100, + kEA80_100 + }; + + // Configuration of an EA interval used for the TTRef rho shift + struct EARhoShift { + EAInterval interval; + const char* label; + float rhoShift; + }; + + // Hard-coded rho-shift values obtained from train 701159 + std::array eaRhoShifts{{{.interval = kMB, .label = "EA_MB", .rhoShift = 0.283998f}, + {.interval = kEA0_20, .label = "EA_Perc_0_20", .rhoShift = 0.199249f}, + {.interval = kEA0_10, .label = "EA_Perc_0_10", .rhoShift = 0.186661f}, + {.interval = kEA20_40, .label = "EA_Perc_20_40", .rhoShift = 0.137945f}, + {.interval = kEA60_80, .label = "EA_Perc_60_80", .rhoShift = 0.0962535f}, + {.interval = kEA50_100, .label = "EA_Perc_50_100", .rhoShift = 0.106919f}, + {.interval = kEA80_100, .label = "EA_Perc_80_100", .rhoShift = 0.0871301f}}}; + void init(InitContext const&) { // Initialize histogram axes: configurable AxisSpec pT{hist.jetPtMax, 0.0, hist.jetPtMax * 1., "#it{p}_{T} (GeV/#it{c})"}; AxisSpec jetPTcorr{hist.jetPtMax + 20, -20., hist.jetPtMax * 1.0, "#it{p}_{T, jet}^{ch, corr} (GeV/#it{c})"}; + AxisSpec jetPTcorrFinnerBin{4800, -20., 100., "#it{p}_{T, jet}^{ch, corr} (GeV/#it{c})"}; AxisSpec scaledFT0A{hist.multNBins, 0.0, 20., "FT0A / #LT FT0A #GT"}; AxisSpec scaledFT0C{hist.multNBins, 0.0, 20., "FT0C / #LT FT0C #GT"}; AxisSpec scaledFT0M{hist.multNBins, 0.0, 20., "FT0M^{*}"}; @@ -279,7 +323,7 @@ struct RecoilJets { AxisSpec pseudorap{40, -1., 1., "#it{#eta}"}; AxisSpec pseudorapJets{20, -0.5, 0.5, "#it{#eta}_{jet}"}; AxisSpec jetArea{50, 0.0, 5., "Area_{jet}"}; - AxisSpec rho{50, 0.0, 50., "#it{#rho}"}; + AxisSpec rho{100, 0.0, 50., "#it{#rho}"}; std::string nameFT0Caxis = "FT0C / #LT FT0C #GT"; std::string nameFT0Maxis = "FT0M^{*}"; @@ -299,10 +343,8 @@ struct RecoilJets { std::array, 2> arrConfigurableAxisEAPartLevel = {{{"FT0C", hist.multFT0CThreshPartLevel, nameFT0Caxis}, {"FT0M", hist.multFT0MThreshPartLevel, nameFT0Maxis}}}; - std::array, 4> arrConfigurableAxisCentrality = {{{"CentFT0A", hist.axisCentrality, nameCentralityAxis}, - {"CentFT0C", hist.axisCentrality, nameCentralityAxis}, - {"CentFT0M", hist.axisCentrality, nameCentralityAxis}, - {"CentFT0CVar1", hist.axisCentrality, nameCentralityAxis}}}; + std::array, 2> arrConfigurableAxisCentrality = {{{"CentFT0C", hist.axisCentrality, nameCentralityAxis}, + {"CentFT0M", hist.axisCentrality, nameCentralityAxis}}}; // Zero-degree calorimeter std::array, 3> arrAxisSpecZDCNeutron = {{{"ZNA", zdcNeutronA}, @@ -384,11 +426,11 @@ struct RecoilJets { spectra.add(Form("hScaled%s_Recoil_JetPt_Corr_TTRef", eaAxis.label), Form("Events w. TT_{Ref}: scaled %s & #it{p}_{T} of recoil jets", eaAxis.label), - kTH2F, {{eaAxis.axis, eaAxis.axisName}, jetPTcorr}, hist.sumw2); + kTH2F, {{eaAxis.axis, eaAxis.axisName}, jetPTcorrFinnerBin}, hist.sumw2); spectra.add(Form("hScaled%s_Recoil_JetPt_Corr_TTSig", eaAxis.label), Form("Events w. TT_{Sig}: scaled %s & #it{p}_{T} of recoil jets", eaAxis.label), - kTH2F, {{eaAxis.axis, eaAxis.axisName}, jetPTcorr}, hist.sumw2); + kTH2F, {{eaAxis.axis, eaAxis.axisName}, jetPTcorrFinnerBin}, hist.sumw2); spectra.add(Form("hScaled%s_Recoil_JetPt_TTRef", eaAxis.label), Form("Events w. TT_{Ref}: scaled %s & #it{p}_{T} of recoil jets", eaAxis.label), @@ -417,11 +459,11 @@ struct RecoilJets { spectra.add(Form("hScaled%s_Recoil_JetPt_Corr_TTRef_RestrictedPhi", eaAxis.label), Form("Events w. TT_{Ref}: scaled %s & #it{p}_{T} of recoil jets", eaAxis.label), - kTH2F, {{eaAxis.axis, eaAxis.axisName}, jetPTcorr}, hist.sumw2); + kTH2F, {{eaAxis.axis, eaAxis.axisName}, jetPTcorrFinnerBin}, hist.sumw2); spectra.add(Form("hScaled%s_Recoil_JetPt_Corr_TTSig_RestrictedPhi", eaAxis.label), Form("Events w. TT_{Sig}: scaled %s & #it{p}_{T} of recoil jets", eaAxis.label), - kTH2F, {{eaAxis.axis, eaAxis.axisName}, jetPTcorr}, hist.sumw2); + kTH2F, {{eaAxis.axis, eaAxis.axisName}, jetPTcorrFinnerBin}, hist.sumw2); spectra.add(Form("hScaled%s_DPhi_JetPt_Corr_TTRef_RestrictedPhi", eaAxis.label), Form("Events w. TT_{Ref} #in #varphi (%.2f, %.2f): scaled %s & #Delta#varphi & #it{p}_{T, jet}^{ch}", phiMin, phiMax, eaAxis.label), @@ -502,11 +544,26 @@ struct RecoilJets { spectra.add(Form("h%s_Recoil_JetPt_Corr_TTRef", centAxis.label), Form("Events w. TT_{Ref}: %s & #it{p}_{T} of recoil jets", centAxis.label), - kTH2F, {{centAxis.axis, centAxis.axisName}, jetPTcorr}, hist.sumw2); + kTH2F, {{centAxis.axis, centAxis.axisName}, jetPTcorrFinnerBin}, hist.sumw2); spectra.add(Form("h%s_Recoil_JetPt_Corr_TTSig", centAxis.label), Form("Events w. TT_{Sig}: %s & #it{p}_{T} of recoil jets", centAxis.label), - kTH2F, {{centAxis.axis, centAxis.axisName}, jetPTcorr}, hist.sumw2); + kTH2F, {{centAxis.axis, centAxis.axisName}, jetPTcorrFinnerBin}, hist.sumw2); + + spectra.add(Form("h%s_FT0MStar", centAxis.label), + Form("Correlation of %s vs. FT0M^{*}", centAxis.label), + kTH2F, {{centAxis.axis, centAxis.axisName}, scaledFT0M}, hist.sumw2); + } + + // Register TTRef recoil spectra with rho-shift correction. + for (const auto& ea : eaRhoShifts) { + spectra.add(Form("h%s_Recoil_JetPt_Corr_RhoShifted_TTRef", ea.label), + Form("%s: recoil jet p_{T} (#rho shifted)", ea.label), + kTH1F, {jetPTcorrFinnerBin}, hist.sumw2); + + spectra.add(Form("h%s_RhoShifted_TTRef", ea.label), + Form("%s: #rho shifted in events w. TT_{Ref}", ea.label), + kTH1F, {rho}, hist.sumw2); } } @@ -565,11 +622,11 @@ struct RecoilJets { spectra.add(Form("hScaled%s_Recoil_JetPt_Corr_TTRef_Part", eaAxis.label), Form("MC events w. TT_{Ref}: scaled %s & #it{p}_{T} of recoil jets", eaAxis.label), - kTH2F, {{eaAxis.axis, eaAxis.axisName}, jetPTcorr}, hist.sumw2); + kTH2F, {{eaAxis.axis, eaAxis.axisName}, jetPTcorrFinnerBin}, hist.sumw2); spectra.add(Form("hScaled%s_Recoil_JetPt_Corr_TTSig_Part", eaAxis.label), Form("MC events w. TT_{Sig}: scaled %s & #it{p}_{T} of recoil jets", eaAxis.label), - kTH2F, {{eaAxis.axis, eaAxis.axisName}, jetPTcorr}, hist.sumw2); + kTH2F, {{eaAxis.axis, eaAxis.axisName}, jetPTcorrFinnerBin}, hist.sumw2); spectra.add(Form("hScaled%s_Recoil_JetPt_TTRef_Part", eaAxis.label), Form("MC events w. TT_{Ref}: scaled %s & #it{p}_{T} of recoil jets", eaAxis.label), @@ -598,11 +655,11 @@ struct RecoilJets { spectra.add(Form("hScaled%s_Recoil_JetPt_Corr_TTRef_RestrictedPhi_Part", eaAxis.label), Form("Events w. TT_{Ref}: scaled %s & #it{p}_{T} of recoil jets", eaAxis.label), - kTH2F, {{eaAxis.axis, eaAxis.axisName}, jetPTcorr}, hist.sumw2); + kTH2F, {{eaAxis.axis, eaAxis.axisName}, jetPTcorrFinnerBin}, hist.sumw2); spectra.add(Form("hScaled%s_Recoil_JetPt_Corr_TTSig_RestrictedPhi_Part", eaAxis.label), Form("Events w. TT_{Sig}: scaled %s & #it{p}_{T} of recoil jets", eaAxis.label), - kTH2F, {{eaAxis.axis, eaAxis.axisName}, jetPTcorr}, hist.sumw2); + kTH2F, {{eaAxis.axis, eaAxis.axisName}, jetPTcorrFinnerBin}, hist.sumw2); spectra.add(Form("hScaled%s_DPhi_JetPt_Corr_TTRef_RestrictedPhi_Part", eaAxis.label), Form("MC events w. TT_{Ref} #in #varphi (%.2f, %.2f): scaled %s & #Delta#varphi & #it{p}_{T, jet}^{ch}", phiMin, phiMax, eaAxis.label), @@ -628,15 +685,8 @@ struct RecoilJets { } // Centrality dependence - uint8_t iCent = 0; for (const auto& centAxis : arrConfigurableAxisCentrality) { - int positionCentFT0CVar1 = 3; - if (iCent == positionCentFT0CVar1) { // case of "CentFT0CVar1" - ++iCent; - continue; - } - // Vertex Z position spectra.add(Form("h%s_vertexZ_Part", centAxis.label), "Centrality vs Z vertex of MC collisions", @@ -690,13 +740,11 @@ struct RecoilJets { spectra.add(Form("h%s_Recoil_JetPt_Corr_TTRef_Part", centAxis.label), Form("MC events w. TT_{Ref}: %s & #it{p}_{T} of recoil jets", centAxis.label), - kTH2F, {{centAxis.axis, centAxis.axisName}, jetPTcorr}, hist.sumw2); + kTH2F, {{centAxis.axis, centAxis.axisName}, jetPTcorrFinnerBin}, hist.sumw2); spectra.add(Form("h%s_Recoil_JetPt_Corr_TTSig_Part", centAxis.label), Form("MC events w. TT_{Sig}: %s & #it{p}_{T} of recoil jets", centAxis.label), - kTH2F, {{centAxis.axis, centAxis.axisName}, jetPTcorr}, hist.sumw2); - - ++iCent; + kTH2F, {{centAxis.axis, centAxis.axisName}, jetPTcorrFinnerBin}, hist.sumw2); } } @@ -772,16 +820,8 @@ struct RecoilJets { //==================================================================================== // Centrality-differential corrected spectra and response QA - uint8_t iCent = 0; for (const auto& centAxis : arrConfigurableAxisCentrality) { - int positionCentFT0A = 0; - int positionCentFT0CVar1 = 3; - if (iCent == positionCentFT0A || iCent == positionCentFT0CVar1) { // case of "CentFT0A" & "CentFT0CVar1" - ++iCent; - continue; - } - spectra.add(Form("h%s_ResponseMatrixInclusiveJetsPtCorr", centAxis.label), Form("%s: correlation inclusive #it{p}_{T, det.}^{corr.} vs. #it{p}_{T, part}^{corr.}", centAxis.label), kTH3F, {{centAxis.axis, centAxis.axisName}, detJetPtCorr, partJetPtCorr}, hist.sumw2); @@ -837,8 +877,6 @@ struct RecoilJets { spectra.add(Form("h%s_MissedRecoilJetsPtCorr", centAxis.label), Form("%s: part. level recoil jets w/o matched pair", centAxis.label), kTH2F, {{centAxis.axis, centAxis.axisName}, partJetPtCorr}, hist.sumw2); - - ++iCent; } //==================================================================================== @@ -1254,12 +1292,15 @@ struct RecoilJets { //============================================================================= // Fill histograms with raw or MC det. level data - template + template void fillHistograms(JCollision const& collision, - Jets const& jets, + ChJets const& chJets, JTracks const& tracks, float weight = 1.) { + // Get the configured scaled FT0M percentile boundaries + const std::vector ft0mEdges = hist.multFT0MThresh; + bool bSigEv = false; std::vector vPhiOfTT; double phiTT = 0.; @@ -1274,48 +1315,46 @@ struct RecoilJets { const auto scaledFT0C = ft0Metrics.scaledFT0C; const auto scaledFT0M = ft0Metrics.scaledFT0M; - const auto centFT0A = collision.centFT0A(); const auto centFT0C = collision.centFT0C(); const auto centFT0M = collision.centFT0M(); - const auto centFT0CVar1 = collision.centFT0CVariant1(); const auto dice = randGen->Rndm(); - if (dice < tt.fracSig) + if (dice < tt.fracSig) { bSigEv = true; + } // EA distributions spectra.fill(HIST("hScaledFT0C"), scaledFT0C, weight); spectra.fill(HIST("hScaledFT0M"), scaledFT0M, weight); // Centrality distribution - spectra.fill(HIST("hCentFT0A"), centFT0A, weight); spectra.fill(HIST("hCentFT0C"), centFT0C, weight); spectra.fill(HIST("hCentFT0M"), centFT0M, weight); - spectra.fill(HIST("hCentFT0CVar1"), centFT0CVar1, weight); + + // Correlation: centrality vs FT0M* + spectra.fill(HIST("hCentFT0C_FT0MStar"), centFT0C, scaledFT0M, weight); + spectra.fill(HIST("hCentFT0M_FT0MStar"), centFT0M, scaledFT0M, weight); // Z vertex position vs EA / centrality spectra.fill(HIST("hScaledFT0C_vertexZ"), scaledFT0C, vertexZ, weight); spectra.fill(HIST("hScaledFT0M_vertexZ"), scaledFT0M, vertexZ, weight); - spectra.fill(HIST("hCentFT0A_vertexZ"), centFT0A, vertexZ, weight); spectra.fill(HIST("hCentFT0C_vertexZ"), centFT0C, vertexZ, weight); spectra.fill(HIST("hCentFT0M_vertexZ"), centFT0M, vertexZ, weight); - spectra.fill(HIST("hCentFT0CVar1_vertexZ"), centFT0CVar1, vertexZ, weight); // Rho vs EA / centrality spectra.fill(HIST("hScaledFT0C_Rho"), scaledFT0C, rho, weight); spectra.fill(HIST("hScaledFT0M_Rho"), scaledFT0M, rho, weight); - spectra.fill(HIST("hCentFT0A_Rho"), centFT0A, rho, weight); spectra.fill(HIST("hCentFT0C_Rho"), centFT0C, rho, weight); spectra.fill(HIST("hCentFT0M_Rho"), centFT0M, rho, weight); - spectra.fill(HIST("hCentFT0CVar1_Rho"), centFT0CVar1, rho, weight); for (const auto& track : tracks) { spectra.fill(HIST("hTrackSelectionCount"), 0.5); - if (skipTrack(track)) + if (skipTrack(track)) { continue; + } const float trackPt = track.pt(); const float trackPhi = track.phi(); @@ -1325,10 +1364,8 @@ struct RecoilJets { spectra.fill(HIST("hScaledFT0CTrackPtEtaPhi"), scaledFT0C, trackPt, trackEta, trackPhi, weight); spectra.fill(HIST("hScaledFT0MTrackPtEtaPhi"), scaledFT0M, trackPt, trackEta, trackPhi, weight); - spectra.fill(HIST("hCentFT0ATrackPtEtaPhi"), centFT0A, trackPt, trackEta, trackPhi, weight); spectra.fill(HIST("hCentFT0CTrackPtEtaPhi"), centFT0C, trackPt, trackEta, trackPhi, weight); spectra.fill(HIST("hCentFT0MTrackPtEtaPhi"), centFT0M, trackPt, trackEta, trackPhi, weight); - spectra.fill(HIST("hCentFT0CVar1TrackPtEtaPhi"), centFT0CVar1, trackPt, trackEta, trackPhi, weight); // Search for TT candidate const auto ptTTsigMin = tt.sigPtRange->at(0); @@ -1375,20 +1412,14 @@ struct RecoilJets { //_____________________________________________________ // Centrality - spectra.fill(HIST("hCentFT0A_Ntrig"), centFT0A, addCountToTTSig, weight); spectra.fill(HIST("hCentFT0C_Ntrig"), centFT0C, addCountToTTSig, weight); spectra.fill(HIST("hCentFT0M_Ntrig"), centFT0M, addCountToTTSig, weight); - spectra.fill(HIST("hCentFT0CVar1_Ntrig"), centFT0CVar1, addCountToTTSig, weight); - spectra.fill(HIST("hCentFT0A_TTSig"), centFT0A, weight); spectra.fill(HIST("hCentFT0C_TTSig"), centFT0C, weight); spectra.fill(HIST("hCentFT0M_TTSig"), centFT0M, weight); - spectra.fill(HIST("hCentFT0CVar1_TTSig"), centFT0CVar1, weight); - spectra.fill(HIST("hCentFT0A_Rho_TTSig"), centFT0A, rho, weight); spectra.fill(HIST("hCentFT0C_Rho_TTSig"), centFT0C, rho, weight); spectra.fill(HIST("hCentFT0M_Rho_TTSig"), centFT0M, rho, weight); - spectra.fill(HIST("hCentFT0CVar1_Rho_TTSig"), centFT0CVar1, rho, weight); } else { spectra.fill(HIST("hScaledFT0C_Ntrig"), scaledFT0C, addCountToTTRef, weight); @@ -1410,39 +1441,34 @@ struct RecoilJets { //_____________________________________________________ // Centrality - spectra.fill(HIST("hCentFT0A_Ntrig"), centFT0A, addCountToTTRef, weight); spectra.fill(HIST("hCentFT0C_Ntrig"), centFT0C, addCountToTTRef, weight); spectra.fill(HIST("hCentFT0M_Ntrig"), centFT0M, addCountToTTRef, weight); - spectra.fill(HIST("hCentFT0CVar1_Ntrig"), centFT0CVar1, addCountToTTRef, weight); - spectra.fill(HIST("hCentFT0A_TTRef"), centFT0A, weight); spectra.fill(HIST("hCentFT0C_TTRef"), centFT0C, weight); spectra.fill(HIST("hCentFT0M_TTRef"), centFT0M, weight); - spectra.fill(HIST("hCentFT0CVar1_TTRef"), centFT0CVar1, weight); - spectra.fill(HIST("hCentFT0A_Rho_TTRef"), centFT0A, rho, weight); spectra.fill(HIST("hCentFT0C_Rho_TTRef"), centFT0C, rho, weight); spectra.fill(HIST("hCentFT0M_Rho_TTRef"), centFT0M, rho, weight); - spectra.fill(HIST("hCentFT0CVar1_Rho_TTRef"), centFT0CVar1, rho, weight); } } - for (const auto& jet : jets) { + for (const auto& chJet : chJets) { // skip jets which have a constituent with pT above specified cut - if (isJetWithHighPtConstituent(jet)) + if (isJetWithHighPtConstituent(chJet)) { continue; + } - const auto jetPt = jet.pt(); - const auto jetArea = jet.area(); + const auto jetPt = chJet.pt(); + const auto jetArea = chJet.area(); const auto jetPtCorr = jetPt - rho * jetArea; - spectra.fill(HIST("hJetPtEtaPhiRhoArea"), jetPt, jet.eta(), jet.phi(), rho, jetArea, weight); + spectra.fill(HIST("hJetPtEtaPhiRhoArea"), jetPt, chJet.eta(), chJet.phi(), rho, jetArea, weight); if (nTT > 0) { const auto phiMin = tt.phiRestr->at(0); const auto phiMax = tt.phiRestr->at(1); - auto [dphi, bRecoilJet] = isRecoilJet(jet, phiTT); + auto [dphi, bRecoilJet] = isRecoilJet(chJet, phiTT); if (bSigEv) { @@ -1460,10 +1486,8 @@ struct RecoilJets { } // Centrality dependence - spectra.fill(HIST("hCentFT0A_DPhi_JetPt_Corr_TTSig"), centFT0A, dphi, jetPtCorr, weight); spectra.fill(HIST("hCentFT0C_DPhi_JetPt_Corr_TTSig"), centFT0C, dphi, jetPtCorr, weight); spectra.fill(HIST("hCentFT0M_DPhi_JetPt_Corr_TTSig"), centFT0M, dphi, jetPtCorr, weight); - spectra.fill(HIST("hCentFT0CVar1_DPhi_JetPt_Corr_TTSig"), centFT0CVar1, dphi, jetPtCorr, weight); if (bRecoilJet) { @@ -1474,10 +1498,8 @@ struct RecoilJets { spectra.fill(HIST("hScaledFT0M_Recoil_JetPt_TTSig"), scaledFT0M, jetPt, weight); // Centrality dependence - spectra.fill(HIST("hCentFT0A_Recoil_JetPt_Corr_TTSig"), centFT0A, jetPtCorr, weight); spectra.fill(HIST("hCentFT0C_Recoil_JetPt_Corr_TTSig"), centFT0C, jetPtCorr, weight); spectra.fill(HIST("hCentFT0M_Recoil_JetPt_Corr_TTSig"), centFT0M, jetPtCorr, weight); - spectra.fill(HIST("hCentFT0CVar1_Recoil_JetPt_Corr_TTSig"), centFT0CVar1, jetPtCorr, weight); if (phiTT > phiMin && phiTT < phiMax) { spectra.fill(HIST("hScaledFT0C_Recoil_JetPt_Corr_TTSig_RestrictedPhi"), scaledFT0C, jetPtCorr, weight); @@ -1500,13 +1522,59 @@ struct RecoilJets { } // Centrality dependence - spectra.fill(HIST("hCentFT0A_DPhi_JetPt_Corr_TTRef"), centFT0A, dphi, jetPtCorr, weight); spectra.fill(HIST("hCentFT0C_DPhi_JetPt_Corr_TTRef"), centFT0C, dphi, jetPtCorr, weight); spectra.fill(HIST("hCentFT0M_DPhi_JetPt_Corr_TTRef"), centFT0M, dphi, jetPtCorr, weight); - spectra.fill(HIST("hCentFT0CVar1_DPhi_JetPt_Corr_TTRef"), centFT0CVar1, dphi, jetPtCorr, weight); if (bRecoilJet) { + // Fill EA-dependent TTRef recoil spectra using the corresponding rho shift + for (const auto& ea : eaRhoShifts) { + + if (!isInEAInterval(scaledFT0M, ft0mEdges, ea.interval)) { + continue; + } + + const float rhoRefShifted = rho + ea.rhoShift; + const float jetPtCorrShifted = jetPt - rhoRefShifted * jetArea; + + switch (ea.interval) { + case kMB: + spectra.fill(HIST("hEA_MB_Recoil_JetPt_Corr_RhoShifted_TTRef"), jetPtCorrShifted, weight); + spectra.fill(HIST("hEA_MB_RhoShifted_TTRef"), rhoRefShifted, weight); + break; + + case kEA0_20: + spectra.fill(HIST("hEA_Perc_0_20_Recoil_JetPt_Corr_RhoShifted_TTRef"), jetPtCorrShifted, weight); + spectra.fill(HIST("hEA_Perc_0_20_RhoShifted_TTRef"), rhoRefShifted, weight); + break; + + case kEA0_10: + spectra.fill(HIST("hEA_Perc_0_10_Recoil_JetPt_Corr_RhoShifted_TTRef"), jetPtCorrShifted, weight); + spectra.fill(HIST("hEA_Perc_0_10_RhoShifted_TTRef"), rhoRefShifted, weight); + break; + + case kEA20_40: + spectra.fill(HIST("hEA_Perc_20_40_Recoil_JetPt_Corr_RhoShifted_TTRef"), jetPtCorrShifted, weight); + spectra.fill(HIST("hEA_Perc_20_40_RhoShifted_TTRef"), rhoRefShifted, weight); + break; + + case kEA60_80: + spectra.fill(HIST("hEA_Perc_60_80_Recoil_JetPt_Corr_RhoShifted_TTRef"), jetPtCorrShifted, weight); + spectra.fill(HIST("hEA_Perc_60_80_RhoShifted_TTRef"), rhoRefShifted, weight); + break; + + case kEA50_100: + spectra.fill(HIST("hEA_Perc_50_100_Recoil_JetPt_Corr_RhoShifted_TTRef"), jetPtCorrShifted, weight); + spectra.fill(HIST("hEA_Perc_50_100_RhoShifted_TTRef"), rhoRefShifted, weight); + break; + + case kEA80_100: + spectra.fill(HIST("hEA_Perc_80_100_Recoil_JetPt_Corr_RhoShifted_TTRef"), jetPtCorrShifted, weight); + spectra.fill(HIST("hEA_Perc_80_100_RhoShifted_TTRef"), rhoRefShifted, weight); + break; + } + } + // EA dependence spectra.fill(HIST("hScaledFT0C_Recoil_JetPt_Corr_TTRef"), scaledFT0C, jetPtCorr, weight); spectra.fill(HIST("hScaledFT0M_Recoil_JetPt_Corr_TTRef"), scaledFT0M, jetPtCorr, weight); @@ -1514,10 +1582,8 @@ struct RecoilJets { spectra.fill(HIST("hScaledFT0M_Recoil_JetPt_TTRef"), scaledFT0M, jetPt, weight); // Centrality dependence - spectra.fill(HIST("hCentFT0A_Recoil_JetPt_Corr_TTRef"), centFT0A, jetPtCorr, weight); spectra.fill(HIST("hCentFT0C_Recoil_JetPt_Corr_TTRef"), centFT0C, jetPtCorr, weight); spectra.fill(HIST("hCentFT0M_Recoil_JetPt_Corr_TTRef"), centFT0M, jetPtCorr, weight); - spectra.fill(HIST("hCentFT0CVar1_Recoil_JetPt_Corr_TTRef"), centFT0CVar1, jetPtCorr, weight); if (phiTT > phiMin && phiTT < phiMax) { spectra.fill(HIST("hScaledFT0C_Recoil_JetPt_Corr_TTRef_RestrictedPhi"), scaledFT0C, jetPtCorr, weight); @@ -1529,9 +1595,9 @@ struct RecoilJets { } } - template + template void fillHistogramsMCPartLevel(JCollision const& collision, - Jets const& jets, + ChJets const& chJets, JParticles const& particles, float weight = 1.) { @@ -1549,20 +1615,19 @@ struct RecoilJets { const auto scaledFT0C = ft0Metrics.scaledFT0C; const auto scaledFT0M = ft0Metrics.scaledFT0M; - const auto centFT0A = collision.centFT0A(); const auto centFT0C = collision.centFT0C(); const auto centFT0M = collision.centFT0M(); const auto dice = randGen->Rndm(); - if (dice < tt.fracSig) + if (dice < tt.fracSig) { bSigEv = true; + } // EA distributions spectra.fill(HIST("hScaledFT0C_Part"), scaledFT0C, weight); spectra.fill(HIST("hScaledFT0M_Part"), scaledFT0M, weight); // Centrality distribution - spectra.fill(HIST("hCentFT0A_Part"), centFT0A, weight); spectra.fill(HIST("hCentFT0C_Part"), centFT0C, weight); spectra.fill(HIST("hCentFT0M_Part"), centFT0M, weight); @@ -1570,7 +1635,6 @@ struct RecoilJets { spectra.fill(HIST("hScaledFT0C_vertexZ_Part"), scaledFT0C, vertexZ, weight); spectra.fill(HIST("hScaledFT0M_vertexZ_Part"), scaledFT0M, vertexZ, weight); - spectra.fill(HIST("hCentFT0A_vertexZ_Part"), centFT0A, vertexZ, weight); spectra.fill(HIST("hCentFT0C_vertexZ_Part"), centFT0C, vertexZ, weight); spectra.fill(HIST("hCentFT0M_vertexZ_Part"), centFT0M, vertexZ, weight); @@ -1579,8 +1643,9 @@ struct RecoilJets { spectra.fill(HIST("hScaledFT0M_Rho_Part"), scaledFT0M, rho, weight); for (const auto& particle : particles) { - if (skipParticle(particle)) + if (skipParticle(particle)) { continue; + } const auto particlePt = particle.pt(); const auto particlePhi = particle.phi(); @@ -1589,7 +1654,6 @@ struct RecoilJets { spectra.fill(HIST("hScaledFT0CTrackPtEtaPhi_Part"), scaledFT0C, particlePt, particleEta, particlePhi, weight); spectra.fill(HIST("hScaledFT0MTrackPtEtaPhi_Part"), scaledFT0M, particlePt, particleEta, particlePhi, weight); - spectra.fill(HIST("hCentFT0ATrackPtEtaPhi_Part"), centFT0A, particlePt, particleEta, particlePhi, weight); spectra.fill(HIST("hCentFT0CTrackPtEtaPhi_Part"), centFT0C, particlePt, particleEta, particlePhi, weight); spectra.fill(HIST("hCentFT0MTrackPtEtaPhi_Part"), centFT0M, particlePt, particleEta, particlePhi, weight); @@ -1638,15 +1702,12 @@ struct RecoilJets { //_____________________________________________________ // Centrality - spectra.fill(HIST("hCentFT0A_Ntrig_Part"), centFT0A, addCountToTTSig, weight); spectra.fill(HIST("hCentFT0C_Ntrig_Part"), centFT0C, addCountToTTSig, weight); spectra.fill(HIST("hCentFT0M_Ntrig_Part"), centFT0M, addCountToTTSig, weight); - spectra.fill(HIST("hCentFT0A_TTSig_Part"), centFT0A, weight); spectra.fill(HIST("hCentFT0C_TTSig_Part"), centFT0C, weight); spectra.fill(HIST("hCentFT0M_TTSig_Part"), centFT0M, weight); - spectra.fill(HIST("hCentFT0A_Rho_TTSig_Part"), centFT0A, rho, weight); spectra.fill(HIST("hCentFT0C_Rho_TTSig_Part"), centFT0C, rho, weight); spectra.fill(HIST("hCentFT0M_Rho_TTSig_Part"), centFT0M, rho, weight); @@ -1670,32 +1731,29 @@ struct RecoilJets { //_____________________________________________________ // Centrality - spectra.fill(HIST("hCentFT0A_Ntrig_Part"), centFT0A, addCountToTTRef, weight); spectra.fill(HIST("hCentFT0C_Ntrig_Part"), centFT0C, addCountToTTRef, weight); spectra.fill(HIST("hCentFT0M_Ntrig_Part"), centFT0M, addCountToTTRef, weight); - spectra.fill(HIST("hCentFT0A_TTRef_Part"), centFT0A, weight); spectra.fill(HIST("hCentFT0C_TTRef_Part"), centFT0C, weight); spectra.fill(HIST("hCentFT0M_TTRef_Part"), centFT0M, weight); - spectra.fill(HIST("hCentFT0A_Rho_TTRef_Part"), centFT0A, rho, weight); spectra.fill(HIST("hCentFT0C_Rho_TTRef_Part"), centFT0C, rho, weight); spectra.fill(HIST("hCentFT0M_Rho_TTRef_Part"), centFT0M, rho, weight); } } - for (const auto& jet : jets) { - const auto jetPt = jet.pt(); - const auto jetArea = jet.area(); + for (const auto& chJet : chJets) { + const auto jetPt = chJet.pt(); + const auto jetArea = chJet.area(); const auto jetPtCorr = jetPt - rho * jetArea; - spectra.fill(HIST("hJetPtEtaPhiRhoArea_Part"), jetPt, jet.eta(), jet.phi(), rho, jetArea, weight); + spectra.fill(HIST("hJetPtEtaPhiRhoArea_Part"), jetPt, chJet.eta(), chJet.phi(), rho, jetArea, weight); if (nTT > 0) { const auto phiMin = tt.phiRestr->at(0); const auto phiMax = tt.phiRestr->at(1); - auto [dphi, bRecoilJet] = isRecoilJet(jet, phiTT); + auto [dphi, bRecoilJet] = isRecoilJet(chJet, phiTT); if (bSigEv) { @@ -1713,7 +1771,6 @@ struct RecoilJets { } // Centrality dependence - spectra.fill(HIST("hCentFT0A_DPhi_JetPt_Corr_TTSig_Part"), centFT0A, dphi, jetPtCorr, weight); spectra.fill(HIST("hCentFT0C_DPhi_JetPt_Corr_TTSig_Part"), centFT0C, dphi, jetPtCorr, weight); spectra.fill(HIST("hCentFT0M_DPhi_JetPt_Corr_TTSig_Part"), centFT0M, dphi, jetPtCorr, weight); @@ -1726,7 +1783,6 @@ struct RecoilJets { spectra.fill(HIST("hScaledFT0M_Recoil_JetPt_TTSig_Part"), scaledFT0M, jetPt, weight); // Centrality dependence - spectra.fill(HIST("hCentFT0A_Recoil_JetPt_Corr_TTSig_Part"), centFT0A, jetPtCorr, weight); spectra.fill(HIST("hCentFT0C_Recoil_JetPt_Corr_TTSig_Part"), centFT0C, jetPtCorr, weight); spectra.fill(HIST("hCentFT0M_Recoil_JetPt_Corr_TTSig_Part"), centFT0M, jetPtCorr, weight); @@ -1751,7 +1807,6 @@ struct RecoilJets { } // Centrality dependence - spectra.fill(HIST("hCentFT0A_DPhi_JetPt_Corr_TTRef_Part"), centFT0A, dphi, jetPtCorr, weight); spectra.fill(HIST("hCentFT0C_DPhi_JetPt_Corr_TTRef_Part"), centFT0C, dphi, jetPtCorr, weight); spectra.fill(HIST("hCentFT0M_DPhi_JetPt_Corr_TTRef_Part"), centFT0M, dphi, jetPtCorr, weight); @@ -1764,7 +1819,6 @@ struct RecoilJets { spectra.fill(HIST("hScaledFT0M_Recoil_JetPt_TTRef_Part"), scaledFT0M, jetPt, weight); // Centrality dependence - spectra.fill(HIST("hCentFT0A_Recoil_JetPt_Corr_TTRef_Part"), centFT0A, jetPtCorr, weight); spectra.fill(HIST("hCentFT0C_Recoil_JetPt_Corr_TTRef_Part"), centFT0C, jetPtCorr, weight); spectra.fill(HIST("hCentFT0M_Recoil_JetPt_Corr_TTRef_Part"), centFT0M, jetPtCorr, weight); @@ -1781,10 +1835,10 @@ struct RecoilJets { //============================================================================= // Construction of response matrix //============================================================================= - template + template void fillMatchedGeoHistograms(JCollision const& collision, - JetsPart const& jetsPart, - JetsDet const& jetsDet, + ChJetsPart const& chJetsPart, + ChJetsDet const& chJetsDet, JTracks const& tracks, const float partLevelCollRho, float weight = 1.) @@ -1813,10 +1867,10 @@ struct RecoilJets { //==================================================================================== // Particle-level jets - for (const auto& jetPart : jetsPart) { - const auto partJetPt = jetPart.pt(); - const auto partJetPtCorr = partJetPt - partLevelCollRho * jetPart.area(); - const bool bIsPartJetRecoil = bIsThereTTSig && get<1>(isRecoilJet(jetPart, phiTTSig)); + for (const auto& chJetPart : chJetsPart) { + const auto partJetPt = chJetPart.pt(); + const auto partJetPtCorr = partJetPt - partLevelCollRho * chJetPart.area(); + const bool bIsPartJetRecoil = bIsThereTTSig && get<1>(isRecoilJet(chJetPart, phiTTSig)); //================================================================================== // All particle-level inclusive jets @@ -1840,8 +1894,8 @@ struct RecoilJets { //================================================================================== // Matched detector-level jets - if (jetPart.has_matchedJetGeo()) { - const auto& jetsDetMatched = jetPart.template matchedJetGeo_as(); + if (chJetPart.has_matchedJetGeo()) { + const auto& jetsDetMatched = chJetPart.template matchedJetGeo_as(); for (const auto& jetDetMatched : jetsDetMatched) { if (isJetWithHighPtConstituent(jetDetMatched)) { @@ -1852,7 +1906,7 @@ struct RecoilJets { const auto detJetPtCorr = detJetPt - detLevelCollRho * jetDetMatched.area(); const auto relJetSmearPt = (partJetPt - detJetPt) / partJetPt; const auto relJetSmearPtCorr = (partJetPtCorr - detJetPtCorr) / partJetPtCorr; - const auto deltaPhi = RecoDecay::constrainAngle(jetPart.phi() - jetDetMatched.phi(), -constants::math::PI); + const auto deltaPhi = RecoDecay::constrainAngle(chJetPart.phi() - jetDetMatched.phi(), -constants::math::PI); hasAcceptedDetMatch = true; //================================================================================ @@ -1926,15 +1980,15 @@ struct RecoilJets { //==================================================================================== // Detector-level jets - for (const auto& jetDet : jetsDet) { - if (isJetWithHighPtConstituent(jetDet)) { + for (const auto& chJetDet : chJetsDet) { + if (isJetWithHighPtConstituent(chJetDet)) { continue; } - const auto detJetPt = jetDet.pt(); - const auto detJetPtCorr = detJetPt - detLevelCollRho * jetDet.area(); + const auto detJetPt = chJetDet.pt(); + const auto detJetPtCorr = detJetPt - detLevelCollRho * chJetDet.area(); - const bool bIsJetRecoil = bIsThereTTSig && get<1>(isRecoilJet(jetDet, phiTTSig)); + const bool bIsJetRecoil = bIsThereTTSig && get<1>(isRecoilJet(chJetDet, phiTTSig)); //================================================================================== // All detector-level inclusive jets @@ -1956,7 +2010,7 @@ struct RecoilJets { //================================================================================== // Fake detector-level jets: no matched particle level jet - if (!jetDet.has_matchedJetGeo()) { + if (!chJetDet.has_matchedJetGeo()) { spectra.fill(HIST("hFakeInclusiveJetsPt"), detJetPt, weight); spectra.fill(HIST("hScaledFT0C_FakeInclusiveJetsPtCorr"), scaledFT0C, detJetPtCorr, weight); spectra.fill(HIST("hScaledFT0M_FakeInclusiveJetsPtCorr"), scaledFT0M, detJetPtCorr, weight); @@ -1974,10 +2028,10 @@ struct RecoilJets { } } - template + template void fillMatchedGeoPtHistograms(JCollision const& collision, - JetsPart const& jetsPart, - JetsDet const& jetsDet, + ChJetsPart const& chJetsPart, + ChJetsDet const& chJetsDet, JTracks const& tracks, const float partLevelCollRho, float weight = 1.) @@ -2005,10 +2059,10 @@ struct RecoilJets { //==================================================================================== // Particle-level jets - for (const auto& jetPart : jetsPart) { - const auto partJetPt = jetPart.pt(); - const auto partJetPtCorr = partJetPt - partLevelCollRho * jetPart.area(); - const bool bIsPartJetRecoil = bIsThereTTSig && get<1>(isRecoilJet(jetPart, phiTTSig)); + for (const auto& chJetPart : chJetsPart) { + const auto partJetPt = chJetPart.pt(); + const auto partJetPtCorr = partJetPt - partLevelCollRho * chJetPart.area(); + const bool bIsPartJetRecoil = bIsThereTTSig && get<1>(isRecoilJet(chJetPart, phiTTSig)); //================================================================================== // All particle-level inclusive jets @@ -2032,9 +2086,9 @@ struct RecoilJets { //================================================================================== // Matched detector-level jets: geo + pT matching - if (jetPart.has_matchedJetGeo() && jetPart.has_matchedJetPt()) { - const auto& jetsDetMatched = jetPart.template matchedJetGeo_as(); - const auto both = intersectMatchIds(jetPart.matchedJetGeoIds(), jetPart.matchedJetPtIds()); + if (chJetPart.has_matchedJetGeo() && chJetPart.has_matchedJetPt()) { + const auto& jetsDetMatched = chJetPart.template matchedJetGeo_as(); + const auto both = intersectMatchIds(chJetPart.matchedJetGeoIds(), chJetPart.matchedJetPtIds()); for (const auto& jetDetMatched : jetsDetMatched) { if (!both.contains(jetDetMatched.globalIndex())) { @@ -2049,7 +2103,7 @@ struct RecoilJets { const auto detJetPtCorr = detJetPt - detLevelCollRho * jetDetMatched.area(); const auto relJetSmearPt = (partJetPt - detJetPt) / partJetPt; const auto relJetSmearPtCorr = (partJetPtCorr - detJetPtCorr) / partJetPtCorr; - const auto deltaPhi = RecoDecay::constrainAngle(jetPart.phi() - jetDetMatched.phi(), -constants::math::PI); + const auto deltaPhi = RecoDecay::constrainAngle(chJetPart.phi() - jetDetMatched.phi(), -constants::math::PI); hasAcceptedDetMatch = true; //================================================================================ @@ -2123,15 +2177,15 @@ struct RecoilJets { //==================================================================================== // Detector-level jets - for (const auto& jetDet : jetsDet) { - if (isJetWithHighPtConstituent(jetDet)) { + for (const auto& chJetDet : chJetsDet) { + if (isJetWithHighPtConstituent(chJetDet)) { continue; } - const auto detJetPt = jetDet.pt(); - const auto detJetPtCorr = detJetPt - detLevelCollRho * jetDet.area(); + const auto detJetPt = chJetDet.pt(); + const auto detJetPtCorr = detJetPt - detLevelCollRho * chJetDet.area(); - const bool bIsJetRecoil = bIsThereTTSig && get<1>(isRecoilJet(jetDet, phiTTSig)); + const bool bIsJetRecoil = bIsThereTTSig && get<1>(isRecoilJet(chJetDet, phiTTSig)); //================================================================================== // All detector-level inclusive jets @@ -2155,8 +2209,8 @@ struct RecoilJets { // Fake detector-level jets: no matched particle-level jet accepted by both geo and pT matching bool hasAcceptedPartMatch = false; - if (jetDet.has_matchedJetGeo() && jetDet.has_matchedJetPt()) { - const auto both = intersectMatchIds(jetDet.matchedJetGeoIds(), jetDet.matchedJetPtIds()); + if (chJetDet.has_matchedJetGeo() && chJetDet.has_matchedJetPt()) { + const auto both = intersectMatchIds(chJetDet.matchedJetGeoIds(), chJetDet.matchedJetPtIds()); hasAcceptedPartMatch = !both.empty(); } @@ -2322,8 +2376,9 @@ struct RecoilJets { spectra.fill(HIST("hEventSelectionCountQA"), 3.5); } - if (!(isGoodZvtxFT0vsPV && isNoSameBunchPileup && isNoCollInTimeRangeStandard)) + if (!(isGoodZvtxFT0vsPV && isNoSameBunchPileup && isNoCollInTimeRangeStandard)) { return; + } spectra.fill(HIST("hEventSelectionCountQA"), 4.5); // All accepted events after 4 flags cut @@ -2387,21 +2442,23 @@ struct RecoilJets { spectra.fill(HIST("hScaledFT0MsignalWithoutBC"), scaledFT0M, fillNumber, weight); } - if (collision.foundBCId() != -1) // -1 if does not + if (collision.foundBCId() != -1) { // -1 if does not spectra.fill(HIST("hIsFT0SignalComeFromCollPerRun"), 2.5, fillNumber, weight); - else + } else { spectra.fill(HIST("hIsFT0SignalComeFromCollPerRun"), 3.5, fillNumber, weight); + } } //============================================================================= // Di-hadron azimuthal correlation in raw and MC det. level (no weight; MB events) data //============================================================================= template - void fillLeadingAndAssociatedTracksTask(JCollision const& collision, JTracks const& tracks, float weight = 1.) + void fillLeadingAndAssociatedTracksTask(JCollision const& collision, + JTracks const& tracks, + float weight = 1.) { std::vector vPhiOfLeadingTracks; std::vector vPtOfLeadingTracks; - std::vector vPhiOfAssociatedTracks; const auto ft0Metrics = getFT0Metrics(collision, ft0a.mean, ft0c.mean); const auto scaledFT0C = ft0Metrics.scaledFT0C; @@ -2409,8 +2466,9 @@ struct RecoilJets { // Search for leading tracks for (const auto& track : tracks) { - if (skipTrack(track)) + if (skipTrack(track)) { continue; + } const auto trackPt = track.pt(); @@ -2432,8 +2490,9 @@ struct RecoilJets { spectra.fill(HIST("hScaledFT0M_NleadTracks"), scaledFT0M, 0.5, weight); for (const auto& track : tracks) { - if (skipTrack(track)) + if (skipTrack(track)) { continue; + } const auto trackPt = track.pt(); const auto trackPhi = track.phi(); @@ -2453,9 +2512,9 @@ struct RecoilJets { //============================================================================= // Background fluctuations in raw data and MC det. level - template + template void fillBkgdFluctuations(JCollision const& collision, - Jets const& jets, + ChJets const& chJets, JTracks const& tracks, float weight = 1.) { @@ -2480,8 +2539,9 @@ struct RecoilJets { uint64_t index = 0; for (const auto& track : tracks) { ++index; - if (skipTrack(track)) + if (skipTrack(track)) { continue; + } float dEta = std::pow(randomConeEta - track.eta(), 2); float dPhi = std::pow(RecoDecay::constrainAngle(randomConePhi - track.phi(), -constants::math::PI), 2); @@ -2505,15 +2565,15 @@ struct RecoilJets { // Avoid leading jet (JE jet reconstruction sorts jets by pT) // square of distance to accept RC placement in events with leading jet - float dMinR2 = std::pow(jet.radius + bkgd.randomConeR + bkgd.minDeltaRToJet, 2); + float dMinR2 = std::pow(jets.radius + bkgd.randomConeR + bkgd.minDeltaRToJet, 2); // max # of attempts to find a place for RC; to avoid possibility with infinite loop in While cycle const int maxAttempts = 15000; - if (jets.size() > 0) // at least 1 jet + if (chJets.size() > 0) // at least 1 jet { - float leadJetEta = jets.iteratorAt(0).eta(); - float leadJetPhi = jets.iteratorAt(0).phi(); + float leadJetEta = chJets.iteratorAt(0).eta(); + float leadJetPhi = chJets.iteratorAt(0).phi(); float dEtaLeadJet = std::pow(leadJetEta - randomConeEta, 2); float dPhiLeadJet = std::pow(RecoDecay::constrainAngle(leadJetPhi - randomConePhi, -constants::math::PI), 2); @@ -2535,8 +2595,9 @@ struct RecoilJets { if (isTherePlaceForRC) { randomConePt = 0.0; for (const auto& track : tracks) { - if (skipTrack(track)) + if (skipTrack(track)) { continue; + } float dEta = std::pow(randomConeEta - track.eta(), 2); float dPhi = std::pow(RecoDecay::constrainAngle(randomConePhi - track.phi(), -constants::math::PI), 2); @@ -2557,8 +2618,9 @@ struct RecoilJets { float perpConePt = 0.0; for (const auto& track : tracks) { - if (skipTrack(track)) + if (skipTrack(track)) { continue; + } float dEta = std::pow(perpConeEta - track.eta(), 2); float dPhi = std::pow(RecoDecay::constrainAngle(perpConePhi - track.phi(), -constants::math::PI), 2); @@ -2574,18 +2636,18 @@ struct RecoilJets { //---------------------------------------------------------- // Avoid leading and subleading jets - if (jets.size() > 1) // at least 2 jets in an event + if (chJets.size() > 1) // at least 2 jets in an event { // Leading jet - float leadJetEta = jets.iteratorAt(0).eta(); - float leadJetPhi = jets.iteratorAt(0).phi(); + float leadJetEta = chJets.iteratorAt(0).eta(); + float leadJetPhi = chJets.iteratorAt(0).phi(); float dEtaLeadJet = std::pow(leadJetEta - randomConeEta, 2); float dPhiLeadJet = std::pow(RecoDecay::constrainAngle(leadJetPhi - randomConePhi, -constants::math::PI), 2); // Subleading jet - float subleadJetEta = jets.iteratorAt(1).eta(); - float subleadJetPhi = jets.iteratorAt(1).phi(); + float subleadJetEta = chJets.iteratorAt(1).eta(); + float subleadJetPhi = chJets.iteratorAt(1).phi(); float dEtaSubleadJet = std::pow(subleadJetEta - randomConeEta, 2); float dPhiSubleadJet = std::pow(RecoDecay::constrainAngle(subleadJetPhi - randomConePhi, -constants::math::PI), 2); @@ -2610,8 +2672,9 @@ struct RecoilJets { if (isTherePlaceForRC) { randomConePt = 0.0; for (const auto& track : tracks) { - if (skipTrack(track)) + if (skipTrack(track)) { continue; + } float dEta = std::pow(randomConeEta - track.eta(), 2); float dPhi = std::pow(RecoDecay::constrainAngle(randomConePhi - track.phi(), -constants::math::PI), 2); @@ -2628,7 +2691,7 @@ struct RecoilJets { //---------------------------------------------------------- // Place cone perpendicular to TTSig candidate - if (vCandForTT.size() > 0) // at least 1 TT + if (!vCandForTT.empty()) // at least 1 TT { auto randIndexTrack = randGen->Integer(vCandForTT.size()); auto objTT = tracks.iteratorAt(vCandForTT[randIndexTrack]); @@ -2640,8 +2703,9 @@ struct RecoilJets { if (std::abs(perpTTConeEta) < (trk.etaCut - bkgd.randomConeR)) { float perpTTConePt = 0.0; for (const auto& track : tracks) { - if (skipTrack(track)) + if (skipTrack(track)) { continue; + } float dEta = std::pow(perpTTConeEta - track.eta(), 2); float dPhi = std::pow(RecoDecay::constrainAngle(perpTTConePhi - track.phi(), -constants::math::PI), 2); @@ -2658,9 +2722,9 @@ struct RecoilJets { } } - template + template void fillBkgdFluctuationsMCPartLevel(JCollision const& collision, - Jets const& jets, + ChJets const& chJets, JParticles const& particles, float weight = 1.) { @@ -2684,8 +2748,9 @@ struct RecoilJets { uint64_t index = 0; for (const auto& particle : particles) { - if (skipParticle(particle)) + if (skipParticle(particle)) { continue; + } float dEta = std::pow(randomConeEta - particle.eta(), 2); float dPhi = std::pow(RecoDecay::constrainAngle(randomConePhi - particle.phi(), -constants::math::PI), 2); @@ -2710,15 +2775,15 @@ struct RecoilJets { // Avoid leading jet (JE jet reconstruction sorts jets by pT) // square of distance to accept RC placement in events with leading jet - float dMinR2 = std::pow(jet.radius + bkgd.randomConeR + bkgd.minDeltaRToJet, 2); + float dMinR2 = std::pow(jets.radius + bkgd.randomConeR + bkgd.minDeltaRToJet, 2); // max # of attempts to find a place for RC; to avoid possibility with infinite loop in While cycle const int maxAttempts = 15000; - if (jets.size() > 0) // at least 1 jet + if (chJets.size() > 0) // at least 1 jet { - float leadJetEta = jets.iteratorAt(0).eta(); - float leadJetPhi = jets.iteratorAt(0).phi(); + float leadJetEta = chJets.iteratorAt(0).eta(); + float leadJetPhi = chJets.iteratorAt(0).phi(); float dEtaLeadJet = std::pow(leadJetEta - randomConeEta, 2); float dPhiLeadJet = std::pow(RecoDecay::constrainAngle(leadJetPhi - randomConePhi, -constants::math::PI), 2); @@ -2740,14 +2805,14 @@ struct RecoilJets { if (isTherePlaceForRC) { randomConePt = 0.0; for (const auto& particle : particles) { - if (skipParticle(particle)) + if (skipParticle(particle)) { continue; + } float dEta = std::pow(randomConeEta - particle.eta(), 2); float dPhi = std::pow(RecoDecay::constrainAngle(randomConePhi - particle.phi(), -constants::math::PI), 2); - if ((dEta + dPhi) < radiusRC2) // inside RC - { + if ((dEta + dPhi) < radiusRC2) { // inside RC randomConePt += particle.pt(); } } @@ -2762,8 +2827,9 @@ struct RecoilJets { float perpConePt = 0.0; for (const auto& particle : particles) { - if (skipParticle(particle)) + if (skipParticle(particle)) { continue; + } float dEta = std::pow(perpConeEta - particle.eta(), 2); float dPhi = std::pow(RecoDecay::constrainAngle(perpConePhi - particle.phi(), -constants::math::PI), 2); @@ -2779,30 +2845,30 @@ struct RecoilJets { //---------------------------------------------------------- // Avoid leading and subleading jets - if (jets.size() > 1) // at least 2 jets in an event + if (chJets.size() > 1) // at least 2 jets in an event { // Leading jet - float leadJetEta = jets.iteratorAt(0).eta(); - float leadJetPhi = jets.iteratorAt(0).phi(); + float leadJetEta = chJets.iteratorAt(0).eta(); + float leadJetPhi = chJets.iteratorAt(0).phi(); float dEtaLeadJet = std::pow(leadJetEta - randomConeEta, 2); float dPhiLeadJet = std::pow(RecoDecay::constrainAngle(leadJetPhi - randomConePhi, -constants::math::PI), 2); // Subleading jet - float subleadJetEta = jets.iteratorAt(1).eta(); - float subleadJetPhi = jets.iteratorAt(1).phi(); + float subleadJetEta = chJets.iteratorAt(1).eta(); + float subleadJetPhi = chJets.iteratorAt(1).phi(); float dEtaSubleadJet = std::pow(subleadJetEta - randomConeEta, 2); float dPhiSubleadJet = std::pow(RecoDecay::constrainAngle(subleadJetPhi - randomConePhi, -constants::math::PI), 2); // Try to add events with TTsig bool keepEventWithTT = false; - if (vCandForTT.size() > 0) // at least 1 TT + if (!vCandForTT.empty()) // at least 1 TT { auto randIndexParticle = randGen->Integer(vCandForTT.size()); auto objTT = particles.iteratorAt(vCandForTT[randIndexParticle]); // Skip events where TT is not a part of leading or subleading jets (mutlijet event, difficult to place RC and avoid hard jets) - if (isTrackInJet(jets.iteratorAt(0), objTT) || isTrackInJet(jets.iteratorAt(1), objTT)) { + if (isTrackInJet(chJets.iteratorAt(0), objTT) || isTrackInJet(chJets.iteratorAt(1), objTT)) { keepEventWithTT = true; } } @@ -2828,8 +2894,9 @@ struct RecoilJets { if (isTherePlaceForRC) { randomConePt = 0.0; for (const auto& particle : particles) { - if (skipParticle(particle)) + if (skipParticle(particle)) { continue; + } float dEta = std::pow(randomConeEta - particle.eta(), 2); float dPhi = std::pow(RecoDecay::constrainAngle(randomConePhi - particle.phi(), -constants::math::PI), 2); @@ -2864,8 +2931,9 @@ struct RecoilJets { bool bSigEv = false; auto dice = randGen->Rndm(); - if (dice < tt.fracSig) + if (dice < tt.fracSig) { bSigEv = true; + } float ptTTMin = 0.0, ptTTMax = 0.0; if (bSigEv) { @@ -2902,12 +2970,14 @@ struct RecoilJets { } // Skip if no TT - if (chosenTTPos < 0) + if (chosenTTPos < 0) { return; + } bool bHasAssocMcPart = tracks.iteratorAt(chosenTTPos).has_mcParticle(); - if (!bHasAssocMcPart) + if (!bHasAssocMcPart) { return; + } // No filter on Particles, it can be outside of |eta| acceptance auto particle = tracks.iteratorAt(chosenTTPos).template mcParticle_as(); @@ -2938,8 +3008,9 @@ struct RecoilJets { spectra.fill(HIST("hScaledFT0M_FractionOfPartTTSigSatisfCond"), scaledFT0M, 1.5, 1.5); } } - if (!bPartWithinEta) + if (!bPartWithinEta) { return; + } //============================================================================= // Fill histograms @@ -2987,8 +3058,9 @@ struct RecoilJets { // Numerator candidates: reconstructed tracks in this detector collision for (const auto& trackPerColl : tracks) { - if (skipTrack(trackPerColl)) + if (skipTrack(trackPerColl)) { continue; + } auto trackPt = trackPerColl.pt(); @@ -3033,8 +3105,9 @@ struct RecoilJets { float particlePt = particleInMcCollision.pt(); bool bKeepParticle = std::fabs(particleEta) < trk.etaCut && !skipParticle(particleInMcCollision); - if (!bKeepParticle) + if (!bKeepParticle) { continue; + } spectra.fill(HIST("hScaledMultFT0C_AllGeneratedParticlesPt"), scaledFT0C, particlePt, weight); spectra.fill(HIST("hScaledMultFT0M_AllGeneratedParticlesPt"), scaledFT0M, particlePt, weight); @@ -3049,32 +3122,34 @@ struct RecoilJets { //============================================================================= void processData(CollRhoDataIt const& collision, TrackTbl const& tracksPerColl, - JetsDataTbl const& jetsPerColl) + ChJetsDataTbl const& chJetsPerColl) { spectra.fill(HIST("hEventSelectionCountReco"), 0.5); - if (skipEvent(collision)) + if (skipEvent(collision)) { return; + } spectra.fill(HIST("hEventSelectionCountReco"), 1.5); // number of events selected for analysis - fillHistograms(collision, jetsPerColl, tracksPerColl); + fillHistograms(collision, chJetsPerColl, tracksPerColl); } PROCESS_SWITCH(RecoilJets, processData, "process raw data", true); //____________________ void processMCDetLevel(CollRhoDataIt const& collision, TrackTbl const& tracksPerColl, - JetsDetTbl const& jetsPerColl) + ChJetsDetTbl const& chJetsPerColl) { spectra.fill(HIST("hEventSelectionCountReco"), 0.5); - if (skipEvent(collision)) + if (skipEvent(collision)) { return; + } spectra.fill(HIST("hEventSelectionCountReco"), 1.5); spectra.fill(HIST("hEventSelectionCountReco"), 4.5); // number of events selected for analysis - fillHistograms(collision, jetsPerColl, tracksPerColl); + fillHistograms(collision, chJetsPerColl, tracksPerColl); } PROCESS_SWITCH(RecoilJets, processMCDetLevel, "process MC det. level data (no weight; MB events)", false); @@ -3082,11 +3157,12 @@ struct RecoilJets { void processMCDetLevelWeighted(CollRhoOutlierDetIt const& collision, aod::JetMcCollisions const&, TrackTbl const& tracksPerColl, - JetsDetTbl const& jetsPerColl) + ChJetsDetTbl const& chJetsPerColl) { spectra.fill(HIST("hEventSelectionCountReco"), 0.5); - if (skipEvent(collision)) + if (skipEvent(collision)) { return; + } spectra.fill(HIST("hEventSelectionCountReco"), 1.5); @@ -3102,14 +3178,14 @@ struct RecoilJets { spectra.fill(HIST("hEventSelectionCountReco"), 4.5); // number of events selected for analysis auto weight = collision.mcCollision().weight(); - fillHistograms(collision, jetsPerColl, tracksPerColl, weight); + fillHistograms(collision, chJetsPerColl, tracksPerColl, weight); } PROCESS_SWITCH(RecoilJets, processMCDetLevelWeighted, "process MC det. level data (weighted JJ)", false); //_____________________ void processMCPartLevel(CollCentRhoPartIt const& collision, PartTbl const& particlesPerColl, - JetsPartTbl const& jetsPerColl) + ChJetsPartTbl const& chJetsPerColl) { spectra.fill(HIST("hEventSelectionCountPartLevel"), 0.5); @@ -3119,14 +3195,14 @@ struct RecoilJets { } spectra.fill(HIST("hEventSelectionCountPartLevel"), 3.5); // number of events selected for analysis - fillHistogramsMCPartLevel(collision, jetsPerColl, particlesPerColl); + fillHistogramsMCPartLevel(collision, chJetsPerColl, particlesPerColl); } PROCESS_SWITCH(RecoilJets, processMCPartLevel, "process MC part. level data (no weight; MB events)", false); //_____________________________ void processMCPartLevelWeighted(CollCentRhoOutlierPartIt const& collision, PartTbl const& particlesPerColl, - JetsPartTbl const& jetsPerColl) + ChJetsPartTbl const& chJetsPerColl) { spectra.fill(HIST("hEventSelectionCountPartLevel"), 0.5); @@ -3144,7 +3220,7 @@ struct RecoilJets { auto weight = collision.weight(); spectra.fill(HIST("ptHat"), collision.ptHard(), weight); - fillHistogramsMCPartLevel(collision, jetsPerColl, particlesPerColl, weight); + fillHistogramsMCPartLevel(collision, chJetsPerColl, particlesPerColl, weight); } PROCESS_SWITCH(RecoilJets, processMCPartLevelWeighted, "process MC part. level data (weighted JJ)", false); @@ -3154,24 +3230,26 @@ struct RecoilJets { void processJetsGeoMatching(CollRhoDetIt const& collision, CollRhoPartTbl const&, TrackTbl const& tracksPerColl, - MatchedJetsDetToPartTbl const& mcDetJetsPerColl, - MatchedJetsPartToDetTbl const& mcPartJets) + MatchedChJetsDetToPartTbl const& mcChJetsDetPerColl, + MatchedChJetsPartToDetTbl const& mcChJetsPart) { - if (skipEvent(collision) || !collision.has_mcCollision()) + if (skipEvent(collision) || !collision.has_mcCollision()) { return; + } const auto mcColl = collision.mcCollision_as(); - if (std::fabs(mcColl.posZ()) > ev.vertexZCut) + if (std::fabs(mcColl.posZ()) > ev.vertexZCut) { return; + } const auto partLevelCollRho = mcColl.rho(); // Slice for mc part level jets associated to a given mcCollisionId const auto mcCollisionId = collision.mcCollisionId(); - const auto mcPartJetsPerMcCollision = mcPartJets.sliceBy(partJetsByMcCollision, mcCollisionId); // signature: (__column to slice___, __index__) + const auto mcChJetsPartPerMcCollision = mcChJetsPart.sliceBy(partJetsByMcCollision, mcCollisionId); // signature: (__column to slice___, __index__) - fillMatchedGeoHistograms(collision, mcPartJetsPerMcCollision, mcDetJetsPerColl, tracksPerColl, partLevelCollRho); + fillMatchedGeoHistograms(collision, mcChJetsPartPerMcCollision, mcChJetsDetPerColl, tracksPerColl, partLevelCollRho); } PROCESS_SWITCH(RecoilJets, processJetsGeoMatching, "process matching of MC jets using Geo criterion (no weight; MB events)", false); @@ -3179,24 +3257,26 @@ struct RecoilJets { void processJetsGeoPtMatching(CollRhoDetIt const& collision, CollRhoPartTbl const&, TrackTbl const& tracksPerColl, - MatchedJetsDetToPartTbl const& mcDetJetsPerColl, - MatchedJetsPartToDetTbl const& mcPartJets) + MatchedChJetsDetToPartTbl const& mcChJetsDetPerColl, + MatchedChJetsPartToDetTbl const& mcChJetsPart) { // Skip detector level collisions - if (skipEvent(collision) || !collision.has_mcCollision()) + if (skipEvent(collision) || !collision.has_mcCollision()) { return; + } const auto mcColl = collision.mcCollision_as(); - if (std::fabs(mcColl.posZ()) > ev.vertexZCut) + if (std::fabs(mcColl.posZ()) > ev.vertexZCut) { return; + } const auto partLevelCollRho = mcColl.rho(); // Slice for mc part level jets associated to a given mcCollisionId const auto mcCollisionId = collision.mcCollisionId(); - const auto mcPartJetsPerMcCollision = mcPartJets.sliceBy(partJetsByMcCollision, mcCollisionId); // signature: (__column to slice___, __index__) + const auto mcChJetsPartPerMcCollision = mcChJetsPart.sliceBy(partJetsByMcCollision, mcCollisionId); // signature: (__column to slice___, __index__) - fillMatchedGeoPtHistograms(collision, mcPartJetsPerMcCollision, mcDetJetsPerColl, tracksPerColl, partLevelCollRho); + fillMatchedGeoPtHistograms(collision, mcChJetsPartPerMcCollision, mcChJetsDetPerColl, tracksPerColl, partLevelCollRho); } PROCESS_SWITCH(RecoilJets, processJetsGeoPtMatching, "process matching of MC jets using Geo+Pt criteria (no weight; MB events)", false); @@ -3204,25 +3284,27 @@ struct RecoilJets { void processJetsGeoMatchingWeighted(CollRhoOutlierDetIt const& collision, CollRhoOutlierPartTbl const&, TrackTbl const& tracksPerColl, - MatchedJetsDetToPartTbl const& mcDetJetsPerColl, - MatchedJetsPartToDetTbl const& mcPartJets) + MatchedChJetsDetToPartTbl const& mcChJetsDetPerColl, + MatchedChJetsPartToDetTbl const& mcChJetsPart) { // Skip detector level collisions - if (skipEvent(collision) || collision.isOutlier() || !collision.has_mcCollision()) + if (skipEvent(collision) || collision.isOutlier() || !collision.has_mcCollision()) { return; + } const auto mcColl = collision.mcCollision_as(); - if (std::fabs(mcColl.posZ()) > ev.vertexZCut) + if (std::fabs(mcColl.posZ()) > ev.vertexZCut) { return; + } const auto partLevelCollRho = mcColl.rho(); const auto weight = mcColl.weight(); // Slice for mc part level jets associated to a given mcCollisionId const auto mcCollisionId = collision.mcCollisionId(); - const auto mcPartJetsPerMcCollision = mcPartJets.sliceBy(partJetsByMcCollision, mcCollisionId); // signature: (__column to slice___, __index__) + const auto mcChJetsPartPerMcCollision = mcChJetsPart.sliceBy(partJetsByMcCollision, mcCollisionId); // signature: (__column to slice___, __index__) - fillMatchedGeoHistograms(collision, mcPartJetsPerMcCollision, mcDetJetsPerColl, tracksPerColl, partLevelCollRho, weight); + fillMatchedGeoHistograms(collision, mcChJetsPartPerMcCollision, mcChJetsDetPerColl, tracksPerColl, partLevelCollRho, weight); } PROCESS_SWITCH(RecoilJets, processJetsGeoMatchingWeighted, "process matching of MC jets using Geo criterion (weighted JJ)", false); @@ -3230,25 +3312,27 @@ struct RecoilJets { void processJetsGeoPtMatchingWeighted(CollRhoOutlierDetIt const& collision, CollRhoOutlierPartTbl const&, TrackTbl const& tracksPerColl, - MatchedJetsDetToPartTbl const& mcDetJetsPerColl, - MatchedJetsPartToDetTbl const& mcPartJets) + MatchedChJetsDetToPartTbl const& mcChJetsDetPerColl, + MatchedChJetsPartToDetTbl const& mcChJetsPart) { // Skip detector level collisions - if (skipEvent(collision) || collision.isOutlier() || !collision.has_mcCollision()) + if (skipEvent(collision) || collision.isOutlier() || !collision.has_mcCollision()) { return; + } auto mcColl = collision.mcCollision_as(); - if (std::fabs(mcColl.posZ()) > ev.vertexZCut) + if (std::fabs(mcColl.posZ()) > ev.vertexZCut) { return; + } auto partLevelCollRho = mcColl.rho(); auto weight = mcColl.weight(); // Slice for mc part level jets associated to a given mcCollisionId auto mcCollisionId = collision.mcCollisionId(); - auto mcPartJetsPerMcCollision = mcPartJets.sliceBy(partJetsByMcCollision, mcCollisionId); // signature: (__column to slice___, __index__) + auto mcChJetsPartPerMcCollision = mcChJetsPart.sliceBy(partJetsByMcCollision, mcCollisionId); // signature: (__column to slice___, __index__) - fillMatchedGeoPtHistograms(collision, mcPartJetsPerMcCollision, mcDetJetsPerColl, tracksPerColl, partLevelCollRho, weight); + fillMatchedGeoPtHistograms(collision, mcChJetsPartPerMcCollision, mcChJetsDetPerColl, tracksPerColl, partLevelCollRho, weight); } PROCESS_SWITCH(RecoilJets, processJetsGeoPtMatchingWeighted, "process matching of MC jets using Geo+Pt criteria (weighted JJ)", false); @@ -3257,8 +3341,9 @@ struct RecoilJets { //============================================================================= void processEventActivity(EvMultZDCDataIt const& collision) { - if (skipEvent(collision)) + if (skipEvent(collision)) { return; + } fillMultiplicityHistograms(collision); } @@ -3268,8 +3353,9 @@ struct RecoilJets { void processEventActivityMCDetLevelWeighted(EvMultOutlierZDCDetIt const& collision, aod::JetMcCollisions const&) { - if (skipEvent(collision) || collision.isOutlier() || !collision.has_mcCollision()) + if (skipEvent(collision) || collision.isOutlier() || !collision.has_mcCollision()) { return; + } auto weight = collision.mcCollision().weight(); fillMultiplicityHistograms(collision, weight); @@ -3281,8 +3367,9 @@ struct RecoilJets { //============================================================================= void processEventActivityMCPartLevel(CollPartIt const& collision) { - if (skipMCEvent(collision)) + if (skipMCEvent(collision)) { return; + } fillMultiplicityHistogramsMCPartLevel(collision); } @@ -3291,8 +3378,9 @@ struct RecoilJets { //__________________________________________ void processEventActivityMCPartLevelWeighted(EvMultOutlierPartIt const& collision) { - if (skipMCEvent(collision) || collision.isOutlier()) + if (skipMCEvent(collision) || collision.isOutlier()) { return; + } auto weight = collision.weight(); fillMultiplicityHistogramsMCPartLevel(collision, weight); @@ -3307,8 +3395,9 @@ struct RecoilJets { aod::Zdcs const&) { // Base flag for event selection - if (!collision.sel8()) + if (!collision.sel8()) { return; + } fillEventActivitySelectionQAHistograms(collision); } @@ -3320,8 +3409,10 @@ struct RecoilJets { void processLeadingAndAssociatedTracksTask(CollDataIt const& collision, TrackTbl const& tracksPerColl) { - if (skipEvent(collision)) + if (skipEvent(collision)) { return; + } + fillLeadingAndAssociatedTracksTask(collision, tracksPerColl); } PROCESS_SWITCH(RecoilJets, processLeadingAndAssociatedTracksTask, "process di-hadron azimuthal correlation in raw and MC det. level (no weight; MB events) data", false); @@ -3331,24 +3422,26 @@ struct RecoilJets { //============================================================================= void processBkgdFluctuations(CollRhoDataIt const& collision, TrackTbl const& tracksPerColl, - JetsDataTbl const& jetsPerColl) + ChJetsDataTbl const& chJetsPerColl) { - if (skipEvent(collision)) + if (skipEvent(collision)) { return; + } - fillBkgdFluctuations(collision, jetsPerColl, tracksPerColl); + fillBkgdFluctuations(collision, chJetsPerColl, tracksPerColl); } PROCESS_SWITCH(RecoilJets, processBkgdFluctuations, "process raw data to estimate bkgd fluctuations", false); //____________________________________ void processBkgdFluctuationsMCDetLevel(CollRhoDataIt const& collision, TrackTbl const& tracksPerColl, - JetsDetTbl const& jetsPerColl) + ChJetsDetTbl const& chJetsPerColl) { - if (skipEvent(collision)) + if (skipEvent(collision)) { return; + } - fillBkgdFluctuations(collision, jetsPerColl, tracksPerColl); + fillBkgdFluctuations(collision, chJetsPerColl, tracksPerColl); } PROCESS_SWITCH(RecoilJets, processBkgdFluctuationsMCDetLevel, "process MC det. level (no weight; MB events) data to estimate bkgd fluctuations", false); @@ -3356,38 +3449,41 @@ struct RecoilJets { void processBkgdFluctuationsMCDetLevelWeighted(CollRhoOutlierDetIt const& collision, aod::JetMcCollisions const&, TrackTbl const& tracksPerColl, - JetsDetTbl const& jetsPerColl) + ChJetsDetTbl const& chJetsPerColl) { - if (skipEvent(collision) || collision.isOutlier() || !collision.has_mcCollision()) + if (skipEvent(collision) || collision.isOutlier() || !collision.has_mcCollision()) { return; + } auto weight = collision.mcCollision().weight(); - fillBkgdFluctuations(collision, jetsPerColl, tracksPerColl, weight); + fillBkgdFluctuations(collision, chJetsPerColl, tracksPerColl, weight); } PROCESS_SWITCH(RecoilJets, processBkgdFluctuationsMCDetLevelWeighted, "process MC det. level (weighted JJ) data to estimate bkgd fluctuations", false); //_____________________________________ void processBkgdFluctuationsMCPartLevel(CollRhoPartIt const& collision, PartTbl const& particlesPerColl, - JetsPartTbl const& jetsPerColl) + ChJetsPartTbl const& chJetsPerColl) { - if (skipMCEvent(collision)) + if (skipMCEvent(collision)) { return; + } - fillBkgdFluctuationsMCPartLevel(collision, jetsPerColl, particlesPerColl); + fillBkgdFluctuationsMCPartLevel(collision, chJetsPerColl, particlesPerColl); } PROCESS_SWITCH(RecoilJets, processBkgdFluctuationsMCPartLevel, "process MC part. level (no weight; MB events) data to estimate bkgd fluctuations", false); //_____________________________________________ void processBkgdFluctuationsMCPartLevelWeighted(CollRhoOutlierPartIt const& collision, PartTbl const& particlesPerColl, - JetsPartTbl const& jetsPerColl) + ChJetsPartTbl const& chJetsPerColl) { - if (skipMCEvent(collision) || collision.isOutlier()) + if (skipMCEvent(collision) || collision.isOutlier()) { return; + } auto weight = collision.weight(); - fillBkgdFluctuationsMCPartLevel(collision, jetsPerColl, particlesPerColl, weight); + fillBkgdFluctuationsMCPartLevel(collision, chJetsPerColl, particlesPerColl, weight); } PROCESS_SWITCH(RecoilJets, processBkgdFluctuationsMCPartLevelWeighted, "process MC part. level (weighted JJ) data to estimate bkgd fluctuations", false); @@ -3400,8 +3496,9 @@ struct RecoilJets { { // Skip detector level collisions - if (skipEvent(collision)) + if (skipEvent(collision)) { return; + } fillTTSmearingPtPhi(collision, tracksPerColl); } @@ -3415,8 +3512,9 @@ struct RecoilJets { { // Skip detector level collisions - if (skipEvent(collision) || collision.isOutlier() || !collision.has_mcCollision()) + if (skipEvent(collision) || collision.isOutlier() || !collision.has_mcCollision()) { return; + } auto weight = collision.mcCollision().weight(); fillTTSmearingPtPhi(collision, tracksPerColl, weight); @@ -3430,8 +3528,9 @@ struct RecoilJets { aod::JetParticles const& particles) { // Skip detector level collision - if (skipEvent(collision) || !collision.has_mcCollision()) + if (skipEvent(collision) || !collision.has_mcCollision()) { return; + } fillTrackRecoEff(collision, tracksPerColl, particles); } @@ -3444,8 +3543,9 @@ struct RecoilJets { aod::JetParticles const& particles) { // Skip detector level collision - if (skipEvent(collision) || !collision.has_mcCollision() || collision.isOutlier()) + if (skipEvent(collision) || !collision.has_mcCollision() || collision.isOutlier()) { return; + } const auto weight = collision.mcCollision().weight(); fillTrackRecoEff(collision, tracksPerColl, particles, weight); @@ -3485,9 +3585,9 @@ struct RecoilJets { } template - std::tuple isRecoilJet(const Jet& jet, double phiTT) + std::tuple isRecoilJet(const Jet& chJet, double phiTT) { - double dphi = std::fabs(RecoDecay::constrainAngle(jet.phi() - phiTT, -constants::math::PI)); + double dphi = std::fabs(RecoDecay::constrainAngle(chJet.phi() - phiTT, -constants::math::PI)); return {dphi, (constants::math::PI - tt.recoilRegion) < dphi}; } @@ -3538,7 +3638,7 @@ struct RecoilJets { { bool bIsJetWithHighPtConstituent = false; for (const auto& chJetConstituent : chJet.template tracks_as()) { - if (chJetConstituent.pt() > jet.constituentPtMax) { + if (chJetConstituent.pt() > jets.constituentPtMax) { bIsJetWithHighPtConstituent = true; break; } @@ -3547,9 +3647,9 @@ struct RecoilJets { } template - bool isTrackInJet(Jet const& jet, Track const& track) + bool isTrackInJet(Jet const& chJet, Track const& track) { - for (auto const& constituentId : jet.tracksIds()) { + for (auto const& constituentId : chJet.tracksIds()) { if (constituentId == track.globalIndex()) { return true; } @@ -3571,8 +3671,9 @@ struct RecoilJets { } } - if (binNumber == -1) // No bin found + if (binNumber == -1) { // No bin found return 0; + } return binNumber; } @@ -3583,8 +3684,9 @@ struct RecoilJets { std::vector vPhiOfTT; for (const auto& track : tracks) { - if (skipTrack(track)) + if (skipTrack(track)) { continue; + } // Search for TT_Sig candidate const auto ptTTsigMin = tt.sigPtRange->at(0); @@ -3609,8 +3711,9 @@ struct RecoilJets { const int dimOf2DHist = 2; bool isHist2D = histPointer->GetDimension() == dimOf2DHist; - if (isHist2D) + if (isHist2D) { return; + } for (int iFlag = 0; iFlag < nEvSelFlags; ++iFlag) { histPointer->GetZaxis()->SetBinLabel(iFlag + 1, zAxis[iFlag]); @@ -3637,6 +3740,42 @@ struct RecoilJets { } return bothSet; } + + // Check whether scaled FT0 belongs to the requested EA interval + bool isInEAInterval(float scaledFT0, + const std::vector& ft0Edges, + EAInterval eaInterval) + { + switch (eaInterval) { + case kMB: + return true; + + case kEA0_20: + return scaledFT0 >= ft0Edges[kEA20] && + scaledFT0 < ft0Edges[kEA0]; + + case kEA0_10: + return scaledFT0 >= ft0Edges[kEA10] && + scaledFT0 < ft0Edges[kEA0]; + + case kEA20_40: + return scaledFT0 >= ft0Edges[kEA40] && + scaledFT0 < ft0Edges[kEA20]; + + case kEA60_80: + return scaledFT0 >= ft0Edges[kEA80] && + scaledFT0 < ft0Edges[kEA60]; + + case kEA50_100: + return scaledFT0 >= ft0Edges[kEA100] && + scaledFT0 < ft0Edges[kEA50]; + + case kEA80_100: + return scaledFT0 >= ft0Edges[kEA100] && + scaledFT0 < ft0Edges[kEA80]; + } + return false; + } }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)