From c84cbb2e13071ea84c1f3b9812f5c5df7650dd0b Mon Sep 17 00:00:00 2001
From: Dominik Smith <d.smith@gsi.de>
Date: Wed, 21 Oct 2020 13:12:21 +0200
Subject: [PATCH] Cleaned up CbmMuchTransportQa. Added CbmMuchDigitizerQa.
---
macro/run/run_qa.C | 1 +
reco/detectors/much/CMakeLists.txt | 4 +-
reco/detectors/much/CbmMuchRecoLinkDef.h | 1 -
sim/detectors/much/CMakeLists.txt | 3 +
sim/detectors/much/CbmMuchSimLinkDef.h | 2 +
sim/detectors/much/qa/CbmMuchDigitizerQa.cxx | 1012 +++++++++++++++++
sim/detectors/much/qa/CbmMuchDigitizerQa.h | 199 ++++
.../detectors/much/qa}/CbmMuchPointInfo.cxx | 0
.../detectors/much/qa}/CbmMuchPointInfo.h | 0
sim/detectors/much/qa/CbmMuchTransportQa.cxx | 473 ++++----
sim/detectors/much/qa/CbmMuchTransportQa.h | 62 +-
11 files changed, 1507 insertions(+), 250 deletions(-)
create mode 100644 sim/detectors/much/qa/CbmMuchDigitizerQa.cxx
create mode 100644 sim/detectors/much/qa/CbmMuchDigitizerQa.h
rename {reco/detectors/much => sim/detectors/much/qa}/CbmMuchPointInfo.cxx (100%)
rename {reco/detectors/much => sim/detectors/much/qa}/CbmMuchPointInfo.h (100%)
diff --git a/macro/run/run_qa.C b/macro/run/run_qa.C
index dffdca3721..3aa12ef371 100644
--- a/macro/run/run_qa.C
+++ b/macro/run/run_qa.C
@@ -157,6 +157,7 @@ void run_qa(Int_t nEvents = 0,
if (CbmSetup::Instance()->IsActive(ECbmModuleId::kMuch)) {
run->AddTask(new CbmMuchTransportQa());
+ run->AddTask(new CbmMuchDigitizerQa());
}
// ------------------------------------------------------------------------
diff --git a/reco/detectors/much/CMakeLists.txt b/reco/detectors/much/CMakeLists.txt
index e9ef35f9ec..b900b0c7a0 100644
--- a/reco/detectors/much/CMakeLists.txt
+++ b/reco/detectors/much/CMakeLists.txt
@@ -4,6 +4,7 @@ ${CMAKE_CURRENT_SOURCE_DIR}
${CBMDETECTORBASE_DIR}/much
${CBMROOT_SOURCE_DIR}/reco/base
+${CBMROOT_SOURCE_DIR}/sim/detectors/much/qa
${CBMBASE_DIR}
@@ -33,7 +34,6 @@ set(SRCS
CbmMuchFindHitsGem.cxx
CbmMuchHitFinderQa.cxx
CbmMuchHitProducerIdeal.cxx
-CbmMuchPointInfo.cxx
CbmMuchFindTracks.cxx
CbmMuchMatchTracks.cxx
@@ -43,7 +43,7 @@ CbmMuchTrackFinderIdeal.cxx
set(LINKDEF CbmMuchRecoLinkDef.h)
Set(LIBRARY_NAME CbmMuchReco)
Set(DEPENDENCIES
- CbmMuchBase CbmRecoBase CbmBase CbmData Base
+ CbmMuchBase CbmMuchSim CbmRecoBase CbmBase CbmData Base
)
GENERATE_LIBRARY()
diff --git a/reco/detectors/much/CbmMuchRecoLinkDef.h b/reco/detectors/much/CbmMuchRecoLinkDef.h
index bb4cd67f13..e6c3fb6852 100644
--- a/reco/detectors/much/CbmMuchRecoLinkDef.h
+++ b/reco/detectors/much/CbmMuchRecoLinkDef.h
@@ -9,7 +9,6 @@
#pragma link C++ class CbmMuchHitFinderQa + ;
#pragma link C++ class CbmMuchHitProducerIdeal + ;
#pragma link C++ class CbmMuchMatchTracks + ;
-#pragma link C++ class CbmMuchPointInfo + ;
#pragma link C++ class CbmMuchTrackFinderIdeal + ;
#endif
diff --git a/sim/detectors/much/CMakeLists.txt b/sim/detectors/much/CMakeLists.txt
index 3b2bb90cea..6c83a88aa2 100644
--- a/sim/detectors/much/CMakeLists.txt
+++ b/sim/detectors/much/CMakeLists.txt
@@ -9,6 +9,7 @@ ${CBMBASE_DIR}
${CBMBASE_DIR}/utils
${CBMROOT_SOURCE_DIR}/sim/transport/base
+${CBMROOT_SOURCE_DIR}/core/data/base
${CBMROOT_SOURCE_DIR}/core/qa
${CBMDATA_DIR}
@@ -40,6 +41,8 @@ CbmMuchReadoutBuffer.cxx
CbmMuchSignal.cxx
qa/CbmMuchTransportQa.cxx
+qa/CbmMuchPointInfo.cxx
+qa/CbmMuchDigitizerQa.cxx
)
set(LINKDEF CbmMuchSimLinkDef.h)
diff --git a/sim/detectors/much/CbmMuchSimLinkDef.h b/sim/detectors/much/CbmMuchSimLinkDef.h
index f2b5951cf5..67b350f261 100644
--- a/sim/detectors/much/CbmMuchSimLinkDef.h
+++ b/sim/detectors/much/CbmMuchSimLinkDef.h
@@ -13,5 +13,7 @@
#pragma link C++ class CbmMuchReadoutBuffer + ;
#pragma link C++ class CbmMuchSignal + ;
#pragma link C++ class CbmMuchTransportQa + ;
+#pragma link C++ class CbmMuchPointInfo + ;
+#pragma link C++ class CbmMuchDigitizerQa + ;
#endif
diff --git a/sim/detectors/much/qa/CbmMuchDigitizerQa.cxx b/sim/detectors/much/qa/CbmMuchDigitizerQa.cxx
new file mode 100644
index 0000000000..2df89e438c
--- /dev/null
+++ b/sim/detectors/much/qa/CbmMuchDigitizerQa.cxx
@@ -0,0 +1,1012 @@
+/// \file CbmMuchDigitizerQa.cxx
+/// \brief Implementation of the CbmMuchDigitizerQa class
+/// \author Sergey Gorbunov <se.gorbunov@gsi.de>
+/// \author Eugeny Kryshen
+/// \author Vikas Singhal
+/// \author Ekata Nandy
+/// \author Dominik Smith
+/// \date 21.10.2020
+
+#include "CbmMuchDigitizerQa.h"
+using std::cout;
+using std::endl;
+using std::vector;
+
+#define BINNING_CHARGE 100, 0, 300.0
+#define BINNING_LENGTH 100, 0, 2.5
+#define BINNING_CHARGE_LOG 100, 3, 10
+#define BINNING_ENERGY_LOG 100, -0.5, 4.5
+#define BINNING_ENERGY_LOG_EL 100, -0.5, 4.5
+
+// -------------------------------------------------------------------------
+CbmMuchDigitizerQa::CbmMuchDigitizerQa(const char* name, Int_t verbose)
+ : FairTask(name, verbose)
+ , fGeoScheme(CbmMuchGeoScheme::Instance())
+ , fDigiManager(CbmDigiManager::Instance())
+ , fPointInfos(new TClonesArray("CbmMuchPointInfo", 10))
+ , fOutFolder("MuchDigiQA", "MuchDigitizerQA")
+ , fvUsPadsFiredR()
+ , fvUsPadsFiredXY()
+ , fvMcPointCharge()
+ , fvPadsTotalR()
+ , fvPadsFiredR()
+ , fvPadOccupancyR() {}
+// -------------------------------------------------------------------------
+
+// -------------------------------------------------------------------------
+CbmMuchDigitizerQa::~CbmMuchDigitizerQa() { DeInit(); }
+// -------------------------------------------------------------------------
+
+void CbmMuchDigitizerQa::DeInit() {
+ for (int i = 0; i < fNstations; i++) {
+ delete fvUsPadsFiredR[i];
+ delete fvUsPadsFiredXY[i];
+ delete fvMcPointCharge[i];
+ delete fvPadsTotalR[i];
+ delete fvPadsFiredR[i];
+ delete fvPadOccupancyR[i];
+ }
+ fvUsPadsFiredR.clear();
+ fvUsPadsFiredXY.clear();
+ fvMcPointCharge.clear();
+ fvPadsTotalR.clear();
+ fvPadsFiredR.clear();
+ fvPadOccupancyR.clear();
+
+ fNstations = 0;
+ fOutFolder.Clear();
+}
+
+// -------------------------------------------------------------------------
+InitStatus CbmMuchDigitizerQa::Init() {
+
+ TDirectory* oldDirectory = gDirectory;
+ FairRootManager* fManager = FairRootManager::Instance();
+ fMCTracks = (TClonesArray*) fManager->GetObject("MCTrack");
+ fPoints = (TClonesArray*) fManager->GetObject("MuchPoint");
+ // Reading Much Digis from CbmMuchDigiManager which are stored as vector
+ fDigiManager = CbmDigiManager::Instance();
+ fDigiManager->Init();
+
+ histFolder = fOutFolder.AddFolder("hist", "Histogramms");
+ fNstations = fGeoScheme->GetNStations();
+ printf("Init: fNstations = %i\n", fNstations);
+
+ //fVerbose = 3;
+ InitCounters();
+ InitCanvases();
+ InitChargeHistos();
+ InitLengthHistos();
+ InitPadHistos();
+ FillTotalPadsHistos();
+ InitChannelPadInfo();
+ InitFits();
+
+ gDirectory = oldDirectory;
+ return kSUCCESS;
+}
+
+void CbmMuchDigitizerQa::InitCounters() {
+
+ fNall = new Int_t[fNstations];
+ fNpr = new Int_t[fNstations];
+ fNpi = new Int_t[fNstations];
+ fNel = new Int_t[fNstations];
+ fNmu = new Int_t[fNstations];
+ fNka = new Int_t[fNstations];
+ fNprimary = new Int_t[fNstations];
+ fNsecondary = new Int_t[fNstations];
+
+ for (Int_t i = 0; i < fNstations; i++) {
+ fNall[i] = 0;
+ fNpr[i] = 0;
+ fNpi[i] = 0;
+ fNel[i] = 0;
+ fNmu[i] = 0;
+ fNka[i] = 0;
+ fNprimary[i] = 0;
+ fNsecondary[i] = 0;
+ }
+ fPointsTotal = 0;
+ fPointsUnderCounted = 0;
+ fPointsOverCounted = 0;
+}
+
+void CbmMuchDigitizerQa::InitChannelPadInfo() {
+
+ fPadMinLx = std::numeric_limits<Double_t>::max();
+ fPadMinLy = std::numeric_limits<Double_t>::max();
+ fPadMaxLx = std::numeric_limits<Double_t>::min();
+ fPadMaxLy = std::numeric_limits<Double_t>::min();
+
+ if (fVerbose > 2) {
+ printf("=========================================================\n");
+ printf(" Station Nr.\t| Sectors\t| Channels\t| Pad min size\t\t| Pad max"
+ "length\t \n");
+ printf("---------------------------------------------------------\n");
+ }
+
+ Int_t nTotSectors = 0;
+ Int_t nTotChannels = 0;
+ for (Int_t iStation = 0; iStation < fNstations; iStation++) {
+ Int_t nChannels = GetNChannels(iStation);
+ Int_t nSectors = GetNSectors(iStation);
+ Double_t padMinLx = GetMinPadSize(iStation).X();
+ Double_t padMinLy = GetMinPadSize(iStation).Y();
+ Double_t padMaxLx = GetMaxPadSize(iStation).X();
+ Double_t padMaxLy = GetMaxPadSize(iStation).Y();
+ if (fPadMinLx > padMinLx) fPadMinLx = padMinLx;
+ if (fPadMinLy > padMinLy) fPadMinLy = padMinLy;
+ if (fPadMaxLx < padMaxLx) fPadMaxLx = padMaxLx;
+ if (fPadMaxLy < padMaxLy) fPadMaxLy = padMaxLy;
+ nTotSectors += nSectors;
+ nTotChannels += nChannels;
+
+ if (fVerbose > 2) {
+ printf("%i\t\t| %i\t\t| %i\t| %5.4fx%5.4f\t\t| %5.4fx%5.4f\n",
+ iStation + 1,
+ nSectors,
+ nChannels,
+ padMinLx,
+ padMinLy,
+ padMaxLx,
+ padMaxLy);
+ printf("%i\t\t| %i\t\t\n", iStation + 1, nChannels);
+ }
+ }
+ printf("-----------------------------------------------------------\n");
+ printf(" Total:\t\t| %i\t\t\n", nTotChannels);
+ printf("===========================================================\n");
+}
+
+void CbmMuchDigitizerQa::InitCanvases() {
+
+ fCanvCharge =
+ new CbmQaCanvas("cMcPointCharge", "MC point charge", 2 * 800, 2 * 400);
+ fCanvCharge->Divide2D(3);
+
+ fCanvStationCharge = new CbmQaCanvas(
+ "cMcPointChargeVsStation", "MC point charge per station", 2 * 800, 2 * 400);
+ fCanvStationCharge->Divide2D(fNstations);
+
+ fCanvChargeVsEnergy = new CbmQaCanvas("cMcPointChargeVsEnergy",
+ "MC point charge vs particle Energy",
+ 2 * 800,
+ 2 * 400);
+ fCanvChargeVsEnergy->Divide2D(4);
+
+ fCanvChargeVsLength = new CbmQaCanvas("cMcPointChargeVsLength",
+ "MC point charge vs track length",
+ 2 * 800,
+ 2 * 400);
+ fCanvChargeVsLength->Divide2D(4);
+
+ fCanvTrackLength =
+ new CbmQaCanvas("cTrackLength", "track length", 2 * 800, 2 * 400);
+ fCanvTrackLength->Divide2D(4);
+
+ fCanvNpadsVsArea =
+ new CbmQaCanvas("cNpadsVsArea", "N pads Vs Area", 2 * 800, 2 * 400);
+
+ fCanvUsPadsFiredXY = new CbmQaCanvas(
+ "cPadsFiredXY", "Number of pads fired vs XY", 2 * 800, 2 * 400);
+ fCanvUsPadsFiredXY->Divide2D(fNstations);
+
+ fCanvPadOccupancyR = new CbmQaCanvas(
+ "cPadOccupancyR", "Pad occupancy [%] vs radius", 2 * 800, 2 * 400);
+ fCanvPadOccupancyR->Divide2D(fNstations);
+
+ fCanvPadsTotalR =
+ new CbmQaCanvas("cPadsTotalR", "Total pads vs radius", 2 * 800, 2 * 400);
+ fCanvPadsTotalR->Divide2D(fNstations);
+
+ fOutFolder.Add(fCanvCharge);
+ fOutFolder.Add(fCanvStationCharge);
+ fOutFolder.Add(fCanvChargeVsEnergy);
+ fOutFolder.Add(fCanvChargeVsLength);
+ fOutFolder.Add(fCanvTrackLength);
+ fOutFolder.Add(fCanvNpadsVsArea);
+ fOutFolder.Add(fCanvUsPadsFiredXY);
+ fOutFolder.Add(fCanvPadOccupancyR);
+ fOutFolder.Add(fCanvPadsTotalR);
+}
+
+void CbmMuchDigitizerQa::InitChargeHistos() {
+
+ fhMcPointCharge =
+ new TH1F("hCharge", "Charge distribution from tracks", BINNING_CHARGE);
+ fhMcPointCharge->GetXaxis()->SetTitle("Charge [10^{4} electrons]");
+
+ fhMcPointChargeLog = new TH1F(
+ "hChargeLog", "Charge (log.) distribution from tracks", BINNING_CHARGE_LOG);
+ fhMcPointChargeLog->GetXaxis()->SetTitle("Charge [Lg(Number of electrons)]");
+
+ fhMcPointChargePr_1GeV_3mm =
+ new TH1F("hChargePr_1GeV_3mm", "Charge for 1 GeV protons", BINNING_CHARGE);
+ fhMcPointChargePr_1GeV_3mm->GetXaxis()->SetTitle("Charge [10^{4} electrons]");
+
+ fhMcPointChargeVsTrackEnergyLog =
+ new TH2F("hChargeEnergyLog",
+ "Charge vs energy (log.) for all tracks",
+ BINNING_ENERGY_LOG,
+ BINNING_CHARGE);
+
+ fhMcPointChargeVsTrackEnergyLogPi =
+ new TH2F("hChargeEnergyLogPi",
+ "Charge vs energy (log.) for pions",
+ BINNING_ENERGY_LOG,
+ BINNING_CHARGE);
+
+ fhMcPointChargeVsTrackEnergyLogPr =
+ new TH2F("hChargeEnergyLogPr",
+ "Charge vs energy (log.) for protons",
+ BINNING_ENERGY_LOG,
+ BINNING_CHARGE);
+
+ fhMcPointChargeVsTrackEnergyLogEl =
+ new TH2F("hChargeEnergyLogEl",
+ "Charge vs energy (log.) for electrons",
+ BINNING_ENERGY_LOG_EL,
+ BINNING_CHARGE);
+
+ fhMcPointChargeVsTrackLength = new TH2F("hChargeTrackLength",
+ "Charge vs length for all tracks",
+ BINNING_LENGTH,
+ BINNING_CHARGE);
+
+ fhMcPointChargeVsTrackLengthPi = new TH2F("hChargeTrackLengthPi",
+ "Charge vs length for pions",
+ BINNING_LENGTH,
+ BINNING_CHARGE);
+
+ fhMcPointChargeVsTrackLengthPr = new TH2F("hChargeTrackLengthPr",
+ "Charge vs length for proton",
+ BINNING_LENGTH,
+ BINNING_CHARGE);
+
+ fhMcPointChargeVsTrackLengthEl = new TH2F("hChargeTrackLengthEl",
+ "Charge vs length for electrons",
+ BINNING_LENGTH,
+ BINNING_CHARGE);
+ std::vector<TH2F*> vChargeHistos;
+ vChargeHistos.push_back(fhMcPointChargeVsTrackEnergyLog);
+ vChargeHistos.push_back(fhMcPointChargeVsTrackEnergyLogPi);
+ vChargeHistos.push_back(fhMcPointChargeVsTrackEnergyLogPr);
+ vChargeHistos.push_back(fhMcPointChargeVsTrackEnergyLogEl);
+ vChargeHistos.push_back(fhMcPointChargeVsTrackLength);
+ vChargeHistos.push_back(fhMcPointChargeVsTrackLengthPi);
+ vChargeHistos.push_back(fhMcPointChargeVsTrackLengthPr);
+ vChargeHistos.push_back(fhMcPointChargeVsTrackLengthEl);
+
+ for (UInt_t i = 0; i < vChargeHistos.size(); i++) {
+ TH2F* histo = vChargeHistos[i];
+ histo->SetStats(0);
+ histo->GetYaxis()->SetDecimals(kTRUE);
+ histo->GetYaxis()->SetTitleOffset(1.4);
+ histo->GetYaxis()->CenterTitle();
+ histo->GetYaxis()->SetTitle("Charge [10^{4} electrons]");
+ if (i < 4) {
+ histo->GetXaxis()->SetTitle("Lg(E_{kin} [MeV])");
+ } else {
+ histo->GetXaxis()->SetTitle("Track length [cm]");
+ }
+ histFolder->Add(histo);
+ }
+}
+
+void CbmMuchDigitizerQa::InitLengthHistos() {
+
+ fhTrackLength = new TH1F("hTrackLength", "Track length", BINNING_LENGTH);
+
+ fhTrackLengthPi =
+ new TH1F("hTrackLengthPi", "Track length for pions", BINNING_LENGTH);
+
+ fhTrackLengthPr =
+ new TH1F("hTrackLengthPr", "Track length for protons", BINNING_LENGTH);
+
+ fhTrackLengthEl =
+ new TH1F("hTrackLengthEl", "Track length for electrons", BINNING_LENGTH);
+
+ std::vector<TH1F*> vLengthHistos;
+ vLengthHistos.push_back(fhTrackLength);
+ vLengthHistos.push_back(fhTrackLengthPi);
+ vLengthHistos.push_back(fhTrackLengthPr);
+ vLengthHistos.push_back(fhTrackLengthEl);
+
+ for (UInt_t i = 0; i < vLengthHistos.size(); i++) {
+ TH1F* histo = vLengthHistos[i];
+ histo->GetXaxis()->SetTitle("Track length [cm]");
+ histFolder->Add(histo);
+ }
+}
+
+void CbmMuchDigitizerQa::InitPadHistos() {
+
+ fvMcPointCharge.resize(fNstations);
+ fvPadsTotalR.resize(fNstations);
+ fvUsPadsFiredR.resize(fNstations);
+ fvUsPadsFiredXY.resize(fNstations);
+ fvPadsFiredR.resize(fNstations);
+ fvPadOccupancyR.resize(fNstations);
+
+ for (Int_t i = 0; i < fNstations; i++) {
+ CbmMuchStation* station = fGeoScheme->GetStation(i);
+ Double_t rMax = station->GetRmax();
+ Double_t rMin = station->GetRmin();
+
+ fvMcPointCharge[i] = new TH1F(
+ Form("hMcPointCharge%i", i + 1),
+ Form("MC point charge : Station %i; Charge [10^4 e]; Count", i + 1),
+ BINNING_CHARGE);
+
+ fvPadsTotalR[i] =
+ new TH1F(Form("hPadsTotal%i", i + 1),
+ Form("Number of pads vs radius: Station %i;Radius [cm]", i + 1),
+ 100,
+ 0.6 * rMin,
+ 1.2 * rMax);
+
+ fvUsPadsFiredR[i] = new TH1F(
+ Form("hUsPadsFired%i", i + 1),
+ Form("Number of fired pads vs radius: Station %i;Radius [cm]", i + 1),
+ 100,
+ 0.6 * rMin,
+ 1.2 * rMax);
+
+ fvUsPadsFiredXY[i] =
+ new TH2F(Form("hUsPadsFiredXY%i", i + 1),
+ Form("Pads fired XY : Station %i; X; Y", i + 1),
+ 100,
+ -1.2 * rMax,
+ 1.2 * rMax,
+ 100,
+ -1.2 * rMax,
+ 1.2 * rMax);
+
+ fvPadsFiredR[i] = new TH1F(
+ Form("hPadsFired%i", i + 1),
+ Form("Number of fired pads vs radius: Station %i;Radius [cm]", i + 1),
+ 100,
+ 0.6 * rMin,
+ 1.2 * rMax);
+
+ fvPadOccupancyR[i] = new TH1F(
+ Form("hOccupancy%i", i + 1),
+ Form("Pad occupancy vs radius: Station %i;Radius [cm];Occupancy, %%",
+ i + 1),
+ 100,
+ 0.6 * rMin,
+ 1.2 * rMax);
+
+ histFolder->Add(fvMcPointCharge[i]);
+ histFolder->Add(fvPadsTotalR[i]);
+ histFolder->Add(fvUsPadsFiredXY[i]);
+ histFolder->Add(fvPadsFiredR[i]);
+ histFolder->Add(fvPadOccupancyR[i]);
+ }
+ fhNpadsVsS = new TH2F("hNpadsVsS",
+ "Number of fired pads vs pad area:area:n pads",
+ 10,
+ -5,
+ 0,
+ 10,
+ 0.5,
+ 10.5);
+ histFolder->Add(fhNpadsVsS);
+}
+
+void CbmMuchDigitizerQa::FillTotalPadsHistos() {
+
+ vector<CbmMuchModule*> modules = fGeoScheme->GetModules();
+ for (vector<CbmMuchModule*>::iterator im = modules.begin();
+ im != modules.end();
+ im++) {
+ CbmMuchModule* mod = (*im);
+ if (mod->GetDetectorType() == 4
+ || mod->GetDetectorType() == 3) { // modified for rpc
+ CbmMuchModuleGem* module = (CbmMuchModuleGem*) mod;
+ vector<CbmMuchPad*> pads = module->GetPads();
+ for (UInt_t ip = 0; ip < pads.size(); ip++) {
+ CbmMuchPad* pad = pads[ip];
+ Int_t stationId = CbmMuchAddress::GetStationIndex(pad->GetAddress());
+ Double_t x0 = pad->GetX();
+ Double_t y0 = pad->GetY();
+ Double_t r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
+ fvPadsTotalR[stationId]->Fill(r0);
+ }
+ }
+ }
+
+ Double_t errors[100];
+ for (Int_t i = 0; i < 100; i++) {
+ errors[i] = 0;
+ }
+ for (Int_t i = 0; i < fNstations; i++) {
+ //fvPadsTotalR[i]->Sumw2();
+ fvPadsTotalR[i]->SetError(errors);
+ }
+}
+
+void CbmMuchDigitizerQa::InitFits() {
+
+ fFitEl = new TF1("fit_el", LandauMPV, -0.5, 4.5, 1);
+ fFitEl->SetParameter(0, 0.511);
+ fFitEl->SetLineWidth(2);
+ fFitEl->SetLineColor(kBlack);
+
+ fFitPi = new TF1("fit_pi", LandauMPV, -0.5, 4.5, 1);
+ fFitPi->SetParameter(0, 139.57);
+ fFitPi->SetLineWidth(2);
+ fFitPi->SetLineColor(kBlack);
+
+ fFitPr = new TF1("fit_pr", LandauMPV, -0.5, 4.5, 1);
+ fFitPr->SetParameter(0, 938.27);
+ fFitPr->SetLineWidth(2);
+ fFitPr->SetLineColor(kBlack);
+}
+
+// -------------------------------------------------------------------------
+
+// -------------------------------------------------------------------------
+void CbmMuchDigitizerQa::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(); cout<<"\n\n SG: stations:
+ // "<<stations->GetEntriesFast()<<endl;
+ // TString geoTag;
+ // CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kMuch, geoTag);
+ // bool mcbmFlag = geoTag.Contains("mcbm", TString::kIgnoreCase);
+ // CbmMuchGeoScheme::Instance()->Init(stations, mcbmFlag);
+}
+// -------------------------------------------------------------------------
+
+// -------------------------------------------------------------------------x
+void CbmMuchDigitizerQa::Exec(Option_t*) {
+
+ fNevents++;
+ LOG(info) << "Event: " << fNevents;
+ TDirectory* oldDirectory = gDirectory;
+
+ OccupancyQa();
+ DigitizerQa();
+ DrawCanvases();
+ if (fVerbose > 1) {
+ PrintFrontLayerPoints();
+ PrintFrontLayerDigis();
+ }
+ gDirectory = oldDirectory;
+ return;
+}
+// -------------------------------------------------------------------------
+
+void CbmMuchDigitizerQa::PrintFrontLayerPoints() {
+ for (int i = 0; i < fPoints->GetEntriesFast(); i++) {
+ CbmMuchPoint* point = (CbmMuchPoint*) fPoints->At(i);
+ Int_t stId = CbmMuchAddress::GetStationIndex(point->GetDetectorID());
+ if (stId != 0) continue;
+ Int_t layerId = CbmMuchAddress::GetLayerIndex(point->GetDetectorID());
+ if (layerId != 0) continue;
+ printf("point %4i xin=%6.2f yin=%6.2f xout=%6.2f yout=%6.2f zin=%6.2f\n",
+ i,
+ point->GetXIn(),
+ point->GetYIn(),
+ point->GetXOut(),
+ point->GetYOut(),
+ point->GetZIn());
+ }
+}
+// -------------------------------------------------------------------------
+
+void CbmMuchDigitizerQa::PrintFrontLayerDigis() {
+ for (Int_t i = 0; i < fDigiManager->GetNofDigis(ECbmModuleId::kMuch); i++) {
+ CbmMuchDigi* digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(i);
+ UInt_t address = digi->GetAddress();
+ Int_t stId = CbmMuchAddress::GetStationIndex(address);
+ if (stId != 0) continue;
+ Int_t layerId = CbmMuchAddress::GetLayerIndex(address);
+ if (layerId != 0) continue;
+ CbmMuchModuleGem* module =
+ (CbmMuchModuleGem*) fGeoScheme->GetModuleByDetId(address);
+ if (!module) continue;
+ CbmMuchPad* pad = module->GetPad(address);
+ Double_t x0 = pad->GetX();
+ Double_t y0 = pad->GetY();
+ UInt_t charge = digi->GetAdc();
+ printf("digi %4i x0=%5.1f y0=%5.1f charge=%3i\n", i, x0, y0, charge);
+ }
+}
+// -------------------------------------------------------------------------
+
+void CbmMuchDigitizerQa::DrawCanvases() {
+
+ fCanvCharge->cd(1);
+ fhMcPointCharge->DrawCopy("", "");
+ fCanvCharge->cd(2);
+ fhMcPointChargeLog->DrawCopy("", "");
+ fCanvCharge->cd(3);
+ fhMcPointChargePr_1GeV_3mm->DrawCopy("", "");
+
+ for (Int_t i = 0; i < 4; i++) {
+ fCanvTrackLength->cd(i + 1);
+ gPad->SetLogy();
+ gStyle->SetOptStat(1110);
+ }
+ fCanvTrackLength->cd(1);
+ fhTrackLength->DrawCopy("", "");
+ fCanvTrackLength->cd(2);
+ fhTrackLengthPi->DrawCopy("", "");
+ fCanvTrackLength->cd(3);
+ fhTrackLengthPr->DrawCopy("", "");
+ fCanvTrackLength->cd(4);
+ fhTrackLengthEl->DrawCopy("", "");
+
+ for (Int_t i = 0; i < 4; i++) {
+ fCanvChargeVsEnergy->cd(i + 1);
+ gPad->Range(0, 0, 200, 200);
+ gPad->SetBottomMargin(0.11);
+ gPad->SetRightMargin(0.14);
+ gPad->SetLeftMargin(0.12);
+ gPad->SetLogz();
+ }
+ fCanvChargeVsEnergy->cd(1);
+ fhMcPointChargeVsTrackEnergyLog->DrawCopy("colz", "");
+ fCanvChargeVsEnergy->cd(2);
+ fhMcPointChargeVsTrackEnergyLogPi->DrawCopy("colz", "");
+ fFitPi->DrawCopy("same");
+ fCanvChargeVsEnergy->cd(3);
+ fhMcPointChargeVsTrackEnergyLogPr->DrawCopy("colz", "");
+ fFitPr->DrawCopy("same");
+ fCanvChargeVsEnergy->cd(4);
+ fhMcPointChargeVsTrackEnergyLogEl->DrawCopy("colz", "");
+ fFitEl->DrawCopy("same");
+
+ for (Int_t i = 0; i < 4; i++) {
+ fCanvChargeVsLength->cd(i + 1);
+ gPad->Range(0, 0, 200, 200);
+ gPad->SetBottomMargin(0.11);
+ gPad->SetRightMargin(0.14);
+ gPad->SetLeftMargin(0.12);
+ gPad->SetLogz();
+ }
+ fCanvChargeVsLength->cd(1);
+ fhMcPointChargeVsTrackLength->DrawCopy("colz", "");
+ fCanvChargeVsLength->cd(2);
+ fhMcPointChargeVsTrackLengthPi->DrawCopy("colz", "");
+ fCanvChargeVsLength->cd(3);
+ fhMcPointChargeVsTrackLengthPr->DrawCopy("colz", "");
+ fCanvChargeVsLength->cd(4);
+ fhMcPointChargeVsTrackLengthEl->DrawCopy("colz", "");
+
+ for (Int_t i = 0; i < fNstations; i++) {
+ *fvPadsFiredR[i] = *fvUsPadsFiredR[i];
+ //fvPadsFiredR[i]->Sumw2();
+ fvPadsFiredR[i]->Scale(1. / fNevents);
+ fvPadOccupancyR[i]->Divide(fvPadsFiredR[i], fvPadsTotalR[i]);
+ fvPadOccupancyR[i]->Scale(100.);
+
+ fCanvStationCharge->cd(i + 1);
+ fvMcPointCharge[i]->DrawCopy("", "");
+ fCanvPadOccupancyR->cd(i + 1);
+ fvPadOccupancyR[i]->DrawCopy("", "");
+ fCanvPadsTotalR->cd(i + 1);
+ fvPadsTotalR[i]->DrawCopy("", "");
+ fCanvUsPadsFiredXY->cd(i + 1);
+ fvUsPadsFiredXY[i]->DrawCopy("colz", "");
+ }
+ fCanvNpadsVsArea->cd();
+ fhNpadsVsS->DrawCopy("colz", "");
+}
+
+// -------------------------------------------------------------------------
+void CbmMuchDigitizerQa::FinishTask() {
+
+ printf("FinishTask\n");
+ cout << "\n\n SG: Finish task!" << endl;
+
+ TDirectory* oldDirectory = gDirectory;
+ bool oldBatchMode = gROOT->IsBatch();
+
+ if (!FairRootManager::Instance() || !FairRootManager::Instance()->GetSink()) {
+ cout << "No sink found" << endl;
+ return;
+ }
+ // gROOT->SetBatch(kTRUE);
+ gROOT->SetBatch(kFALSE);
+ gStyle->SetPaperSize(20, 20);
+
+ if (fVerbose > 1) {
+ OutputNvsS();
+ PrintCounters();
+ }
+ FairSink* sink = FairRootManager::Instance()->GetSink();
+ sink->WriteObject(&fOutFolder, nullptr);
+ gDirectory = oldDirectory;
+ gROOT->SetBatch(oldBatchMode);
+}
+// -------------------------------------------------------------------------
+
+void CbmMuchDigitizerQa::OutputNvsS() {
+ TCanvas* c = new TCanvas("nMeanVsS", "nMeanVsS", 2 * 800, 2 * 400);
+ printf("===================================\n");
+ printf("DigitizerQa:\n");
+
+ Double_t nMean[11];
+ Double_t s[11];
+ for (Int_t iS = 1; iS <= 10; iS++) {
+ nMean[iS] = 0;
+ s[iS] = -5.25 + 0.5 * iS;
+ Double_t total = 0;
+ for (Int_t iN = 1; iN <= 10; iN++) {
+ nMean[iS] += iN * fhNpadsVsS->GetBinContent(iS, iN);
+ total += fhNpadsVsS->GetBinContent(iS, iN);
+ }
+ if (total > 0) nMean[iS] /= total;
+ printf("%f %f\n", s[iS], nMean[iS]);
+ }
+ c->cd();
+
+ TGraph* gNvsS = new TGraph(11, s, nMean);
+ //gNvsS->SetLineColor(2);
+ //gNvsS->SetLineWidth(4);
+ gNvsS->SetMarkerColor(4);
+ gNvsS->SetMarkerSize(1.5);
+ gNvsS->SetMarkerStyle(21);
+ gNvsS->SetTitle("nMeanVsS");
+ gNvsS->GetYaxis()->SetTitle("nMean");
+ gNvsS->GetYaxis()->SetTitle("nMean");
+ //gNvsS->DrawClone("ALP");
+ gNvsS->DrawClone("AP");
+ fOutFolder.Add(c);
+}
+
+void CbmMuchDigitizerQa::PrintCounters() {
+
+ printf("All tracks: ;");
+ for (Int_t i = 0; i < fNstations; i++)
+ printf("%8i;", fNall[i]);
+ printf("\n");
+ printf("------------;");
+ for (Int_t i = 0; i < fNstations; i++)
+ printf("---------");
+ printf("\n");
+ printf("Primary: ;");
+ for (Int_t i = 0; i < fNstations; i++)
+ printf("%8i;", fNprimary[i]);
+ printf("\n");
+ printf("Secondary: ;");
+ for (Int_t i = 0; i < fNstations; i++)
+ printf("%8i;", fNsecondary[i]);
+ printf("\n");
+ printf("-------------");
+ for (Int_t i = 0; i < fNstations; i++)
+ printf("---------");
+ printf("\n");
+ printf("Protons: ;");
+ for (Int_t i = 0; i < fNstations; i++)
+ printf("%8i;", fNpr[i]);
+ printf("\n");
+ printf("Pions: ;");
+ for (Int_t i = 0; i < fNstations; i++)
+ printf("%8i;", fNpi[i]);
+ printf("\n");
+ printf("Electrons: ;");
+ for (Int_t i = 0; i < fNstations; i++)
+ printf("%8i;", fNel[i]);
+ printf("\n");
+ printf("Muons: ;");
+ for (Int_t i = 0; i < fNstations; i++)
+ printf("%8i;", fNmu[i]);
+ printf("\n");
+ printf("Kaons: ;");
+ for (Int_t i = 0; i < fNstations; i++)
+ printf("%8i;", fNka[i]);
+ printf("\n");
+}
+
+// -------------------------------------------------------------------------
+Double_t CbmMuchDigitizerQa::LandauMPV(Double_t* lg_x, Double_t* par) {
+ Double_t gaz_gain_mean = 1.7e+4;
+ Double_t scale = 1.e+6;
+ gaz_gain_mean /= scale;
+ Double_t mass = par[0]; // mass in MeV
+ Double_t x = TMath::Power(10, lg_x[0]);
+ return gaz_gain_mean * MPV_n_e(x, mass);
+}
+// -------------------------------------------------------------------------
+
+// -------------------------------------------------------------------------
+void CbmMuchDigitizerQa::DigitizerQa() {
+ TVector3 vIn; // in position of the track
+ TVector3 vOut; // out position of the track
+ TVector3 p; // track momentum
+ fPointInfos->Clear();
+
+ for (Int_t i = 0; i < fPoints->GetEntriesFast(); i++) {
+ CbmMuchPoint* point = (CbmMuchPoint*) fPoints->At(i);
+ Int_t stId = CbmMuchAddress::GetStationIndex(point->GetDetectorID());
+
+ // Check if the point corresponds to a certain MC Track
+ Int_t trackId = point->GetTrackID();
+ if (trackId < 0) {
+ new ((*fPointInfos)[i]) CbmMuchPointInfo(0, 0, 0, 0, 0);
+ continue;
+ }
+ CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(trackId);
+ if (!mcTrack) {
+ new ((*fPointInfos)[i]) CbmMuchPointInfo(0, 0, 0, 0, 0);
+ continue;
+ }
+
+ Int_t motherId = mcTrack->GetMotherId();
+ Int_t pdgCode = mcTrack->GetPdgCode();
+ TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pdgCode);
+ if (!particle || pdgCode == 22 || // photons
+ pdgCode == 2112) // neutrons
+ {
+ new ((*fPointInfos)[i]) CbmMuchPointInfo(0, 0, 0, 0, 0);
+ continue;
+ }
+
+ if (CbmMuchAddress::GetLayerIndex(point->GetDetectorID()) == 0) {
+ UpdateParticleCounters(stId, pdgCode, motherId);
+ }
+
+ Double_t mass = particle->Mass();
+ point->PositionIn(vIn);
+ point->PositionOut(vOut);
+ point->Momentum(p);
+ Double_t length = (vOut - vIn).Mag(); // Track length
+ Double_t kine = sqrt(p.Mag2() + mass * mass) - mass; // Kinetic energy
+ // Get mother pdg code
+ Int_t motherPdg = 0;
+ CbmMCTrack* motherTrack = NULL;
+ if (motherId != -1) motherTrack = (CbmMCTrack*) fMCTracks->At(motherId);
+ if (motherTrack) motherPdg = motherTrack->GetPdgCode();
+ new ((*fPointInfos)[i])
+ CbmMuchPointInfo(pdgCode, motherPdg, kine, length, stId);
+ }
+ FillChargePerPoint();
+ FillDigitizerPerformancePlots();
+}
+// -------------------------------------------------------------------------
+
+void CbmMuchDigitizerQa::UpdateParticleCounters(Int_t stId,
+ Int_t pdgCode,
+ Int_t motherId) {
+
+ fNall[stId]++;
+ if (pdgCode == 2212)
+ fNpr[stId]++;
+ else if (pdgCode == -211 || pdgCode == 211)
+ fNpi[stId]++;
+ else if (pdgCode == -11 || pdgCode == 11)
+ fNel[stId]++;
+ else if (pdgCode == -13 || pdgCode == 13)
+ fNmu[stId]++;
+ else if (pdgCode == -321 || pdgCode == 321)
+ fNka[stId]++;
+
+ if (motherId == -1)
+ fNprimary[stId]++;
+ else
+ fNsecondary[stId]++;
+}
+
+void CbmMuchDigitizerQa::FillChargePerPoint() {
+ for (Int_t i = 0; i < fDigiManager->GetNofDigis(ECbmModuleId::kMuch); i++) {
+ CbmMuchDigi* digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(i);
+ assert(digi);
+ CbmMatch* match =
+ (CbmMatch*) fDigiManager->GetMatch(ECbmModuleId::kMuch, i);
+ assert(match);
+ CbmMuchModule* module = fGeoScheme->GetModuleByDetId(digi->GetAddress());
+ assert(module);
+ if (!module) continue;
+ Double_t area = 0;
+ if (module->GetDetectorType() != 4 && module->GetDetectorType() != 3)
+ continue; /// rpc
+ LOG(debug) << GetName() << " Processing MuchDigi " << i << " at address "
+ << digi->GetAddress() << " Module number "
+ << module->GetDetectorType();
+
+ CbmMuchModuleGem* module1 = (CbmMuchModuleGem*) module;
+ assert(module1);
+ CbmMuchPad* pad = module1->GetPad(digi->GetAddress());
+ assert(pad);
+ area = pad->GetDx() * pad->GetDy();
+ Int_t nofLinks = match->GetNofLinks();
+ for (Int_t pt = 0; pt < nofLinks; pt++) {
+ Int_t pointId = match->GetLink(pt).GetIndex();
+ Int_t charge = match->GetLink(pt).GetWeight();
+ CbmMuchPointInfo* info = (CbmMuchPointInfo*) fPointInfos->At(pointId);
+ if (info->GetPdgCode() == 0) continue;
+ info->AddCharge(charge);
+ info->AddArea(area);
+ }
+ }
+}
+
+void CbmMuchDigitizerQa::FillDigitizerPerformancePlots() {
+ Bool_t verbose = (fVerbose > 2);
+ for (Int_t i = 0; i < fPointInfos->GetEntriesFast(); i++) {
+ CbmMuchPointInfo* info = (CbmMuchPointInfo*) fPointInfos->At(i);
+ if (verbose) {
+ printf("%i", i);
+ info->Print();
+ }
+ Double_t length = info->GetLength();
+ Double_t kine = 1000 * (info->GetKine());
+ Double_t charge = info->GetCharge();
+ Int_t pdg = info->GetPdgCode();
+ if (pdg == 0) continue;
+ if (charge <= 0) continue;
+ Double_t log_kine = TMath::Log10(kine);
+ Double_t log_charge = TMath::Log10(charge);
+ charge = charge / 1.e+4; // measure charge in 10^4 electrons
+
+ Int_t nPads = info->GetNPads();
+ Double_t area = info->GetArea() / nPads;
+ // printf("%f %i\n",TMath::Log2(area),nPads);
+ fhNpadsVsS->Fill(TMath::Log2(area), nPads);
+ fhMcPointCharge->Fill(charge);
+ fvMcPointCharge[info->GetStationId()]->Fill(charge);
+ fhMcPointChargeLog->Fill(log_charge);
+ fhMcPointChargeVsTrackEnergyLog->Fill(log_kine, charge);
+ fhMcPointChargeVsTrackLength->Fill(length, charge);
+ fhTrackLength->Fill(length);
+
+ if (pdg == 2212) {
+ fhMcPointChargeVsTrackEnergyLogPr->Fill(log_kine, charge);
+ fhMcPointChargeVsTrackLengthPr->Fill(length, charge);
+ fhTrackLengthPr->Fill(length);
+ } else if (pdg == 211 || pdg == -211) {
+ fhMcPointChargeVsTrackEnergyLogPi->Fill(log_kine, charge);
+ fhMcPointChargeVsTrackLengthPi->Fill(length, charge);
+ fhTrackLengthPi->Fill(length);
+ } else if (pdg == 11 || pdg == -11) {
+ fhMcPointChargeVsTrackEnergyLogEl->Fill(log_kine, charge);
+ fhMcPointChargeVsTrackLengthEl->Fill(length, charge);
+ fhTrackLengthEl->Fill(length);
+ }
+ if (pdg == 2212 && length > 0.3 && length < 0.32 && kine > 1000
+ && kine < 1020)
+ fhMcPointChargePr_1GeV_3mm->Fill(charge);
+ }
+}
+
+// -------------------------------------------------------------------------
+void CbmMuchDigitizerQa::OccupancyQa() {
+ // Filling occupancy plots
+ for (Int_t i = 0; i < fDigiManager->GetNofDigis(ECbmModuleId::kMuch); i++) {
+ CbmMuchDigi* digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(i);
+ assert(digi);
+ UInt_t address = digi->GetAddress();
+ CbmMuchModule* module = fGeoScheme->GetModuleByDetId(address);
+ if (!module) continue;
+ Double_t r0 = 0;
+ if (module->GetDetectorType() != 4 && module->GetDetectorType() != 3)
+ continue;
+ CbmMuchModuleGem* module1 = (CbmMuchModuleGem*) module;
+ CbmMuchPad* pad = module1->GetPad(
+ address); // fGeoScheme->GetPadByDetId(detectorId, channelId);
+ assert(pad);
+ Double_t x0 = pad->GetX();
+ Double_t y0 = pad->GetY();
+ r0 = TMath::Sqrt(x0 * x0 + y0 * y0);
+ fvUsPadsFiredR[CbmMuchAddress::GetStationIndex(address)]->Fill(r0);
+ fvUsPadsFiredXY[CbmMuchAddress::GetStationIndex(address)]->Fill(x0, y0);
+ }
+}
+// -------------------------------------------------------------------------
+
+// -------------------------------------------------------------------------
+Int_t CbmMuchDigitizerQa::GetNChannels(Int_t iStation) {
+ Int_t nChannels = 0;
+ vector<CbmMuchModule*> modules = fGeoScheme->GetModules(iStation);
+ for (UInt_t im = 0; im < modules.size(); im++) {
+ CbmMuchModule* mod = modules[im];
+ if (mod->GetDetectorType() != 4 && mod->GetDetectorType() != 3) continue;
+ CbmMuchModuleGem* module = (CbmMuchModuleGem*) mod;
+ if (!module) continue;
+ nChannels += module->GetNPads();
+ }
+ return nChannels;
+}
+
+// -------------------------------------------------------------------------
+Int_t CbmMuchDigitizerQa::GetNSectors(Int_t iStation) {
+ Int_t nSectors = 0;
+ vector<CbmMuchModule*> modules = fGeoScheme->GetModules(iStation);
+ for (UInt_t im = 0; im < modules.size(); im++) {
+ CbmMuchModule* mod = modules[im];
+ if (mod->GetDetectorType() != 4 && mod->GetDetectorType() != 3) continue;
+ CbmMuchModuleGem* module = (CbmMuchModuleGem*) mod;
+ if (!module) continue;
+ nSectors += module->GetNSectors();
+ }
+ return nSectors;
+}
+
+// -------------------------------------------------------------------------
+TVector2 CbmMuchDigitizerQa::GetMinPadSize(Int_t iStation) {
+ Double_t padMinLx = std::numeric_limits<Double_t>::max();
+ Double_t padMinLy = std::numeric_limits<Double_t>::max();
+ vector<CbmMuchModule*> modules = fGeoScheme->GetModules(iStation);
+ for (UInt_t im = 0; im < modules.size(); im++) {
+ CbmMuchModule* mod = modules[im];
+ if (mod->GetDetectorType() != 1) continue;
+ CbmMuchModuleGem* module = (CbmMuchModuleGem*) mod;
+ vector<CbmMuchPad*> pads = module->GetPads();
+ for (UInt_t ip = 0; ip < pads.size(); ip++) {
+ CbmMuchPad* pad = pads[ip];
+ if (pad->GetDx() < padMinLx) padMinLx = pad->GetDx();
+ if (pad->GetDy() < padMinLy) padMinLy = pad->GetDy();
+ }
+ }
+ return TVector2(padMinLx, padMinLy);
+}
+// -------------------------------------------------------------------------
+
+// -------------------------------------------------------------------------
+TVector2 CbmMuchDigitizerQa::GetMaxPadSize(Int_t iStation) {
+ Double_t padMaxLx = std::numeric_limits<Double_t>::min();
+ Double_t padMaxLy = std::numeric_limits<Double_t>::min();
+ vector<CbmMuchModule*> modules = fGeoScheme->GetModules(iStation);
+ for (UInt_t im = 0; im < modules.size(); im++) {
+ CbmMuchModule* mod = modules[im];
+ if (mod->GetDetectorType() != 1) continue;
+ CbmMuchModuleGem* module = (CbmMuchModuleGem*) mod;
+ vector<CbmMuchPad*> pads = module->GetPads();
+ for (UInt_t ip = 0; ip < pads.size(); ip++) {
+ CbmMuchPad* pad = pads[ip];
+ if (pad->GetDx() > padMaxLx) padMaxLx = pad->GetDx();
+ if (pad->GetDy() > padMaxLy) padMaxLy = pad->GetDy();
+ }
+ }
+ return TVector2(padMaxLx, padMaxLy);
+}
+// -------------------------------------------------------------------------
+
+// -------------------------------------------------------------------------
+Double_t CbmMuchDigitizerQa::MPV_n_e(Double_t Tkin, Double_t mass) {
+ Double_t logT;
+ TF1 fPol6("fPol6", "pol6", -5, 10);
+ if (mass < 0.1) {
+ logT = log(Tkin * 0.511 / mass);
+ if (logT > 9.21034) logT = 9.21034;
+ if (logT < min_logT_e) logT = min_logT_e;
+ return fPol6.EvalPar(&logT, mpv_e);
+ } else if (mass >= 0.1 && mass < 0.2) {
+ logT = log(Tkin * 105.658 / mass);
+ if (logT > 9.21034) logT = 9.21034;
+ if (logT < min_logT_mu) logT = min_logT_mu;
+ return fPol6.EvalPar(&logT, mpv_mu);
+ } else {
+ logT = log(Tkin * 938.272 / mass);
+ if (logT > 9.21034) logT = 9.21034;
+ if (logT < min_logT_p) logT = min_logT_p;
+ return fPol6.EvalPar(&logT, mpv_p);
+ }
+}
+
+void CbmMuchDigitizerQa::DivideCanvas2D(TCanvas* c, int nPads) {
+ // divide canvas into nPads in 2D
+ if (!c || nPads < 1) { return; }
+ int rows = (int) sqrt(nPads);
+ int cols = nPads / rows;
+ if (cols * rows < nPads) { cols++; }
+ c->Divide(cols, rows);
+}
+
+// -------------------------------------------------------------------------
+
+ClassImp(CbmMuchDigitizerQa)
diff --git a/sim/detectors/much/qa/CbmMuchDigitizerQa.h b/sim/detectors/much/qa/CbmMuchDigitizerQa.h
new file mode 100644
index 0000000000..aad04643d6
--- /dev/null
+++ b/sim/detectors/much/qa/CbmMuchDigitizerQa.h
@@ -0,0 +1,199 @@
+/// \file CbmMuchDigitizerQa.h
+/// \brief Definition of the CbmMuchDigitizerQa class
+/// \author Sergey Gorbunov <se.gorbunov@gsi.de>
+/// \author Eugeny Kryshen
+/// \author Vikas Singhal
+/// \author Ekata Nandy
+/// \date 25.09.2020
+
+#ifndef CbmMuchDigitizerQa_H
+#define CbmMuchDigitizerQa_H
+
+#include "CbmDigiManager.h"
+#include "CbmGeoMuchPar.h"
+#include "CbmMCTrack.h"
+#include "CbmMatch.h"
+#include "CbmMuchAddress.h"
+#include "CbmMuchDigi.h"
+#include "CbmMuchGeoScheme.h"
+#include "CbmMuchModuleGem.h"
+#include "CbmMuchPad.h"
+#include "CbmMuchPoint.h"
+#include "CbmMuchPointInfo.h"
+#include "CbmMuchRecoDefs.h"
+#include "CbmMuchSector.h"
+#include "CbmMuchStation.h"
+#include "CbmQaCanvas.h"
+#include "FairLogger.h"
+#include "FairRootFileSink.h"
+#include "FairRootManager.h"
+#include "FairRun.h"
+#include "FairRuntimeDb.h"
+#include "FairTask.h"
+#include "Riostream.h"
+#include "TArrayI.h"
+#include "TCanvas.h"
+#include "TClonesArray.h"
+#include "TDatabasePDG.h"
+#include "TF1.h"
+#include "TFile.h"
+#include "TGraph.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TObjArray.h"
+#include "TParticlePDG.h"
+#include "TString.h"
+#include "TStyle.h"
+#include <algorithm>
+#include <cassert>
+#include <map>
+#include <vector>
+
+class CbmMuchGeoScheme;
+
+/// QA for the MUCH detector after a "digitization" step of the simulation.
+/// The class reimplements corresponding QA checks from old CbmMuchHitFinderQa class
+/// made by E. Kryshen & V. Singhal & E. Nandy
+///
+class CbmMuchDigitizerQa : public FairTask {
+
+public:
+ CbmMuchDigitizerQa(const char* name = "MuchHitFinderQa", Int_t verbose = 1);
+ virtual ~CbmMuchDigitizerQa();
+ virtual InitStatus Init();
+ virtual void Exec(Option_t* option);
+ virtual void FinishTask();
+ virtual void SetParContainers();
+ static void DivideCanvas2D(TCanvas* c, int nPads);
+
+protected:
+ /* DigitizerQa - analysis of digitizer performance - charge distributions
+ * for tracks. Track length distrivutions. Statistics on particle types
+ */
+ void DigitizerQa();
+
+ /* Occupance analysis - all pads,fired pads,
+ * and fired/all distributions as functions of radius
+ */
+ void OccupancyQa();
+
+private:
+ Int_t GetNChannels(Int_t iStation);
+ Int_t GetNSectors(Int_t iStation);
+ TVector2 GetMinPadSize(Int_t iStation);
+ TVector2 GetMaxPadSize(Int_t iStation);
+ static Double_t LandauMPV(Double_t* x, Double_t* par);
+ static Double_t MPV_n_e(Double_t Tkin, Double_t mass);
+
+ CbmMuchDigitizerQa(const CbmMuchDigitizerQa&);
+ CbmMuchDigitizerQa& operator=(const CbmMuchDigitizerQa&);
+
+ TFolder* histFolder;
+ void InitCounters();
+ void InitCanvases();
+ void InitChargeHistos();
+ void InitLengthHistos();
+ void InitPadHistos();
+ void InitChannelPadInfo();
+ void InitFits();
+ void DeInit();
+
+ void FillTotalPadsHistos();
+ void DrawCanvases();
+ void PrintCounters();
+ void OutputNvsS();
+
+ void UpdateParticleCounters(Int_t stId, Int_t pdgCode, Int_t motherId);
+ void FillChargePerPoint();
+ void FillDigitizerPerformancePlots();
+
+ void PrintFrontLayerPoints();
+ void PrintFrontLayerDigis();
+
+ // geometry
+ CbmMuchGeoScheme* fGeoScheme = nullptr;
+ Int_t fNstations = 0;
+ CbmDigiManager* fDigiManager = nullptr;
+
+ // containers
+ TClonesArray* fPoints = nullptr;
+ TClonesArray* fDigis = nullptr;
+ TClonesArray* fDigiMatches = nullptr;
+ TClonesArray* fMCTracks = nullptr;
+ TClonesArray* fPointInfos = nullptr; /// temporary additional information
+
+ TFolder fOutFolder; /// output folder with histos and canvases
+ Int_t fNevents = 0; /// number of processed events
+
+ // internal unscaled histograms, need to be scaled at the output
+ std::vector<TH1F*> fvUsPadsFiredR; // fired pads vs R, per station
+ std::vector<TH2F*> fvUsPadsFiredXY; // fired pads vs XY, per station
+
+ // output histograms
+ TH1F* fhMcPointCharge = nullptr; /// MC point charge
+ TH1F* fhMcPointChargeLog = nullptr; /// MC point charge log scale
+ TH1F* fhMcPointChargePr_1GeV_3mm =
+ nullptr; /// MC point charge for selected protons
+
+ TH1F* fhTrackLength = nullptr;
+ TH1F* fhTrackLengthPi = nullptr;
+ TH1F* fhTrackLengthPr = nullptr;
+ TH1F* fhTrackLengthEl = nullptr;
+
+ TH2F* fhMcPointChargeVsTrackEnergyLog = nullptr;
+ TH2F* fhMcPointChargeVsTrackEnergyLogPi = nullptr;
+ TH2F* fhMcPointChargeVsTrackEnergyLogPr = nullptr;
+ TH2F* fhMcPointChargeVsTrackEnergyLogEl = nullptr;
+ TH2F* fhMcPointChargeVsTrackLength = nullptr;
+ TH2F* fhMcPointChargeVsTrackLengthPi = nullptr;
+ TH2F* fhMcPointChargeVsTrackLengthPr = nullptr;
+ TH2F* fhMcPointChargeVsTrackLengthEl = nullptr;
+ TH2F* fhNpadsVsS = nullptr;
+
+ std::vector<TH1F*> fvMcPointCharge; // MC point charge per station
+ std::vector<TH1F*> fvPadsTotalR; // number of pads vs R, per station
+ std::vector<TH1F*> fvPadsFiredR; // fired pads vs R, per station
+ std::vector<TH1F*> fvPadOccupancyR; // pad occupancy vs R, per station
+
+ // output canvaces with histogramm collections
+ CbmQaCanvas* fCanvCharge = nullptr;
+ CbmQaCanvas* fCanvStationCharge = nullptr;
+ CbmQaCanvas* fCanvChargeVsEnergy = nullptr;
+ CbmQaCanvas* fCanvChargeVsLength = nullptr;
+ CbmQaCanvas* fCanvTrackLength = nullptr;
+ CbmQaCanvas* fCanvNpadsVsArea = nullptr;
+ CbmQaCanvas* fCanvUsPadsFiredXY = nullptr;
+ CbmQaCanvas* fCanvPadOccupancyR = nullptr;
+ CbmQaCanvas* fCanvPadsTotalR = nullptr;
+
+ TF1* fFitEl = nullptr;
+ TF1* fFitPi = nullptr;
+ TF1* fFitPr = nullptr;
+
+ Int_t fSignalPoints = 0; // Number of signal MC points
+ Int_t fnPadSizesX = 0;
+ Int_t fnPadSizesY = 0;
+
+ Int_t* fNall = nullptr; // number of all tracks at the first station
+ Int_t* fNpr = nullptr; // number of protons at the first station
+ Int_t* fNpi = nullptr; // number of pions at the first station
+ Int_t* fNel = nullptr; // number of electrons at the first station
+ Int_t* fNmu = nullptr; // number of muons at the first station
+ Int_t* fNka = nullptr; // number of kaons at the first station
+ Int_t* fNprimary = nullptr; // number of primary tracks at the first station
+ Int_t* fNsecondary =
+ nullptr; // number of secondary tracks at the first station
+
+ Int_t fPointsTotal = 0;
+ Int_t fPointsUnderCounted = 0;
+ Int_t fPointsOverCounted = 0;
+
+ Double_t fPadMinLx = 0.;
+ Double_t fPadMinLy = 0.;
+ Double_t fPadMaxLx = 0.;
+ Double_t fPadMaxLy = 0.;
+
+ ClassDef(CbmMuchDigitizerQa, 0)
+};
+
+#endif
diff --git a/reco/detectors/much/CbmMuchPointInfo.cxx b/sim/detectors/much/qa/CbmMuchPointInfo.cxx
similarity index 100%
rename from reco/detectors/much/CbmMuchPointInfo.cxx
rename to sim/detectors/much/qa/CbmMuchPointInfo.cxx
diff --git a/reco/detectors/much/CbmMuchPointInfo.h b/sim/detectors/much/qa/CbmMuchPointInfo.h
similarity index 100%
rename from reco/detectors/much/CbmMuchPointInfo.h
rename to sim/detectors/much/qa/CbmMuchPointInfo.h
diff --git a/sim/detectors/much/qa/CbmMuchTransportQa.cxx b/sim/detectors/much/qa/CbmMuchTransportQa.cxx
index 772ea6349a..1db69c4ca2 100644
--- a/sim/detectors/much/qa/CbmMuchTransportQa.cxx
+++ b/sim/detectors/much/qa/CbmMuchTransportQa.cxx
@@ -4,41 +4,10 @@
/// \author Eugeny Kryshen
/// \author Vikas Singhal
/// \author Ekata Nandy
-/// \date 25.09.2020
+/// \author Dominik Smith
+/// \date 21.10.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 "TLegend.h"
-#include "TPie.h"
-#include "TPieSlice.h"
-#include "TStyle.h"
-
-#include <cassert>
-#include <vector>
-
ClassImp(CbmMuchTransportQa);
// -------------------------------------------------------------------------
@@ -55,18 +24,15 @@ CbmMuchTransportQa::CbmMuchTransportQa(const char* name, Int_t verbose)
, fvMcPointPrimRatio() {}
// -------------------------------------------------------------------------
-
// -------------------------------------------------------------------------
CbmMuchTransportQa::~CbmMuchTransportQa() { DeInit(); }
// -------------------------------------------------------------------------
-
// -------------------------------------------------------------------------
void CbmMuchTransportQa::DeInit() {
fPoints = nullptr;
fMcTracks = nullptr;
-
fOutFolder.Clear();
fhNevents.SetVal(0);
@@ -82,14 +48,14 @@ void CbmMuchTransportQa::DeInit() {
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]);
}
+ fvMcPointXY.clear();
+ fvMcPointPhiZ.clear();
fvMcPointRZ.clear();
for (uint i = 0; i < fvMcPointPRatio.size(); i++) {
@@ -98,7 +64,6 @@ void CbmMuchTransportQa::DeInit() {
for (uint i = 0; i < fvMcPointPrimRatio.size(); i++) {
SafeDelete(fvMcPointPrimRatio[i]);
}
-
SafeDelete(fhNtracks);
SafeDelete(fhFractionPrim);
SafeDelete(fhFractionSec);
@@ -110,14 +75,13 @@ void CbmMuchTransportQa::DeInit() {
fvUsNtra.clear();
fvFraction.clear();
-
fvMcPointPRatio.clear();
fvMcPointPrimRatio.clear();
SafeDelete(fCanvStationXY);
SafeDelete(fCanvStationPhiZ);
SafeDelete(fCanvStationRZ);
-
+ SafeDelete(fCanvUsNtra);
SafeDelete(fCanvStationPRatio);
SafeDelete(fCanvStationPrimRatio);
@@ -131,109 +95,118 @@ void CbmMuchTransportQa::DeInit() {
InitStatus CbmMuchTransportQa::Init() {
TDirectory* oldDirectory = gDirectory;
-
FairRootManager* manager = FairRootManager::Instance();
+ fMcTracks = (TClonesArray*) manager->GetObject("MCTrack");
+ fPoints = (TClonesArray*) manager->GetObject("MuchPoint");
+ fNstations = CbmMuchGeoScheme::Instance()->GetNStations();
+ histFolder = fOutFolder.AddFolder("hist", "Histogramms");
+
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");
-
+ 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;
+ }
+ }
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);
- }
+ InitCountingHistos();
+ InitFractionHistos();
+ Init2dSpatialDistributionHistos();
+ InitRatioPieCharts();
+ InitCanvases();
+
+ gDirectory = oldDirectory;
+ return kSUCCESS;
+}
+// -------------------------------------------------------------------------
+
+void CbmMuchTransportQa::InitCountingHistos() {
+
+ 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);
}
+}
+
+void CbmMuchTransportQa::InitFractionHistos() {
- {
- 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");
+ 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");
+ if (i == 0) {
+ h->GetYaxis()->SetTitle("N tracks");
+ } else {
h->GetYaxis()->SetTitle("%");
- histFolder->Add(h);
}
+ histFolder->Add(h);
}
+}
+
+void CbmMuchTransportQa::Init2dSpatialDistributionHistos() {
fvMcPointXY.resize(fNstations);
fvMcPointPhiZ.resize(fNstations);
fvMcPointRZ.resize(fNstations);
- fvMcPointPRatio.resize(fNstations);
- fvMcPointPrimRatio.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();
+ 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),
@@ -243,8 +216,6 @@ InitStatus CbmMuchTransportQa::Init() {
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),
@@ -254,7 +225,6 @@ InitStatus CbmMuchTransportQa::Init() {
100,
-200.,
200.);
- histFolder->Add(fvMcPointPhiZ[i]);
float dR = rMax - rMin;
fvMcPointRZ[i] = new TH2F(Form("hMcPointRZ%i", i + 1),
@@ -265,42 +235,53 @@ InitStatus CbmMuchTransportQa::Init() {
100,
rMin - 0.1 * dR,
rMax + 0.1 * dR);
+ histFolder->Add(fvMcPointXY[i]);
+ histFolder->Add(fvMcPointPhiZ[i]);
histFolder->Add(fvMcPointRZ[i]);
+ }
+}
+
+void CbmMuchTransportQa::InitRatioPieCharts() {
+ fvMcPointPRatio.resize(fNstations);
+ fvMcPointPrimRatio.resize(fNstations);
+ for (Int_t i = 0; i < fNstations; i++) {
fvMcPointPRatio[i] =
new TPie(Form("fvMcPointPRatio%i", i + 1),
Form("McPoint Particle Ratios: Station %i", i + 1),
5);
- //histFolder->Add(fvMcPointPRatio[i]);
-
fvMcPointPrimRatio[i] =
new TPie(Form("fvMcPointPrimRatio%i", i + 1),
Form("McPoint Primary/Secondary Track: Station %i", i + 1),
2);
- //histFolder->Add(fvMcPointPrimRatio[i]);
+ histFolder->Add(fvMcPointPRatio[i]);
+ histFolder->Add(fvMcPointPrimRatio[i]);
}
+}
+
+void CbmMuchTransportQa::InitCanvases() {
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);
+
+ fCanvUsNtra =
+ new CbmQaCanvas("cUsNtra", "Much: MC unscaled counts", 3 * 400, 3 * 400);
+ fCanvUsNtra->Divide2D(9);
fCanvStationPRatio = new CbmQaCanvas(
"cMcPointPRatios", "Much: MC particle ratios", 2 * 400, 2 * 400);
fCanvStationPRatio->Divide2D(fNstations);
- fOutFolder.Add(fCanvStationPRatio);
fCanvStationPrimRatio =
new CbmQaCanvas("cMcPointPrimRatios",
@@ -308,14 +289,14 @@ InitStatus CbmMuchTransportQa::Init() {
2 * 400,
2 * 400);
fCanvStationPrimRatio->Divide2D(fNstations);
- fOutFolder.Add(fCanvStationPrimRatio);
- gDirectory = oldDirectory;
-
- return kSUCCESS;
+ fOutFolder.Add(fCanvStationXY);
+ fOutFolder.Add(fCanvStationPhiZ);
+ fOutFolder.Add(fCanvStationRZ);
+ fOutFolder.Add(fCanvUsNtra);
+ fOutFolder.Add(fCanvStationPRatio);
+ fOutFolder.Add(fCanvStationPrimRatio);
}
-// -------------------------------------------------------------------------
-
// -------------------------------------------------------------------------
InitStatus CbmMuchTransportQa::ReInit() {
@@ -324,7 +305,6 @@ InitStatus CbmMuchTransportQa::ReInit() {
}
// -------------------------------------------------------------------------
-
// -------------------------------------------------------------------------
void CbmMuchTransportQa::SetParContainers() {
// Get run and runtime database
@@ -350,26 +330,20 @@ void CbmMuchTransportQa::SetParContainers() {
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 (!point) {
+ LOG(fatal) << "Much point " << i << " doesn't exist";
+ break;
+ } // Check if the point corresponds to a certain MC Track
if (trackId < 0 || trackId >= fMcTracks->GetEntriesFast()) {
LOG(fatal) << "Much point " << i << ": trackId " << trackId
<< " doesn't belong to [0," << fMcTracks->GetEntriesFast() - 1
@@ -384,12 +358,8 @@ void CbmMuchTransportQa::Exec(Option_t*) {
}
Int_t motherId = mcTrack->GetMotherId();
-
- // Get mass
- Int_t pdgCode = mcTrack->GetPdgCode();
-
+ 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";
@@ -404,77 +374,162 @@ void CbmMuchTransportQa::Exec(Option_t*) {
}
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;
- }
+ FillCountingHistos(stId, motherId, pdgCode);
}
-
trackStaCross[trackId] |= stMask;
+ Fill2dSpatialDistributionHistos(point, stId);
+ }
+}
- 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());
+void CbmMuchTransportQa::FillCountingHistos(Int_t stId,
+ Int_t motherId,
+ Int_t pdgCode) {
+ 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;
}
}
// -------------------------------------------------------------------------
+void CbmMuchTransportQa::Fill2dSpatialDistributionHistos(CbmMuchPoint* point,
+ Int_t stId) {
+
+ 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());
- }
+ 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.);
}
+ MakePRatioPieCharts();
+ MakePrimRatioPieCharts();
+ DrawCanvases();
+
+ gDirectory = oldDirectory;
+ return fOutFolder;
+}
+// -------------------------------------------------------------------------
+
+void CbmMuchTransportQa::DrawCanvases() {
+
+ 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", "");
+
+ fCanvStationPRatio->cd(i + 1);
+ gStyle->SetOptStat(0);
+ fvMcPointPRatio[i]->DrawClone("nol <");
+
+ TLegend* PRatioPieLeg = fvMcPointPRatio[i]->MakeLegend();
+ PRatioPieLeg->SetY1(.56);
+ PRatioPieLeg->SetY2(.86);
+
+ fCanvStationPrimRatio->cd(i + 1);
+ gStyle->SetOptStat(0);
+ fvMcPointPrimRatio[i]->DrawClone("nol <");
+
+ TLegend* PrimRatioPieLeg = fvMcPointPrimRatio[i]->MakeLegend();
+ PrimRatioPieLeg->SetY1(.71);
+ PrimRatioPieLeg->SetY2(.86);
+ PrimRatioPieLeg->SetX1(.40);
+ PrimRatioPieLeg->SetX2(.90);
+
+ gStyle->SetOptStat(1110);
+ }
+
+ fCanvUsNtra->cd(1);
+ gStyle->SetOptStat(0);
+ fhUsNtraAll->DrawCopy("colz", "");
+
+ fCanvUsNtra->cd(2);
+ gStyle->SetOptStat(0);
+ fhNtracks->DrawCopy("colz", "");
+
+ fCanvUsNtra->cd(3);
+ gStyle->SetOptStat(0);
+ fhUsNtraPrim->DrawCopy("colz", "");
+
+ fCanvUsNtra->cd(4);
+ gStyle->SetOptStat(0);
+ fhUsNtraSec->DrawCopy("colz", "");
+
+ fCanvUsNtra->cd(5);
+ gStyle->SetOptStat(0);
+ fhUsNtraPr->DrawCopy("colz", "");
+
+ fCanvUsNtra->cd(6);
+ gStyle->SetOptStat(0);
+ fhUsNtraPi->DrawCopy("colz", "");
+
+ fCanvUsNtra->cd(7);
+ gStyle->SetOptStat(0);
+ fhUsNtraEl->DrawCopy("colz", "");
+
+ fCanvUsNtra->cd(8);
+ gStyle->SetOptStat(0);
+ fhUsNtraMu->DrawCopy("colz", "");
+
+ fCanvUsNtra->cd(9);
+ gStyle->SetOptStat(0);
+ fhUsNtraKa->DrawCopy("colz", "");
+}
+
+void CbmMuchTransportQa::MakePRatioPieCharts() {
for (Int_t i = 0; i < fNstations; i++) {
- Double_t PRatios[] = {fhFractionEl->GetBinContent(i + 1),
+ Double_t PRatios[] = {fhFractionEl->GetBinContent(i + 1),
fhFractionPr->GetBinContent(i + 1),
fhFractionPi->GetBinContent(i + 1),
fhFractionMu->GetBinContent(i + 1),
fhFractionKa->GetBinContent(i + 1)};
-
Int_t PRatiosColors[] = {4, 3, 2, 5, 6};
fvMcPointPRatio[i]->SetEntryVal(0, PRatios[0]);
@@ -497,10 +552,14 @@ TFolder& CbmMuchTransportQa::GetQa() {
fvMcPointPRatio[i]->SetRadius(.33);
fvMcPointPRatio[i]->SetLabelsOffset(-.1);
fvMcPointPRatio[i]->SetLabelFormat("");
+ }
+}
- Double_t PrimRatios[] = {fhFractionPrim->GetBinContent(i + 1),
- fhFractionSec->GetBinContent(i + 1)};
+void CbmMuchTransportQa::MakePrimRatioPieCharts() {
+ for (Int_t i = 0; i < fNstations; i++) {
+ Double_t PrimRatios[] = {fhFractionPrim->GetBinContent(i + 1),
+ fhFractionSec->GetBinContent(i + 1)};
Int_t PrimRatiosColors[] = {6, 4};
fvMcPointPrimRatio[i]->SetEntryVal(0, PrimRatios[0]);
@@ -515,47 +574,7 @@ TFolder& CbmMuchTransportQa::GetQa() {
fvMcPointPrimRatio[i]->SetLabelsOffset(-.1);
fvMcPointPrimRatio[i]->SetLabelFormat("");
}
-
- 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", "");
-
- fCanvStationPRatio->cd(i + 1);
- gStyle->SetOptStat(0);
- fvMcPointPRatio[i]->DrawClone("nol <");
-
- TLegend* PRatioPieLeg = fvMcPointPRatio[i]->MakeLegend();
- PRatioPieLeg->SetY1(.56);
- PRatioPieLeg->SetY2(.86);
-
- fCanvStationPrimRatio->cd(i + 1);
- gStyle->SetOptStat(0);
- fvMcPointPrimRatio[i]->DrawClone("nol <");
-
- TLegend* PrimRatioPieLeg = fvMcPointPrimRatio[i]->MakeLegend();
- PrimRatioPieLeg->SetY1(.71);
- PrimRatioPieLeg->SetY2(.86);
- PrimRatioPieLeg->SetX1(.40);
- PrimRatioPieLeg->SetX2(.90);
-
- gStyle->SetOptStat(1110);
- }
-
- gDirectory = oldDirectory;
- return fOutFolder;
}
-// -------------------------------------------------------------------------
-
// -------------------------------------------------------------------------
void CbmMuchTransportQa::Finish() {
diff --git a/sim/detectors/much/qa/CbmMuchTransportQa.h b/sim/detectors/much/qa/CbmMuchTransportQa.h
index f6232dfd20..3805faa5a1 100644
--- a/sim/detectors/much/qa/CbmMuchTransportQa.h
+++ b/sim/detectors/much/qa/CbmMuchTransportQa.h
@@ -4,26 +4,42 @@
/// \author Eugeny Kryshen
/// \author Vikas Singhal
/// \author Ekata Nandy
-/// \date 25.09.2020
+/// \author Dominik Smith
+/// \date 21.10.2020
#ifndef CbmMuchTransportQa_H
#define CbmMuchTransportQa_H
+#include "CbmGeoMuchPar.h"
+#include "CbmMCTrack.h"
+#include "CbmMuchGeoScheme.h"
+#include "CbmMuchPoint.h"
+#include "CbmMuchStation.h"
+#include "CbmQaCanvas.h"
+#include "FairLogger.h"
+#include "FairRootFileSink.h"
+#include "FairRootManager.h"
+#include "FairRun.h"
+#include "FairRuntimeDb.h"
#include "FairTask.h"
+#include "TClonesArray.h"
+#include "TDatabasePDG.h"
#include "TFolder.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TLegend.h"
#include "TParameter.h"
+#include "TParticlePDG.h"
+#include "TPie.h"
+#include "TPieSlice.h"
+#include "TString.h"
+#include "TStyle.h"
+#include <cassert>
+#include <vector>
-class CbmMuchGeoScheme;
-class TClonesArray;
-class TH1F;
-class TH2F;
-class TPie;
-class TPieSlice;
-class TLegend;
-class CbmQaCanvas;
+#define BINS_STA fNstations, 0, fNstations
/// 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
///
@@ -32,7 +48,6 @@ 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;
@@ -52,10 +67,22 @@ public:
TFolder& GetQa();
private:
+ void InitCountingHistos();
+ void InitFractionHistos();
+ void Init2dSpatialDistributionHistos();
+ void InitRatioPieCharts();
+ void InitCanvases();
+ void FillCountingHistos(Int_t stId, Int_t motherId, Int_t pdgCode);
+ void Fill2dSpatialDistributionHistos(CbmMuchPoint* point, Int_t stId);
+
/// Reset varibles & deallocate memory.
/// When not called by destructor, need to be folloed by Init().
void DeInit();
+ void MakePRatioPieCharts();
+ void MakePrimRatioPieCharts();
+ void DrawCanvases();
+
/// geometry
Int_t fNstations = 0;
@@ -64,11 +91,11 @@ private:
TClonesArray* fMcTracks = nullptr;
///
+ TFolder* histFolder; /// subfolder for histograms
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
@@ -77,11 +104,9 @@ private:
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
+ 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]
@@ -94,24 +119,21 @@ private:
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
+ std::vector<TH1F*> fvFraction; /// pointers to the above histos
/// output pie charts
-
std::vector<TPie*>
fvMcPointPRatio; /// MC point particle ratio pie charts [N stations]
std::vector<TPie*>
fvMcPointPrimRatio; /// MC point particle ratio pie charts [N stations]
// output canvaces with histogramm collections
-
CbmQaCanvas* fCanvStationXY = nullptr;
CbmQaCanvas* fCanvStationPhiZ = nullptr;
CbmQaCanvas* fCanvStationRZ = nullptr;
+ CbmQaCanvas* fCanvUsNtra = nullptr;
// output canvaces with pie chart collections
-
CbmQaCanvas* fCanvStationPRatio = nullptr;
CbmQaCanvas* fCanvStationPrimRatio = nullptr;
--
GitLab