diff --git a/PWGDQ/Core/VarManager.cxx b/PWGDQ/Core/VarManager.cxx index 0ae36125960..6f1c204152e 100644 --- a/PWGDQ/Core/VarManager.cxx +++ b/PWGDQ/Core/VarManager.cxx @@ -69,7 +69,7 @@ o2::vertexing::FwdDCAFitterN<2> VarManager::fgFitterTwoProngFwd; o2::vertexing::FwdDCAFitterN<3> VarManager::fgFitterThreeProngFwd; o2::globaltracking::MatchGlobalFwd VarManager::mMatching; std::map VarManager::fgCalibs; -bool VarManager::fgRunTPCPostCalibration[4] = {false, false, false, false}; +std::array VarManager::fgRunTPCPostCalibration = {false, false, false, false}; int VarManager::fgCalibrationType = 0; // 0 - no calibration, 1 - calibration vs (TPCncls,pIN,eta) typically for pp, 2 - calibration vs (eta,nPV,nLong,tLong) typically for PbPb bool VarManager::fgUseInterpolatedCalibration = true; // use interpolated calibration histograms (default: true) int VarManager::fgEfficiencyType = 0; // type of efficiency to be applied, default is no efficiency diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index 7079269d946..9bcbd1eb352 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -84,8 +84,6 @@ #include #include -#include - using SMatrix55 = ROOT::Math::SMatrix>; using SMatrix5 = ROOT::Math::SVector; using Vec3D = ROOT::Math::SVector; @@ -179,10 +177,10 @@ class VarManager : public TObject enum BarrelTrackFilteringBits { kIsConversionLeg = 0, // electron from conversions - kIsK0sLeg, // pion from K0s - kIsLambdaLeg, // proton or pion from Lambda - kIsALambdaLeg, // proton or pion from anti-Lambda - kIsOmegaLeg, // kaon from Omega baryon decay + kIsK0sLeg = 1, // pion from K0s + kIsLambdaLeg = 2, // proton or pion from Lambda + kIsALambdaLeg = 3, // proton or pion from anti-Lambda + kIsOmegaLeg = 4, // kaon from Omega baryon decay kDalitzBits = 5, // first bit for Dalitz tagged tracks kBarrelUserCutsBits = 13, // first bit for user track cuts kIsTPCPostcalibrated = 63 // tracks were postcalibrated for the TPC PID @@ -193,7 +191,7 @@ class VarManager : public TObject kMuonIsPropagated = 7 // whether the muon was propagated already }; - public: + // NOLINTNEXTLINE(readability-enum-initial-value) enum Variables { kNothing = -1, // Run wise variables @@ -1215,9 +1213,9 @@ class VarManager : public TObject } SetVariableDependencies(); } - static void SetUseVars(const std::vector usedVars) + static void SetUseVars(const std::vector& usedVars) { - for (auto& var : usedVars) { + for (auto const& var : usedVars) { fgUsedVars[var] = true; } } @@ -1357,17 +1355,9 @@ class VarManager : public TObject return (1.0 / harm) * TMath::ATan2(qnya, qnxa); }; - static float getDeltaPsiInRange(float psi1, float psi2, float harmonic) + static float getDeltaPsiInRange(float psi1, float psi2, unsigned int harmonic) { - float deltaPsi = psi1 - psi2; - if (std::abs(deltaPsi) > o2::constants::math::PI / harmonic) { - if (deltaPsi > 0.) { - deltaPsi -= o2::constants::math::TwoPI / harmonic; - } else { - deltaPsi += o2::constants::math::TwoPI / harmonic; - } - } - return deltaPsi; + return RecoDecay::constrainAngle(psi1 - psi2, -o2::constants::math::PI / harmonic, harmonic); } template static o2::dataformats::VertexBase RecalculatePrimaryVertex(T const& track0, T const& track1, const T1& collision); @@ -1396,17 +1386,17 @@ class VarManager : public TObject template static void FillEventFlowResoFactor(T const& hs_sp, T const& hs_ep, float* values = nullptr); template - static void FillTwoEvents(T const& event1, T const& event2, float* values = nullptr); + static void FillTwoEvents(T const& ev1, T const& ev2, float* values = nullptr); template - static void FillTwoMixEvents(T1 const& event1, T1 const& event2, T2 const& tracks1, T2 const& tracks2, float* values = nullptr); + static void FillTwoMixEvents(T1 const& ev1, T1 const& ev2, T2 const& tracks1, T2 const& tracks2, float* values = nullptr); template static void FillTwoMixEventsFlowResoFactor(T const& hs_sp, T const& hs_ep, float* values = nullptr); template - static void FillTwoMixEventsCumulants(T const& h_v22m, T const& h_v24m, T const& h_v22p, T const& h_v24p, T1 const& t1, T2 const& t2, float* values = nullptr); + static void FillTwoMixEventsCumulants(T const& h_v22ev1, T const& h_v24ev1, T const& h_v22ev2, T const& h_v24ev2, T1 const& t1, T2 const& t2, float* values = nullptr); template static void FillTrack(T const& track, float* values = nullptr); template - static void FillPhoton(T const& photon, float* values = nullptr); + static void FillPhoton(T const& track, float* values = nullptr); template static void FillTrackCollision(T const& track, C const& collision, float* values = nullptr); template @@ -1444,11 +1434,11 @@ class VarManager : public TObject template static void FillTripleMC(T1 const& t1, T2 const& t2, T3 const& t3, float* values = nullptr); template - static void FillQuadMC(T1 const& t1, T2 const& t2, T2 const& t3, float* values = nullptr); + static void FillQuadMC(T1 const& dilepton, T2 const& track1, T2 const& track2, float* values = nullptr); template static void FillPairVertexing(C const& collision, T const& t1, T const& t2, bool propToSV = false, float* values = nullptr); template - static void FillPairVertexingRecomputePV(C const& /*collision*/, T const& t1, T const& t2, o2::dataformats::VertexBase pvRefitted, float* values = nullptr); + static void FillPairVertexingRecomputePV(C const& /*collision*/, T const& t1, T const& t2, const o2::dataformats::VertexBase& pvRefitted, float* values = nullptr); template static void FillTripletVertexing(C const& collision, T const& t1, T const& t2, T const& t3, PairCandidateType tripletType, float* values = nullptr); template @@ -1500,19 +1490,19 @@ class VarManager : public TObject { fgCalibs[calib] = obj; // Check whether all the needed objects for TPC postcalibration are available - if (fgCalibs.find(kTPCElectronMean) != fgCalibs.end() && fgCalibs.find(kTPCElectronSigma) != fgCalibs.end()) { + if (fgCalibs.contains(kTPCElectronMean) && fgCalibs.contains(kTPCElectronSigma)) { fgRunTPCPostCalibration[0] = true; fgUsedVars[kTPCnSigmaEl_Corr] = true; } - if (fgCalibs.find(kTPCPionMean) != fgCalibs.end() && fgCalibs.find(kTPCPionSigma) != fgCalibs.end()) { + if (fgCalibs.contains(kTPCPionMean) && fgCalibs.contains(kTPCPionSigma)) { fgRunTPCPostCalibration[1] = true; fgUsedVars[kTPCnSigmaPi_Corr] = true; } - if (fgCalibs.find(kTPCKaonMean) != fgCalibs.end() && fgCalibs.find(kTPCKaonSigma) != fgCalibs.end()) { + if (fgCalibs.contains(kTPCKaonMean) && fgCalibs.contains(kTPCKaonSigma)) { fgRunTPCPostCalibration[2] = true; fgUsedVars[kTPCnSigmaKa_Corr] = true; } - if (fgCalibs.find(kTPCProtonMean) != fgCalibs.end() && fgCalibs.find(kTPCProtonSigma) != fgCalibs.end()) { + if (fgCalibs.contains(kTPCProtonMean) && fgCalibs.contains(kTPCProtonSigma)) { fgRunTPCPostCalibration[3] = true; fgUsedVars[kTPCnSigmaPr_Corr] = true; } @@ -1534,10 +1524,9 @@ class VarManager : public TObject { auto obj = fgCalibs.find(calib); if (obj == fgCalibs.end()) { - return 0x0; - } else { - return obj->second; + return nullptr; } + return obj->second; } static void SetTPCInterSectorBoundary(float boundarySize) { @@ -1557,7 +1546,6 @@ class VarManager : public TObject fgEOR = eor; } - public: VarManager(); ~VarManager() override; @@ -1599,7 +1587,7 @@ class VarManager : public TObject template static float calculatePhiV(const T1& t1, const T2& t2); template - static float LorentzTransformJpsihadroncosChi(TString Option, const T1& v1, const T2& v2); + static float LorentzTransformJpsihadroncosChi(const TString& Option, const T1& v1, const T2& v2); static o2::vertexing::DCAFitterN<2> fgFitterTwoProngBarrel; static o2::vertexing::DCAFitterN<3> fgFitterThreeProngBarrel; @@ -1609,7 +1597,7 @@ class VarManager : public TObject static o2::globaltracking::MatchGlobalFwd mMatching; static std::map fgCalibs; // map of calibration histograms - static bool fgRunTPCPostCalibration[4]; // 0-electron, 1-pion, 2-kaon, 3-proton + static std::array fgRunTPCPostCalibration; // 0-electron, 1-pion, 2-kaon, 3-proton static int fgCalibrationType; // 0 - no calibration, 1 - calibration vs (TPCncls,pIN,eta) typically for pp, 2 - calibration vs (eta,nPV,nLong,tLong) typically for PbPb static bool fgUseInterpolatedCalibration; // use interpolated calibration histograms (default: true) @@ -1657,21 +1645,18 @@ KFPTrack VarManager::createKFPTrackFromTrack(const T& track) track.cSnpSnp(), track.cTglY(), track.cTglZ(), track.cTglSnp(), track.cTglTgl(), track.c1PtY(), track.c1PtZ(), track.c1PtSnp(), track.c1PtTgl(), track.c1Pt21Pt2()}; o2::track::TrackParametrizationWithError trackparCov{track.x(), track.alpha(), std::move(trackpars), std::move(trackcovs)}; - std::array trkpos_par; - std::array trkmom_par; - std::array trk_cov; + std::array trkpos_par{}; + std::array trkmom_par{}; + std::array trk_cov{}; trackparCov.getXYZGlo(trkpos_par); trackparCov.getPxPyPzGlo(trkmom_par); trackparCov.getCovXYZPxPyPzGlo(trk_cov); - float trkpar_KF[6] = {trkpos_par[0], trkpos_par[1], trkpos_par[2], - trkmom_par[0], trkmom_par[1], trkmom_par[2]}; - float trkcov_KF[21]; - for (int i = 0; i < 21; i++) { - trkcov_KF[i] = trk_cov[i]; - } + std::array trkpar_KF = {trkpos_par[0], trkpos_par[1], trkpos_par[2], + trkmom_par[0], trkmom_par[1], trkmom_par[2]}; + std::array trkcov_KF = trk_cov; KFPTrack kfpTrack; - kfpTrack.SetParameters(trkpar_KF); - kfpTrack.SetCovarianceMatrix(trkcov_KF); + kfpTrack.SetParameters(trkpar_KF.data()); + kfpTrack.SetCovarianceMatrix(trkcov_KF.data()); kfpTrack.SetCharge(track.sign()); kfpTrack.SetNDF(track.tpcNClsFound() - 5); kfpTrack.SetChi2(track.tpcChi2NCl() * track.tpcNClsFound()); @@ -1683,17 +1668,14 @@ KFPTrack VarManager::createKFPFwdTrackFromFwdTrack(const T& muon) { o2::track::TrackParCovFwd trackparCov = FwdToTrackPar(muon, muon); - std::array trk_cov; + std::array trk_cov{}; trackparCov.getCovXYZPxPyPzGlo(trk_cov); - double trkpar_KF[6] = {trackparCov.getX(), trackparCov.getY(), trackparCov.getZ(), - trackparCov.getPx(), trackparCov.getPy(), trackparCov.getPz()}; - float trkcov_KF[21]; - for (int i = 0; i < 21; i++) { - trkcov_KF[i] = trk_cov[i]; - } + std::array trkpar_KF = {trackparCov.getX(), trackparCov.getY(), trackparCov.getZ(), + trackparCov.getPx(), trackparCov.getPy(), trackparCov.getPz()}; + std::array trkcov_KF = trk_cov; KFPTrack kfpTrack; - kfpTrack.SetParameters(trkpar_KF); - kfpTrack.SetCovarianceMatrix(trkcov_KF); + kfpTrack.SetParameters(trkpar_KF.data()); + kfpTrack.SetCovarianceMatrix(trkcov_KF.data()); kfpTrack.SetCharge(muon.sign()); kfpTrack.SetNDF(muon.nClusters() - 5); kfpTrack.SetChi2(muon.chi2()); @@ -1746,7 +1728,7 @@ o2::dataformats::VertexBase VarManager::RecalculatePrimaryVertex(T const& track0 o2::math_utils::Point3D vtxXYZ(kVtx.Parameter(0), kVtx.Parameter(1), kVtx.Parameter(2)); std::array vtxCov{kVtx.Covariance(0), kVtx.Covariance(1), kVtx.Covariance(2), kVtx.Covariance(3), kVtx.Covariance(4), kVtx.Covariance(5)}; - o2::dataformats::VertexBase primaryVertexRec = {std::move(vtxXYZ), std::move(vtxCov)}; + o2::dataformats::VertexBase primaryVertexRec = {vtxXYZ, vtxCov}; return primaryVertexRec; } @@ -1854,7 +1836,7 @@ void VarManager::FillPropagateMuon(const T& muon, const C& collision, float* val values[kMuonRAtAbsorberEnd] = std::sqrt(xAbs * xAbs + yAbs * yAbs); } - SMatrix55 cov = propmuon.getCovariances(); + const SMatrix55& cov = propmuon.getCovariances(); values[kMuonCXX] = cov(0, 0); values[kMuonCXY] = cov(1, 0); values[kMuonCYY] = cov(1, 1); @@ -1881,9 +1863,9 @@ void VarManager::FillGlobalMuonRefit(T1 const& muontrack, T2 const& mfttrack, co } if constexpr ((fillMap & MuonCov) > 0 || (fillMap & ReducedMuonCov) > 0) { o2::dataformats::GlobalFwdTrack propmuon = PropagateMuon(muontrack, collision); - double px = propmuon.getP() * sin(M_PI / 2 - atan(mfttrack.tgl())) * cos(mfttrack.phi()); - double py = propmuon.getP() * sin(M_PI / 2 - atan(mfttrack.tgl())) * sin(mfttrack.phi()); - double pz = propmuon.getP() * cos(M_PI / 2 - atan(mfttrack.tgl())); + double px = propmuon.getP() * std::sin(o2::constants::math::PIHalf - std::atan(mfttrack.tgl())) * std::cos(mfttrack.phi()); + double py = propmuon.getP() * std::sin(o2::constants::math::PIHalf - std::atan(mfttrack.tgl())) * std::sin(mfttrack.phi()); + double pz = propmuon.getP() * std::cos(o2::constants::math::PIHalf - std::atan(mfttrack.tgl())); double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2)); auto mftprop = o2::aod::fwdtrackutils::getTrackParCovFwdShift(mfttrack, fgzShiftFwd); values[kX] = mftprop.getX(); @@ -2002,7 +1984,7 @@ void VarManager::FillEvent(T const& event, float* values) values[kIsNoTFBorder] = event.selection_bit(o2::aod::evsel::kNoTimeFrameBorder); } if (fgUsedVars[kIsTriggerZNAZNC]) { - values[kIsTriggerZNAZNC] = event.selection_bit(o2::aod::evsel::kIsBBZNA) && event.selection_bit(o2::aod::evsel::kIsBBZNC); + values[kIsTriggerZNAZNC] = static_cast(event.selection_bit(o2::aod::evsel::kIsBBZNA) && event.selection_bit(o2::aod::evsel::kIsBBZNC)); } if (fgUsedVars[kIsNoSameBunch]) { values[kIsNoSameBunch] = event.selection_bit(o2::aod::evsel::kNoSameBunchPileup); @@ -2017,7 +1999,7 @@ void VarManager::FillEvent(T const& event, float* values) values[kIsVertexTOFmatched] = event.selection_bit(o2::aod::evsel::kIsVertexTOFmatched); } if (fgUsedVars[kIsSel8]) { - values[kIsSel8] = event.selection_bit(o2::aod::evsel::kIsTriggerTVX) && event.selection_bit(o2::aod::evsel::kNoITSROFrameBorder) && event.selection_bit(o2::aod::evsel::kNoTimeFrameBorder); + values[kIsSel8] = static_cast(event.selection_bit(o2::aod::evsel::kIsTriggerTVX) && event.selection_bit(o2::aod::evsel::kNoITSROFrameBorder) && event.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)); } if (fgUsedVars[kIsGoodITSLayer3]) { values[kIsGoodITSLayer3] = event.selection_bit(o2::aod::evsel::kIsGoodITSLayer3); @@ -2029,37 +2011,37 @@ void VarManager::FillEvent(T const& event, float* values) values[kIsGoodITSLayersAll] = event.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll); } if (fgUsedVars[kIsINT7]) { - values[kIsINT7] = (event.alias_bit(kINT7) > 0); + values[kIsINT7] = static_cast(event.alias_bit(kINT7) > 0); } if (fgUsedVars[kIsEMC7]) { - values[kIsEMC7] = (event.alias_bit(kEMC7) > 0); + values[kIsEMC7] = static_cast(event.alias_bit(kEMC7) > 0); } if (fgUsedVars[kIsINT7inMUON]) { - values[kIsINT7inMUON] = (event.alias_bit(kINT7inMUON) > 0); + values[kIsINT7inMUON] = static_cast(event.alias_bit(kINT7inMUON) > 0); } if (fgUsedVars[kIsMuonSingleLowPt7]) { - values[kIsMuonSingleLowPt7] = (event.alias_bit(kMuonSingleLowPt7) > 0); + values[kIsMuonSingleLowPt7] = static_cast(event.alias_bit(kMuonSingleLowPt7) > 0); } if (fgUsedVars[kIsMuonSingleHighPt7]) { - values[kIsMuonSingleHighPt7] = (event.alias_bit(kMuonSingleHighPt7) > 0); + values[kIsMuonSingleHighPt7] = static_cast(event.alias_bit(kMuonSingleHighPt7) > 0); } if (fgUsedVars[kIsMuonUnlikeLowPt7]) { - values[kIsMuonUnlikeLowPt7] = (event.alias_bit(kMuonUnlikeLowPt7) > 0); + values[kIsMuonUnlikeLowPt7] = static_cast(event.alias_bit(kMuonUnlikeLowPt7) > 0); } if (fgUsedVars[kIsMuonLikeLowPt7]) { - values[kIsMuonLikeLowPt7] = (event.alias_bit(kMuonLikeLowPt7) > 0); + values[kIsMuonLikeLowPt7] = static_cast(event.alias_bit(kMuonLikeLowPt7) > 0); } if (fgUsedVars[kIsCUP8]) { - values[kIsCUP8] = (event.alias_bit(kCUP8) > 0); + values[kIsCUP8] = static_cast(event.alias_bit(kCUP8) > 0); } if (fgUsedVars[kIsCUP9]) { - values[kIsCUP9] = (event.alias_bit(kCUP9) > 0); + values[kIsCUP9] = static_cast(event.alias_bit(kCUP9) > 0); } if (fgUsedVars[kIsMUP10]) { - values[kIsMUP10] = (event.alias_bit(kMUP10) > 0); + values[kIsMUP10] = static_cast(event.alias_bit(kMUP10) > 0); } if (fgUsedVars[kIsMUP11]) { - values[kIsMUP11] = (event.alias_bit(kMUP11) > 0); + values[kIsMUP11] = static_cast(event.alias_bit(kMUP11) > 0); } values[kVtxX] = event.posX(); values[kVtxY] = event.posY(); @@ -2140,14 +2122,14 @@ void VarManager::FillEvent(T const& event, float* values) values[kVtxZ] = event.posZ(); values[kVtxNcontrib] = event.numContrib(); if (fgUsedVars[kIsDoubleGap] || fgUsedVars[kIsSingleGap] || fgUsedVars[kIsSingleGapA] || fgUsedVars[kIsSingleGapC] || fgUsedVars[kIsNoGap]) { - values[kIsDoubleGap] = (event.tag_bit(56 + kDoubleGap) > 0); - values[kIsSingleGapA] = (event.tag_bit(56 + kSingleGapA) > 0); - values[kIsSingleGapC] = (event.tag_bit(56 + kSingleGapC) > 0); - values[kIsSingleGap] = values[kIsSingleGapA] || values[kIsSingleGapC]; - values[kIsNoGap] = !values[kIsDoubleGap] && !values[kIsSingleGap]; + values[kIsDoubleGap] = static_cast(event.tag_bit(56 + kDoubleGap) > 0); + values[kIsSingleGapA] = static_cast(event.tag_bit(56 + kSingleGapA) > 0); + values[kIsSingleGapC] = static_cast(event.tag_bit(56 + kSingleGapC) > 0); + values[kIsSingleGap] = static_cast(values[kIsSingleGapA] != 0.f || values[kIsSingleGapC] != 0.f); + values[kIsNoGap] = static_cast(values[kIsDoubleGap] == 0.f && values[kIsSingleGap] == 0.f); } if (fgUsedVars[kIsITSUPCMode]) { - values[kIsITSUPCMode] = (event.tag_bit(56 + kITSUPCMode) > 0); + values[kIsITSUPCMode] = static_cast(event.tag_bit(56 + kITSUPCMode) > 0); } values[kCollisionTime] = event.collisionTime(); values[kCollisionTimeRes] = event.collisionTimeRes(); @@ -2165,34 +2147,34 @@ void VarManager::FillEvent(T const& event, float* values) values[kIsNoITSROFBorderRecomputed] = bcInITSROF > fgITSROFBorderMarginLow && bcInITSROF < fgITSROFlength - fgITSROFBorderMarginHigh ? 1.0 : 0.0; } if (fgUsedVars[kIsNoITSROFBorder]) { - values[kIsNoITSROFBorder] = (event.selection_bit(o2::aod::evsel::kNoITSROFrameBorder) > 0); + values[kIsNoITSROFBorder] = static_cast(event.selection_bit(o2::aod::evsel::kNoITSROFrameBorder) > 0); } if (fgUsedVars[kIsTVXTriggered]) { - values[kIsTVXTriggered] = (event.selection_bit(o2::aod::evsel::kIsTriggerTVX) > 0); + values[kIsTVXTriggered] = static_cast(event.selection_bit(o2::aod::evsel::kIsTriggerTVX) > 0); } if (fgUsedVars[kIsNoTFBorder]) { - values[kIsNoTFBorder] = (event.selection_bit(o2::aod::evsel::kNoTimeFrameBorder) > 0); + values[kIsNoTFBorder] = static_cast(event.selection_bit(o2::aod::evsel::kNoTimeFrameBorder) > 0); } if (fgUsedVars[kNoCollInTimeRangeStandard]) { - values[kNoCollInTimeRangeStandard] = (event.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard) > 0); + values[kNoCollInTimeRangeStandard] = static_cast(event.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard) > 0); } if (fgUsedVars[kIsNoSameBunch]) { - values[kIsNoSameBunch] = (event.selection_bit(o2::aod::evsel::kNoSameBunchPileup) > 0); + values[kIsNoSameBunch] = static_cast(event.selection_bit(o2::aod::evsel::kNoSameBunchPileup) > 0); } if (fgUsedVars[kIsTriggerZNAZNC]) { - values[kIsTriggerZNAZNC] = event.selection_bit(o2::aod::evsel::kIsBBZNA) && event.selection_bit(o2::aod::evsel::kIsBBZNC); + values[kIsTriggerZNAZNC] = static_cast(event.selection_bit(o2::aod::evsel::kIsBBZNA) && event.selection_bit(o2::aod::evsel::kIsBBZNC)); } if (fgUsedVars[kIsGoodZvtxFT0vsPV]) { - values[kIsGoodZvtxFT0vsPV] = (event.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV) > 0); + values[kIsGoodZvtxFT0vsPV] = static_cast(event.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV) > 0); } if (fgUsedVars[kIsVertexITSTPC]) { - values[kIsVertexITSTPC] = (event.selection_bit(o2::aod::evsel::kIsVertexITSTPC) > 0); + values[kIsVertexITSTPC] = static_cast(event.selection_bit(o2::aod::evsel::kIsVertexITSTPC) > 0); } if (fgUsedVars[kIsVertexTOFmatched]) { - values[kIsVertexTOFmatched] = (event.selection_bit(o2::aod::evsel::kIsVertexTOFmatched) > 0); + values[kIsVertexTOFmatched] = static_cast(event.selection_bit(o2::aod::evsel::kIsVertexTOFmatched) > 0); } if (fgUsedVars[kIsSel8]) { - values[kIsSel8] = event.selection_bit(o2::aod::evsel::kIsTriggerTVX) && event.selection_bit(o2::aod::evsel::kNoTimeFrameBorder) && event.selection_bit(o2::aod::evsel::kNoITSROFrameBorder); + values[kIsSel8] = static_cast(event.selection_bit(o2::aod::evsel::kIsTriggerTVX) && event.selection_bit(o2::aod::evsel::kNoTimeFrameBorder) && event.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)); } if (fgUsedVars[kIsGoodITSLayer3]) { values[kIsGoodITSLayer3] = event.selection_bit(o2::aod::evsel::kIsGoodITSLayer3); @@ -2204,37 +2186,37 @@ void VarManager::FillEvent(T const& event, float* values) values[kIsGoodITSLayersAll] = event.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll); } if (fgUsedVars[kIsINT7]) { - values[kIsINT7] = (event.alias_bit(kINT7) > 0); + values[kIsINT7] = static_cast(event.alias_bit(kINT7) > 0); } if (fgUsedVars[kIsEMC7]) { - values[kIsEMC7] = (event.alias_bit(kEMC7) > 0); + values[kIsEMC7] = static_cast(event.alias_bit(kEMC7) > 0); } if (fgUsedVars[kIsINT7inMUON]) { - values[kIsINT7inMUON] = (event.alias_bit(kINT7inMUON) > 0); + values[kIsINT7inMUON] = static_cast(event.alias_bit(kINT7inMUON) > 0); } if (fgUsedVars[kIsMuonSingleLowPt7]) { - values[kIsMuonSingleLowPt7] = (event.alias_bit(kMuonSingleLowPt7) > 0); + values[kIsMuonSingleLowPt7] = static_cast(event.alias_bit(kMuonSingleLowPt7) > 0); } if (fgUsedVars[kIsMuonSingleHighPt7]) { - values[kIsMuonSingleHighPt7] = (event.alias_bit(kMuonSingleHighPt7) > 0); + values[kIsMuonSingleHighPt7] = static_cast(event.alias_bit(kMuonSingleHighPt7) > 0); } if (fgUsedVars[kIsMuonUnlikeLowPt7]) { - values[kIsMuonUnlikeLowPt7] = (event.alias_bit(kMuonUnlikeLowPt7) > 0); + values[kIsMuonUnlikeLowPt7] = static_cast(event.alias_bit(kMuonUnlikeLowPt7) > 0); } if (fgUsedVars[kIsMuonLikeLowPt7]) { - values[kIsMuonLikeLowPt7] = (event.alias_bit(kMuonLikeLowPt7) > 0); + values[kIsMuonLikeLowPt7] = static_cast(event.alias_bit(kMuonLikeLowPt7) > 0); } if (fgUsedVars[kIsCUP8]) { - values[kIsCUP8] = (event.alias_bit(kCUP8) > 0); + values[kIsCUP8] = static_cast(event.alias_bit(kCUP8) > 0); } if (fgUsedVars[kIsCUP9]) { - values[kIsCUP9] = (event.alias_bit(kCUP9) > 0); + values[kIsCUP9] = static_cast(event.alias_bit(kCUP9) > 0); } if (fgUsedVars[kIsMUP10]) { - values[kIsMUP10] = (event.alias_bit(kMUP10) > 0); + values[kIsMUP10] = static_cast(event.alias_bit(kMUP10) > 0); } if (fgUsedVars[kIsMUP11]) { - values[kIsMUP11] = (event.alias_bit(kMUP11) > 0); + values[kIsMUP11] = static_cast(event.alias_bit(kMUP11) > 0); } } @@ -2550,12 +2532,12 @@ void VarManager::FillEvent(T const& event, float* values) } if constexpr ((fillMap & EventFilter) > 0 || (fillMap & RapidityGapFilter) > 0) { - values[kIsDoubleGap] = (event.eventFilter() & (static_cast(1) << kDoubleGap)) > 0; - values[kIsSingleGapA] = (event.eventFilter() & (static_cast(1) << kSingleGapA)) > 0; - values[kIsSingleGapC] = (event.eventFilter() & (static_cast(1) << kSingleGapC)) > 0; - values[kIsSingleGap] = values[kIsSingleGapA] || values[kIsSingleGapC]; - values[kIsNoGap] = !values[kIsDoubleGap] && !values[kIsSingleGap]; - values[kIsITSUPCMode] = (event.eventFilter() & (static_cast(1) << kITSUPCMode)) > 0; + values[kIsDoubleGap] = static_cast((event.eventFilter() & (static_cast(1) << kDoubleGap)) > 0); + values[kIsSingleGapA] = static_cast((event.eventFilter() & (static_cast(1) << kSingleGapA)) > 0); + values[kIsSingleGapC] = static_cast((event.eventFilter() & (static_cast(1) << kSingleGapC)) > 0); + values[kIsSingleGap] = static_cast(values[kIsSingleGapA] != 0.f || values[kIsSingleGapC] != 0.f); + values[kIsNoGap] = static_cast(values[kIsDoubleGap] == 0.f && values[kIsSingleGap] == 0.f); + values[kIsITSUPCMode] = static_cast((event.eventFilter() & (static_cast(1) << kITSUPCMode)) > 0); } if constexpr ((fillMap & ReducedZdc) > 0) { @@ -2729,7 +2711,7 @@ void VarManager::FillEventTracks(T const& tracks, float* values) int counter2mm = 0; int counter5mm = 0; int counter10mm = 0; - for (auto& d : dcazValues) { + for (auto const& d : dcazValues) { double absD = std::abs(d); if (absD > 0.01) { counter100um++; @@ -2810,10 +2792,10 @@ void VarManager::FillTwoMixEventsCumulants(T const& h_v22ev1, T const& h_v24ev1, values = fgValues; } - int idx_v22ev1; - int idx_v24ev1; - int idx_v22ev2; - int idx_v24ev2; + int idx_v22ev1 = 0; + int idx_v24ev1 = 0; + int idx_v22ev2 = 0; + int idx_v24ev2 = 0; if (values[kTwoEvCentFT0C1] >= 0.) { if (t1.sign() < 0) { @@ -2884,8 +2866,8 @@ void VarManager::FillTwoMixEvents(T1 const& ev1, T1 const& ev2, T2 const& /*trac /* uint32_t Track1Filter = 0; uint32_t Track2Filter = 0; - for (auto& track1 : tracks1) { Track1Filter = uint32_t(track1.isMuonSelected());} - for (auto& track2 : tracks2) { Track2Filter = uint32_t(track2.isMuonSelected());} + for (auto const& track1 : tracks1) { Track1Filter = uint32_t(track1.isMuonSelected());} + for (auto const& track2 : tracks2) { Track2Filter = uint32_t(track2.isMuonSelected());} */ if constexpr ((fillMap & CollisionCent) > 0 || (fillMap & ReducedEventExtended) > 0) { values[kTwoEvCentFT0C1] = ev1.centFT0C(); @@ -2956,23 +2938,17 @@ void VarManager::FillTrack(T const& track, float* values) values[kPhi] = track.phi(); values[kCharge] = track.sign(); if (fgUsedVars[kPhiTPCOuter]) { - values[kPhiTPCOuter] = track.phi() - (track.sign() > 0 ? 1.0 : -1.0) * (TMath::PiOver2() - TMath::ACos(0.22 * fgMagField / track.pt())); - if (values[kPhiTPCOuter] > TMath::TwoPi()) { - values[kPhiTPCOuter] -= TMath::TwoPi(); - } - if (values[kPhiTPCOuter] < 0.0) { - values[kPhiTPCOuter] += TMath::TwoPi(); - } + values[kPhiTPCOuter] = RecoDecay::constrainAngle(track.phi() - (track.sign() > 0 ? 1.0 : -1.0) * (o2::constants::math::PIHalf - TMath::ACos(0.22 * fgMagField / track.pt()))); } if (fgUsedVars[kTrackIsInsideTPCModule]) { - float localSectorPhi = values[kPhiTPCOuter] - TMath::Floor(18.0 * values[kPhiTPCOuter] / TMath::TwoPi()) * (TMath::TwoPi() / 18.0); + float localSectorPhi = values[kPhiTPCOuter] - TMath::Floor(18.0 * values[kPhiTPCOuter] / o2::constants::math::TwoPI) * (o2::constants::math::TwoPI / 18.0); float edge = fgTPCInterSectorBoundary / 2.0 / 246.6; // minimal inter-sector boundary as angle float curvature = 3.0 * 3.33 * track.pt() / fgMagField * (1.0 - TMath::Sin(TMath::ACos(0.22 * fgMagField / track.pt()))); if (curvature / 2.466 > edge) { edge = curvature / 2.466; } double min = edge; - double max = TMath::TwoPi() / 18.0 - edge; + double max = o2::constants::math::TwoPI / 18.0 - edge; values[kTrackIsInsideTPCModule] = (localSectorPhi > min && localSectorPhi < max ? 1.0 : 0.0); } @@ -2987,7 +2963,7 @@ void VarManager::FillTrack(T const& track, float* values) values[kIsLegFromAntiLambda] = track.filteringFlags_bit(VarManager::kIsALambdaLeg); values[kIsLegFromOmega] = track.filteringFlags_bit(VarManager::kIsOmegaLeg); - values[kIsProtonFromLambdaAndAntiLambda] = static_cast((values[kIsLegFromLambda] * track.sign() > 0) || (values[kIsLegFromAntiLambda] * (-track.sign()) > 0)); + values[kIsProtonFromLambdaAndAntiLambda] = static_cast((values[kIsLegFromLambda] * track.sign() > 0) || (values[kIsLegFromAntiLambda] * (-track.sign()) > 0)); for (int i = 0; i < 8; i++) { values[kIsDalitzLeg + i] = track.filteringFlags_bit(VarManager::kDalitzBits + i); @@ -3024,44 +3000,44 @@ void VarManager::FillTrack(T const& track, float* values) values[kPin] = track.tpcInnerParam(); values[kSignedPin] = track.tpcInnerParam() * track.sign(); if (fgUsedVars[kIsITSrefit]) { - values[kIsITSrefit] = (track.flags() & o2::aod::track::ITSrefit) > 0; // NOTE: This is just for Run-2 + values[kIsITSrefit] = static_cast((track.flags() & o2::aod::track::ITSrefit) > 0); // NOTE: This is just for Run-2 } if (fgUsedVars[kTrackTimeResIsRange]) { - values[kTrackTimeResIsRange] = (track.flags() & o2::aod::track::TrackTimeResIsRange) > 0; // NOTE: This is NOT for Run-2 + values[kTrackTimeResIsRange] = static_cast((track.flags() & o2::aod::track::TrackTimeResIsRange) > 0); // NOTE: This is NOT for Run-2 } if (fgUsedVars[kIsTPCrefit]) { - values[kIsTPCrefit] = (track.flags() & o2::aod::track::TPCrefit) > 0; // NOTE: This is just for Run-2 + values[kIsTPCrefit] = static_cast((track.flags() & o2::aod::track::TPCrefit) > 0); // NOTE: This is just for Run-2 } if (fgUsedVars[kPVContributor]) { - values[kPVContributor] = (track.flags() & o2::aod::track::PVContributor) > 0; // NOTE: This is NOT for Run-2 + values[kPVContributor] = static_cast((track.flags() & o2::aod::track::PVContributor) > 0); // NOTE: This is NOT for Run-2 } if (fgUsedVars[kIsGoldenChi2]) { - values[kIsGoldenChi2] = (track.flags() & o2::aod::track::GoldenChi2) > 0; // NOTE: This is just for Run-2 + values[kIsGoldenChi2] = static_cast((track.flags() & o2::aod::track::GoldenChi2) > 0); // NOTE: This is just for Run-2 } if (fgUsedVars[kOrphanTrack]) { - values[kOrphanTrack] = (track.flags() & o2::aod::track::OrphanTrack) > 0; // NOTE: This is NOT for Run-2 + values[kOrphanTrack] = static_cast((track.flags() & o2::aod::track::OrphanTrack) > 0); // NOTE: This is NOT for Run-2 } if (fgUsedVars[kIsSPDfirst]) { - values[kIsSPDfirst] = (track.itsClusterMap() & uint8_t(1)) > 0; + values[kIsSPDfirst] = static_cast((track.itsClusterMap() & uint8_t(1)) > 0); } if (fgUsedVars[kIsSPDboth]) { - values[kIsSPDboth] = (track.itsClusterMap() & uint8_t(3)) > 0; + values[kIsSPDboth] = static_cast((track.itsClusterMap() & uint8_t(3)) > 0); } if (fgUsedVars[kIsSPDany]) { - values[kIsSPDany] = (track.itsClusterMap() & uint8_t(1)) || (track.itsClusterMap() & uint8_t(2)); + values[kIsSPDany] = static_cast((track.itsClusterMap() & uint8_t(1)) || (track.itsClusterMap() & uint8_t(2))); } if (fgUsedVars[kITSClusterMap]) { values[kITSClusterMap] = track.itsClusterMap(); } if (fgUsedVars[kIsITSibFirst]) { - values[kIsITSibFirst] = (track.itsClusterMap() & uint8_t(1)) > 0; + values[kIsITSibFirst] = static_cast((track.itsClusterMap() & uint8_t(1)) > 0); } if (fgUsedVars[kIsITSibAny]) { - values[kIsITSibAny] = (track.itsClusterMap() & (1 << uint8_t(0))) > 0 || (track.itsClusterMap() & (1 << uint8_t(1))) > 0 || (track.itsClusterMap() & (1 << uint8_t(2))) > 0; + values[kIsITSibAny] = static_cast((track.itsClusterMap() & (1 << uint8_t(0))) > 0 || (track.itsClusterMap() & (1 << uint8_t(1))) > 0 || (track.itsClusterMap() & (1 << uint8_t(2))) > 0); } if (fgUsedVars[kIsITSibAll]) { - values[kIsITSibAll] = (track.itsClusterMap() & (1 << uint8_t(0))) > 0 && (track.itsClusterMap() & (1 << uint8_t(1))) > 0 && (track.itsClusterMap() & (1 << uint8_t(2))) > 0; + values[kIsITSibAll] = static_cast((track.itsClusterMap() & (1 << uint8_t(0))) > 0 && (track.itsClusterMap() & (1 << uint8_t(1))) > 0 && (track.itsClusterMap() & (1 << uint8_t(2))) > 0); } values[kTrackTime] = track.trackTime(); @@ -3170,18 +3146,18 @@ void VarManager::FillTrack(T const& track, float* values) // Quantities based on the dalitz selections if constexpr ((fillMap & DalitzBits) > 0) { for (int i = 0; i < 8; i++) { - values[kIsDalitzLeg + i] = static_cast(track.dalitzBits() & (uint8_t(1) << i)); + values[kIsDalitzLeg + i] = static_cast(track.dalitzBits() & (uint8_t(1) << i)); } } // Quantities based on the V0 selections if constexpr ((fillMap & TrackV0Bits) > 0) { - values[kIsLegFromGamma] = static_cast(track.pidbit() & (uint8_t(1) << VarManager::kIsConversionLeg)); - values[kIsLegFromK0S] = static_cast(track.pidbit() & (uint8_t(1) << VarManager::kIsK0sLeg)); - values[kIsLegFromLambda] = static_cast(track.pidbit() & (uint8_t(1) << VarManager::kIsLambdaLeg)); - values[kIsLegFromAntiLambda] = static_cast(track.pidbit() & (uint8_t(1) << VarManager::kIsALambdaLeg)); - values[kIsLegFromOmega] = static_cast(track.pidbit() & (uint8_t(1) << VarManager::kIsOmegaLeg)); - values[kIsProtonFromLambdaAndAntiLambda] = static_cast((values[kIsLegFromLambda] * track.sign() > 0) || (values[kIsLegFromAntiLambda] * (-track.sign()) > 0)); + values[kIsLegFromGamma] = static_cast(track.pidbit() & (uint8_t(1) << VarManager::kIsConversionLeg)); + values[kIsLegFromK0S] = static_cast(track.pidbit() & (uint8_t(1) << VarManager::kIsK0sLeg)); + values[kIsLegFromLambda] = static_cast(track.pidbit() & (uint8_t(1) << VarManager::kIsLambdaLeg)); + values[kIsLegFromAntiLambda] = static_cast(track.pidbit() & (uint8_t(1) << VarManager::kIsALambdaLeg)); + values[kIsLegFromOmega] = static_cast(track.pidbit() & (uint8_t(1) << VarManager::kIsOmegaLeg)); + values[kIsProtonFromLambdaAndAntiLambda] = static_cast((values[kIsLegFromLambda] * track.sign() > 0) || (values[kIsLegFromAntiLambda] * (-track.sign()) > 0)); } // Quantities based on the barrel PID tables @@ -3483,13 +3459,13 @@ void VarManager::FillTrackCollisionMC(T1 const& track, T2 const& MotherTrack, C float m = o2::constants::physics::MassBPlus; float pdgLifetime = 1.638e-12; // s - if (std::abs(MotherTrack.pdgCode()) == 511) { + if (std::abs(MotherTrack.pdgCode()) == o2::constants::physics::Pdg::kB0) { m = o2::constants::physics::MassB0; pdgLifetime = 1.517e-12; // s } // Extract the collision primary vertex position using constexpr, since the collision type may be CollisionMC or ReducedMCEvent - double collPos[3] = {0.0, 0.0, 0.0}; + std::array collPos = {0.0, 0.0, 0.0}; if constexpr (fillMap & ObjTypes::CollisionMC) { collPos[0] = collision.posX(); collPos[1] = collision.posY(); @@ -3536,8 +3512,9 @@ void VarManager::FillTrackCollisionMC(T1 const& track, const std::array void VarManager::FillEnergyCorrelatorsMC(T const& track, T1 const& t1, float* values, float Translow, float Transhigh, float Accweight) { // energy correlators - float MassHadron; + float MassHadron = o2::constants::physics::MassPionCharged; if constexpr (pairType == kJpsiHadronMass) { MassHadron = TMath::Sqrt(t1.e() * t1.e() - t1.p() * t1.p()); } - if constexpr (pairType == kJpsiPionMass) { - MassHadron = o2::constants::physics::MassPionCharged; - } ROOT::Math::PtEtaPhiMVector v1(track.pt(), track.eta(), track.phi(), o2::constants::physics::MassJPsi); float deltaphi = RecoDecay::constrainAngle(track.phi() - t1.phi(), -o2::constants::math::PIHalf); float deltaeta = t1.eta() - track.eta(); @@ -3602,7 +3576,7 @@ void VarManager::FillEnergyCorrelatorsMC(T const& track, T1 const& t1, float* va float randomPhi_away = -o2::constants::math::PIHalf; float deltaphitrans = RecoDecay::constrainAngle(track.phi() - t1.phi(), -o2::constants::math::PI); - if ((deltaphitrans > -Transhigh * TMath::Pi() && deltaphitrans < -Translow * TMath::Pi()) || (deltaphitrans > Translow * TMath::Pi() && deltaphitrans < Transhigh * TMath::Pi())) { + if ((deltaphitrans > -Transhigh * o2::constants::math::PI && deltaphitrans < -Translow * o2::constants::math::PI) || (deltaphitrans > Translow * o2::constants::math::PI && deltaphitrans < Transhigh * o2::constants::math::PI)) { randomPhi_trans = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf); randomPhi_toward = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf); randomPhi_away = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf); @@ -3692,7 +3666,7 @@ void VarManager::FillPair(T1 const& t1, T2 const& t2, float* values) values[kPt] = v12.Pt(); values[kEta] = v12.Eta(); // values[kPhi] = v12.Phi(); - values[kPhi] = v12.Phi() > 0 ? v12.Phi() : v12.Phi() + 2. * M_PI; + values[kPhi] = RecoDecay::constrainAngle(v12.Phi()); values[kRap] = -v12.Rapidity(); double Ptot1 = TMath::Sqrt(v1.Px() * v1.Px() + v1.Py() * v1.Py() + v1.Pz() * v1.Pz()); double Ptot2 = TMath::Sqrt(v2.Px() * v2.Px() + v2.Py() * v2.Py() + v2.Pz() * v2.Pz()); @@ -3706,14 +3680,7 @@ void VarManager::FillPair(T1 const& t1, T2 const& t2, float* values) values[kPhi2] = t2.phi(); if (fgUsedVars[kDeltaPhiPair2]) { - double phipair2 = v1.Phi() - v2.Phi(); - if (phipair2 > 3 * TMath::Pi() / 2) { - values[kDeltaPhiPair2] = phipair2 - 2 * TMath::Pi(); - } else if (phipair2 < -TMath::Pi() / 2) { - values[kDeltaPhiPair2] = phipair2 + 2 * TMath::Pi(); - } else { - values[kDeltaPhiPair2] = phipair2; - } + values[kDeltaPhiPair2] = RecoDecay::constrainAngle(v1.Phi() - v2.Phi(), -o2::constants::math::PIHalf); } if (fgUsedVars[kDeltaEtaPair2]) { @@ -3733,10 +3700,12 @@ void VarManager::FillPair(T1 const& t1, T2 const& t2, float* values) values[kOpeningAngle] = 0.; } else { double arg = scalar / Ptot12; - if (arg > 1.) + if (arg > 1.) { arg = 1.; - if (arg < -1) + } + if (arg < -1) { arg = -1; + } values[kOpeningAngle] = TMath::ACos(arg); } } @@ -3761,26 +3730,18 @@ void VarManager::FillPair(T1 const& t1, T2 const& t2, float* values) ROOT::Math::XYZVectorF zaxis_HE{(v12.Vect()).Unit()}; ROOT::Math::XYZVectorF yaxis_HE{(Beam1_CM.Cross(Beam2_CM)).Unit()}; ROOT::Math::XYZVectorF xaxis_HE{(yaxis_HE.Cross(zaxis_HE)).Unit()}; - if (fgUsedVars[kCosThetaHE]) + if (fgUsedVars[kCosThetaHE]) { values[kCosThetaHE] = zaxis_HE.Dot(v_CM); + } if (fgUsedVars[kPhiHE]) { - values[kPhiHE] = TMath::ATan2(yaxis_HE.Dot(v_CM), xaxis_HE.Dot(v_CM)); - if (values[kPhiHE] < 0) { - values[kPhiHE] += 2 * TMath::Pi(); // ensure phi is in [0, 2pi] - } + values[kPhiHE] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_HE.Dot(v_CM), xaxis_HE.Dot(v_CM))); } if (fgUsedVars[kPhiTildeHE]) { if (fgUsedVars[kCosThetaHE] && fgUsedVars[kPhiHE]) { if (values[kCosThetaHE] > 0) { - values[kPhiTildeHE] = values[kPhiHE] - 0.25 * TMath::Pi(); // phi_tilde = phi - pi/4 - if (values[kPhiTildeHE] < 0) { - values[kPhiTildeHE] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kPhiTildeHE] = RecoDecay::constrainAngle(values[kPhiHE] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 } else { - values[kPhiTildeHE] = values[kPhiHE] - 0.75 * TMath::Pi(); // phi_tilde = phi - 3pi/4 - if (values[kPhiTildeHE] < 0) { - values[kPhiTildeHE] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kPhiTildeHE] = RecoDecay::constrainAngle(values[kPhiHE] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 } } else { values[kPhiTildeHE] = -999; // not computable @@ -3792,26 +3753,18 @@ void VarManager::FillPair(T1 const& t1, T2 const& t2, float* values) ROOT::Math::XYZVectorF zaxis_CS{(Beam1_CM - Beam2_CM).Unit()}; ROOT::Math::XYZVectorF yaxis_CS{(Beam1_CM.Cross(Beam2_CM)).Unit()}; ROOT::Math::XYZVectorF xaxis_CS{(yaxis_CS.Cross(zaxis_CS)).Unit()}; - if (fgUsedVars[kCosThetaCS]) + if (fgUsedVars[kCosThetaCS]) { values[kCosThetaCS] = zaxis_CS.Dot(v_CM); + } if (fgUsedVars[kPhiCS]) { - values[kPhiCS] = TMath::ATan2(yaxis_CS.Dot(v_CM), xaxis_CS.Dot(v_CM)); - if (values[kPhiCS] < 0) { - values[kPhiCS] += 2 * TMath::Pi(); // ensure phi is in [0, 2pi] - } + values[kPhiCS] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_CS.Dot(v_CM), xaxis_CS.Dot(v_CM))); } if (fgUsedVars[kPhiTildeCS]) { if (fgUsedVars[kCosThetaCS] && fgUsedVars[kPhiCS]) { if (values[kCosThetaCS] > 0) { - values[kPhiTildeCS] = values[kPhiCS] - 0.25 * TMath::Pi(); // phi_tilde = phi - pi/4 - if (values[kPhiTildeCS] < 0) { - values[kPhiTildeCS] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kPhiTildeCS] = RecoDecay::constrainAngle(values[kPhiCS] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 } else { - values[kPhiTildeCS] = values[kPhiCS] - 0.75 * TMath::Pi(); // phi_tilde = phi - 3pi/4 - if (values[kPhiTildeCS] < 0) { - values[kPhiTildeCS] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kPhiTildeCS] = RecoDecay::constrainAngle(values[kPhiCS] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 } } else { values[kPhiTildeCS] = -999; // not computable @@ -3827,23 +3780,14 @@ void VarManager::FillPair(T1 const& t1, T2 const& t2, float* values) values[kCosThetaPP] = zaxis_PP.Dot(v_CM) / std::sqrt(zaxis_PP.Mag2()); } if (fgUsedVars[kPhiPP]) { - values[kPhiPP] = TMath::ATan2(yaxis_PP.Dot(v_CM), xaxis_PP.Dot(v_CM)); - if (values[kPhiPP] < 0) { - values[kPhiPP] += 2 * TMath::Pi(); // ensure phi is in [0, 2pi] - } + values[kPhiPP] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_PP.Dot(v_CM), xaxis_PP.Dot(v_CM))); } if (fgUsedVars[kPhiTildePP]) { if (fgUsedVars[kCosThetaPP] && fgUsedVars[kPhiPP]) { if (values[kCosThetaPP] > 0) { - values[kPhiTildePP] = values[kPhiPP] - 0.25 * TMath::Pi(); // phi_tilde = phi - pi/4 - if (values[kPhiTildePP] < 0) { - values[kPhiTildePP] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kPhiTildePP] = RecoDecay::constrainAngle(values[kPhiPP] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 } else { - values[kPhiTildePP] = values[kPhiPP] - 0.75 * TMath::Pi(); // phi_tilde = phi - 3pi/4 - if (values[kPhiTildePP] < 0) { - values[kPhiTildePP] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kPhiTildePP] = RecoDecay::constrainAngle(values[kPhiPP] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 } } else { values[kPhiTildePP] = -999; // not computable @@ -3853,10 +3797,11 @@ void VarManager::FillPair(T1 const& t1, T2 const& t2, float* values) if (useRM) { double randomCostheta = gRandom->Uniform(-1., 1.); - double randomPhi = gRandom->Uniform(0., 2. * TMath::Pi()); + double randomPhi = gRandom->Uniform(0., o2::constants::math::TwoPI); ROOT::Math::XYZVectorF zaxis_RM(randomCostheta, std::sqrt(1 - randomCostheta * randomCostheta) * std::cos(randomPhi), std::sqrt(1 - randomCostheta * randomCostheta) * std::sin(randomPhi)); - if (fgUsedVars[kCosThetaRM]) + if (fgUsedVars[kCosThetaRM]) { values[kCosThetaRM] = zaxis_RM.Dot(v_CM); + } } } @@ -3981,11 +3926,8 @@ void VarManager::FillPairRotation(T1 const& t1, T2 const& t2, float* values) m2 = o2::constants::physics::MassMuon; } - double dphi = gRandom->Uniform(0., 2. * TMath::Pi()); - double rotationphi2 = t2.phi() + dphi; - - if (rotationphi2 > 2. * TMath::Pi()) - rotationphi2 -= 2. * TMath::Pi(); + double dphi = gRandom->Uniform(0., o2::constants::math::TwoPI); + double rotationphi2 = RecoDecay::constrainAngle(t2.phi() + dphi); values[kCharge] = t1.sign() + t2.sign(); values[kCharge1] = t1.sign(); @@ -3997,7 +3939,7 @@ void VarManager::FillPairRotation(T1 const& t1, T2 const& t2, float* values) values[kPt] = v12.Pt(); values[kEta] = v12.Eta(); // values[kPhi] = v12.Phi(); - values[kPhi] = v12.Phi() > 0 ? v12.Phi() : v12.Phi() + 2. * M_PI; + values[kPhi] = RecoDecay::constrainAngle(v12.Phi()); values[kRap] = -v12.Rapidity(); double Ptot1 = TMath::Sqrt(v1.Px() * v1.Px() + v1.Py() * v1.Py() + v1.Pz() * v1.Pz()); double Ptot2 = TMath::Sqrt(v2.Px() * v2.Px() + v2.Py() * v2.Py() + v2.Pz() * v2.Pz()); @@ -4226,7 +4168,7 @@ void VarManager::FillPairME(T1 const& t1, T2 const& t2, float* values) values[kPt] = v12.Pt(); values[kEta] = v12.Eta(); // values[kPhi] = v12.Phi(); - values[kPhi] = v12.Phi() > 0 ? v12.Phi() : v12.Phi() + 2. * M_PI; + values[kPhi] = RecoDecay::constrainAngle(v12.Phi()); values[kRap] = -v12.Rapidity(); // Per-track quantities so ME histograms can use kPt1/kPt2/kEta1/kEta2/kPhi1/kPhi2 just like SE FillPair does. @@ -4238,14 +4180,7 @@ void VarManager::FillPairME(T1 const& t1, T2 const& t2, float* values) values[kPhi2] = t2.phi(); if (fgUsedVars[kDeltaPhiPair2]) { - double phipair2ME = v1.Phi() - v2.Phi(); - if (phipair2ME > 3 * TMath::Pi() / 2) { - values[kDeltaPhiPair2] = phipair2ME - 2 * TMath::Pi(); - } else if (phipair2ME < -TMath::Pi() / 2) { - values[kDeltaPhiPair2] = phipair2ME + 2 * TMath::Pi(); - } else { - values[kDeltaPhiPair2] = phipair2ME; - } + values[kDeltaPhiPair2] = RecoDecay::constrainAngle(v1.Phi() - v2.Phi(), -o2::constants::math::PIHalf); } if (fgUsedVars[kDeltaEtaPair2]) { @@ -4272,26 +4207,18 @@ void VarManager::FillPairME(T1 const& t1, T2 const& t2, float* values) ROOT::Math::XYZVectorF zaxis_HE{(v12.Vect()).Unit()}; ROOT::Math::XYZVectorF yaxis_HE{(Beam1_CM.Cross(Beam2_CM)).Unit()}; ROOT::Math::XYZVectorF xaxis_HE{(yaxis_HE.Cross(zaxis_HE)).Unit()}; - if (fgUsedVars[kCosThetaHE]) + if (fgUsedVars[kCosThetaHE]) { values[kCosThetaHE] = zaxis_HE.Dot(v_CM); + } if (fgUsedVars[kPhiHE]) { - values[kPhiHE] = TMath::ATan2(yaxis_HE.Dot(v_CM), xaxis_HE.Dot(v_CM)); - if (values[kPhiHE] < 0) { - values[kPhiHE] += 2 * TMath::Pi(); // ensure phi is in [0, 2pi] - } + values[kPhiHE] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_HE.Dot(v_CM), xaxis_HE.Dot(v_CM))); } if (fgUsedVars[kPhiTildeHE]) { if (fgUsedVars[kCosThetaHE] && fgUsedVars[kPhiHE]) { if (values[kCosThetaHE] > 0) { - values[kPhiTildeHE] = values[kPhiHE] - 0.25 * TMath::Pi(); // phi_tilde = phi - pi/4 - if (values[kPhiTildeHE] < 0) { - values[kPhiTildeHE] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kPhiTildeHE] = RecoDecay::constrainAngle(values[kPhiHE] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 } else { - values[kPhiTildeHE] = values[kPhiHE] - 0.75 * TMath::Pi(); // phi_tilde = phi - 3pi/4 - if (values[kPhiTildeHE] < 0) { - values[kPhiTildeHE] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kPhiTildeHE] = RecoDecay::constrainAngle(values[kPhiHE] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 } } else { values[kPhiTildeHE] = -999; // not computable @@ -4303,26 +4230,18 @@ void VarManager::FillPairME(T1 const& t1, T2 const& t2, float* values) ROOT::Math::XYZVectorF zaxis_CS{(Beam1_CM - Beam2_CM).Unit()}; ROOT::Math::XYZVectorF yaxis_CS{(Beam1_CM.Cross(Beam2_CM)).Unit()}; ROOT::Math::XYZVectorF xaxis_CS{(yaxis_CS.Cross(zaxis_CS)).Unit()}; - if (fgUsedVars[kCosThetaCS]) + if (fgUsedVars[kCosThetaCS]) { values[kCosThetaCS] = zaxis_CS.Dot(v_CM); + } if (fgUsedVars[kPhiCS]) { - values[kPhiCS] = TMath::ATan2(yaxis_CS.Dot(v_CM), xaxis_CS.Dot(v_CM)); - if (values[kPhiCS] < 0) { - values[kPhiCS] += 2 * TMath::Pi(); // ensure phi is in [0, 2pi] - } + values[kPhiCS] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_CS.Dot(v_CM), xaxis_CS.Dot(v_CM))); } if (fgUsedVars[kPhiTildeCS]) { if (fgUsedVars[kCosThetaCS] && fgUsedVars[kPhiCS]) { if (values[kCosThetaCS] > 0) { - values[kPhiTildeCS] = values[kPhiCS] - 0.25 * TMath::Pi(); // phi_tilde = phi - pi/4 - if (values[kPhiTildeCS] < 0) { - values[kPhiTildeCS] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kPhiTildeCS] = RecoDecay::constrainAngle(values[kPhiCS] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 } else { - values[kPhiTildeCS] = values[kPhiCS] - 0.75 * TMath::Pi(); // phi_tilde = phi - 3pi/4 - if (values[kPhiTildeCS] < 0) { - values[kPhiTildeCS] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kPhiTildeCS] = RecoDecay::constrainAngle(values[kPhiCS] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 } } else { values[kPhiTildeCS] = -999; // not computable @@ -4338,23 +4257,14 @@ void VarManager::FillPairME(T1 const& t1, T2 const& t2, float* values) values[kCosThetaPP] = zaxis_PP.Dot(v_CM) / std::sqrt(zaxis_PP.Mag2()); } if (fgUsedVars[kPhiPP]) { - values[kPhiPP] = TMath::ATan2(yaxis_PP.Dot(v_CM), xaxis_PP.Dot(v_CM)); - if (values[kPhiPP] < 0) { - values[kPhiPP] += 2 * TMath::Pi(); // ensure phi is in [0, 2pi] - } + values[kPhiPP] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_PP.Dot(v_CM), xaxis_PP.Dot(v_CM))); } if (fgUsedVars[kPhiTildePP]) { if (fgUsedVars[kCosThetaPP] && fgUsedVars[kPhiPP]) { if (values[kCosThetaPP] > 0) { - values[kPhiTildePP] = values[kPhiPP] - 0.25 * TMath::Pi(); // phi_tilde = phi - pi/4 - if (values[kPhiTildePP] < 0) { - values[kPhiTildePP] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kPhiTildePP] = RecoDecay::constrainAngle(values[kPhiPP] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 } else { - values[kPhiTildePP] = values[kPhiPP] - 0.75 * TMath::Pi(); // phi_tilde = phi - 3pi/4 - if (values[kPhiTildePP] < 0) { - values[kPhiTildePP] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kPhiTildePP] = RecoDecay::constrainAngle(values[kPhiPP] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 } } else { values[kPhiTildePP] = -999; // not computable @@ -4364,10 +4274,11 @@ void VarManager::FillPairME(T1 const& t1, T2 const& t2, float* values) if (useRM) { double randomCostheta = gRandom->Uniform(-1., 1.); - double randomPhi = gRandom->Uniform(0., 2. * TMath::Pi()); + double randomPhi = gRandom->Uniform(0., o2::constants::math::TwoPI); ROOT::Math::XYZVectorF zaxis_RM(randomCostheta, std::sqrt(1 - randomCostheta * randomCostheta) * std::cos(randomPhi), std::sqrt(1 - randomCostheta * randomCostheta) * std::sin(randomPhi)); - if (fgUsedVars[kCosThetaRM]) + if (fgUsedVars[kCosThetaRM]) { values[kCosThetaRM] = zaxis_RM.Dot(v_CM); + } } } @@ -4455,7 +4366,7 @@ void VarManager::FillPairMEAcrossTFs(T const& t1, T const& t2, float* values) values[kMass] = v12.M(); values[kPt] = v12.Pt(); values[kEta] = v12.Eta(); - values[kPhi] = v12.Phi() > 0 ? v12.Phi() : v12.Phi() + 2. * M_PI; + values[kPhi] = RecoDecay::constrainAngle(v12.Phi()); values[kRap] = -v12.Rapidity(); if (fgUsedVars[kCosThetaStarRandom] || fgUsedVars[kCosThetaStarFT0C]) { @@ -4543,26 +4454,18 @@ void VarManager::FillPairMC(T1 const& t1, T2 const& t2, float* values) ROOT::Math::XYZVectorF zaxis_HE{(v12.Vect()).Unit()}; ROOT::Math::XYZVectorF yaxis_HE{(Beam1_CM.Cross(Beam2_CM)).Unit()}; ROOT::Math::XYZVectorF xaxis_HE{(yaxis_HE.Cross(zaxis_HE)).Unit()}; - if (fgUsedVars[kMCCosThetaHE]) + if (fgUsedVars[kMCCosThetaHE]) { values[kMCCosThetaHE] = zaxis_HE.Dot(v_CM); + } if (fgUsedVars[kMCPhiHE]) { - values[kMCPhiHE] = TMath::ATan2(yaxis_HE.Dot(v_CM), xaxis_HE.Dot(v_CM)); - if (values[kMCPhiHE] < 0) { - values[kMCPhiHE] += 2 * TMath::Pi(); // ensure phi is in [0, 2pi] - } + values[kMCPhiHE] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_HE.Dot(v_CM), xaxis_HE.Dot(v_CM))); } if (fgUsedVars[kMCPhiTildeHE]) { if (fgUsedVars[kMCCosThetaHE] && fgUsedVars[kMCPhiHE]) { if (values[kMCCosThetaHE] > 0) { - values[kMCPhiTildeHE] = values[kMCPhiHE] - 0.25 * TMath::Pi(); // phi_tilde = phi - pi/4 - if (values[kMCPhiTildeHE] < 0) { - values[kMCPhiTildeHE] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kMCPhiTildeHE] = RecoDecay::constrainAngle(values[kMCPhiHE] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 } else { - values[kMCPhiTildeHE] = values[kMCPhiHE] - 0.75 * TMath::Pi(); // phi_tilde = phi - 3pi/4 - if (values[kMCPhiTildeHE] < 0) { - values[kMCPhiTildeHE] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kMCPhiTildeHE] = RecoDecay::constrainAngle(values[kMCPhiHE] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 } } else { values[kMCPhiTildeHE] = -999; // not computable @@ -4574,26 +4477,18 @@ void VarManager::FillPairMC(T1 const& t1, T2 const& t2, float* values) ROOT::Math::XYZVectorF zaxis_CS{(Beam1_CM - Beam2_CM).Unit()}; ROOT::Math::XYZVectorF yaxis_CS{(Beam1_CM.Cross(Beam2_CM)).Unit()}; ROOT::Math::XYZVectorF xaxis_CS{(yaxis_CS.Cross(zaxis_CS)).Unit()}; - if (fgUsedVars[kMCCosThetaCS]) + if (fgUsedVars[kMCCosThetaCS]) { values[kMCCosThetaCS] = zaxis_CS.Dot(v_CM); + } if (fgUsedVars[kMCPhiCS]) { - values[kMCPhiCS] = TMath::ATan2(yaxis_CS.Dot(v_CM), xaxis_CS.Dot(v_CM)); - if (values[kMCPhiCS] < 0) { - values[kMCPhiCS] += 2 * TMath::Pi(); // ensure phi is in [0, 2pi] - } + values[kMCPhiCS] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_CS.Dot(v_CM), xaxis_CS.Dot(v_CM))); } if (fgUsedVars[kMCPhiTildeCS]) { if (fgUsedVars[kMCCosThetaCS] && fgUsedVars[kMCPhiCS]) { if (values[kMCCosThetaCS] > 0) { - values[kMCPhiTildeCS] = values[kMCPhiCS] - 0.25 * TMath::Pi(); // phi_tilde = phi - pi/4 - if (values[kMCPhiTildeCS] < 0) { - values[kMCPhiTildeCS] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kMCPhiTildeCS] = RecoDecay::constrainAngle(values[kMCPhiCS] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 } else { - values[kMCPhiTildeCS] = values[kMCPhiCS] - 0.75 * TMath::Pi(); // phi_tilde = phi - 3pi/4 - if (values[kMCPhiTildeCS] < 0) { - values[kMCPhiTildeCS] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kMCPhiTildeCS] = RecoDecay::constrainAngle(values[kMCPhiCS] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 } } else { values[kMCPhiTildeCS] = -999; // not computable @@ -4609,23 +4504,14 @@ void VarManager::FillPairMC(T1 const& t1, T2 const& t2, float* values) values[kMCCosThetaPP] = zaxis_PP.Dot(v_CM); } if (fgUsedVars[kMCPhiPP]) { - values[kMCPhiPP] = TMath::ATan2(yaxis_PP.Dot(v_CM), xaxis_PP.Dot(v_CM)); - if (values[kMCPhiPP] < 0) { - values[kMCPhiPP] += 2 * TMath::Pi(); // ensure phi is in [0, 2pi] - } + values[kMCPhiPP] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_PP.Dot(v_CM), xaxis_PP.Dot(v_CM))); } if (fgUsedVars[kMCPhiTildePP]) { if (fgUsedVars[kMCCosThetaPP] && fgUsedVars[kMCPhiPP]) { if (values[kMCCosThetaPP] > 0) { - values[kMCPhiTildePP] = values[kMCPhiPP] - 0.25 * TMath::Pi(); // phi_tilde = phi - pi/4 - if (values[kMCPhiTildePP] < 0) { - values[kMCPhiTildePP] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kMCPhiTildePP] = RecoDecay::constrainAngle(values[kMCPhiPP] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 } else { - values[kMCPhiTildePP] = values[kMCPhiPP] - 0.75 * TMath::Pi(); // phi_tilde = phi - 3pi/4 - if (values[kMCPhiTildePP] < 0) { - values[kMCPhiTildePP] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kMCPhiTildePP] = RecoDecay::constrainAngle(values[kMCPhiPP] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 } } else { values[kMCPhiTildePP] = -999; // not computable @@ -4635,10 +4521,11 @@ void VarManager::FillPairMC(T1 const& t1, T2 const& t2, float* values) if (useRM) { double randomCostheta = gRandom->Uniform(-1., 1.); - double randomPhi = gRandom->Uniform(0., 2. * TMath::Pi()); + double randomPhi = gRandom->Uniform(0., o2::constants::math::TwoPI); ROOT::Math::XYZVectorF zaxis_RM(randomCostheta, std::sqrt(1 - randomCostheta * randomCostheta) * std::cos(randomPhi), std::sqrt(1 - randomCostheta * randomCostheta) * std::sin(randomPhi)); - if (fgUsedVars[kMCCosThetaRM]) + if (fgUsedVars[kMCCosThetaRM]) { values[kMCCosThetaRM] = zaxis_RM.Dot(v_CM); + } } } @@ -4665,14 +4552,14 @@ void VarManager::FillPairMC(T1 const& t1, T2 const& t2, float* values) ROOT::Math::PtEtaPhiMVector v_daughter = boostv12(t1.pdgCode() > 0 ? v1 : v2); // reaction plane - float phi = v_daughter.Phi() > TMath::Pi() ? v_daughter.Phi() - 2. * TMath::Pi() : v_daughter.Phi(); + float phi = RecoDecay::constrainAngle(v_daughter.Phi(), -o2::constants::math::PI); values[kDeltaPhiRP_Random] = phi - values[kRandomPsi2]; - values[kDeltaPhiRP_Random] = values[kDeltaPhiRP_Random] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiRP_Random] : values[kDeltaPhiRP_Random]; + values[kDeltaPhiRP_Random] = values[kDeltaPhiRP_Random] > o2::constants::math::PI ? o2::constants::math::TwoPI - values[kDeltaPhiRP_Random] : values[kDeltaPhiRP_Random]; values[kDeltaPhiRP_MC] = phi - values[kMCEventPlaneAngle]; - values[kDeltaPhiRP_MC] = values[kDeltaPhiRP_MC] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiRP_MC] : values[kDeltaPhiRP_MC]; + values[kDeltaPhiRP_MC] = values[kDeltaPhiRP_MC] > o2::constants::math::PI ? o2::constants::math::TwoPI - values[kDeltaPhiRP_MC] : values[kDeltaPhiRP_MC]; // fold delta phi into [-pi/2, pi/2] - values[kDeltaPhiRP_Random] = values[kDeltaPhiRP_Random] > TMath::Pi() / 2. ? TMath::Pi() - values[kDeltaPhiRP_Random] : values[kDeltaPhiRP_Random]; - values[kDeltaPhiRP_MC] = values[kDeltaPhiRP_MC] > TMath::Pi() / 2. ? TMath::Pi() - values[kDeltaPhiRP_MC] : values[kDeltaPhiRP_MC]; + values[kDeltaPhiRP_Random] = values[kDeltaPhiRP_Random] > o2::constants::math::PIHalf ? o2::constants::math::PI - values[kDeltaPhiRP_Random] : values[kDeltaPhiRP_Random]; + values[kDeltaPhiRP_MC] = values[kDeltaPhiRP_MC] > o2::constants::math::PIHalf ? o2::constants::math::PI - values[kDeltaPhiRP_MC] : values[kDeltaPhiRP_MC]; values[kCos2DeltaPhiRP_Random] = TMath::Cos(2. * (phi - values[kRandomPsi2])); values[kCos2DeltaPhiRP_MC] = TMath::Cos(2. * (phi - values[kMCEventPlaneAngle])); } @@ -4788,7 +4675,7 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2, ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), m2); ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; - values[kUsedKF] = fgUsedKF; + values[kUsedKF] = static_cast(fgUsedKF); if (!fgUsedKF) { int procCode = 0; @@ -4847,7 +4734,7 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2, // This modifies track momenta! o2::math_utils::Point3D vtxXYZ(collision.posX(), collision.posY(), collision.posZ()); std::array vtxCov{collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()}; - o2::dataformats::VertexBase primaryVertex = {std::move(vtxXYZ), std::move(vtxCov)}; + o2::dataformats::VertexBase primaryVertex = {vtxXYZ, vtxCov}; // auto primaryVertex = getPrimaryVertex(collision); auto covMatrixPV = primaryVertex.getCov(); @@ -4862,8 +4749,9 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2, v1 = {trackParVar0.getPt(), trackParVar0.getEta(), trackParVar0.getPhi(), m1}; v2 = {trackParVar1.getPt(), trackParVar1.getEta(), trackParVar1.getPhi(), m2}; v12 = v1 + v2; - if (fgPVrecalKF) + if (fgPVrecalKF) { primaryVertexNew = RecalculatePrimaryVertex(t1, t2, collision); + } } else if constexpr (pairType == kDecayToMuMu && muonHasCov) { // Get pca candidate from forward DCA fitter @@ -4928,8 +4816,9 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2, values[kVertexingTauxyProjected] = values[kVertexingLxyProjected] * v12.M() / (v12.Pt()); values[kVertexingTauxyProjectedPoleJPsiMass] = values[kVertexingLxyProjected] * o2::constants::physics::MassJPsi / (v12.Pt()); values[kVertexingTauxyProjectedNs] = values[kVertexingTauxyProjected] / o2::constants::physics::LightSpeedCm2NS; - if (fgPVrecalKF) + if (fgPVrecalKF) { values[kVertexingTauxyProjectedPoleJPsiMassRecalculatePV] = values[kVertexingLxyProjectedRecalculatePV] * o2::constants::physics::MassJPsi / (v12.Pt()); + } values[kVertexingTauzProjected] = values[kVertexingLzProjected] * v12.M() / TMath::Abs(v12.Pz()); values[kVertexingTauxyzProjected] = values[kVertexingLxyzProjected] * v12.M() / (v12.P()); } @@ -4969,10 +4858,11 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2, } if (fgUsedVars[kKFMass]) { float mass = 0., massErr = 0.; - if (!KFGeoTwoProng.GetMass(mass, massErr)) + if (!KFGeoTwoProng.GetMass(mass, massErr)) { values[kKFMass] = mass; - else + } else { values[kKFMass] = -999.; + } } if constexpr (eventHasVtxCov) { @@ -4989,14 +4879,17 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2, values[kVertexingLxyErr] = (KFPV.GetCovariance(0) + KFGeoTwoProng.GetCovariance(0)) * dxPair2PV * dxPair2PV + (KFPV.GetCovariance(2) + KFGeoTwoProng.GetCovariance(2)) * dyPair2PV * dyPair2PV + 2 * ((KFPV.GetCovariance(1) + KFGeoTwoProng.GetCovariance(1)) * dxPair2PV * dyPair2PV); values[kVertexingLzErr] = (KFPV.GetCovariance(5) + KFGeoTwoProng.GetCovariance(5)) * dzPair2PV * dzPair2PV; values[kVertexingLxyzErr] = (KFPV.GetCovariance(0) + KFGeoTwoProng.GetCovariance(0)) * dxPair2PV * dxPair2PV + (KFPV.GetCovariance(2) + KFGeoTwoProng.GetCovariance(2)) * dyPair2PV * dyPair2PV + (KFPV.GetCovariance(5) + KFGeoTwoProng.GetCovariance(5)) * dzPair2PV * dzPair2PV + 2 * ((KFPV.GetCovariance(1) + KFGeoTwoProng.GetCovariance(1)) * dxPair2PV * dyPair2PV + (KFPV.GetCovariance(3) + KFGeoTwoProng.GetCovariance(3)) * dxPair2PV * dzPair2PV + (KFPV.GetCovariance(4) + KFGeoTwoProng.GetCovariance(4)) * dyPair2PV * dzPair2PV); - if (fabs(values[kVertexingLxy]) < 1.e-8f) + if (std::fabs(values[kVertexingLxy]) < 1.e-8f) { values[kVertexingLxy] = 1.e-8f; + } values[kVertexingLxyErr] = values[kVertexingLxyErr] < 0. ? 1.e8f : std::sqrt(values[kVertexingLxyErr]) / values[kVertexingLxy]; - if (fabs(values[kVertexingLz]) < 1.e-8f) + if (std::fabs(values[kVertexingLz]) < 1.e-8f) { values[kVertexingLz] = 1.e-8f; + } values[kVertexingLzErr] = values[kVertexingLzErr] < 0. ? 1.e8f : std::sqrt(values[kVertexingLzErr]) / values[kVertexingLz]; - if (fabs(values[kVertexingLxyz]) < 1.e-8f) + if (std::fabs(values[kVertexingLxyz]) < 1.e-8f) { values[kVertexingLxyz] = 1.e-8f; + } values[kVertexingLxyzErr] = values[kVertexingLxyzErr] < 0. ? 1.e8f : std::sqrt(values[kVertexingLxyzErr]) / values[kVertexingLxyz]; values[kVertexingTauxy] = KFGeoTwoProng.GetPseudoProperDecayTime(KFPV, KFGeoTwoProng.GetMass()) / (o2::constants::physics::LightSpeedCm2NS); values[kVertexingTauz] = -1 * dzPair2PV * KFGeoTwoProng.GetMass() / (TMath::Abs(KFGeoTwoProng.GetPz()) * o2::constants::physics::LightSpeedCm2NS); @@ -5028,10 +4921,12 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2, values[kVertexingLxyzOverErr] = values[kVertexingLxyz] / values[kVertexingLxyzErr]; } - if (fgUsedVars[kKFChi2OverNDFGeo]) + if (fgUsedVars[kKFChi2OverNDFGeo]) { values[kKFChi2OverNDFGeo] = KFGeoTwoProng.GetChi2() / KFGeoTwoProng.GetNDF(); - if (fgUsedVars[kKFCosPA]) + } + if (fgUsedVars[kKFCosPA]) { values[kKFCosPA] = calculateCosPA(KFGeoTwoProng, KFPV); + } // in principle, they should be in FillTrack if (fgUsedVars[kKFTrack0DCAxyz] || fgUsedVars[kKFTrack1DCAxyz]) { @@ -5042,10 +4937,12 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2, values[kKFTrack0DCAxy] = trk0KF.GetDistanceFromVertexXY(KFPV); values[kKFTrack1DCAxy] = trk1KF.GetDistanceFromVertexXY(KFPV); } - if (fgUsedVars[kKFDCAxyzBetweenProngs]) + if (fgUsedVars[kKFDCAxyzBetweenProngs]) { values[kKFDCAxyzBetweenProngs] = trk0KF.GetDistanceFromParticle(trk1KF); - if (fgUsedVars[kKFDCAxyBetweenProngs]) + } + if (fgUsedVars[kKFDCAxyBetweenProngs]) { values[kKFDCAxyBetweenProngs] = trk0KF.GetDistanceFromParticleXY(trk1KF); + } if (fgUsedVars[kKFTracksDCAxyzMax]) { values[kKFTracksDCAxyzMax] = values[kKFTrack0DCAxyz] > values[kKFTrack1DCAxyz] ? values[kKFTrack0DCAxyz] : values[kKFTrack1DCAxyz]; @@ -5076,10 +4973,11 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2, KFGeoTopTwoProngBarrel.SetProductionVertex(KFPV); values[kKFChi2OverNDFGeoTop] = KFGeoTopTwoProngBarrel.GetChi2() / KFGeoTopTwoProngBarrel.GetNDF(); float mass = 0., massErr = 0.; - if (!KFGeoTopTwoProngBarrel.GetMass(mass, massErr)) + if (!KFGeoTopTwoProngBarrel.GetMass(mass, massErr)) { values[kKFMassGeoTop] = mass; - else + } else { values[kKFMassGeoTop] = -999.; + } } if (propToSV) { if constexpr ((pairType == kDecayToMuMu) && muonHasCov) { @@ -5090,9 +4988,9 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2, auto geoMan2 = o2::base::GeometryManager::meanMaterialBudget(t2.x(), t2.y(), t2.z(), KFGeoTwoProng.GetX(), KFGeoTwoProng.GetY(), KFGeoTwoProng.GetZ()); auto x2x01 = static_cast(geoMan1.meanX2X0); auto x2x02 = static_cast(geoMan2.meanX2X0); - float B[3]; - float xyz[3] = {0, 0, 0}; - KFGeoTwoProng.GetFieldValue(xyz, B); + std::array B{}; + std::array xyz = {0.f, 0.f, 0.f}; + KFGeoTwoProng.GetFieldValue(xyz.data(), B.data()); // TODO: find better soluton to handle cases where KF outputs negative variances /*float covXX = 0.1; float covYY = 0.1; @@ -5135,7 +5033,7 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2, values[kPt] = v12.Pt(); values[kEta] = v12.Eta(); // values[kPhi] = v12.Phi(); - values[kPhi] = v12.Phi() > 0 ? v12.Phi() : v12.Phi() + 2. * M_PI; + values[kPhi] = RecoDecay::constrainAngle(v12.Phi()); } else { values[kPt1] = t1.pt(); values[kEta1] = t1.eta(); @@ -5148,7 +5046,7 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2, } template -void VarManager::FillPairVertexingRecomputePV(C const& /*collision*/, T const& t1, T const& t2, o2::dataformats::VertexBase pvRefitted, float* values) +void VarManager::FillPairVertexingRecomputePV(C const& /*collision*/, T const& t1, T const& t2, const o2::dataformats::VertexBase& pvRefitted, float* values) { // recompute decay lenght variables using updated primary vertex @@ -5175,8 +5073,9 @@ void VarManager::FillPairVertexingRecomputePV(C const& /*collision*/, T const& t ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), m2); ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; - if (fgFitterTwoProngBarrel.getNCandidates() == 0) + if (fgFitterTwoProngBarrel.getNCandidates() == 0) { return; + } Vec3D secondaryVertex; if (!fgUsedKF) { // to be updated when seconday vertex is computed with KF @@ -5221,13 +5120,10 @@ void VarManager::FillTripletVertexing(C const& collision, T const& t1, T const& values = fgValues; } - float m1, m2, m3; + float m1 = o2::constants::physics::MassKaonCharged; + float m2 = o2::constants::physics::MassPionCharged; + float m3 = o2::constants::physics::MassPionCharged; - if (tripletType == kTripleCandidateToKPiPi) { - m1 = o2::constants::physics::MassKaonCharged; - m2 = o2::constants::physics::MassPionCharged; - m3 = o2::constants::physics::MassPionCharged; - } if (tripletType == kTripleCandidateToPKPi) { m1 = o2::constants::physics::MassProton; m2 = o2::constants::physics::MassKaonCharged; @@ -5238,7 +5134,7 @@ void VarManager::FillTripletVertexing(C const& collision, T const& t1, T const& ROOT::Math::PtEtaPhiMVector v3(t3.pt(), t3.eta(), t3.phi(), m3); ROOT::Math::PtEtaPhiMVector v123 = v1 + v2 + v3; - values[kUsedKF] = fgUsedKF; + values[kUsedKF] = static_cast(fgUsedKF); if (!fgUsedKF) { int procCode = 0; @@ -5298,7 +5194,7 @@ void VarManager::FillTripletVertexing(C const& collision, T const& t1, T const& o2::math_utils::Point3D vtxXYZ(collision.posX(), collision.posY(), collision.posZ()); std::array vtxCov{collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()}; - o2::dataformats::VertexBase primaryVertex = {std::move(vtxXYZ), std::move(vtxCov)}; + o2::dataformats::VertexBase primaryVertex = {vtxXYZ, vtxCov}; auto covMatrixPV = primaryVertex.getCov(); if (fgUsedVars[kVertexingChi2PCA]) { @@ -5378,10 +5274,11 @@ void VarManager::FillTripletVertexing(C const& collision, T const& t1, T const& } if (fgUsedVars[kKFMass]) { float mass = 0., massErr = 0.; - if (!KFGeoThreeProng.GetMass(mass, massErr)) + if (!KFGeoThreeProng.GetMass(mass, massErr)) { values[kKFMass] = mass; - else + } else { values[kKFMass] = -999.; + } } if constexpr (eventHasVtxCov) { @@ -5399,14 +5296,17 @@ void VarManager::FillTripletVertexing(C const& collision, T const& t1, T const& values[kVertexingLxyErr] = (KFPV.GetCovariance(0) + KFGeoThreeProng.GetCovariance(0)) * dxTriplet2PV * dxTriplet2PV + (KFPV.GetCovariance(2) + KFGeoThreeProng.GetCovariance(2)) * dyTriplet2PV * dyTriplet2PV + 2 * ((KFPV.GetCovariance(1) + KFGeoThreeProng.GetCovariance(1)) * dxTriplet2PV * dyTriplet2PV); values[kVertexingLzErr] = (KFPV.GetCovariance(5) + KFGeoThreeProng.GetCovariance(5)) * dzTriplet2PV * dzTriplet2PV; values[kVertexingLxyzErr] = (KFPV.GetCovariance(0) + KFGeoThreeProng.GetCovariance(0)) * dxTriplet2PV * dxTriplet2PV + (KFPV.GetCovariance(2) + KFGeoThreeProng.GetCovariance(2)) * dyTriplet2PV * dyTriplet2PV + (KFPV.GetCovariance(5) + KFGeoThreeProng.GetCovariance(5)) * dzTriplet2PV * dzTriplet2PV + 2 * ((KFPV.GetCovariance(1) + KFGeoThreeProng.GetCovariance(1)) * dxTriplet2PV * dyTriplet2PV + (KFPV.GetCovariance(3) + KFGeoThreeProng.GetCovariance(3)) * dxTriplet2PV * dzTriplet2PV + (KFPV.GetCovariance(4) + KFGeoThreeProng.GetCovariance(4)) * dyTriplet2PV * dzTriplet2PV); - if (fabs(values[kVertexingLxy]) < 1.e-8f) + if (std::fabs(values[kVertexingLxy]) < 1.e-8f) { values[kVertexingLxy] = 1.e-8f; + } values[kVertexingLxyErr] = values[kVertexingLxyErr] < 0. ? 1.e8f : std::sqrt(values[kVertexingLxyErr]) / values[kVertexingLxy]; - if (fabs(values[kVertexingLz]) < 1.e-8f) + if (std::fabs(values[kVertexingLz]) < 1.e-8f) { values[kVertexingLz] = 1.e-8f; + } values[kVertexingLzErr] = values[kVertexingLzErr] < 0. ? 1.e8f : std::sqrt(values[kVertexingLzErr]) / values[kVertexingLz]; - if (fabs(values[kVertexingLxyz]) < 1.e-8f) + if (std::fabs(values[kVertexingLxyz]) < 1.e-8f) { values[kVertexingLxyz] = 1.e-8f; + } values[kVertexingLxyzErr] = values[kVertexingLxyzErr] < 0. ? 1.e8f : std::sqrt(values[kVertexingLxyzErr]) / values[kVertexingLxyz]; values[kVertexingTauxy] = KFGeoThreeProng.GetPseudoProperDecayTime(KFPV, KFGeoThreeProng.GetMass()) / (o2::constants::physics::LightSpeedCm2NS); @@ -5443,18 +5343,16 @@ void VarManager::FillDileptonTrackVertexing(C const& collision, T1 const& lepton values = fgValues; } - float mtrack; - float mlepton1, mlepton2; + float mtrack = o2::constants::physics::MassMuon; + float mlepton1 = o2::constants::physics::MassMuon; + float mlepton2 = o2::constants::physics::MassMuon; int procCode = 0; int procCodeJpsi = 0; - values[kUsedKF] = fgUsedKF; + values[kUsedKF] = static_cast(fgUsedKF); if (!fgUsedKF) { if constexpr ((candidateType == kBcToThreeMuons) && muonHasCov) { - mlepton1 = o2::constants::physics::MassMuon; - mlepton2 = o2::constants::physics::MassMuon; - mtrack = o2::constants::physics::MassMuon; o2::track::TrackParCovFwd pars1 = FwdToTrackPar(lepton1, lepton1); o2::track::TrackParCovFwd pars2 = FwdToTrackPar(lepton2, lepton2); @@ -5463,11 +5361,11 @@ void VarManager::FillDileptonTrackVertexing(C const& collision, T1 const& lepton procCode = VarManager::fgFitterThreeProngFwd.process(pars1, pars2, pars3); procCodeJpsi = VarManager::fgFitterTwoProngFwd.process(pars1, pars2); } else if constexpr ((candidateType == kBtoJpsiEEK || candidateType == kDstarToD0KPiPi) && trackHasCov) { - if constexpr ((candidateType == kBtoJpsiEEK) && trackHasCov) { + if constexpr (candidateType == kBtoJpsiEEK) { mlepton1 = o2::constants::physics::MassElectron; mlepton2 = o2::constants::physics::MassElectron; mtrack = o2::constants::physics::MassKaonCharged; - } else if constexpr ((candidateType == kDstarToD0KPiPi) && trackHasCov) { + } else if constexpr (candidateType == kDstarToD0KPiPi) { mlepton1 = o2::constants::physics::MassKaonCharged; mlepton2 = o2::constants::physics::MassPionCharged; mtrack = o2::constants::physics::MassPionCharged; @@ -5537,13 +5435,13 @@ void VarManager::FillDileptonTrackVertexing(C const& collision, T1 const& lepton Vec3D secondaryVertex; if constexpr (eventHasVtxCov) { - std::array covMatrixPCA; + std::array covMatrixPCA{}; o2::dataformats::DCA impactParameter0; o2::dataformats::DCA impactParameter1; o2::math_utils::Point3D vtxXYZ(collision.posX(), collision.posY(), collision.posZ()); std::array vtxCov{collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()}; - o2::dataformats::VertexBase primaryVertex = {std::move(vtxXYZ), std::move(vtxCov)}; + o2::dataformats::VertexBase primaryVertex = {vtxXYZ, vtxCov}; auto covMatrixPV = primaryVertex.getCov(); if constexpr ((candidateType == kBtoJpsiEEK || candidateType == kDstarToD0KPiPi) && trackHasCov) { @@ -5596,7 +5494,7 @@ void VarManager::FillDileptonTrackVertexing(C const& collision, T1 const& lepton (v123.P() * values[VarManager::kVertexingLxyz]); } // run 2 definitions: Lxy projected onto the momentum vector of the candidate - if (fgUsedVars[kVertexingLxyProjected] || fgUsedVars[kVertexingLxyzProjected] || values[kVertexingTauxyProjected]) { + if (fgUsedVars[kVertexingLxyProjected] || fgUsedVars[kVertexingLxyzProjected] || fgUsedVars[kVertexingTauxyProjected]) { values[kVertexingLzProjected] = (secondaryVertex[2] - collision.posZ()) * v123.Pz(); values[kVertexingLzProjected] = values[kVertexingLzProjected] / TMath::Sqrt(v123.Pz() * v123.Pz()); values[kVertexingLxyProjected] = ((secondaryVertex[0] - collision.posX()) * v123.Px()) + ((secondaryVertex[1] - collision.posY()) * v123.Py()); @@ -5632,15 +5530,17 @@ void VarManager::FillDileptonTrackVertexing(C const& collision, T1 const& lepton } // Quantities between 3rd prong and candidate - if (fgUsedVars[kKFDCAxyzBetweenProngs]) + if (fgUsedVars[kKFDCAxyzBetweenProngs]) { values[kKFDCAxyzBetweenProngs] = KFGeoTwoLeptons.GetDistanceFromParticle(hadronKF); + } KFGeoThreeProng.SetConstructMethod(2); KFGeoThreeProng.AddDaughter(KFGeoTwoLeptons); KFGeoThreeProng.AddDaughter(hadronKF); - if (fgUsedVars[kKFMass]) + if (fgUsedVars[kKFMass]) { values[kKFMass] = KFGeoThreeProng.GetMass(); + } if constexpr (eventHasVtxCov) { KFPVertex kfpVertex = createKFPVertexFromCollision(collision); @@ -5656,26 +5556,32 @@ void VarManager::FillDileptonTrackVertexing(C const& collision, T1 const& lepton values[kVertexingLxyErr] = (KFPV.GetCovariance(0) + KFGeoThreeProng.GetCovariance(0)) * dxTriplet3PV * dxTriplet3PV + (KFPV.GetCovariance(2) + KFGeoThreeProng.GetCovariance(2)) * dyTriplet3PV * dyTriplet3PV + 2 * ((KFPV.GetCovariance(1) + KFGeoThreeProng.GetCovariance(1)) * dxTriplet3PV * dyTriplet3PV); values[kVertexingLzErr] = (KFPV.GetCovariance(5) + KFGeoThreeProng.GetCovariance(5)) * dzTriplet3PV * dzTriplet3PV; values[kVertexingLxyzErr] = (KFPV.GetCovariance(0) + KFGeoThreeProng.GetCovariance(0)) * dxTriplet3PV * dxTriplet3PV + (KFPV.GetCovariance(2) + KFGeoThreeProng.GetCovariance(2)) * dyTriplet3PV * dyTriplet3PV + (KFPV.GetCovariance(5) + KFGeoThreeProng.GetCovariance(5)) * dzTriplet3PV * dzTriplet3PV + 2 * ((KFPV.GetCovariance(1) + KFGeoThreeProng.GetCovariance(1)) * dxTriplet3PV * dyTriplet3PV + (KFPV.GetCovariance(3) + KFGeoThreeProng.GetCovariance(3)) * dxTriplet3PV * dzTriplet3PV + (KFPV.GetCovariance(4) + KFGeoThreeProng.GetCovariance(4)) * dyTriplet3PV * dzTriplet3PV); - if (fabs(values[kVertexingLxy]) < 1.e-8f) + if (std::fabs(values[kVertexingLxy]) < 1.e-8f) { values[kVertexingLxy] = 1.e-8f; + } values[kVertexingLxyErr] = values[kVertexingLxyErr] < 0. ? 1.e8f : std::sqrt(values[kVertexingLxyErr]) / values[kVertexingLxy]; - if (fabs(values[kVertexingLz]) < 1.e-8f) + if (std::fabs(values[kVertexingLz]) < 1.e-8f) { values[kVertexingLz] = 1.e-8f; + } values[kVertexingLzErr] = values[kVertexingLzErr] < 0. ? 1.e8f : std::sqrt(values[kVertexingLzErr]) / values[kVertexingLz]; - if (fabs(values[kVertexingLxyz]) < 1.e-8f) + if (std::fabs(values[kVertexingLxyz]) < 1.e-8f) { values[kVertexingLxyz] = 1.e-8f; + } values[kVertexingLxyzErr] = values[kVertexingLxyzErr] < 0. ? 1.e8f : std::sqrt(values[kVertexingLxyzErr]) / values[kVertexingLxyz]; - if (fgUsedVars[kVertexingTauxy]) + if (fgUsedVars[kVertexingTauxy]) { values[kVertexingTauxy] = KFGeoThreeProng.GetPseudoProperDecayTime(KFPV, KFGeoThreeProng.GetMass()) / (o2::constants::physics::LightSpeedCm2NS); - if (fgUsedVars[kVertexingTauxyErr]) + } + if (fgUsedVars[kVertexingTauxyErr]) { values[kVertexingTauxyErr] = values[kVertexingLxyErr] * KFGeoThreeProng.GetMass() / (KFGeoThreeProng.GetPt() * o2::constants::physics::LightSpeedCm2NS); + } - if (fgUsedVars[kCosPointingAngle]) + if (fgUsedVars[kCosPointingAngle]) { values[VarManager::kCosPointingAngle] = (dxTriplet3PV * KFGeoThreeProng.GetPx() + dyTriplet3PV * KFGeoThreeProng.GetPy() + dzTriplet3PV * KFGeoThreeProng.GetPz()) / (KFGeoThreeProng.GetP() * values[VarManager::kVertexingLxyz]); + } } // end calculate vertex variables // As defined in Run 2 (projected onto momentum) @@ -5759,7 +5665,7 @@ void VarManager::FillQVectorFromGFW(C const& /*collision*/, A const& compA11, A values[kM11REF] = S21A - S12A; values[kM1111REF] = S41A - 6. * S12A * S21A + 8. * S13A * S11A + 3. * S22A - 6. * S14A; values[kCORR2REF] = (norm(compA21) - S12A) / values[kM11REF]; - values[kCORR4REF] = (pow(norm(compA21), 2) + norm(compA42) - 2. * (compA42 * conj(compA21) * conj(compA21)).real() + 8. * (compA23 * conj(compA21)).real() - 4. * S12A * norm(compA21) - 6. * S14A + 2. * S22A) / values[kM1111REF]; + values[kCORR4REF] = (std::pow(norm(compA21), 2) + norm(compA42) - 2. * (compA42 * conj(compA21) * conj(compA21)).real() + 8. * (compA23 * conj(compA21)).real() - 4. * S12A * norm(compA21) - 6. * S14A + 2. * S22A) / values[kM1111REF]; values[kCORR2REF] = std::isnan(values[kM11REF]) || std::isinf(values[kM11REF]) || std::isnan(values[kCORR2REF]) || std::isinf(values[kCORR2REF]) ? 0 : values[kCORR2REF]; values[kM11REF] = std::isnan(values[kM11REF]) || std::isinf(values[kM11REF]) || std::isnan(values[kCORR2REF]) || std::isinf(values[kCORR2REF]) ? 0 : values[kM11REF]; values[kCORR4REF] = std::isnan(values[kM1111REF]) || std::isinf(values[kM1111REF]) || std::isnan(values[kCORR4REF]) || std::isinf(values[kCORR4REF]) ? 0 : values[kCORR4REF]; @@ -5935,9 +5841,9 @@ void VarManager::FillSpectatorPlane(C const& collision, float* values) values[kTimeZPA] = collision.timeZPA(); values[kTimeZPC] = collision.timeZPC(); - constexpr float beamEne = 5.36 * 0.5; - constexpr float x[4] = {-1.75, 1.75, -1.75, 1.75}; - constexpr float y[4] = {-1.75, -1.75, 1.75, 1.75}; + // constexpr float beamEne = 5.36 * 0.5; + constexpr std::array x = {-1.75f, 1.75f, -1.75f, 1.75f}; + constexpr std::array y = {-1.75f, -1.75f, 1.75f, 1.75f}; // constexpr float intcalibZNA[4] = {0.7997028, 0.8453715, 0.7879917, 0.7695486}; // constexpr float intcalibZNC[4] = {0.7631577, 0.8408003, 0.7083920, 0.7731769}; // constexpr float alpha = 0.395; // WARNING: Run 2 coorection, to be checked @@ -5968,9 +5874,9 @@ void VarManager::FillSpectatorPlane(C const& collision, float* values) } if (denZNC != 0.) { - float nSpecnC = zncCommon / beamEne; // WARNING: Run 2 coorection, to be checked - float cZNC = 1.89358 - 0.71262 / (nSpecnC + 0.71789); // WARNING: Run 2 coorection, to be checked - cZNC = 1.; + // float nSpecnC = zncCommon / beamEne; // WARNING: Run 2 correction, to be checked + // float cZNC = 1.89358 - 0.71262 / (nSpecnC + 0.71789); // WARNING: Run 2 correction, to be checked + float cZNC = 1.; values[kQ1ZNCX] = cZNC * numXZNC / denZNC; values[kQ1ZNCY] = cZNC * numYZNC / denZNC; } else { @@ -5978,9 +5884,9 @@ void VarManager::FillSpectatorPlane(C const& collision, float* values) } if (denZNA != 0.) { - float nSpecnA = znaCommon / beamEne; // WARNING: Run 2 coorection, to be checked - float cZNA = 1.89358 - 0.71262 / (nSpecnA + 0.71789); // WARNING: Run 2 coorection, to be checked - cZNA = 1.; + // float nSpecnA = znaCommon / beamEne; // WARNING: Run 2 correction, to be checked + // float cZNA = 1.89358 - 0.71262 / (nSpecnA + 0.71789); // WARNING: Run 2 correction, to be checked + float cZNA = 1.; values[kQ1ZNAX] = cZNA * numXZNA / denZNA; values[kQ1ZNAY] = cZNA * numYZNA / denZNA; } else { @@ -5996,7 +5902,7 @@ void VarManager::FillSpectatorPlane(C const& collision, float* values) values[kQ1ZNACXX] = values[kQ1ZNACYY] = values[kQ1ZNACYX] = values[kQ1ZNACXY] = 999.; } - if (znaCommon != 0 && sumZNA != 0 && zncCommon != 0 && sumZNC) { + if (znaCommon != 0 && sumZNA != 0 && zncCommon != 0 && sumZNC != 0) { values[KIntercalibZNA] = znaCommon - sumZNA; values[KIntercalibZNC] = zncCommon - sumZNC; } @@ -6069,23 +5975,23 @@ void VarManager::FillPairVn(T1 const& t1, T2 const& t2, float* values) values[kV2EP] = std::isnan(V2EP) || std::isinf(V2EP) ? 0. : V2EP; values[kWV2EP] = std::isnan(V2EP) || std::isinf(V2EP) ? 0. : 1.0; - if (std::isnan(VarManager::fgValues[VarManager::kU2Q2]) == true) { + if (std::isnan(VarManager::fgValues[VarManager::kU2Q2])) { values[kU2Q2] = -999.; values[kR2SP_AB] = -999.; values[kR2SP_AC] = -999.; values[kR2SP_BC] = -999.; } - if (std::isnan(VarManager::fgValues[VarManager::kU3Q3]) == true) { + if (std::isnan(VarManager::fgValues[VarManager::kU3Q3])) { values[kU3Q3] = -999.; values[kR3SP] = -999.; } - if (std::isnan(VarManager::fgValues[VarManager::kCos2DeltaPhi]) == true) { + if (std::isnan(VarManager::fgValues[VarManager::kCos2DeltaPhi])) { values[kCos2DeltaPhi] = -999.; values[kR2EP_AB] = -999.; values[kR2EP_AC] = -999.; values[kR2EP_BC] = -999.; } - if (std::isnan(VarManager::fgValues[VarManager::kCos3DeltaPhi]) == true) { + if (std::isnan(VarManager::fgValues[VarManager::kCos3DeltaPhi])) { values[kCos3DeltaPhi] = -999.; values[kR3EP] = -999.; } @@ -6147,8 +6053,8 @@ void VarManager::FillPairVn(T1 const& t1, T2 const& t2, float* values) Q2X0A = nNorm > 0 ? Q2X0A / nNorm : NAN; Q2Y0A = nNorm > 0 ? Q2Y0A / nNorm : NAN; - float Psi2A = getEventPlane(2, Q2X0A, Q2Y0A); - values[kPsi2A] = Psi2A; + float Psi2ATPC = getEventPlane(2, Q2X0A, Q2Y0A); + values[kPsi2A] = Psi2ATPC; // pT ~ 0.2, non-relativistic // ROOT::Math::PtEtaPhiMVector v_daughter = t1.sign() > 0 ? v1 - v2 : v2 - v1; @@ -6162,7 +6068,7 @@ void VarManager::FillPairVn(T1 const& t1, T2 const& t2, float* values) ROOT::Math::XYZVectorF yAxis_RF = zAxis_RF.Cross(zAxis).Unit(); ROOT::Math::XYZVectorF xAxis_RF = yAxis_RF.Cross(zAxis_RF).Unit(); ROOT::Math::XYZVectorF daughterVec_RF{(v_daughter.Vect()).Unit()}; - ROOT::Math::XYZVectorF b_TPC_RF = ROOT::Math::XYZVectorF(std::cos(Psi2A), std::sin(Psi2A), 0.f); + ROOT::Math::XYZVectorF b_TPC_RF = ROOT::Math::XYZVectorF(std::cos(Psi2ATPC), std::sin(Psi2ATPC), 0.f); ROOT::Math::XYZVectorF b_FT0A_RF = ROOT::Math::XYZVectorF(std::cos(Psi2B), std::sin(Psi2B), 0.f); ROOT::Math::XYZVectorF b_FT0C_RF = ROOT::Math::XYZVectorF(std::cos(Psi2C), std::sin(Psi2C), 0.f); float cosPhi = yAxis_RF.Dot(zAxis_RF.Cross(daughterVec_RF)); @@ -6178,16 +6084,16 @@ void VarManager::FillPairVn(T1 const& t1, T2 const& t2, float* values) float sinPsi2CPP = b_FT0C_RF.Dot(yAxis_RF.Cross(daughterVec_RF)); float Psi2CPP = sinPsi2CPP > 0 ? TMath::ACos(cosPsi2CPP) : -1. * TMath::ACos(cosPsi2CPP); values[kDeltaPhiPP_TPC] = phi_PP > Psi2APP ? phi_PP - Psi2APP : Psi2APP - phi_PP; - values[kDeltaPhiPP_TPC] = values[kDeltaPhiPP_TPC] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiPP_TPC] : values[kDeltaPhiPP_TPC]; + values[kDeltaPhiPP_TPC] = values[kDeltaPhiPP_TPC] > o2::constants::math::PI ? o2::constants::math::TwoPI - values[kDeltaPhiPP_TPC] : values[kDeltaPhiPP_TPC]; values[kDeltaPhiPP_FT0A] = phi_PP > Psi2BPP ? phi_PP - Psi2BPP : Psi2BPP - phi_PP; - values[kDeltaPhiPP_FT0A] = values[kDeltaPhiPP_FT0A] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiPP_FT0A] : values[kDeltaPhiPP_FT0A]; + values[kDeltaPhiPP_FT0A] = values[kDeltaPhiPP_FT0A] > o2::constants::math::PI ? o2::constants::math::TwoPI - values[kDeltaPhiPP_FT0A] : values[kDeltaPhiPP_FT0A]; values[kDeltaPhiPP_FT0C] = phi_PP > Psi2CPP ? phi_PP - Psi2CPP : Psi2CPP - phi_PP; - values[kDeltaPhiPP_FT0C] = values[kDeltaPhiPP_FT0C] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiPP_FT0C] : values[kDeltaPhiPP_FT0C]; + values[kDeltaPhiPP_FT0C] = values[kDeltaPhiPP_FT0C] > o2::constants::math::PI ? o2::constants::math::TwoPI - values[kDeltaPhiPP_FT0C] : values[kDeltaPhiPP_FT0C]; // fold delta phi to [0, pi/2] - values[kDeltaPhiPP_TPC] = values[kDeltaPhiPP_TPC] > TMath::Pi() / 2. ? TMath::Pi() - values[kDeltaPhiPP_TPC] : values[kDeltaPhiPP_TPC]; - values[kDeltaPhiPP_FT0A] = values[kDeltaPhiPP_FT0A] > TMath::Pi() / 2. ? TMath::Pi() - values[kDeltaPhiPP_FT0A] : values[kDeltaPhiPP_FT0A]; - values[kDeltaPhiPP_FT0C] = values[kDeltaPhiPP_FT0C] > TMath::Pi() / 2. ? TMath::Pi() - values[kDeltaPhiPP_FT0C] : values[kDeltaPhiPP_FT0C]; - values[kCos2DeltaPhiPP_TPC] = TMath::Cos(2. * (phi_PP - Psi2A)); + values[kDeltaPhiPP_TPC] = values[kDeltaPhiPP_TPC] > o2::constants::math::PIHalf ? o2::constants::math::PI - values[kDeltaPhiPP_TPC] : values[kDeltaPhiPP_TPC]; + values[kDeltaPhiPP_FT0A] = values[kDeltaPhiPP_FT0A] > o2::constants::math::PIHalf ? o2::constants::math::PI - values[kDeltaPhiPP_FT0A] : values[kDeltaPhiPP_FT0A]; + values[kDeltaPhiPP_FT0C] = values[kDeltaPhiPP_FT0C] > o2::constants::math::PIHalf ? o2::constants::math::PI - values[kDeltaPhiPP_FT0C] : values[kDeltaPhiPP_FT0C]; + values[kCos2DeltaPhiPP_TPC] = TMath::Cos(2. * (phi_PP - Psi2ATPC)); values[kCos2DeltaPhiPP_FT0A] = TMath::Cos(2. * (phi_PP - Psi2B)); values[kCos2DeltaPhiPP_FT0C] = TMath::Cos(2. * (phi_PP - Psi2C)); @@ -6199,18 +6105,18 @@ void VarManager::FillPairVn(T1 const& t1, T2 const& t2, float* values) values[kA2EP_PP_FT0C] = std::isnan(A2PP_FT0C) || std::isinf(A2PP_FT0C) ? -999. : A2PP_FT0C; // reaction plane - float phi = v_daughter.Phi() > TMath::Pi() ? 2. * TMath::Pi() - v_daughter.Phi() : v_daughter.Phi(); - values[kDeltaPhiRP_TPC] = phi > Psi2A ? phi - Psi2A : Psi2A - phi; - values[kDeltaPhiRP_TPC] = values[kDeltaPhiRP_TPC] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiRP_TPC] : values[kDeltaPhiRP_TPC]; + float phi = v_daughter.Phi() > o2::constants::math::PI ? o2::constants::math::TwoPI - v_daughter.Phi() : v_daughter.Phi(); + values[kDeltaPhiRP_TPC] = phi > Psi2ATPC ? phi - Psi2ATPC : Psi2ATPC - phi; + values[kDeltaPhiRP_TPC] = values[kDeltaPhiRP_TPC] > o2::constants::math::PI ? o2::constants::math::TwoPI - values[kDeltaPhiRP_TPC] : values[kDeltaPhiRP_TPC]; values[kDeltaPhiRP_FT0A] = phi > Psi2B ? phi - Psi2B : Psi2B - phi; - values[kDeltaPhiRP_FT0A] = values[kDeltaPhiRP_FT0A] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiRP_FT0A] : values[kDeltaPhiRP_FT0A]; + values[kDeltaPhiRP_FT0A] = values[kDeltaPhiRP_FT0A] > o2::constants::math::PI ? o2::constants::math::TwoPI - values[kDeltaPhiRP_FT0A] : values[kDeltaPhiRP_FT0A]; values[kDeltaPhiRP_FT0C] = phi > Psi2C ? phi - Psi2C : Psi2C - phi; - values[kDeltaPhiRP_FT0C] = values[kDeltaPhiRP_FT0C] > TMath::Pi() ? 2. * TMath::Pi() - values[kDeltaPhiRP_FT0C] : values[kDeltaPhiRP_FT0C]; + values[kDeltaPhiRP_FT0C] = values[kDeltaPhiRP_FT0C] > o2::constants::math::PI ? o2::constants::math::TwoPI - values[kDeltaPhiRP_FT0C] : values[kDeltaPhiRP_FT0C]; // fold delta phi to [0, pi/2] - values[kDeltaPhiRP_TPC] = values[kDeltaPhiRP_TPC] > TMath::Pi() / 2. ? TMath::Pi() - values[kDeltaPhiRP_TPC] : values[kDeltaPhiRP_TPC]; - values[kDeltaPhiRP_FT0A] = values[kDeltaPhiRP_FT0A] > TMath::Pi() / 2. ? TMath::Pi() - values[kDeltaPhiRP_FT0A] : values[kDeltaPhiRP_FT0A]; - values[kDeltaPhiRP_FT0C] = values[kDeltaPhiRP_FT0C] > TMath::Pi() / 2. ? TMath::Pi() - values[kDeltaPhiRP_FT0C] : values[kDeltaPhiRP_FT0C]; - values[kCos2DeltaPhiRP_TPC] = TMath::Cos(2. * (phi - Psi2A)); + values[kDeltaPhiRP_TPC] = values[kDeltaPhiRP_TPC] > o2::constants::math::PIHalf ? o2::constants::math::PI - values[kDeltaPhiRP_TPC] : values[kDeltaPhiRP_TPC]; + values[kDeltaPhiRP_FT0A] = values[kDeltaPhiRP_FT0A] > o2::constants::math::PIHalf ? o2::constants::math::PI - values[kDeltaPhiRP_FT0A] : values[kDeltaPhiRP_FT0A]; + values[kDeltaPhiRP_FT0C] = values[kDeltaPhiRP_FT0C] > o2::constants::math::PIHalf ? o2::constants::math::PI - values[kDeltaPhiRP_FT0C] : values[kDeltaPhiRP_FT0C]; + values[kCos2DeltaPhiRP_TPC] = TMath::Cos(2. * (phi - Psi2ATPC)); values[kCos2DeltaPhiRP_FT0A] = TMath::Cos(2. * (phi - Psi2B)); values[kCos2DeltaPhiRP_FT0C] = TMath::Cos(2. * (phi - Psi2C)); @@ -6331,23 +6237,16 @@ void VarManager::FillDileptonHadron(T1 const& dilepton, T2 const& hadron, float* ROOT::Math::PtEtaPhiMVector v12_Qvect = v1 - v2; double Pinv = v12.M(); double Q1 = (dilepton.mass() * dilepton.mass() - hadronMass * hadronMass) / Pinv; - values[kDileptonHadronKstar] = sqrt(Q1 * Q1 - v12_Qvect.M2()) / 2.0; + values[kDileptonHadronKstar] = std::sqrt(Q1 * Q1 - v12_Qvect.M2()) / 2.0; } if (fgUsedVars[kDeltaPhi]) { - double delta = dilepton.phi() - hadron.phi(); - if (delta > 3.0 / 2.0 * M_PI) { - delta -= 2.0 * M_PI; - } - if (delta < -0.5 * M_PI) { - delta += 2.0 * M_PI; - } - values[kDeltaPhi] = delta; + values[kDeltaPhi] = RecoDecay::constrainAngle(dilepton.phi() - hadron.phi(), -o2::constants::math::PIHalf); } if (fgUsedVars[kDeltaPhiSym]) { double delta = std::abs(dilepton.phi() - hadron.phi()); - if (delta > M_PI) { - delta = 2 * M_PI - delta; + if (delta > o2::constants::math::PI) { + delta = o2::constants::math::TwoPI - delta; } values[kDeltaPhiSym] = delta; } @@ -6404,7 +6303,7 @@ void VarManager::FillEnergyCorrelatorTriple(T1 const& lepton1, T2 const& lepton2 float randomPhi_toward = -o2::constants::math::PIHalf; float randomPhi_away = -o2::constants::math::PIHalf; - if ((deltaphi > -Transhigh * TMath::Pi() && deltaphi < -Translow * TMath::Pi()) || (deltaphi > Translow * TMath::Pi() && deltaphi < Transhigh * TMath::Pi())) { + if ((deltaphi > -Transhigh * o2::constants::math::PI && deltaphi < -Translow * o2::constants::math::PI) || (deltaphi > Translow * o2::constants::math::PI && deltaphi < Transhigh * o2::constants::math::PI)) { randomPhi_trans = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf); randomPhi_toward = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf); randomPhi_away = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf); @@ -6445,13 +6344,10 @@ void VarManager::FillEnergyCorrelatorsUnfoldingTriple(T1 const& lepton1, T2 cons dileptonmass = dilepton.mass(); } - float MassHadron; + float MassHadron = o2::constants::physics::MassPionCharged; if constexpr (pairType == kJpsiHadronMass) { MassHadron = TMath::Sqrt(t1.e() * t1.e() - t1.p() * t1.p()); } - if constexpr (pairType == kJpsiPionMass) { - MassHadron = o2::constants::physics::MassPionCharged; - } ROOT::Math::PtEtaPhiMVector v1_gen(track.pt(), track.eta(), track.phi(), o2::constants::physics::MassJPsi); ROOT::Math::PtEtaPhiMVector v2_gen(t1.pt(), t1.eta(), t1.phi(), MassHadron); float E_boost_gen = LorentzTransformJpsihadroncosChi("weight_boost", v1_gen, v2_gen); @@ -6471,7 +6367,7 @@ void VarManager::FillEnergyCorrelatorsUnfoldingTriple(T1 const& lepton1, T2 cons values[kMCCosChi_randomPhi_trans_gen] = -999.9f; float deltaphi_rec = RecoDecay::constrainAngle(dilepton.phi() - hadron.phi(), -o2::constants::math::PIHalf); - if ((deltaphi_rec > -Transhigh * TMath::Pi() && deltaphi_rec < -Translow * TMath::Pi()) || (deltaphi_rec > Translow * TMath::Pi() && deltaphi_rec < Transhigh * TMath::Pi())) { + if ((deltaphi_rec > -Transhigh * o2::constants::math::PI && deltaphi_rec < -Translow * o2::constants::math::PI) || (deltaphi_rec > Translow * o2::constants::math::PI && deltaphi_rec < Transhigh * o2::constants::math::PI)) { float randomPhi_trans_rec = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf); ROOT::Math::PtEtaPhiMVector v2_randomPhi_trans_rec(v2_rec.pt(), v2_rec.eta(), randomPhi_trans_rec, o2::constants::physics::MassPionCharged); values[kMCCosChi_randomPhi_trans_rec] = LorentzTransformJpsihadroncosChi("coschi", v1_rec, v2_randomPhi_trans_rec); @@ -6568,16 +6464,16 @@ void VarManager::FillDileptonTrackTrack(T1 const& dilepton, T2 const& hadron1, T values = fgValues; } - double defaultDileptonMass = 3.096; + double defaultDileptonMass = o2::constants::physics::MassJPsi; double hadronMass1 = o2::constants::physics::MassPionCharged; double hadronMass2 = o2::constants::physics::MassPionCharged; if (candidateType == kXtoJpsiPiPi) { - defaultDileptonMass = 3.096; + defaultDileptonMass = o2::constants::physics::MassJPsi; hadronMass1 = o2::constants::physics::MassPionCharged; hadronMass2 = o2::constants::physics::MassPionCharged; } if (candidateType == kChictoJpsiEE) { - defaultDileptonMass = 3.096; + defaultDileptonMass = o2::constants::physics::MassJPsi; hadronMass1 = o2::constants::physics::MassElectron; hadronMass2 = o2::constants::physics::MassElectron; } @@ -6610,7 +6506,7 @@ void VarManager::FillDileptonTrackTrack(T1 const& dilepton, T2 const& hadron1, T values[kQ] = v123.M() - defaultDileptonMass - v23.M(); values[kDeltaR1] = ROOT::Math::VectorUtil::DeltaR(v1, v2); values[kDeltaR2] = ROOT::Math::VectorUtil::DeltaR(v1, v3); - values[kDeltaR] = sqrt(pow(values[kDeltaR1], 2) + pow(values[kDeltaR2], 2)); + values[kDeltaR] = std::sqrt(std::pow(values[kDeltaR1], 2) + std::pow(values[kDeltaR2], 2)); values[kRap] = v123.Rapidity(); } } @@ -6630,15 +6526,10 @@ void VarManager::FillDileptonTrackTrackVertexing(C const& collision, T1 const& l values = fgValues; } - float mtrack1, mtrack2; - float mlepton1, mlepton2; - - if constexpr (candidateType == kXtoJpsiPiPi || candidateType == kPsi2StoJpsiPiPi) { - mlepton1 = o2::constants::physics::MassElectron; - mlepton2 = o2::constants::physics::MassElectron; - mtrack1 = o2::constants::physics::MassPionCharged; - mtrack2 = o2::constants::physics::MassPionCharged; - } + float mlepton1 = o2::constants::physics::MassElectron; + float mlepton2 = o2::constants::physics::MassElectron; + float mtrack1 = o2::constants::physics::MassPionCharged; + float mtrack2 = o2::constants::physics::MassPionCharged; ROOT::Math::PtEtaPhiMVector v1(lepton1.pt(), lepton1.eta(), lepton1.phi(), mlepton1); ROOT::Math::PtEtaPhiMVector v2(lepton2.pt(), lepton2.eta(), lepton2.phi(), mlepton2); @@ -6649,7 +6540,7 @@ void VarManager::FillDileptonTrackTrackVertexing(C const& collision, T1 const& l int procCodeDilepton = 0; int procCodeDileptonTrackTrack = 0; - values[kUsedKF] = fgUsedKF; + values[kUsedKF] = static_cast(fgUsedKF); if (!fgUsedKF) { // create covariance matrix std::array lepton1pars = {lepton1.y(), lepton1.z(), lepton1.snp(), lepton1.tgl(), lepton1.signed1Pt()}; @@ -6699,57 +6590,57 @@ void VarManager::FillDileptonTrackTrackVertexing(C const& collision, T1 const& l values[kVertexingTauxyProjected] = -999.; values[kVertexingTauxyzProjected] = -999.; return; - } else { - Vec3D secondaryVertex; - std::array covMatrixPCA; - secondaryVertex = fgFitterFourProngBarrel.getPCACandidate(); - covMatrixPCA = fgFitterFourProngBarrel.calcPCACovMatrixFlat(); - - o2::math_utils::Point3D vtxXYZ(collision.posX(), collision.posY(), collision.posZ()); - std::array vtxCov{collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()}; - o2::dataformats::VertexBase primaryVertex = {std::move(vtxXYZ), std::move(vtxCov)}; - auto covMatrixPV = primaryVertex.getCov(); - - double phi = std::atan2(secondaryVertex[1] - collision.posY(), secondaryVertex[0] - collision.posX()); - double theta = std::atan2(secondaryVertex[2] - collision.posZ(), - std::sqrt((secondaryVertex[0] - collision.posX()) * (secondaryVertex[0] - collision.posX()) + - (secondaryVertex[1] - collision.posY()) * (secondaryVertex[1] - collision.posY()))); - - values[kVertexingLxy] = (collision.posX() - secondaryVertex[0]) * (collision.posX() - secondaryVertex[0]) + - (collision.posY() - secondaryVertex[1]) * (collision.posY() - secondaryVertex[1]); - values[kVertexingLz] = (collision.posZ() - secondaryVertex[2]) * (collision.posZ() - secondaryVertex[2]); - values[kVertexingLxyz] = values[kVertexingLxy] + values[kVertexingLz]; - values[kVertexingLxy] = std::sqrt(values[kVertexingLxy]); - values[kVertexingLz] = std::sqrt(values[kVertexingLz]); - values[kVertexingLxyz] = std::sqrt(values[kVertexingLxyz]); - - values[kVertexingLxyzErr] = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); - values[kVertexingLxyErr] = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); - values[kVertexingLzErr] = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, 0, theta) + getRotatedCovMatrixXX(covMatrixPCA, 0, theta)); + } - values[kVertexingTauz] = (collision.posZ() - secondaryVertex[2]) * v1234.M() / (TMath::Abs(v1234.Pz()) * o2::constants::physics::LightSpeedCm2NS); - values[kVertexingTauxy] = values[kVertexingLxy] * v1234.M() / (v1234.Pt() * o2::constants::physics::LightSpeedCm2NS); - - values[kVertexingTauzErr] = values[kVertexingLzErr] * v1234.M() / (TMath::Abs(v1234.Pz()) * o2::constants::physics::LightSpeedCm2NS); - values[kVertexingTauxyErr] = values[kVertexingLxyErr] * v1234.M() / (v1234.Pt() * o2::constants::physics::LightSpeedCm2NS); - - values[kCosPointingAngle] = ((secondaryVertex[0] - collision.posX()) * v1234.Px() + - (secondaryVertex[1] - collision.posY()) * v1234.Py() + - (secondaryVertex[2] - collision.posZ()) * v1234.Pz()) / - (v1234.P() * values[VarManager::kVertexingLxyz]); - // // run 2 definitions: Decay length projected onto the momentum vector of the candidate - values[kVertexingLzProjected] = (secondaryVertex[2] - collision.posZ()) * v1234.Pz(); - values[kVertexingLzProjected] = values[kVertexingLzProjected] / TMath::Sqrt(v1234.Pz() * v1234.Pz()); - values[kVertexingLxyProjected] = ((secondaryVertex[0] - collision.posX()) * v1234.Px()) + ((secondaryVertex[1] - collision.posY()) * v1234.Py()); - values[kVertexingLxyProjected] = values[kVertexingLxyProjected] / TMath::Sqrt((v1234.Px() * v1234.Px()) + (v1234.Py() * v1234.Py())); - values[kVertexingLxyzProjected] = ((secondaryVertex[0] - collision.posX()) * v1234.Px()) + ((secondaryVertex[1] - collision.posY()) * v1234.Py()) + ((secondaryVertex[2] - collision.posZ()) * v1234.Pz()); - values[kVertexingLxyzProjected] = values[kVertexingLxyzProjected] / TMath::Sqrt((v1234.Px() * v1234.Px()) + (v1234.Py() * v1234.Py()) + (v1234.Pz() * v1234.Pz())); - - values[kVertexingTauzProjected] = values[kVertexingLzProjected] * v1234.M() / TMath::Abs(v1234.Pz()); - values[kVertexingTauxyProjected] = values[kVertexingLxyProjected] * v1234.M() / (v1234.Pt()); - values[kVertexingTauxyzProjected] = values[kVertexingLxyzProjected] * v1234.M() / (v1234.P()); - } - } else if (fgUsedKF) { + Vec3D secondaryVertex; + std::array covMatrixPCA{}; + secondaryVertex = fgFitterFourProngBarrel.getPCACandidate(); + covMatrixPCA = fgFitterFourProngBarrel.calcPCACovMatrixFlat(); + + o2::math_utils::Point3D vtxXYZ(collision.posX(), collision.posY(), collision.posZ()); + std::array vtxCov{collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()}; + o2::dataformats::VertexBase primaryVertex = {vtxXYZ, vtxCov}; + auto covMatrixPV = primaryVertex.getCov(); + + double phi = std::atan2(secondaryVertex[1] - collision.posY(), secondaryVertex[0] - collision.posX()); + double theta = std::atan2(secondaryVertex[2] - collision.posZ(), + std::sqrt((secondaryVertex[0] - collision.posX()) * (secondaryVertex[0] - collision.posX()) + + (secondaryVertex[1] - collision.posY()) * (secondaryVertex[1] - collision.posY()))); + + values[kVertexingLxy] = (collision.posX() - secondaryVertex[0]) * (collision.posX() - secondaryVertex[0]) + + (collision.posY() - secondaryVertex[1]) * (collision.posY() - secondaryVertex[1]); + values[kVertexingLz] = (collision.posZ() - secondaryVertex[2]) * (collision.posZ() - secondaryVertex[2]); + values[kVertexingLxyz] = values[kVertexingLxy] + values[kVertexingLz]; + values[kVertexingLxy] = std::sqrt(values[kVertexingLxy]); + values[kVertexingLz] = std::sqrt(values[kVertexingLz]); + values[kVertexingLxyz] = std::sqrt(values[kVertexingLxyz]); + + values[kVertexingLxyzErr] = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); + values[kVertexingLxyErr] = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); + values[kVertexingLzErr] = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, 0, theta) + getRotatedCovMatrixXX(covMatrixPCA, 0, theta)); + + values[kVertexingTauz] = (collision.posZ() - secondaryVertex[2]) * v1234.M() / (TMath::Abs(v1234.Pz()) * o2::constants::physics::LightSpeedCm2NS); + values[kVertexingTauxy] = values[kVertexingLxy] * v1234.M() / (v1234.Pt() * o2::constants::physics::LightSpeedCm2NS); + + values[kVertexingTauzErr] = values[kVertexingLzErr] * v1234.M() / (TMath::Abs(v1234.Pz()) * o2::constants::physics::LightSpeedCm2NS); + values[kVertexingTauxyErr] = values[kVertexingLxyErr] * v1234.M() / (v1234.Pt() * o2::constants::physics::LightSpeedCm2NS); + + values[kCosPointingAngle] = ((secondaryVertex[0] - collision.posX()) * v1234.Px() + + (secondaryVertex[1] - collision.posY()) * v1234.Py() + + (secondaryVertex[2] - collision.posZ()) * v1234.Pz()) / + (v1234.P() * values[VarManager::kVertexingLxyz]); + // // run 2 definitions: Decay length projected onto the momentum vector of the candidate + values[kVertexingLzProjected] = (secondaryVertex[2] - collision.posZ()) * v1234.Pz(); + values[kVertexingLzProjected] = values[kVertexingLzProjected] / TMath::Sqrt(v1234.Pz() * v1234.Pz()); + values[kVertexingLxyProjected] = ((secondaryVertex[0] - collision.posX()) * v1234.Px()) + ((secondaryVertex[1] - collision.posY()) * v1234.Py()); + values[kVertexingLxyProjected] = values[kVertexingLxyProjected] / TMath::Sqrt((v1234.Px() * v1234.Px()) + (v1234.Py() * v1234.Py())); + values[kVertexingLxyzProjected] = ((secondaryVertex[0] - collision.posX()) * v1234.Px()) + ((secondaryVertex[1] - collision.posY()) * v1234.Py()) + ((secondaryVertex[2] - collision.posZ()) * v1234.Pz()); + values[kVertexingLxyzProjected] = values[kVertexingLxyzProjected] / TMath::Sqrt((v1234.Px() * v1234.Px()) + (v1234.Py() * v1234.Py()) + (v1234.Pz() * v1234.Pz())); + + values[kVertexingTauzProjected] = values[kVertexingLzProjected] * v1234.M() / TMath::Abs(v1234.Pz()); + values[kVertexingTauxyProjected] = values[kVertexingLxyProjected] * v1234.M() / (v1234.Pt()); + values[kVertexingTauxyzProjected] = values[kVertexingLxyzProjected] * v1234.M() / (v1234.P()); + } else { KFParticle lepton1KF; // lepton1 KFParticle lepton2KF; // lepton2 KFParticle KFGeoTwoLeptons; @@ -6817,10 +6708,11 @@ void VarManager::FillDileptonTrackTrackVertexing(C const& collision, T1 const& l if (fgUsedVars[kKFMass]) { float mass = 0., massErr = 0.; - if (!KFGeoFourProng.GetMass(mass, massErr)) + if (!KFGeoFourProng.GetMass(mass, massErr)) { values[kKFMass] = mass; - else + } else { values[kKFMass] = -999.; + } } KFPVertex kfpVertex = createKFPVertexFromCollision(collision); @@ -6838,14 +6730,17 @@ void VarManager::FillDileptonTrackTrackVertexing(C const& collision, T1 const& l values[kVertexingLzErr] = (KFPV.GetCovariance(5) + KFGeoFourProng.GetCovariance(5)) * dzQuadlet2PV * dzQuadlet2PV; values[kVertexingLxyzErr] = (KFPV.GetCovariance(0) + KFGeoFourProng.GetCovariance(0)) * dxQuadlet2PV * dxQuadlet2PV + (KFPV.GetCovariance(2) + KFGeoFourProng.GetCovariance(2)) * dyQuadlet2PV * dyQuadlet2PV + (KFPV.GetCovariance(5) + KFGeoFourProng.GetCovariance(5)) * dzQuadlet2PV * dzQuadlet2PV + 2 * ((KFPV.GetCovariance(1) + KFGeoFourProng.GetCovariance(1)) * dxQuadlet2PV * dyQuadlet2PV + (KFPV.GetCovariance(3) + KFGeoFourProng.GetCovariance(3)) * dxQuadlet2PV * dzQuadlet2PV + (KFPV.GetCovariance(4) + KFGeoFourProng.GetCovariance(4)) * dyQuadlet2PV * dzQuadlet2PV); - if (fabs(values[kVertexingLxy]) < 1.e-8f) + if (std::fabs(values[kVertexingLxy]) < 1.e-8f) { values[kVertexingLxy] = 1.e-8f; + } values[kVertexingLxyErr] = values[kVertexingLxyErr] < 0. ? 1.e8f : std::sqrt(values[kVertexingLxyErr]) / values[kVertexingLxy]; - if (fabs(values[kVertexingLz]) < 1.e-8f) + if (std::fabs(values[kVertexingLz]) < 1.e-8f) { values[kVertexingLz] = 1.e-8f; + } values[kVertexingLzErr] = values[kVertexingLzErr] < 0. ? 1.e8f : std::sqrt(values[kVertexingLzErr]) / values[kVertexingLz]; - if (fabs(values[kVertexingLxyz]) < 1.e-8f) + if (std::fabs(values[kVertexingLxyz]) < 1.e-8f) { values[kVertexingLxyz] = 1.e-8f; + } values[kVertexingLxyzErr] = values[kVertexingLxyzErr] < 0. ? 1.e8f : std::sqrt(values[kVertexingLxyzErr]) / values[kVertexingLxyz]; values[kVertexingTauxy] = KFGeoFourProng.GetPseudoProperDecayTime(KFPV, KFGeoFourProng.GetMass()) / (o2::constants::physics::LightSpeedCm2NS); @@ -6870,8 +6765,6 @@ void VarManager::FillDileptonTrackTrackVertexing(C const& collision, T1 const& l values[kVertexingTauxyProjectedNs] = values[kVertexingTauxyProjected] / o2::constants::physics::LightSpeedCm2NS; values[kVertexingTauzProjected] = values[kVertexingLzProjected] * KFGeoFourProng.GetMass() / TMath::Abs(KFGeoFourProng.GetPz()); values[kKFChi2OverNDFGeo] = KFGeoFourProng.GetChi2() / KFGeoFourProng.GetNDF(); - } else { - return; } } @@ -6883,11 +6776,11 @@ void VarManager::FillQuadMC(T1 const& dilepton, T2 const& track1, T2 const& trac values = fgValues; } - double defaultDileptonMass = 3.096; + double defaultDileptonMass = o2::constants::physics::MassJPsi; double hadronMass1 = o2::constants::physics::MassPionCharged; double hadronMass2 = o2::constants::physics::MassPionCharged; if (candidateType == kXtoJpsiPiPi) { - defaultDileptonMass = 3.096; + defaultDileptonMass = o2::constants::physics::MassJPsi; hadronMass1 = o2::constants::physics::MassPionCharged; hadronMass2 = o2::constants::physics::MassPionCharged; } @@ -6905,7 +6798,7 @@ void VarManager::FillQuadMC(T1 const& dilepton, T2 const& track1, T2 const& trac values[kQ] = v123.M() - defaultDileptonMass - v23.M(); values[kDeltaR1] = ROOT::Math::VectorUtil::DeltaR(v1, v2); values[kDeltaR2] = ROOT::Math::VectorUtil::DeltaR(v1, v3); - values[kDeltaR] = sqrt(pow(values[kDeltaR1], 2) + pow(values[kDeltaR2], 2)); + values[kDeltaR] = std::sqrt(std::pow(values[kDeltaR1], 2) + std::pow(values[kDeltaR2], 2)); values[kDitrackMass] = v23.M(); values[kDitrackPt] = v23.Pt(); } @@ -7046,7 +6939,7 @@ void VarManager::FillBdtScore(T1 const& bdtScore, float* values) } //__________________________________________________________________ template -float VarManager::LorentzTransformJpsihadroncosChi(TString Option, T1 const& v1, T2 const& v2) +float VarManager::LorentzTransformJpsihadroncosChi(const TString& Option, T1 const& v1, T2 const& v2) { float value = -999.0f; auto beta_v1 = v1.BoostToCM(); @@ -7226,7 +7119,7 @@ void VarManager::FillTrackAlice3(T const& track, float* values) values[kCharge] = track.sign(); if (fgUsedVars[kPVContributor]) { - values[kPVContributor] = (track.flags() & o2::aod::track::PVContributor) > 0; + values[kPVContributor] = static_cast((track.flags() & o2::aod::track::PVContributor) > 0); } if (fgUsedVars[kITSClusterMap]) { @@ -7387,7 +7280,7 @@ void VarManager::FillPairAlice3(T1 const& t1, T2 const& t2, float* values) values[kPt] = v12.Pt(); values[kEta] = v12.Eta(); // values[kPhi] = v12.Phi(); - values[kPhi] = v12.Phi() > 0 ? v12.Phi() : v12.Phi() + 2. * M_PI; + values[kPhi] = RecoDecay::constrainAngle(v12.Phi()); values[kRap] = -v12.Rapidity(); double Ptot1 = TMath::Sqrt(v1.Px() * v1.Px() + v1.Py() * v1.Py() + v1.Pz() * v1.Pz()); double Ptot2 = TMath::Sqrt(v2.Px() * v2.Px() + v2.Py() * v2.Py() + v2.Pz() * v2.Pz()); @@ -7401,14 +7294,7 @@ void VarManager::FillPairAlice3(T1 const& t1, T2 const& t2, float* values) values[kPhi2] = t2.phi(); if (fgUsedVars[kDeltaPhiPair2]) { - double phipair2 = v1.Phi() - v2.Phi(); - if (phipair2 > 3 * TMath::Pi() / 2) { - values[kDeltaPhiPair2] = phipair2 - 2 * TMath::Pi(); - } else if (phipair2 < -TMath::Pi() / 2) { - values[kDeltaPhiPair2] = phipair2 + 2 * TMath::Pi(); - } else { - values[kDeltaPhiPair2] = phipair2; - } + values[kDeltaPhiPair2] = RecoDecay::constrainAngle(v1.Phi() - v2.Phi(), -o2::constants::math::PIHalf); } if (fgUsedVars[kDeltaEtaPair2]) { @@ -7428,10 +7314,12 @@ void VarManager::FillPairAlice3(T1 const& t1, T2 const& t2, float* values) values[kOpeningAngle] = 0.; } else { double arg = scalar / Ptot12; - if (arg > 1.) + if (arg > 1.) { arg = 1.; - if (arg < -1) + } + if (arg < -1) { arg = -1; + } values[kOpeningAngle] = TMath::ACos(arg); } } @@ -7456,26 +7344,18 @@ void VarManager::FillPairAlice3(T1 const& t1, T2 const& t2, float* values) ROOT::Math::XYZVectorF zaxis_HE{(v12.Vect()).Unit()}; ROOT::Math::XYZVectorF yaxis_HE{(Beam1_CM.Cross(Beam2_CM)).Unit()}; ROOT::Math::XYZVectorF xaxis_HE{(yaxis_HE.Cross(zaxis_HE)).Unit()}; - if (fgUsedVars[kCosThetaHE]) + if (fgUsedVars[kCosThetaHE]) { values[kCosThetaHE] = zaxis_HE.Dot(v_CM); + } if (fgUsedVars[kPhiHE]) { - values[kPhiHE] = TMath::ATan2(yaxis_HE.Dot(v_CM), xaxis_HE.Dot(v_CM)); - if (values[kPhiHE] < 0) { - values[kPhiHE] += 2 * TMath::Pi(); // ensure phi is in [0, 2pi] - } + values[kPhiHE] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_HE.Dot(v_CM), xaxis_HE.Dot(v_CM))); } if (fgUsedVars[kPhiTildeHE]) { if (fgUsedVars[kCosThetaHE] && fgUsedVars[kPhiHE]) { if (values[kCosThetaHE] > 0) { - values[kPhiTildeHE] = values[kPhiHE] - 0.25 * TMath::Pi(); // phi_tilde = phi - pi/4 - if (values[kPhiTildeHE] < 0) { - values[kPhiTildeHE] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kPhiTildeHE] = RecoDecay::constrainAngle(values[kPhiHE] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 } else { - values[kPhiTildeHE] = values[kPhiHE] - 0.75 * TMath::Pi(); // phi_tilde = phi - 3pi/4 - if (values[kPhiTildeHE] < 0) { - values[kPhiTildeHE] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kPhiTildeHE] = RecoDecay::constrainAngle(values[kPhiHE] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 } } else { values[kPhiTildeHE] = -999; // not computable @@ -7487,26 +7367,18 @@ void VarManager::FillPairAlice3(T1 const& t1, T2 const& t2, float* values) ROOT::Math::XYZVectorF zaxis_CS{(Beam1_CM - Beam2_CM).Unit()}; ROOT::Math::XYZVectorF yaxis_CS{(Beam1_CM.Cross(Beam2_CM)).Unit()}; ROOT::Math::XYZVectorF xaxis_CS{(yaxis_CS.Cross(zaxis_CS)).Unit()}; - if (fgUsedVars[kCosThetaCS]) + if (fgUsedVars[kCosThetaCS]) { values[kCosThetaCS] = zaxis_CS.Dot(v_CM); + } if (fgUsedVars[kPhiCS]) { - values[kPhiCS] = TMath::ATan2(yaxis_CS.Dot(v_CM), xaxis_CS.Dot(v_CM)); - if (values[kPhiCS] < 0) { - values[kPhiCS] += 2 * TMath::Pi(); // ensure phi is in [0, 2pi] - } + values[kPhiCS] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_CS.Dot(v_CM), xaxis_CS.Dot(v_CM))); } if (fgUsedVars[kPhiTildeCS]) { if (fgUsedVars[kCosThetaCS] && fgUsedVars[kPhiCS]) { if (values[kCosThetaCS] > 0) { - values[kPhiTildeCS] = values[kPhiCS] - 0.25 * TMath::Pi(); // phi_tilde = phi - pi/4 - if (values[kPhiTildeCS] < 0) { - values[kPhiTildeCS] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kPhiTildeCS] = RecoDecay::constrainAngle(values[kPhiCS] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 } else { - values[kPhiTildeCS] = values[kPhiCS] - 0.75 * TMath::Pi(); // phi_tilde = phi - 3pi/4 - if (values[kPhiTildeCS] < 0) { - values[kPhiTildeCS] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kPhiTildeCS] = RecoDecay::constrainAngle(values[kPhiCS] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 } } else { values[kPhiTildeCS] = -999; // not computable @@ -7522,23 +7394,14 @@ void VarManager::FillPairAlice3(T1 const& t1, T2 const& t2, float* values) values[kCosThetaPP] = zaxis_PP.Dot(v_CM) / std::sqrt(zaxis_PP.Mag2()); } if (fgUsedVars[kPhiPP]) { - values[kPhiPP] = TMath::ATan2(yaxis_PP.Dot(v_CM), xaxis_PP.Dot(v_CM)); - if (values[kPhiPP] < 0) { - values[kPhiPP] += 2 * TMath::Pi(); // ensure phi is in [0, 2pi] - } + values[kPhiPP] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_PP.Dot(v_CM), xaxis_PP.Dot(v_CM))); } if (fgUsedVars[kPhiTildePP]) { if (fgUsedVars[kCosThetaPP] && fgUsedVars[kPhiPP]) { if (values[kCosThetaPP] > 0) { - values[kPhiTildePP] = values[kPhiPP] - 0.25 * TMath::Pi(); // phi_tilde = phi - pi/4 - if (values[kPhiTildePP] < 0) { - values[kPhiTildePP] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kPhiTildePP] = RecoDecay::constrainAngle(values[kPhiPP] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 } else { - values[kPhiTildePP] = values[kPhiPP] - 0.75 * TMath::Pi(); // phi_tilde = phi - 3pi/4 - if (values[kPhiTildePP] < 0) { - values[kPhiTildePP] += 2 * TMath::Pi(); // ensure phi_tilde is in [0, 2pi] - } + values[kPhiTildePP] = RecoDecay::constrainAngle(values[kPhiPP] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 } } else { values[kPhiTildePP] = -999; // not computable @@ -7548,10 +7411,11 @@ void VarManager::FillPairAlice3(T1 const& t1, T2 const& t2, float* values) if (useRM) { double randomCostheta = gRandom->Uniform(-1., 1.); - double randomPhi = gRandom->Uniform(0., 2. * TMath::Pi()); + double randomPhi = gRandom->Uniform(0., o2::constants::math::TwoPI); ROOT::Math::XYZVectorF zaxis_RM(randomCostheta, std::sqrt(1 - randomCostheta * randomCostheta) * std::cos(randomPhi), std::sqrt(1 - randomCostheta * randomCostheta) * std::sin(randomPhi)); - if (fgUsedVars[kCosThetaRM]) + if (fgUsedVars[kCosThetaRM]) { values[kCosThetaRM] = zaxis_RM.Dot(v_CM); + } } }