tracker input QA for TRD

core/data/CbmLink.h

@@ -61,6 +61,8 @@ public:
     return (lhs.GetFile() == rhs.GetFile() && lhs.GetEntry() == rhs.GetEntry() && lhs.GetIndex() == rhs.GetIndex());
   }
 
+  friend Bool_t operator!=(const CbmLink& lhs, const CbmLink& rhs) { return !(lhs == rhs); }
+
   /** Comparison operators by //Dr.Sys **/
   friend Bool_t operator<(const CbmLink& l, const CbmLink& r)
   {

core/data/CbmMCEventList.cxx

@@ -97,6 +97,17 @@ Double_t CbmMCEventList::GetEventTime(UInt_t eventId, UInt_t fileId)
 // ----------------------------------------------------------------------------
 
 
+// -----   Get event index in the list   -------------------------------------
+Int_t CbmMCEventList::GetEventIndex(UInt_t eventId, UInt_t fileId)
+{
+  if (!fIsSorted) Sort();
+  auto it = Find(fileId, eventId);
+  if (it == fEvents.end()) return -1.;
+  return (it - fEvents.begin());
+}
+// ----------------------------------------------------------------------------
+
+
 // -----   Get event time for event at index in list   ------------------------
 Double_t CbmMCEventList::GetEventTimeByIndex(UInt_t index)
 {

core/data/CbmMCEventList.h

@@ -10,6 +10,7 @@
 #ifndef CBMMCEVENTLIST_H
 #define CBMMCEVENTLIST_H 1
 
+#include "CbmLink.h"
 #include "CbmMCEventInfo.h"  // for CbmMCEventInfo
 
 #include <Rtypes.h>      // for THashConsistencyHolder, ClassDef
@@ -72,6 +73,10 @@ public:
      **/
   Double_t GetEventTime(UInt_t event, UInt_t file);
 
+  /** @brief Event start time
+   * * @param link link to the MC event
+   **/
+  Double_t GetEventTime(const CbmLink& link) { return GetEventTime(link.GetEntry(), link.GetFile()); }
 
   /** @brief Event time by index
      ** @value Event time for event at given index in list
@@ -88,6 +93,19 @@ public:
      **/
   Int_t GetFileIdByIndex(UInt_t index);
 
+  /** @brief Event index
+     ** @param event  MC event number
+     ** @param file   MC input file number
+     ** @value eventindex in list
+     **
+     ** Returns -1. if the event is not present in the list.
+     **/
+  Int_t GetEventIndex(UInt_t event, UInt_t file);
+
+  /** @brief Event index
+     **/
+  Int_t GetEventIndex(const CbmLink& link) { return GetEventIndex(link.GetEntry(), link.GetFile()); }
+
 
   /** @brief Number of events in the list
      ** @value Number of events

core/qa/CMakeLists.txt

@@ -24,6 +24,12 @@ CbmQaPie.cxx
 CbmQaHist.cxx
 )
 
+set(HEADERS
+CbmQaCanvas.h
+CbmQaPie.h
+CbmQaHist.h
+)
+
 set(LINKDEF CbmQaBaseLinkDef.h)
 Set(LIBRARY_NAME CbmQaBase)
 Set(DEPENDENCIES Hist Gpad)

core/qa/CbmQaBaseLinkDef.h

@@ -11,10 +11,16 @@
 #pragma link C++ nestedclasses;
 #pragma link C++ nestedtypedef;
 
+// custom streamers are already implemented in the class body
 #pragma link C++ class CbmQaCanvas - ;
-#pragma link C++ class CbmQaPieSlice + ;
 #pragma link C++ class CbmQaPie - ;
-#pragma link C++ class CbmQaH1F + ;
-#pragma link C++ class CbmQaH1D + ;
+
+// create streamers automatically
+#pragma link C++ class CbmQaPieSlice + ;
+#pragma link C++ class CbmQaHist < TH1F> + ;
+#pragma link C++ class CbmQaHist < TH1D> + ;
+#pragma link C++ class CbmQaHist < TH1I> + ;
+#pragma link C++ class CbmQaHist < TProfile> + ;
+#pragma link C++ class CbmQaHist < TProfile2D> + ;
 
 #endif

core/qa/CbmQaHist.cxx

@@ -13,6 +13,9 @@ templateClassImp(CbmQaHist);
 
 // create definitions of specific CbmQaHist classes,
 // otherwise they will not be linked by the ROOT linker
-//
+
 template class CbmQaHist<TH1F>;
 template class CbmQaHist<TH1D>;
+template class CbmQaHist<TH1I>;
+template class CbmQaHist<TProfile>;
+template class CbmQaHist<TProfile2D>;

core/qa/CbmQaHist.h

@@ -12,13 +12,24 @@
 
 #include "CbmQaCanvas.h"
 
-#include "TFitResultPtr.h"
-#include "TH1D.h"
-#include "TH1F.h"
-#include "TStyle.h"
-#include "TVirtualPad.h"
-/// A modification of TH* classes that keeps statistics & fit drawing options
-/// and resizes the stat window accordingly without actual drawing.
+#include <FairLogger.h>
+
+#include <TFitResultPtr.h>
+#include <TPaveStats.h>
+#include <TStyle.h>
+#include <TVirtualPad.h>
+
+// The TH* headers are needed here for the ROOT linker
+
+#include <TF1.h>
+#include <TH1D.h>
+#include <TH1F.h>
+#include <TH1I.h>
+#include <TProfile.h>
+#include <TProfile2D.h>
+
+/// A modification of TH* and TProfile* classes that keeps statistics & fit drawing options
+/// and resizes the stat window accordingly without the actual drawing.
 /// In the original classes, hist->Draw() & canv->Update() must be called
 /// for resetting Stat/Fit window.
 ///
@@ -34,6 +45,7 @@ public:
       fOptStat = gStyle->GetOptStat();
       fOptFit  = gStyle->GetOptFit();
     }
+    this->SetLineWidth(2);
   }
 
   /// Copy constructor
@@ -45,6 +57,7 @@ public:
   CbmQaHist(Types... args) : HistTypeT(args...)
   {
     if (gStyle) { SetOptStatFit(gStyle->GetOptStat(), gStyle->GetOptFit()); }
+    this->SetLineWidth(2);
   }
 
   /// Destructor
@@ -57,54 +70,63 @@ public:
   TFitResultPtr Fit(Types... args)
   {
     TVirtualPad* padsav = gPad;
-    GetCanvas().cd();
+    GetDummyCanvas().cd();
+    this->Sumw2();
     auto ret = HistTypeT::Fit(args...);
-    GetCanvas().Clear();
+    GetDummyCanvas().Clear();
+
+    // make the output look nice
+
     if (padsav) padsav->cd();
     auto* f = this->GetFunction("gaus");
     if (f) {
-      f->SetParName(0, "Peak");
-      f->SetParName(1, "#mu");
-      f->SetParName(2, "#sigma");
+      f->SetParNames("Peak", "#mu", "#sigma");
       f->SetLineColor(kRed);
       f->SetLineWidth(3);
+      TPaveStats* st = (TPaveStats*) this->FindObject("stats");
+      if (!st) { LOG(fatal) << "CbmQaHist: can not access histogram statistics"; }
+      else {
+        st->SetX1NDC(0.6);
+        st->SetX2NDC(0.940);
+        st->SetY1NDC(0.5);
+        st->SetY2NDC(0.930);
+        st->SetOptStat(111110);
+        st->SetOptFit(10001);
+      }
     }
     return ret;
   }
 
-  /// Set stat drawing options and resize the stat window
+  /// Set stat drawing options and autoresize the stat window
   void SetOptStat(Int_t stat = 1) { SetOptStatFit(stat, fOptFit); }
 
-  /// Set fit drawing options and resize the stat window
+  /// Set fit drawing options and autoresize the stat window
   void SetOptFit(Int_t fit = 1) { SetOptStatFit(fOptStat, fit); }
 
-  /// Set stat & fit drawing options and resize the stat window
+  /// Set stat & fit drawing options and autoresize the stat window
   void SetOptStatFit(int stat, int fit)
   {
+    // the only way to create and auto-size the stat window is to draw the histogram
+
     fOptStat = stat;
     fOptFit  = fit;
-    if (!gStyle) { return; }
+    if (!gStyle) { return; }  // should not happen
+
     TVirtualPad* savePad = gPad;
     int saveStat         = gStyle->GetOptStat();
     int saveFit          = gStyle->GetOptFit();
 
-    this->SetStats(0);  // remove the old stat window
-    this->SetStats(1);  // set the flag to create thes stat window during Draw()
-    CbmQaHist* tmp = (CbmQaHist*) this->Clone("myClone");
-    GetCanvas().cd();
+    GetDummyCanvas().cd();
     gStyle->SetOptStat(fOptStat);
     gStyle->SetOptFit(fOptFit);
-    // the only way to create and auto-size the stat window
-    tmp->Draw();
-    GetCanvas().Update();
-    // move the stat window to *this
-    TObject* obj = tmp->FindObject("stats");
-    // the list of functions seems always to exist, but let's check it
-    if (obj && this->GetListOfFunctions()) {
-      tmp->GetListOfFunctions()->Remove(obj);
-      this->GetListOfFunctions()->Add(obj);
-    }
-    GetCanvas().Clear();
+
+    this->SetStats(0);  // remove the old stat window
+    this->SetStats(1);  // set the flag to create the stat window during Draw()
+    this->Draw();
+    GetDummyCanvas().Update();
+    GetDummyCanvas().Clear();
+
+    // restore the environment
     gStyle->SetOptStat(saveStat);
     gStyle->SetOptFit(saveFit);
     if (savePad) savePad->cd();
@@ -112,12 +134,12 @@ public:
 
 private:
   /// a static canvas for temporary drawing
-  static CbmQaCanvas& GetCanvas()
+  static CbmQaCanvas& GetDummyCanvas()
   {
     /// the static variable will be initialised at the first call;
     /// deleted at the application end (c++11)
-    static CbmQaCanvas tmp("CbmQaTempCanvas", "CbmQaTempCanvas", 1, 1);
-    return tmp;
+    static CbmQaCanvas dummy("CbmQaTempCanvas", "CbmQaTempCanvas", 1, 1);
+    return dummy;
   }
 
   int fOptStat = 1;
@@ -126,8 +148,4 @@ private:
   ClassDefNV(CbmQaHist, 1);
 };
 
-// shortcuts
-typedef CbmQaHist<TH1F> CbmQaH1F;
-typedef CbmQaHist<TH1D> CbmQaH1D;
-
 #endif

macro/run/run_qa.C

@@ -48,7 +48,9 @@ void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "d
 
   // -----   Logger settings   ----------------------------------------------
   FairLogger::GetLogger()->SetLogScreenLevel("INFO");
-  FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
+  fair::Logger::DefineVerbosity(
+    "user1", fair::VerbositySpec::Make(fair::VerbositySpec::Info::severity, fair::VerbositySpec::Info::file_line));
+  FairLogger::GetLogger()->SetLogVerbosityLevel("user1");
   // ------------------------------------------------------------------------
 
   // -----   Environment   --------------------------------------------------
@@ -173,6 +175,7 @@ void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "d
     //run->AddTask(new CbmTrdDigitizerPRFQa()); //works put currently doesn't do anything
     //run->AddTask(new CbmTrdHitRateFastQa());  //opens lots of windows
     run->AddTask(new CbmTrdHitProducerQa());
+    run->AddTask(new CbmTrackerInputQaTrd());  // Tracker requirements to TRD
   }
   // ------------------------------------------------------------------------
 

reco/L1/CMakeLists.txt

@@ -41,6 +41,7 @@ ${CBMROOT_SOURCE_DIR}/reco/L1
 ${CBMROOT_SOURCE_DIR}/reco/L1/L1Algo
 ${CBMROOT_SOURCE_DIR}/reco/L1/OffLineInterface
 ${CBMROOT_SOURCE_DIR}/reco/L1/ParticleFinder
+${CBMROOT_SOURCE_DIR}/reco/L1/qa
 
 ${CBMROOT_SOURCE_DIR}/reco/KF
 ${CBMROOT_SOURCE_DIR}/reco/KF/KFQA
@@ -52,7 +53,7 @@ ${CBMROOT_SOURCE_DIR}/sim/transport/geosetup
 ${CBMDATA_DIR}
 ${CBMDATA_DIR}/base
 
-  ${CBMDETECTORBASE_DIR}/sts
+${CBMDETECTORBASE_DIR}/sts
 
 ${CBMROOT_SOURCE_DIR}/mvd
   ${CBMDETECTORBASE_DIR}/trd
@@ -144,6 +145,8 @@ L1Algo/utils/L1AlgoPulls.cxx
 
 ParticleFinder/CbmL1PFFitter.cxx
 ParticleFinder/CbmL1PFMCParticle.cxx
+
+qa/CbmTrackerInputQaTrd.cxx
 )
 
 set(HEADERS
@@ -175,6 +178,8 @@ OffLineInterface/CbmL1GlobalFindTracksEvents.h
 #OffLineInterface / CbmL1SttTrack.h
 L1Algo/L1Def.h
 L1Algo/L1Vector.h
+
+qa/CbmTrackerInputQaTrd.h
 )
 
 
@@ -231,6 +236,7 @@ Set(DEPENDENCIES
     CbmStsBase
     CbmRecoBase
     CbmRecoSts
+    CbmQaBase
     boost_regex
    )
 

reco/L1/L1LinkDef.h

@@ -29,5 +29,7 @@
 //#pragma link C++ class  CbmL1SttHit+;
 //#pragma link C++ class  CbmL1SttTrackFinder+;
 //#pragma link C++ class  CbmL1SttTrack+;
+#pragma link C++ class CbmTrackerInputQaTrd + ;
+
 
 #endif

reco/L1/qa/CbmTrackerInputQaTrd.h

+/* Copyright (C) 2021 Facility for Antiproton and Ion Research in Europe, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergey Gorbunov[committer]*/
+
+/// @file CbmTrackerInputQaTrd.h
+/// @author Sergey Gorbunov
+/// @date 02.11.2021
+
+#ifndef CbmTrackerInputQaTrd_H
+#define CbmTrackerInputQaTrd_H
+
+#include "CbmQaCanvas.h"
+#include "CbmQaHist.h"
+
+#include <FairTask.h>
+
+#include <TFolder.h>
+#include <TH1D.h>
+#include <TH1F.h>
+#include <TH1I.h>
+#include <TH2F.h>
+#include <TParameter.h>
+#include <TProfile.h>
+#include <TProfile2D.h>
+
+class CbmDigiManager;
+class CbmMCDataManager;
+class CbmMCEventList;
+class CbmMCDataArray;
+class CbmTimeSlice;
+class CbmTrdParSetGeo;
+class CbmTrdParSetDigi;
+
+class TClonesArray;
+
+///
+/// CbmTrackerInputQaTrd class represents all the CA tracker requirements for the TRD detector.
+/// When this QA test is passed, the tracker must work (at least that is the idea).
+///
+/// The class ensures that the tracker has the correct understanding of the TRD geometry and interfaces
+/// and checks the quality of all tracking-related aspects of the TRD data.
+///
+class CbmTrackerInputQaTrd : public FairTask {
+
+public:
+  /// Constructor
+  CbmTrackerInputQaTrd(const char* name = "TrackerInputQaTrd", Int_t verbose = 1);
+
+  /// Destructor
+  ~CbmTrackerInputQaTrd();
+
+  /// FairTask: Intialisation at begin of run.
+  InitStatus Init() { return ReInit(); }
+
+  /// FairTask: Reinitialisation.
+  InitStatus ReInit();
+
+  /// FairTask:  Intialise parameter containers.
+  void SetParContainers();
+
+  /// FairTask: Action at end of run. For this task and all of the subtasks.
+  void Finish();
+
+  /// TTask: Clear all data structures created by a previous execution of a task.
+  void Clear(Option_t* /*option*/ = "") {}
+
+  /// TTask: Process a timeslice
+  void Exec(Option_t*);
+
+  /// Prepare the Qa output and return it as a reference to an internal TFolder.
+  /// The reference is non-const, because FairSink can not write const objects
+  TFolder& GetQa();
+
+private:
+  /// Check the geometry
+  InitStatus GeometryQa();
+
+  /// Analysis of hit uncertainty (pull) distributions
+  void ResolutionQa();
+
+  /// Free the memory etc.
+  void DeInit();
+
+
+  // Setup
+
+  Bool_t fIsTrdInSetup {false};
+  Bool_t fIsMcPresent {false};
+
+  Int_t fNtrackingStations {0};
+
+  CbmTimeSlice* fTimeSlice {nullptr};
+  CbmTrdParSetGeo* fTrdGeoPar {nullptr};
+  CbmTrdParSetDigi* fTrdDigiPar {nullptr};
+  CbmDigiManager* fDigiManager {nullptr};
+
+  /// MC data
+  CbmMCEventList* fMcEventList {nullptr};  // list of MC events connected to the current time slice
+  CbmMCDataManager* fMcManager {nullptr};
+  CbmMCDataArray* fMcTracks {nullptr};
+  CbmMCDataArray* fMcPoints {nullptr};
+
+  /// Data
+  TClonesArray* fClusters {nullptr};
+  TClonesArray* fHits {nullptr};
+  TClonesArray* fHitMatches {nullptr};
+
+  /// Output
+
+
+  TFolder fOutFolder {"TrackerInputQaTrd", "TrackerInputQaTrd"};  /// output folder with histos and canvases
+  TFolder* fHistFolder {nullptr};                                 /// subfolder for histograms
+
+  TParameter<int> fNevents {"nEvents", 0};  /// number of processed events
+
+  /// Histogram for Residual Distribution
+  CbmQaHist<TH1D> fh1DresidualU {"h1DresidualU", "Trd1D: Residual U;(U_{reco} - U_{MC})(cm)", 100, -.5, .5};
+  CbmQaHist<TH1D> fh1DresidualV {"h1DresidualV", "Trd1D: Residual V;(V_{reco} - V_{MC})(cm)", 100, -10, 10};
+  CbmQaHist<TH1D> fh1DresidualT {"h1DresidualT", "Trd1D: Residual T;(T_{reco} - T_{MC})(ns)", 100, -50, 50};
+
+  CbmQaHist<TH1D> fh2DresidualX {"h2DresidualX", "Trd2D: Residual X;(X_{reco} - X_{MC})(cm)", 100, -5, 5};
+  CbmQaHist<TH1D> fh2DresidualY {"h2DresidualY", "Trd2D: Residual Y;(Y_{reco} - Y_{MC})(cm)", 100, -5, 5};
+  CbmQaHist<TH1D> fh2DresidualT {"h2DresidualT", "Trd2D: Residual T;(T_{reco} - T_{MC})(ns)", 100, -1000, 1000};
+
+  /// Histogram for PULL Distribution
+  CbmQaHist<TH1D> fh1DpullU {"h1DpullU", "Trd1D: Pull U;(U_{reco} - U_{MC}) / #sigmaU_{reco}", 100, -5, 5};
+  CbmQaHist<TH1D> fh1DpullV {"h1DpullV", "Trd1D: Pull V;(V_{reco} - V_{MC}) / #sigmaV_{reco}", 100, -5, 5};
+  CbmQaHist<TH1D> fh1DpullT {"h1DpullT", "Trd1D: Pull T;(T_{reco} - T_{MC}) / #sigmaT_{reco}", 100, -5, 5};
+
+  CbmQaHist<TH1D> fh2DpullX {"h2DpullX", "Trd2D: Pull X;(X_{reco} - X_{MC}) / #sigmaX_{reco}", 100, -5, 5};
+  CbmQaHist<TH1D> fh2DpullY {"h2DpullY", "Trd2D: Pull Y;(Y_{reco} - Y_{MC}) / #sigmaY_{reco}", 100, -5, 5};
+  CbmQaHist<TH1D> fh2DpullT {"h2DpullT", "Trd2D: Pull T;(T_{reco} - T_{MC}) / #sigmaT_{reco}", 100, -5, 5};
+
+  /// List of the above histogramms
+  std::vector<CbmQaHist<TH1D>*> fHistList;
+
+  /// hits purity
+  std::vector<CbmQaHist<TH1I>> fhPointsPerHit;
+
+  /// hits efficiency
+  std::vector<CbmQaHist<TH1I>> fhHitsPerPoint;
+
+  /// hits efficiency
+  std::vector<CbmQaHist<TProfile2D>> fhEfficiencyXY;
+  std::vector<CbmQaHist<TProfile>> fhEfficiencyR;
+
+  /// Canvaces: collection of histogramms
+
+  CbmQaCanvas fCanvResidual {"cResidual", "Residual Distribution", 3 * 500, 2 * 500};
+  CbmQaCanvas fCanvPull {"cPull", "Pull Distribution", 3 * 500, 2 * 500};
+  CbmQaCanvas fCanvEfficiencyXY {"cEfficiencyXY", "Efficiency XY: % reconstructed McPoint", 2 * 500, 2 * 500};
+  CbmQaCanvas fCanvEfficiencyR {"cEfficiencyR", "Efficiency R: % reconstructed McPoint", 2 * 500, 2 * 500};
+  CbmQaCanvas fCanvHitsPerPoint {"cHitsPerMcPoint", "Efficiency: Hits Per McPoint", 2 * 500, 2 * 500};
+  CbmQaCanvas fCanvPointsPerHit {"cMcPointsPerHit", "Purity: McPoints per Hit", 2 * 500, 2 * 500};
+
+private:
+  /// Suppressed copy constructor
+  CbmTrackerInputQaTrd(const CbmTrackerInputQaTrd&) = delete;
+
+  /// Suppressed assignment operator
+  CbmTrackerInputQaTrd& operator=(const CbmTrackerInputQaTrd&) = delete;
+
+  ClassDef(CbmTrackerInputQaTrd, 0);
+};
+
+#endif