diff --git a/PWGLF/DataModel/LFStrangenessPIDTables.h b/PWGLF/DataModel/LFStrangenessPIDTables.h index e689dc28141..0b7e95ce302 100644 --- a/PWGLF/DataModel/LFStrangenessPIDTables.h +++ b/PWGLF/DataModel/LFStrangenessPIDTables.h @@ -47,6 +47,10 @@ DECLARE_SOA_TABLE(DauTrackTOFPIDs, "AOD", "DAUTRACKTOFPID", // raw table (for po namespace v0data { +// define constants for NSigma operation +const float kNoTOFValue = -1e+6; +const float kEpsilon = 1e-4; + // ==== TOF INFORMATION === // lengths as stored in the AO2D for TOF calculations DECLARE_SOA_COLUMN(PosTOFLengthToPV, posTOFLengthToPV, float); //! positive track length to PV @@ -79,6 +83,45 @@ DECLARE_SOA_COLUMN(TOFNSigmaALaPi, tofNSigmaALaPi, float); //! positive DECLARE_SOA_COLUMN(TOFNSigmaK0PiPlus, tofNSigmaK0PiPlus, float); //! positive track NSigma from pion <- k0short expectation DECLARE_SOA_COLUMN(TOFNSigmaK0PiMinus, tofNSigmaK0PiMinus, float); //! negative track NSigma from pion <- k0short expectation +// dynamics based on n-sigmas with use-only-if-tof-present logic +DECLARE_SOA_DYNAMIC_COLUMN(TofLambdaCompatibility, tofLambdaCompatibility, //! compatibility with being lambda, checked only if TOF present. Argument: number of sigmas + [](float tofNSigmaLaPr, float tofNSigmaLaPi, float nsigma) -> float { + bool compatible = true; + if (std::abs(tofNSigmaLaPr - kNoTOFValue) > kEpsilon && std::abs(tofNSigmaLaPr) > nsigma) { + compatible = false; // reject only if info present and incompatible + } + if (std::abs(tofNSigmaLaPi - kNoTOFValue) > kEpsilon && std::abs(tofNSigmaLaPi) > nsigma) { + compatible = false; // reject only if info present and incompatible + } + return compatible; + }); + +// dynamics based on n-sigmas with use-only-if-tof-present logic +DECLARE_SOA_DYNAMIC_COLUMN(TofAntiLambdaCompatibility, tofAntiLambdaCompatibility, //! compatibility with being lambda, checked only if TOF present. Argument: number of sigmas + [](float tofNSigmaALaPr, float tofNSigmaALaPi, float nsigma) -> float { + bool compatible = true; + if (std::abs(tofNSigmaALaPr - kNoTOFValue) > kEpsilon && std::abs(tofNSigmaALaPr) > nsigma) { + compatible = false; // reject only if info present and incompatible + } + if (std::abs(tofNSigmaALaPi - kNoTOFValue) > kEpsilon && std::abs(tofNSigmaALaPi) > nsigma) { + compatible = false; // reject only if info present and incompatible + } + return compatible; + }); + +// dynamics based on n-sigmas with use-only-if-tof-present logic +DECLARE_SOA_DYNAMIC_COLUMN(TofK0ShortCompatibility, tofK0ShortCompatibility, //! compatibility with being lambda, checked only if TOF present. Argument: number of sigmas + [](float tofNSigmaK0PiPlus, float tofNSigmaK0PiMinus, float nsigma) -> float { + bool compatible = true; + if (std::abs(tofNSigmaK0PiPlus - kNoTOFValue) > kEpsilon && std::abs(tofNSigmaK0PiPlus) > nsigma) { + compatible = false; // reject only if info present and incompatible + } + if (std::abs(tofNSigmaK0PiMinus - kNoTOFValue) > kEpsilon && std::abs(tofNSigmaK0PiMinus) > nsigma) { + compatible = false; // reject only if info present and incompatible + } + return compatible; + }); + // beta values DECLARE_SOA_COLUMN(TofBetaLambda, tofBetaLambda, float); //! beta value with Lambda hypothesis DECLARE_SOA_COLUMN(TofBetaAntiLambda, tofBetaAntiLambda, float); //! beta value with AntiLambda hypothesis @@ -120,10 +163,17 @@ DECLARE_SOA_TABLE(V0TOFBetas, "AOD", "V0TOFBETA", // processed info table (for a DECLARE_SOA_TABLE(V0TOFNSigmas, "AOD", "V0TOFNSIGMA", // processed NSigma table (for analysis) v0data::TOFNSigmaLaPr, v0data::TOFNSigmaLaPi, v0data::TOFNSigmaALaPr, v0data::TOFNSigmaALaPi, - v0data::TOFNSigmaK0PiPlus, v0data::TOFNSigmaK0PiMinus); + v0data::TOFNSigmaK0PiPlus, v0data::TOFNSigmaK0PiMinus, + v0data::TofLambdaCompatibility, + v0data::TofAntiLambdaCompatibility, + v0data::TofK0ShortCompatibility); namespace cascdata { +// define constants for NSigma operation +const float kNoTOFValue = -1e+6; +const float kEpsilon = 1e-4; + // lengths as stored in the AO2D for TOF calculations DECLARE_SOA_COLUMN(PosTOFLengthToPV, posTOFLengthToPV, float); //! positive track length DECLARE_SOA_COLUMN(NegTOFLengthToPV, negTOFLengthToPV, float); //! negative track length @@ -154,6 +204,38 @@ DECLARE_SOA_COLUMN(TOFNSigmaXiPi, tofNSigmaXiPi, float); //! bachelor track DECLARE_SOA_COLUMN(TOFNSigmaOmLaPi, tofNSigmaOmLaPi, float); //! meson track NSigma from pion <- lambda <- om expectation DECLARE_SOA_COLUMN(TOFNSigmaOmLaPr, tofNSigmaOmLaPr, float); //! baryon track NSigma from proton <- lambda <- om expectation DECLARE_SOA_COLUMN(TOFNSigmaOmKa, tofNSigmaOmKa, float); //! bachelor track NSigma from kaon <- om expectation + +// dynamics based on n-sigmas with use-only-if-tof-present logic +DECLARE_SOA_DYNAMIC_COLUMN(TofXiCompatibility, tofXiCompatibility, //! compatibility with being lambda, checked only if TOF present. Argument: number of sigmas + [](float tofNSigmaXiLaPr, float tofNSigmaXiLaPi, float tofNSigmaXiPi, float nsigma) -> float { + bool compatible = true; + if (std::abs(tofNSigmaXiLaPr - kNoTOFValue) > kEpsilon && std::abs(tofNSigmaXiLaPr) > nsigma) { + compatible = false; // reject only if info present and incompatible + } + if (std::abs(tofNSigmaXiLaPi - kNoTOFValue) > kEpsilon && std::abs(tofNSigmaXiLaPi) > nsigma) { + compatible = false; // reject only if info present and incompatible + } + if (std::abs(tofNSigmaXiPi - kNoTOFValue) > kEpsilon && std::abs(tofNSigmaXiPi) > nsigma) { + compatible = false; // reject only if info present and incompatible + } + return compatible; + }); + +DECLARE_SOA_DYNAMIC_COLUMN(TofOmegaCompatibility, tofOmegaCompatibility, //! compatibility with being lambda, checked only if TOF present. Argument: number of sigmas + [](float tofNSigmaOmLaPr, float tofNSigmaOmLaPi, float tofNSigmaOmKa, float nsigma) -> float { + bool compatible = true; + if (std::abs(tofNSigmaOmLaPr - kNoTOFValue) > kEpsilon && std::abs(tofNSigmaOmLaPr) > nsigma) { + compatible = false; // reject only if info present and incompatible + } + if (std::abs(tofNSigmaOmLaPi - kNoTOFValue) > kEpsilon && std::abs(tofNSigmaOmLaPi) > nsigma) { + compatible = false; // reject only if info present and incompatible + } + if (std::abs(tofNSigmaOmKa - kNoTOFValue) > kEpsilon && std::abs(tofNSigmaOmKa) > nsigma) { + compatible = false; // reject only if info present and incompatible + } + return compatible; + }); + } // namespace cascdata // /-|-\-|-/-|-\-|-/-|-\-|-/-|-\-|-/-|-\-|-/-|-\-|-/-|-\-|-/-|-\-|-/-|-\-|-/-|-\-|-/ @@ -173,7 +255,9 @@ DECLARE_SOA_TABLE(CascTOFPIDs, "AOD", "CASCTOFPID", // processed information for cascdata::BachTOFDeltaTOmKa); DECLARE_SOA_TABLE(CascTOFNSigmas, "AOD", "CascTOFNSigmas", // Nsigmas for cascades cascdata::TOFNSigmaXiLaPi, cascdata::TOFNSigmaXiLaPr, cascdata::TOFNSigmaXiPi, - cascdata::TOFNSigmaOmLaPi, cascdata::TOFNSigmaOmLaPr, cascdata::TOFNSigmaOmKa); + cascdata::TOFNSigmaOmLaPi, cascdata::TOFNSigmaOmLaPr, cascdata::TOFNSigmaOmKa, + cascdata::TofXiCompatibility, + cascdata::TofOmegaCompatibility); } // namespace o2::aod #endif // PWGLF_DATAMODEL_LFSTRANGENESSPIDTABLES_H_ diff --git a/PWGLF/TableProducer/Strangeness/CMakeLists.txt b/PWGLF/TableProducer/Strangeness/CMakeLists.txt index 9e4a3fa04f8..b9ae33532f2 100644 --- a/PWGLF/TableProducer/Strangeness/CMakeLists.txt +++ b/PWGLF/TableProducer/Strangeness/CMakeLists.txt @@ -36,8 +36,8 @@ o2physics_add_dpl_workflow(cascademcfinder PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(cascadepid - SOURCES cascadepid.cxx +o2physics_add_dpl_workflow(strangenesstofpid + SOURCES strangenesstofpid.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DetectorsBase COMPONENT_NAME Analysis) @@ -81,11 +81,6 @@ o2physics_add_dpl_workflow(lambdakzeromcfinder PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(lambdakzeropid - SOURCES lambdakzeropid.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DetectorsBase - COMPONENT_NAME Analysis) - o2physics_add_dpl_workflow(lambdakzerospawner SOURCES lambdakzerospawner.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore diff --git a/PWGLF/TableProducer/Strangeness/cascadepid.cxx b/PWGLF/TableProducer/Strangeness/cascadepid.cxx deleted file mode 100644 index 3c7c82dacdd..00000000000 --- a/PWGLF/TableProducer/Strangeness/cascadepid.cxx +++ /dev/null @@ -1,632 +0,0 @@ -// 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. -// -// *+-+*+-+*+-+*+-+*+-+*+-+* -// Cascade PID tables -// *+-+*+-+*+-+*+-+*+-+*+-+* -// -/// \author Nicolò Jacazio -/// \author David Dobrigkeit Chinellato -/// \since 22/11/2023 -/// \brief Table producer for Casc daughter PID info -// -// This task produces daughter PID information for strange daughters -// taking into account the (candidate-by-candidate) time spent as a heavier -// (strange, weakly-decaying) particle. This task is meant to be a test, as -// it hasn't been fully tested yet! Use at your own peril for now :-) - -#include -#include -#include -#include -#include -#include - -#include "Framework/runDataProcessing.h" -#include "Framework/RunningWorkflowInfo.h" -#include "Framework/AnalysisTask.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/ASoAHelpers.h" -#include "DCAFitter/DCAFitterN.h" -#include "ReconstructionDataFormats/Track.h" -#include "Common/Core/RecoDecay.h" -#include "Common/Core/trackUtilities.h" -#include "PWGLF/DataModel/LFStrangenessTables.h" -#include "PWGLF/DataModel/LFStrangenessPIDTables.h" -#include "PWGLF/DataModel/LFParticleIdentification.h" -#include "Common/Core/TrackSelection.h" -#include "Common/DataModel/TrackSelectionTables.h" -#include "DetectorsBase/Propagator.h" -#include "DetectorsBase/GeometryManager.h" -#include "DataFormatsParameters/GRPObject.h" -#include "DataFormatsParameters/GRPMagField.h" -#include "CCDB/BasicCCDBManager.h" -#include "CommonConstants/PhysicsConstants.h" -#include "Common/TableProducer/PID/pidTOFBase.h" -#include "Common/DataModel/PIDResponse.h" - -using namespace o2; -using namespace o2::framework; -using namespace o2::framework::expressions; -using std::array; - -// For original data loops -using CascOriginalDatas = soa::Join; -using TracksWithAllExtras = soa::Join; - -// Cores with references and TOF pid -using dauTracks = soa::Join; -using CascDerivedDatas = soa::Join; - -struct cascadepid { - // TOF pid for strangeness (recalculated with topology) - Produces casctofpids; // table with base info - Produces casctofnsigmas; // table with Nsigmas - - Service ccdb; - - // For manual sliceBy - Preslice perCollisionOriginal = o2::aod::cascdata::collisionId; - ; - Preslice perCollisionDerived = o2::aod::cascdata::straCollisionId; - - HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - - // Operation and minimisation criteria - Configurable d_bz_input{"d_bz", -999, "bz field, -999 is automatic"}; - Configurable tofPosition{"tofPosition", 377.934f, "TOF effective (inscribed) radius"}; - Configurable doQA{"doQA", true, "create QA histos"}; - Configurable qaV0DCADau{"qaV0DCADau", 0.5, "DCA daughters (cm) for QA plots"}; - Configurable qaCascDCADau{"qaCascDCADau", 0.5, "DCA daughters (cm) for QA plots"}; - Configurable qaV0CosPA{"qaV0CosPA", 0.995, "CosPA for QA plots"}; - Configurable qaCascCosPA{"qaCascCosPA", 0.995, "CosPA for QA plots"}; - Configurable qaMassWindow{"qaMassWindow", 0.005, "Mass window around expected (in GeV/c2) for QA plots"}; - Configurable qaTPCNSigma{"qaTPCNSigma", 5, "TPC N-sigma to apply for qa plots"}; - Configurable doNSigmas{"doNSigmas", false, "calculate TOF N-sigma"}; - Configurable doQANSigma{"doQANSigma", false, "create QA of Nsigma histos"}; - - // CCDB options - Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; - Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; - Configurable lutPath{"lutPath", "GLO/Param/MatLUT", "Path of the Lut parametrization"}; - Configurable geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; - Configurable nSigmaPath{"nSigmaPath", "Users/d/ddobrigk/stratof", "Path of information for n-sigma calculation"}; - - // manual - Configurable useCustomRunNumber{"useCustomRunNumber", false, "Use custom timestamp"}; - Configurable manualRunNumber{"manualRunNumber", 544122, "manual run number if no collisions saved"}; - - ConfigurableAxis axisEta{"axisEta", {20, -1.0f, +1.0f}, "#eta"}; - ConfigurableAxis axisDeltaTime{"axisDeltaTime", {2000, -1000.0f, +1000.0f}, "delta-time (ps)"}; - ConfigurableAxis axisNSigma{"axisNSigma", {200, -10.0f, +10.0f}, "N(#sigma)"}; - ConfigurableAxis axisTime{"axisTime", {200, 0.0f, +20000.0f}, "T (ps)"}; - ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "p_{T} (GeV/c)"}; - - bool nSigmaCalibLoaded; - TList* nSigmaCalibObjects; - TH1 *hMeanPosXiPi, *hSigmaPosXiPi; - TH1 *hMeanPosXiPr, *hSigmaPosXiPr; - TH1 *hMeanNegXiPi, *hSigmaNegXiPi; - TH1 *hMeanNegXiPr, *hSigmaNegXiPr; - TH1 *hMeanBachXiPi, *hSigmaBachXiPi; - TH1 *hMeanPosOmPi, *hSigmaPosOmPi; - TH1 *hMeanPosOmPr, *hSigmaPosOmPr; - TH1 *hMeanNegOmPi, *hSigmaNegOmPi; - TH1 *hMeanNegOmPr, *hSigmaNegOmPr; - TH1 *hMeanBachOmKa, *hSigmaBachOmKa; - - int mRunNumber; - float d_bz; - float maxSnp; // max sine phi for propagation - float maxStep; // max step size (cm) for propagation - - /// function to calculate track length of this track up to a certain segment of a detector - /// to be used internally in another funcrtion that calculates length until it finds the proper one - /// warning: this could be optimised further for speed - /// \param track the input track - /// \param x1 x of the first point of the detector segment - /// \param y1 y of the first point of the detector segment - /// \param x2 x of the first point of the detector segment - /// \param y2 y of the first point of the detector segment - /// \param magneticField the magnetic field to use when propagating - float trackLengthToSegment(o2::track::TrackPar track, float x1, float y1, float x2, float y2, float magneticField) - { - // don't make use of the track parametrization - float length = -104; - - // causality protection - std::array mom; - track.getPxPyPzGlo(mom); - // get start point - std::array startPoint; - track.getXYZGlo(startPoint); - - // better replaced with scalar momentum check later - // if (((x1 + x2) * mom[0] + (y1 + y2) * mom[1]) < 0.0f) - // return -101; - - // get circle X, Y please - o2::math_utils::CircleXYf_t trcCircle; - float sna, csa; - track.getCircleParams(magneticField, trcCircle, sna, csa); - - // Calculate necessary inner product - float segmentModulus = std::hypot(x2 - x1, y2 - y1); - float alongSegment = ((trcCircle.xC - x1) * (x2 - x1) + (trcCircle.yC - y1) * (y2 - y1)) / segmentModulus; - - // find point of closest approach between segment and circle center - float pcaX = (x2 - x1) * alongSegment / segmentModulus + x1; - float pcaY = (y2 - y1) * alongSegment / segmentModulus + y1; - - float centerDistToPC = std::hypot(pcaX - trcCircle.xC, pcaY - trcCircle.yC); - - // distance pca-to-intercept in multiples of segment modulus (for convenience) - if (centerDistToPC > trcCircle.rC) - return -103; - - float pcaToIntercept = TMath::Sqrt(TMath::Abs(trcCircle.rC * trcCircle.rC - centerDistToPC * centerDistToPC)); - - float interceptX1 = pcaX + (x2 - x1) / segmentModulus * pcaToIntercept; - float interceptY1 = pcaY + (y2 - y1) / segmentModulus * pcaToIntercept; - float interceptX2 = pcaX - (x2 - x1) / segmentModulus * pcaToIntercept; - float interceptY2 = pcaY - (y2 - y1) / segmentModulus * pcaToIntercept; - - float scalarCheck1 = ((x2 - x1) * (interceptX1 - x1) + (y2 - y1) * (interceptY1 - y1)) / segmentModulus; - float scalarCheck2 = ((x2 - x1) * (interceptX2 - x1) + (y2 - y1) * (interceptY2 - y1)) / segmentModulus; - - float cosAngle1 = -1000, sinAngle1 = -1000, modulus1 = -1000; - float cosAngle2 = -1000, sinAngle2 = -1000, modulus2 = -1000; - float length1 = -1000, length2 = -1000; - - modulus1 = std::hypot(interceptX1 - trcCircle.xC, interceptY1 - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC); - cosAngle1 = (interceptX1 - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (interceptY1 - trcCircle.yC) * (startPoint[1] - trcCircle.yC); - sinAngle1 = (interceptX1 - trcCircle.xC) * (startPoint[1] - trcCircle.yC) - (interceptY1 - trcCircle.yC) * (startPoint[0] - trcCircle.xC); - cosAngle1 /= modulus1; - sinAngle1 /= modulus1; - length1 = trcCircle.rC * TMath::ACos(cosAngle1); - length1 *= sqrt(1.0f + track.getTgl() * track.getTgl()); - - modulus2 = std::hypot(interceptX2 - trcCircle.xC, interceptY2 - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC); - cosAngle2 = (interceptX2 - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (interceptY2 - trcCircle.yC) * (startPoint[1] - trcCircle.yC); - sinAngle2 = (interceptX2 - trcCircle.xC) * (startPoint[1] - trcCircle.yC) - (interceptY2 - trcCircle.yC) * (startPoint[0] - trcCircle.xC); - cosAngle2 /= modulus2; - sinAngle2 /= modulus2; - length2 = trcCircle.rC * TMath::ACos(cosAngle2); - length2 *= sqrt(1.0f + track.getTgl() * track.getTgl()); - - // rotate transverse momentum vector such that it is at intercepts - float angle1 = TMath::ACos(cosAngle1); - if (sinAngle1 < 0) - angle1 *= -1.0f; - float px1 = +TMath::Cos(angle1) * mom[0] + TMath::Sin(angle1) * mom[1]; - float py1 = -TMath::Sin(angle1) * mom[0] + TMath::Cos(angle1) * mom[1]; - - float angle2 = TMath::ACos(cosAngle2); - if (sinAngle2 < 0) - angle2 *= -1.0f; - float px2 = +TMath::Cos(angle2) * mom[0] + TMath::Sin(angle2) * mom[1]; - float py2 = -TMath::Sin(angle2) * mom[0] + TMath::Cos(angle2) * mom[1]; - - float midSegX = 0.5f * (x2 + x1); - float midSegY = 0.5f * (y2 + y1); - - float scalarMomentumCheck1 = px1 * midSegX + py1 * midSegY; - float scalarMomentumCheck2 = px2 * midSegX + py2 * midSegY; - - if (scalarCheck1 > 0.0f && scalarCheck1 < segmentModulus && scalarMomentumCheck1 > 0.0f) { - length = length1; - // X = interceptX1; Y = interceptY1; - } - if (scalarCheck2 > 0.0f && scalarCheck2 < segmentModulus && scalarMomentumCheck2 > 0.0f) { - length = length2; - // X = interceptX2; Y = interceptY2; - } - return length; - } - - /// function to calculate track length of this track up to a certain segmented detector - /// \param track the input track - /// \param magneticField the magnetic field to use when propagating - float findInterceptLength(o2::track::TrackPar track, float magneticField) - { - float length = 1e+6; - for (int iSeg = 0; iSeg < 18; iSeg++) { - // Detector segmentation loop - float segmentAngle = 20.0f / 180.0f * TMath::Pi(); - float theta = static_cast(iSeg) * 20.0f / 180.0f * TMath::Pi(); - float halfWidth = tofPosition * TMath::Tan(0.5f * segmentAngle); - float x1 = TMath::Cos(theta) * (-halfWidth) + TMath::Sin(theta) * tofPosition; - float y1 = -TMath::Sin(theta) * (-halfWidth) + TMath::Cos(theta) * tofPosition; - float x2 = TMath::Cos(theta) * (+halfWidth) + TMath::Sin(theta) * tofPosition; - float y2 = -TMath::Sin(theta) * (+halfWidth) + TMath::Cos(theta) * tofPosition; - float thisLength = trackLengthToSegment(track, x1, y1, x2, y2, magneticField); - if (thisLength < length && thisLength > 0) - length = thisLength; - } - if (length > 1e+5) - length = -100; // force negative to avoid misunderstandings - return length; - } - - void init(InitContext&) - { - mRunNumber = 0; - d_bz = 0; - maxSnp = 0.85f; // could be changed later - maxStep = 2.00f; // could be changed later - - ccdb->setURL(ccdburl); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - ccdb->setFatalWhenNull(false); - - nSigmaCalibLoaded = false; - nSigmaCalibObjects = nullptr; - - // measured vs expected total time QA - if (doQA) { - // standard deltaTime values - histos.add("hArcDebug", "hArcDebug", kTH2F, {axisPt, {500, -5.0f, 10.0f}}); - histos.add("h2dposDeltaTimeAsXiPi", "h2dposDeltaTimeAsXiPi", {HistType::kTH3F, {axisPt, axisEta, axisDeltaTime}}); - histos.add("h2dposDeltaTimeAsXiPr", "h2dposDeltaTimeAsXiPr", {HistType::kTH3F, {axisPt, axisEta, axisDeltaTime}}); - histos.add("h2dnegDeltaTimeAsXiPi", "h2dnegDeltaTimeAsXiPi", {HistType::kTH3F, {axisPt, axisEta, axisDeltaTime}}); - histos.add("h2dnegDeltaTimeAsXiPr", "h2dnegDeltaTimeAsXiPr", {HistType::kTH3F, {axisPt, axisEta, axisDeltaTime}}); - histos.add("h2dbachDeltaTimeAsXiPi", "h2dbachDeltaTimeAsXiPi", {HistType::kTH3F, {axisPt, axisEta, axisDeltaTime}}); - - histos.add("h2dposDeltaTimeAsOmPi", "h2dposDeltaTimeAsOmPi", {HistType::kTH3F, {axisPt, axisEta, axisDeltaTime}}); - histos.add("h2dposDeltaTimeAsOmPr", "h2dposDeltaTimeAsOmPr", {HistType::kTH3F, {axisPt, axisEta, axisDeltaTime}}); - histos.add("h2dnegDeltaTimeAsOmPi", "h2dnegDeltaTimeAsOmPi", {HistType::kTH3F, {axisPt, axisEta, axisDeltaTime}}); - histos.add("h2dnegDeltaTimeAsOmPr", "h2dnegDeltaTimeAsOmPr", {HistType::kTH3F, {axisPt, axisEta, axisDeltaTime}}); - histos.add("h2dbachDeltaTimeAsOmKa", "h2dbachDeltaTimeAsOmKa", {HistType::kTH3F, {axisPt, axisEta, axisDeltaTime}}); - } - - if (doQANSigma) { - // standard NSigma values - histos.add("h2dNSigmaXiLaPi", "h2dNSigmaXiLaPi", {HistType::kTH2F, {axisPt, axisNSigma}}); - histos.add("h2dNSigmaXiLaPr", "h2dNSigmaXiLaPr", {HistType::kTH2F, {axisPt, axisNSigma}}); - histos.add("h2dNSigmaXiPi", "h2dNSigmaXiPi", {HistType::kTH2F, {axisPt, axisNSigma}}); - histos.add("h2dNSigmaOmLaPi", "h2dNSigmaOmLaPi", {HistType::kTH2F, {axisPt, axisNSigma}}); - histos.add("h2dNSigmaOmLaPr", "h2dNSigmaOmLaPr", {HistType::kTH2F, {axisPt, axisNSigma}}); - histos.add("h2dNSigmaOmKa", "h2dNSigmaOmKa", {HistType::kTH2F, {axisPt, axisNSigma}}); - } - } - - void initCCDB(int runNumber) - { - if (mRunNumber == runNumber) { - return; - } - - // In case override, don't proceed, please - no CCDB access required - if (d_bz_input > -990) { - d_bz = d_bz_input; - o2::parameters::GRPMagField grpmag; - if (fabs(d_bz) > 1e-5) { - grpmag.setL3Current(30000.f / (d_bz / 5.0f)); - } - o2::base::Propagator::initFieldFromGRP(&grpmag); - mRunNumber = runNumber; - return; - } - - o2::parameters::GRPObject* grpo = ccdb->getForRun(grpPath, runNumber); - o2::parameters::GRPMagField* grpmag = 0x0; - if (grpo) { - o2::base::Propagator::initFieldFromGRP(grpo); - // Fetch magnetic field from ccdb for current collision - d_bz = grpo->getNominalL3Field(); - LOG(info) << "Retrieved GRP for run " << runNumber << " with magnetic field of " << d_bz << " kZG"; - } else { - grpmag = ccdb->getForRun(grpmagPath, runNumber); - if (!grpmag) { - LOG(fatal) << "Got nullptr from CCDB for path " << grpmagPath << " of object GRPMagField and " << grpPath << " of object GRPObject for run " << runNumber; - } - o2::base::Propagator::initFieldFromGRP(grpmag); - // Fetch magnetic field from ccdb for current collision - d_bz = std::lround(5.f * grpmag->getL3Current() / 30000.f); - LOG(info) << "Retrieved GRP for run " << runNumber << " with magnetic field of " << d_bz << " kZG"; - } - - // if TOF Nsigma desired - if (doNSigmas) { - nSigmaCalibObjects = ccdb->getForRun(nSigmaPath, runNumber); - if (nSigmaCalibObjects) { - LOGF(info, "loaded TList with this many objects: %i", nSigmaCalibObjects->GetEntries()); - - hMeanPosXiPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanPosXiPi")); - hMeanPosXiPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanPosXiPr")); - hMeanNegXiPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanNegXiPi")); - hMeanNegXiPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanNegXiPr")); - hMeanBachXiPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanBachXiPi")); - hMeanPosOmPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanPosOmPi")); - hMeanPosOmPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanPosOmPr")); - hMeanNegOmPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanNegOmPi")); - hMeanNegOmPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanNegOmPr")); - hMeanBachOmKa = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanBachOmKa")); - - hSigmaPosXiPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaPosXiPi")); - hSigmaPosXiPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaPosXiPr")); - hSigmaNegXiPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaNegXiPi")); - hSigmaNegXiPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaNegXiPr")); - hSigmaBachXiPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaBachXiPi")); - hSigmaPosOmPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaPosOmPi")); - hSigmaPosOmPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaPosOmPr")); - hSigmaNegOmPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaNegOmPi")); - hSigmaNegOmPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaNegOmPr")); - hSigmaBachOmKa = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaBachOmKa")); - - if (!hMeanPosXiPi || !hMeanPosXiPr || !hMeanNegXiPi || !hMeanNegXiPr || !hMeanBachXiPi) - LOG(info) << "Problems finding xi mean histograms!"; - if (!hMeanPosOmPi || !hMeanPosOmPr || !hMeanNegOmPi || !hMeanNegOmPr || !hMeanBachOmKa) - LOG(info) << "Problems finding omega sigma histograms!"; - if (!hSigmaPosXiPi || !hSigmaPosXiPr || !hSigmaNegXiPi || !hSigmaNegXiPr || !hSigmaBachXiPi) - LOG(info) << "Problems finding xi sigma histograms!"; - if (!hSigmaPosOmPi || !hSigmaPosOmPr || !hSigmaNegOmPi || !hSigmaNegOmPr || !hSigmaBachOmKa) - LOG(info) << "Problems finding omega sigma histograms!"; - } - } - mRunNumber = runNumber; - } - - float velocity(float lMomentum, float lMass) - { - // Momentum p and mass m -> returns speed in centimeters per picosecond - // Useful for TOF calculations - float lA = (lMomentum / lMass) * (lMomentum / lMass); - return 0.0299792458 * TMath::Sqrt(lA / (1 + lA)); - } - - template - void processCascadeCandidate(TCollision const& collision, TCascade const& cascade, TTrack const& pTra, TTrack const& nTra, TTrack const& bTra) - { - // initialize from positions and momenta as needed - o2::track::TrackPar posTrack = o2::track::TrackPar({cascade.xlambda(), cascade.ylambda(), cascade.zlambda()}, {cascade.pxpos(), cascade.pypos(), cascade.pzpos()}, +1); - o2::track::TrackPar negTrack = o2::track::TrackPar({cascade.xlambda(), cascade.ylambda(), cascade.zlambda()}, {cascade.pxneg(), cascade.pyneg(), cascade.pzneg()}, -1); - o2::track::TrackPar bachTrack = o2::track::TrackPar({cascade.x(), cascade.y(), cascade.z()}, {cascade.pxbach(), cascade.pybach(), cascade.pzbach()}, cascade.sign()); - o2::track::TrackPar cascTrack = o2::track::TrackPar({cascade.x(), cascade.y(), cascade.z()}, {cascade.px(), cascade.py(), cascade.pz()}, cascade.sign()); - - // start calculation: calculate velocities - 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 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 daughter length to TOF intercept - 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? - - // calculate mother lengths - float lengthV0 = std::hypot(cascade.xlambda() - cascade.x(), cascade.ylambda() - cascade.y(), cascade.zlambda() - cascade.z()); - float lengthCascade = -1e+6; - 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, d3d = 0.0f; - float linearToPV = std::hypot(cascade.x() - collision.posX(), cascade.y() - collision.posY(), cascade.z() - collision.posZ()); - if (successPropag) { - std::array cascCloseToPVPosition; - cascTrack.getXYZGlo(cascCloseToPVPosition); - o2::math_utils::CircleXYf_t trcCircleCascade; - float sna, csa; - cascTrack.getCircleParams(d_bz, trcCircleCascade, sna, csa); - - // calculate 2D distance between two points - d = std::hypot(cascade.x() - cascCloseToPVPosition[0], cascade.y() - cascCloseToPVPosition[1]); - d3d = std::hypot(cascade.x() - cascCloseToPVPosition[0], cascade.y() - cascCloseToPVPosition[1], cascade.z() - cascCloseToPVPosition[2]); // cross-check variable - float sinThetaOverTwo = d / (2.0f * trcCircleCascade.rC); - lengthCascade = 2.0f * trcCircleCascade.rC * TMath::ASin(sinThetaOverTwo); - lengthCascade *= sqrt(1.0f + cascTrack.getTgl() * cascTrack.getTgl()); - } - - if (!successPropag) { - lengthCascade = linearToPV; // if propagation failed, use linear estimate (optional: actually do not define?) - } - - // lambda, xi and omega flight time is always defined - float lambdaFlight = lengthV0 / velocityLa; - float xiFlight = lengthCascade / velocityXi; - float omFlight = lengthCascade / velocityOm; - float posFlightPi = lengthPositive / velocityPositivePi; - float posFlightPr = lengthPositive / velocityPositivePr; - float negFlightPi = lengthNegative / velocityNegativePi; - float negFlightPr = lengthNegative / velocityNegativePr; - float bachFlightPi = lengthBachelor / velocityBachelorPi; - float bachFlightKa = lengthBachelor / velocityBachelorKa; - - // initialize delta-times (actual PID variables) - float posDeltaTimeAsXiPi = -1e+6, posDeltaTimeAsXiPr = -1e+6; - float negDeltaTimeAsXiPi = -1e+6, negDeltaTimeAsXiPr = -1e+6; - float bachDeltaTimeAsXiPi = -1e+6; - float posDeltaTimeAsOmPi = -1e+6, posDeltaTimeAsOmPr = -1e+6; - float negDeltaTimeAsOmPi = -1e+6, negDeltaTimeAsOmPr = -1e+6; - float bachDeltaTimeAsOmKa = -1e+6; - - if (pTra.hasTOF()) { - 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()) { - 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()) { - 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 = -1e+6; - float nSigmaXiLaPi = -1e+6; - float nSigmaXiPi = -1e+6; - float nSigmaOmLaPr = -1e+6; - float nSigmaOmLaPi = -1e+6; - float nSigmaOmKa = -1e+6; - - // go for Nsigma values if requested - if (doNSigmas) { - // Xi hypothesis ________________________ - if (cascade.sign() < 0) { // XiMinus - if (posDeltaTimeAsXiPr > -1e+5) // proton from Lambda from XiMinus has signal - nSigmaXiLaPr = (posDeltaTimeAsXiPr - hMeanPosXiPr->Interpolate(cascade.pt())) / hSigmaPosXiPr->Interpolate(cascade.pt()); - if (negDeltaTimeAsXiPi > -1e+5) // pion from Lambda from XiMinus has signal - nSigmaXiLaPi = (negDeltaTimeAsXiPi - hMeanNegXiPi->Interpolate(cascade.pt())) / hSigmaNegXiPi->Interpolate(cascade.pt()); - if (bachDeltaTimeAsXiPi > -1e+5) // pion from XiMinus has signal - nSigmaXiPi = (bachDeltaTimeAsXiPi - hMeanBachXiPi->Interpolate(cascade.pt())) / hSigmaBachXiPi->Interpolate(cascade.pt()); - if (posDeltaTimeAsOmPr > -1e+5) // proton from Lambda from OmegaMinus has signal - nSigmaOmLaPr = (posDeltaTimeAsOmPr - hMeanPosOmPr->Interpolate(cascade.pt())) / hSigmaPosOmPr->Interpolate(cascade.pt()); - if (negDeltaTimeAsOmPi > -1e+5) // pion from Lambda from OmegaMinus has signal - nSigmaOmLaPi = (negDeltaTimeAsOmPi - hMeanNegOmPi->Interpolate(cascade.pt())) / hSigmaNegOmPi->Interpolate(cascade.pt()); - if (bachDeltaTimeAsOmKa > -1e+5) // kaon from OmegaMinus has signal - nSigmaOmKa = (bachDeltaTimeAsOmKa - hMeanBachOmKa->Interpolate(cascade.pt())) / hSigmaBachOmKa->Interpolate(cascade.pt()); - } else { - if (posDeltaTimeAsXiPi > -1e+5) // proton from Lambda from XiMinus has signal - nSigmaXiLaPi = (posDeltaTimeAsXiPi - hMeanPosXiPi->Interpolate(cascade.pt())) / hSigmaPosXiPi->Interpolate(cascade.pt()); - if (negDeltaTimeAsXiPr > -1e+5) // pion from Lambda from XiMinus has signal - nSigmaXiLaPr = (negDeltaTimeAsXiPr - hMeanNegXiPr->Interpolate(cascade.pt())) / hSigmaNegXiPr->Interpolate(cascade.pt()); - if (bachDeltaTimeAsXiPi > -1e+5) // pion from XiMinus has signal - nSigmaXiPi = (bachDeltaTimeAsXiPi - hMeanBachXiPi->Interpolate(cascade.pt())) / hSigmaBachXiPi->Interpolate(cascade.pt()); - if (posDeltaTimeAsOmPi > -1e+5) // proton from Lambda from OmegaMinus has signal - nSigmaOmLaPi = (posDeltaTimeAsOmPi - hMeanPosOmPi->Interpolate(cascade.pt())) / hSigmaPosOmPi->Interpolate(cascade.pt()); - if (negDeltaTimeAsOmPr > -1e+5) // pion from Lambda from OmegaMinus has signal - nSigmaOmLaPr = (negDeltaTimeAsOmPr - hMeanNegOmPr->Interpolate(cascade.pt())) / hSigmaNegOmPr->Interpolate(cascade.pt()); - if (bachDeltaTimeAsOmKa > -1e+5) // kaon from OmegaMinus has signal - nSigmaOmKa = (bachDeltaTimeAsOmKa - hMeanBachOmKa->Interpolate(cascade.pt())) / hSigmaBachOmKa->Interpolate(cascade.pt()); - } - casctofnsigmas(nSigmaXiLaPi, nSigmaXiLaPr, nSigmaXiPi, nSigmaOmLaPi, nSigmaOmLaPr, nSigmaOmKa); - } - - if (doQA) { - // fill QA histograms for cross-checking - histos.fill(HIST("hArcDebug"), cascade.pt(), lengthCascade - d3d); // for debugging purposes - - if (cascade.dcaV0daughters() < qaV0DCADau && cascade.dcacascdaughters() < qaCascDCADau && cascade.v0cosPA(collision.posX(), collision.posY(), collision.posZ()) > qaV0CosPA && cascade.casccosPA(collision.posX(), collision.posY(), collision.posZ()) > qaCascCosPA) { - if (cascade.sign() < 0) { - if (std::abs(cascade.mXi() - 1.32171) < qaMassWindow && fabs(pTra.tpcNSigmaPr()) < qaTPCNSigma && fabs(nTra.tpcNSigmaPi()) < qaTPCNSigma && fabs(bTra.tpcNSigmaPi()) < qaTPCNSigma) { - histos.fill(HIST("h2dposDeltaTimeAsXiPr"), cascade.pt(), cascade.eta(), posDeltaTimeAsXiPr); - histos.fill(HIST("h2dnegDeltaTimeAsXiPi"), cascade.pt(), cascade.eta(), negDeltaTimeAsXiPi); - histos.fill(HIST("h2dbachDeltaTimeAsXiPi"), cascade.pt(), cascade.eta(), bachDeltaTimeAsXiPi); - if (doQANSigma) { - histos.fill(HIST("h2dNSigmaXiLaPi"), cascade.pt(), nSigmaXiLaPi); - histos.fill(HIST("h2dNSigmaXiLaPr"), cascade.pt(), nSigmaXiLaPr); - histos.fill(HIST("h2dNSigmaXiPi"), cascade.pt(), nSigmaXiPi); - } - } - if (std::abs(cascade.mOmega() - 1.67245) < qaMassWindow && fabs(pTra.tpcNSigmaPr()) < qaTPCNSigma && fabs(nTra.tpcNSigmaPi()) < qaTPCNSigma && fabs(bTra.tpcNSigmaKa()) < qaTPCNSigma) { - histos.fill(HIST("h2dposDeltaTimeAsOmPr"), cascade.pt(), cascade.eta(), posDeltaTimeAsOmPr); - histos.fill(HIST("h2dnegDeltaTimeAsOmPi"), cascade.pt(), cascade.eta(), negDeltaTimeAsOmPi); - histos.fill(HIST("h2dbachDeltaTimeAsOmKa"), cascade.pt(), cascade.eta(), bachDeltaTimeAsOmKa); - if (doQANSigma) { - histos.fill(HIST("h2dNSigmaOmLaPi"), cascade.pt(), nSigmaOmLaPi); - histos.fill(HIST("h2dNSigmaOmLaPr"), cascade.pt(), nSigmaOmLaPr); - histos.fill(HIST("h2dNSigmaOmKa"), cascade.pt(), nSigmaOmKa); - } - } - } else { - if (std::abs(cascade.mXi() - 1.32171) < qaMassWindow && fabs(pTra.tpcNSigmaPi()) < qaTPCNSigma && fabs(nTra.tpcNSigmaPr()) < qaTPCNSigma && fabs(bTra.tpcNSigmaPi()) < qaTPCNSigma) { - histos.fill(HIST("h2dposDeltaTimeAsXiPi"), cascade.pt(), cascade.eta(), posDeltaTimeAsXiPi); - histos.fill(HIST("h2dnegDeltaTimeAsXiPr"), cascade.pt(), cascade.eta(), negDeltaTimeAsXiPr); - histos.fill(HIST("h2dbachDeltaTimeAsXiPi"), cascade.pt(), cascade.eta(), bachDeltaTimeAsXiPi); - if (doQANSigma) { - histos.fill(HIST("h2dNSigmaXiLaPi"), cascade.pt(), nSigmaXiLaPi); - histos.fill(HIST("h2dNSigmaXiLaPr"), cascade.pt(), nSigmaXiLaPr); - histos.fill(HIST("h2dNSigmaXiPi"), cascade.pt(), nSigmaXiPi); - } - } - if (std::abs(cascade.mOmega() - 1.67245) < qaMassWindow && fabs(pTra.tpcNSigmaPi()) < qaTPCNSigma && fabs(nTra.tpcNSigmaPr()) < qaTPCNSigma && fabs(bTra.tpcNSigmaKa()) < qaTPCNSigma) { - histos.fill(HIST("h2dposDeltaTimeAsOmPi"), cascade.pt(), cascade.eta(), posDeltaTimeAsOmPi); - histos.fill(HIST("h2dnegDeltaTimeAsOmPr"), cascade.pt(), cascade.eta(), negDeltaTimeAsOmPr); - histos.fill(HIST("h2dbachDeltaTimeAsOmKa"), cascade.pt(), cascade.eta(), bachDeltaTimeAsOmKa); - if (doQANSigma) { - histos.fill(HIST("h2dNSigmaOmLaPi"), cascade.pt(), nSigmaOmLaPi); - histos.fill(HIST("h2dNSigmaOmLaPr"), cascade.pt(), nSigmaOmLaPr); - histos.fill(HIST("h2dNSigmaOmKa"), cascade.pt(), nSigmaOmKa); - } - } - } - } - } - } - - void processStandardData(aod::Collisions const& collisions, CascOriginalDatas const& Cascades, TracksWithAllExtras const&, aod::BCsWithTimestamps const& /*bcs*/) - { - // Fire up CCDB with first collision in record. If no collisions, bypass - if (useCustomRunNumber || collisions.size() < 1) { - initCCDB(manualRunNumber); - } else { - auto collision = collisions.begin(); - auto bc = collision.bc_as(); - initCCDB(bc.runNumber()); - } - - for (const auto& collision : collisions) { - // Do analysis with collision-grouped V0s, retain full collision information - const uint64_t collIdx = collision.globalIndex(); - auto CascTable_thisCollision = Cascades.sliceBy(perCollisionOriginal, collIdx); - // cascade table sliced - for (auto const& cascade : CascTable_thisCollision) { - // de-reference interlinks by hand for derived data - auto pTra = cascade.posTrack_as(); - auto nTra = cascade.negTrack_as(); - auto bTra = cascade.bachelor_as(); - - processCascadeCandidate(collision, cascade, pTra, nTra, bTra); - } - } - } - - void processDerivedData(soa::Join const& collisions, CascDerivedDatas const& Cascades, dauTracks const&) - { - // 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()); - } - - for (const auto& collision : collisions) { - // Do analysis with collision-grouped V0s, retain full collision information - const uint64_t collIdx = collision.globalIndex(); - auto CascTable_thisCollision = Cascades.sliceBy(perCollisionDerived, collIdx); - // cascade table sliced - for (auto const& cascade : CascTable_thisCollision) { - // de-reference interlinks by hand for derived data - auto pTra = cascade.posTrackExtra_as(); - auto nTra = cascade.negTrackExtra_as(); - auto bTra = cascade.bachTrackExtra_as(); - - processCascadeCandidate(collision, cascade, pTra, nTra, bTra); - } - } - } - - PROCESS_SWITCH(cascadepid, processStandardData, "Process standard data", true); - PROCESS_SWITCH(cascadepid, processDerivedData, "Process derived data", false); -}; - -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; -} diff --git a/PWGLF/TableProducer/Strangeness/lambdakzeropid.cxx b/PWGLF/TableProducer/Strangeness/lambdakzeropid.cxx deleted file mode 100644 index c6041aa0e02..00000000000 --- a/PWGLF/TableProducer/Strangeness/lambdakzeropid.cxx +++ /dev/null @@ -1,656 +0,0 @@ -// 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. -// -// *+-+*+-+*+-+*+-+*+-+*+-+* -// Lambdakzero PID -// *+-+*+-+*+-+*+-+*+-+*+-+* -// -/// \author Nicolò Jacazio -/// \author David Dobrigkeit Chinellato -/// \since 11/05/2023 -/// \brief Table producer for V0 daughter PID info -// -// This task produces daughter PID information for strange daughters -// taking into account the (candidate-by-candidate) time spent as a heavier -// (strange, weakly-decaying) particle. This task is meant to be a test, as -// it hasn't been fully tested yet! Use at your own peril for now :-) - -#include -#include -#include -#include -#include -#include - -#include "Framework/runDataProcessing.h" -#include "Framework/RunningWorkflowInfo.h" -#include "Framework/AnalysisTask.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/ASoAHelpers.h" -#include "DCAFitter/DCAFitterN.h" -#include "ReconstructionDataFormats/Track.h" -#include "Common/Core/RecoDecay.h" -#include "Common/Core/trackUtilities.h" -#include "PWGLF/DataModel/LFStrangenessTables.h" -#include "PWGLF/DataModel/LFStrangenessPIDTables.h" -#include "PWGLF/DataModel/LFParticleIdentification.h" -#include "Common/Core/TrackSelection.h" -#include "Common/DataModel/TrackSelectionTables.h" -#include "DetectorsBase/Propagator.h" -#include "DetectorsBase/GeometryManager.h" -#include "DataFormatsParameters/GRPObject.h" -#include "DataFormatsParameters/GRPMagField.h" -#include "CCDB/BasicCCDBManager.h" -#include "DataFormatsCalibration/MeanVertexObject.h" -#include "CommonConstants/PhysicsConstants.h" -#include "Common/TableProducer/PID/pidTOFBase.h" -#include "Common/DataModel/PIDResponse.h" - -using namespace o2; -using namespace o2::framework; -using namespace o2::framework::expressions; -using std::array; - -// For original data loops -using V0OriginalDatas = soa::Join; -using TracksWithAllExtras = soa::Join; - -// For derived data analysis -using dauTracks = soa::Join; -using V0DerivedDatas = soa::Join; - -struct lambdakzeropid { - // TOF pid for strangeness (recalculated with topology) - Produces v0tofpid; // table with Nsigmas - Produces v0tofbeta; // table with betas - Produces v0tofdebugs; // table with extra debug information - Produces v0tofnsigmas; // table with nsigmas - - Service ccdb; - - // mean vertex position to be used if no collision associated - o2::dataformats::MeanVertexObject* mVtx = nullptr; - - // For manual sliceBy - Preslice perCollisionOriginal = o2::aod::v0data::collisionId; - ; - Preslice perCollisionDerived = o2::aod::v0data::straCollisionId; - - HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - - // Operation and minimisation criteria - Configurable d_bz_input{"d_bz", -999, "bz field, -999 is automatic"}; - Configurable tofPosition{"tofPosition", 377.934f, "TOF effective (inscribed) radius"}; - Configurable doQA{"doQA", true, "create QA histos"}; - Configurable doQANSigma{"doQANSigma", true, "create QA of Nsigma histos"}; - Configurable qaDCADau{"qaDCADau", 0.5, "DCA daughters (cm) for QA plots"}; - Configurable qaCosPA{"qaCosPA", 0.999, "CosPA for QA plots"}; - Configurable qaMassWindow{"qaMassWindow", 0.005, "Mass window around expected (in GeV/c2) for QA plots"}; - Configurable qaTPCNSigma{"qaTPCNSigma", 5, "TPC N-sigma to apply for qa plots"}; - Configurable doNSigmas{"doNSigmas", false, "calculate TOF N-sigma"}; - - // CCDB options - Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; - Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; - Configurable lutPath{"lutPath", "GLO/Param/MatLUT", "Path of the Lut parametrization"}; - Configurable geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; - Configurable nSigmaPath{"nSigmaPath", "Users/d/ddobrigk/stratof", "Path of information for n-sigma calculation"}; - Configurable mVtxPath{"mVtxPath", "GLO/Calib/MeanVertex", "Path of the mean vertex file"}; - - // manual - Configurable useCustomRunNumber{"useCustomRunNumber", false, "Use custom timestamp"}; - Configurable manualRunNumber{"manualRunNumber", 544122, "manual run number if no collisions saved"}; - - ConfigurableAxis axisEta{"axisEta", {20, -1.0f, +1.0f}, "#eta"}; - ConfigurableAxis axisDeltaTime{"axisDeltaTime", {2000, -1000.0f, +1000.0f}, "delta-time (ps)"}; - ConfigurableAxis axisTime{"axisTime", {200, 0.0f, +20000.0f}, "T (ps)"}; - ConfigurableAxis axisNSigma{"axisNSigma", {200, -10.0f, +10.0f}, "N(#sigma)"}; - ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "p_{T} (GeV/c)"}; - - // for n-sigma calibration - bool nSigmaCalibLoaded; - TList* nSigmaCalibObjects; - TH1 *hMeanPosLaPi, *hSigmaPosLaPi; - TH1 *hMeanPosLaPr, *hSigmaPosLaPr; - TH1 *hMeanNegLaPi, *hSigmaNegLaPi; - TH1 *hMeanNegLaPr, *hSigmaNegLaPr; - TH1 *hMeanPosK0Pi, *hSigmaPosK0Pi; - TH1 *hMeanNegK0Pi, *hSigmaNegK0Pi; - - int mRunNumber; - float d_bz; - float maxSnp; // max sine phi for propagation - float maxStep; // max step size (cm) for propagation - - // enum to keep track of the TOF-related properties for V0s - enum tofEnum { kLength = 0, - kHasTOF, - kNEnums }; - - /// function to calculate track length of this track up to a certain segment of a detector - /// to be used internally in another funcrtion that calculates length until it finds the proper one - /// warning: this could be optimised further for speed - /// \param track the input track - /// \param x1 x of the first point of the detector segment - /// \param y1 y of the first point of the detector segment - /// \param x2 x of the first point of the detector segment - /// \param y2 y of the first point of the detector segment - /// \param magneticField the magnetic field to use when propagating - float trackLengthToSegment(o2::track::TrackPar track, float x1, float y1, float x2, float y2, float magneticField) - { - // don't make use of the track parametrization - float length = -104; - - // causality protection - std::array mom; - track.getPxPyPzGlo(mom); - // get start point - std::array startPoint; - track.getXYZGlo(startPoint); - - // better replaced with scalar momentum check later - // if (((x1 + x2) * mom[0] + (y1 + y2) * mom[1]) < 0.0f) - // return -101; - - // get circle X, Y please - o2::math_utils::CircleXYf_t trcCircle; - float sna, csa; - track.getCircleParams(magneticField, trcCircle, sna, csa); - - // Calculate necessary inner product - float segmentModulus = std::hypot(x2 - x1, y2 - y1); - float alongSegment = ((trcCircle.xC - x1) * (x2 - x1) + (trcCircle.yC - y1) * (y2 - y1)) / segmentModulus; - - // find point of closest approach between segment and circle center - float pcaX = (x2 - x1) * alongSegment / segmentModulus + x1; - float pcaY = (y2 - y1) * alongSegment / segmentModulus + y1; - - float centerDistToPC = std::hypot(pcaX - trcCircle.xC, pcaY - trcCircle.yC); - - // distance pca-to-intercept in multiples of segment modulus (for convenience) - if (centerDistToPC > trcCircle.rC) - return -103; - - float pcaToIntercept = TMath::Sqrt(TMath::Abs(trcCircle.rC * trcCircle.rC - centerDistToPC * centerDistToPC)); - - float interceptX1 = pcaX + (x2 - x1) / segmentModulus * pcaToIntercept; - float interceptY1 = pcaY + (y2 - y1) / segmentModulus * pcaToIntercept; - float interceptX2 = pcaX - (x2 - x1) / segmentModulus * pcaToIntercept; - float interceptY2 = pcaY - (y2 - y1) / segmentModulus * pcaToIntercept; - - float scalarCheck1 = ((x2 - x1) * (interceptX1 - x1) + (y2 - y1) * (interceptY1 - y1)) / segmentModulus; - float scalarCheck2 = ((x2 - x1) * (interceptX2 - x1) + (y2 - y1) * (interceptY2 - y1)) / segmentModulus; - - float cosAngle1 = -1000, sinAngle1 = -1000, modulus1 = -1000; - float cosAngle2 = -1000, sinAngle2 = -1000, modulus2 = -1000; - float length1 = -1000, length2 = -1000; - - modulus1 = std::hypot(interceptX1 - trcCircle.xC, interceptY1 - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC); - cosAngle1 = (interceptX1 - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (interceptY1 - trcCircle.yC) * (startPoint[1] - trcCircle.yC); - sinAngle1 = (interceptX1 - trcCircle.xC) * (startPoint[1] - trcCircle.yC) - (interceptY1 - trcCircle.yC) * (startPoint[0] - trcCircle.xC); - cosAngle1 /= modulus1; - sinAngle1 /= modulus1; - length1 = trcCircle.rC * TMath::ACos(cosAngle1); - length1 *= sqrt(1.0f + track.getTgl() * track.getTgl()); - - modulus2 = std::hypot(interceptX2 - trcCircle.xC, interceptY2 - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC); - cosAngle2 = (interceptX2 - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (interceptY2 - trcCircle.yC) * (startPoint[1] - trcCircle.yC); - sinAngle2 = (interceptX2 - trcCircle.xC) * (startPoint[1] - trcCircle.yC) - (interceptY2 - trcCircle.yC) * (startPoint[0] - trcCircle.xC); - cosAngle2 /= modulus2; - sinAngle2 /= modulus2; - length2 = trcCircle.rC * TMath::ACos(cosAngle2); - length2 *= sqrt(1.0f + track.getTgl() * track.getTgl()); - - // rotate transverse momentum vector such that it is at intercepts - float angle1 = TMath::ACos(cosAngle1); - if (sinAngle1 < 0) - angle1 *= -1.0f; - float px1 = +TMath::Cos(angle1) * mom[0] + TMath::Sin(angle1) * mom[1]; - float py1 = -TMath::Sin(angle1) * mom[0] + TMath::Cos(angle1) * mom[1]; - - float angle2 = TMath::ACos(cosAngle2); - if (sinAngle2 < 0) - angle2 *= -1.0f; - float px2 = +TMath::Cos(angle2) * mom[0] + TMath::Sin(angle2) * mom[1]; - float py2 = -TMath::Sin(angle2) * mom[0] + TMath::Cos(angle2) * mom[1]; - - float midSegX = 0.5f * (x2 + x1); - float midSegY = 0.5f * (y2 + y1); - - float scalarMomentumCheck1 = px1 * midSegX + py1 * midSegY; - float scalarMomentumCheck2 = px2 * midSegX + py2 * midSegY; - - if (scalarCheck1 > 0.0f && scalarCheck1 < segmentModulus && scalarMomentumCheck1 > 0.0f) { - length = length1; - // X = interceptX1; Y = interceptY1; - } - if (scalarCheck2 > 0.0f && scalarCheck2 < segmentModulus && scalarMomentumCheck2 > 0.0f) { - length = length2; - // X = interceptX2; Y = interceptY2; - } - return length; - } - - /// function to calculate track length of this track up to a certain segmented detector - /// \param track the input track - /// \param magneticField the magnetic field to use when propagating - float findInterceptLength(o2::track::TrackPar track, float magneticField) - { - float length = 1e+6; - for (int iSeg = 0; iSeg < 18; iSeg++) { - // Detector segmentation loop - float segmentAngle = 20.0f / 180.0f * TMath::Pi(); - float theta = static_cast(iSeg) * 20.0f / 180.0f * TMath::Pi(); - float halfWidth = tofPosition * TMath::Tan(0.5f * segmentAngle); - float x1 = TMath::Cos(theta) * (-halfWidth) + TMath::Sin(theta) * tofPosition; - float y1 = -TMath::Sin(theta) * (-halfWidth) + TMath::Cos(theta) * tofPosition; - float x2 = TMath::Cos(theta) * (+halfWidth) + TMath::Sin(theta) * tofPosition; - float y2 = -TMath::Sin(theta) * (+halfWidth) + TMath::Cos(theta) * tofPosition; - float thisLength = trackLengthToSegment(track, x1, y1, x2, y2, magneticField); - if (thisLength < length && thisLength > 0) - length = thisLength; - } - if (length > 1e+5) - length = -100; // force negative to avoid misunderstandings - return length; - } - - void init(InitContext&) - { - nSigmaCalibLoaded = false; - nSigmaCalibObjects = nullptr; - - // for n-sigma calibration - hMeanPosLaPi = nullptr; - hSigmaPosLaPi = nullptr; - hMeanPosLaPr = nullptr; - hSigmaPosLaPr = nullptr; - hMeanNegLaPi = nullptr; - hSigmaNegLaPi = nullptr; - hMeanNegLaPr = nullptr; - hSigmaNegLaPr = nullptr; - hMeanPosK0Pi = nullptr; - hSigmaNegK0Pi = nullptr; - hMeanNegK0Pi = nullptr; - hSigmaNegK0Pi = nullptr; - - mRunNumber = 0; - d_bz = 0; - maxSnp = 0.85f; // could be changed later - maxStep = 2.00f; // could be changed later - - ccdb->setURL(ccdburl); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - ccdb->setFatalWhenNull(false); - - // per event - histos.add("hCandidateCounter", "hCandidateCounter", kTH1F, {{500, -0.5f, 499.5f}}); - - // measured vs expected total time QA - if (doQA) { - histos.add("h2dProtonMeasuredVsExpected", "h2dProtonMeasuredVsExpected", {HistType::kTH2F, {axisTime, axisTime}}); - histos.add("h2dPionMeasuredVsExpected", "h2dPionMeasuredVsExpected", {HistType::kTH2F, {axisTime, axisTime}}); - - // standard deltaTime values - histos.add("h2dDeltaTimePositiveLambdaPi", "h2dDeltaTimePositiveLambdaPi", {HistType::kTH3F, {axisPt, axisEta, axisDeltaTime}}); - histos.add("h2dDeltaTimeNegativeLambdaPi", "h2dDeltaTimeNegativeLambdaPi", {HistType::kTH3F, {axisPt, axisEta, axisDeltaTime}}); - histos.add("h2dDeltaTimePositiveLambdaPr", "h2dDeltaTimePositiveLambdaPr", {HistType::kTH3F, {axisPt, axisEta, axisDeltaTime}}); - histos.add("h2dDeltaTimeNegativeLambdaPr", "h2dDeltaTimeNegativeLambdaPr", {HistType::kTH3F, {axisPt, axisEta, axisDeltaTime}}); - histos.add("h2dDeltaTimePositiveK0ShortPi", "h2dDeltaTimePositiveK0ShortPi", {HistType::kTH3F, {axisPt, axisEta, axisDeltaTime}}); - histos.add("h2dDeltaTimeNegativeK0ShortPi", "h2dDeltaTimeNegativeK0ShortPi", {HistType::kTH3F, {axisPt, axisEta, axisDeltaTime}}); - - histos.add("h2dPositiveTOFProperties", "h2dPositiveTOFProperties", {HistType::kTH2F, {axisPt, {4, -0.5, 3.5f}}}); - histos.add("h2dNegativeTOFProperties", "h2dNegativeTOFProperties", {HistType::kTH2F, {axisPt, {4, -0.5, 3.5f}}}); - - if (doQANSigma) { - // standard NSigma values - histos.add("h2dNSigmaPositiveLambdaPi", "h2dNSigmaPositiveLambdaPi", {HistType::kTH2F, {axisPt, axisNSigma}}); - histos.add("h2dNSigmaNegativeLambdaPi", "h2dNSigmaNegativeLambdaPi", {HistType::kTH2F, {axisPt, axisNSigma}}); - histos.add("h2dNSigmaPositiveLambdaPr", "h2dNSigmaPositiveLambdaPr", {HistType::kTH2F, {axisPt, axisNSigma}}); - histos.add("h2dNSigmaNegativeLambdaPr", "h2dNSigmaNegativeLambdaPr", {HistType::kTH2F, {axisPt, axisNSigma}}); - histos.add("h2dNSigmaPositiveK0ShortPi", "h2dNSigmaPositiveK0ShortPi", {HistType::kTH2F, {axisPt, axisNSigma}}); - histos.add("h2dNSigmaNegativeK0ShortPi", "h2dNSigmaNegativeK0ShortPi", {HistType::kTH2F, {axisPt, axisNSigma}}); - } - - // delta lambda decay time - histos.add("h2dLambdaDeltaDecayTime", "h2dLambdaDeltaDecayTime", {HistType::kTH2F, {axisPt, axisDeltaTime}}); - } - } - - void initCCDB(int runNumber) - { - if (mRunNumber == runNumber) { - return; - } - - // In case override, don't proceed, please - no CCDB access required - if (d_bz_input > -990) { - d_bz = d_bz_input; - o2::parameters::GRPMagField grpmag; - if (fabs(d_bz) > 1e-5) { - grpmag.setL3Current(30000.f / (d_bz / 5.0f)); - } - o2::base::Propagator::initFieldFromGRP(&grpmag); - mVtx = ccdb->getForRun(mVtxPath, runNumber); - mRunNumber = runNumber; - return; - } - - o2::parameters::GRPObject* grpo = ccdb->getForRun(grpPath, runNumber); - o2::parameters::GRPMagField* grpmag = 0x0; - if (grpo) { - o2::base::Propagator::initFieldFromGRP(grpo); - // Fetch magnetic field from ccdb for current collision - d_bz = grpo->getNominalL3Field(); - LOG(info) << "Retrieved GRP for run " << runNumber << " with magnetic field of " << d_bz << " kZG"; - } else { - grpmag = ccdb->getForRun(grpmagPath, runNumber); - if (!grpmag) { - LOG(fatal) << "Got nullptr from CCDB for path " << grpmagPath << " of object GRPMagField and " << grpPath << " of object GRPObject for run " << runNumber; - } - o2::base::Propagator::initFieldFromGRP(grpmag); - // Fetch magnetic field from ccdb for current collision - d_bz = std::lround(5.f * grpmag->getL3Current() / 30000.f); - mVtx = ccdb->getForRun(mVtxPath, runNumber); - LOG(info) << "Retrieved GRP for run " << runNumber << " with magnetic field of " << d_bz << " kZG"; - } - - // if TOF Nsigma desired - if (doNSigmas) { - nSigmaCalibObjects = ccdb->getForRun(nSigmaPath, runNumber); - if (nSigmaCalibObjects) { - LOGF(info, "loaded TList with this many objects: %i", nSigmaCalibObjects->GetEntries()); - - hMeanPosLaPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanPosLaPi")); - hMeanPosLaPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanPosLaPr")); - hMeanNegLaPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanNegLaPi")); - hMeanNegLaPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanNegLaPr")); - hMeanPosK0Pi = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanPosK0Pi")); - hMeanNegK0Pi = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanNegK0Pi")); - - hSigmaPosLaPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaPosLaPi")); - hSigmaPosLaPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaPosLaPr")); - hSigmaNegLaPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaNegLaPi")); - hSigmaNegLaPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaNegLaPr")); - hSigmaPosK0Pi = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaPosK0Pi")); - hSigmaNegK0Pi = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaNegK0Pi")); - - if (!hMeanPosLaPi) - LOG(info) << "Problems finding mean histogram hMeanPosLaPi!"; - if (!hMeanPosLaPr) - LOG(info) << "Problems finding mean histogram hMeanPosLaPr!"; - if (!hMeanNegLaPi) - LOG(info) << "Problems finding mean histogram hMeanNegLaPi!"; - if (!hMeanNegLaPr) - LOG(info) << "Problems finding mean histogram hMeanNegLaPr!"; - if (!hMeanPosK0Pi) - LOG(info) << "Problems finding mean histogram hMeanPosK0Pi!"; - if (!hMeanNegK0Pi) - LOG(info) << "Problems finding mean histogram hMeanNegK0Pi!"; - if (!hSigmaPosK0Pi || !hSigmaNegK0Pi || !hSigmaPosLaPi || !hSigmaPosLaPr || !hSigmaNegLaPi || !hSigmaNegLaPr) { - LOG(info) << "Problems finding sigma histograms!"; - } - } - } - mRunNumber = runNumber; - } - - float velocity(float lMomentum, float lMass) - { - // Momentum p and mass m -> returns speed in centimeters per picosecond - // Useful for TOF calculations - float lA = (lMomentum / lMass) * (lMomentum / lMass); - 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) - { - // 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); - o2::track::TrackPar negTrack = o2::track::TrackPar({v0.x(), v0.y(), v0.z()}, {v0.pxneg(), v0.pyneg(), v0.pzneg()}, -1); - - float deltaTimePositiveLambdaPi = -1e+6; - float deltaTimeNegativeLambdaPi = -1e+6; - float deltaTimePositiveLambdaPr = -1e+6; - float deltaTimeNegativeLambdaPr = -1e+6; - float deltaTimePositiveK0ShortPi = -1e+6; - float deltaTimeNegativeK0ShortPi = -1e+6; - - float nSigmaPositiveLambdaPi = -1e+3; - float nSigmaPositiveLambdaPr = -1e+3; - float nSigmaNegativeLambdaPi = -1e+3; - float nSigmaNegativeLambdaPr = -1e+3; - float nSigmaPositiveK0ShortPi = -1e+3; - float nSigmaNegativeK0ShortPi = -1e+3; - - 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 lengthPositive = findInterceptLength(posTrack, d_bz); // FIXME: tofPosition ok? adjust? - float lengthNegative = findInterceptLength(negTrack, d_bz); // FIXME: tofPosition ok? adjust? - float timePositivePr = lengthPositive / velocityPositivePr; - float timePositivePi = lengthPositive / velocityPositivePi; - float timeNegativePr = lengthNegative / velocityNegativePr; - float timeNegativePi = lengthNegative / velocityNegativePi; - - if (pTra.hasTOF() && lengthPositive > 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() && lengthNegative > 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 (lengthPositive > 0) - posProperties = posProperties | (static_cast(1) << kLength); - if (pTra.hasTOF()) - posProperties = posProperties | (static_cast(1) << kHasTOF); - int negProperties = 0; - if (lengthNegative > 0) - negProperties = negProperties | (static_cast(1) << kLength); - if (nTra.hasTOF()) - negProperties = negProperties | (static_cast(1) << kHasTOF); - - histos.fill(HIST("h2dPositiveTOFProperties"), v0.pt(), posProperties); - histos.fill(HIST("h2dNegativeTOFProperties"), v0.pt(), negProperties); - } - - float deltaDecayTimeLambda = -10e+4; - float deltaDecayTimeAntiLambda = -10e+4; - float deltaDecayTimeK0Short = -10e+4; - if (nTra.hasTOF() && pTra.hasTOF() > 0 && lengthPositive > 0 && lengthNegative > 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); - } - - // calculate betas - - 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 betaLambda = -1e+6; - float betaAntiLambda = -1e+6; - float betaK0Short = -1e+6; - - if (nTra.hasTOF() && pTra.hasTOF()) { - betaLambda = (lengthV0 / decayTimeLambda) / 0.0299792458; - betaAntiLambda = (lengthV0 / decayTimeAntiLambda) / 0.0299792458; - betaK0Short = (lengthV0 / decayTimeK0Short) / 0.0299792458; - } - - 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) { - // sweep through all viable hypotheses and produce N-sigma - - if (deltaTimePositiveLambdaPi > -1e+5) - nSigmaPositiveLambdaPi = (deltaTimePositiveLambdaPi - hMeanPosLaPi->Interpolate(v0.pt())) / hSigmaPosLaPi->Interpolate(v0.pt()); - if (deltaTimePositiveLambdaPr > -1e+5) - nSigmaPositiveLambdaPr = (deltaTimePositiveLambdaPr - hMeanPosLaPr->Interpolate(v0.pt())) / hSigmaPosLaPr->Interpolate(v0.pt()); - if (deltaTimeNegativeLambdaPi > -1e+5) - nSigmaNegativeLambdaPi = (deltaTimeNegativeLambdaPi - hMeanNegLaPi->Interpolate(v0.pt())) / hSigmaNegLaPi->Interpolate(v0.pt()); - if (deltaTimeNegativeLambdaPr > -1e+5) - nSigmaNegativeLambdaPr = (deltaTimeNegativeLambdaPr - hMeanNegLaPr->Interpolate(v0.pt())) / hSigmaNegLaPr->Interpolate(v0.pt()); - if (deltaTimePositiveK0ShortPi > -1e+5) - nSigmaPositiveK0ShortPi = (deltaTimePositiveK0ShortPi - hMeanPosK0Pi->Interpolate(v0.pt())) / hSigmaPosK0Pi->Interpolate(v0.pt()); - if (deltaTimeNegativeK0ShortPi > -1e+5) - nSigmaNegativeK0ShortPi = (deltaTimeNegativeK0ShortPi - hMeanNegK0Pi->Interpolate(v0.pt())) / hSigmaNegK0Pi->Interpolate(v0.pt()); - - v0tofnsigmas( - nSigmaPositiveLambdaPr, nSigmaNegativeLambdaPi, - nSigmaNegativeLambdaPr, nSigmaPositiveLambdaPi, - nSigmaPositiveK0ShortPi, nSigmaNegativeK0ShortPi); - } - - if (doQA) { - if (pTra.hasTOF()) { - histos.fill(HIST("h2dProtonMeasuredVsExpected"), - (timeLambda + timePositivePr), - (pTra.tofSignal() - pTra.tofEvTime())); - if (v0.v0cosPA() > qaCosPA && v0.dcaV0daughters() < qaDCADau) { - if (std::abs(v0.mLambda() - 1.115683) < qaMassWindow && fabs(pTra.tpcNSigmaPr()) < qaTPCNSigma && fabs(nTra.tpcNSigmaPi()) < qaTPCNSigma) { - histos.fill(HIST("h2dDeltaTimePositiveLambdaPr"), v0.pt(), v0.eta(), deltaTimePositiveLambdaPr); - if (doQANSigma) - histos.fill(HIST("h2dNSigmaPositiveLambdaPr"), v0.pt(), nSigmaPositiveLambdaPr); - } - if (std::abs(v0.mAntiLambda() - 1.115683) < qaMassWindow && fabs(pTra.tpcNSigmaPi()) < qaTPCNSigma && fabs(nTra.tpcNSigmaPr()) < qaTPCNSigma) { - histos.fill(HIST("h2dDeltaTimePositiveLambdaPi"), v0.pt(), v0.eta(), deltaTimePositiveLambdaPi); - if (doQANSigma) - histos.fill(HIST("h2dNSigmaPositiveLambdaPi"), v0.pt(), nSigmaPositiveLambdaPi); - } - if (std::abs(v0.mK0Short() - 0.497) < qaMassWindow && fabs(pTra.tpcNSigmaPi()) < qaTPCNSigma && fabs(nTra.tpcNSigmaPi()) < qaTPCNSigma) { - histos.fill(HIST("h2dDeltaTimePositiveK0ShortPi"), v0.pt(), v0.eta(), deltaTimePositiveK0ShortPi); - if (doQANSigma) - histos.fill(HIST("h2dNSigmaPositiveK0ShortPi"), v0.pt(), nSigmaPositiveK0ShortPi); - } - } - } - - if (nTra.hasTOF()) { - histos.fill(HIST("h2dPionMeasuredVsExpected"), - (timeLambda + timeNegativePi), - (nTra.tofSignal() - nTra.tofEvTime())); - if (v0.v0cosPA() > qaCosPA && v0.dcaV0daughters() < qaDCADau) { - if (std::abs(v0.mLambda() - 1.115683) < qaMassWindow && fabs(pTra.tpcNSigmaPr()) < qaTPCNSigma && fabs(nTra.tpcNSigmaPi()) < qaTPCNSigma) { - histos.fill(HIST("h2dDeltaTimeNegativeLambdaPi"), v0.pt(), v0.eta(), deltaTimeNegativeLambdaPi); - if (doQANSigma) - histos.fill(HIST("h2dNSigmaNegativeLambdaPi"), v0.pt(), nSigmaNegativeLambdaPi); - } - if (std::abs(v0.mAntiLambda() - 1.115683) < qaMassWindow && fabs(pTra.tpcNSigmaPi()) < qaTPCNSigma && fabs(nTra.tpcNSigmaPr()) < qaTPCNSigma) { - histos.fill(HIST("h2dDeltaTimeNegativeLambdaPr"), v0.pt(), v0.eta(), deltaTimeNegativeLambdaPr); - if (doQANSigma) - histos.fill(HIST("h2dNSigmaNegativeLambdaPr"), v0.pt(), nSigmaNegativeLambdaPr); - } - if (std::abs(v0.mK0Short() - 0.497) < qaMassWindow && fabs(pTra.tpcNSigmaPi()) < qaTPCNSigma && fabs(nTra.tpcNSigmaPi()) < qaTPCNSigma) { - histos.fill(HIST("h2dDeltaTimeNegativeK0ShortPi"), v0.pt(), v0.eta(), deltaTimeNegativeK0ShortPi); - if (doQANSigma) - histos.fill(HIST("h2dNSigmaNegativeK0ShortPi"), v0.pt(), nSigmaNegativeK0ShortPi); - } - } - } - // delta lambda decay time - histos.fill(HIST("h2dLambdaDeltaDecayTime"), v0.pt(), deltaDecayTimeLambda); - } - } - - void processStandardData(aod::Collisions const& collisions, V0OriginalDatas const& V0s, TracksWithAllExtras const&, aod::BCsWithTimestamps const& /*bcs*/) - { - // Fire up CCDB with first collision in record. If no collisions, bypass - if (useCustomRunNumber || collisions.size() < 1) { - initCCDB(manualRunNumber); - } else { - auto collision = collisions.begin(); - auto bc = collision.bc_as(); - initCCDB(bc.runNumber()); - } - - 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()}); - } - - auto pTra = V0.posTrack_as(); - auto nTra = V0.negTrack_as(); - processV0Candidate(primaryVertex, V0, pTra, nTra); - } - } - - void processDerivedData(soa::Join const& collisions, V0DerivedDatas const& V0s, dauTracks const&) - { - // 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()); - } - - 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()}); - } - - auto pTra = V0.posTrackExtra_as(); - auto nTra = V0.negTrackExtra_as(); - processV0Candidate(primaryVertex, V0, pTra, nTra); - } - } - - PROCESS_SWITCH(lambdakzeropid, processStandardData, "Process standard data", true); - PROCESS_SWITCH(lambdakzeropid, processDerivedData, "Process derived data", false); -}; - -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; -} diff --git a/PWGLF/TableProducer/Strangeness/strangenesstofpid.cxx b/PWGLF/TableProducer/Strangeness/strangenesstofpid.cxx new file mode 100644 index 00000000000..75cb771dbea --- /dev/null +++ b/PWGLF/TableProducer/Strangeness/strangenesstofpid.cxx @@ -0,0 +1,1012 @@ +// 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. +// +// *+-+*+-+*+-+*+-+*+-+*+-+* +// Lambdakzero PID +// *+-+*+-+*+-+*+-+*+-+*+-+* +// +/// \author Nicolò Jacazio +/// \author David Dobrigkeit Chinellato +/// \since 11/05/2023 +/// \brief Table producer for V0 daughter PID info +// +// This task produces daughter PID information for strange daughters +// taking into account the (candidate-by-candidate) time spent as a heavier +// (strange, weakly-decaying) particle. This task is meant to be a test, as +// it hasn't been fully tested yet! Use at your own peril for now :-) + +#include "TableHelper.h" + +#include "PWGLF/DataModel/LFParticleIdentification.h" +#include "PWGLF/DataModel/LFStrangenessPIDTables.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" + +#include "Common/Core/RecoDecay.h" +#include "Common/Core/TrackSelection.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/PIDResponse.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Common/TableProducer/PID/pidTOFBase.h" + +#include "CCDB/BasicCCDBManager.h" +#include "CommonConstants/PhysicsConstants.h" +#include "DCAFitter/DCAFitterN.h" +#include "DataFormatsCalibration/MeanVertexObject.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/Propagator.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/RunningWorkflowInfo.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/Track.h" + +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using std::array; + +// For original data loops +using V0OriginalDatas = soa::Join; +using CascOriginalDatas = soa::Join; +using TracksWithAllExtras = soa::Join; + +// For derived data analysis +using dauTracks = soa::Join; +using V0DerivedDatas = soa::Join; +using CascDerivedDatas = soa::Join; + +struct strangenesstofpid { + // TOF pid for strangeness (recalculated with topology) + Produces v0tofpid; // table with Nsigmas + Produces v0tofbeta; // table with betas + Produces v0tofdebugs; // table with extra debug information + Produces v0tofnsigmas; // table with nsigmas + Produces casctofpids; // cascades: table with base info + Produces casctofnsigmas; // cascades: table with Nsigmas + + Service ccdb; + + // mean vertex position to be used if no collision associated + o2::dataformats::MeanVertexObject* mVtx = nullptr; + + HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + // master switches + Configurable calculateV0s{"calculateV0s", -1, "calculate V0-related TOF PID (0: no, 1: yes, -1: auto)"}; + Configurable calculateCascades{"calculateCascades", -1, "calculate cascade-related TOF PID (0: no, 1: yes, -1: auto)"}; + + // Operation and minimisation criteria + Configurable d_bz_input{"d_bz", -999, "bz field, -999 is automatic"}; + Configurable tofPosition{"tofPosition", 377.934f, "TOF effective (inscribed) radius"}; + Configurable doQA{"doQA", true, "create QA histos"}; + Configurable doNSigmas{"doNSigmas", false, "calculate TOF N-sigma"}; + Configurable doQANSigma{"doQANSigma", true, "create QA of Nsigma histos"}; + + // configurables related to V0s + struct : ConfigurableGroup { + std::string prefix = "v0Calibration"; + Configurable qaDCADau{"qaDCADau", 0.5, "DCA daughters (cm) for QA plots"}; + Configurable qaCosPA{"qaCosPA", 0.999, "CosPA for QA plots"}; + Configurable qaMassWindow{"qaMassWindow", 0.005, "Mass window around expected (in GeV/c2) for QA plots"}; + Configurable qaTPCNSigma{"qaTPCNSigma", 5, "TPC N-sigma to apply for qa plots"}; + } v0Group; + + // configurables related to V0s + struct : ConfigurableGroup { + std::string prefix = "cascadeCalibration"; + Configurable qaV0DCADau{"qaV0DCADau", 0.5, "DCA daughters (cm) for QA plots"}; + Configurable qaCascDCADau{"qaCascDCADau", 0.5, "DCA daughters (cm) for QA plots"}; + Configurable qaV0CosPA{"qaV0CosPA", 0.995, "CosPA for QA plots"}; + Configurable qaCascCosPA{"qaCascCosPA", 0.995, "CosPA for QA plots"}; + Configurable qaMassWindow{"qaMassWindow", 0.005, "Mass window around expected (in GeV/c2) for QA plots"}; + Configurable qaTPCNSigma{"qaTPCNSigma", 5, "TPC N-sigma to apply for qa plots"}; + } cascadeGroup; + + // CCDB options + Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; + Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; + Configurable lutPath{"lutPath", "GLO/Param/MatLUT", "Path of the Lut parametrization"}; + Configurable nSigmaPath{"nSigmaPath", "Users/d/ddobrigk/stratof", "Path of information for n-sigma calculation"}; + Configurable mVtxPath{"mVtxPath", "GLO/Calib/MeanVertex", "Path of the mean vertex file"}; + + // manual + Configurable useCustomRunNumber{"useCustomRunNumber", false, "Use custom timestamp"}; + Configurable manualRunNumber{"manualRunNumber", 544122, "manual run number if no collisions saved"}; + + ConfigurableAxis axisEta{"axisEta", {20, -1.0f, +1.0f}, "#eta"}; + ConfigurableAxis axisDeltaTime{"axisDeltaTime", {2000, -1000.0f, +1000.0f}, "delta-time (ps)"}; + ConfigurableAxis axisTime{"axisTime", {200, 10000.0f, +20000.0f}, "T (ps)"}; + ConfigurableAxis axisNSigma{"axisNSigma", {200, -10.0f, +10.0f}, "N(#sigma)"}; + ConfigurableAxis axisExpectedOverMeasured{"axisExpectedOverMeasured", {200, 0.9f, 2.9f}, "T_{exp}/T_{meas}"}; + + // master p axis + ConfigurableAxis axisP{"axisP", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "p_{T} (GeV/c)"}; + + // for zooming in at low values only (e-loss studies and effective correction) + ConfigurableAxis axisSmallP{"axisSmallP", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f}, "p_{T} (GeV/c)"}; + + // for n-sigma calibration + bool nSigmaCalibLoaded; + TList* nSigmaCalibObjects; + TH1 *hMeanPosLaPi, *hSigmaPosLaPi; + TH1 *hMeanPosLaPr, *hSigmaPosLaPr; + TH1 *hMeanNegLaPi, *hSigmaNegLaPi; + TH1 *hMeanNegLaPr, *hSigmaNegLaPr; + TH1 *hMeanPosK0Pi, *hSigmaPosK0Pi; + TH1 *hMeanNegK0Pi, *hSigmaNegK0Pi; + TH1 *hMeanPosXiPi, *hSigmaPosXiPi; + TH1 *hMeanPosXiPr, *hSigmaPosXiPr; + TH1 *hMeanNegXiPi, *hSigmaNegXiPi; + TH1 *hMeanNegXiPr, *hSigmaNegXiPr; + TH1 *hMeanBachXiPi, *hSigmaBachXiPi; + TH1 *hMeanPosOmPi, *hSigmaPosOmPi; + TH1 *hMeanPosOmPr, *hSigmaPosOmPr; + TH1 *hMeanNegOmPi, *hSigmaNegOmPi; + TH1 *hMeanNegOmPr, *hSigmaNegOmPr; + TH1 *hMeanBachOmKa, *hSigmaBachOmKa; + + int mRunNumber; + float d_bz; + float maxSnp; // max sine phi for propagation + float maxStep; // max step size (cm) for propagation + + // enum to keep track of the TOF-related properties for V0s + enum tofEnum { kLength = 0, + kHasTOF, + kNEnums }; + + /// function to calculate track length of this track up to a certain segment of a detector + /// to be used internally in another funcrtion that calculates length until it finds the proper one + /// warning: this could be optimised further for speed + /// \param track the input track + /// \param x1 x of the first point of the detector segment + /// \param y1 y of the first point of the detector segment + /// \param x2 x of the first point of the detector segment + /// \param y2 y of the first point of the detector segment + /// \param magneticField the magnetic field to use when propagating + float trackLengthToSegment(o2::track::TrackPar track, float x1, float y1, float x2, float y2, float magneticField) + { + // don't make use of the track parametrization + float length = -104; + + // causality protection + std::array mom; + track.getPxPyPzGlo(mom); + // get start point + std::array startPoint; + track.getXYZGlo(startPoint); + + // better replaced with scalar momentum check later + // if (((x1 + x2) * mom[0] + (y1 + y2) * mom[1]) < 0.0f) + // return -101; + + // get circle X, Y please + o2::math_utils::CircleXYf_t trcCircle; + float sna, csa; + track.getCircleParams(magneticField, trcCircle, sna, csa); + + // Calculate necessary inner product + float segmentModulus = std::hypot(x2 - x1, y2 - y1); + float alongSegment = ((trcCircle.xC - x1) * (x2 - x1) + (trcCircle.yC - y1) * (y2 - y1)) / segmentModulus; + + // find point of closest approach between segment and circle center + float pcaX = (x2 - x1) * alongSegment / segmentModulus + x1; + float pcaY = (y2 - y1) * alongSegment / segmentModulus + y1; + + float centerDistToPC = std::hypot(pcaX - trcCircle.xC, pcaY - trcCircle.yC); + + // distance pca-to-intercept in multiples of segment modulus (for convenience) + if (centerDistToPC > trcCircle.rC) + return -103; + + float pcaToIntercept = TMath::Sqrt(TMath::Abs(trcCircle.rC * trcCircle.rC - centerDistToPC * centerDistToPC)); + + float interceptX1 = pcaX + (x2 - x1) / segmentModulus * pcaToIntercept; + float interceptY1 = pcaY + (y2 - y1) / segmentModulus * pcaToIntercept; + float interceptX2 = pcaX - (x2 - x1) / segmentModulus * pcaToIntercept; + float interceptY2 = pcaY - (y2 - y1) / segmentModulus * pcaToIntercept; + + float scalarCheck1 = ((x2 - x1) * (interceptX1 - x1) + (y2 - y1) * (interceptY1 - y1)) / segmentModulus; + float scalarCheck2 = ((x2 - x1) * (interceptX2 - x1) + (y2 - y1) * (interceptY2 - y1)) / segmentModulus; + + float cosAngle1 = -1000, sinAngle1 = -1000, modulus1 = -1000; + float cosAngle2 = -1000, sinAngle2 = -1000, modulus2 = -1000; + float length1 = -1000, length2 = -1000; + + modulus1 = std::hypot(interceptX1 - trcCircle.xC, interceptY1 - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC); + cosAngle1 = (interceptX1 - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (interceptY1 - trcCircle.yC) * (startPoint[1] - trcCircle.yC); + sinAngle1 = (interceptX1 - trcCircle.xC) * (startPoint[1] - trcCircle.yC) - (interceptY1 - trcCircle.yC) * (startPoint[0] - trcCircle.xC); + cosAngle1 /= modulus1; + sinAngle1 /= modulus1; + length1 = trcCircle.rC * TMath::ACos(cosAngle1); + length1 *= sqrt(1.0f + track.getTgl() * track.getTgl()); + + modulus2 = std::hypot(interceptX2 - trcCircle.xC, interceptY2 - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC); + cosAngle2 = (interceptX2 - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (interceptY2 - trcCircle.yC) * (startPoint[1] - trcCircle.yC); + sinAngle2 = (interceptX2 - trcCircle.xC) * (startPoint[1] - trcCircle.yC) - (interceptY2 - trcCircle.yC) * (startPoint[0] - trcCircle.xC); + cosAngle2 /= modulus2; + sinAngle2 /= modulus2; + length2 = trcCircle.rC * TMath::ACos(cosAngle2); + length2 *= sqrt(1.0f + track.getTgl() * track.getTgl()); + + // rotate transverse momentum vector such that it is at intercepts + float angle1 = TMath::ACos(cosAngle1); + if (sinAngle1 < 0) + angle1 *= -1.0f; + float px1 = +TMath::Cos(angle1) * mom[0] + TMath::Sin(angle1) * mom[1]; + float py1 = -TMath::Sin(angle1) * mom[0] + TMath::Cos(angle1) * mom[1]; + + float angle2 = TMath::ACos(cosAngle2); + if (sinAngle2 < 0) + angle2 *= -1.0f; + float px2 = +TMath::Cos(angle2) * mom[0] + TMath::Sin(angle2) * mom[1]; + float py2 = -TMath::Sin(angle2) * mom[0] + TMath::Cos(angle2) * mom[1]; + + float midSegX = 0.5f * (x2 + x1); + float midSegY = 0.5f * (y2 + y1); + + float scalarMomentumCheck1 = px1 * midSegX + py1 * midSegY; + float scalarMomentumCheck2 = px2 * midSegX + py2 * midSegY; + + if (scalarCheck1 > 0.0f && scalarCheck1 < segmentModulus && scalarMomentumCheck1 > 0.0f) { + length = length1; + // X = interceptX1; Y = interceptY1; + } + if (scalarCheck2 > 0.0f && scalarCheck2 < segmentModulus && scalarMomentumCheck2 > 0.0f) { + length = length2; + // X = interceptX2; Y = interceptY2; + } + return length; + } + + /// function to calculate track length of this track up to a certain segmented detector + /// \param track the input track + /// \param magneticField the magnetic field to use when propagating + float findInterceptLength(o2::track::TrackPar track, float magneticField) + { + float length = 1e+6; + for (int iSeg = 0; iSeg < 18; iSeg++) { + // Detector segmentation loop + float segmentAngle = 20.0f / 180.0f * TMath::Pi(); + float theta = static_cast(iSeg) * 20.0f / 180.0f * TMath::Pi(); + float halfWidth = tofPosition * TMath::Tan(0.5f * segmentAngle); + float x1 = TMath::Cos(theta) * (-halfWidth) + TMath::Sin(theta) * tofPosition; + float y1 = -TMath::Sin(theta) * (-halfWidth) + TMath::Cos(theta) * tofPosition; + float x2 = TMath::Cos(theta) * (+halfWidth) + TMath::Sin(theta) * tofPosition; + float y2 = -TMath::Sin(theta) * (+halfWidth) + TMath::Cos(theta) * tofPosition; + float thisLength = trackLengthToSegment(track, x1, y1, x2, y2, magneticField); + if (thisLength < length && thisLength > 0) + length = thisLength; + } + if (length > 1e+5) + length = -100; // force negative to avoid misunderstandings + return length; + } + + void init(InitContext& initContext) + { + if (calculateV0s.value < 0) { + // check if TOF information is required, enable if so + calculateV0s.value = isTableRequiredInWorkflow(initContext, "V0TOFNSigmas"); + if (calculateV0s.value > 0) { + LOGF(info, "Strangeness TOF PID: V0 calculations enabled automatically"); + } + } + if (calculateCascades.value < 0) { + // check if TOF information is required, enable if so + calculateCascades.value = isTableRequiredInWorkflow(initContext, "CascTOFNSigmas"); + if (calculateCascades.value > 0) { + LOGF(info, "Strangeness TOF PID: Cascade calculations enabled automatically"); + } + } + + nSigmaCalibLoaded = false; + nSigmaCalibObjects = nullptr; + + // for n-sigma calibration + hMeanPosLaPi = nullptr; + hSigmaPosLaPi = nullptr; + hMeanPosLaPr = nullptr; + hSigmaPosLaPr = nullptr; + hMeanNegLaPi = nullptr; + hSigmaNegLaPi = nullptr; + hMeanNegLaPr = nullptr; + hSigmaNegLaPr = nullptr; + hMeanPosK0Pi = nullptr; + hSigmaNegK0Pi = nullptr; + hMeanNegK0Pi = nullptr; + hSigmaNegK0Pi = nullptr; + + mRunNumber = 0; + d_bz = 0; + maxSnp = 0.85f; // could be changed later + maxStep = 2.00f; // could be changed later + + ccdb->setURL(ccdburl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); + + // per event + histos.add("hCandidateCounter", "hCandidateCounter", kTH1F, {{500, -0.5f, 499.5f}}); + + // measured vs expected total time QA + if (doQA) { + // plots for effective eloss corrections - Lambda case + histos.add("h2dProtonMeasuredVsExpected", "h2dProtonMeasuredVsExpected", {HistType::kTH3F, {axisSmallP, axisExpectedOverMeasured, axisEta}}); + histos.add("h2dPionMeasuredVsExpected", "h2dPionMeasuredVsExpected", {HistType::kTH3F, {axisSmallP, axisExpectedOverMeasured, axisEta}}); + + histos.add("hArcDebug", "hArcDebug", kTH2F, {axisP, {50, -5.0f, 10.0f}}); + + // standard deltaTime values + if (calculateV0s.value > 0) { + histos.add("h2dDeltaTimePositiveLambdaPi", "h2dDeltaTimePositiveLambdaPi", {HistType::kTH3F, {axisP, axisEta, axisDeltaTime}}); + histos.add("h2dDeltaTimeNegativeLambdaPi", "h2dDeltaTimeNegativeLambdaPi", {HistType::kTH3F, {axisP, axisEta, axisDeltaTime}}); + histos.add("h2dDeltaTimePositiveLambdaPr", "h2dDeltaTimePositiveLambdaPr", {HistType::kTH3F, {axisP, axisEta, axisDeltaTime}}); + histos.add("h2dDeltaTimeNegativeLambdaPr", "h2dDeltaTimeNegativeLambdaPr", {HistType::kTH3F, {axisP, axisEta, axisDeltaTime}}); + histos.add("h2dDeltaTimePositiveK0ShortPi", "h2dDeltaTimePositiveK0ShortPi", {HistType::kTH3F, {axisP, axisEta, axisDeltaTime}}); + histos.add("h2dDeltaTimeNegativeK0ShortPi", "h2dDeltaTimeNegativeK0ShortPi", {HistType::kTH3F, {axisP, axisEta, axisDeltaTime}}); + } + + if (calculateCascades.value > 0) { + histos.add("h2dposDeltaTimeAsXiPi", "h2dposDeltaTimeAsXiPi", {HistType::kTH3F, {axisP, axisEta, axisDeltaTime}}); + histos.add("h2dposDeltaTimeAsXiPr", "h2dposDeltaTimeAsXiPr", {HistType::kTH3F, {axisP, axisEta, axisDeltaTime}}); + histos.add("h2dnegDeltaTimeAsXiPi", "h2dnegDeltaTimeAsXiPi", {HistType::kTH3F, {axisP, axisEta, axisDeltaTime}}); + histos.add("h2dnegDeltaTimeAsXiPr", "h2dnegDeltaTimeAsXiPr", {HistType::kTH3F, {axisP, axisEta, axisDeltaTime}}); + histos.add("h2dbachDeltaTimeAsXiPi", "h2dbachDeltaTimeAsXiPi", {HistType::kTH3F, {axisP, axisEta, axisDeltaTime}}); + + histos.add("h2dposDeltaTimeAsOmPi", "h2dposDeltaTimeAsOmPi", {HistType::kTH3F, {axisP, axisEta, axisDeltaTime}}); + histos.add("h2dposDeltaTimeAsOmPr", "h2dposDeltaTimeAsOmPr", {HistType::kTH3F, {axisP, axisEta, axisDeltaTime}}); + histos.add("h2dnegDeltaTimeAsOmPi", "h2dnegDeltaTimeAsOmPi", {HistType::kTH3F, {axisP, axisEta, axisDeltaTime}}); + histos.add("h2dnegDeltaTimeAsOmPr", "h2dnegDeltaTimeAsOmPr", {HistType::kTH3F, {axisP, axisEta, axisDeltaTime}}); + histos.add("h2dbachDeltaTimeAsOmKa", "h2dbachDeltaTimeAsOmKa", {HistType::kTH3F, {axisP, axisEta, axisDeltaTime}}); + } + + histos.add("h2dPositiveTOFProperties", "h2dPositiveTOFProperties", {HistType::kTH2F, {axisP, {4, -0.5, 3.5f}}}); + histos.add("h2dNegativeTOFProperties", "h2dNegativeTOFProperties", {HistType::kTH2F, {axisP, {4, -0.5, 3.5f}}}); + + if (doQANSigma) { + if (calculateV0s.value > 0) { + histos.add("h2dNSigmaPositiveLambdaPi", "h2dNSigmaPositiveLambdaPi", {HistType::kTH2F, {axisP, axisNSigma}}); + histos.add("h2dNSigmaNegativeLambdaPi", "h2dNSigmaNegativeLambdaPi", {HistType::kTH2F, {axisP, axisNSigma}}); + histos.add("h2dNSigmaPositiveLambdaPr", "h2dNSigmaPositiveLambdaPr", {HistType::kTH2F, {axisP, axisNSigma}}); + histos.add("h2dNSigmaNegativeLambdaPr", "h2dNSigmaNegativeLambdaPr", {HistType::kTH2F, {axisP, axisNSigma}}); + histos.add("h2dNSigmaPositiveK0ShortPi", "h2dNSigmaPositiveK0ShortPi", {HistType::kTH2F, {axisP, axisNSigma}}); + histos.add("h2dNSigmaNegativeK0ShortPi", "h2dNSigmaNegativeK0ShortPi", {HistType::kTH2F, {axisP, axisNSigma}}); + } + + if (calculateCascades.value > 0) { + histos.add("h2dNSigmaXiLaPi", "h2dNSigmaXiLaPi", {HistType::kTH2F, {axisP, axisNSigma}}); + histos.add("h2dNSigmaXiLaPr", "h2dNSigmaXiLaPr", {HistType::kTH2F, {axisP, axisNSigma}}); + histos.add("h2dNSigmaXiPi", "h2dNSigmaXiPi", {HistType::kTH2F, {axisP, axisNSigma}}); + histos.add("h2dNSigmaOmLaPi", "h2dNSigmaOmLaPi", {HistType::kTH2F, {axisP, axisNSigma}}); + histos.add("h2dNSigmaOmLaPr", "h2dNSigmaOmLaPr", {HistType::kTH2F, {axisP, axisNSigma}}); + histos.add("h2dNSigmaOmKa", "h2dNSigmaOmKa", {HistType::kTH2F, {axisP, axisNSigma}}); + } + } + + // delta lambda decay time + histos.add("h2dLambdaDeltaDecayTime", "h2dLambdaDeltaDecayTime", {HistType::kTH2F, {axisP, axisDeltaTime}}); + } + } + + void initCCDB(int runNumber) + { + if (mRunNumber == runNumber) { + return; + } + + // In case override, don't proceed, please - no CCDB access required + if (d_bz_input > -990) { + d_bz = d_bz_input; + o2::parameters::GRPMagField grpmag; + if (fabs(d_bz) > 1e-5) { + grpmag.setL3Current(30000.f / (d_bz / 5.0f)); + } + o2::base::Propagator::initFieldFromGRP(&grpmag); + mVtx = ccdb->getForRun(mVtxPath, runNumber); + mRunNumber = runNumber; + return; + } + + o2::parameters::GRPObject* grpo = ccdb->getForRun(grpPath, runNumber); + o2::parameters::GRPMagField* grpmag = 0x0; + if (grpo) { + o2::base::Propagator::initFieldFromGRP(grpo); + // Fetch magnetic field from ccdb for current collision + d_bz = grpo->getNominalL3Field(); + LOG(info) << "Retrieved GRP for run " << runNumber << " with magnetic field of " << d_bz << " kZG"; + } else { + grpmag = ccdb->getForRun(grpmagPath, runNumber); + if (!grpmag) { + LOG(fatal) << "Got nullptr from CCDB for path " << grpmagPath << " of object GRPMagField and " << grpPath << " of object GRPObject for run " << runNumber; + } + o2::base::Propagator::initFieldFromGRP(grpmag); + // Fetch magnetic field from ccdb for current collision + d_bz = std::lround(5.f * grpmag->getL3Current() / 30000.f); + mVtx = ccdb->getForRun(mVtxPath, runNumber); + LOG(info) << "Retrieved GRP for run " << runNumber << " with magnetic field of " << d_bz << " kZG"; + } + + // if TOF Nsigma desired + if (doNSigmas) { + nSigmaCalibObjects = ccdb->getForRun(nSigmaPath, runNumber); + if (nSigmaCalibObjects) { + LOGF(info, "loaded TList with this many objects: %i", nSigmaCalibObjects->GetEntries()); + nSigmaCalibLoaded = true; // made it thus far, mark loaded + + if (calculateV0s.value) { + hMeanPosLaPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanPosLaPi")); + hMeanPosLaPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanPosLaPr")); + hMeanNegLaPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanNegLaPi")); + hMeanNegLaPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanNegLaPr")); + hMeanPosK0Pi = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanPosK0Pi")); + hMeanNegK0Pi = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanNegK0Pi")); + + hSigmaPosLaPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaPosLaPi")); + hSigmaPosLaPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaPosLaPr")); + hSigmaNegLaPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaNegLaPi")); + hSigmaNegLaPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaNegLaPr")); + hSigmaPosK0Pi = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaPosK0Pi")); + hSigmaNegK0Pi = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaNegK0Pi")); + + if (!hMeanPosLaPi) + LOG(info) << "Problems finding mean histogram hMeanPosLaPi!"; + if (!hMeanPosLaPr) + LOG(info) << "Problems finding mean histogram hMeanPosLaPr!"; + if (!hMeanNegLaPi) + LOG(info) << "Problems finding mean histogram hMeanNegLaPi!"; + if (!hMeanNegLaPr) + LOG(info) << "Problems finding mean histogram hMeanNegLaPr!"; + if (!hMeanPosK0Pi) + LOG(info) << "Problems finding mean histogram hMeanPosK0Pi!"; + if (!hMeanNegK0Pi) + LOG(info) << "Problems finding mean histogram hMeanNegK0Pi!"; + if (!hSigmaPosK0Pi || !hSigmaNegK0Pi || !hSigmaPosLaPi || !hSigmaPosLaPr || !hSigmaNegLaPi || !hSigmaNegLaPr) { + LOG(info) << "Problems finding sigma histograms!"; + } + + if (calculateCascades.value) { + hMeanPosXiPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanPosXiPi")); + hMeanPosXiPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanPosXiPr")); + hMeanNegXiPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanNegXiPi")); + hMeanNegXiPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanNegXiPr")); + hMeanBachXiPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanBachXiPi")); + hMeanPosOmPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanPosOmPi")); + hMeanPosOmPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanPosOmPr")); + hMeanNegOmPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanNegOmPi")); + hMeanNegOmPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanNegOmPr")); + hMeanBachOmKa = reinterpret_cast(nSigmaCalibObjects->FindObject("hMeanBachOmKa")); + + hSigmaPosXiPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaPosXiPi")); + hSigmaPosXiPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaPosXiPr")); + hSigmaNegXiPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaNegXiPi")); + hSigmaNegXiPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaNegXiPr")); + hSigmaBachXiPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaBachXiPi")); + hSigmaPosOmPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaPosOmPi")); + hSigmaPosOmPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaPosOmPr")); + hSigmaNegOmPi = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaNegOmPi")); + hSigmaNegOmPr = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaNegOmPr")); + hSigmaBachOmKa = reinterpret_cast(nSigmaCalibObjects->FindObject("hSigmaBachOmKa")); + + if (!hMeanPosXiPi || !hMeanPosXiPr || !hMeanNegXiPi || !hMeanNegXiPr || !hMeanBachXiPi) + LOG(info) << "Problems finding xi mean histograms!"; + if (!hMeanPosOmPi || !hMeanPosOmPr || !hMeanNegOmPi || !hMeanNegOmPr || !hMeanBachOmKa) + LOG(info) << "Problems finding omega sigma histograms!"; + if (!hSigmaPosXiPi || !hSigmaPosXiPr || !hSigmaNegXiPi || !hSigmaNegXiPr || !hSigmaBachXiPi) + LOG(info) << "Problems finding xi sigma histograms!"; + if (!hSigmaPosOmPi || !hSigmaPosOmPr || !hSigmaNegOmPi || !hSigmaNegOmPr || !hSigmaBachOmKa) + LOG(info) << "Problems finding omega sigma histograms!"; + } + } + } + } + mRunNumber = runNumber; + } + + float velocity(float lMomentum, float lMass) + { + // Momentum p and mass m -> returns speed in centimeters per picosecond + // Useful for TOF calculations + float lA = (lMomentum / lMass) * (lMomentum / lMass); + 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) + { + // 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); + o2::track::TrackPar negTrack = o2::track::TrackPar({v0.x(), v0.y(), v0.z()}, {v0.pxneg(), v0.pyneg(), v0.pzneg()}, -1); + + float deltaTimePositiveLambdaPi = o2::aod::v0data::kNoTOFValue; + float deltaTimeNegativeLambdaPi = o2::aod::v0data::kNoTOFValue; + float deltaTimePositiveLambdaPr = o2::aod::v0data::kNoTOFValue; + float deltaTimeNegativeLambdaPr = o2::aod::v0data::kNoTOFValue; + float deltaTimePositiveK0ShortPi = o2::aod::v0data::kNoTOFValue; + float deltaTimeNegativeK0ShortPi = o2::aod::v0data::kNoTOFValue; + + float nSigmaPositiveLambdaPi = o2::aod::v0data::kNoTOFValue; + float nSigmaPositiveLambdaPr = o2::aod::v0data::kNoTOFValue; + float nSigmaNegativeLambdaPi = o2::aod::v0data::kNoTOFValue; + float nSigmaNegativeLambdaPr = o2::aod::v0data::kNoTOFValue; + float nSigmaPositiveK0ShortPi = o2::aod::v0data::kNoTOFValue; + float nSigmaNegativeK0ShortPi = 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); + + float lengthPositive = findInterceptLength(posTrack, d_bz); // FIXME: tofPosition ok? adjust? + float lengthNegative = findInterceptLength(negTrack, d_bz); // FIXME: tofPosition ok? adjust? + float timePositivePr = lengthPositive / velocityPositivePr; + float timePositivePi = lengthPositive / velocityPositivePi; + float timeNegativePr = lengthNegative / velocityNegativePr; + float timeNegativePi = lengthNegative / velocityNegativePi; + + if (pTra.hasTOF() && lengthPositive > 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() && lengthNegative > 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 (lengthPositive > 0) + posProperties = posProperties | (static_cast(1) << kLength); + if (pTra.hasTOF()) + posProperties = posProperties | (static_cast(1) << kHasTOF); + int negProperties = 0; + if (lengthNegative > 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() > 0 && lengthPositive > 0 && lengthNegative > 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); + } + + // calculate betas + + 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 betaLambda = o2::aod::cascdata::kNoTOFValue; + ; + float betaAntiLambda = o2::aod::cascdata::kNoTOFValue; + ; + float betaK0Short = o2::aod::cascdata::kNoTOFValue; + ; + + if (nTra.hasTOF() && pTra.hasTOF()) { + betaLambda = (lengthV0 / decayTimeLambda) / 0.0299792458; + betaAntiLambda = (lengthV0 / decayTimeAntiLambda) / 0.0299792458; + betaK0Short = (lengthV0 / decayTimeK0Short) / 0.0299792458; + } + + 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); + } + + float positiveP = std::hypot(v0.pxpos(), v0.pypos(), v0.pzpos()); + float negativeP = std::hypot(v0.pxneg(), v0.pyneg(), v0.pzneg()); + + if (doQA) { + if (pTra.hasTOF()) { + 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) { + histos.fill(HIST("h2dDeltaTimePositiveLambdaPr"), v0.p(), v0.eta(), deltaTimePositiveLambdaPr); + histos.fill(HIST("h2dProtonMeasuredVsExpected"), + (timeLambda + timePositivePr) / (pTra.tofSignal() - pTra.tofEvTime()), + positiveP, v0.positiveeta()); + 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) { + histos.fill(HIST("h2dDeltaTimePositiveLambdaPi"), v0.p(), v0.eta(), 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) { + histos.fill(HIST("h2dDeltaTimePositiveK0ShortPi"), v0.p(), v0.eta(), deltaTimePositiveK0ShortPi); + if (doQANSigma) + histos.fill(HIST("h2dNSigmaPositiveK0ShortPi"), v0.p(), nSigmaPositiveK0ShortPi); + } + } + } + + if (nTra.hasTOF()) { + 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) { + histos.fill(HIST("h2dDeltaTimeNegativeLambdaPi"), v0.p(), v0.eta(), deltaTimeNegativeLambdaPi); + histos.fill(HIST("h2dPionMeasuredVsExpected"), + (timeLambda + timeNegativePi) / (nTra.tofSignal() - nTra.tofEvTime()), + negativeP, v0.negativeeta()); + 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) { + histos.fill(HIST("h2dDeltaTimeNegativeLambdaPr"), v0.p(), v0.eta(), 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) { + histos.fill(HIST("h2dDeltaTimeNegativeK0ShortPi"), v0.p(), v0.eta(), deltaTimeNegativeK0ShortPi); + if (doQANSigma) + histos.fill(HIST("h2dNSigmaNegativeK0ShortPi"), v0.p(), nSigmaNegativeK0ShortPi); + } + } + } + // delta lambda decay time + histos.fill(HIST("h2dLambdaDeltaDecayTime"), v0.p(), deltaDecayTimeLambda); + } + } + + template + void processCascadeCandidate(TCollision const& collision, TCascade const& cascade, TTrack const& pTra, TTrack const& nTra, TTrack const& bTra) + { + // initialize from positions and momenta as needed + o2::track::TrackPar posTrack = o2::track::TrackPar({cascade.xlambda(), cascade.ylambda(), cascade.zlambda()}, {cascade.pxpos(), cascade.pypos(), cascade.pzpos()}, +1); + o2::track::TrackPar negTrack = o2::track::TrackPar({cascade.xlambda(), cascade.ylambda(), cascade.zlambda()}, {cascade.pxneg(), cascade.pyneg(), cascade.pzneg()}, -1); + o2::track::TrackPar bachTrack = o2::track::TrackPar({cascade.x(), cascade.y(), cascade.z()}, {cascade.pxbach(), cascade.pybach(), cascade.pzbach()}, cascade.sign()); + o2::track::TrackPar cascTrack = o2::track::TrackPar({cascade.x(), cascade.y(), cascade.z()}, {cascade.px(), cascade.py(), cascade.pz()}, cascade.sign()); + + // start calculation: calculate velocities + 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 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 daughter length to TOF intercept + 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? + + // 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()}; + bool successPropag = o2::base::Propagator::Instance()->propagateToDCA(collVtx, cascTrack, d_bz, 2.f, o2::base::Propagator::MatCorrType::USEMatCorrNONE); + float d = -1.0f, d3d = 0.0f; + float linearToPV = std::hypot(cascade.x() - collision.getX(), cascade.y() - collision.getY(), cascade.z() - collision.getZ()); + if (successPropag) { + std::array cascCloseToPVPosition; + cascTrack.getXYZGlo(cascCloseToPVPosition); + o2::math_utils::CircleXYf_t trcCircleCascade; + float sna, csa; + cascTrack.getCircleParams(d_bz, trcCircleCascade, sna, csa); + + // calculate 2D distance between two points + d = std::hypot(cascade.x() - cascCloseToPVPosition[0], cascade.y() - cascCloseToPVPosition[1]); + d3d = std::hypot(cascade.x() - cascCloseToPVPosition[0], cascade.y() - cascCloseToPVPosition[1], cascade.z() - cascCloseToPVPosition[2]); // cross-check variable + float sinThetaOverTwo = d / (2.0f * trcCircleCascade.rC); + lengthCascade = 2.0f * trcCircleCascade.rC * TMath::ASin(sinThetaOverTwo); + lengthCascade *= sqrt(1.0f + cascTrack.getTgl() * cascTrack.getTgl()); + } + + if (!successPropag) { + lengthCascade = linearToPV; // if propagation failed, use linear estimate (optional: actually do not define?) + } + + // lambda, xi and omega flight time is always defined + float lambdaFlight = lengthV0 / velocityLa; + float xiFlight = lengthCascade / velocityXi; + float omFlight = lengthCascade / velocityOm; + float posFlightPi = lengthPositive / velocityPositivePi; + float posFlightPr = lengthPositive / velocityPositivePr; + float negFlightPi = lengthNegative / velocityNegativePi; + float negFlightPr = lengthNegative / velocityNegativePr; + float bachFlightPi = lengthBachelor / velocityBachelorPi; + float bachFlightKa = lengthBachelor / velocityBachelorKa; + + // 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()) { + 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()) { + 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()) { + 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) { + // fill QA histograms for cross-checking + histos.fill(HIST("hArcDebug"), cascade.p(), lengthCascade - d3d); // for debugging purposes + + 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) { + 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 (doQANSigma) { + histos.fill(HIST("h2dNSigmaXiLaPi"), cascade.p(), nSigmaXiLaPi); + histos.fill(HIST("h2dNSigmaXiLaPr"), cascade.p(), nSigmaXiLaPr); + histos.fill(HIST("h2dNSigmaXiPi"), cascade.p(), nSigmaXiPi); + } + } + if (std::abs(cascade.mOmega() - 1.67245) < cascadeGroup.qaMassWindow && fabs(pTra.tpcNSigmaPr()) < cascadeGroup.qaTPCNSigma && fabs(nTra.tpcNSigmaPi()) < cascadeGroup.qaTPCNSigma && fabs(bTra.tpcNSigmaKa()) < cascadeGroup.qaTPCNSigma) { + 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 (doQANSigma) { + histos.fill(HIST("h2dNSigmaOmLaPi"), cascade.p(), nSigmaOmLaPi); + histos.fill(HIST("h2dNSigmaOmLaPr"), cascade.p(), nSigmaOmLaPr); + histos.fill(HIST("h2dNSigmaOmKa"), cascade.p(), nSigmaOmKa); + } + } + } 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) { + 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 (doQANSigma) { + histos.fill(HIST("h2dNSigmaXiLaPi"), cascade.p(), nSigmaXiLaPi); + histos.fill(HIST("h2dNSigmaXiLaPr"), cascade.p(), nSigmaXiLaPr); + histos.fill(HIST("h2dNSigmaXiPi"), cascade.p(), nSigmaXiPi); + } + } + if (std::abs(cascade.mOmega() - 1.67245) < cascadeGroup.qaMassWindow && fabs(pTra.tpcNSigmaPi()) < cascadeGroup.qaTPCNSigma && fabs(nTra.tpcNSigmaPr()) < cascadeGroup.qaTPCNSigma && fabs(bTra.tpcNSigmaKa()) < cascadeGroup.qaTPCNSigma) { + 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 (doQANSigma) { + histos.fill(HIST("h2dNSigmaOmLaPi"), cascade.p(), nSigmaOmLaPi); + histos.fill(HIST("h2dNSigmaOmLaPr"), cascade.p(), nSigmaOmLaPr); + histos.fill(HIST("h2dNSigmaOmKa"), cascade.p(), nSigmaOmKa); + } + } + } + } + } + } + + void processStandardData(aod::Collisions const& collisions, V0OriginalDatas const& V0s, CascOriginalDatas const& cascades, TracksWithAllExtras const&, aod::BCsWithTimestamps const& /*bcs*/) + { + // Fire up CCDB with first collision in record. If no collisions, bypass + if (useCustomRunNumber || collisions.size() < 1) { + initCCDB(manualRunNumber); + } else { + auto collision = collisions.begin(); + auto bc = collision.bc_as(); + initCCDB(bc.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_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()}); + } + + auto pTra = V0.posTrack_as(); + auto nTra = V0.negTrack_as(); + processV0Candidate(primaryVertex, V0, pTra, nTra); + } + } + + 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()}); + } + + auto pTra = cascade.posTrack_as(); + auto nTra = cascade.negTrack_as(); + auto bTra = cascade.bachelor_as(); + processCascadeCandidate(primaryVertex, cascade, pTra, nTra, bTra); + } + } + } + + void processDerivedData(soa::Join const& collisions, V0DerivedDatas const& V0s, CascDerivedDatas const& cascades, dauTracks const&) + { + // 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()}); + } + + auto pTra = V0.posTrackExtra_as(); + auto nTra = V0.negTrackExtra_as(); + processV0Candidate(primaryVertex, V0, pTra, nTra); + } + } + + 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()}); + } + + auto pTra = cascade.posTrackExtra_as(); + auto nTra = cascade.negTrackExtra_as(); + auto bTra = cascade.bachTrackExtra_as(); + processCascadeCandidate(primaryVertex, cascade, pTra, nTra, bTra); + } + } + } + + PROCESS_SWITCH(strangenesstofpid, processStandardData, "Process standard data", true); + PROCESS_SWITCH(strangenesstofpid, processDerivedData, "Process derived data", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +} diff --git a/PWGLF/Tasks/QC/strangepidqa.cxx b/PWGLF/Tasks/QC/strangepidqa.cxx index c2bd9d89774..2dd9c58b4e7 100644 --- a/PWGLF/Tasks/QC/strangepidqa.cxx +++ b/PWGLF/Tasks/QC/strangepidqa.cxx @@ -54,11 +54,13 @@ struct strangepidqa { ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "p_{T} (GeV/c)"}; ConfigurableAxis axisRadius{"axisRadius", {200, 0.0f, 100.0f}, "V0 radius (cm)"}; - AxisSpec centAxis = {100, 0.0f, 100.0f, "mult percentile"}; + ConfigurableAxis centAxis{"centAxis", {VARIABLE_WIDTH, 0.0f, 5.0f, 10.0f, 20.0f, 30.0f, 50.0f, 70.0f, 100.0f}, "FT0C centrality"}; + AxisSpec massAxisXi = {200, 1.222f, 1.422f, "Inv. Mass (GeV/c^{2})"}; AxisSpec massAxisOmega = {200, 1.572f, 1.772f, "Inv. Mass (GeV/c^{2})"}; // Invariant Mass + ConfigurableAxis axisK0ShortMass{"axisK0ShortMass", {200, 0.497f - 0.050f, 0.497f + 0.050f}, "M_{K0s} (GeV/c^{2})"}; ConfigurableAxis axisLambdaMass{"axisLambdaMass", {200, 1.08f, 1.16f}, "M_{#Lambda} (GeV/c^{2})"}; // time axes @@ -117,6 +119,9 @@ struct strangepidqa { Configurable tofNsigmaOmLaPr{"tpcNsigmaOmLaPr", 1e+5, "TOF NSigma proton <- lambda <- Omega (>10 is no cut)"}; Configurable tofNsigmaOmLaPi{"tpcNsigmaOmLaPi", 1e+5, "TOF NSigma pion <- lambda <- Omega (>10 is no cut)"}; Configurable tofNsigmaOmKa{"tpcNsigmaOmKa", 1e+5, "TOF NSigma Kaon <- Omega (>10 is no cut)"}; + + Configurable tofNsigmaCompatibility{"tofNsigmaCompatibility", 4, "compatibility check for V0s"}; + Configurable tofNsigmaCompatibilityCascades{"tofNsigmaCompatibilityCascades", 4, "compatibility check for cascades"}; //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* void init(InitContext const&) @@ -215,6 +220,21 @@ struct strangepidqa { histos.add("h2dTimeVsLengthProtonProng", "h2dTimeVsLengthProtonProng", {HistType::kTH2F, {axisLength, axisTime}}); histos.add("h2dTimeVsLengthPionProng", "h2dTimeVsLengthPionProng", {HistType::kTH2F, {axisLength, axisTime}}); + histos.add("h1dMassK0Short", "h1dMassK0Short", {HistType::kTH1F, {axisK0ShortMass}}); + histos.add("h1dMassLambda", "h1dMassLambda", {HistType::kTH1F, {axisLambdaMass}}); + histos.add("h1dMassAntiLambda", "h1dMassAntiLambda", {HistType::kTH1F, {axisLambdaMass}}); + histos.add("h1dMassCompatibleK0Short", "h1dMassCompatibleK0Short", {HistType::kTH1F, {axisK0ShortMass}}); + histos.add("h1dMassCompatibleLambda", "h1dMassCompatibleLambda", {HistType::kTH1F, {axisLambdaMass}}); + histos.add("h1dMassCompatibleAntiLambda", "h1dMassCompatibleAntiLambda", {HistType::kTH1F, {axisLambdaMass}}); + + histos.add("h3dMassK0Short", "h3dMassK0Short", {HistType::kTH3F, {centAxis, axisPt, axisK0ShortMass}}); + histos.add("h3dMassLambda", "h3dMassLambda", {HistType::kTH3F, {centAxis, axisPt, axisLambdaMass}}); + histos.add("h3dMassAntiLambda", "h3dMassAntiLambda", {HistType::kTH3F, {centAxis, axisPt, axisLambdaMass}}); + + histos.add("h3dMassCompatibleK0Short", "h3dMassCompatibleK0Short", {HistType::kTH3F, {centAxis, axisPt, axisK0ShortMass}}); + histos.add("h3dMassCompatibleLambda", "h3dMassCompatibleLambda", {HistType::kTH3F, {centAxis, axisPt, axisLambdaMass}}); + histos.add("h3dMassCompatibleAntiLambda", "h3dMassCompatibleAntiLambda", {HistType::kTH3F, {centAxis, axisPt, axisLambdaMass}}); + // --- ASSOCIATED --- // V0 Radius if (doprocessSim) { @@ -243,17 +263,29 @@ struct strangepidqa { } if (doprocessCascades) { + histos.add("h1dMassXiMinus", "h1dMassXiMinus", {HistType::kTH1F, {massAxisXi}}); + histos.add("h1dMassXiPlus", "h1dMassXiPlus", {HistType::kTH1F, {massAxisXi}}); + histos.add("h1dMassOmegaMinus", "h1dMassOmegaMinus", {HistType::kTH1F, {massAxisOmega}}); + histos.add("h1dMassOmegaPlus", "h1dMassOmegaPlus", {HistType::kTH1F, {massAxisOmega}}); + histos.add("h1dMassCompatibleXiMinus", "h1dMassCompatibleXiMinus", {HistType::kTH1F, {massAxisXi}}); + histos.add("h1dMassCompatibleXiPlus", "h1dMassCompatibleXiPlus", {HistType::kTH1F, {massAxisXi}}); + histos.add("h1dMassCompatibleOmegaMinus", "h1dMassCompatibleOmegaMinus", {HistType::kTH1F, {massAxisOmega}}); + histos.add("h1dMassCompatibleOmegaPlus", "h1dMassCompatibleOmegaPlus", {HistType::kTH1F, {massAxisOmega}}); + histos.add("h3dMassXiMinus", "h3dMassXiMinus", {HistType::kTH3F, {centAxis, axisPt, massAxisXi}}); histos.add("h3dMassXiPlus", "h3dMassXiPlus", {HistType::kTH3F, {centAxis, axisPt, massAxisXi}}); histos.add("h3dMassOmegaMinus", "h3dMassOmegaMinus", {HistType::kTH3F, {centAxis, axisPt, massAxisOmega}}); histos.add("h3dMassOmegaPlus", "h3dMassOmegaPlus", {HistType::kTH3F, {centAxis, axisPt, massAxisOmega}}); + + histos.add("h3dMassCompatibleXiMinus", "h3dMassCompatibleXiMinus", {HistType::kTH3F, {centAxis, axisPt, massAxisXi}}); + histos.add("h3dMassCompatibleXiPlus", "h3dMassCompatibleXiPlus", {HistType::kTH3F, {centAxis, axisPt, massAxisXi}}); + histos.add("h3dMassCompatibleOmegaMinus", "h3dMassCompatibleOmegaMinus", {HistType::kTH3F, {centAxis, axisPt, massAxisOmega}}); + histos.add("h3dMassCompatibleOmegaPlus", "h3dMassCompatibleOmegaPlus", {HistType::kTH3F, {centAxis, axisPt, massAxisOmega}}); } } - void processReal(soa::Join::iterator const& coll, soa::Join const& v0s) + void processReal(soa::Join::iterator const& coll, soa::Join const& v0s, soa::Join const&) { - histos.fill(HIST("hEventVertexZ"), coll.posZ()); - if (coll.centFT0C() > maxCentrality || coll.centFT0C() < minCentrality) return; @@ -262,6 +294,39 @@ struct strangepidqa { if (TMath::Abs(lambda.eta()) > 0.5) continue; + auto negExtra = lambda.negTrackExtra_as>(); + auto posExtra = lambda.posTrackExtra_as>(); + + if (TMath::Abs(posExtra.tpcNSigmaPr()) < tpcNsigmaProton && TMath::Abs(negExtra.tpcNSigmaPi()) < tpcNsigmaPion) { + // lambda case + histos.fill(HIST("h3dMassLambda"), coll.centFT0C(), lambda.pt(), lambda.mLambda()); + histos.fill(HIST("h1dMassLambda"), lambda.mLambda()); + if (lambda.tofLambdaCompatibility(tofNsigmaCompatibility.value)) { + histos.fill(HIST("h3dMassCompatibleLambda"), coll.centFT0C(), lambda.pt(), lambda.mLambda()); + histos.fill(HIST("h1dMassCompatibleLambda"), lambda.mLambda()); + } + } + + if (TMath::Abs(posExtra.tpcNSigmaPi()) < tpcNsigmaProton && TMath::Abs(negExtra.tpcNSigmaPr()) < tpcNsigmaPion) { + // lambda case + histos.fill(HIST("h3dMassAntiLambda"), coll.centFT0C(), lambda.pt(), lambda.mAntiLambda()); + histos.fill(HIST("h1dMassAntiLambda"), lambda.mAntiLambda()); + if (lambda.tofAntiLambdaCompatibility(tofNsigmaCompatibility.value)) { + histos.fill(HIST("h3dMassCompatibleAntiLambda"), coll.centFT0C(), lambda.pt(), lambda.mAntiLambda()); + histos.fill(HIST("h1dMassCompatibleAntiLambda"), lambda.mAntiLambda()); + } + } + + if (TMath::Abs(posExtra.tpcNSigmaPi()) < tpcNsigmaPion && TMath::Abs(negExtra.tpcNSigmaPr()) < tpcNsigmaPion) { + // lambda case + histos.fill(HIST("h3dMassK0Short"), coll.centFT0C(), lambda.pt(), lambda.mK0Short()); + histos.fill(HIST("h1dMassK0Short"), lambda.mK0Short()); + if (lambda.tofK0ShortCompatibility(tofNsigmaCompatibility.value)) { + histos.fill(HIST("h3dMassCompatibleK0Short"), coll.centFT0C(), lambda.pt(), lambda.mK0Short()); + histos.fill(HIST("h1dMassCompatibleK0Short"), lambda.mK0Short()); + } + } + histos.fill(HIST("h2dLambdaMassVsTOFCut"), lambda.mLambda(), TMath::Abs(lambda.posTOFDeltaTLaPr())); histos.fill(HIST("h2dLambdaMassVsTOFCutMeson"), lambda.mLambda(), TMath::Abs(lambda.negTOFDeltaTLaPi())); @@ -378,6 +443,8 @@ struct strangepidqa { void processCascades(soa::Join::iterator const& col, soa::Filtered> const& Cascades, soa::Join const&) { histos.fill(HIST("hEventCentrality"), col.centFT0C()); + if (col.centFT0C() > maxCentrality || col.centFT0C() < minCentrality) + return; for (auto& casc : Cascades) { // major selections here @@ -397,42 +464,38 @@ struct strangepidqa { if (casc.sign() < 0) { if (TMath::Abs(posExtra.tpcNSigmaPr()) < tpcNsigmaProton && TMath::Abs(negExtra.tpcNSigmaPi()) < tpcNsigmaPion && TMath::Abs(bachExtra.tpcNSigmaPi()) < tpcNsigmaBachelor) { - if (tofNsigmaXiLaPr < 100 && fabs(casc.tofNSigmaXiLaPr()) > tofNsigmaXiLaPr) - continue; - if (tofNsigmaXiLaPi < 100 && fabs(casc.tofNSigmaXiLaPi()) > tofNsigmaXiLaPi) - continue; - if (tofNsigmaXiPi < 100 && fabs(casc.tofNSigmaXiPi()) > tofNsigmaXiPi) - continue; histos.fill(HIST("h3dMassXiMinus"), col.centFT0C(), casc.pt(), casc.mXi()); + histos.fill(HIST("h1dMassXiMinus"), casc.mXi()); + if (casc.tofXiCompatibility(tofNsigmaCompatibilityCascades.value)) { + histos.fill(HIST("h3dMassCompatibleXiMinus"), col.centFT0C(), casc.pt(), casc.mXi()); + histos.fill(HIST("h1dMassCompatibleXiMinus"), casc.mXi()); + } } if (TMath::Abs(posExtra.tpcNSigmaPr()) < tpcNsigmaProton && TMath::Abs(negExtra.tpcNSigmaPi()) < tpcNsigmaPion && TMath::Abs(bachExtra.tpcNSigmaKa()) < tpcNsigmaBachelor) { - if (tofNsigmaOmLaPr < 100 && fabs(casc.tofNSigmaOmLaPr()) > tofNsigmaOmLaPr) - continue; - if (tofNsigmaOmLaPi < 100 && fabs(casc.tofNSigmaOmLaPi()) > tofNsigmaOmLaPi) - continue; - if (tofNsigmaOmKa < 100 && fabs(casc.tofNSigmaOmKa()) > tofNsigmaOmKa) - continue; histos.fill(HIST("h3dMassOmegaMinus"), col.centFT0C(), casc.pt(), casc.mOmega()); + histos.fill(HIST("h1dMassOmegaMinus"), casc.mOmega()); + if (casc.tofOmegaCompatibility(tofNsigmaCompatibilityCascades.value)) { + histos.fill(HIST("h3dMassCompatibleOmegaMinus"), col.centFT0C(), casc.pt(), casc.mOmega()); + histos.fill(HIST("h1dMassCompatibleOmegaMinus"), casc.mOmega()); + } } } else { if (TMath::Abs(posExtra.tpcNSigmaPi()) < tpcNsigmaPion && TMath::Abs(negExtra.tpcNSigmaPr()) < tpcNsigmaProton && TMath::Abs(bachExtra.tpcNSigmaPi()) < tpcNsigmaBachelor) { - if (tofNsigmaXiLaPr < 100 && fabs(casc.tofNSigmaXiLaPr()) > tofNsigmaXiLaPr) - continue; - if (tofNsigmaXiLaPi < 100 && fabs(casc.tofNSigmaXiLaPi()) > tofNsigmaXiLaPi) - continue; - if (tofNsigmaXiPi < 100 && fabs(casc.tofNSigmaXiPi()) > tofNsigmaXiPi) - continue; histos.fill(HIST("h3dMassXiPlus"), col.centFT0C(), casc.pt(), casc.mXi()); + histos.fill(HIST("h1dMassXiPlus"), casc.mXi()); + if (casc.tofXiCompatibility(tofNsigmaCompatibilityCascades.value)) { + histos.fill(HIST("h3dMassCompatibleXiPlus"), col.centFT0C(), casc.pt(), casc.mXi()); + histos.fill(HIST("h1dMassCompatibleXiPlus"), casc.mXi()); + } } if (TMath::Abs(posExtra.tpcNSigmaPi()) < tpcNsigmaPion && TMath::Abs(negExtra.tpcNSigmaPr()) < tpcNsigmaProton && TMath::Abs(bachExtra.tpcNSigmaKa()) < tpcNsigmaBachelor) { - if (tofNsigmaOmLaPr < 100 && fabs(casc.tofNSigmaOmLaPr()) > tofNsigmaOmLaPr) - continue; - if (tofNsigmaOmLaPi < 100 && fabs(casc.tofNSigmaOmLaPi()) > tofNsigmaOmLaPi) - continue; - if (tofNsigmaOmKa < 100 && fabs(casc.tofNSigmaOmKa()) > tofNsigmaOmKa) - continue; histos.fill(HIST("h3dMassOmegaPlus"), col.centFT0C(), casc.pt(), casc.mOmega()); + histos.fill(HIST("h1dMassOmegaPlus"), casc.mOmega()); + if (casc.tofOmegaCompatibility(tofNsigmaCompatibilityCascades.value)) { + histos.fill(HIST("h3dMassCompatibleOmegaPlus"), col.centFT0C(), casc.pt(), casc.mOmega()); + histos.fill(HIST("h1dMassCompatibleOmegaPlus"), casc.mOmega()); + } } } }