diff --git a/Common/Core/PID/TPCPIDResponse.h b/Common/Core/PID/TPCPIDResponse.h index a99c0bb5cea..58dfebcf58b 100644 --- a/Common/Core/PID/TPCPIDResponse.h +++ b/Common/Core/PID/TPCPIDResponse.h @@ -23,6 +23,7 @@ #include "Framework/Logger.h" // O2 includes #include "ReconstructionDataFormats/PID.h" +#include "Framework/DataTypes.h" #include "DataFormatsTPC/BetheBlochAleph.h" namespace o2::pid::tpc @@ -44,6 +45,7 @@ class Response void SetMIP(const float mip) { mMIP = mip; } void SetChargeFactor(const float chargeFactor) { mChargeFactor = chargeFactor; } void SetMultiplicityNormalization(const float multNormalization) { mMultNormalization = multNormalization; } + void SetNClNormalization(const float nclnorm) { nClNorm = nclnorm; } void SetUseDefaultResolutionParam(const bool useDefault) { mUseDefaultResolutionParam = useDefault; } void SetParameters(const Response* response) { @@ -51,6 +53,7 @@ class Response mResolutionParamsDefault = response->GetResolutionParamsDefault(); mResolutionParams = response->GetResolutionParams(); mMIP = response->GetMIP(); + nClNorm = response->GetNClNormalization(); mChargeFactor = response->GetChargeFactor(); mMultNormalization = response->GetMultiplicityNormalization(); mUseDefaultResolutionParam = response->GetUseDefaultResolutionParam(); @@ -60,6 +63,7 @@ class Response const std::array GetResolutionParamsDefault() const { return mResolutionParamsDefault; } const std::vector GetResolutionParams() const { return mResolutionParams; } const float GetMIP() const { return mMIP; } + const float GetNClNormalization() const { return nClNorm; } const float GetChargeFactor() const { return mChargeFactor; } const float GetMultiplicityNormalization() const { return mMultNormalization; } const bool GetUseDefaultResolutionParam() const { return mUseDefaultResolutionParam; } @@ -89,8 +93,9 @@ class Response float mChargeFactor = 2.299999952316284f; float mMultNormalization = 11000.; bool mUseDefaultResolutionParam = true; + float nClNorm = 152.f; - ClassDefNV(Response, 2); + ClassDefNV(Response, 3); }; // class Response @@ -98,21 +103,27 @@ class Response template inline float Response::GetExpectedSignal(const TrackType& track, const o2::track::PID::ID id) const { + if (!track.hasTPC()) { + return -999.f; + } 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; + return bethe >= 0.f ? bethe : -999.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 { + if (!track.hasTPC()) { + return -999.f; + } 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; + reso >= 0.f ? resolution = reso : resolution = -999.f; } else { - const double ncl = 159. / track.tpcNClsFound(); // + const double ncl = nClNorm / track.tpcNClsFound(); // const double p = track.tpcInnerParam(); const double mass = o2::track::pid_constants::sMasses[id]; const double bg = p / mass; @@ -122,7 +133,7 @@ inline float Response::GetExpectedSigma(const CollisionType& collision, const Tr 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; + reso >= 0.f ? resolution = reso : resolution = -999.f; } return resolution; } @@ -131,6 +142,15 @@ inline float Response::GetExpectedSigma(const CollisionType& collision, const Tr template inline float Response::GetNumberOfSigma(const CollisionType& collision, const TrackType& trk, const o2::track::PID::ID id) const { + if (GetExpectedSigma(collision, trk, id) < 0.) { + return -999.f; + } + if (GetExpectedSignal(trk, id) < 0.) { + return -999.f; + } + if (!trk.hasTPC()) { + return -999.f; + } return ((trk.tpcSignal() - GetExpectedSignal(trk, id)) / GetExpectedSigma(collision, trk, id)); } @@ -138,6 +158,12 @@ inline float Response::GetNumberOfSigma(const CollisionType& collision, const Tr template inline float Response::GetSignalDelta(const TrackType& trk, const o2::track::PID::ID id) const { + if (GetExpectedSignal(trk, id) < 0.) { + return -999.f; + } + if (!trk.hasTPC()) { + return -999.f; + } return (trk.tpcSignal() - GetExpectedSignal(trk, id)); } @@ -173,6 +199,7 @@ inline void Response::PrintAll() const LOGP(info, "mMIP = {}", mMIP); LOGP(info, "mChargeFactor = {}", mChargeFactor); LOGP(info, "mMultNormalization = {}", mMultNormalization); + LOGP(info, "nClNorm = {}", nClNorm); } } // namespace o2::pid::tpc diff --git a/Common/DataModel/PIDResponse.h b/Common/DataModel/PIDResponse.h index 511af101cd5..3865c1a1c8a 100644 --- a/Common/DataModel/PIDResponse.h +++ b/Common/DataModel/PIDResponse.h @@ -888,6 +888,25 @@ DECLARE_SOA_TABLE(pidTOFAl, "AOD", "pidTOFAl", //! Table of the TOF response wit namespace pidtpc { // Expected signals +DECLARE_SOA_DYNAMIC_COLUMN(TPCExpSignalEl, tpcExpSignalEl, //! Expected signal with the TPC detector for electron + [](float nsigma, float sigma, float tpcsignal) -> float { return tpcsignal - nsigma * sigma; }); +DECLARE_SOA_DYNAMIC_COLUMN(TPCExpSignalMu, tpcExpSignalMu, //! Expected signal with the TPC detector for muon + [](float nsigma, float sigma, float tpcsignal) -> float { return tpcsignal - nsigma * sigma; }); +DECLARE_SOA_DYNAMIC_COLUMN(TPCExpSignalPi, tpcExpSignalPi, //! Expected signal with the TPC detector for pion + [](float nsigma, float sigma, float tpcsignal) -> float { return tpcsignal - nsigma * sigma; }); +DECLARE_SOA_DYNAMIC_COLUMN(TPCExpSignalKa, tpcExpSignalKa, //! Expected signal with the TPC detector for kaon + [](float nsigma, float sigma, float tpcsignal) -> float { return tpcsignal - nsigma * sigma; }); +DECLARE_SOA_DYNAMIC_COLUMN(TPCExpSignalPr, tpcExpSignalPr, //! Expected signal with the TPC detector for proton + [](float nsigma, float sigma, float tpcsignal) -> float { return tpcsignal - nsigma * sigma; }); +DECLARE_SOA_DYNAMIC_COLUMN(TPCExpSignalDe, tpcExpSignalDe, //! Expected signal with the TPC detector for deuteron + [](float nsigma, float sigma, float tpcsignal) -> float { return tpcsignal - nsigma * sigma; }); +DECLARE_SOA_DYNAMIC_COLUMN(TPCExpSignalTr, tpcExpSignalTr, //! Expected signal with the TPC detector for triton + [](float nsigma, float sigma, float tpcsignal) -> float { return tpcsignal - nsigma * sigma; }); +DECLARE_SOA_DYNAMIC_COLUMN(TPCExpSignalHe, tpcExpSignalHe, //! Expected signal with the TPC detector for helium3 + [](float nsigma, float sigma, float tpcsignal) -> float { return tpcsignal - nsigma * sigma; }); +DECLARE_SOA_DYNAMIC_COLUMN(TPCExpSignalAl, tpcExpSignalAl, //! Expected signal with the TPC detector for alpha + [](float nsigma, float sigma, float tpcsignal) -> float { return tpcsignal - nsigma * sigma; }); +// Expected signals difference DECLARE_SOA_DYNAMIC_COLUMN(TPCExpSignalDiffEl, tpcExpSignalDiffEl, //! Difference between signal and expected for electron [](float nsigma, float sigma) -> float { return nsigma * sigma; }); DECLARE_SOA_DYNAMIC_COLUMN(TPCExpSignalDiffMu, tpcExpSignalDiffMu, //! Difference between signal and expected for muon @@ -968,23 +987,23 @@ DEFINE_UNWRAP_NSIGMA_COLUMN(TPCNSigmaAl, tpcNSigmaAl); //! Unwrapped (float) nsi // Per particle tables DECLARE_SOA_TABLE(pidTPCFullEl, "AOD", "pidTPCFullEl", //! Table of the TPC (full) response with expected signal, expected resolution and Nsigma for electron - pidtpc::TPCExpSignalDiffEl, pidtpc::TPCExpSigmaEl, pidtpc::TPCNSigmaEl); + pidtpc::TPCExpSignalEl, pidtpc::TPCExpSignalDiffEl, pidtpc::TPCExpSigmaEl, pidtpc::TPCNSigmaEl); DECLARE_SOA_TABLE(pidTPCFullMu, "AOD", "pidTPCFullMu", //! Table of the TPC (full) response with expected signal, expected resolution and Nsigma for muon - pidtpc::TPCExpSignalDiffMu, pidtpc::TPCExpSigmaMu, pidtpc::TPCNSigmaMu); + pidtpc::TPCExpSignalMu, pidtpc::TPCExpSignalDiffMu, pidtpc::TPCExpSigmaMu, pidtpc::TPCNSigmaMu); DECLARE_SOA_TABLE(pidTPCFullPi, "AOD", "pidTPCFullPi", //! Table of the TPC (full) response with expected signal, expected resolution and Nsigma for pion - pidtpc::TPCExpSignalDiffPi, pidtpc::TPCExpSigmaPi, pidtpc::TPCNSigmaPi); + pidtpc::TPCExpSignalPi, pidtpc::TPCExpSignalDiffPi, pidtpc::TPCExpSigmaPi, pidtpc::TPCNSigmaPi); DECLARE_SOA_TABLE(pidTPCFullKa, "AOD", "pidTPCFullKa", //! Table of the TPC (full) response with expected signal, expected resolution and Nsigma for kaon - pidtpc::TPCExpSignalDiffKa, pidtpc::TPCExpSigmaKa, pidtpc::TPCNSigmaKa); + pidtpc::TPCExpSignalKa, pidtpc::TPCExpSignalDiffKa, pidtpc::TPCExpSigmaKa, pidtpc::TPCNSigmaKa); DECLARE_SOA_TABLE(pidTPCFullPr, "AOD", "pidTPCFullPr", //! Table of the TPC (full) response with expected signal, expected resolution and Nsigma for proton - pidtpc::TPCExpSignalDiffPr, pidtpc::TPCExpSigmaPr, pidtpc::TPCNSigmaPr); + pidtpc::TPCExpSignalPr, pidtpc::TPCExpSignalDiffPr, pidtpc::TPCExpSigmaPr, pidtpc::TPCNSigmaPr); DECLARE_SOA_TABLE(pidTPCFullDe, "AOD", "pidTPCFullDe", //! Table of the TPC (full) response with expected signal, expected resolution and Nsigma for deuteron - pidtpc::TPCExpSignalDiffDe, pidtpc::TPCExpSigmaDe, pidtpc::TPCNSigmaDe); + pidtpc::TPCExpSignalDe, pidtpc::TPCExpSignalDiffDe, pidtpc::TPCExpSigmaDe, pidtpc::TPCNSigmaDe); DECLARE_SOA_TABLE(pidTPCFullTr, "AOD", "pidTPCFullTr", //! Table of the TPC (full) response with expected signal, expected resolution and Nsigma for triton - pidtpc::TPCExpSignalDiffTr, pidtpc::TPCExpSigmaTr, pidtpc::TPCNSigmaTr); + pidtpc::TPCExpSignalTr, pidtpc::TPCExpSignalDiffTr, pidtpc::TPCExpSigmaTr, pidtpc::TPCNSigmaTr); DECLARE_SOA_TABLE(pidTPCFullHe, "AOD", "pidTPCFullHe", //! Table of the TPC (full) response with expected signal, expected resolution and Nsigma for helium3 - pidtpc::TPCExpSignalDiffHe, pidtpc::TPCExpSigmaHe, pidtpc::TPCNSigmaHe); + pidtpc::TPCExpSignalHe, pidtpc::TPCExpSignalDiffHe, pidtpc::TPCExpSigmaHe, pidtpc::TPCNSigmaHe); DECLARE_SOA_TABLE(pidTPCFullAl, "AOD", "pidTPCFullAl", //! Table of the TPC (full) response with expected signal, expected resolution and Nsigma for alpha - pidtpc::TPCExpSignalDiffAl, pidtpc::TPCExpSigmaAl, pidtpc::TPCNSigmaAl); + pidtpc::TPCExpSignalAl, pidtpc::TPCExpSignalDiffAl, pidtpc::TPCExpSigmaAl, pidtpc::TPCNSigmaAl); // Tiny size tables DECLARE_SOA_TABLE(pidTPCEl, "AOD", "pidTPCEl", //! Table of the TPC response with binned Nsigma for electron diff --git a/Common/TableProducer/PID/pidTPC.cxx b/Common/TableProducer/PID/pidTPC.cxx index c60fc12b4c2..6e2cb9a391d 100644 --- a/Common/TableProducer/PID/pidTPC.cxx +++ b/Common/TableProducer/PID/pidTPC.cxx @@ -158,6 +158,7 @@ struct tpcPid { enableNetworkOptimizations.value); network = temp_net; network.evalNetwork(std::vector(network.getInputDimensions(), 1.)); // This is an initialisation and might reduce the overhead of the model + network.SetNClNormalization(response.GetNClNormalization()); } } } @@ -186,6 +187,7 @@ struct tpcPid { reserveTable(pidAl, tablePIDAl); std::vector network_prediction; + const float nNclNormalization = response.GetNClNormalization(); if (useNetworkCorrection) { @@ -230,7 +232,7 @@ struct tpcPid { track_properties[counter_track_props + 2] = trk.signed1Pt(); track_properties[counter_track_props + 3] = o2::track::pid_constants::sMasses[i]; track_properties[counter_track_props + 4] = collisions.iteratorAt(trk.collisionId()).multTPC() / 11000.; - track_properties[counter_track_props + 5] = std::sqrt(159. / trk.tpcNClsFound()); + track_properties[counter_track_props + 5] = std::sqrt(nNclNormalization / trk.tpcNClsFound()); counter_track_props += input_dimensions; } diff --git a/Common/TableProducer/PID/pidTPCFull.cxx b/Common/TableProducer/PID/pidTPCFull.cxx index 4dc004b00d4..8d339caf38c 100644 --- a/Common/TableProducer/PID/pidTPCFull.cxx +++ b/Common/TableProducer/PID/pidTPCFull.cxx @@ -158,6 +158,7 @@ struct tpcPidFull { enableNetworkOptimizations.value); network = temp_net; network.evalNetwork(std::vector(network.getInputDimensions(), 1.)); // This is an initialisation and might reduce the overhead of the model + network.SetNClNormalization(response.GetNClNormalization()); } } } @@ -214,7 +215,7 @@ struct tpcPidFull { const unsigned long prediction_size = output_dimensions * tracks_size; network_prediction = std::vector(prediction_size * 9); // For each mass hypotheses - + const float nNclNormalization = response.GetNClNormalization(); float duration_network = 0; std::vector track_properties(track_prop_size); @@ -230,7 +231,7 @@ struct tpcPidFull { track_properties[counter_track_props + 2] = trk.signed1Pt(); track_properties[counter_track_props + 3] = o2::track::pid_constants::sMasses[i]; track_properties[counter_track_props + 4] = collisions.iteratorAt(trk.collisionId()).multTPC() / 11000.; - track_properties[counter_track_props + 5] = std::sqrt(159. / trk.tpcNClsFound()); + track_properties[counter_track_props + 5] = std::sqrt(nNclNormalization / trk.tpcNClsFound()); counter_track_props += input_dimensions; } diff --git a/Common/TableProducer/PID/pidTPCML.cxx b/Common/TableProducer/PID/pidTPCML.cxx index 6207065f1ec..060ff0a5a8e 100644 --- a/Common/TableProducer/PID/pidTPCML.cxx +++ b/Common/TableProducer/PID/pidTPCML.cxx @@ -124,7 +124,7 @@ std::array Network::createInputFromTrack(const C& collision_it, const const float signed1Pt = track.signed1Pt(); const float mass = o2::track::pid_constants::sMasses[id]; const float multTPC = collision_it.multTPC() / 11000.; - const float ncl = std::sqrt(159. / track.tpcNClsFound()); + const float ncl = std::sqrt(nClNorm / track.tpcNClsFound()); return {p, tgl, signed1Pt, mass, multTPC, ncl}; } diff --git a/Common/TableProducer/PID/pidTPCML.h b/Common/TableProducer/PID/pidTPCML.h index 245989b80a7..8ea794346fc 100644 --- a/Common/TableProducer/PID/pidTPCML.h +++ b/Common/TableProducer/PID/pidTPCML.h @@ -46,12 +46,14 @@ class Network template std::array createInputFromTrack(const C&, const T&, const uint8_t) const; // create a std::vector with all the inputs for the network std::vector createTensor(std::array) const; // create a std::vector (= ONNX tensor) for model input - float* evalNetwork(std::vector); // evaluate the network on a std::vector (= ONNX tensor) - float* evalNetwork(std::vector); // evaluate the network on a std::vector + float* evalNetwork(std::vector); // evaluate the network on a std::vector (= ONNX tensor) + float* evalNetwork(std::vector); // evaluate the network on a std::vector // Getters & Setters int getInputDimensions() const { return mInputShapes[0][1]; }; int getOutputDimensions() const { return mOutputShapes[0][1]; }; + void SetNClNormalization(const float nclnorm) { nClNorm = nclnorm; } + const float GetNClNormalization() const { return nClNorm; } private: // Environment variables for the ONNX runtime @@ -65,6 +67,8 @@ class Network std::vector mOutputNames; std::vector> mOutputShapes; + float nClNorm = 152.f; + // Internal function for printing the shape of tensors: See https://github.com/saganatt/PID_ML_in_O2 or O2Physics/Tools/PIDML/simpleApplyPidOnnxModel.cxx std::string printShape(const std::vector& v); @@ -75,4 +79,4 @@ class Network } // namespace o2::pid::tpc -#endif // O2_PID_TPC_ML_H_ +#endif // O2_PID_TPC_ML_H_ \ No newline at end of file diff --git a/Common/Tools/CMakeLists.txt b/Common/Tools/CMakeLists.txt index 3e8898ca1a4..d7285a428a1 100644 --- a/Common/Tools/CMakeLists.txt +++ b/Common/Tools/CMakeLists.txt @@ -24,4 +24,5 @@ o2physics_add_executable(pidparam-tof-reso 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/handleParamTPCResponse.cxx b/Common/Tools/handleParamTPCResponse.cxx index f72c72b1d60..d7430dd9c46 100644 --- a/Common/Tools/handleParamTPCResponse.cxx +++ b/Common/Tools/handleParamTPCResponse.cxx @@ -41,6 +41,7 @@ bool initOptionsAndParse(bpo::options_description& options, int argc, char* argv "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")( + "paramnClNormalization", bpo::value()->default_value(152.), "Maximum nClusters for normalisation (159 for run 2, 152 for run 3)")( "useDefaultParam", bpo::value()->default_value(true), "Use default sigma parametrisation")( "mode", bpo::value()->default_value(""), "Running mode ('read' from file, 'write' to file, 'pull' from CCDB, 'push' to CCDB)")( "help,h", "Produce help message."); @@ -93,6 +94,7 @@ int main(int argc, char* argv[]) const float mipval = arguments["paramMIP"].as(); const float chargefacval = arguments["paramChargeFactor"].as(); const float multNormval = arguments["paramMultNormalization"].as(); + const float nClNormval = arguments["paramnClNormalization"].as(); const bool useDefaultParam = arguments["useDefaultParam"].as(); const std::string optMode = arguments["mode"].as(); if (optMode.empty()) { @@ -172,6 +174,7 @@ int main(int argc, char* argv[]) tpc->SetMIP(mipval); tpc->SetChargeFactor(chargefacval); tpc->SetMultiplicityNormalization(multNormval); + tpc->SetNClNormalization(nClNormval); tpc->SetUseDefaultResolutionParam(useDefaultParam); tpc->PrintAll(); } diff --git a/DPG/Tasks/CMakeLists.txt b/DPG/Tasks/CMakeLists.txt index bc2663ea4ca..3aac657b75e 100644 --- a/DPG/Tasks/CMakeLists.txt +++ b/DPG/Tasks/CMakeLists.txt @@ -10,5 +10,5 @@ # or submit itself to any jurisdiction. add_subdirectory(AOTTrack) +add_subdirectory(TPC) add_subdirectory(FDD) - diff --git a/DPG/Tasks/TPC/CMakeLists.txt b/DPG/Tasks/TPC/CMakeLists.txt new file mode 100644 index 00000000000..80e0ec33a05 --- /dev/null +++ b/DPG/Tasks/TPC/CMakeLists.txt @@ -0,0 +1,15 @@ +# 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. + +o2physics_add_dpl_workflow(pid-tpc-skimscreation + SOURCES tpcSkimsTableCreator.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore + COMPONENT_NAME Analysis) \ No newline at end of file diff --git a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx new file mode 100644 index 00000000000..28e78a0c348 --- /dev/null +++ b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx @@ -0,0 +1,558 @@ +// 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 tpcSkimsTableCreator.cxx +/// \brief Task to produce table with clean selections for TPC PID calibration +/// +/// \author Annalena Kalteyer +/// \author Christian Sonnabend +/// \author Jeremy Wilkinson +/// O2 +#include +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/runDataProcessing.h" +/// O2Physics +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/PIDResponse.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Common/DataModel/StrangenessTables.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/EventSelection.h" + +#include "tpcSkimsTableCreator.h" +/// ROOT +#include "TRandom3.h" +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::track; +using namespace o2::dataformats; + +struct TreeWriterTpcV0 { + + using Trks = soa::Join; + using Coll = soa::Join; + + /// Tables to be produced + Produces rowTPCTree; + + /// Configurables + Configurable cutPAV0{"cutPAV0", 0., "Cut on the cos(pointing angle) of the decay"}; + Configurable nSigmaTOFdautrack{"nSigmaTOFdautrack", 5., "n-sigma TOF cut on the daughter tracks. Set 0 to switch it off."}; + Configurable nClNorm{"nClNorm", 152., "Number of cluster normalization. Run 2: 159, Run 3 152"}; + Configurable applyEvSel{"applyEvSel", 2, "Flag to apply rapidity cut: 0 -> no event selection, 1 -> Run 2 event selection, 2 -> Run 3 event selection"}; + /// Configurables downsampling + Configurable dwnSmplFactor_Pi{"dwnSmplFactor_Pi", 1., "downsampling factor for pions, default fraction to keep is 1."}; + Configurable dwnSmplFactor_Pr{"dwnSmplFactor_Pr", 1., "downsampling factor for protons, default fraction to keep is 1."}; + Configurable dwnSmplFactor_El{"dwnSmplFactor_El", 1., "downsampling factor for electrons, default fraction to keep is 1."}; + Configurable sqrtSNN{"sqrt_s_NN", 0., "sqrt(s_NN), used for downsampling with the Tsallis distribution"}; + Configurable downsamplingTsalisPions{"downsamplingTsalisPions", -1., "Downsampling factor to reduce the number of pions"}; + Configurable downsamplingTsalisProtons{"downsamplingTsalisProtons", -1., "Downsampling factor to reduce the number of protons"}; + Configurable downsamplingTsalisElectrons{"downsamplingTsalisElectrons", -1., "Downsampling factor to reduce the number of electrons"}; + /// Configurables kaon + Configurable invariantMassCutK0Short{"invariantMassCutK0Short", 0.5, "Mass cut for K0short"}; + Configurable cutQTK0min{"cutQTK0min", 0.1075, "Minimum qt for K0short"}; + Configurable cutQTK0max{"cutQTK0max", 0.215, "Minimum qt for K0short"}; + Configurable cutAPK0pinner{"cutAPK0pinner", 0.199, "inner momentum for K0short"}; + Configurable cutAPK0pouter{"cutAPK0pouter", 0.215, "outer momentum for K0short"}; + Configurable cutAPK0alinner{"cutAPK0alinner", 0.8, "inner alpha for K0short"}; + Configurable cutAPK0alouter{"cutAPK0alouter", 0.9, "outer alpha for K0short"}; + /// Configurables lambda/antilambda + Configurable invariantMassCutLambda{"invariantMassCutLambda", 0.5, "Mass cut for Lambda and Anti-Lambda"}; + Configurable cutQTmin{"cutQTmin", 0.03, "the minimum value of qt for Lambda"}; + Configurable cutAlphaminL{"cutAlphaminL", 0.35, "minimum alpha for Lambda"}; + Configurable cutAlphamaxL{"cutAlphamaxL", 0.7, "maximum alpha for Lambda"}; + Configurable cutAlphaminAL{"cutAlphaminAL", -0.7, "minimum alpha for Anti-Lambda"}; + Configurable cutAlphamaxAL{"cutAlphamaxAL", -0.35, "maximum alpha for Anti-Lambda"}; + Configurable cutAPL0{"cutAPL0", 0.10, "cutAPL par0"}; + Configurable cutAPL1{"cutAPL1", 0.69, "cutAPL par1"}; + Configurable cutAPL2{"cutAPL2", 0.5, "cutAPL par2"}; + /// Configurables gamma + Configurable gammaAsymmetryMax{"gammaAsymmetryMax", 0.95, "maximum photon asymetry."}; + Configurable gammaPsiPairMax{"gammaPsiPairMax", 0.1, "maximum psi angle of the track pair. "}; + Configurable gammaQtPtMultiplicator{"gammaQtPtMultiplicator", 0.11, "Multiply pt of V0s by this value to get the 2nd denominator in the armenteros cut. The products maximum value is fV0QtMax."}; + Configurable gammaQtMax{"gammaQtMax", 0.040, "the maximum value of the product, that is the maximum qt"}; + Configurable gammaV0RadiusMin{"gammaV0RadiusMin", 5., "the minimum V0 radius of the gamma decay"}; + Configurable gammaV0RadiusMax{"gammaV0RadiusMax", 180., "the maximum V0 radius of the gamma decay"}; + Configurable cutAlphaG1{"cutAlphaG1", 0.35, "maximum alpha gamma decay cut 1"}; + Configurable cutAlphaG2min{"cutAlphaG2min", 0.6, "minimum alpha gamma decay cut 2"}; + Configurable cutAlphaG2max{"cutAlphaG2max", 0.8, "maximum alpha gamma decay cut 2"}; + Configurable cutQTG{"cutQTG", 0.04, "maximum qt gamma decay"}; + + /// Kaon selection + template + bool selectionKaon(C const& collision, V0 const& v0) + { + // initialise dynamic variables + const float alpha = v0.alpha(); + const float qt = v0.qtarm(); + const float cosPA = v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); + const float q_K = {cutAPK0pinner * std::sqrt(std::abs(1 - ((alpha * alpha) / (cutAPK0alinner * cutAPK0alinner))))}; + /// Armenteros-Podolanski cut + if ((qt < cutQTK0min) || (qt > cutQTK0max) || (qt < q_K)) { + return false; + } + /// Cut on cosine pointing angle + if (cosPA < cutPAV0) { + return false; + } + /// Invariant mass cut + if (std::abs(v0.mK0Short() - 0.497611) > invariantMassCutK0Short) { + return false; + } + + return true; + }; + + /// Pion daughter selection + template + bool selectionPion(T const& track) + { + /// TOF cut daughter tracks + if (nSigmaTOFdautrack != 0. && std::abs(track.tofNSigmaPi()) > nSigmaTOFdautrack) { + return false; + } + /// Pion downsampling + if (downsamplingTsalisPions > 0. && downsampleTsalisCharged(track.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion]) == 0) { + return false; + } + return true; + } + + /// Lambda selection + template + bool selectionLambda(C const& collision, V0 const& v0) + { + // initialise dynamic variables + + const float alpha = v0.alpha(); + const float qt = v0.qtarm(); + const float cosPA = v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); + + const float q_L = cutAPL0 * std::sqrt(std::abs(1 - (((alpha - cutAPL1) * (alpha - cutAPL1)) / (cutAPL2 * cutAPL2)))); + /// Armenteros-Podolanski cut + if ((alpha < cutAlphaminL) || (alpha > cutAlphamaxL) || (qt < cutQTmin) || (qt > q_L)) { + return false; + } + /// Cut on cosine pointing angle + if (cosPA < cutPAV0) { + return false; + } + /// Invariant mass cut + if (std::abs(v0.mLambda() - 1.115683) > invariantMassCutLambda) { + return false; + } + return true; + }; + + /// Proton daughter selection + template + bool selectionProton(T const& track) + { + /// TOF cut daughter tracks + if (nSigmaTOFdautrack != 0. && std::abs(track.tofNSigmaPr()) > nSigmaTOFdautrack) { + return false; + } + /// Proton downsampling + if (downsamplingTsalisProtons > 0. && downsampleTsalisCharged(track.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton]) == 0) { + return false; + } + return true; + } + + /// Antilambda selection + template + bool selectionAntiLambda(C const& collision, V0 const& v0) + { + // initialise dynamic variables + const float alpha = v0.alpha(); + const float qt = v0.qtarm(); + const float cosPA = v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); + + const float q_L = cutAPL0 * std::sqrt(std::abs(1 - (((alpha + cutAPL1) * (alpha + cutAPL1)) / (cutAPL2 * cutAPL2)))); + /// Armenteros-Podolanski cut + if ((alpha < cutAlphaminAL) || (alpha > cutAlphamaxAL) || (qt < cutQTmin) || (qt > q_L)) { + return false; + } + /// Cut on cosine pointing angle + if (cosPA < cutPAV0) { + return false; + } + /// Invariant mass cut + if (std::abs(v0.mAntiLambda() - 1.115683) > invariantMassCutLambda) { + return false; + } + return true; + }; + + /// Gamma selection + template + bool selectionGamma(C const& collision, V0 const& v0) + { + // initialise dynamic variables + const float alpha = v0.alpha(); + const float qt = v0.qtarm(); + const float cosPA = v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); + const float pT = v0.pt(); + float lQtMaxPtDep = gammaQtPtMultiplicator * pT; + if (lQtMaxPtDep > gammaQtMax) { + lQtMaxPtDep = gammaQtMax; + } + + /// Armenteros-Podolanski cut + if ((((qt > gammaQtMax) || (std::abs(alpha) > cutAlphaG1)) || ((qt > cutQTG) || (std::abs(alpha) < cutAlphaG2min) || (std::abs(alpha) > cutAlphaG2max)))) { + return false; + } + if ((std::pow(alpha / gammaAsymmetryMax, 2) + std::pow(qt / lQtMaxPtDep, 2)) < 1) { + return false; + } + if (v0.v0radius() < gammaV0RadiusMin || v0.v0radius() > gammaV0RadiusMax) { + return false; + } + /// Cut on cosine pointing angle + if (cosPA < cutPAV0) { + return false; + } + if ((std::abs(v0.psipair()) < gammaPsiPairMax)) { + return false; + } + return true; + }; + + /// Electron daughter selection + template + bool selectionElectron(T const& track) + { + /// Electron downsampling + if (downsamplingTsalisElectrons > 0. && downsampleTsalisCharged(track.pt(), downsamplingTsalisElectrons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Electron]) == 0) { + return false; + } + return true; + } + + /// Funktion to fill skimmed tables + template + void fillSkimmedV0Table(V0 const& v0, T const& track, C const& collision, const float nSigmaTPC, const float nSigmaTOF, const float dEdxExp, const o2::track::PID::ID id, double dwnSmplFactor) + { + + const double ncl = track.tpcNClsFound(); + const double p = track.tpcInnerParam(); + const double mass = o2::track::pid_constants::sMasses[id]; + const double bg = p / mass; + const int multTPC = collision.multTPC(); + + const float alpha = v0.alpha(); + const float qt = v0.qtarm(); + const float cosPA = v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); + const float pT = v0.pt(); + const float v0radius = v0.v0radius(); + const float gammapsipair = v0.psipair(); + + const double pseudoRndm = track.pt() * 1000. - (long)(track.pt() * 1000); + if (pseudoRndm < dwnSmplFactor) { + rowTPCTree(track.tpcSignal(), + 1. / dEdxExp, + track.tpcInnerParam(), + track.tgl(), + track.signed1Pt(), + track.eta(), + track.phi(), + track.y(), + mass, + bg, + multTPC / 11000., + std::sqrt(nClNorm / ncl), + id, + nSigmaTPC, + nSigmaTOF, + alpha, + qt, + cosPA, + pT, + v0radius, + gammapsipair); + } + }; + + double tsalisCharged(double pt, double mass, double sqrts) + { + const double a = 6.81, b = 59.24; + const double c = 0.082, d = 0.151; + const double mt = std::sqrt(mass * mass + pt * pt); + const double n = a + b / sqrts; + const double T = c + d / sqrts; + const double p0 = n * T; + const double result = std::pow((1. + mt / p0), -n); + return result; + }; + + /// Random downsampling trigger function using Tsalis/Hagedorn spectra fit (sqrt(s) = 62.4 GeV to 13 TeV) + /// as in https://iopscience.iop.org/article/10.1088/2399-6528/aab00f/pdf + TRandom3* fRndm = new TRandom3(0); + int downsampleTsalisCharged(double pt, double factor1Pt, double sqrts, double mass) + { + const double prob = tsalisCharged(pt, mass, sqrts) * pt; + const double probNorm = tsalisCharged(1., mass, sqrts); + int triggerMask = 0; + if ((fRndm->Rndm() * ((prob / probNorm) * pt * pt)) < factor1Pt) { + triggerMask = 1; + } + + return triggerMask; + }; + + /// Event selection + template + bool isEventSelected(const CollisionType& collision, const TrackType& tracks) + { + if (applyEvSel == 1) { + if (!collision.sel7()) { + return false; + } + } else if (applyEvSel == 2) { + if (!collision.sel8()) { + return false; + } + } + return true; + }; + + void init(o2::framework::InitContext& initContext) + { + } + + void process(Coll::iterator const& collision, Trks const& tracks, aod::V0Datas const& v0s) + { + /// Check event slection + if (!isEventSelected(collision, tracks)) { + return; + } + + rowTPCTree.reserve(tracks.size()); + + /// Loop over v0 candidates + for (auto v0 : v0s) { + auto posTrack = v0.posTrack_as(); + auto negTrack = v0.negTrack_as(); + /// Fill table for Kaons + if (selectionKaon(collision, v0)) { + if (selectionPion(posTrack)) { + fillSkimmedV0Table(v0, posTrack, collision, posTrack.tpcNSigmaPi(), posTrack.tofNSigmaPi(), posTrack.tpcExpSignalPi(posTrack.tpcSignal()), o2::track::PID::Pion, dwnSmplFactor_Pi); + } + if (selectionPion(negTrack)) { + fillSkimmedV0Table(v0, negTrack, collision, negTrack.tpcNSigmaPi(), negTrack.tofNSigmaPi(), negTrack.tpcExpSignalPi(negTrack.tpcSignal()), o2::track::PID::Pion, dwnSmplFactor_Pi); + } + } + /// Fill table for Lambdas + if (selectionLambda(collision, v0)) { + if (selectionProton(posTrack)) { + fillSkimmedV0Table(v0, posTrack, collision, posTrack.tpcNSigmaPr(), posTrack.tofNSigmaPr(), posTrack.tpcExpSignalPr(posTrack.tpcSignal()), o2::track::PID::Proton, dwnSmplFactor_Pr); + } + if (selectionPion(negTrack)) { + fillSkimmedV0Table(v0, negTrack, collision, negTrack.tpcNSigmaPi(), negTrack.tofNSigmaPi(), negTrack.tpcExpSignalPi(negTrack.tpcSignal()), o2::track::PID::Pion, dwnSmplFactor_Pi); + } + } + /// Fill table for Antilambdas + if (selectionAntiLambda(collision, v0)) { + if (selectionPion(posTrack)) { + fillSkimmedV0Table(v0, posTrack, collision, posTrack.tpcNSigmaPi(), posTrack.tofNSigmaPi(), posTrack.tpcExpSignalPi(posTrack.tpcSignal()), o2::track::PID::Pion, dwnSmplFactor_Pi); + } + if (selectionProton(negTrack)) { + fillSkimmedV0Table(v0, negTrack, collision, negTrack.tpcNSigmaPr(), negTrack.tofNSigmaPr(), negTrack.tpcExpSignalPr(negTrack.tpcSignal()), o2::track::PID::Proton, dwnSmplFactor_Pr); + } + } + /// Fill table for Gammas + if (selectionGamma(collision, v0)) { + if (selectionElectron(posTrack)) { + fillSkimmedV0Table(v0, posTrack, collision, posTrack.tpcNSigmaEl(), posTrack.tofNSigmaEl(), posTrack.tpcExpSignalEl(posTrack.tpcSignal()), o2::track::PID::Electron, dwnSmplFactor_El); + } + if (selectionElectron(negTrack)) { + fillSkimmedV0Table(v0, negTrack, collision, negTrack.tpcNSigmaEl(), negTrack.tofNSigmaEl(), negTrack.tpcExpSignalEl(negTrack.tpcSignal()), o2::track::PID::Electron, dwnSmplFactor_El); + } + } + } /// Loop V0 candidates + } /// process +}; /// struct TreeWriterTpcV0 + +struct TreeWriterTPCTOF { + using Trks = soa::Join; + using Coll = soa::Join; + + /// Tables to be produced + Produces rowTPCTOFTree; + + /// Configurables + Configurable nClNorm{"nClNorm", 152., "Number of cluster normalization. Run 2: 159, Run 3 152"}; + Configurable applyEvSel{"applyEvSel", 2, "Flag to apply rapidity cut: 0 -> no event selection, 1 -> Run 2 event selection, 2 -> Run 3 event selection"}; + Configurable applyTrkSel{"applyTrkSel", 1, "Flag to apply track selection: 0 -> no track selection, 1 -> track selection"}; + /// Proton + Configurable maxMomTPCOnlyPr{"maxMomTPCOnlyPr", 0.6, "Maximum momentum for TPC only cut proton"}; + Configurable nSigmaTPCOnlyPr{"nSigmaTPCOnlyPr", 4., "number of sigma for TPC only cut proton"}; + Configurable nSigmaTPC_TPCTOF_Pr{"nSigmaTPC_TPCTOF_Pr", 4., "number of sigma for TPC cut for TPC and TOF combined proton"}; + Configurable nSigmaTOF_TPCTOF_Pr{"nSigmaTOF_TPCTOF_Pr", 3., "number of sigma for TOF cut for TPC and TOF combined proton"}; + Configurable dwnSmplFactor_Pr{"dwnSmplFactor_Pr", 1., "downsampling factor for protons, default fraction to keep is 1."}; + /// Kaon + Configurable maxMomTPCOnlyKa{"maxMomTPCOnlyKa", 0.3, "Maximum momentum for TPC only cut kaon"}; + Configurable nSigmaTPCOnlyKa{"nSigmaTPCOnlyKa", 4., "number of sigma for TPC only cut kaon"}; + Configurable nSigmaTPC_TPCTOF_Ka{"nSigmaTPC_TPCTOF_Ka", 4., "number of sigma for TPC cut for TPC and TOF combined kaon"}; + Configurable nSigmaTOF_TPCTOF_Ka{"nSigmaTOF_TPCTOF_Ka", 3., "number of sigma for TOF cut for TPC and TOF combined kaon"}; + Configurable dwnSmplFactor_Ka{"dwnSmplFactor_Ka", 1., "downsampling factor for kaons, default fraction to keep is 1."}; + /// Pion + Configurable maxMomTPCOnlyPi{"maxMomTPCOnlyPi", 0.5, "Maximum momentum for TPC only cut pion"}; + Configurable nSigmaTPCOnlyPi{"nSigmaTPCOnlyPi", 4., "number of sigma for TPC only cut pion"}; + Configurable nSigmaTPC_TPCTOF_Pi{"nSigmaTPC_TPCTOF_Pi", 4., "number of sigma for TPC cut for TPC and TOF combined pion"}; + Configurable nSigmaTOF_TPCTOF_Pi{"nSigmaTOF_TPCTOF_Pi", 4., "number of sigma for TOF cut for TPC and TOF combined pion"}; + Configurable dwnSmplFactor_Pi{"dwnSmplFactor_Pi", 1., "downsampling factor for pions, default fraction to keep is 1."}; + /// pT dependent downsampling + Configurable sqrtSNN{"sqrt_s_NN", 0., "sqrt(s_NN), used for downsampling with the Tsallis distribution"}; + Configurable downsamplingTsalisProtons{"downsamplingTsalisProtons", -1., "Downsampling factor to reduce the number of protons"}; + Configurable downsamplingTsalisKaons{"downsamplingTsalisKaons", -1., "Downsampling factor to reduce the number of kaons"}; + Configurable downsamplingTsalisPions{"downsamplingTsalisPions", -1., "Downsampling factor to reduce the number of pions"}; + + double tsalisCharged(double pt, double mass, double sqrts) + { + const double a = 6.81, b = 59.24; + const double c = 0.082, d = 0.151; + double mt = std::sqrt(mass * mass + pt * pt); + double n = a + b / sqrts; + double T = c + d / sqrts; + double p0 = n * T; + double result = pow((1. + mt / p0), -n); + return result; + }; + + /// Random downsampling trigger function using Tsalis/Hagedorn spectra fit (sqrt(s) = 62.4 GeV to 13 TeV) + /// as in https://iopscience.iop.org/article/10.1088/2399-6528/aab00f/pdf + TRandom3* fRndm = new TRandom3(0); + bool downsampleTsalisCharged(double pt, float factor1Pt, double sqrts, double mass) + { + if (factor1Pt < 0.) { + return true; + } + const double prob = tsalisCharged(pt, mass, sqrts) * pt; + const double probNorm = tsalisCharged(1., mass, sqrts); + if ((fRndm->Rndm() * ((prob / probNorm) * pt * pt)) > factor1Pt) { + return false; + } else { + return true; + } + }; + + /// Function to fill trees + template + void fillSkimmedTPCTOFTable(T const& track, C const& collision, const float nSigmaTPC, const float nSigmaTOF, const float dEdxExp, const o2::track::PID::ID id, double dwnSmplFactor) + { + + const double ncl = track.tpcNClsFound(); + const double p = track.tpcInnerParam(); + const double mass = o2::track::pid_constants::sMasses[id]; + const double bg = p / mass; + const int multTPC = collision.multTPC(); + + const double pseudoRndm = track.pt() * 1000. - (long)(track.pt() * 1000); + if (pseudoRndm < dwnSmplFactor) { + rowTPCTOFTree(track.tpcSignal(), + 1. / dEdxExp, + track.tpcInnerParam(), + track.tgl(), + track.signed1Pt(), + track.eta(), + track.phi(), + track.y(), + mass, + bg, + multTPC / 11000., + std::sqrt(nClNorm / ncl), + id, + nSigmaTPC, + nSigmaTOF); + } + }; + + /// Event selection + template + bool isEventSelected(const CollisionType& collision, const TrackType& tracks) + { + if (applyEvSel == 1) { + if (!collision.sel7()) { + return false; + } + } else if (applyEvSel == 2) { + if (!collision.sel8()) { + return false; + } + } + return true; + }; + + /// Track selection + template + bool isTrackSelected(const CollisionType& collision, const TrackType& track) + { + if (!track.isGlobalTrack()) { // Skipping non global tracks + return false; + } + if (!track.hasITS()) { // Skipping tracks without ITS + return false; + } + if (!track.hasTPC()) { // Skipping tracks without TPC + return false; + } + return true; + }; + + void init(o2::framework::InitContext& initContext) + { + } + void process(Coll::iterator const& collision, Trks const& tracks) + { + /// Check event selection + if (!isEventSelected(collision, tracks)) { + return; + } + rowTPCTOFTree.reserve(tracks.size()); + for (auto const& trk : tracks) { + + /// Check track selection + if (applyTrkSel == 1 && !isTrackSelected(collision, trk)) { + continue; + } + /// Fill tree for protons + if (trk.tpcInnerParam() < maxMomTPCOnlyPr && std::abs(trk.tpcNSigmaPr()) < nSigmaTPCOnlyPr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton])) { + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaPr(), trk.tofNSigmaPr(), trk.tpcExpSignalPr(trk.tpcSignal()), o2::track::PID::Proton, dwnSmplFactor_Pr); + } else if (trk.tpcInnerParam() > maxMomTPCOnlyPr && std::abs(trk.tofNSigmaPr()) < nSigmaTOF_TPCTOF_Pr && std::abs(trk.tpcNSigmaPr()) < nSigmaTPC_TPCTOF_Pr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton])) { + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaPr(), trk.tofNSigmaPr(), trk.tpcExpSignalPr(trk.tpcSignal()), o2::track::PID::Proton, dwnSmplFactor_Pr); + } + /// Fill tree for kaons + if (trk.tpcInnerParam() < maxMomTPCOnlyKa && std::abs(trk.tpcNSigmaKa()) < nSigmaTPCOnlyKa && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon])) { + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaKa(), trk.tofNSigmaKa(), trk.tpcExpSignalKa(trk.tpcSignal()), o2::track::PID::Kaon, dwnSmplFactor_Ka); + } else if (trk.tpcInnerParam() > maxMomTPCOnlyKa && std::abs(trk.tofNSigmaKa()) < nSigmaTOF_TPCTOF_Ka && std::abs(trk.tpcNSigmaKa()) < nSigmaTPC_TPCTOF_Ka && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon])) { + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaKa(), trk.tofNSigmaKa(), trk.tpcExpSignalKa(trk.tpcSignal()), o2::track::PID::Kaon, dwnSmplFactor_Ka); + } + /// Fill tree pions + if (trk.tpcInnerParam() < maxMomTPCOnlyPi && std::abs(trk.tpcNSigmaPi()) < nSigmaTPCOnlyPi && downsampleTsalisCharged(trk.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion])) { + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaPi(), trk.tofNSigmaPi(), trk.tpcExpSignalPi(trk.tpcSignal()), o2::track::PID::Pion, dwnSmplFactor_Pi); + } else if (trk.tpcInnerParam() > maxMomTPCOnlyPi && std::abs(trk.tofNSigmaPi()) < nSigmaTOF_TPCTOF_Pi && std::abs(trk.tpcNSigmaPi()) < nSigmaTPC_TPCTOF_Pi && downsampleTsalisCharged(trk.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion])) { + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaPi(), trk.tofNSigmaPi(), trk.tpcExpSignalPi(trk.tpcSignal()), o2::track::PID::Pion, dwnSmplFactor_Pi); + } + } /// Loop tracks + } /// process +}; /// struct TreeWriterTPCTOF + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + auto workflow = WorkflowSpec{adaptAnalysisTask(cfgc)}; + workflow.push_back(adaptAnalysisTask(cfgc)); + return workflow; +} \ No newline at end of file diff --git a/DPG/Tasks/TPC/tpcSkimsTableCreator.h b/DPG/Tasks/TPC/tpcSkimsTableCreator.h new file mode 100644 index 00000000000..c06efc4f6e6 --- /dev/null +++ b/DPG/Tasks/TPC/tpcSkimsTableCreator.h @@ -0,0 +1,81 @@ +// 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 tpcSkimsTableCreator.h +/// \author Annalena Kalteyer +/// \author Christian Sonnabend +/// \author Jeremy Wilkinson + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::track; +using namespace o2::dataformats; + +namespace o2::aod +{ +namespace tpcskims +{ +DECLARE_SOA_COLUMN(InvDeDxExpTPC, invdEdxExpTPC, float); +DECLARE_SOA_COLUMN(Mass, mass, float); +DECLARE_SOA_COLUMN(BetaGamma, bg, float); +DECLARE_SOA_COLUMN(NormMultTPC, normMultTPC, float); +DECLARE_SOA_COLUMN(NormNClustersTPC, normNClustersTPC, float); +DECLARE_SOA_COLUMN(PidIndex, pidIndexTPC, uint8_t); +DECLARE_SOA_COLUMN(NSigTPC, nsigTPC, float); +DECLARE_SOA_COLUMN(NSigTOF, nsigTOF, float); +DECLARE_SOA_COLUMN(AlphaV0, alphaV0, float); +DECLARE_SOA_COLUMN(QtV0, qtV0, float); +DECLARE_SOA_COLUMN(CosPAV0, cosPAV0, float); +DECLARE_SOA_COLUMN(PtV0, ptV0, float); +DECLARE_SOA_COLUMN(RadiusV0, radiusV0, float); +DECLARE_SOA_COLUMN(GammaPsiPair, gammaPsiPair, float); +} // namespace tpcskims +DECLARE_SOA_TABLE(SkimmedTPCV0Tree, "AOD", "TPCSKIMV0TREE", + o2::aod::track::TPCSignal, + tpcskims::InvDeDxExpTPC, + o2::aod::track::TPCInnerParam, + o2::aod::track::Tgl, + o2::aod::track::Signed1Pt, + o2::aod::track::Eta, + o2::aod::track::Phi, + o2::aod::track::Y, + tpcskims::Mass, + tpcskims::BetaGamma, + tpcskims::NormMultTPC, + tpcskims::NormNClustersTPC, + tpcskims::PidIndex, + tpcskims::NSigTPC, + tpcskims::NSigTOF, + tpcskims::AlphaV0, + tpcskims::QtV0, + tpcskims::CosPAV0, + tpcskims::PtV0, + tpcskims::RadiusV0, + tpcskims::GammaPsiPair); + +DECLARE_SOA_TABLE(SkimmedTPCTOFTree, "AOD", "TPCTOFSKIMTREE", + o2::aod::track::TPCSignal, + tpcskims::InvDeDxExpTPC, + o2::aod::track::TPCInnerParam, + o2::aod::track::Tgl, + o2::aod::track::Signed1Pt, + o2::aod::track::Eta, + o2::aod::track::Phi, + o2::aod::track::Y, + tpcskims::Mass, + tpcskims::BetaGamma, + tpcskims::NormMultTPC, + tpcskims::NormNClustersTPC, + tpcskims::PidIndex, + tpcskims::NSigTPC, + tpcskims::NSigTOF); +} // namespace o2::aod \ No newline at end of file