BmnHypNuclPairFinder.cxx

Go to the documentation of this file.

 
#include "BmnHypNuclPairFinder.h"
 
#include "TH1F.h"
 
BmnHypNuclPairFinder::BmnHypNuclPairFinder() {}
 
BmnHypNuclPairFinder::BmnHypNuclPairFinder(TString hyp, Int_t run, Bool_t exp)
{
    if (hyp.Contains("H3L") || hyp.Contains("h3l")) {
        massPos = fMassHe3;
        massNeg = fMassPi;
    } else if (hyp.Contains("H4L") || hyp.Contains("h4l")) {
        massPos = fMassHe4;
        massNeg = fMassPi;
    } else {
        massPos = fMassHe3;
        massNeg = fMassPi;
    }
    fEventCounter = 0;
    fKalman = nullptr;
    fParticlePair = nullptr;
    fVertex = nullptr;
    fTrigger = nullptr;
    fStsHits = nullptr;
    fUpperGemClusters = nullptr;
    fLowerGemClusters = nullptr;
    fUpperFsdClusters = nullptr;
    fLowerFsdClusters = nullptr;
    fStsTracks = nullptr;
    fGlobTracks = nullptr;
    fTof700Hits = nullptr;
    fDstHeader = nullptr;
    isExp = exp;
    runId = run;
 
    nSt = 11;
    nSect = {6, 10, 14, 18, 8, 8, 8, 8, 8, 8, 8};
    nMomSlices = 6;
    minMomBound = {0.0, 0.5, 1.0, 1.5, 2.0, 3.0};
    maxMomBound = {0.5, 1.0, 1.5, 2.0, 3.0, 20.0};
 
    a = new Double_t**[nSt];
    b = new Double_t**[nSt];
    for (Int_t i = 0; i < nSt; ++i) {
        a[i] = new Double_t*[nSect[i]];
        b[i] = new Double_t*[nSect[i]];
        for (Int_t j = 0; j < nSect[i]; ++j) {
            a[i][j] = new Double_t[nMomSlices];
            b[i][j] = new Double_t[nMomSlices];
            for (Int_t k = 0; k < nMomSlices; ++k) {
                a[i][j][k] = 0.0;
                b[i][j][k] = 0.0;
            }
        }
    }
 
    m2CutLowKplsTof700 = new TF1("m2CutLowKplsTof400", "-0.00970*x*x-0.0070*x+0.239", 0, 3.5);
    m2CutUpKplsTof700 = new TF1("m2CutUpKplsTof400", "0.00462*x*x+0.0357*x+0.279", 0, 3.5);
    m2CutLowKplsTof400 = new TF1("m2CutLowKplsTof700", "-0.02108*x*x-0.0278*x+0.255", 0, 3.5);
    m2CutUpKplsTof400 = new TF1("m2CutUpKplsTof700", "0.00284*x*x+0.0434*x+0.256", 0, 3.5);
 
    if (massPos > 3) {   // He4
        fit_dP_He = new TF1("fit_dP_He", "[0]+[1]*(1+(3.728/x)^2)^[2]");
        fit_dP_He->SetParameters(0.2010, -0.2875, 0.6414);
    } else {   // He3
        fit_dP_He = new TF1("fit_dP_He", "[0]+[1]*(1+(2.809/x)^2)^[2]");
        fit_dP_He->SetParameters(0.0731, -0.1687, 0.8723);
    }
 
    alephFit_d = new TF1("alephFit_d", "[0]/(x^([3]))*([1]-x^([3])-TMath::Log([2]+1/(x^([4]))))", 0, 14.5);
    const Float_t fpar1[5] = {-591.879, -0.537493, 0.338613, -0.00845327,
                              3.16787};   // parameters from fitting hdEdxHe4
    alephFit_d->SetParameters(fpar1[0], fpar1[1], fpar1[2], fpar1[3], fpar1[4]);
 
    if (massPos > 3) {   // He4
        m2CutLow = new TF1("m2CutLow", "-0.0236*x*x+0.0240*x+3.1896", 0, 30);
        m2CutUp = new TF1("m2CutUp", "0.0169*x*x+0.1306*x+4.0525", 0, 30);
    } else {   // He3
        m2CutLow = new TF1("m2CutLow", "-0.0034*x*x-0.0360*x+1.7426", 0, 30);
        m2CutUp = new TF1("m2CutUp", "0.0128*x*x-0.1060*x+2.6261", 0, 30);
    }
}
 
BmnHypNuclPairFinder::~BmnHypNuclPairFinder()
{
    delete m2CutLowKplsTof400;
    delete m2CutUpKplsTof400;
    delete m2CutLowKplsTof700;
    delete m2CutUpKplsTof700;
    delete fit_dP_He;
    delete alephFit_d;
    delete m2CutLow;
    delete m2CutUp;
    if (a) {
        for (Int_t i = 0; i < nSt; ++i) {
            for (Int_t j = 0; j < nSect[i]; ++j) {
                delete[] a[i][j];
            }
            delete[] a[i];
        }
        delete[] a;
        a = nullptr;
    }
    if (b) {
        for (Int_t i = 0; i < nSt; ++i) {
            for (Int_t j = 0; j < nSect[i]; ++j) {
                delete[] b[i][j];
            }
            delete[] b[i];
        }
        delete[] b;
        b = nullptr;
    }
}
 
void BmnHypNuclPairFinder::Analysis()
{
    Float_t VPX = fVertex->GetX();
    Float_t VPY = fVertex->GetY();
    Float_t VPZ = fVertex->GetZ();
 
    fParInVPPos.clear();
    fParInVPNeg.clear();
    fM2Pos700.clear();
    fChi2NDFPos.clear();
    fChi2NDFNeg.clear();
    fDedxPos.clear();
    fNGemHitsNeg.clear();
    fNGemHitsPos.clear();
    fNFsdHitsNeg.clear();
    fNFsdHitsPos.clear();
 
    Bool_t strangeEvent = kFALSE;
 
    for (Int_t i = 0; i < fGlobTracks->GetEntriesFast(); i++) {
        BmnGlobalTrack* gl = (BmnGlobalTrack*)fGlobTracks->UncheckedAt(i);
        if (IsKaPlus(gl)) {
            strangeEvent = kTRUE;
            break;
        }
    }
 
    if (!strangeEvent)
        return;
 
    for (Int_t i = 0; i < fGlobTracks->GetEntriesFast(); i++) {
        BmnGlobalTrack* gl = (BmnGlobalTrack*)fGlobTracks->UncheckedAt(i);
        CbmStsTrack* sts = (CbmStsTrack*)fStsTracks->At(gl->GetGemTrackIndex());
 
        if (sts->GetB() < -0.1)   // TMVA cut for good tracks
            continue;
 
        if (sts->GetNStsHits() < 4)
            continue;
 
        FairTrackParam par = *(sts->GetParamFirst());
        const Float_t tZ = -0.01942;   // z-coordinate of target for RUN-8
        if (fKalman->TGeoTrackPropagate(&par, tZ, fPDG1, nullptr, nullptr) == kBMNERROR)
            continue;
 
        Int_t nGEM = 0;
        Int_t nFSD = 0;
        GetNHits(sts, nGEM, nFSD);
 
        Float_t mom = gl->GetP();
        if (mom > 0) {   //                        He3/He4
 
            Float_t dedx = GetdEdx(sts);   // gl->GetdQdNLower();
            Float_t m2 = gl->GetMass2(2);
            // ============================================================
            // Отбор по dE/dx.
            //
            // 1. dedx < 0 — GetdEdx() не смог измерить (мало хитов,
            //    импульс вне диапазона и т.п.). Отбрасываем трек
            //    явно.
            //
            // 2. alephFit_d определена на [0, 14.5]. Вне этого
            //    диапазона Eval() возвращает 0 или мусор, что даёт
            //    inf/nan при делении. Отбрасываем такие треки явно.
            //
            // 3. Eval() может вернуть <= 0 для некоторых mom —
            //    защищаемся от деления на ноль/отрицательное.
            // ============================================================
            if (dedx < 0.0)
                continue;
 
            if (mom > 14.5)
                continue;
 
            const Float_t expectedDedx = alephFit_d->Eval(mom);
            if (expectedDedx <= 0.0)
                continue;
 
            if (dedx / expectedDedx < 2.25)
                continue;
 
            fDedxPos.push_back(dedx);
            fM2Pos700.push_back(m2);
 
            fParInVPPos.push_back(par);
            fChi2NDFPos.push_back(sts->GetChi2() / sts->GetNDF());
            fNGemHitsPos.push_back(nGEM);
            fNFsdHitsPos.push_back(nFSD);
        } else {   //                               pi-
            fParInVPNeg.push_back(par);
            fChi2NDFNeg.push_back(sts->GetChi2() / sts->GetNDF());
            fNGemHitsNeg.push_back(nGEM);
            fNFsdHitsNeg.push_back(nFSD);
        }
    }
 
    for (size_t i = 0; i < fParInVPPos.size(); ++i) {
        for (size_t j = 0; j < fParInVPNeg.size(); ++j) {
 
            FairTrackParam VP1 = fParInVPPos[i];
            FairTrackParam VP2 = fParInVPNeg[j];
            FairTrackParam V01 = fParInVPPos[i];
            FairTrackParam V02 = fParInVPNeg[j];
 
            const Float_t kFindV0Range = 30.0f;   // см
 
            Float_t V0Z = FindV0ByVirtualPlanes(&V01, &V02, VPZ, kFindV0Range);
            if (V0Z > 15 || V0Z < 0)
                continue;
 
            V01 = fParInVPPos[i];
            V02 = fParInVPNeg[j];
 
            fKalman->TGeoTrackPropagate(&V01, V0Z, fPDG1, nullptr, nullptr, kTRUE);
            fKalman->TGeoTrackPropagate(&V02, V0Z, fPDG2, nullptr, nullptr, kTRUE);
 
            Float_t V0X1 = V01.GetX();
            Float_t V0X2 = V02.GetX();
            Float_t V0Y1 = V01.GetY();
            Float_t V0Y2 = V02.GetY();
            Float_t VPX1 = VP1.GetX();
            Float_t VPX2 = VP2.GetX();
            Float_t VPY1 = VP1.GetY();
            Float_t VPY2 = VP2.GetY();
 
            Float_t V0X = .5 * (V0X1 + V0X2);
            Float_t V0Y = .5 * (V0Y1 + V0Y2);
            Float_t path = Sqrt(Sq(V0X - VPX) + Sq(V0Y - VPY) + Sq(V0Z - VPZ));
 
            BmnHypNuclPair partPair;
 
            partPair.SetM2(fM2Pos700[i]);
            partPair.SetChi2NDF(fChi2NDFPos[i], fChi2NDFNeg[j]);
 
            Float_t cV01_00 = V01.GetCovariance(0, 0);
            Float_t cV01_01 = V01.GetCovariance(0, 1);
            Float_t cV01_11 = V01.GetCovariance(1, 1);
 
            Float_t cV02_00 = V02.GetCovariance(0, 0);
            Float_t cV02_01 = V02.GetCovariance(0, 1);
            Float_t cV02_11 = V02.GetCovariance(1, 1);
 
            Float_t cVP1_00 = VP1.GetCovariance(0, 0);
            Float_t cVP1_01 = VP1.GetCovariance(0, 1);
            Float_t cVP1_11 = VP1.GetCovariance(1, 1);
 
            Float_t cVP2_00 = VP2.GetCovariance(0, 0);
            Float_t cVP2_01 = VP2.GetCovariance(0, 1);
            Float_t cVP2_11 = VP2.GetCovariance(1, 1);
 
            Float_t dx12 = V0X2 - V0X1;
            Float_t dy12 = V0Y2 - V0Y1;
            Float_t chi2_DCA12 = ComputeChi2_2D(dx12, dy12, cV01_00, cV01_01, cV01_11, cV02_00, cV02_01, cV02_11);
            if (chi2_DCA12 < 0.0)
                continue;
 
            Float_t dx1 = VPX1 - VPX;
            Float_t dy1 = VPY1 - VPY;
 
            Float_t cVtx_00 = fVertex->GetCovariance(0, 0);
            Float_t cVtx_01 = 0.0f;
            Float_t cVtx_11 = fVertex->GetCovariance(1, 1);
 
            Float_t chi2_DCA1 = ComputeChi2_2D(dx1, dy1, cVP1_00, cVP1_01, cVP1_11, cVtx_00, cVtx_01, cVtx_11);
            if (chi2_DCA1 < 0.0)
                continue;
            Float_t dx2 = VPX2 - VPX;
            Float_t dy2 = VPY2 - VPY;
            Float_t chi2_DCA2 = ComputeChi2_2D(dx2, dy2, cVP2_00, cVP2_01, cVP2_11, cVtx_00, cVtx_01, cVtx_11);
            if (chi2_DCA2 < 0.0)
                continue;
            partPair.SetPath(path);
 
            partPair.SetNHitsPair(fNFsdHitsPos[i], fNGemHitsPos[i], fNFsdHitsNeg[j], fNGemHitsNeg[j]);
 
            Float_t p1 = Abs(1. / V01.GetQp());
            p1 *= 2.0;   // p1 is actually p1/Q, Q is 2 for helium
            p1 -= fit_dP_He->Eval(p1);
            Float_t p2 = Abs(1. / V02.GetQp());
            Float_t Tx1 = V01.GetTx();
            Float_t Ty1 = V01.GetTy();
            Float_t Tx2 = V02.GetTx();
            Float_t Ty2 = V02.GetTy();
            Float_t A1 = 1. / Sqrt(Tx1 * Tx1 + Ty1 * Ty1 + 1);
            Float_t A2 = 1. / Sqrt(Tx2 * Tx2 + Ty2 * Ty2 + 1);
            Float_t pz1 = p1 * A1;
            Float_t pz2 = p2 * A2;
            Float_t px1 = Tx1 * pz1;
            Float_t px2 = Tx2 * pz2;
            Float_t py1 = Ty1 * pz1;
            Float_t py2 = Ty2 * pz2;
 
            partPair.SetMomXYZ1(px1, py1, pz1);
            partPair.SetMomXYZ2(px2, py2, pz2);
 
            Float_t PzPart0 = pz1 + pz2;
            Float_t PxPart0 = px1 + px2;
            Float_t PyPart0 = py1 + py2;
            Float_t TxPart0 = PxPart0 / PzPart0;
            Float_t TyPart0 = PyPart0 / PzPart0;
            Float_t xPart0 = TxPart0 * (VPZ - V0Z) + V0X;
            Float_t yPart0 = TyPart0 * (VPZ - V0Z) + V0Y;
            Float_t dx0 = xPart0 - VPX;
            Float_t dy0 = yPart0 - VPY;
 
            Float_t chi2_DCA0 = ComputeChi2_2D(dx0, dy0, cVtx_00, cVtx_01, cVtx_11, cVtx_00, cVtx_01, cVtx_11);
            if (chi2_DCA0 < 0.0)
                continue;
            partPair.SetChi2DCA1(chi2_DCA1);
            partPair.SetChi2DCA2(chi2_DCA2);
            partPair.SetChi2DCA12(chi2_DCA12);
            partPair.SetChi2DCA0(chi2_DCA0);
 
            partPair.SetDedx(fDedxPos[i]);
 
            const Double_t E1 =
                TMath::Sqrt(static_cast<Double_t>(px1) * px1 + static_cast<Double_t>(py1) * py1
                            + static_cast<Double_t>(pz1) * pz1 + static_cast<Double_t>(massPos) * massPos);
            const Double_t E2 =
                TMath::Sqrt(static_cast<Double_t>(px2) * px2 + static_cast<Double_t>(py2) * py2
                            + static_cast<Double_t>(pz2) * pz2 + static_cast<Double_t>(massNeg) * massNeg);
 
            const Double_t E = E1 + E2;
            const Double_t Px = static_cast<Double_t>(px1) + px2;
            const Double_t Py = static_cast<Double_t>(py1) + py2;
            const Double_t Pz = static_cast<Double_t>(pz1) + pz2;
 
            const Double_t m2 = E * E - Px * Px - Py * Py - Pz * Pz;
            const Double_t m = (m2 > 0.0) ? TMath::Sqrt(m2) : 0.0;
 
            partPair.SetInvMass(static_cast<Float_t>(m));
 
            partPair.SetNTrInEvent(fVertex->GetNTracks());
 
            new ((*fParticlePair)[fParticlePair->GetEntriesFast()]) BmnHypNuclPair(partPair);
        }
    }
}
 
InitStatus BmnHypNuclPairFinder::Init()
{
    cout << "\nBmnHypNuclPairFinder::Init()" << endl;
 
    FairRootManager* ioman = FairRootManager::Instance();
 
    fStsTracks = (TClonesArray*)ioman->GetObject("StsVector");
    fGlobTracks = (TClonesArray*)ioman->GetObject("BmnGlobalTrack");
    fStsHits = (TClonesArray*)ioman->GetObject("StsHit");
    fVertex = (CbmVertex*)ioman->GetObject("MpdVertex.");
    fTrigger = (BmnTrigInfoDst*)ioman->GetObject("BmnTrigInfo.");
    fUpperGemClusters = (TClonesArray*)ioman->GetObject("BmnGemUpperCluster");
    fLowerGemClusters = (TClonesArray*)ioman->GetObject("BmnGemLowerCluster");
    fUpperFsdClusters = (TClonesArray*)ioman->GetObject("BmnSiliconUpperCluster");
    fLowerFsdClusters = (TClonesArray*)ioman->GetObject("BmnSiliconLowerCluster");
    fTof700Hits = (TClonesArray*)ioman->GetObject("BmnTof700Hit");
    fDstHeader = (DstEventHeader*)ioman->GetObject("DstEventHeader.");
 
    fParticlePair = new TClonesArray("BmnHypNuclPair");
    ioman->Register("pairs", "HYP", fParticlePair, kTRUE);
 
    fKalman = new BmnKalmanFilter();
 
    TString dir = getenv("VMCWORKDIR");
    if (dir.IsNull()) {
        Error("Init", "VMCWORKDIR is not set — cannot locate calibration.root");
        return kERROR;
    }
 
    TString calibPath = dir + "/input/calibration.root";
    TFile* calibFile = TFile::Open(calibPath);
    if (!calibFile || calibFile->IsZombie()) {
        Error("Init", "Cannot open calibration file: %s", calibPath.Data());
        return kERROR;
    }
 
    for (Int_t i = 0; i < nSt; ++i) {
        for (Int_t j = 0; j < nSect[i]; ++j) {
            for (Int_t k = 0; k < nMomSlices; ++k) {
                TH1F* hA = (TH1F*)calibFile->Get(Form("hA_run_%d_%d_%d", i, j, k));
                TH1F* hB = (TH1F*)calibFile->Get(Form("hB_run_%d_%d_%d", i, j, k));
 
                if (!hA || !hB) {
                    Warning("Init",
                            "Missing calibration histograms for (i=%d, j=%d, k=%d) — "
                            "setting a=b=0",
                            i, j, k);
                    a[i][j][k] = 0.0;
                    b[i][j][k] = 0.0;
                    continue;
                }
 
                a[i][j][k] = hA->GetBinContent(hA->FindBin(runId));
                b[i][j][k] = hB->GetBinContent(hB->FindBin(runId));
            }
        }
    }
 
    calibFile->Close();
    delete calibFile;
    calibFile = nullptr;
 
    constexpr size_t kReserveSize = 128;
 
    fParInVPPos.reserve(kReserveSize);
    fParInVPNeg.reserve(kReserveSize);
    fM2Pos700.reserve(kReserveSize);
    fChi2NDFPos.reserve(kReserveSize);
    fChi2NDFNeg.reserve(kReserveSize);
    fDedxPos.reserve(kReserveSize);
    fNGemHitsNeg.reserve(kReserveSize);
    fNGemHitsPos.reserve(kReserveSize);
    fNFsdHitsNeg.reserve(kReserveSize);
    fNFsdHitsPos.reserve(kReserveSize);
 
    return kSUCCESS;
}
 
// -------------------------------------------------------------------
 
void BmnHypNuclPairFinder::Exec(Option_t* option)
{
    fParticlePair->Delete();
 
    fEventCounter++;
 
    if (fEventCounter % 1000 == 0)
        cout << "Event# " << fEventCounter << endl;
 
    if (isExp) {
        if (!(fTrigger->IsTriggerBitTrueAR(7)))   // bit #7 is for CCT2
            return;
        if (fVertex->GetNTracks() < 2)
            return;
        if (fVertex->GetZ() < -0.5 || fVertex->GetZ() > +0.5)
            return;
 
        Float_t centerOfTargetX = 0.40;
        Float_t centerOfTargetY = 0.15;
        Float_t radiusOfTarget = 1.2;
 
        if (TMath::Sqrt(TMath::Sq(fVertex->GetX() - centerOfTargetX) + TMath::Sq(fVertex->GetY() - centerOfTargetY))
            > radiusOfTarget)
            return;
    }
    Analysis();
}
// -------------------------------------------------------------------
 
void BmnHypNuclPairFinder::Finish()
{
    delete fKalman;
    // delete fMagField;
    cout << "\n-I- [BmnHypNuclPairFinder::Finish] " << endl;
}
 
Float_t BmnHypNuclPairFinder::FindV0ByVirtualPlanes(FairTrackParam* par1,
                                                    FairTrackParam* par2,
                                                    Float_t z_0,
                                                    Float_t range)
{
    const FairTrackParam par1_initial = *par1;
    const FairTrackParam par2_initial = *par2;
 
    const Int_t nPlanes = 5;   // FIXME: можно адаптивно (5 → 3)
 
    Float_t minDist2 = numeric_limits<Float_t>::max();
    Float_t minZ = z_0;
 
    while (range >= 0.1) {
        const Float_t zMax = z_0 + range;
        const Float_t zMin = z_0 - range;
        const Float_t zStep = (zMax - zMin) / nPlanes;
 
        *par1 = par1_initial;
        *par2 = par2_initial;
 
        minDist2 = numeric_limits<Float_t>::max();
        minZ = z_0;
 
        for (Int_t iPlane = 0; iPlane < nPlanes; ++iPlane) {
            const Float_t z = zMax - iPlane * zStep;
 
            FairTrackParam p1 = par1_initial;
            FairTrackParam p2 = par2_initial;
 
            BmnStatus stat1 = fKalman->TGeoTrackPropagate(&p1, z, fPDG1, nullptr, nullptr, kTRUE);
            BmnStatus stat2 = fKalman->TGeoTrackPropagate(&p2, z, fPDG2, nullptr, nullptr, kTRUE);
 
            if (stat1 == kBMNERROR || stat2 == kBMNERROR)
                return DBL_MAX;
 
            const Float_t dx = p1.GetX() - p2.GetX();
            const Float_t dy = p1.GetY() - p2.GetY();
            const Float_t d2 = dx * dx + dy * dy;
 
            if (d2 < minDist2) {
                minDist2 = d2;
                minZ = z;
            }
        }
 
        if (TMath::Abs(minZ - z_0) < 0.01)
            break;
 
        z_0 = minZ;
        range *= 0.5;
    }
 
    return z_0;
}
 
Float_t BmnHypNuclPairFinder::GetdEdx(CbmStsTrack* tr)
{
    if (!isExp)
        return -1.0f;
 
    // TODO: когда Q появится в BmnGlobalTrack,
    // добавить проверку if (Q != 2) return -1.0f;
    // Это сразу отсечёт π⁻, K⁺, p и оставит только He3/He4. Этот метод будет не нужен.
 
    Float_t p = 1.0 / tr->GetParamFirst()->GetQp();
    if (p < 0.0 || p > 20)
        return -1.0f;
 
    Int_t momIdx = -1;
    for (Int_t iP = 0; iP < nMomSlices; ++iP) {
        if (p > minMomBound[iP] && p < maxMomBound[iP]) {
            momIdx = iP;
            break;
        }
    }
    if (momIdx == -1)
        return -1.0f;
 
    // ============================================================
    // Собираем сигналы в статический массив.
    //
    // std::set (старая версия) даёт автоматическую сортировку,
    // но платит за это аллокацией на каждый элемент и O(log n)
    // на вставку. std::array + std::nth_element — zero аллокаций
    // и O(n) в среднем.
    //
    // Логика truncated mean сохранена: после nth_element первые
    // usedHits элементов — наименьшие, они и усредняются.
    // Ландау-хвост (наибольшие) отбрасывается.
    //
    // ВАЖНО: в отличие от std::set, std::array НЕ удаляет
    // дубликаты. Если два хита дают одинаковый dE/dx, оба
    // учитываются. Для truncated mean это не критично —
    // важны значения, а не число повторов.
    // ============================================================
    constexpr Int_t kMaxHits = 64;
    std::array<Float_t, kMaxHits> signalsLow;
    Int_t nSignals = 0;
 
    Int_t nHits = tr->GetNStsHits();
    for (Int_t hitIdx = 0; hitIdx < nHits; ++hitIdx) {
        CbmStsHit* hit = (CbmStsHit*)fStsHits->At(tr->GetStsHitIndex(hitIdx));
 
        Int_t sect = hit->GetSectorNr() - 1;
        Int_t st = hit->GetStationNr() - 1;
 
        StripCluster low;
        if (hit->GetSystemId() == kGEM) {
            low = *((StripCluster*)fLowerGemClusters->At(hit->fDigiF % 1000000));
        } else if (hit->GetSystemId() == kSILICON) {
            low = *((StripCluster*)fLowerFsdClusters->At(hit->fDigiF % 1000000));
        } else {
            continue;
        }
 
        if (nSignals >= kMaxHits)
            break;   // защита от переполнения
 
        signalsLow[nSignals++] = low.TotalSignal * a[st][sect][momIdx] + b[st][sect][momIdx];
    }
 
    if (nSignals == 0)
        return -1.0f;
 
    Int_t usedHits = (nHits == 3)    ? 2
                     : (nHits == 4)  ? 2
                     : (nHits == 5)  ? 3
                     : (nHits == 6)  ? 4
                     : (nHits == 7)  ? 4
                     : (nHits == 8)  ? 5
                     : (nHits == 9)  ? 6
                     : (nHits == 10) ? 6
                     : (nHits > 10)  ? 7
                                     : 0;
 
    if (usedHits <= 0)
        return -1.0f;
 
    if (usedHits > nSignals)
        usedHits = nSignals;
 
    // nth_element: первые usedHits — наименьшие (не отсортированы,
    // но это и не нужно — нам важна только сумма)
    std::nth_element(signalsLow.begin(), signalsLow.begin() + usedHits, signalsLow.begin() + nSignals);
 
    Float_t totSigLow = 0.0f;
    for (Int_t i = 0; i < usedHits; ++i)
        totSigLow += signalsLow[i];
 
    return totSigLow / usedHits;
}
 
void BmnHypNuclPairFinder::GetNHits(CbmStsTrack* tr, Int_t& nGEM, Int_t& nFSD)
{
    nGEM = 0;
    nFSD = 0;
    for (Int_t hitIdx = 0; hitIdx < tr->GetNStsHits(); ++hitIdx) {
        CbmStsHit* hit = (CbmStsHit*)fStsHits->At(tr->GetStsHitIndex(hitIdx));
 
        if (hit->GetSystemId() == kGEM)
            nGEM++;
        else if (hit->GetSystemId() == kSILICON)
            nFSD++;
    }
}
 
Bool_t BmnHypNuclPairFinder::IsKaPlus(BmnGlobalTrack* gl)
{
 
    Float_t m2tof400 = gl->GetMass2(1);   // M2 by tof400
    Float_t m2tof700 = gl->GetMass2(2);   // M2 by tof700
    Float_t p = gl->GetP();
 
    if (p > 3.5 || p < 0.0)
        return kFALSE;
 
    // BOTH TOFs
    return (m2tof400 < m2CutUpKplsTof400->Eval(p) && m2tof400 > m2CutLowKplsTof400->Eval(p))
           || (m2tof700 < m2CutUpKplsTof700->Eval(p) && m2tof700 > m2CutLowKplsTof700->Eval(p));
 
    // TOF400
    // return (m2tof400 < m2CutUpKplsTof400->Eval(p) && m2tof400 > m2CutLowKplsTof400->Eval(p));
 
    // TOF700
    // return (m2tof700 < m2CutUpKplsTof700->Eval(p) && m2tof700 > m2CutLowKplsTof700->Eval(p));
}
 
ClassImp(BmnHypNuclPairFinder);