From da169da5bba64fab04e60a4aa4c625639ba62ec0 Mon Sep 17 00:00:00 2001 From: Adrian Nassirpour Date: Thu, 9 May 2024 14:51:15 +0900 Subject: [PATCH] Final functionality required for hyperloop usage --- PWGJE/Tasks/phiInJets.cxx | 94 ++++++++++++++++++++++++++------------- 1 file changed, 63 insertions(+), 31 deletions(-) diff --git a/PWGJE/Tasks/phiInJets.cxx b/PWGJE/Tasks/phiInJets.cxx index e3530bddfd0..52611da091d 100644 --- a/PWGJE/Tasks/phiInJets.cxx +++ b/PWGJE/Tasks/phiInJets.cxx @@ -50,7 +50,7 @@ struct phiInJets { HistogramRegistry JEhistos{"JEhistos", {}, OutputObjHandlingPolicy::AnalysisObject}; HistogramRegistry registry{"registry", - {{"h_jet_pt", "jet pT;#it{p}_{T,jet} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}}, + {{"h_jet_pt", "jet pT;#it{p}_{T,jet} (GeV/#it{c});entries", {HistType::kTH1F, {{4000, 0., 200.}}}}, {"h_jet_eta", "jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}}, {"h_jet_phi", "jet #phi;#phi_{jet};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}}, {"h_matched_REC_jet_pt", "matched_REC level jet pT;#it{p}_{T,jet part} (GeV/#it{c});Delta", {HistType::kTH2F, {{200, 0., 200.}, {400, -20., 20.}}}}, @@ -59,7 +59,7 @@ struct phiInJets { {"h_matched_GEN_jet_pt", "matched_GEN level jet pT;#it{p}_{T,jet part} (GeV/#it{c});Delta", {HistType::kTH2F, {{200, 0., 200.}, {400, -20., 20.}}}}, {"h_matched_GEN_jet_eta", "matched_GEN level jet #eta;#eta_{jet part};Delta", {HistType::kTH2F, {{100, -1.0, 1.0}, {400, -20., 20.}}}}, {"h_matched_GEN_jet_phi", "matched_GEN level jet #phi;#phi_{jet part};Delta", {HistType::kTH2F, {{80, -1.0, 7.}, {400, -20., 20.}}}}, - {"h_part_jet_pt", "particle level jet pT;#it{p}_{T,jet part} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}}, + {"h_part_jet_pt", "particle level jet pT;#it{p}_{T,jet part} (GeV/#it{c});entries", {HistType::kTH1F, {{4000, 0., 200.}}}}, {"h_part_jet_eta", "particle level jet #eta;#eta_{jet part};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}}, {"h_part_jet_phi", "particle level jet #phi;#phi_{jet part};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}}}}; @@ -67,6 +67,7 @@ struct phiInJets { Configurable cfgtrackSelections{"cfgtrackSelections", "globalTracks", "set track selections"}; Configurable cfgtrkMinPt{"cfgtrkMinPt", 0.15, "set track min pT"}; + Configurable cfgtrkMaxEta{"cfgtrkMaxEta", 0.9, "set track max Eta"}; Configurable cfgMaxDCArToPVcut{"cfgMaxDCArToPVcut", 0.5, "Track DCAr cut to PV Maximum"}; Configurable cfgMaxDCAzToPVcut{"cfgMaxDCAzToPVcut", 2.0, "Track DCAz cut to PV Maximum"}; Configurable cfgPrimaryTrack{"cfgPrimaryTrack", true, "Primary track selection"}; // kGoldenChi2 | kDCAxy | kDCAz @@ -131,7 +132,10 @@ struct phiInJets { 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("JetVsPhi_GEN","JetVsPhi_GEN",kTH2F,{{4000, 0., 200.},{200, 0, 20.0}}); + JEhistos.add("JetVsPhi_REC","JetVsPhi_REC",kTH2F,{{4000, 0., 200.},{200, 0, 20.0}}); + JEhistos.add("nJetsPerEvent", "nJetsPerEvent", kTH1F, {{10, 0.0, 10.0}}); JEhistos.add("hDCArToPv", "DCArToPv", kTH1F, {{300, 0.0, 3.0}}); @@ -234,6 +238,9 @@ struct phiInJets { if (track.pt() < cfgtrkMinPt) return false; + if (std::abs(track.eta()) > cfgtrkMaxEta) + return false; + if (std::abs(track.dcaXY()) > cfgMaxDCArToPVcut) return false; @@ -306,7 +313,7 @@ struct phiInJets { lDecayDaughter2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massKa); lResonance = lDecayDaughter1 + lDecayDaughter2; - if (std::abs(lResonance.Rapidity()) > 0.5) + if (std::abs(lResonance.Eta()) > cfgtrkMaxEta) return; ///////////////////////////////////////////////////////////////////////////// @@ -455,7 +462,7 @@ struct phiInJets { if (!trackSelection(originalTrack)) continue; - + JEhistos.fill(HIST("etaHistogram"), trackC.eta()); JEhistos.fill(HIST("phiHistogram"), trackC.phi()); } // JTrack Loop @@ -489,7 +496,7 @@ struct phiInJets { bool INELgt0 = false; for (const auto& track : tracks) { - if (TMath::Abs(track.eta()) < 0.8) { + if (fabs(track.eta()) < cfgtrkMaxEta) { INELgt0 = true; break; } @@ -523,7 +530,7 @@ struct phiInJets { continue; if (track.has_mcParticle()) { auto mcParticle = track.mcParticle(); - if (mcParticle.isPhysicalPrimary() && fabs(mcParticle.y()) <= 0.5) { + if (mcParticle.isPhysicalPrimary() && fabs(mcParticle.eta()) <= cfgtrkMaxEta) { if (abs(mcParticle.pdgCode()) == 211) JEhistos.fill(HIST("ptJEHistogramPion"), mcParticle.pt()); if (abs(mcParticle.pdgCode()) == 321) @@ -540,7 +547,7 @@ struct phiInJets { if (originalTrack.index() >= originalTrack2.index()) continue; - if (fabs(originalTrack.eta()) > 0.8 || fabs(originalTrack2.eta()) > 0.8) + if (fabs(originalTrack.eta()) > cfgtrkMaxEta || fabs(originalTrack2.eta()) > cfgtrkMaxEta) continue; // check PID @@ -582,13 +589,13 @@ struct phiInJets { 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) + if (fabs(lResonance.Eta()) > cfgtrkMaxEta) 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 phidiff = TVector2::Phi_mpi_pi(mcd_phi[i] - lResonance.Phi()); double etadiff = mcd_eta[i] - lResonance.Eta(); double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); if (R < cfgjetR) @@ -618,8 +625,13 @@ struct phiInJets { JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); } //! jetflag - if (hasJets) + if (hasJets){ JEhistos.fill(HIST("ptJEHistogramPhi_JetTrigger"), lResonance.Pt()); + auto triggerjet = std::min_element(mcd_pt.begin(), mcd_pt.end()); + double triggerjet_pt = *triggerjet; + JEhistos.fill(HIST("JetVsPhi_REC"), triggerjet_pt, lResonance.Pt()); + + } JEhistos.fill(HIST("minvJEHistogramPhi"), lResonance.M()); } // mcpart check } // tracks2 @@ -631,8 +643,9 @@ struct phiInJets { //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// int nprocessSimEvents = 0; - void processSim(o2::aod::JMcCollision const& collision, aod::JMcParticles const& mcParticles, soa::Filtered const& mcpjets) - { + // Preslice slice = o2::aod::JCollision::collisionId; + void processSim(o2::aod::JMcCollision const& collision, soa::SmallGroups> const& recocolls, aod::JMcParticles const& mcParticles, soa::Filtered const& mcpjets) + { if (cDebugLevel > 0) { nprocessSimEvents++; if ((nprocessSimEvents + 1) % 10000 == 0) @@ -641,18 +654,29 @@ struct phiInJets { JEhistos.fill(HIST("nEvents_MCGen"), 0.5); - if (fabs(collision.posZ()) > 10) + if(recocolls.size()<=0)//not reconstructed + return; + + for(auto& recocoll : recocolls){ //poorly reconstructed + if (!jetderiveddatautilities::selectCollision(recocoll, jetderiveddatautilities::JCollisionSel::sel8)) + return; + } + + + + if (fabs(collision.posZ()) > 10)//bad vertex return; bool INELgt0 = false; for (const auto& mcParticle : mcParticles) { - if (TMath::Abs(mcParticle.eta()) < 0.8) { + if (fabs(mcParticle.eta()) < cfgtrkMaxEta) { INELgt0 = true; break; } } - if (!INELgt0) + if (!INELgt0) //not INEL return; + std::vector mcp_pt{}; std::vector mcp_phi{}; std::vector mcp_eta{}; @@ -675,7 +699,7 @@ struct phiInJets { // Check pikp and phi for (const auto& mcParticle : mcParticles) { - if (mcParticle.isPhysicalPrimary() && fabs(mcParticle.y()) <= 0.5) { // watch out for context!!! + if (mcParticle.isPhysicalPrimary() && fabs(mcParticle.eta()) <= cfgtrkMaxEta) { // watch out for context!!! if (abs(mcParticle.pdgCode()) == 211) JEhistos.fill(HIST("ptGeneratedPion"), mcParticle.pt()); if (abs(mcParticle.pdgCode()) == 321) @@ -683,7 +707,7 @@ struct phiInJets { if (abs(mcParticle.pdgCode()) == 2212) JEhistos.fill(HIST("ptGeneratedProton"), mcParticle.pt()); } - if (fabs(mcParticle.y()) <= 0.5) { // watch out for context!!! + if (fabs(mcParticle.eta()) <= cfgtrkMaxEta) { // watch out for context!!! TLorentzVector lResonance; lResonance.SetPxPyPzE(mcParticle.px(), mcParticle.py(), mcParticle.pz(), mcParticle.e()); if (abs(mcParticle.pdgCode()) == 333) { @@ -695,7 +719,7 @@ struct phiInJets { 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 phidiff = TVector2::Phi_mpi_pi(mcp_phi[i] - lResonance.Phi()); double etadiff = mcp_eta[i] - lResonance.Eta(); double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); if (R < cfgjetR) @@ -727,9 +751,11 @@ struct phiInJets { } //! jetflag ////////////////////////////Phi found - if (hasJets) { JEhistos.fill(HIST("ptGeneratedPhi_JetTrigger"), mcParticle.pt()); + auto triggerjet = std::min_element(mcp_pt.begin(), mcp_pt.end()); + double triggerjet_pt = *triggerjet; + JEhistos.fill(HIST("JetVsPhi_GEN"), triggerjet_pt, mcParticle.pt()); } // check for jets } // check for phi @@ -747,6 +773,7 @@ struct phiInJets { // void processMatchedGen(o2::aod::JMcCollision const& collision, aod::JMcParticles const& mcParticles, soa::Filtered const& mcpjets) int nprocessSimJEEvents = 0; void processMatchedGen(aod::JMcCollision const& collision, + soa::SmallGroups> const& recocolls, JetMCDTable const& mcdjets, JetMCPTable const& mcpjets, aod::JMcParticles const& mcParticles) @@ -763,11 +790,17 @@ struct phiInJets { if (fabs(collision.posZ()) > 10) return; - if (fabs(collision.posZ()) > 10) + if(recocolls.size()<=0)//not reconstructed return; + for(auto& recocoll : recocolls){ //poorly reconstructed + if (!jetderiveddatautilities::selectCollision(recocoll, jetderiveddatautilities::JCollisionSel::sel8)) + return; + } + + bool INELgt0 = false; for (const auto& mcParticle : mcParticles) { - if (TMath::Abs(mcParticle.eta()) < 0.8) { + if (TMath::Abs(mcParticle.eta()) < cfgtrkMaxEta) { INELgt0 = true; break; } @@ -812,17 +845,16 @@ struct phiInJets { // First we do GEN part for (const auto& mcParticle : mcParticles) { - if (fabs(mcParticle.y()) > 0.5) - continue; - if (fabs(mcParticle.eta() > 0.8)) + if (fabs(mcParticle.eta()) > cfgtrkMaxEta) continue; - if (abs(mcParticle.pdgCode()) == 333) { + + if (fabs(mcParticle.pdgCode()) == 333) { 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 phidiff = TVector2::Phi_mpi_pi(mcp_phi[i] - lResonance.Phi()); double etadiff = mcp_eta[i] - lResonance.Eta(); double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); if (R < cfgjetR) @@ -883,7 +915,7 @@ struct phiInJets { bool INELgt0 = false; for (const auto& track : tracks) { - if (TMath::Abs(track.eta()) < 0.8) { + if (fabs(track.eta()) < cfgtrkMaxEta) { INELgt0 = true; break; } @@ -931,7 +963,7 @@ struct phiInJets { auto trk2 = track2.track_as(); if (trk1.index() >= trk2.index()) continue; - if (fabs(trk1.eta()) > 0.8 || fabs(trk2.eta()) > 0.8) + if (fabs(trk1.eta()) > cfgtrkMaxEta || fabs(trk2.eta()) > cfgtrkMaxEta) continue; if ((trk1.sign() * trk2.sign()) > 0) continue; // Not K+K- @@ -972,12 +1004,12 @@ struct phiInJets { 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) + if (fabs(lResonance.Eta()) > cfgtrkMaxEta) 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 phidiff = TVector2::Phi_mpi_pi(mcd_phi[i] - lResonance.Phi()); double etadiff = mcd_eta[i] - lResonance.Eta(); double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); if (R < cfgjetR)