diff --git a/PWGLF/Tasks/Resonances/kstarInOO.cxx b/PWGLF/Tasks/Resonances/kstarInOO.cxx index 47259ec84ed..d17c9e87110 100644 --- a/PWGLF/Tasks/Resonances/kstarInOO.cxx +++ b/PWGLF/Tasks/Resonances/kstarInOO.cxx @@ -70,7 +70,6 @@ struct kstarInOO { // Event Selection Configurable cfgEventVtxCut{"cfgEventVtxCut", 10.0, "V_z cut selection"}; - 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, 110.0}, "Binning of the centrality axis"}; // Track Selection @@ -81,6 +80,7 @@ struct kstarInOO { Configurable cfgTrackMaxDCAzToPVcut{"cfgTrackMaxDCAzToPVcut", 2.0, "Track DCAz cut to PV Maximum"}; Configurable cfgTrackPrimaryTrack{"cfgTrackPrimaryTrack", true, "Primary track selection"}; // kGoldenChi2 | kDCAxy | kDCAz Configurable cfgTrackConnectedToPV{"cfgTrackConnectedToPV", true, "PV contributor track selection"}; // PV Contriuibutor + Configurable cfgGlobalTrack{"cfgGlobalTrack", false, "Global track selection"}; // kGoldenChi2 | kDCAxy | kDCAz Configurable cfgTrackGlobalWoDCATrack{"cfgTrackGlobalWoDCATrack", true, "Global track selection without DCA"}; // kQualityTracks (kTrackType | kTPCNCls | kTPCCrossedRows | kTPCCrossedRowsOverNCls | kTPCChi2NDF | kTPCRefit | kITSNCls | kITSChi2NDF | kITSRefit | kITSHits) | kInAcceptanceTracks (kPtRange | kEtaRange) // TPC Configurable cfgTracknFindableTPCClusters{"cfgTrackFindableTPCClusters", 50, "nFindable TPC Clusters"}; @@ -88,7 +88,7 @@ struct kstarInOO { Configurable cfgTracknRowsOverFindable{"cfgTrackRowsOverFindable", 1.2, "nRowsOverFindable TPC CLusters"}; Configurable cfgTracknTPCChi2{"cfgTrackTPCChi2", 4.0, "nTPC Chi2 per Cluster"}; - // IT + // ITS Configurable cfgTracknITSChi2{"cfgTrackITSChi2", 36.0, "nITS Chi2 per Cluster"}; // PID @@ -103,6 +103,12 @@ struct kstarInOO { ConfigurableAxis cfgBinsMixVtx{"cfgBinsMixVtx", {VARIABLE_WIDTH, -10.0f, -5.f, 0.f, 5.f, 10.f}, "Mixing bins - z-vertex"}; Configurable cfgMixNMixedEvents{"cfgMixNMixedEvents", 10, "Number of mixed events per event"}; + // MCGen + Configurable cfgForceGenReco{"cfgForceGenReco", false, "Only consider events which are reconstructed (neglect event-loss)"}; + // Configurable cfgForceBR{"cfgForceBR", false, "Only consider K*0->K(pm)pi(mp)"}; + // Configurable cfgForceKaonAcceptence{"cfgForceKaonAcceptence", false, "Only consider K*0's whose daughters decay inside acceptence (no signal loss)"}; + // Configurable cfgForcePionAcceptence{"cfgForcePionAcceptence", false, "Only consider K*0's whose daughters decay inside acceptence (no signal loss)"}; + // Pair Configurable cfgMinvNBins{"cfgMinvNBins", 300, "Number of bins for Minv axis"}; Configurable cfgMinvMin{"cfgMinvMin", 0.60, "Minimum Minv value"}; @@ -111,8 +117,8 @@ struct kstarInOO { // Histogram Configurable cfgEventCutQA{"cfgEventCutsQA", false, "Enable Event QA Hists"}; Configurable cfgTrackCutQA{"cfgTrackCutQA", false, "Enable Track QA Hists"}; - - // std::vector eventSelectionBits; + Configurable cfgDataHistos{"cfgDataHistos", false, "Enable Data Hists"}; + Configurable cfgMcHistos{"cfgMcHistos", false, "Enable MC Hists"}; void init(o2::framework::InitContext&) { @@ -129,9 +135,6 @@ struct kstarInOO { } if (cfgTrackCutQA) { - // histos.add("h_rawpT", "h_rawpT", kTH1F, {{1000, 0.0, 10.0}}); - // histos.add("h_rawpT_Kaon", "h_rawpT_Kaon", kTH1F, {{1000, 0.0, 10.0}}); - // histos.add("h_rawpT_Pion", "h_rawpT_Pion", kTH1F, {{1000, 0.0, 10.0}}); // histos.add("h_eta", "h_eta", kTH1F, {axisEta}); // histos.add("h_phi", "h_phi", kTH1F, {axisPhi}); @@ -152,23 +155,34 @@ struct kstarInOO { histos.add("QA_kaon_TPC_TOF_AC", "QA_kaon_TPC_TOF_AC", {HistType::kTH2F, {pidAxis, pidAxis}}); } - // MC histos - histos.add("hMC_USS", "hMC_USS", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); - histos.add("hMC_LSS", "hMC_LSS", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); - histos.add("hMC_USS_Mix", "hMC_USS_Mix", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); - histos.add("hMC_LSS_Mix", "hMC_LSS_Mix", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + if (cfgDataHistos) { + histos.add("nEvents", "nEvents", kTH1F, {{4, 0.0, 4.0}}); + histos.add("nEvents_Mix", "nEvents_Mix", kTH1F, {{4, 0.0, 4.0}}); - // histos.add("hMC_pt_Pion", "hMC_pt_Pion", kTH1F, {ptAxis}); - // histos.add("hMC_pt_Kaon", "hMC_pt_Kaon", kTH1F, {ptAxis}); - // histos.add("hMC_pt_Proton", "hMC_pt_Proton", kTH1F, {ptAxis}); + histos.add("hUSS", "hUSS", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + histos.add("hLSS", "hLSS", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + histos.add("hUSS_Mix", "hUSS_Mix", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + histos.add("hLSS_Mix", "hLSS_Mix", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + } - // Event Histograms - histos.add("nEvents_MC", "nEvents_MC", kTH1F, {{4, 0.0, 4.0}}); - histos.add("nEvents_MC_Mix", "nEvents_MC_Mix", kTH1F, {{4, 0.0, 4.0}}); + if (cfgMcHistos) { + histos.add("nEvents_MC", "nEvents_MC", kTH1F, {{4, 0.0, 4.0}}); + histos.add("nEvents_MC_Mix", "nEvents_MC_Mix", kTH1F, {{4, 0.0, 4.0}}); + histos.add("nEvents_MC_True", "nEvents_MC_True", kTH1F, {{4, 0.0, 4.0}}); + + histos.add("hMC_USS", "hMC_USS", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + histos.add("hMC_LSS", "hMC_LSS", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + histos.add("hMC_USS_Mix", "hMC_USS_Mix", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + histos.add("hMC_LSS_Mix", "hMC_LSS_Mix", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + histos.add("hMC_USS_True", "hMC_USS_True", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + histos.add("hMC_kstar_True", "hMC_kstar_True", kTHnSparseF, {cfgCentAxis, ptAxis}); + } } // end of init using EventCandidates = soa::Join; //, aod::CentFT0Ms, aod::CentFT0As + using EventCandidatesTrue = aod::McCollisions; + using TrackCandidates = soa::Join; using TrackCandidatesMC = soa::Join; + Partition kaon = !cfgTrackTPCPID || (nabs(aod::pidtpc::tpcNSigmaKa) <= cfgTrackTPCPIDnSig); + Partition pion = !cfgTrackTPCPID || (nabs(aod::pidtpc::tpcNSigmaPi) <= cfgTrackTPCPIDnSig); + Partition kaonMC = !cfgTrackTPCPID || (nabs(aod::pidtpc::tpcNSigmaKa) <= cfgTrackTPCPIDnSig); Partition pionMC = !cfgTrackTPCPID || (nabs(aod::pidtpc::tpcNSigmaPi) <= cfgTrackTPCPIDnSig); @@ -224,6 +241,9 @@ struct kstarInOO { if (std::abs(track.eta()) > cfgTrackMaxEta) return false; + if (cfgGlobalTrack && !track.isGlobalTrack()) + return false; + if (std::abs(track.dcaXY()) > cfgTrackMaxDCArToPVcut) return false; @@ -324,6 +344,43 @@ struct kstarInOO { return false; } + template + void TrackSlicing(const CollisionType& collision1, const TracksType&, const CollisionType& collision2, const TracksType&, const bool IsMix) + { + auto tracks1 = kaon->sliceByCached(aod::track::collisionId, collision1.globalIndex(), cache); + auto tracks2 = pion->sliceByCached(aod::track::collisionId, collision2.globalIndex(), cache); + auto centrality = collision1.centFT0C(); + + for (const auto& [trk1, trk2] : combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { + + if (!trackSelection(trk1) || !trackSelection(trk2)) + continue; + if (!trackPIDKaon(trk1) || !trackPIDPion(trk2)) + continue; + + auto [KstarPt, Minv] = minvReconstruction(trk1, trk2); + if (Minv < 0) + continue; + + double conjugate = trk1.sign() * trk2.sign(); + if (cfgDataHistos) { + if (!IsMix) { + if (conjugate < 0) { + histos.fill(HIST("hUSS"), centrality, KstarPt, Minv); + } else if (conjugate > 0) { + histos.fill(HIST("hLSS"), centrality, KstarPt, Minv); + } + } else { + if (conjugate < 0) { + histos.fill(HIST("hUSS_Mix"), centrality, KstarPt, Minv); + } else if (conjugate > 0) { + histos.fill(HIST("hLSS_Mix"), centrality, KstarPt, Minv); + } + } + } // cfgDataHistos + } // for + } // TrackSlicing + template void TrackSlicingMC(const CollisionType& collision1, const TracksType&, const CollisionType& collision2, const TracksType&, const bool IsMix) { @@ -343,21 +400,66 @@ struct kstarInOO { continue; double conjugate = trk1.sign() * trk2.sign(); - if (!IsMix) { - if (conjugate < 0) { - histos.fill(HIST("hMC_USS"), centrality, KstarPt, Minv); - } else if (conjugate > 0) { - histos.fill(HIST("hMC_LSS"), centrality, KstarPt, Minv); - } - } else { - if (conjugate < 0) { - histos.fill(HIST("hMC_USS_Mix"), centrality, KstarPt, Minv); - } else if (conjugate > 0) { - histos.fill(HIST("hMC_LSS_Mix"), centrality, KstarPt, Minv); + if (cfgMcHistos) { + if (!IsMix) { + if (conjugate < 0) { + histos.fill(HIST("hMC_USS"), centrality, KstarPt, Minv); + } else if (conjugate > 0) { + histos.fill(HIST("hMC_LSS"), centrality, KstarPt, Minv); + } + } else { + if (conjugate < 0) { + histos.fill(HIST("hMC_USS_Mix"), centrality, KstarPt, Minv); + } else if (conjugate > 0) { + histos.fill(HIST("hMC_LSS_Mix"), centrality, KstarPt, Minv); + } } + } // cfgMcHistos + + //====================== + // Gen MC + if (!trk1.has_mcParticle() || !trk2.has_mcParticle()) + continue; + auto particle1 = trk1.mcParticle(); + auto particle2 = trk2.mcParticle(); + if (std::fabs(particle1.pdgCode()) != 321) + continue; // Not Kaon + if (std::fabs(particle2.pdgCode()) != 211) + continue; // Not Pion + + if (!particle1.has_mothers()) { + continue; } - } - } + if (!particle2.has_mothers()) + continue; + + std::vector mothers1{}; + std::vector mothers1PDG{}; + for (auto& particle1_mom : particle1.template mothers_as()) { + mothers1.push_back(particle1_mom.globalIndex()); + mothers1PDG.push_back(particle1_mom.pdgCode()); + } + + std::vector mothers2{}; + std::vector mothers2PDG{}; + for (auto& particle2_mom : particle2.template mothers_as()) { + mothers2.push_back(particle2_mom.globalIndex()); + mothers2PDG.push_back(particle2_mom.pdgCode()); + } + + if (mothers1PDG[0] != 313) + continue; // mother not K*0 + if (mothers2PDG[0] != 313) + continue; // mothers not K*0 + + if (mothers1[0] != mothers2[0]) + continue; // Kaon and pion not from the same K*0 + + if (cfgMcHistos) { + histos.fill(HIST("hMC_USS_True"), centrality, KstarPt, Minv); + } + } // for + } // TrackSlicingMC template std::pair minvReconstruction(const TracksType& trk1, const TracksType& trk2) @@ -384,6 +486,82 @@ struct kstarInOO { return {lResonance.Pt(), lResonance.M()}; } + //======================================================= + //| + //| DATA STUFF (SE) + //| + //======================================================= + + int nEvents = 0; + void processDataSameEvent(EventCandidates::iterator const& collision, TrackCandidates const& tracks) + { + if (cDebugLevel > 0) { + nEvents++; + if ((nEvents + 1) % 10000 == 0) { + std::cout << "Processed Data Events: " << nEvents << std::endl; + } + } + + auto goodEv = eventSelection(collision); + if (cfgDataHistos) { + histos.fill(HIST("nEvents"), 0.5); + } + if (!goodEv) + return; + + bool INELgt0 = false; + for (const auto& track : tracks) { + if (std::fabs(track.eta()) < cfgTrackMaxEta) { + INELgt0 = true; + break; + } + } + if (!INELgt0) + return; + if (cfgDataHistos) { + histos.fill(HIST("nEvents"), 1.5); + } + TrackSlicing(collision, tracks, collision, tracks, false); + + } // processSameEvents + PROCESS_SWITCH(kstarInOO, processDataSameEvent, "process Data Same Event", false); + + //======================================================= + //| + //| DATA STUFF (ME) + //| + //======================================================= + + int nEventsMix = 0; + void processDataMixedEvent(EventCandidates const& collisions, TrackCandidates const& tracks) + { + auto tracksTuple = std::make_tuple(tracks); + BinningType colBinning{{cfgBinsMixVtx, cfgBinsMixMult}, true}; // true is for 'ignore overflows' (true by default) + SameKindPair pairs{colBinning, cfgMixNMixedEvents, -1, collisions, tracksTuple, &cache}; + for (const auto& [collision1, tracks1, collision2, tracks2] : pairs) { + if (cDebugLevel > 0) { + nEventsMix++; + if ((nEventsMix + 1) % 10000 == 0) { + std::cout << "Processed DATA Mixed Events : " << nEventsMix << std::endl; + } + } + auto goodEv1 = eventSelection(collision1); + auto goodEv2 = eventSelection(collision2); + if (cfgDataHistos) { + histos.fill(HIST("nEvents_Mix"), 0.5); + } + + if (!goodEv1 || !goodEv2) + continue; + + if (cfgDataHistos) { + histos.fill(HIST("nEvents_Mix"), 1.5); + } + TrackSlicing(collision1, tracks1, collision2, tracks2, true); + } + } + PROCESS_SWITCH(kstarInOO, processDataMixedEvent, "process DATA Mixed Event", false); + //======================================================= //| //| MC STUFF (SE) @@ -391,19 +569,21 @@ struct kstarInOO { //======================================================= int nEventsMC = 0; - void processSameEventMC(EventCandidates::iterator const& collision, TrackCandidatesMC const& tracks, aod::McParticles const&) + void processMCSameEvent(EventCandidates::iterator const& collision, TrackCandidatesMC const& tracks, aod::McParticles const&) { if (cDebugLevel > 0) { nEventsMC++; if ((nEventsMC + 1) % 10000 == 0) { double histmem = histos.getSize(); std::cout << histmem << std::endl; - std::cout << "process_SameEvent_MC: " << nEventsMC << std::endl; + std::cout << "process MC Same Event : " << nEventsMC << std::endl; } } auto goodEv = eventSelection(collision); - histos.fill(HIST("nEvents_MC"), 0.5); + if (cfgMcHistos) { + histos.fill(HIST("nEvents_MC"), 0.5); + } if (!goodEv) return; @@ -417,11 +597,13 @@ struct kstarInOO { if (!INELgt0) return; - histos.fill(HIST("nEvents_MC"), 1.5); + if (cfgMcHistos) { + histos.fill(HIST("nEvents_MC"), 1.5); + } TrackSlicingMC(collision, tracks, collision, tracks, false); } // processSameEvents_MC - PROCESS_SWITCH(kstarInOO, processSameEventMC, "process Same Event MC", true); + PROCESS_SWITCH(kstarInOO, processMCSameEvent, "process MC Same Event", false); //======================================================= //| @@ -430,7 +612,7 @@ struct kstarInOO { //======================================================= int nEventsMCMix = 0; - void processMixedEventMC(EventCandidates const& collisions, TrackCandidatesMC const& tracks, aod::McParticles const&) + void processMCMixedEvent(EventCandidates const& collisions, TrackCandidatesMC const& tracks, aod::McParticles const&) { auto tracksTuple = std::make_tuple(tracks); BinningType colBinning{{cfgBinsMixVtx, cfgBinsMixMult}, true}; // true is for 'ignore overflows' (true by default) @@ -439,22 +621,77 @@ struct kstarInOO { if (cDebugLevel > 0) { nEventsMCMix++; if ((nEventsMCMix + 1) % 10000 == 0) { - std::cout << "Processed Mixed Events: " << nEventsMCMix << std::endl; + std::cout << "Processed MC Mixed Events : " << nEventsMCMix << std::endl; } } auto goodEv1 = eventSelection(collision1); auto goodEv2 = eventSelection(collision2); - histos.fill(HIST("nEvents_MC_Mix"), 0.5); + if (cfgMcHistos) { + histos.fill(HIST("nEvents_MC_Mix"), 0.5); + } if (!goodEv1 || !goodEv2) continue; - histos.fill(HIST("nEvents_MC_Mix"), 1.5); - + if (cfgMcHistos) { + histos.fill(HIST("nEvents_MC_Mix"), 1.5); + } TrackSlicingMC(collision1, tracks1, collision2, tracks2, true); } // mixing } // processMixedEvent_MC - PROCESS_SWITCH(kstarInOO, processMixedEventMC, "process Mixed Event MC", false); + PROCESS_SWITCH(kstarInOO, processMCMixedEvent, "process MC Mixed Event", false); + + //======================================================= + //| + //| GENERATED MC STUFF (TRUE) + //| + //======================================================= + + int nEventsTrue = 0; + void processMCTrue(EventCandidatesTrue::iterator const& collision, soa::SmallGroups> const& recocolls, aod::McParticles const& particles) + { + if (cDebugLevel > 0) { + ++nEventsTrue; + if ((nEventsTrue & 10000) == 0) { + std::cout << "Processed MC True Events : " << nEventsTrue << std::endl; + } + } + + if (fabs(collision.posZ()) > cfgEventVtxCut) + return; + + if (recocolls.size() <= 0) { // not reconstructed + if (cfgForceGenReco) { + return; + } + } + + double centrality = -1; + for (auto& recocoll : recocolls) { + centrality = recocoll.centFT0C(); + auto goodEv = eventSelection(recocoll); + + if (cfgMcHistos) { + histos.fill(HIST("nEvents_MC_True"), 0.5); + } + if (!goodEv) + continue; + } // for + + for (auto& particle : particles) { + if (particle.pdgCode() != 313) + continue; // Not K*0 + if (std::fabs(particle.eta()) > cfgTrackMaxEta) + continue; + + if (cfgMcHistos) { + histos.fill(HIST("hMC_kstar_True"), centrality, particle.pt()); + } + + } // loop over particles + + } // processMCTrue + PROCESS_SWITCH(kstarInOO, processMCTrue, "process MC True", false); void processEventsDummy(EventCandidates::iterator const&, TrackCandidates const&) {