diff --git a/PWGJE/Tasks/phiInJets.cxx b/PWGJE/Tasks/phiInJets.cxx index 6e1ffb0ab97..2c3d13f37e5 100644 --- a/PWGJE/Tasks/phiInJets.cxx +++ b/PWGJE/Tasks/phiInJets.cxx @@ -82,9 +82,9 @@ struct phiInJets { Configurable cfgnTOFPID{"cfgnTOFPID", 4, "nTOF PID"}; Configurable cfgjetPtMin{"cfgjetPtMin", 5.0, "minimum jet pT cut"}; Configurable cfgjetR{"cfgjetR", 0.4, "jet resolution parameter"}; - Configurable cfgVtxCut{"cfgVtxCut", 10.0, "V_z cut selection"}; + Configurable cfgVtxCut{"cfgVtxCut", 10.0, "V_z cut selection"}; Configurable cDebugLevel{"cDebugLevel", 0, "Resolution of Debug"}; - Configurable cfgBR{"cfgBR", false, "Forces Gen. Charged BR Only"}; + Configurable cfgBR{"cfgBR", false, "Forces Gen. Charged BR Only"}; // CONFIG DONE ///////////////////////////////////////// //INIT @@ -108,9 +108,9 @@ struct phiInJets { JEhistos.add("ptJEHistogramPhi_JetTrigger", "ptJEHistogramPhi_JetTrigger", kTH1F, {PtAxis}); JEhistos.add("minvJEHistogramPhi", "minvJEHistogramPhi", kTH1F, {MinvAxis}); - JEhistos.add("Resp_Matrix", "Resp_Matrix", HistType::kTHnSparseD , {PtAxis,axisPt,PtAxis,axisPt});//REC(Phi,Jet), GEN(Phi,Jet) - JEhistos.add("Resp_Matrix_MATCHED", "Resp_Matrix_MATCHED", HistType::kTHnSparseD , {PtAxis,axisPt,PtAxis,axisPt});//REC(Phi,Jet), GEN(Phi,Jet) - + JEhistos.add("Resp_Matrix", "Resp_Matrix", HistType::kTHnSparseD, {PtAxis, axisPt, PtAxis, axisPt}); // REC(Phi,Jet), GEN(Phi,Jet) + JEhistos.add("Resp_Matrix_MATCHED", "Resp_Matrix_MATCHED", HistType::kTHnSparseD, {PtAxis, axisPt, PtAxis, axisPt}); // REC(Phi,Jet), GEN(Phi,Jet) + JEhistos.add("ptGeneratedPion", "ptGeneratedPion", kTH1F, {PtAxis}); JEhistos.add("ptGeneratedKaon", "ptGeneratedKaon", kTH1F, {PtAxis}); JEhistos.add("ptGeneratedProton", "ptGeneratedProton", kTH1F, {PtAxis}); @@ -140,8 +140,10 @@ struct phiInJets { 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("hMCRec_nonmatch_hUSS_INSID_pt_v_eta", "hMCRec_nonmatch_hUSS_INSID_pt_v_eta", kTH2F, {PtAxis,axisEta}); - JEhistos.add("hMCGen_nonmatch_hUSS_INSID_pt_v_eta", "hMCGen_nonmatch_hUSS_INSID_pt_v_eta", kTH2F, {PtAxis,axisEta}); + JEhistos.add("hMCTrue_nonmatch_hUSS_Kangle_v_pt", "hMCTrue_nonmatch_hUSS_Kangle_v_pt", kTH2F, {axisEta, PtAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_Kangle_v_pt", "hMCRec_nonmatch_hUSS_Kangle_v_pt", kTH2F, {axisEta, PtAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_INSIDE_pt_v_eta", "hMCRec_nonmatch_hUSS_INSIDE_pt_v_eta", kTH2F, {PtAxis, axisEta}); + JEhistos.add("hMCTrue_nonmatch_hUSS_INSIDE_pt_v_eta", "hMCTrue_nonmatch_hUSS_INSIDE_pt_v_eta", kTH2F, {PtAxis, axisEta}); 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}}); @@ -488,7 +490,7 @@ struct phiInJets { using myCompleteJetTracks = soa::Join; 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) + void processRec(o2::aod::JCollision const& collision, myCompleteJetTracks const& tracks, soa::Filtered const& mcdjets, aod::McParticles const&, myCompleteTracks const& /*originalTracks*/) { if (cDebugLevel > 0) { nprocessRecEvents++; @@ -607,12 +609,23 @@ struct phiInJets { 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)); + + double phidiff_K1 = TVector2::Phi_mpi_pi(mcd_phi[i] - lDecayDaughter1.Phi()); + double etadiff_K1 = mcd_eta[i] - lDecayDaughter1.Eta(); + double R_K1 = TMath::Sqrt((etadiff_K1 * etadiff_K1) + (phidiff_K1 * phidiff_K1)); + double phidiff_K2 = TVector2::Phi_mpi_pi(mcd_phi[i] - lDecayDaughter2.Phi()); + double etadiff_K2 = mcd_eta[i] - lDecayDaughter2.Eta(); + double R_K2 = TMath::Sqrt((etadiff_K2 * etadiff_K2) + (phidiff_K2 * phidiff_K2)); + if(R < cfgjetR){ + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_Kangle_v_pt"), R_K1, lResonance.Pt()); + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_Kangle_v_pt"), R_K2, lResonance.Pt()); + } if (R < cfgjetR) jetFlag = true; } if (jetFlag) { - JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_INSID_pt_v_eta"), lResonance.Pt(), lResonance.Eta()); + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_INSIDE_pt_v_eta"), lResonance.Pt(), lResonance.Eta()); 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()); @@ -704,8 +717,8 @@ struct phiInJets { JEhistos.fill(HIST("nEvents_MCGen"), 2.5); // Check pikp and phi - for (const auto& mcParticle : mcParticles) { - if (mcParticle.isPhysicalPrimary() && fabs(mcParticle.eta()) <= cfgtrkMaxEta) { // watch out for context!!! + for (const auto& mcParticle : mcParticles) { + 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) @@ -717,20 +730,20 @@ struct phiInJets { TLorentzVector lResonance; lResonance.SetPxPyPzE(mcParticle.px(), mcParticle.py(), mcParticle.pz(), mcParticle.e()); if (abs(mcParticle.pdgCode()) == 333) { - JEhistos.fill(HIST("ptGeneratedPhi_ALLBR"), mcParticle.pt()); - - bool skip =false; - //First we check for Forced BR - if(mcParticle.has_daughters()) - for (auto& dgth : mcParticle.daughters_as()) - if(fabs(dgth.pdgCode())!= 321) - skip=true; - - if(skip && cfgBR) - continue; - JEhistos.fill(HIST("ptGeneratedPhi"), mcParticle.pt()); + JEhistos.fill(HIST("ptGeneratedPhi_ALLBR"), mcParticle.pt()); + + bool skip = false; + // First we check for Forced BR + if (mcParticle.has_daughters()) + for (auto& dgth : mcParticle.daughters_as()) + if (fabs(dgth.pdgCode()) != 321) + skip = true; + + if (skip && cfgBR) + continue; + JEhistos.fill(HIST("ptGeneratedPhi"), mcParticle.pt()); JEhistos.fill(HIST("mGeneratedPhi"), lResonance.M()); - + ////////////////////////////Implementation of phi finding TLorentzVector lResonance; @@ -740,13 +753,22 @@ struct phiInJets { 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 (mcParticle.has_daughters()){ + for (auto& dgth : mcParticle.daughters_as()){ + double phidiff_K = TVector2::Phi_mpi_pi(mcp_phi[i] - dgth.phi()); + double etadiff_K = mcp_eta[i] - dgth.eta(); + double R_K = TMath::Sqrt((etadiff_K * etadiff_K) + (phidiff_K * phidiff_K)); + if(R < cfgjetR) + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_Kangle_v_pt"), R_K, lResonance.Pt()); + } + } if (R < cfgjetR) jetFlag = true; } if (jetFlag) { - JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_INSID_pt_v_eta"), lResonance.Pt(), lResonance.Eta()); - JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_INSIDE_1D"), lResonance.M()); + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_INSIDE_pt_v_eta"), lResonance.Pt(), lResonance.Eta()); + 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()); @@ -793,7 +815,7 @@ struct phiInJets { int nprocessSimJEEvents = 0; void processMatchedGen(aod::JMcCollision const& collision, soa::SmallGroups> const& recocolls, - JetMCDTable const& mcdjets, + JetMCDTable const& /*mcdjets*/, JetMCPTable const& mcpjets, aod::JMcParticles const& mcParticles) @@ -867,15 +889,15 @@ struct phiInJets { continue; if (fabs(mcParticle.pdgCode()) == 333) { - bool skip =false; - //First we check for Forced BR - if(mcParticle.has_daughters()) - for (auto& dgth : mcParticle.daughters_as()) - if(fabs(dgth.pdgCode())!= 321) - skip=true; - - if(skip && cfgBR) - continue; + bool skip = false; + // First we check for Forced BR + if (mcParticle.has_daughters()) + for (auto& dgth : mcParticle.daughters_as()) + if (fabs(dgth.pdgCode()) != 321) + skip = true; + + if (skip && cfgBR) + continue; TLorentzVector lResonance; lResonance.SetPxPyPzE(mcParticle.px(), mcParticle.py(), mcParticle.pz(), mcParticle.e()); bool jetFlag = false; @@ -892,7 +914,7 @@ struct phiInJets { if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) 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) { JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_TRIG_1D"), lResonance.M()); @@ -1007,20 +1029,20 @@ struct phiInJets { std::vector mothers1{}; std::vector mothers1PDG{}; - std::vector mothers1Pt{}; + std::vector mothers1Pt{}; for (auto& part1_mom : part1.mothers_as()) { mothers1.push_back(part1_mom.globalIndex()); mothers1PDG.push_back(part1_mom.pdgCode()); - mothers1Pt.push_back(part1_mom.pt()); + mothers1Pt.push_back(part1_mom.pt()); } std::vector mothers2{}; std::vector mothers2PDG{}; - std::vector mothers2Pt{}; + std::vector mothers2Pt{}; for (auto& part2_mom : part2.mothers_as()) { mothers2.push_back(part2_mom.globalIndex()); mothers2PDG.push_back(part2_mom.pdgCode()); - mothers2Pt.push_back(part2_mom.pt()); + mothers2Pt.push_back(part2_mom.pt()); } if (mothers1PDG[0] != 333) continue; // mother not phi @@ -1044,21 +1066,21 @@ struct phiInJets { double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); if (R < cfgjetR) jetFlag = true; - - if(jetFlag){ //Fill Resp. Matrix - if(cDebugLevel > 0 ){ - std::cout<<"******************************************"< 0) { + std::cout << "******************************************" << std::endl; + std::cout << "Rec. Phi Pt: " << lResonance.Pt() << std::endl; + std::cout << "Rec. Jet Pt: " << mcd_pt[i] << std::endl; + std::cout << "Gen. Phi Pt: " << mothers1Pt[0] << std::endl; + std::cout << "Gen. Jet Pt: " << mcp_pt[i] << std::endl; + std::cout << "******************************************" << std::endl; + } + JEhistos.fill(HIST("Resp_Matrix_MATCHED"), lResonance.Pt(), mcd_pt[i], mothers1Pt[0], mcp_pt[i]); + } } - if (jetFlag) { + 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());