Skip to content

Update TPC QA #378

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
jgrosseo merged 1 commit into from
Dec 3, 2021
Merged

Update TPC QA #378

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
34 changes: 34 additions & 0 deletions Common/Core/TableHelper.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,34 @@
// 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 TableHelper.h
/// \author Nicolò Jacazio nicolo.jacazio@cern.ch
/// \brief Base of utilities to build advanced tasks
///

#include "Framework/InitContext.h"
#include "Framework/RunningWorkflowInfo.h"
#include <string>

/// Function to check if a table is required in a workflow
bool isTableRequiredInWorkflow(o2::framework::InitContext& initContext, const std::string table)
{
auto& workflows = initContext.services().get<o2::framework::RunningWorkflowInfo const>();
for (auto device : workflows.devices) {
for (auto input : device.inputs) {
if (input.matcher.binding == table) {
return true;
}
}
}
return false;
}
190 changes: 127 additions & 63 deletions Common/TableProducer/PID/pidTPC.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@

///
/// \file pidTPC.cxx
/// \author Nicolo' Jacazio
/// \author Nicolò Jacazio
/// \brief Task to produce PID tables for TPC split for each particle with only the Nsigma information.
/// Only the tables for the mass hypotheses requested are filled, the others are sent empty.
///
Expand All @@ -25,6 +25,9 @@
#include "Common/Core/PID/PIDResponse.h"
#include "Common/Core/PID/PIDTPC.h"
#include "Common/DataModel/TrackSelectionTables.h"
#include "Common/DataModel/EventSelection.h"
#include "TableHelper.h"
#include "Framework/StaticFor.h"

using namespace o2;
using namespace o2::framework;
Expand All @@ -40,6 +43,7 @@ void customize(std::vector<o2::framework::ConfigParamSpec>& workflowOptions)

#include "Framework/runDataProcessing.h"

/// Task to produce the TPC response table
struct tpcPid {
using Trks = soa::Join<aod::Tracks, aod::TracksExtra>;
using Coll = aod::Collisions;
Expand Down Expand Up @@ -76,34 +80,31 @@ struct tpcPid {
void init(o2::framework::InitContext& initContext)
{
// Checking the tables are requested in the workflow and enabling them
auto& workflows = initContext.services().get<RunningWorkflowInfo const>();
for (DeviceSpec device : workflows.devices) {
for (auto input : device.inputs) {
auto enableFlag = [&input](const std::string particle, Configurable<int>& flag) {
const std::string table = "pidTPC" + particle;
if (input.matcher.binding == table) {
if (flag < 0) {
flag.value = 1;
LOG(info) << "Auto-enabling table: " + table;
} else if (flag > 0) {
flag.value = 1;
LOG(info) << "Table enabled: " + table;
} else {
LOG(info) << "Table disabled: " + table;
}
}
};
enableFlag("El", pidEl);
enableFlag("Mu", pidMu);
enableFlag("Pi", pidPi);
enableFlag("Ka", pidKa);
enableFlag("Pr", pidPr);
enableFlag("De", pidDe);
enableFlag("Tr", pidTr);
enableFlag("He", pidHe);
enableFlag("Al", pidAl);
auto enableFlag = [&](const std::string particle, Configurable<int>& flag) {
const std::string table = "pidTPC" + particle;
if (isTableRequiredInWorkflow(initContext, table)) {
if (flag < 0) {
flag.value = 1;
LOG(info) << "Auto-enabling table: " + table;
} else if (flag > 0) {
flag.value = 1;
LOG(info) << "Table enabled: " + table;
} else {
LOG(info) << "Table disabled: " + table;
}
}
}
};

enableFlag("El", pidEl);
enableFlag("Mu", pidMu);
enableFlag("Pi", pidPi);
enableFlag("Ka", pidKa);
enableFlag("Pr", pidPr);
enableFlag("De", pidDe);
enableFlag("Tr", pidTr);
enableFlag("He", pidHe);
enableFlag("Al", pidAl);

// Getting the parametrization parameters
ccdb->setURL(url.value);
ccdb->setTimestamp(timestamp.value);
Expand Down Expand Up @@ -172,12 +173,23 @@ struct tpcPid {
}
};

/// Task to produce the TPC QA plots
struct tpcPidQa {
static constexpr int Np = 9;
static constexpr const char* pT[Np] = {"e", "#mu", "#pi", "K", "p", "d", "t", "^{3}He", "#alpha"};
static constexpr std::string_view hnsigma[Np] = {"nsigma/El", "nsigma/Mu", "nsigma/Pi",
"nsigma/Ka", "nsigma/Pr", "nsigma/De",
"nsigma/Tr", "nsigma/He", "nsigma/Al"};
static constexpr std::string_view hnsigmapt[Np] = {"nsigmapt/El", "nsigmapt/Mu", "nsigmapt/Pi",
"nsigmapt/Ka", "nsigmapt/Pr", "nsigmapt/De",
"nsigmapt/Tr", "nsigmapt/He", "nsigmapt/Al"};
static constexpr std::string_view hnsigmapospt[Np] = {"nsigmapospt/El", "nsigmapospt/Mu", "nsigmapospt/Pi",
"nsigmapospt/Ka", "nsigmapospt/Pr", "nsigmapospt/De",
"nsigmapospt/Tr", "nsigmapospt/He", "nsigmapospt/Al"};
static constexpr std::string_view hnsigmanegpt[Np] = {"nsigmanegpt/El", "nsigmanegpt/Mu", "nsigmanegpt/Pi",
"nsigmanegpt/Ka", "nsigmanegpt/Pr", "nsigmanegpt/De",
"nsigmanegpt/Tr", "nsigmanegpt/He", "nsigmanegpt/Al"};

HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::QAObject};

Configurable<int> logAxis{"logAxis", 0, "Flag to use a log momentum axis"};
Expand All @@ -187,73 +199,125 @@ struct tpcPidQa {
Configurable<int> nBinsNSigma{"nBinsNSigma", 200, "Number of bins for the NSigma"};
Configurable<float> minNSigma{"minNSigma", -10.f, "Minimum NSigma in range"};
Configurable<float> maxNSigma{"maxNSigma", 10.f, "Maximum NSigma in range"};
Configurable<int> applyEvSel{"applyEvSel", 0, "Flag to apply rapidity cut: 0 -> no event selection, 1 -> Run 2 event selection, 2 -> Run 3 event selection"};
Configurable<bool> applyTrackCut{"applyTrackCut", false, "Flag to apply standard track cuts"};
Configurable<bool> applyRapidityCut{"applyRapidityCut", false, "Flag to apply rapidity cut"};

template <uint8_t i>
void addParticleHistos()
void addParticleHistos(const AxisSpec& pAxis, const AxisSpec& ptAxis)
{
AxisSpec pAxis{nBinsP, minP, maxP, "#it{p} (GeV/#it{c})"};
if (logAxis) {
pAxis.makeLogaritmic();
}
const AxisSpec nSigmaAxis{nBinsNSigma, minNSigma, maxNSigma, Form("N_{#sigma}^{TPC}(%s)", pT[i])};

// NSigma
const AxisSpec nSigmaAxis{nBinsNSigma, minNSigma, maxNSigma, Form("N_{#sigma}^{TPC}(%s)", pT[i])};
histos.add(hnsigma[i].data(), Form("N_{#sigma}^{TPC}(%s)", pT[i]), kTH2F, {pAxis, nSigmaAxis});
histos.add(hnsigmapt[i].data(), Form("N_{#sigma}^{TPC}(%s)", pT[i]), kTH2F, {ptAxis, nSigmaAxis});
histos.add(hnsigmapospt[i].data(), Form("N_{#sigma}^{TPC}(%s)", pT[i]), kTH2F, {ptAxis, nSigmaAxis});
histos.add(hnsigmanegpt[i].data(), Form("N_{#sigma}^{TPC}(%s)", pT[i]), kTH2F, {ptAxis, nSigmaAxis});
}

void init(o2::framework::InitContext&)
{

const AxisSpec multAxis{1000, 0.f, 1000.f, "Track multiplicity"};
const AxisSpec vtxZAxis{100, -20, 20, "Vtx_{z} (cm)"};
const AxisSpec pAxisPosNeg{nBinsP, -maxP, maxP, "Signed #it{p} (GeV/#it{c})"};
AxisSpec pAxis{nBinsP, minP, maxP, "#it{p} (GeV/#it{c})"};
AxisSpec ptAxis{nBinsP, minP, maxP, "#it{p}_{T} (GeV/#it{c})"};
if (logAxis) {
ptAxis.makeLogaritmic();
pAxis.makeLogaritmic();
}
const AxisSpec vtxZAxis{100, -20, 20, "Vtx_{z} (cm)"};
const AxisSpec dedxAxis{1000, 0, 1000, "d#it{E}/d#it{x} A.U."};

// Event properties
auto h = histos.add<TH1>("event/evsel", "", kTH1F, {{10, 0.5, 10.5, "Ev. Sel."}});
h->GetXaxis()->SetBinLabel(1, "Events read");
h->GetXaxis()->SetBinLabel(2, "Passed ev. sel.");
h->GetXaxis()->SetBinLabel(3, "Passed mult.");
h->GetXaxis()->SetBinLabel(4, "Passed vtx Z");

histos.add("event/vertexz", "", kTH1F, {vtxZAxis});
histos.add("event/multiplicity", "", kTH1F, {multAxis});
histos.add("event/tpcsignal", "", kTH2F, {pAxis, dedxAxis});
histos.add("event/signedtpcsignal", "", kTH2F, {pAxisPosNeg, dedxAxis});

addParticleHistos<0>();
addParticleHistos<1>();
addParticleHistos<2>();
addParticleHistos<3>();
addParticleHistos<4>();
addParticleHistos<5>();
addParticleHistos<6>();
addParticleHistos<7>();
addParticleHistos<8>();
static_for<0, 8>([&](auto i) {
addParticleHistos<i>(pAxis, ptAxis);
});
}

template <uint8_t i, typename T>
void fillParticleHistos(const T& t, const float nsigma)
template <o2::track::PID::ID id, typename T>
void fillParticleHistos(const T& t, const float& nsigma)
{
histos.fill(HIST(hnsigma[i]), t.p(), nsigma);
if (applyRapidityCut) {
const float y = TMath::ASinH(t.pt() / TMath::Sqrt(PID::getMass2(id) + t.pt() * t.pt()) * TMath::SinH(t.eta()));
if (abs(y) > 0.5) {
return;
}
}
// Fill histograms
histos.fill(HIST(hnsigma[id]), t.p(), nsigma);
histos.fill(HIST(hnsigmapt[id]), t.pt(), nsigma);
if (t.sign() > 0) {
histos.fill(HIST(hnsigmapospt[id]), t.pt(), nsigma);
} else {
histos.fill(HIST(hnsigmanegpt[id]), t.pt(), nsigma);
}
}

void process(aod::Collision const& collision, soa::Join<aod::Tracks, aod::TracksExtra,
aod::pidTPCEl, aod::pidTPCMu, aod::pidTPCPi,
aod::pidTPCKa, aod::pidTPCPr, aod::pidTPCDe,
aod::pidTPCTr, aod::pidTPCHe, aod::pidTPCAl,
aod::TrackSelection> const& tracks)
void process(soa::Join<aod::Collisions, aod::EvSels>::iterator const& collision,
soa::Join<aod::Tracks, aod::TracksExtra,
aod::pidTPCEl, aod::pidTPCMu, aod::pidTPCPi,
aod::pidTPCKa, aod::pidTPCPr, aod::pidTPCDe,
aod::pidTPCTr, aod::pidTPCHe, aod::pidTPCAl,
aod::TrackSelection> const& tracks)
{
histos.fill(HIST("event/evsel"), 1);
if (applyEvSel == 1) {
if (!collision.sel7()) {
return;
}
} else if (applyEvSel == 2) {
if (!collision.sel8()) {
return;
}
}
histos.fill(HIST("event/evsel"), 2);

float ntracks = 0;
for (auto t : tracks) {
if (applyTrackCut && !t.isGlobalTrack()) {
continue;
}
ntracks += 1;
}
// if (0 && ntracks < 1) {
// return;
// }
histos.fill(HIST("event/evsel"), 3);
if (abs(collision.posZ()) > 10.f) {
return;
}
histos.fill(HIST("event/evsel"), 4);
histos.fill(HIST("event/vertexz"), collision.posZ());
histos.fill(HIST("event/multiplicity"), ntracks);

for (auto t : tracks) {
if (applyTrackCut && !t.isGlobalTrack()) {
continue;
}
// const float mom = t.p();
const float mom = t.tpcInnerParam();
histos.fill(HIST("event/tpcsignal"), mom, t.tpcSignal());
histos.fill(HIST("event/signedtpcsignal"), mom * t.sign(), t.tpcSignal());
//
fillParticleHistos<0>(t, t.tpcNSigmaEl());
fillParticleHistos<1>(t, t.tpcNSigmaMu());
fillParticleHistos<2>(t, t.tpcNSigmaPi());
fillParticleHistos<3>(t, t.tpcNSigmaKa());
fillParticleHistos<4>(t, t.tpcNSigmaPr());
fillParticleHistos<5>(t, t.tpcNSigmaDe());
fillParticleHistos<6>(t, t.tpcNSigmaTr());
fillParticleHistos<7>(t, t.tpcNSigmaHe());
fillParticleHistos<8>(t, t.tpcNSigmaAl());
fillParticleHistos<PID::Electron>(t, t.tpcNSigmaEl());
fillParticleHistos<PID::Muon>(t, t.tpcNSigmaMu());
fillParticleHistos<PID::Pion>(t, t.tpcNSigmaPi());
fillParticleHistos<PID::Kaon>(t, t.tpcNSigmaKa());
fillParticleHistos<PID::Proton>(t, t.tpcNSigmaPr());
fillParticleHistos<PID::Deuteron>(t, t.tpcNSigmaDe());
fillParticleHistos<PID::Triton>(t, t.tpcNSigmaTr());
fillParticleHistos<PID::Helium3>(t, t.tpcNSigmaHe());
fillParticleHistos<PID::Alpha>(t, t.tpcNSigmaAl());
}
}
};
Expand Down
Loading