Skip to content

[PWGJE] fix a few errors. introduced dcaxy cut. added event counter for subsamples #17974

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

Merged
nzardosh merged 1 commit into from
Sep 21, 2026
+52 −61
Merged

[PWGJE] fix a few errors. introduced dcaxy cut. added event counter for subsamples #17974

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
113 changes: 52 additions & 61 deletions PWGJE/Tasks/jetSpectraEseTask.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,8 +70,11 @@ using namespace o2::framework::expressions;

struct JetSpectraEseTask {
Configurable<int> cfgSystFlag{"cfgSystFlag", -1, "-1 = use configurable cuts, 0 = nominal defaults, >0 = systematic cut variation"};
Configurable<std::string> cfgEfficiency{"cfgEfficiency", "", "CCDB path to efficiency"};
Configurable<std::string> cfgEfficiency3D{"cfgEfficiency3D", "", "CCDB path to 3D efficiency"};

struct : ConfigurableGroup {
Configurable<std::string> cfgEfficiency{"cfgEfficiency", "", "CCDB path to efficiency"};
Configurable<std::string> cfgEfficiency3D{"cfgEfficiency3D", "", "CCDB path to 3D efficiency"};
} efficiencyFiles;
Configurable<float> jetPtMin{"jetPtMin", 5.0, "minimum jet pT cut"};
Configurable<float> jetR{"jetR", 0.2, "jet resolution parameter"};
Configurable<float> randomConeR{"randomConeR", 0.4, "size of random Cone for estimating background fluctuations"};
Expand Down Expand Up @@ -101,8 +104,9 @@ struct JetSpectraEseTask {
Configurable<float> trackDCAzMax{"trackDCAzMax", 2.0, "maximum DCAz for tracks"};
Configurable<float> cfgNTPCXrows{"cfgNTPCXrows", 70, "minimum number of crossed rows in the TPC for tracks"};
Configurable<float> cfgNTPCCls{"cfgNTPCCls", 50, "minimum number of clusters in the TPC for tracks"};
Configurable<float> cfgChi2PrTPCcls{"cfgChi2PrTPCcls", 2.5, "cut for chi2 per TPC cluster for tracks"};
Configurable<float> cfgChi2PrTPCcls{"cfgChi2PrTPCcls", 4.0, "cut for chi2 per TPC cluster for tracks"};
Configurable<float> cfgChi2PrITScls{"cfgChi2PrITScls", 36, "cut for chi2 per ITS cluster for tracks"};
Configurable<float> dcaXYSigmaMax{"dcaXYSigmaMax", 7.0, "maximum DCAxy sigma significance for tracks"};
} trackCuts;

Configurable<std::vector<float>> trackPtRhoPhi{"trackPtRhoPhi", {0.2, 5.0}, "pT range for tracks used in rho(phi) calculation"};
Expand Down Expand Up @@ -162,22 +166,15 @@ struct JetSpectraEseTask {
static constexpr float RScale = 100.0f;
static constexpr float EtaAcceptance = 0.9f;
static constexpr float LowFT0Cut = 1e-8;
static constexpr float NominalVertexZCut = 10.0f;
static constexpr float NominalTrackDCAzMax = 2.0f;
static constexpr float NominalNTPCXrows = 70.0f;
static constexpr float NominalNTPCCls = 50.0f;
static constexpr float NominalChi2PrTPCcls = 2.5f;
static constexpr float NominalChi2PrITScls = 36.0f;
static constexpr float NoDCAxySigmaCut = -1.0f;
static constexpr float LooseVertexZFilter = 999.0f;

struct SystematicCuts {
float vertexZCut = 10.0f;
float trackDCAzMax = 2.0f;
float dcaXYSigmaMax = -1.0f;
float dcaXYSigmaMax = 7.0f;
float nTPCXrows = 70.0f;
float nTPCCls = 50.0f;
float chi2PrTPCcls = 2.5f;
float chi2PrTPCcls = 4.0f;
float chi2PrITScls = 36.0f;
} systCuts;

Expand Down Expand Up @@ -252,12 +249,13 @@ struct JetSpectraEseTask {
static constexpr int NumSubSmpl = 5;
static constexpr int NumSavedRhoFitEvents = 5;
std::array<std::shared_ptr<THnSparse>, NumSubSmpl> hSameSub;
std::array<std::shared_ptr<TH1>, NumSubSmpl> hCentralityAnalyzedSub;

void applySystematicPreset()
{
systCuts.vertexZCut = vertexZCut.value;
systCuts.trackDCAzMax = trackCuts.trackDCAzMax.value;
systCuts.dcaXYSigmaMax = NoDCAxySigmaCut;
systCuts.dcaXYSigmaMax = trackCuts.dcaXYSigmaMax.value;
systCuts.nTPCXrows = trackCuts.cfgNTPCXrows.value;
systCuts.nTPCCls = trackCuts.cfgNTPCCls.value;
systCuts.chi2PrTPCcls = trackCuts.cfgChi2PrTPCcls.value;
Expand All @@ -268,15 +266,6 @@ struct JetSpectraEseTask {
LOGF(info, "JetSpectraEseTask::applySystematicPreset() - using configurable cut values");
return;
}

systCuts.vertexZCut = NominalVertexZCut;
systCuts.trackDCAzMax = NominalTrackDCAzMax;
systCuts.dcaXYSigmaMax = NoDCAxySigmaCut;
systCuts.nTPCXrows = NominalNTPCXrows;
systCuts.nTPCCls = NominalNTPCCls;
systCuts.chi2PrTPCcls = NominalChi2PrTPCcls;
systCuts.chi2PrITScls = NominalChi2PrITScls;

switch (flag) {
case 0:
break;
Expand Down Expand Up @@ -305,7 +294,7 @@ struct JetSpectraEseTask {
systCuts.nTPCCls = 70.0f;
break;
case 9:
systCuts.chi2PrTPCcls = 4.0f;
systCuts.chi2PrTPCcls = 2.5f;
break;
case 10:
systCuts.nTPCXrows = 80.0f;
Expand Down Expand Up @@ -435,6 +424,8 @@ struct JetSpectraEseTask {
for (int i = 0; i < NumSubSmpl; ++i) {
std::string n = fmt::format("subsamples/thn_jethad_corr_same_subsample{}", i);
hSameSub[i] = registry.add<THnSparse>(n, "same;...", o2::framework::HistType::kTHnSparseF, axes);
std::string n2 = fmt::format("subsamples/hCentralityAnalyzed_subsample{}", i);
hCentralityAnalyzedSub[i] = registry.add<TH1>(n2, ";Centrality;entries", o2::framework::HistType::kTH1F, {{centAxis}});
}
registry.add("hNtrig", "", {HistType::kTHnSparseF, {{centAxis}, {jetPtAxis}, {dPhiAxis}, {eseAxis}}});

Expand All @@ -448,8 +439,10 @@ struct JetSpectraEseTask {
registry.add("trackQA/after/hNClsCrossed", "", {HistType::kTH1F, {{100, 40, 140}}});
registry.add("trackQA/before/hNClsMin", "", {HistType::kTH1F, {{100, 40, 140}}});
registry.add("trackQA/after/hNClsMin", "", {HistType::kTH1F, {{100, 40, 140}}});
registry.add("trackQA/before/hDCAz", "", {HistType::kTH1F, {{50, 0, 3}}});
registry.add("trackQA/after/hDCAz", "", {HistType::kTH1F, {{50, 0, 3}}});
registry.add("trackQA/before/hDCAz", "", {HistType::kTH1F, {{50, -3, 3}}});
registry.add("trackQA/after/hDCAz", "", {HistType::kTH1F, {{50, -3, 3}}});
registry.add("trackQA/before/hDCAxy", "", {HistType::kTH1F, {{50, -1, 1}}});
registry.add("trackQA/after/hDCAxy", "", {HistType::kTH1F, {{50, -1, 1}}});
registry.add("trackQA/before/hChi2TPC", "", {HistType::kTH1F, {{50, 0, 6}}});
registry.add("trackQA/after/hChi2TPC", "", {HistType::kTH1F, {{50, 0, 6}}});
registry.add("trackQA/before/hChi2ITS", "", {HistType::kTH1F, {{50, 0, 38}}});
Expand Down Expand Up @@ -641,28 +634,27 @@ struct JetSpectraEseTask {
return;
}
const auto efficiencyName = fmt::format("efficiency_{}", std::max(cfgSystFlag.value, 0));
if (!cfgEfficiency.value.empty()) {
cfg.hEffList = ccdb->getForTimeStamp<TList>(cfgEfficiency, timestamp);
if (cfg.hEffList == nullptr) {
LOGF(fatal, "Could not load track efficiency list from %s", cfgEfficiency.value.c_str());
}
cfg.hEff = dynamic_cast<TH1F*>(cfg.hEffList->FindObject(efficiencyName.c_str()));
if (cfg.hEff == nullptr) {
LOGF(fatal, "Could not find %s as TH1F in track efficiency list %s", efficiencyName.c_str(), cfgEfficiency.value.c_str());
}
LOGF(info, "Loaded tracking efficiency %s from %s", efficiencyName.c_str(), cfgEfficiency.value.c_str());
}
if (!cfgEfficiency3D.value.empty()) {
cfg.h3EffList = ccdb->getForTimeStamp<TList>(cfgEfficiency3D, timestamp);
if (!efficiencyFiles.cfgEfficiency3D.value.empty()) {
cfg.h3EffList = ccdb->getForTimeStamp<TList>(efficiencyFiles.cfgEfficiency3D, timestamp);
if (cfg.h3EffList == nullptr) {
LOGF(fatal, "Could not load 3D track efficiency list from %s", cfgEfficiency3D.value.c_str());
LOGF(fatal, "Could not load 3D track efficiency list from %s", efficiencyFiles.cfgEfficiency3D.value.c_str());
}
cfg.h3Eff = dynamic_cast<TH3F*>(cfg.h3EffList->FindObject(efficiencyName.c_str()));
if (cfg.h3Eff == nullptr) {
LOGF(fatal, "Could not find %s as TH3F in 3D track efficiency list %s", efficiencyName.c_str(), cfgEfficiency3D.value.c_str());
LOGF(fatal, "Could not find %s as TH3F in 3D track efficiency list %s", efficiencyName.c_str(), efficiencyFiles.cfgEfficiency3D.value.c_str());
}
LOGF(info, "Loaded 3D tracking efficiency %s from %s", efficiencyName.c_str(), cfgEfficiency3D.value.c_str());
LOGF(info, "Loaded 3D tracking efficiency %s from %s", efficiencyName.c_str(), efficiencyFiles.cfgEfficiency3D.value.c_str());
cfg.is3D = true;
} else if (!efficiencyFiles.cfgEfficiency.value.empty()) {
cfg.hEffList = ccdb->getForTimeStamp<TList>(efficiencyFiles.cfgEfficiency, timestamp);
if (cfg.hEffList == nullptr) {
LOGF(fatal, "Could not load track efficiency list from %s", efficiencyFiles.cfgEfficiency.value.c_str());
}
cfg.hEff = dynamic_cast<TH1F*>(cfg.hEffList->FindObject(efficiencyName.c_str()));
if (cfg.hEff == nullptr) {
LOGF(fatal, "Could not find %s as TH1F in track efficiency list %s", efficiencyName.c_str(), efficiencyFiles.cfgEfficiency.value.c_str());
}
LOGF(info, "Loaded tracking efficiency %s from %s", efficiencyName.c_str(), efficiencyFiles.cfgEfficiency.value.c_str());
}
cfg.isLoaded = true;
}
Expand Down Expand Up @@ -729,6 +721,7 @@ struct JetSpectraEseTask {
registry.fill(HIST("eventQA/hRho"), centrality, collision.rho());
registry.fill(HIST("eventQA/hCentralityAnalyzed"), centrality);
int lRndInd = fRndm->Integer(NumSubSmpl);
hCentralityAnalyzedSub[lRndInd]->Fill(centrality);

auto corrL = [&](const auto& j) { return j.pt() - evalRho(rhoFit.get(), jetR, j.phi(), collision.rho()) * j.area(); };
for (const auto& jet : jets) {
Expand Down Expand Up @@ -790,6 +783,7 @@ struct JetSpectraEseTask {
registry.fill(HIST("trackQA/before/hNClsCrossed"), trk.tpcNClsCrossedRows());
registry.fill(HIST("trackQA/before/hNClsMin"), trk.tpcNClsFound());
registry.fill(HIST("trackQA/before/hDCAz"), trk.dcaZ());
registry.fill(HIST("trackQA/before/hDCAxy"), trk.dcaXY());
registry.fill(HIST("trackQA/before/hChi2TPC"), trk.tpcChi2NCl());
registry.fill(HIST("trackQA/before/hChi2ITS"), trk.itsChi2NCl());
if (!jetderiveddatautilities::selectTrack(track, trackSelection)) {
Expand All @@ -804,6 +798,7 @@ struct JetSpectraEseTask {
registry.fill(HIST("trackQA/after/hNClsCrossed"), trk.tpcNClsCrossedRows());
registry.fill(HIST("trackQA/after/hNClsMin"), trk.tpcNClsFound());
registry.fill(HIST("trackQA/after/hDCAz"), trk.dcaZ());
registry.fill(HIST("trackQA/after/hDCAxy"), trk.dcaXY());
registry.fill(HIST("trackQA/after/hChi2TPC"), trk.tpcChi2NCl());
registry.fill(HIST("trackQA/after/hChi2ITS"), trk.itsChi2NCl());
registry.fill(HIST("h3CenttrPhiPsi2"), centrality, RecoDecay::constrainAngle(track.phi() - psi.psi2, -o2::constants::math::PI), qPerc[0]);
Expand Down Expand Up @@ -873,7 +868,7 @@ struct JetSpectraEseTask {
registry.fill(HIST("eventQA/hEventCounterMixed"), kRhoLocal);
if (fLeadJetPtCut) {
if (!isAcceptedLeadingJet<false>(c1, jets1, centrality)) {
return;
continue;
}
}
registry.fill(HIST("eventQA/hEventCounterMixed"), kLeadJetCut);
Expand Down Expand Up @@ -1271,7 +1266,7 @@ struct JetSpectraEseTask {
if (mcCollision.size() < 1) {
return;
}
if (collisions.size() < 1) {
if (collisions.size() != 1) {
return;
}
if (!(std::abs(mcCollision.posZ()) < systCuts.vertexZCut)) {
Expand Down Expand Up @@ -1589,7 +1584,10 @@ struct JetSpectraEseTask {
modulationFit->FixParameter(2, (ep.psi2 < 0) ? RecoDecay::constrainAngle(ep.psi2) : ep.psi2);
modulationFit->FixParameter(4, (ep.psi3 < 0) ? RecoDecay::constrainAngle(ep.psi3) : ep.psi3);

hPhiPt->Fit(modulationFit.get(), "QN", "", 0, o2::constants::math::TwoPI);
const int fitStatus = hPhiPt->Fit(modulationFit.get(), "QN", "", 0, o2::constants::math::TwoPI);
if (fitStatus != 0) {
return nullptr;
}

if constexpr (fillHist) {
registry.fill(HIST("eventQA/hfitPar0"), getCentrality(col), modulationFit->GetParameter(0));
Expand All @@ -1603,22 +1601,13 @@ struct JetSpectraEseTask {
return nullptr;
}

double chi2{0.};
for (int i{0}; i < hPhiPt->GetXaxis()->GetNbins(); i++) {
if (hPhiPt->GetBinContent(i + 1) <= 0.) {
continue;
}
chi2 += std::pow((hPhiPt->GetBinContent(i + 1) - modulationFit->Eval(hPhiPt->GetXaxis()->GetBinCenter(1 + i))), 2) / hPhiPt->GetBinContent(i + 1);
}

int nDF{1};
int numParams{2};
nDF = static_cast<int>(modulationFit->GetXaxis()->GetNbins()) - numParams;
const double chi2 = modulationFit->GetChisquare();
const int nDF = modulationFit->GetNDF();
if (nDF <= 0) {
return nullptr;
}

auto cDF = 1. - TMath::Gamma(nDF, chi2);
const double cDF = TMath::Prob(chi2, nDF);
if constexpr (fillHist) {
registry.fill(HIST("eventQA/hRhoPhiCheck"), 0.5);
}
Expand Down Expand Up @@ -1798,11 +1787,11 @@ struct JetSpectraEseTask {
}
}
}
registry.fill(HIST("hCentRhoRandomConeRndTrackDirwoOneLeadingJet"), getCentrality(collision), randomConePtWithoutOneLeadJet - o2::constants::math::PI * randomConeR * randomConeR * rho, dPhiRC, qPerc[0]);
registry.fill(HIST("hCentRhoRandomConeRndTrackDirwoTwoLeadingJet"), getCentrality(collision), randomConePtWithoutTwoLeadJet - o2::constants::math::PI * randomConeR * randomConeR * rho, dPhiRC, qPerc[0]);
}
registry.fill(HIST("h3CentdeltapTRndmConePhi_rhovsphi"), getCentrality(collision), randomConePt - o2::constants::math::PI * randomConeR * randomConeR * collision.rho(), dPhiRC);
registry.fill(HIST("h3CentdeltapTRndmConePhi_localrhovsphi"), getCentrality(collision), randomConePt - o2::constants::math::PI * randomConeR * randomConeR * rho, dPhiRC);
registry.fill(HIST("hCentRhoRandomConeRndTrackDirwoOneLeadingJet"), getCentrality(collision), randomConePtWithoutOneLeadJet - o2::constants::math::PI * randomConeR * randomConeR * rho, dPhiRC, qPerc[0]);
registry.fill(HIST("hCentRhoRandomConeRndTrackDirwoTwoLeadingJet"), getCentrality(collision), randomConePtWithoutTwoLeadJet - o2::constants::math::PI * randomConeR * randomConeR * rho, dPhiRC, qPerc[0]);
}
template <typename TTracks, typename TJets>
bool isTrackInJet(TTracks const& track, TJets const& jet)
Expand Down Expand Up @@ -1970,16 +1959,18 @@ struct JetSpectraEseTask {
if (track.eta() <= trackCuts.trackEtaMin || track.eta() >= trackCuts.trackEtaMax) {
return false;
}
if (track.tpcNClsCrossedRows() <= systCuts.nTPCXrows) {
if (track.tpcNClsCrossedRows() < systCuts.nTPCXrows) {
return false;
}
if (track.tpcNClsFound() <= systCuts.nTPCCls) {
if (track.tpcNClsFound() < systCuts.nTPCCls) {
return false;
}
if (std::fabs(track.dcaZ()) >= systCuts.trackDCAzMax) {
return false;
}
// cfgSystFlag=6 reserves systCuts.dcaXYSigmaMax for a future DCAxy significance cut.
if (std::fabs(track.dcaXY()) > (systCuts.dcaXYSigmaMax / 7.f) * (0.0105f + 0.0350f / std::pow(track.pt(), 1.1f))) {
return false;
}
if (track.tpcChi2NCl() > systCuts.chi2PrTPCcls) {
return false;
}
Expand Down
Loading