CbmKFTrErrMCPoints.cxx

Go to the documentation of this file.

/*
 *====================================================================
 *
 *  CBM KF Track Quality
 *  
 *  Authors: M.Zyzak
 *
 *  e-mail : 
 *
 *====================================================================
 *
 *  KF Fit performance 
 *
 *====================================================================
 */
#include "CbmKFTrErrMCPoints.h"
#include "CbmKF.h"
#include "CbmMCTrack.h"
#include "CbmStsAddress.h"
#include "CbmStsSetup.h"
#include <algorithm>
 
using std::vector;
 
ClassImp(CbmKFTrErrMCPoints)
 
  CbmKFTrErrMCPoints::CbmKFTrErrMCPoints()
  : StsArray(), MvdArray(), TofArray(), StsHitsArray(), MvdHitsArray() {}
 
int CbmKFTrErrMCPoints::GetNConsMCStations() {
  if ((GetNMvdPoints() + GetNStsPoints()) < 1) return 0;
  // TODO get station number of the point using methods of the point class!
  float zStation[8] = {30., 40., 50., 60., 70., 80., 90., 100.};
  vector<int> iStations;
  for (int iMvd = 0; iMvd < GetNMvdPoints(); ++iMvd) {
    iStations.push_back(GetMvdPoint(iMvd)->GetStationNr() - 1);
    //    std::cout << GetMvdPoint(iMvd)->GetStationNr() << " " << GetMvdPoint(iMvd)->GetZ() << std::endl;
  }
  for (int iSts = 0; iSts < GetNStsPoints(); ++iSts) {
    int stNumber = -1;
    for (int iSt = 0; iSt < 8; iSt++)
      if (TMath::Abs(zStation[iSt] - GetStsPoint(iSts)->GetZ()) < 2.5)
        stNumber = iSt;
    if (stNumber >= 0)
      iStations.push_back(stNumber + CbmKF::Instance()->GetNMvdStations());
  }
 
  std::sort(iStations.begin(), iStations.end());
 
  int nMaxConsStations = 1;
  int nConsStations    = 1;
  if (iStations.size() == 0) return 0;
  int iPrevSt = iStations[0];
  for (unsigned int iP = 1; iP < iStations.size(); iP++) {
    if ((iStations[iP] - iPrevSt) == 1) {
      nConsStations++;
      iPrevSt = iStations[iP];
    } else if ((iStations[iP] - iPrevSt) > 1) {
      if (nConsStations > nMaxConsStations) nMaxConsStations = nConsStations;
      nConsStations = 1;
      iPrevSt       = iStations[iP];
    }
  }
  if (nConsStations > nMaxConsStations) nMaxConsStations = nConsStations;
 
  return nMaxConsStations;
}
 
int CbmKFTrErrMCPoints::GetNConsHitStations() {
  if ((GetNMvdHits() + GetNStsHits()) < 1) return 0;
  // TODO get station number of the point using methods of the point class!
  vector<int> iStations;
  for (int iMvd = 0; iMvd < GetNMvdHits(); ++iMvd) {
    iStations.push_back(GetMvdHit(iMvd)->GetStationNr() - 1);
    //    std::cout << GetMvdHit(iMvd)->GetStationNr() << " " << GetMvdHit(iMvd)->GetZ() << std::endl;
  }
  for (int iSts = 0; iSts < GetNStsHits(); ++iSts)
    iStations.push_back(
      CbmStsSetup::Instance()->GetStationNumber(GetStsHit(iSts)->GetAddress())
      - 1 + CbmKF::Instance()->GetNMvdStations());
 
  std::sort(iStations.begin(), iStations.end());
 
  int nMaxConsStations = 1;
  int nConsStations    = 1;
  int iPrevSt          = iStations[0];
  for (unsigned int iP = 1; iP < iStations.size(); iP++) {
    if ((iStations[iP] - iPrevSt) == 1) {
      nConsStations++;
      iPrevSt = iStations[iP];
    } else if ((iStations[iP] - iPrevSt) > 1) {
      if (nConsStations > nMaxConsStations) nMaxConsStations = nConsStations;
      nConsStations = 1;
      iPrevSt       = iStations[iP];
    }
  }
  if (nConsStations > nMaxConsStations) nMaxConsStations = nConsStations;
 
  return nMaxConsStations;
}
 
int CbmKFTrErrMCPoints::GetNHitStations() {
  if ((GetNMvdHits() + GetNStsHits()) < 1) return 0;
  // TODO get station number of the point using methods of the point class!
  vector<int> iStations;
  for (int iMvd = 0; iMvd < GetNMvdHits(); ++iMvd) {
    iStations.push_back(GetMvdHit(iMvd)->GetStationNr() - 1);
    //    std::cout << GetMvdHit(iMvd)->GetStationNr() << " " << GetMvdHit(iMvd)->GetZ() << std::endl;
  }
  for (int iSts = 0; iSts < GetNStsHits(); ++iSts)
    iStations.push_back(
      CbmStsSetup::Instance()->GetStationNumber(GetStsHit(iSts)->GetAddress())
      - 1 + CbmKF::Instance()->GetNMvdStations());
 
  std::sort(iStations.begin(), iStations.end());
 
  int nStations = 1;
  int iPrevSt   = iStations[0];
  for (unsigned int iP = 1; iP < iStations.size(); iP++) {
    if ((iStations[iP] - iPrevSt) >= 1) {
      nStations++;
      iPrevSt = iStations[iP];
    }
  }
 
  return nStations;
}
 
int CbmKFTrErrMCPoints::GetNMaxMCPointsOnStation() {
  if ((GetNMvdPoints() + GetNStsPoints()) < 1) return 0;
  // TODO get station number of the point using methods of the point class!
  float zStation[8] = {30., 40., 50., 60., 70., 80., 90., 100.};
  vector<int> iStations;
  for (int iMvd = 0; iMvd < GetNMvdPoints(); ++iMvd) {
    iStations.push_back(GetMvdPoint(iMvd)->GetStationNr() - 1);
    //    std::cout << GetMvdPoint(iMvd)->GetStationNr() << " " << GetMvdPoint(iMvd)->GetZ() << std::endl;
  }
  for (int iSts = 0; iSts < GetNStsPoints(); ++iSts) {
    int stNumber = -1;
    for (int iSt = 0; iSt < 8; iSt++)
      if (TMath::Abs(zStation[iSt] - GetStsPoint(iSts)->GetZ()) < 2.5)
        stNumber = iSt;
    if (stNumber >= 0)
      iStations.push_back(stNumber + CbmKF::Instance()->GetNMvdStations());
  }
 
  std::sort(iStations.begin(), iStations.end());
 
  int nMaxMCPointsOnStation = 1;
  int nMCPointsOnStation    = 1;
  if (iStations.size() == 0) return 0;
  int iPrevSt = iStations[0];
  for (unsigned int iP = 1; iP < iStations.size(); iP++) {
    if ((iStations[iP] - iPrevSt) == 0) {
      nMCPointsOnStation++;
      iPrevSt = iStations[iP];
    } else {
      if (nMCPointsOnStation > nMaxMCPointsOnStation)
        nMaxMCPointsOnStation = nMCPointsOnStation;
      nMCPointsOnStation = 1;
      iPrevSt            = iStations[iP];
    }
  }
 
  return nMaxMCPointsOnStation;
}
 
Bool_t CbmKFTrErrMCPoints::IsReconstructable(CbmMCTrack* mcTr,
                                             int MinNStations,
                                             int PerformanceMode,
                                             float MinRecoMom) {
  Bool_t f = 1;
 
  // reject very slow tracks from analysis
  f &= (mcTr->GetP() > MinRecoMom);
  // detected at least in 4 stations
  if (PerformanceMode == 3)
    f &= (GetNConsMCStations() >= MinNStations);  // L1-MC
  if (PerformanceMode == 2)
    f &= (GetNHitStations() >= MinNStations);  // QA definition
  if (PerformanceMode == 1)
    f &= (GetNConsHitStations() >= MinNStations);  // L1 definition
 
  // maximul 4 layers for a station.
  f &= (GetNMaxMCPointsOnStation() <= 4);
 
  return f;
}