From 39e26b86d1f29cd3498afb7301f4b646772c6d13 Mon Sep 17 00:00:00 2001 From: Adrian Nassirpour Date: Fri, 26 Apr 2024 20:10:09 +0900 Subject: [PATCH 1/2] Full rework of the MC, now fully functional --- PWGJE/Tasks/phiInJets.cxx | 647 +++++++++++++++++++++++++------------- 1 file changed, 436 insertions(+), 211 deletions(-) diff --git a/PWGJE/Tasks/phiInJets.cxx b/PWGJE/Tasks/phiInJets.cxx index f13033fe533..98e828432c9 100644 --- a/PWGJE/Tasks/phiInJets.cxx +++ b/PWGJE/Tasks/phiInJets.cxx @@ -102,12 +102,15 @@ struct phiInJets { JEhistos.add("ptJEHistogramKaon", "ptJEHistogramKaon", kTH1F, {PtAxis}); JEhistos.add("ptJEHistogramProton", "ptJEHistogramProton", kTH1F, {PtAxis}); JEhistos.add("ptJEHistogramPhi", "ptJEHistogramPhi", kTH1F, {PtAxis}); + JEhistos.add("ptJEHistogramPhi_JetTrigger", "ptJEHistogramPhi_JetTrigger", kTH1F, {PtAxis}); JEhistos.add("minvJEHistogramPhi", "minvJEHistogramPhi", kTH1F, {MinvAxis}); JEhistos.add("ptGeneratedPion", "ptGeneratedPion", kTH1F, {PtAxis}); JEhistos.add("ptGeneratedKaon", "ptGeneratedKaon", kTH1F, {PtAxis}); JEhistos.add("ptGeneratedProton", "ptGeneratedProton", kTH1F, {PtAxis}); JEhistos.add("ptGeneratedPhi", "ptGeneratedPhi", kTH1F, {PtAxis}); + JEhistos.add("ptGeneratedPhi_JetTrigger", "ptGeneratedPhi_JetTrigger", kTH1F, {PtAxis}); + JEhistos.add("mGeneratedPhi", "mGeneratedPhi", kTH1F, {MinvAxis}); JEhistos.add("etaHistogram", "etaHistogram", kTH1F, {axisEta}); @@ -126,6 +129,8 @@ struct phiInJets { JEhistos.add("nEvents", "nEvents", kTH1F, {{4, 0.0, 4.0}}); JEhistos.add("nEvents_MCRec", "nEvents_MCRec", kTH1F, {{4, 0.0, 4.0}}); JEhistos.add("nEvents_MCGen", "nEvents_MCGen", kTH1F, {{4, 0.0, 4.0}}); + JEhistos.add("nEvents_MCRec_MATCHED", "nEvents_MCRec_MATCHED", kTH1F, {{4, 0.0, 4.0}}); + JEhistos.add("nEvents_MCGen_MATCHED", "nEvents_MCGen_MATCHED", kTH1F, {{4, 0.0, 4.0}}); JEhistos.add("nJetsPerEvent", "nJetsPerEvent", kTH1F, {{10, 0.0, 10.0}}); @@ -167,6 +172,20 @@ struct phiInJets { JEhistos.add("hMCTrue_hUSS_OUTSIDE", "hMCTrue_hUSS_OUTSIDE", kTH3F, {dRAxis, PtAxis, MinvAxis}); JEhistos.add("hMCTrue_hUSS_OUTSIDE_1D", "hMCTrue_hUSS_OUTSIDE_1D", kTH1F, {MinvAxis}); JEhistos.add("hMCTrue_hUSS_OUTSIDE_1D_2_3", "hMCTrue_hUSS_OUTSIDE_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hMCTrue_hUSS_OUTSIDE_TRIG", "hMCTrue_hUSS_OUTSIDE_TRIG", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCTrue_hUSS_OUTSIDE_TRIG_1D", "hMCTrue_hUSS_OUTSIDE_TRIG_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCTrue_hUSS_OUTSIDE_TRIG_1D_2_3", "hMCTrue_hUSS_OUTSIDE_TRIG_1D_2_3", kTH1F, {MinvAxis}); + + JEhistos.add("hMCTrue_nonmatch_hUSS_INSIDE", "hMCTrue_nonmatch_hUSS_INSIDE", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCTrue_nonmatch_hUSS_INSIDE_1D", "hMCTrue_nonmatch_hUSS_INSIDE_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCTrue_nonmatch_hUSS_INSIDE_1D_2_3", "hMCTrue_nonmatch_hUSS_INSIDE_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hMCTrue_nonmatch_hUSS_OUTSIDE", "hMCTrue_nonmatch_hUSS_OUTSIDE", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCTrue_nonmatch_hUSS_OUTSIDE_1D", "hMCTrue_nonmatch_hUSS_OUTSIDE_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCTrue_nonmatch_hUSS_OUTSIDE_1D_2_3", "hMCTrue_nonmatch_hUSS_OUTSIDE_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG", "hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D", "hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3", "hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_hUSS", "hMCRec_hUSS", kTH3F, {dRAxis, PtAxis, MinvAxis}); JEhistos.add("hMCRec_hUSS_1D", "hMCRec_hUSS_1D", kTH1F, {MinvAxis}); @@ -177,6 +196,24 @@ struct phiInJets { JEhistos.add("hMCRec_hUSS_OUTSIDE", "hMCRec_hUSS_OUTSIDE", kTH3F, {dRAxis, PtAxis, MinvAxis}); JEhistos.add("hMCRec_hUSS_OUTSIDE_1D", "hMCRec_hUSS_OUTSIDE_1D", kTH1F, {MinvAxis}); JEhistos.add("hMCRec_hUSS_OUTSIDE_1D_2_3", "hMCRec_hUSS_OUTSIDE_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_hUSS_OUTSIDE_TRIG", "hMCRec_hUSS_OUTSIDE_TRIG", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCRec_hUSS_OUTSIDE_TRIG_1D", "hMCRec_hUSS_OUTSIDE_TRIG_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_hUSS_OUTSIDE_TRIG_1D_2_3", "hMCRec_hUSS_OUTSIDE_TRIG_1D_2_3", kTH1F, {MinvAxis}); + + JEhistos.add("hMCRec_nonmatch_hUSS", "hMCRec_nonmatch_hUSS", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_1D", "hMCRec_nonmatch_hUSS_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_1D_2_3", "hMCRec_nonmatch_hUSS_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_INSIDE", "hMCRec_nonmatch_hUSS_INSIDE", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_INSIDE_1D", "hMCRec_nonmatch_hUSS_INSIDE_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_INSIDE_1D_2_3", "hMCRec_nonmatch_hUSS_INSIDE_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_OUTSIDE", "hMCRec_nonmatch_hUSS_OUTSIDE", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_OUTSIDE_1D", "hMCRec_nonmatch_hUSS_OUTSIDE_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_OUTSIDE_1D_2_3", "hMCRec_nonmatch_hUSS_OUTSIDE_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG", "hMCRec_nonmatch_hUSS_OUTSIDE_TRIG", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D", "hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3", "hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3", kTH1F, {MinvAxis}); + + // EVENT SELECTION eventSelection = jetderiveddatautilities::initialiseEventSelection(static_cast(cfgeventSelections)); @@ -188,7 +225,7 @@ struct phiInJets { using EventCandidates = soa::Join; // , aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs using TrackCandidates = soa::Join; - Filter jetCuts = aod::jet::pt > cfgjetPtMin&& aod::jet::r == nround(cfgjetR.node() * 100.0f); + Filter jetCuts = aod::jet::pt > cfgjetPtMin && aod::jet::r == nround(cfgjetR.node() * 100.0f); // Function for track quality cuts ///////////////////////////////////////////////////////////////////////////// @@ -197,7 +234,7 @@ struct phiInJets { bool trackSelection(const TrackType track) { // basic track cuts - if (std::abs(track.pt()) < cfgtrkMinPt) + if (track.pt() < cfgtrkMinPt) return false; if (std::abs(track.dcaXY()) > cfgMaxDCArToPVcut) @@ -212,20 +249,21 @@ struct phiInJets { if (track.tpcNClsFindable() < cfgnFindableTPCClusters) return false; - if (std::abs(track.tpcNClsCrossedRows()) < cfgnTPCCrossedRows) + if (track.tpcNClsCrossedRows() < cfgnTPCCrossedRows) return false; - if (std::abs(track.tpcCrossedRowsOverFindableCls()) > cfgnRowsOverFindable) + if (track.tpcCrossedRowsOverFindableCls() > cfgnRowsOverFindable) return false; - if (std::abs(track.tpcChi2NCl()) > cfgnTPCChi2) + if (track.tpcChi2NCl() > cfgnTPCChi2) return false; - if (std::abs(track.itsChi2NCl()) > cfgnITSChi2) + if (track.itsChi2NCl() > cfgnITSChi2) return false; if (cfgConnectedToPV && !track.isPVContributor()) return false; + return true; }; ///////////////////////////////////////////////////////////////////////////// @@ -369,21 +407,20 @@ struct phiInJets { ///////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////// - int nEvents = 0; void processJetTracks(aod::JCollision const& collision, soa::Filtered> const& chargedjets, soa::Join const& tracks, TrackCandidates const&) { if (cDebugLevel > 0) { nEvents++; if ((nEvents + 1) % 10000 == 0) - std::cout << nEvents << std::endl; + std::cout <<"Processed Data Events: "<< nEvents << std::endl; } JEhistos.fill(HIST("nEvents"), 0.5); - - if (!jetderiveddatautilities::selectCollision(collision, eventSelection)) - return; + if (fabs(collision.posZ()) > 10) return; + if (!jetderiveddatautilities::selectCollision(collision, jetderiveddatautilities::JCollisionSel::sel8)) + return; for (auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(tracks, tracks))) { auto trk1 = track1.track_as>(); @@ -395,19 +432,21 @@ struct phiInJets { JEhistos.fill(HIST("FJetaHistogram"), chargedjet.eta()); JEhistos.fill(HIST("FJphiHistogram"), chargedjet.phi()); JEhistos.fill(HIST("FJptHistogram"), chargedjet.pt()); - JEhistos.fill(HIST("FJnchHistogram"), chargedjet.tracks().size()); + // JEhistos.fill(HIST("FJnchHistogram"), chargedjet.tracks().size()); nJets++; } + JEhistos.fill(HIST("nJetsPerEvent"), nJets); JEhistos.fill(HIST("nEvents"), 1.5); + // return; - for (auto const& track : tracks) { - auto originalTrack = track.track_as>(); + for (auto& trackC : tracks) { + auto originalTrack = trackC.track_as>(); JEhistos.fill(HIST("hDCArToPv"), originalTrack.dcaXY()); JEhistos.fill(HIST("hDCAzToPv"), originalTrack.dcaZ()); JEhistos.fill(HIST("rawpT"), originalTrack.pt()); - JEhistos.fill(HIST("rawDpT"), track.pt(), track.pt() - originalTrack.pt()); + JEhistos.fill(HIST("rawDpT"), trackC.pt(), trackC.pt() - originalTrack.pt()); JEhistos.fill(HIST("hIsPrim"), originalTrack.isPrimaryTrack()); JEhistos.fill(HIST("hIsGood"), originalTrack.isGlobalTrackWoDCA()); JEhistos.fill(HIST("hIsPrimCont"), originalTrack.isPVContributor()); @@ -420,8 +459,8 @@ struct phiInJets { if (!trackSelection(originalTrack)) continue; - JEhistos.fill(HIST("etaHistogram"), track.eta()); - JEhistos.fill(HIST("phiHistogram"), track.phi()); + JEhistos.fill(HIST("etaHistogram"), trackC.eta()); + JEhistos.fill(HIST("phiHistogram"), trackC.phi()); } // JTrack Loop }; // Process Switch @@ -436,20 +475,50 @@ struct phiInJets { using myCompleteTracks = soa::Join; using myCompleteJetTracks = soa::Join; int nJEEvents = 0; - void processRec(o2::aod::JCollision const& collision, myCompleteJetTracks const& tracks, soa::Filtered const& mcdjets, aod::McParticles const&, myCompleteTracks const&) - { + int nprocessRecEvents = 0; + void processRec(o2::aod::JCollision const& collision, myCompleteJetTracks const& tracks, soa::Filtered const& mcdjets, aod::McParticles const&, myCompleteTracks const& originalTracks) + { if (cDebugLevel > 0) { - nJEEvents++; - if ((nJEEvents + 1) % 10000 == 0) - std::cout << "JEvents: " << nJEEvents << std::endl; + nprocessRecEvents++; + if ((nprocessRecEvents + 1) % 10000 == 0) + std::cout << "processRec: " << nprocessRecEvents << std::endl; } JEhistos.fill(HIST("nEvents_MCRec"), 0.5); if (fabs(collision.posZ()) > 10) return; - if (!jetderiveddatautilities::selectCollision(collision, eventSelection)) + if (!jetderiveddatautilities::selectCollision(collision, jetderiveddatautilities::JCollisionSel::sel8)) return; + + bool INELgt0 = false; + for(const auto& track : tracks){ + if(TMath::Abs(track.eta())<0.8){ + INELgt0 = true; + break; + } + } + if(!INELgt0) + return; + JEhistos.fill(HIST("nEvents_MCRec"), 1.5); + + std::vector mcd_pt{}; + std::vector mcd_phi{}; + std::vector mcd_eta{}; + + + bool hasJets = kFALSE; + for (auto& mcdjet : mcdjets) { + hasJets=kTRUE; + mcd_pt.push_back(mcdjet.pt()); + mcd_eta.push_back(mcdjet.eta()); + mcd_phi.push_back(mcdjet.phi()); + registry.fill(HIST("h_jet_pt"), mcdjet.pt()); + registry.fill(HIST("h_jet_eta"), mcdjet.eta()); + registry.fill(HIST("h_jet_phi"), mcdjet.phi()); + } + if(hasJets) + JEhistos.fill(HIST("nEvents_MCRec"), 2.5); // Track Eff for (const auto& track : tracks) { @@ -458,7 +527,7 @@ struct phiInJets { continue; if (track.has_mcParticle()) { auto mcParticle = track.mcParticle(); - if (mcParticle.isPhysicalPrimary() && fabs(mcParticle.y()) <= 0.5) { // do this in the context of the track ! (context matters!!!) + if (mcParticle.isPhysicalPrimary() && fabs(mcParticle.y()) <= 0.5) { if (abs(mcParticle.pdgCode()) == 211) JEhistos.fill(HIST("ptJEHistogramPion"), mcParticle.pt()); if (abs(mcParticle.pdgCode()) == 321) @@ -467,93 +536,149 @@ struct phiInJets { JEhistos.fill(HIST("ptJEHistogramProton"), mcParticle.pt()); } } - } - + for (const auto& track2 : tracks) { + auto originalTrack2 = track2.track_as(); + if (!trackSelection(originalTrack2)) + continue; + + if (originalTrack.index() >= originalTrack2.index()) + continue; + + if (fabs(originalTrack.eta()) > 0.8 || fabs(originalTrack2.eta()) > 0.8) + continue; + + //check PID + + if(track.has_mcParticle() && track2.has_mcParticle()){ + auto part1 = track.mcParticle(); + auto part2 = track2.mcParticle(); + if (fabs(part1.pdgCode()) != 321) + continue; // Not Kaon + if (fabs(part2.pdgCode()) != 321) + continue; // Not Kaon + if (!part1.has_mothers()) + continue; // Not decaying Kaon + if (!part2.has_mothers()) + continue; // Not decaying Kaon + + + std::vector mothers1{}; + std::vector mothers1PDG{}; + for (auto& part1_mom : part1.mothers_as()) { + mothers1.push_back(part1_mom.globalIndex()); + mothers1PDG.push_back(part1_mom.pdgCode()); + } + + std::vector mothers2{}; + std::vector mothers2PDG{}; + for (auto& part2_mom : part2.mothers_as()) { + mothers2.push_back(part2_mom.globalIndex()); + mothers2PDG.push_back(part2_mom.pdgCode()); + } + + + if (mothers1PDG[0] != 333) + continue; // mother not phi + if (mothers2PDG[0] != 333) + continue; // mother not phi + if (mothers1[0] != mothers2[0]) + continue; // Kaons not from the same phi + + TLorentzVector lDecayDaughter1, lDecayDaughter2, lResonance; + lDecayDaughter1.SetXYZM(originalTrack.px(), originalTrack.py(), originalTrack.pz(), massKa); + lDecayDaughter2.SetXYZM(originalTrack2.px(), originalTrack2.py(), originalTrack2.pz(), massKa); + lResonance = lDecayDaughter1 + lDecayDaughter2; + if (lResonance.Rapidity() > 0.5) + continue; + JEhistos.fill(HIST("ptJEHistogramPhi"), lResonance.Pt()); + + bool jetFlag = false; + for (int i = 0; i < mcd_pt.size(); i++) { + double phidiff = TVector2::Phi_mpi_pi(mcd_pt[i] - lResonance.Phi()); + double etadiff = mcd_eta[i] - lResonance.Eta(); + double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); + if (R < cfgjetR) + jetFlag = true; + } + + if (jetFlag) { + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_INSIDE_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_INSIDE_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_INSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag && mcd_pt.size()>0) { + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag) { + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + } //! jetflag + + if(hasJets) + JEhistos.fill(HIST("ptJEHistogramPhi_JetTrigger"), lResonance.Pt()); + JEhistos.fill(HIST("minvJEHistogramPhi"), lResonance.M()); + }//mcpart check + }//tracks2 + }//tracks1 // Jet Eff - for (auto& mcdjet : mcdjets) { - registry.fill(HIST("h_jet_pt"), mcdjet.pt()); - registry.fill(HIST("h_jet_eta"), mcdjet.eta()); - registry.fill(HIST("h_jet_phi"), mcdjet.phi()); - } - // Phi Eff - - for (auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(tracks, tracks))) { - auto trk1 = track1.track_as(); - auto trk2 = track2.track_as(); - if (trk1.index() >= trk2.index()) - continue; - if (fabs(trk1.eta()) > 0.8 || fabs(trk2.eta()) > 0.8) - continue; - - if (!trackSelection(trk1)) - continue; - if (!trackSelection(trk2)) - continue; - if ((trk1.sign() * trk2.sign()) > 0) - continue; // Not K+K- - if (!track1.has_mcParticle()) - continue; - if (!track2.has_mcParticle()) - continue; - - auto part1 = track1.mcParticle(); - auto part2 = track2.mcParticle(); - if (fabs(part1.pdgCode()) != 321) - continue; // Not Kaon - if (fabs(part2.pdgCode()) != 321) - continue; // Not Kaon - if (!part1.has_mothers()) - continue; // Not decaying Kaon - if (!part2.has_mothers()) - continue; // Not decaying Kaon - - std::vector mothers1{}; - std::vector mothers1PDG{}; - for (auto& part1_mom : part1.mothers_as()) { - mothers1.push_back(part1_mom.globalIndex()); - mothers1PDG.push_back(part1_mom.pdgCode()); - } - if (mothers1.size() > 1) - continue; // Strange multi-mother decay - - std::vector mothers2{}; - std::vector mothers2PDG{}; - for (auto& part2_mom : part2.mothers_as()) { - mothers2.push_back(part2_mom.globalIndex()); - mothers2PDG.push_back(part2_mom.pdgCode()); - } - - if (mothers2.size() > 1) - continue; // Strange multi-mother decay - if (mothers1PDG[0] != 333) - continue; // mother not phi - if (mothers2PDG[0] != 333) - continue; // mother not phi - if (mothers1[0] != mothers2[0]) - continue; // Kaons not from the same phi - TLorentzVector lDecayDaughter1, lDecayDaughter2, lResonance; - lDecayDaughter1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), massKa); - lDecayDaughter2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massKa); - lResonance = lDecayDaughter1 + lDecayDaughter2; - if (lResonance.Rapidity() > 0.5) - continue; - JEhistos.fill(HIST("ptJEHistogramPhi"), lResonance.Pt()); - JEhistos.fill(HIST("minvJEHistogramPhi"), lResonance.M()); - - } // end of phi check } PROCESS_SWITCH(phiInJets, processRec, "pikp detector level MC JE", true); //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + int nprocessSimEvents = 0; void processSim(o2::aod::JMcCollision const& collision, aod::JMcParticles const& mcParticles, soa::Filtered const& mcpjets) { + if (cDebugLevel > 0) { + nprocessSimEvents++; + if ((nprocessSimEvents + 1) % 10000 == 0) + std::cout << "processSim: " << nprocessSimEvents << std::endl; + } + JEhistos.fill(HIST("nEvents_MCGen"), 0.5); if (fabs(collision.posZ()) > 10) return; + bool INELgt0 = false; + for(const auto& mcParticle : mcParticles){ + if(TMath::Abs(mcParticle.eta())<0.8){ + INELgt0 = true; + break; + } + } + if(!INELgt0) + return; + + std::vector mcp_pt{}; + std::vector mcp_phi{}; + std::vector mcp_eta{}; JEhistos.fill(HIST("nEvents_MCGen"), 1.5); + bool hasJets = kFALSE; + // Check jets + for (auto& mcpjet : mcpjets) { + hasJets=kTRUE; + mcp_pt.push_back(mcpjet.pt()); + mcp_eta.push_back(mcpjet.eta()); + mcp_phi.push_back(mcpjet.phi()); + + registry.fill(HIST("h_part_jet_pt"), mcpjet.pt()); + registry.fill(HIST("h_part_jet_eta"), mcpjet.eta()); + registry.fill(HIST("h_part_jet_phi"), mcpjet.phi()); + } + if(hasJets) + JEhistos.fill(HIST("nEvents_MCGen"), 2.5); + // Check pikp and phi for (const auto& mcParticle : mcParticles) { if (mcParticle.isPhysicalPrimary() && fabs(mcParticle.y()) <= 0.5) { // watch out for context!!! @@ -567,19 +692,57 @@ struct phiInJets { if (fabs(mcParticle.y()) <= 0.5) { // watch out for context!!! TLorentzVector lResonance; lResonance.SetPxPyPzE(mcParticle.px(), mcParticle.py(), mcParticle.pz(), mcParticle.e()); - if (abs(mcParticle.pdgCode()) == 333) + if (abs(mcParticle.pdgCode()) == 333){ JEhistos.fill(HIST("ptGeneratedPhi"), mcParticle.pt()); - if (abs(mcParticle.pdgCode()) == 333) JEhistos.fill(HIST("mGeneratedPhi"), lResonance.M()); - } - } - // Check jets - for (auto& mcpjet : mcpjets) { - registry.fill(HIST("h_part_jet_pt"), mcpjet.pt()); - registry.fill(HIST("h_part_jet_eta"), mcpjet.eta()); - registry.fill(HIST("h_part_jet_phi"), mcpjet.phi()); - } - } + ////////////////////////////Implementation of phi finding + + TLorentzVector lResonance; + lResonance.SetPxPyPzE(mcParticle.px(), mcParticle.py(), mcParticle.pz(), mcParticle.e()); + bool jetFlag = false; + for (int i = 0; i < mcp_pt.size(); i++) { + double phidiff = TVector2::Phi_mpi_pi(mcp_pt[i] - lResonance.Phi()); + double etadiff = mcp_eta[i] - lResonance.Eta(); + double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); + if (R < cfgjetR) + jetFlag = true; + } + + if (jetFlag) { + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_INSIDE_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_INSIDE_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_INSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag && mcp_pt.size()>0) { + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag) { + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + + } //! jetflag + + ////////////////////////////Phi found + + + if(hasJets){ + JEhistos.fill(HIST("ptGeneratedPhi_JetTrigger"), mcParticle.pt()); + }//check for jets + + }//check for phi + }//check for rapidity + }//loop over particles + }//process switch PROCESS_SWITCH(phiInJets, processSim, "pikp particle level MC", true); //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -589,14 +752,38 @@ struct phiInJets { using JetMCDTable = soa::Filtered>; // void processMatchedGen(o2::aod::JMcCollision const& collision, aod::JMcParticles const& mcParticles, soa::Filtered const& mcpjets) - void processMatchedGen(aod::JMcCollision const&, - JetMCDTable const&, + int nprocessSimJEEvents=0; + void processMatchedGen(aod::JMcCollision const& collision, + JetMCDTable const& mcdjets, JetMCPTable const& mcpjets, aod::JMcParticles const& mcParticles) { - // if(fabs(collision.posZ())>10) - // return; + if (cDebugLevel > 0) { + nprocessSimJEEvents++; + if ((nprocessSimJEEvents + 1) % 10000 == 0) + std::cout << "JEvents: " << nprocessSimJEEvents << std::endl; + } + + JEhistos.fill(HIST("nEvents_MCGen_MATCHED"), 0.5); + + if(fabs(collision.posZ())>10) + return; + + if (fabs(collision.posZ()) > 10) + return; + bool INELgt0 = false; + for(const auto& mcParticle : mcParticles){ + if(TMath::Abs(mcParticle.eta())<0.8){ + INELgt0 = true; + break; + } + } + if(!INELgt0) + return; + + + JEhistos.fill(HIST("nEvents_MCGen_MATCHED"), 1.5); std::vector mcd_pt{}; std::vector mcd_phi{}; @@ -605,6 +792,8 @@ struct phiInJets { std::vector mcp_phi{}; std::vector mcp_eta{}; + bool hasJets=kFALSE; + for (auto& mcpjet : mcpjets) { if (!mcpjet.has_matchedJetGeo()) continue; @@ -612,7 +801,7 @@ struct phiInJets { for (auto& mcdjet : mcpjet.template matchedJetGeo_as()) { if (!mcdjet.has_matchedJetGeo()) continue; - + hasJets=kTRUE; registry.fill(HIST("h_matched_GEN_jet_pt"), mcpjet.pt(), mcdjet.pt() - mcpjet.pt()); registry.fill(HIST("h_matched_GEN_jet_phi"), mcpjet.phi(), mcdjet.phi() - mcpjet.phi()); registry.fill(HIST("h_matched_GEN_jet_eta"), mcpjet.eta(), mcdjet.eta() - mcpjet.eta()); @@ -626,6 +815,10 @@ struct phiInJets { } // mcpjets } // mcdjets + if(hasJets) + JEhistos.fill(HIST("nEvents_MCGen_MATCHED"), 2.5); + + // First we do GEN part for (const auto& mcParticle : mcParticles) { if (fabs(mcParticle.y()) > 0.5) @@ -651,16 +844,23 @@ struct phiInJets { JEhistos.fill(HIST("hMCTrue_hUSS_INSIDE_1D_2_3"), lResonance.M()); JEhistos.fill(HIST("hMCTrue_hUSS_INSIDE"), 1.0, lResonance.Pt(), lResonance.M()); - } // jetflag - if (!jetFlag) { - JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_1D"), lResonance.M()); + } else if (!jetFlag && mcp_pt.size()>0) { + JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_TRIG_1D"), lResonance.M()); if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) - JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_TRIG_1D_2_3"), lResonance.M()); - JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_TRIG"), 1.0, lResonance.Pt(), lResonance.M()); - } //! jetflag + } else if (!jetFlag) { + JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + + } //! jetflag } // chech for phi } // MC Particles } // main fcn @@ -668,7 +868,7 @@ struct phiInJets { //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - + int nprocessRecJEEvents=0; // void processMatchedRec(o2::aod::JCollision const& collision, myCompleteJetTracks const& tracks, soa::Filtered const& mcdjets, aod::McParticles const&, myCompleteTracks const& originalTracks) void processMatchedRec(aod::JCollision const& collision, JetMCDTable const& mcdjets, @@ -677,11 +877,31 @@ struct phiInJets { myCompleteTracks const&, aod::McParticles const&) { + + if (cDebugLevel > 0) { + nprocessRecJEEvents++; + if ((nprocessRecJEEvents + 1) % 10000 == 0) + std::cout << "JEvents: " << nprocessRecJEEvents << std::endl; + } + JEhistos.fill(HIST("nEvents_MCRec_MATCHED"), 0.5); + if (fabs(collision.posZ()) > 10) return; - if (!jetderiveddatautilities::selectCollision(collision, eventSelection)) + if (!jetderiveddatautilities::selectCollision(collision, jetderiveddatautilities::JCollisionSel::sel8)) + return; + + bool INELgt0 = false; + for(const auto& track : tracks){ + if(TMath::Abs(track.eta())<0.8){ + INELgt0 = true; + break; + } + } + if(!INELgt0) return; + JEhistos.fill(HIST("nEvents_MCRec_MATCHED"), 1.5); + std::vector mcd_pt{}; std::vector mcd_phi{}; std::vector mcd_eta{}; @@ -689,16 +909,14 @@ struct phiInJets { std::vector mcp_phi{}; std::vector mcp_eta{}; - // std::vector mcp_ID{}; - // std::vector mcp_ID{}; - + bool hasJets=kFALSE; for (auto& mcdjet : mcdjets) { if (!mcdjet.has_matchedJetGeo()) continue; for (auto& mcpjet : mcdjet.template matchedJetGeo_as()) { if (!mcpjet.has_matchedJetGeo()) continue; - + hasJets=kTRUE; registry.fill(HIST("h_matched_REC_jet_pt"), mcpjet.pt(), mcdjet.pt() - mcpjet.pt()); registry.fill(HIST("h_matched_REC_jet_phi"), mcpjet.phi(), mcdjet.phi() - mcpjet.phi()); registry.fill(HIST("h_matched_REC_jet_eta"), mcpjet.eta(), mcdjet.eta() - mcpjet.eta()); @@ -711,97 +929,104 @@ struct phiInJets { mcp_phi.push_back(mcpjet.phi()); } // mcpjets } // mcdjets - // Now we do REC part - for (auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(tracks, tracks))) { - auto trk1 = track1.track_as(); - auto trk2 = track2.track_as(); - if (trk1.index() >= trk2.index()) - continue; - if (fabs(trk1.eta()) > 0.8 || fabs(trk2.eta()) > 0.8) - continue; - if (!trackSelection(trk1)) - continue; - if (!trackSelection(trk2)) - continue; - if ((trk1.sign() * trk2.sign()) > 0) - continue; // Not K+K- - if (!track1.has_mcParticle()) - continue; - if (!track2.has_mcParticle()) - continue; + if(hasJets) + JEhistos.fill(HIST("nEvents_MCRec_MATCHED"), 2.5); - auto part1 = track1.mcParticle(); - auto part2 = track2.mcParticle(); - if (fabs(part1.pdgCode()) != 321) - continue; // Not Kaon - if (fabs(part2.pdgCode()) != 321) - continue; // Not Kaon - if (!part1.has_mothers()) - continue; // Not decaying Kaon - if (!part2.has_mothers()) - continue; // Not decaying Kaon - - std::vector mothers1{}; - std::vector mothers1PDG{}; - for (auto& part1_mom : part1.mothers_as()) { - mothers1.push_back(part1_mom.globalIndex()); - mothers1PDG.push_back(part1_mom.pdgCode()); - } - if (mothers1.size() > 1) - continue; // Strange multi-mother decay - - std::vector mothers2{}; - std::vector mothers2PDG{}; - for (auto& part2_mom : part2.mothers_as()) { - mothers2.push_back(part2_mom.globalIndex()); - mothers2PDG.push_back(part2_mom.pdgCode()); - } - if (mothers2.size() > 1) - continue; // Strange multi-mother decay - if (mothers1PDG[0] != 333) - continue; // mother not phi - if (mothers2PDG[0] != 333) - continue; // mother not phi - if (mothers1[0] != mothers2[0]) - continue; // Kaons not from the same phi - - TLorentzVector lDecayDaughter1, lDecayDaughter2, lResonance; - lDecayDaughter1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), massKa); - lDecayDaughter2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massKa); - lResonance = lDecayDaughter1 + lDecayDaughter2; - if (lResonance.Rapidity() > 0.5) - continue; - - bool jetFlag = false; - for (int i = 0; i < mcd_pt.size(); i++) { - double phidiff = TVector2::Phi_mpi_pi(mcd_pt[i] - lResonance.Phi()); - double etadiff = mcd_eta[i] - lResonance.Eta(); - double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); - if (R < cfgjetR) - jetFlag = true; - } - - if (jetFlag) { - JEhistos.fill(HIST("hMCRec_hUSS_INSIDE_1D"), lResonance.M()); - if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) - JEhistos.fill(HIST("hMCRec_hUSS_INSIDE_1D_2_3"), lResonance.M()); - JEhistos.fill(HIST("hMCRec_hUSS_INSIDE"), 1.0, lResonance.Pt(), lResonance.M()); - - } // jetflag - if (!jetFlag) { - JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_1D"), lResonance.M()); - - if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) - JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); - - JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); - - } //! jetflag - } // tracks + // for (auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(tracks, tracks))) { + for (const auto& track1 : tracks) { + auto trk1 = track1.track_as(); + for (const auto& track2 : tracks) { + auto trk2 = track2.track_as(); + if (trk1.index() >= trk2.index()) + continue; + if (fabs(trk1.eta()) > 0.8 || fabs(trk2.eta()) > 0.8) + continue; + if ((trk1.sign() * trk2.sign()) > 0) + continue; // Not K+K- + if(trackSelection(trk1) && trackSelection(trk2)){ + if(track1.has_mcParticle() && track2.has_mcParticle()){ + auto part1 = track1.mcParticle(); + auto part2 = track2.mcParticle(); + if (fabs(part1.pdgCode()) != 321) + continue; // Not Kaon + if (fabs(part2.pdgCode()) != 321) + continue; // Not Kaon + if (!part1.has_mothers()) + continue; // Not decaying Kaon + if (!part2.has_mothers()) + continue; // Not decaying Kaon + + std::vector mothers1{}; + std::vector mothers1PDG{}; + for (auto& part1_mom : part1.mothers_as()) { + mothers1.push_back(part1_mom.globalIndex()); + mothers1PDG.push_back(part1_mom.pdgCode()); + } + + std::vector mothers2{}; + std::vector mothers2PDG{}; + for (auto& part2_mom : part2.mothers_as()) { + mothers2.push_back(part2_mom.globalIndex()); + mothers2PDG.push_back(part2_mom.pdgCode()); + } + if (mothers1PDG[0] != 333) + continue; // mother not phi + if (mothers2PDG[0] != 333) + continue; // mother not phi + if (mothers1[0] != mothers2[0]) + continue; // Kaons not from the same phi + + TLorentzVector lDecayDaughter1, lDecayDaughter2, lResonance; + lDecayDaughter1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), massKa); + lDecayDaughter2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massKa); + lResonance = lDecayDaughter1 + lDecayDaughter2; + if (lResonance.Rapidity() > 0.5) + continue; + + bool jetFlag = false; + for (int i = 0; i < mcd_pt.size(); i++) { + double phidiff = TVector2::Phi_mpi_pi(mcd_pt[i] - lResonance.Phi()); + double etadiff = mcd_eta[i] - lResonance.Eta(); + double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); + if (R < cfgjetR) + jetFlag = true; + } + + if (jetFlag) { + JEhistos.fill(HIST("hMCRec_hUSS_INSIDE_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_hUSS_INSIDE_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hMCRec_hUSS_INSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag && mcd_pt.size()>0) { + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_TRIG_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_TRIG_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_TRIG"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag) { + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + + } //! jetflag + + }//pass track cut + }//has mc particle + + }//tracks 2 + }//tracks 1 + // } // tracks } // main fcn PROCESS_SWITCH(phiInJets, processMatchedRec, "phi matched Rec level MC", true); + }; // end of main struct WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From 9402316da7d1d19ff7801b570a90025d2df3a4ea Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Fri, 26 Apr 2024 11:13:27 +0000 Subject: [PATCH 2/2] Please consider the following formatting changes --- PWGJE/Tasks/phiInJets.cxx | 596 +++++++++++++++++++------------------- 1 file changed, 293 insertions(+), 303 deletions(-) diff --git a/PWGJE/Tasks/phiInJets.cxx b/PWGJE/Tasks/phiInJets.cxx index 98e828432c9..e3530bddfd0 100644 --- a/PWGJE/Tasks/phiInJets.cxx +++ b/PWGJE/Tasks/phiInJets.cxx @@ -186,7 +186,6 @@ struct phiInJets { JEhistos.add("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D", "hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D", kTH1F, {MinvAxis}); JEhistos.add("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3", "hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3", kTH1F, {MinvAxis}); - JEhistos.add("hMCRec_hUSS", "hMCRec_hUSS", kTH3F, {dRAxis, PtAxis, MinvAxis}); JEhistos.add("hMCRec_hUSS_1D", "hMCRec_hUSS_1D", kTH1F, {MinvAxis}); JEhistos.add("hMCRec_hUSS_1D_2_3", "hMCRec_hUSS_1D_2_3", kTH1F, {MinvAxis}); @@ -213,8 +212,6 @@ struct phiInJets { JEhistos.add("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D", "hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D", kTH1F, {MinvAxis}); JEhistos.add("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3", "hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3", kTH1F, {MinvAxis}); - - // EVENT SELECTION eventSelection = jetderiveddatautilities::initialiseEventSelection(static_cast(cfgeventSelections)); @@ -225,7 +222,7 @@ struct phiInJets { using EventCandidates = soa::Join; // , aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs using TrackCandidates = soa::Join; - Filter jetCuts = aod::jet::pt > cfgjetPtMin && aod::jet::r == nround(cfgjetR.node() * 100.0f); + Filter jetCuts = aod::jet::pt > cfgjetPtMin&& aod::jet::r == nround(cfgjetR.node() * 100.0f); // Function for track quality cuts ///////////////////////////////////////////////////////////////////////////// @@ -263,7 +260,7 @@ struct phiInJets { if (cfgConnectedToPV && !track.isPVContributor()) return false; - + return true; }; ///////////////////////////////////////////////////////////////////////////// @@ -413,10 +410,10 @@ struct phiInJets { if (cDebugLevel > 0) { nEvents++; if ((nEvents + 1) % 10000 == 0) - std::cout <<"Processed Data Events: "<< nEvents << std::endl; + std::cout << "Processed Data Events: " << nEvents << std::endl; } JEhistos.fill(HIST("nEvents"), 0.5); - + if (fabs(collision.posZ()) > 10) return; if (!jetderiveddatautilities::selectCollision(collision, jetderiveddatautilities::JCollisionSel::sel8)) @@ -435,11 +432,11 @@ struct phiInJets { // JEhistos.fill(HIST("FJnchHistogram"), chargedjet.tracks().size()); nJets++; } - + JEhistos.fill(HIST("nJetsPerEvent"), nJets); JEhistos.fill(HIST("nEvents"), 1.5); - // return; + // return; for (auto& trackC : tracks) { auto originalTrack = trackC.track_as>(); @@ -477,11 +474,11 @@ struct phiInJets { int nJEEvents = 0; int nprocessRecEvents = 0; void processRec(o2::aod::JCollision const& collision, myCompleteJetTracks const& tracks, soa::Filtered const& mcdjets, aod::McParticles const&, myCompleteTracks const& originalTracks) - { + { if (cDebugLevel > 0) { nprocessRecEvents++; if ((nprocessRecEvents + 1) % 10000 == 0) - std::cout << "processRec: " << nprocessRecEvents << std::endl; + std::cout << "processRec: " << nprocessRecEvents << std::endl; } JEhistos.fill(HIST("nEvents_MCRec"), 0.5); @@ -490,26 +487,25 @@ struct phiInJets { if (!jetderiveddatautilities::selectCollision(collision, jetderiveddatautilities::JCollisionSel::sel8)) return; - bool INELgt0 = false; - for(const auto& track : tracks){ - if(TMath::Abs(track.eta())<0.8){ - INELgt0 = true; - break; - } - } - if(!INELgt0) - return; - + bool INELgt0 = false; + for (const auto& track : tracks) { + if (TMath::Abs(track.eta()) < 0.8) { + INELgt0 = true; + break; + } + } + if (!INELgt0) + return; + JEhistos.fill(HIST("nEvents_MCRec"), 1.5); - + std::vector mcd_pt{}; std::vector mcd_phi{}; std::vector mcd_eta{}; - bool hasJets = kFALSE; for (auto& mcdjet : mcdjets) { - hasJets=kTRUE; + hasJets = kTRUE; mcd_pt.push_back(mcdjet.pt()); mcd_eta.push_back(mcdjet.eta()); mcd_phi.push_back(mcdjet.phi()); @@ -517,7 +513,7 @@ struct phiInJets { registry.fill(HIST("h_jet_eta"), mcdjet.eta()); registry.fill(HIST("h_jet_phi"), mcdjet.phi()); } - if(hasJets) + if (hasJets) JEhistos.fill(HIST("nEvents_MCRec"), 2.5); // Track Eff @@ -537,99 +533,97 @@ struct phiInJets { } } for (const auto& track2 : tracks) { - auto originalTrack2 = track2.track_as(); - if (!trackSelection(originalTrack2)) - continue; - - if (originalTrack.index() >= originalTrack2.index()) - continue; - - if (fabs(originalTrack.eta()) > 0.8 || fabs(originalTrack2.eta()) > 0.8) - continue; - - //check PID - - if(track.has_mcParticle() && track2.has_mcParticle()){ - auto part1 = track.mcParticle(); - auto part2 = track2.mcParticle(); - if (fabs(part1.pdgCode()) != 321) - continue; // Not Kaon - if (fabs(part2.pdgCode()) != 321) - continue; // Not Kaon - if (!part1.has_mothers()) - continue; // Not decaying Kaon - if (!part2.has_mothers()) - continue; // Not decaying Kaon - - - std::vector mothers1{}; - std::vector mothers1PDG{}; - for (auto& part1_mom : part1.mothers_as()) { - mothers1.push_back(part1_mom.globalIndex()); - mothers1PDG.push_back(part1_mom.pdgCode()); - } - - std::vector mothers2{}; - std::vector mothers2PDG{}; - for (auto& part2_mom : part2.mothers_as()) { - mothers2.push_back(part2_mom.globalIndex()); - mothers2PDG.push_back(part2_mom.pdgCode()); - } - - - if (mothers1PDG[0] != 333) - continue; // mother not phi - if (mothers2PDG[0] != 333) - continue; // mother not phi - if (mothers1[0] != mothers2[0]) - continue; // Kaons not from the same phi - - TLorentzVector lDecayDaughter1, lDecayDaughter2, lResonance; - lDecayDaughter1.SetXYZM(originalTrack.px(), originalTrack.py(), originalTrack.pz(), massKa); - lDecayDaughter2.SetXYZM(originalTrack2.px(), originalTrack2.py(), originalTrack2.pz(), massKa); - lResonance = lDecayDaughter1 + lDecayDaughter2; - if (lResonance.Rapidity() > 0.5) - continue; - JEhistos.fill(HIST("ptJEHistogramPhi"), lResonance.Pt()); - - bool jetFlag = false; - for (int i = 0; i < mcd_pt.size(); i++) { - double phidiff = TVector2::Phi_mpi_pi(mcd_pt[i] - lResonance.Phi()); - double etadiff = mcd_eta[i] - lResonance.Eta(); - double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); - if (R < cfgjetR) - jetFlag = true; - } - - if (jetFlag) { - JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_INSIDE_1D"), lResonance.M()); - if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) - JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_INSIDE_1D_2_3"), lResonance.M()); - JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_INSIDE"), 1.0, lResonance.Pt(), lResonance.M()); - - } else if (!jetFlag && mcd_pt.size()>0) { - JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D"), lResonance.M()); - - if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) - JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3"), lResonance.M()); - - JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG"), 1.0, lResonance.Pt(), lResonance.M()); - - } else if (!jetFlag) { - JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_1D"), lResonance.M()); - - if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) - JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); - - JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); - } //! jetflag - - if(hasJets) - JEhistos.fill(HIST("ptJEHistogramPhi_JetTrigger"), lResonance.Pt()); - JEhistos.fill(HIST("minvJEHistogramPhi"), lResonance.M()); - }//mcpart check - }//tracks2 - }//tracks1 + auto originalTrack2 = track2.track_as(); + if (!trackSelection(originalTrack2)) + continue; + + if (originalTrack.index() >= originalTrack2.index()) + continue; + + if (fabs(originalTrack.eta()) > 0.8 || fabs(originalTrack2.eta()) > 0.8) + continue; + + // check PID + + if (track.has_mcParticle() && track2.has_mcParticle()) { + auto part1 = track.mcParticle(); + auto part2 = track2.mcParticle(); + if (fabs(part1.pdgCode()) != 321) + continue; // Not Kaon + if (fabs(part2.pdgCode()) != 321) + continue; // Not Kaon + if (!part1.has_mothers()) + continue; // Not decaying Kaon + if (!part2.has_mothers()) + continue; // Not decaying Kaon + + std::vector mothers1{}; + std::vector mothers1PDG{}; + for (auto& part1_mom : part1.mothers_as()) { + mothers1.push_back(part1_mom.globalIndex()); + mothers1PDG.push_back(part1_mom.pdgCode()); + } + + std::vector mothers2{}; + std::vector mothers2PDG{}; + for (auto& part2_mom : part2.mothers_as()) { + mothers2.push_back(part2_mom.globalIndex()); + mothers2PDG.push_back(part2_mom.pdgCode()); + } + + if (mothers1PDG[0] != 333) + continue; // mother not phi + if (mothers2PDG[0] != 333) + continue; // mother not phi + if (mothers1[0] != mothers2[0]) + continue; // Kaons not from the same phi + + TLorentzVector lDecayDaughter1, lDecayDaughter2, lResonance; + lDecayDaughter1.SetXYZM(originalTrack.px(), originalTrack.py(), originalTrack.pz(), massKa); + lDecayDaughter2.SetXYZM(originalTrack2.px(), originalTrack2.py(), originalTrack2.pz(), massKa); + lResonance = lDecayDaughter1 + lDecayDaughter2; + if (lResonance.Rapidity() > 0.5) + continue; + JEhistos.fill(HIST("ptJEHistogramPhi"), lResonance.Pt()); + + bool jetFlag = false; + for (int i = 0; i < mcd_pt.size(); i++) { + double phidiff = TVector2::Phi_mpi_pi(mcd_pt[i] - lResonance.Phi()); + double etadiff = mcd_eta[i] - lResonance.Eta(); + double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); + if (R < cfgjetR) + jetFlag = true; + } + + if (jetFlag) { + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_INSIDE_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_INSIDE_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_INSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag && mcd_pt.size() > 0) { + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag) { + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + } //! jetflag + + if (hasJets) + JEhistos.fill(HIST("ptJEHistogramPhi_JetTrigger"), lResonance.Pt()); + JEhistos.fill(HIST("minvJEHistogramPhi"), lResonance.M()); + } // mcpart check + } // tracks2 + } // tracks1 // Jet Eff } PROCESS_SWITCH(phiInJets, processRec, "pikp detector level MC JE", true); @@ -642,32 +636,32 @@ struct phiInJets { if (cDebugLevel > 0) { nprocessSimEvents++; if ((nprocessSimEvents + 1) % 10000 == 0) - std::cout << "processSim: " << nprocessSimEvents << std::endl; + std::cout << "processSim: " << nprocessSimEvents << std::endl; } JEhistos.fill(HIST("nEvents_MCGen"), 0.5); if (fabs(collision.posZ()) > 10) return; - bool INELgt0 = false; - for(const auto& mcParticle : mcParticles){ - if(TMath::Abs(mcParticle.eta())<0.8){ - INELgt0 = true; - break; - } - } - if(!INELgt0) - return; - + bool INELgt0 = false; + for (const auto& mcParticle : mcParticles) { + if (TMath::Abs(mcParticle.eta()) < 0.8) { + INELgt0 = true; + break; + } + } + if (!INELgt0) + return; + std::vector mcp_pt{}; std::vector mcp_phi{}; std::vector mcp_eta{}; JEhistos.fill(HIST("nEvents_MCGen"), 1.5); - bool hasJets = kFALSE; + bool hasJets = kFALSE; // Check jets for (auto& mcpjet : mcpjets) { - hasJets=kTRUE; + hasJets = kTRUE; mcp_pt.push_back(mcpjet.pt()); mcp_eta.push_back(mcpjet.eta()); mcp_phi.push_back(mcpjet.phi()); @@ -676,7 +670,7 @@ struct phiInJets { registry.fill(HIST("h_part_jet_eta"), mcpjet.eta()); registry.fill(HIST("h_part_jet_phi"), mcpjet.phi()); } - if(hasJets) + if (hasJets) JEhistos.fill(HIST("nEvents_MCGen"), 2.5); // Check pikp and phi @@ -692,57 +686,56 @@ struct phiInJets { if (fabs(mcParticle.y()) <= 0.5) { // watch out for context!!! TLorentzVector lResonance; lResonance.SetPxPyPzE(mcParticle.px(), mcParticle.py(), mcParticle.pz(), mcParticle.e()); - if (abs(mcParticle.pdgCode()) == 333){ + if (abs(mcParticle.pdgCode()) == 333) { JEhistos.fill(HIST("ptGeneratedPhi"), mcParticle.pt()); JEhistos.fill(HIST("mGeneratedPhi"), lResonance.M()); - ////////////////////////////Implementation of phi finding - - TLorentzVector lResonance; - lResonance.SetPxPyPzE(mcParticle.px(), mcParticle.py(), mcParticle.pz(), mcParticle.e()); - bool jetFlag = false; - for (int i = 0; i < mcp_pt.size(); i++) { - double phidiff = TVector2::Phi_mpi_pi(mcp_pt[i] - lResonance.Phi()); - double etadiff = mcp_eta[i] - lResonance.Eta(); - double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); - if (R < cfgjetR) - jetFlag = true; - } - - if (jetFlag) { - JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_INSIDE_1D"), lResonance.M()); - if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) - JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_INSIDE_1D_2_3"), lResonance.M()); - JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_INSIDE"), 1.0, lResonance.Pt(), lResonance.M()); - - } else if (!jetFlag && mcp_pt.size()>0) { - JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D"), lResonance.M()); - - if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) - JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3"), lResonance.M()); - - JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG"), 1.0, lResonance.Pt(), lResonance.M()); - - } else if (!jetFlag) { - JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_1D"), lResonance.M()); - - if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) - JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); - - JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); - - } //! jetflag - - ////////////////////////////Phi found - - - if(hasJets){ - JEhistos.fill(HIST("ptGeneratedPhi_JetTrigger"), mcParticle.pt()); - }//check for jets - - }//check for phi - }//check for rapidity - }//loop over particles - }//process switch + ////////////////////////////Implementation of phi finding + + TLorentzVector lResonance; + lResonance.SetPxPyPzE(mcParticle.px(), mcParticle.py(), mcParticle.pz(), mcParticle.e()); + bool jetFlag = false; + for (int i = 0; i < mcp_pt.size(); i++) { + double phidiff = TVector2::Phi_mpi_pi(mcp_pt[i] - lResonance.Phi()); + double etadiff = mcp_eta[i] - lResonance.Eta(); + double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); + if (R < cfgjetR) + jetFlag = true; + } + + if (jetFlag) { + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_INSIDE_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_INSIDE_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_INSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag && mcp_pt.size() > 0) { + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag) { + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + + } //! jetflag + + ////////////////////////////Phi found + + if (hasJets) { + JEhistos.fill(HIST("ptGeneratedPhi_JetTrigger"), mcParticle.pt()); + } // check for jets + + } // check for phi + } // check for rapidity + } // loop over particles + } // process switch PROCESS_SWITCH(phiInJets, processSim, "pikp particle level MC", true); //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -752,7 +745,7 @@ struct phiInJets { using JetMCDTable = soa::Filtered>; // void processMatchedGen(o2::aod::JMcCollision const& collision, aod::JMcParticles const& mcParticles, soa::Filtered const& mcpjets) - int nprocessSimJEEvents=0; + int nprocessSimJEEvents = 0; void processMatchedGen(aod::JMcCollision const& collision, JetMCDTable const& mcdjets, JetMCPTable const& mcpjets, @@ -762,27 +755,26 @@ struct phiInJets { if (cDebugLevel > 0) { nprocessSimJEEvents++; if ((nprocessSimJEEvents + 1) % 10000 == 0) - std::cout << "JEvents: " << nprocessSimJEEvents << std::endl; + std::cout << "JEvents: " << nprocessSimJEEvents << std::endl; } JEhistos.fill(HIST("nEvents_MCGen_MATCHED"), 0.5); - if(fabs(collision.posZ())>10) + if (fabs(collision.posZ()) > 10) return; if (fabs(collision.posZ()) > 10) return; - bool INELgt0 = false; - for(const auto& mcParticle : mcParticles){ - if(TMath::Abs(mcParticle.eta())<0.8){ - INELgt0 = true; - break; - } - } - if(!INELgt0) - return; - - + bool INELgt0 = false; + for (const auto& mcParticle : mcParticles) { + if (TMath::Abs(mcParticle.eta()) < 0.8) { + INELgt0 = true; + break; + } + } + if (!INELgt0) + return; + JEhistos.fill(HIST("nEvents_MCGen_MATCHED"), 1.5); std::vector mcd_pt{}; @@ -792,7 +784,7 @@ struct phiInJets { std::vector mcp_phi{}; std::vector mcp_eta{}; - bool hasJets=kFALSE; + bool hasJets = kFALSE; for (auto& mcpjet : mcpjets) { if (!mcpjet.has_matchedJetGeo()) @@ -801,7 +793,7 @@ struct phiInJets { for (auto& mcdjet : mcpjet.template matchedJetGeo_as()) { if (!mcdjet.has_matchedJetGeo()) continue; - hasJets=kTRUE; + hasJets = kTRUE; registry.fill(HIST("h_matched_GEN_jet_pt"), mcpjet.pt(), mcdjet.pt() - mcpjet.pt()); registry.fill(HIST("h_matched_GEN_jet_phi"), mcpjet.phi(), mcdjet.phi() - mcpjet.phi()); registry.fill(HIST("h_matched_GEN_jet_eta"), mcpjet.eta(), mcdjet.eta() - mcpjet.eta()); @@ -815,10 +807,9 @@ struct phiInJets { } // mcpjets } // mcdjets - if(hasJets) + if (hasJets) JEhistos.fill(HIST("nEvents_MCGen_MATCHED"), 2.5); - // First we do GEN part for (const auto& mcParticle : mcParticles) { if (fabs(mcParticle.y()) > 0.5) @@ -844,7 +835,7 @@ struct phiInJets { JEhistos.fill(HIST("hMCTrue_hUSS_INSIDE_1D_2_3"), lResonance.M()); JEhistos.fill(HIST("hMCTrue_hUSS_INSIDE"), 1.0, lResonance.Pt(), lResonance.M()); - } else if (!jetFlag && mcp_pt.size()>0) { + } else if (!jetFlag && mcp_pt.size() > 0) { JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_TRIG_1D"), lResonance.M()); if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) @@ -853,14 +844,14 @@ struct phiInJets { JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_TRIG"), 1.0, lResonance.Pt(), lResonance.M()); } else if (!jetFlag) { - JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_1D"), lResonance.M()); + JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_1D"), lResonance.M()); - if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) - JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); - JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); - } //! jetflag + } //! jetflag } // chech for phi } // MC Particles } // main fcn @@ -868,7 +859,7 @@ struct phiInJets { //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - int nprocessRecJEEvents=0; + int nprocessRecJEEvents = 0; // void processMatchedRec(o2::aod::JCollision const& collision, myCompleteJetTracks const& tracks, soa::Filtered const& mcdjets, aod::McParticles const&, myCompleteTracks const& originalTracks) void processMatchedRec(aod::JCollision const& collision, JetMCDTable const& mcdjets, @@ -881,7 +872,7 @@ struct phiInJets { if (cDebugLevel > 0) { nprocessRecJEEvents++; if ((nprocessRecJEEvents + 1) % 10000 == 0) - std::cout << "JEvents: " << nprocessRecJEEvents << std::endl; + std::cout << "JEvents: " << nprocessRecJEEvents << std::endl; } JEhistos.fill(HIST("nEvents_MCRec_MATCHED"), 0.5); @@ -889,15 +880,15 @@ struct phiInJets { return; if (!jetderiveddatautilities::selectCollision(collision, jetderiveddatautilities::JCollisionSel::sel8)) return; - + bool INELgt0 = false; - for(const auto& track : tracks){ - if(TMath::Abs(track.eta())<0.8){ - INELgt0 = true; - break; - } + for (const auto& track : tracks) { + if (TMath::Abs(track.eta()) < 0.8) { + INELgt0 = true; + break; + } } - if(!INELgt0) + if (!INELgt0) return; JEhistos.fill(HIST("nEvents_MCRec_MATCHED"), 1.5); @@ -909,14 +900,14 @@ struct phiInJets { std::vector mcp_phi{}; std::vector mcp_eta{}; - bool hasJets=kFALSE; + bool hasJets = kFALSE; for (auto& mcdjet : mcdjets) { if (!mcdjet.has_matchedJetGeo()) continue; for (auto& mcpjet : mcdjet.template matchedJetGeo_as()) { if (!mcpjet.has_matchedJetGeo()) continue; - hasJets=kTRUE; + hasJets = kTRUE; registry.fill(HIST("h_matched_REC_jet_pt"), mcpjet.pt(), mcdjet.pt() - mcpjet.pt()); registry.fill(HIST("h_matched_REC_jet_phi"), mcpjet.phi(), mcdjet.phi() - mcpjet.phi()); registry.fill(HIST("h_matched_REC_jet_eta"), mcpjet.eta(), mcdjet.eta() - mcpjet.eta()); @@ -930,103 +921,102 @@ struct phiInJets { } // mcpjets } // mcdjets // Now we do REC part - if(hasJets) + if (hasJets) JEhistos.fill(HIST("nEvents_MCRec_MATCHED"), 2.5); // for (auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(tracks, tracks))) { for (const auto& track1 : tracks) { auto trk1 = track1.track_as(); for (const auto& track2 : tracks) { - auto trk2 = track2.track_as(); - if (trk1.index() >= trk2.index()) - continue; - if (fabs(trk1.eta()) > 0.8 || fabs(trk2.eta()) > 0.8) - continue; - if ((trk1.sign() * trk2.sign()) > 0) - continue; // Not K+K- - if(trackSelection(trk1) && trackSelection(trk2)){ - if(track1.has_mcParticle() && track2.has_mcParticle()){ - auto part1 = track1.mcParticle(); - auto part2 = track2.mcParticle(); - if (fabs(part1.pdgCode()) != 321) - continue; // Not Kaon - if (fabs(part2.pdgCode()) != 321) - continue; // Not Kaon - if (!part1.has_mothers()) - continue; // Not decaying Kaon - if (!part2.has_mothers()) - continue; // Not decaying Kaon - - std::vector mothers1{}; - std::vector mothers1PDG{}; - for (auto& part1_mom : part1.mothers_as()) { - mothers1.push_back(part1_mom.globalIndex()); - mothers1PDG.push_back(part1_mom.pdgCode()); - } - - std::vector mothers2{}; - std::vector mothers2PDG{}; - for (auto& part2_mom : part2.mothers_as()) { - mothers2.push_back(part2_mom.globalIndex()); - mothers2PDG.push_back(part2_mom.pdgCode()); - } - if (mothers1PDG[0] != 333) - continue; // mother not phi - if (mothers2PDG[0] != 333) - continue; // mother not phi - if (mothers1[0] != mothers2[0]) - continue; // Kaons not from the same phi - - TLorentzVector lDecayDaughter1, lDecayDaughter2, lResonance; - lDecayDaughter1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), massKa); - lDecayDaughter2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massKa); - lResonance = lDecayDaughter1 + lDecayDaughter2; - if (lResonance.Rapidity() > 0.5) - continue; - - bool jetFlag = false; - for (int i = 0; i < mcd_pt.size(); i++) { - double phidiff = TVector2::Phi_mpi_pi(mcd_pt[i] - lResonance.Phi()); - double etadiff = mcd_eta[i] - lResonance.Eta(); - double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); - if (R < cfgjetR) - jetFlag = true; - } - - if (jetFlag) { - JEhistos.fill(HIST("hMCRec_hUSS_INSIDE_1D"), lResonance.M()); - if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) - JEhistos.fill(HIST("hMCRec_hUSS_INSIDE_1D_2_3"), lResonance.M()); - JEhistos.fill(HIST("hMCRec_hUSS_INSIDE"), 1.0, lResonance.Pt(), lResonance.M()); - - } else if (!jetFlag && mcd_pt.size()>0) { - JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_TRIG_1D"), lResonance.M()); - - if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) - JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_TRIG_1D_2_3"), lResonance.M()); - - JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_TRIG"), 1.0, lResonance.Pt(), lResonance.M()); - - } else if (!jetFlag) { - JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_1D"), lResonance.M()); - - if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) - JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); - - JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); - - } //! jetflag - - }//pass track cut - }//has mc particle - - }//tracks 2 - }//tracks 1 + auto trk2 = track2.track_as(); + if (trk1.index() >= trk2.index()) + continue; + if (fabs(trk1.eta()) > 0.8 || fabs(trk2.eta()) > 0.8) + continue; + if ((trk1.sign() * trk2.sign()) > 0) + continue; // Not K+K- + if (trackSelection(trk1) && trackSelection(trk2)) { + if (track1.has_mcParticle() && track2.has_mcParticle()) { + auto part1 = track1.mcParticle(); + auto part2 = track2.mcParticle(); + if (fabs(part1.pdgCode()) != 321) + continue; // Not Kaon + if (fabs(part2.pdgCode()) != 321) + continue; // Not Kaon + if (!part1.has_mothers()) + continue; // Not decaying Kaon + if (!part2.has_mothers()) + continue; // Not decaying Kaon + + std::vector mothers1{}; + std::vector mothers1PDG{}; + for (auto& part1_mom : part1.mothers_as()) { + mothers1.push_back(part1_mom.globalIndex()); + mothers1PDG.push_back(part1_mom.pdgCode()); + } + + std::vector mothers2{}; + std::vector mothers2PDG{}; + for (auto& part2_mom : part2.mothers_as()) { + mothers2.push_back(part2_mom.globalIndex()); + mothers2PDG.push_back(part2_mom.pdgCode()); + } + if (mothers1PDG[0] != 333) + continue; // mother not phi + if (mothers2PDG[0] != 333) + continue; // mother not phi + if (mothers1[0] != mothers2[0]) + continue; // Kaons not from the same phi + + TLorentzVector lDecayDaughter1, lDecayDaughter2, lResonance; + lDecayDaughter1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), massKa); + lDecayDaughter2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massKa); + lResonance = lDecayDaughter1 + lDecayDaughter2; + if (lResonance.Rapidity() > 0.5) + continue; + + bool jetFlag = false; + for (int i = 0; i < mcd_pt.size(); i++) { + double phidiff = TVector2::Phi_mpi_pi(mcd_pt[i] - lResonance.Phi()); + double etadiff = mcd_eta[i] - lResonance.Eta(); + double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); + if (R < cfgjetR) + jetFlag = true; + } + + if (jetFlag) { + JEhistos.fill(HIST("hMCRec_hUSS_INSIDE_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_hUSS_INSIDE_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hMCRec_hUSS_INSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag && mcd_pt.size() > 0) { + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_TRIG_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_TRIG_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_TRIG"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag) { + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + + } //! jetflag + + } // pass track cut + } // has mc particle + + } // tracks 2 + } // tracks 1 // } // tracks } // main fcn PROCESS_SWITCH(phiInJets, processMatchedRec, "phi matched Rec level MC", true); - }; // end of main struct WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)