diff --git a/DPG/Tasks/AOTEvent/CMakeLists.txt b/DPG/Tasks/AOTEvent/CMakeLists.txt index a7e69449b34..73adb8c9344 100644 --- a/DPG/Tasks/AOTEvent/CMakeLists.txt +++ b/DPG/Tasks/AOTEvent/CMakeLists.txt @@ -58,3 +58,8 @@ o2physics_add_dpl_workflow(bc-occupancy-qa SOURCES bcOccupancyQa.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(mc-event-validation + SOURCES mcEventValidation.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/DPG/Tasks/AOTEvent/lightIonsEvSelQa.cxx b/DPG/Tasks/AOTEvent/lightIonsEvSelQa.cxx index 3a5579ebe51..88e360a88a8 100644 --- a/DPG/Tasks/AOTEvent/lightIonsEvSelQa.cxx +++ b/DPG/Tasks/AOTEvent/lightIonsEvSelQa.cxx @@ -1007,7 +1007,7 @@ struct LightIonsEvSelQa { // Ne-Ne if (lastRunNumber == 564468) meanDiff = -0.51; - if (lastRunNumber == 564468) + if (lastRunNumber == 564472) meanDiff = -0.60; float stdDev = (multT0M > 10) ? 0.144723 + 13.5345 / sqrt(multT0M) : 1.5; // cm diff --git a/DPG/Tasks/AOTEvent/macros/plotMcEventValidation.C b/DPG/Tasks/AOTEvent/macros/plotMcEventValidation.C new file mode 100644 index 00000000000..7b1bb83e92e --- /dev/null +++ b/DPG/Tasks/AOTEvent/macros/plotMcEventValidation.C @@ -0,0 +1,539 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file plotMcEventValidation.C +/// \brief Plots and summary of the MC event validation (output of o2-analysis-mc-event-validation): +/// Nch_gen vs FT0C, BC association, vertex reconstruction efficiency and purity +/// +/// Usage: +/// root -l -b -q 'plotMcEventValidation.C("AnalysisResults.root")' +/// root -l -b -q 'plotMcEventValidation.C("AnalysisResults.root", "mcEventValidation.pdf", "mc-event-validation-task", 4)' +/// (the last argument: rebinning of the distance axis for the efficiency and purity vs distance to the nearest other MC collision) +/// Output: multi-page PDF, and a ROOT file with the canvases and the 2D projections (name of the PDF with .root) +/// +/// \author Igor Altsybeev + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace mcevval +{ +// histogram from the task directory (e.g. "Reco/hNchVsFT0C"), detached from the file +template +T* getHist(TDirectory* dir, const char* name) +{ + auto* h = dynamic_cast(dir->Get(name)); + if (!h) { + std::printf("histogram %s not found\n", name); + return nullptr; + } + h->SetDirectory(nullptr); + h->SetStats(false); + return h; +} + +// 2D projection (x: FT0C, y: multiplicity) of a TH3 for the categories [catFirst, catLast] (bins of the z axis, from 1) +TH2* project(TH3* h3, int catFirst, int catLast, const char* name, const char* title) +{ + h3->GetZaxis()->SetRange(catFirst, catLast); + auto* h2 = dynamic_cast(h3->Project3D("yx")); + h3->GetZaxis()->SetRange(0, 0); + if (!h2) { + return nullptr; + } + h2->SetName(name); + h2->SetTitle(title); + h2->SetStats(false); + return h2; +} + +// draw a 2D histogram with log z in pad i +void draw2D(TCanvas* c, int i, TH2* h) +{ + if (!h) { + return; + } + c->cd(i)->SetLogz(); + c->cd(i)->SetRightMargin(0.13); + h->Draw("colz"); +} + +// efficiency reco / gen (Bayesian 68% intervals; bins without generated collisions are skipped) +TGraphAsymmErrors* efficiency(TH1* reco, TH1* gen, const char* name, int color) +{ + auto* eff = new TGraphAsymmErrors(); + const int errorLevel = gErrorIgnoreLevel; + gErrorIgnoreLevel = kError; // no warning about skipped empty bins + eff->Divide(reco, gen, "cl=0.683 b(1,1) mode"); + gErrorIgnoreLevel = errorLevel; + eff->SetName(name); + eff->SetLineColor(color); + eff->SetMarkerColor(color); + eff->SetMarkerStyle(kFullCircle); + eff->SetMarkerSize(0.6); + return eff; +} + +// mean of a profile over all bins (weighted with the bin entries) +double profileMean(const TProfile* p) +{ + double sum = 0, n = 0; + for (int i = 1; i <= p->GetNbinsX(); i++) { + sum += p->GetBinContent(i) * p->GetBinEntries(i); + n += p->GetBinEntries(i); + } + return n > 0 ? sum / n : -1.; +} +// profile as a graph: bins without entries are skipped +TGraphErrors* profileGraph(const TProfile* p, const char* name, int color, int marker) +{ + auto* g = new TGraphErrors(); + g->SetName(name); + for (int i = 1; i <= p->GetNbinsX(); i++) { + if (p->GetBinEntries(i) > 0) { + const int n = g->GetN(); + g->SetPoint(n, p->GetBinCenter(i), p->GetBinContent(i)); + g->SetPointError(n, p->GetBinWidth(i) / 2., p->GetBinError(i)); + } + } + g->SetLineColor(color); + g->SetMarkerColor(color); + g->SetMarkerStyle(marker); + g->SetMarkerSize(0.6); + return g; +} +} // namespace mcevval + +void plotMcEventValidation(const char* fileName = "AnalysisResults.root", const char* outPdf = "mcEventValidation.pdf", const char* dirName = "mc-event-validation-task", int rebinDist = 4) +{ + using namespace mcevval; + TH1::AddDirectory(false); + gStyle->SetOptStat(0); + gStyle->SetPalette(kViridis); + gStyle->SetPadLeftMargin(0.14); + + std::unique_ptr f(TFile::Open(fileName, "READ")); + if (!f || f->IsZombie()) { + std::printf("cannot open %s\n", fileName); + return; + } + TDirectory* dir = f->GetDirectory(dirName); + if (!dir) { + std::printf("directory %s not found in %s\n", dirName, fileName); + return; + } + + // MC truth + auto* hNchVsB = getHist(dir, "MC/hNch"); + auto* hNchFromVertexVsB = getHist(dir, "MC/hNchFromVertex"); + auto* pFracDisplaced = getHist(dir, "MC/pFractionDisplacedVsB"); + auto* hNchVsNchFT0C = getHist(dir, "MC/hNchVsNchFT0C"); + auto* hNchVsNchFT0CFromVertex = getHist(dir, "MC/hNchVsNchFT0CFromVertex"); + auto* hDispPt = getHist(dir, "MC/Displaced/hPt"); + auto* hDispDz = getHist(dir, "MC/Displaced/hDeltaZ"); + auto* hDispR = getHist(dir, "MC/Displaced/hR"); + auto* hDispSpecies = getHist(dir, "MC/Displaced/hSpecies"); + auto* hDispGenStatus = getHist(dir, "MC/Displaced/hGenStatus"); + auto* hVtxGenStatus = getHist(dir, "MC/FromVertex/hGenStatus"); + auto* hNumMcCollsInBC = getHist(dir, "MC/hNumMcCollsInBC"); + auto* hNumRecoPerMcColl = getHist(dir, "MC/hNumRecoPerMcColl"); + // reconstructed collisions + auto* hCounter = getHist(dir, "Reco/hCounter"); + auto* hCategory = getHist(dir, "Reco/hCategory"); + auto* hDeltaBC = getHist(dir, "Reco/hDeltaBC"); + auto* hPurity = getHist(dir, "Reco/hPurity"); + auto* h3NchVsFT0C = getHist(dir, "Reco/hNchVsFT0C"); + auto* h3NchFromVertexVsFT0C = getHist(dir, "Reco/hNchFromVertexVsFT0C"); + auto* h3NGlobalVsFT0C = getHist(dir, "Reco/hNGlobalVsFT0C"); + auto* h3NGlobalVsNch = getHist(dir, "Reco/hNGlobalVsNchFromVertex"); + // vertex reconstruction efficiency and purity + auto* hNchGen = getHist(dir, "Efficiency/hNchGen"); + auto* hNchGenReco = getHist(dir, "Efficiency/hNchGenReco"); + auto* hNchGenWide = getHist(dir, "Efficiency/hNchGenWide"); + auto* hNchGenRecoWide = getHist(dir, "Efficiency/hNchGenRecoWide"); + auto* hDistGen = getHist(dir, "Efficiency/hDistGen"); + auto* hDistGenReco = getHist(dir, "Efficiency/hDistGenReco"); + auto* pPurityPVVsNch = getHist(dir, "Purity/pPurityPVVsNch"); + auto* pPurityAssocVsNch = getHist(dir, "Purity/pPurityAssocVsNch"); + auto* hPurityPVVsNch = getHist(dir, "Purity/hPurityPVVsNch"); + auto* hPurityAssocVsNch = getHist(dir, "Purity/hPurityAssocVsNch"); + auto* pPurityPVVsDist = getHist(dir, "Purity/pPurityPVVsDist"); + auto* pPurityAssocVsDist = getHist(dir, "Purity/pPurityAssocVsDist"); + + // ### summary + std::printf("\n=== MC event validation: %s ===\n", fileName); + if (hNchVsB && hNchFromVertexVsB && hNchVsB->GetEntries() > 0) { + const double meanAll = hNchVsB->GetMean(2), meanVtx = hNchFromVertexVsB->GetMean(2); + std::printf("MC collisions: %.0f; mean Nch(|eta|<0.5): all charged physical primaries %.1f, produced at the vertex %.1f -> %.1f%% displaced\n", + hNchVsB->GetEntries(), meanAll, meanVtx, meanAll > 0 ? 100. * (1. - meanVtx / meanAll) : 0.); + } + if (pFracDisplaced) { + for (const auto& [bMin, bMax] : {std::array{0., 5.}, std::array{5., 10.}, std::array{10., 20.}}) { + double sum = 0, n = 0; + for (int i = pFracDisplaced->FindBin(bMin + 1e-6); i <= pFracDisplaced->FindBin(bMax - 1e-6); i++) { + sum += pFracDisplaced->GetBinContent(i) * pFracDisplaced->GetBinEntries(i); + n += pFracDisplaced->GetBinEntries(i); + } + std::printf(" fraction of Nch(|eta|<0.5) not produced at the vertex, %4.1f < b < %4.1f fm: %.3f\n", bMin, bMax, n > 0 ? sum / n : 0.); + } + } + if (hNchVsNchFT0C && hNchVsNchFT0CFromVertex) { + std::printf("correlation factor Nch(|eta|<0.5) vs Nch(FT0C acceptance): all primaries %.4f, produced at the vertex %.4f\n", + hNchVsNchFT0C->GetCorrelationFactor(), hNchVsNchFT0CFromVertex->GetCorrelationFactor()); + } + if (hDispSpecies && hDispSpecies->GetEntries() > 0) { + std::printf("displaced charged physical primaries: e- %.1f%%, e+ %.1f%%, other %.1f%%; mean pT %.4f GeV/c\n", 100. * hDispSpecies->GetBinContent(1) / hDispSpecies->GetEntries(), + 100. * hDispSpecies->GetBinContent(2) / hDispSpecies->GetEntries(), 100. * hDispSpecies->GetBinContent(3) / hDispSpecies->GetEntries(), hDispPt ? hDispPt->GetMean() : 0.); + } + if (hNumMcCollsInBC && hNumMcCollsInBC->GetEntries() > 0) { + std::printf("MC collisions sharing their BC with other MC collisions (same-bunch pile-up in MC): %.2f%%\n", 100. * (1. - hNumMcCollsInBC->GetBinContent(1) / hNumMcCollsInBC->GetEntries())); + } + if (hNumRecoPerMcColl && hNumRecoPerMcColl->GetEntries() > 0) { + const double n = hNumRecoPerMcColl->GetEntries(); + std::printf("MC collisions with 0 / 1 / >=2 reconstructed collisions: %.2f%% / %.2f%% / %.2f%%\n", 100. * hNumRecoPerMcColl->GetBinContent(1) / n, 100. * hNumRecoPerMcColl->GetBinContent(2) / n, + 100. * hNumRecoPerMcColl->Integral(3, hNumRecoPerMcColl->GetNbinsX() + 1) / n); + } + if (hCounter) { + std::printf("reconstructed collisions:"); + for (int i = 1; i <= hCounter->GetNbinsX(); i++) { + std::printf(" [%s] %.0f", hCounter->GetXaxis()->GetBinLabel(i), hCounter->GetBinContent(i)); + } + std::printf("\n"); + } + if (hCategory && hCategory->GetEntries() > 0) { + std::printf("selected reconstructed collisions by found BC:"); + for (int i = 1; i <= hCategory->GetNbinsX(); i++) { + std::printf(" [%s] %.2f%%", hCategory->GetXaxis()->GetBinLabel(i), 100. * hCategory->GetBinContent(i) / hCategory->GetEntries()); + } + std::printf("\n"); + } + if (hPurity) { + for (int i = 1; i <= hPurity->GetNbinsX(); i++) { + std::unique_ptr p(hPurity->ProjectionY(Form("purity_%d", i), i, i)); + if (p->GetEntries() > 0) { + std::printf(" [%s] mean fraction of PV contributors from the labelled MC collision: %.3f\n", hPurity->GetXaxis()->GetBinLabel(i), p->GetMean()); + } + } + } + if (hNchGenWide && hNchGenRecoWide && hNchGenWide->GetEntries() > 0) { + std::printf("vertex reconstruction efficiency (matched reconstructed collision): all MC collisions %.4f (%.0f / %.0f)\n", hNchGenRecoWide->GetEntries() / hNchGenWide->GetEntries(), + hNchGenRecoWide->GetEntries(), hNchGenWide->GetEntries()); + } + if (hNchGen && hNchGenReco && hNchGenWide && hNchGenRecoWide) { + // classes by Nch (|eta|<0.8), as in the task: [0, 10), [10, 20), [20, 30), >= 30 + const std::array classEdges = {0, 10, 20, 30, -1}; // o2-linter: disable=magic-number (Nch classes, as in the task) + for (size_t i = 0; i + 1 < classEdges.size(); i++) { + double nGen = 0, nReco = 0; + if (classEdges[i + 1] > 0) { + nGen = hNchGen->Integral(hNchGen->FindBin(classEdges[i]), hNchGen->FindBin(classEdges[i + 1] - 1)); + nReco = hNchGenReco->Integral(hNchGenReco->FindBin(classEdges[i]), hNchGenReco->FindBin(classEdges[i + 1] - 1)); + std::printf(" Nch(|eta|<0.8) in [%d, %d): ", classEdges[i], classEdges[i + 1]); + } else { // last class: all MC collisions minus those below the edge + nGen = hNchGenWide->GetEntries() - hNchGen->Integral(1, hNchGen->FindBin(classEdges[i] - 1)); + nReco = hNchGenRecoWide->GetEntries() - hNchGenReco->Integral(1, hNchGenReco->FindBin(classEdges[i] - 1)); + std::printf(" Nch(|eta|<0.8) >= %d: ", classEdges[i]); + } + std::printf("efficiency %.4f (%.0f / %.0f)\n", nGen > 0 ? nReco / nGen : 0., nReco, nGen); + } + } + if (hPurityPVVsNch && hPurityAssocVsNch) { + std::printf("mean purity of reconstructed collisions wrt the labelled MC collision: PV contributors %.4f, associated tracks %.4f\n", hPurityPVVsNch->GetMean(2), hPurityAssocVsNch->GetMean(2)); + } + if (pPurityPVVsDist && pPurityAssocVsDist) { + for (int i = 1; i <= pPurityPVVsDist->GetNbinsY(); i++) { + std::unique_ptr pPV(pPurityPVVsDist->ProfileX(Form("pPV_%d", i), i, i)); + std::unique_ptr pAssoc(pPurityAssocVsDist->ProfileX(Form("pAssoc_%d", i), i, i)); + std::printf(" [%s] mean purity: PV contributors %.4f, associated tracks %.4f (%.0f collisions)\n", pPurityPVVsDist->GetYaxis()->GetBinLabel(i), profileMean(pPV.get()), + profileMean(pAssoc.get()), pPV->GetEntries()); + } + } + + // ### plots + std::vector canvases; + auto newPage = [&canvases](const char* name, const char* title, int width, int height, int nx, int ny) { + auto* canvas = new TCanvas(name, title, width, height); + canvas->Divide(nx, ny); + canvases.push_back(canvas); + return canvas; + }; + const int nCat = h3NchVsFT0C ? h3NchVsFT0C->GetNbinsZ() : 0; + constexpr int NCategories = 3; // categories of the found BC: correct, shared with other MC collisions, wrong + + // page 1: MC truth, Nch(|eta|<0.5) vs Nch in the FT0C acceptance and vs impact parameter + TCanvas* c = newPage("cMcTruth", "MC truth", 1500, 1000, 2, 2); + if (hNchVsNchFT0C) { + hNchVsNchFT0C->SetTitle("MC: all charged physical primaries (as in MultModule)"); + } + if (hNchVsNchFT0CFromVertex) { + hNchVsNchFT0CFromVertex->SetTitle("MC: charged physical primaries produced at the vertex"); + } + draw2D(c, 1, hNchVsNchFT0C); + draw2D(c, 2, hNchVsNchFT0CFromVertex); + draw2D(c, 3, hNchVsB); + if (pFracDisplaced) { + c->cd(4); + pFracDisplaced->SetMinimum(0.); + pFracDisplaced->SetMaximum(1.); + pFracDisplaced->Draw(); + } + + // page 2: displaced charged physical primaries + c = newPage("cDisplaced", "displaced charged physical primaries", 1500, 1000, 3, 2); + int pad = 1; + for (const auto& h : {hDispPt, hDispDz, hDispR}) { + if (h) { + c->cd(pad)->SetLogy(); + h->Draw("hist"); + } + pad++; + } + if (hDispSpecies) { + c->cd(4); + hDispSpecies->SetMinimum(0.); + hDispSpecies->Draw("hist text"); + } + if (hDispGenStatus && hVtxGenStatus) { + c->cd(5)->SetLogy(); + hVtxGenStatus->SetTitle("generator status of charged physical primaries;generator status code;particles"); + hVtxGenStatus->SetLineColor(kBlue + 1); + hDispGenStatus->SetLineColor(kRed + 1); + hVtxGenStatus->SetMinimum(0.5); + hVtxGenStatus->SetMaximum(10. * std::max(hVtxGenStatus->GetMaximum(), hDispGenStatus->GetMaximum())); + hVtxGenStatus->Draw("hist"); + hDispGenStatus->Draw("hist same"); + auto* leg = new TLegend(0.45, 0.75, 0.89, 0.89); + leg->AddEntry(hVtxGenStatus, "produced at the vertex", "l"); + leg->AddEntry(hDispGenStatus, "displaced", "l"); + leg->Draw(); + } + + // page 3: reconstructed collisions, Nch_gen vs FT0C + c = newPage("cNchVsFT0C", "Nch_gen vs FT0C", 1500, 1000, 2, 2); + std::vector derived; // projections, profiles and efficiencies saved in the output file + if (h3NchVsFT0C && h3NchFromVertexVsFT0C && nCat >= NCategories) { + std::array hs = {project(h3NchVsFT0C, 1, nCat, "hNchVsFT0C_all", "all primaries (as in MultModule), all collisions;FT0C amplitude (found BC);N_{ch} (|#eta|<0.5)"), + project(h3NchFromVertexVsFT0C, 1, nCat, "hNchFromVertexVsFT0C_all", "primaries from the vertex, all collisions;FT0C amplitude (found BC);N_{ch} (|#eta|<0.5)"), + project(h3NchFromVertexVsFT0C, 1, 1, "hNchFromVertexVsFT0C_correctBC", "primaries from the vertex, correct BC;FT0C amplitude (found BC);N_{ch} (|#eta|<0.5)"), + project(h3NchFromVertexVsFT0C, 2, 3, "hNchFromVertexVsFT0C_sharedOrWrongBC", "primaries from the vertex, shared or wrong BC;FT0C amplitude (found BC);N_{ch} (|#eta|<0.5)")}; + for (int i = 0; i < 4; i++) { // o2-linter: disable=magic-number (number of pads) + draw2D(c, i + 1, hs[i]); + derived.push_back(hs[i]); + } + } + + // page 4: reconstructed collisions, global-like tracks vs FT0C and vs Nch_gen + c = newPage("cTracksVsFT0C", "tracks vs FT0C", 1500, 1000, 3, 2); + if (h3NGlobalVsFT0C && nCat >= NCategories) { + std::array hs = {project(h3NGlobalVsFT0C, 1, nCat, "hNGlobalVsFT0C_all", "all collisions;FT0C amplitude (found BC);global-like tracks"), + project(h3NGlobalVsFT0C, 1, 1, "hNGlobalVsFT0C_correctBC", "correct BC;FT0C amplitude (found BC);global-like tracks"), + project(h3NGlobalVsFT0C, 2, 2, "hNGlobalVsFT0C_sharedBC", "BC shared with other MC collisions;FT0C amplitude (found BC);global-like tracks"), + project(h3NGlobalVsFT0C, 3, 3, "hNGlobalVsFT0C_wrongBC", "wrong BC;FT0C amplitude (found BC);global-like tracks")}; + for (int i = 0; i < 4; i++) { // o2-linter: disable=magic-number (number of pads) + draw2D(c, i + 1, hs[i]); + derived.push_back(hs[i]); + } + } + if (h3NGlobalVsNch && nCat >= NCategories) { + auto* h = project(h3NGlobalVsNch, 1, nCat, "hNGlobalVsNchFromVertex_all", "all collisions;N_{ch} (|#eta|<0.5) from the vertex;global-like tracks"); + draw2D(c, 5, h); + derived.push_back(h); + } + + // page 5: BC association, same-bunch pile-up in MC, vertex purity and splitting + c = newPage("cAssociation", "BC association and pile-up", 1500, 1000, 3, 2); + if (hCategory) { + c->cd(1)->SetBottomMargin(0.25); + hCategory->SetMinimum(0.); + hCategory->Draw("hist text"); + } + if (hDeltaBC) { + c->cd(2)->SetLogy(); + hDeltaBC->Draw("hist"); + } + if (hNumMcCollsInBC) { + c->cd(3)->SetLogy(); + hNumMcCollsInBC->Draw("hist text"); + } + if (hPurity) { + c->cd(4)->SetLogy(); + auto* leg = new TLegend(0.15, 0.65, 0.6, 0.89); + const std::array colors = {kBlue + 1, kOrange + 7, kRed + 1}; + for (int i = 1; i <= hPurity->GetNbinsX(); i++) { + auto* p = hPurity->ProjectionY(Form("hPurity_%d", i), i, i); + p->SetTitle("fraction of PV contributors from the labelled MC collision;purity;collisions"); + p->SetLineColor(colors[(i - 1) % 3]); + p->SetStats(false); + p->Draw(i == 1 ? "hist" : "hist same"); + leg->AddEntry(p, hPurity->GetXaxis()->GetBinLabel(i), "l"); + derived.push_back(p); + } + leg->Draw(); + } + if (hNumRecoPerMcColl) { + c->cd(5)->SetLogy(); + hNumRecoPerMcColl->Draw("hist text"); + } + + // page 6: vertex reconstruction efficiency + c = newPage("cEfficiency", "vertex reconstruction efficiency", 1500, 1000, 2, 2); + const std::array classColors = {kRed + 1, kOrange + 7, kGreen + 2, kBlue + 1}; + pad = 1; + for (const auto& [reco, gen] : {std::pair{hNchGenReco, hNchGen}, std::pair{hNchGenRecoWide, hNchGenWide}}) { + if (reco && gen) { + c->cd(pad)->SetGridy(); + gPad->DrawFrame(gen->GetXaxis()->GetXmin(), 0., gen->GetXaxis()->GetXmax(), 1.1, + pad == 1 ? "vertex reconstruction efficiency;N_{ch} (|#eta|<0.8) from the vertex;efficiency" : "vertex reconstruction efficiency (full range);N_{ch} (|#eta|<0.8) from the vertex;efficiency"); + auto* eff = efficiency(reco, gen, pad == 1 ? "gEffVsNch" : "gEffVsNchWide", kBlack); + eff->Draw("p"); + derived.push_back(eff); + } + pad++; + } + if (hDistGen && hDistGenReco) { + auto* leg = new TLegend(0.55, 0.15, 0.89, 0.4); + auto* legGen = new TLegend(0.55, 0.65, 0.89, 0.89); + TH1* genFirst = nullptr; // first drawn class: its maximum sets the y range for all classes + double genMax = 0; + for (int i = 1; i <= hDistGen->GetNbinsY(); i++) { + auto* gen = hDistGen->ProjectionX(Form("hDistGen_%d", i), i, i); + auto* reco = hDistGenReco->ProjectionX(Form("hDistGenReco_%d", i), i, i); + gen->Rebin(rebinDist); + reco->Rebin(rebinDist); + const int color = classColors[(i - 1) % classColors.size()]; + c->cd(3)->SetGridy(); + if (i == 1) { + gPad->DrawFrame(gen->GetXaxis()->GetXmin(), 0., gen->GetXaxis()->GetXmax(), 1.1, "vertex reconstruction efficiency;distance to the nearest other MC collision (BCs);efficiency"); + } + auto* eff = efficiency(reco, gen, Form("gEffVsDist_%d", i), color); + eff->Draw("p"); + leg->AddEntry(eff, hDistGen->GetYaxis()->GetBinLabel(i), "lp"); + c->cd(4)->SetLogy(); + gen->SetTitle("MC collisions;distance to the nearest other MC collision (BCs);MC collisions"); + gen->SetStats(false); + gen->SetLineColor(color); + gen->SetMinimum(0.5); + gen->Draw(i == 1 ? "hist" : "hist same"); + genFirst = genFirst ? genFirst : gen; + genMax = std::max(genMax, gen->GetMaximum()); + legGen->AddEntry(gen, hDistGen->GetYaxis()->GetBinLabel(i), "l"); + derived.push_back(eff); + derived.push_back(gen); + derived.push_back(reco); + } + c->cd(3); + leg->Draw(); + c->cd(4); + if (genFirst) { + genFirst->SetMaximum(5. * genMax); + } + legGen->Draw(); + } + + // page 7: purity of the reconstructed collisions wrt the labelled MC collision + c = newPage("cPurity", "purity of reconstructed collisions", 1500, 1000, 3, 2); + if (pPurityPVVsNch && pPurityAssocVsNch) { + c->cd(1)->SetGridy(); + gPad->DrawFrame(pPurityPVVsNch->GetXaxis()->GetXmin(), 0., pPurityPVVsNch->GetXaxis()->GetXmax(), 1.1, "purity of reconstructed collisions;N_{ch} (|#eta|<0.8) of the labelled MC collision;purity"); + auto* gPV = profileGraph(pPurityPVVsNch, "gPurityPVVsNch", kBlue + 1, kFullCircle); + auto* gAssoc = profileGraph(pPurityAssocVsNch, "gPurityAssocVsNch", kRed + 1, kOpenCircle); + gPV->Draw("p"); + gAssoc->Draw("p"); + auto* leg = new TLegend(0.45, 0.15, 0.89, 0.3); + leg->AddEntry(gPV, "PV contributors", "lp"); + leg->AddEntry(gAssoc, "all associated tracks", "lp"); + leg->Draw(); + derived.push_back(gPV); + derived.push_back(gAssoc); + } + if (hPurityPVVsNch) { + hPurityPVVsNch->SetTitle("purity: PV contributors;N_{ch} (|#eta|<0.8) of the labelled MC collision;purity"); + } + if (hPurityAssocVsNch) { + hPurityAssocVsNch->SetTitle("purity: all associated tracks;N_{ch} (|#eta|<0.8) of the labelled MC collision;purity"); + } + draw2D(c, 2, hPurityPVVsNch); + draw2D(c, 3, hPurityAssocVsNch); + pad = 4; + for (const auto& [p2, what] : {std::pair{pPurityPVVsDist, "PV contributors"}, std::pair{pPurityAssocVsDist, "all associated tracks"}}) { + if (p2) { + c->cd(pad)->SetGridy(); + gPad->DrawFrame(p2->GetXaxis()->GetXmin(), 0., p2->GetXaxis()->GetXmax(), 1.1, Form("purity: %s;distance to the nearest other MC collision (BCs);purity", what)); + auto* leg = new TLegend(0.42, 0.15, 0.89, 0.4); + for (int i = 1; i <= p2->GetNbinsY(); i++) { + std::unique_ptr p(p2->ProfileX(Form("%s_%d", p2->GetName(), i), i, i)); + p->Rebin(rebinDist); + auto* g = profileGraph(p.get(), Form("g%s_%d", p2->GetName(), i), classColors[(i - 1) % classColors.size()], kFullCircle); + g->Draw("p"); + leg->AddEntry(g, p2->GetYaxis()->GetBinLabel(i), "lp"); + derived.push_back(g); + } + leg->Draw(); + } + pad++; + } + + // PDF: each page with the aspect ratio of its canvas + for (size_t i = 0; i < canvases.size(); i++) { + gStyle->SetPaperSize(30., 30. * canvases[i]->GetWh() / canvases[i]->GetWw()); + if (i == 0) { + canvases[i]->Print(Form("%s[", outPdf)); + } + canvases[i]->Print(outPdf); + if (i == canvases.size() - 1) { + canvases[i]->Print(Form("%s]", outPdf)); + } + } + std::printf("plots saved in %s\n", outPdf); + + // ROOT file with the canvases + TString rootName = outPdf; + if (rootName.EndsWith(".pdf")) { + rootName.Resize(rootName.Length() - 4); + } + rootName += ".root"; + std::unique_ptr fOut(TFile::Open(rootName, "RECREATE")); + if (fOut && !fOut->IsZombie()) { + fOut->mkdir("canvases")->cd(); + for (const auto& canvas : canvases) { + canvas->Write(); + } + fOut->mkdir("projections")->cd(); + for (const auto& h : derived) { + h->Write(); + } + std::printf("canvases and projections saved in %s\n", rootName.Data()); + } +} diff --git a/DPG/Tasks/AOTEvent/mcEventValidation.cxx b/DPG/Tasks/AOTEvent/mcEventValidation.cxx new file mode 100644 index 00000000000..d2e56d02b26 --- /dev/null +++ b/DPG/Tasks/AOTEvent/mcEventValidation.cxx @@ -0,0 +1,448 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file mcEventValidation.cxx +/// \brief MC validation of the event reconstruction: generated multiplicity vs FT0C, vertex reconstruction efficiency and purity +/// +/// Investigates the off-diagonal populations in Nch_gen vs FT0C and in Ntracks vs FT0C (O2-6499): +/// - MC truth: charged physical primaries (as in MultModule) split into those produced at the MC collision vertex and the displaced ones; +/// - reconstructed collisions: FT0C from the found BC, categorised by the BC association (correct BC, BC shared with other MC collisions, wrong BC); +/// - number of MC collisions per BC (same-bunch pile-up in MC), PV contributors from other MC collisions, number of reconstructed collisions per MC collision; +/// - vertex reconstruction efficiency and purity vs the generated multiplicity (|eta|<0.8) and vs the distance (in BCs) to the nearest other MC collision. +/// +/// \author Igor Altsybeev + +#include "Common/CCDB/EventSelectionParams.h" +#include "Common/DataModel/EventSelection.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::aod::evsel; + +using ColsMC = soa::Join; +using TracksMC = soa::Join; + +struct McEventValidationTask { + Service pdg; + Preslice perCollision = aod::track::collisionId; + + // event selection (as in the Nch_gen vs FT0C studies: sel8, |vZ| < 10 cm, kNoSameBunchPileup) + Configurable requireSel8{"requireSel8", true, "require sel8 for reconstructed collisions"}; + Configurable requireNoSameBunchPileup{"requireNoSameBunchPileup", true, "require kNoSameBunchPileup for reconstructed collisions"}; + Configurable maxAbsVtxZ{"maxAbsVtxZ", 10.f, "max |vZ| of reconstructed collisions (cm)"}; + // MC particles + Configurable maxDistFromVertex{"maxDistFromVertex", 1.f, "particles produced within this distance (cm, in r and in z) from the MC collision vertex are counted as 'from the vertex'"}; + // tracks: PV contributors with ITS, and "global-like" tracks with TPC + Configurable minPt{"minPt", 0.2f, "min pT of tracks (GeV/c)"}; + Configurable maxAbsEta{"maxAbsEta", 0.8f, "max |eta| of tracks"}; + Configurable minITSNCls{"minITSNCls", 5, "min number of ITS clusters of tracks"}; + Configurable minTPCNCls{"minTPCNCls", 50, "min number of TPC clusters of global-like tracks"}; + // vertex reconstruction efficiency: an MC collision is reconstructed if a reconstructed collision has more than minFractionMatch of its PV contributors from it + // and the reconstructed vZ agrees with the generated one within maxDeltaVzMatch + Configurable minFractionMatch{"minFractionMatch", 0.5f, "min fraction of PV contributors from the MC collision for the matching"}; + Configurable maxDeltaVzMatch{"maxDeltaVzMatch", 0.3f, "max |vZ(reco) - vZ(gen)| for the matching (cm)"}; + Configurable maxDistBC{"maxDistBC", 800, "upper edge of the axis for the distance (in BCs) to the nearest other MC collision"}; + Configurable nBinsDistBC{"nBinsDistBC", 160, "number of bins for the distance (in BCs) to the nearest other MC collision"}; + // axes + Configurable nBinsNch{"nBinsNch", 250, "number of bins for Nch and number of tracks"}; + Configurable maxNch{"maxNch", 5000.f, "upper edge of the axes for Nch and number of tracks"}; + Configurable nBinsFT0C{"nBinsFT0C", 300, "number of bins for the FT0C amplitude"}; + Configurable maxFT0C{"maxFT0C", 60000.f, "upper edge of the axis for the FT0C amplitude"}; + + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + // categories of reconstructed collisions by the association of the found BC to the true BC of the MC collision + enum BcCategory { kCorrectBC = 0, // found BC = BC of the MC collision, no other MC collisions in this BC + kSharedBC, // found BC = BC of the MC collision, other MC collisions in this BC (same-bunch pile-up in MC) + kWrongBC, // found BC != BC of the MC collision + kNBcCategories }; + + static constexpr float MaxEtaNch = 0.5f; // |eta| range for Nch, as in MultModule + static constexpr float MinEtaFT0C = -3.3f; // FT0C acceptance, as in MultModule + static constexpr float MaxEtaFT0C = -2.1f; // FT0C acceptance, as in MultModule + static constexpr float MinCharge = 1e-3f; // charged particles + static constexpr int MaxDeltaBC = 50; // axis range for found BC - true BC + static constexpr int MaxNumMcCollsInBC = 10; // axis range for the number of MC collisions in a BC + static constexpr int MaxNumRecoPerMcColl = 10; // axis range for the number of reconstructed collisions per MC collision + static constexpr float MaxEtaNch08 = 0.8f; // |eta| range for Nch in the efficiency and purity studies + static constexpr int MaxNchLow = 200; // axis range for Nch (|eta|<0.8) at low multiplicity, one bin per particle + static constexpr int NNchClasses = 4; // classes of MC collisions by Nch (|eta|<0.8): [0, 10), [10, 20), [20, 30), >= 30 + static constexpr std::array NchClassLimits = {10, 20, 30}; + + // class of an MC collision by Nch (|eta|<0.8): index of the first class limit above nch + static int nchClass(int nch) + { + return static_cast(std::upper_bound(NchClassLimits.begin(), NchClassLimits.end(), nch) - NchClassLimits.begin()); + } + + void init(InitContext&) + { + const AxisSpec axisNch{nBinsNch.value, 0., maxNch.value, "N_{ch} (|#eta|<0.5)"}; + const AxisSpec axisNchFT0C{nBinsNch.value, 0., maxNch.value, "N_{ch} in FT0C acceptance (-3.3<#eta<-2.1)"}; + const AxisSpec axisNtracks{nBinsNch.value, 0., maxNch.value, "number of tracks"}; + const AxisSpec axisFT0C{nBinsFT0C.value, 0., maxFT0C.value, "FT0C amplitude (found BC)"}; + const AxisSpec axisB{200, 0., 20., "impact parameter b (fm)"}; + const AxisSpec axisCategory{kNBcCategories, -0.5, static_cast(kNBcCategories) - 0.5, ""}; + + // MC truth + histos.add("MC/hNch", "MC collisions: charged physical primaries in |#eta|<0.5", kTH2F, {axisB, axisNch}); + histos.add("MC/hNchFromVertex", "MC collisions: charged physical primaries in |#eta|<0.5 produced at the vertex", kTH2F, {axisB, axisNch}); + histos.add("MC/pFractionDisplacedVsB", "MC collisions: fraction of N_{ch} (|#eta|<0.5) not produced at the vertex;impact parameter b (fm);fraction", kTProfile, {axisB}); + histos.add("MC/hNchVsNchFT0C", "MC collisions: all charged physical primaries", kTH2F, {axisNchFT0C, axisNch}); + histos.add("MC/hNchVsNchFT0CFromVertex", "MC collisions: charged physical primaries produced at the vertex", kTH2F, {axisNchFT0C, axisNch}); + // displaced charged physical primaries (all eta) + histos.add("MC/Displaced/hPt", "displaced charged physical primaries;p_{T} (GeV/c)", kTH1F, {{500, 0., 0.05}}); + histos.add("MC/Displaced/hDeltaZ", "displaced charged physical primaries;z_{prod} - z_{vertex} (cm)", kTH1F, {{500, -500., 500.}}); + histos.add("MC/Displaced/hR", "displaced charged physical primaries;r_{prod} (cm)", kTH1F, {{500, 0., 500.}}); + histos.add("MC/Displaced/hSpecies", "displaced charged physical primaries", kTH1F, {{3, -0.5, 2.5}}); + histos.add("MC/Displaced/hGenStatus", "displaced charged physical primaries;generator status code", kTH1F, {{201, -1.5, 199.5}}); + histos.add("MC/FromVertex/hGenStatus", "charged physical primaries produced at the vertex;generator status code", kTH1F, {{201, -1.5, 199.5}}); + // same-bunch pile-up in MC and splitting of MC collisions into several reconstructed ones + histos.add("MC/hNumMcCollsInBC", "number of MC collisions in the BC of a given MC collision;number of MC collisions in the BC", kTH1F, {{MaxNumMcCollsInBC, 0.5, MaxNumMcCollsInBC + 0.5}}); + histos.add("MC/hNumRecoPerMcColl", "MC collisions;number of reconstructed collisions", kTH1F, {{MaxNumRecoPerMcColl + 1, -0.5, MaxNumRecoPerMcColl + 0.5}}); + histos.add("MC/hNumRecoPerMcCollVsB", "MC collisions", kTH2F, {axisB, {MaxNumRecoPerMcColl + 1, -0.5, MaxNumRecoPerMcColl + 0.5, "number of reconstructed collisions"}}); + + // reconstructed collisions + histos.add("Reco/hCounter", "reconstructed collisions", kTH1F, {{6, -0.5, 5.5}}); + histos.add("Reco/hCategory", "selected reconstructed collisions", kTH1F, {axisCategory}); + histos.add("Reco/hDeltaBC", "selected reconstructed collisions;found BC - BC of the MC collision", kTH1F, {{2 * MaxDeltaBC + 1, -MaxDeltaBC - 0.5, MaxDeltaBC + 0.5}}); + histos.add("Reco/hNumMcCollsInFoundBC", "selected reconstructed collisions;number of MC collisions in the found BC", kTH1F, {{MaxNumMcCollsInBC + 1, -0.5, MaxNumMcCollsInBC + 0.5}}); + histos.add("Reco/hPurity", "selected reconstructed collisions;fraction of PV contributors from the labelled MC collision", kTH2F, {axisCategory, {110, 0., 1.1, "purity"}}); + histos.add("Reco/hNchVsFT0C", "selected reconstructed collisions: all charged physical primaries", kTH3F, {axisFT0C, axisNch, axisCategory}); + histos.add("Reco/hNchFromVertexVsFT0C", "selected reconstructed collisions: charged physical primaries produced at the vertex", kTH3F, {axisFT0C, axisNch, axisCategory}); + histos.add("Reco/hNchFT0CFromVertexVsFT0C", "selected reconstructed collisions: N_{ch} in FT0C acceptance produced at the vertex", kTH3F, {axisFT0C, axisNchFT0C, axisCategory}); + histos.add("Reco/hNPVVsFT0C", "selected reconstructed collisions: PV contributors with ITS", kTH3F, {axisFT0C, axisNtracks, axisCategory}); + histos.add("Reco/hNGlobalVsFT0C", "selected reconstructed collisions: global-like tracks", kTH3F, {axisFT0C, axisNtracks, axisCategory}); + histos.add("Reco/hNGlobalVsNchFromVertex", "selected reconstructed collisions: global-like tracks", kTH3F, {axisNch, axisNtracks, axisCategory}); + + // vertex reconstruction efficiency and purity vs Nch (|eta|<0.8, produced at the vertex) and vs the distance to the nearest other MC collision + const AxisSpec axisNch08Low{MaxNchLow, -0.5, MaxNchLow - 0.5, "N_{ch} (|#eta|<0.8)"}; + const AxisSpec axisNch08{nBinsNch.value, 0., maxNch.value, "N_{ch} (|#eta|<0.8)"}; + const AxisSpec axisDist{nBinsDistBC.value, 0., static_cast(maxDistBC.value), "distance to the nearest other MC collision (BCs)"}; + const AxisSpec axisNchClass{NNchClasses, -0.5, NNchClasses - 0.5, ""}; + const AxisSpec axisPurity{110, 0., 1.1, "purity"}; + histos.add("Efficiency/hNchGen", "MC collisions", kTH1F, {axisNch08Low}); + histos.add("Efficiency/hNchGenReco", "MC collisions with a matched reconstructed collision", kTH1F, {axisNch08Low}); + histos.add("Efficiency/hNchGenWide", "MC collisions", kTH1F, {axisNch08}); + histos.add("Efficiency/hNchGenRecoWide", "MC collisions with a matched reconstructed collision", kTH1F, {axisNch08}); + histos.add("Efficiency/hDistGen", "MC collisions", kTH2F, {axisDist, axisNchClass}); + histos.add("Efficiency/hDistGenReco", "MC collisions with a matched reconstructed collision", kTH2F, {axisDist, axisNchClass}); + histos.add("Purity/pPurityPVVsNch", "reconstructed collisions;N_{ch} (|#eta|<0.8) of the labelled MC collision;purity of PV contributors", kTProfile, {axisNch08Low}); + histos.add("Purity/pPurityAssocVsNch", "reconstructed collisions;N_{ch} (|#eta|<0.8) of the labelled MC collision;purity of associated tracks", kTProfile, {axisNch08Low}); + histos.add("Purity/hPurityPVVsNch", "reconstructed collisions: PV contributors", kTH2F, {axisNch08, axisPurity}); + histos.add("Purity/hPurityAssocVsNch", "reconstructed collisions: associated tracks", kTH2F, {axisNch08, axisPurity}); + histos.add("Purity/pPurityPVVsDist", "reconstructed collisions: PV contributors;distance to the nearest other MC collision (BCs);N_{ch} class", kTProfile2D, {axisDist, axisNchClass}); + histos.add("Purity/pPurityAssocVsDist", "reconstructed collisions: associated tracks;distance to the nearest other MC collision (BCs);N_{ch} class", kTProfile2D, {axisDist, axisNchClass}); + for (int i = 0; i < NNchClasses; i++) { + const char* label = nullptr; + if (i == 0) { + label = Form("N_{ch} < %d", NchClassLimits[i]); + } else if (i == NNchClasses - 1) { + label = Form("N_{ch} #geq %d", NchClassLimits[i - 1]); + } else { + label = Form("%d #leq N_{ch} < %d", NchClassLimits[i - 1], NchClassLimits[i]); + } + histos.get(HIST("Efficiency/hDistGen"))->GetYaxis()->SetBinLabel(i + 1, label); + histos.get(HIST("Efficiency/hDistGenReco"))->GetYaxis()->SetBinLabel(i + 1, label); + histos.get(HIST("Purity/pPurityPVVsDist"))->GetYaxis()->SetBinLabel(i + 1, label); + histos.get(HIST("Purity/pPurityAssocVsDist"))->GetYaxis()->SetBinLabel(i + 1, label); + } + + // bin labels + const std::array categoryNames = {"correct BC", "BC shared with other MC colls", "wrong BC"}; + for (int i = 0; i < kNBcCategories; i++) { + histos.get(HIST("Reco/hCategory"))->GetXaxis()->SetBinLabel(i + 1, categoryNames[i]); + histos.get(HIST("Reco/hPurity"))->GetXaxis()->SetBinLabel(i + 1, categoryNames[i]); + for (const auto& h : {histos.get(HIST("Reco/hNchVsFT0C")).get(), histos.get(HIST("Reco/hNchFromVertexVsFT0C")).get(), histos.get(HIST("Reco/hNchFT0CFromVertexVsFT0C")).get(), + histos.get(HIST("Reco/hNPVVsFT0C")).get(), histos.get(HIST("Reco/hNGlobalVsFT0C")).get(), histos.get(HIST("Reco/hNGlobalVsNchFromVertex")).get()}) { + h->GetZaxis()->SetBinLabel(i + 1, categoryNames[i]); + } + } + auto hSpecies = histos.get(HIST("MC/Displaced/hSpecies")); + hSpecies->GetXaxis()->SetBinLabel(1, "e^{-}"); + hSpecies->GetXaxis()->SetBinLabel(2, "e^{+}"); + hSpecies->GetXaxis()->SetBinLabel(3, "other"); + auto hCounter = histos.get(HIST("Reco/hCounter")); + const std::array counterNames = {"all", "with MC label", "with found FT0", "sel8", "|vZ| cut", "kNoSameBunchPileup"}; + for (size_t i = 0; i < counterNames.size(); i++) { + hCounter->GetXaxis()->SetBinLabel(i + 1, counterNames[i]); + } + } + + void process(ColsMC const& cols, aod::McCollisions const& mcCols, aod::McParticles const& mcParts, TracksMC const& tracks, aod::BCs const& bcs, aod::FT0s const&) + { + // ### MC truth: charged physical primaries (as in MultModule), split into produced at the MC collision vertex and displaced + const size_t nMcColls = mcCols.size(); + std::vector nch(nMcColls, 0), nchFromVertex(nMcColls, 0), nchFT0C(nMcColls, 0), nchFT0CFromVertex(nMcColls, 0), nch08FromVertex(nMcColls, 0); + for (const auto& particle : mcParts) { + if (!particle.isPhysicalPrimary() || !particle.has_mcCollision()) { + continue; + } + const auto* pdgParticle = pdg->GetParticle(particle.pdgCode()); + if (pdgParticle == nullptr || std::abs(pdgParticle->Charge()) < MinCharge) { + continue; + } + const auto mcCol = particle.mcCollision(); + const float dz = particle.vz() - mcCol.posZ(); + const float dr = std::hypot(particle.vx() - mcCol.posX(), particle.vy() - mcCol.posY()); + const bool isFromVertex = std::abs(dz) < maxDistFromVertex && dr < maxDistFromVertex; + const float eta = particle.eta(); + const bool isInNchAcc = std::abs(eta) < MaxEtaNch; + const bool isInFT0CAcc = eta > MinEtaFT0C && eta < MaxEtaFT0C; + const auto iMc = particle.mcCollisionId(); + nch[iMc] += isInNchAcc; + nchFromVertex[iMc] += isInNchAcc && isFromVertex; + nchFT0C[iMc] += isInFT0CAcc; + nchFT0CFromVertex[iMc] += isInFT0CAcc && isFromVertex; + nch08FromVertex[iMc] += std::abs(eta) < MaxEtaNch08 && isFromVertex; + if (isFromVertex) { + histos.fill(HIST("MC/FromVertex/hGenStatus"), particle.getGenStatusCode()); + continue; + } + histos.fill(HIST("MC/Displaced/hPt"), particle.pt()); + histos.fill(HIST("MC/Displaced/hDeltaZ"), dz); + histos.fill(HIST("MC/Displaced/hR"), std::hypot(particle.vx(), particle.vy())); + histos.fill(HIST("MC/Displaced/hSpecies"), particle.pdgCode() == PDG_t::kElectron ? 0 : (particle.pdgCode() == PDG_t::kPositron ? 1 : 2)); + histos.fill(HIST("MC/Displaced/hGenStatus"), particle.getGenStatusCode()); + } + + // number of MC collisions per BC (same-bunch pile-up in MC) + std::unordered_map numMcCollsInBC; + for (const auto& mcCol : mcCols) { + numMcCollsInBC[mcCol.bc_as().globalBC()]++; + } + // number of reconstructed collisions per MC collision + std::vector numRecoPerMcColl(nMcColls, 0); + for (const auto& col : cols) { + if (col.has_mcCollision()) { + numRecoPerMcColl[col.mcCollisionId()]++; + } + } + + for (const auto& mcCol : mcCols) { + const auto iMc = mcCol.globalIndex(); + const float b = mcCol.impactParameter(); + histos.fill(HIST("MC/hNch"), b, nch[iMc]); + histos.fill(HIST("MC/hNchFromVertex"), b, nchFromVertex[iMc]); + if (nch[iMc] > 0) { + histos.fill(HIST("MC/pFractionDisplacedVsB"), b, 1. - static_cast(nchFromVertex[iMc]) / nch[iMc]); + } + histos.fill(HIST("MC/hNchVsNchFT0C"), nchFT0C[iMc], nch[iMc]); + histos.fill(HIST("MC/hNchVsNchFT0CFromVertex"), nchFT0CFromVertex[iMc], nchFromVertex[iMc]); + histos.fill(HIST("MC/hNumMcCollsInBC"), numMcCollsInBC[mcCol.bc_as().globalBC()]); + histos.fill(HIST("MC/hNumRecoPerMcColl"), numRecoPerMcColl[iMc]); + histos.fill(HIST("MC/hNumRecoPerMcCollVsB"), b, numRecoPerMcColl[iMc]); + } + + // ### vertex reconstruction efficiency and purity + // distance (in BCs) from each MC collision to the nearest other MC collision in this DF; -1 if it can be outside the DF (closer to the DF edges) + int64_t firstBCinDF = std::numeric_limits::max(), lastBCinDF = 0; + for (const auto& bc : bcs) { + firstBCinDF = std::min(firstBCinDF, static_cast(bc.globalBC())); + lastBCinDF = std::max(lastBCinDF, static_cast(bc.globalBC())); + } + std::vector> mcCollBCs; // (global BC, MC collision index) + for (const auto& mcCol : mcCols) { + mcCollBCs.emplace_back(static_cast(mcCol.bc_as().globalBC()), mcCol.globalIndex()); + } + std::sort(mcCollBCs.begin(), mcCollBCs.end()); + std::vector distToNearest(nMcColls, -1); + for (size_t k = 0; k < mcCollBCs.size(); k++) { + const int64_t bc = mcCollBCs[k].first; + int64_t dist = std::numeric_limits::max(); + if (k > 0) { + dist = std::min(dist, bc - mcCollBCs[k - 1].first); + } + if (k + 1 < mcCollBCs.size()) { + dist = std::min(dist, mcCollBCs[k + 1].first - bc); + } + if (dist <= std::min(bc - firstBCinDF, lastBCinDF - bc)) { // the nearest one cannot be outside the DF + distToNearest[mcCollBCs[k].second] = dist; + } + } + + // matching of reconstructed collisions to MC collisions (PV contributors and vZ), and purity of the reconstructed collisions + std::vector isReconstructed(nMcColls, false); + std::unordered_map nPVPerMcColl, nAssocPerMcColl; + for (const auto& col : cols) { + nPVPerMcColl.clear(); + nAssocPerMcColl.clear(); + int nPVAll = 0, nAssocAll = 0; + for (const auto& track : tracks.sliceBy(perCollision, col.globalIndex())) { + const int64_t iMcOfTrack = track.has_mcParticle() ? mcParts.iteratorAt(track.mcParticleId()).mcCollisionId() : -1; + nAssocAll++; + nAssocPerMcColl[iMcOfTrack]++; + if (track.isPVContributor()) { + nPVAll++; + nPVPerMcColl[iMcOfTrack]++; + } + } + // matching: MC collision with most PV contributors + int64_t iMcDominant = -1; + int nPVDominant = 0; + for (const auto& [iMc, n] : nPVPerMcColl) { + if (iMc >= 0 && n > nPVDominant) { + iMcDominant = iMc; + nPVDominant = n; + } + } + if (iMcDominant >= 0 && nPVDominant > minFractionMatch.value * nPVAll && std::abs(col.posZ() - mcCols.iteratorAt(iMcDominant).posZ()) < maxDeltaVzMatch.value) { + isReconstructed[iMcDominant] = true; + } + // purity wrt the labelled MC collision + if (!col.has_mcCollision()) { + continue; + } + const auto iMc = col.mcCollisionId(); + const int nch08 = nch08FromVertex[iMc]; + const int64_t dist = distToNearest[iMc]; + const int iClass = nchClass(nch08); + if (nPVAll > 0) { + const double purityPV = static_cast(nPVPerMcColl[iMc]) / nPVAll; + histos.fill(HIST("Purity/pPurityPVVsNch"), nch08, purityPV); + histos.fill(HIST("Purity/hPurityPVVsNch"), nch08, purityPV); + if (dist >= 0) { + histos.fill(HIST("Purity/pPurityPVVsDist"), dist, iClass, purityPV); + } + } + if (nAssocAll > 0) { + const double purityAssoc = static_cast(nAssocPerMcColl[iMc]) / nAssocAll; + histos.fill(HIST("Purity/pPurityAssocVsNch"), nch08, purityAssoc); + histos.fill(HIST("Purity/hPurityAssocVsNch"), nch08, purityAssoc); + if (dist >= 0) { + histos.fill(HIST("Purity/pPurityAssocVsDist"), dist, iClass, purityAssoc); + } + } + } + // efficiency + for (size_t iMc = 0; iMc < nMcColls; iMc++) { + const int nch08 = nch08FromVertex[iMc]; + histos.fill(HIST("Efficiency/hNchGen"), nch08); + histos.fill(HIST("Efficiency/hNchGenWide"), nch08); + if (isReconstructed[iMc]) { + histos.fill(HIST("Efficiency/hNchGenReco"), nch08); + histos.fill(HIST("Efficiency/hNchGenRecoWide"), nch08); + } + if (distToNearest[iMc] < 0) { + continue; + } + histos.fill(HIST("Efficiency/hDistGen"), distToNearest[iMc], nchClass(nch08)); + if (isReconstructed[iMc]) { + histos.fill(HIST("Efficiency/hDistGenReco"), distToNearest[iMc], nchClass(nch08)); + } + } + + // ### reconstructed collisions + for (const auto& col : cols) { + histos.fill(HIST("Reco/hCounter"), 0); + if (!col.has_mcCollision()) { + continue; + } + histos.fill(HIST("Reco/hCounter"), 1); + if (!col.has_foundFT0()) { + continue; + } + histos.fill(HIST("Reco/hCounter"), 2); + if (requireSel8 && !col.sel8()) { + continue; + } + histos.fill(HIST("Reco/hCounter"), 3); + if (std::abs(col.posZ()) > maxAbsVtxZ) { + continue; + } + histos.fill(HIST("Reco/hCounter"), 4); + if (requireNoSameBunchPileup && !col.selection_bit(kNoSameBunchPileup)) { + continue; + } + histos.fill(HIST("Reco/hCounter"), 5); + + const float ampFT0C = col.foundFT0().sumAmpC(); // as in MultModule: FT0 of the found BC + const auto mcCol = col.mcCollision(); + const auto iMc = col.mcCollisionId(); + const auto foundBC = col.foundBC_as().globalBC(); + const auto trueBC = mcCol.bc_as().globalBC(); + const int numMcInFoundBC = numMcCollsInBC.contains(foundBC) ? numMcCollsInBC[foundBC] : 0; + int category = kCorrectBC; + if (foundBC != trueBC) { + category = kWrongBC; + } else if (numMcInFoundBC > 1) { + category = kSharedBC; + } + + // tracks: PV contributors with ITS, global-like tracks, and PV contributors from the labelled MC collision + int nPV = 0, nGlobal = 0, nPVAll = 0, nPVFromLabelledColl = 0; + for (const auto& track : tracks.sliceBy(perCollision, col.globalIndex())) { + if (!track.isPVContributor()) { + continue; + } + nPVAll++; + if (track.has_mcParticle() && mcParts.iteratorAt(track.mcParticleId()).mcCollisionId() == iMc) { + nPVFromLabelledColl++; + } + if (track.pt() < minPt || std::abs(track.eta()) > maxAbsEta || track.itsNCls() < minITSNCls) { + continue; + } + nPV++; + if (track.hasTPC() && track.tpcNClsFound() >= minTPCNCls) { + nGlobal++; + } + } + + histos.fill(HIST("Reco/hCategory"), category); + histos.fill(HIST("Reco/hDeltaBC"), static_cast(foundBC) - static_cast(trueBC)); + histos.fill(HIST("Reco/hNumMcCollsInFoundBC"), numMcInFoundBC); + if (nPVAll > 0) { + histos.fill(HIST("Reco/hPurity"), category, static_cast(nPVFromLabelledColl) / nPVAll); + } + histos.fill(HIST("Reco/hNchVsFT0C"), ampFT0C, nch[iMc], category); + histos.fill(HIST("Reco/hNchFromVertexVsFT0C"), ampFT0C, nchFromVertex[iMc], category); + histos.fill(HIST("Reco/hNchFT0CFromVertexVsFT0C"), ampFT0C, nchFT0CFromVertex[iMc], category); + histos.fill(HIST("Reco/hNPVVsFT0C"), ampFT0C, nPV, category); + histos.fill(HIST("Reco/hNGlobalVsFT0C"), ampFT0C, nGlobal, category); + histos.fill(HIST("Reco/hNGlobalVsNchFromVertex"), nchFromVertex[iMc], nGlobal, category); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +}