diff --git a/PWGJE/Tasks/trackJetSpectra.cxx b/PWGJE/Tasks/trackJetSpectra.cxx index 7331e20c985..c0248a74c7a 100644 --- a/PWGJE/Tasks/trackJetSpectra.cxx +++ b/PWGJE/Tasks/trackJetSpectra.cxx @@ -32,12 +32,12 @@ #include #include #include -#include #include "TRandom3.h" +#include #include -#include #include +#include #include #include @@ -47,8 +47,6 @@ #include #include - - using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; @@ -79,12 +77,6 @@ using ParticlesFullTable = soa::Join; using TracksMCDetFullTable = soa::Join; - - - - - - using FilteredTracks = soa::Filtered; using FilteredMatchedJetsDetLevel = soa::Filtered>; @@ -117,27 +109,27 @@ struct TrackJetSpectra { Configurable dPhiCut{"dPhiCut", 0.6, ""}; // List of configurable parameters for histograms Configurable histJetPt{"histJetPt", 150, "Maximum value of jet pT shown in histograms"}; - + Configurable meanFT0A{"meanFT0A", -1., "Mean value of FT0A signal"}; Configurable meanFT0C{"meanFT0C", -1., "Mean value of FT0C signal"}; - + Configurable pionMass{"pionMass", o2::constants::physics::MassPiPlus, "Mass of a pion"}; - + Configurable isMCJJProd{"isMCJJProd", false, "Flag to select MC production: MB (false) or JJ(true)"}, - skipMBGapEvents{"skipMBGapEvents", false, - "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"}; + skipMBGapEvents{"skipMBGapEvents", false, + "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"}; //------------------------------------------------------------ - + TRandom3* rand = new TRandom3(0); // Declare filter on collision Z vertex Filter collisionFilter = nabs(aod::jcollision::posZ) < vertexZCut; Filter collisionFilterMC = nabs(aod::jmccollision::posZ) < vertexZCut; // Declare filters on accepted tracks and MC particles (settings for jet reco are provided in the jet finder wagon) - Filter trackFilter = aod::jtrack::pt > trkPtMin && aod::jtrack::pt < trkPtMax && nabs(aod::jtrack::eta) < trkEtaCut; + Filter trackFilter = aod::jtrack::pt > trkPtMin&& aod::jtrack::pt < trkPtMax&& nabs(aod::jtrack::eta) < trkEtaCut; Filter partFilter = nabs(aod::jmcparticle::eta) < trkEtaCut; // Declare filter on jets @@ -147,23 +139,16 @@ struct TrackJetSpectra { std::vector eventSelectionBits; int trackSelection = -1; - + Service pdg; Preslice partJetsPerCollision = aod::jet::mcCollisionId; - - Preslice tracksPerJCollision = o2::aod::jtrack::collisionId; - - - - - - + Preslice tracksPerJCollision = o2::aod::jtrack::collisionId; void init(InitContext const&) { // Initialize histogram axes - AxisSpec pT{histJetPt+20, -20.0, histJetPt * 1., "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec pT{histJetPt + 20, -20.0, histJetPt * 1., "#it{p}_{T} (GeV/#it{c})"}; AxisSpec pTRes{100, -5.0, 5.0, "#it{p}_{T} Resolution"}; AxisSpec pTTrack{100, 0.0, 100.0, "#it{p}_{T,Track} (GeV/#it{c})"}; AxisSpec phiAngle{40, 0.0, constants::math::TwoPI, "#it{#varphi} (rad)"}; @@ -171,300 +156,274 @@ struct TrackJetSpectra { AxisSpec etaTracks{100, -1.0, 1.0, "#it{#eta}_{trk}"}; AxisSpec detaTracks{400, -0.9, 0.9, "#it{#eta}_{trk}"}; AxisSpec etaJets{100, -0.6, 0.6, "#it{#eta}_{jet}"}; - - + AxisSpec jetArea{50, 0.0, 3., "Area_{jet}"}; AxisSpec rho{50, 0.0, 50., "#it{#rho}"}; AxisSpec rhoArea{60, 0.0, 60., "#it{#rho} #times Area_{jet}"}; - + AxisSpec centrality{10, 0.0, 100., "centrality"}; - + // Convert configurable strings to std::string std::string evSelToString = static_cast(evSel); std::string trkSelToString = static_cast(trkSel); eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(evSelToString); trackSelection = jetderiveddatautilities::initialiseTrackSelection(trkSelToString); - - if(doprocessFilteredCollisions) - { - // Event selection - spectra.add("hEventSelectionCount", "Count # of events in the analysis", kTH1F, {{3, 0.0, 3.}}); - spectra.get(HIST("hEventSelectionCount"))->GetXaxis()->SetBinLabel(1, "Total # of events w/o cuts"); - spectra.get(HIST("hEventSelectionCount"))->GetXaxis()->SetBinLabel(2, Form("# of events after sel. %s", evSelToString.data())); - spectra.get(HIST("hEventSelectionCount"))->GetXaxis()->SetBinLabel(3, "# of events after z cut"); //? - - // Z coordinate of collision vertex - spectra.add("hVertexZ_NoCut", "z vertex of collisions w/o cut", kTH1F, {{100, -20., 20., "#it{z}_{vertex}"}}); - spectra.add("hVertexZ_Cut", "z vertex of collisions w. cut", kTH1F, {{100, -20., 20., "#it{z}_{vertex}"}}); - spectra.add("hVertexZ_EventFiltering", "z vertex of collisions w. event filtering", kTH1F, {{100, -20., 20., "#it{z}_{vertex}"}}); + + if (doprocessFilteredCollisions) { + // Event selection + spectra.add("hEventSelectionCount", "Count # of events in the analysis", kTH1F, {{3, 0.0, 3.}}); + spectra.get(HIST("hEventSelectionCount"))->GetXaxis()->SetBinLabel(1, "Total # of events w/o cuts"); + spectra.get(HIST("hEventSelectionCount"))->GetXaxis()->SetBinLabel(2, Form("# of events after sel. %s", evSelToString.data())); + spectra.get(HIST("hEventSelectionCount"))->GetXaxis()->SetBinLabel(3, "# of events after z cut"); //? + + // Z coordinate of collision vertex + spectra.add("hVertexZ_NoCut", "z vertex of collisions w/o cut", kTH1F, {{100, -20., 20., "#it{z}_{vertex}"}}); + spectra.add("hVertexZ_Cut", "z vertex of collisions w. cut", kTH1F, {{100, -20., 20., "#it{z}_{vertex}"}}); + spectra.add("hVertexZ_EventFiltering", "z vertex of collisions w. event filtering", kTH1F, {{100, -20., 20., "#it{z}_{vertex}"}}); } // Disitribution of tracks - if(doprocessJetsMCDet||doprocessJetsMCDetWeighted||doprocessJets) - { - spectra.add("hPocetTracku", "Pocet", kTH1F, {{1000, 0., 1000., "Pocet Eventu"}}); - spectra.add("hPocetSignalTriggeru", "Pocet", kTH1F, {{10, 0., 10., "Pocet Eventu"}}); - spectra.add("hPocetReferenceTriggeru", "Pocet", kTH1F, {{10, 0., 10., "Pocet Eventu"}}); - spectra.add("hRecoilJetRefPt", "#it{p}_{T} distribution of recoil jets", kTH2F, {{200, 0., 200., "Pocet Eventu"}, centrality}); - spectra.add("hRecoilJetSigPt", "#it{p}_{T} distribution of recoil jets", kTH2F, {{200, 0., 200., "Pocet Eventu"}, centrality}); - spectra.add("hRecoilJetRefCorrPt", "Corrected #it{p}_{T} distribution of recoil jets", kTH2F, {{200, 0., 200., "Pocet Eventu"}, centrality}); - spectra.add("hRecoilJetSigCorrPt", "Corrected #it{p}_{T} distribution of recoil jets", kTH2F, {{200, 0., 200., "Pocet Eventu"}, centrality}); - spectra.add("hTrackPt", "#it{p}_{T} distribution of tracks", kTH2F, {pT, centrality}); - spectra.add("hTrackPhi", "#varphi distribution of tracks", kTH1F, {phiAngle}); - spectra.add("hTrackEta", "#eta distribution of tracks", kTH1F, {etaTracks}); - - spectra.add("hTrackPtPhi", "#it{p}_{T} vs. #varphi distribution of tracks", kTH2F, {pT, phiAngle}); - spectra.add("hTrackPtEta", "#it{p}_{T} vs. #eta distribution of tracks", kTH2F, {pT, etaTracks}); - - spectra.add("hTrackPtPhiEta", "#it{p}_{T} vs. #varphi vs. #eta distribution of tracks", kTH3F, {pT, phiAngle, etaTracks}); - - // Distribution of jets - spectra.add("hJetPt", "#it{p}_{T} distribution of jets", kTH1F, {pT}); - spectra.add("hJetPhi", "#varphi distribution of jets", kTH1F, {phiAngle}); - spectra.add("hJetEta", "#eta distribution of jets", kTH1F, {etaJets}); - - spectra.add("hJetPtPhi", "#it{p}_{T} vs. #varphi distribution of jets", kTH2F, {pT, phiAngle}); - spectra.add("hJetPtEta", "#it{p}_{T} vs. #eta distribution of jets", kTH2F, {pT, etaJets}); - - spectra.add("hJetPtPhiEta", "#it{p}_{T} vs. #varphi vs. #eta distribution of jets", kTH3F, {pT, phiAngle, etaJets}); - - spectra.add("hMultFT0A", "Mult. signal from FTOA", kTH2F, {{2000, 0.0, 40000., "FT0A"}, centrality}); - spectra.add("hMultFT0C", "Mult. signal from FTOC", kTH2F, {{2000, 0.0, 40000., "FT0C"}, centrality}); - - spectra.add("hJetPtCorr", "#it{p}_{T} distribution of jets", kTH1F, {pT}); - - spectra.add("hRho", "distribution of rho", kTH1F, {rho}); - spectra.add("hjetArea", "distribution of jet area", kTH1F, {jetArea}); - - spectra.add("JetAreavsPt", "#it{p}_{T} vs. jet area distribution of jets", kTH2F, {pT, jetArea}); - spectra.add("hRhovsNTracks", "Rho vs number of tracks", kTH2F, {{2000, 0.0, 40000., "NTracks"}, rho}); - spectra.add("hRhovsCentrality", "Rho vs number of tracks", kTH2F, {centrality, rho}); - - spectra.add("hSigCorrConePt", "distribution of Signal Corrected perpendicular Cone Pt", kTH1F, {pT}); - spectra.add("hRefCorrConePt", "distribution of Reference Corrected perpendicular Cone Pt", kTH1F, {pT}); - - spectra.add("hSigRandCorrConePt", "distribution of Signal Corrected random Cone Pt", kTH2F, {pT, centrality}); - spectra.add("hRefRandCorrConePt", "distribution of Reference Corrected random Cone Pt", kTH2F, {pT, centrality}); - - spectra.add("hTTCountData", "Pocet", kTH2F, {{2, 0., 2., "Pocet TT Data"}, centrality}); - } - - if(doprocessJetsMatched) - { - spectra.add("hMissedJets_pT", "Part. level jets w/o matched pair", kTH1F, {{200, 0.0, 200.}}, setSumw2); - spectra.add("hMissedJets_pT_RecoilJets", "Part. level jets w/o matched pair", kTH1F, {{200, 0.0, 200.}}, setSumw2); - - spectra.add("hJetPt_resolution", "Jet p_{T} relative resolution as a func. of jet #it{p}_{T, part}", kTH2F, {{100, -5., 5.}, pT}, setSumw2); - spectra.add("hJetPhi_resolution", "#varphi resolution as a func. of jet #it{p}_{T, part}", kTH2F, {{40, -1., 1.}, pT}, setSumw2); - - spectra.add("hJetPt_DetLevel_vs_PartLevel_RecoilJets", "Correlation recoil jet pT at part. vs. det. levels", kTH2F, {{200, 0.0, 200.}, {200, 0.0, 200.}}, setSumw2); - spectra.add("hJetPt_resolution_RecoilJets", "Jet p_{T} relative resolution as a func. of jet #it{p}_{T, part}", kTH2F, {{100, -5., 5.}, pT}, setSumw2); - spectra.add("hJetPhi_resolution_RecoilJets", "#varphi resolution as a func. of jet #it{p}_{T, part}", kTH2F, {{40, -1., 1.}, pT}, setSumw2); - - + if (doprocessJetsMCDet || doprocessJetsMCDetWeighted || doprocessJets) { + spectra.add("hPocetTracku", "Pocet", kTH1F, {{1000, 0., 1000., "Pocet Eventu"}}); + spectra.add("hPocetSignalTriggeru", "Pocet", kTH1F, {{10, 0., 10., "Pocet Eventu"}}); + spectra.add("hPocetReferenceTriggeru", "Pocet", kTH1F, {{10, 0., 10., "Pocet Eventu"}}); + spectra.add("hRecoilJetRefPt", "#it{p}_{T} distribution of recoil jets", kTH2F, {{200, 0., 200., "Pocet Eventu"}, centrality}); + spectra.add("hRecoilJetSigPt", "#it{p}_{T} distribution of recoil jets", kTH2F, {{200, 0., 200., "Pocet Eventu"}, centrality}); + spectra.add("hRecoilJetRefCorrPt", "Corrected #it{p}_{T} distribution of recoil jets", kTH2F, {{200, 0., 200., "Pocet Eventu"}, centrality}); + spectra.add("hRecoilJetSigCorrPt", "Corrected #it{p}_{T} distribution of recoil jets", kTH2F, {{200, 0., 200., "Pocet Eventu"}, centrality}); + spectra.add("hTrackPt", "#it{p}_{T} distribution of tracks", kTH2F, {pT, centrality}); + spectra.add("hTrackPhi", "#varphi distribution of tracks", kTH1F, {phiAngle}); + spectra.add("hTrackEta", "#eta distribution of tracks", kTH1F, {etaTracks}); + + spectra.add("hTrackPtPhi", "#it{p}_{T} vs. #varphi distribution of tracks", kTH2F, {pT, phiAngle}); + spectra.add("hTrackPtEta", "#it{p}_{T} vs. #eta distribution of tracks", kTH2F, {pT, etaTracks}); + + spectra.add("hTrackPtPhiEta", "#it{p}_{T} vs. #varphi vs. #eta distribution of tracks", kTH3F, {pT, phiAngle, etaTracks}); + + // Distribution of jets + spectra.add("hJetPt", "#it{p}_{T} distribution of jets", kTH1F, {pT}); + spectra.add("hJetPhi", "#varphi distribution of jets", kTH1F, {phiAngle}); + spectra.add("hJetEta", "#eta distribution of jets", kTH1F, {etaJets}); + + spectra.add("hJetPtPhi", "#it{p}_{T} vs. #varphi distribution of jets", kTH2F, {pT, phiAngle}); + spectra.add("hJetPtEta", "#it{p}_{T} vs. #eta distribution of jets", kTH2F, {pT, etaJets}); + + spectra.add("hJetPtPhiEta", "#it{p}_{T} vs. #varphi vs. #eta distribution of jets", kTH3F, {pT, phiAngle, etaJets}); + + spectra.add("hMultFT0A", "Mult. signal from FTOA", kTH2F, {{2000, 0.0, 40000., "FT0A"}, centrality}); + spectra.add("hMultFT0C", "Mult. signal from FTOC", kTH2F, {{2000, 0.0, 40000., "FT0C"}, centrality}); + + spectra.add("hJetPtCorr", "#it{p}_{T} distribution of jets", kTH1F, {pT}); + + spectra.add("hRho", "distribution of rho", kTH1F, {rho}); + spectra.add("hjetArea", "distribution of jet area", kTH1F, {jetArea}); + + spectra.add("JetAreavsPt", "#it{p}_{T} vs. jet area distribution of jets", kTH2F, {pT, jetArea}); + spectra.add("hRhovsNTracks", "Rho vs number of tracks", kTH2F, {{2000, 0.0, 40000., "NTracks"}, rho}); + spectra.add("hRhovsCentrality", "Rho vs number of tracks", kTH2F, {centrality, rho}); + + spectra.add("hSigCorrConePt", "distribution of Signal Corrected perpendicular Cone Pt", kTH1F, {pT}); + spectra.add("hRefCorrConePt", "distribution of Reference Corrected perpendicular Cone Pt", kTH1F, {pT}); + + spectra.add("hSigRandCorrConePt", "distribution of Signal Corrected random Cone Pt", kTH2F, {pT, centrality}); + spectra.add("hRefRandCorrConePt", "distribution of Reference Corrected random Cone Pt", kTH2F, {pT, centrality}); + + spectra.add("hTTCountData", "Pocet", kTH2F, {{2, 0., 2., "Pocet TT Data"}, centrality}); } - - if(doprocessMCPartLevel||doprocessMCPartLevelWeighted) - { - spectra.add("hRecoilJetRefPtMCP", "#it{p}_{T} distribution of recoil jets", kTH1F, {{200, 0., 200., "p_{T,ch jet} (GeV/c)"}}); - spectra.add("hRecoilJetSigPtMCP", "#it{p}_{T} distribution of recoil jets", kTH1F, {{200, 0., 200., "p_{T,ch jet} (GeV/c)"}}); - spectra.add("hRecoilJetRefCorrPtMCP", "Corrected #it{p}_{T} distribution of recoil jets", kTH1F, {{200, 0., 200., "p_{T,ch jet} (GeV/c)"}}); - spectra.add("hRecoilJetSigCorrPtMCP", "Corrected #it{p}_{T} distribution of recoil jets", kTH1F, {{200, 0., 200., "p_{T,ch jet} (GeV/c)"}}); - - spectra.add("hEventSelectionCountPartLevel", "Pocet", kTH1F, {{10, 0., 10., "Pocet Eventu"}}); - spectra.add("hMCPTrackPt", "#it{p}_{T} distribution of tracks", kTH1F, {pT}); - - spectra.add("hJetPtEtaPhiRhoArea_Part", "Charact. of inclusive part. level jets", kTHnSparseF, {pT, etaJets, phiAngle, rho, jetArea}); - spectra.add("hJetPtMCP", "#it{p}_{T} distribution of jets", kTH1F, {{200, 0., 200., "p_{T,ch jet} (GeV/c)"}}); - - spectra.add("hTTCountPartLevel", "Pocet", kTH1F, {{2, 0., 2., "Pocet TT MCP"}}); - + + if (doprocessJetsMatched) { + spectra.add("hMissedJets_pT", "Part. level jets w/o matched pair", kTH1F, {{200, 0.0, 200.}}, setSumw2); + spectra.add("hMissedJets_pT_RecoilJets", "Part. level jets w/o matched pair", kTH1F, {{200, 0.0, 200.}}, setSumw2); + + spectra.add("hJetPt_resolution", "Jet p_{T} relative resolution as a func. of jet #it{p}_{T, part}", kTH2F, {{100, -5., 5.}, pT}, setSumw2); + spectra.add("hJetPhi_resolution", "#varphi resolution as a func. of jet #it{p}_{T, part}", kTH2F, {{40, -1., 1.}, pT}, setSumw2); + + spectra.add("hJetPt_DetLevel_vs_PartLevel_RecoilJets", "Correlation recoil jet pT at part. vs. det. levels", kTH2F, {{200, 0.0, 200.}, {200, 0.0, 200.}}, setSumw2); + spectra.add("hJetPt_resolution_RecoilJets", "Jet p_{T} relative resolution as a func. of jet #it{p}_{T, part}", kTH2F, {{100, -5., 5.}, pT}, setSumw2); + spectra.add("hJetPhi_resolution_RecoilJets", "#varphi resolution as a func. of jet #it{p}_{T, part}", kTH2F, {{40, -1., 1.}, pT}, setSumw2); } - - - if(doprocessTrackEfficiency||doprocessTrackEfficiencyMB) - { - spectra.add("hTrackingEfficiencyDenominatorPt", "#it{p}_{T} distribution", kTH2F, {{200, 0., 200., "p_{T,ch} (GeV/c)"}, centrality}); - spectra.add("hNonprimaryParticlesPt", "#it{p}_{T} distribution", kTH1F, {{200, 0., 200., "p_{T,ch} (GeV/c)"}}); - spectra.add("hTrackingEfficiencyNumeratorPt", "#it{p}_{T} distribution", kTH2F, {{200, 0., 200., "p_{T,ch} (GeV/c)"}, centrality}); - - spectra.add("hTrackingPhiSmearingPt", "#it{p}_{T} distribution", kTH2F, {pT, dphiAngle}); - spectra.add("hTrackingEtaSmearingPt", "#it{p}_{T} distribution", kTH2F, {pT, detaTracks}); - spectra.add("hTrackingPtResolutionPt", "#it{p}_{T} distribution", kTH2F, {pTTrack, pTRes}); - - spectra.add("hNTracksPerCollision", "Pocet", kTH1F, {{1000, 0., 1000., "Pocet Eventu"}}); - - spectra.add("CentralityC", "#it{p}_{T} distribution", kTH1F, {{100, 0., 100., "p_{T,ch} (GeV/c)"}}); - spectra.add("CentralityM", "#it{p}_{T} distribution", kTH1F, {{100, 0., 100., "p_{T,ch} (GeV/c)"}}); - spectra.add("AllTracksPtCent", "#it{p}_{T} distribution", kTH2F, {{200, 0., 200., "p_{T,ch} (GeV/c)"}, centrality}); - - + + if (doprocessMCPartLevel || doprocessMCPartLevelWeighted) { + spectra.add("hRecoilJetRefPtMCP", "#it{p}_{T} distribution of recoil jets", kTH1F, {{200, 0., 200., "p_{T,ch jet} (GeV/c)"}}); + spectra.add("hRecoilJetSigPtMCP", "#it{p}_{T} distribution of recoil jets", kTH1F, {{200, 0., 200., "p_{T,ch jet} (GeV/c)"}}); + spectra.add("hRecoilJetRefCorrPtMCP", "Corrected #it{p}_{T} distribution of recoil jets", kTH1F, {{200, 0., 200., "p_{T,ch jet} (GeV/c)"}}); + spectra.add("hRecoilJetSigCorrPtMCP", "Corrected #it{p}_{T} distribution of recoil jets", kTH1F, {{200, 0., 200., "p_{T,ch jet} (GeV/c)"}}); + + spectra.add("hEventSelectionCountPartLevel", "Pocet", kTH1F, {{10, 0., 10., "Pocet Eventu"}}); + spectra.add("hMCPTrackPt", "#it{p}_{T} distribution of tracks", kTH1F, {pT}); + + spectra.add("hJetPtEtaPhiRhoArea_Part", "Charact. of inclusive part. level jets", kTHnSparseF, {pT, etaJets, phiAngle, rho, jetArea}); + spectra.add("hJetPtMCP", "#it{p}_{T} distribution of jets", kTH1F, {{200, 0., 200., "p_{T,ch jet} (GeV/c)"}}); + + spectra.add("hTTCountPartLevel", "Pocet", kTH1F, {{2, 0., 2., "Pocet TT MCP"}}); } + if (doprocessTrackEfficiency || doprocessTrackEfficiencyMB) { + spectra.add("hTrackingEfficiencyDenominatorPt", "#it{p}_{T} distribution", kTH2F, {{200, 0., 200., "p_{T,ch} (GeV/c)"}, centrality}); + spectra.add("hNonprimaryParticlesPt", "#it{p}_{T} distribution", kTH1F, {{200, 0., 200., "p_{T,ch} (GeV/c)"}}); + spectra.add("hTrackingEfficiencyNumeratorPt", "#it{p}_{T} distribution", kTH2F, {{200, 0., 200., "p_{T,ch} (GeV/c)"}, centrality}); + + spectra.add("hTrackingPhiSmearingPt", "#it{p}_{T} distribution", kTH2F, {pT, dphiAngle}); + spectra.add("hTrackingEtaSmearingPt", "#it{p}_{T} distribution", kTH2F, {pT, detaTracks}); + spectra.add("hTrackingPtResolutionPt", "#it{p}_{T} distribution", kTH2F, {pTTrack, pTRes}); + + spectra.add("hNTracksPerCollision", "Pocet", kTH1F, {{1000, 0., 1000., "Pocet Eventu"}}); + + spectra.add("CentralityC", "#it{p}_{T} distribution", kTH1F, {{100, 0., 100., "p_{T,ch} (GeV/c)"}}); + spectra.add("CentralityM", "#it{p}_{T} distribution", kTH1F, {{100, 0., 100., "p_{T,ch} (GeV/c)"}}); + spectra.add("AllTracksPtCent", "#it{p}_{T} distribution", kTH2F, {{200, 0., 200., "p_{T,ch} (GeV/c)"}, centrality}); + } } //------------------------------------------------------------------------------ - - - - template + + template void fillHistograms(Collision const& collision, Jets const& jets, Tracks const& tracks, float weight = 1.) - + { const float randomConeDistanceFactor = 2.; const int maxRandomConeIterations = 1000; - + std::vector sigTrackPhi; std::vector sigTrackEta; std::vector refTrackPhi; std::vector refTrackEta; - + float multFT0A = collision.multFT0A(); float multFT0C = collision.multFT0C(); float centralityM = collision.centFT0M(); - + spectra.fill(HIST("hMultFT0A"), multFT0A, centralityM, weight); spectra.fill(HIST("hMultFT0C"), multFT0C, centralityM, weight); - - - - float jetR2 = jetR*jetR; + float jetR2 = jetR * jetR; float rho = collision.rho(); - + spectra.fill(HIST("hRho"), rho, weight); - + int i = 0; - for (const auto& track : tracks) - { - // check whether track passes the selection flags - if (skipTrack(track)) - continue; - i++; - auto trackPt = track.pt(); - auto trackPhi = track.phi(); - auto trackEta = track.eta(); - - spectra.fill(HIST("hTrackPt"), trackPt, centralityM, weight); - spectra.fill(HIST("hTrackPhi"), trackPhi, weight); - spectra.fill(HIST("hTrackEta"), trackEta, weight); - - spectra.fill(HIST("hTrackPtPhi"), trackPt, trackPhi, weight); - spectra.fill(HIST("hTrackPtEta"), trackPt, trackEta, weight); - - spectra.fill(HIST("hTrackPtPhiEta"), trackPt, trackPhi, trackEta, weight); - - - - if((signalTriggerMin<=trackPt)&&(trackPt<=signalTriggerMax)){ - sigTrackPhi.push_back(trackPhi); - sigTrackEta.push_back(trackEta); - } - - if((referenceTriggerMin<=trackPt)&&(trackPt<=referenceTriggerMax)){ - refTrackPhi.push_back(trackPhi); - refTrackEta.push_back(trackEta); - } + for (const auto& track : tracks) { + // check whether track passes the selection flags + if (skipTrack(track)) + continue; + i++; + auto trackPt = track.pt(); + auto trackPhi = track.phi(); + auto trackEta = track.eta(); + + spectra.fill(HIST("hTrackPt"), trackPt, centralityM, weight); + spectra.fill(HIST("hTrackPhi"), trackPhi, weight); + spectra.fill(HIST("hTrackEta"), trackEta, weight); + + spectra.fill(HIST("hTrackPtPhi"), trackPt, trackPhi, weight); + spectra.fill(HIST("hTrackPtEta"), trackPt, trackEta, weight); + + spectra.fill(HIST("hTrackPtPhiEta"), trackPt, trackPhi, trackEta, weight); + + if ((signalTriggerMin <= trackPt) && (trackPt <= signalTriggerMax)) { + sigTrackPhi.push_back(trackPhi); + sigTrackEta.push_back(trackEta); + } + + if ((referenceTriggerMin <= trackPt) && (trackPt <= referenceTriggerMax)) { + refTrackPhi.push_back(trackPhi); + refTrackEta.push_back(trackEta); + } } spectra.fill(HIST("hRhovsNTracks"), i, rho, weight); spectra.fill(HIST("hRhovsCentrality"), centralityM, rho, weight); - - - + float perpConePhi = -99; float phiTT = -999; // this will trigger track phi float etaTT = -999; // trigger track eta float rnd = rand->Rndm(); - bool analyzeSignal =0; // 0= reference TT and 1= signal TT - - if (rnd < sigToRefFraction){ - analyzeSignal=1; + bool analyzeSignal = 0; // 0= reference TT and 1= signal TT + + if (rnd < sigToRefFraction) { + analyzeSignal = 1; } - if (analyzeSignal==0&&refTrackPhi.size()>0){ - - int ii = rand->Integer(refTrackPhi.size()); - spectra.fill(HIST("hTTCountData"), 0.5, centralityM, weight); - spectra.fill(HIST("hPocetReferenceTriggeru"), static_cast(refTrackPhi.size()), weight); - phiTT = refTrackPhi[ii]; - etaTT = refTrackEta[ii]; - - }else if(analyzeSignal==1&&sigTrackPhi.size()>0){ - - int ii = rand->Integer(sigTrackPhi.size()); - spectra.fill(HIST("hTTCountData"), 1.5, centralityM, weight); - spectra.fill(HIST("hPocetSignalTriggeru"), static_cast(sigTrackPhi.size()), weight); - phiTT = sigTrackPhi[ii]; - etaTT = sigTrackEta[ii]; + if (analyzeSignal == 0 && refTrackPhi.size() > 0) { + + int ii = rand->Integer(refTrackPhi.size()); + spectra.fill(HIST("hTTCountData"), 0.5, centralityM, weight); + spectra.fill(HIST("hPocetReferenceTriggeru"), static_cast(refTrackPhi.size()), weight); + phiTT = refTrackPhi[ii]; + etaTT = refTrackEta[ii]; + + } else if (analyzeSignal == 1 && sigTrackPhi.size() > 0) { + + int ii = rand->Integer(sigTrackPhi.size()); + spectra.fill(HIST("hTTCountData"), 1.5, centralityM, weight); + spectra.fill(HIST("hPocetSignalTriggeru"), static_cast(sigTrackPhi.size()), weight); + phiTT = sigTrackPhi[ii]; + etaTT = sigTrackEta[ii]; } - - spectra.fill(HIST("hPocetTracku"), i, weight); - - - for (const auto& jet : jets){ - //jet spectra before TT selection - auto jetPt = jet.pt(); - - auto jetPhi = jet.phi(); - auto jetEta = jet.eta(); - float jetArea = jet.area(); - float jetPtCorr = jetPt - rho * jetArea; - - spectra.fill(HIST("hJetPtCorr"), jetPtCorr, weight); - spectra.fill(HIST("hJetPt"), jetPt, weight); - spectra.fill(HIST("hJetPhi"), jetPhi, weight); - spectra.fill(HIST("hJetEta"), jetEta, weight); - - spectra.fill(HIST("hJetPtPhi"), jetPt, jetPhi, weight); - spectra.fill(HIST("hJetPtEta"), jetPt, jetEta, weight); - - spectra.fill(HIST("hJetPtPhiEta"), jetPt, jetPhi, jetEta, weight); - - spectra.fill(HIST("hjetArea"), jetArea, weight); - spectra.fill(HIST("JetAreavsPt"), jetPt, jetArea, weight); - + + spectra.fill(HIST("hPocetTracku"), i, weight); + + for (const auto& jet : jets) { + // jet spectra before TT selection + auto jetPt = jet.pt(); + + auto jetPhi = jet.phi(); + auto jetEta = jet.eta(); + float jetArea = jet.area(); + float jetPtCorr = jetPt - rho * jetArea; + + spectra.fill(HIST("hJetPtCorr"), jetPtCorr, weight); + spectra.fill(HIST("hJetPt"), jetPt, weight); + spectra.fill(HIST("hJetPhi"), jetPhi, weight); + spectra.fill(HIST("hJetEta"), jetEta, weight); + + spectra.fill(HIST("hJetPtPhi"), jetPt, jetPhi, weight); + spectra.fill(HIST("hJetPtEta"), jetPt, jetEta, weight); + + spectra.fill(HIST("hJetPtPhiEta"), jetPt, jetPhi, jetEta, weight); + + spectra.fill(HIST("hjetArea"), jetArea, weight); + spectra.fill(HIST("JetAreavsPt"), jetPt, jetArea, weight); } - - const int minPhiTT = -100; //phi for TT not found - if(phiTT < minPhiTT){ - return; //trigger track was not found + + const int minPhiTT = -100; // phi for TT not found + if (phiTT < minPhiTT) { + return; // trigger track was not found } - - - perpConePhi = TVector2::Phi_mpi_pi(phiTT- constants::math::PIHalf); - - + + perpConePhi = TVector2::Phi_mpi_pi(phiTT - constants::math::PIHalf); + float conePt = 0.; ROOT::Math::PtEtaPhiMVector sumlv(0., 0., 0., 0.); - - for (const auto& track : tracks){ - // check whether track passes the selection flags - if (skipTrack(track)) continue; - - auto trackPt = track.pt(); - auto trackPhi = track.phi(); - auto trackEta = track.eta(); - - float dPhiPerp = TVector2::Phi_mpi_pi(trackPhi - perpConePhi); - float dEtaPerp = trackEta - etaTT; - - if((dPhiPerp)*(dPhiPerp) + (dEtaPerp)*(dEtaPerp) < jetR2){ - ROOT::Math::PtEtaPhiMVector lv(trackPt, trackEta, trackPhi, pionMass); - sumlv += lv; - } + + for (const auto& track : tracks) { + // check whether track passes the selection flags + if (skipTrack(track)) + continue; + + auto trackPt = track.pt(); + auto trackPhi = track.phi(); + auto trackEta = track.eta(); + + float dPhiPerp = TVector2::Phi_mpi_pi(trackPhi - perpConePhi); + float dEtaPerp = trackEta - etaTT; + + if ((dPhiPerp) * (dPhiPerp) + (dEtaPerp) * (dEtaPerp) < jetR2) { + ROOT::Math::PtEtaPhiMVector lv(trackPt, trackEta, trackPhi, pionMass); + sumlv += lv; + } } conePt = static_cast(sumlv.Pt()); - //calculate delta pT from perpendicular cone - float corrConeDeltaPt = conePt - constants::math::PI*jetR2*rho; - if (analyzeSignal==0){ - spectra.fill(HIST("hRefCorrConePt"), corrConeDeltaPt, weight); - } else{ - spectra.fill(HIST("hSigCorrConePt"), corrConeDeltaPt, weight); + // calculate delta pT from perpendicular cone + float corrConeDeltaPt = conePt - constants::math::PI * jetR2 * rho; + if (analyzeSignal == 0) { + spectra.fill(HIST("hRefCorrConePt"), corrConeDeltaPt, weight); + } else { + spectra.fill(HIST("hSigCorrConePt"), corrConeDeltaPt, weight); } - + int jj = 0; const int jet1 = 1; const int jet2 = 2; @@ -472,427 +431,419 @@ struct TrackJetSpectra { float jet1Eta = 0.; float jet2Phi = 0.; float jet2Eta = 0.; - - for (const auto& jet : jets){ - - jj ++; - - auto jetPhi = jet.phi(); - auto jetEta = jet.eta(); - - if(jj == jet1){ - jet1Phi = jetPhi; - jet1Eta = jetEta; - } - - if(jj == jet2){ - jet2Phi = jetPhi; - jet2Eta = jetEta; - } - } - - bool bCloseJet = true; + + for (const auto& jet : jets) { + + jj++; + + auto jetPhi = jet.phi(); + auto jetEta = jet.eta(); + + if (jj == jet1) { + jet1Phi = jetPhi; + jet1Eta = jetEta; + } + + if (jj == jet2) { + jet2Phi = jetPhi; + jet2Eta = jetEta; + } + } + + bool bCloseJet = true; float randJetPhi = -999.; float randJetEta = -999.; int iii = 0; - - while(bCloseJet){ - iii++; - if(iii > maxRandomConeIterations) { - break; - } - - randJetPhi = rand->Uniform(-constants::math::PI, constants::math::PI); - randJetEta = rand->Uniform(-0.5, 0.5); - - float vphi1 = TVector2::Phi_mpi_pi(randJetPhi - jet1Phi); - float vphi2 = TVector2::Phi_mpi_pi(randJetPhi - jet2Phi); - float vphiTT = TVector2::Phi_mpi_pi(randJetPhi - phiTT); - float veta1 = randJetEta - jet1Eta; - float veta2 = randJetEta - jet2Eta; - float vetaTT = randJetEta - etaTT; - float dist1 = std::sqrt(vphi1*vphi1 + veta1*veta1); - float dist2 = std::sqrt(vphi2*vphi2 + veta2*veta2); - float distTT = std::sqrt(vphiTT*vphiTT + vetaTT*vetaTT); - - - if(dist1 > randomConeDistanceFactor * jetR && + + while (bCloseJet) { + iii++; + if (iii > maxRandomConeIterations) { + break; + } + + randJetPhi = rand->Uniform(-constants::math::PI, constants::math::PI); + randJetEta = rand->Uniform(-0.5, 0.5); + + float vphi1 = TVector2::Phi_mpi_pi(randJetPhi - jet1Phi); + float vphi2 = TVector2::Phi_mpi_pi(randJetPhi - jet2Phi); + float vphiTT = TVector2::Phi_mpi_pi(randJetPhi - phiTT); + float veta1 = randJetEta - jet1Eta; + float veta2 = randJetEta - jet2Eta; + float vetaTT = randJetEta - etaTT; + float dist1 = std::sqrt(vphi1 * vphi1 + veta1 * veta1); + float dist2 = std::sqrt(vphi2 * vphi2 + veta2 * veta2); + float distTT = std::sqrt(vphiTT * vphiTT + vetaTT * vetaTT); + + if (dist1 > randomConeDistanceFactor * jetR && dist2 > randomConeDistanceFactor * jetR && - distTT > randomConeDistanceFactor * jetR) bCloseJet = false; + distTT > randomConeDistanceFactor * jetR) + bCloseJet = false; } - + float randConePt = 0.; sumlv = ROOT::Math::PtEtaPhiMVector(0., 0., 0., 0.); - for (const auto& track : tracks){ - if (skipTrack(track)) continue; - - auto trackPt = track.pt(); - auto trackPhi = track.phi(); - auto trackEta = track.eta(); - float dPhiRand = TVector2::Phi_mpi_pi(trackPhi - randJetPhi); - float dEtaRand = trackEta - randJetEta; - - - if((dPhiRand)*(dPhiRand) + (dEtaRand)*(dEtaRand) < jetR2){ - ROOT::Math::PtEtaPhiMVector lv(trackPt, trackEta, trackPhi, pionMass); - sumlv += lv; - } + for (const auto& track : tracks) { + if (skipTrack(track)) + continue; + + auto trackPt = track.pt(); + auto trackPhi = track.phi(); + auto trackEta = track.eta(); + float dPhiRand = TVector2::Phi_mpi_pi(trackPhi - randJetPhi); + float dEtaRand = trackEta - randJetEta; + + if ((dPhiRand) * (dPhiRand) + (dEtaRand) * (dEtaRand) < jetR2) { + ROOT::Math::PtEtaPhiMVector lv(trackPt, trackEta, trackPhi, pionMass); + sumlv += lv; + } } randConePt = static_cast(sumlv.Pt()); - - float corrRandConeDeltaPt = randConePt - constants::math::PI*jetR2*rho; - - if (analyzeSignal==0){ - spectra.fill(HIST("hRefRandCorrConePt"), corrRandConeDeltaPt, centralityM, weight); - }else{ - spectra.fill(HIST("hSigRandCorrConePt"), corrRandConeDeltaPt, centralityM,weight); - } - - - // Recoil jet spectra - for (const auto& jet : jets){ - auto jetPt = jet.pt(); - auto jetPhi = jet.phi(); - float jetArea = jet.area(); - float jetPtCorr = jetPt - rho * jetArea; - - float dPhi = -100; - dPhi = std::fabs(TVector2::Phi_mpi_pi(jetPhi- phiTT)); - - if(dPhi > constants::math::PI - dPhiCut){ - if(analyzeSignal == 0){ - spectra.fill(HIST("hRecoilJetRefPt"), jetPt, centralityM, weight); - spectra.fill(HIST("hRecoilJetRefCorrPt"), jetPtCorr, centralityM, weight); - }else{ - spectra.fill(HIST("hRecoilJetSigPt"), jetPt, centralityM, weight); - spectra.fill(HIST("hRecoilJetSigCorrPt"), jetPtCorr, centralityM, weight); - } - } - } - } - //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - //template + + float corrRandConeDeltaPt = randConePt - constants::math::PI * jetR2 * rho; + + if (analyzeSignal == 0) { + spectra.fill(HIST("hRefRandCorrConePt"), corrRandConeDeltaPt, centralityM, weight); + } else { + spectra.fill(HIST("hSigRandCorrConePt"), corrRandConeDeltaPt, centralityM, weight); + } + + // Recoil jet spectra + for (const auto& jet : jets) { + auto jetPt = jet.pt(); + auto jetPhi = jet.phi(); + float jetArea = jet.area(); + float jetPtCorr = jetPt - rho * jetArea; + + float dPhi = -100; + dPhi = std::fabs(TVector2::Phi_mpi_pi(jetPhi - phiTT)); + + if (dPhi > constants::math::PI - dPhiCut) { + if (analyzeSignal == 0) { + spectra.fill(HIST("hRecoilJetRefPt"), jetPt, centralityM, weight); + spectra.fill(HIST("hRecoilJetRefCorrPt"), jetPtCorr, centralityM, weight); + } else { + spectra.fill(HIST("hRecoilJetSigPt"), jetPt, centralityM, weight); + spectra.fill(HIST("hRecoilJetSigCorrPt"), jetPtCorr, centralityM, weight); + } + } + } + } + //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + // template template void fillMatchedHistograms(JTracksTable const& tracks, JetsBase const& jetsBase, JetsTag const& jetsTag, float weight = 1.) { - - - - std::vector sigTrackPhi; - double phiTT = 0.0; - for (const auto& track : tracks) { - if (skipTrack(track)) continue; - if((signalTriggerMin<=track.pt())&&(track.pt()<=signalTriggerMax)) sigTrackPhi.push_back(track.phi()); - } - - bool bIsThereTTSig = sigTrackPhi.size() > 0; + std::vector sigTrackPhi; + double phiTT = 0.0; - if (bIsThereTTSig) phiTT = sigTrackPhi[rand->Integer(sigTrackPhi.size())]; + for (const auto& track : tracks) { + if (skipTrack(track)) + continue; + if ((signalTriggerMin <= track.pt()) && (track.pt() <= signalTriggerMax)) + sigTrackPhi.push_back(track.phi()); + } + + bool bIsThereTTSig = sigTrackPhi.size() > 0; + + if (bIsThereTTSig) + phiTT = sigTrackPhi[rand->Integer(sigTrackPhi.size())]; for (const auto& jetBase : jetsBase) { - bool bIsBaseJetRecoil = false; /////////////////////////////////////////////////////////////////////////////////////////////////// - - if(bIsThereTTSig){ - float dPhi = std::fabs(TVector2::Phi_mpi_pi(jetBase.phi()- phiTT)); - if(dPhi > constants::math::PI - dPhiCut) bIsBaseJetRecoil = true; - } - - dataForUnfolding(jetBase, jetsTag, bIsBaseJetRecoil, tracks, weight); + bool bIsBaseJetRecoil = false; /////////////////////////////////////////////////////////////////////////////////////////////////// + + if (bIsThereTTSig) { + float dPhi = std::fabs(TVector2::Phi_mpi_pi(jetBase.phi() - phiTT)); + if (dPhi > constants::math::PI - dPhiCut) + bIsBaseJetRecoil = true; + } + + dataForUnfolding(jetBase, jetsTag, bIsBaseJetRecoil, tracks, weight); } } - - //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - + + //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + void processFilteredCollisions(soa::Filtered::iterator const& collision) { // check whether event contains the required flag for event selection (sel8 in this case) - if (skipEvent(collision)) return; + if (skipEvent(collision)) + return; spectra.fill(HIST("hVertexZ_EventFiltering"), collision.posZ()); // compare with "hVertexZ_Cut" histogram } - + PROCESS_SWITCH(TrackJetSpectra, processFilteredCollisions, "process filtered collisions", false); - //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - void processJets(FilteredColl const& collision, soa::Filtered const& jets, soa::Filtered const& tracks) { - if (skipEvent(collision)){ + if (skipEvent(collision)) { return; - } + } // check whether event contains the required flag for event selection (sel8 in this case) - + fillHistograms(collision, jets, tracks); - - } - + } + PROCESS_SWITCH(TrackJetSpectra, processJets, "process inclusive jets", false); - + //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - + void processJetsMCDet(FilteredColl const& collision, - FilteredJetsDetLevel const& jets, - soa::Filtered const& tracks) + FilteredJetsDetLevel const& jets, + soa::Filtered const& tracks) { // check whether event contains the required flag for event selection (sel8 in this case) - if (skipEvent(collision)){ + if (skipEvent(collision)) { return; - } + } fillHistograms(collision, jets, tracks); - - } - + } + PROCESS_SWITCH(TrackJetSpectra, processJetsMCDet, "process MCDet jets", false); - + //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - + void processJetsMCDetWeighted(FilteredCollDetLevelGetWeight const& collision, - FilteredJetsDetLevel const& jets, - soa::Filtered const& tracks) + FilteredJetsDetLevel const& jets, + soa::Filtered const& tracks) { - if (skipEvent(collision)|| collision.isOutlier()){ + if (skipEvent(collision) || collision.isOutlier()) { return; - } + } // check whether event contains the required flag for event selection (sel8 in this case) auto weight = collision.weight(); fillHistograms(collision, jets, tracks, weight); - - } - + } + PROCESS_SWITCH(TrackJetSpectra, processJetsMCDetWeighted, "process MCDet jets", false); - + //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - + void processMCPartLevel(FilteredCollPartLevelMB const& collision, FilteredParticles const& particles, FilteredJetsPartLevel const& jets) { // check whether event contains the required flag for event selection (sel8 in this case) - if (skipEvent(collision)) return; - - fillMCPHistograms(collision, jets, particles); - - - } - + if (skipEvent(collision)) + return; + + fillMCPHistograms(collision, jets, particles); + } + PROCESS_SWITCH(TrackJetSpectra, processMCPartLevel, "process MCDet jets", false); - + //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - + void processMCPartLevelWeighted(FilteredCollPartLevel const& collision, - FilteredParticles const& particles, - FilteredJetsPartLevel const& jets) + FilteredParticles const& particles, + FilteredJetsPartLevel const& jets) { // check whether event contains the required flag for event selection (sel8 in this case) - if (skipEvent(collision)|| collision.isOutlier()) return; - + if (skipEvent(collision) || collision.isOutlier()) + return; + auto weight = collision.weight(); fillMCPHistograms(collision, jets, particles, weight); - - } - + } + PROCESS_SWITCH(TrackJetSpectra, processMCPartLevelWeighted, "process MCDet jets", false); - + //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - + void processTrackEfficiency(CollPartLevelEfficiencyWeighted const& collisionMC, CollisionMCD const& collision, TracksMCDetFullTable const& tracks, ParticlesFullTable const& mcParticles) - { - if (collisionMC.isOutlier()) return; - auto weight = collisionMC.weight(); - - fillEfficiency(collisionMC, - collision, - tracks, - mcParticles, - weight); + { + if (collisionMC.isOutlier()) + return; + auto weight = collisionMC.weight(); + + fillEfficiency(collisionMC, + collision, + tracks, + mcParticles, + weight); } - -PROCESS_SWITCH(TrackJetSpectra, processTrackEfficiency, "process efficiency of tracks", false); - - //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - void processTrackEfficiencyMB(CollPartLevelMBEfficiency const& collisionMC, - CollisionMCD const& collision, - TracksMCDetFullTable const& tracks, - ParticlesFullTable const& mcParticles) -{ - // if (skipEvent(collision)) return; - - - fillEfficiency(collisionMC, - collision, - tracks, - mcParticles); - + + PROCESS_SWITCH(TrackJetSpectra, processTrackEfficiency, "process efficiency of tracks", false); + + //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + void processTrackEfficiencyMB(CollPartLevelMBEfficiency const& collisionMC, + CollisionMCD const& collision, + TracksMCDetFullTable const& tracks, + ParticlesFullTable const& mcParticles) + { + // if (skipEvent(collision)) return; + + fillEfficiency(collisionMC, + collision, + tracks, + mcParticles); } -PROCESS_SWITCH(TrackJetSpectra, processTrackEfficiencyMB, "process efficiency of tracks", false); - - - //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - - template + PROCESS_SWITCH(TrackJetSpectra, processTrackEfficiencyMB, "process efficiency of tracks", false); + + //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + template void fillEfficiency(JCollisionMC const& collisionMC, JCollision const& collisions, JTracks const& tracks, JMcParticles const& mcParticles, float weight = 1.) - { if (!(std::fabs(collisionMC.posZ()) < vertexZCut)) { return; } - + if (collisions.size() < 1) { return; } - + if (!jetderiveddatautilities::selectCollision(collisions.begin(), eventSelectionBits, skipMBGapEvents)) { // Skipping MC events that have their first associated collision not reconstructed return; } - + float centralityC = collisions.begin().centFT0C(); - + float centralityM = collisions.begin().centFT0M(); - - spectra.fill(HIST("CentralityC"), centralityC, weight); - + spectra.fill(HIST("CentralityM"), centralityM, weight); - - for (auto const& particle : mcParticles){ - auto pdgParticle = pdg->GetParticle(particle.pdgCode()); - if (!pdgParticle) continue; - - if (static_cast(pdgParticle->Charge()) == 0) continue; - if (!particle.isPhysicalPrimary()) continue; - - if (std::fabs(particle.eta()) > trkEtaCut) continue; - if (particle.pt() < trkPtMin) continue; - - - - spectra.fill(HIST("hTrackingEfficiencyDenominatorPt"), particle.pt(), centralityM, weight); - } - - for (auto const& collision : collisions){ - if (skipEvent(collision)) continue; - if (!(std::fabs(collision.posZ()) < vertexZCut)) continue; - - - auto collTracks = tracks.sliceBy(tracksPerJCollision, collision.globalIndex()); - spectra.fill(HIST("hNTracksPerCollision"), collTracks.size()); - //for (auto const& particle : particles){ - - - for (auto const& track : collTracks) { - if (skipTrack(track)) continue; - - if (std::fabs(track.eta()) > trkEtaCut) continue; - if (track.pt() < trkPtMin) continue; - - spectra.fill(HIST("AllTracksPtCent"), track.pt(), centralityM, weight); - - - if (!track.has_mcParticle()) continue; - - - auto mcPart = track.template mcParticle_as(); - - auto pdgParticle = pdg->GetParticle(mcPart.pdgCode()); - if (static_cast(pdgParticle->Charge()) == 0) continue; - - if (std::fabs(mcPart.eta()) > trkEtaCut) continue; - if (mcPart.pt() < trkPtMin) continue; - - - if (!mcPart.isPhysicalPrimary()) { - spectra.fill(HIST("hNonprimaryParticlesPt"), track.pt(), weight); - continue; - } - spectra.fill(HIST("hTrackingEfficiencyNumeratorPt"), mcPart.pt(), centralityM, weight); - - float dPhi = TVector2::Phi_mpi_pi(mcPart.phi() - track.phi()); - float dEta = mcPart.eta() - track.eta(); - float dPtResolution = (mcPart.pt() - track.pt())/mcPart.pt(); - - spectra.fill(HIST("hTrackingPhiSmearingPt"), mcPart.pt(), dPhi, weight); - spectra.fill(HIST("hTrackingEtaSmearingPt"),mcPart.pt(), dEta, weight); - spectra.fill(HIST("hTrackingPtResolutionPt"),mcPart.pt(), dPtResolution, weight); - + + for (auto const& particle : mcParticles) { + auto pdgParticle = pdg->GetParticle(particle.pdgCode()); + if (!pdgParticle) + continue; + + if (static_cast(pdgParticle->Charge()) == 0) + continue; + if (!particle.isPhysicalPrimary()) + continue; + + if (std::fabs(particle.eta()) > trkEtaCut) + continue; + if (particle.pt() < trkPtMin) + continue; + + spectra.fill(HIST("hTrackingEfficiencyDenominatorPt"), particle.pt(), centralityM, weight); + } + + for (auto const& collision : collisions) { + if (skipEvent(collision)) + continue; + if (!(std::fabs(collision.posZ()) < vertexZCut)) + continue; + + auto collTracks = tracks.sliceBy(tracksPerJCollision, collision.globalIndex()); + spectra.fill(HIST("hNTracksPerCollision"), collTracks.size()); + // for (auto const& particle : particles){ + + for (auto const& track : collTracks) { + if (skipTrack(track)) + continue; + + if (std::fabs(track.eta()) > trkEtaCut) + continue; + if (track.pt() < trkPtMin) + continue; + + spectra.fill(HIST("AllTracksPtCent"), track.pt(), centralityM, weight); + + if (!track.has_mcParticle()) + continue; + + auto mcPart = track.template mcParticle_as(); + + auto pdgParticle = pdg->GetParticle(mcPart.pdgCode()); + if (static_cast(pdgParticle->Charge()) == 0) + continue; + + if (std::fabs(mcPart.eta()) > trkEtaCut) + continue; + if (mcPart.pt() < trkPtMin) + continue; + + if (!mcPart.isPhysicalPrimary()) { + spectra.fill(HIST("hNonprimaryParticlesPt"), track.pt(), weight); + continue; + } + spectra.fill(HIST("hTrackingEfficiencyNumeratorPt"), mcPart.pt(), centralityM, weight); + + float dPhi = TVector2::Phi_mpi_pi(mcPart.phi() - track.phi()); + float dEta = mcPart.eta() - track.eta(); + float dPtResolution = (mcPart.pt() - track.pt()) / mcPart.pt(); + + spectra.fill(HIST("hTrackingPhiSmearingPt"), mcPart.pt(), dPhi, weight); + spectra.fill(HIST("hTrackingEtaSmearingPt"), mcPart.pt(), dEta, weight); + spectra.fill(HIST("hTrackingPtResolutionPt"), mcPart.pt(), dPtResolution, weight); + } + } } - } -} -//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - - - template + //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + template void fillMCPHistograms(Collision const& collision, Jets const& jets, - Particles const& particles, float weight = 1.) - + Particles const& particles, float weight = 1.) + { - + spectra.fill(HIST("hEventSelectionCountPartLevel"), 0.5, weight); - - std::vector vPhiOfTT; - + + std::vector vPhiOfTT; + float rnd = rand->Rndm(); bool analyzeSignal = 0; float rho = collision.rho(); - - if (rnd < sigToRefFraction) analyzeSignal=1; - - + + if (rnd < sigToRefFraction) + analyzeSignal = 1; + for (const auto& particle : particles) { auto pdgParticle = pdg->GetParticle(particle.pdgCode()); - if (!pdgParticle) continue; + if (!pdgParticle) + continue; // Need charge and physical primary particles bool bParticleNeutral = (static_cast(pdgParticle->Charge()) == 0); - if (bParticleNeutral || !particle.isPhysicalPrimary()) continue; - + if (bParticleNeutral || !particle.isPhysicalPrimary()) + continue; + float particlePt = particle.pt(); float particlePhi = particle.phi(); - + if (analyzeSignal && (particlePt > signalTriggerMin && particlePt < signalTriggerMax)) { - vPhiOfTT.push_back(particlePhi); + vPhiOfTT.push_back(particlePhi); } if (!analyzeSignal && (particlePt > referenceTriggerMin && particlePt < referenceTriggerMax)) { - vPhiOfTT.push_back(particlePhi); + vPhiOfTT.push_back(particlePhi); } } - - - int ii = -1; //index of selected particle level trigger track - if (vPhiOfTT.size() > 0){ - ii = rand->Integer(vPhiOfTT.size()); - if(analyzeSignal == 0){ - spectra.fill(HIST("hTTCountPartLevel"), 0.5, weight); - }else{ - spectra.fill(HIST("hTTCountPartLevel"), 1.5, weight); - } + + int ii = -1; // index of selected particle level trigger track + if (vPhiOfTT.size() > 0) { + ii = rand->Integer(vPhiOfTT.size()); + if (analyzeSignal == 0) { + spectra.fill(HIST("hTTCountPartLevel"), 0.5, weight); + } else { + spectra.fill(HIST("hTTCountPartLevel"), 1.5, weight); + } } - - + for (const auto& jet : jets) { float jetPt = jet.pt(); float jetPhi = jet.phi(); @@ -901,35 +852,35 @@ PROCESS_SWITCH(TrackJetSpectra, processTrackEfficiencyMB, "process efficiency of spectra.fill(HIST("hJetPtEtaPhiRhoArea_Part"), jetPt, jet.eta(), jet.phi(), rho, jetArea, weight); spectra.fill(HIST("hJetPtMCP"), jetPt, weight); - - if (ii > -1) { //only events with TT particle level - - float phiTT = vPhiOfTT[ii]; - - float deltaPhi = std::fabs(TVector2::Phi_mpi_pi(phiTT- jetPhi)); - if (deltaPhi > constants::math::PI - dPhiCut){ - if(analyzeSignal == 0){ - spectra.fill(HIST("hRecoilJetRefPtMCP"), jetPt, weight); - spectra.fill(HIST("hRecoilJetRefCorrPtMCP"), jetPtCorr, weight); - }else{ - spectra.fill(HIST("hRecoilJetSigPtMCP"), jetPt, weight); - spectra.fill(HIST("hRecoilJetSigCorrPtMCP"), jetPtCorr, weight); - } - } + + if (ii > -1) { // only events with TT particle level + + float phiTT = vPhiOfTT[ii]; + + float deltaPhi = std::fabs(TVector2::Phi_mpi_pi(phiTT - jetPhi)); + if (deltaPhi > constants::math::PI - dPhiCut) { + if (analyzeSignal == 0) { + spectra.fill(HIST("hRecoilJetRefPtMCP"), jetPt, weight); + spectra.fill(HIST("hRecoilJetRefCorrPtMCP"), jetPtCorr, weight); + } else { + spectra.fill(HIST("hRecoilJetSigPtMCP"), jetPt, weight); + spectra.fill(HIST("hRecoilJetSigCorrPtMCP"), jetPtCorr, weight); + } + } } } } - - + //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - + void processJetsMatched(FilteredCollDetLevelGetWeight const& collision, aod::JetMcCollisions const&, FilteredTracks const& tracks, FilteredMatchedJetsDetLevel const& mcdjets, FilteredMatchedJetsPartLevel const& mcpjets) { - if (skipEvent(collision)|| collision.isOutlier()) return; + if (skipEvent(collision) || collision.isOutlier()) + return; auto mcpjetsPerMCCollision = mcpjets.sliceBy(partJetsPerCollision, collision.mcCollisionId()); @@ -939,39 +890,38 @@ PROCESS_SWITCH(TrackJetSpectra, processTrackEfficiencyMB, "process efficiency of //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // Auxiliary functions - + template bool skipEvent(const Collision& coll) { /// \brief: trigger cut is needed for pp data return !jetderiveddatautilities::selectCollision(coll, eventSelectionBits, skipMBGapEvents); } - - //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - -template -bool skipMCEvent(const Collision& coll) -{ - return !jetderiveddatautilities::selectCollision(coll, eventSelectionBits); -} - + //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - + + template + bool skipMCEvent(const Collision& coll) + { + return !jetderiveddatautilities::selectCollision(coll, eventSelectionBits); + } + + //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + template bool skipTrack(const Track& track) { - - //check track quality select global tracks - return !jetderiveddatautilities::selectTrack(track, trackSelection); + // check track quality select global tracks + return !jetderiveddatautilities::selectTrack(track, trackSelection); } - + //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - + template bool isJetWithHighPtConstituent(Jet const& jet, Tracks const&) { - //kill all jets which have a constituent track with pT > 100 GeV + // kill all jets which have a constituent track with pT > 100 GeV bool bIsJetWithHighPtConstituent = false; for (const auto& jetConstituent : jet.template tracks_as()) { if (jetConstituent.pt() > maxJetConstituentPt) { @@ -981,45 +931,45 @@ bool skipMCEvent(const Collision& coll) } return bIsJetWithHighPtConstituent; } - + //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - + template void dataForUnfolding(PartJet const& partJet, DetJet const& detJets, bool bIsBaseJetRecoil, TracksTable const& tracks, float weight = 1.) { - //function which matches particle jet with detector level jet + // function which matches particle jet with detector level jet float partJetPt = partJet.pt(); bool bIsThereMatchedJet = partJet.has_matchedJetGeo(); - + if (bIsThereMatchedJet) { - const auto& jetsMatched = partJet.template matchedJetGeo_as>(); - - for (const auto& jetMatched : jetsMatched) { - // skip matches where detector level jets have a constituent with pT above specified cut - bool skipMatchedDetJet = isJetWithHighPtConstituent(jetMatched, tracks); - - if (skipMatchedDetJet) { - // Miss jets - spectra.fill(HIST("hMissedJets_pT"), partJetPt, weight); - - if (bIsBaseJetRecoil) spectra.fill(HIST("hMissedJets_pT_RecoilJets"), partJetPt, weight); - - } else { - float detJetPt = jetMatched.pt(); - - spectra.fill(HIST("hJetPt_resolution"), (partJetPt - detJetPt) / partJetPt, partJetPt, weight); - spectra.fill(HIST("hJetPhi_resolution"), partJet.phi() - jetMatched.phi(), partJetPt, weight); - - if (bIsBaseJetRecoil) { - spectra.fill(HIST("hJetPt_DetLevel_vs_PartLevel_RecoilJets"), detJetPt, partJetPt, weight); - spectra.fill(HIST("hJetPt_resolution_RecoilJets"), (partJetPt - detJetPt) / partJetPt, partJetPt, weight); - spectra.fill(HIST("hJetPhi_resolution_RecoilJets"), partJet.phi() - jetMatched.phi(), partJetPt, weight); - } - } + const auto& jetsMatched = partJet.template matchedJetGeo_as>(); + + for (const auto& jetMatched : jetsMatched) { + // skip matches where detector level jets have a constituent with pT above specified cut + bool skipMatchedDetJet = isJetWithHighPtConstituent(jetMatched, tracks); + + if (skipMatchedDetJet) { + // Miss jets + spectra.fill(HIST("hMissedJets_pT"), partJetPt, weight); + + if (bIsBaseJetRecoil) + spectra.fill(HIST("hMissedJets_pT_RecoilJets"), partJetPt, weight); + + } else { + float detJetPt = jetMatched.pt(); + + spectra.fill(HIST("hJetPt_resolution"), (partJetPt - detJetPt) / partJetPt, partJetPt, weight); + spectra.fill(HIST("hJetPhi_resolution"), partJet.phi() - jetMatched.phi(), partJetPt, weight); + + if (bIsBaseJetRecoil) { + spectra.fill(HIST("hJetPt_DetLevel_vs_PartLevel_RecoilJets"), detJetPt, partJetPt, weight); + spectra.fill(HIST("hJetPt_resolution_RecoilJets"), (partJetPt - detJetPt) / partJetPt, partJetPt, weight); + spectra.fill(HIST("hJetPhi_resolution_RecoilJets"), partJet.phi() - jetMatched.phi(), partJetPt, weight); + } + } } } } - }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)