From 73c04eff98a08ecb04944459d5e34b623ea2bcc3 Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Wed, 26 Apr 2023 14:14:51 +0200
Subject: [PATCH] CA: added input QA to the mCBM macro

---
 macro/mcbm/mcbm_qa.C            | 53 +++++++++++++++++++++++++++++----
 macro/run/run_qa.C              |  2 --
 reco/L1/qa/CbmCaInputQaMuch.cxx | 26 ++++++++--------
 reco/L1/qa/CbmCaInputQaSts.cxx  | 26 ++++++++--------
 reco/L1/qa/CbmCaInputQaTof.cxx  | 28 ++++++++---------
 reco/L1/qa/CbmCaInputQaTrd.cxx  | 30 ++++++++-----------
 6 files changed, 98 insertions(+), 67 deletions(-)

diff --git a/macro/mcbm/mcbm_qa.C b/macro/mcbm/mcbm_qa.C
index 310a3e3b62..6c71367ca9 100644
--- a/macro/mcbm/mcbm_qa.C
+++ b/macro/mcbm/mcbm_qa.C
@@ -175,21 +175,62 @@ void mcbm_qa(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_2020_03_test",
   }
   // ------------------------------------------------------------------------
 
-  // ----- MUCH QA  ---------------------------------
-  if (setup->IsActive(ECbmModuleId::kMuch)) {
+  // ----- Tracking detector interface --------------------------------------
+  run->AddTask(new CbmTrackingDetectorInterfaceInit());
+  // ------------------------------------------------------------------------
+
+  // NOTE (FIXME) SZh, 26.04.2023
+  // For some reason QA tasks interfere one with another. It leads to the
+  // segmentational violation in CbmMuchHitFinderQa::DrawCanvases() function,
+  // if this task is called after the CbmCaInputQaSts task. The problem dis-
+  // appears, if the CbmMuchHitFinderQa task runs before the CbmCaInputQaSts
+  // task.
+
+  // ----- MUCH QA  ---------------------------------------------------------
+  if (bUseMuch) {
     run->AddTask(new CbmMuchTransportQa());
     run->AddTask(new CbmMuchDigitizerQa());
+
+
     CbmMuchHitFinderQa* muchHitFinderQa = new CbmMuchHitFinderQa();
     muchHitFinderQa->SetGeoFileName(muchParFile);
-    run->AddTask(muchHitFinderQa);
+    //run->AddTask(muchHitFinderQa);
+
+    // CA Input QA
+    auto* pCaInputMuch = new CbmCaInputQaMuch(verbose, bUseMC);
+    pCaInputMuch->SetEfficiencyThrsh(0.5, 0, 100);
+    run->AddTask(pCaInputMuch);
   }
   // ------------------------------------------------------------------------
 
+  // ----- STS QA -----------------------------------------------------------
+  if (bUseSts) {
+    // CA Input QA
+    auto* pCaInputSts = new CbmCaInputQaSts(verbose, bUseMC);
+    pCaInputSts->SetEfficiencyThrsh(0.5, 0, 100);
+    run->AddTask(pCaInputSts);
+  }
+  // ------------------------------------------------------------------------
 
-  // ----- CA tracking QA ---------------------------------------------------
-  // Tracking detector interface initialization
-  run->AddTask(new CbmTrackingDetectorInterfaceInit());
+  // ----- TRD QA -----------------------------------------------------------
+  if (bUseTrd) {
+    // CA Input QA
+    auto* pCaInputTrd = new CbmCaInputQaTrd(verbose, bUseMC);
+    pCaInputTrd->SetEfficiencyThrsh(0.5, 0, 100);
+    run->AddTask(pCaInputTrd);
+  }
+  // ------------------------------------------------------------------------
 
+  // ----- TOF QA -----------------------------------------------------------
+  if (bUseTof) {
+    // CA Input QA
+    auto* pCaInputTof = new CbmCaInputQaTof(verbose, bUseMC);
+    pCaInputTof->SetEfficiencyThrsh(0.5, 0, 100);
+    run->AddTask(pCaInputTof);
+  }
+  // ------------------------------------------------------------------------
+
+  // ----- CA tracking QA ---------------------------------------------------
   // Kalman Filter (currently needed to access the magnetic filed, to be
   // removed soon)
   run->AddTask(new CbmKF());
diff --git a/macro/run/run_qa.C b/macro/run/run_qa.C
index aa3b57db2d..93afc8071e 100644
--- a/macro/run/run_qa.C
+++ b/macro/run/run_qa.C
@@ -217,7 +217,6 @@ void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "d
     }
     run->AddTask(trdHitProducerQa);
     run->AddTask(new CbmTrdCalibTracker());
-    run->AddTask(new CbmTrackerInputQaTrd());  // Tracker requirements to TRD
     auto* pCaInputQaTrd = new CbmCaInputQaTrd(verbose, bUseMC);
     pCaInputQaTrd->SetEfficiencyThrsh(0.5, 0, 100);
     run->AddTask(pCaInputQaTrd);
@@ -226,7 +225,6 @@ void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "d
 
   // ----- TOF QA  ---------------------------------
   if (CbmSetup::Instance()->IsActive(ECbmModuleId::kTof)) {
-    //run->AddTask(new CbmTrackerInputQaTof());  // Tracker requirements to TOF
     auto* pCaInputQaTof = new CbmCaInputQaTof(verbose, bUseMC);
     pCaInputQaTof->SetEfficiencyThrsh(0.5, 0, 100);
     run->AddTask(pCaInputQaTof);
diff --git a/reco/L1/qa/CbmCaInputQaMuch.cxx b/reco/L1/qa/CbmCaInputQaMuch.cxx
index 272a448535..7e2a1d651e 100644
--- a/reco/L1/qa/CbmCaInputQaMuch.cxx
+++ b/reco/L1/qa/CbmCaInputQaMuch.cxx
@@ -144,19 +144,18 @@ bool CbmCaInputQaMuch::Check()
 
     // Function to fit a residual distribution, returns a structure
     auto FitResidualsAndGetMean = [&](TH1* pHist) {
-      auto* pFit = (TF1*) gROOT->FindObject("gausn");
-      LOG_IF(fatal, !pFit) << fName << ": residuals fit function is not found";
+      auto fit        = TF1("fitres", "gausn");
       double statMean = pHist->GetMean();
       double statSigm = pHist->GetStdDev();
-      pFit->SetParameters(100., statMean, statSigm);
-      pHist->Fit(pFit, "MQ");
-      pHist->Fit(pFit, "MQ");
-      pHist->Fit(pFit, "MQ");
+      fit.SetParameters(100., statMean, statSigm);
+      pHist->Fit("fitres", "MQ");
+      pHist->Fit("fitres", "MQ");
+      pHist->Fit("fitres", "MQ");
       // NOTE: Several fit procedures are made to avoid empty fit results
       std::array<double, 3> result;
-      result[0] = pFit->GetParameter(1);
-      result[1] = -pFit->GetParameter(2) * fResMeanThrsh;
-      result[2] = +pFit->GetParameter(2) * fResMeanThrsh;
+      result[0] = fit.GetParameter(1);
+      result[1] = -fit.GetParameter(2) * fResMeanThrsh;
+      result[2] = +fit.GetParameter(2) * fResMeanThrsh;
       return result;
     };
 
@@ -197,11 +196,10 @@ bool CbmCaInputQaMuch::Check()
     //TF1* pPullFit = new TF1("pullFitGausn", "[0] * Expk(-[2] * (x - [1]) * (x - [1]), [3])", -10., 10.);
 
     auto FitPullsAndGetSigma = [&](TH1* pHist) {
-      auto* pFit = (TF1*) gROOT->FindObject("gausn");
-      LOG_IF(fatal, !pFit) << fName << ": pulls fit function is not found";
-      pFit->SetParameters(100, 0.001, 1.000);
-      pHist->Fit(pFit, "Q");
-      return pFit->GetParameter(2);
+      auto fit = TF1("fitpull", "gausn(0)");
+      fit.SetParameters(100, 0.001, 1.000);
+      pHist->Fit("fitpull", "Q");
+      return fit.GetParameter(2);
     };
 
     auto* pPullsTable = MakeTable("pulls_rms", "Pulls std. dev. values in different stations", nSt, 4);
diff --git a/reco/L1/qa/CbmCaInputQaSts.cxx b/reco/L1/qa/CbmCaInputQaSts.cxx
index 26c4b9dd42..181950dd31 100644
--- a/reco/L1/qa/CbmCaInputQaSts.cxx
+++ b/reco/L1/qa/CbmCaInputQaSts.cxx
@@ -145,19 +145,18 @@ bool CbmCaInputQaSts::Check()
 
     // Function to fit a residual distribution, returns a structure
     auto FitResidualsAndGetMean = [&](TH1* pHist) {
-      auto* pFit = (TF1*) gROOT->FindObject("gausn");
-      LOG_IF(fatal, !pFit) << fName << ": residuals fit function is not found";
+      auto fit        = TF1("fitres", "gausn");
       double statMean = pHist->GetMean();
       double statSigm = pHist->GetStdDev();
-      pFit->SetParameters(100., statMean, statSigm);
-      pHist->Fit(pFit, "MQ");
-      pHist->Fit(pFit, "MQ");
-      pHist->Fit(pFit, "MQ");
+      fit.SetParameters(100., statMean, statSigm);
+      pHist->Fit("fitres", "MQ");
+      pHist->Fit("fitres", "MQ");
+      pHist->Fit("fitres", "MQ");
       // NOTE: Several fit procedures are made to avoid empty fit results
       std::array<double, 3> result;
-      result[0] = pFit->GetParameter(1);
-      result[1] = -pFit->GetParameter(2) * fResMeanThrsh;
-      result[2] = +pFit->GetParameter(2) * fResMeanThrsh;
+      result[0] = fit.GetParameter(1);
+      result[1] = -fit.GetParameter(2) * fResMeanThrsh;
+      result[2] = +fit.GetParameter(2) * fResMeanThrsh;
       return result;
     };
 
@@ -198,11 +197,10 @@ bool CbmCaInputQaSts::Check()
     //TF1* pPullFit = new TF1("pullFitGausn", "[0] * Expk(-[2] * (x - [1]) * (x - [1]), [3])", -10., 10.);
 
     auto FitPullsAndGetSigma = [&](TH1* pHist) {
-      auto* pFit = (TF1*) gROOT->FindObject("gausn");
-      LOG_IF(fatal, !pFit) << fName << ": pulls fit function is not found";
-      pFit->SetParameters(100, 0.001, 1.000);
-      pHist->Fit(pFit, "Q");
-      return pFit->GetParameter(2);
+      auto fit = TF1("fitpull", "gausn(0)");
+      fit.SetParameters(100, 0.001, 1.000);
+      pHist->Fit("fitpull", "Q");
+      return fit.GetParameter(2);
     };
 
     auto* pPullsTable = MakeTable("pulls_rms", "Pulls std. dev. values in different stations", nSt, 4);
diff --git a/reco/L1/qa/CbmCaInputQaTof.cxx b/reco/L1/qa/CbmCaInputQaTof.cxx
index db984e9c1d..e8903f8e8b 100644
--- a/reco/L1/qa/CbmCaInputQaTof.cxx
+++ b/reco/L1/qa/CbmCaInputQaTof.cxx
@@ -147,19 +147,18 @@ bool CbmCaInputQaTof::Check()
 
     // Function to fit a residual distribution, returns a structure
     auto FitResidualsAndGetMean = [&](TH1* pHist) {
-      auto* pFit = (TF1*) gROOT->FindObject("gausn");
-      LOG_IF(fatal, !pFit) << fName << ": residuals fit function is not found";
+      auto fit        = TF1("fitres", "gausn");
       double statMean = pHist->GetMean();
       double statSigm = pHist->GetStdDev();
-      pFit->SetParameters(100., statMean, statSigm);
-      pHist->Fit(pFit, "MQ");
-      pHist->Fit(pFit, "MQ");
-      pHist->Fit(pFit, "MQ");
+      fit.SetParameters(100., statMean, statSigm);
+      pHist->Fit("fitres", "MQ");
+      pHist->Fit("fitres", "MQ");
+      pHist->Fit("fitres", "MQ");
       // NOTE: Several fit procedures are made to avoid empty fit results
       std::array<double, 3> result;
-      result[0] = pFit->GetParameter(1);
-      result[1] = -pFit->GetParameter(2) * fResMeanThrsh;
-      result[2] = +pFit->GetParameter(2) * fResMeanThrsh;
+      result[0] = fit.GetParameter(1);
+      result[1] = -fit.GetParameter(2) * fResMeanThrsh;
+      result[2] = +fit.GetParameter(2) * fResMeanThrsh;
       return result;
     };
 
@@ -200,11 +199,10 @@ bool CbmCaInputQaTof::Check()
     //TF1* pPullFit = new TF1("pullFitGausn", "[0] * Expk(-[2] * (x - [1]) * (x - [1]), [3])", -10., 10.);
 
     auto FitPullsAndGetSigma = [&](TH1* pHist) {
-      auto* pFit = (TF1*) gROOT->FindObject("gausn");
-      LOG_IF(fatal, !pFit) << fName << ": pulls fit function is not found";
-      pFit->SetParameters(100, 0.001, 1.000);
-      pHist->Fit(pFit, "Q");
-      return pFit->GetParameter(2);
+      auto fit = TF1("fitpull", "gausn(0)");
+      fit.SetParameters(100, 0.001, 1.000);
+      pHist->Fit("fitpull", "Q");
+      return fit.GetParameter(2);
     };
 
     auto* pPullsTable = MakeTable("pulls_rms", "Pulls std. dev. values in different stations", nSt, 4);
@@ -305,6 +303,8 @@ void CbmCaInputQaTof::FillHistograms()
     int32_t address = pHit->GetAddress();
     if (0x00202806 == address || 0x00002806 == address) { continue; }
 
+    if (5 == CbmTofAddress::GetSmType(pHit->GetAddress())) { continue; }
+
     // *************************
     // ** Reconstructed hit QA
 
diff --git a/reco/L1/qa/CbmCaInputQaTrd.cxx b/reco/L1/qa/CbmCaInputQaTrd.cxx
index 827f066b38..d9a85050ed 100644
--- a/reco/L1/qa/CbmCaInputQaTrd.cxx
+++ b/reco/L1/qa/CbmCaInputQaTrd.cxx
@@ -143,19 +143,18 @@ bool CbmCaInputQaTrd::Check()
 
     // Function to fit a residual distribution, returns a structure
     auto FitResidualsAndGetMean = [&](TH1* pHist) {
-      auto* pFit = (TF1*) gROOT->FindObject("gausn");
-      LOG_IF(fatal, !pFit) << fName << ": residuals fit function is not found";
+      auto fit        = TF1("fitres", "gausn");
       double statMean = pHist->GetMean();
       double statSigm = pHist->GetStdDev();
-      pFit->SetParameters(100., statMean, statSigm);
-      pHist->Fit(pFit, "MQ");
-      pHist->Fit(pFit, "MQ");
-      pHist->Fit(pFit, "MQ");
+      fit.SetParameters(100., statMean, statSigm);
+      pHist->Fit("fitres", "MQ");
+      pHist->Fit("fitres", "MQ");
+      pHist->Fit("fitres", "MQ");
       // NOTE: Several fit procedures are made to avoid empty fit results
       std::array<double, 3> result;
-      result[0] = pFit->GetParameter(1);
-      result[1] = -pFit->GetParameter(2) * fResMeanThrsh;
-      result[2] = +pFit->GetParameter(2) * fResMeanThrsh;
+      result[0] = fit.GetParameter(1);
+      result[1] = -fit.GetParameter(2) * fResMeanThrsh;
+      result[2] = +fit.GetParameter(2) * fResMeanThrsh;
       return result;
     };
 
@@ -196,11 +195,10 @@ bool CbmCaInputQaTrd::Check()
     //TF1* pPullFit = new TF1("pullFitGausn", "[0] * Expk(-[2] * (x - [1]) * (x - [1]), [3])", -10., 10.);
 
     auto FitPullsAndGetSigma = [&](TH1* pHist) {
-      auto* pFit = (TF1*) gROOT->FindObject("gausn");
-      LOG_IF(fatal, !pFit) << fName << ": pulls fit function is not found";
-      pFit->SetParameters(100, 0.001, 1.000);
-      pHist->Fit(pFit, "Q");
-      return pFit->GetParameter(2);
+      auto fit = TF1("fitpull", "gausn(0)");
+      fit.SetParameters(100, 0.001, 1.000);
+      pHist->Fit("fitpull", "Q");
+      return fit.GetParameter(2);
     };
 
     auto* pPullsTable = MakeTable("pulls_rms", "Pulls std. dev. values in different stations", nSt, 4);
@@ -393,7 +391,6 @@ void CbmCaInputQaTrd::FillHistograms()
       double t0MC = fpMCEventList->GetEventTime(bestPointLink);
       LOG_IF(fatal, t0MC < 0) << fName << ": MC time zero is lower then 0 ns: " << t0MC;
 
-
       // ----- MC point properties
       //
       double mass = pMCTrack->GetMass();
@@ -434,7 +431,6 @@ void CbmCaInputQaTrd::FillHistograms()
       if (std::fabs(pMCPoint->GetPzOut()) < fMinMomentum) { continue; }  // CUT ON MINIMUM MOMENTUM
       //if (pMCo < cuts::kMinP) { continue; }  // Momentum threshold
 
-
       fvph_point_ypos_vs_xpos[iSt]->Fill(xMC, yMC);
       fvph_point_xpos_vs_zpos[iSt]->Fill(zMC, xMC);
       fvph_point_ypos_vs_zpos[iSt]->Fill(zMC, yMC);
@@ -932,7 +928,7 @@ InitStatus CbmCaInputQaTrd::InitHistograms()
     fvpe_reco_eff_vs_xy.resize(nSt + 2, nullptr);
     fvpe_reco_eff_vs_r.resize(nSt + 2, nullptr);
 
-    for (int iSt = 0; iSt <= nSt; ++iSt) {
+    for (int iSt = 0; iSt < nSt + 2; ++iSt) {
       TString nsuff = "";  // Histogram name suffix
       TString tsuff = "";  // Histogram title suffix
 
-- 
GitLab