diff --git a/Analysis/Tasks/CMakeLists.txt b/Analysis/Tasks/CMakeLists.txt index 44ce185982209..37006fd177484 100644 --- a/Analysis/Tasks/CMakeLists.txt +++ b/Analysis/Tasks/CMakeLists.txt @@ -32,10 +32,6 @@ o2_add_dpl_workflow(trackqa PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore COMPONENT_NAME Analysis) -o2_add_dpl_workflow(qatask - SOURCES qaTask.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore - COMPONENT_NAME Analysis) o2_add_dpl_workflow(pid-tof SOURCES pidTOF.cxx diff --git a/Analysis/Tasks/PWGHF/CMakeLists.txt b/Analysis/Tasks/PWGHF/CMakeLists.txt index 934b0800e2bb2..16e4180937401 100644 --- a/Analysis/Tasks/PWGHF/CMakeLists.txt +++ b/Analysis/Tasks/PWGHF/CMakeLists.txt @@ -8,6 +8,11 @@ # granted to it by virtue of its status as an Intergovernmental Organization or # submit itself to any jurisdiction. +o2_add_dpl_workflow(qatask + SOURCES qaTask.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore + COMPONENT_NAME Analysis) + o2_add_dpl_workflow(hf-track-index-skims-creator SOURCES HFTrackIndexSkimsCreator.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing ROOT::EG diff --git a/Analysis/Tasks/PWGHF/qaTask.cxx b/Analysis/Tasks/PWGHF/qaTask.cxx new file mode 100644 index 0000000000000..b819dcacd6118 --- /dev/null +++ b/Analysis/Tasks/PWGHF/qaTask.cxx @@ -0,0 +1,420 @@ +// 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 Peter Hristov , CERN +/// \author Gian Michele Innocenti , CERN +/// \author Henrique J C Zanoli , Utrecht University +/// \author Nicolo' Jacazio , CERN + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/HistogramRegistry.h" +#include "Analysis/trackUtilities.h" +#include "ReconstructionDataFormats/DCA.h" + +#include "TH1D.h" + +#include +#include +#include +#include "boost/algorithm/string.hpp" + +namespace o2fw = o2::framework; +namespace o2exp = o2::framework::expressions; +namespace o2df = o2::dataformats; + +namespace track_utils +{ +/// Converts a angle phi to the -pi to pi range. +double ConvertPhiRange(double phi) +{ + if (phi > M_PI) { + phi -= 2 * M_PI; + } else if (phi < -M_PI) { + phi += 2 * M_PI; + } + + return phi; +} + +/// Determines the impact parameter and its error for a given track. +/// \param track the track to get the impact parameter from. +/// \param primaryVertex the primary vertex of th collision. +/// \param impactParameterRPhi variable to save the impact parameter (in r phi) in micrometers. +/// \param impactParameterRPhiError variable to save the impact parameter (in r phi) error in micrometers. +/// \param impactParameterZ variable to save the impact parameter (in Z) in micrometers. +/// \param impactParameterZError variable to save the impact parameter (in Z) error in micrometers. +template +bool GetImpactParameterAndError(const Track& track, const o2df::VertexBase& primaryVertex, double& impactParameterRPhi, + double& impactParameterRPhiError, double& impactParameterZ, + double& impactParameterErrorZ) +{ + impactParameterRPhi = -999.; + impactParameterRPhiError = -999.; + impactParameterZ = -999; + impactParameterErrorZ = -999; + + o2df::DCA dca; + // FIXME: get this from CCDB + constexpr float magneticField{5.0}; // in kG + auto trackParameter = getTrackParCov(track); + bool propagate = trackParameter.propagateToDCA(primaryVertex, magneticField, &dca); + + constexpr float conversion_to_micrometer = 1000; + if (propagate) { + impactParameterRPhi = conversion_to_micrometer * dca.getY(); + impactParameterRPhiError = conversion_to_micrometer * std::sqrt(dca.getSigmaY2()); + impactParameterZ = conversion_to_micrometer * dca.getZ(); + impactParameterErrorZ = conversion_to_micrometer * std::sqrt(dca.getSigmaZ2()); + } + return propagate; +} +} // namespace track_utils + +namespace o2::qa::features +{ + +/// Class to abstract the naming of a particular feature. It can help you to build the histogram +/// labels in a consistent way and generate the titles. +class Feature +{ + public: + Feature() = default; + Feature(std::string name, std::string unit = "") : mName(std::move(name)), mUnit(std::move(unit)){}; + + /// Returns the feature tagged as MC with no unit. Example: p_{T}^{MC}. + std::string MCRaw() const { return mName + "^{MC}"; }; + + /// Returns the feature tagged as Reconstructed with no unit. Example: p_{T}^{Rec}. + std::string RecRaw() const { return mName + "^{Rec}"; }; + + /// Returns the unit with no brackets. Example: GeV/c. + std::string UnitRaw() const { return mUnit; }; + + /// Returns the name with no units. Example: p_{T}. + std::string NameRaw() const { return mName; }; + + /// Returns the name with the unit. Example: p_{T} [Gev/c]. + std::string Name() const { return mName + Unit(); }; + + /// Returns the name tagged as MC with the unit. Example: p_{T}^{MC} [Gev/c]. + std::string MC() const { return MCRaw() + Unit(); }; + + /// Returns the name tagged as Reconstructed with the unit. Example: p_{T}^{Rec} [Gev/c] + std::string Rec() const { return RecRaw() + Unit(); }; + + /// Returns the name difference between the MC and reconstructed with the unit. + /// Example: p_{T}^{MC} - p_{T}^{Rec} [Gev/c] + std::string MCRecDiff() const { return MCRaw() + " - " + RecRaw() + Unit(); }; + + /// Returns the name difference between the MC and reconstructed divided by the reconstructed. + /// For example: (p_{T}^{MC} - p_{T}^{Rec})/p_{T}^{Rec} + std::string RelativeMCRecDiff() const { return "(" + MCRaw() + " - " + RecRaw() + ")/" + RecRaw(); }; + + /// Returns the unit formatted to be used in the axis. If no unit is present, returns an empty value. + /// Example: [GeV/c] + std::string Unit() const + { + if (UnitRaw().empty()) { + return ""; + } + return " [" + UnitRaw() + "]"; + }; + + operator std::string() { return Name(); } + + private: + const std::string mName; + const std::string mUnit; +}; + +/// Makes a title for an histogram +std::string MakeTitle(std::vector axisTitles, const std::string& counts = "Counts") +{ + axisTitles.push_back(counts); + return ";" + boost::algorithm::join(axisTitles, ";"); +} + +Feature Eta("#eta"); +Feature TrackMultiplicity("Track Multiplicity"); +Feature Phi{"#varphi", "rad"}; +Feature Pt("p_{T}", "GeV/c"); +Feature VertexX("X", "cm"); +Feature VertexY("Y", "cm"); +Feature VertexZ("Z", "cm"); +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"); +} // namespace o2::qa::features + +namespace qafeat = o2::qa::features; + +/// Task to QA global observables of the event +struct QAGlobalObservables { + o2fw::Configurable nBinsNumberOfTracks{"nBinsNumberOfTracks", 2000, "Number of bins for the Number of Tracks"}; + o2fw::Configurable nBinsVertexPosition{"nBinsPt", 100, "Number of bins for the Vertex Position"}; + + std::array collisionZRange = {-20., 20.}; + std::array collisionXYRange = {-0.01, 0.01}; + + std::array numberOfTracksRange = {0, 400}; + + o2fw::OutputObj 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]}}); + + 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, + {{nBinsVertexPosition, collisionZRange[0], collisionZRange[1]}}); + + histograms.add("multiplicity/numberOfTracks", qafeat::MakeTitle({qafeat::TrackMultiplicity}).c_str(), o2fw::kTH1D, + {{nBinsNumberOfTracks, numberOfTracksRange[0], numberOfTracksRange[1]}}); + } + + void process(const o2::aod::Collision& collision, const o2::aod::Tracks& tracks) + { + eventCount->Fill(0); + histograms.fill("collision/collisionX", collision.posX()); + histograms.fill("collision/collisionY", collision.posY()); + histograms.fill("collision/collisionZ", collision.posZ()); + + int nTracks(0); + for (const auto& track : tracks) { + nTracks++; + } + + histograms.fill("multiplicity/numberOfTracks", nTracks); + } +}; + +/// Task to QA the kinematic properties of the tracks +struct QATrackingKine { + o2fw::Configurable nBinsPt{"nBinsPt", 100, "Number of bins for Pt"}; + std::array ptRange = {0, 10.}; + + o2fw::Configurable nBinsPhi{"nBinsPhi", 100, "Number of bins for Phi"}; + + o2fw::Configurable nBinsEta{"nBinsEta", 100, "Number of bins for the eta histogram."}; + std::array etaRange = {-6, 6}; + + o2fw::HistogramRegistry histos{"HistogramsKineQA"}; + + void init(o2fw::InitContext&) + { + histos.add("tracking/pt", qafeat::MakeTitle({qafeat::Pt}).c_str(), o2fw::kTH1D, + {{nBinsPt, ptRange[0], ptRange[1]}}); + 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}}); + } + + void process(const o2::aod::Track& track) + { + histos.fill("tracking/eta", track.eta()); + histos.fill("tracking/pt", track.pt()); + histos.fill("tracking/phi", track.phi()); + } +}; + +/// Task to evaluate the tracking resolution (Pt, Eta, Phi and impact parameter) +struct QATrackingResolution { + std::vector 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 nBinsEta{"nBinsEta", 60, "Number of bins for the pseudorapidity"}; + std::array etaRange = {-3, 3}; + + o2fw::Configurable nBinsPhi{"nBinsPhi", 50, "Number of bins for Phi"}; + std::array phiRange = {0, 2 * M_PI}; + + o2fw::Configurable nBinsDeltaPt{"nBinsDeltaPt", 100, "Number of bins for the transverse momentum differences"}; + std::array deltaPtRange = {-0.5, 0.5}; + + o2fw::Configurable nBinsDeltaPhi{"nBinsDeltaPhi", 100, "Number of bins for the azimuthal angle differences"}; + std::array deltaPhiRange = {-0.1, 0.1}; + + o2fw::Configurable nBinsDeltaEta{"nBinsDeltaEta", 100, "Number of bins for the pseudorapidity differences"}; + std::array deltaEtaRange = {-0.1, 0.1}; + + o2fw::Configurable nBinsImpactParameter{"nBinsImpactParameter", 1000, "Number of bins for the Impact parameter"}; + + std::array impactParameterRange = {-1500, 1500}; // micrometer + std::array impactParameterResolutionRange = {0, 1000}; // micrometer + + // Registry of histograms + o2fw::HistogramRegistry histos{"HistogramsTrackingResolutionQA"}; + + void init(o2fw::InitContext&) + { + // Histogram axis definitions + + o2fw::AxisSpec ptAxis{ptBins}; + o2fw::AxisSpec deltaPtAxis{nBinsDeltaPt, deltaPtRange[0], deltaPtRange[1]}; + o2fw::AxisSpec deltaPtRelativeAxis{nBinsDeltaPt, deltaPtRange[0], deltaPtRange[1]}; + o2fw::AxisSpec deltaPtAbsoluteRelativeAxis{nBinsDeltaPt, 0., deltaPtRange[1]}; + + o2fw::AxisSpec etaAxis{nBinsEta, etaRange[0], etaRange[1]}; + o2fw::AxisSpec deltaEtaAxis{nBinsDeltaEta, deltaEtaRange[0], deltaEtaRange[1]}; + + o2fw::AxisSpec phiAxis{nBinsPhi, phiRange[0], phiRange[1]}; + o2fw::AxisSpec deltaPhiAxis{nBinsDeltaPhi, deltaPhiRange[0], deltaPhiRange[1]}; + + o2fw::AxisSpec impactParRPhiAxis{nBinsImpactParameter, impactParameterRange[0], impactParameterRange[1]}; + o2fw::AxisSpec impactParRPhiErrorAxis{nBinsImpactParameter, impactParameterResolutionRange[0], + impactParameterResolutionRange[1]}; + + o2fw::AxisSpec impactParZAxis{nBinsImpactParameter, impactParameterRange[0], impactParameterRange[1]}; + o2fw::AxisSpec impactParZErrorAxis{nBinsImpactParameter, impactParameterResolutionRange[0], + impactParameterResolutionRange[1]}; + + // Eta + histos.add("eta/etaDiffMCReco", qafeat::MakeTitle({qafeat::Eta.MCRecDiff()}).c_str(), o2fw::kTH1D, {deltaEtaAxis}); + + histos.add("eta/etaDiffMCRecoVsEtaMC", qafeat::MakeTitle({qafeat::Eta.MCRecDiff(), qafeat::Eta.MC()}).c_str(), + o2fw::kTH2D, {deltaEtaAxis, etaAxis}); + + histos.add("eta/etaDiffMCRecoVsEtaReco", qafeat::MakeTitle({qafeat::Eta.MCRecDiff(), qafeat::Eta.Rec()}).c_str(), + o2fw::kTH2D, {deltaEtaAxis, etaAxis}); + + // Phi + histos.add("phi/phiDiffMCRec", qafeat::MakeTitle({qafeat::Phi.MCRecDiff()}).c_str(), o2fw::kTH1D, {deltaPhiAxis}); + + // Pt + histos.add("pt/ptDiffMCRec", qafeat::MakeTitle({qafeat::Pt.MCRecDiff()}).c_str(), o2fw::kTH1D, {deltaPtAxis}); + + histos.add("pt/ptResolution", qafeat::MakeTitle({qafeat::Pt.RelativeMCRecDiff()}).c_str(), o2fw::kTH1D, + {deltaPtRelativeAxis}); + + histos.add("pt/ptResolutionVsPt", qafeat::MakeTitle({qafeat::Pt.Rec(), qafeat::Pt.RelativeMCRecDiff()}).c_str(), + o2fw::kTH2D, {ptAxis, deltaPtAbsoluteRelativeAxis}); + + histos.add("pt/ptResolutionVsEta", qafeat::MakeTitle({qafeat::Eta.Rec(), qafeat::Pt.RelativeMCRecDiff()}).c_str(), + o2fw::kTH2D, {etaAxis, deltaPtAbsoluteRelativeAxis}); + + histos.add("pt/ptResolutionVsPhi", qafeat::MakeTitle({qafeat::Phi.Rec(), qafeat::Pt.RelativeMCRecDiff()}).c_str(), + o2fw::kTH2D, {phiAxis, deltaPtAbsoluteRelativeAxis}); + + // Impact parameters + histos.add("impactParameter/impactParameterRPhiVsPt", + qafeat::MakeTitle({qafeat::Pt.Rec(), qafeat::ImpactParameterRPhi}).c_str(), o2fw::kTH2D, + {ptAxis, impactParRPhiAxis}); + + histos.add("impactParameter/impactParameterRPhiVsEta", + qafeat::MakeTitle({qafeat::Eta.Rec(), qafeat::ImpactParameterRPhi}).c_str(), o2fw::kTH2D, + {etaAxis, impactParRPhiAxis}); + + histos.add("impactParameter/impactParameterRPhiVsPhi", + qafeat::MakeTitle({qafeat::Phi.Rec(), qafeat::ImpactParameterRPhi}).c_str(), o2fw::kTH2D, + {phiAxis, impactParRPhiAxis}); + + histos.add("impactParameter/impactParameterErrorRPhiVsPt", + qafeat::MakeTitle({qafeat::Pt.Rec(), qafeat::ImpactParameterRPhiError}).c_str(), o2fw::kTH2D, + {ptAxis, impactParRPhiErrorAxis}); + + histos.add("impactParameter/impactParameterErrorRPhiVsEta", + qafeat::MakeTitle({qafeat::Eta.Rec(), qafeat::ImpactParameterRPhiError}).c_str(), o2fw::kTH2D, + {etaAxis, impactParRPhiErrorAxis}); + + histos.add("impactParameter/impactParameterErrorRPhiVsPhi", + qafeat::MakeTitle({qafeat::Phi.Rec(), qafeat::ImpactParameterRPhiError}).c_str(), o2fw::kTH2D, + {phiAxis, impactParRPhiErrorAxis}); + + histos.add("impactParameter/impactParameterZVsPt", + qafeat::MakeTitle({qafeat::Pt.Rec(), qafeat::ImpactParameterZ}).c_str(), o2fw::kTH2D, + {ptAxis, impactParZAxis}); + + histos.add("impactParameter/impactParameterZVsEta", + qafeat::MakeTitle({qafeat::Eta.Rec(), qafeat::ImpactParameterZ}).c_str(), o2fw::kTH2D, + {etaAxis, impactParZAxis}); + + histos.add("impactParameter/impactParameterZVsPhi", + qafeat::MakeTitle({qafeat::Phi.Rec(), qafeat::ImpactParameterZ}).c_str(), o2fw::kTH2D, + {phiAxis, impactParZAxis}); + + histos.add("impactParameter/impactParameterErrorZVsPt", + qafeat::MakeTitle({qafeat::Pt.Rec(), qafeat::ImpactParameterZError}).c_str(), o2fw::kTH2D, + {ptAxis, impactParZErrorAxis}); + + histos.add("impactParameter/impactParameterErrorZVsEta", + qafeat::MakeTitle({qafeat::Eta.Rec(), qafeat::ImpactParameterZError}).c_str(), o2fw::kTH2D, + {etaAxis, impactParZErrorAxis}); + + histos.add("impactParameter/impactParameterErrorZVsPhi", + qafeat::MakeTitle({qafeat::Phi.Rec(), qafeat::ImpactParameterZError}).c_str(), o2fw::kTH2D, + {phiAxis, impactParZErrorAxis}); + } + + void process(const o2::soa::Join::iterator& collision, + const o2::soa::Join& tracks, + const o2::aod::McParticles& mcParticles, const o2::aod::McCollisions& mcCollisions) + { + const o2df::VertexBase primaryVertex = getPrimaryVertex(collision); + + for (const auto& track : tracks) { + const double deltaPt = track.label().pt() - track.pt(); + histos.fill("pt/ptDiffMCRec", deltaPt); + + const double deltaPtOverPt = deltaPt / track.pt(); + + histos.fill("pt/ptResolution", deltaPtOverPt); + histos.fill("pt/ptResolutionVsPt", track.pt(), std::abs(deltaPtOverPt)); + histos.fill("pt/ptResolutionVsEta", track.eta(), std::abs(deltaPtOverPt)); + histos.fill("pt/ptResolutionVsPhi", track.phi(), std::abs(deltaPtOverPt)); + + const double deltaEta = track.label().eta() - track.eta(); + histos.fill("eta/etaDiffMCReco", deltaEta); + histos.fill("eta/etaDiffMCRecoVsEtaMC", deltaEta, track.label().eta()); + histos.fill("eta/etaDiffMCRecoVsEtaReco", deltaEta, track.eta()); + + const double deltaPhi = track_utils::ConvertPhiRange(track.label().phi() - track.phi()); + histos.fill("phi/phiDiffMCRec", deltaPhi); + + double impactParameterRPhi{-999.}, impactParameterRPhiError{-999.}; + double impactParameterZ{-999.}, impactParameterErrorZ{-999.}; + + const bool propagate = track_utils::GetImpactParameterAndError( + track, primaryVertex, impactParameterRPhi, impactParameterRPhiError, impactParameterZ, impactParameterErrorZ); + + if (propagate) { + histos.fill("impactParameter/impactParameterRPhiVsPt", track.pt(), impactParameterRPhi); + histos.fill("impactParameter/impactParameterRPhiVsEta", track.eta(), impactParameterRPhi); + histos.fill("impactParameter/impactParameterRPhiVsPhi", track.phi(), impactParameterRPhi); + + histos.fill("impactParameter/impactParameterZVsPt", track.pt(), impactParameterZ); + histos.fill("impactParameter/impactParameterZVsEta", track.eta(), impactParameterZ); + histos.fill("impactParameter/impactParameterZVsPhi", track.phi(), impactParameterZ); + + histos.fill("impactParameter/impactParameterErrorRPhiVsPt", track.pt(), impactParameterRPhiError); + histos.fill("impactParameter/impactParameterErrorRPhiVsEta", track.eta(), impactParameterRPhiError); + histos.fill("impactParameter/impactParameterErrorRPhiVsPhi", track.phi(), impactParameterRPhiError); + + histos.fill("impactParameter/impactParameterErrorZVsPt", track.pt(), impactParameterErrorZ); + histos.fill("impactParameter/impactParameterErrorZVsEta", track.eta(), impactParameterErrorZ); + histos.fill("impactParameter/impactParameterErrorZVsPhi", track.phi(), impactParameterErrorZ); + } + } + } +}; + +o2fw::WorkflowSpec defineDataProcessing(o2fw::ConfigContext const&) +{ + return o2fw::WorkflowSpec{o2fw::adaptAnalysisTask("qa-global-observables"), + o2fw::adaptAnalysisTask("qa-tracking-kine"), + o2fw::adaptAnalysisTask("qa-tracking-resolution")}; +} diff --git a/Analysis/Tasks/qaTask.cxx b/Analysis/Tasks/qaTask.cxx deleted file mode 100644 index b6c2b0a641a16..0000000000000 --- a/Analysis/Tasks/qaTask.cxx +++ /dev/null @@ -1,224 +0,0 @@ -// 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 Peter Hristov , CERN -/// \author Gian Michele Innocenti , CERN -/// \author Henrique J C Zanoli , Utrecht University -/// \author Nicolo' Jacazio , CERN - -#include "Framework/runDataProcessing.h" -#include "Framework/AnalysisTask.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/HistogramRegistry.h" -#include "Analysis/trackUtilities.h" -#include "ReconstructionDataFormats/DCA.h" -#include "Analysis/MC.h" - -#include -#include -#include - -namespace o2fw = o2::framework; -namespace o2exp = o2::framework::expressions; -namespace o2df = o2::dataformats; - -namespace track_utils -{ -/// Converts a angle phi to the -pi to pi range. -double ConvertPhiRange(double phi) -{ - if (phi > M_PI) { - phi -= 2 * M_PI; - } else if (phi < -M_PI) { - phi += 2 * M_PI; - } - - return phi; -} - -/// Determines the impact parameter and its error for a given track. -/// \param track the track to get the impact parameter from. -/// \param primaryVertex the primary vertex of th collision. -/// \param impactParameter variable to save the impact parameter in micrometers. -/// \param impactParameterError variable to save the impact parameter error in micrometers. -template -bool GetImpactParameterAndError(const Track& track, const o2df::VertexBase& primaryVertex, double& impactParameter, - double& impactParameterError) -{ - impactParameter = -999.; - impactParameterError = -999.; - - o2df::DCA dca; - // FIXME: get this from CCDB - constexpr float magneticField{5.0}; // in kG - auto trackParameter = getTrackParCov(track); - bool propagate = trackParameter.propagateToDCA(primaryVertex, magneticField, &dca); - - if (propagate) { - impactParameter = 1000 * dca.getY(); - impactParameterError = 1000 * std::sqrt(dca.getSigmaY2()); - } - return propagate; -} -} // namespace track_utils - -/// Task to QA global observables of the event -struct QAGlobalObservables { - std::array collisionPositionRange = {-20., 20.}; - std::array numberOfTracksRange = {0, 2000}; - - o2fw::Configurable nBinsNumberOfTracks{"nBinsNumberOfTracks", 2000, "Number of bins fot the Number of Tracks"}; - - o2fw::Configurable nBinsVertexPosition{"nBinsPt", 100, "Number of bins for the Vertex Position"}; - - o2fw::OutputObj hCollisionX{ - TH1F("collisionX", "; X [cm];Counts", nBinsVertexPosition, collisionPositionRange[0], collisionPositionRange[1])}; - - o2fw::OutputObj hCollisionY{ - TH1F("collisionY", "; Y [cm];Counts", nBinsVertexPosition, collisionPositionRange[0], collisionPositionRange[1])}; - - o2fw::OutputObj hCollisionZ{ - TH1F("collisionZ", "; Z [cm];Counts", nBinsVertexPosition, collisionPositionRange[0], collisionPositionRange[1])}; - - o2fw::OutputObj hNumberOfTracks{TH1F("NumberOfTracks", "; Number of Tracks;Counts", nBinsNumberOfTracks, - numberOfTracksRange[0], numberOfTracksRange[1])}; - - void process(const o2::aod::Collision& collision, const o2::aod::Tracks& tracks) - { - hCollisionX->Fill(collision.posX()); - hCollisionY->Fill(collision.posY()); - hCollisionZ->Fill(collision.posZ()); - - int nTracks(0); - for (const auto& track : tracks) { - nTracks++; - } - hNumberOfTracks->Fill(nTracks); - } -}; - -/// Task to QA the kinematic properties of the tracks -struct QATrackingKine { - o2fw::Configurable nBinsPt{"nBinsPt", 100, "Number of bins for Pt"}; - std::array ptRange = {0, 10.}; - o2fw::Configurable nBinsPhi{"nBinsPhi", 100, "Number of bins for Phi"}; - o2fw::Configurable nBinsEta{"nBinsEta", 100, "Number of bins for the eta histogram."}; - std::array etaRange = {-6, 6}; - - o2fw::OutputObj hPt{TH1F("pt", ";p_{T} [GeV];Counts", nBinsPt, ptRange[0], ptRange[1])}; - o2fw::OutputObj hEta{TH1F("eta", ";#eta;Counts", nBinsEta, etaRange[0], etaRange[1])}; - o2fw::OutputObj hPhi{TH1F("phi", ";#varphi [rad];Counts", nBinsPhi, 0, 2 * M_PI)}; - - void process(const o2::aod::Track& track) - { - hEta->Fill(track.eta()); - hPhi->Fill(track.phi()); - hPt->Fill(track.pt()); - } -}; - -/// Task to evaluate the tracking resolution (Pt, Eta, Phi and impact parameter) -struct QATrackingResolution { - o2fw::Configurable nBinsPtTrack{"nBinsPtTrack", 100, "Number of bins for the transverse momentum"}; - std::array ptRange = {0, 10.}; - - o2fw::Configurable nBinsEta{"nBinsEta", 400, "Number of bins for the pseudorapidity"}; - std::array etaRange = {-6, 6}; - - o2fw::Configurable nBinsDeltaPt{"nBinsDeltaPt", 400, "Number of bins for the transverse momentum differences"}; - - o2fw::Configurable nBinsDeltaPhi{"nBinsPhi", 100, "Number of bins for the azimuthal angle differences"}; - - o2fw::Configurable nBinsDeltaEta{"nBinsDeltaEta", 100, "Number of bins for the pseudorapidity differences"}; - - o2fw::Configurable nBinsImpactParameter{"nBinsImpactParameter", 1000, "Number of bins for the Impact parameter"}; - - std::array impactParameterRange = {-1500, 1500}; // micrometer - std::array impactParameterResolutionRange = {0, 1000}; // micrometer - - // Registry of histograms - o2fw::HistogramRegistry histos{"Histos", {}, o2fw::OutputObjHandlingPolicy::AnalysisObject}; - - void init(o2fw::InitContext&) - { - // Eta - histos.add("Eta/etaDiffMCReco", ";#eta_{MC} - #eta_{Rec}", o2fw::kTH1F, - {{nBinsDeltaEta, -2, 2}}); - histos.add("Eta/etaDiffMCRecoVsEtaMC", ";#eta_{MC} - #eta_{Rec};#eta_{MC}", o2fw::kTH2F, - {{nBinsDeltaEta, -2, 2}, {nBinsEta, etaRange[0], etaRange[1]}}); - histos.add("Eta/etaDiffMCRecoVsEtaReco", ";#eta_{MC} - #eta_{Rec};#eta_{Rec}", o2fw::kTH2F, - {{nBinsDeltaEta, -2, 2}, {nBinsEta, etaRange[0], etaRange[1]}}); - // Phi - histos.add("Phi/phiDiffMCRec", ";#varphi_{MC} - #varphi_{Rec} [rad]", o2fw::kTH1F, - {{nBinsDeltaPhi, -M_PI, M_PI}}); - // Pt - histos.add("Pt/ptDiffMCRec", ";p_{T}_{MC} - p_{T}_{Rec} [GeV/c]", o2fw::kTH1F, - {{nBinsDeltaPt, -2., 2.}}); - histos.add("Pt/ptResolution", ";(p_{T}_{MC} - p_{T}_{Rec})/(p_{T}_{Rec})", o2fw::kTH1F, - {{nBinsDeltaPt, -2., 2.}}); - histos.add("Pt/ptResolutionVsPt", ";p_{T} [GeV/c];(p_{T}_{MC} - p_{T}_{Rec})/(p_{T}_{Rec})", o2fw::kTH2F, - {{nBinsPtTrack, ptRange[0], ptRange[1]}, {nBinsDeltaPt, 0, 2.}}); - histos.add("Pt/ptResolutionVsEta", ";#eta;(p_{T}_{MC} - p_{T}_{Rec})/(p_{T}_{Rec})", o2fw::kTH2F, - {{nBinsEta, -4., 4.}, {nBinsDeltaPt, -2., 2.}}); - // Impact parameters - histos.add("ImpactParameter/impactParameterVsPt", ";p_{T} [GeV/c];Impact Parameter [{#mu}m]", o2fw::kTH2F, - {{nBinsPtTrack, ptRange[0], ptRange[1]}, {nBinsImpactParameter, impactParameterRange[0], impactParameterRange[1]}}); - histos.add("ImpactParameter/impactParameterVsEta", "#eta;Impact Parameter [{#mu}m]", o2fw::kTH2F, - {{nBinsEta, etaRange[0], etaRange[1]}, {nBinsImpactParameter, impactParameterRange[0], impactParameterRange[1]}}); - histos.add("ImpactParameter/impactParameterErrorVsPt", ";p_{T} [GeV/c];Impact Parameter Error [#mum]", o2fw::kTH2F, - {{nBinsPtTrack, ptRange[0], ptRange[1]}, {nBinsImpactParameter, impactParameterRange[0], impactParameterRange[1]}}); - histos.add("ImpactParameter/impactParameterErrorVsEta", ";#eta;Impact Parameter Error [#mum]", o2fw::kTH2F, - {{nBinsEta, etaRange[0], etaRange[1]}, {nBinsImpactParameter, 0, impactParameterRange[1]}}); - } - - void process(const o2::soa::Join::iterator& collision, - const o2::soa::Join& tracks, - const o2::aod::McParticles& mcParticles, const o2::aod::McCollisions& mcCollisions) - { - const o2df::VertexBase primaryVertex = getPrimaryVertex(collision); - - for (const auto& track : tracks) { - const double deltaPt = track.label().pt() - track.pt(); - histos.fill("Pt/ptDiffMCRec", deltaPt); - - const double deltaPtOverPt = deltaPt / track.pt(); - histos.fill("Pt/ptResolution", deltaPtOverPt); - histos.fill("Pt/ptResolutionVsPt", track.pt(), abs(deltaPtOverPt)); - histos.fill("Pt/ptResolutionVsEta", track.eta(), abs(deltaPtOverPt)); - - const double deltaEta = track.label().eta() - track.eta(); - histos.fill("Eta/etaDiffMCReco", deltaEta); - histos.fill("Eta/etaDiffMCRecoVsEtaMC", deltaEta, track.label().eta()); - histos.fill("Eta/etaDiffMCRecoVsEtaReco", deltaEta, track.eta()); - - const auto deltaPhi = track_utils::ConvertPhiRange(track.label().phi() - track.phi()); - histos.fill("Phi/phiDiffMCRec", deltaPhi); - - double impactParameter{-999.}; - double impactParameterError{-999.}; - - const bool propagate = - track_utils::GetImpactParameterAndError(track, primaryVertex, impactParameter, impactParameterError); - - if (propagate) { - histos.fill("ImpactParameter/impactParameterVsPt", track.pt(), impactParameter); - histos.fill("ImpactParameter/impactParameterVsEta", track.eta(), impactParameter); - histos.fill("ImpactParameter/impactParameterErrorVsPt", track.pt(), impactParameterError); - histos.fill("ImpactParameter/impactParameterErrorVsEta", track.eta(), impactParameterError); - } - } - } -}; - -o2fw::WorkflowSpec defineDataProcessing(o2fw::ConfigContext const&) -{ - return o2fw::WorkflowSpec{o2fw::adaptAnalysisTask("qa-global-observables"), - o2fw::adaptAnalysisTask("qa-tracking-kine"), - o2fw::adaptAnalysisTask("qa-tracking-resolution")}; -}