CbmL1ReadEvent.cxx

Go to the documentation of this file.

/*
 *====================================================================
 *
 *  CBM Level 1 Reconstruction 
 *  
 *  Authors: I.Kisel,  S.Gorbunov, I. Rostovtseva (2016)
 *  
 *
 *  e-mail : ikisel@kip.uni-heidelberg.de 
 *
 *====================================================================
 *
 *  Read Event information to L1
 *
 *====================================================================
 */
 
#include "CbmKF.h"
#include "CbmL1.h"
#include "CbmMatch.h"
#include "CbmStsAddress.h"
#include "CbmStsSetup.h"
#include "L1Algo/L1Algo.h"
#include "L1AlgoInputData.h"
 
#include "CbmMatch.h"
#include "CbmMuchGeoScheme.h"
#include "CbmMuchPixelHit.h"
#include "CbmMuchPoint.h"
#include "CbmTofHit.h"
#include "CbmTofPoint.h"
#include "CbmTrdHit.h"
#include "CbmTrdPoint.h"
 
//#include "CbmMvdHitMatch.h"
 
 
#include "TDatabasePDG.h"
#include "TRandom.h"
 
#include <iostream>
#include <vector>
 
using std::cout;
using std::endl;
using std::find;
using std::vector;
 
//#define MVDIDEALHITS
//#define STSIDEALHITS
 
 
static bool compareZ(const int& a, const int& b) {
  //        return (CbmL1::Instance()->vMCPoints[a].z < CbmL1::Instance()->vMCPoints[b].z);
  const CbmL1* l1 = CbmL1::Instance();
  return l1->Get_Z_vMCPoint(a) < l1->Get_Z_vMCPoint(b);
}
 
 
struct TmpHit {  // used for sort Hits before writing in the normal arrays
  int iStripF,
    iStripB;  // indices of real-strips, sts-strips (one got from detector. They consist from parts with differen positions, so will be divided before using)
  int indStripF, indStripB;  // indices of L1-strips, indices in TmpStrip arrays
  int iStation;
  int ExtIndex;  // index of hit in the TClonesArray array ( negative for MVD )
  bool isStrip;
  double u_front, u_back;  // positions of strips
  double x, y;             // position of hit
  double dx, dy, dxy;
  double du, dv;
  int iMC;  // index of MCPoint in the vMCPoints array
  double time, t_er, z;
  int Det;
  int id;
  int track;
  static bool Compare(const TmpHit& a, const TmpHit& b) {
    return (a.iStation < b.iStation)
           || ((a.iStation == b.iStation) && (a.y < b.y));
  }
};
 
struct TmpStrip {
  fscal u;
  fscal time;
  int iStation;
  int iStrip;
  bool isStrip;
  int effIndex;  // used for unefficiency
};
 
void CbmL1::ReadEvent(L1AlgoInputData* fData_, CbmEvent* event) {
  if (fVerbose >= 10) cout << "ReadEvent: start." << endl;
 
  // clear arrays for next event
  vMCPoints.clear();
  vMCPoints_in_Time_Slice.clear();
  vMCTracks.clear();
  vStsHits.clear();
  vRTracks.clear();
  vHitMCRef.clear();
  vHitStore.clear();
  dFEI2vMCPoints.clear();
  dFEI2vMCTracks.clear();
  if (fVerbose >= 10) cout << "ReadEvent: clear is done." << endl;
 
  vector<TmpHit> tmpHits;
  vector<TmpStrip> tmpStrips;
  vector<TmpStrip> tmpStripsB;
 
  // -- produce Sts hits from space points --
 
  for (int i = 0; i < NStation; i++) {
 
    fData_->StsHitsStartIndex[i] = static_cast<THitI>(-1);
    fData_->StsHitsStopIndex[i]  = 0;
  }
 
  //Int_t nMvdPoints = 0;
  nMvdPoints = 0;
  // get MVD hits
  Int_t nMvdHits  = 0;
  Int_t nMuchHits = 0;
  Int_t nTrdHits  = 0;
  Int_t nTofHits  = 0;
  // get STS hits
  int nStsHits = 0;
 
  int nStsPoints = 0;
  int nTrdPoints = 0;
  nMuchPoints    = 0;
  int nTofPoints = 0;
 
  vector<CbmLink*> ToFPointsMatch;
 
  if (fPerformance) {
    Fill_vMCTracks();
 
    for (DFSET::iterator set_it = vFileEvent.begin();
         set_it != vFileEvent.end();
         ++set_it) {
      Int_t iFile  = set_it->first;
      Int_t iEvent = set_it->second;
 
 
      if (fMvdPoints && listMvdPts) {
        Int_t nMvdPointsInEvent = fMvdPoints->Size(iFile, iEvent);
        for (Int_t iMC = 0; iMC < nMvdPointsInEvent; iMC++) {
          CbmL1MCPoint MC;
 
          MC.event = iEvent;
 
          if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 1)) {
            MC.iStation    = -1;
            L1Station* sta = algo->vStations;
            for (Int_t iSt = 0; iSt < NStsStations; iSt++)
              MC.iStation = (MC.z > sta[iSt].z[0] - 1) ? iSt : MC.iStation;
 
            Double_t dtrck          = dFEI(iFile, iEvent, MC.ID);
            DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck);
            if (trk_it == dFEI2vMCTracks.end()) continue;
            Int_t IND_Track = trk_it->second;
            vMCTracks[IND_Track].Points.push_back(vMCPoints.size());
 
            MC.ID = trk_it->second;
 
            //           vMCTracks[MC.ID].Points.push_back(vMCPoints.size());
 
            vMCPoints.push_back(MC);
            vMCPoints_in_Time_Slice.push_back(0);
 
            dFEI2vMCPoints.insert(DFEI2I::value_type(dFEI(iFile, iEvent, iMC),
                                                     vMCPoints.size() - 1));
            nMvdPoints++;
          }
        }
      }
 
      Int_t nMC = fStsPoints->Size(iFile, iEvent);
 
      for (Int_t iMC = 0; iMC < nMC; iMC++) {
        CbmL1MCPoint MC;
 
        MC.event = iEvent;
 
        if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 0)) {
          MC.iStation    = -1;
          L1Station* sta = algo->vStations + NMvdStations;
          for (Int_t iSt = 0; iSt < NStsStations; iSt++)
            MC.iStation =
              (MC.z > sta[iSt].z[0] - 2.5) ? (NMvdStations + iSt) : MC.iStation;
 
          Double_t dtrck          = dFEI(iFile, iEvent, MC.ID);
          DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck);
          if (trk_it == dFEI2vMCTracks.end()) continue;
          Int_t IND_Track = trk_it->second;
          vMCTracks[IND_Track].Points.push_back(vMCPoints.size());
 
          MC.ID = trk_it->second;
          vMCPoints.push_back(MC);
          vMCPoints_in_Time_Slice.push_back(0);
 
          dFEI2vMCPoints.insert(DFEI2I::value_type(
            dFEI(iFile, iEvent, iMC + nMvdPoints), vMCPoints.size() - 1));
          nStsPoints++;
        }
      }
 
      if (fMuchPoints) {
 
        for (Int_t iMC = 0; iMC < fMuchPoints->Size(iFile, iEvent); iMC++) {
          CbmL1MCPoint MC;
 
          MC.event = iEvent;
 
          if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 2)) {
            MC.iStation    = -1;
            L1Station* sta = algo->vStations + NMvdStations + NStsStations;
            for (Int_t iSt = 0; iSt < NMuchStations; iSt++)
              MC.iStation = (MC.z > sta[iSt].z[0] - 2.5)
                              ? (NMvdStations + NStsStations + iSt)
                              : MC.iStation;
 
            Double_t dtrck          = dFEI(iFile, iEvent, MC.ID);
            DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck);
            if (trk_it == dFEI2vMCTracks.end()) continue;
            Int_t IND_Track = trk_it->second;
 
            vMCTracks[IND_Track].Points.push_back(vMCPoints.size());
 
            MC.ID = trk_it->second;
            vMCPoints.push_back(MC);
            vMCPoints_in_Time_Slice.push_back(0);
 
            dFEI2vMCPoints.insert(DFEI2I::value_type(
              dFEI(iFile, iEvent, iMC + nMvdPoints + nStsPoints),
              vMCPoints.size() - 1));
            nMuchPoints++;
          }
        }
      }
 
 
      if (fTrdPoints)
        for (Int_t iMC = 0; iMC < fTrdPoints->Size(iFile, iEvent); iMC++) {
          CbmL1MCPoint MC;
 
          MC.event = iEvent;
 
 
          if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 3)) {
 
            MC.iStation = -1;
            L1Station* sta =
              algo->vStations + NMvdStations + NStsStations + NMuchStations;
            for (Int_t iSt = 0; iSt < NTrdStations; iSt++)
              MC.iStation =
                (MC.z > sta[iSt].z[0] - 4.0)
                  ? (NMvdStations + NStsStations + NMuchStations + iSt)
                  : MC.iStation;
 
 
            Double_t dtrck          = dFEI(iFile, iEvent, MC.ID);
            DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck);
            if (trk_it == dFEI2vMCTracks.end()) continue;
            Int_t IND_Track = trk_it->second;
 
            vMCTracks[IND_Track].Points.push_back(vMCPoints.size());
 
            MC.ID = trk_it->second;
            vMCPoints.push_back(MC);
            vMCPoints_in_Time_Slice.push_back(0);
 
            dFEI2vMCPoints.insert(DFEI2I::value_type(
              dFEI(iFile, iEvent, iMC + nMvdPoints + nStsPoints + nMuchPoints),
              vMCPoints.size() - 1));
            nTrdPoints++;
          }
        }
 
      ToFPointsMatch.resize(0);
 
      if (fTofPoints) {
        for (int j = 0; j < fTofHits->GetEntries(); j++) {
 
          CbmTofHit* mh = L1_DYNAMIC_CAST<CbmTofHit*>(fTofHits->At(j));
 
          CbmMatch* matchHitMatch =
            L1_DYNAMIC_CAST<CbmMatch*>(fTofHitDigiMatches->At(j));
 
          if (matchHitMatch->GetNofLinks() > 0) {
 
            CbmLink* link   = (CbmLink*) &matchHitMatch->GetLink(0);
            CbmTofPoint* pt = (CbmTofPoint*) fTofPoints->Get(
              link->GetFile(), link->GetEntry(), link->GetIndex());
 
            for (int iLink = 1; iLink < matchHitMatch->GetNofLinks(); iLink++) {
 
              CbmLink* link1 = (CbmLink*) &(matchHitMatch->GetLink(iLink));
 
              CbmTofPoint* pt_cur = (CbmTofPoint*) fTofPoints->Get(
                link1->GetFile(), link1->GetEntry(), link1->GetIndex());
 
              TVector3 pos_cur, pos_old, pos_hit;
 
 
              pt_cur->Position(pos_cur);
              pt->Position(pos_old);
              mh->Position(pos_hit);
 
              if (fabs(pos_cur.Z() - pos_hit.Z())
                  < fabs(pos_old.Z() - pos_hit.Z())) {
                pt   = pt_cur;
                link = link1;
              }
            }
 
            ToFPointsMatch.push_back(link);
          }  // for j
        }    // if listTrdHits
 
 
        for (UInt_t iMC = 0; iMC < ToFPointsMatch.size(); iMC++) {
          CbmL1MCPoint MC;
 
          MC.event = iEvent;
 
          int eventNr = ToFPointsMatch[iMC]->GetEntry();
 
          if (eventNr != iEvent) continue;
 
          if (!ReadMCPoint(&MC,
                           ToFPointsMatch[iMC]->GetIndex(),
                           ToFPointsMatch[iMC]->GetFile(),
                           ToFPointsMatch[iMC]->GetEntry(),
                           4)) {
 
            MC.iStation =
              NMvdStations + NStsStations + NMuchStations + NTrdStations;
 
            Double_t dtrck          = dFEI(iFile, eventNr, MC.ID);
            DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck);
            if (trk_it == dFEI2vMCTracks.end()) continue;
 
            Int_t IND_Track = trk_it->second;
 
            vMCTracks[IND_Track].Points.push_back(vMCPoints.size());
 
            MC.ID = trk_it->second;
 
            vMCPoints.push_back(MC);
            vMCPoints_in_Time_Slice.push_back(0);
 
            dFEI2vMCPoints.insert(DFEI2I::value_type(
              dFEI(iFile,
                   eventNr,
                   ToFPointsMatch[iMC]->GetIndex() + nMvdPoints + nStsPoints
                     + nMuchPoints + nTrdPoints),
              vMCPoints.size() - 1));
            nTofPoints++;
          }
        }
      }
 
    }  //time_slice
 
    for (unsigned int iTr = 0; iTr < vMCTracks.size(); iTr++) {
 
      sort(
        vMCTracks[iTr].Points.begin(), vMCTracks[iTr].Points.end(), compareZ);
 
      if (vMCTracks[iTr].mother_ID >= 0) {
        Double_t dtrck          = dFEI(vMCTracks[iTr].iFile,
                              vMCTracks[iTr].iEvent,
                              vMCTracks[iTr].mother_ID);
        DFEI2I::iterator map_it = dFEI2vMCTracks.find(dtrck);
        if (map_it == dFEI2vMCTracks.end())
          vMCTracks[iTr].mother_ID = -2;
        else
          vMCTracks[iTr].mother_ID = map_it->second;
      }
    }  //iTr
    if (fVerbose >= 10) cout << "Points in vMCTracks are sorted." << endl;
 
  }  //fPerformance
 
 
  // get MVD hits
 
  if (listMvdHits) {
 
    //     if(listMvdPts)
    //     {
    //       isUsedMvdPoint.resize(nMC);
    //       for(int iMc=0; iMc<nMC; iMc++) isUsedMvdPoint[iMc]=0;
    //     }
 
    for (int j = 0; j < listMvdHits->GetEntries(); j++) {
      TmpHit th;
      {
        CbmMvdHit* mh = L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(j));
        th.ExtIndex   = -(1 + j);
        th.iStation   = mh->GetStationNr();
        //   th.iSector  = 0;
        th.isStrip = 0;
        th.iStripF = j;
        th.iStripB = -1;
        if (th.iStripF < 0) continue;
        if (th.iStripF >= 0 && th.iStripB >= 0) th.isStrip = 1;
        if (th.iStripB < 0) th.iStripB = th.iStripF;
 
        TVector3 pos, err;
        mh->Position(pos);
        mh->PositionError(err);
 
        th.dx = mh->GetDx();
        th.dy = mh->GetDy();
 
        th.du  = mh->GetDx();
        th.dv  = mh->GetDy();
        th.dxy = mh->GetDxy();
 
        th.x = pos.X();
        th.y = pos.Y();
 
        L1Station& st = algo->vStations[th.iStation];
        th.u_front =
          th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
        th.u_back =
          th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
      }
      th.Det  = 0;
      th.iMC  = -1;
      int iMC = -1;
      //   int iMCTr = -1;
      //   if( listMvdHitMatches ){
      //     CbmMvdHitMatch *match = (CbmMvdHitMatch*) listMvdHitMatches->At(j);
      //     if( match){
      //       iMC = match->GetPointId();
      //       iMCTr = match->GetTrackId();
      //     }
      //   }
      if (fPerformance) {
        if (listMvdHitMatches) {
          CbmMatch* mvdHitMatch =
            L1_DYNAMIC_CAST<CbmMatch*>(listMvdHitMatches->At(j));
 
 
          if (mvdHitMatch->GetNofLinks() > 0)
            if (mvdHitMatch->GetLink(0).GetIndex() < nMvdPoints) {
              iMC    = mvdHitMatch->GetLink(0).GetIndex();
              th.iMC = iMC;
              //    th.track = iMC;
            }
        }
        //       if( listMvdPts && iMC>=0 ){ // TODO1: don't need this with FairLinks
        //         CbmL1MCPoint MC;
        //         if( ! ReadMCPoint( &MC, iMC, 1 ) ){
        //           MC.iStation = th.iStation;
        //       //    isUsedMvdPoint[iMC] = 1;
        //
        // //       MC.ID = iMCTr; // because atch->GetPointId() == 0 !!! and ReadMCPoint don't work
        // //       MC.z = th.iStation; // for sort in right order
        //
        //           vMCPoints.push_back( MC );
        //           nMvdPoints++;
        //           vMCPoints_in_Time_Slice.push_back(0);
        //           th.iMC = vMCPoints.size()-1;
        //         }
        //         Int_t iFile =-1;
        //         Int_t iEvent=-1;
        //
        //
        //                if ( !fTimesliceMode )
        //         {
        //           iFile  = vFileEvent.begin()->first;
        //           iEvent = vFileEvent.begin()->second;
        //         }
        //
        //       Double_t dtrck =dFEI(iFile, iEvent, iMC);
        //       DFEI2I::iterator trk_it = dFEI2vMCPoints.find(dtrck);
        //       if (trk_it==dFEI2vMCPoints.end()) continue;
        //       iMC = trk_it->second;
        //
        //       } // if listStsPts
      }
      //if(  h.MC_Point >=0 ) // DEBUG !!!!
      {
        tmpHits.push_back(th);
        nMvdHits++;
      }
    }  // for j
  }    // if listMvdHits
  if (fVerbose >= 10) cout << "ReadEvent: mvd hits are gotten." << endl;
 
 
  if (listStsHits) {
    Int_t nEnt = 0;
    if (fTimesliceMode)
      nEnt = listStsHits->GetEntries();
    else
      nEnt = (event ? event->GetNofData(ECbmDataType::kStsHit)
                    : listStsHits->GetEntries());
 
    for (Int_t j = 0; j < nEnt; j++) {
      Int_t hitIndex = 0;
      if (fTimesliceMode)
        hitIndex = j;
      else
        hitIndex = (event ? event->GetIndex(ECbmDataType::kStsHit, j) : j);
 
      CbmStsHit* sh = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndex));
 
      TmpHit th;
      {
        CbmStsHit* mh = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndex));
        th.ExtIndex   = hitIndex;
        th.Det        = 1;
        th.iStation   = NMvdStations
                      + CbmStsSetup::Instance()->GetStationNumber(
                        mh->GetAddress());  //mh->GetStationNr() - 1;
        th.isStrip = 0;
        th.iStripF = 0;  //mh->GetFrontDigiId();
        th.iStripB = 0;  //mh->GetBackDigiId();
 
        //Get time
        th.time = mh->GetTime();
        th.t_er = mh->GetTimeError();
 
 
        th.iStripF += nMvdHits;
        th.iStripB += nMvdHits;
 
        TVector3 pos, err;
        mh->Position(pos);
        mh->PositionError(err);
 
        th.x = pos.X();
        th.y = pos.Y();
 
        th.dx  = mh->GetDx();
        th.dy  = mh->GetDy();
        th.dxy = mh->GetDxy();
 
        th.du = mh->GetDu();
        th.dv = mh->GetDv();
 
        L1Station& st = algo->vStations[th.iStation];
        th.u_front =
          th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
        th.u_back =
          th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
      }
      th.iMC = -1;
 
      Int_t iMC = -1;
 
      if (fPerformance) {
        if (listStsClusterMatch) {
          const CbmMatch* frontClusterMatch = static_cast<const CbmMatch*>(
            listStsClusterMatch->At(sh->GetFrontClusterId()));
          const CbmMatch* backClusterMatch = static_cast<const CbmMatch*>(
            listStsClusterMatch->At(sh->GetBackClusterId()));
          CbmMatch stsHitMatch;
 
          for (Int_t iFrontLink = 0;
               iFrontLink < frontClusterMatch->GetNofLinks();
               iFrontLink++) {
            const CbmLink& frontLink = frontClusterMatch->GetLink(iFrontLink);
 
            for (Int_t iBackLink = 0;
                 iBackLink < backClusterMatch->GetNofLinks();
                 iBackLink++) {
              const CbmLink& backLink = backClusterMatch->GetLink(iBackLink);
              if (frontLink == backLink) {
                stsHitMatch.AddLink(frontLink);
                stsHitMatch.AddLink(backLink);
              }
            }
          }
 
          if (stsHitMatch.GetNofLinks() > 0) {
            Float_t bestWeight = 0.f;
            for (Int_t iLink = 0; iLink < stsHitMatch.GetNofLinks(); iLink++) {
              if (stsHitMatch.GetLink(iLink).GetWeight() > bestWeight) {
                bestWeight   = stsHitMatch.GetLink(iLink).GetWeight();
                Int_t iFile  = stsHitMatch.GetLink(iLink).GetFile();
                Int_t iEvent = stsHitMatch.GetLink(iLink).GetEntry();
                //                 if(!fTimesliceMode) //TODO Fix the event number in links
                //                   iEvent+=1;
                Int_t iIndex =
                  stsHitMatch.GetLink(iLink).GetIndex() + nMvdPoints;
                Double_t dtrck          = dFEI(iFile, iEvent, iIndex);
                DFEI2I::iterator trk_it = dFEI2vMCPoints.find(dtrck);
                if (trk_it == dFEI2vMCPoints.end()) continue;
                iMC = trk_it->second;
              }
            }
          }
        } else
          iMC = sh->GetRefId();  // TODO1: don't need this with FairLinks
      }                          //fPerformance
 
      if (iMC > -1) {
        th.iMC = iMC;
        // th.track = iMC;
      }
 
      tmpHits.push_back(th);
      nStsHits++;
 
    }  // for j
  }    // if listStsHits
 
  if (fMuchPixelHits) {
 
    Int_t nEnt = fMuchPixelHits->GetEntries();
 
 
    for (int j = 0; j < nEnt; j++) {
      TmpHit th;
      {
 
        CbmMuchPixelHit* mh =
          static_cast<CbmMuchPixelHit*>(fMuchPixelHits->At(j));
 
        th.ExtIndex = j;
        th.Det      = 2;
 
 
        Int_t stationNumber =
          CbmMuchGeoScheme::GetStationIndex(mh->GetAddress());
        Int_t layerNumber = CbmMuchGeoScheme::GetLayerIndex(mh->GetAddress());
 
        int DetId = stationNumber * 3 + layerNumber;
 
 
        th.iStation = DetId + NMvdStations + NStsStations;
        //Get time
        th.time = mh->GetTime() - 17;
        th.t_er = mh->GetTimeError();
 
 
        //   th.iSector  = 0;
        th.isStrip = 0;
        th.iStripF = j;
        th.iStripB = -1;
        if (th.iStripF < 0) continue;
        if (th.iStripF >= 0 && th.iStripB >= 0) th.isStrip = 1;
        if (th.iStripB < 0) th.iStripB = th.iStripF;
 
        th.iStripF += nMvdHits + nStsHits;
        th.iStripB += nMvdHits + nStsHits;
 
        th.x = mh->GetX();
        th.y = mh->GetY();
        th.z = mh->GetZ();
 
        th.dx  = mh->GetDx();
        th.dy  = mh->GetDy();
        th.dxy = 0;
 
        th.du = mh->GetDx();
        th.dv = mh->GetDy();
 
 
        L1Station& st = algo->vStations[th.iStation];
        th.u_front =
          th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
        th.u_back =
          th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
      }
      th.iMC  = -1;
      int iMC = -1;
 
 
      if (fPerformance) {
        if (listMuchHitMatches) {
          CbmMatch* matchHitMatch =
            L1_DYNAMIC_CAST<CbmMatch*>(listMuchHitMatches->At(j));
 
 
          for (Int_t iLink = 0; iLink < matchHitMatch->GetNofLinks(); iLink++) {
            if (matchHitMatch->GetLink(iLink).GetIndex() < nMuchPoints) {
              iMC          = matchHitMatch->GetLink(iLink).GetIndex();
              Int_t iIndex = iMC + nMvdPoints + nStsPoints;
 
              Int_t iFile  = matchHitMatch->GetLink(0).GetFile();
              Int_t iEvent = matchHitMatch->GetLink(0).GetEntry();
 
              Double_t dtrck          = dFEI(iFile, iEvent, iIndex);
              DFEI2I::iterator trk_it = dFEI2vMCPoints.find(dtrck);
              if (trk_it == dFEI2vMCPoints.end()) continue;
              th.iMC = trk_it->second;
              //             th.track = vMCPoints[th.iMC].ID;
              //              th.qp = vMCPoints[iMC].p;
              //    if(matchHitMatch -> GetNofLinks() == 0) continue;
              //      Float_t bestWeight = 0.f;
              //      Float_t totalWeight = 0.f;
              //     int iMCPoint = -1;
              //   CbmLink link;
 
              //              CbmMuchPoint* pt = (CbmMuchPoint*) fMuchPoints->Get(
              //                matchHitMatch->GetLink(0).GetFile(),
              //                matchHitMatch->GetLink(0).GetEntry(),
              //                matchHitMatch->GetLink(0).GetIndex());
              // double mom = sqrt(pt->GetPxOut()*pt->GetPxOut()+pt->GetPyOut()*pt->GetPyOut()+pt->GetPzOut()*pt->GetPzOut());
              //     th.p = mom;
              //     th.q = pt->GetTrackID();//(L1_DYNAMIC_CAST<CbmMCTrack*>(fMCTracks->Get(link.GetFile(),link.GetEntry(),  pt->GetTrackID()) ))->GetCharge();
 
              //                th.x = pt->GetX( th.z ) + gRandom->Gaus(0,th.dx);
              //
              //                th.y = pt->GetY(th.z) + gRandom->Gaus(0,th.dy);
              //     th.time = pt->GetTime();  //+ gRandom->Gaus(0,th.t_er);
 
              //              L1Station& st = algo->vStations[th.iStation];
              //              th.u_front =
              //                th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
              //              th.u_back =
              //                th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
            }
          }
        }
      }
 
      tmpHits.push_back(th);
      nMuchHits++;
 
    }  // for j
  }    // if listMvdHits
 
  if (listTrdHits) {
 
    for (int j = 0; j < listTrdHits->GetEntries(); j++) {
      TmpHit th;
 
      CbmTrdHit* mh = L1_DYNAMIC_CAST<CbmTrdHit*>(listTrdHits->At(j));
      th.ExtIndex   = j;
      th.Det        = 3;
      /*  num variable not used
        int num =mh->GetPlaneId();
        
        if ((mh->GetPlaneId())==0) num = 0;
        if ((mh->GetPlaneId())==1) num = 1;
        if ((mh->GetPlaneId())==2) num = 2;
        if ((mh->GetPlaneId())==3) num = 3;
        if ((mh->GetPlaneId())==4) num = 4;
//        if (num == 1) continue;
//        if (num == 3) continue;
        */
 
      //    for (int k=0; k <TrdHitsOnStation[num+1].size(); k++ ){
 
      th.iStation =
        NMvdStations + mh->GetPlaneId() + NStsStations + NMuchStations;
 
 
      th.time = mh->GetTime();
      th.t_er = mh->GetTimeError();
 
      //   th.iSector  = 0;
      th.isStrip = 0;
      th.iStripF = j;
      th.iStripB = -1;  //TrdHitsOnStationIndex[num+1][k];
 
      th.iStripF += nMvdHits + nStsHits + nMuchHits;
      th.iStripB += nMvdHits + nStsHits + nMuchHits;
 
      TVector3 pos, err;
      mh->Position(pos);
      mh->PositionError(err);
 
      th.x = pos.X();
      th.y = pos.Y();
 
 
      th.dx  = fabs(mh->GetDx());
      th.dy  = fabs(mh->GetDy());
      th.dxy = 0;
 
      th.du = fabs(mh->GetDx());
      th.dv = fabs(mh->GetDy());
 
      //         CbmTrdHit mh1 = TrdHitsOnStation[num+1][k];
      //
      //         TVector3 pos1, err1;
      //         mh1.Position(pos1);
      //         mh1.PositionError(err1);
      //
      //         if (th.dx>th.dy){
      //
      //         th.x = pos1.X();
      //
      //
      //         th.dx = mh1.GetDx();
      //
      //
      //         th.du = mh1.GetDx();
      //
      //         }
      //
      //         if (th.dy>th.dx){
      //
      //         th.y = pos1.Y();
      //
      //         th.dy = mh1.GetDy();
      //
      //         th.dv = mh1.GetDy();
      //         }
 
      L1Station& st = algo->vStations[th.iStation];
      th.u_front =
        th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
      th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
 
      th.iMC  = -1;
      int iMC = -1;
 
      if (fPerformance) {
 
        if (fTrdHitMatches) {
 
          CbmMatch* trdHitMatch =
            L1_DYNAMIC_CAST<CbmMatch*>(fTrdHitMatches->At(j));
 
          //  CbmMatch *trdHitMatch1 = L1_DYNAMIC_CAST<CbmMatch*>( fTrdHitMatches->At(TrdHitsOnStationIndex[num+1][k]) );
 
 
          if (trdHitMatch->GetNofLinks() > 0)
            if (trdHitMatch->GetLink(0).GetIndex() < nTrdPoints) {
              iMC    = trdHitMatch->GetLink(0).GetIndex();
              th.iMC = iMC + nMvdPoints + nStsPoints + nMuchPoints;
              //      th.track = vMCPoints[th.iMC].ID;
 
              //              CbmTrdPoint* pt = (CbmTrdPoint*) fTrdPoints->Get(
              //                trdHitMatch->GetLink(0).GetFile(),
              //                trdHitMatch->GetLink(0).GetEntry(),
              //                trdHitMatch->GetLink(0).GetIndex());
 
              //  float min = 0.1;
 
 
              //          if( trdHitMatch1->GetNofLinks() > 0 )
              //           if( trdHitMatch1->GetLink(0).GetIndex() < nTrdPoints )
              //           {
              //
              //             CbmTrdPoint* pt1 = (CbmTrdPoint*) fTrdPoints->Get(trdHitMatch1->GetLink(0).GetFile(),trdHitMatch1->GetLink(0).GetEntry(),trdHitMatch1->GetLink(0).GetIndex());
              //
              //            if (mh->GetDx()>mh->GetDy()){
              //
              //             th.dx = mh1.GetDx();
              //             th.du = mh1.GetDx();
              //             th.x = pt1->GetXOut()+ gRandom->Gaus(0,th.dx);
              //             }
              //
              //            if (mh->GetDy()>mh->GetDx()){
              //
              //             th.dy = mh1.GetDy();
              //             th.dv = mh1.GetDy();
              //             th.y = pt1->GetYOut()+ gRandom->Gaus(0,th.dy);
              //            }
              //
              //
              //          L1Station& st = algo->vStations[th.iStation];
              //          th.u_front =
              //            th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
              //          th.u_back =
              //            th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
              //
            }
        }
      }
      tmpHits.push_back(th);
      nTrdHits++;
 
      //}
    }  // for j
  }    // if listTrdHits
 
  if (fTofHits) {
 
    for (int j = 0; j < fTofHits->GetEntries(); j++) {
      TmpHit th;
 
      CbmTofHit* mh = L1_DYNAMIC_CAST<CbmTofHit*>(fTofHits->At(j));
      th.ExtIndex   = j;
      th.Det        = 4;
      th.iStation   = NMvdStations + mh->GetPlaneId() + NStsStations
                    + NMuchStations + NTrdStations;
 
      th.time = mh->GetTime();
 
      th.t_er = mh->GetTimeError();
 
      th.dx  = mh->GetDx();
      th.dy  = mh->GetDy();
      th.dxy = 0;
 
      th.du = mh->GetDx();
      th.dv = mh->GetDy();
 
      //   th.iSector  = 0;
      th.isStrip = 0;
      th.iStripF = j;
      th.iStripB = -1;
      if (th.iStripF < 0) continue;
      if (th.iStripF >= 0 && th.iStripB >= 0) th.isStrip = 1;
      if (th.iStripB < 0) th.iStripB = th.iStripF;
 
      th.iStripF += nMvdHits + nStsHits + nMuchHits + nTrdHits;
      th.iStripB += nMvdHits + nStsHits + nMuchHits + nTrdHits;
 
      TVector3 pos, err;
      mh->Position(pos);
      mh->PositionError(err);
 
      th.x = pos.X();
      th.y = pos.Y();
 
      L1Station& st = algo->vStations[th.iStation];
      th.u_front =
        th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
      th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
 
 
      th.iMC = -1;
      //      int iMC = -1;
 
 
      if (fPerformance) {
 
        //   if (ToFPointsMatch[j]!=0)
 
        //  th.iMC = j+nMvdPoints+nStsPoints+nTrdPoints+nMuchPoints;
 
        Int_t iFile  = ToFPointsMatch[j]->GetFile();
        Int_t iEvent = ToFPointsMatch[j]->GetEntry();
 
        Int_t iIndex = ToFPointsMatch[j]->GetIndex() + nMvdPoints + nStsPoints
                       + nMuchPoints + nTrdPoints;
 
        Double_t dtrck          = dFEI(iFile, iEvent, iIndex);
        DFEI2I::iterator trk_it = dFEI2vMCPoints.find(dtrck);
        if (trk_it == dFEI2vMCPoints.end()) continue;
        th.iMC = trk_it->second;
        //         th.track = vMCPoints[th.iMC].ID;
        //         th.qp = vMCPoints[th.iMC].p;
 
        //        CbmTofPoint* pt = (CbmTofPoint*) fTofPoints->Get(ToFPointsMatch[j]->GetFile(),ToFPointsMatch[j]->GetEntry(),ToFPointsMatch[j]->GetIndex());
 
 
        //                 th.x = pt->GetX() + gRandom->Gaus(0,th.dx);
        //
        //         th.y = pt->GetY()+ gRandom->Gaus(0,th.dy);
        //         th.time =  pt->GetTime()+ gRandom->Gaus(0,th.t_er);
        //
        //                 L1Station &st = algo->vStations[th.iStation];
        //         th.u_front = th.x*st.frontInfo.cos_phi[0] + th.y*st.frontInfo.sin_phi[0];
        //         th.u_back  = th.x* st.backInfo.cos_phi[0] + th.y*st.backInfo.sin_phi[0];
      }
 
 
      tmpHits.push_back(th);
      nTofHits++;
 
    }  // for j
  }    // if listTrdHits
 
  if (fVerbose >= 10) cout << "ReadEvent: sts hits are gotten." << endl;
 
  // sort hits
  int nHits = nMvdHits + nStsHits + nMuchHits + nTrdHits + nTofHits;
 
  sort(tmpHits.begin(), tmpHits.end(), TmpHit::Compare);
 
  // -- create strips --
  int NStrips = 0, NStripsB = 0;
  for (int ih = 0; ih < nHits; ih++) {
    TmpHit& th = tmpHits[ih];
 
    // try to find the respective front and back strip from the already created
    th.indStripF = -1;
    th.indStripB = -1;
    for (int is = 0; is < NStrips; is++) {
      TmpStrip& s = tmpStrips[is];
      if (s.iStation != th.iStation) continue;
      if (s.iStrip != th.iStripF) continue;
      if (fabs(s.time - th.time) > 30) continue;
      if (fabs(s.u - th.u_front) > 1.e-4) continue;
      th.indStripF = is;
    }
    for (int is = 0; is < NStripsB; is++) {
      TmpStrip& s = tmpStripsB[is];
      if (s.iStation != th.iStation) continue;
      if (s.iStrip != th.iStripB) continue;
      if (fabs(s.time - th.time) > 30) continue;
      if (fabs(s.u - th.u_back) > 1.e-4) continue;
      th.indStripB = is;
    }
    // create new strips
    if (th.indStripF < 0) {
      TmpStrip tmp;
      tmp.iStation = th.iStation;
      tmp.iStrip   = th.iStripF;
      tmp.u        = th.u_front;
      tmp.time     = th.time;
      tmp.isStrip  = th.isStrip;
      tmpStrips.push_back(tmp);
      th.indStripF = NStrips++;
    }
    if (th.indStripB < 0) {
      TmpStrip tmp1;
      tmp1.iStation = th.iStation;
      tmp1.iStrip   = th.iStripB;
      tmp1.isStrip  = th.isStrip;
      tmp1.u        = th.u_back;
      tmp1.time     = th.time;
      tmpStripsB.push_back(tmp1);
      th.indStripB = NStripsB++;
    }
  }  // ih
 
  // take into account unefficiency and save strips in L1Algo
  Int_t NEffStrips = 0, NEffStripsB = 0;
  for (int i = 0; i < NStrips; i++) {
    TmpStrip& ts = tmpStrips[i];
    //     if( ts.effIndex == -1 ){
    ts.effIndex = NEffStrips++;
    char flag   = ts.iStation * 4;
 
    fData_->vStsStrips.push_back(ts.u);
    fData_->vSFlag.push_back(flag);
  }
  for (int i = 0; i < NStripsB; i++) {
    TmpStrip& ts = tmpStripsB[i];
    //     if( ts.effIndex == -1 ){
    ts.effIndex = NEffStripsB++;
    char flag   = ts.iStation * 4;
 
    fData_->vStsStripsB.push_back(ts.u);
    fData_->vSFlagB.push_back(flag);
  }
 
 
  if (fVerbose >= 10) cout << "ReadEvent: strips are read." << endl;
 
  // -- save hits --
  int nEffHits = 0;
  vector<float>
    vStsZPos_temp;  // temp array for unsorted z positions of detectors segments
  for (int i = 0; i < nHits; i++) {
    TmpHit& th = tmpHits[i];
 
    CbmL1HitStore s;
    s.Det      = th.Det;
    s.ExtIndex = th.ExtIndex;
    s.iStation = th.iStation;
    s.x        = th.x;
    s.y        = th.y;
    s.dx       = th.dx;
    s.dy       = th.dy;
    s.dxy      = th.dxy;
    s.time     = th.time;
 
    if (th.indStripF < 0 || th.indStripF >= NStrips) continue;
    if (th.indStripB < 0 || th.indStripB >= NStripsB) continue;
 
    TmpStrip& tsF = tmpStrips[th.indStripF];
    TmpStrip& tsB = tmpStripsB[th.indStripB];
 
    L1StsHit h;
    h.f = tsF.effIndex;
    h.b = tsB.effIndex;
 
    h.t_reco = th.time;
    h.t_er   = th.t_er;
    //  h.track = th.track;
    //    h.dx  = th.dx;
    //    h.dy  = th.dy;
    h.du = th.du;
    h.dv = th.dv;
    //    h.dxy = th.dxy;
    //     h.p = th.p;
    //     h.q = th.q;
    // h.ista = th.iStation;
 
    // find and save z positions
    float z_tmp = -111.;
    int ist     = th.iStation;
    if (ist < NMvdStations) {
#ifdef MVDIDEALHITS
      CbmMvdPoint* point =
        L1_DYNAMIC_CAST<CbmMvdPoint*>(listMvdPts->At(-s.ExtIndex - 1));
      z_tmp = 0.5 * (point->GetZOut() + point->GetZ());
#else
      CbmMvdHit* mh_m =
        L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(-s.ExtIndex - 1));
      z_tmp = mh_m->GetZ();
#endif
    }
 
    if ((ist >= NMvdStations) && (ist < (NStsStations + NMvdStations))) {
#ifdef STSIDEALHITS
      CbmStsPoint* point =
        L1_DYNAMIC_CAST<CbmStsPoint*>(listStsPts->At(s.ExtIndex));
      z_tmp = 0.5 * (point->GetZOut() + point->GetZIn());
#else
      CbmStsHit* mh_m =
        L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(s.ExtIndex));
      z_tmp = mh_m->GetZ();
#endif
    }
 
    if ((ist >= NStsStations + NMvdStations)
        && (ist < (NStsStations + NMvdStations + NMuchStations))) {
      //#ifdef STSIDEALHITS
      //  CbmStsPoint* point = L1_DYNAMIC_CAST<CbmStsPoint*>(listStsPts->At(s.ExtIndex));
      //   z_tmp = 0.5 * ( point->GetZOut() + point->GetZIn() );
      //#else
      CbmMuchPixelHit* mh =
        static_cast<CbmMuchPixelHit*>(fMuchPixelHits->At(s.ExtIndex));
      z_tmp = mh->GetZ();
      //#endif
    }
 
    if ((ist >= NStsStations + NMvdStations + NMuchStations)
        && (ist
            < (NStsStations + NMvdStations + NMuchStations + NTrdStations))) {
      //#ifdef STSIDEALHITS
      //  CbmStsPoint* point = L1_DYNAMIC_CAST<CbmStsPoint*>(listStsPts->At(s.ExtIndex));
      //   z_tmp = 0.5 * ( point->GetZOut() + point->GetZIn() );
      //#else
      CbmTrdHit* mh = static_cast<CbmTrdHit*>(listTrdHits->At(s.ExtIndex));
      z_tmp         = mh->GetZ();
      //#endif
    }
 
    if ((ist >= NStsStations + NMvdStations + NMuchStations + NTrdStations
                  + NTOFStation)
        && (ist
            < (NStsStations + NMvdStations + NMuchStations + NTrdStations))) {
      //#ifdef STSIDEALHITS
      //  CbmStsPoint* point = L1_DYNAMIC_CAST<CbmStsPoint*>(listStsPts->At(s.ExtIndex));
      //   z_tmp = 0.5 * ( point->GetZOut() + point->GetZIn() );
      //#else
      CbmTofHit* mh = static_cast<CbmTofHit*>(fTofHits->At(s.ExtIndex));
      z_tmp         = mh->GetZ();
      //#endif
    }
 
    unsigned int k;
    for (k = 0; k < vStsZPos_temp.size(); k++) {
      if (vStsZPos_temp[k] == z_tmp) {
        h.iz = k;
        break;
      }
    }
    if (k == vStsZPos_temp.size()) {
      h.iz = vStsZPos_temp.size();
      vStsZPos_temp.push_back(z_tmp);
    }
 
    // save hit
    vStsHits.push_back(
      CbmL1StsHit(fData->vStsHits.size(), th.ExtIndex, th.Det));
 
    vStsHits[vStsHits.size() - 1].x = th.x;
    vStsHits[vStsHits.size() - 1].y = th.y;
    vStsHits[vStsHits.size() - 1].t = th.time;
 
 
    fData_->vStsHits.push_back(h);
 
    int sta = th.iStation;
 
    if (fData_->StsHitsStartIndex[sta] == static_cast<THitI>(-1))
      fData_->StsHitsStartIndex[sta] = nEffHits;
    nEffHits++;
 
    fData_->StsHitsStopIndex[sta] = nEffHits;
 
    vHitStore.push_back(s);
    vHitMCRef.push_back(th.iMC);
  }
 
 
  for (int i = 0; i < NStation; i++) {
 
    if (fData_->StsHitsStartIndex[i] == static_cast<THitI>(-1))
      fData_->StsHitsStartIndex[i] = fData_->StsHitsStopIndex[i];
  }
 
  if (fVerbose >= 10) cout << "ReadEvent: mvd and sts are saved." << endl;
 
 
  // sort z-pos
  if (vStsZPos_temp.size() != 0) {
    vector<float> vStsZPos_temp2;
    vStsZPos_temp2.clear();
    vStsZPos_temp2.push_back(vStsZPos_temp[0]);
    vector<int> newToOldIndex;
    newToOldIndex.clear();
    newToOldIndex.push_back(0);
 
    for (unsigned int k = 1; k < vStsZPos_temp.size(); k++) {
      vector<float>::iterator itpos = vStsZPos_temp2.begin() + 1;
      vector<int>::iterator iti     = newToOldIndex.begin() + 1;
      for (; itpos < vStsZPos_temp2.end(); itpos++, iti++) {
        if (vStsZPos_temp[k] < *itpos) {
          vStsZPos_temp2.insert(itpos, vStsZPos_temp[k]);
          newToOldIndex.insert(iti, k);
          break;
        }
      }
      if (itpos == vStsZPos_temp2.end()) {
        vStsZPos_temp2.push_back(vStsZPos_temp[k]);
        newToOldIndex.push_back(k);
      }
    }  // k
 
 
    if (fVerbose >= 10) cout << "ReadEvent: z-pos are sorted." << endl;
 
    for (unsigned int k = 0; k < vStsZPos_temp2.size(); k++)
      fData_->vStsZPos.push_back(vStsZPos_temp2[k]);
 
    int size_nto_tmp = newToOldIndex.size();
    vector<int> oldToNewIndex;
    oldToNewIndex.clear();
    oldToNewIndex.resize(size_nto_tmp);
    for (int k = 0; k < size_nto_tmp; k++)
      oldToNewIndex[newToOldIndex[k]] = k;
 
    int size_hs_tmp = vHitStore.size();
    for (int k = 0; k < size_hs_tmp; k++)
      fData_->vStsHits[k].iz = oldToNewIndex[fData_->vStsHits[k].iz];
  }
 
  if (fVerbose >= 10) cout << "ReadEvent: z-pos are saved." << endl;
 
  algo->SetData(fData_->GetStsHits(),
                fData_->GetStsStrips(),
                fData_->GetStsStripsB(),
                fData_->GetStsZPos(),
                fData_->GetSFlag(),
                fData_->GetSFlagB(),
                fData_->GetStsHitsStartIndex(),
                fData_->GetStsHitsStopIndex());
 
  if (fPerformance) {
 
    if (fVerbose >= 10) cout << "HitMatch is done." << endl;
    if (fVerbose >= 10) cout << "MCPoints and MCTracks are saved." << endl;
  }
 
  if (fVerbose >= 10) cout << "ReadEvent is done." << endl;
 
 
}  // void CbmL1::ReadEvent()
 
 
void CbmL1::Fill_vMCTracks() {
  PrimVtx.MC_ID = 999;
  {
    CbmL1Vtx Vtxcurr;
    int nvtracks = 0, nvtrackscurr = 0;
 
    vMCTracks.clear();
 
    for (DFSET::iterator set_it = vFileEvent.begin();
         set_it != vFileEvent.end();
         ++set_it) {
      Int_t iFile  = set_it->first;
      Int_t iEvent = set_it->second;
 
 
      Int_t nMCTrack = fMCTracks->Size(iFile, iEvent);
 
 
      for (Int_t iMCTrack = 0; iMCTrack < nMCTrack; iMCTrack++) {
        CbmMCTrack* MCTrack =
          L1_DYNAMIC_CAST<CbmMCTrack*>(fMCTracks->Get(iFile, iEvent, iMCTrack));
        if (!MCTrack) continue;
 
        int mother_ID = MCTrack->GetMotherId();
 
 
        if (mother_ID < 0 && mother_ID != -2) mother_ID = -iEvent - 1;
        TVector3 vr;
        TLorentzVector vp;
        MCTrack->GetStartVertex(vr);
        MCTrack->Get4Momentum(vp);
 
        Int_t pdg  = MCTrack->GetPdgCode();
        Double_t q = 1, mass = 0.;
        if (pdg < 9999999
            && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) {
          q    = TDatabasePDG::Instance()->GetParticle(pdg)->Charge() / 3.0;
          mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
        } else
          q = 0;
        Int_t IND_Track = vMCTracks.size();  //or iMCTrack?
        CbmL1MCTrack T(mass, q, vr, vp, IND_Track, mother_ID, pdg);
        //        CbmL1MCTrack T(mass, q, vr, vp, iMCTrack, mother_ID, pdg);
        T.time   = MCTrack->GetStartT();
        T.iFile  = iFile;
        T.iEvent = iEvent;
 
        vMCTracks.push_back(T);
        //    Double_t dtrck =dFEI(iFile,iEvent,iMCTrack);
        dFEI2vMCTracks.insert(DFEI2I::value_type(dFEI(iFile, iEvent, iMCTrack),
                                                 vMCTracks.size() - 1));
 
        if (T.mother_ID < 0) {  // vertex track
          if (PrimVtx.MC_ID == 999
              || fabs(T.z - Vtxcurr.MC_z) > 1.e-7) {  // new vertex
            if (nvtrackscurr > nvtracks) {
              PrimVtx  = Vtxcurr;
              nvtracks = nvtrackscurr;
            }
            Vtxcurr.MC_x  = T.x;
            Vtxcurr.MC_y  = T.y;
            Vtxcurr.MC_z  = T.z;
            Vtxcurr.MC_ID = T.mother_ID;
            nvtrackscurr  = 1;
          } else
            nvtrackscurr++;
        }
      }  //iTracks
    }    //Links
    if (nvtrackscurr > nvtracks) PrimVtx = Vtxcurr;
  }  //PrimVtx
 
  if (fVerbose && PrimVtx.MC_ID == 999) cout << "No primary vertex !!!" << endl;
}  //Fill_vMCTracks
 
bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC,
                        int iPoint,
                        int file,
                        int event,
                        int MVD) {
  TVector3 xyzI, PI, xyzO, PO;
  Int_t mcID    = -1;
  Double_t time = 0.f;
  if (MVD == 1) {
    CbmMvdPoint* pt = L1_DYNAMIC_CAST<CbmMvdPoint*>(
      fMvdPoints->Get(file, event, iPoint));  // file, event, object
    //CbmMvdPoint *pt = L1_DYNAMIC_CAST<CbmMvdPoint*> (Point);
 
    if (!pt) return 1;
    pt->Position(xyzI);
    pt->Momentum(PI);
    pt->PositionOut(xyzO);
    pt->MomentumOut(PO);
    mcID = pt->GetTrackID();
    time = pt->GetTime();
  }
  if (MVD == 0) {
    CbmStsPoint* pt = L1_DYNAMIC_CAST<CbmStsPoint*>(
      fStsPoints->Get(file, event, iPoint));  // file, event, object
    if (!pt) return 1;
    //     if ( fTimesliceMode )
    //     {
    //       Double_t StartTime = fTimeSlice->GetStartTime();
    //       Double_t EndTime = fTimeSlice->GetEndTime();
    //       Double_t Time_MC_point =  pt->GetTime() + fEventList->GetEventTime(event, file);
    //       if (Time_MC_point < StartTime )
    //         return 1;
    //
    //       if (Time_MC_point > EndTime )
    //         return 1;
    //     } //if ( fTimesliceMode )
 
    pt->Position(xyzI);
    pt->Momentum(PI);
    pt->PositionOut(xyzO);
    pt->MomentumOut(PO);
    mcID = pt->GetTrackID();
    time = pt->GetTime();
  }
 
 
  if (MVD == 2) {
    CbmMuchPoint* pt = L1_DYNAMIC_CAST<CbmMuchPoint*>(
      fMuchPoints->Get(file, event, iPoint));  // file, event, object
    if (!pt) return 1;
    if (fTimesliceMode) {
      Double_t StartTime = fTimeSlice->GetStartTime();
      Double_t EndTime   = fTimeSlice->GetEndTime();
      Double_t Time_MC_point =
        pt->GetTime() + fEventList->GetEventTime(event, file);
      if (Time_MC_point < StartTime) return 1;
 
      if (Time_MC_point > EndTime) return 1;
    }  //if ( fTimesliceMode )
 
    pt->Position(xyzI);
    pt->Momentum(PI);
    pt->PositionOut(xyzO);
    pt->Momentum(PO);
    mcID = pt->GetTrackID();
    time = pt->GetTime();
  }
 
 
  if (MVD == 3) {
    CbmTrdPoint* pt = L1_DYNAMIC_CAST<CbmTrdPoint*>(
      fTrdPoints->Get(file, event, iPoint));  // file, event, object
 
    if (!pt) return 1;
    if (fTimesliceMode) {
      //Double_t StartTime = fTimeSlice->GetStartTime();  // not used
      //Double_t EndTime = fTimeSlice->GetEndTime();      // not used
      //Double_t Time_MC_point =  pt->GetTime() + fEventList->GetEventTime(event, file); // not used
      //       if (Time_MC_point < StartTime )
      //         return 1;
      //
      //       if (Time_MC_point > EndTime )
      //         return 1;
    }  //if ( fTimesliceMode )
 
    pt->Position(xyzI);
    pt->Momentum(PI);
    pt->PositionOut(xyzO);
    pt->MomentumOut(PO);
    mcID = pt->GetTrackID();
 
    time = pt->GetTime();
  }
 
  if (MVD == 4) {
    CbmTofPoint* pt = L1_DYNAMIC_CAST<CbmTofPoint*>(
      fTofPoints->Get(file, event, iPoint));  // file, event, object
    if (!pt) return 1;
    if (fTimesliceMode) {
      Double_t StartTime = fTimeSlice->GetStartTime();
      Double_t EndTime   = fTimeSlice->GetEndTime();
      Double_t Time_MC_point =
        pt->GetTime() + fEventList->GetEventTime(event, file);
      if (Time_MC_point < StartTime) return 1;
 
      if (Time_MC_point > EndTime) return 1;
    }  //if ( fTimesliceMode )
 
    pt->Position(xyzI);
    pt->Momentum(PI);
    pt->Position(xyzO);
    pt->Momentum(PO);
    mcID = pt->GetTrackID();
    time = pt->GetTime();
  }
 
  TVector3 xyz = .5 * (xyzI + xyzO);
  TVector3 P   = .5 * (PI + PO);
  MC->x        = xyz.X();
  MC->y        = xyz.Y();
  MC->z        = xyz.Z();
  MC->px       = P.X();
  MC->py       = P.Y();
  MC->pz       = P.Z();
  MC->xIn      = xyzI.X();
  MC->yIn      = xyzI.Y();
  MC->zIn      = xyzI.Z();
  MC->pxIn     = PI.X();
  MC->pyIn     = PI.Y();
  MC->pzIn     = PI.Z();
  MC->xOut     = xyzO.X();
  MC->yOut     = xyzO.Y();
  MC->zOut     = xyzO.Z();
  MC->pxOut    = PO.X();
  MC->pyOut    = PO.Y();
  MC->pzOut    = PO.Z();
  MC->p       = sqrt(fabs(MC->px * MC->px + MC->py * MC->py + MC->pz * MC->pz));
  MC->ID      = mcID;
  MC->pointId = iPoint;
  MC->file    = file;
  MC->event   = event;
 
  MC->time = time;
 
  if (MC->ID < 0) return 1;
  CbmMCTrack* MCTrack =
    L1_DYNAMIC_CAST<CbmMCTrack*>(fMCTracks->Get(file, event, MC->ID));
  if (!MCTrack) return 1;
  MC->pdg       = MCTrack->GetPdgCode();
  MC->mother_ID = MCTrack->GetMotherId();
 
  TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(MC->pdg);
  if (particlePDG) {
    MC->q    = particlePDG->Charge() / 3.0;
    MC->mass = particlePDG->Mass();
  }
 
  return 0;
}
 
bool CbmL1::ReadMCPoint(CbmL1MCPoint* /*MC*/, int /*iPoint*/, int /*MVD*/) {
  return 0;
}
 
void CbmL1::HitMatch() {
  const int NHits = vStsHits.size();
  for (int iH = 0; iH < NHits; iH++) {
    CbmL1StsHit& hit = vStsHits[iH];
 
    if (hit.Det == 1) {
      CbmStsHit* sh =
        L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(vStsHits[iH].extIndex));
 
      int iP = -1;
 
      vector<int> iEvent1;
 
 
      if (listStsClusterMatch) {
 
        const CbmMatch* frontClusterMatch = static_cast<const CbmMatch*>(
          listStsClusterMatch->At(sh->GetFrontClusterId()));
        const CbmMatch* backClusterMatch = static_cast<const CbmMatch*>(
          listStsClusterMatch->At(sh->GetBackClusterId()));
        CbmMatch stsHitMatch;
 
        Float_t Sum_of_front = 0;
        Float_t Sum_of_back  = 0;
 
 
        for (Int_t iFrontLink = 0;
             iFrontLink < frontClusterMatch->GetNofLinks();
             iFrontLink++) {
          const CbmLink& frontLink = frontClusterMatch->GetLink(iFrontLink);
          Sum_of_front             = Sum_of_front + frontLink.GetWeight();
        }
 
        for (Int_t iBackLink = 0; iBackLink < backClusterMatch->GetNofLinks();
             iBackLink++) {
          const CbmLink& backLink = backClusterMatch->GetLink(iBackLink);
          Sum_of_back             = Sum_of_back + backLink.GetWeight();
        }
 
        for (Int_t iFrontLink = 0;
             iFrontLink < frontClusterMatch->GetNofLinks();
             iFrontLink++) {
          const CbmLink& frontLink = frontClusterMatch->GetLink(iFrontLink);
 
          // Float_t Fraction_front = frontLink.GetWeight()/Sum_of_front;
 
          for (Int_t iBackLink = 0; iBackLink < backClusterMatch->GetNofLinks();
               iBackLink++) {
            const CbmLink& backLink = backClusterMatch->GetLink(iBackLink);
 
            // Float_t  Fraction_back = backLink.GetWeight()/Sum_of_back;
 
            if (frontLink == backLink) {
              stsHitMatch.AddLink(frontLink);
              stsHitMatch.AddLink(backLink);
            }
          }
        }
 
        Float_t bestWeight  = 0.f;
        Float_t totalWeight = 0.f;
        for (Int_t iLink = 0; iLink < stsHitMatch.GetNofLinks(); iLink++) {
          Int_t iFile  = stsHitMatch.GetLink(iLink).GetFile();
          Int_t iEvent = stsHitMatch.GetLink(iLink).GetEntry();
          Int_t iIndex = stsHitMatch.GetLink(iLink).GetIndex();
 
          iEvent1.push_back(iEvent);
 
 
          if (!fTimesliceMode) {
            iFile  = vFileEvent.begin()->first;
            iEvent = vFileEvent.begin()->second;
          }
          int nMvdPoints_ = 0;
 
          if (!fTimesliceMode)
            if (listMvdPts) nMvdPoints_ = listMvdPts->GetEntriesFast();
 
          Double_t dtrck          = dFEI(iFile, iEvent, iIndex + nMvdPoints_);
          DFEI2I::iterator trk_it = dFEI2vMCPoints.find(dtrck);
 
          if (trk_it == dFEI2vMCPoints.end()) continue;
 
          totalWeight += stsHitMatch.GetLink(iLink).GetWeight();
          if (stsHitMatch.GetLink(iLink).GetWeight() > bestWeight) {
            bestWeight = stsHitMatch.GetLink(iLink).GetWeight();
            iP         = trk_it->second;
          }
        }
      }  //mach cluster
 
 
      if (iP >= 0) {
 
        for (unsigned int i = 0; i < iEvent1.size(); i++) {
          hit.event = iEvent1[i];
        }
 
        hit.event = vMCPoints[iP].event;
        hit.mcPointIds.push_back(iP);
        vMCPoints[iP].hitIds.push_back(iH);
      } else {
        hit.event   = -1;
        int idPoint = vHitMCRef[iH];
        if (idPoint >= 0) {
          hit.mcPointIds.push_back(idPoint);
          vMCPoints[idPoint].hitIds.push_back(iH);
        }
      }  // if no clusters
    }
 
 
    if (hit.Det == 0) {  // if no use Links or this is mvd hit
 
      int iP = vHitMCRef[iH];
      if (iP >= 0) {
        hit.mcPointIds.push_back(iP);
        vMCPoints[iP].hitIds.push_back(iH);
      }
    }
 
    if (hit.Det == 2) {  // if no use Links or this is mvd hit
 
      int iP = vHitMCRef[iH];
 
      if (iP >= 0) {
        hit.mcPointIds.push_back(iP);
        vMCPoints[iP].hitIds.push_back(iH);
      }
    }
 
    if (hit.Det == 3) {  // if no use Links or this is mvd hit
 
 
      int iP = vHitMCRef[iH];
      if (iP >= 0) {
        hit.mcPointIds.push_back(iP);
        vMCPoints[iP].hitIds.push_back(iH);
      }
    }
 
    if (hit.Det == 4) {  // if no use Links or this is mvd hit
 
 
      int iP = vHitMCRef[iH];
 
      if (iP >= 0) {
        hit.mcPointIds.push_back(iP);
        vMCPoints[iP].hitIds.push_back(iH);
      }
    }
 
 
  }  // for hits
}