Skip to content

Commit 2d26511

Browse files
LuzhiyonggLuzhiyongg
committed
debug pt-dependence efficiency
1 parent 3199b1d commit 2d26511

File tree

2 files changed

+71
-21
lines changed

2 files changed

+71
-21
lines changed

‎PWGCF/Flow/Tasks/flowMc.cxx‎

Lines changed: 49 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -78,7 +78,7 @@ struct FlowMc {
7878
O2_DEFINE_CONFIGURABLE(cfgFlowEfficiency, std::string, "", "CCDB path to efficiency object")
7979
O2_DEFINE_CONFIGURABLE(cfgCentVsIPTruth, std::string, "", "CCDB path to centrality vs IP truth")
8080
O2_DEFINE_CONFIGURABLE(cfgIsGlobalTrack, bool, false, "Use global tracks instead of hasTPC&&hasITS")
81-
O2_DEFINE_CONFIGURABLE(cfgK0Lambda0Enabled, bool, false, "Add K0 and Lambda0")
81+
O2_DEFINE_CONFIGURABLE(cfgK0Lambda0Enabled, bool, true, "Add K0 and Lambda0")
8282
O2_DEFINE_CONFIGURABLE(cfgFlowCumulantEnabled, bool, false, "switch of calculating flow")
8383
O2_DEFINE_CONFIGURABLE(cfgFlowCumulantNbootstrap, int, 30, "Number of subsamples")
8484
O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrUse, bool, false, "Use track density efficiency correction")
@@ -96,7 +96,6 @@ struct FlowMc {
9696

9797
Configurable<std::vector<double>> cfgTrackDensityP0{"cfgTrackDensityP0", std::vector<double>{0.6003720411, 0.6152630970, 0.6288860646, 0.6360694031, 0.6409494798, 0.6450540203, 0.6482117301, 0.6512592056, 0.6640008690, 0.6862631416, 0.7005738691, 0.7106567432, 0.7170728333}, "parameter 0 for track density efficiency correction"};
9898
Configurable<std::vector<double>> cfgTrackDensityP1{"cfgTrackDensityP1", std::vector<double>{-1.007592e-05, -8.932635e-06, -9.114538e-06, -1.054818e-05, -1.220212e-05, -1.312304e-05, -1.376433e-05, -1.412813e-05, -1.289562e-05, -1.050065e-05, -8.635725e-06, -7.380821e-06, -6.201250e-06}, "parameter 1 for track density efficiency correction"};
99-
float maxEta = 0.8;
10099

101100
ConfigurableAxis axisB{"axisB", {100, 0.0f, 20.0f}, ""};
102101
ConfigurableAxis axisCentrality{"axisCentrality", {VARIABLE_WIDTH, 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90}, "X axis for histograms"};
@@ -162,9 +161,10 @@ struct FlowMc {
162161
const AxisSpec axisEta{20, -1., 1., "#eta"};
163162
const AxisSpec axisCounter{1, 0, +1, ""};
164163
// QA histograms
165-
histos.add<TH1>("mcEventCounter", "Monte Carlo Truth EventCounter", HistType::kTH1F, {axisCounter});
166-
histos.add<TH1>("numberOfRecoCollisions", "numberOfRecoCollisions", HistType::kTH1F, {{10, -0.5f, 9.5f}});
167-
histos.add<TH1>("RecoEventCounter", "Reconstruction EventCounter", HistType::kTH1F, {axisCounter});
164+
histos.add<TH1>("mcEventCounter", "Monte Carlo Truth EventCounter", HistType::kTH1F, {{5, 0, 5}});
165+
histos.add<TH1>("numberOfRecoCollisions", "numberOfRecoCollisions", HistType::kTH1F, {{10, -1.5f, 8.5f}});
166+
histos.add<TH1>("RecoEventCounter", "Reconstruction EventCounter", HistType::kTH1F, {{5, 0, 5}});
167+
histos.add<TH1>("RecoProcessEventCounter", "Reconstruction EventCounter", HistType::kTH1F, {{5, 0, 5}});
168168
histos.add<TH1>("hnTPCClu", "Number of found TPC clusters", HistType::kTH1D, {{100, 40, 180}});
169169
histos.add<TH1>("hnITSClu", "Number of found ITS clusters", HistType::kTH1D, {{100, 0, 20}});
170170
// pT histograms
@@ -209,6 +209,8 @@ struct FlowMc {
209209
histos.add<TH1>("hPtMCGen", "Monte Carlo Truth; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
210210
histos.add<TH3>("hEtaPtVtxzMCGen", "Monte Carlo Truth; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}});
211211
histos.add<TH1>("hPtMCGlobal", "Monte Carlo Global; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
212+
histos.add<TH1>("hPtReco", "Monte Carlo Reco Global; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
213+
histos.add<TH1>("hPtReco_PhysicalPrimary", "Monte Carlo Reco Global; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
212214
histos.add<TH3>("hEtaPtVtxzMCGlobal", "Monte Carlo Global; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}});
213215
histos.add<TH1>("hPhiWeightedTrDen", "corrected #phi distribution, considering track density", {HistType::kTH1D, {axisPhi}});
214216

@@ -447,9 +449,9 @@ struct FlowMc {
447449
return myTrackSel.IsSelected(track);
448450
}
449451

450-
void process(FilteredMcCollisions::iterator const& mcCollision, aod::BCsWithTimestamps const&, soa::SmallGroups<soa::Join<aod::McCollisionLabels, aod::Collisions, aod::EvSels>> const& collisions, FilteredMcParticles const& mcParticles, FilteredTracks const&)
452+
void processMCTrue(FilteredMcCollisions::iterator const& mcCollision, aod::BCsWithTimestamps const&, soa::SmallGroups<soa::Join<aod::McCollisionLabels, aod::Collisions, aod::EvSels>> const& collisions, FilteredMcParticles const& mcParticles, FilteredTracks const&)
451453
{
452-
454+
histos.fill(HIST("mcEventCounter"), 0.5);
453455
float imp = mcCollision.impactParameter();
454456
float evPhi = mcCollision.eventPlaneAngle();
455457
float vtxz = mcCollision.posZ();
@@ -465,21 +467,23 @@ struct FlowMc {
465467
loadCorrections(bc.timestamp());
466468

467469
if (collisions.size() > -1) {
468-
histos.fill(HIST("mcEventCounter"), 0.5);
469470
histos.fill(HIST("numberOfRecoCollisions"), collisions.size()); // number of times coll was reco-ed
471+
histos.fill(HIST("RecoEventCounter"), 0.5);
470472
if (cfgRecoEvRejectMC) {
471473
if (collisions.size() != 1) { // only pass those have one reconstruction event
472474
return;
473475
}
476+
histos.fill(HIST("RecoEventCounter"), 1.5);
474477
for (auto const& collision : collisions) {
475478
if (!eventSelected(collision))
476479
return;
477480
}
481+
histos.fill(HIST("RecoEventCounter"), 2.5);
478482
}
479-
histos.fill(HIST("RecoEventCounter"), 0.5);
480483
}
484+
histos.fill(HIST("RecoEventCounter"), 3.5);
481485

482-
if (imp > minB && imp < maxB) {
486+
if (imp >= minB && imp <= maxB) {
483487
// event within range
484488
histos.fill(HIST("hImpactParameter"), imp);
485489
histos.fill(HIST("hEventPlaneAngle"), evPhi);
@@ -502,7 +506,7 @@ struct FlowMc {
502506
continue;
503507
if (!mcParticle.isPhysicalPrimary())
504508
continue;
505-
if (std::fabs(mcParticle.eta()) > maxEta) // main acceptance
509+
if (std::fabs(mcParticle.eta()) > cfgCutEta) // main acceptance
506510
continue;
507511
if (mcParticle.has_tracks()) {
508512
auto const& tracks = mcParticle.tracks_as<FilteredTracks>();
@@ -551,7 +555,7 @@ struct FlowMc {
551555

552556
if (!mcParticle.isPhysicalPrimary())
553557
continue;
554-
if (std::fabs(mcParticle.eta()) > maxEta) // main acceptance
558+
if (std::fabs(mcParticle.eta()) > cfgCutEta) // main acceptance
555559
continue;
556560

557561
float deltaPhi = mcParticle.phi() - mcCollision.eventPlaneAngle();
@@ -696,6 +700,39 @@ struct FlowMc {
696700
}
697701
histos.fill(HIST("hNchVsImpactParameter"), imp, nCh);
698702
}
703+
PROCESS_SWITCH(FlowMc, processMCTrue, "process pure simulation information", true);
704+
705+
void processReco(soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels>::iterator const& collision, aod::BCsWithTimestamps const&, FilteredTracks const& tracks, aod::McParticles const&, aod::McCollisions const&)
706+
{
707+
histos.fill(HIST("RecoProcessEventCounter"), 0.5);
708+
if (!eventSelected(collision))
709+
return;
710+
histos.fill(HIST("RecoProcessEventCounter"), 1.5);
711+
if (tracks.size() < 1)
712+
return;
713+
auto bc = collision.bc_as<aod::BCsWithTimestamps>();
714+
int runNumber = bc.runNumber();
715+
histos.fill(HIST("RecoProcessEventCounter"), 2.5);
716+
for (const auto& track : tracks) {
717+
if (!trackSelected(track))
718+
continue;
719+
histos.fill(HIST("hPtReco"), track.pt());
720+
auto mcParticle = track.mcParticle();
721+
int pdgCode = std::abs(mcParticle.pdgCode());
722+
bool extraPDGType = true;
723+
if (cfgK0Lambda0Enabled) {
724+
extraPDGType = (pdgCode != PDG_t::kK0Short && pdgCode != PDG_t::kLambda0);
725+
}
726+
if (extraPDGType && pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != kKPlus && pdgCode != PDG_t::kProton)
727+
continue;
728+
if (!mcParticle.isPhysicalPrimary())
729+
continue;
730+
if (std::fabs(mcParticle.eta()) > cfgCutEta) // main acceptance
731+
continue;
732+
histos.fill(HIST("hPtReco_PhysicalPrimary"), track.pt());
733+
}
734+
}
735+
PROCESS_SWITCH(FlowMc, processReco, "process reconstructed information", true);
699736
};
700737

701738
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

‎PWGCF/Flow/Tasks/flowTask.cxx‎

Lines changed: 22 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,7 @@
4141
#include "TList.h"
4242
#include <TF1.h>
4343
#include <TObjArray.h>
44+
#include <TPDGCode.h>
4445
#include <TProfile.h>
4546
#include <TRandom3.h>
4647

@@ -300,6 +301,7 @@ struct FlowTask {
300301
std::string hCentTitle = "Centrality distribution, Estimator " + std::to_string(cfgCentEstimator);
301302
registry.add("hCent", hCentTitle.c_str(), {HistType::kTH1D, {{100, 0, 100}}});
302303
if (doprocessMCGen) {
304+
registry.add("MCGen/hMCEventCount", "Number of Event;; Count", {HistType::kTH1D, {{5, 0, 5}}});
303305
registry.add("MCGen/MChVtxZ", "Vexter Z distribution", {HistType::kTH1D, {axisVertex}});
304306
registry.add("MCGen/MChMult", "Multiplicity distribution", {HistType::kTH1D, {{3000, 0.5, 3000.5}}});
305307
registry.add("MCGen/MChCent", hCentTitle.c_str(), {HistType::kTH1D, {{100, 0, 100}}});
@@ -360,7 +362,7 @@ struct FlowTask {
360362
if (doprocessMCGen) {
361363
registry.add("MCGen/MChPhi", "#phi distribution", {HistType::kTH1D, {axisPhi}});
362364
registry.add("MCGen/MChEta", "#eta distribution", {HistType::kTH1D, {axisEta}});
363-
registry.add("MCGen/MChPtRef", "p_{T} distribution after cut", {HistType::kTH1D, {axisPtHist}});
365+
registry.add("MCGen/MChPt", "p_{T} distribution after cut", {HistType::kTH1D, {axisPtHist}});
364366
registry.add("hMeanPtWithinGap08_MC", "mean p_{T}", {HistType::kTProfile, {axisIndependent}});
365367
for (auto i = 0; i < cfgNbootstrap; i++) {
366368
bootstrapArray[i][kMeanPtWithinGap08_MC] = registry.add<TProfile>(Form("BootstrapContainer_%d/hMeanPtWithinGap08_MC", i), "", {HistType::kTProfile, {axisIndependent}});
@@ -1242,14 +1244,26 @@ struct FlowTask {
12421244

12431245
void processMCGen(FilteredMcCollisions::iterator const& mcCollision, FilteredSmallGroupMcCollisions const& collisions, FilteredMcParticles const& mcParticles)
12441246
{
1247+
registry.fill(HIST("MCGen/hMCEventCount"), 0.5);
12451248
if (collisions.size() != 1)
12461249
return;
1250+
registry.fill(HIST("MCGen/hMCEventCount"), 1.5);
12471251

12481252
float cent = -1.;
12491253
for (const auto& collision : collisions) {
12501254
cent = getCentrality(collision);
12511255
}
12521256

1257+
if (cfgUseAdditionalEventCut) {
1258+
for (auto const& collision : collisions) {
1259+
if (!collision.sel8())
1260+
return;
1261+
if (!eventSelected(collision, mcParticles.size(), cent))
1262+
return;
1263+
}
1264+
}
1265+
registry.fill(HIST("MCGen/hMCEventCount"), 2.5);
1266+
12531267
float lRandom = fRndm->Rndm();
12541268
float vtxz = mcCollision.posZ();
12551269
registry.fill(HIST("MCGen/MChVtxZ"), vtxz);
@@ -1262,26 +1276,25 @@ struct FlowTask {
12621276
fGFW->Clear();
12631277
fFCptgen->clearVector();
12641278

1265-
if (cfgUseAdditionalEventCut) {
1266-
for (auto const& collision : collisions) {
1267-
if (!eventSelected(collision, mcParticles.size(), cent))
1268-
return;
1269-
}
1270-
}
1271-
12721279
double ptSum_Gap08 = 0.;
12731280
double count_Gap08 = 0.;
12741281
for (const auto& mcParticle : mcParticles) {
1282+
int pdgCode = std::abs(mcParticle.pdgCode());
1283+
bool extraPDGType = (pdgCode != PDG_t::kK0Short && pdgCode != PDG_t::kLambda0);
1284+
if (extraPDGType && pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != kKPlus && pdgCode != PDG_t::kProton)
1285+
continue;
12751286
if (!mcParticle.isPhysicalPrimary())
12761287
continue;
1288+
if (std::fabs(mcParticle.eta()) > cfgCutEta) // main acceptance
1289+
continue;
12771290
bool withinPtPOI = (cfgCutPtPOIMin < mcParticle.pt()) && (mcParticle.pt() < cfgCutPtPOIMax); // within POI pT range
12781291
bool withinPtRef = (cfgCutPtRefMin < mcParticle.pt()) && (mcParticle.pt() < cfgCutPtRefMax); // within RF pT range
12791292
bool withinEtaGap08 = (std::abs(mcParticle.eta()) < cfgCutEta);
12801293

1294+
registry.fill(HIST("MCGen/MChPt"), mcParticle.pt());
12811295
if (withinPtRef) {
12821296
registry.fill(HIST("MCGen/MChPhi"), mcParticle.phi());
12831297
registry.fill(HIST("MCGen/MChEta"), mcParticle.eta());
1284-
registry.fill(HIST("MCGen/MChPtRef"), mcParticle.pt());
12851298
if (withinEtaGap08) {
12861299
ptSum_Gap08 += mcParticle.pt();
12871300
count_Gap08 += 1.;

0 commit comments

Comments
 (0)