From f7f9b769958ab74019a37ae3b6fbda46d55551d5 Mon Sep 17 00:00:00 2001 From: David Dobrigkeit Chinellato Date: Fri, 18 Jul 2025 13:31:06 -0300 Subject: [PATCH 01/12] [PWGLF/Common] Add test version of propagation service task --- .../TableProducer/trackPropagationTester.cxx | 2 +- .../TableProducer/Strangeness/CMakeLists.txt | 5 + .../Strangeness/propagationService.cxx | 146 + PWGLF/Utils/strangenessBuilderModule.h | 2536 +++++++++++++++++ 4 files changed, 2688 insertions(+), 1 deletion(-) create mode 100644 PWGLF/TableProducer/Strangeness/propagationService.cxx create mode 100644 PWGLF/Utils/strangenessBuilderModule.h diff --git a/Common/TableProducer/trackPropagationTester.cxx b/Common/TableProducer/trackPropagationTester.cxx index 694ff77d3c0..817c3f59db5 100644 --- a/Common/TableProducer/trackPropagationTester.cxx +++ b/Common/TableProducer/trackPropagationTester.cxx @@ -105,4 +105,4 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { WorkflowSpec workflow{adaptAnalysisTask(cfgc)}; return workflow; -} +} \ No newline at end of file diff --git a/PWGLF/TableProducer/Strangeness/CMakeLists.txt b/PWGLF/TableProducer/Strangeness/CMakeLists.txt index 55bf4550601..bc8ba1da673 100644 --- a/PWGLF/TableProducer/Strangeness/CMakeLists.txt +++ b/PWGLF/TableProducer/Strangeness/CMakeLists.txt @@ -166,3 +166,8 @@ o2physics_add_dpl_workflow(lambdaspincorrelation SOURCES lambdaspincorrelation.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(propagationservice + SOURCES propagationService.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter KFParticle::KFParticle O2Physics::TPCDriftManager + COMPONENT_NAME Analysis) diff --git a/PWGLF/TableProducer/Strangeness/propagationService.cxx b/PWGLF/TableProducer/Strangeness/propagationService.cxx new file mode 100644 index 00000000000..11522ff1525 --- /dev/null +++ b/PWGLF/TableProducer/Strangeness/propagationService.cxx @@ -0,0 +1,146 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file propagationService.cxx +/// \brief +/// \author ALICE + +//=============================================================== +// +// Merged track propagation + strangeness building task +// +// Provides a common task to deal with track propagation and +// strangeness building in a single DPL device that is particularly +// adequate for pipelining. +// +// Currently meant for testing and performance evaluation +// +//=============================================================== + +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" +#include "Framework/RunningWorkflowInfo.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Common/Core/trackUtilities.h" +#include "ReconstructionDataFormats/DCA.h" +#include "DetectorsBase/Propagator.h" +#include "DetectorsBase/GeometryManager.h" +#include "CommonUtils/NameConf.h" +#include "CCDB/CcdbApi.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "CCDB/BasicCCDBManager.h" +#include "Framework/HistogramRegistry.h" +#include "DataFormatsCalibration/MeanVertexObject.h" +#include "CommonConstants/GeomConstants.h" +#include "PWGLF/Utils/strangenessBuilderModule.h" +#include "Common/Tools/TrackPropagationModule.h" +#include "Common/Tools/StandardCCDBLoader.h" + +using namespace o2; +using namespace o2::framework; +// using namespace o2::framework::expressions; + +// use parameters + cov mat non-propagated, aux info + (extension propagated) +using FullTracksExt = soa::Join; +using FullTracksExtIU = soa::Join; +using FullTracksExtWithPID = soa::Join; +using FullTracksExtIUWithPID = soa::Join; +using FullTracksExtLabeled = soa::Join; +using FullTracksExtLabeledIU = soa::Join; +using FullTracksExtLabeledWithPID = soa::Join; +using FullTracksExtLabeledIUWithPID = soa::Join; +using TracksWithExtra = soa::Join; + +// For dE/dx association in pre-selection +using TracksExtraWithPID = soa::Join; + +struct propagationService { + // CCDB boilerplate declarations + o2::framework::Configurable ccdburl{"ccdburl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Service ccdb; + + // propagation stuff + o2::common::StandardCCDBLoaderConfigurables standardCCDBLoaderConfigurables; + o2::common::StandardCCDBLoader ccdbLoader; + + // boilerplate: strangeness builder stuff + o2::pwglf::strangenessbuilder::products products; + o2::pwglf::strangenessbuilder::coreConfigurables baseOpts; + o2::pwglf::strangenessbuilder::v0Configurables v0BuilderOpts; + o2::pwglf::strangenessbuilder::cascadeConfigurables cascadeBuilderOpts; + o2::pwglf::strangenessbuilder::preSelectOpts preSelectOpts; + o2::pwglf::strangenessbuilder::BuilderModule strangenessBuilderModule; + + // track propagation + o2::common::TrackPropagationProducts trackPropagationProducts; + o2::common::TrackPropagationConfigurables trackPropagationConfigurables; + o2::common::TrackPropagationModule trackPropagation; + + // registry + HistogramRegistry histos{"histos"}; + + void init(o2::framework::InitContext& initContext) + { + // CCDB boilerplate init + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setURL(ccdburl.value); + + // task-specific + trackPropagation.init(trackPropagationConfigurables, histos, initContext); + strangenessBuilderModule.init(baseOpts, v0BuilderOpts, cascadeBuilderOpts, preSelectOpts, histos, initContext); + } + + void processRealData(soa::Join const& collisions, aod::V0s const& v0s, aod::Cascades const& cascades, aod::TrackedCascades const& trackedCascades, FullTracksExtIU const& tracks, aod::BCsWithTimestamps const& bcs) + { + ccdbLoader.initCCDBfromBCs(standardCCDBLoaderConfigurables, ccdb, bcs); + trackPropagation.fillTrackTables(trackPropagationConfigurables, ccdbLoader, collisions, tracks, trackPropagationProducts, histos); + strangenessBuilderModule.dataProcess(ccdb, histos, collisions, static_cast(nullptr), v0s, cascades, trackedCascades, tracks, bcs, static_cast(nullptr), products); + } + + void processMonteCarlo(soa::Join const& collisions, aod::McCollisions const& mccollisions, aod::V0s const& v0s, aod::Cascades const& cascades, aod::TrackedCascades const& trackedCascades, FullTracksExtLabeledIU const& tracks, aod::BCsWithTimestamps const& bcs, aod::McParticles const& mcParticles) + { + ccdbLoader.initCCDBfromBCs(standardCCDBLoaderConfigurables, ccdb, bcs); + trackPropagation.fillTrackTables(trackPropagationConfigurables, ccdbLoader, collisions, tracks, trackPropagationProducts, histos); + strangenessBuilderModule.dataProcess(ccdb, histos, collisions, mccollisions, v0s, cascades, trackedCascades, tracks, bcs, mcParticles, products); + } + + void processRealDataWithPID(soa::Join const& collisions, aod::V0s const& v0s, aod::Cascades const& cascades, aod::TrackedCascades const& trackedCascades, FullTracksExtIUWithPID const& tracks, aod::BCsWithTimestamps const& bcs) + { + ccdbLoader.initCCDBfromBCs(standardCCDBLoaderConfigurables, ccdb, bcs); + trackPropagation.fillTrackTables(trackPropagationConfigurables, ccdbLoader, collisions, tracks, trackPropagationProducts, histos); + strangenessBuilderModule.dataProcess(ccdb, histos, collisions, static_cast(nullptr), v0s, cascades, trackedCascades, tracks, bcs, static_cast(nullptr), products); + } + + void processMonteCarloWithPID(soa::Join const& collisions, aod::McCollisions const& mccollisions, aod::V0s const& v0s, aod::Cascades const& cascades, aod::TrackedCascades const& trackedCascades, FullTracksExtLabeledIUWithPID const& tracks, aod::BCsWithTimestamps const& bcs, aod::McParticles const& mcParticles) + { + ccdbLoader.initCCDBfromBCs(standardCCDBLoaderConfigurables, ccdb, bcs); + trackPropagation.fillTrackTables(trackPropagationConfigurables, ccdbLoader, collisions, tracks, trackPropagationProducts, histos); + strangenessBuilderModule.dataProcess(ccdb, histos, collisions, mccollisions, v0s, cascades, trackedCascades, tracks, bcs, mcParticles, products); + } + + PROCESS_SWITCH(propagationService, processRealData, "process real data", true); + PROCESS_SWITCH(propagationService, processMonteCarlo, "process monte carlo", false); + PROCESS_SWITCH(propagationService, processRealDataWithPID, "process real data", false); + PROCESS_SWITCH(propagationService, processMonteCarloWithPID, "process monte carlo", false); +}; + +//**************************************************************************************** +/** + * Workflow definition. + */ +//**************************************************************************************** +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{adaptAnalysisTask(cfgc)}; + return workflow; +} \ No newline at end of file diff --git a/PWGLF/Utils/strangenessBuilderModule.h b/PWGLF/Utils/strangenessBuilderModule.h new file mode 100644 index 00000000000..2c982acddee --- /dev/null +++ b/PWGLF/Utils/strangenessBuilderModule.h @@ -0,0 +1,2536 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file strangenessBuilderModule.h +/// \brief strangeness builder module +/// \author ALICE + +// simple checkers, but ensure 8 bit integers +#define BITSET(var, nbit) ((var) |= (static_cast(1) << static_cast(nbit))) + +#ifndef PWGLF_UTILS_STRANGENESSBUILDERMODULE_H_ +#define PWGLF_UTILS_STRANGENESSBUILDERMODULE_H_ + +#include +#include +#include +#include +#include +#include "Framework/AnalysisDataModel.h" +#include "Framework/Configurable.h" +#include "Framework/HistogramSpec.h" +#include "TableHelper.h" +#include "Common/Core/TPCVDriftManager.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/Utils/strangenessBuilderHelper.h" + +//__________________________________________ +// strangeness builder module + +namespace o2 +{ +namespace pwglf +{ +namespace strangenessbuilder // avoid polluting other namespaces +{ + +// statics necessary for the configurables in this namespace +static constexpr int nParameters = 1; +static const std::vector tableNames{ + "V0Indices", //.0 (standard analysis: V0Cores) + "V0CoresBase", //.1 (standard analyses: main table) + "V0Covs", //.2 (joinable with V0Cores) + "CascIndices", //.3 (standard analyses: CascData) + "KFCascIndices", //.4 (standard analyses: KFCascData) + "TraCascIndices", //.5 (standard analyses: TraCascData) + "StoredCascCores", //.6 (standard analyses: CascData, main table) + "StoredKFCascCores", //.7 (standard analyses: KFCascData, main table) + "StoredTraCascCores", //.8 (standard analyses: TraCascData, main table) + "CascCovs", //.9 (joinable with CascData) + "KFCascCovs", //.10 (joinable with KFCascData) + "TraCascCovs", //.11 (joinable with TraCascData) + "V0TrackXs", //.12 (joinable with V0Data) + "CascTrackXs", //.13 (joinable with CascData) + "CascBBs", //.14 (standard, bachelor-baryon vars) + "V0DauCovs", //.15 (requested: tracking studies) + "V0DauCovIUs", //.16 (requested: tracking studies) + "V0TraPosAtDCAs", //.17 (requested: tracking studies) + "V0TraPosAtIUs", //.18 (requested: tracking studies) + "V0Ivanovs", //.19 (requested: tracking studies) + "McV0Labels", //.20 (MC/standard analysis) + "V0MCCores", //.21 (MC, all generated desired V0s) + "V0CoreMCLabels", //.22 (MC, refs V0Cores to V0MCCores) + "V0MCCollRefs", //.23 (MC, refs V0MCCores to McCollisions) + "McCascLabels", //.24 (MC/standard analysis) + "McKFCascLabels", //.25 (MC, refs KFCascCores to CascMCCores) + "McTraCascLabels", //.26 (MC, refs TraCascCores to CascMCCores) + "McCascBBTags", //.27 (MC, joinable with CascCores, tags reco-ed) + "CascMCCores", //.28 (MC, all generated desired cascades) + "CascCoreMCLabels", //.29 (MC, refs CascCores to CascMCCores) + "CascMCCollRefs", // 30 (MC, refs CascMCCores to McCollisions) + "CascToTraRefs", //.31 (interlink CascCores -> TraCascCores) + "CascToKFRefs", //.32 (interlink CascCores -> KFCascCores) + "TraToCascRefs", //.33 (interlink TraCascCores -> CascCores) + "KFToCascRefs", //.34 (interlink KFCascCores -> CascCores) + "V0FoundTags", //.35 (tags found vs findable V0s in findable mode) + "CascFoundTags" //.36 (tags found vs findable Cascades in findable mode) +}; + +static constexpr int nTablesConst = 37; + +static const std::vector parameterNames{"enable"}; +static const int defaultParameters[nTablesConst][nParameters]{ + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, // index 9 + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, // index 19 + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, // index 29 + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}, + {-1}}; + +// table index : match order above +enum tableIndex { kV0Indices = 0, + kV0CoresBase, + kV0Covs, + kCascIndices, + kKFCascIndices, + kTraCascIndices, + kStoredCascCores, + kStoredKFCascCores, + kStoredTraCascCores, + kCascCovs, + kKFCascCovs, + kTraCascCovs, + kV0TrackXs, + kCascTrackXs, + kCascBBs, + kV0DauCovs, + kV0DauCovIUs, + kV0TraPosAtDCAs, + kV0TraPosAtIUs, + kV0Ivanovs, + kMcV0Labels, + kV0MCCores, + kV0CoreMCLabels, + kV0MCCollRefs, + kMcCascLabels, + kMcKFCascLabels, + kMcTraCascLabels, + kMcCascBBTags, + kCascMCCores, + kCascCoreMCLabels, + kCascMCCollRefs, + kCascToTraRefs, + kCascToKFRefs, + kTraToCascRefs, + kKFToCascRefs, + kV0FoundTags, + kCascFoundTags, + nTables }; + +enum V0PreSelection : uint8_t { selGamma = 0, + selK0Short, + selLambda, + selAntiLambda }; + +enum CascPreSelection : uint8_t { selXiMinus = 0, + selXiPlus, + selOmegaMinus, + selOmegaPlus }; + +static constexpr float defaultK0MassWindowParameters[1][4] = {{2.81882e-03, 1.14057e-03, 1.72138e-03, 5.00262e-01}}; +static constexpr float defaultLambdaWindowParameters[1][4] = {{1.17518e-03, 1.24099e-04, 5.47937e-03, 3.08009e-01}}; +static constexpr float defaultXiMassWindowParameters[1][4] = {{1.43210e-03, 2.03561e-04, 2.43187e-03, 7.99668e-01}}; +static constexpr float defaultOmMassWindowParameters[1][4] = {{1.43210e-03, 2.03561e-04, 2.43187e-03, 7.99668e-01}}; + +static constexpr float defaultLifetimeCuts[1][4] = {{20, 60, 40, 20}}; + +struct products : o2::framework::ProducesGroup { + //__________________________________________________ + // V0 tables + o2::framework::Produces v0indices; // standard part of V0Datas + o2::framework::Produces v0cores; // standard part of V0Datas + o2::framework::Produces v0covs; // for decay chain reco + + //__________________________________________________ + // cascade tables + o2::framework::Produces cascidx; // standard part of CascDatas + o2::framework::Produces kfcascidx; // standard part of KFCascDatas + o2::framework::Produces tracascidx; // standard part of TraCascDatas + o2::framework::Produces cascdata; // standard part of CascDatas + o2::framework::Produces kfcascdata; // standard part of KFCascDatas + o2::framework::Produces tracascdata; // standard part of TraCascDatas + o2::framework::Produces casccovs; // for decay chain reco + o2::framework::Produces kfcasccovs; // for decay chain reco + o2::framework::Produces tracasccovs; // for decay chain reco + + //__________________________________________________ + // interlink tables + o2::framework::Produces v0dataLink; // de-refs V0s -> V0Data + o2::framework::Produces cascdataLink; // de-refs Cascades -> CascData + o2::framework::Produces kfcascdataLink; // de-refs Cascades -> KFCascData + o2::framework::Produces tracascdataLink; // de-refs Cascades -> TraCascData + + //__________________________________________________ + // secondary auxiliary tables + o2::framework::Produces v0trackXs; // for decay chain reco + o2::framework::Produces cascTrackXs; // for decay chain reco + + //__________________________________________________ + // further auxiliary / optional if desired + o2::framework::Produces cascbb; + o2::framework::Produces v0daucovs; // covariances of daughter tracks + o2::framework::Produces v0daucovIUs; // covariances of daughter tracks + o2::framework::Produces v0dauPositions; // auxiliary debug information + o2::framework::Produces v0dauPositionsIU; // auxiliary debug information + o2::framework::Produces v0ivanovs; // information for Marian's tests + + //__________________________________________________ + // MC information: V0 + o2::framework::Produces v0labels; // MC labels for V0s + o2::framework::Produces v0mccores; // mc info storage + o2::framework::Produces v0CoreMCLabels; // interlink V0Cores -> V0MCCores + o2::framework::Produces v0mccollref; // references collisions from V0MCCores + + // MC information: Cascades + o2::framework::Produces casclabels; // MC labels for cascades + o2::framework::Produces kfcasclabels; // MC labels for KF cascades + o2::framework::Produces tracasclabels; // MC labels for tracked cascades + o2::framework::Produces bbtags; // bb tags (inv structure tagging in mc) + o2::framework::Produces cascmccores; // mc info storage + o2::framework::Produces cascCoreMClabels; // interlink CascCores -> CascMCCores + o2::framework::Produces cascmccollrefs; // references MC collisions from MC cascades + + //__________________________________________________ + // cascade interlinks + // FIXME: commented out until strangederivedbuilder adjusted accordingly + // o2::framework::Produces cascToTraRefs; // cascades -> tracked + // o2::framework::Produces cascToKFRefs; // cascades -> KF + // o2::framework::Produces traToCascRefs; // tracked -> cascades + // o2::framework::Produces kfToCascRefs; // KF -> cascades + + //__________________________________________________ + // Findable tags + o2::framework::Produces v0FoundTag; + o2::framework::Produces cascFoundTag; +}; + +// strangenessBuilder: 1st-order configurables +struct coreConfigurables : o2::framework::ConfigurableGroup { + o2::framework::Configurable> enabledTables{"enabledTables", + {defaultParameters[0], nTables, nParameters, tableNames, parameterNames}, + "Produce this table: -1 for autodetect; otherwise, 0/1 is false/true"}; + std::vector mEnabledTables; // Vector of enabled tables + + // first order deduplication implementation + // more algorithms to be added as necessary + o2::framework::Configurable deduplicationAlgorithm{"deduplicationAlgorithm", 1, "0: disabled; 1: best pointing angle wins; 2: best DCA daughters wins; 3: best pointing and best DCA wins"}; + + // V0 buffer for V0s used in cascades: master switch + // exchanges CPU (generate V0s again) with memory (save pre-generated V0s) + o2::framework::Configurable useV0BufferForCascades{"useV0BufferForCascades", false, "store array of V0s for cascades or not. False (default): save RAM, use more CPU; true: save CPU, use more RAM"}; + + o2::framework::Configurable mc_findableMode{"mc_findableMode", 0, "0: disabled; 1: add findable-but-not-found to existing V0s from AO2D; 2: reset V0s and generate only findable-but-not-found"}; +}; + +// strangenessBuilder: V0 building options +struct v0Configurables : o2::framework::ConfigurableGroup { + std::string prefix = "v0BuilderOpts"; + o2::framework::Configurable generatePhotonCandidates{"generatePhotonCandidates", false, "generate gamma conversion candidates (V0s using TPC-only tracks)"}; + o2::framework::Configurable moveTPCOnlyTracks{"moveTPCOnlyTracks", true, "if dealing with TPC-only tracks, move them according to TPC drift / time info"}; + + // baseline conditionals of V0 building + o2::framework::Configurable minCrossedRows{"minCrossedRows", 50, "minimum TPC crossed rows for daughter tracks"}; + o2::framework::Configurable dcanegtopv{"dcanegtopv", .1, "DCA Neg To PV"}; + o2::framework::Configurable dcapostopv{"dcapostopv", .1, "DCA Pos To PV"}; + o2::framework::Configurable v0cospa{"v0cospa", 0.95, "V0 CosPA"}; // double -> N.B. dcos(x)/dx = 0 at x=0) + o2::framework::Configurable dcav0dau{"dcav0dau", 1.0, "DCA V0 Daughters"}; + o2::framework::Configurable v0radius{"v0radius", 0.9, "v0radius"}; + o2::framework::Configurable maxDaughterEta{"maxDaughterEta", 5.0, "Maximum daughter eta (in abs value)"}; + + // MC builder options + o2::framework::Configurable mc_populateV0MCCoresSymmetric{"mc_populateV0MCCoresSymmetric", false, "populate V0MCCores table for derived data analysis, keep V0MCCores joinable with V0Cores"}; + o2::framework::Configurable mc_populateV0MCCoresAsymmetric{"mc_populateV0MCCoresAsymmetric", true, "populate V0MCCores table for derived data analysis, create V0Cores -> V0MCCores interlink. Saves only labeled V0s."}; + o2::framework::Configurable mc_treatPiToMuDecays{"mc_treatPiToMuDecays", true, "if true, will correctly capture pi -> mu and V0 label will still point to originating V0 decay in those cases. Nota bene: prong info will still be for the muon!"}; + o2::framework::Configurable mc_rapidityWindow{"mc_rapidityWindow", 0.5, "rapidity window to save non-recoed candidates"}; + o2::framework::Configurable mc_keepOnlyPhysicalPrimary{"mc_keepOnlyPhysicalPrimary", false, "Keep only physical primary generated V0s if not recoed"}; + o2::framework::Configurable mc_addGeneratedK0Short{"mc_addGeneratedK0Short", true, "add V0MCCore entry for generated, not-recoed K0Short"}; + o2::framework::Configurable mc_addGeneratedLambda{"mc_addGeneratedLambda", true, "add V0MCCore entry for generated, not-recoed Lambda"}; + o2::framework::Configurable mc_addGeneratedAntiLambda{"mc_addGeneratedAntiLambda", true, "add V0MCCore entry for generated, not-recoed AntiLambda"}; + o2::framework::Configurable mc_addGeneratedGamma{"mc_addGeneratedGamma", false, "add V0MCCore entry for generated, not-recoed Gamma"}; + o2::framework::Configurable mc_addGeneratedGammaMakeCollinear{"mc_addGeneratedGammaMakeCollinear", true, "when adding findable gammas, mark them as collinear"}; + o2::framework::Configurable mc_findableDetachedV0{"mc_findableDetachedV0", false, "if true, generate findable V0s that have collisionId -1. Caution advised."}; +}; + +// strangenessBuilder: cascade building options +struct cascadeConfigurables : o2::framework::ConfigurableGroup { + std::string prefix = "cascadeBuilderOpts"; + o2::framework::Configurable useCascadeMomentumAtPrimVtx{"useCascadeMomentumAtPrimVtx", false, "use cascade momentum at PV"}; + + // conditionals + o2::framework::Configurable minCrossedRows{"minCrossedRows", 50, "minimum TPC crossed rows for daughter tracks"}; + o2::framework::Configurable dcabachtopv{"dcabachtopv", .05, "DCA Bach To PV"}; + o2::framework::Configurable cascradius{"cascradius", 0.9, "cascradius"}; + o2::framework::Configurable casccospa{"casccospa", 0.95, "casccospa"}; + o2::framework::Configurable dcacascdau{"dcacascdau", 1.0, "DCA cascade Daughters"}; + o2::framework::Configurable lambdaMassWindow{"lambdaMassWindow", .010, "Distance from Lambda mass (does not apply to KF path)"}; + o2::framework::Configurable maxDaughterEta{"maxDaughterEta", 5.0, "Maximum daughter eta (in abs value)"}; + + // KF building specific + o2::framework::Configurable kfTuneForOmega{"kfTuneForOmega", false, "if enabled, take main cascade properties from Omega fit instead of Xi fit (= default)"}; + o2::framework::Configurable kfConstructMethod{"kfConstructMethod", 2, "KF Construct Method"}; + o2::framework::Configurable kfUseV0MassConstraint{"kfUseV0MassConstraint", true, "KF: use Lambda mass constraint"}; + o2::framework::Configurable kfUseCascadeMassConstraint{"kfUseCascadeMassConstraint", false, "KF: use Cascade mass constraint - WARNING: not adequate for inv mass analysis of Xi"}; + o2::framework::Configurable kfDoDCAFitterPreMinimV0{"kfDoDCAFitterPreMinimV0", true, "KF: do DCAFitter pre-optimization before KF fit to include material corrections for V0"}; + o2::framework::Configurable kfDoDCAFitterPreMinimCasc{"kfDoDCAFitterPreMinimCasc", true, "KF: do DCAFitter pre-optimization before KF fit to include material corrections for Xi"}; + + // MC builder options + o2::framework::Configurable mc_populateCascMCCoresSymmetric{"mc_populateCascMCCoresSymmetric", false, "populate CascMCCores table for derived data analysis, keep CascMCCores joinable with CascCores"}; + o2::framework::Configurable mc_populateCascMCCoresAsymmetric{"mc_populateCascMCCoresAsymmetric", true, "populate CascMCCores table for derived data analysis, create CascCores -> CascMCCores interlink. Saves only labeled Cascades."}; + o2::framework::Configurable mc_addGeneratedXiMinus{"mc_addGeneratedXiMinus", true, "add CascMCCore entry for generated, not-recoed XiMinus"}; + o2::framework::Configurable mc_addGeneratedXiPlus{"mc_addGeneratedXiPlus", true, "add CascMCCore entry for generated, not-recoed XiPlus"}; + o2::framework::Configurable mc_addGeneratedOmegaMinus{"mc_addGeneratedOmegaMinus", true, "add CascMCCore entry for generated, not-recoed OmegaMinus"}; + o2::framework::Configurable mc_addGeneratedOmegaPlus{"mc_addGeneratedOmegaPlus", true, "add CascMCCore entry for generated, not-recoed OmegaPlus"}; + o2::framework::Configurable mc_treatPiToMuDecays{"mc_treatPiToMuDecays", true, "if true, will correctly capture pi -> mu and V0 label will still point to originating V0 decay in those cases. Nota bene: prong info will still be for the muon!"}; + o2::framework::Configurable mc_rapidityWindow{"mc_rapidityWindow", 0.5, "rapidity window to save non-recoed candidates"}; + o2::framework::Configurable mc_keepOnlyPhysicalPrimary{"mc_keepOnlyPhysicalPrimary", false, "Keep only physical primary generated cascades if not recoed"}; + o2::framework::Configurable mc_findableDetachedCascade{"mc_findableDetachedCascade", false, "if true, generate findable cascades that have collisionId -1. Caution advised."}; +}; + +// preselection options +struct preSelectOpts : o2::framework::ConfigurableGroup { + std::string prefix = "preSelectOpts"; + o2::framework::Configurable preselectOnlyDesiredV0s{"preselectOnlyDesiredV0s", false, "preselect only V0s with compatible TPC PID and mass info"}; + o2::framework::Configurable preselectOnlyDesiredCascades{"preselectOnlyDesiredCascades", false, "preselect only Cascades with compatible TPC PID and mass info"}; + + // lifetime preselection options + // apply lifetime cuts to V0 and cascade candidates + // unit of measurement: centimeters + // lifetime of K0Short ~2.6844 cm, no feeddown and plenty to cut + // lifetime of Lambda ~7.9 cm but keep in mind feeddown from cascades + // lifetime of Xi ~4.91 cm + // lifetime of Omega ~2.461 cm + o2::framework::Configurable> lifetimeCut{"lifetimeCut", {defaultLifetimeCuts[0], 4, {"lifetimeCutK0S", "lifetimeCutLambda", "lifetimeCutXi", "lifetimeCutOmega"}}, "Lifetime cut for V0s and cascades (cm)"}; + + // mass preselection options + o2::framework::Configurable massCutPhoton{"massCutPhoton", 0.3, "Photon max mass"}; + o2::framework::Configurable> massCutK0{"massCutK0", {defaultK0MassWindowParameters[0], 4, {"constant", "linear", "expoConstant", "expoRelax"}}, "mass parameters for K0"}; + o2::framework::Configurable> massCutLambda{"massCutLambda", {defaultLambdaWindowParameters[0], 4, {"constant", "linear", "expoConstant", "expoRelax"}}, "mass parameters for Lambda"}; + o2::framework::Configurable> massCutXi{"massCutXi", {defaultXiMassWindowParameters[0], 4, {"constant", "linear", "expoConstant", "expoRelax"}}, "mass parameters for Xi"}; + o2::framework::Configurable> massCutOm{"massCutOm", {defaultOmMassWindowParameters[0], 4, {"constant", "linear", "expoConstant", "expoRelax"}}, "mass parameters for Omega"}; + o2::framework::Configurable massWindownumberOfSigmas{"massWindownumberOfSigmas", 20, "number of sigmas around mass peaks to keep"}; + o2::framework::Configurable massWindowSafetyMargin{"massWindowSafetyMargin", 0.001, "Extra mass window safety margin (in GeV/c2)"}; + + // TPC PID preselection options + o2::framework::Configurable maxTPCpidNsigma{"maxTPCpidNsigma", 5.0, "Maximum TPC PID N sigma (in abs value)"}; +}; + +class BuilderModule +{ + public: + BuilderModule() + { + // constructor + } + + // mass windows + float getMassSigmaK0Short(float pt) + { + return preSelectOpts.massCutK0->get("constant") + pt * preSelectOpts.massCutK0->get("linear") + preSelectOpts.massCutK0->get("expoConstant") * TMath::Exp(-pt / preSelectOpts.massCutK0->get("expoRelax")); + } + float getMassSigmaLambda(float pt) + { + return preSelectOpts.massCutLambda->get("constant") + pt * preSelectOpts.massCutLambda->get("linear") + preSelectOpts.massCutLambda->get("expoConstant") * TMath::Exp(-pt / preSelectOpts.massCutLambda->get("expoRelax")); + } + float getMassSigmaXi(float pt) + { + return preSelectOpts.massCutXi->get("constant") + pt * preSelectOpts.massCutXi->get("linear") + preSelectOpts.massCutXi->get("expoConstant") * TMath::Exp(-pt / preSelectOpts.massCutXi->get("expoRelax")); + } + float getMassSigmaOmega(float pt) + { + return preSelectOpts.massCutOm->get("constant") + pt * preSelectOpts.massCutOm->get("linear") + preSelectOpts.massCutOm->get("expoConstant") * TMath::Exp(-pt / preSelectOpts.massCutOm->get("expoRelax")); + } + + // helper object + o2::pwglf::strangenessBuilderHelper straHelper; + + // for handling TPC-only tracks (photons) + int mRunNumber; + o2::aod::common::TPCVDriftManager mVDriftMgr; + + // for establishing CascData/KFData/TraCascData interlinks + struct { + std::vector cascCoreToCascades; + std::vector kfCascCoreToCascades; + std::vector traCascCoreToCascades; + std::vector cascadeToCascCores; + std::vector cascadeToKFCascCores; + std::vector cascadeToTraCascCores; + } interlinks; + + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + // struct to add abstraction layer between V0s, Cascades and build indices + // desirable for adding extra (findable, etc) V0s, Cascades to built list + struct trackEntry { + int globalId = -1; + int originId = -1; + int mcCollisionId = -1; + int pdgCode = -1; + }; + struct v0Entry { + int globalId = -1; + int collisionId = -1; + int posTrackId = -1; + int negTrackId = -1; + int v0Type = 0; + int pdgCode = 0; // undefined if not MC - useful for faster finding + int particleId = -1; // de-reference the V0 particle if necessary + bool isCollinearV0 = false; + bool found = false; + }; + struct cascadeEntry { + int globalId = -1; + int collisionId = -1; + int v0Id = -1; + int posTrackId = -1; + int negTrackId = -1; + int bachTrackId = -1; + bool found = false; + }; + + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + // Helper struct to contain V0MCCore information prior to filling + struct mcV0info { + int label = -1; + int motherLabel = -1; + int pdgCode = 0; + int pdgCodeMother = 0; + int pdgCodePositive = 0; + int pdgCodeNegative = 0; + int mcCollision = -1; + bool isPhysicalPrimary = false; + int processPositive = -1; + int processNegative = -1; + std::array xyz; + std::array posP; + std::array negP; + std::array momentum; + }; + mcV0info thisInfo; + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + // Helper struct to contain CascMCCore information prior to filling + struct mcCascinfo { + int label; + int motherLabel; + int mcCollision; + int pdgCode; + int pdgCodeMother; + int pdgCodeV0; + int pdgCodePositive; + int pdgCodeNegative; + int pdgCodeBachelor; + bool isPhysicalPrimary; + int processPositive = -1; + int processNegative = -1; + int processBachelor = -1; + std::array xyz; + std::array lxyz; + std::array posP; + std::array negP; + std::array bachP; + std::array momentum; + int mcParticlePositive; + int mcParticleNegative; + int mcParticleBachelor; + }; + mcCascinfo thisCascInfo; + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + + std::vector v0List; + std::vector cascadeList; + std::vector sorted_v0; + std::vector sorted_cascade; + + // for tagging V0s used in cascades + std::vector v0sFromCascades; // Vector of v0 candidates used in cascades + std::vector ao2dV0toV0List; // index to relate v0s -> v0List + std::vector v0Map; // index to relate v0List -> v0sFromCascades + + // declaration of structs here + // (N.B.: will be invisible to the outside, create your own copies) + o2::pwglf::strangenessbuilder::coreConfigurables baseOpts; + o2::pwglf::strangenessbuilder::v0Configurables v0BuilderOpts; + o2::pwglf::strangenessbuilder::cascadeConfigurables cascadeBuilderOpts; + o2::pwglf::strangenessbuilder::preSelectOpts preSelectOpts; + + template + void init(TBaseConfigurables const& inputBaseOpts, TV0Configurables const& inputV0BuilderOpts, TCascadeConfigurables const& inputCascadeBuilderOpts, TPreSelOpts const& inputPreSelectOpts, THistoRegistry& histos, TInitContext& context) + { + // read in configurations from the task where it's used + // could be grouped even further, but should work + baseOpts = inputBaseOpts; + v0BuilderOpts = inputV0BuilderOpts; + cascadeBuilderOpts = inputCascadeBuilderOpts; + preSelectOpts = inputPreSelectOpts; + + // setup bookkeeping histogram + auto h = histos.template add("hTableBuildingStatistics", "hTableBuildingStatistics", o2::framework::kTH1D, {{nTablesConst, -0.5f, static_cast(nTablesConst)}}); + auto h2 = histos.template add("hInputStatistics", "hInputStatistics", o2::framework::kTH1D, {{nTablesConst, -0.5f, static_cast(nTablesConst)}}); + h2->SetTitle("Input table sizes"); + + if (v0BuilderOpts.generatePhotonCandidates.value == true) { + auto hDeduplicationStatistics = histos.template add("hDeduplicationStatistics", "hDeduplicationStatistics", o2::framework::kTH1D, {{2, -0.5f, 1.5f}}); + hDeduplicationStatistics->GetXaxis()->SetBinLabel(1, "AO2D V0s"); + hDeduplicationStatistics->GetXaxis()->SetBinLabel(2, "Deduplicated V0s"); + } + + if (preSelectOpts.preselectOnlyDesiredV0s.value == true) { + auto hPreselectionV0s = histos.template add("hPreselectionV0s", "hPreselectionV0s", o2::framework::kTH1D, {{16, -0.5f, 15.5f}}); + hPreselectionV0s->GetXaxis()->SetBinLabel(1, "Not preselected"); + hPreselectionV0s->GetXaxis()->SetBinLabel(2, "#gamma"); + hPreselectionV0s->GetXaxis()->SetBinLabel(3, "K^{0}_{S}"); + hPreselectionV0s->GetXaxis()->SetBinLabel(4, "#gamma, K^{0}_{S}"); + hPreselectionV0s->GetXaxis()->SetBinLabel(5, "#Lambda"); + hPreselectionV0s->GetXaxis()->SetBinLabel(6, "#gamma, #Lambda"); + hPreselectionV0s->GetXaxis()->SetBinLabel(7, "K^{0}_{S}, #Lambda"); + hPreselectionV0s->GetXaxis()->SetBinLabel(8, "#gamma, K^{0}_{S}, #Lambda"); + hPreselectionV0s->GetXaxis()->SetBinLabel(9, "#bar{#Lambda}"); + hPreselectionV0s->GetXaxis()->SetBinLabel(10, "#gamma, #bar{#Lambda}"); + hPreselectionV0s->GetXaxis()->SetBinLabel(11, "K^{0}_{S}, #bar{#Lambda}"); + hPreselectionV0s->GetXaxis()->SetBinLabel(12, "#gamma, K^{0}_{S}, #bar{#Lambda}"); + hPreselectionV0s->GetXaxis()->SetBinLabel(13, "#Lambda, #bar{#Lambda}"); + hPreselectionV0s->GetXaxis()->SetBinLabel(14, "#gamma, #Lambda, #bar{#Lambda}"); + hPreselectionV0s->GetXaxis()->SetBinLabel(15, "K^{0}_{S}, #Lambda, #bar{#Lambda}"); + hPreselectionV0s->GetXaxis()->SetBinLabel(16, "#gamma, K^{0}_{S}, #Lambda, #bar{#Lambda}"); + } + + if (preSelectOpts.preselectOnlyDesiredCascades.value == true) { + auto hPreselectionCascades = histos.template add("hPreselectionCascades", "hPreselectionCascades", o2::framework::kTH1D, {{16, -0.5f, 15.5f}}); + hPreselectionCascades->GetXaxis()->SetBinLabel(1, "Not preselected"); + hPreselectionCascades->GetXaxis()->SetBinLabel(2, "#Xi^{-}"); + hPreselectionCascades->GetXaxis()->SetBinLabel(3, "#Xi^{+}"); + hPreselectionCascades->GetXaxis()->SetBinLabel(4, "#Xi^{-}, #Xi^{+}"); + hPreselectionCascades->GetXaxis()->SetBinLabel(5, "#Omega^{-}"); + hPreselectionCascades->GetXaxis()->SetBinLabel(6, "#Xi^{-}, #Omega^{-}"); + hPreselectionCascades->GetXaxis()->SetBinLabel(7, "#Xi^{+}, #Omega^{-}"); + hPreselectionCascades->GetXaxis()->SetBinLabel(8, "#Xi^{-}, #Xi^{+}, #Omega^{-}"); + hPreselectionCascades->GetXaxis()->SetBinLabel(9, "#Omega^{+}"); + hPreselectionCascades->GetXaxis()->SetBinLabel(10, "#Xi^{-}, #Omega^{+}"); + hPreselectionCascades->GetXaxis()->SetBinLabel(11, "#Xi^{+}, #Omega^{+}"); + hPreselectionCascades->GetXaxis()->SetBinLabel(12, "#Xi^{-}, #Xi^{+}, #Omega^{+}"); + hPreselectionCascades->GetXaxis()->SetBinLabel(13, "#Omega^{-}, #Omega^{+}"); + hPreselectionCascades->GetXaxis()->SetBinLabel(14, "#Xi^{-}, #Omega^{-}, #Omega^{+}"); + hPreselectionCascades->GetXaxis()->SetBinLabel(15, "#Xi^{+}, #Omega^{-}, #Omega^{+}"); + hPreselectionCascades->GetXaxis()->SetBinLabel(16, "#Xi^{-}, #Xi^{+}, #Omega^{-}, #Omega^{+}"); + } + + if (baseOpts.mc_findableMode.value > 0) { + // save statistics of findable candidate processing + auto hFindable = histos.template add("hFindableStatistics", "hFindableStatistics", o2::framework::kTH1D, {{6, -0.5f, 5.5f}}); + hFindable->SetTitle(Form("Findable mode: %i", static_cast(baseOpts.mc_findableMode.value))); + hFindable->GetXaxis()->SetBinLabel(1, "AO2D V0s"); + hFindable->GetXaxis()->SetBinLabel(2, "V0s to be built"); + hFindable->GetXaxis()->SetBinLabel(3, "V0s with collId -1"); + hFindable->GetXaxis()->SetBinLabel(4, "AO2D Cascades"); + hFindable->GetXaxis()->SetBinLabel(5, "Cascades to be built"); + hFindable->GetXaxis()->SetBinLabel(6, "Cascades with collId -1"); + } + + auto hPrimaryV0s = histos.template add("hPrimaryV0s", "hPrimaryV0s", o2::framework::kTH1D, {{2, -0.5f, 1.5f}}); + hPrimaryV0s->GetXaxis()->SetBinLabel(1, "All V0s"); + hPrimaryV0s->GetXaxis()->SetBinLabel(2, "Primary V0s"); + + mRunNumber = 0; + + baseOpts.mEnabledTables.resize(nTables, 0); + + LOGF(info, "Configuring tables to generate"); + auto& workflows = context.services().template get(); + + TString listOfRequestors[nTables]; + for (int i = 0; i < nTables; i++) { + // adjust bookkeeping histogram + h->GetXaxis()->SetBinLabel(i + 1, tableNames[i].c_str()); + h2->GetXaxis()->SetBinLabel(i + 1, tableNames[i].c_str()); + h->SetBinContent(i + 1, -1); // mark all as disabled to start + + int f = baseOpts.enabledTables->get(tableNames[i].c_str(), "enable"); + if (f == 1) { + baseOpts.mEnabledTables[i] = 1; + listOfRequestors[i] = "manual enabling"; + } + if (f == -1) { + // autodetect this table in other devices + for (o2::framework::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("strangenessbuilder-initializer") == 0) + continue; // don't listen to the initializer + if (o2::framework::DataSpecUtils::partialMatch(input.matcher, o2::header::DataOrigin("AOD"))) { + auto&& [origin, description, version] = o2::framework::DataSpecUtils::asConcreteDataMatcher(input.matcher); + std::string tableNameWithVersion = tableNames[i]; + if (version > 0) { + tableNameWithVersion += Form("_%03d", version); + } + if (input.matcher.binding == tableNameWithVersion) { + LOGF(info, "Device %s has subscribed to %s (version %i)", device.name, tableNames[i], version); + listOfRequestors[i].Append(Form("%s ", device.name.c_str())); + baseOpts.mEnabledTables[i] = 1; + } + } + } + } + } + } + + LOGF(info, "*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*"); + LOGF(info, " Strangeness builder: basic configuration listing"); + LOGF(info, "*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*"); + + if (baseOpts.mc_findableMode.value == 1) { + LOGF(info, " ===> findable mode 1 is enabled: complement reco-ed with non-found findable"); + } + if (baseOpts.mc_findableMode.value == 2) { + LOGF(info, " ===> findable mode 2 is enabled: re-generate all findable from scratch"); + } + + // list enabled tables + + for (int i = 0; i < nTables; i++) { + // printout to be improved in the future + if (baseOpts.mEnabledTables[i]) { + LOGF(info, " -~> Table enabled: %s, requested by %s", tableNames[i], listOfRequestors[i].Data()); + h->SetBinContent(i + 1, 0); // mark enabled + } + } + LOGF(info, "*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*"); + // print base cuts + LOGF(info, "-~> V0 | min crossed rows ..............: %i", v0BuilderOpts.minCrossedRows.value); + LOGF(info, "-~> V0 | DCA pos track to PV ...........: %f", v0BuilderOpts.dcapostopv.value); + LOGF(info, "-~> V0 | DCA neg track to PV ...........: %f", v0BuilderOpts.dcanegtopv.value); + LOGF(info, "-~> V0 | V0 cosine of PA ...............: %f", v0BuilderOpts.v0cospa.value); + LOGF(info, "-~> V0 | DCA between V0 daughters ......: %f", v0BuilderOpts.dcav0dau.value); + LOGF(info, "-~> V0 | V0 2D decay radius ............: %f", v0BuilderOpts.v0radius.value); + LOGF(info, "-~> V0 | Maximum daughter eta ..........: %f", v0BuilderOpts.maxDaughterEta.value); + + LOGF(info, "-~> Cascade | min crossed rows .........: %i", cascadeBuilderOpts.minCrossedRows.value); + LOGF(info, "-~> Cascade | DCA bach track to PV .....: %f", cascadeBuilderOpts.dcabachtopv.value); + LOGF(info, "-~> Cascade | Cascade cosine of PA .....: %f", cascadeBuilderOpts.casccospa.value); + LOGF(info, "-~> Cascade | Cascade daughter DCA .....: %f", cascadeBuilderOpts.dcacascdau.value); + LOGF(info, "-~> Cascade | Cascade radius ...........: %f", cascadeBuilderOpts.cascradius.value); + LOGF(info, "-~> Cascade | Lambda mass window .......: %f", cascadeBuilderOpts.lambdaMassWindow.value); + LOGF(info, "-~> Cascade | Maximum daughter eta .....: %f", cascadeBuilderOpts.maxDaughterEta.value); + LOGF(info, "*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*"); + + // set V0 parameters in the helper + straHelper.v0selections.minCrossedRows = v0BuilderOpts.minCrossedRows; + straHelper.v0selections.dcanegtopv = v0BuilderOpts.dcanegtopv; + straHelper.v0selections.dcapostopv = v0BuilderOpts.dcapostopv; + straHelper.v0selections.v0cospa = v0BuilderOpts.v0cospa; + straHelper.v0selections.dcav0dau = v0BuilderOpts.dcav0dau; + straHelper.v0selections.v0radius = v0BuilderOpts.v0radius; + straHelper.v0selections.maxDaughterEta = v0BuilderOpts.maxDaughterEta; + + // set cascade parameters in the helper + straHelper.cascadeselections.minCrossedRows = cascadeBuilderOpts.minCrossedRows; + straHelper.cascadeselections.dcabachtopv = cascadeBuilderOpts.dcabachtopv; + straHelper.cascadeselections.cascradius = cascadeBuilderOpts.cascradius; + straHelper.cascadeselections.casccospa = cascadeBuilderOpts.casccospa; + straHelper.cascadeselections.dcacascdau = cascadeBuilderOpts.dcacascdau; + straHelper.cascadeselections.lambdaMassWindow = cascadeBuilderOpts.lambdaMassWindow; + straHelper.cascadeselections.maxDaughterEta = cascadeBuilderOpts.maxDaughterEta; + } + + // for sorting + template + std::vector sort_indices(const std::vector& v, bool doSorting = false) + { + std::vector idx(v.size()); + std::iota(idx.begin(), idx.end(), 0); + if (doSorting) { + // do sorting only if requested (not always necessary) + std::stable_sort(idx.begin(), idx.end(), + [&v](std::size_t i1, std::size_t i2) { return v[i1].collisionId < v[i2].collisionId; }); + } + return idx; + } + + template + bool initCCDB(TCCDB& ccdb, aod::BCsWithTimestamps const& bcs, TCollisions const& collisions) + { + auto bc = collisions.size() ? collisions.begin().template bc_as() : bcs.begin(); + if (!bcs.size()) { + LOGF(warn, "No BC found, skipping this DF."); + return false; // signal to skip this DF + } + + if (mRunNumber == bc.runNumber()) { + return true; + } + + auto timestamp = bc.timestamp(); + + // Fetch magnetic field from ccdb for current collision + auto magneticField = o2::base::Propagator::Instance()->getNominalBz(); + LOG(info) << "Configuring for timestamp " << timestamp << " with magnetic field of " << magneticField << " kG"; + + // Set magnetic field value once known + straHelper.fitter.setBz(magneticField); + + LOG(info) << "Fully configured for run: " << bc.runNumber(); + // mmark this run as configured + mRunNumber = bc.runNumber(); + + if (v0BuilderOpts.generatePhotonCandidates.value && v0BuilderOpts.moveTPCOnlyTracks.value) { + // initialize only if needed, avoid unnecessary CCDB calls + mVDriftMgr.init(&ccdb->instance()); + mVDriftMgr.update(timestamp); + } + + return true; + } + + //__________________________________________________ + void resetInterlinks() + { + interlinks.cascCoreToCascades.clear(); + interlinks.kfCascCoreToCascades.clear(); + interlinks.traCascCoreToCascades.clear(); + interlinks.cascadeToCascCores.clear(); + interlinks.cascadeToKFCascCores.clear(); + interlinks.cascadeToTraCascCores.clear(); + } + + //__________________________________________________ + void populateCascadeInterlinks() + { + // if (mEnabledTables[kCascToKFRefs]) { + // for (const auto& cascCore : interlinks.cascCoreToCascades) { + // cascToKFRefs(interlinks.cascadeToKFCascCores[cascCore]); + // histos.fill(HIST("hTableBuildingStatistics"), kCascToKFRefs); + // } + // } + // if (mEnabledTables[kCascToTraRefs]) { + // for (const auto& cascCore : interlinks.cascCoreToCascades) { + // cascToTraRefs(interlinks.cascadeToTraCascCores[cascCore]); + // histos.fill(HIST("hTableBuildingStatistics"), kCascToTraRefs); + // } + // } + // if (mEnabledTables[kKFToCascRefs]) { + // for (const auto& kfCascCore : interlinks.kfCascCoreToCascades) { + // kfToCascRefs(interlinks.cascadeToCascCores[kfCascCore]); + // histos.fill(HIST("hTableBuildingStatistics"), kKFToCascRefs); + // } + // } + // if (mEnabledTables[kTraToCascRefs]) { + // for (const auto& traCascCore : interlinks.traCascCoreToCascades) { + // traToCascRefs(interlinks.cascadeToCascCores[traCascCore]); + // histos.fill(HIST("hTableBuildingStatistics"), kTraToCascRefs); + // } + // } + } + + //__________________________________________________ + template + void prepareBuildingLists(THistoRegistry& histos, TCollisions const& collisions, TMCCollisions const& mcCollisions, TV0s const& v0s, TCascades const& cascades, TTracks const& tracks, TMCParticles const& mcParticles) + { + // this function prepares the v0List and cascadeList depending on + // how the task has been set up. Standard operation simply uses + // the existing V0s and Cascades from AO2D, while findable MC + // operation either complements with all findable-but-not-found + // or resets and fills with all findable. + // + // Whenever using findable candidates, they will be appropriately + // marked for posterior analysis using 'tag' variables. + // + // findable mode legend: + // 0: simple passthrough of V0s, Cascades in AO2Ds + // (in data, this is the only mode possible!) + // 1: add extra findable that haven't been found + // 2: generate only findable (no background) + + // redo lists from scratch + v0List.clear(); + cascadeList.clear(); + sorted_v0.clear(); + sorted_cascade.clear(); + ao2dV0toV0List.clear(); + + trackEntry currentTrackEntry; + v0Entry currentV0Entry; + cascadeEntry currentCascadeEntry; + + std::vector bestCollisionArray; // stores McCollision -> Collision map + std::vector bestCollisionNContribsArray; // stores Ncontribs for biggest coll assoc to mccoll + + int collisionLessV0s = 0; + int collisionLessCascades = 0; + + if (baseOpts.mc_findableMode.value > 0) { + if constexpr (soa::is_table) { + // if mcCollisions exist, assemble mcColl -> bestRecoColl map here + bestCollisionArray.clear(); + bestCollisionNContribsArray.clear(); + bestCollisionArray.resize(mcCollisions.size(), -1); // marks not reconstructed + bestCollisionNContribsArray.resize(mcCollisions.size(), -1); // marks not reconstructed + + // single loop over double loop at a small cost in memory for extra array + for (const auto& collision : collisions) { + if (collision.has_mcCollision()) { + if (collision.numContrib() > bestCollisionNContribsArray[collision.mcCollisionId()]) { + bestCollisionArray[collision.mcCollisionId()] = collision.globalIndex(); + bestCollisionNContribsArray[collision.mcCollisionId()] = collision.numContrib(); + } + } + } // end collision loop + } // end is_table + } // end findable mode check + + if (baseOpts.mc_findableMode.value < 2) { + // keep all unless de-duplication active + ao2dV0toV0List.resize(v0s.size(), -1); // -1 means keep, -2 means do not keep + + if (baseOpts.deduplicationAlgorithm > 0 && v0BuilderOpts.generatePhotonCandidates) { + // handle duplicates explicitly: group V0s according to (p,n) indices + // will provide a list of collisionIds (in V0group), allowing for + // easy de-duplication when passing to the v0List + std::vector v0tableGrouped = o2::pwglf::groupDuplicates(v0s); + histos.fill(HIST("hDeduplicationStatistics"), 0.0, v0s.size()); + histos.fill(HIST("hDeduplicationStatistics"), 1.0, v0tableGrouped.size()); + + // process grouped duplicates, remove 'bad' ones + for (size_t iV0 = 0; iV0 < v0tableGrouped.size(); iV0++) { + auto pTrack = tracks.rawIteratorAt(v0tableGrouped[iV0].posTrackId); + auto nTrack = tracks.rawIteratorAt(v0tableGrouped[iV0].negTrackId); + + bool isPosTPCOnly = (pTrack.hasTPC() && !pTrack.hasITS() && !pTrack.hasTRD() && !pTrack.hasTOF()); + bool isNegTPCOnly = (nTrack.hasTPC() && !nTrack.hasITS() && !nTrack.hasTRD() && !nTrack.hasTOF()); + + // skip single copy V0s + if (v0tableGrouped[iV0].collisionIds.size() == 1) { + continue; + } + + // don't try to de-duplicate if no track is TPC only + if (!isPosTPCOnly && !isNegTPCOnly) { + continue; + } + + // fitness criteria defined here + float bestPointingAngle = 10; // a nonsense angle, anything's better + size_t bestPointingAngleIndex = -1; + + float bestDCADaughters = 1e+3; // an excessively large DCA + size_t bestDCADaughtersIndex = -1; + + for (size_t ic = 0; ic < v0tableGrouped[iV0].collisionIds.size(); ic++) { + // get track parametrizations, collisions + auto posTrackPar = getTrackParCov(pTrack); + auto negTrackPar = getTrackParCov(nTrack); + auto const& collision = collisions.rawIteratorAt(v0tableGrouped[iV0].collisionIds[ic]); + + // handle TPC-only tracks properly (photon conversions) + if (v0BuilderOpts.moveTPCOnlyTracks) { + if (isPosTPCOnly) { + // Nota bene: positive is TPC-only -> this entire V0 merits treatment as photon candidate + posTrackPar.setPID(o2::track::PID::Electron); + negTrackPar.setPID(o2::track::PID::Electron); + if (!mVDriftMgr.moveTPCTrack(collision, pTrack, posTrackPar)) { + return; + } + } + if (isNegTPCOnly) { + // Nota bene: negative is TPC-only -> this entire V0 merits treatment as photon candidate + posTrackPar.setPID(o2::track::PID::Electron); + negTrackPar.setPID(o2::track::PID::Electron); + if (!mVDriftMgr.moveTPCTrack(collision, nTrack, negTrackPar)) { + return; + } + } + } // end TPC drift treatment + + // process candidate with helper, generate properties for consulting + // : do not apply selections: do as much as possible to preserve + // candidate at this level and do not select with topo selections + if (straHelper.buildV0Candidate(v0tableGrouped[iV0].collisionIds[ic], collision.posX(), collision.posY(), collision.posZ(), pTrack, nTrack, posTrackPar, negTrackPar, true, false, true)) { + // candidate built, check pointing angle + if (straHelper.v0.pointingAngle < bestPointingAngle) { + bestPointingAngle = straHelper.v0.pointingAngle; + bestPointingAngleIndex = ic; + } + if (straHelper.v0.daughterDCA < bestDCADaughters) { + bestDCADaughters = straHelper.v0.daughterDCA; + bestDCADaughtersIndex = ic; + } + } // end build V0 + } // end candidate loop + + // mark de-duplicated candidates + for (size_t ic = 0; ic < v0tableGrouped[iV0].collisionIds.size(); ic++) { + ao2dV0toV0List[v0tableGrouped[iV0].V0Ids[ic]] = -2; + // algorithm 1: best pointing angle + if (bestPointingAngleIndex == ic && baseOpts.deduplicationAlgorithm.value == 1) { + ao2dV0toV0List[v0tableGrouped[iV0].V0Ids[ic]] = -1; // keep best only + } + if (bestDCADaughtersIndex == ic && baseOpts.deduplicationAlgorithm.value == 2) { + ao2dV0toV0List[v0tableGrouped[iV0].V0Ids[ic]] = -1; // keep best only + } + if (bestDCADaughtersIndex == ic && bestPointingAngleIndex == ic && baseOpts.deduplicationAlgorithm.value == 3) { + ao2dV0toV0List[v0tableGrouped[iV0].V0Ids[ic]] = -1; // keep best only + } + } + } // end V0 loop + } // end de-duplication process + + for (const auto& v0 : v0s) { + if (ao2dV0toV0List[v0.globalIndex()] == -1) { // keep only de-duplicated + ao2dV0toV0List[v0.globalIndex()] = v0List.size(); // maps V0s to the corresponding v0List entry + currentV0Entry.globalId = v0.globalIndex(); + currentV0Entry.collisionId = v0.collisionId(); + currentV0Entry.posTrackId = v0.posTrackId(); + currentV0Entry.negTrackId = v0.negTrackId(); + currentV0Entry.v0Type = v0.v0Type(); + currentV0Entry.pdgCode = 0; + currentV0Entry.particleId = -1; + currentV0Entry.isCollinearV0 = v0.isCollinearV0(); + currentV0Entry.found = true; + v0List.push_back(currentV0Entry); + } + } + } + // any mode other than 0 will require mcParticles + if constexpr (soa::is_table) { + if (baseOpts.mc_findableMode.value > 0) { + // for search if existing or not + int v0ListReconstructedSize = v0List.size(); + + // find extra candidates, step 1: find subset of tracks that interest + std::vector positiveTrackArray; + std::vector negativeTrackArray; + // vector elements: track index, origin index [, mc collision id, pdg code] + int dummy = -1; // unnecessary in this path + for (const auto& track : tracks) { + if (!track.has_mcParticle()) { + continue; // skip this, it's trouble + } + auto particle = track.template mcParticle_as(); + int originParticleIndex = getOriginatingParticle(particle, dummy, v0BuilderOpts.mc_treatPiToMuDecays); + if (originParticleIndex < 0) { + continue; // skip this, it's trouble (2) + } + auto originParticle = mcParticles.rawIteratorAt(originParticleIndex); + + bool trackIsInteresting = false; + if ( + (originParticle.pdgCode() == 310 && v0BuilderOpts.mc_addGeneratedK0Short.value > 0) || + (originParticle.pdgCode() == 3122 && v0BuilderOpts.mc_addGeneratedLambda.value > 0) || + (originParticle.pdgCode() == -3122 && v0BuilderOpts.mc_addGeneratedAntiLambda.value > 0) || + (originParticle.pdgCode() == 22 && v0BuilderOpts.mc_addGeneratedGamma.value > 0)) { + trackIsInteresting = true; + } + if (!trackIsInteresting) { + continue; // skip this, it's uninteresting + } + + currentTrackEntry.globalId = static_cast(track.globalIndex()); + currentTrackEntry.originId = originParticleIndex; + currentTrackEntry.mcCollisionId = originParticle.mcCollisionId(); + currentTrackEntry.pdgCode = originParticle.pdgCode(); + + // now separate according to particle species + if (track.sign() < 0) { + negativeTrackArray.push_back(currentTrackEntry); + } else { + positiveTrackArray.push_back(currentTrackEntry); + } + } + + // Nested loop only with valuable tracks + for (const auto& positiveTrackIndex : positiveTrackArray) { + for (const auto& negativeTrackIndex : negativeTrackArray) { + if (positiveTrackIndex.originId != negativeTrackIndex.originId) { + continue; // not the same originating particle + } + // findable mode 1: add non-reconstructed as v0Type 8 + if (baseOpts.mc_findableMode.value == 1) { + bool detected = false; + for (int ii = 0; ii < v0ListReconstructedSize; ii++) { + // check if this particular combination already exists in v0List + if (v0List[ii].posTrackId == positiveTrackIndex.globalId && + v0List[ii].negTrackId == negativeTrackIndex.globalId) { + detected = true; + // override pdg code with something useful for cascade findable math + v0List[ii].pdgCode = positiveTrackIndex.pdgCode; + break; + } + } + if (detected == false) { + // collision index: from best-version-of-this-mcCollision + // nota bene: this could be negative, caution advised + currentV0Entry.globalId = -1; + currentV0Entry.collisionId = bestCollisionArray[positiveTrackIndex.mcCollisionId]; + currentV0Entry.posTrackId = positiveTrackIndex.globalId; + currentV0Entry.negTrackId = negativeTrackIndex.globalId; + currentV0Entry.v0Type = 1; + currentV0Entry.pdgCode = positiveTrackIndex.pdgCode; + currentV0Entry.particleId = positiveTrackIndex.originId; + currentV0Entry.isCollinearV0 = false; + if (v0BuilderOpts.mc_addGeneratedGammaMakeCollinear.value && currentV0Entry.pdgCode == 22) { + currentV0Entry.isCollinearV0 = true; + } + currentV0Entry.found = false; + if (bestCollisionArray[positiveTrackIndex.mcCollisionId] < 0) { + collisionLessV0s++; + } + if (v0BuilderOpts.mc_findableDetachedV0.value || currentV0Entry.collisionId >= 0) { + v0List.push_back(currentV0Entry); + } + } + } + // findable mode 2 + if (baseOpts.mc_findableMode.value == 2) { + currentV0Entry.globalId = -1; + currentV0Entry.collisionId = bestCollisionArray[positiveTrackIndex.mcCollisionId]; + currentV0Entry.posTrackId = positiveTrackIndex.globalId; + currentV0Entry.negTrackId = negativeTrackIndex.globalId; + currentV0Entry.v0Type = 1; + currentV0Entry.pdgCode = positiveTrackIndex.pdgCode; + currentV0Entry.particleId = positiveTrackIndex.originId; + currentV0Entry.isCollinearV0 = false; + if (v0BuilderOpts.mc_addGeneratedGammaMakeCollinear.value && currentV0Entry.pdgCode == 22) { + currentV0Entry.isCollinearV0 = true; + } + currentV0Entry.found = false; + for (const auto& v0 : v0s) { + if (v0.posTrackId() == positiveTrackIndex.globalId && + v0.negTrackId() == negativeTrackIndex.globalId) { + // this will override type, but not collision index + // N.B.: collision index checks still desirable! + currentV0Entry.globalId = v0.globalIndex(); + currentV0Entry.v0Type = v0.v0Type(); + currentV0Entry.isCollinearV0 = v0.isCollinearV0(); + currentV0Entry.found = true; + break; + } + } + if (v0BuilderOpts.mc_findableDetachedV0.value || currentV0Entry.collisionId >= 0) { + v0List.push_back(currentV0Entry); + } + } + } + } // end positive / negative track loops + + // fill findable statistics table + histos.fill(HIST("hFindableStatistics"), 0.0, v0s.size()); + histos.fill(HIST("hFindableStatistics"), 1.0, v0List.size()); + histos.fill(HIST("hFindableStatistics"), 2.0, collisionLessV0s); + + } // end findableMode > 0 check + } // end soa::is_table + + // determine properly collision-id-sorted index array for later use + // N.B.: necessary also before cascade part + sorted_v0.clear(); + sorted_v0 = sort_indices(v0List, (baseOpts.mc_findableMode.value > 0)); + + // Cascade part if cores are requested, skip otherwise + if (baseOpts.mEnabledTables[kStoredCascCores] || baseOpts.mEnabledTables[kStoredKFCascCores]) { + if (baseOpts.mc_findableMode.value < 2) { + // simple passthrough: copy existing cascades to build list + for (const auto& cascade : cascades) { + auto const& v0 = cascade.v0(); + currentCascadeEntry.globalId = cascade.globalIndex(); + currentCascadeEntry.collisionId = cascade.collisionId(); + currentCascadeEntry.v0Id = ao2dV0toV0List[v0.globalIndex()]; + currentCascadeEntry.posTrackId = v0.posTrackId(); + currentCascadeEntry.negTrackId = v0.negTrackId(); + currentCascadeEntry.bachTrackId = cascade.bachelorId(); + currentCascadeEntry.found = true; + cascadeList.push_back(currentCascadeEntry); + } + } + + // any mode other than 0 will require mcParticles + if constexpr (soa::is_table) { + if (baseOpts.mc_findableMode.value > 0) { + // for search if existing or not + size_t cascadeListReconstructedSize = cascadeList.size(); + + // determine which tracks are of interest + std::vector bachelorTrackArray; + // vector elements: track index, origin index, mc collision id, pdg code] + int dummy = -1; // unnecessary in this path + for (const auto& track : tracks) { + if (!track.has_mcParticle()) { + continue; // skip this, it's trouble + } + auto particle = track.template mcParticle_as(); + int originParticleIndex = getOriginatingParticle(particle, dummy, cascadeBuilderOpts.mc_treatPiToMuDecays); + if (originParticleIndex < 0) { + continue; // skip this, it's trouble (2) + } + auto originParticle = mcParticles.rawIteratorAt(originParticleIndex); + + bool trackIsInteresting = false; + if ( + (originParticle.pdgCode() == 3312 && cascadeBuilderOpts.mc_addGeneratedXiMinus.value > 0) || + (originParticle.pdgCode() == -3312 && cascadeBuilderOpts.mc_addGeneratedXiPlus.value > 0) || + (originParticle.pdgCode() == 3334 && cascadeBuilderOpts.mc_addGeneratedOmegaMinus.value > 0) || + (originParticle.pdgCode() == -3334 && cascadeBuilderOpts.mc_addGeneratedOmegaPlus.value > 0)) { + trackIsInteresting = true; + } + if (!trackIsInteresting) { + continue; // skip this, it's uninteresting + } + + currentTrackEntry.globalId = static_cast(track.globalIndex()); + currentTrackEntry.originId = originParticleIndex; + currentTrackEntry.mcCollisionId = originParticle.mcCollisionId(); + currentTrackEntry.pdgCode = originParticle.pdgCode(); + + // populate list of bachelor tracks to pair + bachelorTrackArray.push_back(currentTrackEntry); + } + + // determine which V0s are of interest to pair and do pairing + for (size_t v0i = 0; v0i < v0List.size(); v0i++) { + auto v0 = v0List[sorted_v0[v0i]]; + + if (std::abs(v0.pdgCode) != 3122) { + continue; // this V0 isn't a lambda, can't come from a cascade: skip + } + if (v0.particleId < 0) { + continue; // no de-referencing possible (e.g. background, ...) + } + auto v0Particle = mcParticles.rawIteratorAt(v0.particleId); + + int v0OriginParticleIndex = -1; + if (v0Particle.has_mothers()) { + auto const& motherList = v0Particle.template mothers_as(); + if (motherList.size() == 1) { + for (const auto& mother : motherList) { + v0OriginParticleIndex = mother.globalIndex(); + } + } + } + if (v0OriginParticleIndex < 0) { + continue; + } + auto v0OriginParticle = mcParticles.rawIteratorAt(v0OriginParticleIndex); + + if (std::abs(v0OriginParticle.pdgCode()) != 3312 && std::abs(v0OriginParticle.pdgCode()) != 3334) { + continue; // this V0 does not come from any particle of interest, don't try + } + for (const auto& bachelorTrackIndex : bachelorTrackArray) { + if (bachelorTrackIndex.originId != v0OriginParticle.globalIndex()) { + continue; + } + // if we are here: v0 origin is 3312 or 3334, bachelor origin matches V0 origin + // findable mode 1: add non-reconstructed as cascadeType 1 + if (baseOpts.mc_findableMode.value == 1) { + bool detected = false; + for (size_t ii = 0; ii < cascadeListReconstructedSize; ii++) { + // check if this particular combination already exists in cascadeList + // caution: use track indices (immutable) but not V0 indices (re-indexing) + if (cascadeList[ii].posTrackId == v0.posTrackId && + cascadeList[ii].negTrackId == v0.negTrackId && + cascadeList[ii].bachTrackId == bachelorTrackIndex.globalId) { + detected = true; + break; + } + } + if (detected == false) { + // collision index: from best-version-of-this-mcCollision + // nota bene: this could be negative, caution advised + currentCascadeEntry.globalId = -1; + currentCascadeEntry.collisionId = bestCollisionArray[bachelorTrackIndex.mcCollisionId]; + currentCascadeEntry.v0Id = v0i; // correct information here + currentCascadeEntry.posTrackId = v0.posTrackId; + currentCascadeEntry.negTrackId = v0.negTrackId; + currentCascadeEntry.bachTrackId = bachelorTrackIndex.globalId; + currentCascadeEntry.found = false; + cascadeList.push_back(currentCascadeEntry); + if (bestCollisionArray[bachelorTrackIndex.mcCollisionId] < 0) { + collisionLessCascades++; + } + if (cascadeBuilderOpts.mc_findableDetachedCascade.value || currentCascadeEntry.collisionId >= 0) { + cascadeList.push_back(currentCascadeEntry); + } + } + } + + // findable mode 2: determine type based on cascade table, + // with type 1 being reserved to findable-but-not-found + if (baseOpts.mc_findableMode.value == 2) { + currentCascadeEntry.globalId = -1; + currentCascadeEntry.collisionId = bestCollisionArray[bachelorTrackIndex.mcCollisionId]; + currentCascadeEntry.v0Id = v0i; // fill this in one go later + currentCascadeEntry.posTrackId = v0.posTrackId; + currentCascadeEntry.negTrackId = v0.negTrackId; + currentCascadeEntry.bachTrackId = bachelorTrackIndex.globalId; + currentCascadeEntry.found = false; + if (bestCollisionArray[bachelorTrackIndex.mcCollisionId] < 0) { + collisionLessCascades++; + } + for (const auto& cascade : cascades) { + auto const& v0fromAOD = cascade.v0(); + if (v0fromAOD.posTrackId() == v0.posTrackId && + v0fromAOD.negTrackId() == v0.negTrackId && + cascade.bachelorId() == bachelorTrackIndex.globalId) { + // this will override type, but not collision index + // N.B.: collision index checks still desirable! + currentCascadeEntry.found = true; + currentCascadeEntry.globalId = cascade.globalIndex(); + break; + } + } + if (cascadeBuilderOpts.mc_findableDetachedCascade.value || currentCascadeEntry.collisionId >= 0) { + cascadeList.push_back(currentCascadeEntry); + } + } + } // end bachelorTrackArray loop + } // end v0List loop + + // at this stage, cascadeList is alright, but the v0 indices are still not + // correct. We'll have to loop over all V0s and find the appropriate matches + // ---> but only in mode 1, and only for AO2D-native V0s + if (baseOpts.mc_findableMode.value == 1) { + for (size_t casci = 0; casci < cascadeListReconstructedSize; casci++) { + // loop over v0List to find corresponding v0 index, but do it in sorted way + for (size_t v0i = 0; v0i < v0List.size(); v0i++) { + auto v0 = v0List[sorted_v0[v0i]]; + if (cascadeList[casci].posTrackId == v0.posTrackId && + cascadeList[casci].negTrackId == v0.negTrackId) { + cascadeList[casci].v0Id = v0i; // fix, point to correct V0 index + break; + } + } + } + } + // we should now be done! collect statistics + histos.fill(HIST("hFindableStatistics"), 3.0, cascades.size()); + histos.fill(HIST("hFindableStatistics"), 4.0, cascadeList.size()); + histos.fill(HIST("hFindableStatistics"), 5.0, collisionLessCascades); + + } // end findable mode check + } // end soa::is_table + + // we need to allow for sorted use of cascadeList + sorted_cascade.clear(); + sorted_cascade = sort_indices(cascadeList, (baseOpts.mc_findableMode.value > 0)); + } + + LOGF(info, "AO2D input: %i V0s, %i cascades. Building list sizes: %i V0s, %i cascades", v0s.size(), cascades.size(), v0List.size(), cascadeList.size()); + } + + //__________________________________________________ + template + void markV0sUsedInCascades(TV0s const& v0s, TCascades const& cascades, TTrackedCascades const& trackedCascades) + { + int v0sUsedInCascades = 0; + v0sFromCascades.clear(); + v0Map.clear(); + v0Map.resize(v0List.size(), -2); // marks not used + if (baseOpts.useV0BufferForCascades.value == false) { + return; // don't attempt to mark needlessly + } + if (baseOpts.mEnabledTables[kStoredCascCores]) { + for (const auto& cascade : cascadeList) { + if (cascade.v0Id < 0) + continue; + if (v0Map[cascade.v0Id] == -2) { + v0sUsedInCascades++; + } + v0Map[cascade.v0Id] = -1; // marks used (but isn't the index of a properly built V0, which would be >= 0) + } + } + int trackedCascadeCount = 0; + if constexpr (soa::is_table) { + // tracked only created outside of findable mode + if (baseOpts.mEnabledTables[kStoredTraCascCores] && baseOpts.mc_findableMode.value == 0) { + trackedCascadeCount = trackedCascades.size(); + for (const auto& trackedCascade : trackedCascades) { + auto const& cascade = trackedCascade.cascade(); + if (v0Map[ao2dV0toV0List[cascade.v0Id()]] == -2) { + v0sUsedInCascades++; + } + v0Map[ao2dV0toV0List[cascade.v0Id()]] = -1; // marks used (but isn't the index of a built V0, which would be >= 0) + } + } + } + LOGF(debug, "V0 total %i, Cascade total %i, Tracked cascade total %i, V0s flagged used in cascades: %i", v0s.size(), cascades.size(), trackedCascadeCount, v0sUsedInCascades); + } + + //__________________________________________________ + template + void buildV0s(THistoRegistry& histos, TCollisions const& collisions, TV0s const& v0s, TTracks const& tracks, TMCParticles const& mcParticles, TProducts& products) + { + // prepare MC containers (not necessarily used) + std::vector mcV0infos; // V0MCCore information + std::vector mcParticleIsReco; + + if constexpr (soa::is_table) { + // do this if provided with a mcParticle table as well + mcParticleIsReco.resize(mcParticles.size(), false); + } + + int nV0s = 0; + // Loops over all V0s in the time frame + histos.fill(HIST("hInputStatistics"), kV0CoresBase, v0s.size()); + for (size_t iv0 = 0; iv0 < v0List.size(); iv0++) { + const auto& v0 = v0List[sorted_v0[iv0]]; + + if (!baseOpts.mEnabledTables[kV0CoresBase] && v0Map[iv0] == -2) { + // this v0 hasn't been used by cascades and we're not generating V0s, so skip it + products.v0dataLink(-1, -1); + continue; + } + + // Get tracks and generate candidate + // if collisionId positive: get vertex, negative: origin + // could be replaced by mean vertex (but without much benefit...) + float pvX = 0.0f, pvY = 0.0f, pvZ = 0.0f; + if (v0.collisionId >= 0) { + auto const& collision = collisions.rawIteratorAt(v0.collisionId); + pvX = collision.posX(); + pvY = collision.posY(); + pvZ = collision.posZ(); + } + auto const& posTrack = tracks.rawIteratorAt(v0.posTrackId); + auto const& negTrack = tracks.rawIteratorAt(v0.negTrackId); + + auto posTrackPar = getTrackParCov(posTrack); + auto negTrackPar = getTrackParCov(negTrack); + + // handle TPC-only tracks properly (photon conversions) + if (v0BuilderOpts.moveTPCOnlyTracks) { + bool isPosTPCOnly = (posTrack.hasTPC() && !posTrack.hasITS() && !posTrack.hasTRD() && !posTrack.hasTOF()); + if (isPosTPCOnly) { + // Nota bene: positive is TPC-only -> this entire V0 merits treatment as photon candidate + posTrackPar.setPID(o2::track::PID::Electron); + negTrackPar.setPID(o2::track::PID::Electron); + + auto const& collision = collisions.rawIteratorAt(v0.collisionId); + if (!mVDriftMgr.moveTPCTrack(collision, posTrack, posTrackPar)) { + return; + } + } + + bool isNegTPCOnly = (negTrack.hasTPC() && !negTrack.hasITS() && !negTrack.hasTRD() && !negTrack.hasTOF()); + if (isNegTPCOnly) { + // Nota bene: negative is TPC-only -> this entire V0 merits treatment as photon candidate + posTrackPar.setPID(o2::track::PID::Electron); + negTrackPar.setPID(o2::track::PID::Electron); + + auto const& collision = collisions.rawIteratorAt(v0.collisionId); + if (!mVDriftMgr.moveTPCTrack(collision, negTrack, negTrackPar)) { + return; + } + } + } + + if (!straHelper.buildV0Candidate(v0.collisionId, pvX, pvY, pvZ, posTrack, negTrack, posTrackPar, negTrackPar, v0.isCollinearV0, baseOpts.mEnabledTables[kV0Covs], true)) { + products.v0dataLink(-1, -1); + continue; + } + if constexpr (requires { posTrack.tpcNSigmaEl(); }) { + if (preSelectOpts.preselectOnlyDesiredV0s) { + float lPt = RecoDecay::sqrtSumOfSquares( + straHelper.v0.positiveMomentum[0] + straHelper.v0.negativeMomentum[0], + straHelper.v0.positiveMomentum[1] + straHelper.v0.negativeMomentum[1]); + + float lPtot = RecoDecay::sqrtSumOfSquares( + straHelper.v0.positiveMomentum[0] + straHelper.v0.negativeMomentum[0], + straHelper.v0.positiveMomentum[1] + straHelper.v0.negativeMomentum[1], + straHelper.v0.positiveMomentum[2] + straHelper.v0.negativeMomentum[2]); + + float lLengthTraveled = RecoDecay::sqrtSumOfSquares( + straHelper.v0.position[0] - pvX, + straHelper.v0.position[1] - pvY, + straHelper.v0.position[2] - pvZ); + + uint8_t maskV0Preselection = 0; + + if ( // photon PID, mass, lifetime selection + std::abs(posTrack.tpcNSigmaEl()) < preSelectOpts.maxTPCpidNsigma && + std::abs(negTrack.tpcNSigmaEl()) < preSelectOpts.maxTPCpidNsigma && + std::abs(straHelper.v0.massGamma) < preSelectOpts.massCutPhoton) { + BITSET(maskV0Preselection, selGamma); + } + + if ( // K0Short PID, mass, lifetime selection + std::abs(posTrack.tpcNSigmaPi()) < preSelectOpts.maxTPCpidNsigma && + std::abs(negTrack.tpcNSigmaPi()) < preSelectOpts.maxTPCpidNsigma && + o2::constants::physics::MassKaonNeutral * lLengthTraveled / (lPtot + 1e-13) < preSelectOpts.lifetimeCut->get("lifetimeCutK0S") && + std::abs(straHelper.v0.massK0Short - o2::constants::physics::MassKaonNeutral) < preSelectOpts.massWindownumberOfSigmas * getMassSigmaK0Short(lPt) + preSelectOpts.massWindowSafetyMargin) { + BITSET(maskV0Preselection, selK0Short); + } + + if ( // Lambda PID, mass, lifetime selection + std::abs(posTrack.tpcNSigmaPr()) < preSelectOpts.maxTPCpidNsigma && + std::abs(negTrack.tpcNSigmaPi()) < preSelectOpts.maxTPCpidNsigma && + o2::constants::physics::MassLambda * lLengthTraveled / (lPtot + 1e-13) < preSelectOpts.lifetimeCut->get("lifetimeCutLambda") && + std::abs(straHelper.v0.massLambda - o2::constants::physics::MassLambda) < preSelectOpts.massWindownumberOfSigmas * getMassSigmaLambda(lPt) + preSelectOpts.massWindowSafetyMargin) { + BITSET(maskV0Preselection, selLambda); + } + + if ( // antiLambda PID, mass, lifetime selection + std::abs(posTrack.tpcNSigmaPi()) < preSelectOpts.maxTPCpidNsigma && + std::abs(negTrack.tpcNSigmaPr()) < preSelectOpts.maxTPCpidNsigma && + o2::constants::physics::MassLambda * lLengthTraveled / (lPtot + 1e-13) < preSelectOpts.lifetimeCut->get("lifetimeCutLambda") && + std::abs(straHelper.v0.massAntiLambda - o2::constants::physics::MassLambda) < preSelectOpts.massWindownumberOfSigmas * getMassSigmaLambda(lPt) + preSelectOpts.massWindowSafetyMargin) { + BITSET(maskV0Preselection, selAntiLambda); + } + + histos.fill(HIST("hPreselectionV0s"), maskV0Preselection); + + if (maskV0Preselection == 0) { + products.v0dataLink(-1, -1); + continue; + } + } + } + if (v0Map[iv0] == -1 && baseOpts.useV0BufferForCascades) { + v0Map[iv0] = v0sFromCascades.size(); // provide actual valid index in buffer + v0sFromCascades.push_back(straHelper.v0); + } + // fill requested cursors only if type is not 0 + if (v0.v0Type == 1 || (v0.v0Type > 1 && v0BuilderOpts.generatePhotonCandidates)) { + nV0s++; + if (baseOpts.mEnabledTables[kV0Indices]) { + // for referencing (especially - but not only - when using derived data) + products.v0indices(v0.posTrackId, v0.negTrackId, + v0.collisionId, iv0); + histos.fill(HIST("hTableBuildingStatistics"), kV0Indices); + } + if (baseOpts.mEnabledTables[kV0TrackXs]) { + // further decay chains may need this + products.v0trackXs(straHelper.v0.positiveTrackX, straHelper.v0.negativeTrackX); + histos.fill(HIST("hTableBuildingStatistics"), kV0TrackXs); + } + if (baseOpts.mEnabledTables[kV0CoresBase]) { + // standard analysis + products.v0cores(straHelper.v0.position[0], straHelper.v0.position[1], straHelper.v0.position[2], + straHelper.v0.positiveMomentum[0], straHelper.v0.positiveMomentum[1], straHelper.v0.positiveMomentum[2], + straHelper.v0.negativeMomentum[0], straHelper.v0.negativeMomentum[1], straHelper.v0.negativeMomentum[2], + straHelper.v0.daughterDCA, + straHelper.v0.positiveDCAxy, + straHelper.v0.negativeDCAxy, + TMath::Cos(straHelper.v0.pointingAngle), + straHelper.v0.dcaToPV, + v0.v0Type); + products.v0dataLink(products.v0cores.lastIndex(), -1); + histos.fill(HIST("hTableBuildingStatistics"), kV0CoresBase); + } + if (baseOpts.mEnabledTables[kV0TraPosAtDCAs]) { + // for tracking studies + products.v0dauPositions(straHelper.v0.positivePosition[0], straHelper.v0.positivePosition[1], straHelper.v0.positivePosition[2], + straHelper.v0.negativePosition[0], straHelper.v0.negativePosition[1], straHelper.v0.negativePosition[2]); + histos.fill(HIST("hTableBuildingStatistics"), kV0TraPosAtDCAs); + } + if (baseOpts.mEnabledTables[kV0TraPosAtIUs]) { + // for tracking studies + std::array positivePositionIU; + std::array negativePositionIU; + o2::track::TrackPar positiveTrackParam = getTrackPar(posTrack); + o2::track::TrackPar negativeTrackParam = getTrackPar(negTrack); + positiveTrackParam.getXYZGlo(positivePositionIU); + negativeTrackParam.getXYZGlo(negativePositionIU); + products.v0dauPositionsIU(positivePositionIU[0], positivePositionIU[1], positivePositionIU[2], + negativePositionIU[0], negativePositionIU[1], negativePositionIU[2]); + histos.fill(HIST("hTableBuildingStatistics"), kV0TraPosAtIUs); + } + if (baseOpts.mEnabledTables[kV0Covs]) { + products.v0covs(straHelper.v0.positionCovariance, straHelper.v0.momentumCovariance); + histos.fill(HIST("hTableBuildingStatistics"), kV0Covs); + } + + //_________________________________________________________ + // MC handling part + if constexpr (soa::is_table) { + // only worry about this if someone else worried about this + if ((baseOpts.mEnabledTables[kV0MCCores] || baseOpts.mEnabledTables[kMcV0Labels] || baseOpts.mEnabledTables[kV0MCCollRefs])) { + thisInfo.label = -1; + thisInfo.motherLabel = -1; + thisInfo.pdgCode = 0; + thisInfo.pdgCodeMother = 0; + thisInfo.pdgCodePositive = 0; + thisInfo.pdgCodeNegative = 0; + thisInfo.mcCollision = -1; + thisInfo.xyz[0] = thisInfo.xyz[1] = thisInfo.xyz[2] = 0.0f; + thisInfo.posP[0] = thisInfo.posP[1] = thisInfo.posP[2] = 0.0f; + thisInfo.negP[0] = thisInfo.negP[1] = thisInfo.negP[2] = 0.0f; + thisInfo.momentum[0] = thisInfo.momentum[1] = thisInfo.momentum[2] = 0.0f; + + // Association check + // There might be smarter ways of doing this in the future + if (negTrack.has_mcParticle() && posTrack.has_mcParticle()) { + auto lMCNegTrack = negTrack.template mcParticle_as(); + auto lMCPosTrack = posTrack.template mcParticle_as(); + + thisInfo.pdgCodePositive = lMCPosTrack.pdgCode(); + thisInfo.pdgCodeNegative = lMCNegTrack.pdgCode(); + thisInfo.processPositive = lMCPosTrack.getProcess(); + thisInfo.processNegative = lMCNegTrack.getProcess(); + thisInfo.posP[0] = lMCPosTrack.px(); + thisInfo.posP[1] = lMCPosTrack.py(); + thisInfo.posP[2] = lMCPosTrack.pz(); + thisInfo.negP[0] = lMCNegTrack.px(); + thisInfo.negP[1] = lMCNegTrack.py(); + thisInfo.negP[2] = lMCNegTrack.pz(); + + // check for pi -> mu + antineutrino decay + // if present, de-reference original V0 correctly and provide label to original object + // NOTA BENE: the prong info will still correspond to a muon, treat carefully! + int negOriginating = -1, posOriginating = -1, particleForDecayPositionIdx = -1; + negOriginating = getOriginatingParticle(lMCNegTrack, particleForDecayPositionIdx, v0BuilderOpts.mc_treatPiToMuDecays); + posOriginating = getOriginatingParticle(lMCPosTrack, particleForDecayPositionIdx, v0BuilderOpts.mc_treatPiToMuDecays); + + if (negOriginating > -1 && negOriginating == posOriginating) { + auto originatingV0 = mcParticles.rawIteratorAt(negOriginating); + auto particleForDecayPosition = mcParticles.rawIteratorAt(particleForDecayPositionIdx); + + thisInfo.label = originatingV0.globalIndex(); + thisInfo.xyz[0] = particleForDecayPosition.vx(); + thisInfo.xyz[1] = particleForDecayPosition.vy(); + thisInfo.xyz[2] = particleForDecayPosition.vz(); + + if (originatingV0.has_mcCollision()) { + thisInfo.mcCollision = originatingV0.mcCollisionId(); // save this reference, please + } + + // acquire information + thisInfo.pdgCode = originatingV0.pdgCode(); + thisInfo.isPhysicalPrimary = originatingV0.isPhysicalPrimary(); + thisInfo.momentum[0] = originatingV0.px(); + thisInfo.momentum[1] = originatingV0.py(); + thisInfo.momentum[2] = originatingV0.pz(); + + if (originatingV0.has_mothers()) { + for (const auto& lV0Mother : originatingV0.template mothers_as()) { + thisInfo.pdgCodeMother = lV0Mother.pdgCode(); + thisInfo.motherLabel = lV0Mother.globalIndex(); + } + } + } + + } // end association check + // Construct label table (note: this will be joinable with V0Datas!) + if (baseOpts.mEnabledTables[kMcV0Labels]) { + products.v0labels(thisInfo.label, thisInfo.motherLabel); + histos.fill(HIST("hTableBuildingStatistics"), kMcV0Labels); + } + + // Construct found tag + if (baseOpts.mEnabledTables[kV0FoundTags]) { + products.v0FoundTag(v0.found); + histos.fill(HIST("hTableBuildingStatistics"), kV0FoundTags); + } + + // Mark mcParticle as recoed (no searching necessary afterwards) + if (thisInfo.label > -1) { + mcParticleIsReco[thisInfo.label] = true; + } + + // ---] Symmetric populate [--- + // in this approach, V0Cores will be joinable with V0MCCores. + // this is the most pedagogical approach, but it is also more limited + // and it might use more disk space unnecessarily. + if (v0BuilderOpts.mc_populateV0MCCoresSymmetric) { + if (baseOpts.mEnabledTables[kV0MCCores]) { + products.v0mccores( + thisInfo.label, thisInfo.pdgCode, + thisInfo.pdgCodeMother, thisInfo.pdgCodePositive, thisInfo.pdgCodeNegative, + thisInfo.isPhysicalPrimary, thisInfo.xyz[0], thisInfo.xyz[1], thisInfo.xyz[2], + thisInfo.posP[0], thisInfo.posP[1], thisInfo.posP[2], + thisInfo.negP[0], thisInfo.negP[1], thisInfo.negP[2], + thisInfo.momentum[0], thisInfo.momentum[1], thisInfo.momentum[2]); + histos.fill(HIST("hTableBuildingStatistics"), kV0MCCores); + histos.fill(HIST("hPrimaryV0s"), 0); + if (thisInfo.isPhysicalPrimary) + histos.fill(HIST("hPrimaryV0s"), 1); + } + if (baseOpts.mEnabledTables[kV0MCCollRefs]) { + products.v0mccollref(thisInfo.mcCollision); + histos.fill(HIST("hTableBuildingStatistics"), kV0MCCollRefs); + } + + // n.b. placing the interlink index here allows for the writing of + // code that is agnostic with respect to the joinability of + // V0Cores and V0MCCores (always dereference -> safe) + if (baseOpts.mEnabledTables[kV0CoreMCLabels]) { + products.v0CoreMCLabels(iv0); // interlink index + histos.fill(HIST("hTableBuildingStatistics"), kV0CoreMCLabels); + } + } + // ---] Asymmetric populate [--- + // in this approach, V0Cores will NOT be joinable with V0MCCores. + // an additional reference to V0MCCore that IS joinable with V0Cores + // will be provided to the user. + if (v0BuilderOpts.mc_populateV0MCCoresAsymmetric) { + int thisV0MCCoreIndex = -1; + // step 1: check if this element is already provided in the table + // using the packedIndices variable calculated above + for (uint32_t ii = 0; ii < mcV0infos.size(); ii++) { + if (thisInfo.label == mcV0infos[ii].label && mcV0infos[ii].label > -1) { + thisV0MCCoreIndex = ii; + break; // this exists already in list + } + } + if (thisV0MCCoreIndex < 0 && thisInfo.label > -1) { + // this V0MCCore does not exist yet. Create it and reference it + thisV0MCCoreIndex = mcV0infos.size(); + mcV0infos.push_back(thisInfo); + } + if (baseOpts.mEnabledTables[kV0CoreMCLabels]) { + products.v0CoreMCLabels(thisV0MCCoreIndex); // interlink index + histos.fill(HIST("hTableBuildingStatistics"), kV0CoreMCLabels); + } + } + } // enabled tables check + } // constexpr requires check + } + } + + // finish populating V0MCCores if in asymmetric mode + if constexpr (soa::is_table) { + if (v0BuilderOpts.mc_populateV0MCCoresAsymmetric && (baseOpts.mEnabledTables[kV0MCCores] || baseOpts.mEnabledTables[kV0MCCollRefs])) { + // first step: add any un-recoed v0mmcores that were requested + for (const auto& mcParticle : mcParticles) { + thisInfo.label = -1; + thisInfo.motherLabel = -1; + thisInfo.pdgCode = 0; + thisInfo.pdgCodeMother = -1; + thisInfo.pdgCodePositive = -1; + thisInfo.pdgCodeNegative = -1; + thisInfo.mcCollision = -1; + thisInfo.xyz[0] = thisInfo.xyz[1] = thisInfo.xyz[2] = 0.0f; + thisInfo.posP[0] = thisInfo.posP[1] = thisInfo.posP[2] = 0.0f; + thisInfo.negP[0] = thisInfo.negP[1] = thisInfo.negP[2] = 0.0f; + thisInfo.momentum[0] = thisInfo.momentum[1] = thisInfo.momentum[2] = 0.0f; + + if (mcParticleIsReco[mcParticle.globalIndex()] == true) + continue; // skip if already created in list + + if (std::fabs(mcParticle.y()) > v0BuilderOpts.mc_rapidityWindow) + continue; // skip outside midrapidity + + if (v0BuilderOpts.mc_keepOnlyPhysicalPrimary && !mcParticle.isPhysicalPrimary()) + continue; // skip secondary MC V0s + + if ( + (v0BuilderOpts.mc_addGeneratedK0Short && mcParticle.pdgCode() == 310) || + (v0BuilderOpts.mc_addGeneratedLambda && mcParticle.pdgCode() == 3122) || + (v0BuilderOpts.mc_addGeneratedAntiLambda && mcParticle.pdgCode() == -3122) || + (v0BuilderOpts.mc_addGeneratedGamma && mcParticle.pdgCode() == 22)) { + thisInfo.pdgCode = mcParticle.pdgCode(); + thisInfo.isPhysicalPrimary = mcParticle.isPhysicalPrimary(); + thisInfo.label = mcParticle.globalIndex(); + + if (mcParticle.has_mcCollision()) { + thisInfo.mcCollision = mcParticle.mcCollisionId(); // save this reference, please + } + + // + thisInfo.momentum[0] = mcParticle.px(); + thisInfo.momentum[1] = mcParticle.py(); + thisInfo.momentum[2] = mcParticle.pz(); + + if (mcParticle.has_mothers()) { + auto const& mother = mcParticle.template mothers_first_as(); + thisInfo.pdgCodeMother = mother.pdgCode(); + thisInfo.motherLabel = mother.globalIndex(); + } + if (mcParticle.has_daughters()) { + auto const& daughters = mcParticle.template daughters_as(); + + for (const auto& dau : daughters) { + if (dau.getProcess() != 4) + continue; + + if (dau.pdgCode() > 0) { + thisInfo.pdgCodePositive = dau.pdgCode(); + thisInfo.processPositive = dau.getProcess(); + thisInfo.posP[0] = dau.px(); + thisInfo.posP[1] = dau.py(); + thisInfo.posP[2] = dau.pz(); + thisInfo.xyz[0] = dau.vx(); + thisInfo.xyz[1] = dau.vy(); + thisInfo.xyz[2] = dau.vz(); + } + if (dau.pdgCode() < 0) { + thisInfo.pdgCodeNegative = dau.pdgCode(); + thisInfo.processNegative = dau.getProcess(); + thisInfo.negP[0] = dau.px(); + thisInfo.negP[1] = dau.py(); + thisInfo.negP[2] = dau.pz(); + } + } + } + + // if I got here, it means this MC particle was not recoed and is of interest. Add it please + mcV0infos.push_back(thisInfo); + } + } + + for (const auto& info : mcV0infos) { + if (baseOpts.mEnabledTables[kV0MCCores]) { + products.v0mccores( + info.label, info.pdgCode, + info.pdgCodeMother, info.pdgCodePositive, info.pdgCodeNegative, + info.isPhysicalPrimary, info.xyz[0], info.xyz[1], info.xyz[2], + info.posP[0], info.posP[1], info.posP[2], + info.negP[0], info.negP[1], info.negP[2], + info.momentum[0], info.momentum[1], info.momentum[2]); + histos.fill(HIST("hTableBuildingStatistics"), kV0MCCores); + histos.fill(HIST("hPrimaryV0s"), 0); + if (info.isPhysicalPrimary) + histos.fill(HIST("hPrimaryV0s"), 1); + } + if (baseOpts.mEnabledTables[kV0MCCollRefs]) { + products.v0mccollref(info.mcCollision); + histos.fill(HIST("hTableBuildingStatistics"), kV0MCCollRefs); + } + } + } // end V0MCCores filling in case of MC + } // end constexpr requires mcParticles + + LOGF(debug, "V0s in DF: %i, V0s built: %i, V0s built and buffered for cascades: %i.", v0s.size(), nV0s, v0sFromCascades.size()); + } + + //__________________________________________________ + template + void extractMonteCarloProperties(TTrack const& posTrack, TTrack const& negTrack, TTrack const& bachTrack, TMCParticles const& mcParticles) + { + // encapsulates acquisition of MC properties from MC + thisCascInfo.pdgCode = -1, thisCascInfo.pdgCodeMother = -1; + thisCascInfo.pdgCodePositive = -1, thisCascInfo.pdgCodeNegative = -1; + thisCascInfo.pdgCodeBachelor = -1, thisCascInfo.pdgCodeV0 = -1; + thisCascInfo.isPhysicalPrimary = false; + thisCascInfo.xyz[0] = -999.0f, thisCascInfo.xyz[1] = -999.0f, thisCascInfo.xyz[2] = -999.0f; + thisCascInfo.lxyz[0] = -999.0f, thisCascInfo.lxyz[1] = -999.0f, thisCascInfo.lxyz[2] = -999.0f; + thisCascInfo.posP[0] = -999.0f, thisCascInfo.posP[1] = -999.0f, thisCascInfo.posP[2] = -999.0f; + thisCascInfo.negP[0] = -999.0f, thisCascInfo.negP[1] = -999.0f, thisCascInfo.negP[2] = -999.0f; + thisCascInfo.bachP[0] = -999.0f, thisCascInfo.bachP[1] = -999.0f, thisCascInfo.bachP[2] = -999.0f; + thisCascInfo.momentum[0] = -999.0f, thisCascInfo.momentum[1] = -999.0f, thisCascInfo.momentum[2] = -999.0f; + thisCascInfo.label = -1, thisCascInfo.motherLabel = -1; + thisCascInfo.mcParticlePositive = -1; + thisCascInfo.mcParticleNegative = -1; + thisCascInfo.mcParticleBachelor = -1; + + // Association check + // There might be smarter ways of doing this in the future + if (negTrack.has_mcParticle() && posTrack.has_mcParticle() && bachTrack.has_mcParticle()) { + auto lMCBachTrack = bachTrack.template mcParticle_as(); + auto lMCNegTrack = negTrack.template mcParticle_as(); + auto lMCPosTrack = posTrack.template mcParticle_as(); + + thisCascInfo.mcParticlePositive = lMCPosTrack.globalIndex(); + thisCascInfo.mcParticleNegative = lMCNegTrack.globalIndex(); + thisCascInfo.mcParticleBachelor = lMCBachTrack.globalIndex(); + thisCascInfo.pdgCodePositive = lMCPosTrack.pdgCode(); + thisCascInfo.pdgCodeNegative = lMCNegTrack.pdgCode(); + thisCascInfo.pdgCodeBachelor = lMCBachTrack.pdgCode(); + thisCascInfo.posP[0] = lMCPosTrack.px(); + thisCascInfo.posP[1] = lMCPosTrack.py(); + thisCascInfo.posP[2] = lMCPosTrack.pz(); + thisCascInfo.negP[0] = lMCNegTrack.px(); + thisCascInfo.negP[1] = lMCNegTrack.py(); + thisCascInfo.negP[2] = lMCNegTrack.pz(); + thisCascInfo.bachP[0] = lMCBachTrack.px(); + thisCascInfo.bachP[1] = lMCBachTrack.py(); + thisCascInfo.bachP[2] = lMCBachTrack.pz(); + thisCascInfo.processPositive = lMCPosTrack.getProcess(); + thisCascInfo.processNegative = lMCNegTrack.getProcess(); + thisCascInfo.processBachelor = lMCBachTrack.getProcess(); + + // Step 0: treat pi -> mu + antineutrino + // if present, de-reference original V0 correctly and provide label to original object + // NOTA BENE: the prong info will still correspond to a muon, treat carefully! + int negOriginating = -1, posOriginating = -1, bachOriginating = -1; + int particleForLambdaDecayPositionIdx = -1, particleForCascadeDecayPositionIdx = -1; + negOriginating = getOriginatingParticle(lMCNegTrack, particleForLambdaDecayPositionIdx, cascadeBuilderOpts.mc_treatPiToMuDecays); + posOriginating = getOriginatingParticle(lMCPosTrack, particleForLambdaDecayPositionIdx, cascadeBuilderOpts.mc_treatPiToMuDecays); + bachOriginating = getOriginatingParticle(lMCBachTrack, particleForCascadeDecayPositionIdx, cascadeBuilderOpts.mc_treatPiToMuDecays); + + if (negOriginating > -1 && negOriginating == posOriginating) { + auto originatingV0 = mcParticles.rawIteratorAt(negOriginating); + auto particleForLambdaDecayPosition = mcParticles.rawIteratorAt(particleForLambdaDecayPositionIdx); + + thisCascInfo.label = originatingV0.globalIndex(); + thisCascInfo.lxyz[0] = particleForLambdaDecayPosition.vx(); + thisCascInfo.lxyz[1] = particleForLambdaDecayPosition.vy(); + thisCascInfo.lxyz[2] = particleForLambdaDecayPosition.vz(); + thisCascInfo.pdgCodeV0 = originatingV0.pdgCode(); + + if (originatingV0.has_mothers()) { + for (const auto& lV0Mother : originatingV0.template mothers_as()) { + if (lV0Mother.globalIndex() == bachOriginating) { // found mother particle + thisCascInfo.label = lV0Mother.globalIndex(); + + if (lV0Mother.has_mcCollision()) { + thisCascInfo.mcCollision = lV0Mother.mcCollisionId(); // save this reference, please + } + + thisCascInfo.pdgCode = lV0Mother.pdgCode(); + thisCascInfo.isPhysicalPrimary = lV0Mother.isPhysicalPrimary(); + thisCascInfo.xyz[0] = originatingV0.vx(); + thisCascInfo.xyz[1] = originatingV0.vy(); + thisCascInfo.xyz[2] = originatingV0.vz(); + thisCascInfo.momentum[0] = lV0Mother.px(); + thisCascInfo.momentum[1] = lV0Mother.py(); + thisCascInfo.momentum[2] = lV0Mother.pz(); + if (lV0Mother.has_mothers()) { + for (const auto& lV0GrandMother : lV0Mother.template mothers_as()) { + thisCascInfo.pdgCodeMother = lV0GrandMother.pdgCode(); + thisCascInfo.motherLabel = lV0GrandMother.globalIndex(); + } + } + } + } // end v0 mother loop + } // end has_mothers check for V0 + } // end conditional of pos/neg originating being the same + } // end association check + } + + //__________________________________________________ + template + void buildCascades(THistoRegistry& histos, TCollisions const& collisions, TCascades const& cascades, TTracks const& tracks, TMCParticles const& mcParticles, TProducts& products) + { + // prepare MC containers (not necessarily used) + std::vector mcCascinfos; // V0MCCore information + std::vector mcParticleIsReco; + + if constexpr (soa::is_table) { + // do this if provided with a mcParticle table as well + mcParticleIsReco.resize(mcParticles.size(), false); + } + + if (!baseOpts.mEnabledTables[kStoredCascCores]) { + return; // don't do if no request for cascades in place + } + int nCascades = 0; + // Loops over all cascades in the time frame + histos.fill(HIST("hInputStatistics"), kStoredCascCores, cascades.size()); + for (size_t icascade = 0; icascade < cascades.size(); icascade++) { + // Get tracks and generate candidate + auto const& cascade = cascades[sorted_cascade[icascade]]; + // if collisionId positive: get vertex, negative: origin + // could be replaced by mean vertex (but without much benefit...) + float pvX = 0.0f, pvY = 0.0f, pvZ = 0.0f; + if (cascade.collisionId >= 0) { + auto const& collision = collisions.rawIteratorAt(cascade.collisionId); + pvX = collision.posX(); + pvY = collision.posY(); + pvZ = collision.posZ(); + } + auto const& posTrack = tracks.rawIteratorAt(cascade.posTrackId); + auto const& negTrack = tracks.rawIteratorAt(cascade.negTrackId); + auto const& bachTrack = tracks.rawIteratorAt(cascade.bachTrackId); + if (baseOpts.useV0BufferForCascades) { + // this processing path uses a buffer of V0s so that no + // additional minimization step is redone. It consumes less + // CPU at the cost of more memory. Since memory is a more + // limited commodity, this isn't the default option. + + // check if cached - if not, skip + if (cascade.v0Id < 0 || v0Map[cascade.v0Id] < 0) { + // this V0 hasn't been stored / cached + products.cascdataLink(-1); + interlinks.cascadeToCascCores.push_back(-1); + continue; // didn't work out, skip + } + + if (!straHelper.buildCascadeCandidate(cascade.collisionId, pvX, pvY, pvZ, + v0sFromCascades[v0Map[cascade.v0Id]], + posTrack, + negTrack, + bachTrack, + baseOpts.mEnabledTables[kCascBBs], + cascadeBuilderOpts.useCascadeMomentumAtPrimVtx, + baseOpts.mEnabledTables[kCascCovs])) { + products.cascdataLink(-1); + interlinks.cascadeToCascCores.push_back(-1); + continue; // didn't work out, skip + } + } else { + // this processing path generates the entire cascade + // from tracks, without any need to have V0s generated. + if (!straHelper.buildCascadeCandidate(cascade.collisionId, pvX, pvY, pvZ, + posTrack, + negTrack, + bachTrack, + baseOpts.mEnabledTables[kCascBBs], + cascadeBuilderOpts.useCascadeMomentumAtPrimVtx, + baseOpts.mEnabledTables[kCascCovs])) { + products.cascdataLink(-1); + interlinks.cascadeToCascCores.push_back(-1); + continue; // didn't work out, skip + } + } + nCascades++; + + if constexpr (requires { posTrack.tpcNSigmaEl(); }) { + if (preSelectOpts.preselectOnlyDesiredCascades) { + float lPt = RecoDecay::sqrtSumOfSquares( + straHelper.cascade.bachelorMomentum[0] + straHelper.cascade.positiveMomentum[0] + straHelper.cascade.negativeMomentum[0], + straHelper.cascade.bachelorMomentum[1] + straHelper.cascade.positiveMomentum[1] + straHelper.cascade.negativeMomentum[1]); + + float lPtot = RecoDecay::sqrtSumOfSquares( + straHelper.cascade.bachelorMomentum[0] + straHelper.cascade.positiveMomentum[0] + straHelper.cascade.negativeMomentum[0], + straHelper.cascade.bachelorMomentum[1] + straHelper.cascade.positiveMomentum[1] + straHelper.cascade.negativeMomentum[1], + straHelper.cascade.bachelorMomentum[2] + straHelper.cascade.positiveMomentum[2] + straHelper.cascade.negativeMomentum[2]); + + float lV0Ptot = RecoDecay::sqrtSumOfSquares( + straHelper.cascade.positiveMomentum[0] + straHelper.cascade.negativeMomentum[0], + straHelper.cascade.positiveMomentum[1] + straHelper.cascade.negativeMomentum[1], + straHelper.cascade.positiveMomentum[2] + straHelper.cascade.negativeMomentum[2]); + + float lLengthTraveled = RecoDecay::sqrtSumOfSquares( + straHelper.cascade.cascadePosition[0] - pvX, + straHelper.cascade.cascadePosition[1] - pvY, + straHelper.cascade.cascadePosition[2] - pvZ); + + float lV0LengthTraveled = RecoDecay::sqrtSumOfSquares( + straHelper.cascade.v0Position[0] - straHelper.cascade.cascadePosition[0], + straHelper.cascade.v0Position[1] - straHelper.cascade.cascadePosition[1], + straHelper.cascade.v0Position[2] - straHelper.cascade.cascadePosition[2]); + + uint8_t maskCascadePreselection = 0; + + if ( // XiMinus PID and mass selection + straHelper.cascade.charge < 0 && + std::abs(posTrack.tpcNSigmaPr()) < preSelectOpts.maxTPCpidNsigma && + std::abs(negTrack.tpcNSigmaPi()) < preSelectOpts.maxTPCpidNsigma && + std::abs(bachTrack.tpcNSigmaPi()) < preSelectOpts.maxTPCpidNsigma && + o2::constants::physics::MassLambda * lV0LengthTraveled / (lV0Ptot + 1e-13) < preSelectOpts.lifetimeCut->get("lifetimeCutLambda") && + o2::constants::physics::MassXiMinus * lLengthTraveled / (lPtot + 1e-13) < preSelectOpts.lifetimeCut->get("lifetimeCutXi") && + std::abs(straHelper.cascade.massXi - o2::constants::physics::MassXiMinus) < preSelectOpts.massWindownumberOfSigmas * getMassSigmaXi(lPt) + preSelectOpts.massWindowSafetyMargin) { + BITSET(maskCascadePreselection, selXiMinus); + } + + if ( // XiPlus PID and mass selection + straHelper.cascade.charge > 0 && + std::abs(posTrack.tpcNSigmaPi()) < preSelectOpts.maxTPCpidNsigma && + std::abs(negTrack.tpcNSigmaPr()) < preSelectOpts.maxTPCpidNsigma && + std::abs(bachTrack.tpcNSigmaPi()) < preSelectOpts.maxTPCpidNsigma && + o2::constants::physics::MassLambda * lV0LengthTraveled / (lV0Ptot + 1e-13) < preSelectOpts.lifetimeCut->get("lifetimeCutLambda") && + o2::constants::physics::MassXiMinus * lLengthTraveled / (lPtot + 1e-13) < preSelectOpts.lifetimeCut->get("lifetimeCutXi") && + std::abs(straHelper.cascade.massXi - o2::constants::physics::MassXiMinus) < preSelectOpts.massWindownumberOfSigmas * getMassSigmaXi(lPt) + preSelectOpts.massWindowSafetyMargin) { + BITSET(maskCascadePreselection, selXiPlus); + } + + if ( // OmegaMinus PID and mass selection + straHelper.cascade.charge < 0 && + std::abs(posTrack.tpcNSigmaPr()) < preSelectOpts.maxTPCpidNsigma && + std::abs(negTrack.tpcNSigmaPi()) < preSelectOpts.maxTPCpidNsigma && + std::abs(bachTrack.tpcNSigmaKa()) < preSelectOpts.maxTPCpidNsigma && + o2::constants::physics::MassLambda * lV0LengthTraveled / (lV0Ptot + 1e-13) < preSelectOpts.lifetimeCut->get("lifetimeCutLambda") && + o2::constants::physics::MassOmegaMinus * lLengthTraveled / (lPtot + 1e-13) < preSelectOpts.lifetimeCut->get("lifetimeCutOmega") && + std::abs(straHelper.cascade.massOmega - o2::constants::physics::MassOmegaMinus) < preSelectOpts.massWindownumberOfSigmas * getMassSigmaOmega(lPt) + preSelectOpts.massWindowSafetyMargin) { + BITSET(maskCascadePreselection, selOmegaMinus); + } + + if ( // OmegaPlus PID and mass selection + straHelper.cascade.charge > 0 && + std::abs(posTrack.tpcNSigmaPi()) < preSelectOpts.maxTPCpidNsigma && + std::abs(negTrack.tpcNSigmaPr()) < preSelectOpts.maxTPCpidNsigma && + std::abs(bachTrack.tpcNSigmaKa()) < preSelectOpts.maxTPCpidNsigma && + o2::constants::physics::MassLambda * lV0LengthTraveled / (lV0Ptot + 1e-13) < preSelectOpts.lifetimeCut->get("lifetimeCutLambda") && + o2::constants::physics::MassOmegaMinus * lLengthTraveled / (lPtot + 1e-13) < preSelectOpts.lifetimeCut->get("lifetimeCutOmega") && + std::abs(straHelper.cascade.massOmega - o2::constants::physics::MassOmegaMinus) < preSelectOpts.massWindownumberOfSigmas * getMassSigmaOmega(lPt) + preSelectOpts.massWindowSafetyMargin) { + BITSET(maskCascadePreselection, selOmegaPlus); + } + + histos.fill(HIST("hPreselectionCascades"), maskCascadePreselection); + + if (maskCascadePreselection == 0) { + products.cascdataLink(-1); + interlinks.cascadeToCascCores.push_back(-1); + continue; + } + } + } + + // generate analysis tables as required + if (baseOpts.mEnabledTables[kCascIndices]) { + products.cascidx(cascade.globalId, + straHelper.cascade.positiveTrack, straHelper.cascade.negativeTrack, + straHelper.cascade.bachelorTrack, straHelper.cascade.collisionId); + histos.fill(HIST("hTableBuildingStatistics"), kCascIndices); + } + if (baseOpts.mEnabledTables[kStoredCascCores]) { + products.cascdata(straHelper.cascade.charge, straHelper.cascade.massXi, straHelper.cascade.massOmega, + straHelper.cascade.cascadePosition[0], straHelper.cascade.cascadePosition[1], straHelper.cascade.cascadePosition[2], + straHelper.cascade.v0Position[0], straHelper.cascade.v0Position[1], straHelper.cascade.v0Position[2], + straHelper.cascade.positiveMomentum[0], straHelper.cascade.positiveMomentum[1], straHelper.cascade.positiveMomentum[2], + straHelper.cascade.negativeMomentum[0], straHelper.cascade.negativeMomentum[1], straHelper.cascade.negativeMomentum[2], + straHelper.cascade.bachelorMomentum[0], straHelper.cascade.bachelorMomentum[1], straHelper.cascade.bachelorMomentum[2], + straHelper.cascade.cascadeMomentum[0], straHelper.cascade.cascadeMomentum[1], straHelper.cascade.cascadeMomentum[2], + straHelper.cascade.v0DaughterDCA, straHelper.cascade.cascadeDaughterDCA, + straHelper.cascade.positiveDCAxy, straHelper.cascade.negativeDCAxy, + straHelper.cascade.bachelorDCAxy, straHelper.cascade.cascadeDCAxy, straHelper.cascade.cascadeDCAz); + histos.fill(HIST("hTableBuildingStatistics"), kStoredCascCores); + + // interlink always produced if cascades generated + products.cascdataLink(products.cascdata.lastIndex()); + interlinks.cascCoreToCascades.push_back(cascade.globalId); + interlinks.cascadeToCascCores.push_back(products.cascdata.lastIndex()); + } + + if (baseOpts.mEnabledTables[kCascTrackXs]) { + products.cascTrackXs(straHelper.cascade.positiveTrackX, straHelper.cascade.negativeTrackX, straHelper.cascade.bachelorTrackX); + histos.fill(HIST("hTableBuildingStatistics"), kCascTrackXs); + } + if (baseOpts.mEnabledTables[kCascBBs]) { + products.cascbb(straHelper.cascade.bachBaryonCosPA, straHelper.cascade.bachBaryonDCAxyToPV); + histos.fill(HIST("hTableBuildingStatistics"), kCascBBs); + } + if (baseOpts.mEnabledTables[kCascCovs]) { + products.casccovs(straHelper.cascade.covariance); + histos.fill(HIST("hTableBuildingStatistics"), kCascCovs); + } + + //_________________________________________________________ + // MC handling part + if constexpr (soa::is_table) { + // only worry about this if someone else worried about this + if ((baseOpts.mEnabledTables[kCascMCCores] || baseOpts.mEnabledTables[kMcCascLabels] || baseOpts.mEnabledTables[kCascMCCollRefs])) { + extractMonteCarloProperties(posTrack, negTrack, bachTrack, mcParticles); + + // Construct label table (note: this will be joinable with CascDatas) + if (baseOpts.mEnabledTables[kMcCascLabels]) { + products.casclabels( + thisCascInfo.label, thisCascInfo.motherLabel); + histos.fill(HIST("hTableBuildingStatistics"), kMcCascLabels); + } + + // Construct found tag + if (baseOpts.mEnabledTables[kCascFoundTags]) { + products.cascFoundTag(cascade.found); + histos.fill(HIST("hTableBuildingStatistics"), kCascFoundTags); + } + + // Mark mcParticle as recoed (no searching necessary afterwards) + if (thisCascInfo.label > -1) { + mcParticleIsReco[thisCascInfo.label] = true; + } + + if (cascadeBuilderOpts.mc_populateCascMCCoresSymmetric) { + if (baseOpts.mEnabledTables[kCascMCCores]) { + products.cascmccores( + thisCascInfo.pdgCode, thisCascInfo.pdgCodeMother, thisCascInfo.pdgCodeV0, thisCascInfo.isPhysicalPrimary, + thisCascInfo.pdgCodePositive, thisCascInfo.pdgCodeNegative, thisCascInfo.pdgCodeBachelor, + thisCascInfo.xyz[0], thisCascInfo.xyz[1], thisCascInfo.xyz[2], + thisCascInfo.lxyz[0], thisCascInfo.lxyz[1], thisCascInfo.lxyz[2], + thisCascInfo.posP[0], thisCascInfo.posP[1], thisCascInfo.posP[2], + thisCascInfo.negP[0], thisCascInfo.negP[1], thisCascInfo.negP[2], + thisCascInfo.bachP[0], thisCascInfo.bachP[1], thisCascInfo.bachP[2], + thisCascInfo.momentum[0], thisCascInfo.momentum[1], thisCascInfo.momentum[2]); + histos.fill(HIST("hTableBuildingStatistics"), kCascMCCores); + } + if (baseOpts.mEnabledTables[kCascMCCollRefs]) { + products.cascmccollrefs(thisCascInfo.mcCollision); + histos.fill(HIST("hTableBuildingStatistics"), kCascMCCollRefs); + } + } + + if (cascadeBuilderOpts.mc_populateCascMCCoresAsymmetric) { + int thisCascMCCoreIndex = -1; + // step 1: check if this element is already provided in the table + // using the packedIndices variable calculated above + for (uint32_t ii = 0; ii < mcCascinfos.size(); ii++) { + if (thisCascInfo.label == mcCascinfos[ii].label && mcCascinfos[ii].label > -1) { + thisCascMCCoreIndex = ii; + break; // this exists already in list + } + } + if (thisCascMCCoreIndex < 0) { + // this CascMCCore does not exist yet. Create it and reference it + thisCascMCCoreIndex = mcCascinfos.size(); + mcCascinfos.push_back(thisCascInfo); + } + if (baseOpts.mEnabledTables[kCascCoreMCLabels]) { + products.cascCoreMClabels(thisCascMCCoreIndex); // interlink: reconstructed -> MC index + histos.fill(HIST("hTableBuildingStatistics"), kCascCoreMCLabels); + } + } + + } // enabled tables check + + // if BB tags requested, generate them now + if (baseOpts.mEnabledTables[kMcCascBBTags]) { + bool bbTag = false; + if (bachTrack.has_mcParticle()) { + auto bachelorParticle = bachTrack.template mcParticle_as(); + if (bachelorParticle.pdgCode() == 211) { // pi+, look for antiproton in negative prong + if (negTrack.has_mcParticle()) { + auto baryonParticle = negTrack.template mcParticle_as(); + if (baryonParticle.has_mothers() && bachelorParticle.has_mothers() && baryonParticle.pdgCode() == -2212) { + for (const auto& baryonMother : baryonParticle.template mothers_as()) { + for (const auto& pionMother : bachelorParticle.template mothers_as()) { + if (baryonMother.globalIndex() == pionMother.globalIndex() && baryonMother.pdgCode() == -3122) { + bbTag = true; + } + } + } + } + } + } // end if-pion + if (bachelorParticle.pdgCode() == -211) { // pi-, look for proton in positive prong + if (posTrack.has_mcParticle()) { + auto baryonParticle = posTrack.template mcParticle_as(); + if (baryonParticle.has_mothers() && bachelorParticle.has_mothers() && baryonParticle.pdgCode() == 2212) { + for (const auto& baryonMother : baryonParticle.template mothers_as()) { + for (const auto& pionMother : bachelorParticle.template mothers_as()) { + if (baryonMother.globalIndex() == pionMother.globalIndex() && baryonMother.pdgCode() == 3122) { + bbTag = true; + } + } + } + } + } + } // end if-pion + } // end bachelor has mcparticle + // Construct label table (note: this will be joinable with CascDatas) + products.bbtags(bbTag); + histos.fill(HIST("hTableBuildingStatistics"), kMcCascBBTags); + } // end BB tag table enabled check + + } // constexpr requires mcParticles check + } // cascades loop + + //_________________________________________________________ + // MC handling part + if constexpr (soa::is_table) { + if ((baseOpts.mEnabledTables[kCascMCCores] || baseOpts.mEnabledTables[kMcCascLabels] || baseOpts.mEnabledTables[kCascMCCollRefs])) { + // now populate V0MCCores if in asymmetric mode + if (cascadeBuilderOpts.mc_populateCascMCCoresAsymmetric) { + // first step: add any un-recoed v0mmcores that were requested + for (const auto& mcParticle : mcParticles) { + thisCascInfo.pdgCode = -1, thisCascInfo.pdgCodeMother = -1; + thisCascInfo.pdgCodePositive = -1, thisCascInfo.pdgCodeNegative = -1; + thisCascInfo.pdgCodeBachelor = -1, thisCascInfo.pdgCodeV0 = -1; + thisCascInfo.isPhysicalPrimary = false; + thisCascInfo.xyz[0] = 0.0f, thisCascInfo.xyz[1] = 0.0f, thisCascInfo.xyz[2] = 0.0f; + thisCascInfo.lxyz[0] = 0.0f, thisCascInfo.lxyz[1] = 0.0f, thisCascInfo.lxyz[2] = 0.0f; + thisCascInfo.posP[0] = 0.0f, thisCascInfo.posP[1] = 0.0f, thisCascInfo.posP[2] = 0.0f; + thisCascInfo.negP[0] = 0.0f, thisCascInfo.negP[1] = 0.0f, thisCascInfo.negP[2] = 0.0f; + thisCascInfo.bachP[0] = 0.0f, thisCascInfo.bachP[1] = 0.0f, thisCascInfo.bachP[2] = 0.0f; + thisCascInfo.momentum[0] = 0.0f, thisCascInfo.momentum[1] = 0.0f, thisCascInfo.momentum[2] = 0.0f; + thisCascInfo.label = -1, thisCascInfo.motherLabel = -1; + thisCascInfo.mcParticlePositive = -1; + thisCascInfo.mcParticleNegative = -1; + thisCascInfo.mcParticleBachelor = -1; + + if (mcParticleIsReco[mcParticle.globalIndex()] == true) + continue; // skip if already created in list + + if (std::fabs(mcParticle.y()) > cascadeBuilderOpts.mc_rapidityWindow) + continue; // skip outside midrapidity + + if (cascadeBuilderOpts.mc_keepOnlyPhysicalPrimary && !mcParticle.isPhysicalPrimary()) + continue; // skip secondary MC cascades + + if ( + (cascadeBuilderOpts.mc_addGeneratedXiMinus && mcParticle.pdgCode() == 3312) || + (cascadeBuilderOpts.mc_addGeneratedXiPlus && mcParticle.pdgCode() == -3312) || + (cascadeBuilderOpts.mc_addGeneratedOmegaMinus && mcParticle.pdgCode() == 3334) || + (cascadeBuilderOpts.mc_addGeneratedOmegaPlus && mcParticle.pdgCode() == -3334)) { + thisCascInfo.pdgCode = mcParticle.pdgCode(); + thisCascInfo.isPhysicalPrimary = mcParticle.isPhysicalPrimary(); + + if (mcParticle.has_mcCollision()) { + thisCascInfo.mcCollision = mcParticle.mcCollisionId(); // save this reference, please + } + thisCascInfo.momentum[0] = mcParticle.px(); + thisCascInfo.momentum[1] = mcParticle.py(); + thisCascInfo.momentum[2] = mcParticle.pz(); + thisCascInfo.label = mcParticle.globalIndex(); + + if (mcParticle.has_daughters()) { + auto const& daughters = mcParticle.template daughters_as(); + for (const auto& dau : daughters) { + if (dau.getProcess() != 4) // check whether the daughter comes from a decay + continue; + + if (std::abs(dau.pdgCode()) == 211 || std::abs(dau.pdgCode()) == 321) { + thisCascInfo.pdgCodeBachelor = dau.pdgCode(); + thisCascInfo.bachP[0] = dau.px(); + thisCascInfo.bachP[1] = dau.py(); + thisCascInfo.bachP[2] = dau.pz(); + thisCascInfo.xyz[0] = dau.vx(); + thisCascInfo.xyz[1] = dau.vy(); + thisCascInfo.xyz[2] = dau.vz(); + thisCascInfo.mcParticleBachelor = dau.globalIndex(); + } + if (std::abs(dau.pdgCode()) == 2212) { + thisCascInfo.pdgCodeV0 = dau.pdgCode(); + + for (const auto& v0Dau : dau.template daughters_as()) { + if (v0Dau.getProcess() != 4) + continue; + + if (v0Dau.pdgCode() > 0) { + thisCascInfo.pdgCodePositive = v0Dau.pdgCode(); + thisCascInfo.processPositive = v0Dau.getProcess(); + thisCascInfo.posP[0] = v0Dau.px(); + thisCascInfo.posP[1] = v0Dau.py(); + thisCascInfo.posP[2] = v0Dau.pz(); + thisCascInfo.lxyz[0] = v0Dau.vx(); + thisCascInfo.lxyz[1] = v0Dau.vy(); + thisCascInfo.lxyz[2] = v0Dau.vz(); + thisCascInfo.mcParticlePositive = v0Dau.globalIndex(); + } + if (v0Dau.pdgCode() < 0) { + thisCascInfo.pdgCodeNegative = v0Dau.pdgCode(); + thisCascInfo.processNegative = v0Dau.getProcess(); + thisCascInfo.negP[0] = v0Dau.px(); + thisCascInfo.negP[1] = v0Dau.py(); + thisCascInfo.negP[2] = v0Dau.pz(); + thisCascInfo.mcParticleNegative = v0Dau.globalIndex(); + } + } + } + } + } + + // if I got here, it means this MC particle was not recoed and is of interest. Add it please + mcCascinfos.push_back(thisCascInfo); + } + } + + for (const auto& thisInfoToFill : mcCascinfos) { + if (baseOpts.mEnabledTables[kCascMCCores]) { + products.cascmccores( // a lot of the info below will be compressed in case of not-recoed MC (good!) + thisInfoToFill.pdgCode, thisInfoToFill.pdgCodeMother, thisInfoToFill.pdgCodeV0, thisInfoToFill.isPhysicalPrimary, + thisInfoToFill.pdgCodePositive, thisInfoToFill.pdgCodeNegative, thisInfoToFill.pdgCodeBachelor, + thisInfoToFill.xyz[0], thisInfoToFill.xyz[1], thisInfoToFill.xyz[2], + thisInfoToFill.lxyz[0], thisInfoToFill.lxyz[1], thisInfoToFill.lxyz[2], + thisInfoToFill.posP[0], thisInfoToFill.posP[1], thisInfoToFill.posP[2], + thisInfoToFill.negP[0], thisInfoToFill.negP[1], thisInfoToFill.negP[2], + thisInfoToFill.bachP[0], thisInfoToFill.bachP[1], thisInfoToFill.bachP[2], + thisInfoToFill.momentum[0], thisInfoToFill.momentum[1], thisInfoToFill.momentum[2]); + histos.fill(HIST("hTableBuildingStatistics"), kCascMCCores); + } + if (baseOpts.mEnabledTables[kCascMCCollRefs]) { + products.cascmccollrefs(thisInfoToFill.mcCollision); + histos.fill(HIST("hTableBuildingStatistics"), kCascMCCollRefs); + } + } + } + } // enabled tables check + } // constexpr requires mcParticles check + + LOGF(debug, "Cascades in DF: %i, cascades built: %i", cascades.size(), nCascades); + } + + //__________________________________________________ + template + void buildKFCascades(THistoRegistry& histos, TCollisions const& collisions, TCascades const& cascades, TTracks const& tracks, TMCParticles const& mcParticles, TProducts& products) + { + if (!baseOpts.mEnabledTables[kStoredKFCascCores]) { + return; // don't do if no request for cascades in place + } + int nCascades = 0; + // Loops over all cascades in the time frame + histos.fill(HIST("hInputStatistics"), kStoredKFCascCores, cascades.size()); + for (size_t icascade = 0; icascade < cascades.size(); icascade++) { + // Get tracks and generate candidate + auto const& cascade = cascades[sorted_cascade[icascade]]; + // if collisionId positive: get vertex, negative: origin + // could be replaced by mean vertex (but without much benefit...) + float pvX = 0.0f, pvY = 0.0f, pvZ = 0.0f; + if (cascade.collisionId >= 0) { + auto const& collision = collisions.rawIteratorAt(cascade.collisionId); + pvX = collision.posX(); + pvY = collision.posY(); + pvZ = collision.posZ(); + } + auto const& posTrack = tracks.rawIteratorAt(cascade.posTrackId); + auto const& negTrack = tracks.rawIteratorAt(cascade.negTrackId); + auto const& bachTrack = tracks.rawIteratorAt(cascade.bachTrackId); + if (!straHelper.buildCascadeCandidateWithKF(cascade.collisionId, pvX, pvY, pvZ, + posTrack, + negTrack, + bachTrack, + baseOpts.mEnabledTables[kCascBBs], + cascadeBuilderOpts.kfConstructMethod, + cascadeBuilderOpts.kfTuneForOmega, + cascadeBuilderOpts.kfUseV0MassConstraint, + cascadeBuilderOpts.kfUseCascadeMassConstraint, + cascadeBuilderOpts.kfDoDCAFitterPreMinimV0, + cascadeBuilderOpts.kfDoDCAFitterPreMinimCasc)) { + products.kfcascdataLink(-1); + interlinks.cascadeToKFCascCores.push_back(-1); + continue; // didn't work out, skip + } + nCascades++; + + // generate analysis tables as required + if (baseOpts.mEnabledTables[kKFCascIndices]) { + products.kfcascidx(cascade.globalId, + straHelper.cascade.positiveTrack, straHelper.cascade.negativeTrack, + straHelper.cascade.bachelorTrack, straHelper.cascade.collisionId); + histos.fill(HIST("hTableBuildingStatistics"), kKFCascIndices); + } + if (baseOpts.mEnabledTables[kStoredKFCascCores]) { + products.kfcascdata(straHelper.cascade.charge, straHelper.cascade.massXi, straHelper.cascade.massOmega, + straHelper.cascade.cascadePosition[0], straHelper.cascade.cascadePosition[1], straHelper.cascade.cascadePosition[2], + straHelper.cascade.v0Position[0], straHelper.cascade.v0Position[1], straHelper.cascade.v0Position[2], + straHelper.cascade.positivePosition[0], straHelper.cascade.positivePosition[1], straHelper.cascade.positivePosition[2], + straHelper.cascade.negativePosition[0], straHelper.cascade.negativePosition[1], straHelper.cascade.negativePosition[2], + straHelper.cascade.positiveMomentum[0], straHelper.cascade.positiveMomentum[1], straHelper.cascade.positiveMomentum[2], + straHelper.cascade.negativeMomentum[0], straHelper.cascade.negativeMomentum[1], straHelper.cascade.negativeMomentum[2], + straHelper.cascade.bachelorMomentum[0], straHelper.cascade.bachelorMomentum[1], straHelper.cascade.bachelorMomentum[2], + straHelper.cascade.v0Momentum[0], straHelper.cascade.v0Momentum[1], straHelper.cascade.v0Momentum[2], + straHelper.cascade.cascadeMomentum[0], straHelper.cascade.cascadeMomentum[1], straHelper.cascade.cascadeMomentum[2], + straHelper.cascade.v0DaughterDCA, straHelper.cascade.cascadeDaughterDCA, + straHelper.cascade.positiveDCAxy, straHelper.cascade.negativeDCAxy, + straHelper.cascade.bachelorDCAxy, straHelper.cascade.cascadeDCAxy, straHelper.cascade.cascadeDCAz, + straHelper.cascade.kfMLambda, straHelper.cascade.kfV0Chi2, straHelper.cascade.kfCascadeChi2); + histos.fill(HIST("hTableBuildingStatistics"), kStoredKFCascCores); + + // interlink always produced if cascades generated + products.kfcascdataLink(products.kfcascdata.lastIndex()); + interlinks.kfCascCoreToCascades.push_back(cascade.globalId); + interlinks.cascadeToKFCascCores.push_back(products.kfcascdata.lastIndex()); + } + if (baseOpts.mEnabledTables[kKFCascCovs]) { + products.kfcasccovs(straHelper.cascade.covariance, straHelper.cascade.kfTrackCovarianceV0, straHelper.cascade.kfTrackCovariancePos, straHelper.cascade.kfTrackCovarianceNeg); + histos.fill(HIST("hTableBuildingStatistics"), kKFCascCovs); + } + + //_________________________________________________________ + // MC handling part (labels only) + if constexpr (soa::is_table) { + // only worry about this if someone else worried about this + if ((baseOpts.mEnabledTables[kMcKFCascLabels])) { + extractMonteCarloProperties(posTrack, negTrack, bachTrack, mcParticles); + + // Construct label table (note: this will be joinable with KFCascDatas) + products.kfcasclabels(thisCascInfo.label); + histos.fill(HIST("hTableBuildingStatistics"), kMcKFCascLabels); + } // enabled tables check + } // constexpr requires mcParticles check + } // end loop over cascades + + LOGF(debug, "KF Cascades in DF: %i, KF cascades built: %i", cascades.size(), nCascades); + } + + //__________________________________________________ + template + void buildTrackedCascades(THistoRegistry& histos, TCollisions const& collisions, TStrangeTracks const& cascadeTracks, TMCParticles const& mcParticles, TProducts& products) + { + if (!baseOpts.mEnabledTables[kStoredTraCascCores] || baseOpts.mc_findableMode.value != 0) { + return; // don't do if no request for cascades in place or findable mode used + } + int nCascades = 0; + // Loops over all V0s in the time frame + histos.fill(HIST("hInputStatistics"), kStoredTraCascCores, cascadeTracks.size()); + for (const auto& cascadeTrack : cascadeTracks) { + // Get tracks and generate candidate + if (!cascadeTrack.has_track()) + continue; // safety (should be fine but depends on future stratrack dev) + + auto const& strangeTrack = cascadeTrack.template track_as(); + + // if collisionId positive: get vertex, negative: origin + // could be replaced by mean vertex (but without much benefit...) + float pvX = 0.0f, pvY = 0.0f, pvZ = 0.0f; + if (strangeTrack.has_collision()) { + auto const& collision = collisions.rawIteratorAt(strangeTrack.collisionId()); + pvX = collision.posX(); + pvY = collision.posY(); + pvZ = collision.posZ(); + } + auto const& cascade = cascadeTrack.cascade(); + auto const& v0 = cascade.v0(); + auto const& posTrack = v0.template posTrack_as(); + auto const& negTrack = v0.template negTrack_as(); + auto const& bachTrack = cascade.template bachelor_as(); + if (!straHelper.buildCascadeCandidate(strangeTrack.collisionId(), pvX, pvY, pvZ, + posTrack, + negTrack, + bachTrack, + baseOpts.mEnabledTables[kCascBBs], + cascadeBuilderOpts.useCascadeMomentumAtPrimVtx, + baseOpts.mEnabledTables[kCascCovs])) { + products.tracascdataLink(-1); + interlinks.cascadeToTraCascCores.push_back(-1); + continue; // didn't work out, skip + } + + // recalculate DCAxy, DCAz with strange track + auto strangeTrackParCov = getTrackParCov(strangeTrack); + std::array dcaInfo; + strangeTrackParCov.setPID(o2::track::PID::XiMinus); // FIXME: not OK for omegas + o2::base::Propagator::Instance()->propagateToDCABxByBz({pvX, pvY, pvZ}, strangeTrackParCov, 2.f, straHelper.fitter.getMatCorrType(), &dcaInfo); + straHelper.cascade.cascadeDCAxy = dcaInfo[0]; + straHelper.cascade.cascadeDCAz = dcaInfo[1]; + + // get momentum from strange track (should not be very different) + strangeTrackParCov.getPxPyPzGlo(straHelper.cascade.cascadeMomentum); + + // accounting + nCascades++; + + // generate analysis tables as required + if (baseOpts.mEnabledTables[kTraCascIndices]) { + products.tracascidx(cascade.globalIndex(), + straHelper.cascade.positiveTrack, straHelper.cascade.negativeTrack, + straHelper.cascade.bachelorTrack, cascadeTrack.trackId(), straHelper.cascade.collisionId); + histos.fill(HIST("hTableBuildingStatistics"), kTraCascIndices); + } + if (baseOpts.mEnabledTables[kStoredTraCascCores]) { + products.tracascdata(straHelper.cascade.charge, cascadeTrack.xiMass(), cascadeTrack.omegaMass(), + cascadeTrack.decayX(), cascadeTrack.decayY(), cascadeTrack.decayZ(), + straHelper.cascade.v0Position[0], straHelper.cascade.v0Position[1], straHelper.cascade.v0Position[2], + straHelper.cascade.positiveMomentum[0], straHelper.cascade.positiveMomentum[1], straHelper.cascade.positiveMomentum[2], + straHelper.cascade.negativeMomentum[0], straHelper.cascade.negativeMomentum[1], straHelper.cascade.negativeMomentum[2], + straHelper.cascade.bachelorMomentum[0], straHelper.cascade.bachelorMomentum[1], straHelper.cascade.bachelorMomentum[2], + straHelper.cascade.cascadeMomentum[0], straHelper.cascade.cascadeMomentum[1], straHelper.cascade.cascadeMomentum[2], + straHelper.cascade.v0DaughterDCA, straHelper.cascade.cascadeDaughterDCA, + straHelper.cascade.positiveDCAxy, straHelper.cascade.negativeDCAxy, + straHelper.cascade.bachelorDCAxy, straHelper.cascade.cascadeDCAxy, straHelper.cascade.cascadeDCAz, + cascadeTrack.matchingChi2(), cascadeTrack.topologyChi2(), cascadeTrack.itsClsSize()); + histos.fill(HIST("hTableBuildingStatistics"), kStoredTraCascCores); + + // interlink always produced if base core table generated + products.tracascdataLink(products.tracascdata.lastIndex()); + interlinks.traCascCoreToCascades.push_back(cascade.globalIndex()); + interlinks.cascadeToTraCascCores.push_back(products.tracascdata.lastIndex()); + } + if (baseOpts.mEnabledTables[kCascCovs]) { + std::array traCovMat = {0.}; + strangeTrackParCov.getCovXYZPxPyPzGlo(traCovMat); + float traCovMatArray[21]; + for (int ii = 0; ii < 21; ii++) { + traCovMatArray[ii] = traCovMat[ii]; + } + products.tracasccovs(traCovMatArray); + histos.fill(HIST("hTableBuildingStatistics"), kCascCovs); + } + + //_________________________________________________________ + // MC handling part (labels only) + if constexpr (soa::is_table) { + // only worry about this if someone else worried about this + if ((baseOpts.mEnabledTables[kMcTraCascLabels])) { + extractMonteCarloProperties(posTrack, negTrack, bachTrack, mcParticles); + + // Construct label table (note: this will be joinable with KFCascDatas) + products.tracasclabels(thisCascInfo.label); + histos.fill(HIST("hTableBuildingStatistics"), kMcTraCascLabels); + } // enabled tables check + } // constexpr requires mcParticles check + } // end loop over cascades + LOGF(debug, "Tracked cascades in DF: %i, tracked cascades built: %i", cascadeTracks.size(), nCascades); + } + + //__________________________________________________ + // MC kink handling + template + int getOriginatingParticle(mcpart const& part, int& indexForPositionOfDecay, bool treatPiToMuDecays) + { + int returnValue = -1; + if (part.has_mothers()) { + auto const& motherList = part.template mothers_as(); + if (motherList.size() == 1) { + for (const auto& mother : motherList) { + if (std::abs(part.pdgCode()) == 13 && treatPiToMuDecays) { + // muon decay, de-ref mother twice + if (mother.has_mothers()) { + auto grandMotherList = mother.template mothers_as(); + if (grandMotherList.size() == 1) { + for (const auto& grandMother : grandMotherList) { + returnValue = grandMother.globalIndex(); + indexForPositionOfDecay = mother.globalIndex(); // for V0 decay position: grab muon + } + } + } + } else { + returnValue = mother.globalIndex(); + indexForPositionOfDecay = part.globalIndex(); + } + } + } + } + return returnValue; + } + + //__________________________________________________ + template + void dataProcess(TCCDB& ccdb, THistoRegistry& histos, TCollisions const& collisions, TMCCollisions const& mccollisions, TV0s const& v0s, TCascades const& cascades, TTrackedCascades const& trackedCascades, TTracks const& tracks, TBCs const& bcs, TMCParticles const& mcParticles, TProducts& products) + { + if (!initCCDB(ccdb, bcs, collisions)) + return; + + // reset vectors for cascade interlinks + resetInterlinks(); + + // prepare v0List, cascadeList + prepareBuildingLists(histos, collisions, mccollisions, v0s, cascades, tracks, mcParticles); + + // mark V0s that will be buffered for the cascade building + markV0sUsedInCascades(v0List, cascadeList, trackedCascades); + + // build V0s + buildV0s(histos, collisions, v0List, tracks, mcParticles, products); + + // build cascades + buildCascades(histos, collisions, cascadeList, tracks, mcParticles, products); + buildKFCascades(histos, collisions, cascadeList, tracks, mcParticles, products); + + // build tracked cascades only if subscription is Run 3 like (doesn't exist in Run 2) + if constexpr (soa::is_table) { + buildTrackedCascades(histos, collisions, trackedCascades, mcParticles, products); + } + + populateCascadeInterlinks(); + } +}; // end BuilderModule + +} // namespace strangenessbuilder +} // namespace pwglf +} // namespace o2 + +#endif // PWGLF_UTILS_STRANGENESSBUILDERMODULE_H_ From c18f387dbc6f7d6862c26e4e284094df13a08877 Mon Sep 17 00:00:00 2001 From: ALICE Builder Date: Fri, 18 Jul 2025 18:37:16 +0200 Subject: [PATCH 02/12] Please consider the following formatting changes (#440) --- .../TableProducer/trackPropagationTester.cxx | 3 +- .../Strangeness/propagationService.cxx | 49 ++++++++++--------- PWGLF/Utils/strangenessBuilderModule.h | 22 +++++---- 3 files changed, 41 insertions(+), 33 deletions(-) diff --git a/Common/TableProducer/trackPropagationTester.cxx b/Common/TableProducer/trackPropagationTester.cxx index 817c3f59db5..c0a84ee8bad 100644 --- a/Common/TableProducer/trackPropagationTester.cxx +++ b/Common/TableProducer/trackPropagationTester.cxx @@ -105,4 +105,5 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { WorkflowSpec workflow{adaptAnalysisTask(cfgc)}; return workflow; -} \ No newline at end of file +} + \ No newline at end of file diff --git a/PWGLF/TableProducer/Strangeness/propagationService.cxx b/PWGLF/TableProducer/Strangeness/propagationService.cxx index 11522ff1525..1b42d60409d 100644 --- a/PWGLF/TableProducer/Strangeness/propagationService.cxx +++ b/PWGLF/TableProducer/Strangeness/propagationService.cxx @@ -10,40 +10,42 @@ // or submit itself to any jurisdiction. /// \file propagationService.cxx -/// \brief +/// \brief /// \author ALICE //=============================================================== // -// Merged track propagation + strangeness building task -// -// Provides a common task to deal with track propagation and -// strangeness building in a single DPL device that is particularly -// adequate for pipelining. +// Merged track propagation + strangeness building task // -// Currently meant for testing and performance evaluation +// Provides a common task to deal with track propagation and +// strangeness building in a single DPL device that is particularly +// adequate for pipelining. +// +// Currently meant for testing and performance evaluation // //=============================================================== -#include "Framework/AnalysisDataModel.h" -#include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" -#include "Framework/RunningWorkflowInfo.h" -#include "Common/DataModel/TrackSelectionTables.h" +#include "PWGLF/Utils/strangenessBuilderModule.h" + #include "Common/Core/trackUtilities.h" -#include "ReconstructionDataFormats/DCA.h" -#include "DetectorsBase/Propagator.h" -#include "DetectorsBase/GeometryManager.h" -#include "CommonUtils/NameConf.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Common/Tools/StandardCCDBLoader.h" +#include "Common/Tools/TrackPropagationModule.h" + +#include "CCDB/BasicCCDBManager.h" #include "CCDB/CcdbApi.h" +#include "CommonConstants/GeomConstants.h" +#include "CommonUtils/NameConf.h" +#include "DataFormatsCalibration/MeanVertexObject.h" #include "DataFormatsParameters/GRPMagField.h" -#include "CCDB/BasicCCDBManager.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/Propagator.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" #include "Framework/HistogramRegistry.h" -#include "DataFormatsCalibration/MeanVertexObject.h" -#include "CommonConstants/GeomConstants.h" -#include "PWGLF/Utils/strangenessBuilderModule.h" -#include "Common/Tools/TrackPropagationModule.h" -#include "Common/Tools/StandardCCDBLoader.h" +#include "Framework/RunningWorkflowInfo.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/DCA.h" using namespace o2; using namespace o2::framework; @@ -143,4 +145,5 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { WorkflowSpec workflow{adaptAnalysisTask(cfgc)}; return workflow; -} \ No newline at end of file +} + \ No newline at end of file diff --git a/PWGLF/Utils/strangenessBuilderModule.h b/PWGLF/Utils/strangenessBuilderModule.h index 2c982acddee..ce57c8f23e6 100644 --- a/PWGLF/Utils/strangenessBuilderModule.h +++ b/PWGLF/Utils/strangenessBuilderModule.h @@ -19,19 +19,23 @@ #ifndef PWGLF_UTILS_STRANGENESSBUILDERMODULE_H_ #define PWGLF_UTILS_STRANGENESSBUILDERMODULE_H_ -#include -#include -#include -#include -#include -#include "Framework/AnalysisDataModel.h" -#include "Framework/Configurable.h" -#include "Framework/HistogramSpec.h" #include "TableHelper.h" -#include "Common/Core/TPCVDriftManager.h" + #include "PWGLF/DataModel/LFStrangenessTables.h" #include "PWGLF/Utils/strangenessBuilderHelper.h" +#include "Common/Core/TPCVDriftManager.h" + +#include "Framework/AnalysisDataModel.h" +#include "Framework/Configurable.h" +#include "Framework/HistogramSpec.h" + +#include +#include +#include +#include +#include + //__________________________________________ // strangeness builder module From 1357fee5552311e95e1229da8e6139c1c45e8ab8 Mon Sep 17 00:00:00 2001 From: ALICE Builder Date: Fri, 18 Jul 2025 18:41:29 +0200 Subject: [PATCH 03/12] Please consider the following formatting changes (#441) --- Common/TableProducer/trackPropagationTester.cxx | 1 - PWGLF/TableProducer/Strangeness/propagationService.cxx | 1 - 2 files changed, 2 deletions(-) diff --git a/Common/TableProducer/trackPropagationTester.cxx b/Common/TableProducer/trackPropagationTester.cxx index c0a84ee8bad..694ff77d3c0 100644 --- a/Common/TableProducer/trackPropagationTester.cxx +++ b/Common/TableProducer/trackPropagationTester.cxx @@ -106,4 +106,3 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) WorkflowSpec workflow{adaptAnalysisTask(cfgc)}; return workflow; } - \ No newline at end of file diff --git a/PWGLF/TableProducer/Strangeness/propagationService.cxx b/PWGLF/TableProducer/Strangeness/propagationService.cxx index 1b42d60409d..3455cfa64a7 100644 --- a/PWGLF/TableProducer/Strangeness/propagationService.cxx +++ b/PWGLF/TableProducer/Strangeness/propagationService.cxx @@ -146,4 +146,3 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) WorkflowSpec workflow{adaptAnalysisTask(cfgc)}; return workflow; } - \ No newline at end of file From 6665a49322990ddb5e7499bb52a592abf9ebbbc7 Mon Sep 17 00:00:00 2001 From: David Dobrigkeit Chinellato Date: Fri, 18 Jul 2025 15:40:08 -0300 Subject: [PATCH 04/12] Fix headers --- PWGLF/TableProducer/Strangeness/propagationService.cxx | 1 - PWGLF/Utils/strangenessBuilderModule.h | 3 +++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/PWGLF/TableProducer/Strangeness/propagationService.cxx b/PWGLF/TableProducer/Strangeness/propagationService.cxx index 3455cfa64a7..dd7353a594d 100644 --- a/PWGLF/TableProducer/Strangeness/propagationService.cxx +++ b/PWGLF/TableProducer/Strangeness/propagationService.cxx @@ -36,7 +36,6 @@ #include "CCDB/CcdbApi.h" #include "CommonConstants/GeomConstants.h" #include "CommonUtils/NameConf.h" -#include "DataFormatsCalibration/MeanVertexObject.h" #include "DataFormatsParameters/GRPMagField.h" #include "DetectorsBase/GeometryManager.h" #include "DetectorsBase/Propagator.h" diff --git a/PWGLF/Utils/strangenessBuilderModule.h b/PWGLF/Utils/strangenessBuilderModule.h index ce57c8f23e6..70f87e73a10 100644 --- a/PWGLF/Utils/strangenessBuilderModule.h +++ b/PWGLF/Utils/strangenessBuilderModule.h @@ -27,8 +27,11 @@ #include "Common/Core/TPCVDriftManager.h" #include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisHelpers.h" #include "Framework/Configurable.h" #include "Framework/HistogramSpec.h" +#include "Framework/HistogramRegistry.h" +#include "DataFormatsCalibration/MeanVertexObject.h" #include #include From 0434fdd31192fab4656aa7d986f98177ee4c9d84 Mon Sep 17 00:00:00 2001 From: David Dobrigkeit Chinellato Date: Fri, 18 Jul 2025 16:24:12 -0300 Subject: [PATCH 05/12] One more header fix --- PWGLF/Utils/strangenessBuilderModule.h | 1 + 1 file changed, 1 insertion(+) diff --git a/PWGLF/Utils/strangenessBuilderModule.h b/PWGLF/Utils/strangenessBuilderModule.h index 70f87e73a10..b83ca91c146 100644 --- a/PWGLF/Utils/strangenessBuilderModule.h +++ b/PWGLF/Utils/strangenessBuilderModule.h @@ -32,6 +32,7 @@ #include "Framework/HistogramSpec.h" #include "Framework/HistogramRegistry.h" #include "DataFormatsCalibration/MeanVertexObject.h" +#include "DataFormatsParameters/GRPMagField.h" #include #include From 68e906ea0a4d9adc8eda30bbbbb456d46af2792c Mon Sep 17 00:00:00 2001 From: ALICE Builder Date: Fri, 18 Jul 2025 21:28:33 +0200 Subject: [PATCH 06/12] Please consider the following formatting changes (#442) --- PWGLF/Utils/strangenessBuilderModule.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/PWGLF/Utils/strangenessBuilderModule.h b/PWGLF/Utils/strangenessBuilderModule.h index b83ca91c146..f9f370001ab 100644 --- a/PWGLF/Utils/strangenessBuilderModule.h +++ b/PWGLF/Utils/strangenessBuilderModule.h @@ -26,13 +26,13 @@ #include "Common/Core/TPCVDriftManager.h" +#include "DataFormatsCalibration/MeanVertexObject.h" +#include "DataFormatsParameters/GRPMagField.h" #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisHelpers.h" #include "Framework/Configurable.h" -#include "Framework/HistogramSpec.h" #include "Framework/HistogramRegistry.h" -#include "DataFormatsCalibration/MeanVertexObject.h" -#include "DataFormatsParameters/GRPMagField.h" +#include "Framework/HistogramSpec.h" #include #include From c3a9667eeb409ee243108e2fed0ad715e0851bd8 Mon Sep 17 00:00:00 2001 From: David Dobrigkeit Chinellato Date: Sat, 19 Jul 2025 16:21:28 -0300 Subject: [PATCH 07/12] Fully suppress strangenes not required --- PWGLF/Utils/strangenessBuilderModule.h | 85 +++++++++++++++----------- 1 file changed, 51 insertions(+), 34 deletions(-) diff --git a/PWGLF/Utils/strangenessBuilderModule.h b/PWGLF/Utils/strangenessBuilderModule.h index f9f370001ab..27948b692e1 100644 --- a/PWGLF/Utils/strangenessBuilderModule.h +++ b/PWGLF/Utils/strangenessBuilderModule.h @@ -397,6 +397,8 @@ class BuilderModule return preSelectOpts.massCutOm->get("constant") + pt * preSelectOpts.massCutOm->get("linear") + preSelectOpts.massCutOm->get("expoConstant") * TMath::Exp(-pt / preSelectOpts.massCutOm->get("expoRelax")); } + int nEnabledTables = 0; + // helper object o2::pwglf::strangenessBuilderHelper straHelper; @@ -519,6 +521,49 @@ class BuilderModule cascadeBuilderOpts = inputCascadeBuilderOpts; preSelectOpts = inputPreSelectOpts; + baseOpts.mEnabledTables.resize(nTables, 0); + + LOGF(info, "Checking if strangeness building is required"); + auto& workflows = context.services().template get(); + + nEnabledTables = 0; + + TString listOfRequestors[nTables]; + for (int i = 0; i < nTables; i++) { + int f = baseOpts.enabledTables->get(tableNames[i].c_str(), "enable"); + if (f == 1) { + baseOpts.mEnabledTables[i] = 1; + listOfRequestors[i] = "manual enabling"; + } + if (f == -1) { + // autodetect this table in other devices + for (o2::framework::DeviceSpec const& device : workflows.devices) { + // Step 1: check if this device subscribed to the V0data table + for (auto const& input : device.inputs) { + if (o2::framework::DataSpecUtils::partialMatch(input.matcher, o2::header::DataOrigin("AOD"))) { + auto&& [origin, description, version] = o2::framework::DataSpecUtils::asConcreteDataMatcher(input.matcher); + std::string tableNameWithVersion = tableNames[i]; + if (version > 0) { + tableNameWithVersion += Form("_%03d", version); + } + if (input.matcher.binding == tableNameWithVersion) { + LOGF(info, "Device %s has subscribed to %s (version %i)", device.name, tableNames[i], version); + listOfRequestors[i].Append(Form("%s ", device.name.c_str())); + baseOpts.mEnabledTables[i] = 1; + nEnabledTables++; + } + } + } + } + } + } + + if(nEnabledTables==0){ + LOGF(info, "Strangeness building not required. Will suppress all functionality, including logs, from this point forward."); + return; + } + + // setup bookkeeping histogram auto h = histos.template add("hTableBuildingStatistics", "hTableBuildingStatistics", o2::framework::kTH1D, {{nTablesConst, -0.5f, static_cast(nTablesConst)}}); auto h2 = histos.template add("hInputStatistics", "hInputStatistics", o2::framework::kTH1D, {{nTablesConst, -0.5f, static_cast(nTablesConst)}}); @@ -588,44 +633,12 @@ class BuilderModule mRunNumber = 0; - baseOpts.mEnabledTables.resize(nTables, 0); - - LOGF(info, "Configuring tables to generate"); - auto& workflows = context.services().template get(); - - TString listOfRequestors[nTables]; for (int i = 0; i < nTables; i++) { // adjust bookkeeping histogram h->GetXaxis()->SetBinLabel(i + 1, tableNames[i].c_str()); h2->GetXaxis()->SetBinLabel(i + 1, tableNames[i].c_str()); - h->SetBinContent(i + 1, -1); // mark all as disabled to start - - int f = baseOpts.enabledTables->get(tableNames[i].c_str(), "enable"); - if (f == 1) { - baseOpts.mEnabledTables[i] = 1; - listOfRequestors[i] = "manual enabling"; - } - if (f == -1) { - // autodetect this table in other devices - for (o2::framework::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("strangenessbuilder-initializer") == 0) - continue; // don't listen to the initializer - if (o2::framework::DataSpecUtils::partialMatch(input.matcher, o2::header::DataOrigin("AOD"))) { - auto&& [origin, description, version] = o2::framework::DataSpecUtils::asConcreteDataMatcher(input.matcher); - std::string tableNameWithVersion = tableNames[i]; - if (version > 0) { - tableNameWithVersion += Form("_%03d", version); - } - if (input.matcher.binding == tableNameWithVersion) { - LOGF(info, "Device %s has subscribed to %s (version %i)", device.name, tableNames[i], version); - listOfRequestors[i].Append(Form("%s ", device.name.c_str())); - baseOpts.mEnabledTables[i] = 1; - } - } - } - } + if(baseOpts.mEnabledTables[i] == false){ + h->SetBinContent(i + 1, -1); // mark disabled tables, distinguish from zero counts } } @@ -2509,6 +2522,10 @@ class BuilderModule template void dataProcess(TCCDB& ccdb, THistoRegistry& histos, TCollisions const& collisions, TMCCollisions const& mccollisions, TV0s const& v0s, TCascades const& cascades, TTrackedCascades const& trackedCascades, TTracks const& tracks, TBCs const& bcs, TMCParticles const& mcParticles, TProducts& products) { + if(nEnabledTables == 0){ + return; // fully suppressed + } + if (!initCCDB(ccdb, bcs, collisions)) return; From 9dbefd49ac282d36453850ed7a1b500266bd4858 Mon Sep 17 00:00:00 2001 From: ALICE Builder Date: Sun, 20 Jul 2025 09:36:11 +0200 Subject: [PATCH 08/12] Please consider the following formatting changes (#443) --- PWGLF/Utils/strangenessBuilderModule.h | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/PWGLF/Utils/strangenessBuilderModule.h b/PWGLF/Utils/strangenessBuilderModule.h index 27948b692e1..41d995fbb44 100644 --- a/PWGLF/Utils/strangenessBuilderModule.h +++ b/PWGLF/Utils/strangenessBuilderModule.h @@ -558,12 +558,11 @@ class BuilderModule } } - if(nEnabledTables==0){ + if (nEnabledTables == 0) { LOGF(info, "Strangeness building not required. Will suppress all functionality, including logs, from this point forward."); return; } - // setup bookkeeping histogram auto h = histos.template add("hTableBuildingStatistics", "hTableBuildingStatistics", o2::framework::kTH1D, {{nTablesConst, -0.5f, static_cast(nTablesConst)}}); auto h2 = histos.template add("hInputStatistics", "hInputStatistics", o2::framework::kTH1D, {{nTablesConst, -0.5f, static_cast(nTablesConst)}}); @@ -637,7 +636,7 @@ class BuilderModule // adjust bookkeeping histogram h->GetXaxis()->SetBinLabel(i + 1, tableNames[i].c_str()); h2->GetXaxis()->SetBinLabel(i + 1, tableNames[i].c_str()); - if(baseOpts.mEnabledTables[i] == false){ + if (baseOpts.mEnabledTables[i] == false) { h->SetBinContent(i + 1, -1); // mark disabled tables, distinguish from zero counts } } @@ -2522,7 +2521,7 @@ class BuilderModule template void dataProcess(TCCDB& ccdb, THistoRegistry& histos, TCollisions const& collisions, TMCCollisions const& mccollisions, TV0s const& v0s, TCascades const& cascades, TTrackedCascades const& trackedCascades, TTracks const& tracks, TBCs const& bcs, TMCParticles const& mcParticles, TProducts& products) { - if(nEnabledTables == 0){ + if (nEnabledTables == 0) { return; // fully suppressed } From 3ccb42dda5c0a0701b0baafe18bcd4a8e257f459 Mon Sep 17 00:00:00 2001 From: ddobrigk Date: Mon, 21 Jul 2025 12:25:45 +0200 Subject: [PATCH 09/12] Update propagationService.cxx --- PWGLF/TableProducer/Strangeness/propagationService.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGLF/TableProducer/Strangeness/propagationService.cxx b/PWGLF/TableProducer/Strangeness/propagationService.cxx index dd7353a594d..fa88d6fb47a 100644 --- a/PWGLF/TableProducer/Strangeness/propagationService.cxx +++ b/PWGLF/TableProducer/Strangeness/propagationService.cxx @@ -21,7 +21,7 @@ // strangeness building in a single DPL device that is particularly // adequate for pipelining. // -// Currently meant for testing and performance evaluation +// Currently meant for testing and performance checks // //=============================================================== From dd83a33431bdf58747fcc52dee7ab55dd13e788b Mon Sep 17 00:00:00 2001 From: David Dobrigkeit Chinellato Date: Tue, 22 Jul 2025 04:51:34 -0300 Subject: [PATCH 10/12] Comment out filtering part --- .../Strangeness/propagationService.cxx | 33 ++++++++++--------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/propagationService.cxx b/PWGLF/TableProducer/Strangeness/propagationService.cxx index fa88d6fb47a..c470f163533 100644 --- a/PWGLF/TableProducer/Strangeness/propagationService.cxx +++ b/PWGLF/TableProducer/Strangeness/propagationService.cxx @@ -115,24 +115,27 @@ struct propagationService { strangenessBuilderModule.dataProcess(ccdb, histos, collisions, mccollisions, v0s, cascades, trackedCascades, tracks, bcs, mcParticles, products); } - void processRealDataWithPID(soa::Join const& collisions, aod::V0s const& v0s, aod::Cascades const& cascades, aod::TrackedCascades const& trackedCascades, FullTracksExtIUWithPID const& tracks, aod::BCsWithTimestamps const& bcs) - { - ccdbLoader.initCCDBfromBCs(standardCCDBLoaderConfigurables, ccdb, bcs); - trackPropagation.fillTrackTables(trackPropagationConfigurables, ccdbLoader, collisions, tracks, trackPropagationProducts, histos); - strangenessBuilderModule.dataProcess(ccdb, histos, collisions, static_cast(nullptr), v0s, cascades, trackedCascades, tracks, bcs, static_cast(nullptr), products); - } - - void processMonteCarloWithPID(soa::Join const& collisions, aod::McCollisions const& mccollisions, aod::V0s const& v0s, aod::Cascades const& cascades, aod::TrackedCascades const& trackedCascades, FullTracksExtLabeledIUWithPID const& tracks, aod::BCsWithTimestamps const& bcs, aod::McParticles const& mcParticles) - { - ccdbLoader.initCCDBfromBCs(standardCCDBLoaderConfigurables, ccdb, bcs); - trackPropagation.fillTrackTables(trackPropagationConfigurables, ccdbLoader, collisions, tracks, trackPropagationProducts, histos); - strangenessBuilderModule.dataProcess(ccdb, histos, collisions, mccollisions, v0s, cascades, trackedCascades, tracks, bcs, mcParticles, products); - } + // FIXME: the part below is only viable if TPC PID + // switches to using TracksIU (circular dependency) + // + // void processRealDataWithPID(soa::Join const& collisions, aod::V0s const& v0s, aod::Cascades const& cascades, aod::TrackedCascades const& trackedCascades, FullTracksExtIUWithPID const& tracks, aod::BCsWithTimestamps const& bcs) + // { + // ccdbLoader.initCCDBfromBCs(standardCCDBLoaderConfigurables, ccdb, bcs); + // trackPropagation.fillTrackTables(trackPropagationConfigurables, ccdbLoader, collisions, tracks, trackPropagationProducts, histos); + // strangenessBuilderModule.dataProcess(ccdb, histos, collisions, static_cast(nullptr), v0s, cascades, trackedCascades, tracks, bcs, static_cast(nullptr), products); + // } + + // void processMonteCarloWithPID(soa::Join const& collisions, aod::McCollisions const& mccollisions, aod::V0s const& v0s, aod::Cascades const& cascades, aod::TrackedCascades const& trackedCascades, FullTracksExtLabeledIUWithPID const& tracks, aod::BCsWithTimestamps const& bcs, aod::McParticles const& mcParticles) + // { + // ccdbLoader.initCCDBfromBCs(standardCCDBLoaderConfigurables, ccdb, bcs); + // trackPropagation.fillTrackTables(trackPropagationConfigurables, ccdbLoader, collisions, tracks, trackPropagationProducts, histos); + // strangenessBuilderModule.dataProcess(ccdb, histos, collisions, mccollisions, v0s, cascades, trackedCascades, tracks, bcs, mcParticles, products); + // } PROCESS_SWITCH(propagationService, processRealData, "process real data", true); PROCESS_SWITCH(propagationService, processMonteCarlo, "process monte carlo", false); - PROCESS_SWITCH(propagationService, processRealDataWithPID, "process real data", false); - PROCESS_SWITCH(propagationService, processMonteCarloWithPID, "process monte carlo", false); + // PROCESS_SWITCH(propagationService, processRealDataWithPID, "process real data", false); + // PROCESS_SWITCH(propagationService, processMonteCarloWithPID, "process monte carlo", false); }; //**************************************************************************************** From 1310637ed222c34ab4a6345ebda2d905d2375222 Mon Sep 17 00:00:00 2001 From: ALICE Builder Date: Tue, 22 Jul 2025 10:03:25 +0200 Subject: [PATCH 11/12] Please consider the following formatting changes (#445) --- PWGLF/TableProducer/Strangeness/propagationService.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGLF/TableProducer/Strangeness/propagationService.cxx b/PWGLF/TableProducer/Strangeness/propagationService.cxx index c470f163533..3febe10fc6e 100644 --- a/PWGLF/TableProducer/Strangeness/propagationService.cxx +++ b/PWGLF/TableProducer/Strangeness/propagationService.cxx @@ -115,9 +115,9 @@ struct propagationService { strangenessBuilderModule.dataProcess(ccdb, histos, collisions, mccollisions, v0s, cascades, trackedCascades, tracks, bcs, mcParticles, products); } - // FIXME: the part below is only viable if TPC PID + // FIXME: the part below is only viable if TPC PID // switches to using TracksIU (circular dependency) - // + // // void processRealDataWithPID(soa::Join const& collisions, aod::V0s const& v0s, aod::Cascades const& cascades, aod::TrackedCascades const& trackedCascades, FullTracksExtIUWithPID const& tracks, aod::BCsWithTimestamps const& bcs) // { // ccdbLoader.initCCDBfromBCs(standardCCDBLoaderConfigurables, ccdb, bcs); From aa751187675cf84abfcb732ef70e23f8ddcdc9bf Mon Sep 17 00:00:00 2001 From: ddobrigk Date: Tue, 22 Jul 2025 18:46:20 +0200 Subject: [PATCH 12/12] Update propagationService.cxx --- PWGLF/TableProducer/Strangeness/propagationService.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGLF/TableProducer/Strangeness/propagationService.cxx b/PWGLF/TableProducer/Strangeness/propagationService.cxx index 3febe10fc6e..6140dd2cc8b 100644 --- a/PWGLF/TableProducer/Strangeness/propagationService.cxx +++ b/PWGLF/TableProducer/Strangeness/propagationService.cxx @@ -21,7 +21,7 @@ // strangeness building in a single DPL device that is particularly // adequate for pipelining. // -// Currently meant for testing and performance checks +// Currently meant for testing and performance check // //===============================================================