CbmEventVertexDca.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 "CbmEventVertexDca.h"
 
#include "CbmStsAddress.h"
 
 
void CbmEventVertexDca::BookHistograms()
{
  std::string h_name;
 
  // ---- Binning configuration ----
  // 3D Vertex
  double vertex_dx    = 0.01;
  double vertex_dy    = 0.01;
  double vertex_dz    = 0.01;
  double vertex_x_min = -35 + 0.5 * vertex_dx;
  double vertex_x_max = +35 + 0.5 * vertex_dx;
  double vertex_y_min = -35 + 0.5 * vertex_dy;
  double vertex_y_max = +35 + 0.5 * vertex_dy;
  double vertex_z_min = -65 + 0.5 * vertex_dz;
  double vertex_z_max = +65 + 0.5 * vertex_dz;
 
  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};
  // ---- --------------------- ----
 
  // ---- 3D Vertex 2D projections ----
  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]");
  // ---- ------------------------ ----
 
  // ---- PCA 2D projections ----
  if (fVertexFinder) {
    h_name = "pca_y_vs_x";
    fH2DShared[h_name] =
      std::make_shared<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);
    fH2DShared[h_name]->GetXaxis()->SetTitle("Vertex_{X} [cm]");
    fH2DShared[h_name]->GetYaxis()->SetTitle("Vertex_{Y} [cm]");
 
    h_name = "pca_x_vs_z";
    fH2DShared[h_name] =
      std::make_shared<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);
    fH2DShared[h_name]->GetXaxis()->SetTitle("Vertex_{Z} [cm]");
    fH2DShared[h_name]->GetYaxis()->SetTitle("Vertex_{X} [cm]");
 
    h_name = "pca_y_vs_z";
    fH2DShared[h_name] =
      std::make_shared<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);
    fH2DShared[h_name]->GetXaxis()->SetTitle("Vertex_{Z} [cm]");
    fH2DShared[h_name]->GetYaxis()->SetTitle("Vertex_{Y} [cm]");
 
    std::cout << " - DcaVertexFinder - Enabling QA" << std::endl;
    fVertexFinder->EnableQa({fH2DShared["pca_y_vs_x"], fH2DShared["pca_x_vs_z"], fH2DShared["pca_y_vs_z"]});
  }
  // ---- ------------------------ ----
 
  // ----        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));
  }
  // ---- ------------------ ----
 
  // ----     Residual with MC      ----
  if (input_has_mc) {
    LOG(info) << " - INFO - Enabling MC comparison";
    for (const char* di : {"x", "y", "z"}) {
      h_name = Form("res_%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("%s_{rec} - %s_{MC} [cm]", di, di));
 
      h_name = Form("res_%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("%s_{rec} - %s_{MC} [cm]", di, di));
 
      h_name = Form("res_%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("%s_{rec} - %s_{MC} [cm]", di, di));
    }
  }
  // ---- ------------------------- ----
 
  // Track multiplicity
  h_name       = "ca_track_multiplicity_all";
  fH1D[h_name] = std::make_unique<TH1D>(h_name.c_str(), h_name.c_str(), 50, 0, 50);
  fH1D[h_name]->GetXaxis()->SetTitle("N_{CATrack}");
  fH1D[h_name]->GetYaxis()->SetTitle("Entries");
 
  h_name       = "ca_track_multiplicity_sel";
  fH1D[h_name] = std::make_unique<TH1D>(h_name.c_str(), h_name.c_str(), 50, 0, 50);
  fH1D[h_name]->GetXaxis()->SetTitle("N_{CATrack}");
  fH1D[h_name]->GetYaxis()->SetTitle("Entries");
}
 
 
void CbmEventVertexDca::FinishTask() { SaveToFile(); }
 
void CbmEventVertexDca::CheckVertex()
{
  fH1D["ca_track_multiplicity_sel"]->Fill(fGlbTrks.size());
 
  fVertexFinder->SetTracks(fGlbTrks);
  const std::optional<CbmVertex> vertex = fVertexFinder->FindVertex();
 
  if (vertex.has_value()) {
 
    TVector3 vtx(vertex->GetX(), vertex->GetY(), vertex->GetZ());
 
    fH2D["vertex_y_vs_x"]->Fill(vtx[0], vtx[1]);
    fH2D["vertex_x_vs_z"]->Fill(vtx[2], vtx[0]);
    fH2D["vertex_y_vs_z"]->Fill(vtx[2], vtx[1]);
 
    for (const CbmGlobalTrack* trk : fGlbTrks) {
      float trk_tx = trk->GetParamFirst()->GetTx();
      float trk_ty = trk->GetParamFirst()->GetTy();
      float trk_x0 = trk->GetParamFirst()->GetX();
      float trk_y0 = trk->GetParamFirst()->GetY();
      float trk_z0 = trk->GetParamFirst()->GetZ();
 
      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 (input_has_mc) {
      TVector3 mc_vertex;
      int n_of_mc_trks = fMCTrkArray->GetEntriesFast();
      for (int mc_trk_idx = 0; mc_trk_idx < n_of_mc_trks; mc_trk_idx++) {
        auto mc_trk = (CbmMCTrack*) fMCTrkArray->UncheckedAt(mc_trk_idx);
        if (mc_trk->GetMotherId() == -1) {
          mc_trk->GetStartVertex(mc_vertex);
          break;
        }
      }
 
      fH2D["res_x_vs_x"]->Fill(vtx[0] - mc_vertex[0], vtx[0]);
      fH2D["res_x_vs_y"]->Fill(vtx[0] - mc_vertex[0], vtx[1]);
      fH2D["res_x_vs_z"]->Fill(vtx[0] - mc_vertex[0], vtx[2]);
 
      fH2D["res_y_vs_x"]->Fill(vtx[1] - mc_vertex[1], vtx[0]);
      fH2D["res_y_vs_y"]->Fill(vtx[1] - mc_vertex[1], vtx[1]);
      fH2D["res_y_vs_z"]->Fill(vtx[1] - mc_vertex[1], vtx[2]);
 
      fH2D["res_z_vs_x"]->Fill(vtx[2] - mc_vertex[2], vtx[0]);
      fH2D["res_z_vs_y"]->Fill(vtx[2] - mc_vertex[2], vtx[1]);
      fH2D["res_z_vs_z"]->Fill(vtx[2] - mc_vertex[2], vtx[2]);
    }
  }
 
  // Clear track containers
  fGlbTrks.clear();
}
 
void CbmEventVertexDca::ProcessEvent(CbmEvent* evt)
{
  int nb_of_glob_trk = evt->GetNofData(ECbmDataType::kGlobalTrack);
  fH1D["ca_track_multiplicity_all"]->Fill(nb_of_glob_trk);
 
  if (fAnalysisCuts != nullptr && !fAnalysisCuts->CheckEvent(evt)) return;
 
 
  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);
    ProcessGlobalTrack((CbmGlobalTrack*) fGlbTrkArray->UncheckedAt(glob_trk_idx));
  }
 
  CheckVertex();
}
 
void CbmEventVertexDca::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;
  }
 
  fGlbTrks.push_back(trk);
}
 
void CbmEventVertexDca::Exec(Option_t*)
{
  LOG(info) << "Running CbmEventVertexDca;\t"
            << "Entry: " << entry_;
 
  // Check if CbmEvent
  if (fCbmEvtArray != nullptr) {
    uint32_t nb_events = fCbmEvtArray->GetEntriesFast();
    LOG(info) << Form(" ... number of CbmEvent: %u\n", nb_events);
    for (uint32_t evt_idx = 0; evt_idx < nb_events; evt_idx++) {
      ProcessEvent((CbmEvent*) fCbmEvtArray->UncheckedAt(evt_idx));
    }
  }
  else {
 
    LOG(debug) << " ... running time-like*** Using all tracks in TS - risk of high combinatorial";
 
    uint32_t nb_of_glob_trk = fGlbTrkArray->GetEntriesFast();
    LOG(info) << Form("Number of GlobalTracks: %u\n", 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);
    }
 
    CheckVertex();
  }
  entry_++;
}
 
InitStatus CbmEventVertexDca::Init()
{
  LOG(debug) << "Init CbmEventVertexDca ...";
 
  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");
 
    fMCTrkArray  = (TClonesArray*) ioman->GetObject("MCTrack");
    input_has_mc = fMCTrkArray != nullptr;
  }
 
  LoadSetup();
 
  BookHistograms();
 
  fGlbTrks.reserve(100);
 
  return kSUCCESS;
}
 
void CbmEventVertexDca::SetVertexFinder(std::shared_ptr<CbmDcaVertexFinder> vtx_finder) { fVertexFinder = vtx_finder; }