diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index 21e6b6a63f8..fcf19d4f16e 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -146,6 +146,7 @@ struct OnTheFlyTracker { struct : ConfigurableGroup { std::string prefix = "fastTrackerSettings"; // JSON group name Configurable minSiliconHits{"minSiliconHits", 6, "minimum number of silicon hits to accept track"}; + Configurable minSiliconHitsForKinkReco{"minSiliconHitsForKinkReco", 4, "minimum number of silicon hits to accept track"}; Configurable minSiliconHitsIfTPCUsed{"minSiliconHitsIfTPCUsed", 2, "minimum number of silicon hits to accept track in case TPC info is present"}; Configurable minTPCClusters{"minTPCClusters", 70, "minimum number of TPC hits necessary to consider minSiliconHitsIfTPCUsed"}; Configurable applyZacceptance{"applyZacceptance", false, "apply z limits to detector layers or not"}; @@ -168,6 +169,7 @@ struct OnTheFlyTracker { Configurable findXi{"findXi", false, "if decayXi on, find Xi and fill Tracks table also with Xi"}; Configurable trackXi{"trackXi", false, "if findXi on, attempt to track Xi"}; Configurable doXiQA{"doXiQA", false, "QA plots for when treating Xi"}; + Configurable doKinkReco{"doKinkReco", 0, "Flag for kink reco setting: 0 - disabled, 1 - complementary, 2 - only"}; } cascadeDecaySettings; struct : ConfigurableGroup { @@ -866,7 +868,10 @@ struct OnTheFlyTracker { std::vector nHits(kCascProngs); // total std::vector nSiliconHits(kCascProngs); // silicon type std::vector nTPCHits(kCascProngs); // TPC type + + bool tryKinkReco = false; if (cascadeDecaySettings.decayXi && mcParticle.pdgCode() == kXiMinus) { + bool reconstructedCascade = false; if (cascadeDecaySettings.doXiQA) { getHist(TH1, histPath + "hXiBuilding")->Fill(0.0f); } @@ -914,23 +919,18 @@ struct OnTheFlyTracker { } } - if (cascadeDecaySettings.doXiQA && mcParticle.pdgCode() == kXiMinus) { - if (isReco[0] && isReco[1] && isReco[2]) { - getHist(TH1, histPath + "hXiBuilding")->Fill(2.0f); - getHist(TH2, histPath + "hRecoXi")->Fill(xiDecayRadius2D, mcParticle.pt()); - } - if (isReco[0]) - getHist(TH2, histPath + "hRecoPiFromXi")->Fill(xiDecayRadius2D, cascadeDecayProducts[0].Pt()); - if (isReco[1]) - getHist(TH2, histPath + "hRecoPiFromLa")->Fill(laDecayRadius2D, cascadeDecayProducts[1].Pt()); - if (isReco[2]) - getHist(TH2, histPath + "hRecoPrFromLa")->Fill(laDecayRadius2D, cascadeDecayProducts[2].Pt()); + if (!isReco[1] || !isReco[2]) { + tryKinkReco = true; // Lambda outside acceptance, set flag for kink reco to be used if mode 1 + } + + if (isReco[0] && isReco[1] && isReco[2]) { + reconstructedCascade = true; } // +-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+ // combine particles into actual Xi candidate // cascade building starts here - if (cascadeDecaySettings.findXi && mcParticle.pdgCode() == kXiMinus && isReco[0] && isReco[1] && isReco[2]) { + if (cascadeDecaySettings.findXi && isReco[0] && isReco[1] && isReco[2] && cascadeDecaySettings.doKinkReco != 2) { if (cascadeDecaySettings.doXiQA) { getHist(TH1, histPath + "hXiBuilding")->Fill(3.0f); } @@ -1111,24 +1111,111 @@ struct OnTheFlyTracker { 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}); - if (cascadeDecaySettings.doXiQA) { - getHist(TH1, histPath + "hXiBuilding")->Fill(6.0f); - getHist(TH2, histPath + "h2dDeltaPtVsPt")->Fill(trackParCov.getPt(), cascadeTrack.getPt() - trackParCov.getPt()); - getHist(TH2, histPath + "h2dDeltaEtaVsPt")->Fill(trackParCov.getPt(), cascadeTrack.getEta() - trackParCov.getEta()); - - getHist(TH1, histPath + "hMassLambda")->Fill(thisCascade.mLambda); - getHist(TH1, histPath + "hMassXi")->Fill(thisCascade.mXi); - getHist(TH2, histPath + "hFoundVsFindable")->Fill(thisCascade.findableClusters, thisCascade.foundClusters); - } - // add this cascade to vector (will fill cursor later with collision ID) cascadesAlice3.push_back(thisCascade); } } } // end cascade building + + if (isReco[0] && ((cascadeDecaySettings.doKinkReco == 1 && tryKinkReco) || cascadeDecaySettings.doKinkReco == 2)) { // mode 1 or 2 + o2::track::TrackParCov prefectCascadeTrack, trackedCasc; + const o2::track::TrackParCov& trackedBach = xiDaughterTrackParCovsTracked[0]; + o2::upgrade::convertMCParticleToO2Track(mcParticle, prefectCascadeTrack, pdgDB); + + // back track is already smeared + int nCascHits = fastTracker[icfg]->FastTrack(prefectCascadeTrack, trackedCasc, dNdEta); + reconstructedCascade = (fastTrackerSettings.minSiliconHitsForKinkReco < nCascHits) ? false : true; + + if (reconstructedCascade) { + std::array pCasc; + std::array pBach; + std::array pV0; + trackedCasc.getPxPyPzGlo(pCasc); + trackedBach.getPxPyPzGlo(pBach); + for (size_t i = 0; i < pCasc.size(); ++i) { + pV0[i] = pCasc[i] - pBach[i]; + } + + 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; + } + + int nCand = 0; + bool kinkFitterOK = true; + try { + nCand = fitter.process(trackedCasc, trackedBach); + } catch (...) { + kinkFitterOK = false; + } + + if (nCand == 0) { + kinkFitterOK = false; + } + + if (kinkFitterOK) { + if (cascadeDecaySettings.doXiQA) { + getHist(TH1, histPath + "hXiBuilding")->Fill(6.0f); + } + } + + 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.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}); + + tracksAlice3.push_back(TrackAlice3{trackedCasc, 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 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 + "h2dDeltaPtVsPt")->Fill(trackParCov.getPt(), cascadeTrack.getPt() - trackParCov.getPt()); + // getHist(TH2, histPath + "h2dDeltaEtaVsPt")->Fill(trackParCov.getPt(), cascadeTrack.getEta() - trackParCov.getEta()); + 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()); + } + if (isReco[2]) { + getHist(TH2, histPath + "hRecoPrFromLa")->Fill(laDecayRadius2D, cascadeDecayProducts[2].Pt()); + } + } + continue; // Cascade handling done, should not be considered anymore } + // V0 handling std::vector v0DaughterTrackParCovsPerfect(2); std::vector v0DaughterTrackParCovsTracked(2); diff --git a/ALICE3/TableProducer/alice3MulticharmFinder.cxx b/ALICE3/TableProducer/alice3MulticharmFinder.cxx index dd799366d18..9b636b03bb4 100644 --- a/ALICE3/TableProducer/alice3MulticharmFinder.cxx +++ b/ALICE3/TableProducer/alice3MulticharmFinder.cxx @@ -200,6 +200,13 @@ struct Alice3MulticharmFinder { int nTPCHitsPiCC; } thisXiCCcandidate; + struct ProngInfo { + float pt = 1e+10; + float eta = 1e+10; + float dcaXY = 1e+10; + float dcaZ = 1e+10; + }; + template bool buildDecayCandidateTwoBody(TTrackType const& t0, TTrackType const& t1, float mass0, float mass1) { @@ -557,10 +564,6 @@ struct Alice3MulticharmFinder { } uint32_t nCombinationsC = 0; - auto bach = xiCand.bachTrack_as(); // de-reference bach track - auto neg = xiCand.negTrack_as(); // de-reference neg track - auto pos = xiCand.posTrack_as(); // de-reference pos track - if (!BIT_CHECK(xi.decayMap(), kTrueXiFromXiC)) { continue; } @@ -822,13 +825,38 @@ struct Alice3MulticharmFinder { picc.hasSigPi(), picc.nSigmaPionRich(), getPdgCodeForTrack(picc)); + ProngInfo bachelor, positive, negative; + if (xiCand.has_bachTrack()) { + auto bach = xiCand.bachTrack_as(); // de-reference bach track + bachelor.pt = bach.pt(); + bachelor.eta = bach.eta(); + bachelor.dcaXY = bach.dcaXY(); + bachelor.dcaZ = bach.dcaZ(); + } + + if (xiCand.has_negTrack()) { + auto neg = xiCand.negTrack_as(); // de-reference neg track + negative.pt = neg.pt(); + negative.eta = neg.eta(); + negative.dcaXY = neg.dcaXY(); + negative.dcaZ = neg.dcaZ(); + } + + if (xiCand.has_posTrack()) { + auto pos = xiCand.posTrack_as(); // de-reference pos track + positive.pt = pos.pt(); + positive.eta = pos.eta(); + positive.dcaXY = pos.dcaXY(); + positive.dcaZ = pos.dcaZ(); + } + multiCharmExtra( - bach.pt(), bach.eta(), - bach.dcaXY(), bach.dcaZ(), - pos.pt(), pos.eta(), - pos.dcaXY(), pos.dcaZ(), - neg.pt(), neg.eta(), - neg.dcaXY(), neg.dcaZ(), + bachelor.pt, bachelor.eta, + bachelor.dcaXY, bachelor.dcaZ, + positive.pt, positive.eta, + positive.dcaXY, positive.dcaZ, + negative.pt, negative.eta, + negative.dcaXY, negative.dcaZ, pi1c.eta(), pi2c.eta(), picc.eta()); } }