From 4771eb94b3f312e2553347b2983b81258087c325 Mon Sep 17 00:00:00 2001 From: ddobrigk Date: Tue, 7 Feb 2023 09:08:23 -0300 Subject: [PATCH 1/2] change subscription --- PWGLF/TableProducer/cascadelabelbuilder.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/PWGLF/TableProducer/cascadelabelbuilder.cxx b/PWGLF/TableProducer/cascadelabelbuilder.cxx index 0e3df9bfa57..3e5eb7899e0 100644 --- a/PWGLF/TableProducer/cascadelabelbuilder.cxx +++ b/PWGLF/TableProducer/cascadelabelbuilder.cxx @@ -50,7 +50,7 @@ struct cascadeLabelBuilder { //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* // build cascade labels - void process(aod::Collision const& collision, aod::CascDataExt const& casctable, aod::V0sLinked const&, aod::V0Datas const& v0table, LabeledTracks const&, aod::McParticles const&) + void process(aod::Collision const& collision, aod::CascDataExt const& casctable, aod::V0sLinked const&, aod::V0Datas const& v0table, aod::McTrackLabels const&, aod::McParticles const&) { for (auto& casc : casctable) { // Loop over those that actually have the corresponding V0 associated to them @@ -63,9 +63,9 @@ struct cascadeLabelBuilder { int lLabel = -1; // Acquire all three daughter tracks, please - auto lBachTrack = casc.bachelor_as(); - auto lNegTrack = v0data.negTrack_as(); - auto lPosTrack = v0data.posTrack_as(); + auto lBachTrack = casc.bachelor_as(); + auto lNegTrack = v0data.negTrack_as(); + auto lPosTrack = v0data.posTrack_as(); // Association check // There might be smarter ways of doing this in the future From 822c36b28072fa44716d23dc163f431c9860991c Mon Sep 17 00:00:00 2001 From: ddobrigk Date: Tue, 7 Feb 2023 09:08:53 -0300 Subject: [PATCH 2/2] change subscription --- .../TableProducer/lambdakzerolabelbuilder.cxx | 882 +++++++++++++++++- 1 file changed, 851 insertions(+), 31 deletions(-) diff --git a/PWGLF/TableProducer/lambdakzerolabelbuilder.cxx b/PWGLF/TableProducer/lambdakzerolabelbuilder.cxx index bb3b61ecb74..0a962b3d739 100644 --- a/PWGLF/TableProducer/lambdakzerolabelbuilder.cxx +++ b/PWGLF/TableProducer/lambdakzerolabelbuilder.cxx @@ -9,9 +9,27 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. // -// *+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* -// Lambdakzero label builder task -// *+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* +// *+-+*+-+*+-+*+-+*+-+*+-+* +// Lambdakzero builder task +// *+-+*+-+*+-+*+-+*+-+*+-+* +// +// This task loops over a set of V0 indices and +// creates the corresponding analysis tables that contain +// the typical information required for analysis. +// +// PERFORMANCE WARNING: this task includes several track +// propagation calls that are intrinsically heavy. Please +// also be cautious when adjusting selections: these can +// increase / decrease CPU consumption quite significantly. +// +// IDEAL USAGE: if you are interested in taking V0s and +// cascades and propagating TrackParCovs based on these, +// please do not re-propagate the daughters. Instead, +// the tables generated by this builder task can be used +// to instantiate a TrackPar object (default operation) +// or even a TrackParCov object (for which you will +// need to enable the option of producing the V0Cov and +// CascCov tables too). // // Comments, questions, complaints, suggestions? // Please write to: @@ -30,59 +48,861 @@ #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/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" using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; using std::array; +namespace o2::aod +{ +namespace v0tag +{ +// Global bool +DECLARE_SOA_COLUMN(IsInteresting, isInteresting, bool); //! will this be built or not? + +// MC association bools +DECLARE_SOA_COLUMN(IsTrueGamma, isTrueGamma, bool); //! PDG checked correctly in MC +DECLARE_SOA_COLUMN(IsTrueK0Short, isTrueK0Short, bool); //! PDG checked correctly in MC +DECLARE_SOA_COLUMN(IsTrueLambda, isTrueLambda, bool); //! PDG checked correctly in MC +DECLARE_SOA_COLUMN(IsTrueAntiLambda, isTrueAntiLambda, bool); //! PDG checked correctly in MC +DECLARE_SOA_COLUMN(IsTrueHypertriton, isTrueHypertriton, bool); //! PDG checked correctly in MC +DECLARE_SOA_COLUMN(IsTrueAntiHypertriton, isTrueAntiHypertriton, bool); //! PDG checked correctly in MC + +// dE/dx compatibility bools +DECLARE_SOA_COLUMN(IsGammaCandidate, isGammaCandidate, bool); //! compatible with dE/dx hypotheses +DECLARE_SOA_COLUMN(IsK0ShortCandidate, isK0ShortCandidate, bool); //! compatible with dE/dx hypotheses +DECLARE_SOA_COLUMN(IsLambdaCandidate, isLambdaCandidate, bool); //! compatible with dE/dx hypotheses +DECLARE_SOA_COLUMN(IsAntiLambdaCandidate, isAntiLambdaCandidate, bool); //! compatible with dE/dx hypotheses +DECLARE_SOA_COLUMN(IsHypertritonCandidate, isHypertritonCandidate, bool); //! compatible with dE/dx hypotheses +DECLARE_SOA_COLUMN(IsAntiHypertritonCandidate, isAntiHypertritonCandidate, bool); //! compatible with dE/dx hypotheses +} +DECLARE_SOA_TABLE(V0Tags, "AOD", "V0TAGS", + v0tag::IsInteresting, + v0tag::IsTrueGamma, + v0tag::IsTrueK0Short, + v0tag::IsTrueLambda, + v0tag::IsTrueAntiLambda, + v0tag::IsTrueHypertriton, + v0tag::IsTrueAntiHypertriton, + v0tag::IsGammaCandidate, + v0tag::IsK0ShortCandidate, + v0tag::IsLambdaCandidate, + v0tag::IsAntiLambdaCandidate, + v0tag::IsHypertritonCandidate, + v0tag::IsAntiHypertritonCandidate); +} // namespace o2::aod + +// use parameters + cov mat non-propagated, aux info + (extension propagated) +using FullTracksExt = soa::Join; +using FullTracksExtIU = soa::Join; +using TracksWithExtra = soa::Join; + +// For dE/dx association in pre-selection +using TracksExtraWithPID = soa::Join; + +// For MC and dE/dx association +using TracksExtraWithPIDandLabels = soa::Join; + +// Pre-selected V0s +using TaggedV0s = soa::Join; + // For MC association in pre-selection -using LabeledTracks = soa::Join; +using LabeledTracksExtra = soa::Join; + +struct lambdakzeroBuilder { + Produces v0data; + Produces v0covs; // covariances + Service ccdb; + + // Configurables related to table creation + Configurable createV0CovMats{"createV0CovMats", -1, {"Produces V0 cov matrices. -1: auto, 0: don't, 1: yes. Default: auto (-1)"}}; + + // use auto-detect configuration + Configurable d_UseAutodetectMode{"d_UseAutodetectMode", true, "Autodetect requested topo sels"}; + + Configurable dcanegtopv{"dcanegtopv", .1, "DCA Neg To PV"}; + Configurable dcapostopv{"dcapostopv", .1, "DCA Pos To PV"}; + Configurable v0cospa{"v0cospa", 0.995, "V0 CosPA"}; // double -> N.B. dcos(x)/dx = 0 at x=0) + Configurable dcav0dau{"dcav0dau", 1.0, "DCA V0 Daughters"}; + Configurable v0radius{"v0radius", 0.9, "v0radius"}; + + Configurable tpcrefit{"tpcrefit", 0, "demand TPC refit"}; + + // Operation and minimisation criteria + Configurable d_bz_input{"d_bz", -999, "bz field, -999 is automatic"}; + Configurable d_UseAbsDCA{"d_UseAbsDCA", true, "Use Abs DCAs"}; + Configurable d_UseWeightedPCA{"d_UseWeightedPCA", false, "Vertices use cov matrices"}; + Configurable useMatCorrType{"useMatCorrType", 0, "0: none, 1: TGeo, 2: LUT"}; + Configurable rejDiffCollTracks{"rejDiffCollTracks", 0, "rejDiffCollTracks"}; + Configurable d_doTrackQA{"d_doTrackQA", false, "do track QA"}; + + // 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"}; + + int mRunNumber; + float d_bz; + float maxSnp; // max sine phi for propagation + float maxStep; // max step size (cm) for propagation + o2::base::MatLayerCylSet* lut = nullptr; + + // Define o2 fitter, 2-prong, active memory (no need to redefine per event) + o2::vertexing::DCAFitterN<2> fitter; + + Filter taggedFilter = aod::v0tag::isInteresting == true; + + // For manual sliceBy + Preslice perCollision = o2::aod::v0::collisionId; + + enum v0step { kV0All = 0, + kV0TPCrefit, + kV0DCAxy, + kV0DCADau, + kV0CosPA, + kV0Radius, + kNV0Steps }; + + // Helper struct to pass V0 information + struct { + float posTrackX; + float negTrackX; + std::array pos; + std::array posP; + std::array negP; + float dcaV0dau; + float posDCAxy; + float negDCAxy; + float cosPA; + float V0radius; + float lambdaMass; + float antilambdaMass; + } v0candidate; + + // Helper struct to do bookkeeping of building parameters + struct { + std::array v0stats; + std::array posITSclu; + std::array negITSclu; + long exceptions; + long eventCounter; + } statisticsRegistry; + + HistogramRegistry registry{ + "registry", + {{"hEventCounter", "hEventCounter", {HistType::kTH1F, {{1, 0.0f, 1.0f}}}}, + {"hCaughtExceptions", "hCaughtExceptions", {HistType::kTH1F, {{1, 0.0f, 1.0f}}}}, + {"hPositiveITSClusters", "hPositiveITSClusters", {HistType::kTH1F, {{10, -0.5f, 9.5f}}}}, + {"hNegativeITSClusters", "hNegativeITSClusters", {HistType::kTH1F, {{10, -0.5f, 9.5f}}}}, + {"hV0Criteria", "hV0Criteria", {HistType::kTH1F, {{10, -0.5f, 9.5f}}}}}}; + + void resetHistos() + { + statisticsRegistry.exceptions = 0; + statisticsRegistry.eventCounter = 0; + for (Int_t ii = 0; ii < kNV0Steps; ii++) + statisticsRegistry.v0stats[ii] = 0; + for (Int_t ii = 0; ii < 10; ii++) { + statisticsRegistry.posITSclu[ii] = 0; + statisticsRegistry.negITSclu[ii] = 0; + } + } + + void fillHistos() + { + registry.fill(HIST("hEventCounter"), 0.0, statisticsRegistry.eventCounter); + registry.fill(HIST("hCaughtExceptions"), 0.0, statisticsRegistry.exceptions); + for (Int_t ii = 0; ii < kNV0Steps; ii++) + registry.fill(HIST("hV0Criteria"), ii, statisticsRegistry.v0stats[ii]); + if (d_doTrackQA) { + for (Int_t ii = 0; ii < 10; ii++) { + registry.fill(HIST("hPositiveITSClusters"), ii, statisticsRegistry.posITSclu[ii]); + registry.fill(HIST("hNegativeITSClusters"), ii, statisticsRegistry.negITSclu[ii]); + } + } + } + + o2::track::TrackParCov lPositiveTrack; + o2::track::TrackParCov lNegativeTrack; + + void init(InitContext& context) + { + resetHistos(); + + 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); + + if (useMatCorrType == 1) { + LOGF(info, "TGeo correction requested, loading geometry"); + if (!o2::base::GeometryManager::isGeometryLoaded()) { + ccdb->get(geoPath); + } + } + if (useMatCorrType == 2) { + LOGF(info, "LUT correction requested, loading LUT"); + lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get(lutPath)); + } + + if (doprocessRun2 == false && doprocessRun3 == false) { + LOGF(fatal, "Neither processRun2 nor processRun3 enabled. Please choose one."); + } + if (doprocessRun2 == true && doprocessRun3 == true) { + LOGF(fatal, "Cannot enable processRun2 and processRun3 at the same time. Please choose one."); + } + + if (d_UseAutodetectMode) { + double loosest_v0cospa = 100; + float loosest_dcav0dau = -100; + float loosest_dcapostopv = 100; + float loosest_dcanegtopv = 100; + float loosest_radius = 100; + + double detected_v0cospa = -100; + float detected_dcav0dau = -100; + float detected_dcapostopv = 100; + float detected_dcanegtopv = 100; + float detected_radius = 100; + + LOGF(info, "*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*"); + LOGF(info, " Single-strange builder self-configuration"); + LOGF(info, "*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*"); + auto& workflows = context.services().get(); + for (DeviceSpec const& device : workflows.devices) { + // Step 1: check if this device subscribed to the V0data table + for (auto const& input : device.inputs) { + if (device.name.compare("lambdakzero-initializer") == 0) + continue; // don't listen to the initializer, it's just to extend stuff + const std::string v0DataName = "V0Datas"; + if (input.matcher.binding == v0DataName && device.name.compare("multistrange-builder") != 0) { + LOGF(info, "Device named %s has subscribed to V0datas table! Will now scan for desired settings...", device.name); + for (auto const& option : device.options) { + // 5 V0 topological selections + if (option.name.compare("v0setting_cospa") == 0) { + detected_v0cospa = option.defaultValue.get(); + LOGF(info, "%s requested V0 cospa = %f", device.name, detected_v0cospa); + if (detected_v0cospa < loosest_v0cospa) + loosest_v0cospa = detected_v0cospa; + } + if (option.name.compare("v0setting_dcav0dau") == 0) { + detected_dcav0dau = option.defaultValue.get(); + LOGF(info, "%s requested DCA V0 daughters = %f", device.name, detected_dcav0dau); + if (detected_dcav0dau > loosest_dcav0dau) + loosest_dcav0dau = detected_dcav0dau; + } + if (option.name.compare("v0setting_dcapostopv") == 0) { + detected_dcapostopv = option.defaultValue.get(); + LOGF(info, "%s requested DCA positive daughter to PV = %f", device.name, detected_dcapostopv); + if (detected_dcapostopv < loosest_dcapostopv) + loosest_dcapostopv = detected_dcapostopv; + } + if (option.name.compare("v0setting_dcanegtopv") == 0) { + detected_dcanegtopv = option.defaultValue.get(); + LOGF(info, "%s requested DCA negative daughter to PV = %f", device.name, detected_dcanegtopv); + if (detected_dcanegtopv < loosest_dcanegtopv) + loosest_dcanegtopv = detected_dcanegtopv; + } + if (option.name.compare("v0setting_radius") == 0) { + detected_radius = option.defaultValue.get(); + LOGF(info, "%s requested minimum V0 radius = %f", device.name, detected_radius); + if (detected_radius < loosest_radius) + loosest_radius = detected_radius; + } + } + } + const std::string V0CovsName = "V0Covs"; + if (input.matcher.binding == V0CovsName) { + LOGF(info, "Device named %s has subscribed to V0Covs table! Enabling.", device.name); + createV0CovMats.value = 1; + } + } + } + LOGF(info, "Self-configuration finished! Decided on selections:"); + LOGF(info, " -+*> V0 cospa ..............: %.6f", loosest_v0cospa); + LOGF(info, " -+*> DCA V0 daughters ......: %.6f", loosest_dcav0dau); + LOGF(info, " -+*> DCA positive daughter .: %.6f", loosest_dcapostopv); + LOGF(info, " -+*> DCA negative daughter .: %.6f", loosest_dcanegtopv); + LOGF(info, " -+*> Minimum V0 radius .....: %.6f", loosest_radius); + + dcanegtopv.value = loosest_dcanegtopv; + dcapostopv.value = loosest_dcapostopv; + v0cospa.value = loosest_v0cospa; + dcav0dau.value = loosest_dcav0dau; + v0radius.value = loosest_radius; + } + + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + LOGF(info, " -+*> process call configuration:"); + if (doprocessRun2 == true) { + LOGF(info, " ---+*> Run 2 processing enabled. Will subscribe to Tracks table."); + }; + if (doprocessRun3 == true) { + LOGF(info, " ---+*> Run 3 processing enabled. Will subscribe to TracksIU table."); + }; + if (createV0CovMats > 0) { + LOGF(info, " ---+*> Will produce V0 cov mat table"); + }; + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + + // initialize O2 2-prong fitter (only once) + fitter.setPropagateToPCA(true); + fitter.setMaxR(200.); + fitter.setMinParamChange(1e-3); + fitter.setMinRelChi2Change(0.9); + fitter.setMaxDZIni(1e9); + fitter.setMaxChi2(1e9); + fitter.setUseAbsDCA(d_UseAbsDCA); + fitter.setWeightedFinalPCA(d_UseWeightedPCA); + + // Material correction in the DCA fitter + o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE; + if (useMatCorrType == 1) + matCorr = o2::base::Propagator::MatCorrType::USEMatCorrTGeo; + if (useMatCorrType == 2) + matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT; + fitter.setMatCorrType(matCorr); + } + + void initCCDB(aod::BCsWithTimestamps::iterator const& bc) + { + if (mRunNumber == bc.runNumber()) { + return; + } + + // In case override, don't proceed, please - no CCDB access required + if (d_bz_input > -990) { + d_bz = d_bz_input; + fitter.setBz(d_bz); + o2::parameters::GRPMagField grpmag; + if (fabs(d_bz) > 1e-5) { + grpmag.setL3Current(30000.f / (d_bz / 5.0f)); + } + o2::base::Propagator::initFieldFromGRP(&grpmag); + mRunNumber = bc.runNumber(); + return; + } + + auto run3grp_timestamp = bc.timestamp(); + o2::parameters::GRPObject* grpo = ccdb->getForTimeStamp(grpPath, run3grp_timestamp); + 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 timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG"; + } else { + grpmag = ccdb->getForTimeStamp(grpmagPath, run3grp_timestamp); + if (!grpmag) { + LOG(fatal) << "Got nullptr from CCDB for path " << grpmagPath << " of object GRPMagField and " << grpPath << " of object GRPObject for timestamp " << run3grp_timestamp; + } + 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 timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG"; + } + mRunNumber = bc.runNumber(); + // Set magnetic field value once known + fitter.setBz(d_bz); + + if (useMatCorrType == 2) { + // setMatLUT only after magfield has been initalized + // (setMatLUT has implicit and problematic init field call if not) + o2::base::Propagator::Instance()->setMatLUT(lut); + } + } + + template + bool buildV0Candidate(TV0Object const& V0) + { + // Get tracks + auto const& posTrack = V0.template posTrack_as(); + auto const& negTrack = V0.template negTrack_as(); + auto const& collision = V0.collision(); + + // value 0.5: any considered V0 + statisticsRegistry.v0stats[kV0All]++; + if (tpcrefit) { + if (!(posTrack.trackType() & o2::aod::track::TPCrefit)) { + return false; + } + if (!(negTrack.trackType() & o2::aod::track::TPCrefit)) { + return false; + } + } + + // Passes TPC refit + statisticsRegistry.v0stats[kV0TPCrefit]++; + + // Calculate DCA with respect to the collision associated to the V0, not individual tracks + gpu::gpustd::array dcaInfo; + + auto posTrackPar = getTrackPar(posTrack); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, posTrackPar, 2.f, fitter.getMatCorrType(), &dcaInfo); + auto posTrackdcaXY = dcaInfo[0]; + + auto negTrackPar = getTrackPar(negTrack); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, negTrackPar, 2.f, fitter.getMatCorrType(), &dcaInfo); + auto negTrackdcaXY = dcaInfo[0]; + + if (fabs(posTrackdcaXY) < dcapostopv || fabs(negTrackdcaXY) < dcanegtopv) { + return false; + } + + // Initialize properly, please + v0candidate.posDCAxy = posTrackdcaXY; + v0candidate.negDCAxy = negTrackdcaXY; + + // passes DCAxy + statisticsRegistry.v0stats[kV0DCAxy]++; + + // Change strangenessBuilder tracks + lPositiveTrack = getTrackParCov(posTrack); + lNegativeTrack = getTrackParCov(negTrack); + + //---/---/---/ + // Move close to minima + int nCand = 0; + try { + nCand = fitter.process(lPositiveTrack, lNegativeTrack); + } catch (...) { + statisticsRegistry.exceptions++; + LOG(error) << "Exception caught in DCA fitter process call!"; + return false; + } + if (nCand == 0) { + return false; + } + + v0candidate.posTrackX = fitter.getTrack(0).getX(); + v0candidate.negTrackX = fitter.getTrack(1).getX(); + + lPositiveTrack = fitter.getTrack(0); + lNegativeTrack = fitter.getTrack(1); + lPositiveTrack.getPxPyPzGlo(v0candidate.posP); + lNegativeTrack.getPxPyPzGlo(v0candidate.negP); + + // get decay vertex coordinates + const auto& vtx = fitter.getPCACandidate(); + for (int i = 0; i < 3; i++) { + v0candidate.pos[i] = vtx[i]; + } + + v0candidate.dcaV0dau = TMath::Sqrt(fitter.getChi2AtPCACandidate()); + + // Apply selections so a skimmed table is created only + if (v0candidate.dcaV0dau > dcav0dau) { + return false; + } + + // Passes DCA between daughters check + statisticsRegistry.v0stats[kV0DCADau]++; + + v0candidate.cosPA = RecoDecay::cpa(array{collision.posX(), collision.posY(), collision.posZ()}, array{v0candidate.pos[0], v0candidate.pos[1], v0candidate.pos[2]}, array{v0candidate.posP[0] + v0candidate.negP[0], v0candidate.posP[1] + v0candidate.negP[1], v0candidate.posP[2] + v0candidate.negP[2]}); + if (v0candidate.cosPA < v0cospa) { + return false; + } + + // Passes CosPA check + statisticsRegistry.v0stats[kV0CosPA]++; + + v0candidate.V0radius = RecoDecay::sqrtSumOfSquares(v0candidate.pos[0], v0candidate.pos[1]); + if (v0candidate.V0radius < v0radius) { + return false; + } + + // Passes radius check + statisticsRegistry.v0stats[kV0Radius]++; + // Return OK: passed all v0 candidate selecton criteria + if (d_doTrackQA) { + if (posTrack.itsNCls() < 10) + statisticsRegistry.posITSclu[posTrack.itsNCls()]++; + if (negTrack.itsNCls() < 10) + statisticsRegistry.negITSclu[negTrack.itsNCls()]++; + } + return true; + } + + template + void buildStrangenessTables(TV0Table const& V0s) + { + statisticsRegistry.eventCounter++; + + // Loops over all V0s in the time frame + for (auto& V0 : V0s) { + // populates v0candidate struct declared inside strangenessbuilder + bool validCandidate = buildV0Candidate(V0); + + if (!validCandidate) { + continue; // doesn't pass selections + } + + // populates table for V0 analysis + v0data(V0.posTrackId(), + V0.negTrackId(), + V0.collisionId(), + V0.globalIndex(), + v0candidate.posTrackX, v0candidate.negTrackX, + v0candidate.pos[0], v0candidate.pos[1], v0candidate.pos[2], + v0candidate.posP[0], v0candidate.posP[1], v0candidate.posP[2], + v0candidate.negP[0], v0candidate.negP[1], v0candidate.negP[2], + v0candidate.dcaV0dau, + v0candidate.posDCAxy, + v0candidate.negDCAxy); + + // populate V0 covariance matrices if required by any other task + if (createV0CovMats) { + // Calculate position covariance matrix + auto covVtxV = fitter.calcPCACovMatrix(0); + // std::array positionCovariance; + float positionCovariance[6]; + positionCovariance[0] = covVtxV(0, 0); + positionCovariance[1] = covVtxV(1, 0); + positionCovariance[2] = covVtxV(1, 1); + positionCovariance[3] = covVtxV(2, 0); + positionCovariance[4] = covVtxV(2, 1); + positionCovariance[5] = covVtxV(2, 2); + // store momentum covariance matrix + std::array covTpositive = {0.}; + std::array covTnegative = {0.}; + // std::array momentumCovariance; + float momentumCovariance[6]; + lPositiveTrack.getCovXYZPxPyPzGlo(covTpositive); + lNegativeTrack.getCovXYZPxPyPzGlo(covTnegative); + constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component + for (int i = 0; i < 6; i++) { + momentumCovariance[i] = covTpositive[MomInd[i]] + covTnegative[MomInd[i]]; + } + v0covs(positionCovariance, momentumCovariance); + } + } + // En masse histo filling at end of process call + fillHistos(); + resetHistos(); + } + + void processRun2(aod::Collisions const& collisions, soa::Filtered const& V0s, FullTracksExt const&, aod::BCsWithTimestamps const&) + { + for (const auto& collision : collisions) { + // Fire up CCDB + auto bc = collision.bc_as(); + initCCDB(bc); + // Do analysis with collision-grouped V0s, retain full collision information + const uint64_t collIdx = collision.globalIndex(); + auto V0Table_thisCollision = V0s.sliceBy(perCollision, collIdx); + buildStrangenessTables(V0Table_thisCollision); + } + } + PROCESS_SWITCH(lambdakzeroBuilder, processRun2, "Produce Run 2 V0 tables", false); + + void processRun3(aod::Collisions const& collisions, soa::Filtered const& V0s, FullTracksExtIU const&, aod::BCsWithTimestamps const&) + { + for (const auto& collision : collisions) { + // Fire up CCDB + auto bc = collision.bc_as(); + initCCDB(bc); + // Do analysis with collision-grouped V0s, retain full collision information + const uint64_t collIdx = collision.globalIndex(); + auto V0Table_thisCollision = V0s.sliceBy(perCollision, collIdx); + buildStrangenessTables(V0Table_thisCollision); + } + } + PROCESS_SWITCH(lambdakzeroBuilder, processRun3, "Produce Run 3 V0 tables", true); +}; //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* -struct lambdakzeroLabelBuilder { - Produces v0labels; // MC labels for V0s +struct lambdakzeroPreselector { + Produces v0tags; // MC tags + + Configurable dIfMCgenerateK0Short{"dIfMCgenerateK0Short", true, "if MC, generate MC true K0Short (yes/no)"}; + Configurable dIfMCgenerateLambda{"dIfMCgenerateLambda", true, "if MC, generate MC true Lambda (yes/no)"}; + Configurable dIfMCgenerateAntiLambda{"dIfMCgenerateAntiLambda", true, "if MC, generate MC true AntiLambda (yes/no)"}; + Configurable dIfMCgenerateGamma{"dIfMCgenerateGamma", false, "if MC, generate MC true gamma (yes/no)"}; + Configurable dIfMCgenerateHypertriton{"dIfMCgenerateHypertriton", false, "if MC, generate MC true hypertritons (yes/no)"}; + Configurable dIfMCgenerateAntiHypertriton{"dIfMCgenerateAntiHypertriton", false, "if MC, generate MC true antihypertritons (yes/no)"}; + + Configurable ddEdxPreSelectK0Short{"ddEdxPreSelectK0Short", true, "pre-select dE/dx compatibility with K0Short (yes/no)"}; + Configurable ddEdxPreSelectLambda{"ddEdxPreSelectLambda", true, "pre-select dE/dx compatibility with Lambda (yes/no)"}; + Configurable ddEdxPreSelectAntiLambda{"ddEdxPreSelectAntiLambda", true, "pre-select dE/dx compatibility with AntiLambda (yes/no)"}; + Configurable ddEdxPreSelectGamma{"ddEdxPreSelectGamma", false, "pre-select dE/dx compatibility with gamma (yes/no)"}; + Configurable ddEdxPreSelectHypertriton{"ddEdxPreSelectHypertriton", false, "pre-select dE/dx compatibility with hypertritons (yes/no)"}; + Configurable ddEdxPreSelectAntiHypertriton{"ddEdxPreSelectAntiHypertriton", false, "pre-select dE/dx compatibility with antihypertritons (yes/no)"}; + + // dEdx pre-selection compatibility + Configurable ddEdxPreSelectionWindow{"ddEdxPreSelectionWindow", 7, "Nsigma window for dE/dx preselection"}; + + // tpc quality pre-selection + Configurable dTPCNCrossedRows{"dTPCNCrossedRows", 50, "Minimum TPC crossed rows"}; + + // context-aware selections + Configurable dPreselectOnlyBaryons{"dPreselectOnlyBaryons", false, "apply TPC dE/dx and quality only to baryon daughters"}; + void init(InitContext const&) {} //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* - // build V0 labels - void process(aod::Collision const& collision, aod::V0Datas const& v0table, LabeledTracks const&, aod::McParticles const& particlesMC) + /// function to check track quality + template + void checkTrackQuality(TV0Object const& lV0Candidate, bool& lIsInteresting, bool lIsGamma, bool lIsK0Short, bool lIsLambda, bool lIsAntiLambda, bool lIsHypertriton, bool lIsAntiHypertriton) { - for (auto& v0 : v0table) { - int lLabel = -1; - - auto lNegTrack = v0.negTrack_as(); - auto lPosTrack = v0.posTrack_as(); - - // Association check - // There might be smarter ways of doing this in the future - if (lNegTrack.has_mcParticle() && lPosTrack.has_mcParticle()) { - auto lMCNegTrack = lNegTrack.mcParticle_as(); - auto lMCPosTrack = lPosTrack.mcParticle_as(); - if (lMCNegTrack.has_mothers() && lMCPosTrack.has_mothers()) { - - for (auto& lNegMother : lMCNegTrack.mothers_as()) { - for (auto& lPosMother : lMCPosTrack.mothers_as()) { - if (lNegMother.globalIndex() == lPosMother.globalIndex()) { - lLabel = lNegMother.globalIndex(); - } + lIsInteresting = false; + auto lNegTrack = lV0Candidate.template negTrack_as(); + auto lPosTrack = lV0Candidate.template posTrack_as(); + + // No baryons in decay + if ((lIsGamma || lIsK0Short) && (lPosTrack.tpcNClsCrossedRows() >= dTPCNCrossedRows && lNegTrack.tpcNClsCrossedRows() >= dTPCNCrossedRows)) + lIsInteresting = true; + // With baryons in decay + if (lIsLambda && (lPosTrack.tpcNClsCrossedRows() >= dTPCNCrossedRows && (lNegTrack.tpcNClsCrossedRows() >= dTPCNCrossedRows || dPreselectOnlyBaryons))) + lIsInteresting = true; + if (lIsAntiLambda && (lNegTrack.tpcNClsCrossedRows() >= dTPCNCrossedRows && (lPosTrack.tpcNClsCrossedRows() >= dTPCNCrossedRows || dPreselectOnlyBaryons))) + lIsInteresting = true; + if (lIsHypertriton && (lPosTrack.tpcNClsCrossedRows() >= dTPCNCrossedRows && (lNegTrack.tpcNClsCrossedRows() >= dTPCNCrossedRows || dPreselectOnlyBaryons))) + lIsInteresting = true; + if (lIsHypertriton && (lNegTrack.tpcNClsCrossedRows() >= dTPCNCrossedRows && (lPosTrack.tpcNClsCrossedRows() >= dTPCNCrossedRows || dPreselectOnlyBaryons))) + lIsInteresting = true; + } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + /// function to check PDG association + template + void checkPDG(TV0Object const& lV0Candidate, bool& lIsInteresting, bool& lIsGamma, bool& lIsK0Short, bool& lIsLambda, bool& lIsAntiLambda, bool& lIsHypertriton, bool& lIsAntiHypertriton) + { + int lPDG = -1; + auto lNegTrack = lV0Candidate.template negTrack_as(); + auto lPosTrack = lV0Candidate.template posTrack_as(); + + // Association check + // There might be smarter ways of doing this in the future + if (lNegTrack.has_mcParticle() && lPosTrack.has_mcParticle()) { + auto lMCNegTrack = lNegTrack.template mcParticle_as(); + auto lMCPosTrack = lPosTrack.template mcParticle_as(); + if (lMCNegTrack.has_mothers() && lMCPosTrack.has_mothers()) { + + for (auto& lNegMother : lMCNegTrack.template mothers_as()) { + for (auto& lPosMother : lMCPosTrack.template mothers_as()) { + if (lNegMother.globalIndex() == lPosMother.globalIndex()) { + lPDG = lNegMother.pdgCode(); } } } - } // end association check - // Construct label table (note: this will be joinable with V0Datas!) - v0labels( - lLabel); + } + } // end association check + if (lPDG == 310 && dIfMCgenerateK0Short) { + lIsK0Short = true; + lIsInteresting = true; + } + if (lPDG == 3122 && dIfMCgenerateLambda) { + lIsLambda = true; + lIsInteresting = true; + } + if (lPDG == -3122 && dIfMCgenerateAntiLambda) { + lIsAntiLambda = true; + lIsInteresting = true; + } + if (lPDG == 22 && dIfMCgenerateGamma) { + lIsGamma = true; + lIsInteresting = true; + } + if (lPDG == 1010010030 && dIfMCgenerateHypertriton) { + lIsHypertriton = true; + lIsInteresting = true; + } + if (lPDG == -1010010030 && dIfMCgenerateAntiHypertriton) { + lIsAntiHypertriton = true; + lIsInteresting = true; + } + } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + template + void checkdEdx(TV0Object const& lV0Candidate, bool& lIsInteresting, bool& lIsGamma, bool& lIsK0Short, bool& lIsLambda, bool& lIsAntiLambda, bool& lIsHypertriton, bool& lIsAntiHypertriton) + { + auto lNegTrack = lV0Candidate.template negTrack_as(); + auto lPosTrack = lV0Candidate.template posTrack_as(); + + // dEdx check with LF PID + if (TMath::Abs(lNegTrack.tpcNSigmaEl()) < ddEdxPreSelectionWindow && + TMath::Abs(lPosTrack.tpcNSigmaEl()) < ddEdxPreSelectionWindow && + ddEdxPreSelectGamma) { + lIsGamma = 1; + lIsInteresting = 1; + } + if (TMath::Abs(lNegTrack.tpcNSigmaPi()) < ddEdxPreSelectionWindow && + TMath::Abs(lPosTrack.tpcNSigmaPi()) < ddEdxPreSelectionWindow && + ddEdxPreSelectK0Short) { + lIsK0Short = 1; + lIsInteresting = 1; + } + if ((TMath::Abs(lNegTrack.tpcNSigmaPi()) < ddEdxPreSelectionWindow || dPreselectOnlyBaryons) && + TMath::Abs(lPosTrack.tpcNSigmaPr()) < ddEdxPreSelectionWindow && + ddEdxPreSelectLambda) { + lIsLambda = 1; + lIsInteresting = 1; + } + if (TMath::Abs(lNegTrack.tpcNSigmaPr()) < ddEdxPreSelectionWindow && + (TMath::Abs(lPosTrack.tpcNSigmaPi()) < ddEdxPreSelectionWindow || dPreselectOnlyBaryons) && + ddEdxPreSelectAntiLambda) { + lIsAntiLambda = 1; + lIsInteresting = 1; + } + if (TMath::Abs(lNegTrack.tpcNSigmaPi()) < ddEdxPreSelectionWindow && + (TMath::Abs(lPosTrack.tpcNSigmaHe()) < ddEdxPreSelectionWindow || dPreselectOnlyBaryons) && + ddEdxPreSelectHypertriton) { + lIsHypertriton = 1; + lIsInteresting = 1; + } + if ((TMath::Abs(lNegTrack.tpcNSigmaHe()) < ddEdxPreSelectionWindow || dPreselectOnlyBaryons) && + TMath::Abs(lPosTrack.tpcNSigmaPi()) < ddEdxPreSelectionWindow && + ddEdxPreSelectAntiHypertriton) { + lIsAntiHypertriton = 1; + lIsInteresting = 1; + } + } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + /// This process function ensures that all V0s are built. It will simply tag everything as true. + void processBuildAll(aod::V0s const& v0table, aod::TracksExtra const&) + { + for (auto& v0 : v0table) { + bool lIsQualityInteresting = false; + checkTrackQuality(v0, lIsQualityInteresting, true, true, true, true, true, true); + v0tags(lIsQualityInteresting, + true, true, true, true, true, true, + true, true, true, true, true, true); + } + } + PROCESS_SWITCH(lambdakzeroPreselector, processBuildAll, "Switch to build all V0s", true); + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + void processBuildMCAssociated(aod::Collision const& collision, aod::V0s const& v0table, LabeledTracksExtra const&, aod::McParticles const& particlesMC) + { + for (auto& v0 : v0table) { + bool lIsInteresting = false; + bool lIsTrueGamma = false; + bool lIsTrueK0Short = false; + bool lIsTrueLambda = false; + bool lIsTrueAntiLambda = false; + bool lIsTrueHypertriton = false; + bool lIsTrueAntiHypertriton = false; + + bool lIsQualityInteresting = false; + + checkPDG(v0, lIsInteresting, lIsTrueGamma, lIsTrueK0Short, lIsTrueLambda, lIsTrueAntiLambda, lIsTrueHypertriton, lIsTrueAntiHypertriton); + checkTrackQuality(v0, lIsQualityInteresting, true, true, true, true, true, true); + v0tags(lIsInteresting * lIsQualityInteresting, + lIsTrueGamma, lIsTrueK0Short, lIsTrueLambda, lIsTrueAntiLambda, lIsTrueHypertriton, lIsTrueAntiHypertriton, + true, true, true, true, true, true); + } + } + PROCESS_SWITCH(lambdakzeroPreselector, processBuildMCAssociated, "Switch to build MC-associated V0s", false); + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + void processBuildValiddEdx(aod::Collision const& collision, aod::V0s const& v0table, TracksExtraWithPID const&) + { + for (auto& v0 : v0table) { + bool lIsInteresting = false; + bool lIsdEdxGamma = false; + bool lIsdEdxK0Short = false; + bool lIsdEdxLambda = false; + bool lIsdEdxAntiLambda = false; + bool lIsdEdxHypertriton = false; + bool lIsdEdxAntiHypertriton = false; + + bool lIsQualityInteresting = false; + + checkdEdx(v0, lIsInteresting, lIsdEdxGamma, lIsdEdxK0Short, lIsdEdxLambda, lIsdEdxAntiLambda, lIsdEdxHypertriton, lIsdEdxAntiHypertriton); + checkTrackQuality(v0, lIsQualityInteresting, lIsdEdxGamma, lIsdEdxK0Short, lIsdEdxLambda, lIsdEdxAntiLambda, lIsdEdxHypertriton, lIsdEdxAntiHypertriton); + v0tags(lIsInteresting * lIsQualityInteresting, + true, true, true, true, true, true, + lIsdEdxGamma, lIsdEdxK0Short, lIsdEdxLambda, lIsdEdxAntiLambda, lIsdEdxHypertriton, lIsdEdxAntiHypertriton); } } + PROCESS_SWITCH(lambdakzeroPreselector, processBuildValiddEdx, "Switch to build V0s with dE/dx preselection", false); + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + void processBuildValiddEdxMCAssociated(aod::Collision const& collision, aod::V0s const& v0table, TracksExtraWithPIDandLabels const&) + { + for (auto& v0 : v0table) { + bool lIsTrueInteresting = false; + bool lIsTrueGamma = false; + bool lIsTrueK0Short = false; + bool lIsTrueLambda = false; + bool lIsTrueAntiLambda = false; + bool lIsTrueHypertriton = false; + bool lIsTrueAntiHypertriton = false; + + bool lIsdEdxInteresting = false; + bool lIsdEdxGamma = false; + bool lIsdEdxK0Short = false; + bool lIsdEdxLambda = false; + bool lIsdEdxAntiLambda = false; + bool lIsdEdxHypertriton = false; + bool lIsdEdxAntiHypertriton = false; + + bool lIsQualityInteresting = false; + + checkPDG(v0, lIsTrueInteresting, lIsTrueGamma, lIsTrueK0Short, lIsTrueLambda, lIsTrueAntiLambda, lIsTrueHypertriton, lIsTrueAntiHypertriton); + checkdEdx(v0, lIsdEdxInteresting, lIsdEdxGamma, lIsdEdxK0Short, lIsdEdxLambda, lIsdEdxAntiLambda, lIsdEdxHypertriton, lIsdEdxAntiHypertriton); + checkTrackQuality(v0, lIsQualityInteresting, lIsdEdxGamma, lIsdEdxK0Short, lIsdEdxLambda, lIsdEdxAntiLambda, lIsdEdxHypertriton, lIsdEdxAntiHypertriton); + v0tags(lIsTrueInteresting * lIsdEdxInteresting * lIsQualityInteresting, + lIsTrueGamma, lIsTrueK0Short, lIsTrueLambda, lIsTrueAntiLambda, lIsTrueHypertriton, lIsTrueAntiHypertriton, + lIsdEdxGamma, lIsdEdxK0Short, lIsdEdxLambda, lIsdEdxAntiLambda, lIsdEdxHypertriton, lIsdEdxAntiHypertriton); + } + } + PROCESS_SWITCH(lambdakzeroPreselector, processBuildValiddEdxMCAssociated, "Switch to build MC-associated V0s with dE/dx preselection", false); + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* +}; + +//*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* +struct lambdakzeroV0DataLinkBuilder { + Produces v0dataLink; + + void init(InitContext const&) {} + + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + // build V0 -> V0Data link table + void process(aod::V0s const& v0table, aod::V0Datas const& v0datatable) + { + std::vector lIndices; + lIndices.reserve(v0table.size()); + for (int ii = 0; ii < v0table.size(); ii++) + lIndices[ii] = -1; + for (auto& v0data : v0datatable) { + lIndices[v0data.v0Id()] = v0data.globalIndex(); + } + for (int ii = 0; ii < v0table.size(); ii++) { + v0dataLink(lIndices[ii]); + } + } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* +}; + +// Extends the v0data table with expression columns +struct lambdakzeroInitializer { + Spawns v0datas; + void init(InitContext const&) {} }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc)}; }