HalCbmV0Builder.cxx

Go to the documentation of this file.

/* Copyright (C) 2025-2025 Warsaw University of Technology, Warsaw
    SPDX-License-Identifier: GPL-3.0-only
    Authors: Daniel Wielanek [committer] */
#include "HalCbmV0Builder.h"
 
#include "HalCbmAnalysisManager.h"
#include "HalCbmDetectorID.h"
#include "HalCbmTrack.h"
#include "HalCbmV0TempTrack.h"
#include "HalCbmV0Track.h"
 
#include <TClonesArray.h>
#include <TDatabasePDG.h>
#include <TParticlePDG.h>
 
#include <Hal/Cout.h>
#include <Hal/DataManager.h>
 
namespace Hal
{
 
} /* namespace Hal */
 
HalCbmV0Builder::v0rawInfo HalCbmV0Builder::Convert(const AnalysisTree::Particle& p)
{
  HalCbmV0Builder::v0rawInfo v0Info;
  v0Info.chi2geo          = p.GetField<float>(fV0Container->GetFieldIds().chi2geo);
  v0Info.chi2_prim_first  = p.GetField<float>(fV0Container->GetFieldIds().chi2_prim_first);
  v0Info.chi2_prim_second = p.GetField<float>(fV0Container->GetFieldIds().chi2_prim_second);
  v0Info.chi2_topo        = p.GetField<float>(fV0Container->GetFieldIds().chi2_topo);
  v0Info.cosine_first     = p.GetField<float>(fV0Container->GetFieldIds().cosine_first);
  v0Info.cosine_second    = p.GetField<float>(fV0Container->GetFieldIds().cosine_second);
  v0Info.cosine_topo      = p.GetField<float>(fV0Container->GetFieldIds().cosine_topo);
  v0Info.distance         = p.GetField<float>(fV0Container->GetFieldIds().distance);
  v0Info.l                = p.GetField<float>(fV0Container->GetFieldIds().l);
  v0Info.l_dl             = p.GetField<float>(fV0Container->GetFieldIds().l_over_dl);
  v0Info.mass             = p.GetField<float>(fV0Container->GetFieldIds().mass);
  v0Info.px               = p.GetField<float>(fV0Container->GetFieldIds().px);
  v0Info.py               = p.GetField<float>(fV0Container->GetFieldIds().py);
  v0Info.pz               = p.GetField<float>(fV0Container->GetFieldIds().pz);
  v0Info.x                = p.GetField<float>(fV0Container->GetFieldIds().x);
  v0Info.y                = p.GetField<float>(fV0Container->GetFieldIds().y);
  v0Info.z                = p.GetField<float>(fV0Container->GetFieldIds().z);
  v0Info.pid              = p.GetField<int>(fV0Container->GetFieldIds().pdg);
  v0Info.dau1id           = p.GetField<int>(fV0Container->GetFieldIds().dau1id);
  v0Info.dau2id           = p.GetField<int>(fV0Container->GetFieldIds().dau2id);
 
#ifdef DEBUG_EXTR
  v0Info.dau1px = p.GetField<float>(fV0Container->GetFieldIds().dau1px);
  v0Info.dau1py = p.GetField<float>(fV0Container->GetFieldIds().dau1py);
  v0Info.dau1pz = p.GetField<float>(fV0Container->GetFieldIds().dau1pz);
  v0Info.dau1x  = p.GetField<float>(fV0Container->GetFieldIds().dau1x);
  v0Info.dau1y  = p.GetField<float>(fV0Container->GetFieldIds().dau1y);
  v0Info.dau1z  = p.GetField<float>(fV0Container->GetFieldIds().dau1z);
 
  v0Info.dau2px = p.GetField<float>(fV0Container->GetFieldIds().dau2px);
  v0Info.dau2py = p.GetField<float>(fV0Container->GetFieldIds().dau2py);
  v0Info.dau2pz = p.GetField<float>(fV0Container->GetFieldIds().dau2pz);
  v0Info.dau2x  = p.GetField<float>(fV0Container->GetFieldIds().dau2x);
  v0Info.dau2y  = p.GetField<float>(fV0Container->GetFieldIds().dau2y);
  v0Info.dau2z  = p.GetField<float>(fV0Container->GetFieldIds().dau2z);
#endif
  return v0Info;
}
 
HalCbmV0Builder::HalCbmV0Builder(Int_t motherPid, Int_t posdau, Int_t negdau)
  : fPidMom(motherPid)
  , fPidPosDau(posdau)
  , fPidNegDau(negdau)
{
}
 
HalCbmV0Builder::EInitFlag HalCbmV0Builder::Init()
{
  fSimContainer  = (CbmAnaTreeMcContainer*) Hal::DataManager::Instance()->GetObject(HalCbm::GetContainerName("mc"));
  fRecoContainer = (CbmAnaTreeRecoContainer*) Hal::DataManager::Instance()->GetObject(HalCbm::GetContainerName("reco"));
  fV0Container   = (CbmAnaTreeV0Container*) Hal::DataManager::Instance()->GetObject(HalCbm::GetContainerName("v0"));
  auto trk       = TDatabasePDG::Instance()->GetParticle(fPidMom);
  if (!trk) {
    Hal::Cout::PrintInfo(Form("HalCbmV0Builder::Init cannot find selected mother PID=%i", fPidMom), Hal::EInfo::kError);
    return EInitFlag::kERROR;
  }
  fMass = trk->Mass();
  if (!fV0Container) {
    Hal::Cout::PrintInfo("HalCbmV0Builder::Init cannot find V0 data", Hal::EInfo::kError);
    return EInitFlag::kERROR;
  }
  if (!fRecoContainer) {
    Hal::Cout::PrintInfo("HalCbmV0Builder::Init cannot find reco data", Hal::EInfo::kError);
    return EInitFlag::kERROR;
  }
  fV0Temp = (TClonesArray*) Hal::DataManager::Instance()->GetObject("TempV0.");
  if (!fV0Temp) {
    fV0Temp = new TClonesArray("HalCbmV0TempTrack");
    Hal::DataManager::Instance()->Register("TempV0.", "TempV0", fV0Temp, kFALSE);
    fOwner = kTRUE;
  }
  return EInitFlag::kSUCCESS;
}
 
void HalCbmV0Builder::Exec(Option_t* /*option*/)
{
  auto recoCont    = fRecoContainer->GetVtxTracks();
  auto recoMatch   = fRecoContainer->GetVtx2Sim();
  auto sim         = fSimContainer->GetParticles();
  int charge_field = fRecoContainer->GetFieldIds().vtx_q;
  int nv0s         = fV0Container->GetParticles()->GetNumberOfChannels();
  int mother_field = fSimContainer->GetFieldIds().motherId;
  int ent          = 0;
  if (fOwner) {
    fV0Temp->Clear("");
    fV0Temp->ExpandCreateFast(nv0s);
  }
  else {
    ent = fV0Temp->GetEntriesFast();
    fV0Temp->ExpandCreateFast(ent + nv0s);
  }
  for (int i = 0; i < nv0s; i++) {
    auto particleV0            = fV0Container->GetParticles()->GetChannel(i);
    auto v0struct              = Convert(particleV0);
    auto dau1                  = recoCont->GetChannel(v0struct.dau1id);
    auto dau2                  = recoCont->GetChannel(v0struct.dau2id);
    int ch1                    = dau1.GetField<int>(charge_field);
    int ch2                    = dau2.GetField<int>(charge_field);
    AnalysisTree::Track daupos = dau1, dauneg = dau2;
    if (ch1 < 0 && ch2 > 0) {  // daughters swap
      daupos          = dau2;
      dauneg          = dau1;
      int temp        = v0struct.dau1id;
      v0struct.dau1id = v0struct.dau2id;
      v0struct.dau2id = temp;
    }
    auto v0Track = (HalCbmV0TempTrack*) fV0Temp->ConstructedAt(ent + i);
    v0Track->GetTrack().SetVertexChi2(v0struct.chi2_topo);
    double e = TMath::Sqrt(v0struct.px * v0struct.px + v0struct.py * v0struct.py + v0struct.pz * v0struct.pz
                           + v0struct.mass * v0struct.mass);
    v0Track->GetTrack().SetMomentum(v0struct.px, v0struct.py, v0struct.pz, e);
    v0Track->GetTrack().SetDCA(v0struct.x, v0struct.y, v0struct.z);
    int field_hit = fRecoContainer->GetFieldIds().vtx_mvdhits;
    int field_sts = fRecoContainer->GetFieldIds().vtx_nhits;
    v0Track->GetTrack().SetNMvdHits(daupos.GetField<int>(field_hit) + dauneg.GetField<int>(field_hit));
    v0Track->GetTrack().SetNStsHits(daupos.GetField<int>(field_sts) + dauneg.GetField<int>(field_sts)
                                    - daupos.GetField<int>(field_hit) - dauneg.GetField<int>(field_hit));
    v0Track->GetV0().SetPdgDaughters(fPidPosDau, fPidNegDau);
    v0Track->GetV0().SetPosId(v0struct.dau1id);
    v0Track->GetV0().SetNegId(v0struct.dau2id);
    v0Track->GetV0().SetDecLenght(v0struct.l);
    v0Track->GetV0().SetDecDl(v0struct.l_dl);
    v0Track->GetV0().SetDauDist(v0struct.distance);
 
#ifdef DEBUG_EXTR
    v0Track->GetV0().SetMomentumPos(TVector3(v0struct.dau1px, v0struct.dau1py, v0struct.dau1pz));
    v0Track->GetV0().SetMomentumNeg(TVector3(v0struct.dau2px, v0struct.dau2py, v0struct.dau2pz));
#else
    v0Track->GetV0().SetMomentumPos(TVector3(daupos.GetPx(), daupos.GetPy(), daupos.GetPz()));
    v0Track->GetV0().SetMomentumNeg(TVector3(dauneg.GetPx(), dauneg.GetPy(), dauneg.GetPz()));
#endif
    v0Track->GetV0().SetMcId(v0struct.pid);
    v0Track->SetMcId(-1);
    //mc matching
    if (!recoMatch) continue;
    int mc1 = recoMatch->GetMatchDirect(v0struct.dau1id);
    int mc2 = recoMatch->GetMatchDirect(v0struct.dau2id);
    if (fBeFriendly) {
      if (mc1 >= 0 && mc2 >= 0) {
        auto mctrack1 = sim->GetChannel(mc1);
        auto mctrack2 = sim->GetChannel(mc2);
        int mom1      = mctrack1.GetField<int>(mother_field);
        v0Track->SetMcId(mom1);
      }
      else if (mc1 >= 0) {
        auto mctrack = sim->GetChannel(mc1);
        v0Track->SetMcId(mctrack.GetField<int>(mother_field));
      }
      else if (mc2 >= 0) {
        auto mctrack = sim->GetChannel(mc2);
        v0Track->SetMcId(mctrack.GetField<int>(mother_field));
      }
    }
    else {
      if (mc1 < 0) continue;
      if (mc2 < 0) continue;
      auto mctrack1 = sim->GetChannel(mc1);
      auto mctrack2 = sim->GetChannel(mc2);
      int mom1      = mctrack1.GetField<int>(mother_field);
      int mom2      = mctrack2.GetField<int>(mother_field);
      if (mom1 == mom2) {
        v0Track->SetMcId(mom1);
      }
    }
  }
}
 
HalCbmV0Builder::~HalCbmV0Builder() {}