diff --git a/PWGLF/Tasks/Resonances/lambdastarproxy.cxx b/PWGLF/Tasks/Resonances/lambdastarproxy.cxx index b923950622e..dc677ad52da 100644 --- a/PWGLF/Tasks/Resonances/lambdastarproxy.cxx +++ b/PWGLF/Tasks/Resonances/lambdastarproxy.cxx @@ -45,6 +45,7 @@ #include #include #include +#include #include #include #include @@ -148,6 +149,22 @@ struct Lambdastarproxy { Configurable lstarEnableExpSignalTOFDe{"lstarEnableExpSignalTOFDe", 0, "Enable deuteron-only TOF expected-signal-difference cut when lstarTOFExpSignalDiffDeMax > 0"}; Configurable lstarTOFExpSignalDiffDeMax{"lstarTOFExpSignalDiffDeMax", -1.f, "Maximum |tofExpSignalDiffDe| for deuterons; <=0 disables the numeric cut"}; + // Opt-in deuteron-only refinements. Disabled defaults preserve the old + // selections, histogram schema and K/p PID, including legacy strategy 2. + struct : ConfigurableGroup { + Configurable lstarEnableITSPIDDe{"lstarEnableITSPIDDe", 0, "Apply one-sided ITS cluster-size PID to deuteron candidates"}; + Configurable lstarCutNsigmaITSDeMin{"lstarCutNsigmaITSDeMin", -3.f, "Require nSigmaITS_d greater than this lower bound (not an absolute-value cut)"}; + Configurable lstarITSNClustersDe{"lstarITSNClustersDe", 0, "Additional deuteron-only minimum ITS clusters; 0 disables"}; + Configurable lstarTPCNClsFoundDeMin{"lstarTPCNClsFoundDeMin", 0, "Additional deuteron-only minimum TPC found clusters (not crossed rows); 0 disables"}; + Configurable lstarRequireTOFAbovePDe{"lstarRequireTOFAbovePDe", -1.f, "Strategy 2: require TOF match and TOF-d nSigma cut at full track p >= threshold in GeV/c; negative disables"}; + Configurable lstarEnableDeuteronPurityQA{"lstarEnableDeuteronPurityQA", 0, "Book compact charge-separated pre-PID/selected d response and unique-track MC species controls; independent of legacy PID QA"}; + // Z=1 defaults and formula follow Common/DataModel/PIDResponseITS.h. + // Keep parameters per task instance: shared ITSResponse initialization would + // introduce an its-pid dependency and mutable state across cloned wagons. + Configurable> lstarITSDeParamsData{"lstarITSDeParamsData", {1.18941f, 1.53792f, 1.69961f, 0.194669f, -0.208616f, 1.30753f}, "Data ITS d response: A,B,C,r0,r1,r2; expected=A/(p/m)^B+C; relative sigma=r0*erf((p/m-r1)/r2)"}; + Configurable> lstarITSDeParamsMC{"lstarITSDeParamsMC", {1.63806f, 1.58847f, 2.52275f, 0.140487f, -0.431078f, 1.50052f}, "Reconstructed-MC ITS d response: A,B,C,r0,r1,r2; override for production-specific calibration"}; + } deuteronPID; + // PID strategy for final K/p/d candidate selection. // 0 = pT-ref dependent rectangular cuts: // pT < pTref: require |TPC| < TPC cut only @@ -179,6 +196,16 @@ struct Lambdastarproxy { Configurable lstarTPCCrossedRowsOverFindableMin{"lstarTPCCrossedRowsOverFindableMin", 0.8f, "Minimum TPC crossed rows over findable clusters"}; Configurable lstarDcaXYMax{"lstarDcaXYMax", float{DcaXYMaxDefault}, "Max |DCA_{xy}| to PV (cm)"}; Configurable lstarDcaZMax{"lstarDcaZMax", float{DcaZMaxDefault}, "Max |DCA_{z}| to PV (cm)"}; + // Optional common-track DCAz model. The fixed cut above is unchanged when + // disabled. pT is the original measured track pT, before proxy scaling. + struct : ConfigurableGroup { + Configurable lstarEnablePtDepDcaZ{"lstarEnablePtDepDcaZ", 0, "Replace the explicit fixed DCAz cut by a pT-dependent model for K/p/d tracks; requires lstarOnlyGlobalTrackCuts=0"}; + Configurable> lstarPtDepDcaZParams{"lstarPtDepDcaZParams", {0.004f, 0.013f, 1.f, 1.f}, "DCAz parameters {offset,coefficient,power,scale}: max|DCAz|=scale*(offset+coefficient/pT^power) cm, pT in GeV/c; not a calibrated DCAz resolution"}; + } dcaZSelection; + struct : ConfigurableGroup { + Configurable lstarEnableTrackSpectra{"lstarEnableTrackSpectra", 0, "Book charge-separated raw selected d/p/K track pT and pT-versus-FT0M spectra, independent of PID QA and pair partners"}; + Configurable lstarTrackSpectraAbsYMax{"lstarTrackSpectraAbsYMax", 0.5f, "Spectra-only |y| window under each species mass; does not change pair candidates; raw spectra need corrections before published-yield comparison"}; + } trackSpectraConfig; Configurable lstarChi2PerTPCCluster{"lstarChi2PerTPCCluster", float{Chi2PerTPCClusterDefault}, "Maximum Chi2/cluster for the TPC track segment"}; Configurable lstarChi2PerITSCluster{"lstarChi2PerITSCluster", float{Chi2PerITSClusterDefault}, "Maximum Chi2/cluster for the ITS track segment"}; Configurable lstarITSNClusters{"lstarITSNClusters", int{ITSNClustersDefault}, "Minimum number of ITS clusters"}; @@ -441,6 +468,42 @@ struct Lambdastarproxy { aod::pidTPCFullDe, aod::pidTOFFullDe, aod::pidTOFbeta>; + static constexpr size_t DcaZPtParameterCount = 4; + static constexpr size_t DcaZPtPowerIndex = 2; + static constexpr size_t DcaZPtScaleIndex = 3; + + void validatePtDependentDcaZ() const + { + if (dcaZSelection.lstarEnablePtDepDcaZ.value == 0) { + return; + } + if (lstarOnlyGlobalTrackCuts.value != 0) { + LOGF(fatal, "lstarEnablePtDepDcaZ requires lstarOnlyGlobalTrackCuts=0; global/primary-track flags still apply independently"); + } + auto const& parameters = dcaZSelection.lstarPtDepDcaZParams.value; + if (parameters.size() != DcaZPtParameterCount || + !std::all_of(parameters.begin(), parameters.end(), [](float value) { return std::isfinite(value); }) || + parameters[0] < 0.f || parameters[1] < 0.f || (parameters[0] == 0.f && parameters[1] == 0.f) || + parameters[DcaZPtPowerIndex] < 0.f || parameters[DcaZPtScaleIndex] <= 0.f) { + LOGF(fatal, "DCAz pT model requires four finite {offset,coefficient,power,scale} values; offset/coefficient nonnegative with a nonzero term, power nonnegative, scale positive"); + } + } + + bool passDcaZSelection(float pt, float dcaZ) const + { + if (dcaZSelection.lstarEnablePtDepDcaZ.value == 0) { + // Preserve the legacy fixed-cut comparison, including its disabled-cut + // convention and handling of nonfinite values, for exact old selection. + return !(lstarDcaZMax.value > 0.f && std::abs(dcaZ) > lstarDcaZMax.value); + } + if (!std::isfinite(pt) || pt <= 0.f || !std::isfinite(dcaZ)) { + return false; + } + auto const& parameters = dcaZSelection.lstarPtDepDcaZParams.value; + const float maximum = parameters[DcaZPtScaleIndex] * (parameters[0] + parameters[1] / std::pow(pt, parameters[DcaZPtPowerIndex])); + return std::isfinite(maximum) && maximum > 0.f && std::abs(dcaZ) <= maximum; + } + template bool passTrackQuality(const TTrack& trk) const { @@ -513,7 +576,7 @@ struct Lambdastarproxy { } if constexpr (requires { trk.dcaZ(); }) { - if (lstarDcaZMax.value > 0.f && std::abs(trk.dcaZ()) > lstarDcaZMax.value) { + if (!passDcaZSelection(trk.pt(), trk.dcaZ())) { return false; } } @@ -554,6 +617,12 @@ struct Lambdastarproxy { if (lstarNoMixedEvents.value < 0 || lstarSparseMultBins.value <= 0) { LOGF(fatal, "lstarNoMixedEvents must be nonnegative; lstarSparseMultBins must be positive"); } + validateDeuteronRefinements(); + validatePtDependentDcaZ(); + if (trackSpectraConfig.lstarEnableTrackSpectra.value != 0 && + ((!doprocess.value && !doprocessMCReco.value) || !std::isfinite(trackSpectraConfig.lstarTrackSpectraAbsYMax.value) || trackSpectraConfig.lstarTrackSpectraAbsYMax.value <= 0.f)) { + LOGF(fatal, "Track spectra require data or reconstructed MC and a positive finite spectra rapidity window"); + } AxisSpec massAxis{200, 1.4, 1.9, "M_{pK} (GeV/c^{2})"}; AxisSpec ptAxis{100, 0., 10., "p_{T} (GeV/c)"}; AxisSpec nsAxis{500, -50., 50., "n#sigma"}; @@ -1028,10 +1097,117 @@ struct Lambdastarproxy { } bookMC(); bookProxyDiagnostics(); + bookDeuteronPurityQA(); + bookTrackSpectra(); } static constexpr int PdgLambdaStar = static_cast(PDG_t::kLambda1520); static constexpr int PdgDeuteron = static_cast(o2::constants::physics::Pdg::kDeuteron); + enum class TrackSpectrumSpecies : size_t { Deuteron, + Proton, + Kaon }; + static constexpr size_t TrackSpectrumSpeciesCount = 3; + static constexpr size_t TrackSpectrumChargeCount = 2; + static constexpr int TrackSpectrumPtBins = 200; + static constexpr float TrackSpectrumPtMax = 10.f; + static constexpr std::array, TrackSpectrumSpeciesCount> TrackSpectrumNames{{{"Deuteron", "AntiDeuteron"}, {"Proton", "AntiProton"}, {"KaonPlus", "KaonMinus"}}}; + static constexpr std::array TrackSpectrumMasses{o2::constants::physics::MassDeuteron, o2::constants::physics::MassProton, o2::constants::physics::MassKaonCharged}; + static constexpr std::array TrackSpectrumPdgs{PdgDeuteron, PdgProton, PdgKaon}; + enum class TrackSpectrumTruth : int { Unmatched, + WrongIdentity, + WrongCollision, + NonPrimary, + Primary }; + static constexpr int TrackSpectrumTruthClasses = 5; + struct TrackSpectraHistograms { + std::array, TrackSpectrumChargeCount>, TrackSpectrumSpeciesCount> pt{}; + std::array, TrackSpectrumChargeCount>, TrackSpectrumSpeciesCount> ptFT0M{}; + std::shared_ptr events, eventsFT0M; + std::shared_ptr mcTruth; + } trackSpectra; + + void bookTrackSpectra() + { + if (trackSpectraConfig.lstarEnableTrackSpectra.value == 0) { + return; + } + const AxisSpec pt{TrackSpectrumPtBins, 0., TrackSpectrumPtMax, "Full measured track pT (GeV/c)"}; + const AxisSpec ft0m{FT0MPercentileBins, 0., FT0MPercentileMax, "Raw FT0M percentile"}; + for (size_t species = 0; species < TrackSpectrumSpeciesCount; ++species) { + for (size_t charge = 0; charge < TrackSpectrumChargeCount; ++charge) { + const std::string name = TrackSpectrumNames[species][charge]; + trackSpectra.pt[species][charge] = histos.add(("hTrackPt" + name).c_str(), "Raw selected single-track candidates;full measured track pT (GeV/c);Counts;species-mass rapidity window", HistType::kTH1D, {pt}, true); + trackSpectra.ptFT0M[species][charge] = histos.add(("hTrackPtVsFT0M" + name).c_str(), "Raw selected single-track candidates;regular raw FT0M only;species-mass rapidity window", HistType::kTH2D, {pt, ft0m}, true); + } + } + trackSpectra.events = histos.add("hTrackSpectraEventSelection", "Selected reconstructed events, including zero candidates;1=all,2=regular FT0M,3=invalid or unavailable FT0M;Events", HistType::kTH1D, {AxisSpec{3, 0.5, 3.5, "Event scope"}}, true); + trackSpectra.eventsFT0M = histos.add("hTrackSpectraEventsVsFT0M", "Selected events for pT-versus-FT0M normalization, including zero candidates;regular raw FT0M percentile;Events", HistType::kTH1D, {ft0m}, true); + if (doprocessMCReco) { + trackSpectra.mcTruth = histos.add("hMCTrackSpectraTruth", "Same raw selected tracks:0=unmatched,1=wrong species/charge,2=wrong MC collision,3=correct nonprimary,4=correct primary;full pT;role=2*species+negativeCharge (d=0,p=1,K=2);truth class", + HistType::kTH3D, {pt, AxisSpec{TrackSpectrumSpeciesCount * TrackSpectrumChargeCount, -0.5, TrackSpectrumSpeciesCount * TrackSpectrumChargeCount - 0.5, "Species/charge role"}, AxisSpec{TrackSpectrumTruthClasses, -0.5, TrackSpectrumTruthClasses - 0.5, "Truth class"}}, true); + } + } + + template + float fillTrackSpectraEvent(TCollision const& collision) + { + float rawFT0M = std::numeric_limits::quiet_NaN(); + if (trackSpectraConfig.lstarEnableTrackSpectra.value == 0) { + return rawFT0M; + } + if constexpr (requires { collision.centFT0M(); }) { + rawFT0M = collision.centFT0M(); + } + trackSpectra.events->Fill(1.); + if (regularMultiplicity(rawFT0M)) { + trackSpectra.events->Fill(2.); + trackSpectra.eventsFT0M->Fill(rawFT0M); + } else { + trackSpectra.events->Fill(3.); + } + return rawFT0M; + } + + template + void fillTrackSpectrum(TTrack const& track, TrackSpectrumSpecies role, float rawFT0M, int mcCollisionId) + { + if (trackSpectraConfig.lstarEnableTrackSpectra.value == 0) { + return; + } + const size_t species = static_cast(role); + const size_t charge = track.sign() < 0 ? 1 : 0; + const double pt = track.pt(); + const double eta = track.eta(); + if (track.sign() == 0 || !std::isfinite(pt) || pt <= 0. || !std::isfinite(eta)) { + return; + } + const double mass = TrackSpectrumMasses[species]; + const double y = std::asinh(pt * std::sinh(eta) / std::sqrt(mass * mass + pt * pt)); + if (!std::isfinite(y) || std::abs(y) >= trackSpectraConfig.lstarTrackSpectraAbsYMax.value) { + return; + } + // Called once in each selected-species loop, never inside pair/mixing loops. + // The spectra-only rapidity cut does not remove candidates from those loops. + trackSpectra.pt[species][charge]->Fill(pt); + if (regularMultiplicity(rawFT0M)) { + trackSpectra.ptFT0M[species][charge]->Fill(pt, rawFT0M); + } + if constexpr (IsMC) { + auto truth = TrackSpectrumTruth::Unmatched; + if (track.has_mcParticle()) { + const auto& particle = track.mcParticle(); + const int expected = track.sign() > 0 ? TrackSpectrumPdgs[species] : -TrackSpectrumPdgs[species]; + if (particle.pdgCode() != expected) { + truth = TrackSpectrumTruth::WrongIdentity; + } else if (mcCollisionId < 0 || particle.mcCollisionId() != mcCollisionId) { + truth = TrackSpectrumTruth::WrongCollision; + } else { + truth = particle.isPhysicalPrimary() ? TrackSpectrumTruth::Primary : TrackSpectrumTruth::NonPrimary; + } + } + trackSpectra.mcTruth->Fill(pt, TrackSpectrumChargeCount * species + charge, static_cast(truth)); + } + } // MC outputs are 2D mass/pT histograms. They add no high-dimensional maps. std::array, 5> mcGenSpectra{}; std::array, 2>, 4> mcRecoProxy{}; @@ -1075,6 +1251,200 @@ struct Lambdastarproxy { std::shared_ptr probeTrackStatus; std::shared_ptr probeEventScope; + static constexpr int ITSClusterLayers = 7; + static constexpr int ITSClusterBits = 4; + static constexpr uint32_t ITSClusterMask = 0xf; + static constexpr int MinITSClustersForPID = 2; + static constexpr size_t ITSResponseParameters = 6; + static constexpr size_t ITSResolutionScale = 3; + static constexpr size_t ITSResolutionShift = 4; + static constexpr size_t ITSResolutionWidth = 5; + enum class DeuteronPIDStage : int { BeforePID, + AfterTPC, + Selected }; + static constexpr int DeuteronPIDStages = 3; + static constexpr int DeuteronPIDStatesPerCharge = 2 * DeuteronPIDStages; + static constexpr int DeuteronTruthUnmatched = SpeciesCount; + static constexpr int DeuteronTruthWrongCharge = SpeciesCount + 1; + static constexpr std::array DeuteronCharges{"D", "AntiD"}; + struct DeuteronPurityHistograms { + std::array, 2> dePurityTPCBefore{}, dePurityITSBefore{}, dePurityTPCSelected{}, dePurityITSSelected{}; + std::array, 2> dePurityMass2AfterTPC{}, dePurityMass2Selected{}, dePurityTrackStages{}; + std::shared_ptr dePurityMCSpecies; + } dePurity; + + void validateDeuteronRefinements() const + { + if (deuteronPID.lstarITSNClustersDe.value < 0 || deuteronPID.lstarITSNClustersDe.value > ITSClusterLayers || deuteronPID.lstarTPCNClsFoundDeMin.value < 0) { + LOGF(fatal, "Deuteron-only ITS minimum must be 0..7; TPC-found minimum must be nonnegative"); + } + if (!std::isfinite(deuteronPID.lstarRequireTOFAbovePDe.value)) { + LOGF(fatal, "lstarRequireTOFAbovePDe must be finite; a negative value disables it"); + } + if (deuteronPID.lstarRequireTOFAbovePDe.value >= 0.f && + (deuteronPID.lstarRequireTOFAbovePDe.value == 0.f || lstarPidStrategy.value != PidStrategyNucleiDeuteronTPC || + !std::isfinite(lstarCutNsigmaTOFDe.value) || lstarCutNsigmaTOFDe.value <= 0.f)) { + LOGF(fatal, "Momentum-dependent deuteron TOF requires strategy 2, a positive p threshold and a positive finite TOF-d cut"); + } + if (deuteronPID.lstarEnableITSPIDDe.value != 0 || deuteronPID.lstarEnableDeuteronPurityQA.value != 0) { + if (!std::isfinite(deuteronPID.lstarCutNsigmaITSDeMin.value)) { + LOGF(fatal, "The ITS-d lower nSigma bound must be finite"); + } + auto const& parameters = doprocessMCReco.value ? deuteronPID.lstarITSDeParamsMC.value : deuteronPID.lstarITSDeParamsData.value; + if (parameters.size() != ITSResponseParameters || + !std::all_of(parameters.begin(), parameters.end(), [](float value) { return std::isfinite(value); }) || + parameters[0] <= 0.f || parameters[1] <= 0.f || parameters[2] <= 0.f || + parameters[ITSResolutionScale] <= 0.f || parameters[ITSResolutionWidth] <= 0.f) { + LOGF(fatal, "ITS-d response requires six finite A,B,C,r0,r1,r2 parameters with A,B,C,r0,r2 positive"); + } + } + if (deuteronPID.lstarEnableDeuteronPurityQA.value != 0 && !doprocess.value && !doprocessMCReco.value) { + LOGF(fatal, "Deuteron purity QA requires data or reconstructed MC"); + } + } + + template + double deuteronITSNSigma(TTrack const& track) const + { + const double invalid = std::numeric_limits::quiet_NaN(); + const double momentum = track.p(); + if (!std::isfinite(momentum) || momentum <= 0. || !std::isfinite(track.eta())) { + return invalid; + } + int clusters = 0; + int sum = 0; + int largest = 0; + const uint32_t packed = track.itsClusterSizes(); + for (int layer = 0; layer < ITSClusterLayers; ++layer) { + const int size = (packed >> (layer * ITSClusterBits)) & ITSClusterMask; + if (size > 0) { + ++clusters; + sum += size; + largest = std::max(largest, size); + } + } + if (clusters < MinITSClustersForPID) { + return invalid; + } + auto const& parameters = doprocessMCReco.value ? deuteronPID.lstarITSDeParamsMC.value : deuteronPID.lstarITSDeParamsData.value; + const double betaGamma = momentum / o2::constants::physics::MassDeuteron; + const double expected = parameters[0] / std::pow(betaGamma, parameters[1]) + parameters[2]; + const double sigma = expected * parameters[ITSResolutionScale] * std::erf((betaGamma - parameters[ITSResolutionShift]) / parameters[ITSResolutionWidth]); + if (!std::isfinite(expected) || !std::isfinite(sigma) || expected <= 0. || sigma <= 0.) { + return invalid; + } + const double truncatedMean = static_cast(sum - largest) / (clusters - 1); + const double response = (truncatedMean / std::cosh(static_cast(track.eta())) - expected) / sigma; + return std::isfinite(response) ? response : invalid; + } + + template + bool passDeuteronRefinements(TTrack const& track) const + { + if (deuteronPID.lstarITSNClustersDe.value > 0 && track.itsNCls() < deuteronPID.lstarITSNClustersDe.value) { + return false; + } + if (deuteronPID.lstarTPCNClsFoundDeMin.value > 0 && track.tpcNClsFound() < deuteronPID.lstarTPCNClsFoundDeMin.value) { + return false; + } + if (deuteronPID.lstarEnableITSPIDDe.value != 0) { + const double response = deuteronITSNSigma(track); + if (!std::isfinite(response) || response <= deuteronPID.lstarCutNsigmaITSDeMin.value) { + return false; + } + } + if (deuteronPID.lstarRequireTOFAbovePDe.value >= 0.f) { + // Full measured total p, before proxy division by two; equality belongs + // to the TOF-required branch. The existing global requirement is stronger. + const double momentum = track.p(); + if (!std::isfinite(momentum) || momentum <= 0.) { + return false; + } + if (momentum >= deuteronPID.lstarRequireTOFAbovePDe.value && + (!track.hasTOF() || !std::isfinite(track.tofNSigmaDe()) || std::abs(track.tofNSigmaDe()) >= lstarCutNsigmaTOFDe.value)) { + return false; + } + } + return true; + } + + void bookDeuteronPurityQA() + { + if (deuteronPID.lstarEnableDeuteronPurityQA.value == 0) { + return; + } + const AxisSpec momentum{100, 0., 10., "Full measured d-track p (GeV/c)"}; + const AxisSpec tpc{160, -12., 12., "TPC nSigma_d"}; + const AxisSpec its{100, -10., 10., "ITS nSigma_d"}; + const AxisSpec mass2{200, -4., 4., "TOF mass2 minus deuteron mass2 (GeV2/c4)"}; + uint64_t payload = 0; + auto add = [&](std::string const& name, std::string const& title, AxisSpec const& response) { + auto histogram = histos.add(name.c_str(), title.c_str(), HistType::kTH2D, {momentum, response}, true); + payload += static_cast(histogram->GetNcells()) * 2 * sizeof(double); + return histogram; + }; + for (size_t sign = 0; sign < DeuteronCharges.size(); ++sign) { + const std::string suffix = DeuteronCharges[sign]; + dePurity.dePurityTPCBefore[sign] = add("hDeuteronPurityTPCBefore" + suffix, "After common track quality, before all d-only cuts;all TOF states", tpc); + dePurity.dePurityITSBefore[sign] = add("hDeuteronPurityITSBefore" + suffix, "Same pre-PID tracks;valid ITS cluster response only", its); + dePurity.dePurityTPCSelected[sign] = add("hDeuteronPurityTPCSelected" + suffix, "Final selected proxy tracks;all TOF states", tpc); + dePurity.dePurityITSSelected[sign] = add("hDeuteronPurityITSSelected" + suffix, "Final selected proxy tracks;valid ITS response only", its); + dePurity.dePurityMass2AfterTPC[sign] = add("hDeuteronPurityMass2AfterTPC" + suffix, "After only TPC-d window and common track quality;valid TOF;before d-only and TOF PID cuts", mass2); + dePurity.dePurityMass2Selected[sign] = add("hDeuteronPurityMass2Selected" + suffix, "Final selected proxy tracks;valid TOF beta only", mass2); + dePurity.dePurityTrackStages[sign] = add("hDeuteronPurityTrackStages" + suffix, "Unique track rows per stage;full p;hasTOF+2*stage (0=common quality,1=TPC window,2=final proxy selection)", + AxisSpec{DeuteronPIDStatesPerCharge, -0.5, DeuteronPIDStatesPerCharge - 0.5, "hasTOF+2*stage"}); + } + if (doprocessMCReco) { + dePurity.dePurityMCSpecies = histos.add("hMCDeuteronPurityTrackSpecies", "Unique track rows: e,mu,pi,K,p,d,other,unmatched,wrong-charge;full p;species;hasTOF+2*stage+6*negativeCharge", + HistType::kTH3D, {momentum, AxisSpec{SpeciesCount + 2, -0.5, SpeciesCount + 1.5, "True single-track species"}, AxisSpec{2 * DeuteronPIDStatesPerCharge, -0.5, 2 * DeuteronPIDStatesPerCharge - 0.5, "Detector, stage and charge state"}}, true); + payload += static_cast(dePurity.dePurityMCSpecies->GetNcells()) * 2 * sizeof(double); + } + LOGF(info, "Deuteron purity QA content+Sumw2 payload: %.2f MiB per wagon (fixed bins; no THnSparse)", static_cast(payload) / (1024. * 1024.)); + } + + template + void fillDeuteronPurityQA(TTrack const& track, DeuteronPIDStage stage) + { + if (deuteronPID.lstarEnableDeuteronPurityQA.value == 0) { + return; + } + const double momentum = track.p(); + if (!std::isfinite(momentum) || momentum <= 0.) { + return; + } + const size_t sign = track.sign() < 0 ? 1 : 0; + const int state = static_cast(track.hasTOF()) + 2 * static_cast(stage); + dePurity.dePurityTrackStages[sign]->Fill(momentum, state); + if constexpr (IsMC) { + int species = DeuteronTruthUnmatched; + if (track.has_mcParticle()) { + const auto& particle = track.mcParticle(); + species = chargeFromPdg(particle.pdgCode()) == track.sign() ? speciesGroup(particle.pdgCode()) : DeuteronTruthWrongCharge; + } + dePurity.dePurityMCSpecies->Fill(momentum, species, state + DeuteronPIDStatesPerCharge * sign); + } + if (stage == DeuteronPIDStage::BeforePID || stage == DeuteronPIDStage::Selected) { + const double tpc = track.tpcNSigmaDe(); + const double its = deuteronITSNSigma(track); + if (std::isfinite(tpc)) { + (stage == DeuteronPIDStage::BeforePID ? dePurity.dePurityTPCBefore[sign] : dePurity.dePurityTPCSelected[sign])->Fill(momentum, tpc); + } + if (std::isfinite(its)) { + (stage == DeuteronPIDStage::BeforePID ? dePurity.dePurityITSBefore[sign] : dePurity.dePurityITSSelected[sign])->Fill(momentum, its); + } + } + if (stage != DeuteronPIDStage::BeforePID && track.hasTOF()) { + const double beta = track.beta(); + if (std::isfinite(beta) && beta > TofBetaMin && beta < TofBetaMax) { + constexpr double MassDeuteronForQA = o2::constants::physics::MassDeuteron; + const double delta = momentum * momentum * (1. / (beta * beta) - 1.) - MassDeuteronForQA * MassDeuteronForQA; + if (std::isfinite(delta)) { + (stage == DeuteronPIDStage::AfterTPC ? dePurity.dePurityMass2AfterTPC[sign] : dePurity.dePurityMass2Selected[sign])->Fill(momentum, delta); + } + } + } + } + struct PIDProbeCand { float px = 0.f, py = 0.f, pz = 0.f; double momentum = 0., mass2 = 0.; @@ -1332,7 +1702,7 @@ struct Lambdastarproxy { const bool passProton = passFinalCandidatePID(track.pt(), track.tpcNSigmaPr(), track.tofNSigmaPr(), true, lstarCutNsigmaTPCPr.value, lstarCutNsigmaTOFPr.value, lstarPidCircularCutPr.value, lstarPidPtRefPr.value); probe.nominal = passFinalCandidatePID(track.pt(), tpc[role], tof[role], true, lstarCutNsigmaTPCDe.value, lstarCutNsigmaTOFDe.value, lstarPidCircularCutDe.value, lstarPidPtRefDe.value, true) && - passOptionalDeuteronTOFExtras(track) && (lstarProxyUseProtonPIDAsDeuteron.value == 0 || passProton); + passOptionalDeuteronTOFExtras(track) && passDeuteronRefinements(track) && (lstarProxyUseProtonPIDAsDeuteron.value == 0 || passProton); } else { probe.nominal = passFinalCandidatePID(track.pt(), tpc[role], tof[role], true, lstarCutNsigmaTPCKaon.value, lstarCutNsigmaTOFKaon.value, lstarPidCircularCutKaon.value, lstarPidPtRefKaon.value); } @@ -1979,6 +2349,13 @@ struct Lambdastarproxy { const float eventMult = eventMultiplicityFT0MOrFallback(collision, eventMultFallback); // Exactly once after event cuts, before PID probes or candidate requirements. fillEventMultiplicityCounters(collision, eventMult); + const float spectraRawFT0M = fillTrackSpectraEvent(collision); + int spectraMCCollisionId = -1; + if constexpr (IsMC) { + if (trackSpectraConfig.lstarEnableTrackSpectra.value != 0) { + spectraMCCollisionId = collision.mcCollisionId(); + } + } if (lstarEnableProxyPIDProbes.value != 0) { if constexpr (IsMC) { fillProxyPIDProbes(tracks, eventMult, collision.mcCollisionId()); @@ -1999,6 +2376,13 @@ struct Lambdastarproxy { continue; } + if (deuteronPID.lstarEnableDeuteronPurityQA.value != 0) { + fillDeuteronPurityQA(trkD, DeuteronPIDStage::BeforePID); + if (std::isfinite(trkD.tpcNSigmaDe()) && std::abs(trkD.tpcNSigmaDe()) < lstarCutNsigmaTPCDe.value) { + fillDeuteronPurityQA(trkD, DeuteronPIDStage::AfterTPC); + } + } + // Deuteron kinematics needed before PID because the PID strategy can depend on pT const float ptD = trkD.pt(); const float etaD = trkD.eta(); @@ -2027,6 +2411,9 @@ struct Lambdastarproxy { if (!passesDeuteronSelection) { continue; } + if (!passDeuteronRefinements(trkD)) { + continue; + } const float nsTPCPrAsProxy = trkD.tpcNSigmaPr(); @@ -2051,6 +2438,9 @@ struct Lambdastarproxy { continue; } + fillDeuteronPurityQA(trkD, DeuteronPIDStage::Selected); + fillTrackSpectrum(trkD, TrackSpectrumSpecies::Deuteron, spectraRawFT0M, spectraMCCollisionId); + const double pD = static_cast(ptD) * std::cosh(static_cast(etaD)); // Candidate QA after final deuteron PID cut @@ -2179,6 +2569,7 @@ struct Lambdastarproxy { if (!isProton) { continue; } + fillTrackSpectrum(trkP, TrackSpectrumSpecies::Proton, spectraRawFT0M, spectraMCCollisionId); const double pP = static_cast(ptP) * std::cosh(static_cast(etaP)); @@ -2250,6 +2641,7 @@ struct Lambdastarproxy { if (!isKaon) { continue; } + fillTrackSpectrum(trkK, TrackSpectrumSpecies::Kaon, spectraRawFT0M, spectraMCCollisionId); const double pK = static_cast(ptK) * std::cosh(static_cast(etaK));