BmnUpstreamTracking.h

Go to the documentation of this file.

// Author: Vasilisa Lenivenko (VBLHEP) <vasilisa@jinr.ru> 2019-07-18
 
////////////////////////////////////////////////////////////////////////////////
//                                                                            //
// BmnUpstreamTracking                                                      //
//                                                                            //
// Implementation of an algorithm developed by                                //
// V.Lenivenko and V.Paltchik                                                 //
// to the BmnRoot software                                                    //
//                                                                            //
// The algorithm serves for searching for tracks                              //
// in the Silicon detectors & MWPCs                                           //
// of the configuration SRC of BM@N experiment                                //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////
 
#ifndef BMNUPSTREAMTRACKING_H
#define BMNUPSTREAMTRACKING_H 1
 
#include "BmnMwpcHit.h"
#include "BmnMwpcSegment.h"
#include "BmnMwpcTrack.h"
#include "BmnSiliconHit.h"
#include "BmnSiliconTrackFinder.h"
#include "BmnTrack.h"
#include "FairMCPoint.h"
#include "FairTask.h"
#include "FairTrackParam.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TList.h"
 
#include <TClonesArray.h>
#include <TMath.h>
#include <TString.h>
#include <TVector2.h>
#include <fstream>
#include <vector>
using namespace std;
using namespace TMath;
 
struct Smatch
{
    Int_t Ind1 = -1;    // Si
    Int_t Ind2 = -1;    // Pctr
    Int_t Inds2 = -1;   // seg2
    Int_t Inds3 = -1;   // seg3
    Int_t Nhits1 = -1;
    Int_t Nhits2 = -1;
    Double_t Z1 = -999.;
    Double_t Z2 = -999.;
    Double_t Zs2 = -999.;
    Double_t Zs3 = -999.;
    Double_t distX = -1.;
    Double_t distY = -1.;
    Double_t Chi2m = 999.;
    Double_t param1[4] = {999., 999., 999., 999.};
    Double_t param2[4] = {999., 999., 999., 999.};
    Double_t params2[4] = {999., 999., 999., 999.};
    Double_t params3[4] = {999., 999., 999., 999.};
    // Silicon
    Double_t SiCoordX[4] = {-999., -999., -999., -999.};
    Double_t SiCoordY[4] = {-999., -999., -999., -999.};
    Double_t SiCoordZ[4] = {-999., -999., -999., -999.};
    Double_t SiSigmaX[4] = {-999., -999., -999., -999.};
    Double_t SiSigmaY[4] = {-999., -999., -999., -999.};
    Double_t SiMod[4] = {-999., -999., -999., -999.};
};
 
struct UpTracks
{
    Double_t param[4] = {999., 999., 999., 999.};
    Double_t CoordX[5] = {-999., -999., -999., -999., -999.};
    Double_t CoordY[5] = {-999., -999., -999., -999., -999.};
    Double_t CoordZ[5] = {-999., -999., -999., -999., -999.};
    Double_t Chi2 = 999.;
    Int_t Nhits = -1;
    Int_t Pdg = -1;
};
 
class BmnUpstreamTracking : public FairTask
{
  public:
    BmnUpstreamTracking();
    BmnUpstreamTracking(Bool_t, Int_t);
    virtual ~BmnUpstreamTracking();
 
    virtual InitStatus Init();
    virtual void Exec(Option_t* opt);
    virtual void Finish();
 
    static bool compareSegments(const Smatch& a, const Smatch& b) { return a.Chi2m < b.Chi2m; }
 
    struct MC_points
    {
        Int_t Id = -1;    //
        Int_t Pdg = -1;   //
        Int_t Np_P = 0;
        Int_t Np_Si = 0;
        Int_t Np_ch2 = 0;
        Int_t Np_ch3 = 0;
        Int_t Np = 0;
 
        Bool_t xWas2 = 0;
        Bool_t uWas2 = 0;
        Bool_t vWas2 = 0;
        Bool_t xWas3 = 0;
        Bool_t uWas3 = 0;
        Bool_t vWas3 = 0;
        Bool_t np_3si = 0;
 
        Double_t param[4] = {999., 999., 999., 999.};
        Double_t x_si[4] = {-999., -999., -999., -999.};
        Double_t y_si[4] = {-999., -999., -999., -999.};
        Double_t z_si[4] = {-999., -999., -999., -999.};
    };
    struct twotracks
    {
        Double_t track1[5] = {-999., -999., -999., -999., -999.};   // Ax, x, Ay, Y
        Int_t track1_pdg = -999;
        Double_t track2[5] = {-999., -999., -999., -999., -999.};   // Ax, x, Ay, Y
        Int_t track2_pdg = -999;
    };
 
  private:
    Bool_t fDebug = 0;
    UInt_t fEventNo = 0;   // event counter
    Int_t fRunNumber;
    Bool_t expData;
    TList fList;
    TClonesArray* fSiTracks;       //!
    TClonesArray* fSilHits;        //!
    TClonesArray* fMWPCTracks;     //!
    TClonesArray* fMWPCSegments;   //!
    TString fInputBranchName1;
    TString fInputBranchName2;
    TString fInputBranchName3;
    TString fInputBranchHits;
    //--mc--
    TClonesArray* fSiTracksSim;   //!
    TString fInputBranchNameSimTrue;
    TClonesArray* fBmnHitsArray;   //!
    TString fInputBranchName;
 
    /** Output array of Silicon Hits **/
    TClonesArray* fBmnUpstreamTracksArray;   //!
    TString fOutputTracksBranchName;
    TString fOutputFileName;
    //--------
    Double_t*** par_ab_Ch;     //!
    Double_t*** par_ab_tr;     //!
    Double_t*** SiXYhits;      //!
    Double_t*** Points;        //!
    Double_t** par_ab_SiTr;    //!
    Double_t** par_Seg_z;      //!
    Double_t** par_ab_trz;     //!
    Double_t** Xforglfit;      //!
    Double_t** Amatr;          //!
    Double_t* par_SiTr_z;      //!
    Double_t* X_shift_seg;     //!
    Double_t* X_shiftUp;       //!
    Double_t* Y_shiftUp;       //!
    Double_t* Z_pair;          //!
    Double_t* Z_Chamber;       //!
    Double_t* AmatrInverted;   //!
    Double_t* rhs;             //!
    Double_t* AmatrArray;      //!
    Double_t* line;            //!
    Int_t* Nseg_Ch;            //!
    Int_t* NPCTracks;          //!
    Int_t* NSiXYhits;          //!
    Int_t NSiTracks;
    Int_t NumPoints;
    Double_t Shift_toCenterOfMagnetX;    //  = 0.039;
    Double_t Shift_toCenterOfMagnetY;    //  = 0.13;
    Double_t Shift_toCenterOfMagnetAX;   // = 0.;
    Double_t Shift_toCenterOfMagnetAY;   // = 0.0019;
    Double_t X_shift;
    Double_t Y_shift;
 
    vector<Smatch> vtmpSeg;
    vector<Smatch> OutVector;
    vector<UpTracks> vecUpTracks;
    vector<MC_points> vec_points;
    vector<twotracks> vec_twotracks;
    //--------
    void PrepareArraysToProcessEvent();
    void ReadSiliconTracks(Double_t**, Double_t*, Int_t&, vector<MC_points>&);
    void ReadSiliconHits(Double_t***, Int_t*);
    void ReadMWPCSegments(Double_t***, Double_t**, Int_t*, vector<MC_points>&);
    void ReadMWPCTracks(Double_t***, Double_t**, Int_t*);
    void PCTracksSiTracksMatching(Double_t**,
                                  Double_t*,
                                  Int_t&,
                                  Double_t***,
                                  Double_t**,
                                  Int_t*,
                                  vector<Smatch>&,
                                  vector<Smatch>&);
    void PCSegmentsSiTracksMatching(Double_t**,
                                    Double_t*,
                                    Int_t&,
                                    Double_t***,
                                    Double_t**,
                                    Int_t*,
                                    vector<Smatch>&,
                                    vector<Smatch>&);
    void RecordingBest(Int_t&, vector<Smatch>&, vector<Smatch>&);
    void RecordCandidateSeg(Int_t&,
                            Int_t&,
                            Int_t&,
                            Double_t**,
                            Double_t***,
                            Double_t*,
                            Double_t**,
                            Double_t&,
                            Double_t&,
                            vector<Smatch>&);
    void GetAddSiHits(vector<Smatch>&, Double_t***, Int_t*);
    void GetHitsToPoints(vector<Smatch>&, Double_t***, Int_t*, Double_t***, Int_t*, Double_t***, Int_t&);
    void GlobalFit(Double_t**, Double_t**, Double_t*);
    bool InvertMatrix(Double_t*, Double_t*);
    void CalculationOfChi2(Double_t***, Int_t&, vector<UpTracks>&);
    void PrintAllTracks(vector<UpTracks>&);
    void TrackRecording(vector<UpTracks>&);
    void MCefficiencyCalculation(vector<MC_points>&, vector<UpTracks>&);
 
    //--------
    const Int_t kBig = 200;
    const Int_t kPoints = 5;        // kNumDets
    const Int_t kPoin_Par = 5;      // parameters
    const Int_t kNumPars = 10;      // parameters
    const Int_t kNumPairs = 2;      // number of MWPC pairs
    const Int_t kNumChambers = 4;   // number of MWPC stat
    const Int_t kNumStSi = 4;       // number of Si stat
    const Int_t kMaxMC = 100;
    const Int_t fNstations = 4;
    const Int_t kNPlanes = 6;
    const Double_t kChi2_Max = 50.;
    const Double_t dw = 0.25;                       // [cm]//MWPC wire step
    const Double_t SigmXmwpc = sqrt(2) * dw / 6.;   // dw/6.;//dw / sq12;
    const Double_t SigmYmwpc = sqrt(2) * dw / 3.;   // dw/3.
    // matching consts
    const Double_t cutAx = 0.02;   // 0.012;
    const Double_t cutAy = 0.02;   // 0.015;
    const Double_t cutX = 1.5;     // cm
    const Double_t cutY = 2.0;     // cm
    // important Z
    const Double_t Zcentr = -350.;         // cm
    const Double_t kZ_target = -645.191;   // cm
    const Double_t kZSi_cent = -392.5;
    // MC
    const Int_t PDG_Be7 = 1000040070;   // Be7
    const Int_t PDG_Li6 = 1000030060;   // Li6
    const Int_t PDG_Li7 = 1000030070;   // Li7
    const Int_t PDG_He3 = 1000020030;   // He3
    const Int_t PDG_He4 = 1000020040;   // He4
    const Int_t PDG_H2 = 1000010020;    // H2
    const Int_t PDG_p = 2212;           // proton
 
    // some hists
    TH1D *hAx_fitUp, *hAy_fitUp, *hx_fitUp, *hy_fitUp, *hchi2_fitUp, *hNhitsUp, *hdAx_tr_mc_comb, *hdX_tr_mc_comb,
        *hdAy_tr_mc_comb, *hdY_tr_mc_comb, *hdAx_uptr_mc, *hdX_uptr_mc, *hdAy_uptr_mc, *hdY_uptr_mc, *hDen_mcuptr,
        *hNum_mcuptr, *hEff_mcuptr, *hAx_upmc, *hAy_upmc, *hX_upmc, *hY_upmc, *hNtr_reco, *hNtr_mc, *hNrecoTrif2mc,
        *hAngle_reco, *hAngle_recoifNmc2, *hAngle_recoifNmc2Cases;                 //!
    TH2D *hNtr_mc_vs_reco, *hy_vs_x_Up, *hY_vs_Xmctrue, *hvertexXYUp, *hTyTx_Up;   //!
    vector<TH1D*> hResXst, hResYst;
 
    ClassDef(BmnUpstreamTracking, 1)
};
 
#endif