Skip to content

Commit 91b9298

Browse files
lietavalietava
authored
[PWGLF] Omegahm: thnsparse for nult studies and for response of dN/deta (AliceO2Group#14110)
1 parent 34a77ad commit 91b9298

File tree

1 file changed

+153
-56
lines changed

1 file changed

+153
-56
lines changed

‎PWGLF/Tasks/Strangeness/nonPromptCascade.cxx‎

Lines changed: 153 additions & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -37,7 +37,9 @@
3737

3838
#include "Math/Vector4D.h"
3939
#include "TDatabasePDG.h"
40+
#include "THnSparse.h"
4041
#include "TParticlePDG.h"
42+
#include "TTree.h"
4143

4244
#include <cmath>
4345
#include <memory>
@@ -185,6 +187,7 @@ struct NonPromptCascadeTask {
185187
using CollisionCandidatesRun3 = soa::Join<aod::Collisions, aod::EvSels, aod::FT0Mults, aod::FV0Mults, aod::CentFT0Cs, aod::CentFV0As, aod::CentFT0Ms, aod::MultsGlobal>;
186188
using CollisionCandidatesRun3MC = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::FT0Mults, aod::FV0Mults, aod::CentFT0Cs, aod::CentFV0As, aod::CentFT0Ms, aod::MultsGlobal>;
187189
using CollisionsWithLabel = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::MultsGlobal>;
190+
using TracksWithLabel = soa::Join<aod::Tracks, aod::McTrackLabels>;
188191

189192
Preslice<TracksExtData> perCollision = aod::track::collisionId;
190193
Preslice<TracksExtMC> perCollisionMC = aod::track::collisionId;
@@ -212,8 +215,7 @@ struct NonPromptCascadeTask {
212215

213216
Configurable<float> cfgMaxMult{"cfgMaxMult", 8000.f, "Upper range of multiplicty histo"};
214217
Configurable<float> cfgMaxMultFV0{"cfgMaxMultFV0", 10000.f, "Upper range of multiplicty FV0 histo"};
215-
Configurable<float> cfgMinMult{"cfgMinMult", 3000.f, "Lower range of FT0M histo in zoomed histo"};
216-
Configurable<float> cfgMaxCent{"cfgMaxCent", 8.0025f, "Upper range of FT0M histo"};
218+
Configurable<std::string> cfgPtEdgesdNdeta{"ptEdges", "0,0.2,0.4,0.6,0.8,1,1.2,1.6,2.0,2.4,2.8,3.2,3.6,4,4.5,5,5.5,6,7,8,10", "Pt bin edges (comma-separated)"};
217219

218220
Zorro mZorro;
219221
OutputObj<ZorroSummary> mZorroSummary{"ZorroSummary"};
@@ -225,23 +227,11 @@ struct NonPromptCascadeTask {
225227
o2::vertexing::DCAFitterN<2> mDCAFitter;
226228
std::array<int, 2> mProcessCounter = {0, 0}; // {Tracked, All}
227229
std::map<uint64_t, uint32_t> mToiMap;
228-
std::unordered_map<std::string, std::shared_ptr<TH2>> mHistsPerRunMultVsCent;
229-
std::unordered_map<std::string, std::shared_ptr<TH2>> mHistsPerRunMultVsCentZoom;
230-
std::unordered_map<std::string, std::shared_ptr<TH2>> mHistsPerRunNtracktVsCent;
231-
std::unordered_map<std::string, std::shared_ptr<TH2>> mHistsPerRunNtracktVsCentZoom;
232-
233-
int nBinsMult = cfgMaxMult;
234-
int nBinsMultFV0 = cfgMaxMultFV0;
235-
int nBinsMultZoom = cfgMaxMult - cfgMinMult;
236-
int nBinsCentZoom = (cfgMaxCent + 0.0025) / 0.005;
237-
238-
AxisSpec multAxis = {nBinsMult, 0, cfgMaxMult, "Multiplicity FT0M"};
239-
AxisSpec multAxisFV0 = {nBinsMultFV0, 0, cfgMaxMultFV0, "Multiplicity FT0M"};
240-
AxisSpec centAxis = {101, -0.025, 101.025, "Centrality"};
241-
AxisSpec centAxisZoom = {nBinsCentZoom, -0.0025, cfgMaxCent, "Centrality"};
242-
AxisSpec multAxisZoom = {nBinsMultZoom, cfgMinMult, cfgMaxMult, "Multiplicity FT0M"};
243-
AxisSpec nTracksAxis = {100, 0., 100., "NTracksGlobal"};
244-
AxisSpec nTracksAxisMC = {100, 0., 100., "NTracksMC"};
230+
//
231+
HistogramRegistry mRegistryMults{"Multhistos"};
232+
HistogramRegistry mRegistrydNdeta{"dNdetahistos"};
233+
234+
//
245235

246236
void initCCDB(aod::BCsWithTimestamps::iterator const& bc)
247237
{
@@ -262,7 +252,7 @@ struct NonPromptCascadeTask {
262252
}
263253
}
264254

265-
void init(InitContext const&)
255+
void init(InitContext& context)
266256
{
267257
mZorroSummary.setObject(mZorro.getZorroSummary());
268258
mCCDB->setURL(ccdbUrl);
@@ -283,18 +273,77 @@ struct NonPromptCascadeTask {
283273
std::array<std::string, 7> cutsNames{"# candidates", "hasTOF", "nClusTPC", "nSigmaTPCbach", "nSigmaTPCprotontrack", "nSigmaTPCpiontrack", "cosPA"};
284274
auto cutsOmega{std::get<std::shared_ptr<TH2>>(mRegistry.add("h_PIDcutsOmega", ";;Invariant mass (GeV/#it{c}^{2})", HistType::kTH2D, {{cutsNames.size(), -0.5, -0.5 + cutsNames.size()}, {125, 1.650, 1.700}}))};
285275
auto cutsXi{std::get<std::shared_ptr<TH2>>(mRegistry.add("h_PIDcutsXi", ";;Invariant mass (GeV/#it{c}^{2})", HistType::kTH2D, {{6, -0.5, 5.5}, {125, 1.296, 1.346}}))};
286-
mRegistry.add("hMultVsCent", "hMultVsCent", HistType::kTH2F, {centAxis, multAxis});
287-
mRegistry.add("hMultVsCentZoom", "hMultVsCentZoom", HistType::kTH2F, {centAxisZoom, multAxisZoom});
288-
mRegistry.add("hNTracksVsCent", "hNTracksVsCent", HistType::kTH2F, {centAxis, nTracksAxis});
289-
mRegistry.add("hNTracksVsCentZoom", "hNTracksVsCentZoom", HistType::kTH2F, {centAxisZoom, nTracksAxis});
290-
mRegistry.add("hMultFV0VshNTracks", "hMultFV0VshNTracks", HistType::kTH2F, {nTracksAxis, multAxisFV0});
291-
mRegistry.add("hNTracksVsCentFV0A", "hNTracksVsCentFV0A", HistType::kTH2F, {nTracksAxis, centAxis});
292-
mRegistry.add("hNTracksMCVsTracksReco", "hNTracksMCVsTracksReco", HistType::kTH2F, {nTracksAxisMC, nTracksAxis});
293-
mRegistry.add("hNTracksMCNotInReco", "hNTracksMCNotInReco", HistType::kTH1F, {nTracksAxisMC});
294276
for (size_t iBin{0}; iBin < cutsNames.size(); ++iBin) {
295277
cutsOmega->GetYaxis()->SetBinLabel(iBin + 1, cutsNames[iBin].c_str());
296278
cutsXi->GetYaxis()->SetBinLabel(iBin + 1, cutsNames[iBin].c_str());
297279
}
280+
//
281+
int nBinsMult = cfgMaxMult;
282+
int nBinsMultFV0 = cfgMaxMultFV0;
283+
284+
AxisSpec multAxis = {nBinsMult, 0, cfgMaxMult, "Multiplicity FT0M"};
285+
AxisSpec multAxisFV0 = {nBinsMultFV0, 0, cfgMaxMultFV0, "Multiplicity FT0M"};
286+
AxisSpec nTracksAxis = {100, 0., 100., "NTracksGlobal"};
287+
AxisSpec nTracksAxisMC = {100, 0., 100., "NTracksMC"};
288+
std::vector<double> centBinning;
289+
std::vector<double> multBinning;
290+
std::vector<double> trackBinning;
291+
std::vector<double> runsBinning;
292+
double cent = -0.0;
293+
while (cent < 100) {
294+
if (cent < 0.1) {
295+
centBinning.push_back(cent);
296+
cent += 0.01;
297+
} else if (cent < 1) {
298+
centBinning.push_back(cent);
299+
cent += 0.1;
300+
} else {
301+
centBinning.push_back(cent);
302+
cent += 1;
303+
}
304+
}
305+
double ntracks = 0;
306+
while (ntracks < 100) {
307+
trackBinning.push_back(ntracks);
308+
ntracks++;
309+
}
310+
double run = 550367 - 0.5;
311+
for (int i = 0; i < 9000; i++) {
312+
runsBinning.push_back(run);
313+
run++;
314+
}
315+
AxisSpec centAxis{centBinning, "Centrality (%)"};
316+
// AxisSpec multAxis{multBinning, "Multiplicity"};
317+
AxisSpec trackAxisMC{trackBinning, "NTracks MC"};
318+
AxisSpec trackAxis{trackBinning, "NTracks Global Reco"};
319+
AxisSpec runsAxis{runsBinning, "Run Number"};
320+
321+
mRegistryMults.add("hCentMultsRuns", "hCentMultsRuns", HistType::kTHnSparseF, {centAxis, multAxis, centAxis, multAxisFV0, nTracksAxis, runsAxis});
322+
//
323+
// dN/deta
324+
//
325+
bool runMCdNdeta = context.options().get<bool>("processdNdetaMC");
326+
// std::cout << "runMCdNdeta: " << runMCdNdeta << std::endl;
327+
if (runMCdNdeta) {
328+
std::vector<double> ptBins;
329+
std::vector<std::string> tokens = o2::utils::Str::tokenize(cfgPtEdgesdNdeta, ',');
330+
for (auto const& pts : tokens) {
331+
double pt = 0;
332+
try {
333+
pt = std::stof(pts);
334+
} catch (...) {
335+
LOG(error) << "Wrong cfgPtEdgesdNdeta string:" << cfgPtEdgesdNdeta << std::endl;
336+
}
337+
ptBins.push_back(pt);
338+
}
339+
AxisSpec ptAxisMC{ptBins, "pT MC"};
340+
AxisSpec ptAxisReco{ptBins, "pT Reco"};
341+
342+
// multMeasured, multMC, ptMeasured, ptMC
343+
mRegistrydNdeta.add("hdNdetaRM/hdNdetaRM", "hdNdetaRM", HistType::kTHnSparseF, {nTracksAxisMC, nTracksAxis, ptAxisMC, ptAxisReco});
344+
mRegistrydNdeta.add("hdNdetaRM/hdNdetaRMNotInRecoCol", "hdNdetaRMNotInRecoCol", HistType::kTHnSparseF, {nTracksAxisMC, ptAxisMC});
345+
mRegistrydNdeta.add("hdNdetaRM/hdNdetaRMNotInRecoTrk", "hdNdetaRMNotInRecoTrk", HistType::kTHnSparseF, {nTracksAxisMC, ptAxisMC});
346+
}
298347
}
299348

300349
template <typename CollisionType, typename TrackType>
@@ -361,27 +410,13 @@ struct NonPromptCascadeTask {
361410
{
362411
// std::cout << "Filling mult histos" << std::endl;
363412
for (const auto& coll : collisions) {
364-
std::string histNameMvC = "mult/hMultVsCent_run" + std::to_string(mRunNumber);
365-
std::string histNameMvCZ = "mult/hMultVsCentZoom_run" + std::to_string(mRunNumber);
366-
std::string histNameTvC = "mult/hNTracksVsCent_run" + std::to_string(mRunNumber);
367-
std::string histNameTvCZ = "mult/hNTracksVsCentZoom_run" + std::to_string(mRunNumber);
368-
if (!mHistsPerRunMultVsCent.contains(histNameMvC)) {
369-
mHistsPerRunMultVsCent[histNameMvC] = std::get<std::shared_ptr<TH2>>(mRegistry.add(histNameMvC.c_str(), histNameMvC.c_str(), HistType::kTH2F, {centAxis, multAxis}));
370-
mHistsPerRunMultVsCentZoom[histNameMvCZ] = std::get<std::shared_ptr<TH2>>(mRegistry.add(histNameMvCZ.c_str(), histNameMvCZ.c_str(), HistType::kTH2F, {centAxisZoom, multAxisZoom}));
371-
mHistsPerRunNtracktVsCent[histNameTvC] = std::get<std::shared_ptr<TH2>>(mRegistry.add(histNameTvC.c_str(), histNameTvC.c_str(), HistType::kTH2F, {centAxis, nTracksAxis}));
372-
mHistsPerRunNtracktVsCentZoom[histNameTvCZ] = std::get<std::shared_ptr<TH2>>(mRegistry.add(histNameTvCZ.c_str(), histNameTvCZ.c_str(), HistType::kTH2F, {centAxisZoom, nTracksAxis}));
373-
}
374-
mHistsPerRunMultVsCent[histNameMvC]->Fill(coll.centFT0M(), coll.multFT0M());
375-
mHistsPerRunMultVsCentZoom[histNameMvCZ]->Fill(coll.centFT0M(), coll.multFT0M());
376-
mHistsPerRunNtracktVsCent[histNameTvC]->Fill(coll.centFT0M(), coll.multNTracksGlobal());
377-
mHistsPerRunNtracktVsCentZoom[histNameTvCZ]->Fill(coll.centFT0M(), coll.multNTracksGlobal());
378-
// run integrated histos
379-
mRegistry.fill(HIST("hMultVsCent"), coll.centFT0M(), coll.multFT0M());
380-
mRegistry.fill(HIST("hMultVsCentZoom"), coll.centFT0M(), coll.multFT0M());
381-
mRegistry.fill(HIST("hNTracksVsCent"), coll.centFT0M(), (float)coll.multNTracksGlobal());
382-
mRegistry.fill(HIST("hNTracksVsCentZoom"), coll.centFT0M(), coll.multNTracksGlobal());
383-
mRegistry.fill(HIST("hMultFV0VshNTracks"), coll.multNTracksGlobal(), coll.multFV0A());
384-
mRegistry.fill(HIST("hNTracksVsCentFV0A"), coll.multNTracksGlobal(), coll.centFV0A());
413+
float centFT0M = coll.centFT0M();
414+
float multFT0M = coll.multFT0M();
415+
float centFV0A = coll.centFV0A();
416+
float multFV0A = coll.multFV0A();
417+
float multNTracks = coll.multNTracksGlobal();
418+
float run = mRunNumber;
419+
mRegistryMults.fill(HIST("hCentMultsRuns"), centFT0M, multFT0M, centFV0A, multFV0A, multNTracks, run);
385420
}
386421
};
387422

@@ -758,7 +793,7 @@ struct NonPromptCascadeTask {
758793
}
759794
PROCESS_SWITCH(NonPromptCascadeTask, processCascadesData, "process cascades: Data analysis", false);
760795

761-
int getMCMult(aod::McParticles const& mcParticles, int mcCollId)
796+
int getMCMult(aod::McParticles const& mcParticles, int mcCollId, std::vector<float>& ptMCvec)
762797
{
763798
int mult = 0;
764799
for (auto const& mcp : mcParticles) {
@@ -776,35 +811,97 @@ struct NonPromptCascadeTask {
776811
accept = accept && (q != 0);
777812
if (accept) {
778813
++mult;
814+
ptMCvec.push_back(mcp.pt());
779815
}
780816
}
781817
}
782818
return mult;
783819
}
784-
void processNegMC(CollisionsWithLabel const& colls, aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles)
820+
void processdNdetaMC(CollisionsWithLabel const& colls, aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles, TracksWithLabel const& tracks)
785821
{
786822
// std::cout << "ProcNegMC" << std::endl;
823+
// Map: collision index -> reco multiplicity
824+
std::vector<int> recoMult(colls.size(), 0);
825+
for (auto const& trk : tracks) {
826+
int collId = trk.collisionId();
827+
// Here you can impose same track cuts as for MC (eta, primaries, etc.)
828+
if (std::abs(trk.eta()) > 0.5f) {
829+
continue;
830+
}
831+
++recoMult[collId];
832+
}
833+
std::vector<int> isReco(mcParticles.size(), 0);
834+
std::vector<int> isRecoMult(mcParticles.size(), 0);
787835
std::vector<int> mcReconstructed(mcCollisions.size(), 0);
836+
// std::cout << " reco cols with mc:" << colls.size() << " tracks:" << tracks.size() << " mccols:" << mcCollisions.size() << "mcParticles:" << mcParticles.size() << std::endl;
788837
for (auto const& col : colls) {
789838
int mcCollId = col.mcCollisionId(); // col.template mcCollision_as<aod::McCollisions>();
839+
int collId = col.globalIndex();
790840
// auto mc = col.mcCollision();
791841
// int mcId = mc.globalIndex();
792842
// std::cout << "globalIndex:" << mcId << " colID:" << mcCollId << std::endl;
793-
int mult = getMCMult(mcParticles, mcCollId);
843+
std::vector<float> mcptvec;
844+
float mult = getMCMult(mcParticles, mcCollId, mcptvec);
794845
mcReconstructed[mcCollId] = 1;
795-
mRegistry.fill(HIST("hNTracksMCVsTracksReco"), mult, col.multNTracksGlobal());
846+
for (auto const& trk : tracks) {
847+
int mcPid = trk.mcParticleId(); // depends on your label table
848+
if (mcPid >= 0 && mcPid < (int)mcParticles.size()) {
849+
isReco[mcPid] = 1;
850+
isRecoMult[mcPid] = mult;
851+
} else {
852+
continue;
853+
}
854+
if (trk.collisionId() != collId) {
855+
continue; // different event
856+
}
857+
auto mcPar = mcParticles.rawIteratorAt(mcPid);
858+
// Apply same acceptance as in MC multiplicity
859+
if (!mcPar.isPhysicalPrimary()) {
860+
continue;
861+
}
862+
if (std::abs(mcPar.eta()) > 0.5f) {
863+
continue;
864+
}
865+
int q = 0;
866+
auto pdgEntry = TDatabasePDG::Instance()->GetParticle(mcPar.pdgCode());
867+
if (pdgEntry) {
868+
q = int(std::round(pdgEntry->Charge() / 3.0));
869+
}
870+
if (q == 0) {
871+
continue;
872+
}
873+
// float multReco = recoMult[collId];
874+
float multReco = col.multNTracksGlobal();
875+
float ptReco = trk.pt();
876+
float ptMC = mcPar.pt();
877+
mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRM"), mult, multReco, ptMC, ptReco);
878+
}
879+
880+
// mRegistry.fill(HIST("hNTracksMCVsTracksReco"), mult, col.multNTracksGlobal());
796881
}
882+
// count mc particles with no reco tracks
883+
for (auto const& mcp : mcParticles) {
884+
int mcPidG = mcp.globalIndex();
885+
// std::cout << "mcPidG:" << mcPidG << std::endl;
886+
if (!isReco[mcPidG]) {
887+
mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoTrk"), isRecoMult[mcPidG], mcp.pt());
888+
}
889+
}
890+
// count unreconstructed mc collisions
797891
for (auto const& mc : mcCollisions) {
798892
int gindex = mc.globalIndex();
799893
// std::cout << "mc globalIndex:" << gindex << std::endl;
800894
if (!mcReconstructed[gindex]) {
801-
int mult = getMCMult(mcParticles, gindex);
895+
std::vector<float> mcptvec;
896+
int mult = getMCMult(mcParticles, gindex, mcptvec);
802897
// std::cout << "===> unreconstructed:" << mult << std::endl;
803-
mRegistry.fill(HIST("hNTracksMCNotInReco"), mult);
898+
for (auto const& pt : mcptvec) {
899+
mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoCol"), mult, pt);
900+
}
804901
}
805902
}
806903
}
807-
PROCESS_SWITCH(NonPromptCascadeTask, processNegMC, "process mc", false);
904+
PROCESS_SWITCH(NonPromptCascadeTask, processdNdetaMC, "process mc dN/deta", false);
808905
};
809906

810907
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)