diff --git a/PWGLF/DataModel/LFStrangenessPIDTables.h b/PWGLF/DataModel/LFStrangenessPIDTables.h index 0b7e95ce302..b157097a605 100644 --- a/PWGLF/DataModel/LFStrangenessPIDTables.h +++ b/PWGLF/DataModel/LFStrangenessPIDTables.h @@ -15,36 +15,156 @@ #ifndef PWGLF_DATAMODEL_LFSTRANGENESSPIDTABLES_H_ #define PWGLF_DATAMODEL_LFSTRANGENESSPIDTABLES_H_ -#include -#include "Framework/AnalysisDataModel.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" + #include "Common/Core/RecoDecay.h" + #include "CommonConstants/PhysicsConstants.h" +#include "Framework/AnalysisDataModel.h" + +#include namespace o2::aod { namespace dautrack { -// ==== TPC INFORMATION === -DECLARE_SOA_COLUMN(TPCSignal, tpcSignal, float); //! track TPC signal +// ==== define packing helpers === +namespace packing +{ +// define variables for packing +static constexpr int nbins = (1 << 8 * sizeof(int8_t)) - 2; +static constexpr int8_t overflowBin = nbins >> 1; +static constexpr int8_t underflowBin = -(nbins >> 1); +static constexpr float binned_max = 6.35; +static constexpr float binned_min = -6.35; +static constexpr float bin_width = (binned_max - binned_min) / nbins; +static constexpr float underflow_return = -100.0f; +static constexpr float overflow_return = +100.0f; + +// define helper function to do packing +int8_t packInInt8(float nSigma) +{ + // calculate + if (nSigma <= binned_min) + return underflowBin; + if (nSigma >= binned_max) + return overflowBin; + if (nSigma >= 0) { + return static_cast((nSigma / bin_width) + 0.5f); + } + // automatic: this is the case in which nSigma < 0 + return static_cast((nSigma / bin_width) - 0.5f); +} + +// define helper function to do unpacking +float unpackInt8(int8_t nSigma) +{ + if (nSigma == underflowBin) { + return underflow_return; + } + if (nSigma == overflowBin) { + return overflow_return; + } + return bin_width * nSigma; +} + +} // namespace packing +} // namespace dautrack + +namespace dautrack_legacy +{ +// ==== LEGACY TPC INFORMATION (full size tables) === DECLARE_SOA_COLUMN(TPCNSigmaEl, tpcNSigmaEl, float); //! Nsigma proton DECLARE_SOA_COLUMN(TPCNSigmaPi, tpcNSigmaPi, float); //! Nsigma proton DECLARE_SOA_COLUMN(TPCNSigmaKa, tpcNSigmaKa, float); //! Nsigma proton DECLARE_SOA_COLUMN(TPCNSigmaPr, tpcNSigmaPr, float); //! Nsigma proton DECLARE_SOA_COLUMN(TPCNSigmaHe, tpcNSigmaHe, float); //! Nsigma proton +} // namespace dautrack_legacy + +namespace dautrack +{ +// ==== COMPACT TPC INFORMATION (full size tables) === +DECLARE_SOA_COLUMN(TPCSignal, tpcSignal, float); //! track TPC signal +DECLARE_SOA_COLUMN(PackedTPCNSigmaEl, packedTpcNSigmaEl, int8_t); //! Nsigma proton +DECLARE_SOA_COLUMN(PackedTPCNSigmaPi, packedTpcNSigmaPi, int8_t); //! Nsigma proton +DECLARE_SOA_COLUMN(PackedTPCNSigmaKa, packedTpcNSigmaKa, int8_t); //! Nsigma proton +DECLARE_SOA_COLUMN(PackedTPCNSigmaPr, packedTpcNSigmaPr, int8_t); //! Nsigma proton + +DECLARE_SOA_DYNAMIC_COLUMN(TPCNSigmaEl, tpcNSigmaEl, //! unpacked TPC nsigma + [](int8_t nsigma_packed) -> float { return o2::aod::dautrack::packing::unpackInt8(nsigma_packed); }); +DECLARE_SOA_DYNAMIC_COLUMN(TPCNSigmaPi, tpcNSigmaPi, //! unpacked TPC nsigma + [](int8_t nsigma_packed) -> float { return o2::aod::dautrack::packing::unpackInt8(nsigma_packed); }); +DECLARE_SOA_DYNAMIC_COLUMN(TPCNSigmaKa, tpcNSigmaKa, //! unpacked TPC nsigma + [](int8_t nsigma_packed) -> float { return o2::aod::dautrack::packing::unpackInt8(nsigma_packed); }); +DECLARE_SOA_DYNAMIC_COLUMN(TPCNSigmaPr, tpcNSigmaPr, //! unpacked TPC nsigma + [](int8_t nsigma_packed) -> float { return o2::aod::dautrack::packing::unpackInt8(nsigma_packed); }); // ==== TOF INFORMATION === +DECLARE_SOA_INDEX_COLUMN(DauTrackExtra, dauTrackExtra); //! point to daughter this TOF info belongs to +DECLARE_SOA_INDEX_COLUMN(StraCollision, straCollision); //! point to collision associated with this track (not the V0/Casc) DECLARE_SOA_COLUMN(TOFSignal, tofSignal, float); //! track TOF signal -DECLARE_SOA_COLUMN(TOFEvTime, tofEvTime, float); //! track TOF signal -DECLARE_SOA_COLUMN(Length, length, float); //! track TOF signal +DECLARE_SOA_COLUMN(TOFEvTime, tofEvTime, float); //! event time +DECLARE_SOA_COLUMN(Length, length, float); //! track length (to assigned PV) +DECLARE_SOA_COLUMN(TOFExpMom, tofExpMom, float); //! tof Exp Mom (to assigned PV) + +// dynamics with expected times +DECLARE_SOA_DYNAMIC_COLUMN(TOFExpTimeEl, tofExpTimeEl, //! Expected time for the track to reach the TOF under the electron hypothesis + [](float length, float tofExpMom) -> float { + constexpr float massSquared = o2::constants::physics::MassElectron * o2::constants::physics::MassElectron; + return o2::framework::pid::tof::MassToExpTime(tofExpMom, length, massSquared); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(TOFExpTimePi, tofExpTimePi, //! Expected time for the track to reach the TOF under the pion hypothesis + [](float length, float tofExpMom) -> float { + constexpr float massSquared = o2::constants::physics::MassPionCharged * o2::constants::physics::MassPionCharged; + return o2::framework::pid::tof::MassToExpTime(tofExpMom, length, massSquared); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(TOFExpTimeKa, tofExpTimeKa, //! Expected time for the track to reach the TOF under the kaon hypothesis + [](float length, float tofExpMom) -> float { + constexpr float massSquared = o2::constants::physics::MassKaonCharged * o2::constants::physics::MassKaonCharged; + return o2::framework::pid::tof::MassToExpTime(tofExpMom, length, massSquared); + }); + +DECLARE_SOA_DYNAMIC_COLUMN(TOFExpTimePr, tofExpTimePr, //! Expected time for the track to reach the TOF under the proton hypothesis + [](float length, float tofExpMom) -> float { + constexpr float massSquared = o2::constants::physics::MassProton * o2::constants::physics::MassProton; + return o2::framework::pid::tof::MassToExpTime(tofExpMom, length, massSquared); + }); + } // namespace dautrack -DECLARE_SOA_TABLE(DauTrackTPCPIDs, "AOD", "DAUTRACKTPCPID", // nsigma table (for analysis) - dautrack::TPCSignal, dautrack::TPCNSigmaEl, - dautrack::TPCNSigmaPi, dautrack::TPCNSigmaKa, - dautrack::TPCNSigmaPr, dautrack::TPCNSigmaHe); -DECLARE_SOA_TABLE(DauTrackTOFPIDs, "AOD", "DAUTRACKTOFPID", // raw table (for posterior TOF calculation) +DECLARE_SOA_TABLE(DauTrackTPCPIDs_000, "AOD", "DAUTRACKTPCPID", // nsigma table (for analysis) + dautrack::TPCSignal, dautrack_legacy::TPCNSigmaEl, + dautrack_legacy::TPCNSigmaPi, dautrack_legacy::TPCNSigmaKa, + dautrack_legacy::TPCNSigmaPr, dautrack_legacy::TPCNSigmaHe); + +DECLARE_SOA_TABLE_VERSIONED(DauTrackTPCPIDs_001, "AOD", "DAUTRACKTPCPID", 1, // nsigma table (for analysis) + dautrack::TPCSignal, + dautrack::PackedTPCNSigmaEl, dautrack::PackedTPCNSigmaPi, + dautrack::PackedTPCNSigmaKa, dautrack::PackedTPCNSigmaPr, + dautrack::TPCNSigmaEl, + dautrack::TPCNSigmaPi, + dautrack::TPCNSigmaKa, + dautrack::TPCNSigmaPr); + +using DauTrackTPCPIDs = DauTrackTPCPIDs_001; // second gen: packed Nsigma, no He + +DECLARE_SOA_TABLE(DauTrackTOFPIDs_000, "AOD", "DAUTRACKTOFPID", // raw table (for posterior TOF calculation) dautrack::TOFSignal, dautrack::TOFEvTime, dautrack::Length); +DECLARE_SOA_TABLE_VERSIONED(DauTrackTOFPIDs_001, "AOD", "DAUTRACKTOFPID", 1, // raw table (for posterior TOF calculation) + o2::soa::Index<>, + dautrack::StraCollisionId, dautrack::DauTrackExtraId, + dautrack::TOFSignal, dautrack::TOFEvTime, + dautrack::Length, dautrack::TOFExpMom, + dautrack::TOFExpTimeEl, + dautrack::TOFExpTimePi, + dautrack::TOFExpTimeKa, + dautrack::TOFExpTimePr); + +using DauTrackTOFPIDs = DauTrackTOFPIDs_001; // second gen: with collision Id, with TOFExpMom + namespace v0data { // define constants for NSigma operation diff --git a/PWGLF/DataModel/LFStrangenessTables.h b/PWGLF/DataModel/LFStrangenessTables.h index ca58973dff8..6da514fa64e 100644 --- a/PWGLF/DataModel/LFStrangenessTables.h +++ b/PWGLF/DataModel/LFStrangenessTables.h @@ -63,6 +63,9 @@ DECLARE_SOA_DYNAMIC_COLUMN(EnergyCommonZNA, energyCommonZNA, //! get the total s [](float value) -> float { return value; }); DECLARE_SOA_DYNAMIC_COLUMN(EnergyCommonZNC, energyCommonZNC, //! get the total sum of the ZN A amplitudes [](float value) -> float { return value; }); + +// event time +DECLARE_SOA_COLUMN(EventTime, eventTime, float); //! event time (FT0, TOF) for TOF PID (stored once per event) } // namespace stracollision //______________________________________________________ @@ -308,6 +311,8 @@ DECLARE_SOA_TABLE(StraStamps_000, "AOD", "STRASTAMPS", //! information for ID-in bc::RunNumber, timestamp::Timestamp); DECLARE_SOA_TABLE_VERSIONED(StraStamps_001, "AOD", "STRASTAMPS", 1, //! information for ID-ing mag field if needed bc::RunNumber, timestamp::Timestamp, bc::GlobalBC); +DECLARE_SOA_TABLE(StraEvTimes, "AOD", "STRAEVTIMES", //! event time (FT0, TOF) + stracollision::EventTime); using StraRawCents = StraRawCents_004; using StraCents = StraCents_001; diff --git a/PWGLF/TableProducer/Strangeness/Converters/CMakeLists.txt b/PWGLF/TableProducer/Strangeness/Converters/CMakeLists.txt index 842137dafe1..e0ebdb862f9 100644 --- a/PWGLF/TableProducer/Strangeness/Converters/CMakeLists.txt +++ b/PWGLF/TableProducer/Strangeness/Converters/CMakeLists.txt @@ -9,11 +9,21 @@ # granted to it by virtue of its status as an Intergovernmental Organization # or submit itself to any jurisdiction. +o2physics_add_dpl_workflow(stradautrackstpcpidconverter + SOURCES stradautrackstpcpidconverter.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(stradautrackstofpidconverter SOURCES stradautrackstofpidconverter.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(stradautrackstofpidconverter2 + SOURCES stradautrackstofpidconverter2.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(stradautracksextraconverter2 SOURCES stradautracksextraconverter2.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore diff --git a/PWGLF/TableProducer/Strangeness/Converters/stradautrackstofpidconverter.cxx b/PWGLF/TableProducer/Strangeness/Converters/stradautrackstofpidconverter.cxx index 8731939ca53..26cefdf2485 100644 --- a/PWGLF/TableProducer/Strangeness/Converters/stradautrackstofpidconverter.cxx +++ b/PWGLF/TableProducer/Strangeness/Converters/stradautrackstofpidconverter.cxx @@ -53,7 +53,7 @@ struct stradautrackstofpidconverter { lTOFEvTimes[casc.bachTrackExtraId()] = casc.bachTOFEventTime(); } for (int ii = 0; ii < dauTracks.size(); ii++) { - dautracktofpids(lTOFSignals[ii], lTOFEvTimes[ii], lLengths[ii]); + dautracktofpids(-1, -1, lTOFSignals[ii], lTOFEvTimes[ii], lLengths[ii], 0.0f); } } }; diff --git a/PWGLF/TableProducer/Strangeness/Converters/stradautrackstofpidconverter2.cxx b/PWGLF/TableProducer/Strangeness/Converters/stradautrackstofpidconverter2.cxx new file mode 100644 index 00000000000..371c5c133d3 --- /dev/null +++ b/PWGLF/TableProducer/Strangeness/Converters/stradautrackstofpidconverter2.cxx @@ -0,0 +1,66 @@ +// 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. +#include "PWGLF/DataModel/LFStrangenessPIDTables.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" + +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" + +using namespace o2; +using namespace o2::framework; + +// converts DauTrackTOFPIDs_000 to _001 +struct stradautrackstofpidconverter2 { + Produces dautracktofpids; + Produces straEvTimes; + + void process(aod::StraCollisions const& collisions, soa::Join const& dauTracks, soa::Join const& v0s) + { + // create new TOFPIDs + for (int ii = 0; ii < dauTracks.size(); ii++) { + auto dauTrack = dauTracks.rawIteratorAt(ii); + dautracktofpids(-1, -1, dauTrack.tofSignal(), dauTrack.tofEvTime(), dauTrack.length(), 0.0f); + } + + // fill EvTimes (created simultaneously, so done in the same converter) + // recover as much as possible from the previous format + std::vector collisionEventTime(collisions.size(), 0.0); + std::vector collisionNtracks(collisions.size(), 0); + + for (const auto& v0 : v0s) { + auto posTrackTOF = dauTracks.rawIteratorAt(v0.posTrackExtraId()); + auto negTrackTOF = dauTracks.rawIteratorAt(v0.negTrackExtraId()); + if (posTrackTOF.hasTOF()) { + collisionEventTime[v0.straCollisionId()] += posTrackTOF.tofEvTime(); + collisionNtracks[v0.straCollisionId()]++; + } + if (negTrackTOF.hasTOF()) { + collisionEventTime[v0.straCollisionId()] += negTrackTOF.tofEvTime(); + collisionNtracks[v0.straCollisionId()]++; + } + } + for (const auto& collision : collisions) { + if (collisionNtracks[collision.globalIndex()] > 0) { + collisionEventTime[collision.globalIndex()] /= static_cast(collisionNtracks[collision.globalIndex()]); + } else { + collisionEventTime[collision.globalIndex()] = -1e+6; // undefined + } + straEvTimes(collisionEventTime[collision.globalIndex()]); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +} diff --git a/PWGLF/TableProducer/Strangeness/Converters/stradautrackstpcpidconverter.cxx b/PWGLF/TableProducer/Strangeness/Converters/stradautrackstpcpidconverter.cxx new file mode 100644 index 00000000000..627872555bc --- /dev/null +++ b/PWGLF/TableProducer/Strangeness/Converters/stradautrackstpcpidconverter.cxx @@ -0,0 +1,42 @@ +// 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. +#include "PWGLF/DataModel/LFStrangenessPIDTables.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" + +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" + +using namespace o2; +using namespace o2::framework; + +// converts DauTrackTOFPIDs_000 to _001 +struct stradautrackstpcpidconverter { + Produces dautrackpcpids; + + void process(aod::DauTrackTPCPIDs_000 const& v000s) + { + for (int ii = 0; ii < v000s.size(); ii++) { + auto dauTrackTPCPID = v000s.rawIteratorAt(ii); + dautrackpcpids(dauTrackTPCPID.tpcSignal(), + aod::dautrack::packing::packInInt8(dauTrackTPCPID.tpcNSigmaEl()), + aod::dautrack::packing::packInInt8(dauTrackTPCPID.tpcNSigmaPi()), + aod::dautrack::packing::packInInt8(dauTrackTPCPID.tpcNSigmaKa()), + aod::dautrack::packing::packInInt8(dauTrackTPCPID.tpcNSigmaPr())); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +} diff --git a/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx b/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx index 8140448576a..5a563a10a3d 100644 --- a/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx +++ b/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx @@ -72,74 +72,77 @@ using BCsWithTimestampsAndRun2Infos = soa::Join strangeColl; // characterises collisions - Produces strangeCollLabels; // characterises collisions - Produces strangeMCColl; // characterises collisions / MC - Produces strangeMCMults; // characterises collisions / MC mults - Produces strangeCents; // characterises collisions / centrality in Run 3 - Produces strangeCentsRun2; // characterises collisions / centrality in Run 2 - Produces strangeEvSels; // characterises collisions / centrality / sel8 selection in Run 3 - Produces strangeEvSelsRun2; // characterises collisions / centrality / sel8 selection in Run 2 - Produces strangeStamps; // provides timestamps, run numbers - Produces v0collref; // references collisions from V0s - Produces casccollref; // references collisions from cascades - Produces kfcasccollref; // references collisions from KF cascades - Produces tracasccollref; // references collisions from tracked cascades - - //__________________________________________________ - // track extra references - Produces dauTrackExtras; // daughter track detector properties - Produces dauTrackMCIds; // daughter track MC Particle ID - Produces dauTrackTPCPIDs; // daughter track TPC PID - Produces dauTrackTOFPIDs; // daughter track TOF PID - Produces v0Extras; // references DauTracks from V0s - Produces cascExtras; // references DauTracks from cascades - Produces straTrackExtras; // references DauTracks from tracked cascades - - //__________________________________________________ - // cascade interlinks - Produces cascToTraRefs; // cascades -> tracked - Produces cascToKFRefs; // cascades -> KF - Produces traToCascRefs; // tracked -> cascades - Produces kfToCascRefs; // KF -> cascades - - //__________________________________________________ - // mother information - Produces v0mothers; // V0 mother references - Produces cascmothers; // casc mother references - Produces motherMCParts; // mc particles for mothers - - //__________________________________________________ - // UPC specific information - Produces zdcNeutrons; // Primary neutrons within ZDC acceptance - Produces zdcNeutronsMCCollRefs; // references collisions from ZDCNeutrons - - //__________________________________________________ - // Q-vectors - Produces StraFT0AQVs; // FT0A Q-vector - Produces StraFT0CQVs; // FT0C Q-vector - Produces StraFT0MQVs; // FT0M Q-vector - Produces StraFV0AQVs; // FV0A Q-vector - Produces StraTPCQVs; // TPC Q-vector - Produces StraFT0CQVsEv; // events used to compute FT0C Q-vector (LF) - Produces StraZDCSP; // ZDC Sums and Products - - //__________________________________________________ - // Generated binned data - // this is a hack while the system does not do better - Produces geK0Short; - Produces geLambda; - Produces geAntiLambda; - Produces geXiMinus; - Produces geXiPlus; - Produces geOmegaMinus; - Produces geOmegaPlus; - - //__________________________________________________ - // Debug - Produces straOrigin; + struct : ProducesGroup { + //__________________________________________________ + // fundamental building blocks of derived data + Produces strangeColl; // characterises collisions + Produces strangeCollLabels; // characterises collisions + Produces strangeMCColl; // characterises collisions / MC + Produces strangeMCMults; // characterises collisions / MC mults + Produces strangeCents; // characterises collisions / centrality in Run 3 + Produces strangeCentsRun2; // characterises collisions / centrality in Run 2 + Produces strangeEvSels; // characterises collisions / centrality / sel8 selection in Run 3 + Produces strangeEvSelsRun2; // characterises collisions / centrality / sel8 selection in Run 2 + Produces strangeStamps; // provides timestamps, run numbers + Produces straEvTimes; // provides event times (FT0, TOF) + Produces v0collref; // references collisions from V0s + Produces casccollref; // references collisions from cascades + Produces kfcasccollref; // references collisions from KF cascades + Produces tracasccollref; // references collisions from tracked cascades + + //__________________________________________________ + // track extra references + Produces dauTrackExtras; // daughter track detector properties + Produces dauTrackMCIds; // daughter track MC Particle ID + Produces dauTrackTPCPIDs; // daughter track TPC PID + Produces dauTrackTOFPIDs; // daughter track TOF PID + Produces v0Extras; // references DauTracks from V0s + Produces cascExtras; // references DauTracks from cascades + Produces straTrackExtras; // references DauTracks from tracked cascades + + //__________________________________________________ + // cascade interlinks + Produces cascToTraRefs; // cascades -> tracked + Produces cascToKFRefs; // cascades -> KF + Produces traToCascRefs; // tracked -> cascades + Produces kfToCascRefs; // KF -> cascades + + //__________________________________________________ + // mother information + Produces v0mothers; // V0 mother references + Produces cascmothers; // casc mother references + Produces motherMCParts; // mc particles for mothers + + //__________________________________________________ + // UPC specific information + Produces zdcNeutrons; // Primary neutrons within ZDC acceptance + Produces zdcNeutronsMCCollRefs; // references collisions from ZDCNeutrons + + //__________________________________________________ + // Q-vectors + Produces StraFT0AQVs; // FT0A Q-vector + Produces StraFT0CQVs; // FT0C Q-vector + Produces StraFT0MQVs; // FT0M Q-vector + Produces StraFV0AQVs; // FV0A Q-vector + Produces StraTPCQVs; // TPC Q-vector + Produces StraFT0CQVsEv; // events used to compute FT0C Q-vector (LF) + Produces StraZDCSP; // ZDC Sums and Products + + //__________________________________________________ + // Generated binned data + // this is a hack while the system does not do better + Produces geK0Short; + Produces geLambda; + Produces geAntiLambda; + Produces geXiMinus; + Produces geXiPlus; + Produces geOmegaMinus; + Produces geOmegaPlus; + + //__________________________________________________ + // Debug + Produces straOrigin; + } products; // histogram registry for bookkeeping HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -482,10 +485,10 @@ struct strangederivedbuilder { // +-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+ // fill collision tables if (strange || fillEmptyCollisions) { - strangeStamps(bc.runNumber(), bc.timestamp(), bc.globalBC()); - strangeColl(collision.posX(), collision.posY(), collision.posZ()); + products.strangeStamps(bc.runNumber(), bc.timestamp(), bc.globalBC()); + products.strangeColl(collision.posX(), collision.posY(), collision.posZ()); if constexpr (requires { collision.mcCollisionId(); }) { // check if MC information is available and if so fill labels - strangeCollLabels(collision.mcCollisionId()); + products.strangeCollLabels(collision.mcCollisionId()); } if constexpr (requires { collision.centFT0C(); }) { // check if we are in Run 3 @@ -495,89 +498,111 @@ struct strangederivedbuilder { centrality = hRawCentrality->GetBinContent(hRawCentrality->FindBin(collision.multFT0C())); } - strangeCents(collision.centFT0M(), collision.centFT0A(), - centrality, collision.centFV0A(), collision.centFT0CVariant1(), - collision.centMFT(), collision.centNGlobal()); - strangeEvSels(collision.sel8(), collision.selection_raw(), - collision.multFT0A() * static_cast(fillTruncationOptions.fillRawFT0A), - collision.multFT0C() * static_cast(fillTruncationOptions.fillRawFT0C), - collision.multFV0A() * static_cast(fillTruncationOptions.fillRawFV0A), - collision.multFDDA() * static_cast(fillTruncationOptions.fillRawFDDA), - collision.multFDDC() * static_cast(fillTruncationOptions.fillRawFDDC), - collision.multNTracksPVeta1() * static_cast(fillTruncationOptions.fillRawNTracksEta1), - collision.multPVTotalContributors() * static_cast(fillTruncationOptions.fillRawNTracksForCorrelation), - collision.multNTracksGlobal() * static_cast(fillTruncationOptions.fillRawNTracksForCorrelation), - collision.multNTracksITSTPC() * static_cast(fillTruncationOptions.fillRawNTracksForCorrelation), - collision.multAllTracksTPCOnly() * static_cast(fillTruncationOptions.fillRawNTracksForCorrelation), - collision.multAllTracksITSTPC() * static_cast(fillTruncationOptions.fillRawNTracksForCorrelation), - collision.multZNA() * static_cast(fillTruncationOptions.fillRawZDC), - collision.multZNC() * static_cast(fillTruncationOptions.fillRawZDC), - collision.multZEM1() * static_cast(fillTruncationOptions.fillRawZDC), - collision.multZEM2() * static_cast(fillTruncationOptions.fillRawZDC), - collision.multZPA() * static_cast(fillTruncationOptions.fillRawZDC), - collision.multZPC() * static_cast(fillTruncationOptions.fillRawZDC), - collision.trackOccupancyInTimeRange(), - collision.ft0cOccupancyInTimeRange(), - // UPC info - gapSide, - totalFT0AmplitudeA, totalFT0AmplitudeC, totalFV0AmplitudeA, - totalFDDAmplitudeA, totalFDDAmplitudeC, - energyCommonZNA, energyCommonZNC, - // Collision flags - collision.flags(), - collision.alias_raw(), - collision.rct_raw()); + products.strangeCents(collision.centFT0M(), collision.centFT0A(), + centrality, collision.centFV0A(), collision.centFT0CVariant1(), + collision.centMFT(), collision.centNGlobal()); + products.strangeEvSels(collision.sel8(), collision.selection_raw(), + collision.multFT0A() * static_cast(fillTruncationOptions.fillRawFT0A), + collision.multFT0C() * static_cast(fillTruncationOptions.fillRawFT0C), + collision.multFV0A() * static_cast(fillTruncationOptions.fillRawFV0A), + collision.multFDDA() * static_cast(fillTruncationOptions.fillRawFDDA), + collision.multFDDC() * static_cast(fillTruncationOptions.fillRawFDDC), + collision.multNTracksPVeta1() * static_cast(fillTruncationOptions.fillRawNTracksEta1), + collision.multPVTotalContributors() * static_cast(fillTruncationOptions.fillRawNTracksForCorrelation), + collision.multNTracksGlobal() * static_cast(fillTruncationOptions.fillRawNTracksForCorrelation), + collision.multNTracksITSTPC() * static_cast(fillTruncationOptions.fillRawNTracksForCorrelation), + collision.multAllTracksTPCOnly() * static_cast(fillTruncationOptions.fillRawNTracksForCorrelation), + collision.multAllTracksITSTPC() * static_cast(fillTruncationOptions.fillRawNTracksForCorrelation), + collision.multZNA() * static_cast(fillTruncationOptions.fillRawZDC), + collision.multZNC() * static_cast(fillTruncationOptions.fillRawZDC), + collision.multZEM1() * static_cast(fillTruncationOptions.fillRawZDC), + collision.multZEM2() * static_cast(fillTruncationOptions.fillRawZDC), + collision.multZPA() * static_cast(fillTruncationOptions.fillRawZDC), + collision.multZPC() * static_cast(fillTruncationOptions.fillRawZDC), + collision.trackOccupancyInTimeRange(), + collision.ft0cOccupancyInTimeRange(), + // UPC info + gapSide, + totalFT0AmplitudeA, totalFT0AmplitudeC, totalFV0AmplitudeA, + totalFDDAmplitudeA, totalFDDAmplitudeC, + energyCommonZNA, energyCommonZNC, + // Collision flags + collision.flags(), + collision.alias_raw(), + collision.rct_raw()); } else { // We are in Run 2 - strangeCentsRun2(collision.centRun2V0M(), collision.centRun2V0A(), - collision.centRun2SPDTracklets(), collision.centRun2SPDClusters()); - strangeEvSelsRun2(collision.sel8(), collision.sel7(), collision.selection_raw(), - collision.multFT0A() * static_cast(fillTruncationOptions.fillRawFT0A), - collision.multFT0C() * static_cast(fillTruncationOptions.fillRawFT0C), - collision.multFV0A() * static_cast(fillTruncationOptions.fillRawFV0A), - collision.multFV0C() * static_cast(fillTruncationOptions.fillRawFV0C), - collision.multFDDA() * static_cast(fillTruncationOptions.fillRawFDDA), - collision.multFDDC() * static_cast(fillTruncationOptions.fillRawFDDC), - bc.spdClustersL0() * static_cast(fillTruncationOptions.fillRawSPDclsL0Run2), - bc.spdClustersL1() * static_cast(fillTruncationOptions.fillRawSPDclsL1Run2), - collision.multNTracksPVeta1() * static_cast(fillTruncationOptions.fillRawNTracksEta1), - collision.multTracklets() * static_cast(fillTruncationOptions.fillRawTrackletsRun2), - -1, /* dummy number of PV contribs total while waiting for the multiplicity task to produce it */ - -1, /* dummy global track multiplicities while waiting for the multiplicity task to produce it */ - -1, /* dummy track multiplicities, PV contribs, no eta cut while waiting for the multiplicity task to produce it */ - -1, /* dummy TPConly track multiplicities, all, no eta cut while waiting for the multiplicity task to produce it */ - -1, /* dummy ITSTPC track multiplicities, all, no eta cut waiting for the multiplicity task to produce it */ - collision.multZNA() * static_cast(fillTruncationOptions.fillRawZDC), - collision.multZNC() * static_cast(fillTruncationOptions.fillRawZDC), - collision.multZEM1() * static_cast(fillTruncationOptions.fillRawZDC), - collision.multZEM2() * static_cast(fillTruncationOptions.fillRawZDC), - collision.multZPA() * static_cast(fillTruncationOptions.fillRawZDC), - collision.multZPC() * static_cast(fillTruncationOptions.fillRawZDC), - collision.alias_raw()); + products.strangeCentsRun2(collision.centRun2V0M(), collision.centRun2V0A(), + collision.centRun2SPDTracklets(), collision.centRun2SPDClusters()); + products.strangeEvSelsRun2(collision.sel8(), collision.sel7(), collision.selection_raw(), + collision.multFT0A() * static_cast(fillTruncationOptions.fillRawFT0A), + collision.multFT0C() * static_cast(fillTruncationOptions.fillRawFT0C), + collision.multFV0A() * static_cast(fillTruncationOptions.fillRawFV0A), + collision.multFV0C() * static_cast(fillTruncationOptions.fillRawFV0C), + collision.multFDDA() * static_cast(fillTruncationOptions.fillRawFDDA), + collision.multFDDC() * static_cast(fillTruncationOptions.fillRawFDDC), + bc.spdClustersL0() * static_cast(fillTruncationOptions.fillRawSPDclsL0Run2), + bc.spdClustersL1() * static_cast(fillTruncationOptions.fillRawSPDclsL1Run2), + collision.multNTracksPVeta1() * static_cast(fillTruncationOptions.fillRawNTracksEta1), + collision.multTracklets() * static_cast(fillTruncationOptions.fillRawTrackletsRun2), + -1, /* dummy number of PV contribs total while waiting for the multiplicity task to produce it */ + -1, /* dummy global track multiplicities while waiting for the multiplicity task to produce it */ + -1, /* dummy track multiplicities, PV contribs, no eta cut while waiting for the multiplicity task to produce it */ + -1, /* dummy TPConly track multiplicities, all, no eta cut while waiting for the multiplicity task to produce it */ + -1, /* dummy ITSTPC track multiplicities, all, no eta cut waiting for the multiplicity task to produce it */ + collision.multZNA() * static_cast(fillTruncationOptions.fillRawZDC), + collision.multZNC() * static_cast(fillTruncationOptions.fillRawZDC), + collision.multZEM1() * static_cast(fillTruncationOptions.fillRawZDC), + collision.multZEM2() * static_cast(fillTruncationOptions.fillRawZDC), + collision.multZPA() * static_cast(fillTruncationOptions.fillRawZDC), + collision.multZPC() * static_cast(fillTruncationOptions.fillRawZDC), + collision.alias_raw()); } } for (const auto& v0 : V0Table_thisColl) - V0CollIndices[v0.globalIndex()] = strangeColl.lastIndex(); + V0CollIndices[v0.globalIndex()] = products.strangeColl.lastIndex(); for (const auto& casc : CascTable_thisColl) - CascadeCollIndices[casc.globalIndex()] = strangeColl.lastIndex(); + CascadeCollIndices[casc.globalIndex()] = products.strangeColl.lastIndex(); for (const auto& casc : KFCascTable_thisColl) - KFCascadeCollIndices[casc.globalIndex()] = strangeColl.lastIndex(); + KFCascadeCollIndices[casc.globalIndex()] = products.strangeColl.lastIndex(); for (const auto& casc : TraCascTable_thisColl) - TraCascadeCollIndices[casc.globalIndex()] = strangeColl.lastIndex(); + TraCascadeCollIndices[casc.globalIndex()] = products.strangeColl.lastIndex(); } // +-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+-<*>-+ // populate references, including those that might not be assigned for (const auto& v0 : V0s) { - v0collref(V0CollIndices[v0.globalIndex()]); + products.v0collref(V0CollIndices[v0.globalIndex()]); } for (const auto& casc : Cascades) { - casccollref(CascadeCollIndices[casc.globalIndex()]); + products.casccollref(CascadeCollIndices[casc.globalIndex()]); } for (const auto& casc : KFCascades) { - kfcasccollref(KFCascadeCollIndices[casc.globalIndex()]); + products.kfcasccollref(KFCascadeCollIndices[casc.globalIndex()]); } for (const auto& casc : TraCascades) { - tracasccollref(TraCascadeCollIndices[casc.globalIndex()]); + products.tracasccollref(TraCascadeCollIndices[casc.globalIndex()]); + } + } + + // helper function to estimate collision time + template + void populateEventTimes(coll const& collisions, TTracks const& tracks) + { + std::vector collisionEventTime(collisions.size(), 0.0); + std::vector collisionNtracks(collisions.size(), 0); + for (const auto& track : tracks) { + if (track.hasTOF() && track.collisionId() >= 0) { + collisionEventTime[track.collisionId()] += track.tofEvTime(); + collisionNtracks[track.collisionId()]++; + } + } + for (const auto& collision : collisions) { + if (collisionNtracks[collision.globalIndex()] > 0) { + collisionEventTime[collision.globalIndex()] /= static_cast(collisionNtracks[collision.globalIndex()]); + } else { + collisionEventTime[collision.globalIndex()] = -1e+6; // undefined + } + products.straEvTimes(collisionEventTime[collision.globalIndex()]); } } @@ -611,13 +636,13 @@ struct strangederivedbuilder { totalMult++; } - strangeMCColl(mccollision.posX(), mccollision.posY(), mccollision.posZ(), - mccollision.impactParameter(), mccollision.eventPlaneAngle(), mccollision.generatorsID()); - strangeMCMults(mccollision.multMCFT0A(), mccollision.multMCFT0C(), - mccollision.multMCNParticlesEta05(), - mccollision.multMCNParticlesEta08(), - mccollision.multMCNParticlesEta10(), - totalMult); + products.strangeMCColl(mccollision.posX(), mccollision.posY(), mccollision.posZ(), + mccollision.impactParameter(), mccollision.eventPlaneAngle(), mccollision.generatorsID()); + products.strangeMCMults(mccollision.multMCFT0A(), mccollision.multMCFT0C(), + mccollision.multMCNParticlesEta05(), + mccollision.multMCNParticlesEta08(), + mccollision.multMCNParticlesEta10(), + totalMult); } } @@ -680,21 +705,21 @@ struct strangederivedbuilder { for (auto const& v0 : V0s) { auto const& posTrack = v0.posTrack_as(); auto const& negTrack = v0.negTrack_as(); - v0Extras(trackMap[posTrack.globalIndex()], - trackMap[negTrack.globalIndex()]); // joinable with V0Datas + products.v0Extras(trackMap[posTrack.globalIndex()], + trackMap[negTrack.globalIndex()]); // joinable with V0Datas } //__________________________________________________ // circle back and populate actual DauTrackExtra table for (auto const& tr : tracksExtra) { if (trackMap[tr.globalIndex()] >= 0) { - dauTrackExtras(tr.itsChi2NCl(), - tr.tpcChi2NCl(), - tr.detectorMap(), - tr.itsClusterSizes(), - tr.tpcNClsFindable(), - tr.tpcNClsFindableMinusFound(), - tr.tpcNClsFindableMinusCrossedRows(), - tr.tpcNClsShared()); + products.dauTrackExtras(tr.itsChi2NCl(), + tr.tpcChi2NCl(), + tr.detectorMap(), + tr.itsClusterSizes(), + tr.tpcNClsFindable(), + tr.tpcNClsFindableMinusFound(), + tr.tpcNClsFindableMinusCrossedRows(), + tr.tpcNClsShared()); } } // done! @@ -760,8 +785,8 @@ struct strangederivedbuilder { for (auto const& v0 : V0s) { auto const& posTrack = v0.template posTrack_as(); auto const& negTrack = v0.template negTrack_as(); - v0Extras(trackMap[posTrack.globalIndex()], - trackMap[negTrack.globalIndex()]); // joinable with V0Datas + products.v0Extras(trackMap[posTrack.globalIndex()], + trackMap[negTrack.globalIndex()]); // joinable with V0Datas } //__________________________________________________ // populate track references @@ -769,61 +794,60 @@ struct strangederivedbuilder { auto bachTrack = casc.template bachelor_as(); auto posTrack = casc.template posTrack_as(); auto negTrack = casc.template negTrack_as(); - cascExtras(trackMap[posTrack.globalIndex()], - trackMap[negTrack.globalIndex()], - trackMap[bachTrack.globalIndex()]); // joinable with CascDatas + products.cascExtras(trackMap[posTrack.globalIndex()], + trackMap[negTrack.globalIndex()], + trackMap[bachTrack.globalIndex()]); // joinable with CascDatas } //__________________________________________________ // populate track references for (auto const& casc : TraCascades) { auto strangeTrack = casc.template strangeTrack_as(); - straTrackExtras(trackMap[strangeTrack.globalIndex()]); // joinable with TraCascDatas + products.straTrackExtras(trackMap[strangeTrack.globalIndex()]); // joinable with TraCascDatas } //__________________________________________________ // circle back and populate actual DauTrackExtra table for (auto const& tr : tracksExtra) { if (trackMap[tr.globalIndex()] >= 0) { - dauTrackExtras(tr.itsChi2NCl(), - tr.tpcChi2NCl(), - tr.detectorMap(), - tr.itsClusterSizes(), - tr.tpcNClsFindable(), - tr.tpcNClsFindableMinusFound(), - tr.tpcNClsFindableMinusCrossedRows(), - tr.tpcNClsShared()); + products.dauTrackExtras(tr.itsChi2NCl(), + tr.tpcChi2NCl(), + tr.detectorMap(), + tr.itsClusterSizes(), + tr.tpcNClsFindable(), + tr.tpcNClsFindableMinusFound(), + tr.tpcNClsFindableMinusCrossedRows(), + tr.tpcNClsShared()); // _________________________________________ // if the table has MC info if constexpr (requires { tr.mcParticle(); }) { // do your thing with the mcParticleIds only in case the table has the MC info - dauTrackMCIds(tr.mcParticleId()); // joinable with dauTrackExtras + products.dauTrackMCIds(tr.mcParticleId()); // joinable with dauTrackExtras } if constexpr (requires { tr.tpcNSigmaEl(); }) { - if (roundNSigmaVariables) { // round if requested - dauTrackTPCPIDs(tr.tpcSignal(), - roundToPrecision(tr.tpcNSigmaEl(), precisionNSigmas), - roundToPrecision(tr.tpcNSigmaPi(), precisionNSigmas), - roundToPrecision(tr.tpcNSigmaKa(), precisionNSigmas), - roundToPrecision(tr.tpcNSigmaPr(), precisionNSigmas), - roundToPrecision(tr.tpcNSigmaHe(), precisionNSigmas)); - } else { - dauTrackTPCPIDs(tr.tpcSignal(), tr.tpcNSigmaEl(), - tr.tpcNSigmaPi(), tr.tpcNSigmaKa(), - tr.tpcNSigmaPr(), tr.tpcNSigmaHe()); - } + products.dauTrackTPCPIDs(tr.tpcSignal(), + aod::dautrack::packing::packInInt8(tr.tpcNSigmaEl()), + aod::dautrack::packing::packInInt8(tr.tpcNSigmaPi()), + aod::dautrack::packing::packInInt8(tr.tpcNSigmaKa()), + aod::dautrack::packing::packInInt8(tr.tpcNSigmaPr())); // populate daughter-level TOF information - dauTrackTOFPIDs(tr.tofSignal(), tr.tofEvTime(), tr.length()); + if (tr.hasTOF()) { + products.dauTrackTOFPIDs(tr.collisionId(), products.dauTrackExtras.lastIndex(), tr.tofSignal(), tr.tofEvTime(), tr.length(), tr.tofExpMom()); + } } else { // populate with empty fully-compatible Nsigmas if no corresponding table available - dauTrackTPCPIDs(0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f); - dauTrackTOFPIDs(0.0f, 0.0f, 0.0f); + products.dauTrackTPCPIDs(0.0f, 0, 0, 0, 0); } } } // done! } + void processPopulateEventTimes(aod::Collisions const& collisions, soa::Join const& tracks) + { + populateEventTimes(collisions, tracks); + } + void processTrackExtras(aod::V0Datas const& V0s, aod::CascDatas const& Cascades, aod::KFCascDatas const& KFCascades, aod::TraCascDatas const& TraCascades, soa::Join const& tracksExtra, aod::V0s const&) { fillTrackExtras(V0s, Cascades, KFCascades, TraCascades, tracksExtra); @@ -872,23 +896,23 @@ struct strangederivedbuilder { // populate track references for (auto const& v0 : V0s) { if (v0.mcMotherParticleId() > -1) { - v0mothers(motherReference[v0.mcMotherParticleId()]); // joinable with V0Datas + products.v0mothers(motherReference[v0.mcMotherParticleId()]); // joinable with V0Datas } else { - v0mothers(-1); // joinable with V0Datas + products.v0mothers(-1); // joinable with V0Datas } } for (auto const& ca : Cascades) { if (ca.mcMotherParticleId() > -1) { - cascmothers(motherReference[ca.mcMotherParticleId()]); // joinable with CascDatas + products.cascmothers(motherReference[ca.mcMotherParticleId()]); // joinable with CascDatas } else { - cascmothers(-1); // joinable with CascDatas + products.cascmothers(-1); // joinable with CascDatas } } //__________________________________________________ // populate motherMCParticles for (auto const& tr : mcParticles) { if (motherReference[tr.globalIndex()] >= 0) { - motherMCParts(tr.px(), tr.py(), tr.pz(), tr.pdgCode(), tr.isPhysicalPrimary()); + products.motherMCParts(tr.px(), tr.py(), tr.pz(), tr.pdgCode(), tr.isPhysicalPrimary()); } } } @@ -904,7 +928,7 @@ struct strangederivedbuilder { auto cascade = c.cascade_as(); indexTracked = cascade.traCascDataId(); } - cascToTraRefs(indexTracked); + products.cascToTraRefs(indexTracked); } // Tracked to standard for (auto const& c : TraCascades) { @@ -913,7 +937,7 @@ struct strangederivedbuilder { auto cascade = c.cascade_as(); index = cascade.cascDataId(); } - traToCascRefs(index); + products.traToCascRefs(index); } } @@ -926,7 +950,7 @@ struct strangederivedbuilder { auto cascade = c.cascade_as(); indexKF = cascade.kfCascDataId(); } - cascToKFRefs(indexKF); + products.cascToKFRefs(indexKF); } // KF to standard for (auto const& c : KFCascades) { @@ -935,7 +959,7 @@ struct strangederivedbuilder { auto cascade = c.cascade_as(); index = cascade.cascDataId(); } - kfToCascRefs(index); + products.kfToCascRefs(index); } } @@ -1024,48 +1048,48 @@ struct strangederivedbuilder { } // at end of data frame // -> pack information from this DF into a generated histogram, once / DF - geK0Short(genK0Short); - geLambda(genLambda); - geAntiLambda(genAntiLambda); - geXiMinus(genXiMinus); - geXiPlus(genXiPlus); - geOmegaMinus(genOmegaMinus); - geOmegaPlus(genOmegaPlus); + products.geK0Short(genK0Short); + products.geLambda(genLambda); + products.geAntiLambda(genAntiLambda); + products.geXiMinus(genXiMinus); + products.geXiPlus(genXiPlus); + products.geOmegaMinus(genOmegaMinus); + products.geOmegaPlus(genOmegaPlus); } void processFT0AQVectors(soa::Join::iterator const& collision) { - StraFT0AQVs(collision.qvecFT0ARe(), collision.qvecFT0AIm(), collision.sumAmplFT0A()); + products.StraFT0AQVs(collision.qvecFT0ARe(), collision.qvecFT0AIm(), collision.sumAmplFT0A()); } void processFT0CQVectors(soa::Join::iterator const& collision) { - StraFT0CQVs(collision.qvecFT0CRe(), collision.qvecFT0CIm(), collision.sumAmplFT0C()); + products.StraFT0CQVs(collision.qvecFT0CRe(), collision.qvecFT0CIm(), collision.sumAmplFT0C()); } void processFT0CQVectorsLF(soa::Join::iterator const& collision) { - StraFT0CQVs(collision.qFT0C() * std::cos(2 * collision.psiFT0C()), collision.qFT0C() * std::sin(2 * collision.psiFT0C()), collision.qFT0C()); - StraFT0CQVsEv(collision.triggereventep()); + products.StraFT0CQVs(collision.qFT0C() * std::cos(2 * collision.psiFT0C()), collision.qFT0C() * std::sin(2 * collision.psiFT0C()), collision.qFT0C()); + products.StraFT0CQVsEv(collision.triggereventep()); } void processZDCSP(soa::Join::iterator const& collision) { - StraZDCSP(collision.triggereventsp(), - collision.psiZDCA(), collision.psiZDCC(), collision.qxZDCA(), collision.qxZDCC(), collision.qyZDCA(), collision.qyZDCC()); + products.StraZDCSP(collision.triggereventsp(), + collision.psiZDCA(), collision.psiZDCC(), collision.qxZDCA(), collision.qxZDCC(), collision.qyZDCA(), collision.qyZDCC()); } void processFT0MQVectors(soa::Join::iterator const& collision) { - StraFT0MQVs(collision.qvecFT0MRe(), collision.qvecFT0MIm(), collision.sumAmplFT0M()); + products.StraFT0MQVs(collision.qvecFT0MRe(), collision.qvecFT0MIm(), collision.sumAmplFT0M()); } void processFV0AQVectors(soa::Join::iterator const& collision) { - StraFV0AQVs(collision.qvecFV0ARe(), collision.qvecFV0AIm(), collision.sumAmplFV0A()); + products.StraFV0AQVs(collision.qvecFV0ARe(), collision.qvecFV0AIm(), collision.sumAmplFV0A()); } void processTPCQVectors(soa::Join::iterator const& collision) { - StraTPCQVs(collision.qvecBNegRe(), collision.qvecBNegIm(), collision.nTrkBNeg(), collision.qvecBPosRe(), collision.qvecBPosIm(), collision.nTrkBPos()); + products.StraTPCQVs(collision.qvecBNegRe(), collision.qvecBNegIm(), collision.nTrkBNeg(), collision.qvecBPosRe(), collision.qvecBPosIm(), collision.nTrkBPos()); } void processTPCQVectorsLF(soa::Join::iterator const& collision) { - StraTPCQVs(collision.qTPCL() * std::cos(2 * collision.psiTPCL()), collision.qTPCL() * std::sin(2 * collision.psiTPCL()), collision.qTPCL(), collision.qTPCR() * std::cos(2 * collision.psiTPCR()), collision.qTPCR() * std::sin(2 * collision.psiTPCR()), collision.qTPCR()); + products.StraTPCQVs(collision.qTPCL() * std::cos(2 * collision.psiTPCL()), collision.qTPCL() * std::sin(2 * collision.psiTPCL()), collision.qTPCL(), collision.qTPCR() * std::cos(2 * collision.psiTPCR()), collision.qTPCR() * std::sin(2 * collision.psiTPCR()), collision.qTPCR()); } using uint128_t = __uint128_t; @@ -1077,7 +1101,7 @@ struct strangederivedbuilder { void processDataframeIDs(aod::Origins const& origins) { auto origin = origins.begin(); - straOrigin(origin.dataframeID()); + products.straOrigin(origin.dataframeID()); } void processSimulatedZDCNeutrons(soa::Join const& mcCollisions, aod::McParticles const& mcParticlesEntireTable) @@ -1089,11 +1113,11 @@ struct strangederivedbuilder { for (const auto& mcPart : mcParticles) { if (std::abs(mcPart.pdgCode()) == kNeutron) { // check if it is a neutron or anti-neutron if (std::abs(mcPart.eta()) > 8.7) { // check if it is within ZDC acceptance - zdcNeutrons(mcPart.pdgCode(), mcPart.statusCode(), mcPart.flags(), - mcPart.vx(), mcPart.vy(), mcPart.vz(), mcPart.vt(), - mcPart.px(), mcPart.py(), mcPart.pz(), mcPart.e()); + products.zdcNeutrons(mcPart.pdgCode(), mcPart.statusCode(), mcPart.flags(), + mcPart.vx(), mcPart.vy(), mcPart.vz(), mcPart.vt(), + mcPart.px(), mcPart.py(), mcPart.pz(), mcPart.e()); - zdcNeutronsMCCollRefs(mcPart.mcCollisionId()); + products.zdcNeutronsMCCollRefs(mcPart.mcCollisionId()); } } } @@ -1111,6 +1135,8 @@ struct strangederivedbuilder { // Run 2: collision processing PROCESS_SWITCH(strangederivedbuilder, processCollisionsRun2, "Produce collisions (Run2)", false); PROCESS_SWITCH(strangederivedbuilder, processCollisionsRun2WithMC, "Produce collisions (Run 2) with MC info", false); + // Event times + PROCESS_SWITCH(strangederivedbuilder, processPopulateEventTimes, "Populate event times", true); // detailed information processing PROCESS_SWITCH(strangederivedbuilder, processTrackExtrasV0sOnly, "Produce track extra information (V0s only)", true); diff --git a/PWGLF/TableProducer/Strangeness/strangenesstofpid.cxx b/PWGLF/TableProducer/Strangeness/strangenesstofpid.cxx index f74a6d95d1d..b19d7ba04c8 100644 --- a/PWGLF/TableProducer/Strangeness/strangenesstofpid.cxx +++ b/PWGLF/TableProducer/Strangeness/strangenesstofpid.cxx @@ -70,7 +70,7 @@ using CascOriginalDatas = soa::Join; using TracksWithAllExtras = soa::Join; // For derived data analysis -using dauTracks = soa::Join; +using dauTracks = soa::Join; using V0DerivedDatas = soa::Join; using V0DerivedDatasMC = soa::Join; using CascDerivedDatas = soa::Join; @@ -104,22 +104,6 @@ struct strangenesstofpid { struct : ConfigurableGroup { Configurable d_bz_input{"d_bz", -999, "bz field, -999 is automatic"}; Configurable tofPosition{"tofPosition", 377.934f, "TOF effective (inscribed) radius"}; - - Configurable correctELossInclination{"correctELossInclination", false, "factor out inclination when doing effective e-loss correction (0: no, 1: yes)"}; - Configurable numberOfStepsFirstStage{"numberOfStepsFirstStage", 500, "Max number of alpha rotations to attempt in first stage"}; - Configurable numberOfStepsSecondStage{"numberOfStepsSecondStage", 500, "Max number of steps rotations to attempt in second stage"}; - Configurable stepSizeFirstStage{"stepSizeFirstStage", 2.0f, "Max number of alpha rotations to attempt in first stage"}; - Configurable firstApproximationThreshold{"firstApproximationThreshold", 4.0f, "be satisfied if first approach to TOF radius is OK within this threshold (cm)"}; - - // regulate e-loss calculation to save CPU - Configurable maxPionMomentumForEloss{"maxPionMomentumForEloss", 1.5f, "above this momentum, do fast analytical TOF calculation for pions"}; - Configurable maxKaonMomentumForEloss{"maxKaonMomentumForEloss", 2.0f, "above this momentum, do fast analytical TOF calculation for kaons"}; - Configurable maxProtonMomentumForEloss{"maxProtonMomentumForEloss", 2.5f, "above this momentum, do fast analytical TOF calculation for protons"}; - - Configurable rejectKaonMomentumForEloss{"rejectKaonMomentumForEloss", 0.2f, "below this momentum, reject kaon hypothesis (won't reach TOF)"}; - Configurable rejectProtonMomentumForEloss{"rejectProtonMomentumForEloss", 0.2f, "below this momentum, reject proton hypothesis (won't reach TOF)"}; - - Configurable tpcNsigmaThreshold{"tpcNsigmaThreshold", 6.0f, "require TPC compatibility to attempt eloss propagation (otherwise, don't calculate)"}; } propagationConfiguration; Configurable doQA{"doQA", false, "create QA histos"}; @@ -335,138 +319,6 @@ struct strangenesstofpid { return length; } - /// O2 Propagator + TrackLTIntegral approach helpers - - /// function to calculate segmented (truncated) radius based on a certain x, y position - float segmentedRadius(float x, float y) - { - float atAngle = std::atan2(y, x); - float roundedAngle = TMath::Pi() / 9; // 18 segments = use 9 here - float angleSegmentAxis = 0.5f * roundedAngle + roundedAngle * static_cast(std::floor(atAngle / roundedAngle)); - float xSegmentAxis = TMath::Cos(angleSegmentAxis); - float ySegmentAxis = TMath::Sin(angleSegmentAxis); - return xSegmentAxis * x + ySegmentAxis * y; // inner product - } - - /// function to calculate track length of this track up to a certain segmented detector - /// \param track the input track - /// \param time returned time (with PID given by track PID) - void calculateTOF(o2::track::TrackPar track, float& time) - { - time = -1e+6; - - if (track.getPID() == o2::track::PID::Proton && track.getP() < propagationConfiguration.rejectProtonMomentumForEloss.value) { - return; // don't attempt to calculate below-threshold protons (will stop anyway) - } - if (track.getPID() == o2::track::PID::Kaon && track.getP() < propagationConfiguration.rejectKaonMomentumForEloss.value) { - return; // don't attempt to calculate below-threshold kaons (will stop anyway) - } - - o2::track::TrackLTIntegral ltIntegral; - - static constexpr float MAX_SIN_PHI = 0.85f; - static constexpr float MAX_STEP = 2.0f; - static constexpr float MAX_STEP_FINAL_STAGE = 0.5f; - static constexpr float MAX_FINAL_X = 450.0f; // maximum extra X on top of TOF X for correcting value - - //____________________________________________________________ - // stage 1: propagate to TOF-inscribed circle at tofPosition - - // define standard variables - std::array xyz; - std::array pxpypz; - track.getXYZGlo(xyz); - float segmentedRstart = segmentedRadius(xyz[0], xyz[1]); - - bool firstPropag = true; - for (int iRot = 0; iRot < propagationConfiguration.numberOfStepsFirstStage.value; iRot++) { - track.rotateParam(track.getPhi()); // start rotated - float currentX = track.getX(); - float tofX = currentX; - track.getXatLabR(propagationConfiguration.tofPosition, tofX, d_bz, o2::track::DirOutward); - if (std::abs(tofX - currentX) < propagationConfiguration.firstApproximationThreshold.value) { - // signal conclusion - if (calculationMethod.value == 2) { - histos.fill(HIST("hInitialPropagationSteps"), iRot); // store number of steps - } - break; - } - float nextX = std::min(currentX + propagationConfiguration.stepSizeFirstStage.value, tofX); - firstPropag = o2::base::Propagator::Instance()->propagateToX(track, nextX, d_bz, MAX_SIN_PHI, MAX_STEP, o2::base::Propagator::MatCorrType::USEMatCorrLUT, <Integral); - } - - // mark start position of next step - track.getXYZGlo(xyz); - float snp = track.getSnp(); - float segmentedR = segmentedRadius(xyz[0], xyz[1]); - float segmentedRintermediate = segmentedR; - float currentTime = ltIntegral.getTOF(track.getPID()); - if (calculationMethod.value == 2) { - histos.fill(HIST("hSegRadiusFirstPropagVsStart"), segmentedRstart, segmentedR); // for debugging purposes - histos.fill(HIST("hSnp"), track.getSnp()); // for debugging purposes - histos.fill(HIST("hTOFPosition"), xyz[0], xyz[1]); // for debugging purposes - histos.fill(HIST("hSegRadius"), segmentedR); // for debugging purposes - if (!firstPropag) { - histos.fill(HIST("hTOFPositionFirstPropagFail"), xyz[0], xyz[1]); // for debugging purposes - histos.fill(HIST("hSegRadiusFirstPropagFail"), segmentedR); // for debugging purposes - histos.fill(HIST("hSnpFirstPropagFail"), track.getSnp()); // for debugging purposes - } - } - - // correct for TOF segmentation - bool trackOKextra = true; - float trackXextra = track.getX(); - int propagationSteps = 0; - int maxPropagationSteps = propagationConfiguration.numberOfStepsSecondStage.value; - while ((trackXextra < MAX_FINAL_X) && (propagationSteps < maxPropagationSteps)) { - // continue with alpha aligned with pT - track.getPxPyPzGlo(pxpypz); - track.rotateParam(std::atan2(pxpypz[1], pxpypz[0])); - trackXextra = track.getX() + MAX_STEP_FINAL_STAGE; - trackOKextra = o2::base::Propagator::Instance()->propagateToX(track, trackXextra, d_bz, MAX_SIN_PHI, MAX_STEP, o2::base::Propagator::MatCorrType::USEMatCorrLUT, <Integral); - if (!trackOKextra) { - if (calculationMethod.value == 2) { - track.getXYZGlo(xyz); - histos.fill(HIST("hTOFPositionStopped"), xyz[0], xyz[1]); // for debugging purposes - histos.fill(HIST("hSnpStopped"), snp); // for debugging purposes - histos.fill(HIST("hSegRadiusStopped"), segmentedRadius(xyz[0], xyz[1])); // for debugging purposes - histos.fill(HIST("hSegRadiusStoppedVsFirstPropag"), segmentedRintermediate, segmentedRadius(xyz[0], xyz[1])); // for debugging purposes - } - time = -1e+6; - return; // propagation failed, skip, won't look reasonable - } - - // re-evaluate - did we cross? if yes break - float previousX = xyz[0], previousY = xyz[1]; - track.getXYZGlo(xyz); - if (segmentedRadius(xyz[0], xyz[1]) > propagationConfiguration.tofPosition) { - // crossed boundary -> do proportional scaling with how much we actually crossed the boundary - float segmentedRFinal = segmentedRadius(xyz[0], xyz[1]); - float timeFinal = ltIntegral.getTOF(track.getPID()); - float fraction = (propagationConfiguration.tofPosition - segmentedR) / (segmentedRFinal - segmentedR + 1e-6); // proportional fraction - time = currentTime + (timeFinal - currentTime) * fraction; - if (calculationMethod.value == 2) { - histos.fill(HIST("hTOFPositionFinal"), previousX + fraction * (xyz[0] - previousX), previousY + fraction * (xyz[1] - previousY)); // for debugging purposes - histos.fill(HIST("hSegRadiusFinal"), segmentedRadius(previousX + fraction * (xyz[0] - previousX), previousY + fraction * (xyz[1] - previousY))); // for debugging purposes - histos.fill(HIST("hSnpFinal"), track.getSnp()); // for debugging purposes - histos.fill(HIST("hSegRadiusFinalVsFirstPropag"), segmentedRintermediate, segmentedRadius(previousX + fraction * (xyz[0] - previousX), previousY + fraction * (xyz[1] - previousY))); // for debugging purposes - histos.fill(HIST("hRefinedPropagationSteps"), propagationSteps, 1.0f); // for debugging purposes - } - return; // get out of the entire function and return (don't just break) - } - - // prepare for next step by setting current position and desired variables - segmentedR = segmentedRadius(xyz[0], xyz[1]); - currentTime = ltIntegral.getTOF(track.getPID()); - propagationSteps++; - } - if (calculationMethod.value == 2) { - histos.fill(HIST("hRefinedPropagationSteps"), propagationSteps, 0.0f); // for debugging purposes - track.getXYZGlo(xyz); - histos.fill(HIST("hSegRadiusGotLost"), segmentedRadius(xyz[0], xyz[1])); // for debugging purposes - } - } - void init(InitContext& initContext) { if (calculateV0s.value < 0) { @@ -566,93 +418,6 @@ struct strangenesstofpid { // delta lambda decay time histos.add("h2dLambdaDeltaDecayTime", "h2dLambdaDeltaDecayTime", {HistType::kTH2F, {axisP, axisDeltaTime}}); } - - if (calculationMethod.value == 2) { - //_____________________________________________________________________ - // special mode in which comparison histograms are required - - //_____________________________________________________________________ - // histograms for debugging modes 0 vs 1 - // encoded success rates in each hypothesis and method vs prong p - histos.add("h2dSucessRatePion", "h2dSucessRatePion", kTH2F, {axisSmallP, {4, -0.5f, 3.5f}}); - histos.add("h2dSucessRateKaon", "h2dSucessRateKaon", kTH2F, {axisSmallP, {4, -0.5f, 3.5f}}); - histos.add("h2dSucessRateProton", "h2dSucessRateProton", kTH2F, {axisSmallP, {4, -0.5f, 3.5f}}); - - histos.add("hInitialPropagationSteps", "hInitialPropagationSteps", kTH1F, {{500, -0.5f, 499.5f}}); - histos.add("hRefinedPropagationSteps", "hRefinedPropagationSteps", kTH2F, {{1000, -0.5f, 999.5f}, {2, -0.5f, 1.5f}}); - - // base ArcDebug: comparison between times of arrival in different methods - histos.add("hArcDebug", "hArcDebug", kTH2F, {axisTime, axisTime}); - - // Position of TrackLTIntegral method: intermediate (getXatLabR) and final (reach segmented detector) - histos.add("hTOFPosition", "hTOFPosition", kTH2F, {axisPosition, axisPosition}); - histos.add("hTOFPositionFirstPropagFail", "hTOFPositionFirstPropagFail", kTH2F, {axisPosition, axisPosition}); - histos.add("hTOFPositionFinal", "hTOFPositionFinal", kTH2F, {axisPosition, axisPosition}); - histos.add("hTOFPositionGetXAtLabFail", "hTOFPositionGetXAtLabFail", kTH2F, {axisPosition, axisPosition}); - histos.add("hTOFPositionStopped", "hTOFPositionStopped", kTH2F, {axisPosition, axisPosition}); - - // Snp cross-check - histos.add("hSnp", "hSnp", kTH1F, {axisSnp}); - histos.add("hSnpFirstPropagFail", "hSnpFirstPropagFail", kTH1F, {axisSnp}); - histos.add("hSnpFinal", "hSnpFinal", kTH1F, {axisSnp}); - histos.add("hSnpGetXAtLabFail", "hSnpGetXAtLabFail", kTH1F, {axisSnp}); - histos.add("hSnpStopped", "hSnpStopped", kTH1F, {axisSnp}); - - // segmented radius: positions - histos.add("hSegRadius", "hSegRadius", kTH1F, {{400, 0.0f, 400.0f}}); - histos.add("hSegRadiusFirstPropagFail", "hSegRadiusFirstPropagFail", kTH1F, {{400, 0.0f, 400.0f}}); - histos.add("hSegRadiusFinal", "hSegRadiusFinal", kTH1F, {{400, 0.0f, 400.0f}}); - histos.add("hSegRadiusStopped", "hSegRadiusStopped", kTH1F, {{400, 0.0f, 400.0f}}); - histos.add("hSegRadiusGotLost", "hSegRadiusGotLost", kTH1F, {{400, 0.0f, 400.0f}}); - - histos.add("hSegRadiusFirstPropagVsStart", "hSegRadiusFirstPropagVsStart", kTH2F, {{400, 0.0f, 400.0f}, {400, 0.0f, 400.0f}}); - histos.add("hSegRadiusStoppedVsFirstPropag", "hSegRadiusStoppedVsFirstPropag", kTH2F, {{400, 0.0f, 400.0f}, {400, 0.0f, 400.0f}}); - histos.add("hSegRadiusFinalVsFirstPropag", "hSegRadiusFinalVsFirstPropag", kTH2F, {{400, 0.0f, 400.0f}, {400, 0.0f, 400.0f}}); - - // Delta-times of each method for the various species - histos.add("hDeltaTimeMethodsVsP_posLaPr", "hDeltaTimeMethodsVsP_posLaPr", kTH3F, {axisSmallP, axisEta, axisDeltaTime}); - histos.add("hDeltaTimeMethodsVsP_posLaPi", "hDeltaTimeMethodsVsP_posLaPi", kTH3F, {axisSmallP, axisEta, axisDeltaTime}); - histos.add("hDeltaTimeMethodsVsP_posK0Pi", "hDeltaTimeMethodsVsP_posK0Pi", kTH3F, {axisSmallP, axisEta, axisDeltaTime}); - histos.add("hDeltaTimeMethodsVsP_negLaPr", "hDeltaTimeMethodsVsP_negLaPr", kTH3F, {axisSmallP, axisEta, axisDeltaTime}); - histos.add("hDeltaTimeMethodsVsP_negLaPi", "hDeltaTimeMethodsVsP_negLaPi", kTH3F, {axisSmallP, axisEta, axisDeltaTime}); - histos.add("hDeltaTimeMethodsVsP_negK0Pi", "hDeltaTimeMethodsVsP_negK0Pi", kTH3F, {axisSmallP, axisEta, axisDeltaTime}); - - histos.add("hDeltaTimeMethodsVsP_posXiPi", "hDeltaTimeMethodsVsP_posXiPi", kTH3F, {axisSmallP, axisEta, axisDeltaTime}); - histos.add("hDeltaTimeMethodsVsP_posXiPr", "hDeltaTimeMethodsVsP_posXiPr", kTH3F, {axisSmallP, axisEta, axisDeltaTime}); - histos.add("hDeltaTimeMethodsVsP_negXiPi", "hDeltaTimeMethodsVsP_negXiPi", kTH3F, {axisSmallP, axisEta, axisDeltaTime}); - histos.add("hDeltaTimeMethodsVsP_negXiPr", "hDeltaTimeMethodsVsP_negXiPr", kTH3F, {axisSmallP, axisEta, axisDeltaTime}); - histos.add("hDeltaTimeMethodsVsP_bachXiPi", "hDeltaTimeMethodsVsP_bachXiPi", kTH3F, {axisSmallP, axisEta, axisDeltaTime}); - - histos.add("hDeltaTimeMethodsVsP_posOmPi", "hDeltaTimeMethodsVsP_posOmPi", kTH3F, {axisSmallP, axisEta, axisDeltaTime}); - histos.add("hDeltaTimeMethodsVsP_posOmPr", "hDeltaTimeMethodsVsP_posOmPr", kTH3F, {axisSmallP, axisEta, axisDeltaTime}); - histos.add("hDeltaTimeMethodsVsP_negOmPi", "hDeltaTimeMethodsVsP_negOmPi", kTH3F, {axisSmallP, axisEta, axisDeltaTime}); - histos.add("hDeltaTimeMethodsVsP_negOmPr", "hDeltaTimeMethodsVsP_negOmPr", kTH3F, {axisSmallP, axisEta, axisDeltaTime}); - histos.add("hDeltaTimeMethodsVsP_bachOmKa", "hDeltaTimeMethodsVsP_bachOmKa", kTH3F, {axisSmallP, axisEta, axisDeltaTime}); - - histos.add("hRatioTimeMethodsVsP_posLaPr", "hRatioTimeMethodsVsP_posLaPr", kTH3F, {axisSmallP, axisEta, axisRatioMethods}); - histos.add("hRatioTimeMethodsVsP_posLaPi", "hRatioTimeMethodsVsP_posLaPi", kTH3F, {axisSmallP, axisEta, axisRatioMethods}); - histos.add("hRatioTimeMethodsVsP_posK0Pi", "hRatioTimeMethodsVsP_posK0Pi", kTH3F, {axisSmallP, axisEta, axisRatioMethods}); - histos.add("hRatioTimeMethodsVsP_negLaPr", "hRatioTimeMethodsVsP_negLaPr", kTH3F, {axisSmallP, axisEta, axisRatioMethods}); - histos.add("hRatioTimeMethodsVsP_negLaPi", "hRatioTimeMethodsVsP_negLaPi", kTH3F, {axisSmallP, axisEta, axisRatioMethods}); - histos.add("hRatioTimeMethodsVsP_negK0Pi", "hRatioTimeMethodsVsP_negK0Pi", kTH3F, {axisSmallP, axisEta, axisRatioMethods}); - - histos.add("hRatioTimeMethodsVsP_posXiPi", "hRatioTimeMethodsVsP_posXiPi", kTH3F, {axisSmallP, axisEta, axisRatioMethods}); - histos.add("hRatioTimeMethodsVsP_posXiPr", "hRatioTimeMethodsVsP_posXiPr", kTH3F, {axisSmallP, axisEta, axisRatioMethods}); - histos.add("hRatioTimeMethodsVsP_negXiPi", "hRatioTimeMethodsVsP_negXiPi", kTH3F, {axisSmallP, axisEta, axisRatioMethods}); - histos.add("hRatioTimeMethodsVsP_negXiPr", "hRatioTimeMethodsVsP_negXiPr", kTH3F, {axisSmallP, axisEta, axisRatioMethods}); - histos.add("hRatioTimeMethodsVsP_bachXiPi", "hRatioTimeMethodsVsP_bachXiPi", kTH3F, {axisSmallP, axisEta, axisRatioMethods}); - - histos.add("hRatioTimeMethodsVsP_posOmPi", "hRatioTimeMethodsVsP_posOmPi", kTH3F, {axisSmallP, axisEta, axisRatioMethods}); - histos.add("hRatioTimeMethodsVsP_posOmPr", "hRatioTimeMethodsVsP_posOmPr", kTH3F, {axisSmallP, axisEta, axisRatioMethods}); - histos.add("hRatioTimeMethodsVsP_negOmPi", "hRatioTimeMethodsVsP_negOmPi", kTH3F, {axisSmallP, axisEta, axisRatioMethods}); - histos.add("hRatioTimeMethodsVsP_negOmPr", "hRatioTimeMethodsVsP_negOmPr", kTH3F, {axisSmallP, axisEta, axisRatioMethods}); - histos.add("hRatioTimeMethodsVsP_bachOmKa", "hRatioTimeMethodsVsP_bachOmKa", kTH3F, {axisSmallP, axisEta, axisRatioMethods}); - } - - // list memory consumption at start if running in modes with more output - if (calculationMethod.value == 2 || doQA) { - histos.print(); - } } void initCCDB(int runNumber) @@ -767,15 +532,15 @@ struct strangenesstofpid { } } - if (calculationMethod.value > 0 && !lut) { - // setMatLUT only after magfield has been initalized - // (setMatLUT has implicit and problematic init field call if not) - LOG(info) << "Loading full (all-radius) material look-up table for run number: " << runNumber; - lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->getForRun(ccdbConfigurations.lutPath, runNumber)); - o2::base::Propagator::Instance()->setMatLUT(lut); - o2::base::Propagator::Instance()->setTGeoFallBackAllowed(false); - LOG(info) << "Material look-up table loaded!"; - } + // if (calculationMethod.value > 0 && !lut) { + // // setMatLUT only after magfield has been initalized + // // (setMatLUT has implicit and problematic init field call if not) + // LOG(info) << "Loading full (all-radius) material look-up table for run number: " << runNumber; + // lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->getForRun(ccdbConfigurations.lutPath, runNumber)); + // o2::base::Propagator::Instance()->setMatLUT(lut); + // o2::base::Propagator::Instance()->setTGeoFallBackAllowed(false); + // LOG(info) << "Material look-up table loaded!"; + // } mRunNumber = runNumber; } @@ -787,24 +552,12 @@ struct strangenesstofpid { return 0.0299792458 * TMath::Sqrt(lA / (1 + lA)); } - // templatized process function for symmetric operation in derived and original AO2D - template - void processV0Candidate(TCollision const& collision, TV0 const& v0, TTrack const& pTra, TTrack const& nTra, int v0pdg) - { - // time of V0 segment - float lengthV0 = std::hypot(v0.x() - collision.getX(), v0.y() - collision.getY(), v0.z() - collision.getZ()); - float velocityK0Short = velocity(v0.p(), o2::constants::physics::MassKaonNeutral); - float velocityLambda = velocity(v0.p(), o2::constants::physics::MassLambda); - float timeK0Short = lengthV0 / velocityK0Short; // in picoseconds - float timeLambda = lengthV0 / velocityLambda; // in picoseconds - - // initialize from V0 position and momenta - o2::track::TrackPar posTrack = o2::track::TrackPar({v0.x(), v0.y(), v0.z()}, {v0.pxpos(), v0.pypos(), v0.pzpos()}, +1, false); - o2::track::TrackPar negTrack = o2::track::TrackPar({v0.x(), v0.y(), v0.z()}, {v0.pxneg(), v0.pyneg(), v0.pzneg()}, -1, false); - - // at minimum - float positiveP = std::hypot(v0.pxpos(), v0.pypos(), v0.pzpos()); - float negativeP = std::hypot(v0.pxneg(), v0.pyneg(), v0.pzneg()); + // structs to hold information + struct v0TofInfo { // holds processed information regarding TOF for V0s + float timePositivePr = o2::aod::v0data::kNoTOFValue; + float timePositivePi = o2::aod::v0data::kNoTOFValue; + float timeNegativePr = o2::aod::v0data::kNoTOFValue; + float timeNegativePi = o2::aod::v0data::kNoTOFValue; float deltaTimePositiveLambdaPi = o2::aod::v0data::kNoTOFValue; float deltaTimeNegativeLambdaPi = o2::aod::v0data::kNoTOFValue; @@ -819,304 +572,229 @@ struct strangenesstofpid { float nSigmaNegativeLambdaPr = o2::aod::v0data::kNoTOFValue; float nSigmaPositiveK0ShortPi = o2::aod::v0data::kNoTOFValue; float nSigmaNegativeK0ShortPi = o2::aod::v0data::kNoTOFValue; + }; - float timePositivePr = o2::aod::v0data::kNoTOFValue; - float timePositivePi = o2::aod::v0data::kNoTOFValue; - float timeNegativePr = o2::aod::v0data::kNoTOFValue; - float timeNegativePi = o2::aod::v0data::kNoTOFValue; - - float timePositivePr_Method0 = o2::aod::v0data::kNoTOFValue; - float timePositivePi_Method0 = o2::aod::v0data::kNoTOFValue; - float timeNegativePr_Method0 = o2::aod::v0data::kNoTOFValue; - float timeNegativePi_Method0 = o2::aod::v0data::kNoTOFValue; - - float timePositivePr_Method1 = o2::aod::v0data::kNoTOFValue; - float timePositivePi_Method1 = o2::aod::v0data::kNoTOFValue; - float timeNegativePr_Method1 = o2::aod::v0data::kNoTOFValue; - float timeNegativePi_Method1 = o2::aod::v0data::kNoTOFValue; - - float velocityPositivePr = velocity(posTrack.getP(), o2::constants::physics::MassProton); - float velocityPositivePi = velocity(posTrack.getP(), o2::constants::physics::MassPionCharged); - float velocityNegativePr = velocity(negTrack.getP(), o2::constants::physics::MassProton); - float velocityNegativePi = velocity(negTrack.getP(), o2::constants::physics::MassPionCharged); - - if (pTra.hasTOF() && pTra.hasITS()) { - float lengthPositive = findInterceptLength(posTrack, d_bz); // FIXME: tofPosition ok? adjust? - if (lengthPositive > 0) { - timePositivePr_Method0 = lengthPositive / velocityPositivePr; - timePositivePi_Method0 = lengthPositive / velocityPositivePi; - } - } - if (nTra.hasTOF() && nTra.hasITS()) { - float lengthNegative = findInterceptLength(negTrack, d_bz); // FIXME: tofPosition ok? adjust? - if (lengthNegative > 0) { - timeNegativePr_Method0 = lengthNegative / velocityNegativePr; - timeNegativePi_Method0 = lengthNegative / velocityNegativePi; - } - } - - if (calculationMethod.value > 0) { - // method to calculate the time and length via Propagator TrackLTIntegral - if (pTra.hasTOF() && pTra.hasITS()) { // calculate if signal present, otherwise skip - if (posTrack.getP() < propagationConfiguration.maxProtonMomentumForEloss.value && std::abs(pTra.tpcNSigmaPr()) < propagationConfiguration.tpcNsigmaThreshold) { - o2::track::TrackPar posTrackAsProton(posTrack); - posTrackAsProton.setPID(o2::track::PID::Proton); - calculateTOF(posTrackAsProton, timePositivePr_Method1); - } else { - timePositivePr_Method1 = timePositivePr_Method0; - } - if (posTrack.getP() < propagationConfiguration.maxPionMomentumForEloss.value && std::abs(pTra.tpcNSigmaPi()) < propagationConfiguration.tpcNsigmaThreshold) { - o2::track::TrackPar posTrackAsPion(posTrack); - posTrackAsPion.setPID(o2::track::PID::Pion); - calculateTOF(posTrackAsPion, timePositivePi_Method1); - } else { - timePositivePi_Method1 = timePositivePi_Method0; - } - } - if (nTra.hasTOF() && nTra.hasITS()) { // calculate if signal present, otherwise skip - if (negTrack.getP() < propagationConfiguration.maxProtonMomentumForEloss.value && std::abs(nTra.tpcNSigmaPr()) < propagationConfiguration.tpcNsigmaThreshold) { - o2::track::TrackPar negTrackAsProton(negTrack); - negTrackAsProton.setPID(o2::track::PID::Proton); - calculateTOF(negTrackAsProton, timeNegativePr_Method1); - } else { - timeNegativePr_Method1 = timeNegativePr_Method0; - } - - if (negTrack.getP() < propagationConfiguration.maxPionMomentumForEloss.value && std::abs(nTra.tpcNSigmaPi()) < propagationConfiguration.tpcNsigmaThreshold) { - o2::track::TrackPar negTrackAsPion(negTrack); - negTrackAsPion.setPID(o2::track::PID::Pion); - calculateTOF(negTrackAsPion, timeNegativePi_Method1); - } else { - timeNegativePi_Method1 = timeNegativePi_Method0; - } - } - } - - // assign values to be used in main calculation - if (calculationMethod.value == 0) { - timePositivePr = timePositivePr_Method0; - timePositivePi = timePositivePi_Method0; - timeNegativePr = timeNegativePr_Method0; - timeNegativePi = timeNegativePi_Method0; - } else { - timePositivePr = timePositivePr_Method1; - timePositivePi = timePositivePi_Method1; - timeNegativePr = timeNegativePr_Method1; - timeNegativePi = timeNegativePi_Method1; - } - - if (calculationMethod.value == 2) { - // do analysis of successes and failures - bool positiveSuccessMethod0Pr = std::abs(timePositivePr_Method0 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon; - bool negativeSuccessMethod0Pr = std::abs(timeNegativePr_Method0 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon; - bool positiveSuccessMethod1Pr = std::abs(timePositivePr_Method1 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon; - bool negativeSuccessMethod1Pr = std::abs(timeNegativePr_Method1 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon; - bool positiveSuccessMethod0Pi = std::abs(timePositivePi_Method0 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon; - bool negativeSuccessMethod0Pi = std::abs(timeNegativePi_Method0 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon; - bool positiveSuccessMethod1Pi = std::abs(timePositivePi_Method1 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon; - bool negativeSuccessMethod1Pi = std::abs(timeNegativePi_Method1 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon; - - int encodedPositiveSuccessPi = positiveSuccessMethod0Pi + 2 * positiveSuccessMethod1Pi; - int encodedPositiveSuccessPr = positiveSuccessMethod0Pr + 2 * positiveSuccessMethod1Pr; - int encodedNegativeSuccessPi = negativeSuccessMethod0Pi + 2 * negativeSuccessMethod1Pi; - int encodedNegativeSuccessPr = negativeSuccessMethod0Pr + 2 * negativeSuccessMethod1Pr; - - if (pTra.hasTOF()) { - histos.fill(HIST("h2dSucessRateProton"), positiveP, encodedPositiveSuccessPr); - histos.fill(HIST("h2dSucessRatePion"), positiveP, encodedPositiveSuccessPi); - } - if (nTra.hasTOF()) { - histos.fill(HIST("h2dSucessRateProton"), negativeP, encodedNegativeSuccessPr); - histos.fill(HIST("h2dSucessRatePion"), negativeP, encodedNegativeSuccessPi); - } - } - - if (pTra.hasTOF() && pTra.hasITS() && timePositivePr > 0) { - deltaTimePositiveLambdaPr = (pTra.tofSignal() - pTra.tofEvTime()) - (timeLambda + timePositivePr); - deltaTimePositiveLambdaPi = (pTra.tofSignal() - pTra.tofEvTime()) - (timeLambda + timePositivePi); - deltaTimePositiveK0ShortPi = (pTra.tofSignal() - pTra.tofEvTime()) - (timeK0Short + timePositivePi); - } - if (nTra.hasTOF() && nTra.hasITS() && timeNegativePr > 0) { - deltaTimeNegativeLambdaPr = (nTra.tofSignal() - nTra.tofEvTime()) - (timeLambda + timeNegativePr); - deltaTimeNegativeLambdaPi = (nTra.tofSignal() - nTra.tofEvTime()) - (timeLambda + timeNegativePi); - deltaTimeNegativeK0ShortPi = (nTra.tofSignal() - nTra.tofEvTime()) - (timeK0Short + timeNegativePi); - } - - if (doQA) { - // calculate and pack properties for QA purposes - int posProperties = 0; - if (timePositivePr > 0) - posProperties = posProperties | (static_cast(1) << kLength); - if (pTra.hasTOF()) - posProperties = posProperties | (static_cast(1) << kHasTOF); - int negProperties = 0; - if (timeNegativePr > 0) - negProperties = negProperties | (static_cast(1) << kLength); - if (nTra.hasTOF()) - negProperties = negProperties | (static_cast(1) << kHasTOF); - - histos.fill(HIST("h2dPositiveTOFProperties"), v0.p(), posProperties); - histos.fill(HIST("h2dNegativeTOFProperties"), v0.p(), negProperties); - } - - float deltaDecayTimeLambda = -10e+4; - float deltaDecayTimeAntiLambda = -10e+4; - float deltaDecayTimeK0Short = -10e+4; - if (nTra.hasTOF() && pTra.hasTOF() && timePositivePr > 0 && timeNegativePr > 0) { // does not depend on event time - deltaDecayTimeLambda = (pTra.tofSignal() - timePositivePr) - (nTra.tofSignal() - timeNegativePi); - deltaDecayTimeAntiLambda = (pTra.tofSignal() - timePositivePi) - (nTra.tofSignal() - timeNegativePr); - deltaDecayTimeK0Short = (pTra.tofSignal() - timePositivePi) - (nTra.tofSignal() - timeNegativePi); - } + // structs to hold information + struct cascTofInfo { // holds processed information regarding TOF for Cascades + float posFlightPi = o2::aod::cascdata::kNoTOFValue; + float posFlightPr = o2::aod::cascdata::kNoTOFValue; + float negFlightPi = o2::aod::cascdata::kNoTOFValue; + float negFlightPr = o2::aod::cascdata::kNoTOFValue; + float bachFlightPi = o2::aod::cascdata::kNoTOFValue; + float bachFlightKa = o2::aod::cascdata::kNoTOFValue; - // calculate betas + float posDeltaTimeAsXiPi = o2::aod::cascdata::kNoTOFValue, posDeltaTimeAsXiPr = o2::aod::cascdata::kNoTOFValue; + float negDeltaTimeAsXiPi = o2::aod::cascdata::kNoTOFValue, negDeltaTimeAsXiPr = o2::aod::cascdata::kNoTOFValue; + float bachDeltaTimeAsXiPi = o2::aod::cascdata::kNoTOFValue; + float posDeltaTimeAsOmPi = o2::aod::cascdata::kNoTOFValue, posDeltaTimeAsOmPr = o2::aod::cascdata::kNoTOFValue; + float negDeltaTimeAsOmPi = o2::aod::cascdata::kNoTOFValue, negDeltaTimeAsOmPr = o2::aod::cascdata::kNoTOFValue; + float bachDeltaTimeAsOmKa = o2::aod::cascdata::kNoTOFValue; - float evTimeMean = 0.5f * (pTra.tofEvTime() + nTra.tofEvTime()); - float decayTimeLambda = 0.5f * ((pTra.tofSignal() - timePositivePr) + (nTra.tofSignal() - timeNegativePi)) - evTimeMean; - float decayTimeAntiLambda = 0.5f * ((pTra.tofSignal() - timePositivePi) + (nTra.tofSignal() - timeNegativePr)) - evTimeMean; - float decayTimeK0Short = 0.5f * ((pTra.tofSignal() - timePositivePi) + (nTra.tofSignal() - timeNegativePi)) - evTimeMean; + float nSigmaXiLaPr = o2::aod::cascdata::kNoTOFValue; + float nSigmaXiLaPi = o2::aod::cascdata::kNoTOFValue; + float nSigmaXiPi = o2::aod::cascdata::kNoTOFValue; + float nSigmaOmLaPr = o2::aod::cascdata::kNoTOFValue; + float nSigmaOmLaPi = o2::aod::cascdata::kNoTOFValue; + float nSigmaOmKa = o2::aod::cascdata::kNoTOFValue; + }; + + struct trackTofInfo { // holds input track info + bool hasITS = false; + bool hasTPC = false; + bool hasTOF = false; + int collisionId = -1; + float tofExpMom = 0.0f; + float tofSignal = 0.0f; + float tofEvTime = 0.0f; + float length = 0.0f; + + // save TPC PID here for completeness too + float tpcNSigmaPi = 0.0f; + float tpcNSigmaKa = 0.0f; + float tpcNSigmaPr = 0.0f; + }; - float betaLambda = o2::aod::cascdata::kNoTOFValue; - float betaAntiLambda = o2::aod::cascdata::kNoTOFValue; - float betaK0Short = o2::aod::cascdata::kNoTOFValue; + // templatized process function for symmetric operation in derived and original AO2D + /// \param collisions the collisions table (needed for de-referencing V0 and progns) + /// \param v0 the V0 being processed + /// \param pTof the TOF information for the positive track + /// \param nTof the TOF information for the negative track + template + v0TofInfo calculateTofInfoV0(TCollisions const& collisions, int const& collisionId, TV0 const& v0, TTOFInfo const& pTof, TTOFInfo const& nTof) + { + v0TofInfo v0tof; // return this struct + auto collision = collisions.rawIteratorAt(collisionId); - if (nTra.hasTOF() && pTra.hasTOF()) { - betaLambda = (lengthV0 / decayTimeLambda) / 0.0299792458; - betaAntiLambda = (lengthV0 / decayTimeAntiLambda) / 0.0299792458; - betaK0Short = (lengthV0 / decayTimeK0Short) / 0.0299792458; - } + //_____________________________________________________________________________________________ + // daughter tracks: initialize from V0 position and momenta + o2::track::TrackPar posTrack = o2::track::TrackPar({v0.x(), v0.y(), v0.z()}, {v0.pxpos(), v0.pypos(), v0.pzpos()}, +1, false); + o2::track::TrackPar negTrack = o2::track::TrackPar({v0.x(), v0.y(), v0.z()}, {v0.pxneg(), v0.pyneg(), v0.pzneg()}, -1, false); - v0tofpid(deltaTimePositiveLambdaPi, deltaTimePositiveLambdaPr, - deltaTimeNegativeLambdaPi, deltaTimeNegativeLambdaPr, - deltaTimePositiveK0ShortPi, deltaTimeNegativeK0ShortPi, - deltaDecayTimeLambda, deltaDecayTimeAntiLambda, deltaDecayTimeK0Short); - v0tofbeta(betaLambda, betaAntiLambda, betaK0Short); - v0tofdebugs(timeLambda, timeK0Short, timePositivePr, timePositivePi, timeNegativePr, timeNegativePi); - - // do Nsigmas if requested - if (doNSigmas && nSigmaCalibLoaded) { - // sweep through all viable hypotheses and produce N-sigma - - if (deltaTimePositiveLambdaPi > -1e+5) - nSigmaPositiveLambdaPi = (deltaTimePositiveLambdaPi - hMeanPosLaPi->Interpolate(v0.p())) / hSigmaPosLaPi->Interpolate(v0.p()); - if (deltaTimePositiveLambdaPr > -1e+5) - nSigmaPositiveLambdaPr = (deltaTimePositiveLambdaPr - hMeanPosLaPr->Interpolate(v0.p())) / hSigmaPosLaPr->Interpolate(v0.p()); - if (deltaTimeNegativeLambdaPi > -1e+5) - nSigmaNegativeLambdaPi = (deltaTimeNegativeLambdaPi - hMeanNegLaPi->Interpolate(v0.p())) / hSigmaNegLaPi->Interpolate(v0.p()); - if (deltaTimeNegativeLambdaPr > -1e+5) - nSigmaNegativeLambdaPr = (deltaTimeNegativeLambdaPr - hMeanNegLaPr->Interpolate(v0.p())) / hSigmaNegLaPr->Interpolate(v0.p()); - if (deltaTimePositiveK0ShortPi > -1e+5) - nSigmaPositiveK0ShortPi = (deltaTimePositiveK0ShortPi - hMeanPosK0Pi->Interpolate(v0.p())) / hSigmaPosK0Pi->Interpolate(v0.p()); - if (deltaTimeNegativeK0ShortPi > -1e+5) - nSigmaNegativeK0ShortPi = (deltaTimeNegativeK0ShortPi - hMeanNegK0Pi->Interpolate(v0.p())) / hSigmaNegK0Pi->Interpolate(v0.p()); - - v0tofnsigmas( - nSigmaPositiveLambdaPr, nSigmaNegativeLambdaPi, - nSigmaNegativeLambdaPr, nSigmaPositiveLambdaPi, - nSigmaPositiveK0ShortPi, nSigmaNegativeK0ShortPi); - } + //_____________________________________________________________________________________________ + // time of V0 segment + float lengthV0 = std::hypot(v0.x() - collision.posX(), v0.y() - collision.posY(), v0.z() - collision.posZ()); + float velocityK0Short = velocity(v0.p(), o2::constants::physics::MassKaonNeutral); + float velocityLambda = velocity(v0.p(), o2::constants::physics::MassLambda); + float timeK0Short = lengthV0 / velocityK0Short; // in picoseconds + float timeLambda = lengthV0 / velocityLambda; // in picoseconds - if (doQA) { - // length factor due to eta (to offset e-loss) - float positiveCosine = 1.0f / sqrt(1.0f + posTrack.getTgl() * posTrack.getTgl()); - float negativeCosine = 1.0f / sqrt(1.0f + negTrack.getTgl() * negTrack.getTgl()); - if (propagationConfiguration.correctELossInclination.value == false) { - negativeCosine = positiveCosine = 1.0f; + //_____________________________________________________________________________________________ + // define simple checks + bool passesQAcuts = (v0.v0cosPA() > v0Group.qaCosPA && v0.dcaV0daughters() < v0Group.qaDCADau); + bool lambdaCandidate = std::abs(v0.mLambda() - o2::constants::physics::MassLambda) < v0Group.qaMassWindow && + std::abs(pTof.tpcNSigmaPr) < v0Group.qaTPCNSigma && + std::abs(nTof.tpcNSigmaPi) < v0Group.qaTPCNSigma; + bool antiLambdaCandidate = std::abs(v0.mAntiLambda() - o2::constants::physics::MassLambda) < v0Group.qaMassWindow && + std::abs(pTof.tpcNSigmaPi) < v0Group.qaTPCNSigma && + std::abs(nTof.tpcNSigmaPr) < v0Group.qaTPCNSigma; + bool k0ShortCandidate = std::abs(v0.mK0Short() - o2::constants::physics::MassKaonNeutral) < v0Group.qaMassWindow && + std::abs(pTof.tpcNSigmaPi) < v0Group.qaTPCNSigma && + std::abs(nTof.tpcNSigmaPi) < v0Group.qaTPCNSigma; + + //_____________________________________________________________________________________________ + // Actual calculation + if (pTof.hasTOF && pTof.hasITS) { + float velocityPositivePr, velocityPositivePi, lengthPositive; + velocityPositivePr = velocityPositivePi = lengthPositive = o2::aod::v0data::kNoTOFValue; + // method 0: legacy standalone without use of primary particle TOF + if (calculationMethod.value == 0) { + velocityPositivePr = velocity(posTrack.getP(), o2::constants::physics::MassProton); + velocityPositivePi = velocity(posTrack.getP(), o2::constants::physics::MassPionCharged); + lengthPositive = findInterceptLength(posTrack, d_bz); + v0tof.timePositivePr = lengthPositive / velocityPositivePr; + v0tof.timePositivePi = lengthPositive / velocityPositivePi; } + // method 1: correct primary particle TOF information + // length -> revise by removing travel length to primary vertex + // expected momentum -> kept as is for now, could correct at second stage + // use main method from TOF to calculate expected time + if (calculationMethod.value == 1) { + if (pTof.collisionId >= 0) { + auto trackCollision = collisions.rawIteratorAt(pTof.collisionId); + const o2::math_utils::Point3D trackVertex{trackCollision.posX(), trackCollision.posY(), trackCollision.posZ()}; + o2::track::TrackLTIntegral ltIntegral; + bool successPropag = o2::base::Propagator::Instance()->propagateToDCA(trackVertex, posTrack, d_bz, 2.f, o2::base::Propagator::MatCorrType::USEMatCorrNONE, nullptr, <Integral); + if (successPropag) { + lengthPositive = pTof.length - ltIntegral.getL(); + v0tof.timePositivePr = o2::framework::pid::tof::MassToExpTime(pTof.tofExpMom, lengthPositive, o2::constants::physics::MassProton * o2::constants::physics::MassProton); + v0tof.timePositivePi = o2::framework::pid::tof::MassToExpTime(pTof.tofExpMom, lengthPositive, o2::constants::physics::MassPionCharged * o2::constants::physics::MassPionCharged); + } + } + } + if (lengthPositive > 0.0f) { + v0tof.deltaTimePositiveLambdaPr = (pTof.tofSignal - pTof.tofEvTime) - (timeLambda + v0tof.timePositivePr); + v0tof.deltaTimePositiveLambdaPi = (pTof.tofSignal - pTof.tofEvTime) - (timeLambda + v0tof.timePositivePi); + v0tof.deltaTimePositiveK0ShortPi = (pTof.tofSignal - pTof.tofEvTime) - (timeK0Short + v0tof.timePositivePi); + + // de facto nsigma + if (nSigmaCalibLoaded) { + v0tof.nSigmaPositiveLambdaPi = (v0tof.deltaTimePositiveLambdaPi - hMeanPosLaPi->Interpolate(v0.p())) / hSigmaPosLaPi->Interpolate(v0.p()); + v0tof.nSigmaPositiveLambdaPr = (v0tof.deltaTimePositiveLambdaPr - hMeanPosLaPr->Interpolate(v0.p())) / hSigmaPosLaPr->Interpolate(v0.p()); + v0tof.nSigmaPositiveK0ShortPi = (v0tof.deltaTimePositiveK0ShortPi - hMeanPosK0Pi->Interpolate(v0.p())) / hSigmaPosK0Pi->Interpolate(v0.p()); + } - if (pTra.hasTOF() && pTra.hasITS()) { - if (v0.v0cosPA() > v0Group.qaCosPA && v0.dcaV0daughters() < v0Group.qaDCADau) { - if (std::abs(v0.mLambda() - 1.115683) < v0Group.qaMassWindow && fabs(pTra.tpcNSigmaPr()) < v0Group.qaTPCNSigma && fabs(nTra.tpcNSigmaPi()) < v0Group.qaTPCNSigma && ((v0pdg == 0) || (v0pdg == 3122))) { - histos.fill(HIST("h2dDeltaTimePositiveLambdaPr"), v0.p(), v0.eta(), deltaTimePositiveLambdaPr); - if (calculationMethod.value == 2 && std::abs(timePositivePr_Method0 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon && std::abs(timePositivePr_Method1 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon) { - histos.fill(HIST("hDeltaTimeMethodsVsP_posLaPr"), positiveP, v0.positiveeta(), (timePositivePr_Method0 - timePositivePr_Method1) * positiveCosine); - histos.fill(HIST("hRatioTimeMethodsVsP_posLaPr"), positiveP, v0.positiveeta(), (timePositivePr_Method1 / timePositivePr_Method0) * positiveCosine); + // do QA histograms (calibration / QC) + if (doQA) { + if (passesQAcuts) { + if (lambdaCandidate) { + histos.fill(HIST("h2dDeltaTimePositiveLambdaPr"), v0.p(), v0.eta(), v0tof.deltaTimePositiveLambdaPr); } - if (doQANSigma) - histos.fill(HIST("h2dNSigmaPositiveLambdaPr"), v0.p(), nSigmaPositiveLambdaPr); - } - if (std::abs(v0.mAntiLambda() - 1.115683) < v0Group.qaMassWindow && fabs(pTra.tpcNSigmaPi()) < v0Group.qaTPCNSigma && fabs(nTra.tpcNSigmaPr()) < v0Group.qaTPCNSigma && ((v0pdg == 0) || (v0pdg == -3122))) { - histos.fill(HIST("h2dDeltaTimePositiveLambdaPi"), v0.p(), v0.eta(), deltaTimePositiveLambdaPi); - if (calculationMethod.value == 2 && std::abs(timePositivePi_Method0 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon && std::abs(timePositivePi_Method1 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon) { - histos.fill(HIST("hDeltaTimeMethodsVsP_posLaPi"), positiveP, v0.positiveeta(), (timePositivePi_Method0 - timePositivePi_Method1) * positiveCosine); - histos.fill(HIST("hRatioTimeMethodsVsP_posLaPi"), positiveP, v0.positiveeta(), (timePositivePi_Method1 / timePositivePi_Method0) * positiveCosine); + if (antiLambdaCandidate) { + histos.fill(HIST("h2dDeltaTimePositiveLambdaPi"), v0.p(), v0.eta(), v0tof.deltaTimePositiveLambdaPi); } - if (doQANSigma) - histos.fill(HIST("h2dNSigmaPositiveLambdaPi"), v0.p(), nSigmaPositiveLambdaPi); - } - if (std::abs(v0.mK0Short() - 0.497) < v0Group.qaMassWindow && fabs(pTra.tpcNSigmaPi()) < v0Group.qaTPCNSigma && fabs(nTra.tpcNSigmaPi()) < v0Group.qaTPCNSigma && ((v0pdg == 0) || (v0pdg == 310))) { - histos.fill(HIST("h2dDeltaTimePositiveK0ShortPi"), v0.p(), v0.eta(), deltaTimePositiveK0ShortPi); - if (calculationMethod.value == 2 && std::abs(timePositivePi_Method0 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon && std::abs(timePositivePi_Method1 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon) { - histos.fill(HIST("hDeltaTimeMethodsVsP_posK0Pi"), positiveP, v0.positiveeta(), (timePositivePi_Method0 - timePositivePi_Method1) * positiveCosine); - histos.fill(HIST("hRatioTimeMethodsVsP_posK0Pi"), positiveP, v0.positiveeta(), (timePositivePi_Method1 / timePositivePi_Method0) * positiveCosine); + if (k0ShortCandidate) { + histos.fill(HIST("h2dDeltaTimePositiveK0ShortPi"), v0.p(), v0.eta(), v0tof.deltaTimePositiveK0ShortPi); } - if (doQANSigma) - histos.fill(HIST("h2dNSigmaPositiveK0ShortPi"), v0.p(), nSigmaPositiveK0ShortPi); } } } + } + if (nTof.hasTOF && nTof.hasITS) { + float velocityNegativePr, velocityNegativePi, lengthNegative; + velocityNegativePr = velocityNegativePi = lengthNegative = o2::aod::v0data::kNoTOFValue; + // method 0: legacy standalone without use of primary particle TOF + if (calculationMethod.value == 0) { + velocityNegativePr = velocity(negTrack.getP(), o2::constants::physics::MassProton); + velocityNegativePi = velocity(negTrack.getP(), o2::constants::physics::MassPionCharged); + lengthNegative = findInterceptLength(negTrack, d_bz); + v0tof.timeNegativePr = lengthNegative / velocityNegativePr; + v0tof.timeNegativePi = lengthNegative / velocityNegativePi; + } + // method 1: correct primary particle TOF information + // length -> revise by removing travel length to primary vertex + // expected momentum -> kept as is for now, could correct at second stage + // use main method from TOF to calculate expected time + if (calculationMethod.value == 1) { + if (nTof.collisionId >= 0) { + auto trackCollision = collisions.rawIteratorAt(nTof.collisionId); + const o2::math_utils::Point3D trackVertex{trackCollision.posX(), trackCollision.posY(), trackCollision.posZ()}; + o2::track::TrackLTIntegral ltIntegral; + bool successPropag = o2::base::Propagator::Instance()->propagateToDCA(trackVertex, negTrack, d_bz, 2.f, o2::base::Propagator::MatCorrType::USEMatCorrNONE, nullptr, <Integral); + if (successPropag) { + lengthNegative = nTof.length - ltIntegral.getL(); + v0tof.timeNegativePr = o2::framework::pid::tof::MassToExpTime(nTof.tofExpMom, nTof.length - ltIntegral.getL(), o2::constants::physics::MassProton * o2::constants::physics::MassProton); + v0tof.timeNegativePi = o2::framework::pid::tof::MassToExpTime(nTof.tofExpMom, nTof.length - ltIntegral.getL(), o2::constants::physics::MassPionCharged * o2::constants::physics::MassPionCharged); + } + } + } + if (lengthNegative > 0.0f) { + v0tof.deltaTimeNegativeLambdaPr = (nTof.tofSignal - nTof.tofEvTime) - (timeLambda + v0tof.timeNegativePr); + v0tof.deltaTimeNegativeLambdaPi = (nTof.tofSignal - nTof.tofEvTime) - (timeLambda + v0tof.timeNegativePi); + v0tof.deltaTimeNegativeK0ShortPi = (nTof.tofSignal - nTof.tofEvTime) - (timeK0Short + v0tof.timeNegativePi); + + // de facto nsigma + if (nSigmaCalibLoaded) { + v0tof.nSigmaNegativeLambdaPi = (v0tof.deltaTimeNegativeLambdaPi - hMeanNegLaPi->Interpolate(v0.p())) / hSigmaNegLaPi->Interpolate(v0.p()); + v0tof.nSigmaNegativeLambdaPr = (v0tof.deltaTimeNegativeLambdaPr - hMeanNegLaPr->Interpolate(v0.p())) / hSigmaNegLaPr->Interpolate(v0.p()); + v0tof.nSigmaNegativeK0ShortPi = (v0tof.deltaTimeNegativeK0ShortPi - hMeanNegK0Pi->Interpolate(v0.p())) / hSigmaNegK0Pi->Interpolate(v0.p()); + } - if (nTra.hasTOF() && nTra.hasITS()) { - if (v0.v0cosPA() > v0Group.qaCosPA && v0.dcaV0daughters() < v0Group.qaDCADau) { - if (std::abs(v0.mLambda() - 1.115683) < v0Group.qaMassWindow && fabs(pTra.tpcNSigmaPr()) < v0Group.qaTPCNSigma && fabs(nTra.tpcNSigmaPi()) < v0Group.qaTPCNSigma && ((v0pdg == 0) || (v0pdg == 3122))) { - histos.fill(HIST("h2dDeltaTimeNegativeLambdaPi"), v0.p(), v0.eta(), deltaTimeNegativeLambdaPi); - if (calculationMethod.value == 2 && std::abs(timeNegativePi_Method0 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon && std::abs(timeNegativePi_Method1 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon) { - histos.fill(HIST("hDeltaTimeMethodsVsP_negLaPi"), negativeP, v0.negativeeta(), (timeNegativePi_Method0 - timeNegativePi_Method1) * negativeCosine); - histos.fill(HIST("hRatioTimeMethodsVsP_negLaPi"), negativeP, v0.negativeeta(), (timeNegativePi_Method1 / timeNegativePi_Method0)); + // do QA histograms (calibration / QC) + if (doQA) { + if (passesQAcuts) { + if (lambdaCandidate) { + histos.fill(HIST("h2dDeltaTimeNegativeLambdaPi"), v0.p(), v0.eta(), v0tof.deltaTimeNegativeLambdaPi); } - // delta lambda decay time - histos.fill(HIST("h2dLambdaDeltaDecayTime"), v0.p(), deltaDecayTimeLambda); - if (doQANSigma) - histos.fill(HIST("h2dNSigmaNegativeLambdaPi"), v0.p(), nSigmaNegativeLambdaPi); - } - if (std::abs(v0.mAntiLambda() - 1.115683) < v0Group.qaMassWindow && fabs(pTra.tpcNSigmaPi()) < v0Group.qaTPCNSigma && fabs(nTra.tpcNSigmaPr()) < v0Group.qaTPCNSigma && ((v0pdg == 0) || (v0pdg == -3122))) { - histos.fill(HIST("h2dDeltaTimeNegativeLambdaPr"), v0.p(), v0.eta(), deltaTimeNegativeLambdaPr); - if (calculationMethod.value == 2 && std::abs(timeNegativePr_Method0 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon && std::abs(timeNegativePr_Method1 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon) { - histos.fill(HIST("hDeltaTimeMethodsVsP_negLaPr"), negativeP, v0.negativeeta(), (timeNegativePr_Method0 - timeNegativePr_Method1) * negativeCosine); - histos.fill(HIST("hRatioTimeMethodsVsP_negLaPr"), negativeP, v0.negativeeta(), (timeNegativePr_Method1 / timeNegativePr_Method0)); + if (antiLambdaCandidate) { + histos.fill(HIST("h2dDeltaTimeNegativeLambdaPr"), v0.p(), v0.eta(), v0tof.deltaTimeNegativeLambdaPr); } - if (doQANSigma) - histos.fill(HIST("h2dNSigmaNegativeLambdaPr"), v0.p(), nSigmaNegativeLambdaPr); - } - if (std::abs(v0.mK0Short() - 0.497) < v0Group.qaMassWindow && fabs(pTra.tpcNSigmaPi()) < v0Group.qaTPCNSigma && fabs(nTra.tpcNSigmaPi()) < v0Group.qaTPCNSigma && ((v0pdg == 0) || (v0pdg == 310))) { - histos.fill(HIST("h2dDeltaTimeNegativeK0ShortPi"), v0.p(), v0.eta(), deltaTimeNegativeK0ShortPi); - if (calculationMethod.value == 2 && std::abs(timeNegativePi_Method0 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon && std::abs(timeNegativePi_Method1 - o2::aod::v0data::kNoTOFValue) > o2::aod::v0data::kEpsilon) { - histos.fill(HIST("hDeltaTimeMethodsVsP_negK0Pi"), negativeP, v0.negativeeta(), (timeNegativePi_Method0 - timeNegativePi_Method1) * negativeCosine); - histos.fill(HIST("hRatioTimeMethodsVsP_negK0Pi"), negativeP, v0.negativeeta(), (timeNegativePi_Method1 / timeNegativePi_Method0)); + if (k0ShortCandidate) { + histos.fill(HIST("h2dDeltaTimeNegativeK0ShortPi"), v0.p(), v0.eta(), v0tof.deltaTimeNegativeK0ShortPi); } - if (doQANSigma) - histos.fill(HIST("h2dNSigmaNegativeK0ShortPi"), v0.p(), nSigmaNegativeK0ShortPi); } } } } - } - template - void processCascadeCandidate(TCollision const& collision, TCascade const& cascade, TTrack const& pTra, TTrack const& nTra, TTrack const& bTra, int cascpdg) + return v0tof; + } // end calculation altogether + + template + cascTofInfo calculateTofInfoCascade(TCollisions const& collisions, int const& collisionId, TCascade const& cascade, TTOFInfo const& pTof, TTOFInfo const& nTof, TTOFInfo const& bTof) { - // initialize from positions and momenta as needed + cascTofInfo casctof; // return this struct + auto collision = collisions.rawIteratorAt(collisionId); + + //_____________________________________________________________________________________________ + // daughter tracks: initialize from V0 position and momenta o2::track::TrackPar posTrack = o2::track::TrackPar({cascade.xlambda(), cascade.ylambda(), cascade.zlambda()}, {cascade.pxpos(), cascade.pypos(), cascade.pzpos()}, +1, false); o2::track::TrackPar negTrack = o2::track::TrackPar({cascade.xlambda(), cascade.ylambda(), cascade.zlambda()}, {cascade.pxneg(), cascade.pyneg(), cascade.pzneg()}, -1, false); o2::track::TrackPar bachTrack = o2::track::TrackPar({cascade.x(), cascade.y(), cascade.z()}, {cascade.pxbach(), cascade.pybach(), cascade.pzbach()}, cascade.sign(), false); o2::track::TrackPar cascTrack = o2::track::TrackPar({cascade.x(), cascade.y(), cascade.z()}, {cascade.px(), cascade.py(), cascade.pz()}, cascade.sign(), false); - float positiveP = std::hypot(cascade.pxpos(), cascade.pypos(), cascade.pzpos()); - float negativeP = std::hypot(cascade.pxneg(), cascade.pyneg(), cascade.pzneg()); - float bachelorP = std::hypot(cascade.pxbach(), cascade.pybach(), cascade.pzbach()); - - // start calculation: calculate velocities + //_____________________________________________________________________________________________ + // time of V0 segment float velocityXi = velocity(cascTrack.getP(), o2::constants::physics::MassXiMinus); float velocityOm = velocity(cascTrack.getP(), o2::constants::physics::MassOmegaMinus); float velocityLa = velocity(std::hypot(cascade.pxlambda(), cascade.pylambda(), cascade.pzlambda()), o2::constants::physics::MassLambda); - - // calculate mother lengths float lengthV0 = std::hypot(cascade.xlambda() - cascade.x(), cascade.ylambda() - cascade.y(), cascade.zlambda() - cascade.z()); float lengthCascade = o2::aod::cascdata::kNoTOFValue; - ; - const o2::math_utils::Point3D collVtx{collision.getX(), collision.getY(), collision.getZ()}; + + // cascade length (N.B. could be simpler via trackLTIntegral, kept with legacy calculation) + const o2::math_utils::Point3D collVtx{collision.posX(), collision.posY(), collision.posZ()}; bool successPropag = o2::base::Propagator::Instance()->propagateToDCA(collVtx, cascTrack, d_bz, 2.f, o2::base::Propagator::MatCorrType::USEMatCorrNONE); float d = -1.0f; - float linearToPV = std::hypot(cascade.x() - collision.getX(), cascade.y() - collision.getY(), cascade.z() - collision.getZ()); + float linearToPV = std::hypot(cascade.x() - collision.posX(), cascade.y() - collision.posY(), cascade.z() - collision.posZ()); if (successPropag) { std::array cascCloseToPVPosition; cascTrack.getXYZGlo(cascCloseToPVPosition); @@ -1136,349 +814,234 @@ struct strangenesstofpid { lengthCascade = linearToPV; // if propagation failed, use linear estimate (optional: actually do not define?) } - // lambda, xi and omega flight time is always defined + // flight times of decaying particles float lambdaFlight = lengthV0 / velocityLa; float xiFlight = lengthCascade / velocityXi; float omFlight = lengthCascade / velocityOm; - float posFlightPi = o2::aod::cascdata::kNoTOFValue; - float posFlightPr = o2::aod::cascdata::kNoTOFValue; - float negFlightPi = o2::aod::cascdata::kNoTOFValue; - float negFlightPr = o2::aod::cascdata::kNoTOFValue; - float bachFlightPi = o2::aod::cascdata::kNoTOFValue; - float bachFlightKa = o2::aod::cascdata::kNoTOFValue; - float posFlightPi_Method0 = o2::aod::cascdata::kNoTOFValue; - float posFlightPr_Method0 = o2::aod::cascdata::kNoTOFValue; - float negFlightPi_Method0 = o2::aod::cascdata::kNoTOFValue; - float negFlightPr_Method0 = o2::aod::cascdata::kNoTOFValue; - float bachFlightPi_Method0 = o2::aod::cascdata::kNoTOFValue; - float bachFlightKa_Method0 = o2::aod::cascdata::kNoTOFValue; - - float posFlightPi_Method1 = o2::aod::cascdata::kNoTOFValue; - float posFlightPr_Method1 = o2::aod::cascdata::kNoTOFValue; - float negFlightPi_Method1 = o2::aod::cascdata::kNoTOFValue; - float negFlightPr_Method1 = o2::aod::cascdata::kNoTOFValue; - float bachFlightPi_Method1 = o2::aod::cascdata::kNoTOFValue; - float bachFlightKa_Method1 = o2::aod::cascdata::kNoTOFValue; - - // actual time-of-flight of daughter calculation - if (calculationMethod.value == 0 || calculationMethod.value == 2) { - float velocityPositivePr = velocity(posTrack.getP(), o2::constants::physics::MassProton); - float velocityPositivePi = velocity(posTrack.getP(), o2::constants::physics::MassPionCharged); - float velocityNegativePr = velocity(negTrack.getP(), o2::constants::physics::MassProton); - float velocityNegativePi = velocity(negTrack.getP(), o2::constants::physics::MassPionCharged); - float velocityBachelorPi = velocity(bachTrack.getP(), o2::constants::physics::MassPionCharged); - float velocityBachelorKa = velocity(bachTrack.getP(), o2::constants::physics::MassKaonCharged); - - float lengthPositive = findInterceptLength(posTrack, d_bz); // FIXME: tofPosition ok? adjust? - float lengthNegative = findInterceptLength(negTrack, d_bz); // FIXME: tofPosition ok? adjust? - float lengthBachelor = findInterceptLength(bachTrack, d_bz); // FIXME: tofPosition ok? adjust? - - if (lengthPositive > 0 && pTra.hasTOF() && pTra.hasITS()) { - posFlightPi_Method0 = lengthPositive / velocityPositivePi; - posFlightPr_Method0 = lengthPositive / velocityPositivePr; - } - if (lengthNegative > 0 && nTra.hasTOF() && nTra.hasITS()) { - negFlightPi_Method0 = lengthNegative / velocityNegativePi; - negFlightPr_Method0 = lengthNegative / velocityNegativePr; - } - if (lengthBachelor > 0 && bTra.hasTOF() && bTra.hasITS()) { - bachFlightPi_Method0 = lengthBachelor / velocityBachelorPi; - bachFlightKa_Method0 = lengthBachelor / velocityBachelorKa; - } - } - - if (calculationMethod.value > 0) { - if (pTra.hasTOF() && pTra.hasITS()) { // calculate if signal present, otherwise skip - if (posTrack.getP() < propagationConfiguration.maxProtonMomentumForEloss.value && std::abs(pTra.tpcNSigmaPr()) < propagationConfiguration.tpcNsigmaThreshold) { - o2::track::TrackPar posTrackAsProton(posTrack); - posTrackAsProton.setPID(o2::track::PID::Proton); - calculateTOF(posTrackAsProton, posFlightPr_Method1); - } else { - posFlightPr_Method1 = posFlightPr_Method0; - } - - if (posTrack.getP() < propagationConfiguration.maxPionMomentumForEloss.value && std::abs(pTra.tpcNSigmaPi()) < propagationConfiguration.tpcNsigmaThreshold) { - o2::track::TrackPar posTrackAsPion(posTrack); - posTrackAsPion.setPID(o2::track::PID::Pion); - calculateTOF(posTrackAsPion, posFlightPi_Method1); - } else { - posFlightPi_Method1 = posFlightPi_Method0; - } + //_____________________________________________________________________________________________ + // define simple checks + bool passesQAcuts = (cascade.dcaV0daughters() < cascadeGroup.qaV0DCADau && cascade.dcacascdaughters() < cascadeGroup.qaCascDCADau && cascade.v0cosPA(collision.posX(), collision.posY(), collision.posZ()) > cascadeGroup.qaV0CosPA && cascade.casccosPA(collision.posX(), collision.posY(), collision.posZ()) > cascadeGroup.qaCascCosPA); + bool xiMinusCandidate = cascade.sign() < 0 && + std::abs(cascade.mXi() - o2::constants::physics::MassXiMinus) < cascadeGroup.qaMassWindow && + std::abs(pTof.tpcNSigmaPr) < cascadeGroup.qaTPCNSigma && + std::abs(nTof.tpcNSigmaPi) < cascadeGroup.qaTPCNSigma && + std::abs(bTof.tpcNSigmaPi) < cascadeGroup.qaTPCNSigma; + bool xiPlusCandidate = cascade.sign() > 0 && + std::abs(cascade.mXi() - o2::constants::physics::MassXiMinus) < cascadeGroup.qaMassWindow && + std::abs(pTof.tpcNSigmaPi) < cascadeGroup.qaTPCNSigma && + std::abs(nTof.tpcNSigmaPr) < cascadeGroup.qaTPCNSigma && + std::abs(bTof.tpcNSigmaPi) < cascadeGroup.qaTPCNSigma; + bool omegaMinusCandidate = cascade.sign() < 0 && + std::abs(cascade.mOmega() - o2::constants::physics::MassOmegaMinus) < cascadeGroup.qaMassWindow && + std::abs(pTof.tpcNSigmaPr) < cascadeGroup.qaTPCNSigma && + std::abs(nTof.tpcNSigmaPi) < cascadeGroup.qaTPCNSigma && + std::abs(bTof.tpcNSigmaKa) < cascadeGroup.qaTPCNSigma; + bool omegaPlusCandidate = cascade.sign() > 0 && + std::abs(cascade.mOmega() - o2::constants::physics::MassOmegaMinus) < cascadeGroup.qaMassWindow && + std::abs(pTof.tpcNSigmaPi) < cascadeGroup.qaTPCNSigma && + std::abs(nTof.tpcNSigmaPr) < cascadeGroup.qaTPCNSigma && + std::abs(bTof.tpcNSigmaKa) < cascadeGroup.qaTPCNSigma; + + //_____________________________________________________________________________________________ + // Actual calculation + if (pTof.hasTOF && pTof.hasITS) { + float velocityPositivePr, velocityPositivePi, lengthPositive; + velocityPositivePr = velocityPositivePi = lengthPositive = o2::aod::v0data::kNoTOFValue; + if (calculationMethod.value == 0) { + velocityPositivePr = velocity(posTrack.getP(), o2::constants::physics::MassProton); + velocityPositivePi = velocity(posTrack.getP(), o2::constants::physics::MassPionCharged); + lengthPositive = findInterceptLength(posTrack, d_bz); + casctof.posFlightPr = lengthPositive / velocityPositivePr; + casctof.posFlightPi = lengthPositive / velocityPositivePi; } - if (nTra.hasTOF() && nTra.hasITS()) { // calculate if signal present, otherwise skip - if (negTrack.getP() < propagationConfiguration.maxProtonMomentumForEloss.value && std::abs(nTra.tpcNSigmaPr()) < propagationConfiguration.tpcNsigmaThreshold) { - o2::track::TrackPar negTrackAsProton(negTrack); - negTrackAsProton.setPID(o2::track::PID::Proton); - calculateTOF(negTrackAsProton, negFlightPr_Method1); - } else { - negFlightPr_Method1 = negFlightPr_Method0; - } - - if (negTrack.getP() < propagationConfiguration.maxProtonMomentumForEloss.value && std::abs(nTra.tpcNSigmaPi()) < propagationConfiguration.tpcNsigmaThreshold) { - o2::track::TrackPar negTrackAsPion(negTrack); - negTrackAsPion.setPID(o2::track::PID::Pion); - calculateTOF(negTrackAsPion, negFlightPi_Method1); - } else { - negFlightPi_Method1 = negFlightPi_Method0; + // method 1: correct primary particle TOF information + // length -> revise by removing travel length to primary vertex + // expected momentum -> kept as is for now, could correct at second stage + // use main method from TOF to calculate expected time + if (calculationMethod.value == 1) { + if (pTof.collisionId >= 0) { + auto trackCollision = collisions.rawIteratorAt(pTof.collisionId); + const o2::math_utils::Point3D trackVertex{trackCollision.posX(), trackCollision.posY(), trackCollision.posZ()}; + o2::track::TrackLTIntegral ltIntegral; + bool successPropag = o2::base::Propagator::Instance()->propagateToDCA(trackVertex, posTrack, d_bz, 2.f, o2::base::Propagator::MatCorrType::USEMatCorrNONE, nullptr, <Integral); + if (successPropag) { + lengthPositive = pTof.length - ltIntegral.getL(); + casctof.posFlightPr = o2::framework::pid::tof::MassToExpTime(pTof.tofExpMom, pTof.length - ltIntegral.getL(), o2::constants::physics::MassProton * o2::constants::physics::MassProton); + casctof.posFlightPi = o2::framework::pid::tof::MassToExpTime(pTof.tofExpMom, pTof.length - ltIntegral.getL(), o2::constants::physics::MassPionCharged * o2::constants::physics::MassPionCharged); + } } } - if (bTra.hasTOF() && bTra.hasITS()) { // calculate if signal present, otherwise skip - if (bachTrack.getP() < propagationConfiguration.maxPionMomentumForEloss.value && std::abs(bTra.tpcNSigmaPi()) < propagationConfiguration.tpcNsigmaThreshold) { - o2::track::TrackPar bachTrackAsPion(bachTrack); - bachTrackAsPion.setPID(o2::track::PID::Pion); - calculateTOF(bachTrackAsPion, bachFlightPi_Method1); - } else { - bachFlightPi_Method1 = bachFlightPi_Method0; + if (lengthPositive > 0.0f) { + casctof.posDeltaTimeAsXiPi = (pTof.tofSignal - pTof.tofEvTime) - (xiFlight + lambdaFlight + casctof.posFlightPi); + casctof.posDeltaTimeAsXiPr = (pTof.tofSignal - pTof.tofEvTime) - (xiFlight + lambdaFlight + casctof.posFlightPr); + casctof.posDeltaTimeAsOmPi = (pTof.tofSignal - pTof.tofEvTime) - (omFlight + lambdaFlight + casctof.posFlightPi); + casctof.posDeltaTimeAsOmPr = (pTof.tofSignal - pTof.tofEvTime) - (omFlight + lambdaFlight + casctof.posFlightPr); + + // de facto nsigma + if (nSigmaCalibLoaded) { + if (cascade.sign() < 0) { + casctof.nSigmaXiLaPr = (casctof.posDeltaTimeAsXiPr - hMeanPosXiPr->Interpolate(cascade.p())) / hSigmaPosXiPr->Interpolate(cascade.p()); + casctof.nSigmaOmLaPr = (casctof.posDeltaTimeAsOmPr - hMeanPosOmPr->Interpolate(cascade.p())) / hSigmaPosOmPr->Interpolate(cascade.p()); + } else { + casctof.nSigmaXiLaPi = (casctof.posDeltaTimeAsXiPi - hMeanPosXiPi->Interpolate(cascade.p())) / hSigmaPosXiPi->Interpolate(cascade.p()); + casctof.nSigmaOmLaPi = (casctof.posDeltaTimeAsOmPi - hMeanPosOmPi->Interpolate(cascade.p())) / hSigmaPosOmPi->Interpolate(cascade.p()); + } } - if (bachTrack.getP() < propagationConfiguration.maxKaonMomentumForEloss.value && std::abs(bTra.tpcNSigmaKa()) < propagationConfiguration.tpcNsigmaThreshold) { - o2::track::TrackPar bachTrackAsKaon(bachTrack); - bachTrackAsKaon.setPID(o2::track::PID::Kaon); - calculateTOF(bachTrackAsKaon, bachFlightKa_Method1); - } else { - bachFlightKa_Method1 = bachFlightKa_Method0; + // do QA histograms (calibration / QC) + if (doQA) { + if (passesQAcuts) { + if (xiMinusCandidate) { + histos.fill(HIST("h2dposDeltaTimeAsXiPr"), cascade.p(), cascade.eta(), casctof.posDeltaTimeAsXiPr); + } + if (xiPlusCandidate) { + histos.fill(HIST("h2dposDeltaTimeAsXiPi"), cascade.p(), cascade.eta(), casctof.posDeltaTimeAsXiPi); + } + if (omegaMinusCandidate) { + histos.fill(HIST("h2dposDeltaTimeAsOmPr"), cascade.p(), cascade.eta(), casctof.posDeltaTimeAsOmPr); + } + if (omegaPlusCandidate) { + histos.fill(HIST("h2dposDeltaTimeAsOmPi"), cascade.p(), cascade.eta(), casctof.posDeltaTimeAsOmPi); + } + } } } - } - - // assign values to be used in main calculation - if (calculationMethod.value == 0) { - posFlightPi = posFlightPi_Method0; - posFlightPr = posFlightPr_Method0; - negFlightPi = negFlightPi_Method0; - negFlightPr = negFlightPr_Method0; - bachFlightPi = bachFlightPi_Method0; - bachFlightKa = bachFlightKa_Method0; - } else { - posFlightPi = posFlightPi_Method1; - posFlightPr = posFlightPr_Method1; - negFlightPi = negFlightPi_Method1; - negFlightPr = negFlightPr_Method1; - bachFlightPi = bachFlightPi_Method1; - bachFlightKa = bachFlightKa_Method1; - } - - if (calculationMethod.value == 2) { - // do analysis of successes and failures - bool positiveSuccessMethod0Pr = std::abs(posFlightPr_Method0 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon; - bool negativeSuccessMethod0Pr = std::abs(negFlightPr_Method0 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon; - bool positiveSuccessMethod1Pr = std::abs(posFlightPr_Method1 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon; - bool negativeSuccessMethod1Pr = std::abs(negFlightPr_Method1 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon; - bool positiveSuccessMethod0Pi = std::abs(posFlightPi_Method0 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon; - bool negativeSuccessMethod0Pi = std::abs(negFlightPi_Method0 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon; - bool positiveSuccessMethod1Pi = std::abs(posFlightPi_Method1 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon; - bool negativeSuccessMethod1Pi = std::abs(negFlightPi_Method1 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon; - - bool bachelorSuccessMethod0Pi = std::abs(bachFlightPi_Method0 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon; - bool bachelorSuccessMethod0Ka = std::abs(bachFlightKa_Method0 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon; - bool bachelorSuccessMethod1Pi = std::abs(bachFlightPi_Method1 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon; - bool bachelorSuccessMethod1Ka = std::abs(bachFlightKa_Method1 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon; - - int encodedPositiveSuccessPi = positiveSuccessMethod0Pi + 2 * positiveSuccessMethod1Pi; - int encodedPositiveSuccessPr = positiveSuccessMethod0Pr + 2 * positiveSuccessMethod1Pr; - int encodedNegativeSuccessPi = negativeSuccessMethod0Pi + 2 * negativeSuccessMethod1Pi; - int encodedNegativeSuccessPr = negativeSuccessMethod0Pr + 2 * negativeSuccessMethod1Pr; - int encodedBachelorSuccessPi = bachelorSuccessMethod0Pi + 2 * bachelorSuccessMethod1Pi; - int encodedBachelorSuccessKa = bachelorSuccessMethod0Ka + 2 * bachelorSuccessMethod1Ka; - - if (pTra.hasTOF()) { - histos.fill(HIST("h2dSucessRateProton"), positiveP, encodedPositiveSuccessPr); - histos.fill(HIST("h2dSucessRatePion"), positiveP, encodedPositiveSuccessPi); - } - if (nTra.hasTOF()) { - histos.fill(HIST("h2dSucessRateProton"), negativeP, encodedNegativeSuccessPr); - histos.fill(HIST("h2dSucessRatePion"), negativeP, encodedNegativeSuccessPi); + } // end positive + + if (nTof.hasTOF && nTof.hasITS) { + float velocityNegativePr, velocityNegativePi, lengthNegative; + velocityNegativePr = velocityNegativePi = lengthNegative = o2::aod::v0data::kNoTOFValue; + // method 0: legacy standalone without use of primary particle TOF + if (calculationMethod.value == 0) { + velocityNegativePr = velocity(negTrack.getP(), o2::constants::physics::MassProton); + velocityNegativePi = velocity(negTrack.getP(), o2::constants::physics::MassPionCharged); + lengthNegative = findInterceptLength(negTrack, d_bz); + casctof.negFlightPr = lengthNegative / velocityNegativePr; + casctof.negFlightPi = lengthNegative / velocityNegativePi; } - if (bTra.hasTOF()) { - histos.fill(HIST("h2dSucessRateKaon"), bachelorP, encodedBachelorSuccessKa); - histos.fill(HIST("h2dSucessRatePion"), bachelorP, encodedBachelorSuccessPi); - } - } - - // initialize delta-times (actual PID variables) - float posDeltaTimeAsXiPi = o2::aod::cascdata::kNoTOFValue, posDeltaTimeAsXiPr = o2::aod::cascdata::kNoTOFValue; - float negDeltaTimeAsXiPi = o2::aod::cascdata::kNoTOFValue, negDeltaTimeAsXiPr = o2::aod::cascdata::kNoTOFValue; - float bachDeltaTimeAsXiPi = o2::aod::cascdata::kNoTOFValue; - float posDeltaTimeAsOmPi = o2::aod::cascdata::kNoTOFValue, posDeltaTimeAsOmPr = o2::aod::cascdata::kNoTOFValue; - float negDeltaTimeAsOmPi = o2::aod::cascdata::kNoTOFValue, negDeltaTimeAsOmPr = o2::aod::cascdata::kNoTOFValue; - float bachDeltaTimeAsOmKa = o2::aod::cascdata::kNoTOFValue; - - if (pTra.hasTOF() && pTra.hasITS()) { - posDeltaTimeAsXiPi = (pTra.tofSignal() - pTra.tofEvTime()) - (xiFlight + lambdaFlight + posFlightPi); - posDeltaTimeAsXiPr = (pTra.tofSignal() - pTra.tofEvTime()) - (xiFlight + lambdaFlight + posFlightPr); - posDeltaTimeAsOmPi = (pTra.tofSignal() - pTra.tofEvTime()) - (omFlight + lambdaFlight + posFlightPi); - posDeltaTimeAsOmPr = (pTra.tofSignal() - pTra.tofEvTime()) - (omFlight + lambdaFlight + posFlightPr); - } - if (nTra.hasTOF() && nTra.hasITS()) { - negDeltaTimeAsXiPi = (nTra.tofSignal() - nTra.tofEvTime()) - (xiFlight + lambdaFlight + negFlightPi); - negDeltaTimeAsXiPr = (nTra.tofSignal() - nTra.tofEvTime()) - (xiFlight + lambdaFlight + negFlightPr); - negDeltaTimeAsOmPi = (nTra.tofSignal() - nTra.tofEvTime()) - (omFlight + lambdaFlight + negFlightPi); - negDeltaTimeAsOmPr = (nTra.tofSignal() - nTra.tofEvTime()) - (omFlight + lambdaFlight + negFlightPr); - } - if (bTra.hasTOF() && bTra.hasITS()) { - bachDeltaTimeAsXiPi = (bTra.tofSignal() - bTra.tofEvTime()) - (xiFlight + bachFlightPi); - bachDeltaTimeAsOmKa = (bTra.tofSignal() - bTra.tofEvTime()) - (omFlight + bachFlightKa); - } - - casctofpids( - posDeltaTimeAsXiPi, posDeltaTimeAsXiPr, negDeltaTimeAsXiPi, negDeltaTimeAsXiPr, bachDeltaTimeAsXiPi, - posDeltaTimeAsOmPi, posDeltaTimeAsOmPr, negDeltaTimeAsOmPi, negDeltaTimeAsOmPr, bachDeltaTimeAsOmKa); - - float nSigmaXiLaPr = o2::aod::cascdata::kNoTOFValue; - float nSigmaXiLaPi = o2::aod::cascdata::kNoTOFValue; - float nSigmaXiPi = o2::aod::cascdata::kNoTOFValue; - float nSigmaOmLaPr = o2::aod::cascdata::kNoTOFValue; - float nSigmaOmLaPi = o2::aod::cascdata::kNoTOFValue; - float nSigmaOmKa = o2::aod::cascdata::kNoTOFValue; - - // go for Nsigma values if requested - if (doNSigmas && nSigmaCalibLoaded) { - // Xi hypothesis ________________________ - if (cascade.sign() < 0) { // XiMinus - if (posDeltaTimeAsXiPr > -1e+5) // proton from Lambda from XiMinus has signal - nSigmaXiLaPr = (posDeltaTimeAsXiPr - hMeanPosXiPr->Interpolate(cascade.p())) / hSigmaPosXiPr->Interpolate(cascade.p()); - if (negDeltaTimeAsXiPi > -1e+5) // pion from Lambda from XiMinus has signal - nSigmaXiLaPi = (negDeltaTimeAsXiPi - hMeanNegXiPi->Interpolate(cascade.p())) / hSigmaNegXiPi->Interpolate(cascade.p()); - if (bachDeltaTimeAsXiPi > -1e+5) // pion from XiMinus has signal - nSigmaXiPi = (bachDeltaTimeAsXiPi - hMeanBachXiPi->Interpolate(cascade.p())) / hSigmaBachXiPi->Interpolate(cascade.p()); - if (posDeltaTimeAsOmPr > -1e+5) // proton from Lambda from OmegaMinus has signal - nSigmaOmLaPr = (posDeltaTimeAsOmPr - hMeanPosOmPr->Interpolate(cascade.p())) / hSigmaPosOmPr->Interpolate(cascade.p()); - if (negDeltaTimeAsOmPi > -1e+5) // pion from Lambda from OmegaMinus has signal - nSigmaOmLaPi = (negDeltaTimeAsOmPi - hMeanNegOmPi->Interpolate(cascade.p())) / hSigmaNegOmPi->Interpolate(cascade.p()); - if (bachDeltaTimeAsOmKa > -1e+5) // kaon from OmegaMinus has signal - nSigmaOmKa = (bachDeltaTimeAsOmKa - hMeanBachOmKa->Interpolate(cascade.p())) / hSigmaBachOmKa->Interpolate(cascade.p()); - } else { - if (posDeltaTimeAsXiPi > -1e+5) // proton from Lambda from XiMinus has signal - nSigmaXiLaPi = (posDeltaTimeAsXiPi - hMeanPosXiPi->Interpolate(cascade.p())) / hSigmaPosXiPi->Interpolate(cascade.p()); - if (negDeltaTimeAsXiPr > -1e+5) // pion from Lambda from XiMinus has signal - nSigmaXiLaPr = (negDeltaTimeAsXiPr - hMeanNegXiPr->Interpolate(cascade.p())) / hSigmaNegXiPr->Interpolate(cascade.p()); - if (bachDeltaTimeAsXiPi > -1e+5) // pion from XiMinus has signal - nSigmaXiPi = (bachDeltaTimeAsXiPi - hMeanBachXiPi->Interpolate(cascade.p())) / hSigmaBachXiPi->Interpolate(cascade.p()); - if (posDeltaTimeAsOmPi > -1e+5) // proton from Lambda from OmegaMinus has signal - nSigmaOmLaPi = (posDeltaTimeAsOmPi - hMeanPosOmPi->Interpolate(cascade.p())) / hSigmaPosOmPi->Interpolate(cascade.p()); - if (negDeltaTimeAsOmPr > -1e+5) // pion from Lambda from OmegaMinus has signal - nSigmaOmLaPr = (negDeltaTimeAsOmPr - hMeanNegOmPr->Interpolate(cascade.p())) / hSigmaNegOmPr->Interpolate(cascade.p()); - if (bachDeltaTimeAsOmKa > -1e+5) // kaon from OmegaMinus has signal - nSigmaOmKa = (bachDeltaTimeAsOmKa - hMeanBachOmKa->Interpolate(cascade.p())) / hSigmaBachOmKa->Interpolate(cascade.p()); - } - casctofnsigmas(nSigmaXiLaPi, nSigmaXiLaPr, nSigmaXiPi, nSigmaOmLaPi, nSigmaOmLaPr, nSigmaOmKa); - } - - if (doQA) { - // length factor due to eta (to offset e-loss) - float positiveCosine = 1.0f / sqrt(1.0f + posTrack.getTgl() * posTrack.getTgl()); - float negativeCosine = 1.0f / sqrt(1.0f + negTrack.getTgl() * negTrack.getTgl()); - float bachelorCosine = 1.0f / sqrt(1.0f + bachTrack.getTgl() * bachTrack.getTgl()); - if (propagationConfiguration.correctELossInclination.value == false) { - negativeCosine = positiveCosine = bachelorCosine = 1.0f; + // method 1: correct primary particle TOF information + // length -> revise by removing travel length to primary vertex + // expected momentum -> kept as is for now, could correct at second stage + // use main method from TOF to calculate expected time + if (calculationMethod.value == 1) { + if (nTof.collisionId >= 0) { + auto trackCollision = collisions.rawIteratorAt(nTof.collisionId); + const o2::math_utils::Point3D trackVertex{trackCollision.posX(), trackCollision.posY(), trackCollision.posZ()}; + o2::track::TrackLTIntegral ltIntegral; + bool successPropag = o2::base::Propagator::Instance()->propagateToDCA(trackVertex, negTrack, d_bz, 2.f, o2::base::Propagator::MatCorrType::USEMatCorrNONE, nullptr, <Integral); + if (successPropag) { + lengthNegative = nTof.length - ltIntegral.getL(); + casctof.negFlightPr = o2::framework::pid::tof::MassToExpTime(nTof.tofExpMom, nTof.length - ltIntegral.getL(), o2::constants::physics::MassProton * o2::constants::physics::MassProton); + casctof.negFlightPi = o2::framework::pid::tof::MassToExpTime(nTof.tofExpMom, nTof.length - ltIntegral.getL(), o2::constants::physics::MassPionCharged * o2::constants::physics::MassPionCharged); + } + } } + if (lengthNegative > 0.0f) { + casctof.negDeltaTimeAsXiPi = (nTof.tofSignal - nTof.tofEvTime) - (xiFlight + lambdaFlight + casctof.negFlightPi); + casctof.negDeltaTimeAsXiPr = (nTof.tofSignal - nTof.tofEvTime) - (xiFlight + lambdaFlight + casctof.negFlightPr); + casctof.negDeltaTimeAsOmPi = (nTof.tofSignal - nTof.tofEvTime) - (omFlight + lambdaFlight + casctof.negFlightPi); + casctof.negDeltaTimeAsOmPr = (nTof.tofSignal - nTof.tofEvTime) - (omFlight + lambdaFlight + casctof.negFlightPr); + + // de facto nsigma + if (nSigmaCalibLoaded) { + if (cascade.sign() < 0) { + casctof.nSigmaXiLaPr = (casctof.negDeltaTimeAsXiPr - hMeanPosXiPr->Interpolate(cascade.p())) / hSigmaPosXiPr->Interpolate(cascade.p()); + casctof.nSigmaOmLaPr = (casctof.negDeltaTimeAsOmPr - hMeanPosOmPr->Interpolate(cascade.p())) / hSigmaPosOmPr->Interpolate(cascade.p()); + } else { + casctof.nSigmaXiLaPi = (casctof.negDeltaTimeAsXiPi - hMeanPosXiPi->Interpolate(cascade.p())) / hSigmaPosXiPi->Interpolate(cascade.p()); + casctof.nSigmaOmLaPi = (casctof.negDeltaTimeAsOmPi - hMeanPosOmPi->Interpolate(cascade.p())) / hSigmaPosOmPi->Interpolate(cascade.p()); + } + } - if (cascade.dcaV0daughters() < cascadeGroup.qaV0DCADau && cascade.dcacascdaughters() < cascadeGroup.qaCascDCADau && cascade.v0cosPA(collision.getX(), collision.getY(), collision.getZ()) > cascadeGroup.qaV0CosPA && cascade.casccosPA(collision.getX(), collision.getY(), collision.getZ()) > cascadeGroup.qaCascCosPA) { - if (cascade.sign() < 0) { - if (std::abs(cascade.mXi() - 1.32171) < cascadeGroup.qaMassWindow && fabs(pTra.tpcNSigmaPr()) < cascadeGroup.qaTPCNSigma && fabs(nTra.tpcNSigmaPi()) < cascadeGroup.qaTPCNSigma && fabs(bTra.tpcNSigmaPi()) < cascadeGroup.qaTPCNSigma && ((cascpdg == 0) || (cascpdg == 3312))) { - histos.fill(HIST("h2dposDeltaTimeAsXiPr"), cascade.p(), cascade.eta(), posDeltaTimeAsXiPr); - histos.fill(HIST("h2dnegDeltaTimeAsXiPi"), cascade.p(), cascade.eta(), negDeltaTimeAsXiPi); - histos.fill(HIST("h2dbachDeltaTimeAsXiPi"), cascade.p(), cascade.eta(), bachDeltaTimeAsXiPi); - if (calculationMethod.value == 2) { - if (std::abs(posFlightPr_Method0 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon && std::abs(posFlightPr_Method1 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon) { - histos.fill(HIST("hDeltaTimeMethodsVsP_posXiPr"), positiveP, cascade.positiveeta(), (posFlightPr_Method0 - posFlightPr_Method1) * positiveCosine); - histos.fill(HIST("hRatioTimeMethodsVsP_posXiPr"), positiveP, cascade.positiveeta(), (posFlightPr_Method1 / posFlightPr_Method0)); - } - if (std::abs(negFlightPi_Method0 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon && std::abs(negFlightPi_Method1 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon) { - histos.fill(HIST("hDeltaTimeMethodsVsP_negXiPi"), negativeP, cascade.negativeeta(), (negFlightPi_Method0 - negFlightPi_Method1) * negativeCosine); - histos.fill(HIST("hRatioTimeMethodsVsP_negXiPi"), negativeP, cascade.negativeeta(), (negFlightPi_Method1 / negFlightPi_Method0)); - } - if (std::abs(bachFlightPi_Method0 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon && std::abs(bachFlightPi_Method1 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon) { - histos.fill(HIST("hDeltaTimeMethodsVsP_bachXiPi"), bachelorP, cascade.bacheloreta(), (bachFlightPi_Method0 - bachFlightPi_Method1) * bachelorCosine); - histos.fill(HIST("hRatioTimeMethodsVsP_bachXiPi"), bachelorP, cascade.bacheloreta(), (bachFlightPi_Method1 / bachFlightPi_Method0)); - } + // do QA histograms (calibration / QC) + if (doQA) { + if (passesQAcuts) { + if (xiMinusCandidate) { + histos.fill(HIST("h2dnegDeltaTimeAsXiPi"), cascade.p(), cascade.eta(), casctof.negDeltaTimeAsXiPi); } - if (doQANSigma) { - histos.fill(HIST("h2dNSigmaXiLaPi"), cascade.p(), nSigmaXiLaPi); - histos.fill(HIST("h2dNSigmaXiLaPr"), cascade.p(), nSigmaXiLaPr); - histos.fill(HIST("h2dNSigmaXiPi"), cascade.p(), nSigmaXiPi); + if (xiPlusCandidate) { + histos.fill(HIST("h2dnegDeltaTimeAsXiPr"), cascade.p(), cascade.eta(), casctof.negDeltaTimeAsXiPr); } - } - if (std::abs(cascade.mOmega() - 1.67245) < cascadeGroup.qaMassWindow && fabs(pTra.tpcNSigmaPr()) < cascadeGroup.qaTPCNSigma && fabs(nTra.tpcNSigmaPi()) < cascadeGroup.qaTPCNSigma && fabs(bTra.tpcNSigmaKa()) < cascadeGroup.qaTPCNSigma && ((cascpdg == 0) || (cascpdg == 3334))) { - histos.fill(HIST("h2dposDeltaTimeAsOmPr"), cascade.p(), cascade.eta(), posDeltaTimeAsOmPr); - histos.fill(HIST("h2dnegDeltaTimeAsOmPi"), cascade.p(), cascade.eta(), negDeltaTimeAsOmPi); - histos.fill(HIST("h2dbachDeltaTimeAsOmKa"), cascade.p(), cascade.eta(), bachDeltaTimeAsOmKa); - if (calculationMethod.value == 2) { - if (std::abs(posFlightPr_Method0 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon && std::abs(posFlightPr_Method1 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon) { - histos.fill(HIST("hDeltaTimeMethodsVsP_posOmPr"), positiveP, cascade.positiveeta(), (posFlightPr_Method0 - posFlightPr_Method1) * positiveCosine); - histos.fill(HIST("hRatioTimeMethodsVsP_posOmPr"), positiveP, cascade.positiveeta(), (posFlightPr_Method1 / posFlightPr_Method0)); - } - if (std::abs(negFlightPi_Method0 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon && std::abs(negFlightPi_Method1 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon) { - histos.fill(HIST("hDeltaTimeMethodsVsP_negOmPi"), negativeP, cascade.negativeeta(), (negFlightPi_Method0 - negFlightPi_Method1) * negativeCosine); - histos.fill(HIST("hRatioTimeMethodsVsP_negOmPi"), negativeP, cascade.negativeeta(), (negFlightPi_Method1 / negFlightPi_Method0)); - } - if (std::abs(bachFlightKa_Method0 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon && std::abs(bachFlightKa_Method1 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon) { - histos.fill(HIST("hDeltaTimeMethodsVsP_bachOmKa"), bachelorP, cascade.bacheloreta(), (bachFlightKa_Method0 - bachFlightKa_Method1) * bachelorCosine); - histos.fill(HIST("hRatioTimeMethodsVsP_bachOmKa"), bachelorP, cascade.bacheloreta(), (bachFlightKa_Method1 / bachFlightKa_Method0)); - } + if (omegaMinusCandidate) { + histos.fill(HIST("h2dnegDeltaTimeAsOmPi"), cascade.p(), cascade.eta(), casctof.negDeltaTimeAsOmPi); } - if (doQANSigma) { - histos.fill(HIST("h2dNSigmaOmLaPi"), cascade.p(), nSigmaOmLaPi); - histos.fill(HIST("h2dNSigmaOmLaPr"), cascade.p(), nSigmaOmLaPr); - histos.fill(HIST("h2dNSigmaOmKa"), cascade.p(), nSigmaOmKa); + if (omegaPlusCandidate) { + histos.fill(HIST("h2dnegDeltaTimeAsOmPr"), cascade.p(), cascade.eta(), casctof.negDeltaTimeAsOmPr); } } - } else { - if (std::abs(cascade.mXi() - 1.32171) < cascadeGroup.qaMassWindow && fabs(pTra.tpcNSigmaPi()) < cascadeGroup.qaTPCNSigma && fabs(nTra.tpcNSigmaPr()) < cascadeGroup.qaTPCNSigma && fabs(bTra.tpcNSigmaPi()) < cascadeGroup.qaTPCNSigma && ((cascpdg == 0) || (cascpdg == -3312))) { - histos.fill(HIST("h2dposDeltaTimeAsXiPi"), cascade.p(), cascade.eta(), posDeltaTimeAsXiPi); - histos.fill(HIST("h2dnegDeltaTimeAsXiPr"), cascade.p(), cascade.eta(), negDeltaTimeAsXiPr); - histos.fill(HIST("h2dbachDeltaTimeAsXiPi"), cascade.p(), cascade.eta(), bachDeltaTimeAsXiPi); - if (calculationMethod.value == 2) { - if (std::abs(posFlightPi_Method0 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon && std::abs(posFlightPi_Method1 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon) { - histos.fill(HIST("hDeltaTimeMethodsVsP_posXiPi"), positiveP, cascade.positiveeta(), (posFlightPi_Method0 - posFlightPi_Method1) * positiveCosine); - histos.fill(HIST("hRatioTimeMethodsVsP_posXiPi"), positiveP, cascade.positiveeta(), (posFlightPi_Method1 / posFlightPi_Method1)); - } - if (std::abs(negFlightPr_Method0 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon && std::abs(negFlightPr_Method1 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon) { - histos.fill(HIST("hDeltaTimeMethodsVsP_negXiPr"), negativeP, cascade.negativeeta(), (negFlightPr_Method0 - negFlightPr_Method1) * negativeCosine); - histos.fill(HIST("hRatioTimeMethodsVsP_negXiPr"), negativeP, cascade.negativeeta(), (negFlightPr_Method1 / negFlightPr_Method0)); - } - if (std::abs(bachFlightPi_Method0 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon && std::abs(bachFlightPi_Method1 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon) { - histos.fill(HIST("hDeltaTimeMethodsVsP_bachXiPi"), bachelorP, cascade.bacheloreta(), (bachFlightPi_Method0 - bachFlightPi_Method1) * bachelorCosine); - histos.fill(HIST("hRatioTimeMethodsVsP_bachXiPi"), bachelorP, cascade.bacheloreta(), (bachFlightPi_Method1 / bachFlightPi_Method0)); - } + } + } + } // end negative + + if (bTof.hasTOF && bTof.hasITS) { + float velocityBachelorKa, velocityBachelorPi, lengthBachelor; + velocityBachelorKa = velocityBachelorPi = lengthBachelor = o2::aod::v0data::kNoTOFValue; + // method 0: legacy standalone without use of primary particle TOF + if (calculationMethod.value == 0) { + velocityBachelorPi = velocity(bachTrack.getP(), o2::constants::physics::MassPionCharged); + velocityBachelorKa = velocity(bachTrack.getP(), o2::constants::physics::MassKaonCharged); + lengthBachelor = findInterceptLength(bachTrack, d_bz); + casctof.bachFlightPi = lengthBachelor / velocityBachelorPi; + casctof.bachFlightKa = lengthBachelor / velocityBachelorKa; + } + // method 1: correct primary particle TOF information + // length -> revise by removing travel length to primary vertex + // expected momentum -> kept as is for now, could correct at second stage + // use main method from TOF to calculate expected time + if (calculationMethod.value == 1) { + if (bTof.collisionId >= 0) { + auto trackCollision = collisions.rawIteratorAt(bTof.collisionId); + const o2::math_utils::Point3D trackVertex{trackCollision.posX(), trackCollision.posY(), trackCollision.posZ()}; + o2::track::TrackLTIntegral ltIntegral; + bool successPropag = o2::base::Propagator::Instance()->propagateToDCA(trackVertex, bachTrack, d_bz, 2.f, o2::base::Propagator::MatCorrType::USEMatCorrNONE, nullptr, <Integral); + if (successPropag) { + lengthBachelor = bTof.length - ltIntegral.getL(); + casctof.bachFlightPi = o2::framework::pid::tof::MassToExpTime(bTof.tofExpMom, bTof.length - ltIntegral.getL(), o2::constants::physics::MassPionCharged * o2::constants::physics::MassPionCharged); + casctof.bachFlightKa = o2::framework::pid::tof::MassToExpTime(bTof.tofExpMom, bTof.length - ltIntegral.getL(), o2::constants::physics::MassKaonCharged * o2::constants::physics::MassKaonCharged); + } + } + } + if (lengthBachelor > 0.0f) { + casctof.bachDeltaTimeAsXiPi = (bTof.tofSignal - bTof.tofEvTime) - (xiFlight + casctof.bachFlightPi); + casctof.bachDeltaTimeAsOmKa = (bTof.tofSignal - bTof.tofEvTime) - (omFlight + casctof.bachFlightKa); + + // de facto nsigma + if (nSigmaCalibLoaded) { + if (cascade.sign() < 0) { + casctof.nSigmaXiPi = (casctof.bachDeltaTimeAsXiPi - hMeanBachXiPi->Interpolate(cascade.p())) / hSigmaBachXiPi->Interpolate(cascade.p()); + casctof.nSigmaOmKa = (casctof.bachDeltaTimeAsOmKa - hMeanBachOmKa->Interpolate(cascade.p())) / hSigmaBachOmKa->Interpolate(cascade.p()); + } else { + casctof.nSigmaXiPi = (casctof.bachDeltaTimeAsXiPi - hMeanBachXiPi->Interpolate(cascade.p())) / hSigmaBachXiPi->Interpolate(cascade.p()); + casctof.nSigmaOmKa = (casctof.bachDeltaTimeAsOmKa - hMeanBachOmKa->Interpolate(cascade.p())) / hSigmaBachOmKa->Interpolate(cascade.p()); + } + } + + // do QA histograms (calibration / QC) + if (doQA) { + if (passesQAcuts) { + if (xiMinusCandidate) { + histos.fill(HIST("h2dbachDeltaTimeAsXiPi"), cascade.p(), cascade.eta(), casctof.bachDeltaTimeAsXiPi); } - if (doQANSigma) { - histos.fill(HIST("h2dNSigmaXiLaPi"), cascade.p(), nSigmaXiLaPi); - histos.fill(HIST("h2dNSigmaXiLaPr"), cascade.p(), nSigmaXiLaPr); - histos.fill(HIST("h2dNSigmaXiPi"), cascade.p(), nSigmaXiPi); + if (xiPlusCandidate) { + histos.fill(HIST("h2dbachDeltaTimeAsXiPi"), cascade.p(), cascade.eta(), casctof.bachDeltaTimeAsXiPi); } - } - if (std::abs(cascade.mOmega() - 1.67245) < cascadeGroup.qaMassWindow && fabs(pTra.tpcNSigmaPi()) < cascadeGroup.qaTPCNSigma && fabs(nTra.tpcNSigmaPr()) < cascadeGroup.qaTPCNSigma && fabs(bTra.tpcNSigmaKa()) < cascadeGroup.qaTPCNSigma && ((cascpdg == 0) || (cascpdg == -3334))) { - histos.fill(HIST("h2dposDeltaTimeAsOmPi"), cascade.p(), cascade.eta(), posDeltaTimeAsOmPi); - histos.fill(HIST("h2dnegDeltaTimeAsOmPr"), cascade.p(), cascade.eta(), negDeltaTimeAsOmPr); - histos.fill(HIST("h2dbachDeltaTimeAsOmKa"), cascade.p(), cascade.eta(), bachDeltaTimeAsOmKa); - if (calculationMethod.value == 2) { - if (std::abs(posFlightPi_Method0 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon && std::abs(posFlightPi_Method1 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon) { - histos.fill(HIST("hDeltaTimeMethodsVsP_posOmPi"), positiveP, cascade.positiveeta(), (posFlightPi_Method0 - posFlightPi_Method1) * positiveCosine); - histos.fill(HIST("hRatioTimeMethodsVsP_posOmPi"), positiveP, cascade.positiveeta(), (posFlightPi_Method1 / posFlightPi_Method1)); - } - if (std::abs(negFlightPr_Method0 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon && std::abs(negFlightPr_Method1 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon) { - histos.fill(HIST("hDeltaTimeMethodsVsP_negOmPr"), negativeP, cascade.negativeeta(), (negFlightPr_Method0 - negFlightPr_Method1) * negativeCosine); - histos.fill(HIST("hRatioTimeMethodsVsP_negOmPr"), negativeP, cascade.negativeeta(), (negFlightPr_Method1 / negFlightPr_Method0)); - } - if (std::abs(bachFlightKa_Method0 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon && std::abs(bachFlightKa_Method1 - o2::aod::cascdata::kNoTOFValue) > o2::aod::cascdata::kEpsilon) { - histos.fill(HIST("hDeltaTimeMethodsVsP_bachOmKa"), bachelorP, cascade.bacheloreta(), (bachFlightKa_Method0 - bachFlightKa_Method1) * bachelorCosine); - histos.fill(HIST("hRatioTimeMethodsVsP_bachOmKa"), bachelorP, cascade.bacheloreta(), (bachFlightKa_Method1 / bachFlightKa_Method1)); - } + if (omegaMinusCandidate) { + histos.fill(HIST("h2dbachDeltaTimeAsOmKa"), cascade.p(), cascade.eta(), casctof.bachDeltaTimeAsOmKa); } - if (doQANSigma) { - histos.fill(HIST("h2dNSigmaOmLaPi"), cascade.p(), nSigmaOmLaPi); - histos.fill(HIST("h2dNSigmaOmLaPr"), cascade.p(), nSigmaOmLaPr); - histos.fill(HIST("h2dNSigmaOmKa"), cascade.p(), nSigmaOmKa); + if (omegaPlusCandidate) { + histos.fill(HIST("h2dbachDeltaTimeAsOmKa"), cascade.p(), cascade.eta(), casctof.bachDeltaTimeAsOmKa); } } } } - } + } // end bachelor + + // don't forget to give feedback + return casctof; } - void processStandardData(aod::Collisions const& collisions, V0OriginalDatas const& V0s, CascOriginalDatas const& cascades, TracksWithAllExtras const&, aod::BCsWithTimestamps const& /*bcs*/) + void processStandardData(aod::Collisions const& collisions, V0OriginalDatas const& V0s, CascOriginalDatas const& cascades, TracksWithAllExtras const& tracks, aod::BCsWithTimestamps const& /*bcs*/) { // Fire up CCDB with first collision in record. If no collisions, bypass if (useCustomRunNumber || collisions.size() < 1) { @@ -1489,46 +1052,125 @@ struct strangenesstofpid { initCCDB(bc.runNumber()); } + //________________________________________________________________________ + // estimate event times (only necessary for original data) + std::vector collisionEventTime(collisions.size(), 0.0); + std::vector collisionNtracks(collisions.size(), 0); + for (const auto& track : tracks) { + if (track.hasTOF()) { + collisionEventTime[track.collisionId()] += track.tofEvTime(); + collisionNtracks[track.collisionId()]++; + } + } + for (const auto& collision : collisions) { + if (collisionNtracks[collision.globalIndex()] > 0) { + collisionEventTime[collision.globalIndex()] /= static_cast(collisionNtracks[collision.globalIndex()]); + } else { + collisionEventTime[collision.globalIndex()] = -1e+6; // undefined + } + } + if (calculateV0s.value) { for (const auto& V0 : V0s) { - // for storing whatever is the relevant quantity for the PV - o2::dataformats::VertexBase primaryVertex; - if (V0.has_collision()) { - auto const& collision = V0.collision(); - primaryVertex.setPos({collision.posX(), collision.posY(), collision.posZ()}); - primaryVertex.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()); - } else { - primaryVertex.setPos({mVtx->getX(), mVtx->getY(), mVtx->getZ()}); - } + trackTofInfo pTof, nTof; // information storage auto pTra = V0.posTrack_as(); auto nTra = V0.negTrack_as(); - processV0Candidate(primaryVertex, V0, pTra, nTra, 0); + + pTof.collisionId = pTra.collisionId(); + pTof.hasITS = pTra.hasITS(); + pTof.hasTPC = pTra.hasTPC(); + pTof.hasTOF = pTra.hasTOF(); + pTof.tofExpMom = pTra.tofExpMom(); + pTof.tofEvTime = collisionEventTime[V0.collisionId()]; + pTof.tofSignal = pTra.tofSignal(); + pTof.length = pTra.length(); + pTof.tpcNSigmaPi = pTra.tpcNSigmaPi(); + pTof.tpcNSigmaPr = pTra.tpcNSigmaPr(); + + nTof.collisionId = nTra.collisionId(); + nTof.hasITS = nTra.hasITS(); + nTof.hasTPC = nTra.hasTPC(); + nTof.hasTOF = nTra.hasTOF(); + nTof.tofExpMom = nTra.tofExpMom(); + nTof.tofEvTime = collisionEventTime[V0.collisionId()]; + nTof.tofSignal = nTra.tofSignal(); + nTof.length = nTra.length(); + nTof.tpcNSigmaPi = nTra.tpcNSigmaPi(); + nTof.tpcNSigmaPr = nTra.tpcNSigmaPr(); + + v0TofInfo v0tof = calculateTofInfoV0(collisions, V0.collisionId(), V0, pTof, nTof); + + if (doNSigmas) { + v0tofnsigmas( + v0tof.nSigmaPositiveLambdaPr, v0tof.nSigmaNegativeLambdaPi, + v0tof.nSigmaNegativeLambdaPr, v0tof.nSigmaPositiveLambdaPi, + v0tof.nSigmaPositiveK0ShortPi, v0tof.nSigmaNegativeK0ShortPi); + } } } if (calculateCascades.value) { for (const auto& cascade : cascades) { - // for storing whatever is the relevant quantity for the PV - o2::dataformats::VertexBase primaryVertex; - if (cascade.has_collision()) { - auto const& collision = cascade.collision(); - primaryVertex.setPos({collision.posX(), collision.posY(), collision.posZ()}); - primaryVertex.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()); - } else { - primaryVertex.setPos({mVtx->getX(), mVtx->getY(), mVtx->getZ()}); - } + trackTofInfo pTof, nTof, bTof; // information storage auto pTra = cascade.posTrack_as(); auto nTra = cascade.negTrack_as(); auto bTra = cascade.bachelor_as(); - processCascadeCandidate(primaryVertex, cascade, pTra, nTra, bTra, 0); + + pTof.collisionId = pTra.collisionId(); + pTof.hasITS = pTra.hasITS(); + pTof.hasTPC = pTra.hasTPC(); + pTof.hasTOF = pTra.hasTOF(); + pTof.tofExpMom = pTra.tofExpMom(); + pTof.tofEvTime = collisionEventTime[cascade.collisionId()]; + pTof.tofSignal = pTra.tofSignal(); + pTof.length = pTra.length(); + pTof.tpcNSigmaPi = pTra.tpcNSigmaPi(); + pTof.tpcNSigmaPr = pTra.tpcNSigmaPr(); + + nTof.collisionId = nTra.collisionId(); + nTof.hasITS = nTra.hasITS(); + nTof.hasTPC = nTra.hasTPC(); + nTof.hasTOF = nTra.hasTOF(); + nTof.tofExpMom = nTra.tofExpMom(); + nTof.tofEvTime = collisionEventTime[cascade.collisionId()]; + nTof.tofSignal = nTra.tofSignal(); + nTof.length = nTra.length(); + nTof.tpcNSigmaPi = nTra.tpcNSigmaPi(); + nTof.tpcNSigmaPr = nTra.tpcNSigmaPr(); + + bTof.collisionId = bTra.collisionId(); + bTof.hasITS = bTra.hasITS(); + bTof.hasTPC = bTra.hasTPC(); + bTof.hasTOF = bTra.hasTOF(); + bTof.tofExpMom = bTra.tofExpMom(); + bTof.tofEvTime = collisionEventTime[cascade.collisionId()]; + bTof.tofSignal = bTra.tofSignal(); + bTof.length = bTra.length(); + bTof.tpcNSigmaPi = bTra.tpcNSigmaPi(); + bTof.tpcNSigmaKa = bTra.tpcNSigmaKa(); + + cascTofInfo casctof = calculateTofInfoCascade(collisions, cascade.collisionId(), cascade, pTof, nTof, bTof); + + if (doNSigmas) { + casctofnsigmas( + casctof.nSigmaXiLaPi, casctof.nSigmaXiLaPr, casctof.nSigmaXiPi, + casctof.nSigmaOmLaPi, casctof.nSigmaOmLaPr, casctof.nSigmaOmKa); + } } } } - void processDerivedData(soa::Join const& collisions, V0DerivedDatas const& V0s, CascDerivedDatas const& cascades, dauTracks const&) + void processDerivedData(soa::Join const& collisions, V0DerivedDatas const& V0s, CascDerivedDatas const& cascades, dauTracks const& dauTrackTable, aod::DauTrackTOFPIDs const& dauTrackTOFPIDs) { + // auto-determine if current or old generation of dauTrackTOFPIDs + if (dauTrackTOFPIDs.size() == 0) { + return; + } + auto firstTOFPID = dauTrackTOFPIDs.rawIteratorAt(0); + bool isNewTOFFOrmat = firstTOFPID.straCollisionId() < 0 ? false : true; + // Fire up CCDB with first collision in record. If no collisions, bypass if (useCustomRunNumber || collisions.size() < 1) { initCCDB(manualRunNumber); @@ -1537,117 +1179,130 @@ struct strangenesstofpid { initCCDB(collision.runNumber()); } + // hold indices + std::vector tofIndices(dauTrackTable.size(), -1); + + if (isNewTOFFOrmat) { + // re-index + for (const auto& dauTrackTOFPID : dauTrackTOFPIDs) { + tofIndices[dauTrackTOFPID.dauTrackExtraId()] = dauTrackTOFPID.globalIndex(); + } + } else { + // they are actually joinable + std::iota(tofIndices.begin(), tofIndices.end(), 0); + } + if (calculateV0s.value) { for (const auto& V0 : V0s) { - // for storing whatever is the relevant quantity for the PV - o2::dataformats::VertexBase primaryVertex; - if (V0.has_straCollision()) { - auto const& collision = V0.straCollision_as>(); - primaryVertex.setPos({collision.posX(), collision.posY(), collision.posZ()}); - // cov: won't be used anyways, all fine - primaryVertex.setCov(1e-6, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6); - } else { - primaryVertex.setPos({mVtx->getX(), mVtx->getY(), mVtx->getZ()}); - } + trackTofInfo pTof, nTof; // information storage + auto collision = collisions.rawIteratorAt(V0.straCollisionId()); auto pTra = V0.posTrackExtra_as(); auto nTra = V0.negTrackExtra_as(); - processV0Candidate(primaryVertex, V0, pTra, nTra, 0); + + pTof.hasITS = pTra.hasITS(); + pTof.hasTPC = pTra.hasTPC(); + pTof.hasTOF = pTra.hasTOF(); + pTof.tpcNSigmaPi = pTra.tpcNSigmaPi(); + pTof.tpcNSigmaPr = pTra.tpcNSigmaPr(); + if (tofIndices[V0.posTrackExtraId()] >= 0 && collision.eventTime() > -1e+5) { + auto pTofExt = dauTrackTOFPIDs.rawIteratorAt(tofIndices[V0.posTrackExtraId()]); + pTof.collisionId = pTofExt.straCollisionId(); + pTof.tofExpMom = pTofExt.tofExpMom(); + pTof.tofEvTime = collision.eventTime(); + pTof.tofSignal = pTofExt.tofSignal(); + pTof.length = pTofExt.length(); + } + + nTof.hasITS = nTra.hasITS(); + nTof.hasTPC = nTra.hasTPC(); + nTof.hasTOF = nTra.hasTOF(); + nTof.tpcNSigmaPi = nTra.tpcNSigmaPi(); + nTof.tpcNSigmaPr = nTra.tpcNSigmaPr(); + if (tofIndices[V0.negTrackExtraId()] >= 0 && collision.eventTime() > -1e+5) { + auto nTofExt = dauTrackTOFPIDs.rawIteratorAt(tofIndices[V0.negTrackExtraId()]); + nTof.collisionId = nTofExt.straCollisionId(); + nTof.tofExpMom = nTofExt.tofExpMom(); + nTof.tofEvTime = collision.eventTime(); + nTof.tofSignal = nTofExt.tofSignal(); + nTof.length = nTofExt.length(); + } + + v0TofInfo v0tof = calculateTofInfoV0(collisions, V0.straCollisionId(), V0, pTof, nTof); + + if (doNSigmas) { + v0tofnsigmas( + v0tof.nSigmaPositiveLambdaPr, v0tof.nSigmaNegativeLambdaPi, + v0tof.nSigmaNegativeLambdaPr, v0tof.nSigmaPositiveLambdaPi, + v0tof.nSigmaPositiveK0ShortPi, v0tof.nSigmaNegativeK0ShortPi); + } } } if (calculateCascades.value) { for (const auto& cascade : cascades) { - // for storing whatever is the relevant quantity for the PV - o2::dataformats::VertexBase primaryVertex; - if (cascade.has_straCollision()) { - auto const& collision = cascade.straCollision_as>(); - primaryVertex.setPos({collision.posX(), collision.posY(), collision.posZ()}); - primaryVertex.setCov(1e-6, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6); - } else { - primaryVertex.setPos({mVtx->getX(), mVtx->getY(), mVtx->getZ()}); - } + trackTofInfo pTof, nTof, bTof; // information storage + auto collision = collisions.rawIteratorAt(cascade.straCollisionId()); auto pTra = cascade.posTrackExtra_as(); auto nTra = cascade.negTrackExtra_as(); auto bTra = cascade.bachTrackExtra_as(); - processCascadeCandidate(primaryVertex, cascade, pTra, nTra, bTra, 0); - } - } - } - - void processDerivedDataMCTest(soa::Join const& collisions, V0DerivedDatasMC const& V0s, CascDerivedDatasMC const& cascades, dauTracks const&, aod::V0MCCores const& v0mcs, aod::CascMCCores const& cascmcs) - { - // Fire up CCDB with first collision in record. If no collisions, bypass - if (useCustomRunNumber || collisions.size() < 1) { - initCCDB(manualRunNumber); - } else { - auto collision = collisions.begin(); - initCCDB(collision.runNumber()); - } - if (calculateV0s.value) { - for (const auto& V0 : V0s) { - // for storing whatever is the relevant quantity for the PV - o2::dataformats::VertexBase primaryVertex; - if (V0.has_straCollision()) { - auto const& collision = V0.straCollision_as>(); - primaryVertex.setPos({collision.posX(), collision.posY(), collision.posZ()}); - // cov: won't be used anyways, all fine - primaryVertex.setCov(1e-6, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6); - } else { - primaryVertex.setPos({mVtx->getX(), mVtx->getY(), mVtx->getZ()}); + pTof.hasITS = pTra.hasITS(); + pTof.hasTPC = pTra.hasTPC(); + pTof.hasTOF = pTra.hasTOF(); + pTof.tpcNSigmaPi = pTra.tpcNSigmaPi(); + pTof.tpcNSigmaPr = pTra.tpcNSigmaPr(); + if (tofIndices[cascade.posTrackExtraId()] >= 0 && collision.eventTime() > -1e+5) { + auto pTofExt = dauTrackTOFPIDs.rawIteratorAt(tofIndices[cascade.posTrackExtraId()]); + pTof.collisionId = pTofExt.straCollisionId(); + pTof.tofExpMom = pTofExt.tofExpMom(); + pTof.tofEvTime = collision.eventTime(); + pTof.tofSignal = pTofExt.tofSignal(); + pTof.length = pTofExt.length(); } - // check association - int v0pdg = 0; - if (V0.v0MCCoreId() > -1) { - auto v0mc = v0mcs.rawIteratorAt(V0.v0MCCoreId()); - v0pdg = v0mc.pdgCode(); - if (std::abs(v0pdg) != 3122 && v0pdg != 310) { - continue; // only associated from this point on - } + nTof.hasITS = nTra.hasITS(); + nTof.hasTPC = nTra.hasTPC(); + nTof.hasTOF = nTra.hasTOF(); + nTof.tpcNSigmaPi = nTra.tpcNSigmaPi(); + nTof.tpcNSigmaPr = nTra.tpcNSigmaPr(); + if (tofIndices[cascade.negTrackExtraId()] >= 0 && collision.eventTime() > -1e+5) { + auto nTofExt = dauTrackTOFPIDs.rawIteratorAt(tofIndices[cascade.negTrackExtraId()]); + nTof.collisionId = nTofExt.straCollisionId(); + nTof.tofExpMom = nTofExt.tofExpMom(); + nTof.tofEvTime = collision.eventTime(); + nTof.tofSignal = nTofExt.tofSignal(); + nTof.length = nTofExt.length(); } - auto pTra = V0.posTrackExtra_as(); - auto nTra = V0.negTrackExtra_as(); - processV0Candidate(primaryVertex, V0, pTra, nTra, v0pdg); - } - } - - if (calculateCascades.value) { - for (const auto& cascade : cascades) { - // for storing whatever is the relevant quantity for the PV - o2::dataformats::VertexBase primaryVertex; - if (cascade.has_straCollision()) { - auto const& collision = cascade.straCollision_as>(); - primaryVertex.setPos({collision.posX(), collision.posY(), collision.posZ()}); - primaryVertex.setCov(1e-6, 1e-6, 1e-6, 1e-6, 1e-6, 1e-6); - } else { - primaryVertex.setPos({mVtx->getX(), mVtx->getY(), mVtx->getZ()}); + bTof.hasITS = bTra.hasITS(); + bTof.hasTPC = bTra.hasTPC(); + bTof.hasTOF = bTra.hasTOF(); + bTof.tpcNSigmaPi = bTra.tpcNSigmaPi(); + bTof.tpcNSigmaKa = bTra.tpcNSigmaKa(); + if (tofIndices[cascade.bachTrackExtraId()] >= 0 && collision.eventTime() > -1e+5) { + auto bTofExt = dauTrackTOFPIDs.rawIteratorAt(tofIndices[cascade.bachTrackExtraId()]); + bTof.collisionId = bTofExt.straCollisionId(); + bTof.tofExpMom = bTofExt.tofExpMom(); + bTof.tofEvTime = collision.eventTime(); + bTof.tofSignal = bTofExt.tofSignal(); + bTof.length = bTofExt.length(); } - // check association - int cascpdg = 0; - if (cascade.cascMCCoreId() > -1) { - auto cascmc = cascmcs.rawIteratorAt(cascade.cascMCCoreId()); - cascpdg = cascmc.pdgCode(); - if (std::abs(cascpdg) != 3312 && std::abs(cascpdg) != 3334) { - continue; // only associated from this point on - } - } + cascTofInfo casctof = calculateTofInfoCascade(collisions, cascade.straCollisionId(), cascade, pTof, nTof, bTof); - auto pTra = cascade.posTrackExtra_as(); - auto nTra = cascade.negTrackExtra_as(); - auto bTra = cascade.bachTrackExtra_as(); - processCascadeCandidate(primaryVertex, cascade, pTra, nTra, bTra, cascpdg); + if (doNSigmas) { + casctofnsigmas( + casctof.nSigmaXiLaPi, casctof.nSigmaXiLaPr, casctof.nSigmaXiPi, + casctof.nSigmaOmLaPi, casctof.nSigmaOmLaPr, casctof.nSigmaOmKa); + } } } } PROCESS_SWITCH(strangenesstofpid, processStandardData, "Process standard data", false); PROCESS_SWITCH(strangenesstofpid, processDerivedData, "Process derived data", true); - PROCESS_SWITCH(strangenesstofpid, processDerivedDataMCTest, "Process derived data / MC with assoc / not for analysis", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)