CbmStsCorrelation.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 "CbmStsCorrelation.h"
 
#include "CbmPixelHit.h"
#include "CbmStsAddress.h"
#include "TPaveStats.h"
#include "TStyle.h"
 
void CbmStsCorrelation::BookHistograms()
{
 
  for (auto& [element_i, geo_i] : fStsGeoInfo) {
    for (auto& [element_j, geo_j] : fStsGeoInfo) {
 
      std::string h_name_base = "";
      std::string obj_x, obj_y;
      if (element_i < 8 && element_j < 8) {
        if (element_i >= element_j) continue;
        obj_x       = Form("Sts%d", element_i);
        obj_y       = Form("Sts%d", element_j);
        h_name_base = Form("Sts%d:Sts%d", element_i, element_j);
      }
      else if (element_i > 8 && element_j > 8) {
        int32_t unit_i = CbmStsAddress::GetElementId(element_i, kStsUnit);
        int32_t unit_j = CbmStsAddress::GetElementId(element_j, kStsUnit);
 
        if (unit_i >= unit_j) continue;
 
        obj_x       = Form("Sts0x%x", element_i);
        obj_y       = Form("Sts0x%x", element_j);
        h_name_base = Form("Sts0x%x:Sts0x%x", element_i, element_j);
      }
 
      if (!h_name_base.length()) continue;
      LOG(debug) << Form("Booking for %s", h_name_base.c_str());
 
      double x_min_i = geo_i[0];
      double x_max_i = geo_i[1];
      double y_min_i = geo_i[2];
      double y_max_i = geo_i[3];
 
      unsigned int nb_bins_x_i = (x_max_i - x_min_i) / cbm_sts_utils::kStsDx;
      unsigned int nb_bins_y_i = (y_max_i - y_min_i) / cbm_sts_utils::kStsDy;
 
      double x_min_j = geo_j[0];
      double x_max_j = geo_j[1];
      double y_min_j = geo_j[2];
      double y_max_j = geo_j[3];
 
      LOG(debug) << Form("XX_[%0.3f, %0.3f]:[%0.3f, %0.3f]", x_min_i, x_max_i, x_min_j, x_max_j);
      LOG(debug) << Form("YY_[%0.3f, %0.3f]:[%0.3f, %0.3f]", y_min_i, y_max_i, y_min_j, y_max_j);
 
      unsigned int nb_bins_x_j = (x_max_j - x_min_j) / cbm_sts_utils::kStsDx;
      unsigned int nb_bins_y_j = (y_max_j - y_min_j) / cbm_sts_utils::kStsDy;
 
      std::string h_name;
      h_name = Form("%s_XX", h_name_base.c_str());
      fH2D[h_name] =
        std::make_unique<TH2D>(h_name.c_str(), "", nb_bins_x_i, x_min_i, x_max_i, nb_bins_x_j, x_min_j, x_max_j);
      fH2D[h_name]->GetXaxis()->SetTitle(Form("X_{%s} [cm]", obj_x.c_str()));
      fH2D[h_name]->GetYaxis()->SetTitle(Form("X_{%s} [cm]", obj_y.c_str()));
      fH2D[h_name]->GetXaxis()->CenterTitle(true);
      fH2D[h_name]->GetYaxis()->CenterTitle(true);
 
      h_name = Form("%s_YY", h_name_base.c_str());
      fH2D[h_name] =
        std::make_unique<TH2D>(h_name.c_str(), "", nb_bins_y_i, y_min_i, y_max_i, nb_bins_y_j, y_min_j, y_max_j);
      fH2D[h_name]->GetXaxis()->SetTitle(Form("Y_{%s} [cm]", obj_x.c_str()));
      fH2D[h_name]->GetYaxis()->SetTitle(Form("Y_{%s} [cm]", obj_y.c_str()));
      fH2D[h_name]->GetXaxis()->CenterTitle(true);
      fH2D[h_name]->GetYaxis()->CenterTitle(true);
    }
  }
}
 
void CbmStsCorrelation::BuildCorrelation()
{
  bool is_mc_data = false;
  if (!fStsHits.size()) return;
  for (int sts_unit_i = 0; sts_unit_i < nb_sts_station_ - 1; sts_unit_i++) {
    for (int sts_unit_j = sts_unit_i + 1; sts_unit_j < nb_sts_station_; sts_unit_j++) {
      LOG(debug) << Form("Building correlations for STS%d:STS%d: %ld\t%ld", sts_unit_i, sts_unit_j,
                         fStsHits[sts_unit_i].size(), fStsHits[sts_unit_j].size());
      for (auto sts_hit_i : fStsHits[sts_unit_i]) {
        int32_t sts_address_i = sts_hit_i->GetAddress();
        double sts_i_x        = sts_hit_i->GetX();
        double sts_i_y        = sts_hit_i->GetY();
 
        CbmStsHit* closest_hit = nullptr;
        double abs_time_diff   = 999;
        for (auto sts_hit_j : fStsHits[sts_unit_j]) {
 
          if (is_mc_data) {
            int32_t sts_address_j = sts_hit_j->GetAddress();
            double sts_j_x        = sts_hit_j->GetX();
            double sts_j_y        = sts_hit_j->GetY();
 
            fH2D[Form("Sts%d:Sts%d_YY", sts_unit_i, sts_unit_j)]->Fill(sts_i_y, sts_j_y);
            fH2D[Form("Sts%d:Sts%d_XX", sts_unit_i, sts_unit_j)]->Fill(sts_i_x, sts_j_x);
            fH2D[Form("Sts0x%x:Sts0x%x_YY", sts_address_i, sts_address_j)]->Fill(sts_i_y, sts_j_y);
            fH2D[Form("Sts0x%x:Sts0x%x_XX", sts_address_i, sts_address_j)]->Fill(sts_i_x, sts_j_x);
          }
          else {
 
            double dt = std::abs(sts_hit_i->GetTime() - sts_hit_j->GetTime());
            if (abs_time_diff > dt) {
              abs_time_diff = dt;
              closest_hit   = sts_hit_j;
            }
          }
        }
        if (!is_mc_data && closest_hit != nullptr) {
          int32_t sts_address_j = closest_hit->GetAddress();
          double sts_j_x        = closest_hit->GetX();
          double sts_j_y        = closest_hit->GetY();
 
          fH2D[Form("Sts%d:Sts%d_YY", sts_unit_i, sts_unit_j)]->Fill(sts_i_y, sts_j_y);
          fH2D[Form("Sts%d:Sts%d_XX", sts_unit_i, sts_unit_j)]->Fill(sts_i_x, sts_j_x);
          fH2D[Form("Sts0x%x:Sts0x%x_YY", sts_address_i, sts_address_j)]->Fill(sts_i_y, sts_j_y);
          fH2D[Form("Sts0x%x:Sts0x%x_XX", sts_address_i, sts_address_j)]->Fill(sts_i_x, sts_j_x);
        }
      }
    }
  }
  fStsHits.clear();
}
 
void CbmStsCorrelation::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);
  }
}
 
void CbmStsCorrelation::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 CbmStsCorrelation::Exec(Option_t*)
{
  // Check if CbmEvent
  if (fCbmEvtArray != nullptr) {
    uint32_t nb_events = fCbmEvtArray->GetEntriesFast();
    LOG(info) << Form("Running event-like Entry: %d\t%d CbmEvents", entry_, nb_events);
    for (uint32_t evt_idx = 0; evt_idx < nb_events; evt_idx++) {
      CbmEvent* event = (CbmEvent*) fCbmEvtArray->UncheckedAt(evt_idx);
      ProcessEvent(event);
      BuildCorrelation();
    }
  }
  else {
 
    LOG(debug) << "Running time-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);
    }
    BuildCorrelation();
  }
  entry_++;
}
 
InitStatus CbmStsCorrelation::Init()
{
  LOG(debug) << "Init CbmStsCorrelation ...";
 
  FairRootManager* ioman = FairRootManager::Instance();
  if (ioman != nullptr) {
    fCbmEvtArray = (TClonesArray*) ioman->GetObject("CbmEvent");
    fStsHitArray = (TClonesArray*) ioman->GetObject("StsHit");
    fStsCluArray = (TClonesArray*) ioman->GetObject("StsCluster");
  }
 
  LoadSetup();
 
  BookHistograms();
 
  return kSUCCESS;
}
 
 
void CbmStsCorrelation::Finish()
{
  SaveToFile();
  fH1D.clear();
  fH2D.clear();
}