From abba92eb5e782ad230412e4978e49d9f75e53c7a Mon Sep 17 00:00:00 2001 From: ddobrigk Date: Tue, 6 Aug 2024 10:14:23 +0200 Subject: [PATCH 1/8] ALICE 3: Add OTF strangeness tracking --- ALICE3/DataModel/OTFStrangeness.h | 68 ++++ ALICE3/TableProducer/OTF/CMakeLists.txt | 2 +- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 329 +++++++++++++++---- 3 files changed, 326 insertions(+), 73 deletions(-) create mode 100644 ALICE3/DataModel/OTFStrangeness.h diff --git a/ALICE3/DataModel/OTFStrangeness.h b/ALICE3/DataModel/OTFStrangeness.h new file mode 100644 index 00000000000..7a6b462ccfc --- /dev/null +++ b/ALICE3/DataModel/OTFStrangeness.h @@ -0,0 +1,68 @@ +// 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_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::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/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..2f18b491681 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 { @@ -140,6 +160,28 @@ struct OnTheFlyTracker { bool isDecayDau; }; + // 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,12 +204,6 @@ 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; @@ -339,7 +375,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 +386,12 @@ 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}}); } LOGF(info, "Initializing magnetic field to value: %.3f kG", static_cast(magneticField)); @@ -388,6 +424,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 @@ -402,7 +450,7 @@ struct OnTheFlyTracker { std::vector l0Velocity(3); std::vector xiMomentum = {particle.px(), particle.py(), particle.pz()}; std::vector xiProductionVertex = {particle.vx(), particle.vy(), particle.vz()}; - double bz = magneticField * 0.1; // To tesla + double bz = magneticField * (-0.1); // To tesla (sign?!) TRandom3 rand; rand.SetSeed(seed); double u = rand.Uniform(0, 1); @@ -410,7 +458,7 @@ struct OnTheFlyTracker { xiVelocity[i] = xiMomentum[i] / particle.e(); } double xi_v_tot = particle.p() / particle.e(); - double xi_ctau = 4.91; // xi + double xi_ctau = 4.91 / 100; // xi double xi_charge = -1.6022e-19; double speedOfLight = 3e+8; double xi_v_xy = speedOfLight * sqrt(xiVelocity[0] * xiVelocity[0] + xiVelocity[1] * xiVelocity[1]); @@ -432,7 +480,7 @@ struct OnTheFlyTracker { xiDecay.Generate(); decayDaughters.push_back(*xiDecay.GetDecay(1)); - double l0_ctau = 7.89; // lambda + double l0_ctau = 7.89 / 100; // lambda double l0_v_tot = xiDecay.GetDecay(0)->P() / xiDecay.GetDecay(0)->E(); std::vector l0Daughters = {0.139, 0.938}; l0Velocity[0] = xiDecay.GetDecay(0)->Px() / xiDecay.GetDecay(0)->E(); @@ -509,12 +557,6 @@ struct OnTheFlyTracker { { tracksAlice3.clear(); ghostTracksAlice3.clear(); - trackPdg.clear(); - trackIsFromXi.clear(); - trackIsFromL0.clear(); - ghostTrackPdg.clear(); - ghostTrackIsFromXi.clear(); - ghostTrackIsFromL0.clear(); bcData.clear(); o2::dataformats::DCA dcaInfo; @@ -565,8 +607,8 @@ struct OnTheFlyTracker { if (treatXi) { if (mcParticle.pdgCode() == 3312) { decayParticle(mcParticle, decayProducts, xiDecayVertex, l0DecayVertex); - xiDecayRadius2D = sqrt(xiDecayVertex[0] * xiDecayVertex[0] + xiDecayVertex[1] * xiDecayVertex[1]); - l0DecayRadius2D = sqrt(l0DecayVertex[0] * l0DecayVertex[0] + l0DecayVertex[1] * l0DecayVertex[1]); + xiDecayRadius2D = sqrt(xiDecayVertex[0] * xiDecayVertex[0] + xiDecayVertex[1] * xiDecayVertex[1]) * 100; + l0DecayRadius2D = sqrt(l0DecayVertex[0] * l0DecayVertex[0] + l0DecayVertex[1] * l0DecayVertex[1]) * 100; } } @@ -599,7 +641,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 +652,9 @@ struct OnTheFlyTracker { continue; } + o2::track::TrackParCov trackParCov; + convertMCParticleToO2Track(mcParticle, trackParCov); + bool isDecayDaughter = false; if (mcParticle.getProcess() == 4) isDecayDaughter = true; @@ -618,9 +663,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) { @@ -671,18 +713,12 @@ struct OnTheFlyTracker { } 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]); } 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]); } } - 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 +728,192 @@ 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, true}); + + if(doXiQA){ + histos.fill(HIST("hMassLambda"), thisCascade.mLambda); + histos.fill(HIST("hMassXi"), thisCascade.mXi); + histos.fill(HIST("hFoundVsFindable"), thisCascade.findableClusters, thisCascade.foundClusters); + } + + // fill table + upgradeCascades( + thisCascade.cascadeTrackId, + thisCascade.positiveId, + thisCascade.negativeId, + thisCascade.bachelorId, + thisCascade.dcaV0dau, + thisCascade.dcacascdau, + thisCascade.v0radius, + thisCascade.cascradius, + thisCascade.cascradiusMC, + thisCascade.mLambda, + thisCascade.mXi, + thisCascade.findableClusters, + thisCascade.foundClusters + ); + } + } + } // 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 +945,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 +1023,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,20 +1038,6 @@ 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); - } - } tracksDCA(dcaXY, dcaZ); if (populateTracksDCACov) { tracksDCACov(dcaInfo.getSigmaY2(), dcaInfo.getSigmaZ2()); @@ -872,14 +1067,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 +1114,6 @@ 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++; } } }; From c1d491ea865081011e53c6bb04f25eaf08411a89 Mon Sep 17 00:00:00 2001 From: ddobrigk Date: Tue, 6 Aug 2024 13:11:46 +0200 Subject: [PATCH 2/8] Modifications for QA of DCAxy of secondaries --- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 30 +++++++++++++++++--- 1 file changed, 26 insertions(+), 4 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index 2f18b491681..56340354c7d 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -153,11 +153,18 @@ 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 @@ -392,6 +399,11 @@ struct OnTheFlyTracker { // 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)); @@ -712,9 +724,9 @@ struct OnTheFlyTracker { continue; } if (isReco[i]) { - tracksAlice3.push_back(TrackAlice3{xiDaughterTrackParCovs[i], mcParticle.globalIndex(), t, 100.f * 1e-3, true}); + 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}); + ghostTracksAlice3.push_back(TrackAlice3{xiDaughterTrackParCovs[i], mcParticle.globalIndex(), t, 100.f * 1e-3, true, true, i+2}); } } @@ -882,7 +894,7 @@ struct OnTheFlyTracker { // 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, true}); + tracksAlice3.push_back(TrackAlice3{cascadeTrack, mcParticle.globalIndex(), t, 100.f * 1e-3, false, false, 1}); if(doXiQA){ histos.fill(HIST("hMassLambda"), thisCascade.mLambda); @@ -1038,6 +1050,16 @@ struct OnTheFlyTracker { histos.fill(HIST("h2dDCAxy"), trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please histos.fill(HIST("hTrackXatDCA"), trackParametrization.getX()); } + 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) { tracksDCACov(dcaInfo.getSigmaY2(), dcaInfo.getSigmaZ2()); From bc40276f96c170181a1868bc46a9a1b315b804ee Mon Sep 17 00:00:00 2001 From: ALICE Builder Date: Tue, 6 Aug 2024 13:12:29 +0200 Subject: [PATCH 3/8] Please consider the following formatting changes (#316) --- ALICE3/DataModel/OTFStrangeness.h | 6 +- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 164 +++++++++---------- 2 files changed, 85 insertions(+), 85 deletions(-) diff --git a/ALICE3/DataModel/OTFStrangeness.h b/ALICE3/DataModel/OTFStrangeness.h index 7a6b462ccfc..9b1da51a26f 100644 --- a/ALICE3/DataModel/OTFStrangeness.h +++ b/ALICE3/DataModel/OTFStrangeness.h @@ -31,7 +31,7 @@ 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 +// topo vars DECLARE_SOA_COLUMN(DCAV0Daughters, dcaV0Daughters, float); DECLARE_SOA_COLUMN(DCACascadeDaughters, dcaCascadeDaughters, float); DECLARE_SOA_COLUMN(V0Radius, v0Radius, float); @@ -46,7 +46,7 @@ DECLARE_SOA_COLUMN(FoundClusters, foundClusters, int); } // namespace otfcascade DECLARE_SOA_TABLE(UpgradeCascades, "AOD", "UPGRADECASCADES", - o2::soa::Index<>, + o2::soa::Index<>, otfcascade::CascadeTrackId, otfcascade::PosTrackId, otfcascade::NegTrackId, @@ -57,7 +57,7 @@ DECLARE_SOA_TABLE(UpgradeCascades, "AOD", "UPGRADECASCADES", otfcascade::CascRadius, otfcascade::CascRadiusMC, otfcascade::MLambda, - otfcascade::MXi, + otfcascade::MXi, otfcascade::FindableClusters, otfcascade::FoundClusters); diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index 56340354c7d..de466cca6f7 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -141,7 +141,7 @@ struct OnTheFlyTracker { using PVertex = o2::dataformats::PrimaryVertex; - // for secondary vertex finding + // for secondary vertex finding o2::vertexing::DCAFitterN<2> fitter; // Class to hold the track information for the O2 vertexing @@ -153,35 +153,35 @@ 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, bool weakDecayDauInput = false, int isUsedInCascadingInput = 0) : o2::track::TrackParCov(src), - mcLabel{label}, - timeEst{t, te}, - isDecayDau(decayDauInput), - isWeakDecayDau(weakDecayDauInput), - isUsedInCascading(isUsedInCascadingInput) {} + 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 + 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 + // 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 - + 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; + // for tracking + int findableClusters; int foundClusters; float mLambda; @@ -724,9 +724,9 @@ struct OnTheFlyTracker { continue; } if (isReco[i]) { - tracksAlice3.push_back(TrackAlice3{xiDaughterTrackParCovs[i], mcParticle.globalIndex(), t, 100.f * 1e-3, true, true, i+2}); + 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, true, i+2}); + ghostTracksAlice3.push_back(TrackAlice3{xiDaughterTrackParCovs[i], mcParticle.globalIndex(), t, 100.f * 1e-3, true, true, i + 2}); } } @@ -742,61 +742,61 @@ struct OnTheFlyTracker { } // +-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+ - // combine particles into actual Xi candidate + // combine particles into actual Xi candidate // cascade building starts here - if(findXi && mcParticle.pdgCode() == 3312 && isReco[0] && isReco[1] && isReco[2]){ + 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; + // 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 + // 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!"; + // LOG(error) << "Exception caught in DCA fitter process call!"; dcaFitterOK_V0 = false; } if (nCand == 0) { dcaFitterOK_V0 = false; } // V0 found successfully - if(dcaFitterOK_V0){ + 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) + 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 + // 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 + // 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; + covV[i] = 1e-6; } o2::track::TrackParCov v0Track = o2::track::TrackParCov( {pos[0], pos[1], pos[2]}, @@ -805,58 +805,59 @@ struct OnTheFlyTracker { v0Track.setAbsCharge(0); v0Track.setPID(o2::track::PID::Lambda); - // dca fitter step + // 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!"; + // LOG(error) << "Exception caught in DCA fitter process call!"; dcaFitterOK_Cascade = false; } - if (nCand == 0){ + if (nCand == 0) { dcaFitterOK_Cascade = false; } // Cascade found successfully - if(dcaFitterOK_Cascade){ - o2::track::TrackParCov bachelorTrackAtPCA = fitter.getTrack(1); + 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 + // 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 + // 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 + 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.findableClusters = 0; thisCascade.foundClusters = 0; - if(trackXi){ - // optionally, add the points in the layers before the decay of the Xi + 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--) { + 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)){ + 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 } @@ -865,28 +866,28 @@ struct OnTheFlyTracker { 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...) + + 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 + // 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], + 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; + if (!(cascadeTrack.propagateTo(xyz1[0], magneticField))) + continue; const o2::track::TrackParametrization::dim2_t hitpoint = { static_cast(xyz1[1]), - static_cast(xyz1[2]) - }; + static_cast(xyz1[2])}; const o2::track::TrackParametrization::dim3_t hitpointcov = {pixelResolution * pixelResolution, 0.f, pixelResolution * pixelResolution}; - cascadeTrack.update(hitpoint,hitpointcov); + cascadeTrack.update(hitpoint, hitpointcov); thisCascade.foundClusters++; // add to findable } } @@ -896,31 +897,30 @@ struct OnTheFlyTracker { 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){ + if (doXiQA) { histos.fill(HIST("hMassLambda"), thisCascade.mLambda); histos.fill(HIST("hMassXi"), thisCascade.mXi); histos.fill(HIST("hFoundVsFindable"), thisCascade.findableClusters, thisCascade.foundClusters); } - // fill table - upgradeCascades( - thisCascade.cascadeTrackId, + // fill table + upgradeCascades( + thisCascade.cascadeTrackId, thisCascade.positiveId, - thisCascade.negativeId, - thisCascade.bachelorId, - thisCascade.dcaV0dau, - thisCascade.dcacascdau, - thisCascade.v0radius, - thisCascade.cascradius, - thisCascade.cascradiusMC, - thisCascade.mLambda, + thisCascade.negativeId, + thisCascade.bachelorId, + thisCascade.dcaV0dau, + thisCascade.dcacascdau, + thisCascade.v0radius, + thisCascade.cascradius, + thisCascade.cascradiusMC, + thisCascade.mLambda, thisCascade.mXi, - thisCascade.findableClusters, - thisCascade.foundClusters - ); + thisCascade.findableClusters, + thisCascade.foundClusters); } } - } // end cascade building + } // end cascade building // +-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+ continue; // Not filling the tables with the xi itself @@ -1050,14 +1050,14 @@ struct OnTheFlyTracker { histos.fill(HIST("h2dDCAxy"), trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please histos.fill(HIST("hTrackXatDCA"), trackParametrization.getX()); } - if (doXiQA){ - if(trackParCov.isUsedInCascading==1) + if (doXiQA) { + if (trackParCov.isUsedInCascading == 1) histos.fill(HIST("h2dDCAxyCascade"), trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please - if(trackParCov.isUsedInCascading==2) + if (trackParCov.isUsedInCascading == 2) histos.fill(HIST("h2dDCAxyCascadeBachelor"), trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please - if(trackParCov.isUsedInCascading==3) + if (trackParCov.isUsedInCascading == 3) histos.fill(HIST("h2dDCAxyCascadeNegative"), trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please - if(trackParCov.isUsedInCascading==4) + if (trackParCov.isUsedInCascading == 4) histos.fill(HIST("h2dDCAxyCascadePositive"), trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please } tracksDCA(dcaXY, dcaZ); From 68b45d7843a15c7434711518e7efded50915a414 Mon Sep 17 00:00:00 2001 From: ddobrigk Date: Tue, 6 Aug 2024 15:02:32 +0200 Subject: [PATCH 4/8] Add machinery to do multi-charm --- ALICE3/DataModel/A3DecayFinderTables.h | 7 +- ALICE3/DataModel/OTFMulticharm.h | 57 +++ ALICE3/DataModel/OTFStrangeness.h | 2 + ALICE3/TableProducer/CMakeLists.txt | 9 + ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 36 +- ALICE3/TableProducer/alice3-decayfinder.cxx | 144 -------- .../TableProducer/alice3-decaypreselector.cxx | 222 ++++++++++++ ALICE3/TableProducer/alice3-multicharm.cxx | 331 ++++++++++++++++++ 8 files changed, 647 insertions(+), 161 deletions(-) create mode 100644 ALICE3/DataModel/OTFMulticharm.h create mode 100644 ALICE3/TableProducer/alice3-decaypreselector.cxx create mode 100644 ALICE3/TableProducer/alice3-multicharm.cxx diff --git a/ALICE3/DataModel/A3DecayFinderTables.h b/ALICE3/DataModel/A3DecayFinderTables.h index cc7e8f3d3c7..03183136f49 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 index 9b1da51a26f..861ec3a7af8 100644 --- a/ALICE3/DataModel/OTFStrangeness.h +++ b/ALICE3/DataModel/OTFStrangeness.h @@ -26,6 +26,7 @@ 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"); //! @@ -47,6 +48,7 @@ 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, diff --git a/ALICE3/TableProducer/CMakeLists.txt b/ALICE3/TableProducer/CMakeLists.txt index 52300bfc644..89196ce280a 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-decayfinder.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/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index de466cca6f7..3f4dd0de8eb 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -214,6 +214,7 @@ struct OnTheFlyTracker { std::vector bcData; o2::steer::InteractionSampler irSampler; o2::vertexing::PVertexer vertexer; + std::vector cascadesAlice3; void init(o2::framework::InitContext&) { @@ -903,21 +904,8 @@ struct OnTheFlyTracker { histos.fill(HIST("hFoundVsFindable"), thisCascade.findableClusters, thisCascade.foundClusters); } - // fill table - upgradeCascades( - thisCascade.cascadeTrackId, - thisCascade.positiveId, - thisCascade.negativeId, - thisCascade.bachelorId, - thisCascade.dcaV0dau, - thisCascade.dcacascdau, - thisCascade.v0radius, - thisCascade.cascradius, - thisCascade.cascradiusMC, - thisCascade.mLambda, - thisCascade.mXi, - thisCascade.findableClusters, - thisCascade.foundClusters); + // add this cascade to vector (will fill cursor later with collision ID) + cascadesAlice3.push_back(thisCascade); } } } // end cascade building @@ -1137,6 +1125,24 @@ struct OnTheFlyTracker { } TracksAlice3(false); } + + 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..3079deb3a21 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; @@ -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..c256524241d --- /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, 4232, 211)) //4422 -> 4232 -211 + bitoff(selectionMap[track.globalIndex()], kTruePiFromXiC); + if (!checkPDG(track, 4232, 3312)) //4232 -> 3312 211 211 + bitoff(selectionMap[track.globalIndex()], kTrueXiFromXiC); + if (!checkPDG(track, 4232, 211)) //4232 -> 3312 211 211 + 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..0ac23ba61f0 --- /dev/null +++ b/ALICE3/TableProducer/alice3-multicharm.cxx @@ -0,0 +1,331 @@ +// 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})"}; + + 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 trackSelectionXiFromXiC = 1 << kTrueXiFromXiC; + 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; + } thisCandidate; + + template + bool buildDecayCandidateTwoBody(TTrackType const& t0, TTrackType const& t1, float mass0, float mass1, o2::track::TrackParCov& parentTrack) + { + //}-{}-{}-{}-{}-{}-{}-{}-{}-{} + // 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); + std::array P0; + std::array P1; + t0new.getPxPyPzGlo(P0); + t1new.getPxPyPzGlo(P1); + + // set relevant values + thisCandidate.dca = TMath::Sqrt(fitter.getChi2AtPCACandidate()); + thisCandidate.mass = RecoDecay::m(array{array{P0[0], P0[1], P0[2]}, array{P1[0], P1[1], P1[2]}}, array{mass0, mass1}); + thisCandidate.pt = std::hypot(P0[0] + P1[0], P0[1] + P1[1]); + thisCandidate.eta = RecoDecay::eta(array{P0[0] + P1[0], P0[1] + P1[1], P0[2] + P1[2]}); + parentTrack = fitter.createParentTrackParCov(); + return true; + } + + template + bool buildDecayCandidateThreeBody(TTrackType1 const& prong0, TTrackType2 const& prong1, TTrackType3 const& prong2, float p0mass, float p1mass, float p2mass, o2::track::TrackParCov& parentTrack) + { + 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 = fitter.getTrack(0); + t1 = fitter.getTrack(1); + t2 = fitter.getTrack(2); + std::array P0; + std::array P1; + std::array P2; + t0.getPxPyPzGlo(P0); + t1.getPxPyPzGlo(P1); + t2.getPxPyPzGlo(P2); + + // set relevant values + thisCandidate.dca = TMath::Sqrt(fitter3.getChi2AtPCACandidate()); + thisCandidate.mass = RecoDecay::m(array{array{P0[0], P0[1], P0[2]}, array{P1[0], P1[1], P1[2]}, array{P2[0], P2[1], P2[2]}}, array{p0mass, p1mass, p2mass}); + thisCandidate.pt = std::hypot(P0[0] + P1[0] + P2[0], P0[1] + P1[1] + P2[1]); + thisCandidate.eta = RecoDecay::eta(array{P0[0] + P1[0] + P2[0], P0[1] + P1[1] + P2[1], P0[2] + P1[2] + P2[2]}); + parentTrack = fitter3.createParentTrackParCov(); + 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}); + + 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){ + auto track = cascade.cascadeTrack_as(); // de-reference cascade track + histos.fill(HIST("h2dDCAxyVsPtXiFromXiC"), track.pt(), track.dcaXY() * 1e+4); + } + 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); + } + + // D mesons + for (auto const& xiCand : cascades) { + histos.fill(HIST("hMassXi"), xiCand.mXi()); + auto xi = xiCand.cascadeTrack_as(); // de-reference cascade track + for (auto const& pi1c : tracksPiFromXiCgrouped) { + if (mcSameMotherCheck && !checkSameMother(xi, pi1c)) + continue; + for (auto const& pi2c : tracksPiFromXiCgrouped) { + if (mcSameMotherCheck && !checkSameMother(xi, pi2c)) + continue; + // 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. + o2::track::TrackParCov xicTrack; + if(!buildDecayCandidateThreeBody(xi, pi1c, pi2c, 1.32171, 0.139570, 0.139570, xicTrack)) + continue; // failed at building candidate + + histos.fill(HIST("hMassXiC"), thisCandidate.mass); + + // attempt XiCC finding + for (auto const& picc : tracksPiFromXiCCgrouped) { + o2::track::TrackParCov piccTrack = getTrackParCov(picc); + + o2::track::TrackParCov xiccTrack; + if(!buildDecayCandidateTwoBody(xicTrack, piccTrack, 2.46793, 0.139570, xiccTrack)) + continue; // failed at building candidate + + histos.fill(HIST("hMassXiCC"), thisCandidate.mass); + } + } + } + } + } + //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* + + + //*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<* + PROCESS_SWITCH(alice3multicharm, processGenerated, "fill MC-only histograms", true); + PROCESS_SWITCH(alice3multicharm, processFindXiCC, "find XiCC mesons", true); + //*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<* +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +} From 31ebb16703bb3c41187273f32d05d04d58cd0902 Mon Sep 17 00:00:00 2001 From: ddobrigk Date: Tue, 6 Aug 2024 15:56:56 +0200 Subject: [PATCH 5/8] Latest additions --- ALICE3/TableProducer/CMakeLists.txt | 2 +- ALICE3/TableProducer/alice3-decayfinder.cxx | 6 +- ALICE3/TableProducer/alice3-multicharm.cxx | 118 +++++++++++++++----- 3 files changed, 94 insertions(+), 32 deletions(-) diff --git a/ALICE3/TableProducer/CMakeLists.txt b/ALICE3/TableProducer/CMakeLists.txt index 89196ce280a..9c3d710d358 100644 --- a/ALICE3/TableProducer/CMakeLists.txt +++ b/ALICE3/TableProducer/CMakeLists.txt @@ -32,7 +32,7 @@ o2physics_add_dpl_workflow(alice3-centrality COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(alice3-decaypreselector - SOURCES alice3-decayfinder.cxx + SOURCES alice3-decaypreselector.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter COMPONENT_NAME Analysis) diff --git a/ALICE3/TableProducer/alice3-decayfinder.cxx b/ALICE3/TableProducer/alice3-decayfinder.cxx index 3079deb3a21..4b88305fb60 100644 --- a/ALICE3/TableProducer/alice3-decayfinder.cxx +++ b/ALICE3/TableProducer/alice3-decayfinder.cxx @@ -212,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; diff --git a/ALICE3/TableProducer/alice3-multicharm.cxx b/ALICE3/TableProducer/alice3-multicharm.cxx index 0ac23ba61f0..33f5721fe3c 100644 --- a/ALICE3/TableProducer/alice3-multicharm.cxx +++ b/ALICE3/TableProducer/alice3-multicharm.cxx @@ -117,10 +117,15 @@ struct alice3multicharm { float mass; float pt; float eta; + std::array xyz; + std::array prong0mom; + std::array prong1mom; + std::array prong2mom; + std::array parentTrackCovMatrix = {0}; } thisCandidate; template - bool buildDecayCandidateTwoBody(TTrackType const& t0, TTrackType const& t1, float mass0, float mass1, o2::track::TrackParCov& parentTrack) + bool buildDecayCandidateTwoBody(TTrackType const& t0, TTrackType const& t1, float mass0, float mass1) { //}-{}-{}-{}-{}-{}-{}-{}-{}-{} // Move close to minima @@ -137,22 +142,48 @@ struct alice3multicharm { o2::track::TrackParCov t0new = fitter.getTrack(0); o2::track::TrackParCov t1new = fitter.getTrack(1); - std::array P0; - std::array P1; - t0new.getPxPyPzGlo(P0); - t1new.getPxPyPzGlo(P1); + 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{P0[0], P0[1], P0[2]}, array{P1[0], P1[1], P1[2]}}, array{mass0, mass1}); - thisCandidate.pt = std::hypot(P0[0] + P1[0], P0[1] + P1[1]); - thisCandidate.eta = RecoDecay::eta(array{P0[0] + P1[0], P0[1] + P1[1], P0[2] + P1[2]}); - parentTrack = fitter.createParentTrackParCov(); + 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& parentTrack) + 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); @@ -171,22 +202,49 @@ struct alice3multicharm { } //}-{}-{}-{}-{}-{}-{}-{}-{}-{} - t0 = fitter.getTrack(0); - t1 = fitter.getTrack(1); - t2 = fitter.getTrack(2); - std::array P0; - std::array P1; - std::array P2; - t0.getPxPyPzGlo(P0); - t1.getPxPyPzGlo(P1); - t2.getPxPyPzGlo(P2); + 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{P0[0], P0[1], P0[2]}, array{P1[0], P1[1], P1[2]}, array{P2[0], P2[1], P2[2]}}, array{p0mass, p1mass, p2mass}); - thisCandidate.pt = std::hypot(P0[0] + P1[0] + P2[0], P0[1] + P1[1] + P2[1]); - thisCandidate.eta = RecoDecay::eta(array{P0[0] + P1[0] + P2[0], P0[1] + P1[1] + P2[1], P0[2] + P1[2] + P2[2]}); - parentTrack = fitter3.createParentTrackParCov(); + 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; } @@ -295,18 +353,22 @@ struct alice3multicharm { continue; // 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. - o2::track::TrackParCov xicTrack; - if(!buildDecayCandidateThreeBody(xi, pi1c, pi2c, 1.32171, 0.139570, 0.139570, xicTrack)) + 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 for (auto const& picc : tracksPiFromXiCCgrouped) { o2::track::TrackParCov piccTrack = getTrackParCov(picc); - - o2::track::TrackParCov xiccTrack; - if(!buildDecayCandidateTwoBody(xicTrack, piccTrack, 2.46793, 0.139570, xiccTrack)) + if(!buildDecayCandidateTwoBody(xicTrack, piccTrack, 2.46793, 0.139570)) continue; // failed at building candidate histos.fill(HIST("hMassXiCC"), thisCandidate.mass); From e467e96a5b111ca4db69da24f588cd34918e88c4 Mon Sep 17 00:00:00 2001 From: ALICE Builder Date: Tue, 6 Aug 2024 15:57:37 +0200 Subject: [PATCH 6/8] Please consider the following formatting changes (#317) --- ALICE3/DataModel/A3DecayFinderTables.h | 4 +-- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 2 +- .../TableProducer/alice3-decaypreselector.cxx | 6 ++-- ALICE3/TableProducer/alice3-multicharm.cxx | 29 +++++++++---------- 4 files changed, 20 insertions(+), 21 deletions(-) diff --git a/ALICE3/DataModel/A3DecayFinderTables.h b/ALICE3/DataModel/A3DecayFinderTables.h index 03183136f49..05eb5ad9307 100644 --- a/ALICE3/DataModel/A3DecayFinderTables.h +++ b/ALICE3/DataModel/A3DecayFinderTables.h @@ -44,8 +44,8 @@ enum a3selectionBit : uint32_t { kDCAxy = 0, kTruePiMinusFromLc, kTrueKaMinusFromLc, kTruePrMinusFromLc, - kTrueXiFromXiC, - kTruePiFromXiC, + kTrueXiFromXiC, + kTruePiFromXiC, kTruePiFromXiCC }; namespace o2::aod diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index 3f4dd0de8eb..c524000915a 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -1126,7 +1126,7 @@ struct OnTheFlyTracker { TracksAlice3(false); } - for(const auto& cascade : cascadesAlice3){ + for (const auto& cascade : cascadesAlice3) { upgradeCascades( collisions.lastIndex(), // now we know the collision index -> populate table cascade.cascadeTrackId, diff --git a/ALICE3/TableProducer/alice3-decaypreselector.cxx b/ALICE3/TableProducer/alice3-decaypreselector.cxx index c256524241d..6b41d361468 100644 --- a/ALICE3/TableProducer/alice3-decaypreselector.cxx +++ b/ALICE3/TableProducer/alice3-decaypreselector.cxx @@ -187,11 +187,11 @@ struct alice3decaypreselector { bitoff(selectionMap[track.globalIndex()], kTruePiMinusFromLc); // XiCC daughters - if (!checkPDG(track, 4232, 211)) //4422 -> 4232 -211 + if (!checkPDG(track, 4232, 211)) // 4422 -> 4232 -211 bitoff(selectionMap[track.globalIndex()], kTruePiFromXiC); - if (!checkPDG(track, 4232, 3312)) //4232 -> 3312 211 211 + if (!checkPDG(track, 4232, 3312)) // 4232 -> 3312 211 211 bitoff(selectionMap[track.globalIndex()], kTrueXiFromXiC); - if (!checkPDG(track, 4232, 211)) //4232 -> 3312 211 211 + if (!checkPDG(track, 4232, 211)) // 4232 -> 3312 211 211 bitoff(selectionMap[track.globalIndex()], kTruePiFromXiC); } } diff --git a/ALICE3/TableProducer/alice3-multicharm.cxx b/ALICE3/TableProducer/alice3-multicharm.cxx index 33f5721fe3c..458f97da85f 100644 --- a/ALICE3/TableProducer/alice3-multicharm.cxx +++ b/ALICE3/TableProducer/alice3-multicharm.cxx @@ -151,9 +151,9 @@ struct alice3multicharm { thisCandidate.xyz[i] = vtx[i]; } - // compute cov mat - for (int ii=0; ii<21; ii++) - thisCandidate.parentTrackCovMatrix[ii] = 0.0f; + // compute cov mat + for (int ii = 0; ii < 21; ii++) + thisCandidate.parentTrackCovMatrix[ii] = 0.0f; std::array covA = {0}; std::array covB = {0}; @@ -163,7 +163,7 @@ struct alice3multicharm { 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]; + thisCandidate.parentTrackCovMatrix[j] = covA[j] + covB[j]; } auto covVtx = fitter.calcPCACovMatrix(); @@ -215,9 +215,9 @@ struct alice3multicharm { thisCandidate.xyz[i] = vtx[i]; } - // compute cov mat - for (int ii=0; ii<21; ii++) - thisCandidate.parentTrackCovMatrix[ii] = 0.0f; + // compute cov mat + for (int ii = 0; ii < 21; ii++) + thisCandidate.parentTrackCovMatrix[ii] = 0.0f; std::array covA = {0}; std::array covB = {0}; @@ -229,7 +229,7 @@ struct alice3multicharm { 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]; + thisCandidate.parentTrackCovMatrix[j] = covA[j] + covB[j] + covC[j]; } auto covVtx = fitter3.calcPCACovMatrix(); @@ -331,7 +331,7 @@ struct alice3multicharm { auto tracksPiFromXiCCgrouped = tracksPiFromXiCC->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); if (doDCAplots) { - for (auto const& cascade : cascades){ + for (auto const& cascade : cascades) { auto track = cascade.cascadeTrack_as(); // de-reference cascade track histos.fill(HIST("h2dDCAxyVsPtXiFromXiC"), track.pt(), track.dcaXY() * 1e+4); } @@ -351,9 +351,9 @@ struct alice3multicharm { for (auto const& pi2c : tracksPiFromXiCgrouped) { if (mcSameMotherCheck && !checkSameMother(xi, pi2c)) continue; - // 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. - if(!buildDecayCandidateThreeBody(xi, pi1c, pi2c, 1.32171, 0.139570, 0.139570)) + // 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. + if (!buildDecayCandidateThreeBody(xi, pi1c, pi2c, 1.32171, 0.139570, 0.139570)) continue; // failed at building candidate const std::array momentumC = { @@ -365,10 +365,10 @@ struct alice3multicharm { histos.fill(HIST("hMassXiC"), thisCandidate.mass); - // attempt XiCC finding + // attempt XiCC finding for (auto const& picc : tracksPiFromXiCCgrouped) { o2::track::TrackParCov piccTrack = getTrackParCov(picc); - if(!buildDecayCandidateTwoBody(xicTrack, piccTrack, 2.46793, 0.139570)) + if (!buildDecayCandidateTwoBody(xicTrack, piccTrack, 2.46793, 0.139570)) continue; // failed at building candidate histos.fill(HIST("hMassXiCC"), thisCandidate.mass); @@ -379,7 +379,6 @@ struct alice3multicharm { } //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+* - //*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<*>-~-<* PROCESS_SWITCH(alice3multicharm, processGenerated, "fill MC-only histograms", true); PROCESS_SWITCH(alice3multicharm, processFindXiCC, "find XiCC mesons", true); From 1c972164080b3c8ec25fc738f94af749b658c345 Mon Sep 17 00:00:00 2001 From: Jesper Gumprecht <113693781+jesgum@users.noreply.github.com> Date: Tue, 6 Aug 2024 16:00:58 +0200 Subject: [PATCH 7/8] Minor fix --- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index c524000915a..d400343eb71 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -463,7 +463,7 @@ struct OnTheFlyTracker { std::vector l0Velocity(3); std::vector xiMomentum = {particle.px(), particle.py(), particle.pz()}; std::vector xiProductionVertex = {particle.vx(), particle.vy(), particle.vz()}; - double bz = magneticField * (-0.1); // To tesla (sign?!) + double bz = magneticField * 0.1; // To tesla TRandom3 rand; rand.SetSeed(seed); double u = rand.Uniform(0, 1); @@ -620,8 +620,8 @@ struct OnTheFlyTracker { if (treatXi) { if (mcParticle.pdgCode() == 3312) { decayParticle(mcParticle, decayProducts, xiDecayVertex, l0DecayVertex); - xiDecayRadius2D = sqrt(xiDecayVertex[0] * xiDecayVertex[0] + xiDecayVertex[1] * xiDecayVertex[1]) * 100; - l0DecayRadius2D = sqrt(l0DecayVertex[0] * l0DecayVertex[0] + l0DecayVertex[1] * l0DecayVertex[1]) * 100; + xiDecayRadius2D = sqrt(xiDecayVertex[0] * xiDecayVertex[0] + xiDecayVertex[1] * xiDecayVertex[1]); + l0DecayRadius2D = sqrt(l0DecayVertex[0] * l0DecayVertex[0] + l0DecayVertex[1] * l0DecayVertex[1]); } } From f9bf86a9873f588abbcd198bd4723020ce9c82da Mon Sep 17 00:00:00 2001 From: Jesper Gumprecht <113693781+jesgum@users.noreply.github.com> Date: Tue, 6 Aug 2024 16:05:21 +0200 Subject: [PATCH 8/8] Add Xi QA plots --- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index d400343eb71..08b591a3593 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -394,6 +394,13 @@ 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}); @@ -1039,6 +1046,18 @@ struct OnTheFlyTracker { histos.fill(HIST("hTrackXatDCA"), trackParametrization.getX()); } if (doXiQA) { + 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 (trackParCov.isUsedInCascading == 1) histos.fill(HIST("h2dDCAxyCascade"), trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please if (trackParCov.isUsedInCascading == 2)