Skip to content

Commit bbaf53c

Browse files
joachimckhjoachimckh
committed
fix few bugs, introduce dcaxy cut, event counter for subsamples
1 parent abd5072 commit bbaf53c

1 file changed

Lines changed: 52 additions & 61 deletions

File tree

‎PWGJE/Tasks/jetSpectraEseTask.cxx‎

Lines changed: 52 additions & 61 deletions
Original file line numberDiff line numberDiff line change
@@ -70,8 +70,11 @@ using namespace o2::framework::expressions;
7070

7171
struct JetSpectraEseTask {
7272
Configurable<int> cfgSystFlag{"cfgSystFlag", -1, "-1 = use configurable cuts, 0 = nominal defaults, >0 = systematic cut variation"};
73-
Configurable<std::string> cfgEfficiency{"cfgEfficiency", "", "CCDB path to efficiency"};
74-
Configurable<std::string> cfgEfficiency3D{"cfgEfficiency3D", "", "CCDB path to 3D efficiency"};
73+
74+
struct : ConfigurableGroup {
75+
Configurable<std::string> cfgEfficiency{"cfgEfficiency", "", "CCDB path to efficiency"};
76+
Configurable<std::string> cfgEfficiency3D{"cfgEfficiency3D", "", "CCDB path to 3D efficiency"};
77+
} efficiencyFiles;
7578
Configurable<float> jetPtMin{"jetPtMin", 5.0, "minimum jet pT cut"};
7679
Configurable<float> jetR{"jetR", 0.2, "jet resolution parameter"};
7780
Configurable<float> randomConeR{"randomConeR", 0.4, "size of random Cone for estimating background fluctuations"};
@@ -101,8 +104,9 @@ struct JetSpectraEseTask {
101104
Configurable<float> trackDCAzMax{"trackDCAzMax", 2.0, "maximum DCAz for tracks"};
102105
Configurable<float> cfgNTPCXrows{"cfgNTPCXrows", 70, "minimum number of crossed rows in the TPC for tracks"};
103106
Configurable<float> cfgNTPCCls{"cfgNTPCCls", 50, "minimum number of clusters in the TPC for tracks"};
104-
Configurable<float> cfgChi2PrTPCcls{"cfgChi2PrTPCcls", 2.5, "cut for chi2 per TPC cluster for tracks"};
107+
Configurable<float> cfgChi2PrTPCcls{"cfgChi2PrTPCcls", 4.0, "cut for chi2 per TPC cluster for tracks"};
105108
Configurable<float> cfgChi2PrITScls{"cfgChi2PrITScls", 36, "cut for chi2 per ITS cluster for tracks"};
109+
Configurable<float> dcaXYSigmaMax{"dcaXYSigmaMax", 7.0, "maximum DCAxy sigma significance for tracks"};
106110
} trackCuts;
107111

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

174171
struct SystematicCuts {
175172
float vertexZCut = 10.0f;
176173
float trackDCAzMax = 2.0f;
177-
float dcaXYSigmaMax = -1.0f;
174+
float dcaXYSigmaMax = 7.0f;
178175
float nTPCXrows = 70.0f;
179176
float nTPCCls = 50.0f;
180-
float chi2PrTPCcls = 2.5f;
177+
float chi2PrTPCcls = 4.0f;
181178
float chi2PrITScls = 36.0f;
182179
} systCuts;
183180

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

256254
void applySystematicPreset()
257255
{
258256
systCuts.vertexZCut = vertexZCut.value;
259257
systCuts.trackDCAzMax = trackCuts.trackDCAzMax.value;
260-
systCuts.dcaXYSigmaMax = NoDCAxySigmaCut;
258+
systCuts.dcaXYSigmaMax = trackCuts.dcaXYSigmaMax.value;
261259
systCuts.nTPCXrows = trackCuts.cfgNTPCXrows.value;
262260
systCuts.nTPCCls = trackCuts.cfgNTPCCls.value;
263261
systCuts.chi2PrTPCcls = trackCuts.cfgChi2PrTPCcls.value;
@@ -268,15 +266,6 @@ struct JetSpectraEseTask {
268266
LOGF(info, "JetSpectraEseTask::applySystematicPreset() - using configurable cut values");
269267
return;
270268
}
271-
272-
systCuts.vertexZCut = NominalVertexZCut;
273-
systCuts.trackDCAzMax = NominalTrackDCAzMax;
274-
systCuts.dcaXYSigmaMax = NoDCAxySigmaCut;
275-
systCuts.nTPCXrows = NominalNTPCXrows;
276-
systCuts.nTPCCls = NominalNTPCCls;
277-
systCuts.chi2PrTPCcls = NominalChi2PrTPCcls;
278-
systCuts.chi2PrITScls = NominalChi2PrITScls;
279-
280269
switch (flag) {
281270
case 0:
282271
break;
@@ -305,7 +294,7 @@ struct JetSpectraEseTask {
305294
systCuts.nTPCCls = 70.0f;
306295
break;
307296
case 9:
308-
systCuts.chi2PrTPCcls = 4.0f;
297+
systCuts.chi2PrTPCcls = 2.5f;
309298
break;
310299
case 10:
311300
systCuts.nTPCXrows = 80.0f;
@@ -435,6 +424,8 @@ struct JetSpectraEseTask {
435424
for (int i = 0; i < NumSubSmpl; ++i) {
436425
std::string n = fmt::format("subsamples/thn_jethad_corr_same_subsample{}", i);
437426
hSameSub[i] = registry.add<THnSparse>(n, "same;...", o2::framework::HistType::kTHnSparseF, axes);
427+
std::string n2 = fmt::format("subsamples/hCentralityAnalyzed_subsample{}", i);
428+
hCentralityAnalyzedSub[i] = registry.add<TH1>(n2, ";Centrality;entries", o2::framework::HistType::kTH1F, {{centAxis}});
438429
}
439430
registry.add("hNtrig", "", {HistType::kTHnSparseF, {{centAxis}, {jetPtAxis}, {dPhiAxis}, {eseAxis}}});
440431

@@ -448,8 +439,10 @@ struct JetSpectraEseTask {
448439
registry.add("trackQA/after/hNClsCrossed", "", {HistType::kTH1F, {{100, 40, 140}}});
449440
registry.add("trackQA/before/hNClsMin", "", {HistType::kTH1F, {{100, 40, 140}}});
450441
registry.add("trackQA/after/hNClsMin", "", {HistType::kTH1F, {{100, 40, 140}}});
451-
registry.add("trackQA/before/hDCAz", "", {HistType::kTH1F, {{50, 0, 3}}});
452-
registry.add("trackQA/after/hDCAz", "", {HistType::kTH1F, {{50, 0, 3}}});
442+
registry.add("trackQA/before/hDCAz", "", {HistType::kTH1F, {{50, -3, 3}}});
443+
registry.add("trackQA/after/hDCAz", "", {HistType::kTH1F, {{50, -3, 3}}});
444+
registry.add("trackQA/before/hDCAxy", "", {HistType::kTH1F, {{50, -1, 1}}});
445+
registry.add("trackQA/after/hDCAxy", "", {HistType::kTH1F, {{50, -1, 1}}});
453446
registry.add("trackQA/before/hChi2TPC", "", {HistType::kTH1F, {{50, 0, 6}}});
454447
registry.add("trackQA/after/hChi2TPC", "", {HistType::kTH1F, {{50, 0, 6}}});
455448
registry.add("trackQA/before/hChi2ITS", "", {HistType::kTH1F, {{50, 0, 38}}});
@@ -641,28 +634,27 @@ struct JetSpectraEseTask {
641634
return;
642635
}
643636
const auto efficiencyName = fmt::format("efficiency_{}", std::max(cfgSystFlag.value, 0));
644-
if (!cfgEfficiency.value.empty()) {
645-
cfg.hEffList = ccdb->getForTimeStamp<TList>(cfgEfficiency, timestamp);
646-
if (cfg.hEffList == nullptr) {
647-
LOGF(fatal, "Could not load track efficiency list from %s", cfgEfficiency.value.c_str());
648-
}
649-
cfg.hEff = dynamic_cast<TH1F*>(cfg.hEffList->FindObject(efficiencyName.c_str()));
650-
if (cfg.hEff == nullptr) {
651-
LOGF(fatal, "Could not find %s as TH1F in track efficiency list %s", efficiencyName.c_str(), cfgEfficiency.value.c_str());
652-
}
653-
LOGF(info, "Loaded tracking efficiency %s from %s", efficiencyName.c_str(), cfgEfficiency.value.c_str());
654-
}
655-
if (!cfgEfficiency3D.value.empty()) {
656-
cfg.h3EffList = ccdb->getForTimeStamp<TList>(cfgEfficiency3D, timestamp);
637+
if (!efficiencyFiles.cfgEfficiency3D.value.empty()) {
638+
cfg.h3EffList = ccdb->getForTimeStamp<TList>(efficiencyFiles.cfgEfficiency3D, timestamp);
657639
if (cfg.h3EffList == nullptr) {
658-
LOGF(fatal, "Could not load 3D track efficiency list from %s", cfgEfficiency3D.value.c_str());
640+
LOGF(fatal, "Could not load 3D track efficiency list from %s", efficiencyFiles.cfgEfficiency3D.value.c_str());
659641
}
660642
cfg.h3Eff = dynamic_cast<TH3F*>(cfg.h3EffList->FindObject(efficiencyName.c_str()));
661643
if (cfg.h3Eff == nullptr) {
662-
LOGF(fatal, "Could not find %s as TH3F in 3D track efficiency list %s", efficiencyName.c_str(), cfgEfficiency3D.value.c_str());
644+
LOGF(fatal, "Could not find %s as TH3F in 3D track efficiency list %s", efficiencyName.c_str(), efficiencyFiles.cfgEfficiency3D.value.c_str());
663645
}
664-
LOGF(info, "Loaded 3D tracking efficiency %s from %s", efficiencyName.c_str(), cfgEfficiency3D.value.c_str());
646+
LOGF(info, "Loaded 3D tracking efficiency %s from %s", efficiencyName.c_str(), efficiencyFiles.cfgEfficiency3D.value.c_str());
665647
cfg.is3D = true;
648+
} else if (!efficiencyFiles.cfgEfficiency.value.empty()) {
649+
cfg.hEffList = ccdb->getForTimeStamp<TList>(efficiencyFiles.cfgEfficiency, timestamp);
650+
if (cfg.hEffList == nullptr) {
651+
LOGF(fatal, "Could not load track efficiency list from %s", efficiencyFiles.cfgEfficiency.value.c_str());
652+
}
653+
cfg.hEff = dynamic_cast<TH1F*>(cfg.hEffList->FindObject(efficiencyName.c_str()));
654+
if (cfg.hEff == nullptr) {
655+
LOGF(fatal, "Could not find %s as TH1F in track efficiency list %s", efficiencyName.c_str(), efficiencyFiles.cfgEfficiency.value.c_str());
656+
}
657+
LOGF(info, "Loaded tracking efficiency %s from %s", efficiencyName.c_str(), efficiencyFiles.cfgEfficiency.value.c_str());
666658
}
667659
cfg.isLoaded = true;
668660
}
@@ -729,6 +721,7 @@ struct JetSpectraEseTask {
729721
registry.fill(HIST("eventQA/hRho"), centrality, collision.rho());
730722
registry.fill(HIST("eventQA/hCentralityAnalyzed"), centrality);
731723
int lRndInd = fRndm->Integer(NumSubSmpl);
724+
hCentralityAnalyzedSub[lRndInd]->Fill(centrality);
732725

733726
auto corrL = [&](const auto& j) { return j.pt() - evalRho(rhoFit.get(), jetR, j.phi(), collision.rho()) * j.area(); };
734727
for (const auto& jet : jets) {
@@ -790,6 +783,7 @@ struct JetSpectraEseTask {
790783
registry.fill(HIST("trackQA/before/hNClsCrossed"), trk.tpcNClsCrossedRows());
791784
registry.fill(HIST("trackQA/before/hNClsMin"), trk.tpcNClsFound());
792785
registry.fill(HIST("trackQA/before/hDCAz"), trk.dcaZ());
786+
registry.fill(HIST("trackQA/before/hDCAxy"), trk.dcaXY());
793787
registry.fill(HIST("trackQA/before/hChi2TPC"), trk.tpcChi2NCl());
794788
registry.fill(HIST("trackQA/before/hChi2ITS"), trk.itsChi2NCl());
795789
if (!jetderiveddatautilities::selectTrack(track, trackSelection)) {
@@ -804,6 +798,7 @@ struct JetSpectraEseTask {
804798
registry.fill(HIST("trackQA/after/hNClsCrossed"), trk.tpcNClsCrossedRows());
805799
registry.fill(HIST("trackQA/after/hNClsMin"), trk.tpcNClsFound());
806800
registry.fill(HIST("trackQA/after/hDCAz"), trk.dcaZ());
801+
registry.fill(HIST("trackQA/after/hDCAxy"), trk.dcaXY());
807802
registry.fill(HIST("trackQA/after/hChi2TPC"), trk.tpcChi2NCl());
808803
registry.fill(HIST("trackQA/after/hChi2ITS"), trk.itsChi2NCl());
809804
registry.fill(HIST("h3CenttrPhiPsi2"), centrality, RecoDecay::constrainAngle(track.phi() - psi.psi2, -o2::constants::math::PI), qPerc[0]);
@@ -873,7 +868,7 @@ struct JetSpectraEseTask {
873868
registry.fill(HIST("eventQA/hEventCounterMixed"), kRhoLocal);
874869
if (fLeadJetPtCut) {
875870
if (!isAcceptedLeadingJet<false>(c1, jets1, centrality)) {
876-
return;
871+
continue;
877872
}
878873
}
879874
registry.fill(HIST("eventQA/hEventCounterMixed"), kLeadJetCut);
@@ -1271,7 +1266,7 @@ struct JetSpectraEseTask {
12711266
if (mcCollision.size() < 1) {
12721267
return;
12731268
}
1274-
if (collisions.size() < 1) {
1269+
if (collisions.size() != 1) {
12751270
return;
12761271
}
12771272
if (!(std::abs(mcCollision.posZ()) < systCuts.vertexZCut)) {
@@ -1589,7 +1584,10 @@ struct JetSpectraEseTask {
15891584
modulationFit->FixParameter(2, (ep.psi2 < 0) ? RecoDecay::constrainAngle(ep.psi2) : ep.psi2);
15901585
modulationFit->FixParameter(4, (ep.psi3 < 0) ? RecoDecay::constrainAngle(ep.psi3) : ep.psi3);
15911586

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

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

1606-
double chi2{0.};
1607-
for (int i{0}; i < hPhiPt->GetXaxis()->GetNbins(); i++) {
1608-
if (hPhiPt->GetBinContent(i + 1) <= 0.) {
1609-
continue;
1610-
}
1611-
chi2 += std::pow((hPhiPt->GetBinContent(i + 1) - modulationFit->Eval(hPhiPt->GetXaxis()->GetBinCenter(1 + i))), 2) / hPhiPt->GetBinContent(i + 1);
1612-
}
1613-
1614-
int nDF{1};
1615-
int numParams{2};
1616-
nDF = static_cast<int>(modulationFit->GetXaxis()->GetNbins()) - numParams;
1604+
const double chi2 = modulationFit->GetChisquare();
1605+
const int nDF = modulationFit->GetNDF();
16171606
if (nDF <= 0) {
16181607
return nullptr;
16191608
}
16201609

1621-
auto cDF = 1. - TMath::Gamma(nDF, chi2);
1610+
const double cDF = TMath::Prob(chi2, nDF);
16221611
if constexpr (fillHist) {
16231612
registry.fill(HIST("eventQA/hRhoPhiCheck"), 0.5);
16241613
}
@@ -1798,11 +1787,11 @@ struct JetSpectraEseTask {
17981787
}
17991788
}
18001789
}
1790+
registry.fill(HIST("hCentRhoRandomConeRndTrackDirwoOneLeadingJet"), getCentrality(collision), randomConePtWithoutOneLeadJet - o2::constants::math::PI * randomConeR * randomConeR * rho, dPhiRC, qPerc[0]);
1791+
registry.fill(HIST("hCentRhoRandomConeRndTrackDirwoTwoLeadingJet"), getCentrality(collision), randomConePtWithoutTwoLeadJet - o2::constants::math::PI * randomConeR * randomConeR * rho, dPhiRC, qPerc[0]);
18011792
}
18021793
registry.fill(HIST("h3CentdeltapTRndmConePhi_rhovsphi"), getCentrality(collision), randomConePt - o2::constants::math::PI * randomConeR * randomConeR * collision.rho(), dPhiRC);
18031794
registry.fill(HIST("h3CentdeltapTRndmConePhi_localrhovsphi"), getCentrality(collision), randomConePt - o2::constants::math::PI * randomConeR * randomConeR * rho, dPhiRC);
1804-
registry.fill(HIST("hCentRhoRandomConeRndTrackDirwoOneLeadingJet"), getCentrality(collision), randomConePtWithoutOneLeadJet - o2::constants::math::PI * randomConeR * randomConeR * rho, dPhiRC, qPerc[0]);
1805-
registry.fill(HIST("hCentRhoRandomConeRndTrackDirwoTwoLeadingJet"), getCentrality(collision), randomConePtWithoutTwoLeadJet - o2::constants::math::PI * randomConeR * randomConeR * rho, dPhiRC, qPerc[0]);
18061795
}
18071796
template <typename TTracks, typename TJets>
18081797
bool isTrackInJet(TTracks const& track, TJets const& jet)
@@ -1970,16 +1959,18 @@ struct JetSpectraEseTask {
19701959
if (track.eta() <= trackCuts.trackEtaMin || track.eta() >= trackCuts.trackEtaMax) {
19711960
return false;
19721961
}
1973-
if (track.tpcNClsCrossedRows() <= systCuts.nTPCXrows) {
1962+
if (track.tpcNClsCrossedRows() < systCuts.nTPCXrows) {
19741963
return false;
19751964
}
1976-
if (track.tpcNClsFound() <= systCuts.nTPCCls) {
1965+
if (track.tpcNClsFound() < systCuts.nTPCCls) {
19771966
return false;
19781967
}
19791968
if (std::fabs(track.dcaZ()) >= systCuts.trackDCAzMax) {
19801969
return false;
19811970
}
1982-
// cfgSystFlag=6 reserves systCuts.dcaXYSigmaMax for a future DCAxy significance cut.
1971+
if (std::fabs(track.dcaXY()) > (systCuts.dcaXYSigmaMax / 7.f) * (0.0105f + 0.0350f / std::pow(track.pt(), 1.1f))) {
1972+
return false;
1973+
}
19831974
if (track.tpcChi2NCl() > systCuts.chi2PrTPCcls) {
19841975
return false;
19851976
}

0 commit comments

Comments
 (0)