Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
10 changes: 10 additions & 0 deletions Analysis/ALICE3/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,16 @@ o2_add_dpl_workflow(alice3-qa-multiplicity
PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore
COMPONENT_NAME Analysis)

o2_add_dpl_workflow(alice3-qa-singleparticle
SOURCES src/alice3-qa-singleparticle.cxx
PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore
COMPONENT_NAME Analysis)

o2_add_dpl_workflow(alice3-lutmaker
SOURCES src/alice3-lutmaker.cxx
PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore
COMPONENT_NAME Analysis)

o2_add_dpl_workflow(alice3-pid-rich-qa
SOURCES src/pidRICHqa.cxx
PUBLIC_LINK_LIBRARIES O2::AnalysisDataModel O2::AnalysisCore O2::ALICE3Analysis
Expand Down
176 changes: 176 additions & 0 deletions Analysis/ALICE3/src/alice3-lutmaker.cxx
Original file line number Diff line number Diff line change
@@ -0,0 +1,176 @@
// Copyright CERN and copyright holders of ALICE O2. This software is
// distributed under the terms of the GNU General Public License v3 (GPL
// Version 3), copied verbatim in the file "COPYING".
//
// See http://alice-o2.web.cern.ch/license for full licensing information.
//
// 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.

/// \author Nicolo' Jacazio <nicolo.jacazio@cern.ch>, CERN
/// \brief Task to extract LUTs for the fast simulation from full simulation
/// \since 27/04/2021

// O2 includes
#include "Framework/AnalysisTask.h"
#include "AnalysisCore/MC.h"
#include "ReconstructionDataFormats/Track.h"

using namespace o2;
using namespace framework;
using namespace framework::expressions;

void customize(std::vector<o2::framework::ConfigParamSpec>& workflowOptions)
{
std::vector<ConfigParamSpec> options{
{"lut-el", VariantType::Int, 1, {"LUT input for the Electron PDG code"}},
{"lut-mu", VariantType::Int, 1, {"LUT input for the Muon PDG code"}},
{"lut-pi", VariantType::Int, 1, {"LUT input for the Pion PDG code"}},
{"lut-ka", VariantType::Int, 1, {"LUT input for the Kaon PDG code"}},
{"lut-pr", VariantType::Int, 1, {"LUT input for the Proton PDG code"}},
{"lut-tr", VariantType::Int, 0, {"LUT input for the Triton PDG code"}},
{"lut-de", VariantType::Int, 0, {"LUT input for the Deuteron PDG code"}},
{"lut-he", VariantType::Int, 0, {"LUT input for the Helium3 PDG code"}}};
std::swap(workflowOptions, options);
}

#include "Framework/runDataProcessing.h"

template <o2::track::pid_constants::ID particle>
struct Alice3LutMaker {
static constexpr int nSpecies = 8;
static constexpr int PDGs[nSpecies] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton, 1000010020, 1000010030, 1000020030};
static_assert(particle < nSpecies && "Maximum of particles reached");
static constexpr int pdg = PDGs[particle];
Configurable<bool> selPrim{"sel-prim", false, "If true selects primaries, if not select all particles"};
Configurable<int> etaBins{"eta-bins", 80, "Number of eta bins"};
Configurable<float> etaMin{"eta-min", -4.f, "Lower limit in eta"};
Configurable<float> etaMax{"eta-max", 4.f, "Upper limit in eta"};
Configurable<int> ptBins{"pt-bins", 200, "Number of pT bins"};
Configurable<float> ptMin{"pt-min", -2.f, "Lower limit in pT"};
Configurable<float> ptMax{"pt-max", 2.f, "Upper limit in pT"};
Configurable<int> logPt{"log-pt", 0, "Flag to use a logarithmic pT axis, in this case the pT limits are the expontents"};
HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject};

void init(InitContext&)
{
const TString commonTitle = Form(" PDG %i", pdg);
AxisSpec axisPt{ptBins, ptMin, ptMax, "#it{p}_{T} GeV/#it{c}"};
if (logPt) {
const double min = axisPt.binEdges[0];
const double width = (axisPt.binEdges[1] - axisPt.binEdges[0]) / axisPt.nBins.value();
axisPt.binEdges.clear();
axisPt.binEdges.resize(0);
for (int i = 0; i < axisPt.nBins.value() + 1; i++) {
axisPt.binEdges.push_back(std::pow(10., min + i * width));
}
axisPt.nBins = std::nullopt;
}
const AxisSpec axisEta{etaBins, etaMin, etaMax, "#it{#eta}"};

histos.add("multiplicity", "Track multiplicity;Tracks per event;Events", kTH1F, {{100, 0, 2000}});
histos.add("pt", "pt" + commonTitle, kTH1F, {axisPt});
histos.add("eta", "eta" + commonTitle, kTH1F, {axisEta});
histos.add("CovMat_cYY", "cYY" + commonTitle, kTProfile2D, {axisPt, axisEta});
histos.add("CovMat_cZY", "cZY" + commonTitle, kTProfile2D, {axisPt, axisEta});
histos.add("CovMat_cZZ", "cZZ" + commonTitle, kTProfile2D, {axisPt, axisEta});
histos.add("CovMat_cSnpY", "cSnpY" + commonTitle, kTProfile2D, {axisPt, axisEta});
histos.add("CovMat_cSnpZ", "cSnpZ" + commonTitle, kTProfile2D, {axisPt, axisEta});
histos.add("CovMat_cSnpSnp", "cSnpSnp" + commonTitle, kTProfile2D, {axisPt, axisEta});
histos.add("CovMat_cTglY", "cTglY" + commonTitle, kTProfile2D, {axisPt, axisEta});
histos.add("CovMat_cTglZ", "cTglZ" + commonTitle, kTProfile2D, {axisPt, axisEta});
histos.add("CovMat_cTglSnp", "cTglSnp" + commonTitle, kTProfile2D, {axisPt, axisEta});
histos.add("CovMat_cTglTgl", "cTglTgl" + commonTitle, kTProfile2D, {axisPt, axisEta});
histos.add("CovMat_c1PtY", "c1PtY" + commonTitle, kTProfile2D, {axisPt, axisEta});
histos.add("CovMat_c1PtZ", "c1PtZ" + commonTitle, kTProfile2D, {axisPt, axisEta});
histos.add("CovMat_c1PtSnp", "c1PtSnp" + commonTitle, kTProfile2D, {axisPt, axisEta});
histos.add("CovMat_c1PtTgl", "c1PtTgl" + commonTitle, kTProfile2D, {axisPt, axisEta});
histos.add("CovMat_c1Pt21Pt2", "c1Pt21Pt2" + commonTitle, kTProfile2D, {axisPt, axisEta});

histos.add("Efficiency", "Efficiency" + commonTitle, kTProfile2D, {axisPt, axisEta});
}

void process(const soa::Join<aod::Tracks, aod::TracksCov, aod::McTrackLabels>& tracks,
const aod::McParticles& mcParticles)
{
std::vector<int64_t> recoTracks(tracks.size());
int ntrks = 0;

for (const auto& track : tracks) {
const auto mcParticle = track.mcParticle();
if (mcParticle.pdgCode() != pdg) {
continue;
}
if (selPrim.value && !MC::isPhysicalPrimary(mcParticles, mcParticle)) { // Requiring is physical primary
continue;
}

recoTracks[ntrks++] = mcParticle.globalIndex();

histos.fill(HIST("pt"), mcParticle.pt());
histos.fill(HIST("eta"), mcParticle.eta());
histos.fill(HIST("CovMat_cYY"), mcParticle.pt(), mcParticle.eta(), track.cYY());
histos.fill(HIST("CovMat_cZY"), mcParticle.pt(), mcParticle.eta(), track.cZY());
histos.fill(HIST("CovMat_cZZ"), mcParticle.pt(), mcParticle.eta(), track.cZZ());
histos.fill(HIST("CovMat_cSnpY"), mcParticle.pt(), mcParticle.eta(), track.cSnpY());
histos.fill(HIST("CovMat_cSnpZ"), mcParticle.pt(), mcParticle.eta(), track.cSnpZ());
histos.fill(HIST("CovMat_cSnpSnp"), mcParticle.pt(), mcParticle.eta(), track.cSnpSnp());
histos.fill(HIST("CovMat_cTglY"), mcParticle.pt(), mcParticle.eta(), track.cTglY());
histos.fill(HIST("CovMat_cTglZ"), mcParticle.pt(), mcParticle.eta(), track.cTglZ());
histos.fill(HIST("CovMat_cTglSnp"), mcParticle.pt(), mcParticle.eta(), track.cTglSnp());
histos.fill(HIST("CovMat_cTglTgl"), mcParticle.pt(), mcParticle.eta(), track.cTglTgl());
histos.fill(HIST("CovMat_c1PtY"), mcParticle.pt(), mcParticle.eta(), track.c1PtY());
histos.fill(HIST("CovMat_c1PtZ"), mcParticle.pt(), mcParticle.eta(), track.c1PtZ());
histos.fill(HIST("CovMat_c1PtSnp"), mcParticle.pt(), mcParticle.eta(), track.c1PtSnp());
histos.fill(HIST("CovMat_c1PtTgl"), mcParticle.pt(), mcParticle.eta(), track.c1PtTgl());
histos.fill(HIST("CovMat_c1Pt21Pt2"), mcParticle.pt(), mcParticle.eta(), track.c1Pt21Pt2());
}
histos.fill(HIST("multiplicity"), ntrks);

for (const auto& mcParticle : mcParticles) {
if (mcParticle.pdgCode() != pdg) {
continue;
}
if (!MC::isPhysicalPrimary(mcParticles, mcParticle)) { // Requiring is physical primary
continue;
}

if (std::find(recoTracks.begin(), recoTracks.end(), mcParticle.globalIndex()) != recoTracks.end()) {
histos.fill(HIST("Efficiency"), mcParticle.pt(), mcParticle.eta(), 1.);
} else {
histos.fill(HIST("Efficiency"), mcParticle.pt(), mcParticle.eta(), 0.);
}
}
}
};

WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
{
WorkflowSpec w;
if (cfgc.options().get<int>("lut-el")) {
w.push_back(adaptAnalysisTask<Alice3LutMaker<o2::track::PID::Electron>>(cfgc, TaskName{"alice3-lutmaker-electron"}));
}
if (cfgc.options().get<int>("lut-mu")) {
w.push_back(adaptAnalysisTask<Alice3LutMaker<o2::track::PID::Muon>>(cfgc, TaskName{"alice3-lutmaker-muon"}));
}
if (cfgc.options().get<int>("lut-pi")) {
w.push_back(adaptAnalysisTask<Alice3LutMaker<o2::track::PID::Pion>>(cfgc, TaskName{"alice3-lutmaker-pion"}));
}
if (cfgc.options().get<int>("lut-ka")) {
w.push_back(adaptAnalysisTask<Alice3LutMaker<o2::track::PID::Kaon>>(cfgc, TaskName{"alice3-lutmaker-kaon"}));
}
if (cfgc.options().get<int>("lut-pr")) {
w.push_back(adaptAnalysisTask<Alice3LutMaker<o2::track::PID::Proton>>(cfgc, TaskName{"alice3-lutmaker-proton"}));
}
if (cfgc.options().get<int>("lut-de")) {
w.push_back(adaptAnalysisTask<Alice3LutMaker<o2::track::PID::Deuteron>>(cfgc, TaskName{"alice3-lutmaker-deuteron"}));
}
if (cfgc.options().get<int>("lut-tr")) {
w.push_back(adaptAnalysisTask<Alice3LutMaker<o2::track::PID::Triton>>(cfgc, TaskName{"alice3-lutmaker-triton"}));
}
if (cfgc.options().get<int>("lut-he")) {
w.push_back(adaptAnalysisTask<Alice3LutMaker<o2::track::PID::Helium3>>(cfgc, TaskName{"alice3-lutmaker-helium3"}));
}
return w;
}
135 changes: 135 additions & 0 deletions Analysis/ALICE3/src/alice3-qa-singleparticle.cxx
Original file line number Diff line number Diff line change
@@ -0,0 +1,135 @@
// Copyright CERN and copyright holders of ALICE O2. This software is
// distributed under the terms of the GNU General Public License v3 (GPL
// Version 3), copied verbatim in the file "COPYING".
//
// See http://alice-o2.web.cern.ch/license for full licensing information.
//
// 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.
/// \author Nicolo' Jacazio <nicolo.jacazio@cern.ch>, CERN

// O2 includes
#include "Framework/AnalysisTask.h"
#include "Framework/runDataProcessing.h"
#include "AnalysisCore/MC.h"
#include "Framework/HistogramRegistry.h"

using namespace o2;
using namespace o2::framework;
using namespace o2::framework::expressions;

struct Alice3SingleParticle {
Configurable<int> PDG{"PDG", 2212, "PDG code of the particle of interest"};
Configurable<int> IsStable{"IsStable", 0, "Flag to check stable particles"};
HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject};
Configurable<int> ptBins{"pt-bins", 500, "Number of pT bins"};
Configurable<float> ptMin{"pt-min", 0.f, "Lower limit in pT"};
Configurable<float> ptMax{"pt-max", 5.f, "Upper limit in pT"};
Configurable<int> etaBins{"eta-bins", 500, "Number of eta bins"};
Configurable<float> etaMin{"eta-min", -3.f, "Lower limit in eta"};
Configurable<float> etaMax{"eta-max", 3.f, "Upper limit in eta"};
Configurable<float> yMin{"y-min", -3.f, "Lower limit in y"};
Configurable<float> yMax{"y-max", 3.f, "Upper limit in y"};
Configurable<int> prodBins{"prod-bins", 100, "Number of production vertex bins"};
Configurable<float> prodMin{"prod-min", -1.f, "Lower limit in production vertex"};
Configurable<float> prodMax{"prod-max", 1.f, "Upper limit in production vertex"};
Configurable<float> charge{"charge", 1.f, "Particle charge to scale the reconstructed momentum"};
Configurable<bool> doPrint{"doPrint", false, "Flag to print debug messages"};

void init(InitContext&)
{
const TString tit = Form("%i", PDG.value);
AxisSpec axisPt{ptBins, ptMin, ptMax};
AxisSpec axisEta{etaBins, etaMin, etaMax};
AxisSpec axisProd{prodBins, prodMin, prodMax};

histos.add("particlePt", "Particle Pt " + tit + ";#it{p}_{T} (GeV/#it{c})", kTH1D, {axisPt});
histos.add("prodVx", "Particle Prod. Vertex X " + tit + ";Prod. Vertex X (cm)", kTH1D, {axisProd});
histos.add("prodVy", "Particle Prod. Vertex Y " + tit + ";Prod. Vertex Y (cm)", kTH1D, {axisProd});
histos.add("prodVz", "Particle Prod. Vertex Z " + tit + ";Prod. Vertex Z (cm)", kTH1D, {axisProd});
histos.add("prodRadius", "Particle Prod. Vertex Radius " + tit + ";Prod. Vertex Radius (cm)", kTH1D, {axisProd});
histos.add("prodVxVsPt", "Particle Prod. Vertex X " + tit + ";#it{p}_{T} (GeV/#it{c});Prod. Vertex X (cm)", kTH2D, {axisPt, axisProd});
histos.add("prodVyVsPt", "Particle Prod. Vertex Y " + tit + ";#it{p}_{T} (GeV/#it{c});Prod. Vertex Y (cm)", kTH2D, {axisPt, axisProd});
histos.add("prodVzVsPt", "Particle Prod. Vertex Z " + tit + ";#it{p}_{T} (GeV/#it{c});Prod. Vertex Z (cm)", kTH2D, {axisPt, axisProd});
histos.add("prodRadiusVsPt", "Particle Prod. Vertex Radius " + tit + ";#it{p}_{T} (GeV/#it{c});Prod. Vertex Radius (cm)", kTH2D, {axisPt, axisProd});
histos.add("prodRadius3DVsPt", "Particle Prod. Vertex Radius XYZ " + tit + ";#it{p}_{T} (GeV/#it{c});Prod. Vertex Radius XYZ (cm)", kTH2D, {axisPt, axisProd});
histos.add("trackPt", "Track Pt " + tit + ";#it{p}_{T} (GeV/#it{c})", kTH1D, {axisPt});
histos.add("particleEta", "Particle Eta " + tit + ";#it{#eta}", kTH1D, {axisEta});
histos.add("trackEta", "Track Eta " + tit + ";#it{#eta}", kTH1D, {axisEta});
histos.add("particleY", "Particle Y " + tit + ";#it{y}", kTH1D, {axisEta});
histos.add("primaries", "Source for primaries " + tit + ";PDG Code", kTH1D, {{100, 0.f, 100.f}});
histos.add("secondaries", "Source for secondaries " + tit + ";PDG Code", kTH1D, {{100, 0.f, 100.f}});
}

void process(const soa::Join<o2::aod::Tracks, o2::aod::McTrackLabels>& tracks,
const aod::McParticles& mcParticles)
{

std::vector<int64_t> ParticlesOfInterest;
for (const auto& mcParticle : mcParticles) {
if (mcParticle.pdgCode() != PDG) {
continue;
}
if (mcParticle.y() < yMin || mcParticle.y() > yMax) {
continue;
}
histos.fill(HIST("particlePt"), mcParticle.pt());
histos.fill(HIST("particleEta"), mcParticle.eta());
histos.fill(HIST("particleY"), mcParticle.y());
histos.fill(HIST("prodVx"), mcParticle.vx());
histos.fill(HIST("prodVy"), mcParticle.vy());
histos.fill(HIST("prodVz"), mcParticle.vz());
histos.fill(HIST("prodRadius"), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy()));
histos.fill(HIST("prodVxVsPt"), mcParticle.pt(), mcParticle.vx());
histos.fill(HIST("prodVyVsPt"), mcParticle.pt(), mcParticle.vy());
histos.fill(HIST("prodVzVsPt"), mcParticle.pt(), mcParticle.vz());
histos.fill(HIST("prodRadiusVsPt"), mcParticle.pt(), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy()));
histos.fill(HIST("prodRadius3DVsPt"), mcParticle.pt(), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy() + mcParticle.vz() * mcParticle.vz()));
ParticlesOfInterest.push_back(mcParticle.globalIndex());
}

for (const auto& track : tracks) {
const auto mcParticle = track.mcParticle();
if (!IsStable) {
if (mcParticle.mother0() < 0) {
continue;
}
auto mother = mcParticles.iteratorAt(mcParticle.mother0());
const auto ParticleIsInteresting = std::find(ParticlesOfInterest.begin(), ParticlesOfInterest.end(), mother.globalIndex()) != ParticlesOfInterest.end();
if (!ParticleIsInteresting) {
continue;
}
if (doPrint) {
LOG(INFO) << "Track " << track.globalIndex() << " comes from a " << mother.pdgCode() << " and is a " << mcParticle.pdgCode();
}
} else {
if (mcParticle.pdgCode() != PDG) {
continue;
}
histos.fill(HIST("trackPt"), track.pt() * charge);
histos.fill(HIST("trackEta"), track.eta());
if (mcParticle.mother0() < 0) {
if (doPrint) {
LOG(INFO) << "Track " << track.globalIndex() << " is a " << mcParticle.pdgCode();
}
continue;
}
auto mother = mcParticles.iteratorAt(mcParticle.mother0());
if (MC::isPhysicalPrimary(mcParticles, mcParticle)) {
histos.get<TH1>(HIST("primaries"))->Fill(Form("%i", mother.pdgCode()), 1.f);
} else {
histos.get<TH1>(HIST("secondaries"))->Fill(Form("%i", mother.pdgCode()), 1.f);
}
if (doPrint) {
LOG(INFO) << "Track " << track.globalIndex() << " is a " << mcParticle.pdgCode() << " and comes from a " << mother.pdgCode() << " and is " << (MC::isPhysicalPrimary(mcParticles, mcParticle) ? "" : "not") << " a primary";
}
}
}
}
};

WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
{
return WorkflowSpec{adaptAnalysisTask<Alice3SingleParticle>(cfgc)};
}
Loading