CbmRichRingRecoQa.cxx

Go to the documentation of this file.

/* Copyright (C) 2025 Justus-Liebig-Universitaet Giessen, Giessen
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Martin Beyer [committer] */

 
#include "CbmRichRingRecoQa.h"
 
#include "CbmDigiManager.h"
#include "CbmEvent.h"
#include "CbmMCDataArray.h"
#include "CbmMCDataManager.h"
#include "CbmMCTrack.h"
#include "CbmMatch.h"
#include "CbmMatchRecoToMC.h"
#include "CbmQaTable.h"
#include "CbmRichHit.h"
#include "CbmRichRing.h"
#include "CbmTrackMatchNew.h"
 
#include <FairRootManager.h>
 
#include <TClonesArray.h>
#include <TMCProcess.h>
 
#include <array>
#include <vector>
 
 
InitStatus CbmRichRingRecoQa::InitQa()
{
  FairRootManager* manager = FairRootManager::Instance();
  if (!manager) LOG(fatal) << GetName() << ": No FairRootManager!";
 
  fDigiManager = CbmDigiManager::Instance();
  fDigiManager->Init();
  if (!fDigiManager->IsPresent(ECbmModuleId::kRich)) {
    LOG(fatal) << GetName() << ": No RichDigi array found in CbmDigiManager!";
  }
 
  fCbmEvents = dynamic_cast<TClonesArray*>(manager->GetObject("CbmEvent"));
  LOG(info) << GetName() << ": CbmEvent" << (fCbmEvents ? " " : " NOT ") << "found.\n";
 
  fRichHits = dynamic_cast<TClonesArray*>(manager->GetObject("RichHit"));
  if (!fRichHits) LOG(fatal) << GetName() << ": No RichHit array found in FairRootManager!";
 
  fRichRings = dynamic_cast<TClonesArray*>(manager->GetObject("RichRing"));
  if (!fRichRings) LOG(fatal) << GetName() << ": No RichRing array found in FairRootManager!";
 
  fRichProjections = dynamic_cast<TClonesArray*>(manager->GetObject("RichProjection"));
  if (!fRichProjections) LOG(fatal) << GetName() << ": No RichProjection array found in FairRootManager!";
 
 
  if (IsMCUsed()) {
    auto* mcManager = static_cast<CbmMCDataManager*>(manager->GetObject("MCDataManager"));
    if (!mcManager) LOG(fatal) << GetName() << ": No MCDataManager!";
 
    fMcTracks = mcManager->InitBranch("MCTrack");
    if (!fMcTracks) LOG(fatal) << GetName() << ": No MCTrack array found in CbmMCDataManager!";
 
    fRichPoints = mcManager->InitBranch("RichPoint");
    if (!fRichPoints) LOG(fatal) << GetName() << ": No RichPoint array found in CbmMCDataManager!";
 
    fRichRingMatches = dynamic_cast<TClonesArray*>(manager->GetObject("RichRingMatch"));
    if (!fRichRingMatches) LOG(fatal) << GetName() << ": No RichRingMatch array found in FairRootManager!";
  }
 
  InitHistograms();
 
  return kSUCCESS;
}
 
void CbmRichRingRecoQa::InitHistograms()
{
  // clang-format off
  fhRichRingsPerEvent =
    MakeQaObject<TH1D>("RichRingsPerEvent", "RichRings per event; #RichRings / Event; Entries", 71, -0.5, 70.5);
  fhRingTimeInEvent =
    MakeQaObject<TH1D>("RichRingTimeEvent", "Rich ring time in event; Ring time - Event time [ns]; Entires", 30, 12., 24.);
  fhRingHitsTime =
    MakeQaObject<TH1D>("RingHitsTime", "Rich ring time - Hit times; Ring time - Hit time [ns]; Entires", 100, -5., 5.);
  fhRingCenterXY =
    MakeQaObject<TH2D>("RingXY", "Rich ring positions; x [cm]; y [cm]; Entries", 50, -130., 130.0, 20, 110., 200.);
 
  // 1D feature histograms 
  fhRingFeatures.at(ERingFeature::Radius) =
    MakeQaObject<TH1D>("RingRadius", "Rich ring radius; Radius [cm]; Entries", 50, 0., 10.);
  fhRingFeatures.at(ERingFeature::NHits) =
    MakeQaObject<TH1D>("RingNumHits", "Rich ring number of hits; #Ring hits; Entries", 51, -0.5, 50.5);
  fhRingFeatures.at(ERingFeature::CenterX) =
    MakeQaObject<TH1D>("RingCenterX", "Rich ring center x position; Center x [cm]; Entries", 100, -130., 130.0);
  fhRingFeatures.at(ERingFeature::CenterY) =
    MakeQaObject<TH1D>("RingCenterY", "Rich ring center y position; |Center y| [cm]; Entries", 50, 110., 200.);
  fhRingFeatures.at(ERingFeature::Aaxis) =
    MakeQaObject<TH1D>("RingAaxis", "Rich ring A axis; A axis [cm]; Entries", 50, 0., 10.);
  fhRingFeatures.at(ERingFeature::Baxis) =
    MakeQaObject<TH1D>("RingBaxis", "Rich ring B axis; B axis [cm]; Entries", 50, 0., 10.);
  fhRingFeatures.at(ERingFeature::BoA) =
    MakeQaObject<TH1D>("RingBoA", "Rich ring B over A axis; B/A; Entries", 50, 0., 1.1);
  fhRingFeatures.at(ERingFeature::Phi) =
    MakeQaObject<TH1D>("RingPhi", "Rich ring phi; Phi [rad]; Entries", 40, -2., 2.);
  fhRingFeatures.at(ERingFeature::Chi2oNDF) =
    MakeQaObject<TH1D>("RingChi2oNDF", "Rich ring Chi2 over NDF; Chi2 / NDF; Entries", 40, -0.5, 3.);
  // ! Angle histograms current no filled, since angle value is not passed to CbmRichRing objects
  // fhRingFeatures.at(ERingFeature::Angle) =
  //   MakeQaObject<TH1D>("RingAngle", "Rich ring sum of 3 consecutive angles between hits in ring; Angle [rad]; Entries", 50, 0., 3.1415);
  fhRingFeatures.at(ERingFeature::RadialPosition) =
    MakeQaObject<TH1D>("RingRadialPosition", "Rich ring radial position; RadialPosition [cm]; Entries", 50, 0., 150);
  fhRingFeatures.at(ERingFeature::RadialAngle) =
    MakeQaObject<TH1D>("RingRadialAngle", "Rich ring radial angle; RadialAngle [rad]; Entries", 50, 0., 3.1415);
 
  // * MC part starting here
  if (!IsMCUsed()) return;
  fhRingTypeCounts =
    MakeQaObject<TH1D>("RingTypes", "Reconstructed rings for types; Ring types; Entries", 8, -1.5, 6.5);
  const std::array<std::string, 8> tickLabels{"Reco", "TrueReco", "PrimEl", "SecEl", "Pion", "Other", "Fake", "Clone"};
  for (std::size_t i = 0; i < tickLabels.size(); ++i) {
    fhRingTypeCounts->GetXaxis()->SetBinLabel(i + 1, tickLabels[i].c_str());
  }
 
  fhRingRecoPurityXY  =
    MakeQaObject<TProfile2D>(
      "RingRecoPurityXY",
      "Rich ring reconstruction purity; x [cm]; y [cm]; Purity",
      50, -130., 130.0, 20, 110., 200.);
 
  const std::array<TString, 6> typesStr{"TrueReco", "PrimEl", "SecEl", "Pion", "Other", "Fake"};
  const std::array<TString, 6> typesStrFull{"all correct reconstructed rings", "primary electrons",
                                            "seconday electrons", "pions", "others", "fakes"};
  
  // Efficiency vs MC momentum Profile histograms
  for (int type = ERecoRingType::TrueReco; type <= ERecoRingType::Other; type++) {
    fhRingRecoEff.at(type) =
      MakeQaObject<TProfile>(
        "RingRecoEff" + typesStr.at(type),
        "Rich ring reconstruction efficiency" + typesStrFull.at(type) +"; p_{MC} [GeV/c]; Ring reconstruction efficiency",
        60, 0., 12.);
  }
 
  // 1D feature histograms MC
  for (int type = ERecoRingType::TrueReco; type <= ERecoRingType::Fake; type++) {
    auto& featureHists = fhRingTypeFeatures.at(type);
    featureHists.at(ERingFeature::Radius) =
      MakeQaObject<TH1D>(
        "RingFeatures" + typesStr.at(type) + "/RingRadius" + typesStr.at(type),
        "Rich ring radius " + typesStrFull.at(type) + "; Radius [cm]; Entries",
        50, 0., 10.);
    featureHists.at(ERingFeature::NHits) =
      MakeQaObject<TH1D>(
        "RingFeatures" + typesStr.at(type) + "/RingNumHits" + typesStr.at(type),
        "Rich ring number of hits " + typesStrFull.at(type) + "; #Ring hits; Entries",
        51, -0.5, 50.5);
    featureHists.at(ERingFeature::CenterX) =
      MakeQaObject<TH1D>(
        "RingFeatures" + typesStr.at(type) + "/RingCenterX" + typesStr.at(type),
        "Rich ring center x position " + typesStrFull.at(type) + "; Center x [cm]; Entries",
        100, -130., 130.0);
    featureHists.at(ERingFeature::CenterY) =
      MakeQaObject<TH1D>(
        "RingFeatures" + typesStr.at(type) + "/RingCenterY" + typesStr.at(type),
        "Rich ring center y position " + typesStrFull.at(type) + "; |Center y| [cm]; Entries",
        50, 110., 200.);
    featureHists.at(ERingFeature::Aaxis) =
      MakeQaObject<TH1D>(
        "RingFeatures" + typesStr.at(type) + "/RingAaxis" + typesStr.at(type),
        "Rich ring A axis " + typesStrFull.at(type) + "; A axis [cm]; Entries",
        50, 0., 10.);
    featureHists.at(ERingFeature::Baxis) =
      MakeQaObject<TH1D>(
        "RingFeatures" + typesStr.at(type) + "/RingBaxis" + typesStr.at(type),
        "Rich ring B axis " + typesStrFull.at(type) + "; B axis [cm]; Entries",
        50, 0., 10.);
    featureHists.at(ERingFeature::BoA) =
      MakeQaObject<TH1D>(
        "RingFeatures" + typesStr.at(type) + "/RingBoA" + typesStr.at(type),
        "Rich ring B over A axis " + typesStrFull.at(type) + "; B/A; Entries",
        50, 0., 1.1);
    featureHists.at(ERingFeature::Phi) =
      MakeQaObject<TH1D>(
        "RingFeatures" + typesStr.at(type) + "/RingPhi" + typesStr.at(type),
        "Rich ring phi " + typesStrFull.at(type) + "; Phi [rad]; Entries",
        40, -2., 2.);
    featureHists.at(ERingFeature::Chi2oNDF) =
      MakeQaObject<TH1D>(
        "RingFeatures" + typesStr.at(type) + "/RingChi2oNDF" + typesStr.at(type),
        "Rich ring Chi2 over NDF " + typesStrFull.at(type) + "; Chi2 / NDF; Entries",
        40, -0.5, 3.);
    // featureHists.at(ERingFeature::Angle) =
    //   MakeQaObject<TH1D>(
    //     "RingFeatures" + typesStr.at(type) + "/RingAngle" + typesStr.at(type),
    //     "Rich ring sum of 3 consecutive angles between hits in ring " + typesStrFull.at(type) + "; Angle [rad]; Entries",
    //     50, 0., 3.1415);
    featureHists.at(ERingFeature::RadialPosition) =
      MakeQaObject<TH1D>(
        "RingFeatures" + typesStr.at(type) + "/RingRadialPosition" + typesStr.at(type),
        "Rich ring radial position " + typesStrFull.at(type) + "; RadialPosition [cm]; Entries",
        50, 0., 150);
    featureHists.at(ERingFeature::RadialAngle) =
      MakeQaObject<TH1D>(
        "RingFeatures" + typesStr.at(type) + "/RingRadialAngle" + typesStr.at(type),
        "Rich ring radial angle " + typesStrFull.at(type) + "; RadialAngle [rad]; Entries",
        50, 0., 3.1415);
  }
 
  // 2D feature vs MC momentum histograms (no fakes here!)
  for (int type = ERecoRingType::TrueReco; type <= ERecoRingType::Other; type++) {
    auto& featureHists = fhRingTypeFeaturesMomentum.at(type);
    featureHists.at(ERingFeature::Radius) =
      MakeQaObject<TH2D>(
        "RingFeatures" + typesStr.at(type) + "/RingRadiusMomentum" + typesStr.at(type),
        "Rich ring radius vs momentum " + typesStrFull.at(type) + "; p_{MC} [GeV/c]; Radius [cm]; Entries",
        60, 0., 12., 50, 0., 10.);
    featureHists.at(ERingFeature::NHits) =
      MakeQaObject<TH2D>(
        "RingFeatures" + typesStr.at(type) + "/RingNumHitsMomentum" + typesStr.at(type),
        "Rich ring number of hits vs momentum " + typesStrFull.at(type) + "; p_{MC} [GeV/c]; #Ring hits; Entries",
        60, 0., 12., 51, -0.5, 50.5);
    featureHists.at(ERingFeature::CenterX) =
      MakeQaObject<TH2D>(
        "RingFeatures" + typesStr.at(type) + "/RingCenterXMomentum" + typesStr.at(type),
        "Rich ring center x position vs momentum " + typesStrFull.at(type) + "; p_{MC} [GeV/c]; Center x [cm]; Entries",
        60, 0., 12., 100, -130., 130.0);
    featureHists.at(ERingFeature::CenterY) =
      MakeQaObject<TH2D>(
        "RingFeatures" + typesStr.at(type) + "/RingCenterYMomentum" + typesStr.at(type),
        "Rich ring center y position vs momentum " + typesStrFull.at(type) + "; p_{MC} [GeV/c]; |Center y| [cm]; Entries",
        60, 0., 12., 50, 110., 200.);
    featureHists.at(ERingFeature::Aaxis) =
      MakeQaObject<TH2D>(
        "RingFeatures" + typesStr.at(type) + "/RingAaxisMomentum" + typesStr.at(type),
        "Rich ring A axis vs momentum " + typesStrFull.at(type) + "; p_{MC} [GeV/c]; A axis [cm]; Entries",
        60, 0., 12., 50, 0., 10.);
    featureHists.at(ERingFeature::Baxis) =
      MakeQaObject<TH2D>(
        "RingFeatures" + typesStr.at(type) + "/RingBaxisMomentum" + typesStr.at(type),
        "Rich ring B axis vs momentum " + typesStrFull.at(type) + "; p_{MC} [GeV/c]; B axis [cm]; Entries",
        60, 0., 12., 50, 0., 10.);
    featureHists.at(ERingFeature::BoA) =
      MakeQaObject<TH2D>(
        "RingFeatures" + typesStr.at(type) + "/RingBoAMomentum" + typesStr.at(type),
        "Rich ring B over A axis vs momentum " + typesStrFull.at(type) + "; p_{MC} [GeV/c]; B/A; Entries",
        60, 0., 12., 50, 0., 1.1);
    featureHists.at(ERingFeature::Phi) =
      MakeQaObject<TH2D>(
        "RingFeatures" + typesStr.at(type) + "/RingPhiMomentum" + typesStr.at(type),
        "Rich ring phi vs momentum " + typesStrFull.at(type) + "; p_{MC} [GeV/c]; Phi [rad]; Entries",
        60, 0., 12., 40, -2., 2.);
    featureHists.at(ERingFeature::Chi2oNDF) =
      MakeQaObject<TH2D>(
        "RingFeatures" + typesStr.at(type) + "/RingChi2oNDFMomentum" + typesStr.at(type),
        "Rich ring Chi2 over NDF vs momentum " + typesStrFull.at(type) + "; p_{MC} [GeV/c]; Chi2 / NDF; Entries",
        60, 0., 12., 40, -0.5, 3.);
    // featureHists.at(ERingFeature::Angle) =
    //   MakeQaObject<TH2D>(
    //     "RingFeatures" + typesStr.at(type) + "/RingAngleMomentum" + typesStr.at(type),
    //     "Rich ring sum of 3 consecutive angles between hits in ring vs momentum " + typesStrFull.at(type) + "; p_{MC} [GeV/c]; Angle [rad]; Entries",
    //     60, 0., 12., 50, 0., 3.1415);
    featureHists.at(ERingFeature::RadialPosition) =
      MakeQaObject<TH2D>(
        "RingFeatures" + typesStr.at(type) + "/RingRadialPositionMomentum" + typesStr.at(type),
        "Rich ring radial position vs momentum " + typesStrFull.at(type) + "; p_{MC} [GeV/c]; RadialPosition [cm]; Entries",
        60, 0., 12., 50, 0., 150);
    featureHists.at(ERingFeature::RadialAngle) =
      MakeQaObject<TH2D>(
        "RingFeatures" + typesStr.at(type) + "/RingRadialAngleMomentum" + typesStr.at(type),
        "Rich ring radial angle vs momentum " + typesStrFull.at(type) + "; p_{MC} [GeV/c]; RadialAngle [rad]; Entries",
        60, 0., 12., 50, 0., 3.1415);
  }
  // clang-format on
}
 
void CbmRichRingRecoQa::ExecQa()
{
  const CbmEvent* event = GetCurrentEvent();
 
  const int nRichRings = event ? event->GetNofData(ECbmDataType::kRichRing) : fRichRings->GetEntriesFast();
  if (event) fhRichRingsPerEvent->Fill(nRichRings);
  for (int iRing = 0; iRing < nRichRings; ++iRing) {
    const int iRingIndex = event ? event->GetIndex(ECbmDataType::kRichRing, iRing) : iRing;
    const auto* ring     = static_cast<CbmRichRing*>(fRichRings->At(iRingIndex));
    if (!ring) continue;
 
    if (event) fhRingTimeInEvent->Fill(ring->GetTime() - event->GetTzero());
    fhRingCenterXY->Fill(ring->GetCenterX(), ring->GetCenterY());
    const int nHits = ring->GetNofHits();
    for (int iHit = 0; iHit < nHits; ++iHit) {
      const auto* hit = static_cast<CbmRichHit*>(fRichHits->At(ring->GetHit(iHit)));
      if (!hit) continue;
      fhRingHitsTime->Fill(ring->GetTime() - hit->GetTime());
    }
 
    fhRingFeatures.at(ERingFeature::Radius)->Fill(ring->GetRadius());
    fhRingFeatures.at(ERingFeature::NHits)->Fill(ring->GetNofHits());
    fhRingFeatures.at(ERingFeature::CenterX)->Fill(ring->GetCenterX());
    fhRingFeatures.at(ERingFeature::CenterY)->Fill(ring->GetCenterY());
    fhRingFeatures.at(ERingFeature::Aaxis)->Fill(ring->GetAaxis());
    fhRingFeatures.at(ERingFeature::Baxis)->Fill(ring->GetBaxis());
    fhRingFeatures.at(ERingFeature::BoA)->Fill(ring->GetBaxis() / ring->GetAaxis());
    fhRingFeatures.at(ERingFeature::Phi)->Fill(ring->GetPhi());
    fhRingFeatures.at(ERingFeature::Chi2oNDF)->Fill(ring->GetChi2() / static_cast<double>(ring->GetNDF()));
    // fhRingFeatures.at(ERingFeature::Angle)->Fill(ring->GetAngle());
    fhRingFeatures.at(ERingFeature::RadialPosition)->Fill(ring->GetRadialPosition());
    fhRingFeatures.at(ERingFeature::RadialAngle)->Fill(ring->GetRadialAngle());
  }
 
  // * MC part starting here
  if (!IsMCUsed()) return;
 
  // ! Accepted rings are searched for in each camera seperately.
  //   Reason: It happens that Cherenkov cones are split due to the 2 segment mirror.
  // Accepted rings for upper [0] and lower [1] RICH camera
  std::array<std::map<CbmLink, int>, 2> accRingsCam{};
 
  // Count RichHits for each Cherenkov mother particle
  const int nRichHits = event ? event->GetNofData(ECbmDataType::kRichHit) : fRichHits->GetEntriesFast();
  for (int iHit = 0; iHit < nRichHits; iHit++) {
    const int iHitIndex = event ? event->GetIndex(ECbmDataType::kRichHit, iHit) : iHit;
    const auto* hit     = static_cast<CbmRichHit*>(fRichHits->At(iHitIndex));
    if (!hit) continue;
    const std::vector<CbmLink> motherIds =
      CbmMatchRecoToMC::GetMcTrackMotherIdsForRichHit(fDigiManager, hit, fRichPoints, fMcTracks);
    const bool isUpperCamera = hit->GetY() > 0.;
    for (const auto& motherId : motherIds)
      (isUpperCamera ? accRingsCam[0] : accRingsCam[1])[motherId]++;
  }
 
  // Remove mother particles below minimum amount of hits for accepted ring.
  // Set counters to 0 and reuse them later to count how many reco rings match with a specific acc ring.
  for (auto& accRings : accRingsCam) {
    for (auto it = accRings.begin(); it != accRings.end();) {
      auto& [link, nHits] = *it;
      if (nHits < fAccRingMinNumHits) {
        it = accRings.erase(it);
        continue;
      }
      const auto* mcTrack = static_cast<const CbmMCTrack*>(fMcTracks->Get(link));
      if (!mcTrack) {
        it = accRings.erase(it);
        continue;
      }
      // TODO: implement FillHistsAccRing(mcTrack), e.g. count split Cherenkov cones
      nHits = 0;
      ++it;
    }
  }
 
  // Loop reconstruced rings, fill histograms for specific ring types (E.g. primary, pion, fake, etc.)
  fhRingTypeCounts->Fill(-1., static_cast<double>(nRichRings));
  for (Int_t iRing = 0; iRing < nRichRings; ++iRing) {
    const int iRingIndex = event ? event->GetIndex(ECbmDataType::kRichRing, iRing) : iRing;
    const auto* ring     = static_cast<CbmRichRing*>(fRichRings->At(iRingIndex));
    if (!ring) continue;
 
    const auto* ringMatch = static_cast<CbmTrackMatchNew*>(fRichRingMatches->At(iRingIndex));
    if (!ringMatch) {
      LOG(warn) << GetName() << ": ringMatch not found for ring index " << iRingIndex << " !";
      FillHistsFakeRing(ring);
      continue;
    }
 
    // Full digitizer noise ring
    if (ringMatch->GetNofLinks() <= 0) {
      LOG(warn) << GetName() << ": Full digitizer noise ring for index " << iRingIndex
                << " , this should be extremely rare!";
      FillHistsFakeRing(ring);
      continue;
    }
 
    const CbmLink matchedLink{ringMatch->GetMatchedLink()};
    const auto* mcTrack = static_cast<const CbmMCTrack*>(fMcTracks->Get(matchedLink));
    if (!mcTrack) {
      LOG(warn) << GetName() << ": mcTrack not found for matched link " << matchedLink.ToString() << " !";
      FillHistsFakeRing(ring);
      continue;
    }
 
    // Normal fake rings
    const double quota{ringMatch->GetTrueOverAllHitsRatio()};
    if (quota < fQuotaRich) {
      FillHistsFakeRing(ring);
      continue;
    }
 
    const bool isUpperCamera{ring->GetCenterY() > 0.};
    std::map<CbmLink, int>& accRings = isUpperCamera ? accRingsCam[0] : accRingsCam[1];
    if (accRings.count(matchedLink) > 0) {
      accRings[matchedLink]++;
      if (accRings[matchedLink] > 1) {
        fhRingTypeCounts->Fill(static_cast<double>(ERecoRingType::Clone));
      }
      else {
        FillHistsTrueRecoRing(ring, mcTrack);
      }
    }
    else {
      // True reco, but no accepted, example for this case:
      // 5 Hits from one MC mother particle. A ring is reconstruced
      // from those 5 Hits + 2 other. The quota will be above threshold,
      // hence a true ring. But no acc ring counterpart exists,
      // because we require >= 7 hits for a accepted ring usually.
      FillHistsFakeRing(ring);
      continue;
    }
  }
 
  for (const auto& accRings : accRingsCam) {
    for (const auto& [link, nMatchedRings] : accRings) {
      const auto* mcTrack = static_cast<const CbmMCTrack*>(fMcTracks->Get(link));
      if (!mcTrack) continue;
 
      fhRingRecoEff.at(ERecoRingType::TrueReco)->Fill(mcTrack->GetP(), static_cast<double>(nMatchedRings > 0), 1.);
 
      if (IsMcPrimaryElectron(mcTrack)) {
        fhRingRecoEff.at(ERecoRingType::PrimEl)->Fill(mcTrack->GetP(), static_cast<double>(nMatchedRings > 0), 1.);
      }
      else if (IsMcSecondaryElectron(mcTrack)) {
        fhRingRecoEff.at(ERecoRingType::SecEl)->Fill(mcTrack->GetP(), static_cast<double>(nMatchedRings > 0), 1.);
      }
      else if (IsMcPion(mcTrack)) {
        fhRingRecoEff.at(ERecoRingType::Pion)->Fill(mcTrack->GetP(), static_cast<double>(nMatchedRings > 0), 1.);
      }
      else {
        fhRingRecoEff.at(ERecoRingType::Other)->Fill(mcTrack->GetP(), static_cast<double>(nMatchedRings > 0), 1.);
      }
    }
  }
}
 
void CbmRichRingRecoQa::Check()
{
  // * MC part starting here
  if (!IsMCUsed()) return;
  // Sanity check for amount of ringTypes
  const double nRecoRings{fhRingTypeCounts->GetBinContent(ERecoRingType::Reco + 2)};
  const double nTrueRecoRings{fhRingTypeCounts->GetBinContent(ERecoRingType::TrueReco + 2)};
  const double nSumParticleTypes{fhRingTypeCounts->GetBinContent(ERecoRingType::PrimEl + 2)
                                 + fhRingTypeCounts->GetBinContent(ERecoRingType::SecEl + 2)
                                 + fhRingTypeCounts->GetBinContent(ERecoRingType::Pion + 2)
                                 + fhRingTypeCounts->GetBinContent(ERecoRingType::Other + 2)};
  const double nSumFakeClone{fhRingTypeCounts->GetBinContent(ERecoRingType::Fake + 2)
                             + fhRingTypeCounts->GetBinContent(ERecoRingType::Clone + 2)};
  if (static_cast<int>(std::round(nTrueRecoRings)) != static_cast<int>(std::round(nSumParticleTypes))) {
    LOG(error) << "nTrueRecoRings != nSumParticleTypes. This should never happen!";
  }
  if (static_cast<int>(std::round(nRecoRings)) != static_cast<int>(std::round(nTrueRecoRings + nSumFakeClone))) {
    LOG(error) << "nRecoRingsTotal != nTrueRecoRingsTotal + nSumFakeClone. This should never happen!";
  }
 
  const std::array<std::string, 4> flagNames{"Mean efficiency all", "Mean efficiency primary electrons",
                                             "Mean efficiency secondary electrons", "Mean efficiency pions"};
  for (int type = ERecoRingType::TrueReco; type <= ERecoRingType::Pion; type++) {
    double val      = fhRingRecoEff.at(type)->GetMean(2);
    bool status     = val >= fRingRecoEffsMin.at(type);
    std::string msg = status ? "" : Form("value %f < %f", val, fRingRecoEffsMin.at(type));
    StoreCheckResult(flagNames.at(type), status, msg);
  }
  {
    double purity = fhRingTypeCounts->GetBinContent(ERecoRingType::TrueReco + 2)
                    / fhRingTypeCounts->GetBinContent(ERecoRingType::Reco + 2);
    bool status     = purity >= fRingPurityMin;
    std::string msg = status ? "" : Form("value %f < %f", purity, fRingPurityMin);
    StoreCheckResult("Mean ring reco purity", status, msg);
  }
 
  // Efficiency and purity summary
  auto GetErrorFraction = [&](const double meanNom, const double errNom, const double meanDenom,
                              const double errDenom) {
    const double sqrtFactor =
      std::sqrt((errNom / meanNom) * (errNom / meanNom) + (errDenom / meanDenom) * (errDenom / meanDenom));
    return (meanNom / meanDenom) * sqrtFactor;
  };
 
  auto* pTable = MakeQaObject<CbmQaTable>("table_efficiency", "RICH ring reconstruction efficiencies and purity", 8, 1);
  pTable->SetColWidth(30);
  pTable->SetNamesOfCols({" mean +/- error"});
  pTable->SetNamesOfRows({{"All particle efficiency"},
                          {"Primary electron efficiency"},
                          {"Secondary electron efficiency"},
                          {"Pion efficiency"},
                          {"Other efficiency"},
                          {"Fake rate"},
                          {"Clone rate"},
                          {"Purity"}});
  for (int type = ERecoRingType::TrueReco; type <= ERecoRingType::Other; type++) {
    pTable->SetCell(type, 0, fhRingRecoEff.at(type)->GetMean(2), fhRingRecoEff.at(type)->GetMeanError(2));
  }
  const std::array<ERecoRingType, 3> ringTypes{ERecoRingType::Fake, ERecoRingType::Clone, ERecoRingType::TrueReco};
  const std::array<int, 3> rowIndex{5, 6, 7};
  for (int i = 0; i < 3; ++i) {
    pTable->SetCell(rowIndex.at(i), 0,
                    fhRingTypeCounts->GetBinContent(ringTypes.at(i) + 2)
                      / fhRingTypeCounts->GetBinContent(ERecoRingType::Reco + 2),
                    GetErrorFraction(fhRingTypeCounts->GetBinContent(ringTypes.at(i) + 2),
                                     fhRingTypeCounts->GetBinError(ringTypes.at(i) + 2),
                                     fhRingTypeCounts->GetBinContent(ERecoRingType::Reco + 2),
                                     fhRingTypeCounts->GetBinError(ERecoRingType::Reco + 2)));
  }
  LOG(info) << '\n' << pTable->ToString(4, true);
}
 
void CbmRichRingRecoQa::FillHistsTrueRecoRing(const CbmRichRing* ring, const CbmMCTrack* mcTrack)
{
  fhRingTypeCounts->Fill(static_cast<double>(ERecoRingType::TrueReco));
  FillHistsRingFeatures(ring, mcTrack, ERecoRingType::TrueReco);
 
  if (IsMcPrimaryElectron(mcTrack)) {
    fhRingTypeCounts->Fill(static_cast<double>(ERecoRingType::PrimEl));
    FillHistsRingFeatures(ring, mcTrack, ERecoRingType::PrimEl);
  }
  else if (IsMcSecondaryElectron(mcTrack)) {
    fhRingTypeCounts->Fill(static_cast<double>(ERecoRingType::SecEl));
    FillHistsRingFeatures(ring, mcTrack, ERecoRingType::SecEl);
  }
  else if (IsMcPion(mcTrack)) {
    fhRingTypeCounts->Fill(static_cast<double>(ERecoRingType::Pion));
    FillHistsRingFeatures(ring, mcTrack, ERecoRingType::Pion);
  }
  else {
    fhRingTypeCounts->Fill(static_cast<double>(ERecoRingType::Other));
    FillHistsRingFeatures(ring, mcTrack, ERecoRingType::Other);
  }
}
 
void CbmRichRingRecoQa::FillHistsFakeRing(const CbmRichRing* ring)
{
  fhRingTypeCounts->Fill(static_cast<double>(ERecoRingType::Fake));
  FillHistsRingFeatures(ring, nullptr, ERecoRingType::Fake);
}
 
void CbmRichRingRecoQa::FillHistsRingFeatures(const CbmRichRing* ring, const CbmMCTrack* mcTrack,
                                              const ERecoRingType ringType)
{
  fhRingRecoPurityXY->Fill(ring->GetCenterX(), std::abs(ring->GetCenterY()), ringType == ERecoRingType::Fake ? 0. : 1.);
 
  // Fill 1D ring feature histograms
  auto& featureHists = fhRingTypeFeatures.at(ringType);
  featureHists.at(ERingFeature::Radius)->Fill(ring->GetRadius());
  featureHists.at(ERingFeature::NHits)->Fill(ring->GetNofHits());
  featureHists.at(ERingFeature::CenterX)->Fill(ring->GetCenterX());
  featureHists.at(ERingFeature::CenterY)->Fill(ring->GetCenterY());
  featureHists.at(ERingFeature::Aaxis)->Fill(ring->GetAaxis());
  featureHists.at(ERingFeature::Baxis)->Fill(ring->GetBaxis());
  featureHists.at(ERingFeature::BoA)->Fill(ring->GetBaxis() / ring->GetAaxis());
  featureHists.at(ERingFeature::Phi)->Fill(ring->GetPhi());
  featureHists.at(ERingFeature::Chi2oNDF)->Fill(ring->GetChi2() / static_cast<double>(ring->GetNDF()));
  // featureHists.at(ERingFeature::Angle)->Fill(ring->GetAngle());
  featureHists.at(ERingFeature::RadialPosition)->Fill(ring->GetRadialPosition());
  featureHists.at(ERingFeature::RadialAngle)->Fill(ring->GetRadialAngle());
 
  // Fill 2D ring feature vs MC momentum histograms
  if (ringType == ERecoRingType::Fake || !mcTrack) return;  // Skip ring feature vs momentum for fake rings
  auto& featureHistsMomentum = fhRingTypeFeaturesMomentum.at(ringType);
  featureHistsMomentum.at(ERingFeature::Radius)->Fill(mcTrack->GetP(), ring->GetRadius());
  featureHistsMomentum.at(ERingFeature::NHits)->Fill(mcTrack->GetP(), ring->GetNofHits());
  featureHistsMomentum.at(ERingFeature::CenterX)->Fill(mcTrack->GetP(), ring->GetCenterX());
  featureHistsMomentum.at(ERingFeature::CenterY)->Fill(mcTrack->GetP(), ring->GetCenterY());
  featureHistsMomentum.at(ERingFeature::Aaxis)->Fill(mcTrack->GetP(), ring->GetAaxis());
  featureHistsMomentum.at(ERingFeature::Baxis)->Fill(mcTrack->GetP(), ring->GetBaxis());
  featureHistsMomentum.at(ERingFeature::BoA)->Fill(mcTrack->GetP(), ring->GetBaxis() / ring->GetAaxis());
  featureHistsMomentum.at(ERingFeature::Phi)->Fill(mcTrack->GetP(), ring->GetPhi());
  featureHistsMomentum.at(ERingFeature::Chi2oNDF)
    ->Fill(mcTrack->GetP(), ring->GetChi2() / static_cast<double>(ring->GetNDF()));
  // featureHistsMomentum.at(ERingFeature::Angle)->Fill(mcTrack->GetP(), ring->GetAngle());
  featureHistsMomentum.at(ERingFeature::RadialPosition)->Fill(mcTrack->GetP(), ring->GetRadialPosition());
  featureHistsMomentum.at(ERingFeature::RadialAngle)->Fill(mcTrack->GetP(), ring->GetRadialAngle());
}
 
 
bool CbmRichRingRecoQa::IsMcPrimaryElectron(const CbmMCTrack* mcTrack)
{
  if (!mcTrack) return false;
  return (std::abs(mcTrack->GetPdgCode()) == 11
          && (mcTrack->GetMotherId() == -1 || mcTrack->GetGeantProcessId() == kPPrimary));
}
 
bool CbmRichRingRecoQa::IsMcSecondaryElectron(const CbmMCTrack* mcTrack)
{
  if (!mcTrack) return false;
  return (std::abs(mcTrack->GetPdgCode()) == 11 && !IsMcPrimaryElectron(mcTrack));
}
 
bool CbmRichRingRecoQa::IsMcPion(const CbmMCTrack* mcTrack)
{
  if (!mcTrack) return false;
  return std::abs(mcTrack->GetPdgCode()) == 211;
}