From 8802460561b8f590243f8d1ab52ebf6715d53f03 Mon Sep 17 00:00:00 2001 From: marcellocosti Date: Tue, 28 Apr 2026 14:11:38 +0200 Subject: [PATCH 1/5] Refactor OTF tracker --- ALICE3/DataModel/A3DecayFinderTables.h | 10 + ALICE3/DataModel/tracksAlice3.h | 4 +- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 1930 ++++++++--------- ALICE3/TableProducer/alice3-decayfinder.cxx | 8 +- .../alice3TrackingTranslator.cxx | 3 +- ALICE3/Tasks/alice3HfTask3Prong.cxx | 39 +- Common/TableProducer/qVectorsTable.cxx | 6 +- PWGJE/Tasks/CMakeLists.txt | 18 +- 8 files changed, 979 insertions(+), 1039 deletions(-) diff --git a/ALICE3/DataModel/A3DecayFinderTables.h b/ALICE3/DataModel/A3DecayFinderTables.h index bd9aa5592b3..2ed4796045c 100644 --- a/ALICE3/DataModel/A3DecayFinderTables.h +++ b/ALICE3/DataModel/A3DecayFinderTables.h @@ -202,6 +202,11 @@ DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterXY, impactParameterXY, //! DECLARE_SOA_COLUMN(MlScore0, mlScore0, float); //! DECLARE_SOA_COLUMN(MlScore1, mlScore1, float); //! DECLARE_SOA_COLUMN(MlScore2, mlScore2, float); //! + +// Secondary vertex gen-reco position +DECLARE_SOA_COLUMN(DeltaXSecVtx, deltaXSecVtx, float); //! +DECLARE_SOA_COLUMN(DeltaYSecVtx, deltaYSecVtx, float); //! +DECLARE_SOA_COLUMN(DeltaZSecVtx, deltaZSecVtx, float); //! } // namespace a3_hf_cand #define HFCAND_COLUMNS \ @@ -235,6 +240,11 @@ DECLARE_SOA_COLUMN(MlScore2, mlScore2, float); //! a3_hf_cand::Pt2Prong2, \ a3_hf_cand::PVectorProng2 +DECLARE_SOA_TABLE(Alice3SVResos, "AOD", "ALICE3SVRESO", //! + a3_hf_cand::DeltaXSecVtx, + a3_hf_cand::DeltaYSecVtx, + a3_hf_cand::DeltaZSecVtx); + namespace a3DecayMap { DECLARE_SOA_COLUMN(DecayMap, decayMap, uint32_t); //! simple map to process passing / not passing criteria diff --git a/ALICE3/DataModel/tracksAlice3.h b/ALICE3/DataModel/tracksAlice3.h index 04cb2d1a9d4..374069de4cb 100644 --- a/ALICE3/DataModel/tracksAlice3.h +++ b/ALICE3/DataModel/tracksAlice3.h @@ -30,6 +30,7 @@ DECLARE_SOA_COLUMN(IsReconstructed, isReconstructed, bool); //! is reconstructed DECLARE_SOA_COLUMN(NSiliconHits, nSiliconHits, int); //! number of silicon hits DECLARE_SOA_COLUMN(NTPCHits, nTPCHits, int); //! number of tpc hits DECLARE_SOA_COLUMN(PdgCode, pdgCode, int); //! PDG code of the linked truth MC particle +DECLARE_SOA_COLUMN(TrackType, trackType, int); //! Type of the track } // namespace track_alice3 DECLARE_SOA_TABLE(TracksAlice3, "AOD", "TRACKSALICE3", track_alice3::IsReconstructed); @@ -41,7 +42,8 @@ using TrackAlice3Pdg = TracksAlice3Pdg::iterator; DECLARE_SOA_TABLE(TracksExtraA3, "AOD", "TracksExtraA3", track_alice3::NSiliconHits, - track_alice3::NTPCHits); + track_alice3::NTPCHits, + track_alice3::TrackType); using TrackExtraA3 = TracksExtraA3::iterator; namespace mcparticle_alice3 diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index 7cb3139aaf1..428d5d4bf33 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -116,6 +116,19 @@ using std::array; LOG(fatal) << "Histogram " << name << " not found!"; #define insertHist(name, ...) histPointers[name] = histos.add((name).c_str(), __VA_ARGS__); +enum TrackType { + kNone = 0, + kRecoPrimary, + kGhostPrimary, + kRecoV0Daug, + kGhostV0Daug, + kRecoCasc, + kGenCasc, + kRecoCascDaug, + kGenCascDaug, + kNTrackTypes, +}; + struct OnTheFlyTracker { Produces tableCollisions; Produces tableMcCollisionLabels; @@ -261,14 +274,17 @@ struct OnTheFlyTracker { bool weakDecayDauInput = false, int isUsedInCascadingInput = 0, int nSiliconHitsInput = 0, - int nTPCHitsInput = 0) : o2::track::TrackParCov(src), + int nTPCHitsInput = 0, + TrackType trackTypeInput = TrackType::kRecoPrimary + ) : o2::track::TrackParCov(src), mcLabel{label}, timeEst{time, timeError}, isDecayDau(decayDauInput), isWeakDecayDau(weakDecayDauInput), isUsedInCascading(isUsedInCascadingInput), nSiliconHits(nSiliconHitsInput), - nTPCHits(nTPCHitsInput) {} + nTPCHits(nTPCHitsInput), + trackType(trackTypeInput) {} const TimeEst& getTimeMUS() const { return timeEst; } int64_t mcLabel; TimeEst timeEst; ///< time estimate in ns @@ -277,6 +293,7 @@ struct OnTheFlyTracker { int isUsedInCascading; // 0: not at all, 1: is a cascade, 2: is a bachelor, 3: is a pion, 4: is a proton int nSiliconHits; int nTPCHits; + TrackType trackType; }; // Helper struct to pass cascade information @@ -336,6 +353,9 @@ struct OnTheFlyTracker { o2::constants::physics::kHelium3, o2::constants::physics::kAlpha}; + // Primary vertexing QA + std::pair vertexReconstructionEfficiencyCounters = {0, 0}; // {nVerticesWithMoreThan2Contributors, nVerticesReconstructed} + // necessary for particle charges Service pdgDB; @@ -349,6 +369,10 @@ struct OnTheFlyTracker { std::vector> mSmearer; // For processing and vertexing + std::vector recoPrimaries; + std::vector ghostPrimaries; + std::vector recoV0Daugs; + std::vector recoCascDaugs; std::vector tracksAlice3; std::vector ghostTracksAlice3; std::vector bcData; @@ -364,6 +388,14 @@ struct OnTheFlyTracker { // Configuration defined at init time o2::fastsim::GeometryContainer mGeoContainer; float mMagneticField = 0.0f; + // Time resolution constants + const float timeResolutionNs = 100.f; // ns + const float nsToMus = 1e-3f; + const float timeResolutionUs = timeResolutionNs * nsToMus; // us + + o2::dataformats::DCA dcaInfo; + o2::dataformats::VertexBase vtx; + void init(o2::framework::InitContext& initContext) { LOG(info) << "Initializing OnTheFlyTracker task"; @@ -452,9 +484,9 @@ struct OnTheFlyTracker { insertHist(histPath + "hRecoMultiplicity", "hRecoMultiplicity;Reco Multiplicity;Counts", {kTH1D, {{axes.axisMultiplicity}}}); if (enablePrimaryVertexing) { - insertHist(histPath + "hDeltaXPVRecoGen", "hDeltaXPVRecoGen;Delta X (reco - gen), cm", {kTH1D, {{axes.axisDeltaVtxCoord}}}); - insertHist(histPath + "hDeltaYPVRecoGen", "hDeltaYPVRecoGen;Delta Y (reco - gen), cm", {kTH1D, {{axes.axisDeltaVtxCoord}}}); - insertHist(histPath + "hDeltaZPVRecoGen", "hDeltaZPVRecoGen;Delta Z (reco - gen), cm", {kTH1D, {{axes.axisDeltaVtxCoord}}}); + insertHist(histPath + "hDeltaXPVRecoGen", "hDeltaXPVRecoGen;Delta X (reco - gen), cm", {kTH2D, {{axes.axisDeltaVtxCoord, axes.axisMultiplicity}}}); + insertHist(histPath + "hDeltaYPVRecoGen", "hDeltaYPVRecoGen;Delta Y (reco - gen), cm", {kTH2D, {{axes.axisDeltaVtxCoord, axes.axisMultiplicity}}}); + insertHist(histPath + "hDeltaZPVRecoGen", "hDeltaZPVRecoGen;Delta Z (reco - gen), cm", {kTH2D, {{axes.axisDeltaVtxCoord, axes.axisMultiplicity}}}); insertHist(histPath + "hDeltaMultPVRecoGen", "hDeltaMultPVRecoGen;Delta Multiplicity (reco - gen)", {kTH1D, {{axes.axisDeltaMultPVRecoGen}}}); insertHist(histPath + "hVtxMultGen", "hVtxMultGen;Generated Vertex Multiplicity", {kTH1D, {{axes.axisVtxMult}}}); insertHist(histPath + "hVtxMultReco", "hVtxMultReco;Reconstructed Vertex Multiplicity", {kTH1D, {{axes.axisVtxMult}}}); @@ -464,6 +496,7 @@ struct OnTheFlyTracker { getHist(TH1, histPath + "hVtxTrials")->GetXaxis()->SetBinLabel(2, "Succeeded"); } + LOG(info) << "Enable secondary smearing for configuration " << icfg << ": " << enableSecondarySmearing.value; if (enableSecondarySmearing) { fastTracker.emplace_back(std::make_unique()); fastTracker[icfg]->SetMagneticField(mMagneticField); @@ -474,6 +507,7 @@ struct OnTheFlyTracker { fastTracker[icfg]->Print(); // print fastTracker settings if (cascadeDecaySettings.doXiQA) { + LOG(info) << "Enabling Xi QA histograms for configuration " << icfg; insertHist(histPath + "hXiBuilding", "hXiBuilding", {kTH1F, {{10, -0.5f, 9.5f}}}); insertHist(histPath + "hGenXi", "hGenXi", {kTH2F, {axes.axisDecayRadius, axes.axisMomentum}}); @@ -523,6 +557,7 @@ struct OnTheFlyTracker { if (doExtraQA) { insertHist(histPath + "h2dPtRes", "h2dPtRes;Gen p_{T};#Delta p_{T} / Reco p_{T}", {kTH2D, {{axes.axisMomentum, axes.axisPtRes}}}); + insertHist(histPath + "h2dPtResAbs", "h2dPtResAbs;Gen p_{T};#Delta p_{T}", {kTH2D, {{axes.axisMomentum, axes.axisPtRes}}}); insertHist(histPath + "h2dDCAxy", "h2dDCAxy;p_{T};DCA_{xy}", {kTH2D, {{axes.axisMomentum, axes.axisDCA}}}); insertHist(histPath + "h2dDCAz", "h2dDCAz;p_{T};DCA_{z}", {kTH2D, {{axes.axisMomentum, axes.axisDCA}}}); } @@ -599,7 +634,6 @@ struct OnTheFlyTracker { insertHist(v0histPath + "K0/hRecoNegDaughterFromV0", "hRecoNegDaughterFromV0", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); insertHist(v0histPath + "K0/hRecoPosDaughterFromV0", "hRecoPosDaughterFromV0", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); insertHist(v0histPath + "K0/hMass", "hMass", kTH2F, {axes.axisK0Mass, axes.axisMomentum}); - insertHist(v0histPath + "K0/hFinalMass", "hFinalMass", kTH1F, {axes.axisK0Mass}); // Lambda insertHist(v0histPath + "Lambda/hGen", "hGen", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); insertHist(v0histPath + "Lambda/hReco", "hReco", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); @@ -608,8 +642,6 @@ struct OnTheFlyTracker { insertHist(v0histPath + "Lambda/hRecoNegDaughterFromV0", "hRecoNegDaughterFromV0", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); insertHist(v0histPath + "Lambda/hRecoPosDaughterFromV0", "hRecoPosDaughterFromV0", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); insertHist(v0histPath + "Lambda/hMass", "hMass", kTH2F, {axes.axisLambdaMass, axes.axisMomentum}); - insertHist(v0histPath + "Lambda/hFinalMass", "hFinalMass", kTH1F, {axes.axisLambdaMass}); - // AntiLambda insertHist(v0histPath + "AntiLambda/hGen", "hGen", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); insertHist(v0histPath + "AntiLambda/hReco", "hReco", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); @@ -618,7 +650,6 @@ struct OnTheFlyTracker { insertHist(v0histPath + "AntiLambda/hRecoNegDaughterFromV0", "hRecoNegDaughterFromV0", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); insertHist(v0histPath + "AntiLambda/hRecoPosDaughterFromV0", "hRecoPosDaughterFromV0", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); insertHist(v0histPath + "AntiLambda/hMass", "hMass", kTH2F, {axes.axisLambdaMass, axes.axisMomentum}); - insertHist(v0histPath + "AntiLambda/hFinalMass", "hFinalMass", kTH1F, {axes.axisLambdaMass}); } } @@ -806,30 +837,12 @@ struct OnTheFlyTracker { decayDaughters.push_back(*v0Decay.GetDecay(1)); } - float dNdEta = 0.f; // Charged particle multiplicity to use in the efficiency evaluation - std::pair vertexReconstructionEfficiencyCounters = {0, 0}; // {nVerticesWithMoreThan2Contributors, nVerticesReconstructed} - void processWithLUTs(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles, const int icfg) - { - vertexReconstructionEfficiencyCounters.first += 1; - const int lastTrackIndex = tableStoredTracksCov.lastIndex() + 1; // bookkeep the last added track - const std::string histPath = "Configuration_" + std::to_string(icfg) + "/"; - - tracksAlice3.clear(); - ghostTracksAlice3.clear(); - bcData.clear(); - cascadesAlice3.clear(); - v0sAlice3.clear(); - - o2::dataformats::DCA dcaInfo; - o2::dataformats::VertexBase vtx; - - // generate collision time - auto ir = irSampler.generateCollisionTime(); - const float eventCollisionTimeNS = ir.timeInBCNS; - - // First we compute the number of charged particles in the event - dNdEta = 0.f; - LOG(debug) << "Processing " << mcParticles.size() << " MC particles to compute dNch/deta"; + /// Function to compute dN/deta for a given set of MC particles + /// \param dNdEta the address of the variable to fill with the computed dN/deta value + /// \param mcParticles the set of MC particles to compute dN/deta from + /// \param histPath the path to the histogram where the computed dN/deta value will be stored for QA purposes + template + void computeDNDEta(float& dNdEta, McParticleType const& mcParticles, const std::string histPath) { for (const auto& mcParticle : mcParticles) { if (std::abs(mcParticle.eta()) > multEtaRange) { continue; @@ -860,852 +873,810 @@ struct OnTheFlyTracker { LOG(debug) << "Computed dNch/deta before normalization: " << dNdEta; dNdEta /= (multEtaRange * 2.0f); - uint32_t multiplicityCounter = 0; getHist(TH1, histPath + "hLUTMultiplicity")->Fill(dNdEta); - // Now that the multiplicity is known, we can process the particles to smear them - double xiDecayRadius2D = 0; - double laDecayRadius2D = 0; - double v0DecayRadius2D = 0; + } + + /// Function to study the cascade decay and fill the relevant histograms and output track vector + /// \param trackTableOffset the offset to apply to the global track indices when filling the output track vector + /// \param mcParticle the cascade MC particle to study + /// \param tracksCascadeProngs the address of the vector of output tracks to fill with the cascade daughters + /// \param primaryVertex the primary vertex of the event, needed for the track propagation in the cascade decay + /// \param icfg the index of the current configuration, needed for histogram filling + /// \param dNdEta the dN/deta of the event, needed for the track time smearing + /// \param eventCollisionTimeNS the collision time of the event in ns, needed for the track time smearing + template + void studyCascade(const int trackTableOffset, + const McParticleType& mcParticle, + std::vector& tracksCascadeProngs, + o2::vertexing::PVertex& primaryVertex, + int icfg, + int dNdEta, + float eventCollisionTimeNS) + { + LOG(info) << ""; + LOG(info) << "-------------------------------"; + LOG(info) << "Studying cascade with PDG " << mcParticle.pdgCode() << " and pT " << mcParticle.pt(); + o2::track::TrackParCov trackParCov; + o2::upgrade::convertMCParticleToO2Track(mcParticle, trackParCov, pdgDB); + const std::string histPath = "Configuration_" + std::to_string(icfg) + "/"; + std::vector cascadeDecayProducts; - std::vector v0DecayProducts; - std::vector xiDecayVertex, laDecayVertex, v0DecayVertex; - for (const auto& mcParticle : mcParticles) { - xiDecayRadius2D = 0; - laDecayRadius2D = 0; - v0DecayRadius2D = 0; - xiDecayVertex.clear(); - laDecayVertex.clear(); - v0DecayVertex.clear(); - - cascadeDecayProducts.clear(); - v0DecayProducts.clear(); - const bool isCascade = mcParticle.pdgCode() == kXiMinus; - if (cascadeDecaySettings.decayXi && isCascade) { - o2::track::TrackParCov xiTrackParCov; - o2::upgrade::convertMCParticleToO2Track(mcParticle, xiTrackParCov, pdgDB); - decayCascade(mcParticle, xiTrackParCov, cascadeDecayProducts, xiDecayVertex, laDecayVertex); - if (cascadeDecayProducts.size() != 3) { - LOG(fatal) << "Xi decay did not produce 3 daughters as expected!"; + std::vector xiDecayVertex, laDecayVertex; + static constexpr int kCascProngs = 3; + std::array xiDaughterTrackParCovsPerfect; + std::array xiDaughterTrackParCovsTracked; + std::array isReco; + std::array nHitsCascadeProngs; // total + std::array nSiliconHitsCascadeProngs; // silicon type + std::array nTPCHitsCascadeProngs; // TPC type + + o2::track::TrackParCov xiTrackParCov; + o2::upgrade::convertMCParticleToO2Track(mcParticle, xiTrackParCov, pdgDB); + LOG(info) << "Decaying cascade..."; + decayCascade(mcParticle, xiTrackParCov, cascadeDecayProducts, xiDecayVertex, laDecayVertex); + if (cascadeDecayProducts.size() != 3) { + LOG(fatal) << "Xi decay did not produce 3 daughters as expected!"; + } + double xiDecayRadius2D = std::hypot(xiDecayVertex[0], xiDecayVertex[1]); + double laDecayRadius2D = std::hypot(laDecayVertex[0], laDecayVertex[1]); + + if (cascadeDecaySettings.doXiQA) { + getHist(TH2, histPath + "hGenXi")->Fill(xiDecayRadius2D, mcParticle.pt()); + getHist(TH2, histPath + "hGenPiFromXi")->Fill(xiDecayRadius2D, cascadeDecayProducts[0].Pt()); + getHist(TH2, histPath + "hGenPiFromLa")->Fill(laDecayRadius2D, cascadeDecayProducts[1].Pt()); + getHist(TH2, histPath + "hGenPrFromLa")->Fill(laDecayRadius2D, cascadeDecayProducts[2].Pt()); + } + + if (cascadeDecaySettings.doXiQA) { + getHist(TH1, histPath + "hXiBuilding")->Fill(0.0f); + } + + o2::upgrade::convertTLorentzVectorToO2Track(PDG_t::kPiMinus, cascadeDecayProducts[0], xiDecayVertex, xiDaughterTrackParCovsPerfect[0], pdgDB); + xiDaughterTrackParCovsPerfect[0].setPID(pdgCodeToPID(PDG_t::kPiMinus)); + o2::upgrade::convertTLorentzVectorToO2Track(PDG_t::kPiMinus, cascadeDecayProducts[1], laDecayVertex, xiDaughterTrackParCovsPerfect[1], pdgDB); + xiDaughterTrackParCovsPerfect[1].setPID(pdgCodeToPID(PDG_t::kPiMinus)); + o2::upgrade::convertTLorentzVectorToO2Track(PDG_t::kProton, cascadeDecayProducts[2], laDecayVertex, xiDaughterTrackParCovsPerfect[2], pdgDB); + xiDaughterTrackParCovsPerfect[2].setPID(pdgCodeToPID(PDG_t::kProton)); + o2::track::TrackParCov perfectCascadeTrack; + o2::upgrade::convertMCParticleToO2Track(mcParticle, perfectCascadeTrack, pdgDB); + perfectCascadeTrack.setPID(pdgCodeToPID(PDG_t::kXiMinus)); // FIXME: not OK for omegas + + float trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; + tracksCascadeProngs.push_back(TrackAlice3{xiDaughterTrackParCovsPerfect[0], mcParticle.globalIndex(), trackTime, timeResolutionUs, true, true, 1, -1, -1, TrackType::kGenCascDaug}); + trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; + tracksCascadeProngs.push_back(TrackAlice3{xiDaughterTrackParCovsPerfect[1], mcParticle.globalIndex(), trackTime, timeResolutionUs, true, true, 2, -1, -1, TrackType::kGenCascDaug}); + trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; + tracksCascadeProngs.push_back(TrackAlice3{xiDaughterTrackParCovsPerfect[2], mcParticle.globalIndex(), trackTime, timeResolutionUs, true, true, 3, -1, -1, TrackType::kGenCascDaug}); + trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; + tracksCascadeProngs.push_back(TrackAlice3{perfectCascadeTrack, mcParticle.globalIndex(), trackTime, timeResolutionUs, true, true, 3, -1, -1, TrackType::kGenCascDaug}); + LOG(info) << "tracksCascadeProngs.size() after adding gen tracks: " << tracksCascadeProngs.size(); + + for (int i = 0; i < kCascProngs; i++) { + isReco[i] = false; + nHitsCascadeProngs[i] = 0; + nSiliconHitsCascadeProngs[i] = 0; + nTPCHitsCascadeProngs[i] = 0; + if (enableSecondarySmearing) { + nHitsCascadeProngs[i] = fastTracker[icfg]->FastTrack(xiDaughterTrackParCovsPerfect[i], xiDaughterTrackParCovsTracked[i], dNdEta); + nSiliconHitsCascadeProngs[i] = fastTracker[icfg]->GetNSiliconPoints(); + nTPCHitsCascadeProngs[i] = fastTracker[icfg]->GetNGasPoints(); + + if (nHitsCascadeProngs[i] < 0 && cascadeDecaySettings.doXiQA) { // QA + getHist(TH1, histPath + "hFastTrackerQA")->Fill(o2::math_utils::abs(nHitsCascadeProngs[i])); } - xiDecayRadius2D = std::hypot(xiDecayVertex[0], xiDecayVertex[1]); - laDecayRadius2D = std::hypot(laDecayVertex[0], laDecayVertex[1]); - } - const bool isV0 = std::find(v0PDGs.begin(), v0PDGs.end(), mcParticle.pdgCode()) != v0PDGs.end(); - if (v0DecaySettings.decayV0 && isV0) { - decayV0Particle(mcParticle, v0DecayProducts, v0DecayVertex, mcParticle.pdgCode()); - if (v0DecayProducts.size() != 2) { - LOG(fatal) << "V0 decay did not produce 2 daughters as expected!"; + if (nSiliconHitsCascadeProngs[i] < fastTrackerSettings.minSiliconHits) { + LOG(info) << "Prong " << i << " failed fast tracker silicon hit requirement with " << nSiliconHitsCascadeProngs[i] << " hits, required: " << fastTrackerSettings.minSiliconHits; + continue; } - v0DecayRadius2D = std::hypot(v0DecayVertex[0], v0DecayVertex[1]); - } + if (nSiliconHitsCascadeProngs[i] < fastTrackerSettings.minSiliconHitsIfTPCUsed) { + LOG(info) << "Prong " << i << " failed fast tracker silicon hit requirement with " << nSiliconHitsCascadeProngs[i] << " hits, required: " << fastTrackerSettings.minSiliconHitsIfTPCUsed; + continue; + } + if (nTPCHitsCascadeProngs[i] < fastTrackerSettings.minTPCClusters) { + LOG(info) << "Prong " << i << " failed fast tracker TPC hit requirement with " << nTPCHitsCascadeProngs[i] << " hits, required: " << fastTrackerSettings.minTPCClusters; + continue; + } + isReco[i] = true; - if (!mcParticle.isPhysicalPrimary()) { - continue; + for (uint32_t ih = 0; ih < fastTracker[icfg]->GetNHits() && cascadeDecaySettings.doXiQA; ih++) { + getHist(TH2, histPath + "hFastTrackerHits")->Fill(fastTracker[icfg]->GetHitZ(ih), + std::hypot(fastTracker[icfg]->GetHitX(ih), fastTracker[icfg]->GetHitY(ih))); + } + } else { + isReco[i] = true; + xiDaughterTrackParCovsTracked[i] = xiDaughterTrackParCovsPerfect[i]; } - const bool longLivedToBeHandled = std::find(longLivedHandledPDGs.begin(), longLivedHandledPDGs.end(), std::abs(mcParticle.pdgCode())) != longLivedHandledPDGs.end(); - const bool nucleiToBeHandled = std::find(nucleiPDGs.begin(), nucleiPDGs.end(), std::abs(mcParticle.pdgCode())) != nucleiPDGs.end(); - const bool pdgsToBeHandled = longLivedToBeHandled || (enableNucleiSmearing && nucleiToBeHandled) || (cascadeDecaySettings.decayXi && isCascade) || (v0DecaySettings.decayV0 && isV0); - if (!pdgsToBeHandled) { + if (TMath::IsNaN(xiDaughterTrackParCovsTracked[i].getZ())) { continue; + } else { + histos.fill(HIST("hNaNBookkeeping"), i + 1, 1.0f); } + trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; + // TODO: add flag for whether it's a ghost track or not, currently assuming all are reconstructed tracks if they pass the fast tracker requirements + TrackType trackType = isReco[i] ? TrackType::kRecoCascDaug : TrackType::kGenCascDaug; + tracksCascadeProngs[i] = TrackAlice3{xiDaughterTrackParCovsTracked[i], mcParticle.globalIndex(), trackTime, timeResolutionUs, true, true, i + 2, nSiliconHitsCascadeProngs[i], nTPCHitsCascadeProngs[i], trackType}; + } - if (std::fabs(mcParticle.eta()) > maxEta) { - continue; - } - if (enablePrimarySmearing) { - getHist(TH1, histPath + "hPtGenerated")->Fill(mcParticle.pt()); - getHist(TH1, histPath + "hPhiGenerated")->Fill(mcParticle.phi()); - switch (std::abs(mcParticle.pdgCode())) { - case kElectron: - getHist(TH1, histPath + "hPtGeneratedEl")->Fill(mcParticle.pt()); - break; - case kPiPlus: - getHist(TH1, histPath + "hPtGeneratedPi")->Fill(mcParticle.pt()); - break; - case kKPlus: - getHist(TH1, histPath + "hPtGeneratedKa")->Fill(mcParticle.pt()); - break; - case kProton: - getHist(TH1, histPath + "hPtGeneratedPr")->Fill(mcParticle.pt()); - break; - } - } - if (cascadeDecaySettings.doXiQA && isCascade) { - getHist(TH2, histPath + "hGenXi")->Fill(xiDecayRadius2D, mcParticle.pt()); - getHist(TH2, histPath + "hGenPiFromXi")->Fill(xiDecayRadius2D, cascadeDecayProducts[0].Pt()); - getHist(TH2, histPath + "hGenPiFromLa")->Fill(laDecayRadius2D, cascadeDecayProducts[1].Pt()); - getHist(TH2, histPath + "hGenPrFromLa")->Fill(laDecayRadius2D, cascadeDecayProducts[2].Pt()); + bool tryKinkReco = false; + if (!isReco[1] || !isReco[2]) { + tryKinkReco = true; // Lambda outside acceptance, set flag for kink reco to be used if mode 1 + } + + bool reconstructedCascade = false; + if (isReco[0] && isReco[1] && isReco[2]) { + LOG(info) << "All cascade daughters reconstructed, attempting cascade reconstruction with DCA finder..."; + reconstructedCascade = true; + } + + // +-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+ + // combine particles into actual Xi candidate + // cascade building starts here + if (cascadeDecaySettings.findXi && reconstructedCascade && cascadeDecaySettings.doKinkReco != 2) { + if (cascadeDecaySettings.doXiQA) { + getHist(TH1, histPath + "hXiBuilding")->Fill(3.0f); } - if (v0DecaySettings.doV0QA && isV0) { - for (size_t indexV0 = 0; indexV0 < v0PDGs.size(); indexV0++) { - if (mcParticle.pdgCode() != v0PDGs[indexV0]) { - continue; - } - for (int indexDetector = 0; indexDetector < mGeoContainer.getNumberOfConfigurations(); indexDetector++) { - std::string path = Form("V0Building_Configuration_%i/%s/", indexDetector, NameV0s[indexV0].data()); - fillHist(TH2, path + "hGen", v0DecayRadius2D, mcParticle.pt()); - fillHist(TH2, path + "hGenNegDaughterFromV0", v0DecayRadius2D, v0DecayProducts[0].Pt()); - fillHist(TH2, path + "hGenPosDaughterFromV0", v0DecayRadius2D, v0DecayProducts[1].Pt()); - } - } + // 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.bachelorId = trackTableOffset + 1; + thisCascade.positiveId = trackTableOffset + 2; // Lambda daughters + thisCascade.negativeId = trackTableOffset + 3; // Lambda daughters + thisCascade.cascadeTrackId = trackTableOffset + 4; + + // use DCA fitters + int nCand = 0; + bool dcaFitterOK_V0 = true; + try { + nCand = fitter.process(xiDaughterTrackParCovsTracked[1], xiDaughterTrackParCovsTracked[2]); + } catch (...) { + // LOG(error) << "Exception caught in DCA fitter process call!"; + dcaFitterOK_V0 = false; } - if (mcParticle.pt() < minPt) { - continue; + if (nCand == 0) { + dcaFitterOK_V0 = false; } - o2::track::TrackParCov trackParCov; - o2::upgrade::convertMCParticleToO2Track(mcParticle, trackParCov, pdgDB); - - const bool isDecayDaughter = (mcParticle.getProcess() == TMCProcess::kPDecay); + fitter.propagateTracksToVertex(); + if (!fitter.isPropagateTracksToVertexDone()) { + dcaFitterOK_V0 = false; + } - multiplicityCounter++; - const float timeResolutionNs = 100.f; // ns - const float nsToMus = 1e-3f; - const float timeResolutionUs = timeResolutionNs * nsToMus; // us - const float time = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; - static constexpr int kCascProngs = 3; - std::array xiDaughterTrackParCovsPerfect; - std::array xiDaughterTrackParCovsTracked; - std::array isReco; - std::array nHitsCascadeProngs; // total - std::array nSiliconHitsCascadeProngs; // silicon type - std::array nTPCHitsCascadeProngs; // TPC type - - bool tryKinkReco = false; - if (cascadeDecaySettings.decayXi && isCascade) { - bool reconstructedCascade = false; + // V0 found successfully + if (dcaFitterOK_V0) { + LOG(info) << "V0 candidate found successfully with DCA fitter, now trying to combine with bachelor track to form cascade..."; if (cascadeDecaySettings.doXiQA) { - getHist(TH1, histPath + "hXiBuilding")->Fill(0.0f); + getHist(TH1, histPath + "hXiBuilding")->Fill(4.0f); } - o2::upgrade::convertTLorentzVectorToO2Track(PDG_t::kPiMinus, cascadeDecayProducts[0], xiDecayVertex, xiDaughterTrackParCovsPerfect[0], pdgDB); - xiDaughterTrackParCovsPerfect[0].setPID(pdgCodeToPID(PDG_t::kPiMinus)); - o2::upgrade::convertTLorentzVectorToO2Track(PDG_t::kPiMinus, cascadeDecayProducts[1], laDecayVertex, xiDaughterTrackParCovsPerfect[1], pdgDB); - xiDaughterTrackParCovsPerfect[1].setPID(pdgCodeToPID(PDG_t::kPiMinus)); - o2::upgrade::convertTLorentzVectorToO2Track(PDG_t::kProton, cascadeDecayProducts[2], laDecayVertex, xiDaughterTrackParCovsPerfect[2], pdgDB); - xiDaughterTrackParCovsPerfect[2].setPID(pdgCodeToPID(PDG_t::kProton)); - - for (int i = 0; i < kCascProngs; i++) { - isReco[i] = false; - nHitsCascadeProngs[i] = 0; - nSiliconHitsCascadeProngs[i] = 0; - nTPCHitsCascadeProngs[i] = 0; - if (enableSecondarySmearing) { - nHitsCascadeProngs[i] = fastTracker[icfg]->FastTrack(xiDaughterTrackParCovsPerfect[i], xiDaughterTrackParCovsTracked[i], dNdEta); - nSiliconHitsCascadeProngs[i] = fastTracker[icfg]->GetNSiliconPoints(); - nTPCHitsCascadeProngs[i] = fastTracker[icfg]->GetNGasPoints(); - - if (nHitsCascadeProngs[i] < 0 && cascadeDecaySettings.doXiQA) { // QA - getHist(TH1, histPath + "hFastTrackerQA")->Fill(o2::math_utils::abs(nHitsCascadeProngs[i])); - } - - if (nSiliconHitsCascadeProngs[i] >= fastTrackerSettings.minSiliconHits || (nSiliconHitsCascadeProngs[i] >= fastTrackerSettings.minSiliconHitsIfTPCUsed && nTPCHitsCascadeProngs[i] >= fastTrackerSettings.minTPCClusters)) { - isReco[i] = true; - } else { - continue; // extra sure - } - for (uint32_t ih = 0; ih < fastTracker[icfg]->GetNHits() && cascadeDecaySettings.doXiQA; ih++) { - getHist(TH2, histPath + "hFastTrackerHits")->Fill(fastTracker[icfg]->GetHitZ(ih), std::hypot(fastTracker[icfg]->GetHitX(ih), fastTracker[icfg]->GetHitY(ih))); - } - } else { - isReco[i] = true; - xiDaughterTrackParCovsTracked[i] = xiDaughterTrackParCovsPerfect[i]; - } + 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]; + } - if (TMath::IsNaN(xiDaughterTrackParCovsTracked[i].getZ())) { - continue; - } else { - histos.fill(HIST("hNaNBookkeeping"), i + 1, 1.0f); - } - if (isReco[i]) { - tracksAlice3.push_back(TrackAlice3{xiDaughterTrackParCovsTracked[i], mcParticle.globalIndex(), time, timeResolutionUs, true, true, i + 2, nSiliconHitsCascadeProngs[i], nTPCHitsCascadeProngs[i]}); - } else { - ghostTracksAlice3.push_back(TrackAlice3{xiDaughterTrackParCovsTracked[i], mcParticle.globalIndex(), time, timeResolutionUs, true, true, i + 2}); - } + // 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(std::array{std::array{posP[0], posP[1], posP[2]}, + std::array{negP[0], negP[1], negP[2]}}, + std::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; } - if (!isReco[1] || !isReco[2]) { - tryKinkReco = true; // Lambda outside acceptance, set flag for kink reco to be used if mode 1 + 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(pdgCodeToPID(PDG_t::kLambda0)); + + // dca fitter step + nCand = 0; + bool dcaFitterOK_Cascade = true; + try { + nCand = fitter.process(v0Track, xiDaughterTrackParCovsTracked[0]); + } catch (...) { + // LOG(error) << "Exception caught in DCA fitter process call!"; + dcaFitterOK_Cascade = false; + } + if (nCand == 0) { + dcaFitterOK_Cascade = false; } - if (isReco[0] && isReco[1] && isReco[2]) { - reconstructedCascade = true; + fitter.propagateTracksToVertex(); + if (!fitter.isPropagateTracksToVertexDone()) { + dcaFitterOK_Cascade = false; } - // +-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+ - // combine particles into actual Xi candidate - // cascade building starts here - if (cascadeDecaySettings.findXi && isReco[0] && isReco[1] && isReco[2] && cascadeDecaySettings.doKinkReco != 2) { + const u_int8_t fitterStatusCode = fitter.getFitStatus(); + histos.fill(HIST("hFitterStatusCode"), fitterStatusCode); + // Cascade found successfully + if (dcaFitterOK_Cascade) { if (cascadeDecaySettings.doXiQA) { - getHist(TH1, histPath + "hXiBuilding")->Fill(3.0f); - } - // 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 = lastTrackIndex + tracksAlice3.size() - 1; - thisCascade.negativeId = lastTrackIndex + tracksAlice3.size() - 2; - thisCascade.bachelorId = lastTrackIndex + tracksAlice3.size() - 3; - - // use DCA fitters - int nCand = 0; - bool dcaFitterOK_V0 = true; - try { - nCand = fitter.process(xiDaughterTrackParCovsTracked[1], xiDaughterTrackParCovsTracked[2]); - } catch (...) { - // LOG(error) << "Exception caught in DCA fitter process call!"; - dcaFitterOK_V0 = false; - } - if (nCand == 0) { - dcaFitterOK_V0 = false; + getHist(TH1, histPath + "hXiBuilding")->Fill(5.0f); } - fitter.propagateTracksToVertex(); - if (!fitter.isPropagateTracksToVertexDone()) { - dcaFitterOK_V0 = false; + o2::track::TrackParCov bachelorTrackAtPCA = fitter.getTrack(1); + const auto& vtxCascade = fitter.getPCACandidate(); + for (int i = 0; i < 3; i++) { + posCascade[i] = vtxCascade[i]; } - // V0 found successfully - if (dcaFitterOK_V0) { - if (cascadeDecaySettings.doXiQA) { - getHist(TH1, histPath + "hXiBuilding")->Fill(4.0f); - } - - 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]; - } + // 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(std::array{std::array{bachP[0], bachP[1], bachP[2]}, + std::array{posP[0] + negP[0], posP[1] + negP[1], posP[2] + negP[2]}}, + std::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(pdgCodeToPID(PDG_t::kXiMinus)); // FIXME: not OK for omegas + + thisCascade.cascradiusMC = xiDecayRadius2D; + thisCascade.pt = cascadeTrack.getPt(); + thisCascade.eta = cascadeTrack.getEta(); + thisCascade.findableClusters = 0; + thisCascade.foundClusters = 0; + + if (cascadeDecaySettings.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 = fastTracker[icfg]->GetLayers().size() - 1; i >= 0; --i) { + o2::fastsim::DetLayer layer = fastTracker[icfg]->GetLayer(i); + if (layer.isInert()) { + continue; // Not an active tracking layer + } - // 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(std::array{std::array{posP[0], posP[1], posP[2]}, - std::array{negP[0], negP[1], negP[2]}}, - std::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; - } + if (thisCascade.cascradiusMC < layer.getRadius()) { + continue; // Cascade did not reach this layer + } - 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(pdgCodeToPID(PDG_t::kLambda0)); - - // dca fitter step - nCand = 0; - bool dcaFitterOK_Cascade = true; - try { - nCand = fitter.process(v0Track, xiDaughterTrackParCovsTracked[0]); - } catch (...) { - // LOG(error) << "Exception caught in DCA fitter process call!"; - dcaFitterOK_Cascade = false; - } - if (nCand == 0) { - dcaFitterOK_Cascade = false; - } + // cascade decayed after the corresponding radius + thisCascade.findableClusters++; // add to findable - fitter.propagateTracksToVertex(); - if (!fitter.isPropagateTracksToVertexDone()) { - dcaFitterOK_Cascade = false; - } + // find perfect intercept XYZ + float targetX = 1e+3; + trackParCov.getXatLabR(layer.getRadius(), targetX, mMagneticField); + if (targetX > 999) { + continue; // failed to find intercept + } - const u_int8_t fitterStatusCode = fitter.getFitStatus(); - histos.fill(HIST("hFitterStatusCode"), fitterStatusCode); - // Cascade found successfully - if (dcaFitterOK_Cascade) { - if (cascadeDecaySettings.doXiQA) { - getHist(TH1, histPath + "hXiBuilding")->Fill(5.0f); + if (!trackParCov.propagateTo(targetX, mMagneticField)) { + continue; // failed to propagate } - o2::track::TrackParCov bachelorTrackAtPCA = fitter.getTrack(1); - const auto& vtxCascade = fitter.getPCACandidate(); - for (int i = 0; i < 3; i++) { - posCascade[i] = vtxCascade[i]; + // 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::constants::math::PI}; + + if (layer.getResolutionRPhi() > 1e-8 && layer.getResolutionZ() > 1e-8) { // catch zero (though should not really happen...) + phi = gRandom->Gaus(phi, std::asin(layer.getResolutionRPhi() / r)); + posClusterCandidate[0] = r * std::cos(phi); + posClusterCandidate[1] = r * std::sin(phi); + posClusterCandidate[2] = gRandom->Gaus(posClusterCandidate[2], layer.getResolutionZ()); } - // 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(std::array{std::array{bachP[0], bachP[1], bachP[2]}, - std::array{posP[0] + negP[0], posP[1] + negP[1], posP[2] + negP[2]}}, - std::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(pdgCodeToPID(PDG_t::kXiMinus)); // FIXME: not OK for omegas - - thisCascade.cascradiusMC = xiDecayRadius2D; - thisCascade.pt = cascadeTrack.getPt(); - thisCascade.eta = cascadeTrack.getEta(); - thisCascade.findableClusters = 0; - thisCascade.foundClusters = 0; - - if (cascadeDecaySettings.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 = fastTracker[icfg]->GetLayers().size() - 1; i >= 0; --i) { - o2::fastsim::DetLayer layer = fastTracker[icfg]->GetLayer(i); - if (layer.isInert()) { - continue; // Not an active tracking layer - } - - if (thisCascade.cascradiusMC < layer.getRadius()) { - continue; // Cascade did not reach this layer - } - - // cascade decayed after the corresponding radius - thisCascade.findableClusters++; // add to findable - - // find perfect intercept XYZ - float targetX = 1e+3; - trackParCov.getXatLabR(layer.getRadius(), targetX, mMagneticField); - if (targetX > 999) { - continue; // failed to find intercept - } - - if (!trackParCov.propagateTo(targetX, mMagneticField)) { - 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::constants::math::PI}; - - if (layer.getResolutionRPhi() > 1e-8 && layer.getResolutionZ() > 1e-8) { // catch zero (though should not really happen...) - phi = gRandom->Gaus(phi, std::asin(layer.getResolutionRPhi() / r)); - posClusterCandidate[0] = r * std::cos(phi); - posClusterCandidate[1] = r * std::sin(phi); - posClusterCandidate[2] = gRandom->Gaus(posClusterCandidate[2], layer.getResolutionZ()); - } - - if (std::isnan(phi)) { - continue; // Catch when getXatLabR misses layer[i] - } - - // 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], mMagneticField))) { - continue; - } - - const o2::track::TrackParametrization::dim2_t hitpoint = {static_cast(xyz1[1]), static_cast(xyz1[2])}; - const o2::track::TrackParametrization::dim3_t hitpointcov = {layer.getResolutionRPhi() * layer.getResolutionRPhi(), 0.f, layer.getResolutionZ() * layer.getResolutionZ()}; - if (layer.isInDeadPhiRegion(phi)) { - continue; // No hit for strangeness tracking update - } - - cascadeTrack.update(hitpoint, hitpointcov); - thisCascade.foundClusters++; // add to findable - } - - if (thisCascade.foundClusters < cascadeDecaySettings.minStraTrackHits) { - continue; // We didn't find enough hits for strangeness tracking - } + if (std::isnan(phi)) { + continue; // Catch when getXatLabR misses layer[i] } - // add cascade track - thisCascade.cascadeTrackId = lastTrackIndex + tracksAlice3.size(); // this is the next index to be filled -> should be it - tracksAlice3.push_back(TrackAlice3{cascadeTrack, mcParticle.globalIndex(), time, timeResolutionUs, false, false, 1, thisCascade.foundClusters}); + // 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]}; - // add this cascade to vector (will fill cursor later with collision ID) - cascadesAlice3.push_back(thisCascade); - } - } - } // end cascade building + if (!(cascadeTrack.propagateTo(xyz1[0], mMagneticField))) { + continue; + } - if (isReco[0] && ((cascadeDecaySettings.doKinkReco == 1 && tryKinkReco) || cascadeDecaySettings.doKinkReco == 2)) { // mode 1 or 2 - o2::track::TrackParCov prefectCascadeTrack, trackedCascade; - const o2::track::TrackParCov& trackedBach = xiDaughterTrackParCovsTracked[0]; - o2::upgrade::convertMCParticleToO2Track(mcParticle, prefectCascadeTrack, pdgDB); + const o2::track::TrackParametrization::dim2_t hitpoint = {static_cast(xyz1[1]), static_cast(xyz1[2])}; + const o2::track::TrackParametrization::dim3_t hitpointcov = {layer.getResolutionRPhi() * layer.getResolutionRPhi(), 0.f, layer.getResolutionZ() * layer.getResolutionZ()}; + if (layer.isInDeadPhiRegion(phi)) { + continue; // No hit for strangeness tracking update + } - // back track is already smeared - prefectCascadeTrack.setPID(pdgCodeToPID(PDG_t::kXiMinus)); // FIXME: not OK for omegas - const int nCascHits = fastTracker[icfg]->FastTrack(prefectCascadeTrack, trackedCascade, dNdEta, xiDecayRadius2D); - reconstructedCascade = (fastTrackerSettings.minSiliconHitsForKinkReco < nCascHits) ? true : false; - if (reconstructedCascade) { - std::array pCasc; - std::array pBach; - std::array pV0; - trackedCascade.getPxPyPzGlo(pCasc); - trackedBach.getPxPyPzGlo(pBach); - for (size_t i = 0; i < pCasc.size(); ++i) { - pV0[i] = pCasc[i] - pBach[i]; + cascadeTrack.update(hitpoint, hitpointcov); + thisCascade.foundClusters++; // add to findable } - if (isReco[1] && !isReco[2]) { - thisCascade.negativeId = lastTrackIndex + tracksAlice3.size() - 1; - thisCascade.positiveId = -1; - } else if (!isReco[1] && isReco[2]) { - thisCascade.negativeId = -1; - thisCascade.positiveId = lastTrackIndex + tracksAlice3.size() - 1; - } else if (isReco[1] && isReco[2]) { - thisCascade.positiveId = lastTrackIndex + tracksAlice3.size() - 1; - thisCascade.negativeId = lastTrackIndex + tracksAlice3.size() - 2; - } else { - thisCascade.positiveId = -1; - thisCascade.negativeId = -1; + if (thisCascade.foundClusters < cascadeDecaySettings.minStraTrackHits) { + return; // We didn't find enough hits for strangeness tracking } + } + tracksCascadeProngs[kCascProngs + 1] = TrackAlice3{cascadeTrack, mcParticle.globalIndex(), trackTime, timeResolutionUs, false, false, 1, thisCascade.foundClusters, TrackType::kGenCascDaug}; + } + } + } // end cascade building + + if (isReco[0] && ((cascadeDecaySettings.doKinkReco == 1 && tryKinkReco) || cascadeDecaySettings.doKinkReco == 2)) { // mode 1 or 2 + o2::track::TrackParCov trackedCascade; + const o2::track::TrackParCov& trackedBach = xiDaughterTrackParCovsTracked[0]; + const int nCascHits = fastTracker[icfg]->FastTrack(perfectCascadeTrack, trackedCascade, dNdEta, xiDecayRadius2D); + reconstructedCascade = (fastTrackerSettings.minSiliconHitsForKinkReco < nCascHits) ? true : false; + if (reconstructedCascade) { + std::array pCasc; + std::array pBach; + std::array pV0; + trackedCascade.getPxPyPzGlo(pCasc); + trackedBach.getPxPyPzGlo(pBach); + for (size_t i = 0; i < pCasc.size(); ++i) { + pV0[i] = pCasc[i] - pBach[i]; + } - int nCand = 0; - bool kinkFitterOK = true; - try { - nCand = fitter.process(trackedCascade, trackedBach); - } catch (...) { - kinkFitterOK = false; - } + // Reset indices + if (!isReco[1]) { thisCascade.negativeId = -1; } + if (!isReco[2]) { thisCascade.positiveId = -1; } - if (nCand == 0) { - kinkFitterOK = false; - } + int nCand = 0; + bool kinkFitterOK = true; + try { + nCand = fitter.process(trackedCascade, trackedBach); + } catch (...) { + kinkFitterOK = false; + } - fitter.propagateTracksToVertex(); - if (!fitter.isPropagateTracksToVertexDone()) { - kinkFitterOK = false; - } + if (nCand == 0) { + kinkFitterOK = false; + } - const u_int8_t fitterStatusCode = fitter.getFitStatus(); - histos.fill(HIST("hFitterStatusCode"), fitterStatusCode); - if (kinkFitterOK) { - if (cascadeDecaySettings.doXiQA) { - getHist(TH1, histPath + "hXiBuilding")->Fill(6.0f); - } + fitter.propagateTracksToVertex(); + if (!fitter.isPropagateTracksToVertexDone()) { + kinkFitterOK = false; + } - o2::track::TrackParCov newCascadeTrack = fitter.getTrack(0); // (cascade) - std::array kinkVtx = {-999, -999, -999}; - kinkVtx = fitter.getPCACandidatePos(); - thisCascade.bachelorId = lastTrackIndex + tracksAlice3.size() - isReco.size(); - thisCascade.cascadeTrackId = lastTrackIndex + tracksAlice3.size(); // this should be ok - thisCascade.dcaV0dau = -1.f; // unknown - thisCascade.v0radius = -1.f; // unknown - thisCascade.dcacascdau = std::sqrt(fitter.getChi2AtPCACandidate()); - thisCascade.cascradius = std::hypot(kinkVtx[0], kinkVtx[1]); - thisCascade.cascradiusMC = xiDecayRadius2D; - thisCascade.mLambda = o2::constants::physics::MassLambda; - thisCascade.findableClusters = nCascHits; - thisCascade.foundClusters = nCascHits; - thisCascade.pt = newCascadeTrack.getPt(); - thisCascade.eta = newCascadeTrack.getEta(); - thisCascade.mXi = RecoDecay::m(std::array{std::array{pBach[0], pBach[1], pBach[2]}, std::array{pV0[0], pV0[1], pV0[2]}}, - std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassLambda}); - - newCascadeTrack.setPID(pdgCodeToPID(PDG_t::kXiMinus)); // FIXME: not OK for omegas - tracksAlice3.push_back(TrackAlice3{newCascadeTrack, mcParticle.globalIndex(), time, timeResolutionUs, false, false, 1, thisCascade.foundClusters}); - - // add this cascade to vector (will fill cursor later with collision ID) - cascadesAlice3.push_back(thisCascade); - } // end fitter OK - } // end cascade found - } // end cascade kink building - - // +-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+ - if (cascadeDecaySettings.doXiQA) { - if (reconstructedCascade) { - getHist(TH2, histPath + "hRecoXi")->Fill(xiDecayRadius2D, mcParticle.pt()); - getHist(TH1, histPath + "hMassLambda")->Fill(thisCascade.mLambda); - getHist(TH1, histPath + "hMassXi")->Fill(thisCascade.mXi); - getHist(TH2, histPath + "h2dMassXi")->Fill(thisCascade.mXi, thisCascade.pt); - getHist(TH2, histPath + "h2dDeltaPtVsPt")->Fill(thisCascade.pt, (mcParticle.pt() - thisCascade.pt) / thisCascade.pt); - getHist(TH2, histPath + "h2dDeltaEtaVsPt")->Fill(thisCascade.pt, mcParticle.eta() - thisCascade.eta); - getHist(TH2, histPath + "hFoundVsFindable")->Fill(thisCascade.findableClusters, thisCascade.foundClusters); - } - if (isReco[0]) { - getHist(TH2, histPath + "hRecoPiFromXi")->Fill(xiDecayRadius2D, cascadeDecayProducts[0].Pt()); - } - if (isReco[1]) { - getHist(TH2, histPath + "hRecoPiFromLa")->Fill(laDecayRadius2D, cascadeDecayProducts[1].Pt()); + const u_int8_t fitterStatusCode = fitter.getFitStatus(); + histos.fill(HIST("hFitterStatusCode"), fitterStatusCode); + if (kinkFitterOK) { + if (cascadeDecaySettings.doXiQA) { + getHist(TH1, histPath + "hXiBuilding")->Fill(6.0f); } - if (isReco[2]) { - getHist(TH2, histPath + "hRecoPrFromLa")->Fill(laDecayRadius2D, cascadeDecayProducts[2].Pt()); + + o2::track::TrackParCov newCascadeTrack = fitter.getTrack(0); // (cascade) + std::array kinkVtx = {-999, -999, -999}; + kinkVtx = fitter.getPCACandidatePos(); + thisCascade.dcaV0dau = -1.f; // unknown + thisCascade.v0radius = -1.f; // unknown + thisCascade.dcacascdau = std::sqrt(fitter.getChi2AtPCACandidate()); + thisCascade.cascradius = std::hypot(kinkVtx[0], kinkVtx[1]); + thisCascade.cascradiusMC = xiDecayRadius2D; + thisCascade.mLambda = o2::constants::physics::MassLambda; + thisCascade.findableClusters = nCascHits; + thisCascade.foundClusters = nCascHits; + thisCascade.pt = newCascadeTrack.getPt(); + thisCascade.eta = newCascadeTrack.getEta(); + thisCascade.mXi = RecoDecay::m(std::array{std::array{pBach[0], pBach[1], pBach[2]}, std::array{pV0[0], pV0[1], pV0[2]}}, + std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassLambda}); + + newCascadeTrack.setPID(pdgCodeToPID(PDG_t::kXiMinus)); // FIXME: not OK for omegas + float trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; + if (reconstructedCascade) { + tracksCascadeProngs[kCascProngs + 1] = TrackAlice3{newCascadeTrack, mcParticle.globalIndex(), trackTime, timeResolutionUs, false, false, 1, thisCascade.foundClusters, TrackType::kRecoCascDaug}; } - } + } // end fitter OK + } // end cascade found + } // end cascade kink building + + // +-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+ + if (cascadeDecaySettings.doXiQA) { + double dcaXY{-1.}, dcaZ{-1.}; + if (reconstructedCascade) { + getHist(TH2, histPath + "hRecoXi")->Fill(xiDecayRadius2D, mcParticle.pt()); + getHist(TH1, histPath + "hMassLambda")->Fill(thisCascade.mLambda); + getHist(TH1, histPath + "hMassXi")->Fill(thisCascade.mXi); + getHist(TH2, histPath + "h2dMassXi")->Fill(thisCascade.mXi, thisCascade.pt); + getHist(TH2, histPath + "h2dDeltaPtVsPt")->Fill(thisCascade.pt, (mcParticle.pt() - thisCascade.pt) / thisCascade.pt); + getHist(TH2, histPath + "h2dDeltaEtaVsPt")->Fill(thisCascade.pt, mcParticle.eta() - thisCascade.eta); + getHist(TH2, histPath + "hFoundVsFindable")->Fill(thisCascade.findableClusters, thisCascade.foundClusters); - continue; // Cascade handling done, should not be considered anymore + o2::track::TrackParCov trackParametrization(trackParCov); + trackParametrization.propagateToDCA(primaryVertex, mMagneticField, &dcaInfo); + getHist(TH2, histPath + "h2dDCAxyCascade")->Fill(trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please + getHist(TH2, histPath + "h2dDCAzCascade")->Fill(trackParametrization.getPt(), dcaZ * 1e+4); // in microns, please + } + if (isReco[0]) { + getHist(TH2, histPath + "hRecoPiFromXi")->Fill(xiDecayRadius2D, cascadeDecayProducts[0].Pt()); + o2::track::TrackParCov trackParametrizationCascProng0(tracksCascadeProngs[0]); + if (populateTracksDCA && tracksCascadeProngs[0].propagateToDCA(primaryVertex, mMagneticField, &dcaInfo)) { // FIXME: this is not the right trackParametrization, need to propagate the bachelor track + dcaXY = dcaInfo.getY(); + dcaZ = dcaInfo.getZ(); + getHist(TH2, histPath + "h2dDCAxyCascadeBachelor")->Fill(trackParametrizationCascProng0.getPt(), dcaXY * 1e+4); // in microns, please + getHist(TH2, histPath + "h2dDCAzCascadeBachelor")->Fill(trackParametrizationCascProng0.getPt(), dcaZ * 1e+4); // in microns, please + } + } + if (isReco[1]) { + getHist(TH2, histPath + "hRecoPiFromLa")->Fill(laDecayRadius2D, cascadeDecayProducts[1].Pt()); + o2::track::TrackParCov trackParametrizationCascProng1(tracksCascadeProngs[1]); + if (populateTracksDCA && tracksCascadeProngs[1].propagateToDCA(primaryVertex, mMagneticField, &dcaInfo)) { // FIXME: this is not the right trackParametrization, need to propagate the negative pion track + dcaXY = dcaInfo.getY(); + dcaZ = dcaInfo.getZ(); + getHist(TH2, histPath + "h2dDCAxyCascadeNegative")->Fill(trackParametrizationCascProng1.getPt(), dcaXY * 1e+4); // in microns, please + getHist(TH2, histPath + "h2dDCAzCascadeNegative")->Fill(trackParametrizationCascProng1.getPt(), dcaZ * 1e+4); // in microns, please + } + } + if (isReco[2]) { + getHist(TH2, histPath + "hRecoPrFromLa")->Fill(laDecayRadius2D, cascadeDecayProducts[2].Pt()); + o2::track::TrackParCov trackParametrizationCascProng2(tracksCascadeProngs[2]); + if (populateTracksDCA && tracksCascadeProngs[2].propagateToDCA(primaryVertex, mMagneticField, &dcaInfo)) { // FIXME: this is not the right trackParametrization, need to propagate the positive proton track + dcaXY = dcaInfo.getY(); + dcaZ = dcaInfo.getZ(); + getHist(TH2, histPath + "h2dDCAxyCascadePositive")->Fill(trackParametrizationCascProng2.getPt(), dcaXY * 1e+4); // in microns, please + getHist(TH2, histPath + "h2dDCAzCascadePositive")->Fill(trackParametrizationCascProng2.getPt(), dcaZ * 1e+4); // in microns, please + } } + } - // V0 handling - std::vector v0DaughterTrackParCovsPerfect(2); - std::vector v0DaughterTrackParCovsTracked(2); - std::vector isV0Reco(kv0Prongs); - std::vector nV0Hits(kv0Prongs); // total - std::vector nV0SiliconHits(kv0Prongs); // silicon type - std::vector nV0TPCHits(kv0Prongs); // TPC type - if (v0DecaySettings.decayV0 && isV0) { - if (v0DecaySettings.doV0QA) { - fillHist(TH1, Form("V0Building_Configuration_%i/hV0Building", icfg), 0.0f); + // populate Cascades + tableUpgradeCascades(tableCollisions.lastIndex(), + 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); + } + + /// Function to study V0s and fill the relevant histograms + /// \param trackTableOffset offset to be added to the indices of the daughter tracks when filling the V0 table + /// \param mcParticle the MC particle corresponding to the V0 to be studied + /// \param tracksV0Daugs vector of tracks to which the V0 daughters will be added + /// \param icfg index of the current configuration, used for histogram filling + /// \param dNdEta multiplicity of the event, used for fast tracker smearing + /// \param eventCollisionTimeNS collision time of the event in nanoseconds, used for track time smearing + template + void studyV0(const int trackTableOffset, + const McParticleType& mcParticle, + std::vector& tracksV0Daugs, + int icfg, + int dNdEta, + float eventCollisionTimeNS) { + + o2::track::TrackParCov trackParCov; + o2::upgrade::convertMCParticleToO2Track(mcParticle, trackParCov, pdgDB); + const std::string histPath = "Configuration_" + std::to_string(icfg) + "/"; + + std::vector v0DecayProducts; + std::vector laDecayVertex, v0DecayVertex; + decayV0Particle(mcParticle, v0DecayProducts, v0DecayVertex, mcParticle.pdgCode()); + if (v0DecayProducts.size() != 2) { + LOG(fatal) << "V0 decay did not produce 2 daughters as expected!"; + } + double v0DecayRadius2D = std::hypot(v0DecayVertex[0], v0DecayVertex[1]); + + if (v0DecaySettings.doV0QA) { + for (size_t indexV0 = 0; indexV0 < v0PDGs.size(); indexV0++) { + if (mcParticle.pdgCode() != v0PDGs[indexV0]) { + continue; } - switch (mcParticle.pdgCode()) { - case kK0Short: - o2::upgrade::convertTLorentzVectorToO2Track(kPiMinus, v0DecayProducts[0], v0DecayVertex, v0DaughterTrackParCovsPerfect[0], pdgDB); - o2::upgrade::convertTLorentzVectorToO2Track(kPiPlus, v0DecayProducts[1], v0DecayVertex, v0DaughterTrackParCovsPerfect[1], pdgDB); - break; - case kLambda0: - o2::upgrade::convertTLorentzVectorToO2Track(kPiMinus, v0DecayProducts[0], v0DecayVertex, v0DaughterTrackParCovsPerfect[0], pdgDB); - o2::upgrade::convertTLorentzVectorToO2Track(kProton, v0DecayProducts[1], v0DecayVertex, v0DaughterTrackParCovsPerfect[1], pdgDB); - break; - case kLambda0Bar: - o2::upgrade::convertTLorentzVectorToO2Track(kProtonBar, v0DecayProducts[0], v0DecayVertex, v0DaughterTrackParCovsPerfect[0], pdgDB); - o2::upgrade::convertTLorentzVectorToO2Track(kPiPlus, v0DecayProducts[1], v0DecayVertex, v0DaughterTrackParCovsPerfect[1], pdgDB); - break; - default: - LOG(fatal) << "Unhandled V0 PDG code " << mcParticle.pdgCode(); + for (int indexDetector = 0; indexDetector < mGeoContainer.getNumberOfConfigurations(); indexDetector++) { + std::string path = Form("V0Building_Configuration_%i/%s/", indexDetector, NameV0s[indexV0].data()); + fillHist(TH2, path + "hGen", v0DecayRadius2D, mcParticle.pt()); + fillHist(TH2, path + "hGenNegDaughterFromV0", v0DecayRadius2D, v0DecayProducts[0].Pt()); + fillHist(TH2, path + "hGenPosDaughterFromV0", v0DecayRadius2D, v0DecayProducts[1].Pt()); } - for (int i = 0; i < kv0Prongs; i++) { - isV0Reco[i] = false; - nV0Hits[i] = 0; - nV0SiliconHits[i] = 0; - nV0TPCHits[i] = 0; - if (enableSecondarySmearing) { - nV0Hits[i] = fastTracker[icfg]->FastTrack(v0DaughterTrackParCovsPerfect[i], v0DaughterTrackParCovsTracked[i], dNdEta); - nV0SiliconHits[i] = fastTracker[icfg]->GetNSiliconPoints(); - nV0TPCHits[i] = fastTracker[icfg]->GetNGasPoints(); - - if (nV0Hits[i] < 0 && v0DecaySettings.doV0QA) { // QA - fillHist(TH1, Form("V0Building_Configuration_%i/hFastTrackerQA", icfg), o2::math_utils::abs(nV0Hits[i])); - } + } + } - if (nV0SiliconHits[i] >= fastTrackerSettings.minSiliconHits || (nV0SiliconHits[i] >= fastTrackerSettings.minSiliconHitsIfTPCUsed && nV0TPCHits[i] >= fastTrackerSettings.minTPCClusters)) { - isReco[i] = true; - } else { - continue; // extra sure - } - for (uint32_t ih = 0; ih < fastTracker[icfg]->GetNHits() && v0DecaySettings.doV0QA; ih++) { - fillHist(TH2, Form("V0Building_Configuration_%i/hFastTrackerHits", icfg), fastTracker[icfg]->GetHitZ(ih), std::hypot(fastTracker[icfg]->GetHitX(ih), fastTracker[icfg]->GetHitY(ih))); - } - } else { - isReco[i] = true; - v0DaughterTrackParCovsTracked[i] = v0DaughterTrackParCovsPerfect[i]; - } + // V0 handling + std::vector v0DaughterTrackParCovsPerfect(2); + std::vector v0DaughterTrackParCovsTracked(2); + std::vector isV0Reco(kv0Prongs); + std::vector nV0Hits(kv0Prongs); // total + std::vector nV0SiliconHits(kv0Prongs); // silicon type + std::vector nV0TPCHits(kv0Prongs); // TPC type + if (v0DecaySettings.doV0QA) { + fillHist(TH1, Form("V0Building_Configuration_%i/hV0Building", icfg), 0.0f); + } + switch (mcParticle.pdgCode()) { + case kK0Short: + o2::upgrade::convertTLorentzVectorToO2Track(kPiMinus, v0DecayProducts[0], v0DecayVertex, v0DaughterTrackParCovsPerfect[0], pdgDB); + o2::upgrade::convertTLorentzVectorToO2Track(kPiPlus, v0DecayProducts[1], v0DecayVertex, v0DaughterTrackParCovsPerfect[1], pdgDB); + break; + case kLambda0: + o2::upgrade::convertTLorentzVectorToO2Track(kPiMinus, v0DecayProducts[0], v0DecayVertex, v0DaughterTrackParCovsPerfect[0], pdgDB); + o2::upgrade::convertTLorentzVectorToO2Track(kProton, v0DecayProducts[1], v0DecayVertex, v0DaughterTrackParCovsPerfect[1], pdgDB); + break; + case kLambda0Bar: + o2::upgrade::convertTLorentzVectorToO2Track(kProtonBar, v0DecayProducts[0], v0DecayVertex, v0DaughterTrackParCovsPerfect[0], pdgDB); + o2::upgrade::convertTLorentzVectorToO2Track(kPiPlus, v0DecayProducts[1], v0DecayVertex, v0DaughterTrackParCovsPerfect[1], pdgDB); + break; + default: + LOG(fatal) << "Unhandled V0 PDG code " << mcParticle.pdgCode(); + } - // if (TMath::IsNaN(v0DaughterTrackParCovsTracked[i].getZ())) { - // continue; - // } else { - // histos.fill(HIST("hNaNBookkeeping"), i + 1, 1.0f); - // } - if (isReco[i]) { - tracksAlice3.push_back(TrackAlice3{v0DaughterTrackParCovsTracked[i], mcParticle.globalIndex(), time, timeResolutionUs, true, true, i + 2, nSiliconHitsCascadeProngs[i], nTPCHitsCascadeProngs[i]}); - } else { - ghostTracksAlice3.push_back(TrackAlice3{v0DaughterTrackParCovsTracked[i], mcParticle.globalIndex(), time, timeResolutionUs, true, true, i + 2}); + // Store not reconstructed daughters, will update them in case reconstruction is successful + float trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; + tracksV0Daugs.push_back(TrackAlice3{v0DaughterTrackParCovsPerfect[0], mcParticle.globalIndex(), 0.f, timeResolutionUs, true, true, 1, TrackType::kGhostV0Daug}); + trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; + tracksV0Daugs.push_back(TrackAlice3{v0DaughterTrackParCovsPerfect[1], mcParticle.globalIndex(), 0.f, timeResolutionUs, true, true, 1, TrackType::kGhostV0Daug}); + + for (int i = 0; i < kv0Prongs; i++) { + isV0Reco[i] = false; + nV0Hits[i] = 0; + nV0SiliconHits[i] = 0; + nV0TPCHits[i] = 0; + if (enableSecondarySmearing) { + nV0Hits[i] = fastTracker[icfg]->FastTrack(v0DaughterTrackParCovsPerfect[i], v0DaughterTrackParCovsTracked[i], dNdEta); + nV0SiliconHits[i] = fastTracker[icfg]->GetNSiliconPoints(); + nV0TPCHits[i] = fastTracker[icfg]->GetNGasPoints(); + + if (nV0SiliconHits[i] < fastTrackerSettings.minSiliconHits) continue; + if (nV0SiliconHits[i] < fastTrackerSettings.minSiliconHitsIfTPCUsed) continue; + if (nV0TPCHits[i] < fastTrackerSettings.minTPCClusters) continue; + isV0Reco[i] = true; + + if (v0DecaySettings.doV0QA) { // QA + if (nV0Hits[i] < 0) { + fillHist(TH1, Form("V0Building_Configuration_%i/hFastTrackerQA", icfg), o2::math_utils::abs(nV0Hits[i])); + } + for (uint32_t ih = 0; ih < fastTracker[icfg]->GetNHits(); ih++) { + fillHist(TH2, Form("V0Building_Configuration_%i/hFastTrackerHits", icfg), fastTracker[icfg]->GetHitZ(ih), std::hypot(fastTracker[icfg]->GetHitX(ih), fastTracker[icfg]->GetHitY(ih))); } } - if (v0DecaySettings.doV0QA) { - if (isReco[0] && isReco[1]) { - fillHist(TH1, Form("V0Building_Configuration_%i/hV0Building", icfg), 1.0f); - for (size_t indexV0 = 0; indexV0 < v0PDGs.size(); indexV0++) { - if (mcParticle.pdgCode() == v0PDGs[indexV0]) { - fillHist(TH2, Form("V0Building_Configuration_%i/%s/hReco", icfg, NameV0s[indexV0].data()), v0DecayRadius2D, mcParticle.pt()); - } - } + } else { + isV0Reco[i] = true; + v0DaughterTrackParCovsTracked[i] = v0DaughterTrackParCovsPerfect[i]; + } + + trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; + // TODO: flag to separate ghost and reco tracks + TrackType trackType = isV0Reco[i] ? TrackType::kRecoV0Daug : TrackType::kGhostV0Daug; + tracksV0Daugs[i] = TrackAlice3{v0DaughterTrackParCovsTracked[i], mcParticle.globalIndex(), trackTime, timeResolutionUs, true, true, i + 2, trackType}; + } + if (v0DecaySettings.doV0QA) { + if (isV0Reco[0] && isV0Reco[1]) { + fillHist(TH1, Form("V0Building_Configuration_%i/hV0Building", icfg), 1.0f); + for (size_t indexV0 = 0; indexV0 < v0PDGs.size(); indexV0++) { + if (mcParticle.pdgCode() == v0PDGs[indexV0]) { + fillHist(TH2, Form("V0Building_Configuration_%i/%s/hReco", icfg, NameV0s[indexV0].data()), v0DecayRadius2D, mcParticle.pt()); } - if (isReco[0]) { - for (size_t indexV0 = 0; indexV0 < v0PDGs.size(); indexV0++) { - if (mcParticle.pdgCode() == v0PDGs[indexV0]) { - fillHist(TH2, Form("V0Building_Configuration_%i/%s/hRecoNegDaughterFromV0", icfg, NameV0s[indexV0].data()), v0DecayRadius2D, v0DecayProducts[0].Pt()); - } - } + } + } + if (isV0Reco[0]) { + for (size_t indexV0 = 0; indexV0 < v0PDGs.size(); indexV0++) { + if (mcParticle.pdgCode() == v0PDGs[indexV0]) { + fillHist(TH2, Form("V0Building_Configuration_%i/%s/hRecoNegDaughterFromV0", icfg, NameV0s[indexV0].data()), v0DecayRadius2D, v0DecayProducts[0].Pt()); } - if (isReco[1]) { - for (size_t indexV0 = 0; indexV0 < v0PDGs.size(); indexV0++) { - if (mcParticle.pdgCode() == v0PDGs[indexV0]) { - fillHist(TH2, Form("V0Building_Configuration_%i/%s/hRecoPosDaughterFromV0", icfg, NameV0s[indexV0].data()), v0DecayRadius2D, v0DecayProducts[1].Pt()); - } - } + } + } + if (isV0Reco[1]) { + for (size_t indexV0 = 0; indexV0 < v0PDGs.size(); indexV0++) { + if (mcParticle.pdgCode() == v0PDGs[indexV0]) { + fillHist(TH2, Form("V0Building_Configuration_%i/%s/hRecoPosDaughterFromV0", icfg, NameV0s[indexV0].data()), v0DecayRadius2D, v0DecayProducts[1].Pt()); } } + } + } - // +-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+ - // combine particles into actual V0 candidate - // V0 building starts here - if (v0DecaySettings.findV0 && isReco[0] && isReco[1]) { - if (v0DecaySettings.doV0QA) { - fillHist(TH1, Form("V0Building_Configuration_%i/hV0Building", icfg), 2.0f); - } + // +-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+ + // combine particles into actual V0 candidate + // V0 building starts here + if (v0DecaySettings.findV0 && isV0Reco[0] && isV0Reco[1]) { + if (v0DecaySettings.doV0QA) { + fillHist(TH1, Form("V0Building_Configuration_%i/hV0Building", icfg), 2.0f); + } - // assign indices of the particles we've used - // they should be the last ones to be filled, in order: - // n-1: positive Track from V0 - // n-2: negative Track from V0 - thisV0.positiveId = lastTrackIndex + tracksAlice3.size() - 1; - thisV0.negativeId = lastTrackIndex + tracksAlice3.size() - 2; - thisV0.mcParticleId = mcParticle.globalIndex(); - // use DCA fitters - int nCand = 0; - bool dcaFitterOK_V0 = true; - try { - nCand = fitter.process(v0DaughterTrackParCovsTracked[0], v0DaughterTrackParCovsTracked[1]); - } catch (...) { - // LOG(error) << "Exception caught in DCA fitter process call!"; - dcaFitterOK_V0 = false; - } - if (nCand == 0) { - dcaFitterOK_V0 = false; - } + // assign indices of the particles we've used + // they should be the last ones to be filled, in order: + // n-1: positive Track from V0 + // n-2: negative Track from V0 + thisV0.mcParticleId = mcParticle.globalIndex(); + thisV0.negativeId = trackTableOffset + 1; + thisV0.positiveId = trackTableOffset + 2; + + // use DCA fitters + int nCand = 0; + bool dcaFitterOK_V0 = true; + try { + nCand = fitter.process(v0DaughterTrackParCovsTracked[0], v0DaughterTrackParCovsTracked[1]); + } catch (...) { + // LOG(error) << "Exception caught in DCA fitter process call!"; + dcaFitterOK_V0 = false; + } + if (nCand == 0) { + dcaFitterOK_V0 = false; + } - fitter.propagateTracksToVertex(); - if (!fitter.isPropagateTracksToVertexDone()) { - dcaFitterOK_V0 = false; - } + fitter.propagateTracksToVertex(); + if (!fitter.isPropagateTracksToVertexDone()) { + dcaFitterOK_V0 = false; + } - const u_int8_t fitterStatusCode = fitter.getFitStatus(); - histos.fill(HIST("hFitterStatusCode"), fitterStatusCode); - // V0 found successfully - if (dcaFitterOK_V0) { - if (v0DecaySettings.doV0QA) { - fillHist(TH1, Form("V0Building_Configuration_%i/hV0Building", icfg), 3.0f); - } + const u_int8_t fitterStatusCode = fitter.getFitStatus(); + histos.fill(HIST("hFitterStatusCode"), fitterStatusCode); + // V0 found successfully + if (dcaFitterOK_V0) { + if (v0DecaySettings.doV0QA) { + fillHist(TH1, Form("V0Building_Configuration_%i/hV0Building", icfg), 3.0f); + } - std::array pos; - std::array posP; - std::array negP; + std::array pos; + std::array posP; + std::array negP; - o2::track::TrackParCov pTrackAtPCA = fitter.getTrack(1); // (positive) - o2::track::TrackParCov nTrackAtPCA = fitter.getTrack(0); // (negative) - pTrackAtPCA.getPxPyPzGlo(posP); - nTrackAtPCA.getPxPyPzGlo(negP); + o2::track::TrackParCov pTrackAtPCA = fitter.getTrack(1); // (positive) + o2::track::TrackParCov nTrackAtPCA = fitter.getTrack(0); // (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]; - } + // 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 - thisV0.dcaV0dau = TMath::Sqrt(fitter.getChi2AtPCACandidate()); - thisV0.v0radius = std::hypot(pos[0], pos[1]); - thisV0.pt = std::hypot(std::cos(v0DaughterTrackParCovsTracked[0].getPhi()) * v0DaughterTrackParCovsTracked[0].getPt() + std::cos(v0DaughterTrackParCovsTracked[1].getPhi()) * v0DaughterTrackParCovsTracked[1].getPt(), - std::sin(v0DaughterTrackParCovsTracked[0].getPhi()) * v0DaughterTrackParCovsTracked[0].getPt() + std::sin(v0DaughterTrackParCovsTracked[1].getPhi()) * v0DaughterTrackParCovsTracked[1].getPt()); - if (std::abs(mcParticle.pdgCode()) == kK0Short) { - thisV0.mK0 = RecoDecay::m(std::array{std::array{posP[0], posP[1], posP[2]}, + // calculate basic V0 properties here + // DCA to PV taken care of in daughter tracks already, not necessary + thisV0.dcaV0dau = TMath::Sqrt(fitter.getChi2AtPCACandidate()); + thisV0.v0radius = std::hypot(pos[0], pos[1]); + thisV0.pt = std::hypot(std::cos(v0DaughterTrackParCovsTracked[0].getPhi()) * v0DaughterTrackParCovsTracked[0].getPt() + std::cos(v0DaughterTrackParCovsTracked[1].getPhi()) * v0DaughterTrackParCovsTracked[1].getPt(), + std::sin(v0DaughterTrackParCovsTracked[0].getPhi()) * v0DaughterTrackParCovsTracked[0].getPt() + std::sin(v0DaughterTrackParCovsTracked[1].getPhi()) * v0DaughterTrackParCovsTracked[1].getPt()); + + thisV0.mK0 = -1, thisV0.mLambda = -1, thisV0.mAntiLambda = -1; // initialize to unphysical values + if (std::abs(mcParticle.pdgCode()) == kK0Short) { + thisV0.mK0 = RecoDecay::m(std::array{std::array{posP[0], posP[1], posP[2]}, + std::array{negP[0], negP[1], negP[2]}}, + std::array{o2::constants::physics::MassPionCharged, + o2::constants::physics::MassPionCharged}); + } + + if (mcParticle.pdgCode() == kLambda0) { + thisV0.mLambda = RecoDecay::m(std::array{std::array{posP[0], posP[1], posP[2]}, std::array{negP[0], negP[1], negP[2]}}, - std::array{o2::constants::physics::MassPionCharged, + std::array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged}); - } else { - thisV0.mK0 = -1; - } + } - if (mcParticle.pdgCode() == kLambda0) { - thisV0.mLambda = RecoDecay::m(std::array{std::array{posP[0], posP[1], posP[2]}, + if (mcParticle.pdgCode() == kLambda0Bar) { + thisV0.mAntiLambda = RecoDecay::m(std::array{std::array{posP[0], posP[1], posP[2]}, std::array{negP[0], negP[1], negP[2]}}, - std::array{o2::constants::physics::MassProton, - o2::constants::physics::MassPionCharged}); - } else { - thisV0.mLambda = -1; - } - - if (mcParticle.pdgCode() == kLambda0Bar) { - thisV0.mAntiLambda = RecoDecay::m(std::array{std::array{posP[0], posP[1], posP[2]}, - std::array{negP[0], negP[1], negP[2]}}, - std::array{o2::constants::physics::MassPionCharged, - o2::constants::physics::MassProton}); - } else { - thisV0.mAntiLambda = -1; - } - - if (v0DecaySettings.doV0QA) { - fillHist(TH1, Form("V0Building_Configuration_%i/hV0Building", icfg), 4.0f); - if (std::abs(mcParticle.pdgCode()) == kK0Short) { - fillHist(TH2, Form("V0Building_Configuration_%i/K0/hMass", icfg), thisV0.mK0, thisV0.pt); - } - if (mcParticle.pdgCode() == kLambda0) { - fillHist(TH2, Form("V0Building_Configuration_%i/Lambda/hMass", icfg), thisV0.mLambda, thisV0.pt); - } - if (mcParticle.pdgCode() == kLambda0Bar) { - fillHist(TH2, Form("V0Building_Configuration_%i/AntiLambda/hMass", icfg), thisV0.mAntiLambda, thisV0.pt); - } - } + std::array{o2::constants::physics::MassPionCharged, + o2::constants::physics::MassProton}); + } - // add this V0 to vector (will fill cursor later with collision ID) - v0sAlice3.push_back(thisV0); + if (v0DecaySettings.doV0QA) { + fillHist(TH1, Form("V0Building_Configuration_%i/hV0Building", icfg), 4.0f); + if (std::abs(mcParticle.pdgCode()) == kK0Short) { + fillHist(TH2, Form("V0Building_Configuration_%i/K0/hMass", icfg), thisV0.mK0, thisV0.pt); + } + if (mcParticle.pdgCode() == kLambda0) { + fillHist(TH2, Form("V0Building_Configuration_%i/Lambda/hMass", icfg), thisV0.mLambda, thisV0.pt); + } + if (mcParticle.pdgCode() == kLambda0Bar) { + fillHist(TH2, Form("V0Building_Configuration_%i/AntiLambda/hMass", icfg), thisV0.mAntiLambda, thisV0.pt); } } } + } + + // populate V0s + tableUpgradeV0s(tableCollisions.lastIndex(), // now we know the collision index -> populate table + thisV0.mcParticleId, + thisV0.positiveId, + thisV0.negativeId, + thisV0.dcaV0dau, + thisV0.v0radius, + thisV0.mLambda, + thisV0.mAntiLambda, + thisV0.mK0, + thisV0.pt); + } + + /// Function to compute the primary vertex position using the provided tracks and MC collision information + /// \param mcCollision the MC collision information, used to get the true vertex position for comparison + /// \param prmTrks the vector of tracks to be used for vertex reconstruction + /// \param primaryVertex the output variable where the computed primary vertex will be stored + /// \param icfg index of the current configuration, used for histogram filling + template + void computeVertex(McCollisionType& mcCollision, const std::vector& prmTrks, o2::vertexing::PVertex& primaryVertex, const int icfg) { + + if (!enablePrimaryVertexing) { + primaryVertex.setXYZ(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ()); + return; + } + vertexReconstructionEfficiencyCounters.first += 1; + const std::string histPath = "Configuration_" + std::to_string(icfg) + "/"; + fillHist(TH1, histPath + "hVtxMultGen", prmTrks.size()); + std::vector lblTracks; + std::vector vertices; + std::vector vertexTrackIDs; + std::vector v2tRefs; + std::vector lblVtx; + lblVtx.emplace_back(mcCollision.globalIndex(), 1); + std::vector idxVec; // store IDs + + idxVec.reserve(prmTrks.size()); + for (unsigned i = 0; i < prmTrks.size(); i++) { + lblTracks.emplace_back(prmTrks[i].mcLabel, mcCollision.globalIndex(), 1, false); + idxVec.emplace_back(i, o2::dataformats::GlobalTrackID::ITS); // let's say ITS + } + + getHist(TH1, histPath + "hVtxTrials")->Fill(0); // Tried vertexing + // Calculate vertices + const int n_vertices = vertexer.process(prmTrks, // track array + idxVec, + gsl::span{bcData}, + vertices, + vertexTrackIDs, + v2tRefs, + gsl::span{lblTracks}, + lblVtx); + + if (n_vertices < 1) { if (doExtraQA) { - histos.fill(HIST("hSimTrackX"), trackParCov.getX()); + histos.fill(HIST("h1dVerticesNotReco"), prmTrks.size()); } - if (isV0) { - continue; // V0 handling done, should not be considered anymore - } - - bool reconstructed = true; - int nTrkHits = 0; - if (enablePrimarySmearing && !fastPrimaryTrackerSettings.fastTrackPrimaries) { - reconstructed = mSmearer[icfg]->smearTrack(trackParCov, mcParticle.pdgCode(), dNdEta); - nTrkHits = fastTrackerSettings.minSiliconHits; - } else if (fastPrimaryTrackerSettings.fastTrackPrimaries) { - o2::track::TrackParCov o2Track; - o2::upgrade::convertMCParticleToO2Track(mcParticle, o2Track, pdgDB); - o2Track.setPID(pdgCodeToPID(mcParticle.pdgCode())); - nTrkHits = fastTracker[icfg]->FastTrack(o2Track, trackParCov, dNdEta); - if (nTrkHits < fastPrimaryTrackerSettings.minSiliconHits) { - reconstructed = false; - } - } - - if (!reconstructed && !processUnreconstructedTracks) { - continue; - } - if (TMath::IsNaN(trackParCov.getZ())) { - // capture rare smearing mistakes / corrupted tracks - histos.fill(HIST("hNaNBookkeeping"), 0.0f, 0.0f); - continue; - } else { - histos.fill(HIST("hNaNBookkeeping"), 0.0f, 1.0f); // ok! - } - - // Base QA (note: reco pT here) - if (enablePrimarySmearing) { - getHist(TH1, histPath + "hPtReconstructed")->Fill(trackParCov.getPt()); - if (std::abs(mcParticle.pdgCode()) == kElectron) - getHist(TH1, histPath + "hPtReconstructedEl")->Fill(trackParCov.getPt()); - if (std::abs(mcParticle.pdgCode()) == kPiPlus) - getHist(TH1, histPath + "hPtReconstructedPi")->Fill(trackParCov.getPt()); - if (std::abs(mcParticle.pdgCode()) == kKPlus) - getHist(TH1, histPath + "hPtReconstructedKa")->Fill(trackParCov.getPt()); - if (std::abs(mcParticle.pdgCode()) == kProton) - getHist(TH1, histPath + "hPtReconstructedPr")->Fill(trackParCov.getPt()); - } - if (doExtraQA) { - getHist(TH2, histPath + "h2dPtRes")->Fill(trackParCov.getPt(), (trackParCov.getPt() - mcParticle.pt()) / trackParCov.getPt()); - histos.fill(HIST("hRecoTrackX"), trackParCov.getX()); - } - - // populate vector with track if we reco-ed it - if (reconstructed) { - tracksAlice3.push_back(TrackAlice3{trackParCov, mcParticle.globalIndex(), time, timeResolutionUs, isDecayDaughter, false, 0, nTrkHits}); - } else { - ghostTracksAlice3.push_back(TrackAlice3{trackParCov, mcParticle.globalIndex(), time, timeResolutionUs, isDecayDaughter}); + return; // primary vertex not reconstructed + } + vertexReconstructionEfficiencyCounters.second += 1; + getHist(TH1, histPath + "hVtxTrials")->Fill(1); // Succeeded vertexing + + // Find largest vertex + int largestVertex = 0; + for (Int_t iv = 1; iv < n_vertices; iv++) { + if (vertices[iv].getNContributors() > vertices[largestVertex].getNContributors()) { + largestVertex = iv; } } - - // *+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+* - // Calculate primary vertex with tracks from this collision - // data preparation - o2::vertexing::PVertex primaryVertex; - if (enablePrimaryVertexing) { - LOG(debug) << "Starting primary vertexing with " << tracksAlice3.size() << " tracks."; - fillHist(TH1, histPath + "hVtxMultGen", tracksAlice3.size()); - std::vector lblTracks; - std::vector vertices; - std::vector vertexTrackIDs; - std::vector v2tRefs; - std::vector lblVtx; - lblVtx.emplace_back(mcCollision.globalIndex(), 1); - std::vector idxVec; // store IDs - - idxVec.reserve(tracksAlice3.size()); - for (unsigned i = 0; i < tracksAlice3.size(); i++) { - lblTracks.emplace_back(tracksAlice3[i].mcLabel, mcCollision.globalIndex(), 1, false); - idxVec.emplace_back(i, o2::dataformats::GlobalTrackID::ITS); // let's say ITS - } - - getHist(TH1, histPath + "hVtxTrials")->Fill(0); // Tried vertexing - - // Calculate vertices - const int n_vertices = vertexer.process(tracksAlice3, // track array - idxVec, - gsl::span{bcData}, - vertices, - vertexTrackIDs, - v2tRefs, - gsl::span{lblTracks}, - lblVtx); - - LOG(debug) << "Vertex reconstruction efficiency " << vertexReconstructionEfficiencyCounters.second << "/" << vertexReconstructionEfficiencyCounters.first << "=" << vertexReconstructionEfficiencyCounters.second / vertexReconstructionEfficiencyCounters.first; - if (n_vertices < 1) { - LOG(debug) << "Vertexing completed, no vtx found."; - if (doExtraQA) { - histos.fill(HIST("h1dVerticesNotReco"), tracksAlice3.size()); - } - return; // primary vertex not reconstructed - } - vertexReconstructionEfficiencyCounters.second += 1; - getHist(TH1, histPath + "hVtxTrials")->Fill(1); // Succeeded vertexing - LOG(debug) << "Vertexing completed with " << n_vertices << " vertices found."; - - // Find largest vertex - int largestVertex = 0; - for (Int_t iv = 1; iv < n_vertices; iv++) { - if (vertices[iv].getNContributors() > vertices[largestVertex].getNContributors()) { - largestVertex = iv; - } - } - primaryVertex = vertices[largestVertex]; - if (doExtraQA) { - histos.fill(HIST("h2dVerticesVsContributors"), primaryVertex.getNContributors(), n_vertices); - } - fillHist(TH1, histPath + "hVtxMultReco", primaryVertex.getNContributors()); - fillHist(TH1, histPath + "hDeltaMultPVRecoGen", static_cast(primaryVertex.getNContributors()) - static_cast(tracksAlice3.size())); - fillHist(TH1, histPath + "hDeltaXPVRecoGen", primaryVertex.getX() - mcCollision.posX()); - fillHist(TH1, histPath + "hDeltaYPVRecoGen", primaryVertex.getY() - mcCollision.posY()); - fillHist(TH1, histPath + "hDeltaZPVRecoGen", primaryVertex.getZ() - mcCollision.posZ()); - } else { - primaryVertex.setXYZ(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ()); - } - // *+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+* - - // debug / informational - getHist(TH1, histPath + "hSimMultiplicity")->Fill(multiplicityCounter); - getHist(TH1, histPath + "hRecoMultiplicity")->Fill(tracksAlice3.size()); - getHist(TH1, histPath + "hPVz")->Fill(primaryVertex.getZ()); - + primaryVertex = vertices[largestVertex]; if (doExtraQA) { - histos.fill(HIST("hRecoVsSimMultiplicity"), multiplicityCounter, tracksAlice3.size()); + histos.fill(HIST("h2dVerticesVsContributors"), primaryVertex.getNContributors(), n_vertices); } + fillHist(TH1, histPath + "hVtxMultReco", primaryVertex.getNContributors()); + fillHist(TH1, histPath + "hDeltaMultPVRecoGen", static_cast(primaryVertex.getNContributors()) - static_cast(prmTrks.size())); + fillHist(TH2, histPath + "hDeltaXPVRecoGen", primaryVertex.getX() - mcCollision.posX(), primaryVertex.getNContributors()); + fillHist(TH2, histPath + "hDeltaYPVRecoGen", primaryVertex.getY() - mcCollision.posY(), primaryVertex.getNContributors()); + fillHist(TH2, histPath + "hDeltaZPVRecoGen", primaryVertex.getZ() - mcCollision.posZ(), primaryVertex.getNContributors()); + } - // *+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+* - // populate collisions - tableCollisions(-1, // BC is irrelevant in synthetic MC tests for now, could be adjusted in future - primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ(), - primaryVertex.getSigmaX2(), primaryVertex.getSigmaXY(), primaryVertex.getSigmaY2(), - primaryVertex.getSigmaXZ(), primaryVertex.getSigmaYZ(), primaryVertex.getSigmaZ2(), - 0, primaryVertex.getChi2(), primaryVertex.getNContributors(), - eventCollisionTimeNS, 0.f); // For the moment the event collision time is taken as the "GEANT" time, the computation of the event time is done a posteriori from the tracks in the OTF TOF PID task - tableMcCollisionLabels(mcCollision.globalIndex(), 0); - tableCollisionsAlice3(dNdEta); - tableOTFLUTConfigId(icfg); // Populate OTF LUT configuration ID - // *+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+* + /// Function to fill track information into the relevant tables and histograms + /// \param tracks the vector of tracks to be processed + /// \param primaryVertex the primary vertex position, used for DCA calculation + /// \param icfg index of the current configuration, used for histogram filling + void fillTracksInfo(std::vector const& tracks, o2::vertexing::PVertex const& primaryVertex, const int icfg) { + const std::string histPath = "Configuration_" + std::to_string(icfg) + "/"; - // *+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+* - // populate tracks - LOG(debug) << "Populating " << tracksAlice3.size() << " tracks."; - for (const auto& trackParCov : tracksAlice3) { + for (const auto& trackParCov : tracks) { // Fixme: collision index could be changeable aod::track::TrackTypeEnum trackType = aod::track::Track; @@ -1721,24 +1692,6 @@ struct OnTheFlyTracker { getHist(TH2, histPath + "h2dDCAz")->Fill(trackParametrization.getPt(), dcaZ * 1e+4); histos.fill(HIST("hTrackXatDCA"), trackParametrization.getX()); } - if (cascadeDecaySettings.doXiQA) { - if (trackParCov.isUsedInCascading == 1) { - getHist(TH2, histPath + "h2dDCAxyCascade")->Fill(trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please - getHist(TH2, histPath + "h2dDCAzCascade")->Fill(trackParametrization.getPt(), dcaZ * 1e+4); // in microns, please - } - if (trackParCov.isUsedInCascading == 2) { - getHist(TH2, histPath + "h2dDCAxyCascadeBachelor")->Fill(trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please - getHist(TH2, histPath + "h2dDCAzCascadeBachelor")->Fill(trackParametrization.getPt(), dcaZ * 1e+4); // in microns, please - } - if (trackParCov.isUsedInCascading == 3) { - getHist(TH2, histPath + "h2dDCAxyCascadeNegative")->Fill(trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please - getHist(TH2, histPath + "h2dDCAzCascadeNegative")->Fill(trackParametrization.getPt(), dcaZ * 1e+4); // in microns, please - } - if (trackParCov.isUsedInCascading == 4) { - getHist(TH2, histPath + "h2dDCAxyCascadePositive")->Fill(trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please - getHist(TH2, histPath + "h2dDCAzCascadePositive")->Fill(trackParametrization.getPt(), dcaZ * 1e+4); // in microns, please - } - } tableTracksDCA(dcaXY, dcaZ); if (populateTracksDCACov) { tableTracksDCACov(dcaInfo.getSigmaY2(), dcaInfo.getSigmaZ2()); @@ -1755,7 +1708,7 @@ struct OnTheFlyTracker { trackParCov.getSigmaTgl2(), trackParCov.getSigma1PtY(), trackParCov.getSigma1PtZ(), trackParCov.getSigma1PtSnp(), trackParCov.getSigma1PtTgl(), trackParCov.getSigma1Pt2()); tableMcTrackLabels(trackParCov.mcLabel, 0); - tableTracksExtraA3(trackParCov.nSiliconHits, trackParCov.nTPCHits); + tableTracksExtraA3(trackParCov.nSiliconHits, trackParCov.nTPCHits, trackParCov.trackType); // populate extra tables if required to do so if (populateTracksExtra) { @@ -1770,97 +1723,178 @@ struct OnTheFlyTracker { } tableTracksAlice3(true); } + } - // populate ghost tracks - LOG(debug) << "Populating " << ghostTracksAlice3.size() << " ghost tracks."; - for (const auto& trackParCov : ghostTracksAlice3) { - // Fixme: collision index could be changeable - aod::track::TrackTypeEnum trackType = aod::track::Track; + void processWithLUTs(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles, const int icfg) + { + LOG(info) << "Processing collision " << mcCollision.globalIndex() << " with " << mcParticles.size() << " particles"; + const std::string histPath = "Configuration_" + std::to_string(icfg) + "/"; - if (populateTracksDCA) { - float dcaXY = 1e+10, dcaZ = 1e+10; - o2::track::TrackParCov trackParametrization(trackParCov); - if (trackParametrization.propagateToDCA(primaryVertex, mMagneticField, &dcaInfo)) { - dcaXY = dcaInfo.getY(); - dcaZ = dcaInfo.getZ(); - } - if (doExtraQA && (!extraQAwithoutDecayDaughters || (extraQAwithoutDecayDaughters && !trackParCov.isDecayDau))) { - histos.fill(HIST("h2dDCAxy"), trackParametrization.getPt(), dcaXY * 1e+4); // in microns, please - histos.fill(HIST("h2dDCAz"), trackParametrization.getPt(), dcaZ * 1e+4); // in microns, please - histos.fill(HIST("hTrackXatDCA"), trackParametrization.getX()); - } - tableTracksDCA(dcaXY, dcaZ); - if (populateTracksDCACov) { - tableTracksDCACov(dcaInfo.getSigmaY2(), dcaInfo.getSigmaZ2()); - } + std::vector genCascades; + std::vector genV0s; + bcData.clear(); + + o2::dataformats::DCA dcaInfo; + // o2::dataformats::VertexBase vtx; + + // *+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+* + // Study collision and perform vertexing + float dNdEta{0.f}; // Charged particle multiplicity to use in the efficiency evaluation + computeDNDEta(dNdEta, mcParticles, histPath); + auto ir = irSampler.generateCollisionTime(); + const float eventCollisionTimeNS = ir.timeInBCNS; + + uint32_t multiplicityCounter = 0; + // Now that the multiplicity is known, we can process the particles to smear them + for (const auto& mcParticle : mcParticles) { + + if ((std::fabs(mcParticle.eta()) > maxEta) || (mcParticle.pt() < minPt)) { + continue; + } + const bool isCascade = mcParticle.pdgCode() == kXiMinus; + if (isCascade && cascadeDecaySettings.decayXi) { + genCascades.push_back(mcParticle.globalIndex()); + } + const bool isV0 = std::find(v0PDGs.begin(), v0PDGs.end(), mcParticle.pdgCode()) != v0PDGs.end(); + if (isV0 && v0DecaySettings.decayV0) { + genV0s.push_back(mcParticle.globalIndex()); + } + if (!mcParticle.isPhysicalPrimary()) { + continue; } - tableStoredTracks(tableCollisions.lastIndex(), trackType, trackParCov.getX(), trackParCov.getAlpha(), trackParCov.getY(), trackParCov.getZ(), trackParCov.getSnp(), trackParCov.getTgl(), trackParCov.getQ2Pt()); - tableTracksExtension(trackParCov.getPt(), trackParCov.getP(), trackParCov.getEta(), trackParCov.getPhi()); + const bool longLivedToBeHandled = std::find(longLivedHandledPDGs.begin(), longLivedHandledPDGs.end(), std::abs(mcParticle.pdgCode())) != longLivedHandledPDGs.end(); + const bool nucleiToBeHandled = std::find(nucleiPDGs.begin(), nucleiPDGs.end(), std::abs(mcParticle.pdgCode())) != nucleiPDGs.end(); + const bool pdgsToBeHandled = longLivedToBeHandled || + (enableNucleiSmearing && nucleiToBeHandled) || + (cascadeDecaySettings.decayXi && isCascade) || + (v0DecaySettings.decayV0 && isV0); + if (!pdgsToBeHandled) { + // LOG(info) << "Particle with PDG " << mcParticle.pdgCode() << " is not in the list of particles to be handled, skipping."; + continue; + } - // TODO do we keep the rho as 0? Also the sigma's are duplicated information - tableStoredTracksCov(std::sqrt(trackParCov.getSigmaY2()), std::sqrt(trackParCov.getSigmaZ2()), std::sqrt(trackParCov.getSigmaSnp2()), - std::sqrt(trackParCov.getSigmaTgl2()), std::sqrt(trackParCov.getSigma1Pt2()), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); - tableTracksCovExtension(trackParCov.getSigmaY2(), trackParCov.getSigmaZY(), trackParCov.getSigmaZ2(), trackParCov.getSigmaSnpY(), - trackParCov.getSigmaSnpZ(), trackParCov.getSigmaSnp2(), trackParCov.getSigmaTglY(), trackParCov.getSigmaTglZ(), trackParCov.getSigmaTglSnp(), - trackParCov.getSigmaTgl2(), trackParCov.getSigma1PtY(), trackParCov.getSigma1PtZ(), trackParCov.getSigma1PtSnp(), trackParCov.getSigma1PtTgl(), - trackParCov.getSigma1Pt2()); - tableMcTrackLabels(trackParCov.mcLabel, 0); - tableTracksExtraA3(trackParCov.nSiliconHits, trackParCov.nTPCHits); + multiplicityCounter++; + o2::track::TrackParCov trackParCov; + const bool isDecayDaughter = (mcParticle.getProcess() == TMCProcess::kPDecay); + if (doExtraQA) { + histos.fill(HIST("hSimTrackX"), trackParCov.getX()); + } - // populate extra tables if required to do so - if (populateTracksExtra) { - tableStoredTracksExtra(0.0f, static_cast(0), static_cast(0), static_cast(0), static_cast(0), - static_cast(0), static_cast(0), static_cast(0), static_cast(0), - 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, - 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f); + bool reconstructed = true; + int nTrkHits = 0; + if (enablePrimarySmearing) { + if (fastPrimaryTrackerSettings.fastTrackPrimaries) { + o2::track::TrackParCov perfectTrackParCov; + o2::upgrade::convertMCParticleToO2Track(mcParticle, perfectTrackParCov, pdgDB); + perfectTrackParCov.setPID(pdgCodeToPID(mcParticle.pdgCode())); + nTrkHits = fastTracker[icfg]->FastTrack(perfectTrackParCov, trackParCov, dNdEta); + if (nTrkHits < fastPrimaryTrackerSettings.minSiliconHits) { + reconstructed = false; + } + } else { + o2::upgrade::convertMCParticleToO2Track(mcParticle, trackParCov, pdgDB); + reconstructed = mSmearer[icfg]->smearTrack(trackParCov, mcParticle.pdgCode(), dNdEta); + nTrkHits = fastTrackerSettings.minSiliconHits; + } + getHist(TH1, histPath + "hPtGenerated")->Fill(mcParticle.pt()); + getHist(TH1, histPath + "hPhiGenerated")->Fill(mcParticle.phi()); + if (std::abs(mcParticle.pdgCode()) == kElectron) getHist(TH1, histPath + "hPtGeneratedEl")->Fill(mcParticle.pt()); + if (std::abs(mcParticle.pdgCode()) == kPiPlus) getHist(TH1, histPath + "hPtGeneratedPi")->Fill(mcParticle.pt()); + if (std::abs(mcParticle.pdgCode()) == kKPlus) getHist(TH1, histPath + "hPtGeneratedKa")->Fill(mcParticle.pt()); + if (std::abs(mcParticle.pdgCode()) == kProton) getHist(TH1, histPath + "hPtGeneratedPr")->Fill(mcParticle.pt()); + + if (!reconstructed && !processUnreconstructedTracks) { + continue; + } + + // LOG(info) << "Particle with PDG " << mcParticle.pdgCode() << (reconstructed ? " was reconstructed." : " was not reconstructed."); + getHist(TH1, histPath + "hPtReconstructed")->Fill(trackParCov.getPt()); + if (std::abs(mcParticle.pdgCode()) == kElectron) getHist(TH1, histPath + "hPtReconstructedEl")->Fill(trackParCov.getPt()); + if (std::abs(mcParticle.pdgCode()) == kPiPlus) getHist(TH1, histPath + "hPtReconstructedPi")->Fill(trackParCov.getPt()); + if (std::abs(mcParticle.pdgCode()) == kKPlus) getHist(TH1, histPath + "hPtReconstructedKa")->Fill(trackParCov.getPt()); + if (std::abs(mcParticle.pdgCode()) == kProton) getHist(TH1, histPath + "hPtReconstructedPr")->Fill(trackParCov.getPt()); } - if (populateTrackSelection) { - tableTrackSelection(static_cast(0), false, false, false, false, false, false); - tableTrackSelectionExtension(false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false); + if (doExtraQA) { + getHist(TH2, histPath + "h2dPtRes")->Fill(trackParCov.getPt(), (trackParCov.getPt() - mcParticle.pt()) / trackParCov.getPt()); + getHist(TH2, histPath + "h2dPtResAbs")->Fill(trackParCov.getPt(), trackParCov.getPt() - mcParticle.pt()); + histos.fill(HIST("hRecoTrackX"), trackParCov.getX()); + } + + if (TMath::IsNaN(trackParCov.getZ())) { + histos.fill(HIST("hNaNBookkeeping"), 0.0f, 0.0f); + continue; // capture smearing mistakes + } + histos.fill(HIST("hNaNBookkeeping"), 0.0f, 1.0f); // ok! + + // LOG(info) << "Particle with PDG " << mcParticle.pdgCode() << " has " << nTrkHits << " hits."; + // Time associated to the mcParticle: collision time + smearing + float trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; + TrackType trackType = reconstructed ? TrackType::kRecoPrimary : TrackType::kGhostPrimary; + if (reconstructed) { + recoPrimaries.push_back(TrackAlice3{trackParCov, mcParticle.globalIndex(), trackTime, timeResolutionUs, isDecayDaughter, false, 0, nTrkHits, trackType}); + } else { + ghostPrimaries.push_back(TrackAlice3{trackParCov, mcParticle.globalIndex(), trackTime, timeResolutionUs, isDecayDaughter, false, 0, nTrkHits, trackType}); } - tableTracksAlice3(false); } - // populate Cascades - LOG(debug) << "Populating " << cascadesAlice3.size() << " cascades."; - for (const auto& cascade : cascadesAlice3) { - tableUpgradeCascades(tableCollisions.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); + // Compute primary vertex with primary tracks + o2::vertexing::PVertex primaryVertex; + if (recoPrimaries.size() <= 2) { + LOG(info) << "Not enough primary tracks (" << recoPrimaries.size() << ") to compute vertex, skipping vertexing and collision population."; + return; } + computeVertex(mcCollision, recoPrimaries, primaryVertex, icfg); + // LOG(info) << "Primary vertex position: (" << primaryVertex.getX() << ", " << primaryVertex.getY() << ", " << primaryVertex.getZ() << ")"; + getHist(TH1, histPath + "hPVz")->Fill(primaryVertex.getZ()); + // populate collisions + tableCollisions(-1, // BC is irrelevant in synthetic MC tests for now, could be adjusted in future + primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ(), + primaryVertex.getSigmaX2(), primaryVertex.getSigmaXY(), primaryVertex.getSigmaY2(), + primaryVertex.getSigmaXZ(), primaryVertex.getSigmaYZ(), primaryVertex.getSigmaZ2(), + 0, primaryVertex.getChi2(), primaryVertex.getNContributors(), + eventCollisionTimeNS, 0.f); // For the moment the event collision time is taken as the "GEANT" time, the computation of the event time is done a posteriori from the tracks in the OTF TOF PID task + tableMcCollisionLabels(mcCollision.globalIndex(), 0); + tableCollisionsAlice3(dNdEta); + tableOTFLUTConfigId(icfg); // Populate OTF LUT configuration ID - // populate V0s - LOG(debug) << "Populating " << v0sAlice3.size() << " V0s."; - for (const auto& v0 : v0sAlice3) { - tableUpgradeV0s(tableCollisions.lastIndex(), // now we know the collision index -> populate table - v0.mcParticleId, - v0.positiveId, - v0.negativeId, - v0.dcaV0dau, - v0.v0radius, - v0.mLambda, - v0.mAntiLambda, - v0.mK0, - v0.pt); + // Study V0s and cascades + int trackTableOffset = tableStoredTracks.lastIndex(); + std::vector tracksCascadeProngs; + for (size_t iCasc = 0; iCasc < genCascades.size(); iCasc++) { + auto genCasc = mcParticles.rawIteratorAt(genCascades[iCasc] - mcParticles.offset()); + studyCascade(trackTableOffset, genCasc, tracksCascadeProngs, primaryVertex, icfg, dNdEta, eventCollisionTimeNS); + LOG(info) << "[CASCADE] Populating number of cascade tracks: " << tracksCascadeProngs.size(); + fillTracksInfo(tracksCascadeProngs, primaryVertex, icfg); + trackTableOffset += tracksCascadeProngs.size(); // each cascade takes 4 tracks in the table (cascade + 3 daughters) + tracksCascadeProngs.clear(); + } + std::vector tracksV0Daugs; + for (size_t iV0 = 0; iV0 < genV0s.size(); iV0++) { + auto genV0 = mcParticles.rawIteratorAt(genV0s[iV0] - mcParticles.offset()); + studyV0(trackTableOffset, genV0, tracksV0Daugs, icfg, dNdEta, eventCollisionTimeNS); + LOG(info) << "Populating " << tracksV0Daugs.size() << " tracks associated to V0 candidate"; + fillTracksInfo(tracksV0Daugs, primaryVertex, icfg); + trackTableOffset += tracksV0Daugs.size(); // each V0 takes 2 tracks in the table (2 daughters) + tracksV0Daugs.clear(); } + // *+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+* + // populate tracks + fillTracksInfo(recoPrimaries, primaryVertex, icfg); + fillTracksInfo(ghostPrimaries, primaryVertex, icfg); + // do bookkeeping of fastTracker tracking if (enableSecondarySmearing) { histos.fill(HIST("hCovMatOK"), 0.0f, fastTracker[icfg]->GetCovMatNotOK()); histos.fill(HIST("hCovMatOK"), 1.0f, fastTracker[icfg]->GetCovMatOK()); } + if (doExtraQA) { + histos.fill(HIST("hRecoVsSimMultiplicity"), multiplicityCounter, recoPrimaries.size()); + getHist(TH1, histPath + "hSimMultiplicity")->Fill(multiplicityCounter); + getHist(TH1, histPath + "hRecoMultiplicity")->Fill(recoPrimaries.size()); + } + LOG(debug) << " <- Finished processing OTF tracking with LUT configuration ID " << icfg; } // end process @@ -1871,7 +1905,6 @@ struct OnTheFlyTracker { processWithLUTs(mcCollision, mcParticles, static_cast(icfg)); } } - void processConfigurationDev(aod::McCollision const& mcCollision, aod::McPartWithDaus const& mcParticles, const int icfg) { // const int lastTrackIndex = tableStoredTracksCov.lastIndex() + 1; // bookkeep the last added track @@ -1890,39 +1923,10 @@ struct OnTheFlyTracker { const float eventCollisionTimeNS = ir.timeInBCNS; // First we compute the number of charged particles in the event - dNdEta = 0.f; - for (const auto& mcParticle : mcParticles) { - if (std::abs(mcParticle.eta()) > multEtaRange) { - continue; - } - - if (!mcParticle.isPhysicalPrimary()) { - continue; - } - - const auto pdg = std::abs(mcParticle.pdgCode()); - const bool longLivedToBeHandled = std::find(longLivedHandledPDGs.begin(), longLivedHandledPDGs.end(), pdg) != longLivedHandledPDGs.end(); - const bool nucleiToBeHandled = std::find(nucleiPDGs.begin(), nucleiPDGs.end(), pdg) != nucleiPDGs.end(); - const bool pdgsToBeHandled = longLivedToBeHandled || (enableNucleiSmearing && nucleiToBeHandled); - if (!pdgsToBeHandled) { - continue; - } - - const auto& pdgInfo = pdgDB->GetParticle(mcParticle.pdgCode()); - if (!pdgInfo) { - LOG(warning) << "PDG code " << mcParticle.pdgCode() << " not found in the database"; - continue; - } - if (pdgInfo->Charge() == 0) { - continue; - } - dNdEta += 1.f; - } + float dNdEta{0.f}; + computeDNDEta(dNdEta, mcParticles, histPath); - dNdEta /= (multEtaRange * 2.0f); uint32_t multiplicityCounter = 0; - getHist(TH1, histPath + "hLUTMultiplicity")->Fill(dNdEta); - // Now that the multiplicity is known, we can process the particles to smear them for (const auto& mcParticle : mcParticles) { const bool longLivedToBeHandled = std::find(longLivedHandledPDGs.begin(), longLivedHandledPDGs.end(), std::abs(mcParticle.pdgCode())) != longLivedHandledPDGs.end(); @@ -1962,9 +1966,6 @@ struct OnTheFlyTracker { multiplicityCounter++; o2::track::TrackParCov trackParCov; const bool isDecayDaughter = (mcParticle.getProcess() == TMCProcess::kPDecay); - const float timeResolutionNs = 100.f; // ns - const float nsToMus = 1e-3f; - const float timeResolutionUs = timeResolutionNs * nsToMus; // us const float time = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; bool reconstructed = false; @@ -2008,18 +2009,18 @@ struct OnTheFlyTracker { } if (reconstructed) { - tracksAlice3.push_back(TrackAlice3{trackParCov, mcParticle.globalIndex(), time, timeResolutionUs, isDecayDaughter, false, 0, nTrkHits}); + tracksAlice3.push_back(TrackAlice3{trackParCov, mcParticle.globalIndex(), time, timeResolutionUs, isDecayDaughter, false, 0, nTrkHits, kRecoPrimary}); } else { - ghostTracksAlice3.push_back(TrackAlice3{trackParCov, mcParticle.globalIndex(), time, timeResolutionUs, isDecayDaughter}); + ghostTracksAlice3.push_back(TrackAlice3{trackParCov, mcParticle.globalIndex(), time, timeResolutionUs, isDecayDaughter, false, 0, nTrkHits, kGhostPrimary}); } } o2::vertexing::PVertex primaryVertex; - if (enablePrimaryVertexing) { // disabled for now - primaryVertex.setXYZ(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ()); - } else { - primaryVertex.setXYZ(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ()); + if (recoPrimaries.size() <= 2) { + LOG(info) << "Not enough primary tracks (" << recoPrimaries.size() << ") to compute vertex, skipping vertexing and collision population."; + return; } + computeVertex(mcCollision, tracksAlice3, primaryVertex, icfg); getHist(TH1, histPath + "hSimMultiplicity")->Fill(multiplicityCounter); getHist(TH1, histPath + "hRecoMultiplicity")->Fill(tracksAlice3.size()); @@ -2039,95 +2040,8 @@ struct OnTheFlyTracker { // *+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+* // populate tracks - for (const auto& trackParCov : tracksAlice3) { - aod::track::TrackTypeEnum trackType = aod::track::Track; - - if (populateTracksDCA) { - float dcaXY = 1e+10, dcaZ = 1e+10; - o2::track::TrackParCov trackParametrization(trackParCov); - if (trackParametrization.propagateToDCA(primaryVertex, mMagneticField, &dcaInfo)) { - dcaXY = dcaInfo.getY(); - dcaZ = dcaInfo.getZ(); - } - - tableTracksDCA(dcaXY, dcaZ); - if (populateTracksDCACov) { - tableTracksDCACov(dcaInfo.getSigmaY2(), dcaInfo.getSigmaZ2()); - } - } - - tableStoredTracks(tableCollisions.lastIndex(), trackType, trackParCov.getX(), trackParCov.getAlpha(), trackParCov.getY(), trackParCov.getZ(), trackParCov.getSnp(), trackParCov.getTgl(), trackParCov.getQ2Pt()); - tableTracksExtension(trackParCov.getPt(), trackParCov.getP(), trackParCov.getEta(), trackParCov.getPhi()); - - // TODO do we keep the rho as 0? Also the sigma's are duplicated information - tableStoredTracksCov(std::sqrt(trackParCov.getSigmaY2()), std::sqrt(trackParCov.getSigmaZ2()), std::sqrt(trackParCov.getSigmaSnp2()), - std::sqrt(trackParCov.getSigmaTgl2()), std::sqrt(trackParCov.getSigma1Pt2()), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); - tableTracksCovExtension(trackParCov.getSigmaY2(), trackParCov.getSigmaZY(), trackParCov.getSigmaZ2(), trackParCov.getSigmaSnpY(), - trackParCov.getSigmaSnpZ(), trackParCov.getSigmaSnp2(), trackParCov.getSigmaTglY(), trackParCov.getSigmaTglZ(), trackParCov.getSigmaTglSnp(), - trackParCov.getSigmaTgl2(), trackParCov.getSigma1PtY(), trackParCov.getSigma1PtZ(), trackParCov.getSigma1PtSnp(), trackParCov.getSigma1PtTgl(), - trackParCov.getSigma1Pt2()); - tableMcTrackWithDauLabels(trackParCov.mcLabel, 0); - tableTracksExtraA3(trackParCov.nSiliconHits, trackParCov.nTPCHits); - - // populate extra tables if required to do so - if (populateTracksExtra) { - tableStoredTracksExtra(0.0f, static_cast(0), static_cast(0), static_cast(0), static_cast(0), - static_cast(0), static_cast(0), static_cast(0), static_cast(0), - 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, - 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f); - } - if (populateTrackSelection) { - tableTrackSelection(static_cast(0), false, false, false, false, false, false); - tableTrackSelectionExtension(false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false); - } - tableTracksAlice3(true); - } - - // *+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+* - // populate ghost tracks - for (const auto& trackParCov : ghostTracksAlice3) { - aod::track::TrackTypeEnum trackType = aod::track::Track; - - if (populateTracksDCA) { - float dcaXY = 1e+10, dcaZ = 1e+10; - o2::track::TrackParCov trackParametrization(trackParCov); - if (trackParametrization.propagateToDCA(primaryVertex, mMagneticField, &dcaInfo)) { - dcaXY = dcaInfo.getY(); - dcaZ = dcaInfo.getZ(); - } - - tableTracksDCA(dcaXY, dcaZ); - if (populateTracksDCACov) { - tableTracksDCACov(dcaInfo.getSigmaY2(), dcaInfo.getSigmaZ2()); - } - } - - tableStoredTracks(tableCollisions.lastIndex(), trackType, trackParCov.getX(), trackParCov.getAlpha(), trackParCov.getY(), trackParCov.getZ(), trackParCov.getSnp(), trackParCov.getTgl(), trackParCov.getQ2Pt()); - tableTracksExtension(trackParCov.getPt(), trackParCov.getP(), trackParCov.getEta(), trackParCov.getPhi()); - - // TODO do we keep the rho as 0? Also the sigma's are duplicated information - tableStoredTracksCov(std::sqrt(trackParCov.getSigmaY2()), std::sqrt(trackParCov.getSigmaZ2()), std::sqrt(trackParCov.getSigmaSnp2()), - std::sqrt(trackParCov.getSigmaTgl2()), std::sqrt(trackParCov.getSigma1Pt2()), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); - tableTracksCovExtension(trackParCov.getSigmaY2(), trackParCov.getSigmaZY(), trackParCov.getSigmaZ2(), trackParCov.getSigmaSnpY(), - trackParCov.getSigmaSnpZ(), trackParCov.getSigmaSnp2(), trackParCov.getSigmaTglY(), trackParCov.getSigmaTglZ(), trackParCov.getSigmaTglSnp(), - trackParCov.getSigmaTgl2(), trackParCov.getSigma1PtY(), trackParCov.getSigma1PtZ(), trackParCov.getSigma1PtSnp(), trackParCov.getSigma1PtTgl(), - trackParCov.getSigma1Pt2()); - tableMcTrackWithDauLabels(trackParCov.mcLabel, 0); - tableTracksExtraA3(trackParCov.nSiliconHits, trackParCov.nTPCHits); - - // populate extra tables if required to do so - if (populateTracksExtra) { - tableStoredTracksExtra(0.0f, static_cast(0), static_cast(0), static_cast(0), static_cast(0), - static_cast(0), static_cast(0), static_cast(0), static_cast(0), - 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, - 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f); - } - if (populateTrackSelection) { - tableTrackSelection(static_cast(0), false, false, false, false, false, false); - tableTrackSelectionExtension(false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false, false); - } - tableTracksAlice3(false); - } + fillTracksInfo(tracksAlice3, primaryVertex, icfg); + fillTracksInfo(ghostTracksAlice3, primaryVertex, icfg); } void processDecayer(aod::McCollision const& mcCollision, aod::McPartWithDaus const& mcParticles) @@ -2152,4 +2066,4 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) return WorkflowSpec{ adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc)}; -} +} \ No newline at end of file diff --git a/ALICE3/TableProducer/alice3-decayfinder.cxx b/ALICE3/TableProducer/alice3-decayfinder.cxx index f828ee5fbff..b61ddb1a421 100644 --- a/ALICE3/TableProducer/alice3-decayfinder.cxx +++ b/ALICE3/TableProducer/alice3-decayfinder.cxx @@ -79,6 +79,7 @@ struct alice3decayFinder { Produces mcRecFlags; // contains MC truth info (is true Lc, or bkg) Produces pidInfoLcDaugs; // contains PID info for Lc candidates Produces mcGenFlags; // contains MC gen info for 3-prong candidates + Produces coordsSVResos; // contains MC gen and reco coordinates of the secondary vertex info for 3-prong candidates // Vertexing struct : ConfigurableGroup { @@ -231,6 +232,7 @@ struct alice3decayFinder { std::array Pdaug2; // pion track std::array primaryVertex; // primary vertex coordinates std::array secondaryVertex; // secondary vertex coordinates + std::array deltaRecoGenSV; // difference between reconstructed and generated secondary vertex coordinates float impactParameterY0; // impact parameters float errorImpactParameterY0; // impact parameters error float impactParameterY1; // impact parameters @@ -434,6 +436,9 @@ struct alice3decayFinder { mCandidate3Prong.origin = 0; // Default: unknown origin if (mCandidate3Prong.particleMcRec > -1) { const auto& motherParticle = mcParticles.rawIteratorAt(mCandidate3Prong.particleMcRec); + mCandidate3Prong.deltaRecoGenSV[0] = mCandidate3Prong.secondaryVertex[0] - motherParticle.vx(); + mCandidate3Prong.deltaRecoGenSV[1] = mCandidate3Prong.secondaryVertex[1] - motherParticle.vy(); + mCandidate3Prong.deltaRecoGenSV[2] = mCandidate3Prong.secondaryVertex[2] - motherParticle.vz(); mCandidate3Prong.flagMc = motherParticle.pdgCode() > 0 ? charmHadFlag : -charmHadFlag; // Particle std::vector idxBhadMothers{}; mCandidate3Prong.origin = RecoDecay::getCharmHadronOrigin(mcParticles, motherParticle, false, &idxBhadMothers); @@ -1054,8 +1059,9 @@ struct alice3decayFinder { std::sqrt(mCandidate3Prong.errorImpactParameterZ1), std::sqrt(mCandidate3Prong.errorImpactParameterZ2), candPx, candPy, candPz); - mcRecFlags(mCandidate3Prong.origin, mCandidate3Prong.ptBMotherRec, mCandidate3Prong.flagMc, mCandidate3Prong.particleMcRec); // placeholder for prompt/non-prompt + mcRecFlags(mCandidate3Prong.origin, mCandidate3Prong.ptBMotherRec, mCandidate3Prong.flagMc, mCandidate3Prong.particleMcRec); fillPidTable(prong0, prong1, prong2); + coordsSVResos(mCandidate3Prong.deltaRecoGenSV[0], mCandidate3Prong.deltaRecoGenSV[1], mCandidate3Prong.deltaRecoGenSV[2]); } } } diff --git a/ALICE3/TableProducer/alice3TrackingTranslator.cxx b/ALICE3/TableProducer/alice3TrackingTranslator.cxx index c68cb2adbe3..d9b2d3874f6 100644 --- a/ALICE3/TableProducer/alice3TrackingTranslator.cxx +++ b/ALICE3/TableProducer/alice3TrackingTranslator.cxx @@ -582,7 +582,8 @@ struct Alice3TrackingTranslator { tableTracksAlice3Pdg(pdgCode); // PdgCode to the linked MC truth particle tableTracksExtraA3(m_nMeasurements, // nSiliconHits (using m_nMeasurements as proxy) - 0); // nTPCHits + 0, // nTPCHits + 0); // trackType // Fill extra track info tableStoredTracksExtra(0.f, // TPCInnerParam diff --git a/ALICE3/Tasks/alice3HfTask3Prong.cxx b/ALICE3/Tasks/alice3HfTask3Prong.cxx index 4f21bada994..ead69052c9e 100644 --- a/ALICE3/Tasks/alice3HfTask3Prong.cxx +++ b/ALICE3/Tasks/alice3HfTask3Prong.cxx @@ -64,8 +64,8 @@ struct Alice3HfTask3Prong { int selectedPdg{-1}; - using Cands3PReco = soa::Filtered>; - using Cands3PRecoWMl = soa::Filtered>; + using Cands3PReco = soa::Filtered>; + using Cands3PRecoWMl = soa::Filtered>; using Cands3PGen = soa::Join; Filter filterSelectCandidates = (aod::a3_hf_sel_3prong::isSelMassHypo0 == true || aod::a3_hf_sel_3prong::isSelMassHypo1 == true); @@ -79,6 +79,7 @@ struct Alice3HfTask3Prong { ConfigurableAxis thnConfigAxisCanType{"thnConfigAxisCanType", {5, 0., 5.}, ""}; ConfigurableAxis thnAxisRapidity{"thnAxisRapidity", {20, -1, 1}, "Cand. rapidity bins"}; ConfigurableAxis thnConfigAxisGenPtB{"thnConfigAxisGenPtB", {1000, 0, 100}, "Gen Pt B"}; + ConfigurableAxis thnConfigAxisDeltaSecVtxCoord{"thnConfigAxisDeltaSecVtxCoord", {1000, -10, 10}, ""}; HistogramRegistry registry{"registry", {}}; @@ -142,6 +143,9 @@ struct Alice3HfTask3Prong { addHistogramsRec("hCPAxy", "cosine of pointing angle xy", "entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}); addHistogramsRec("hDca2", "prong Chi2PCA to sec. vertex (cm)", "entries", {HistType::kTH1F, {{400, 0., 20.}}}); addHistogramsRec("hEta", "#it{#eta}", "entries", {HistType::kTH1F, {{100, -2., 2.}}}); + addHistogramsRec("hDeltaXRecoGenSecVtx", "p_{T} (GeV/#it{c})", "#Delta X Sec. Vtx. (cm)", {HistType::kTH2F, {thnConfigAxisPt, thnConfigAxisDeltaSecVtxCoord}}); + addHistogramsRec("hDeltaYRecoGenSecVtx", "p_{T} (GeV/#it{c})", "#Delta Y Sec. Vtx. (cm)", {HistType::kTH2F, {thnConfigAxisPt, thnConfigAxisDeltaSecVtxCoord}}); + addHistogramsRec("hDeltaZRecoGenSecVtx", "p_{T} (GeV/#it{c})", "#Delta Z Sec. Vtx. (cm)", {HistType::kTH2F, {thnConfigAxisPt, thnConfigAxisDeltaSecVtxCoord}}); addHistogramsRec("hd0VsPtProng0", "prong 0 DCAxy to prim. vertex (cm)", "#it{p}_{T} (GeV/#it{c})", {HistType::kTH2F, {{600, -0.4, 0.4}, {vbins}}}); addHistogramsRec("hd0VsPtProng1", "prong 1 DCAxy to prim. vertex (cm)", "#it{p}_{T} (GeV/#it{c})", {HistType::kTH2F, {{600, -0.4, 0.4}, {vbins}}}); addHistogramsRec("hd0VsPtProng2", "prong 2 DCAxy to prim. vertex (cm)", "#it{p}_{T} (GeV/#it{c})", {HistType::kTH2F, {{600, -0.4, 0.4}, {vbins}}}); @@ -257,6 +261,9 @@ struct Alice3HfTask3Prong { registry.fill(histoPrefix + HIST("hImpParErrProng1VsPt") + histoSuffix, candidate.errorImpactParameterY1(), candidate.pt()); registry.fill(histoPrefix + HIST("hImpParErrProng2VsPt") + histoSuffix, candidate.errorImpactParameterY2(), candidate.pt()); registry.fill(histoPrefix + HIST("hDecLenErrVsPt") + histoSuffix, candidate.errorDecayLength(), candidate.pt()); + registry.fill(histoPrefix + HIST("hDeltaXRecoGenSecVtx") + histoSuffix, candidate.pt(), candidate.deltaXSecVtx()); + registry.fill(histoPrefix + HIST("hDeltaYRecoGenSecVtx") + histoSuffix, candidate.pt(), candidate.deltaYSecVtx()); + registry.fill(histoPrefix + HIST("hDeltaZRecoGenSecVtx") + histoSuffix, candidate.pt(), candidate.deltaZSecVtx()); } /// Fill MC histograms at reconstruction level @@ -264,8 +271,8 @@ struct Alice3HfTask3Prong { /// \tparam SaveMl indicates whether ML scores are saved in the THnSparse /// \tparam CandsRec is the type of the reconstructed candidates collection /// \param candidates is the collection of reconstructed candidates - template - void fillHistosMcRec(CandsRec const& candidates, AllParticles const& allParticles) + template + void fillHistosMcRec(CandsRec const& candidates, GenParticles const& genParticles) { registry.fill(HIST("MC/rec/hCandidateCounter"), 0.); for (const auto& candidate : candidates) { @@ -277,13 +284,13 @@ struct Alice3HfTask3Prong { registry.fill(HIST("MC/rec/hCandidateCounter"), 2.); if (candidate.particleMcRec() >= 0) { registry.fill(HIST("MC/rec/hCandidateCounter"), 3.); - auto mcParticle = allParticles.iteratorAt(candidate.particleMcRec()); + auto mcParticle = genParticles.iteratorAt(candidate.particleMcRec()); if (mcParticle.has_daughters()) { auto daughters = mcParticle.daughtersIds(); LOG(debug) << "Reco candidate matched to MC particle with PDG " << mcParticle.pdgCode() << " daughters: " << daughters.size(); int prongIdx = 0; for (int dauId = daughters[0]; dauId <= daughters[1]; dauId++) { - auto dau = allParticles.iteratorAt(dauId); + auto dau = genParticles.iteratorAt(dauId); LOG(debug) << " dauId: " << dauId << " PDG: " << dau.pdgCode() << " with pT: " << dau.pt(); switch (prongIdx) { case 0: @@ -376,9 +383,9 @@ struct Alice3HfTask3Prong { /// \tparam ParticleType is the type of the generated particle /// \tparam TableType is the type of the full table (non-Partition) /// \param particle is a generated particle - /// \param allParticles is the full table of particles for iteratorAt access + /// \param genParticles is the full table of particles for iteratorAt access template - void fillHistogramsGen(ParticleType const& particle, TableType const& allParticles) + void fillHistogramsGen(ParticleType const& particle, TableType const& genParticles) { LOG(debug) << "Filling generated histograms for signal type " << SignalType; static constexpr auto histoPrefix = HIST("MC/gen/") + HIST(SignalFolders[SignalType]) + HIST("/"); @@ -400,7 +407,7 @@ struct Alice3HfTask3Prong { float pz = 0.f; float e = 0.f; for (int iDau = firstDauIdx; iDau <= lastDauIdx && iDau > 0; iDau++) { - const auto& dau = allParticles.iteratorAt(iDau); + const auto& dau = genParticles.iteratorAt(iDau); e += dau.e(); px += dau.px(); py += dau.py(); @@ -415,11 +422,11 @@ struct Alice3HfTask3Prong { /// Fill MC histograms at generated level /// \tparam CharmHad is the charm hadron species /// \tparam CandsGen is the type of the generated candidates collection - /// \tparam AllParticles is the type of the full particle table + /// \tparam GenParticles is the type of the full particle table /// \param mcParticles is the collection of generated particles (can be a Partition) - /// \param allParticles is the full table of particles - template - void fillHistosMcGen(CandsGen const& mcParticles, AllParticles const& allParticles) + /// \param genParticles is the full table of particles + template + void fillHistosMcGen(CandsGen const& mcParticles, GenParticles const& genParticles) { // MC gen. for (const auto& particle : mcParticles) { @@ -431,14 +438,14 @@ struct Alice3HfTask3Prong { const auto ptGen = particle.pt(); const auto originType = particle.originMcGen(); - fillHistogramsGen(particle, allParticles); + fillHistogramsGen(particle, genParticles); float ptGenB = -1.f; if (originType == RecoDecay::OriginType::Prompt) { - fillHistogramsGen(particle, allParticles); + fillHistogramsGen(particle, genParticles); } else if (particle.originMcGen() == RecoDecay::OriginType::NonPrompt) { ptGenB = particle.bHadMotherPtGen(); - fillHistogramsGen(particle, allParticles); + fillHistogramsGen(particle, genParticles); } if (fillThn) { diff --git a/Common/TableProducer/qVectorsTable.cxx b/Common/TableProducer/qVectorsTable.cxx index ce5152aace6..677453eb6ca 100644 --- a/Common/TableProducer/qVectorsTable.cxx +++ b/Common/TableProducer/qVectorsTable.cxx @@ -590,15 +590,15 @@ struct qVectorsTable { float sumAmplFT0C = 0.; float sumAmplFT0M = 0.; float sumAmplFV0A = 0.; - + if (coll.has_foundFT0() && (useDetector["QvectorFT0As"] || useDetector["QvectorFT0Cs"] || useDetector["QvectorFT0Ms"])) { auto ft0 = coll.foundFT0(); - + if (useDetector["QvectorFT0As"]) { for (std::size_t iChA = 0; iChA < ft0.channelA().size(); iChA++) { float ampl = ft0.amplitudeA()[iChA]; int FT0AchId = ft0.channelA()[iChA]; - + histosQA.fill(HIST("FT0Amp"), ampl, FT0AchId); histosQA.fill(HIST("FT0AmpCor"), ampl / FT0RelGainConst[FT0AchId], FT0AchId); diff --git a/PWGJE/Tasks/CMakeLists.txt b/PWGJE/Tasks/CMakeLists.txt index a1b103a528c..63e8fe99f69 100644 --- a/PWGJE/Tasks/CMakeLists.txt +++ b/PWGJE/Tasks/CMakeLists.txt @@ -188,10 +188,10 @@ if(FastJet_FOUND) SOURCES v0QA.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore COMPONENT_NAME Analysis) - o2physics_add_dpl_workflow(jet-finder-charged-qa - SOURCES jetFinderQA.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore - COMPONENT_NAME Analysis) +# o2physics_add_dpl_workflow(jet-finder-charged-qa +# SOURCES jetFinderQA.cxx +# PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore +# COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(jet-outlier-qa SOURCES jetOutlierQA.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore @@ -366,11 +366,11 @@ if(FastJet_FOUND) SOURCES dijetFinderQA.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore COMPONENT_NAME Analysis) - o2physics_add_dpl_workflow(bjet-tagging-gnn - SOURCES bjetTaggingGnn.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore O2Physics::EventFilteringUtils - O2Physics::AnalysisCCDB - COMPONENT_NAME Analysis) + # o2physics_add_dpl_workflow(bjet-tagging-gnn + # SOURCES bjetTaggingGnn.cxx + # PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore O2Physics::EventFilteringUtils + # O2Physics::AnalysisCCDB + # COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(jet-shape SOURCES jetShape.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore From 3148b35a065613b31b0fd7293de4f2fd2d459135 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Tue, 12 May 2026 06:14:35 +0000 Subject: [PATCH 2/5] Please consider the following formatting changes --- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 98 ++++++++++++-------- Common/TableProducer/qVectorsTable.cxx | 6 +- 2 files changed, 60 insertions(+), 44 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index 428d5d4bf33..0cede1f5fa0 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -275,16 +275,15 @@ struct OnTheFlyTracker { int isUsedInCascadingInput = 0, int nSiliconHitsInput = 0, int nTPCHitsInput = 0, - TrackType trackTypeInput = TrackType::kRecoPrimary - ) : o2::track::TrackParCov(src), - mcLabel{label}, - timeEst{time, timeError}, - isDecayDau(decayDauInput), - isWeakDecayDau(weakDecayDauInput), - isUsedInCascading(isUsedInCascadingInput), - nSiliconHits(nSiliconHitsInput), - nTPCHits(nTPCHitsInput), - trackType(trackTypeInput) {} + TrackType trackTypeInput = TrackType::kRecoPrimary) : o2::track::TrackParCov(src), + mcLabel{label}, + timeEst{time, timeError}, + isDecayDau(decayDauInput), + isWeakDecayDau(weakDecayDauInput), + isUsedInCascading(isUsedInCascadingInput), + nSiliconHits(nSiliconHitsInput), + nTPCHits(nTPCHitsInput), + trackType(trackTypeInput) {} const TimeEst& getTimeMUS() const { return timeEst; } int64_t mcLabel; TimeEst timeEst; ///< time estimate in ns @@ -842,7 +841,8 @@ struct OnTheFlyTracker { /// \param mcParticles the set of MC particles to compute dN/deta from /// \param histPath the path to the histogram where the computed dN/deta value will be stored for QA purposes template - void computeDNDEta(float& dNdEta, McParticleType const& mcParticles, const std::string histPath) { + void computeDNDEta(float& dNdEta, McParticleType const& mcParticles, const std::string histPath) + { for (const auto& mcParticle : mcParticles) { if (std::abs(mcParticle.eta()) > multEtaRange) { continue; @@ -874,7 +874,6 @@ struct OnTheFlyTracker { dNdEta /= (multEtaRange * 2.0f); getHist(TH1, histPath + "hLUTMultiplicity")->Fill(dNdEta); - } /// Function to study the cascade decay and fill the relevant histograms and output track vector @@ -981,8 +980,7 @@ struct OnTheFlyTracker { isReco[i] = true; for (uint32_t ih = 0; ih < fastTracker[icfg]->GetNHits() && cascadeDecaySettings.doXiQA; ih++) { - getHist(TH2, histPath + "hFastTrackerHits")->Fill(fastTracker[icfg]->GetHitZ(ih), - std::hypot(fastTracker[icfg]->GetHitX(ih), fastTracker[icfg]->GetHitY(ih))); + getHist(TH2, histPath + "hFastTrackerHits")->Fill(fastTracker[icfg]->GetHitZ(ih), std::hypot(fastTracker[icfg]->GetHitX(ih), fastTracker[icfg]->GetHitY(ih))); } } else { isReco[i] = true; @@ -1024,8 +1022,8 @@ struct OnTheFlyTracker { // n-2: pion from lambda // n-3: pion from xi thisCascade.bachelorId = trackTableOffset + 1; - thisCascade.positiveId = trackTableOffset + 2; // Lambda daughters - thisCascade.negativeId = trackTableOffset + 3; // Lambda daughters + thisCascade.positiveId = trackTableOffset + 2; // Lambda daughters + thisCascade.negativeId = trackTableOffset + 3; // Lambda daughters thisCascade.cascadeTrackId = trackTableOffset + 4; // use DCA fitters @@ -1237,8 +1235,12 @@ struct OnTheFlyTracker { } // Reset indices - if (!isReco[1]) { thisCascade.negativeId = -1; } - if (!isReco[2]) { thisCascade.positiveId = -1; } + if (!isReco[1]) { + thisCascade.negativeId = -1; + } + if (!isReco[2]) { + thisCascade.positiveId = -1; + } int nCand = 0; bool kinkFitterOK = true; @@ -1267,8 +1269,8 @@ struct OnTheFlyTracker { o2::track::TrackParCov newCascadeTrack = fitter.getTrack(0); // (cascade) std::array kinkVtx = {-999, -999, -999}; kinkVtx = fitter.getPCACandidatePos(); - thisCascade.dcaV0dau = -1.f; // unknown - thisCascade.v0radius = -1.f; // unknown + thisCascade.dcaV0dau = -1.f; // unknown + thisCascade.v0radius = -1.f; // unknown thisCascade.dcacascdau = std::sqrt(fitter.getChi2AtPCACandidate()); thisCascade.cascradius = std::hypot(kinkVtx[0], kinkVtx[1]); thisCascade.cascradiusMC = xiDecayRadius2D; @@ -1319,7 +1321,7 @@ struct OnTheFlyTracker { if (isReco[1]) { getHist(TH2, histPath + "hRecoPiFromLa")->Fill(laDecayRadius2D, cascadeDecayProducts[1].Pt()); o2::track::TrackParCov trackParametrizationCascProng1(tracksCascadeProngs[1]); - if (populateTracksDCA && tracksCascadeProngs[1].propagateToDCA(primaryVertex, mMagneticField, &dcaInfo)) { // FIXME: this is not the right trackParametrization, need to propagate the negative pion track + if (populateTracksDCA && tracksCascadeProngs[1].propagateToDCA(primaryVertex, mMagneticField, &dcaInfo)) { // FIXME: this is not the right trackParametrization, need to propagate the negative pion track dcaXY = dcaInfo.getY(); dcaZ = dcaInfo.getZ(); getHist(TH2, histPath + "h2dDCAxyCascadeNegative")->Fill(trackParametrizationCascProng1.getPt(), dcaXY * 1e+4); // in microns, please @@ -1329,7 +1331,7 @@ struct OnTheFlyTracker { if (isReco[2]) { getHist(TH2, histPath + "hRecoPrFromLa")->Fill(laDecayRadius2D, cascadeDecayProducts[2].Pt()); o2::track::TrackParCov trackParametrizationCascProng2(tracksCascadeProngs[2]); - if (populateTracksDCA && tracksCascadeProngs[2].propagateToDCA(primaryVertex, mMagneticField, &dcaInfo)) { // FIXME: this is not the right trackParametrization, need to propagate the positive proton track + if (populateTracksDCA && tracksCascadeProngs[2].propagateToDCA(primaryVertex, mMagneticField, &dcaInfo)) { // FIXME: this is not the right trackParametrization, need to propagate the positive proton track dcaXY = dcaInfo.getY(); dcaZ = dcaInfo.getZ(); getHist(TH2, histPath + "h2dDCAxyCascadePositive")->Fill(trackParametrizationCascProng2.getPt(), dcaXY * 1e+4); // in microns, please @@ -1365,10 +1367,11 @@ struct OnTheFlyTracker { template void studyV0(const int trackTableOffset, const McParticleType& mcParticle, - std::vector& tracksV0Daugs, + std::vector& tracksV0Daugs, int icfg, int dNdEta, - float eventCollisionTimeNS) { + float eventCollisionTimeNS) + { o2::track::TrackParCov trackParCov; o2::upgrade::convertMCParticleToO2Track(mcParticle, trackParCov, pdgDB); @@ -1439,9 +1442,12 @@ struct OnTheFlyTracker { nV0SiliconHits[i] = fastTracker[icfg]->GetNSiliconPoints(); nV0TPCHits[i] = fastTracker[icfg]->GetNGasPoints(); - if (nV0SiliconHits[i] < fastTrackerSettings.minSiliconHits) continue; - if (nV0SiliconHits[i] < fastTrackerSettings.minSiliconHitsIfTPCUsed) continue; - if (nV0TPCHits[i] < fastTrackerSettings.minTPCClusters) continue; + if (nV0SiliconHits[i] < fastTrackerSettings.minSiliconHits) + continue; + if (nV0SiliconHits[i] < fastTrackerSettings.minSiliconHitsIfTPCUsed) + continue; + if (nV0TPCHits[i] < fastTrackerSettings.minTPCClusters) + continue; isV0Reco[i] = true; if (v0DecaySettings.doV0QA) { // QA @@ -1607,7 +1613,8 @@ struct OnTheFlyTracker { /// \param primaryVertex the output variable where the computed primary vertex will be stored /// \param icfg index of the current configuration, used for histogram filling template - void computeVertex(McCollisionType& mcCollision, const std::vector& prmTrks, o2::vertexing::PVertex& primaryVertex, const int icfg) { + void computeVertex(McCollisionType& mcCollision, const std::vector& prmTrks, o2::vertexing::PVertex& primaryVertex, const int icfg) + { if (!enablePrimaryVertexing) { primaryVertex.setXYZ(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ()); @@ -1673,7 +1680,8 @@ struct OnTheFlyTracker { /// \param tracks the vector of tracks to be processed /// \param primaryVertex the primary vertex position, used for DCA calculation /// \param icfg index of the current configuration, used for histogram filling - void fillTracksInfo(std::vector const& tracks, o2::vertexing::PVertex const& primaryVertex, const int icfg) { + void fillTracksInfo(std::vector const& tracks, o2::vertexing::PVertex const& primaryVertex, const int icfg) + { const std::string histPath = "Configuration_" + std::to_string(icfg) + "/"; for (const auto& trackParCov : tracks) { @@ -1765,8 +1773,8 @@ struct OnTheFlyTracker { const bool longLivedToBeHandled = std::find(longLivedHandledPDGs.begin(), longLivedHandledPDGs.end(), std::abs(mcParticle.pdgCode())) != longLivedHandledPDGs.end(); const bool nucleiToBeHandled = std::find(nucleiPDGs.begin(), nucleiPDGs.end(), std::abs(mcParticle.pdgCode())) != nucleiPDGs.end(); - const bool pdgsToBeHandled = longLivedToBeHandled || - (enableNucleiSmearing && nucleiToBeHandled) || + const bool pdgsToBeHandled = longLivedToBeHandled || + (enableNucleiSmearing && nucleiToBeHandled) || (cascadeDecaySettings.decayXi && isCascade) || (v0DecaySettings.decayV0 && isV0); if (!pdgsToBeHandled) { @@ -1799,10 +1807,14 @@ struct OnTheFlyTracker { } getHist(TH1, histPath + "hPtGenerated")->Fill(mcParticle.pt()); getHist(TH1, histPath + "hPhiGenerated")->Fill(mcParticle.phi()); - if (std::abs(mcParticle.pdgCode()) == kElectron) getHist(TH1, histPath + "hPtGeneratedEl")->Fill(mcParticle.pt()); - if (std::abs(mcParticle.pdgCode()) == kPiPlus) getHist(TH1, histPath + "hPtGeneratedPi")->Fill(mcParticle.pt()); - if (std::abs(mcParticle.pdgCode()) == kKPlus) getHist(TH1, histPath + "hPtGeneratedKa")->Fill(mcParticle.pt()); - if (std::abs(mcParticle.pdgCode()) == kProton) getHist(TH1, histPath + "hPtGeneratedPr")->Fill(mcParticle.pt()); + if (std::abs(mcParticle.pdgCode()) == kElectron) + getHist(TH1, histPath + "hPtGeneratedEl")->Fill(mcParticle.pt()); + if (std::abs(mcParticle.pdgCode()) == kPiPlus) + getHist(TH1, histPath + "hPtGeneratedPi")->Fill(mcParticle.pt()); + if (std::abs(mcParticle.pdgCode()) == kKPlus) + getHist(TH1, histPath + "hPtGeneratedKa")->Fill(mcParticle.pt()); + if (std::abs(mcParticle.pdgCode()) == kProton) + getHist(TH1, histPath + "hPtGeneratedPr")->Fill(mcParticle.pt()); if (!reconstructed && !processUnreconstructedTracks) { continue; @@ -1810,10 +1822,14 @@ struct OnTheFlyTracker { // LOG(info) << "Particle with PDG " << mcParticle.pdgCode() << (reconstructed ? " was reconstructed." : " was not reconstructed."); getHist(TH1, histPath + "hPtReconstructed")->Fill(trackParCov.getPt()); - if (std::abs(mcParticle.pdgCode()) == kElectron) getHist(TH1, histPath + "hPtReconstructedEl")->Fill(trackParCov.getPt()); - if (std::abs(mcParticle.pdgCode()) == kPiPlus) getHist(TH1, histPath + "hPtReconstructedPi")->Fill(trackParCov.getPt()); - if (std::abs(mcParticle.pdgCode()) == kKPlus) getHist(TH1, histPath + "hPtReconstructedKa")->Fill(trackParCov.getPt()); - if (std::abs(mcParticle.pdgCode()) == kProton) getHist(TH1, histPath + "hPtReconstructedPr")->Fill(trackParCov.getPt()); + if (std::abs(mcParticle.pdgCode()) == kElectron) + getHist(TH1, histPath + "hPtReconstructedEl")->Fill(trackParCov.getPt()); + if (std::abs(mcParticle.pdgCode()) == kPiPlus) + getHist(TH1, histPath + "hPtReconstructedPi")->Fill(trackParCov.getPt()); + if (std::abs(mcParticle.pdgCode()) == kKPlus) + getHist(TH1, histPath + "hPtReconstructedKa")->Fill(trackParCov.getPt()); + if (std::abs(mcParticle.pdgCode()) == kProton) + getHist(TH1, histPath + "hPtReconstructedPr")->Fill(trackParCov.getPt()); } if (doExtraQA) { getHist(TH2, histPath + "h2dPtRes")->Fill(trackParCov.getPt(), (trackParCov.getPt() - mcParticle.pt()) / trackParCov.getPt()); @@ -1881,7 +1897,7 @@ struct OnTheFlyTracker { // *+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+* // populate tracks - fillTracksInfo(recoPrimaries, primaryVertex, icfg); + fillTracksInfo(recoPrimaries, primaryVertex, icfg); fillTracksInfo(ghostPrimaries, primaryVertex, icfg); // do bookkeeping of fastTracker tracking @@ -2066,4 +2082,4 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) return WorkflowSpec{ adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc)}; -} \ No newline at end of file +} diff --git a/Common/TableProducer/qVectorsTable.cxx b/Common/TableProducer/qVectorsTable.cxx index 677453eb6ca..ce5152aace6 100644 --- a/Common/TableProducer/qVectorsTable.cxx +++ b/Common/TableProducer/qVectorsTable.cxx @@ -590,15 +590,15 @@ struct qVectorsTable { float sumAmplFT0C = 0.; float sumAmplFT0M = 0.; float sumAmplFV0A = 0.; - + if (coll.has_foundFT0() && (useDetector["QvectorFT0As"] || useDetector["QvectorFT0Cs"] || useDetector["QvectorFT0Ms"])) { auto ft0 = coll.foundFT0(); - + if (useDetector["QvectorFT0As"]) { for (std::size_t iChA = 0; iChA < ft0.channelA().size(); iChA++) { float ampl = ft0.amplitudeA()[iChA]; int FT0AchId = ft0.channelA()[iChA]; - + histosQA.fill(HIST("FT0Amp"), ampl, FT0AchId); histosQA.fill(HIST("FT0AmpCor"), ampl / FT0RelGainConst[FT0AchId], FT0AchId); From d87a246cf9bf10e054dd87a847bd6f259c7b0373 Mon Sep 17 00:00:00 2001 From: marcellocosti Date: Tue, 12 May 2026 08:20:52 +0200 Subject: [PATCH 3/5] Restore unwanted changes --- ALICE3/DataModel/A3DecayFinderTables.h | 10 ------ ALICE3/TableProducer/alice3-decayfinder.cxx | 8 +---- ALICE3/Tasks/alice3HfTask3Prong.cxx | 39 +++++++++------------ PWGJE/Tasks/CMakeLists.txt | 18 +++++----- 4 files changed, 26 insertions(+), 49 deletions(-) diff --git a/ALICE3/DataModel/A3DecayFinderTables.h b/ALICE3/DataModel/A3DecayFinderTables.h index 2ed4796045c..bd9aa5592b3 100644 --- a/ALICE3/DataModel/A3DecayFinderTables.h +++ b/ALICE3/DataModel/A3DecayFinderTables.h @@ -202,11 +202,6 @@ DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterXY, impactParameterXY, //! DECLARE_SOA_COLUMN(MlScore0, mlScore0, float); //! DECLARE_SOA_COLUMN(MlScore1, mlScore1, float); //! DECLARE_SOA_COLUMN(MlScore2, mlScore2, float); //! - -// Secondary vertex gen-reco position -DECLARE_SOA_COLUMN(DeltaXSecVtx, deltaXSecVtx, float); //! -DECLARE_SOA_COLUMN(DeltaYSecVtx, deltaYSecVtx, float); //! -DECLARE_SOA_COLUMN(DeltaZSecVtx, deltaZSecVtx, float); //! } // namespace a3_hf_cand #define HFCAND_COLUMNS \ @@ -240,11 +235,6 @@ DECLARE_SOA_COLUMN(DeltaZSecVtx, deltaZSecVtx, float); //! a3_hf_cand::Pt2Prong2, \ a3_hf_cand::PVectorProng2 -DECLARE_SOA_TABLE(Alice3SVResos, "AOD", "ALICE3SVRESO", //! - a3_hf_cand::DeltaXSecVtx, - a3_hf_cand::DeltaYSecVtx, - a3_hf_cand::DeltaZSecVtx); - namespace a3DecayMap { DECLARE_SOA_COLUMN(DecayMap, decayMap, uint32_t); //! simple map to process passing / not passing criteria diff --git a/ALICE3/TableProducer/alice3-decayfinder.cxx b/ALICE3/TableProducer/alice3-decayfinder.cxx index b61ddb1a421..f828ee5fbff 100644 --- a/ALICE3/TableProducer/alice3-decayfinder.cxx +++ b/ALICE3/TableProducer/alice3-decayfinder.cxx @@ -79,7 +79,6 @@ struct alice3decayFinder { Produces mcRecFlags; // contains MC truth info (is true Lc, or bkg) Produces pidInfoLcDaugs; // contains PID info for Lc candidates Produces mcGenFlags; // contains MC gen info for 3-prong candidates - Produces coordsSVResos; // contains MC gen and reco coordinates of the secondary vertex info for 3-prong candidates // Vertexing struct : ConfigurableGroup { @@ -232,7 +231,6 @@ struct alice3decayFinder { std::array Pdaug2; // pion track std::array primaryVertex; // primary vertex coordinates std::array secondaryVertex; // secondary vertex coordinates - std::array deltaRecoGenSV; // difference between reconstructed and generated secondary vertex coordinates float impactParameterY0; // impact parameters float errorImpactParameterY0; // impact parameters error float impactParameterY1; // impact parameters @@ -436,9 +434,6 @@ struct alice3decayFinder { mCandidate3Prong.origin = 0; // Default: unknown origin if (mCandidate3Prong.particleMcRec > -1) { const auto& motherParticle = mcParticles.rawIteratorAt(mCandidate3Prong.particleMcRec); - mCandidate3Prong.deltaRecoGenSV[0] = mCandidate3Prong.secondaryVertex[0] - motherParticle.vx(); - mCandidate3Prong.deltaRecoGenSV[1] = mCandidate3Prong.secondaryVertex[1] - motherParticle.vy(); - mCandidate3Prong.deltaRecoGenSV[2] = mCandidate3Prong.secondaryVertex[2] - motherParticle.vz(); mCandidate3Prong.flagMc = motherParticle.pdgCode() > 0 ? charmHadFlag : -charmHadFlag; // Particle std::vector idxBhadMothers{}; mCandidate3Prong.origin = RecoDecay::getCharmHadronOrigin(mcParticles, motherParticle, false, &idxBhadMothers); @@ -1059,9 +1054,8 @@ struct alice3decayFinder { std::sqrt(mCandidate3Prong.errorImpactParameterZ1), std::sqrt(mCandidate3Prong.errorImpactParameterZ2), candPx, candPy, candPz); - mcRecFlags(mCandidate3Prong.origin, mCandidate3Prong.ptBMotherRec, mCandidate3Prong.flagMc, mCandidate3Prong.particleMcRec); + mcRecFlags(mCandidate3Prong.origin, mCandidate3Prong.ptBMotherRec, mCandidate3Prong.flagMc, mCandidate3Prong.particleMcRec); // placeholder for prompt/non-prompt fillPidTable(prong0, prong1, prong2); - coordsSVResos(mCandidate3Prong.deltaRecoGenSV[0], mCandidate3Prong.deltaRecoGenSV[1], mCandidate3Prong.deltaRecoGenSV[2]); } } } diff --git a/ALICE3/Tasks/alice3HfTask3Prong.cxx b/ALICE3/Tasks/alice3HfTask3Prong.cxx index ead69052c9e..4f21bada994 100644 --- a/ALICE3/Tasks/alice3HfTask3Prong.cxx +++ b/ALICE3/Tasks/alice3HfTask3Prong.cxx @@ -64,8 +64,8 @@ struct Alice3HfTask3Prong { int selectedPdg{-1}; - using Cands3PReco = soa::Filtered>; - using Cands3PRecoWMl = soa::Filtered>; + using Cands3PReco = soa::Filtered>; + using Cands3PRecoWMl = soa::Filtered>; using Cands3PGen = soa::Join; Filter filterSelectCandidates = (aod::a3_hf_sel_3prong::isSelMassHypo0 == true || aod::a3_hf_sel_3prong::isSelMassHypo1 == true); @@ -79,7 +79,6 @@ struct Alice3HfTask3Prong { ConfigurableAxis thnConfigAxisCanType{"thnConfigAxisCanType", {5, 0., 5.}, ""}; ConfigurableAxis thnAxisRapidity{"thnAxisRapidity", {20, -1, 1}, "Cand. rapidity bins"}; ConfigurableAxis thnConfigAxisGenPtB{"thnConfigAxisGenPtB", {1000, 0, 100}, "Gen Pt B"}; - ConfigurableAxis thnConfigAxisDeltaSecVtxCoord{"thnConfigAxisDeltaSecVtxCoord", {1000, -10, 10}, ""}; HistogramRegistry registry{"registry", {}}; @@ -143,9 +142,6 @@ struct Alice3HfTask3Prong { addHistogramsRec("hCPAxy", "cosine of pointing angle xy", "entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}); addHistogramsRec("hDca2", "prong Chi2PCA to sec. vertex (cm)", "entries", {HistType::kTH1F, {{400, 0., 20.}}}); addHistogramsRec("hEta", "#it{#eta}", "entries", {HistType::kTH1F, {{100, -2., 2.}}}); - addHistogramsRec("hDeltaXRecoGenSecVtx", "p_{T} (GeV/#it{c})", "#Delta X Sec. Vtx. (cm)", {HistType::kTH2F, {thnConfigAxisPt, thnConfigAxisDeltaSecVtxCoord}}); - addHistogramsRec("hDeltaYRecoGenSecVtx", "p_{T} (GeV/#it{c})", "#Delta Y Sec. Vtx. (cm)", {HistType::kTH2F, {thnConfigAxisPt, thnConfigAxisDeltaSecVtxCoord}}); - addHistogramsRec("hDeltaZRecoGenSecVtx", "p_{T} (GeV/#it{c})", "#Delta Z Sec. Vtx. (cm)", {HistType::kTH2F, {thnConfigAxisPt, thnConfigAxisDeltaSecVtxCoord}}); addHistogramsRec("hd0VsPtProng0", "prong 0 DCAxy to prim. vertex (cm)", "#it{p}_{T} (GeV/#it{c})", {HistType::kTH2F, {{600, -0.4, 0.4}, {vbins}}}); addHistogramsRec("hd0VsPtProng1", "prong 1 DCAxy to prim. vertex (cm)", "#it{p}_{T} (GeV/#it{c})", {HistType::kTH2F, {{600, -0.4, 0.4}, {vbins}}}); addHistogramsRec("hd0VsPtProng2", "prong 2 DCAxy to prim. vertex (cm)", "#it{p}_{T} (GeV/#it{c})", {HistType::kTH2F, {{600, -0.4, 0.4}, {vbins}}}); @@ -261,9 +257,6 @@ struct Alice3HfTask3Prong { registry.fill(histoPrefix + HIST("hImpParErrProng1VsPt") + histoSuffix, candidate.errorImpactParameterY1(), candidate.pt()); registry.fill(histoPrefix + HIST("hImpParErrProng2VsPt") + histoSuffix, candidate.errorImpactParameterY2(), candidate.pt()); registry.fill(histoPrefix + HIST("hDecLenErrVsPt") + histoSuffix, candidate.errorDecayLength(), candidate.pt()); - registry.fill(histoPrefix + HIST("hDeltaXRecoGenSecVtx") + histoSuffix, candidate.pt(), candidate.deltaXSecVtx()); - registry.fill(histoPrefix + HIST("hDeltaYRecoGenSecVtx") + histoSuffix, candidate.pt(), candidate.deltaYSecVtx()); - registry.fill(histoPrefix + HIST("hDeltaZRecoGenSecVtx") + histoSuffix, candidate.pt(), candidate.deltaZSecVtx()); } /// Fill MC histograms at reconstruction level @@ -271,8 +264,8 @@ struct Alice3HfTask3Prong { /// \tparam SaveMl indicates whether ML scores are saved in the THnSparse /// \tparam CandsRec is the type of the reconstructed candidates collection /// \param candidates is the collection of reconstructed candidates - template - void fillHistosMcRec(CandsRec const& candidates, GenParticles const& genParticles) + template + void fillHistosMcRec(CandsRec const& candidates, AllParticles const& allParticles) { registry.fill(HIST("MC/rec/hCandidateCounter"), 0.); for (const auto& candidate : candidates) { @@ -284,13 +277,13 @@ struct Alice3HfTask3Prong { registry.fill(HIST("MC/rec/hCandidateCounter"), 2.); if (candidate.particleMcRec() >= 0) { registry.fill(HIST("MC/rec/hCandidateCounter"), 3.); - auto mcParticle = genParticles.iteratorAt(candidate.particleMcRec()); + auto mcParticle = allParticles.iteratorAt(candidate.particleMcRec()); if (mcParticle.has_daughters()) { auto daughters = mcParticle.daughtersIds(); LOG(debug) << "Reco candidate matched to MC particle with PDG " << mcParticle.pdgCode() << " daughters: " << daughters.size(); int prongIdx = 0; for (int dauId = daughters[0]; dauId <= daughters[1]; dauId++) { - auto dau = genParticles.iteratorAt(dauId); + auto dau = allParticles.iteratorAt(dauId); LOG(debug) << " dauId: " << dauId << " PDG: " << dau.pdgCode() << " with pT: " << dau.pt(); switch (prongIdx) { case 0: @@ -383,9 +376,9 @@ struct Alice3HfTask3Prong { /// \tparam ParticleType is the type of the generated particle /// \tparam TableType is the type of the full table (non-Partition) /// \param particle is a generated particle - /// \param genParticles is the full table of particles for iteratorAt access + /// \param allParticles is the full table of particles for iteratorAt access template - void fillHistogramsGen(ParticleType const& particle, TableType const& genParticles) + void fillHistogramsGen(ParticleType const& particle, TableType const& allParticles) { LOG(debug) << "Filling generated histograms for signal type " << SignalType; static constexpr auto histoPrefix = HIST("MC/gen/") + HIST(SignalFolders[SignalType]) + HIST("/"); @@ -407,7 +400,7 @@ struct Alice3HfTask3Prong { float pz = 0.f; float e = 0.f; for (int iDau = firstDauIdx; iDau <= lastDauIdx && iDau > 0; iDau++) { - const auto& dau = genParticles.iteratorAt(iDau); + const auto& dau = allParticles.iteratorAt(iDau); e += dau.e(); px += dau.px(); py += dau.py(); @@ -422,11 +415,11 @@ struct Alice3HfTask3Prong { /// Fill MC histograms at generated level /// \tparam CharmHad is the charm hadron species /// \tparam CandsGen is the type of the generated candidates collection - /// \tparam GenParticles is the type of the full particle table + /// \tparam AllParticles is the type of the full particle table /// \param mcParticles is the collection of generated particles (can be a Partition) - /// \param genParticles is the full table of particles - template - void fillHistosMcGen(CandsGen const& mcParticles, GenParticles const& genParticles) + /// \param allParticles is the full table of particles + template + void fillHistosMcGen(CandsGen const& mcParticles, AllParticles const& allParticles) { // MC gen. for (const auto& particle : mcParticles) { @@ -438,14 +431,14 @@ struct Alice3HfTask3Prong { const auto ptGen = particle.pt(); const auto originType = particle.originMcGen(); - fillHistogramsGen(particle, genParticles); + fillHistogramsGen(particle, allParticles); float ptGenB = -1.f; if (originType == RecoDecay::OriginType::Prompt) { - fillHistogramsGen(particle, genParticles); + fillHistogramsGen(particle, allParticles); } else if (particle.originMcGen() == RecoDecay::OriginType::NonPrompt) { ptGenB = particle.bHadMotherPtGen(); - fillHistogramsGen(particle, genParticles); + fillHistogramsGen(particle, allParticles); } if (fillThn) { diff --git a/PWGJE/Tasks/CMakeLists.txt b/PWGJE/Tasks/CMakeLists.txt index 63e8fe99f69..a1b103a528c 100644 --- a/PWGJE/Tasks/CMakeLists.txt +++ b/PWGJE/Tasks/CMakeLists.txt @@ -188,10 +188,10 @@ if(FastJet_FOUND) SOURCES v0QA.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore COMPONENT_NAME Analysis) -# o2physics_add_dpl_workflow(jet-finder-charged-qa -# SOURCES jetFinderQA.cxx -# PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore -# COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(jet-finder-charged-qa + SOURCES jetFinderQA.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore + COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(jet-outlier-qa SOURCES jetOutlierQA.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore @@ -366,11 +366,11 @@ if(FastJet_FOUND) SOURCES dijetFinderQA.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore COMPONENT_NAME Analysis) - # o2physics_add_dpl_workflow(bjet-tagging-gnn - # SOURCES bjetTaggingGnn.cxx - # PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore O2Physics::EventFilteringUtils - # O2Physics::AnalysisCCDB - # COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(bjet-tagging-gnn + SOURCES bjetTaggingGnn.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore O2Physics::EventFilteringUtils + O2Physics::AnalysisCCDB + COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(jet-shape SOURCES jetShape.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore From 112dc168936feabb35f27ac37e48d234d75d6d83 Mon Sep 17 00:00:00 2001 From: marcellocosti Date: Tue, 12 May 2026 14:00:58 +0200 Subject: [PATCH 4/5] Updates --- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index 0cede1f5fa0..5aef7ca3f26 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -1741,9 +1741,8 @@ struct OnTheFlyTracker { std::vector genCascades; std::vector genV0s; bcData.clear(); - - o2::dataformats::DCA dcaInfo; - // o2::dataformats::VertexBase vtx; + recoPrimaries.clear(); + ghostPrimaries.clear(); // *+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+*+~+* // Study collision and perform vertexing @@ -1856,7 +1855,7 @@ struct OnTheFlyTracker { // Compute primary vertex with primary tracks o2::vertexing::PVertex primaryVertex; - if (recoPrimaries.size() <= 2) { + if (enablePrimaryVertexing && recoPrimaries.size() <= 2) { LOG(info) << "Not enough primary tracks (" << recoPrimaries.size() << ") to compute vertex, skipping vertexing and collision population."; return; } @@ -1921,6 +1920,7 @@ struct OnTheFlyTracker { processWithLUTs(mcCollision, mcParticles, static_cast(icfg)); } } + void processConfigurationDev(aod::McCollision const& mcCollision, aod::McPartWithDaus const& mcParticles, const int icfg) { // const int lastTrackIndex = tableStoredTracksCov.lastIndex() + 1; // bookkeep the last added track @@ -2032,7 +2032,7 @@ struct OnTheFlyTracker { } o2::vertexing::PVertex primaryVertex; - if (recoPrimaries.size() <= 2) { + if (enablePrimaryVertexing && recoPrimaries.size() <= 2) { LOG(info) << "Not enough primary tracks (" << recoPrimaries.size() << ") to compute vertex, skipping vertexing and collision population."; return; } From 7e50cf4b5c461898020621aa01bb11af0fc7b25b Mon Sep 17 00:00:00 2001 From: marcellocosti Date: Tue, 12 May 2026 22:26:12 +0200 Subject: [PATCH 5/5] Fix TPC cluster condition --- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 149 ++++++++++--------- 1 file changed, 76 insertions(+), 73 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index 5aef7ca3f26..62444181c49 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -495,7 +495,6 @@ struct OnTheFlyTracker { getHist(TH1, histPath + "hVtxTrials")->GetXaxis()->SetBinLabel(2, "Succeeded"); } - LOG(info) << "Enable secondary smearing for configuration " << icfg << ": " << enableSecondarySmearing.value; if (enableSecondarySmearing) { fastTracker.emplace_back(std::make_unique()); fastTracker[icfg]->SetMagneticField(mMagneticField); @@ -506,36 +505,43 @@ struct OnTheFlyTracker { fastTracker[icfg]->Print(); // print fastTracker settings if (cascadeDecaySettings.doXiQA) { - LOG(info) << "Enabling Xi QA histograms for configuration " << icfg; insertHist(histPath + "hXiBuilding", "hXiBuilding", {kTH1F, {{10, -0.5f, 9.5f}}}); - - insertHist(histPath + "hGenXi", "hGenXi", {kTH2F, {axes.axisDecayRadius, axes.axisMomentum}}); - insertHist(histPath + "hRecoXi", "hRecoXi", {kTH2F, {axes.axisDecayRadius, axes.axisMomentum}}); - - insertHist(histPath + "hGenPiFromXi", "hGenPiFromXi", {kTH2F, {axes.axisDecayRadius, axes.axisMomentum}}); - insertHist(histPath + "hGenPiFromLa", "hGenPiFromLa", {kTH2F, {axes.axisDecayRadius, axes.axisMomentum}}); - insertHist(histPath + "hGenPrFromLa", "hGenPrFromLa", {kTH2F, {axes.axisDecayRadius, axes.axisMomentum}}); - insertHist(histPath + "hRecoPiFromXi", "hRecoPiFromXi", {kTH2F, {axes.axisDecayRadius, axes.axisMomentum}}); - insertHist(histPath + "hRecoPiFromLa", "hRecoPiFromLa", {kTH2F, {axes.axisDecayRadius, axes.axisMomentum}}); - insertHist(histPath + "hRecoPrFromLa", "hRecoPrFromLa", {kTH2F, {axes.axisDecayRadius, axes.axisMomentum}}); + getHist(TH1, histPath + "hXiBuilding")->GetXaxis()->SetBinLabel(1, "Generated"); + getHist(TH1, histPath + "hXiBuilding")->GetXaxis()->SetBinLabel(2, "Secondary smearing prong 0"); + getHist(TH1, histPath + "hXiBuilding")->GetXaxis()->SetBinLabel(3, "Secondary smearing prong 1"); + getHist(TH1, histPath + "hXiBuilding")->GetXaxis()->SetBinLabel(4, "Secondary smearing prong 2"); + getHist(TH1, histPath + "hXiBuilding")->GetXaxis()->SetBinLabel(5, "Not Nan"); + getHist(TH1, histPath + "hXiBuilding")->GetXaxis()->SetBinLabel(6, "Start Reco"); + getHist(TH1, histPath + "hXiBuilding")->GetXaxis()->SetBinLabel(7, "V0 fitter ok"); + getHist(TH1, histPath + "hXiBuilding")->GetXaxis()->SetBinLabel(8, "Kink fitter ok"); + + insertHist(histPath + "hGenXi", "hGenXi;Decay Radius;#it{p}_{T}", {kTH2F, {axes.axisDecayRadius, axes.axisMomentum}}); + insertHist(histPath + "hRecoXi", "hRecoXi;Decay Radius;#it{p}_{T}", {kTH2F, {axes.axisDecayRadius, axes.axisMomentum}}); + + insertHist(histPath + "hGenPiFromXi", "hGenPiFromXi;Decay Radius;#it{p}_{T}", {kTH2F, {axes.axisDecayRadius, axes.axisMomentum}}); + insertHist(histPath + "hGenPiFromLa", "hGenPiFromLa;Decay Radius;#it{p}_{T}", {kTH2F, {axes.axisDecayRadius, axes.axisMomentum}}); + insertHist(histPath + "hGenPrFromLa", "hGenPrFromLa;Decay Radius;#it{p}_{T}", {kTH2F, {axes.axisDecayRadius, axes.axisMomentum}}); + insertHist(histPath + "hRecoPiFromXi", "hRecoPiFromXi;Decay Radius;#it{p}_{T}", {kTH2F, {axes.axisDecayRadius, axes.axisMomentum}}); + insertHist(histPath + "hRecoPiFromLa", "hRecoPiFromLa;Decay Radius;#it{p}_{T}", {kTH2F, {axes.axisDecayRadius, axes.axisMomentum}}); + insertHist(histPath + "hRecoPrFromLa", "hRecoPrFromLa;Decay Radius;#it{p}_{T}", {kTH2F, {axes.axisDecayRadius, axes.axisMomentum}}); // basic mass histograms to see if we're in business - insertHist(histPath + "hMassLambda", "hMassLambda", {kTH1F, {{axes.axisLambdaMass}}}); - insertHist(histPath + "hMassXi", "hMassXi", {kTH1F, {{axes.axisXiMass}}}); - insertHist(histPath + "h2dMassXi", "h2dMassXi", {kTH2F, {{axes.axisXiMass}, {axes.axisMomentum}}}); + insertHist(histPath + "hMassLambda", "hMassLambda;Mass;Counts", {kTH1F, {{axes.axisLambdaMass}}}); + insertHist(histPath + "hMassXi", "hMassXi;Mass;Counts", {kTH1F, {{axes.axisXiMass}}}); + insertHist(histPath + "h2dMassXi", "h2dMassXi;Mass;Counts", {kTH2F, {{axes.axisXiMass}, {axes.axisMomentum}}}); // OTF strangeness tracking QA - insertHist(histPath + "hFoundVsFindable", "hFoundVsFindable", {kTH2F, {{10, -0.5f, 9.5f}, {10, -0.5f, 9.5f}}}); + insertHist(histPath + "hFoundVsFindable", "hFoundVsFindable;Found;Findable", {kTH2F, {{10, -0.5f, 9.5f}, {10, -0.5f, 9.5f}}}); - insertHist(histPath + "h2dDCAxyCascade", "h2dDCAxyCascade", {kTH2F, {axes.axisMomentum, axes.axisDCA}}); - insertHist(histPath + "h2dDCAxyCascadeBachelor", "h2dDCAxyCascadeBachelor", {kTH2F, {axes.axisMomentum, axes.axisDCA}}); - insertHist(histPath + "h2dDCAxyCascadeNegative", "h2dDCAxyCascadeNegative", {kTH2F, {axes.axisMomentum, axes.axisDCA}}); - insertHist(histPath + "h2dDCAxyCascadePositive", "h2dDCAxyCascadePositive", {kTH2F, {axes.axisMomentum, axes.axisDCA}}); + insertHist(histPath + "h2dDCAxyCascade", "h2dDCAxyCascade;#it{p}_{T};DCA_{xy}", {kTH2F, {axes.axisMomentum, axes.axisDCA}}); + insertHist(histPath + "h2dDCAxyCascadeBachelor", "h2dDCAxyCascadeBachelor;#it{p}_{T};DCA_{xy}", {kTH2F, {axes.axisMomentum, axes.axisDCA}}); + insertHist(histPath + "h2dDCAxyCascadeNegative", "h2dDCAxyCascadeNegative;#it{p}_{T};DCA_{xy}", {kTH2F, {axes.axisMomentum, axes.axisDCA}}); + insertHist(histPath + "h2dDCAxyCascadePositive", "h2dDCAxyCascadePositive;#it{p}_{T};DCA_{xy}", {kTH2F, {axes.axisMomentum, axes.axisDCA}}); - insertHist(histPath + "h2dDCAzCascade", "h2dDCAzCascade", {kTH2F, {axes.axisMomentum, axes.axisDCA}}); - insertHist(histPath + "h2dDCAzCascadeBachelor", "h2dDCAzCascadeBachelor", {kTH2F, {axes.axisMomentum, axes.axisDCA}}); - insertHist(histPath + "h2dDCAzCascadeNegative", "h2dDCAzCascadeNegative", {kTH2F, {axes.axisMomentum, axes.axisDCA}}); - insertHist(histPath + "h2dDCAzCascadePositive", "h2dDCAzCascadePositive", {kTH2F, {axes.axisMomentum, axes.axisDCA}}); + insertHist(histPath + "h2dDCAzCascade", "h2dDCAzCascade;#it{p}_{T};DCA_{z}", {kTH2F, {axes.axisMomentum, axes.axisDCA}}); + insertHist(histPath + "h2dDCAzCascadeBachelor", "h2dDCAzCascadeBachelor;#it{p}_{T};DCA_{z}", {kTH2F, {axes.axisMomentum, axes.axisDCA}}); + insertHist(histPath + "h2dDCAzCascadeNegative", "h2dDCAzCascadeNegative;#it{p}_{T};DCA_{z}", {kTH2F, {axes.axisMomentum, axes.axisDCA}}); + insertHist(histPath + "h2dDCAzCascadePositive", "h2dDCAzCascadePositive;#it{p}_{T};DCA_{z}", {kTH2F, {axes.axisMomentum, axes.axisDCA}}); insertHist(histPath + "h2dDeltaPtVsPt", "h2dDeltaPtVsPt;Gen p_{T};#Delta p_{T}", {kTH2F, {axes.axisMomentum, axes.axisDeltaPt}}); insertHist(histPath + "h2dDeltaEtaVsPt", "h2dDeltaEtaVsPt;Gen p_{T};#Delta #eta", {kTH2F, {axes.axisMomentum, axes.axisDeltaEta}}); @@ -550,7 +556,6 @@ struct OnTheFlyTracker { getHist(TH1, histPath + "hFastTrackerQA")->GetXaxis()->SetBinLabel(6, "multiple scattering"); getHist(TH1, histPath + "hFastTrackerQA")->GetXaxis()->SetBinLabel(7, "energy loss"); getHist(TH1, histPath + "hFastTrackerQA")->GetXaxis()->SetBinLabel(8, "efficiency"); - getHist(TH1, histPath + "hFastTrackerQA")->GetXaxis()->SetBinLabel(9, "no layers hit"); } } @@ -893,9 +898,6 @@ struct OnTheFlyTracker { int dNdEta, float eventCollisionTimeNS) { - LOG(info) << ""; - LOG(info) << "-------------------------------"; - LOG(info) << "Studying cascade with PDG " << mcParticle.pdgCode() << " and pT " << mcParticle.pt(); o2::track::TrackParCov trackParCov; o2::upgrade::convertMCParticleToO2Track(mcParticle, trackParCov, pdgDB); const std::string histPath = "Configuration_" + std::to_string(icfg) + "/"; @@ -912,7 +914,6 @@ struct OnTheFlyTracker { o2::track::TrackParCov xiTrackParCov; o2::upgrade::convertMCParticleToO2Track(mcParticle, xiTrackParCov, pdgDB); - LOG(info) << "Decaying cascade..."; decayCascade(mcParticle, xiTrackParCov, cascadeDecayProducts, xiDecayVertex, laDecayVertex); if (cascadeDecayProducts.size() != 3) { LOG(fatal) << "Xi decay did not produce 3 daughters as expected!"; @@ -941,6 +942,11 @@ struct OnTheFlyTracker { o2::upgrade::convertMCParticleToO2Track(mcParticle, perfectCascadeTrack, pdgDB); perfectCascadeTrack.setPID(pdgCodeToPID(PDG_t::kXiMinus)); // FIXME: not OK for omegas + thisCascade.bachelorId = trackTableOffset + 1; + thisCascade.negativeId = trackTableOffset + 2; // Lambda daughters + thisCascade.positiveId = trackTableOffset + 3; // Lambda daughters + thisCascade.cascadeTrackId = trackTableOffset + 4; + float trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; tracksCascadeProngs.push_back(TrackAlice3{xiDaughterTrackParCovsPerfect[0], mcParticle.globalIndex(), trackTime, timeResolutionUs, true, true, 1, -1, -1, TrackType::kGenCascDaug}); trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; @@ -949,7 +955,6 @@ struct OnTheFlyTracker { tracksCascadeProngs.push_back(TrackAlice3{xiDaughterTrackParCovsPerfect[2], mcParticle.globalIndex(), trackTime, timeResolutionUs, true, true, 3, -1, -1, TrackType::kGenCascDaug}); trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; tracksCascadeProngs.push_back(TrackAlice3{perfectCascadeTrack, mcParticle.globalIndex(), trackTime, timeResolutionUs, true, true, 3, -1, -1, TrackType::kGenCascDaug}); - LOG(info) << "tracksCascadeProngs.size() after adding gen tracks: " << tracksCascadeProngs.size(); for (int i = 0; i < kCascProngs; i++) { isReco[i] = false; @@ -965,17 +970,15 @@ struct OnTheFlyTracker { getHist(TH1, histPath + "hFastTrackerQA")->Fill(o2::math_utils::abs(nHitsCascadeProngs[i])); } - if (nSiliconHitsCascadeProngs[i] < fastTrackerSettings.minSiliconHits) { - LOG(info) << "Prong " << i << " failed fast tracker silicon hit requirement with " << nSiliconHitsCascadeProngs[i] << " hits, required: " << fastTrackerSettings.minSiliconHits; - continue; - } - if (nSiliconHitsCascadeProngs[i] < fastTrackerSettings.minSiliconHitsIfTPCUsed) { - LOG(info) << "Prong " << i << " failed fast tracker silicon hit requirement with " << nSiliconHitsCascadeProngs[i] << " hits, required: " << fastTrackerSettings.minSiliconHitsIfTPCUsed; - continue; + if (nSiliconHitsCascadeProngs[i] >= fastTrackerSettings.minSiliconHits || + (nSiliconHitsCascadeProngs[i] >= fastTrackerSettings.minSiliconHitsIfTPCUsed && + nTPCHitsCascadeProngs[i] >= fastTrackerSettings.minTPCClusters)) { + isReco[i] = true; + } else { + continue; // extra sure } - if (nTPCHitsCascadeProngs[i] < fastTrackerSettings.minTPCClusters) { - LOG(info) << "Prong " << i << " failed fast tracker TPC hit requirement with " << nTPCHitsCascadeProngs[i] << " hits, required: " << fastTrackerSettings.minTPCClusters; - continue; + if (cascadeDecaySettings.doXiQA) { + getHist(TH1, histPath + "hXiBuilding")->Fill(static_cast(i+1)); } isReco[i] = true; @@ -985,11 +988,15 @@ struct OnTheFlyTracker { } else { isReco[i] = true; xiDaughterTrackParCovsTracked[i] = xiDaughterTrackParCovsPerfect[i]; + if (cascadeDecaySettings.doXiQA) { + getHist(TH1, histPath + "hXiBuilding")->Fill(static_cast(i+1)); + } } if (TMath::IsNaN(xiDaughterTrackParCovsTracked[i].getZ())) { continue; } else { + getHist(TH1, histPath + "hXiBuilding")->Fill(4.0f); histos.fill(HIST("hNaNBookkeeping"), i + 1, 1.0f); } trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; @@ -1005,10 +1012,11 @@ struct OnTheFlyTracker { bool reconstructedCascade = false; if (isReco[0] && isReco[1] && isReco[2]) { - LOG(info) << "All cascade daughters reconstructed, attempting cascade reconstruction with DCA finder..."; reconstructedCascade = true; } + bool fillCascadeTable{false}; + // +-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+ // combine particles into actual Xi candidate // cascade building starts here @@ -1016,15 +1024,6 @@ struct OnTheFlyTracker { if (cascadeDecaySettings.doXiQA) { getHist(TH1, histPath + "hXiBuilding")->Fill(3.0f); } - // 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.bachelorId = trackTableOffset + 1; - thisCascade.positiveId = trackTableOffset + 2; // Lambda daughters - thisCascade.negativeId = trackTableOffset + 3; // Lambda daughters - thisCascade.cascadeTrackId = trackTableOffset + 4; // use DCA fitters int nCand = 0; @@ -1046,7 +1045,6 @@ struct OnTheFlyTracker { // V0 found successfully if (dcaFitterOK_V0) { - LOG(info) << "V0 candidate found successfully with DCA fitter, now trying to combine with bachelor track to form cascade..."; if (cascadeDecaySettings.doXiQA) { getHist(TH1, histPath + "hXiBuilding")->Fill(4.0f); } @@ -1115,7 +1113,7 @@ struct OnTheFlyTracker { // Cascade found successfully if (dcaFitterOK_Cascade) { if (cascadeDecaySettings.doXiQA) { - getHist(TH1, histPath + "hXiBuilding")->Fill(5.0f); + getHist(TH1, histPath + "hXiBuilding")->Fill(6.0f); } o2::track::TrackParCov bachelorTrackAtPCA = fitter.getTrack(1); @@ -1214,7 +1212,8 @@ struct OnTheFlyTracker { return; // We didn't find enough hits for strangeness tracking } } - tracksCascadeProngs[kCascProngs + 1] = TrackAlice3{cascadeTrack, mcParticle.globalIndex(), trackTime, timeResolutionUs, false, false, 1, thisCascade.foundClusters, TrackType::kGenCascDaug}; + tracksCascadeProngs[kCascProngs] = TrackAlice3{cascadeTrack, mcParticle.globalIndex(), trackTime, timeResolutionUs, false, false, 1, thisCascade.foundClusters, TrackType::kGenCascDaug}; + fillCascadeTable = true; } } } // end cascade building @@ -1263,7 +1262,7 @@ struct OnTheFlyTracker { histos.fill(HIST("hFitterStatusCode"), fitterStatusCode); if (kinkFitterOK) { if (cascadeDecaySettings.doXiQA) { - getHist(TH1, histPath + "hXiBuilding")->Fill(6.0f); + getHist(TH1, histPath + "hXiBuilding")->Fill(7.0f); } o2::track::TrackParCov newCascadeTrack = fitter.getTrack(0); // (cascade) @@ -1281,12 +1280,12 @@ struct OnTheFlyTracker { thisCascade.eta = newCascadeTrack.getEta(); thisCascade.mXi = RecoDecay::m(std::array{std::array{pBach[0], pBach[1], pBach[2]}, std::array{pV0[0], pV0[1], pV0[2]}}, std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassLambda}); - newCascadeTrack.setPID(pdgCodeToPID(PDG_t::kXiMinus)); // FIXME: not OK for omegas float trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; if (reconstructedCascade) { tracksCascadeProngs[kCascProngs + 1] = TrackAlice3{newCascadeTrack, mcParticle.globalIndex(), trackTime, timeResolutionUs, false, false, 1, thisCascade.foundClusters, TrackType::kRecoCascDaug}; } + fillCascadeTable = true; } // end fitter OK } // end cascade found } // end cascade kink building @@ -1340,6 +1339,11 @@ struct OnTheFlyTracker { } } + if (!fillCascadeTable) { + tracksCascadeProngs.clear(); // clear the tracks added for this cascade since we won't be filling the table + return; + } + // populate Cascades tableUpgradeCascades(tableCollisions.lastIndex(), thisCascade.cascadeTrackId, @@ -1432,6 +1436,7 @@ struct OnTheFlyTracker { trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; tracksV0Daugs.push_back(TrackAlice3{v0DaughterTrackParCovsPerfect[1], mcParticle.globalIndex(), 0.f, timeResolutionUs, true, true, 1, TrackType::kGhostV0Daug}); + bool fillV0Table{false}; for (int i = 0; i < kv0Prongs; i++) { isV0Reco[i] = false; nV0Hits[i] = 0; @@ -1442,13 +1447,13 @@ struct OnTheFlyTracker { nV0SiliconHits[i] = fastTracker[icfg]->GetNSiliconPoints(); nV0TPCHits[i] = fastTracker[icfg]->GetNGasPoints(); - if (nV0SiliconHits[i] < fastTrackerSettings.minSiliconHits) - continue; - if (nV0SiliconHits[i] < fastTrackerSettings.minSiliconHitsIfTPCUsed) - continue; - if (nV0TPCHits[i] < fastTrackerSettings.minTPCClusters) - continue; - isV0Reco[i] = true; + if (nV0SiliconHits[i] >= fastTrackerSettings.minSiliconHits || + (nV0SiliconHits[i] >= fastTrackerSettings.minSiliconHitsIfTPCUsed && + nV0TPCHits[i] >= fastTrackerSettings.minTPCClusters)) { + isV0Reco[i] = true; + } else { + continue; // extra sure + } if (v0DecaySettings.doV0QA) { // QA if (nV0Hits[i] < 0) { @@ -1591,9 +1596,17 @@ struct OnTheFlyTracker { fillHist(TH2, Form("V0Building_Configuration_%i/AntiLambda/hMass", icfg), thisV0.mAntiLambda, thisV0.pt); } } + + LOG(info) << "Will fill V0 table for MC particle with PDG " << mcParticle.pdgCode() << " and global index " << mcParticle.globalIndex(); + fillV0Table = true; } } + if (!fillV0Table) { + tracksV0Daugs.clear(); // clear the tracks added for this V0 since we won't be filling the table + return; // don't fill the table if we didn't find a V0 candidate + } + // populate V0s tableUpgradeV0s(tableCollisions.lastIndex(), // now we know the collision index -> populate table thisV0.mcParticleId, @@ -1735,7 +1748,6 @@ struct OnTheFlyTracker { void processWithLUTs(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles, const int icfg) { - LOG(info) << "Processing collision " << mcCollision.globalIndex() << " with " << mcParticles.size() << " particles"; const std::string histPath = "Configuration_" + std::to_string(icfg) + "/"; std::vector genCascades; @@ -1777,7 +1789,6 @@ struct OnTheFlyTracker { (cascadeDecaySettings.decayXi && isCascade) || (v0DecaySettings.decayV0 && isV0); if (!pdgsToBeHandled) { - // LOG(info) << "Particle with PDG " << mcParticle.pdgCode() << " is not in the list of particles to be handled, skipping."; continue; } @@ -1819,7 +1830,6 @@ struct OnTheFlyTracker { continue; } - // LOG(info) << "Particle with PDG " << mcParticle.pdgCode() << (reconstructed ? " was reconstructed." : " was not reconstructed."); getHist(TH1, histPath + "hPtReconstructed")->Fill(trackParCov.getPt()); if (std::abs(mcParticle.pdgCode()) == kElectron) getHist(TH1, histPath + "hPtReconstructedEl")->Fill(trackParCov.getPt()); @@ -1842,7 +1852,6 @@ struct OnTheFlyTracker { } histos.fill(HIST("hNaNBookkeeping"), 0.0f, 1.0f); // ok! - // LOG(info) << "Particle with PDG " << mcParticle.pdgCode() << " has " << nTrkHits << " hits."; // Time associated to the mcParticle: collision time + smearing float trackTime = (eventCollisionTimeNS + gRandom->Gaus(0., timeResolutionNs)) * nsToMus; TrackType trackType = reconstructed ? TrackType::kRecoPrimary : TrackType::kGhostPrimary; @@ -1860,7 +1869,6 @@ struct OnTheFlyTracker { return; } computeVertex(mcCollision, recoPrimaries, primaryVertex, icfg); - // LOG(info) << "Primary vertex position: (" << primaryVertex.getX() << ", " << primaryVertex.getY() << ", " << primaryVertex.getZ() << ")"; getHist(TH1, histPath + "hPVz")->Fill(primaryVertex.getZ()); // populate collisions tableCollisions(-1, // BC is irrelevant in synthetic MC tests for now, could be adjusted in future @@ -1879,16 +1887,15 @@ struct OnTheFlyTracker { for (size_t iCasc = 0; iCasc < genCascades.size(); iCasc++) { auto genCasc = mcParticles.rawIteratorAt(genCascades[iCasc] - mcParticles.offset()); studyCascade(trackTableOffset, genCasc, tracksCascadeProngs, primaryVertex, icfg, dNdEta, eventCollisionTimeNS); - LOG(info) << "[CASCADE] Populating number of cascade tracks: " << tracksCascadeProngs.size(); fillTracksInfo(tracksCascadeProngs, primaryVertex, icfg); trackTableOffset += tracksCascadeProngs.size(); // each cascade takes 4 tracks in the table (cascade + 3 daughters) tracksCascadeProngs.clear(); } + std::vector tracksV0Daugs; for (size_t iV0 = 0; iV0 < genV0s.size(); iV0++) { auto genV0 = mcParticles.rawIteratorAt(genV0s[iV0] - mcParticles.offset()); studyV0(trackTableOffset, genV0, tracksV0Daugs, icfg, dNdEta, eventCollisionTimeNS); - LOG(info) << "Populating " << tracksV0Daugs.size() << " tracks associated to V0 candidate"; fillTracksInfo(tracksV0Daugs, primaryVertex, icfg); trackTableOffset += tracksV0Daugs.size(); // each V0 takes 2 tracks in the table (2 daughters) tracksV0Daugs.clear(); @@ -1923,7 +1930,6 @@ struct OnTheFlyTracker { void processConfigurationDev(aod::McCollision const& mcCollision, aod::McPartWithDaus const& mcParticles, const int icfg) { - // const int lastTrackIndex = tableStoredTracksCov.lastIndex() + 1; // bookkeep the last added track const std::string histPath = "Configuration_" + std::to_string(icfg) + "/"; tracksAlice3.clear(); ghostTracksAlice3.clear(); @@ -1931,9 +1937,6 @@ struct OnTheFlyTracker { cascadesAlice3.clear(); v0sAlice3.clear(); - o2::dataformats::DCA dcaInfo; - o2::dataformats::VertexBase vtx; - // generate collision time auto ir = irSampler.generateCollisionTime(); const float eventCollisionTimeNS = ir.timeInBCNS;