From d3035c60597f1d676df737867627cd5ed4d568f8 Mon Sep 17 00:00:00 2001 From: ariffero Date: Fri, 30 Aug 2024 17:28:16 +0200 Subject: [PATCH 1/6] Add phi to the variables in the tree --- PWGUD/Tasks/fwdMuonsUPC.cxx | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/PWGUD/Tasks/fwdMuonsUPC.cxx b/PWGUD/Tasks/fwdMuonsUPC.cxx index 34075b7a862..cb18088572b 100644 --- a/PWGUD/Tasks/fwdMuonsUPC.cxx +++ b/PWGUD/Tasks/fwdMuonsUPC.cxx @@ -20,6 +20,7 @@ #include "TLorentzVector.h" #include "TSystem.h" #include "TMath.h" +#include "TRandom3.h" namespace dimu { @@ -28,6 +29,8 @@ DECLARE_SOA_COLUMN(M, m, float); DECLARE_SOA_COLUMN(Pt, pt, float); DECLARE_SOA_COLUMN(Rap, rap, float); DECLARE_SOA_COLUMN(Phi, phi, float); +DECLARE_SOA_COLUMN(PhiAv, phiAv, float); +DECLARE_SOA_COLUMN(PhiCh, phiCh, float); // tracks positive (p) and negative (n) DECLARE_SOA_COLUMN(Ptp, ptp, float); DECLARE_SOA_COLUMN(Etap, etap, float); @@ -47,6 +50,7 @@ namespace o2::aod { DECLARE_SOA_TABLE(DiMu, "AOD", "DIMU", dimu::M, dimu::Pt, dimu::Rap, dimu::Phi, + dimu::PhiAv, dimu::PhiCh, dimu::Ptp, dimu::Etap, dimu::Phip, dimu::Ptn, dimu::Etan, dimu::Phin, dimu::Tzna, dimu::Ezna, dimu::Tznc, dimu::Eznc, dimu::Nclass); @@ -167,6 +171,8 @@ struct fwdMuonsUPC { registry.add("hTimeZNC", "ZNC Times;;#counts", kTH1D, {axisTimeZN}); registry.add("hEnergyZN", "ZNA vs ZNC energy", kTH2D, {axisEnergyZNA, axisEnergyZNC}); registry.add("hSameSign", "hSameSign;;#counts", kTH1D, {{6, -0.5, 5.5}}); + registry.add("hPhiCharge", "#phi #it{charge}", kTH1D,{axisPhi}); + registry.add("hPhiAverage", "#phi #it{average}", kTH1D,{axisPhi}); reg0n0n.add("hMass", "Ivariant mass of muon pairs - 0n0n;;#counts", kTH1D, {axisMass}); reg0n0n.add("hPt", "Transverse momentum mass of muon pairs - 0n0n;;#counts", kTH1D, {axisPt}); @@ -314,6 +320,24 @@ struct fwdMuonsUPC { if (p.Rapidity() > yCutUp) return; + //compute phi for azimuth anisotropy + TLorentzVector tSum, tDiffAv, tDiffCh; + tSum = p1 + p2; + if(tr1.sign()>0){ + tDiffCh = p1 - p2; + if(gRandom->Rndm()>0.5) + tDiffAv = p1 - p2; + else tDiffAv = p2 - p1; + }else{ + tDiffCh = p2 - p1; + if(gRandom->Rndm()>0.5) + tDiffAv = p2 - p1; + else tDiffAv = p1 - p2; + } + //average + float phiAverage = tSum.DeltaPhi(tDiffAv); + float phiCharge = tSum.DeltaPhi(tDiffCh); + // zdc info if (TMath::Abs(zdc.timeA) < 10) registry.fill(HIST("hTimeZNA"), zdc.timeA); @@ -381,15 +405,19 @@ struct fwdMuonsUPC { registry.fill(HIST("hPhi"), p.Phi()); registry.fill(HIST("hCharge"), tr1.sign()); registry.fill(HIST("hCharge"), tr2.sign()); + registry.fill(HIST("hPhiAverage"), phiAverage); + registry.fill(HIST("hPhiCharge"), phiCharge); // store the event to save it into a tree if (tr1.sign() > 0) { dimuSel(p.M(), p.Pt(), p.Rapidity(), p.Phi(), + phiAverage, phiCharge, p1.Pt(), p1.PseudoRapidity(), p1.Phi(), p2.Pt(), p2.PseudoRapidity(), p2.Phi(), zdc.timeA, zdc.enA, zdc.timeC, zdc.enC, znClass); } else { dimuSel(p.M(), p.Pt(), p.Rapidity(), p.Phi(), + phiAverage, phiCharge, p2.Pt(), p2.PseudoRapidity(), p2.Phi(), p1.Pt(), p1.PseudoRapidity(), p1.Phi(), zdc.timeA, zdc.enA, zdc.timeC, zdc.enC, znClass); From 90ffe9f08e50c3a0cd91f3e33828f6a93ebaeb29 Mon Sep 17 00:00:00 2001 From: ariffero Date: Sun, 20 Oct 2024 17:10:38 +0200 Subject: [PATCH 2/6] Change values of znClass for the different classes --- PWGUD/Tasks/fwdMuonsUPC.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/PWGUD/Tasks/fwdMuonsUPC.cxx b/PWGUD/Tasks/fwdMuonsUPC.cxx index cb18088572b..d9558d5cd8f 100644 --- a/PWGUD/Tasks/fwdMuonsUPC.cxx +++ b/PWGUD/Tasks/fwdMuonsUPC.cxx @@ -363,7 +363,7 @@ struct fwdMuonsUPC { // fill the histos in neutron classes and assign neutron class label // 0n0n if (neutron_C == false && neutron_A == false) { - znClass = 0; + znClass = 1; reg0n0n.fill(HIST("hMass"), p.M()); reg0n0n.fill(HIST("hPt"), p.Pt()); reg0n0n.fill(HIST("hPtFit"), p.Pt()); @@ -371,16 +371,16 @@ struct fwdMuonsUPC { reg0n0n.fill(HIST("hRapidity"), p.Rapidity()); } else if (neutron_A ^ neutron_C) { // Xn0n + 0nXn if (neutron_A) - znClass = 1; - else if (neutron_C) znClass = 2; + else if (neutron_C) + znClass = 3; regXn0n.fill(HIST("hMass"), p.M()); regXn0n.fill(HIST("hPt"), p.Pt()); regXn0n.fill(HIST("hPtFit"), p.Pt()); regXn0n.fill(HIST("hEta"), p.Eta()); regXn0n.fill(HIST("hRapidity"), p.Rapidity()); } else if (neutron_A && neutron_C) { // XnXn - znClass = 3; + znClass = 4; regXnXn.fill(HIST("hMass"), p.M()); regXnXn.fill(HIST("hPt"), p.Pt()); regXnXn.fill(HIST("hPtFit"), p.Pt()); From 50a3b2f91803db5ca176e42cf82d5796b10ab71f Mon Sep 17 00:00:00 2001 From: ariffero Date: Tue, 22 Oct 2024 09:24:47 +0200 Subject: [PATCH 3/6] Add cuts on pt, mass, and Y selected from config. --- PWGUD/Tasks/fwdMuonsUPC.cxx | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/PWGUD/Tasks/fwdMuonsUPC.cxx b/PWGUD/Tasks/fwdMuonsUPC.cxx index d9558d5cd8f..2a2be73d8d2 100644 --- a/PWGUD/Tasks/fwdMuonsUPC.cxx +++ b/PWGUD/Tasks/fwdMuonsUPC.cxx @@ -304,8 +304,21 @@ struct fwdMuonsUPC { if (nMIDs != 2) return; - // cuts on pair kinematics - if (!(p.M() > 2 && p.M() < 6 && p.Pt() < 5)) + //cut on pair kinematics + //select mass + if(p.M() < lowMass) + return; + if(p.M() > highMass) + return; + //select pt + if(p.Pt() < lowPt) + return; + if(p.Pt() > highPt) + return; + // select rapidity + if (p.Rapidity() < yCutLow) + return; + if (p.Rapidity() > yCutUp) return; // select opposite charge events only @@ -314,12 +327,6 @@ struct fwdMuonsUPC { return; } - // select rapidity ranges - if (p.Rapidity() < yCutLow) - return; - if (p.Rapidity() > yCutUp) - return; - //compute phi for azimuth anisotropy TLorentzVector tSum, tDiffAv, tDiffCh; tSum = p1 + p2; From 37bc94145acd17ecc58ba5a9cf0d8279c6406fa0 Mon Sep 17 00:00:00 2001 From: ariffero Date: Wed, 20 Nov 2024 18:22:42 +0100 Subject: [PATCH 4/6] Modify fwdMuonUPC --- PWGUD/Tasks/CMakeLists.txt | 10 ++ PWGUD/Tasks/fwdMuonsUPC.cxx | 346 ++++++++++++++++++++++++++++++++---- 2 files changed, 323 insertions(+), 33 deletions(-) diff --git a/PWGUD/Tasks/CMakeLists.txt b/PWGUD/Tasks/CMakeLists.txt index a0f1eb17747..a1b9c1b90e5 100644 --- a/PWGUD/Tasks/CMakeLists.txt +++ b/PWGUD/Tasks/CMakeLists.txt @@ -212,3 +212,13 @@ o2physics_add_dpl_workflow(event-by-event SOURCES eventByevent.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGPIDSelector COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(mumu-forward + SOURCES mumuForward.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGPIDSelector + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(fwd-muons-mc + SOURCES fwdMuonsMC.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGPIDSelector + COMPONENT_NAME Analysis) diff --git a/PWGUD/Tasks/fwdMuonsUPC.cxx b/PWGUD/Tasks/fwdMuonsUPC.cxx index 2a2be73d8d2..770648a3149 100644 --- a/PWGUD/Tasks/fwdMuonsUPC.cxx +++ b/PWGUD/Tasks/fwdMuonsUPC.cxx @@ -17,6 +17,7 @@ #include "DataFormatsParameters/GRPECSObject.h" #include "PWGUD/DataModel/UDTables.h" +#include "TDatabasePDG.h" #include "TLorentzVector.h" #include "TSystem.h" #include "TMath.h" @@ -26,15 +27,27 @@ namespace dimu { // dimuon DECLARE_SOA_COLUMN(M, m, float); +DECLARE_SOA_COLUMN(E, energy, float); +DECLARE_SOA_COLUMN(Px, px, float); +DECLARE_SOA_COLUMN(Py, py, float); +DECLARE_SOA_COLUMN(Pz, pz, float); DECLARE_SOA_COLUMN(Pt, pt, float); DECLARE_SOA_COLUMN(Rap, rap, float); DECLARE_SOA_COLUMN(Phi, phi, float); DECLARE_SOA_COLUMN(PhiAv, phiAv, float); DECLARE_SOA_COLUMN(PhiCh, phiCh, float); // tracks positive (p) and negative (n) +DECLARE_SOA_COLUMN(Ep, energyp, float); +DECLARE_SOA_COLUMN(Pxp, pxp, float); +DECLARE_SOA_COLUMN(Pyp, pyp, float); +DECLARE_SOA_COLUMN(Pzp, pzp, float); DECLARE_SOA_COLUMN(Ptp, ptp, float); DECLARE_SOA_COLUMN(Etap, etap, float); DECLARE_SOA_COLUMN(Phip, phip, float); +DECLARE_SOA_COLUMN(En, energyn, float); +DECLARE_SOA_COLUMN(Pxn, pxn, float); +DECLARE_SOA_COLUMN(Pyn, pyn, float); +DECLARE_SOA_COLUMN(Pzn, pzn, float); DECLARE_SOA_COLUMN(Ptn, ptn, float); DECLARE_SOA_COLUMN(Etan, etan, float); DECLARE_SOA_COLUMN(Phin, phin, float); @@ -49,10 +62,10 @@ DECLARE_SOA_COLUMN(Nclass, nclass, int); namespace o2::aod { DECLARE_SOA_TABLE(DiMu, "AOD", "DIMU", - dimu::M, dimu::Pt, dimu::Rap, dimu::Phi, + dimu::M, dimu::E, dimu::Px, dimu::Py, dimu::Pz, dimu::Pt, dimu::Rap, dimu::Phi, dimu::PhiAv, dimu::PhiCh, - dimu::Ptp, dimu::Etap, dimu::Phip, - dimu::Ptn, dimu::Etan, dimu::Phin, + dimu::Ep, dimu::Pxp, dimu::Pyp, dimu::Pzp, dimu::Ptp, dimu::Etap, dimu::Phip, + dimu::En, dimu::Pxn, dimu::Pyn, dimu::Pzn, dimu::Ptn, dimu::Etan, dimu::Phin, dimu::Tzna, dimu::Ezna, dimu::Tznc, dimu::Eznc, dimu::Nclass); } @@ -63,6 +76,7 @@ using namespace o2::framework::expressions; // defining constants double mMu = 0.10566; // mass of muon + // constants used in the track selection const float kRAbsMin = 17.6; const float kRAbsMid = 26.5; @@ -77,20 +91,27 @@ struct fwdMuonsUPC { using CandidatesFwd = soa::Join; using ForwardTracks = soa::Join; + using CompleteCandidatesFwd = soa::Join; + using CompleteFwdTracks = soa::Join; Produces dimuSel; - // defining histograms using histogram registry + // defining histograms using histogram registry: different histos for the different process functions HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; HistogramRegistry reg0n0n{"reg0n0n", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; HistogramRegistry regXn0n{"regXn0n", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; HistogramRegistry regXnXn{"regXnXn", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + HistogramRegistry McGenRegistry{"McGenRegistry", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + HistogramRegistry McCorrReg{"McCorrReg", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + + + // CONFIGURABLES // pT of muon pairs Configurable nBinsPt{"nBinsPt", 250, "N bins in pT histo"}; Configurable lowPt{"lowPt", 0., "lower limit in pT histo"}; - Configurable highPt{"highPt", 0.5, "upper limit in pT histo"}; + Configurable highPt{"highPt", 2, "upper limit in pT histo"}; // mass of muon pairs Configurable nBinsMass{"nBinsMass", 500, "N bins in mass histo"}; Configurable lowMass{"lowMass", 0., "lower limit in mass histo"}; @@ -125,10 +146,10 @@ struct fwdMuonsUPC { Configurable highEnZN{"highEnZN", 250., "upper limit in ZN energy histo"}; // configuarble rapidity cuts - Configurable yCutLow{"yCutLow", -4, "Lower cut in pair rapidity"}; - Configurable yCutUp{"yCutUp", -2.5, "Upper cut in pair rapidity"}; + //Configurable yCutLow{"yCutLow", -4, "Lower cut in pair rapidity"}; + //Configurable yCutUp{"yCutUp", -2.5, "Upper cut in pair rapidity"}; - void init(InitContext&) + void init(InitContext& context) { // binning of pT axis fr fit std::vector ptFitBinning = { @@ -152,9 +173,10 @@ struct fwdMuonsUPC { const AxisSpec axisPhiSingle{nBinsPhiSingle, lowPhiSingle, highPhiSingle, "#varphi_{trk}"}; // histos + // data and reco MC registry.add("hMass", "Ivariant mass of muon pairs;;#counts", kTH1D, {axisMass}); - registry.add("hPt", "Transverse momentum mass of muon pairs;;#counts", kTH1D, {axisPt}); - registry.add("hPtFit", "Transverse momentum mass of muon pairs;;#counts", kTH1D, {axisPtFit}); + registry.add("hPt", "Transverse momentum of muon pairs;;#counts", kTH1D, {axisPt}); + registry.add("hPtFit", "Transverse momentum of muon pairs;;#counts", kTH1D, {axisPtFit}); registry.add("hEta", "Pseudorapidty of muon pairs;;#counts", kTH1D, {axisEta}); registry.add("hRapidity", "Rapidty of muon pairs;;#counts", kTH1D, {axisRapidity}); registry.add("hPhi", "#varphi of muon pairs;;#counts", kTH1D, {axisPhi}); @@ -167,30 +189,50 @@ struct fwdMuonsUPC { registry.add("hEtaTrkNeg", "#eta of negative muons;;#counts", kTH1D, {axisEtaSingle}); registry.add("hPhiTrkPos", "#varphi of positive muons;;#counts", kTH1D, {axisPhiSingle}); registry.add("hPhiTrkNeg", "#varphi of negative muons;;#counts", kTH1D, {axisPhiSingle}); - registry.add("hTimeZNA", "ZNA Times;;#counts", kTH1D, {axisTimeZN}); - registry.add("hTimeZNC", "ZNC Times;;#counts", kTH1D, {axisTimeZN}); - registry.add("hEnergyZN", "ZNA vs ZNC energy", kTH2D, {axisEnergyZNA, axisEnergyZNC}); registry.add("hSameSign", "hSameSign;;#counts", kTH1D, {{6, -0.5, 5.5}}); registry.add("hPhiCharge", "#phi #it{charge}", kTH1D,{axisPhi}); registry.add("hPhiAverage", "#phi #it{average}", kTH1D,{axisPhi}); + // data + registry.add("hTimeZNA", "ZNA Times;;#counts", kTH1D, {axisTimeZN}); + registry.add("hTimeZNC", "ZNC Times;;#counts", kTH1D, {axisTimeZN}); + registry.add("hEnergyZN", "ZNA vs ZNC energy", kTH2D, {axisEnergyZNA, axisEnergyZNC}); + reg0n0n.add("hMass", "Ivariant mass of muon pairs - 0n0n;;#counts", kTH1D, {axisMass}); - reg0n0n.add("hPt", "Transverse momentum mass of muon pairs - 0n0n;;#counts", kTH1D, {axisPt}); + reg0n0n.add("hPt", "Transverse momentum of muon pairs - 0n0n;;#counts", kTH1D, {axisPt}); reg0n0n.add("hEta", "Pseudorapidty of muon pairs - 0n0n;;#counts", kTH1D, {axisEta}); reg0n0n.add("hRapidity", "Rapidty of muon pairs - 0n0n;;#counts", kTH1D, {axisRapidity}); - reg0n0n.add("hPtFit", "Transverse momentum mass of muon pairs - 0n0n;;#counts", kTH1D, {axisPtFit}); + reg0n0n.add("hPtFit", "Transverse momentum of muon pairs - 0n0n;;#counts", kTH1D, {axisPtFit}); regXn0n.add("hMass", "Ivariant mass of muon pairs - Xn0n;;#counts", kTH1D, {axisMass}); - regXn0n.add("hPt", "Transverse momentum mass of muon pairs - Xn0n;;#counts", kTH1D, {axisPt}); + regXn0n.add("hPt", "Transverse momentum of muon pairs - Xn0n;;#counts", kTH1D, {axisPt}); regXn0n.add("hEta", "Pseudorapidty of muon pairs - Xn0n;;#counts", kTH1D, {axisEta}); regXn0n.add("hRapidity", "Rapidty of muon pairs - Xn0n;;#counts", kTH1D, {axisRapidity}); - regXn0n.add("hPtFit", "Transverse momentum mass of muon pairs - Xn0n;;#counts", kTH1D, {axisPtFit}); + regXn0n.add("hPtFit", "Transverse momentum of muon pairs - Xn0n;;#counts", kTH1D, {axisPtFit}); regXnXn.add("hMass", "Ivariant mass of muon pairs - XnXn;;#counts", kTH1D, {axisMass}); - regXnXn.add("hPt", "Transverse momentum mass of muon pairs - XnXn;;#counts", kTH1D, {axisPt}); + regXnXn.add("hPt", "Transverse momentum of muon pairs - XnXn;;#counts", kTH1D, {axisPt}); regXnXn.add("hEta", "Pseudorapidty of muon pairs - XnXn;;#counts", kTH1D, {axisEta}); regXnXn.add("hRapidity", "Rapidty of muon pairs - XnXn;;#counts", kTH1D, {axisRapidity}); - regXnXn.add("hPtFit", "Transverse momentum mass of muon pairs - XnXn;;#counts", kTH1D, {axisPtFit}); + regXnXn.add("hPtFit", "Transverse momentum of muon pairs - XnXn;;#counts", kTH1D, {axisPtFit}); + + // gen MC + McGenRegistry.add("hMass", "Ivariant mass of muon pairs;;#counts", kTH1D, {axisMass}); + McGenRegistry.add("hPt", "Transverse momentum of muon pairs;;#counts", kTH1D, {axisPt}); + McGenRegistry.add("hEta", "Pseudorapidty of muon pairs;;#counts", kTH1D, {axisEta}); + McGenRegistry.add("hRapidity", "Rapidty of muon pairs;;#counts", kTH1D, {axisRapidity}); + McGenRegistry.add("hPhi", "#varphi of muon pairs;;#counts", kTH1D, {axisPhi}); + McGenRegistry.add("hPtTrkPos", "Pt of positive muons;;#counts", kTH1D, {axisPtSingle}); + McGenRegistry.add("hPtTrkNeg", "Pt of negative muons;;#counts", kTH1D, {axisPtSingle}); + McGenRegistry.add("hEtaTrkPos", "#eta of positive muons;;#counts", kTH1D, {axisEtaSingle}); + McGenRegistry.add("hEtaTrkNeg", "#eta of negative muons;;#counts", kTH1D, {axisEtaSingle}); + McGenRegistry.add("hPhiTrkPos", "#varphi of positive muons;;#counts", kTH1D, {axisPhiSingle}); + McGenRegistry.add("hPhiTrkNeg", "#varphi of negative muons;;#counts", kTH1D, {axisPhiSingle}); + McGenRegistry.add("hPhiCharge", "#phi #it{charge}", kTH1D,{axisPhi}); + McGenRegistry.add("hPhiAverage", "#phi #it{average}", kTH1D,{axisPhi}); + + //corr gen-reco + McCorrReg.add("hPtcorr", "gen pT vs reco pT", kTH2D, {axisPt, axisPt}); } // FUNCTIONS @@ -210,6 +252,38 @@ struct fwdMuonsUPC { } } + // template function that fills a map with the collision id of each udmccollision as key + // and a vector with the tracks + // map == (key, element) == (udMcCollisionId, vector of mc particles) + template + void collectMcCandIDs(std::unordered_map>& tracksPerCand, TTracks& tracks) + { + for (const auto& tr : tracks) { + int32_t candId = tr.udMcCollisionId(); + if (candId < 0) { + continue; + } + tracksPerCand[candId].push_back(tr.globalIndex()); + } + } + + template + void collectMcGenRecoCandIDs(std::unordered_map>& tracksPerCand, TTracks& tracks) + { + for (const auto& tr : tracks) { + if(tr.has_mcParticle()){ + int32_t candId = tr.udCollisionId(); + if (candId < 0) { + continue; + } + tracksPerCand[candId].push_back(tr.globalIndex()); + + auto mcPart = tr.mcParticle(); + tracksPerCand[candId].push_back(mcPart.globalIndex()); + } + } + } + // struct used to store the ZDC info in a map struct ZDCinfo { float timeA; @@ -269,6 +343,30 @@ struct fwdMuonsUPC { return true; } + //compute phi for azimuth anisotropy + void computePhiAnis(TLorentzVector p1, TLorentzVector p2, int sign1, float &phiAverage, float &phiCharge){ + + TLorentzVector tSum, tDiffAv, tDiffCh; + tSum = p1 + p2; + if(sign1>0){ + tDiffCh = p1 - p2; + if(gRandom->Rndm()>0.5) + tDiffAv = p1 - p2; + else tDiffAv = p2 - p1; + }else{ + tDiffCh = p2 - p1; + if(gRandom->Rndm()>0.5) + tDiffAv = p2 - p1; + else tDiffAv = p1 - p2; + } + + //average + phiAverage = tSum.DeltaPhi(tDiffAv); + //charge + phiCharge = tSum.DeltaPhi(tDiffCh); + } + + // function that processes the candidates: // it applies V0 selection, trk selection, and fills the histograms void processCand(CandidatesFwd::iterator const& cand, @@ -295,6 +393,12 @@ struct fwdMuonsUPC { if (!isMuonSelected(tr2)) return; + //select MCH-MID tracks only + //if(tr1.trackType() != 3) + // return; + //if(tr2.trackType() != 3) + // return; + // MCH-MID match selection int nMIDs = 0; if (tr1.chi2MatchMCHMID() > 0) @@ -316,9 +420,9 @@ struct fwdMuonsUPC { if(p.Pt() > highPt) return; // select rapidity - if (p.Rapidity() < yCutLow) + if (p.Rapidity() < lowRapidity) return; - if (p.Rapidity() > yCutUp) + if (p.Rapidity() > highRapidity) return; // select opposite charge events only @@ -377,9 +481,11 @@ struct fwdMuonsUPC { reg0n0n.fill(HIST("hEta"), p.Eta()); reg0n0n.fill(HIST("hRapidity"), p.Rapidity()); } else if (neutron_A ^ neutron_C) { // Xn0n + 0nXn - if (neutron_A) - znClass = 2; - else if (neutron_C) + if (neutron_A)// select opposite charge events only + if (cand.netCharge() != 0) { + registry.fill(HIST("hSameSign"), cand.numContrib()); + return; + } znClass = 3; regXn0n.fill(HIST("hMass"), p.M()); regXn0n.fill(HIST("hPt"), p.Pt()); @@ -417,22 +523,144 @@ struct fwdMuonsUPC { // store the event to save it into a tree if (tr1.sign() > 0) { - dimuSel(p.M(), p.Pt(), p.Rapidity(), p.Phi(), + dimuSel(p.M(), p.E(), p.Px(), p.Py(), p.Pz(), p.Pt(), p.Rapidity(), p.Phi(), phiAverage, phiCharge, - p1.Pt(), p1.PseudoRapidity(), p1.Phi(), - p2.Pt(), p2.PseudoRapidity(), p2.Phi(), - zdc.timeA, zdc.enA, zdc.timeC, zdc.enC, znClass); + p1.E(), p1.Px(), p1.Py(), p1.Pz(), p1.Pt(), p1.PseudoRapidity(), p1.Phi(), + p2.E(), p2.Px(), p2.Py(), p2.Pz(), p2.Pt(), p2.PseudoRapidity(), p2.Phi(), + zdc.timeA, zdc.enA, zdc.timeC, zdc.enC, znClass);TLorentzVector p1, p2; + p1.SetXYZM(tr1.px(), tr1.py(), tr1.pz(), mMu); + p2.SetXYZM(tr2.px(), tr2.py(), tr2.pz(), mMu); + TLorentzVector p = p1 + p2; } else { - dimuSel(p.M(), p.Pt(), p.Rapidity(), p.Phi(), + dimuSel(p.M(), p.E(), p.Px(), p.Py(), p.Pz(), p.Pt(), p.Rapidity(), p.Phi(), phiAverage, phiCharge, - p2.Pt(), p2.PseudoRapidity(), p2.Phi(), - p1.Pt(), p1.PseudoRapidity(), p1.Phi(), + p2.E(), p2.Px(), p2.Py(), p2.Pz(), p2.Pt(), p2.PseudoRapidity(), p2.Phi(), + p1.E(), p1.Px(), p1.Py(), p1.Pz(), p1.Pt(), p1.PseudoRapidity(), p1.Phi(), zdc.timeA, zdc.enA, zdc.timeC, zdc.enC, znClass); } + + } + + // function that processes the MC gen candidates: + // it applies some kinematics cut and fills the histograms + void processMcCand(aod::UDMcCollisions::iterator const& mcCand, + aod::UDMcParticles::iterator const& McPart1, aod::UDMcParticles::iterator const& McPart2){ + + //check that all pairs are mu+mu- + if(McPart1.pdgCode() + McPart2.pdgCode() != 0) LOGF(info,"PDG codes: %d | %d" ,McPart1.pdgCode(),McPart2.pdgCode()); + + // create Lorentz vectors + TLorentzVector p1, p2; + p1.SetXYZM(McPart1.px(), McPart1.py(), McPart1.pz(), mMu); + p2.SetXYZM(McPart2.px(), McPart2.py(), McPart2.pz(), mMu); + TLorentzVector p = p1 + p2; + + //cut on pair kinematics + //select mass + if(p.M() < lowMass) + return; + if(p.M() > highMass) + return; + //select pt + if(p.Pt() < lowPt) + return; + if(p.Pt() > highPt) + return; + // select rapidity + if (p.Rapidity() < lowRapidity) + return; + if (p.Rapidity() > highRapidity) + return; + + //compute phi for azimuth anisotropy + float phiAverage = 0; + float phiCharge = 0; + computePhiAnis(p1,p2,McPart1.pdgCode(),phiAverage, phiCharge); + + // fill the histos + McGenRegistry.fill(HIST("hPtTrkPos"), p1.Pt()); + McGenRegistry.fill(HIST("hPtTrkNeg"), p2.Pt()); + McGenRegistry.fill(HIST("hEtaTrkPos"), p1.Eta()); + McGenRegistry.fill(HIST("hEtaTrkNeg"), p2.Eta()); + McGenRegistry.fill(HIST("hPhiTrkPos"), p1.Phi()); + McGenRegistry.fill(HIST("hPhiTrkNeg"), p2.Phi()); + McGenRegistry.fill(HIST("hMass"), p.M()); + McGenRegistry.fill(HIST("hPt"), p.Pt()); + McGenRegistry.fill(HIST("hEta"), p.Eta()); + McGenRegistry.fill(HIST("hRapidity"), p.Rapidity()); + McGenRegistry.fill(HIST("hPhi"), p.Phi()); + McGenRegistry.fill(HIST("hPhiAverage"), phiAverage); + McGenRegistry.fill(HIST("hPhiCharge"), phiCharge); + + } + + void processMcCandGenReco(CandidatesFwd::iterator const& cand, + const ForwardTracks::iterator& tr1, const ForwardTracks::iterator& tr2, + aod::UDMcParticles::iterator const& McPart1, aod::UDMcParticles::iterator const& McPart2) + { + // V0 selection + const auto& ampsV0A = cand.amplitudesV0A(); + const auto& ampsRelBCsV0A = cand.ampRelBCsV0A(); + for (unsigned int i = 0; i < ampsV0A.size(); ++i) { + if (std::abs(ampsRelBCsV0A[i]) <= 1) { + if (ampsV0A[i] > 100.) + return; + } + } + + // track selection + TLorentzVector p1, p2; + p1.SetXYZM(tr1.px(), tr1.py(), tr1.pz(), mMu); + p2.SetXYZM(tr2.px(), tr2.py(), tr2.pz(), mMu); + TLorentzVector p = p1 + p2; + if (!isMuonSelected(tr1)) + return; + if (!isMuonSelected(tr2)) + return; + + // MCH-MID match selection + int nMIDs = 0; + if (tr1.chi2MatchMCHMID() > 0) + nMIDs++; + if (tr2.chi2MatchMCHMID() > 0) + nMIDs++; + if (nMIDs != 2) + return; + + // select opposite charge events only + if (cand.netCharge() != 0) { + registry.fill(HIST("hSameSign"), cand.numContrib()); + return; + } + + // mc particle + TLorentzVector p1Mc, p2Mc; + p1Mc.SetXYZM(McPart1.px(), McPart1.py(), McPart1.pz(), mMu); + p2Mc.SetXYZM(McPart2.px(), McPart2.py(), McPart2.pz(), mMu); + TLorentzVector pMc = p1Mc + p2Mc; + + //cut on pair kinematics + //select mass + if(p.M() < lowMass) + return; + if(p.M() > highMass) + return; + //select pt + if(p.Pt() < lowPt) + return; + if(p.Pt() > highPt) + return; + // select rapidity + if (p.Rapidity() < lowRapidity) + return; + if (p.Rapidity() > highRapidity) + return; + + McCorrReg.fill(HIST("hPtcorr"), p.Pt(),pMc.Pt()); } // PROCESS FUNCTION - void process(CandidatesFwd const& eventCandidates, + void processDataAndReco(CandidatesFwd const& eventCandidates, o2::aod::UDZdcsReduced& ZDCs, ForwardTracks const& fwdTracks) { @@ -469,7 +697,59 @@ struct fwdMuonsUPC { } } - PROCESS_SWITCH(fwdMuonsUPC, process, "", false); + PROCESS_SWITCH(fwdMuonsUPC, processDataAndReco, "", true); + + + // process MC Truth + void processGen(aod::UDMcCollisions const& mccollisions, aod::UDMcParticles const& McParts){ + + // map with the tracks + std::unordered_map> tracksPerCand; + collectMcCandIDs(tracksPerCand, McParts); + + // loop over the candidates + for (const auto& item : tracksPerCand) { + int32_t trId1 = item.second[0]; + int32_t trId2 = item.second[1]; + int32_t candID = item.first; + auto cand = mccollisions.iteratorAt(candID); + auto tr1 = McParts.iteratorAt(trId1); + auto tr2 = McParts.iteratorAt(trId2); + + processMcCand(cand, tr1, tr2); + } + } + PROCESS_SWITCH(fwdMuonsUPC, processGen, "", false); + +/* + void processGenAndReco(CandidatesFwd const& eventCandidates, + ForwardTracks const& fwdTracks, + aod::UDMcParticles const&) + { + // map with the tracks + std::unordered_map> tracksPerCand; + collectMcGenRecoCandIDs(tracksPerCand, fwdTracks); + + // loop over the candidates + for (const auto& item : tracksPerCand) { + int32_t trId1 = item.second[0]; + int32_t trId2 = item.second[1]; + + int32_t candID = item.first; + auto cand = eventCandidates.iteratorAt(candID); + auto tr1 = fwdTracks.iteratorAt(trId1); + auto tr2 = fwdTracks.iteratorAt(trId2); + + int32_t trMcId1 = item.second[2]; + int32_t trMcId2 = item.second[3]; + auto trMc1 = fwdTracks.iteratorAt(trMcId1); + auto trMc2 = fwdTracks.iteratorAt(trMcId2); + + processMcCandGenReco(cand, tr1, tr2, trMc1, trMc2); + } + } + PROCESS_SWITCH(fwdMuonsUPC, processGenAndReco, "", false); +*/ }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From f578c615d6bb32f04bdfddf587c98514a3ae548a Mon Sep 17 00:00:00 2001 From: ariffero Date: Wed, 20 Nov 2024 18:35:44 +0100 Subject: [PATCH 5/6] delete fwdMuonsMC and mumuForward and update CMakeList --- PWGUD/Tasks/CMakeLists.txt | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/PWGUD/Tasks/CMakeLists.txt b/PWGUD/Tasks/CMakeLists.txt index a1b9c1b90e5..a0f1eb17747 100644 --- a/PWGUD/Tasks/CMakeLists.txt +++ b/PWGUD/Tasks/CMakeLists.txt @@ -212,13 +212,3 @@ o2physics_add_dpl_workflow(event-by-event SOURCES eventByevent.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGPIDSelector COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(mumu-forward - SOURCES mumuForward.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGPIDSelector - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(fwd-muons-mc - SOURCES fwdMuonsMC.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGPIDSelector - COMPONENT_NAME Analysis) From ce094ee3f87b8223dbd438b092a7c8a35ae2fa93 Mon Sep 17 00:00:00 2001 From: ariffero Date: Wed, 27 Nov 2024 17:34:20 +0100 Subject: [PATCH 6/6] Modify fwdMuonsUPC.cxx - add process functions to analyze MC (pure MC and reconstructed MC) - correct typos - add info on the file - retrieve muon mass from PDG - add new variables to the tables to save as trees - compute phi for azimuth anisotropy --- PWGUD/Tasks/CMakeLists.txt | 2 +- PWGUD/Tasks/fwdMuonsUPC.cxx | 426 +++++++++++++++++++++++++----------- 2 files changed, 297 insertions(+), 131 deletions(-) diff --git a/PWGUD/Tasks/CMakeLists.txt b/PWGUD/Tasks/CMakeLists.txt index a0f1eb17747..b144b7becb8 100644 --- a/PWGUD/Tasks/CMakeLists.txt +++ b/PWGUD/Tasks/CMakeLists.txt @@ -211,4 +211,4 @@ o2physics_add_dpl_workflow(upc-rho-analysis o2physics_add_dpl_workflow(event-by-event SOURCES eventByevent.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGPIDSelector - COMPONENT_NAME Analysis) + COMPONENT_NAME Analysis) \ No newline at end of file diff --git a/PWGUD/Tasks/fwdMuonsUPC.cxx b/PWGUD/Tasks/fwdMuonsUPC.cxx index 770648a3149..eb2dfd073a1 100644 --- a/PWGUD/Tasks/fwdMuonsUPC.cxx +++ b/PWGUD/Tasks/fwdMuonsUPC.cxx @@ -9,6 +9,13 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +/// \file fwdMuonsUPC.cxx +/// \brief perform some selections on fwd events and saves the results + +/// executable name o2-analysis-ud-fwd-muon-upc + +/// \author Andrea Giovanni Riffero + #include "Framework/runDataProcessing.h" #include "Framework/AnalysisTask.h" #include "Framework/AnalysisDataModel.h" @@ -23,6 +30,8 @@ #include "TMath.h" #include "TRandom3.h" + +// table for saving tree with info on data namespace dimu { // dimuon @@ -69,14 +78,52 @@ DECLARE_SOA_TABLE(DiMu, "AOD", "DIMU", dimu::Tzna, dimu::Ezna, dimu::Tznc, dimu::Eznc, dimu::Nclass); } + +// for saving tree with info on reco MC +namespace recodimu +{ +// dimuon +DECLARE_SOA_COLUMN(Pt, pt, float); +DECLARE_SOA_COLUMN(Rap, rap, float); +DECLARE_SOA_COLUMN(Phi, phi, float); +DECLARE_SOA_COLUMN(PhiAv, phiAv, float); +DECLARE_SOA_COLUMN(PhiCh, phiCh, float); +// tracks positive (p) and negative (n) +DECLARE_SOA_COLUMN(Ptp, ptp, float); +DECLARE_SOA_COLUMN(Etap, etap, float); +DECLARE_SOA_COLUMN(Phip, phip, float); +DECLARE_SOA_COLUMN(Ptn, ptn, float); +DECLARE_SOA_COLUMN(Etan, etan, float); +DECLARE_SOA_COLUMN(Phin, phin, float); +//gen info dimuon +DECLARE_SOA_COLUMN(GenPt, gen_pt, float); +DECLARE_SOA_COLUMN(GenRap, gen_rap, float); +DECLARE_SOA_COLUMN(GenPhi, gen_phi, float); +//gen info trks +DECLARE_SOA_COLUMN(GenPtp, gen_ptp, float); +DECLARE_SOA_COLUMN(GenEtap, gen_etap, float); +DECLARE_SOA_COLUMN(GenPhip, gen_phip, float); +DECLARE_SOA_COLUMN(GenPtn, gen_ptn, float); +DECLARE_SOA_COLUMN(GenEtan, gen_etan, float); +DECLARE_SOA_COLUMN(GenPhin, gen_phin, float); +} // namespace recodimu + +namespace o2::aod +{ +DECLARE_SOA_TABLE(recoDiMu, "AOD", "RECODIMU", + recodimu::Pt, recodimu::Rap, recodimu::Phi, + recodimu::PhiAv, recodimu::PhiCh, + recodimu::Ptp, recodimu::Etap, recodimu::Phip, + recodimu::Ptn, recodimu::Etan, recodimu::Phin, + recodimu::GenPt, recodimu::GenRap, recodimu::GenPhi, + recodimu::GenPtp, recodimu::GenEtap, recodimu::GenPhip, + recodimu::GenPtn, recodimu::GenEtan, recodimu::GenPhin); +} + using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; -// defining constants -double mMu = 0.10566; // mass of muon - - // constants used in the track selection const float kRAbsMin = 17.6; const float kRAbsMid = 26.5; @@ -89,12 +136,15 @@ const float kPtMin = 0.; struct fwdMuonsUPC { + // a pdg object + TDatabasePDG* pdg = nullptr; + using CandidatesFwd = soa::Join; using ForwardTracks = soa::Join; - using CompleteCandidatesFwd = soa::Join; - using CompleteFwdTracks = soa::Join; + using CompleteFwdTracks = soa::Join; Produces dimuSel; + Produces dimuReco; // defining histograms using histogram registry: different histos for the different process functions HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; @@ -102,10 +152,7 @@ struct fwdMuonsUPC { HistogramRegistry regXn0n{"regXn0n", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; HistogramRegistry regXnXn{"regXnXn", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; HistogramRegistry McGenRegistry{"McGenRegistry", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; - HistogramRegistry McCorrReg{"McCorrReg", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; - - - + HistogramRegistry McRecoRegistry{"McRecoRegistry", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; // CONFIGURABLES // pT of muon pairs @@ -145,12 +192,12 @@ struct fwdMuonsUPC { Configurable lowEnZN{"lowEnZN", -50., "lower limit in ZN energy histo"}; Configurable highEnZN{"highEnZN", 250., "upper limit in ZN energy histo"}; - // configuarble rapidity cuts - //Configurable yCutLow{"yCutLow", -4, "Lower cut in pair rapidity"}; - //Configurable yCutUp{"yCutUp", -2.5, "Upper cut in pair rapidity"}; - void init(InitContext& context) + void init(InitContext&) { + // PDG + pdg = TDatabasePDG::Instance(); + // binning of pT axis fr fit std::vector ptFitBinning = { 0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10, @@ -174,13 +221,13 @@ struct fwdMuonsUPC { // histos // data and reco MC - registry.add("hMass", "Ivariant mass of muon pairs;;#counts", kTH1D, {axisMass}); + registry.add("hMass", "Invariant mass of muon pairs;;#counts", kTH1D, {axisMass}); registry.add("hPt", "Transverse momentum of muon pairs;;#counts", kTH1D, {axisPt}); registry.add("hPtFit", "Transverse momentum of muon pairs;;#counts", kTH1D, {axisPtFit}); registry.add("hEta", "Pseudorapidty of muon pairs;;#counts", kTH1D, {axisEta}); registry.add("hRapidity", "Rapidty of muon pairs;;#counts", kTH1D, {axisRapidity}); registry.add("hPhi", "#varphi of muon pairs;;#counts", kTH1D, {axisPhi}); - registry.add("hCharge", "Charge;#it{charge};;#counts", kTH1D, {{5, -2.5, 2.5}}); + registry.add("hCharge", "Charge;;;#counts", kTH1D, {{5, -2.5, 2.5}}); registry.add("hContrib", "hContrib;;#counts", kTH1D, {{6, -0.5, 5.5}}); registry.add("hEvSign", "Sum of the charges of all the tracks in each event;;#counts", kTH1D, {{5, -2.5, 2.5}}); registry.add("hPtTrkPos", "Pt of positive muons;;#counts", kTH1D, {axisPtSingle}); @@ -198,26 +245,26 @@ struct fwdMuonsUPC { registry.add("hTimeZNC", "ZNC Times;;#counts", kTH1D, {axisTimeZN}); registry.add("hEnergyZN", "ZNA vs ZNC energy", kTH2D, {axisEnergyZNA, axisEnergyZNC}); - reg0n0n.add("hMass", "Ivariant mass of muon pairs - 0n0n;;#counts", kTH1D, {axisMass}); + reg0n0n.add("hMass", "Invariant mass of muon pairs - 0n0n;;#counts", kTH1D, {axisMass}); reg0n0n.add("hPt", "Transverse momentum of muon pairs - 0n0n;;#counts", kTH1D, {axisPt}); reg0n0n.add("hEta", "Pseudorapidty of muon pairs - 0n0n;;#counts", kTH1D, {axisEta}); reg0n0n.add("hRapidity", "Rapidty of muon pairs - 0n0n;;#counts", kTH1D, {axisRapidity}); reg0n0n.add("hPtFit", "Transverse momentum of muon pairs - 0n0n;;#counts", kTH1D, {axisPtFit}); - regXn0n.add("hMass", "Ivariant mass of muon pairs - Xn0n;;#counts", kTH1D, {axisMass}); + regXn0n.add("hMass", "Invariant mass of muon pairs - Xn0n;;#counts", kTH1D, {axisMass}); regXn0n.add("hPt", "Transverse momentum of muon pairs - Xn0n;;#counts", kTH1D, {axisPt}); regXn0n.add("hEta", "Pseudorapidty of muon pairs - Xn0n;;#counts", kTH1D, {axisEta}); regXn0n.add("hRapidity", "Rapidty of muon pairs - Xn0n;;#counts", kTH1D, {axisRapidity}); regXn0n.add("hPtFit", "Transverse momentum of muon pairs - Xn0n;;#counts", kTH1D, {axisPtFit}); - regXnXn.add("hMass", "Ivariant mass of muon pairs - XnXn;;#counts", kTH1D, {axisMass}); + regXnXn.add("hMass", "Invariant mass of muon pairs - XnXn;;#counts", kTH1D, {axisMass}); regXnXn.add("hPt", "Transverse momentum of muon pairs - XnXn;;#counts", kTH1D, {axisPt}); regXnXn.add("hEta", "Pseudorapidty of muon pairs - XnXn;;#counts", kTH1D, {axisEta}); regXnXn.add("hRapidity", "Rapidty of muon pairs - XnXn;;#counts", kTH1D, {axisRapidity}); regXnXn.add("hPtFit", "Transverse momentum of muon pairs - XnXn;;#counts", kTH1D, {axisPtFit}); // gen MC - McGenRegistry.add("hMass", "Ivariant mass of muon pairs;;#counts", kTH1D, {axisMass}); + McGenRegistry.add("hMass", "Invariant mass of muon pairs;;#counts", kTH1D, {axisMass}); McGenRegistry.add("hPt", "Transverse momentum of muon pairs;;#counts", kTH1D, {axisPt}); McGenRegistry.add("hEta", "Pseudorapidty of muon pairs;;#counts", kTH1D, {axisEta}); McGenRegistry.add("hRapidity", "Rapidty of muon pairs;;#counts", kTH1D, {axisRapidity}); @@ -231,12 +278,45 @@ struct fwdMuonsUPC { McGenRegistry.add("hPhiCharge", "#phi #it{charge}", kTH1D,{axisPhi}); McGenRegistry.add("hPhiAverage", "#phi #it{average}", kTH1D,{axisPhi}); + // reco MC + McRecoRegistry.add("hMass", "Invariant mass of muon pairs;;#counts", kTH1D, {axisMass}); + McRecoRegistry.add("hPt", "Transverse momentum of muon pairs;;#counts", kTH1D, {axisPt}); + McRecoRegistry.add("hPtFit", "Transverse momentum of muon pairs;;#counts", kTH1D, {axisPtFit}); + McRecoRegistry.add("hEta", "Pseudorapidty of muon pairs;;#counts", kTH1D, {axisEta}); + McRecoRegistry.add("hRapidity", "Rapidty of muon pairs;;#counts", kTH1D, {axisRapidity}); + McRecoRegistry.add("hPhi", "#varphi of muon pairs;;#counts", kTH1D, {axisPhi}); + McRecoRegistry.add("hCharge", "Charge;;;#counts", kTH1D, {{5, -2.5, 2.5}}); + McRecoRegistry.add("hContrib", "hContrib;;#counts", kTH1D, {{6, -0.5, 5.5}}); + McRecoRegistry.add("hEvSign", "Sum of the charges of all the tracks in each event;;#counts", kTH1D, {{5, -2.5, 2.5}}); + McRecoRegistry.add("hPtTrkPos", "Pt of positive muons;;#counts", kTH1D, {axisPtSingle}); + McRecoRegistry.add("hPtTrkNeg", "Pt of negative muons;;#counts", kTH1D, {axisPtSingle}); + McRecoRegistry.add("hEtaTrkPos", "#eta of positive muons;;#counts", kTH1D, {axisEtaSingle}); + McRecoRegistry.add("hEtaTrkNeg", "#eta of negative muons;;#counts", kTH1D, {axisEtaSingle}); + McRecoRegistry.add("hPhiTrkPos", "#varphi of positive muons;;#counts", kTH1D, {axisPhiSingle}); + McRecoRegistry.add("hPhiTrkNeg", "#varphi of negative muons;;#counts", kTH1D, {axisPhiSingle}); + McRecoRegistry.add("hSameSign", "hSameSign;;#counts", kTH1D, {{6, -0.5, 5.5}}); + McRecoRegistry.add("hPhiCharge", "#phi #it{charge}", kTH1D,{axisPhi}); + McRecoRegistry.add("hPhiAverage", "#phi #it{average}", kTH1D,{axisPhi}); + //corr gen-reco - McCorrReg.add("hPtcorr", "gen pT vs reco pT", kTH2D, {axisPt, axisPt}); + McRecoRegistry.add("hPtcorr", "gen pT vs reco pT", kTH2D, {axisPt, axisPt}); + McRecoRegistry.add("hRapcorr", "gen rapidity vs reco rapidity", kTH2D, {axisRapidity, axisRapidity}); + McRecoRegistry.add("hPhicorr", "gen #phi vs reco #phi", kTH2D, {axisPhi, axisPhi}); } // FUNCTIONS + // retrieve particle mass (GeV/c^2) from TDatabasePDG + float particleMass(TDatabasePDG* pdg, int pid) + { + auto mass = 0.; + TParticlePDG* pdgparticle = pdg->GetParticle(pid); + if (pdgparticle != nullptr) { + mass = pdgparticle->Mass(); + } + return mass; + } + // template function that fills a map with the collision id of each udcollision as key // and a vector with the tracks // map == (key, element) == (udCollisionId, vector of trks) @@ -267,20 +347,27 @@ struct fwdMuonsUPC { } } + // template function that fills a map with the collision id of each udmccollision as key + // and a vector with the tracks + // map == (key, element) == (udCollisionId, vector(track1, mcPart1, udCollisionId1, track2, mcPart2, udCollisionId2)) template - void collectMcGenRecoCandIDs(std::unordered_map>& tracksPerCand, TTracks& tracks) + void collectRecoCandID(std::unordered_map>& tracksPerCand, TTracks& tracks) { for (const auto& tr : tracks) { - if(tr.has_mcParticle()){ - int32_t candId = tr.udCollisionId(); - if (candId < 0) { - continue; - } - tracksPerCand[candId].push_back(tr.globalIndex()); - - auto mcPart = tr.mcParticle(); - tracksPerCand[candId].push_back(mcPart.globalIndex()); + int32_t candId = tr.udCollisionId(); + if (candId < 0) + continue; + + if(!tr.has_udMcParticle()){ + //LOGF(info,"tr does not have mc part"); + continue; } + //retrieve mc particle from the reco track + auto mcPart = tr.udMcParticle(); + + tracksPerCand[candId].push_back(tr.globalIndex()); + tracksPerCand[candId].push_back(mcPart.globalIndex()); + tracksPerCand[candId].push_back(mcPart.udMcCollisionId()); } } @@ -322,11 +409,13 @@ struct fwdMuonsUPC { } // function to select muon tracks - bool isMuonSelected(const ForwardTracks::iterator& fwdTrack) + template + bool isMuonSelected(const TTracks& fwdTrack) { float rAbs = fwdTrack.rAtAbsorberEnd(); float pDca = fwdTrack.pDca(); TLorentzVector p; + auto mMu = particleMass(pdg, 13); p.SetXYZM(fwdTrack.px(), fwdTrack.py(), fwdTrack.pz(), mMu); float eta = p.Eta(); float pt = p.Pt(); @@ -342,8 +431,8 @@ struct fwdMuonsUPC { return false; return true; } - - //compute phi for azimuth anisotropy + + // function to compute phi for azimuth anisotropy void computePhiAnis(TLorentzVector p1, TLorentzVector p2, int sign1, float &phiAverage, float &phiCharge){ TLorentzVector tSum, tDiffAv, tDiffCh; @@ -360,17 +449,19 @@ struct fwdMuonsUPC { else tDiffAv = p1 - p2; } - //average + // average phiAverage = tSum.DeltaPhi(tDiffAv); - //charge + // charge phiCharge = tSum.DeltaPhi(tDiffCh); } // function that processes the candidates: - // it applies V0 selection, trk selection, and fills the histograms + // it applies V0 selection, trk selection, kine selection, and fills the histograms + // it also divides the data in neutron classes + // used for real data void processCand(CandidatesFwd::iterator const& cand, - const ForwardTracks::iterator& tr1, const ForwardTracks::iterator& tr2, + ForwardTracks::iterator const& tr1, ForwardTracks::iterator const& tr2, ZDCinfo& zdc) { // V0 selection @@ -382,23 +473,19 @@ struct fwdMuonsUPC { return; } } + + // select opposite charge events only + if (cand.netCharge() != 0) { + registry.fill(HIST("hSameSign"), cand.numContrib()); + return; + } // track selection - TLorentzVector p1, p2; - p1.SetXYZM(tr1.px(), tr1.py(), tr1.pz(), mMu); - p2.SetXYZM(tr2.px(), tr2.py(), tr2.pz(), mMu); - TLorentzVector p = p1 + p2; - if (!isMuonSelected(tr1)) + if (!isMuonSelected(*tr1)) return; - if (!isMuonSelected(tr2)) + if (!isMuonSelected(*tr2)) return; - //select MCH-MID tracks only - //if(tr1.trackType() != 3) - // return; - //if(tr2.trackType() != 3) - // return; - // MCH-MID match selection int nMIDs = 0; if (tr1.chi2MatchMCHMID() > 0) @@ -408,6 +495,13 @@ struct fwdMuonsUPC { if (nMIDs != 2) return; + //form Lorentz vectors + TLorentzVector p1, p2; + auto mMu = particleMass(pdg, 13); + p1.SetXYZM(tr1.px(), tr1.py(), tr1.pz(), mMu); + p2.SetXYZM(tr2.px(), tr2.py(), tr2.pz(), mMu); + TLorentzVector p = p1 + p2; + //cut on pair kinematics //select mass if(p.M() < lowMass) @@ -425,30 +519,11 @@ struct fwdMuonsUPC { if (p.Rapidity() > highRapidity) return; - // select opposite charge events only - if (cand.netCharge() != 0) { - registry.fill(HIST("hSameSign"), cand.numContrib()); - return; - } - - //compute phi for azimuth anisotropy - TLorentzVector tSum, tDiffAv, tDiffCh; - tSum = p1 + p2; - if(tr1.sign()>0){ - tDiffCh = p1 - p2; - if(gRandom->Rndm()>0.5) - tDiffAv = p1 - p2; - else tDiffAv = p2 - p1; - }else{ - tDiffCh = p2 - p1; - if(gRandom->Rndm()>0.5) - tDiffAv = p2 - p1; - else tDiffAv = p1 - p2; - } - //average - float phiAverage = tSum.DeltaPhi(tDiffAv); - float phiCharge = tSum.DeltaPhi(tDiffCh); - + // compute phi for azimuth anisotropy + float phiAverage = 0; + float phiCharge = 0; + computePhiAnis(p1,p2,tr1.sign(),phiAverage, phiCharge); + // zdc info if (TMath::Abs(zdc.timeA) < 10) registry.fill(HIST("hTimeZNA"), zdc.timeA); @@ -481,11 +556,9 @@ struct fwdMuonsUPC { reg0n0n.fill(HIST("hEta"), p.Eta()); reg0n0n.fill(HIST("hRapidity"), p.Rapidity()); } else if (neutron_A ^ neutron_C) { // Xn0n + 0nXn - if (neutron_A)// select opposite charge events only - if (cand.netCharge() != 0) { - registry.fill(HIST("hSameSign"), cand.numContrib()); - return; - } + if (neutron_A) + znClass = 2; + else if (neutron_C) znClass = 3; regXn0n.fill(HIST("hMass"), p.M()); regXn0n.fill(HIST("hPt"), p.Pt()); @@ -527,10 +600,7 @@ struct fwdMuonsUPC { phiAverage, phiCharge, p1.E(), p1.Px(), p1.Py(), p1.Pz(), p1.Pt(), p1.PseudoRapidity(), p1.Phi(), p2.E(), p2.Px(), p2.Py(), p2.Pz(), p2.Pt(), p2.PseudoRapidity(), p2.Phi(), - zdc.timeA, zdc.enA, zdc.timeC, zdc.enC, znClass);TLorentzVector p1, p2; - p1.SetXYZM(tr1.px(), tr1.py(), tr1.pz(), mMu); - p2.SetXYZM(tr2.px(), tr2.py(), tr2.pz(), mMu); - TLorentzVector p = p1 + p2; + zdc.timeA, zdc.enA, zdc.timeC, zdc.enC, znClass); } else { dimuSel(p.M(), p.E(), p.Px(), p.Py(), p.Pz(), p.Pt(), p.Rapidity(), p.Phi(), phiAverage, phiCharge, @@ -538,19 +608,19 @@ struct fwdMuonsUPC { p1.E(), p1.Px(), p1.Py(), p1.Pz(), p1.Pt(), p1.PseudoRapidity(), p1.Phi(), zdc.timeA, zdc.enA, zdc.timeC, zdc.enC, znClass); } - } // function that processes the MC gen candidates: // it applies some kinematics cut and fills the histograms - void processMcCand(aod::UDMcCollisions::iterator const& mcCand, - aod::UDMcParticles::iterator const& McPart1, aod::UDMcParticles::iterator const& McPart2){ + void processMcGenCand(aod::UDMcCollisions::iterator const& mcCand, + aod::UDMcParticles::iterator const& McPart1, aod::UDMcParticles::iterator const& McPart2){ //check that all pairs are mu+mu- if(McPart1.pdgCode() + McPart2.pdgCode() != 0) LOGF(info,"PDG codes: %d | %d" ,McPart1.pdgCode(),McPart2.pdgCode()); // create Lorentz vectors TLorentzVector p1, p2; + auto mMu = particleMass(pdg, 13); p1.SetXYZM(McPart1.px(), McPart1.py(), McPart1.pz(), mMu); p2.SetXYZM(McPart2.px(), McPart2.py(), McPart2.pz(), mMu); TLorentzVector p = p1 + p2; @@ -591,12 +661,13 @@ struct fwdMuonsUPC { McGenRegistry.fill(HIST("hPhi"), p.Phi()); McGenRegistry.fill(HIST("hPhiAverage"), phiAverage); McGenRegistry.fill(HIST("hPhiCharge"), phiCharge); - } - void processMcCandGenReco(CandidatesFwd::iterator const& cand, - const ForwardTracks::iterator& tr1, const ForwardTracks::iterator& tr2, - aod::UDMcParticles::iterator const& McPart1, aod::UDMcParticles::iterator const& McPart2) + // function that processes MC reco candidates + // it applies V0 selection, trk selection, kine selection, and fills the histograms + void processMcRecoCand(CandidatesFwd::iterator const& cand, + CompleteFwdTracks::iterator const& tr1, aod::UDMcParticles::iterator const& McPart1, + CompleteFwdTracks::iterator const& tr2, aod::UDMcParticles::iterator const& McPart2) { // V0 selection const auto& ampsV0A = cand.amplitudesV0A(); @@ -607,15 +678,17 @@ struct fwdMuonsUPC { return; } } + + // select opposite charge events only + if (cand.netCharge() != 0) { + registry.fill(HIST("hSameSign"), cand.numContrib()); + return; + } // track selection - TLorentzVector p1, p2; - p1.SetXYZM(tr1.px(), tr1.py(), tr1.pz(), mMu); - p2.SetXYZM(tr2.px(), tr2.py(), tr2.pz(), mMu); - TLorentzVector p = p1 + p2; - if (!isMuonSelected(tr1)) + if (!isMuonSelected(*tr1)) return; - if (!isMuonSelected(tr2)) + if (!isMuonSelected(*tr2)) return; // MCH-MID match selection @@ -627,19 +700,14 @@ struct fwdMuonsUPC { if (nMIDs != 2) return; - // select opposite charge events only - if (cand.netCharge() != 0) { - registry.fill(HIST("hSameSign"), cand.numContrib()); - return; - } - - // mc particle - TLorentzVector p1Mc, p2Mc; - p1Mc.SetXYZM(McPart1.px(), McPart1.py(), McPart1.pz(), mMu); - p2Mc.SetXYZM(McPart2.px(), McPart2.py(), McPart2.pz(), mMu); - TLorentzVector pMc = p1Mc + p2Mc; + // form Lorentz vectors + TLorentzVector p1, p2; + auto mMu = particleMass(pdg, 13); + p1.SetXYZM(tr1.px(), tr1.py(), tr1.pz(), mMu); + p2.SetXYZM(tr2.px(), tr2.py(), tr2.pz(), mMu); + TLorentzVector p = p1 + p2; - //cut on pair kinematics + //cut on pair kinematics (reco candidates) //select mass if(p.M() < lowMass) return; @@ -656,11 +724,95 @@ struct fwdMuonsUPC { if (p.Rapidity() > highRapidity) return; - McCorrReg.fill(HIST("hPtcorr"), p.Pt(),pMc.Pt()); + // compute phi for azimuth anisotropy + float phiAverage = 0; + float phiCharge = 0; + computePhiAnis(p1,p2,tr1.sign(),phiAverage, phiCharge); + + // gen particle + TLorentzVector p1Mc, p2Mc; + p1Mc.SetXYZM(McPart1.px(), McPart1.py(), McPart1.pz(), mMu); + p2Mc.SetXYZM(McPart2.px(), McPart2.py(), McPart2.pz(), mMu); + TLorentzVector pMc = p1Mc + p2Mc; + + //compute gen phi for azimuth anisotropy + float phiGenAverage = 0; + float phiGenCharge = 0; + computePhiAnis(p1,p2,McPart1.pdgCode(),phiGenAverage, phiGenCharge); + + // print info in case of problems + if(tr1.sign()*McPart1.pdgCode()>0 || tr2.sign()*McPart2.pdgCode()>0){ + LOGF(info, "Problem: "); + LOGF(info, "real: %d | %d",(int)tr1.sign(), (int)tr2.sign()); + LOGF(info, "mc : %i | %i",(int)McPart1.pdgCode(), (int)McPart2.pdgCode()); + LOGF(info, "contrib: %d", (int)cand.numContrib()); + } + + // fill the histos + // reco info + McRecoRegistry.fill(HIST("hContrib"), cand.numContrib()); + McRecoRegistry.fill(HIST("hPtTrkPos"), p1.Pt()); + McRecoRegistry.fill(HIST("hPtTrkNeg"), p2.Pt()); + McRecoRegistry.fill(HIST("hEtaTrkPos"), p1.Eta()); + McRecoRegistry.fill(HIST("hEtaTrkNeg"), p2.Eta()); + McRecoRegistry.fill(HIST("hPhiTrkPos"), p1.Phi()); + McRecoRegistry.fill(HIST("hPhiTrkNeg"), p2.Phi()); + McRecoRegistry.fill(HIST("hEvSign"), cand.netCharge()); + McRecoRegistry.fill(HIST("hMass"), p.M()); + McRecoRegistry.fill(HIST("hPt"), p.Pt()); + McRecoRegistry.fill(HIST("hPtFit"), p.Pt()); + McRecoRegistry.fill(HIST("hEta"), p.Eta()); + McRecoRegistry.fill(HIST("hRapidity"), p.Rapidity()); + McRecoRegistry.fill(HIST("hPhi"), p.Phi()); + McRecoRegistry.fill(HIST("hCharge"), tr1.sign()); + McRecoRegistry.fill(HIST("hCharge"), tr2.sign()); + McRecoRegistry.fill(HIST("hPhiAverage"), phiAverage); + McRecoRegistry.fill(HIST("hPhiCharge"), phiCharge); + + //gen info (of reco events) + McGenRegistry.fill(HIST("hPtTrkPos"), p1Mc.Pt()); + McGenRegistry.fill(HIST("hPtTrkNeg"), p2Mc.Pt()); + McGenRegistry.fill(HIST("hEtaTrkPos"), p1Mc.Eta()); + McGenRegistry.fill(HIST("hEtaTrkNeg"), p2Mc.Eta()); + McGenRegistry.fill(HIST("hPhiTrkPos"), p1Mc.Phi()); + McGenRegistry.fill(HIST("hPhiTrkNeg"), p2Mc.Phi()); + McGenRegistry.fill(HIST("hMass"), pMc.M()); + McGenRegistry.fill(HIST("hPt"), pMc.Pt()); + McGenRegistry.fill(HIST("hEta"), pMc.Eta()); + McGenRegistry.fill(HIST("hRapidity"), pMc.Rapidity()); + McGenRegistry.fill(HIST("hPhi"), pMc.Phi()); + McGenRegistry.fill(HIST("hPhiAverage"), phiGenAverage); + McGenRegistry.fill(HIST("hPhiCharge"), phiGenCharge); + + // reco-gen correlations + McRecoRegistry.fill(HIST("hPtcorr"), p.Pt(),pMc.Pt()); + McRecoRegistry.fill(HIST("hRapcorr"), p.Eta(),pMc.Eta()); + McRecoRegistry.fill(HIST("hPhicorr"), p.Phi(),pMc.Phi()); + + //store the event to save it into a tree + if(tr1.sign() > 0){ + dimuReco(p.Pt(), p.Rapidity(), p.Phi(), + phiAverage, phiCharge, + p1.Pt(), p1.PseudoRapidity(), p1.Phi(), + p2.Pt(), p2.PseudoRapidity(), p2.Phi(), + //gen info + pMc.Pt(), pMc.Rapidity(), pMc.Phi(), + p1Mc.Pt(), p1Mc.PseudoRapidity(), p1Mc.Phi(), + p2Mc.Pt(), p2Mc.PseudoRapidity(), p2Mc.Phi()); + } else { + dimuReco(p.Pt(), p.Rapidity(), p.Phi(), + phiAverage, phiCharge, + p2.Pt(), p2.PseudoRapidity(), p2.Phi(), + p1.Pt(), p1.PseudoRapidity(), p1.Phi(), + //gen info + pMc.Pt(), pMc.Rapidity(), pMc.Phi(), + p2Mc.Pt(), p2Mc.PseudoRapidity(), p2Mc.Phi(), + p1Mc.Pt(), p1Mc.PseudoRapidity(), p1Mc.Phi()); + } } // PROCESS FUNCTION - void processDataAndReco(CandidatesFwd const& eventCandidates, + void processData(CandidatesFwd const& eventCandidates, o2::aod::UDZdcsReduced& ZDCs, ForwardTracks const& fwdTracks) { @@ -697,11 +849,11 @@ struct fwdMuonsUPC { } } - PROCESS_SWITCH(fwdMuonsUPC, processDataAndReco, "", true); + PROCESS_SWITCH(fwdMuonsUPC, processData, "", true); // process MC Truth - void processGen(aod::UDMcCollisions const& mccollisions, aod::UDMcParticles const& McParts){ + void processMcGen(aod::UDMcCollisions const& mccollisions, aod::UDMcParticles const& McParts){ // map with the tracks std::unordered_map> tracksPerCand; @@ -716,40 +868,54 @@ struct fwdMuonsUPC { auto tr1 = McParts.iteratorAt(trId1); auto tr2 = McParts.iteratorAt(trId2); - processMcCand(cand, tr1, tr2); + processMcGenCand(cand, tr1, tr2); } } - PROCESS_SWITCH(fwdMuonsUPC, processGen, "", false); + PROCESS_SWITCH(fwdMuonsUPC, processMcGen, "", false); -/* - void processGenAndReco(CandidatesFwd const& eventCandidates, - ForwardTracks const& fwdTracks, - aod::UDMcParticles const&) + + // process reco MC (gen info included) + void processMcReco(CandidatesFwd const& eventCandidates, + CompleteFwdTracks const& fwdTracks, + aod::UDMcCollisions const& mcCandidates, + aod::UDMcParticles const& McParts) { - // map with the tracks - std::unordered_map> tracksPerCand; - collectMcGenRecoCandIDs(tracksPerCand, fwdTracks); + std::unordered_map> tracksPerCandAll; + collectRecoCandID(tracksPerCandAll, fwdTracks); // loop over the candidates - for (const auto& item : tracksPerCand) { + for (const auto& item : tracksPerCandAll) { + if(item.second.size()!=6){ + //LOGF(info, "error: reco track(s) not gen"); + continue; + }; int32_t trId1 = item.second[0]; - int32_t trId2 = item.second[1]; + int32_t trId2 = item.second[3]; //[2] int32_t candID = item.first; auto cand = eventCandidates.iteratorAt(candID); auto tr1 = fwdTracks.iteratorAt(trId1); auto tr2 = fwdTracks.iteratorAt(trId2); - int32_t trMcId1 = item.second[2]; - int32_t trMcId2 = item.second[3]; - auto trMc1 = fwdTracks.iteratorAt(trMcId1); - auto trMc2 = fwdTracks.iteratorAt(trMcId2); + auto trMcId1 = item.second[1]; + auto trMcId2 = item.second[4]; + auto trMc1 = McParts.iteratorAt(trMcId1); + auto trMc2 = McParts.iteratorAt(trMcId2); - processMcCandGenReco(cand, tr1, tr2, trMc1, trMc2); + auto mcCandID1 = mcCandidates.iteratorAt(item.second[2]); + auto mcCandID2 = mcCandidates.iteratorAt(item.second[5]); + + if(mcCandID1 != mcCandID2){ + //LOGF(info, "mc tracks belong to different collisions"); + } + + //auto mcTr1 = McParts.iteratorAt(tr1.udMcParticleId()); + //auto mcTr2 = McParts.iteratorAt(tr2.udMcParticleId()); + processMcRecoCand(cand, tr1, trMc1, tr2, trMc2); } } - PROCESS_SWITCH(fwdMuonsUPC, processGenAndReco, "", false); -*/ + PROCESS_SWITCH(fwdMuonsUPC, processMcReco, "", false); + }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)