From aef0a1321ce195c6a28af4a84c47447e8159aaca Mon Sep 17 00:00:00 2001 From: AnnalenaKa Date: Mon, 14 Mar 2022 15:25:36 +0100 Subject: [PATCH 1/2] New TPC PID Response after BugFix --- Common/Core/AnalysisCoreLinkDef.h | 3 +- Common/Core/CMakeLists.txt | 6 +- Common/Core/PID/BetheBloch.h | 55 ---- Common/Core/PID/PIDTPC.h | 70 ---- Common/Core/PID/TPCPIDResponse.h | 180 ++++++++++ Common/Core/PID/TPCReso.h | 54 --- Common/TableProducer/PID/CMakeLists.txt | 4 +- Common/TableProducer/PID/pidBayes.cxx | 94 +++--- Common/TableProducer/PID/pidTPC.cxx | 225 ++++++------- Common/TableProducer/PID/pidTPCFull.cxx | 385 +++++++++++++++++----- Common/TableProducer/PID/qaTPCMC.cxx | 4 +- Common/Tools/CMakeLists.txt | 10 +- Common/Tools/handleParamTPCBetheBloch.cxx | 166 ---------- Common/Tools/handleParamTPCResponse.cxx | 300 +++++++++++++++++ PWGEM/Tasks/gammaConversions.cxx | 1 - PWGEM/Tasks/gammaConversionsmc.cxx | 1 - 16 files changed, 938 insertions(+), 620 deletions(-) delete mode 100644 Common/Core/PID/BetheBloch.h delete mode 100644 Common/Core/PID/PIDTPC.h create mode 100644 Common/Core/PID/TPCPIDResponse.h delete mode 100644 Common/Core/PID/TPCReso.h delete mode 100644 Common/Tools/handleParamTPCBetheBloch.cxx create mode 100644 Common/Tools/handleParamTPCResponse.cxx diff --git a/Common/Core/AnalysisCoreLinkDef.h b/Common/Core/AnalysisCoreLinkDef.h index 754a5819ce2..16ae21e0bc5 100644 --- a/Common/Core/AnalysisCoreLinkDef.h +++ b/Common/Core/AnalysisCoreLinkDef.h @@ -19,6 +19,5 @@ #pragma link C++ class o2::pid::Parametrization + ; #pragma link C++ class o2::pid::tof::TOFReso + ; -#pragma link C++ class o2::pid::tpc::BetheBloch + ; -#pragma link C++ class o2::pid::tpc::TPCReso + ; #pragma link C++ class OrbitRange + ; +#pragma link C++ class o2::pid::tpc::Response + ; \ No newline at end of file diff --git a/Common/Core/CMakeLists.txt b/Common/Core/CMakeLists.txt index c9202cc0877..9f3491dcc89 100644 --- a/Common/Core/CMakeLists.txt +++ b/Common/Core/CMakeLists.txt @@ -23,12 +23,10 @@ o2physics_target_root_dictionary(AnalysisCore PID/DetectorResponse.h PID/PIDTOF.h PID/TOFReso.h - PID/PIDTPC.h - PID/BetheBloch.h - PID/TPCReso.h + PID/TPCPIDResponse.h LINKDEF AnalysisCoreLinkDef.h) o2physics_add_executable(merger COMPONENT_NAME aod SOURCES aodMerger.cxx - PUBLIC_LINK_LIBRARIES ROOT::Hist ROOT::Core ROOT::Net) + PUBLIC_LINK_LIBRARIES ROOT::Hist ROOT::Core ROOT::Net) \ No newline at end of file diff --git a/Common/Core/PID/BetheBloch.h b/Common/Core/PID/BetheBloch.h deleted file mode 100644 index d59112fa84a..00000000000 --- a/Common/Core/PID/BetheBloch.h +++ /dev/null @@ -1,55 +0,0 @@ -// 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 BetheBloch.h -/// \author Nicolo' Jacazio -/// \since 07/08/2020 -/// \brief Implementation for the TPC PID response of the BB parametrization -/// - -#ifndef O2_ANALYSIS_PID_BETHEBLOCH_H_ -#define O2_ANALYSIS_PID_BETHEBLOCH_H_ - -#include "DataFormatsTPC/BetheBlochAleph.h" -#include "Common/Core/PID/ParamBase.h" -#include "ReconstructionDataFormats/PID.h" -#include - -namespace o2::pid::tpc -{ - -class BetheBloch : public Parametrization -{ - public: - BetheBloch() : Parametrization("BetheBloch", 7){}; - ~BetheBloch() override = default; - float operator()(const float* x) const override - { - return mParameters[5] * o2::tpc::BetheBlochAleph(x[0], mParameters[0], mParameters[1], mParameters[2], mParameters[3], mParameters[4]) * std::pow(x[1], mParameters[6]); - } - ClassDef(BetheBloch, 1); -}; - -float BetheBlochParam(const float& momentum, const float& mass, const float& charge, const Parameters& parameters) -{ - return parameters[5] * o2::tpc::BetheBlochAleph(momentum / mass, parameters[0], parameters[1], parameters[2], parameters[3], parameters[4]) * std::pow(charge, parameters[6]); -} - -template -float BetheBlochParamTrack(const T& track, const Parameters& parameters) -{ - return BetheBlochParam(track.tpcInnerParam(), o2::track::pid_constants::sMasses[id], (float)o2::track::pid_constants::sCharges[id], parameters); -} - -} // namespace o2::pid::tpc - -#endif diff --git a/Common/Core/PID/PIDTPC.h b/Common/Core/PID/PIDTPC.h deleted file mode 100644 index 74531bbaa97..00000000000 --- a/Common/Core/PID/PIDTPC.h +++ /dev/null @@ -1,70 +0,0 @@ -// 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 PIDTPC.h -/// \author Nicolo' Jacazio -/// \since 2020-07-24 -/// \brief Implementation of the TPC detector response for PID -/// - -#ifndef O2_FRAMEWORK_PIDTPC_H_ -#define O2_FRAMEWORK_PIDTPC_H_ - -// ROOT includes -#include "Rtypes.h" -#include "TMath.h" - -// O2 includes -#include "Framework/Logger.h" -#include "ReconstructionDataFormats/PID.h" -#include "Common/Core/PID/DetectorResponse.h" - -namespace o2::pid::tpc -{ - -/// \brief Class to handle the the TPC detector response -template -class ELoss -{ - public: - ELoss() = default; - ~ELoss() = default; - - /// Gets the expected signal of the measurement - float GetExpectedSignal(const DetectorResponse& response, const TrackType& trk) const; - - /// Gets the expected resolution of the measurement - float GetExpectedSigma(const DetectorResponse& response, const TrackType& trk) const; - - /// Gets the number of sigmas with respect the expected value - float GetSeparation(const DetectorResponse& response, const TrackType& trk) const { return (trk.tpcSignal() - GetExpectedSignal(response, trk)) / GetExpectedSigma(response, trk); } -}; - -template -float ELoss::GetExpectedSignal(const DetectorResponse& response, const TrackType& trk) const -{ - const float x[2] = {trk.tpcInnerParam() / o2::track::PID::getMass(id), (float)o2::track::PID::getCharge(id)}; - const float bethe = response(DetectorResponse::kSignal, x); - return bethe >= 0.f ? bethe : 0.f; -} - -template -float ELoss::GetExpectedSigma(const DetectorResponse& response, const TrackType& trk) const -{ - const float x[2] = {trk.tpcSignal(), (float)trk.tpcNClsFound()}; - const float reso = response(DetectorResponse::kSigma, x); - return reso >= 0.f ? reso : 0.f; -} - -} // namespace o2::pid::tpc - -#endif // O2_FRAMEWORK_PIDTPC_H_ diff --git a/Common/Core/PID/TPCPIDResponse.h b/Common/Core/PID/TPCPIDResponse.h new file mode 100644 index 00000000000..a99c0bb5cea --- /dev/null +++ b/Common/Core/PID/TPCPIDResponse.h @@ -0,0 +1,180 @@ +// 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 TPCPIDResponse.h +/// \author Annalena Kalteyer +/// \author Christian Sonnabend +/// \brief + +#ifndef O2_PID_TPC_RESPONSE_H_ +#define O2_PID_TPC_RESPONSE_H_ + +#include +#include +#include +#include "Framework/Logger.h" +// O2 includes +#include "ReconstructionDataFormats/PID.h" +#include "DataFormatsTPC/BetheBlochAleph.h" + +namespace o2::pid::tpc +{ + +/// \brief Class to handle the TPC PID response + +class Response +{ + + public: + Response() = default; + ~Response() = default; + + /// Setter and Getter for the private parameters + void SetBetheBlochParams(const std::array& betheBlochParams) { mBetheBlochParams = betheBlochParams; } + void SetResolutionParamsDefault(const std::array& resolutionParamsDefault) { mResolutionParamsDefault = resolutionParamsDefault; } + void SetResolutionParams(const std::vector& resolutionParams) { mResolutionParams = resolutionParams; } + void SetMIP(const float mip) { mMIP = mip; } + void SetChargeFactor(const float chargeFactor) { mChargeFactor = chargeFactor; } + void SetMultiplicityNormalization(const float multNormalization) { mMultNormalization = multNormalization; } + void SetUseDefaultResolutionParam(const bool useDefault) { mUseDefaultResolutionParam = useDefault; } + void SetParameters(const Response* response) + { + mBetheBlochParams = response->GetBetheBlochParams(); + mResolutionParamsDefault = response->GetResolutionParamsDefault(); + mResolutionParams = response->GetResolutionParams(); + mMIP = response->GetMIP(); + mChargeFactor = response->GetChargeFactor(); + mMultNormalization = response->GetMultiplicityNormalization(); + mUseDefaultResolutionParam = response->GetUseDefaultResolutionParam(); + } + + const std::array GetBetheBlochParams() const { return mBetheBlochParams; } + const std::array GetResolutionParamsDefault() const { return mResolutionParamsDefault; } + const std::vector GetResolutionParams() const { return mResolutionParams; } + const float GetMIP() const { return mMIP; } + const float GetChargeFactor() const { return mChargeFactor; } + const float GetMultiplicityNormalization() const { return mMultNormalization; } + const bool GetUseDefaultResolutionParam() const { return mUseDefaultResolutionParam; } + + /// Gets the expected signal of the track + template + float GetExpectedSignal(const TrackType& track, const o2::track::PID::ID id) const; + /// Gets the expected resolution of the track + template + float GetExpectedSigma(const CollisionType& collision, const TrackType& trk, const o2::track::PID::ID id) const; + /// Gets the number of sigmas with respect the expected value + template + float GetNumberOfSigma(const CollisionType& collision, const TrackType& trk, const o2::track::PID::ID id) const; + /// Gets the deviation to the expected signal + template + float GetSignalDelta(const TrackType& trk, const o2::track::PID::ID id) const; + /// Gets relative dEdx resolution contribution due to relative pt resolution + float GetRelativeResolutiondEdx(const float p, const float mass, const float charge, const float resol) const; + + void PrintAll() const; + + private: + std::array mBetheBlochParams = {0.03209809958934784, 19.9768009185791, 2.5266601063857674e-16, 2.7212300300598145, 6.080920219421387}; + std::array mResolutionParamsDefault = {0.07, 0.0}; + std::vector mResolutionParams = {5.43799e-7, 0.053044, 0.667584, 0.0142667, 0.00235175, 1.22482, 2.3501e-7, 0.031585}; + float mMIP = 50.f; + float mChargeFactor = 2.299999952316284f; + float mMultNormalization = 11000.; + bool mUseDefaultResolutionParam = true; + + ClassDefNV(Response, 2); + +}; // class Response + +/// Get expected Signal of the measurement +template +inline float Response::GetExpectedSignal(const TrackType& track, const o2::track::PID::ID id) const +{ + const float bethe = mMIP * o2::tpc::BetheBlochAleph(track.tpcInnerParam() / o2::track::pid_constants::sMasses[id], mBetheBlochParams[0], mBetheBlochParams[1], mBetheBlochParams[2], mBetheBlochParams[3], mBetheBlochParams[4]) * std::pow((float)o2::track::pid_constants::sCharges[id], mChargeFactor); + return bethe >= 0.f ? bethe : 0.f; +} + +/// Gets the expected resolution of the measurement +template +inline float Response::GetExpectedSigma(const CollisionType& collision, const TrackType& track, const o2::track::PID::ID id) const +{ + float resolution = 0.; + if (mUseDefaultResolutionParam) { + const float reso = track.tpcSignal() * mResolutionParamsDefault[0] * ((float)track.tpcNClsFound() > 0 ? std::sqrt(1. + mResolutionParamsDefault[1] / (float)track.tpcNClsFound()) : 1.f); + reso >= 0.f ? resolution = reso : resolution = 0.f; + } else { + + const double ncl = 159. / track.tpcNClsFound(); // + const double p = track.tpcInnerParam(); + const double mass = o2::track::pid_constants::sMasses[id]; + const double bg = p / mass; + const double dEdx = o2::tpc::BetheBlochAleph((float)bg, mBetheBlochParams[0], mBetheBlochParams[1], mBetheBlochParams[2], mBetheBlochParams[3], mBetheBlochParams[4]) * std::pow((float)o2::track::pid_constants::sCharges[id], mChargeFactor); + const double relReso = GetRelativeResolutiondEdx(p, mass, o2::track::pid_constants::sCharges[id], mResolutionParams[3]); + + const std::vector values{1.f / dEdx, track.tgl(), std::sqrt(ncl), relReso, track.signed1Pt(), collision.multTPC() / mMultNormalization}; + + const float reso = sqrt(pow(mResolutionParams[0], 2) * values[0] + pow(mResolutionParams[1], 2) * (values[2] * mResolutionParams[5]) * pow(values[0] / sqrt(1 + pow(values[1], 2)), mResolutionParams[2]) + values[2] * pow(values[3], 2) + pow(mResolutionParams[4] * values[4], 2) + pow(values[5] * mResolutionParams[6], 2) + pow(values[5] * (values[0] / sqrt(1 + pow(values[1], 2))) * mResolutionParams[7], 2)) * dEdx * mMIP; + reso >= 0.f ? resolution = reso : resolution = 0.f; + } + return resolution; +} + +/// Gets the number of sigma between the actual signal and the expected signal +template +inline float Response::GetNumberOfSigma(const CollisionType& collision, const TrackType& trk, const o2::track::PID::ID id) const +{ + return ((trk.tpcSignal() - GetExpectedSignal(trk, id)) / GetExpectedSigma(collision, trk, id)); +} + +/// Gets the deviation between the actual signal and the expected signal +template +inline float Response::GetSignalDelta(const TrackType& trk, const o2::track::PID::ID id) const +{ + return (trk.tpcSignal() - GetExpectedSignal(trk, id)); +} + +//// Gets relative dEdx resolution contribution due relative pt resolution +inline float Response::GetRelativeResolutiondEdx(const float p, const float mass, const float charge, const float resol) const +{ + const float bg = p / mass; + const float dEdx = o2::tpc::BetheBlochAleph(bg, mBetheBlochParams[0], mBetheBlochParams[1], mBetheBlochParams[2], mBetheBlochParams[3], mBetheBlochParams[4]) * std::pow(charge, mChargeFactor); + const float deltaP = resol * std::sqrt(dEdx); + const float bgDelta = p * (1 + deltaP) / mass; + const float dEdx2 = o2::tpc::BetheBlochAleph(bgDelta, mBetheBlochParams[0], mBetheBlochParams[1], mBetheBlochParams[2], mBetheBlochParams[3], mBetheBlochParams[4]) * std::pow(charge, mChargeFactor); + const float deltaRel = std::abs(dEdx2 - dEdx) / dEdx; + return deltaRel; +} + +inline void Response::PrintAll() const +{ + LOGP(info, "==== TPC PID response parameters: ===="); + for (int i = 0; i < int(mBetheBlochParams.size()); i++) { + LOGP(info, "BB param [{}] = {}", i, mBetheBlochParams[i]); + } + LOGP(info, "use default resolution parametrization = {}", mUseDefaultResolutionParam); + if (mUseDefaultResolutionParam) { + LOGP(info, "Default Resolution parametrization: "); + for (int i = 0; i < int(mResolutionParamsDefault.size()); i++) { + LOGP(info, "Resolution param [{}] = {}", i, mResolutionParamsDefault[i]); + } + } else { + for (int i = 0; i < int(mResolutionParams.size()); i++) { + LOGP(info, "Resolution param [{}] = {}", i, mResolutionParams[i]); + } + } + LOGP(info, "mMIP = {}", mMIP); + LOGP(info, "mChargeFactor = {}", mChargeFactor); + LOGP(info, "mMultNormalization = {}", mMultNormalization); +} + +} // namespace o2::pid::tpc + +#endif // O2_FRAMEWORK_PIDRESPONSE_H_ \ No newline at end of file diff --git a/Common/Core/PID/TPCReso.h b/Common/Core/PID/TPCReso.h deleted file mode 100644 index 231f6ecdf4e..00000000000 --- a/Common/Core/PID/TPCReso.h +++ /dev/null @@ -1,54 +0,0 @@ -// 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 TPCReso.h -/// \author Nicolo' Jacazio -/// \since 07/08/2020 -/// \brief Implementation for the TPC PID response of the BB resolution -/// - -#ifndef O2_ANALYSIS_PID_TPCRESO_H_ -#define O2_ANALYSIS_PID_TPCRESO_H_ - -#include "Common/Core/PID/ParamBase.h" -#include "ReconstructionDataFormats/PID.h" - -namespace o2::pid::tpc -{ - -class TPCReso : public Parametrization -{ - public: - TPCReso() : Parametrization("TPCReso", 2){}; - ~TPCReso() override = default; - float operator()(const float* x) const override - { - // relative dEdx resolution rel sigma = fRes0*sqrt(1+fResN2/npoint) - return x[0] * mParameters[0] * (x[1] > 0 ? sqrt(1. + mParameters[1] / x[1]) : 1.f); - } - ClassDef(TPCReso, 1); -}; - -float TPCResoParam(const float& signal, const float& npoints, const Parameters& parameters) -{ - return signal * parameters[0] * (npoints > 0 ? sqrt(1. + parameters[1] / npoints) : 1.f); -} - -template -float TPCResoParamTrack(const T& track, const Parameters& parameters) -{ - return TPCResoParam(track.tpcSignal(), (float)track.tpcNClsFound(), parameters); -} - -} // namespace o2::pid::tpc - -#endif diff --git a/Common/TableProducer/PID/CMakeLists.txt b/Common/TableProducer/PID/CMakeLists.txt index f2c542e2796..7ca41886b92 100644 --- a/Common/TableProducer/PID/CMakeLists.txt +++ b/Common/TableProducer/PID/CMakeLists.txt @@ -37,7 +37,7 @@ o2physics_add_dpl_workflow(pid-tpc SOURCES pidTPC.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore COMPONENT_NAME Analysis) - + o2physics_add_dpl_workflow(pid-tpc-full SOURCES pidTPCFull.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore @@ -60,4 +60,4 @@ o2physics_add_dpl_workflow(pid-hmpid-qa o2physics_add_dpl_workflow(pid-bayes SOURCES pidBayes.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore - COMPONENT_NAME Analysis) + COMPONENT_NAME Analysis) \ No newline at end of file diff --git a/Common/TableProducer/PID/pidBayes.cxx b/Common/TableProducer/PID/pidBayes.cxx index 6c0ad4a5859..ecc663bbb7e 100644 --- a/Common/TableProducer/PID/pidBayes.cxx +++ b/Common/TableProducer/PID/pidBayes.cxx @@ -24,7 +24,8 @@ #include #include "Common/Core/PID/PIDResponse.h" #include "Common/Core/PID/PIDTOF.h" -#include "Common/Core/PID/PIDTPC.h" +#include "Common/Core/PID/TPCPIDResponse.h" +#include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/TrackSelectionTables.h" using namespace o2; @@ -43,7 +44,7 @@ void customize(std::vector& workflowOptions) struct bayesPid { using Trks = soa::Join; - using Coll = aod::Collisions; + using Coll = soa::Join; // Tables to produce Produces tablePIDEl; /// Table for the Electron @@ -69,11 +70,13 @@ struct bayesPid { // Detector response and input parameters std::array Response; static constexpr const char* detectorName[kNDet] = {"TOF", "TPC"}; + // TPC PID Response + o2::pid::tpc::Response responseTPC; + o2::pid::tpc::Response* responseTPCptr = nullptr; Service ccdb; - Configurable paramfile{"param-file", "", "Path to the parametrization object, if emtpy the parametrization is not taken from file"}; + Configurable paramfileTOF{"param-file-TOF", "", "Path to the TOF parametrization object, if emtpy the parametrization is not taken from file"}; + Configurable paramfileTPC{"param-file-TPC", "", "Path to the TPC parametrization object, if emtpy the parametrization is not taken from file"}; Configurable TOFsigmaname{"param-tof-sigma", "TOFReso", "Name of the parametrization for the expected sigma, used in both file and CCDB mode"}; - Configurable TPCsignalname{"param-tpc-signal", "BetheBloch", "Name of the parametrization for the expected signal, used in both file and CCDB mode"}; - Configurable TPCsigmaname{"param-tpc-sigma", "TPCReso", "Name of the parametrization for the expected sigma, used in both file and CCDB mode"}; Configurable enableTOF{"enableTOF", false, "Enabling TOF"}; Configurable enableTPC{"enableTPC", false, "Enabling TPC"}; /// Ordering has to respect the one in ProbType @@ -82,7 +85,7 @@ struct bayesPid { Configurable url{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable ccdbPathTOF{"ccdbPathTOF", "Analysis/PID/TOF", "Path of the TOF parametrization on the CCDB"}; - Configurable ccdbPathTPC{"ccdbPathTPC", "Analysis/PID/TPC", "Path of the TPC parametrization on the CCDB"}; + Configurable ccdbPathTPC{"ccdbPathTPC", "Analysis/PID/TPC/Response", "Path of the TPC parametrization on the CCDB"}; Configurable timestamp{"ccdb-timestamp", -1, "timestamp of the object"}; // Configuration flags to include and exclude particle hypotheses // Configurable> pid{"pid", @@ -219,58 +222,54 @@ struct bayesPid { // const std::vector p = {0.008, 0.008, 0.002, 40.0}; Response[kTOF].SetParameters(DetectorResponse::kSigma, p); - const std::string fname = paramfile.value; - if (!fname.empty()) { // Loading the parametrization from file - LOG(info) << "Loading exp. sigma parametrization from file" << fname << ", using param: " << TOFsigmaname.value; - Response[kTOF].LoadParamFromFile(fname.data(), TOFsigmaname.value, DetectorResponse::kSigma); - - LOG(info) << "Loading exp. signal parametrization from file" << fname << ", using param: " << TPCsignalname.value; - Response[kTPC].LoadParamFromFile(fname.data(), TPCsignalname.value, DetectorResponse::kSignal); - - LOG(info) << "Loading exp. sigma parametrization from file" << fname << ", using param: " << TPCsigmaname.value; - Response[kTPC].LoadParamFromFile(fname.data(), TPCsigmaname.value, DetectorResponse::kSigma); - + const std::string fnameTOF = paramfileTOF.value; + const TString fnameTPC = paramfileTPC.value; + if (!fnameTOF.empty()) { // Loading the parametrization from file + LOG(info) << "Loading exp. sigma parametrization from file" << fnameTOF << ", using param: " << TOFsigmaname.value; + Response[kTOF].LoadParamFromFile(fnameTOF.data(), TOFsigmaname.value, DetectorResponse::kSigma); } else { // Loading it from CCDB std::string path = ccdbPathTOF.value + "/" + TOFsigmaname.value; LOG(info) << "Loading exp. sigma parametrization from CCDB, using path: " << path << " for timestamp " << timestamp.value; Response[kTOF].LoadParam(DetectorResponse::kSigma, ccdb->getForTimeStamp(path, timestamp.value)); - - path = ccdbPathTPC.value + "/" + TPCsignalname.value; - LOG(info) << "Loading exp. signal parametrization from CCDB, using path: " << path << " for timestamp " << timestamp.value; - Response[kTPC].LoadParam(DetectorResponse::kSignal, ccdb->getForTimeStamp(path, timestamp.value)); - - path = ccdbPathTPC.value + "/" + TPCsigmaname.value; - LOG(info) << "Loading exp. sigma parametrization from CCDB, using path: " << path << " for timestamp " << timestamp.value; - Response[kTPC].LoadParam(DetectorResponse::kSigma, ccdb->getForTimeStamp(path, timestamp.value)); + } + if (fnameTPC != "") { // Loading the parametrization from file + LOGP(info, "Loading TPC response from file {}", fnameTPC); + try { + std::unique_ptr f(TFile::Open(fnameTPC, "READ")); + f->GetObject("Response", responseTPCptr); + responseTPC.SetParameters(responseTPCptr); + } catch (...) { + LOGP(info, "Loading the TPC PID Response from file {} failed!", fnameTPC); + }; + } else { + const std::string pathTPC = ccdbPathTPC.value; + const auto time = timestamp.value; + responseTPC.SetParameters(ccdb->getForTimeStamp(pathTPC, time)); + LOGP(info, "Loading TPC response from CCDB, using path: {} for timestamp {}", pathTPC, time); + responseTPC.PrintAll(); } } - template - using respTPC = tpc::ELoss; - /// Computes PID probabilities for the TPC template - void ComputeTPCProbability(const Trks::iterator& track) + void ComputeTPCProbability(const Coll::iterator& collision, const Trks::iterator& track) { if (!checkEnabled()) { return; } - constexpr respTPC responseTPCPID; - const float dedx = track.tpcSignal(); bool mismatch = true; // if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)){ // dedx = GetTPCsignalTunedOnData(track); // } - - const float bethe = responseTPCPID.GetExpectedSignal(Response[kTPC], track); - const float sigma = responseTPCPID.GetExpectedSigma(Response[kTPC], track); - LOG(info) << "For " << pid_constants::sNames[pid] << " computing bethe " << bethe << " and sigma " << sigma; - // bethe = fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection, fUseTPCPileupCorrection); - // sigma = fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection, fUseTPCPileupCorrection); + const float bethe = responseTPC.GetExpectedSignal(track, pid); + const float sigma = responseTPC.GetExpectedSigma(collision, track, pid); + // LOG(info) << "For " << pid_constants::sNames[pid] << " computing bethe " << bethe << " and sigma " << sigma; + // bethe = fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection, fUseTPCPileupCorrection); + // sigma = fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection, fUseTPCPileupCorrection); if (abs(dedx - bethe) > fRange * sigma) { // Probability[kTPC][pid] = exp(-0.5 * fRange * fRange) / sigma; // BUG fix @@ -453,15 +452,16 @@ struct bayesPid { for (auto const& trk : tracks) { // Loop on Tracks - ComputeTPCProbability(trk); - ComputeTPCProbability(trk); - ComputeTPCProbability(trk); - ComputeTPCProbability(trk); - ComputeTPCProbability(trk); - ComputeTPCProbability(trk); - ComputeTPCProbability(trk); - ComputeTPCProbability(trk); - ComputeTPCProbability(trk); + auto collision = collisions.iteratorAt(trk.collisionId()); + ComputeTPCProbability(collision, trk); + ComputeTPCProbability(collision, trk); + ComputeTPCProbability(collision, trk); + ComputeTPCProbability(collision, trk); + ComputeTPCProbability(collision, trk); + ComputeTPCProbability(collision, trk); + ComputeTPCProbability(collision, trk); + ComputeTPCProbability(collision, trk); + ComputeTPCProbability(collision, trk); ComputeTOFProbability(trk); ComputeTOFProbability(trk); @@ -646,4 +646,4 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) workflow.push_back(adaptAnalysisTask(cfgc)); } return workflow; -} +} \ No newline at end of file diff --git a/Common/TableProducer/PID/pidTPC.cxx b/Common/TableProducer/PID/pidTPC.cxx index 06b30eea891..cf430a07157 100644 --- a/Common/TableProducer/PID/pidTPC.cxx +++ b/Common/TableProducer/PID/pidTPC.cxx @@ -11,22 +11,27 @@ /// /// \file pidTPC.cxx -/// \author Nicolò Jacazio +/// \author Annalena Kalteyer annalena.sophie.kalteyer@cern.ch +/// \author Christian Sonnabend christian.sonnabend@cern.ch +/// \author Nicolò Jacazio nicolo.jacazio@cern.ch /// \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. /// +#include "TFile.h" + // O2 includes #include "Framework/AnalysisTask.h" #include "Framework/HistogramRegistry.h" -#include "Framework/RunningWorkflowInfo.h" #include "ReconstructionDataFormats/Track.h" #include #include "Common/Core/PID/PIDResponse.h" -#include "Common/Core/PID/PIDTPC.h" +#include "Common/Core/PID/TPCPIDResponse.h" #include "Common/DataModel/TrackSelectionTables.h" -#include "Common/DataModel/EventSelection.h" +#include "Framework/AnalysisDataModel.h" +#include "Common/DataModel/Multiplicity.h" #include "TableHelper.h" +#include "Common/DataModel/EventSelection.h" #include "Framework/StaticFor.h" using namespace o2; @@ -46,7 +51,7 @@ void customize(std::vector& workflowOptions) /// Task to produce the TPC response table struct tpcPid { using Trks = soa::Join; - using Coll = aod::Collisions; + using Coll = soa::Join; // Tables to produce Produces tablePIDEl; Produces tablePIDMu; @@ -57,14 +62,14 @@ struct tpcPid { Produces tablePIDTr; Produces tablePIDHe; Produces tablePIDAl; - // Detector response and input parameters - DetectorResponse response; + // TPC PID Response + o2::pid::tpc::Response response; + o2::pid::tpc::Response* responseptr = nullptr; + // Input parameters Service ccdb; Configurable paramfile{"param-file", "", "Path to the parametrization object, if emtpy the parametrization is not taken from file"}; - Configurable signalname{"param-signal", "BetheBloch", "Name of the parametrization for the expected signal, used in both file and CCDB mode"}; - Configurable sigmaname{"param-sigma", "TPCReso", "Name of the parametrization for the expected sigma, used in both file and CCDB mode"}; Configurable url{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - Configurable ccdbPath{"ccdbPath", "Analysis/PID/TPC", "Path of the TPC parametrization on the CCDB"}; + Configurable ccdbPath{"ccdbPath", "Analysis/PID/TPC/Response", "Path of the TPC parametrization on the CCDB"}; Configurable ccdbTimestamp{"ccdb-timestamp", 0, "timestamp of the object used to query in CCDB the detector response. Exceptions: -1 gets the latest object, 0 gets the run dependent timestamp"}; // Configuration flags to include and exclude particle hypotheses Configurable pidEl{"pid-el", -1, {"Produce PID information for the Electron mass hypothesis, overrides the automatic setup: the corresponding table can be set off (0) or on (1)"}}; @@ -79,8 +84,6 @@ struct tpcPid { // Paramatrization configuration bool useCCDBParam = false; - std::string ccdbPathSignal = ""; - std::string ccdbPathSigma = ""; void init(o2::framework::InitContext& initContext) { @@ -99,7 +102,6 @@ struct tpcPid { } } }; - enableFlag("El", pidEl); enableFlag("Mu", pidMu); enableFlag("Pi", pidPi); @@ -110,102 +112,66 @@ struct tpcPid { enableFlag("He", pidHe); enableFlag("Al", pidAl); - // Getting the parametrization parameters - ccdb->setURL(url.value); - ccdb->setTimestamp(ccdbTimestamp.value); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - // Not later than now objects - ccdb->setCreatedNotAfter(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count()); - // - const std::string fname = paramfile.value; - if (!fname.empty()) { // Loading the parametrization from file - LOG(info) << "Loading exp. signal parametrization from file" << fname << ", using param: " << signalname.value; - response.LoadParamFromFile(fname.data(), signalname.value, DetectorResponse::kSignal); - - LOG(info) << "Loading exp. sigma parametrization from file" << fname << ", using param: " << sigmaname.value; - response.LoadParamFromFile(fname.data(), sigmaname.value, DetectorResponse::kSigma); - } else { // Loading it from CCDB + const TString fname = paramfile.value; + if (fname != "") { // Loading the parametrization from file + LOGP(info, "Loading TPC response from file {}", fname); + try { + std::unique_ptr f(TFile::Open(fname, "READ")); + f->GetObject("Response", responseptr); + response.SetParameters(responseptr); + } catch (...) { + LOGP(info, "Loading the TPC PID Response from file {} failed!", fname); + }; + } else { useCCDBParam = true; - ccdbPathSignal = ccdbPath.value + "/" + signalname.value; - LOG(info) << "Loading exp. signal parametrization from CCDB, using path: " << ccdbPathSignal << " for timestamp " << ccdbTimestamp.value; - response.LoadParam(DetectorResponse::kSignal, ccdb->getForTimeStamp(ccdbPathSignal, ccdbTimestamp.value)); - - ccdbPathSigma = ccdbPath.value + "/" + sigmaname.value; - LOG(info) << "Loading exp. sigma parametrization from CCDB, using path: " << ccdbPathSigma << " for timestamp " << ccdbTimestamp.value; - response.LoadParam(DetectorResponse::kSigma, ccdb->getForTimeStamp(ccdbPathSigma, ccdbTimestamp.value)); + const std::string path = ccdbPath.value; + const auto time = ccdbTimestamp.value; + ccdb->setURL(url.value); + ccdb->setTimestamp(time); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setCreatedNotAfter(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count()); + response.SetParameters(ccdb->getForTimeStamp(path, time)); + LOGP(info, "Loading TPC response from CCDB, using path: {} for ccdbTimestamp {}", path, time); + response.PrintAll(); } } - template - using ResponseImplementation = o2::pid::tpc::ELoss; - void process(Trks const& tracks, - aod::BCsWithTimestamps const&, - Coll const&) + void process(Coll const& collisions, Trks const& tracks, + aod::BCsWithTimestamps const&) { - constexpr auto responseEl = ResponseImplementation(); - constexpr auto responseMu = ResponseImplementation(); - constexpr auto responsePi = ResponseImplementation(); - constexpr auto responseKa = ResponseImplementation(); - constexpr auto responsePr = ResponseImplementation(); - constexpr auto responseDe = ResponseImplementation(); - constexpr auto responseTr = ResponseImplementation(); - constexpr auto responseHe = ResponseImplementation(); - constexpr auto responseAl = ResponseImplementation(); - - auto reserveTable = [&tracks](const Configurable& flag, auto& table) { - if (flag.value != 1) { - return; + // Check and fill enabled tables + auto makeTable = [&tracks, &collisions, this](const Configurable& flag, auto& table, const o2::track::PID::ID pid) { + if (flag.value == 1) { + // Prepare memory for enabled tables + table.reserve(tracks.size()); + int lastCollisionId = -1; // Last collision ID analysed + for (auto const& trk : tracks) { // Loop on Tracks + auto collision = collisions.iteratorAt(trk.collisionId()); + if (useCCDBParam && ccdbTimestamp.value == 0 && trk.has_collision() && trk.collisionId() != lastCollisionId) { // Updating parametrization only if the initial timestamp is 0 + lastCollisionId = trk.collisionId(); + const auto& bc = collision.bc_as(); + response.SetParameters(ccdb->getForTimeStamp(ccdbPath.value, bc.timestamp())); + } + const float numbersigma = response.GetNumberOfSigma(collision, trk, pid); + aod::pidutils::packInTable(numbersigma, table, + aod::pidtpc_tiny::binned_min, + aod::pidtpc_tiny::binned_max, + aod::pidtpc_tiny::bin_width); + } } - table.reserve(tracks.size()); }; - - reserveTable(pidEl, tablePIDEl); - reserveTable(pidMu, tablePIDMu); - reserveTable(pidPi, tablePIDPi); - reserveTable(pidKa, tablePIDKa); - reserveTable(pidPr, tablePIDPr); - reserveTable(pidDe, tablePIDDe); - reserveTable(pidTr, tablePIDTr); - reserveTable(pidHe, tablePIDHe); - reserveTable(pidAl, tablePIDAl); - - int lastCollisionId = -1; // Last collision ID analysed - - for (auto const& t : tracks) { // Loop on tracks - - if (useCCDBParam && ccdbTimestamp.value == 0 && t.has_collision() && t.collisionId() != lastCollisionId) { // Updating parametrization only if the initial timestamp is 0 - - lastCollisionId = t.collisionId(); /// Cache last collision ID - - const auto& bc = t.collision().bc_as(); - response.LoadParam(DetectorResponse::kSignal, ccdb->getForTimeStamp(ccdbPathSignal, bc.timestamp())); - response.LoadParam(DetectorResponse::kSigma, ccdb->getForTimeStamp(ccdbPathSigma, bc.timestamp())); - } - // Check and fill enabled tables - auto makeTable = [&t](const Configurable& flag, auto& table, const DetectorResponse& response, const auto& responsePID) { - if (flag.value != 1) { - return; - } - const float& separation = responsePID.GetSeparation(response, t); - aod::pidutils::packInTable(separation, table, - aod::pidtpc_tiny::binned_min, - aod::pidtpc_tiny::binned_max, - aod::pidtpc_tiny::bin_width); - }; - - makeTable(pidEl, tablePIDEl, response, responseEl); - makeTable(pidMu, tablePIDMu, response, responseMu); - makeTable(pidPi, tablePIDPi, response, responsePi); - makeTable(pidKa, tablePIDKa, response, responseKa); - makeTable(pidPr, tablePIDPr, response, responsePr); - makeTable(pidDe, tablePIDDe, response, responseDe); - makeTable(pidTr, tablePIDTr, response, responseTr); - makeTable(pidHe, tablePIDHe, response, responseHe); - makeTable(pidAl, tablePIDAl, response, responseAl); - } + makeTable(pidEl, tablePIDEl, o2::track::PID::Electron); + makeTable(pidMu, tablePIDMu, o2::track::PID::Muon); + makeTable(pidPi, tablePIDPi, o2::track::PID::Pion); + makeTable(pidKa, tablePIDKa, o2::track::PID::Kaon); + makeTable(pidPr, tablePIDPr, o2::track::PID::Proton); + makeTable(pidDe, tablePIDDe, o2::track::PID::Deuteron); + makeTable(pidTr, tablePIDTr, o2::track::PID::Triton); + makeTable(pidHe, tablePIDHe, o2::track::PID::Helium3); + makeTable(pidAl, tablePIDAl, o2::track::PID::Alpha); } }; @@ -228,14 +194,14 @@ struct tpcPidQa { HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::QAObject}; - Configurable logAxis{"logAxis", 0, "Flag to use a log momentum axis"}; - Configurable nBinsP{"nBinsP", 400, "Number of bins for the momentum"}; - Configurable minP{"minP", 0, "Minimum momentum in range"}; + Configurable logAxis{"logAxis", 1, "Flag to use a log momentum axis"}; + Configurable nBinsP{"nBinsP", 3000, "Number of bins for the momentum"}; + Configurable minP{"minP", 0.01, "Minimum momentum in range"}; Configurable maxP{"maxP", 20, "Maximum momentum 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 applyEvSel{"applyEvSel", 2, "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"}; @@ -243,20 +209,24 @@ struct tpcPidQa { void addParticleHistos(const AxisSpec& pAxis, const AxisSpec& ptAxis) { // 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}); + const char* axisTitle = Form("N_{#sigma}^{TPC}(%s)", pT[i]); + const AxisSpec nSigmaAxis{nBinsNSigma, minNSigma, maxNSigma, axisTitle}; + histos.add(hnsigma[i].data(), axisTitle, kTH2F, {pAxis, nSigmaAxis}); + histos.add(hnsigmapt[i].data(), axisTitle, kTH2F, {ptAxis, nSigmaAxis}); + histos.add(hnsigmapospt[i].data(), axisTitle, kTH2F, {ptAxis, nSigmaAxis}); + histos.add(hnsigmanegpt[i].data(), axisTitle, 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 etaAxis{100, -2, 2, "#it{#eta}"}; + const AxisSpec lAxis{100, 0, 500, "Track length (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})"}; + AxisSpec pAxis{nBinsP, minP, maxP, "#it{p} (GeV/#it{c})"}; if (logAxis) { ptAxis.makeLogaritmic(); pAxis.makeLogaritmic(); @@ -271,9 +241,17 @@ struct tpcPidQa { h->GetXaxis()->SetBinLabel(4, "Passed vtx Z"); histos.add("event/vertexz", "", kTH1F, {vtxZAxis}); - histos.add("event/multiplicity", "", kTH1F, {multAxis}); + h = histos.add("event/particlehypo", "", kTH1F, {{10, 0, 10, "PID in tracking"}}); + for (int i = 0; i < 9; i++) { + h->GetXaxis()->SetBinLabel(i + 1, PID::getName(i)); + } + histos.add("event/trackmultiplicity", "", kTH1F, {multAxis}); histos.add("event/tpcsignal", "", kTH2F, {pAxis, dedxAxis}); histos.add("event/signedtpcsignal", "", kTH2F, {pAxisPosNeg, dedxAxis}); + histos.add("event/eta", "", kTH1F, {etaAxis}); + histos.add("event/length", "", kTH1F, {lAxis}); + histos.add("event/pt", "", kTH1F, {ptAxis}); + histos.add("event/p", "", kTH1F, {pAxis}); static_for<0, 8>([&](auto i) { addParticleHistos(pAxis, ptAxis); @@ -289,8 +267,9 @@ struct tpcPidQa { return; } } + // Fill histograms - const auto& nsigma = o2::aod::pidutils::tpcNSigma(id, t); + const auto& nsigma = o2::aod::pidutils::tpcNSigma(t); histos.fill(HIST(hnsigma[id]), t.p(), nsigma); histos.fill(HIST(hnsigmapt[id]), t.pt(), nsigma); if (t.sign() > 0) { @@ -307,6 +286,7 @@ struct tpcPidQa { aod::pidTPCTr, aod::pidTPCHe, aod::pidTPCAl, aod::TrackSelection> const& tracks) { + histos.fill(HIST("event/evsel"), 1); if (applyEvSel == 1) { if (!collision.sel7()) { @@ -317,8 +297,10 @@ struct tpcPidQa { return; } } + histos.fill(HIST("event/evsel"), 2); + // Computing Multiplicity first float ntracks = 0; for (auto t : tracks) { if (applyTrackCut && !t.isGlobalTrack()) { @@ -335,16 +317,19 @@ struct tpcPidQa { } histos.fill(HIST("event/evsel"), 4); histos.fill(HIST("event/vertexz"), collision.posZ()); - histos.fill(HIST("event/multiplicity"), ntracks); + histos.fill(HIST("event/trackmultiplicity"), 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()); + histos.fill(HIST("event/particlehypo"), t.pidForTracking()); + histos.fill(HIST("event/tpcsignal"), t.tpcInnerParam(), t.tpcSignal()); + histos.fill(HIST("event/signedtpcsignal"), t.tpcInnerParam() * t.sign(), t.tpcSignal()); + histos.fill(HIST("event/eta"), t.eta()); + histos.fill(HIST("event/length"), t.length()); + histos.fill(HIST("event/pt"), t.pt()); + histos.fill(HIST("event/p"), t.p()); // fillParticleHistos(t); fillParticleHistos(t); @@ -366,4 +351,4 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) workflow.push_back(adaptAnalysisTask(cfgc)); } return workflow; -} +} \ No newline at end of file diff --git a/Common/TableProducer/PID/pidTPCFull.cxx b/Common/TableProducer/PID/pidTPCFull.cxx index 5702bb839b7..127cf05b8ce 100644 --- a/Common/TableProducer/PID/pidTPCFull.cxx +++ b/Common/TableProducer/PID/pidTPCFull.cxx @@ -11,22 +11,27 @@ /// /// \file pidTPCFull.cxx +/// \author Annalena Kalteyer annalena.sophie.kalteyer@cern.ch +/// \author Christian Sonnabend christian.sonnabend@cern.ch /// \author Nicolò Jacazio nicolo.jacazio@cern.ch /// \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. /// +#include "TFile.h" + // O2 includes #include "Framework/AnalysisTask.h" #include "Framework/HistogramRegistry.h" -#include "Framework/RunningWorkflowInfo.h" #include "ReconstructionDataFormats/Track.h" #include #include "Common/Core/PID/PIDResponse.h" -#include "Common/Core/PID/PIDTPC.h" +#include "Common/Core/PID/TPCPIDResponse.h" #include "Common/DataModel/TrackSelectionTables.h" -#include "Common/DataModel/EventSelection.h" +#include "Framework/AnalysisDataModel.h" +#include "Common/DataModel/Multiplicity.h" #include "TableHelper.h" +#include "Common/DataModel/EventSelection.h" #include "Framework/StaticFor.h" using namespace o2; @@ -46,7 +51,8 @@ void customize(std::vector& workflowOptions) /// Task to produce the TPC response table struct tpcPidFull { using Trks = soa::Join; - using Coll = aod::Collisions; + using Coll = soa::Join; + // Tables to produce Produces tablePIDEl; Produces tablePIDMu; @@ -57,15 +63,15 @@ struct tpcPidFull { Produces tablePIDTr; Produces tablePIDHe; Produces tablePIDAl; - // Detector response and input parameters - DetectorResponse response; + // TPC PID Response + o2::pid::tpc::Response response; + o2::pid::tpc::Response* responseptr = nullptr; + // Input parameters Service ccdb; Configurable paramfile{"param-file", "", "Path to the parametrization object, if emtpy the parametrization is not taken from file"}; - Configurable signalname{"param-signal", "BetheBloch", "Name of the parametrization for the expected signal, used in both file and CCDB mode"}; - Configurable sigmaname{"param-sigma", "TPCReso", "Name of the parametrization for the expected sigma, used in both file and CCDB mode"}; Configurable url{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - Configurable ccdbPath{"ccdbPath", "Analysis/PID/TPC", "Path of the TPC parametrization on the CCDB"}; - Configurable timestamp{"ccdb-timestamp", -1, "timestamp of the object"}; + Configurable ccdbPath{"ccdbPath", "Analysis/PID/TPC/Response", "Path of the TPC parametrization on the CCDB"}; + Configurable ccdbTimestamp{"ccdb-timestamp", 0, "timestamp of the object used to query in CCDB the detector response. Exceptions: -1 gets the latest object, 0 gets the run dependent timestamp"}; // Configuration flags to include and exclude particle hypotheses Configurable pidEl{"pid-el", -1, {"Produce PID information for the Electron mass hypothesis, overrides the automatic setup: the corresponding table can be set off (0) or on (1)"}}; Configurable pidMu{"pid-mu", -1, {"Produce PID information for the Muon mass hypothesis, overrides the automatic setup: the corresponding table can be set off (0) or on (1)"}}; @@ -77,6 +83,9 @@ struct tpcPidFull { Configurable pidHe{"pid-he", -1, {"Produce PID information for the Helium3 mass hypothesis, overrides the automatic setup: the corresponding table can be set off (0) or on (1)"}}; Configurable pidAl{"pid-al", -1, {"Produce PID information for the Alpha mass hypothesis, overrides the automatic setup: the corresponding table can be set off (0) or on (1)"}}; + // Paramatrization configuration + bool useCCDBParam = false; + void init(o2::framework::InitContext& initContext) { // Checking the tables are requested in the workflow and enabling them @@ -94,7 +103,6 @@ struct tpcPidFull { } } }; - enableFlag("El", pidEl); enableFlag("Mu", pidMu); enableFlag("Pi", pidPi); @@ -105,66 +113,62 @@ struct tpcPidFull { enableFlag("He", pidHe); enableFlag("Al", pidAl); - // Getting the parametrization parameters - ccdb->setURL(url.value); - ccdb->setTimestamp(timestamp.value); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - // Not later than now objects - ccdb->setCreatedNotAfter(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count()); - // - const std::string fname = paramfile.value; - if (!fname.empty()) { // Loading the parametrization from file - LOG(info) << "Loading exp. signal parametrization from file" << fname << ", using param: " << signalname.value; - response.LoadParamFromFile(fname.data(), signalname.value, DetectorResponse::kSignal); - - LOG(info) << "Loading exp. sigma parametrization from file" << fname << ", using param: " << sigmaname.value; - response.LoadParamFromFile(fname.data(), sigmaname.value, DetectorResponse::kSigma); - } else { // Loading it from CCDB - std::string path = ccdbPath.value + "/" + signalname.value; - LOG(info) << "Loading exp. signal parametrization from CCDB, using path: " << path << " for timestamp " << timestamp.value; - response.LoadParam(DetectorResponse::kSignal, ccdb->getForTimeStamp(path, timestamp.value)); - - path = ccdbPath.value + "/" + sigmaname.value; - LOG(info) << "Loading exp. sigma parametrization from CCDB, using path: " << path << " for timestamp " << timestamp.value; - response.LoadParam(DetectorResponse::kSigma, ccdb->getForTimeStamp(path, timestamp.value)); + const TString fname = paramfile.value; + if (fname != "") { // Loading the parametrization from file + LOGP(info, "Loading TPC response from file {}", fname); + try { + std::unique_ptr f(TFile::Open(fname, "READ")); + f->GetObject("Response", responseptr); + response.SetParameters(responseptr); + } catch (...) { + LOGP(info, "Loading the TPC PID Response from file {} failed!", fname); + }; + } else { + useCCDBParam = true; + const std::string path = ccdbPath.value; + const auto time = ccdbTimestamp.value; + ccdb->setURL(url.value); + ccdb->setTimestamp(time); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setCreatedNotAfter(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count()); + response.SetParameters(ccdb->getForTimeStamp(path, time)); + LOGP(info, "Loading TPC response from CCDB, using path: {} for ccdbTimestamp {}", path, time); + response.PrintAll(); } } - template - using ResponseImplementation = o2::pid::tpc::ELoss; - void process(Coll const& collisions, Trks const& tracks) + void process(Coll const& collisions, Trks const& tracks, + aod::BCsWithTimestamps const&) { - constexpr auto responseEl = ResponseImplementation(); - constexpr auto responseMu = ResponseImplementation(); - constexpr auto responsePi = ResponseImplementation(); - constexpr auto responseKa = ResponseImplementation(); - constexpr auto responsePr = ResponseImplementation(); - constexpr auto responseDe = ResponseImplementation(); - constexpr auto responseTr = ResponseImplementation(); - constexpr auto responseHe = ResponseImplementation(); - constexpr auto responseAl = ResponseImplementation(); - // Check and fill enabled tables - auto makeTable = [&tracks](const Configurable& flag, auto& table, const DetectorResponse& response, const auto& responsePID) { - if (flag.value == 1) { + auto makeTable = [&tracks, &collisions, this](const Configurable& flag, auto& table, const o2::track::PID::ID pid) { + if (flag.value == 1) { // Loop on Tracks // Prepare memory for enabled tables table.reserve(tracks.size()); - for (auto const& trk : tracks) { // Loop on Tracks - table(responsePID.GetExpectedSigma(response, trk), - responsePID.GetSeparation(response, trk)); + int lastCollisionId = -1; // Last collision ID analysed + for (auto const& trk : tracks) { + auto collision = collisions.iteratorAt(trk.collisionId()); // Loop on Tracks + if (useCCDBParam && ccdbTimestamp.value == 0 && trk.has_collision() && trk.collisionId() != lastCollisionId) { // Updating parametrization only if the initial timestamp is 0 + lastCollisionId = trk.collisionId(); + const auto& bc = collision.bc_as(); + response.SetParameters(ccdb->getForTimeStamp(ccdbPath.value, bc.timestamp())); + } + table(response.GetExpectedSigma(collision, trk, pid), + response.GetNumberOfSigma(collision, trk, pid)); } } }; - makeTable(pidEl, tablePIDEl, response, responseEl); - makeTable(pidMu, tablePIDMu, response, responseMu); - makeTable(pidPi, tablePIDPi, response, responsePi); - makeTable(pidKa, tablePIDKa, response, responseKa); - makeTable(pidPr, tablePIDPr, response, responsePr); - makeTable(pidDe, tablePIDDe, response, responseDe); - makeTable(pidTr, tablePIDTr, response, responseTr); - makeTable(pidHe, tablePIDHe, response, responseHe); - makeTable(pidAl, tablePIDAl, response, responseAl); + // const o2::pid::tpc::Response& response; + makeTable(pidEl, tablePIDEl, o2::track::PID::Electron); + makeTable(pidMu, tablePIDMu, o2::track::PID::Muon); + makeTable(pidPi, tablePIDPi, o2::track::PID::Pion); + makeTable(pidKa, tablePIDKa, o2::track::PID::Kaon); + makeTable(pidPr, tablePIDPr, o2::track::PID::Proton); + makeTable(pidDe, tablePIDDe, o2::track::PID::Deuteron); + makeTable(pidTr, tablePIDTr, o2::track::PID::Triton); + makeTable(pidHe, tablePIDHe, o2::track::PID::Helium3); + makeTable(pidAl, tablePIDAl, o2::track::PID::Alpha); } }; @@ -196,9 +200,9 @@ struct tpcPidFullQa { HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::QAObject}; - Configurable logAxis{"logAxis", 0, "Flag to use a log momentum axis"}; - Configurable nBinsP{"nBinsP", 400, "Number of bins for the momentum"}; - Configurable minP{"minP", 0, "Minimum momentum in range"}; + Configurable logAxis{"logAxis", 1, "Flag to use a log momentum axis"}; + Configurable nBinsP{"nBinsP", 3000, "Number of bins for the momentum"}; + Configurable minP{"minP", 0.01, "Minimum momentum in range"}; Configurable maxP{"maxP", 20, "Maximum momentum in range"}; Configurable nBinsDelta{"nBinsDelta", 200, "Number of bins for the Delta"}; Configurable minDelta{"minDelta", -1000.f, "Minimum Delta in range"}; @@ -209,7 +213,7 @@ struct tpcPidFullQa { 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 applyEvSel{"applyEvSel", 2, "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"}; @@ -229,20 +233,24 @@ struct tpcPidFullQa { 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(), 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}); + const char* axisTitle = Form("N_{#sigma}^{TPC}(%s)", pT[i]); + const AxisSpec nSigmaAxis{nBinsNSigma, minNSigma, maxNSigma, axisTitle}; + histos.add(hnsigma[i].data(), axisTitle, kTH2F, {pAxis, nSigmaAxis}); + histos.add(hnsigmapt[i].data(), axisTitle, kTH2F, {ptAxis, nSigmaAxis}); + histos.add(hnsigmapospt[i].data(), axisTitle, kTH2F, {ptAxis, nSigmaAxis}); + histos.add(hnsigmanegpt[i].data(), axisTitle, 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 etaAxis{100, -2, 2, "#it{#eta}"}; + const AxisSpec lAxis{100, 0, 500, "Track length (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})"}; + AxisSpec pAxis{nBinsP, minP, maxP, "#it{p} (GeV/#it{c})"}; if (logAxis) { ptAxis.makeLogaritmic(); pAxis.makeLogaritmic(); @@ -257,9 +265,17 @@ struct tpcPidFullQa { h->GetXaxis()->SetBinLabel(4, "Passed vtx Z"); histos.add("event/vertexz", "", kTH1F, {vtxZAxis}); - histos.add("event/multiplicity", "", kTH1F, {multAxis}); + h = histos.add("event/particlehypo", "", kTH1F, {{10, 0, 10, "PID in tracking"}}); + for (int i = 0; i < 9; i++) { + h->GetXaxis()->SetBinLabel(i + 1, PID::getName(i)); + } + histos.add("event/trackmultiplicity", "", kTH1F, {multAxis}); histos.add("event/tpcsignal", "", kTH2F, {pAxis, dedxAxis}); histos.add("event/signedtpcsignal", "", kTH2F, {pAxisPosNeg, dedxAxis}); + histos.add("event/eta", "", kTH1F, {etaAxis}); + histos.add("event/length", "", kTH1F, {lAxis}); + histos.add("event/pt", "", kTH1F, {ptAxis}); + histos.add("event/p", "", kTH1F, {pAxis}); static_for<0, 8>([&](auto i) { addParticleHistos(pAxis, ptAxis); @@ -267,7 +283,7 @@ struct tpcPidFullQa { } template - void fillParticleHistos(const T& t, const float& mom, const float& exp_diff, const float& expsigma) + void fillParticleHistos(const T& t) { if (applyRapidityCut) { const float y = TMath::ASinH(t.pt() / TMath::Sqrt(PID::getMass2(id) + t.pt() * t.pt()) * TMath::SinH(t.eta())); @@ -275,11 +291,14 @@ struct tpcPidFullQa { return; } } + + const auto& nsigma = o2::aod::pidutils::tpcNSigma(t); + const auto& diff = o2::aod::pidutils::tpcExpSignalDiff(t); + // 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); - const auto& nsigma = o2::aod::pidutils::tpcNSigma(id, t); + histos.fill(HIST(hexpected[id]), t.tpcInnerParam(), t.tpcSignal() - diff); + histos.fill(HIST(hexpected_diff[id]), t.tpcInnerParam(), diff); + histos.fill(HIST(hexpsigma[id]), t.tpcInnerParam(), o2::aod::pidutils::tpcExpSigma(t)); histos.fill(HIST(hnsigma[id]), t.p(), nsigma); histos.fill(HIST(hnsigmapt[id]), t.pt(), nsigma); if (t.sign() > 0) { @@ -296,6 +315,7 @@ struct tpcPidFullQa { aod::pidTPCFullTr, aod::pidTPCFullHe, aod::pidTPCFullAl, aod::TrackSelection> const& tracks) { + histos.fill(HIST("event/evsel"), 1); if (applyEvSel == 1) { if (!collision.sel7()) { @@ -306,8 +326,10 @@ struct tpcPidFullQa { return; } } + histos.fill(HIST("event/evsel"), 2); + // Computing Multiplicity first float ntracks = 0; for (auto t : tracks) { if (applyTrackCut && !t.isGlobalTrack()) { @@ -324,7 +346,190 @@ struct tpcPidFullQa { } histos.fill(HIST("event/evsel"), 4); histos.fill(HIST("event/vertexz"), collision.posZ()); - histos.fill(HIST("event/multiplicity"), ntracks); + histos.fill(HIST("event/trackmultiplicity"), ntracks); + + for (auto t : tracks) { + if (applyTrackCut && !t.isGlobalTrack()) { + continue; + } + histos.fill(HIST("event/particlehypo"), t.pidForTracking()); + histos.fill(HIST("event/tpcsignal"), t.tpcInnerParam(), t.tpcSignal()); + histos.fill(HIST("event/signedtpcsignal"), t.tpcInnerParam() * t.sign(), t.tpcSignal()); + histos.fill(HIST("event/eta"), t.eta()); + histos.fill(HIST("event/length"), t.length()); + histos.fill(HIST("event/pt"), t.pt()); + histos.fill(HIST("event/p"), t.p()); + // + fillParticleHistos(t); + fillParticleHistos(t); + fillParticleHistos(t); + fillParticleHistos(t); + fillParticleHistos(t); + fillParticleHistos(t); + fillParticleHistos(t); + fillParticleHistos(t); + fillParticleHistos(t); + } + } +}; + +/// Task to produce the TPC QA plots with TOF PID +struct tpcPidFullQaWTof { + 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 hexpected[Np] = {"expected/El", "expected/Mu", "expected/Pi", + "expected/Ka", "expected/Pr", "expected/De", + "expected/Tr", "expected/He", "expected/Al"}; + 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 hsignal[Np] = {"signal/El", "signal/Mu", "signal/Pi", + "signal/Ka", "signal/Pr", "signal/De", + "signal/Tr", "signal/He", "signal/Al"}; + + HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::QAObject}; + + Configurable logAxis{"logAxis", 1, "Flag to use a log momentum axis"}; + Configurable nBinsP{"nBinsP", 3000, "Number of bins for the momentum"}; + Configurable minP{"minP", 0.01, "Minimum momentum in range"}; + Configurable maxP{"maxP", 20, "Maximum momentum in range"}; + 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", 2, "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(const AxisSpec& pAxis, const AxisSpec& ptAxis) + { + // Exp signal + const AxisSpec expAxis{1000, 0, 1000, Form("d#it{E}/d#it{x}_(%s) A.U.", pT[id])}; + histos.add(hexpected[id].data(), "", kTH2F, {pAxis, expAxis}); + + // Signal - Expected signal + const AxisSpec deltaAxis{nBinsDelta, minDelta, maxDelta, Form("d#it{E}/d#it{x} - d#it{E}/d#it{x}(%s)", pT[id])}; + histos.add(hexpected_diff[id].data(), "", kTH2F, {pAxis, deltaAxis}); + + // Exp Sigma + const AxisSpec expSigmaAxis{nBinsExpSigma, minExpSigma, maxExpSigma, Form("Exp_{#sigma}^{TPC}(%s)", pT[id])}; + histos.add(hexpsigma[id].data(), "", kTH2F, {pAxis, expSigmaAxis}); + + // NSigma + const char* axisTitle = Form("N_{#sigma}^{TPC}(%s)", pT[id]); + const AxisSpec nSigmaAxis{nBinsNSigma, minNSigma, maxNSigma, axisTitle}; + histos.add(hnsigma[id].data(), axisTitle, kTH2F, {pAxis, nSigmaAxis}); + histos.add(hnsigmapt[id].data(), axisTitle, kTH2F, {ptAxis, nSigmaAxis}); + const AxisSpec dedxAxis{1000, 0, 1000, "d#it{E}/d#it{x} A.U."}; + histos.add(hsignal[id].data(), "", kTH2F, {pAxis, dedxAxis}); + } + + 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(); + } + + // 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}); + h = histos.add("event/particlehypo", "", kTH1F, {{10, 0, 10, "PID in tracking"}}); + for (int i = 0; i < 9; i++) { + h->GetXaxis()->SetBinLabel(i + 1, PID::getName(i)); + } + histos.add("event/trackmultiplicity", "", kTH1F, {multAxis}); + + static_for<0, 8>([&](auto i) { + addParticleHistos(pAxis, ptAxis); + }); + } + + template + void fillParticleHistos(const T& t, const float& mom) + { + 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 + const auto& nsigma = o2::aod::pidutils::tpcNSigma(t); + const auto& nsigmatof = o2::aod::pidutils::tofNSigma(t); + const auto& diff = o2::aod::pidutils::tpcExpSignalDiff(t); + if (std::abs(nsigmatof) < 3.f) { + histos.fill(HIST(hexpected[id]), mom, t.tpcSignal() - diff); + histos.fill(HIST(hexpected_diff[id]), mom, diff); + histos.fill(HIST(hexpsigma[id]), t.p(), o2::aod::pidutils::tpcExpSigma(t)); + histos.fill(HIST(hnsigma[id]), t.p(), nsigma); + histos.fill(HIST(hnsigmapt[id]), t.pt(), nsigma); + histos.fill(HIST(hsignal[id]), mom, t.tpcSignal()); + // histos.fill(HIST("event/signedtpcsignal"), mom * t.sign(), t.tpcSignal()); + } + } + + 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; + } + 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/trackmultiplicity"), ntracks); for (auto t : tracks) { if (applyTrackCut && !t.isGlobalTrack()) { @@ -332,18 +537,17 @@ struct tpcPidFullQa { } // 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()); + histos.fill(HIST("event/particlehypo"), t.pidForTracking()); // - fillParticleHistos(t, mom, t.tpcExpSignalDiffEl(), t.tpcExpSigmaEl()); - fillParticleHistos(t, mom, t.tpcExpSignalDiffMu(), t.tpcExpSigmaMu()); - fillParticleHistos(t, mom, t.tpcExpSignalDiffPi(), t.tpcExpSigmaPi()); - fillParticleHistos(t, mom, t.tpcExpSignalDiffKa(), t.tpcExpSigmaKa()); - fillParticleHistos(t, mom, t.tpcExpSignalDiffPr(), t.tpcExpSigmaPr()); - fillParticleHistos(t, mom, t.tpcExpSignalDiffDe(), t.tpcExpSigmaDe()); - fillParticleHistos(t, mom, t.tpcExpSignalDiffTr(), t.tpcExpSigmaTr()); - fillParticleHistos(t, mom, t.tpcExpSignalDiffHe(), t.tpcExpSigmaHe()); - fillParticleHistos(t, mom, t.tpcExpSignalDiffAl(), t.tpcExpSigmaAl()); + fillParticleHistos(t, mom); + fillParticleHistos(t, mom); + fillParticleHistos(t, mom); + fillParticleHistos(t, mom); + fillParticleHistos(t, mom); + fillParticleHistos(t, mom); + fillParticleHistos(t, mom); + fillParticleHistos(t, mom); + fillParticleHistos(t, mom); } } }; @@ -353,6 +557,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) auto workflow = WorkflowSpec{adaptAnalysisTask(cfgc)}; if (cfgc.options().get("add-qa")) { workflow.push_back(adaptAnalysisTask(cfgc)); + workflow.push_back(adaptAnalysisTask(cfgc)); } return workflow; -} +} \ No newline at end of file diff --git a/Common/TableProducer/PID/qaTPCMC.cxx b/Common/TableProducer/PID/qaTPCMC.cxx index b8dc5320c01..ab32ccb3443 100644 --- a/Common/TableProducer/PID/qaTPCMC.cxx +++ b/Common/TableProducer/PID/qaTPCMC.cxx @@ -23,11 +23,9 @@ #include "ReconstructionDataFormats/Track.h" #include #include "Common/Core/PID/PIDResponse.h" -#include "Common/Core/PID/PIDTPC.h" using namespace o2; using namespace o2::framework; -using namespace o2::pid; using namespace o2::framework::expressions; using namespace o2::track; @@ -270,4 +268,4 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) workflow.push_back(adaptAnalysisTask>(cfgc, TaskName{"pidTPC-qa-Al"})); } return workflow; -} +} \ No newline at end of file diff --git a/Common/Tools/CMakeLists.txt b/Common/Tools/CMakeLists.txt index f664827c93e..4bcad206458 100644 --- a/Common/Tools/CMakeLists.txt +++ b/Common/Tools/CMakeLists.txt @@ -1,3 +1,4 @@ + # 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. @@ -14,13 +15,12 @@ o2physics_add_executable(aod-data-model-graph PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore ) -o2physics_add_executable(pidparam-tpc-bethe-bloch - SOURCES handleParamTPCBetheBloch.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore - ) - o2physics_add_executable(pidparam-tof-reso SOURCES handleParamTOFReso.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore ) +o2physics_add_executable(pidparam-tpc-response + SOURCES handleParamTPCResponse.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + ) \ No newline at end of file diff --git a/Common/Tools/handleParamTPCBetheBloch.cxx b/Common/Tools/handleParamTPCBetheBloch.cxx deleted file mode 100644 index 8e619f62f85..00000000000 --- a/Common/Tools/handleParamTPCBetheBloch.cxx +++ /dev/null @@ -1,166 +0,0 @@ -// 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 handleParamTPCBetheBloch.cxx -/// \author Nicolò Jacazio nicolo.jacazio@cern.ch -/// \since 2020-06-22 -/// \brief A simple tool to produce Bethe Bloch parametrization objects for the TPC PID Response -/// - -#include "Common/Core/PID/BetheBloch.h" -#include "Common/Core/PID/TPCReso.h" -#include "handleParamBase.h" - -using namespace o2::pid::tpc; - -bool initOptionsAndParse(bpo::options_description& options, int argc, char* argv[]) -{ - options.add_options()( - "url,u", bpo::value()->default_value("http://alice-ccdb.cern.ch"), "URL of the CCDB database e.g. http://ccdb-test.cern.ch:8080 or http://alice-ccdb.cern.ch")( - "ccdb-path,c", bpo::value()->default_value("Analysis/PID/TPC"), "CCDB path for storage/retrieval")( - "rct-path", bpo::value()->default_value("RCT/RunInformation"), "path to the ccdb RCT objects for the SOR/EOR timestamps")( - "start,s", bpo::value()->default_value(0), "Start timestamp of object validity. If 0 and runnumber != 0 it will be set to the run SOR")( - "stop,S", bpo::value()->default_value(0), "Stop timestamp of object validity. If 0 and runnumber != 0 it will be set to the run EOR")( - "timestamp,T", bpo::value()->default_value(-1), "Timestamp of the object to retrieve, used in alternative to the run number")( - "runnumber,R", bpo::value()->default_value(0), "Timestamp of the object to retrieve, used in alternative to the timestamp (if 0 using the timestamp)")( - "delete-previous,delete_previous,d", bpo::value()->default_value(0), "Flag to delete previous versions of converter objects in the CCDB before uploading the new one so as to avoid proliferation on CCDB")( - "save-to-file,file,f,o", bpo::value()->default_value(""), "Option to save parametrization to file instead of uploading to ccdb")( - "read-from-file,i", bpo::value()->default_value(""), "Option to get parametrization from a file")( - "exp-name,n", bpo::value()->default_value("BetheBloch"), "Name of the parametrization object")( - "reso-name,n", bpo::value()->default_value("TPCReso"), "Name of the parametrization object")( - "mode,m", bpo::value()->default_value(1), "Working mode: 0 push 1 pull and test")( - "p0", bpo::value()->default_value(0.0320981), "Parameter 0 of the TPC expected value")( - "p1", bpo::value()->default_value(19.9768), "Parameter 1 of the TPC expected value")( - "p2", bpo::value()->default_value(2.52666e-16), "Parameter 2 of the TPC expected value")( - "p3", bpo::value()->default_value(2.72123), "Parameter 3 of the TPC expected value")( - "p4", bpo::value()->default_value(6.08092), "Parameter 4 of the TPC expected value")( - "p5", bpo::value()->default_value(50.f), "Parameter 5 of the TPC expected value")( - "p6", bpo::value()->default_value(2.3), "Parameter 6 of the TPC expected value")( - "dryrun,D", bpo::value()->default_value(0), "Perform a dryrun check before uploading")( - "verbose,v", bpo::value()->default_value(0), "Verbose level 0, 1")( - "help,h", "Produce help message."); - try { - bpo::store(parse_command_line(argc, argv, options), arguments); - - // help - if (arguments.count("help")) { - LOG(info) << options; - return false; - } - - bpo::notify(arguments); - } catch (const bpo::error& e) { - LOG(error) << e.what() << "\n"; - LOG(error) << "Error parsing command line arguments; Available options:"; - LOG(error) << options; - return false; - } - return true; -} - -int main(int argc, char* argv[]) -{ - bpo::options_description options("Allowed options"); - if (!initOptionsAndParse(options, argc, argv)) { - return 1; - } - - // Fetch options - const auto mode = arguments["mode"].as(); - const auto runnumber = arguments["runnumber"].as(); - auto timestamp = arguments["timestamp"].as(); - const auto path = arguments["ccdb-path"].as(); - auto start = arguments["start"].as(); - auto stop = arguments["stop"].as(); - - // Init CCDB - const std::string url = arguments["url"].as(); - api.init(url); - if (!api.isHostReachable()) { - LOG(warning) << "CCDB host " << url << " is not reacheable, cannot go forward"; - return 1; - } - - // Init timestamps - setupTimestamps(timestamp, start, stop); - - BetheBloch* bb = nullptr; - TPCReso* reso = nullptr; - const std::string exp_name = arguments["exp-name"].as(); - const std::string reso_name = arguments["reso-name"].as(); - if (mode == 0) { // Push mode - LOG(info) << "Handling TPC parametrization in create mode"; - const std::string input_file_name = arguments["read-from-file"].as(); - if (!input_file_name.empty()) { - TFile f(input_file_name.data(), "READ"); - if (!f.IsOpen()) { - LOG(warning) << "Input file " << input_file_name << " is not reacheable, cannot get param from file"; - } - f.GetObject(exp_name.c_str(), bb); - f.GetObject(reso_name.c_str(), reso); - f.Close(); - } - if (!bb) { - bb = new BetheBloch(); - const std::vector bbparams = {arguments["p0"].as(), arguments["p1"].as(), arguments["p2"].as(), arguments["p3"].as(), arguments["p4"].as(), arguments["p5"].as(), arguments["p6"].as()}; - bb->SetParameters(bbparams); - } - bb->Print(); - if (!reso) { - reso = new TPCReso(); - const std::vector resoparams = {0.07, 0.0}; - reso->SetParameters(resoparams); - } - reso->Print(); - const std::string fname = arguments["save-to-file"].as(); - if (!fname.empty()) { // Saving it to file - LOG(info) << "Saving parametrization to file " << fname; - TFile f(fname.data(), "RECREATE"); - bb->Write(); - bb->GetParameters().Write(); - reso->Write(); - reso->GetParameters().Write(); - f.ls(); - f.Close(); - } else { // Saving it to CCDB - LOG(info) << "Saving parametrization to CCDB " << path << " with validity " << start << " " - << stop; - - std::map metadata; - if (runnumber != 0) { - metadata["runnumber"] = Form("%i", runnumber); - } - // Storing parametrization objects - storeOnCCDB(path + "/" + exp_name, metadata, start, stop, bb); - storeOnCCDB(path + "/" + reso_name, metadata, start, stop, reso); - // Storing parametrization parameters - o2::pid::Parameters* params; - bb->GetParameters(params); - storeOnCCDB(path + "/Parameters/" + exp_name, metadata, start, stop, params); - reso->GetParameters(params); - storeOnCCDB(path + "/Parameters/" + reso_name, metadata, start, stop, params); - } - } else { // Pull and test mode - LOG(info) << "Handling TPC parametrization in test mode for timestamp " - << timestamp << " -> " << timeStampToHReadble(timestamp); - const float x[2] = {1, 1}; - bb = retrieveFromCCDB(path + "/" + exp_name, timestamp); - bb->Print(); - LOG(info) << "BetheBloch " << bb->operator()(x); - - reso = retrieveFromCCDB(path + "/" + reso_name, timestamp); - reso->Print(); - LOG(info) << "TPCReso " << reso->operator()(x); - } - - return 0; -} diff --git a/Common/Tools/handleParamTPCResponse.cxx b/Common/Tools/handleParamTPCResponse.cxx new file mode 100644 index 00000000000..14f955806f5 --- /dev/null +++ b/Common/Tools/handleParamTPCResponse.cxx @@ -0,0 +1,300 @@ +// 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 handleParamTPCResponse.cxx +/// \author Jeremy Wilkinson +/// \brief exec for writing and reading Response object + +#include "CCDB/CcdbApi.h" +#include +#include +#include +#include +#include +#include "TFile.h" +#include "Common/Core/PID/PIDResponse.h" +#include "Common/Core/PID/TPCPIDResponse.h" +#include "handleParamBase.h" +#include "Algorithm/RangeTokenizer.h" +using namespace o2::pid::tpc; +namespace bpo = boost::program_options; + +bool initOptionsAndParse(bpo::options_description& options, int argc, char* argv[], bpo::variables_map& arguments) +{ + options.add_options()( + "url,u", bpo::value()->default_value("http://alice-ccdb.cern.ch"), "URL of the CCDB database e.g. http://ccdb-test.cern.ch:8080 or http://alice-ccdb.cern.ch")( + "ccdb-path,c", bpo::value()->default_value("Analysis/PID/TPC"), "CCDB path for storage/retrieval")( + "rct-path", bpo::value()->default_value("RCT/RunInformation"), "path to the ccdb RCT objects for the SOR/EOR timestamps")( + "start,s", bpo::value()->default_value(0), "Start timestamp of object validity. If 0 and runnumber != 0 it will be set to the run SOR")( + "stop,S", bpo::value()->default_value(0), "Stop timestamp of object validity. If 0 and runnumber != 0 it will be set to the run EOR")( + "timestamp,T", bpo::value()->default_value(-1), "Timestamp of the object to retrieve, used in alternative to the run number")( + "runnumber,R", bpo::value()->default_value(0), "Timestamp of the object to retrieve, used in alternative to the timestamp (if 0 using the timestamp)")( + "delete-previous,delete_previous,d", bpo::value()->default_value(0), "Flag to delete previous versions of converter objects in the CCDB before uploading the new one so as to avoid proliferation on CCDB")( + "save-to-file,file,f,o", bpo::value()->default_value(""), "Option to save parametrization to file instead of uploading to ccdb")( + "read-from-file,i", bpo::value()->default_value(""), "Option to get parametrization from a file")( + "objname,n", bpo::value()->default_value("Response"), "Object name to be stored in file")( + "inobjname,n", bpo::value()->default_value("Response"), "Object name to be read from file in 'push' mode")( + "bb0", bpo::value()->default_value(0.03209809958934784f), "Bethe-Bloch parameter 0")( + "bb1", bpo::value()->default_value(19.9768009185791f), "Bethe-Bloch parameter 1")( + "bb2", bpo::value()->default_value(2.5266601063857674e-16f), "Bethe-Bloch parameter 2")( + "bb3", bpo::value()->default_value(2.7212300300598145f), "Bethe-Bloch parameter 3")( + "bb4", bpo::value()->default_value(6.080920219421387f), "Bethe-Bloch parameter 4")( + "reso-param-path", bpo::value()->default_value(""), "Path to resolution parameter file")( + "sigmaGlobal", bpo::value()->default_value("5.43799e-7,0.053044,0.667584,0.0142667,0.00235175,1.22482,2.3501e-7,0.031585"), "Sigma parameters global")( + "paramMIP", bpo::value()->default_value(50.f), "MIP parameter value")( + "paramChargeFactor", bpo::value()->default_value(2.299999952316284f), "Charge factor value")( + "paramMultNormalization", bpo::value()->default_value(11000.), "Multiplicity Normalization")( + "useDefaultParam", bpo::value()->default_value(true), "Use default parametrizatio")( + "dryrun,D", bpo::value()->default_value(0), "Perform a dryrun check before uploading")( + "mode", bpo::value()->default_value(""), "Running mode ('read' from file, 'write' to file, 'pull' from CCDB, 'push' to CCDB)")( + "help,h", "Produce help message."); + try { + bpo::store(parse_command_line(argc, argv, options), arguments); + + // help + if (arguments.count("help")) { + LOG(info) << options; + return false; + } + bpo::notify(arguments); + } catch (const bpo::error& e) { + LOG(error) << e.what() << "\n"; + LOG(error) << "Error parsing command line arguments; Available options:"; + LOG(error) << options; + return false; + } + return true; + +} // initOptionsAndParse + +int main(int argc, char* argv[]) +{ + bpo::options_description options("Allowed options"); + bpo::variables_map arguments; + if (!initOptionsAndParse(options, argc, argv, arguments)) { + return 1; + } + + std::unique_ptr tpc = nullptr; + + const std::string urlCCDB = arguments["url"].as(); + const std::string pathCCDB = arguments["ccdb-path"].as(); + long startTime = arguments["start"].as(); + long endTime = arguments["stop"].as(); + const auto runnumber = arguments["runnumber"].as(); + auto timestamp = arguments["timestamp"].as(); + const int optDelete = arguments["delete-previous"].as(); + + const std::string outFilename = arguments["save-to-file"].as(); + const std::string inFilename = arguments["read-from-file"].as(); + const std::string objname = arguments["objname"].as(); + const std::string inobjname = arguments["inobjname"].as(); + + const float bb0 = arguments["bb0"].as(); + const float bb1 = arguments["bb1"].as(); + const float bb2 = arguments["bb2"].as(); + const float bb3 = arguments["bb3"].as(); + const float bb4 = arguments["bb4"].as(); + const std::string pathResoParam = arguments["reso-param-path"].as(); + const std::string sigmaGlobal = arguments["sigmaGlobal"].as(); + const float mipval = arguments["paramMIP"].as(); + const float chargefacval = arguments["paramChargeFactor"].as(); + const float multNormval = arguments["paramMultNormalization"].as(); + const bool useDefaultParam = arguments["useDefaultParam"].as(); + const std::string optMode = arguments["mode"].as(); + if (optMode.empty()) { + LOG(error) << "--mode must be specified (read, write, pull, push)"; + return 1; + } + // create parameter arrays from commandline options + std::array BBparams = {bb0, bb1, bb2, bb3, bb4}; + + std::vector sigparamsGlobal; + + if (pathResoParam != "") { // Read resolution parameters from file + std::ifstream infile(pathResoParam); + std::string paramstring; + std::vector sigmaParams; + std::getline(infile, paramstring); + std::istringstream ss(paramstring); + sigmaParams.push_back({}); + double param = 0.; + while (ss >> param) { + sigmaParams.push_back(param); + } + sigparamsGlobal = sigmaParams; + } else { + sigparamsGlobal = o2::RangeTokenizer::tokenize(sigmaGlobal); + } + + std::map metadata, headers; + o2::ccdb::CcdbApi api; + if (optMode.compare("push") == 0 || optMode.compare("pull") == 0) { // Initialise CCDB if in push/pull mode + api.init(urlCCDB); + if (!api.isHostReachable()) { + LOG(warning) << "CCDB mode (push/pull) enabled but host " << urlCCDB << " is unreachable."; + return 1; + } + // Init timestamps + if (runnumber != 0) { + const auto rct_path = arguments["rct-path"].as(); + const std::string run_path = Form("%s/%i", rct_path.data(), runnumber); + LOG(info) << "Getting timestamp for run " << runnumber << " from CCDB in path " << run_path; + headers = api.retrieveHeaders(run_path, metadata, -1); + if (headers.count("SOR") == 0) { + LOGF(fatal, "Cannot find run-number to timestamp in path '%s'.", run_path.data()); + } + if (headers.count("SOR") == 0) { + LOGF(fatal, "Cannot find run-number SOR in path '%s'.", run_path.data()); + } + if (headers.count("EOR") == 0) { + LOGF(fatal, "Cannot find run-number EOR in path '%s'.", run_path.data()); + } + timestamp = atol(headers["SOR"].c_str()); // timestamp of the SOR in ms + if (startTime == 0) { + + startTime = atol(headers["SOR"].c_str()); // timestamp of the SOR in ms + LOG(info) << "Setting startTime of object from run number: " << startTime << " -> " << timeStampToHReadble(startTime); + } + if (endTime == 0) { + endTime = atol(headers["EOR"].c_str()); // timestamp of the EOR in ms + LOG(info) << "Setting endTime of object from run number: " << endTime << " -> " << timeStampToHReadble(endTime); + } + } + if (endTime == 0) { // Default value for endTime + endTime = 4108971600000; + } + } + + if (optMode.compare("read") == 0) { // Read existing object from local file + if (inFilename.empty()) { + LOG(error) << "read mode defined with no input file, please set --read-from-file"; + return 1; + } + + TFile fin(inFilename.data(), "READ"); + if (!fin.IsOpen()) { + LOG(error) << "Input file " << inFilename << " could not be read"; + return 1; + } + + tpc.reset(fin.Get(objname.c_str())); + if (!tpc) { + LOG(error) << "Object with name " << objname << " could not be found in file " << inFilename; + return 1; + } + LOG(info) << "Reading existing TPCPIDResponse object " << objname << " from " << inFilename << ":"; + tpc->PrintAll(); + return 0; + } + + else if (optMode.compare("write") == 0 || optMode.compare("push") == 0) // Create new object to write to local file or push to CCDB + { + if (!inFilename.empty()) { // Read from existing file to push to CCDB + LOG(info) << "Reading from existing file to write to CCDB:"; + TFile fin(inFilename.data(), "READ"); + if (!fin.IsOpen()) { + LOG(error) << "Input file " << inFilename << " could not be read"; + return 1; + } + tpc.reset(fin.Get(inobjname.c_str())); + if (!tpc) { + LOG(error) << "Object with name " << objname << " could not be found in file " << inFilename; + return 1; + } + tpc->PrintAll(); + } else { // Create new object if file not specified + LOG(info) << "Creating new TPCPIDResponse object with defined parameters:"; + + tpc.reset(new Response()); + tpc->SetBetheBlochParams(BBparams); + tpc->SetResolutionParams(sigparamsGlobal); + tpc->SetMIP(mipval); + tpc->SetChargeFactor(chargefacval); + tpc->SetMultiplicityNormalization(multNormval); + tpc->SetUseDefaultResolutionParam(useDefaultParam); + tpc->PrintAll(); + } + if (optMode.compare("write") == 0) { + if (outFilename.empty()) { + LOG(error) << "'write' mode specified, but no output filename. Quitting"; + return 1; + } + + LOG(info) << "Writing to output file " << outFilename; + TFile fout(outFilename.data(), "RECREATE"); + if (!fout.IsOpen()) { + LOG(error) << "Output file " << inFilename << " could not be written to"; + return 1; + } + fout.cd(); + // tpc->Print(); + fout.WriteObject(tpc.get(), objname.c_str()); + fout.Close(); + LOG(info) << "File successfully written"; + return 0; + } + + if (optMode.compare("push") == 0) { + LOG(info) << "Attempting to push object to CCDB"; + + if (optDelete) { + api.truncate(pathCCDB); + } + + if (runnumber != 0) { + metadata["runnumber"] = Form("%i", runnumber); + } + // storeOnCCDB(pathCCDB + "/" + objname, metadata, startTime, endTime, tpc.get()); + const auto dryrun = arguments["dryrun"].as(); + if (!dryrun) { + if (arguments["delete-previous"].as()) { + LOG(info) << "Truncating CCDB path " << pathCCDB + "/" + objname; + api.truncate(pathCCDB + "/" + objname); + } + api.storeAsTFileAny(tpc.get(), pathCCDB + "/" + objname, metadata, startTime, endTime); + } else { + LOG(info) << "Dryrunning 'api.storeAsTFileAny(obj," << pathCCDB + "/" + objname << ", metadata, " << startTime << ", " << endTime << ")'"; + } + } + } + + else if (optMode.compare("pull") == 0) { // pull existing from CCDB; write out to file if requested + LOG(info) << "Attempting to pull object from CCDB (" << urlCCDB << "): " << pathCCDB << "/" << objname; + tpc.reset(api.retrieveFromTFileAny(pathCCDB + "/" + objname, metadata, -1, nullptr)); + if (!tpc) { // Quit gracefully if pulling object fails + return 1; + } + + tpc->PrintAll(); + + if (!outFilename.empty()) { + LOG(info) << "Writing pulled object to local file"; + TFile fout(outFilename.data(), "RECREATE"); + if (!fout.IsOpen()) { + LOG(error) << "Output file " << inFilename << " could not be written to"; + return 1; + } + fout.cd(); + // tpc->Print(); + fout.WriteObject(tpc.get(), objname.c_str()); + fout.Close(); + LOG(info) << "File successfully written"; + } + return 0; + } + + else { + LOG(error) << "Invalid mode specified! (must be 'read', 'write', 'pull' or 'push')"; + return 1; + } +} // main \ No newline at end of file diff --git a/PWGEM/Tasks/gammaConversions.cxx b/PWGEM/Tasks/gammaConversions.cxx index 1b5d508f082..aaf9b010c59 100644 --- a/PWGEM/Tasks/gammaConversions.cxx +++ b/PWGEM/Tasks/gammaConversions.cxx @@ -19,7 +19,6 @@ #include "PWGHF/DataModel/HFSecondaryVertex.h" // for BigTracks #include "Common/Core/PID/PIDResponse.h" -#include "Common/Core/PID/PIDTPC.h" #include #include diff --git a/PWGEM/Tasks/gammaConversionsmc.cxx b/PWGEM/Tasks/gammaConversionsmc.cxx index b06a620c801..c54dc2e698c 100644 --- a/PWGEM/Tasks/gammaConversionsmc.cxx +++ b/PWGEM/Tasks/gammaConversionsmc.cxx @@ -19,7 +19,6 @@ #include "PWGHF/DataModel/HFSecondaryVertex.h" // for BigTracks #include "Common/Core/PID/PIDResponse.h" -#include "Common/Core/PID/PIDTPC.h" #include #include From d96e9e26e99407dc944ddbd5bf1fe155e1e2cf1e Mon Sep 17 00:00:00 2001 From: AnnalenaKa Date: Wed, 16 Mar 2022 15:51:42 +0100 Subject: [PATCH 2/2] remove SetParameters from makeTable and use functions handleParamBase --- Common/TableProducer/PID/pidTPC.cxx | 77 +++++++++++++---------- Common/TableProducer/PID/pidTPCFull.cxx | 68 ++++++++++++--------- Common/Tools/CMakeLists.txt | 1 - Common/Tools/handleParamTPCResponse.cxx | 81 +++++-------------------- 4 files changed, 102 insertions(+), 125 deletions(-) diff --git a/Common/TableProducer/PID/pidTPC.cxx b/Common/TableProducer/PID/pidTPC.cxx index cf430a07157..a83798fc4ac 100644 --- a/Common/TableProducer/PID/pidTPC.cxx +++ b/Common/TableProducer/PID/pidTPC.cxx @@ -120,7 +120,7 @@ struct tpcPid { f->GetObject("Response", responseptr); response.SetParameters(responseptr); } catch (...) { - LOGP(info, "Loading the TPC PID Response from file {} failed!", fname); + LOGF(fatal, "Loading the TPC PID Response from file {} failed!", fname); }; } else { useCCDBParam = true; @@ -140,38 +140,53 @@ struct tpcPid { void process(Coll const& collisions, Trks const& tracks, aod::BCsWithTimestamps const&) { - // Check and fill enabled tables - auto makeTable = [&tracks, &collisions, this](const Configurable& flag, auto& table, const o2::track::PID::ID pid) { - if (flag.value == 1) { - // Prepare memory for enabled tables - table.reserve(tracks.size()); - int lastCollisionId = -1; // Last collision ID analysed - for (auto const& trk : tracks) { // Loop on Tracks - auto collision = collisions.iteratorAt(trk.collisionId()); - if (useCCDBParam && ccdbTimestamp.value == 0 && trk.has_collision() && trk.collisionId() != lastCollisionId) { // Updating parametrization only if the initial timestamp is 0 - lastCollisionId = trk.collisionId(); - const auto& bc = collision.bc_as(); - response.SetParameters(ccdb->getForTimeStamp(ccdbPath.value, bc.timestamp())); - } - const float numbersigma = response.GetNumberOfSigma(collision, trk, pid); - aod::pidutils::packInTable(numbersigma, table, - aod::pidtpc_tiny::binned_min, - aod::pidtpc_tiny::binned_max, - aod::pidtpc_tiny::bin_width); - } + auto reserveTable = [&tracks](const Configurable& flag, auto& table) { + if (flag.value != 1) { + return; } + table.reserve(tracks.size()); }; - makeTable(pidEl, tablePIDEl, o2::track::PID::Electron); - makeTable(pidMu, tablePIDMu, o2::track::PID::Muon); - makeTable(pidPi, tablePIDPi, o2::track::PID::Pion); - makeTable(pidKa, tablePIDKa, o2::track::PID::Kaon); - makeTable(pidPr, tablePIDPr, o2::track::PID::Proton); - makeTable(pidDe, tablePIDDe, o2::track::PID::Deuteron); - makeTable(pidTr, tablePIDTr, o2::track::PID::Triton); - makeTable(pidHe, tablePIDHe, o2::track::PID::Helium3); - makeTable(pidAl, tablePIDAl, o2::track::PID::Alpha); + // Prepare memory for enabled tables + reserveTable(pidEl, tablePIDEl); + reserveTable(pidMu, tablePIDMu); + reserveTable(pidPi, tablePIDPi); + reserveTable(pidKa, tablePIDKa); + reserveTable(pidPr, tablePIDPr); + reserveTable(pidDe, tablePIDDe); + reserveTable(pidTr, tablePIDTr); + reserveTable(pidHe, tablePIDHe); + reserveTable(pidAl, tablePIDAl); + int lastCollisionId = -1; // Last collision ID analysed + for (auto const& trk : tracks) { // Loop on Tracks + if (useCCDBParam && ccdbTimestamp.value == 0 && trk.has_collision() && trk.collisionId() != lastCollisionId) { // Updating parametrization only if the initial timestamp is 0 + lastCollisionId = trk.collisionId(); + const auto& bc = collisions.iteratorAt(trk.collisionId()).bc_as(); + response.SetParameters(ccdb->getForTimeStamp(ccdbPath.value, bc.timestamp())); + } + // Check and fill enabled tables + auto makeTable = [&trk, &collisions, this](const Configurable& flag, auto& table, const o2::track::PID::ID pid) { + if (flag.value != 1) { + return; + } + + const float numbersigma = response.GetNumberOfSigma(collisions.iteratorAt(trk.collisionId()), trk, pid); + aod::pidutils::packInTable(numbersigma, table, + aod::pidtpc_tiny::binned_min, + aod::pidtpc_tiny::binned_max, + aod::pidtpc_tiny::bin_width); + }; + makeTable(pidEl, tablePIDEl, o2::track::PID::Electron); + makeTable(pidMu, tablePIDMu, o2::track::PID::Muon); + makeTable(pidPi, tablePIDPi, o2::track::PID::Pion); + makeTable(pidKa, tablePIDKa, o2::track::PID::Kaon); + makeTable(pidPr, tablePIDPr, o2::track::PID::Proton); + makeTable(pidDe, tablePIDDe, o2::track::PID::Deuteron); + makeTable(pidTr, tablePIDTr, o2::track::PID::Triton); + makeTable(pidHe, tablePIDHe, o2::track::PID::Helium3); + makeTable(pidAl, tablePIDAl, o2::track::PID::Alpha); + } } }; diff --git a/Common/TableProducer/PID/pidTPCFull.cxx b/Common/TableProducer/PID/pidTPCFull.cxx index 127cf05b8ce..878e1cbd09e 100644 --- a/Common/TableProducer/PID/pidTPCFull.cxx +++ b/Common/TableProducer/PID/pidTPCFull.cxx @@ -121,7 +121,7 @@ struct tpcPidFull { f->GetObject("Response", responseptr); response.SetParameters(responseptr); } catch (...) { - LOGP(info, "Loading the TPC PID Response from file {} failed!", fname); + LOGF(fatal, "Loading the TPC PID Response from file {} failed!", fname); }; } else { useCCDBParam = true; @@ -141,34 +141,48 @@ struct tpcPidFull { void process(Coll const& collisions, Trks const& tracks, aod::BCsWithTimestamps const&) { - // Check and fill enabled tables - auto makeTable = [&tracks, &collisions, this](const Configurable& flag, auto& table, const o2::track::PID::ID pid) { - if (flag.value == 1) { // Loop on Tracks - // Prepare memory for enabled tables - table.reserve(tracks.size()); - int lastCollisionId = -1; // Last collision ID analysed - for (auto const& trk : tracks) { - auto collision = collisions.iteratorAt(trk.collisionId()); // Loop on Tracks - if (useCCDBParam && ccdbTimestamp.value == 0 && trk.has_collision() && trk.collisionId() != lastCollisionId) { // Updating parametrization only if the initial timestamp is 0 - lastCollisionId = trk.collisionId(); - const auto& bc = collision.bc_as(); - response.SetParameters(ccdb->getForTimeStamp(ccdbPath.value, bc.timestamp())); - } - table(response.GetExpectedSigma(collision, trk, pid), - response.GetNumberOfSigma(collision, trk, pid)); - } + auto reserveTable = [&tracks](const Configurable& flag, auto& table) { + if (flag.value != 1) { + return; } + table.reserve(tracks.size()); }; - // const o2::pid::tpc::Response& response; - makeTable(pidEl, tablePIDEl, o2::track::PID::Electron); - makeTable(pidMu, tablePIDMu, o2::track::PID::Muon); - makeTable(pidPi, tablePIDPi, o2::track::PID::Pion); - makeTable(pidKa, tablePIDKa, o2::track::PID::Kaon); - makeTable(pidPr, tablePIDPr, o2::track::PID::Proton); - makeTable(pidDe, tablePIDDe, o2::track::PID::Deuteron); - makeTable(pidTr, tablePIDTr, o2::track::PID::Triton); - makeTable(pidHe, tablePIDHe, o2::track::PID::Helium3); - makeTable(pidAl, tablePIDAl, o2::track::PID::Alpha); + // Prepare memory for enabled tables + reserveTable(pidEl, tablePIDEl); + reserveTable(pidMu, tablePIDMu); + reserveTable(pidPi, tablePIDPi); + reserveTable(pidKa, tablePIDKa); + reserveTable(pidPr, tablePIDPr); + reserveTable(pidDe, tablePIDDe); + reserveTable(pidTr, tablePIDTr); + reserveTable(pidHe, tablePIDHe); + reserveTable(pidAl, tablePIDAl); + int lastCollisionId = -1; // Last collision ID analysed + for (auto const& trk : tracks) { // Loop on Tracks + if (useCCDBParam && ccdbTimestamp.value == 0 && trk.has_collision() && trk.collisionId() != lastCollisionId) { // Updating parametrization only if the initial timestamp is 0 + lastCollisionId = trk.collisionId(); + const auto& bc = collisions.iteratorAt(trk.collisionId()).bc_as(); + response.SetParameters(ccdb->getForTimeStamp(ccdbPath.value, bc.timestamp())); + } + // Check and fill enabled tables + auto makeTable = [&trk, &collisions, this](const Configurable& flag, auto& table, const o2::track::PID::ID pid) { + if (flag.value != 1) { + return; + } + table(response.GetExpectedSigma(collisions.iteratorAt(trk.collisionId()), trk, pid), + response.GetNumberOfSigma(collisions.iteratorAt(trk.collisionId()), trk, pid)); + }; + // const o2::pid::tpc::Response& response; + makeTable(pidEl, tablePIDEl, o2::track::PID::Electron); + makeTable(pidMu, tablePIDMu, o2::track::PID::Muon); + makeTable(pidPi, tablePIDPi, o2::track::PID::Pion); + makeTable(pidKa, tablePIDKa, o2::track::PID::Kaon); + makeTable(pidPr, tablePIDPr, o2::track::PID::Proton); + makeTable(pidDe, tablePIDDe, o2::track::PID::Deuteron); + makeTable(pidTr, tablePIDTr, o2::track::PID::Triton); + makeTable(pidHe, tablePIDHe, o2::track::PID::Helium3); + makeTable(pidAl, tablePIDAl, o2::track::PID::Alpha); + } } }; diff --git a/Common/Tools/CMakeLists.txt b/Common/Tools/CMakeLists.txt index 4bcad206458..9de705d1fde 100644 --- a/Common/Tools/CMakeLists.txt +++ b/Common/Tools/CMakeLists.txt @@ -1,4 +1,3 @@ - # 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. diff --git a/Common/Tools/handleParamTPCResponse.cxx b/Common/Tools/handleParamTPCResponse.cxx index 14f955806f5..be91239c65c 100644 --- a/Common/Tools/handleParamTPCResponse.cxx +++ b/Common/Tools/handleParamTPCResponse.cxx @@ -14,21 +14,16 @@ /// \author Jeremy Wilkinson /// \brief exec for writing and reading Response object -#include "CCDB/CcdbApi.h" -#include -#include #include #include #include #include "TFile.h" -#include "Common/Core/PID/PIDResponse.h" #include "Common/Core/PID/TPCPIDResponse.h" #include "handleParamBase.h" #include "Algorithm/RangeTokenizer.h" using namespace o2::pid::tpc; -namespace bpo = boost::program_options; -bool initOptionsAndParse(bpo::options_description& options, int argc, char* argv[], bpo::variables_map& arguments) +bool initOptionsAndParse(bpo::options_description& options, int argc, char* argv[]) { options.add_options()( "url,u", bpo::value()->default_value("http://alice-ccdb.cern.ch"), "URL of the CCDB database e.g. http://ccdb-test.cern.ch:8080 or http://alice-ccdb.cern.ch")( @@ -79,17 +74,16 @@ bool initOptionsAndParse(bpo::options_description& options, int argc, char* argv int main(int argc, char* argv[]) { bpo::options_description options("Allowed options"); - bpo::variables_map arguments; - if (!initOptionsAndParse(options, argc, argv, arguments)) { + if (!initOptionsAndParse(options, argc, argv)) { //, arguments return 1; } - std::unique_ptr tpc = nullptr; + Response* tpc = nullptr; const std::string urlCCDB = arguments["url"].as(); - const std::string pathCCDB = arguments["ccdb-path"].as(); - long startTime = arguments["start"].as(); - long endTime = arguments["stop"].as(); + const auto pathCCDB = arguments["ccdb-path"].as(); + auto startTime = arguments["start"].as(); + auto endTime = arguments["stop"].as(); const auto runnumber = arguments["runnumber"].as(); auto timestamp = arguments["timestamp"].as(); const int optDelete = arguments["delete-previous"].as(); @@ -136,43 +130,13 @@ int main(int argc, char* argv[]) sigparamsGlobal = o2::RangeTokenizer::tokenize(sigmaGlobal); } - std::map metadata, headers; - o2::ccdb::CcdbApi api; if (optMode.compare("push") == 0 || optMode.compare("pull") == 0) { // Initialise CCDB if in push/pull mode api.init(urlCCDB); if (!api.isHostReachable()) { LOG(warning) << "CCDB mode (push/pull) enabled but host " << urlCCDB << " is unreachable."; return 1; } - // Init timestamps - if (runnumber != 0) { - const auto rct_path = arguments["rct-path"].as(); - const std::string run_path = Form("%s/%i", rct_path.data(), runnumber); - LOG(info) << "Getting timestamp for run " << runnumber << " from CCDB in path " << run_path; - headers = api.retrieveHeaders(run_path, metadata, -1); - if (headers.count("SOR") == 0) { - LOGF(fatal, "Cannot find run-number to timestamp in path '%s'.", run_path.data()); - } - if (headers.count("SOR") == 0) { - LOGF(fatal, "Cannot find run-number SOR in path '%s'.", run_path.data()); - } - if (headers.count("EOR") == 0) { - LOGF(fatal, "Cannot find run-number EOR in path '%s'.", run_path.data()); - } - timestamp = atol(headers["SOR"].c_str()); // timestamp of the SOR in ms - if (startTime == 0) { - - startTime = atol(headers["SOR"].c_str()); // timestamp of the SOR in ms - LOG(info) << "Setting startTime of object from run number: " << startTime << " -> " << timeStampToHReadble(startTime); - } - if (endTime == 0) { - endTime = atol(headers["EOR"].c_str()); // timestamp of the EOR in ms - LOG(info) << "Setting endTime of object from run number: " << endTime << " -> " << timeStampToHReadble(endTime); - } - } - if (endTime == 0) { // Default value for endTime - endTime = 4108971600000; - } + setupTimestamps(timestamp, startTime, endTime); } if (optMode.compare("read") == 0) { // Read existing object from local file @@ -187,7 +151,7 @@ int main(int argc, char* argv[]) return 1; } - tpc.reset(fin.Get(objname.c_str())); + tpc = fin.Get(objname.c_str()); if (!tpc) { LOG(error) << "Object with name " << objname << " could not be found in file " << inFilename; return 1; @@ -206,7 +170,7 @@ int main(int argc, char* argv[]) LOG(error) << "Input file " << inFilename << " could not be read"; return 1; } - tpc.reset(fin.Get(inobjname.c_str())); + tpc = fin.Get(inobjname.c_str()); if (!tpc) { LOG(error) << "Object with name " << objname << " could not be found in file " << inFilename; return 1; @@ -215,7 +179,7 @@ int main(int argc, char* argv[]) } else { // Create new object if file not specified LOG(info) << "Creating new TPCPIDResponse object with defined parameters:"; - tpc.reset(new Response()); + tpc = new Response(); tpc->SetBetheBlochParams(BBparams); tpc->SetResolutionParams(sigparamsGlobal); tpc->SetMIP(mipval); @@ -237,8 +201,7 @@ int main(int argc, char* argv[]) return 1; } fout.cd(); - // tpc->Print(); - fout.WriteObject(tpc.get(), objname.c_str()); + fout.WriteObject(tpc, objname.c_str()); fout.Close(); LOG(info) << "File successfully written"; return 0; @@ -250,30 +213,17 @@ int main(int argc, char* argv[]) if (optDelete) { api.truncate(pathCCDB); } - + std::map metadata; if (runnumber != 0) { metadata["runnumber"] = Form("%i", runnumber); } - // storeOnCCDB(pathCCDB + "/" + objname, metadata, startTime, endTime, tpc.get()); - const auto dryrun = arguments["dryrun"].as(); - if (!dryrun) { - if (arguments["delete-previous"].as()) { - LOG(info) << "Truncating CCDB path " << pathCCDB + "/" + objname; - api.truncate(pathCCDB + "/" + objname); - } - api.storeAsTFileAny(tpc.get(), pathCCDB + "/" + objname, metadata, startTime, endTime); - } else { - LOG(info) << "Dryrunning 'api.storeAsTFileAny(obj," << pathCCDB + "/" + objname << ", metadata, " << startTime << ", " << endTime << ")'"; - } + storeOnCCDB(pathCCDB + "/" + objname, metadata, startTime, endTime, tpc); } } else if (optMode.compare("pull") == 0) { // pull existing from CCDB; write out to file if requested LOG(info) << "Attempting to pull object from CCDB (" << urlCCDB << "): " << pathCCDB << "/" << objname; - tpc.reset(api.retrieveFromTFileAny(pathCCDB + "/" + objname, metadata, -1, nullptr)); - if (!tpc) { // Quit gracefully if pulling object fails - return 1; - } + tpc = retrieveFromCCDB(pathCCDB + "/" + objname, timestamp); tpc->PrintAll(); @@ -285,8 +235,7 @@ int main(int argc, char* argv[]) return 1; } fout.cd(); - // tpc->Print(); - fout.WriteObject(tpc.get(), objname.c_str()); + fout.WriteObject(tpc, objname.c_str()); fout.Close(); LOG(info) << "File successfully written"; }