BmnHypNuclPairFinder.h

Go to the documentation of this file.

 
#ifndef BMNHYPNUCLPAIRFINDER_H
#define BMNHYPNUCLPAIRFINDER_H 1
 
#include "BmnFieldMap.h"
#include "BmnGlobalTrack.h"
#include "BmnHypNuclPair.h"
#include "BmnKalmanFilter.h"
#include "BmnNewFieldMap.h"
#include "BmnStripData.h"
#include "BmnTofHit.h"
#include "BmnTrigInfoDst.h"
#include "CbmStsHit.h"
#include "CbmStsTrack.h"
#include "CbmVertex.h"
#include "DstEventHeader.h"
#include "FairTask.h"
#include "TFile.h"
#include "TGeoManager.h"
#include "TRandom.h"
#include "UniRun.h"
 
using namespace std;
using namespace TMath;
 
class BmnHypNuclPairFinder : public FairTask
{
  private:
    UInt_t fEventCounter;   //! event counter
 
    TClonesArray* fStsHits;
    TClonesArray* fUpperGemClusters;
    TClonesArray* fLowerGemClusters;
    TClonesArray* fUpperFsdClusters;
    TClonesArray* fLowerFsdClusters;
    TClonesArray* fStsTracks;
    TClonesArray* fGlobTracks;
    TClonesArray* fTof700Hits;
    CbmVertex* fVertex;
    BmnTrigInfoDst* fTrigger;
    DstEventHeader* fDstHeader;
 
    BmnKalmanFilter* fKalman;
 
    TClonesArray* fParticlePair;
 
    Float_t massPos;
    Float_t massNeg;
 
    const Int_t fPDG1 = 2212;   // we use proton insted of he3, because kalman cannot work with ions; 1000020030; //He3
    const Int_t fPDG2 = -211;   // pi-minus
 
    const Float_t fMassHe4 = 3.7284013;
    const Float_t fMassHe3 = 2.8094135;
    const Float_t fMassPi = 0.1395704;
    const Float_t fMassP = 0.9382721;
    const Float_t fMassD = 1.8756129;
 
    Bool_t isExp;
    Int_t runId;
 
    // TF1s for parabolic cuts on K+
    TF1* m2CutLowKplsTof700;
    TF1* m2CutUpKplsTof700;
    TF1* m2CutLowKplsTof400;
    TF1* m2CutUpKplsTof400;
 
    // Parametrization for dP = Pout - Pin over Pout.
    // Needed to correct momentum (KF doesn't work correctly for Q=2, default value is Q=1)
    // This parametrization got in macro ELossTest.C (on VM)
    TF1* fit_dP_He;
 
    // Cut from Matvey Klimenko by 10mln events
    // 4He/d data were divided by rough cut to d and 4He
    // d data were fitted by aleph (medians vs p prepared first)
    TF1* alephFit_d;
 
    // m^2 cut to select He-4 (low)
    TF1* m2CutLow;
    // m^2 cut to select He-4 (up)
    TF1* m2CutUp;
 
    Int_t nSt;
    vector<Int_t> nSect;
    Int_t nMomSlices;
    vector<Float_t> minMomBound;
    vector<Float_t> maxMomBound;
 
    Double_t*** a;
    Double_t*** b;
 
    vector<FairTrackParam> fParInVPPos;
    vector<FairTrackParam> fParInVPNeg;
    vector<Float_t> fM2Pos700;
    vector<Float_t> fChi2NDFPos;
    vector<Float_t> fChi2NDFNeg;
    vector<Float_t> fDedxPos;
    vector<Short_t> fNGemHitsNeg;
    vector<Short_t> fNGemHitsPos;
    vector<Short_t> fNFsdHitsNeg;
    vector<Short_t> fNFsdHitsPos;
 
  public:
    BmnHypNuclPairFinder();
    // hyp = "H3L"/"h3l" or "H4L"/"h4l"
    BmnHypNuclPairFinder(TString hyp, Int_t run, Bool_t exp = kTRUE);
    virtual ~BmnHypNuclPairFinder();
 
    virtual void Exec(Option_t* option);
    virtual InitStatus Init();
    virtual void Finish();
 
  private:
    // ============================================================
    // Sentinel-значение для невалидного χ².
    // ============================================================
    static constexpr Float_t kInvalidChi2 = -1.0f;
    // ============================================================
    // Вспомогательная функция: 2D χ² с учётом корреляций.
    // dx, dy   — разности координат (например, x1 - x2)
    // C1, C2   — ковариационные матрицы (2x2, представленные как
    //            массивы [c00, c01, c11])
    // Возвращает χ² = Δᵀ · (C1 + C2)⁻¹ · Δ, или kInvalidChi2
    // при вырожденной матрице.
    //
    // ВАЖНО: внутренние вычисления — в Double_t. Это критично для
    // треков с хорошим разрешением (σ ~ 0.01 см), где det матрицы
    // ~ 1e-8 и Float_t теряет точность (вплоть до отрицательного det).
    // ============================================================
    static inline Float_t ComputeChi2_2D(Float_t dx,
                                         Float_t dy,
                                         Float_t c1_00,
                                         Float_t c1_01,
                                         Float_t c1_11,
                                         Float_t c2_00,
                                         Float_t c2_01,
                                         Float_t c2_11)
    {
        // Суммарная ковариация — в Double_t
        const Double_t s00 = static_cast<Double_t>(c1_00) + c2_00;
        const Double_t s01 = static_cast<Double_t>(c1_01) + c2_01;
        const Double_t s11 = static_cast<Double_t>(c1_11) + c2_11;
 
        const Double_t det = s00 * s11 - s01 * s01;
        if (det <= 0.0)
            return kInvalidChi2;   // вырожденная матрица
 
        // Обратная матрица 2x2
        const Double_t inv00 = s11 / det;
        const Double_t inv01 = -s01 / det;
        const Double_t inv11 = s00 / det;
 
        // χ² = Δᵀ · Cov⁻¹ · Δ
        const Double_t dxd = dx;
        const Double_t dyd = dy;
        const Double_t chi2 = dxd * dxd * inv00 + 2.0 * dxd * dyd * inv01 + dyd * dyd * inv11;
 
        // Защита от отрицательного χ² из-за численных ошибок
        // (теоретически невозможно, но double round-off может дать -1e-15)
        return (chi2 >= 0.0) ? static_cast<Float_t>(chi2) : 0.0f;
    }
 
    void Analysis();
 
    Float_t FindV0ByVirtualPlanes(FairTrackParam* par1, FairTrackParam* par2, Float_t z_0, Float_t range = 50.);
    Float_t GetdEdx(CbmStsTrack* tr);
    void GetNHits(CbmStsTrack* tr, Int_t& nGEM, Int_t& nFSD);
    Float_t GetSigXTof700He3(FairTrackParam* par);
    Float_t GetSigYTof700He3(FairTrackParam* par);
    Float_t GetDyTof700(FairTrackParam* par);
    Float_t GetDxTof700(FairTrackParam* par);
 
    // Method to check the presence of strangness in event by presence of K+ or K-
    Bool_t IsKaPlus(BmnGlobalTrack* gl);
 
    ClassDef(BmnHypNuclPairFinder, 0)
};
 
#endif