diff --git a/PWGLF/Tasks/Resonances/kstarInOO.cxx b/PWGLF/Tasks/Resonances/kstarInOO.cxx index b2840dca788..0a3e99db80c 100644 --- a/PWGLF/Tasks/Resonances/kstarInOO.cxx +++ b/PWGLF/Tasks/Resonances/kstarInOO.cxx @@ -75,9 +75,6 @@ struct kstarInOO { Configurable cfgEventVtxCut{"cfgEventVtxCut", 10.0, "V_z cut selection"}; Configurable cfgEventMaxEta{"cfgEventMaxEta", 1.0, "set INEL event eta cut"}; ConfigurableAxis cfgCentAxis{"cfgCentAxis", {VARIABLE_WIDTH, 0.0, 1.0, 5.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0}, "Binning of the centrality axis"}; - Configurable cfgOccupancySel{"cfgOccupancySel", false, "Occupancy selection"}; - Configurable cfgOccupancyMax{"cfgOccupancyMax", 999999., "maximum occupancy of tracks in neighbouring collisions in a given time range"}; - Configurable cfgOccupancyMin{"cfgOccupancyMin", -100., "minimum occupancy of tracks in neighbouring collisions in a given time range"}; // Track Selection // General @@ -106,9 +103,12 @@ struct kstarInOO { Configurable cfgTrackCirclePIDCut{"cfgTrackCirclePIDCut", true, "Enables PID cut shape circle switch"}; Configurable cfgTrackCircleValue{"cfgTrackCircleValue", 2, "Enables TOF TPC PID circle cut value"}; Configurable cfgTrackTOFHard{"cfgTrackTOFHard", false, "Enables TOF Hard"}; + Configurable isMIDCut{"isMIDCut", false, "Enables Mis-IDentified PID cut"}; Configurable cfgTrackTPCPIDnSig{"cfgTrackTPCPIDnSig", 4.0, "nTPC PID sigma"}; Configurable cfgTrackTOFPIDnSig{"cfgTrackTOFPIDnSig", 4.0, "nTOF PID sigma"}; + Configurable cfgMIDnSig{"cfgMIDnSig", 1.5, "nTPC Mis-IDentification Particle sigma"}; + Configurable cDebugLevel{"cDebugLevel", 0, "Resolution of Debug"}; // Mixing @@ -133,9 +133,8 @@ struct kstarInOO { Configurable cfgEventCutQA{"cfgEventCutsQA", false, "Enable Event QA Hists"}; Configurable cfgTrackCutQA{"cfgTrackCutQA", false, "Enable Track QA Hists"}; Configurable cfgJetQAHistos{"cfgJetQAHistos", false, "Enable Jet QA Histos"}; - Configurable cfgWoJetQA{"cfgWoJetQA", false, "Enable Without Jet QA Histos"}; - Configurable cfgMCHistos{"cfgMCHistos", false, "Enable MC Hists"}; + Configurable cfgMCHistos{"cfgMCHistos", false, "Enable MC Histos"}; Configurable cfgMixedHistos{"cfgMixedHistos", false, "Enable Mixed Histos"}; Configurable cfgJetHistos{"cfgJetHistos", false, "Enable Jet Histos"}; Configurable cfgJetMCHistos{"cfgJetMCHistos", false, "Enable Jet MC Histos"}; @@ -160,9 +159,6 @@ struct kstarInOO { Configurable cfgTriggerMasksTest3{"cfgTriggerMasksTest3", "", "possible JE Trigger masks Test3"}; Configurable cfgForceTrueINELgt{"cfgForceTrueINELgt", true, "Check generated event with True INELgt0"}; - Configurable cfgIsKstar{"cfgIsKstar", false, "Swaps Phi for Kstar analysis"}; - Configurable cfgBR{"cfgBR", false, "Forces Gen. Charged BR Only"}; - std::vector eventSelectionBits; std::vector RealTriggerMaskBits; std::vector triggerMaskBitsTest1; @@ -199,7 +195,6 @@ struct kstarInOO { "kNoITSROFrameBorder", "kNoCollInTimeRangeStandard", "kIsGoodITSLayersAll", - Form("Occupancy < %.0f", cfgOccupancyMax.value), "All passed events"}; for (size_t i = 0; i < eventCutLabels.size(); ++i) { hCutFlow->GetXaxis()->SetBinLabel(i + 1, eventCutLabels[i].c_str()); @@ -224,6 +219,9 @@ struct kstarInOO { histos.add("QA_kaon_TPC_TOF_BC", "QA_kaon_TPC_TOF_BC", {HistType::kTH2F, {pidAxis, pidAxis}}); histos.add("QA_track_pT_BC", "QA_track_pT_BC", kTH1F, {{13, 0.0, 13.0}}); + histos.add("Possible_MID_Pion_TPC", "Possible_MID_Pion_TPC", {HistType::kTH2F, {pidAxis, pidAxis}}); + histos.add("Possible_MID_Kaon_TPC", "Possible_MID_Kaon_TPC", {HistType::kTH2F, {pidAxis, pidAxis}}); + histos.add("hDCArToPv_AC", "DCArToPv_AC", kTH1F, {axisDCAxy}); histos.add("hDCAzToPv_AC", "DCAzToPv_AC", kTH1F, {axisDCAz}); histos.add("hIsPrim_AC", "hIsPrim_AC", kTH1F, {{2, -0.5, 1.5}}); @@ -252,8 +250,6 @@ struct kstarInOO { histos.add("nGoodJets", "The number of good jets", kTH1F, {{6, -0.5, 5.5}}); histos.add("nJetsPerEvent", "The number of jet per event", kTH1F, {{6, -0.5, 5.5}}); - histos.add("rawDimpT", "rawDimpT", kTH2F, {{1000, 0.0, 10.0}, {100, -0.5, 0.5}}); - histos.add("jetTrackEta", "Jet Track Eta", kTH1F, {{100, -1.0, 1.0}}); histos.add("jetTrackPhi", "Jet Track Phi", kTH1F, {{80, -1.0, 7.0}}); } @@ -262,14 +258,28 @@ struct kstarInOO { histos.add("mcpjet_eta", "MC particle-level Jet Eta", kTH1F, {{100, -1.0, 1.0}}); histos.add("mcpjet_phi", "MC particle-level Jet Phi", kTH1F, {{80, -1.0, 7.0}}); - histos.add("Gen_particle_All_BR", "Generated particle All pT (GeV/c)", kTH1F, {{2000, 0., 100.}}); - histos.add("Gen_particle_BR", "Gen_particle_BR", kTH1F, {minvAxis}); - histos.add("Gen_particle_pT", "Generated particle pT (GeV/c)", kTH1F, {{2000, 0., 100.}}); + histos.add("nJet_6_8", "MC particle-level nJet pT 6 < jet pT < 8 (GeV/c)", kTH1F, {{2000, 0., 100.}}); + + histos.add("Gen_KstarPt_All_BR", "Generated K*0 All pT (GeV/c)", kTH1F, {{2000, 0., 100.}}); + histos.add("Gen_Kstar_BR", "Generated K*0 BR", kTH1F, {minvAxis}); + histos.add("Gen_KstarPt_BR", "Generated K*0 pT (GeV/c)", kTH1F, {{2000, 0., 100.}}); - histos.add("dR_taggedjet_kaon", "dR between tagged jet and kaon wo pT cut", {HistType::kTH2F, {{60, 0., 1.5}, {100, 0., 20.}}}); - histos.add("dR_taggedjet_pion", "dR between tagged jet and pion wo pT cut", {HistType::kTH2F, {{60, 0., 1.5}, {100, 0., 20.}}}); - histos.add("dR_taggedjet_all", "dR between tagged jet and kpi", {HistType::kTH2F, {{60, 0., 1.5}, {100, 0., 20.}}}); - histos.add("dR_taggedjet_all_6_8", "dR between tagged jet and kpi 6 < jet pT < 8", {HistType::kTH2F, {{60, 0., 1.5}, {100, 0., 20.}}}); + histos.add("dR_taggedjet_kaon", "dR between tagged jet and kaon wo pT cut", {HistType::kTH2F, {{100, 0., 3.0}, {150, 0., 30.}}}); + histos.add("dR_taggedjet_pion", "dR between tagged jet and pion wo pT cut", {HistType::kTH2F, {{100, 0., 3.0}, {150, 0., 30.}}}); + histos.add("dR_taggedjet_all", "dR between tagged jet and kpi", {HistType::kTH2F, {{100, 0., 3.0}, {150, 0., 30.}}}); + histos.add("dR_taggedjet_all_6_8", "dR between tagged jet and kpi 6 < jet pT < 8", {HistType::kTH2F, {{100, 0., 3.0}, {150, 0., 30.}}}); + + histos.add("dR_kstar_pion", "dR between kstar and pion", {HistType::kTH2F, {{100, 0., 3.0}, {150, 0., 30.}}}); + histos.add("dR_kstar_kaon", "dR between kstar and kaon", {HistType::kTH2F, {{100, 0., 3.0}, {150, 0., 30.}}}); + + histos.add("Kstar_pT_INJet_6_8", "K*0 pT Inside 6-8 Jet", kTH1F, {{1000, 0., 50.}}); + histos.add("Kstar_pT_INJet_8_10", "K*0 pT Inside 8-10 Jet", kTH1F, {{1000, 0., 50.}}); + histos.add("Kstar_pT_INJet_10_12", "K*0 pT Inside 10-12 Jet", kTH1F, {{1000, 0., 50.}}); + + histos.add("JetpT_8to", "Jet pT 8 to", kTH1F, {{2000, 0., 100.}}); + histos.add("recoveredJetpT_6_8to8", "Recovered Jet pT 6-8 --> 8", kTH1F, {{2000, 0., 100.}}); + histos.add("Kstar_pT_6_8to8_IN_recoveredJet", "K*0 pT In recovered Jet pT 6-8 --> 8", kTH1F, {{1000, 0., 50.}}); + histos.add("Kstar_pT_6_8to8_10_IN_recoveredJet", "K*0 pT In recovered Jet pT 6-8 --> 8-10", kTH1F, {{1000, 0., 50.}}); histos.add("missed_kpi_INJets_6_8", "missed kpi In Jets with 6 < jetPt < 8", {HistType::kTH2F, {{120, 0.0, 1.2}, {100, 0., 20.}}}); histos.add("missed_kpi_INJets_8_10", "missed kpi In Jets with 8 < jetPt < 10", {HistType::kTH2F, {{120, 0.0, 1.2}, {100, 0., 20.}}}); @@ -279,22 +289,14 @@ struct kstarInOO { histos.add("missed_kpi_INJets_25_infinite", "missed kpi In Jets with 25 < jetPt < infinite", {HistType::kTH2F, {{120, 0.0, 1.2}, {100, 0., 20.}}}); histos.add("missed_kpi_INJets_8_infinite", "missed kpi In Jets with 8 < jetPt < infinite", {HistType::kTH2F, {{120, 0.0, 1.2}, {100, 0., 20.}}}); - histos.add("recoveredJetpT_6_8to8_10", "recovered Jet pT", kTH1F, {{2000, 0., 100.}}); - histos.add("recoveredJetpT_6_8to8_10_kstarSpectra", "Kstar pT within the recovered Jet pT", kTH1F, {{2000, 0., 100.}}); - - histos.add("normalJetpT_8_kstarSpectra", "kstar pT in Jet > 8GeV/c", kTH1F, {{2000, 0., 100.}}); - histos.add("normalJetpT_6_8_kstarSpectra", "6 GeV/c < kstar pT in Jet < 8 GeV/c", kTH1F, {{2000, 0., 100.}}); - histos.add("normalJetpT_8_10_kstarSpectra", "8 GeV/c < kstar pT in Jet < 10 GeV/c", kTH1F, {{2000, 0., 100.}}); - histos.add("normalJetpT_10_12_kstarSpectra", "10 GeV/c < kstar pT in Jet < 12 GeV/c", kTH1F, {{2000, 0., 100.}}); - histos.add("JetMigration", "bin to bin migration", {HistType::kTH2F, {{150, 0.0, 15.0, "True jet pT (GeV/c)"}, {150, 0., 15., "Recovered jet pT (GeV/c)"}}}); } //////////////////////////////////// histos.add("nEvents", "nEvents", kTH1F, {{7, 0.0, 7.0}}); histos.add("hUSS_KPi", "hUSS_KPi", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); - histos.add("hUSS_PiK", "hUSS_PiK", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); histos.add("hLSS_KPi", "hLSS_KPi", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + histos.add("hUSS_PiK", "hUSS_PiK", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); histos.add("hLSS_PiK", "hLSS_PiK", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); if (cfgMixedHistos) { @@ -457,57 +459,50 @@ struct kstarInOO { } // pion pid Selection } // fill QA - enum class objectType { MB, + enum class objectType { Inclusive, Jets, - MBRecParticle, - JetsRecParticle }; + RecoParticle, + JetsRecoParticle }; template - void fillMinv(objectType type, const TrackType& trk1, const TrackType& trk2, const ROOT::Math::PxPyPzMVector& lReso, double centrality, double jetpt, bool IsMix, bool flip) + void fillMinv(objectType type, const TrackType& trk1, const TrackType& trk2, const ROOT::Math::PxPyPzMVector& lReso, double centrality, double jetpt, bool flip) { double conjugate = trk1.sign() * trk2.sign(); switch (type) { - case objectType::MB: - if (IsMix && cfgMixedHistos) { - if (conjugate < 0) { - if (!flip) - histos.fill(HIST("hUSS_KPi_Mix"), centrality, lReso.Pt(), lReso.M()); - else - histos.fill(HIST("hUSS_PiK_Mix"), centrality, lReso.Pt(), lReso.M()); + case objectType::Inclusive: + if (conjugate < 0) { + if (!flip) { + histos.fill(HIST("hUSS_KPi"), centrality, lReso.Pt(), lReso.M()); + } else { + histos.fill(HIST("hUSS_PiK"), centrality, lReso.Pt(), lReso.M()); } - } else { - if (conjugate < 0) { - if (!flip) - histos.fill(HIST("hUSS_KPi"), centrality, lReso.Pt(), lReso.M()); - else - histos.fill(HIST("hUSS_PiK"), centrality, lReso.Pt(), lReso.M()); - } else if (conjugate > 0) { - if (!flip) - histos.fill(HIST("hLSS_KPi"), centrality, lReso.Pt(), lReso.M()); - else - histos.fill(HIST("hLSS_PiK"), centrality, lReso.Pt(), lReso.M()); + } else if (conjugate > 0) { + if (!flip) { + histos.fill(HIST("hLSS_KPi"), centrality, lReso.Pt(), lReso.M()); + } else { + histos.fill(HIST("hLSS_PiK"), centrality, lReso.Pt(), lReso.M()); } } break; case objectType::Jets: - if (!IsMix && cfgJetHistos) { + if (cfgJetHistos) { if (conjugate < 0) { - if (!flip) + if (!flip) { histos.fill(HIST("hUSS_KPi_INSIDE"), centrality, jetpt, lReso.Pt(), lReso.M()); - else + } else { histos.fill(HIST("hUSS_PiK_INSIDE"), centrality, jetpt, lReso.Pt(), lReso.M()); - } else { - if (conjugate > 0) { - if (!flip) - histos.fill(HIST("hLSS_KPi_INSIDE"), centrality, jetpt, lReso.Pt(), lReso.M()); - else - histos.fill(HIST("hLSS_PiK_INSIDE"), centrality, jetpt, lReso.Pt(), lReso.M()); + } + } else if (conjugate > 0) { + if (!flip) { + histos.fill(HIST("hLSS_KPi_INSIDE"), centrality, jetpt, lReso.Pt(), lReso.M()); + } else { + histos.fill(HIST("hLSS_PiK_INSIDE"), centrality, jetpt, lReso.Pt(), lReso.M()); } } } break; - case objectType::MBRecParticle: + case objectType::RecoParticle: if (cfgMCHistos) { if (conjugate < 0) { histos.fill(HIST("hUSS_TrueRec"), centrality, lReso.Pt(), lReso.M()); @@ -515,7 +510,7 @@ struct kstarInOO { } break; - case objectType::JetsRecParticle: + case objectType::JetsRecoParticle: if (cfgJetMCHistos) { if (conjugate < 0) { histos.fill(HIST("hUSS_TrueRec_INSIDE"), centrality, jetpt, lReso.Pt(), lReso.M()); @@ -527,68 +522,71 @@ struct kstarInOO { //====================================================================== template - std::pair eventSelection(const EventType event, const bool QA) + std::pair eventSelection(const EventType& event, const bool QA) { if (cfgEventCutQA && QA) { fillQA(false, event, 1); histos.fill(HIST("hEvent_Cut"), 0); } - if (!event.sel8()) + if (!event.sel8()) { return {false, 1}; + } if (cfgEventCutQA && QA) { histos.fill(HIST("hEvent_Cut"), 1); } - if (std::abs(event.posZ()) > cfgEventVtxCut) + if (std::abs(event.posZ()) > cfgEventVtxCut) { return {false, 2}; + } if (cfgEventCutQA && QA) { histos.fill(HIST("hEvent_Cut"), 2); } - if (!event.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) + if (!event.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) { return {false, 3}; + } if (cfgEventCutQA && QA) { histos.fill(HIST("hEvent_Cut"), 3); } - if (!event.selection_bit(aod::evsel::kNoSameBunchPileup)) + if (!event.selection_bit(aod::evsel::kNoSameBunchPileup)) { return {false, 4}; + } if (cfgEventCutQA && QA) { histos.fill(HIST("hEvent_Cut"), 4); } - if (!event.selection_bit(aod::evsel::kNoTimeFrameBorder)) + if (!event.selection_bit(aod::evsel::kNoTimeFrameBorder)) { return {false, 5}; + } if (cfgEventCutQA && QA) { histos.fill(HIST("hEvent_Cut"), 5); } - if (!event.selection_bit(aod::evsel::kNoITSROFrameBorder)) + if (!event.selection_bit(aod::evsel::kNoITSROFrameBorder)) { return {false, 6}; + } if (cfgEventCutQA && QA) { histos.fill(HIST("hEvent_Cut"), 6); } - if (!event.selection_bit(aod::evsel::kNoCollInTimeRangeStandard)) + if (!event.selection_bit(aod::evsel::kNoCollInTimeRangeStandard)) { return {false, 7}; + } if (cfgEventCutQA && QA) { histos.fill(HIST("hEvent_Cut"), 7); } - if (!event.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) + if (!event.selection_bit(aod::evsel::kIsGoodITSLayersAll)) { return {false, 8}; - if (cfgEventCutQA && QA) { - histos.fill(HIST("hEvent_Cut"), 8); } - if (cfgOccupancySel && (event.trackOccupancyInTimeRange() > cfgOccupancyMax || event.trackOccupancyInTimeRange() < cfgOccupancyMin)) - return {false, 9}; if (cfgEventCutQA && QA) { - histos.fill(HIST("hEvent_Cut"), 9); + histos.fill(HIST("hEvent_Cut"), 8); } if (cfgEventCutQA && QA) { fillQA(true, event, 1); - histos.fill(HIST("hEvent_Cut"), 10); + histos.fill(HIST("hEvent_Cut"), 9); } - return {true, 11}; + return {true, 10}; }; template - std::pair JetEventSelection(const EventType event, const bool QA) + std::pair JetEventSelection(const EventType& event, const bool QA) { if (cfgEventCutQA && QA) { fillQA(false, event, 1); @@ -597,8 +595,9 @@ struct kstarInOO { if (!jetderiveddatautilities::selectCollision(event, eventSelectionBits)) { // sel8 return {false, 1}; } - if (std::abs(event.posZ()) > cfgEventVtxCut) + if (std::abs(event.posZ()) > cfgEventVtxCut) { return {false, 2}; + } if (cfgEventCutQA && QA) { fillQA(true, event, 1); @@ -607,38 +606,51 @@ struct kstarInOO { }; template - bool trackSelection(const TracksType track, const bool QA) + bool trackSelection(const TracksType& track, const bool QA) { if (cfgTrackCutQA && QA) { fillQA(false, track, 2); } - if (cfgTrackGlobalSel && !track.isGlobalTrack()) + if (cfgTrackGlobalSel && !track.isGlobalTrack()) { return false; - if (track.pt() < cfgTrackMinPt) + } + if (track.pt() < cfgTrackMinPt) { return false; - if (std::abs(track.eta()) > cfgTrackMaxEta) + } + if (std::abs(track.eta()) > cfgTrackMaxEta) { return false; - if (std::abs(track.dcaXY()) > cfgTrackMaxDCArToPVcut) + } + if (std::abs(track.dcaXY()) > cfgTrackMaxDCArToPVcut) { return false; - if (std::abs(track.dcaZ()) > cfgTrackMaxDCAzToPVcut) + } + if (std::abs(track.dcaZ()) > cfgTrackMaxDCAzToPVcut) { return false; - if (cfgTrackPrimaryTrack && !track.isPrimaryTrack()) + } + if (cfgTrackPrimaryTrack && !track.isPrimaryTrack()) { return false; - if (cfgTrackGlobalWoDCATrack && !track.isGlobalTrackWoDCA()) + } + if (cfgTrackGlobalWoDCATrack && !track.isGlobalTrackWoDCA()) { return false; - if (cfgTrackFindableTPCClusters > 0 && track.tpcNClsFindable() < cfgTrackFindableTPCClusters) + } + if (cfgTrackFindableTPCClusters > 0 && track.tpcNClsFindable() < cfgTrackFindableTPCClusters) { return false; - if (track.tpcNClsCrossedRows() < cfgTrackTPCCrossedRows) + } + if (track.tpcNClsCrossedRows() < cfgTrackTPCCrossedRows) { return false; - if (cfgTrackRowsOverFindable > 0 && track.tpcCrossedRowsOverFindableCls() > cfgTrackRowsOverFindable) + } + if (cfgTrackRowsOverFindable > 0 && track.tpcCrossedRowsOverFindableCls() > cfgTrackRowsOverFindable) { return false; - if (track.tpcChi2NCl() > cfgTrackTPCChi2) + } + if (track.tpcChi2NCl() > cfgTrackTPCChi2) { return false; - if (track.itsChi2NCl() > cfgTrackITSChi2) + } + if (track.itsChi2NCl() > cfgTrackITSChi2) { return false; - if (cfgTrackConnectedToPV && !track.isPVContributor()) + } + if (cfgTrackConnectedToPV && !track.isPVContributor()) { return false; + } if (cfgTrackCutQA && QA) { fillQA(true, track, 2); @@ -646,6 +658,24 @@ struct kstarInOO { return true; }; + template + bool isMIDPion(const TrackPID& candidate, const bool QA) + { + if (cfgTrackCutQA && QA) { + histos.fill(HIST("Possible_MID_Pion_TPC"), candidate.tpcNSigmaKa(), candidate.tofNSigmaKa()); + } + return std::abs(candidate.tpcNSigmaKa()) < cfgMIDnSig && std::abs(candidate.tpcNSigmaPi()) < cfgMIDnSig; + } + + template + bool isMIDKaon(const TrackPID& candidate, const bool QA) + { + if (cfgTrackCutQA && QA) { + histos.fill(HIST("Possible_MID_Kaon_TPC"), candidate.tpcNSigmaPi(), candidate.tofNSigmaPi()); + } + return std::abs(candidate.tpcNSigmaPi()) < cfgMIDnSig && std::abs(candidate.tpcNSigmaKa()) < cfgMIDnSig; + } + template bool trackPIDKaon(const TrackPID& candidate, const bool QA) { @@ -659,23 +689,27 @@ struct kstarInOO { // TPC if (cfgTrackSquarePIDCut) { - if (std::abs(candidate.tpcNSigmaKa()) < cfgTrackTPCPIDnSig) + if (std::abs(candidate.tpcNSigmaKa()) < cfgTrackTPCPIDnSig) { tpcPIDPassed = true; - if (candidate.hasTOF()) { - if (std::abs(candidate.tofNSigmaKa()) < cfgTrackTOFPIDnSig) { - tofPIDPassed = true; - } - } else { - if (!cfgTrackTOFHard) { - tofPIDPassed = true; + + if (candidate.hasTOF()) { + if (std::abs(candidate.tofNSigmaKa()) < cfgTrackTOFPIDnSig) { + tofPIDPassed = true; + } } else { - tofPIDPassed = false; + if (!cfgTrackTOFHard) { + tofPIDPassed = true; + } else { + tofPIDPassed = false; + } } - } + } // TPCPID cut } // end of square cut + if (cfgTrackCirclePIDCut) { - if (std::abs(candidate.tpcNSigmaKa()) < cfgTrackTPCPIDnSig) + if (std::abs(candidate.tpcNSigmaKa()) < cfgTrackTPCPIDnSig) { tpcpid = std::abs(candidate.tpcNSigmaKa()); + } tofpid = 0; if (candidate.hasTOF()) { @@ -693,9 +727,8 @@ struct kstarInOO { // TPC & TOF if (tpcPIDPassed && tofPIDPassed) { - if (cfgTrackCutQA && QA) { - fillQA(true, candidate, 3); - } + // if (cfgTrackCutQA && QA) + // fillQA(true, candidate, 3); return true; } return false; @@ -714,23 +747,28 @@ struct kstarInOO { // TPC if (cfgTrackSquarePIDCut) { - if (std::abs(candidate.tpcNSigmaPi()) < cfgTrackTPCPIDnSig) + if (std::abs(candidate.tpcNSigmaPi()) < cfgTrackTPCPIDnSig) { tpcPIDPassed = true; - if (candidate.hasTOF()) { - if (std::abs(candidate.tofNSigmaPi()) < cfgTrackTOFPIDnSig) { - tofPIDPassed = true; - } - } else { - if (!cfgTrackTOFHard) { - tofPIDPassed = true; + + if (candidate.hasTOF()) { + if (std::abs(candidate.tofNSigmaPi()) < cfgTrackTOFPIDnSig) { + tofPIDPassed = true; + } } else { - tofPIDPassed = false; + if (!cfgTrackTOFHard) { + tofPIDPassed = true; + } else { + tofPIDPassed = false; + } } - } + } // TPCPID cut + } // end of square cut + if (cfgTrackCirclePIDCut) { - if (std::abs(candidate.tpcNSigmaPi()) < cfgTrackTPCPIDnSig) + if (std::abs(candidate.tpcNSigmaPi()) < cfgTrackTPCPIDnSig) { tpcpid = std::abs(candidate.tpcNSigmaPi()); + } tofpid = 0; if (candidate.hasTOF()) { @@ -748,9 +786,8 @@ struct kstarInOO { // TPC & TOF if (tpcPIDPassed && tofPIDPassed) { - if (cfgTrackCutQA && QA) { - fillQA(true, candidate, 4); - } + // if (cfgTrackCutQA && QA) + // fillQA(true, candidate, 4); return true; } return false; @@ -759,12 +796,27 @@ struct kstarInOO { template ROOT::Math::PxPyPzMVector minvReconstruction(const TracksType& trk1, const TracksType& trk2, const bool QA, const bool flip) { - if (!trackSelection(trk1, false) || !trackSelection(trk2, false)) + if (!trackSelection(trk1, false) || !trackSelection(trk2, false)) { return {}; - if (!trackPIDKaon(trk1, QA) || !trackPIDPion(trk2, QA)) + } + if (!trackPIDKaon(trk1, QA) || !trackPIDPion(trk2, QA)) { return {}; - if (trk1.globalIndex() >= trk2.globalIndex()) + } + if (trk1.globalIndex() >= trk2.globalIndex()) { + return {}; + } + + if (isMIDCut && isMIDPion(trk1, QA)) { return {}; + } + if (isMIDCut && isMIDKaon(trk2, QA)) { + return {}; + } + + if (cfgTrackCutQA && QA) { + fillQA(true, trk1, 3); + fillQA(true, trk2, 4); + } ROOT::Math::PxPyPzMVector lDecayDaughter1, lDecayDaughter2, lResonance; if (!flip) { @@ -776,8 +828,9 @@ struct kstarInOO { } lResonance = lDecayDaughter1 + lDecayDaughter2; - if (std::abs(lResonance.Eta()) > cfgTrackMaxEta) + if (std::abs(lResonance.Eta()) > cfgTrackMaxEta) { return {}; + } return {lResonance}; } // minvReconstruction @@ -786,8 +839,9 @@ struct kstarInOO { ROOT::Math::PxPyPzMVector TrueReconstruction(const TracksType& trk1, const TracksType& trk2) { double conjugate = trk1.sign() * trk2.sign(); - if (conjugate > 0) + if (conjugate > 0) { return {}; + } auto particle1 = trk1.mcParticle(); auto particle2 = trk2.mcParticle(); @@ -813,17 +867,21 @@ struct kstarInOO { mothers2PDG.push_back(particle2_mom.pdgCode()); } - if (mothers1PDG[0] != Kstar0PDG) + if (mothers1PDG[0] != Kstar0PDG) { return {}; // mother not K*0 - if (mothers2PDG[0] != Kstar0PDG) + } + if (mothers2PDG[0] != Kstar0PDG) { return {}; // mothers not K*0 - if (mothers1[0] != mothers2[0]) + } + if (mothers1[0] != mothers2[0]) { return {}; // Kaon and pion not from the same K*0 - - if (std::abs(particle1.pdgCode()) != PionPDG && std::abs(particle1.pdgCode()) != KaonPDG) + } + if (std::abs(particle1.pdgCode()) != PionPDG && std::abs(particle1.pdgCode()) != KaonPDG) { return {}; - if (std::abs(particle2.pdgCode()) != PionPDG && std::abs(particle2.pdgCode()) != KaonPDG) + } + if (std::abs(particle2.pdgCode()) != PionPDG && std::abs(particle2.pdgCode()) != KaonPDG) { return {}; + } double track1_mass, track2_mass; if (std::abs(particle1.pdgCode()) == PionPDG) { @@ -847,16 +905,18 @@ struct kstarInOO { lTrueDaughter2 = ROOT::Math::PxPyPzMVector(particle2.px(), particle2.py(), particle2.pz(), track2_mass); lTrueReso = lTrueDaughter1 + lTrueDaughter2; - if (lTrueReso.M() < 0) + if (lTrueReso.M() < 0) { return {}; - if (std::abs(lTrueReso.Eta()) > cfgTrackMaxEta) + } + if (std::abs(lTrueReso.Eta()) > cfgTrackMaxEta) { return {}; + } return {lTrueReso}; } // TrueReconstruction template - void TrackSlicing(const CollisionType& collision1, const TracksType&, const CollisionType& collision2, const TracksType&, const bool IsMix, const bool QA) + void TrackSlicing(const CollisionType& collision1, const TracksType&, const CollisionType& collision2, const TracksType&, const bool QA) { auto tracks1 = kaon->sliceByCached(aod::track::collisionId, collision1.globalIndex(), cache); auto tracks2 = pion->sliceByCached(aod::track::collisionId, collision2.globalIndex(), cache); @@ -865,46 +925,49 @@ struct kstarInOO { for (const auto& [trk1, trk2] : combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { for (bool flip : {false, true}) { auto lReso = minvReconstruction(trk1, trk2, QA, flip); - if (lReso.M() < 0) + if (lReso.M() < 0) { continue; + } - fillMinv(objectType::MB, trk1, trk2, lReso, centrality, -1.0, IsMix, flip); + fillMinv(objectType::Inclusive, trk1, trk2, lReso, centrality, -1.0, flip); } // flip } // for } // TrackSlicing template - void TrackSlicingMC(const CollisionType& collision1, const TracksType&, const CollisionType& collision2, const TracksType&, const bool IsMix, const bool QA) + void TrackSlicingMC(const CollisionType& collision1, const TracksType&, const CollisionType& collision2, const TracksType&, const bool QA) { auto tracks1 = kaonMC->sliceByCached(aod::track::collisionId, collision1.globalIndex(), cache); auto tracks2 = pionMC->sliceByCached(aod::track::collisionId, collision2.globalIndex(), cache); auto centrality = collision1.centFT0C(); for (const auto& [trk1, trk2] : combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { - if (!trk1.has_mcParticle() || !trk2.has_mcParticle()) + if (!trk1.has_mcParticle() || !trk2.has_mcParticle()) { continue; + } for (bool flip : {false, true}) { auto lReso = minvReconstruction(trk1, trk2, QA, flip); - if (lReso.M() < 0) + if (lReso.M() < 0) { continue; - - fillMinv(objectType::MB, trk1, trk2, lReso, centrality, -1.0, IsMix, flip); + } + fillMinv(objectType::Inclusive, trk1, trk2, lReso, centrality, -1.0, flip); } // flip //============================ //| True Reconstruction //============================ auto lTrueReso = TrueReconstruction(trk1, trk2); - fillMinv(objectType::MBRecParticle, trk1, trk2, lTrueReso, centrality, -1.0, IsMix, false); + fillMinv(objectType::RecoParticle, trk1, trk2, lTrueReso, centrality, -1.0, false); } // tracks } // TrackSlicingMC template - double DistinguishJets(const JetType& jets, ROOT::Math::PxPyPzMVector lResonance) + double DistinguishJets(const JetType& jets, ROOT::Math::PxPyPzMVector& lResonance) { - if (cDebugLevel > 0) + if (cDebugLevel > 0) { LOG(info) << "Found multiple jets to the same phi."; + } double bestR = 0; double bestJetpT = 0; @@ -924,7 +987,7 @@ struct kstarInOO { } template - void JetTrackSlicing(aod::JetCollision const& collision, TracksType const& jetTracks, const JetType& chargedjets, const bool IsMix, const bool QA) + void JetTrackSlicing(aod::JetCollision const& collision, TracksType const& jetTracks, const JetType& chargedjets, const bool QA) { //============= //| Inclusive @@ -936,10 +999,11 @@ struct kstarInOO { for (bool flip : {false, true}) { auto lResonance = minvReconstruction(trk1, trk2, QA, flip); - if (lResonance.M() < 0) + if (lResonance.M() < 0) { continue; + } - fillMinv(objectType::MB, trk1, trk2, lResonance, centrality, -1.0, IsMix, flip); + fillMinv(objectType::Inclusive, trk1, trk2, lResonance, centrality, -1.0, flip); //======== //| Jets @@ -947,13 +1011,13 @@ struct kstarInOO { bool jetFlag = false; int goodjets = 0; double jetpt = 0; - double R = 0; + double dR = 0; for (auto const& jet : chargedjets) { - double phidiff = TVector2::Phi_mpi_pi(jet.phi() - lResonance.Phi()); - double etadiff = jet.eta() - lResonance.Eta(); - R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); + double dphi = TVector2::Phi_mpi_pi(jet.phi() - lResonance.Phi()); + double deta = jet.eta() - lResonance.Eta(); + dR = TMath::Sqrt((deta * deta) + (dphi * dphi)); - if (R < cfgJetR) { + if (dR < cfgJetR) { jetFlag = true; jetpt = jet.pt(); goodjets++; @@ -970,14 +1034,14 @@ struct kstarInOO { } if (jetFlag) { - fillMinv(objectType::Jets, trk1, trk2, lResonance, centrality, jetpt, IsMix, flip); + fillMinv(objectType::Jets, trk1, trk2, lResonance, centrality, jetpt, flip); } // jetFlag } // flip } // jetTracks } // JetTrackSlicing template - void JetTrackSlicingMC(aod::JetCollision const& collision, TracksType const& jetTracks, const JetType& mcdjets, const bool IsMix, const bool QA) + void JetTrackSlicingMC(aod::JetCollision const& collision, TracksType const& jetTracks, const JetType& mcdjets, const bool QA) { //============================ //| MB: Track Reconstruction @@ -987,15 +1051,17 @@ struct kstarInOO { auto trk1 = track1.template track_as(); auto trk2 = track2.template track_as(); - if (!trk1.has_mcParticle() || !trk2.has_mcParticle()) + if (!trk1.has_mcParticle() || !trk2.has_mcParticle()) { continue; + } for (bool flip : {false, true}) { auto lResonance = minvReconstruction(trk1, trk2, QA, flip); - if (lResonance.M() < 0) + if (lResonance.M() < 0) { continue; + } - fillMinv(objectType::MB, trk1, trk2, lResonance, centrality, -1.0, IsMix, flip); + fillMinv(objectType::Inclusive, trk1, trk2, lResonance, centrality, -1.0, flip); //============================== //| Jets: Track Reconstruction //============================== @@ -1022,7 +1088,7 @@ struct kstarInOO { } if (jetFlag) { - fillMinv(objectType::Jets, trk1, trk2, lResonance, centrality, jetpt, IsMix, flip); + fillMinv(objectType::Jets, trk1, trk2, lResonance, centrality, jetpt, flip); } // jetFlag } // filp } // Tracks loop @@ -1041,8 +1107,9 @@ struct kstarInOO { { if (cDebugLevel > 0) { nJetEvents++; - if ((nJetEvents + 1) % 10000 == 0) + if ((nJetEvents + 1) % 10000 == 0) { LOG(info) << "Processed Jet Data Events: " << nJetEvents; + } } histos.fill(HIST("nEvents"), 0.5); // Raw event @@ -1053,13 +1120,16 @@ struct kstarInOO { break; } } - if (!INELgt0) + if (!INELgt0) { return; + } histos.fill(HIST("nEvents"), 1.5); // INEL>0 event auto [goodEv, code] = JetEventSelection(collision, true); - if (!goodEv) + // auto [goodEv, code] = eventSelection(collision, true); + if (!goodEv) { return; + } histos.fill(HIST("nEvents"), 2.5); // After selection // Trigger before we start jet finding @@ -1097,32 +1167,29 @@ struct kstarInOO { } if (!RT) { return; - } else if (RT && (VTtest1 || VTtest2 || VTtest3)) { + } + if (VTtest1 || VTtest2 || VTtest3) { return; } } // Trigger cut histos.fill(HIST("nEvents"), 3.5); // After trigger cut - std::vector jetpT{}; - std::vector jetEta{}; - std::vector jetPhi{}; bool HasJets = false; int nJets = 0; for (auto chargedjet : chargedjets) { - if (std::abs(chargedjet.eta()) > cfgJetMaxEta - cfgJetdR) - return; + if (std::abs(chargedjet.eta()) > (cfgJetMaxEta - cfgJetdR)) { + continue; + } - jetpT.push_back(chargedjet.pt()); - jetEta.push_back(chargedjet.eta()); - jetPhi.push_back(chargedjet.phi()); nJets++; if (cfgJetQAHistos) { histos.fill(HIST("JetpT"), chargedjet.pt()); histos.fill(HIST("JetEta"), chargedjet.eta()); histos.fill(HIST("JetPhi"), chargedjet.phi()); } - if (chargedjet.pt() > cfgJetpT) + if (chargedjet.pt() > cfgJetpT) { HasJets = true; + } } if (cfgJetQAHistos) { histos.fill(HIST("nJetsPerEvent"), nJets); @@ -1132,24 +1199,26 @@ struct kstarInOO { //| Has Jets //==================== if (cfgReqJets) { - if (!HasJets) + if (!HasJets) { return; + } } histos.fill(HIST("nEvents"), 4.5); // Has jets + // Just for checking track QA plot for (auto& jetTrack : jetTracks) { auto originTrack = jetTrack.track_as(); - if (!trackSelection(originTrack, true)) + if (!trackSelection(originTrack, true)) { continue; + } if (cfgJetQAHistos) { - histos.fill(HIST("rawDimpT"), jetTrack.pt(), jetTrack.pt() - originTrack.pt()); histos.fill(HIST("jetTrackEta"), jetTrack.eta()); histos.fill(HIST("jetTrackPhi"), jetTrack.phi()); } } // jetTrack loop - JetTrackSlicing(collision, jetTracks, chargedjets, false, true); + JetTrackSlicing(collision, jetTracks, chargedjets, true); } // ProcessDataJets PROCESS_SWITCH(kstarInOO, processDataJets, "process Data Jets", false); @@ -1163,8 +1232,9 @@ struct kstarInOO { { if (cDebugLevel > 0) { nJetMCEvents++; - if ((nJetMCEvents + 1) % 10000 == 0) + if ((nJetMCEvents + 1) % 10000 == 0) { LOG(info) << "Processed Jet MC Events: " << nJetMCEvents; + } } histos.fill(HIST("nEvents"), 0.5); // Gen event @@ -1175,13 +1245,15 @@ struct kstarInOO { break; } } // jetTrack loop - if (!INELgt0) + if (!INELgt0) { return; + } histos.fill(HIST("nEvents"), 1.5); // INEL>0 event auto [goodEv, code] = JetEventSelection(collision, true); - if (!goodEv) + if (!goodEv) { return; + } histos.fill(HIST("nEvents"), 2.5); // Inclusive event // The trigger option doesn't work @@ -1189,19 +1261,13 @@ struct kstarInOO { // return; // histos.fill(HIST("nEvents"), 3.5); // Events for Inclusive - std::vector mcdjetpT{}; - std::vector mcdjetEta{}; - std::vector mcdjetPhi{}; - bool HasJets = false; int nJets = 0; for (auto mcdjet : mcdjets) { - if (std::abs(mcdjet.eta()) > cfgJetMaxEta - cfgJetdR) - return; + if (std::abs(mcdjet.eta()) > cfgJetMaxEta - cfgJetdR) { + continue; + } - mcdjetpT.push_back(mcdjet.pt()); - mcdjetEta.push_back(mcdjet.eta()); - mcdjetPhi.push_back(mcdjet.phi()); nJets++; if (cfgJetQAHistos) { histos.fill(HIST("JetpT"), mcdjet.pt()); @@ -1211,11 +1277,10 @@ struct kstarInOO { if (mcdjet.pt() > cfgJetpT) { HasJets = true; } - } + } // mcdjets if (cfgJetQAHistos) { histos.fill(HIST("nJetsPerEvent"), nJets); } - //==================== //| Has Jets //==================== @@ -1226,7 +1291,7 @@ struct kstarInOO { } histos.fill(HIST("nEvents"), 3.5); // Jet event - // JetTrackSlicingMC(collision, jetTracks, mcdjets, false, true); + // JetTrackSlicingMC(collision, jetTracks, mcdjets, true); // Instead of using "JetTrackSlicingMC", Once we have to test showing the distribution each condition for (auto& [track1, track2] : combinations(o2::soa::CombinationsUpperIndexPolicy(jetTracks, jetTracks))) { @@ -1243,17 +1308,18 @@ struct kstarInOO { histos.fill(HIST("hEffRecTest0_pT"), lResonanceTest1.Pt()); } - if (!trk1.has_mcParticle() || !trk2.has_mcParticle()) + if (!trk1.has_mcParticle() || !trk2.has_mcParticle()) { continue; + } if (cfgJetMCHistos) { histos.fill(HIST("hEffRecTest1_pT"), lResonanceTest1.Pt()); } - ////////////////////////////// ////////////////////////////// if (cfgReqMcEffTrackQA) { // false, true - if (!trackSelection(trk1, true) || !trackSelection(trk2, false)) + if (!trackSelection(trk1, true) || !trackSelection(trk2, false)) { continue; + } } ////////////////////////////// ////////////////////////////// @@ -1261,15 +1327,14 @@ struct kstarInOO { if (cfgJetMCHistos) { histos.fill(HIST("hEffRecTest2_pT"), lResonanceTest1.Pt()); } - if (cfgReqMcEffPID) { // false, true - if (!trackPIDKaon(trk1, true) || !trackPIDPion(trk2, true)) + if (!trackPIDKaon(trk1, true) || !trackPIDPion(trk2, true)) { continue; + } } if (cfgJetMCHistos) { histos.fill(HIST("hEffRecTest3_pT"), lResonanceTest1.Pt()); } - auto particle1 = trk1.mcParticle(); auto particle2 = trk2.mcParticle(); //////////////////////////////////////////////////////////////////////// @@ -1293,6 +1358,7 @@ struct kstarInOO { if (cfgJetMCHistos) { histos.fill(HIST("hMotherPdg1"), std::abs(mothers1PDG[0])); } + std::vector mothers2{}; std::vector mothers2PDG{}; for (auto& particle2_mom : particle2.template mothers_as()) { @@ -1302,30 +1368,30 @@ struct kstarInOO { if (cfgJetMCHistos) { histos.fill(HIST("hMotherPdg2"), std::abs(mothers2PDG[0])); } - - if (mothers1[0] != mothers2[0]) + if (mothers1[0] != mothers2[0]) { continue; // Kaon and pion not from the same K*0 - + } if (cfgJetMCHistos) { histos.fill(HIST("hEffRecTest5_pT"), lResonanceTest1.Pt()); } - - if (std::abs(particle1.pdgCode()) != KaonPDG) // kaon + // kaon + if (std::abs(particle1.pdgCode()) != KaonPDG) { continue; + } if (cfgJetMCHistos) { histos.fill(HIST("hEffRecTest6_pT"), lResonanceTest1.Pt()); } - - if (std::abs(particle2.pdgCode()) != PionPDG) // pion + // pion + if (std::abs(particle2.pdgCode()) != PionPDG) { continue; - + } if (cfgJetMCHistos) { histos.fill(HIST("hEffRecTest7_pT"), lResonanceTest1.Pt()); } - - if (std::abs(mothers1PDG[0]) != Kstar0PDG) + if (std::abs(mothers1PDG[0]) != Kstar0PDG) { continue; // mother not K*0 + } if (cfgJetMCHistos) { histos.fill(HIST("hEffRecTest8_pT"), lResonanceTest1.Pt()); if (lResonanceTest1.Pt() < 0.4) { @@ -1334,13 +1400,12 @@ struct kstarInOO { } } - if (std::abs(mothers2PDG[0]) != Kstar0PDG) + if (std::abs(mothers2PDG[0]) != Kstar0PDG) { continue; // mothers not K*0 - + } if (cfgJetMCHistos) { histos.fill(HIST("hEffRec_pT"), lResonanceTest1.Pt()); } - } // track loop } // process loop PROCESS_SWITCH(kstarInOO, processMCJets, "process MC Jets", false); @@ -1355,8 +1420,9 @@ struct kstarInOO { { if (cDebugLevel > 0) { nEvents++; - if ((nEvents + 1) % 10000 == 0) + if ((nEvents + 1) % 10000 == 0) { LOG(info) << "Processed Data Events: " << nEvents; + } } histos.fill(HIST("nEvents"), 0.5); @@ -1367,23 +1433,25 @@ struct kstarInOO { break; } } - if (!INELgt0) + if (!INELgt0) { return; + } histos.fill(HIST("nEvents"), 1.5); auto [goodEv, code] = eventSelection(collision, true); - if (!goodEv) + if (!goodEv) { return; + } histos.fill(HIST("nEvents"), 2.5); // for trackQA plot for (const auto& track : tracks) { - if (!trackSelection(track, true)) + if (!trackSelection(track, true)) { continue; + } } - TrackSlicing(collision, tracks, collision, tracks, false, true); - + TrackSlicing(collision, tracks, collision, tracks, true); } // processSameEvents PROCESS_SWITCH(kstarInOO, processDataSameEvent, "process Data Same Event", false); @@ -1401,27 +1469,33 @@ struct kstarInOO { for (const auto& [collision1, tracks1, collision2, tracks2] : pairs) { if (cDebugLevel > 0) { nEventsMix++; - if ((nEventsMix + 1) % 10000 == 0) + if ((nEventsMix + 1) % 10000 == 0) { LOG(info) << "Processed DATA Mixed Events : " << nEventsMix; + } } auto [goodEv1, code1] = eventSelection(collision1, false); auto [goodEv2, code2] = eventSelection(collision2, false); bool VtxMixFlag = false; bool CentMixFlag = false; // bool OccupanacyMixFlag = false; - if (std::abs(collision1.posZ() - collision2.posZ()) <= cfgVtxMixCut) // set default to maybe 10 + if (std::abs(collision1.posZ() - collision2.posZ()) <= cfgVtxMixCut) { // set default to maybe 10 VtxMixFlag = true; - if (std::abs(collision1.centFT0C() - collision2.centFT0C()) <= cfgVtxMixCut) // set default to maybe 10 + } + if (std::abs(collision1.centFT0C() - collision2.centFT0C()) <= cfgVtxMixCut) { // set default to maybe 10 CentMixFlag = true; + } - if (!goodEv1 || !goodEv2) + if (!goodEv1 || !goodEv2) { continue; - if (!CentMixFlag) + } + if (!CentMixFlag) { continue; - if (!VtxMixFlag) + } + if (!VtxMixFlag) { continue; + } - TrackSlicing(collision1, tracks1, collision2, tracks2, true, false); + TrackSlicing(collision1, tracks1, collision2, tracks2, false); } } PROCESS_SWITCH(kstarInOO, processDataMixedEvent, "process DATA Mixed Event", false); @@ -1451,22 +1525,25 @@ struct kstarInOO { break; } } - if (!INELgt0) + if (!INELgt0) { return; + } histos.fill(HIST("nEvents"), 1.5); auto [goodEv, code] = eventSelection(collision, true); // sel8 & vtx - if (!goodEv) + if (!goodEv) { return; + } histos.fill(HIST("nEvents"), 2.5); // for trackQA plot for (const auto& track : tracks) { - if (!trackSelection(track, true)) + if (!trackSelection(track, true)) { continue; + } } - TrackSlicingMC(collision, tracks, collision, tracks, false, true); + TrackSlicingMC(collision, tracks, collision, tracks, true); } // processSameEvents_MC PROCESS_SWITCH(kstarInOO, processSameEventMC, "process Same Event MC", false); @@ -1484,8 +1561,9 @@ struct kstarInOO { for (const auto& [collision1, tracks1, collision2, tracks2] : pairs) { if (cDebugLevel > 0) { nEventsMCMix++; - if ((nEventsMCMix + 1) % 10000 == 0) + if ((nEventsMCMix + 1) % 10000 == 0) { LOG(info) << "Processed Mixed Events: " << nEventsMCMix; + } } auto [goodEv1, code1] = eventSelection(collision1, false); auto [goodEv2, code2] = eventSelection(collision2, false); @@ -1493,7 +1571,7 @@ struct kstarInOO { continue; } - TrackSlicingMC(collision1, tracks1, collision2, tracks2, true, false); + TrackSlicingMC(collision1, tracks1, collision2, tracks2, false); } // mixing } // processMixedEvent_MC PROCESS_SWITCH(kstarInOO, processMixedEventMC, "process Mixed Event MC", false); @@ -1508,22 +1586,24 @@ struct kstarInOO { { if (cDebugLevel > 0) { ++nEventsGen; - if (nEventsGen % 10000 == 0) + if (nEventsGen % 10000 == 0) { LOG(info) << "Processed MC (GEN) Events: " << nEventsGen; + } } if (cfgMCHistos) { histos.fill(HIST("nEvents_Gen"), 0.5); // Gen events } - //======================= //| Event & Signal loss //======================= bool INELgt0 = false; for (auto& particle : mcParticles) { - if (std::abs(particle.eta()) > cfgEventMaxEta) + if (std::abs(particle.eta()) > cfgEventMaxEta) { continue; - if (particle.pt() <= 0.0) + } + if (particle.pt() <= 0.0) { continue; + } INELgt0 = true; break; } @@ -1535,15 +1615,17 @@ struct kstarInOO { if (cfgMCHistos) { histos.fill(HIST("nEvents_Gen"), 1.5); // INEL>0 Gen events & EL: denominator } - - if (std::abs(collision.posZ()) > cfgEventVtxCut) + if (std::abs(collision.posZ()) > cfgEventVtxCut) { return; + } for (auto& particle : mcParticles) { - if (std::abs(particle.pdgCode()) != Kstar0PDG) + if (std::abs(particle.pdgCode()) != Kstar0PDG) { continue; - if (std::abs(particle.eta()) > cfgTrackMaxEta) + } + if (std::abs(particle.eta()) > cfgTrackMaxEta) { continue; + } if (cfgMCHistos) { histos.fill(HIST("hGen_pT_Kstar"), particle.pt()); // SL: denominator @@ -1563,33 +1645,37 @@ struct kstarInOO { bool hasGoodEv = false; for (auto& recocoll : recocolls) { // poorly reconstructed auto [goodEv, code] = eventSelection(recocoll, true); - if (recocoll.posZ() > cfgEventVtxCut) + if (recocoll.posZ() > cfgEventVtxCut) { goodEv = false; - if (!goodEv) + } + if (!goodEv) { continue; + } hasGoodEv = true; centrality = recocoll.centFT0C(); break; } // recocolls - if (!hasGoodEv) + if (!hasGoodEv) { return; + } if (cfgMCHistos) { histos.fill(HIST("nEvents_Gen"), 2.5); // EL: numerator } - for (auto& particle : mcParticles) { - if (std::abs(particle.pdgCode()) != Kstar0PDG) + if (std::abs(particle.pdgCode()) != Kstar0PDG) { continue; // Not K*0 - if (std::abs(particle.eta()) > cfgTrackMaxEta) + } + if (std::abs(particle.eta()) > cfgTrackMaxEta) { continue; - if (particle.pt() < cfgTrackMinPt) + } + if (particle.pt() < cfgTrackMinPt) { continue; + } if (cfgMCHistos) { histos.fill(HIST("hRec_pT_Kstar"), centrality, particle.pt()); // SL: numerator // eff: denominator } - } // loop over particles } // processGen PROCESS_SWITCH(kstarInOO, processGen, "process Generated Particles", false); @@ -1604,8 +1690,9 @@ struct kstarInOO { { if (cDebugLevel > 0) { ++nprocessGenEvents; - if (nprocessGenEvents % 10000 == 0) + if (nprocessGenEvents % 10000 == 0) { LOG(info) << "Processed MC (GEN) Events: " << nprocessGenEvents; + } } if (cfgJetMCHistos) { histos.fill(HIST("nEvents_Gen"), 0.5); @@ -1615,10 +1702,12 @@ struct kstarInOO { //======================= bool INELgt0 = false; for (auto& particle : mcParticles) { - if (std::abs(particle.eta()) > cfgEventMaxEta) + if (std::abs(particle.eta()) > cfgEventMaxEta) { continue; - if (particle.pt() <= 0.0) + } + if (particle.pt() <= 0.0) { continue; + } INELgt0 = true; break; } @@ -1630,15 +1719,17 @@ struct kstarInOO { if (cfgJetMCHistos) { histos.fill(HIST("nEvents_Gen"), 1.5); // EL: denominator } - - if (std::abs(collision.posZ()) > cfgEventVtxCut) + if (std::abs(collision.posZ()) > cfgEventVtxCut) { return; + } for (auto& particle : mcParticles) { - if (std::abs(particle.pdgCode()) != Kstar0PDG) + if (std::abs(particle.pdgCode()) != Kstar0PDG) { continue; - if (std::abs(particle.eta()) > cfgTrackMaxEta) + } + if (std::abs(particle.eta()) > cfgTrackMaxEta) { continue; + } if (cfgJetMCHistos) { histos.fill(HIST("hGen_pT_Kstar"), particle.pt()); // SL: denominator @@ -1656,23 +1747,28 @@ struct kstarInOO { if (goodEv) { goodEv = jetderiveddatautilities::selectTrigger(recocoll, RealTriggerMaskBits); } - if (std::abs(recocoll.posZ()) > cfgEventVtxCut) + if (std::abs(recocoll.posZ()) > cfgEventVtxCut) { goodEv = false; + } - if (!goodEv) + if (!goodEv) { return; + } } // reco.coll if (cfgJetMCHistos) { histos.fill(HIST("nEvents_Gen"), 2.5); // EL: numerator } for (auto& particle : mcParticles) { - if (std::abs(particle.pdgCode()) != Kstar0PDG) + if (std::abs(particle.pdgCode()) != Kstar0PDG) { continue; - if (std::abs(particle.eta()) > cfgTrackMaxEta) + } + if (std::abs(particle.eta()) > cfgTrackMaxEta) { continue; - if (particle.pt() < cfgTrackMinPt) + } + if (particle.pt() < cfgTrackMinPt) { continue; + } /* // Not Yet if (cfg_Force_BR) { @@ -1698,7 +1794,6 @@ struct kstarInOO { if (cfgJetMCHistos) { histos.fill(HIST("hRec_pT_Kstar"), particle.pt()); // SL: numerator and eff: denominator } - } // loop over particles } // end of process PROCESS_SWITCH(kstarInOO, processGenJets, "Process Generated Particles Inclusive&Jets", false); @@ -1708,85 +1803,78 @@ struct kstarInOO { { if (cDebugLevel > 0) { ++ndRtest; - if (ndRtest % 10000 == 0) + if (ndRtest % 10000 == 0) { LOG(info) << "Processed dR test: " << ndRtest; + } } bool INELgt0 = false; for (auto& mcpjet : mcpjets) { - if (std::abs(mcpjet.eta()) > cfgEventMaxEta) + if (std::abs(mcpjet.eta()) > cfgEventMaxEta) { continue; - if (mcpjet.pt() <= 0.0) + } + if (mcpjet.pt() <= 0.0) { continue; + } INELgt0 = true; break; } // Selection event through mcpjet loop - if (!INELgt0) + if (!INELgt0) { return; - - if (std::abs(collision.posZ()) > cfgEventVtxCut) + } + if (std::abs(collision.posZ()) > cfgEventVtxCut) { return; + } histos.fill(HIST("nEvents"), 0.5); for (auto& mcParticle : mcParticles) { - if (std::abs(mcParticle.eta()) > cfgTrackMaxEta) + if (std::abs(mcParticle.eta()) > cfgTrackMaxEta) { continue; - if (!mcParticle.has_daughters()) + } + if (!mcParticle.has_daughters()) { continue; - - ROOT::Math::PxPyPzEVector lResonance; - lResonance = ROOT::Math::PxPyPzEVector(mcParticle.px(), mcParticle.py(), mcParticle.pz(), mcParticle.e()); - - int GenPID = 0; - if (!cfgIsKstar) - GenPID = 333; - else - GenPID = 313; - - if (std::abs(mcParticle.pdgCode()) != GenPID) + } + if (std::abs(mcParticle.pdgCode()) != Kstar0PDG) { continue; + } if (cfgJetdRHistos) { - histos.fill(HIST("Gen_particle_All_BR"), mcParticle.pt()); + histos.fill(HIST("Gen_KstarPt_All_BR"), mcParticle.pt()); } - + // 4vector: K*0 + ROOT::Math::PxPyPzEVector lResonance = ROOT::Math::PxPyPzEVector(mcParticle.px(), mcParticle.py(), mcParticle.pz(), mcParticle.e()); bool skip = false; int daughter_kaon = 0; int daughter_pion = 0; - if (!cfgIsKstar) { - for (auto& daughter : mcParticle.daughters_as()) { - if (std::abs(daughter.pdgCode()) != KaonPDG) - skip = true; - } - } else { - for (auto& daughter : mcParticle.daughters_as()) { - if (std::abs(daughter.pdgCode()) == KaonPDG) - ++daughter_kaon; - else if (std::abs(daughter.pdgCode()) == PionPDG) - ++daughter_pion; + for (auto& daughter : mcParticle.daughters_as()) { + if (std::abs(daughter.pdgCode()) == KaonPDG) { + ++daughter_kaon; + } else if (std::abs(daughter.pdgCode()) == PionPDG) { + ++daughter_pion; } - if (daughter_kaon != 1 || daughter_pion != 1) - skip = true; - } // K*(892) - - if (skip && cfgBR) + } + if (daughter_kaon != 1 || daughter_pion != 1) { + skip = true; + } + if (skip) { continue; + } if (cfgJetdRHistos) { - histos.fill(HIST("Gen_particle_BR"), lResonance.M()); - histos.fill(HIST("Gen_particle_pT"), mcParticle.pt()); + histos.fill(HIST("Gen_Kstar_BR"), lResonance.M()); + histos.fill(HIST("Gen_KstarPt_BR"), mcParticle.pt()); } //================== // Distinguish Jets double bestR = 999; - double bestJetpT = 0; - double bestJetPhi = 0; - double bestJetEta = 0; + ROOT::Math::PxPyPzEVector bestJet; + for (auto& mcpjet : mcpjets) { - if (mcpjet.pt() < cfgJetpT) + if (mcpjet.pt() < cfgJetpT) { continue; + } if (cfgJetdRHistos) { histos.fill(HIST("mcpjet_eta"), mcpjet.eta()); @@ -1796,123 +1884,135 @@ struct kstarInOO { double dphi = TVector2::Phi_mpi_pi(mcpjet.phi() - lResonance.Phi()); double deta = mcpjet.eta() - lResonance.Eta(); - double R = TMath::Sqrt((dphi * dphi) + (deta * deta)); - if (R < bestR) { - bestR = R; - bestJetpT = mcpjet.pt(); - bestJetPhi = mcpjet.phi(); - bestJetEta = mcpjet.eta(); - } + double dR = TMath::Sqrt((dphi * dphi) + (deta * deta)); + if (dR < bestR) { + bestR = dR; + bestJet = ROOT::Math::PxPyPzEVector(mcpjet.px(), mcpjet.py(), mcpjet.pz(), mcpjet.energy()); + } // bestR } // mcpJets - if (bestR > cfgJetdR) + if (bestR > cfgJetdR) { continue; + } + + if (bestJet.Pt() > 6.0 && bestJet.Pt() < 8.0) { + if (cfgJetdRHistos) { + histos.fill(HIST("nJet_6_8"), bestJet.Pt()); + } + } //================== // daughters - double missing_pt = 0; - double dR_kaon, dR_pion; - bool kaon_out = false; - bool pion_out = false; - for (auto& daughter : mcParticle.daughters_as()) { - if (cfgIsKstar) { - if (std::abs(daughter.pdgCode()) == KaonPDG) { + double dR_kaon = 0; + double dR_pion = 0; + double dR_k_kstar = 0; + double dR_pi_kstar = 0; + bool isKaonInJet = false; + bool isPionInJet = false; + ROOT::Math::PxPyPzEVector Daughter, missing_daughter; - double dphi_kaon = TVector2::Phi_mpi_pi(bestJetPhi - daughter.phi()); - double deta_kaon = bestJetEta - daughter.eta(); - dR_kaon = TMath::Sqrt((dphi_kaon * dphi_kaon) + (deta_kaon * deta_kaon)); - - if (bestR < cfgJetdR) { - if (cfgJetdRHistos) { - histos.fill(HIST("dR_taggedjet_kaon"), dR_kaon, lResonance.Pt()); - histos.fill(HIST("dR_taggedjet_all"), dR_kaon, lResonance.Pt()); - } - if (dR_kaon > cfgJetdR) { - kaon_out = true; - missing_pt += daughter.pt(); - } - } // INSIDE Jets - } // kaon daughter - if (std::abs(daughter.pdgCode()) == PionPDG) { + for (auto& daughter : mcParticle.daughters_as()) { + Daughter = ROOT::Math::PxPyPzEVector(daughter.px(), daughter.py(), daughter.pz(), daughter.energy()); + + if (std::abs(daughter.pdgCode()) == KaonPDG) { + //=============== + // K*0 and kaon + double dphi_k_kstar = TVector2::Phi_mpi_pi(lResonance.Phi() - Daughter.Phi()); + double deta_k_kstar = lResonance.Eta() - Daughter.Eta(); + dR_k_kstar = TMath::Sqrt((dphi_k_kstar * dphi_k_kstar) + (deta_k_kstar * deta_k_kstar)); + + // Jet and kaon + double dphi_kaon = TVector2::Phi_mpi_pi(bestJet.Phi() - Daughter.Phi()); + double deta_kaon = bestJet.Eta() - Daughter.Eta(); + dR_kaon = TMath::Sqrt((dphi_kaon * dphi_kaon) + (deta_kaon * deta_kaon)); + + if (dR_kaon <= cfgJetdR) { + isKaonInJet = true; + } else { + missing_daughter += Daughter; + } + } else if (std::abs(daughter.pdgCode()) == PionPDG) { + //=============== + // K*0 and pion + double dphi_pi_kstar = TVector2::Phi_mpi_pi(lResonance.Phi() - Daughter.Phi()); + double deta_pi_kstar = lResonance.Eta() - Daughter.Eta(); + dR_pi_kstar = TMath::Sqrt((dphi_pi_kstar * dphi_pi_kstar) + (deta_pi_kstar * deta_pi_kstar)); + + // Jet and pion + double dphi_pion = TVector2::Phi_mpi_pi(bestJet.Phi() - Daughter.Phi()); + double deta_pion = bestJet.Eta() - Daughter.Eta(); + dR_pion = TMath::Sqrt((dphi_pion * dphi_pion) + (deta_pion * deta_pion)); + + if (dR_pion <= cfgJetdR) { + isPionInJet = true; + } else { + missing_daughter += Daughter; + } + } + } // daughter - double dphi_pion = TVector2::Phi_mpi_pi(bestJetPhi - daughter.phi()); - double deta_pion = bestJetEta - daughter.eta(); - dR_pion = TMath::Sqrt((dphi_pion * dphi_pion) + (deta_pion * deta_pion)); + if (cfgJetdRHistos) { + histos.fill(HIST("dR_kstar_kaon"), dR_k_kstar, lResonance.Pt()); + histos.fill(HIST("dR_kstar_pion"), dR_pi_kstar, lResonance.Pt()); + + if (bestR < cfgJetdR) { + histos.fill(HIST("dR_taggedjet_all"), dR_kaon, lResonance.Pt()); + histos.fill(HIST("dR_taggedjet_all"), dR_pion, lResonance.Pt()); + histos.fill(HIST("dR_taggedjet_kaon"), dR_kaon, lResonance.Pt()); + histos.fill(HIST("dR_taggedjet_pion"), dR_pion, lResonance.Pt()); + + if (bestJet.Pt() > 6.0 && bestJet.Pt() < 8.0) { + histos.fill(HIST("dR_taggedjet_all_6_8"), dR_kaon, lResonance.Pt()); + histos.fill(HIST("dR_taggedjet_all_6_8"), dR_pion, lResonance.Pt()); + histos.fill(HIST("Kstar_pT_INJet_6_8"), lResonance.Pt()); + } + if (bestJet.Pt() > 8.0) { + histos.fill(HIST("JetpT_8to"), bestJet.Pt()); + } + if (bestJet.Pt() > 8.0 && bestJet.Pt() < 10.0) { + histos.fill(HIST("Kstar_pT_INJet_8_10"), lResonance.Pt()); + } + if (bestJet.Pt() > 10.0 && bestJet.Pt() < 12.0) { + histos.fill(HIST("Kstar_pT_INJet_10_12"), lResonance.Pt()); + } - if (bestR < cfgJetdR) { - if (cfgJetdRHistos) { - histos.fill(HIST("dR_taggedjet_pion"), dR_pion, lResonance.Pt()); - histos.fill(HIST("dR_taggedjet_all"), dR_pion, lResonance.Pt()); + // Daughters out of jets + if (!isKaonInJet || !isPionInJet) { + ROOT::Math::PxPyPzEVector recoveredJet = bestJet + missing_daughter; + if (bestJet.Pt() > 6.0 && bestJet.Pt() < 8.0) { + histos.fill(HIST("missed_kpi_INJets_6_8"), (bestJet.Pt() - missing_daughter.Pt()) / bestJet.Pt(), lResonance.Pt()); - if (bestJetpT > 6.0 && bestJetpT < 8.0) - histos.fill(HIST("dR_taggedjet_all_6_8"), dR_pion, lResonance.Pt()); - } - if (dR_pion > cfgJetdR) { - pion_out = true; - missing_pt += daughter.pt(); + if (recoveredJet.Pt() > 8.0) { + histos.fill(HIST("recoveredJetpT_6_8to8"), recoveredJet.Pt()); + histos.fill(HIST("Kstar_pT_6_8to8_IN_recoveredJet"), lResonance.Pt()); } - } // INSIDE Jets - } // pion daughter - } else { - if (std::abs(daughter.pdgCode()) == KaonPDG) { - double dphi_kaon = TVector2::Phi_mpi_pi(bestJetPhi - daughter.phi()); - double deta_kaon = bestJetEta - daughter.eta(); - dR_kaon = TMath::Sqrt((dphi_kaon * dphi_kaon) + (deta_kaon * deta_kaon)); - - if (bestR < cfgJetdR) { - if (cfgJetdRHistos) { - histos.fill(HIST("dR_taggedjet_kaon"), dR_kaon, lResonance.Pt()); - histos.fill(HIST("dR_taggedjet_all"), dR_kaon, lResonance.Pt()); - - if (bestJetpT > 6.0 && bestJetpT < 8.0) - histos.fill(HIST("dR_taggedjet_all_6_8"), dR_kaon, lResonance.Pt()); - } - - if (dR_kaon > cfgJetdR) { - kaon_out = true; - missing_pt = daughter.pt(); + if (recoveredJet.Pt() > 8.0 && recoveredJet.Pt() < 10.0) { + histos.fill(HIST("Kstar_pT_6_8to8_10_IN_recoveredJet"), lResonance.Pt()); } + } // 6 < bestJetpT < 8 GeV/c + if (bestJet.Pt() > 8.0 && bestJet.Pt() < 10.0) { + histos.fill(HIST("missed_kpi_INJets_8_10"), (bestJet.Pt() - missing_daughter.Pt()) / bestJet.Pt(), lResonance.Pt()); } - } // kaon daughter - } // phi(1020) - } // daughter - - if (kaon_out || pion_out) { - double recoveredJetpT = bestJetpT + missing_pt; - if (cfgJetdRHistos) { - if (bestJetpT > 6.0 && bestJetpT < 8.0) { - histos.fill(HIST("normalJetpT_6_8_kstarSpectra"), lResonance.Pt()); - histos.fill(HIST("missed_kpi_INJets_6_8"), (bestJetpT - missing_pt) / bestJetpT, lResonance.Pt()); - if (recoveredJetpT > 8.0) { - histos.fill(HIST("recoveredJetpT_6_8to8_10"), recoveredJetpT); - histos.fill(HIST("recoveredJetpT_6_8to8_10_kstarSpectra"), lResonance.Pt()); + if (bestJet.Pt() > 10.0 && bestJet.Pt() < 12.0) { + histos.fill(HIST("missed_kpi_INJets_10_12"), (bestJet.Pt() - missing_daughter.Pt()) / bestJet.Pt(), lResonance.Pt()); } - } - - if (bestJetpT > 8.0 && bestJetpT < 10.0) { - histos.fill(HIST("missed_kpi_INJets_8_10"), (bestJetpT - missing_pt) / bestJetpT, lResonance.Pt()); - histos.fill(HIST("normalJetpT_8_10_kstarSpectra"), lResonance.Pt()); - } - if (bestJetpT > 10.0 && bestJetpT < 12.0) { - histos.fill(HIST("missed_kpi_INJets_10_12"), (bestJetpT - missing_pt) / bestJetpT, lResonance.Pt()); - histos.fill(HIST("normalJetpT_10_12_kstarSpectra"), lResonance.Pt()); - } - if (bestJetpT > 12.0 && bestJetpT < 15.0) - histos.fill(HIST("missed_kpi_INJets_12_15"), (bestJetpT - missing_pt) / bestJetpT, lResonance.Pt()); - if (bestJetpT > 15.0 && bestJetpT < 25.0) - histos.fill(HIST("missed_kpi_INJets_15_25"), (bestJetpT - missing_pt) / bestJetpT, lResonance.Pt()); - if (bestJetpT > 25.0) - histos.fill(HIST("missed_kpi_INJets_25_infinite"), (bestJetpT - missing_pt) / bestJetpT, lResonance.Pt()); - - if (bestJetpT > 8.0) { - histos.fill(HIST("missed_kpi_INJets_8_infinite"), (bestJetpT - missing_pt) / bestJetpT, lResonance.Pt()); - histos.fill(HIST("normalJetpT_8_kstarSpectra"), lResonance.Pt()); - } - histos.fill(HIST("JetMigration"), bestJetpT, recoveredJetpT); - } // cfgJetdRHistos - } // kaon_out || pion_out + if (bestJet.Pt() > 12.0 && bestJet.Pt() < 15.0) { + histos.fill(HIST("missed_kpi_INJets_12_15"), (bestJet.Pt() - missing_daughter.Pt()) / bestJet.Pt(), lResonance.Pt()); + } + if (bestJet.Pt() > 15.0 && bestJet.Pt() < 25.0) { + histos.fill(HIST("missed_kpi_INJets_15_25"), (bestJet.Pt() - missing_daughter.Pt()) / bestJet.Pt(), lResonance.Pt()); + } + if (bestJet.Pt() > 25.0) { + histos.fill(HIST("missed_kpi_INJets_25_infinite"), (bestJet.Pt() - missing_daughter.Pt()) / bestJet.Pt(), lResonance.Pt()); + } + if (bestJet.Pt() > 8.0) { + histos.fill(HIST("missed_kpi_INJets_8_infinite"), (bestJet.Pt() - missing_daughter.Pt()) / bestJet.Pt(), lResonance.Pt()); + } + histos.fill(HIST("JetMigration"), bestJet.Pt(), recoveredJet.Pt()); + } // Daughters out of jet + } // K*0 InJet + } // cfgJetdRHistos } // mcParticles }; - PROCESS_SWITCH(kstarInOO, processJetQA, "Process dR of K*0 Inclusive and Inside jet", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)