From bbe28efdcf731dc6c9387b25f203af145e1f329f Mon Sep 17 00:00:00 2001
From: sgorbuno <se.gorbunov@gsi.de>
Date: Tue, 29 Sep 2020 11:51:54 +0000
Subject: [PATCH] qa task for MUCH transport level
---
macro/run/CMakeLists.txt | 12 +-
macro/run/run_qa.C | 245 +++++++---
macro/run/run_qa_old.C | 121 +++++
sim/detectors/much/CMakeLists.txt | 6 +-
sim/detectors/much/CbmMuchSimLinkDef.h | 1 +
sim/detectors/much/qa/CbmMuchTransportQa.cxx | 457 +++++++++++++++++++
sim/detectors/much/qa/CbmMuchTransportQa.h | 106 +++++
7 files changed, 877 insertions(+), 71 deletions(-)
create mode 100644 macro/run/run_qa_old.C
create mode 100644 sim/detectors/much/qa/CbmMuchTransportQa.cxx
create mode 100644 sim/detectors/much/qa/CbmMuchTransportQa.h
diff --git a/macro/run/CMakeLists.txt b/macro/run/CMakeLists.txt
index f917221899..eedd71c102 100644
--- a/macro/run/CMakeLists.txt
+++ b/macro/run/CMakeLists.txt
@@ -4,6 +4,7 @@ GENERATE_ROOT_TEST_SCRIPT(${CBMROOT_SOURCE_DIR}/macro/run/run_digi.C)
GENERATE_ROOT_TEST_SCRIPT(${CBMROOT_SOURCE_DIR}/macro/run/run_reco_event.C)
GENERATE_ROOT_TEST_SCRIPT(${CBMROOT_SOURCE_DIR}/macro/run/check_overlaps.C)
GENERATE_ROOT_TEST_SCRIPT(${CBMROOT_SOURCE_DIR}/macro/run/run_transport_qa.C)
+GENERATE_ROOT_TEST_SCRIPT(${CBMROOT_SOURCE_DIR}/macro/run/run_qa.C)
GENERATE_ROOT_TEST_SCRIPT(${CBMROOT_SOURCE_DIR}/macro/run/run_reco_tb_track.C)
GENERATE_ROOT_TEST_SCRIPT(${CBMROOT_SOURCE_DIR}/macro/run/run_recoqa_tb_track.C)
GENERATE_ROOT_TEST_SCRIPT(${CBMROOT_SOURCE_DIR}/macro/run/run_reco_analysis.C)
@@ -76,7 +77,16 @@ ForEach(setup IN LISTS cbm_setup)
set_tests_properties(${testname} PROPERTIES FIXTURES_REQUIRED ${fixture_reco_event_mc})
Set(fixture_digi_tb fixture_digi_tb_${testname})
set_tests_properties(${testname} PROPERTIES FIXTURES_SETUP ${fixture_digi_tb})
-
+
+ # --- Quality Assurance
+ Set(testname run_qa_event_mc_${setup})
+ Add_Test(${testname} ${CBMROOT_BINARY_DIR}/macro/run/run_qa.sh ${NumEvents} \"data/${setup}_test\" \"${setup}\")
+ Set_Tests_Properties(${testname} PROPERTIES TIMEOUT ${timeOutTime})
+ Set_Tests_Properties(${testname} PROPERTIES PASS_REGULAR_EXPRESSION "Test Passed;All ok")
+ Set_Tests_Properties(${testname} PROPERTIES FAIL_REGULAR_EXPRESSION "ERROR")
+ Set_Tests_Properties(${testname} PROPERTIES FAIL_REGULAR_EXPRESSION "FATAL")
+ Set_Tests_properties(${testname} PROPERTIES FIXTURES_REQUIRED ${fixture_reco_event_mc})
+
# # --- Reconstruction Quality Assurance
# Set(testname run_reco_event_qa_${setup})
# Add_Test(${testname} ${CBMROOT_BINARY_DIR}/macro/run/run_reco_analysis.sh \"data/${setup}_test\" {{\"MvdCluster\",\"CbmMvdCluster\",\"MVD_Cluster_Finder\"},{\"MvdHit\",\"CbmMvdHit\",\"MVD_Hit_Finder\"},{\"StsCluster\",\"CbmStsCluster\",\"StsFindCluster\"},{\"StsHit\",\"CbmStsHit\",\"StsFindHit\"},{\"StsTrack\",\"CbmStsTrack\",\"StsFindTracksEvents\"}} {{\"sts\",{${pullDevAllowed},1,1,1}},{\"mvd\",{${pullDevAllowed},1,1,1}}} 2500 ${uploadHistJPG})
diff --git a/macro/run/run_qa.C b/macro/run/run_qa.C
index bafa1a89f0..c77ad24d82 100644
--- a/macro/run/run_qa.C
+++ b/macro/run/run_qa.C
@@ -1,47 +1,120 @@
// --------------------------------------------------------------------------
//
-// Macro for STS QA
+// Macro for simulation & reconstruction QA
//
-// V. Friese 13/01/2006
+// The following naming conventions are assumed:
+// Raw data file: [dataset].event.raw.root
+// Transport file: [dataset].tra.root
+// Parameter file: [dataset].par.root
+// Reconstruction file: [dataset].rec.root
+//
+// S. Gorbunov 28/09/2020
//
// --------------------------------------------------------------------------
-void run_qa(Int_t nEvents = 1, const char* setupName = "sis100_electron") {
+
+// Includes needed for IDE
+#if !defined(__CLING__)
+#include <FairFileSource.h>
+#include <FairMonitor.h>
+#include <FairParAsciiFileIo.h>
+#include <FairParRootFileIo.h>
+#include <FairRootFileSink.h>
+#include <FairRunAna.h>
+#include <FairRuntimeDb.h>
+#include <FairSystemInfo.h>
+
+#include "CbmDefs.h"
+#include "CbmMCDataManager.h"
+#include "CbmMuchTransportQa.h"
+#include "CbmSetup.h"
+
+#include <TStopwatch.h>
+#endif
+
+void run_qa(Int_t nEvents = 0,
+ TString dataset = "data/sis100_muon_jpsi_test",
+ TString setup = "sis100_muon_jpsi") {
// ========================================================================
// Adjust this part according to your requirements
- // Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug)
- Int_t iVerbose = 1;
-
- TString outDir = "data/";
-
- // MC file
- TString simFile =
- outDir + setupName + "_test.mc.root"; // Input file (MC events)
+ // ----- Logger settings ----------------------------------------------
- // Reco file
- TString recFile = outDir + setupName + "_test.eds.root"; // Output file
-
- // Parameter file
- TString parFile = outDir + setupName + "_params.root"; // Parameter file
+ FairLogger::GetLogger()->SetLogScreenLevel("INFO");
+ FairLogger::GetLogger()->SetLogVerbosityLevel("LOW");
+ // ------------------------------------------------------------------------
- // Output file
- TString outFile = outDir + setupName + "_test.qa.root"; // Output file
+ // ----- Environment --------------------------------------------------
+ TString myName = "run_qa"; // this macro's name for screen output
+ TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ // ------------------------------------------------------------------------
- TString srcDir = gSystem->Getenv("VMCWORKDIR");
- TString paramDir = srcDir + "/parameters/";
+ // ----- In- and output file names ------------------------------------
+ TString rawFile = dataset + ".event.raw.root";
+ TString traFile = dataset + ".tra.root";
+ TString parFile = dataset + ".par.root";
+ TString recFile = dataset + ".rec.root";
+ TString sinkFile = dataset + ".qa.root";
+ // ------------------------------------------------------------------------
// ----- Load the geometry setup -------------------------------------
std::cout << std::endl;
- TString setupFile = srcDir + "/geometry/setup/setup_" + setupName + ".C";
- TString setupFunct = "setup_";
- setupFunct = setupFunct + setupName + "()";
- std::cout << "-I- "
- << "run_qa"
- << ": Loading macro " << setupFile << std::endl;
- gROOT->LoadMacro(setupFile);
- gROOT->ProcessLine(setupFunct);
+ std::cout << "-I- " << myName << ": Loading setup " << setup << std::endl;
+ CbmSetup::Instance()->LoadSetup(setup);
+ // You can modify the pre-defined setup by using
+ // CbmSetup::Instance()->RemoveModule(ESystemId) or
+ // CbmSetup::Instance()->SetModule(ESystemId, const char*, Bool_t) or
+ // CbmSetup::Instance()->SetActive(ESystemId, Bool_t)
+ // See the class documentation of CbmSetup.
+ // ------------------------------------------------------------------------
+
+ // ----- Parameter files as input to the runtime database -------------
+ std::cout << std::endl;
+ std::cout << "-I- " << myName << ": Defining parameter files " << std::endl;
+ TList* parFileList = new TList();
+ TString geoTag;
+
+ // - MUCH digitisation parameters
+ if (CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kMuch, geoTag)) {
+ bool mcbmFlag = geoTag.Contains("mcbm", TString::kIgnoreCase);
+ TString parFile = srcDir + "/parameters/much/much_";
+ parFile += (mcbmFlag) ? geoTag : geoTag(0, 4);
+ parFile += "_digi_sector.root";
+ { // init geometry from the file
+ TFile* f = new TFile(parFile, "R");
+ TObjArray* stations = (TObjArray*) f->Get("stations");
+ CbmMuchGeoScheme::Instance()->Init(stations, mcbmFlag);
+ }
+ }
+
+ // - TRD digitisation parameters
+ if (CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kTrd, geoTag)) {
+ const Char_t* npar[4] = {"asic", "digi", "gas", "gain"};
+ TObjString* trdParFile(NULL);
+ for (Int_t i(0); i < 4; i++) {
+ trdParFile = new TObjString(srcDir + "/parameters/trd/trd_" + geoTag + "."
+ + npar[i] + ".par");
+ parFileList->Add(trdParFile);
+ std::cout << "-I- " << myName << ": Using parameter file "
+ << trdParFile->GetString() << std::endl;
+ }
+ }
+
+ // - TOF digitisation parameters
+ if (CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kTof, geoTag)) {
+ TObjString* tofFile =
+ new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digi.par");
+ parFileList->Add(tofFile);
+ std::cout << "-I- " << myName << ": Using parameter file "
+ << tofFile->GetString() << std::endl;
+ TObjString* tofBdfFile =
+ new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digibdf.par");
+ parFileList->Add(tofBdfFile);
+ std::cout << "-I- " << myName << ": Using parameter file "
+ << tofBdfFile->GetString() << std::endl;
+ }
+ // ------------------------------------------------------------------------
// In general, the following parts need not be touched
// ========================================================================
@@ -51,71 +124,105 @@ void run_qa(Int_t nEvents = 1, const char* setupName = "sis100_electron") {
timer.Start();
// ------------------------------------------------------------------------
- // ----- Analysis run --------------------------------------------------
- FairRunAna* fRun = new FairRunAna();
- FairFileSource* inputSource = new FairFileSource(recFile);
- fRun->SetSource(inputSource);
- // fRun->SetInputFile(simFile);
- // fRun->AddFriend(recFile);
- // fRun->SetInputFile(recFile);
- // fRun->AddFriend(simFile);
- fRun->SetOutputFile(outFile);
- // Define output file for FairMonitor histograms
- TString monitorFile {outFile};
+ // ---- Debug option -------------------------------------------------
+ gDebug = 0;
+ // ------------------------------------------------------------------------
+
+ // ----- FairRunAna ---------------------------------------------------
+ FairRunAna* run = new FairRunAna();
+
+ FairFileSource* inputSource = new FairFileSource(rawFile);
+ inputSource->AddFriend(traFile);
+ inputSource->AddFriend(recFile);
+ run->SetSource(inputSource);
+ run->SetGenerateRunInfo(kFALSE);
+
+ FairRootFileSink* sink = new FairRootFileSink(sinkFile);
+ run->SetSink(sink);
+
+ TString monitorFile {sinkFile};
monitorFile.ReplaceAll("qa", "qa.monitor");
FairMonitor::GetMonitor()->EnableMonitor(kTRUE, monitorFile);
// ------------------------------------------------------------------------
-
- // ----- Parameter database --------------------------------------------
- FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
- FairParRootFileIo* parIo1 = new FairParRootFileIo();
- parIo1->open(parFile.Data(), "UPDATE");
- rtdb->setFirstInput(parIo1);
+ // ----- MCDataManager (legacy mode) -----------------------------------
+ CbmMCDataManager* mcManager = new CbmMCDataManager("MCDataManager", 1);
+ mcManager->AddFile(traFile);
+ run->AddTask(mcManager);
// ------------------------------------------------------------------------
- // ----- Intialise and run --------------------------------------------
- fRun->Init();
+ // ----- MUCH hit finder QA ---------------------------------
- rtdb->setOutput(parIo1);
- rtdb->saveOutput();
- rtdb->print();
+ if (CbmSetup::Instance()->IsActive(ECbmModuleId::kMuch)) {
+ run->AddTask(new CbmMuchTransportQa());
+ }
- cout << "Starting run" << endl;
- fRun->Run(0, nEvents);
// ------------------------------------------------------------------------
+ // ----- Parameter database --------------------------------------------
+ {
+ std::cout << std::endl << std::endl;
+ std::cout << "-I- " << myName << ": Set runtime DB" << std::endl;
+ FairRuntimeDb* rtdb = run->GetRuntimeDb();
+ FairParRootFileIo* parIo1 = new FairParRootFileIo();
+ parIo1->open(parFile.Data(), "in");
+ rtdb->setFirstInput(parIo1);
+ if (!parFileList->IsEmpty()) {
+ FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
+ parIo2->open(parFileList, "in");
+ rtdb->setSecondInput(parIo2);
+ }
+ rtdb->print();
+ }
+
+ // ----- Run initialisation -------------------------------------------
+ std::cout << std::endl;
+ std::cout << "-I- " << myName << ": Initialise run" << std::endl;
+ run->Init();
+
+ // ----- Start run ----------------------------------------------------
+ std::cout << std::endl << std::endl;
+ std::cout << "-I- " << myName << ": Starting run" << std::endl;
+ run->Run(0, nEvents);
+ // ------------------------------------------------------------------------
+
+ //std::cout<< "Test failed" << std::endl;
// ----- Finish -------------------------------------------------------
timer.Stop();
+
+ FairMonitor::GetMonitor()->Print();
Double_t rtime = timer.RealTime();
Double_t ctime = timer.CpuTime();
- cout << endl << endl;
- cout << "Macro finished succesfully." << endl;
- cout << "Output file is " << outFile << endl;
- cout << "Parameter file is " << parFile << endl;
- cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
- cout << endl;
+ std::cout << std::endl << std::endl;
+ std::cout << "Macro finished successfully." << std::endl;
+ std::cout << "Output file is " << sinkFile << std::endl;
+ std::cout << "Parameter file is " << parFile << std::endl;
+ std::cout << "Real time " << rtime << " s, CPU time " << ctime << " s"
+ << std::endl;
+ std::cout << std::endl;
+ std::cout << " Test passed" << std::endl;
+ std::cout << " All ok " << std::endl;
// ------------------------------------------------------------------------
+
+ // ----- Resource monitoring ------------------------------------------
// Extract the maximal used memory an add is as Dart measurement
// This line is filtered by CTest and the value send to CDash
FairSystemInfo sysInfo;
Float_t maxMemory = sysInfo.GetMaxMemory();
- cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
- cout << maxMemory;
- cout << "</DartMeasurement>" << endl;
+ std::cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
+ std::cout << maxMemory;
+ std::cout << "</DartMeasurement>" << std::endl;
Float_t cpuUsage = ctime / rtime;
- cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
- cout << cpuUsage;
- cout << "</DartMeasurement>" << endl;
-
- FairMonitor* tempMon = FairMonitor::GetMonitor();
- tempMon->Print();
+ std::cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
+ std::cout << cpuUsage;
+ std::cout << "</DartMeasurement>" << std::endl;
+ // ------------------------------------------------------------------------
+ // ----- Function needed for CTest runtime dependency -----------------
RemoveGeoManager();
- cout << " Test passed" << endl;
- cout << " All ok " << endl;
+ // ------------------------------------------------------------------------
}
diff --git a/macro/run/run_qa_old.C b/macro/run/run_qa_old.C
new file mode 100644
index 0000000000..51202c3919
--- /dev/null
+++ b/macro/run/run_qa_old.C
@@ -0,0 +1,121 @@
+// --------------------------------------------------------------------------
+//
+// Macro for STS QA
+//
+// V. Friese 13/01/2006
+//
+// --------------------------------------------------------------------------
+void run_qa_old(Int_t nEvents = 1, const char* setupName = "sis100_electron") {
+
+ // ========================================================================
+ // Adjust this part according to your requirements
+
+ // Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug)
+ Int_t iVerbose = 1;
+
+ TString outDir = "data/";
+
+ // MC file
+ TString simFile =
+ outDir + setupName + "_test.mc.root"; // Input file (MC events)
+
+ // Reco file
+ TString recFile = outDir + setupName + "_test.eds.root"; // Output file
+
+ // Parameter file
+ TString parFile = outDir + setupName + "_params.root"; // Parameter file
+
+ // Output file
+ TString outFile = outDir + setupName + "_test.qa.root"; // Output file
+
+ TString srcDir = gSystem->Getenv("VMCWORKDIR");
+ TString paramDir = srcDir + "/parameters/";
+
+
+ // ----- Load the geometry setup -------------------------------------
+ std::cout << std::endl;
+ TString setupFile = srcDir + "/geometry/setup/setup_" + setupName + ".C";
+ TString setupFunct = "setup_";
+ setupFunct = setupFunct + setupName + "()";
+ std::cout << "-I- "
+ << "run_qa"
+ << ": Loading macro " << setupFile << std::endl;
+ gROOT->LoadMacro(setupFile);
+ gROOT->ProcessLine(setupFunct);
+
+ // In general, the following parts need not be touched
+ // ========================================================================
+
+ // ----- Timer --------------------------------------------------------
+ TStopwatch timer;
+ timer.Start();
+ // ------------------------------------------------------------------------
+
+ // ----- Analysis run --------------------------------------------------
+ FairRunAna* fRun = new FairRunAna();
+ FairFileSource* inputSource = new FairFileSource(recFile);
+ fRun->SetSource(inputSource);
+ // fRun->SetInputFile(simFile);
+ // fRun->AddFriend(recFile);
+ // fRun->SetInputFile(recFile);
+ // fRun->AddFriend(simFile);
+ fRun->SetOutputFile(outFile);
+ // Define output file for FairMonitor histograms
+ TString monitorFile {outFile};
+ monitorFile.ReplaceAll("qa", "qa.monitor");
+ FairMonitor::GetMonitor()->EnableMonitor(kTRUE, monitorFile);
+ // ------------------------------------------------------------------------
+
+
+ // ----- Parameter database --------------------------------------------
+ FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
+ FairParRootFileIo* parIo1 = new FairParRootFileIo();
+ parIo1->open(parFile.Data(), "UPDATE");
+ rtdb->setFirstInput(parIo1);
+ // ------------------------------------------------------------------------
+
+
+ // ----- Intialise and run --------------------------------------------
+ fRun->Init();
+
+ rtdb->setOutput(parIo1);
+ rtdb->saveOutput();
+ rtdb->print();
+
+ cout << "Starting run" << endl;
+ fRun->Run(0, nEvents);
+ // ------------------------------------------------------------------------
+
+
+ // ----- Finish -------------------------------------------------------
+ timer.Stop();
+ Double_t rtime = timer.RealTime();
+ Double_t ctime = timer.CpuTime();
+ cout << endl << endl;
+ cout << "Macro finished succesfully." << endl;
+ cout << "Output file is " << outFile << endl;
+ cout << "Parameter file is " << parFile << endl;
+ cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
+ cout << endl;
+ // ------------------------------------------------------------------------
+
+ // Extract the maximal used memory an add is as Dart measurement
+ // This line is filtered by CTest and the value send to CDash
+ FairSystemInfo sysInfo;
+ Float_t maxMemory = sysInfo.GetMaxMemory();
+ cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
+ cout << maxMemory;
+ cout << "</DartMeasurement>" << endl;
+
+ Float_t cpuUsage = ctime / rtime;
+ cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
+ cout << cpuUsage;
+ cout << "</DartMeasurement>" << endl;
+
+ FairMonitor* tempMon = FairMonitor::GetMonitor();
+ tempMon->Print();
+
+ RemoveGeoManager();
+ cout << " Test passed" << endl;
+ cout << " All ok " << endl;
+}
diff --git a/sim/detectors/much/CMakeLists.txt b/sim/detectors/much/CMakeLists.txt
index ad082df4cd..3b2bb90cea 100644
--- a/sim/detectors/much/CMakeLists.txt
+++ b/sim/detectors/much/CMakeLists.txt
@@ -1,6 +1,7 @@
Set(INCLUDE_DIRECTORIES
${CMAKE_CURRENT_SOURCE_DIR}
+${CMAKE_CURRENT_SOURCE_DIR}/qa
${CBMDETECTORBASE_DIR}/much
@@ -8,6 +9,7 @@ ${CBMBASE_DIR}
${CBMBASE_DIR}/utils
${CBMROOT_SOURCE_DIR}/sim/transport/base
+${CBMROOT_SOURCE_DIR}/core/qa
${CBMDATA_DIR}
${CBMDATA_DIR}/much
@@ -36,12 +38,14 @@ CbmMuch.cxx
CbmMuchDigitizeGem.cxx
CbmMuchReadoutBuffer.cxx
CbmMuchSignal.cxx
+
+qa/CbmMuchTransportQa.cxx
)
set(LINKDEF CbmMuchSimLinkDef.h)
Set(LIBRARY_NAME CbmMuchSim)
Set(DEPENDENCIES
- CbmMuchBase CbmBase CbmData Base
+ CbmMuchBase CbmBase CbmData Base CbmQaBase
)
GENERATE_LIBRARY()
diff --git a/sim/detectors/much/CbmMuchSimLinkDef.h b/sim/detectors/much/CbmMuchSimLinkDef.h
index d46400411f..f2b5951cf5 100644
--- a/sim/detectors/much/CbmMuchSimLinkDef.h
+++ b/sim/detectors/much/CbmMuchSimLinkDef.h
@@ -12,5 +12,6 @@
#pragma link C++ class CbmMuchDigitizeGem + ;
#pragma link C++ class CbmMuchReadoutBuffer + ;
#pragma link C++ class CbmMuchSignal + ;
+#pragma link C++ class CbmMuchTransportQa + ;
#endif
diff --git a/sim/detectors/much/qa/CbmMuchTransportQa.cxx b/sim/detectors/much/qa/CbmMuchTransportQa.cxx
new file mode 100644
index 0000000000..0892cee664
--- /dev/null
+++ b/sim/detectors/much/qa/CbmMuchTransportQa.cxx
@@ -0,0 +1,457 @@
+/// \file CbmMuchTransportQa.cxx
+/// \brief Implementation of the CbmMuchTransportQa class
+/// \author Sergey Gorbunov <se.gorbunov@gsi.de>
+/// \author Eugeny Kryshen
+/// \author Vikas Singhal
+/// \author Ekata Nandy
+/// \date 25.09.2020
+
+#include "CbmMuchTransportQa.h"
+
+#include "FairLogger.h"
+#include "FairRootFileSink.h"
+#include "FairRootManager.h"
+#include "FairRun.h"
+#include "FairRuntimeDb.h"
+
+#include "CbmGeoMuchPar.h"
+#include "CbmMuchGeoScheme.h"
+#include "CbmMuchStation.h"
+
+#include "CbmMuchPoint.h"
+
+#include "CbmMCTrack.h"
+#include "TDatabasePDG.h"
+#include "TParticlePDG.h"
+
+
+#include "TString.h"
+
+#include "TClonesArray.h"
+
+#include "CbmQaCanvas.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TStyle.h"
+
+#include <cassert>
+#include <vector>
+
+ClassImp(CbmMuchTransportQa);
+
+// -------------------------------------------------------------------------
+CbmMuchTransportQa::CbmMuchTransportQa(const char* name, Int_t verbose)
+ : FairTask(name, verbose)
+ , fOutFolder("MuchTransportQA", "Much Transport QA")
+ , fhNevents("nEvents", 0)
+ , fvUsNtra()
+ , fvMcPointXY()
+ , fvMcPointPhiZ()
+ , fvMcPointRZ()
+ , fvFraction() {}
+// -------------------------------------------------------------------------
+
+
+// -------------------------------------------------------------------------
+CbmMuchTransportQa::~CbmMuchTransportQa() { DeInit(); }
+// -------------------------------------------------------------------------
+
+
+// -------------------------------------------------------------------------
+void CbmMuchTransportQa::DeInit() {
+
+ fPoints = nullptr;
+ fMcTracks = nullptr;
+
+ fOutFolder.Clear();
+ fhNevents.SetVal(0);
+
+ SafeDelete(fhUsNtraAll);
+ SafeDelete(fhUsNtraPrim);
+ SafeDelete(fhUsNtraSec);
+ SafeDelete(fhUsNtraPr);
+ SafeDelete(fhUsNtraPi);
+ SafeDelete(fhUsNtraEl);
+ SafeDelete(fhUsNtraMu);
+ SafeDelete(fhUsNtraKa);
+
+ for (uint i = 0; i < fvMcPointXY.size(); i++) {
+ SafeDelete(fvMcPointXY[i]);
+ }
+ fvMcPointXY.clear();
+ for (uint i = 0; i < fvMcPointPhiZ.size(); i++) {
+ SafeDelete(fvMcPointPhiZ[i]);
+ }
+ fvMcPointPhiZ.clear();
+ for (uint i = 0; i < fvMcPointRZ.size(); i++) {
+ SafeDelete(fvMcPointRZ[i]);
+ }
+ fvMcPointRZ.clear();
+
+ SafeDelete(fhNtracks);
+ SafeDelete(fhFractionPrim);
+ SafeDelete(fhFractionSec);
+ SafeDelete(fhFractionPr);
+ SafeDelete(fhFractionPi);
+ SafeDelete(fhFractionEl);
+ SafeDelete(fhFractionMu);
+ SafeDelete(fhFractionKa);
+
+ fvUsNtra.clear();
+ fvFraction.clear();
+
+ SafeDelete(fCanvStationXY);
+ SafeDelete(fCanvStationPhiZ);
+ SafeDelete(fCanvStationRZ);
+
+ fNstations = 0;
+ fOutFolder.Clear();
+}
+// -------------------------------------------------------------------------
+
+
+// -------------------------------------------------------------------------
+InitStatus CbmMuchTransportQa::Init() {
+
+ TDirectory* oldDirectory = gDirectory;
+
+ FairRootManager* manager = FairRootManager::Instance();
+ if (!manager) {
+ LOG(error) << "No FairRootManager found";
+ return kERROR;
+ }
+
+ fMcTracks = (TClonesArray*) manager->GetObject("MCTrack");
+ if (!fMcTracks) {
+ LOG(error) << "No MC tracks found";
+ return kERROR;
+ }
+
+ fPoints = (TClonesArray*) manager->GetObject("MuchPoint");
+ if (!fPoints) {
+ LOG(error) << "No MC points found";
+ return kERROR;
+ }
+
+ if (!CbmMuchGeoScheme::Instance()) {
+ LOG(fatal) << "No CbmMuchGeoScheme found";
+ return kFATAL;
+ }
+
+ fNstations = CbmMuchGeoScheme::Instance()->GetNStations();
+
+ if (fNstations == 0) {
+ LOG(error) << "CbmMuchGeoScheme is not initialized";
+ return kERROR;
+ }
+
+ TFolder* histFolder = fOutFolder.AddFolder("hist", "Histogramms");
+
+ fhNevents.SetVal(0);
+ histFolder->Add(&fhNevents);
+
+#define BINS_STA fNstations, 0, fNstations
+
+ {
+ fvUsNtra.clear();
+ std::vector<TH1F*>& v = fvUsNtra;
+ v.push_back(fhUsNtraAll = new TH1F("hUsNtraAll", "N tracks", BINS_STA));
+ v.push_back(fhUsNtraPrim =
+ new TH1F("hUsNtraPrim", "N primary tracks", BINS_STA));
+ v.push_back(fhUsNtraSec =
+ new TH1F("hUsNtraSec", "N secondary tracks", BINS_STA));
+ v.push_back(fhUsNtraPr = new TH1F("hUsNtraPr", "N protons", BINS_STA));
+ v.push_back(fhUsNtraPi = new TH1F("hUsNtraPi", "N pions", BINS_STA));
+ v.push_back(fhUsNtraEl = new TH1F("hUsNtraEl", "N electrons", BINS_STA));
+ v.push_back(fhUsNtraMu = new TH1F("hUsNtraMu", "N muons", BINS_STA));
+ v.push_back(fhUsNtraKa = new TH1F("hUsNtraKa", "N kaons", BINS_STA));
+ for (uint i = 0; i < fvUsNtra.size(); i++) {
+ TH1F* h = fvUsNtra[i];
+ h->SetStats(0);
+ h->GetXaxis()->SetTitle("Station");
+ histFolder->Add(h);
+ }
+ }
+
+ {
+ fvFraction.clear();
+ std::vector<TH1F*>& v = fvFraction;
+ v.push_back(fhNtracks =
+ new TH1F("hNtracks", "N tracks per event", BINS_STA));
+ v.push_back(fhFractionPrim = new TH1F(
+ "hFractionPrim", "Fraction of primary tracks", BINS_STA));
+ v.push_back(fhFractionSec = new TH1F(
+ "hFractionSec", "Fraction of secondary tracks", BINS_STA));
+ v.push_back(fhFractionPr =
+ new TH1F("hFractionPr", "Fraction of protons", BINS_STA));
+ v.push_back(fhFractionPi =
+ new TH1F("hFractionPi", "Fraction of pions", BINS_STA));
+ v.push_back(fhFractionEl =
+ new TH1F("hFractionEl", "Fraction of electrons", BINS_STA));
+ v.push_back(fhFractionMu =
+ new TH1F("hFractionMu", "Fraction of muons", BINS_STA));
+ v.push_back(fhFractionKa =
+ new TH1F("hFractionKa", "Fraction of kaons", BINS_STA));
+
+ for (uint i = 0; i < fvFraction.size(); i++) {
+ TH1F* h = fvFraction[i];
+ h->SetStats(0);
+ h->GetXaxis()->SetTitle("Station");
+ h->GetYaxis()->SetTitle("%");
+ histFolder->Add(h);
+ }
+ }
+
+ fvMcPointXY.resize(fNstations);
+ fvMcPointPhiZ.resize(fNstations);
+ fvMcPointRZ.resize(fNstations);
+
+ for (Int_t i = 0; i < fNstations; i++) {
+ CbmMuchStation* station = CbmMuchGeoScheme::Instance()->GetStation(i);
+ if (!station) {
+ LOG(fatal) << "Much station " << i << " doesn't exist";
+ return kFATAL;
+ }
+ Double_t rMax = station->GetRmax();
+ Double_t rMin = station->GetRmin();
+
+ fvMcPointXY[i] = new TH2F(Form("hMcPointXY%i", i + 1),
+ Form("MC point XY : Station %i; X; Y", i + 1),
+ 100,
+ -1.2 * rMax,
+ 1.2 * rMax,
+ 100,
+ -1.2 * rMax,
+ 1.2 * rMax);
+ histFolder->Add(fvMcPointXY[i]);
+
+ fvMcPointPhiZ[i] =
+ new TH2F(Form("hMcPointPhiZ%i", i + 1),
+ Form("MC point Phi vs Z : Station %i; Z; Phi", i + 1),
+ 100,
+ station->GetZ() - station->GetTubeDz() - 5.,
+ station->GetZ() + station->GetTubeDz() + 5.,
+ 100,
+ -200.,
+ 200.);
+ histFolder->Add(fvMcPointPhiZ[i]);
+
+ float dR = rMax - rMin;
+ fvMcPointRZ[i] = new TH2F(Form("hMcPointRZ%i", i + 1),
+ Form("MC point R vs Z : Station %i; Z; R", i + 1),
+ 100,
+ station->GetZ() - station->GetTubeDz() - 5.,
+ station->GetZ() + station->GetTubeDz() + 5.,
+ 100,
+ rMin - 0.1 * dR,
+ rMax + 0.1 * dR);
+ histFolder->Add(fvMcPointRZ[i]);
+ }
+
+ fCanvStationXY =
+ new CbmQaCanvas("cMcPointXY", "Much: MC point XY", 2 * 400, 2 * 400);
+ fCanvStationXY->Divide2D(fNstations);
+ fOutFolder.Add(fCanvStationXY);
+
+ fCanvStationPhiZ = new CbmQaCanvas(
+ "cMcPointPhiZ", "Much: MC point Phi vs Z", 2 * 800, 2 * 400);
+ fCanvStationPhiZ->Divide2D(fNstations);
+ fOutFolder.Add(fCanvStationPhiZ);
+
+ fCanvStationRZ =
+ new CbmQaCanvas("cMcPointRZ", "Much: MC point R vs Z", 2 * 800, 2 * 400);
+ fCanvStationRZ->Divide2D(fNstations);
+ fOutFolder.Add(fCanvStationRZ);
+
+ gDirectory = oldDirectory;
+
+ return kSUCCESS;
+}
+// -------------------------------------------------------------------------
+
+
+// -------------------------------------------------------------------------
+InitStatus CbmMuchTransportQa::ReInit() {
+ DeInit();
+ return Init();
+}
+// -------------------------------------------------------------------------
+
+
+// -------------------------------------------------------------------------
+void CbmMuchTransportQa::SetParContainers() {
+ // Get run and runtime database
+
+ // The code currently does not work,
+ // CbmMuchGeoScheme::Instance() must be initialised outside.
+ // - Sergey
+
+ // FairRuntimeDb* db = FairRuntimeDb::instance();
+ // if ( ! db ) Fatal("SetParContainers", "No runtime database");
+ // Get MUCH geometry parameter container
+ // CbmGeoMuchPar *fGeoPar = (CbmGeoMuchPar*)
+ // db->getContainer("CbmGeoMuchPar"); TObjArray *stations =
+ // fGeoPar->GetStations();
+ // TString geoTag;
+ // CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kMuch, geoTag);
+ // bool mcbmFlag = geoTag.Contains("mcbm", TString::kIgnoreCase);
+ // CbmMuchGeoScheme::Instance()->Init(stations, mcbmFlag);
+}
+// -------------------------------------------------------------------------
+
+// -------------------------------------------------------------------------
+void CbmMuchTransportQa::Exec(Option_t*) {
+
+ LOG(info) << "Event: " << fhNevents.GetVal();
+
+ fhNevents.SetVal(fhNevents.GetVal() + 1);
+
+ // bitmask tells which stations were crossed by mc track
+ std::vector<UInt_t> trackStaCross(fMcTracks->GetEntriesFast(), 0);
+
+ for (Int_t i = 0; i < fPoints->GetEntriesFast(); i++) {
+
+ CbmMuchPoint* point = (CbmMuchPoint*) fPoints->At(i);
+
+ if (!point) {
+ LOG(fatal) << "Much point " << i << " doesn't exist";
+ break;
+ }
+
+ Int_t stId = CbmMuchAddress::GetStationIndex(point->GetDetectorID());
+ UInt_t stMask = (1 << stId);
+
+ // Check if the point corresponds to a certain MC Track
+ Int_t trackId = point->GetTrackID();
+ if (trackId < 0 || trackId >= fMcTracks->GetEntriesFast()) {
+ LOG(fatal) << "Much point " << i << ": trackId " << trackId
+ << " doesn't belong to [0," << fMcTracks->GetEntriesFast() - 1
+ << "]";
+ break;
+ }
+
+ CbmMCTrack* mcTrack = (CbmMCTrack*) fMcTracks->At(trackId);
+ if (!mcTrack) {
+ LOG(fatal) << "MC track " << trackId << " doesn't exist";
+ break;
+ }
+
+ Int_t motherId = mcTrack->GetMotherId();
+
+ // Get mass
+ Int_t pdgCode = mcTrack->GetPdgCode();
+
+ //if (pdgCode == 0) continue;
+
+ TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pdgCode);
+ if (!particle) {
+ LOG(warning) << "Particle with pdg code " << pdgCode << " doesn't exist";
+ //continue;
+ }
+ if (pdgCode == 0 || pdgCode == 22 || // photons
+ pdgCode == 2112) // neutrons
+ {
+ LOG(warning) << "Particle with pdg code " << pdgCode
+ << " left an mc point";
+ //continue;
+ }
+
+ if (!(trackStaCross[trackId] & stMask)) {
+ fhUsNtraAll->Fill(stId);
+ if (motherId == -1) {
+ fhUsNtraPrim->Fill(stId);
+ } else {
+ fhUsNtraSec->Fill(stId);
+ }
+ switch (abs(pdgCode)) {
+ case 2212: // proton
+ fhUsNtraPr->Fill(stId);
+ break;
+ case 211: // pion
+ fhUsNtraPi->Fill(stId);
+ break;
+ case 11: // electron
+ fhUsNtraEl->Fill(stId);
+ break;
+ case 13: // muon
+ fhUsNtraMu->Fill(stId);
+ break;
+ case 321: // kaon
+ fhUsNtraKa->Fill(stId);
+ break;
+ }
+ }
+
+ trackStaCross[trackId] |= stMask;
+
+ TVector3 v1; // in position of the track
+ TVector3 v2; // out position of the track
+
+ point->PositionIn(v1);
+ point->PositionOut(v2);
+
+ fvMcPointXY[stId]->Fill(v1.X(), v1.Y());
+ fvMcPointXY[stId]->Fill(v2.X(), v2.Y());
+
+ fvMcPointPhiZ[stId]->Fill(v1.Z(), v1.Phi() * TMath::RadToDeg());
+ fvMcPointPhiZ[stId]->Fill(v2.Z(), v2.Phi() * TMath::RadToDeg());
+
+ fvMcPointRZ[stId]->Fill(v1.Z(), v1.Perp());
+ fvMcPointRZ[stId]->Fill(v2.Z(), v2.Perp());
+ }
+}
+// -------------------------------------------------------------------------
+
+
+// -------------------------------------------------------------------------
+TFolder& CbmMuchTransportQa::GetQa() {
+
+ TDirectory* oldDirectory = gDirectory;
+
+ fhNtracks->Reset();
+ fhNtracks->Add(fhUsNtraAll, 1. / fhNevents.GetVal());
+
+ {
+ std::vector<Double_t> errors(fNstations, 0.);
+ fhUsNtraAll->SetError(errors.data());
+ }
+
+ for (uint i = 1; i < fvFraction.size(); i++) {
+ fvFraction[i]->Divide(fvUsNtra[i], fhUsNtraAll);
+ fvFraction[i]->Scale(100.);
+ }
+
+ for (Int_t i = 0; i < fNstations; i++) {
+
+ fCanvStationXY->cd(i + 1);
+ gStyle->SetOptStat(0);
+ fvMcPointXY[i]->DrawCopy("colz", "");
+
+ fCanvStationPhiZ->cd(i + 1);
+ gStyle->SetOptStat(0);
+ fvMcPointPhiZ[i]->DrawCopy("colz", "");
+
+ fCanvStationRZ->cd(i + 1);
+ gStyle->SetOptStat(0);
+ fvMcPointRZ[i]->DrawCopy("colz", "");
+
+ gStyle->SetOptStat(1110);
+ }
+
+ gDirectory = oldDirectory;
+ return fOutFolder;
+}
+// -------------------------------------------------------------------------
+
+
+// -------------------------------------------------------------------------
+void CbmMuchTransportQa::Finish() {
+
+ if (!FairRootManager::Instance() || !FairRootManager::Instance()->GetSink()) {
+ LOG(error) << "No sink found";
+ return;
+ }
+ FairSink* sink = FairRootManager::Instance()->GetSink();
+ sink->WriteObject(&GetQa(), nullptr);
+}
+// -------------------------------------------------------------------------
diff --git a/sim/detectors/much/qa/CbmMuchTransportQa.h b/sim/detectors/much/qa/CbmMuchTransportQa.h
new file mode 100644
index 0000000000..b6c113ffe9
--- /dev/null
+++ b/sim/detectors/much/qa/CbmMuchTransportQa.h
@@ -0,0 +1,106 @@
+/// \file CbmMuchTransportQa.h
+/// \brief Definition of the CbmMuchTransportQa class
+/// \author Sergey Gorbunov <se.gorbunov@gsi.de>
+/// \author Eugeny Kryshen
+/// \author Vikas Singhal
+/// \author Ekata Nandy
+/// \date 25.09.2020
+
+#ifndef CbmMuchTransportQa_H
+#define CbmMuchTransportQa_H
+
+#include "FairTask.h"
+#include "TFolder.h"
+#include "TParameter.h"
+
+class CbmMuchGeoScheme;
+class TClonesArray;
+class TH1F;
+class TH2F;
+class CbmQaCanvas;
+
+/// QA for the MUCH detector after a "transport" step of the simulation.
+///
+/// The class reimplements corresponding QA checks from old CbmMuchHitFinderQa class
+/// made by E. Kryshen & V. Singhal & E. Nandy
+///
+class CbmMuchTransportQa : public FairTask {
+
+public:
+ /// Constructor
+ CbmMuchTransportQa(const char* name = "MuchHitFinderQa", Int_t verbose = 1);
+
+ /// Deactivated copy constructors
+ CbmMuchTransportQa(const CbmMuchTransportQa&) = delete;
+ CbmMuchTransportQa& operator=(const CbmMuchTransportQa&) = delete;
+
+ /// Destructor
+ virtual ~CbmMuchTransportQa();
+
+ /// FairTask methods
+ void SetParContainers();
+ InitStatus Init();
+ InitStatus ReInit();
+ void Exec(Option_t* option);
+ void Finish();
+
+ /// Prepare Qa output and return it as a reference to an internal folder.
+ /// The reference is non-const, because FairSink can not write const objects
+ TFolder& GetQa();
+
+private:
+ /// Reset varibles & deallocate memory.
+ /// When not called by destructor, need to be folloed by Init().
+ void DeInit();
+
+ /// geometry
+ Int_t fNstations = 0;
+
+ /// containers
+ TClonesArray* fPoints = nullptr;
+ TClonesArray* fMcTracks = nullptr;
+ ///
+
+ TFolder fOutFolder; /// output folder with histos and canvases
+ TParameter<int> fhNevents; /// number of processed events
+
+ /// internal unscaled histogramms
+
+ TH1F* fhUsNtraAll = nullptr; /// number of all tracks
+ TH1F* fhUsNtraPrim = nullptr; /// number of primary tracks
+ TH1F* fhUsNtraSec = nullptr; /// number of secondary tracks
+ TH1F* fhUsNtraPr = nullptr; /// number of protons
+ TH1F* fhUsNtraPi = nullptr; /// number of pions
+ TH1F* fhUsNtraEl = nullptr; /// number of electrons
+ TH1F* fhUsNtraMu = nullptr; /// number of muons
+ TH1F* fhUsNtraKa = nullptr; /// number of kaons
+
+ std::vector<TH1F*> fvUsNtra; /// pointers to the above fhUsNtra* histos
+
+ /// output histograms
+
+ std::vector<TH2F*> fvMcPointXY; /// MC point Y vs X [N stations]
+ std::vector<TH2F*> fvMcPointPhiZ; /// MC point Phi vs Z [N stations]
+ std::vector<TH2F*> fvMcPointRZ; /// MC point R vs Z [N stations]
+
+ TH1F* fhNtracks = nullptr; /// number of all tracks / event
+ TH1F* fhFractionPrim = nullptr; /// fraction of primary tracks
+ TH1F* fhFractionSec = nullptr; /// fraction of secondary tracks
+ TH1F* fhFractionPr = nullptr; /// fraction of protons
+ TH1F* fhFractionPi = nullptr; /// fraction of pions
+ TH1F* fhFractionEl = nullptr; /// fraction of electrons
+ TH1F* fhFractionMu = nullptr; /// fraction of muons
+ TH1F* fhFractionKa = nullptr; /// fraction of kaons
+
+ std::vector<TH1F*> fvFraction; /// pointers to the above histos
+
+ // output canvaces with histogramm collections
+
+ CbmQaCanvas* fCanvStationXY = nullptr;
+ CbmQaCanvas* fCanvStationPhiZ = nullptr;
+ CbmQaCanvas* fCanvStationRZ = nullptr;
+
+ ClassDef(CbmMuchTransportQa, 0)
+};
+
+#endif
--
GitLab