diff --git a/PWGUD/Tasks/CMakeLists.txt b/PWGUD/Tasks/CMakeLists.txt index b30af171280..963307998fb 100644 --- a/PWGUD/Tasks/CMakeLists.txt +++ b/PWGUD/Tasks/CMakeLists.txt @@ -200,7 +200,7 @@ o2physics_add_dpl_workflow(upc-pion-analysis COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(fwd-muons-upc - SOURCES FwdMuonsUPC.cxx + SOURCES fwdMuonsUPC.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 similarity index 58% rename from PWGUD/Tasks/FwdMuonsUPC.cxx rename to PWGUD/Tasks/fwdMuonsUPC.cxx index b3351965727..0f1cc7eb6eb 100644 --- a/PWGUD/Tasks/FwdMuonsUPC.cxx +++ b/PWGUD/Tasks/fwdMuonsUPC.cxx @@ -9,7 +9,7 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// \file FwdMuonsUPC.cxx +/// \file fwdMuonsUPC.cxx /// \brief perform some selections on fwd events and saves the results /// executable name o2-analysis-ud-fwd-muon-upc @@ -18,24 +18,17 @@ #include "PWGUD/DataModel/UDTables.h" -#include +#include +#include #include -#include #include -#include -#include -#include -#include #include -#include #include -#include +#include +#include #include -#include -#include -#include #include #include @@ -183,12 +176,12 @@ const int kReqMatchMIDTracks = 2; const int kReqMatchMFTTracks = 2; const int kMaxChi2MFTMatch = 30; const float kMaxZDCTime = 2.; -const float kMaxZDCTimeHisto = 10.; -const int kMuonPDG = 13; +const int k2Tracks = 2; +const int k4Tracks = 4; -struct FwdMuonsUPC { +struct fwdMuonsUPC { - // a pdg object + // a PDG object Service pdg; using CandidatesFwd = soa::Join; @@ -201,11 +194,6 @@ struct FwdMuonsUPC { // 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 mcRecoRegistry{"mcRecoRegistry", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; // CONFIGURABLES static constexpr double Pi = o2::constants::math::PI; @@ -217,34 +205,10 @@ struct FwdMuonsUPC { Configurable nBinsMass{"nBinsMass", 500, "N bins in mass histo"}; Configurable lowMass{"lowMass", 0., "lower limit in mass histo"}; Configurable highMass{"highMass", 10., "upper limit in mass histo"}; - // eta of muon pairs - Configurable nBinsEta{"nBinsEta", 600, "N bins in eta histo"}; - Configurable lowEta{"lowEta", -10., "lower limit in eta histo"}; - Configurable highEta{"highEta", -2., "upper limit in eta histo"}; // rapidity of muon pairs Configurable nBinsRapidity{"nBinsRapidity", 250, "N bins in rapidity histo"}; Configurable lowRapidity{"lowRapidity", -4.5, "lower limit in rapidity histo"}; Configurable highRapidity{"highRapidity", -2., "upper limit in rapidity histo"}; - // phi of muon pairs - Configurable nBinsPhi{"nBinsPhi", 600, "N bins in phi histo"}; - Configurable lowPhi{"lowPhi", -Pi, "lower limit in phi histo"}; - Configurable highPhi{"highPhi", Pi, "upper limit in phi histo"}; - // pT of single muons - Configurable nBinsPtSingle{"nBinsPtSingle", 500, "N bins in pT histo single muon"}; - Configurable lowPtSingle{"lowPtSingle", 0., "lower limit in pT histo single muon"}; - Configurable highPtSingle{"highPtSingle", 2., "upper limit in pT histo single muon"}; - // eta of single muons - Configurable nBinsEtaSingle{"nBinsEtaSingle", 250, "N bins in eta histo single muon"}; - Configurable lowEtaSingle{"lowEtaSingle", -4.5, "lower limit in eta histo single muon"}; - Configurable highEtaSingle{"highEtaSingle", -2., "upper limit in eta histo single muon"}; - // phi of single muons - Configurable nBinsPhiSingle{"nBinsPhiSingle", 600, "N bins in phi histo single muon"}; - Configurable lowPhiSingle{"lowPhiSingle", -Pi, "lower limit in phi histo single muon"}; - Configurable highPhiSingle{"highPhiSingle", Pi, "upper limit in phi histo single muon"}; - // ZDC - Configurable nBinsZDCen{"nBinsZDCen", 200, "N bins in ZN energy"}; - Configurable lowEnZN{"lowEnZN", -50., "lower limit in ZN energy histo"}; - Configurable highEnZN{"highEnZN", 250., "upper limit in ZN energy histo"}; // my track type // 0 = MCH-MID-MFT // 1 = MCH-MID @@ -271,114 +235,18 @@ struct FwdMuonsUPC { // axis const AxisSpec axisPt{nBinsPt, lowPt, highPt, "#it{p}_{T} GeV/#it{c}"}; const AxisSpec axisPtFit = {ptFitBinning, "#it{p}_{T} (GeV/c)"}; - const AxisSpec axisPtFit2 = {ptFitBinningHalfWidth, "#it{p}_{T} (GeV/c)"}; const AxisSpec axisMass{nBinsMass, lowMass, highMass, "m_{#mu#mu} GeV/#it{c}^{2}"}; - const AxisSpec axisEta{nBinsEta, lowEta, highEta, "#eta"}; const AxisSpec axisRapidity{nBinsRapidity, lowRapidity, highRapidity, "Rapidity"}; - const AxisSpec axisPhi{nBinsPhi, lowPhi, highPhi, "#varphi"}; - const AxisSpec axisPtSingle{nBinsPtSingle, lowPtSingle, highPtSingle, "#it{p}_{T}_{ trk} GeV/#it{c}"}; - const AxisSpec axisTimeZN{200, -10, 10, "ZDC time (ns)"}; - const AxisSpec axisEnergyZNA{nBinsZDCen, lowEnZN, highEnZN, "ZNA energy (TeV)"}; - const AxisSpec axisEnergyZNC{nBinsZDCen, lowEnZN, highEnZN, "ZNC energy (TeV)"}; - const AxisSpec axisEtaSingle{nBinsEtaSingle, lowEtaSingle, highEtaSingle, "#eta_{trk}"}; - const AxisSpec axisPhiSingle{nBinsPhiSingle, lowPhiSingle, highPhiSingle, "#varphi_{trk}"}; // histos - // data and reco MC 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("hPtFit2", "Transverse momentum of muon pairs;;#counts", kTH1D, {axisPtFit2}); - 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;;;#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}); - registry.add("hPtTrkNeg", "Pt of negative muons;;#counts", kTH1D, {axisPtSingle}); - registry.add("hEtaTrkPos", "#eta of positive muons;;#counts", kTH1D, {axisEtaSingle}); - 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("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", "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", "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", "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", "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}); - 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}); - - // 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("hPtFit2", "Transverse momentum of muon pairs;;#counts", kTH1D, {axisPtFit2}); - 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 - 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(int pid) - { - auto mass = pdg->Mass(pid); - 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) @@ -405,6 +273,9 @@ struct FwdMuonsUPC { if (candId < 0) { continue; } + if (std::abs(tr.pdgCode()) != pdg->GetParticle("mu-")->PdgCode()) { + continue; + } tracksPerCand[candId].push_back(tr.globalIndex()); } } @@ -475,9 +346,8 @@ struct FwdMuonsUPC { { float rAbs = fwdTrack.rAtAbsorberEnd(); float pDca = fwdTrack.pDca(); - TLorentzVector p; - auto mMu = particleMass(kMuonPDG); - p.SetXYZM(fwdTrack.px(), fwdTrack.py(), fwdTrack.pz(), mMu); + ROOT::Math::PxPyPzMVector p{fwdTrack.px(), fwdTrack.py(), fwdTrack.pz(), o2::constants::physics::MassMuon}; + float eta = p.Eta(); float pt = p.Pt(); float pDcaMax = rAbs < kRAbsMid ? kPDca1 : kPDca2; @@ -494,9 +364,9 @@ struct FwdMuonsUPC { } // function to compute phi for azimuth anisotropy - void computePhiAnis(TLorentzVector p1, TLorentzVector p2, int sign1, float& phiAverage, float& phiCharge) + void computePhiAnis(ROOT::Math::PxPyPzMVector p1, ROOT::Math::PxPyPzMVector p2, int sign1, float& phiAverage, float& phiCharge) { - TLorentzVector tSum, tDiffAv, tDiffCh; + ROOT::Math::PxPyPzMVector tSum, tDiffAv, tDiffCh; tSum = p1 + p2; float halfUnity = 0.5; if (sign1 > 0) { @@ -514,9 +384,9 @@ struct FwdMuonsUPC { } // average - phiAverage = tSum.DeltaPhi(tDiffAv); + phiAverage = ROOT::Math::VectorUtil::DeltaPhi(tSum, tDiffAv); // charge - phiCharge = tSum.DeltaPhi(tDiffCh); + phiCharge = ROOT::Math::VectorUtil::DeltaPhi(tSum, tDiffCh); } // function that processes the candidates: @@ -537,9 +407,14 @@ struct FwdMuonsUPC { } } + // select events with exactly 2 forward tracks + if (cand.numContrib() != k2Tracks) { + return; + } + // select opposite charge events only if (cand.netCharge() != 0) { - registry.fill(HIST("hSameSign"), cand.numContrib()); + // registry.fill(HIST("hSameSign"), cand.numContrib()); return; } @@ -574,11 +449,9 @@ struct FwdMuonsUPC { return; // form Lorentz vectors - TLorentzVector p1, p2; - auto mMu = particleMass(kMuonPDG); - p1.SetXYZM(tr1.px(), tr1.py(), tr1.pz(), mMu); - p2.SetXYZM(tr2.px(), tr2.py(), tr2.pz(), mMu); - TLorentzVector p = p1 + p2; + ROOT::Math::PxPyPzMVector p1{tr1.px(), tr1.py(), tr1.pz(), o2::constants::physics::MassMuon}; + ROOT::Math::PxPyPzMVector p2{tr2.px(), tr2.py(), tr2.pz(), o2::constants::physics::MassMuon}; + ROOT::Math::PxPyPzMVector p = p1 + p2; // cut on pair kinematics // select mass @@ -602,13 +475,6 @@ struct FwdMuonsUPC { float phiCharge = 0; computePhiAnis(p1, p2, tr1.sign(), phiAverage, phiCharge); - // zdc info - if (std::abs(zdc.timeA) < kMaxZDCTimeHisto) - registry.fill(HIST("hTimeZNA"), zdc.timeA); - if (std::abs(zdc.timeC) < kMaxZDCTimeHisto) - registry.fill(HIST("hTimeZNC"), zdc.timeC); - registry.fill(HIST("hEnergyZN"), zdc.enA, zdc.enC); - // divide the events in neutron classes bool neutronA = false; bool neutronC = false; @@ -628,65 +494,35 @@ struct FwdMuonsUPC { // 0n0n if (neutronC == false && neutronA == false) { znClass = 1; - reg0n0n.fill(HIST("hMass"), p.M()); - reg0n0n.fill(HIST("hPt"), p.Pt()); - reg0n0n.fill(HIST("hPtFit"), p.Pt()); - reg0n0n.fill(HIST("hEta"), p.Eta()); - reg0n0n.fill(HIST("hRapidity"), p.Rapidity()); } else if (neutronA ^ neutronC) { // Xn0n + 0nXn if (neutronA) znClass = 2; else if (neutronC) 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 (neutronA && neutronC) { // XnXn znClass = 4; - regXnXn.fill(HIST("hMass"), p.M()); - regXnXn.fill(HIST("hPt"), p.Pt()); - regXnXn.fill(HIST("hPtFit"), p.Pt()); - regXnXn.fill(HIST("hEta"), p.Eta()); - regXnXn.fill(HIST("hRapidity"), p.Rapidity()); } - // fill the histos without looking at neutron emission - registry.fill(HIST("hContrib"), cand.numContrib()); - registry.fill(HIST("hPtTrkPos"), p1.Pt()); - registry.fill(HIST("hPtTrkNeg"), p2.Pt()); - registry.fill(HIST("hEtaTrkPos"), p1.Eta()); - registry.fill(HIST("hEtaTrkNeg"), p2.Eta()); - registry.fill(HIST("hPhiTrkPos"), p1.Phi()); - registry.fill(HIST("hPhiTrkNeg"), p2.Phi()); - registry.fill(HIST("hEvSign"), cand.netCharge()); + // fill the histos registry.fill(HIST("hMass"), p.M()); registry.fill(HIST("hPt"), p.Pt()); registry.fill(HIST("hPtFit"), p.Pt()); - registry.fill(HIST("hPtFit2"), p.Pt()); - registry.fill(HIST("hEta"), p.Eta()); registry.fill(HIST("hRapidity"), p.Rapidity()); - 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(cand.runNumber(), p.M(), p.E(), p.Px(), p.Py(), p.Pz(), p.Pt(), p.Rapidity(), p.Phi(), phiAverage, phiCharge, - p1.E(), p1.Px(), p1.Py(), p1.Pz(), p1.Pt(), p1.PseudoRapidity(), p1.Phi(), static_cast(myTrackType), - p2.E(), p2.Px(), p2.Py(), p2.Pz(), p2.Pt(), p2.PseudoRapidity(), p2.Phi(), static_cast(myTrackType), + p1.E(), p1.Px(), p1.Py(), p1.Pz(), p1.Pt(), p1.Eta(), p1.Phi(), static_cast(myTrackType), + p2.E(), p2.Px(), p2.Py(), p2.Pz(), p2.Pt(), p2.Eta(), p2.Phi(), static_cast(myTrackType), zdc.timeA, zdc.enA, zdc.timeC, zdc.enC, znClass); } else { dimuSel(cand.runNumber(), p.M(), p.E(), p.Px(), p.Py(), p.Pz(), p.Pt(), p.Rapidity(), p.Phi(), phiAverage, phiCharge, - p2.E(), p2.Px(), p2.Py(), p2.Pz(), p2.Pt(), p2.PseudoRapidity(), p2.Phi(), static_cast(myTrackType), - p1.E(), p1.Px(), p1.Py(), p1.Pz(), p1.Pt(), p1.PseudoRapidity(), p1.Phi(), static_cast(myTrackType), + p2.E(), p2.Px(), p2.Py(), p2.Pz(), p2.Pt(), p2.Eta(), p2.Phi(), static_cast(myTrackType), + p1.E(), p1.Px(), p1.Py(), p1.Pz(), p1.Pt(), p1.Eta(), p1.Phi(), static_cast(myTrackType), zdc.timeA, zdc.enA, zdc.timeC, zdc.enC, znClass); } } @@ -698,15 +534,18 @@ struct FwdMuonsUPC { { // check that all pairs are mu+mu- - if (std::abs(McPart1.pdgCode()) != kMuonPDG && std::abs(McPart2.pdgCode()) != kMuonPDG) + if (std::abs(McPart1.pdgCode()) != pdg->GetParticle("mu-")->PdgCode() || std::abs(McPart2.pdgCode()) != pdg->GetParticle("mu-")->PdgCode()) { LOGF(debug, "PDG codes: %d | %d", McPart1.pdgCode(), McPart2.pdgCode()); + return; + } + if (McPart1.pdgCode() + McPart2.pdgCode() != 0) { + return; + } // create Lorentz vectors - TLorentzVector p1, p2; - auto mMu = particleMass(kMuonPDG); - p1.SetXYZM(McPart1.px(), McPart1.py(), McPart1.pz(), mMu); - p2.SetXYZM(McPart2.px(), McPart2.py(), McPart2.pz(), mMu); - TLorentzVector p = p1 + p2; + ROOT::Math::PxPyPzMVector p1{McPart1.px(), McPart1.py(), McPart1.pz(), o2::constants::physics::MassMuon}; + ROOT::Math::PxPyPzMVector p2{McPart2.px(), McPart2.py(), McPart2.pz(), o2::constants::physics::MassMuon}; + ROOT::Math::PxPyPzMVector p = p1 + p2; // cut on pair kinematics // select mass @@ -731,31 +570,22 @@ struct FwdMuonsUPC { 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); + registry.fill(HIST("hMass"), p.M()); + registry.fill(HIST("hPt"), p.Pt()); + registry.fill(HIST("hPtFit"), p.Pt()); + registry.fill(HIST("hRapidity"), p.Rapidity()); // store the event to save it into a tree if (McPart1.pdgCode() < 0) { dimuGen(p.M(), p.Pt(), p.Rapidity(), p.Phi(), phiAverage, phiCharge, - p1.Pt(), p1.PseudoRapidity(), p1.Phi(), - p2.Pt(), p2.PseudoRapidity(), p2.Phi()); + p1.Pt(), p1.Eta(), p1.Phi(), + p2.Pt(), p2.Eta(), p2.Phi()); } else { dimuGen(p.M(), p.Pt(), p.Rapidity(), p.Phi(), phiAverage, phiCharge, - p2.Pt(), p2.PseudoRapidity(), p2.Phi(), - p1.Pt(), p1.PseudoRapidity(), p1.Phi()); + p2.Pt(), p2.Eta(), p2.Phi(), + p1.Pt(), p1.Eta(), p1.Phi()); } } @@ -767,7 +597,7 @@ struct FwdMuonsUPC { { // check that all pairs are mu+mu- - if (std::abs(McPart1.pdgCode()) != kMuonPDG && std::abs(McPart2.pdgCode()) != kMuonPDG) + if (std::abs(McPart1.pdgCode()) != pdg->GetParticle("mu-")->PdgCode() || std::abs(McPart2.pdgCode()) != pdg->GetParticle("mu-")->PdgCode()) LOGF(debug, "PDG codes: %d | %d", McPart1.pdgCode(), McPart2.pdgCode()); // V0 selection @@ -782,7 +612,7 @@ struct FwdMuonsUPC { // select opposite charge events only if (cand.netCharge() != 0) { - registry.fill(HIST("hSameSign"), cand.numContrib()); + // registry.fill(HIST("hSameSign"), cand.numContrib()); return; } @@ -817,11 +647,9 @@ struct FwdMuonsUPC { return; // form Lorentz vectors - TLorentzVector p1, p2; - auto mMu = particleMass(kMuonPDG); - p1.SetXYZM(tr1.px(), tr1.py(), tr1.pz(), mMu); - p2.SetXYZM(tr2.px(), tr2.py(), tr2.pz(), mMu); - TLorentzVector p = p1 + p2; + ROOT::Math::PxPyPzMVector p1{tr1.px(), tr1.py(), tr1.pz(), o2::constants::physics::MassMuon}; + ROOT::Math::PxPyPzMVector p2{tr2.px(), tr2.py(), tr2.pz(), o2::constants::physics::MassMuon}; + ROOT::Math::PxPyPzMVector p = p1 + p2; // cut on pair kinematics (reco candidates) // select mass @@ -845,16 +673,14 @@ struct FwdMuonsUPC { 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; + ROOT::Math::PxPyPzMVector p1Mc{McPart1.px(), McPart1.py(), McPart1.pz(), o2::constants::physics::MassMuon}; + ROOT::Math::PxPyPzMVector p2Mc{McPart2.px(), McPart2.py(), McPart2.pz(), o2::constants::physics::MassMuon}; + ROOT::Math::PxPyPzMVector pMc = p2Mc + p2Mc; // compute gen phi for azimuth anisotropy float phiGenAverage = 0; float phiGenCharge = 0; - computePhiAnis(p1, p2, -McPart1.pdgCode(), phiGenAverage, phiGenCharge); + computePhiAnis(p1Mc, p2Mc, -McPart1.pdgCode(), phiGenAverage, phiGenCharge); // print info in case of problems if (tr1.sign() * McPart1.pdgCode() > 0 || tr2.sign() * McPart2.pdgCode() > 0) { @@ -865,68 +691,32 @@ struct FwdMuonsUPC { } // 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("hPtFit2"), 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.Rapidity(), pMc.Rapidity()); - mcRecoRegistry.fill(HIST("hPhicorr"), p.Phi(), pMc.Phi()); + registry.fill(HIST("hMass"), p.M()); + registry.fill(HIST("hPt"), p.Pt()); + registry.fill(HIST("hPtFit"), p.Pt()); + registry.fill(HIST("hRapidity"), p.Rapidity()); // store the event to save it into a tree if (tr1.sign() > 0) { dimuReco(cand.runNumber(), p.M(), p.Pt(), p.Rapidity(), p.Phi(), phiAverage, phiCharge, - p1.Pt(), p1.PseudoRapidity(), p1.Phi(), static_cast(myTrackType), - p2.Pt(), p2.PseudoRapidity(), p2.Phi(), static_cast(myTrackType), + p1.Pt(), p1.Eta(), p1.Phi(), static_cast(myTrackType), + p2.Pt(), p2.Eta(), p2.Phi(), static_cast(myTrackType), // gen info pMc.Pt(), pMc.Rapidity(), pMc.Phi(), - p1Mc.Pt(), p1Mc.PseudoRapidity(), p1Mc.Phi(), - p2Mc.Pt(), p2Mc.PseudoRapidity(), p2Mc.Phi()); + p1Mc.Pt(), p1Mc.Eta(), p1Mc.Phi(), + p2Mc.Pt(), p2Mc.Eta(), p2Mc.Phi()); } else { dimuReco(cand.runNumber(), p.M(), p.Pt(), p.Rapidity(), p.Phi(), phiAverage, phiCharge, - p2.Pt(), p2.PseudoRapidity(), p2.Phi(), static_cast(myTrackType), - p1.Pt(), p1.PseudoRapidity(), p1.Phi(), static_cast(myTrackType), + p2.Pt(), p2.Eta(), p2.Phi(), static_cast(myTrackType), + p1.Pt(), p1.Eta(), p1.Phi(), static_cast(myTrackType), // gen info pMc.Pt(), pMc.Rapidity(), pMc.Phi(), - p2Mc.Pt(), p2Mc.PseudoRapidity(), p2Mc.Phi(), - p1Mc.Pt(), p1Mc.PseudoRapidity(), p1Mc.Phi()); + p2Mc.Pt(), p2Mc.Eta(), p2Mc.Phi(), + p1Mc.Pt(), p1Mc.Eta(), p1Mc.Phi()); } } @@ -946,6 +736,10 @@ struct FwdMuonsUPC { // loop over the candidates for (const auto& item : tracksPerCand) { + if (item.second.size() != k2Tracks) { + LOGF(debug, "number track = %d", item.second.size()); + continue; + } int32_t trId1 = item.second[0]; int32_t trId2 = item.second[1]; int32_t candID = item.first; @@ -968,18 +762,28 @@ struct FwdMuonsUPC { } } - PROCESS_SWITCH(FwdMuonsUPC, processData, "", true); + PROCESS_SWITCH(fwdMuonsUPC, processData, "", true); // process MC Truth void processMcGen(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) { + if (item.second.size() != k2Tracks) { + LOGF(debug, "mc parts = %d", item.second.size()); + for (const auto& id : item.second) { + auto p = McParts.iteratorAt(id); + LOGF(debug, + " part %d: pdg=%d status=%d has_mothers=%d has_daughters=%d", + id, p.pdgCode(), p.statusCode(), + p.has_mothers(), p.has_daughters()); + } + continue; + } int32_t trId1 = item.second[0]; int32_t trId2 = item.second[1]; int32_t candID = item.first; @@ -990,12 +794,12 @@ struct FwdMuonsUPC { processMcGenCand(cand, tr1, tr2); } } - PROCESS_SWITCH(FwdMuonsUPC, processMcGen, "", false); + PROCESS_SWITCH(fwdMuonsUPC, processMcGen, "", false); // process reco MC (gen info included) void processMcReco(CandidatesFwd const& eventCandidates, CompleteFwdTracks const& fwdTracks, - aod::UDMcCollisions const&, + aod::UDMcCollisions const& /*mccollisions*/, aod::UDMcParticles const& McParts) { std::unordered_map> tracksPerCandAll; @@ -1003,7 +807,7 @@ struct FwdMuonsUPC { // loop over the candidates for (const auto& item : tracksPerCandAll) { - if (item.second.size() != 4) { + if (item.second.size() != k4Tracks) { LOGF(debug, "number track (reco + gen) = %d", item.second.size()); continue; } @@ -1030,12 +834,12 @@ struct FwdMuonsUPC { processMcRecoCand(cand, tr1, trMc1, tr2, trMc2); } } - PROCESS_SWITCH(FwdMuonsUPC, processMcReco, "", false); + PROCESS_SWITCH(fwdMuonsUPC, processMcReco, "", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), }; }