From 6d87bc1645d24e4f1c46b483087d370a6f71279f Mon Sep 17 00:00:00 2001
From: Dominik Smith <d.smith@gsi.de>
Date: Sun, 22 Nov 2020 18:41:52 +0100
Subject: [PATCH] Much reconstruction QA
---
macro/mcbm/mcbm_qa.C | 1 +
macro/run/run_qa.C | 3 +
reco/detectors/much/CMakeLists.txt | 5 +-
reco/detectors/much/CbmMuchHitFinderQa.cxx | 1642 -----------------
reco/detectors/much/CbmMuchHitFinderQa.h | 178 --
reco/detectors/much/qa/CbmMuchHitFinderQa.cxx | 692 +++++++
reco/detectors/much/qa/CbmMuchHitFinderQa.h | 112 ++
sim/detectors/much/qa/CbmMuchDigitizerQa.cxx | 93 +-
sim/detectors/much/qa/CbmMuchDigitizerQa.h | 10 +-
sim/detectors/much/qa/CbmMuchTransportQa.cxx | 38 +-
sim/detectors/much/qa/CbmMuchTransportQa.h | 6 +-
11 files changed, 911 insertions(+), 1869 deletions(-)
delete mode 100644 reco/detectors/much/CbmMuchHitFinderQa.cxx
delete mode 100644 reco/detectors/much/CbmMuchHitFinderQa.h
create mode 100644 reco/detectors/much/qa/CbmMuchHitFinderQa.cxx
create mode 100644 reco/detectors/much/qa/CbmMuchHitFinderQa.h
diff --git a/macro/mcbm/mcbm_qa.C b/macro/mcbm/mcbm_qa.C
index 7a0283d61e..da953d9e82 100644
--- a/macro/mcbm/mcbm_qa.C
+++ b/macro/mcbm/mcbm_qa.C
@@ -153,6 +153,7 @@ void mcbm_qa(Int_t nEvents = 0,
if (CbmSetup::Instance()->IsActive(ECbmModuleId::kMuch)) {
run->AddTask(new CbmMuchTransportQa());
run->AddTask(new CbmMuchDigitizerQa());
+ run->AddTask(new CbmMuchHitFinderQa());
}
// ------------------------------------------------------------------------
diff --git a/macro/run/run_qa.C b/macro/run/run_qa.C
index 35d1c4dcf7..acdd0c755d 100644
--- a/macro/run/run_qa.C
+++ b/macro/run/run_qa.C
@@ -25,6 +25,8 @@
#include "CbmDefs.h"
#include "CbmMCDataManager.h"
+#include "CbmMuchDigitizerQa.h"
+#include "CbmMuchHitFinderQa.h"
#include "CbmMuchTransportQa.h"
#include "CbmSetup.h"
@@ -153,6 +155,7 @@ void run_qa(Int_t nEvents = 0,
if (CbmSetup::Instance()->IsActive(ECbmModuleId::kMuch)) {
run->AddTask(new CbmMuchTransportQa());
run->AddTask(new CbmMuchDigitizerQa());
+ run->AddTask(new CbmMuchHitFinderQa());
}
// ------------------------------------------------------------------------
diff --git a/reco/detectors/much/CMakeLists.txt b/reco/detectors/much/CMakeLists.txt
index 135c25a423..ecddb5b576 100644
--- a/reco/detectors/much/CMakeLists.txt
+++ b/reco/detectors/much/CMakeLists.txt
@@ -4,7 +4,7 @@ ${CMAKE_CURRENT_SOURCE_DIR}
${CBMDETECTORBASE_DIR}/much
${CBMROOT_SOURCE_DIR}/reco/base
-${CBMROOT_SOURCE_DIR}/sim/detectors/much/qa
+${CBMROOT_SOURCE_DIR}/core/qa
${CBMBASE_DIR}
@@ -32,12 +32,13 @@ link_directories( ${LINK_DIRECTORIES})
set(SRCS
CbmMuchFindHitsGem.cxx
-CbmMuchHitFinderQa.cxx
CbmMuchHitProducerIdeal.cxx
CbmMuchFindTracks.cxx
CbmMuchMatchTracks.cxx
CbmMuchTrackFinderIdeal.cxx
+
+qa/CbmMuchHitFinderQa.cxx
)
set(LINKDEF CbmMuchRecoLinkDef.h)
diff --git a/reco/detectors/much/CbmMuchHitFinderQa.cxx b/reco/detectors/much/CbmMuchHitFinderQa.cxx
deleted file mode 100644
index 5b5ff70d5f..0000000000
--- a/reco/detectors/much/CbmMuchHitFinderQa.cxx
+++ /dev/null
@@ -1,1642 +0,0 @@
-// -------------------------------------------------------------------------
-// ----- CbmMuchHitProducerQa source file -----
-// ----- Modified since 21/06/2019 by Ekata Nandy (ekata@vecc.gov.in) -Inclusion of RPC detector type
-// ----- Modified 02/18 by Vikas Singhal -----
-// ----- Created 16/11/07 by E. Kryshen -----
-// -------------------------------------------------------------------------
-// AUG 18 RELEASE VERSION
-#include "CbmMuchHitFinderQa.h"
-#include "CbmMuchCluster.h"
-#include "CbmMuchDigi.h"
-#include "CbmMuchPixelHit.h"
-#include "CbmMuchPoint.h"
-#include "FairRootManager.h"
-
-#include "CbmDigiManager.h"
-#include "CbmMuchAddress.h"
-#include "CbmMuchDigi.h"
-#include "CbmMuchModuleGem.h"
-#include "CbmMuchPad.h"
-#include "CbmMuchSector.h"
-#include "CbmMuchStation.h"
-
-#include "CbmMCTrack.h"
-#include "CbmMatch.h"
-#include "FairLogger.h"
-#include "TDatabasePDG.h"
-#include "TParticlePDG.h"
-
-#include "CbmMuchPointInfo.h"
-#include "CbmMuchRecoDefs.h"
-#include "TCanvas.h"
-#include "TF1.h"
-#include "TFile.h"
-
-#include "CbmMuchGeoScheme.h"
-#include "TArrayI.h"
-#include "TObjArray.h"
-#include "TStyle.h"
-
-
-#include "CbmGeoMuchPar.h"
-#include "FairRuntimeDb.h"
-#include "TGraph.h"
-
-#include <algorithm>
-#include <cassert>
-#include <map>
-#include <vector>
-
-using std::cout;
-using std::endl;
-using std::map;
-using std::vector;
-// -------------------------------------------------------------------------
-CbmMuchHitFinderQa::CbmMuchHitFinderQa(const char* name, Int_t verbose)
- : FairTask(name, verbose)
- , fGeoScheme(CbmMuchGeoScheme::Instance())
- , fGeoFileName()
- , fFileName()
- , fSignalPoints(0)
- , fSignalHits(0)
- , fVerbose(verbose)
- , fEvent(0)
- , fFlag(0)
- , fPoints(NULL)
- ,
- //fDigis(NULL),
- //fDigiMatches(NULL),
- fClusters(NULL)
- , fHits(NULL)
- , fMCTracks(NULL)
- , fPointInfos(new TClonesArray("CbmMuchPointInfo", 10))
- , fNstations(0)
- , fChargeHistos(NULL)
- , fhChargeEnergyLog(NULL)
- , fhChargeEnergyLogPi(NULL)
- , fhChargeEnergyLogPr(NULL)
- , fhChargeEnergyLogEl(NULL)
- , fhChargeTrackLength(NULL)
- , fhChargeTrackLengthPi(NULL)
- , fhChargeTrackLengthPr(NULL)
- , fhChargeTrackLengthEl(NULL)
- , fhChargeLog(NULL)
- , fhChargePr_1GeV_3mm(NULL)
- , fhNpadsVsS(NULL)
- , fhCharge(NULL)
- , fhOccupancyR(NULL)
- , fhPadsTotalR(NULL)
- , fhPadsFiredR(NULL)
- , fhPullXpads1(NULL)
- , fhPullYpads1(NULL)
- , fhPullXpads2(NULL)
- , fhPullYpads2(NULL)
- , fhPullXpads3(NULL)
- , fhPullYpads3(NULL)
- , fnPadSizesX(0)
- , fnPadSizesY(0)
- , fNTimingPulls(8)
- , fhPullX(NULL)
- , fhPullY(NULL)
- , fhPullT(NULL)
- , fhResidualX(NULL)
- , fhResidualY(NULL)
- , fhResidualT(NULL)
- , fhPointsInCluster(NULL)
- , fhDigisInCluster(NULL)
- , fhHitsInCluster(NULL)
- , fNall(NULL)
- , fNpr(NULL)
- , fNpi(NULL)
- , fNel(NULL)
- , fNmu(NULL)
- , fNka(NULL)
- , fNprimary(NULL)
- , fNsecondary(NULL)
- , fPointsTotal(0)
- , fPointsUnderCounted(0)
- , fPointsOverCounted(0)
- , fOccupancyQaOn(kTRUE)
- , fPullsQaOn(kTRUE)
- , fDigitizerQaOn(kTRUE)
- , fStatisticsQaOn(kTRUE)
- , fClusterDeconvQaOn(kTRUE)
- , fPrintToFileOn(kTRUE)
- , fPadMinLx(0.)
- , fPadMinLy(0.)
- , fPadMaxLx(0.)
- , fPadMaxLy(0.)
- , pointsFile(NULL)
- , padsFile(NULL)
-
-{}
-// -------------------------------------------------------------------------
-
-
-// -------------------------------------------------------------------------
-CbmMuchHitFinderQa::~CbmMuchHitFinderQa() {}
-// -------------------------------------------------------------------------
-
-
-// -------------------------------------------------------------------------
-InitStatus CbmMuchHitFinderQa::Init() {
- FairRootManager* fManager = FairRootManager::Instance();
- fMCTracks = (TClonesArray*) fManager->GetObject("MCTrack");
- fPoints = (TClonesArray*) fManager->GetObject("MuchPoint");
- fHits = (TClonesArray*) fManager->GetObject("MuchPixelHit");
- //fDigis = (TClonesArray*) fManager->GetObject("MuchDigi");
- //fDigiMatches = (TClonesArray*) fManager->GetObject("MuchDigiMatch"); ///added ekata
- fClusters = (TClonesArray*) fManager->GetObject("MuchCluster");
- // Reading Much Digis from CbmMuchDigiManager which are stored as vector
- fDigiManager = CbmDigiManager::Instance();
- fDigiManager->Init();
-
- // printf(" %i",fMCTracks);
- // printf(" %i",fPoints);
- // printf(" %i",fHits);
- // printf(" %i",fDigis);
- // printf(" %i",fDigiMatches);
- // printf(" %i",fClusters);
- // printf("\n");
-
- TFile* f = new TFile(fGeoFileName, "R");
- TObjArray* stations = (TObjArray*) f->Get("stations");
- fGeoScheme->Init(stations, fFlag);
- fNstations = fGeoScheme->GetNStations();
- printf("Init: fNstations = %i\n", fNstations);
-
- 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];
-
- // Reset counters
- 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;
-
-#define BINNING_CHARGE 100, 0, 3.0
-#define BINNING_LENGTH 100, 0, 2.5
-#define BINNING_CHARGE_LOG 100, 4, 8
-#define BINNING_ENERGY_LOG 100, -0.5, 4.5
-#define BINNING_ENERGY_LOG_EL 100, -0.5, 4.5
-
- // fhCharge = new TH1D("hCharge",
- // "Charge distribution from tracks",
- // BINNING_CHARGE);
-
- fhChargeLog = new TH1D(
- "hChargeLog", "Charge (log.) distribution from tracks", BINNING_CHARGE_LOG);
-
- fhChargeEnergyLog = new TH2D("hChargeEnergyLog",
- "Charge vs energy (log.) for all tracks",
- BINNING_ENERGY_LOG,
- BINNING_CHARGE);
-
- fhChargeEnergyLogPi = new TH2D("hChargeEnergyLogPi",
- "Charge vs energy (log.) for pions",
- BINNING_ENERGY_LOG,
- BINNING_CHARGE);
-
- fhChargeEnergyLogPr = new TH2D("hChargeEnergyLogPr",
- "Charge vs energy (log.) for protons",
- BINNING_ENERGY_LOG,
- BINNING_CHARGE);
-
- fhChargeEnergyLogEl = new TH2D("hChargeEnergyLogEl",
- "Charge vs energy (log.) for electrons",
- BINNING_ENERGY_LOG_EL,
- BINNING_CHARGE);
-
- fhChargeTrackLength = new TH2D("hChargeTrackLength",
- "Charge vs length for all tracks",
- BINNING_LENGTH,
- BINNING_CHARGE);
-
- fhChargeTrackLengthPi = new TH2D("hChargeTrackLengthPi",
- "Charge vs length for pions",
- BINNING_LENGTH,
- BINNING_CHARGE);
- fhChargeTrackLengthPr = new TH2D("hChargeTrackLengthPr",
- "Charge vs length for proton",
- BINNING_LENGTH,
- BINNING_CHARGE);
-
- fhChargeTrackLengthEl = new TH2D("hChargeTrackLengthEl",
- "Charge vs length for electrons",
- BINNING_LENGTH,
- BINNING_CHARGE);
-
- fhChargePr_1GeV_3mm =
- new TH1D("hChargePr_1GeV_3mm", "Charge for 1 GeV protons", BINNING_CHARGE);
-
- // fhPointsInCluster = new TH1I("hPointsInCluster",
- // "Number of tracks in cluster",
- // 10,0,10);
-
- // fhDigisInCluster = new TH1I("hDigisInCluster",
- // "Number of digis in cluster",
- // 10,0,10);
-
- // fhHitsInCluster = new TH1I("hHitsInCluster",
- // "Number of hits in cluster",
- // 10,0,10);
-
-
- // fhCharge ->GetXaxis()->SetTitle("Charge [10^{6} electrons]");
- fhChargeLog->GetXaxis()->SetTitle("Lg (Charge) [Number of electrons]");
- fhChargePr_1GeV_3mm->GetXaxis()->SetTitle("Charge [10^{6} electrons]");
-
-
- fChargeHistos = new TObjArray();
- fChargeHistos->Add(fhChargeEnergyLog);
- fChargeHistos->Add(fhChargeEnergyLogEl);
- fChargeHistos->Add(fhChargeEnergyLogPi);
- fChargeHistos->Add(fhChargeEnergyLogPr);
- fChargeHistos->Add(fhChargeTrackLength);
- fChargeHistos->Add(fhChargeTrackLengthEl);
- fChargeHistos->Add(fhChargeTrackLengthPi);
- fChargeHistos->Add(fhChargeTrackLengthPr);
- // fChargeHistos->Add(fhCharge);
- fChargeHistos->Add(fhChargeLog);
-
- for (Int_t i = 0; i < 8; i++) {
- TH2D* histo = (TH2D*) fChargeHistos->At(i);
- histo->SetStats(0);
- histo->GetYaxis()->SetDecimals(kTRUE);
- histo->GetYaxis()->SetTitleOffset(1.4);
- histo->GetYaxis()->CenterTitle();
- histo->GetYaxis()->SetTitle("Charge [10^{6} electrons]");
- if (i < 4)
- histo->GetXaxis()->SetTitle("Lg E_{kin} [MeV]");
- else
- histo->GetXaxis()->SetTitle("Track length [cm]");
- }
-
- fhPointsInCluster = new TH1I*[fNstations];
- fhDigisInCluster = new TH1I*[fNstations];
- fhHitsInCluster = new TH1I*[fNstations];
-
- fhCharge = new TH1D*[fNstations];
- fhOccupancyR = new TH1D*[fNstations];
- fhPadsTotalR = new TH1D*[fNstations];
- fhPadsFiredR = new TH1D*[fNstations];
-
- for (Int_t i = 0; i < fNstations; i++) {
- CbmMuchStation* station = fGeoScheme->GetStation(i);
- Double_t rMax = station->GetRmax();
- Double_t rMin = station->GetRmin();
- fhPadsTotalR[i] =
- new TH1D(Form("hPadsTotal%i", i + 1),
- Form("Number of pads vs radius: station %i;Radius [cm]", i + 1),
- 100,
- 0.6 * rMin,
- 1.2 * rMax);
- fhPadsFiredR[i] = new TH1D(
- 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);
- fhOccupancyR[i] = new TH1D(
- Form("hOccupancy%i", i + 1),
- Form("Occupancy vs radius: station %i;Radius [cm];Occupancy, %%", i + 1),
- 100,
- 0.6 * rMin,
- 1.2 * rMax);
-
- fhCharge[i] =
- new TH1D(Form("hCharge%i", i + 1),
- Form("Charge : station %i; Charge(1e4); Count", i + 1),
- 200,
- 0,
- 500);
- fhPointsInCluster[i] =
- new TH1I(Form("hPointsInCluster%i", i + 1),
- Form("Points in Cluster : Station %i ", i + 1),
- 10,
- 0,
- 10);
- fhDigisInCluster[i] =
- new TH1I(Form("hDigisInCluster%i", i + 1),
- Form("Digis in Cluster : Station %i ", i + 1),
- 10,
- 0,
- 10);
- fhHitsInCluster[i] = new TH1I(Form("hHitsInCluster%i", i + 1),
- Form("Hits in Cluster : Station %i ", i + 1),
- 10,
- 0,
- 10);
- }
-
- 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);
- fhPadsTotalR[stationId]->Fill(r0);
- }
- }
- }
- /*
- 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();
-
- Int_t nTotSectors = 0;
- Int_t nTotChannels = 0;
- printf("=========================================================================================\n");
- printf(" Station Nr.\t| Sectors\t| Channels\t| Pad min size\t\t| Pad max length\t \n");
- printf("-----------------------------------------------------------------------------------------\n");
- */
- Int_t nTotChannels = 0;
- for (Int_t iStation = 0; iStation < fNstations; 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();
- */
- Int_t nChannels = GetNChannels(iStation);
- /*
- Int_t nSectors = GetNSectors(iStation);
- if (fPadMinLx>padMinLx) fPadMinLx = padMinLx;
- if (fPadMinLy>padMinLy) fPadMinLy = padMinLy;
- if (fPadMaxLx<padMaxLx) fPadMaxLx = padMaxLx;
- if (fPadMaxLy<padMaxLy) fPadMaxLy = padMaxLy;
- 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);
- */
- // nTotSectors += nSectors;
- printf("%i\t\t| %i\t\t\n", iStation + 1, nChannels);
- nTotChannels += nChannels;
- }
- printf("---------------------------------------------------------------------"
- "--------------------\n");
- printf(" Total:\t\t| %i\t\t\n", nTotChannels);
- printf("====================================================================="
- "====================\n");
- /*
- fnPadSizesX = TMath::CeilNint(TMath::Log2(fPadMaxLx/fPadMinLx)+1);
- fnPadSizesY = TMath::CeilNint(TMath::Log2(fPadMaxLy/fPadMinLy)+1);
- Info("Init","nPadSizesX=%i",fnPadSizesX);
- Info("Init","nPadSizesY=%i",fnPadSizesY);
- fhPullXpads1 = new TH1D*[fnPadSizesX];
- fhPullYpads1 = new TH1D*[fnPadSizesY];
- fhPullXpads2 = new TH1D*[fnPadSizesX];
- fhPullYpads2 = new TH1D*[fnPadSizesY];
- fhPullXpads3 = new TH1D*[fnPadSizesX];
- fhPullYpads3 = new TH1D*[fnPadSizesY];
- fhPullT = new TH1D*[fNTimingPulls];
-
- for (Int_t i=0;i<fnPadSizesX;i++){
- fhPullXpads1[i] = new TH1D(Form("hPullXpads1%i",i),Form("Pull distribution X. Npads = 1 Size =%i; (x_{RC} - x_{MC}) / dx_{RC}",i),100,-5,5);
- fhPullXpads2[i] = new TH1D(Form("hPullXpads2%i",i),Form("Pull distribution X. Npads = 2 Size =%i; (x_{RC} - x_{MC}) / dx_{RC}",i),100,-5,5);
- fhPullXpads3[i] = new TH1D(Form("hPullXpads3%i",i),Form("Pull distribution X. Npads = 3 Size =%i; (x_{RC} - x_{MC}) / dx_{RC}",i),100,-5,5);
- }
-
- for (Int_t i=0;i<fnPadSizesY;i++){
- fhPullYpads1[i] = new TH1D(Form("hPullYpads1%i",i),Form("Pull distribution Y. Npads = 1 Size =%i; (y_{RC} - y_{MC}) / dy_{RC}",i),100,-5,5);
- fhPullYpads2[i] = new TH1D(Form("hPullYpads2%i",i),Form("Pull distribution Y. Npads = 2 Size =%i; (y_{RC} - y_{MC}) / dy_{RC}",i),100,-5,5);
- fhPullYpads3[i] = new TH1D(Form("hPullYpads3%i",i),Form("Pull distribution Y. Npads = 3 Size =%i; (y_{RC} - y_{MC}) / dy_{RC}",i),100,-5,5);
- }
-
- fhPullT[0] = new TH1D("hPullT","Pull distribution T, all pads; (t_{RC} - t_{MC}) / dt_{RC}",100,-5,5);
- for (Int_t i=1;i<fNTimingPulls;i++){
- fhPullT[i] = new TH1D(Form("hPullT%i",i),Form("Pull distribution T. Npads = %i; (t_{RC} - t_{MC}) / dt_{RC}",i),100,-5,5);
- }
- */
- //Pull Distribution
- fhPullX = new TH1D(
- "hPullX", "Pull distribution X;(x_{RC} - x_{MC}) / dx_{RC}", 500, -5, 5);
- fhPullY = new TH1D(
- "hPullY", "Pull distribution Y;(y_{RC} - y_{MC}) / dy_{RC}", 500, -5, 5);
- fhPullT = new TH1D(
- "hPullT", "Pull distribution T;(t_{RC} - t_{MC}) / dt_{RC}", 120, -10, 50);
-
- //Residual Distribution
- fhResidualX = new TH1D(
- "hResidualX", "Residual distribution X;(x_{RC} - x_{MC})(cm)", 500, -5, 5);
- fhResidualY = new TH1D(
- "hResidualY", "Residual distribution Y;(y_{RC} - y_{MC})(cm)", 500, -5, 5);
- fhResidualT = new TH1D("hResidualT",
- "Residual distribution T;(t_{RC} - t_{MC})(ns)",
- 140,
- -20,
- 50);
-
- fhNpadsVsS = new TH2D("hNpadsVsS",
- "Number of fired pads vs pad area:area:n pads",
- 10,
- -5,
- 0,
- 10,
- 0.5,
- 10.5);
- pointsFile = fopen("points.txt", "w");
- padsFile = fopen("pads.txt", "w");
- return kSUCCESS;
-}
-// -------------------------------------------------------------------------
-
-
-// -------------------------------------------------------------------------
-void CbmMuchHitFinderQa::SetParContainers() {
- // Get run and runtime database
- // FairRuntimeDb* db = FairRuntimeDb::instance();
- // if ( ! db ) Fatal("SetParContainers", "No runtime database");
- // Get MUCH geometry parameter container
- // fGeoPar = (CbmGeoMuchPar*) db->getContainer("CbmGeoMuchPar");
-}
-// -------------------------------------------------------------------------
-
-
-// -------------------------------------------------------------------------x
-void CbmMuchHitFinderQa::Exec(Option_t*) {
- fEvent++;
- LOG(info) << "Event: " << fEvent;
- fprintf(pointsFile, "Event %i\n", fEvent);
- fprintf(padsFile, "Event %i\n", fEvent);
-
- if (fPullsQaOn) PullsQa();
- if (fOccupancyQaOn) OccupancyQa();
- if (fDigitizerQaOn) DigitizerQa();
- if (fStatisticsQaOn) StatisticsQa();
- if (fClusterDeconvQaOn) ClusterDeconvQa();
-
- return;
- 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());
- fprintf(pointsFile,
- "%7.3f %7.3f %7.3f %7.3f\n",
- point->GetXIn(),
- point->GetYIn(),
- point->GetXOut(),
- point->GetYOut());
- }
-
- // old code
- /* for (Int_t i=0;i<fDigis->GetEntriesFast();i++){
- CbmMuchDigi* digi = (CbmMuchDigi*) fDigis->At(i); */
- 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);
- fprintf(padsFile, "%5.1f %5.1f %3i\n", x0, y0, charge);
- }
-}
-// -------------------------------------------------------------------------
-
-
-// -------------------------------------------------------------------------
-void CbmMuchHitFinderQa::FinishTask() {
- printf("FinishTask\n");
- fclose(pointsFile);
- fclose(padsFile);
-
- gStyle->SetPaperSize(20, 20);
- // Int_t NuOfColoums = 1;
- // Int_t NuOfRows = 1;
-
- //Setting Canvas according to the Number of stations
- /*
- if(fNstations>4){ NuOfColoums = 3; NuOfRows = 2;}
- else{
- if(fNstations<=4 && fNstations >=3){ NuOfColoums = 2; NuOfRows = 2;}
- else { NuOfColoums = fNstations; NuOfRows = 1;}
- }
-*/
-
- Double_t errors[100];
- for (Int_t i = 0; i < 100; i++)
- errors[i] = 0;
-
- for (Int_t i = 0; i < fNstations; i++) {
- fhPadsTotalR[i]->Sumw2();
- fhPadsTotalR[i]->SetError(errors);
- fhPadsFiredR[i]->Sumw2();
- fhPadsFiredR[i]->Scale(1. / fEvent);
- fhOccupancyR[i]->Divide(fhPadsFiredR[i], fhPadsTotalR[i]);
- fhOccupancyR[i]->Scale(100);
- fhCharge[i]->Scale(1. / fEvent);
- fhPointsInCluster[i]->Scale(1. / fEvent);
- fhDigisInCluster[i]->Scale(1. / fEvent);
- fhHitsInCluster[i]->Scale(1. / fEvent);
- }
-
- if (fPullsQaOn && fVerbose > 1) {
- printf("===================================\n");
- printf("PullsQa:\n");
-
- // TCanvas* cTiming = new TCanvas("cTiming","Timing pulls",250*4,250*2);
- // cTiming->Divide(4,2);
- // for (Int_t i=0;i<fNTimingPulls;i++){
- // cTiming->cd(i+1);
- // fhPullT[i]->Fit("gaus");
- // fhPullT[i]->Draw("e");
- // }
- //
- /* TCanvas* c4 = new TCanvas("c4","Pulls",800,400);
- c4->Divide(3,1);
- c4->cd(1);*/
- fhPullX->Sumw2();
- fhPullX->Fit("gaus", "Q");
- fhPullX->GetFunction("gaus")->SetLineWidth(3);
- fhPullX->GetFunction("gaus")->SetLineColor(kRed);
- // fhPullX->Draw();
- fhPullX->Write();
- if (fPrintToFileOn) gPad->Print(".gif");
- if (fPrintToFileOn) gPad->Print(".eps");
- // c4->cd(2);
- fhPullY->Sumw2();
- fhPullY->Fit("gaus", "Q");
- fhPullY->GetFunction("gaus")->SetLineWidth(3);
- fhPullY->GetFunction("gaus")->SetLineColor(kRed);
- // fhPullY->Draw();
- fhPullY->Write();
- if (fPrintToFileOn) gPad->Print(".gif");
- if (fPrintToFileOn) gPad->Print(".eps");
- // c4->cd(3);
- fhPullT->Sumw2();
- fhPullT->Fit("gaus", "Q");
- fhPullT->GetFunction("gaus")->SetLineWidth(3);
- fhPullT->GetFunction("gaus")->SetLineColor(kRed);
- // fhPullT->Draw();
- fhPullT->Write();
- if (fPrintToFileOn) gPad->Print(".gif");
- if (fPrintToFileOn) gPad->Print(".eps");
- // c4->cd();
-
- // TCanvas* cR5 = new TCanvas("cR5","Residuals",800,400);
- // cR5->Divide(3,1);
- // cR5->cd(1);
- fhResidualX->Sumw2();
- fhResidualX->Fit("gaus", "Q");
- fhResidualX->GetFunction("gaus")->SetLineWidth(3);
- fhResidualX->GetFunction("gaus")->SetLineColor(kRed);
- // fhResidualX->Draw();
- fhResidualX->Write();
- // if (fPrintToFileOn) gPad->Print(".gif");
- // if (fPrintToFileOn) gPad->Print(".eps");
- // cR5->cd(2);
- fhResidualY->Sumw2();
- fhResidualY->Fit("gaus", "Q");
- fhResidualY->GetFunction("gaus")->SetLineWidth(3);
- fhResidualY->GetFunction("gaus")->SetLineColor(kRed);
- // fhResidualY->Draw();
- fhResidualY->Write();
- // if (fPrintToFileOn) gPad->Print(".gif");
- // if (fPrintToFileOn) gPad->Print(".eps");
- // cR5->cd(3);
- fhResidualT->Sumw2();
- fhResidualT->Fit("gaus", "Q");
- fhResidualT->GetFunction("gaus")->SetLineWidth(3);
- fhResidualT->GetFunction("gaus")->SetLineColor(kRed);
- /// fhResidualT->Draw();
- fhResidualT->Write();
- // if (fPrintToFileOn) gPad->Print(".gif");
- // if (fPrintToFileOn) gPad->Print(".eps");
- // cR5->cd();
-
- // TCanvas* c_alone = new TCanvas("c_alone","Pulls",400,400);
- // new TCanvas("c_alone","Pulls",400,400);
- // fhPullX->Draw();
- // gPad->Print(".gif");
- //
- // TCanvas* c4x = new TCanvas("c4x","X-pulls vs pad size and cluster size",fnPadSizesX*300,3*300);
- // c4x->Divide(fnPadSizesX,3);
- // for (Int_t i=0;i<fnPadSizesX;i++){
- // c4x->cd(i+1);
- // fhPullXpads1[i]->Sumw2();
- // fhPullXpads1[i]->Fit("gaus","Q");
- // fhPullXpads1[i]->GetFunction("gaus")->SetLineWidth(2);
- // fhPullXpads1[i]->GetFunction("gaus")->SetLineColor(kBlue);
- // fhPullXpads1[i]->Draw();
- // if (fPrintToFileOn) gPad->Print(".gif");
- // if (fPrintToFileOn) gPad->Print(".eps");
- //
- // c4x->cd(fnPadSizesX+i+1);
- // fhPullXpads2[i]->Sumw2();
- // fhPullXpads2[i]->Fit("gaus","Q");
- // fhPullXpads2[i]->GetFunction("gaus")->SetLineWidth(2);
- // fhPullXpads2[i]->GetFunction("gaus")->SetLineColor(kBlue);
- // fhPullXpads2[i]->Draw();
- // if (fPrintToFileOn) gPad->Print(".gif");
- // if (fPrintToFileOn) gPad->Print(".eps");
- //
- // c4x->cd(2*fnPadSizesX+i+1);
- // fhPullXpads3[i]->Sumw2();
- // fhPullXpads3[i]->Fit("gaus","Q");
- // fhPullXpads3[i]->GetFunction("gaus")->SetLineWidth(2);
- // fhPullXpads3[i]->GetFunction("gaus")->SetLineColor(kBlue);
- // fhPullXpads3[i]->Draw();
- // if (fPrintToFileOn) gPad->Print(".gif");
- // if (fPrintToFileOn) gPad->Print(".eps");
- // }
-
- // TCanvas* c4y = new TCanvas("c4y","Y-pulls vs pad size and cluster size",fnPadSizesY*300,3*300);
- // c4y->Divide(fnPadSizesY,3);
- // for (Int_t i=0;i<fnPadSizesY;i++){
- // c4y->cd(i+1);
- // fhPullYpads1[i]->Sumw2();
- // fhPullYpads1[i]->Fit("gaus","Q");
- // fhPullYpads1[i]->GetFunction("gaus")->SetLineWidth(2);
- // fhPullYpads1[i]->GetFunction("gaus")->SetLineColor(kBlue);
- // fhPullYpads1[i]->Draw();
- // if (fPrintToFileOn) gPad->Print(".gif");
- // if (fPrintToFileOn) gPad->Print(".eps");
- //
- // c4y->cd(fnPadSizesY+i+1);
- // fhPullYpads2[i]->Sumw2();
- // fhPullYpads2[i]->Fit("gaus","Q");
- // fhPullYpads2[i]->GetFunction("gaus")->SetLineWidth(2);
- // fhPullYpads2[i]->GetFunction("gaus")->SetLineColor(kBlue);
- // fhPullYpads2[i]->Draw();
- // if (fPrintToFileOn) gPad->Print(".gif");
- // if (fPrintToFileOn) gPad->Print(".eps");
- //
- // c4y->cd(2*fnPadSizesY+i+1);
- // fhPullYpads3[i]->Sumw2();
- // fhPullYpads3[i]->Fit("gaus","Q");
- // fhPullYpads3[i]->GetFunction("gaus")->SetLineWidth(2);
- // fhPullYpads3[i]->GetFunction("gaus")->SetLineColor(kBlue);
- // fhPullYpads3[i]->Draw();
- // if (fPrintToFileOn) gPad->Print(".gif");
- // if (fPrintToFileOn) gPad->Print(".eps");
- // }
- }
-
-
- if (fOccupancyQaOn && fVerbose > 1) {
- printf("===================================\n");
- printf("OccupancyQa:\n");
- // TCanvas* c1 = new TCanvas("c1","All pad distributions",1200,800);
- // c1->Divide(NuOfColoums,NuOfRows);
- for (Int_t i = 0; i < fNstations; i++) {
- // c1->cd(i+1);
- // fhPadsTotalR[i]->SetStats(0);
- // fhPadsTotalR[i]->Draw("hist");
- // if (fPrintToFileOn) gPad->Print(".gif");
- // if (fPrintToFileOn) gPad->Print(".eps");
- }
- // c1->cd();
-
- // TCanvas* c2 = new TCanvas("c2","Fired pad distributions",1200,800);
- // c2->Divide(NuOfColoums,NuOfRows);
- for (Int_t i = 0; i < fNstations; i++) {
- // c2->cd(i+1);
- // fhPadsFiredR[i]->SetStats(0);
- // fhPadsFiredR[i]->Draw();
- fhPadsFiredR[i]->Write();
- // if (fPrintToFileOn) gPad->Print(".gif");
- // if (fPrintToFileOn) gPad->Print(".eps");
- }
- // c2->cd();
-
- // TCanvas* c3 = new TCanvas("c3","Occupancy plots",2400,1600);
- // c3->Divide(NuOfColoums,NuOfRows);
- for (Int_t i = 0; i < fNstations; i++) {
- // c3->cd(i+1);
- // fhOccupancyR[i]->SetStats(0);
- // fhOccupancyR[i]->Draw();
- fhOccupancyR[i]->Write();
- // if (fPrintToFileOn) gPad->Print(".gif");
- // if (fPrintToFileOn) gPad->Print(".eps");
- }
- // c3->cd();
-
- // TCanvas *c4 = new TCanvas("c4","Charge plot",2400,1600); ////added for station charge
- // c4->Divide(3,2);
- for (Int_t i = 0; i < fNstations; i++) {
- // c4->cd(i+1);
- // fhCharge[i]->Draw();
- fhCharge[i]->Write();
- // fhPointsInCluster[i]->Draw();
- fhPointsInCluster[i]->Write();
- fhDigisInCluster[i]->Write();
- fhHitsInCluster[i]->Write();
- }
-
-
- /*
- // TCanvas* oc_1 = new TCanvas("oc_1","",1200,1000);
- new TCanvas("oc_1","",1200,1000);
- fhOccupancyR[0]->GetYaxis()->SetTitleOffset(1.6);
- fhOccupancyR[0]->Draw();
- if (fPrintToFileOn) gPad->Print(".gif");
- // TCanvas* oc_2 = new TCanvas("oc_2","",1200,1000);
- new TCanvas("oc_2","",1200,1000);
- fhOccupancyR[1]->GetYaxis()->SetTitleOffset(1.6);
- fhOccupancyR[1]->Draw();
- if (fPrintToFileOn) gPad->Print(".gif");
- // TCanvas* oc_3 = new TCanvas("oc_3","",1200,1000);
- new TCanvas("oc_3","",1200,1000);
- fhOccupancyR[2]->GetYaxis()->SetTitleOffset(1.6);
- fhOccupancyR[2]->Draw();
- if (fPrintToFileOn) gPad->Print(".gif");*/
- }
-
- if (fDigitizerQaOn && fVerbose > 1) {
- printf("===================================\n");
- printf("DigitizerQa:\n");
-
- TF1* fit_el = new TF1("fit_el", LandauMPV, -0.5, 4.5, 1);
- fit_el->SetParameter(0, 0.511);
- fit_el->SetLineWidth(2);
- fit_el->SetLineColor(kBlack);
-
- TF1* fit_pi = new TF1("fit_pi", LandauMPV, -0.5, 4.5, 1);
- fit_pi->SetParameter(0, 139.57);
- fit_pi->SetLineWidth(2);
- fit_pi->SetLineColor(kBlack);
-
- TF1* fit_pr = new TF1("fit_pr", LandauMPV, -0.5, 4.5, 1);
- fit_pr->SetParameter(0, 938.27);
- fit_pr->SetLineWidth(2);
- fit_pr->SetLineColor(kBlack);
-
- TH1D* hLength = fhChargeTrackLength->ProjectionX();
- TH1D* hLengthEl = fhChargeTrackLengthEl->ProjectionX();
- ;
- TH1D* hLengthPi = fhChargeTrackLengthPi->ProjectionX();
- ;
- TH1D* hLengthPr = fhChargeTrackLengthPr->ProjectionX();
- ;
- hLength->SetTitle("Track length for all tracks");
- hLengthEl->SetTitle("Track length for electrons");
- hLengthPi->SetTitle("Track length for pions");
- hLengthPr->SetTitle("Track length for protons");
-
- fChargeHistos->Add(hLength);
- fChargeHistos->Add(hLengthEl);
- fChargeHistos->Add(hLengthPi);
- fChargeHistos->Add(hLengthPr);
-
- // TCanvas* c5 = new TCanvas("c5","Charge vs energy and length",2100,1300);
- // c5->Divide(4,3);
- for (Int_t i = 0; i < 8; i++) {
- // c5->cd(i+1);
- // gPad->Range(0,0,200,200);
- // gPad->SetBottomMargin(0.11);
- // gPad->SetRightMargin(0.14);
- //gPad->SetLeftMargin(0.12);
- //gPad->SetLogz();
- // TH2D* histo = (TH2D*) fChargeHistos->At(i);
- // histo->Draw("colz");
- // if (i==1) fit_el->Draw("same");
- // if (i==2) fit_pi->Draw("same");
- // if (i==3) fit_pr->Draw("same");
- // if (fPrintToFileOn) gPad->Print(".gif");
- // if (fPrintToFileOn) gPad->Print(".eps");
- }
-
- // TCanvas* ch_1 = new TCanvas("ch_1","",1200,1000);
- // new TCanvas("ch_1","",1200,1000);
- // gPad->SetBottomMargin(0.11);
- // gPad->SetRightMargin(0.14);
- // gPad->SetLeftMargin(0.12);
- // gPad->SetLogz();
- // TH2D* histo1 = (TH2D*) fChargeHistos->At(1);
- // histo1->Draw("colz");
- // fit_el->Draw("same");
- // if (fPrintToFileOn) gPad->Print(".gif");
-
- // TCanvas* ch_2 = new TCanvas("ch_2","",1200,1000);
- /* new TCanvas("ch_2","",1200,1000);
- gPad->SetBottomMargin(0.11);
- gPad->SetRightMargin(0.14);
- gPad->SetLeftMargin(0.12);
- gPad->SetLogz();*/
- // TH2D* histo2 = (TH2D*) fChargeHistos->At(2);
- /* histo2->Draw("colz");
- fit_pi->Draw("same");
- if (fPrintToFileOn) gPad->Print(".gif");*/
-
- // TCanvas* ch_3 = new TCanvas("ch_3","",1200,1000);
- /* new TCanvas("ch_3","",1200,1000);
- gPad->SetBottomMargin(0.11);
- gPad->SetRightMargin(0.14);
- gPad->SetLeftMargin(0.12);
- gPad->SetLogz(); */
- // TH2D* histo3 = (TH2D*) fChargeHistos->At(3);
- /* histo3->Draw("colz");
- fit_pr->Draw("same");
- if (fPrintToFileOn) gPad->Print(".gif"); */
-
-
- for (Int_t i = 10; i < 14; i++) {
- // c5->cd(i-1);
- // gPad->SetLogy();
- // gStyle->SetOptStat(1110);
- // TH1D* histo = (TH1D*) fChargeHistos->At(i);
- // histo->Draw();
- //if (fPrintToFileOn) gPad->Print(".gif");
- //if (fPrintToFileOn) gPad->Print(".eps");
- }
-
- // below for 31st CBM Collaboration Meeting
- // TCanvas* c6 = new TCanvas("c6","Charge distribution",1200,400);
- // c6->Divide(2,1);
- // c6->cd(1);
- // fhCharge->Draw();
- // fhCharge->Write();
- // if (fPrintToFileOn) gPad->Print(".gif");
- // if (fPrintToFileOn) gPad->Print(".eps");
-
- //c6->cd(2);
- // fhChargeLog->Draw();
- fhChargeLog->Write();
- // if (fPrintToFileOn) gPad->Print(".gif");
- // if (fPrintToFileOn) gPad->Print(".eps");
-
-
- /* Commented below for 31st CBM Collaboration Meeting
- TCanvas* c6 = new TCanvas("c6","Charge distribution",1200,400);
- c6->Divide(3,1);
- c6->cd(1);
- fhCharge->Draw();
- if (fPrintToFileOn) gPad->Print(".gif");
- if (fPrintToFileOn) gPad->Print(".eps");
-
- c6->cd(2);
- fhChargeLog->Draw();
- if (fPrintToFileOn) gPad->Print(".gif");
- if (fPrintToFileOn) gPad->Print(".eps");
-
- c6->cd(3);
- fhChargePr_1GeV_3mm->Draw();
- if (fPrintToFileOn) gPad->Print(".gif");
- if (fPrintToFileOn) gPad->Print(".eps");
- */
-
- // TCanvas* c8 = new TCanvas("c8","Square vs nPads",800,400);
- // c8->Divide(2,1);
- //c8->cd(1);
- // fhNpadsVsS->Draw("colz");
- fhNpadsVsS->Write();
- 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]);
- }
- // c8->cd(2);
- // TGraph* gNvsS = new TGraph(11,s,nMean);
- //gNvsS->DrawClone();
-
-
- 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");
- }
-
-
- if (fStatisticsQaOn) {
- printf("===================================\n");
- printf("StatisticsQa:\n");
- // TCanvas* c7 = new TCanvas("c7","Cluster statistics",1200,400);
- // c7->Divide(3,1);
- // c7->cd(1);
- // gStyle->SetOptStat(1110);
- // gPad->SetLogy();
- // fhPointsInCluster->Draw();
- // fhPointsInCluster->Write();
- // if (fPrintToFileOn) gPad->Print(".gif");
- // if (fPrintToFileOn) gPad->Print(".eps");
- // c7->cd(2);
- // gStyle->SetOptStat(1110);
- // gPad->SetLogy();
- // fhDigisInCluster->Draw();
- // fhDigisInCluster->Write();
- // if (fPrintToFileOn) gPad->Print(".gif");
- // if (fPrintToFileOn) gPad->Print(".eps");
- // c7->cd(3);
- // gStyle->SetOptStat(1110);
- // gPad->SetLogy();
- // fhHitsInCluster->Draw();
- // fhHitsInCluster->Write();
- // if (fPrintToFileOn) gPad->Print(".gif");
- // if (fPrintToFileOn) gPad->Print(".eps");
-
- printf("Total number of points: %i\n", fPointsTotal);
- printf("Points overcounted: %i\n", fPointsOverCounted);
- printf("Points undercounted: %i\n", fPointsUnderCounted);
- }
-
- if (fClusterDeconvQaOn) {
- printf("Signal points: %i\n", fSignalPoints);
- printf("Signal hits: %i\n", fSignalHits);
- }
-
- // TFile* performanceFile = new TFile(fFileName, "recreate");
- //
- // for (Int_t i=0;i<fNstations;i++) {
- // fhPadsTotalR[i]->Write();
- // fhPadsFiredR[i]->Write();
- // fhOccupancyR[i]->Write();
- // }
- //
- // fhPullX->Write();
- // fhPullY->Write();
- //
- // fhPointsInCluster->Write();
- // fhDigisInCluster->Write();
- // fhHitsInCluster->Write();
- //
- // for (Int_t i=0;i<10;i++) fChargeHistos->At(i)->Write();
- //
- // fhNpadsVsS->Write();
- // performanceFile->Close();
-}
-// -------------------------------------------------------------------------
-
-
-// -------------------------------------------------------------------------
-Double_t 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 CbmMuchHitFinderQa::DigitizerQa() {
- Bool_t verbose = (fVerbose > 2);
- 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();
-
- // Get mass
- 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) {
- 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]++;
- }
-
- 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);
- }
-
-
- // Filling generated charge for each point
- // Changing below loop from DigiMatches to Digis
-
- //for (Int_t i=0;i<fDigiMatches->GetEntriesFast();i++){
- // old code
- /*for (Int_t i=0;i<fDigis->GetEntriesFast();i++){
- // Get pad area
- CbmMuchDigi* digi = (CbmMuchDigi*) fDigis->At(i); */
- for (Int_t i = 0; i < fDigiManager->GetNofDigis(ECbmModuleId::kMuch); i++) {
- CbmMuchDigi* digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(i);
- assert(digi);
- //CbmMatch* match = (CbmMatch*) fDigiMatches->At(i);
- 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);
- }
- }
-
- //Filling digitizer performance plots
- 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);
- fhCharge[info->GetStationId()]->Fill(charge);
- fhChargeLog->Fill(log_charge);
- fhChargeEnergyLog->Fill(log_kine, charge);
- fhChargeTrackLength->Fill(length, charge);
-
- if (pdg == 2212)
- fhChargeEnergyLogPr->Fill(log_kine, charge);
- else if (pdg == 211 || pdg == -211)
- fhChargeEnergyLogPi->Fill(log_kine, charge);
- else if (pdg == 11 || pdg == -11)
- fhChargeEnergyLogEl->Fill(log_kine, charge);
-
- if (pdg == 2212)
- fhChargeTrackLengthPr->Fill(length, charge);
- else if (pdg == 211 || pdg == -211)
- fhChargeTrackLengthPi->Fill(length, charge);
- else if (pdg == 11 || pdg == -11)
- fhChargeTrackLengthEl->Fill(length, charge);
-
- if (pdg == 2212 && length > 0.3 && length < 0.32 && kine > 1000
- && kine < 1020)
- fhChargePr_1GeV_3mm->Fill(charge);
- }
-}
-// -------------------------------------------------------------------------
-
-
-// -------------------------------------------------------------------------
-void CbmMuchHitFinderQa::OccupancyQa() {
- // Bool_t verbose = (fVerbose>2);
- // Filling occupancy plots
- // old code
- /*for (Int_t i=0;i<fDigis->GetEntriesFast();i++){
- CbmMuchDigi* digi = (CbmMuchDigi*) fDigis->At(i); */
- 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);
- fhPadsFiredR[CbmMuchAddress::GetStationIndex(address)]->Fill(r0);
- }
-}
-// -------------------------------------------------------------------------
-
-
-// -------------------------------------------------------------------------
-void CbmMuchHitFinderQa::StatisticsQa() {
- // Bool_t verbose = (fVerbose>2);
- Int_t nClusters = fClusters->GetEntriesFast();
- TArrayI hitsInCluster;
- TArrayI pointsInCluster;
- hitsInCluster.Set(nClusters);
- pointsInCluster.Set(nClusters);
- cout << " start Stat QA " << endl;
- for (Int_t i = 0; i < nClusters; i++)
- hitsInCluster[i] = 0;
- for (Int_t i = 0; i < nClusters; i++)
- pointsInCluster[i] = 0;
-
- for (Int_t i = 0; i < fHits->GetEntriesFast(); i++) {
- CbmMuchPixelHit* hit = (CbmMuchPixelHit*) fHits->At(i);
- //cout<<" hit index "<<i<<" plane id "<<hit->GetPlaneId()<<" x "<<hit->GetX()<<" y "<<hit->GetY()<<" z "<<hit->GetZ()<<" cluster Id "<< hit->GetRefId()<<endl;
-
- // cout<<" hit index "<<i<<" plane id "<<hit->GetPlaneId()<<endl;
- Int_t clusterId = hit->GetRefId();
- hitsInCluster[clusterId]++;
- }
-
- for (Int_t i = 0; i < nClusters; i++) {
- CbmMuchCluster* cluster = (CbmMuchCluster*) fClusters->At(i);
-
- map<Int_t, Int_t> map_points;
- Int_t nDigis = cluster->GetNofDigis();
-
- auto address = cluster->GetAddress();
- auto StationIndex = CbmMuchAddress::GetStationIndex(address);
- //cout<<" station index "<<StationIndex<<endl;
-
- fhDigisInCluster[StationIndex]->Fill(nDigis);
-
- for (Int_t digiId = 0; digiId < nDigis; digiId++) {
- Int_t index = cluster->GetDigi(digiId);
- //Access Match from CbmDigi only
- // CbmMuchDigi* digi= (CbmMuchDigi*) fDigis->At(index);
- //CbmMatch* match = (CbmMatch*) fDigiMatches->At(index);
- CbmMatch* match =
- (CbmMatch*) fDigiManager->GetMatch(ECbmModuleId::kMuch, index);
- Int_t nPoints = match->GetNofLinks();
- for (Int_t iRefPoint = 0; iRefPoint < nPoints; iRefPoint++) {
- Int_t pointId = match->GetLink(iRefPoint).GetIndex();
- map_points[pointId] = 1;
- }
- }
- pointsInCluster[i] = map_points.size();
- map_points.clear();
- }
-
-
- for (Int_t i = 0; i < nClusters; i++) {
- // added
- CbmMuchCluster* cluster = (CbmMuchCluster*) fClusters->At(i);
- auto address = cluster->GetAddress();
- auto StationIndex = CbmMuchAddress::GetStationIndex(address);
- //cout<<" station index "<<StationIndex<<endl;
- /// end add
-
- Int_t nPts = pointsInCluster[i];
- Int_t nHits = hitsInCluster[i];
- fhPointsInCluster[StationIndex]->Fill(nPts);
- // fhPointsInCluster->Fill(nPts);
- fhHitsInCluster[StationIndex]->Fill(nHits);
- if (nPts > nHits) fPointsUnderCounted += (nPts - nHits);
- if (nHits > nPts) fPointsOverCounted += (nHits - nPts);
- fPointsTotal += nPts;
- }
-}
-// -------------------------------------------------------------------------
-
-
-// -------------------------------------------------------------------------
-void CbmMuchHitFinderQa::PullsQa() {
- Bool_t verbose = (fVerbose > 2);
- // Filling residuals and pools for hits at the first layer
- for (Int_t i = 0; i < fHits->GetEntriesFast(); i++) {
- CbmMuchPixelHit* hit = (CbmMuchPixelHit*) fHits->At(i);
- // Select hits from the first station only
- Int_t iStation = CbmMuchAddress::GetStationIndex(hit->GetAddress());
- Int_t iLayer = CbmMuchAddress::GetLayerIndex(hit->GetAddress());
- // if((iStation !=0 && iStation !=1))continue;
- // if((iStation !=0 && iStation !=1))continue;
- // cout<<" PULLS QA STATION INDEX "<<iStation<<endl;
- //Earlier finding for only one station
- //if(!(iStation == 0)) continue;
- // if(!(iStation == 3 && iLayer == 0)) continue;
- if (verbose)
- printf(" Hit %i, station %i, layer %i ", i, iStation, iLayer);
-
- Int_t clusterId = hit->GetRefId();
- Bool_t hit_unique = 1;
- for (Int_t j = i + 1; j < fHits->GetEntriesFast(); j++) {
- CbmMuchPixelHit* hit1 = (CbmMuchPixelHit*) fHits->At(j);
- //if (hit1->GetStationNr()>stationNr) break;
- if (hit1->GetRefId() == clusterId) {
- hit_unique = 0;
- break;
- }
- }
- if (verbose) printf("hit_unique=%i", hit_unique);
- if (!hit_unique) {
- if (verbose) printf("\n");
- continue;
- }
-
- // Select hits with clusters formed by only one MC Point
- CbmMuchCluster* cluster = (CbmMuchCluster*) fClusters->At(clusterId);
- Bool_t point_unique = 1;
- Int_t pointId = -1;
-
- /* Double_t xmax = 0;
- Double_t xmin = 0;
- Double_t ymax = 0;
- Double_t ymin = 0;
- Double_t dxmin = 0;
- Double_t dymin = 0;
- */
- // if (cluster->GetNDigis()>1) {if (verbose) printf("\n"); continue;}
- // cout<<" digi number per cluster "<<cluster->GetNofDigis()<<endl;
- for (Int_t digiId = 0; digiId < cluster->GetNofDigis(); digiId++) {
- Int_t index = cluster->GetDigi(digiId);
- // printf("%i\n",index);
- // CbmMuchDigi* digi= (CbmMuchDigi*) fDigis->At(index);
- CbmMuchDigi* digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(index);
- // cout<<" check 1"<<endl;
- if (!digi) continue;
- if (index < 0) continue;
- //CbmMatch* match = (CbmMatch*) fDigiMatches->At(index);
- CbmMatch* match =
- (CbmMatch*) fDigiManager->GetMatch(ECbmModuleId::kMuch, index);
- if (!match) continue; ///added ekata
- // cout<<" check 2 "<<endl;
- //CbmMuchDigiMatch* match = (CbmMuchDigiMatch*) fDigiMatches->At(index);
- // Not unique if the pad has several mcPoint references
- if (verbose) printf(" n=%i", match->GetNofLinks());
- if (match->GetNofLinks() == 0) {
- printf(" noise hit");
- point_unique = 0;
- break;
- }
- if (match->GetNofLinks() > 1) {
- point_unique = 0;
- break;
- }
- Int_t currentPointId = match->GetLink(0).GetIndex();
- CbmMuchModuleGem* module =
- (CbmMuchModuleGem*) fGeoScheme->GetModuleByDetId(digi->GetAddress());
- if (!module) continue;
- // CbmMuchPad* pad = module->GetPad(digi->GetAddress());//fGeoScheme->GetPadByDetId(digi->GetDetectorId(), digi->GetChannelId());
- /*Double_t x = pad->GetX();
- Double_t y = pad->GetY();
- Double_t dx = pad->GetDx();
- Double_t dy = pad->GetDy();
- if (digiId==0 || dxmin>dx ) dxmin=dx;
- if (digiId==0 || dymin>dy ) dymin=dy;
- if (digiId==0 || xmin>x-dx/2) xmin=x-dx/2;
- if (digiId==0 || xmax<x+dx/2) xmax=x+dx/2;
- if (digiId==0 || ymin>y-dy/2) ymin=y-dy/2;
- if (digiId==0 || ymax<y+dy/2) ymax=y+dy/2;*/
- if (digiId == 0) {
- pointId = currentPointId;
- continue;
- }
- // Not unique if mcPoint references differ for different digis
- if (currentPointId != pointId) {
- point_unique = 0;
- break;
- }
- }
-
-
- if (verbose) printf(" point_unique=%i", point_unique);
- if (!point_unique) {
- if (verbose) printf("\n");
- continue;
- }
- //printf(" %f %f %f %f %f %f\n",xmin,xmax,ymin,ymax,dxmin,dymin);
- // Int_t nPadsX=Int_t((xmax-xmin)/dxmin);
- // Int_t nPadsY=Int_t((ymax-ymin)/dymin);
- //printf("nPadsX=%i nPadsY=%i\n",nPadsX,nPadsY);
-
- if (verbose) printf(" pointId=%i", pointId);
- CbmMuchPoint* point = (CbmMuchPoint*) fPoints->At(pointId);
-
- Double_t xMC = 0.5 * (point->GetXIn() + point->GetXOut());
- Double_t yMC = 0.5 * (point->GetYIn() + point->GetYOut());
- Double_t tMC = point->GetTime();
- // cout<<" MC point time "<<tMC<<" z "<<point->GetZ()<<endl;
- Double_t xRC = hit->GetX();
- Double_t yRC = hit->GetY();
- Double_t dx = hit->GetDx();
- Double_t dy = hit->GetDy();
-
- Double_t tRC = hit->GetTime();
- Double_t dt = hit->GetTimeError();
- // cout<<" Rec Hit time "<<tRC<<endl;
-
- if (dx < 1.e-10) {
- printf("Anomalously small dx\n");
- continue;
- }
- if (dy < 1.e-10) {
- printf("Anomalously small dy\n");
- continue;
- }
- fhPullX->Fill((xRC - xMC) / dx);
- fhPullY->Fill((yRC - yMC) / dy);
- fhPullT->Fill((tRC - tMC) / dt);
- fhResidualX->Fill((xRC - xMC));
- fhResidualY->Fill((yRC - yMC));
- fhResidualT->Fill((tRC - tMC));
- // Int_t nDigis = cluster->GetNDigis();
- // if (nDigis<=fNTimingPulls) fhPullT[nDigis]->Fill((tRC-tMC)/dt);
-
- if (verbose) printf("\n");
-
- // Int_t index = cluster->GetDigiIndex(0);
- //
- // // printf("index=%i\n",index);
- // CbmMuchDigi* digi = (CbmMuchDigi*) fDigis->At(index);
- //
- //
- // Int_t padSizeX = Int_t(TMath::Log2(dxmin/fPadMinLx));
- // Int_t padSizeY = Int_t(TMath::Log2(dymin/fPadMinLy));
- // if (padSizeX>=fnPadSizesX || padSizeX<0) { printf("wrong x pad size\n"); continue; }
- // if (padSizeY>=fnPadSizesY || padSizeY<0) { printf("wrong y pad size\n"); continue; }
- // if (nPadsX==1 && nPadsY==1) fhPullXpads1[padSizeX]->Fill((xRC-xMC)/dx);
- // if (nPadsY==1 && nPadsX==1) fhPullYpads1[padSizeY]->Fill((yRC-yMC)/dy);
- // if (nPadsX==2 && nPadsY==1) fhPullXpads2[padSizeX]->Fill((xRC-xMC)/dx);
- // if (nPadsY==2 && nPadsX==1) fhPullYpads2[padSizeY]->Fill((yRC-yMC)/dy);
- // if (nPadsX==3 && nPadsY==1) fhPullXpads3[padSizeX]->Fill((xRC-xMC)/dx);
- // if (nPadsY==3 && nPadsX==1) fhPullYpads3[padSizeY]->Fill((yRC-yMC)/dy);
- }
-}
-// -------------------------------------------------------------------------
-
-
-// -------------------------------------------------------------------------
-void CbmMuchHitFinderQa::ClusterDeconvQa() {
- Int_t nPoints = fPoints->GetEntriesFast();
- // Int_t nMatches = fDigiMatches->GetEntriesFast();
- Int_t nClusters = fClusters->GetEntriesFast();
- vector<Int_t> pIndices;
- vector<Int_t>::iterator it;
-
- for (Int_t iPoint = 0; iPoint < nPoints; ++iPoint) {
- if (IsSignalPoint(iPoint)) fSignalPoints++;
- }
-
- for (Int_t iCluster = 0; iCluster < nClusters; ++iCluster) {
- CbmMuchCluster* cluster = (CbmMuchCluster*) fClusters->At(iCluster);
- if (!cluster) continue;
- Int_t nDigis = cluster->GetNofDigis();
- for (Int_t id = 0; id < nDigis; ++id) {
- Int_t iDigi = cluster->GetDigi(id);
- // CbmMuchDigi* digi= (CbmMuchDigi*) fDigis->At(iDigi);
- //CbmMatch* match = (CbmMatch*) fDigiMatches->At(iDigi);
- CbmMatch* match =
- (CbmMatch*) fDigiManager->GetMatch(ECbmModuleId::kMuch, iDigi);
-
- //CbmMuchDigiMatch* match = (CbmMuchDigiMatch*) fDigiMatches->At(iDigi);
- if (!match) continue;
- for (Int_t ip = 0; ip < match->GetNofLinks(); ++ip) {
- Int_t iPoint = match->GetLink(ip).GetIndex();
- it = find(pIndices.begin(), pIndices.end(), iPoint);
- if (it != pIndices.end()) continue;
- pIndices.push_back(iPoint);
- if (IsSignalPoint(iPoint)) fSignalHits++;
- }
- }
- }
-}
-// -------------------------------------------------------------------------
-
-Bool_t CbmMuchHitFinderQa::IsSignalPoint(Int_t iPoint) {
- Int_t nPoints = fPoints->GetEntriesFast();
- Int_t nTracks = fMCTracks->GetEntriesFast();
- CbmMuchPoint* point = (iPoint < 0 || iPoint > nPoints - 1)
- ? NULL
- : (CbmMuchPoint*) fPoints->At(iPoint);
- if (!point) return kFALSE;
- Int_t iTrack = point->GetTrackID();
- CbmMCTrack* track = (iTrack < 0 || iTrack > nTracks - 1)
- ? NULL
- : (CbmMCTrack*) fMCTracks->At(iTrack);
- if (!track) return kFALSE;
- if (iTrack != 0 && iTrack != 1)
- return kFALSE; // Signal tracks must be fist ones
- // Verify if it is a signal muon
- if (track->GetMotherId() < 0) { // No mother for signal
- Int_t pdgCode = track->GetPdgCode();
- if (TMath::Abs(pdgCode) == 13) { // If it is a muon
- return kTRUE;
- }
- }
- return kFALSE;
-}
-
-// -------------------------------------------------------------------------
-Int_t CbmMuchHitFinderQa::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 CbmMuchHitFinderQa::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 CbmMuchHitFinderQa::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 CbmMuchHitFinderQa::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 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);
- }
-}
-// -------------------------------------------------------------------------
-
-ClassImp(CbmMuchHitFinderQa)
diff --git a/reco/detectors/much/CbmMuchHitFinderQa.h b/reco/detectors/much/CbmMuchHitFinderQa.h
deleted file mode 100644
index ec6eb87f3c..0000000000
--- a/reco/detectors/much/CbmMuchHitFinderQa.h
+++ /dev/null
@@ -1,178 +0,0 @@
-// -------------------------------------------------------------------------
-// ----- CbmMuchHitFinderQa header file -----
-// ----- Modified 01/18 by Vikas Singhal -----
-// ----- Created 16/11/07 by E. Kryshen -----
-// -------------------------------------------------------------------------
-
-#ifndef CbmMuchHitFinderQa_H
-#define CbmMuchHitFinderQa_H
-
-#include "FairTask.h"
-
-#include "CbmMuchDigi.h"
-#include "TClonesArray.h"
-#include "TH1.h"
-#include "TH2.h"
-#include "TString.h"
-
-class CbmDigiManager;
-class CbmMuchGeoScheme;
-class TObjArray;
-class TVector2;
-
-Double_t LandauMPV(Double_t* x, Double_t* par);
-Double_t MPV_n_e(Double_t Tkin, Double_t mass);
-
-class CbmMuchHitFinderQa : public FairTask {
-
-public:
- CbmMuchHitFinderQa(const char* name = "MuchHitFinderQa", Int_t verbose = 1);
- virtual ~CbmMuchHitFinderQa();
- virtual InitStatus Init();
- virtual void Exec(Option_t* option);
- virtual void FinishTask();
- virtual void SetParContainers();
- void SetGeoFileName(TString fileName) { fGeoFileName = fileName; }
- void SetPerformanceFileName(TString fileName) { fFileName = fileName; }
- void SetGeometryID(Int_t flag) { fFlag = flag; }
-
- void SetPullsQa(Bool_t on) { fPullsQaOn = on; }
- void SetOccupancyQa(Bool_t on) { fOccupancyQaOn = on; }
- void SetDigitizerQa(Bool_t on) { fDigitizerQaOn = on; }
- void SetStatisticsQa(Bool_t on) { fStatisticsQaOn = on; }
- void SetClusterDeconvQa(Bool_t on) { fClusterDeconvQaOn = on; }
-
-
- void SetPrintToFile(Bool_t on) { fPrintToFileOn = on; }
-
-protected:
- /* Analysis of hit uncertainty (pull) distributions
- * as function of pad size and cluster shape
- */
- void PullsQa();
-
- /* Occupance analysis - all pads,fired pads,
- * and fired/all distributions as functions of radius
- */
- void OccupancyQa();
-
- /* DigitizerQa - analysis of digitizer performance - charge distributions
- * for tracks. Track length distrivutions. Statistics on particle types
- */
- void DigitizerQa();
-
- /* Information on clusters - number of pads in a cluster, number of points, contributed to
- * a cluster, number of hits, created from a cluster
- */
- void StatisticsQa();
-
- /* Number of points and hits in a cluster for signal muons (MotherId = 0) */
- void ClusterDeconvQa();
-
-private:
- CbmMuchGeoScheme* fGeoScheme;
- TString fGeoFileName;
- TString fFileName;
- Int_t fSignalPoints; // Number of signal MC points
- Int_t fSignalHits; // Number of signal hits
- Int_t fVerbose;
- Int_t fEvent;
- Int_t fFlag;
- TClonesArray* fPoints;
- TClonesArray* fDigis;
- CbmDigiManager* fDigiManager;
- TClonesArray* fDigiMatches;
- TClonesArray* fClusters;
- TClonesArray* fHits;
- TClonesArray* fMCTracks;
- TClonesArray* fPointInfos;
-
- Int_t fNstations;
-
- TObjArray* fChargeHistos;
- TH2D* fhChargeEnergyLog;
- TH2D* fhChargeEnergyLogPi;
- TH2D* fhChargeEnergyLogPr;
- TH2D* fhChargeEnergyLogEl;
- TH2D* fhChargeTrackLength;
- TH2D* fhChargeTrackLengthPi;
- TH2D* fhChargeTrackLengthPr;
- TH2D* fhChargeTrackLengthEl;
- // TH1D* fhCharge;
- TH1D* fhChargeLog;
- TH1D* fhChargePr_1GeV_3mm;
- TH2D* fhNpadsVsS;
-
- TH1D** fhCharge;
- TH1D** fhOccupancyR;
- TH1D** fhPadsTotalR;
- TH1D** fhPadsFiredR;
- TH1D** fhPullXpads1;
- TH1D** fhPullYpads1;
- TH1D** fhPullXpads2;
- TH1D** fhPullYpads2;
- TH1D** fhPullXpads3;
- TH1D** fhPullYpads3;
- Int_t fnPadSizesX;
- Int_t fnPadSizesY;
-
- //TH1D** fhPullT;
- Int_t fNTimingPulls;
- //1D Histogram for PULL Distribution
- TH1D* fhPullX;
- TH1D* fhPullY;
- TH1D* fhPullT;
-
- //1D Histogram for Residual Distribution
- TH1D* fhResidualX;
- TH1D* fhResidualY;
- TH1D* fhResidualT;
-
- TH1I** fhPointsInCluster;
- TH1I** fhDigisInCluster;
- TH1I** fhHitsInCluster;
-
- Int_t* fNall; // number of all tracks at the first station
- Int_t* fNpr; // number of protons at the first station
- Int_t* fNpi; // number of pions at the first station
- Int_t* fNel; // number of electrons at the first station
- Int_t* fNmu; // number of muons at the first station
- Int_t* fNka; // number of kaons at the first station
- Int_t* fNprimary; // number of primary tracks at the first station
- Int_t* fNsecondary; // number of secondary tracks at the first station
-
- Int_t fPointsTotal;
- Int_t fPointsUnderCounted;
- Int_t fPointsOverCounted;
-
- Bool_t fOccupancyQaOn;
- Bool_t fPullsQaOn;
- Bool_t fDigitizerQaOn;
- Bool_t fStatisticsQaOn;
- Bool_t fClusterDeconvQaOn;
-
- Bool_t fPrintToFileOn;
-
- Double_t fPadMinLx;
- Double_t fPadMinLy;
- Double_t fPadMaxLx;
- Double_t fPadMaxLy;
-
- FILE* pointsFile;
- FILE* padsFile;
-
- /** Defines whether the point with the given index is signal point. **/
- Bool_t IsSignalPoint(Int_t iPoint);
-
- Int_t GetNChannels(Int_t iStation);
- Int_t GetNSectors(Int_t iStation);
- TVector2 GetMinPadSize(Int_t iStation);
- TVector2 GetMaxPadSize(Int_t iStation);
-
- CbmMuchHitFinderQa(const CbmMuchHitFinderQa&);
- CbmMuchHitFinderQa& operator=(const CbmMuchHitFinderQa&);
-
- ClassDef(CbmMuchHitFinderQa, 1)
-};
-
-#endif
diff --git a/reco/detectors/much/qa/CbmMuchHitFinderQa.cxx b/reco/detectors/much/qa/CbmMuchHitFinderQa.cxx
new file mode 100644
index 0000000000..6b397aaeb3
--- /dev/null
+++ b/reco/detectors/much/qa/CbmMuchHitFinderQa.cxx
@@ -0,0 +1,692 @@
+// -------------------------------------------------------------------------
+// ----- CbmMuchHitProducerQa source file -----
+// ----- Modified since 21/06/2019 by Ekata Nandy (ekata@vecc.gov.in) -Inclusion of RPC detector type
+// ----- Modified 02/18 by Vikas Singhal -----
+// ----- Created 16/11/07 by E. Kryshen -----
+// -------------------------------------------------------------------------
+// AUG 18 RELEASE VERSION
+
+#include "CbmMuchHitFinderQa.h"
+
+#include "CbmDefs.h"
+#include "CbmDigiManager.h"
+#include "CbmLink.h"
+#include "CbmMCTrack.h"
+#include "CbmMatch.h"
+#include "CbmMuchAddress.h"
+#include "CbmMuchCluster.h"
+#include "CbmMuchDigi.h"
+#include "CbmMuchGeoScheme.h"
+#include "CbmMuchModuleGem.h"
+#include "CbmMuchPixelHit.h"
+#include "CbmMuchPoint.h"
+#include "CbmQaCanvas.h"
+
+#include "FairRootManager.h"
+#include <FairSink.h>
+#include <FairTask.h>
+#include <Logger.h>
+
+#include "TArrayI.h"
+#include "TF1.h"
+#include "TFile.h"
+#include "TObjArray.h"
+#include "TParameter.h"
+#include "TPaveStats.h"
+#include "TStyle.h"
+#include <TClonesArray.h>
+#include <TGenericClassInfo.h>
+#include <TH1.h>
+#include <TMathBase.h>
+#include <TString.h>
+
+#include <boost/exception/exception.hpp>
+#include <boost/type_index/type_index_facade.hpp>
+
+#include <algorithm>
+#include <iostream>
+#include <map>
+#include <stdio.h>
+#include <vector>
+
+using std::cout;
+using std::endl;
+using std::map;
+using std::vector;
+// -------------------------------------------------------------------------
+CbmMuchHitFinderQa::CbmMuchHitFinderQa(const char* name, Int_t verbose)
+ : FairTask(name, verbose)
+ , fGeoScheme(CbmMuchGeoScheme::Instance())
+ , fGeoFileName()
+ , fFileName()
+ , fVerbose(verbose)
+ , fOutFolder("MuchHitFinderQA", "MuchHitFinderQA")
+ , fhPointsInCluster()
+ , fhDigisInCluster()
+ , fhHitsPerCluster()
+ , fNevents("nEvents", 0)
+ , fSignalPoints("SignalPoints", 0)
+ , fSignalHits("SignalHits", 0)
+ , fPointsTotal("PointsTotal", 0)
+ , fPointsUnderCounted("PointsUnderCounted", 0)
+ , fPointsOverCounted("PointsOverCounted", 0) {}
+// -------------------------------------------------------------------------
+
+
+// -------------------------------------------------------------------------
+CbmMuchHitFinderQa::~CbmMuchHitFinderQa() { DeInit(); }
+// -------------------------------------------------------------------------
+
+// -------------------------------------------------------------------------
+void CbmMuchHitFinderQa::DeInit() {
+ fPoints = nullptr;
+ fDigiManager = nullptr;
+ fClusters = nullptr;
+ fHits = nullptr;
+ fMCTracks = nullptr;
+ histFolder = nullptr;
+ fOutFolder.Clear();
+
+ SafeDelete(fhPullX);
+ SafeDelete(fhPullY);
+ SafeDelete(fhPullT);
+ SafeDelete(fhResidualX);
+ SafeDelete(fhResidualY);
+ SafeDelete(fhResidualT);
+
+ for (uint i = 0; i < fhPointsInCluster.size(); i++) {
+ SafeDelete(fhPointsInCluster[i]);
+ }
+ for (uint i = 0; i < fhDigisInCluster.size(); i++) {
+ SafeDelete(fhDigisInCluster[i]);
+ }
+ for (uint i = 0; i < fhHitsPerCluster.size(); i++) {
+ SafeDelete(fhHitsPerCluster[i]);
+ }
+ fhPointsInCluster.clear();
+ fhDigisInCluster.clear();
+ fhHitsPerCluster.clear();
+
+ SafeDelete(fCanvPointsInCluster);
+ SafeDelete(fCanvDigisInCluster);
+ SafeDelete(fCanvHitsPerCluster);
+ SafeDelete(fCanvPull);
+ SafeDelete(fCanvResidual);
+
+ fVerbose = 0;
+ fFlag = 0;
+ fNstations = 0;
+
+ fNevents.SetVal(0);
+ fSignalPoints.SetVal(0);
+ fSignalHits.SetVal(0);
+ fPointsTotal.SetVal(0);
+ fPointsUnderCounted.SetVal(0);
+ fPointsOverCounted.SetVal(0);
+}
+// -------------------------------------------------------------------------
+
+
+// -------------------------------------------------------------------------
+InitStatus CbmMuchHitFinderQa::Init() {
+ DeInit();
+ FairRootManager* fManager = FairRootManager::Instance();
+ fMCTracks = (TClonesArray*) fManager->GetObject("MCTrack");
+ fPoints = (TClonesArray*) fManager->GetObject("MuchPoint");
+ fHits = (TClonesArray*) fManager->GetObject("MuchPixelHit");
+ fClusters = (TClonesArray*) fManager->GetObject("MuchCluster");
+ // Reading Much Digis from CbmMuchDigiManager which are stored as vector
+ fDigiManager = CbmDigiManager::Instance();
+
+ TFile* f = new TFile(fGeoFileName, "R");
+ TObjArray* stations = (TObjArray*) f->Get("stations");
+ fGeoScheme->Init(stations, fFlag);
+
+ if (!fManager) {
+ LOG(fatal) << "Can not find FairRootManager";
+ return kFATAL;
+ }
+ if (!fDigiManager) {
+ LOG(fatal) << "Can not find Much digi manager";
+ return kFATAL;
+ }
+ if (!fGeoScheme) {
+ LOG(fatal) << "Can not find Much geo scheme";
+ return kFATAL;
+ }
+ if (!fMCTracks) {
+ LOG(error) << "No MC tracks found";
+ return kERROR;
+ }
+ if (!fPoints) {
+ LOG(error) << "No MC points found";
+ return kERROR;
+ }
+ if (!fHits) {
+ LOG(error) << "No hits found";
+ return kERROR;
+ }
+ if (!fClusters) {
+ LOG(error) << "No hits found";
+ return kERROR;
+ }
+ histFolder = fOutFolder.AddFolder("hist", "Histogramms");
+
+ fNevents.SetVal(0);
+ fSignalPoints.SetVal(0);
+ fSignalHits.SetVal(0);
+ fPointsTotal.SetVal(0);
+ fPointsUnderCounted.SetVal(0);
+ fPointsOverCounted.SetVal(0);
+
+ histFolder->Add(&fNevents);
+ histFolder->Add(&fSignalPoints);
+ histFolder->Add(&fSignalHits);
+ histFolder->Add(&fPointsTotal);
+ histFolder->Add(&fPointsUnderCounted);
+ histFolder->Add(&fPointsOverCounted);
+
+ fDigiManager->Init();
+ fNstations = fGeoScheme->GetNStations();
+ LOG(debug) << "Init(): fNstations = " << fNstations;
+
+ fhPointsInCluster.resize(fNstations);
+ fhDigisInCluster.resize(fNstations);
+ fhHitsPerCluster.resize(fNstations);
+
+ for (Int_t i = 0; i < fNstations; i++) {
+
+ fhPointsInCluster[i] =
+ new TH1I(Form("hMCPointsInCluster%i", i + 1),
+ Form("MC Points in Cluster : Station %i ", i + 1),
+ 10,
+ 0,
+ 10);
+ fhDigisInCluster[i] =
+ new TH1I(Form("hDigisInCluster%i", i + 1),
+ Form("Digis in Cluster : Station %i ", i + 1),
+ 10,
+ 0,
+ 10);
+ fhHitsPerCluster[i] =
+ new TH1I(Form("hHitsPerCluster%i", i + 1),
+ Form("Hits per Cluster : Station %i ", i + 1),
+ 10,
+ 0,
+ 10);
+ histFolder->Add(fhPointsInCluster[i]);
+ histFolder->Add(fhDigisInCluster[i]);
+ histFolder->Add(fhHitsPerCluster[i]);
+ }
+ gStyle->SetOptStat(1);
+
+ //Pull Distribution
+ fhPullX = new TH1D(
+ "hPullX", "Pull distribution X;(x_{RC} - x_{MC}) / dx_{RC}", 500, -5, 5);
+ fhPullY = new TH1D(
+ "hPullY", "Pull distribution Y;(y_{RC} - y_{MC}) / dy_{RC}", 500, -5, 5);
+ fhPullT = new TH1D(
+ "hPullT", "Pull distribution T;(t_{RC} - t_{MC}) / dt_{RC}", 120, -5, 5);
+
+ //Residual Distribution
+ fhResidualX = new TH1D(
+ "hResidualX", "Residual distribution X;(x_{RC} - x_{MC})(cm)", 500, -3, 3);
+ fhResidualY = new TH1D(
+ "hResidualY", "Residual distribution Y;(y_{RC} - y_{MC})(cm)", 500, -3, 3);
+ fhResidualT = new TH1D("hResidualT",
+ "Residual distribution T;(t_{RC} - t_{MC})(ns)",
+ 140,
+ -17,
+ 17);
+
+ histFolder->Add(fhPullX);
+ histFolder->Add(fhPullY);
+ histFolder->Add(fhPullT);
+ histFolder->Add(fhResidualX);
+ histFolder->Add(fhResidualY);
+ histFolder->Add(fhResidualT);
+
+ fCanvPointsInCluster = new CbmQaCanvas(
+ "cMCPointsInCluster", "MC Points In Cluster", 2 * 400, 2 * 400);
+ fCanvPointsInCluster->Divide2D(fNstations);
+
+ fCanvDigisInCluster =
+ new CbmQaCanvas("cDigisInCluster", "Digis In Cluster", 2 * 400, 2 * 400);
+ fCanvDigisInCluster->Divide2D(fNstations);
+
+ fCanvHitsPerCluster =
+ new CbmQaCanvas("cHitsPerCluster", "Hits Per Cluster", 2 * 400, 2 * 400);
+ fCanvHitsPerCluster->Divide2D(fNstations);
+
+ fCanvPull = new CbmQaCanvas("cPull", "Pull Distribution", 3 * 600, 1 * 400);
+ fCanvPull->Divide2D(3);
+
+ fCanvResidual =
+ new CbmQaCanvas("cResidual", "Residual Distribution", 3 * 600, 1 * 400);
+ fCanvResidual->Divide2D(3);
+
+ fOutFolder.Add(fCanvPointsInCluster);
+ fOutFolder.Add(fCanvDigisInCluster);
+ fOutFolder.Add(fCanvHitsPerCluster);
+ fOutFolder.Add(fCanvPull);
+ fOutFolder.Add(fCanvResidual);
+
+ return kSUCCESS;
+}
+// -------------------------------------------------------------------------
+
+
+// -------------------------------------------------------------------------
+void CbmMuchHitFinderQa::SetParContainers() {
+ // Get run and runtime database
+ // FairRuntimeDb* db = FairRuntimeDb::instance();
+ // if ( ! db ) Fatal("SetParContainers", "No runtime database");
+ // Get MUCH geometry parameter container
+ // fGeoPar = (CbmGeoMuchPar*) db->getContainer("CbmGeoMuchPar");
+}
+// -------------------------------------------------------------------------
+
+
+// -------------------------------------------------------------------------x
+void CbmMuchHitFinderQa::Exec(Option_t*) {
+ fNevents.SetVal(fNevents.GetVal() + 1);
+ LOG(debug) << "Event: " << fNevents.GetVal();
+
+ PullsQa();
+ StatisticsQa();
+ ClusterDeconvQa();
+}
+// -------------------------------------------------------------------------
+
+
+// -------------------------------------------------------------------------
+void CbmMuchHitFinderQa::FinishTask() {
+ printf("FinishTask\n");
+
+ gStyle->SetPaperSize(20, 20);
+
+ for (Int_t i = 0; i < fNstations; i++) {
+ fhPointsInCluster[i]->Scale(1. / fNevents.GetVal());
+ fhDigisInCluster[i]->Scale(1. / fNevents.GetVal());
+ fhHitsPerCluster[i]->Scale(1. / fNevents.GetVal());
+ }
+
+ if (fVerbose > 1) {
+ printf("===================================\n");
+ printf("PullsQa:\n");
+ }
+
+ std::vector<TH1D*> vResHistos;
+ vResHistos.push_back(fhPullX);
+ vResHistos.push_back(fhPullY);
+ vResHistos.push_back(fhPullT);
+ vResHistos.push_back(fhResidualX);
+ vResHistos.push_back(fhResidualY);
+ vResHistos.push_back(fhResidualT);
+
+ for (UInt_t i = 0; i < vResHistos.size(); i++) {
+ TH1D* histo = vResHistos[i];
+ histo->Sumw2();
+ histo->Fit("gaus", "Q");
+ histo->GetFunction("gaus")->SetLineWidth(3);
+ histo->GetFunction("gaus")->SetLineColor(kRed);
+ // histo->SetStats(kTRUE);
+ }
+
+ printf("===================================\n");
+ printf("StatisticsQa:\n");
+ printf("Total number of points: %i\n", fPointsTotal.GetVal());
+ printf("Points overcounted: %i\n", fPointsOverCounted.GetVal());
+ printf("Points undercounted: %i\n", fPointsUnderCounted.GetVal());
+ printf("Signal points: %i\n", fSignalPoints.GetVal());
+ printf("Signal hits: %i\n", fSignalHits.GetVal());
+
+ DrawCanvases();
+
+ FairSink* sink = FairRootManager::Instance()->GetSink();
+ sink->WriteObject(&fOutFolder, nullptr);
+}
+// -------------------------------------------------------------------------
+
+
+// -------------------------------------------------------------------------
+void CbmMuchHitFinderQa::DrawCanvases() {
+
+ for (Int_t i = 0; i < fNstations; i++) {
+ fCanvPointsInCluster->cd(i + 1);
+ fhPointsInCluster[i]->DrawCopy("", "");
+
+ fCanvDigisInCluster->cd(i + 1);
+ fhDigisInCluster[i]->DrawCopy("", "");
+
+ fCanvHitsPerCluster->cd(i + 1);
+ fhHitsPerCluster[i]->DrawCopy("", "");
+ }
+
+ std::vector<TH1D*> vPullHistos;
+ vPullHistos.push_back(fhPullX);
+ vPullHistos.push_back(fhPullY);
+ vPullHistos.push_back(fhPullT);
+
+ for (UInt_t i = 0; i < vPullHistos.size(); i++) {
+ TH1D* histo = vPullHistos[i];
+ fCanvPull->cd(i + 1);
+ histo->Draw(); //necessary to create stats pointer
+ fCanvPull->Update();
+ TPaveStats* st = (TPaveStats*) histo->FindObject("stats");
+ st->SetX1NDC(0.621);
+ st->SetX2NDC(0.940);
+ st->SetY1NDC(0.657);
+ st->SetY2NDC(0.929);
+ st->SetOptStat(1110);
+ st->SetOptFit(11);
+ //st->SetTextSize(0.04);
+ histo->DrawCopy("", "");
+ //version below only changes canvas but not hist folder
+ //TH1* hClone = histo->DrawCopy("", "");
+ //fCanvPull->Update();
+ //TPaveStats *st = (TPaveStats*)hClone->FindObject("stats");
+ //st->SetOptStat(1110);
+ //st->SetOptFit(11);
+ }
+
+ std::vector<TH1D*> vResHistos;
+ vResHistos.push_back(fhResidualX);
+ vResHistos.push_back(fhResidualY);
+ vResHistos.push_back(fhResidualT);
+
+ for (UInt_t i = 0; i < vResHistos.size(); i++) {
+ TH1D* histo = vResHistos[i];
+ fCanvResidual->cd(i + 1);
+ histo->Draw(); //necessary to create stats pointer
+ fCanvResidual->Update();
+ TPaveStats* st = (TPaveStats*) histo->FindObject("stats");
+ st->SetX1NDC(0.621);
+ st->SetX2NDC(0.940);
+ st->SetY1NDC(0.657);
+ st->SetY2NDC(0.929);
+ st->SetOptStat(1110);
+ st->SetOptFit(11);
+ //st->SetTextSize(0.04);
+ histo->DrawCopy("", "");
+ }
+}
+// -------------------------------------------------------------------------
+
+
+// -------------------------------------------------------------------------
+void CbmMuchHitFinderQa::StatisticsQa() {
+ // Bool_t verbose = (fVerbose>2);
+ Int_t nClusters = fClusters->GetEntriesFast();
+ TArrayI hitsInCluster;
+ TArrayI pointsInCluster;
+ hitsInCluster.Set(nClusters);
+ pointsInCluster.Set(nClusters);
+ LOG(debug) << " start Stat QA ";
+ for (Int_t i = 0; i < nClusters; i++)
+ hitsInCluster[i] = 0;
+ for (Int_t i = 0; i < nClusters; i++)
+ pointsInCluster[i] = 0;
+
+ for (Int_t i = 0; i < fHits->GetEntriesFast(); i++) {
+ CbmMuchPixelHit* hit = (CbmMuchPixelHit*) fHits->At(i);
+ //cout<<" hit index "<<i<<" plane id "<<hit->GetPlaneId()<<" x "<<hit->GetX()<<" y "<<hit->GetY()<<" z "<<hit->GetZ()<<" cluster Id "<< hit->GetRefId()<<endl;
+
+ // cout<<" hit index "<<i<<" plane id "<<hit->GetPlaneId()<<endl;
+ Int_t clusterId = hit->GetRefId();
+ hitsInCluster[clusterId]++;
+ }
+
+ for (Int_t i = 0; i < nClusters; i++) {
+ CbmMuchCluster* cluster = (CbmMuchCluster*) fClusters->At(i);
+
+ map<Int_t, Int_t> map_points;
+ Int_t nDigis = cluster->GetNofDigis();
+
+ auto address = cluster->GetAddress();
+ auto StationIndex = CbmMuchAddress::GetStationIndex(address);
+ //cout<<" station index "<<StationIndex<<endl;
+
+ fhDigisInCluster[StationIndex]->Fill(nDigis);
+
+ if (fDigiManager->IsMatchPresent(ECbmModuleId::kMuch)) {
+ for (Int_t digiId = 0; digiId < nDigis; digiId++) {
+ Int_t index = cluster->GetDigi(digiId);
+ //Access Match from CbmDigi only
+ CbmMatch* match =
+ (CbmMatch*) fDigiManager->GetMatch(ECbmModuleId::kMuch, index);
+ if (!match) {
+ LOG(fatal) << "CbmMuchHitFinderQa::StatisticsQa(): Match should be "
+ "present but is null.";
+ return;
+ }
+ Int_t nPoints = match->GetNofLinks();
+ for (Int_t iRefPoint = 0; iRefPoint < nPoints; iRefPoint++) {
+ Int_t pointId = match->GetLink(iRefPoint).GetIndex();
+ map_points[pointId] = 1;
+ }
+ }
+ }
+ pointsInCluster[i] = map_points.size();
+ map_points.clear();
+ }
+
+ for (Int_t i = 0; i < nClusters; i++) {
+ // added
+ CbmMuchCluster* cluster = (CbmMuchCluster*) fClusters->At(i);
+ auto address = cluster->GetAddress();
+ auto StationIndex = CbmMuchAddress::GetStationIndex(address);
+ //cout<<" station index "<<StationIndex<<endl;
+ /// end add
+
+ Int_t nPts = pointsInCluster[i];
+ Int_t nHits = hitsInCluster[i];
+ fhPointsInCluster[StationIndex]->Fill(nPts);
+ fhHitsPerCluster[StationIndex]->Fill(nHits);
+ if (nPts > nHits)
+ fPointsUnderCounted.SetVal(fPointsUnderCounted.GetVal() + (nPts - nHits));
+ if (nHits > nPts)
+ fPointsOverCounted.SetVal(fPointsOverCounted.GetVal() + (nHits - nPts));
+ fPointsTotal.SetVal(fPointsTotal.GetVal() + nPts);
+ }
+}
+// -------------------------------------------------------------------------
+
+
+// -------------------------------------------------------------------------
+void CbmMuchHitFinderQa::PullsQa() {
+ Bool_t verbose = (fVerbose > 2);
+ // Filling residuals and pools for hits at the first layer
+ for (Int_t i = 0; i < fHits->GetEntriesFast(); i++) {
+ CbmMuchPixelHit* hit = (CbmMuchPixelHit*) fHits->At(i);
+ // Select hits from the first station only
+ Int_t iStation = CbmMuchAddress::GetStationIndex(hit->GetAddress());
+ Int_t iLayer = CbmMuchAddress::GetLayerIndex(hit->GetAddress());
+ // if((iStation !=0 && iStation !=1))continue;
+ // if((iStation !=0 && iStation !=1))continue;
+ // cout<<" PULLS QA STATION INDEX "<<iStation<<endl;
+ //Earlier finding for only one station
+ //if(!(iStation == 0)) continue;
+ // if(!(iStation == 3 && iLayer == 0)) continue;
+ if (verbose)
+ printf(" Hit %i, station %i, layer %i ", i, iStation, iLayer);
+
+ Int_t clusterId = hit->GetRefId();
+ Bool_t hit_unique = 1;
+ for (Int_t j = i + 1; j < fHits->GetEntriesFast(); j++) {
+ CbmMuchPixelHit* hit1 = (CbmMuchPixelHit*) fHits->At(j);
+ //if (hit1->GetStationNr()>stationNr) break;
+ if (hit1->GetRefId() == clusterId) {
+ hit_unique = 0;
+ break;
+ }
+ }
+ if (verbose) printf("hit_unique=%i", hit_unique);
+ if (!hit_unique) {
+ if (verbose) printf("\n");
+ continue;
+ }
+
+ // Select hits with clusters formed by only one MC Point
+ CbmMuchCluster* cluster = (CbmMuchCluster*) fClusters->At(clusterId);
+ Bool_t point_unique = 1;
+ Int_t pointId = -1;
+
+ for (Int_t digiId = 0; digiId < cluster->GetNofDigis(); digiId++) {
+ Int_t index = cluster->GetDigi(digiId);
+ // printf("%i\n",index);
+ CbmMuchDigi* digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(index);
+ // cout<<" check 1"<<endl;
+ if (!digi) {
+ LOG(fatal) << "CbmMuchHitFinderQa::PullsQa(): Error, digi not found.";
+ return;
+ }
+ if (index < 0) {
+ LOG(fatal)
+ << "CbmMuchHitFinderQa::PullsQa(): Error, index out of bounds.";
+ return;
+ }
+ CbmMatch* match =
+ (CbmMatch*) fDigiManager->GetMatch(ECbmModuleId::kMuch, index);
+ if (!match) {
+ LOG(fatal) << "CbmMuchHitFinderQa::PullsQa(): Error, match not found.";
+ return;
+ }
+ // Not unique if the pad has several mcPoint references
+ if (verbose) printf(" n=%i", match->GetNofLinks());
+ if (match->GetNofLinks() == 0) {
+ printf(" noise hit");
+ point_unique = 0;
+ break;
+ }
+ if (match->GetNofLinks() > 1) {
+ point_unique = 0;
+ break;
+ }
+ Int_t currentPointId = match->GetLink(0).GetIndex();
+
+ CbmMuchModuleGem* module =
+ (CbmMuchModuleGem*) fGeoScheme->GetModuleByDetId(digi->GetAddress());
+ if (!module) {
+ LOG(fatal) << "CbmMuchHitFinderQa::PullsQa(): Error, module not found.";
+ return;
+ }
+ if (digiId == 0) {
+ pointId = currentPointId;
+ continue;
+ }
+ // Not unique if mcPoint references differ for different digis
+ if (currentPointId != pointId) {
+ point_unique = 0;
+ break;
+ }
+ }
+
+ if (verbose) printf(" point_unique=%i", point_unique);
+ if (!point_unique) {
+ if (verbose) printf("\n");
+ continue;
+ }
+
+ if (verbose) printf(" pointId=%i", pointId);
+ CbmMuchPoint* point = (CbmMuchPoint*) fPoints->At(pointId);
+
+ Double_t xMC = 0.5 * (point->GetXIn() + point->GetXOut());
+ Double_t yMC = 0.5 * (point->GetYIn() + point->GetYOut());
+ Double_t tMC = point->GetTime();
+ // cout<<" MC point time "<<tMC<<" z "<<point->GetZ()<<endl;
+ Double_t xRC = hit->GetX();
+ Double_t yRC = hit->GetY();
+ Double_t dx = hit->GetDx();
+ Double_t dy = hit->GetDy();
+
+ Double_t tRC = hit->GetTime();
+ Double_t dt = hit->GetTimeError();
+ // cout<<" Rec Hit time "<<tRC<<endl;
+
+ if (dx < 1.e-10) {
+ printf("Anomalously small dx\n");
+ continue;
+ }
+ if (dy < 1.e-10) {
+ printf("Anomalously small dy\n");
+ continue;
+ }
+ fhPullX->Fill((xRC - xMC) / dx);
+ fhPullY->Fill((yRC - yMC) / dy);
+ fhPullT->Fill((tRC - tMC) / dt);
+ fhResidualX->Fill((xRC - xMC));
+ fhResidualY->Fill((yRC - yMC));
+ fhResidualT->Fill((tRC - tMC));
+
+ if (verbose) printf("\n");
+ }
+}
+// -------------------------------------------------------------------------
+
+
+// -------------------------------------------------------------------------
+void CbmMuchHitFinderQa::ClusterDeconvQa() {
+ Int_t nPoints = fPoints->GetEntriesFast();
+ Int_t nClusters = fClusters->GetEntriesFast();
+ vector<Int_t> pIndices;
+ vector<Int_t>::iterator it;
+
+ for (Int_t iPoint = 0; iPoint < nPoints; ++iPoint) {
+ if (IsSignalPoint(iPoint)) fSignalPoints.SetVal(fSignalPoints.GetVal() + 1);
+ }
+
+ for (Int_t iCluster = 0; iCluster < nClusters; ++iCluster) {
+ CbmMuchCluster* cluster = (CbmMuchCluster*) fClusters->At(iCluster);
+ if (!cluster) {
+ LOG(fatal)
+ << "CbmMuchHitFinderQa::ClusterDeconvQa(): Error, cluster not found.";
+ return;
+ }
+ Int_t nDigis = cluster->GetNofDigis();
+ for (Int_t id = 0; id < nDigis; ++id) {
+ Int_t iDigi = cluster->GetDigi(id);
+ CbmMatch* match =
+ (CbmMatch*) fDigiManager->GetMatch(ECbmModuleId::kMuch, iDigi);
+ if (!match) {
+ LOG(fatal)
+ << "CbmMuchHitFinderQa::ClusterDeconvQa(): Error, match not found.";
+ return;
+ }
+ for (Int_t ip = 0; ip < match->GetNofLinks(); ++ip) {
+ Int_t iPoint = match->GetLink(ip).GetIndex();
+ it = find(pIndices.begin(), pIndices.end(), iPoint);
+ if (it != pIndices.end()) continue;
+ pIndices.push_back(iPoint);
+ if (IsSignalPoint(iPoint)) fSignalHits.SetVal(fSignalHits.GetVal() + 1);
+ }
+ }
+ }
+}
+// -------------------------------------------------------------------------
+
+Bool_t CbmMuchHitFinderQa::IsSignalPoint(Int_t iPoint) {
+ Int_t nPoints = fPoints->GetEntriesFast();
+ Int_t nTracks = fMCTracks->GetEntriesFast();
+ CbmMuchPoint* point = (iPoint < 0 || iPoint > nPoints - 1)
+ ? NULL
+ : (CbmMuchPoint*) fPoints->At(iPoint);
+ if (!point) return kFALSE;
+ Int_t iTrack = point->GetTrackID();
+ CbmMCTrack* track = (iTrack < 0 || iTrack > nTracks - 1)
+ ? NULL
+ : (CbmMCTrack*) fMCTracks->At(iTrack);
+ if (!track) return kFALSE;
+ if (iTrack != 0 && iTrack != 1)
+ return kFALSE; // Signal tracks must be fist ones
+ // Verify if it is a signal muon
+ if (track->GetMotherId() < 0) { // No mother for signal
+ Int_t pdgCode = track->GetPdgCode();
+ if (TMath::Abs(pdgCode) == 13) { // If it is a muon
+ return kTRUE;
+ }
+ }
+ return kFALSE;
+}
+
+ClassImp(CbmMuchHitFinderQa)
diff --git a/reco/detectors/much/qa/CbmMuchHitFinderQa.h b/reco/detectors/much/qa/CbmMuchHitFinderQa.h
new file mode 100644
index 0000000000..a61a438f11
--- /dev/null
+++ b/reco/detectors/much/qa/CbmMuchHitFinderQa.h
@@ -0,0 +1,112 @@
+// -------------------------------------------------------------------------
+// ----- CbmMuchHitFinderQa header file -----
+// ----- Modified 01/18 by Vikas Singhal -----
+// ----- Created 16/11/07 by E. Kryshen -----
+// -------------------------------------------------------------------------
+
+#ifndef CbmMuchHitFinderQa_H
+#define CbmMuchHitFinderQa_H
+
+#include "FairTask.h"
+
+#include "TParameter.h"
+#include "TString.h"
+#include <Rtypes.h>
+#include <RtypesCore.h>
+#include <TFolder.h>
+
+class CbmDigiManager;
+class CbmMuchGeoScheme;
+class CbmQaCanvas;
+class TBuffer;
+class TClass;
+class TClonesArray;
+class TH1D;
+class TH1I;
+class TMemberInspector;
+
+class CbmMuchHitFinderQa : public FairTask {
+
+public:
+ CbmMuchHitFinderQa(const char* name = "MuchHitFinderQa", Int_t verbose = 1);
+ virtual ~CbmMuchHitFinderQa();
+ virtual InitStatus Init();
+ virtual void Exec(Option_t* option);
+ virtual void FinishTask();
+ virtual void SetParContainers();
+ void SetGeoFileName(TString fileName) { fGeoFileName = fileName; }
+ void SetPerformanceFileName(TString fileName) { fFileName = fileName; }
+ void SetGeometryID(Int_t flag) { fFlag = flag; }
+
+protected:
+ /* Analysis of hit uncertainty (pull) distributions
+ * as function of pad size and cluster shape
+ */
+ void PullsQa();
+
+ /* Information on clusters - number of pads in a cluster, number of points, contributed to
+ * a cluster, number of hits, created from a cluster
+ */
+ void StatisticsQa();
+
+ /* Number of points and hits in a cluster for signal muons (MotherId = 0) */
+ void ClusterDeconvQa();
+
+private:
+ void DeInit();
+ void DrawCanvases();
+
+ CbmMuchGeoScheme* fGeoScheme;
+ TString fGeoFileName;
+ TString fFileName;
+ Int_t fVerbose = 0;
+ Int_t fFlag = 0;
+ TClonesArray* fPoints = nullptr;
+ CbmDigiManager* fDigiManager = nullptr;
+
+ TFolder fOutFolder; /// output folder with histos and canvases
+ TFolder* histFolder; /// subfolder for histograms
+
+ TClonesArray* fClusters = nullptr;
+ TClonesArray* fHits = nullptr;
+ TClonesArray* fMCTracks = nullptr;
+
+ Int_t fNstations = 0;
+
+ //1D Histogram for PULL Distribution
+ TH1D* fhPullX = nullptr;
+ TH1D* fhPullY = nullptr;
+ TH1D* fhPullT = nullptr;
+
+ //1D Histogram for Residual Distribution
+ TH1D* fhResidualX = nullptr;
+ TH1D* fhResidualY = nullptr;
+ TH1D* fhResidualT = nullptr;
+
+ std::vector<TH1I*> fhPointsInCluster;
+ std::vector<TH1I*> fhDigisInCluster;
+ std::vector<TH1I*> fhHitsPerCluster;
+
+ CbmQaCanvas* fCanvPointsInCluster = nullptr;
+ CbmQaCanvas* fCanvDigisInCluster = nullptr;
+ CbmQaCanvas* fCanvHitsPerCluster = nullptr;
+ CbmQaCanvas* fCanvPull = nullptr;
+ CbmQaCanvas* fCanvResidual = nullptr;
+
+ TParameter<int> fNevents; /// number of processed events
+ TParameter<int> fSignalPoints; // Number of signal MC points
+ TParameter<int> fSignalHits; // Number of signal hits
+ TParameter<int> fPointsTotal;
+ TParameter<int> fPointsUnderCounted;
+ TParameter<int> fPointsOverCounted;
+
+ /** Defines whether the point with the given index is signal point. **/
+ Bool_t IsSignalPoint(Int_t iPoint);
+
+ CbmMuchHitFinderQa(const CbmMuchHitFinderQa&);
+ CbmMuchHitFinderQa& operator=(const CbmMuchHitFinderQa&);
+
+ ClassDef(CbmMuchHitFinderQa, 1)
+};
+
+#endif
diff --git a/sim/detectors/much/qa/CbmMuchDigitizerQa.cxx b/sim/detectors/much/qa/CbmMuchDigitizerQa.cxx
index 5d62efda00..967038b0c4 100644
--- a/sim/detectors/much/qa/CbmMuchDigitizerQa.cxx
+++ b/sim/detectors/much/qa/CbmMuchDigitizerQa.cxx
@@ -29,6 +29,7 @@
#include "TGraph.h"
#include "TH1.h"
#include "TH2.h"
+#include "TParameter.h"
#include "TString.h"
#include "TStyle.h"
#include <FairRootManager.h>
@@ -46,7 +47,6 @@
#include <math.h>
#include <vector>
-using std::cout;
using std::endl;
using std::vector;
@@ -61,10 +61,8 @@ ClassImp(CbmMuchDigitizerQa);
// -------------------------------------------------------------------------
CbmMuchDigitizerQa::CbmMuchDigitizerQa(const char* name, Int_t verbose)
: FairTask(name, verbose)
- , fGeoScheme(nullptr)
- , fDigiManager(nullptr)
- , fPointInfos(new TClonesArray("CbmMuchPointInfo", 10))
, fOutFolder("MuchDigiQA", "MuchDigitizerQA")
+ , fNevents("nEvents", 0)
, fvUsPadsFiredR()
, fvUsPadsFiredXY()
, fvTrackCharge()
@@ -77,13 +75,35 @@ CbmMuchDigitizerQa::~CbmMuchDigitizerQa() { DeInit(); }
// -------------------------------------------------------------------------
void CbmMuchDigitizerQa::DeInit() {
- for (int i = 0; i < fNstations; i++) {
- delete fvUsPadsFiredR[i];
- delete fvUsPadsFiredXY[i];
- delete fvTrackCharge[i];
- delete fvPadsTotalR[i];
- delete fvPadsFiredR[i];
- delete fvPadOccupancyR[i];
+
+ histFolder = nullptr;
+ fGeoScheme = nullptr;
+ fDigiManager = nullptr;
+ fPoints = nullptr;
+ fDigis = nullptr;
+ fDigiMatches = nullptr;
+ fMCTracks = nullptr;
+ fOutFolder.Clear();
+
+ SafeDelete(fPointInfos);
+
+ for (uint i = 0; i < fvUsPadsFiredR.size(); i++) {
+ SafeDelete(fvUsPadsFiredR[i]);
+ }
+ for (uint i = 0; i < fvUsPadsFiredXY.size(); i++) {
+ SafeDelete(fvUsPadsFiredXY[i]);
+ }
+ for (uint i = 0; i < fvTrackCharge.size(); i++) {
+ SafeDelete(fvTrackCharge[i]);
+ }
+ for (uint i = 0; i < fvPadsTotalR.size(); i++) {
+ SafeDelete(fvPadsTotalR[i]);
+ }
+ for (uint i = 0; i < fvPadsFiredR.size(); i++) {
+ SafeDelete(fvPadsFiredR[i]);
+ }
+ for (uint i = 0; i < fvPadOccupancyR.size(); i++) {
+ SafeDelete(fvPadOccupancyR[i]);
}
fvUsPadsFiredR.clear();
fvUsPadsFiredXY.clear();
@@ -92,13 +112,46 @@ void CbmMuchDigitizerQa::DeInit() {
fvPadsFiredR.clear();
fvPadOccupancyR.clear();
- fNstations = 0;
- fOutFolder.Clear();
+ SafeDelete(fhTrackLength);
+ SafeDelete(fhTrackLengthPi);
+ SafeDelete(fhTrackLengthPr);
+ SafeDelete(fhTrackLengthEl);
+ SafeDelete(fhTrackChargeVsTrackEnergyLog);
+ SafeDelete(fhTrackChargeVsTrackEnergyLogPi);
+ SafeDelete(fhTrackChargeVsTrackEnergyLogPr);
+ SafeDelete(fhTrackChargeVsTrackEnergyLogEl);
+ SafeDelete(fhTrackChargeVsTrackLength);
+ SafeDelete(fhTrackChargeVsTrackLengthPi);
+ SafeDelete(fhTrackChargeVsTrackLengthPr);
+ SafeDelete(fhTrackChargeVsTrackLengthEl);
+ SafeDelete(fhNpadsVsS);
+ SafeDelete(fCanvCharge);
+ SafeDelete(fCanvStationCharge);
+ SafeDelete(fCanvChargeVsEnergy);
+ SafeDelete(fCanvChargeVsLength);
+ SafeDelete(fCanvTrackLength);
+ SafeDelete(fCanvNpadsVsArea);
+ SafeDelete(fCanvUsPadsFiredXY);
+ SafeDelete(fCanvPadOccupancyR);
+ SafeDelete(fCanvPadsTotalR);
+ SafeDelete(fFitEl);
+ SafeDelete(fFitPi);
+ SafeDelete(fFitPr);
+
+ fNevents.SetVal(0);
+ fNstations = 0;
+ fnPadSizesX = 0;
+ fnPadSizesY = 0;
+ fPadMinLx = 0.;
+ fPadMinLy = 0.;
+ fPadMaxLx = 0.;
+ fPadMaxLy = 0.;
}
// -------------------------------------------------------------------------
InitStatus CbmMuchDigitizerQa::Init() {
-
+ DeInit();
+ fPointInfos = new TClonesArray("CbmMuchPointInfo", 10);
TDirectory* oldDirectory = gDirectory;
FairRootManager* fManager = FairRootManager::Instance();
@@ -134,8 +187,9 @@ InitStatus CbmMuchDigitizerQa::Init() {
fMCTracks = nullptr;
fPoints = nullptr;
}
-
histFolder = fOutFolder.AddFolder("hist", "Histogramms");
+ fNevents.SetVal(0);
+ histFolder->Add(&fNevents);
//fVerbose = 3;
InitCanvases();
@@ -570,9 +624,8 @@ void CbmMuchDigitizerQa::SetParContainers() {
// -------------------------------------------------------------------------x
void CbmMuchDigitizerQa::Exec(Option_t*) {
-
- fNevents++;
- LOG(info) << "Event: " << fNevents;
+ fNevents.SetVal(fNevents.GetVal() + 1);
+ LOG(debug) << "Event: " << fNevents.GetVal();
if (CheckConsistency() != 0) { return; }
@@ -910,7 +963,7 @@ void CbmMuchDigitizerQa::DrawPadCanvases() {
for (Int_t i = 0; i < fNstations; i++) {
*fvPadsFiredR[i] = *fvUsPadsFiredR[i];
//fvPadsFiredR[i]->Sumw2();
- fvPadsFiredR[i]->Scale(1. / fNevents);
+ fvPadsFiredR[i]->Scale(1. / fNevents.GetVal());
fvPadOccupancyR[i]->Divide(fvPadsFiredR[i], fvPadsTotalR[i]);
fvPadOccupancyR[i]->Scale(100.);
@@ -1075,4 +1128,4 @@ Double_t CbmMuchDigitizerQa::MPV_n_e(Double_t Tkin, Double_t mass) {
if (logT < min_logT_p) logT = min_logT_p;
return fPol6.EvalPar(&logT, mpv_p);
}
-}
\ No newline at end of file
+}
diff --git a/sim/detectors/much/qa/CbmMuchDigitizerQa.h b/sim/detectors/much/qa/CbmMuchDigitizerQa.h
index 5d180849e1..1ef0a78cd1 100644
--- a/sim/detectors/much/qa/CbmMuchDigitizerQa.h
+++ b/sim/detectors/much/qa/CbmMuchDigitizerQa.h
@@ -10,6 +10,7 @@
#define CbmMuchDigitizerQa_H
#include "FairTask.h"
+#include "TParameter.h"
#include <Rtypes.h>
#include <RtypesCore.h>
#include <TFolder.h>
@@ -96,8 +97,8 @@ private:
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
+ TFolder fOutFolder; /// output folder with histos and canvases
+ TParameter<int> fNevents; /// number of processed events
// internal unscaled histograms, need to be scaled at the output
std::vector<TH1F*> fvUsPadsFiredR; // fired pads vs R, per station
@@ -144,9 +145,8 @@ private:
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 fnPadSizesX = 0;
+ Int_t fnPadSizesY = 0;
Double_t fPadMinLx = 0.;
Double_t fPadMinLy = 0.;
diff --git a/sim/detectors/much/qa/CbmMuchTransportQa.cxx b/sim/detectors/much/qa/CbmMuchTransportQa.cxx
index 5556f27e3a..032bdd84f8 100644
--- a/sim/detectors/much/qa/CbmMuchTransportQa.cxx
+++ b/sim/detectors/much/qa/CbmMuchTransportQa.cxx
@@ -41,7 +41,7 @@ ClassImp(CbmMuchTransportQa);
CbmMuchTransportQa::CbmMuchTransportQa(const char* name, Int_t verbose)
: FairTask(name, verbose)
, fOutFolder("MuchTransportQA", "Much Transport QA")
- , fhNevents("nEvents", 0)
+ , fNevents("nEvents", 0)
, fvUsNtra()
, fvMcPointXY()
, fvMcPointPhiZ()
@@ -59,7 +59,8 @@ void CbmMuchTransportQa::DeInit() {
fPoints = nullptr;
fMcTracks = nullptr;
fOutFolder.Clear();
- fhNevents.SetVal(0);
+ histFolder = nullptr;
+ fNevents.SetVal(0);
SafeDelete(fhUsNtraAll);
SafeDelete(fhUsNtraPrim);
@@ -69,6 +70,7 @@ void CbmMuchTransportQa::DeInit() {
SafeDelete(fhUsNtraEl);
SafeDelete(fhUsNtraMu);
SafeDelete(fhUsNtraKa);
+ fvUsNtra.clear();
for (uint i = 0; i < fvMcPointXY.size(); i++) {
SafeDelete(fvMcPointXY[i]);
@@ -83,12 +85,6 @@ void CbmMuchTransportQa::DeInit() {
fvMcPointPhiZ.clear();
fvMcPointRZ.clear();
- for (uint i = 0; i < fvMcPointPRatio.size(); i++) {
- SafeDelete(fvMcPointPRatio[i]);
- }
- for (uint i = 0; i < fvMcPointPrimRatio.size(); i++) {
- SafeDelete(fvMcPointPrimRatio[i]);
- }
SafeDelete(fhNtracks);
SafeDelete(fhFractionPrim);
SafeDelete(fhFractionSec);
@@ -97,9 +93,14 @@ void CbmMuchTransportQa::DeInit() {
SafeDelete(fhFractionEl);
SafeDelete(fhFractionMu);
SafeDelete(fhFractionKa);
-
- fvUsNtra.clear();
fvFraction.clear();
+
+ for (uint i = 0; i < fvMcPointPRatio.size(); i++) {
+ SafeDelete(fvMcPointPRatio[i]);
+ }
+ for (uint i = 0; i < fvMcPointPrimRatio.size(); i++) {
+ SafeDelete(fvMcPointPrimRatio[i]);
+ }
fvMcPointPRatio.clear();
fvMcPointPrimRatio.clear();
@@ -109,13 +110,12 @@ void CbmMuchTransportQa::DeInit() {
SafeDelete(fCanvNtra);
SafeDelete(fCanvStationPRatio);
SafeDelete(fCanvStationPrimRatio);
-
fNstations = 0;
- fOutFolder.Clear();
}
// -------------------------------------------------------------------------
InitStatus CbmMuchTransportQa::Init() {
+ DeInit();
TDirectory* oldDirectory = gDirectory;
FairRootManager* manager = FairRootManager::Instance();
@@ -151,8 +151,8 @@ InitStatus CbmMuchTransportQa::Init() {
return kFATAL;
}
}
- fhNevents.SetVal(0);
- histFolder->Add(&fhNevents);
+ fNevents.SetVal(0);
+ histFolder->Add(&fNevents);
InitCountingHistos();
InitFractionHistos();
@@ -355,8 +355,8 @@ void CbmMuchTransportQa::SetParContainers() {
// -------------------------------------------------------------------------
void CbmMuchTransportQa::Exec(Option_t*) {
- LOG(info) << "Event: " << fhNevents.GetVal();
- fhNevents.SetVal(fhNevents.GetVal() + 1);
+ fNevents.SetVal(fNevents.GetVal() + 1);
+ LOG(debug) << "Event: " << fNevents.GetVal();
// bitmask tells which stations were crossed by mc track
std::vector<UInt_t> trackStaCross(fMcTracks->GetEntriesFast(), 0);
@@ -450,7 +450,7 @@ TFolder& CbmMuchTransportQa::GetQa() {
TDirectory* oldDirectory = gDirectory;
fhNtracks->Reset();
- fhNtracks->Add(fhUsNtraAll, 1. / fhNevents.GetVal());
+ fhNtracks->Add(fhUsNtraAll, 1. / fNevents.GetVal());
std::vector<Double_t> errors(fNstations, 0.);
fhUsNtraAll->SetError(errors.data());
@@ -497,7 +497,7 @@ void CbmMuchTransportQa::DrawCanvases() {
PrimRatioPieLeg->SetX2(.90);
}
- double scale = (fhNevents.GetVal() > 0) ? 1. / fhNevents.GetVal() : 0;
+ double scale = (fNevents.GetVal() > 0) ? 1. / fNevents.GetVal() : 0;
int i = 1;
fCanvNtra->cd(i++);
@@ -589,4 +589,4 @@ void CbmMuchTransportQa::Finish() {
}
FairSink* sink = FairRootManager::Instance()->GetSink();
sink->WriteObject(&GetQa(), nullptr);
-}
\ No newline at end of file
+}
diff --git a/sim/detectors/much/qa/CbmMuchTransportQa.h b/sim/detectors/much/qa/CbmMuchTransportQa.h
index bf1f6e78b5..73cc74e190 100644
--- a/sim/detectors/much/qa/CbmMuchTransportQa.h
+++ b/sim/detectors/much/qa/CbmMuchTransportQa.h
@@ -78,9 +78,9 @@ private:
TClonesArray* fMcTracks = nullptr;
///
- TFolder* histFolder; /// subfolder for histograms
- TFolder fOutFolder; /// output folder with histos and canvases
- TParameter<int> fhNevents; /// number of processed events
+ TFolder* histFolder; /// subfolder for histograms
+ TFolder fOutFolder; /// output folder with histos and canvases
+ TParameter<int> fNevents; /// number of processed events
/// internal unscaled histogramms
TH1F* fhUsNtraAll = nullptr; /// number of all tracks
--
GitLab