CbmStsEfficiency.cxx

Go to the documentation of this file.

/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Dario Ramirez [committer] */
 
#include "CbmStsEfficiency.h"
 
#include "CbmStsAddress.h"
#include "CbmStsUtils.h"
 
CbmStsEfficiency::CbmStsEfficiency(uint32_t layer_idx) : fTestLayer(layer_idx)
{
  LOG(debug) << "Creating an instance of CbmStsEfficiency ...";
}
CbmStsEfficiency::CbmStsEfficiency(uint32_t test_layer, double max_dis_trg_vtx, double max_dca_trk_vtx,
                                   double default_residual)
  : fTestLayer(test_layer)
  , fMaxDisTrgVtx(max_dis_trg_vtx)
  , fMaxDCATrkVtx(max_dca_trk_vtx)
  , fDefaultResidual(default_residual)
{
  LOG(debug) << "Creating an instance of CbmStsEfficiency ...";
}
 
void CbmStsEfficiency::BookHistograms()
{
  // ---- Binning configuration ----
  // Vertex
  double vertex_x_min = -15;
  double vertex_x_max = +15;
  double vertex_y_min = -15;
  double vertex_y_max = +15;
  double vertex_z_min = -15;
  double vertex_z_max = +15;
  double vertex_dx    = 0.1;
  double vertex_dy    = 0.1;
  double vertex_dz    = 0.1;
 
  auto x_vtx_binning =
    cbm_sts_utils::HBinning{uint32_t((vertex_x_max - vertex_x_min) / vertex_dx), vertex_x_min, vertex_x_max};
  auto y_vtx_binning =
    cbm_sts_utils::HBinning{uint32_t((vertex_y_max - vertex_y_min) / vertex_dy), vertex_y_min, vertex_y_max};
  auto z_vtx_binning =
    cbm_sts_utils::HBinning{uint32_t((vertex_z_max - vertex_z_min) / vertex_dz), vertex_z_min, vertex_z_max};
 
  // DCA
  double dca_dx      = .001;  // 10 um
  double dca_min     = -2.5 - 0.5 * dca_dx;
  double dca_max     = +2.5 + 0.5 * dca_dx;
  auto r_dca_binning = cbm_sts_utils::HBinning{uint32_t((dca_max - dca_min) / dca_dx), dca_min, dca_max};
 
  // Residuals track - hit
  double res_dx      = +0.001;  // 10 um
  double res_min     = -2.50 - 0.5 * res_dx;
  double res_max     = +2.50 + 0.5 * res_dx;
  auto r_sen_binning = cbm_sts_utils::HBinning{uint32_t((res_max - res_min) / res_dx), res_min, res_max};
  // ---- --------------------- ----
 
  std::string h_name;
  for (auto& [element, geo] : fStsGeoInfo) {
    std::string h_name_base = "";
    if (element == static_cast<int32_t>(fTestLayer)) {
      h_name_base = Form("Sts%d", element);
    }
    else if (element > 8) {
      if (CbmStsAddress::GetElementId(element, kStsUnit) != fTestLayer) continue;
      fDUTList.push_back(element);
      h_name_base = Form("Sts0x%x", element);
    }
 
    if (!h_name_base.length()) continue;
 
    auto x_sen_binning = cbm_sts_utils::HBinning{uint32_t((geo[1] - geo[0]) / cbm_sts_utils::kStsDx), geo[0], geo[1]};
    auto y_sen_binning = cbm_sts_utils::HBinning{uint32_t((geo[3] - geo[2]) / cbm_sts_utils::kStsDy), geo[2], geo[3]};
 
    // ---- efficiency ----
    h_name = Form("%s_found", h_name_base.c_str());
    fH2D[h_name] =
      std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), x_sen_binning.n_of_bins, x_sen_binning.x_min,
                             x_sen_binning.x_max, y_sen_binning.n_of_bins, y_sen_binning.x_min, y_sen_binning.x_max);
    fH2D[h_name]->GetXaxis()->SetTitle("StsHit^{rec}_{X} [cm]");
    fH2D[h_name]->GetYaxis()->SetTitle("StsHit^{rec}_{Y} [cm]");
 
    h_name = Form("%s_track", h_name_base.c_str());
    fH2D[h_name] =
      std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), x_sen_binning.n_of_bins, x_sen_binning.x_min,
                             x_sen_binning.x_max, y_sen_binning.n_of_bins, y_sen_binning.x_min, y_sen_binning.x_max);
    fH2D[h_name]->GetXaxis()->SetTitle("StsHit^{ext}_{X} [cm]");
    fH2D[h_name]->GetYaxis()->SetTitle("StsHit^{ext}_{Y} [cm]");
 
    h_name = Form("%s_hit_map", h_name_base.c_str());
    fH2D[h_name] =
      std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), x_sen_binning.n_of_bins, x_sen_binning.x_min,
                             x_sen_binning.x_max, y_sen_binning.n_of_bins, y_sen_binning.x_min, y_sen_binning.x_max);
    fH2D[h_name]->GetXaxis()->SetTitle("StsHit_{X} [cm]");
    fH2D[h_name]->GetYaxis()->SetTitle("StsHit_{Y} [cm]");
    // ---- ---------- ----
 
    // ---- residuals ----
    for (const char* di : {"x", "y"}) {
      h_name = Form("%s_d%s_vs_x", h_name_base.c_str(), di);
      fH2D[h_name] =
        std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), r_sen_binning.n_of_bins, r_sen_binning.x_min,
                               r_sen_binning.x_max, x_sen_binning.n_of_bins, x_sen_binning.x_min, x_sen_binning.x_max);
      fH2D[h_name]->GetYaxis()->SetTitle("X [cm]");
      fH2D[h_name]->GetXaxis()->SetTitle(Form("\\rho_{%s} [cm]", di));
 
      h_name = Form("%s_d%s_vs_y", h_name_base.c_str(), di);
      fH2D[h_name] =
        std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), r_sen_binning.n_of_bins, r_sen_binning.x_min,
                               r_sen_binning.x_max, y_sen_binning.n_of_bins, y_sen_binning.x_min, y_sen_binning.x_max);
      fH2D[h_name]->GetYaxis()->SetTitle("Y [cm]");
      fH2D[h_name]->GetXaxis()->SetTitle(Form("\\rho_{%s} [cm]", di));
    }
    // ---- --------- ----
  }  // end geo loop
 
  // ----     3D vertex      ----
  h_name = "vertex_y_vs_x";
  fH2D[h_name] =
    std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), x_vtx_binning.n_of_bins, x_vtx_binning.x_min,
                           x_vtx_binning.x_max, y_vtx_binning.n_of_bins, y_vtx_binning.x_min, y_vtx_binning.x_max);
  fH2D[h_name]->GetXaxis()->SetTitle("Vertex_{X} [cm]");
  fH2D[h_name]->GetYaxis()->SetTitle("Vertex_{Y} [cm]");
 
  h_name = "vertex_x_vs_z";
  fH2D[h_name] =
    std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), z_vtx_binning.n_of_bins, z_vtx_binning.x_min,
                           z_vtx_binning.x_max, x_vtx_binning.n_of_bins, x_vtx_binning.x_min, x_vtx_binning.x_max);
  fH2D[h_name]->GetXaxis()->SetTitle("Vertex_{Z} [cm]");
  fH2D[h_name]->GetYaxis()->SetTitle("Vertex_{X} [cm]");
 
  h_name = "vertex_y_vs_z";
  fH2D[h_name] =
    std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), z_vtx_binning.n_of_bins, z_vtx_binning.x_min,
                           z_vtx_binning.x_max, y_vtx_binning.n_of_bins, y_vtx_binning.x_min, y_vtx_binning.x_max);
  fH2D[h_name]->GetXaxis()->SetTitle("Vertex_{Z} [cm]");
  fH2D[h_name]->GetYaxis()->SetTitle("Vertex_{Y} [cm]");
  // ---- ------------------ ----
 
  // ----        DCA         ----
  for (const char* di : {"x", "y", "z"}) {
    h_name = Form("dca_d%s_vs_x", di);
    fH2D[h_name] =
      std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), r_dca_binning.n_of_bins, r_dca_binning.x_min,
                             r_dca_binning.x_max, x_vtx_binning.n_of_bins, x_vtx_binning.x_min, x_vtx_binning.x_max);
    fH2D[h_name]->GetYaxis()->SetTitle("X [cm]");
    fH2D[h_name]->GetXaxis()->SetTitle(Form("DCA_{%s} [cm]", di));
 
    h_name = Form("dca_d%s_vs_y", di);
    fH2D[h_name] =
      std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), r_dca_binning.n_of_bins, r_dca_binning.x_min,
                             r_dca_binning.x_max, y_vtx_binning.n_of_bins, y_vtx_binning.x_min, y_vtx_binning.x_max);
    fH2D[h_name]->GetYaxis()->SetTitle("Y [cm]");
    fH2D[h_name]->GetXaxis()->SetTitle(Form("DCA_{%s} [cm]", di));
 
    h_name = Form("dca_d%s_vs_z", di);
    fH2D[h_name] =
      std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), r_dca_binning.n_of_bins, r_dca_binning.x_min,
                             r_dca_binning.x_max, z_vtx_binning.n_of_bins, z_vtx_binning.x_min, z_vtx_binning.x_max);
    fH2D[h_name]->GetYaxis()->SetTitle(Form("Z [cm]"));
    fH2D[h_name]->GetXaxis()->SetTitle(Form("DCA_{%s} [cm]", di));
  }
  // ---- ------------------ ----
 
  // ---- Track multiplicity ----
  for (const char* mod : {":all", ":sel"}) {
    h_name       = Form("ca_track_multiplicity%s", mod);
    fH1D[h_name] = std::make_unique<TH1D>(h_name.c_str(), h_name.c_str(), 50, 0, 50);  // HARDCODE
    fH1D[h_name]->GetXaxis()->SetTitle("N_{CATrack}");
    fH1D[h_name]->GetYaxis()->SetTitle("Entries");
 
    h_name       = Form("ca_track_chi2%s", mod);
    fH1D[h_name] = std::make_unique<TH1D>(h_name.c_str(), h_name.c_str(), 10000, 0, 10);  // HARDCODE
    fH1D[h_name]->GetXaxis()->SetTitle("Chi2");
    fH1D[h_name]->GetYaxis()->SetTitle("Entries");
  }
  // ---- ------------------ ----
}
 
void CbmStsEfficiency::CheckEfficiency()
{
  LOG(debug) << "Checking efficiency of STS station " << fTestLayer << " using " << fGlbTrks.size()
             << " CbmGlobalTracks";
  fVertexFinder->SetTracks(fGlbTrks);
  const std::optional<CbmVertex> vertex = fVertexFinder->FindVertex();
 
  if (!vertex.has_value()) return;
 
  if (std::abs(vertex->GetZ()) > fMaxDisTrgVtx) return;
 
  fH2D["vertex_y_vs_x"]->Fill(vertex->GetX(), vertex->GetY());
  fH2D["vertex_x_vs_z"]->Fill(vertex->GetZ(), vertex->GetX());
  fH2D["vertex_y_vs_z"]->Fill(vertex->GetZ(), vertex->GetY());
 
  TVector3 vtx(vertex->GetX(), vertex->GetY(), vertex->GetZ());
  std::map<int32_t, std::vector<bool>> hit_used;
  for (auto trk : fGlbTrks) {
    double trk_tx = trk->GetParamFirst()->GetTx();
    double trk_ty = trk->GetParamFirst()->GetTy();
    double trk_x0 = trk->GetParamFirst()->GetX();
    double trk_y0 = trk->GetParamFirst()->GetY();
    double trk_z0 = trk->GetParamFirst()->GetZ();
 
    // Check if tack comes from vertex
    TVector3 p(trk_x0, trk_y0, trk_z0);
    TVector3 e(trk_tx, trk_ty, 1);
    TVector3 trk_to_vtx = ((vtx - p).Cross(e)) * (1. / e.Mag());
 
    fH2D["dca_dx_vs_x"]->Fill(trk_to_vtx.Px(), vtx.Px());
    fH2D["dca_dx_vs_y"]->Fill(trk_to_vtx.Px(), vtx.Py());
    fH2D["dca_dx_vs_z"]->Fill(trk_to_vtx.Px(), vtx.Pz());
 
    fH2D["dca_dy_vs_x"]->Fill(trk_to_vtx.Py(), vtx.Px());
    fH2D["dca_dy_vs_y"]->Fill(trk_to_vtx.Py(), vtx.Py());
    fH2D["dca_dy_vs_z"]->Fill(trk_to_vtx.Py(), vtx.Pz());
 
    fH2D["dca_dz_vs_x"]->Fill(trk_to_vtx.Pz(), vtx.Px());
    fH2D["dca_dz_vs_y"]->Fill(trk_to_vtx.Pz(), vtx.Py());
    fH2D["dca_dz_vs_z"]->Fill(trk_to_vtx.Pz(), vtx.Pz());
 
    if (trk_to_vtx.Mag() > fMaxDCATrkVtx) continue;
 
    for (int32_t& sensor : fDUTList) {
      auto hits        = fStsHits[sensor];
      size_t n_of_hits = hits.size();
 
      // Extrapolate to the current sensor z-plane
      double sensor_z    = 0.5 * (fStsGeoInfo[sensor][4] + fStsGeoInfo[sensor][5]);
      double predicted_x = trk_x0 + trk_tx * (sensor_z - trk_z0);
      double predicted_y = trk_y0 + trk_ty * (sensor_z - trk_z0);
 
      // Check if the sensor contain the space point
      if (predicted_x < fStsGeoInfo[sensor][0] || predicted_x > fStsGeoInfo[sensor][1]
          || predicted_y < fStsGeoInfo[sensor][2] || predicted_y > fStsGeoInfo[sensor][3]) {
        continue;
      }
 
      int unit = CbmStsAddress::GetElementId(sensor, kStsUnit);
      fH2D[Form("Sts%d_track", unit)]->Fill(predicted_x, predicted_y);
      fH2D[Form("Sts0x%x_track", sensor)]->Fill(predicted_x, predicted_y);
 
      if (n_of_hits) {
        /* Block to search in available hits */
        double closest_id  = 0;
        double dx          = hits[0]->GetX() - predicted_x - fResiduals[sensor].mean.Px();
        double dy          = hits[0]->GetY() - predicted_y - fResiduals[sensor].mean.Py();
        double closest_rho = dx * dx + dy * dy;
 
        for (size_t idx = 1; idx < n_of_hits; idx++) {
          dx         = hits[idx]->GetX() - predicted_x - fResiduals[sensor].mean.Px();
          dy         = hits[idx]->GetY() - predicted_y - fResiduals[sensor].mean.Py();
          double rho = dx * dx + dy * dy;
 
          if (rho < closest_rho) {
            closest_id  = idx;
            closest_rho = rho;
          }
        }
 
        dx = hits[closest_id]->GetX() - predicted_x - fResiduals[sensor].mean.Px();
        dy = hits[closest_id]->GetY() - predicted_y - fResiduals[sensor].mean.Py();
 
        fH2D[Form("Sts%d_dx_vs_x", unit)]->Fill(dx, hits[closest_id]->GetX());
        fH2D[Form("Sts%d_dy_vs_x", unit)]->Fill(dy, hits[closest_id]->GetX());
        fH2D[Form("Sts%d_dx_vs_y", unit)]->Fill(dx, hits[closest_id]->GetY());
        fH2D[Form("Sts%d_dy_vs_y", unit)]->Fill(dy, hits[closest_id]->GetY());
 
        fH2D[Form("Sts0x%x_dx_vs_x", sensor)]->Fill(dx, hits[closest_id]->GetX());
        fH2D[Form("Sts0x%x_dy_vs_x", sensor)]->Fill(dy, hits[closest_id]->GetX());
        fH2D[Form("Sts0x%x_dx_vs_y", sensor)]->Fill(dx, hits[closest_id]->GetY());
        fH2D[Form("Sts0x%x_dy_vs_y", sensor)]->Fill(dy, hits[closest_id]->GetY());
 
        double sensor_resolution =
          fResiduals[sensor].sigma.Mag() != 0 ? fResiduals[sensor].sigma.Mag() : fDefaultResidual;
        if (closest_rho <= 3.0 * sensor_resolution) {
          fH2D[Form("Sts%d_found", unit)]->Fill(predicted_x, predicted_y);
          fH2D[Form("Sts0x%x_found", sensor)]->Fill(predicted_x, predicted_y);
        }
      }
 
 
      /* Block to search in avaliable hit */
 
    }  // end sensor under test lopp
  }    // end track loop
}
 
TH2D* FindInactiveArea(TH2D* h, int dead_size_x = 10, int dead_size_y = 10)
{
  TH2D* dead = (TH2D*) h->Clone();
  dead->Reset();
 
  for (int bin_y = 1; bin_y < dead->GetNbinsY(); bin_y++) {
    for (int bin_x = dead_size_x; bin_x < dead->GetNbinsX() - dead_size_x; bin_x++) {
 
      bool all_dead = true;
      for (int test_bin = -dead_size_x; test_bin <= +dead_size_x; test_bin++) {
        if (h->GetBinContent(bin_x + test_bin, bin_y)) {
          all_dead = false;
          break;
        }
      }
      if (all_dead) {
        for (int test_bin = -dead_size_x; test_bin <= +dead_size_x; test_bin++)
          dead->SetBinContent(bin_x + test_bin, bin_y, 1);
      }
    }
  }
 
  for (int bin_x = 1; bin_x < dead->GetNbinsX(); bin_x++) {
    for (int bin_y = dead_size_y; bin_y < dead->GetNbinsY() - dead_size_y; bin_y++) {
 
      bool all_dead = true;
      for (int test_bin = -dead_size_y; test_bin <= +dead_size_y; test_bin++) {
        if (h->GetBinContent(bin_x, bin_y + test_bin)) {
          all_dead = false;
          break;
        }
      }
      if (all_dead) {
        for (int test_bin = -dead_size_y; test_bin <= +dead_size_y; test_bin++)
          dead->SetBinContent(bin_x, bin_y + test_bin, 1);
      }
    }
  }
 
  return dead;
}
 
void CbmStsEfficiency::Efficiency(uint32_t sensor)
{
  std::string h_name_base = "";
  if (sensor == fTestLayer) {
    h_name_base = Form("Sts%d", sensor);
  }
  else if (sensor > 8) {
    if (CbmStsAddress::GetElementId(sensor, kStsUnit) != fTestLayer) return;
    h_name_base = Form("Sts0x%x", sensor);
  }
  auto found = fH2D[Form("%s_found", h_name_base.c_str())].get();
  auto track = fH2D[Form("%s_track", h_name_base.c_str())].get();
 
  if (found == nullptr || track == nullptr) return;
 
  TH2D* inactive_area = FindInactiveArea(fH2D[Form("%s_hit_map", h_name_base.c_str())].get(), 50, 50);
 
  // Manual integration
  int n_of_found = 0;
  int n_of_track = 0;
  std::vector<double> samples;
 
  for (int bin_y = 50; bin_y < inactive_area->GetNbinsY() - 50; bin_y++) {
    for (int bin_x = 50; bin_x < inactive_area->GetNbinsX() - 50; bin_x++) {
      if (!inactive_area->GetBinContent(bin_x, bin_y)) {
        n_of_found += found->GetBinContent(bin_x, bin_y);
        n_of_track += track->GetBinContent(bin_x, bin_y);
      }
    }
  }
  LOG(debug) << n_of_found << "\t" << n_of_track;
  double efficiency = n_of_track != 0 ? double(n_of_found) / n_of_track : -1;
  if (efficiency != -1) {
    samples.push_back(efficiency);
  }
 
  double mean  = 0;
  double sigma = 0;
  int n        = samples.size();
  if (n == 0) {
    return;
  }
 
  // Take efficiency as mean value
  for (double& val : samples) {
    mean += val;
  }
  mean = double(mean / n);
 
  // Take error as RMS of the values
  for (double& val : samples) {
    sigma += (val - mean) * (val - mean);
  }
  sigma = sqrt(double(sigma / n));
 
  std::cout << Form("0x%x:\t%0.4f\t%0.4f\n", sensor, mean, sigma);
}
 
void CbmStsEfficiency::FinishTask()
{
  for (auto& [element, geo] : fStsGeoInfo) {
    Efficiency(element);
  }
 
  SaveToFile();
}
 
 
void CbmStsEfficiency::ProcessEvent(CbmEvent* evt)
{
  int nb_sts_hits = evt->GetNofData(ECbmDataType::kStsHit);
 
  for (int hit_evt_idx = 0; hit_evt_idx < nb_sts_hits; hit_evt_idx++) {
    int sts_hit_idx = evt->GetIndex(ECbmDataType::kStsHit, hit_evt_idx);
    CbmStsHit* hit  = (CbmStsHit*) fStsHitArray->UncheckedAt(sts_hit_idx);
    ProcessHit(hit);
  }
 
  int nb_of_glob_trk = evt->GetNofData(ECbmDataType::kGlobalTrack);
  for (int evt_glob_trk_idx = 0; evt_glob_trk_idx < nb_of_glob_trk; evt_glob_trk_idx++) {
    int glob_trk_idx         = evt->GetIndex(ECbmDataType::kGlobalTrack, evt_glob_trk_idx);
    CbmGlobalTrack* glob_trk = (CbmGlobalTrack*) fGlbTrkArray->UncheckedAt(glob_trk_idx);
    ProcessGlobalTrack(glob_trk);
  }
 
  fH1D["ca_track_multiplicity:all"]->Fill(nb_of_glob_trk);
  fH1D["ca_track_multiplicity:sel"]->Fill(fGlbTrks.size());
 
  if (fGlbTrks.size()) CheckEfficiency();
 
  fGlbTrks.clear();
  fStsHits.clear();
}
 
void CbmStsEfficiency::ProcessGlobalTrack(CbmGlobalTrack* trk)
{
  auto sts_trk_idx  = trk->GetStsTrackIndex();
  auto rich_trk_idx = trk->GetRichRingIndex();
  auto much_trk_idx = trk->GetMuchTrackIndex();
  auto trd_trk_idx  = trk->GetTrdTrackIndex();
  auto tof_trk_idx  = trk->GetTofTrackIndex();
  double trk_p_val  = TMath::Prob(trk->GetChi2(), trk->GetNDF());
  fH1D["ca_track_chi2:all"]->Fill(trk->GetChi2());
 
  // Apply GlobalTracks cuts
  int32_t mvd_trk_size = sts_trk_idx != -1 && fStsTrkArray != nullptr
                           ? ((CbmStsTrack*) fStsTrkArray->UncheckedAt(sts_trk_idx))->GetNofMvdHits()
                           : 0;
 
  int32_t sts_trk_size = sts_trk_idx != -1 && fStsTrkArray != nullptr
                           ? ((CbmStsTrack*) fStsTrkArray->UncheckedAt(sts_trk_idx))->GetNofStsHits()
                           : 0;
 
  int32_t rich_trk_size = rich_trk_idx != -1 && fRchTrkArray != nullptr
                            ? ((CbmTrack*) fRchTrkArray->UncheckedAt(rich_trk_idx))->GetNofHits()
                            : 0;
 
  int32_t much_trk_size = much_trk_idx != -1 && fMchTrkArray != nullptr
                            ? ((CbmTrack*) fMchTrkArray->UncheckedAt(much_trk_idx))->GetNofHits()
                            : 0;
 
  int32_t trd_trk_size = trd_trk_idx != -1 && fTrdTrkArray != nullptr
                           ? ((CbmTrack*) fTrdTrkArray->UncheckedAt(trd_trk_idx))->GetNofHits()
                           : 0;
 
  int32_t tof_trk_size = tof_trk_idx != -1 && fTofTrkArray != nullptr
                           ? ((CbmTrack*) fTofTrkArray->UncheckedAt(tof_trk_idx))->GetNofHits()
                           : 0;
 
  if (fAnalysisCuts != nullptr
      && (!fAnalysisCuts->Check(CbmCutId::kGlobalTrackMvdSize, mvd_trk_size)
          || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackStsSize, sts_trk_size)
          || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackRichSize, rich_trk_size)
          || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackMuchSize, much_trk_size)
          || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackTrdSize, trd_trk_size)
          || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackTofSize, tof_trk_size)
          || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackChi2, trk->GetChi2())
          || !fAnalysisCuts->Check(CbmCutId::kGlobalTrackPval, trk_p_val))) {
    return;
  }
 
  fH1D["ca_track_chi2:sel"]->Fill(trk->GetChi2());
  fGlbTrks.push_back(trk);
}
 
void CbmStsEfficiency::ProcessHit(CbmStsHit* hit)
{
  if (hit == nullptr) return;
  uint32_t address = hit->GetAddress();
  uint32_t unit    = CbmStsAddress::GetElementId(address, kStsUnit);
 
  if (unit != fTestLayer) return;
 
  if (fAnalysisCuts != nullptr && !fAnalysisCuts->CheckStsHit(hit, fStsCluArray)) return;
 
  fStsHits[address].push_back(hit);
  fH2D[Form("Sts%u_hit_map", unit)]->Fill(hit->GetX(), hit->GetY());
  fH2D[Form("Sts0x%x_hit_map", address)]->Fill(hit->GetX(), hit->GetY());
}
 
void CbmStsEfficiency::Exec(Option_t*)
{
  // Check if CbmEvent
  if (fCbmEvtArray != nullptr) {
    LOG(info) << "Running CbmStsEfficiency - Event-like - " << entry_;
 
    uint32_t nb_events = fCbmEvtArray->GetEntriesFast();
    for (uint32_t evt_idx = 0; evt_idx < nb_events; evt_idx++) {
      CbmEvent* event = (CbmEvent*) fCbmEvtArray->UncheckedAt(evt_idx);
      if (fAnalysisCuts != nullptr && !fAnalysisCuts->CheckEvent(event)) continue;
      LOG(debug) << Form("Process CbmEvent %d / %d", evt_idx, nb_events);
      ProcessEvent(event);
    }
  }
  else {
 
    LOG(info) << "Running CbmStsEfficiency - TS-like - ";
    uint32_t n_of_hits = fStsHitArray->GetEntriesFast();
    for (uint32_t hit_idx = 0; hit_idx < n_of_hits; hit_idx++) {
      CbmStsHit* sts_hit = (CbmStsHit*) fStsHitArray->UncheckedAt(hit_idx);
      ProcessHit(sts_hit);
    }
    uint32_t nb_of_glob_trk = fGlbTrkArray->GetEntriesFast();
    LOG(debug) << Form("Number of GlobalTracks: %u", nb_of_glob_trk);
    for (uint32_t glob_trk_idx = 0; glob_trk_idx < nb_of_glob_trk; glob_trk_idx++) {
      CbmGlobalTrack* glob_trk = (CbmGlobalTrack*) fGlbTrkArray->UncheckedAt(glob_trk_idx);
      ProcessGlobalTrack(glob_trk);
    }
 
    if (fGlbTrks.size()) CheckEfficiency();
 
    fGlbTrks.clear();
    fStsHits.clear();
  }
  entry_++;
}
 
InitStatus CbmStsEfficiency::Init()
{
  LOG(debug) << "Init CbmStsEfficiency ...";
 
  FairRootManager* ioman = FairRootManager::Instance();
  if (ioman != nullptr) {
    fCbmEvtArray = (TClonesArray*) ioman->GetObject("CbmEvent");
 
    fGlbTrkArray = (TClonesArray*) ioman->GetObject("GlobalTrack");
 
    fStsTrkArray = (TClonesArray*) ioman->GetObject("StsTrack");
    fRchTrkArray = (TClonesArray*) ioman->GetObject("RichTrack");
    fMchTrkArray = (TClonesArray*) ioman->GetObject("MuchTrack");
    fTrdTrkArray = (TClonesArray*) ioman->GetObject("TrdTrack");
    fTofTrkArray = (TClonesArray*) ioman->GetObject("TofTrack");
 
    fStsHitArray = (TClonesArray*) ioman->GetObject("StsHit");
 
    fStsCluArray = (TClonesArray*) ioman->GetObject("StsCluster");
  }
 
  LoadSetup();
 
  BookHistograms();
 
  std::stringstream ss;
  for (auto s : fDUTList) {
    ss << Form("\t\t0x%x - U%d L%d M%d\n", s, CbmStsAddress::GetElementId(s, kStsUnit),
               CbmStsAddress::GetElementId(s, kStsLadder), CbmStsAddress::GetElementId(s, kStsModule));
  }
 
  LOG(info) << "Max vtx target dz: " << fMaxDisTrgVtx;
  LOG(info) << "Max trk vtx DCA: " << fMaxDCATrkVtx;
  LOG(info) << "Default residual: " << fDefaultResidual;
  LOG(info) << "Test layer: " << fTestLayer;
  LOG(info) << "Layer modules:\n" << ss.str();
 
  return kSUCCESS;
}
 
void CbmStsEfficiency::SetSensorResidual(int32_t address, TVector3 m, TVector3 s)
{
  fResiduals[address] = Residual(m, s);
}
 
void CbmStsEfficiency::SetResidual(std::string file_name)
{
  std::ifstream in(file_name);
  int32_t sensor;
  double mean_x, mean_y, sigma_x, sigma_y;
  while (in >> std::hex >> sensor >> std::dec >> mean_x >> mean_y >> sigma_x >> sigma_y) {
    SetSensorResidual(sensor, TVector3(mean_x, mean_y, 0), TVector3(sigma_x, sigma_y, 0));
  }
}
 
 
void CbmStsEfficiency::SetVertexFinder(CbmDcaVertexFinder* vtx_finder) { fVertexFinder = vtx_finder; }