diff --git a/PWGLF/DataModel/Reduced3BodyTables.h b/PWGLF/DataModel/Reduced3BodyTables.h new file mode 100644 index 00000000000..a11081cd848 --- /dev/null +++ b/PWGLF/DataModel/Reduced3BodyTables.h @@ -0,0 +1,269 @@ +// 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. + +#ifndef PWGLF_DATAMODEL_REDUCED3BODYTABLES_H_ +#define PWGLF_DATAMODEL_REDUCED3BODYTABLES_H_ + +#include +#include "Framework/AnalysisDataModel.h" +#include "Common/Core/RecoDecay.h" +#include "CommonConstants/PhysicsConstants.h" +#include "PWGLF/DataModel/Vtx3BodyTables.h" + +namespace o2::aod +{ + +DECLARE_SOA_TABLE(ReducedCollisions, "AOD", "REDUCEDCOLLISION", //! reduced collision table (same structure as the original collision table) + o2::soa::Index<>, collision::BCId, + collision::PosX, collision::PosY, collision::PosZ, + collision::CovXX, collision::CovXY, collision::CovYY, collision::CovXZ, collision::CovYZ, collision::CovZZ, + collision::Flags, collision::Chi2, collision::NumContrib, + collision::CollisionTime, collision::CollisionTimeRes, + bc::RunNumber); + +DECLARE_SOA_TABLE(ReducedPVMults, "AOD", "REDUCEDPVMULT", //! Multiplicity from the PV contributors, joinable with reducedCollisions + mult::MultNTracksPV); + +DECLARE_SOA_TABLE(ReducedCentFT0Cs, "AOD", "REDUCEDCENTFT0C", //! Reduced Run 3 FT0C centrality table, joinable with reducedCollisions + cent::CentFT0C); + +namespace reducedtracks3body +{ +// track parameter definition +DECLARE_SOA_INDEX_COLUMN(Collision, collision); //! Collision to which this track belongs +DECLARE_SOA_COLUMN(X, x, float); //! +DECLARE_SOA_COLUMN(Alpha, alpha, float); //! +DECLARE_SOA_COLUMN(Y, y, float); //! +DECLARE_SOA_COLUMN(Z, z, float); //! +DECLARE_SOA_COLUMN(Snp, snp, float); //! +DECLARE_SOA_COLUMN(Tgl, tgl, float); //! +DECLARE_SOA_COLUMN(Signed1Pt, signed1Pt, float); //! (sign of charge)/Pt in c/GeV. Use pt() and sign() instead +DECLARE_SOA_EXPRESSION_COLUMN(Phi, phi, float, //! Phi of the track, in radians within [0, 2pi) + ifnode(nasin(aod::track::snp) + aod::track::alpha < 0.0f, nasin(aod::track::snp) + aod::track::alpha + o2::constants::math::TwoPI, + ifnode(nasin(aod::track::snp) + aod::track::alpha >= o2::constants::math::TwoPI, nasin(aod::track::snp) + aod::track::alpha - o2::constants::math::TwoPI, + nasin(aod::track::snp) + aod::track::alpha))); +DECLARE_SOA_EXPRESSION_COLUMN(Eta, eta, float, //! Pseudorapidity + -1.f * nlog(ntan(o2::constants::math::PIQuarter - 0.5f * natan(aod::track::tgl)))); +DECLARE_SOA_EXPRESSION_COLUMN(Pt, pt, float, //! Transverse momentum of the track in GeV/c + ifnode(nabs(aod::track::signed1Pt) <= o2::constants::math::Almost0, o2::constants::math::VeryBig, nabs(1.f / aod::track::signed1Pt))); +DECLARE_SOA_DYNAMIC_COLUMN(Sign, sign, //! Charge: positive: 1, negative: -1 + [](float signed1Pt) -> short { return (signed1Pt > 0) ? 1 : -1; }); +DECLARE_SOA_DYNAMIC_COLUMN(Px, px, //! Momentum in x-direction in GeV/c + [](float signed1Pt, float snp, float alpha) -> float { + auto pt = 1.f / std::abs(signed1Pt); + // FIXME: GCC & clang should optimize to sincosf + float cs = cosf(alpha), sn = sinf(alpha); + auto r = std::sqrt((1.f - snp) * (1.f + snp)); + return pt * (r * cs - snp * sn); + }); +DECLARE_SOA_DYNAMIC_COLUMN(Py, py, //! Momentum in y-direction in GeV/c + [](float signed1Pt, float snp, float alpha) -> float { + auto pt = 1.f / std::abs(signed1Pt); + // FIXME: GCC & clang should optimize to sincosf + float cs = cosf(alpha), sn = sinf(alpha); + auto r = std::sqrt((1.f - snp) * (1.f + snp)); + return pt * (snp * cs + r * sn); + }); +DECLARE_SOA_DYNAMIC_COLUMN(Pz, pz, //! Momentum in z-direction in GeV/c + [](float signed1Pt, float tgl) -> float { + auto pt = 1.f / std::abs(signed1Pt); + return pt * tgl; + }); +DECLARE_SOA_EXPRESSION_COLUMN(P, p, float, //! Momentum in Gev/c + ifnode(nabs(aod::track::signed1Pt) <= o2::constants::math::Almost0, o2::constants::math::VeryBig, 0.5f * (ntan(o2::constants::math::PIQuarter - 0.5f * natan(aod::track::tgl)) + 1.f / ntan(o2::constants::math::PIQuarter - 0.5f * natan(aod::track::tgl))) / nabs(aod::track::signed1Pt))); +DECLARE_SOA_DYNAMIC_COLUMN(Rapidity, rapidity, //! Track rapidity, computed under the mass assumption given as input + [](float signed1Pt, float tgl, float mass) -> float { + const auto pt = 1.f / std::abs(signed1Pt); + const auto pz = pt * tgl; + const auto p = 0.5f * (std::tan(o2::constants::math::PIQuarter - 0.5f * std::atan(tgl)) + 1.f / std::tan(o2::constants::math::PIQuarter - 0.5f * std::atan(tgl))) * pt; + const auto energy = std::sqrt(p * p + mass * mass); + return 0.5f * std::log((energy + pz) / (energy - pz)); + }); + +// tracks cov matrix parameter definition +DECLARE_SOA_COLUMN(SigmaY, sigmaY, float); //! Covariance matrix +DECLARE_SOA_COLUMN(SigmaZ, sigmaZ, float); //! Covariance matrix +DECLARE_SOA_COLUMN(SigmaSnp, sigmaSnp, float); //! Covariance matrix +DECLARE_SOA_COLUMN(SigmaTgl, sigmaTgl, float); //! Covariance matrix +DECLARE_SOA_COLUMN(Sigma1Pt, sigma1Pt, float); //! Covariance matrix +DECLARE_SOA_COLUMN(RhoZY, rhoZY, int8_t); //! Covariance matrix in compressed form +DECLARE_SOA_COLUMN(RhoSnpY, rhoSnpY, int8_t); //! Covariance matrix in compressed form +DECLARE_SOA_COLUMN(RhoSnpZ, rhoSnpZ, int8_t); //! Covariance matrix in compressed form +DECLARE_SOA_COLUMN(RhoTglY, rhoTglY, int8_t); //! Covariance matrix in compressed form +DECLARE_SOA_COLUMN(RhoTglZ, rhoTglZ, int8_t); //! Covariance matrix in compressed form +DECLARE_SOA_COLUMN(RhoTglSnp, rhoTglSnp, int8_t); //! Covariance matrix in compressed form +DECLARE_SOA_COLUMN(Rho1PtY, rho1PtY, int8_t); //! Covariance matrix in compressed form +DECLARE_SOA_COLUMN(Rho1PtZ, rho1PtZ, int8_t); //! Covariance matrix in compressed form +DECLARE_SOA_COLUMN(Rho1PtSnp, rho1PtSnp, int8_t); //! Covariance matrix in compressed form +DECLARE_SOA_COLUMN(Rho1PtTgl, rho1PtTgl, int8_t); //! Covariance matrix in compressed form + +DECLARE_SOA_EXPRESSION_COLUMN(CYY, cYY, float, //! + aod::track::sigmaY* aod::track::sigmaY); +DECLARE_SOA_EXPRESSION_COLUMN(CZY, cZY, float, //! + (aod::track::rhoZY / 128.f) * (aod::track::sigmaZ * aod::track::sigmaY)); +DECLARE_SOA_EXPRESSION_COLUMN(CZZ, cZZ, float, //! + aod::track::sigmaZ* aod::track::sigmaZ); +DECLARE_SOA_EXPRESSION_COLUMN(CSnpY, cSnpY, float, //! + (aod::track::rhoSnpY / 128.f) * (aod::track::sigmaSnp * aod::track::sigmaY)); +DECLARE_SOA_EXPRESSION_COLUMN(CSnpZ, cSnpZ, float, //! + (aod::track::rhoSnpZ / 128.f) * (aod::track::sigmaSnp * aod::track::sigmaZ)); +DECLARE_SOA_EXPRESSION_COLUMN(CSnpSnp, cSnpSnp, float, //! + aod::track::sigmaSnp* aod::track::sigmaSnp); +DECLARE_SOA_EXPRESSION_COLUMN(CTglY, cTglY, float, //! + (aod::track::rhoTglY / 128.f) * (aod::track::sigmaTgl * aod::track::sigmaY)); +DECLARE_SOA_EXPRESSION_COLUMN(CTglZ, cTglZ, float, //! + (aod::track::rhoTglZ / 128.f) * (aod::track::sigmaTgl * aod::track::sigmaZ)); +DECLARE_SOA_EXPRESSION_COLUMN(CTglSnp, cTglSnp, float, //! + (aod::track::rhoTglSnp / 128.f) * (aod::track::sigmaTgl * aod::track::sigmaSnp)); +DECLARE_SOA_EXPRESSION_COLUMN(CTglTgl, cTglTgl, float, //! + aod::track::sigmaTgl* aod::track::sigmaTgl); +DECLARE_SOA_EXPRESSION_COLUMN(C1PtY, c1PtY, float, //! + (aod::track::rho1PtY / 128.f) * (aod::track::sigma1Pt * aod::track::sigmaY)); +DECLARE_SOA_EXPRESSION_COLUMN(C1PtZ, c1PtZ, float, //! + (aod::track::rho1PtZ / 128.f) * (aod::track::sigma1Pt * aod::track::sigmaZ)); +DECLARE_SOA_EXPRESSION_COLUMN(C1PtSnp, c1PtSnp, float, //! + (aod::track::rho1PtSnp / 128.f) * (aod::track::sigma1Pt * aod::track::sigmaSnp)); +DECLARE_SOA_EXPRESSION_COLUMN(C1PtTgl, c1PtTgl, float, //! + (aod::track::rho1PtTgl / 128.f) * (aod::track::sigma1Pt * aod::track::sigmaTgl)); +DECLARE_SOA_EXPRESSION_COLUMN(C1Pt21Pt2, c1Pt21Pt2, float, //! + aod::track::sigma1Pt* aod::track::sigma1Pt); + +// tracks extra parameter definition +DECLARE_SOA_COLUMN(TPCInnerParam, tpcInnerParam, float); //! Momentum at inner wall of the TPC +DECLARE_SOA_COLUMN(Flags, flags, uint32_t); //! Track flags. Run 2: see TrackFlagsRun2Enum | Run 3: see TrackFlags +DECLARE_SOA_COLUMN(ITSClusterSizes, itsClusterSizes, uint32_t); //! Clusters sizes, four bits per a layer, starting from the innermost +DECLARE_SOA_COLUMN(TPCNClsFindable, tpcNClsFindable, uint8_t); //! Findable TPC clusters for this track geometry +DECLARE_SOA_COLUMN(TPCNClsFindableMinusFound, tpcNClsFindableMinusFound, int8_t); //! TPC Clusters: Findable - Found +DECLARE_SOA_COLUMN(TPCNClsFindableMinusCrossedRows, tpcNClsFindableMinusCrossedRows, int8_t); //! TPC Clusters: Findable - crossed rows +DECLARE_SOA_COLUMN(TRDPattern, trdPattern, uint8_t); //! Contributor to the track on TRD layer in bits 0-5, starting from the innermost, bit 6 indicates a potentially split tracklet, bit 7 if the track crossed a padrow +DECLARE_SOA_COLUMN(TPCChi2NCl, tpcChi2NCl, float); //! Chi2 / cluster for the TPC track segment +DECLARE_SOA_COLUMN(TOFChi2, tofChi2, float); //! Chi2 for the TOF track segment +DECLARE_SOA_COLUMN(TPCSignal, tpcSignal, float); //! dE/dx signal in the TPC +DECLARE_SOA_COLUMN(TOFExpMom, tofExpMom, float); //! TOF expected momentum obtained in tracking, used to compute the expected times + +DECLARE_SOA_EXPRESSION_COLUMN(DetectorMap, detectorMap, uint8_t, //! Detector map version 1, see enum DetectorMapEnum + ifnode(aod::track::itsClusterSizes > (uint32_t)0, static_cast(o2::aod::track::ITS), (uint8_t)0x0) | + ifnode(aod::track::tpcNClsFindable > (uint8_t)0, static_cast(o2::aod::track::TPC), (uint8_t)0x0) | + ifnode(aod::track::trdPattern > (uint8_t)0, static_cast(o2::aod::track::TRD), (uint8_t)0x0) | + ifnode((aod::track::tofChi2 >= 0.f) && (aod::track::tofExpMom > 0.f), static_cast(o2::aod::track::TOF), (uint8_t)0x0)); + +DECLARE_SOA_DYNAMIC_COLUMN(ITSClsSizeInLayer, itsClsSizeInLayer, //! Size of the ITS cluster in a given layer + [](uint32_t itsClusterSizes, int layer) -> uint8_t { + if (layer >= 7 || layer < 0) { + return 0; + } + return (itsClusterSizes >> (layer * 4)) & 0xf; + }); + +DECLARE_SOA_DYNAMIC_COLUMN(HasTPC, hasTPC, //! Flag to check if track has a TPC match + [](uint8_t detectorMap) -> bool { return detectorMap & o2::aod::track::TPC; }); +DECLARE_SOA_DYNAMIC_COLUMN(HasTOF, hasTOF, //! Flag to check if track has a TOF measurement + [](uint8_t detectorMap) -> bool { return detectorMap & o2::aod::track::TOF; }); +DECLARE_SOA_DYNAMIC_COLUMN(IsPVContributor, isPVContributor, //! Run 3: Has this track contributed to the collision vertex fit + [](uint8_t flags) -> bool { return (flags & o2::aod::track::PVContributor) == o2::aod::track::PVContributor; }); +DECLARE_SOA_DYNAMIC_COLUMN(PIDForTracking, pidForTracking, //! PID hypothesis used during tracking. See the constants in the class PID in PID.h + [](uint32_t flags) -> uint32_t { return flags >> 28; }); +DECLARE_SOA_DYNAMIC_COLUMN(TPCNClsFound, tpcNClsFound, //! Number of found TPC clusters + [](uint8_t tpcNClsFindable, int8_t tpcNClsFindableMinusFound) -> int16_t { return (int16_t)tpcNClsFindable - tpcNClsFindableMinusFound; }); +DECLARE_SOA_DYNAMIC_COLUMN(TPCNClsCrossedRows, tpcNClsCrossedRows, //! Number of crossed TPC Rows + [](uint8_t tpcNClsFindable, int8_t TPCNClsFindableMinusCrossedRows) -> int16_t { return (int16_t)tpcNClsFindable - TPCNClsFindableMinusCrossedRows; }); +DECLARE_SOA_DYNAMIC_COLUMN(TPCCrossedRowsOverFindableCls, tpcCrossedRowsOverFindableCls, //! Ratio crossed rows over findable clusters + [](uint8_t tpcNClsFindable, int8_t tpcNClsFindableMinusCrossedRows) -> float { + int16_t tpcNClsCrossedRows = (int16_t)tpcNClsFindable - tpcNClsFindableMinusCrossedRows; + return (float)tpcNClsCrossedRows / (float)tpcNClsFindable; + }); + +// track PID definition +DECLARE_SOA_COLUMN(TPCNSigmaPr, tpcNSigmaPr, float); //! Nsigma separation with the TPC detector for proton +DECLARE_SOA_COLUMN(TPCNSigmaDe, tpcNSigmaDe, float); //! Nsigma separation with the TPC detector for deuteron +DECLARE_SOA_COLUMN(TPCNSigmaPi, tpcNSigmaPi, float); //! Nsigma separation with the TPC detector for pion +DECLARE_SOA_COLUMN(TOFNSigmaDe, tofNSigmaDe, float); //! Nsigma separation with the TOF detector for deuteron (recalculated) + +} // namespace reducedtracks3body + +DECLARE_SOA_TABLE_FULL(StoredReducedTracksIU, "ReducedTracks_IU", "AOD", "REDUCEDTRACK_IU", //! On disk version of the track parameters at inner most update (e.g. ITS) as it comes from the tracking + o2::soa::Index<>, reducedtracks3body::CollisionId, + reducedtracks3body::X, reducedtracks3body::Alpha, + reducedtracks3body::Y, reducedtracks3body::Z, reducedtracks3body::Snp, reducedtracks3body::Tgl, + reducedtracks3body::Signed1Pt, + // cov matrix + reducedtracks3body::SigmaY, reducedtracks3body::SigmaZ, reducedtracks3body::SigmaSnp, reducedtracks3body::SigmaTgl, reducedtracks3body::Sigma1Pt, + reducedtracks3body::RhoZY, reducedtracks3body::RhoSnpY, reducedtracks3body::RhoSnpZ, reducedtracks3body::RhoTglY, reducedtracks3body::RhoTglZ, + reducedtracks3body::RhoTglSnp, reducedtracks3body::Rho1PtY, reducedtracks3body::Rho1PtZ, reducedtracks3body::Rho1PtSnp, reducedtracks3body::Rho1PtTgl, + // tracks extra + reducedtracks3body::TPCInnerParam, reducedtracks3body::Flags, reducedtracks3body::ITSClusterSizes, + reducedtracks3body::TPCNClsFindable, reducedtracks3body::TPCNClsFindableMinusFound, reducedtracks3body::TPCNClsFindableMinusCrossedRows, + reducedtracks3body::TRDPattern, reducedtracks3body::TPCChi2NCl, reducedtracks3body::TOFChi2, + reducedtracks3body::TPCSignal, reducedtracks3body::TOFExpMom, + // TPC PID + reducedtracks3body::TPCNSigmaPr, reducedtracks3body::TPCNSigmaPi, reducedtracks3body::TPCNSigmaDe, + reducedtracks3body::TOFNSigmaDe, + + // ----------- dynmaic columns ------------ + // tracks IU + reducedtracks3body::Px, + reducedtracks3body::Py, + reducedtracks3body::Pz, + reducedtracks3body::Rapidity, + reducedtracks3body::Sign, + // tracks extra + reducedtracks3body::PIDForTracking, + reducedtracks3body::IsPVContributor, + reducedtracks3body::HasTPC, + reducedtracks3body::HasTOF, + reducedtracks3body::TPCNClsFound, + reducedtracks3body::TPCNClsCrossedRows, + reducedtracks3body::ITSClsSizeInLayer, + reducedtracks3body::TPCCrossedRowsOverFindableCls); + +DECLARE_SOA_EXTENDED_TABLE(ReducedTracksIU, StoredReducedTracksIU, "EXREDUCEDTRACK_IU", 0, //! Track parameters at inner most update (e.g. ITS) as it comes from the tracking + reducedtracks3body::Pt, + reducedtracks3body::P, + reducedtracks3body::Eta, + reducedtracks3body::Phi, + // cov matrix + reducedtracks3body::CYY, + reducedtracks3body::CZY, + reducedtracks3body::CZZ, + reducedtracks3body::CSnpY, + reducedtracks3body::CSnpZ, + reducedtracks3body::CSnpSnp, + reducedtracks3body::CTglY, + reducedtracks3body::CTglZ, + reducedtracks3body::CTglSnp, + reducedtracks3body::CTglTgl, + reducedtracks3body::C1PtY, + reducedtracks3body::C1PtZ, + reducedtracks3body::C1PtSnp, + reducedtracks3body::C1PtTgl, + reducedtracks3body::C1Pt21Pt2, + // tracks extra + reducedtracks3body::DetectorMap); + +using ReducedTrackIU = ReducedTracksIU::iterator; + +namespace reduceddecay3body +{ +DECLARE_SOA_INDEX_COLUMN_FULL(Track0, track0, int, ReducedTracksIU, "_0"); //! Track 0 index +DECLARE_SOA_INDEX_COLUMN_FULL(Track1, track1, int, ReducedTracksIU, "_1"); //! Track 1 index +DECLARE_SOA_INDEX_COLUMN_FULL(Track2, track2, int, ReducedTracksIU, "_2"); //! Track 2 index +DECLARE_SOA_INDEX_COLUMN_FULL(Collision, collision, int, ReducedCollisions, ""); //! Collision index +} // namespace reduceddecay3body + +DECLARE_SOA_TABLE(ReducedDecay3Bodys, "AOD", "REDUCEDDECAY3BODY", //! reduced 3-body decay table + o2::soa::Index<>, reduceddecay3body::CollisionId, reduceddecay3body::Track0Id, reduceddecay3body::Track1Id, reduceddecay3body::Track2Id); + +using ReducedDecay3BodysLinked = soa::Join; +using ReducedDecay3BodyLinked = ReducedDecay3BodysLinked::iterator; + +} // namespace o2::aod + +#endif // PWGLF_DATAMODEL_REDUCED3BODYTABLES_H_ diff --git a/PWGLF/DataModel/pidTOFGeneric.h b/PWGLF/DataModel/pidTOFGeneric.h index b43d4ecc9ff..a0d287a7643 100644 --- a/PWGLF/DataModel/pidTOFGeneric.h +++ b/PWGLF/DataModel/pidTOFGeneric.h @@ -88,19 +88,76 @@ class TofPidNewCollision template float GetTOFNSigma(TTrack const& track, TCollision const& originalcol, TCollision const& correctedcol, bool EnableBCAO2D = true); + + float GetTOFNSigma(TTrack const& track); + float GetTOFNSigma(o2::track::PID::ID pidId, TTrack const& track); + + float CalculateTOFNSigma(o2::track::PID::ID pidId, TTrack const& track, double tofsignal, double evTime, double evTimeErr) + { + + float expSigma, tofNsigma = -999; + + switch (pidId) { + case 0: + expSigma = responseEl.GetExpectedSigma(mRespParamsV2, track, tofsignal, evTimeErr); + tofNsigma = (tofsignal - evTime - responseEl.GetCorrectedExpectedSignal(mRespParamsV2, track)) / expSigma; + break; + case 1: + expSigma = responseMu.GetExpectedSigma(mRespParamsV2, track, tofsignal, evTimeErr); + tofNsigma = (tofsignal - evTime - responseMu.GetCorrectedExpectedSignal(mRespParamsV2, track)) / expSigma; + break; + case 2: + expSigma = responsePi.GetExpectedSigma(mRespParamsV2, track, tofsignal, evTimeErr); + tofNsigma = (tofsignal - evTime - responsePi.GetCorrectedExpectedSignal(mRespParamsV2, track)) / expSigma; + break; + case 3: + expSigma = responseKa.GetExpectedSigma(mRespParamsV2, track, tofsignal, evTimeErr); + tofNsigma = (tofsignal - evTime - responseKa.GetCorrectedExpectedSignal(mRespParamsV2, track)) / expSigma; + break; + case 4: + expSigma = responsePr.GetExpectedSigma(mRespParamsV2, track, tofsignal, evTimeErr); + tofNsigma = (tofsignal - evTime - responsePr.GetCorrectedExpectedSignal(mRespParamsV2, track)) / expSigma; + break; + case 5: + expSigma = responseDe.GetExpectedSigma(mRespParamsV2, track, tofsignal, evTimeErr); + tofNsigma = (tofsignal - evTime - responseDe.GetCorrectedExpectedSignal(mRespParamsV2, track)) / expSigma; + break; + case 6: + expSigma = responseTr.GetExpectedSigma(mRespParamsV2, track, tofsignal, evTimeErr); + tofNsigma = (tofsignal - evTime - responseTr.GetCorrectedExpectedSignal(mRespParamsV2, track)) / expSigma; + break; + case 7: + expSigma = responseHe.GetExpectedSigma(mRespParamsV2, track, tofsignal, evTimeErr); + tofNsigma = (tofsignal - evTime - responseHe.GetCorrectedExpectedSignal(mRespParamsV2, track)) / expSigma; + break; + case 8: + expSigma = responseAl.GetExpectedSigma(mRespParamsV2, track, tofsignal, evTimeErr); + tofNsigma = (tofsignal - evTime - responseAl.GetCorrectedExpectedSignal(mRespParamsV2, track)) / expSigma; + break; + default: + LOG(fatal) << "Wrong particle ID in TofPidSecondary class"; + return -999; + } + + return tofNsigma; + } }; template template float TofPidNewCollision::GetTOFNSigma(o2::track::PID::ID pidId, TTrack const& track, TCollision const& originalcol, TCollision const& correctedcol, bool EnableBCAO2D) { + + if (!track.has_collision() || !track.hasTOF()) { + return -999; + } + float mMassHyp = o2::track::pid_constants::sMasses2Z[track.pidForTracking()]; float expTime = track.length() * sqrt((mMassHyp * mMassHyp) + (track.tofExpMom() * track.tofExpMom())) / (kCSPEED * track.tofExpMom()); // L*E/(p*c) = L/v float evTime = correctedcol.evTime(); float evTimeErr = correctedcol.evTimeErr(); float tofsignal = track.trackTime() * 1000 + expTime; // in ps - float expSigma, tofNsigma; if (originalcol.globalIndex() == correctedcol.globalIndex()) { evTime = track.evTimeForTrack(); @@ -121,48 +178,7 @@ float TofPidNewCollision::GetTOFNSigma(o2::track::PID::ID pidId, TTrack } } - switch (pidId) { - case 0: - expSigma = responseEl.GetExpectedSigma(mRespParamsV2, track, tofsignal, evTimeErr); - tofNsigma = (tofsignal - evTime - responseEl.GetCorrectedExpectedSignal(mRespParamsV2, track)) / expSigma; - break; - case 1: - expSigma = responseMu.GetExpectedSigma(mRespParamsV2, track, tofsignal, evTimeErr); - tofNsigma = (tofsignal - evTime - responseMu.GetCorrectedExpectedSignal(mRespParamsV2, track)) / expSigma; - break; - case 2: - expSigma = responsePi.GetExpectedSigma(mRespParamsV2, track, tofsignal, evTimeErr); - tofNsigma = (tofsignal - evTime - responsePi.GetCorrectedExpectedSignal(mRespParamsV2, track)) / expSigma; - break; - case 3: - expSigma = responseKa.GetExpectedSigma(mRespParamsV2, track, tofsignal, evTimeErr); - tofNsigma = (tofsignal - evTime - responseKa.GetCorrectedExpectedSignal(mRespParamsV2, track)) / expSigma; - break; - case 4: - expSigma = responsePr.GetExpectedSigma(mRespParamsV2, track, tofsignal, evTimeErr); - tofNsigma = (tofsignal - evTime - responsePr.GetCorrectedExpectedSignal(mRespParamsV2, track)) / expSigma; - break; - case 5: - expSigma = responseDe.GetExpectedSigma(mRespParamsV2, track, tofsignal, evTimeErr); - tofNsigma = (tofsignal - evTime - responseDe.GetCorrectedExpectedSignal(mRespParamsV2, track)) / expSigma; - break; - case 6: - expSigma = responseTr.GetExpectedSigma(mRespParamsV2, track, tofsignal, evTimeErr); - tofNsigma = (tofsignal - evTime - responseTr.GetCorrectedExpectedSignal(mRespParamsV2, track)) / expSigma; - break; - case 7: - expSigma = responseHe.GetExpectedSigma(mRespParamsV2, track, tofsignal, evTimeErr); - tofNsigma = (tofsignal - evTime - responseHe.GetCorrectedExpectedSignal(mRespParamsV2, track)) / expSigma; - break; - case 8: - expSigma = responseAl.GetExpectedSigma(mRespParamsV2, track, tofsignal, evTimeErr); - tofNsigma = (tofsignal - evTime - responseAl.GetCorrectedExpectedSignal(mRespParamsV2, track)) / expSigma; - break; - default: - LOG(fatal) << "Wrong particle ID in TofPidSecondary class"; - return -999; - } - + float tofNsigma = CalculateTOFNSigma(pidId, track, tofsignal, evTime, evTimeErr); return tofNsigma; } @@ -173,6 +189,31 @@ float TofPidNewCollision::GetTOFNSigma(TTrack const& track, TCollision c return GetTOFNSigma(pidType, track, originalcol, correctedcol, EnableBCAO2D); } +template +float TofPidNewCollision::GetTOFNSigma(o2::track::PID::ID pidId, TTrack const& track) +{ + + if (!track.has_collision() || !track.hasTOF()) { + return -999; + } + + float mMassHyp = o2::track::pid_constants::sMasses2Z[track.pidForTracking()]; + float expTime = track.length() * sqrt((mMassHyp * mMassHyp) + (track.tofExpMom() * track.tofExpMom())) / (kCSPEED * track.tofExpMom()); // L*E/(p*c) = L/v + + float evTime = track.evTimeForTrack(); + float evTimeErr = track.evTimeErrForTrack(); + float tofsignal = track.trackTime() * 1000 + expTime; // in ps + + float tofNsigma = CalculateTOFNSigma(pidId, track, tofsignal, evTime, evTimeErr); + return tofNsigma; +} + +template +float TofPidNewCollision::GetTOFNSigma(TTrack const& track) +{ + return GetTOFNSigma(pidType, track); +} + } // namespace pidtofgeneric } // namespace o2::aod diff --git a/PWGLF/TableProducer/Nuspex/CMakeLists.txt b/PWGLF/TableProducer/Nuspex/CMakeLists.txt index be724d080ad..076349f3af5 100644 --- a/PWGLF/TableProducer/Nuspex/CMakeLists.txt +++ b/PWGLF/TableProducer/Nuspex/CMakeLists.txt @@ -98,3 +98,8 @@ o2physics_add_dpl_workflow(tr-he-analysis SOURCES trHeAnalysis.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(reduced3body-creator + SOURCES reduced3bodyCreator.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::EventFilteringUtils O2::TOFBase + COMPONENT_NAME Analysis) diff --git a/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx b/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx index 80b27744557..693668581f2 100644 --- a/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx +++ b/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx @@ -32,6 +32,7 @@ #include "Common/Core/RecoDecay.h" #include "Common/Core/trackUtilities.h" #include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/DataModel/Reduced3BodyTables.h" #include "PWGLF/DataModel/Vtx3BodyTables.h" #include "PWGLF/DataModel/pidTOFGeneric.h" #include "Common/Core/TrackSelection.h" @@ -69,7 +70,7 @@ using namespace o2::framework::expressions; using std::array; using FullTracksExtIU = soa::Join; -using FullTracksExtPIDIU = soa::Join; +using FullTracksExtPIDIU = soa::Join; using ColwithEvTimes = o2::soa::Join; using FullCols = o2::soa::Join; @@ -324,6 +325,7 @@ struct decay3bodyBuilder { // Filters and slices // Filter collisionFilter = (aod::evsel::sel8 == true && nabs(aod::collision::posZ) < 10.f); Preslice perCollision = o2::aod::decay3body::collisionId; + Preslice perReducedCollision = o2::aod::reduceddecay3body::collisionId; int mRunNumber; float d_bz; @@ -399,7 +401,7 @@ struct decay3bodyBuilder { fitter3body.setMatCorrType(matCorr); // Add histograms separately for different process functions - if (doprocessRun3 == true || doprocessRun3EM == true || doprocessRun3EMLikeSign == true) { + if (doprocessRun3 == true || doprocessRun3Reduced || doprocessRun3EM == true || doprocessRun3EMLikeSign == true) { registry.add("hEventCounter", "hEventCounter", HistType::kTH1F, {{1, 0.0f, 1.0f}}); auto hVtx3BodyCounter = registry.add("hVtx3BodyCounter", "hVtx3BodyCounter", HistType::kTH1F, {{6, 0.0f, 6.0f}}); hVtx3BodyCounter->GetXaxis()->SetBinLabel(1, "Total"); @@ -511,7 +513,7 @@ struct decay3bodyBuilder { KFParticle::SetField(d_bz); #endif o2::parameters::GRPMagField grpmag; - if (fabs(d_bz) > 1e-5) { + if (std::fabs(d_bz) > 1e-5) { grpmag.setL3Current(30000.f / (d_bz / 5.0f)); } o2::base::Propagator::initFieldFromGRP(&grpmag); @@ -644,7 +646,7 @@ struct decay3bodyBuilder { //------------------------------------------------------------------ // 3body candidate builder template - void fillVtxCand(TCollisionTable const& collision, TTrackTable const& t0, TTrackTable const& t1, TTrackTable const& t2, int64_t decay3bodyId, int bachelorcharge = 1) + void fillVtxCand(TCollisionTable const& collision, TTrackTable const& t0, TTrackTable const& t1, TTrackTable const& t2, int64_t decay3bodyId, int bachelorcharge = 1, double tofNSigmaBach = -999) { registry.fill(HIST("hVtx3BodyCounter"), kVtxAll); @@ -654,18 +656,6 @@ struct decay3bodyBuilder { } registry.fill(HIST("hVtx3BodyCounter"), kVtxTPCNcls); - // Recalculate the TOF PID - double tofNSigmaBach = -999; - if (t2.has_collision() && t2.hasTOF()) { - if (decay3bodyId == -1) { - // for event-mixing, the collisionId of tracks not equal to global index - tofNSigmaBach = bachelorTOFPID.GetTOFNSigma(t2, collision, collision); - } else { - auto originalcol = t2.template collision_as(); - tofNSigmaBach = bachelorTOFPID.GetTOFNSigma(t2, originalcol, collision); - } - } - if (enablePidCut) { if (t2.sign() > 0) { if (!checkPIDH3L(t0, t1, t2, tofNSigmaBach)) @@ -728,7 +718,7 @@ struct decay3bodyBuilder { } registry.fill(HIST("hVtx3BodyCounter"), kVtxDcaDau); - float VtxcosPA = RecoDecay::cpa(array{collision.posX(), collision.posY(), collision.posZ()}, array{pos[0], pos[1], pos[2]}, array{p3B[0], p3B[1], p3B[2]}); + float VtxcosPA = RecoDecay::cpa(array{collision.posX(), collision.posY(), collision.posZ()}, std::array{pos[0], pos[1], pos[2]}, std::array{p3B[0], p3B[1], p3B[2]}); if (VtxcosPA < minCosPA3body) { return; } @@ -1452,21 +1442,52 @@ struct decay3bodyBuilder { initCCDB(bc); registry.fill(HIST("hEventCounter"), 0.5); - const auto& d3bodysInCollision = decay3bodys.sliceBy(perCollision, collision.globalIndex()); - for (auto& d3body : d3bodysInCollision) { + const auto& d3bodys_thisCollision = decay3bodys.sliceBy(perCollision, collision.globalIndex()); + for (const auto& d3body : d3bodys_thisCollision) { auto t0 = d3body.template track0_as(); auto t1 = d3body.template track1_as(); auto t2 = d3body.template track2_as(); - fillVtxCand(collision, t0, t1, t2, d3body.globalIndex(), bachelorcharge); + + // Recalculate the TOF PID + double tofNSigmaBach = -999; + if (t2.has_collision() && t2.hasTOF()) { + auto originalcol = t2.template collision_as(); + tofNSigmaBach = bachelorTOFPID.GetTOFNSigma(t2, originalcol, collision); + } + + fillVtxCand(collision, t0, t1, t2, d3body.globalIndex(), bachelorcharge, tofNSigmaBach); } } - for (auto& candVtx : vtxCandidates) { + for (const auto& candVtx : vtxCandidates) { fillVtx3BodyTable(candVtx); } } PROCESS_SWITCH(decay3bodyBuilder, processRun3, "Produce DCA fitter decay3body tables", true); + //------------------------------------------------------------------ + void processRun3Reduced(aod::ReducedCollisions const& collisions, aod::ReducedTracksIU const&, aod::ReducedDecay3Bodys const& decay3bodys, aod::BCsWithTimestamps const&) + { + vtxCandidates.clear(); + + registry.fill(HIST("hEventCounter"), 0.5, collisions.size()); + + for (const auto& d3body : decay3bodys) { + auto t0 = d3body.template track0_as(); + auto t1 = d3body.template track1_as(); + auto t2 = d3body.template track2_as(); + auto collision = d3body.template collision_as(); + auto bc = collision.bc_as(); + initCCDB(bc); + fillVtxCand(collision, t0, t1, t2, d3body.globalIndex(), bachelorcharge, t2.tofNSigmaDe()); + } + + for (const auto& candVtx : vtxCandidates) { + fillVtx3BodyTable(candVtx); + } + } + PROCESS_SWITCH(decay3bodyBuilder, processRun3Reduced, "Produce DCA fitter decay3body tables with reduced data", false); + //------------------------------------------------------------------ // Event-mixing background void processRun3EM(FullCols const& collisions, TrackExtPIDIUwithEvTimes const& tracksIU, aod::BCsWithTimestamps const&) @@ -1491,7 +1512,7 @@ struct decay3bodyBuilder { candBachelors.bindTable(tracksIU); candAntiBachelors.bindTable(tracksIU); - for (auto& [c1, tracks1, c2, tracks2] : pair) { + for (const auto& [c1, tracks1, c2, tracks2] : pair) { if (EMTrackSel.em_event_sel8_selection && (!c1.sel8() || !c2.sel8())) { continue; } @@ -1505,10 +1526,12 @@ struct decay3bodyBuilder { auto antibachelors = candAntiBachelors->sliceByCached(aod::track::collisionId, c2.globalIndex(), cache); for (auto const& [tpos, tneg, tbach] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(protons, pionsminus, bachelors))) { - fillVtxCand(c1, tpos, tneg, tbach, -1, bachelorcharge); + double tofNSigmaBach = bachelorTOFPID.GetTOFNSigma(tbach); // Recalculate the TOF PID + fillVtxCand(c1, tpos, tneg, tbach, -1, bachelorcharge, tofNSigmaBach); } for (auto const& [tpos, tneg, tbach] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(pionsplus, antiprotons, antibachelors))) { - fillVtxCand(c1, tpos, tneg, tbach, -1, bachelorcharge); + double tofNSigmaBach = bachelorTOFPID.GetTOFNSigma(tbach); // Recalculate the TOF PID + fillVtxCand(c1, tpos, tneg, tbach, -1, bachelorcharge, tofNSigmaBach); } } @@ -1517,7 +1540,7 @@ struct decay3bodyBuilder { return a.collisionId < b.collisionId; }); - for (auto& candVtx : vtxCandidates) { + for (const auto& candVtx : vtxCandidates) { fillVtx3BodyTable(candVtx); } } @@ -1547,7 +1570,7 @@ struct decay3bodyBuilder { candBachelors.bindTable(tracksIU); candAntiBachelors.bindTable(tracksIU); - for (auto& [c1, tracks1, c2, tracks2] : pair) { + for (const auto& [c1, tracks1, c2, tracks2] : pair) { if (EMTrackSel.em_event_sel8_selection && (!c1.sel8() || !c2.sel8())) { continue; } @@ -1561,10 +1584,12 @@ struct decay3bodyBuilder { auto antibachelors = candAntiBachelors->sliceByCached(aod::track::collisionId, c2.globalIndex(), cache); for (auto const& [tpos, tneg, tbach] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(protons, pionsminus, antibachelors))) { - fillVtxCand(c1, tpos, tneg, tbach, -1, bachelorcharge); + double tofNSigmaBach = bachelorTOFPID.GetTOFNSigma(tbach); // Recalculate the TOF PID + fillVtxCand(c1, tpos, tneg, tbach, -1, bachelorcharge, tofNSigmaBach); } for (auto const& [tpos, tneg, tbach] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(pionsplus, antiprotons, bachelors))) { - fillVtxCand(c1, tpos, tneg, tbach, -1, bachelorcharge); + double tofNSigmaBach = bachelorTOFPID.GetTOFNSigma(tbach); // Recalculate the TOF PID + fillVtxCand(c1, tpos, tneg, tbach, -1, bachelorcharge, tofNSigmaBach); } } @@ -1573,7 +1598,7 @@ struct decay3bodyBuilder { return a.collisionId < b.collisionId; }); - for (auto& candVtx : vtxCandidates) { + for (const auto& candVtx : vtxCandidates) { fillVtx3BodyTable(candVtx); } } @@ -1727,13 +1752,14 @@ struct decay3bodyDataLinkBuilder { void init(InitContext const&) {} - void process(aod::Decay3Bodys const& decay3bodytable, aod::Vtx3BodyDatas const& vtxdatatable) + template + void buildDecay3BodyLabel(TDecay3Bodys const& decay3bodytable, TVtx3BodyDatas const& vtxdatatable) { std::vector lIndices; lIndices.reserve(decay3bodytable.size()); for (int ii = 0; ii < decay3bodytable.size(); ii++) lIndices[ii] = -1; - for (auto& vtxdata : vtxdatatable) { + for (const auto& vtxdata : vtxdatatable) { if (vtxdata.decay3bodyId() != -1) { lIndices[vtxdata.decay3bodyId()] = vtxdata.globalIndex(); } @@ -1742,6 +1768,18 @@ struct decay3bodyDataLinkBuilder { vtxdataLink(lIndices[ii]); } } + + void processStandard(aod::Decay3Bodys const& decay3bodytable, aod::Vtx3BodyDatas const& vtxdatatable) + { + buildDecay3BodyLabel(decay3bodytable, vtxdatatable); + } + PROCESS_SWITCH(decay3bodyDataLinkBuilder, processStandard, "Produce label from decay3body to vtx3body", true); + + void processReduced(aod::ReducedDecay3Bodys const& decay3bodytable, aod::Vtx3BodyDatas const& vtxdatatable) + { + buildDecay3BodyLabel(decay3bodytable, vtxdatatable); + } + PROCESS_SWITCH(decay3bodyDataLinkBuilder, processReduced, "Produce label from reducedDecay3body to vtx3body", false); }; struct kfdecay3bodyDataLinkBuilder { @@ -1814,7 +1852,7 @@ struct decay3bodyLabelBuilder { lIndices[ii] = -1; } - for (auto& decay3body : decay3bodys) { + for (const auto& decay3body : decay3bodys) { int lLabel = -1; int lPDG = -1; @@ -1842,9 +1880,9 @@ struct decay3bodyLabelBuilder { continue; } - for (auto& lMother0 : lMCTrack0.mothers_as()) { - for (auto& lMother1 : lMCTrack1.mothers_as()) { - for (auto& lMother2 : lMCTrack2.mothers_as()) { + for (const auto& lMother0 : lMCTrack0.mothers_as()) { + for (const auto& lMother1 : lMCTrack1.mothers_as()) { + for (const auto& lMother2 : lMCTrack2.mothers_as()) { if (lMother0.globalIndex() == lMother1.globalIndex() && lMother0.globalIndex() == lMother2.globalIndex()) { lGlobalIndex = lMother1.globalIndex(); lPt = lMother1.pt(); @@ -1864,7 +1902,7 @@ struct decay3bodyLabelBuilder { // Intended for hypertriton cross-checks only if (lPDG == 1010010030 && lMCTrack0.pdgCode() == 2212 && lMCTrack1.pdgCode() == -211 && lMCTrack2.pdgCode() == 1000010020) { lLabel = lGlobalIndex; - double hypertritonMCMass = RecoDecay::m(array{array{lMCTrack0.px(), lMCTrack0.py(), lMCTrack0.pz()}, array{lMCTrack1.px(), lMCTrack1.py(), lMCTrack1.pz()}, array{lMCTrack2.px(), lMCTrack2.py(), lMCTrack2.pz()}}, array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged, o2::constants::physics::MassDeuteron}); + double hypertritonMCMass = RecoDecay::m(std::array{std::array{lMCTrack0.px(), lMCTrack0.py(), lMCTrack0.pz()}, std::array{lMCTrack1.px(), lMCTrack1.py(), lMCTrack1.pz()}, std::array{lMCTrack2.px(), lMCTrack2.py(), lMCTrack2.pz()}}, std::array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged, o2::constants::physics::MassDeuteron}); registry.fill(HIST("hLabelCounter"), 2.5); registry.fill(HIST("hHypertritonMCPt"), lPt); registry.fill(HIST("hHypertritonMCLifetime"), MClifetime); @@ -1872,7 +1910,7 @@ struct decay3bodyLabelBuilder { } if (lPDG == -1010010030 && lMCTrack0.pdgCode() == 211 && lMCTrack1.pdgCode() == -2212 && lMCTrack2.pdgCode() == -1000010020) { lLabel = lGlobalIndex; - double antiHypertritonMCMass = RecoDecay::m(array{array{lMCTrack0.px(), lMCTrack0.py(), lMCTrack0.pz()}, array{lMCTrack1.px(), lMCTrack1.py(), lMCTrack1.pz()}, array{lMCTrack2.px(), lMCTrack2.py(), lMCTrack2.pz()}}, array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassProton, o2::constants::physics::MassDeuteron}); + double antiHypertritonMCMass = RecoDecay::m(std::array{std::array{lMCTrack0.px(), lMCTrack0.py(), lMCTrack0.pz()}, std::array{lMCTrack1.px(), lMCTrack1.py(), lMCTrack1.pz()}, std::array{lMCTrack2.px(), lMCTrack2.py(), lMCTrack2.pz()}}, std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassProton, o2::constants::physics::MassDeuteron}); registry.fill(HIST("hLabelCounter"), 2.5); registry.fill(HIST("hAntiHypertritonMCPt"), lPt); registry.fill(HIST("hAntiHypertritonMCLifetime"), MClifetime); diff --git a/PWGLF/TableProducer/Nuspex/reduced3bodyCreator.cxx b/PWGLF/TableProducer/Nuspex/reduced3bodyCreator.cxx new file mode 100644 index 00000000000..e6d38c83bc9 --- /dev/null +++ b/PWGLF/TableProducer/Nuspex/reduced3bodyCreator.cxx @@ -0,0 +1,318 @@ +// 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 Task to produce reduced AO2Ds for use in the hypertriton 3body reconstruction with the decay3bodybuilder.cxx +/// \author Yuanzhe Wang +/// \author Carolina Reetz + +#include +#include +#include +#include +#include +#include + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoAHelpers.h" +#include "ReconstructionDataFormats/Track.h" +#include "Common/Core/RecoDecay.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/Centrality.h" +#include "PWGLF/DataModel/pidTOFGeneric.h" +#include "PWGLF/DataModel/Reduced3BodyTables.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/PIDResponse.h" +#include "Common/Core/PID/PIDTOF.h" +#include "TableHelper.h" + +#include "EventFiltering/Zorro.h" +#include "EventFiltering/ZorroSummary.h" + +#include "DetectorsBase/GeometryManager.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "CCDB/BasicCCDBManager.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +using FullTracksExtIU = soa::Join; +using FullTracksExtPIDIU = soa::Join; + +using ColwithEvTimes = o2::soa::Join; +using ColwithEvTimesMultsCents = o2::soa::Join; +using TrackExtIUwithEvTimes = soa::Join; +using TrackExtPIDIUwithEvTimes = soa::Join; + +struct reduced3bodyCreator { + + Produces reducedCollisions; + Produces reducedPVMults; + Produces reducedCentFTOCs; + Produces reducedDecay3Bodys; + Produces reducedFullTracksPIDIU; + + Service ccdb; + Zorro zorro; + OutputObj zorroSummary{"zorroSummary"}; + + o2::aod::pidtofgeneric::TofPidNewCollision bachelorTOFPID; + + std::vector daughterTracks; + + Configurable event_sel8_selection{"event_sel8_selection", true, "event selection count post sel8 cut"}; + Configurable mc_event_selection{"mc_event_selection", true, "mc event selection count post kIsTriggerTVX and kNoTimeFrameBorder"}; + Configurable event_posZ_selection{"event_posZ_selection", true, "event selection count post poZ cut"}; + // CCDB options + Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + // CCDB TOF PID paras + Configurable timestamp{"ccdb-timestamp", -1, "timestamp of the object"}; + Configurable paramFileName{"paramFileName", "", "Path to the parametrization object. If empty the parametrization is not taken from file"}; + Configurable parametrizationPath{"parametrizationPath", "TOF/Calib/Params", "Path of the TOF parametrization on the CCDB or in the file, if the paramFileName is not empty"}; + Configurable passName{"passName", "", "Name of the pass inside of the CCDB parameter collection. If empty, the automatically deceted from metadata (to be implemented!!!)"}; + Configurable timeShiftCCDBPath{"timeShiftCCDBPath", "", "Path of the TOF time shift vs eta. If empty none is taken"}; + Configurable loadResponseFromCCDB{"loadResponseFromCCDB", false, "Flag to load the response from the CCDB"}; + Configurable fatalOnPassNotAvailable{"fatalOnPassNotAvailable", true, "Flag to throw a fatal if the pass is not available in the retrieved CCDB object"}; + // Zorro counting + Configurable cfgSkimmedProcessing{"cfgSkimmedProcessing", false, "Skimmed dataset processing"}; + + int mRunNumber; + o2::pid::tof::TOFResoParamsV2 mRespParamsV2; + + HistogramRegistry registry{"registry", {}}; + + void init(InitContext&) + { + mRunNumber = 0; + zorroSummary.setObject(zorro.getZorroSummary()); + bachelorTOFPID.SetPidType(o2::track::PID::Deuteron); + + ccdb->setURL(ccdburl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); + + registry.add("hAllSelEventsVtxZ", "hAllSelEventsVtxZ", HistType::kTH1F, {{500, -15.0f, 15.0f, "PV Z (cm)"}}); + + auto hEventCounter = registry.add("hEventCounter", "hEventCounter", HistType::kTH1F, {{4, 0.0f, 4.0f}}); + hEventCounter->GetXaxis()->SetBinLabel(1, "total"); + hEventCounter->GetXaxis()->SetBinLabel(2, "sel8"); + hEventCounter->GetXaxis()->SetBinLabel(3, "vertexZ"); + hEventCounter->GetXaxis()->SetBinLabel(4, "reduced"); + hEventCounter->LabelsOption("v"); + + auto hEventCounterZorro = registry.add("hEventCounterZorro", "hEventCounterZorro", HistType::kTH1D, {{2, -0.5, 1.5}}); + hEventCounterZorro->GetXaxis()->SetBinLabel(1, "Zorro before evsel"); + hEventCounterZorro->GetXaxis()->SetBinLabel(2, "Zorro after evsel"); + } + + void initCCDB(aod::BCsWithTimestamps::iterator const& bc) + { + // In case override, don't proceed, please - no CCDB access required + if (mRunNumber == bc.runNumber()) { + return; + } + + mRunNumber = bc.runNumber(); + if (cfgSkimmedProcessing) { + zorro.initCCDB(ccdb.service, bc.runNumber(), bc.timestamp(), "fH3L3Body"); + zorro.populateHistRegistry(registry, bc.runNumber()); + } + + // Initial TOF PID Paras, copied from PIDTOF.h + timestamp.value = bc.timestamp(); + ccdb->setTimestamp(timestamp.value); + // Not later than now objects + ccdb->setCreatedNotAfter(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count()); + // TODO: implement the automatic pass name detection from metadata + if (passName.value == "") { + passName.value = "unanchored"; // temporary default + LOG(warning) << "Passed autodetect mode for pass, not implemented yet, waiting for metadata. Taking '" << passName.value << "'"; + } + LOG(info) << "Using parameter collection, starting from pass '" << passName.value << "'"; + + const std::string fname = paramFileName.value; + if (!fname.empty()) { // Loading the parametrization from file + LOG(info) << "Loading exp. sigma parametrization from file " << fname << ", using param: " << parametrizationPath.value; + if (1) { + o2::tof::ParameterCollection paramCollection; + paramCollection.loadParamFromFile(fname, parametrizationPath.value); + LOG(info) << "+++ Loaded parameter collection from file +++"; + if (!paramCollection.retrieveParameters(mRespParamsV2, passName.value)) { + if (fatalOnPassNotAvailable) { + LOGF(fatal, "Pass '%s' not available in the retrieved CCDB object", passName.value.data()); + } else { + LOGF(warning, "Pass '%s' not available in the retrieved CCDB object", passName.value.data()); + } + } else { + mRespParamsV2.setShiftParameters(paramCollection.getPars(passName.value)); + mRespParamsV2.printShiftParameters(); + } + } else { + mRespParamsV2.loadParamFromFile(fname.data(), parametrizationPath.value); + } + } else if (loadResponseFromCCDB) { // Loading it from CCDB + LOG(info) << "Loading exp. sigma parametrization from CCDB, using path: " << parametrizationPath.value << " for timestamp " << timestamp.value; + o2::tof::ParameterCollection* paramCollection = ccdb->getForTimeStamp(parametrizationPath.value, timestamp.value); + paramCollection->print(); + if (!paramCollection->retrieveParameters(mRespParamsV2, passName.value)) { // Attempt at loading the parameters with the pass defined + if (fatalOnPassNotAvailable) { + LOGF(fatal, "Pass '%s' not available in the retrieved CCDB object", passName.value.data()); + } else { + LOGF(warning, "Pass '%s' not available in the retrieved CCDB object", passName.value.data()); + } + } else { // Pass is available, load non standard parameters + mRespParamsV2.setShiftParameters(paramCollection->getPars(passName.value)); + mRespParamsV2.printShiftParameters(); + } + } + mRespParamsV2.print(); + if (timeShiftCCDBPath.value != "") { + if (timeShiftCCDBPath.value.find(".root") != std::string::npos) { + mRespParamsV2.setTimeShiftParameters(timeShiftCCDBPath.value, "gmean_Pos", true); + mRespParamsV2.setTimeShiftParameters(timeShiftCCDBPath.value, "gmean_Neg", false); + } else { + mRespParamsV2.setTimeShiftParameters(ccdb->getForTimeStamp(Form("%s/pos", timeShiftCCDBPath.value.c_str()), timestamp.value), true); + mRespParamsV2.setTimeShiftParameters(ccdb->getForTimeStamp(Form("%s/neg", timeShiftCCDBPath.value.c_str()), timestamp.value), false); + } + } + + bachelorTOFPID.SetParams(mRespParamsV2); + } + + void process(ColwithEvTimesMultsCents const& collisions, TrackExtPIDIUwithEvTimes const&, aod::Decay3Bodys const& decay3bodys, aod::BCsWithTimestamps const&) + { + + int lastCollisionID = -1; // collisionId of last analysed decay3body. Table is sorted. + + // Event counting + for (const auto& collision : collisions) { + // Zorro event counting + bool isZorroSelected = false; + if (cfgSkimmedProcessing) { + isZorroSelected = zorro.isSelected(collision.bc_as().globalBC()); + if (isZorroSelected) { + registry.fill(HIST("hEventCounterZorro"), 0.5); + } + } + + // Event selection + registry.fill(HIST("hEventCounter"), 0.5); + if (event_sel8_selection && !collision.sel8()) { + continue; + } + registry.fill(HIST("hEventCounter"), 1.5); + if (event_posZ_selection && (collision.posZ() >= 10.0f || collision.posZ() <= -10.0f)) { // 10cm + continue; + } + registry.fill(HIST("hEventCounter"), 2.5); + registry.fill(HIST("hAllSelEventsVtxZ"), collision.posZ()); + + if (cfgSkimmedProcessing && isZorroSelected) { + registry.fill(HIST("hEventCounterZorro"), 1.5); + } + } + + // Creat reduced table + for (const auto& d3body : decay3bodys) { + + daughterTracks.clear(); + + auto collision = d3body.template collision_as(); + + if (event_sel8_selection && !collision.sel8()) { + continue; + } + if (event_posZ_selection && (collision.posZ() >= 10.0f || collision.posZ() <= -10.0f)) { // 10cm + continue; + } + + auto bc = collision.bc_as(); + initCCDB(bc); + + // Save the collision + if (collision.globalIndex() != lastCollisionID) { + int runNumber = bc.runNumber(); + reducedCollisions( + collision.bcId(), + collision.posX(), collision.posY(), collision.posZ(), + collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ(), + collision.flags(), collision.chi2(), collision.numContrib(), + collision.collisionTime(), collision.collisionTimeRes(), + runNumber); + reducedPVMults(collision.multNTracksPV()); + reducedCentFTOCs(collision.centFT0C()); + + lastCollisionID = collision.globalIndex(); + } + + // Save daughter tracks + const auto daughter0 = d3body.template track0_as(); + const auto daughter1 = d3body.template track1_as(); + const auto daughter2 = d3body.template track2_as(); + + // TOF PID of bachelor must be calcualted here + // ---------------------------------------------- + auto originalcol = daughter2.template collision_as(); + double tofNSigmaBach = bachelorTOFPID.GetTOFNSigma(daughter2, originalcol, collision); + // ---------------------------------------------- + + // save reduced track table with decay3body daughters + daughterTracks.push_back(daughter0); + daughterTracks.push_back(daughter1); + daughterTracks.push_back(daughter2); + for (int i = 0; i < 3; i++) { + double tofNSigmaTrack = (i == 2) ? tofNSigmaBach : -999.; + reducedFullTracksPIDIU( + // TrackIU + // reducedTrackID + i, + reducedCollisions.lastIndex(), + daughterTracks[i].x(), daughterTracks[i].alpha(), + daughterTracks[i].y(), daughterTracks[i].z(), daughterTracks[i].snp(), daughterTracks[i].tgl(), + daughterTracks[i].signed1Pt(), + // TracksCovIU + daughterTracks[i].sigmaY(), daughterTracks[i].sigmaZ(), daughterTracks[i].sigmaSnp(), daughterTracks[i].sigmaTgl(), daughterTracks[i].sigma1Pt(), + daughterTracks[i].rhoZY(), daughterTracks[i].rhoSnpY(), daughterTracks[i].rhoSnpZ(), daughterTracks[i].rhoTglY(), daughterTracks[i].rhoTglZ(), + daughterTracks[i].rhoTglSnp(), daughterTracks[i].rho1PtY(), daughterTracks[i].rho1PtZ(), daughterTracks[i].rho1PtSnp(), daughterTracks[i].rho1PtTgl(), + // TracksExtra + daughterTracks[i].tpcInnerParam(), daughterTracks[i].flags(), daughterTracks[i].itsClusterSizes(), + daughterTracks[i].tpcNClsFindable(), daughterTracks[i].tpcNClsFindableMinusFound(), daughterTracks[i].tpcNClsFindableMinusCrossedRows(), + daughterTracks[i].trdPattern(), daughterTracks[i].tpcChi2NCl(), daughterTracks[i].tofChi2(), + daughterTracks[i].tpcSignal(), daughterTracks[i].tofExpMom(), + // PID + daughterTracks[i].tpcNSigmaPr(), daughterTracks[i].tpcNSigmaPi(), daughterTracks[i].tpcNSigmaDe(), + tofNSigmaTrack); + } + + // save reduced decay3body table + reducedDecay3Bodys(reducedCollisions.lastIndex(), reducedFullTracksPIDIU.lastIndex() - 2, reducedFullTracksPIDIU.lastIndex() - 1, reducedFullTracksPIDIU.lastIndex()); + } + + registry.fill(HIST("hEventCounter"), 3.5, reducedCollisions.lastIndex() + 1); + } +}; + +struct reduced3bodyInitializer { + Spawns reducedTracksIU; + void init(InitContext const&) {} +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), + }; +}