CbmStsHitAna.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 "CbmStsHitAna.h"
 
#include <CbmStsTrack.h>
 
#include <ctime>
#include <iostream>
#include <typeinfo>
 
CbmStsHitAna::CbmStsHitAna() : fHitModifier({":all"}) {}
 
CbmStsHitAna::CbmStsHitAna(std::string cal_par_file) : fHitModifier({":all"}), fCalibrationFile(cal_par_file) {}
 
InitStatus CbmStsHitAna::Init()
{
  if (fModuleParSet == nullptr) {
    LOG(debug) << "Loading charge calibration ...";
    fModuleParSet = std::make_unique<CbmStsParSetModule>("CbmParSetModule", "STS parameters", "Default");
    fModuleParSet->LoadParASCII(fCalibrationFile);
  }
 
  FairRootManager* ioman = FairRootManager::Instance();
  if (ioman == nullptr) return kERROR;
 
  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");
 
  if (fGlbTrkArray) fHitModifier.push_back(":trk");
 
  LoadSetup();
 
  return kSUCCESS;
}
 
void CbmStsHitAna::BookHistograms(int32_t address)
{
  LOG(debug) << Form("Booking histograms for module: 0x%x", address);
 
  const CbmStsParModule& par_module           = fModuleParSet->GetParModule(address);
  const auto [n_side_binning, p_side_binning] = cbm_sts_utils::ChargeBinning(par_module, fMaxCluSize);
 
  std::string h_name;
 
  double x_min = fStsGeoInfo.count(address) ? fStsGeoInfo[address][0] : -15;
  double x_max = fStsGeoInfo.count(address) ? fStsGeoInfo[address][1] : +15;
  double y_min = fStsGeoInfo.count(address) ? fStsGeoInfo[address][2] : -15;
  double y_max = fStsGeoInfo.count(address) ? fStsGeoInfo[address][3] : +15;
 
  double t_min = -150 - 0.5 * cbm_sts_utils::kStsClock;
  double t_max = +150 + 0.5 * cbm_sts_utils::kStsClock;
 
  auto x_binning = cbm_sts_utils::HBinning{uint32_t((x_max - x_min) / cbm_sts_utils::kStsDx), x_min, x_max};
  auto y_binning = cbm_sts_utils::HBinning{uint32_t((y_max - y_min) / cbm_sts_utils::kStsDy), y_min, y_max};
  auto t_binning = cbm_sts_utils::HBinning{uint32_t((t_max - t_min) / cbm_sts_utils::kStsClock), t_min, t_max};
 
  for (const char* mod : fHitModifier) {
    h_name       = Form("0x%x_Q_asymmetry%s", address, mod);
    fH1D[h_name] = std::make_unique<TH1D>(h_name.c_str(), h_name.c_str(), 1000, -1.01, +1.01);
 
    h_name       = Form("0x%x_Qp_vs_Qn%s", address, mod);
    fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), n_side_binning.n_of_bins,
                                          n_side_binning.x_min, n_side_binning.x_max, p_side_binning.n_of_bins,
                                          p_side_binning.x_min, p_side_binning.x_max);
 
    h_name       = Form("0x%x_Qp_vs_size%s", address, mod);
    fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), p_side_binning.n_of_bins,
                                          p_side_binning.x_min, p_side_binning.x_max, fMaxCluSize, 1, fMaxCluSize + 1);
 
    h_name       = Form("0x%x_Qn_vs_size%s", address, mod);
    fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), n_side_binning.n_of_bins,
                                          n_side_binning.x_min, n_side_binning.x_max, fMaxCluSize, 1, fMaxCluSize + 1);
 
    h_name       = Form("0x%x_psize_vs_nsize%s", address, mod);
    fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), fMaxCluSize, 1, fMaxCluSize + 1, fMaxCluSize,
                                          1, fMaxCluSize + 1);
 
    h_name       = Form("0x%x_y_vs_x%s", address, mod);
    fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), x_binning.n_of_bins, x_binning.x_min,
                                          x_binning.x_max, y_binning.n_of_bins, y_binning.x_min, y_binning.x_max);
 
    h_name = Form("0x%x_cluster_dt%s", address, mod);
    fH1D[h_name] =
      std::make_unique<TH1D>(h_name.c_str(), h_name.c_str(), t_binning.n_of_bins, t_binning.x_min, t_binning.x_max);
  }
}
 
void CbmStsHitAna::ProcessEvent(CbmEvent* evt)
{
  if (fAnalysisCuts != nullptr && !fAnalysisCuts->CheckEvent(evt)) return;
 
  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, false);
  }
 
  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);
  }
}
 
void CbmStsHitAna::ProcessHit(CbmStsHit* hit, bool belong_to_trk)
{
  if (hit == nullptr) return;
  int32_t address = hit->GetAddress();
 
  if (fAnalysisCuts != nullptr && !fAnalysisCuts->CheckStsHit(hit, fStsCluArray)) {
    return;
  }
 
  if (!fAddressBook.count(address)) {
    fAddressBook.insert(address);
    BookHistograms(address);
  }
 
  double q_n = cbm_sts_utils::GetHitChargeF(hit, fStsCluArray);
  double q_p = cbm_sts_utils::GetHitChargeB(hit, fStsCluArray);
 
  int32_t size_n = cbm_sts_utils::GetHitCluSizeF(hit, fStsCluArray);
  int32_t size_p = cbm_sts_utils::GetHitCluSizeB(hit, fStsCluArray);
 
  const char* mod = fHitModifier[belong_to_trk];
 
  fH1D[Form("0x%x_Q_asymmetry%s", address, mod)]->Fill(cbm_sts_utils::GetHitChargeAsy(hit, fStsCluArray));
  fH2D[Form("0x%x_Qp_vs_Qn%s", address, mod)]->Fill(q_n, q_p);
  fH2D[Form("0x%x_Qp_vs_size%s", address, mod)]->Fill(q_p, size_p);
  fH2D[Form("0x%x_Qn_vs_size%s", address, mod)]->Fill(q_n, size_n);
  fH2D[Form("0x%x_psize_vs_nsize%s", address, mod)]->Fill(size_p, size_n);
  fH2D[Form("0x%x_y_vs_x%s", address, mod)]->Fill(hit->GetX(), hit->GetY());
  fH1D[Form("0x%x_cluster_dt%s", address, mod)]->Fill(cbm_sts_utils::GetHitTimeF(hit, fStsCluArray)
                                                      - cbm_sts_utils::GetHitTimeB(hit, fStsCluArray));
}
 
void CbmStsHitAna::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();
  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;
  }
  LOG(debug) << Form("Processing %d StsHit that were attached to a CATrack", sts_trk_size);
  CbmStsTrack* sts_track = (CbmStsTrack*) fStsTrkArray->UncheckedAt(sts_trk_idx);
 
  // ------------------------------------
  // Process StsHit attached to a CATrack
  // ------------------------------------
  for (int hit_idx = 0; hit_idx < sts_trk_size; hit_idx++) {
    uint32_t sts_hit_idx = sts_track->GetStsHitIndex(hit_idx);
    CbmStsHit* sts_hit   = (CbmStsHit*) fStsHitArray->UncheckedAt(sts_hit_idx);
    ProcessHit(sts_hit, true);
  }
  // ------------------------------------
}
 
void CbmStsHitAna::Exec(Option_t*)
{
  if (fCbmEvtArray != nullptr) {
    uint32_t nb_events        = fCbmEvtArray->GetEntriesFast();
    double avg_nb_sts_hits    = 0;
    double avg_nb_of_glob_trk = 0;
    for (uint32_t evt_idx = 0; evt_idx < nb_events; evt_idx++) {
      ProcessEvent((CbmEvent*) fCbmEvtArray->UncheckedAt(evt_idx));
      avg_nb_sts_hits += ((CbmEvent*) fCbmEvtArray->UncheckedAt(evt_idx))->GetNofData(ECbmDataType::kStsHit);
      avg_nb_of_glob_trk += ((CbmEvent*) fCbmEvtArray->UncheckedAt(evt_idx))->GetNofData(ECbmDataType::kGlobalTrack);
    }
    LOG_IF(info, nb_events) << "Running CbmStsHitAna - Event-like - Entry: " << entry_ << "\tEvents: " << nb_events
                            << "\tAvg StsHit/Events: " << avg_nb_sts_hits / nb_events
                            << "\tAvg GlobalTrk/Events: " << avg_nb_of_glob_trk / nb_events;
    LOG_IF(info, !nb_events) << "Running CbmStsHitAna - Event-like - Entry: " << entry_ << "\tEvents: " << nb_events;
  }
  else {
    uint32_t n_of_hits = fStsHitArray->GetEntriesFast();
    LOG(info) << "Running CbmStsHitAna - Time-like - Entry: " << entry_ << "\tStsHits: " << n_of_hits;
    for (uint32_t hit_idx = 0; hit_idx < n_of_hits; hit_idx++) {
      CbmStsHit* sts_hit = (CbmStsHit*) fStsHitArray->UncheckedAt(hit_idx);
      ProcessHit(sts_hit, false);
    }
  }
  entry_++;
}
 
void CbmStsHitAna::Finish()
{
  SaveToFile();
  fH1D.clear();
  fH2D.clear();
}