From 6167c3495cceaadb66f18034b9d0591d2db63326 Mon Sep 17 00:00:00 2001 From: glromane <95305986+glromane@users.noreply.github.com> Date: Fri, 22 Dec 2023 18:03:41 +0100 Subject: [PATCH] adding binning for some variables and adding an MC table (#4225) * removing lines with table reservation and fixing linter bugs * Update femto3dPairTask.h fixing double-track cut * Update femto3dPairTask.h * Update femto3dPairTask.h * Update femto3dPairTask.cxx * Update femto3dQA.cxx * Update femto3dPairTask.cxx * Update femto3dPairTaskMC.cxx * Update femto3dPairTask.cxx --- PWGCF/Femto3D/CMakeLists.txt | 4 +- PWGCF/Femto3D/Core/CMakeLists.txt | 10 + PWGCF/Femto3D/Core/femto3dPairTask.h | 230 +++++++++++ PWGCF/Femto3D/DataModel/CMakeLists.txt | 0 PWGCF/Femto3D/DataModel/singletrackselector.h | 240 ++++++++--- PWGCF/Femto3D/TableProducer/CMakeLists.txt | 0 .../TableProducer/singleTrackSelector.cxx | 113 ++++-- PWGCF/Femto3D/Tasks/CMakeLists.txt | 25 ++ PWGCF/Femto3D/Tasks/femto3dPairTask.cxx | 382 ++++++++++++++++++ PWGCF/Femto3D/Tasks/femto3dPairTaskMC.cxx | 375 +++++++++++++++++ PWGCF/Femto3D/Tasks/femto3dQA.cxx | 221 ++++++++++ 11 files changed, 1503 insertions(+), 97 deletions(-) mode change 100644 => 100755 PWGCF/Femto3D/CMakeLists.txt create mode 100755 PWGCF/Femto3D/Core/CMakeLists.txt create mode 100755 PWGCF/Femto3D/Core/femto3dPairTask.h mode change 100644 => 100755 PWGCF/Femto3D/DataModel/CMakeLists.txt mode change 100644 => 100755 PWGCF/Femto3D/TableProducer/CMakeLists.txt create mode 100755 PWGCF/Femto3D/Tasks/CMakeLists.txt create mode 100755 PWGCF/Femto3D/Tasks/femto3dPairTask.cxx create mode 100755 PWGCF/Femto3D/Tasks/femto3dPairTaskMC.cxx create mode 100755 PWGCF/Femto3D/Tasks/femto3dQA.cxx diff --git a/PWGCF/Femto3D/CMakeLists.txt b/PWGCF/Femto3D/CMakeLists.txt old mode 100644 new mode 100755 index 7c83285fff0..e57cd178dd8 --- a/PWGCF/Femto3D/CMakeLists.txt +++ b/PWGCF/Femto3D/CMakeLists.txt @@ -9,7 +9,7 @@ # granted to it by virtue of its status as an Intergovernmental Organization # or submit itself to any jurisdiction. -#add_subdirectory(Core) +add_subdirectory(Core) add_subdirectory(DataModel) add_subdirectory(TableProducer) -#add_subdirectory(Tasks) +add_subdirectory(Tasks) \ No newline at end of file diff --git a/PWGCF/Femto3D/Core/CMakeLists.txt b/PWGCF/Femto3D/Core/CMakeLists.txt new file mode 100755 index 00000000000..62f88c324cf --- /dev/null +++ b/PWGCF/Femto3D/Core/CMakeLists.txt @@ -0,0 +1,10 @@ +# 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. \ No newline at end of file diff --git a/PWGCF/Femto3D/Core/femto3dPairTask.h b/PWGCF/Femto3D/Core/femto3dPairTask.h new file mode 100755 index 00000000000..c1ff8ecd3ed --- /dev/null +++ b/PWGCF/Femto3D/Core/femto3dPairTask.h @@ -0,0 +1,230 @@ +// 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. +// +/// \brief utility functions for femto task +/// \author Sofia Tomassini, Gleb Romanenko, Nicolò Jacazio +/// \since 30 May 2023 + +#ifndef PWGCF_FEMTO3D_CORE_FEMTO3DPAIRTASK_H_ +#define PWGCF_FEMTO3D_CORE_FEMTO3DPAIRTASK_H_ + +// #include "Framework/ASoA.h" +// #include "Framework/DataTypes.h" +// #include "Framework/AnalysisDataModel.h" +// #include "Common/DataModel/PIDResponse.h" +// #include "Framework/Logger.h" +// #include "Common/DataModel/Multiplicity.h" + +#include +#include "TLorentzVector.h" +#include "TVector3.h" +#include "TDatabasePDG.h" + +double particle_mass(int PDGcode) +{ + // if(PDGcode == 2212) return TDatabasePDG::Instance()->GetParticle(2212)->Mass(); + if (PDGcode == 1000010020) + return 1.87561294257; + else + return TDatabasePDG::Instance()->GetParticle(PDGcode)->Mass(); +} + +namespace o2::aod::singletrackselector +{ +template +Type getBinIndex(float const& value, std::vector const& binning, int const& NsubBins = 1) +{ + Type res = -100; + for (int i = 0; i < binning.size() - 1; i++) { + if (value >= binning[i] && binning[i + 1] > value) { + if (NsubBins < 2) { + res = (Type)i; + break; + } else { + float subBinWidth = (binning[i + 1] - binning[i]) / NsubBins; + int subBin = std::floor((value - binning[i]) / subBinWidth); + int delimeter = std::pow(10, std::to_string(NsubBins).size()); + + res = (Type)i + (Type)subBin / delimeter; + break; + } + } + } + return res; +} + +//==================================================================================== + +float GetKstarFrom4vectors(TLorentzVector& first4momentum, TLorentzVector& second4momentum, bool isIdentical) +{ + if (isIdentical) { + TLorentzVector fourmomentadiff = first4momentum - second4momentum; + return 0.5 * abs(fourmomentadiff.Mag()); + } else { + TLorentzVector fourmomentasum = first4momentum + second4momentum; + + first4momentum.Boost((-1) * fourmomentasum.BoostVector()); + second4momentum.Boost((-1) * fourmomentasum.BoostVector()); + + TVector3 qinv = first4momentum.Vect() - second4momentum.Vect(); + return 0.5 * abs(qinv.Mag()); + } +} + +//==================================================================================== + +template +class FemtoPair +{ + public: + FemtoPair() {} + FemtoPair(TrackType const& first, TrackType const& second) + { + _first = first; + _second = second; + } + FemtoPair(TrackType const& first, TrackType const& second, const bool& isidentical) + { + _first = first; + _second = second; + _isidentical = isidentical; + } + + FemtoPair(const FemtoPair& obj) + { + SetFirstParticle(obj.GetFirstParticle()); + SetSecondParticle(obj.GetSecondParticle()); + } + explicit FemtoPair(const FemtoPair* obj) + { + SetFirstParticle(obj->GetFirstParticle()); + SetSecondParticle(obj->GetSecondParticle()); + } + ~FemtoPair() {} + FemtoPair& operator=(const FemtoPair& obj) + { + if (this != &obj) { + SetFirstParticle(obj.GetFirstParticle()); + SetSecondParticle(obj.GetSecondParticle()); + } + return *this; + } + + void SetPair(TrackType const& first, TrackType const& second) + { + _first = first; + _second = second; + } + void SetFirstParticle(TrackType const& first) { _first = first; } + void SetSecondParticle(TrackType const& second) { _second = second; } + void SetIdentical(const bool& isidentical) { _isidentical = isidentical; } + void SetMagField1(const float& magfield1) { _magfield1 = magfield1; } + void SetMagField2(const float& magfield2) { _magfield2 = magfield2; } + void SetPDG1(const int& PDG1) { _PDG1 = PDG1; } + void SetPDG2(const int& PDG2) { _PDG2 = PDG2; } + void ResetPair(); + void ResetAll(); + + TrackType* GetFirstParticle() const { return _first; } + TrackType* GetSecondParticle() const { return _second; } + bool IsIdentical() { return _isidentical; } + + bool IsClosePair(const float& deta = 0.01, const float& dphi = 0.01, const float& radius = 1.2); + float GetEtaDiff() const + { + if (_first != NULL && _second != NULL) + return _first->eta() - _second->eta(); + else + return 1000; + } + float GetPhiStarDiff(const float& radius = 1.2) const + { + if (_first != NULL && _second != NULL) + return _first->phiStar(_magfield1, radius) - _second->phiStar(_magfield2, radius); + else + return 1000; + } + float GetKstar() const; + float GetKt() const; + + private: + TrackType _first = NULL; + TrackType _second = NULL; + float _magfield1 = 0.0, _magfield2 = 0.0; + int _PDG1 = 0, _PDG2 = 0; + bool _isidentical = true; +}; + +template +void FemtoPair::ResetPair() +{ + _first = NULL; + _second = NULL; +} + +template +void FemtoPair::ResetAll() +{ + _first = NULL; + _second = NULL; + _magfield1 = 0.0; + _magfield2 = 0.0; + _PDG1 = 0; + _PDG2 = 0; + _isidentical = true; +} + +template +bool FemtoPair::IsClosePair(const float& deta, const float& dphi, const float& radius) +{ + if (_first == NULL || _second == NULL) + return true; + if (!(_magfield1 * _magfield2)) + return true; + if (abs(GetEtaDiff()) < deta && abs(GetPhiStarDiff(radius)) < dphi) + return true; + + return false; +} + +template +float FemtoPair::GetKstar() const +{ + if (_first == NULL || _second == NULL) + return -1000; + if (!(_magfield1 * _magfield2)) + return -1000; + if (!(_PDG1 * _PDG2)) + return -1000; + + TLorentzVector first4momentum; + first4momentum.SetPtEtaPhiM(_first->pt(), _first->eta(), _first->phi(), particle_mass(_PDG1)); + TLorentzVector second4momentum; + second4momentum.SetPtEtaPhiM(_second->pt(), _second->eta(), _second->phi(), particle_mass(_PDG2)); + + return GetKstarFrom4vectors(first4momentum, second4momentum, _isidentical); +} + +template +float FemtoPair::GetKt() const +{ + if (_first == NULL || _second == NULL) + return -1000; + if (!(_magfield1 * _magfield2)) + return -1000; + if (!(_PDG1 * _PDG2)) + return -1000; + + return 0.5 * std::sqrt(std::pow(_first->px() + _second->px(), 2) + std::pow(_first->py() + _second->py(), 2)); +} +} // namespace o2::aod::singletrackselector + +#endif // PWGCF_FEMTO3D_CORE_FEMTO3DPAIRTASK_H_ diff --git a/PWGCF/Femto3D/DataModel/CMakeLists.txt b/PWGCF/Femto3D/DataModel/CMakeLists.txt old mode 100644 new mode 100755 diff --git a/PWGCF/Femto3D/DataModel/singletrackselector.h b/PWGCF/Femto3D/DataModel/singletrackselector.h index 91dfb63e856..3db909f2911 100755 --- a/PWGCF/Femto3D/DataModel/singletrackselector.h +++ b/PWGCF/Femto3D/DataModel/singletrackselector.h @@ -39,29 +39,65 @@ inline typename binningType::binned_t packInTable(const float& valueToBin) } else if (valueToBin >= binningType::binned_max) { return binningType::overflowBin; } else { - return static_cast(valueToBin / binningType::bin_width); + return static_cast((valueToBin - binningType::binned_center) / binningType::bin_width); } } template inline float unPack(const typename binningType::binned_t& b) { - return static_cast(binningType::bin_width * b); + return binningType::bin_width * b + binningType::binned_center; } -namespace nsigma +template +inline o2::framework::expressions::Node unPack(const T& b) { -struct binning { + return binningType::bin_width * b + binningType::binned_center; +} + +namespace binning +{ +struct binningParent { public: typedef int8_t binned_t; static constexpr int nbins = (1 << 8 * sizeof(binned_t)) - 2; static constexpr binned_t overflowBin = nbins >> 1; static constexpr binned_t underflowBin = -(nbins >> 1); - static constexpr float binned_max = 10.0; +}; + +struct nsigma : binningParent { + public: static constexpr float binned_min = -10.0; + static constexpr float binned_max = 10.0; + static constexpr float binned_center = 0.5 * (binned_min + binned_max); + static constexpr float bin_width = (binned_max - binned_min) / nbins; +}; + +struct dca : binningParent { + public: + static constexpr float binned_min = -1.0; + static constexpr float binned_max = 1.0; + static constexpr float binned_center = 0.5 * (binned_min + binned_max); static constexpr float bin_width = (binned_max - binned_min) / nbins; }; -} // namespace nsigma + +struct chi2 : binningParent { + public: + static constexpr float binned_min = 0.0; + static constexpr float binned_max = 10.0; + static constexpr float binned_center = 0.5 * (binned_min + binned_max); + static constexpr float bin_width = (binned_max - binned_min) / nbins; +}; + +struct rowsOverFindable : binningParent { + public: + static constexpr float binned_min = 0.0; + static constexpr float binned_max = 3.0; + static constexpr float binned_center = 0.5 * (binned_min + binned_max); + static constexpr float bin_width = (binned_max - binned_min) / nbins; +}; + +} // namespace binning DECLARE_SOA_COLUMN(Mult, mult, int); // Multiplicity of the collision DECLARE_SOA_COLUMN(MultPercentile, multPerc, float); // Percentiles of multiplicity of the collision @@ -79,27 +115,29 @@ DECLARE_SOA_TABLE(SingleCollSels, "AOD", "SINGLECOLLSEL", // Table of the variab namespace singletrackselector { - -DECLARE_SOA_INDEX_COLUMN(SingleCollSel, singleCollSel); // Index to the collision -DECLARE_SOA_COLUMN(P, p, float); // Momentum of the track -DECLARE_SOA_COLUMN(DcaXY, dcaXY, float); // impact parameter of the track -DECLARE_SOA_COLUMN(DcaZ, dcaZ, float); // impact parameter of the track -DECLARE_SOA_COLUMN(TPCInnerParam, tpcInnerParam, float); // Momentum at inner wall of the TPC -DECLARE_SOA_COLUMN(TPCSignal, tpcSignal, float); // dE/dx TPC -DECLARE_SOA_COLUMN(Beta, beta, float); // TOF beta -DECLARE_SOA_COLUMN(TPCNClsFound, tpcNClsFound, int16_t); // Number of TPC clusters -DECLARE_SOA_COLUMN(TPCChi2NCl, tpcChi2NCl, float); // TPC chi2 -DECLARE_SOA_COLUMN(TPCCrossedRowsOverFindableCls, tpcCrossedRowsOverFindableCls, float); // Ratio of found over findable clusters -DECLARE_SOA_COLUMN(TPCNClsShared, tpcNClsShared, uint8_t); // Number of shared TPC clusters -DECLARE_SOA_COLUMN(ITSNCls, itsNCls, uint8_t); // Number of ITS clusters -DECLARE_SOA_COLUMN(ITSChi2NCl, itsChi2NCl, float); // ITS chi2 -DECLARE_SOA_COLUMN(Sign, sign, int8_t); +DECLARE_SOA_INDEX_COLUMN(SingleCollSel, singleCollSel); // Index to the collision +DECLARE_SOA_COLUMN(P, p, float); // Momentum of the track DECLARE_SOA_COLUMN(Eta, eta, float); DECLARE_SOA_COLUMN(Phi, phi, float); -DECLARE_SOA_COLUMN(StoredTOFNSigmaPr, storedTofNSigmaPr, nsigma::binning::binned_t); -DECLARE_SOA_COLUMN(StoredTPCNSigmaPr, storedTpcNSigmaPr, nsigma::binning::binned_t); -DECLARE_SOA_COLUMN(StoredTOFNSigmaDe, storedTofNSigmaDe, nsigma::binning::binned_t); -DECLARE_SOA_COLUMN(StoredTPCNSigmaDe, storedTpcNSigmaDe, nsigma::binning::binned_t); +DECLARE_SOA_COLUMN(Sign, sign, int8_t); +DECLARE_SOA_COLUMN(TPCNClsFound, tpcNClsFound, int16_t); // Number of TPC clusters +DECLARE_SOA_COLUMN(TPCNClsShared, tpcNClsShared, uint8_t); // Number of shared TPC clusters +DECLARE_SOA_COLUMN(ITSNCls, itsNCls, uint8_t); // Number of ITS clusters + +DECLARE_SOA_COLUMN(StoredDcaXY, storedDcaXY, binning::dca::binned_t); // impact parameter of the track +DECLARE_SOA_COLUMN(StoredDcaZ, storedDcaZ, binning::dca::binned_t); // impact parameter of the track +DECLARE_SOA_COLUMN(StoredTPCChi2NCl, storedTpcChi2NCl, binning::chi2::binned_t); // TPC chi2 +DECLARE_SOA_COLUMN(StoredITSChi2NCl, storedItsChi2NCl, binning::chi2::binned_t); // ITS chi2 +DECLARE_SOA_COLUMN(StoredTPCCrossedRowsOverFindableCls, storedTpcCrossedRowsOverFindableCls, binning::rowsOverFindable::binned_t); // Ratio of found over findable clusters + +DECLARE_SOA_COLUMN(StoredTOFNSigmaPi, storedTofNSigmaPi, binning::nsigma::binned_t); +DECLARE_SOA_COLUMN(StoredTPCNSigmaPi, storedTpcNSigmaPi, binning::nsigma::binned_t); +DECLARE_SOA_COLUMN(StoredTOFNSigmaKa, storedTofNSigmaKa, binning::nsigma::binned_t); +DECLARE_SOA_COLUMN(StoredTPCNSigmaKa, storedTpcNSigmaKa, binning::nsigma::binned_t); +DECLARE_SOA_COLUMN(StoredTOFNSigmaPr, storedTofNSigmaPr, binning::nsigma::binned_t); +DECLARE_SOA_COLUMN(StoredTPCNSigmaPr, storedTpcNSigmaPr, binning::nsigma::binned_t); +DECLARE_SOA_COLUMN(StoredTOFNSigmaDe, storedTofNSigmaDe, binning::nsigma::binned_t); +DECLARE_SOA_COLUMN(StoredTPCNSigmaDe, storedTpcNSigmaDe, binning::nsigma::binned_t); DECLARE_SOA_DYNAMIC_COLUMN(Energy, energy, [](float p, float mass) -> float { return sqrt(p * p + mass * mass); }); DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, [](float p, float eta) -> float { return p / std::cosh(eta); }); @@ -115,22 +153,39 @@ DECLARE_SOA_DYNAMIC_COLUMN(PhiStar, phiStar, } }); +DECLARE_SOA_DYNAMIC_COLUMN(DcaXY, dcaXY, + [](binning::dca::binned_t dca_binned) -> float { return singletrackselector::unPack(dca_binned); }); +DECLARE_SOA_DYNAMIC_COLUMN(DcaZ, dcaZ, + [](binning::dca::binned_t dca_binned) -> float { return singletrackselector::unPack(dca_binned); }); +DECLARE_SOA_DYNAMIC_COLUMN(TPCChi2NCl, tpcChi2NCl, + [](binning::chi2::binned_t chi2_binned) -> float { return singletrackselector::unPack(chi2_binned); }); +DECLARE_SOA_DYNAMIC_COLUMN(ITSChi2NCl, itsChi2NCl, + [](binning::chi2::binned_t chi2_binned) -> float { return singletrackselector::unPack(chi2_binned); }); + +DECLARE_SOA_DYNAMIC_COLUMN(TPCCrossedRowsOverFindableCls, tpcCrossedRowsOverFindableCls, + [](binning::rowsOverFindable::binned_t rowsOverFindable_binned) -> float { return singletrackselector::unPack(rowsOverFindable_binned); }); + +DECLARE_SOA_DYNAMIC_COLUMN(TOFNSigmaPi, tofNSigmaPi, + [](binning::nsigma::binned_t nsigma_binned) -> float { return singletrackselector::unPack(nsigma_binned); }); +DECLARE_SOA_DYNAMIC_COLUMN(TPCNSigmaPi, tpcNSigmaPi, + [](binning::nsigma::binned_t nsigma_binned) -> float { return singletrackselector::unPack(nsigma_binned); }); +DECLARE_SOA_DYNAMIC_COLUMN(TOFNSigmaKa, tofNSigmaKa, + [](binning::nsigma::binned_t nsigma_binned) -> float { return singletrackselector::unPack(nsigma_binned); }); +DECLARE_SOA_DYNAMIC_COLUMN(TPCNSigmaKa, tpcNSigmaKa, + [](binning::nsigma::binned_t nsigma_binned) -> float { return singletrackselector::unPack(nsigma_binned); }); + DECLARE_SOA_DYNAMIC_COLUMN(TOFNSigmaPr, tofNSigmaPr, - [](nsigma::binning::binned_t nsigma_binned) -> float { - return singletrackselector::unPack(nsigma_binned); - }); + [](binning::nsigma::binned_t nsigma_binned) -> float { return singletrackselector::unPack(nsigma_binned); }); DECLARE_SOA_DYNAMIC_COLUMN(TPCNSigmaPr, tpcNSigmaPr, - [](nsigma::binning::binned_t nsigma_binned) -> float { - return singletrackselector::unPack(nsigma_binned); - }); + [](binning::nsigma::binned_t nsigma_binned) -> float { return singletrackselector::unPack(nsigma_binned); }); DECLARE_SOA_DYNAMIC_COLUMN(TOFNSigmaDe, tofNSigmaDe, - [](nsigma::binning::binned_t nsigma_binned) -> float { - return singletrackselector::unPack(nsigma_binned); - }); + [](binning::nsigma::binned_t nsigma_binned) -> float { return singletrackselector::unPack(nsigma_binned); }); DECLARE_SOA_DYNAMIC_COLUMN(TPCNSigmaDe, tpcNSigmaDe, - [](nsigma::binning::binned_t nsigma_binned) -> float { - return singletrackselector::unPack(nsigma_binned); - }); + [](binning::nsigma::binned_t nsigma_binned) -> float { return singletrackselector::unPack(nsigma_binned); }); + +DECLARE_SOA_COLUMN(TPCInnerParam, tpcInnerParam, float); // Momentum at inner wall of the TPC +DECLARE_SOA_COLUMN(TPCSignal, tpcSignal, float); // dE/dx TPC +DECLARE_SOA_COLUMN(Beta, beta, float); // TOF beta } // namespace singletrackselector @@ -138,36 +193,82 @@ DECLARE_SOA_TABLE_FULL(SingleTrackSels, "SelTracks", "AOD", "SINGLETRACKSEL", // o2::soa::Index<>, singletrackselector::SingleCollSelId, singletrackselector::P, - singletrackselector::DcaXY, - singletrackselector::DcaZ, - singletrackselector::TPCInnerParam, - singletrackselector::TPCSignal, - singletrackselector::Beta, + singletrackselector::Eta, + singletrackselector::Phi, + singletrackselector::Sign, singletrackselector::TPCNClsFound, - singletrackselector::TPCChi2NCl, - singletrackselector::TPCCrossedRowsOverFindableCls, singletrackselector::TPCNClsShared, singletrackselector::ITSNCls, - singletrackselector::ITSChi2NCl, - singletrackselector::Sign, - singletrackselector::Eta, - singletrackselector::Phi, + singletrackselector::StoredDcaXY, + singletrackselector::StoredDcaZ, + singletrackselector::StoredTPCChi2NCl, + singletrackselector::StoredITSChi2NCl, + singletrackselector::StoredTPCCrossedRowsOverFindableCls, + + singletrackselector::StoredTOFNSigmaPi, + singletrackselector::StoredTPCNSigmaPi, + singletrackselector::StoredTOFNSigmaKa, + singletrackselector::StoredTPCNSigmaKa, singletrackselector::StoredTOFNSigmaPr, singletrackselector::StoredTPCNSigmaPr, singletrackselector::StoredTOFNSigmaDe, singletrackselector::StoredTPCNSigmaDe, + + singletrackselector::DcaXY, + singletrackselector::DcaZ, + singletrackselector::TPCChi2NCl, + singletrackselector::ITSChi2NCl, + singletrackselector::TPCCrossedRowsOverFindableCls, + + singletrackselector::TOFNSigmaPi, + singletrackselector::TPCNSigmaPi, + singletrackselector::TOFNSigmaKa, + singletrackselector::TPCNSigmaKa, + singletrackselector::TOFNSigmaPr, + singletrackselector::TPCNSigmaPr, + singletrackselector::TOFNSigmaDe, + singletrackselector::TPCNSigmaDe, + singletrackselector::Energy, singletrackselector::Pt, singletrackselector::Px, singletrackselector::Py, singletrackselector::Pz, - singletrackselector::PhiStar, - singletrackselector::TOFNSigmaPr, - singletrackselector::TPCNSigmaPr, - singletrackselector::TOFNSigmaDe, - singletrackselector::TPCNSigmaDe); + singletrackselector::PhiStar); + +DECLARE_SOA_TABLE(SingleTrkExtras, "AOD", "SINGLETRKEXTRA", + singletrackselector::TPCInnerParam, + singletrackselector::TPCSignal, + singletrackselector::Beta); + +namespace singletrackselector +{ +DECLARE_SOA_COLUMN(PdgCode, pdgCode, int); +DECLARE_SOA_COLUMN(Origin, origin, int); // 0 - prymary; 1 - weak decay; 2 - material +DECLARE_SOA_COLUMN(P_MC, p_MC, float); +DECLARE_SOA_COLUMN(Eta_MC, eta_MC, float); +DECLARE_SOA_COLUMN(Phi_MC, phi_MC, float); + +DECLARE_SOA_DYNAMIC_COLUMN(Pt_MC, pt_MC, [](float p, float eta) -> float { return p / std::cosh(eta); }); +DECLARE_SOA_DYNAMIC_COLUMN(Px_MC, px_MC, [](float p, float eta, float phi) -> float { return (p / std::cosh(eta)) * std::sin(phi); }); +DECLARE_SOA_DYNAMIC_COLUMN(Py_MC, py_MC, [](float p, float eta, float phi) -> float { return (p / std::cosh(eta)) * std::cos(phi); }); +DECLARE_SOA_DYNAMIC_COLUMN(Pz_MC, pz_MC, [](float p, float eta) -> float { return p * std::tanh(eta); }); + +} // namespace singletrackselector + +DECLARE_SOA_TABLE(SingleTrkMCs, "AOD", "SINGLETRKMC", // Table with generatad info from MC + singletrackselector::PdgCode, + singletrackselector::Origin, + singletrackselector::P_MC, + singletrackselector::Eta_MC, + singletrackselector::Phi_MC, + singletrackselector::Pt_MC, + singletrackselector::Px_MC, + singletrackselector::Py_MC, + singletrackselector::Pz_MC); } // namespace o2::aod + #endif // PWGCF_FEMTO3D_DATAMODEL_SINGLETRACKSELECTOR_H_ namespace o2::aod::singletrackselector @@ -176,15 +277,23 @@ namespace o2::aod::singletrackselector template inline bool TPCselection(TrackType const& track, std::pair> const& PIDcuts) { + int PDG = PIDcuts.first; float Nsigma = -1000; - - switch (PIDcuts.first) { + switch (PDG) { case 2212: Nsigma = track.tpcNSigmaPr(); break; case 1000010020: Nsigma = track.tpcNSigmaDe(); break; + case 211: + Nsigma = track.tpcNSigmaPi(); + break; + case 321: + Nsigma = track.tpcNSigmaKa(); + break; + case 0: + return false; default: LOG(fatal) << "Cannot interpret PDG for TPC selection: " << PIDcuts.first; } @@ -196,11 +305,14 @@ inline bool TPCselection(TrackType const& track, std::pair -inline bool TOFselection(TrackType const& track, std::pair> const& PIDcuts) +inline bool TOFselection(TrackType const& track, std::pair> const& PIDcuts, float const& TPCresidualCut = 5.0f) { - float Nsigma = -1000; + int PDG = PIDcuts.first; + if (!TPCselection(track, std::make_pair(PDG, std::vector{-TPCresidualCut, TPCresidualCut}))) + return false; - switch (PIDcuts.first) { + float Nsigma = -1000; + switch (PDG) { case 2212: Nsigma = track.tofNSigmaPr(); break; @@ -208,13 +320,13 @@ inline bool TOFselection(TrackType const& track, std::pair::value) { - Nsigma = track.tofNSigmaPi(); - } + Nsigma = track.tofNSigmaPi(); + break; case 321: - if constexpr (std::experimental::is_detected::value) { - Nsigma = track.tofNSigmaKa(); - } + Nsigma = track.tofNSigmaKa(); + break; + case 0: + return false; default: LOG(fatal) << "Cannot interpret PDG for TOF selection: " << PIDcuts.first; } diff --git a/PWGCF/Femto3D/TableProducer/CMakeLists.txt b/PWGCF/Femto3D/TableProducer/CMakeLists.txt old mode 100644 new mode 100755 diff --git a/PWGCF/Femto3D/TableProducer/singleTrackSelector.cxx b/PWGCF/Femto3D/TableProducer/singleTrackSelector.cxx index 4ebc3723246..89b21efb5ff 100755 --- a/PWGCF/Femto3D/TableProducer/singleTrackSelector.cxx +++ b/PWGCF/Femto3D/TableProducer/singleTrackSelector.cxx @@ -55,29 +55,34 @@ struct singleTrackSelector { Configurable> _particlesToReject{"particlesToRejectPDGs", std::vector{}, "PDG codes of perticles that will be rejected with TOF (only pion, kaon, proton and deurton are supported now)"}; Configurable> rejectWithinNsigmaTOF{"rejectWithinNsigmaTOF", std::vector{-5.0f, 5.0f}, "TOF rejection Nsigma range for particles specified with PDG to be rejected"}; - Configurable _tpcNClsFound{"tpcNClsFound", 70, "Minimun number of cluster found in TPC"}; - Configurable _itsNCls{"itsNCls", 0, "Minimun number of cluster found in ITS"}; - Configurable _tpcCrossedRowsOverFindableCls{"tpcCrossedRowsOverFindableCls", 0.f, "Minimun number of crossed rows over findable clusters"}; + Configurable _min_P{"min_P", 0.f, "lower mometum limit"}; + Configurable _max_P{"max_P", 100.f, "upper mometum limit"}; + Configurable _eta{"eta", 100.f, "abs eta value limit"}; Configurable _dcaXY{"dcaXY", 1000.f, "Maximum dca of track in xy"}; Configurable _dcaZ{"dcaZ", 1000.f, "Maximum dca of track in xy"}; using Trks = soa::Join; - using Coll = soa::Join; + + using Coll = soa::Join; Produces tableRow; Produces tableRowColl; + Produces tableRowExtra; + Produces tableRowMC; Filter eventFilter = (applyEvSel.node() == 0) || ((applyEvSel.node() == 1) && (aod::evsel::sel7 == true)) || ((applyEvSel.node() == 2) && (aod::evsel::sel8 == true)); Filter vertexFilter = ((o2::aod::collision::posZ < 15.f) && (o2::aod::collision::posZ > -15.f)); Filter trackFilter = ((o2::aod::track::itsChi2NCl <= 36.f) && (o2::aod::track::itsChi2NCl >= 0.f) && (o2::aod::track::tpcChi2NCl >= 0.f) && (o2::aod::track::tpcChi2NCl <= 4.f)); - Filter tpcFilter = ((o2::aod::singletrackselector::tpcNClsFound >= _tpcNClsFound) && (o2::aod::singletrackselector::tpcCrossedRowsOverFindableCls >= _tpcCrossedRowsOverFindableCls)); - Filter itsFilter = (o2::aod::singletrackselector::itsNCls >= _itsNCls); + + Filter pFilter = o2::aod::track::p > _min_P&& o2::aod::track::p < _max_P; + Filter etaFilter = nabs(o2::aod::track::eta) < _eta; Filter dcaFilter = (o2::aod::track::dcaXY <= _dcaXY) && (o2::aod::track::dcaZ <= _dcaZ); int mRunNumber = 0; @@ -127,14 +132,10 @@ struct singleTrackSelector { d_bz = 0.1 * d_bz; } - void process(soa::Filtered::iterator const& collision, soa::Filtered const& tracks, aod::BCsWithTimestamps const&) + template + inline void fillTheTables(Col collision, Trks const& tracks) { - bool skip_track = false; // flag used for track rejection - tableRow.reserve(tracks.size()); - - auto bc = collision.bc_as(); - initCCDB(bc); tableRowColl(collision.multTPC(), collision.centFT0M(), @@ -142,8 +143,14 @@ struct singleTrackSelector { d_bz); for (auto& track : tracks) { + if constexpr (isMC) { + if (!track.has_mcParticle()) + continue; + } skip_track = false; - + if (Configurable{"rejectNotPropagatedTrks", true, "rejects tracks that are not propagated to the primary vertex"} && track.trackType() != aod::track::Track) { + continue; + } for (auto i : particlesToReject) { // if satisfied, want to continue in the upper loop (over tracks) -- skip the current track // cannot use simple 'continue' since it will be applied to the current loop, so have to use a flag @@ -161,29 +168,73 @@ struct singleTrackSelector { tableRow(tableRowColl.lastIndex(), track.p(), - track.dcaXY(), - track.dcaZ(), - track.tpcInnerParam(), - track.tpcSignal(), - track.beta(), + track.eta(), + track.phi(), + track.sign(), track.tpcNClsFound(), - track.tpcChi2NCl(), - track.tpcCrossedRowsOverFindableCls(), track.tpcNClsShared(), track.itsNCls(), - track.itsChi2NCl(), - track.sign(), - track.eta(), - track.phi(), - singletrackselector::packInTable(track.tofNSigmaPr()), - singletrackselector::packInTable(track.tpcNSigmaPr()), - singletrackselector::packInTable(track.tofNSigmaDe()), - singletrackselector::packInTable(track.tpcNSigmaDe())); + singletrackselector::packInTable(track.dcaXY()), + singletrackselector::packInTable(track.dcaZ()), + singletrackselector::packInTable(track.tpcChi2NCl()), + singletrackselector::packInTable(track.itsChi2NCl()), + singletrackselector::packInTable(track.tpcCrossedRowsOverFindableCls()), + singletrackselector::packInTable(track.tofNSigmaPi()), + singletrackselector::packInTable(track.tpcNSigmaPi()), + singletrackselector::packInTable(track.tofNSigmaKa()), + singletrackselector::packInTable(track.tpcNSigmaKa()), + singletrackselector::packInTable(track.tofNSigmaPr()), + singletrackselector::packInTable(track.tpcNSigmaPr()), + singletrackselector::packInTable(track.tofNSigmaDe()), + singletrackselector::packInTable(track.tpcNSigmaDe())); + + tableRowExtra(track.tpcInnerParam(), + track.tpcSignal(), + track.beta()); + + if constexpr (isMC) { + int origin = -1; + if (track.mcParticle().isPhysicalPrimary()) + origin = 0; + if (!track.mcParticle().isPhysicalPrimary() && track.mcParticle().producedByGenerator()) + origin = 1; + if (!track.mcParticle().isPhysicalPrimary() && !track.mcParticle().producedByGenerator()) + origin = 2; + + if (origin == -1) + LOGF(fatal, "Could not define the origin (primary/weak decay/material) of the track!!!"); + + tableRowMC(track.mcParticle().pdgCode(), + origin, + track.mcParticle().p(), + track.mcParticle().eta(), + track.mcParticle().phi()); + } break; // break the loop with particlesToKeep after the 'if' condition is satisfied -- don't want double entries } } } + + void processData(soa::Filtered::iterator const& collision, soa::Filtered const& tracks, aod::BCsWithTimestamps const&) + { + + auto bc = collision.bc_as(); + initCCDB(bc); + + fillTheTables(collision, tracks); + } + PROCESS_SWITCH(singleTrackSelector, processData, "process data", true); + + void processMC(soa::Filtered::iterator const& collision, soa::Filtered> const& tracks, aod::McParticles const&, aod::BCsWithTimestamps const&) + { + + auto bc = collision.bc_as(); + initCCDB(bc); + + fillTheTables(collision, tracks); + } + PROCESS_SWITCH(singleTrackSelector, processMC, "process MC", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGCF/Femto3D/Tasks/CMakeLists.txt b/PWGCF/Femto3D/Tasks/CMakeLists.txt new file mode 100755 index 00000000000..a0338dce678 --- /dev/null +++ b/PWGCF/Femto3D/Tasks/CMakeLists.txt @@ -0,0 +1,25 @@ +# 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(femto3d-qa + SOURCES femto3dQA.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(femto3d-pair-task + SOURCES femto3dPairTask.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(femto3d-pair-task-mc + SOURCES femto3dPairTaskMC.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore + COMPONENT_NAME Analysis) \ No newline at end of file diff --git a/PWGCF/Femto3D/Tasks/femto3dPairTask.cxx b/PWGCF/Femto3D/Tasks/femto3dPairTask.cxx new file mode 100755 index 00000000000..d469eb78bf9 --- /dev/null +++ b/PWGCF/Femto3D/Tasks/femto3dPairTask.cxx @@ -0,0 +1,382 @@ +// 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. +// +/// \brief Femto3D pair mixing task +/// \author Sofia Tomassini, Gleb Romanenko, Nicolò Jacazio +/// \since 31 May 2023 + +#include +#include +#include + +#include "PWGCF/Femto3D/Core/femto3dPairTask.h" +#include "PWGCF/Femto3D/DataModel/singletrackselector.h" +#include "TLorentzVector.h" + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/ASoA.h" +#include "Framework/DataTypes.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/Expressions.h" +#include "Framework/StaticFor.h" +#include "MathUtils/Utils.h" +#include "Common/DataModel/Multiplicity.h" + +using namespace o2; +using namespace o2::soa; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct FemtoCorrelations { + // using allinfo = soa::Join; // aod::pidTPCPr + /// Construct a registry object with direct declaration + HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; + + Configurable _min_P{"min_P", 0.0, "lower mometum limit"}; + Configurable _max_P{"max_P", 100.0, "upper mometum limit"}; + Configurable _eta{"eta", 100.0, "abs eta value limit"}; + Configurable _dcaXY{"dcaXY", 10.0, "abs dcaXY value limit"}; + Configurable _dcaZ{"dcaZ", 10.0, "abs dcaZ value limit"}; + Configurable _tpcNClsFound{"minTpcNClsFound", 0, "minimum allowed number of TPC clasters"}; + Configurable _tpcChi2NCl{"tpcChi2NCl", 100.0, "upper limit for chi2 value of a fit over TPC clasters"}; + Configurable _tpcCrossedRowsOverFindableCls{"tpcCrossedRowsOverFindableCls", 0, "lower limit of TPC CrossedRows/FindableCls value"}; + Configurable _tpcNClsShared{"maxTpcNClsShared", 100, "maximum allowed number of TPC shared clasters"}; + Configurable _itsNCls{"minItsNCls", 0, "minimum allowed number of ITS clasters"}; + Configurable _itsChi2NCl{"itsChi2NCl", 100.0, "upper limit for chi2 value of a fit over ITS clasters"}; + Configurable _vertexZ{"VertexZ", 10.0, "abs vertexZ value limit"}; + + Configurable _sign_1{"sign_1", 1, "sign of the first particle in a pair"}; + Configurable _particlePDG_1{"particlePDG_1", 2212, "PDG code of the first particle in a pair to perform PID for (only proton and deurton are supported now)"}; + Configurable> _tpcNSigma_1{"tpcNSigma_1", std::vector{-3.0f, 3.0f}, "first particle PID: Nsigma range in TPC before the TOF is used"}; + Configurable _PIDtrshld_1{"PIDtrshld_1", 10.0, "first particle PID: value of momentum from which the PID is done with TOF (before that only TPC is used)"}; + Configurable> _tofNSigma_1{"tofNSigma_1", std::vector{-3.0f, 3.0f}, "first particle PID: Nsigma range in TOF"}; + + Configurable _sign_2{"sign_2", 1, "sign of the second particle in a pair"}; + Configurable _particlePDG_2{"particlePDG_2", 2212, "PDG code of the second particle in a pair to perform PID for (only proton and deurton are supported now)"}; + Configurable> _tpcNSigma_2{"tpcNSigma_2", std::vector{-3.0f, 3.0f}, "second particle PID: Nsigma range in TPC before the TOF is used"}; + Configurable _PIDtrshld_2{"PIDtrshld_2", 10.0, "second particle PID: value of momentum from which the PID is done with TOF (before that only TPC is used)"}; + Configurable> _tofNSigma_2{"tofNSigma_2", std::vector{-3.0f, 3.0f}, "second particle PID: Nsigma range in TOF"}; + + Configurable _particlePDGtoReject{"particlePDGtoRejectFromSecond", 0, "applied only if the particles are non-identical and only to the second particle in the pair!!!"}; + Configurable> _rejectWithinNsigmaTOF{"rejectWithinNsigmaTOF", std::vector{-0.0f, 0.0f}, "TOF rejection Nsigma range for the particle specified with PDG to be rejected"}; + + Configurable _deta{"deta", 0.01, "minimum allowed defference in eta between two tracks in a pair"}; + Configurable _dphi{"dphi", 0.01, "minimum allowed defference in phi_star between two tracks in a pair"}; + Configurable _radiusTPC{"radiusTPC", 1.2, "TPC radius to calculate phi_star for"}; + + Configurable _vertexNbinsToMix{"vertexNbinsToMix", 10, "Number of vertexZ bins for the mixing"}; + Configurable> _centBins{"multBins", std::vector{0.0f, 100.0f}, "multiplicity percentile/centrality binning (min:0, max:100)"}; + Configurable _multNsubBins{"multSubBins", 1, "number of sub-bins to perform the mixing within"}; + Configurable> _kTbins{"kTbins", std::vector{0.0f, 100.0f}, "pair transverse momentum kT binning"}; + ConfigurableAxis CFkStarBinning{"CFkStarBinning", {500, 0.005, 5.005}, "k* binning of the CF (Nbins, lowlimit, uplimit)"}; + + bool IsIdentical; + + std::pair> TPCcuts_1; + std::pair> TOFcuts_1; + + std::pair> TPCcuts_2; + std::pair> TOFcuts_2; + + using FilteredCollisions = aod::SingleCollSels; + using FilteredTracks = aod::SingleTrackSels; + + typedef std::shared_ptr::iterator> trkType; + typedef std::shared_ptr::iterator> colType; + + std::map> selectedtracks_1; + std::map> selectedtracks_2; + std::map, std::vector> mixbins; + + std::unique_ptr> Pair = std::make_unique>(); + + Filter pFilter = o2::aod::singletrackselector::p > _min_P&& o2::aod::singletrackselector::p < _max_P; + Filter etaFilter = nabs(o2::aod::singletrackselector::eta) < _eta; + + Filter tpcTrkFilter = o2::aod::singletrackselector::tpcNClsFound >= _tpcNClsFound && + o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedTpcChi2NCl) < _tpcChi2NCl && + o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedTpcCrossedRowsOverFindableCls) > _tpcCrossedRowsOverFindableCls; + + Filter dcaFilter = nabs(o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedDcaXY)) < _dcaXY && + nabs(o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedDcaZ)) < _dcaZ; + + Filter itsTrkFilter = o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedItsChi2NCl) < _itsChi2NCl; + + Filter vertexFilter = nabs(o2::aod::singletrackselector::posZ) < _vertexZ; + + std::vector> MultHistos; + std::vector>> kThistos; + std::vector>> SEhistos; + std::vector>> MEhistos; + + void init(o2::framework::InitContext&) + { + + if (_centBins.value.size() < 2) + LOGF(fatal, "The configured number of multiplicity/centrality bins in the array is less than 2 !!!"); + if (_kTbins.value.size() < 2) + LOGF(fatal, "The configured number of kT bins in the array is less than 2 !!!"); + if (_vertexNbinsToMix.value < 1) + LOGF(fatal, "The configured number of VertexZ bins is less than 1 !!!"); + + IsIdentical = (_sign_1 * _particlePDG_1 == _sign_2 * _particlePDG_2); + + Pair->SetIdentical(IsIdentical); + Pair->SetPDG1(_particlePDG_1); + Pair->SetPDG2(_particlePDG_2); + + TPCcuts_1 = std::make_pair(_particlePDG_1, _tpcNSigma_1); + TOFcuts_1 = std::make_pair(_particlePDG_1, _tofNSigma_1); + TPCcuts_2 = std::make_pair(_particlePDG_2, _tpcNSigma_2); + TOFcuts_2 = std::make_pair(_particlePDG_2, _tofNSigma_2); + + const AxisSpec kStarAxis{CFkStarBinning, "k* (GeV/c)"}; + + for (int i = 0; i < _centBins.value.size() - 1; i++) { + std::vector> SEperMult; + std::vector> MEperMult; + std::vector> kTperMult; + + auto hMult = registry.add(Form("Cent%i/Mult_vs_cent%i", i, i), Form("Mult_vs_cent%i", i), kTH1F, {{5001, -0.5, 5000.5, "Nch"}}); + MultHistos.push_back(std::move(hMult)); + + for (int j = 0; j < _kTbins.value.size() - 1; j++) { + auto hSE = registry.add(Form("Cent%i/SE_cent%i_kT%i", i, i, j), Form("SE_cent%i_kT%i", i, j), kTH1F, {kStarAxis}); + auto hME = registry.add(Form("Cent%i/ME_cent%i_kT%i", i, i, j), Form("ME_cent%i_kT%i", i, j), kTH1F, {kStarAxis}); + auto hkT = registry.add(Form("Cent%i/kT_cent%i_kT%i", i, i, j), Form("kT_cent%i_kT%i", i, j), kTH1F, {{500, 0., 5., "kT"}}); + SEperMult.push_back(std::move(hSE)); + MEperMult.push_back(std::move(hME)); + kTperMult.push_back(std::move(hkT)); + } + + SEhistos.push_back(std::move(SEperMult)); + MEhistos.push_back(std::move(MEperMult)); + kThistos.push_back(std::move(kTperMult)); + } + + registry.add("p_first", Form("p_%i", static_cast(_particlePDG_1)), kTH1F, {{100, 0., 5., "p"}}); + registry.add("nsigmaTOF_first", Form("nsigmaTOF_%i", static_cast(_particlePDG_1)), kTH2F, {{100, 0., 5.}, {100, -10., 10.}}); + registry.add("nsigmaTPC_first", Form("nsigmaTPC_%i", static_cast(_particlePDG_1)), kTH2F, {{100, 0., 5.}, {100, -10., 10.}}); + if (!IsIdentical) { + registry.add("p_second", Form("p_%i", static_cast(_particlePDG_2)), kTH1F, {{100, 0., 5., "p"}}); + registry.add("nsigmaTOF_second", Form("nsigmaTOF_%i", static_cast(_particlePDG_2)), kTH2F, {{100, 0., 5.}, {100, -10., 10.}}); + registry.add("nsigmaTPC_second", Form("nsigmaTPC_%i", static_cast(_particlePDG_2)), kTH2F, {{100, 0., 5.}, {100, -10., 10.}}); + } + } + + template + void mixTracks(Type const& tracks, int multBin) + { // template for identical particles from the same collision + if (multBin < 0 && multBin > SEhistos.size()) + LOGF(fatal, "multBin value passed to the mixTracks function is less than 0 or exceeds the configured number of Cent. bins"); + + for (int ii = 0; ii < tracks.size(); ii++) { // nested loop for all the combinations + for (int iii = ii + 1; iii < tracks.size(); iii++) { + + Pair->SetPair(tracks[ii], tracks[iii]); + float pair_kT = Pair->GetKt(); + int kTbin = o2::aod::singletrackselector::getBinIndex(pair_kT, _kTbins); + if (kTbin < 0 && kTbin > SEhistos[multBin].size()) + LOGF(fatal, "kTbin value obtained for a pair is less than 0 or exceeds the configured number of kT bins"); + + if (!Pair->IsClosePair(_deta, _dphi, _radiusTPC)) { + kThistos[multBin][kTbin]->Fill(pair_kT); + SEhistos[multBin][kTbin]->Fill(Pair->GetKstar()); // close pair rejection and fillig the SE histo + } + Pair->ResetPair(); + } + } + } + + template + void mixTracks(Type const& tracks1, Type const& tracks2, int multBin) + { // last value: 0 -- SE; 1 -- ME + if (multBin < 0 && multBin > SEhistos.size()) + LOGF(fatal, "multBin value passed to the mixTracks function is less than 0 or exceeds the configured number of Cent. bins"); + + for (auto ii : tracks1) { + for (auto iii : tracks2) { + + Pair->SetPair(ii, iii); + float pair_kT = Pair->GetKt(); + int kTbin = o2::aod::singletrackselector::getBinIndex(pair_kT, _kTbins); + if (kTbin < 0 && kTbin > SEhistos[multBin].size()) + LOGF(fatal, "kTbin value obtained for a pair is less than 0 or exceeds the configured number of kT bins"); + + if (!Pair->IsClosePair(_deta, _dphi, _radiusTPC)) { + if (!SE_or_ME) { + SEhistos[multBin][kTbin]->Fill(Pair->GetKstar()); + kThistos[multBin][kTbin]->Fill(pair_kT); + } else { + MEhistos[multBin][kTbin]->Fill(Pair->GetKstar()); + } + } + Pair->ResetPair(); + } + } + } + + void process(soa::Filtered const& collisions, soa::Filtered const& tracks) + { + if (_particlePDG_1 == 0 || _particlePDG_2 == 0) + LOGF(fatal, "One of passed PDG is 0!!!"); + + for (auto track : tracks) { + if (abs(track.singleCollSel().posZ()) > _vertexZ) + continue; + if (track.tpcNClsShared() > _tpcNClsShared || track.itsNCls() < _itsNCls) + continue; + if (track.singleCollSel().multPerc() < *_centBins.value.begin() || track.singleCollSel().multPerc() > *(_centBins.value.end() - 1)) + continue; + + if (track.sign() == _sign_1 && (track.p() < _PIDtrshld_1 ? o2::aod::singletrackselector::TPCselection(track, TPCcuts_1) : o2::aod::singletrackselector::TOFselection(track, TOFcuts_1))) { // filling the map: eventID <-> selected particles1 + selectedtracks_1[track.singleCollSelId()].push_back(std::make_shared(track)); + + registry.fill(HIST("p_first"), track.p()); + if (_particlePDG_1 == 211) { + registry.fill(HIST("nsigmaTOF_first"), track.p(), track.tofNSigmaPi()); + registry.fill(HIST("nsigmaTPC_first"), track.p(), track.tpcNSigmaPi()); + } + if (_particlePDG_1 == 321) { + registry.fill(HIST("nsigmaTOF_first"), track.p(), track.tofNSigmaKa()); + registry.fill(HIST("nsigmaTPC_first"), track.p(), track.tpcNSigmaKa()); + } + if (_particlePDG_1 == 2212) { + registry.fill(HIST("nsigmaTOF_first"), track.p(), track.tofNSigmaPr()); + registry.fill(HIST("nsigmaTPC_first"), track.p(), track.tpcNSigmaPr()); + } + if (_particlePDG_1 == 1000010020) { + registry.fill(HIST("nsigmaTOF_first"), track.p(), track.tofNSigmaDe()); + registry.fill(HIST("nsigmaTPC_first"), track.p(), track.tpcNSigmaDe()); + } + } + + if (IsIdentical) { + continue; + } else if (track.sign() != _sign_2 && !TOFselection(track, std::make_pair(_particlePDGtoReject, _rejectWithinNsigmaTOF)) && (track.p() < _PIDtrshld_2 ? o2::aod::singletrackselector::TPCselection(track, TPCcuts_2) : o2::aod::singletrackselector::TOFselection(track, TOFcuts_2))) { // filling the map: eventID <-> selected particles2 if (see condition above ^) + selectedtracks_2[track.singleCollSelId()].push_back(std::make_shared(track)); + + registry.fill(HIST("p_second"), track.p()); + if (_particlePDG_2 == 211) { + registry.fill(HIST("nsigmaTOF_second"), track.p(), track.tofNSigmaPi()); + registry.fill(HIST("nsigmaTPC_second"), track.p(), track.tpcNSigmaPi()); + } + if (_particlePDG_2 == 321) { + registry.fill(HIST("nsigmaTOF_second"), track.p(), track.tofNSigmaKa()); + registry.fill(HIST("nsigmaTPC_second"), track.p(), track.tpcNSigmaKa()); + } + if (_particlePDG_2 == 2212) { + registry.fill(HIST("nsigmaTOF_second"), track.p(), track.tofNSigmaPr()); + registry.fill(HIST("nsigmaTPC_second"), track.p(), track.tpcNSigmaPr()); + } + if (_particlePDG_2 == 1000010020) { + registry.fill(HIST("nsigmaTOF_second"), track.p(), track.tofNSigmaDe()); + registry.fill(HIST("nsigmaTPC_second"), track.p(), track.tpcNSigmaDe()); + } + } + } + + for (auto collision : collisions) { + if (collision.multPerc() < *_centBins.value.begin() || collision.multPerc() > *(_centBins.value.end() - 1)) + continue; + + if (selectedtracks_1.find(collision.globalIndex()) == selectedtracks_1.end()) { + if (IsIdentical) + continue; + else if (selectedtracks_2.find(collision.globalIndex()) == selectedtracks_2.end()) + continue; + } + int vertexBinToMix = round(collision.posZ() / (2 * _vertexZ / _vertexNbinsToMix)); + float centBinToMix = o2::aod::singletrackselector::getBinIndex(collision.multPerc(), _centBins, _multNsubBins); + + mixbins[std::pair{vertexBinToMix, centBinToMix}].push_back(std::make_shared(collision)); + } + + //====================================== mixing starts here ====================================== + + if (IsIdentical) { //====================================== mixing identical ====================================== + + for (auto i = mixbins.begin(); i != mixbins.end(); i++) { // iterating over all vertex&mult bins + + for (int indx1 = 0; indx1 < (i->second).size(); indx1++) { // loop over all the events in each vertex&mult bin + + auto col1 = (i->second)[indx1]; + + Pair->SetMagField1(col1->magField()); + Pair->SetMagField2(col1->magField()); + + int centBin = std::floor((i->first).second); + MultHistos[centBin]->Fill(col1->mult()); + + mixTracks(selectedtracks_1[col1->index()], centBin); // mixing SE identical + + for (int indx2 = indx1 + 1; indx2 < (i->second).size(); indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin + + auto col2 = (i->second)[indx2]; + + Pair->SetMagField2(col2->magField()); + mixTracks<1>(selectedtracks_1[col1->index()], selectedtracks_1[col2->index()], centBin); // mixing ME identical, in <> brackets: 0 -- SE; 1 -- ME + } + } + } + + } else { //====================================== mixing non-identical ====================================== + + for (auto i = mixbins.begin(); i != mixbins.end(); i++) { // iterating over all vertex&mult bins + + for (int indx1 = 0; indx1 < (i->second).size(); indx1++) { // loop over all the events in each vertex&mult bin + + auto col1 = (i->second)[indx1]; + + Pair->SetMagField1(col1->magField()); + Pair->SetMagField2(col1->magField()); + + int centBin = std::floor((i->first).second); + MultHistos[centBin]->Fill(col1->mult()); + + mixTracks<0>(selectedtracks_1[col1->index()], selectedtracks_2[col1->index()], centBin); // mixing SE non-identical, in <> brackets: 0 -- SE; 1 -- ME + + for (int indx2 = indx1 + 1; indx2 < (i->second).size(); indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin + + auto col2 = (i->second)[indx2]; + + Pair->SetMagField2(col2->magField()); + mixTracks<1>(selectedtracks_1[col1->index()], selectedtracks_2[col2->index()], centBin); // mixing ME non-identical, in <> brackets: 0 -- SE; 1 -- ME + } + } + } + + } //====================================== end of mixing non-identical ====================================== + + // clearing up + for (auto i = selectedtracks_1.begin(); i != selectedtracks_1.end(); i++) + (i->second).clear(); + selectedtracks_1.clear(); + + if (!IsIdentical) { + for (auto i = selectedtracks_2.begin(); i != selectedtracks_2.end(); i++) + (i->second).clear(); + selectedtracks_2.clear(); + } + + for (auto i = mixbins.begin(); i != mixbins.end(); i++) + (i->second).clear(); + mixbins.clear(); + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} diff --git a/PWGCF/Femto3D/Tasks/femto3dPairTaskMC.cxx b/PWGCF/Femto3D/Tasks/femto3dPairTaskMC.cxx new file mode 100755 index 00000000000..882a1e2221f --- /dev/null +++ b/PWGCF/Femto3D/Tasks/femto3dPairTaskMC.cxx @@ -0,0 +1,375 @@ +// 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. +// +/// \brief Femto3D pair mixing task +/// \author Sofia Tomassini, Gleb Romanenko, Nicolò Jacazio +/// \since 31 May 2023 + +#include +#include +#include + +#include "PWGCF/Femto3D/Core/femto3dPairTask.h" +#include "PWGCF/Femto3D/DataModel/singletrackselector.h" +#include "TLorentzVector.h" + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/ASoA.h" +#include "Framework/DataTypes.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/Expressions.h" +#include "Framework/StaticFor.h" +#include "MathUtils/Utils.h" +#include "Common/DataModel/Multiplicity.h" + +using namespace o2; +using namespace o2::soa; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct FemtoCorrelationsMC { + // using allinfo = soa::Join; // aod::pidTPCPr + /// Construct a registry object with direct declaration + HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; + + Configurable _min_P{"min_P", 0.0, "lower mometum limit"}; + Configurable _max_P{"max_P", 100.0, "upper mometum limit"}; + Configurable _eta{"eta", 100.0, "abs eta value limit"}; + Configurable _dcaXY{"dcaXY", 10.0, "abs dcaXY value limit"}; + Configurable _dcaZ{"dcaZ", 10.0, "abs dcaZ value limit"}; + Configurable _tpcNClsFound{"minTpcNClsFound", 0, "minimum allowed number of TPC clasters"}; + Configurable _tpcChi2NCl{"tpcChi2NCl", 100.0, "upper limit for chi2 value of a fit over TPC clasters"}; + Configurable _tpcCrossedRowsOverFindableCls{"tpcCrossedRowsOverFindableCls", 0, "lower limit of TPC CrossedRows/FindableCls value"}; + Configurable _tpcNClsShared{"maxTpcNClsShared", 100, "maximum allowed number of TPC shared clasters"}; + Configurable _itsNCls{"minItsNCls", 0, "minimum allowed number of ITS clasters"}; + Configurable _itsChi2NCl{"itsChi2NCl", 100.0, "upper limit for chi2 value of a fit over ITS clasters"}; + Configurable _vertexZ{"VertexZ", 10.0, "abs vertexZ value limit"}; + + Configurable _sign_1{"sign_1", 1, "sign of the first particle in a pair"}; + Configurable _particlePDG_1{"particlePDG_1", 2212, "PDG code of the first particle in a pair to perform PID for (only proton and deurton are supported now)"}; + Configurable> _tpcNSigma_1{"tpcNSigma_1", std::vector{-3.0f, 3.0f}, "first particle PID: Nsigma range in TPC before the TOF is used"}; + Configurable _PIDtrshld_1{"PIDtrshld_1", 10.0, "first particle PID: value of momentum from which the PID is done with TOF (before that only TPC is used)"}; + Configurable> _tofNSigma_1{"tofNSigma_1", std::vector{-3.0f, 3.0f}, "first particle PID: Nsigma range in TOF"}; + + Configurable _sign_2{"sign_2", 1, "sign of the second particle in a pair"}; + Configurable _particlePDG_2{"particlePDG_2", 2212, "PDG code of the second particle in a pair to perform PID for (only proton and deurton are supported now)"}; + Configurable> _tpcNSigma_2{"tpcNSigma_2", std::vector{-3.0f, 3.0f}, "second particle PID: Nsigma range in TPC before the TOF is used"}; + Configurable _PIDtrshld_2{"PIDtrshld_2", 10.0, "second particle PID: value of momentum from which the PID is done with TOF (before that only TPC is used)"}; + Configurable> _tofNSigma_2{"tofNSigma_2", std::vector{-3.0f, 3.0f}, "second particle PID: Nsigma range in TOF"}; + + Configurable _particlePDGtoReject{"particlePDGtoRejectFromSecond", 0, "applied only if the particles are non-identical and only to the second particle in the pair!!!"}; + Configurable> _rejectWithinNsigmaTOF{"rejectWithinNsigmaTOF", std::vector{-0.0f, 0.0f}, "TOF rejection Nsigma range for the particle specified with PDG to be rejected"}; + + Configurable _radiusTPC{"radiusTPC", 1.2, "TPC radius to calculate phi_star for"}; + + ConfigurableAxis CFkStarBinning{"CFkStarBinning", {500, 0.005, 5.005}, "k* binning of the res. matrix (Nbins, lowlimit, uplimit)"}; + + bool IsIdentical; + + std::pair> TPCcuts_1; + std::pair> TOFcuts_1; + + std::pair> TPCcuts_2; + std::pair> TOFcuts_2; + + using FilteredCollisions = aod::SingleCollSels; + using FilteredTracks = soa::Join; + + typedef std::shared_ptr::iterator> trkType; + typedef std::shared_ptr::iterator> colType; + + std::map> selectedtracks_1; + std::map> selectedtracks_2; + + std::unique_ptr> Pair = std::make_unique>(); + + Filter pFilter = o2::aod::singletrackselector::p > _min_P&& o2::aod::singletrackselector::p < _max_P; + Filter etaFilter = nabs(o2::aod::singletrackselector::eta) < _eta; + + Filter tpcTrkFilter = o2::aod::singletrackselector::tpcNClsFound >= _tpcNClsFound && + o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedTpcChi2NCl) < _tpcChi2NCl && + o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedTpcCrossedRowsOverFindableCls) > _tpcCrossedRowsOverFindableCls; + + Filter itsTrkFilter = o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedItsChi2NCl) < _itsChi2NCl; + + Filter vertexFilter = nabs(o2::aod::singletrackselector::posZ) < _vertexZ; + + void init(o2::framework::InitContext&) + { + + IsIdentical = (_sign_1 * _particlePDG_1 == _sign_2 * _particlePDG_2); + + Pair->SetIdentical(IsIdentical); + Pair->SetPDG1(_particlePDG_1); + Pair->SetPDG2(_particlePDG_2); + + TPCcuts_1 = std::make_pair(_particlePDG_1, _tpcNSigma_1); + TOFcuts_1 = std::make_pair(_particlePDG_1, _tofNSigma_1); + TPCcuts_2 = std::make_pair(_particlePDG_2, _tpcNSigma_2); + TOFcuts_2 = std::make_pair(_particlePDG_2, _tofNSigma_2); + + registry.add("FirstParticle/dcaxyz_vs_pt_primary", "dcaxyz_vs_pt_primary", kTH3F, {{100, 0., 5., "pt"}, {250, -1., 1., "DCA_XY(pt) primary"}, {250, -1., 1., "DCA_Z(pt) primary"}}); + registry.add("FirstParticle/dcaxyz_vs_pt_weakdecay", "dcaxyz_vs_pt_weakdecay", kTH3F, {{100, 0., 5., "pt"}, {250, -1., 1., "DCA_XY(pt) weakdecay"}, {250, -1., 1., "DCA_Z(pt) weakdecay"}}); + registry.add("FirstParticle/dcaxyz_vs_pt_material", "dcaxyz_vs_pt_material", kTH3F, {{100, 0., 5., "pt"}, {250, -1., 1., "DCA_XY(pt) material"}, {250, -1., 1., "DCA_Z(pt) material"}}); + + registry.add("FirstParticle/3dmomentumEl", "3dmomentumEl", kTH3F, {{100, 0., 5., "px"}, {100, 0., 5., "py"}, {100, 0., 5., "pz"}}); + registry.add("FirstParticle/3dmomentumMu", "3dmomentumMu", kTH3F, {{100, 0., 5., "px"}, {100, 0., 5., "py"}, {100, 0., 5., "pz"}}); + registry.add("FirstParticle/3dmomentumPi", "3dmomentumPi", kTH3F, {{100, 0., 5., "px"}, {100, 0., 5., "py"}, {100, 0., 5., "pz"}}); + registry.add("FirstParticle/3dmomentumKa", "3dmomentumKa", kTH3F, {{100, 0., 5., "px"}, {100, 0., 5., "py"}, {100, 0., 5., "pz"}}); + registry.add("FirstParticle/3dmomentumPr", "3dmomentumPr", kTH3F, {{100, 0., 5., "px"}, {100, 0., 5., "py"}, {100, 0., 5., "pz"}}); + registry.add("FirstParticle/3dmomentumDe", "3dmomentumDe", kTH3F, {{100, 0., 5., "px"}, {100, 0., 5., "py"}, {100, 0., 5., "pz"}}); + registry.add("FirstParticle/3dmomentumAll", "3dmomentumAll", kTH3F, {{100, 0., 5., "px"}, {100, 0., 5., "py"}, {100, 0., 5., "pz"}}); + registry.add("FirstParticle/PtSpectraEl", "PtSpectraEl", kTH1F, {{100, 0., 5.}}); + registry.add("FirstParticle/PtSpectraMu", "PtSpectraMu", kTH1F, {{100, 0., 5.}}); + registry.add("FirstParticle/PtSpectraPi", "PtSpectraPi", kTH1F, {{100, 0., 5.}}); + registry.add("FirstParticle/PtSpectraKa", "PtSpectraKa", kTH1F, {{100, 0., 5.}}); + registry.add("FirstParticle/PtSpectraPr", "PtSpectraPr", kTH1F, {{100, 0., 5.}}); + registry.add("FirstParticle/PtSpectraDe", "PtSpectraDe", kTH1F, {{100, 0., 5.}}); + registry.add("FirstParticle/PtSpectraAll", "PtSpectrAll", kTH1F, {{100, 0., 5.}}); + + if (!IsIdentical) { + registry.add("SecondParticle/dcaxyz_vs_pt_primary", "dcaxyz_vs_pt_primary", kTH3F, {{100, 0., 5., "pt"}, {200, -1., 1., "DCA_XY(pt) primary"}, {200, -1., 1., "DCA_Z(pt) primary"}}); + registry.add("SecondParticle/dcaxyz_vs_pt_weakdecay", "dcaxyz_vs_pt_weakdecay", kTH3F, {{100, 0., 5., "pt"}, {200, -1., 1., "DCA_XY(pt) weakdecay"}, {200, -1., 1., "DCA_Z(pt) weakdecay"}}); + registry.add("SecondParticle/dcaxyz_vs_pt_material", "dcaxyz_vs_pt_material", kTH3F, {{100, 0., 5., "pt"}, {200, -1., 1., "DCA_XY(pt) material"}, {200, -1., 1., "DCA_Z(pt) material"}}); + + registry.add("SecondParticle/3dmomentumEl", "3dmomentumEl", kTH3F, {{100, 0., 5., "px"}, {100, 0., 5., "py"}, {100, 0., 5., "pz"}}); + registry.add("SecondParticle/3dmomentumMu", "3dmomentumMu", kTH3F, {{100, 0., 5., "px"}, {100, 0., 5., "py"}, {100, 0., 5., "pz"}}); + registry.add("SecondParticle/3dmomentumPi", "3dmomentumPi", kTH3F, {{100, 0., 5., "px"}, {100, 0., 5., "py"}, {100, 0., 5., "pz"}}); + registry.add("SecondParticle/3dmomentumKa", "3dmomentumKa", kTH3F, {{100, 0., 5., "px"}, {100, 0., 5., "py"}, {100, 0., 5., "pz"}}); + registry.add("SecondParticle/3dmomentumPr", "3dmomentumPr", kTH3F, {{100, 0., 5., "px"}, {100, 0., 5., "py"}, {100, 0., 5., "pz"}}); + registry.add("SecondParticle/3dmomentumDe", "3dmomentumDe", kTH3F, {{100, 0., 5., "px"}, {100, 0., 5., "py"}, {100, 0., 5., "pz"}}); + registry.add("SecondParticle/3dmomentumAll", "3dmomentumAll", kTH3F, {{100, 0., 5., "px"}, {100, 0., 5., "py"}, {100, 0., 5., "pz"}}); + registry.add("SecondParticle/PtSpectraEl", "PtSpectraEl", kTH1F, {{100, 0., 5.}}); + registry.add("SecondParticle/PtSpectraMu", "PtSpectraMu", kTH1F, {{100, 0., 5.}}); + registry.add("SecondParticle/PtSpectraPi", "PtSpectraPi", kTH1F, {{100, 0., 5.}}); + registry.add("SecondParticle/PtSpectraKa", "PtSpectraKa", kTH1F, {{100, 0., 5.}}); + registry.add("SecondParticle/PtSpectraPr", "PtSpectraPr", kTH1F, {{100, 0., 5.}}); + registry.add("SecondParticle/PtSpectraDe", "PtSpectraDe", kTH1F, {{100, 0., 5.}}); + registry.add("SecondParticle/PtSpectraAll", "PtSpectrAll", kTH1F, {{100, 0., 5.}}); + } + + registry.add("ResolutionMatrix", "ResolutionMatrix_rec(gen)", kTH2F, {{CFkStarBinning, "k*_gen (GeV/c)"}, {CFkStarBinning, "k*_rec (GeV/c)"}}); + registry.add("DoubleTrackEffects", "DoubleTrackEffects_deta(dphi*)", kTH2F, {{200, -M_PI, M_PI, "dphi*"}, {200, -0.5, 0.5, "deta"}}); + } + + template + void fillEtaPhi(Type const& tracks) + { // template for particles from the same collision identical + for (int ii = 0; ii < tracks.size(); ii++) { // nested loop for all the combinations + for (int iii = ii + 1; iii < tracks.size(); iii++) { + + Pair->SetPair(tracks[ii], tracks[iii]); + Pair->SetMagField1((tracks[ii]->singleCollSel()).magField()); + Pair->SetMagField2((tracks[iii]->singleCollSel()).magField()); + + registry.fill(HIST("DoubleTrackEffects"), Pair->GetPhiStarDiff(_radiusTPC), Pair->GetEtaDiff()); + Pair->ResetPair(); + } + } + } + + template + void fillEtaPhi(Type const& tracks1, Type const& tracks2) + { // template for particles from the same collision non-identical + for (auto ii : tracks1) { + for (auto iii : tracks2) { + + Pair->SetPair(ii, iii); + Pair->SetMagField1((ii->singleCollSel()).magField()); + Pair->SetMagField2((iii->singleCollSel()).magField()); + + registry.fill(HIST("DoubleTrackEffects"), Pair->GetPhiStarDiff(_radiusTPC), Pair->GetEtaDiff()); + Pair->ResetPair(); + } + } + } + + template + void fillResMatrix(Type const& tracks1, Type const& tracks2) + { // template for ME + for (auto ii : tracks1) { + for (auto iii : tracks2) { + Pair->SetPair(ii, iii); + + TLorentzVector first4momentumGen; + first4momentumGen.SetPtEtaPhiM(ii->pt_MC(), ii->eta_MC(), ii->phi_MC(), particle_mass(_particlePDG_1)); + TLorentzVector second4momentumGen; + second4momentumGen.SetPtEtaPhiM(iii->pt_MC(), iii->eta_MC(), iii->phi_MC(), particle_mass(_particlePDG_2)); + + registry.fill(HIST("ResolutionMatrix"), o2::aod::singletrackselector::GetKstarFrom4vectors(first4momentumGen, second4momentumGen, IsIdentical), Pair->GetKstar()); + Pair->ResetPair(); + } + } + } + + void process(soa::Filtered const& collisions, soa::Filtered const& tracks) + { + if (_particlePDG_1 == 0 || _particlePDG_2 == 0) + LOGF(fatal, "One of passed PDG is 0!!!"); + + int trackPDG, trackOrigin; + + for (auto track : tracks) { + if (abs(track.singleCollSel().posZ()) > _vertexZ) + continue; + if (track.tpcNClsShared() > _tpcNClsShared || track.itsNCls() < _itsNCls) + continue; + + if (track.sign() == _sign_1 && (track.p() < _PIDtrshld_1 ? o2::aod::singletrackselector::TPCselection(track, TPCcuts_1) : o2::aod::singletrackselector::TOFselection(track, TOFcuts_1))) { + trackOrigin = track.origin(); + switch (trackOrigin) { + case 0: + registry.fill(HIST("FirstParticle/dcaxyz_vs_pt_primary"), track.pt(), track.dcaXY(), track.dcaZ()); + break; + case 1: + registry.fill(HIST("FirstParticle/dcaxyz_vs_pt_weakdecay"), track.pt(), track.dcaXY(), track.dcaZ()); + break; + case 2: + registry.fill(HIST("FirstParticle/dcaxyz_vs_pt_material"), track.pt(), track.dcaXY(), track.dcaZ()); + break; + } + if (abs(track.dcaXY()) > _dcaXY || abs(track.dcaZ()) > _dcaZ) + continue; + + selectedtracks_1[track.singleCollSelId()].push_back(std::make_shared(track)); // filling the map: eventID <-> selected particles1 + + trackPDG = abs(track.pdgCode()); + + registry.fill(HIST("FirstParticle/3dmomentumAll"), track.px_MC(), track.py_MC(), track.pz_MC()); + registry.fill(HIST("FirstParticle/PtSpectraAll"), track.pt_MC()); + + switch (trackPDG) { + case 11: + registry.fill(HIST("FirstParticle/3dmomentumEl"), track.px_MC(), track.py_MC(), track.pz_MC()); + registry.fill(HIST("FirstParticle/PtSpectraEl"), track.pt_MC()); + break; + case 13: + registry.fill(HIST("FirstParticle/3dmomentumMu"), track.px_MC(), track.py_MC(), track.pz_MC()); + registry.fill(HIST("FirstParticle/PtSpectraMu"), track.pt_MC()); + break; + case 211: + registry.fill(HIST("FirstParticle/3dmomentumPi"), track.px_MC(), track.py_MC(), track.pz_MC()); + registry.fill(HIST("FirstParticle/PtSpectraPi"), track.pt_MC()); + break; + case 321: + registry.fill(HIST("FirstParticle/3dmomentumKa"), track.px_MC(), track.py_MC(), track.pz_MC()); + registry.fill(HIST("FirstParticle/PtSpectraKa"), track.pt_MC()); + break; + case 2212: + registry.fill(HIST("FirstParticle/3dmomentumPr"), track.px_MC(), track.py_MC(), track.pz_MC()); + registry.fill(HIST("FirstParticle/PtSpectraPr"), track.pt_MC()); + break; + case 1000010020: + registry.fill(HIST("FirstParticle/3dmomentumDe"), track.px_MC(), track.py_MC(), track.pz_MC()); + registry.fill(HIST("FirstParticle/PtSpectraDe"), track.pt_MC()); + break; + } + } + + if (IsIdentical) { + continue; + } else if (track.sign() != _sign_2 && !TOFselection(track, std::make_pair(_particlePDGtoReject, _rejectWithinNsigmaTOF)) && (track.p() < _PIDtrshld_2 ? o2::aod::singletrackselector::TPCselection(track, TPCcuts_2) : o2::aod::singletrackselector::TOFselection(track, TOFcuts_2))) { // filling the map: eventID <-> selected particles2 if (see condition above ^) + trackOrigin = track.origin(); + switch (trackOrigin) { + case 0: + registry.fill(HIST("SecondParticle/dcaxyz_vs_pt_primary"), track.pt(), track.dcaXY(), track.dcaZ()); + break; + case 1: + registry.fill(HIST("SecondParticle/dcaxyz_vs_pt_weakdecay"), track.pt(), track.dcaXY(), track.dcaZ()); + break; + case 2: + registry.fill(HIST("SecondParticle/dcaxyz_vs_pt_material"), track.pt(), track.dcaXY(), track.dcaZ()); + break; + } + if (abs(track.dcaXY()) > _dcaXY || abs(track.dcaZ()) > _dcaZ) + continue; + + selectedtracks_2[track.singleCollSelId()].push_back(std::make_shared(track)); // filling the map: eventID <-> selected particles2 + + trackPDG = abs(track.pdgCode()); + + registry.fill(HIST("SecondParticle/3dmomentumAll"), track.px_MC(), track.py_MC(), track.pz_MC()); + registry.fill(HIST("SecondParticle/PtSpectraAll"), track.pt_MC()); + + switch (trackPDG) { + case 11: + registry.fill(HIST("SecondParticle/3dmomentumEl"), track.px_MC(), track.py_MC(), track.pz_MC()); + registry.fill(HIST("SecondParticle/PtSpectraEl"), track.pt_MC()); + break; + case 13: + registry.fill(HIST("SecondParticle/3dmomentumMu"), track.px_MC(), track.py_MC(), track.pz_MC()); + registry.fill(HIST("SecondParticle/PtSpectraMu"), track.pt_MC()); + break; + case 211: + registry.fill(HIST("SecondParticle/3dmomentumPi"), track.px_MC(), track.py_MC(), track.pz_MC()); + registry.fill(HIST("SecondParticle/PtSpectraPi"), track.pt_MC()); + break; + case 321: + registry.fill(HIST("SecondParticle/3dmomentumKa"), track.px_MC(), track.py_MC(), track.pz_MC()); + registry.fill(HIST("SecondParticle/PtSpectraKa"), track.pt_MC()); + break; + case 2212: + registry.fill(HIST("SecondParticle/3dmomentumPr"), track.px_MC(), track.py_MC(), track.pz_MC()); + registry.fill(HIST("SecondParticle/PtSpectraPr"), track.pt_MC()); + break; + case 1000010020: + registry.fill(HIST("SecondParticle/3dmomentumDe"), track.px_MC(), track.py_MC(), track.pz_MC()); + registry.fill(HIST("SecondParticle/PtSpectraDe"), track.pt_MC()); + break; + } + } + } + + //====================================== filling deta(dphi*) & res. matrix starts here ====================================== + + if (IsIdentical) { //====================================== identical ====================================== + + for (auto i = selectedtracks_1.begin(); i != selectedtracks_1.end(); i++) { // iterating over all selected collisions with selected tracks + fillEtaPhi(i->second); // filling deta(dphi*) -- SE identical + auto j = i; + + for (++j; j != selectedtracks_1.end(); j++) { // nested loop to do all the ME identical combinations + fillResMatrix(i->second, j->second); // filling res. matrix -- ME identical + } + } + } else { //====================================== non-identical ====================================== + + for (auto i = selectedtracks_1.begin(); i != selectedtracks_1.end(); i++) { // iterating over all selected collisions with selected tracks1 + auto ii = selectedtracks_2.find(i->first); + if (ii != selectedtracks_2.end()) + fillEtaPhi(i->second, ii->second); // checking if there are tracks2 for the choosen collision and filling deta(dphi*) -- SE non-identical + + for (auto j = selectedtracks_2.begin(); j != selectedtracks_2.end(); j++) { // nested loop to do all the ME non-identical combinations + fillResMatrix(i->second, j->second); // filling res. matrix -- ME non-identical + } + } + } //====================================== end of mixing non-identical ====================================== + + // clearing up + for (auto i = selectedtracks_1.begin(); i != selectedtracks_1.end(); i++) + (i->second).clear(); + selectedtracks_1.clear(); + + if (!IsIdentical) { + for (auto i = selectedtracks_2.begin(); i != selectedtracks_2.end(); i++) + (i->second).clear(); + selectedtracks_2.clear(); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} diff --git a/PWGCF/Femto3D/Tasks/femto3dQA.cxx b/PWGCF/Femto3D/Tasks/femto3dQA.cxx new file mode 100755 index 00000000000..d652d4d0d3f --- /dev/null +++ b/PWGCF/Femto3D/Tasks/femto3dQA.cxx @@ -0,0 +1,221 @@ +// 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. +// +/// \brief Femto3D QA task +/// \author Sofia Tomassini, Gleb Romanenko, Nicolò Jacazio +/// \since 31 May 2023 + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include +#include + +#include "Framework/ASoA.h" +#include "MathUtils/Utils.h" +#include "Framework/DataTypes.h" +#include "Common/DataModel/Multiplicity.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/Expressions.h" + +#include "Framework/StaticFor.h" +#include "PWGCF/Femto3D/DataModel/singletrackselector.h" + +using namespace o2; +using namespace o2::soa; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct QAHistograms { + // using allinfo = soa::Join; // aod::pidTPCPr + /// Construct a registry object with direct declaration + HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; + + Configurable _sign{"sign", 1, "sign of a track"}; + Configurable _vertexZ{"VertexZ", 10.0, "abs vertexZ value limit"}; + Configurable _min_P{"min_P", 0.0, "lower mometum limit"}; + Configurable _max_P{"max_P", 100.0, "upper mometum limit"}; + Configurable _eta{"eta", 100.0, "abs eta value limit"}; + Configurable _dcaXY{"dcaXY", 10.0, "abs dcaXY value limit"}; + Configurable _dcaZ{"dcaZ", 10.0, "abs dcaZ value limit"}; + Configurable _tpcNClsFound{"minTpcNClsFound", 0, "minimum allowed number of TPC clasters"}; + Configurable _tpcChi2NCl{"tpcChi2NCl", 100.0, "upper limit for chi2 value of a fit over TPC clasters"}; + Configurable _tpcCrossedRowsOverFindableCls{"tpcCrossedRowsOverFindableCls", 0, "lower limit of TPC CrossedRows/FindableCls value"}; + Configurable _tpcNClsShared{"maxTpcNClsShared", 0, "maximum allowed number of TPC shared clasters"}; + Configurable _itsNCls{"minItsNCls", 0, "minimum allowed number of ITS clasters for a track"}; + Configurable _itsChi2NCl{"itsChi2NCl", 100.0, "upper limit for chi2 value of a fit over ITS clasters for a track"}; + Configurable _particlePDG{"particlePDG", 2212, "PDG code of a particle to perform PID for (only proton and deurton are supported now)"}; + Configurable> _tpcNSigma{"tpcNSigma", std::vector{-4.0f, 4.0f}, "Nsigma range in TPC before the TOF is used"}; + Configurable _PIDtrshld{"PIDtrshld", 10.0, "value of momentum from which the PID is done with TOF (before that only TPC is used)"}; + Configurable> _tofNSigma{"tofNSigma", std::vector{-4.0f, 4.0f}, "Nsigma range in TOF"}; + + Configurable _particlePDGtoReject{"particlePDGtoReject", 211, "PDG codes of perticles that will be rejected with TOF (only proton and deurton are supported now)"}; + Configurable> _rejectWithinNsigmaTOF{"rejectWithinNsigmaTOF", std::vector{-0.0f, 0.0f}, "TOF rejection Nsigma range for particles specified with PDG to be rejected"}; + + std::pair> TPCcuts; + std::pair> TOFcuts; + + Filter signFilter = o2::aod::singletrackselector::sign == _sign; + Filter pFilter = o2::aod::singletrackselector::p > _min_P&& o2::aod::singletrackselector::p < _max_P; + Filter etaFilter = nabs(o2::aod::singletrackselector::eta) < _eta; + + Filter tpcTrkFilter = o2::aod::singletrackselector::tpcNClsFound >= _tpcNClsFound && + o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedTpcChi2NCl) < _tpcChi2NCl && + o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedTpcCrossedRowsOverFindableCls) > _tpcCrossedRowsOverFindableCls; + + Filter dcaFilter = nabs(o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedDcaXY)) < _dcaXY && + nabs(o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedDcaZ)) < _dcaZ; + + Filter itsTrkFilter = o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedItsChi2NCl) < _itsChi2NCl; + + Filter vertexFilter = nabs(o2::aod::singletrackselector::posZ) < _vertexZ; + + void init(o2::framework::InitContext&) + { + TPCcuts = std::make_pair(_particlePDG, _tpcNSigma); + TOFcuts = std::make_pair(_particlePDG, _tofNSigma); + + registry.add("TPCSignal_nocuts", "TPC signal without cuts", kTH2F, {{{200, 0., 5.0, "#it{p}_{inner} (GeV/#it{c})"}, {1000, 0., 1000.0, "dE/dx in TPC (arbitrary units)"}}}); + registry.add("TOFSignal_nocuts", "TOF signal without cuts", kTH2F, {{{200, 0., 5.0, "#it{p} (GeV/#it{c})"}, {100, 0., 1.5, "#beta"}}}); + + registry.add("eta", "eta", kTH1F, {{200, -2.5, 2.5, "eta"}}); + registry.add("phi", "phi", kTH1F, {{200, 0., 2. * M_PI, "phi"}}); + registry.add("px", "px", kTH1F, {{100, 0., 5., "px"}}); + registry.add("py", "py", kTH1F, {{100, 0., 5., "py"}}); + registry.add("pz", "pz", kTH1F, {{100, 0., 5., "pz"}}); + registry.add("p", "p", kTH1F, {{100, 0., 5., "p"}}); + registry.add("pt", "pt", kTH1F, {{100, 0., 5., "pt"}}); + registry.add("sign", "sign", kTH1F, {{3, -1.5, 1.5, "sign"}}); + registry.add("dcaxy_to_p", "dcaxy_to_p", kTH2F, {{100, 0., 5.0, "p"}, {200, -1., 1., "dcaxy"}}); + registry.add("dcaxy_to_pt", "dcaxy_to_pt", kTH2F, {{100, 0., 5., "pt"}, {200, -1., 1., "dcaxy"}}); + registry.add("dcaz_to_p", "dcaz_to_p", kTH2F, {{100, 0., 5., "p"}, {200, -1., 1., "dcaz"}}); + registry.add("dcaz_to_pt", "dcaz_to_pt", kTH2F, {{100, 0., 5., "pt"}, {200, -1., 1., "dcaz"}}); + registry.add("TPCClusters", "TPCClusters", kTH1F, {{163, -0.5, 162.5, "NTPCClust"}}); + registry.add("TPCCrossedRowsOverFindableCls", "TPCCrossedRowsOverFindableCls", kTH1F, {{100, 0.0, 10.0, "NcrossedRowsOverFindable"}}); + registry.add("TPCNClsShared", "TPCNClsShared", kTH1F, {{160, -0.5, 159.5, "TPCNshared"}}); + registry.add("ITSClusters", "ITSClusters", kTH1F, {{10, -0.5, 9.5, "NITSClust"}}); + registry.add("ITSchi2", "ITSchi2", kTH1F, {{100, 0.0, 40., "ITSchi2"}}); + registry.add("TPCchi2", "TPCchi2", kTH1F, {{100, 0.0, 6., "TPCchi2"}}); + + registry.add("TPCSignal", "TPC Signal", kTH2F, {{{200, 0., 5.0, "#it{p}_{inner} (GeV/#it{c})"}, {1000, 0., 1000.0, "dE/dx in TPC (arbitrary units)"}}}); + registry.add("TOFSignal", "TOF Signal", kTH2F, {{200, 0., 5.0, "#it{p} (GeV/#it{c})"}, {100, 0., 1.5, "#beta"}}); + + switch (_particlePDG) { + case 211: + registry.add("nsigmaTOFPi", "nsigmaTOFPi", kTH2F, {{100, 0., 5.}, {100, -10., 10.}}); + registry.add("nsigmaTPCPi", "nsigmaTPCPi", kTH2F, {{100, 0., 5.}, {100, -10., 10.}}); + break; + case 321: + registry.add("nsigmaTOFKa", "nsigmaTOFKa", kTH2F, {{100, 0., 5.}, {100, -10., 10.}}); + registry.add("nsigmaTPCKa", "nsigmaTPCKa", kTH2F, {{100, 0., 5.}, {100, -10., 10.}}); + break; + case 2212: + registry.add("nsigmaTOFPr", "nsigmaTOFPr", kTH2F, {{100, 0., 5.}, {100, -10., 10.}}); + registry.add("nsigmaTPCPr", "nsigmaTPCPr", kTH2F, {{100, 0., 5.}, {100, -10., 10.}}); + break; + case 1000010020: + registry.add("nsigmaTOFDe", "nsigmaTOFDe", kTH2F, {{100, 0., 5.}, {100, -10., 10.}}); + registry.add("nsigmaTPCDe", "nsigmaTPCDe", kTH2F, {{100, 0., 5.}, {100, -10., 10.}}); + break; + default: + break; + } + + registry.add("posZ", "posZ", kTH1F, {{300, -16., 16., "posZ"}}); + registry.add("mult", "mult", kTH1F, {{500, 0., 500., "mult"}}); + } + + template + void fillHistograms(ColsType const& collisions, TracksType const& tracks) + { + for (auto& collision : collisions) { + registry.fill(HIST("posZ"), collision.posZ()); + registry.fill(HIST("mult"), collision.mult()); + } + + for (auto& track : tracks) { + if (abs(track.singleCollSel().posZ()) > _vertexZ) + continue; + if ((track.tpcNClsShared()) > _tpcNClsShared || (track.itsNCls()) < _itsNCls) + continue; + + if constexpr (FillExtra) { + registry.fill(HIST("TPCSignal_nocuts"), track.tpcInnerParam(), track.tpcSignal()); + registry.fill(HIST("TOFSignal_nocuts"), track.p(), track.beta()); + } + + if (!TOFselection(track, std::make_pair(_particlePDGtoReject, _rejectWithinNsigmaTOF)) && (track.p() < _PIDtrshld ? o2::aod::singletrackselector::TPCselection(track, TPCcuts) : o2::aod::singletrackselector::TOFselection(track, TOFcuts))) { + registry.fill(HIST("eta"), track.eta()); + registry.fill(HIST("phi"), track.phi()); + registry.fill(HIST("px"), track.px()); + registry.fill(HIST("py"), track.py()); + registry.fill(HIST("pz"), track.pz()); + registry.fill(HIST("p"), track.p()); + registry.fill(HIST("pt"), track.pt()); + registry.fill(HIST("sign"), track.sign()); + registry.fill(HIST("dcaxy_to_p"), track.p(), track.dcaXY()); + registry.fill(HIST("dcaxy_to_pt"), track.pt(), track.dcaXY()); + registry.fill(HIST("dcaz_to_p"), track.p(), track.dcaZ()); + registry.fill(HIST("dcaz_to_pt"), track.pt(), track.dcaZ()); + registry.fill(HIST("TPCClusters"), track.tpcNClsFound()); + registry.fill(HIST("TPCCrossedRowsOverFindableCls"), track.tpcCrossedRowsOverFindableCls()); + registry.fill(HIST("TPCNClsShared"), track.tpcNClsShared()); + registry.fill(HIST("ITSClusters"), track.itsNCls()); + registry.fill(HIST("ITSchi2"), track.itsChi2NCl()); + registry.fill(HIST("TPCchi2"), track.tpcChi2NCl()); + + switch (_particlePDG) { + case 211: + registry.fill(HIST("nsigmaTOFPi"), track.p(), track.tofNSigmaPi()); + registry.fill(HIST("nsigmaTPCPi"), track.p(), track.tpcNSigmaPi()); + break; + case 321: + registry.fill(HIST("nsigmaTOFKa"), track.p(), track.tofNSigmaKa()); + registry.fill(HIST("nsigmaTPCKa"), track.p(), track.tpcNSigmaKa()); + break; + case 2212: + registry.fill(HIST("nsigmaTOFPr"), track.p(), track.tofNSigmaPr()); + registry.fill(HIST("nsigmaTPCPr"), track.p(), track.tpcNSigmaPr()); + break; + case 1000010020: + registry.fill(HIST("nsigmaTOFDe"), track.p(), track.tofNSigmaDe()); + registry.fill(HIST("nsigmaTPCDe"), track.p(), track.tpcNSigmaDe()); + break; + default: + break; + } + + if constexpr (FillExtra) { + registry.fill(HIST("TPCSignal"), track.tpcInnerParam(), track.tpcSignal()); + registry.fill(HIST("TOFSignal"), track.p(), track.beta()); + } + } + } + } + + void processDefault(soa::Filtered const& collisions, soa::Filtered const& tracks) + { + fillHistograms(collisions, tracks); + } + PROCESS_SWITCH(QAHistograms, processDefault, "process default", true); + + void processExtra(soa::Filtered const& collisions, soa::Filtered> const& tracks) + { + fillHistograms(collisions, tracks); + } + PROCESS_SWITCH(QAHistograms, processExtra, "process extra", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}