From af30ac4aa39020ce5eb876b22c23736079c4b303 Mon Sep 17 00:00:00 2001 From: dmallick2 Date: Tue, 13 Oct 2026 14:21:07 +0200 Subject: [PATCH 1/4] Add GJ axis definition and fill histograms --- PWGDQ/Core/HistogramsLibrary.cxx | 30 ++ PWGDQ/Core/VarManager.cxx | 36 +++ PWGDQ/Core/VarManager.h | 284 +++++++++++++++++- PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx | 32 ++ 4 files changed, 378 insertions(+), 4 deletions(-) diff --git a/PWGDQ/Core/HistogramsLibrary.cxx b/PWGDQ/Core/HistogramsLibrary.cxx index 662a89f514e..fda0e671aad 100644 --- a/PWGDQ/Core/HistogramsLibrary.cxx +++ b/PWGDQ/Core/HistogramsLibrary.cxx @@ -1308,6 +1308,20 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h std::array xmax = {4., 3., 100, 1., 6.28}; hm->AddHistogram(histClass, "MC_Dielectron_Mass_Pt_Cent_cosThetaHE", "", 5, varspTHE.data(), bins.data(), xmin.data(), xmax.data(), 0, -1, kFALSE); } + if (groupStr.CompareTo("polarization-dielectron-pbpb-midy-gjb1-gen") == 0) { + std::array varspTGJB1 = {VarManager::kMCMass, VarManager::kMCPt, VarManager::kMCEventCentrFT0C, VarManager::kMCCosThetaGJB1, VarManager::kMCPhiGJB1}; + std::array bins = {50, 30, 10, 20, 10}; + std::array xmin = {2., 0., 0, -1., 0.0}; + std::array xmax = {4., 3., 100, 1., 6.28}; + hm->AddHistogram(histClass, "MC_Dielectron_Mass_Pt_Cent_cosThetaGJB1", "", 5, varspTGJB1.data(), bins.data(), xmin.data(), xmax.data(), 0, -1, kFALSE); + } + if (groupStr.CompareTo("polarization-dielectron-pbpb-midy-gjb2-gen") == 0) { + std::array varspTGJB2 = {VarManager::kMCMass, VarManager::kMCPt, VarManager::kMCEventCentrFT0C, VarManager::kMCCosThetaGJB2, VarManager::kMCPhiGJB2}; + std::array bins = {50, 30, 10, 20, 10}; + std::array xmin = {2., 0., 0, -1., 0.0}; + std::array xmax = {4., 3., 100, 1., 6.28}; + hm->AddHistogram(histClass, "MC_Dielectron_Mass_Pt_Cent_cosThetaGJB2", "", 5, varspTGJB2.data(), bins.data(), xmin.data(), xmax.data(), 0, -1, kFALSE); + } if (groupStr.CompareTo("polarization-dielectron-pbpb-midy-cs-gen") == 0) { std::array varspTCS = {VarManager::kMCMass, VarManager::kMCPt, VarManager::kMCEventCentrFT0C, VarManager::kMCCosThetaCS, VarManager::kMCPhiCS}; std::array bins = {50, 30, 10, 20, 10}; @@ -1422,6 +1436,22 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h hm->AddHistogram(histClass, "Dielectron_Mass_Pt_Cent_cosThetaHE", "", 5, varsHEpbpb.data(), binspT.data(), xminpT.data(), xmaxpT.data(), nullptr, -1, kFALSE); } + if (subGroupStr.Contains("dielectron-polarization-gjb1-pbpb")) { + std::array varsGJB1pbpb = {VarManager::kMass, VarManager::kPt, VarManager::kCentFT0C, VarManager::kCosThetaGJB1, VarManager::kPhiGJB1}; + std::array binspT = {100, 30, 10, 10, 10}; + std::array xminpT = {2., 0., 0, -1., 0.0}; + std::array xmaxpT = {4.5, 3., 100, 1., 6.28}; + hm->AddHistogram(histClass, "Dielectron_Mass_Pt_Cent_cosThetaGJB1", "", 5, varsGJB1pbpb.data(), binspT.data(), xminpT.data(), xmaxpT.data(), nullptr, -1, kFALSE); + } + + if (subGroupStr.Contains("dielectron-polarization-gjb2-pbpb")) { + std::array varsGJB2pbpb = {VarManager::kMass, VarManager::kPt, VarManager::kCentFT0C, VarManager::kCosThetaGJB2, VarManager::kPhiGJB2}; + std::array binspT = {100, 30, 10, 10, 10}; + std::array xminpT = {2., 0., 0, -1., 0.0}; + std::array xmaxpT = {4.5, 3., 100, 1., 6.28}; + hm->AddHistogram(histClass, "Dielectron_Mass_Pt_Cent_cosThetaGJB2", "", 5, varsGJB2pbpb.data(), binspT.data(), xminpT.data(), xmaxpT.data(), nullptr, -1, kFALSE); + } + if (subGroupStr.Contains("dielectron-polarization-cs-pbpb")) { std::array varsCSpbpb = {VarManager::kMass, VarManager::kPt, VarManager::kCentFT0C, VarManager::kCosThetaCS, VarManager::kPhiCS}; std::array binspT = {100, 30, 10, 10, 10}; diff --git a/PWGDQ/Core/VarManager.cxx b/PWGDQ/Core/VarManager.cxx index d1b193ec12c..ce0d01b55c3 100644 --- a/PWGDQ/Core/VarManager.cxx +++ b/PWGDQ/Core/VarManager.cxx @@ -1273,6 +1273,18 @@ void VarManager::SetDefaultVarNames() fgVariableUnits[kMCPhiCS] = "rad"; fgVariableNames[kMCPhiTildeCS] = "MC #tilde{#varphi}_{CS}"; fgVariableUnits[kMCPhiTildeCS] = "rad"; + fgVariableNames[kMCCosThetaGJB1] = "MC cos#it{#theta}_{GJ1}"; + fgVariableUnits[kMCCosThetaGJB1] = ""; + fgVariableNames[kMCPhiGJB1] = "MC #varphi_{GJ1}"; + fgVariableUnits[kMCPhiGJB1] = "rad"; + fgVariableNames[kMCPhiTildeGJB1] = "MC #tilde{#varphi}_{GJ1}"; + fgVariableUnits[kMCPhiTildeGJB1] = "rad"; + fgVariableNames[kMCCosThetaGJB2] = "MC cos#it{#theta}_{GJ2}"; + fgVariableUnits[kMCCosThetaGJB2] = ""; + fgVariableNames[kMCPhiGJB2] = "MC #varphi_{GJ2}"; + fgVariableUnits[kMCPhiGJB2] = "rad"; + fgVariableNames[kMCPhiTildeGJB2] = "MC #tilde{#varphi}_{GJ2}"; + fgVariableUnits[kMCPhiTildeGJB2] = "rad"; fgVariableNames[kMCCosThetaPP] = "MC cos(#theta_{PP})"; fgVariableUnits[kMCCosThetaPP] = ""; fgVariableNames[kMCPhiPP] = "MC #varphi_{PP}"; @@ -1715,6 +1727,18 @@ void VarManager::SetDefaultVarNames() fgVariableUnits[kPhiCS] = "rad."; fgVariableNames[kPhiTildeCS] = "#tilde{#varphi}_{CS}"; fgVariableUnits[kPhiTildeCS] = "rad."; + fgVariableNames[kCosThetaGJB1] = "cos#it{#theta}_{GJ1}"; + fgVariableUnits[kCosThetaGJB1] = ""; + fgVariableNames[kPhiGJB1] = "#varphi_{GJ1}"; + fgVariableUnits[kPhiGJB1] = "rad."; + fgVariableNames[kPhiTildeGJB1] = "#tilde{#varphi}_{GJ1}"; + fgVariableUnits[kPhiTildeGJB1] = "rad."; + fgVariableNames[kCosThetaGJB2] = "cos#it{#theta}_{GJ2}"; + fgVariableUnits[kCosThetaGJB2] = ""; + fgVariableNames[kPhiGJB2] = "#varphi_{GJ2}"; + fgVariableUnits[kPhiGJB2] = "rad."; + fgVariableNames[kPhiTildeGJB2] = "#tilde{#varphi}_{GJ2}"; + fgVariableUnits[kPhiTildeGJB2] = "rad."; fgVariableNames[kCosThetaPP] = "cos#it{#theta}_{PP}"; fgVariableUnits[kCosThetaPP] = ""; fgVariableNames[kPhiPP] = "#varphi_{PP}"; @@ -2559,6 +2583,12 @@ void VarManager::SetDefaultVarNames() fgVarNamesMap["kMCCosThetaCS"] = kMCCosThetaCS; fgVarNamesMap["kMCPhiCS"] = kMCPhiCS; fgVarNamesMap["kMCPhiTildeCS"] = kMCPhiTildeCS; + fgVarNamesMap["kMCCosThetaGJB1"] = kMCCosThetaGJB1; + fgVarNamesMap["kMCPhiGJB1"] = kMCPhiGJB1; + fgVarNamesMap["kMCPhiTildeGJB1"] = kMCPhiTildeGJB1; + fgVarNamesMap["kMCCosThetaGJB2"] = kMCCosThetaGJB2; + fgVarNamesMap["kMCPhiGJB2"] = kMCPhiGJB2; + fgVarNamesMap["kMCPhiTildeGJB2"] = kMCPhiTildeGJB2; fgVarNamesMap["kMCCosThetaPP"] = kMCCosThetaPP; fgVarNamesMap["kMCPhiPP"] = kMCPhiPP; fgVarNamesMap["kMCPhiTildePP"] = kMCPhiTildePP; @@ -2605,6 +2635,12 @@ void VarManager::SetDefaultVarNames() fgVarNamesMap["kCosThetaCS"] = kCosThetaCS; fgVarNamesMap["kPhiCS"] = kPhiCS; fgVarNamesMap["kPhiTildeCS"] = kPhiTildeCS; + fgVarNamesMap["kCosThetaGJB1"] = kCosThetaGJB1; + fgVarNamesMap["kPhiGJB1"] = kPhiGJB1; + fgVarNamesMap["kPhiTildeGJB1"] = kPhiTildeGJB1; + fgVarNamesMap["kCosThetaGJB2"] = kCosThetaGJB2; + fgVarNamesMap["kPhiGJB2"] = kPhiGJB2; + fgVarNamesMap["kPhiTildeGJB2"] = kPhiTildeGJB2; fgVarNamesMap["kCosThetaPP"] = kCosThetaPP; fgVarNamesMap["kPhiPP"] = kPhiPP; fgVarNamesMap["kPhiTildePP"] = kPhiTildePP; diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index ed09b758655..8e41ab60b9b 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -776,6 +776,12 @@ class VarManager : public TObject kMCCosThetaCS, kMCPhiCS, kMCPhiTildeCS, + kMCCosThetaGJB1, + kMCPhiGJB1, + kMCPhiTildeGJB1, + kMCCosThetaGJB2, + kMCPhiGJB2, + kMCPhiTildeGJB2, kMCCosThetaPP, kMCPhiPP, kMCPhiTildePP, @@ -831,6 +837,12 @@ class VarManager : public TObject kCosThetaCS, kPhiCS, kPhiTildeCS, + kCosThetaGJB1, + kPhiGJB1, + kPhiTildeGJB1, + kCosThetaGJB2, + kPhiGJB2, + kPhiTildeGJB2, kCosThetaPP, kPhiPP, kPhiTildePP, @@ -3889,10 +3901,12 @@ void VarManager::FillPair(T1 const& t1, T2 const& t2, float* values) // polarization parameters bool useHE = fgUsedVars[kCosThetaHE] || fgUsedVars[kPhiHE]; // helicity frame bool useCS = fgUsedVars[kCosThetaCS] || fgUsedVars[kPhiCS]; // Collins-Soper frame + bool useGJB1 = fgUsedVars[kCosThetaGJB1] || fgUsedVars[kPhiGJB1]; // Gottfried-Jackson frame, beam 1 (+z) + bool useGJB2 = fgUsedVars[kCosThetaGJB2] || fgUsedVars[kPhiGJB2]; // Gottfried-Jackson frame, beam 2 (-z) bool usePP = fgUsedVars[kCosThetaPP]; // production plane frame bool useRM = fgUsedVars[kCosThetaRM]; // Random frame - if (useHE || useCS || usePP || useRM) { + if (useHE || useCS || useGJB1 || useGJB2 || usePP || useRM) { ROOT::Math::Boost boostv12{v12.BoostToCM()}; ROOT::Math::XYZVectorF v1_CM{(boostv12(v1).Vect()).Unit()}; ROOT::Math::XYZVectorF v2_CM{(boostv12(v2).Vect()).Unit()}; @@ -3948,6 +3962,70 @@ void VarManager::FillPair(T1 const& t1, T2 const& t2, float* values) } } + if (useGJB1) { + // z-axis: beam direction in the mother rest frame + ROOT::Math::XYZVectorF zaxis_GJB1{Beam1_CM.Unit()}; + // mother momentum in LAB (defines the production plane) + ROOT::Math::XYZVectorF motherLab_GJB1{(ROOT::Math::XYZVectorF(v12.Vect())).Unit()}; + // y-axis: normal to the production plane (fallback to (0,1,0) if degenerate) + ROOT::Math::XYZVectorF yaxis_GJB1{zaxis_GJB1.Cross(motherLab_GJB1)}; + if (yaxis_GJB1.Mag2() > 1.e-6f) { + yaxis_GJB1 = yaxis_GJB1.Unit(); + } else { + yaxis_GJB1 = ROOT::Math::XYZVectorF(0.f, 1.f, 0.f); + } + ROOT::Math::XYZVectorF xaxis_GJB1{(yaxis_GJB1.Cross(zaxis_GJB1)).Unit()}; + if (fgUsedVars[kCosThetaGJB1]) { + values[kCosThetaGJB1] = zaxis_GJB1.Dot(v_CM); + } + if (fgUsedVars[kPhiGJB1]) { + values[kPhiGJB1] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_GJB1.Dot(v_CM), xaxis_GJB1.Dot(v_CM))); + } + if (fgUsedVars[kPhiTildeGJB1]) { + if (fgUsedVars[kCosThetaGJB1] && fgUsedVars[kPhiGJB1]) { + if (values[kCosThetaGJB1] > 0) { + values[kPhiTildeGJB1] = RecoDecay::constrainAngle(values[kPhiGJB1] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 + } else { + values[kPhiTildeGJB1] = RecoDecay::constrainAngle(values[kPhiGJB1] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 + } + } else { + values[kPhiTildeGJB1] = -999; // not computable + } + } + } + + if (useGJB2) { + // z-axis: beam direction in the mother rest frame + ROOT::Math::XYZVectorF zaxis_GJB2{Beam2_CM.Unit()}; + // mother momentum in LAB (defines the production plane) + ROOT::Math::XYZVectorF motherLab_GJB2{(ROOT::Math::XYZVectorF(v12.Vect())).Unit()}; + // y-axis: normal to the production plane (fallback to (0,1,0) if degenerate) + ROOT::Math::XYZVectorF yaxis_GJB2{zaxis_GJB2.Cross(motherLab_GJB2)}; + if (yaxis_GJB2.Mag2() > 1.e-6f) { + yaxis_GJB2 = yaxis_GJB2.Unit(); + } else { + yaxis_GJB2 = ROOT::Math::XYZVectorF(0.f, 1.f, 0.f); + } + ROOT::Math::XYZVectorF xaxis_GJB2{(yaxis_GJB2.Cross(zaxis_GJB2)).Unit()}; + if (fgUsedVars[kCosThetaGJB2]) { + values[kCosThetaGJB2] = zaxis_GJB2.Dot(v_CM); + } + if (fgUsedVars[kPhiGJB2]) { + values[kPhiGJB2] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_GJB2.Dot(v_CM), xaxis_GJB2.Dot(v_CM))); + } + if (fgUsedVars[kPhiTildeGJB2]) { + if (fgUsedVars[kCosThetaGJB2] && fgUsedVars[kPhiGJB2]) { + if (values[kCosThetaGJB2] > 0) { + values[kPhiTildeGJB2] = RecoDecay::constrainAngle(values[kPhiGJB2] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 + } else { + values[kPhiTildeGJB2] = RecoDecay::constrainAngle(values[kPhiGJB2] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 + } + } else { + values[kPhiTildeGJB2] = -999; // not computable + } + } + } + if (usePP) { ROOT::Math::XYZVector zaxis_PP = ROOT::Math::XYZVector(v12.Py(), -v12.Px(), 0.f); ROOT::Math::XYZVector yaxis_PP{(v12.Vect()).Unit()}; @@ -4384,10 +4462,12 @@ void VarManager::FillPairME(T1 const& t1, T2 const& t2, float* values) // polarization parameters bool useHE = fgUsedVars[kCosThetaHE] || fgUsedVars[kPhiHE]; // helicity frame bool useCS = fgUsedVars[kCosThetaCS] || fgUsedVars[kPhiCS]; // Collins-Soper frame + bool useGJB1 = fgUsedVars[kCosThetaGJB1] || fgUsedVars[kPhiGJB1]; // Gottfried-Jackson frame, beam 1 (+z) + bool useGJB2 = fgUsedVars[kCosThetaGJB2] || fgUsedVars[kPhiGJB2]; // Gottfried-Jackson frame, beam 2 (-z) bool usePP = fgUsedVars[kCosThetaPP]; // production plane frame bool useRM = fgUsedVars[kCosThetaRM]; // Random frame - if (useHE || useCS || usePP || useRM) { + if (useHE || useCS || useGJB1 || useGJB2 || usePP || useRM) { ROOT::Math::Boost boostv12{v12.BoostToCM()}; ROOT::Math::XYZVectorF v1_CM{(boostv12(v1).Vect()).Unit()}; ROOT::Math::XYZVectorF v2_CM{(boostv12(v2).Vect()).Unit()}; @@ -4443,6 +4523,70 @@ void VarManager::FillPairME(T1 const& t1, T2 const& t2, float* values) } } + if (useGJB1) { + // z-axis: beam direction in the mother rest frame + ROOT::Math::XYZVectorF zaxis_GJB1{Beam1_CM.Unit()}; + // mother momentum in LAB (defines the production plane) + ROOT::Math::XYZVectorF motherLab_GJB1{(ROOT::Math::XYZVectorF(v12.Vect())).Unit()}; + // y-axis: normal to the production plane + ROOT::Math::XYZVectorF yaxis_GJB1{zaxis_GJB1.Cross(motherLab_GJB1)}; + if (yaxis_GJB1.Mag2() > 1.e-6f) { + yaxis_GJB1 = yaxis_GJB1.Unit(); + } else { + yaxis_GJB1 = ROOT::Math::XYZVectorF(0.f, 1.f, 0.f); + } + ROOT::Math::XYZVectorF xaxis_GJB1{(yaxis_GJB1.Cross(zaxis_GJB1)).Unit()}; + if (fgUsedVars[kCosThetaGJB1]) { + values[kCosThetaGJB1] = zaxis_GJB1.Dot(v_CM); + } + if (fgUsedVars[kPhiGJB1]) { + values[kPhiGJB1] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_GJB1.Dot(v_CM), xaxis_GJB1.Dot(v_CM))); + } + if (fgUsedVars[kPhiTildeGJB1]) { + if (fgUsedVars[kCosThetaGJB1] && fgUsedVars[kPhiGJB1]) { + if (values[kCosThetaGJB1] > 0) { + values[kPhiTildeGJB1] = RecoDecay::constrainAngle(values[kPhiGJB1] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 + } else { + values[kPhiTildeGJB1] = RecoDecay::constrainAngle(values[kPhiGJB1] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 + } + } else { + values[kPhiTildeGJB1] = -999; // not computable + } + } + } + + if (useGJB2) { + // z-axis: beam direction in the mother rest frame + ROOT::Math::XYZVectorF zaxis_GJB2{Beam2_CM.Unit()}; + // mother momentum in LAB (defines the production plane) + ROOT::Math::XYZVectorF motherLab_GJB2{(ROOT::Math::XYZVectorF(v12.Vect())).Unit()}; + // y-axis: normal to the production plane (fallback to (0,1,0) if degenerate) + ROOT::Math::XYZVectorF yaxis_GJB2{zaxis_GJB2.Cross(motherLab_GJB2)}; + if (yaxis_GJB2.Mag2() > 1.e-6f) { + yaxis_GJB2 = yaxis_GJB2.Unit(); + } else { + yaxis_GJB2 = ROOT::Math::XYZVectorF(0.f, 1.f, 0.f); + } + ROOT::Math::XYZVectorF xaxis_GJB2{(yaxis_GJB2.Cross(zaxis_GJB2)).Unit()}; + if (fgUsedVars[kCosThetaGJB2]) { + values[kCosThetaGJB2] = zaxis_GJB2.Dot(v_CM); + } + if (fgUsedVars[kPhiGJB2]) { + values[kPhiGJB2] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_GJB2.Dot(v_CM), xaxis_GJB2.Dot(v_CM))); + } + if (fgUsedVars[kPhiTildeGJB2]) { + if (fgUsedVars[kCosThetaGJB2] && fgUsedVars[kPhiGJB2]) { + if (values[kCosThetaGJB2] > 0) { + values[kPhiTildeGJB2] = RecoDecay::constrainAngle(values[kPhiGJB2] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 + } else { + values[kPhiTildeGJB2] = RecoDecay::constrainAngle(values[kPhiGJB2] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 + } + } else { + values[kPhiTildeGJB2] = -999; // not computable + } + } + } + if (usePP) { ROOT::Math::XYZVector zaxis_PP = ROOT::Math::XYZVector(v12.Py(), -v12.Px(), 0.f); ROOT::Math::XYZVector yaxis_PP{(v12.Vect()).Unit()}; @@ -4746,10 +4890,12 @@ void VarManager::FillPairMC(T1 const& t1, T2 const& t2, float* values) // polarization parameters bool useHE = fgUsedVars[kMCCosThetaHE] || fgUsedVars[kMCPhiHE]; // helicity frame bool useCS = fgUsedVars[kMCCosThetaCS] || fgUsedVars[kMCPhiCS]; // Collins-Soper frame + bool useGJB1 = fgUsedVars[kMCCosThetaGJB1] || fgUsedVars[kMCPhiGJB1]; // Gottfried-Jackson frame, beam 1 (+z) + bool useGJB2 = fgUsedVars[kMCCosThetaGJB2] || fgUsedVars[kMCPhiGJB2]; // Gottfried-Jackson frame, beam 2 (-z) bool usePP = fgUsedVars[kMCCosThetaPP]; // production plane frame bool useRM = fgUsedVars[kMCCosThetaRM]; // Random frame - if (useHE || useCS || usePP || useRM) { + if (useHE || useCS || useGJB1 || useGJB2 || usePP || useRM) { ROOT::Math::Boost boostv12{v12.BoostToCM()}; ROOT::Math::XYZVectorF v1_CM{(boostv12(v1).Vect()).Unit()}; ROOT::Math::XYZVectorF v2_CM{(boostv12(v2).Vect()).Unit()}; @@ -4805,6 +4951,70 @@ void VarManager::FillPairMC(T1 const& t1, T2 const& t2, float* values) } } + if (useGJB1) { + // z-axis: beam direction in the mother rest frame + ROOT::Math::XYZVectorF zaxis_GJB1{Beam1_CM.Unit()}; + // mother momentum in LAB (defines the production plane) + ROOT::Math::XYZVectorF motherLab_GJB1{(ROOT::Math::XYZVectorF(v12.Vect())).Unit()}; + // y-axis: normal to the production plane (fallback to (0,1,0) if degenerate) + ROOT::Math::XYZVectorF yaxis_GJB1{zaxis_GJB1.Cross(motherLab_GJB1)}; + if (yaxis_GJB1.Mag2() > 1.e-6f) { + yaxis_GJB1 = yaxis_GJB1.Unit(); + } else { + yaxis_GJB1 = ROOT::Math::XYZVectorF(0.f, 1.f, 0.f); + } + ROOT::Math::XYZVectorF xaxis_GJB1{(yaxis_GJB1.Cross(zaxis_GJB1)).Unit()}; + if (fgUsedVars[kMCCosThetaGJB1]) { + values[kMCCosThetaGJB1] = zaxis_GJB1.Dot(v_CM); + } + if (fgUsedVars[kMCPhiGJB1]) { + values[kMCPhiGJB1] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_GJB1.Dot(v_CM), xaxis_GJB1.Dot(v_CM))); + } + if (fgUsedVars[kMCPhiTildeGJB1]) { + if (fgUsedVars[kMCCosThetaGJB1] && fgUsedVars[kMCPhiGJB1]) { + if (values[kMCCosThetaGJB1] > 0) { + values[kMCPhiTildeGJB1] = RecoDecay::constrainAngle(values[kMCPhiGJB1] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 + } else { + values[kMCPhiTildeGJB1] = RecoDecay::constrainAngle(values[kMCPhiGJB1] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 + } + } else { + values[kMCPhiTildeGJB1] = -999; // not computable + } + } + } + + if (useGJB2) { + // z-axis: beam direction in the mother rest frame + ROOT::Math::XYZVectorF zaxis_GJB2{Beam2_CM.Unit()}; + // mother momentum in LAB (defines the production plane) + ROOT::Math::XYZVectorF motherLab_GJB2{(ROOT::Math::XYZVectorF(v12.Vect())).Unit()}; + // y-axis: normal to the production plane (fallback to (0,1,0) if degenerate) + ROOT::Math::XYZVectorF yaxis_GJB2{zaxis_GJB2.Cross(motherLab_GJB2)}; + if (yaxis_GJB2.Mag2() > 1.e-6f) { + yaxis_GJB2 = yaxis_GJB2.Unit(); + } else { + yaxis_GJB2 = ROOT::Math::XYZVectorF(0.f, 1.f, 0.f); + } + ROOT::Math::XYZVectorF xaxis_GJB2{(yaxis_GJB2.Cross(zaxis_GJB2)).Unit()}; + if (fgUsedVars[kMCCosThetaGJB2]) { + values[kMCCosThetaGJB2] = zaxis_GJB2.Dot(v_CM); + } + if (fgUsedVars[kMCPhiGJB2]) { + values[kMCPhiGJB2] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_GJB2.Dot(v_CM), xaxis_GJB2.Dot(v_CM))); + } + if (fgUsedVars[kMCPhiTildeGJB2]) { + if (fgUsedVars[kMCCosThetaGJB2] && fgUsedVars[kMCPhiGJB2]) { + if (values[kMCCosThetaGJB2] > 0) { + values[kMCPhiTildeGJB2] = RecoDecay::constrainAngle(values[kMCPhiGJB2] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 + } else { + values[kMCPhiTildeGJB2] = RecoDecay::constrainAngle(values[kMCPhiGJB2] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 + } + } else { + values[kMCPhiTildeGJB2] = -999; // not computable + } + } + } + if (usePP) { ROOT::Math::XYZVector zaxis_PP = ROOT::Math::XYZVector(v12.Py(), -v12.Px(), 0.f); ROOT::Math::XYZVector yaxis_PP{v12.Vect().Unit()}; @@ -7559,10 +7769,12 @@ void VarManager::FillPairAlice3(T1 const& t1, T2 const& t2, float* values) // polarization parameters bool useHE = fgUsedVars[kCosThetaHE] || fgUsedVars[kPhiHE]; // helicity frame bool useCS = fgUsedVars[kCosThetaCS] || fgUsedVars[kPhiCS]; // Collins-Soper frame + bool useGJB1 = fgUsedVars[kCosThetaGJB1] || fgUsedVars[kPhiGJB1]; // Gottfried-Jackson frame, beam 1 (+z) + bool useGJB2 = fgUsedVars[kCosThetaGJB2] || fgUsedVars[kPhiGJB2]; // Gottfried-Jackson frame, beam 2 (-z) bool usePP = fgUsedVars[kCosThetaPP]; // production plane frame bool useRM = fgUsedVars[kCosThetaRM]; // Random frame - if (useHE || useCS || usePP || useRM) { + if (useHE || useCS || useGJB1 || useGJB2 || usePP || useRM) { ROOT::Math::Boost boostv12{v12.BoostToCM()}; ROOT::Math::XYZVectorF v1_CM{(boostv12(v1).Vect()).Unit()}; ROOT::Math::XYZVectorF v2_CM{(boostv12(v2).Vect()).Unit()}; @@ -7618,6 +7830,70 @@ void VarManager::FillPairAlice3(T1 const& t1, T2 const& t2, float* values) } } + if (useGJB1) { + // z-axis: beam direction in the mother rest frame + ROOT::Math::XYZVectorF zaxis_GJB1{Beam1_CM.Unit()}; + // mother momentum in LAB (defines the production plane) + ROOT::Math::XYZVectorF motherLab_GJB1{(ROOT::Math::XYZVectorF(v12.Vect())).Unit()}; + // y-axis: normal to the production plane (fallback to (0,1,0) if degenerate) + ROOT::Math::XYZVectorF yaxis_GJB1{zaxis_GJB1.Cross(motherLab_GJB1)}; + if (yaxis_GJB1.Mag2() > 1.e-6f) { + yaxis_GJB1 = yaxis_GJB1.Unit(); + } else { + yaxis_GJB1 = ROOT::Math::XYZVectorF(0.f, 1.f, 0.f); + } + ROOT::Math::XYZVectorF xaxis_GJB1{(yaxis_GJB1.Cross(zaxis_GJB1)).Unit()}; + if (fgUsedVars[kCosThetaGJB1]) { + values[kCosThetaGJB1] = zaxis_GJB1.Dot(v_CM); + } + if (fgUsedVars[kPhiGJB1]) { + values[kPhiGJB1] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_GJB1.Dot(v_CM), xaxis_GJB1.Dot(v_CM))); + } + if (fgUsedVars[kPhiTildeGJB1]) { + if (fgUsedVars[kCosThetaGJB1] && fgUsedVars[kPhiGJB1]) { + if (values[kCosThetaGJB1] > 0) { + values[kPhiTildeGJB1] = RecoDecay::constrainAngle(values[kPhiGJB1] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 + } else { + values[kPhiTildeGJB1] = RecoDecay::constrainAngle(values[kPhiGJB1] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 + } + } else { + values[kPhiTildeGJB1] = -999; // not computable + } + } + } + + if (useGJB2) { + // z-axis: beam direction in the mother rest frame + ROOT::Math::XYZVectorF zaxis_GJB2{Beam2_CM.Unit()}; + // mother momentum in LAB (defines the production plane) + ROOT::Math::XYZVectorF motherLab_GJB2{(ROOT::Math::XYZVectorF(v12.Vect())).Unit()}; + // y-axis: normal to the production plane (fallback to (0,1,0) if degenerate) + ROOT::Math::XYZVectorF yaxis_GJB2{zaxis_GJB2.Cross(motherLab_GJB2)}; + if (yaxis_GJB2.Mag2() > 1.e-6f) { + yaxis_GJB2 = yaxis_GJB2.Unit(); + } else { + yaxis_GJB2 = ROOT::Math::XYZVectorF(0.f, 1.f, 0.f); + } + ROOT::Math::XYZVectorF xaxis_GJB2{(yaxis_GJB2.Cross(zaxis_GJB2)).Unit()}; + if (fgUsedVars[kCosThetaGJB2]) { + values[kCosThetaGJB2] = zaxis_GJB2.Dot(v_CM); + } + if (fgUsedVars[kPhiGJB2]) { + values[kPhiGJB2] = RecoDecay::constrainAngle(TMath::ATan2(yaxis_GJB2.Dot(v_CM), xaxis_GJB2.Dot(v_CM))); + } + if (fgUsedVars[kPhiTildeGJB2]) { + if (fgUsedVars[kCosThetaGJB2] && fgUsedVars[kPhiGJB2]) { + if (values[kCosThetaGJB2] > 0) { + values[kPhiTildeGJB2] = RecoDecay::constrainAngle(values[kPhiGJB2] - o2::constants::math::PIQuarter); // phi_tilde = phi - pi/4 + } else { + values[kPhiTildeGJB2] = RecoDecay::constrainAngle(values[kPhiGJB2] - 3. * o2::constants::math::PIQuarter); // phi_tilde = phi - 3pi/4 + } + } else { + values[kPhiTildeGJB2] = -999; // not computable + } + } + } + if (usePP) { ROOT::Math::XYZVector zaxis_PP = ROOT::Math::XYZVector(v12.Py(), -v12.Px(), 0.f); ROOT::Math::XYZVector yaxis_PP{(v12.Vect()).Unit()}; diff --git a/PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx b/PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx index c35430d09dc..6cfbd56c2e4 100644 --- a/PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx +++ b/PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx @@ -1160,6 +1160,8 @@ struct AnalysisSameEventPairing { Configurable fConfigPseudoRMQA{"cfgPseudoRMQA", false, "If true, fill QA histograms with pseudo-proper decay length analysis with polarization RM"}; Configurable fConfigTruthPbPbMIDYHE{"cfgTruthPbPbmidyHE", false, "If true, fill QA histograms with dielectron pairs in helicity frame"}; Configurable fConfigTruthPbPbMIDYCS{"cfgTruthPbPbmidyCS", false, "If true, fill QA histograms with dielectron pairs in collion-soper frame"}; + Configurable fConfigTruthPbPbMIDYGJB1{"cfgTruthPbPbmidyGJB1", false, "If true, fill QA histograms with dielectron pairs in Gottfried-Jackson frame, beam 1 (+z)"}; + Configurable fConfigTruthPbPbMIDYGJB2{"cfgTruthPbPbmidyGJB2", false, "If true, fill QA histograms with dielectron pairs in Gottfried-Jackson frame, beam 2 (-z)"}; Configurable fConfigAddSEPHistogram{"cfgAddSEPHistogram", "", "Comma separated list of histograms"}; Configurable fConfigAddJSONHistograms{"cfgAddJSONHistograms", "", "Histograms in JSON format"}; Configurable useRemoteField{"cfgUseRemoteField", false, "Chose whether to fetch the magnetic field from ccdb or set it manually"}; @@ -1537,6 +1539,12 @@ struct AnalysisSameEventPairing { if (fConfigOptions.fConfigTruthPbPbMIDYCS.value) { histNames += Form("MCTruthGenPoldielectronPbPbPairCSSel_%s;", sig->GetName()); } + if (fConfigOptions.fConfigTruthPbPbMIDYGJB1.value) { + histNames += Form("MCTruthGenPoldielectronPbPbPairGJB1Sel_%s;", sig->GetName()); + } + if (fConfigOptions.fConfigTruthPbPbMIDYGJB2.value) { + histNames += Form("MCTruthGenPoldielectronPbPbPairGJB2Sel_%s;", sig->GetName()); + } fHasTwoProngGenMCsignals = true; } // for these pair level signals, also add histograms for each MCgenAcc cut if specified @@ -1560,6 +1568,12 @@ struct AnalysisSameEventPairing { if (fConfigOptions.fConfigTruthPbPbMIDYCS.value) { histNames += Form("MCTruthGenPoldielectronPbPbPairCSSel_%s_%s;", sig->GetName(), cut->GetName()); } + if (fConfigOptions.fConfigTruthPbPbMIDYGJB1.value) { + histNames += Form("MCTruthGenPoldielectronPbPbPairGJB1Sel_%s_%s;", sig->GetName(), cut->GetName()); + } + if (fConfigOptions.fConfigTruthPbPbMIDYGJB2.value) { + histNames += Form("MCTruthGenPoldielectronPbPbPairGJB2Sel_%s_%s;", sig->GetName(), cut->GetName()); + } } } } @@ -2239,6 +2253,12 @@ struct AnalysisSameEventPairing { if (fConfigOptions.fConfigTruthPbPbMIDYCS.value) { fHistMan->FillHistClass(Form("MCTruthGenPoldielectronPbPbPairCSSel_%s", sig->GetName()), static_cast(VarManager::fgValues)); } + if (fConfigOptions.fConfigTruthPbPbMIDYGJB1.value) { + fHistMan->FillHistClass(Form("MCTruthGenPoldielectronPbPbPairGJB1Sel_%s", sig->GetName()), static_cast(VarManager::fgValues)); + } + if (fConfigOptions.fConfigTruthPbPbMIDYGJB2.value) { + fHistMan->FillHistClass(Form("MCTruthGenPoldielectronPbPbPairGJB2Sel_%s", sig->GetName()), static_cast(VarManager::fgValues)); + } if (fUseMCGenAccCut) { for (const auto& cut : fMCGenAccCuts) { @@ -2262,6 +2282,12 @@ struct AnalysisSameEventPairing { if (fConfigOptions.fConfigTruthPbPbMIDYCS.value) { fHistMan->FillHistClass(Form("MCTruthGenPoldielectronPbPbPairCSSel_%s_%s", sig->GetName(), cut->GetName()), static_cast(VarManager::fgValues)); } + if (fConfigOptions.fConfigTruthPbPbMIDYGJB1.value) { + fHistMan->FillHistClass(Form("MCTruthGenPoldielectronPbPbPairGJB1Sel_%s_%s", sig->GetName(), cut->GetName()), static_cast(VarManager::fgValues)); + } + if (fConfigOptions.fConfigTruthPbPbMIDYGJB2.value) { + fHistMan->FillHistClass(Form("MCTruthGenPoldielectronPbPbPairGJB2Sel_%s_%s", sig->GetName(), cut->GetName()), static_cast(VarManager::fgValues)); + } } } } @@ -3825,6 +3851,12 @@ void DefineHistograms(HistogramManager* histMan, const TString& histClasses, con if (classStr.Contains("MCTruthGenPoldielectronPbPbPairCS")) { dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "polarization-dielectron-pbpb-midy-cs-gen", histName); } + if (classStr.Contains("MCTruthGenPoldielectronPbPbPairGJB1")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "polarization-dielectron-pbpb-midy-gjb1-gen", histName); + } + if (classStr.Contains("MCTruthGenPoldielectronPbPbPairGJB2")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "polarization-dielectron-pbpb-midy-gjb2-gen", histName); + } if (classStr.Contains("MCTruthGenSelBR")) { dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "mctruth_triple"); } else if (classStr.Contains("MCTruthGen")) { From ae3b9d64d9f77fe51b2ff24b91ce754f790663b6 Mon Sep 17 00:00:00 2001 From: dmallick2 Date: Tue, 13 Oct 2026 14:36:52 +0200 Subject: [PATCH 2/4] Add GJ axis definition and fill histograms --- PWGDQ/Core/VarManager.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index 8e41ab60b9b..c83f57daa7c 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -4528,7 +4528,7 @@ void VarManager::FillPairME(T1 const& t1, T2 const& t2, float* values) ROOT::Math::XYZVectorF zaxis_GJB1{Beam1_CM.Unit()}; // mother momentum in LAB (defines the production plane) ROOT::Math::XYZVectorF motherLab_GJB1{(ROOT::Math::XYZVectorF(v12.Vect())).Unit()}; - // y-axis: normal to the production plane + // y-axis: normal to the production plane ROOT::Math::XYZVectorF yaxis_GJB1{zaxis_GJB1.Cross(motherLab_GJB1)}; if (yaxis_GJB1.Mag2() > 1.e-6f) { yaxis_GJB1 = yaxis_GJB1.Unit(); From 446c05f4909cc0832e9667a3c11a8c19561d4882 Mon Sep 17 00:00:00 2001 From: dmallick2 Date: Tue, 13 Oct 2026 15:36:55 +0200 Subject: [PATCH 3/4] Fix uninitialized members and variable shadowing --- PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx b/PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx index 6cfbd56c2e4..3d48aee1ae9 100644 --- a/PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx +++ b/PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx @@ -776,7 +776,7 @@ struct AnalysisTrackSelection { // compute quantities which depend on the associated collision, such as DCA if (track.collisionId() != event.globalIndex()) { VarManager::FillTrackCollision(track, event); - // cout << "Filled track observables for association" << endl; + // cout << "Filled track observables for association" << endl; } bool isCorrectAssoc = false; @@ -1990,9 +1990,9 @@ struct AnalysisSameEventPairing { if (twoTrackFilter & (static_cast(1) << icut)) { isAmbiInBunch = (twoTrackFilter & (static_cast(1) << 28)) || (twoTrackFilter & (static_cast(1) << 29)); isAmbiOutOfBunch = (twoTrackFilter & (static_cast(1) << 30)) || (twoTrackFilter & (static_cast(1) << 31)); - if (sign1 * sign2 < 0) { // +- pairs + if (sign1 * sign2 < 0) { // +- pairs fHistMan->FillHistClass(histNames[icut][0].Data(), static_cast(VarManager::fgValues)); // reconstructed, unmatched - for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals if (mcDecision & (static_cast(1) << isig)) { PromptNonPromptSepTable(VarManager::fgValues[VarManager::kMass], VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kRap], VarManager::fgValues[VarManager::kPhi], VarManager::fgValues[VarManager::kVertexingTauxyProjected], VarManager::fgValues[VarManager::kVertexingTauxyProjectedPoleJPsiMass], VarManager::fgValues[VarManager::kVertexingTauzProjected], VarManager::fgValues[VarManager::kVertexingTauxyProjectedPoleJPsiMassRecalculatePV], @@ -3250,14 +3250,14 @@ struct AnalysisDileptonTrackTrack { Produces DileptonTrackTrackTable; - int fCurrentRun; // needed to detect if the run changed and trigger update of calibrations etc. + int fCurrentRun = -1; // needed to detect if the run changed and trigger update of calibrations etc. // uint32_t fTrackCutBitMap; // track cut bit mask to be used in the selection of tracks associated with dileptons // cut name setting TString fTrackCutName1; TString fTrackCutName2; std::vector fTrackCutNames; - uint32_t fTrackCutBitMap1; - uint32_t fTrackCutBitMap2; + uint32_t fTrackCutBitMap1 = 0; + uint32_t fTrackCutBitMap2 = 0; bool fIsSameTrackCut = false; AnalysisCompositeCut fDileptonCut; std::vector fQuadrupletCutNames; @@ -3272,8 +3272,8 @@ struct AnalysisDileptonTrackTrack { constexpr static uint32_t fgDileptonFillMap = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::Pair; // fill map - float* fValuesQuadruplet; - HistogramManager* fHistMan; + float* fValuesQuadruplet = nullptr; + HistogramManager* fHistMan = nullptr; // o2::analysis::QuadMlResponse fQuadMlResponse; // std::vector fOutputMlQuad; @@ -3473,7 +3473,7 @@ struct AnalysisDileptonTrackTrack { } uint32_t mcDecision = static_cast(0); - size_t isig = 0; + // size_t isig = 0; for (auto dilepton : dileptons) { // get full track info of tracks based on the index @@ -3554,7 +3554,7 @@ struct AnalysisDileptonTrackTrack { auto track1MC = track1.mcParticle(); auto track2MC = track2.mcParticle(); mcDecision = 0; - isig = 0; + size_t isig = 0; for (auto sig = fRecMCSignals.begin(); sig != fRecMCSignals.end(); sig++, isig++) { if ((*sig)->CheckSignal(true, lepton1MC, lepton2MC, track1MC, track2MC)) { mcDecision |= (static_cast(1) << isig); From d0e1da4fa40240a5bdcd5b81e774f7c5dc645bab Mon Sep 17 00:00:00 2001 From: dmallick2 Date: Tue, 13 Oct 2026 15:48:19 +0200 Subject: [PATCH 4/4] Fix uninitialized members and variable shadowing --- PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx b/PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx index 3d48aee1ae9..424d809dbe1 100644 --- a/PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx +++ b/PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx @@ -3554,13 +3554,12 @@ struct AnalysisDileptonTrackTrack { auto track1MC = track1.mcParticle(); auto track2MC = track2.mcParticle(); mcDecision = 0; - size_t isig = 0; - for (auto sig = fRecMCSignals.begin(); sig != fRecMCSignals.end(); sig++, isig++) { - if ((*sig)->CheckSignal(true, lepton1MC, lepton2MC, track1MC, track2MC)) { - mcDecision |= (static_cast(1) << isig); + // size_t isig = 0; + for (size_t iSignal = 0; iSignal < fRecMCSignals.size(); ++iSignal) { + if (fRecMCSignals[iSignal]->CheckSignal(true, lepton1MC, lepton2MC, track1MC, track2MC)) { + mcDecision |= (static_cast(1) << iSignal); } } - int iCut = 0; uint32_t CutDecision = 0; for (auto cutname = fQuadrupletCutNames.begin(); cutname != fQuadrupletCutNames.end(); cutname++, iCut++) {