CbmStsResolution.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 "CbmStsResolution.h"
 
#include "CbmPixelHit.h"
#include "CbmStsAddress.h"
#include "CbmStsTrack.h"
#include "TPaveStats.h"
#include "TStyle.h"
 
CbmStsResolution::CbmStsResolution(double d_max) : fImpactParMax(d_max)
{
  LOG(debug) << "Creating an instance of CbmStsResolution ...";
}
 
void CbmStsResolution::BookHistograms()
{
  for (auto& [element, geo] : fStsGeoInfo) {
    std::string h_name_base = element < 8 ? Form("Sts%d", element) : Form("Sts0x%x", element);
 
    LOG(debug) << Form("Booking for %s", h_name_base.c_str());
 
    double x_min = geo[0];
    double x_max = geo[1];
    double y_min = geo[2];
    double y_max = geo[3];
 
    unsigned int nb_bins_x = (x_max - x_min) / cbm_sts_utils::kStsDx;
    unsigned int nb_bins_y = (y_max - y_min) / cbm_sts_utils::kStsDy;
 
 
    double dr              = 5e-4;
    double r_min           = -2.7 - 0.5 * dr;
    double r_max           = +2.7 + 0.5 * dr;
    unsigned int nb_bins_r = (r_max - r_min) / dr;
 
    for (auto mod : {":all", ":trk"}) {
      std::string h_name;
      h_name       = Form("%s_dx_vs_x%s", h_name_base.c_str(), mod);
      fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), "", nb_bins_r, r_min, r_max, nb_bins_x, x_min, x_max);
      fH2D[h_name]->GetXaxis()->SetTitle(Form("X_{trk} - X_{hit} [cm]"));
      fH2D[h_name]->GetYaxis()->SetTitle("X_{hit} [cm]");
      fH2D[h_name]->GetXaxis()->CenterTitle(true);
      fH2D[h_name]->GetYaxis()->CenterTitle(true);
 
      h_name       = Form("%s_dx_vs_y%s", h_name_base.c_str(), mod);
      fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), "", nb_bins_r, r_min, r_max, nb_bins_y, y_min, y_max);
      fH2D[h_name]->GetXaxis()->SetTitle(Form("X_{trk} - X_{hit} [cm]"));
      fH2D[h_name]->GetYaxis()->SetTitle("Y_{hit} [cm]");
      fH2D[h_name]->GetXaxis()->CenterTitle(true);
      fH2D[h_name]->GetYaxis()->CenterTitle(true);
 
      h_name       = Form("%s_dy_vs_x%s", h_name_base.c_str(), mod);
      fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), "", nb_bins_r, r_min, r_max, nb_bins_x, x_min, x_max);
      fH2D[h_name]->GetXaxis()->SetTitle(Form("Y_{trk} - Y_{hit} [cm]"));
      fH2D[h_name]->GetYaxis()->SetTitle("X_{hit} [cm]");
      fH2D[h_name]->GetXaxis()->CenterTitle(true);
      fH2D[h_name]->GetYaxis()->CenterTitle(true);
 
      h_name       = Form("%s_dy_vs_y%s", h_name_base.c_str(), mod);
      fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), "", nb_bins_r, r_min, r_max, nb_bins_y, y_min, y_max);
      fH2D[h_name]->GetXaxis()->SetTitle(Form("Y_{trk} - Y_{hit} [cm]"));
      fH2D[h_name]->GetYaxis()->SetTitle("Y_{hit} [cm]");
      fH2D[h_name]->GetXaxis()->CenterTitle(true);
      fH2D[h_name]->GetYaxis()->CenterTitle(true);
    }
    if (element < 8) {
      std::string h_name = Form("Sts%d_ref_hits", element);
      fH2D[h_name]       = std::make_unique<TH2D>(h_name.c_str(), "", nb_bins_x, x_min, x_max, nb_bins_y, y_min, y_max);
      fH2D[h_name]->GetXaxis()->SetTitle(Form("X_{TrkHit} [cm]"));
      fH2D[h_name]->GetYaxis()->SetTitle("Y_{TrkHit} [cm]");
      fH2D[h_name]->GetXaxis()->CenterTitle(true);
      fH2D[h_name]->GetYaxis()->CenterTitle(true);
    }
  }
}
 
void CbmStsResolution::ProcessEvent(CbmEvent* evt)
{
  if (fAnalysisCuts != nullptr && !fAnalysisCuts->CheckEvent(evt)) return;
 
  int nb_sts_hits = evt->GetNofData(ECbmDataType::kStsHit);
  LOG(debug) << Form("Processing event with %d StsHit", nb_sts_hits);
  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);
  }
  BuildResidual();
}
 
void CbmStsResolution::BuildResidual()
{
  int n_of_units = fStsHits.size();
 
  for (auto trk : fGlbTrks) {
    auto sts_trk_idx = trk->GetStsTrackIndex();
    if (sts_trk_idx == -1) continue;
 
    CbmStsTrack* sts_track = (CbmStsTrack*) fStsTrkArray->UncheckedAt(sts_trk_idx);
    if (((CbmStsTrack*) fStsTrkArray->UncheckedAt(sts_trk_idx))->GetNofStsHits() != 2) continue;
 
    CbmStsHit* hit_i = (CbmStsHit*) fStsHitArray->UncheckedAt(sts_track->GetStsHitIndex(0));
    CbmStsHit* hit_j = (CbmStsHit*) fStsHitArray->UncheckedAt(sts_track->GetStsHitIndex(1));
 
    // Charge cut to Sts track hits
    if (fTrkHitQmin > 0
        && (cbm_sts_utils::GetHitCharge(hit_i, fStsCluArray) < fTrkHitQmin
            || cbm_sts_utils::GetHitCharge(hit_j, fStsCluArray) < fTrkHitQmin))
      continue;
 
    int unit_i = CbmStsAddress::GetElementId(hit_i->GetAddress(), kStsUnit);
    int unit_j = CbmStsAddress::GetElementId(hit_j->GetAddress(), kStsUnit);
 
    // Checkk track crossed eneabled sensors
    if (fActiveRefSensors.size()) {
      if (std::find(fActiveRefSensors.begin(), fActiveRefSensors.end(), hit_i->GetAddress()) == fActiveRefSensors.end())
        continue;
      if (std::find(fActiveRefSensors.begin(), fActiveRefSensors.end(), hit_j->GetAddress()) == fActiveRefSensors.end())
        continue;
    }
 
    double x_i = hit_i->GetX();
    double y_i = hit_i->GetY();
    double z_i = hit_i->GetZ();
 
    double x_j = hit_j->GetX();
    double y_j = hit_j->GetY();
    double z_j = hit_j->GetZ();
 
    fH2D[Form("Sts%d_ref_hits", unit_i)]->Fill(x_i, y_i);
    fH2D[Form("Sts%d_ref_hits", unit_j)]->Fill(x_j, y_j);
 
    double x_0, y_0, t_x, t_y;
    if (tracking_use_two_sts_) {
      t_x = (x_i - x_j) / (z_i - z_j);
      t_y = (y_i - y_j) / (z_i - z_j);
      x_0 = x_i - t_x * z_i;
      y_0 = y_i - t_y * z_i;
    }
    else {
      t_x = trk->GetParamFirst()->GetTx();
      t_y = trk->GetParamFirst()->GetTy();
      x_0 = x_i - t_x * z_i;
      y_0 = y_i - t_y * z_i;
    }
 
    double d_trk_trg = TMath::Sqrt(x_0 * x_0 + y_0 * y_0);
 
    if (fImpactParMax > 0 && fImpactParMax < d_trk_trg) continue;
 
    for (int uut = 0; uut < n_of_units; uut++) {
      if (uut == unit_i || uut == unit_j) continue;
 
      CbmStsHit* closest_hit = nullptr;
      double min_dist        = 1e+6;
      for (auto hit_dut : fStsHits[uut]) {
        double x_dut = hit_dut->GetX();
        double y_dut = hit_dut->GetY();
        double z_dut = hit_dut->GetZ();
 
        double x_trk = x_i + t_x * (z_dut - z_i);
        double y_trk = y_i + t_y * (z_dut - z_i);
 
        double dx = x_trk - x_dut;
        double dy = y_trk - y_dut;
 
        fH2D[Form("Sts0x%x_dx_vs_x:all", hit_dut->GetAddress())]->Fill(dx, x_dut);
        fH2D[Form("Sts0x%x_dx_vs_y:all", hit_dut->GetAddress())]->Fill(dx, y_dut);
        fH2D[Form("Sts0x%x_dy_vs_x:all", hit_dut->GetAddress())]->Fill(dy, x_dut);
        fH2D[Form("Sts0x%x_dy_vs_y:all", hit_dut->GetAddress())]->Fill(dy, y_dut);
 
        fH2D[Form("Sts%d_dx_vs_x:all", uut)]->Fill(dx, x_dut);
        fH2D[Form("Sts%d_dx_vs_y:all", uut)]->Fill(dx, y_dut);
        fH2D[Form("Sts%d_dy_vs_x:all", uut)]->Fill(dy, x_dut);
        fH2D[Form("Sts%d_dy_vs_y:all", uut)]->Fill(dy, y_dut);
 
        if (min_dist > dx * dx + dy * dy) {
          min_dist    = dx * dx + dy * dy;
          closest_hit = hit_dut;
        }
      }  // hit_dut loop
 
      if (closest_hit) {
        double x_dut = closest_hit->GetX();
        double y_dut = closest_hit->GetY();
        double z_dut = closest_hit->GetZ();
 
        double x_trk = x_i + t_x * (z_dut - z_i);
        double y_trk = y_i + t_y * (z_dut - z_i);
 
        double dx = x_trk - x_dut;
        double dy = y_trk - y_dut;
 
        fH2D[Form("Sts0x%x_dx_vs_x:trk", closest_hit->GetAddress())]->Fill(dx, x_dut);
        fH2D[Form("Sts0x%x_dx_vs_y:trk", closest_hit->GetAddress())]->Fill(dx, y_dut);
        fH2D[Form("Sts0x%x_dy_vs_x:trk", closest_hit->GetAddress())]->Fill(dy, x_dut);
        fH2D[Form("Sts0x%x_dy_vs_y:trk", closest_hit->GetAddress())]->Fill(dy, y_dut);
 
        fH2D[Form("Sts%d_dx_vs_x:trk", uut)]->Fill(dx, x_dut);
        fH2D[Form("Sts%d_dx_vs_y:trk", uut)]->Fill(dx, y_dut);
        fH2D[Form("Sts%d_dy_vs_x:trk", uut)]->Fill(dy, x_dut);
        fH2D[Form("Sts%d_dy_vs_y:trk", uut)]->Fill(dy, y_dut);
      }
    }  // uut loop
  }    // track loop
 
  fStsHits.clear();
  fGlbTrks.clear();
}
 
void CbmStsResolution::ProcessHit(CbmStsHit* hit)
{
  if (hit == nullptr) return;
 
  if (fAnalysisCuts != nullptr && !fAnalysisCuts->CheckStsHit(hit, fStsCluArray)) return;
 
  int32_t address = hit->GetAddress();
  int32_t unit    = CbmStsAddress::GetElementId(address, kStsUnit);
 
  fStsHits[unit].push_back(hit);
}
 
void CbmStsResolution::EnableRefSensor(int32_t address) { fActiveRefSensors.push_back(address); }
 
void CbmStsResolution::ProcessGlobalTrack(CbmGlobalTrack* trk)
{
  if (trk == nullptr) return;
 
  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();
  float trk_p_val   = TMath::Prob(trk->GetChi2(), trk->GetNDF());
 
  // 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;
  }
 
  fGlbTrks.push_back(trk);
}
 
void CbmStsResolution::Exec(Option_t*)
{
  switch (fInputType) {
    case kEventMode: {
      if (fCbmEvtArray == nullptr) {
        LOG(error) << "Invalid branch for event-mode: Branch CbmEvent -> nullptr";
        break;
      }
 
      uint32_t nb_events = fCbmEvtArray->GetEntriesFast();
      LOG(info) << "Processing entry " << entry_ << "\tCbmEvent: " << nb_events;
 
      for (uint32_t evt_idx = 0; evt_idx < nb_events; evt_idx++) {
        CbmEvent* event = (CbmEvent*) fCbmEvtArray->UncheckedAt(evt_idx);
        ProcessEvent(event);
      }
      break;
    }
    case kMCTimeMode: {
      break;
    }
    case kMCEventMode:
    case kTimeMode:
    default: {
      uint32_t n_of_hits    = fStsHitArray->GetEntriesFast();
      uint32_t n_of_glb_trk = fGlbTrkArray->GetEntriesFast();
      LOG(info) << "Processing entry " << entry_ << "\tStsHit: " << n_of_hits << "\tGlobalTrack: " << n_of_glb_trk;
 
      for (uint32_t hit_idx = 0; hit_idx < n_of_hits; hit_idx++) {
        ProcessHit((CbmStsHit*) fStsHitArray->UncheckedAt(hit_idx));
      }
 
      for (uint32_t glb_trk_idx = 0; glb_trk_idx < n_of_glb_trk; glb_trk_idx++) {
        ProcessGlobalTrack((CbmGlobalTrack*) fGlbTrkArray->UncheckedAt(glb_trk_idx));
      }
 
      BuildResidual();
      break;
    }
  }
  entry_++;
}
 
InitStatus CbmStsResolution::Init()
{
  LOG(debug) << "Init CbmStsResolution ...";
 
  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();
 
  return kSUCCESS;
}
 
void CbmStsResolution::Finish() { SaveToFile(); }