CbmCaMCModule.h

Go to the documentation of this file.

/* Copyright (C) 2022-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Sergei Zharko [committer] */
 
 
#ifndef CbmCaMCModule_h
#define CbmCaMCModule_h 1
 
#include "CaMonitor.h"
#include "CaParameters.h"
#include "CaToolsMCData.h"
#include "CaToolsMCPoint.h"
#include "CaVector.h"
#include "CbmL1DetectorID.h"
#include "CbmL1Hit.h"
#include "CbmL1Track.h"
#include "CbmLink.h"
#include "CbmMCDataArray.h"
#include "CbmMCEventList.h"
#include "CbmMCTrack.h"
#include "CbmMatch.h"
#include "CbmMuchPoint.h"
#include "CbmMuchTrackingInterface.h"
#include "CbmMvdPoint.h"
#include "CbmMvdTrackingInterface.h"
#include "CbmStsPoint.h"
#include "CbmStsTrackingInterface.h"
#include "CbmTimeSlice.h"
#include "CbmTofPoint.h"
#include "CbmTofTrackingInterface.h"
#include "CbmTrdPoint.h"
#include "CbmTrdTrackingInterface.h"
#include "TClonesArray.h"
#include "TDatabasePDG.h"
 
#include <numeric>
#include <optional>
#include <string>
#include <string_view>
#include <type_traits>
 
class CbmEvent;
class CbmMCDataObject;
class CbmL1HitDebugInfo;
class CbmL1Track;
 
namespace cbm::ca
{
  namespace ca = cbm::algo::ca;
 
  class MCModule {
 
 
   public:
    MCModule(int verb = 1, int perfMode = 1) : fVerbose(verb), fPerformanceMode(perfMode)
    {
      LOG(info) << "cbm::ca::MCModule: performance mode = " << fPerformanceMode;
    }
 
    ~MCModule() = default;
 
    MCModule(const MCModule&) = delete;
 
    MCModule(MCModule&&) = delete;
 
    MCModule& operator=(const MCModule&) = delete;
 
    MCModule& operator=(MCModule&&) = delete;
 
 
    void Finish();
 
    const tools::MCData* GetMCData() const { return fpMCData; }
 
    void InitEvent(CbmEvent* pEvent);
 
    bool InitRun();
 
    template<ca::EDetectorID DetId>
    std::tuple<int, std::vector<int>> MatchHitWithMc(int iHitExt);
 
    void MatchHits();
 
    void MatchTracks();
 
    void ProcessEvent(CbmEvent* pEvent);
 
    void RegisterFirstHitIndexes(const std::array<int, constants::size::MaxNdetectors + 1>& source)
    {
      fpvFstHitId = &source;
    }
 
    void SetDetector(ca::EDetectorID detID, bool flag) { fvbUseDet[detID] = flag; }
 
    void RegisterMCData(tools::MCData& mcData) { fpMCData = &mcData; }
 
    void RegisterRecoTrackContainer(ca::Vector<CbmL1Track>& vRecoTracks) { fpvRecoTracks = &vRecoTracks; }
 
    void RegisterHitIndexContainer(ca::Vector<CbmL1HitId>& vHitIds) { fpvHitIds = &vHitIds; }
 
    void RegisterParameters(std::shared_ptr<ca::Parameters<float>>& pParameters) { fpParameters = pParameters; }
 
    void RegisterQaHitContainer(ca::Vector<CbmL1HitDebugInfo>& vQaHits) { fpvQaHits = &vQaHits; }
 
    int GetVerbosity() const { return fVerbose; }
 
    void CheckInit() const;
 
    void InitTrackInfo();
 
    template<ca::EDetectorID DetID>
    void MatchPointsAndHits();
 
    void MatchRecoAndMCTracks();
 
   private:
    void ReadMCTracks();
 
    void ReadMCPoints();
 
    template<ca::EDetectorID DetID>
    void ReadMCPointsForDetector();
 
    template<ca::EDetectorID DetID>
    std::optional<tools::MCPoint> FillMCPoint(int iExtId, int iEvent, int iFile);
 
    enum class EMonitorKey
    {
      kMcTrack,                 
      kMcTrackReconstructable,  
      kMcPoint,                 
      kRecoNevents,             
      kMissedMatchesMvd,        
      kMissedMatchesSts,        
      kMissedMatchesMuch,       
      kMissedMatchesTrd,        
      kMissedMatchesTof,        
      END
    };
    ca::Monitor<EMonitorKey> fMonitor{"CA MC Module"};  
 
    // ------ Flags
    DetIdArr_t<bool> fvbUseDet = {{false}};  
    int fVerbose               = 1;          
    int fPerformanceMode       = -1;         
 
    std::shared_ptr<ca::Parameters<float>> fpParameters = nullptr;  
    // ------ Input data branches
    const CbmTimeSlice* fpTimeSlice  = nullptr;  
    CbmMCEventList* fpMCEventList    = nullptr;  
    CbmMCDataObject* fpMCEventHeader = nullptr;  
    CbmMCDataArray* fpMCTracks       = nullptr;  
 
    DetIdArr_t<CbmTrackingDetectorInterfaceBase*> fvpDetInterface = {{nullptr}};  
 
    DetIdArr_t<CbmMCDataArray*> fvpBrPoints   = {{nullptr}};  
    DetIdArr_t<TClonesArray*> fvpBrHitMatches = {{nullptr}};  
 
    // Matching information
    std::set<std::pair<int, int>> fFileEventIDs;  
    int fBestMcFile  = -1;                        
    int fBestMcEvent = -1;                        
 
 
    // ----- Internal MC data
    tools::MCData* fpMCData = nullptr;  
 
    // ----- Internal reconstructed data
    ca::Vector<CbmL1Track>* fpvRecoTracks    = nullptr;  
    ca::Vector<CbmL1HitId>* fpvHitIds        = nullptr;  
    ca::Vector<CbmL1HitDebugInfo>* fpvQaHits = nullptr;  
 
    const std::array<int, constants::size::MaxNdetectors + 1>* fpvFstHitId = nullptr;
  };
 
 
  // **********************************************
  // **     Template function implementation     **
  // **********************************************
 
  // -------------------------------------------------------------------------------------------------------------------
  //
  template<ca::EDetectorID DetID>
  std::tuple<int, std::vector<int>> MCModule::MatchHitWithMc(int iHitExt)
  {
    int iPoint = -1;
    std::vector<int> vPoints;
    const auto* pHitMatch = dynamic_cast<CbmMatch*>(fvpBrHitMatches[DetID]->At(iHitExt));
    if (!pHitMatch) {
      LOG(warn) << "Hit match with index " << iHitExt << " is missing for " << kDetName[DetID];
      if constexpr (ca::EDetectorID::kMvd == DetID) {
        fMonitor.IncrementCounter(EMonitorKey::kMissedMatchesMvd);
      }
      else if constexpr (ca::EDetectorID::kSts == DetID) {
        fMonitor.IncrementCounter(EMonitorKey::kMissedMatchesSts);
      }
      else if constexpr (ca::EDetectorID::kMuch == DetID) {
        fMonitor.IncrementCounter(EMonitorKey::kMissedMatchesMuch);
      }
      else if constexpr (ca::EDetectorID::kTrd == DetID) {
        fMonitor.IncrementCounter(EMonitorKey::kMissedMatchesTrd);
      }
      else if constexpr (ca::EDetectorID::kTof == DetID) {
        fMonitor.IncrementCounter(EMonitorKey::kMissedMatchesTof);
      }
      return std::tuple(iPoint, vPoints);
    }
 
    for (int iLink = 0; iLink < pHitMatch->GetNofLinks(); ++iLink) {
      const auto& link = pHitMatch->GetLink(iLink);
      int iPointExt    = link.GetIndex();
      int iEvent       = link.GetEntry();
      int iFile        = link.GetFile();
      if (iPointExt < 0) continue;
      int id = fpMCData->FindInternalPointIndex(DetID, iPointExt, iEvent, iFile);
      vPoints.push_back(id);
    }
 
    if (pHitMatch->GetNofLinks() > 0) {
      const auto& link = pHitMatch->GetMatchedLink();
      if (link.GetIndex() > -1) {
        int index = link.GetIndex();
        int event = link.GetEntry();
        int file  = link.GetFile();
        iPoint    = fpMCData->FindInternalPointIndex(DetID, index, event, file);
      }
    }
 
    return std::tuple(iPoint, vPoints);
  }
 
  // -------------------------------------------------------------------------------------------------------------------
  //
  template<ca::EDetectorID DetID>
  std::optional<tools::MCPoint> MCModule::FillMCPoint(int iExtId, int iEvent, int iFile)
  {
    auto oPoint = std::make_optional<tools::MCPoint>();
 
    // ----- Get positions, momenta, time and track ID
    TVector3 posIn;   // Position at entrance to station [cm]
    TVector3 posOut;  // Position at exist of station [cm]
    TVector3 momIn;   // 3-momentum at entrance to station [GeV/c]
    TVector3 momOut;  // 3-momentum at exit of station [GeV/c]
 
    auto* pBrPoints = fvpBrPoints[DetID];
 
    PointTypes_t::at<DetID>* pExtPoint = dynamic_cast<PointTypes_t::at<DetID>*>(pBrPoints->Get(iFile, iEvent, iExtId));
    if (!pExtPoint) {
      LOG(warn) << "CbmCaMCModule: " << kDetName[DetID] << " MC point with iExtId = " << iExtId
                << ", iEvent = " << iEvent << ", iFile = " << iFile << " does not exist";
      return std::nullopt;
    }
    if constexpr (ca::EDetectorID::kMvd == DetID) {
      pExtPoint->Position(posIn);
      pExtPoint->PositionOut(posOut);
      pExtPoint->Momentum(momIn);
      pExtPoint->MomentumOut(momOut);
    }
    // STS
    else if constexpr (ca::EDetectorID::kSts == DetID) {
      pExtPoint->Position(posIn);
      pExtPoint->PositionOut(posOut);
      pExtPoint->Momentum(momIn);
      pExtPoint->MomentumOut(momOut);
    }
    // MuCh
    else if constexpr (ca::EDetectorID::kMuch == DetID) {
      pExtPoint->Position(posIn);
      pExtPoint->PositionOut(posOut);
      pExtPoint->Momentum(momIn);
      pExtPoint->Momentum(momOut);
    }
    // TRD
    else if constexpr (ca::EDetectorID::kTrd == DetID) {
      pExtPoint->Position(posIn);
      pExtPoint->PositionOut(posOut);
      pExtPoint->Momentum(momIn);
      pExtPoint->MomentumOut(momOut);
    }
    // TOF
    else if constexpr (ca::EDetectorID::kTof == DetID) {
      pExtPoint->Position(posIn);
      pExtPoint->Position(posOut);
      pExtPoint->Momentum(momIn);
      pExtPoint->Momentum(momOut);
    }
    double time = pExtPoint->GetTime();
    int iTmcExt = pExtPoint->GetTrackID();
 
    if (iTmcExt < 0) {
      LOG(warn) << "CbmCaMCModule: For MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
                << ", iFile = " << iFile << " MC track is undefined (has ID = " << iTmcExt << ')';
      return std::nullopt;
    }
    TVector3 posMid = 0.5 * (posIn + posOut);
    TVector3 momMid = 0.5 * (momIn + momOut);
 
    // // ----- Get station index
    int iStLoc = fvpDetInterface[DetID]->GetTrackingStationIndex(pExtPoint);
    if (iStLoc < 0) {
      return std::nullopt;
    }
 
    int stationID = fpParameters->GetStationIndexActive(iStLoc, DetID);
    if (stationID == -1) {
      return std::nullopt;
    }  // Skip points from inactive stations
 
    // Update point time with event time
    time += fpMCEventList->GetEventTime(iEvent, iFile);
 
    // ----- Reject MC points falling out of the time slice
    // STS, MuCh, TRD, TOF
    if constexpr (DetID != ca::EDetectorID::kMvd) {
      // TODO: SZh 18.06.2024: Avoid dependency from CbmTimeSlice, if possible (no TimeSlice branch in the online unpack)
      double startT = fpTimeSlice->GetStartTime();
      double endT   = fpTimeSlice->GetEndTime();
 
      if ((startT > 0. && time < startT) || (endT > 0. && time > endT)) {
        LOG(warn) << "CbmCaMCModule: MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
                  << ", iFile = " << iFile << " and det id " << int(DetID) << " fell out of the TS duration [" << startT
                  << ", " << endT << "] with measured time = " << time << " [ns]";
        return std::nullopt;
      }
    }
 
    // ----- Fill MC point
    oPoint->SetExternalId(fpMCData->GetPointGlobExtIndex(DetID, iExtId));
    oPoint->SetEventId(iEvent);
    oPoint->SetFileId(iFile);
    oPoint->SetTime(time);
    oPoint->SetX(posMid.X());
    oPoint->SetY(posMid.Y());
    oPoint->SetZ(posMid.Z());
    oPoint->SetXIn(posIn.X());
    oPoint->SetYIn(posIn.Y());
    oPoint->SetZIn(posIn.Z());
    oPoint->SetXOut(posOut.X());
    oPoint->SetYOut(posOut.Y());
    oPoint->SetZOut(posOut.Z());
    oPoint->SetPx(momMid.X());
    oPoint->SetPy(momMid.Y());
    oPoint->SetPz(momMid.Z());
    oPoint->SetPxIn(momIn.X());
    oPoint->SetPyIn(momIn.Y());
    oPoint->SetPzIn(momIn.Z());
    oPoint->SetPxOut(momOut.X());
    oPoint->SetPyOut(momOut.Y());
    oPoint->SetPzOut(momOut.Z());
 
    // Select current number of points as a local id of points
    oPoint->SetId(fpMCData->GetNofPoints());
 
    // Match MC track and point to each other
    int iTmcInt = fpMCData->FindInternalTrackIndex(iTmcExt, iEvent, iFile);
 
    oPoint->SetTrackId(iTmcInt);
    if (iTmcInt > -1) {
      fpMCData->GetTrack(iTmcInt).AddPointIndex(oPoint->GetId());
    }
 
    oPoint->SetStationId(stationID);
    oPoint->SetDetectorId(DetID);
 
    auto* pExtTrk = dynamic_cast<CbmMCTrack*>(fpMCTracks->Get(iFile, iEvent, iTmcExt));
    if (!pExtTrk) {
      LOG(warn) << "CbmCaMCModule: MC track with iTmcExt = " << iTmcExt << ", iEvent = " << iEvent
                << ", iFile = " << iFile << " MC track is undefined (nullptr)";
      return std::nullopt;
    }
    oPoint->SetPdgCode(pExtTrk->GetPdgCode());
    oPoint->SetMotherId(pExtTrk->GetMotherId());
 
    auto* pPdgDB = TDatabasePDG::Instance()->GetParticle(oPoint->GetPdgCode());
    oPoint->SetMass(pPdgDB ? pPdgDB->Mass() : 0.);  
    oPoint->SetCharge(pPdgDB ? pPdgDB->Charge() / 3. : 0.);
 
    return oPoint;
  }
 
  // -------------------------------------------------------------------------------------------------------------------
  //
  template<ca::EDetectorID DetID>
  void MCModule::MatchPointsAndHits()
  {
    if (!fvbUseDet[DetID]) {
      return;
    }
 
    for (int iH = (*fpvFstHitId)[static_cast<int>(DetID)]; iH < (*fpvFstHitId)[static_cast<int>(DetID) + 1]; ++iH) {
      auto& hit            = (*fpvQaHits)[iH];
      auto [iBestP, vAllP] = MatchHitWithMc<DetID>(hit.ExtIndex);
      if (iBestP >= 0) {
        hit.SetBestMcPointId(iBestP);
      }
      for (auto iP : vAllP) {
        if (iP >= 0) {
          hit.AddMcPointId(iP);
          fpMCData->GetPoint(iP).AddHitID(iH);
        }
      }
    }
  }
  // -------------------------------------------------------------------------------------------------------------------
  // NOTE: template is used, because another template function FillMCPoint is used inside
  template<ca::EDetectorID DetID>
  void MCModule::ReadMCPointsForDetector()
  {
    if (!fvbUseDet[DetID]) {
      return;
    }
    for (const auto& [iFile, iEvent] : fFileEventIDs) {
      int nPointsEvent = fvpBrPoints[DetID]->Size(iFile, iEvent);
      for (int iP = 0; iP < nPointsEvent; ++iP) {
        std::optional<tools::MCPoint> oPoint = FillMCPoint<DetID>(iP, iEvent, iFile);
        if (oPoint) {
          fpMCData->AddPoint(*oPoint);
        }
      }  // iP: end
    }    // key: end
  }
}  // namespace cbm::ca
 
#endif  // CbmCaMCModule_h