diff --git a/PWGLF/DataModel/LFStrangenessPIDTables.h b/PWGLF/DataModel/LFStrangenessPIDTables.h index af9dd2e44db..ab8e0ccd37b 100644 --- a/PWGLF/DataModel/LFStrangenessPIDTables.h +++ b/PWGLF/DataModel/LFStrangenessPIDTables.h @@ -43,12 +43,16 @@ namespace v0data { // ==== TOF INFORMATION === // lengths as stored in the AO2D for TOF calculations -DECLARE_SOA_COLUMN(PosTOFLength, posTOFLength, float); //! positive track length -DECLARE_SOA_COLUMN(NegTOFLength, negTOFLength, float); //! negative track length -DECLARE_SOA_COLUMN(PosTOFSignal, posTOFSignal, float); //! positive track signal -DECLARE_SOA_COLUMN(NegTOFSignal, negTOFSignal, float); //! negative track signal -DECLARE_SOA_COLUMN(PosTOFEventTime, posTOFEventTime, float); //! positive track event time -DECLARE_SOA_COLUMN(NegTOFEventTime, negTOFEventTime, float); //! negative track event time +DECLARE_SOA_COLUMN(PosTOFLengthToPV, posTOFLengthToPV, float); //! positive track length to PV +DECLARE_SOA_COLUMN(NegTOFLengthToPV, negTOFLengthToPV, float); //! negative track length to PV +DECLARE_SOA_COLUMN(PosTOFSignal, posTOFSignal, float); //! positive track signal +DECLARE_SOA_COLUMN(NegTOFSignal, negTOFSignal, float); //! negative track signal +DECLARE_SOA_COLUMN(PosTOFEventTime, posTOFEventTime, float); //! positive track event time +DECLARE_SOA_COLUMN(NegTOFEventTime, negTOFEventTime, float); //! negative track event time + +// recalculated lengths +DECLARE_SOA_COLUMN(PosTOFLength, posTOFLength, float); //! positive track length, recalculated +DECLARE_SOA_COLUMN(NegTOFLength, negTOFLength, float); //! negative track length, recalculated // delta-times DECLARE_SOA_COLUMN(PosTOFDeltaTLaPi, posTOFDeltaTLaPi, float); //! positive track TOFDeltaT from pion <- lambda expectation @@ -68,13 +72,14 @@ DECLARE_SOA_COLUMN(NegNSigmaK0Pi, negNSigmaK0Pi, float); //! positive track NSig } // namespace v0data DECLARE_SOA_TABLE(V0TOFs, "AOD", "V0TOF", // raw information table (for debug, etc) - v0data::PosTOFLength, v0data::NegTOFLength, + v0data::PosTOFLengthToPV, v0data::NegTOFLengthToPV, v0data::PosTOFSignal, v0data::NegTOFSignal, - v0data::PosTOFEventTime, v0data::NegTOFEventTime, + v0data::PosTOFEventTime, v0data::NegTOFEventTime); +DECLARE_SOA_TABLE(V0TOFPIDs, "AOD", "V0TOFPID", // processed info table (for analysis) + v0data::PosTOFLength, v0data::NegTOFLength, v0data::PosTOFDeltaTLaPi, v0data::PosTOFDeltaTLaPr, v0data::NegTOFDeltaTLaPi, v0data::NegTOFDeltaTLaPr, - v0data::PosTOFDeltaTK0Pi, v0data::NegTOFDeltaTK0Pi); -DECLARE_SOA_TABLE(V0TOFPIDs, "AOD", "V0TOFPID", // nsigma table (for analysis) + v0data::PosTOFDeltaTK0Pi, v0data::NegTOFDeltaTK0Pi, v0data::PosNSigmaLaPi, v0data::PosNSigmaLaPr, v0data::NegNSigmaLaPi, v0data::NegNSigmaLaPr, v0data::PosNSigmaK0Pi, v0data::NegNSigmaK0Pi); @@ -83,15 +88,15 @@ namespace cascdata { // ==== TOF INFORMATION === // lengths as stored in the AO2D for TOF calculations -DECLARE_SOA_COLUMN(PosTOFLength, posTOFLength, float); //! positive track length -DECLARE_SOA_COLUMN(NegTOFLength, negTOFLength, float); //! negative track length -DECLARE_SOA_COLUMN(BachTOFLength, bachTOFLength, float); //! bachelor track length -DECLARE_SOA_COLUMN(PosTOFSignal, posTOFSignal, float); //! positive track signal -DECLARE_SOA_COLUMN(NegTOFSignal, negTOFSignal, float); //! negative track signal -DECLARE_SOA_COLUMN(BachTOFSignal, bachTOFSignal, float); //! bachelor track signal -DECLARE_SOA_COLUMN(PosTOFEventTime, posTOFEventTime, float); //! positive track event time -DECLARE_SOA_COLUMN(NegTOFEventTime, negTOFEventTime, float); //! negative track event time -DECLARE_SOA_COLUMN(BachTOFEventTime, bachTOFEventTime, float); //! bachelor track event time +DECLARE_SOA_COLUMN(PosTOFLengthToPV, posTOFLengthToPV, float); //! positive track length +DECLARE_SOA_COLUMN(NegTOFLengthToPV, negTOFLengthToPV, float); //! negative track length +DECLARE_SOA_COLUMN(BachTOFLengthToPV, bachTOFLengthToPV, float); //! bachelor track length +DECLARE_SOA_COLUMN(PosTOFSignal, posTOFSignal, float); //! positive track signal +DECLARE_SOA_COLUMN(NegTOFSignal, negTOFSignal, float); //! negative track signal +DECLARE_SOA_COLUMN(BachTOFSignal, bachTOFSignal, float); //! bachelor track signal +DECLARE_SOA_COLUMN(PosTOFEventTime, posTOFEventTime, float); //! positive track event time +DECLARE_SOA_COLUMN(NegTOFEventTime, negTOFEventTime, float); //! negative track event time +DECLARE_SOA_COLUMN(BachTOFEventTime, bachTOFEventTime, float); //! bachelor track event time // delta-times DECLARE_SOA_COLUMN(PosTOFDeltaTXiPi, posTOFDeltaTXiPi, float); //! positive track TOFDeltaT from pion <- lambda <- xi expectation @@ -106,29 +111,29 @@ DECLARE_SOA_COLUMN(NegTOFDeltaTOmPr, negTOFDeltaTOmPr, float); //! negative tr DECLARE_SOA_COLUMN(BachTOFDeltaTOmPi, bachTOFDeltaTOmPi, float); //! bachelor track TOFDeltaT from pion <- omega expectation // n-sigmas -DECLARE_SOA_COLUMN(PosNSigmaXiPi, posNSigmaXiPi, float); //! positive track NSigma from pion <- lambda <- xi expectation -DECLARE_SOA_COLUMN(PosNSigmaXiPr, posNSigmaXiPr, float); //! positive track NSigma from proton <- lambda <- xi expectation -DECLARE_SOA_COLUMN(NegNSigmaXiPi, negNSigmaXiPi, float); //! negative track NSigma from pion <- lambda <- xi expectation -DECLARE_SOA_COLUMN(NegNSigmaXiPr, negNSigmaXiPr, float); //! negative track NSigma from proton <- lambda <- xi expectation -DECLARE_SOA_COLUMN(BachNSigmaXiPi, bachNSigmaXiPi, float); //! bachelor track NSigma from pion <- xi expectation -DECLARE_SOA_COLUMN(PosNSigmaOmPi, posNSigmaOmPi, float); //! positive track NSigma from pion <- lambda <- omega expectation -DECLARE_SOA_COLUMN(PosNSigmaOmPr, posNSigmaOmPr, float); //! positive track NSigma from proton <- lambda <- omega expectation -DECLARE_SOA_COLUMN(NegNSigmaOmPi, negNSigmaOmPi, float); //! negative track NSigma from pion <- lambda <- omega expectation -DECLARE_SOA_COLUMN(NegNSigmaOmPr, negNSigmaOmPr, float); //! negative track NSigma from proton <- lambda <- omega expectation -DECLARE_SOA_COLUMN(BachNSigmaOmKa, bachNSigmaOmKa, float); //! bachelor track NSigma from kaon <- omega expectation +DECLARE_SOA_COLUMN(PosNSigmaXiPi, posNSigmaXiPi, float); //! positive track NSigma from pion <- lambda <- xi expectation +DECLARE_SOA_COLUMN(PosNSigmaXiPr, posNSigmaXiPr, float); //! positive track NSigma from proton <- lambda <- xi expectation +DECLARE_SOA_COLUMN(NegNSigmaXiPi, negNSigmaXiPi, float); //! negative track NSigma from pion <- lambda <- xi expectation +DECLARE_SOA_COLUMN(NegNSigmaXiPr, negNSigmaXiPr, float); //! negative track NSigma from proton <- lambda <- xi expectation +DECLARE_SOA_COLUMN(BachNSigmaXiPi, bachNSigmaXiPi, float); //! bachelor track NSigma from pion <- xi expectation +DECLARE_SOA_COLUMN(PosNSigmaOmPi, posNSigmaOmPi, float); //! positive track NSigma from pion <- lambda <- omega expectation +DECLARE_SOA_COLUMN(PosNSigmaOmPr, posNSigmaOmPr, float); //! positive track NSigma from proton <- lambda <- omega expectation +DECLARE_SOA_COLUMN(NegNSigmaOmPi, negNSigmaOmPi, float); //! negative track NSigma from pion <- lambda <- omega expectation +DECLARE_SOA_COLUMN(NegNSigmaOmPr, negNSigmaOmPr, float); //! negative track NSigma from proton <- lambda <- omega expectation +DECLARE_SOA_COLUMN(BachNSigmaOmKa, bachNSigmaOmKa, float); //! bachelor track NSigma from kaon <- omega expectation } // namespace cascdata DECLARE_SOA_TABLE(CascTOFs, "AOD", "CascTOF", // raw information table (for debug, etc) - cascdata::PosTOFLength, cascdata::NegTOFLength, cascdata::BachTOFLength, + cascdata::PosTOFLengthToPV, cascdata::NegTOFLengthToPV, cascdata::BachTOFLengthToPV, cascdata::PosTOFSignal, cascdata::NegTOFSignal, cascdata::BachTOFSignal, - cascdata::PosTOFEventTime, cascdata::NegTOFEventTime, cascdata::BachTOFEventTime, + cascdata::PosTOFEventTime, cascdata::NegTOFEventTime, cascdata::BachTOFEventTime); +DECLARE_SOA_TABLE(CascTOFPIDs, "AOD", "CASCTOFPID", // processed information for analysis cascdata::PosTOFDeltaTXiPi, cascdata::PosTOFDeltaTXiPr, cascdata::NegTOFDeltaTXiPi, cascdata::NegTOFDeltaTXiPr, cascdata::BachTOFDeltaTXiPi, cascdata::PosTOFDeltaTOmPi, cascdata::PosTOFDeltaTOmPr, cascdata::NegTOFDeltaTOmPi, cascdata::NegTOFDeltaTOmPr, - cascdata::BachTOFDeltaTOmPi); -DECLARE_SOA_TABLE(CascTOFPIDs, "AOD", "CASCTOFPID", // nsigma table (for analysis) + cascdata::BachTOFDeltaTOmPi, cascdata::PosNSigmaXiPi, cascdata::PosNSigmaXiPr, cascdata::NegNSigmaXiPi, cascdata::NegNSigmaXiPr, cascdata::BachNSigmaXiPi, diff --git a/PWGLF/DataModel/LFStrangenessTables.h b/PWGLF/DataModel/LFStrangenessTables.h index 641e500669f..d3d1ef73813 100644 --- a/PWGLF/DataModel/LFStrangenessTables.h +++ b/PWGLF/DataModel/LFStrangenessTables.h @@ -34,6 +34,8 @@ DECLARE_SOA_TABLE(StraEPs, "AOD", "STRAEPS", //! centrality percentiles qvec::QvecFT0CRe, qvec::QvecFT0CIm, qvec::SumAmplFT0C, qvec::QvecFT0MRe, qvec::QvecFT0MIm, qvec::SumAmplFT0M, qvec::QvecFV0ARe, qvec::QvecFV0AIm, qvec::SumAmplFV0A); +DECLARE_SOA_TABLE(StraStamps, "AOD", "STRASTAMPS", //! information for ID-ing mag field if needed + bc::RunNumber, timestamp::Timestamp); using StraCollision = StraCollisions::iterator; using StraCent = StraCents::iterator; diff --git a/PWGLF/TableProducer/cascadepid.cxx b/PWGLF/TableProducer/cascadepid.cxx index 93beff60367..eb73ff9f635 100644 --- a/PWGLF/TableProducer/cascadepid.cxx +++ b/PWGLF/TableProducer/cascadepid.cxx @@ -77,8 +77,7 @@ using LabeledTracksExtra = soa::Join; struct cascadepid { // TOF pid for strangeness (recalculated with topology) - Produces casctof; // raw table for checks - Produces casctofpid; // table with Nsigmas + Produces casctofpids; // table with Nsigmas Service ccdb; @@ -246,14 +245,11 @@ struct cascadepid { auto negTrack = cascade.negTrack_as(); // FIXME: TOF calculation: under construction, to follow - - if (fillRawPID) { - casctof(posTrack.length(), negTrack.length(), bachTrack.length(), - posTrack.tofSignal(), negTrack.tofSignal(), bachTrack.tofSignal(), - posTrack.tofEvTime(), negTrack.tofEvTime(), bachTrack.tofEvTime(), - 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, - 0.0f, 0.0f, 0.0f, 0.0f, 0.0f); - } + casctofpids(0.0f, 0.0f, + posTrack.length(), negTrack.length(), bachTrack.length(), + 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, + 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, + 0.0f, 0.0f, 0.0f, 0.0f, 0.0f); } } } diff --git a/PWGLF/TableProducer/lambdakzeropid.cxx b/PWGLF/TableProducer/lambdakzeropid.cxx index 0d353634411..626b1adc5c8 100644 --- a/PWGLF/TableProducer/lambdakzeropid.cxx +++ b/PWGLF/TableProducer/lambdakzeropid.cxx @@ -77,7 +77,6 @@ using LabeledTracksExtra = soa::Join; struct lambdakzeropid { // TOF pid for strangeness (recalculated with topology) - Produces v0tof; // raw table for checks Produces v0tofpid; // table with Nsigmas Service ccdb; @@ -89,9 +88,8 @@ struct lambdakzeropid { // Operation and minimisation criteria Configurable d_bz_input{"d_bz", -999, "bz field, -999 is automatic"}; - Configurable tofPosition{"tofPosition", 370, "TOF position for tests"}; + Configurable tofPosition{"tofPosition", 377.934f, "TOF effective (inscribed) radius"}; Configurable checkTPCCompatibility{"checkTPCCompatibility", true, "check compatibility with dE/dx in QA plots"}; - Configurable fillRawPID{"fillRawPID", true, "fill raw PID tables for debug/x-check"}; // CCDB options Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; @@ -110,71 +108,132 @@ struct lambdakzeropid { 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 radius + /// 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 radius the radius of the layer you're calculating the length to + /// \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 trackLength(o2::track::TrackParCov track, float radius, float magneticField) + float trackLengthToSegment(o2::track::TrackParCov track, float x1, float y1, float x2, float y2, float magneticField) { // don't make use of the track parametrization - float length = -100; + float length = -104; + // causality protection + std::array mom; + track.getPxPyPzGlo(mom); + // get start point + std::array startPoint; + track.getXYZGlo(startPoint); + + 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); - // distance between circle centers (one circle is at origin -> easy) - float centerDistance = std::hypot(trcCircle.xC, trcCircle.yC); - - // condition of circles touching - if not satisfied returned length will be -100 - if (centerDistance < trcCircle.rC + radius && centerDistance > fabs(trcCircle.rC - radius)) { - length = 0.0f; - - // base radical direction - float ux = trcCircle.xC / centerDistance; - float uy = trcCircle.yC / centerDistance; - // calculate perpendicular vector (normalized) for +/- displacement - float vx = -uy; - float vy = +ux; - // calculate coordinate for radical line - float radical = (centerDistance * centerDistance - trcCircle.rC * trcCircle.rC + radius * radius) / (2.0f * centerDistance); - // calculate absolute displacement from center-to-center axis - float displace = (0.5f / centerDistance) * TMath::Sqrt( - (-centerDistance + trcCircle.rC - radius) * - (-centerDistance - trcCircle.rC + radius) * - (-centerDistance + trcCircle.rC + radius) * - (centerDistance + trcCircle.rC + radius)); - - // possible intercept points of track and TOF layer in 2D plane - float point1[2] = {radical * ux + displace * vx, radical * uy + displace * vy}; - float point2[2] = {radical * ux - displace * vx, radical * uy - displace * vy}; - - // decide on correct intercept point - std::array mom; - track.getPxPyPzGlo(mom); - float scalarProduct1 = point1[0] * mom[0] + point1[1] * mom[1]; - float scalarProduct2 = point2[0] * mom[0] + point2[1] * mom[1]; - - // get start point - std::array startPoint; - track.getXYZGlo(startPoint); - - float cosAngle = -1000, modulus = -1000; - - if (scalarProduct1 > scalarProduct2) { - modulus = std::hypot(point1[0] - trcCircle.xC, point1[1] - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC); - cosAngle = (point1[0] - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (point1[1] - trcCircle.yC) * (startPoint[0] - trcCircle.yC); - } else { - modulus = std::hypot(point2[0] - trcCircle.xC, point2[1] - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC); - cosAngle = (point2[0] - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (point2[1] - trcCircle.yC) * (startPoint[0] - trcCircle.yC); - } - cosAngle /= modulus; - length = trcCircle.rC * TMath::ACos(cosAngle); - length *= sqrt(1.0f + track.getTgl() * track.getTgl()); + // 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; + + float halfPerimeter = TMath::Pi() * trcCircle.rC; // perimeter check. Length should not pass this ever + + if (scalarCheck1 > 0.0f && scalarCheck1 < segmentModulus && length1 < halfPerimeter && scalarMomentumCheck1 > 0.0f) { + length = length1; + // X = interceptX1; Y = interceptY1; + } + if (scalarCheck2 > 0.0f && scalarCheck2 < segmentModulus && length2 < halfPerimeter && 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::TrackParCov track, float magneticField) + { + 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 length = trackLengthToSegment(track, x1, y1, x2, y2, magneticField); + if (length > 0) + return length; + } + return -100; // not reached / not found + } + void init(InitContext& context) { mRunNumber = 0; @@ -282,8 +341,8 @@ struct lambdakzeropid { posTrack.propagateTo(v0.posX(), d_bz); negTrack.propagateTo(v0.negX(), d_bz); - float lengthPositive = trackLength(posTrack, tofPosition, d_bz); // FIXME: tofPosition ok? adjust? - float lengthNegative = trackLength(negTrack, tofPosition, d_bz); // FIXME: tofPosition ok? adjust? + 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; @@ -296,14 +355,11 @@ struct lambdakzeropid { deltaTimePositiveK0ShortPi = (posTrackRow.tofSignal() - posTrackRow.tofEvTime()) - (timeK0Short + timeNegativePi); deltaTimeNegativeK0ShortPi = (negTrackRow.tofSignal() - negTrackRow.tofEvTime()) - (timeK0Short + timeNegativePi); - if (fillRawPID) { - v0tof(posTrackRow.length(), negTrackRow.length(), - posTrackRow.tofSignal(), negTrackRow.tofSignal(), - posTrackRow.tofEvTime(), negTrackRow.tofEvTime(), - deltaTimePositiveLambdaPi, deltaTimePositiveLambdaPr, - deltaTimeNegativeLambdaPi, deltaTimeNegativeLambdaPr, - deltaTimePositiveK0ShortPi, deltaTimeNegativeK0ShortPi); - } + v0tofpid(lengthPositive, lengthNegative, + deltaTimePositiveLambdaPi, deltaTimePositiveLambdaPr, + deltaTimeNegativeLambdaPi, deltaTimeNegativeLambdaPr, + deltaTimePositiveK0ShortPi, deltaTimeNegativeK0ShortPi, + 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f); // FIXME auto originalV0 = v0.v0_as(); // this could look confusing, so: // the first v0 is the v0data row; the getter de-references the v0 (stored indices) row diff --git a/PWGLF/TableProducer/strangederivedbuilder.cxx b/PWGLF/TableProducer/strangederivedbuilder.cxx index 05196d577be..ef9dc048c5d 100644 --- a/PWGLF/TableProducer/strangederivedbuilder.cxx +++ b/PWGLF/TableProducer/strangederivedbuilder.cxx @@ -52,6 +52,7 @@ using namespace o2::framework::expressions; using std::array; using TracksWithExtra = soa::Join; +using FullTracksExtIUTOF = soa::Join; // simple checkers #define bitset(var, nbit) ((var) |= (1 << (nbit))) @@ -62,6 +63,7 @@ struct strangederivedbuilder { // fundamental building blocks of derived data Produces strangeColl; // characterises collisions Produces strangeCents; // characterises collisions / centrality + Produces strangeStamps; // provides timestamps, run numbers Produces v0collref; // references collisions from V0s Produces casccollref; // references collisions from cascades Produces kfcasccollref; // references collisions from KF cascades @@ -84,9 +86,14 @@ struct strangederivedbuilder { //__________________________________________________ // mother information - Produces v0mothers; // V0 mother references - Produces cascmothers; // casc mother references - Produces motherMCParts; // mc particles for mothers + Produces v0mothers; // V0 mother references + Produces cascmothers; // casc mother references + Produces motherMCParts; // mc particles for mothers + + //__________________________________________________ + // raw TOF PID for posterior use if requested + Produces v0tofs; // V0 part + Produces casctofs; // cascade part // histogram registry for bookkeeping HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -119,8 +126,11 @@ struct strangederivedbuilder { Preslice KFCascperCollision = o2::aod::cascdata::collisionId; Preslice TraCascperCollision = o2::aod::cascdata::collisionId; + int64_t currentCollIdx; + void init(InitContext& context) { + currentCollIdx = -1; // setup map for fast checking if enabled static_for<0, nSpecies - 1>([&](auto i) { constexpr int index = i.value; @@ -135,30 +145,27 @@ struct strangederivedbuilder { histos.add(Form("hGen%s", particleNames[i].data()), Form("hGen%s", particleNames[i].data()), kTH1D, {axisPt}); } - void processCollisionsV0sOnly(soa::Join const& collisions, aod::V0Datas const& V0s) + void processCollisionsV0sOnly(soa::Join const& collisions, aod::V0Datas const& V0s, aod::BCsWithTimestamps const&) { - int currentCollIdx = -1; for (const auto& collision : collisions) { const uint64_t collIdx = collision.globalIndex(); auto V0Table_thisColl = V0s.sliceBy(V0perCollision, collIdx); bool strange = V0Table_thisColl.size() > 0; // casc table sliced if (strange || fillEmptyCollisions) { - if (currentCollIdx != collIdx) { - strangeColl(collision.posX(), collision.posY(), collision.posZ()); - strangeCents(collision.centFT0M(), collision.centFT0A(), - collision.centFT0C(), collision.centFV0A()); - currentCollIdx = collIdx; - } + strangeColl(collision.posX(), collision.posY(), collision.posZ()); + strangeCents(collision.centFT0M(), collision.centFT0A(), + collision.centFT0C(), collision.centFV0A()); + auto bc = collision.bc_as(); + strangeStamps(bc.timestamp(), bc.runNumber()); } for (int i = 0; i < V0Table_thisColl.size(); i++) v0collref(strangeColl.lastIndex()); } } - void processCollisions(soa::Join const& collisions, aod::V0Datas const& V0s, aod::CascDatas const& Cascades, aod::KFCascDatas const& KFCascades, aod::TraCascDatas const& TraCascades) + void processCollisions(soa::Join const& collisions, aod::V0Datas const& V0s, aod::CascDatas const& Cascades, aod::KFCascDatas const& KFCascades, aod::TraCascDatas const& TraCascades, aod::BCsWithTimestamps const&) { - int currentCollIdx = -1; for (const auto& collision : collisions) { const uint64_t collIdx = collision.globalIndex(); auto V0Table_thisColl = V0s.sliceBy(V0perCollision, collIdx); @@ -171,12 +178,11 @@ struct strangederivedbuilder { TraCascTable_thisColl.size() > 0; // casc table sliced if (strange || fillEmptyCollisions) { - if (currentCollIdx != collIdx) { - strangeColl(collision.posX(), collision.posY(), collision.posZ()); - strangeCents(collision.centFT0M(), collision.centFT0A(), - collision.centFT0C(), collision.centFV0A()); - currentCollIdx = collIdx; - } + strangeColl(collision.posX(), collision.posY(), collision.posZ()); + strangeCents(collision.centFT0M(), collision.centFT0A(), + collision.centFT0C(), collision.centFV0A()); + auto bc = collision.bc_as(); + strangeStamps(bc.timestamp(), bc.runNumber()); } for (int i = 0; i < V0Table_thisColl.size(); i++) v0collref(strangeColl.lastIndex()); @@ -418,6 +424,28 @@ struct strangederivedbuilder { } } + void processProduceV0TOFs(aod::Collision const& collision, aod::V0Datas const& V0s, FullTracksExtIUTOF const&) + { + for (auto const& v0 : V0s) { + auto const& posTrackRow = v0.posTrack_as(); + auto const& negTrackRow = v0.negTrack_as(); + v0tofs(posTrackRow.length(), negTrackRow.length(), + posTrackRow.tofSignal(), negTrackRow.tofSignal(), + posTrackRow.tofEvTime(), negTrackRow.tofEvTime()); + } + } + void processProduceCascTOFs(aod::Collision const& collision, aod::CascDatas const& Cascades, FullTracksExtIUTOF const&) + { + for (auto const& cascade : Cascades) { + auto const& posTrackRow = cascade.posTrack_as(); + auto const& negTrackRow = cascade.negTrack_as(); + auto const& bachTrackRow = cascade.bachelor_as(); + casctofs(posTrackRow.length(), negTrackRow.length(), bachTrackRow.length(), + posTrackRow.tofSignal(), negTrackRow.tofSignal(), bachTrackRow.tofSignal(), + posTrackRow.tofEvTime(), negTrackRow.tofEvTime(), bachTrackRow.tofEvTime()); + } + } + PROCESS_SWITCH(strangederivedbuilder, processCollisionsV0sOnly, "Produce collisions (V0s only)", true); PROCESS_SWITCH(strangederivedbuilder, processCollisions, "Produce collisions (V0s + casc)", true); PROCESS_SWITCH(strangederivedbuilder, processTrackExtrasV0sOnly, "Produce track extra information (V0s only)", true); @@ -426,6 +454,8 @@ struct strangederivedbuilder { PROCESS_SWITCH(strangederivedbuilder, processCascadeInterlinkTracked, "Produce tables interconnecting cascades", false); PROCESS_SWITCH(strangederivedbuilder, processCascadeInterlinkKF, "Produce tables interconnecting cascades", false); PROCESS_SWITCH(strangederivedbuilder, processSimulation, "Produce simulated information", true); + PROCESS_SWITCH(strangederivedbuilder, processProduceV0TOFs, "Produce V0TOFs table", true); + PROCESS_SWITCH(strangederivedbuilder, processProduceCascTOFs, "Produce CascTOFs table", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)