Skip to content

[PWGDQ] add weighted method #16201

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
zjxiongOvO wants to merge 2 commits into from
Closed

[PWGDQ] add weighted method #16201

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8ec221c
update dqefficiency for PbPb analysis
Apr 29, 2026
40a7845
add weighted method
May 11, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 33 additions & 11 deletions PWGDQ/Core/HistogramsLibrary.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,6 +111,9 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h
hm->AddHistogram(histClass, "CentFT0C_MultTPC", "CentFT0C vs MultTPC", false, 100, 0., 100., VarManager::kCentFT0C, 100, 0., 50000., VarManager::kMultTPC);
hm->AddHistogram(histClass, "CentFT0C_Run", "Cent FT0C", true, 1, -0.5, 0.5, VarManager::kRunNo, 100, 0., 100., VarManager::kCentFT0C, 1, 0, 1, VarManager::kNothing, "", "", "", VarManager::kNothing, VarManager::kNothing, false, true);
}
if (subGroupStr.Contains("randomplane")) {
hm->AddHistogram(histClass, "random plane angle", "random Psi", false, 100, -TMath::Pi() / 2, TMath::Pi() / 2, VarManager::kRandomPsi2);
}
if (subGroupStr.Contains("mult")) {
if (subGroupStr.Contains("pp")) {
hm->AddHistogram(histClass, "MultTPC", "MultTPC", false, 250, 0.0, 500.0, VarManager::kMultTPC);
Expand Down Expand Up @@ -246,6 +249,8 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h
hm->AddHistogram(histClass, "MCVtxZ_VtxZ", "Vtx Z (MC vs rec)", false, 75, -15.0, 15.0, VarManager::kVtxZ, 75, -15.0, 15.0, VarManager::kMCVtxZ);
hm->AddHistogram(histClass, "MCVtxZ", "Vtx Z (MC)", false, 75, -15.0, 15.0, VarManager::kMCVtxZ);
hm->AddHistogram(histClass, "MCImpPar_CentVZERO", "MC impact param vs CentVZERO", false, 50, 0.0, 100.0, VarManager::kCentVZERO, 20, 0.0, 20.0, VarManager::kMCEventImpParam);
hm->AddHistogram(histClass, "MCEventPlaneAngle", "MC event plane angle", false, 100, -TMath::Pi(), TMath::Pi(), VarManager::kMCEventPlaneAngle);
hm->AddHistogram(histClass, "kMCEventCentrFT0C_CentFT0C", "MC Centrality FT0C vs Centrality FT0C", false, 100, 0.0, 100.0, VarManager::kCentFT0C, 100, 0.0, 100.0, VarManager::kMCEventCentrFT0C);
}
if (subGroupStr.Contains("generator")) {
hm->AddHistogram(histClass, "MCVtxX", "Vtx X", false, 1000, -0.5, 0.5, VarManager::kMCVtxX);
Expand Down Expand Up @@ -1161,13 +1166,21 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h
hm->AddHistogram(histClass, "Eta_Pt", "", false, 40, -2.0, 2.0, VarManager::kMCEta, 200, 0.0, 20.0, VarManager::kMCPt);
hm->AddHistogram(histClass, "Phi_Eta", "#phi vs #eta distribution", false, 200, -5.0, 5.0, VarManager::kMCEta, 200, -2. * o2::constants::math::PI, 2. * o2::constants::math::PI, VarManager::kMCPhi);
if (subGroupStr.Contains("polarization")) {
int varspTHE[4] = {VarManager::kMCPt, VarManager::kMCCosThetaHE, VarManager::kMCPhiHE, VarManager::kMCPhiTildeHE};
int varspTCS[4] = {VarManager::kMCPt, VarManager::kMCCosThetaCS, VarManager::kMCPhiCS, VarManager::kMCPhiTildeCS};
int bins[4] = {20, 20, 20, 20};
double xmin[4] = {0., -1., 0., 0.};
double xmax[4] = {20., 1., 2. * o2::constants::math::PI, 2. * o2::constants::math::PI};
hm->AddHistogram(histClass, "Pt_cosThetaHE_phiHE_phiTildeHE", "", 4, varspTHE, bins, xmin, xmax, 0, -1, kFALSE);
hm->AddHistogram(histClass, "Pt_cosThetaCS_phiCS_phiTildeCS", "", 4, varspTCS, bins, xmin, xmax, 0, -1, kFALSE);
if (subGroupStr.Contains("pp")) {
int varspTHE[4] = {VarManager::kMCPt, VarManager::kMCCosThetaHE, VarManager::kMCPhiHE, VarManager::kMCPhiTildeHE};
int varspTCS[4] = {VarManager::kMCPt, VarManager::kMCCosThetaCS, VarManager::kMCPhiCS, VarManager::kMCPhiTildeCS};
int bins[4] = {20, 20, 20, 20};
double xmin[4] = {0., -1., 0., 0.};
double xmax[4] = {20., 1., 2. * o2::constants::math::PI, 2. * o2::constants::math::PI};
hm->AddHistogram(histClass, "Pt_cosThetaHE_phiHE_phiTildeHE", "", 4, varspTHE, bins, xmin, xmax, 0, -1, kFALSE);
hm->AddHistogram(histClass, "Pt_cosThetaCS_phiCS_phiTildeCS", "", 4, varspTCS, bins, xmin, xmax, 0, -1, kFALSE);
}
if (subGroupStr.Contains("pbpb")) {
hm->AddHistogram(histClass, "CosThetaStarRandom", "", false, 100, -1.0, 1.0, VarManager::kCosThetaStarRandom);
hm->AddHistogram(histClass, "Pt_Cent_CosThetaStarRandom", "", false, 20, 0.0, 20.0, VarManager::kMCPt, 9, 0., 90., VarManager::kMCEventCentrFT0C, 20, -1.0, 1.0, VarManager::kCosThetaStarRandom);
hm->AddHistogram(histClass, "Cos2ThetaStarRandom_Mass", "", true, 50, 2.0, 4.0, VarManager::kMCMass, 100, -1.0, 1.0, VarManager::kCos2ThetaStarRandom);
hm->AddHistogram(histClass, "CosThetaStar", "", false, 100, -1.0, 1.0, VarManager::kMCCosThetaStar);
}
}
}
if (!groupStr.CompareTo("mctruth_quad")) {
Expand Down Expand Up @@ -1348,12 +1361,21 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h
}
}
if (subGroupStr.Contains("globalpolarization")) {
hm->AddHistogram(histClass, "CosThetaStarTPC", "", false, 100, -1.0, 1.0, VarManager::kCosThetaStarTPC);
hm->AddHistogram(histClass, "CosThetaStarFT0A", "", false, 100, -1.0, 1.0, VarManager::kCosThetaStarFT0A);
hm->AddHistogram(histClass, "CosThetaStarFT0C", "", false, 100, -1.0, 1.0, VarManager::kCosThetaStarFT0C);
hm->AddHistogram(histClass, "CosThetaStarRandom", "", false, 100, -1.0, 1.0, VarManager::kCosThetaStarRandom);
hm->AddHistogram(histClass, "Cos2ThetaStarFT0C_Mass", "", true, 50, 2.0, 4.0, VarManager::kMass, 100, -1.0, 1.0, VarManager::kCos2ThetaStarFT0C);
hm->AddHistogram(histClass, "Cos2ThetaStarRandom_Mass", "", true, 50, 2.0, 4.0, VarManager::kMass, 100, -1.0, 1.0, VarManager::kCos2ThetaStarRandom);
if (subGroupStr.Contains("tpc")) {
hm->AddHistogram(histClass, "CosThetaStarTPC", "", false, 100, -1.0, 1.0, VarManager::kCosThetaStarTPC);
}
if (subGroupStr.Contains("ft0")) {
hm->AddHistogram(histClass, "CosThetaStarFT0A", "", false, 100, -1.0, 1.0, VarManager::kCosThetaStarFT0A);
hm->AddHistogram(histClass, "CosThetaStarFT0C", "", false, 100, -1.0, 1.0, VarManager::kCosThetaStarFT0C);
hm->AddHistogram(histClass, "Cos2ThetaStarFT0C_Mass", "", true, 50, 2.0, 4.0, VarManager::kMass, 100, -1.0, 1.0, VarManager::kCos2ThetaStarFT0C);
hm->AddHistogram(histClass, "Cos2ThetaStarFT0C_Mass_Weight", "", true, 50, 2.0, 4.0, VarManager::kMass, 100, -1.0, 1.0, VarManager::kCos2ThetaStarFT0C, 1, 0, 1, VarManager::kNothing, "", "", "", -1, VarManager::kPairWeight, false, false);
}
if (subGroupStr.Contains("mc")) {
hm->AddHistogram(histClass, "EventPlaneAngle", "", false, 100, -TMath::Pi(), TMath::Pi(), VarManager::kMCEventPlaneAngle);
hm->AddHistogram(histClass, "CosThetaStarMC", "", false, 100, -1.0, 1.0, VarManager::kMCCosThetaStar);
}
}
if (subGroupStr.Contains("upsilon")) {
hm->AddHistogram(histClass, "MassUpsilon_Pt", "", false, 500, 7.0, 12.0, VarManager::kMass, 400, 0.0, 40.0, VarManager::kPt);
Expand Down
60 changes: 60 additions & 0 deletions PWGDQ/Core/VarManager.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,8 @@ std::map<VarManager::CalibObjects, TObject*> VarManager::fgCalibs;
bool VarManager::fgRunTPCPostCalibration[4] = {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
TObject* VarManager::fgEfficiencyHist = nullptr; // histogram for efficiency

//__________________________________________________________________
VarManager::VarManager() : TObject()
Expand Down Expand Up @@ -380,6 +382,58 @@ double VarManager::ComputePIDcalibration(int species, double nSigmaValue)
}
}

//__________________________________________________________________
void VarManager::SetEfficiencyObject(int efficiencyType, TObject* obj)
{
// check the type of the efficiency object and set it accordingly
if (efficiencyType >= kNEfficiencyTypes || efficiencyType < 0) {
LOG(warning) << "SetEfficiencyObject: unknown efficiency type " << efficiencyType;
return;
}

// set the efficiency type
fgEfficiencyType = efficiencyType;
// set the efficiency object
fgEfficiencyHist = obj;
}

void VarManager::FillEfficiency(float* values)
{
// depending on the efficiency type, we use different types of efficiency histograms and different variables to get the efficiency value
if (!values) {
values = fgValues;
}

if (fgEfficiencyType == kNone) {
values[kPairEfficiency] = 1.0; // if no efficiency is to be applied, set the efficiency value to 1
values[kPairWeight] = 1.0; // set the weight to 1
} else if (fgEfficiencyType == kPairPtCentFT0cCosThetaStarFT0c) {
if (!fgEfficiencyHist) {
LOG(fatal) << "efficiency histogram not set";
return;
}
TH3F* efficiencyHist = reinterpret_cast<TH3F*>(fgEfficiencyHist);
// Get the bin indices for the efficiency histogram
int binPt = efficiencyHist->GetXaxis()->FindBin(values[kPairPt]);
binPt = (binPt == 0 ? 1 : binPt);
binPt = (binPt > efficiencyHist->GetXaxis()->GetNbins() ? efficiencyHist->GetXaxis()->GetNbins() : binPt);
int binCent = efficiencyHist->GetYaxis()->FindBin(values[kCentFT0C]);
binCent = (binCent == 0 ? 1 : binCent);
binCent = (binCent > efficiencyHist->GetYaxis()->GetNbins() ? efficiencyHist->GetYaxis()->GetNbins() : binCent);
int binCosThetaStarFT0c = efficiencyHist->GetZaxis()->FindBin(values[kCosThetaStarFT0C]);
binCosThetaStarFT0c = (binCosThetaStarFT0c == 0 ? 1 : binCosThetaStarFT0c);
binCosThetaStarFT0c = (binCosThetaStarFT0c > efficiencyHist->GetZaxis()->GetNbins() ? efficiencyHist->GetZaxis()->GetNbins() : binCosThetaStarFT0c);

// get the efficiency value from the histogram
values[kPairEfficiency] = efficiencyHist->GetBinContent(binPt, binCent, binCosThetaStarFT0c);
values[kPairWeight] = 1.0 / (values[kPairEfficiency] > 0 ? values[kPairEfficiency] : 1.0); // set the weight as the inverse of the efficiency, but avoid division by zero
} else {
LOG(warning) << "FillEfficiency: unknown efficiency type " << fgEfficiencyType << ", using default efficiency = 1";
values[kPairEfficiency] = 1;
values[kPairWeight] = 1;
}
}

//__________________________________________________________________
std::tuple<float, float, float, float, float> VarManager::BimodalityCoefficientUnbinned(const std::vector<float>& data)
{
Expand Down Expand Up @@ -672,6 +726,7 @@ void VarManager::SetDefaultVarNames()
fgVariableNames[kMCEventTime] = "MC event time";
fgVariableNames[kMCEventWeight] = "MC event weight";
fgVariableNames[kMCEventImpParam] = "MC impact parameter";
fgVariableNames[kMCEventPlaneAngle] = "MC event plane angle";
fgVariableNames[kMCEventCentrFT0C] = "MC Centrality FT0C";
fgVariableNames[kMultMCNParticlesEta05] = "MC Multiplicity Central Barrel for |eta| < 0.5";
fgVariableNames[kMultMCNParticlesEta08] = "MC Multiplicity Central Barrel for |eta| < 0.8";
Expand All @@ -684,6 +739,7 @@ void VarManager::SetDefaultVarNames()
fgVariableUnits[kMCEventTime] = ""; // TODO: add proper unit
fgVariableUnits[kMCEventWeight] = "";
fgVariableUnits[kMCEventImpParam] = "b";
fgVariableUnits[kMCEventPlaneAngle] = "";
fgVariableUnits[kMCEventCentrFT0C] = "%";
fgVariableUnits[kMultMCNParticlesEta05] = "Multiplicity_eta05";
fgVariableUnits[kMultMCNParticlesEta08] = "Multiplicity_eta08";
Expand Down Expand Up @@ -1531,6 +1587,8 @@ void VarManager::SetDefaultVarNames()
fgVariableUnits[kCosThetaStarRandom] = "";
fgVariableNames[kCos2ThetaStarRandom] = "cos^{2}#it{#theta}^{*}_{Random}";
fgVariableUnits[kCos2ThetaStarRandom] = "";
fgVariableNames[kMCCosThetaStar] = "cos#it{#theta}^{*}_{MC}";
fgVariableUnits[kMCCosThetaStar] = "";
fgVariableNames[kCosPhiVP] = "cos#it{#varphi}_{VP}";
fgVariableUnits[kCosPhiVP] = "";
fgVariableNames[kPhiVP] = "#varphi_{VP} - #Psi_{2}";
Expand Down Expand Up @@ -1948,6 +2006,7 @@ void VarManager::SetDefaultVarNames()
fgVarNamesMap["kMCEventTime"] = kMCEventTime;
fgVarNamesMap["kMCEventWeight"] = kMCEventWeight;
fgVarNamesMap["kMCEventImpParam"] = kMCEventImpParam;
fgVarNamesMap["kMCEventPlaneAngle"] = kMCEventPlaneAngle;
fgVarNamesMap["kQ1ZNAX"] = kQ1ZNAX;
fgVarNamesMap["kQ1ZNAY"] = kQ1ZNAY;
fgVarNamesMap["kQ1ZNCX"] = kQ1ZNCX;
Expand Down Expand Up @@ -2357,6 +2416,7 @@ void VarManager::SetDefaultVarNames()
fgVarNamesMap["kCos2ThetaStarFT0C"] = kCos2ThetaStarFT0C;
fgVarNamesMap["kCosThetaStarRandom"] = kCosThetaStarRandom;
fgVarNamesMap["kCos2ThetaStarRandom"] = kCos2ThetaStarRandom;
fgVarNamesMap["kMCCosThetaStar"] = kMCCosThetaStar;
fgVarNamesMap["kCosPhiVP"] = kCosPhiVP;
fgVarNamesMap["kPhiVP"] = kPhiVP;
fgVarNamesMap["kDeltaPhiPair2"] = kDeltaPhiPair2;
Expand Down
72 changes: 64 additions & 8 deletions PWGDQ/Core/VarManager.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,8 @@
#include <KFPTrack.h>
#include <KFPVertex.h>
#include <KFParticle.h>
#include <rapidjson/document.h>
#include <rapidjson/error/error.h>

#include <GPUROOTCartesianFwd.h>
#include <Rtypes.h>
Expand Down Expand Up @@ -338,6 +340,7 @@ class VarManager : public TObject
kMCEventWeight,
kMCEventImpParam,
kMCEventCentrFT0C,
kMCEventPlaneAngle,
kMultMCNParticlesEta10,
kMultMCNParticlesEta08,
kMultMCNParticlesEta05,
Expand Down Expand Up @@ -726,6 +729,7 @@ class VarManager : public TObject
kMCPhiPP,
kMCPhiTildePP,
kMCCosThetaRM,
kMCCosThetaStar,

// Pair variables
kCandidateId,
Expand Down Expand Up @@ -924,6 +928,8 @@ class VarManager : public TObject
kS12,
kS13,
kS23,
kPairEfficiency,
kPairWeight,
kNPairVariables,

// Candidate-track correlation variables
Expand Down Expand Up @@ -1109,6 +1115,13 @@ class VarManager : public TObject
kNCalibObjects
};

enum EfficiencyType {
kNone = 0,
kPairPtCentFT0cCosThetaStarFT0c,
// Add more efficiency types as needed
kNEfficiencyTypes
};

enum DileptonCharmHadronTypes {
kJPsi = 0,
kD0ToPiK,
Expand Down Expand Up @@ -1459,7 +1472,8 @@ class VarManager : public TObject
fgUseInterpolatedCalibration = useInterpolation;
}
static double ComputePIDcalibration(int species, double nSigmaValue);

static void SetEfficiencyObject(int type, TObject* obj);
static void FillEfficiency(float* values = nullptr);
static TObject* GetCalibrationObject(CalibObjects calib)
{
auto obj = fgCalibs.find(calib);
Expand Down Expand Up @@ -1517,6 +1531,7 @@ class VarManager : public TObject

// static void FillEventDerived(float* values = nullptr);
static void FillTrackDerived(float* values = nullptr);
// static void FillPairDerived(float* values = nullptr);
template <typename T, typename U, typename V>
static auto getRotatedCovMatrixXX(const T& matrix, U phi, V theta);
template <typename T>
Expand All @@ -1543,6 +1558,9 @@ class VarManager : public TObject
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)

static int fgEfficiencyType; // type of efficiency correction to apply
static TObject* fgEfficiencyHist; // histogram for efficiency correction

VarManager& operator=(const VarManager& c);
VarManager(const VarManager& c);

Expand Down Expand Up @@ -1904,6 +1922,10 @@ void VarManager::FillEvent(T const& event, float* values)
values[kCollisionRandom] = gRandom->Rndm();
}

if (fgUsedVars[kRandomPsi2]) {
values[kRandomPsi2] = gRandom->Uniform(-o2::constants::math::PIHalf, o2::constants::math::PIHalf);
}

if constexpr ((fillMap & Collision) > 0) {
// TODO: trigger info from the event selection requires a separate flag
// so that it can be switched off independently of the rest of Collision variables (e.g. if event selection is not available)
Expand Down Expand Up @@ -2210,7 +2232,6 @@ void VarManager::FillEvent(T const& event, float* values)
values[VarManager::kPsi2A] = Psi2A;
values[VarManager::kPsi2B] = Psi2B;
values[VarManager::kPsi2C] = Psi2C;
values[VarManager::kRandomPsi2] = gRandom->Uniform(-o2::constants::math::PIHalf, o2::constants::math::PIHalf);

if constexpr ((fillMap & ReducedEventQvectorExtra) > 0) {
values[kQ42XA] = event.q42xa();
Expand Down Expand Up @@ -2332,7 +2353,8 @@ void VarManager::FillEvent(T const& event, float* values)
values[kMCEventTime] = event.t();
values[kMCEventWeight] = event.weight();
values[kMCEventImpParam] = event.impactParameter();
if constexpr ((fillMap & CollisionCent) > 0) {
values[kMCEventPlaneAngle] = event.eventPlaneAngle();
if constexpr (requires { event.bestCollisionCentFT0C(); }) {
// WARNING: temporary solution, ongoing work to provide proper MC gen. centrality
values[kMCEventCentrFT0C] = event.bestCollisionCentFT0C();
values[kMultMCNParticlesEta05] = event.multMCNParticlesEta05();
Expand Down Expand Up @@ -3661,6 +3683,27 @@ void VarManager::FillPair(T1 const& t1, T2 const& t2, float* values)
}
}

if (fgUsedVars[kCosThetaStarRandom]) {
ROOT::Math::Boost boostv12{v12.BoostToCM()};
ROOT::Math::XYZVectorF v1_CM{(boostv12(v1).Vect()).Unit()};
ROOT::Math::XYZVectorF v2_CM{(boostv12(v2).Vect()).Unit()};

// using positive sign convention for the first track
ROOT::Math::XYZVectorF v_CM = (t1.sign() > 0 ? v1_CM : v2_CM);

// Randomize the event plane angle to check the unpolarized contribution
ROOT::Math::XYZVector zaxisRandom = ROOT::Math::XYZVector(TMath::Cos(values[kRandomPsi2]), TMath::Sin(values[kRandomPsi2]), 0).Unit();
values[kCosThetaStarRandom] = v_CM.Dot(zaxisRandom);
values[kCos2ThetaStarRandom] = values[kCosThetaStarRandom] * values[kCosThetaStarRandom];

// if the truth event plane angle is available, calculate the cos(theta*) with respect to the true event plane angle for comparison
if (fgUsedVars[kMCCosThetaStar] && fgUsedVars[kMCEventPlaneAngle]) {
// truth event plane angle
ROOT::Math::XYZVector zaxisTrue = ROOT::Math::XYZVector(TMath::Cos(values[kMCEventPlaneAngle]), TMath::Sin(values[kMCEventPlaneAngle]), 0).Unit();
values[kMCCosThetaStar] = v_CM.Dot(zaxisTrue);
}
}

if constexpr ((pairType == kDecayToEE) && ((fillMap & TrackCov) > 0 || (fillMap & ReducedTrackBarrelCov) > 0)) {

if (fgUsedVars[kQuadDCAabsXY] || fgUsedVars[kQuadDCAsigXY] || fgUsedVars[kQuadDCAabsZ] || fgUsedVars[kQuadDCAsigZ] || fgUsedVars[kQuadDCAsigXYZ] || fgUsedVars[kSignQuadDCAsigXY]) {
Expand Down Expand Up @@ -4293,6 +4336,24 @@ void VarManager::FillPairMC(T1 const& t1, T2 const& t2, float* values)
values[kMCCosThetaRM] = zaxis_RM.Dot(v_CM);
}
}

if (fgUsedVars[kCosThetaStarRandom] || fgUsedVars[kMCCosThetaStar]) {
ROOT::Math::Boost boostv12{v12.BoostToCM()};
ROOT::Math::XYZVectorF v1_CM{(boostv12(v1).Vect()).Unit()};
ROOT::Math::XYZVectorF v2_CM{(boostv12(v2).Vect()).Unit()};

// using positive sign convention for the first track
ROOT::Math::XYZVectorF v_CM = (t1.pdgCode() > 0 ? v1_CM : v2_CM);

// Randomize the event plane angle to check the unpolarized contribution
ROOT::Math::XYZVector zaxisRandom = ROOT::Math::XYZVector(TMath::Cos(values[kRandomPsi2]), TMath::Sin(values[kRandomPsi2]), 0).Unit();
values[kCosThetaStarRandom] = v_CM.Dot(zaxisRandom);
values[kCos2ThetaStarRandom] = values[kCosThetaStarRandom] * values[kCosThetaStarRandom];

// truth event plane angle
ROOT::Math::XYZVector zaxisTrue = ROOT::Math::XYZVector(TMath::Cos(values[kMCEventPlaneAngle]), TMath::Sin(values[kMCEventPlaneAngle]), 0).Unit();
values[kMCCosThetaStar] = v_CM.Dot(zaxisTrue);
}
}

template <int candidateType, typename T1, typename T2, typename T3>
Expand Down Expand Up @@ -5721,11 +5782,6 @@ void VarManager::FillPairVn(T1 const& t1, T2 const& t2, float* values)
values[kCosThetaStarFT0C] = v_CM.Dot(zaxisFT0C);
values[kAbsCosThetaStarFT0C] = std::abs(values[kCosThetaStarFT0C]);
values[kCos2ThetaStarFT0C] = values[kCosThetaStarFT0C] * values[kCosThetaStarFT0C];

// Randomize the event plane angle to check the unpolarized contribution
ROOT::Math::XYZVector zaxisRandom = ROOT::Math::XYZVector(TMath::Cos(values[kRandomPsi2]), TMath::Sin(values[kRandomPsi2]), 0).Unit();
values[kCosThetaStarRandom] = v_CM.Dot(zaxisRandom);
values[kCos2ThetaStarRandom] = values[kCosThetaStarRandom] * values[kCosThetaStarRandom];
}

// kV4, kC4POI, kC4REF etc.
Expand Down
Loading
Loading