From b6be33b8dc7a960d557f876a9a5747b64731c2f1 Mon Sep 17 00:00:00 2001 From: Mario Ciacco Date: Fri, 20 Oct 2023 02:00:55 +0200 Subject: [PATCH] add analysis of split vertices vs. time + update configurables (#3649) * add analysis of split vertices vs. time + update configurables * implement clang formatting * clang formatting in function to compute time difference * add 'namespace' after the scope --- PWGLF/Tasks/QC/vertexQA.cxx | 95 ++++++++++++++++++++++++++++++++++--- 1 file changed, 89 insertions(+), 6 deletions(-) diff --git a/PWGLF/Tasks/QC/vertexQA.cxx b/PWGLF/Tasks/QC/vertexQA.cxx index 961b8949f0e..9fc60a68436 100644 --- a/PWGLF/Tasks/QC/vertexQA.cxx +++ b/PWGLF/Tasks/QC/vertexQA.cxx @@ -11,6 +11,8 @@ #include #include +#include +#include #include "Framework/runDataProcessing.h" #include "Framework/AnalysisDataModel.h" @@ -20,6 +22,24 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; +using BCcoll = std::pair; + +namespace +{ +constexpr double LHCRFFreq = 400.789e6; +constexpr double LHCBunchSpacingNS = 10 * 1.e9 / LHCRFFreq; +double deltaTimeColl(BCcoll const bccoll1, BCcoll const bccoll2) +{ + auto coll1 = std::get(bccoll1); + auto coll2 = std::get(bccoll2); + auto bc1 = std::get(bccoll1); + auto bc2 = std::get(bccoll2); + int64_t tmpDT = int64_t(bc1.globalBC()) - int64_t(bc2.globalBC()); + double deltaT = tmpDT * LHCBunchSpacingNS + coll1.collisionTime() - coll2.collisionTime(); + return deltaT; +} +} // namespace + namespace o2::aod { DECLARE_SOA_TABLE(VtxQAtable, "AOD", "VTXQATABLE", @@ -40,16 +60,25 @@ struct vertexQA { Configurable storeTree{"storeTree", 1000, "Store in tree collisions from BC's with more than 'storeTree' vertices, for in-depth analysis"}; + Configurable nCollMax{"nCollMax", 20, "Maximum size of collision buffer"}; + Configurable nSigmaZ{"nSigmaZ", 1000., "Number of sigmas for z of vertices"}; + Configurable nSigmaR{"nSigmaR", 1000., "Number of sigmas for transverse displacement of vertices"}; + Configurable nSigmaT{"nSigmaT", 1000., "Number of sigmas for time of vertices"}; + Configurable maxTime{"maxTime", 100000., "Maximum time difference between split vertices in ns"}; + ConfigurableAxis xVtxAxis{"xVtxBins", {100, -0.1f, 0.1f}, "Binning for the vertex x (y) in cm"}; - ConfigurableAxis zVtxAxis{"zVtxBins", {100, -20.f, 20.f}, "Binning for the vertex z in cm"}; - ConfigurableAxis tVtxAxis{"tVtxBins", {100, -50.f, 50.f}, "Binning for the vertex t in ns"}; + ConfigurableAxis zVtxAxis{"zVtxBins", {100, -10.f, 10.f}, "Binning for the vertex z in cm"}; + ConfigurableAxis tVtxAxis{"tVtxBins", {100, -20.f, 20.f}, "Binning for the vertex t in ns"}; - ConfigurableAxis xDiffVtxAxis{"xDiffVtxBins", {100, -0.1f, 0.1f}, "Binning for the vertex x (y) distance in cm"}; - ConfigurableAxis xyDiffVtxAxis{"xyDiffVtxBins", {100, 0.f, 0.1f}, "Binning for the vertex xy distance in cm"}; - ConfigurableAxis zDiffVtxAxis{"zDiffVtxBins", {100, -20.f, 20.f}, "Binning for the vertex z distance in cm"}; - ConfigurableAxis tDiffVtxAxis{"tDiffVtxBins", {100, -50.f, 50.f}, "Binning for the vertex t distance in ns"}; + ConfigurableAxis xDiffVtxAxis{"xDiffVtxBins", {200, -0.1f, 0.1f}, "Binning for the vertex x (y) distance in cm"}; + ConfigurableAxis xyDiffVtxAxis{"xyDiffVtxBins", {200, 0.f, 0.1f}, "Binning for the vertex xy distance in cm"}; + ConfigurableAxis zDiffVtxAxis{"zDiffVtxBins", {200, -10.f, 10.f}, "Binning for the vertex z distance in cm"}; + ConfigurableAxis tDiffVtxAxis{"tDiffVtxBins", {200, 0.f, 50.f}, "Binning for the vertex t distance in ns"}; + ConfigurableAxis tDiffVtxAxisExtend{"tDiffVtxBinsExtend", {1000, 0.f, 100000.f}, "Binning for the vertex t distance in ns, extended range"}; ConfigurableAxis nVtxAxis{"nVtxBins", {11, -0.5, 10.5}, "Binning for the number of reconstructed vertices per BC"}; + ConfigurableAxis nVtxTwoAxis{"nVtxTwoBins", {2, 0.5, 2.5}, "Binning for the number of reconstructed vertices vs. time"}; + ConfigurableAxis nContribAxis{"nContribBins", {1000, 0, 5000}, "Binning for number of contributors to PV"}; HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -72,8 +101,17 @@ struct vertexQA { histos.add("yTwoVtxHistogram", ";#it{y}_{vtx}^{1} (cm);#it{y}_{vtx}^{2} (cm)", HistType::kTH2F, {xVtxAxis, xVtxAxis}); histos.add("zTwoVtxHistogram", ";#it{z}_{vtx}^{1} (cm);#it{z}_{vtx}^{2} (cm)", HistType::kTH2F, {zVtxAxis, zVtxAxis}); histos.add("tTwoVtxHistogram", ";#it{t}_{vtx}^{1} (ns);#it{t}_{vtx}^{2} (ns)", HistType::kTH2F, {tVtxAxis, tVtxAxis}); + + histos.add("nVtxTimeSeriesHistogram", ";#it{N}_{vtx}^{rec};Entries", HistType::kTH1F, {nVtxTwoAxis}); + histos.add("zDistVtxTimeSeriesHistogram", ";#Delta#it{z}_{vtx} (cm);Entries", HistType::kTH1F, {zDiffVtxAxis}); + histos.add("tDistVtxTimeSeriesHistogram", ";#Delta#it{t}_{vtx} (ns);Entries", HistType::kTH1F, {tDiffVtxAxisExtend}); + histos.add("tCollTwoVtxTimeSeriesHistogram", ";#Delta#it{t}^{coll}_{vtx,1} (ns);#Delta#it{t}^{coll}_{vtx,2} (ns)", HistType::kTH2F, {tVtxAxis, tVtxAxis}); + histos.add("zDistVsTDistVtxTimeSeriesHistogram", ";#Delta#it{t}_{vtx} (ns);#Delta#it{z}_{vtx} (cm)", HistType::kTH2F, {tDiffVtxAxisExtend, zDiffVtxAxis}); + histos.add("nContribTwoVtxTimeSeriesHistogram", ";#it{N}_{vtx,1};#it{N}_{vtx,2}", HistType::kTH2F, {nContribAxis, nContribAxis}); } + std::deque colls; + void process(aod::BC const& bc, aod::Collisions const& collisions) { auto collSize{collisions.size()}; @@ -124,6 +162,51 @@ struct vertexQA { histos.fill(HIST("tTwoVtxHistogram"), collPosT[0], collPosT[1]); } + for (auto col : collisions) { + colls.emplace_back(bc, col); + } + + if (colls.size() > nCollMax) { + auto delta = colls.size() - nCollMax; + for (long unsigned int iC{0}; iC < delta; ++iC) { + histos.fill(HIST("nVtxTimeSeriesHistogram"), 1); + colls.pop_front(); + } + } + + auto compareCollisions = [&](BCcoll other) { + auto coll1 = std::get(colls.front()); + auto coll2 = std::get(other); + double testZ = std::abs(coll1.posZ() - coll2.posZ()) / std::sqrt(coll1.covZZ() + coll2.covZZ()); + double distR = std::hypot(coll1.posX() - coll2.posX(), coll1.posY() - coll2.posY()); + double drdxCol = (coll1.posX() - coll2.posX()) / distR; + double drdxOther = -(coll1.posX() - coll2.posX()) / distR; + double drdyCol = (coll1.posY() - coll2.posY()) / distR; + double drdyOther = -(coll1.posY() - coll2.posY()) / distR; + double covR = std::pow(drdxCol, 2) * coll1.covXX() + std::pow(drdxOther, 2) * coll2.covXX() + std::pow(drdyCol, 2) * coll1.covYY() + std::pow(drdxOther, 2) * coll2.covYY() + 2 * drdxCol * drdyCol * coll1.covXY() + 2 * drdxOther * drdyOther * coll2.covXY(); + double testR = distR / std::sqrt(covR); + double deltaT = deltaTimeColl(colls.front(), other); + double testT = std::abs(deltaT) / std::sqrt(std::pow(coll1.collisionTimeRes(), 2) + std::pow(coll2.collisionTimeRes(), 2)); + return (testT < nSigmaT && testZ < nSigmaZ && testR < nSigmaR && std::abs(deltaT) < maxTime); + }; + + if (colls.size() > 1) { + auto id = std::find_if(colls.begin() + 1, colls.end(), compareCollisions); + if (id != colls.end()) { + auto coll1 = std::get(colls.front()); + auto coll2 = std::get(*id); + double deltaT = deltaTimeColl(colls.front(), *id); + histos.fill(HIST("zDistVtxTimeSeriesHistogram"), coll1.posZ() - coll2.posZ()); + histos.fill(HIST("tDistVtxTimeSeriesHistogram"), std::abs(deltaT)); + histos.fill(HIST("tCollTwoVtxTimeSeriesHistogram"), coll1.collisionTime(), coll2.collisionTime()); + histos.fill(HIST("zDistVsTDistVtxTimeSeriesHistogram"), std::abs(deltaT), coll1.posZ() - coll2.posZ()); + histos.fill(HIST("nContribTwoVtxTimeSeriesHistogram"), coll1.numContrib(), coll2.numContrib()); + histos.fill(HIST("nVtxTimeSeriesHistogram"), 2); + colls.erase(id); + colls.pop_front(); + } + } + histos.fill(HIST("nVtxHistogram"), collSize); if (collSize > storeTree) { for (int i{0}; i < collSize; ++i) {