diff --git a/ALICE3/DataModel/A3DecayFinderTables.h b/ALICE3/DataModel/A3DecayFinderTables.h index cc7e8f3d3c7..05eb5ad9307 100644 --- a/ALICE3/DataModel/A3DecayFinderTables.h +++ b/ALICE3/DataModel/A3DecayFinderTables.h @@ -43,7 +43,10 @@ enum a3selectionBit : uint32_t { kDCAxy = 0, kTruePrPlusFromLc, kTruePiMinusFromLc, kTrueKaMinusFromLc, - kTruePrMinusFromLc }; + kTruePrMinusFromLc, + kTrueXiFromXiC, + kTruePiFromXiC, + kTruePiFromXiCC }; namespace o2::aod { @@ -51,7 +54,7 @@ namespace a3DecayMap { DECLARE_SOA_COLUMN(DecayMap, decayMap, uint32_t); //! simple map to process passing / not passing criteria } // namespace a3DecayMap -DECLARE_SOA_TABLE(Alice3DecayMaps, "AOD", "ALICE3DECAYMAP", +DECLARE_SOA_TABLE(Alice3DecayMaps, "AOD", "ALICE3DECAYMAPS", a3DecayMap::DecayMap); using Alice3DecayMap = Alice3DecayMaps::iterator; diff --git a/ALICE3/DataModel/OTFMulticharm.h b/ALICE3/DataModel/OTFMulticharm.h new file mode 100644 index 00000000000..021288c0426 --- /dev/null +++ b/ALICE3/DataModel/OTFMulticharm.h @@ -0,0 +1,57 @@ +// 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 OTFStrangeness.h +/// \author David Dobrigkeit Chinellato +/// \since 05/08/2024 +/// \brief Set of tables for the ALICE3 strangeness information +/// + +#ifndef ALICE3_DATAMODEL_OTFMULTICHARM_H_ +#define ALICE3_DATAMODEL_OTFMULTICHARM_H_ + +// O2 includes +#include "Framework/AnalysisDataModel.h" + +namespace o2::aod +{ +namespace otfmulticharm +{ +DECLARE_SOA_INDEX_COLUMN_FULL(Cascade, cascade, int, UpgradeCascades, "_Cascade"); +DECLARE_SOA_INDEX_COLUMN_FULL(XiCPion1, xiCPion1, int, Tracks, "_Pi1XiC"); +DECLARE_SOA_INDEX_COLUMN_FULL(XiCPion2, xiCPion2, int, Tracks, "_Pi2XiC"); +DECLARE_SOA_INDEX_COLUMN_FULL(XiCCPion, xiCCPion, int, Tracks, "_PiXiCC"); + +// topo vars +DECLARE_SOA_COLUMN(DCAXiCDaughters, dcaXiCDaughters, float); +DECLARE_SOA_COLUMN(DCAXiCCDaughters, dcaXiCCDaughters, float); + +DECLARE_SOA_COLUMN(MXiC, mXiC, float); +DECLARE_SOA_COLUMN(MXiCC, mXiCC, float); + +} // namespace otfmulticharm +DECLARE_SOA_TABLE(MultiCharmStates, "AOD", "MultiCharmStates", + o2::soa::Index<>, + otfcascade::CascadeId, + otfcascade::XiCPion1Id, + otfcascade::XiCPion2Id, + otfcascade::XiCCPionId, + otfcascade::DCAXiCDaughters, + otfcascade::DCAXiCCDaughters, + otfcascade::MXiC, + otfcascade::MXiCC); + +using MultiCharmState = MultiCharmState::iterator; + +} // namespace o2::aod + +#endif // ALICE3_DATAMODEL_OTFMULTICHARM_H_ diff --git a/ALICE3/DataModel/OTFStrangeness.h b/ALICE3/DataModel/OTFStrangeness.h new file mode 100644 index 00000000000..861ec3a7af8 --- /dev/null +++ b/ALICE3/DataModel/OTFStrangeness.h @@ -0,0 +1,70 @@ +// 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 OTFStrangeness.h +/// \author David Dobrigkeit Chinellato +/// \since 05/08/2024 +/// \brief Set of tables for the ALICE3 strangeness information +/// + +#ifndef ALICE3_DATAMODEL_OTFSTRANGENESS_H_ +#define ALICE3_DATAMODEL_OTFSTRANGENESS_H_ + +// O2 includes +#include "Framework/AnalysisDataModel.h" + +namespace o2::aod +{ +namespace otfcascade +{ +DECLARE_SOA_INDEX_COLUMN(Collision, collision); //! +DECLARE_SOA_INDEX_COLUMN_FULL(CascadeTrack, cascadeTrack, int, Tracks, "_Cascade"); //! +DECLARE_SOA_INDEX_COLUMN_FULL(PosTrack, posTrack, int, Tracks, "_Pos"); //! +DECLARE_SOA_INDEX_COLUMN_FULL(NegTrack, negTrack, int, Tracks, "_Neg"); //! +DECLARE_SOA_INDEX_COLUMN_FULL(BachTrack, bachTrack, int, Tracks, "_Bach"); //! + +// topo vars +DECLARE_SOA_COLUMN(DCAV0Daughters, dcaV0Daughters, float); +DECLARE_SOA_COLUMN(DCACascadeDaughters, dcaCascadeDaughters, float); +DECLARE_SOA_COLUMN(V0Radius, v0Radius, float); +DECLARE_SOA_COLUMN(CascRadius, cascRadius, float); +DECLARE_SOA_COLUMN(CascRadiusMC, cascRadiusMC, float); +DECLARE_SOA_COLUMN(MLambda, mLambda, float); +DECLARE_SOA_COLUMN(MXi, mXi, float); + +// strangeness tracking +DECLARE_SOA_COLUMN(FindableClusters, findableClusters, int); +DECLARE_SOA_COLUMN(FoundClusters, foundClusters, int); + +} // namespace otfcascade +DECLARE_SOA_TABLE(UpgradeCascades, "AOD", "UPGRADECASCADES", + o2::soa::Index<>, + otfcascade::CollisionId, + otfcascade::CascadeTrackId, + otfcascade::PosTrackId, + otfcascade::NegTrackId, + otfcascade::BachTrackId, + otfcascade::DCAV0Daughters, + otfcascade::DCACascadeDaughters, + otfcascade::V0Radius, + otfcascade::CascRadius, + otfcascade::CascRadiusMC, + otfcascade::MLambda, + otfcascade::MXi, + otfcascade::FindableClusters, + otfcascade::FoundClusters); + +using UpgradeCascade = UpgradeCascades::iterator; + +} // namespace o2::aod + +#endif // ALICE3_DATAMODEL_OTFSTRANGENESS_H_ diff --git a/ALICE3/TableProducer/CMakeLists.txt b/ALICE3/TableProducer/CMakeLists.txt index 52300bfc644..9c3d710d358 100644 --- a/ALICE3/TableProducer/CMakeLists.txt +++ b/ALICE3/TableProducer/CMakeLists.txt @@ -31,8 +31,17 @@ o2physics_add_dpl_workflow(alice3-centrality PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(alice3-decaypreselector + SOURCES alice3-decaypreselector.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(alice3-decayfinder SOURCES alice3-decayfinder.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(alice3-multicharm + SOURCES alice3-multicharm.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter + COMPONENT_NAME Analysis) diff --git a/ALICE3/TableProducer/OTF/CMakeLists.txt b/ALICE3/TableProducer/OTF/CMakeLists.txt index 305fcce4b66..55427691f4f 100644 --- a/ALICE3/TableProducer/OTF/CMakeLists.txt +++ b/ALICE3/TableProducer/OTF/CMakeLists.txt @@ -11,7 +11,7 @@ o2physics_add_dpl_workflow(onthefly-tracker SOURCES onTheFlyTracker.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2::DetectorsVertexing O2Physics::ALICE3Core + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2::DetectorsVertexing O2::DCAFitter O2Physics::ALICE3Core COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(onthefly-tofpid diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index a73fabdbc80..1b5682daad4 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -24,6 +24,7 @@ /// #include +#include #include #include @@ -35,6 +36,8 @@ #include "Framework/runDataProcessing.h" #include "Framework/HistogramRegistry.h" #include +#include "DCAFitter/DCAFitterN.h" +#include "Common/Core/RecoDecay.h" #include "Framework/O2DatabasePDGPlugin.h" #include "Common/DataModel/TrackSelectionTables.h" #include "ReconstructionDataFormats/DCA.h" @@ -45,13 +48,21 @@ #include "SimulationDataFormat/InteractionSampler.h" #include "Field/MagneticField.h" +#include "ITSMFTSimulation/Hit.h" +#include "ITStracking/Configuration.h" +#include "ITStracking/IOUtils.h" +#include "ITStracking/Tracker.h" +#include "ITStracking/Vertexer.h" +#include "ITStracking/VertexerTraits.h" + #include "ALICE3/Core/DelphesO2TrackSmearer.h" #include "ALICE3/DataModel/collisionAlice3.h" #include "ALICE3/DataModel/tracksAlice3.h" -#include "ALICE3/DataModel/OTFMcTrackExtra.h" +#include "ALICE3/DataModel/OTFStrangeness.h" using namespace o2; using namespace o2::framework; +using std::array; struct OnTheFlyTracker { Produces collisions; @@ -65,7 +76,7 @@ struct OnTheFlyTracker { Produces tracksDCACov; Produces collisionsAlice3; Produces TracksAlice3; - Produces OTFMcExtra; + Produces upgradeCascades; // optionally produced, empty (to be tuned later) Produces tracksExtra; // base table, extend later @@ -81,7 +92,10 @@ struct OnTheFlyTracker { Configurable enableNucleiSmearing{"enableNucleiSmearing", false, "Enable smearing of nuclei"}; Configurable enablePrimaryVertexing{"enablePrimaryVertexing", true, "Enable primary vertexing"}; Configurable interpolateLutEfficiencyVsNch{"interpolateLutEfficiencyVsNch", true, "interpolate LUT efficiency as f(Nch)"}; - Configurable treatXi{"treatXi", false, "Manually decay Xi^{-} and fill tables with daughters"}; + Configurable treatXi{"treatXi", false, "Manually decay Xi and fill tables with daughters"}; + Configurable findXi{"findXi", false, "if treatXi on, find Xi and fill Tracks table also with Xi"}; + Configurable trackXi{"trackXi", false, "if findXi on, attempt to track Xi"}; + Configurable pixelResolution{"pixelResolution", 2.5e-4, "pixel resolution in centimeters (for OTF strangeness tracking)"}; Configurable populateTracksDCA{"populateTracksDCA", true, "populate TracksDCA table"}; Configurable populateTracksDCACov{"populateTracksDCACov", false, "populate TracksDCACov table"}; @@ -91,6 +105,7 @@ struct OnTheFlyTracker { Configurable processUnreconstructedTracks{"processUnreconstructedTracks", false, "process (smear) unreco-ed tracks"}; Configurable doExtraQA{"doExtraQA", false, "do extra 2D QA plots"}; Configurable extraQAwithoutDecayDaughters{"extraQAwithoutDecayDaughters", false, "remove decay daughters from qa plots (yes/no)"}; + Configurable doXiQA{"doXiQA", false, "QA plots for when treating Xi"}; Configurable lutEl{"lutEl", "lutCovm.el.dat", "LUT for electrons"}; Configurable lutMu{"lutMu", "lutCovm.mu.dat", "LUT for muons"}; @@ -121,9 +136,14 @@ struct OnTheFlyTracker { ConfigurableAxis axisDCA{"axisDCA", {400, -200, 200}, "DCA (#mum)"}; ConfigurableAxis axisX{"axisX", {250, -50, 200}, "track X (cm)"}; ConfigurableAxis axisRadius{"axisRadius", {55, 0.01, 100}, "decay radius"}; + ConfigurableAxis axisLambdaMass{"axisLambdaMass", {200, 1.101f, 1.131f}, ""}; + ConfigurableAxis axisXiMass{"axisXiMass", {200, 1.22f, 1.42f}, ""}; using PVertex = o2::dataformats::PrimaryVertex; + // for secondary vertex finding + o2::vertexing::DCAFitterN<2> fitter; + // Class to hold the track information for the O2 vertexing class TrackAlice3 : public o2::track::TrackParCov { @@ -133,12 +153,41 @@ struct OnTheFlyTracker { TrackAlice3() = default; ~TrackAlice3() = default; TrackAlice3(const TrackAlice3& src) = default; - TrackAlice3(const o2::track::TrackParCov& src, const int64_t label, const float t = 0, const float te = 1, bool decayDauInput = false) : o2::track::TrackParCov(src), mcLabel{label}, timeEst{t, te}, isDecayDau(decayDauInput) {} + TrackAlice3(const o2::track::TrackParCov& src, const int64_t label, const float t = 0, const float te = 1, bool decayDauInput = false, bool weakDecayDauInput = false, int isUsedInCascadingInput = 0) : o2::track::TrackParCov(src), + mcLabel{label}, + timeEst{t, te}, + isDecayDau(decayDauInput), + isWeakDecayDau(weakDecayDauInput), + isUsedInCascading(isUsedInCascadingInput) {} const TimeEst& getTimeMUS() const { return timeEst; } int64_t mcLabel; TimeEst timeEst; ///< time estimate in ns bool isDecayDau; + bool isWeakDecayDau; + int isUsedInCascading; // 0: not at all, 1: is a cascade, 2: is a bachelor, 3: is a pion, 4: is a proton + }; + + // Helper struct to pass cascade information + struct cascadecandidate { + int cascadeTrackId; // track index in the Tracks table + int positiveId; // track index in the Tracks table + int negativeId; // track index in the Tracks table + int bachelorId; // track index in the Tracks table + + float dcaV0dau; + float dcacascdau; + float v0radius; + float cascradius; + float cascradiusMC; + + // for tracking + int findableClusters; + int foundClusters; + + float mLambda; + float mXi; }; + cascadecandidate thisCascade; // necessary for particle charges Service pdgDB; @@ -162,15 +211,10 @@ struct OnTheFlyTracker { // For processing and vertexing std::vector tracksAlice3; std::vector ghostTracksAlice3; - std::vector trackPdg; - std::vector trackIsFromXi; - std::vector trackIsFromL0; - std::vector ghostTrackPdg; - std::vector ghostTrackIsFromXi; - std::vector ghostTrackIsFromL0; std::vector bcData; o2::steer::InteractionSampler irSampler; o2::vertexing::PVertexer vertexer; + std::vector cascadesAlice3; void init(o2::framework::InitContext&) { @@ -339,7 +383,7 @@ struct OnTheFlyTracker { histos.add("hTrackXatDCA", "hTrackXatDCA", kTH1F, {axisX}); } - if (treatXi) { + if (doXiQA) { histos.add("hGenXi", "hGenXi", kTH2F, {axisRadius, axisMomentum}); histos.add("hRecoXi", "hRecoXi", kTH2F, {axisRadius, axisMomentum}); @@ -350,12 +394,17 @@ struct OnTheFlyTracker { histos.add("hRecoPiFromL0", "hRecoPiFromL0", kTH2F, {axisRadius, axisMomentum}); histos.add("hRecoPrFromL0", "hRecoPrFromL0", kTH2F, {axisRadius, axisMomentum}); - histos.add("hPiFromXiDCAxy", "hPiFromXiDCAxy", kTH1F, {axisDCA}); - histos.add("hPiFromL0DCAxy", "hPiFromL0DCAxy", kTH1F, {axisDCA}); - histos.add("hPrFromL0DCAxy", "hPrFromL0DCAxy", kTH1F, {axisDCA}); - histos.add("hPiFromXiDCAxyVsPt", "hPiFromXiDCAxyVsPt", kTH2F, {axisMomentum, axisDCA}); - histos.add("hPiFromL0DCAxyVsPt", "hPiFromL0DCAxyVsPt", kTH2F, {axisMomentum, axisDCA}); - histos.add("hPrFromL0DCAxyVsPt", "hPrFromL0DCAxyVsPt", kTH2F, {axisMomentum, axisDCA}); + // basic mass histograms to see if we're in business + histos.add("hMassLambda", "hMassLambda", kTH1F, {axisLambdaMass}); + histos.add("hMassXi", "hMassXi", kTH1F, {axisXiMass}); + + // OTF strangeness tracking QA + histos.add("hFoundVsFindable", "hFoundVsFindable", kTH2F, {{10, -0.5f, 9.5f}, {10, -0.5f, 9.5f}}); + + histos.add("h2dDCAxyCascade", "h2dDCAxyCascade", kTH2F, {axisMomentum, axisDCA}); + histos.add("h2dDCAxyCascadeBachelor", "h2dDCAxyCascadeBachelor", kTH2F, {axisMomentum, axisDCA}); + histos.add("h2dDCAxyCascadeNegative", "h2dDCAxyCascadeNegative", kTH2F, {axisMomentum, axisDCA}); + histos.add("h2dDCAxyCascadePositive", "h2dDCAxyCascadePositive", kTH2F, {axisMomentum, axisDCA}); } LOGF(info, "Initializing magnetic field to value: %.3f kG", static_cast(magneticField)); @@ -388,6 +437,18 @@ struct OnTheFlyTracker { vertexer.setValidateWithIR(kFALSE); vertexer.setBunchFilling(irSampler.getBunchFilling()); vertexer.init(); + + // initialize O2 2-prong fitter + fitter.setPropagateToPCA(true); + fitter.setMaxR(200.); + fitter.setMinParamChange(1e-3); + fitter.setMinRelChi2Change(0.9); + fitter.setMaxDZIni(1e9); + fitter.setMaxDXYIni(4); + fitter.setMaxChi2(1e9); + fitter.setUseAbsDCA(true); + fitter.setWeightedFinalPCA(false); + fitter.setMatCorrType(o2::base::Propagator::MatCorrType::USEMatCorrNONE); // such a light detector here } /// Function to decay the xi @@ -509,13 +570,8 @@ struct OnTheFlyTracker { { tracksAlice3.clear(); ghostTracksAlice3.clear(); - trackPdg.clear(); - trackIsFromXi.clear(); - trackIsFromL0.clear(); - ghostTrackPdg.clear(); - ghostTrackIsFromXi.clear(); - ghostTrackIsFromL0.clear(); bcData.clear(); + cascadesAlice3.clear(); o2::dataformats::DCA dcaInfo; o2::dataformats::VertexBase vtx; @@ -599,7 +655,7 @@ struct OnTheFlyTracker { if (TMath::Abs(mcParticle.pdgCode()) == 2212) histos.fill(HIST("hPtGeneratedPr"), mcParticle.pt()); - if (treatXi && mcParticle.pdgCode() == 3312) { + if (doXiQA && mcParticle.pdgCode() == 3312) { histos.fill(HIST("hGenXi"), xiDecayRadius2D, mcParticle.pt()); histos.fill(HIST("hGenPiFromXi"), xiDecayRadius2D, decayProducts[0].Pt()); histos.fill(HIST("hGenPiFromL0"), l0DecayRadius2D, decayProducts[1].Pt()); @@ -610,6 +666,9 @@ struct OnTheFlyTracker { continue; } + o2::track::TrackParCov trackParCov; + convertMCParticleToO2Track(mcParticle, trackParCov); + bool isDecayDaughter = false; if (mcParticle.getProcess() == 4) isDecayDaughter = true; @@ -618,9 +677,6 @@ struct OnTheFlyTracker { const float t = (ir.timeInBCNS + gRandom->Gaus(0., 100.)) * 1e-3; std::vector xiDaughterTrackParCovs(3); std::vector isReco(3); - std::vector dauPdg = {-211, -211, 2212}; - std::vector fromXi = {true, false, false}; - std::vector fromL0 = {false, true, true}; std::vector smearer = {mSmearer0, mSmearer1, mSmearer2, mSmearer3, mSmearer4, mSmearer5}; if (treatXi && mcParticle.pdgCode() == 3312) { if (xiDecayRadius2D > 20) { @@ -670,19 +726,13 @@ struct OnTheFlyTracker { continue; } if (isReco[i]) { - tracksAlice3.push_back(TrackAlice3{xiDaughterTrackParCovs[i], mcParticle.globalIndex(), t, 100.f * 1e-3, true}); - trackPdg.push_back(dauPdg[i]); - trackIsFromXi.push_back(fromXi[i]); - trackIsFromL0.push_back(fromL0[i]); + tracksAlice3.push_back(TrackAlice3{xiDaughterTrackParCovs[i], mcParticle.globalIndex(), t, 100.f * 1e-3, true, true, i + 2}); } else { - ghostTracksAlice3.push_back(TrackAlice3{xiDaughterTrackParCovs[i], mcParticle.globalIndex(), t, 100.f * 1e-3, true}); - ghostTrackPdg.push_back(dauPdg[i]); - ghostTrackIsFromXi.push_back(fromXi[i]); - ghostTrackIsFromL0.push_back(fromL0[i]); + ghostTracksAlice3.push_back(TrackAlice3{xiDaughterTrackParCovs[i], mcParticle.globalIndex(), t, 100.f * 1e-3, true, true, i + 2}); } } - if (treatXi && mcParticle.pdgCode() == 3312) { + if (doXiQA && mcParticle.pdgCode() == 3312) { if (isReco[0] && isReco[1] && isReco[2]) histos.fill(HIST("hRecoXi"), xiDecayRadius2D, mcParticle.pt()); if (isReco[0]) @@ -692,12 +742,179 @@ struct OnTheFlyTracker { if (isReco[2]) histos.fill(HIST("hRecoPrFromL0"), l0DecayRadius2D, decayProducts[2].Pt()); } + + // +-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+ + // combine particles into actual Xi candidate + // cascade building starts here + if (findXi && mcParticle.pdgCode() == 3312 && isReco[0] && isReco[1] && isReco[2]) { + // assign indices of the particles we've used + // they should be the last ones to be filled, in order: + // n-1: proton from lambda + // n-2: pion from lambda + // n-3: pion from xi + thisCascade.positiveId = tracksAlice3.size() - 1; + thisCascade.negativeId = tracksAlice3.size() - 2; + thisCascade.bachelorId = tracksAlice3.size() - 3; + + // use DCA fitters + int nCand = 0; + bool dcaFitterOK_V0 = true; + try { + nCand = fitter.process(xiDaughterTrackParCovs[1], xiDaughterTrackParCovs[2]); + } catch (...) { + // LOG(error) << "Exception caught in DCA fitter process call!"; + dcaFitterOK_V0 = false; + } + if (nCand == 0) { + dcaFitterOK_V0 = false; + } + // V0 found successfully + if (dcaFitterOK_V0) { + std::array pos; + std::array posCascade; + std::array posP; + std::array negP; + std::array bachP; + + o2::track::TrackParCov pTrackAtPCA = fitter.getTrack(1); // proton (positive) + o2::track::TrackParCov nTrackAtPCA = fitter.getTrack(0); // pion (negative) + pTrackAtPCA.getPxPyPzGlo(posP); + nTrackAtPCA.getPxPyPzGlo(negP); + + // get decay vertex coordinates + const auto& vtx = fitter.getPCACandidate(); + for (int i = 0; i < 3; i++) { + pos[i] = vtx[i]; + } + + // calculate basic V0 properties here + // DCA to PV taken care of in daughter tracks already, not necessary + thisCascade.dcaV0dau = TMath::Sqrt(fitter.getChi2AtPCACandidate()); + thisCascade.v0radius = std::hypot(pos[0], pos[1]); + thisCascade.mLambda = RecoDecay::m(array{array{posP[0], posP[1], posP[2]}, array{negP[0], negP[1], negP[2]}}, array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged}); + + // go for cascade: create V0 (pseudo)track from reconstructed V0 + std::array covV = {0.}; + constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component + for (int i = 0; i < 6; i++) { + covV[MomInd[i]] = 1e-6; + covV[i] = 1e-6; + } + o2::track::TrackParCov v0Track = o2::track::TrackParCov( + {pos[0], pos[1], pos[2]}, + {posP[0] + negP[0], posP[1] + negP[1], posP[2] + negP[2]}, + covV, 0, true); + v0Track.setAbsCharge(0); + v0Track.setPID(o2::track::PID::Lambda); + + // dca fitter step + nCand = 0; + bool dcaFitterOK_Cascade = true; + try { + nCand = fitter.process(v0Track, xiDaughterTrackParCovs[0]); + } catch (...) { + // LOG(error) << "Exception caught in DCA fitter process call!"; + dcaFitterOK_Cascade = false; + } + if (nCand == 0) { + dcaFitterOK_Cascade = false; + } + + // Cascade found successfully + if (dcaFitterOK_Cascade) { + o2::track::TrackParCov bachelorTrackAtPCA = fitter.getTrack(1); + + const auto& vtxCascade = fitter.getPCACandidate(); + for (int i = 0; i < 3; i++) { + posCascade[i] = vtxCascade[i]; + } + + // basic properties of the cascade + thisCascade.dcacascdau = TMath::Sqrt(fitter.getChi2AtPCACandidate()); + thisCascade.cascradius = std::hypot(posCascade[0], posCascade[1]); + bachelorTrackAtPCA.getPxPyPzGlo(bachP); + + thisCascade.mXi = RecoDecay::m(array{array{bachP[0], bachP[1], bachP[2]}, array{posP[0] + negP[0], posP[1] + negP[1], posP[2] + negP[2]}}, array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassLambda}); + + // initialize cascade track + o2::track::TrackParCov cascadeTrack = fitter.createParentTrackParCov(); + cascadeTrack.setAbsCharge(-1); // may require more adjustments + cascadeTrack.setPID(o2::track::PID::XiMinus); // FIXME: not OK for omegas + + thisCascade.cascradiusMC = xiDecayRadius2D; + thisCascade.findableClusters = 0; + thisCascade.foundClusters = 0; + + if (trackXi) { + // optionally, add the points in the layers before the decay of the Xi + // will back-track the perfect MC cascade to relevant layers, find hit, smear and add to smeared cascade + for (int i = layers.size() - 1; i >= 0; i--) { + if (thisCascade.cascradiusMC > layers[i]) { + // will add this layer, since cascade decayed after the corresponding radius + thisCascade.findableClusters++; // add to findable + + // find perfect intercept XYZ + float targetX = 1e+3; + trackParCov.getXatLabR(layers[i], targetX, magneticField); + if (targetX > 999) + continue; // failed to find intercept + + if (!trackParCov.propagateTo(targetX, magneticField)) { + continue; // failed to propagate + } + + // get potential cluster position + std::array posClusterCandidate; + trackParCov.getXYZGlo(posClusterCandidate); + float r{std::hypot(posClusterCandidate[0], posClusterCandidate[1])}; + float phi{std::atan2(-posClusterCandidate[1], -posClusterCandidate[0]) + o2::its::constants::math::Pi}; + + if (pixelResolution > 1e-8) { // catch zero (though should not really happen...) + phi = gRandom->Gaus(phi, std::asin(pixelResolution / r)); + posClusterCandidate[0] = r * std::cos(phi); + posClusterCandidate[1] = r * std::sin(phi); + posClusterCandidate[2] = gRandom->Gaus(posClusterCandidate[2], pixelResolution); + } + + // towards adding cluster: move to track alpha + double alpha = cascadeTrack.getAlpha(); + double xyz1[3]{ + TMath::Cos(alpha) * posClusterCandidate[0] + TMath::Sin(alpha) * posClusterCandidate[1], + -TMath::Sin(alpha) * posClusterCandidate[0] + TMath::Cos(alpha) * posClusterCandidate[1], + posClusterCandidate[2]}; + + if (!(cascadeTrack.propagateTo(xyz1[0], magneticField))) + continue; + const o2::track::TrackParametrization::dim2_t hitpoint = { + static_cast(xyz1[1]), + static_cast(xyz1[2])}; + const o2::track::TrackParametrization::dim3_t hitpointcov = {pixelResolution * pixelResolution, 0.f, pixelResolution * pixelResolution}; + cascadeTrack.update(hitpoint, hitpointcov); + thisCascade.foundClusters++; // add to findable + } + } + } + + // add cascade track + thisCascade.cascadeTrackId = tracksAlice3.size(); // this is the next index to be filled -> should be it + tracksAlice3.push_back(TrackAlice3{cascadeTrack, mcParticle.globalIndex(), t, 100.f * 1e-3, false, false, 1}); + + if (doXiQA) { + histos.fill(HIST("hMassLambda"), thisCascade.mLambda); + histos.fill(HIST("hMassXi"), thisCascade.mXi); + histos.fill(HIST("hFoundVsFindable"), thisCascade.findableClusters, thisCascade.foundClusters); + } + + // add this cascade to vector (will fill cursor later with collision ID) + cascadesAlice3.push_back(thisCascade); + } + } + } // end cascade building + // +-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+ + continue; // Not filling the tables with the xi itself } - o2::track::TrackParCov trackParCov; - convertMCParticleToO2Track(mcParticle, trackParCov); - if (doExtraQA) { histos.fill(HIST("hSimTrackX"), trackParCov.getX()); } @@ -729,14 +946,8 @@ struct OnTheFlyTracker { // populate vector with track if we reco-ed it if (reconstructed) { tracksAlice3.push_back(TrackAlice3{trackParCov, mcParticle.globalIndex(), t, 100.f * 1e-3, isDecayDaughter}); - trackPdg.push_back(mcParticle.pdgCode()); - trackIsFromXi.push_back(false); - trackIsFromL0.push_back(false); } else { ghostTracksAlice3.push_back(TrackAlice3{trackParCov, mcParticle.globalIndex(), t, 100.f * 1e-3, isDecayDaughter}); - ghostTrackPdg.push_back(mcParticle.pdgCode()); - ghostTrackIsFromXi.push_back(false); - ghostTrackIsFromL0.push_back(false); } } @@ -813,7 +1024,6 @@ struct OnTheFlyTracker { // *+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+* // populate tracks - unsigned trackCounter = 0; for (const auto& trackParCov : tracksAlice3) { // Fixme: collision index could be changeable aod::track::TrackTypeEnum trackType = aod::track::Track; @@ -829,19 +1039,15 @@ struct OnTheFlyTracker { histos.fill(HIST("h2dDCAxy"), trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please histos.fill(HIST("hTrackXatDCA"), trackParametrization.getX()); } - if (treatXi) { - if (trackPdg[trackCounter] == -211 && trackIsFromXi[trackCounter]) { - histos.fill(HIST("hPiFromXiDCAxy"), dcaXY); - histos.fill(HIST("hPiFromXiDCAxyVsPt"), trackParametrization.getPt(), dcaXY); - } - if (trackPdg[trackCounter] == -211 && trackIsFromL0[trackCounter]) { - histos.fill(HIST("hPiFromL0DCAxy"), dcaXY); - histos.fill(HIST("hPiFromL0DCAxyVsPt"), trackParametrization.getPt(), dcaXY); - } - if (trackPdg[trackCounter] == 2212 && trackIsFromL0[trackCounter]) { - histos.fill(HIST("hPrFromL0DCAxy"), dcaXY); - histos.fill(HIST("hPrFromL0DCAxyVsPt"), trackParametrization.getPt(), dcaXY); - } + if (doXiQA) { + if (trackParCov.isUsedInCascading == 1) + histos.fill(HIST("h2dDCAxyCascade"), trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please + if (trackParCov.isUsedInCascading == 2) + histos.fill(HIST("h2dDCAxyCascadeBachelor"), trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please + if (trackParCov.isUsedInCascading == 3) + histos.fill(HIST("h2dDCAxyCascadeNegative"), trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please + if (trackParCov.isUsedInCascading == 4) + histos.fill(HIST("h2dDCAxyCascadePositive"), trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please } tracksDCA(dcaXY, dcaZ); if (populateTracksDCACov) { @@ -872,14 +1078,8 @@ struct OnTheFlyTracker { trackSelectionExtension(false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false); } TracksAlice3(true); - - // populate with xi daughter info - OTFMcExtra(trackPdg[trackCounter], trackIsFromXi[trackCounter], trackIsFromL0[trackCounter]); - trackCounter++; } - // populate ghost tracks - unsigned ghostTrackCounter = 0; for (const auto& trackParCov : ghostTracksAlice3) { // Fixme: collision index could be changeable aod::track::TrackTypeEnum trackType = aod::track::Track; @@ -925,10 +1125,24 @@ struct OnTheFlyTracker { trackSelectionExtension(false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false); } TracksAlice3(false); + } - // populate table with xi daughter info - OTFMcExtra(ghostTrackPdg[ghostTrackCounter], ghostTrackIsFromXi[trackCounter], ghostTrackIsFromL0[trackCounter]); - ghostTrackCounter++; + for (const auto& cascade : cascadesAlice3) { + upgradeCascades( + collisions.lastIndex(), // now we know the collision index -> populate table + cascade.cascadeTrackId, + cascade.positiveId, + cascade.negativeId, + cascade.bachelorId, + cascade.dcaV0dau, + cascade.dcacascdau, + cascade.v0radius, + cascade.cascradius, + cascade.cascradiusMC, + cascade.mLambda, + cascade.mXi, + cascade.findableClusters, + cascade.foundClusters); } } }; diff --git a/ALICE3/TableProducer/alice3-decayfinder.cxx b/ALICE3/TableProducer/alice3-decayfinder.cxx index 80b9de80238..4b88305fb60 100644 --- a/ALICE3/TableProducer/alice3-decayfinder.cxx +++ b/ALICE3/TableProducer/alice3-decayfinder.cxx @@ -65,149 +65,6 @@ using tofTracks = soa::Join; using richTracks = soa::Join; using alice3tracks = soa::Join; -struct alice3decayPreselector { - Produces a3decayMaps; - - // Operation and minimisation criteria - Configurable nSigmaTOF{"nSigmaTOF", 4.0f, "Nsigma for TOF PID (if enabled)"}; - Configurable nSigmaRICH{"nSigmaRICH", 4.0f, "Nsigma for RICH PID (if enabled)"}; - - // Define o2 fitter, 2-prong, active memory (no need to redefine per event) - o2::vertexing::DCAFitterN<2> fitter; - - // for bit-packed maps - std::vector selectionMap; - - HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - - //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* - /// function to check PDG + PDG mother - template - bool checkPDG(TTrack const& track, int pdgMother, int pdg) - { - bool returnValue = false; - // Association check - if (track.has_mcParticle()) { - auto mcParticle = track.template mcParticle_as(); - if (mcParticle.has_mothers()) { - for (auto& mcParticleMother : mcParticle.template mothers_as()) { - if (mcParticle.pdgCode() == pdg && mcParticleMother.pdgCode() == pdgMother) - returnValue = true; - } - } - } // end association check - return returnValue; - } - //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* - - void init(InitContext&) - { - // future dev if needed - } - - // go declarative: use partitions instead of "if", then just toggle bits to allow for mask selection later - Partition pInnerTOFPi = nabs(aod::upgrade_tof::nSigmaPionInnerTOF) > nSigmaTOF; - Partition pInnerTOFKa = nabs(aod::upgrade_tof::nSigmaKaonInnerTOF) > nSigmaTOF; - Partition pInnerTOFPr = nabs(aod::upgrade_tof::nSigmaProtonInnerTOF) > nSigmaTOF; - Partition pOuterTOFPi = nabs(aod::upgrade_tof::nSigmaPionOuterTOF) > nSigmaTOF; - Partition pOuterTOFKa = nabs(aod::upgrade_tof::nSigmaKaonOuterTOF) > nSigmaTOF; - Partition pOuterTOFPr = nabs(aod::upgrade_tof::nSigmaProtonOuterTOF) > nSigmaTOF; - Partition pRICHPi = nabs(aod::alice3rich::richNsigmaPi) > nSigmaRICH; - Partition pRICHKa = nabs(aod::alice3rich::richNsigmaKa) > nSigmaRICH; - Partition pRICHPr = nabs(aod::alice3rich::richNsigmaPr) > nSigmaRICH; - - //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* - /// Initialization of mask vectors if uninitialized - void initializeMasks(int size) - { - selectionMap.clear(); - selectionMap.resize(size, 0xFFFFFFFF); // all bits 1, please - } - //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* - /// This process function ensures that all V0s are built. It will simply tag everything as true. - void processInitialize(aod::Tracks const& tracks) - { - initializeMasks(tracks.size()); - } - //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* - void processFilterInnerTOF(tofTracks const&) - { - for (auto const& track : pInnerTOFPi) - bitoff(selectionMap[track.globalIndex()], kInnerTOFPion); - for (auto const& track : pInnerTOFKa) - bitoff(selectionMap[track.globalIndex()], kInnerTOFKaon); - for (auto const& track : pInnerTOFPr) - bitoff(selectionMap[track.globalIndex()], kInnerTOFProton); - } - //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* - void processFilterOuterTOF(tofTracks const&) - { - for (auto const& track : pOuterTOFPi) - bitoff(selectionMap[track.globalIndex()], kOuterTOFPion); - for (auto const& track : pOuterTOFKa) - bitoff(selectionMap[track.globalIndex()], kOuterTOFKaon); - for (auto const& track : pOuterTOFPr) - bitoff(selectionMap[track.globalIndex()], kOuterTOFProton); - } - //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* - void processFilterRICH(richTracks const&) - { - for (auto const& track : pRICHPi) - bitoff(selectionMap[track.globalIndex()], kRICHPion); - for (auto const& track : pRICHKa) - bitoff(selectionMap[track.globalIndex()], kRICHKaon); - for (auto const& track : pRICHPr) - bitoff(selectionMap[track.globalIndex()], kRICHProton); - } - //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* - void processFilterOnMonteCarloTruth(labeledTracks const& tracks, aod::McParticles const&) - { - for (auto const& track : tracks) { - // D mesons - if (!checkPDG(track, 421, -321)) //+421 -> -321 +211 - bitoff(selectionMap[track.globalIndex()], kTrueKaMinusFromD); - if (!checkPDG(track, -421, +321)) //-421 -> +321 -211 - bitoff(selectionMap[track.globalIndex()], kTrueKaPlusFromD); - if (!checkPDG(track, 421, +211)) //+421 -> -321 +211 - bitoff(selectionMap[track.globalIndex()], kTruePiPlusFromD); - if (!checkPDG(track, -421, -211)) //-421 -> +321 -211 - bitoff(selectionMap[track.globalIndex()], kTruePiMinusFromD); - - // Lambdac baryons - if (!checkPDG(track, +4122, +2212)) //+4122 -> +2212 -321 +211 - bitoff(selectionMap[track.globalIndex()], kTruePrPlusFromLc); - if (!checkPDG(track, +4122, -321)) //+4122 -> +2212 -321 +211 - bitoff(selectionMap[track.globalIndex()], kTrueKaMinusFromLc); - if (!checkPDG(track, +4122, +211)) //+4122 -> +2212 -321 +211 - bitoff(selectionMap[track.globalIndex()], kTruePiPlusFromLc); - if (!checkPDG(track, -4122, -2212)) //-4122 -> -2212 +321 -211 - bitoff(selectionMap[track.globalIndex()], kTruePrMinusFromLc); - if (!checkPDG(track, -4122, +321)) //-4122 -> -2212 +321 -211 - bitoff(selectionMap[track.globalIndex()], kTrueKaPlusFromLc); - if (!checkPDG(track, -4122, -211)) //-4122 -> -2212 +321 -211 - bitoff(selectionMap[track.globalIndex()], kTruePiMinusFromLc); - } - } - //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* - void processPublishDecision(aod::Tracks const& tracks) - { - for (uint32_t i = 0; i < tracks.size(); i++) { - a3decayMaps(selectionMap[i]); - } - selectionMap.clear(); - } - //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* - - //*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<* - PROCESS_SWITCH(alice3decayPreselector, processInitialize, "Initialize (MUST be on)", true); - PROCESS_SWITCH(alice3decayPreselector, processFilterInnerTOF, "Switch to use inner TOF PID", false); - PROCESS_SWITCH(alice3decayPreselector, processFilterOuterTOF, "Switch to use outer TOF PID", false); - PROCESS_SWITCH(alice3decayPreselector, processFilterRICH, "Switch to use RICH", false); - PROCESS_SWITCH(alice3decayPreselector, processFilterOnMonteCarloTruth, "Switch to use MC truth", false); - PROCESS_SWITCH(alice3decayPreselector, processPublishDecision, "Fill decision mask table (MUST be on)", true); - //*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<* -}; - struct alice3decayFinder { SliceCache cache; @@ -355,9 +212,9 @@ struct alice3decayFinder { } //}-{}-{}-{}-{}-{}-{}-{}-{}-{} - t0 = fitter.getTrack(0); - t1 = fitter.getTrack(1); - t2 = fitter.getTrack(2); + t0 = fitter3.getTrack(0); + t1 = fitter3.getTrack(1); + t2 = fitter3.getTrack(2); std::array P0; std::array P1; std::array P2; @@ -593,6 +450,5 @@ struct alice3decayFinder { WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc)}; } diff --git a/ALICE3/TableProducer/alice3-decaypreselector.cxx b/ALICE3/TableProducer/alice3-decaypreselector.cxx new file mode 100644 index 00000000000..b92b2c996be --- /dev/null +++ b/ALICE3/TableProducer/alice3-decaypreselector.cxx @@ -0,0 +1,222 @@ +// 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. +// +// *+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* +// Decay finder task for ALICE 3 +// *+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* +// +// Uses specific ALICE 3 PID and performance for studying +// HF decays. Work in progress: use at your own risk! +// + +#include +#include +#include +#include +#include +#include + +#include "Framework/runDataProcessing.h" +#include "Framework/RunningWorkflowInfo.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoAHelpers.h" +#include "DCAFitter/DCAFitterN.h" +#include "ReconstructionDataFormats/Track.h" +#include "Common/Core/RecoDecay.h" +#include "Common/Core/trackUtilities.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/DataModel/LFParticleIdentification.h" +#include "Common/Core/TrackSelection.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "DetectorsBase/Propagator.h" +#include "DetectorsBase/GeometryManager.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "CCDB/BasicCCDBManager.h" +#include "DataFormatsCalibration/MeanVertexObject.h" +#include "ALICE3/DataModel/OTFTOF.h" +#include "ALICE3/DataModel/RICH.h" +#include "ALICE3/DataModel/A3DecayFinderTables.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using std::array; + +// simple checkers +// #define biton(var, nbit) ((var) |= (static_cast(1) << (nbit))) +#define bitoff(var, nbit) ((var) &= ~(static_cast(1) << (nbit))) //((a) &= ~(1ULL<<(b))) +// #define bitcheck(var, nbit) ((var) & (static_cast(1) << (nbit))) + +using FullTracksExt = soa::Join; + +// For MC association in pre-selection +using labeledTracks = soa::Join; +using tofTracks = soa::Join; +using richTracks = soa::Join; + +struct alice3decaypreselector { + Produces a3decayMaps; + + // Operation and minimisation criteria + Configurable nSigmaTOF{"nSigmaTOF", 4.0f, "Nsigma for TOF PID (if enabled)"}; + Configurable nSigmaRICH{"nSigmaRICH", 4.0f, "Nsigma for RICH PID (if enabled)"}; + + // Define o2 fitter, 2-prong, active memory (no need to redefine per event) + o2::vertexing::DCAFitterN<2> fitter; + + // for bit-packed maps + std::vector selectionMap; + + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + /// function to check PDG + PDG mother + template + bool checkPDG(TTrack const& track, int pdgMother, int pdg) + { + bool returnValue = false; + // Association check + if (track.has_mcParticle()) { + auto mcParticle = track.template mcParticle_as(); + if (mcParticle.has_mothers()) { + for (auto& mcParticleMother : mcParticle.template mothers_as()) { + if (mcParticle.pdgCode() == pdg && mcParticleMother.pdgCode() == pdgMother) + returnValue = true; + } + } + } // end association check + return returnValue; + } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + + void init(InitContext&) + { + // future dev if needed + } + + // go declarative: use partitions instead of "if", then just toggle bits to allow for mask selection later + Partition pInnerTOFPi = nabs(aod::upgrade_tof::nSigmaPionInnerTOF) > nSigmaTOF; + Partition pInnerTOFKa = nabs(aod::upgrade_tof::nSigmaKaonInnerTOF) > nSigmaTOF; + Partition pInnerTOFPr = nabs(aod::upgrade_tof::nSigmaProtonInnerTOF) > nSigmaTOF; + Partition pOuterTOFPi = nabs(aod::upgrade_tof::nSigmaPionOuterTOF) > nSigmaTOF; + Partition pOuterTOFKa = nabs(aod::upgrade_tof::nSigmaKaonOuterTOF) > nSigmaTOF; + Partition pOuterTOFPr = nabs(aod::upgrade_tof::nSigmaProtonOuterTOF) > nSigmaTOF; + Partition pRICHPi = nabs(aod::alice3rich::richNsigmaPi) > nSigmaRICH; + Partition pRICHKa = nabs(aod::alice3rich::richNsigmaKa) > nSigmaRICH; + Partition pRICHPr = nabs(aod::alice3rich::richNsigmaPr) > nSigmaRICH; + + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + /// Initialization of mask vectors if uninitialized + void initializeMasks(int size) + { + selectionMap.clear(); + selectionMap.resize(size, 0xFFFFFFFF); // all bits 1, please + } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + /// This process function ensures that all V0s are built. It will simply tag everything as true. + void processInitialize(aod::Tracks const& tracks) + { + initializeMasks(tracks.size()); + } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + void processFilterInnerTOF(tofTracks const&) + { + for (auto const& track : pInnerTOFPi) + bitoff(selectionMap[track.globalIndex()], kInnerTOFPion); + for (auto const& track : pInnerTOFKa) + bitoff(selectionMap[track.globalIndex()], kInnerTOFKaon); + for (auto const& track : pInnerTOFPr) + bitoff(selectionMap[track.globalIndex()], kInnerTOFProton); + } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + void processFilterOuterTOF(tofTracks const&) + { + for (auto const& track : pOuterTOFPi) + bitoff(selectionMap[track.globalIndex()], kOuterTOFPion); + for (auto const& track : pOuterTOFKa) + bitoff(selectionMap[track.globalIndex()], kOuterTOFKaon); + for (auto const& track : pOuterTOFPr) + bitoff(selectionMap[track.globalIndex()], kOuterTOFProton); + } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + void processFilterRICH(richTracks const&) + { + for (auto const& track : pRICHPi) + bitoff(selectionMap[track.globalIndex()], kRICHPion); + for (auto const& track : pRICHKa) + bitoff(selectionMap[track.globalIndex()], kRICHKaon); + for (auto const& track : pRICHPr) + bitoff(selectionMap[track.globalIndex()], kRICHProton); + } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + void processFilterOnMonteCarloTruth(labeledTracks const& tracks, aod::McParticles const&) + { + for (auto const& track : tracks) { + // D mesons + if (!checkPDG(track, 421, -321)) //+421 -> -321 +211 + bitoff(selectionMap[track.globalIndex()], kTrueKaMinusFromD); + if (!checkPDG(track, -421, +321)) //-421 -> +321 -211 + bitoff(selectionMap[track.globalIndex()], kTrueKaPlusFromD); + if (!checkPDG(track, 421, +211)) //+421 -> -321 +211 + bitoff(selectionMap[track.globalIndex()], kTruePiPlusFromD); + if (!checkPDG(track, -421, -211)) //-421 -> +321 -211 + bitoff(selectionMap[track.globalIndex()], kTruePiMinusFromD); + + // Lambdac baryons + if (!checkPDG(track, +4122, +2212)) //+4122 -> +2212 -321 +211 + bitoff(selectionMap[track.globalIndex()], kTruePrPlusFromLc); + if (!checkPDG(track, +4122, -321)) //+4122 -> +2212 -321 +211 + bitoff(selectionMap[track.globalIndex()], kTrueKaMinusFromLc); + if (!checkPDG(track, +4122, +211)) //+4122 -> +2212 -321 +211 + bitoff(selectionMap[track.globalIndex()], kTruePiPlusFromLc); + if (!checkPDG(track, -4122, -2212)) //-4122 -> -2212 +321 -211 + bitoff(selectionMap[track.globalIndex()], kTruePrMinusFromLc); + if (!checkPDG(track, -4122, +321)) //-4122 -> -2212 +321 -211 + bitoff(selectionMap[track.globalIndex()], kTrueKaPlusFromLc); + if (!checkPDG(track, -4122, -211)) //-4122 -> -2212 +321 -211 + bitoff(selectionMap[track.globalIndex()], kTruePiMinusFromLc); + + // XiCC daughters + if (!checkPDG(track, 4422, 211)) // 4422 -> 4232 -211, pi from xicc + bitoff(selectionMap[track.globalIndex()], kTruePiFromXiCC); + if (!checkPDG(track, 4232, 3312)) // 4232 -> 3312 211 211, xi from xic + bitoff(selectionMap[track.globalIndex()], kTrueXiFromXiC); + if (!checkPDG(track, 4232, 211)) // 4232 -> 3312 211 211, pi from xic + bitoff(selectionMap[track.globalIndex()], kTruePiFromXiC); + } + } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + void processPublishDecision(aod::Tracks const& tracks) + { + for (uint32_t i = 0; i < tracks.size(); i++) { + a3decayMaps(selectionMap[i]); + } + selectionMap.clear(); + } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + + //*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<* + PROCESS_SWITCH(alice3decaypreselector, processInitialize, "Initialize (MUST be on)", true); + PROCESS_SWITCH(alice3decaypreselector, processFilterInnerTOF, "Switch to use inner TOF PID", false); + PROCESS_SWITCH(alice3decaypreselector, processFilterOuterTOF, "Switch to use outer TOF PID", false); + PROCESS_SWITCH(alice3decaypreselector, processFilterRICH, "Switch to use RICH", false); + PROCESS_SWITCH(alice3decaypreselector, processFilterOnMonteCarloTruth, "Switch to use MC truth", false); + PROCESS_SWITCH(alice3decaypreselector, processPublishDecision, "Fill decision mask table (MUST be on)", true); + //*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<* +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +} diff --git a/ALICE3/TableProducer/alice3-multicharm.cxx b/ALICE3/TableProducer/alice3-multicharm.cxx new file mode 100644 index 00000000000..3cfdea25c00 --- /dev/null +++ b/ALICE3/TableProducer/alice3-multicharm.cxx @@ -0,0 +1,411 @@ +// 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. +// +// *+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* +// Decay finder task for ALICE 3 +// *+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* +// +// Uses specific ALICE 3 PID and performance for studying +// HF decays. Work in progress: use at your own risk! +// + +#include +#include +#include +#include +#include +#include + +#include "Framework/runDataProcessing.h" +#include "Framework/RunningWorkflowInfo.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoAHelpers.h" +#include "DCAFitter/DCAFitterN.h" +#include "ReconstructionDataFormats/Track.h" +#include "Common/Core/RecoDecay.h" +#include "Common/Core/trackUtilities.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/DataModel/LFParticleIdentification.h" +#include "Common/Core/TrackSelection.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "DetectorsBase/Propagator.h" +#include "DetectorsBase/GeometryManager.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "CCDB/BasicCCDBManager.h" +#include "DataFormatsCalibration/MeanVertexObject.h" +#include "ALICE3/DataModel/OTFTOF.h" +#include "ALICE3/DataModel/RICH.h" +#include "ALICE3/DataModel/A3DecayFinderTables.h" +#include "ALICE3/DataModel/OTFStrangeness.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using std::array; + +// simple checkers +// #define biton(var, nbit) ((var) |= (static_cast(1) << (nbit))) +#define bitoff(var, nbit) ((var) &= ~(static_cast(1) << (nbit))) //((a) &= ~(1ULL<<(b))) +// #define bitcheck(var, nbit) ((var) & (static_cast(1) << (nbit))) + +using FullTracksExt = soa::Join; + +// For MC association in pre-selection +using labeledTracks = soa::Join; +using tofTracks = soa::Join; +using richTracks = soa::Join; +using alice3tracks = soa::Join; + +struct alice3multicharm { + SliceCache cache; + + // Operation and minimisation criteria + Configurable magneticField{"magneticField", 20.0f, "Magnetic field (in kilogauss)"}; + Configurable doDCAplots{"doDCAplots", true, "do daughter prong DCA plots for D mesons"}; + Configurable mcSameMotherCheck{"mcSameMotherCheck", true, "check if tracks come from the same MC mother"}; + Configurable dcaXiCDaughtersSelection{"dcaXiCDaughtersSelection", 1000.0f, "DCA between XiC daughters (cm)"}; + Configurable dcaXiCCDaughtersSelection{"dcaXiCCDaughtersSelection", 1000.0f, "DCA between XiCC daughters (cm)"}; + + Configurable piFromXiC_dcaXYconstant{"piFromXiC_dcaXYconstant", -1.0f, "[0] in |DCAxy| > [0]+[1]/pT"}; + Configurable piFromXiC_dcaXYpTdep{"piFromXiC_dcaXYpTdep", 0.0, "[1] in |DCAxy| > [0]+[1]/pT"}; + Configurable piFromXiCC_dcaXYconstant{"piFromXiCC_dcaXYconstant", -1.0f, "[0] in |DCAxy| > [0]+[1]/pT"}; + Configurable piFromXiCC_dcaXYpTdep{"piFromXiCC_dcaXYpTdep", 0.0, "[1] in |DCAxy| > [0]+[1]/pT"}; + Configurable xiFromXiC_dcaXYconstant{"xiFromXiC_dcaXYconstant", -1.0f, "[0] in |DCAxy| > [0]+[1]/pT"}; + Configurable xiFromXiC_dcaXYpTdep{"xiFromXiC_dcaXYpTdep", 0.0, "[1] in |DCAxy| > [0]+[1]/pT"}; + + ConfigurableAxis axisEta{"axisEta", {8, -4.0f, +4.0f}, "#eta"}; + ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "pt axis for QA histograms"}; + ConfigurableAxis axisDCA{"axisDCA", {200, -100, 100}, "DCA (#mum)"}; + + ConfigurableAxis axisXiMass{"axisXiMass", {200, 1.221f, 1.421f}, "Xi Inv Mass (GeV/c^{2})"}; + ConfigurableAxis axisXiCMass{"axisXiCMass", {200, 2.368f, 2.568f}, "XiC Inv Mass (GeV/c^{2})"}; + ConfigurableAxis axisXiCCMass{"axisXiCCMass", {200, 3.521f, 3.721f}, "XiCC Inv Mass (GeV/c^{2})"}; + + ConfigurableAxis axisNConsidered{"axisNConsidered", {200, -0.5f, 199.5f}, "Number of considered track combinations"}; + + o2::vertexing::DCAFitterN<2> fitter; + o2::vertexing::DCAFitterN<3> fitter3; + + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + Partition trueXi = aod::mcparticle::pdgCode == 3312; + Partition trueXiC = aod::mcparticle::pdgCode == 4232; + Partition trueXiCC = aod::mcparticle::pdgCode == 4422; + + // filter expressions for D mesons + static constexpr uint32_t trackSelectionPiFromXiC = 1 << kInnerTOFPion | 1 << kOuterTOFPion | 1 << kRICHPion | 1 << kTruePiFromXiC; + static constexpr uint32_t trackSelectionPiFromXiCC = 1 << kInnerTOFPion | 1 << kOuterTOFPion | 1 << kRICHPion | 1 << kTruePiFromXiCC; + + // partitions for Xi daughters + Partition tracksPiFromXiC = + ((aod::a3DecayMap::decayMap & trackSelectionPiFromXiC) == trackSelectionPiFromXiC) && aod::track::signed1Pt > 0.0f && nabs(aod::track::dcaXY) > piFromXiC_dcaXYconstant + piFromXiC_dcaXYpTdep* nabs(aod::track::signed1Pt); + Partition tracksPiFromXiCC = + ((aod::a3DecayMap::decayMap & trackSelectionPiFromXiCC) == trackSelectionPiFromXiCC) && aod::track::signed1Pt > 0.0f && nabs(aod::track::dcaXY) > piFromXiCC_dcaXYconstant + piFromXiCC_dcaXYpTdep* nabs(aod::track::signed1Pt); + + // Helper struct to pass candidate information + struct { + float dca; + float mass; + float pt; + float eta; + std::array xyz; + std::array prong0mom; + std::array prong1mom; + std::array prong2mom; + std::array parentTrackCovMatrix; + } thisCandidate; + + template + bool buildDecayCandidateTwoBody(TTrackType const& t0, TTrackType const& t1, float mass0, float mass1) + { + //}-{}-{}-{}-{}-{}-{}-{}-{}-{} + // Move close to minima + int nCand = 0; + try { + nCand = fitter.process(t0, t1); + } catch (...) { + return false; + } + if (nCand == 0) { + return false; + } + //}-{}-{}-{}-{}-{}-{}-{}-{}-{} + + o2::track::TrackParCov t0new = fitter.getTrack(0); + o2::track::TrackParCov t1new = fitter.getTrack(1); + t0new.getPxPyPzGlo(thisCandidate.prong0mom); + t1new.getPxPyPzGlo(thisCandidate.prong1mom); + + // get decay vertex coordinates + const auto& vtx = fitter.getPCACandidate(); + for (int i = 0; i < 3; i++) { + thisCandidate.xyz[i] = vtx[i]; + } + + // compute cov mat + for (int ii = 0; ii < 21; ii++) + thisCandidate.parentTrackCovMatrix[ii] = 0.0f; + + std::array covA = {0}; + std::array covB = {0}; + fitter.getTrack(0).getCovXYZPxPyPzGlo(covA); + fitter.getTrack(1).getCovXYZPxPyPzGlo(covB); + + const int momInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component + for (int i = 0; i < 6; i++) { + int j = momInd[i]; + thisCandidate.parentTrackCovMatrix[j] = covA[j] + covB[j]; + } + + auto covVtx = fitter.calcPCACovMatrix(); + thisCandidate.parentTrackCovMatrix[0] = covVtx(0, 0); + thisCandidate.parentTrackCovMatrix[1] = covVtx(1, 0); + thisCandidate.parentTrackCovMatrix[2] = covVtx(1, 1); + thisCandidate.parentTrackCovMatrix[3] = covVtx(2, 0); + thisCandidate.parentTrackCovMatrix[4] = covVtx(2, 1); + thisCandidate.parentTrackCovMatrix[5] = covVtx(2, 2); + + // set relevant values + thisCandidate.dca = TMath::Sqrt(fitter.getChi2AtPCACandidate()); + thisCandidate.mass = RecoDecay::m(array{array{thisCandidate.prong0mom[0], thisCandidate.prong0mom[1], thisCandidate.prong0mom[2]}, array{thisCandidate.prong1mom[0], thisCandidate.prong1mom[1], thisCandidate.prong1mom[2]}}, array{mass0, mass1}); + thisCandidate.pt = std::hypot(thisCandidate.prong0mom[0] + thisCandidate.prong1mom[0], thisCandidate.prong0mom[1] + thisCandidate.prong1mom[1]); + thisCandidate.eta = RecoDecay::eta(array{thisCandidate.prong0mom[0] + thisCandidate.prong1mom[0], thisCandidate.prong0mom[1] + thisCandidate.prong1mom[1], thisCandidate.prong0mom[2] + thisCandidate.prong1mom[2]}); + return true; + } + + template + bool buildDecayCandidateThreeBody(TTrackType1 const& prong0, TTrackType2 const& prong1, TTrackType3 const& prong2, float p0mass, float p1mass, float p2mass) + { + o2::track::TrackParCov t0 = getTrackParCov(prong0); + o2::track::TrackParCov t1 = getTrackParCov(prong1); + o2::track::TrackParCov t2 = getTrackParCov(prong2); + + //}-{}-{}-{}-{}-{}-{}-{}-{}-{} + // Move close to minima + int nCand = 0; + try { + nCand = fitter3.process(t0, t1, t2); + } catch (...) { + return false; + } + if (nCand == 0) { + return false; + } + //}-{}-{}-{}-{}-{}-{}-{}-{}-{} + + t0 = fitter3.getTrack(0); + t1 = fitter3.getTrack(1); + t2 = fitter3.getTrack(2); + t0.getPxPyPzGlo(thisCandidate.prong0mom); + t1.getPxPyPzGlo(thisCandidate.prong1mom); + t2.getPxPyPzGlo(thisCandidate.prong2mom); + + // get decay vertex coordinates + const auto& vtx = fitter3.getPCACandidate(); + for (int i = 0; i < 3; i++) { + thisCandidate.xyz[i] = vtx[i]; + } + + // compute cov mat + for (int ii = 0; ii < 21; ii++) + thisCandidate.parentTrackCovMatrix[ii] = 0.0f; + + std::array covA = {0}; + std::array covB = {0}; + std::array covC = {0}; + fitter3.getTrack(0).getCovXYZPxPyPzGlo(covA); + fitter3.getTrack(1).getCovXYZPxPyPzGlo(covB); + fitter3.getTrack(2).getCovXYZPxPyPzGlo(covC); + + const int momInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component + for (int i = 0; i < 6; i++) { + int j = momInd[i]; + thisCandidate.parentTrackCovMatrix[j] = covA[j] + covB[j] + covC[j]; + } + + auto covVtx = fitter3.calcPCACovMatrix(); + thisCandidate.parentTrackCovMatrix[0] = covVtx(0, 0); + thisCandidate.parentTrackCovMatrix[1] = covVtx(1, 0); + thisCandidate.parentTrackCovMatrix[2] = covVtx(1, 1); + thisCandidate.parentTrackCovMatrix[3] = covVtx(2, 0); + thisCandidate.parentTrackCovMatrix[4] = covVtx(2, 1); + thisCandidate.parentTrackCovMatrix[5] = covVtx(2, 2); + + // set relevant values + thisCandidate.dca = TMath::Sqrt(fitter3.getChi2AtPCACandidate()); + thisCandidate.mass = RecoDecay::m(array{array{thisCandidate.prong0mom[0], thisCandidate.prong0mom[1], thisCandidate.prong0mom[2]}, array{thisCandidate.prong1mom[0], thisCandidate.prong1mom[1], thisCandidate.prong1mom[2]}, array{thisCandidate.prong2mom[0], thisCandidate.prong2mom[1], thisCandidate.prong2mom[2]}}, array{p0mass, p1mass, p2mass}); + thisCandidate.pt = std::hypot(thisCandidate.prong0mom[0] + thisCandidate.prong1mom[0] + thisCandidate.prong2mom[0], thisCandidate.prong0mom[1] + thisCandidate.prong1mom[1] + thisCandidate.prong2mom[1]); + thisCandidate.eta = RecoDecay::eta(array{thisCandidate.prong0mom[0] + thisCandidate.prong1mom[0] + thisCandidate.prong2mom[0], thisCandidate.prong0mom[1] + thisCandidate.prong1mom[1] + thisCandidate.prong2mom[1], thisCandidate.prong0mom[2] + thisCandidate.prong1mom[2] + thisCandidate.prong2mom[2]}); + return true; + } + + /// function to check if tracks have the same mother in MC + template + bool checkSameMother(TTrackType1 const& track1, TTrackType2 const& track2) + { + bool returnValue = false; + // Association check + // There might be smarter ways of doing this in the future + if (track1.has_mcParticle() && track2.has_mcParticle()) { + auto mcParticle1 = track1.template mcParticle_as(); + auto mcParticle2 = track2.template mcParticle_as(); + if (mcParticle1.has_mothers() && mcParticle2.has_mothers()) { + for (auto& mcParticleMother1 : mcParticle1.template mothers_as()) { + for (auto& mcParticleMother2 : mcParticle2.template mothers_as()) { + if (mcParticleMother1.globalIndex() == mcParticleMother2.globalIndex()) { + returnValue = true; + } + } + } + } + } // end association check + return returnValue; + } + + void init(InitContext&) + { + // initialize O2 2-prong fitter (only once) + fitter.setPropagateToPCA(true); + fitter.setMaxR(200.); + fitter.setMinParamChange(1e-3); + fitter.setMinRelChi2Change(0.9); + fitter.setMaxDZIni(1e9); + fitter.setMaxChi2(1e9); + fitter.setUseAbsDCA(true); + fitter.setWeightedFinalPCA(false); + fitter.setBz(magneticField); + fitter.setMatCorrType(o2::base::Propagator::MatCorrType::USEMatCorrNONE); + + fitter3.setPropagateToPCA(true); + fitter3.setMaxR(200.); + fitter3.setMinParamChange(1e-3); + fitter3.setMinRelChi2Change(0.9); + fitter3.setMaxDZIni(1e9); + fitter3.setMaxChi2(1e9); + fitter3.setUseAbsDCA(true); + fitter3.setWeightedFinalPCA(false); + fitter3.setBz(magneticField); + fitter3.setMatCorrType(o2::base::Propagator::MatCorrType::USEMatCorrNONE); + + histos.add("h2dGenXi", "h2dGenXi", kTH2F, {axisPt, axisEta}); + histos.add("h2dGenXiC", "h2dGenXiC", kTH2F, {axisPt, axisEta}); + histos.add("h2dGenXiCC", "h2dGenXiCC", kTH2F, {axisPt, axisEta}); + + histos.add("hMassXi", "hMassXi", kTH1F, {axisXiMass}); + histos.add("hMassXiC", "hMassXiC", kTH1F, {axisXiCMass}); + histos.add("hMassXiCC", "hMassXiCC", kTH1F, {axisXiCCMass}); + + histos.add("hCombinationsXiC", "hCombinationsXiC", kTH1F, {axisNConsidered}); + histos.add("hCombinationsXiCC", "hCombinationsXiCC", kTH1F, {axisNConsidered}); + + if (doDCAplots) { + histos.add("h2dDCAxyVsPtXiFromXiC", "h2dDCAxyVsPtXiFromXiC", kTH2F, {axisPt, axisDCA}); + histos.add("h2dDCAxyVsPtPiFromXiC", "h2dDCAxyVsPtPiFromXiC", kTH2F, {axisPt, axisDCA}); + histos.add("h2dDCAxyVsPtPiFromXiCC", "h2dDCAxyVsPtPiFromXiCC", kTH2F, {axisPt, axisDCA}); + } + } + + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + void processGenerated(aod::McParticles const&) + { + for (auto const& mcParticle : trueXi) + histos.fill(HIST("h2dGenXi"), mcParticle.pt(), mcParticle.eta()); + for (auto const& mcParticle : trueXiC) + histos.fill(HIST("h2dGenXiC"), mcParticle.pt(), mcParticle.eta()); + for (auto const& mcParticle : trueXiCC) + histos.fill(HIST("h2dGenXiCC"), mcParticle.pt(), mcParticle.eta()); + } + + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + void processFindXiCC(aod::Collision const& collision, alice3tracks const&, aod::McParticles const&, aod::UpgradeCascades const& cascades) + { + // group with this collision + // n.b. cascades do not need to be grouped, being used directly in iterator-grouping + auto tracksPiFromXiCgrouped = tracksPiFromXiC->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + auto tracksPiFromXiCCgrouped = tracksPiFromXiCC->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + + if (doDCAplots) { + for (auto const& cascade : cascades) { + if (cascade.has_cascadeTrack()) { + auto track = cascade.cascadeTrack_as(); // de-reference cascade track + histos.fill(HIST("h2dDCAxyVsPtXiFromXiC"), track.pt(), track.dcaXY() * 1e+4); + } else { + LOGF(info, "Damn, something is wrong"); + } + } + for (auto const& track : tracksPiFromXiCgrouped) + histos.fill(HIST("h2dDCAxyVsPtPiFromXiC"), track.pt(), track.dcaXY() * 1e+4); + for (auto const& track : tracksPiFromXiCCgrouped) + histos.fill(HIST("h2dDCAxyVsPtPiFromXiCC"), track.pt(), track.dcaXY() * 1e+4); + } + + for (auto const& xiCand : cascades) { + histos.fill(HIST("hMassXi"), xiCand.mXi()); + auto xi = xiCand.cascadeTrack_as(); // de-reference cascade track + uint32_t nCombinationsC = 0; + for (auto const& pi1c : tracksPiFromXiCgrouped) { + if (mcSameMotherCheck && !checkSameMother(xi, pi1c)) + continue; + for (auto const& pi2c : tracksPiFromXiCgrouped) { + if (mcSameMotherCheck && !checkSameMother(xi, pi2c)) + continue; // keep only if same mother + if (pi1c.globalIndex() >= pi2c.globalIndex()) + continue; // avoid same-mother, avoid double-counting + + // if I am here, it means this is a triplet to be considered for XiC vertexing. + // will now attempt to build a three-body decay candidate with these three track rows. + + nCombinationsC++; + if (!buildDecayCandidateThreeBody(xi, pi1c, pi2c, 1.32171, 0.139570, 0.139570)) + continue; // failed at building candidate + + const std::array momentumC = { + thisCandidate.prong0mom[0] + thisCandidate.prong1mom[0] + thisCandidate.prong2mom[0], + thisCandidate.prong0mom[1] + thisCandidate.prong1mom[1] + thisCandidate.prong2mom[1], + thisCandidate.prong0mom[2] + thisCandidate.prong1mom[2] + thisCandidate.prong2mom[2]}; + + o2::track::TrackParCov xicTrack(thisCandidate.xyz, momentumC, thisCandidate.parentTrackCovMatrix, +1); + + histos.fill(HIST("hMassXiC"), thisCandidate.mass); + + // attempt XiCC finding + uint32_t nCombinationsCC = 0; + for (auto const& picc : tracksPiFromXiCCgrouped) { + // to-do: check same mother here + + nCombinationsCC++; + o2::track::TrackParCov piccTrack = getTrackParCov(picc); + if (!buildDecayCandidateTwoBody(xicTrack, piccTrack, 2.46793, 0.139570)) + continue; // failed at building candidate + + histos.fill(HIST("hMassXiCC"), thisCandidate.mass); + } + histos.fill(HIST("hCombinationsXiCC"), nCombinationsCC); + } + } + histos.fill(HIST("hCombinationsXiC"), nCombinationsC); + } + } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + + //*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<* + PROCESS_SWITCH(alice3multicharm, processGenerated, "fill MC-only histograms", true); + PROCESS_SWITCH(alice3multicharm, processFindXiCC, "find XiCC baryons", true); + //*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<* +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +}