From f88c945d4529cc128019387e7e623e22fdebb40d Mon Sep 17 00:00:00 2001 From: Sara Haidlova Date: Tue, 12 Mar 2024 16:30:53 +0100 Subject: [PATCH 1/6] PWGUD: new task for central barrel J/Psi --- PWGUD/Core/UDFSParser.h | 2 +- PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h | 221 +++++ PWGUD/Tasks/CMakeLists.txt | 4 + PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx | 846 ++++++++++++++++++++ 4 files changed, 1072 insertions(+), 1 deletion(-) create mode 100644 PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h create mode 100644 PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx diff --git a/PWGUD/Core/UDFSParser.h b/PWGUD/Core/UDFSParser.h index fb3d36ea524..6ceddafa5df 100644 --- a/PWGUD/Core/UDFSParser.h +++ b/PWGUD/Core/UDFSParser.h @@ -12,7 +12,7 @@ #ifndef PWGUD_CORE_UDFSPARSER_H_ #define PWGUD_CORE_UDFSPARSER_H_ -//#include +// #include #include #include diff --git a/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h b/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h new file mode 100644 index 00000000000..a8d0fa2af0a --- /dev/null +++ b/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h @@ -0,0 +1,221 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \brief +/// \author Sara Haidlova, sara.haidlova@cern.ch + +#ifndef PWGUD_CORE_UPCJPSICENTRALBARRELCORR_H_ +#define PWGUD_CORE_UPCJPSICENTRALBARRELCORR_H_ + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace std; + +enum ParticleType { + P_ELECTRON = 0, + P_MUON = 1, + P_PROTON = 2 +}; + +template +int testPIDhypo(T trackPID) +{ + float nSigmaTPC[3]; + nSigmaTPC[P_ELECTRON] = std::abs(trackPID.tpcNSigmaEl()); + nSigmaTPC[P_MUON] = std::abs(trackPID.tpcNSigmaMu()); + nSigmaTPC[P_PROTON] = std::abs(trackPID.tpcNSigmaPr()); + int enumChoiceTPC = std::distance(std::begin(nSigmaTPC), + std::min_element(std::begin(nSigmaTPC), std::end(nSigmaTPC))); + if (trackPID.hasTPC()) { + return enumChoiceTPC; + } else + return -1; +} + +float mom(float px, float py, float pz) +{ + return std::sqrt(px * px + py * py + pz * pz); +} + +float phi(float px, float py) +{ + if (px != 0) { + return std::atan2(py, px); + } else + return -999; +} + +float eta(float px, float py, float pz) +{ + float mom = std::sqrt(px * px + py * py + pz * pz); + if (mom != 0) { + return std::atanh(pz / mom); + } else + return -999; +} + +/*float energy(float mass, float px, float py, float pz) +// Just a simple function to return track energy +{ + return std::sqrt(mass * mass + px * px + py * py + pz * pz); +} + +float rap(float mass, float px, float py, float pz) +// Just a simple function to return track rapidity +{ + return 0.5 * std::log((energy(mass, px, py, pz) + pz) / (energy(mass, px, py, pz) - pz)); +}*/ + +float* correlation(TLorentzVector* lv1, TLorentzVector* lv2, TLorentzVector* lv) +{ + TLorentzVector pa(1., 0., 0, 1); // projectile + TLorentzVector pb(1., 0., 0, -1); // target + + float* q = new float[3]; // to save values + + // Accoplanarity angle + Float_t deltaPhi; + deltaPhi = lv1->Phi() - lv2->Phi(); + float accoplCut = 1. - TMath::Abs(deltaPhi) / TMath::Pi(); + // z + TLorentzVector z; + Float_t part1, part2; + + // Dot product: v1*v2 = t1*t2-x1*x2-y1*y2-z1*z2 + + part1 = lv->Dot(pb); + part2 = lv->Dot(pa); + + Float_t part3x = pa.X() * part1; + Float_t part3y = pa.Y() * part1; + Float_t part3z = pa.Z() * part1; + Float_t part3e = pa.T() * part1; + + Float_t part4x = pb.X() * part2; + Float_t part4y = pb.Y() * part2; + Float_t part4z = pb.Z() * part2; + Float_t part4e = pb.T() * part2; + + TLorentzVector part3(TVector3(part3x, part3y, part3z), part3e); + TLorentzVector part4(TVector3(part4x, part4y, part4z), part4e); + + // Un-normalized Z + z = part3 - part4; + + // Normalized z + Float_t normz = TMath::Sqrt(-z * z); + Float_t znx = z.X() / normz; + Float_t zny = z.Y() / normz; + Float_t znz = z.Z() / normz; + Float_t zne = z.E() / normz; + + // Normalized z + TLorentzVector zhat(TVector3(znx, zny, znz), zne); + + // calculate x + TLorentzVector x; + + Float_t constant1 = (lv->Dot(*lv)) / (2 * (lv->Dot(pa))); + Float_t constant2 = (lv->Dot(*lv)) / (2 * (lv->Dot(pb))); + + Float_t comp1x = pa.X() * constant1; + Float_t comp1y = pa.Y() * constant1; + Float_t comp1z = pa.Z() * constant1; + Float_t comp1e = pa.T() * constant1; + + TLorentzVector comp1(TVector3(comp1x, comp1y, comp1z), comp1e); + + Float_t comp2x = pb.X() * constant2; + Float_t comp2y = pb.Y() * constant2; + Float_t comp2z = pb.Z() * constant2; + Float_t comp2e = pb.T() * constant2; + + TLorentzVector comp2(TVector3(comp2x, comp2y, comp2z), comp2e); + + // Un-normalized x + x = *lv - comp1 - comp2; + // normalize x + Float_t normx = TMath::Sqrt(-x * x); + Float_t xnx = x.X() / normx; + Float_t xny = x.Y() / normx; + Float_t xnz = x.Z() / normx; + Float_t xne = x.E() / normx; + + // Normalized x + TLorentzVector xhat(TVector3(xnx, xny, xnz), xne); + + // calculate y + // TLorentzVector y; + Float_t yone = pa.Y() * pb.Z() * lv->E() - pa.Z() * pb.Y() * lv->E() + pa.Z() * pb.E() * lv->Y() + pa.E() * pb.Y() * lv->Z() - pa.Y() * pb.E() * lv->Z() - pa.E() * pb.Z() * lv->Y(); + Float_t ytwo = -pa.Z() * pb.E() * lv->X() + pa.Z() * pb.X() * lv->E() - pa.X() * pb.Z() * lv->E() + pa.X() * pb.E() * lv->Z() - pa.E() * pb.X() * lv->Z() + pa.E() * pb.Z() * lv->X(); + Float_t ythree = pa.X() * pb.Y() * lv->E() - pa.Y() * pb.X() * lv->E() + pa.Y() * pb.E() * lv->X() - pa.X() * pb.E() * lv->Y() + pa.E() * pb.X() * lv->Y() - pa.E() * pb.Y() * lv->X(); + Float_t yfour = -pa.X() * pb.Y() * lv->Z() + pa.X() * pb.Z() * lv->Y() - pa.Z() * pb.X() * lv->Y() + pa.Z() * pb.Y() * lv->X() - pa.Y() * pb.Z() * lv->X() + pa.Y() * pb.X() * lv->Z(); + + // Un-normalized y + TLorentzVector y(TVector3(yone, ytwo, ythree), yfour); + + // normalize y + Float_t normy = TMath::Sqrt(-y * y); + Float_t ynx = y.X() / normy; + Float_t yny = y.Y() / normy; + Float_t ynz = y.Z() / normy; + Float_t yne = y.E() / normy; + + // normalized y + TLorentzVector yhat(TVector3(ynx, yny, ynz), yne); + + // Lepton momentum difference + TLorentzVector diff; + diff = (*lv1 - *lv2); + Float_t diff2x = diff.X() / 2.; + Float_t diff2y = diff.Y() / 2.; + Float_t diff2z = diff.Z() / 2.; + Float_t diff2e = diff.E() / 2.; + TLorentzVector diff2(TVector3(diff2x, diff2y, diff2z), diff2e); + + // Normalize diff2 + Float_t norm2 = TMath::Sqrt(-diff2 * diff2); + Float_t diff3x = diff2.X() / norm2; + Float_t diff3y = diff2.Y() / norm2; + Float_t diff3z = diff2.Z() / norm2; + Float_t diff3e = diff2.E() / norm2; + + TLorentzVector diff3(TVector3(diff3x, diff3y, diff3z), diff3e); + + // computing the angles + float cosThetaCS = zhat * diff3; + Double_t SinThetaCosPhiCS = xhat * diff3; + Double_t SinThetaSinPhiCS = yhat * diff3; + //************************************** + + float phi = atan2(SinThetaSinPhiCS, SinThetaCosPhiCS); + // if (phi>=0) phi = phi; + // if (phi<0) phi = phi + 2*TMath::Pi(); + + q[0] = accoplCut; + q[1] = phi; + q[2] = cosThetaCS; + + return q; +} + +double DeltaPhi(TLorentzVector lv1, TLorentzVector lv2) +{ + TLorentzVector lv_sum = lv1 + lv2; + TLorentzVector lv_diff = lv1 - lv2; + + double dp = lv_sum.DeltaPhi(lv_diff); + + return dp; +} + +#endif // PWGUD_CORE_UPCJPSICENTRALBARRELCORR_H_ diff --git a/PWGUD/Tasks/CMakeLists.txt b/PWGUD/Tasks/CMakeLists.txt index edf72cb8984..bc672b3e2d8 100644 --- a/PWGUD/Tasks/CMakeLists.txt +++ b/PWGUD/Tasks/CMakeLists.txt @@ -71,4 +71,8 @@ o2physics_add_dpl_workflow(upc-tau-rl o2physics_add_dpl_workflow(polarisation-rho SOURCES PolarisationRho.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGPIDSelector + COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(upc-jpsi-corr + SOURCES upcJpsiCentralBarrelCorr.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGPIDSelector COMPONENT_NAME Analysis) \ No newline at end of file diff --git a/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx b/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx new file mode 100644 index 00000000000..4ff739baaf0 --- /dev/null +++ b/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx @@ -0,0 +1,846 @@ + +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \brief +/// \author Sara Haidlova, sara.haidlova@cern.ch + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" + +#include "TDatabasePDG.h" +#include "TLorentzVector.h" +#include "TMath.h" +#include "TVector3.h" +#include "Common/DataModel/PIDResponse.h" +#include "PWGUD/DataModel/UDTables.h" +#include "PWGUD/Core/UDHelpers.h" +#include +#include "UPCJpsiCentralBarrelCorrHelper.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::constants::math; + +struct upcJpsiCentralBarrel { + // configurable axes + ConfigurableAxis IVMAxis{"IVMAxis", {350.0f, 0.0f, 4.5f}, "M_#it{inv} (GeV/#it{c}^{2})"}; + ConfigurableAxis ptAxis{"ptAxis", {250.0f, 0.1f, 3.0f}, "#it{p}_T (GeV/#it{c})"}; + ConfigurableAxis pAxis{"pAxis", {250.0f, 0.1f, 3.0f}, "#it{p} (GeV/#it{c})"}; + ConfigurableAxis etaAxis{"etaAxis", {250.0f, -1.5f, 1.5f}, "#eta (-)"}; + ConfigurableAxis countAxis{"countAxis", {10.0f, 0.0f, 10.0f}, "Number of events (-)"}; + ConfigurableAxis phiAxis{"phiAxis", {250.0f, -PI, PI}, "#phi (rad)"}; + ConfigurableAxis accoplAxis{"accoplAxis", {250.0f, -0.2f, 0.2f}, "accAngle"}; + ConfigurableAxis thetaAxis{"thetaAxis", {250.0f, -1.5f, 1.5f}, "cos #theta (-)"}; + ConfigurableAxis sigTPCAxis{"sigTPCAxis", {100.0f, 0, 200.0f}, "TPC d#it{E}/d#it{x}"}; + ConfigurableAxis sigTOFAxis{"sigTOFAxis", {100.0f, 0, 200.0f}, "TOF d#it{E}/d#it{x}"}; + + // configurable cuts (modify in json) + Configurable TPCNClsCrossedRows{"TPCNClsCrossedRows", 70, "number of crossed rows in TPC"}; + Configurable TPCNSigmaMu{"TPCNSigmaMu", 5, "PID for TPC Mu track"}; + Configurable EtaCut{"EtaCut", 0.9f, "acceptance cut per track"}; + Configurable RapCut{"RapCut", 0.8f, "choose event in midrapidity"}; + Configurable dcaZCut{"dcaZCut", 2, "cut on the impact parameter in z of the track to the PV"}; + Configurable dcaXYCut{"dcaXYCut", 1e10, "cut on the impact parameter in xy of the track to the PV"}; + Configurable ITSNClsCut{"ITSNClsCut", 1, "minimal number of ITS clusters"}; + Configurable ITSChi2NClsCut{"ITSChi2NClsCut", 36, "minimal Chi2/cluster for the ITS track"}; + Configurable TPCNClsCrossedRowsCut{"TPCNClsCrossedRowsCut", 70, "minimal number of crossed TPC rows"}; + Configurable TPCChi2NCls{"TPCChi2NCls", 4, "minimal Chi2/cluster for the TPC track"}; + + TDatabasePDG* pdg = nullptr; + + // initialize histogram registry + HistogramRegistry Statistics{ + "Statistics", + {}}; + + HistogramRegistry RawData{ + "RawData", + {}}; + + HistogramRegistry PVContributors{ + "PVContributors", + {}}; + + HistogramRegistry TGmu{ + "TGmu", + {}}; + + HistogramRegistry TGmuCand{ + "TGmuCand", + {}}; + + HistogramRegistry TGel{ + "TGel", + {}}; + + HistogramRegistry TGelCand{ + "TGelCand", + {}}; + + HistogramRegistry TGp{ + "TGp", + {}}; + + HistogramRegistry TGpCand{ + "TGpCand", + {}}; + + HistogramRegistry JPsiToEl{ + "JPsiToEl", + {}}; + + HistogramRegistry JPsiToMu{ + "JPsiToMu", + {}}; + + HistogramRegistry JPsiToP{ + "JPsiToP", + {}}; + + HistogramRegistry Correlation{ + "Correlation", + {}}; + + HistogramRegistry Asymmetry{ + "Asymmetry", + {}}; + + using UDCollisionsFull = soa::Join; + using UDCollisionFull = UDCollisionsFull::iterator; + using UDTracksFull = soa::Join; + using UDTrackFull = UDTracksFull::iterator; + + void init(InitContext&) + { + pdg = TDatabasePDG::Instance(); + + const AxisSpec axisIVM{IVMAxis, "IVM axis"}; + const AxisSpec axispt{ptAxis, "pt axis"}; + const AxisSpec axiseta{etaAxis, "eta axis"}; + const AxisSpec axisCounter{countAxis, "counter axis"}; + const AxisSpec axisAccAngle{accoplAxis, "accAngle"}; + const AxisSpec axisAngTheta{thetaAxis, "cosTheta"}; + const AxisSpec axisPhi{phiAxis, "phi"}; + const AxisSpec axisp{pAxis, "p axis"}; + const AxisSpec axisTPC{sigTPCAxis, ""}; + const AxisSpec axisTOF{sigTOFAxis, ""}; + + // statistics histograms for counters + Statistics.add("Statistics/hNumberOfCollisions", "hNumberOfCollisions", {HistType::kTH1F, {axisCounter}}); + Statistics.add("Statistics/hNumberOfTracks", "hNumberOfTracks", {HistType::kTH1F, {axisCounter}}); + Statistics.add("Statistics/hNumberGT", "hNumberGT", {HistType::kTH1F, {axisCounter}}); + Statistics.add("Statistics/hNumberGTselected", "hNumberGTselected", {HistType::kTH1F, {axisCounter}}); + Statistics.add("Statistics/hNumberGTel", "hNumberGTel", {HistType::kTH1F, {axisCounter}}); + Statistics.add("Statistics/hNumberGTmu", "hNumberGTmu", {HistType::kTH1F, {axisCounter}}); + Statistics.add("Statistics/hNumberGTp", "hNumberGTp", {HistType::kTH1F, {axisCounter}}); + + // raw data histograms + RawData.add("RawData/hTrackPt", "hTrackPt", {HistType::kTH1F, {axispt}}); + RawData.add("RawData/hTrackEta", "hTrackEta", {HistType::kTH1F, {axiseta}}); + RawData.add("RawData/hTrackPhi", "hTrackPhi", {HistType::kTH1F, {axisPhi}}); + RawData.add("RawData/PID/hTPCVsP", "hTPCVsP", {HistType::kTH2F, {axisp, axisTPC}}); + RawData.add("RawData/PID/hTPCVsPt", "hTPCVsPt", {HistType::kTH2F, {axispt, axisTPC}}); + RawData.add("RawData/PID/hTPCVsPhi", "hTPCVsPhi", {HistType::kTH2F, {axisPhi, axisTPC}}); + RawData.add("RawData/PID/hTPCVsEta", "hTPCVsEta", {HistType::kTH2F, {axiseta, axisTPC}}); + RawData.add("RawData/PID/hTOFVsP", "hTOFVsP", {HistType::kTH2F, {axisp, axisTOF}}); + RawData.add("RawData/PID/hTOFVsPt", "hTOFVsPt", {HistType::kTH2F, {axispt, axisTOF}}); + RawData.add("RawData/PID/hTOFVsPhi", "hTOFVsPhi", {HistType::kTH2F, {axisPhi, axisTOF}}); + RawData.add("RawData/PID/hTOFVsEta", "hTOFVsEta", {HistType::kTH2F, {axiseta, axisTOF}}); + + // PVContributors histograms + PVContributors.add("PVContributors/hTrackPt", "hTrackPt", {HistType::kTH1F, {axispt}}); + PVContributors.add("PVContributors/hTrackEta", "hTrackEta", {HistType::kTH1F, {axiseta}}); + PVContributors.add("PVContributors/hTrackPhi", "hTrackPhi", {HistType::kTH1F, {axisPhi}}); + PVContributors.add("PVContributors/PID/hTPCVsP", "hTPCVsP", {HistType::kTH2F, {axisp, axisTPC}}); + PVContributors.add("PVContributors/PID/hTPCVsPt", "hTPCVsPt", {HistType::kTH2F, {axispt, axisTPC}}); + PVContributors.add("PVContributors/PID/hTPCVsPhi", "hTPCVsPhi", {HistType::kTH2F, {axisPhi, axisTPC}}); + PVContributors.add("PVContributors/PID/hTPCVsEta", "hTPCVsEta", {HistType::kTH2F, {axiseta, axisTPC}}); + PVContributors.add("PVContributors/PID/hTOFVsP", "hTOFVsP", {HistType::kTH2F, {axisp, axisTOF}}); + PVContributors.add("PVContributors/PID/hTOFVsPt", "hTOFVsPt", {HistType::kTH2F, {axispt, axisTOF}}); + PVContributors.add("PVContributors/PID/hTOFVsPhi", "hTOFVsPhi", {HistType::kTH2F, {axisPhi, axisTOF}}); + PVContributors.add("PVContributors/PID/hTOFVsEta", "hTOFVsEta", {HistType::kTH2F, {axiseta, axisTOF}}); + + // TGmu histograms + TGmu.add("TGmu/hTrackPt", "hTrackPt", {HistType::kTH1F, {axispt}}); + TGmu.add("TGmu/hTrackEta", "hTrackEta", {HistType::kTH1F, {axiseta}}); + TGmu.add("TGmu/hTrackPhi", "hTrackPhi", {HistType::kTH1F, {axisPhi}}); + + // TGmuCand histograms + TGmuCand.add("TGmuCand/hTrackPt", "hTrackPt", {HistType::kTH1F, {axispt}}); + TGmuCand.add("TGmuCand/hTrackEta", "hTrackEta", {HistType::kTH1F, {axiseta}}); + TGmuCand.add("TGmuCand/hTrackPhi", "hTrackPhi", {HistType::kTH1F, {axisPhi}}); + TGmuCand.add("TGmuCand/hJpsiPt", "hJpsiPt", {HistType::kTH1F, {axispt}}); + + // TGel histograms + TGel.add("TGel/hTrackPt", "hTrackPt", {HistType::kTH1F, {axispt}}); + TGel.add("TGel/hTrackEta", "hTrackEta", {HistType::kTH1F, {axiseta}}); + TGel.add("TGel/hTrackPhi", "hTrackPhi", {HistType::kTH1F, {axisPhi}}); + + // TGelCand histograms + TGelCand.add("TGelCand/hTrackPt", "hTrackPt", {HistType::kTH1F, {axispt}}); + TGelCand.add("TGelCand/hTrackEta", "hTrackEta", {HistType::kTH1F, {axiseta}}); + TGelCand.add("TGelCand/hTrackPhi", "hTrackPhi", {HistType::kTH1F, {axisPhi}}); + TGelCand.add("TGelCand/hJpsiPt", "hJpsiPt", {HistType::kTH1F, {axispt}}); + + // TGp histograms + TGp.add("TGp/hTrackPt", "hTrackPt", {HistType::kTH1F, {axispt}}); + TGp.add("TGp/hTrackEta", "hTrackEta", {HistType::kTH1F, {axiseta}}); + TGp.add("TGp/hTrackPhi", "hTrackPhi", {HistType::kTH1F, {axisPhi}}); + + // TGpCand histograms + TGpCand.add("TGpCand/hTrackPt", "hTrackPt", {HistType::kTH1F, {axispt}}); + TGpCand.add("TGpCand/hTrackEta", "hTrackEta", {HistType::kTH1F, {axiseta}}); + TGpCand.add("TGpCand/hTrackPhi", "hTrackPhi", {HistType::kTH1F, {axisPhi}}); + TGpCand.add("TGpCand/hJpsiPt", "hJpsiPt", {HistType::kTH1F, {axispt}}); + + // JPsiToEl histograms + JPsiToEl.add("JPsiToEl/Coherent/hPt", "Pt of J/Psi ; p_{T} {GeV/c]", {HistType::kTH1F, {axispt}}); + JPsiToEl.add("JPsiToEl/Coherent/hPt1", "pT of track 1 ; p_{T} {GeV/c]", {HistType::kTH1F, {axispt}}); + JPsiToEl.add("JPsiToEl/Coherent/hPt2", "pT of track 2 ; p_{T} {GeV/c]", {HistType::kTH1F, {axispt}}); + JPsiToEl.add("JPsiToEl/Coherent/hEta1", "eta of track 1 ; #eta {-]", {HistType::kTH1F, {axiseta}}); + JPsiToEl.add("JPsiToEl/Coherent/hEta2", "eta of track 2 ; #eta {-]", {HistType::kTH1F, {axiseta}}); + JPsiToEl.add("JPsiToEl/Coherent/hPhi1", "phi of track 1 ; #phi {-]", {HistType::kTH1F, {axisPhi}}); + JPsiToEl.add("JPsiToEl/Coherent/hPhi2", "phi of track 2 ; #phi {-]", {HistType::kTH1F, {axisPhi}}); + JPsiToEl.add("JPsiToEl/Coherent/hIVM", "J/Psi Invariant Mass ; m {GeV]", {HistType::kTH1F, {axisIVM}}); + JPsiToEl.add("JPsiToEl/Coherent/hRap", "Rap of J/Psi ; y {-]", {HistType::kTH1F, {axiseta}}); + JPsiToEl.add("JPsiToEl/Coherent/hEta", "Eta of J/Psi ; #eta {-]", {HistType::kTH1F, {axiseta}}); + JPsiToEl.add("JPsiToEl/Coherent/hPhi", "Phi of J/Psi ; #phi {-]", {HistType::kTH1F, {axisPhi}}); + JPsiToEl.add("JPsiToEl/Coherent/PID/hTPCVsP", "hTPCVsP", {HistType::kTH2F, {axisp, axisTPC}}); + JPsiToEl.add("JPsiToEl/Coherent/PID/hTPCVsPt", "hTPCVsPt", {HistType::kTH2F, {axispt, axisTPC}}); + JPsiToEl.add("JPsiToEl/Coherent/PID/hTPCVsPhi", "hTPCVsPhi", {HistType::kTH2F, {axisPhi, axisTPC}}); + JPsiToEl.add("JPsiToEl/Coherent/PID/hTPCVsEta", "hTPCVsEta", {HistType::kTH2F, {axiseta, axisTPC}}); + JPsiToEl.add("JPsiToEl/Coherent/PID/hTOFVsP", "hTOFVsP", {HistType::kTH2F, {axisp, axisTOF}}); + JPsiToEl.add("JPsiToEl/Coherent/PID/hTOFVsPt", "hTOFVsPt", {HistType::kTH2F, {axispt, axisTOF}}); + JPsiToEl.add("JPsiToEl/Coherent/PID/hTOFVsPhi", "hTOFVsPhi", {HistType::kTH2F, {axisPhi, axisTOF}}); + JPsiToEl.add("JPsiToEl/Coherent/PID/hTOFVsEta", "hTOFVsEta", {HistType::kTH2F, {axiseta, axisTOF}}); + + JPsiToEl.add("JPsiToEl/Incoherent/hPt", "Pt of J/Psi ; p_{T} {GeV/c]", {HistType::kTH1F, {axispt}}); + JPsiToEl.add("JPsiToEl/Incoherent/hPt1", "pT of track 1 ; p_{T} {GeV/c]", {HistType::kTH1F, {axispt}}); + JPsiToEl.add("JPsiToEl/Incoherent/hPt2", "pT of track 2 ; p_{T} {GeV/c]", {HistType::kTH1F, {axispt}}); + JPsiToEl.add("JPsiToEl/Incoherent/hEta1", "eta of track 1 ; #eta {-]", {HistType::kTH1F, {axiseta}}); + JPsiToEl.add("JPsiToEl/Incoherent/hEta2", "eta of track 2 ; #eta {-]", {HistType::kTH1F, {axiseta}}); + JPsiToEl.add("JPsiToEl/Incoherent/hPhi1", "phi of track 1 ; #phi {-]", {HistType::kTH1F, {axisPhi}}); + JPsiToEl.add("JPsiToEl/Incoherent/hPhi2", "phi of track 2 ; #phi {-]", {HistType::kTH1F, {axisPhi}}); + JPsiToEl.add("JPsiToEl/Incoherent/hIVM", "J/Psi Invariant Mass ; m {GeV]", {HistType::kTH1F, {axisIVM}}); + JPsiToEl.add("JPsiToEl/Incoherent/hRap", "Rap of J/Psi ; y {-]", {HistType::kTH1F, {axiseta}}); + JPsiToEl.add("JPsiToEl/Incoherent/hEta", "Eta of J/Psi ; #eta {-]", {HistType::kTH1F, {axiseta}}); + JPsiToEl.add("JPsiToEl/Incoherent/hPhi", "Phi of J/Psi ; #phi {-]", {HistType::kTH1F, {axisPhi}}); + JPsiToEl.add("JPsiToEl/Incoherent/PID/hTPCVsP", "hTPCVsP", {HistType::kTH2F, {axisp, axisTPC}}); + JPsiToEl.add("JPsiToEl/Incoherent/PID/hTPCVsPt", "hTPCVsPt", {HistType::kTH2F, {axispt, axisTPC}}); + JPsiToEl.add("JPsiToEl/Incoherent/PID/hTPCVsPhi", "hTPCVsPhi", {HistType::kTH2F, {axisPhi, axisTPC}}); + JPsiToEl.add("JPsiToEl/Incoherent/PID/hTPCVsEta", "hTPCVsEta", {HistType::kTH2F, {axiseta, axisTPC}}); + JPsiToEl.add("JPsiToEl/Incoherent/PID/hTOFVsP", "hTOFVsP", {HistType::kTH2F, {axisp, axisTOF}}); + JPsiToEl.add("JPsiToEl/Incoherent/PID/hTOFVsPt", "hTOFVsPt", {HistType::kTH2F, {axispt, axisTOF}}); + JPsiToEl.add("JPsiToEl/Incoherent/PID/hTOFVsPhi", "hTOFVsPhi", {HistType::kTH2F, {axisPhi, axisTOF}}); + JPsiToEl.add("JPsiToEl/Incoherent/PID/hTOFVsEta", "hTOFVsEta", {HistType::kTH2F, {axiseta, axisTOF}}); + + // JPsiToMu histograms + JPsiToMu.add("JPsiToMu/Coherent/hPt", "Pt of J/Psi ; p_{T} {GeV/c]", {HistType::kTH1F, {axispt}}); + JPsiToMu.add("JPsiToMu/Coherent/hPt1", "pT of track 1 ; p_{T} {GeV/c]", {HistType::kTH1F, {axispt}}); + JPsiToMu.add("JPsiToMu/Coherent/hPt2", "pT of track 2 ; p_{T} {GeV/c]", {HistType::kTH1F, {axispt}}); + JPsiToMu.add("JPsiToMu/Coherent/hEta1", "eta of track 1 ; #eta {-]", {HistType::kTH1F, {axiseta}}); + JPsiToMu.add("JPsiToMu/Coherent/hEta2", "eta of track 2 ; #eta {-]", {HistType::kTH1F, {axiseta}}); + JPsiToMu.add("JPsiToMu/Coherent/hPhi1", "phi of track 1 ; #phi {-]", {HistType::kTH1F, {axisPhi}}); + JPsiToMu.add("JPsiToMu/Coherent/hPhi2", "phi of track 2 ; #phi {-]", {HistType::kTH1F, {axisPhi}}); + JPsiToMu.add("JPsiToMu/Coherent/hIVM", "J/Psi Invariant Mass ; m {GeV]", {HistType::kTH1F, {axisIVM}}); + JPsiToMu.add("JPsiToMu/Coherent/hRap", "Rap of J/Psi ; y {-]", {HistType::kTH1F, {axiseta}}); + JPsiToMu.add("JPsiToMu/Coherent/hEta", "Eta of J/Psi ; #eta {-]", {HistType::kTH1F, {axiseta}}); + JPsiToMu.add("JPsiToMu/Coherent/hPhi", "Phi of J/Psi ; #phi {-]", {HistType::kTH1F, {axisPhi}}); + JPsiToMu.add("JPsiToMu/Coherent/PID/hTPCVsP", "hTPCVsP", {HistType::kTH2F, {axisp, axisTPC}}); + JPsiToMu.add("JPsiToMu/Coherent/PID/hTPCVsPt", "hTPCVsPt", {HistType::kTH2F, {axispt, axisTPC}}); + JPsiToMu.add("JPsiToMu/Coherent/PID/hTPCVsPhi", "hTPCVsPhi", {HistType::kTH2F, {axisPhi, axisTPC}}); + JPsiToMu.add("JPsiToMu/Coherent/PID/hTPCVsEta", "hTPCVsEta", {HistType::kTH2F, {axiseta, axisTPC}}); + JPsiToMu.add("JPsiToMu/Coherent/PID/hTOFVsP", "hTOFVsP", {HistType::kTH2F, {axisp, axisTOF}}); + JPsiToMu.add("JPsiToMu/Coherent/PID/hTOFVsPt", "hTOFVsPt", {HistType::kTH2F, {axispt, axisTOF}}); + JPsiToMu.add("JPsiToMu/Coherent/PID/hTOFVsPhi", "hTOFVsPhi", {HistType::kTH2F, {axisPhi, axisTOF}}); + JPsiToMu.add("JPsiToMu/Coherent/PID/hTOFVsEta", "hTOFVsEta", {HistType::kTH2F, {axiseta, axisTOF}}); + + JPsiToMu.add("JPsiToMu/Incoherent/hPt", "Pt of J/Psi ; p_{T} {GeV/c]", {HistType::kTH1F, {axispt}}); + JPsiToMu.add("JPsiToMu/Incoherent/hPt1", "pT of track 1 ; p_{T} {GeV/c]", {HistType::kTH1F, {axispt}}); + JPsiToMu.add("JPsiToMu/Incoherent/hPt2", "pT of track 2 ; p_{T} {GeV/c]", {HistType::kTH1F, {axispt}}); + JPsiToMu.add("JPsiToMu/Incoherent/hEta1", "eta of track 1 ; #eta {-]", {HistType::kTH1F, {axiseta}}); + JPsiToMu.add("JPsiToMu/Incoherent/hEta2", "eta of track 2 ; #eta {-]", {HistType::kTH1F, {axiseta}}); + JPsiToMu.add("JPsiToMu/Incoherent/hPhi1", "phi of track 1 ; #phi {-]", {HistType::kTH1F, {axisPhi}}); + JPsiToMu.add("JPsiToMu/Incoherent/hPhi2", "phi of track 2 ; #phi {-]", {HistType::kTH1F, {axisPhi}}); + JPsiToMu.add("JPsiToMu/Incoherent/hIVM", "J/Psi Invariant Mass ; m {GeV]", {HistType::kTH1F, {axisIVM}}); + JPsiToMu.add("JPsiToMu/Incoherent/hRap", "Rap of J/Psi ; y {-]", {HistType::kTH1F, {axiseta}}); + JPsiToMu.add("JPsiToMu/Incoherent/hEta", "Eta of J/Psi ; #eta {-]", {HistType::kTH1F, {axiseta}}); + JPsiToMu.add("JPsiToMu/Incoherent/hPhi", "Phi of J/Psi ; #phi {-]", {HistType::kTH1F, {axisPhi}}); + JPsiToMu.add("JPsiToMu/Incoherent/PID/hTPCVsP", "hTPCVsP", {HistType::kTH2F, {axisp, axisTPC}}); + JPsiToMu.add("JPsiToMu/Incoherent/PID/hTPCVsPt", "hTPCVsPt", {HistType::kTH2F, {axispt, axisTPC}}); + JPsiToMu.add("JPsiToMu/Incoherent/PID/hTPCVsPhi", "hTPCVsPhi", {HistType::kTH2F, {axisPhi, axisTPC}}); + JPsiToMu.add("JPsiToMu/Incoherent/PID/hTPCVsEta", "hTPCVsEta", {HistType::kTH2F, {axiseta, axisTPC}}); + JPsiToMu.add("JPsiToMu/Incoherent/PID/hTOFVsP", "hTOFVsP", {HistType::kTH2F, {axisp, axisTOF}}); + JPsiToMu.add("JPsiToMu/Incoherent/PID/hTOFVsPt", "hTOFVsPt", {HistType::kTH2F, {axispt, axisTOF}}); + JPsiToMu.add("JPsiToMu/Incoherent/PID/hTOFVsPhi", "hTOFVsPhi", {HistType::kTH2F, {axisPhi, axisTOF}}); + JPsiToMu.add("JPsiToMu/Incoherent/PID/hTOFVsEta", "hTOFVsEta", {HistType::kTH2F, {axiseta, axisTOF}}); + + // JPsiToP histograms + JPsiToP.add("JPsiToP/Coherent/hPt", "Pt of J/Psi ; p_{T} {GeV/c]", {HistType::kTH1F, {axispt}}); + JPsiToP.add("JPsiToP/Coherent/hPt1", "pT of track 1 ; p_{T} {GeV/c]", {HistType::kTH1F, {axispt}}); + JPsiToP.add("JPsiToP/Coherent/hPt2", "pT of track 2 ; p_{T} {GeV/c]", {HistType::kTH1F, {axispt}}); + JPsiToP.add("JPsiToP/Coherent/hEta1", "eta of track 1 ; #eta {-]", {HistType::kTH1F, {axiseta}}); + JPsiToP.add("JPsiToP/Coherent/hEta2", "eta of track 2 ; #eta {-]", {HistType::kTH1F, {axiseta}}); + JPsiToP.add("JPsiToP/Coherent/hPhi1", "phi of track 1 ; #phi {-]", {HistType::kTH1F, {axisPhi}}); + JPsiToP.add("JPsiToP/Coherent/hPhi2", "phi of track 2 ; #phi {-]", {HistType::kTH1F, {axisPhi}}); + JPsiToP.add("JPsiToP/Coherent/hIVM", "J/Psi Invariant Mass ; m {GeV]", {HistType::kTH1F, {axisIVM}}); + JPsiToP.add("JPsiToP/Coherent/hRap", "Rap of J/Psi ; y {-]", {HistType::kTH1F, {axiseta}}); + JPsiToP.add("JPsiToP/Coherent/hEta", "Eta of J/Psi ; #eta {-]", {HistType::kTH1F, {axiseta}}); + JPsiToP.add("JPsiToP/Coherent/hPhi", "Phi of J/Psi ; #phi {-]", {HistType::kTH1F, {axisPhi}}); + JPsiToP.add("JPsiToP/Coherent/PID/hTPCVsP", "hTPCVsP", {HistType::kTH2F, {axisp, axisTPC}}); + JPsiToP.add("JPsiToP/Coherent/PID/hTPCVsPt", "hTPCVsPt", {HistType::kTH2F, {axispt, axisTPC}}); + JPsiToP.add("JPsiToP/Coherent/PID/hTPCVsPhi", "hTPCVsPhi", {HistType::kTH2F, {axisPhi, axisTPC}}); + JPsiToP.add("JPsiToP/Coherent/PID/hTPCVsEta", "hTPCVsEta", {HistType::kTH2F, {axiseta, axisTPC}}); + JPsiToP.add("JPsiToP/Coherent/PID/hTOFVsP", "hTOFVsP", {HistType::kTH2F, {axisp, axisTOF}}); + JPsiToP.add("JPsiToP/Coherent/PID/hTOFVsPt", "hTOFVsPt", {HistType::kTH2F, {axispt, axisTOF}}); + JPsiToP.add("JPsiToP/Coherent/PID/hTOFVsPhi", "hTOFVsPhi", {HistType::kTH2F, {axisPhi, axisTOF}}); + JPsiToP.add("JPsiToP/Coherent/PID/hTOFVsEta", "hTOFVsEta", {HistType::kTH2F, {axiseta, axisTOF}}); + + JPsiToP.add("JPsiToP/Incoherent/hPt", "Pt of J/Psi ; p_{T} {GeV/c]", {HistType::kTH1F, {axispt}}); + JPsiToP.add("JPsiToP/Incoherent/hPt1", "pT of track 1 ; p_{T} {GeV/c]", {HistType::kTH1F, {axispt}}); + JPsiToP.add("JPsiToP/Incoherent/hPt2", "pT of track 2 ; p_{T} {GeV/c]", {HistType::kTH1F, {axispt}}); + JPsiToP.add("JPsiToP/Incoherent/hEta1", "eta of track 1 ; #eta {-]", {HistType::kTH1F, {axiseta}}); + JPsiToP.add("JPsiToP/Incoherent/hEta2", "eta of track 2 ; #eta {-]", {HistType::kTH1F, {axiseta}}); + JPsiToP.add("JPsiToP/Incoherent/hPhi1", "phi of track 1 ; #phi {-]", {HistType::kTH1F, {axisPhi}}); + JPsiToP.add("JPsiToP/Incoherent/hPhi2", "phi of track 2 ; #phi {-]", {HistType::kTH1F, {axisPhi}}); + JPsiToP.add("JPsiToP/Incoherent/hIVM", "J/Psi Invariant Mass ; m {GeV]", {HistType::kTH1F, {axisIVM}}); + JPsiToP.add("JPsiToP/Incoherent/hRap", "Rap of J/Psi ; y {-]", {HistType::kTH1F, {axiseta}}); + JPsiToP.add("JPsiToP/Incoherent/hEta", "Eta of J/Psi ; #eta {-]", {HistType::kTH1F, {axiseta}}); + JPsiToP.add("JPsiToP/Incoherent/hPhi", "Phi of J/Psi ; #phi {-]", {HistType::kTH1F, {axisPhi}}); + JPsiToP.add("JPsiToP/Incoherent/PID/hTPCVsP", "hTPCVsP", {HistType::kTH2F, {axisp, axisTPC}}); + JPsiToP.add("JPsiToP/Incoherent/PID/hTPCVsPt", "hTPCVsPt", {HistType::kTH2F, {axispt, axisTPC}}); + JPsiToP.add("JPsiToP/Incoherent/PID/hTPCVsPhi", "hTPCVsPhi", {HistType::kTH2F, {axisPhi, axisTPC}}); + JPsiToP.add("JPsiToP/Incoherent/PID/hTPCVsEta", "hTPCVsEta", {HistType::kTH2F, {axiseta, axisTPC}}); + JPsiToP.add("JPsiToP/Incoherent/PID/hTOFVsP", "hTOFVsP", {HistType::kTH2F, {axisp, axisTOF}}); + JPsiToP.add("JPsiToP/Incoherent/PID/hTOFVsPt", "hTOFVsPt", {HistType::kTH2F, {axispt, axisTOF}}); + JPsiToP.add("JPsiToP/Incoherent/PID/hTOFVsPhi", "hTOFVsPhi", {HistType::kTH2F, {axisPhi, axisTOF}}); + JPsiToP.add("JPsiToP/Incoherent/PID/hTOFVsEta", "hTOFVsEta", {HistType::kTH2F, {axiseta, axisTOF}}); + + // Correlation histograms + Correlation.add("Correlation/Muon/Coherent/AccoplAngle", "AccoplAngle", {HistType::kTH1F, {axisAccAngle}}); + Correlation.add("Correlation/Muon/Coherent/CosTheta", "CosTheta", {HistType::kTH1F, {axisAngTheta}}); + Correlation.add("Correlation/Muon/Coherent/Phi", "Phi", {HistType::kTH1F, {axisPhi}}); + Correlation.add("Correlation/Muon/Coherent/Phi1", "Phi1", {HistType::kTH1F, {axisPhi}}); + Correlation.add("Correlation/Muon/Coherent/Phi2", "Phi2", {HistType::kTH1F, {axisPhi}}); + Correlation.add("Correlation/Muon/Coherent/CosThetaPhi", "CosThetaPhi", {HistType::kTH2F, {{axisAngTheta}, {axisPhi}}}); + + Correlation.add("Correlation/Muon/Incoherent/AccoplAngle", "AccoplAngle", {HistType::kTH1F, {axisAccAngle}}); + Correlation.add("Correlation/Muon/Incoherent/CosTheta", "CosTheta", {HistType::kTH1F, {axisAngTheta}}); + Correlation.add("Correlation/Muon/Incoherent/Phi", "Phi", {HistType::kTH1F, {axisPhi}}); + Correlation.add("Correlation/Muon/Incoherent/Phi1", "Phi1", {HistType::kTH1F, {axisPhi}}); + Correlation.add("Correlation/Muon/Incoherent/Phi2", "Phi2", {HistType::kTH1F, {axisPhi}}); + Correlation.add("Correlation/Muon/Incoherent/CosThetaPhi", "CosThetaPhi", {HistType::kTH2F, {{axisAngTheta}, {axisPhi}}}); + + Correlation.add("Correlation/Electron/Coherent/AccoplAngle", "AccoplAngle", {HistType::kTH1F, {axisAccAngle}}); + Correlation.add("Correlation/Electron/Coherent/CosTheta", "CosTheta", {HistType::kTH1F, {axisAngTheta}}); + Correlation.add("Correlation/Electron/Coherent/Phi", "Phi", {HistType::kTH1F, {axisPhi}}); + Correlation.add("Correlation/Electron/Coherent/Phi1", "Phi1", {HistType::kTH1F, {axisPhi}}); + Correlation.add("Correlation/Electron/Coherent/Phi2", "Phi2", {HistType::kTH1F, {axisPhi}}); + Correlation.add("Correlation/Electron/Coherent/CosThetaPhi", "CosThetaPhi", {HistType::kTH2F, {{axisAngTheta}, {axisPhi}}}); + + Correlation.add("Correlation/Electron/Incoherent/AccoplAngle", "AccoplAngle", {HistType::kTH1F, {axisAccAngle}}); + Correlation.add("Correlation/Electron/Incoherent/CosTheta", "CosTheta", {HistType::kTH1F, {axisAngTheta}}); + Correlation.add("Correlation/Electron/Incoherent/Phi", "Phi", {HistType::kTH1F, {axisPhi}}); + Correlation.add("Correlation/Electron/Incoherent/Phi1", "Phi1", {HistType::kTH1F, {axisPhi}}); + Correlation.add("Correlation/Electron/Incoherent/Phi2", "Phi2", {HistType::kTH1F, {axisPhi}}); + Correlation.add("Correlation/Electron/Incoherent/CosThetaPhi", "CosThetaPhi", {HistType::kTH2F, {{axisAngTheta}, {axisPhi}}}); + + // Asymmetry histograms + Asymmetry.add("Asymmetry/Muon/Coherent/DeltaPhi", "DeltaPhi", {HistType::kTH1F, {axisPhi}}); + Asymmetry.add("Asymmetry/Muon/Incoherent/DeltaPhi", "DeltaPhi", {HistType::kTH1F, {axisPhi}}); + Asymmetry.add("Asymmetry/Electron/Coherent/DeltaPhi", "DeltaPhi", {HistType::kTH1F, {axisPhi}}); + Asymmetry.add("Asymmetry/Electron/Incoherent/DeltaPhi", "DeltaPhi", {HistType::kTH1F, {axisPhi}}); + } + + float particleMass(TDatabasePDG* pdg, int pid) + { + auto mass = 0.; + TParticlePDG* pdgparticle = pdg->GetParticle(pid); + if (pdgparticle != nullptr) { + mass = pdgparticle->Mass(); + } + return mass; + } + + template + bool GoodTrackCuts(T const& track) + { + // kinematics + if (eta(track.px(), track.py(), track.pz()) < -EtaCut || eta(track.px(), track.py(), track.pz()) > EtaCut) { + return false; + } + // DCA + if (track.dcaZ() > dcaZCut) { + return false; + } + if (track.dcaXY() > dcaXYCut) { + return false; + } + // ITS + if (!track.hasITS()) { + return false; + } + if (track.itsNCls() < ITSNClsCut) { + return false; + } + if (track.itsChi2NCl() > ITSChi2NClsCut) { + return false; + } + // TPC + if (!track.hasTPC()) { + return false; + } + if (track.tpcNClsCrossedRows() < TPCNClsCrossedRowsCut) { + return false; + } + if (track.tpcChi2NCl() > TPCChi2NCls) { + return false; // TPC chi2 + } + + return true; + } + + template + bool CandidateCuts(C const& mother) + { + if (abs(mother.Rapidity()) > 0.8f) { + return false; + } + + if (mother.M() < 2.5f) { + return false; + } + + return true; + } + + void process(UDCollisionFull const& collision, UDTracksFull const& tracks) + { + Statistics.get(HIST("Statistics/hNumberOfCollisions"))->Fill(0); // number of collisions without any cuts + + // loop over tracks without selections + for (auto& track : tracks) { + float trkPx = track.px(); + float trkPy = track.py(); + float trkPz = track.pz(); + + Statistics.get(HIST("Statistics/hNumberOfTracks"))->Fill(0); + if (track.isPVContributor() == 1) { + Statistics.get(HIST("Statistics/hNumberOfTracks"))->Fill(1); + PVContributors.get(HIST("PVContributors/hTrackPt"))->Fill(track.pt()); + PVContributors.get(HIST("PVContributors/hTrackEta"))->Fill(eta(trkPx, trkPy, trkPz)); + PVContributors.get(HIST("PVContributors/hTrackPhi"))->Fill(phi(trkPx, trkPy)); + + if (track.hasTPC()) { + PVContributors.get(HIST("PVContributors/PID/hTPCVsP"))->Fill(mom(trkPx, trkPy, trkPz), track.tpcSignal()); + PVContributors.get(HIST("PVContributors/PID/hTPCVsPt"))->Fill(track.pt(), track.tpcSignal()); + PVContributors.get(HIST("PVContributors/PID/hTPCVsEta"))->Fill(eta(trkPx, trkPy, trkPz), track.tpcSignal()); + PVContributors.get(HIST("PVContributors/PID/hTPCVsPhi"))->Fill(phi(trkPx, trkPy), track.tpcSignal()); + } + + if (track.hasTOF()) { + PVContributors.get(HIST("PVContributors/PID/hTOFVsP"))->Fill(mom(trkPx, trkPy, trkPz), track.tofSignal()); + PVContributors.get(HIST("PVContributors/PID/hTOFVsPt"))->Fill(track.pt(), track.tofSignal()); + PVContributors.get(HIST("PVContributors/PID/hTOFVsEta"))->Fill(eta(trkPx, trkPy, trkPz), track.tofSignal()); + PVContributors.get(HIST("PVContributors/PID/hTOFVsPhi"))->Fill(phi(trkPx, trkPy), track.tofSignal()); + } + } + + RawData.get(HIST("RawData/hTrackPt"))->Fill(track.pt()); + RawData.get(HIST("RawData/hTrackEta"))->Fill(eta(trkPx, trkPy, trkPz)); + RawData.get(HIST("RawData/hTrackPhi"))->Fill(phi(trkPx, trkPy)); + + if (track.hasTPC()) { + RawData.get(HIST("RawData/PID/hTPCVsP"))->Fill(mom(trkPx, trkPy, trkPz), track.tpcSignal()); + RawData.get(HIST("RawData/PID/hTPCVsPt"))->Fill(track.pt(), track.tpcSignal()); + RawData.get(HIST("RawData/PID/hTPCVsEta"))->Fill(eta(trkPx, trkPy, trkPz), track.tpcSignal()); + RawData.get(HIST("RawData/PID/hTPCVsPhi"))->Fill(phi(trkPx, trkPy), track.tpcSignal()); + } + + if (track.hasTOF()) { + RawData.get(HIST("RawData/PID/hTOFVsP"))->Fill(mom(trkPx, trkPy, trkPz), track.tofSignal()); + RawData.get(HIST("RawData/PID/hTOFVsPt"))->Fill(track.pt(), track.tofSignal()); + RawData.get(HIST("RawData/PID/hTOFVsEta"))->Fill(eta(trkPx, trkPy, trkPz), track.tofSignal()); + RawData.get(HIST("RawData/PID/hTOFVsPhi"))->Fill(phi(trkPx, trkPy), track.tofSignal()); + } + } + + int countGT = 0; + int countGTselected = 0; + int countGTel = 0; + int countGTmu = 0; + int countGTp = 0; + std::vector trkIdx; + // loop over tracks with selections + for (auto& track : tracks) { + // select primary vertex contributors + if (track.isPVContributor() != 1) { + return; + } + // select good tracks + if (GoodTrackCuts(track) != 1) { + return; + } + countGT++; + int hypoID = testPIDhypo(track); + if (hypoID == P_ELECTRON || hypoID == P_MUON) { + countGTselected++; + trkIdx.push_back(track.index()); + if (hypoID == P_ELECTRON) { + countGTel++; + } + if (hypoID == P_MUON) { + countGTmu++; + } + if (hypoID == P_PROTON) { + countGTp++; + } + } + } + + Statistics.get(HIST("Statistics/hNumberOfTracks"))->Fill(2., countGT); + Statistics.get(HIST("Statistics/hNumberOfTracks"))->Fill(3., countGTselected); + Statistics.get(HIST("Statistics/hNumberOfTracks"))->Fill(4., countGTel); + Statistics.get(HIST("Statistics/hNumberOfTracks"))->Fill(5., countGTmu); + Statistics.get(HIST("Statistics/hNumberGT"))->Fill(countGT); + Statistics.get(HIST("Statistics/hNumberGTselected"))->Fill(countGTselected); + Statistics.get(HIST("Statistics/hNumberGTel"))->Fill(countGTel); + Statistics.get(HIST("Statistics/hNumberGTmu"))->Fill(countGTmu); + Statistics.get(HIST("Statistics/hNumberGTp"))->Fill(countGTp); + + if (countGT == 2) { + TLorentzVector mother, daughter[2]; + auto trkDaughter1 = tracks.iteratorAt(trkIdx[0]); + auto trkDaughter2 = tracks.iteratorAt(trkIdx[1]); + if ((trkDaughter1.sign() * trkDaughter2.sign()) > 0) { + return; + } + if (countGTel == 2) { + if (!(pow(trkDaughter1.tpcNSigmaEl(), 2) + pow(trkDaughter2.tpcNSigmaEl(), 2) < pow(trkDaughter1.tpcNSigmaMu(), 2) + pow(trkDaughter2.tpcNSigmaMu(), 2))) { + return; + } + auto ene1 = sqrt(pow(trkDaughter1.px(), 2.) + pow(trkDaughter1.py(), 2.) + pow(trkDaughter1.pz(), 2.) + pow(particleMass(pdg, 11), 2.)); + auto ene2 = sqrt(pow(trkDaughter2.px(), 2.) + pow(trkDaughter2.py(), 2.) + pow(trkDaughter2.pz(), 2.) + pow(particleMass(pdg, 11), 2.)); + daughter[0].SetPxPyPzE(trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz(), ene1); + daughter[1].SetPxPyPzE(trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz(), ene2); + mother = daughter[0] + daughter[1]; + + TGel.get(HIST("TGel/hTrackPt"))->Fill(daughter[0].Pt()); + TGel.get(HIST("TGel/hTrackPt"))->Fill(daughter[1].Pt()); + TGel.get(HIST("TGel/hTrackEta"))->Fill(daughter[0].Eta()); + TGel.get(HIST("TGel/hTrackEta"))->Fill(daughter[1].Eta()); + TGel.get(HIST("TGel/hTrackPhi"))->Fill(daughter[0].Phi()); + TGel.get(HIST("TGel/hTrackPhi"))->Fill(daughter[1].Phi()); + + if (CandidateCuts(mother) != 1) { + return; + } + + TGelCand.get(HIST("TGelCand/hTrackPt"))->Fill(daughter[0].Pt()); + TGelCand.get(HIST("TGelCand/hTrackPt"))->Fill(daughter[1].Pt()); + TGelCand.get(HIST("TGelCand/hTrackEta"))->Fill(daughter[0].Eta()); + TGelCand.get(HIST("TGelCand/hTrackEta"))->Fill(daughter[1].Eta()); + TGelCand.get(HIST("TGelCand/hTrackPhi"))->Fill(daughter[0].Phi()); + TGelCand.get(HIST("TGelCand/hTrackPhi"))->Fill(daughter[1].Phi()); + TGelCand.get(HIST("TGelCand/hJpsiPt"))->Fill(mother.Pt()); + if (mother.Pt() < 0.2f) { + // fill track histos + JPsiToEl.get(HIST("JPsiToEl/Coherent/hPt1"))->Fill(daughter[0].Pt()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hPt2"))->Fill(daughter[1].Pt()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hEta1"))->Fill(daughter[0].Eta()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hEta2"))->Fill(daughter[1].Eta()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hPhi1"))->Fill(daughter[0].Phi()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hPhi2"))->Fill(daughter[1].Phi()); + if (trkDaughter1.hasTPC()) { + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsPt"))->Fill(daughter[0].Pt(), trkDaughter1.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsEta"))->Fill(daughter[0].Eta(), trkDaughter1.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsPhi"))->Fill(daughter[0].Phi(), trkDaughter1.tpcSignal()); + } + if (trkDaughter2.hasTPC()) { + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsPt"))->Fill(daughter[1].Pt(), trkDaughter2.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsEta"))->Fill(daughter[1].Eta(), trkDaughter2.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsPhi"))->Fill(daughter[1].Phi(), trkDaughter2.tpcSignal()); + } + + if (trkDaughter1.hasTOF()) { + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsPt"))->Fill(daughter[0].Pt(), trkDaughter1.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsEta"))->Fill(daughter[0].Eta(), trkDaughter1.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsPhi"))->Fill(daughter[0].Phi(), trkDaughter1.tofSignal()); + } + if (trkDaughter2.hasTOF()) { + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsPt"))->Fill(daughter[1].Pt(), trkDaughter2.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsEta"))->Fill(daughter[1].Eta(), trkDaughter2.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsPhi"))->Fill(daughter[1].Phi(), trkDaughter2.tofSignal()); + } + // fill J/psi histos + JPsiToEl.get(HIST("JPsiToEl/Coherent/hPt"))->Fill(mother.Pt()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hEta"))->Fill(mother.Eta()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hPhi"))->Fill(mother.Phi()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hRap"))->Fill(mother.Rapidity()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hIVM"))->Fill(mother.M()); + + float* q = correlation(&daughter[0], &daughter[1], &mother); + Correlation.get(HIST("Correlation/Electron/Coherent/Phi1"))->Fill(daughter[0].Phi(), 1.); + Correlation.get(HIST("Correlation/Electron/Coherent/Phi2"))->Fill(daughter[1].Phi(), 1.); + Correlation.get(HIST("Correlation/Electron/Coherent/Phi"))->Fill(q[1], 1.); + Correlation.get(HIST("Correlation/Electron/Coherent/CosTheta"))->Fill(q[2], 1.); + Correlation.get(HIST("Correlation/Electron/Coherent/AccoplAngle"))->Fill(q[0], 1.); + Correlation.get(HIST("Correlation/Electron/Coherent/CosThetaPhi"))->Fill(q[2], q[1]); + + double dp = DeltaPhi(daughter[0], daughter[1]); + Asymmetry.get(HIST("Asymmetry/Electron/Coherent/DeltaPhi"))->Fill(dp); + } + if (mother.Pt() > 0.2f) { + // fill track histos + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPt1"))->Fill(daughter[0].Pt()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPt2"))->Fill(daughter[1].Pt()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hEta1"))->Fill(daughter[0].Eta()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hEta2"))->Fill(daughter[1].Eta()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPhi1"))->Fill(daughter[0].Phi()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPhi2"))->Fill(daughter[1].Phi()); + if (trkDaughter1.hasTPC()) { + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsPt"))->Fill(daughter[0].Pt(), trkDaughter1.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsEta"))->Fill(daughter[0].Eta(), trkDaughter1.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsPhi"))->Fill(daughter[0].Phi(), trkDaughter1.tpcSignal()); + } + if (trkDaughter2.hasTPC()) { + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsPt"))->Fill(daughter[1].Pt(), trkDaughter2.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsEta"))->Fill(daughter[1].Eta(), trkDaughter2.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsPhi"))->Fill(daughter[1].Phi(), trkDaughter2.tpcSignal()); + } + + if (trkDaughter1.hasTOF()) { + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsPt"))->Fill(daughter[0].Pt(), trkDaughter1.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsEta"))->Fill(daughter[0].Eta(), trkDaughter1.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsPhi"))->Fill(daughter[0].Phi(), trkDaughter1.tofSignal()); + } + if (trkDaughter2.hasTOF()) { + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsPt"))->Fill(daughter[1].Pt(), trkDaughter2.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsEta"))->Fill(daughter[1].Eta(), trkDaughter2.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsPhi"))->Fill(daughter[1].Phi(), trkDaughter2.tofSignal()); + } + // fill J/psi histos + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPt"))->Fill(mother.Pt()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hEta"))->Fill(mother.Eta()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPhi"))->Fill(mother.Phi()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hRap"))->Fill(mother.Rapidity()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hIVM"))->Fill(mother.M()); + + float* q = correlation(&daughter[0], &daughter[1], &mother); + Correlation.get(HIST("Correlation/Electron/Incoherent/Phi1"))->Fill(daughter[0].Phi(), 1.); + Correlation.get(HIST("Correlation/Electron/Incoherent/Phi2"))->Fill(daughter[1].Phi(), 1.); + Correlation.get(HIST("Correlation/Electron/Incoherent/Phi"))->Fill(q[1], 1.); + Correlation.get(HIST("Correlation/Electron/Incoherent/CosTheta"))->Fill(q[2], 1.); + Correlation.get(HIST("Correlation/Electron/Incoherent/AccoplAngle"))->Fill(q[0], 1.); + Correlation.get(HIST("Correlation/Electron/Incoherent/CosThetaPhi"))->Fill(q[2], q[1]); + + double dp = DeltaPhi(daughter[0], daughter[1]); + Asymmetry.get(HIST("Asymmetry/Electron/Incoherent/DeltaPhi"))->Fill(dp); + } + } // end electrons + if (countGTmu == 2) { + TLorentzVector mother, daughter[2]; + auto trkDaughter1 = tracks.iteratorAt(trkIdx[0]); + auto trkDaughter2 = tracks.iteratorAt(trkIdx[1]); + if (!(pow(trkDaughter1.tpcNSigmaMu(), 2) + pow(trkDaughter2.tpcNSigmaMu(), 2) < pow(trkDaughter1.tpcNSigmaEl(), 2) + pow(trkDaughter2.tpcNSigmaEl(), 2))) { + return; + } + auto ene1 = sqrt(pow(trkDaughter1.px(), 2.) + pow(trkDaughter1.py(), 2.) + pow(trkDaughter1.pz(), 2.) + pow(particleMass(pdg, 13), 2.)); + auto ene2 = sqrt(pow(trkDaughter2.px(), 2.) + pow(trkDaughter2.py(), 2.) + pow(trkDaughter2.pz(), 2.) + pow(particleMass(pdg, 13), 2.)); + daughter[0].SetPxPyPzE(trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz(), ene1); + daughter[1].SetPxPyPzE(trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz(), ene2); + mother = daughter[0] + daughter[1]; + + TGmu.get(HIST("TGmu/hTrackPt"))->Fill(daughter[0].Pt()); + TGmu.get(HIST("TGmu/hTrackPt"))->Fill(daughter[1].Pt()); + TGmu.get(HIST("TGmu/hTrackEta"))->Fill(daughter[0].Eta()); + TGmu.get(HIST("TGmu/hTrackEta"))->Fill(daughter[1].Eta()); + TGmu.get(HIST("TGmu/hTrackPhi"))->Fill(daughter[0].Phi()); + TGmu.get(HIST("TGmu/hTrackPhi"))->Fill(daughter[1].Phi()); + + if (CandidateCuts(mother) != 1) { + return; + } + + TGmuCand.get(HIST("TGmuCand/hTrackPt"))->Fill(daughter[0].Pt()); + TGmuCand.get(HIST("TGmuCand/hTrackPt"))->Fill(daughter[1].Pt()); + TGmuCand.get(HIST("TGmuCand/hTrackEta"))->Fill(daughter[0].Eta()); + TGmuCand.get(HIST("TGmuCand/hTrackEta"))->Fill(daughter[1].Eta()); + TGmuCand.get(HIST("TGmuCand/hTrackPhi"))->Fill(daughter[0].Phi()); + TGmuCand.get(HIST("TGmuCand/hTrackPhi"))->Fill(daughter[1].Phi()); + TGmuCand.get(HIST("TGmuCand/hJpsiPt"))->Fill(mother.Pt()); + + if (mother.Pt() < 0.2f) { + // fill track histos + JPsiToMu.get(HIST("JPsiToMu/Coherent/hPt1"))->Fill(daughter[0].Pt()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hPt2"))->Fill(daughter[1].Pt()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hEta1"))->Fill(daughter[0].Eta()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hEta2"))->Fill(daughter[1].Eta()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hPhi1"))->Fill(daughter[0].Phi()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hPhi2"))->Fill(daughter[1].Phi()); + if (trkDaughter1.hasTPC()) { + JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsPt"))->Fill(daughter[0].Pt(), trkDaughter1.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsEta"))->Fill(daughter[0].Eta(), trkDaughter1.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsPhi"))->Fill(daughter[0].Phi(), trkDaughter1.tpcSignal()); + } + if (trkDaughter2.hasTPC()) { + JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsPt"))->Fill(daughter[1].Pt(), trkDaughter2.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsEta"))->Fill(daughter[1].Eta(), trkDaughter2.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsPhi"))->Fill(daughter[1].Phi(), trkDaughter2.tpcSignal()); + } + // fill J/psi histos + JPsiToMu.get(HIST("JPsiToMu/Coherent/hPt"))->Fill(mother.Pt()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hEta"))->Fill(mother.Eta()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hPhi"))->Fill(mother.Phi()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hRap"))->Fill(mother.Rapidity()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hIVM"))->Fill(mother.M()); + + float* q = correlation(&daughter[0], &daughter[1], &mother); + Correlation.get(HIST("Correlation/Muon/Coherent/Phi1"))->Fill(daughter[0].Phi(), 1.); + Correlation.get(HIST("Correlation/Muon/Coherent/Phi2"))->Fill(daughter[1].Phi(), 1.); + Correlation.get(HIST("Correlation/Muon/Coherent/Phi"))->Fill(q[1], 1.); + Correlation.get(HIST("Correlation/Muon/Coherent/CosTheta"))->Fill(q[2], 1.); + Correlation.get(HIST("Correlation/Muon/Coherent/AccoplAngle"))->Fill(q[0], 1.); + Correlation.get(HIST("Correlation/Muon/Coherent/CosThetaPhi"))->Fill(q[2], q[1]); + + double dp = DeltaPhi(daughter[0], daughter[1]); + Asymmetry.get(HIST("Asymmetry/Muon/Coherent/DeltaPhi"))->Fill(dp); + } + if (mother.Pt() > 0.2f) { + // fill track histos + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPt1"))->Fill(daughter[0].Pt()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPt2"))->Fill(daughter[1].Pt()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hEta1"))->Fill(daughter[0].Eta()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hEta2"))->Fill(daughter[1].Eta()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPhi1"))->Fill(daughter[0].Phi()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPhi2"))->Fill(daughter[1].Phi()); + if (trkDaughter1.hasTPC()) { + JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsPt"))->Fill(daughter[0].Pt(), trkDaughter1.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsEta"))->Fill(daughter[0].Eta(), trkDaughter1.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsPhi"))->Fill(daughter[0].Phi(), trkDaughter1.tpcSignal()); + } + if (trkDaughter2.hasTPC()) { + JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsPt"))->Fill(daughter[1].Pt(), trkDaughter2.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsEta"))->Fill(daughter[1].Eta(), trkDaughter2.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsPhi"))->Fill(daughter[1].Phi(), trkDaughter2.tpcSignal()); + } + // fill J/psi histos + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPt"))->Fill(mother.Pt()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hEta"))->Fill(mother.Eta()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPhi"))->Fill(mother.Phi()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hRap"))->Fill(mother.Rapidity()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hIVM"))->Fill(mother.M()); + + float* q = correlation(&daughter[0], &daughter[1], &mother); + Correlation.get(HIST("Correlation/Muon/Incoherent/Phi1"))->Fill(daughter[0].Phi(), 1.); + Correlation.get(HIST("Correlation/Muon/Incoherent/Phi2"))->Fill(daughter[1].Phi(), 1.); + Correlation.get(HIST("Correlation/Muon/Incoherent/Phi"))->Fill(q[1], 1.); + Correlation.get(HIST("Correlation/Muon/Incoherent/CosTheta"))->Fill(q[2], 1.); + Correlation.get(HIST("Correlation/Muon/Incoherent/AccoplAngle"))->Fill(q[0], 1.); + Correlation.get(HIST("Correlation/Muon/Incoherent/CosThetaPhi"))->Fill(q[2], q[1]); + + double dp = DeltaPhi(daughter[0], daughter[1]); + Asymmetry.get(HIST("Asymmetry/Muon/Incoherent/DeltaPhi"))->Fill(dp); + } + } // end muons + if (countGTp == 2) { + TLorentzVector mother, daughter[2]; + auto trkDaughter1 = tracks.iteratorAt(trkIdx[0]); + auto trkDaughter2 = tracks.iteratorAt(trkIdx[1]); + + auto ene1 = sqrt(pow(trkDaughter1.px(), 2.) + pow(trkDaughter1.py(), 2.) + pow(trkDaughter1.pz(), 2.) + pow(particleMass(pdg, 2212), 2.)); + auto ene2 = sqrt(pow(trkDaughter2.px(), 2.) + pow(trkDaughter2.py(), 2.) + pow(trkDaughter2.pz(), 2.) + pow(particleMass(pdg, 2212), 2.)); + daughter[0].SetPxPyPzE(trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz(), ene1); + daughter[1].SetPxPyPzE(trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz(), ene2); + mother = daughter[0] + daughter[1]; + + TGp.get(HIST("TGp/hTrackPt"))->Fill(daughter[0].Pt()); + TGp.get(HIST("TGp/hTrackPt"))->Fill(daughter[1].Pt()); + TGp.get(HIST("TGp/hTrackEta"))->Fill(daughter[0].Eta()); + TGp.get(HIST("TGp/hTrackEta"))->Fill(daughter[1].Eta()); + TGp.get(HIST("TGp/hTrackPhi"))->Fill(daughter[0].Phi()); + TGp.get(HIST("TGp/hTrackPhi"))->Fill(daughter[1].Phi()); + + if (CandidateCuts(mother) != 1) { + return; + } + + TGpCand.get(HIST("TGpCand/hTrackPt"))->Fill(daughter[0].Pt()); + TGpCand.get(HIST("TGpCand/hTrackPt"))->Fill(daughter[1].Pt()); + TGpCand.get(HIST("TGpCand/hTrackEta"))->Fill(daughter[0].Eta()); + TGpCand.get(HIST("TGpCand/hTrackEta"))->Fill(daughter[1].Eta()); + TGpCand.get(HIST("TGpCand/hTrackPhi"))->Fill(daughter[0].Phi()); + TGpCand.get(HIST("TGpCand/hTrackPhi"))->Fill(daughter[1].Phi()); + TGpCand.get(HIST("TGpCand/hJpsiPt"))->Fill(mother.Pt()); + + if (mother.Pt() < 0.2f) { + // fill track histos + JPsiToP.get(HIST("JPsiToP/Coherent/hPt1"))->Fill(daughter[0].Pt()); + JPsiToP.get(HIST("JPsiToP/Coherent/hPt2"))->Fill(daughter[1].Pt()); + JPsiToP.get(HIST("JPsiToP/Coherent/hEta1"))->Fill(daughter[0].Eta()); + JPsiToP.get(HIST("JPsiToP/Coherent/hEta2"))->Fill(daughter[1].Eta()); + JPsiToP.get(HIST("JPsiToP/Coherent/hPhi1"))->Fill(daughter[0].Phi()); + JPsiToP.get(HIST("JPsiToP/Coherent/hPhi2"))->Fill(daughter[1].Phi()); + if (trkDaughter1.hasTPC()) { + JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsPt"))->Fill(daughter[0].Pt(), trkDaughter1.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsEta"))->Fill(daughter[0].Eta(), trkDaughter1.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsPhi"))->Fill(daughter[0].Phi(), trkDaughter1.tpcSignal()); + } + if (trkDaughter2.hasTPC()) { + JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsPt"))->Fill(daughter[1].Pt(), trkDaughter2.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsEta"))->Fill(daughter[1].Eta(), trkDaughter2.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsPhi"))->Fill(daughter[1].Phi(), trkDaughter2.tpcSignal()); + } + // fill J/psi histos + JPsiToP.get(HIST("JPsiToP/Coherent/hPt"))->Fill(mother.Pt()); + JPsiToP.get(HIST("JPsiToP/Coherent/hEta"))->Fill(mother.Eta()); + JPsiToP.get(HIST("JPsiToP/Coherent/hPhi"))->Fill(mother.Phi()); + JPsiToP.get(HIST("JPsiToP/Coherent/hRap"))->Fill(mother.Rapidity()); + JPsiToP.get(HIST("JPsiToP/Coherent/hIVM"))->Fill(mother.M()); + } + if (mother.Pt() > 0.2f) { + // fill track histos + JPsiToP.get(HIST("JPsiToP/Incoherent/hPt1"))->Fill(daughter[0].Pt()); + JPsiToP.get(HIST("JPsiToP/Incoherent/hPt2"))->Fill(daughter[1].Pt()); + JPsiToP.get(HIST("JPsiToP/Incoherent/hEta1"))->Fill(daughter[0].Eta()); + JPsiToP.get(HIST("JPsiToP/Incoherent/hEta2"))->Fill(daughter[1].Eta()); + JPsiToP.get(HIST("JPsiToP/Incoherent/hPhi1"))->Fill(daughter[0].Phi()); + JPsiToP.get(HIST("JPsiToP/Incoherent/hPhi2"))->Fill(daughter[1].Phi()); + if (trkDaughter1.hasTPC()) { + JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsPt"))->Fill(daughter[0].Pt(), trkDaughter1.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsEta"))->Fill(daughter[0].Eta(), trkDaughter1.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsPhi"))->Fill(daughter[0].Phi(), trkDaughter1.tpcSignal()); + } + if (trkDaughter2.hasTPC()) { + JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsPt"))->Fill(daughter[1].Pt(), trkDaughter2.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsEta"))->Fill(daughter[1].Eta(), trkDaughter2.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsPhi"))->Fill(daughter[1].Phi(), trkDaughter2.tpcSignal()); + } + // fill J/psi histos + JPsiToP.get(HIST("JPsiToP/Incoherent/hPt"))->Fill(mother.Pt()); + JPsiToP.get(HIST("JPsiToP/Incoherent/hEta"))->Fill(mother.Eta()); + JPsiToP.get(HIST("JPsiToP/Incoherent/hPhi"))->Fill(mother.Phi()); + JPsiToP.get(HIST("JPsiToP/Incoherent/hRap"))->Fill(mother.Rapidity()); + JPsiToP.get(HIST("JPsiToP/Incoherent/hIVM"))->Fill(mother.M()); + } + } // end protons + } // end two tracks + } // end process +}; // end struct + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"upc-jpsi-corr"}), + }; +} From 72e632d7677ab9b45c8cb908fd992af44579f588 Mon Sep 17 00:00:00 2001 From: Sara Haidlova Date: Wed, 13 Mar 2024 11:20:35 +0100 Subject: [PATCH 2/6] PWGUD: New central barrel task, additional changes --- PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h | 29 +++++++++------------ PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx | 15 +++++++---- 2 files changed, 22 insertions(+), 22 deletions(-) diff --git a/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h b/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h index a8d0fa2af0a..18c18f7095a 100644 --- a/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h +++ b/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h @@ -11,10 +11,13 @@ /// /// \brief /// \author Sara Haidlova, sara.haidlova@cern.ch +/// \since March 2024 #ifndef PWGUD_CORE_UPCJPSICENTRALBARRELCORR_H_ #define PWGUD_CORE_UPCJPSICENTRALBARRELCORR_H_ +#include + using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; @@ -37,8 +40,9 @@ int testPIDhypo(T trackPID) std::min_element(std::begin(nSigmaTPC), std::end(nSigmaTPC))); if (trackPID.hasTPC()) { return enumChoiceTPC; - } else - return -1; + } else { + return -1; + } } float mom(float px, float py, float pz) @@ -50,8 +54,9 @@ float phi(float px, float py) { if (px != 0) { return std::atan2(py, px); - } else - return -999; + } else { + return -999; + } } float eta(float px, float py, float pz) @@ -59,21 +64,11 @@ float eta(float px, float py, float pz) float mom = std::sqrt(px * px + py * py + pz * pz); if (mom != 0) { return std::atanh(pz / mom); - } else - return -999; + } else { + return -999; + } } -/*float energy(float mass, float px, float py, float pz) -// Just a simple function to return track energy -{ - return std::sqrt(mass * mass + px * px + py * py + pz * pz); -} - -float rap(float mass, float px, float py, float pz) -// Just a simple function to return track rapidity -{ - return 0.5 * std::log((energy(mass, px, py, pz) + pz) / (energy(mass, px, py, pz) - pz)); -}*/ float* correlation(TLorentzVector* lv1, TLorentzVector* lv2, TLorentzVector* lv) { diff --git a/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx b/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx index 4ff739baaf0..483518b2d7d 100644 --- a/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx +++ b/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx @@ -12,19 +12,24 @@ /// /// \brief /// \author Sara Haidlova, sara.haidlova@cern.ch +/// \since March 2024 +//O2 headers #include "Framework/runDataProcessing.h" #include "Framework/AnalysisTask.h" +#include "CommonConstants/MathConstants.h" +//O2Physics headers +#include "Common/DataModel/PIDResponse.h" +#include "PWGUD/DataModel/UDTables.h" +#include "PWGUD/Core/UDHelpers.h" +#include "PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h" + +//ROOT headers #include "TDatabasePDG.h" #include "TLorentzVector.h" #include "TMath.h" #include "TVector3.h" -#include "Common/DataModel/PIDResponse.h" -#include "PWGUD/DataModel/UDTables.h" -#include "PWGUD/Core/UDHelpers.h" -#include -#include "UPCJpsiCentralBarrelCorrHelper.h" using namespace o2; using namespace o2::framework; From 11c905a12e65ce8b3029c79eb0afdce188b46157 Mon Sep 17 00:00:00 2001 From: Sara Haidlova Date: Wed, 13 Mar 2024 12:27:09 +0100 Subject: [PATCH 3/6] PWGUD: New JPsi task, clang corrections --- PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h | 13 ++++++------- PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx | 6 +++--- 2 files changed, 9 insertions(+), 10 deletions(-) diff --git a/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h b/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h index 18c18f7095a..8f1d9a4c691 100644 --- a/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h +++ b/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h @@ -13,8 +13,8 @@ /// \author Sara Haidlova, sara.haidlova@cern.ch /// \since March 2024 -#ifndef PWGUD_CORE_UPCJPSICENTRALBARRELCORR_H_ -#define PWGUD_CORE_UPCJPSICENTRALBARRELCORR_H_ +#ifndef PWGUD_CORE_UPCJPSICENTRALBARRELCORRHELPER_H_ +#define PWGUD_CORE_UPCJPSICENTRALBARRELCORRHELPER_H_ #include @@ -41,7 +41,7 @@ int testPIDhypo(T trackPID) if (trackPID.hasTPC()) { return enumChoiceTPC; } else { - return -1; + return -1; } } @@ -55,7 +55,7 @@ float phi(float px, float py) if (px != 0) { return std::atan2(py, px); } else { - return -999; + return -999; } } @@ -65,11 +65,10 @@ float eta(float px, float py, float pz) if (mom != 0) { return std::atanh(pz / mom); } else { - return -999; + return -999; } } - float* correlation(TLorentzVector* lv1, TLorentzVector* lv2, TLorentzVector* lv) { TLorentzVector pa(1., 0., 0, 1); // projectile @@ -213,4 +212,4 @@ double DeltaPhi(TLorentzVector lv1, TLorentzVector lv2) return dp; } -#endif // PWGUD_CORE_UPCJPSICENTRALBARRELCORR_H_ +#endif // PWGUD_CORE_UPCJPSICENTRALBARRELCORRHELPER_H_ diff --git a/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx b/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx index 483518b2d7d..5cd85494502 100644 --- a/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx +++ b/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx @@ -14,18 +14,18 @@ /// \author Sara Haidlova, sara.haidlova@cern.ch /// \since March 2024 -//O2 headers +// O2 headers #include "Framework/runDataProcessing.h" #include "Framework/AnalysisTask.h" #include "CommonConstants/MathConstants.h" -//O2Physics headers +// O2Physics headers #include "Common/DataModel/PIDResponse.h" #include "PWGUD/DataModel/UDTables.h" #include "PWGUD/Core/UDHelpers.h" #include "PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h" -//ROOT headers +// ROOT headers #include "TDatabasePDG.h" #include "TLorentzVector.h" #include "TMath.h" From a52abf63b604674ee05149100b9f48bf2d27b628 Mon Sep 17 00:00:00 2001 From: Sara Haidlova Date: Thu, 21 Mar 2024 16:02:44 +0100 Subject: [PATCH 4/6] PWGUD: additional code changes including RecoDecay.h --- PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h | 30 +- PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx | 451 ++++++++++---------- 2 files changed, 239 insertions(+), 242 deletions(-) diff --git a/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h b/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h index 8f1d9a4c691..15ea539ebf9 100644 --- a/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h +++ b/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h @@ -17,6 +17,7 @@ #define PWGUD_CORE_UPCJPSICENTRALBARRELCORRHELPER_H_ #include +#include "CommonConstants/MathConstants.h" using namespace o2; using namespace o2::framework; @@ -45,30 +46,6 @@ int testPIDhypo(T trackPID) } } -float mom(float px, float py, float pz) -{ - return std::sqrt(px * px + py * py + pz * pz); -} - -float phi(float px, float py) -{ - if (px != 0) { - return std::atan2(py, px); - } else { - return -999; - } -} - -float eta(float px, float py, float pz) -{ - float mom = std::sqrt(px * px + py * py + pz * pz); - if (mom != 0) { - return std::atanh(pz / mom); - } else { - return -999; - } -} - float* correlation(TLorentzVector* lv1, TLorentzVector* lv2, TLorentzVector* lv) { TLorentzVector pa(1., 0., 0, 1); // projectile @@ -79,7 +56,7 @@ float* correlation(TLorentzVector* lv1, TLorentzVector* lv2, TLorentzVector* lv) // Accoplanarity angle Float_t deltaPhi; deltaPhi = lv1->Phi() - lv2->Phi(); - float accoplCut = 1. - TMath::Abs(deltaPhi) / TMath::Pi(); + float accoplCut = 1. - abs(deltaPhi) / o2::constants::math::PI; // z TLorentzVector z; Float_t part1, part2; @@ -193,7 +170,8 @@ float* correlation(TLorentzVector* lv1, TLorentzVector* lv2, TLorentzVector* lv) float phi = atan2(SinThetaSinPhiCS, SinThetaCosPhiCS); // if (phi>=0) phi = phi; - // if (phi<0) phi = phi + 2*TMath::Pi(); + if (phi < 0) + phi = phi + o2::constants::math::TwoPI; q[0] = accoplCut; q[1] = phi; diff --git a/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx b/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx index 5cd85494502..74fa86f360b 100644 --- a/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx +++ b/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx @@ -1,4 +1,3 @@ - // Copyright 2019-2020 CERN and copyright holders of ALICE O2. // See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. // All rights not expressly granted are reserved. @@ -21,6 +20,7 @@ // O2Physics headers #include "Common/DataModel/PIDResponse.h" +#include "Common/Core/RecoDecay.h" #include "PWGUD/DataModel/UDTables.h" #include "PWGUD/Core/UDHelpers.h" #include "PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h" @@ -28,8 +28,6 @@ // ROOT headers #include "TDatabasePDG.h" #include "TLorentzVector.h" -#include "TMath.h" -#include "TVector3.h" using namespace o2; using namespace o2::framework; @@ -43,7 +41,7 @@ struct upcJpsiCentralBarrel { ConfigurableAxis pAxis{"pAxis", {250.0f, 0.1f, 3.0f}, "#it{p} (GeV/#it{c})"}; ConfigurableAxis etaAxis{"etaAxis", {250.0f, -1.5f, 1.5f}, "#eta (-)"}; ConfigurableAxis countAxis{"countAxis", {10.0f, 0.0f, 10.0f}, "Number of events (-)"}; - ConfigurableAxis phiAxis{"phiAxis", {250.0f, -PI, PI}, "#phi (rad)"}; + ConfigurableAxis phiAxis{"phiAxis", {250.0f, 0, TwoPI}, "#phi (rad)"}; ConfigurableAxis accoplAxis{"accoplAxis", {250.0f, -0.2f, 0.2f}, "accAngle"}; ConfigurableAxis thetaAxis{"thetaAxis", {250.0f, -1.5f, 1.5f}, "cos #theta (-)"}; ConfigurableAxis sigTPCAxis{"sigTPCAxis", {100.0f, 0, 200.0f}, "TPC d#it{E}/d#it{x}"}; @@ -361,10 +359,10 @@ struct upcJpsiCentralBarrel { Correlation.add("Correlation/Electron/Incoherent/CosThetaPhi", "CosThetaPhi", {HistType::kTH2F, {{axisAngTheta}, {axisPhi}}}); // Asymmetry histograms - Asymmetry.add("Asymmetry/Muon/Coherent/DeltaPhi", "DeltaPhi", {HistType::kTH1F, {axisPhi}}); - Asymmetry.add("Asymmetry/Muon/Incoherent/DeltaPhi", "DeltaPhi", {HistType::kTH1F, {axisPhi}}); - Asymmetry.add("Asymmetry/Electron/Coherent/DeltaPhi", "DeltaPhi", {HistType::kTH1F, {axisPhi}}); - Asymmetry.add("Asymmetry/Electron/Incoherent/DeltaPhi", "DeltaPhi", {HistType::kTH1F, {axisPhi}}); + Asymmetry.add("Asymmetry/Muon/Coherent/DeltaPhi", "DeltaPhi", {HistType::kTH1F, {{180, -PI, PI}}}); + Asymmetry.add("Asymmetry/Muon/Incoherent/DeltaPhi", "DeltaPhi", {HistType::kTH1F, {{180, -PI, PI}}}); + Asymmetry.add("Asymmetry/Electron/Coherent/DeltaPhi", "DeltaPhi", {HistType::kTH1F, {{180, -PI, PI}}}); + Asymmetry.add("Asymmetry/Electron/Incoherent/DeltaPhi", "DeltaPhi", {HistType::kTH1F, {{180, -PI, PI}}}); } float particleMass(TDatabasePDG* pdg, int pid) @@ -381,7 +379,7 @@ struct upcJpsiCentralBarrel { bool GoodTrackCuts(T const& track) { // kinematics - if (eta(track.px(), track.py(), track.pz()) < -EtaCut || eta(track.px(), track.py(), track.pz()) > EtaCut) { + if (RecoDecay::eta(std::array{track.px(), track.py(), track.pz()}) < -EtaCut || RecoDecay::eta(std::array{track.px(), track.py(), track.pz()}) > EtaCut) { return false; } // DCA @@ -418,11 +416,12 @@ struct upcJpsiCentralBarrel { template bool CandidateCuts(C const& mother) { - if (abs(mother.Rapidity()) > 0.8f) { + if (abs(RecoDecay::y(std::array{mother[0], mother[1], mother[2]}, mother[3])) > 0.8f) { return false; } - if (mother.M() < 2.5f) { + double energy = RecoDecay::e(mother[0], mother[1], mother[2], mother[3]); + if (RecoDecay::m(std::array{mother[0], mother[1], mother[2]}, energy) < 2.5f) { return false; } @@ -443,40 +442,40 @@ struct upcJpsiCentralBarrel { if (track.isPVContributor() == 1) { Statistics.get(HIST("Statistics/hNumberOfTracks"))->Fill(1); PVContributors.get(HIST("PVContributors/hTrackPt"))->Fill(track.pt()); - PVContributors.get(HIST("PVContributors/hTrackEta"))->Fill(eta(trkPx, trkPy, trkPz)); - PVContributors.get(HIST("PVContributors/hTrackPhi"))->Fill(phi(trkPx, trkPy)); + PVContributors.get(HIST("PVContributors/hTrackEta"))->Fill(RecoDecay::eta(std::array{trkPx, trkPy, trkPz})); + PVContributors.get(HIST("PVContributors/hTrackPhi"))->Fill(RecoDecay::phi(trkPx, trkPy)); if (track.hasTPC()) { - PVContributors.get(HIST("PVContributors/PID/hTPCVsP"))->Fill(mom(trkPx, trkPy, trkPz), track.tpcSignal()); + PVContributors.get(HIST("PVContributors/PID/hTPCVsP"))->Fill(RecoDecay::sqrtSumOfSquares(trkPx, trkPy, trkPz), track.tpcSignal()); PVContributors.get(HIST("PVContributors/PID/hTPCVsPt"))->Fill(track.pt(), track.tpcSignal()); - PVContributors.get(HIST("PVContributors/PID/hTPCVsEta"))->Fill(eta(trkPx, trkPy, trkPz), track.tpcSignal()); - PVContributors.get(HIST("PVContributors/PID/hTPCVsPhi"))->Fill(phi(trkPx, trkPy), track.tpcSignal()); + PVContributors.get(HIST("PVContributors/PID/hTPCVsEta"))->Fill(RecoDecay::eta(std::array{trkPx, trkPy, trkPz}), track.tpcSignal()); + PVContributors.get(HIST("PVContributors/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(trkPx, trkPy), track.tpcSignal()); } if (track.hasTOF()) { - PVContributors.get(HIST("PVContributors/PID/hTOFVsP"))->Fill(mom(trkPx, trkPy, trkPz), track.tofSignal()); + PVContributors.get(HIST("PVContributors/PID/hTOFVsP"))->Fill(RecoDecay::sqrtSumOfSquares(trkPx, trkPy, trkPz), track.tofSignal()); PVContributors.get(HIST("PVContributors/PID/hTOFVsPt"))->Fill(track.pt(), track.tofSignal()); - PVContributors.get(HIST("PVContributors/PID/hTOFVsEta"))->Fill(eta(trkPx, trkPy, trkPz), track.tofSignal()); - PVContributors.get(HIST("PVContributors/PID/hTOFVsPhi"))->Fill(phi(trkPx, trkPy), track.tofSignal()); + PVContributors.get(HIST("PVContributors/PID/hTOFVsEta"))->Fill(RecoDecay::eta(std::array{trkPx, trkPy, trkPz}), track.tofSignal()); + PVContributors.get(HIST("PVContributors/PID/hTOFVsPhi"))->Fill(RecoDecay::phi(trkPx, trkPy), track.tofSignal()); } } RawData.get(HIST("RawData/hTrackPt"))->Fill(track.pt()); - RawData.get(HIST("RawData/hTrackEta"))->Fill(eta(trkPx, trkPy, trkPz)); - RawData.get(HIST("RawData/hTrackPhi"))->Fill(phi(trkPx, trkPy)); + RawData.get(HIST("RawData/hTrackEta"))->Fill(RecoDecay::eta(std::array{trkPx, trkPy, trkPz})); + RawData.get(HIST("RawData/hTrackPhi"))->Fill(RecoDecay::phi(trkPx, trkPy)); if (track.hasTPC()) { - RawData.get(HIST("RawData/PID/hTPCVsP"))->Fill(mom(trkPx, trkPy, trkPz), track.tpcSignal()); + RawData.get(HIST("RawData/PID/hTPCVsP"))->Fill(RecoDecay::sqrtSumOfSquares(trkPx, trkPy, trkPz), track.tpcSignal()); RawData.get(HIST("RawData/PID/hTPCVsPt"))->Fill(track.pt(), track.tpcSignal()); - RawData.get(HIST("RawData/PID/hTPCVsEta"))->Fill(eta(trkPx, trkPy, trkPz), track.tpcSignal()); - RawData.get(HIST("RawData/PID/hTPCVsPhi"))->Fill(phi(trkPx, trkPy), track.tpcSignal()); + RawData.get(HIST("RawData/PID/hTPCVsEta"))->Fill(RecoDecay::eta(std::array{trkPx, trkPy, trkPz}), track.tpcSignal()); + RawData.get(HIST("RawData/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(trkPx, trkPy), track.tpcSignal()); } if (track.hasTOF()) { - RawData.get(HIST("RawData/PID/hTOFVsP"))->Fill(mom(trkPx, trkPy, trkPz), track.tofSignal()); + RawData.get(HIST("RawData/PID/hTOFVsP"))->Fill(RecoDecay::sqrtSumOfSquares(trkPx, trkPy, trkPz), track.tofSignal()); RawData.get(HIST("RawData/PID/hTOFVsPt"))->Fill(track.pt(), track.tofSignal()); - RawData.get(HIST("RawData/PID/hTOFVsEta"))->Fill(eta(trkPx, trkPy, trkPz), track.tofSignal()); - RawData.get(HIST("RawData/PID/hTOFVsPhi"))->Fill(phi(trkPx, trkPy), track.tofSignal()); + RawData.get(HIST("RawData/PID/hTOFVsEta"))->Fill(RecoDecay::eta(std::array{trkPx, trkPy, trkPz}), track.tofSignal()); + RawData.get(HIST("RawData/PID/hTOFVsPhi"))->Fill(RecoDecay::phi(trkPx, trkPy), track.tofSignal()); } } @@ -524,79 +523,85 @@ struct upcJpsiCentralBarrel { Statistics.get(HIST("Statistics/hNumberGTp"))->Fill(countGTp); if (countGT == 2) { - TLorentzVector mother, daughter[2]; + TLorentzVector mom, daughter[2]; auto trkDaughter1 = tracks.iteratorAt(trkIdx[0]); auto trkDaughter2 = tracks.iteratorAt(trkIdx[1]); if ((trkDaughter1.sign() * trkDaughter2.sign()) > 0) { return; } if (countGTel == 2) { - if (!(pow(trkDaughter1.tpcNSigmaEl(), 2) + pow(trkDaughter2.tpcNSigmaEl(), 2) < pow(trkDaughter1.tpcNSigmaMu(), 2) + pow(trkDaughter2.tpcNSigmaMu(), 2))) { + if (!(RecoDecay::sumOfSquares(trkDaughter1.tpcNSigmaMu(), trkDaughter2.tpcNSigmaMu()) < RecoDecay::sumOfSquares(trkDaughter1.tpcNSigmaEl(), trkDaughter2.tpcNSigmaEl()))) { return; } - auto ene1 = sqrt(pow(trkDaughter1.px(), 2.) + pow(trkDaughter1.py(), 2.) + pow(trkDaughter1.pz(), 2.) + pow(particleMass(pdg, 11), 2.)); - auto ene2 = sqrt(pow(trkDaughter2.px(), 2.) + pow(trkDaughter2.py(), 2.) + pow(trkDaughter2.pz(), 2.) + pow(particleMass(pdg, 11), 2.)); + auto ene1 = RecoDecay::e(trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz(), particleMass(pdg, 11)); + auto ene2 = RecoDecay::e(trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz(), particleMass(pdg, 11)); daughter[0].SetPxPyPzE(trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz(), ene1); daughter[1].SetPxPyPzE(trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz(), ene2); - mother = daughter[0] + daughter[1]; + mom = daughter[0] + daughter[1]; + + std::array daughter1 = {trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz()}; + std::array daughter2 = {trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz()}; + + std::array mother = {trkDaughter1.px() + trkDaughter2.px(), trkDaughter1.py() + trkDaughter2.py(), trkDaughter1.pz() + trkDaughter2.pz(), 2 * particleMass(pdg, 11)}; + double energy = RecoDecay::e(mother[0], mother[1], mother[2], mother[3]); - TGel.get(HIST("TGel/hTrackPt"))->Fill(daughter[0].Pt()); - TGel.get(HIST("TGel/hTrackPt"))->Fill(daughter[1].Pt()); - TGel.get(HIST("TGel/hTrackEta"))->Fill(daughter[0].Eta()); - TGel.get(HIST("TGel/hTrackEta"))->Fill(daughter[1].Eta()); - TGel.get(HIST("TGel/hTrackPhi"))->Fill(daughter[0].Phi()); - TGel.get(HIST("TGel/hTrackPhi"))->Fill(daughter[1].Phi()); + TGel.get(HIST("TGel/hTrackPt"))->Fill(trkDaughter1.pt()); + TGel.get(HIST("TGel/hTrackPt"))->Fill(trkDaughter2.pt()); + TGel.get(HIST("TGel/hTrackEta"))->Fill(RecoDecay::eta(daughter1)); + TGel.get(HIST("TGel/hTrackEta"))->Fill(RecoDecay::eta(daughter2)); + TGel.get(HIST("TGel/hTrackPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); + TGel.get(HIST("TGel/hTrackPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); if (CandidateCuts(mother) != 1) { return; } - TGelCand.get(HIST("TGelCand/hTrackPt"))->Fill(daughter[0].Pt()); - TGelCand.get(HIST("TGelCand/hTrackPt"))->Fill(daughter[1].Pt()); - TGelCand.get(HIST("TGelCand/hTrackEta"))->Fill(daughter[0].Eta()); - TGelCand.get(HIST("TGelCand/hTrackEta"))->Fill(daughter[1].Eta()); - TGelCand.get(HIST("TGelCand/hTrackPhi"))->Fill(daughter[0].Phi()); - TGelCand.get(HIST("TGelCand/hTrackPhi"))->Fill(daughter[1].Phi()); - TGelCand.get(HIST("TGelCand/hJpsiPt"))->Fill(mother.Pt()); - if (mother.Pt() < 0.2f) { + TGelCand.get(HIST("TGelCand/hTrackPt"))->Fill(trkDaughter1.pt()); + TGelCand.get(HIST("TGelCand/hTrackPt"))->Fill(trkDaughter2.pt()); + TGelCand.get(HIST("TGelCand/hTrackEta"))->Fill(RecoDecay::eta(daughter1)); + TGelCand.get(HIST("TGelCand/hTrackEta"))->Fill(RecoDecay::eta(daughter2)); + TGelCand.get(HIST("TGelCand/hTrackPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); + TGelCand.get(HIST("TGelCand/hTrackPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); + TGelCand.get(HIST("TGelCand/hJpsiPt"))->Fill(RecoDecay::pt(mother)); + if (RecoDecay::pt(mother) < 0.2f) { // fill track histos - JPsiToEl.get(HIST("JPsiToEl/Coherent/hPt1"))->Fill(daughter[0].Pt()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/hPt2"))->Fill(daughter[1].Pt()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/hEta1"))->Fill(daughter[0].Eta()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/hEta2"))->Fill(daughter[1].Eta()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/hPhi1"))->Fill(daughter[0].Phi()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/hPhi2"))->Fill(daughter[1].Phi()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hPt1"))->Fill(trkDaughter1.pt()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hPt2"))->Fill(trkDaughter2.pt()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hEta1"))->Fill(RecoDecay::eta(daughter1)); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hEta2"))->Fill(RecoDecay::eta(daughter2)); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hPhi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hPhi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); if (trkDaughter1.hasTPC()) { - JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsPt"))->Fill(daughter[0].Pt(), trkDaughter1.tpcSignal()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsEta"))->Fill(daughter[0].Eta(), trkDaughter1.tpcSignal()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsPhi"))->Fill(daughter[0].Phi(), trkDaughter1.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsPt"))->Fill(trkDaughter1.pt(), trkDaughter1.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter1), trkDaughter1.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), trkDaughter1.tpcSignal()); } if (trkDaughter2.hasTPC()) { - JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsPt"))->Fill(daughter[1].Pt(), trkDaughter2.tpcSignal()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsEta"))->Fill(daughter[1].Eta(), trkDaughter2.tpcSignal()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsPhi"))->Fill(daughter[1].Phi(), trkDaughter2.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsPt"))->Fill(trkDaughter2.pt(), trkDaughter2.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter2), trkDaughter2.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), trkDaughter2.tpcSignal()); } if (trkDaughter1.hasTOF()) { - JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsPt"))->Fill(daughter[0].Pt(), trkDaughter1.tofSignal()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsEta"))->Fill(daughter[0].Eta(), trkDaughter1.tofSignal()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsPhi"))->Fill(daughter[0].Phi(), trkDaughter1.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsPt"))->Fill(trkDaughter1.pt(), trkDaughter1.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsEta"))->Fill(RecoDecay::eta(daughter1), trkDaughter1.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), trkDaughter1.tofSignal()); } if (trkDaughter2.hasTOF()) { - JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsPt"))->Fill(daughter[1].Pt(), trkDaughter2.tofSignal()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsEta"))->Fill(daughter[1].Eta(), trkDaughter2.tofSignal()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsPhi"))->Fill(daughter[1].Phi(), trkDaughter2.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsPt"))->Fill(trkDaughter2.pt(), trkDaughter2.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsEta"))->Fill(RecoDecay::eta(daughter2), trkDaughter2.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), trkDaughter2.tofSignal()); } // fill J/psi histos - JPsiToEl.get(HIST("JPsiToEl/Coherent/hPt"))->Fill(mother.Pt()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/hEta"))->Fill(mother.Eta()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/hPhi"))->Fill(mother.Phi()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/hRap"))->Fill(mother.Rapidity()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/hIVM"))->Fill(mother.M()); - - float* q = correlation(&daughter[0], &daughter[1], &mother); - Correlation.get(HIST("Correlation/Electron/Coherent/Phi1"))->Fill(daughter[0].Phi(), 1.); - Correlation.get(HIST("Correlation/Electron/Coherent/Phi2"))->Fill(daughter[1].Phi(), 1.); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hPt"))->Fill(RecoDecay::pt(mother)); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hEta"))->Fill(RecoDecay::eta(std::array{mother[0], mother[1], mother[2]})); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hPhi"))->Fill(RecoDecay::phi(mother[0], mother[1])); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hRap"))->Fill(RecoDecay::y(std::array{mother[0], mother[1], mother[2]}, mother[3])); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hIVM"))->Fill(RecoDecay::m(std::array{mother[0], mother[1], mother[2]}, energy)); + + float* q = correlation(&daughter[0], &daughter[1], &mom); + Correlation.get(HIST("Correlation/Electron/Coherent/Phi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), 1.); + Correlation.get(HIST("Correlation/Electron/Coherent/Phi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), 1.); Correlation.get(HIST("Correlation/Electron/Coherent/Phi"))->Fill(q[1], 1.); Correlation.get(HIST("Correlation/Electron/Coherent/CosTheta"))->Fill(q[2], 1.); Correlation.get(HIST("Correlation/Electron/Coherent/AccoplAngle"))->Fill(q[0], 1.); @@ -605,45 +610,45 @@ struct upcJpsiCentralBarrel { double dp = DeltaPhi(daughter[0], daughter[1]); Asymmetry.get(HIST("Asymmetry/Electron/Coherent/DeltaPhi"))->Fill(dp); } - if (mother.Pt() > 0.2f) { + if (RecoDecay::pt(mother) > 0.2f) { // fill track histos - JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPt1"))->Fill(daughter[0].Pt()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPt2"))->Fill(daughter[1].Pt()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/hEta1"))->Fill(daughter[0].Eta()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/hEta2"))->Fill(daughter[1].Eta()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPhi1"))->Fill(daughter[0].Phi()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPhi2"))->Fill(daughter[1].Phi()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPt1"))->Fill(trkDaughter1.pt()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPt2"))->Fill(trkDaughter2.pt()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hEta1"))->Fill(RecoDecay::eta(daughter1)); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hEta2"))->Fill(RecoDecay::eta(daughter2)); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPhi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPhi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); if (trkDaughter1.hasTPC()) { - JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsPt"))->Fill(daughter[0].Pt(), trkDaughter1.tpcSignal()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsEta"))->Fill(daughter[0].Eta(), trkDaughter1.tpcSignal()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsPhi"))->Fill(daughter[0].Phi(), trkDaughter1.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsPt"))->Fill(trkDaughter1.pt(), trkDaughter1.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter1), trkDaughter1.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), trkDaughter1.tpcSignal()); } if (trkDaughter2.hasTPC()) { - JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsPt"))->Fill(daughter[1].Pt(), trkDaughter2.tpcSignal()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsEta"))->Fill(daughter[1].Eta(), trkDaughter2.tpcSignal()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsPhi"))->Fill(daughter[1].Phi(), trkDaughter2.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsPt"))->Fill(trkDaughter2.pt(), trkDaughter2.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter2), trkDaughter2.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), trkDaughter2.tpcSignal()); } if (trkDaughter1.hasTOF()) { - JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsPt"))->Fill(daughter[0].Pt(), trkDaughter1.tofSignal()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsEta"))->Fill(daughter[0].Eta(), trkDaughter1.tofSignal()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsPhi"))->Fill(daughter[0].Phi(), trkDaughter1.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsPt"))->Fill(trkDaughter1.pt(), trkDaughter1.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsEta"))->Fill(RecoDecay::eta(daughter1), trkDaughter1.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), trkDaughter1.tofSignal()); } if (trkDaughter2.hasTOF()) { - JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsPt"))->Fill(daughter[1].Pt(), trkDaughter2.tofSignal()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsEta"))->Fill(daughter[1].Eta(), trkDaughter2.tofSignal()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsPhi"))->Fill(daughter[1].Phi(), trkDaughter2.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsPt"))->Fill(trkDaughter2.pt(), trkDaughter2.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsEta"))->Fill(RecoDecay::eta(daughter2), trkDaughter2.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), trkDaughter2.tofSignal()); } // fill J/psi histos - JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPt"))->Fill(mother.Pt()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/hEta"))->Fill(mother.Eta()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPhi"))->Fill(mother.Phi()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/hRap"))->Fill(mother.Rapidity()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/hIVM"))->Fill(mother.M()); - - float* q = correlation(&daughter[0], &daughter[1], &mother); - Correlation.get(HIST("Correlation/Electron/Incoherent/Phi1"))->Fill(daughter[0].Phi(), 1.); - Correlation.get(HIST("Correlation/Electron/Incoherent/Phi2"))->Fill(daughter[1].Phi(), 1.); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPt"))->Fill(RecoDecay::pt(mother)); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hEta"))->Fill(RecoDecay::eta(std::array{mother[0], mother[1], mother[2]})); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPhi"))->Fill(RecoDecay::phi(mother[0], mother[1])); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hRap"))->Fill(RecoDecay::y(std::array{mother[0], mother[1], mother[2]}, mother[3])); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hIVM"))->Fill(RecoDecay::m(std::array{mother[0], mother[1], mother[2]}, energy)); + + float* q = correlation(&daughter[0], &daughter[1], &mom); + Correlation.get(HIST("Correlation/Electron/Incoherent/Phi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), 1.); + Correlation.get(HIST("Correlation/Electron/Incoherent/Phi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), 1.); Correlation.get(HIST("Correlation/Electron/Incoherent/Phi"))->Fill(q[1], 1.); Correlation.get(HIST("Correlation/Electron/Incoherent/CosTheta"))->Fill(q[2], 1.); Correlation.get(HIST("Correlation/Electron/Incoherent/AccoplAngle"))->Fill(q[0], 1.); @@ -654,65 +659,72 @@ struct upcJpsiCentralBarrel { } } // end electrons if (countGTmu == 2) { - TLorentzVector mother, daughter[2]; + TLorentzVector mom, daughter[2]; auto trkDaughter1 = tracks.iteratorAt(trkIdx[0]); auto trkDaughter2 = tracks.iteratorAt(trkIdx[1]); - if (!(pow(trkDaughter1.tpcNSigmaMu(), 2) + pow(trkDaughter2.tpcNSigmaMu(), 2) < pow(trkDaughter1.tpcNSigmaEl(), 2) + pow(trkDaughter2.tpcNSigmaEl(), 2))) { + if (!(RecoDecay::sumOfSquares(trkDaughter1.tpcNSigmaEl(), trkDaughter2.tpcNSigmaEl() < RecoDecay::sumOfSquares(trkDaughter1.tpcNSigmaMu(), trkDaughter2.tpcNSigmaMu())))) { return; } - auto ene1 = sqrt(pow(trkDaughter1.px(), 2.) + pow(trkDaughter1.py(), 2.) + pow(trkDaughter1.pz(), 2.) + pow(particleMass(pdg, 13), 2.)); - auto ene2 = sqrt(pow(trkDaughter2.px(), 2.) + pow(trkDaughter2.py(), 2.) + pow(trkDaughter2.pz(), 2.) + pow(particleMass(pdg, 13), 2.)); + auto ene1 = RecoDecay::e(trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz(), particleMass(pdg, 11)); + auto ene2 = RecoDecay::e(trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz(), particleMass(pdg, 11)); daughter[0].SetPxPyPzE(trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz(), ene1); daughter[1].SetPxPyPzE(trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz(), ene2); - mother = daughter[0] + daughter[1]; + mom = daughter[0] + daughter[1]; + + std::array daughter1 = {trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz()}; + std::array daughter2 = {trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz()}; + + std::array mother = {trkDaughter1.px() + trkDaughter2.px(), trkDaughter1.py() + trkDaughter2.py(), trkDaughter1.pz() + trkDaughter2.pz(), 2 * particleMass(pdg, 11)}; - TGmu.get(HIST("TGmu/hTrackPt"))->Fill(daughter[0].Pt()); - TGmu.get(HIST("TGmu/hTrackPt"))->Fill(daughter[1].Pt()); - TGmu.get(HIST("TGmu/hTrackEta"))->Fill(daughter[0].Eta()); - TGmu.get(HIST("TGmu/hTrackEta"))->Fill(daughter[1].Eta()); - TGmu.get(HIST("TGmu/hTrackPhi"))->Fill(daughter[0].Phi()); - TGmu.get(HIST("TGmu/hTrackPhi"))->Fill(daughter[1].Phi()); + double energy = RecoDecay::e(mother[0], mother[1], mother[2], mother[3]); + + TGmu.get(HIST("TGmu/hTrackPt"))->Fill(trkDaughter1.pt()); + TGmu.get(HIST("TGmu/hTrackPt"))->Fill(trkDaughter2.pt()); + TGmu.get(HIST("TGmu/hTrackEta"))->Fill(RecoDecay::eta(daughter1)); + TGmu.get(HIST("TGmu/hTrackEta"))->Fill(RecoDecay::eta(daughter2)); + TGmu.get(HIST("TGmu/hTrackPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); + TGmu.get(HIST("TGmu/hTrackPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); if (CandidateCuts(mother) != 1) { return; } - TGmuCand.get(HIST("TGmuCand/hTrackPt"))->Fill(daughter[0].Pt()); - TGmuCand.get(HIST("TGmuCand/hTrackPt"))->Fill(daughter[1].Pt()); - TGmuCand.get(HIST("TGmuCand/hTrackEta"))->Fill(daughter[0].Eta()); - TGmuCand.get(HIST("TGmuCand/hTrackEta"))->Fill(daughter[1].Eta()); - TGmuCand.get(HIST("TGmuCand/hTrackPhi"))->Fill(daughter[0].Phi()); - TGmuCand.get(HIST("TGmuCand/hTrackPhi"))->Fill(daughter[1].Phi()); - TGmuCand.get(HIST("TGmuCand/hJpsiPt"))->Fill(mother.Pt()); + TGmuCand.get(HIST("TGmuCand/hTrackPt"))->Fill(trkDaughter1.pt()); + TGmuCand.get(HIST("TGmuCand/hTrackPt"))->Fill(trkDaughter2.pt()); + TGmuCand.get(HIST("TGmuCand/hTrackEta"))->Fill(RecoDecay::eta(daughter1)); + TGmuCand.get(HIST("TGmuCand/hTrackEta"))->Fill(RecoDecay::eta(daughter2)); + TGmuCand.get(HIST("TGmuCand/hTrackPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); + TGmuCand.get(HIST("TGmuCand/hTrackPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); + TGmuCand.get(HIST("TGmuCand/hJpsiPt"))->Fill(RecoDecay::pt(mother)); - if (mother.Pt() < 0.2f) { + if (RecoDecay::pt(mother) < 0.2f) { // fill track histos - JPsiToMu.get(HIST("JPsiToMu/Coherent/hPt1"))->Fill(daughter[0].Pt()); - JPsiToMu.get(HIST("JPsiToMu/Coherent/hPt2"))->Fill(daughter[1].Pt()); - JPsiToMu.get(HIST("JPsiToMu/Coherent/hEta1"))->Fill(daughter[0].Eta()); - JPsiToMu.get(HIST("JPsiToMu/Coherent/hEta2"))->Fill(daughter[1].Eta()); - JPsiToMu.get(HIST("JPsiToMu/Coherent/hPhi1"))->Fill(daughter[0].Phi()); - JPsiToMu.get(HIST("JPsiToMu/Coherent/hPhi2"))->Fill(daughter[1].Phi()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hPt1"))->Fill(trkDaughter1.pt()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hPt2"))->Fill(trkDaughter2.pt()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hEta1"))->Fill(RecoDecay::eta(daughter1)); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hEta2"))->Fill(RecoDecay::eta(daughter2)); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hPhi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hPhi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); if (trkDaughter1.hasTPC()) { - JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsPt"))->Fill(daughter[0].Pt(), trkDaughter1.tpcSignal()); - JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsEta"))->Fill(daughter[0].Eta(), trkDaughter1.tpcSignal()); - JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsPhi"))->Fill(daughter[0].Phi(), trkDaughter1.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsPt"))->Fill(trkDaughter1.pt(), trkDaughter1.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter1), trkDaughter1.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), trkDaughter1.tpcSignal()); } if (trkDaughter2.hasTPC()) { - JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsPt"))->Fill(daughter[1].Pt(), trkDaughter2.tpcSignal()); - JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsEta"))->Fill(daughter[1].Eta(), trkDaughter2.tpcSignal()); - JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsPhi"))->Fill(daughter[1].Phi(), trkDaughter2.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsPt"))->Fill(trkDaughter2.pt(), trkDaughter2.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter2), trkDaughter2.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), trkDaughter2.tpcSignal()); } // fill J/psi histos - JPsiToMu.get(HIST("JPsiToMu/Coherent/hPt"))->Fill(mother.Pt()); - JPsiToMu.get(HIST("JPsiToMu/Coherent/hEta"))->Fill(mother.Eta()); - JPsiToMu.get(HIST("JPsiToMu/Coherent/hPhi"))->Fill(mother.Phi()); - JPsiToMu.get(HIST("JPsiToMu/Coherent/hRap"))->Fill(mother.Rapidity()); - JPsiToMu.get(HIST("JPsiToMu/Coherent/hIVM"))->Fill(mother.M()); - - float* q = correlation(&daughter[0], &daughter[1], &mother); - Correlation.get(HIST("Correlation/Muon/Coherent/Phi1"))->Fill(daughter[0].Phi(), 1.); - Correlation.get(HIST("Correlation/Muon/Coherent/Phi2"))->Fill(daughter[1].Phi(), 1.); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hPt"))->Fill(RecoDecay::pt(mother)); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hEta"))->Fill(RecoDecay::eta(std::array{mother[0], mother[1], mother[2]})); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hPhi"))->Fill(RecoDecay::phi(mother[0], mother[1])); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hRap"))->Fill(RecoDecay::y(std::array{mother[0], mother[1], mother[2]}, mother[3])); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hIVM"))->Fill(RecoDecay::m(std::array{mother[0], mother[1], mother[2]}, energy)); + + float* q = correlation(&daughter[0], &daughter[1], &mom); + Correlation.get(HIST("Correlation/Muon/Coherent/Phi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), 1.); + Correlation.get(HIST("Correlation/Muon/Coherent/Phi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), 1.); Correlation.get(HIST("Correlation/Muon/Coherent/Phi"))->Fill(q[1], 1.); Correlation.get(HIST("Correlation/Muon/Coherent/CosTheta"))->Fill(q[2], 1.); Correlation.get(HIST("Correlation/Muon/Coherent/AccoplAngle"))->Fill(q[0], 1.); @@ -721,34 +733,34 @@ struct upcJpsiCentralBarrel { double dp = DeltaPhi(daughter[0], daughter[1]); Asymmetry.get(HIST("Asymmetry/Muon/Coherent/DeltaPhi"))->Fill(dp); } - if (mother.Pt() > 0.2f) { + if (RecoDecay::pt(mother) > 0.2f) { // fill track histos - JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPt1"))->Fill(daughter[0].Pt()); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPt2"))->Fill(daughter[1].Pt()); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/hEta1"))->Fill(daughter[0].Eta()); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/hEta2"))->Fill(daughter[1].Eta()); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPhi1"))->Fill(daughter[0].Phi()); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPhi2"))->Fill(daughter[1].Phi()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPt1"))->Fill(trkDaughter1.pt()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPt2"))->Fill(trkDaughter2.pt()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hEta1"))->Fill(RecoDecay::eta(daughter1)); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hEta2"))->Fill(RecoDecay::eta(daughter2)); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPhi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPhi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); if (trkDaughter1.hasTPC()) { - JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsPt"))->Fill(daughter[0].Pt(), trkDaughter1.tpcSignal()); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsEta"))->Fill(daughter[0].Eta(), trkDaughter1.tpcSignal()); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsPhi"))->Fill(daughter[0].Phi(), trkDaughter1.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsPt"))->Fill(trkDaughter1.pt(), trkDaughter1.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter1), trkDaughter1.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), trkDaughter1.tpcSignal()); } if (trkDaughter2.hasTPC()) { - JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsPt"))->Fill(daughter[1].Pt(), trkDaughter2.tpcSignal()); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsEta"))->Fill(daughter[1].Eta(), trkDaughter2.tpcSignal()); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsPhi"))->Fill(daughter[1].Phi(), trkDaughter2.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsPt"))->Fill(trkDaughter2.pt(), trkDaughter2.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter2), trkDaughter2.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), trkDaughter2.tpcSignal()); } // fill J/psi histos - JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPt"))->Fill(mother.Pt()); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/hEta"))->Fill(mother.Eta()); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPhi"))->Fill(mother.Phi()); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/hRap"))->Fill(mother.Rapidity()); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/hIVM"))->Fill(mother.M()); - - float* q = correlation(&daughter[0], &daughter[1], &mother); - Correlation.get(HIST("Correlation/Muon/Incoherent/Phi1"))->Fill(daughter[0].Phi(), 1.); - Correlation.get(HIST("Correlation/Muon/Incoherent/Phi2"))->Fill(daughter[1].Phi(), 1.); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPt"))->Fill(RecoDecay::pt(mother)); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hEta"))->Fill(RecoDecay::eta(std::array{mother[0], mother[1], mother[2]})); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPhi"))->Fill(RecoDecay::phi(mother[0], mother[1])); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hRap"))->Fill(RecoDecay::y(std::array{mother[0], mother[1], mother[2]}, mother[3])); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hIVM"))->Fill(RecoDecay::m(std::array{mother[0], mother[1], mother[2]}, energy)); + + float* q = correlation(&daughter[0], &daughter[1], &mom); + Correlation.get(HIST("Correlation/Muon/Incoherent/Phi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), 1.); + Correlation.get(HIST("Correlation/Muon/Incoherent/Phi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), 1.); Correlation.get(HIST("Correlation/Muon/Incoherent/Phi"))->Fill(q[1], 1.); Correlation.get(HIST("Correlation/Muon/Incoherent/CosTheta"))->Fill(q[2], 1.); Correlation.get(HIST("Correlation/Muon/Incoherent/AccoplAngle"))->Fill(q[0], 1.); @@ -759,84 +771,91 @@ struct upcJpsiCentralBarrel { } } // end muons if (countGTp == 2) { - TLorentzVector mother, daughter[2]; + TLorentzVector mom, daughter[2]; auto trkDaughter1 = tracks.iteratorAt(trkIdx[0]); auto trkDaughter2 = tracks.iteratorAt(trkIdx[1]); - auto ene1 = sqrt(pow(trkDaughter1.px(), 2.) + pow(trkDaughter1.py(), 2.) + pow(trkDaughter1.pz(), 2.) + pow(particleMass(pdg, 2212), 2.)); - auto ene2 = sqrt(pow(trkDaughter2.px(), 2.) + pow(trkDaughter2.py(), 2.) + pow(trkDaughter2.pz(), 2.) + pow(particleMass(pdg, 2212), 2.)); + auto ene1 = RecoDecay::e(trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz(), particleMass(pdg, 11)); + auto ene2 = RecoDecay::e(trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz(), particleMass(pdg, 11)); daughter[0].SetPxPyPzE(trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz(), ene1); daughter[1].SetPxPyPzE(trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz(), ene2); - mother = daughter[0] + daughter[1]; + mom = daughter[0] + daughter[1]; + + std::array daughter1 = {trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz()}; + std::array daughter2 = {trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz()}; + + std::array mother = {trkDaughter1.px() + trkDaughter2.px(), trkDaughter1.py() + trkDaughter2.py(), trkDaughter1.pz() + trkDaughter2.pz(), 2 * particleMass(pdg, 11)}; + + double energy = RecoDecay::e(mother[0], mother[1], mother[2], mother[3]); - TGp.get(HIST("TGp/hTrackPt"))->Fill(daughter[0].Pt()); - TGp.get(HIST("TGp/hTrackPt"))->Fill(daughter[1].Pt()); - TGp.get(HIST("TGp/hTrackEta"))->Fill(daughter[0].Eta()); - TGp.get(HIST("TGp/hTrackEta"))->Fill(daughter[1].Eta()); - TGp.get(HIST("TGp/hTrackPhi"))->Fill(daughter[0].Phi()); - TGp.get(HIST("TGp/hTrackPhi"))->Fill(daughter[1].Phi()); + TGp.get(HIST("TGp/hTrackPt"))->Fill(trkDaughter1.pt()); + TGp.get(HIST("TGp/hTrackPt"))->Fill(trkDaughter2.pt()); + TGp.get(HIST("TGp/hTrackEta"))->Fill(RecoDecay::eta(daughter1)); + TGp.get(HIST("TGp/hTrackEta"))->Fill(RecoDecay::eta(daughter2)); + TGp.get(HIST("TGp/hTrackPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); + TGp.get(HIST("TGp/hTrackPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); if (CandidateCuts(mother) != 1) { return; } - TGpCand.get(HIST("TGpCand/hTrackPt"))->Fill(daughter[0].Pt()); - TGpCand.get(HIST("TGpCand/hTrackPt"))->Fill(daughter[1].Pt()); - TGpCand.get(HIST("TGpCand/hTrackEta"))->Fill(daughter[0].Eta()); - TGpCand.get(HIST("TGpCand/hTrackEta"))->Fill(daughter[1].Eta()); - TGpCand.get(HIST("TGpCand/hTrackPhi"))->Fill(daughter[0].Phi()); - TGpCand.get(HIST("TGpCand/hTrackPhi"))->Fill(daughter[1].Phi()); - TGpCand.get(HIST("TGpCand/hJpsiPt"))->Fill(mother.Pt()); + TGpCand.get(HIST("TGpCand/hTrackPt"))->Fill(trkDaughter1.pt()); + TGpCand.get(HIST("TGpCand/hTrackPt"))->Fill(trkDaughter2.pt()); + TGpCand.get(HIST("TGpCand/hTrackEta"))->Fill(RecoDecay::eta(daughter1)); + TGpCand.get(HIST("TGpCand/hTrackEta"))->Fill(RecoDecay::eta(daughter2)); + TGpCand.get(HIST("TGpCand/hTrackPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); + TGpCand.get(HIST("TGpCand/hTrackPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); + TGpCand.get(HIST("TGpCand/hJpsiPt"))->Fill(RecoDecay::pt(mother)); - if (mother.Pt() < 0.2f) { + if (RecoDecay::pt(mother) < 0.2f) { // fill track histos - JPsiToP.get(HIST("JPsiToP/Coherent/hPt1"))->Fill(daughter[0].Pt()); - JPsiToP.get(HIST("JPsiToP/Coherent/hPt2"))->Fill(daughter[1].Pt()); - JPsiToP.get(HIST("JPsiToP/Coherent/hEta1"))->Fill(daughter[0].Eta()); - JPsiToP.get(HIST("JPsiToP/Coherent/hEta2"))->Fill(daughter[1].Eta()); - JPsiToP.get(HIST("JPsiToP/Coherent/hPhi1"))->Fill(daughter[0].Phi()); - JPsiToP.get(HIST("JPsiToP/Coherent/hPhi2"))->Fill(daughter[1].Phi()); + JPsiToP.get(HIST("JPsiToP/Coherent/hPt1"))->Fill(trkDaughter1.pt()); + JPsiToP.get(HIST("JPsiToP/Coherent/hPt2"))->Fill(trkDaughter2.pt()); + JPsiToP.get(HIST("JPsiToP/Coherent/hEta1"))->Fill(RecoDecay::eta(daughter1)); + JPsiToP.get(HIST("JPsiToP/Coherent/hEta2"))->Fill(RecoDecay::eta(daughter2)); + JPsiToP.get(HIST("JPsiToP/Coherent/hPhi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); + JPsiToP.get(HIST("JPsiToP/Coherent/hPhi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); if (trkDaughter1.hasTPC()) { - JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsPt"))->Fill(daughter[0].Pt(), trkDaughter1.tpcSignal()); - JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsEta"))->Fill(daughter[0].Eta(), trkDaughter1.tpcSignal()); - JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsPhi"))->Fill(daughter[0].Phi(), trkDaughter1.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsPt"))->Fill(trkDaughter1.pt(), trkDaughter1.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter1), trkDaughter1.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), trkDaughter1.tpcSignal()); } if (trkDaughter2.hasTPC()) { - JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsPt"))->Fill(daughter[1].Pt(), trkDaughter2.tpcSignal()); - JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsEta"))->Fill(daughter[1].Eta(), trkDaughter2.tpcSignal()); - JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsPhi"))->Fill(daughter[1].Phi(), trkDaughter2.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsPt"))->Fill(trkDaughter2.pt(), trkDaughter2.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter2), trkDaughter2.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), trkDaughter2.tpcSignal()); } // fill J/psi histos - JPsiToP.get(HIST("JPsiToP/Coherent/hPt"))->Fill(mother.Pt()); - JPsiToP.get(HIST("JPsiToP/Coherent/hEta"))->Fill(mother.Eta()); - JPsiToP.get(HIST("JPsiToP/Coherent/hPhi"))->Fill(mother.Phi()); - JPsiToP.get(HIST("JPsiToP/Coherent/hRap"))->Fill(mother.Rapidity()); - JPsiToP.get(HIST("JPsiToP/Coherent/hIVM"))->Fill(mother.M()); + JPsiToP.get(HIST("JPsiToP/Coherent/hPt"))->Fill(RecoDecay::pt(mother)); + JPsiToP.get(HIST("JPsiToP/Coherent/hEta"))->Fill(RecoDecay::eta(std::array{mother[0], mother[1], mother[2]})); + JPsiToP.get(HIST("JPsiToP/Coherent/hPhi"))->Fill(RecoDecay::phi(mother[0], mother[1])); + JPsiToP.get(HIST("JPsiToP/Coherent/hRap"))->Fill(RecoDecay::y(std::array{mother[0], mother[1], mother[2]}, mother[3])); + JPsiToP.get(HIST("JPsiToP/Coherent/hIVM"))->Fill(RecoDecay::m(std::array{mother[0], mother[1], mother[2]}, energy)); } - if (mother.Pt() > 0.2f) { + if (RecoDecay::pt(mother) > 0.2f) { // fill track histos - JPsiToP.get(HIST("JPsiToP/Incoherent/hPt1"))->Fill(daughter[0].Pt()); - JPsiToP.get(HIST("JPsiToP/Incoherent/hPt2"))->Fill(daughter[1].Pt()); - JPsiToP.get(HIST("JPsiToP/Incoherent/hEta1"))->Fill(daughter[0].Eta()); - JPsiToP.get(HIST("JPsiToP/Incoherent/hEta2"))->Fill(daughter[1].Eta()); - JPsiToP.get(HIST("JPsiToP/Incoherent/hPhi1"))->Fill(daughter[0].Phi()); - JPsiToP.get(HIST("JPsiToP/Incoherent/hPhi2"))->Fill(daughter[1].Phi()); + JPsiToP.get(HIST("JPsiToP/Incoherent/hPt1"))->Fill(trkDaughter1.pt()); + JPsiToP.get(HIST("JPsiToP/Incoherent/hPt2"))->Fill(trkDaughter2.pt()); + JPsiToP.get(HIST("JPsiToP/Incoherent/hEta1"))->Fill(RecoDecay::eta(daughter1)); + JPsiToP.get(HIST("JPsiToP/Incoherent/hEta2"))->Fill(RecoDecay::eta(daughter2)); + JPsiToP.get(HIST("JPsiToP/Incoherent/hPhi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); + JPsiToP.get(HIST("JPsiToP/Incoherent/hPhi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); if (trkDaughter1.hasTPC()) { - JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsPt"))->Fill(daughter[0].Pt(), trkDaughter1.tpcSignal()); - JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsEta"))->Fill(daughter[0].Eta(), trkDaughter1.tpcSignal()); - JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsPhi"))->Fill(daughter[0].Phi(), trkDaughter1.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsPt"))->Fill(trkDaughter1.pt(), trkDaughter1.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter1), trkDaughter1.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), trkDaughter1.tpcSignal()); } if (trkDaughter2.hasTPC()) { - JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsPt"))->Fill(daughter[1].Pt(), trkDaughter2.tpcSignal()); - JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsEta"))->Fill(daughter[1].Eta(), trkDaughter2.tpcSignal()); - JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsPhi"))->Fill(daughter[1].Phi(), trkDaughter2.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsPt"))->Fill(trkDaughter2.pt(), trkDaughter2.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter2), trkDaughter2.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), trkDaughter2.tpcSignal()); } // fill J/psi histos - JPsiToP.get(HIST("JPsiToP/Incoherent/hPt"))->Fill(mother.Pt()); - JPsiToP.get(HIST("JPsiToP/Incoherent/hEta"))->Fill(mother.Eta()); - JPsiToP.get(HIST("JPsiToP/Incoherent/hPhi"))->Fill(mother.Phi()); - JPsiToP.get(HIST("JPsiToP/Incoherent/hRap"))->Fill(mother.Rapidity()); - JPsiToP.get(HIST("JPsiToP/Incoherent/hIVM"))->Fill(mother.M()); + JPsiToP.get(HIST("JPsiToP/Incoherent/hPt"))->Fill(RecoDecay::pt(mother)); + JPsiToP.get(HIST("JPsiToP/Incoherent/hEta"))->Fill(RecoDecay::eta(std::array{mother[0], mother[1], mother[2]})); + JPsiToP.get(HIST("JPsiToP/Incoherent/hPhi"))->Fill(RecoDecay::phi(mother[0], mother[1])); + JPsiToP.get(HIST("JPsiToP/Incoherent/hRap"))->Fill(RecoDecay::y(std::array{mother[0], mother[1], mother[2]}, mother[3])); + JPsiToP.get(HIST("JPsiToP/Incoherent/hIVM"))->Fill(RecoDecay::m(std::array{mother[0], mother[1], mother[2]}, energy)); } } // end protons } // end two tracks From c65a1bb6a0f564d18e40353d329ca969ff8d1b42 Mon Sep 17 00:00:00 2001 From: Sara Haidlova Date: Tue, 9 Apr 2024 21:05:47 +0200 Subject: [PATCH 5/6] PWGUD: fix for memory leak and mass calculations --- PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h | 10 +- PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx | 204 +++++++++++--------- 2 files changed, 116 insertions(+), 98 deletions(-) diff --git a/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h b/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h index 15ea539ebf9..54b6ec1b276 100644 --- a/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h +++ b/PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h @@ -56,7 +56,7 @@ float* correlation(TLorentzVector* lv1, TLorentzVector* lv2, TLorentzVector* lv) // Accoplanarity angle Float_t deltaPhi; deltaPhi = lv1->Phi() - lv2->Phi(); - float accoplCut = 1. - abs(deltaPhi) / o2::constants::math::PI; + float accoplCut = 1. - std::abs(deltaPhi) / o2::constants::math::PI; // z TLorentzVector z; Float_t part1, part2; @@ -83,7 +83,7 @@ float* correlation(TLorentzVector* lv1, TLorentzVector* lv2, TLorentzVector* lv) z = part3 - part4; // Normalized z - Float_t normz = TMath::Sqrt(-z * z); + Float_t normz = std::sqrt(-z * z); Float_t znx = z.X() / normz; Float_t zny = z.Y() / normz; Float_t znz = z.Z() / normz; @@ -115,7 +115,7 @@ float* correlation(TLorentzVector* lv1, TLorentzVector* lv2, TLorentzVector* lv) // Un-normalized x x = *lv - comp1 - comp2; // normalize x - Float_t normx = TMath::Sqrt(-x * x); + Float_t normx = std::sqrt(-x * x); Float_t xnx = x.X() / normx; Float_t xny = x.Y() / normx; Float_t xnz = x.Z() / normx; @@ -135,7 +135,7 @@ float* correlation(TLorentzVector* lv1, TLorentzVector* lv2, TLorentzVector* lv) TLorentzVector y(TVector3(yone, ytwo, ythree), yfour); // normalize y - Float_t normy = TMath::Sqrt(-y * y); + Float_t normy = std::sqrt(-y * y); Float_t ynx = y.X() / normy; Float_t yny = y.Y() / normy; Float_t ynz = y.Z() / normy; @@ -154,7 +154,7 @@ float* correlation(TLorentzVector* lv1, TLorentzVector* lv2, TLorentzVector* lv) TLorentzVector diff2(TVector3(diff2x, diff2y, diff2z), diff2e); // Normalize diff2 - Float_t norm2 = TMath::Sqrt(-diff2 * diff2); + Float_t norm2 = std::sqrt(-diff2 * diff2); Float_t diff3x = diff2.X() / norm2; Float_t diff3y = diff2.Y() / norm2; Float_t diff3z = diff2.Z() / norm2; diff --git a/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx b/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx index 74fa86f360b..d911ad90dc5 100644 --- a/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx +++ b/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx @@ -379,7 +379,7 @@ struct upcJpsiCentralBarrel { bool GoodTrackCuts(T const& track) { // kinematics - if (RecoDecay::eta(std::array{track.px(), track.py(), track.pz()}) < -EtaCut || RecoDecay::eta(std::array{track.px(), track.py(), track.pz()}) > EtaCut) { + if (std::abs(RecoDecay::eta(std::array{track.px(), track.py(), track.pz()})) > EtaCut) { return false; } // DCA @@ -413,15 +413,14 @@ struct upcJpsiCentralBarrel { return true; } - template - bool CandidateCuts(C const& mother) + // template + bool CandidateCuts(float massJpsi, float rapJpsi) { - if (abs(RecoDecay::y(std::array{mother[0], mother[1], mother[2]}, mother[3])) > 0.8f) { + if (rapJpsi > RapCut) { return false; } - double energy = RecoDecay::e(mother[0], mother[1], mother[2], mother[3]); - if (RecoDecay::m(std::array{mother[0], mother[1], mother[2]}, energy) < 2.5f) { + if (massJpsi < 2.5f) { return false; } @@ -522,6 +521,10 @@ struct upcJpsiCentralBarrel { Statistics.get(HIST("Statistics/hNumberGTmu"))->Fill(countGTmu); Statistics.get(HIST("Statistics/hNumberGTp"))->Fill(countGTp); + float massEl = 0.000510998950; + float massMu = 0.10565837; + float massPr = 0.93827208816; + if (countGT == 2) { TLorentzVector mom, daughter[2]; auto trkDaughter1 = tracks.iteratorAt(trkIdx[0]); @@ -533,8 +536,8 @@ struct upcJpsiCentralBarrel { if (!(RecoDecay::sumOfSquares(trkDaughter1.tpcNSigmaMu(), trkDaughter2.tpcNSigmaMu()) < RecoDecay::sumOfSquares(trkDaughter1.tpcNSigmaEl(), trkDaughter2.tpcNSigmaEl()))) { return; } - auto ene1 = RecoDecay::e(trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz(), particleMass(pdg, 11)); - auto ene2 = RecoDecay::e(trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz(), particleMass(pdg, 11)); + auto ene1 = RecoDecay::e(trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz(), massEl); + auto ene2 = RecoDecay::e(trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz(), massEl); daughter[0].SetPxPyPzE(trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz(), ene1); daughter[1].SetPxPyPzE(trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz(), ene2); mom = daughter[0] + daughter[1]; @@ -542,17 +545,20 @@ struct upcJpsiCentralBarrel { std::array daughter1 = {trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz()}; std::array daughter2 = {trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz()}; - std::array mother = {trkDaughter1.px() + trkDaughter2.px(), trkDaughter1.py() + trkDaughter2.py(), trkDaughter1.pz() + trkDaughter2.pz(), 2 * particleMass(pdg, 11)}; - double energy = RecoDecay::e(mother[0], mother[1], mother[2], mother[3]); + std::array mother = {trkDaughter1.px() + trkDaughter2.px(), trkDaughter1.py() + trkDaughter2.py(), trkDaughter1.pz() + trkDaughter2.pz()}; + + auto arrMom = std::array{daughter1, daughter2}; + float massJpsi = RecoDecay::m(arrMom, std::array{massEl, massEl}); + float rapJpsi = RecoDecay::y(mother, massJpsi); TGel.get(HIST("TGel/hTrackPt"))->Fill(trkDaughter1.pt()); TGel.get(HIST("TGel/hTrackPt"))->Fill(trkDaughter2.pt()); TGel.get(HIST("TGel/hTrackEta"))->Fill(RecoDecay::eta(daughter1)); TGel.get(HIST("TGel/hTrackEta"))->Fill(RecoDecay::eta(daughter2)); - TGel.get(HIST("TGel/hTrackPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); - TGel.get(HIST("TGel/hTrackPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); + TGel.get(HIST("TGel/hTrackPhi"))->Fill(RecoDecay::phi(daughter1)); + TGel.get(HIST("TGel/hTrackPhi"))->Fill(RecoDecay::phi(daughter2)); - if (CandidateCuts(mother) != 1) { + if (CandidateCuts(massJpsi, rapJpsi) != 1) { return; } @@ -560,8 +566,8 @@ struct upcJpsiCentralBarrel { TGelCand.get(HIST("TGelCand/hTrackPt"))->Fill(trkDaughter2.pt()); TGelCand.get(HIST("TGelCand/hTrackEta"))->Fill(RecoDecay::eta(daughter1)); TGelCand.get(HIST("TGelCand/hTrackEta"))->Fill(RecoDecay::eta(daughter2)); - TGelCand.get(HIST("TGelCand/hTrackPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); - TGelCand.get(HIST("TGelCand/hTrackPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); + TGelCand.get(HIST("TGelCand/hTrackPhi"))->Fill(RecoDecay::phi(daughter1)); + TGelCand.get(HIST("TGelCand/hTrackPhi"))->Fill(RecoDecay::phi(daughter2)); TGelCand.get(HIST("TGelCand/hJpsiPt"))->Fill(RecoDecay::pt(mother)); if (RecoDecay::pt(mother) < 0.2f) { // fill track histos @@ -569,39 +575,39 @@ struct upcJpsiCentralBarrel { JPsiToEl.get(HIST("JPsiToEl/Coherent/hPt2"))->Fill(trkDaughter2.pt()); JPsiToEl.get(HIST("JPsiToEl/Coherent/hEta1"))->Fill(RecoDecay::eta(daughter1)); JPsiToEl.get(HIST("JPsiToEl/Coherent/hEta2"))->Fill(RecoDecay::eta(daughter2)); - JPsiToEl.get(HIST("JPsiToEl/Coherent/hPhi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); - JPsiToEl.get(HIST("JPsiToEl/Coherent/hPhi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hPhi1"))->Fill(RecoDecay::phi(daughter1)); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hPhi2"))->Fill(RecoDecay::phi(daughter2)); if (trkDaughter1.hasTPC()) { JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsPt"))->Fill(trkDaughter1.pt(), trkDaughter1.tpcSignal()); JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter1), trkDaughter1.tpcSignal()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), trkDaughter1.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter1), trkDaughter1.tpcSignal()); } if (trkDaughter2.hasTPC()) { JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsPt"))->Fill(trkDaughter2.pt(), trkDaughter2.tpcSignal()); JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter2), trkDaughter2.tpcSignal()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), trkDaughter2.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter2), trkDaughter2.tpcSignal()); } if (trkDaughter1.hasTOF()) { JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsPt"))->Fill(trkDaughter1.pt(), trkDaughter1.tofSignal()); JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsEta"))->Fill(RecoDecay::eta(daughter1), trkDaughter1.tofSignal()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), trkDaughter1.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsPhi"))->Fill(RecoDecay::phi(daughter1), trkDaughter1.tofSignal()); } if (trkDaughter2.hasTOF()) { JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsPt"))->Fill(trkDaughter2.pt(), trkDaughter2.tofSignal()); JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsEta"))->Fill(RecoDecay::eta(daughter2), trkDaughter2.tofSignal()); - JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), trkDaughter2.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Coherent/PID/hTOFVsPhi"))->Fill(RecoDecay::phi(daughter2), trkDaughter2.tofSignal()); } // fill J/psi histos JPsiToEl.get(HIST("JPsiToEl/Coherent/hPt"))->Fill(RecoDecay::pt(mother)); - JPsiToEl.get(HIST("JPsiToEl/Coherent/hEta"))->Fill(RecoDecay::eta(std::array{mother[0], mother[1], mother[2]})); - JPsiToEl.get(HIST("JPsiToEl/Coherent/hPhi"))->Fill(RecoDecay::phi(mother[0], mother[1])); - JPsiToEl.get(HIST("JPsiToEl/Coherent/hRap"))->Fill(RecoDecay::y(std::array{mother[0], mother[1], mother[2]}, mother[3])); - JPsiToEl.get(HIST("JPsiToEl/Coherent/hIVM"))->Fill(RecoDecay::m(std::array{mother[0], mother[1], mother[2]}, energy)); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hEta"))->Fill(RecoDecay::eta(mother)); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hPhi"))->Fill(RecoDecay::phi(mother)); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hRap"))->Fill(rapJpsi); + JPsiToEl.get(HIST("JPsiToEl/Coherent/hIVM"))->Fill(massJpsi); float* q = correlation(&daughter[0], &daughter[1], &mom); - Correlation.get(HIST("Correlation/Electron/Coherent/Phi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), 1.); - Correlation.get(HIST("Correlation/Electron/Coherent/Phi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), 1.); + Correlation.get(HIST("Correlation/Electron/Coherent/Phi1"))->Fill(RecoDecay::phi(daughter1), 1.); + Correlation.get(HIST("Correlation/Electron/Coherent/Phi2"))->Fill(RecoDecay::phi(daughter2), 1.); Correlation.get(HIST("Correlation/Electron/Coherent/Phi"))->Fill(q[1], 1.); Correlation.get(HIST("Correlation/Electron/Coherent/CosTheta"))->Fill(q[2], 1.); Correlation.get(HIST("Correlation/Electron/Coherent/AccoplAngle"))->Fill(q[0], 1.); @@ -609,6 +615,8 @@ struct upcJpsiCentralBarrel { double dp = DeltaPhi(daughter[0], daughter[1]); Asymmetry.get(HIST("Asymmetry/Electron/Coherent/DeltaPhi"))->Fill(dp); + + delete[] q; } if (RecoDecay::pt(mother) > 0.2f) { // fill track histos @@ -616,39 +624,39 @@ struct upcJpsiCentralBarrel { JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPt2"))->Fill(trkDaughter2.pt()); JPsiToEl.get(HIST("JPsiToEl/Incoherent/hEta1"))->Fill(RecoDecay::eta(daughter1)); JPsiToEl.get(HIST("JPsiToEl/Incoherent/hEta2"))->Fill(RecoDecay::eta(daughter2)); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPhi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPhi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPhi1"))->Fill(RecoDecay::phi(daughter1)); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPhi2"))->Fill(RecoDecay::phi(daughter2)); if (trkDaughter1.hasTPC()) { JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsPt"))->Fill(trkDaughter1.pt(), trkDaughter1.tpcSignal()); JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter1), trkDaughter1.tpcSignal()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), trkDaughter1.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter1), trkDaughter1.tpcSignal()); } if (trkDaughter2.hasTPC()) { JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsPt"))->Fill(trkDaughter2.pt(), trkDaughter2.tpcSignal()); JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter2), trkDaughter2.tpcSignal()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), trkDaughter2.tpcSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter2), trkDaughter2.tpcSignal()); } if (trkDaughter1.hasTOF()) { JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsPt"))->Fill(trkDaughter1.pt(), trkDaughter1.tofSignal()); JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsEta"))->Fill(RecoDecay::eta(daughter1), trkDaughter1.tofSignal()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), trkDaughter1.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsPhi"))->Fill(RecoDecay::phi(daughter1), trkDaughter1.tofSignal()); } if (trkDaughter2.hasTOF()) { JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsPt"))->Fill(trkDaughter2.pt(), trkDaughter2.tofSignal()); JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsEta"))->Fill(RecoDecay::eta(daughter2), trkDaughter2.tofSignal()); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), trkDaughter2.tofSignal()); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/PID/hTOFVsPhi"))->Fill(RecoDecay::phi(daughter2), trkDaughter2.tofSignal()); } // fill J/psi histos JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPt"))->Fill(RecoDecay::pt(mother)); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/hEta"))->Fill(RecoDecay::eta(std::array{mother[0], mother[1], mother[2]})); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPhi"))->Fill(RecoDecay::phi(mother[0], mother[1])); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/hRap"))->Fill(RecoDecay::y(std::array{mother[0], mother[1], mother[2]}, mother[3])); - JPsiToEl.get(HIST("JPsiToEl/Incoherent/hIVM"))->Fill(RecoDecay::m(std::array{mother[0], mother[1], mother[2]}, energy)); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hEta"))->Fill(RecoDecay::eta(mother)); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hPhi"))->Fill(RecoDecay::phi(mother)); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hRap"))->Fill(rapJpsi); + JPsiToEl.get(HIST("JPsiToEl/Incoherent/hIVM"))->Fill(massJpsi); float* q = correlation(&daughter[0], &daughter[1], &mom); - Correlation.get(HIST("Correlation/Electron/Incoherent/Phi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), 1.); - Correlation.get(HIST("Correlation/Electron/Incoherent/Phi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), 1.); + Correlation.get(HIST("Correlation/Electron/Incoherent/Phi1"))->Fill(RecoDecay::phi(daughter1), 1.); + Correlation.get(HIST("Correlation/Electron/Incoherent/Phi2"))->Fill(RecoDecay::phi(daughter2), 1.); Correlation.get(HIST("Correlation/Electron/Incoherent/Phi"))->Fill(q[1], 1.); Correlation.get(HIST("Correlation/Electron/Incoherent/CosTheta"))->Fill(q[2], 1.); Correlation.get(HIST("Correlation/Electron/Incoherent/AccoplAngle"))->Fill(q[0], 1.); @@ -656,6 +664,8 @@ struct upcJpsiCentralBarrel { double dp = DeltaPhi(daughter[0], daughter[1]); Asymmetry.get(HIST("Asymmetry/Electron/Incoherent/DeltaPhi"))->Fill(dp); + + delete[] q; } } // end electrons if (countGTmu == 2) { @@ -665,8 +675,8 @@ struct upcJpsiCentralBarrel { if (!(RecoDecay::sumOfSquares(trkDaughter1.tpcNSigmaEl(), trkDaughter2.tpcNSigmaEl() < RecoDecay::sumOfSquares(trkDaughter1.tpcNSigmaMu(), trkDaughter2.tpcNSigmaMu())))) { return; } - auto ene1 = RecoDecay::e(trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz(), particleMass(pdg, 11)); - auto ene2 = RecoDecay::e(trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz(), particleMass(pdg, 11)); + auto ene1 = RecoDecay::e(trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz(), massMu); + auto ene2 = RecoDecay::e(trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz(), massMu); daughter[0].SetPxPyPzE(trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz(), ene1); daughter[1].SetPxPyPzE(trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz(), ene2); mom = daughter[0] + daughter[1]; @@ -674,18 +684,20 @@ struct upcJpsiCentralBarrel { std::array daughter1 = {trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz()}; std::array daughter2 = {trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz()}; - std::array mother = {trkDaughter1.px() + trkDaughter2.px(), trkDaughter1.py() + trkDaughter2.py(), trkDaughter1.pz() + trkDaughter2.pz(), 2 * particleMass(pdg, 11)}; + std::array mother = {trkDaughter1.px() + trkDaughter2.px(), trkDaughter1.py() + trkDaughter2.py(), trkDaughter1.pz() + trkDaughter2.pz()}; - double energy = RecoDecay::e(mother[0], mother[1], mother[2], mother[3]); + auto arrMom = std::array{daughter1, daughter2}; + float massJpsi = RecoDecay::m(arrMom, std::array{massEl, massEl}); + float rapJpsi = RecoDecay::y(mother, massJpsi); TGmu.get(HIST("TGmu/hTrackPt"))->Fill(trkDaughter1.pt()); TGmu.get(HIST("TGmu/hTrackPt"))->Fill(trkDaughter2.pt()); TGmu.get(HIST("TGmu/hTrackEta"))->Fill(RecoDecay::eta(daughter1)); TGmu.get(HIST("TGmu/hTrackEta"))->Fill(RecoDecay::eta(daughter2)); - TGmu.get(HIST("TGmu/hTrackPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); - TGmu.get(HIST("TGmu/hTrackPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); + TGmu.get(HIST("TGmu/hTrackPhi"))->Fill(RecoDecay::phi(daughter1)); + TGmu.get(HIST("TGmu/hTrackPhi"))->Fill(RecoDecay::phi(daughter2)); - if (CandidateCuts(mother) != 1) { + if (CandidateCuts(massJpsi, rapJpsi) != 1) { return; } @@ -693,8 +705,8 @@ struct upcJpsiCentralBarrel { TGmuCand.get(HIST("TGmuCand/hTrackPt"))->Fill(trkDaughter2.pt()); TGmuCand.get(HIST("TGmuCand/hTrackEta"))->Fill(RecoDecay::eta(daughter1)); TGmuCand.get(HIST("TGmuCand/hTrackEta"))->Fill(RecoDecay::eta(daughter2)); - TGmuCand.get(HIST("TGmuCand/hTrackPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); - TGmuCand.get(HIST("TGmuCand/hTrackPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); + TGmuCand.get(HIST("TGmuCand/hTrackPhi"))->Fill(RecoDecay::phi(daughter1)); + TGmuCand.get(HIST("TGmuCand/hTrackPhi"))->Fill(RecoDecay::phi(daughter2)); TGmuCand.get(HIST("TGmuCand/hJpsiPt"))->Fill(RecoDecay::pt(mother)); if (RecoDecay::pt(mother) < 0.2f) { @@ -703,28 +715,28 @@ struct upcJpsiCentralBarrel { JPsiToMu.get(HIST("JPsiToMu/Coherent/hPt2"))->Fill(trkDaughter2.pt()); JPsiToMu.get(HIST("JPsiToMu/Coherent/hEta1"))->Fill(RecoDecay::eta(daughter1)); JPsiToMu.get(HIST("JPsiToMu/Coherent/hEta2"))->Fill(RecoDecay::eta(daughter2)); - JPsiToMu.get(HIST("JPsiToMu/Coherent/hPhi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); - JPsiToMu.get(HIST("JPsiToMu/Coherent/hPhi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hPhi1"))->Fill(RecoDecay::phi(daughter1)); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hPhi2"))->Fill(RecoDecay::phi(daughter2)); if (trkDaughter1.hasTPC()) { JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsPt"))->Fill(trkDaughter1.pt(), trkDaughter1.tpcSignal()); JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter1), trkDaughter1.tpcSignal()); - JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), trkDaughter1.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter1), trkDaughter1.tpcSignal()); } if (trkDaughter2.hasTPC()) { JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsPt"))->Fill(trkDaughter2.pt(), trkDaughter2.tpcSignal()); JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter2), trkDaughter2.tpcSignal()); - JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), trkDaughter2.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Coherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter2), trkDaughter2.tpcSignal()); } // fill J/psi histos JPsiToMu.get(HIST("JPsiToMu/Coherent/hPt"))->Fill(RecoDecay::pt(mother)); - JPsiToMu.get(HIST("JPsiToMu/Coherent/hEta"))->Fill(RecoDecay::eta(std::array{mother[0], mother[1], mother[2]})); - JPsiToMu.get(HIST("JPsiToMu/Coherent/hPhi"))->Fill(RecoDecay::phi(mother[0], mother[1])); - JPsiToMu.get(HIST("JPsiToMu/Coherent/hRap"))->Fill(RecoDecay::y(std::array{mother[0], mother[1], mother[2]}, mother[3])); - JPsiToMu.get(HIST("JPsiToMu/Coherent/hIVM"))->Fill(RecoDecay::m(std::array{mother[0], mother[1], mother[2]}, energy)); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hEta"))->Fill(RecoDecay::eta(mother)); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hPhi"))->Fill(RecoDecay::phi(mother)); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hRap"))->Fill(rapJpsi); + JPsiToMu.get(HIST("JPsiToMu/Coherent/hIVM"))->Fill(massJpsi); float* q = correlation(&daughter[0], &daughter[1], &mom); - Correlation.get(HIST("Correlation/Muon/Coherent/Phi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), 1.); - Correlation.get(HIST("Correlation/Muon/Coherent/Phi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), 1.); + Correlation.get(HIST("Correlation/Muon/Coherent/Phi1"))->Fill(RecoDecay::phi(daughter1), 1.); + Correlation.get(HIST("Correlation/Muon/Coherent/Phi2"))->Fill(RecoDecay::phi(daughter2), 1.); Correlation.get(HIST("Correlation/Muon/Coherent/Phi"))->Fill(q[1], 1.); Correlation.get(HIST("Correlation/Muon/Coherent/CosTheta"))->Fill(q[2], 1.); Correlation.get(HIST("Correlation/Muon/Coherent/AccoplAngle"))->Fill(q[0], 1.); @@ -732,6 +744,8 @@ struct upcJpsiCentralBarrel { double dp = DeltaPhi(daughter[0], daughter[1]); Asymmetry.get(HIST("Asymmetry/Muon/Coherent/DeltaPhi"))->Fill(dp); + + delete[] q; } if (RecoDecay::pt(mother) > 0.2f) { // fill track histos @@ -739,28 +753,28 @@ struct upcJpsiCentralBarrel { JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPt2"))->Fill(trkDaughter2.pt()); JPsiToMu.get(HIST("JPsiToMu/Incoherent/hEta1"))->Fill(RecoDecay::eta(daughter1)); JPsiToMu.get(HIST("JPsiToMu/Incoherent/hEta2"))->Fill(RecoDecay::eta(daughter2)); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPhi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPhi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPhi1"))->Fill(RecoDecay::phi(daughter1)); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPhi2"))->Fill(RecoDecay::phi(daughter2)); if (trkDaughter1.hasTPC()) { JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsPt"))->Fill(trkDaughter1.pt(), trkDaughter1.tpcSignal()); JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter1), trkDaughter1.tpcSignal()); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), trkDaughter1.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter1), trkDaughter1.tpcSignal()); } if (trkDaughter2.hasTPC()) { JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsPt"))->Fill(trkDaughter2.pt(), trkDaughter2.tpcSignal()); JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter2), trkDaughter2.tpcSignal()); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), trkDaughter2.tpcSignal()); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter2), trkDaughter2.tpcSignal()); } // fill J/psi histos JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPt"))->Fill(RecoDecay::pt(mother)); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/hEta"))->Fill(RecoDecay::eta(std::array{mother[0], mother[1], mother[2]})); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPhi"))->Fill(RecoDecay::phi(mother[0], mother[1])); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/hRap"))->Fill(RecoDecay::y(std::array{mother[0], mother[1], mother[2]}, mother[3])); - JPsiToMu.get(HIST("JPsiToMu/Incoherent/hIVM"))->Fill(RecoDecay::m(std::array{mother[0], mother[1], mother[2]}, energy)); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hEta"))->Fill(RecoDecay::eta(mother)); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hPhi"))->Fill(RecoDecay::phi(mother)); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hRap"))->Fill(rapJpsi); + JPsiToMu.get(HIST("JPsiToMu/Incoherent/hIVM"))->Fill(massJpsi); float* q = correlation(&daughter[0], &daughter[1], &mom); - Correlation.get(HIST("Correlation/Muon/Incoherent/Phi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), 1.); - Correlation.get(HIST("Correlation/Muon/Incoherent/Phi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), 1.); + Correlation.get(HIST("Correlation/Muon/Incoherent/Phi1"))->Fill(RecoDecay::phi(daughter1), 1.); + Correlation.get(HIST("Correlation/Muon/Incoherent/Phi2"))->Fill(RecoDecay::phi(daughter2), 1.); Correlation.get(HIST("Correlation/Muon/Incoherent/Phi"))->Fill(q[1], 1.); Correlation.get(HIST("Correlation/Muon/Incoherent/CosTheta"))->Fill(q[2], 1.); Correlation.get(HIST("Correlation/Muon/Incoherent/AccoplAngle"))->Fill(q[0], 1.); @@ -768,6 +782,8 @@ struct upcJpsiCentralBarrel { double dp = DeltaPhi(daughter[0], daughter[1]); Asymmetry.get(HIST("Asymmetry/Muon/Incoherent/DeltaPhi"))->Fill(dp); + + delete[] q; } } // end muons if (countGTp == 2) { @@ -775,8 +791,8 @@ struct upcJpsiCentralBarrel { auto trkDaughter1 = tracks.iteratorAt(trkIdx[0]); auto trkDaughter2 = tracks.iteratorAt(trkIdx[1]); - auto ene1 = RecoDecay::e(trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz(), particleMass(pdg, 11)); - auto ene2 = RecoDecay::e(trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz(), particleMass(pdg, 11)); + auto ene1 = RecoDecay::e(trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz(), massPr); + auto ene2 = RecoDecay::e(trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz(), massPr); daughter[0].SetPxPyPzE(trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz(), ene1); daughter[1].SetPxPyPzE(trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz(), ene2); mom = daughter[0] + daughter[1]; @@ -784,18 +800,20 @@ struct upcJpsiCentralBarrel { std::array daughter1 = {trkDaughter1.px(), trkDaughter1.py(), trkDaughter1.pz()}; std::array daughter2 = {trkDaughter2.px(), trkDaughter2.py(), trkDaughter2.pz()}; - std::array mother = {trkDaughter1.px() + trkDaughter2.px(), trkDaughter1.py() + trkDaughter2.py(), trkDaughter1.pz() + trkDaughter2.pz(), 2 * particleMass(pdg, 11)}; + std::array mother = {trkDaughter1.px() + trkDaughter2.px(), trkDaughter1.py() + trkDaughter2.py(), trkDaughter1.pz() + trkDaughter2.pz()}; - double energy = RecoDecay::e(mother[0], mother[1], mother[2], mother[3]); + auto arrMom = std::array{daughter1, daughter2}; + float massJpsi = RecoDecay::m(arrMom, std::array{massEl, massEl}); + float rapJpsi = RecoDecay::y(mother, massJpsi); TGp.get(HIST("TGp/hTrackPt"))->Fill(trkDaughter1.pt()); TGp.get(HIST("TGp/hTrackPt"))->Fill(trkDaughter2.pt()); TGp.get(HIST("TGp/hTrackEta"))->Fill(RecoDecay::eta(daughter1)); TGp.get(HIST("TGp/hTrackEta"))->Fill(RecoDecay::eta(daughter2)); - TGp.get(HIST("TGp/hTrackPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); - TGp.get(HIST("TGp/hTrackPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); + TGp.get(HIST("TGp/hTrackPhi"))->Fill(RecoDecay::phi(daughter1)); + TGp.get(HIST("TGp/hTrackPhi"))->Fill(RecoDecay::phi(daughter2)); - if (CandidateCuts(mother) != 1) { + if (CandidateCuts(massJpsi, rapJpsi) != 1) { return; } @@ -803,8 +821,8 @@ struct upcJpsiCentralBarrel { TGpCand.get(HIST("TGpCand/hTrackPt"))->Fill(trkDaughter2.pt()); TGpCand.get(HIST("TGpCand/hTrackEta"))->Fill(RecoDecay::eta(daughter1)); TGpCand.get(HIST("TGpCand/hTrackEta"))->Fill(RecoDecay::eta(daughter2)); - TGpCand.get(HIST("TGpCand/hTrackPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); - TGpCand.get(HIST("TGpCand/hTrackPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); + TGpCand.get(HIST("TGpCand/hTrackPhi"))->Fill(RecoDecay::phi(daughter1)); + TGpCand.get(HIST("TGpCand/hTrackPhi"))->Fill(RecoDecay::phi(daughter2)); TGpCand.get(HIST("TGpCand/hJpsiPt"))->Fill(RecoDecay::pt(mother)); if (RecoDecay::pt(mother) < 0.2f) { @@ -813,24 +831,24 @@ struct upcJpsiCentralBarrel { JPsiToP.get(HIST("JPsiToP/Coherent/hPt2"))->Fill(trkDaughter2.pt()); JPsiToP.get(HIST("JPsiToP/Coherent/hEta1"))->Fill(RecoDecay::eta(daughter1)); JPsiToP.get(HIST("JPsiToP/Coherent/hEta2"))->Fill(RecoDecay::eta(daughter2)); - JPsiToP.get(HIST("JPsiToP/Coherent/hPhi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); - JPsiToP.get(HIST("JPsiToP/Coherent/hPhi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); + JPsiToP.get(HIST("JPsiToP/Coherent/hPhi1"))->Fill(RecoDecay::phi(daughter1)); + JPsiToP.get(HIST("JPsiToP/Coherent/hPhi2"))->Fill(RecoDecay::phi(daughter2)); if (trkDaughter1.hasTPC()) { JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsPt"))->Fill(trkDaughter1.pt(), trkDaughter1.tpcSignal()); JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter1), trkDaughter1.tpcSignal()); - JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), trkDaughter1.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter1), trkDaughter1.tpcSignal()); } if (trkDaughter2.hasTPC()) { JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsPt"))->Fill(trkDaughter2.pt(), trkDaughter2.tpcSignal()); JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter2), trkDaughter2.tpcSignal()); - JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), trkDaughter2.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Coherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter2), trkDaughter2.tpcSignal()); } // fill J/psi histos JPsiToP.get(HIST("JPsiToP/Coherent/hPt"))->Fill(RecoDecay::pt(mother)); - JPsiToP.get(HIST("JPsiToP/Coherent/hEta"))->Fill(RecoDecay::eta(std::array{mother[0], mother[1], mother[2]})); - JPsiToP.get(HIST("JPsiToP/Coherent/hPhi"))->Fill(RecoDecay::phi(mother[0], mother[1])); - JPsiToP.get(HIST("JPsiToP/Coherent/hRap"))->Fill(RecoDecay::y(std::array{mother[0], mother[1], mother[2]}, mother[3])); - JPsiToP.get(HIST("JPsiToP/Coherent/hIVM"))->Fill(RecoDecay::m(std::array{mother[0], mother[1], mother[2]}, energy)); + JPsiToP.get(HIST("JPsiToP/Coherent/hEta"))->Fill(RecoDecay::eta(mother)); + JPsiToP.get(HIST("JPsiToP/Coherent/hPhi"))->Fill(RecoDecay::phi(mother)); + JPsiToP.get(HIST("JPsiToP/Coherent/hRap"))->Fill(rapJpsi); + JPsiToP.get(HIST("JPsiToP/Coherent/hIVM"))->Fill(massJpsi); } if (RecoDecay::pt(mother) > 0.2f) { // fill track histos @@ -838,24 +856,24 @@ struct upcJpsiCentralBarrel { JPsiToP.get(HIST("JPsiToP/Incoherent/hPt2"))->Fill(trkDaughter2.pt()); JPsiToP.get(HIST("JPsiToP/Incoherent/hEta1"))->Fill(RecoDecay::eta(daughter1)); JPsiToP.get(HIST("JPsiToP/Incoherent/hEta2"))->Fill(RecoDecay::eta(daughter2)); - JPsiToP.get(HIST("JPsiToP/Incoherent/hPhi1"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1])); - JPsiToP.get(HIST("JPsiToP/Incoherent/hPhi2"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1])); + JPsiToP.get(HIST("JPsiToP/Incoherent/hPhi1"))->Fill(RecoDecay::phi(daughter1)); + JPsiToP.get(HIST("JPsiToP/Incoherent/hPhi2"))->Fill(RecoDecay::phi(daughter2)); if (trkDaughter1.hasTPC()) { JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsPt"))->Fill(trkDaughter1.pt(), trkDaughter1.tpcSignal()); JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter1), trkDaughter1.tpcSignal()); - JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter1[0], daughter1[1]), trkDaughter1.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter1), trkDaughter1.tpcSignal()); } if (trkDaughter2.hasTPC()) { JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsPt"))->Fill(trkDaughter2.pt(), trkDaughter2.tpcSignal()); JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsEta"))->Fill(RecoDecay::eta(daughter2), trkDaughter2.tpcSignal()); - JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter2[0], daughter2[1]), trkDaughter2.tpcSignal()); + JPsiToP.get(HIST("JPsiToP/Incoherent/PID/hTPCVsPhi"))->Fill(RecoDecay::phi(daughter2), trkDaughter2.tpcSignal()); } // fill J/psi histos JPsiToP.get(HIST("JPsiToP/Incoherent/hPt"))->Fill(RecoDecay::pt(mother)); - JPsiToP.get(HIST("JPsiToP/Incoherent/hEta"))->Fill(RecoDecay::eta(std::array{mother[0], mother[1], mother[2]})); - JPsiToP.get(HIST("JPsiToP/Incoherent/hPhi"))->Fill(RecoDecay::phi(mother[0], mother[1])); - JPsiToP.get(HIST("JPsiToP/Incoherent/hRap"))->Fill(RecoDecay::y(std::array{mother[0], mother[1], mother[2]}, mother[3])); - JPsiToP.get(HIST("JPsiToP/Incoherent/hIVM"))->Fill(RecoDecay::m(std::array{mother[0], mother[1], mother[2]}, energy)); + JPsiToP.get(HIST("JPsiToP/Incoherent/hEta"))->Fill(RecoDecay::eta(mother)); + JPsiToP.get(HIST("JPsiToP/Incoherent/hPhi"))->Fill(RecoDecay::phi(mother)); + JPsiToP.get(HIST("JPsiToP/Incoherent/hRap"))->Fill(rapJpsi); + JPsiToP.get(HIST("JPsiToP/Incoherent/hIVM"))->Fill(massJpsi); } } // end protons } // end two tracks From 91096c3b61a284e743c06d627de4edf6e3406911 Mon Sep 17 00:00:00 2001 From: Sara Haidlova Date: Wed, 10 Apr 2024 20:57:10 +0200 Subject: [PATCH 6/6] PWGUD: use of PhysicsConstants.h for masses --- PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx | 21 ++++----------------- 1 file changed, 4 insertions(+), 17 deletions(-) diff --git a/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx b/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx index d911ad90dc5..f95fbf4d286 100644 --- a/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx +++ b/PWGUD/Tasks/upcJpsiCentralBarrelCorr.cxx @@ -17,6 +17,7 @@ #include "Framework/runDataProcessing.h" #include "Framework/AnalysisTask.h" #include "CommonConstants/MathConstants.h" +#include "CommonConstants/PhysicsConstants.h" // O2Physics headers #include "Common/DataModel/PIDResponse.h" @@ -26,7 +27,6 @@ #include "PWGUD/Core/UPCJpsiCentralBarrelCorrHelper.h" // ROOT headers -#include "TDatabasePDG.h" #include "TLorentzVector.h" using namespace o2; @@ -59,8 +59,6 @@ struct upcJpsiCentralBarrel { Configurable TPCNClsCrossedRowsCut{"TPCNClsCrossedRowsCut", 70, "minimal number of crossed TPC rows"}; Configurable TPCChi2NCls{"TPCChi2NCls", 4, "minimal Chi2/cluster for the TPC track"}; - TDatabasePDG* pdg = nullptr; - // initialize histogram registry HistogramRegistry Statistics{ "Statistics", @@ -125,7 +123,6 @@ struct upcJpsiCentralBarrel { void init(InitContext&) { - pdg = TDatabasePDG::Instance(); const AxisSpec axisIVM{IVMAxis, "IVM axis"}; const AxisSpec axispt{ptAxis, "pt axis"}; @@ -365,16 +362,6 @@ struct upcJpsiCentralBarrel { Asymmetry.add("Asymmetry/Electron/Incoherent/DeltaPhi", "DeltaPhi", {HistType::kTH1F, {{180, -PI, PI}}}); } - float particleMass(TDatabasePDG* pdg, int pid) - { - auto mass = 0.; - TParticlePDG* pdgparticle = pdg->GetParticle(pid); - if (pdgparticle != nullptr) { - mass = pdgparticle->Mass(); - } - return mass; - } - template bool GoodTrackCuts(T const& track) { @@ -521,9 +508,9 @@ struct upcJpsiCentralBarrel { Statistics.get(HIST("Statistics/hNumberGTmu"))->Fill(countGTmu); Statistics.get(HIST("Statistics/hNumberGTp"))->Fill(countGTp); - float massEl = 0.000510998950; - float massMu = 0.10565837; - float massPr = 0.93827208816; + float massEl = o2::constants::physics::MassElectron; + float massMu = o2::constants::physics::MassMuonMinus; + float massPr = o2::constants::physics::MassProton; if (countGT == 2) { TLorentzVector mom, daughter[2];