CbmStsRecoBeamSpot.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 "CbmStsRecoBeamSpot.h"
 
#include "CbmStsAddress.h"
 
void CbmStsRecoBeamSpot::AddTarget(CbmTarget* trg) { AddTarget("", trg); }
void CbmStsRecoBeamSpot::AddTarget(std::string trg_name, CbmTarget* trg)
{
  if (trg_name.length() == 0) {
    trg_name = Form("target_%ld", fTargets.size());
  }
  fTargets[trg_name] = trg;
}
 
void CbmStsRecoBeamSpot::BookHistograms()
{
  std::string h_name;
  std::string h_name_base = "";
  for (auto& [element_i, geo_i] : fStsGeoInfo) {
    if (element_i < 8) continue;
 
    for (auto& [element_j, geo_j] : fStsGeoInfo) {
      if (element_j < 8) continue;
      int32_t unit_i = CbmStsAddress::GetElementId(element_i, kStsUnit);
      int32_t unit_j = CbmStsAddress::GetElementId(element_j, kStsUnit);
 
      if (unit_i >= unit_j) continue;
 
      for (auto& [trg_name, trg] : fTargets) {
        h_name_base = Form("0x%x:0x%x_%s", element_i, element_j, trg_name.c_str());
 
        fSampleRangeXmin = trg->GetPosition().Px() - 10.;  // HARDCODE
        fSampleRangeXmax = trg->GetPosition().Px() + 10.;  // HARDCODE
        fSampleRangeYmin = trg->GetPosition().Py() - 10.;  // HARDCODE
        fSampleRangeYmax = trg->GetPosition().Py() + 10.;  // HARDCODE
        fSampleRangeZmin = trg->GetPosition().Pz() - 10.;  // HARDCODE
        fSampleRangeZmax = trg->GetPosition().Pz() + 10.;  // HARDCODE
 
 
        unsigned int nb_bins_x = (fSampleRangeXmax - fSampleRangeXmin) / fSampleBinSizeX;
        unsigned int nb_bins_y = (fSampleRangeYmax - fSampleRangeYmin) / fSampleBinSizeY;
        unsigned int nb_bins_z = (fSampleRangeZmax - fSampleRangeZmin) / fSampleBinSizeZ;
 
        h_name       = Form("%s_Y_vs_X", h_name_base.c_str());
        fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), nb_bins_x, fSampleRangeXmin,
                                              fSampleRangeXmax, nb_bins_y, fSampleRangeYmin, fSampleRangeYmax);
        fH2D[h_name]->GetXaxis()->SetTitle("X [cm]");
        fH2D[h_name]->GetYaxis()->SetTitle("Y [cm]");
 
        h_name       = Form("%s_X_vs_Z", h_name_base.c_str());
        fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), nb_bins_z, fSampleRangeZmin,
                                              fSampleRangeZmax, nb_bins_x, fSampleRangeXmin, fSampleRangeXmax);
        fH2D[h_name]->GetXaxis()->SetTitle("Z [cm]");
        fH2D[h_name]->GetYaxis()->SetTitle("X [cm]");
 
        h_name       = Form("%s_Y_vs_Z", h_name_base.c_str());
        fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), nb_bins_z, fSampleRangeZmin,
                                              fSampleRangeZmax, nb_bins_y, fSampleRangeYmin, fSampleRangeYmax);
        fH2D[h_name]->GetXaxis()->SetTitle("Z [cm]");
        fH2D[h_name]->GetYaxis()->SetTitle("Y [cm]");
 
        h_name       = Form("%s_Y_beam_vs_X_beam", h_name_base.c_str());
        fH2D[h_name] = std::make_unique<TH2D>(h_name.c_str(), h_name.c_str(), nb_bins_x, fSampleRangeXmin,
                                              fSampleRangeXmax, nb_bins_y, fSampleRangeYmin, fSampleRangeYmax);
        fH2D[h_name]->GetXaxis()->SetTitle("X_beam [cm]");
        fH2D[h_name]->GetYaxis()->SetTitle("Y_beam [cm]");
 
      }  // end targets loop
    }    // end sts_i modules loop
  }      // end sts_j modules loop
}
 
TVector3 CbmStsRecoBeamSpot::ExtrapolateTrackTo(CbmPixelHit* hit_0, CbmPixelHit* hit_1, CbmTarget* trg)
{
  TVector3 n0 = trg->GetPosition();
  TVector3 n1(.0, .0, 1.);
  n1.RotateY(trg->GetRotation());
 
  TVector3 l0(hit_0->GetX(), hit_0->GetY(), hit_0->GetZ());
  TVector3 l1(hit_1->GetX() - hit_0->GetX(), hit_1->GetY() - hit_0->GetY(), hit_1->GetZ() - hit_0->GetZ());
 
  double l1_dot_n1 = l1.Dot(n1);
  if (l1_dot_n1 <= 1e-6) {  // tracklet is parallel (or close to it) to target plane
    throw std::invalid_argument("Track-let parallel to target plane");
  }
 
  double t = (n0 - l0).Dot(n1) / l1.Dot(n1);
  return l0 + t * l1;
}
 
void CbmStsRecoBeamSpot::BeamSpotReco()
{
  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++) {
      for (auto sts_hit_i : fStsHits[sts_unit_i]) {
        for (auto sts_hit_j : fStsHits[sts_unit_j]) {
          int32_t sts_address_i = sts_hit_i->GetAddress();
          int32_t sts_address_j = sts_hit_j->GetAddress();
 
          for (auto& [trg_name, trg] : fTargets) {
            TVector3 sol            = ExtrapolateTrackTo(sts_hit_i, sts_hit_j, trg);
            std::string h_name_base = Form("0x%x:0x%x_%s", sts_address_i, sts_address_j, trg_name.c_str());
 
            std::string h_name = Form("%s_X_vs_Z", h_name_base.c_str());
            fH2D[h_name]->Fill(sol.Pz(), sol.Px());
 
            h_name = Form("%s_Y_vs_Z", h_name_base.c_str());
            fH2D[h_name]->Fill(sol.Pz(), sol.Py());
 
            h_name = Form("%s_Y_vs_X", h_name_base.c_str());
            fH2D[h_name]->Fill(sol.Px(), sol.Py());
 
            TVector3 sol_beam = sol - trg->GetPosition();
 
            double x_beam = sol_beam.Px() / cos(trg->GetRotation());
            double y_beam = sol_beam.Py();
            h_name        = Form("%s_Y_beam_vs_X_beam", h_name_base.c_str());
            fH2D[h_name]->Fill(x_beam, y_beam);
          }
        }
      }
    }
  }
  fStsHits.clear();
}
 
 
void CbmStsRecoBeamSpot::FinishTask() { SaveToFile(); }
 
 
void CbmStsRecoBeamSpot::ProcessEvent(CbmEvent* evt)
{
  int nb_sts_hits = evt->GetNofData(ECbmDataType::kStsHit);
  if (nb_sts_hits > 0)
    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 CbmStsRecoBeamSpot::ProcessHit(CbmStsHit* hit)
{
  if (hit == nullptr) return;
  int32_t address = hit->GetAddress();
  int32_t unit    = CbmStsAddress::GetElementId(address, kStsUnit);
 
  if (fAnalysisCuts != nullptr && fAnalysisCuts->CheckStsHit(hit, fStsCluArray)) {
    fStsHits[unit].push_back(hit);
  }
}
 
void CbmStsRecoBeamSpot::Exec(Option_t*)
{
 
  // Check if CbmEvent
  if (fCbmEvtArray != nullptr) {
 
    uint32_t nb_events = fCbmEvtArray->GetEntriesFast();
    LOG(info) << "Running CbmStsRecoBeamSpot - Event like - Entry:" << entry_ << "\tEvents: " << nb_events;
    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)) {
        ProcessEvent(event);
        BeamSpotReco();
      }
    }
  }
  else {
 
    uint32_t n_of_hits = fStsHitArray->GetEntriesFast();
    LOG(info) << "Running CbmStsRecoBeamSpot - Time like - Entry:" << entry_ << "\tStsHit: " << 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);
    }
    BeamSpotReco();
  }
  entry_++;
}
 
InitStatus CbmStsRecoBeamSpot::Init()
{
  LOG(debug) << "Init CbmStsRecoBeamSpot ...";
 
  FairRootManager* ioman = FairRootManager::Instance();
  if (ioman == nullptr) return kERROR;
 
  fCbmEvtArray = (TClonesArray*) ioman->GetObject("CbmEvent");
  fStsHitArray = (TClonesArray*) ioman->GetObject("StsHit");
  fStsCluArray = (TClonesArray*) ioman->GetObject("StsCluster");
 
  for (auto& [trg_name, trg] : fTargets) {
    LOG(debug) << trg_name << ": " << trg->ToString();
  }
 
  LoadSetup();
  BookHistograms();
 
  return kSUCCESS;
}