From f866218fac1326d48e595fc3357156711d4cecdd Mon Sep 17 00:00:00 2001 From: feisenhu Date: Fri, 26 Apr 2024 14:53:11 +0200 Subject: [PATCH] [PWGEM/DQ] seperate calculation of phiV, add pairType as template for FillPairMC --- PWGDQ/Core/VarManager.h | 221 ++++++++++++++++----------- PWGDQ/Tasks/dqEfficiency.cxx | 14 +- PWGEM/Dilepton/Tasks/MCtemplates.cxx | 10 +- 3 files changed, 142 insertions(+), 103 deletions(-) diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index fae39cd2663..cae1cf45886 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -823,8 +823,8 @@ class VarManager : public TObject static void FillTriple(T1 const& t1, T2 const& t2, T3 const& t3, float* values = nullptr, PairCandidateType pairType = kTripleCandidateToEEPhoton); template static void FillPairME(T1 const& t1, T2 const& t2, float* values = nullptr); - template - static void FillPairMC(T1 const& t1, T2 const& t2, float* values = nullptr, PairCandidateType pairType = kDecayToEE); + template + static void FillPairMC(T1 const& t1, T2 const& t2, float* values = nullptr); template static void FillTripleMC(T1 const& t1, T2 const& t2, T3 const& t3, float* values = nullptr, PairCandidateType pairType = kTripleCandidateToEEPhoton); template @@ -912,6 +912,8 @@ class VarManager : public TObject template static KFPVertex createKFPVertexFromCollision(const T& collision); static float calculateCosPA(KFParticle kfp, KFParticle PV); + template + static float calculatePhiV(const T1& t1, const T2& t2); static o2::vertexing::DCAFitterN<2> fgFitterTwoProngBarrel; static o2::vertexing::DCAFitterN<3> fgFitterThreeProngBarrel; @@ -2234,93 +2236,7 @@ void VarManager::FillPair(T1 const& t1, T2 const& t2, float* values) } } if (fgUsedVars[kPairPhiv]) { - // cos(phiv) = w*a /|w||a| - // with w = u x v - // and a = u x z / |u x z| , unit vector perpendicular to v12 and z-direction (magnetic field) - // u = v12 / |v12| , the unit vector of v12 - // v = v1 x v2 / |v1 x v2| , unit vector perpendicular to v1 and v2 - - float bz = fgFitterTwoProngBarrel.getBz(); - - bool swapTracks = false; - if (v1.Pt() < v2.Pt()) { // ordering of track, pt1 > pt2 - ROOT::Math::PtEtaPhiMVector v3 = v1; - v1 = v2; - v2 = v3; - swapTracks = true; - } - - // momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by definition. - // vector product of pep X pem - float vpx = 0, vpy = 0, vpz = 0; - if (t1.sign() * t2.sign() > 0) { // Like Sign - if (!swapTracks) { - if (bz * t1.sign() < 0) { - vpx = v1.Py() * v2.Pz() - v1.Pz() * v2.Py(); - vpy = v1.Pz() * v2.Px() - v1.Px() * v2.Pz(); - vpz = v1.Px() * v2.Py() - v1.Py() * v2.Px(); - } else { - vpx = v2.Py() * v1.Pz() - v2.Pz() * v1.Py(); - vpy = v2.Pz() * v1.Px() - v2.Px() * v1.Pz(); - vpz = v2.Px() * v1.Py() - v2.Py() * v1.Px(); - } - } else { // swaped tracks - if (bz * t2.sign() < 0) { - vpx = v1.Py() * v2.Pz() - v1.Pz() * v2.Py(); - vpy = v1.Pz() * v2.Px() - v1.Px() * v2.Pz(); - vpz = v1.Px() * v2.Py() - v1.Py() * v2.Px(); - } else { - vpx = v2.Py() * v1.Pz() - v2.Pz() * v1.Py(); - vpy = v2.Pz() * v1.Px() - v2.Px() * v1.Pz(); - vpz = v2.Px() * v1.Py() - v2.Py() * v1.Px(); - } - } - } else { // Unlike Sign - if (!swapTracks) { - if (bz * t1.sign() > 0) { - vpx = v1.Py() * v2.Pz() - v1.Pz() * v2.Py(); - vpy = v1.Pz() * v2.Px() - v1.Px() * v2.Pz(); - vpz = v1.Px() * v2.Py() - v1.Py() * v2.Px(); - } else { - vpx = v2.Py() * v1.Pz() - v2.Pz() * v1.Py(); - vpy = v2.Pz() * v1.Px() - v2.Px() * v1.Pz(); - vpz = v2.Px() * v1.Py() - v2.Py() * v1.Px(); - } - } else { // swaped tracks - if (bz * t2.sign() > 0) { - vpx = v1.Py() * v2.Pz() - v1.Pz() * v2.Py(); - vpy = v1.Pz() * v2.Px() - v1.Px() * v2.Pz(); - vpz = v1.Px() * v2.Py() - v1.Py() * v2.Px(); - } else { - vpx = v2.Py() * v1.Pz() - v2.Pz() * v1.Py(); - vpy = v2.Pz() * v1.Px() - v2.Px() * v1.Pz(); - vpz = v2.Px() * v1.Py() - v2.Py() * v1.Px(); - } - } - } - - // unit vector of pep X pem - float vx = vpx / TMath::Sqrt(vpx * vpx + vpy * vpy + vpz * vpz); - float vy = vpy / TMath::Sqrt(vpx * vpx + vpy * vpy + vpz * vpz); - float vz = vpz / TMath::Sqrt(vpx * vpx + vpy * vpy + vpz * vpz); - - float px = v12.Px(); - float py = v12.Py(); - float pz = v12.Pz(); - - // unit vector of (pep+pem) - float ux = px / TMath::Sqrt(px * px + py * py + pz * pz); - float uy = py / TMath::Sqrt(px * px + py * py + pz * pz); - float uz = pz / TMath::Sqrt(px * px + py * py + pz * pz); - float ax = uy / TMath::Sqrt(ux * ux + uy * uy); - float ay = -ux / TMath::Sqrt(ux * ux + uy * uy); - - // The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz) - float wx = uy * vz - uz * vy; - float wy = uz * vx - ux * vz; - // by construction, (wx,wy,wz) must be a unit vector. Measure angle between (wx,wy,wz) and (ax,ay,0). - // The angle between them should be small if the pair is conversion. This function then returns values close to pi! - values[kPairPhiv] = TMath::ACos(wx * ax + wy * ay); // phiv in [0,pi] //cosPhiV = wx * ax + wy * ay; + values[kPairPhiv] = calculatePhiV(t1, t2); } } @@ -2402,10 +2318,14 @@ void VarManager::FillPairME(T1 const& t1, T2 const& t2, float* values) values[kU2Q2Ev2] = values[kQ2X0A2] * std::cos(2 * v2.Phi()) + values[kQ2Y0A2] * std::sin(2 * v2.Phi()); values[kCos2DeltaPhiMu1] = std::cos(2 * (v1.Phi() - v12.Phi())); values[kCos2DeltaPhiMu2] = std::cos(2 * (v2.Phi() - v12.Phi())); + + if (fgUsedVars[kPairPhiv]) { + values[kPairPhiv] = calculatePhiV(t1, t2); + } } -template -void VarManager::FillPairMC(T1 const& t1, T2 const& t2, float* values, PairCandidateType pairType) +template +void VarManager::FillPairMC(T1 const& t1, T2 const& t2, float* values) { if (!values) { values = fgValues; @@ -2436,6 +2356,10 @@ void VarManager::FillPairMC(T1 const& t1, T2 const& t2, float* values, PairCandi values[kEta] = v12.Eta(); values[kPhi] = v12.Phi(); values[kRap] = -v12.Rapidity(); + + if (fgUsedVars[kPairPhiv]) { + values[kPairPhiv] = calculatePhiV(t1, t2); + } } template @@ -3523,4 +3447,119 @@ void VarManager::FillDileptonTrackTrack(T1 const& dilepton, T2 const& hadron1, T } } +//__________________________________________________________________ +template +float VarManager::calculatePhiV(T1 const& t1, T2 const& t2) +{ + // cos(phiv) = w*a /|w||a| + // with w = u x v + // and a = u x z / |u x z| , unit vector perpendicular to v12 and z-direction (magnetic field) + // u = v12 / |v12| , the unit vector of v12 + // v = v1 x v2 / |v1 x v2| , unit vector perpendicular to v1 and v2 + + float m1 = o2::constants::physics::MassElectron; + float m2 = o2::constants::physics::MassElectron; + if constexpr (pairType == kDecayToMuMu) { + m1 = o2::constants::physics::MassMuon; + m2 = o2::constants::physics::MassMuon; + } + + if constexpr (pairType == kDecayToPiPi) { + m1 = o2::constants::physics::MassPionCharged; + m2 = o2::constants::physics::MassPionCharged; + } + + if constexpr (pairType == kElectronMuon) { + m2 = o2::constants::physics::MassMuon; + } + + ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), m1); + ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), m2); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + + float pairPhiV = -999; + float bz = fgFitterTwoProngBarrel.getBz(); + + bool swapTracks = false; + if (v1.Pt() < v2.Pt()) { // ordering of track, pt1 > pt2 + ROOT::Math::PtEtaPhiMVector v3 = v1; + v1 = v2; + v2 = v3; + swapTracks = true; + } + + // momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by definition. + // vector product of pep X pem + float vpx = 0, vpy = 0, vpz = 0; + if (t1.sign() * t2.sign() > 0) { // Like Sign + if (!swapTracks) { + if (bz * t1.sign() < 0) { + vpx = v1.Py() * v2.Pz() - v1.Pz() * v2.Py(); + vpy = v1.Pz() * v2.Px() - v1.Px() * v2.Pz(); + vpz = v1.Px() * v2.Py() - v1.Py() * v2.Px(); + } else { + vpx = v2.Py() * v1.Pz() - v2.Pz() * v1.Py(); + vpy = v2.Pz() * v1.Px() - v2.Px() * v1.Pz(); + vpz = v2.Px() * v1.Py() - v2.Py() * v1.Px(); + } + } else { // swaped tracks + if (bz * t2.sign() < 0) { + vpx = v1.Py() * v2.Pz() - v1.Pz() * v2.Py(); + vpy = v1.Pz() * v2.Px() - v1.Px() * v2.Pz(); + vpz = v1.Px() * v2.Py() - v1.Py() * v2.Px(); + } else { + vpx = v2.Py() * v1.Pz() - v2.Pz() * v1.Py(); + vpy = v2.Pz() * v1.Px() - v2.Px() * v1.Pz(); + vpz = v2.Px() * v1.Py() - v2.Py() * v1.Px(); + } + } + } else { // Unlike Sign + if (!swapTracks) { + if (bz * t1.sign() > 0) { + vpx = v1.Py() * v2.Pz() - v1.Pz() * v2.Py(); + vpy = v1.Pz() * v2.Px() - v1.Px() * v2.Pz(); + vpz = v1.Px() * v2.Py() - v1.Py() * v2.Px(); + } else { + vpx = v2.Py() * v1.Pz() - v2.Pz() * v1.Py(); + vpy = v2.Pz() * v1.Px() - v2.Px() * v1.Pz(); + vpz = v2.Px() * v1.Py() - v2.Py() * v1.Px(); + } + } else { // swaped tracks + if (bz * t2.sign() > 0) { + vpx = v1.Py() * v2.Pz() - v1.Pz() * v2.Py(); + vpy = v1.Pz() * v2.Px() - v1.Px() * v2.Pz(); + vpz = v1.Px() * v2.Py() - v1.Py() * v2.Px(); + } else { + vpx = v2.Py() * v1.Pz() - v2.Pz() * v1.Py(); + vpy = v2.Pz() * v1.Px() - v2.Px() * v1.Pz(); + vpz = v2.Px() * v1.Py() - v2.Py() * v1.Px(); + } + } + } + + // unit vector of pep X pem + float vx = vpx / TMath::Sqrt(vpx * vpx + vpy * vpy + vpz * vpz); + float vy = vpy / TMath::Sqrt(vpx * vpx + vpy * vpy + vpz * vpz); + float vz = vpz / TMath::Sqrt(vpx * vpx + vpy * vpy + vpz * vpz); + + float px = v12.Px(); + float py = v12.Py(); + float pz = v12.Pz(); + + // unit vector of (pep+pem) + float ux = px / TMath::Sqrt(px * px + py * py + pz * pz); + float uy = py / TMath::Sqrt(px * px + py * py + pz * pz); + float uz = pz / TMath::Sqrt(px * px + py * py + pz * pz); + float ax = uy / TMath::Sqrt(ux * ux + uy * uy); + float ay = -ux / TMath::Sqrt(ux * ux + uy * uy); + + // The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz) + float wx = uy * vz - uz * vy; + float wy = uz * vx - ux * vz; + // by construction, (wx,wy,wz) must be a unit vector. Measure angle between (wx,wy,wz) and (ax,ay,0). + // The angle between them should be small if the pair is conversion. This function then returns values close to pi! + pairPhiV = TMath::ACos(wx * ax + wy * ay); // phiv in [0,pi] //cosPhiV = wx * ax + wy * ay; + return pairPhiV; +} + #endif // PWGDQ_CORE_VARMANAGER_H_ diff --git a/PWGDQ/Tasks/dqEfficiency.cxx b/PWGDQ/Tasks/dqEfficiency.cxx index e661f4f7563..21b60728b4f 100644 --- a/PWGDQ/Tasks/dqEfficiency.cxx +++ b/PWGDQ/Tasks/dqEfficiency.cxx @@ -836,7 +836,7 @@ struct AnalysisSameEventPairing { } // end loop over barrel track pairs } // end runPairing - template + template void runMCGen(TTracksMC& groupedMCTracks) { // loop over mc stack and fill histograms for pure MC truth signals @@ -879,7 +879,7 @@ struct AnalysisSameEventPairing { checked = sig.CheckSignal(false, t1, t2); } if (checked) { - VarManager::FillPairMC(t1, t2); + VarManager::FillPairMC(t1, t2); fHistMan->FillHistClass(Form("MCTruthGenPair_%s", sig.GetName()), VarManager::fgValues); } } @@ -901,7 +901,7 @@ struct AnalysisSameEventPairing { runPairing(event, tracks, tracks, eventsMC, tracksMC); auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); groupedMCTracks.bindInternalIndicesTo(&tracksMC); - runMCGen(groupedMCTracks); + runMCGen(groupedMCTracks); } void processDecayToEEVertexingSkimmed(soa::Filtered::iterator const& event, @@ -916,7 +916,7 @@ struct AnalysisSameEventPairing { runPairing(event, tracks, tracks, eventsMC, tracksMC); auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); groupedMCTracks.bindInternalIndicesTo(&tracksMC); - runMCGen(groupedMCTracks); + runMCGen(groupedMCTracks); } void processDecayToMuMuSkimmed(soa::Filtered::iterator const& event, @@ -931,7 +931,7 @@ struct AnalysisSameEventPairing { runPairing(event, muons, muons, eventsMC, tracksMC); auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); groupedMCTracks.bindInternalIndicesTo(&tracksMC); - runMCGen(groupedMCTracks); + runMCGen(groupedMCTracks); } void processDecayToMuMuVertexingSkimmed(soa::Filtered::iterator const& event, @@ -946,7 +946,7 @@ struct AnalysisSameEventPairing { runPairing(event, muons, muons, eventsMC, tracksMC); auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); groupedMCTracks.bindInternalIndicesTo(&tracksMC); - runMCGen(groupedMCTracks); + runMCGen(groupedMCTracks); } /*void processElectronMuonSkimmed(soa::Filtered::iterator const& event, @@ -959,7 +959,7 @@ struct AnalysisSameEventPairing { runPairing(event, tracks, muons, eventsMC, tracksMC); auto groupedMCTracks = tracksMC.sliceBy(aod::reducedtrackMC::reducedMCeventId, event.reducedMCevent().globalIndex()); - runMCGen(groupedMCTracks); + runMCGen(groupedMCTracks); }*/ void processDummy(MyEvents&) { diff --git a/PWGEM/Dilepton/Tasks/MCtemplates.cxx b/PWGEM/Dilepton/Tasks/MCtemplates.cxx index 71a7ee44c85..ef84d8859eb 100644 --- a/PWGEM/Dilepton/Tasks/MCtemplates.cxx +++ b/PWGEM/Dilepton/Tasks/MCtemplates.cxx @@ -577,7 +577,7 @@ struct AnalysisSameEventPairing { } // end loop over barrel track pairs } // end runPairing - template + template void runMCGenPair(TTracksMC const& groupedMCTracks) { // loop over mc stack and fill histograms for pure MC truth signals @@ -602,7 +602,7 @@ struct AnalysisSameEventPairing { checked = sig.CheckSignal(true, t1, t2); } if (checked) { - VarManager::FillPairMC(t1, t2); + VarManager::FillPairMC(t1, t2); fHistMan->FillHistClass(Form("MCTruthGenPair_%s", sig.GetName()), VarManager::fgValues); } } @@ -626,7 +626,7 @@ struct AnalysisSameEventPairing { runPairing(tracks, tracks); auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); groupedMCTracks.bindInternalIndicesTo(&tracksMC); - runMCGenPair(groupedMCTracks); + runMCGenPair(groupedMCTracks); } void processDecayToEESkimmedWithCov(soa::Filtered::iterator const& event, @@ -641,7 +641,7 @@ struct AnalysisSameEventPairing { runPairing(tracks, tracks); auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); groupedMCTracks.bindInternalIndicesTo(&tracksMC); - runMCGenPair(groupedMCTracks); + runMCGenPair(groupedMCTracks); } void processDecayToEEAOD(soa::Filtered::iterator const& event, @@ -656,7 +656,7 @@ struct AnalysisSameEventPairing { runPairing(tracks, tracks); auto groupedMCTracks = tracksMC.sliceBy(perMcCollision, event.mcCollision().globalIndex()); groupedMCTracks.bindInternalIndicesTo(&tracksMC); - runMCGenPair(groupedMCTracks); + runMCGenPair(groupedMCTracks); } void processDummy(MyEvents&)