diff --git a/PWGLF/TableProducer/Nuspex/lnnRecoTask.cxx b/PWGLF/TableProducer/Nuspex/lnnRecoTask.cxx index da509a17295..e5ebc7a704a 100644 --- a/PWGLF/TableProducer/Nuspex/lnnRecoTask.cxx +++ b/PWGLF/TableProducer/Nuspex/lnnRecoTask.cxx @@ -76,6 +76,15 @@ std::shared_ptr hIsMatterGen; std::shared_ptr hIsMatterGenTwoBody; std::shared_ptr hDCAxy3H; +float alphaAP(std::array const& momB, std::array const& momC) +{ + std::array momA = {momB[0] + momC[0], momB[1] + momC[1], momB[2] + momC[2]}; + float momTot = std::sqrt(momA[0] * momA[0] + momA[1] * momA[1] + momA[2] * momA[2]); + float lQlPos = (momB[0] * momA[0] + momB[1] * momA[1] + momB[2] * momA[2]) / momTot; + float lQlNeg = (momC[0] * momA[0] + momC[1] * momA[1] + momC[2] * momA[2]) / momTot; + return (lQlPos - lQlNeg) / (lQlPos + lQlNeg); +} + } // namespace struct lnnCandidate { @@ -142,7 +151,6 @@ struct lnnRecoTask { Configurable nSigmaCutMaxTPC{"nSigmaCutMaxTPC", 5, "triton dEdx cut (n sigma)"}; Configurable nTPCClusMin3H{"nTPCClusMin3H", 80, "triton NTPC clusters cut"}; Configurable ptMinTOF{"ptMinTOF", 0.8, "minimum pt for TOF cut"}; - Configurable ptMaxTOF{"ptMaxTOF", 3.5, "minimum pt for TOF cut"}; Configurable TrTOFMass2Cut{"TrTOFMass2Cut", 5.5, "minimum Triton mass square to TOF"}; Configurable mcSignalOnly{"mcSignalOnly", true, "If true, save only signal in MC"}; @@ -339,16 +347,18 @@ struct lnnRecoTask { hdEdxTot->Fill(-negRigidity, negTrack.tpcSignal()); // ITS only tracks do not have TPC information. TPCnSigma: only lower cut to allow for triton reconstruction - bool is3H = posTrack.hasTPC() && nSigmaTPCpos > nSigmaCutMinTPC && nSigmaTPCpos < nSigmaCutMaxTPC; bool isAnti3H = negTrack.hasTPC() && nSigmaTPCneg > nSigmaCutMinTPC && nSigmaTPCneg < nSigmaCutMaxTPC; - if (!is3H && !isAnti3H) + if (!is3H && !isAnti3H) // discard if both tracks are not 3H candidates continue; - // Describing lnn as matter candidate + // if alphaAP is > 0 the candidate is 3H, if < 0 it is anti-3H + std::array momPos = std::array{posTrack.px(), posTrack.py(), posTrack.pz()}; + std::array momNeg = std::array{negTrack.px(), negTrack.py(), negTrack.pz()}; + float alpha = alphaAP(momPos, momNeg); lnnCandidate lnnCand; - lnnCand.isMatter = is3H && isAnti3H ? std::abs(nSigmaTPCpos) < std::abs(nSigmaTPCneg) : is3H; + lnnCand.isMatter = alpha > 0; auto& h3track = lnnCand.isMatter ? posTrack : negTrack; auto& h3Rigidity = lnnCand.isMatter ? posRigidity : negRigidity; @@ -372,29 +382,25 @@ struct lnnRecoTask { lnnCand.flags |= lnnCand.isMatter ? static_cast((posTrack.pidForTracking() & 0xF) << 4) : static_cast((negTrack.pidForTracking() & 0xF) << 4); lnnCand.flags |= lnnCand.isMatter ? static_cast(negTrack.pidForTracking() & 0xF) : static_cast(posTrack.pidForTracking() & 0xF); - auto h3TrackCov = getTrackParCov(posTrack); + auto posTrackCov = getTrackParCov(posTrack); auto negTrackCov = getTrackParCov(negTrack); - int chargeFactor = -1 + 2 * lnnCand.isMatter; - hdEdx3HSel->Fill(chargeFactor * lnnCand.mom3HTPC, h3track.tpcSignal()); - hNsigma3HSel->Fill(chargeFactor * lnnCand.mom3HTPC, lnnCand.nSigma3H); - hDCAxy3H->Fill(h3track.pt(), h3track.dcaXY()); - if (is3H) { - hdEdx3HPosTrack->Fill(lnnCand.mom3HTPC, h3track.tpcSignal()); - } - if (h3track.hasTOF()) { - float beta = h3track.beta(); + + float beta = -1.f; + if (h3track.pt() >= ptMinTOF) { + if (!h3track.hasTOF()) { + continue; + } + beta = h3track.beta(); lnnCand.mass2TrTOF = h3track.mass() * h3track.mass(); - if (h3track.pt() >= ptMinTOF && h3track.pt() <= ptMaxTOF && lnnCand.mass2TrTOF >= TrTOFMass2Cut) { - h3HSignalPtTOF->Fill(chargeFactor * h3track.pt(), beta); - hNsigma3HSelTOF->Fill(chargeFactor * h3track.pt(), h3track.tofNSigmaTr()); - h3HMassPtTOF->Fill(chargeFactor * h3track.pt(), lnnCand.mass2TrTOF); + if (lnnCand.mass2TrTOF < TrTOFMass2Cut) { + continue; } } int nCand = 0; try { - nCand = fitter.process(h3TrackCov, negTrackCov); + nCand = fitter.process(posTrackCov, negTrackCov); } catch (...) { LOG(error) << "Exception caught in DCA fitter process call!"; continue; @@ -459,12 +465,11 @@ struct lnnRecoTask { // if survived all selections, propagate decay daughters to PV gpu::gpustd::array dcaInfo; + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, h3PropTrack, 2.f, fitter.getMatCorrType(), &dcaInfo); + lnnCand.h3DCAXY = dcaInfo[0]; - o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, h3TrackCov, 2.f, fitter.getMatCorrType(), &dcaInfo); - lnnCand.isMatter ? lnnCand.h3DCAXY = dcaInfo[0] : lnnCand.piDCAXY = dcaInfo[0]; - - o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, negTrackCov, 2.f, fitter.getMatCorrType(), &dcaInfo); - lnnCand.isMatter ? lnnCand.piDCAXY = dcaInfo[0] : lnnCand.h3DCAXY = dcaInfo[0]; + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, piPropTrack, 2.f, fitter.getMatCorrType(), &dcaInfo); + lnnCand.piDCAXY = dcaInfo[0]; // finally, push back the candidate lnnCand.isReco = true; @@ -472,6 +477,16 @@ struct lnnRecoTask { lnnCand.negTrackID = negTrack.globalIndex(); lnnCandidates.push_back(lnnCand); + + // Fill QA histograms + hdEdx3HSel->Fill(chargeFactor * lnnCand.mom3HTPC, h3track.tpcSignal()); + hNsigma3HSel->Fill(chargeFactor * lnnCand.mom3HTPC, lnnCand.nSigma3H); + hDCAxy3H->Fill(h3track.pt(), h3track.dcaXY()); + if (h3track.hasTOF()) { + h3HSignalPtTOF->Fill(chargeFactor * h3track.pt(), beta); + hNsigma3HSelTOF->Fill(chargeFactor * h3track.pt(), h3track.tofNSigmaTr()); + h3HMassPtTOF->Fill(chargeFactor * h3track.pt(), lnnCand.mass2TrTOF); + } } }