Skip to content

Commit

Permalink
add analysis of split vertices vs. time + update configurables (Alice…
Browse files Browse the repository at this point in the history 
…O2Group#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
  • Loading branch information
@maciacco
maciacco authored Oct 20, 2023
1 parent fd6551e commit b6be33b
Showing 1 changed file with 89 additions and 6 deletions.
95 changes: 89 additions & 6 deletions PWGLF/Tasks/QC/vertexQA.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,8 @@

#include <cmath>
#include <vector>
#include <deque>
#include <algorithm>

#include "Framework/runDataProcessing.h"
#include "Framework/AnalysisDataModel.h"
Expand All @@ -20,6 +22,24 @@ using namespace o2;
using namespace o2::framework;
using namespace o2::framework::expressions;

using BCcoll = std::pair<aod::BC, aod::Collision>;

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<aod::Collision>(bccoll1);
auto coll2 = std::get<aod::Collision>(bccoll2);
auto bc1 = std::get<aod::BC>(bccoll1);
auto bc2 = std::get<aod::BC>(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",
Expand All @@ -40,16 +60,25 @@ struct vertexQA {

Configurable<int> storeTree{"storeTree", 1000, "Store in tree collisions from BC's with more than 'storeTree' vertices, for in-depth analysis"};

Configurable<long unsigned int> nCollMax{"nCollMax", 20, "Maximum size of collision buffer"};
Configurable<double> nSigmaZ{"nSigmaZ", 1000., "Number of sigmas for z of vertices"};
Configurable<double> nSigmaR{"nSigmaR", 1000., "Number of sigmas for transverse displacement of vertices"};
Configurable<double> nSigmaT{"nSigmaT", 1000., "Number of sigmas for time of vertices"};
Configurable<double> 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};

Expand All @@ -72,8 +101,17 @@ struct vertexQA {
histos.add<TH2>("yTwoVtxHistogram", ";#it{y}_{vtx}^{1} (cm);#it{y}_{vtx}^{2} (cm)", HistType::kTH2F, {xVtxAxis, xVtxAxis});
histos.add<TH2>("zTwoVtxHistogram", ";#it{z}_{vtx}^{1} (cm);#it{z}_{vtx}^{2} (cm)", HistType::kTH2F, {zVtxAxis, zVtxAxis});
histos.add<TH2>("tTwoVtxHistogram", ";#it{t}_{vtx}^{1} (ns);#it{t}_{vtx}^{2} (ns)", HistType::kTH2F, {tVtxAxis, tVtxAxis});

histos.add<TH1>("nVtxTimeSeriesHistogram", ";#it{N}_{vtx}^{rec};Entries", HistType::kTH1F, {nVtxTwoAxis});
histos.add<TH1>("zDistVtxTimeSeriesHistogram", ";#Delta#it{z}_{vtx} (cm);Entries", HistType::kTH1F, {zDiffVtxAxis});
histos.add<TH1>("tDistVtxTimeSeriesHistogram", ";#Delta#it{t}_{vtx} (ns);Entries", HistType::kTH1F, {tDiffVtxAxisExtend});
histos.add<TH2>("tCollTwoVtxTimeSeriesHistogram", ";#Delta#it{t}^{coll}_{vtx,1} (ns);#Delta#it{t}^{coll}_{vtx,2} (ns)", HistType::kTH2F, {tVtxAxis, tVtxAxis});
histos.add<TH2>("zDistVsTDistVtxTimeSeriesHistogram", ";#Delta#it{t}_{vtx} (ns);#Delta#it{z}_{vtx} (cm)", HistType::kTH2F, {tDiffVtxAxisExtend, zDiffVtxAxis});
histos.add<TH2>("nContribTwoVtxTimeSeriesHistogram", ";#it{N}_{vtx,1};#it{N}_{vtx,2}", HistType::kTH2F, {nContribAxis, nContribAxis});
}

std::deque<BCcoll> colls;

void process(aod::BC const& bc, aod::Collisions const& collisions)
{
auto collSize{collisions.size()};
Expand Down Expand Up @@ -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<aod::Collision>(colls.front());
auto coll2 = std::get<aod::Collision>(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<aod::Collision>(colls.front());
auto coll2 = std::get<aod::Collision>(*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) {
Expand Down

0 comments on commit b6be33b

Please sign in to comment.