From c5e563d49c0181788808593d47c484a0efa5f9bb Mon Sep 17 00:00:00 2001 From: AnnalenaKa Date: Fri, 8 Jul 2022 09:56:55 +0200 Subject: [PATCH 01/11] Add downsampling Tsallis --- Common/DataModel/PIDResponse.h | 28 +- Common/TableProducer/PID/pidTPCFull.cxx | 12 +- Common/Tools/CMakeLists.txt | 6 +- Common/Tools/tpcSkimsTableCreator.cxx | 538 ++++++++++++++++++++++++ 4 files changed, 569 insertions(+), 15 deletions(-) create mode 100644 Common/Tools/tpcSkimsTableCreator.cxx diff --git a/Common/DataModel/PIDResponse.h b/Common/DataModel/PIDResponse.h index 511af101cd5..92f344dfe56 100644 --- a/Common/DataModel/PIDResponse.h +++ b/Common/DataModel/PIDResponse.h @@ -888,6 +888,16 @@ DECLARE_SOA_TABLE(pidTOFAl, "AOD", "pidTOFAl", //! Table of the TOF response wit namespace pidtpc { // Expected signals +DECLARE_SOA_COLUMN(TPCExpSignalEl, tpcExpSignalEl, float); //! Expected resolution with the TPC detector for electron +DECLARE_SOA_COLUMN(TPCExpSignalMu, tpcExpSignalMu, float); //! Expected resolution with the TPC detector for muon +DECLARE_SOA_COLUMN(TPCExpSignalPi, tpcExpSignalPi, float); //! Expected resolution with the TPC detector for pion +DECLARE_SOA_COLUMN(TPCExpSignalKa, tpcExpSignalKa, float); //! Expected resolution with the TPC detector for kaon +DECLARE_SOA_COLUMN(TPCExpSignalPr, tpcExpSignalPr, float); //! Expected resolution with the TPC detector for proton +DECLARE_SOA_COLUMN(TPCExpSignalDe, tpcExpSignalDe, float); //! Expected resolution with the TPC detector for deuteron +DECLARE_SOA_COLUMN(TPCExpSignalTr, tpcExpSignalTr, float); //! Expected resolution with the TPC detector for triton +DECLARE_SOA_COLUMN(TPCExpSignalHe, tpcExpSignalHe, float); //! Expected resolution with the TPC detector for helium3 +DECLARE_SOA_COLUMN(TPCExpSignalAl, tpcExpSignalAl, float); //! Expected resolution with the TPC detector for alpha +// 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 +978,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/pidTPCFull.cxx b/Common/TableProducer/PID/pidTPCFull.cxx index 1131259c07b..769e2e30f2e 100644 --- a/Common/TableProducer/PID/pidTPCFull.cxx +++ b/Common/TableProducer/PID/pidTPCFull.cxx @@ -275,24 +275,26 @@ struct tpcPidFull { // Here comes the application of the network. The output--dimensions of the network dtermine the application: 1: mean, 2: sigma, 3: sigma asymmetric // For now only the option 2: sigma will be used. The other options are kept if there would be demand later on if (network.getOutputDimensions() == 1) { - table(response.GetExpectedSigma(collisions.iteratorAt(trk.collisionId()), trk, pid), + table(response.GetExpectedSignal(trk, pid), + response.GetExpectedSigma(collisions.iteratorAt(trk.collisionId()), trk, pid), (trk.tpcSignal() - network_prediction[count_tracks + tracks_size * pid] * response.GetExpectedSignal(trk, pid)) / response.GetExpectedSigma(collisions.iteratorAt(trk.collisionId()), trk, pid)); } else if (network.getOutputDimensions() == 2) { - table((network_prediction[2 * (count_tracks + tracks_size * pid) + 1] - network_prediction[2 * (count_tracks + tracks_size * pid)]) * response.GetExpectedSignal(trk, pid), + table(response.GetExpectedSignal(trk, pid), (network_prediction[2 * (count_tracks + tracks_size * pid) + 1] - network_prediction[2 * (count_tracks + tracks_size * pid)]) * response.GetExpectedSignal(trk, pid), (trk.tpcSignal() / response.GetExpectedSignal(trk, pid) - network_prediction[2 * (count_tracks + tracks_size * pid)]) / (network_prediction[2 * (count_tracks + tracks_size * pid) + 1] - network_prediction[2 * (count_tracks + tracks_size * pid)])); } else if (network.getOutputDimensions() == 3) { if (trk.tpcSignal() / response.GetExpectedSignal(trk, pid) >= network_prediction[3 * (count_tracks + tracks_size * pid)]) { - table((network_prediction[3 * (count_tracks + tracks_size * pid) + 1] - network_prediction[3 * (count_tracks + tracks_size * pid)]) * response.GetExpectedSignal(trk, pid), + table(response.GetExpectedSignal(trk, pid), (network_prediction[3 * (count_tracks + tracks_size * pid) + 1] - network_prediction[3 * (count_tracks + tracks_size * pid)]) * response.GetExpectedSignal(trk, pid), (trk.tpcSignal() / response.GetExpectedSignal(trk, pid) - network_prediction[3 * (count_tracks + tracks_size * pid)]) / (network_prediction[3 * (count_tracks + tracks_size * pid) + 1] - network_prediction[3 * (count_tracks + tracks_size * pid)])); } else { - table((network_prediction[3 * (count_tracks + tracks_size * pid)] - network_prediction[3 * (count_tracks + tracks_size * pid) + 2]) * response.GetExpectedSignal(trk, pid), + table(response.GetExpectedSignal(trk, pid),(network_prediction[3 * (count_tracks + tracks_size * pid)] - network_prediction[3 * (count_tracks + tracks_size * pid) + 2]) * response.GetExpectedSignal(trk, pid), (trk.tpcSignal() / response.GetExpectedSignal(trk, pid) - network_prediction[3 * (count_tracks + tracks_size * pid)]) / (network_prediction[3 * (count_tracks + tracks_size * pid)] - network_prediction[3 * (count_tracks + tracks_size * pid) + 2])); } } else { LOGF(fatal, "Network output-dimensions incompatible!"); } } else { - table(response.GetExpectedSigma(collisions.iteratorAt(trk.collisionId()), trk, pid), + table(response.GetExpectedSignal(trk, pid), + response.GetExpectedSigma(collisions.iteratorAt(trk.collisionId()), trk, pid), response.GetNumberOfSigma(collisions.iteratorAt(trk.collisionId()), trk, pid)); } }; diff --git a/Common/Tools/CMakeLists.txt b/Common/Tools/CMakeLists.txt index 3e8898ca1a4..fb3d4fff106 100644 --- a/Common/Tools/CMakeLists.txt +++ b/Common/Tools/CMakeLists.txt @@ -24,4 +24,8 @@ 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 + ) +o2physics_add_dpl_workflow(pid-tpc-skimscreation + SOURCES tpcSkimsTableCreator.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/Common/Tools/tpcSkimsTableCreator.cxx b/Common/Tools/tpcSkimsTableCreator.cxx new file mode 100644 index 00000000000..60c8d2f70a5 --- /dev/null +++ b/Common/Tools/tpcSkimsTableCreator.cxx @@ -0,0 +1,538 @@ +// 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" +/// 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; + +// void customize(std::vector& workflowOptions) +// { +// std::vector options{// runtime customisation goes here +// {"useV0", VariantType::Int, 0, {"Use V0 information for QA"}}}; +// std::swap(workflowOptions, options); +// } + +namespace o2::aod +{ +namespace tpcskims +{ +DECLARE_SOA_COLUMN(InvDeDxExpV0, invdEdxExpV0, float); +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); +} // namespace tpcskims +DECLARE_SOA_TABLE(SkimmedTPCV0Tree, "AOD", "TPCSKIMV0TREE", + o2::aod::track::TPCSignal, + tpcskims::InvDeDxExpV0, + 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); + +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 + +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."}; + /// 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 + float alpha = v0.alpha(); + float qt = v0.qtarm(); + float cosPA = v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); + float q_K = {cutAPK0pinner * sqrt(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 (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. && 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 + + float alpha = v0.alpha(); + float qt = v0.qtarm(); + float cosPA = v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); + + float q_L = cutAPL0 * sqrt(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 (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. && 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 + float alpha = v0.alpha(); + float qt = v0.qtarm(); + float cosPA = v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); + + float q_L = cutAPL0 * sqrt(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 (abs(v0.mAntiLambda() - 1.115683) > invariantMassCutLambda) { + return false; + } + return true; + }; + + /// Gamma selection + template + bool selectionGamma(C const& collision, V0 const& v0) + { + // initialise dynamic variables + float alpha = v0.alpha(); + float qt = v0.qtarm(); + float cosPA = v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); + float pT = v0.pt(); + float lQtMaxPtDep = gammaQtPtMultiplicator * pT; + if (lQtMaxPtDep > gammaQtMax) { + lQtMaxPtDep = gammaQtMax; + } + + /// Armenteros-Podolanski cut + if ((((qt > gammaQtMax) || (TMath::Abs(alpha) > cutAlphaG1)) || ((qt > cutQTG) || (TMath::Abs(alpha) < cutAlphaG2min) || (TMath::Abs(alpha) > cutAlphaG2max)))) { + return false; + } + if ((TMath::Power(alpha / gammaAsymmetryMax, 2) + TMath::Power(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 ((TMath::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; + } + + template + void fillSkimmedV0Table(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(); + + 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(159. / ncl), + id, + nSigmaTPC, + nSigmaTOF); + } + }; + + + + 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=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); + int downsampleTsalisCharged(double pt, double factor1Pt, double sqrts, double mass){ + double prob=tsalisCharged(pt,mass,sqrts)*pt; + double probNorm=tsalisCharged(1.,mass,sqrts); + int triggerMask=0; + if ((fRndm->Rndm()*((prob/probNorm)*pt*pt))(); + auto negTrack = v0.negTrack_as(); + /// Fill table for Kaons + if (selectionKaon(collision, v0, posTrack, negTrack)) { + if (selectionPion(posTrack)) + { + fillSkimmedV0Table(posTrack, collision, posTrack.tpcNSigmaPi(), posTrack.tofNSigmaPi(), posTrack.tpcExpSignalPi(), o2::track::PID::Pion, dwnSmplFactor_Pi); + } + if (selectionPion(negTrack)) + { + fillSkimmedV0Table(negTrack, collision, negTrack.tpcNSigmaPi(), negTrack.tofNSigmaPi(), negTrack.tpcExpSignalPi(), o2::track::PID::Pion, dwnSmplFactor_Pi); + } + } + /// Fill table for Lambdas + if (selectionLambda(collision, v0, posTrack, negTrack)) { + if (selectionProton(posTrack)) + { + fillSkimmedV0Table(posTrack, collision, posTrack.tpcNSigmaPr(), posTrack.tofNSigmaPr(), posTrack.tpcExpSignalPr(), o2::track::PID::Proton, dwnSmplFactor_Pr); + } + if (selectionPion(negTrack)) + { + fillSkimmedV0Table(negTrack, collision, negTrack.tpcNSigmaPi(), negTrack.tofNSigmaPi(), negTrack.tpcExpSignalPi(), o2::track::PID::Pion, dwnSmplFactor_Pi); + } + } + /// Fill table for Antilambdas + if (selectionAntiLambda(collision, v0, posTrack, negTrack)) { + if (selectionPion(posTrack)) + { + fillSkimmedV0Table(posTrack, collision, posTrack.tpcNSigmaPi(), posTrack.tofNSigmaPi(), posTrack.tpcExpSignalPi(), o2::track::PID::Pion, dwnSmplFactor_Pi); + } + if (selectionProton(negTrack)) + { + fillSkimmedV0Table(negTrack, collision, negTrack.tpcNSigmaPr(), negTrack.tofNSigmaPr(), negTrack.tpcExpSignalPr(), o2::track::PID::Proton, dwnSmplFactor_Pr); + } + } + /// Fill table for Gammas + if (selectionGamma(collision, v0)) { + if (selectionElectron(posTrack)) + { + fillSkimmedV0Table(posTrack, collision, posTrack.tpcNSigmaEl(), posTrack.tofNSigmaEl(), posTrack.tpcExpSignalEl(), o2::track::PID::Electron, dwnSmplFactor_El); + } + if (selectionElectron(negTrack)) + { + fillSkimmedV0Table(negTrack, collision, negTrack.tpcNSigmaEl(), negTrack.tofNSigmaEl(), negTrack.tpcExpSignalEl(), 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 + /// 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{"downsamplingTsalisKaonsn", -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=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, double factor1Pt, double sqrts, double mass, float dwnsmplFact){ + if (dwnsmplFact<0.) return true; + double prob=tsalisCharged(pt,mass,sqrts)*pt; + double probNorm=tsalisCharged(1.,mass,sqrts); + if ((fRndm->Rndm()*((prob/probNorm)*pt*pt))>factor1Pt) + { + return false; + } + else { + return true; + } + + }; + + 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(); + + 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(159. / ncl), + id, + nSigmaTPC, + nSigmaTOF); + } + }; + + void init(o2::framework::InitContext& initContext) + { + } + void process(Coll::iterator const& collision, Trks const& tracks) + { + rowTPCTOFTree.reserve(tracks.size()); + for (auto const& trk : tracks) { + /// Fill tree for protons + if (trk.tpcInnerParam() < maxMomTPCOnlyPr && 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(), o2::track::PID::Proton, dwnSmplFactor_Pr); + } else if (trk.tpcInnerParam() > maxMomTPCOnlyPr && trk.tofNSigmaPr() < nSigmaTOF_TPCTOF_Pr && 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(), o2::track::PID::Proton, dwnSmplFactor_Pr); + } + /// Fill tree for kaons + if (trk.tpcInnerParam() < maxMomTPCOnlyKa && 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(), o2::track::PID::Kaon, dwnSmplFactor_Ka); + } else if (trk.tpcInnerParam() > maxMomTPCOnlyKa && trk.tofNSigmaKa() < nSigmaTOF_TPCTOF_Ka && 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(), o2::track::PID::Kaon, dwnSmplFactor_Ka); + } + /// Fill tree pions + if (trk.tpcInnerParam() < maxMomTPCOnlyPi && 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(), o2::track::PID::Pion, dwnSmplFactor_Pi); + } else if (trk.tpcInnerParam() > maxMomTPCOnlyPi && trk.tofNSigmaPi() < nSigmaTOF_TPCTOF_Pi && 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(), o2::track::PID::Pion, dwnSmplFactor_Pi); + } + } /// Loop tracks + } /// process +}; /// struct TreeWriterTPCTOF + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + // const int useV0 = cfgc.options().get("useV0"); + auto workflow = WorkflowSpec{adaptAnalysisTask(cfgc)}; + workflow.push_back(adaptAnalysisTask(cfgc)); + // if (useV0) { + + // } + return workflow; +} \ No newline at end of file From 563b159ceae9350f91d58c1c74309a1800c40dbd Mon Sep 17 00:00:00 2001 From: AnnalenaKa Date: Fri, 8 Jul 2022 10:08:36 +0200 Subject: [PATCH 02/11] fix errors --- Common/Tools/tpcSkimsTableCreator.cxx | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/Common/Tools/tpcSkimsTableCreator.cxx b/Common/Tools/tpcSkimsTableCreator.cxx index 60c8d2f70a5..b668aa855ca 100644 --- a/Common/Tools/tpcSkimsTableCreator.cxx +++ b/Common/Tools/tpcSkimsTableCreator.cxx @@ -360,7 +360,7 @@ struct TreeWriterTpcV0 { auto posTrack = v0.posTrack_as(); auto negTrack = v0.negTrack_as(); /// Fill table for Kaons - if (selectionKaon(collision, v0, posTrack, negTrack)) { + if (selectionKaon(collision, v0)) { if (selectionPion(posTrack)) { fillSkimmedV0Table(posTrack, collision, posTrack.tpcNSigmaPi(), posTrack.tofNSigmaPi(), posTrack.tpcExpSignalPi(), o2::track::PID::Pion, dwnSmplFactor_Pi); @@ -371,7 +371,7 @@ struct TreeWriterTpcV0 { } } /// Fill table for Lambdas - if (selectionLambda(collision, v0, posTrack, negTrack)) { + if (selectionLambda(collision, v0)) { if (selectionProton(posTrack)) { fillSkimmedV0Table(posTrack, collision, posTrack.tpcNSigmaPr(), posTrack.tofNSigmaPr(), posTrack.tpcExpSignalPr(), o2::track::PID::Proton, dwnSmplFactor_Pr); @@ -382,7 +382,7 @@ struct TreeWriterTpcV0 { } } /// Fill table for Antilambdas - if (selectionAntiLambda(collision, v0, posTrack, negTrack)) { + if (selectionAntiLambda(collision, v0)) { if (selectionPion(posTrack)) { fillSkimmedV0Table(posTrack, collision, posTrack.tpcNSigmaPi(), posTrack.tofNSigmaPi(), posTrack.tpcExpSignalPi(), o2::track::PID::Pion, dwnSmplFactor_Pi); @@ -505,21 +505,21 @@ struct TreeWriterTPCTOF { rowTPCTOFTree.reserve(tracks.size()); for (auto const& trk : tracks) { /// Fill tree for protons - if (trk.tpcInnerParam() < maxMomTPCOnlyPr && trk.tpcNSigmaPr() < nSigmaTPCOnlyPr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton])) { + if (trk.tpcInnerParam() < maxMomTPCOnlyPr && trk.tpcNSigmaPr() < nSigmaTPCOnlyPr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton], downsamplingTsalisProtons)) { fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaPr(), trk.tofNSigmaPr(), trk.tpcExpSignalPr(), o2::track::PID::Proton, dwnSmplFactor_Pr); - } else if (trk.tpcInnerParam() > maxMomTPCOnlyPr && trk.tofNSigmaPr() < nSigmaTOF_TPCTOF_Pr && trk.tpcNSigmaPr() < nSigmaTPC_TPCTOF_Pr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton])) { + } else if (trk.tpcInnerParam() > maxMomTPCOnlyPr && trk.tofNSigmaPr() < nSigmaTOF_TPCTOF_Pr && trk.tpcNSigmaPr() < nSigmaTPC_TPCTOF_Pr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton], downsamplingTsalisProtons)) { fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaPr(), trk.tofNSigmaPr(), trk.tpcExpSignalPr(), o2::track::PID::Proton, dwnSmplFactor_Pr); } /// Fill tree for kaons - if (trk.tpcInnerParam() < maxMomTPCOnlyKa && trk.tpcNSigmaKa() < nSigmaTPCOnlyKa && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon])) { + if (trk.tpcInnerParam() < maxMomTPCOnlyKa && trk.tpcNSigmaKa() < nSigmaTPCOnlyKa && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon], downsamplingTsalisKaons)) { fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaKa(), trk.tofNSigmaKa(), trk.tpcExpSignalKa(), o2::track::PID::Kaon, dwnSmplFactor_Ka); - } else if (trk.tpcInnerParam() > maxMomTPCOnlyKa && trk.tofNSigmaKa() < nSigmaTOF_TPCTOF_Ka && trk.tpcNSigmaKa() < nSigmaTPC_TPCTOF_Ka && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon])) { + } else if (trk.tpcInnerParam() > maxMomTPCOnlyKa && trk.tofNSigmaKa() < nSigmaTOF_TPCTOF_Ka && trk.tpcNSigmaKa() < nSigmaTPC_TPCTOF_Ka && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon], downsamplingTsalisKaons)) { fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaKa(), trk.tofNSigmaKa(), trk.tpcExpSignalKa(), o2::track::PID::Kaon, dwnSmplFactor_Ka); } /// Fill tree pions - if (trk.tpcInnerParam() < maxMomTPCOnlyPi && trk.tpcNSigmaPi() < nSigmaTPCOnlyPi && downsampleTsalisCharged(trk.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion])) { + if (trk.tpcInnerParam() < maxMomTPCOnlyPi && trk.tpcNSigmaPi() < nSigmaTPCOnlyPi && downsampleTsalisCharged(trk.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion], downsamplingTsalisPions)) { fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaPi(), trk.tofNSigmaPi(), trk.tpcExpSignalPi(), o2::track::PID::Pion, dwnSmplFactor_Pi); - } else if (trk.tpcInnerParam() > maxMomTPCOnlyPi && trk.tofNSigmaPi() < nSigmaTOF_TPCTOF_Pi && trk.tpcNSigmaPi() < nSigmaTPC_TPCTOF_Pi && downsampleTsalisCharged(trk.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion])) { + } else if (trk.tpcInnerParam() > maxMomTPCOnlyPi && trk.tofNSigmaPi() < nSigmaTOF_TPCTOF_Pi && trk.tpcNSigmaPi() < nSigmaTPC_TPCTOF_Pi && downsampleTsalisCharged(trk.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion], downsamplingTsalisPions)) { fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaPi(), trk.tofNSigmaPi(), trk.tpcExpSignalPi(), o2::track::PID::Pion, dwnSmplFactor_Pi); } } /// Loop tracks From a59b72bbbf91d581b2f48be8472fda076c28ce64 Mon Sep 17 00:00:00 2001 From: AnnalenaKa Date: Tue, 12 Jul 2022 09:35:35 +0200 Subject: [PATCH 03/11] Fix configurable name --- Common/Tools/tpcSkimsTableCreator.cxx | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/Common/Tools/tpcSkimsTableCreator.cxx b/Common/Tools/tpcSkimsTableCreator.cxx index b668aa855ca..a7b46d6e215 100644 --- a/Common/Tools/tpcSkimsTableCreator.cxx +++ b/Common/Tools/tpcSkimsTableCreator.cxx @@ -436,7 +436,7 @@ struct TreeWriterTPCTOF { /// 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{"downsamplingTsalisKaonsn", -1., "Downsampling factor to reduce the number of kaons"}; + 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){ @@ -453,8 +453,8 @@ struct TreeWriterTPCTOF { /// 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, double factor1Pt, double sqrts, double mass, float dwnsmplFact){ - if (dwnsmplFact<0.) return true; + bool downsampleTsalisCharged(double pt, float factor1Pt, double sqrts, double mass){ + if (factor1Pt<0.) return true; double prob=tsalisCharged(pt,mass,sqrts)*pt; double probNorm=tsalisCharged(1.,mass,sqrts); if ((fRndm->Rndm()*((prob/probNorm)*pt*pt))>factor1Pt) @@ -505,21 +505,21 @@ struct TreeWriterTPCTOF { rowTPCTOFTree.reserve(tracks.size()); for (auto const& trk : tracks) { /// Fill tree for protons - if (trk.tpcInnerParam() < maxMomTPCOnlyPr && trk.tpcNSigmaPr() < nSigmaTPCOnlyPr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton], downsamplingTsalisProtons)) { + if (trk.tpcInnerParam() < maxMomTPCOnlyPr && 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(), o2::track::PID::Proton, dwnSmplFactor_Pr); - } else if (trk.tpcInnerParam() > maxMomTPCOnlyPr && trk.tofNSigmaPr() < nSigmaTOF_TPCTOF_Pr && trk.tpcNSigmaPr() < nSigmaTPC_TPCTOF_Pr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton], downsamplingTsalisProtons)) { + } else if (trk.tpcInnerParam() > maxMomTPCOnlyPr && abs(trk.tofNSigmaPr()) < nSigmaTOF_TPCTOF_Pr && 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(), o2::track::PID::Proton, dwnSmplFactor_Pr); } /// Fill tree for kaons - if (trk.tpcInnerParam() < maxMomTPCOnlyKa && trk.tpcNSigmaKa() < nSigmaTPCOnlyKa && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon], downsamplingTsalisKaons)) { + if (trk.tpcInnerParam() < maxMomTPCOnlyKa && 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(), o2::track::PID::Kaon, dwnSmplFactor_Ka); - } else if (trk.tpcInnerParam() > maxMomTPCOnlyKa && trk.tofNSigmaKa() < nSigmaTOF_TPCTOF_Ka && trk.tpcNSigmaKa() < nSigmaTPC_TPCTOF_Ka && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon], downsamplingTsalisKaons)) { + } else if (trk.tpcInnerParam() > maxMomTPCOnlyKa && abs(trk.tofNSigmaKa()) < nSigmaTOF_TPCTOF_Ka && 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(), o2::track::PID::Kaon, dwnSmplFactor_Ka); } /// Fill tree pions - if (trk.tpcInnerParam() < maxMomTPCOnlyPi && trk.tpcNSigmaPi() < nSigmaTPCOnlyPi && downsampleTsalisCharged(trk.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion], downsamplingTsalisPions)) { + if (trk.tpcInnerParam() < maxMomTPCOnlyPi && 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(), o2::track::PID::Pion, dwnSmplFactor_Pi); - } else if (trk.tpcInnerParam() > maxMomTPCOnlyPi && trk.tofNSigmaPi() < nSigmaTOF_TPCTOF_Pi && trk.tpcNSigmaPi() < nSigmaTPC_TPCTOF_Pi && downsampleTsalisCharged(trk.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion], downsamplingTsalisPions)) { + } else if (trk.tpcInnerParam() > maxMomTPCOnlyPi && abs(trk.tofNSigmaPi()) < nSigmaTOF_TPCTOF_Pi && 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(), o2::track::PID::Pion, dwnSmplFactor_Pi); } } /// Loop tracks From b5ef5bb13d478b910a28d9d31408330d45208b02 Mon Sep 17 00:00:00 2001 From: AnnalenaKa Date: Fri, 15 Jul 2022 22:07:46 +0200 Subject: [PATCH 04/11] Correct NCl norm and default values tpc signal and reso --- Common/Core/PID/TPCPIDResponse.h | 21 +++- Common/TableProducer/PID/pidTPC.cxx | 4 +- Common/TableProducer/PID/pidTPCFull.cxx | 5 +- Common/TableProducer/PID/pidTPCML.cxx | 2 +- Common/TableProducer/PID/pidTPCML.h | 6 +- Common/Tools/CMakeLists.txt | 5 +- Common/Tools/handleParamTPCResponse.cxx | 3 + DPG/Tasks/TPC/CMakeLists.txt | 15 +++ .../Tasks/TPC}/tpcSkimsTableCreator.cxx | 104 +++++++++--------- 9 files changed, 98 insertions(+), 67 deletions(-) create mode 100644 DPG/Tasks/TPC/CMakeLists.txt rename {Common/Tools => DPG/Tasks/TPC}/tpcSkimsTableCreator.cxx (90%) diff --git a/Common/Core/PID/TPCPIDResponse.h b/Common/Core/PID/TPCPIDResponse.h index a99c0bb5cea..ebc9818f8b6 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,6 +93,7 @@ class Response float mChargeFactor = 2.299999952316284f; float mMultNormalization = 11000.; bool mUseDefaultResolutionParam = true; + float nClNorm = 152.f; ClassDefNV(Response, 2); @@ -98,6 +103,9 @@ class Response template inline float Response::GetExpectedSignal(const TrackType& track, const o2::track::PID::ID id) const { + if (!track.hasTPC()) { + return 0.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; } @@ -106,13 +114,16 @@ inline float Response::GetExpectedSignal(const TrackType& track, const o2::track 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,9 @@ 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 (!trk.hasTPC()) { + return -999.f; + } return ((trk.tpcSignal() - GetExpectedSignal(trk, id)) / GetExpectedSigma(collision, trk, id)); } @@ -173,6 +187,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/TableProducer/PID/pidTPC.cxx b/Common/TableProducer/PID/pidTPC.cxx index 4a69b2b31a9..fe8a3337640 100644 --- a/Common/TableProducer/PID/pidTPC.cxx +++ b/Common/TableProducer/PID/pidTPC.cxx @@ -156,6 +156,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()); } } } @@ -183,6 +184,7 @@ struct tpcPid { reserveTable(pidAl, tablePIDAl); std::vector network_prediction; + float nNclNormalization = response.GetNClNormalization(); if (useNetworkCorrection) { @@ -227,7 +229,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 769e2e30f2e..595e440a563 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 - + 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..96ffc83d57a 100644 --- a/Common/TableProducer/PID/pidTPCML.h +++ b/Common/TableProducer/PID/pidTPCML.h @@ -52,6 +52,8 @@ class Network // 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 @@ -64,6 +66,8 @@ class Network std::vector> mInputShapes; 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 fb3d4fff106..d7285a428a1 100644 --- a/Common/Tools/CMakeLists.txt +++ b/Common/Tools/CMakeLists.txt @@ -25,7 +25,4 @@ o2physics_add_executable(pidparam-tpc-response SOURCES handleParamTPCResponse.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore ) -o2physics_add_dpl_workflow(pid-tpc-skimscreation - SOURCES tpcSkimsTableCreator.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore - COMPONENT_NAME Analysis) + diff --git a/Common/Tools/handleParamTPCResponse.cxx b/Common/Tools/handleParamTPCResponse.cxx index b8a0bccf8c1..933e9076948 100644 --- a/Common/Tools/handleParamTPCResponse.cxx +++ b/Common/Tools/handleParamTPCResponse.cxx @@ -50,6 +50,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.), "Multiplicity Normalization")( "useDefaultParam", bpo::value()->default_value(true), "Use default sigma parametrisation")( "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)")( @@ -107,6 +108,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()) { @@ -190,6 +192,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/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/Common/Tools/tpcSkimsTableCreator.cxx b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx similarity index 90% rename from Common/Tools/tpcSkimsTableCreator.cxx rename to DPG/Tasks/TPC/tpcSkimsTableCreator.cxx index a7b46d6e215..fe3f46eea11 100644 --- a/Common/Tools/tpcSkimsTableCreator.cxx +++ b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx @@ -103,6 +103,7 @@ struct TreeWriterTpcV0 { /// 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"}; /// 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."}; @@ -175,7 +176,7 @@ struct TreeWriterTpcV0 { return false; } /// Pion downsampling - if(downsamplingTsalisPions>0. && downsampleTsalisCharged(track.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion])==0) { + if (downsamplingTsalisPions > 0. && downsampleTsalisCharged(track.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion]) == 0) { return false; } return true; @@ -216,7 +217,7 @@ struct TreeWriterTpcV0 { return false; } /// Proton downsampling - if(downsamplingTsalisProtons>0. && downsampleTsalisCharged(track.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton])==0) { + if (downsamplingTsalisProtons > 0. && downsampleTsalisCharged(track.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton]) == 0) { return false; } return true; @@ -286,7 +287,7 @@ struct TreeWriterTpcV0 { 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) { + if (downsamplingTsalisElectrons > 0. && downsampleTsalisCharged(track.pt(), downsamplingTsalisElectrons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Electron]) == 0) { return false; } return true; @@ -315,34 +316,35 @@ struct TreeWriterTpcV0 { mass, bg, multTPC / 11000., - std::sqrt(159. / ncl), + std::sqrt(nClNorm / ncl), id, nSigmaTPC, nSigmaTOF); } }; - - - 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=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); + 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 = 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); - int downsampleTsalisCharged(double pt, double factor1Pt, double sqrts, double mass){ - double prob=tsalisCharged(pt,mass,sqrts)*pt; - double probNorm=tsalisCharged(1.,mass,sqrts); - int triggerMask=0; - if ((fRndm->Rndm()*((prob/probNorm)*pt*pt))Rndm() * ((prob / probNorm) * pt * pt)) < factor1Pt) + triggerMask = 1; return triggerMask; }; @@ -350,7 +352,6 @@ struct TreeWriterTpcV0 { { } - void process(Coll::iterator const& collision, Trks const& tracks, aod::V0Datas const& v0s) { rowTPCTree.reserve(tracks.size()); @@ -361,45 +362,37 @@ struct TreeWriterTpcV0 { auto negTrack = v0.negTrack_as(); /// Fill table for Kaons if (selectionKaon(collision, v0)) { - if (selectionPion(posTrack)) - { + if (selectionPion(posTrack)) { fillSkimmedV0Table(posTrack, collision, posTrack.tpcNSigmaPi(), posTrack.tofNSigmaPi(), posTrack.tpcExpSignalPi(), o2::track::PID::Pion, dwnSmplFactor_Pi); } - if (selectionPion(negTrack)) - { + if (selectionPion(negTrack)) { fillSkimmedV0Table(negTrack, collision, negTrack.tpcNSigmaPi(), negTrack.tofNSigmaPi(), negTrack.tpcExpSignalPi(), o2::track::PID::Pion, dwnSmplFactor_Pi); } } /// Fill table for Lambdas if (selectionLambda(collision, v0)) { - if (selectionProton(posTrack)) - { + if (selectionProton(posTrack)) { fillSkimmedV0Table(posTrack, collision, posTrack.tpcNSigmaPr(), posTrack.tofNSigmaPr(), posTrack.tpcExpSignalPr(), o2::track::PID::Proton, dwnSmplFactor_Pr); } - if (selectionPion(negTrack)) - { + if (selectionPion(negTrack)) { fillSkimmedV0Table(negTrack, collision, negTrack.tpcNSigmaPi(), negTrack.tofNSigmaPi(), negTrack.tpcExpSignalPi(), o2::track::PID::Pion, dwnSmplFactor_Pi); } } /// Fill table for Antilambdas if (selectionAntiLambda(collision, v0)) { - if (selectionPion(posTrack)) - { + if (selectionPion(posTrack)) { fillSkimmedV0Table(posTrack, collision, posTrack.tpcNSigmaPi(), posTrack.tofNSigmaPi(), posTrack.tpcExpSignalPi(), o2::track::PID::Pion, dwnSmplFactor_Pi); } - if (selectionProton(negTrack)) - { + if (selectionProton(negTrack)) { fillSkimmedV0Table(negTrack, collision, negTrack.tpcNSigmaPr(), negTrack.tofNSigmaPr(), negTrack.tpcExpSignalPr(), o2::track::PID::Proton, dwnSmplFactor_Pr); } } /// Fill table for Gammas if (selectionGamma(collision, v0)) { - if (selectionElectron(posTrack)) - { + if (selectionElectron(posTrack)) { fillSkimmedV0Table(posTrack, collision, posTrack.tpcNSigmaEl(), posTrack.tofNSigmaEl(), posTrack.tpcExpSignalEl(), o2::track::PID::Electron, dwnSmplFactor_El); } - if (selectionElectron(negTrack)) - { + if (selectionElectron(negTrack)) { fillSkimmedV0Table(negTrack, collision, negTrack.tpcNSigmaEl(), negTrack.tofNSigmaEl(), negTrack.tpcExpSignalEl(), o2::track::PID::Electron, dwnSmplFactor_El); } } @@ -415,6 +408,7 @@ struct TreeWriterTPCTOF { Produces rowTPCTOFTree; /// Configurables + Configurable nClNorm{"nClNorm", 152., "Number of cluster normalization. Run 2: 159, Run 3 152"}; /// 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"}; @@ -439,32 +433,32 @@ struct TreeWriterTPCTOF { 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=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); + 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 = 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; - double prob=tsalisCharged(pt,mass,sqrts)*pt; - double probNorm=tsalisCharged(1.,mass,sqrts); - if ((fRndm->Rndm()*((prob/probNorm)*pt*pt))>factor1Pt) - { + bool downsampleTsalisCharged(double pt, float factor1Pt, double sqrts, double mass) + { + if (factor1Pt < 0.) + return true; + double prob = tsalisCharged(pt, mass, sqrts) * pt; + double probNorm = tsalisCharged(1., mass, sqrts); + if ((fRndm->Rndm() * ((prob / probNorm) * pt * pt)) > factor1Pt) { return false; - } - else { + } else { return true; - } - + } }; template @@ -490,7 +484,7 @@ struct TreeWriterTPCTOF { mass, bg, multTPC / 11000., - std::sqrt(159. / ncl), + std::sqrt(nClNorm / ncl), id, nSigmaTPC, nSigmaTOF); From fc19923dc3611e70b820f19ec068ba987bed4e11 Mon Sep 17 00:00:00 2001 From: AnnalenaKa Date: Tue, 19 Jul 2022 00:47:06 +0200 Subject: [PATCH 05/11] default resolution 0 and variables constant --- Common/Core/PID/TPCPIDResponse.h | 11 +++++--- Common/TableProducer/PID/pidTPC.cxx | 2 +- Common/TableProducer/PID/pidTPCFull.cxx | 6 ++--- Common/Tools/handleParamTPCResponse.cxx | 2 +- DPG/Tasks/CMakeLists.txt | 3 +++ DPG/Tasks/TPC/tpcSkimsTableCreator.cxx | 35 +++++++++++++------------ 6 files changed, 33 insertions(+), 26 deletions(-) diff --git a/Common/Core/PID/TPCPIDResponse.h b/Common/Core/PID/TPCPIDResponse.h index ebc9818f8b6..0f5a6c38423 100644 --- a/Common/Core/PID/TPCPIDResponse.h +++ b/Common/Core/PID/TPCPIDResponse.h @@ -95,7 +95,7 @@ class Response bool mUseDefaultResolutionParam = true; float nClNorm = 152.f; - ClassDefNV(Response, 2); + ClassDefNV(Response, 3); }; // class Response @@ -115,12 +115,12 @@ template inline float Response::GetExpectedSigma(const CollisionType& collision, const TrackType& track, const o2::track::PID::ID id) const { if (!track.hasTPC()) { - return -999.f; + return 0.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 = -999.f; + reso >= 0.f ? resolution = reso : resolution = 0.f; } else { const double ncl = nClNorm / track.tpcNClsFound(); // @@ -133,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 = -999.f; + reso >= 0.f ? resolution = reso : resolution = 0.f; } return resolution; } @@ -142,6 +142,9 @@ 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 (!trk.hasTPC()) { return -999.f; } diff --git a/Common/TableProducer/PID/pidTPC.cxx b/Common/TableProducer/PID/pidTPC.cxx index fe8a3337640..c313b780af3 100644 --- a/Common/TableProducer/PID/pidTPC.cxx +++ b/Common/TableProducer/PID/pidTPC.cxx @@ -184,7 +184,7 @@ struct tpcPid { reserveTable(pidAl, tablePIDAl); std::vector network_prediction; - float nNclNormalization = response.GetNClNormalization(); + const float nNclNormalization = response.GetNClNormalization(); if (useNetworkCorrection) { diff --git a/Common/TableProducer/PID/pidTPCFull.cxx b/Common/TableProducer/PID/pidTPCFull.cxx index 595e440a563..5426e105852 100644 --- a/Common/TableProducer/PID/pidTPCFull.cxx +++ b/Common/TableProducer/PID/pidTPCFull.cxx @@ -215,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 - float nNclNormalization = response.GetNClNormalization(); + const float nNclNormalization = response.GetNClNormalization(); float duration_network = 0; std::vector track_properties(track_prop_size); @@ -277,7 +277,7 @@ struct tpcPidFull { // For now only the option 2: sigma will be used. The other options are kept if there would be demand later on if (network.getOutputDimensions() == 1) { table(response.GetExpectedSignal(trk, pid), - response.GetExpectedSigma(collisions.iteratorAt(trk.collisionId()), trk, pid), + response.GetExpectedSigma(collisions.iteratorAt(trk.collisionId()), trk, pid), (trk.tpcSignal() - network_prediction[count_tracks + tracks_size * pid] * response.GetExpectedSignal(trk, pid)) / response.GetExpectedSigma(collisions.iteratorAt(trk.collisionId()), trk, pid)); } else if (network.getOutputDimensions() == 2) { table(response.GetExpectedSignal(trk, pid), (network_prediction[2 * (count_tracks + tracks_size * pid) + 1] - network_prediction[2 * (count_tracks + tracks_size * pid)]) * response.GetExpectedSignal(trk, pid), @@ -287,7 +287,7 @@ struct tpcPidFull { table(response.GetExpectedSignal(trk, pid), (network_prediction[3 * (count_tracks + tracks_size * pid) + 1] - network_prediction[3 * (count_tracks + tracks_size * pid)]) * response.GetExpectedSignal(trk, pid), (trk.tpcSignal() / response.GetExpectedSignal(trk, pid) - network_prediction[3 * (count_tracks + tracks_size * pid)]) / (network_prediction[3 * (count_tracks + tracks_size * pid) + 1] - network_prediction[3 * (count_tracks + tracks_size * pid)])); } else { - table(response.GetExpectedSignal(trk, pid),(network_prediction[3 * (count_tracks + tracks_size * pid)] - network_prediction[3 * (count_tracks + tracks_size * pid) + 2]) * response.GetExpectedSignal(trk, pid), + table(response.GetExpectedSignal(trk, pid), (network_prediction[3 * (count_tracks + tracks_size * pid)] - network_prediction[3 * (count_tracks + tracks_size * pid) + 2]) * response.GetExpectedSignal(trk, pid), (trk.tpcSignal() / response.GetExpectedSignal(trk, pid) - network_prediction[3 * (count_tracks + tracks_size * pid)]) / (network_prediction[3 * (count_tracks + tracks_size * pid)] - network_prediction[3 * (count_tracks + tracks_size * pid) + 2])); } } else { diff --git a/Common/Tools/handleParamTPCResponse.cxx b/Common/Tools/handleParamTPCResponse.cxx index 933e9076948..7544757c5ef 100644 --- a/Common/Tools/handleParamTPCResponse.cxx +++ b/Common/Tools/handleParamTPCResponse.cxx @@ -50,7 +50,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.), "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")( "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)")( diff --git a/DPG/Tasks/CMakeLists.txt b/DPG/Tasks/CMakeLists.txt index 2355a911e38..db6e8b5698a 100644 --- a/DPG/Tasks/CMakeLists.txt +++ b/DPG/Tasks/CMakeLists.txt @@ -71,3 +71,6 @@ o2physics_add_dpl_workflow(pid-tpc-qa-mc o2physics_add_dpl_workflow(fdd-qa SOURCES qaFDD.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase COMPONENT_NAME Analysis) + + +add_subdirectory(TPC) \ No newline at end of file diff --git a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx index fe3f46eea11..a4b85383c3c 100644 --- a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx +++ b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx @@ -147,10 +147,10 @@ struct TreeWriterTpcV0 { bool selectionKaon(C const& collision, V0 const& v0) { // initialise dynamic variables - float alpha = v0.alpha(); - float qt = v0.qtarm(); - float cosPA = v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); - float q_K = {cutAPK0pinner * sqrt(abs(1 - ((alpha * alpha) / (cutAPK0alinner * cutAPK0alinner))))}; + 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 * sqrt(abs(1 - ((alpha * alpha) / (cutAPK0alinner * cutAPK0alinner))))}; /// Armenteros-Podolanski cut if ((qt < cutQTK0min) || (qt > cutQTK0max) || (qt < q_K)) { return false; @@ -188,11 +188,11 @@ struct TreeWriterTpcV0 { { // initialise dynamic variables - float alpha = v0.alpha(); - float qt = v0.qtarm(); - float cosPA = v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); + const float alpha = v0.alpha(); + const float qt = v0.qtarm(); + const float cosPA = v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); - float q_L = cutAPL0 * sqrt(abs(1 - (((alpha - cutAPL1) * (alpha - cutAPL1)) / (cutAPL2 * cutAPL2)))); + const float q_L = cutAPL0 * sqrt(abs(1 - (((alpha - cutAPL1) * (alpha - cutAPL1)) / (cutAPL2 * cutAPL2)))); /// Armenteros-Podolanski cut if ((alpha < cutAlphaminL) || (alpha > cutAlphamaxL) || (qt < cutQTmin) || (qt > q_L)) { return false; @@ -228,11 +228,11 @@ struct TreeWriterTpcV0 { bool selectionAntiLambda(C const& collision, V0 const& v0) { // initialise dynamic variables - float alpha = v0.alpha(); - float qt = v0.qtarm(); - float cosPA = v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); + const float alpha = v0.alpha(); + const float qt = v0.qtarm(); + const float cosPA = v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); - float q_L = cutAPL0 * sqrt(abs(1 - (((alpha + cutAPL1) * (alpha + cutAPL1)) / (cutAPL2 * cutAPL2)))); + const float q_L = cutAPL0 * sqrt(abs(1 - (((alpha + cutAPL1) * (alpha + cutAPL1)) / (cutAPL2 * cutAPL2)))); /// Armenteros-Podolanski cut if ((alpha < cutAlphaminAL) || (alpha > cutAlphamaxAL) || (qt < cutQTmin) || (qt > q_L)) { return false; @@ -253,10 +253,10 @@ struct TreeWriterTpcV0 { bool selectionGamma(C const& collision, V0 const& v0) { // initialise dynamic variables - float alpha = v0.alpha(); - float qt = v0.qtarm(); - float cosPA = v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); - float pT = v0.pt(); + 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; @@ -450,8 +450,9 @@ struct TreeWriterTPCTOF { TRandom3* fRndm = new TRandom3(0); bool downsampleTsalisCharged(double pt, float factor1Pt, double sqrts, double mass) { - if (factor1Pt < 0.) + if (factor1Pt < 0.) { return true; + } double prob = tsalisCharged(pt, mass, sqrts) * pt; double probNorm = tsalisCharged(1., mass, sqrts); if ((fRndm->Rndm() * ((prob / probNorm) * pt * pt)) > factor1Pt) { From ceb25a9697800542201828d9578029c654f76d9e Mon Sep 17 00:00:00 2001 From: AnnalenaKa Date: Thu, 21 Jul 2022 09:49:34 +0200 Subject: [PATCH 06/11] change default values sigma and signal --- Common/Core/PID/TPCPIDResponse.h | 19 ++- Common/DataModel/PIDResponse.h | 45 +++--- Common/TableProducer/PID/pidTPCFull.cxx | 12 +- DPG/Tasks/TPC/tpcSkimsTableCreator.cxx | 176 +++++++++++++----------- DPG/Tasks/TPC/tpcSkimsTableCreator.h | 77 +++++++++++ 5 files changed, 220 insertions(+), 109 deletions(-) create mode 100644 DPG/Tasks/TPC/tpcSkimsTableCreator.h diff --git a/Common/Core/PID/TPCPIDResponse.h b/Common/Core/PID/TPCPIDResponse.h index 0f5a6c38423..58dfebcf58b 100644 --- a/Common/Core/PID/TPCPIDResponse.h +++ b/Common/Core/PID/TPCPIDResponse.h @@ -104,10 +104,10 @@ template inline float Response::GetExpectedSignal(const TrackType& track, const o2::track::PID::ID id) const { if (!track.hasTPC()) { - return 0.f; + 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 @@ -115,12 +115,12 @@ template inline float Response::GetExpectedSigma(const CollisionType& collision, const TrackType& track, const o2::track::PID::ID id) const { if (!track.hasTPC()) { - return 0.f; + 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 = nClNorm / track.tpcNClsFound(); // @@ -133,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; } @@ -145,6 +145,9 @@ inline float Response::GetNumberOfSigma(const CollisionType& collision, const Tr if (GetExpectedSigma(collision, trk, id) < 0.) { return -999.f; } + if (GetExpectedSignal(trk, id) < 0.) { + return -999.f; + } if (!trk.hasTPC()) { return -999.f; } @@ -155,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)); } diff --git a/Common/DataModel/PIDResponse.h b/Common/DataModel/PIDResponse.h index 92f344dfe56..8eebe118a4b 100644 --- a/Common/DataModel/PIDResponse.h +++ b/Common/DataModel/PIDResponse.h @@ -888,15 +888,24 @@ DECLARE_SOA_TABLE(pidTOFAl, "AOD", "pidTOFAl", //! Table of the TOF response wit namespace pidtpc { // Expected signals -DECLARE_SOA_COLUMN(TPCExpSignalEl, tpcExpSignalEl, float); //! Expected resolution with the TPC detector for electron -DECLARE_SOA_COLUMN(TPCExpSignalMu, tpcExpSignalMu, float); //! Expected resolution with the TPC detector for muon -DECLARE_SOA_COLUMN(TPCExpSignalPi, tpcExpSignalPi, float); //! Expected resolution with the TPC detector for pion -DECLARE_SOA_COLUMN(TPCExpSignalKa, tpcExpSignalKa, float); //! Expected resolution with the TPC detector for kaon -DECLARE_SOA_COLUMN(TPCExpSignalPr, tpcExpSignalPr, float); //! Expected resolution with the TPC detector for proton -DECLARE_SOA_COLUMN(TPCExpSignalDe, tpcExpSignalDe, float); //! Expected resolution with the TPC detector for deuteron -DECLARE_SOA_COLUMN(TPCExpSignalTr, tpcExpSignalTr, float); //! Expected resolution with the TPC detector for triton -DECLARE_SOA_COLUMN(TPCExpSignalHe, tpcExpSignalHe, float); //! Expected resolution with the TPC detector for helium3 -DECLARE_SOA_COLUMN(TPCExpSignalAl, tpcExpSignalAl, float); //! Expected resolution with the TPC detector for alpha +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; }); @@ -978,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::TPCExpSignalEl, 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::TPCExpSignalMu, 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::TPCExpSignalPi, 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::TPCExpSignalKa, 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::TPCExpSignalPr, 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::TPCExpSignalDe, 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::TPCExpSignalTr, 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::TPCExpSignalHe, 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::TPCExpSignalAl, 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/pidTPCFull.cxx b/Common/TableProducer/PID/pidTPCFull.cxx index 5426e105852..4baf507e920 100644 --- a/Common/TableProducer/PID/pidTPCFull.cxx +++ b/Common/TableProducer/PID/pidTPCFull.cxx @@ -276,26 +276,24 @@ struct tpcPidFull { // Here comes the application of the network. The output--dimensions of the network dtermine the application: 1: mean, 2: sigma, 3: sigma asymmetric // For now only the option 2: sigma will be used. The other options are kept if there would be demand later on if (network.getOutputDimensions() == 1) { - table(response.GetExpectedSignal(trk, pid), - response.GetExpectedSigma(collisions.iteratorAt(trk.collisionId()), trk, pid), + table(response.GetExpectedSigma(collisions.iteratorAt(trk.collisionId()), trk, pid), (trk.tpcSignal() - network_prediction[count_tracks + tracks_size * pid] * response.GetExpectedSignal(trk, pid)) / response.GetExpectedSigma(collisions.iteratorAt(trk.collisionId()), trk, pid)); } else if (network.getOutputDimensions() == 2) { - table(response.GetExpectedSignal(trk, pid), (network_prediction[2 * (count_tracks + tracks_size * pid) + 1] - network_prediction[2 * (count_tracks + tracks_size * pid)]) * response.GetExpectedSignal(trk, pid), + table((network_prediction[2 * (count_tracks + tracks_size * pid) + 1] - network_prediction[2 * (count_tracks + tracks_size * pid)]) * response.GetExpectedSignal(trk, pid), (trk.tpcSignal() / response.GetExpectedSignal(trk, pid) - network_prediction[2 * (count_tracks + tracks_size * pid)]) / (network_prediction[2 * (count_tracks + tracks_size * pid) + 1] - network_prediction[2 * (count_tracks + tracks_size * pid)])); } else if (network.getOutputDimensions() == 3) { if (trk.tpcSignal() / response.GetExpectedSignal(trk, pid) >= network_prediction[3 * (count_tracks + tracks_size * pid)]) { - table(response.GetExpectedSignal(trk, pid), (network_prediction[3 * (count_tracks + tracks_size * pid) + 1] - network_prediction[3 * (count_tracks + tracks_size * pid)]) * response.GetExpectedSignal(trk, pid), + table((network_prediction[3 * (count_tracks + tracks_size * pid) + 1] - network_prediction[3 * (count_tracks + tracks_size * pid)]) * response.GetExpectedSignal(trk, pid), (trk.tpcSignal() / response.GetExpectedSignal(trk, pid) - network_prediction[3 * (count_tracks + tracks_size * pid)]) / (network_prediction[3 * (count_tracks + tracks_size * pid) + 1] - network_prediction[3 * (count_tracks + tracks_size * pid)])); } else { - table(response.GetExpectedSignal(trk, pid), (network_prediction[3 * (count_tracks + tracks_size * pid)] - network_prediction[3 * (count_tracks + tracks_size * pid) + 2]) * response.GetExpectedSignal(trk, pid), + table((network_prediction[3 * (count_tracks + tracks_size * pid)] - network_prediction[3 * (count_tracks + tracks_size * pid) + 2]) * response.GetExpectedSignal(trk, pid), (trk.tpcSignal() / response.GetExpectedSignal(trk, pid) - network_prediction[3 * (count_tracks + tracks_size * pid)]) / (network_prediction[3 * (count_tracks + tracks_size * pid)] - network_prediction[3 * (count_tracks + tracks_size * pid) + 2])); } } else { LOGF(fatal, "Network output-dimensions incompatible!"); } } else { - table(response.GetExpectedSignal(trk, pid), - response.GetExpectedSigma(collisions.iteratorAt(trk.collisionId()), trk, pid), + table(response.GetExpectedSigma(collisions.iteratorAt(trk.collisionId()), trk, pid), response.GetNumberOfSigma(collisions.iteratorAt(trk.collisionId()), trk, pid)); } }; diff --git a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx index a4b85383c3c..9d8ef9a58ce 100644 --- a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx +++ b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx @@ -26,6 +26,9 @@ #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 @@ -36,66 +39,10 @@ using namespace o2::framework::expressions; using namespace o2::track; using namespace o2::dataformats; -// void customize(std::vector& workflowOptions) -// { -// std::vector options{// runtime customisation goes here -// {"useV0", VariantType::Int, 0, {"Use V0 information for QA"}}}; -// std::swap(workflowOptions, options); -// } - -namespace o2::aod -{ -namespace tpcskims -{ -DECLARE_SOA_COLUMN(InvDeDxExpV0, invdEdxExpV0, float); -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); -} // namespace tpcskims -DECLARE_SOA_TABLE(SkimmedTPCV0Tree, "AOD", "TPCSKIMV0TREE", - o2::aod::track::TPCSignal, - tpcskims::InvDeDxExpV0, - 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); - -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 - struct TreeWriterTpcV0 { using Trks = soa::Join; - using Coll = soa::Join; + using Coll = soa::Join; /// Tables to be produced Produces rowTPCTree; @@ -104,6 +51,7 @@ struct TreeWriterTpcV0 { 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."}; @@ -293,8 +241,9 @@ struct TreeWriterTpcV0 { return true; } - template - void fillSkimmedV0Table(T const& track, C const& collision, const float nSigmaTPC, const float nSigmaTOF, const float dEdxExp, const o2::track::PID::ID id, double dwnSmplFactor) + /// 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(); @@ -303,6 +252,11 @@ struct TreeWriterTpcV0 { 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(); + double pseudoRndm = track.pt() * 1000. - (long)(track.pt() * 1000); if (pseudoRndm < dwnSmplFactor) { rowTPCTree(track.tpcSignal(), @@ -319,7 +273,11 @@ struct TreeWriterTpcV0 { std::sqrt(nClNorm / ncl), id, nSigmaTPC, - nSigmaTOF); + nSigmaTOF, + alpha, + qt, + cosPA, + pT); } }; @@ -348,12 +306,33 @@ struct TreeWriterTpcV0 { 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 @@ -363,37 +342,37 @@ struct TreeWriterTpcV0 { /// Fill table for Kaons if (selectionKaon(collision, v0)) { if (selectionPion(posTrack)) { - fillSkimmedV0Table(posTrack, collision, posTrack.tpcNSigmaPi(), posTrack.tofNSigmaPi(), posTrack.tpcExpSignalPi(), o2::track::PID::Pion, dwnSmplFactor_Pi); + fillSkimmedV0Table(v0, posTrack, collision, posTrack.tpcNSigmaPi(), posTrack.tofNSigmaPi(), posTrack.tpcExpSignalPi(posTrack.tpcSignal()), o2::track::PID::Pion, dwnSmplFactor_Pi); } if (selectionPion(negTrack)) { - fillSkimmedV0Table(negTrack, collision, negTrack.tpcNSigmaPi(), negTrack.tofNSigmaPi(), negTrack.tpcExpSignalPi(), o2::track::PID::Pion, dwnSmplFactor_Pi); + 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(posTrack, collision, posTrack.tpcNSigmaPr(), posTrack.tofNSigmaPr(), posTrack.tpcExpSignalPr(), o2::track::PID::Proton, dwnSmplFactor_Pr); + fillSkimmedV0Table(v0, posTrack, collision, posTrack.tpcNSigmaPr(), posTrack.tofNSigmaPr(), posTrack.tpcExpSignalPr(posTrack.tpcSignal()), o2::track::PID::Proton, dwnSmplFactor_Pr); } if (selectionPion(negTrack)) { - fillSkimmedV0Table(negTrack, collision, negTrack.tpcNSigmaPi(), negTrack.tofNSigmaPi(), negTrack.tpcExpSignalPi(), o2::track::PID::Pion, dwnSmplFactor_Pi); + 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(posTrack, collision, posTrack.tpcNSigmaPi(), posTrack.tofNSigmaPi(), posTrack.tpcExpSignalPi(), o2::track::PID::Pion, dwnSmplFactor_Pi); + fillSkimmedV0Table(v0, posTrack, collision, posTrack.tpcNSigmaPi(), posTrack.tofNSigmaPi(), posTrack.tpcExpSignalPi(posTrack.tpcSignal()), o2::track::PID::Pion, dwnSmplFactor_Pi); } if (selectionProton(negTrack)) { - fillSkimmedV0Table(negTrack, collision, negTrack.tpcNSigmaPr(), negTrack.tofNSigmaPr(), negTrack.tpcExpSignalPr(), o2::track::PID::Proton, dwnSmplFactor_Pr); + 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(posTrack, collision, posTrack.tpcNSigmaEl(), posTrack.tofNSigmaEl(), posTrack.tpcExpSignalEl(), o2::track::PID::Electron, dwnSmplFactor_El); + fillSkimmedV0Table(v0, posTrack, collision, posTrack.tpcNSigmaEl(), posTrack.tofNSigmaEl(), posTrack.tpcExpSignalEl(posTrack.tpcSignal()), o2::track::PID::Electron, dwnSmplFactor_El); } if (selectionElectron(negTrack)) { - fillSkimmedV0Table(negTrack, collision, negTrack.tpcNSigmaEl(), negTrack.tofNSigmaEl(), negTrack.tpcExpSignalEl(), o2::track::PID::Electron, dwnSmplFactor_El); + fillSkimmedV0Table(v0, negTrack, collision, negTrack.tpcNSigmaEl(), negTrack.tofNSigmaEl(), negTrack.tpcExpSignalEl(negTrack.tpcSignal()), o2::track::PID::Electron, dwnSmplFactor_El); } } } /// Loop V0 candidates @@ -402,13 +381,14 @@ struct TreeWriterTpcV0 { struct TreeWriterTPCTOF { using Trks = soa::Join; - using Coll = 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"}; /// 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"}; @@ -462,6 +442,7 @@ struct TreeWriterTPCTOF { } }; + /// 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) { @@ -492,30 +473,71 @@ struct TreeWriterTPCTOF { } }; + /// 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 (!isTrackSelected(collision, trk)) { + continue; + } /// Fill tree for protons if (trk.tpcInnerParam() < maxMomTPCOnlyPr && 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(), o2::track::PID::Proton, dwnSmplFactor_Pr); + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaPr(), trk.tofNSigmaPr(), trk.tpcExpSignalPr(trk.tpcSignal()), o2::track::PID::Proton, dwnSmplFactor_Pr); } else if (trk.tpcInnerParam() > maxMomTPCOnlyPr && abs(trk.tofNSigmaPr()) < nSigmaTOF_TPCTOF_Pr && 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(), o2::track::PID::Proton, dwnSmplFactor_Pr); + 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 && 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(), o2::track::PID::Kaon, dwnSmplFactor_Ka); + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaKa(), trk.tofNSigmaKa(), trk.tpcExpSignalKa(trk.tpcSignal()), o2::track::PID::Kaon, dwnSmplFactor_Ka); } else if (trk.tpcInnerParam() > maxMomTPCOnlyKa && abs(trk.tofNSigmaKa()) < nSigmaTOF_TPCTOF_Ka && 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(), o2::track::PID::Kaon, dwnSmplFactor_Ka); + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaKa(), trk.tofNSigmaKa(), trk.tpcExpSignalKa(trk.tpcSignal()), o2::track::PID::Kaon, dwnSmplFactor_Ka); } /// Fill tree pions if (trk.tpcInnerParam() < maxMomTPCOnlyPi && 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(), o2::track::PID::Pion, dwnSmplFactor_Pi); + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaPi(), trk.tofNSigmaPi(), trk.tpcExpSignalPi(trk.tpcSignal()), o2::track::PID::Pion, dwnSmplFactor_Pi); } else if (trk.tpcInnerParam() > maxMomTPCOnlyPi && abs(trk.tofNSigmaPi()) < nSigmaTOF_TPCTOF_Pi && 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(), o2::track::PID::Pion, dwnSmplFactor_Pi); + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaPi(), trk.tofNSigmaPi(), trk.tpcExpSignalPi(trk.tpcSignal()), o2::track::PID::Pion, dwnSmplFactor_Pi); } } /// Loop tracks } /// process @@ -523,11 +545,7 @@ struct TreeWriterTPCTOF { WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - // const int useV0 = cfgc.options().get("useV0"); auto workflow = WorkflowSpec{adaptAnalysisTask(cfgc)}; workflow.push_back(adaptAnalysisTask(cfgc)); - // if (useV0) { - - // } 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..b69ef4a8e76 --- /dev/null +++ b/DPG/Tasks/TPC/tpcSkimsTableCreator.h @@ -0,0 +1,77 @@ +// 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); +} // 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); + +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 From 6770ec9a4c603b13793f9615f0fb99818bc5d29b Mon Sep 17 00:00:00 2001 From: AnnalenaKa Date: Thu, 21 Jul 2022 11:43:10 +0200 Subject: [PATCH 07/11] Fix formatting errors --- Common/DataModel/PIDResponse.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Common/DataModel/PIDResponse.h b/Common/DataModel/PIDResponse.h index 8eebe118a4b..3865c1a1c8a 100644 --- a/Common/DataModel/PIDResponse.h +++ b/Common/DataModel/PIDResponse.h @@ -999,7 +999,7 @@ DECLARE_SOA_TABLE(pidTPCFullPr, "AOD", "pidTPCFullPr", //! Table of the TPC (ful DECLARE_SOA_TABLE(pidTPCFullDe, "AOD", "pidTPCFullDe", //! Table of the TPC (full) response with expected signal, expected resolution and Nsigma for deuteron 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::TPCExpSignalTr, 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::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 From 46f1e4310ff69805c2d62e7daa2408e9464bb69f Mon Sep 17 00:00:00 2001 From: AnnalenaKa Date: Thu, 21 Jul 2022 15:15:56 +0200 Subject: [PATCH 08/11] declare variables const --- DPG/Tasks/TPC/tpcSkimsTableCreator.cxx | 32 +++++++++++++++----------- DPG/Tasks/TPC/tpcSkimsTableCreator.h | 6 ++++- 2 files changed, 24 insertions(+), 14 deletions(-) diff --git a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx index 9d8ef9a58ce..73819f28fa9 100644 --- a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx +++ b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx @@ -256,8 +256,10 @@ struct TreeWriterTpcV0 { 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(); - double pseudoRndm = track.pt() * 1000. - (long)(track.pt() * 1000); + const double pseudoRndm = track.pt() * 1000. - (long)(track.pt() * 1000); if (pseudoRndm < dwnSmplFactor) { rowTPCTree(track.tpcSignal(), 1. / dEdxExp, @@ -277,7 +279,9 @@ struct TreeWriterTpcV0 { alpha, qt, cosPA, - pT); + pT, + v0radius, + gammapsipair); } }; @@ -285,11 +289,11 @@ struct TreeWriterTpcV0 { { const double a = 6.81, b = 59.24; const double c = 0.082, d = 0.151; - double mt = 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); + 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; }; @@ -298,11 +302,13 @@ struct TreeWriterTpcV0 { TRandom3* fRndm = new TRandom3(0); int downsampleTsalisCharged(double pt, double factor1Pt, double sqrts, double mass) { - double prob = tsalisCharged(pt, mass, sqrts) * pt; - double probNorm = tsalisCharged(1., mass, sqrts); + 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) + if ((fRndm->Rndm() * ((prob / probNorm) * pt * pt)) < factor1Pt) { triggerMask = 1; + } + return triggerMask; }; @@ -433,8 +439,8 @@ struct TreeWriterTPCTOF { if (factor1Pt < 0.) { return true; } - double prob = tsalisCharged(pt, mass, sqrts) * pt; - double probNorm = tsalisCharged(1., mass, sqrts); + 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 { @@ -453,7 +459,7 @@ struct TreeWriterTPCTOF { const double bg = p / mass; const int multTPC = collision.multTPC(); - double pseudoRndm = track.pt() * 1000. - (long)(track.pt() * 1000); + const double pseudoRndm = track.pt() * 1000. - (long)(track.pt() * 1000); if (pseudoRndm < dwnSmplFactor) { rowTPCTOFTree(track.tpcSignal(), 1. / dEdxExp, diff --git a/DPG/Tasks/TPC/tpcSkimsTableCreator.h b/DPG/Tasks/TPC/tpcSkimsTableCreator.h index b69ef4a8e76..c06efc4f6e6 100644 --- a/DPG/Tasks/TPC/tpcSkimsTableCreator.h +++ b/DPG/Tasks/TPC/tpcSkimsTableCreator.h @@ -36,6 +36,8 @@ 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, @@ -56,7 +58,9 @@ DECLARE_SOA_TABLE(SkimmedTPCV0Tree, "AOD", "TPCSKIMV0TREE", tpcskims::AlphaV0, tpcskims::QtV0, tpcskims::CosPAV0, - tpcskims::PtV0); + tpcskims::PtV0, + tpcskims::RadiusV0, + tpcskims::GammaPsiPair); DECLARE_SOA_TABLE(SkimmedTPCTOFTree, "AOD", "TPCTOFSKIMTREE", o2::aod::track::TPCSignal, From e1c78fdd88b81c67ed38213b49f7a8ebc0a76e1e Mon Sep 17 00:00:00 2001 From: AnnalenaKa Date: Thu, 21 Jul 2022 15:46:54 +0200 Subject: [PATCH 09/11] remove TMath and use std --- DPG/Tasks/TPC/tpcSkimsTableCreator.cxx | 36 +++++++++++++------------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx index 73819f28fa9..7eed941b68a 100644 --- a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx +++ b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx @@ -98,7 +98,7 @@ struct TreeWriterTpcV0 { 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 * sqrt(abs(1 - ((alpha * alpha) / (cutAPK0alinner * cutAPK0alinner))))}; + 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; @@ -108,7 +108,7 @@ struct TreeWriterTpcV0 { return false; } /// Invariant mass cut - if (abs(v0.mK0Short() - 0.497611) > invariantMassCutK0Short) { + if (std::abs(v0.mK0Short() - 0.497611) > invariantMassCutK0Short) { return false; } @@ -120,7 +120,7 @@ struct TreeWriterTpcV0 { bool selectionPion(T const& track) { /// TOF cut daughter tracks - if (nSigmaTOFdautrack != 0. && abs(track.tofNSigmaPi()) > nSigmaTOFdautrack) { + if (nSigmaTOFdautrack != 0. && std::abs(track.tofNSigmaPi()) > nSigmaTOFdautrack) { return false; } /// Pion downsampling @@ -140,7 +140,7 @@ struct TreeWriterTpcV0 { const float qt = v0.qtarm(); const float cosPA = v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); - const float q_L = cutAPL0 * sqrt(abs(1 - (((alpha - cutAPL1) * (alpha - cutAPL1)) / (cutAPL2 * cutAPL2)))); + 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; @@ -150,7 +150,7 @@ struct TreeWriterTpcV0 { return false; } /// Invariant mass cut - if (abs(v0.mLambda() - 1.115683) > invariantMassCutLambda) { + if (std::abs(v0.mLambda() - 1.115683) > invariantMassCutLambda) { return false; } return true; @@ -161,7 +161,7 @@ struct TreeWriterTpcV0 { bool selectionProton(T const& track) { /// TOF cut daughter tracks - if (nSigmaTOFdautrack != 0. && abs(track.tofNSigmaPr()) > nSigmaTOFdautrack) { + if (nSigmaTOFdautrack != 0. && std::abs(track.tofNSigmaPr()) > nSigmaTOFdautrack) { return false; } /// Proton downsampling @@ -180,7 +180,7 @@ struct TreeWriterTpcV0 { const float qt = v0.qtarm(); const float cosPA = v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); - const float q_L = cutAPL0 * sqrt(abs(1 - (((alpha + cutAPL1) * (alpha + cutAPL1)) / (cutAPL2 * cutAPL2)))); + 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; @@ -190,7 +190,7 @@ struct TreeWriterTpcV0 { return false; } /// Invariant mass cut - if (abs(v0.mAntiLambda() - 1.115683) > invariantMassCutLambda) { + if (std::abs(v0.mAntiLambda() - 1.115683) > invariantMassCutLambda) { return false; } return true; @@ -211,10 +211,10 @@ struct TreeWriterTpcV0 { } /// Armenteros-Podolanski cut - if ((((qt > gammaQtMax) || (TMath::Abs(alpha) > cutAlphaG1)) || ((qt > cutQTG) || (TMath::Abs(alpha) < cutAlphaG2min) || (TMath::Abs(alpha) > cutAlphaG2max)))) { + if ((((qt > gammaQtMax) || (std::abs(alpha) > cutAlphaG1)) || ((qt > cutQTG) || (std::abs(alpha) < cutAlphaG2min) || (std::abs(alpha) > cutAlphaG2max)))) { return false; } - if ((TMath::Power(alpha / gammaAsymmetryMax, 2) + TMath::Power(qt / lQtMaxPtDep, 2)) < 1) { + if ((std::pow(alpha / gammaAsymmetryMax, 2) + std::pow(qt / lQtMaxPtDep, 2)) < 1) { return false; } if (v0.v0radius() < gammaV0RadiusMin || v0.v0radius() > gammaV0RadiusMax) { @@ -224,7 +224,7 @@ struct TreeWriterTpcV0 { if (cosPA < cutPAV0) { return false; } - if ((TMath::Abs(v0.psipair()) < gammaPsiPairMax)) { + if ((std::abs(v0.psipair()) < gammaPsiPairMax)) { return false; } return true; @@ -423,7 +423,7 @@ struct TreeWriterTPCTOF { { const double a = 6.81, b = 59.24; const double c = 0.082, d = 0.151; - double mt = sqrt(mass * mass + pt * pt); + double mt = std::sqrt(mass * mass + pt * pt); double n = a + b / sqrts; double T = c + d / sqrts; double p0 = n * T; @@ -528,21 +528,21 @@ struct TreeWriterTPCTOF { continue; } /// Fill tree for protons - if (trk.tpcInnerParam() < maxMomTPCOnlyPr && abs(trk.tpcNSigmaPr()) < nSigmaTPCOnlyPr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton])) { + 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 && abs(trk.tofNSigmaPr()) < nSigmaTOF_TPCTOF_Pr && abs(trk.tpcNSigmaPr()) < nSigmaTPC_TPCTOF_Pr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton])) { + } 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 && abs(trk.tpcNSigmaKa()) < nSigmaTPCOnlyKa && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon])) { + 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 && abs(trk.tofNSigmaKa()) < nSigmaTOF_TPCTOF_Ka && abs(trk.tpcNSigmaKa()) < nSigmaTPC_TPCTOF_Ka && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon])) { + } 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 && abs(trk.tpcNSigmaPi()) < nSigmaTPCOnlyPi && downsampleTsalisCharged(trk.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion])) { + 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 && abs(trk.tofNSigmaPi()) < nSigmaTOF_TPCTOF_Pi && abs(trk.tpcNSigmaPi()) < nSigmaTPC_TPCTOF_Pi && downsampleTsalisCharged(trk.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion])) { + } 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 From 7cb17aab99c6fcca22ceabf1295f611261998e3a Mon Sep 17 00:00:00 2001 From: AnnalenaKa Date: Tue, 26 Jul 2022 14:18:20 +0200 Subject: [PATCH 10/11] Add configurable track selection --- DPG/Tasks/TPC/tpcSkimsTableCreator.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx index 7eed941b68a..28e78a0c348 100644 --- a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx +++ b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx @@ -395,6 +395,7 @@ struct TreeWriterTPCTOF { /// 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"}; @@ -524,7 +525,7 @@ struct TreeWriterTPCTOF { for (auto const& trk : tracks) { /// Check track selection - if (!isTrackSelected(collision, trk)) { + if (applyTrkSel == 1 && !isTrackSelected(collision, trk)) { continue; } /// Fill tree for protons From 42b26fc4d15d14aa4a2eb5ae08af865ceaa6dbd2 Mon Sep 17 00:00:00 2001 From: AnnalenaKa Date: Thu, 28 Jul 2022 16:28:33 +0200 Subject: [PATCH 11/11] Fix formatting errors --- Common/TableProducer/PID/pidTPCML.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Common/TableProducer/PID/pidTPCML.h b/Common/TableProducer/PID/pidTPCML.h index 96ffc83d57a..8ea794346fc 100644 --- a/Common/TableProducer/PID/pidTPCML.h +++ b/Common/TableProducer/PID/pidTPCML.h @@ -46,8 +46,8 @@ 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]; }; @@ -66,7 +66,7 @@ class Network std::vector> mInputShapes; 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