From 79eb39d0b38012dcfbe3899c5d62336dde639230 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Thu, 2 Dec 2021 12:06:28 +0100 Subject: [PATCH] Update TPC QA - Use common utilities - Add track cuts in QA - Add event cuts in QA - Extend available QA plots --- Common/Core/TableHelper.h | 34 ++++ Common/TableProducer/PID/pidTPC.cxx | 190 ++++++++++++++-------- Common/TableProducer/PID/pidTPCFull.cxx | 205 ++++++++++++++++-------- 3 files changed, 301 insertions(+), 128 deletions(-) create mode 100644 Common/Core/TableHelper.h diff --git a/Common/Core/TableHelper.h b/Common/Core/TableHelper.h new file mode 100644 index 00000000000..0bbda4b7241 --- /dev/null +++ b/Common/Core/TableHelper.h @@ -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 + +/// 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(); + for (auto device : workflows.devices) { + for (auto input : device.inputs) { + if (input.matcher.binding == table) { + return true; + } + } + } + return false; +} diff --git a/Common/TableProducer/PID/pidTPC.cxx b/Common/TableProducer/PID/pidTPC.cxx index a5e64c711cc..d4ae79c76cd 100644 --- a/Common/TableProducer/PID/pidTPC.cxx +++ b/Common/TableProducer/PID/pidTPC.cxx @@ -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. /// @@ -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; @@ -40,6 +43,7 @@ void customize(std::vector& workflowOptions) #include "Framework/runDataProcessing.h" +/// Task to produce the TPC response table struct tpcPid { using Trks = soa::Join; using Coll = aod::Collisions; @@ -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(); - for (DeviceSpec device : workflows.devices) { - for (auto input : device.inputs) { - auto enableFlag = [&input](const std::string particle, Configurable& 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& 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); @@ -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 logAxis{"logAxis", 0, "Flag to use a log momentum axis"}; @@ -187,73 +199,125 @@ struct tpcPidQa { Configurable nBinsNSigma{"nBinsNSigma", 200, "Number of bins for the NSigma"}; Configurable minNSigma{"minNSigma", -10.f, "Minimum NSigma in range"}; Configurable maxNSigma{"maxNSigma", 10.f, "Maximum NSigma in range"}; + Configurable applyEvSel{"applyEvSel", 0, "Flag to apply rapidity cut: 0 -> no event selection, 1 -> Run 2 event selection, 2 -> Run 3 event selection"}; + Configurable applyTrackCut{"applyTrackCut", false, "Flag to apply standard track cuts"}; + Configurable applyRapidityCut{"applyRapidityCut", false, "Flag to apply rapidity cut"}; template - 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("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(pAxis, ptAxis); + }); } - template - void fillParticleHistos(const T& t, const float nsigma) + template + 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 const& tracks) + void process(soa::Join::iterator const& collision, + soa::Join 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(t, t.tpcNSigmaEl()); + fillParticleHistos(t, t.tpcNSigmaMu()); + fillParticleHistos(t, t.tpcNSigmaPi()); + fillParticleHistos(t, t.tpcNSigmaKa()); + fillParticleHistos(t, t.tpcNSigmaPr()); + fillParticleHistos(t, t.tpcNSigmaDe()); + fillParticleHistos(t, t.tpcNSigmaTr()); + fillParticleHistos(t, t.tpcNSigmaHe()); + fillParticleHistos(t, t.tpcNSigmaAl()); } } }; diff --git a/Common/TableProducer/PID/pidTPCFull.cxx b/Common/TableProducer/PID/pidTPCFull.cxx index 0a814f8bc04..e3fb24f36c7 100644 --- a/Common/TableProducer/PID/pidTPCFull.cxx +++ b/Common/TableProducer/PID/pidTPCFull.cxx @@ -11,7 +11,7 @@ /// /// \file pidTPCFull.cxx -/// \author Nicolo' Jacazio +/// \author Nicolò Jacazio /// \brief Task to produce PID tables for TPC split for each particle. /// Only the tables for the mass hypotheses requested are filled, the others are sent empty. /// @@ -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; @@ -40,6 +43,7 @@ void customize(std::vector& workflowOptions) #include "Framework/runDataProcessing.h" +/// Task to produce the TPC response table struct tpcPidFull { using Trks = soa::Join; using Coll = aod::Collisions; @@ -76,34 +80,31 @@ struct tpcPidFull { void init(o2::framework::InitContext& initContext) { // Checking the tables are requested in the workflow and enabling them - auto& workflows = initContext.services().get(); - for (DeviceSpec device : workflows.devices) { - for (auto input : device.inputs) { - auto enableFlag = [&input](const std::string particle, Configurable& flag) { - const std::string table = "pidTPCFull" + 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& flag) { + const std::string table = "pidTPCFull" + 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); @@ -167,6 +168,7 @@ struct tpcPidFull { } }; +/// Task to produce the TPC QA plots struct tpcPidFullQa { static constexpr int Np = 9; static constexpr const char* pT[Np] = {"e", "#mu", "#pi", "K", "p", "d", "t", "^{3}He", "#alpha"}; @@ -176,9 +178,22 @@ struct tpcPidFullQa { static constexpr std::string_view hexpected_diff[Np] = {"expected_diff/El", "expected_diff/Mu", "expected_diff/Pi", "expected_diff/Ka", "expected_diff/Pr", "expected_diff/De", "expected_diff/Tr", "expected_diff/He", "expected_diff/Al"}; + static constexpr std::string_view hexpsigma[Np] = {"expsigma/El", "expsigma/Mu", "expsigma/Pi", + "expsigma/Ka", "expsigma/Pr", "expsigma/De", + "expsigma/Tr", "expsigma/He", "expsigma/Al"}; 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 logAxis{"logAxis", 0, "Flag to use a log momentum axis"}; @@ -188,18 +203,19 @@ struct tpcPidFullQa { Configurable nBinsDelta{"nBinsDelta", 200, "Number of bins for the Delta"}; Configurable minDelta{"minDelta", -1000.f, "Minimum Delta in range"}; Configurable maxDelta{"maxDelta", 1000.f, "Maximum Delta in range"}; + Configurable nBinsExpSigma{"nBinsExpSigma", 200, "Number of bins for the ExpSigma"}; + Configurable minExpSigma{"minExpSigma", 0.f, "Minimum ExpSigma in range"}; + Configurable maxExpSigma{"maxExpSigma", 200.f, "Maximum ExpSigma in range"}; Configurable nBinsNSigma{"nBinsNSigma", 200, "Number of bins for the NSigma"}; Configurable minNSigma{"minNSigma", -10.f, "Minimum NSigma in range"}; Configurable maxNSigma{"maxNSigma", 10.f, "Maximum NSigma in range"}; + Configurable applyEvSel{"applyEvSel", 0, "Flag to apply rapidity cut: 0 -> no event selection, 1 -> Run 2 event selection, 2 -> Run 3 event selection"}; + Configurable applyTrackCut{"applyTrackCut", false, "Flag to apply standard track cuts"}; + Configurable applyRapidityCut{"applyRapidityCut", false, "Flag to apply rapidity cut"}; template - void addParticleHistos() + void addParticleHistos(const AxisSpec& pAxis, const AxisSpec& ptAxis) { - AxisSpec pAxis{nBinsP, minP, maxP, "#it{p} (GeV/#it{c})"}; - if (logAxis) { - pAxis.makeLogaritmic(); - } - // Exp signal const AxisSpec expAxis{1000, 0, 1000, Form("d#it{E}/d#it{x}_(%s) A.U.", pT[i])}; histos.add(hexpected[i].data(), "", kTH2F, {pAxis, expAxis}); @@ -208,66 +224,125 @@ struct tpcPidFullQa { const AxisSpec deltaAxis{nBinsDelta, minDelta, maxDelta, Form("d#it{E}/d#it{x} - d#it{E}/d#it{x}(%s)", pT[i])}; histos.add(hexpected_diff[i].data(), "", kTH2F, {pAxis, deltaAxis}); + // Exp Sigma + const AxisSpec expSigmaAxis{nBinsExpSigma, minExpSigma, maxExpSigma, Form("Exp_{#sigma}^{TPC}(%s)", pT[i])}; + histos.add(hexpsigma[i].data(), "", kTH2F, {pAxis, expSigmaAxis}); + // NSigma const AxisSpec nSigmaAxis{nBinsNSigma, minNSigma, maxNSigma, Form("N_{#sigma}^{TPC}(%s)", pT[i])}; - histos.add(hnsigma[i].data(), "", kTH2F, {pAxis, nSigmaAxis}); + 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("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(pAxis, ptAxis); + }); } - template - void fillParticleHistos(const T& t, const float mom, const float exp_diff, const float nsigma) + template + void fillParticleHistos(const T& t, const float& mom, const float& exp_diff, const float& expsigma, const float& nsigma) { - histos.fill(HIST(hexpected[i]), mom, t.tpcSignal() - exp_diff); - histos.fill(HIST(hexpected_diff[i]), mom, exp_diff); - 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(hexpected[id]), mom, t.tpcSignal() - exp_diff); + histos.fill(HIST(hexpected_diff[id]), mom, exp_diff); + histos.fill(HIST(hexpsigma[id]), t.p(), expsigma); + 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 const& tracks) + void process(soa::Join::iterator const& collision, + soa::Join 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, mom, t.tpcExpSignalDiffEl(), t.tpcNSigmaEl()); - fillParticleHistos<1>(t, mom, t.tpcExpSignalDiffMu(), t.tpcNSigmaMu()); - fillParticleHistos<2>(t, mom, t.tpcExpSignalDiffPi(), t.tpcNSigmaPi()); - fillParticleHistos<3>(t, mom, t.tpcExpSignalDiffKa(), t.tpcNSigmaKa()); - fillParticleHistos<4>(t, mom, t.tpcExpSignalDiffPr(), t.tpcNSigmaPr()); - fillParticleHistos<5>(t, mom, t.tpcExpSignalDiffDe(), t.tpcNSigmaDe()); - fillParticleHistos<6>(t, mom, t.tpcExpSignalDiffTr(), t.tpcNSigmaTr()); - fillParticleHistos<7>(t, mom, t.tpcExpSignalDiffHe(), t.tpcNSigmaHe()); - fillParticleHistos<8>(t, mom, t.tpcExpSignalDiffAl(), t.tpcNSigmaAl()); + fillParticleHistos(t, mom, t.tpcExpSignalDiffEl(), t.tpcExpSigmaEl(), t.tpcNSigmaEl()); + fillParticleHistos(t, mom, t.tpcExpSignalDiffMu(), t.tpcExpSigmaMu(), t.tpcNSigmaMu()); + fillParticleHistos(t, mom, t.tpcExpSignalDiffPi(), t.tpcExpSigmaPi(), t.tpcNSigmaPi()); + fillParticleHistos(t, mom, t.tpcExpSignalDiffKa(), t.tpcExpSigmaKa(), t.tpcNSigmaKa()); + fillParticleHistos(t, mom, t.tpcExpSignalDiffPr(), t.tpcExpSigmaPr(), t.tpcNSigmaPr()); + fillParticleHistos(t, mom, t.tpcExpSignalDiffDe(), t.tpcExpSigmaDe(), t.tpcNSigmaDe()); + fillParticleHistos(t, mom, t.tpcExpSignalDiffTr(), t.tpcExpSigmaTr(), t.tpcNSigmaTr()); + fillParticleHistos(t, mom, t.tpcExpSignalDiffHe(), t.tpcExpSigmaHe(), t.tpcNSigmaHe()); + fillParticleHistos(t, mom, t.tpcExpSignalDiffAl(), t.tpcExpSigmaAl(), t.tpcNSigmaAl()); } } };