Skip to content

Commit ec38e3e

Browse files
romainschotterromainschotter
authored
Fix mcCentralityModule.h
1 parent 6a2c947 commit ec38e3e

1 file changed

Lines changed: 94 additions & 89 deletions

File tree

‎PWGLF/Utils/mcCentralityModule.h‎

Lines changed: 94 additions & 89 deletions
Original file line numberDiff line numberDiff line change
@@ -18,35 +18,40 @@
1818
#ifndef PWGLF_UTILS_MCCENTRALITYMODULE_H_
1919
#define PWGLF_UTILS_MCCENTRALITYMODULE_H
2020

21-
#include "TableHelper.h"
22-
2321
#include "PWGLF/DataModel/mcCentrality.h"
2422

25-
#include "Common/DataModel/Centrality.h"
26-
#include "Common/DataModel/Multiplicity.h"
27-
#include "Common/DataModel/TrackSelectionTables.h"
28-
2923
#include <CCDB/BasicCCDBManager.h>
30-
#include <CCDB/CcdbApi.h>
3124
#include <Framework/AnalysisDataModel.h>
32-
#include <Framework/AnalysisTask.h>
25+
#include <Framework/AnalysisHelpers.h>
26+
#include <Framework/Configurable.h>
27+
#include <Framework/Array2D.h>
28+
#include <Framework/DeviceSpec.h>
29+
#include <Framework/DataSpecUtils.h>
3330
#include <Framework/HistogramRegistry.h>
34-
#include <Framework/O2DatabasePDGPlugin.h>
31+
#include <Framework/HistogramSpec.h>
3532
#include <Framework/RunningWorkflowInfo.h>
36-
#include <Framework/StaticFor.h>
3733
#include <Framework/runDataProcessing.h>
38-
#include <ReconstructionDataFormats/Track.h>
3934

35+
#include <TH1.h>
36+
#include <TH2.h>
37+
#include <TProfile.h>
4038
#include <TRandom3.h>
41-
39+
#include <RtypesCore.h>
40+
#include <TString.h>
41+
42+
#include <concepts>
43+
#include <cstddef>
44+
#include <cstdio>
45+
#include <cstdlib>
46+
#include <cstdint>
4247
#include <memory>
4348
#include <string>
44-
#include <variant>
49+
#include <type_traits>
50+
#include <vector>
4551

4652
using namespace o2;
4753
using namespace o2::framework;
4854
using namespace o2::framework::expressions;
49-
using namespace o2::track;
5055

5156
//__________________________________________
5257
// strangeness builder module
@@ -61,16 +66,17 @@ namespace mccentrality // avoid polluting other namespaces
6166
// statics necessary for the configurables in this namespace
6267
static constexpr int nParameters = 1;
6368
static const std::vector<std::string> tableNames{
64-
"McCentFV0As",
65-
"McCentFT0Ms",
66-
"McCentFT0As",
67-
"McCentFT0Cs",
68-
"McCentFT0CVariant1s",
69-
"McCentFT0CVariant2s",
70-
"McCentFDDMs",
71-
"McCentNTPVs",
72-
"McCentNGlobals",
73-
"McCentMFTs"};
69+
"McCentFV0As",
70+
"McCentFT0Ms",
71+
"McCentFT0As",
72+
"McCentFT0Cs",
73+
"McCentFT0CVariant1s",
74+
"McCentFT0CVariant2s",
75+
"McCentFDDMs",
76+
"McCentNTPVs",
77+
"McCentNGlobals",
78+
"McCentMFTs"
79+
};
7480

7581
static constexpr int nTablesConst = 10;
7682
static const std::vector<std::string> parameterNames{"enable"};
@@ -97,19 +103,20 @@ enum tableIndex { kFV0A = 0,
97103
kNTPV,
98104
kNGlobal,
99105
kMFT,
100-
kNestimators };
106+
kNestimators};
101107

102108
static constexpr const char* DirList[] = {
103-
"FV0A",
104-
"FT0M",
105-
"FT0A",
106-
"FT0C",
107-
"FT0CVariant1",
108-
"FT0CVariant2",
109-
"FDDM",
110-
"NTPV",
111-
"NGlobal",
112-
"MFT"};
109+
"FV0A",
110+
"FT0M",
111+
"FT0A",
112+
"FT0C",
113+
"FT0CVariant1",
114+
"FT0CVariant2",
115+
"FDDM",
116+
"NTPV",
117+
"NGlobal",
118+
"MFT"
119+
};
113120

114121
// mcCentralityModule: 1st-order configurables
115122
struct coreConfigurables : o2::framework::ConfigurableGroup {
@@ -129,7 +136,7 @@ struct coreConfigurables : o2::framework::ConfigurableGroup {
129136
ConfigurableAxis binsPercentileFine{"binsPercentileFine", {VARIABLE_WIDTH, 0, 0.001, 0.01, 1.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0}, "Binning of the percentile axis"};
130137
ConfigurableAxis binsMultiplicity{"binsMultiplicity", {1000, 0, 5000}, "Binning of the multiplicity axis"};
131138

132-
// ccdb configurables
139+
//ccdb configurables
133140
Configurable<std::string> path{"path", "/tmp/InputCalibMC.root", "path to calib file or ccdb path if begins with ccdb://"};
134141

135142
// debug option
@@ -152,11 +159,11 @@ struct products : o2::framework::ProducesGroup {
152159

153160
template <typename T>
154161
concept HasMcMults = requires(typename T::iterator a) {
155-
{ a.multMCFT0A() } -> std::convertible_to<int>;
156-
{ a.multMCFT0C() } -> std::convertible_to<int>;
157-
{ a.multMCFV0A() } -> std::convertible_to<int>;
158-
{ a.multMCFDDA() } -> std::convertible_to<int>;
159-
{ a.multMCFDDC() } -> std::convertible_to<int>;
162+
{ a.multMCFT0A() } -> std::convertible_to<int>;
163+
{ a.multMCFT0C() } -> std::convertible_to<int>;
164+
{ a.multMCFV0A() } -> std::convertible_to<int>;
165+
{ a.multMCFDDA() } -> std::convertible_to<int>;
166+
{ a.multMCFDDC() } -> std::convertible_to<int>;
160167
{ a.multMCNParticlesEta08() } -> std::convertible_to<int>;
161168
{ a.multMCNParticlesEta05() } -> std::convertible_to<int>;
162169
};
@@ -371,8 +378,7 @@ struct BuilderModule {
371378
}
372379

373380
template <typename THist>
374-
THist* getHist(const char* name)
375-
{
381+
THist* getHist(const char* name) {
376382
if (!this->MCCentralityCalibObjects) {
377383
return (THist*)0x0;
378384
}
@@ -390,18 +396,17 @@ struct BuilderModule {
390396
return hist;
391397
}
392398

393-
TH1D* extractCentralityCalibration(int idx, const char* name, bool reverse = false)
394-
{
399+
TH1D* extractCentralityCalibration(int idx, const char* name, bool reverse = false) {
395400

396-
auto CalibMC = [&](TString estimator) {
397-
printf("\e[1;31mCalibration for %s estimator\e[0;00m -> Starting...\n", estimator.Data());
401+
auto CalibMC = [&](TString estimator){
402+
LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> Starting...\n", estimator.Data());
398403

399404
std::vector<double> percentile_center(nCentBins);
400405
std::vector<double> epercentile_center(nCentBins);
401406

402-
// Histograms
403-
TH2D* h2dMultVsCent_Data = getHist<TH2D>(Form("h2dMultVsCent%s_Data", estimator.Data()));
404-
TH2D* h2dMultRecoVsMultGen_MC = getHist<TH2D>(Form("hMultEta05VsGenMult%s", estimator.Data()));
407+
//Histograms
408+
TH2D *h2dMultVsCent_Data = getHist<TH2D>(Form("h2dMultVsCent%s_Data", estimator.Data()));
409+
TH2D *h2dMultRecoVsMultGen_MC = getHist<TH2D>(Form("hMultEta05VsGenMult%s", estimator.Data()));
405410
if (!h2dMultVsCent_Data || !h2dMultRecoVsMultGen_MC) {
406411
return (TH1D*)0x0;
407412
}
@@ -413,12 +418,12 @@ struct BuilderModule {
413418
TProfile* hPVData = hCalibPVData[idx].get();
414419
TProfile* hPVMC = hCalibPVMC[idx].get();
415420

416-
TH1D* h1dCalib = h2dMultRecoVsMultGen_MC->ProjectionX(Form("h1d%s", estimator.Data()), 1, h2dMultRecoVsMultGen_MC->GetNbinsX());
421+
TH1D *h1dCalib = h2dMultRecoVsMultGen_MC->ProjectionX(Form("h1d%s", estimator.Data()), 1, h2dMultRecoVsMultGen_MC->GetNbinsX());
417422
h1dCalib->Reset();
418423
h1dCalib->SetTitle(Form("%s calibration object", estimator.Data()));
419424
h1dCalib->GetXaxis()->SetTitle(Form("#it{N}_{%s, gen.}", estimator.Data()));
420425
h1dCalib->GetYaxis()->SetTitle(Form("%s percentile (%%)", estimator.Data()));
421-
426+
422427
// NOTE: candidate-by-candidate assignment (assignCentralityPerCandidate) does not go through this
423428
// mean-matching path at all. It samples directly from the "hGenMultEta05VsCentrality<estimator>"
424429
// joint histogram (x = reco centrality of the matched MC collision, y = generated mult |eta|<0.5)
@@ -427,16 +432,16 @@ struct BuilderModule {
427432
if (reverse) {
428433
for (int i = 0; i < nCentBins; i++) {
429434
int irev = i;
430-
percentile_center[i] = (centralityBins[irev] + centralityBins[irev + 1]) / 2;
431-
epercentile_center[i] = (centralityBins[irev] - centralityBins[irev + 1]) / 2;
435+
percentile_center[i] = (centralityBins[irev] + centralityBins[irev+1]) / 2;
436+
epercentile_center[i] = (centralityBins[irev] - centralityBins[irev+1]) / 2;
432437
}
433438

434439
int startBinMc = h1dCalib->GetNbinsX();
435440
for (int i = 0; i < nCentBins; i++) { // Loop over centrality bins
436441
// start from the end (from the high multiplicity collisions)
437442
int irev = i;
438-
TH1D* projData = h2dMultVsCent_Data->ProjectionY(Form("projData_%d", i), h2dMultVsCent_Data->GetXaxis()->FindBin(centralityBins[irev] + 1e-5),
439-
h2dMultVsCent_Data->GetXaxis()->FindBin(centralityBins[irev + 1] - 1e-5));
443+
TH1D* projData = h2dMultVsCent_Data->ProjectionY(Form("projData_%d", i), h2dMultVsCent_Data->GetXaxis()->FindBin(centralityBins[irev]+1e-5),
444+
h2dMultVsCent_Data->GetXaxis()->FindBin(centralityBins[irev+1]-1e-5));
440445

441446
double meanMult_Data = projData->GetMean();
442447

@@ -462,7 +467,7 @@ struct BuilderModule {
462467
endBinMc = j;
463468
}
464469
}
465-
if (i == nCentBins - 1) {
470+
if (i == nCentBins-1) {
466471
endBinMc = 1;
467472
}
468473

@@ -487,30 +492,30 @@ struct BuilderModule {
487492
}
488493
}
489494

490-
printf("\e[1;31mCalibration for %s estimator\e[0;00m -> Data centrality bin \e[1;31m%g-%g%%\e[0;00m\n", estimator.Data(), centralityBins[irev], centralityBins[irev + 1]);
491-
printf("\e[1;31mCalibration for %s estimator\e[0;00m -> MC multiplicity range %g-%g\n", estimator.Data(), projMC->GetBinLowEdge(endBinMc), projMC->GetBinLowEdge(startBinMc));
492-
printf("\e[1;31mCalibration for %s estimator\e[0;00m -> <PV> data = %.4f Vs <PV> MC = %.4f (MC/Data = %.4f%%)\n", estimator.Data(), meanMult_Data, meanMult_MC, (meanMult_MC - meanMult_Data) * 100 / meanMult_Data);
493-
printf("\e[1;31mCalibration for %s estimator\e[0;00m -> N entries Data = %g Vs N entries MC = %g\n", estimator.Data(), projData->Integral(), projMC->Integral());
494-
printf("\n");
495+
LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> Data centrality bin \e[1;31m%g-%g%%\e[0;00m\n", estimator.Data(), centralityBins[irev], centralityBins[irev+1]);
496+
LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> MC multiplicity range %g-%g\n", estimator.Data(), projMC->GetBinLowEdge(endBinMc), projMC->GetBinLowEdge(startBinMc));
497+
LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> <PV> data = %.4f Vs <PV> MC = %.4f (MC/Data = %.4f%%)\n", estimator.Data(), meanMult_Data, meanMult_MC, (meanMult_MC-meanMult_Data) * 100 / meanMult_Data);
498+
LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> N entries Data = %g Vs N entries MC = %g\n", estimator.Data(), projData->Integral(), projMC->Integral());
499+
LOGF(info, "\n");
495500
for (int ibin = 1; ibin <= h1dCalib->GetNbinsX(); ibin++) {
496-
if (ibin <= startBinMc && ibin >= endBinMc) {
497-
h1dCalib->SetBinContent(ibin, percentile_center[i]);
501+
if (ibin <= startBinMc && ibin >= endBinMc){
502+
h1dCalib->SetBinContent(ibin, percentile_center[i]);
498503
}
499504
}
500505
startBinMc = endBinMc;
501506
} // End loop over centrality bins
502507
} else {
503508
for (int i = 0; i < nCentBins; i++) {
504509
int irev = nCentBins - i;
505-
percentile_center[i] = (centralityBins[irev - 1] + centralityBins[irev]) / 2;
506-
epercentile_center[i] = (centralityBins[irev] - centralityBins[irev - 1]) / 2;
510+
percentile_center[i] = (centralityBins[irev-1] + centralityBins[irev]) / 2;
511+
epercentile_center[i] = (centralityBins[irev] - centralityBins[irev-1]) / 2;
507512
}
508513
int startBinMc = 1;
509514
for (int i = 0; i < nCentBins; i++) { // Loop over centrality bins
510515
// start from the end (from the low multiplicity collisions)
511516
int irev = nCentBins - i;
512-
TH1D* projData = h2dMultVsCent_Data->ProjectionY(Form("projData_%d", i), h2dMultVsCent_Data->GetXaxis()->FindBin(centralityBins[irev - 1] + 1e-5),
513-
h2dMultVsCent_Data->GetXaxis()->FindBin(centralityBins[irev] - 1e-5));
517+
TH1D* projData = h2dMultVsCent_Data->ProjectionY(Form("projData_%d", i), h2dMultVsCent_Data->GetXaxis()->FindBin(centralityBins[irev-1]+1e-5),
518+
h2dMultVsCent_Data->GetXaxis()->FindBin(centralityBins[irev]-1e-5));
514519

515520
double meanMult_Data = projData->GetMean();
516521

@@ -536,8 +541,8 @@ struct BuilderModule {
536541
endBinMc = j;
537542
}
538543
}
539-
if (i == nCentBins - 1) {
540-
endBinMc = h2dMultRecoVsMultGen_MC->GetNbinsX();
544+
if (i == nCentBins-1) {
545+
endBinMc = h2dMultRecoVsMultGen_MC->GetNbinsX();
541546
}
542547

543548
TH1D* projMC = h2dMultRecoVsMultGen_MC->ProjectionY(Form("projMC_%d", i), startBinMc, endBinMc);
@@ -561,21 +566,21 @@ struct BuilderModule {
561566
}
562567
}
563568

564-
printf("\e[1;31mCalibration for %s estimator\e[0;00m -> Data centrality bin \e[1;31m%g-%g%%\e[0;00m\n", estimator.Data(), centralityBins[irev - 1], centralityBins[irev]);
565-
printf("\e[1;31mCalibration for %s estimator\e[0;00m -> MC multiplicity range %g-%g\n", estimator.Data(), projMC->GetBinLowEdge(startBinMc), projMC->GetBinLowEdge(endBinMc + 1));
566-
printf("\e[1;31mCalibration for %s estimator\e[0;00m -> <PV> data = %.4f Vs <PV> MC = %.4f\n", estimator.Data(), meanMult_Data, meanMult_MC);
567-
printf("\e[1;31mCalibration for %s estimator\e[0;00m -> N entries Data = %g Vs N entries MC = %g\n", estimator.Data(), projData->Integral(), projMC->Integral());
569+
LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> Data centrality bin \e[1;31m%g-%g%%\e[0;00m\n", estimator.Data(), centralityBins[irev-1], centralityBins[irev]);
570+
LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> MC multiplicity range %g-%g\n", estimator.Data(), projMC->GetBinLowEdge(startBinMc), projMC->GetBinLowEdge(endBinMc+1));
571+
LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> <PV> data = %.4f Vs <PV> MC = %.4f\n", estimator.Data(), meanMult_Data, meanMult_MC);
572+
LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> N entries Data = %g Vs N entries MC = %g\n", estimator.Data(), projData->Integral(), projMC->Integral());
568573
for (int ibin = 1; ibin <= h1dCalib->GetNbinsX(); ibin++) {
569-
if (ibin <= endBinMc && ibin >= startBinMc) {
570-
h1dCalib->SetBinContent(ibin, percentile_center[i]);
574+
if (ibin <= endBinMc && ibin >= startBinMc){
575+
h1dCalib->SetBinContent(ibin, percentile_center[i]);
571576
}
572577
}
573578

574-
startBinMc = endBinMc + 1;
579+
startBinMc = endBinMc+1;
575580
} // End loop over centrality bins
576581
}
577582

578-
printf("\e[1;31mCalibration for %s estimator\e[0;00m -> Done!\n", estimator.Data());
583+
LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> Done!\n", estimator.Data());
579584

580585
return h1dCalib;
581586
};
@@ -657,8 +662,8 @@ struct BuilderModule {
657662
// which has no systematic effect on the sampled distribution.
658663
ULong64_t computeSamplingSeed(int64_t collisionIndex, int tableIdx) const
659664
{
660-
static constexpr ULong64_t kCollisionStride = 131ull; // > kNestimators, so tableIdx can't alias into the collision term
661-
static constexpr ULong64_t kRunStride = 1000003ull; // prime, well above kCollisionStride * (typical collisions per DF)
665+
static constexpr ULong64_t kCollisionStride = 131ull; // > kNestimators, so tableIdx can't alias into the collision term
666+
static constexpr ULong64_t kRunStride = 1000003ull; // prime, well above kCollisionStride * (typical collisions per DF)
662667
return baseOpts.centralitySamplingSeed.value + kRunStride * static_cast<ULong64_t>(mRunNumber) +
663668
kCollisionStride * static_cast<ULong64_t>(collisionIndex) + static_cast<ULong64_t>(tableIdx);
664669
}
@@ -712,15 +717,15 @@ struct BuilderModule {
712717
return percentile;
713718
};
714719

715-
populateTable(products.centFV0A, h1dFV0A, h2dCentVsGenMult[kFV0A], nFV0A, std::integral_constant<int, kFV0A>{});
716-
populateTable(products.centFT0M, h1dFT0M, h2dCentVsGenMult[kFT0M], nFT0M, std::integral_constant<int, kFT0M>{});
717-
populateTable(products.centFT0A, h1dFT0A, h2dCentVsGenMult[kFT0A], nFT0A, std::integral_constant<int, kFT0A>{});
718-
populateTable(products.centFT0C, h1dFT0C, h2dCentVsGenMult[kFT0C], nFT0C, std::integral_constant<int, kFT0C>{});
719-
populateTable(products.centFT0CVariant1, h1dFT0CVariant1, h2dCentVsGenMult[kFT0CVariant1], nFT0C, std::integral_constant<int, kFT0CVariant1>{});
720-
populateTable(products.centFT0CVariant2, h1dFT0CVariant2, h2dCentVsGenMult[kFT0CVariant2], nFT0C, std::integral_constant<int, kFT0CVariant2>{});
721-
populateTable(products.centFDDM, h1dFDDM, h2dCentVsGenMult[kFDDM], nFDDM, std::integral_constant<int, kFDDM>{});
722-
populateTable(products.centNTPV, h1dNTPV, h2dCentVsGenMult[kNTPV], nGlobal, std::integral_constant<int, kNTPV>{});
723-
populateTable(products.centNGlobal, h1dNGlobal, h2dCentVsGenMult[kNGlobal], nGlobal, std::integral_constant<int, kNGlobal>{});
720+
populateTable(products.centFV0A, h1dFV0A, h2dCentVsGenMult[kFV0A], nFV0A, std::integral_constant<int, kFV0A>{});
721+
populateTable(products.centFT0M, h1dFT0M, h2dCentVsGenMult[kFT0M], nFT0M, std::integral_constant<int, kFT0M>{});
722+
populateTable(products.centFT0A, h1dFT0A, h2dCentVsGenMult[kFT0A], nFT0A, std::integral_constant<int, kFT0A>{});
723+
populateTable(products.centFT0C, h1dFT0C, h2dCentVsGenMult[kFT0C], nFT0C, std::integral_constant<int, kFT0C>{});
724+
populateTable(products.centFT0CVariant1, h1dFT0CVariant1, h2dCentVsGenMult[kFT0CVariant1], nFT0C, std::integral_constant<int, kFT0CVariant1>{});
725+
populateTable(products.centFT0CVariant2, h1dFT0CVariant2, h2dCentVsGenMult[kFT0CVariant2], nFT0C, std::integral_constant<int, kFT0CVariant2>{});
726+
populateTable(products.centFDDM, h1dFDDM, h2dCentVsGenMult[kFDDM], nFDDM, std::integral_constant<int, kFDDM>{});
727+
populateTable(products.centNTPV, h1dNTPV, h2dCentVsGenMult[kNTPV], nGlobal, std::integral_constant<int, kNTPV>{});
728+
populateTable(products.centNGlobal, h1dNGlobal, h2dCentVsGenMult[kNGlobal], nGlobal, std::integral_constant<int, kNGlobal>{});
724729
// populateTable(products.centMFT, h1dMFT, h2dCentVsGenMult[kMFT], nMFT, std::integral_constant<int, kMFT>{}); // to be added later
725730

726731
mcCollisionCounter++;

0 commit comments

Comments
 (0)