diff --git a/PWGLF/DataModel/LFStrangenessPIDTables.h b/PWGLF/DataModel/LFStrangenessPIDTables.h index ab8e0ccd37b..fdef30d1d1f 100644 --- a/PWGLF/DataModel/LFStrangenessPIDTables.h +++ b/PWGLF/DataModel/LFStrangenessPIDTables.h @@ -62,13 +62,31 @@ DECLARE_SOA_COLUMN(NegTOFDeltaTLaPr, negTOFDeltaTLaPr, float); //! negative trac DECLARE_SOA_COLUMN(PosTOFDeltaTK0Pi, posTOFDeltaTK0Pi, float); //! positive track TOFDeltaT from pion <- k0short expectation DECLARE_SOA_COLUMN(NegTOFDeltaTK0Pi, negTOFDeltaTK0Pi, float); //! positive track TOFDeltaT from pion <- k0short expectation -// n-sigmas +// delta-decay-times (event-time-independent) +DECLARE_SOA_COLUMN(DeltaDecayTimeLambda, deltaDecayTimeLambda, float); //! delta-decay time estimate from proton/pion from Lambda +DECLARE_SOA_COLUMN(DeltaDecayTimeAntiLambda, deltaDecayTimeAntiLambda, float); //! delta-decay time estimate from pion/proton from AntiLambda +DECLARE_SOA_COLUMN(DeltaDecayTimeK0Short, deltaDecayTimeK0Short, float); //! delta-decay time estimate from pion/pion from K0Short + +// n-sigmas - unused as of now DECLARE_SOA_COLUMN(PosNSigmaLaPi, posNSigmaLaPi, float); //! positive track NSigma from pion <- lambda expectation DECLARE_SOA_COLUMN(PosNSigmaLaPr, posNSigmaLaPr, float); //! positive track NSigma from proton <- lambda expectation DECLARE_SOA_COLUMN(NegNSigmaLaPi, negNSigmaLaPi, float); //! negative track NSigma from pion <- lambda expectation DECLARE_SOA_COLUMN(NegNSigmaLaPr, negNSigmaLaPr, float); //! negative track NSigma from proton <- lambda expectation DECLARE_SOA_COLUMN(PosNSigmaK0Pi, posNSigmaK0Pi, float); //! positive track NSigma from pion <- k0short expectation DECLARE_SOA_COLUMN(NegNSigmaK0Pi, negNSigmaK0Pi, float); //! positive track NSigma from pion <- k0short expectation + +// beta values +DECLARE_SOA_COLUMN(TofBetaLambda, tofBetaLambda, float); //! beta value with Lambda hypothesis +DECLARE_SOA_COLUMN(TofBetaAntiLambda, tofBetaAntiLambda, float); //! beta value with AntiLambda hypothesis +DECLARE_SOA_COLUMN(TofBetaK0Short, tofBetaK0Short, float); //! beta value with K0Short hypothesis + +// debug quantities +DECLARE_SOA_COLUMN(V0LifetimeLambda, v0LifetimeLambda, float); //! lifetime of V0 assuming lambda mass (ps) +DECLARE_SOA_COLUMN(V0LifetimeK0Short, v0LifetimeK0Short, float); //! lifetime of V0 assuming K0 mass (ps) +DECLARE_SOA_COLUMN(PosLifetimePr, posLifetimePr, float); //! lifetime (to TOF) of pos prong assuming proton mass (ps) +DECLARE_SOA_COLUMN(PosLifetimePi, posLifetimePi, float); //! lifetime (to TOF) of pos prong assuming pion mass (ps) +DECLARE_SOA_COLUMN(NegLifetimePr, negLifetimePr, float); //! lifetime (to TOF) of neg prong assuming proton mass (ps) +DECLARE_SOA_COLUMN(NegLifetimePi, negLifetimePi, float); //! lifetime (to TOF) of neg prong assuming pion mass (ps) } // namespace v0data DECLARE_SOA_TABLE(V0TOFs, "AOD", "V0TOF", // raw information table (for debug, etc) @@ -80,9 +98,17 @@ DECLARE_SOA_TABLE(V0TOFPIDs, "AOD", "V0TOFPID", // processed info table (for ana v0data::PosTOFDeltaTLaPi, v0data::PosTOFDeltaTLaPr, v0data::NegTOFDeltaTLaPi, v0data::NegTOFDeltaTLaPr, v0data::PosTOFDeltaTK0Pi, v0data::NegTOFDeltaTK0Pi, - v0data::PosNSigmaLaPi, v0data::PosNSigmaLaPr, - v0data::NegNSigmaLaPi, v0data::NegNSigmaLaPr, - v0data::PosNSigmaK0Pi, v0data::NegNSigmaK0Pi); + v0data::DeltaDecayTimeLambda, + v0data::DeltaDecayTimeAntiLambda, + v0data::DeltaDecayTimeK0Short); + +DECLARE_SOA_TABLE(V0TOFDebugs, "AOD", "V0TOFDEBUG", // table with intermediate information solely for debugging + v0data::V0LifetimeLambda, v0data::V0LifetimeK0Short, + v0data::PosLifetimePr, v0data::PosLifetimePi, + v0data::NegLifetimePr, v0data::NegLifetimePi); + +DECLARE_SOA_TABLE(V0TOFBetas, "AOD", "V0TOFBETA", // processed info table (for analysis) + v0data::TofBetaLambda, v0data::TofBetaAntiLambda, v0data::TofBetaK0Short); namespace cascdata { diff --git a/PWGLF/DataModel/LFStrangenessTables.h b/PWGLF/DataModel/LFStrangenessTables.h index fdd0e6e993c..9d962372322 100644 --- a/PWGLF/DataModel/LFStrangenessTables.h +++ b/PWGLF/DataModel/LFStrangenessTables.h @@ -108,7 +108,7 @@ DECLARE_SOA_COLUMN(PDGCode, pdgCode, int); //! pdg code DECLARE_SOA_COLUMN(IsPhysicalPrimary, isPhysicalPrimary, bool); //! primary criterion } // namespace motherParticle -DECLARE_SOA_TABLE(MotherMCParts, "AOD", "MOTHERMCPARTS", //! mother MC information, abbreviated name due to size limit +DECLARE_SOA_TABLE(MotherMCParts, "AOD", "MOTHERMCPART", //! mother MC information, abbreviated name due to size limit motherParticle::Px, motherParticle::Py, motherParticle::Pz, motherParticle::PDGCode, motherParticle::IsPhysicalPrimary); diff --git a/PWGLF/TableProducer/lambdakzeropid.cxx b/PWGLF/TableProducer/lambdakzeropid.cxx index 2dd4ecb79c7..bdda2efffab 100644 --- a/PWGLF/TableProducer/lambdakzeropid.cxx +++ b/PWGLF/TableProducer/lambdakzeropid.cxx @@ -81,6 +81,8 @@ using V0FullCores = 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 Service ccdb; @@ -92,7 +94,7 @@ struct lambdakzeropid { // 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 checkTPCCompatibility{"checkTPCCompatibility", true, "check compatibility with dE/dx in QA plots"}; + Configurable doQA{"doQA", true, "create QA histos"}; // CCDB options Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; @@ -101,10 +103,9 @@ struct lambdakzeropid { Configurable lutPath{"lutPath", "GLO/Param/MatLUT", "Path of the Lut parametrization"}; Configurable geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; - ConfigurableAxis axisPtQA{"axisPtQA", {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}, "pt axis for QA histograms"}; ConfigurableAxis axisDeltaTime{"axisDeltaTime", {2000, -1000.0f, +1000.0f}, "delta-time (ps)"}; - ConfigurableAxis axisK0ShortMass{"axisK0ShortMass", {200, 0.400f, 0.600f}, "Inv. Mass (GeV/c^{2})"}; - ConfigurableAxis axisLambdaMass{"axisLambdaMass", {200, 1.01f, 1.21f}, "Inv. Mass (GeV/c^{2})"}; + 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)"}; int mRunNumber; float d_bz; @@ -132,8 +133,9 @@ struct lambdakzeropid { std::array startPoint; track.getXYZGlo(startPoint); - if (((x1 + x2) * mom[0] + (y1 + y2) * mom[1]) < 0.0f) - return -101; + // 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; @@ -203,13 +205,11 @@ struct lambdakzeropid { 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) { + if (scalarCheck1 > 0.0f && scalarCheck1 < segmentModulus && scalarMomentumCheck1 > 0.0f) { length = length1; // X = interceptX1; Y = interceptY1; } - if (scalarCheck2 > 0.0f && scalarCheck2 < segmentModulus && length2 < halfPerimeter && scalarMomentumCheck2 > 0.0f) { + if (scalarCheck2 > 0.0f && scalarCheck2 < segmentModulus && scalarMomentumCheck2 > 0.0f) { length = length2; // X = interceptX2; Y = interceptY2; } @@ -221,6 +221,7 @@ struct lambdakzeropid { /// \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(); @@ -230,11 +231,13 @@ struct lambdakzeropid { 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; + float thisLength = trackLengthToSegment(track, x1, y1, x2, y2, magneticField); + if (thisLength < length && thisLength > 0) + length = thisLength; } - return -100; // not reached / not found + if (length > 1e+5) + length = -100; // force negative to avoid misunderstandings + return length; } void init(InitContext& context) @@ -251,6 +254,15 @@ struct lambdakzeropid { // 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}}); + + // delta lambda decay time + histos.add("h2dLambdaDeltaDecayTime", "h2dLambdaDeltaDecayTime", {HistType::kTH2F, {axisPt, axisDeltaTime}}); + } } void initCCDB(soa::Join::iterator const& collision) @@ -341,18 +353,63 @@ struct lambdakzeropid { float timeNegativePr = lengthNegative / velocityNegativePr; float timeNegativePi = lengthNegative / velocityNegativePi; - deltaTimePositiveLambdaPr = (v0.posTOFSignal() - v0.posTOFEventTime()) - (timeLambda + timePositivePr); - deltaTimePositiveLambdaPi = (v0.posTOFSignal() - v0.posTOFEventTime()) - (timeLambda + timePositivePi); - deltaTimeNegativeLambdaPr = (v0.negTOFSignal() - v0.negTOFEventTime()) - (timeLambda + timeNegativePr); - deltaTimeNegativeLambdaPi = (v0.negTOFSignal() - v0.negTOFEventTime()) - (timeLambda + timeNegativePi); - deltaTimePositiveK0ShortPi = (v0.posTOFSignal() - v0.posTOFEventTime()) - (timeK0Short + timeNegativePi); - deltaTimeNegativeK0ShortPi = (v0.negTOFSignal() - v0.negTOFEventTime()) - (timeK0Short + timeNegativePi); + if (v0.posTOFSignal() > 0 && v0.posTOFEventTime() > 0 && lengthPositive > 0) { + deltaTimePositiveLambdaPr = (v0.posTOFSignal() - v0.posTOFEventTime()) - (timeLambda + timePositivePr); + deltaTimePositiveLambdaPi = (v0.posTOFSignal() - v0.posTOFEventTime()) - (timeLambda + timePositivePi); + deltaTimePositiveK0ShortPi = (v0.posTOFSignal() - v0.posTOFEventTime()) - (timeK0Short + timeNegativePi); + } + if (v0.negTOFSignal() > 0 && v0.negTOFEventTime() > 0 && lengthNegative > 0) { + deltaTimeNegativeLambdaPr = (v0.negTOFSignal() - v0.negTOFEventTime()) - (timeLambda + timeNegativePr); + deltaTimeNegativeLambdaPi = (v0.negTOFSignal() - v0.negTOFEventTime()) - (timeLambda + timeNegativePi); + deltaTimeNegativeK0ShortPi = (v0.negTOFSignal() - v0.negTOFEventTime()) - (timeK0Short + timeNegativePi); + } + float deltaDecayTimeLambda = -10e+4; + float deltaDecayTimeAntiLambda = -10e+4; + float deltaDecayTimeK0Short = -10e+4; + if (v0.posTOFSignal() > 0 && v0.negTOFSignal() > 0 && lengthPositive > 0 && lengthNegative > 0) { // does not depend on event time + deltaDecayTimeLambda = (v0.posTOFSignal() - timePositivePr) - (v0.negTOFSignal() - timeNegativePi); + deltaDecayTimeAntiLambda = (v0.posTOFSignal() - timePositivePi) - (v0.negTOFSignal() - timeNegativePr); + deltaDecayTimeK0Short = (v0.posTOFSignal() - timePositivePi) - (v0.negTOFSignal() - timeNegativePi); + } + + // calculate betas + + float evTimeMean = 0.5f * (v0.posTOFEventTime() + v0.negTOFEventTime()); + float decayTimeLambda = 0.5f * ((v0.posTOFSignal() - timePositivePr) + (v0.negTOFSignal() - timeNegativePi)) - evTimeMean; + float decayTimeAntiLambda = 0.5f * ((v0.posTOFSignal() - timePositivePi) + (v0.negTOFSignal() - timeNegativePr)) - evTimeMean; + float decayTimeK0Short = 0.5f * ((v0.posTOFSignal() - timePositivePi) + (v0.negTOFSignal() - timeNegativePi)) - evTimeMean; + + float betaLambda = -1e+6; + float betaAntiLambda = -1e+6; + float betaK0Short = -1e+6; + + if (v0.posTOFSignal() > 0 && v0.negTOFSignal() > 0 && v0.posTOFEventTime() > 0 && v0.negTOFEventTime() > 0) { + betaLambda = (lengthV0 / decayTimeLambda) / 0.0299792458; + betaAntiLambda = (lengthV0 / decayTimeAntiLambda) / 0.0299792458; + betaK0Short = (lengthV0 / decayTimeK0Short) / 0.0299792458; + } v0tofpid(lengthPositive, lengthNegative, deltaTimePositiveLambdaPi, deltaTimePositiveLambdaPr, deltaTimeNegativeLambdaPi, deltaTimeNegativeLambdaPr, deltaTimePositiveK0ShortPi, deltaTimeNegativeK0ShortPi, - 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f); // FIXME + deltaDecayTimeLambda, deltaDecayTimeAntiLambda, deltaDecayTimeK0Short); + v0tofbeta(betaLambda, betaAntiLambda, betaK0Short); + v0tofdebugs(timeLambda, timeK0Short, timePositivePr, timePositivePi, timeNegativePr, timeNegativePi); + + if (doQA) { + if (v0.posTOFSignal() > 0 && v0.posTOFEventTime() > 0) + histos.fill(HIST("h2dProtonMeasuredVsExpected"), + (timeLambda + timePositivePr), + (v0.posTOFSignal() - v0.posTOFEventTime())); + if (v0.negTOFSignal() > 0 && v0.negTOFEventTime() > 0) + histos.fill(HIST("h2dPionMeasuredVsExpected"), + (timeLambda + timeNegativePi), + (v0.negTOFSignal() - v0.negTOFEventTime())); + + // delta lambda decay time + histos.fill(HIST("h2dLambdaDeltaDecayTime"), v0.pt(), deltaDecayTimeLambda); + } } } } diff --git a/PWGLF/Tasks/QC/strangepidqa.cxx b/PWGLF/Tasks/QC/strangepidqa.cxx index a27aa7df1d3..5dc7173b1a3 100644 --- a/PWGLF/Tasks/QC/strangepidqa.cxx +++ b/PWGLF/Tasks/QC/strangepidqa.cxx @@ -57,23 +57,100 @@ struct strangepidqa { // Invariant Mass ConfigurableAxis axisLambdaMass{"axisLambdaMass", {200, 1.08f, 1.16f}, "M_{#Lambda} (GeV/c^{2})"}; - // Delta time + // time axes ConfigurableAxis axisDeltaTime{"axisDeltaTime", {200, -1000.0f, +1000.0f}, "#Delta time (ps)"}; + ConfigurableAxis axisTime{"axisTime", {200, 0.0f, +20000.0f}, "T (ps)"}; + ConfigurableAxis axisBeta{"axisBeta", {1200, 0.0f, +1.2f}, "#Beta"}; // Length axis ConfigurableAxis axisLength{"axisLength", {600, 0.0f, +600.0f}, "track Length (cm)"}; + // TOF cut axis + ConfigurableAxis axisTOFCut{"axisTOFCut", {100, 0.0f, +10000.0f}, "TOF compat. cut (ps)"}; + + // TOF selection matters + Configurable requireTOFsignalPion{"requireTOFsignalPion", true, "require that pion prongs have TOF"}; + Configurable requireTOFsignalProton{"requireTOFsignalProton", true, "require that proton prongs have TOF"}; + Configurable requireTOFEventTimePion{"requireTOFEventTimePion", true, "require that pion prongs have TOF event time"}; + Configurable requireTOFEventTimeProton{"requireTOFEventTimeProton", true, "require that proton prongs have TOF event time"}; + + Configurable maxDeltaTimeProton{"maxDeltaTimeProton", 1e+9, "check maximum allowed time"}; + Configurable maxDeltaTimePion{"maxDeltaTimePion", 1e+9, "check maximum allowed time"}; + Configurable maxDeltaTimeDecay{"maxDeltaTimeDecay", 1e+9, "check maximum allowed delta-time-decay"}; + + Configurable minCentrality{"minCentrality", 60, "max value of centrality allowed"}; + Configurable maxCentrality{"maxCentrality", 100, "max value of centrality allowed"}; + + Configurable minV0Radius{"minV0Radius", 1.5, "min radius"}; + Configurable minCosPA{"minCosPA", .98, "min cosPA"}; + + Configurable minPtV0{"minPtV0", 1.0, "min pT for integrated mass histograms"}; + Configurable maxPtV0{"maxPtV0", 3.0, "max pT for integrated mass histograms"}; + void init(InitContext const&) { // Event counter histos.add("hEventVertexZ", "hEventVertexZ", kTH1F, {vertexZ}); + // Presence of negative and positive signals and event time + auto hPositiveStatus = histos.add("hPositiveStatus", "hPositiveStatus", kTH1F, {{4, -0.5f, 3.5f}}); + auto hNegativeStatus = histos.add("hNegativeStatus", "hNegativeStatus", kTH1F, {{4, -0.5f, 3.5f}}); + auto hPositiveStatusReachedTOF = histos.add("hPositiveStatusReachedTOF", "hPositiveStatusReachedTOF", kTH1F, {{4, -0.5f, 3.5f}}); + auto hNegativeStatusReachedTOF = histos.add("hNegativeStatusReachedTOF", "hNegativeStatusReachedTOF", kTH1F, {{4, -0.5f, 3.5f}}); + hPositiveStatus->GetXaxis()->SetBinLabel(1, "All positive"); + hPositiveStatus->GetXaxis()->SetBinLabel(2, "Has only TOF sig"); + hPositiveStatus->GetXaxis()->SetBinLabel(3, "Has only TOF ev time"); + hPositiveStatus->GetXaxis()->SetBinLabel(4, "Has full time info"); + hNegativeStatus->GetXaxis()->SetBinLabel(1, "All negative"); + hNegativeStatus->GetXaxis()->SetBinLabel(2, "Has only TOF sig"); + hNegativeStatus->GetXaxis()->SetBinLabel(3, "Has only TOF ev time"); + hNegativeStatus->GetXaxis()->SetBinLabel(4, "Has full time info"); + + hPositiveStatusReachedTOF->GetXaxis()->SetBinLabel(1, "All positive"); + hPositiveStatusReachedTOF->GetXaxis()->SetBinLabel(2, "Has only TOF sig"); + hPositiveStatusReachedTOF->GetXaxis()->SetBinLabel(3, "Has only TOF ev time"); + hPositiveStatusReachedTOF->GetXaxis()->SetBinLabel(4, "Has full time info"); + hNegativeStatusReachedTOF->GetXaxis()->SetBinLabel(1, "All negative"); + hNegativeStatusReachedTOF->GetXaxis()->SetBinLabel(2, "Has only TOF sig"); + hNegativeStatusReachedTOF->GetXaxis()->SetBinLabel(3, "Has only TOF ev time"); + hNegativeStatusReachedTOF->GetXaxis()->SetBinLabel(4, "Has full time info"); + + // V0 Radius + histos.add("hLambdaMass", "hLambdaMass", {HistType::kTH1F, {axisLambdaMass}}); + histos.add("hAssocLambdaMass", "hAssocLambdaMass", {HistType::kTH1F, {axisLambdaMass}}); + histos.add("hAssocLambdaMassGoodTime", "hAssocLambdaMassGoodTime", {HistType::kTH1F, {axisLambdaMass}}); + histos.add("hAssocLambdaMassBadTime", "hAssocLambdaMassBadTime", {HistType::kTH1F, {axisLambdaMass}}); + // V0 Radius histos.add("h2dLambdaRadiusVsPt", "hLambdaRadiusVsPt", {HistType::kTH2F, {axisPt, axisRadius}}); // Invariant Mass histos.add("h2dLambdaMassVsPt", "hLambdaMassVsPt", {HistType::kTH2F, {axisPt, axisLambdaMass}}); + // Invariant Mass + histos.add("h2dLambdaMassVsTOFCut", "h2dLambdaMassVsTOFCut", {HistType::kTH2F, {axisLambdaMass, axisTOFCut}}); + histos.add("h2dLambdaMassVsTOFCutWithSignals", "h2dLambdaMassVsTOFCutWithSignals", {HistType::kTH2F, {axisLambdaMass, axisTOFCut}}); + histos.add("h2dLambdaMassVsTOFCutMeson", "h2dLambdaMassVsTOFCutMeson", {HistType::kTH2F, {axisLambdaMass, axisTOFCut}}); + histos.add("h2dLambdaMassVsTOFCutMesonWithSignals", "h2dLambdaMassVsTOFCutMesonWithSignals", {HistType::kTH2F, {axisLambdaMass, axisTOFCut}}); + + // Invariant Mass with TOF selections + histos.add("hLambdaMass_ProtonTOF", "hLambdaMass_ProtonTOF", {HistType::kTH1F, {axisLambdaMass}}); + histos.add("hLambdaMass_PionTOF", "hLambdaMass_PionTOF", {HistType::kTH1F, {axisLambdaMass}}); + histos.add("hLambdaMass_AllTOF", "hLambdaMass_AllTOF", {HistType::kTH1F, {axisLambdaMass}}); + histos.add("hLambdaMass_DeltaDecayTime", "hLambdaMass_DeltaDecayTime", {HistType::kTH1F, {axisLambdaMass}}); + histos.add("h2dLambdaMassVsPt_ProtonTOF", "hLambdaMassVsPtProtonTOF", {HistType::kTH2F, {axisPt, axisLambdaMass}}); + histos.add("h2dLambdaMassVsPt_PionTOF", "hLambdaMassVsPtPionTOF", {HistType::kTH2F, {axisPt, axisLambdaMass}}); + histos.add("h2dLambdaMassVsPt_AllTOF", "hLambdaMassVsPtAllTOF", {HistType::kTH2F, {axisPt, axisLambdaMass}}); + histos.add("h2dLambdaMassVsPt_DeltaDecayTime", "hLambdaMassVsPtDeltaDecayTime", {HistType::kTH2F, {axisPt, axisLambdaMass}}); + + // Invariant Mass with TOF selections + histos.add("hLambdaMass_InvertProtonTOF", "hLambdaMass_InvertProtonTOF", {HistType::kTH1F, {axisLambdaMass}}); + histos.add("hLambdaMass_InvertPionTOF", "hLambdaMass_InvertPionTOF", {HistType::kTH1F, {axisLambdaMass}}); + histos.add("hLambdaMass_InvertAllTOF", "hLambdaMass_InvertAllTOF", {HistType::kTH1F, {axisLambdaMass}}); + histos.add("h2dLambdaMassVsPt_InvertProtonTOF", "hLambdaMassVsPtInvertProtonTOF", {HistType::kTH2F, {axisPt, axisLambdaMass}}); + histos.add("h2dLambdaMassVsPt_InvertPionTOF", "hLambdaMassVsPtInvertPionTOF", {HistType::kTH2F, {axisPt, axisLambdaMass}}); + histos.add("h2dLambdaMassVsPt_InvertAllTOF", "hLambdaMassVsPtInvertAllTOF", {HistType::kTH2F, {axisPt, axisLambdaMass}}); + // radius vs prong length histos.add("h2dProtonLengthVsRadius", "h2dProtonLengthVsRadius", {HistType::kTH2F, {axisRadius, axisLength}}); histos.add("h2dPionLengthVsRadius", "h2dPionLengthVsRadius", {HistType::kTH2F, {axisRadius, axisLength}}); @@ -87,15 +164,132 @@ struct strangepidqa { histos.add("h2dProtonDeltaTimeVsRadius", "h2dProtonDeltaTimeVsRadius", {HistType::kTH2F, {axisRadius, axisDeltaTime}}); histos.add("h2dPionDeltaTimeVsPt", "h2dPionDeltaTimeVsPt", {HistType::kTH2F, {axisPt, axisDeltaTime}}); histos.add("h2dPionDeltaTimeVsRadius", "h2dPionDeltaTimeVsRadius", {HistType::kTH2F, {axisRadius, axisDeltaTime}}); + + // TOF PID testing for prongs + histos.add("h2dProtonDeltaTimeVsPt_MassSelected", "h2dProtonDeltaTimeVsPt_MassSelected", {HistType::kTH2F, {axisPt, axisDeltaTime}}); + histos.add("h2dProtonDeltaTimeVsRadius_MassSelected", "h2dProtonDeltaTimeVsRadius_MassSelected", {HistType::kTH2F, {axisRadius, axisDeltaTime}}); + histos.add("h2dPionDeltaTimeVsPt_MassSelected", "h2dPionDeltaTimeVsPt_MassSelected", {HistType::kTH2F, {axisPt, axisDeltaTime}}); + histos.add("h2dPionDeltaTimeVsRadius_MassSelected", "h2dPionDeltaTimeVsRadius_MassSelected", {HistType::kTH2F, {axisRadius, axisDeltaTime}}); + + // delta lambda decay time + histos.add("h2dLambdaDeltaDecayTime", "h2dLambdaDeltaDecayTime", {HistType::kTH2F, {axisPt, axisDeltaTime}}); + histos.add("h2dLambdaDeltaDecayTime_MassSelected", "h2dLambdaDeltaDecayTime_MassSelected", {HistType::kTH2F, {axisPt, axisDeltaTime}}); + + // beta plots + histos.add("h2dLambdaBeta", "h2dLambdaBeta", {HistType::kTH2F, {axisPt, axisBeta}}); + histos.add("h2dLambdaBeta_MassSelected", "h2dLambdaBeta_MassSelected", {HistType::kTH2F, {axisPt, axisBeta}}); + + // length vs time for prongs / debug + histos.add("h2dTimeVsLengthProtonProng", "h2dTimeVsLengthProtonProng", {HistType::kTH2F, {axisLength, axisTime}}); + histos.add("h2dTimeVsLengthPionProng", "h2dTimeVsLengthPionProng", {HistType::kTH2F, {axisLength, axisTime}}); + + // --- ASSOCIATED --- + // V0 Radius + if (doprocessSim) { + histos.add("h2dAssocLambdaRadiusVsPt", "hLambdaRadiusVsPt", {HistType::kTH2F, {axisPt, axisRadius}}); + + // Invariant Mass + histos.add("h2dAssocLambdaMassVsPt", "hLambdaMassVsPt", {HistType::kTH2F, {axisPt, axisLambdaMass}}); + + // radius vs prong length + histos.add("h2dAssocProtonLengthVsRadius", "h2dAssocProtonLengthVsRadius", {HistType::kTH2F, {axisRadius, axisLength}}); + histos.add("h2dAssocPionLengthVsRadius", "h2dAssocPionLengthVsRadius", {HistType::kTH2F, {axisRadius, axisLength}}); + + // recalculated vs topv lengths + histos.add("h2dAssocProtonLengthVsLengthToPV", "h2dAssocProtonLengthVsRadiusToPV", {HistType::kTH2F, {axisRadius, axisRadius}}); + histos.add("h2dAssocPionLengthVsLengthToPV", "h2dAssocPionLengthVsLengthToPV", {HistType::kTH2F, {axisRadius, axisRadius}}); + + // TOF PID testing for prongs + histos.add("h2dAssocProtonDeltaTimeVsPt", "h2dAssocProtonDeltaTimeVsPt", {HistType::kTH2F, {axisPt, axisDeltaTime}}); + histos.add("h2dAssocProtonDeltaTimeVsRadius", "h2dAssocProtonDeltaTimeVsRadius", {HistType::kTH2F, {axisRadius, axisDeltaTime}}); + histos.add("h2dAssocPionDeltaTimeVsPt", "h2dAssocPionDeltaTimeVsPt", {HistType::kTH2F, {axisPt, axisDeltaTime}}); + histos.add("h2dAssocPionDeltaTimeVsRadius", "h2dAssocPionDeltaTimeVsRadius", {HistType::kTH2F, {axisRadius, axisDeltaTime}}); + + // delta lambda decay time + histos.add("h2dAssocLambdaDeltaDecayTime", "h2dAssocLambdaDeltaDecayTime", {HistType::kTH2F, {axisPt, axisDeltaTime}}); + histos.add("h2dAssocLambdaDeltaDecayTime_MassSelected", "h2dAssocLambdaDeltaDecayTime_MassSelected", {HistType::kTH2F, {axisPt, axisDeltaTime}}); + } } - void process(aod::StraCollision const& coll, soa::Join const& v0s) + void processReal(soa::Join::iterator const& coll, soa::Join const& v0s) { histos.fill(HIST("hEventVertexZ"), coll.posZ()); + if (coll.centFT0C() > maxCentrality || coll.centFT0C() < minCentrality) + return; + for (auto& lambda : v0s) { // selecting photons from Sigma0 - if (lambda.pdgCode() != 3122) + + if (TMath::Abs(lambda.eta()) > 0.5) + continue; + + if (TMath::Abs(lambda.mLambda() - 1.115683) < 0.01 && lambda.v0cosPA() > minCosPA) { + histos.fill(HIST("hPositiveStatus"), 0.0f); + if (lambda.posTOFSignal() > 1e-5 && lambda.posTOFEventTime() < 0.0f) + histos.fill(HIST("hPositiveStatus"), 1.0f); + if (lambda.posTOFSignal() < 0.0f && lambda.posTOFEventTime() > 1e-5) + histos.fill(HIST("hPositiveStatus"), 2.0f); + if (lambda.posTOFSignal() > 1e-5 && lambda.posTOFEventTime() > 1e-5) + histos.fill(HIST("hPositiveStatus"), 3.0f); + + histos.fill(HIST("hNegativeStatus"), 0.0f); + if (lambda.negTOFSignal() > 1e-5 && lambda.negTOFEventTime() < 0.0f) + histos.fill(HIST("hNegativeStatus"), 1.0f); + if (lambda.negTOFSignal() < 0.0f && lambda.negTOFEventTime() > 1e-5) + histos.fill(HIST("hNegativeStatus"), 2.0f); + if (lambda.negTOFSignal() > 1e-5 && lambda.negTOFEventTime() > 1e-5) + histos.fill(HIST("hNegativeStatus"), 3.0f); + + if (lambda.posTOFLength() > 0.0f) { + histos.fill(HIST("hPositiveStatusReachedTOF"), 0.0f); + if (lambda.posTOFSignal() > 1e-5 && lambda.posTOFEventTime() < 0.0f) + histos.fill(HIST("hPositiveStatusReachedTOF"), 1.0f); + if (lambda.posTOFSignal() < 0.0f && lambda.posTOFEventTime() > 1e-5) + histos.fill(HIST("hPositiveStatusReachedTOF"), 2.0f); + if (lambda.posTOFSignal() > 1e-5 && lambda.posTOFEventTime() > 1e-5) + histos.fill(HIST("hPositiveStatusReachedTOF"), 3.0f); + } + + if (lambda.negTOFLength() > 0.0f) { + histos.fill(HIST("hNegativeStatusReachedTOF"), 0.0f); + if (lambda.negTOFSignal() > 1e-5 && lambda.negTOFEventTime() < 0.0f) + histos.fill(HIST("hNegativeStatusReachedTOF"), 1.0f); + if (lambda.negTOFSignal() < 0.0f && lambda.negTOFEventTime() > 1e-5) + histos.fill(HIST("hNegativeStatusReachedTOF"), 2.0f); + if (lambda.negTOFSignal() > 1e-5 && lambda.negTOFEventTime() > 1e-5) + histos.fill(HIST("hNegativeStatusReachedTOF"), 3.0f); + } + + // properties of the positive and negative prongs / debug + histos.fill(HIST("h2dTimeVsLengthProtonProng"), lambda.posTOFLength(), lambda.posLifetimePr()); + histos.fill(HIST("h2dTimeVsLengthPionProng"), lambda.negTOFLength(), lambda.negLifetimePr()); + } + + histos.fill(HIST("h2dLambdaMassVsTOFCut"), lambda.mLambda(), TMath::Abs(lambda.posTOFDeltaTLaPr())); + histos.fill(HIST("h2dLambdaMassVsTOFCutMeson"), lambda.mLambda(), TMath::Abs(lambda.negTOFDeltaTLaPi())); + + if (lambda.v0radius() < minV0Radius) + continue; + if (lambda.posTOFSignal() < 0 && requireTOFsignalProton) continue; + if (lambda.negTOFSignal() < 0 && requireTOFsignalPion) + continue; + + histos.fill(HIST("h2dLambdaDeltaDecayTime"), lambda.pt(), lambda.deltaDecayTimeLambda()); + if (TMath::Abs(lambda.mLambda() - 1.115683) < 0.01 && lambda.v0cosPA() > minCosPA) { + histos.fill(HIST("h2dLambdaDeltaDecayTime_MassSelected"), lambda.pt(), lambda.deltaDecayTimeLambda()); + } + + if (lambda.posTOFEventTime() < 0 && requireTOFEventTimeProton) + continue; + if (lambda.negTOFEventTime() < 0 && requireTOFEventTimePion) + continue; + + histos.fill(HIST("h2dLambdaMassVsTOFCutWithSignals"), lambda.mLambda(), TMath::Abs(lambda.posTOFDeltaTLaPr())); + histos.fill(HIST("h2dLambdaMassVsTOFCutMesonWithSignals"), lambda.mLambda(), TMath::Abs(lambda.negTOFDeltaTLaPi())); + + if (lambda.pt() > minPtV0 && lambda.pt() < maxPtV0) + histos.fill(HIST("hLambdaMass"), lambda.mLambda()); histos.fill(HIST("h2dLambdaRadiusVsPt"), lambda.pt(), lambda.v0radius()); histos.fill(HIST("h2dLambdaMassVsPt"), lambda.pt(), lambda.mLambda()); @@ -104,14 +298,109 @@ struct strangepidqa { histos.fill(HIST("h2dProtonDeltaTimeVsRadius"), lambda.v0radius(), lambda.posTOFDeltaTLaPr()); histos.fill(HIST("h2dPionDeltaTimeVsPt"), lambda.pt(), lambda.negTOFDeltaTLaPi()); histos.fill(HIST("h2dPionDeltaTimeVsRadius"), lambda.v0radius(), lambda.negTOFDeltaTLaPi()); + histos.fill(HIST("h2dLambdaBeta"), lambda.p(), lambda.tofBetaLambda()); + + if (TMath::Abs(lambda.mLambda() - 1.115683) < 0.01 && lambda.v0cosPA() > minCosPA) { + histos.fill(HIST("h2dProtonDeltaTimeVsPt_MassSelected"), lambda.pt(), lambda.posTOFDeltaTLaPr()); + histos.fill(HIST("h2dProtonDeltaTimeVsRadius_MassSelected"), lambda.v0radius(), lambda.posTOFDeltaTLaPr()); + histos.fill(HIST("h2dPionDeltaTimeVsPt_MassSelected"), lambda.pt(), lambda.negTOFDeltaTLaPi()); + histos.fill(HIST("h2dPionDeltaTimeVsRadius_MassSelected"), lambda.v0radius(), lambda.negTOFDeltaTLaPi()); + histos.fill(HIST("h2dLambdaBeta_MassSelected"), lambda.p(), lambda.tofBetaLambda()); + } histos.fill(HIST("h2dProtonLengthVsRadius"), lambda.v0radius(), lambda.posTOFLength()); histos.fill(HIST("h2dPionLengthVsRadius"), lambda.v0radius(), lambda.negTOFLength()); histos.fill(HIST("h2dProtonLengthVsLengthToPV"), lambda.posTOFLengthToPV(), lambda.posTOFLength()); histos.fill(HIST("h2dPionLengthVsLengthToPV"), lambda.negTOFLengthToPV(), lambda.negTOFLength()); + + // Standard selection of time + if (TMath::Abs(lambda.deltaDecayTimeLambda()) < maxDeltaTimeDecay) { + histos.fill(HIST("h2dLambdaMassVsPt_DeltaDecayTime"), lambda.pt(), lambda.mLambda()); + if (lambda.pt() > minPtV0 && lambda.pt() < maxPtV0) + histos.fill(HIST("hLambdaMass_DeltaDecayTime"), lambda.mLambda()); + } + + if (TMath::Abs(lambda.posTOFDeltaTLaPr()) < maxDeltaTimeProton) { + histos.fill(HIST("h2dLambdaMassVsPt_ProtonTOF"), lambda.pt(), lambda.mLambda()); + if (lambda.pt() > minPtV0 && lambda.pt() < maxPtV0) + histos.fill(HIST("hLambdaMass_ProtonTOF"), lambda.mLambda()); + } + if (TMath::Abs(lambda.negTOFDeltaTLaPi()) < maxDeltaTimeProton) { + histos.fill(HIST("h2dLambdaMassVsPt_PionTOF"), lambda.pt(), lambda.mLambda()); + if (lambda.pt() > minPtV0 && lambda.pt() < maxPtV0) + histos.fill(HIST("hLambdaMass_PionTOF"), lambda.mLambda()); + if (TMath::Abs(lambda.posTOFDeltaTLaPr()) < maxDeltaTimeProton) { + histos.fill(HIST("h2dLambdaMassVsPt_AllTOF"), lambda.pt(), lambda.mLambda()); + if (lambda.pt() > minPtV0 && lambda.pt() < maxPtV0) + histos.fill(HIST("hLambdaMass_AllTOF"), lambda.mLambda()); + } + } + // Inversion of time selection for debug + // if(TMath::Abs(lambda.deltaDecayTimeLambda())>maxDeltaTimeDecay && TMath::Abs(lambda.deltaDecayTimeLambda()) < 5e+4){ + // histos.fill(HIST("h2dLambdaMassVsPt_DeltaDecayTime"), lambda.pt(), lambda.mLambda()); + // if(lambda.pt()>minPtV0 && lambda.pt() < maxPtV0) histos.fill(HIST("hLambdaMass_DeltaDecayTime"), lambda.mLambda()); + // } + + if (TMath::Abs(lambda.posTOFDeltaTLaPr()) > maxDeltaTimeProton && TMath::Abs(lambda.posTOFDeltaTLaPr()) < 5e+4) { + histos.fill(HIST("h2dLambdaMassVsPt_InvertProtonTOF"), lambda.pt(), lambda.mLambda()); + if (lambda.pt() > minPtV0 && lambda.pt() < maxPtV0) + histos.fill(HIST("hLambdaMass_InvertProtonTOF"), lambda.mLambda()); + } + if (TMath::Abs(lambda.negTOFDeltaTLaPi()) > maxDeltaTimeProton && TMath::Abs(lambda.negTOFDeltaTLaPi()) < 5e+4) { + histos.fill(HIST("h2dLambdaMassVsPt_InvertPionTOF"), lambda.pt(), lambda.mLambda()); + if (lambda.pt() > minPtV0 && lambda.pt() < maxPtV0) + histos.fill(HIST("hLambdaMass_InvertPionTOF"), lambda.mLambda()); + if (TMath::Abs(lambda.posTOFDeltaTLaPr()) > maxDeltaTimeProton && TMath::Abs(lambda.posTOFDeltaTLaPr()) < 5e+4) { + histos.fill(HIST("h2dLambdaMassVsPt_InvertAllTOF"), lambda.pt(), lambda.mLambda()); + if (lambda.pt() > minPtV0 && lambda.pt() < maxPtV0) + histos.fill(HIST("hLambdaMass_InvertAllTOF"), lambda.mLambda()); + } + } + } + } + + void processSim(aod::StraCollision const& coll, soa::Join const& v0s) + { + for (auto& lambda : v0s) { // selecting photons from Sigma0 + if (lambda.pdgCode() != 3122) + continue; + if (!lambda.isPhysicalPrimary()) + continue; + if (lambda.pdgCodePositive() != 2212) + continue; + if (lambda.pdgCodeNegative() != -211) + continue; + + histos.fill(HIST("hAssocLambdaMass"), lambda.mLambda()); + + if (TMath::Abs(lambda.negTOFDeltaTLaPi()) > 800) { + histos.fill(HIST("hAssocLambdaMassGoodTime"), lambda.mLambda()); + } else { + histos.fill(HIST("hAssocLambdaMassBadTime"), lambda.mLambda()); + } + + histos.fill(HIST("h2dAssocLambdaRadiusVsPt"), lambda.pt(), lambda.v0radius()); + histos.fill(HIST("h2dAssocLambdaMassVsPt"), lambda.pt(), lambda.mLambda()); + + histos.fill(HIST("h2dAssocProtonDeltaTimeVsPt"), lambda.pt(), lambda.posTOFDeltaTLaPr()); + histos.fill(HIST("h2dAssocProtonDeltaTimeVsRadius"), lambda.v0radius(), lambda.posTOFDeltaTLaPr()); + histos.fill(HIST("h2dAssocPionDeltaTimeVsPt"), lambda.pt(), lambda.negTOFDeltaTLaPi()); + histos.fill(HIST("h2dAssocPionDeltaTimeVsRadius"), lambda.v0radius(), lambda.negTOFDeltaTLaPi()); + + histos.fill(HIST("h2dAssocProtonLengthVsRadius"), lambda.v0radius(), lambda.posTOFLength()); + histos.fill(HIST("h2dAssocPionLengthVsRadius"), lambda.v0radius(), lambda.negTOFLength()); + + histos.fill(HIST("h2dAssocProtonLengthVsLengthToPV"), lambda.posTOFLengthToPV(), lambda.posTOFLength()); + histos.fill(HIST("h2dAssocPionLengthVsLengthToPV"), lambda.negTOFLengthToPV(), lambda.negTOFLength()); + + // delta lambda decay time + histos.fill(HIST("h2dAssocLambdaDeltaDecayTime"), lambda.pt(), lambda.deltaDecayTimeLambda()); } } + + PROCESS_SWITCH(strangepidqa, processReal, "Produce real information", true); + PROCESS_SWITCH(strangepidqa, processSim, "Produce simulated information", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)