Skip to content

Commit

Permalink
Add sampling to DCAr calculation, optional pT selection and DCAr vs p…
Browse files Browse the repository at this point in the history 
…T, NCls, eta plots.
  • Loading branch information
Maximilian Korwieser committed Oct 26, 2023
1 parent 05a1811 commit 7fbccb5
Show file tree
Hide file tree
Showing 2 changed files with 72 additions and 19 deletions.
14 changes: 10 additions & 4 deletions Detectors/TPC/qc/include/TPCQC/Tracks.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,11 +69,14 @@ class Tracks

// To set the elementary track cuts
void setTrackCuts(float AbsEta = 1.,
int nClusterCut = 60, float dEdxTot = 20)
int nClusterCut = 60, float dEdxTot = 20, float cutDCAr = 3, float cutPtForDCAr = 1.5, float samplingFractionDCAr = 0.1)
{
mCutAbsEta = AbsEta;
mCutMinnCls = nClusterCut;
mCutMindEdxTot = dEdxTot;
mCutMinDCAr = cutDCAr;
mCutMaxPtDCAr = cutPtForDCAr;
mSamplingFractionDCAr = samplingFractionDCAr;
}

// Just for backward compatibility with crrent QC, temporary, will be removed in the next PR
Expand All @@ -96,9 +99,12 @@ class Tracks
const std::unordered_map<std::string, std::unique_ptr<TH1>>& getMapHist() const { return mMapHist; }

private:
float mCutAbsEta = 1.f; // Eta cut
int mCutMinnCls = 60; // minimum N clusters
float mCutMindEdxTot = 20.f; // dEdxTot min value
float mCutAbsEta = 1.f; // Eta cut
int mCutMinnCls = 60; // minimum N clusters
float mCutMindEdxTot = 20.f; // dEdxTot min value
float mCutMinDCAr = 3.f; // minimum DCAr
float mCutMaxPtDCAr = 1.5f; // maximum pT for DCAr plots
float mSamplingFractionDCAr = 0.1f; // sampling rate for calculation of DCAr
std::unordered_map<std::string, std::unique_ptr<TH1>> mMapHist;
std::vector<TH1F> mHist1D{}; ///< Initialize vector of 1D histograms
std::vector<TH2F> mHist2D{}; ///< Initialize vector of 2D histograms
Expand Down
77 changes: 62 additions & 15 deletions Detectors/TPC/qc/src/Tracks.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@

// root includes
#include "TFile.h"
#include "TRandom3.h"

// o2 includes
#include "DataFormatsTPC/TrackTPC.h"
Expand Down Expand Up @@ -85,6 +86,13 @@ void Tracks::initializeHistograms()
for (const auto type : types) {
mMapHist[fmt::format("hDCAr_{}", type).data()] = std::make_unique<TH2F>(fmt::format("hDCAr_{}", type).data(), fmt::format("DCAr {};phi;DCAr (cm)", type).data(), 360, 0, o2::math_utils::twoPid(), 250, -10., 10.);
}
// DCA vs variables Histograms
mMapHist["hDCArVsPtPos"] = std::make_unique<TH2F>("hDCArVsPtPos", "DCAr; pT; DCAr (cm)", logPtBinning.size() - 1, logPtBinning.data(), 250, -10., 10.);
mMapHist["hDCArVsEtaPos"] = std::make_unique<TH2F>("hDCArVsEtaPos", "DCAr; #eta; DCAr (cm)", 400, -2., 2., 250, -10., 10.);
mMapHist["hDCArVsNClsPos"] = std::make_unique<TH2F>("hDCArVsNClsPos", "DCAr; NClusters; DCAr (cm)", 400, -0.5, 399.5, 250, -10., 10.);
mMapHist["hDCArVsPtNeg"] = std::make_unique<TH2F>("hDCArVsPtNeg", "DCAr; pT; DCAr (cm)", logPtBinning.size() - 1, logPtBinning.data(), 250, -10., 10.);
mMapHist["hDCArVsEtaNeg"] = std::make_unique<TH2F>("hDCArVsEtaNeg", "DCAr; #eta; DCAr (cm)", 400, -2., 2., 250, -10., 10.);
mMapHist["hDCArVsNClsNeg"] = std::make_unique<TH2F>("hDCArVsNClsNeg", "DCAr; NClusters; DCAr (cm)", 400, -0.5, 399.5, 250, -10., 10.);
}
//______________________________________________________________________________
void Tracks::resetHistograms()
Expand Down Expand Up @@ -130,23 +138,62 @@ bool Tracks::processTrack(const o2::tpc::TrackTPC& track)
auto propagator = o2::base::Propagator::Instance(true);
const int type = (track.getQ2Pt() < 0) + 2 * track.hasCSideClustersOnly();
auto dcaHist = mMapHist[fmt::format("hDCAr_{}", types[type]).data()].get();
auto dcaHist_pT_pos = mMapHist["hDCArVsPtPos"].get();
auto dcaHist_eta_pos = mMapHist["hDCArVsEtaPos"].get();
auto dcaHist_nCluster_pos = mMapHist["hDCArVsNClsPos"].get();
auto dcaHist_pT_neg = mMapHist["hDCArVsPtNeg"].get();
auto dcaHist_eta_neg = mMapHist["hDCArVsEtaNeg"].get();
auto dcaHist_nCluster_neg = mMapHist["hDCArVsNClsNeg"].get();

if (propagator->getMatLUT() && propagator->hasMagFieldSet()) {
// ---| fill DCA histos |---
o2::gpu::gpustd::array<float, 2> dca;
const o2::math_utils::Point3D<float> refPoint{0, 0, 0};
o2::track::TrackPar propTrack(track);
if (propagator->propagateToDCABxByBz(refPoint, propTrack, 2.f, o2::base::Propagator::MatCorrType::USEMatCorrLUT, &dca)) {
const auto phi = o2::math_utils::to02PiGen(track.getPhi());
dcaHist->Fill(phi, dca[0]);
}
} else {
static bool reported = false;
if (!reported) {
LOGP(error, "o2::base::Propagator not properly initialized, MatLUT ({}) and / or Field ({}) missing, will not fill DCA histograms", (void*)propagator->getMatLUT(), (void*)propagator->hasMagFieldSet());
reported = true;
// set-up sampling for the DCA calculation
Double_t sampleProb = 1;

if (mSamplingFractionDCAr > 0) { // for now no SEED is given.
TRandom3 randomGenerator;
sampleProb = randomGenerator.Uniform(1);
}

if (sampleProb > (Double_t)(1. - mSamplingFractionDCAr)) {

if (propagator->getMatLUT() && propagator->hasMagFieldSet()) {
// ---| fill DCA histos |---
o2::gpu::gpustd::array<float, 2> dca;
const o2::math_utils::Point3D<float> refPoint{0, 0, 0};
o2::track::TrackPar propTrack(track);
if (propagator->propagateToDCABxByBz(refPoint, propTrack, 2.f, o2::base::Propagator::MatCorrType::USEMatCorrLUT, &dca)) {
const auto phi = o2::math_utils::to02PiGen(track.getPhi());
dcaHist->Fill(phi, dca[0]);
if (sign < 0.) {
dcaHist_pT_neg->Fill(pt, dca[0]);
if (pt < mCutMaxPtDCAr) {
dcaHist_eta_neg->Fill(eta, dca[0]);
dcaHist_nCluster_neg->Fill(nCls, dca[0]);
}
} else {
dcaHist_pT_pos->Fill(pt, dca[0]);
if (pt < mCutMaxPtDCAr) {
dcaHist_eta_pos->Fill(eta, dca[0]);
dcaHist_nCluster_pos->Fill(nCls, dca[0]);
}
}
}
} else {
static bool reported = false;
if (!reported) {
LOGP(error, "o2::base::Propagator not properly initialized, MatLUT ({}) and / or Field ({}) missing, will not fill DCA histograms", (void*)propagator->getMatLUT(), (void*)propagator->hasMagFieldSet());
reported = true;
}
dcaHist->SetTitle(fmt::format("DCAr {} o2::base::Propagator not properly initialized", types[type]).data());
if (sign < 0.) {
dcaHist_pT_neg->SetTitle(fmt::format("DCAr {} o2::base::Propagator not properly initialized", "p_T neg").data());
dcaHist_eta_neg->SetTitle(fmt::format("DCAr {} o2::base::Propagator not properly initialized", "eta neg").data());
dcaHist_nCluster_neg->SetTitle(fmt::format("DCAr {} o2::base::Propagator not properly initialized", "nClusters neg").data());
} else {
dcaHist_pT_pos->SetTitle(fmt::format("DCAr {} o2::base::Propagator not properly initialized", "p_T pos").data());
dcaHist_eta_pos->SetTitle(fmt::format("DCAr {} o2::base::Propagator not properly initialized", "eta pos").data());
dcaHist_nCluster_pos->SetTitle(fmt::format("DCAr {} o2::base::Propagator not properly initialized", "nClusters pos").data());
}
}
dcaHist->SetTitle(fmt::format("DCAr {} o2::base::Propagator not properly initialized", types[type]).data());
}

if (hasASideOnly == 1) {
Expand Down

0 comments on commit 7fbccb5

Please sign in to comment.