Skip to content

[PWGDQ] Add and change the code for the energy correlator analysis #14325

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

Open
YazhenLin wants to merge 9 commits into
base: master
Choose a base branch
from
+76 −33
Open

[PWGDQ] Add and change the code for the energy correlator analysis #14325

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
f7aef3e
Add files via upload
YazhenLin Dec 17, 2025
8f59022
Add files via upload
YazhenLin Dec 17, 2025
2646b10
Add files via upload
YazhenLin Dec 17, 2025
fc5de50
Add files via upload
YazhenLin Dec 17, 2025
b05a247
Add files via upload
YazhenLin Dec 17, 2025
ab05454
Merge branch 'AliceO2Group:master' into master
YazhenLin Dec 22, 2025
b12e9b5
Add files via upload
YazhenLin Dec 22, 2025
f62d934
Add files via upload
YazhenLin Dec 22, 2025
73b8812
Add files via upload
YazhenLin Dec 22, 2025
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
33 changes: 23 additions & 10 deletions PWGDQ/Core/HistogramsLibrary.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1942,23 +1942,36 @@
hm->AddHistogram(histClass, "Delta_Mass_DstarD0region", "", false, 50, 0.14, 0.16, VarManager::kDeltaMass);
}
if (subGroupStr.Contains("energy-correlator")) {
hm->AddHistogram(histClass, "Coschi", "", false, 40, -1.0, 1.0, VarManager::kCosChi, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kECWeight);
hm->AddHistogram(histClass, "CosTheta_woweight", "", false, 40, -1.0, 1.0, VarManager::kCosTheta);
double coschiBins[26];
for (int i = 0; i < 26; i++) {
coschiBins[i] = -1.0 + 2.0 * TMath::Power(0.04 * i, 2.0);
}

double deltaetaBins[21];
for (int i = 0; i < 21; i++) {
deltaetaBins[i] = -2.0 + 0.2 * i;
}

hm->AddHistogram(histClass, "Coschi", "", false, 25, coschiBins, VarManager::kCosChi, 0, nullptr, -1, 0, nullptr, -1, "", "", "", -1, VarManager::kECWeight);
hm->AddHistogram(histClass, "CosTheta", "", false, 40, -1.0, 1.0, VarManager::kCosTheta, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kEWeight_before);
hm->AddHistogram(histClass, "Pt_Hadron", ";P_{T}", false, 120, 0.0, 30.0, VarManager::kPtDau);
hm->AddHistogram(histClass, "Coschi_wo", "", false, 25, coschiBins, VarManager::kCosChi);
hm->AddHistogram(histClass, "Eta_Hadron", ";#eta", false, 120, -2.0, 2.0, VarManager::kEtaDau);
hm->AddHistogram(histClass, "Phi_Hadron", ";#phi", false, 120, -8, 8, VarManager::kPhiDau);
hm->AddHistogram(histClass, "DeltaEta_DeltaPhi", "", false, 20, -2.0, 2.0, VarManager::kDeltaEta, 50, -8.0, 8.0, VarManager::kDeltaPhi);
hm->AddHistogram(histClass, "DeltaEta_DeltaPhi_weight", "", false, 20, -2.0, 2.0, VarManager::kDeltaEta, 50, -2.0, 6.0, VarManager::kDeltaPhi, 0, 0, 0, -1, "", "", "", -1, VarManager::kEWeight_before);
hm->AddHistogram(histClass, "DeltaEta", "", false, 20, -2.0, 2.0, VarManager::kDeltaEta);
hm->AddHistogram(histClass, "DeltaPhi", "", false, 50, -8.0, 8.0, VarManager::kDeltaPhi);
hm->AddHistogram(histClass, "DeltaPhi", "", false, 50, -2.0, 6.0, VarManager::kDeltaPhi);
hm->AddHistogram(histClass, "DeltaPhi_weight", "", false, 50, -2.0, 6.0, VarManager::kDeltaPhi, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kEWeight_before);
hm->AddHistogram(histClass, "Coschi_DeltaEta", "", false, 25, coschiBins, VarManager::kCosChi, 20, deltaetaBins, VarManager::kDeltaEta, 0, nullptr, -1, "", "", "", -1, VarManager::kECWeight);
hm->AddHistogram(histClass, "Coschi_wo_DeltaEta", "", false, 25, coschiBins, VarManager::kCosChi, 20, deltaetaBins, VarManager::kDeltaEta);
// for bkg
hm->AddHistogram(histClass, "DeltaPhi_randomPhi_trans", "", false, 50, -8.0, 8.0, VarManager::kdeltaphi_randomPhi_trans);
hm->AddHistogram(histClass, "DeltaPhi_randomPhi_toward", "", false, 50, -8.0, 8.0, VarManager::kdeltaphi_randomPhi_toward);
hm->AddHistogram(histClass, "DeltaPhi_randomPhi_away", "", false, 50, -8.0, 8.0, VarManager::kdeltaphi_randomPhi_away);
hm->AddHistogram(histClass, "Coschi_randomPhi_trans", "", false, 25, coschiBins, VarManager::kCosChi_randomPhi_trans, 0, nullptr, -1, 0, nullptr, -1, "", "", "", -1, VarManager::kWeight_randomPhi_trans);
hm->AddHistogram(histClass, "Coschi_randomPhi_toward", "", false, 25, coschiBins, VarManager::kCosChi_randomPhi_toward, 0, nullptr, -1, 0, nullptr, -1, "", "", "", -1, VarManager::kWeight_randomPhi_toward);
hm->AddHistogram(histClass, "Coschi_randomPhi_away", "", false, 25, coschiBins, VarManager::kCosChi_randomPhi_away, 0, nullptr, -1, 0, nullptr, -1, "", "", "", -1, VarManager::kWeight_randomPhi_away);

hm->AddHistogram(histClass, "Coschi_randomPhi_trans", "", false, 40, -1.0, 1.0, VarManager::kCosChi_randomPhi_trans, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kWeight_randomPhi_trans);
hm->AddHistogram(histClass, "Coschi_randomPhi_toward", "", false, 40, -1.0, 1.0, VarManager::kCosChi_randomPhi_toward, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kWeight_randomPhi_toward);
hm->AddHistogram(histClass, "Coschi_randomPhi_away", "", false, 40, -1.0, 1.0, VarManager::kCosChi_randomPhi_away, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kWeight_randomPhi_away);
hm->AddHistogram(histClass, "Coschi_wo_randomPhi_trans", "", false, 25, coschiBins, VarManager::kCosChi_randomPhi_trans);
hm->AddHistogram(histClass, "Coschi_wo_randomPhi_toward", "", false, 25, coschiBins, VarManager::kCosChi_randomPhi_toward);
hm->AddHistogram(histClass, "Coschi_wo_randomPhi_away", "", false, 25, coschiBins, VarManager::kCosChi_randomPhi_away);
}
}

Expand Down Expand Up @@ -2131,7 +2144,7 @@
LOG(fatal) << "histClass field should be an array of strings, e.g. [class1, class2]";
return false;
}
for (auto& v : hist->FindMember("histClass")->value.GetArray()) {

Check failure on line 2147 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.

Check failure on line 2147 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
if (!v.IsString()) {
LOG(fatal) << "histClass field should be an array of strings, e.g. [class1, class2]";
return false;
Expand Down Expand Up @@ -2287,7 +2300,7 @@
}
}
if (isTHn) {
for (auto& v : hist->FindMember("vars")->value.GetArray()) {

Check failure on line 2303 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.

Check failure on line 2303 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
if (VarManager::fgVarNamesMap.find(v.GetString()) == VarManager::fgVarNamesMap.end()) {
LOG(fatal) << "Bad variable in vars (" << v.GetString() << ") specified for histogram";
return false;
Expand Down Expand Up @@ -2350,7 +2363,7 @@

// create an array of strings to store the different histogram classes
std::vector<const char*> histClasses;
for (auto& v : hist.FindMember("histClass")->value.GetArray()) {

Check failure on line 2366 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.

Check failure on line 2366 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
histClasses.push_back(v.GetString());
}
const char* title = hist.FindMember("title")->value.GetString();
Expand All @@ -2361,7 +2374,7 @@

int* vars = new int[nDimensions];
int iDim = 0;
for (auto& v : hist.FindMember("vars")->value.GetArray()) {

Check failure on line 2377 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.

Check failure on line 2377 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
LOG(debug) << "iDim " << iDim << ": " << v.GetString();
vars[iDim++] = VarManager::fgVarNamesMap[v.GetString()];
}
Expand All @@ -2375,27 +2388,27 @@
xmin = new double[nDimensions];
xmax = new double[nDimensions];
int iDim = 0;
for (auto& v : hist.FindMember("nBins")->value.GetArray()) {

Check failure on line 2391 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.

Check failure on line 2391 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
nBins[iDim++] = v.GetInt();
LOG(debug) << "nBins " << iDim << ": " << nBins[iDim - 1];
}
iDim = 0;
for (auto& v : hist.FindMember("xmin")->value.GetArray()) {

Check failure on line 2396 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.

Check failure on line 2396 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
xmin[iDim++] = v.GetDouble();
LOG(debug) << "xmin " << iDim << ": " << xmin[iDim - 1];
}
iDim = 0;
for (auto& v : hist.FindMember("xmax")->value.GetArray()) {

Check failure on line 2401 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.

Check failure on line 2401 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
xmax[iDim++] = v.GetDouble();
LOG(debug) << "xmax " << iDim << ": " << xmax[iDim - 1];
}
} else {
int iDim = 0;
binLimits = new TArrayD[nDimensions];
for (auto& v : hist.FindMember("binLimits")->value.GetArray()) {

Check failure on line 2408 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.

Check failure on line 2408 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
double* lims = new double[v.GetArray().Size()];
int iElem = 0;
for (auto& lim : v.GetArray()) {

Check failure on line 2411 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.

Check failure on line 2411 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
lims[iElem++] = lim.GetDouble();
}
binLimits[iDim++] = TArrayD(v.GetArray().Size(), lims);
Expand All @@ -2406,7 +2419,7 @@
if (hist.HasMember("axLabels")) {
axLabels = new TString[hist.FindMember("axLabels")->value.GetArray().Size()];
int iDim = 0;
for (auto& v : hist.FindMember("axLabels")->value.GetArray()) {

Check failure on line 2422 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.

Check failure on line 2422 in PWGDQ/Core/HistogramsLibrary.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
axLabels[iDim++] = v.GetString();
}
}
Expand Down
57 changes: 34 additions & 23 deletions PWGDQ/Core/VarManager.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1231,6 +1231,8 @@ class VarManager : public TObject
template <typename T1, typename T2>
static void FillDileptonHadron(T1 const& dilepton, T2 const& hadron, float* values = nullptr, float hadronMass = 0.0f);
template <typename T1, typename T2>
static void FillEnergyCorrelator(T1 const& dilepton, T2 const& hadron, float* values = nullptr, bool applyFitMass = false, float sidebandMass = 0.0f);
template <typename T1, typename T2>
static void FillDileptonPhoton(T1 const& dilepton, T2 const& photon, float* values = nullptr);
template <typename T>
static void FillHadron(T const& hadron, float* values = nullptr, float hadronMass = 0.0f);
Expand Down Expand Up @@ -5301,8 +5303,39 @@ void VarManager::FillDileptonHadron(T1 const& dilepton, T2 const& hadron, float*
double Q1 = (dilepton.mass() * dilepton.mass() - hadronMass * hadronMass) / Pinv;
values[kDileptonHadronKstar] = sqrt(Q1 * Q1 - v12_Qvect.M2()) / 2.0;
}

if (fgUsedVars[kDeltaPhi]) {
double delta = dilepton.phi() - hadron.phi();
if (delta > 3.0 / 2.0 * M_PI) {
delta -= 2.0 * M_PI;
}
if (delta < -0.5 * M_PI) {
delta += 2.0 * M_PI;
}
values[kDeltaPhi] = delta;
}
if (fgUsedVars[kDeltaPhiSym]) {
double delta = std::abs(dilepton.phi() - hadron.phi());
if (delta > M_PI) {
delta = 2 * M_PI - delta;
}
values[kDeltaPhiSym] = delta;
}
if (fgUsedVars[kDeltaEta]) {
values[kDeltaEta] = dilepton.eta() - hadron.eta();
}
}

template <typename T1, typename T2>
void VarManager::FillEnergyCorrelator(T1 const& dilepton, T2 const& hadron, float* values, bool applyFitMass, float sidebandMass)
{
float dileptonmass = dilepton.mass();
if (applyFitMass) {
dileptonmass = sidebandMass;
}

if (fgUsedVars[kCosChi] || fgUsedVars[kECWeight] || fgUsedVars[kCosTheta] || fgUsedVars[kEWeight_before] || fgUsedVars[kPtDau] || fgUsedVars[kEtaDau] || fgUsedVars[kPhiDau] || fgUsedVars[kCosChi_randomPhi_trans] || fgUsedVars[kCosChi_randomPhi_toward] || fgUsedVars[kCosChi_randomPhi_away]) {
ROOT::Math::PtEtaPhiMVector v1(dilepton.pt(), dilepton.eta(), dilepton.phi(), dilepton.mass());
ROOT::Math::PtEtaPhiMVector v1(dilepton.pt(), dilepton.eta(), dilepton.phi(), dileptonmass);
ROOT::Math::PtEtaPhiMVector v2(hadron.pt(), hadron.eta(), hadron.phi(), o2::constants::physics::MassPionCharged);
values[kCosChi] = LorentzTransformJpsihadroncosChi("coschi", v1, v2);
float E_boost = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2);
Expand Down Expand Up @@ -5348,29 +5381,7 @@ void VarManager::FillDileptonHadron(T1 const& dilepton, T2 const& hadron, float*
values[kdeltaphi_randomPhi_away] = RecoDecay::constrainAngle(v1.phi() - randomPhi_away, -o2::constants::math::PIHalf);
}
}

if (fgUsedVars[kDeltaPhi]) {
double delta = dilepton.phi() - hadron.phi();
if (delta > 3.0 / 2.0 * M_PI) {
delta -= 2.0 * M_PI;
}
if (delta < -0.5 * M_PI) {
delta += 2.0 * M_PI;
}
values[kDeltaPhi] = delta;
}
if (fgUsedVars[kDeltaPhiSym]) {
double delta = std::abs(dilepton.phi() - hadron.phi());
if (delta > M_PI) {
delta = 2 * M_PI - delta;
}
values[kDeltaPhiSym] = delta;
}
if (fgUsedVars[kDeltaEta]) {
values[kDeltaEta] = dilepton.eta() - hadron.eta();
}
}

template <typename T1, typename T2>
void VarManager::FillDileptonPhoton(T1 const& dilepton, T2 const& photon, float* values)
{
Expand Down
19 changes: 19 additions & 0 deletions PWGDQ/Tasks/tableReader_withAssoc.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,11 +42,13 @@
#include "ITSMFTBase/DPLAlpideParam.h"

#include "TGeoGlobalMagField.h"
#include <TF1.h>
#include <TH1F.h>
#include <TH3F.h>
#include <THashList.h>
#include <TList.h>
#include <TObjString.h>
#include <TRandom.h>
#include <TString.h>

#include <algorithm>
Expand Down Expand Up @@ -3133,6 +3135,8 @@ struct AnalysisDileptonTrack {
Configurable<std::string> fConfigGeoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"};
Configurable<bool> fConfigUseRapcut{"cfgUseMCRapcut", false, "Use Rap cut for dileptons used in the triplet vertexing"};
Configurable<bool> fConfigEnergycorrelator{"cfgEnergycorrelator", false, "Add some hist for energy correlator study"};
Configurable<bool> fConfigApplyMassEC{"cfgApplyMassEC", false, "Apply fit mass for sideband for the energy correlator study"};
Configurable<std::vector<float>> fConfigFitmassEC{"cfgTFitmassEC", std::vector<float>{61199.9144, -19417.7588, 1709.4958, 2.8, 3.2}, "parameter from the fit fuction and fit range"};

int fCurrentRun; // needed to detect if the run changed and trigger update of calibrations etc.
int fNCuts; // number of dilepton leg cuts
Expand Down Expand Up @@ -3165,6 +3169,8 @@ struct AnalysisDileptonTrack {

NoBinningPolicy<aod::dqanalysisflags::MixingHash> fHashBin;

TF1* fMassBkg = nullptr;

void init(o2::framework::InitContext& context)
{
bool isBarrel = context.mOptions.get<bool>("processBarrelSkimmed");
Expand Down Expand Up @@ -3389,6 +3395,12 @@ struct AnalysisDileptonTrack {
LOG(info) << "Loading geometry from CCDB in dilepton-track task";
fCCDB->get<TGeoManager>(fConfigGeoPath);
}

// the background mass distribution in signal region
std::vector<float> fMassBkgpar = fConfigFitmassEC;
fMassBkg = new TF1("fMassBkg", "[0]+[1]*x+[2]*x*x", fMassBkgpar[3], fMassBkgpar[4]);
fMassBkg->SetParameters(fMassBkgpar[0], fMassBkgpar[1], fMassBkgpar[2]);
fMassBkg->SetNpx(10000);
}

// init parameters from CCDB
Expand Down Expand Up @@ -3490,6 +3502,10 @@ struct AnalysisDileptonTrack {
// compute needed quantities
VarManager::FillDileptonHadron(dilepton, track, fValuesHadron);
VarManager::FillDileptonTrackVertexing<TCandidateType, TEventFillMap, TTrackFillMap>(event, lepton1, lepton2, track, fValuesHadron);

// for the energy correlator analysis
VarManager::FillEnergyCorrelator(dilepton, track, fValuesHadron, fConfigApplyMassEC, fMassBkg->GetRandom());

// table to be written out for ML analysis
BmesonsTable(event.runNumber(), event.globalIndex(), event.timestamp(), fValuesHadron[VarManager::kPairMass], dilepton.mass(), fValuesHadron[VarManager::kDeltaMass], fValuesHadron[VarManager::kPairPt], fValuesHadron[VarManager::kPairEta], fValuesHadron[VarManager::kPairPhi], fValuesHadron[VarManager::kPairRap],
fValuesHadron[VarManager::kVertexingLxy], fValuesHadron[VarManager::kVertexingLxyz], fValuesHadron[VarManager::kVertexingLz],
Expand Down Expand Up @@ -3685,6 +3701,9 @@ struct AnalysisDileptonTrack {
// compute dilepton - track quantities
VarManager::FillDileptonHadron(dilepton, track, VarManager::fgValues);

// for the energy correlator analysis
VarManager::FillEnergyCorrelator(dilepton, track, VarManager::fgValues, fConfigApplyMassEC, fMassBkg->GetRandom());

// loop over dilepton leg cuts and track cuts and fill histograms separately for each combination
for (int icut = 0; icut < fNCuts; icut++) {
if (!dilepton.filterMap_bit(icut)) {
Expand Down
Loading