Skip to content

PWGHF: Add efficiency and fix impact parameter calculation #5052

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ginnocen merged 8 commits into from
Jan 18, 2021
Merged

PWGHF: Add efficiency and fix impact parameter calculation #5052

Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
a2a6b90
Add physical primary requirements
hzanoli Dec 11, 2020
3e1b1aa
Fix order mc and add particle spefic eff
hzanoli Dec 11, 2020
800fef9
Particle dependent efficiency
hzanoli Dec 14, 2020
d58b0a1
Fix pt bins
hzanoli Dec 17, 2020
99bb72d
Unify Pt bins
hzanoli Dec 18, 2020
7401ce4
Resolve conflics for histogram names
hzanoli Dec 18, 2020
f2e69f9
Fix another histogram name
hzanoli Dec 21, 2020
f9a3250
Merge branch 'dev' into qa-efficiency.
hzanoli Dec 21, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
201 changes: 152 additions & 49 deletions Analysis/Tasks/PWGHF/qaTask.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,15 +18,17 @@
#include "Framework/HistogramRegistry.h"
#include "AnalysisCore/trackUtilities.h"
#include "ReconstructionDataFormats/DCA.h"
#include "AnalysisDataModel/TrackSelectionTables.h"
#include "AnalysisCore/MC.h"
#include "TPDGCode.h"

#include "TH1D.h"

#include <cmath>
#include <string>
#include <sstream>
#include <algorithm>
#include "boost/algorithm/string.hpp"

using namespace o2::framework;
namespace o2fw = o2::framework;
namespace o2exp = o2::framework::expressions;
namespace o2df = o2::dataformats;
Expand Down Expand Up @@ -77,6 +79,7 @@ bool GetImpactParameterAndError(const Track& track, const o2df::VertexBase& prim
}
return propagate;
}

} // namespace track_utils

namespace o2::qa::features
Expand Down Expand Up @@ -140,7 +143,7 @@ class Feature
std::string MakeTitle(std::vector<std::string> axisTitles, const std::string& counts = "Counts")
{
axisTitles.push_back(counts);
return ";" + boost::algorithm::join(axisTitles, ";");
return "; " + boost::algorithm::join(axisTitles, "; ");
}

Feature Eta("#eta");
Expand All @@ -154,48 +157,97 @@ Feature ImpactParameterRPhi("Impact Parameter r#varphi", "#mum");
Feature ImpactParameterRPhiError("Impact Parameter Error r#varphi", "#mum");
Feature ImpactParameterZ("Impact Parameter Z", "#mum");
Feature ImpactParameterZError("Impact Parameter Z Error", "#mum");
Feature NumberOfContributors("Number Of contributors to the PV.");
Feature CovarianceXX("Cov_{xx}", "cm^{2}");
Feature CovarianceXY("Cov_{xy}", "cm^{2}");
Feature CovarianceXZ("Cov_{xz}", "cm^{2}");
Feature CovarianceYY("Cov_{yy}", "cm^{2}");
Feature CovarianceYZ("Cov_{yz}", "cm^{2}");
Feature CovarianceZZ("Cov_{zz}", "cm^{2}");
Feature VertexChi2("#Chi^{2}");

std::vector<double> PtBins = {
0.01, 0.0101, 0.0102, 0.0103, 0.0104, 0.0105, 0.0106, 0.0107, 0.0108, 0.0109, 0.011,
0.0111, 0.0112, 0.0113, 0.0114, 0.0115, 0.0116, 0.0117, 0.0118, 0.0119, 0.012,
0.0121, 0.0122, 0.0123, 0.0124, 0.0125, 0.0126, 0.0127, 0.0128, 0.0129, 0.013,
0.0131, 0.0132, 0.0133, 0.0134, 0.0135, 0.0136, 0.0137, 0.0138, 0.0139, 0.014,
0.0141, 0.0142, 0.0143, 0.0144, 0.0145, 0.0146, 0.0147, 0.0148, 0.0149, 0.015,
0.0151, 0.0152, 0.0153, 0.0154, 0.0155, 0.0156, 0.0157, 0.0158, 0.0159, 0.016,
0.0161, 0.0162, 0.0163, 0.0164, 0.0165, 0.0166, 0.0167, 0.0168, 0.0169, 0.017,
0.0171, 0.0172, 0.0173, 0.0174, 0.0175, 0.0176, 0.0177, 0.0178, 0.0179, 0.018,
0.0181, 0.0182, 0.0183, 0.0184, 0.0185, 0.0186, 0.0187, 0.0188, 0.0189, 0.019,
0.0191, 0.0192, 0.0193, 0.0194, 0.0195, 0.0196, 0.0197, 0.0198, 0.0199, 0.02,
0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.12, 0.14, 0.155, 0.16, 0.165,
0.175, 0.18, 0.185, 0.2, 0.225, 0.25, 0.275, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6,
0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 6.0,
8.0, 10.0, 15.0, 20.0, 30.0, 50.0, 100.0};

} // namespace o2::qa::features

namespace qafeat = o2::qa::features;

/// Task to QA global observables of the event
struct QAGlobalObservables {
o2fw::Configurable<int> nBinsNumberOfTracks{"nBinsNumberOfTracks", 2000, "Number of bins for the Number of Tracks"};
o2fw::Configurable<int> nBinsVertexPosition{"nBinsPt", 100, "Number of bins for the Vertex Position"};

o2fw::Configurable<int> nBinsVertexPosition{"nBinsVertexPosition", 100, "Number of bins for the Vertex Position"};

o2fw::Configurable<int> nBinsNumberOfContributorsVertex{
"nBinsNumberOfContributorsVertex", 200, "Number bins for the number of contributors to the primary vertex"};

o2fw::Configurable<int> numberOfContributorsVertexMax{
"numberOfContributorsVertexMax", 200, "Maximum value for the Number of contributors to the primary vertex"};

o2fw::Configurable<int> nBinsVertexCovarianceMatrix{"nBinsVertexCovarianceMatrix", 100,
"Number bins for the vertex covariance matrix"};

std::array<float, 2> collisionZRange = {-20., 20.};
std::array<float, 2> collisionXYRange = {-0.01, 0.01};

std::array<float, 2> numberOfTracksRange = {0, 400};
std::array<float, 2> vertexCovarianceMatrixRange = {-0.1, 0.1};

o2fw::OutputObj<TH1D> eventCount{TH1D("eventCount", qafeat::MakeTitle({"Selected Events"}).c_str(), 2, 0, 2)};
o2fw::HistogramRegistry histograms{"HistogramsGlobalQA"};

void init(o2fw::InitContext&)
{
histograms.add("collision/collisionX", qafeat::MakeTitle({qafeat::VertexX}).c_str(), o2fw::kTH1D,
{{nBinsVertexPosition, collisionXYRange[0], collisionXYRange[1]}});
o2fw::AxisSpec collisionXAxis{nBinsVertexPosition, collisionXYRange[0], collisionXYRange[1]};
o2fw::AxisSpec collisionYAxis{nBinsVertexPosition, collisionXYRange[0], collisionXYRange[1]};
o2fw::AxisSpec collisionZAxis{nBinsVertexPosition, collisionZRange[0], collisionZRange[1]};

o2fw::AxisSpec numberOfContributorsAxis{nBinsNumberOfContributorsVertex, 0, float(numberOfContributorsVertexMax)};
o2fw::AxisSpec numberOfTrackAxis{nBinsNumberOfTracks, numberOfTracksRange[0], numberOfTracksRange[1]};

o2fw::AxisSpec vertexCovarianceMatrixAxis{nBinsVertexCovarianceMatrix, vertexCovarianceMatrixRange[0],
vertexCovarianceMatrixRange[1]};

histograms.add("collision/collisionX", qafeat::MakeTitle({qafeat::VertexX}).c_str(), o2fw::kTH1D, {collisionXAxis});

histograms.add("collision/collisionY", qafeat::MakeTitle({qafeat::VertexY}).c_str(), o2fw::kTH1D, {collisionYAxis});

histograms.add("collision/collisionY", qafeat::MakeTitle({qafeat::VertexY}).c_str(), o2fw::kTH1D,
{{nBinsVertexPosition, collisionXYRange[0], collisionXYRange[1]}});
histograms.add("collision/collisionZ", qafeat::MakeTitle({qafeat::VertexZ}).c_str(), o2fw::kTH1D, {collisionZAxis});

histograms.add("collision/collisionZ", qafeat::MakeTitle({qafeat::VertexZ}).c_str(), o2fw::kTH1D,
{{nBinsVertexPosition, collisionZRange[0], collisionZRange[1]}});
histograms.add("collision/numberOfContributors", qafeat::MakeTitle({qafeat::NumberOfContributors}).c_str(),
o2fw::kTH1D, {numberOfContributorsAxis});

histograms.add("collision/vertexChi2", qafeat::MakeTitle({qafeat::VertexChi2}).c_str(),
o2fw::kTH1D, {{100, 0, 10}});

histograms.add("multiplicity/numberOfTracks", qafeat::MakeTitle({qafeat::TrackMultiplicity}).c_str(), o2fw::kTH1D,
{{nBinsNumberOfTracks, numberOfTracksRange[0], numberOfTracksRange[1]}});

// covariance histograms
histograms.add("Covariance/xx", "xx;Cov_{xx} [cm^{2}]", kTH1D, {{200, -0.1, 0.1}});
histograms.add("Covariance/xy", "xy;Cov_{xy} [cm^{2}]", kTH1D, {{200, -0.1, 0.1}});
histograms.add("Covariance/xz", "xz;Cov_{xz} [cm^{2}]", kTH1D, {{200, -0.1, 0.1}});
histograms.add("Covariance/yy", "yy;Cov_{yy} [cm^{2}]", kTH1D, {{200, -0.1, 0.1}});
histograms.add("Covariance/yz", "yz;Cov_{yz} [cm^{2}]", kTH1D, {{200, -0.1, 0.1}});
histograms.add("Covariance/zz", "zz;Cov_{zz} [cm^{2}]", kTH1D, {{200, -0.1, 0.1}});

// quality histograms
histograms.add("Quality/Chi2", "#Chi^{2};#Chi^{2}", kTH1D, {{100, 0, 10}});
histograms.add("Quality/Contributors", "Contributors;Contributors", kTH1D, {{100, 0, 100}});
{numberOfTrackAxis});

histograms.add("covariance/xx", qafeat::MakeTitle({qafeat::CovarianceXX}).c_str(), o2fw::kTH1D,
{vertexCovarianceMatrixAxis});
histograms.add("covariance/xy", qafeat::MakeTitle({qafeat::CovarianceXY}).c_str(), o2fw::kTH1D,
{vertexCovarianceMatrixAxis});
histograms.add("covariance/xz", qafeat::MakeTitle({qafeat::CovarianceXZ}).c_str(), o2fw::kTH1D,
{vertexCovarianceMatrixAxis});
histograms.add("covariance/yy", qafeat::MakeTitle({qafeat::CovarianceYY}).c_str(), o2fw::kTH1D,
{vertexCovarianceMatrixAxis});
histograms.add("covariance/yz", qafeat::MakeTitle({qafeat::CovarianceYZ}).c_str(), o2fw::kTH1D,
{vertexCovarianceMatrixAxis});
histograms.add("covariance/zz", qafeat::MakeTitle({qafeat::CovarianceZZ}).c_str(), o2fw::kTH1D,
{vertexCovarianceMatrixAxis});
}

void process(const o2::aod::Collision& collision, const o2::aod::Tracks& tracks)
Expand All @@ -205,32 +257,27 @@ struct QAGlobalObservables {
histograms.fill(HIST("collision/collisionY"), collision.posY());
histograms.fill(HIST("collision/collisionZ"), collision.posZ());

histograms.fill(HIST("collision/numberOfContributors"), collision.numContrib());
histograms.fill(HIST("collision/vertexChi2"), collision.chi2());

histograms.fill(HIST("covariance/xx"), collision.covXX());
histograms.fill(HIST("covariance/xy"), collision.covXY());
histograms.fill(HIST("covariance/xz"), collision.covXZ());
histograms.fill(HIST("covariance/yy"), collision.covYY());
histograms.fill(HIST("covariance/yz"), collision.covYZ());
histograms.fill(HIST("covariance/zz"), collision.covZZ());

int nTracks(0);
for (const auto& track : tracks) {
nTracks++;
}

histograms.fill(HIST("multiplicity/numberOfTracks"), nTracks);

// fill covariance variables
histograms.fill(HIST("Covariance/xx"), collision.covXX());
histograms.fill(HIST("Covariance/xy"), collision.covXY());
histograms.fill(HIST("Covariance/xz"), collision.covXZ());
histograms.fill(HIST("Covariance/yy"), collision.covYY());
histograms.fill(HIST("Covariance/yz"), collision.covYZ());
histograms.fill(HIST("Covariance/zz"), collision.covZZ());

// fill quality variables
histograms.fill(HIST("Quality/Chi2"), collision.chi2());
histograms.fill(HIST("Quality/Contributors"), collision.numContrib());
}
};

/// Task to QA the kinematic properties of the tracks
struct QATrackingKine {
o2fw::Configurable<int> nBinsPt{"nBinsPt", 100, "Number of bins for Pt"};
std::array<double, 2> ptRange = {0, 10.};

o2fw::Configurable<int> nBinsPhi{"nBinsPhi", 100, "Number of bins for Phi"};

o2fw::Configurable<int> nBinsEta{"nBinsEta", 100, "Number of bins for the eta histogram."};
Expand All @@ -241,7 +288,7 @@ struct QATrackingKine {
void init(o2fw::InitContext&)
{
histos.add("tracking/pt", qafeat::MakeTitle({qafeat::Pt}).c_str(), o2fw::kTH1D,
{{nBinsPt, ptRange[0], ptRange[1]}});
{qafeat::PtBins});
histos.add("tracking/eta", qafeat::MakeTitle({qafeat::Eta.NameRaw()}).c_str(), o2fw::kTH1D,
{{nBinsEta, etaRange[0], etaRange[1]}});
histos.add("tracking/phi", qafeat::MakeTitle({qafeat::Phi}).c_str(), o2fw::kTH1D, {{nBinsPhi, 0, 2 * M_PI}});
Expand All @@ -257,10 +304,8 @@ struct QATrackingKine {

/// Task to evaluate the tracking resolution (Pt, Eta, Phi and impact parameter)
struct QATrackingResolution {
std::vector<double> ptBins = {0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10, 0.12,
0.14, 0.16, 0.18, 0.2, 0.225, 0.25, 0.275, 0.3, 0.35, 0.4, 0.45, 0.5,
0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0,
3.5, 4., 5., 6., 8., 10., 15., 20., 30., 50., 100.};
o2fw::Configurable<bool> useOnlyPhysicsPrimary{"useOnlyPhysicsPrimary", true,
"Whether to use only physical primary particles for the resolution."};

o2fw::Configurable<int> nBinsEta{"nBinsEta", 60, "Number of bins for the pseudorapidity"};
std::array<double, 2> etaRange = {-3, 3};
Expand All @@ -277,9 +322,9 @@ struct QATrackingResolution {
o2fw::Configurable<int> nBinsDeltaEta{"nBinsDeltaEta", 100, "Number of bins for the pseudorapidity differences"};
std::array<double, 2> deltaEtaRange = {-0.1, 0.1};

o2fw::Configurable<int> nBinsImpactParameter{"nBinsImpactParameter", 1000, "Number of bins for the Impact parameter"};
o2fw::Configurable<int> nBinsImpactParameter{"nBinsImpactParameter", 2000, "Number of bins for the Impact parameter"};

std::array<double, 2> impactParameterRange = {-1500, 1500}; // micrometer
std::array<double, 2> impactParameterRange = {-500, 500}; // micrometer
std::array<double, 2> impactParameterResolutionRange = {0, 1000}; // micrometer

// Registry of histograms
Expand All @@ -289,7 +334,7 @@ struct QATrackingResolution {
{
// Histogram axis definitions

o2fw::AxisSpec ptAxis{ptBins};
o2fw::AxisSpec ptAxis{qafeat::PtBins};
o2fw::AxisSpec deltaPtAxis{nBinsDeltaPt, deltaPtRange[0], deltaPtRange[1]};
o2fw::AxisSpec deltaPtRelativeAxis{nBinsDeltaPt, deltaPtRange[0], deltaPtRange[1]};
o2fw::AxisSpec deltaPtAbsoluteRelativeAxis{nBinsDeltaPt, 0., deltaPtRange[1]};
Expand Down Expand Up @@ -392,6 +437,13 @@ struct QATrackingResolution {
const o2df::VertexBase primaryVertex = getPrimaryVertex(collision);

for (const auto& track : tracks) {

if (useOnlyPhysicsPrimary) {
const auto mcParticle = track.label();
if (!MC::isPhysicalPrimary(mcParticles, mcParticle)) {
continue;
}
}
const double deltaPt = track.label().pt() - track.pt();
histos.fill(HIST("pt/ptDiffMCRec"), deltaPt);

Expand Down Expand Up @@ -437,9 +489,60 @@ struct QATrackingResolution {
}
};

/// Task to QA the efficiency of a particular particle defined by particlePDG
template <PDG_t particlePDG>
struct QATrackingEfficiency {
o2fw::Configurable<int> nBinsEta{"nBinsEta", 30, "Number of bins for the pseudorapidity"};
std::array<double, 2> etaRange = {-3, 3};

o2fw::Configurable<int> nBinsPhi{"nBinsPhi", 20, "Number of bins for Phi"};
std::array<double, 2> phiRange = {0, 2 * M_PI};

o2fw::HistogramRegistry histos{"histogramsTrackingEfficiencyQA"};

void init(o2fw::InitContext&)
{
o2fw::AxisSpec ptAxis{qafeat::PtBins};
o2fw::AxisSpec phiAxis{nBinsPhi, phiRange[0], phiRange[1]};
o2fw::AxisSpec etaAxis{nBinsEta, etaRange[0], etaRange[1]};

histos.add("reconstructedKinematics",
qafeat::MakeTitle({qafeat::Pt.MC(), qafeat::Eta.MC(), qafeat::Phi.MC()}).c_str(),
o2fw::kTH3D, {ptAxis, etaAxis, phiAxis});

histos.add("generatedKinematics",
qafeat::MakeTitle({qafeat::Pt.MC(), qafeat::Eta.MC(), qafeat::Phi.MC()}).c_str(),
o2fw::kTH3D, {ptAxis, etaAxis, phiAxis});
}

void process(const o2::soa::Join<o2::aod::Tracks, o2::aod::McTrackLabels>& tracks,
const o2::aod::McParticles& mcParticles)
{
for (const auto& track : tracks) {
const auto mcParticle = track.label();
if (MC::isPhysicalPrimary(mcParticles, mcParticle) &&
abs(mcParticle.pdgCode()) == particlePDG) {
histos.fill(HIST("reconstructedKinematics"), mcParticle.pt(), mcParticle.eta(), mcParticle.phi());
}
}

for (const auto& mcParticle : mcParticles) {
if (MC::isPhysicalPrimary(mcParticles, mcParticle) &&
abs(mcParticle.pdgCode()) == particlePDG) {
histos.fill(HIST("generatedKinematics"), mcParticle.pt(), mcParticle.eta(), mcParticle.phi());
}
}
}
};

o2fw::WorkflowSpec defineDataProcessing(o2fw::ConfigContext const&)
{
return o2fw::WorkflowSpec{o2fw::adaptAnalysisTask<QAGlobalObservables>("qa-global-observables"),
o2fw::adaptAnalysisTask<QATrackingKine>("qa-tracking-kine"),
o2fw::adaptAnalysisTask<QATrackingResolution>("qa-tracking-resolution")};
}
o2fw::adaptAnalysisTask<QATrackingResolution>("qa-tracking-resolution"),
o2fw::adaptAnalysisTask<QATrackingEfficiency<PDG_t::kPiPlus>>("qa-tracking-efficiency-pion"),
o2fw::adaptAnalysisTask<QATrackingEfficiency<PDG_t::kProton>>("qa-tracking-efficiency-proton"),
o2fw::adaptAnalysisTask<QATrackingEfficiency<PDG_t::kElectron>>("qa-tracking-efficiency-electron"),
o2fw::adaptAnalysisTask<QATrackingEfficiency<PDG_t::kMuonMinus>>("qa-tracking-efficiency-muon"),
o2fw::adaptAnalysisTask<QATrackingEfficiency<PDG_t::kKPlus>>("qa-tracking-efficiency-kaon")};
};