Skip to content

Commit

Permalink
[PWGEM/DQ] seperate calculation of phiV, add pairType as template for…
Browse files Browse the repository at this point in the history 
… FillPairMC
  • Loading branch information
@feisenhu
feisenhu committed Apr 26, 2024
1 parent 52114c9 commit f866218
Show file tree
Hide file tree
Showing 3 changed files with 142 additions and 103 deletions.
221 changes: 130 additions & 91 deletions PWGDQ/Core/VarManager.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <int pairType, typename T1, typename T2>
static void FillPairME(T1 const& t1, T2 const& t2, float* values = nullptr);
template <typename T1, typename T2>
static void FillPairMC(T1 const& t1, T2 const& t2, float* values = nullptr, PairCandidateType pairType = kDecayToEE);
template <int pairType, typename T1, typename T2>
static void FillPairMC(T1 const& t1, T2 const& t2, float* values = nullptr);
template <typename T1, typename T2, typename T3>
static void FillTripleMC(T1 const& t1, T2 const& t2, T3 const& t3, float* values = nullptr, PairCandidateType pairType = kTripleCandidateToEEPhoton);
template <int pairType, uint32_t collFillMap, uint32_t fillMap, typename C, typename T>
Expand Down Expand Up @@ -912,6 +912,8 @@ class VarManager : public TObject
template <typename T>
static KFPVertex createKFPVertexFromCollision(const T& collision);
static float calculateCosPA(KFParticle kfp, KFParticle PV);
template <int pairType, typename T1, typename T2>
static float calculatePhiV(const T1& t1, const T2& t2);

static o2::vertexing::DCAFitterN<2> fgFitterTwoProngBarrel;
static o2::vertexing::DCAFitterN<3> fgFitterThreeProngBarrel;
Expand Down Expand Up @@ -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<pairType>(t1, t2);
}
}

Expand Down Expand Up @@ -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<pairType>(t1, t2);
}
}

template <typename T1, typename T2>
void VarManager::FillPairMC(T1 const& t1, T2 const& t2, float* values, PairCandidateType pairType)
template <int pairType, typename T1, typename T2>
void VarManager::FillPairMC(T1 const& t1, T2 const& t2, float* values)
{
if (!values) {
values = fgValues;
Expand Down Expand Up @@ -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<pairType>(t1, t2);
}
}

template <typename T1, typename T2, typename T3>
Expand Down Expand Up @@ -3523,4 +3447,119 @@ void VarManager::FillDileptonTrackTrack(T1 const& dilepton, T2 const& hadron1, T
}
}

//__________________________________________________________________
template <int pairType, typename T1, typename T2>
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_
14 changes: 7 additions & 7 deletions PWGDQ/Tasks/dqEfficiency.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -836,7 +836,7 @@ struct AnalysisSameEventPairing {
} // end loop over barrel track pairs
} // end runPairing

template <typename TTracksMC>
template <int TPairType,typename TTracksMC>
void runMCGen(TTracksMC& groupedMCTracks)
{
// loop over mc stack and fill histograms for pure MC truth signals
Expand Down Expand Up @@ -879,7 +879,7 @@ struct AnalysisSameEventPairing {
checked = sig.CheckSignal(false, t1, t2);
}
if (checked) {
VarManager::FillPairMC(t1, t2);
VarManager::FillPairMC<TPairType>(t1, t2);
fHistMan->FillHistClass(Form("MCTruthGenPair_%s", sig.GetName()), VarManager::fgValues);
}
}
Expand All @@ -901,7 +901,7 @@ struct AnalysisSameEventPairing {
runPairing<VarManager::kDecayToEE, gkEventFillMap, gkMCEventFillMap, gkTrackFillMap>(event, tracks, tracks, eventsMC, tracksMC);
auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex());
groupedMCTracks.bindInternalIndicesTo(&tracksMC);
runMCGen(groupedMCTracks);
runMCGen<VarManager::kDecayToEE>(groupedMCTracks);
}

void processDecayToEEVertexingSkimmed(soa::Filtered<MyEventsVtxCovSelected>::iterator const& event,
Expand All @@ -916,7 +916,7 @@ struct AnalysisSameEventPairing {
runPairing<VarManager::kDecayToEE, gkEventFillMapWithCov, gkMCEventFillMap, gkTrackFillMapWithCov>(event, tracks, tracks, eventsMC, tracksMC);
auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex());
groupedMCTracks.bindInternalIndicesTo(&tracksMC);
runMCGen(groupedMCTracks);
runMCGen<VarManager::kDecayToEE>(groupedMCTracks);
}

void processDecayToMuMuSkimmed(soa::Filtered<MyEventsSelected>::iterator const& event,
Expand All @@ -931,7 +931,7 @@ struct AnalysisSameEventPairing {
runPairing<VarManager::kDecayToMuMu, gkEventFillMap, gkMCEventFillMap, gkMuonFillMap>(event, muons, muons, eventsMC, tracksMC);
auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex());
groupedMCTracks.bindInternalIndicesTo(&tracksMC);
runMCGen(groupedMCTracks);
runMCGen<VarManager::kDecayToMuMu>(groupedMCTracks);
}

void processDecayToMuMuVertexingSkimmed(soa::Filtered<MyEventsVtxCovSelected>::iterator const& event,
Expand All @@ -946,7 +946,7 @@ struct AnalysisSameEventPairing {
runPairing<VarManager::kDecayToMuMu, gkEventFillMapWithCov, gkMCEventFillMap, gkMuonFillMapWithCov>(event, muons, muons, eventsMC, tracksMC);
auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex());
groupedMCTracks.bindInternalIndicesTo(&tracksMC);
runMCGen(groupedMCTracks);
runMCGen<VarManager::kDecayToMuMu>(groupedMCTracks);
}

/*void processElectronMuonSkimmed(soa::Filtered<MyEventsSelected>::iterator const& event,
Expand All @@ -959,7 +959,7 @@ struct AnalysisSameEventPairing {
runPairing<VarManager::kElectronMuon, gkEventFillMap, gkMCEventFillMap, gkTrackFillMap>(event, tracks, muons, eventsMC, tracksMC);
auto groupedMCTracks = tracksMC.sliceBy(aod::reducedtrackMC::reducedMCeventId, event.reducedMCevent().globalIndex());
runMCGen(groupedMCTracks);
runMCGen<VarManager::kElectronMuon>(groupedMCTracks);
}*/
void processDummy(MyEvents&)
{
Expand Down
10 changes: 5 additions & 5 deletions PWGEM/Dilepton/Tasks/MCtemplates.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -577,7 +577,7 @@ struct AnalysisSameEventPairing {
} // end loop over barrel track pairs
} // end runPairing

template <typename TTracksMC>
template <int TPairType, typename TTracksMC>
void runMCGenPair(TTracksMC const& groupedMCTracks)
{
// loop over mc stack and fill histograms for pure MC truth signals
Expand All @@ -602,7 +602,7 @@ struct AnalysisSameEventPairing {
checked = sig.CheckSignal(true, t1, t2);
}
if (checked) {
VarManager::FillPairMC(t1, t2);
VarManager::FillPairMC<TPairType>(t1, t2);
fHistMan->FillHistClass(Form("MCTruthGenPair_%s", sig.GetName()), VarManager::fgValues);
}
}
Expand All @@ -626,7 +626,7 @@ struct AnalysisSameEventPairing {
runPairing<VarManager::kDecayToEE, gkTrackFillMap>(tracks, tracks);
auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex());
groupedMCTracks.bindInternalIndicesTo(&tracksMC);
runMCGenPair(groupedMCTracks);
runMCGenPair<VarManager::kDecayToEE>(groupedMCTracks);
}

void processDecayToEESkimmedWithCov(soa::Filtered<MyEventsVtxCovSelected>::iterator const& event,
Expand All @@ -641,7 +641,7 @@ struct AnalysisSameEventPairing {
runPairing<VarManager::kDecayToEE, gkTrackFillMapWithCov>(tracks, tracks);
auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex());
groupedMCTracks.bindInternalIndicesTo(&tracksMC);
runMCGenPair(groupedMCTracks);
runMCGenPair<VarManager::kDecayToEE>(groupedMCTracks);
}

void processDecayToEEAOD(soa::Filtered<MyEventsSelectedAOD>::iterator const& event,
Expand All @@ -656,7 +656,7 @@ struct AnalysisSameEventPairing {
runPairing<VarManager::kDecayToEE, gkTrackFillMapAOD>(tracks, tracks);
auto groupedMCTracks = tracksMC.sliceBy(perMcCollision, event.mcCollision().globalIndex());
groupedMCTracks.bindInternalIndicesTo(&tracksMC);
runMCGenPair(groupedMCTracks);
runMCGenPair<VarManager::kDecayToEE>(groupedMCTracks);
}

void processDummy(MyEvents&)
Expand Down

0 comments on commit f866218

Please sign in to comment.