Skip to content

Commit 8f4460e

Browse files
lietavalietava
committed
dev: use of thnsparse for mults and response
1 parent 8019d0e commit 8f4460e

File tree

1 file changed

+161
-62
lines changed

1 file changed

+161
-62
lines changed

‎PWGLF/Tasks/Strangeness/nonPromptCascade.cxx‎

Lines changed: 161 additions & 62 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,8 @@
3838
#include "Math/Vector4D.h"
3939
#include "TDatabasePDG.h"
4040
#include "TParticlePDG.h"
41+
#include "TTree.h"
42+
#include "THnSparse.h"
4143

4244
#include <cmath>
4345
#include <memory>
@@ -185,6 +187,8 @@ 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>;
191+
188192

189193
Preslice<TracksExtData> perCollision = aod::track::collisionId;
190194
Preslice<TracksExtMC> perCollisionMC = aod::track::collisionId;
@@ -212,8 +216,7 @@ struct NonPromptCascadeTask {
212216

213217
Configurable<float> cfgMaxMult{"cfgMaxMult", 8000.f, "Upper range of multiplicty histo"};
214218
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"};
219+
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)"};
217220

218221
Zorro mZorro;
219222
OutputObj<ZorroSummary> mZorroSummary{"ZorroSummary"};
@@ -225,23 +228,11 @@ struct NonPromptCascadeTask {
225228
o2::vertexing::DCAFitterN<2> mDCAFitter;
226229
std::array<int, 2> mProcessCounter = {0, 0}; // {Tracked, All}
227230
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"};
231+
//
232+
HistogramRegistry mRegistryMults{"Multhistos"};
233+
HistogramRegistry mRegistrydNdeta{"dNdetahistos"};
234+
235+
//
245236

246237
void initCCDB(aod::BCsWithTimestamps::iterator const& bc)
247238
{
@@ -262,7 +253,7 @@ struct NonPromptCascadeTask {
262253
}
263254
}
264255

265-
void init(InitContext const&)
256+
void init(InitContext & context)
266257
{
267258
mZorroSummary.setObject(mZorro.getZorroSummary());
268259
mCCDB->setURL(ccdbUrl);
@@ -283,18 +274,77 @@ struct NonPromptCascadeTask {
283274
std::array<std::string, 7> cutsNames{"# candidates", "hasTOF", "nClusTPC", "nSigmaTPCbach", "nSigmaTPCprotontrack", "nSigmaTPCpiontrack", "cosPA"};
284275
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}}))};
285276
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});
294277
for (size_t iBin{0}; iBin < cutsNames.size(); ++iBin) {
295278
cutsOmega->GetYaxis()->SetBinLabel(iBin + 1, cutsNames[iBin].c_str());
296279
cutsXi->GetYaxis()->SetBinLabel(iBin + 1, cutsNames[iBin].c_str());
297280
}
281+
//
282+
int nBinsMult = cfgMaxMult;
283+
int nBinsMultFV0 = cfgMaxMultFV0;
284+
285+
AxisSpec multAxis = {nBinsMult, 0, cfgMaxMult, "Multiplicity FT0M"};
286+
AxisSpec multAxisFV0 = {nBinsMultFV0, 0, cfgMaxMultFV0, "Multiplicity FT0M"};
287+
AxisSpec nTracksAxis = {100, 0., 100., "NTracksGlobal"};
288+
AxisSpec nTracksAxisMC = {100, 0., 100., "NTracksMC"};
289+
std::vector<double> centBinning;
290+
std::vector<double> multBinning;
291+
std::vector<double> trackBinning;
292+
std::vector<double> runsBinning;
293+
double cent = -0.0;
294+
while(cent < 100) {
295+
if(cent < 0.1) {
296+
centBinning.push_back(cent);
297+
cent += 0.01;
298+
} else if(cent < 1) {
299+
centBinning.push_back(cent);
300+
cent += 0.1;
301+
} else {
302+
centBinning.push_back(cent);
303+
cent += 1;
304+
}
305+
}
306+
double ntracks = 0;
307+
while(ntracks < 100) {
308+
trackBinning.push_back(ntracks);
309+
ntracks++;
310+
}
311+
double run = 550367 - 0.5;
312+
for(int i = 0; i < 9000; i++) {
313+
runsBinning.push_back(run);
314+
run++;
315+
}
316+
AxisSpec centAxis{centBinning, "Centrality (%)"};
317+
//AxisSpec multAxis{multBinning, "Multiplicity"};
318+
AxisSpec trackAxisMC{trackBinning, "NTracks MC"};
319+
AxisSpec trackAxis{trackBinning, "NTracks Global Reco"};
320+
AxisSpec runsAxis{runsBinning, "Run Number"};
321+
322+
mRegistryMults.add("hCentMultsRuns", "hCentMultsRuns", HistType::kTHnSparseF, {centAxis, multAxis, centAxis, multAxisFV0, nTracksAxis, runsAxis});
323+
//
324+
// dN/deta
325+
//
326+
bool runMCdNdeta = context.options().get<bool>("processdNdetaMC");
327+
//std::cout << "runMCdNdeta: " << runMCdNdeta << std::endl;
328+
if(runMCdNdeta) {
329+
std::vector<double> ptBins;
330+
std::vector<std::string> tokens = o2::utils::Str::tokenize(cfgPtEdgesdNdeta, ',');
331+
for(auto const& pts: tokens){
332+
double pt = 0;
333+
try{
334+
pt = std::stof(pts);
335+
} catch (...) {
336+
LOG(error) << "Wrong cfgPtEdgesdNdeta string:" << cfgPtEdgesdNdeta << std::endl;
337+
}
338+
ptBins.push_back(pt);
339+
}
340+
AxisSpec ptAxisMC{ptBins, "pT MC"};
341+
AxisSpec ptAxisReco{ptBins, "pT Reco"};
342+
343+
// multMeasured, multMC, ptMeasured, ptMC
344+
mRegistrydNdeta.add("hdNdetaRM/hdNdetaRM", "hdNdetaRM", HistType::kTHnSparseF, {nTracksAxisMC, nTracksAxis, ptAxisMC, ptAxisReco});
345+
mRegistrydNdeta.add("hdNdetaRM/hdNdetaRMNotInRecoCol", "hdNdetaRMNotInRecoCol", HistType::kTHnSparseF, {nTracksAxisMC, ptAxisMC});
346+
mRegistrydNdeta.add("hdNdetaRM/hdNdetaRMNotInRecoTrk", "hdNdetaRMNotInRecoTrk", HistType::kTHnSparseF, {nTracksAxisMC, ptAxisMC});
347+
}
298348
}
299349

300350
template <typename CollisionType, typename TrackType>
@@ -361,27 +411,13 @@ struct NonPromptCascadeTask {
361411
{
362412
// std::cout << "Filling mult histos" << std::endl;
363413
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());
414+
float centFT0M = coll.centFT0M();
415+
float multFT0M = coll.multFT0M();
416+
float centFV0A = coll.centFV0A();
417+
float multFV0A = coll.multFV0A();
418+
float multNTracks = coll.multNTracksGlobal();
419+
float run = mRunNumber;
420+
mRegistryMults.fill(HIST("hCentMultsRuns"), centFT0M, multFT0M, centFV0A, multFV0A, multNTracks, run);
385421
}
386422
};
387423

@@ -758,7 +794,7 @@ struct NonPromptCascadeTask {
758794
}
759795
PROCESS_SWITCH(NonPromptCascadeTask, processCascadesData, "process cascades: Data analysis", false);
760796

761-
int getMCMult(aod::McParticles const& mcParticles, int mcCollId)
797+
int getMCMult(aod::McParticles const& mcParticles, int mcCollId, std::vector<float>& ptMCvec)
762798
{
763799
int mult = 0;
764800
for (auto const& mcp : mcParticles) {
@@ -775,36 +811,99 @@ struct NonPromptCascadeTask {
775811
}
776812
accept = accept && (q != 0);
777813
if (accept) {
778-
++mult;
814+
++mult;
815+
ptMCvec.push_back(mcp.pt());
779816
}
780817
}
781818
}
782819
return mult;
783820
}
784-
void processNegMC(CollisionsWithLabel const& colls, aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles)
821+
void processdNdetaMC(CollisionsWithLabel const& colls, aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles, TracksWithLabel const& tracks)
785822
{
786823
// std::cout << "ProcNegMC" << std::endl;
824+
// Map: collision index -> reco multiplicity
825+
std::vector<int> recoMult(colls.size(), 0);
826+
for (auto const& trk : tracks) {
827+
int collId = trk.collisionId();
828+
// Here you can impose same track cuts as for MC (eta, primaries, etc.)
829+
if (std::abs(trk.eta()) > 0.5f) {
830+
continue;
831+
}
832+
++recoMult[collId];
833+
}
834+
std::vector<int> isReco(mcParticles.size(), 0);
835+
std::vector<int> isRecoMult(mcParticles.size(), 0);
787836
std::vector<int> mcReconstructed(mcCollisions.size(), 0);
837+
//std::cout << " reco cols with mc:" << colls.size() << " tracks:" << tracks.size() << " mccols:" << mcCollisions.size() << "mcParticles:" << mcParticles.size() << std::endl;
788838
for (auto const& col : colls) {
789839
int mcCollId = col.mcCollisionId(); // col.template mcCollision_as<aod::McCollisions>();
790-
// auto mc = col.mcCollision();
791-
// int mcId = mc.globalIndex();
792-
// std::cout << "globalIndex:" << mcId << " colID:" << mcCollId << std::endl;
793-
int mult = getMCMult(mcParticles, mcCollId);
840+
int collId = col.globalIndex();
841+
//auto mc = col.mcCollision();
842+
//int mcId = mc.globalIndex();
843+
//std::cout << "globalIndex:" << mcId << " colID:" << mcCollId << std::endl;
844+
std::vector<float> mcptvec;
845+
float mult = getMCMult(mcParticles, mcCollId, mcptvec);
794846
mcReconstructed[mcCollId] = 1;
795-
mRegistry.fill(HIST("hNTracksMCVsTracksReco"), mult, col.multNTracksGlobal());
847+
for (auto const& trk : tracks) {
848+
int mcPid = trk.mcParticleId(); // depends on your label table
849+
if (mcPid >= 0 && mcPid < (int)mcParticles.size()) {
850+
isReco[mcPid] = 1;
851+
isRecoMult[mcPid] = mult;
852+
} else {
853+
continue;
854+
}
855+
if (trk.collisionId() != collId) {
856+
continue; // different event
857+
}
858+
auto mcPar = mcParticles.rawIteratorAt(mcPid);
859+
// Apply same acceptance as in MC multiplicity
860+
if (!mcPar.isPhysicalPrimary()) {
861+
continue;
862+
}
863+
if (std::abs(mcPar.eta()) > 0.5f) {
864+
continue;
865+
}
866+
int q = 0;
867+
auto pdgEntry = TDatabasePDG::Instance()->GetParticle(mcPar.pdgCode());
868+
if (pdgEntry) {
869+
q = int(std::round(pdgEntry->Charge() / 3.0));
870+
}
871+
if(q == 0) {
872+
continue;
873+
}
874+
//float multReco = recoMult[collId];
875+
float multReco = col.multNTracksGlobal();
876+
float ptReco = trk.pt();
877+
float ptMC = mcPar.pt();
878+
mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRM"), mult, multReco, ptMC, ptReco);
879+
880+
}
881+
882+
//mRegistry.fill(HIST("hNTracksMCVsTracksReco"), mult, col.multNTracksGlobal());
796883
}
884+
// count mc particles with no reco tracks
885+
for (auto const& mcp : mcParticles) {
886+
int mcPidG = mcp.globalIndex();
887+
//std::cout << "mcPidG:" << mcPidG << std::endl;
888+
if (!isReco[mcPidG]) {
889+
mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoTrk"), isRecoMult[mcPidG], mcp.pt());
890+
}
891+
}
892+
// count unreconstructed mc collisions
797893
for (auto const& mc : mcCollisions) {
798894
int gindex = mc.globalIndex();
799-
// std::cout << "mc globalIndex:" << gindex << std::endl;
895+
//std::cout << "mc globalIndex:" << gindex << std::endl;
800896
if (!mcReconstructed[gindex]) {
801-
int mult = getMCMult(mcParticles, gindex);
802-
// std::cout << "===> unreconstructed:" << mult << std::endl;
803-
mRegistry.fill(HIST("hNTracksMCNotInReco"), mult);
897+
std::vector<float> mcptvec;
898+
int mult = getMCMult(mcParticles, gindex, mcptvec);
899+
//std::cout << "===> unreconstructed:" << mult << std::endl;
900+
for(auto const& pt: mcptvec){
901+
mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoCol"), mult,pt);
902+
}
804903
}
805904
}
806905
}
807-
PROCESS_SWITCH(NonPromptCascadeTask, processNegMC, "process mc", false);
906+
PROCESS_SWITCH(NonPromptCascadeTask, processdNdetaMC, "process mc dN/deta", false);
808907
};
809908

810909
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)