diff --git a/PWGDQ/Core/HistogramsLibrary.cxx b/PWGDQ/Core/HistogramsLibrary.cxx index b65c12de647..96227af31ce 100644 --- a/PWGDQ/Core/HistogramsLibrary.cxx +++ b/PWGDQ/Core/HistogramsLibrary.cxx @@ -972,14 +972,34 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h hm->AddHistogram(histClass, "VtxZ_VtxNcontribReal", "VtxZ vs VtxNcontribReal", false, 240, -12.0, 12.0, VarManager::kVtxZ, 200, 0, 200.0, VarManager::kVtxNcontribReal); } if (subGroupStr.Contains("polarization")) { - hm->AddHistogram(histClass, "cosThetaHE", "", false, 100, -1., 1., VarManager::kCosThetaHE); - hm->AddHistogram(histClass, "cosThetaCS", "", false, 100, -1., 1., VarManager::kCosThetaCS); - hm->AddHistogram(histClass, "PhiHE", "", false, 100, -o2::constants::math::PI, o2::constants::math::PI, VarManager::kPhiHE); - hm->AddHistogram(histClass, "PhiCS", "", false, 100, -o2::constants::math::PI, o2::constants::math::PI, VarManager::kPhiCS); - hm->AddHistogram(histClass, "Mass_Pt_cosThetaHE", "", false, 100, 1.0, 5.0, VarManager::kMass, 250, 0.0, 25.0, VarManager::kPt, 40, -1., 1., VarManager::kCosThetaHE); - hm->AddHistogram(histClass, "Mass_Pt_cosThetaCS", "", false, 100, 1.0, 5.0, VarManager::kMass, 250, 0.0, 25.0, VarManager::kPt, 40, -1., 1., VarManager::kCosThetaCS); - hm->AddHistogram(histClass, "Mass_Pt_PhiHE", "", false, 100, 1.0, 5.0, VarManager::kMass, 250, 0.0, 25.0, VarManager::kPt, 40, -o2::constants::math::PI, o2::constants::math::PI, VarManager::kPhiHE); - hm->AddHistogram(histClass, "Mass_Pt_PhiCS", "", false, 100, 1.0, 5.0, VarManager::kMass, 250, 0.0, 25.0, VarManager::kPt, 40, -o2::constants::math::PI, o2::constants::math::PI, VarManager::kPhiCS); + if (subGroupStr.Contains("helicity")) { + hm->AddHistogram(histClass, "cosThetaHE", "", false, 100, -1., 1., VarManager::kCosThetaHE); + hm->AddHistogram(histClass, "phiHE", "", false, 100, 0, 2 * o2::constants::math::PI, VarManager::kPhiHE); + hm->AddHistogram(histClass, "phitildeHE", "", false, 100, 0, 2 * o2::constants::math::PI, VarManager::kPhiTildeHE); + hm->AddHistogram(histClass, "Mass_Pt_CosThetaHE", "", false, 100, 1.0, 5.0, VarManager::kMass, 40, 0.0, 20.0, VarManager::kPt, 20, -1., 1., VarManager::kCosThetaHE); + hm->AddHistogram(histClass, "Mass_Pt_PhiHE", "", false, 100, 1.0, 5.0, VarManager::kMass, 40, 0.0, 20.0, VarManager::kPt, 20, 0., 2 * o2::constants::math::PI, VarManager::kPhiHE); + hm->AddHistogram(histClass, "Mass_Pt_PhiTildeHE", "", false, 100, 1.0, 5.0, VarManager::kMass, 40, 0.0, 20.0, VarManager::kPt, 20, 0., 2 * o2::constants::math::PI, VarManager::kPhiTildeHE); + } + if (subGroupStr.Contains("collins-soper")) { + hm->AddHistogram(histClass, "cosThetaCS", "", false, 100, -1., 1., VarManager::kCosThetaCS); + hm->AddHistogram(histClass, "phiCS", "", false, 100, 0, 2 * o2::constants::math::PI, VarManager::kPhiCS); + hm->AddHistogram(histClass, "phitildeCS", "", false, 100, 0, 2 * o2::constants::math::PI, VarManager::kPhiTildeCS); + hm->AddHistogram(histClass, "Mass_Pt_CosThetaCS", "", false, 100, 1.0, 5.0, VarManager::kMass, 40, 0.0, 20.0, VarManager::kPt, 20, -1., 1., VarManager::kCosThetaCS); + hm->AddHistogram(histClass, "Mass_Pt_PhiCS", "", false, 100, 1.0, 5.0, VarManager::kMass, 40, 0.0, 20.0, VarManager::kPt, 20, 0., 2 * o2::constants::math::PI, VarManager::kPhiCS); + hm->AddHistogram(histClass, "Mass_Pt_PhiTildeCS", "", false, 100, 1.0, 5.0, VarManager::kMass, 40, 0.0, 20.0, VarManager::kPt, 20, 0., 2 * o2::constants::math::PI, VarManager::kPhiTildeCS); + } + if (subGroupStr.Contains("production")) { + hm->AddHistogram(histClass, "cosThetaPP", "", false, 100, -1., 1., VarManager::kCosThetaPP); + hm->AddHistogram(histClass, "phiPP", "", false, 100, 0, 2 * o2::constants::math::PI, VarManager::kPhiPP); + hm->AddHistogram(histClass, "phitildePP", "", false, 100, 0, 2 * o2::constants::math::PI, VarManager::kPhiTildePP); + hm->AddHistogram(histClass, "Mass_Pt_CosThetaPP", "", false, 100, 1.0, 5.0, VarManager::kMass, 40, 0.0, 20.0, VarManager::kPt, 20, -1., 1., VarManager::kCosThetaPP); + hm->AddHistogram(histClass, "Mass_Pt_PhiPP", "", false, 100, 1.0, 5.0, VarManager::kMass, 40, 0.0, 20.0, VarManager::kPt, 20, 0., 2 * o2::constants::math::PI, VarManager::kPhiPP); + hm->AddHistogram(histClass, "Mass_Pt_PhiTildePP", "", false, 100, 1.0, 5.0, VarManager::kMass, 40, 0.0, 20.0, VarManager::kPt, 20, 0., 2 * o2::constants::math::PI, VarManager::kPhiTildePP); + } + if (subGroupStr.Contains("random")) { + hm->AddHistogram(histClass, "cosThetaRM", "", false, 100, -1., 1., VarManager::kCosThetaRM); + hm->AddHistogram(histClass, "Mass_Pt_CosThetaRM", "", false, 200, 1.0, 5.0, VarManager::kMass, 40, 0.0, 20.0, VarManager::kPt, 20, -1., 1., VarManager::kCosThetaRM); + } } if (subGroupStr.Contains("upsilon")) { hm->AddHistogram(histClass, "MassUpsilon_Pt", "", false, 500, 7.0, 12.0, VarManager::kMass, 400, 0.0, 40.0, VarManager::kPt); diff --git a/PWGDQ/Core/VarManager.cxx b/PWGDQ/Core/VarManager.cxx index dc3c963be43..503101d303f 100644 --- a/PWGDQ/Core/VarManager.cxx +++ b/PWGDQ/Core/VarManager.cxx @@ -542,18 +542,26 @@ void VarManager::SetDefaultVarNames() fgVariableUnits[kMCVy] = "cm"; // TODO: check the unit fgVariableNames[kMCVz] = "MC vz"; fgVariableUnits[kMCVz] = "cm"; // TODO: check the unit - fgVariableNames[kMCCosThetaHE] = "cos#it{#theta}^{MC}_{HE}"; + fgVariableNames[kMCCosThetaHE] = "MC cos(#theta_{HE})"; fgVariableUnits[kMCCosThetaHE] = ""; - fgVariableNames[kMCPhiHE] = "#varphi^{MC}_{HE}"; - fgVariableUnits[kMCPhiHE] = "rad."; - fgVariableNames[kMCCosThetaCS] = "cos#it{#theta}^{MC}_{CS}"; + fgVariableNames[kMCPhiHE] = "MC #varphi_{HE}"; + fgVariableUnits[kMCPhiHE] = "rad"; + fgVariableNames[kMCPhiTildeHE] = "MC #tilde{#varphi}_{HE}"; + fgVariableUnits[kMCPhiTildeHE] = "rad"; + fgVariableNames[kMCCosThetaCS] = "MC cos(#theta_{CS})"; fgVariableUnits[kMCCosThetaCS] = ""; - fgVariableNames[kMCPhiCS] = "#varphi^{MC}_{CS}"; - fgVariableUnits[kMCPhiCS] = "rad."; - fgVariableNames[kMCCosThetaPP] = "cos#it{#theta}^{MC}_{PP}"; + fgVariableNames[kMCPhiCS] = "MC #varphi_{CS}"; + fgVariableUnits[kMCPhiCS] = "rad"; + fgVariableNames[kMCPhiTildeCS] = "MC #tilde{#varphi}_{CS}"; + fgVariableUnits[kMCPhiTildeCS] = "rad"; + fgVariableNames[kMCCosThetaPP] = "MC cos(#theta_{PP})"; fgVariableUnits[kMCCosThetaPP] = ""; - fgVariableNames[kMCPhiPP] = "#varphi^{MC}_{PP}"; - fgVariableUnits[kMCPhiPP] = "rad."; + fgVariableNames[kMCPhiPP] = "MC #varphi_{PP}"; + fgVariableUnits[kMCPhiPP] = "rad"; + fgVariableNames[kMCPhiTildePP] = "MC #tilde{#varphi}_{PP}"; + fgVariableUnits[kMCPhiTildePP] = "rad"; + fgVariableNames[kMCCosThetaRM] = "MC cos(#theta_{RM})"; + fgVariableUnits[kMCCosThetaRM] = ""; fgVariableNames[kCandidateId] = ""; fgVariableUnits[kCandidateId] = ""; fgVariableNames[kPairType] = "Pair type"; @@ -934,18 +942,26 @@ void VarManager::SetDefaultVarNames() fgVariableUnits[kDeltaPhi] = "rad."; fgVariableNames[kDeltaPhiSym] = "#Delta#phi"; fgVariableUnits[kDeltaPhiSym] = "rad."; - fgVariableNames[kCosThetaHE] = "cos#it{#theta}_{HE}"; + fgVariableNames[kCosThetaHE] = "cos#it{#theta}"; fgVariableUnits[kCosThetaHE] = ""; fgVariableNames[kPhiHE] = "#varphi_{HE}"; fgVariableUnits[kPhiHE] = "rad."; + fgVariableNames[kPhiTildeHE] = "#tilde{#varphi}_{HE}"; + fgVariableUnits[kPhiTildeHE] = "rad."; fgVariableNames[kCosThetaCS] = "cos#it{#theta}_{CS}"; fgVariableUnits[kCosThetaCS] = ""; fgVariableNames[kPhiCS] = "#varphi_{CS}"; fgVariableUnits[kPhiCS] = "rad."; + fgVariableNames[kPhiTildeCS] = "#tilde{#varphi}_{CS}"; + fgVariableUnits[kPhiTildeCS] = "rad."; fgVariableNames[kCosThetaPP] = "cos#it{#theta}_{PP}"; fgVariableUnits[kCosThetaPP] = ""; fgVariableNames[kPhiPP] = "#varphi_{PP}"; fgVariableUnits[kPhiPP] = "rad."; + fgVariableNames[kPhiTildePP] = "#tilde{#varphi}_{PP}"; + fgVariableUnits[kPhiTildePP] = "rad."; + fgVariableNames[kCosThetaRM] = "cos#it{#theta}_{RM}"; + fgVariableUnits[kCosThetaRM] = ""; fgVariableNames[kCosPhiVP] = "cos#it{#varphi}_{VP}"; fgVariableUnits[kCosPhiVP] = ""; fgVariableNames[kPhiVP] = "#varphi_{VP} - #Psi_{2}"; @@ -1462,11 +1478,15 @@ void VarManager::SetDefaultVarNames() fgVarNamesMap["kMCEta"] = kMCEta; fgVarNamesMap["kMCY"] = kMCY; fgVarNamesMap["kMCCosThetaHE"] = kMCCosThetaHE; - fgVarNamesMap["kMCCosThetaCS"] = kMCCosThetaCS; - fgVarNamesMap["kMCCosThetaPP"] = kMCCosThetaPP; fgVarNamesMap["kMCPhiHE"] = kMCPhiHE; + fgVarNamesMap["kMCPhiTildeHE"] = kMCPhiTildeHE; + fgVarNamesMap["kMCCosThetaCS"] = kMCCosThetaCS; fgVarNamesMap["kMCPhiCS"] = kMCPhiCS; + fgVarNamesMap["kMCPhiTildeCS"] = kMCPhiTildeCS; + fgVarNamesMap["kMCCosThetaPP"] = kMCCosThetaPP; fgVarNamesMap["kMCPhiPP"] = kMCPhiPP; + fgVarNamesMap["kMCPhiTildePP"] = kMCPhiTildePP; + fgVarNamesMap["kMCCosThetaRM"] = kMCCosThetaRM; fgVarNamesMap["kMCParticleGeneratorId"] = kMCParticleGeneratorId; fgVarNamesMap["kNMCParticleVariables"] = kNMCParticleVariables; fgVarNamesMap["kMCMotherPdgCode"] = kMCMotherPdgCode; @@ -1496,11 +1516,15 @@ void VarManager::SetDefaultVarNames() fgVarNamesMap["kVertexingProcCode"] = kVertexingProcCode; fgVarNamesMap["kVertexingChi2PCA"] = kVertexingChi2PCA; fgVarNamesMap["kCosThetaHE"] = kCosThetaHE; - fgVarNamesMap["kCosThetaCS"] = kCosThetaCS; - fgVarNamesMap["kCosThetaPP"] = kCosThetaPP; fgVarNamesMap["kPhiHE"] = kPhiHE; + fgVarNamesMap["kPhiTildeHE"] = kPhiTildeHE; + fgVarNamesMap["kCosThetaCS"] = kCosThetaCS; fgVarNamesMap["kPhiCS"] = kPhiCS; + fgVarNamesMap["kPhiTildeCS"] = kPhiTildeCS; + fgVarNamesMap["kCosThetaPP"] = kCosThetaPP; fgVarNamesMap["kPhiPP"] = kPhiPP; + fgVarNamesMap["kPhiTildePP"] = kPhiTildePP; + fgVarNamesMap["kCosThetaRM"] = kCosThetaRM; fgVarNamesMap["kCosPhiVP"] = kCosPhiVP; fgVarNamesMap["kPhiVP"] = kPhiVP; fgVarNamesMap["kDeltaPhiPair2"] = kDeltaPhiPair2; diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index b1e18ab7635..1da8fa30562 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -621,11 +621,15 @@ class VarManager : public TObject // MC pair variables kMCCosThetaHE, - kMCCosThetaCS, - kMCCosThetaPP, kMCPhiHE, + kMCPhiTildeHE, + kMCCosThetaCS, kMCPhiCS, + kMCPhiTildeCS, + kMCCosThetaPP, kMCPhiPP, + kMCPhiTildePP, + kMCCosThetaRM, // Pair variables kCandidateId, @@ -654,11 +658,15 @@ class VarManager : public TObject kVertexingProcCode, kVertexingChi2PCA, kCosThetaHE, - kCosThetaCS, - kCosThetaPP, kPhiHE, + kPhiTildeHE, + kCosThetaCS, kPhiCS, + kPhiTildeCS, + kCosThetaPP, kPhiPP, + kPhiTildePP, + kCosThetaRM, kCosPhiVP, kPhiVP, kDeltaPhiPair2, @@ -2814,53 +2822,128 @@ void VarManager::FillPair(T1 const& t1, T2 const& t2, float* values) } } - // TO DO: get the correct values from CCDB - double BeamMomentum = TMath::Sqrt(fgCenterOfMassEnergy * fgCenterOfMassEnergy / 4 - fgMassofCollidingParticle * fgMassofCollidingParticle); // GeV - ROOT::Math::PxPyPzEVector Beam1(0., 0., -BeamMomentum, fgCenterOfMassEnergy / 2); - ROOT::Math::PxPyPzEVector Beam2(0., 0., BeamMomentum, fgCenterOfMassEnergy / 2); - - // Boost to center of mass frame - ROOT::Math::Boost boostv12{v12.BoostToCM()}; - ROOT::Math::XYZVectorF v1_CM{(boostv12(v1).Vect()).Unit()}; - ROOT::Math::XYZVectorF v2_CM{(boostv12(v2).Vect()).Unit()}; - ROOT::Math::XYZVectorF Beam1_CM{(boostv12(Beam1).Vect()).Unit()}; - ROOT::Math::XYZVectorF Beam2_CM{(boostv12(Beam2).Vect()).Unit()}; - - // Helicity frame - 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()}; - - // Collins-Soper frame - ROOT::Math::XYZVectorF zaxis_CS{((Beam1_CM.Unit() - Beam2_CM.Unit()).Unit())}; - ROOT::Math::XYZVectorF yaxis_CS{(Beam1_CM.Cross(Beam2_CM)).Unit()}; - ROOT::Math::XYZVectorF xaxis_CS{(yaxis_CS.Cross(zaxis_CS)).Unit()}; - - // Production frame - ROOT::Math::XYZVector normalVec = ROOT::Math::XYZVector(v12.Py(), -v12.Px(), 0.f); - - if (fgUsedVars[kCosThetaHE]) { - values[kCosThetaHE] = (t1.sign() > 0 ? zaxis_HE.Dot(v1_CM) : zaxis_HE.Dot(v2_CM)); - } - - if (fgUsedVars[kPhiHE]) { - values[kPhiHE] = (t1.sign() > 0 ? TMath::ATan2(yaxis_HE.Dot(v1_CM), xaxis_HE.Dot(v1_CM)) : TMath::ATan2(yaxis_HE.Dot(v2_CM), xaxis_HE.Dot(v2_CM))); - } - - if (fgUsedVars[kCosThetaCS]) { - values[kCosThetaCS] = (t1.sign() > 0 ? zaxis_CS.Dot(v1_CM) : zaxis_CS.Dot(v2_CM)); - } + // polarization parameters + bool useHE = fgUsedVars[kCosThetaHE] || fgUsedVars[kPhiHE]; // helicity frame + bool useCS = fgUsedVars[kCosThetaCS] || fgUsedVars[kPhiCS]; // Collins-Soper frame + bool usePP = fgUsedVars[kCosThetaPP]; // production plane frame + bool useRM = fgUsedVars[kCosThetaRM]; // Random frame + + if (useHE || useCS || usePP || useRM) { + // TO DO: get the correct values from CCDB + double BeamMomentum = TMath::Sqrt(fgCenterOfMassEnergy * fgCenterOfMassEnergy / 4 - fgMassofCollidingParticle * fgMassofCollidingParticle); // GeV + ROOT::Math::PxPyPzEVector Beam1(0., 0., -BeamMomentum, fgCenterOfMassEnergy / 2); + ROOT::Math::PxPyPzEVector Beam2(0., 0., BeamMomentum, fgCenterOfMassEnergy / 2); + + ROOT::Math::Boost boostv12{v12.BoostToCM()}; + ROOT::Math::XYZVectorF v1_CM{(boostv12(v1).Vect()).Unit()}; + ROOT::Math::XYZVectorF v2_CM{(boostv12(v2).Vect()).Unit()}; + ROOT::Math::XYZVectorF Beam1_CM{(boostv12(Beam1).Vect()).Unit()}; + ROOT::Math::XYZVectorF Beam2_CM{(boostv12(Beam2).Vect()).Unit()}; + + // using positive sign convention for the first track + ROOT::Math::XYZVectorF v_CM = (t1.sign() > 0 ? v1_CM : v2_CM); + + if (useHE) { + 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]) + 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] + } + } + 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] + } + } 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] + } + } + } else { + values[kPhiTildeHE] = -999; // not computable + } + } + } - if (fgUsedVars[kPhiCS]) { - values[kPhiCS] = (t1.sign() > 0 ? TMath::ATan2(yaxis_CS.Dot(v1_CM), xaxis_CS.Dot(v1_CM)) : TMath::ATan2(yaxis_CS.Dot(v2_CM), xaxis_CS.Dot(v2_CM))); - } + if (useCS) { + 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]) + 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] + } + } + 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] + } + } 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] + } + } + } else { + values[kPhiTildeCS] = -999; // not computable + } + } + } - if (fgUsedVars[kCosThetaPP]) { - values[kCosThetaPP] = (t1.sign() > 0 ? normalVec.Dot(v1_CM) : normalVec.Dot(v2_CM)); - } + if (usePP) { + ROOT::Math::XYZVector zaxis_PP = ROOT::Math::XYZVector(v12.Py(), -v12.Px(), 0.f); + ROOT::Math::XYZVector yaxis_PP{(v12.Vect()).Unit()}; + ROOT::Math::XYZVector xaxis_PP{(yaxis_PP.Cross(zaxis_PP)).Unit()}; + if (fgUsedVars[kCosThetaPP]) { + values[kCosThetaPP] = zaxis_PP.Dot(v_CM); + } + 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] + } + } + 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] + } + } 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] + } + } + } else { + values[kPhiTildePP] = -999; // not computable + } + } + } - if (fgUsedVars[kPhiPP]) { - values[kPhiPP] = (t1.sign() > 0 ? TMath::ATan2((normalVec.Dot(v1_CM)), zaxis_HE.Dot(v1_CM)) : TMath::ATan2((normalVec.Dot(v2_CM)), zaxis_HE.Dot(v2_CM))); + if (useRM) { + double randomCostheta = gRandom->Uniform(-1., 1.); + double randomPhi = gRandom->Uniform(0., 2. * TMath::Pi()); + 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]) + values[kCosThetaRM] = zaxis_RM.Dot(v_CM); + } } if constexpr ((pairType == kDecayToEE) && ((fillMap & TrackCov) > 0 || (fillMap & ReducedTrackBarrelCov) > 0)) { @@ -3253,38 +3336,129 @@ void VarManager::FillPairMC(T1 const& t1, T2 const& t2, float* values) values[kMCEta] = v12.Eta(); values[kMCPhi] = v12.Phi(); values[kMCY] = -v12.Rapidity(); - double BeamMomentum = TMath::Sqrt(fgCenterOfMassEnergy * fgCenterOfMassEnergy / 4 - fgMassofCollidingParticle * fgMassofCollidingParticle); // GeV - ROOT::Math::PxPyPzEVector Beam1(0., 0., -BeamMomentum, fgCenterOfMassEnergy / 2); - ROOT::Math::PxPyPzEVector Beam2(0., 0., BeamMomentum, fgCenterOfMassEnergy / 2); - - // Boost to center of mass frame - ROOT::Math::Boost boostv12{v12.BoostToCM()}; - ROOT::Math::XYZVectorF v1_CM{(boostv12(v1).Vect()).Unit()}; - ROOT::Math::XYZVectorF v2_CM{(boostv12(v2).Vect()).Unit()}; - ROOT::Math::XYZVectorF Beam1_CM{(boostv12(Beam1).Vect()).Unit()}; - ROOT::Math::XYZVectorF Beam2_CM{(boostv12(Beam2).Vect()).Unit()}; - - if (fgUsedVars[kMCCosThetaHE] || fgUsedVars[kMCPhiHE]) { - 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()}; - values[kMCCosThetaHE] = (t1.pdgCode() < 0 ? zaxis_HE.Dot(v1_CM) : zaxis_HE.Dot(v2_CM)); - values[kMCPhiHE] = (t1.pdgCode() < 0 ? TMath::ATan2(yaxis_HE.Dot(v1_CM), xaxis_HE.Dot(v1_CM)) : TMath::ATan2(yaxis_HE.Dot(v2_CM), xaxis_HE.Dot(v2_CM))); - } - - if (fgUsedVars[kMCCosThetaCS] || fgUsedVars[kMCPhiCS]) { - ROOT::Math::XYZVectorF zaxis_CS{((Beam1_CM.Unit() - Beam2_CM.Unit()).Unit())}; - ROOT::Math::XYZVectorF yaxis_CS{(Beam1_CM.Cross(Beam2_CM)).Unit()}; - ROOT::Math::XYZVectorF xaxis_CS{(yaxis_CS.Cross(zaxis_CS)).Unit()}; - values[kMCCosThetaCS] = (t1.pdgCode() < 0 ? zaxis_CS.Dot(v1_CM) : zaxis_CS.Dot(v2_CM)); - values[kMCPhiCS] = (t1.pdgCode() < 0 ? TMath::ATan2(yaxis_CS.Dot(v1_CM), xaxis_CS.Dot(v1_CM)) : TMath::ATan2(yaxis_CS.Dot(v2_CM), xaxis_CS.Dot(v2_CM))); - } - - if (fgUsedVars[kMCCosThetaPP] || fgUsedVars[kMCPhiPP]) { - ROOT::Math::XYZVectorF zaxis_HE{(v12.Vect()).Unit()}; - ROOT::Math::XYZVector normalVec = ROOT::Math::XYZVector(v12.Py(), -v12.Px(), 0.f); - values[kMCCosThetaPP] = (t1.pdgCode() < 0 ? normalVec.Dot(v1_CM) : normalVec.Dot(v2_CM)); - values[kMCPhiPP] = (t1.pdgCode() < 0 ? TMath::ATan2((normalVec.Dot(v1_CM)), zaxis_HE.Dot(v1_CM)) : TMath::ATan2((normalVec.Dot(v2_CM)), zaxis_HE.Dot(v2_CM))); + + // polarization parameters + bool useHE = fgUsedVars[kMCCosThetaHE] || fgUsedVars[kMCPhiHE]; // helicity frame + bool useCS = fgUsedVars[kMCCosThetaCS] || fgUsedVars[kMCPhiCS]; // Collins-Soper frame + bool usePP = fgUsedVars[kMCCosThetaPP]; // production plane frame + bool useRM = fgUsedVars[kMCCosThetaRM]; // Random frame + + if (useHE || useCS || usePP || useRM) { + // TO DO: get the correct values from CCDB + double BeamMomentum = TMath::Sqrt(fgCenterOfMassEnergy * fgCenterOfMassEnergy / 4 - fgMassofCollidingParticle * fgMassofCollidingParticle); // GeV + ROOT::Math::PxPyPzEVector Beam1(0., 0., -BeamMomentum, fgCenterOfMassEnergy / 2); + ROOT::Math::PxPyPzEVector Beam2(0., 0., BeamMomentum, fgCenterOfMassEnergy / 2); + + ROOT::Math::Boost boostv12{v12.BoostToCM()}; + ROOT::Math::XYZVectorF v1_CM{(boostv12(v1).Vect()).Unit()}; + ROOT::Math::XYZVectorF v2_CM{(boostv12(v2).Vect()).Unit()}; + ROOT::Math::XYZVectorF Beam1_CM{(boostv12(Beam1).Vect()).Unit()}; + ROOT::Math::XYZVectorF Beam2_CM{(boostv12(Beam2).Vect()).Unit()}; + + // using positive sign convention for the first track + ROOT::Math::XYZVectorF v_CM = (t1.pdgCode() > 0 ? v1_CM : v2_CM); + + if (useHE) { + 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]) + 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] + } + } + 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] + } + } 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] + } + } + } else { + values[kMCPhiTildeHE] = -999; // not computable + } + } + } + + if (useCS) { + 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]) + 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] + } + } + 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] + } + } 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] + } + } + } else { + values[kMCPhiTildeCS] = -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()}; + ROOT::Math::XYZVector xaxis_PP{(yaxis_PP.Cross(zaxis_PP)).Unit()}; + if (fgUsedVars[kMCCosThetaPP]) { + 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] + } + } + 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] + } + } 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] + } + } + } else { + values[kMCPhiTildePP] = -999; // not computable + } + } + } + + if (useRM) { + double randomCostheta = gRandom->Uniform(-1., 1.); + double randomPhi = gRandom->Uniform(0., 2. * TMath::Pi()); + 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]) + values[kMCCosThetaRM] = zaxis_RM.Dot(v_CM); + } } }