From 8ea81a6d00318b1359faa4e2c4e6988a9798a616 Mon Sep 17 00:00:00 2001
From: sgorbuno <se.gorbunov@gsi.de>
Date: Wed, 10 Mar 2021 13:53:36 +0000
Subject: [PATCH] make the QA tasks work with timeslices
---
core/detectors/much/CbmMuchPointInfo.h | 16 +-
macro/run/run_qa.C | 20 +-
reco/detectors/much/qa/CbmMuchHitFinderQa.cxx | 121 +--
reco/detectors/much/qa/CbmMuchHitFinderQa.h | 24 +-
reco/detectors/sts/CMakeLists.txt | 4 +
reco/detectors/sts/CbmStsFindTracksQa.cxx | 726 +++++++++---------
reco/detectors/sts/CbmStsFindTracksQa.h | 176 +++--
sim/detectors/much/qa/CbmMuchDigitizerQa.cxx | 234 +++---
sim/detectors/much/qa/CbmMuchDigitizerQa.h | 32 +-
sim/detectors/much/qa/CbmMuchTransportQa.cxx | 122 +--
sim/detectors/much/qa/CbmMuchTransportQa.h | 14 +-
11 files changed, 833 insertions(+), 656 deletions(-)
diff --git a/core/detectors/much/CbmMuchPointInfo.h b/core/detectors/much/CbmMuchPointInfo.h
index 0d38fa7009..83d2a23b5e 100644
--- a/core/detectors/much/CbmMuchPointInfo.h
+++ b/core/detectors/much/CbmMuchPointInfo.h
@@ -31,14 +31,14 @@ public:
}
void Print(Option_t* = "") const;
- Double_t GetKine() { return fKine; }
- Double_t GetLength() { return fLength; }
- Int_t GetPdgCode() { return fPdgCode; }
- Int_t GetMotherPdg() { return fMotherPdg; }
- Int_t GetCharge() { return fCharge; }
- Int_t GetStationId() { return fStationId; }
- Double_t GetArea() { return fS; }
- Int_t GetNPads() { return fNPads; }
+ Double_t GetKine() const { return fKine; }
+ Double_t GetLength() const { return fLength; }
+ Int_t GetPdgCode() const { return fPdgCode; }
+ Int_t GetMotherPdg() const { return fMotherPdg; }
+ Int_t GetCharge() const { return fCharge; }
+ Int_t GetStationId() const { return fStationId; }
+ Double_t GetArea() const { return fS; }
+ Int_t GetNPads() const { return fNPads; }
private:
Double_t fKine;
diff --git a/macro/run/run_qa.C b/macro/run/run_qa.C
index d0104352c3..86bfae43e8 100644
--- a/macro/run/run_qa.C
+++ b/macro/run/run_qa.C
@@ -33,9 +33,8 @@
#include <TStopwatch.h>
#endif
-void run_qa(Int_t nEvents = 0,
- TString dataset = "data/sis100_muon_jpsi_test",
- TString setup = "sis100_muon_jpsi") {
+void run_qa(TString dataset = "data/sis100_muon_jpsi_test", TString setup = "sis100_muon_jpsi", Int_t nEvents = -1)
+{
// ========================================================================
// Adjust this part according to your requirements
@@ -51,10 +50,10 @@ void run_qa(Int_t nEvents = 0,
// ------------------------------------------------------------------------
// ----- In- and output file names ------------------------------------
- TString rawFile = dataset + ".event.raw.root";
+ TString rawFile = dataset + ".raw.root";
TString traFile = dataset + ".tra.root";
TString parFile = dataset + ".par.root";
- TString recFile = dataset + ".rec.root";
+ TString recFile = dataset + ".reco.root";
TString sinkFile = dataset + ".qa.root";
// ------------------------------------------------------------------------
@@ -140,18 +139,25 @@ void run_qa(Int_t nEvents = 0,
FairMonitor::GetMonitor()->EnableMonitor(kTRUE, monitorFile);
// ------------------------------------------------------------------------
- // ----- MCDataManager (legacy mode) -----------------------------------
- CbmMCDataManager* mcManager = new CbmMCDataManager("MCDataManager", 1);
+ // ----- MCDataManager -----------------------------------
+ CbmMCDataManager* mcManager = new CbmMCDataManager("MCDataManager", 0);
mcManager->AddFile(traFile);
run->AddTask(mcManager);
// ------------------------------------------------------------------------
+ // ----- Match reco to MC ------
+ CbmMatchRecoToMC* matchTask = new CbmMatchRecoToMC();
+ run->AddTask(matchTask);
+
// ----- MUCH QA ---------------------------------
if (CbmSetup::Instance()->IsActive(ECbmModuleId::kMuch)) {
run->AddTask(new CbmMuchTransportQa());
run->AddTask(new CbmMuchDigitizerQa());
run->AddTask(new CbmMuchHitFinderQa());
}
+
+ // ----- STS QA ---------------------------------
+ if (CbmSetup::Instance()->IsActive(ECbmModuleId::kSts)) { run->AddTask(new CbmStsFindTracksQa()); }
// ------------------------------------------------------------------------
// ----- Parameter database --------------------------------------------
diff --git a/reco/detectors/much/qa/CbmMuchHitFinderQa.cxx b/reco/detectors/much/qa/CbmMuchHitFinderQa.cxx
index ebe7f84917..ca3e59c9fe 100644
--- a/reco/detectors/much/qa/CbmMuchHitFinderQa.cxx
+++ b/reco/detectors/much/qa/CbmMuchHitFinderQa.cxx
@@ -11,6 +11,8 @@
#include "CbmDefs.h"
#include "CbmDigiManager.h"
#include "CbmLink.h"
+#include "CbmMCDataArray.h"
+#include "CbmMCDataManager.h"
#include "CbmMCTrack.h"
#include "CbmMatch.h"
#include "CbmMuchAddress.h"
@@ -21,6 +23,7 @@
#include "CbmMuchPixelHit.h"
#include "CbmMuchPoint.h"
#include "CbmQaCanvas.h"
+#include "CbmTimeSlice.h"
#include "FairRootManager.h"
#include <FairSink.h>
@@ -46,9 +49,10 @@
#include <algorithm>
#include <iostream>
#include <map>
-#include <stdio.h>
#include <vector>
+#include <stdio.h>
+
using std::cout;
using std::endl;
using std::map;
@@ -129,10 +133,19 @@ void CbmMuchHitFinderQa::DeInit() {
// -------------------------------------------------------------------------
InitStatus CbmMuchHitFinderQa::Init() {
+
DeInit();
- FairRootManager* fManager = FairRootManager::Instance();
- fMCTracks = (TClonesArray*) fManager->GetObject("MCTrack");
- fPoints = (TClonesArray*) fManager->GetObject("MuchPoint");
+
+ fManager = FairRootManager::Instance();
+ fMcManager = dynamic_cast<CbmMCDataManager*>(fManager->GetObject("MCDataManager"));
+
+ fTimeSlice = (CbmTimeSlice*) fManager->GetObject("TimeSlice.");
+
+ if (fMcManager) {
+ fMCTracks = fMcManager->InitBranch("MCTrack");
+ fPoints = fMcManager->InitBranch("MuchPoint");
+ }
+
fHits = (TClonesArray*) fManager->GetObject("MuchPixelHit");
fClusters = (TClonesArray*) fManager->GetObject("MuchCluster");
// Reading Much Digis from CbmMuchDigiManager which are stored as vector
@@ -142,7 +155,7 @@ InitStatus CbmMuchHitFinderQa::Init() {
TFile* oldFile = gFile;
TDirectory* oldDir = gDirectory;
- TFile* f = new TFile(fGeoFileName, "R");
+ TFile* f = new TFile(fGeoFileName, "R");
/// Restore old global file and folder pointer to avoid messing with FairRoot
gFile = oldFile;
@@ -310,7 +323,8 @@ void CbmMuchHitFinderQa::Exec(Option_t*) {
// -------------------------------------------------------------------------
void CbmMuchHitFinderQa::FinishTask() {
- printf("FinishTask\n");
+
+ if (fVerbose > 0) { printf(" CbmMuchHitFinderQa FinishTask\n"); }
gStyle->SetPaperSize(20, 20);
@@ -320,11 +334,6 @@ void CbmMuchHitFinderQa::FinishTask() {
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);
@@ -337,18 +346,23 @@ void CbmMuchHitFinderQa::FinishTask() {
TH1D* histo = vResHistos[i];
histo->Sumw2();
histo->Fit("gaus", "Q");
- histo->GetFunction("gaus")->SetLineWidth(3);
- histo->GetFunction("gaus")->SetLineColor(kRed);
+ auto f = histo->GetFunction("gaus");
+ if (f) {
+ f->SetLineWidth(3);
+ f->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());
+ if (fVerbose > 0) {
+ 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();
@@ -505,7 +519,7 @@ void CbmMuchHitFinderQa::StatisticsQa() {
// -------------------------------------------------------------------------
void CbmMuchHitFinderQa::PullsQa() {
Bool_t verbose = (fVerbose > 2);
- // Filling residuals and pools for hits at the first layer
+ // Filling residuals and pulls 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
@@ -539,7 +553,7 @@ void CbmMuchHitFinderQa::PullsQa() {
// 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;
+ CbmLink link;
for (Int_t digiId = 0; digiId < cluster->GetNofDigis(); digiId++) {
Int_t index = cluster->GetDigi(digiId);
@@ -572,7 +586,7 @@ void CbmMuchHitFinderQa::PullsQa() {
point_unique = 0;
break;
}
- Int_t currentPointId = match->GetLink(0).GetIndex();
+ CbmLink currentLink = match->GetLink(0);
CbmMuchModuleGem* module =
(CbmMuchModuleGem*) fGeoScheme->GetModuleByDetId(digi->GetAddress());
@@ -581,11 +595,11 @@ void CbmMuchHitFinderQa::PullsQa() {
return;
}
if (digiId == 0) {
- pointId = currentPointId;
+ link = currentLink;
continue;
}
// Not unique if mcPoint references differ for different digis
- if (currentPointId != pointId) {
+ if (!(currentLink == link)) {
point_unique = 0;
break;
}
@@ -597,8 +611,8 @@ void CbmMuchHitFinderQa::PullsQa() {
continue;
}
- if (verbose) printf(" pointId=%i", pointId);
- CbmMuchPoint* point = (CbmMuchPoint*) fPoints->At(pointId);
+ if (verbose) { printf(" file %i event %i pointId %i", link.GetFile(), link.GetEntry(), link.GetIndex()); }
+ CbmMuchPoint* point = (CbmMuchPoint*) fPoints->Get(link);
Double_t xMC = 0.5 * (point->GetXIn() + point->GetXOut());
Double_t yMC = 0.5 * (point->GetYIn() + point->GetYOut());
@@ -614,11 +628,11 @@ void CbmMuchHitFinderQa::PullsQa() {
// cout<<" Rec Hit time "<<tRC<<endl;
if (dx < 1.e-10) {
- printf("Anomalously small dx\n");
+ LOG(error) << "Anomalously small dx";
continue;
}
if (dy < 1.e-10) {
- printf("Anomalously small dy\n");
+ LOG(error) << "Anomalously small dy";
continue;
}
fhPullX->Fill((xRC - xMC) / dx);
@@ -636,15 +650,26 @@ void CbmMuchHitFinderQa::PullsQa() {
// -------------------------------------------------------------------------
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);
+ Int_t nClusters = fClusters->GetEntriesFast();
+ vector<CbmLink> vPoints;
+
+ {
+ const CbmMatch& match = fTimeSlice->GetMatch();
+ for (int iLink = 0; iLink < match.GetNofLinks(); iLink++) {
+ CbmLink link = match.GetLink(iLink);
+ int nMuchPoints = fPoints->Size(link);
+ for (Int_t iPoint = 0; iPoint < nMuchPoints; iPoint++) {
+ link.SetIndex(iPoint);
+ link.SetWeight(0.);
+ if (IsSignalPoint(link)) { fSignalPoints.SetVal(fSignalPoints.GetVal() + 1); }
+ vPoints.push_back(link);
+ }
+ }
}
+ std::sort(vPoints.begin(), vPoints.end());
+
for (Int_t iCluster = 0; iCluster < nClusters; ++iCluster) {
CbmMuchCluster* cluster = (CbmMuchCluster*) fClusters->At(iCluster);
if (!cluster) {
@@ -663,29 +688,27 @@ void CbmMuchHitFinderQa::ClusterDeconvQa() {
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);
+ CbmLink pointLink = match->GetLink(ip);
+ auto it = find(vPoints.begin(), vPoints.end(), pointLink);
+ assert(it != vPoints.end());
+ if (it->GetWeight() > 0.) continue;
+ it->SetWeight(1.);
+ if (IsSignalPoint(pointLink)) { 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);
+Bool_t CbmMuchHitFinderQa::IsSignalPoint(CbmLink pointLink)
+{
+
+ CbmMuchPoint* point = (CbmMuchPoint*) fPoints->Get(pointLink);
if (!point) return kFALSE;
Int_t iTrack = point->GetTrackID();
- CbmMCTrack* track = (iTrack < 0 || iTrack > nTracks - 1)
- ? NULL
- : (CbmMCTrack*) fMCTracks->At(iTrack);
+ CbmMCTrack* track = (CbmMCTrack*) fMCTracks->Get(pointLink.GetFile(), pointLink.GetEntry(), 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
diff --git a/reco/detectors/much/qa/CbmMuchHitFinderQa.h b/reco/detectors/much/qa/CbmMuchHitFinderQa.h
index a61a438f11..1225734fc1 100644
--- a/reco/detectors/much/qa/CbmMuchHitFinderQa.h
+++ b/reco/detectors/much/qa/CbmMuchHitFinderQa.h
@@ -7,6 +7,8 @@
#ifndef CbmMuchHitFinderQa_H
#define CbmMuchHitFinderQa_H
+#include "CbmLink.h"
+
#include "FairTask.h"
#include "TParameter.h"
@@ -24,6 +26,9 @@ class TClonesArray;
class TH1D;
class TH1I;
class TMemberInspector;
+class CbmMCDataArray;
+class CbmMCDataManager;
+class CbmTimeSlice;
class CbmMuchHitFinderQa : public FairTask {
@@ -56,12 +61,20 @@ private:
void DeInit();
void DrawCanvases();
+ // setup
+ FairRootManager* fManager = nullptr; //!
+ CbmMCDataManager* fMcManager = nullptr; //!
+ CbmTimeSlice* fTimeSlice = nullptr; //!
+
+ //
CbmMuchGeoScheme* fGeoScheme;
TString fGeoFileName;
TString fFileName;
- Int_t fVerbose = 0;
- Int_t fFlag = 0;
- TClonesArray* fPoints = nullptr;
+ Int_t fVerbose = 0;
+ Int_t fFlag = 0;
+
+ CbmMCDataArray* fMCTracks = nullptr;
+ CbmMCDataArray* fPoints = nullptr;
CbmDigiManager* fDigiManager = nullptr;
TFolder fOutFolder; /// output folder with histos and canvases
@@ -69,7 +82,6 @@ private:
TClonesArray* fClusters = nullptr;
TClonesArray* fHits = nullptr;
- TClonesArray* fMCTracks = nullptr;
Int_t fNstations = 0;
@@ -101,12 +113,12 @@ private:
TParameter<int> fPointsOverCounted;
/** Defines whether the point with the given index is signal point. **/
- Bool_t IsSignalPoint(Int_t iPoint);
+ Bool_t IsSignalPoint(CbmLink pointLink);
CbmMuchHitFinderQa(const CbmMuchHitFinderQa&);
CbmMuchHitFinderQa& operator=(const CbmMuchHitFinderQa&);
- ClassDef(CbmMuchHitFinderQa, 1)
+ ClassDef(CbmMuchHitFinderQa, 0)
};
#endif
diff --git a/reco/detectors/sts/CMakeLists.txt b/reco/detectors/sts/CMakeLists.txt
index 732f085c0b..7a22abb28b 100644
--- a/reco/detectors/sts/CMakeLists.txt
+++ b/reco/detectors/sts/CMakeLists.txt
@@ -41,6 +41,10 @@ ${CBMDATA_DIR}
${CBMDATA_DIR}/base
${CBMDATA_DIR}/global
${CBMDATA_DIR}/sts
+${CBMDATA_DIR}/mvd
+
+# MVD
+${CBMROOT_SOURCE_DIR}/mvd
# STS
${CBMDETECTORBASE_DIR}/sts
diff --git a/reco/detectors/sts/CbmStsFindTracksQa.cxx b/reco/detectors/sts/CbmStsFindTracksQa.cxx
index 038969aeea..afba76850f 100644
--- a/reco/detectors/sts/CbmStsFindTracksQa.cxx
+++ b/reco/detectors/sts/CbmStsFindTracksQa.cxx
@@ -21,16 +21,21 @@
#include "FairRun.h"
// Includes from CbmRoot
-#include "CbmEvent.h"
#include "CbmMCDataArray.h"
#include "CbmMCDataManager.h"
#include "CbmMCTrack.h"
+#include "CbmMvdDetector.h"
+#include "CbmMvdHit.h"
+#include "CbmMvdPoint.h"
+#include "CbmMvdStationPar.h"
#include "CbmStsHit.h"
#include "CbmStsPoint.h"
#include "CbmStsSetup.h"
#include "CbmStsTrack.h"
+#include "CbmTimeSlice.h"
#include "CbmTrackMatchNew.h"
+#include "FairRunAna.h"
using std::fixed;
using std::right;
@@ -39,125 +44,18 @@ using std::setw;
// ----- Default constructor -------------------------------------------
-CbmStsFindTracksQa::CbmStsFindTracksQa(Int_t iVerbose)
- : FairTask("STSFindTracksQA", iVerbose)
- , fHitMap()
- , fMatchMap()
- , fQualiMap()
- , fEvents()
- , fMCTracks(NULL)
- , fStsPoints(NULL)
- , fStsHits(NULL)
- , fStsHitMatch(NULL)
- , fStsTracks(NULL)
- , fMatches(NULL)
- , fLegacy(kFALSE)
- , fTargetPos(0., 0., 0.)
- , fSetup(NULL)
- , fNStations(0)
- , fMinStations(3)
- , fQuota(0.7)
- , fhMomAccAll(new TH1F())
- , fhMomRecAll(new TH1F())
- , fhMomEffAll(new TH1F())
- , fhMomAccPrim(new TH1F())
- , fhMomRecPrim(new TH1F())
- , fhMomEffPrim(new TH1F())
- , fhMomAccSec(new TH1F())
- , fhMomRecSec(new TH1F())
- , fhMomEffSec(new TH1F())
- , fhNpAccAll(new TH1F())
- , fhNpRecAll(new TH1F())
- , fhNpEffAll(new TH1F())
- , fhNpAccPrim(new TH1F())
- , fhNpRecPrim(new TH1F())
- , fhNpEffPrim(new TH1F())
- , fhNpAccSec(new TH1F())
- , fhNpRecSec(new TH1F())
- , fhNpEffSec(new TH1F())
- , fhZAccSec(new TH1F())
- , fhZRecSec(new TH1F())
- , fhZEffSec(new TH1F())
- , fhNhClones(new TH1F())
- , fhNhGhosts(new TH1F())
- , fHistoList(new TList())
- , fNAccAll(0)
- , fNAccPrim(0)
- , fNAccRef(0)
- , fNAccSec(0)
- , fNRecAll(0)
- , fNRecPrim(0)
- , fNRecRef(0)
- , fNRecSec(0)
- , fNGhosts(0)
- , fNClones(0)
- , fNEvents(0)
- , fNEventsFailed(0)
- , fTime(0.)
- , fTimer() {}
+CbmStsFindTracksQa::CbmStsFindTracksQa(Int_t iVerbose) : FairTask("STSFindTracksQA", iVerbose) {}
// -------------------------------------------------------------------------
// ----- Standard constructor ------------------------------------------
-CbmStsFindTracksQa::CbmStsFindTracksQa(Int_t minStations,
- Double_t quota,
- Int_t iVerbose)
+CbmStsFindTracksQa::CbmStsFindTracksQa(Int_t minStations, Double_t quota, Int_t iVerbose)
: FairTask("STSFindTracksQA", iVerbose)
- , fHitMap()
- , fMatchMap()
- , fQualiMap()
- , fEvents(NULL)
- , fMCTracks(NULL)
- , fStsPoints(NULL)
- , fStsHits(NULL)
- , fStsHitMatch(NULL)
- , fStsTracks(NULL)
- , fMatches(NULL)
- , fLegacy(kFALSE)
- , fTargetPos(0., 0., 0.)
- , fSetup(NULL)
- , fNStations(0)
, fMinStations(minStations)
, fQuota(quota)
- , fhMomAccAll(new TH1F())
- , fhMomRecAll(new TH1F())
- , fhMomEffAll(new TH1F())
- , fhMomAccPrim(new TH1F())
- , fhMomRecPrim(new TH1F())
- , fhMomEffPrim(new TH1F())
- , fhMomAccSec(new TH1F())
- , fhMomRecSec(new TH1F())
- , fhMomEffSec(new TH1F())
- , fhNpAccAll(new TH1F())
- , fhNpRecAll(new TH1F())
- , fhNpEffAll(new TH1F())
- , fhNpAccPrim(new TH1F())
- , fhNpRecPrim(new TH1F())
- , fhNpEffPrim(new TH1F())
- , fhNpAccSec(new TH1F())
- , fhNpRecSec(new TH1F())
- , fhNpEffSec(new TH1F())
- , fhZAccSec(new TH1F())
- , fhZRecSec(new TH1F())
- , fhZEffSec(new TH1F())
- , fhNhClones(new TH1F())
- , fhNhGhosts(new TH1F())
- , fHistoList(new TList())
- , fNAccAll(0)
- , fNAccPrim(0)
- , fNAccRef(0)
- , fNAccSec(0)
- , fNRecAll(0)
- , fNRecPrim(0)
- , fNRecRef(0)
- , fNRecSec(0)
- , fNGhosts(0)
- , fNClones(0)
- , fNEvents(0)
- , fNEventsFailed(0)
- , fTime(0.)
- , fTimer() {}
+{
+}
// -------------------------------------------------------------------------
@@ -173,174 +71,82 @@ CbmStsFindTracksQa::~CbmStsFindTracksQa() {
// ----- Task execution ------------------------------------------------
void CbmStsFindTracksQa::Exec(Option_t* /*opt*/) {
- // If there is an event branch: do the event loop
- if (fEvents) {
- Int_t nEvents = fEvents->GetEntriesFast();
- LOG(debug) << GetName() << ": found time slice with " << nEvents
- << " events.";
-
- for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
- CbmEvent* event = dynamic_cast<CbmEvent*>(fEvents->At(iEvent));
- assert(event);
- ProcessEvent(event);
- }
- }
-
- // If there is no event branch, process the entire tree
- else {
- ProcessEvent();
- }
-}
-// -------------------------------------------------------------------------
-
-
-// ----- Public method SetParContainers --------------------------------
-void CbmStsFindTracksQa::SetParContainers() {}
-// -------------------------------------------------------------------------
-
-
-// ----- Public method Init --------------------------------------------
-InitStatus CbmStsFindTracksQa::Init() {
-
- LOG(info) << "\n\n====================================================";
- LOG(info) << GetName() << ": Initialising...";
-
- // Get RootManager
- FairRootManager* ioman = FairRootManager::Instance();
- assert(ioman);
-
- // Get STS setup
- fSetup = CbmStsSetup::Instance();
-
- // Get MCDataManager
- CbmMCDataManager* mcManager =
- dynamic_cast<CbmMCDataManager*>(ioman->GetObject("MCDataManager"));
- assert(mcManager);
-
- // Get MCTrack array
- fMCTracks = mcManager->InitBranch("MCTrack");
- assert(fMCTracks);
-
- // Get StsPoint array
- fStsPoints = mcManager->InitBranch("StsPoint");
- assert(fStsPoints);
-
- // Get Event array
- fEvents = dynamic_cast<TClonesArray*>(ioman->GetObject("CbmEvent"));
-
- // Get StsHit array
- fStsHits = (TClonesArray*) ioman->GetObject("StsHit");
- assert(fStsHits);
-
- // Get StsHitMatch array
- fStsHitMatch = (TClonesArray*) ioman->GetObject("StsHitMatch");
- assert(fStsHitMatch);
-
- // Get StsTrack array
- fStsTracks = (TClonesArray*) ioman->GetObject("StsTrack");
- assert(fStsTracks);
-
- // Get StsTrackMatch array
- fMatches = (TClonesArray*) ioman->GetObject("StsTrackMatch");
- assert(fMatches);
-
-
- // Get the geometry of target and STS
- InitStatus geoStatus = GetGeometry();
- if (geoStatus != kSUCCESS) {
- LOG(error) << GetName() << "::Init: Error in reading geometry!";
- return geoStatus;
- }
-
- // Create histograms
- CreateHistos();
- Reset();
-
- // Output
- LOG(info) << " Number of STS stations : " << fNStations;
- LOG(info) << " Target position ( " << fTargetPos.X() << ", "
- << fTargetPos.Y() << ", " << fTargetPos.Z() << ") cm";
- LOG(info) << " Minimum number of STS stations : " << fMinStations;
- LOG(info) << " Matching quota : " << fQuota;
- LOG(info) << "====================================================";
-
- return geoStatus;
-}
-// -------------------------------------------------------------------------
-
-
-// ----- Public method ReInit ------------------------------------------
-InitStatus CbmStsFindTracksQa::ReInit() {
-
- LOG(info) << "\n\n====================================================";
- LOG(info) << GetName() << ": Re-initialising...";
-
- // Get the geometry of target and STS
- InitStatus geoStatus = GetGeometry();
- if (geoStatus != kSUCCESS) {
- LOG(error) << GetName() << "::Init: Error in reading geometry!";
- return geoStatus;
- }
-
- // --- Screen log
- LOG(info) << " Number of STS stations : " << fNStations;
- LOG(info) << " Target position ( " << fTargetPos.X() << ", "
- << fTargetPos.Y() << ", " << fTargetPos.Z() << ") cm";
- LOG(info) << " Minimum number of STS stations : " << fMinStations;
- LOG(info) << " Matching quota : " << fQuota;
- LOG(info) << "====================================================";
-
- return geoStatus;
-}
-// -------------------------------------------------------------------------
-
-
-// ----- Public method Exec --------------------------------------------
-void CbmStsFindTracksQa::ProcessEvent(CbmEvent* event) {
-
- // --- Event number. Note that the FairRun counting start with 1.
- Int_t eventNumber =
- (event ? event->GetNumber()
- : FairRun::Instance()->GetEventHeader()->GetMCEntryNumber() - 1);
-
- LOG(debug) << GetName() << ": Process event " << eventNumber;
+ LOG(debug) << GetName() << ": Process event ";
// Timer
fTimer.Start();
// Eventwise counters
// Int_t nMCTracks = 0;
- Int_t nTracks = 0;
- Int_t nGhosts = 0;
- Int_t nClones = 0;
- Int_t nAll = 0;
- Int_t nAcc = 0;
- Int_t nRecAll = 0;
- Int_t nPrim = 0;
- Int_t nRecPrim = 0;
- Int_t nRef = 0;
- Int_t nRecRef = 0;
- Int_t nSec = 0;
- Int_t nRecSec = 0;
+ Int_t nTracks = 0;
+ Int_t nGhosts = 0;
+ Int_t nClones = 0;
+ Int_t nAll = 0;
+ Int_t nAcc = 0;
+ Int_t nRecAll = 0;
+ Int_t nPrim = 0;
+ Int_t nRecPrim = 0;
+ Int_t nRef = 0;
+ Int_t nRecRef = 0;
+ Int_t nRefLong = 0;
+ Int_t nRecRefLong = 0;
+ Int_t nSec = 0;
+ Int_t nRecSec = 0;
TVector3 momentum;
TVector3 vertex;
+ // check consistency
+ assert(fStsTracks->GetEntriesFast() == fStsTrackMatches->GetEntriesFast());
+
+ {
+ fMcTrackInfoMap.clear();
+ std::vector<CbmLink> events = fTimeSlice->GetMatch().GetLinks();
+ std::sort(events.begin(), events.end());
+ McTrackInfo info;
+ for (uint iLink = 0; iLink < events.size(); iLink++) {
+ CbmLink link = events[iLink];
+ Int_t nMCTracks = fMCTracks->Size(link);
+ for (Int_t iTr = 0; iTr < nMCTracks; iTr++) {
+ link.SetIndex(iTr);
+ fMcTrackInfoMap.insert({link, info});
+ }
+ }
+ }
+
// Fill hit and track maps
- FillHitMap(event);
- FillMatchMap(event, nTracks, nGhosts, nClones);
+ FillHitMap();
+ FillMatchMap(nTracks, nGhosts, nClones);
+
+ int nMcTracks = fMcTrackInfoMap.size();
// Loop over MCTracks
- Int_t nMcTracks = fMCTracks->Size(0, eventNumber);
- for (Int_t mcTrackId = 0; mcTrackId < nMcTracks; mcTrackId++) {
- CbmMCTrack* mcTrack =
- dynamic_cast<CbmMCTrack*>(fMCTracks->Get(0, eventNumber, mcTrackId));
+ int iMcTrack = 0;
+ for (auto itTrack = fMcTrackInfoMap.begin(); itTrack != fMcTrackInfoMap.end(); ++itTrack, ++iMcTrack) {
+ const CbmLink& link = itTrack->first;
+ McTrackInfo& info = itTrack->second;
+ CbmMCTrack* mcTrack = dynamic_cast<CbmMCTrack*>(fMCTracks->Get(link));
assert(mcTrack);
// Continue only for reconstructible tracks
nAll++;
- if (fHitMap.find(mcTrackId) == fHitMap.end()) continue; // No hits
- Int_t nStations = fHitMap[mcTrackId].size();
+ Int_t nStations = info.fHitMap.size();
if (nStations < fMinStations) continue; // Too few stations
+
+ int nContStations = 0; // Number of continious stations
+ {
+ int istaprev = -1;
+ int len = 0;
+ for (auto itSta = info.fHitMap.begin(); itSta != info.fHitMap.end(); itSta++) {
+ if (len == 0 || itSta->first == istaprev + 1) { len++; }
+ else {
+ len = 1;
+ }
+ if (nContStations < len) { nContStations = len; }
+ istaprev = itSta->first;
+ }
+ }
+ if (nContStations < fMinStations) continue; // Too few stations
+
nAcc++;
// Check origin of MCTrack
@@ -363,6 +169,12 @@ void CbmStsFindTracksQa::ProcessEvent(CbmEvent* event) {
nRef++;
}
+ Bool_t isRefLong = kFALSE;
+ if (isRef && nContStations >= fStsNstations) {
+ isRefLong = kTRUE;
+ nRefLong++;
+ }
+
// Fill histograms for reconstructible tracks
fhMomAccAll->Fill(mom);
fhNpAccAll->Fill(Double_t(nStations));
@@ -376,29 +188,26 @@ void CbmStsFindTracksQa::ProcessEvent(CbmEvent* event) {
}
// Get matched StsTrack
- Int_t trackId = -1;
- Double_t quali = 0.;
+ Int_t trackId = info.fStsTrackMatch;
+ Double_t quali = info.fQuali;
// Bool_t isRec = kFALSE;
- if (fMatchMap.find(mcTrackId) != fMatchMap.end()) {
- trackId = fMatchMap[mcTrackId];
+ if (trackId >= 0) {
// isRec = kTRUE;
CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(trackId);
assert(stsTrack);
- quali = fQualiMap[mcTrackId];
assert(quali >= fQuota);
- CbmTrackMatchNew* match = (CbmTrackMatchNew*) fMatches->At(trackId);
+ CbmTrackMatchNew* match = (CbmTrackMatchNew*) fStsTrackMatches->At(trackId);
assert(match);
Int_t nTrue = match->GetNofTrueHits();
Int_t nWrong = match->GetNofWrongHits();
//Int_t nFake = match->GetNofFakeHits();
Int_t nFake = 0;
- Int_t nAllHits = stsTrack->GetNofStsHits();
+ Int_t nAllHits = stsTrack->GetNofStsHits() + stsTrack->GetNofMvdHits();
+ if (!fIsMvdActive) { assert(stsTrack->GetNofMvdHits() == 0); }
assert(nTrue + nWrong + nFake == nAllHits);
-
// Verbose output
- LOG(debug1) << GetName() << ": MCTrack " << mcTrackId << ", stations "
- << nStations << ", hits " << nAllHits << ", true hits "
- << nTrue;
+ LOG(debug1) << GetName() << ": MCTrack " << iMcTrack << ", stations " << nStations << ", hits " << nAllHits
+ << ", true hits " << nTrue;
// Fill histograms for reconstructed tracks
nRecAll++;
@@ -408,24 +217,25 @@ void CbmStsFindTracksQa::ProcessEvent(CbmEvent* event) {
nRecPrim++;
fhMomRecPrim->Fill(mom);
fhNpRecPrim->Fill(Double_t(nAllHits));
- if (isRef) nRecRef++;
} else {
nRecSec++;
fhMomRecSec->Fill(mom);
fhNpRecSec->Fill(Double_t(nAllHits));
fhZRecSec->Fill(vertex.Z());
}
+ if (isRef) nRecRef++;
+ if (isRefLong) nRecRefLong++;
} // Match found in map?
} // Loop over MCTracks
-
// Calculate efficiencies
- Double_t effAll = Double_t(nRecAll) / Double_t(nAcc);
- Double_t effPrim = Double_t(nRecPrim) / Double_t(nPrim);
- Double_t effRef = Double_t(nRecRef) / Double_t(nRef);
- Double_t effSec = Double_t(nRecSec) / Double_t(nSec);
+ Double_t effAll = Double_t(nRecAll) / Double_t(nAcc);
+ Double_t effPrim = Double_t(nRecPrim) / Double_t(nPrim);
+ Double_t effRef = Double_t(nRecRef) / Double_t(nRef);
+ Double_t effRefLong = Double_t(nRecRefLong) / Double_t(nRefLong);
+ Double_t effSec = Double_t(nRecSec) / Double_t(nSec);
fTimer.Stop();
@@ -446,6 +256,8 @@ void CbmStsFindTracksQa::ProcessEvent(CbmEvent* event) {
LOG(debug) << "Reference : reconstructible: " << nRef
<< ", reconstructed: " << nRecRef << ", efficiency "
<< effRef * 100. << "%";
+ LOG(debug) << "Reference long : reconstructible: " << nRefLong << ", reconstructed: " << nRecRefLong
+ << ", efficiency " << effRefLong * 100. << "%";
LOG(debug) << "Non-vertex : reconstructible: " << nSec
<< ", reconstructed: " << nRecSec << ", efficiency "
<< effSec * 100. << "%";
@@ -457,13 +269,16 @@ void CbmStsFindTracksQa::ProcessEvent(CbmEvent* event) {
// Increase counters
+ fNAll += nAll;
fNAccAll += nAcc;
fNAccPrim += nPrim;
fNAccRef += nRef;
+ fNAccRefLong += nRefLong;
fNAccSec += nSec;
fNRecAll += nRecAll;
fNRecPrim += nRecPrim;
fNRecRef += nRecRef;
+ fNRecRefLong += nRecRefLong;
fNRecSec += nRecSec;
fNGhosts += nGhosts;
fNClones += nClones;
@@ -473,6 +288,136 @@ void CbmStsFindTracksQa::ProcessEvent(CbmEvent* event) {
// -------------------------------------------------------------------------
+// ----- Public method SetParContainers --------------------------------
+void CbmStsFindTracksQa::SetParContainers() {}
+// -------------------------------------------------------------------------
+
+
+// ----- Public method Init --------------------------------------------
+InitStatus CbmStsFindTracksQa::Init()
+{
+
+ LOG(info) << "\n\n====================================================";
+ LOG(info) << GetName() << ": Initialising...";
+
+ // Get STS setup
+ fStsSetup = CbmStsSetup::Instance();
+ assert(fStsSetup);
+ if (!fStsSetup->IsInit()) { fStsSetup->Init(); }
+
+ fManager = FairRootManager::Instance();
+ assert(fManager);
+
+ fMcManager = dynamic_cast<CbmMCDataManager*>(fManager->GetObject("MCDataManager"));
+
+ assert(fMcManager);
+
+ fTimeSlice = static_cast<CbmTimeSlice*>(fManager->GetObject("TimeSlice."));
+
+ if (fTimeSlice == nullptr) { LOG(fatal) << "CbmStsFindTracksQa: No time slice object"; }
+
+ if (fMcManager) {
+ fMCTracks = fMcManager->InitBranch("MCTrack");
+ fStsPoints = fMcManager->InitBranch("StsPoint");
+ }
+
+ assert(fMCTracks);
+ assert(fStsPoints);
+
+ // Get the geometry
+ InitStatus geoStatus = GetGeometry();
+ if (geoStatus != kSUCCESS) {
+ LOG(error) << GetName() << "::Init: Error in reading geometry!";
+ return geoStatus;
+ }
+
+ // MVD
+
+ fMvdPoints = fMcManager->InitBranch("MvdPoint");
+ fMvdCluster = (TClonesArray*) (fManager->GetObject("MvdCluster"));
+ fMvdHits = (TClonesArray*) (fManager->GetObject("MvdHit"));
+ fMvdHitMatch = (TClonesArray*) fManager->GetObject("MvdHitMatch");
+
+ // Currently in the time-based mode MVD is present but not reconstructed
+ // TODO: remove the check once the reconstruction works
+ if (fIsMvdActive && !fMvdHits) {
+ LOG(warning) << "CbmStsFindTracksQa: MVD hits are missing, MVD will not be "
+ "included to the STS track match";
+ fIsMvdActive = false;
+ }
+
+ if (fIsMvdActive) {
+ assert(fMvdPoints);
+ assert(fMvdCluster);
+ assert(fMvdHits);
+ assert(fMvdHitMatch);
+ }
+
+ // STS
+
+ // Get StsHit array
+ fStsHits = (TClonesArray*) fManager->GetObject("StsHit");
+ assert(fStsHits);
+
+ // Get StsHitMatch array
+ fStsHitMatch = (TClonesArray*) fManager->GetObject("StsHitMatch");
+ assert(fStsHitMatch);
+
+ // Get StsHitMatch array
+ fStsClusterMatch = (TClonesArray*) fManager->GetObject("StsClusterMatch");
+ assert(fStsClusterMatch);
+
+ // Get StsTrack array
+ fStsTracks = (TClonesArray*) fManager->GetObject("StsTrack");
+ assert(fStsTracks);
+
+ // Get StsTrackMatch array
+ fStsTrackMatches = (TClonesArray*) fManager->GetObject("StsTrackMatch");
+ assert(fStsTrackMatches);
+
+
+ // Create histograms
+ CreateHistos();
+ Reset();
+
+ // Output
+ LOG(info) << " Number of STS stations : " << fStsNstations;
+ LOG(info) << " Target position ( " << fTargetPos.X() << ", " << fTargetPos.Y() << ", " << fTargetPos.Z() << ") cm";
+ LOG(info) << " Minimum number of STS stations : " << fMinStations;
+ LOG(info) << " Matching quota : " << fQuota;
+ LOG(info) << "====================================================";
+
+ return geoStatus;
+}
+// -------------------------------------------------------------------------
+
+
+// ----- Public method ReInit ------------------------------------------
+InitStatus CbmStsFindTracksQa::ReInit()
+{
+
+ LOG(info) << "\n\n====================================================";
+ LOG(info) << GetName() << ": Re-initialising...";
+
+ // Get the geometry of target and STS
+ InitStatus geoStatus = GetGeometry();
+ if (geoStatus != kSUCCESS) {
+ LOG(error) << GetName() << "::Init: Error in reading geometry!";
+ return geoStatus;
+ }
+
+ // --- Screen log
+ LOG(info) << " Number of STS stations : " << fStsNstations;
+ LOG(info) << " Target position ( " << fTargetPos.X() << ", " << fTargetPos.Y() << ", " << fTargetPos.Z() << ") cm";
+ LOG(info) << " Minimum number of STS stations : " << fMinStations;
+ LOG(info) << " Matching quota : " << fQuota;
+ LOG(info) << "====================================================";
+
+ return geoStatus;
+}
+// -------------------------------------------------------------------------
+
+
// ----- Private method Finish -----------------------------------------
void CbmStsFindTracksQa::Finish() {
@@ -495,9 +440,10 @@ void CbmStsFindTracksQa::Finish() {
Double_t effAll = Double_t(fNRecAll) / Double_t(fNAccAll);
Double_t effPrim = Double_t(fNRecPrim) / Double_t(fNAccPrim);
Double_t effRef = Double_t(fNRecRef) / Double_t(fNAccRef);
+ Double_t effRefLong = Double_t(fNRecRefLong) / Double_t(fNAccRefLong);
Double_t effSec = Double_t(fNRecSec) / Double_t(fNAccSec);
- Double_t rateGhosts = Double_t(fNGhosts) / Double_t(fNEvents);
- Double_t rateClones = Double_t(fNClones) / Double_t(fNEvents);
+ Double_t rateGhosts = Double_t(fNGhosts) / Double_t(fNRecAll);
+ Double_t rateClones = Double_t(fNClones) / Double_t(fNRecAll);
// Run summary to screen
std::cout << std::endl;
@@ -510,34 +456,53 @@ void CbmStsFindTracksQa::Finish() {
<< fNRecPrim << "/" << fNAccPrim << ")";
LOG(info) << "Eff. reference tracks : " << effRef * 100 << " % (" << fNRecRef
<< "/" << fNAccRef << ")";
+ LOG(info) << "Eff. reference long tracks : " << effRefLong * 100 << " % (" << fNRecRefLong << "/" << fNAccRefLong
+ << ")";
LOG(info) << "Eff. secondary tracks : " << effSec * 100 << " % (" << fNRecSec
<< "/" << fNAccSec << ")";
- LOG(info) << "Ghost rate : " << rateGhosts << " per event";
- LOG(info) << "Clone rate : " << rateClones << " per event";
+ LOG(info) << "Ghost rate : " << rateGhosts * 100 << " % (" << fNGhosts << "/" << fNRecAll << ")";
+ LOG(info) << "Clone rate : " << rateClones * 100 << " % (" << fNClones << "/" << fNRecAll << ")";
+ LOG(info) << "mc tracks/event " << fNAll / fNEvents << " accepted " << fNRecAll / fNEvents;
LOG(info) << "Time per event : " << setprecision(6)
<< fTime / Double_t(fNEvents) << " s";
+
+ if (fMvdNstations > 0 && !fIsMvdActive) {
+ LOG(warning) << "CbmStsFindTracksQa: MVD hits are missing, MVD is not "
+ "included to the STS track match";
+ }
+
LOG(info) << "=====================================";
// Write histos to output
+ /*
gDirectory->mkdir("STSFindTracksQA");
gDirectory->cd("STSFindTracksQA");
TIter next(fHistoList);
while (TH1* histo = ((TH1*) next()))
histo->Write();
gDirectory->cd("..");
+*/
+ FairSink* sink = FairRootManager::Instance()->GetSink();
+ sink->WriteObject(&fOutFolder, nullptr);
}
// -------------------------------------------------------------------------
// ----- Private method GetGeometry ------------------------------------
InitStatus CbmStsFindTracksQa::GetGeometry() {
-
// Get target geometry
GetTargetPosition();
-
- fNStations = CbmStsSetup::Instance()->GetNofStations();
-
-
+ fMvdNstations = 0;
+ {
+ CbmMvdDetector* mvdDetector = CbmMvdDetector::Instance();
+ if (mvdDetector) {
+ CbmMvdStationPar* mvdStationPar = mvdDetector->GetParameterFile();
+ assert(mvdStationPar);
+ fMvdNstations = mvdStationPar->GetStationCount();
+ }
+ }
+ fIsMvdActive = (fMvdNstations > 0);
+ fStsNstations = CbmStsSetup::Instance()->GetNofStations();
return kSUCCESS;
}
// -------------------------------------------------------------------------
@@ -601,6 +566,8 @@ void CbmStsFindTracksQa::GetTargetPosition() {
// ----- Private method CreateHistos -----------------------------------
void CbmStsFindTracksQa::CreateHistos() {
+ fOutFolder.Clear();
+
// Histogram list
fHistoList = new TList();
@@ -692,6 +659,11 @@ void CbmStsFindTracksQa::CreateHistos() {
new TH1F("hNhGhosts", "number of hits for ghosts", nBinsNp, minNp, maxNp);
fHistoList->Add(fhNhClones);
fHistoList->Add(fhNhGhosts);
+
+ TIter next(fHistoList);
+ while (TH1* histo = ((TH1*) next())) {
+ fOutFolder.Add(histo);
+ }
}
// -------------------------------------------------------------------------
@@ -703,118 +675,144 @@ void CbmStsFindTracksQa::Reset() {
while (TH1* histo = ((TH1*) next()))
histo->Reset();
- fNAccAll = fNAccPrim = fNAccRef = fNAccSec = 0;
- fNRecAll = fNRecPrim = fNRecRef = fNRecSec = 0;
+ fNAccAll = fNAccPrim = fNAccRef = fNAccRefLong = fNAccSec = 0;
+ fNRecAll = fNRecPrim = fNRecRef = fNRecRefLong = fNRecSec = 0;
fNGhosts = fNClones = fNEvents = 0;
}
// -------------------------------------------------------------------------
// ----- Private method FillHitMap -------------------------------------
-void CbmStsFindTracksQa::FillHitMap(CbmEvent* event) {
-
- // --- Event number. Note that the FairRun counting starts with 1.
- Int_t eventNumber =
- (event ? event->GetNumber()
- : FairRun::Instance()->GetEventHeader()->GetMCEntryNumber() - 1);
+void CbmStsFindTracksQa::FillHitMap()
+{
// --- Fill hit map ( mcTrack -> ( station -> number of hits ) )
- fHitMap.clear();
- Int_t nHits = (event ? event->GetNofData(ECbmDataType::kStsHit)
- : fStsHits->GetEntriesFast());
- for (Int_t iHit = 0; iHit < nHits; iHit++) {
- Int_t hitIndex =
- (event ? event->GetIndex(ECbmDataType::kStsHit, iHit) : iHit);
- CbmStsHit* hit = (CbmStsHit*) fStsHits->At(hitIndex);
- CbmMatch* hitMatch = (CbmMatch*) fStsHitMatch->At(hitIndex);
- Int_t pointIndex = hitMatch->GetMatchedLink().GetIndex();
- assert(pointIndex >= 0);
- CbmStsPoint* stsPoint =
- dynamic_cast<CbmStsPoint*>(fStsPoints->Get(0, eventNumber, pointIndex));
- assert(stsPoint);
- Int_t mcTrackIndex = stsPoint->GetTrackID();
- Int_t station = fSetup->GetStationNumber(hit->GetAddress());
- fHitMap[mcTrackIndex][station]++;
+
+ // pocess MVD hits
+
+ if (fIsMvdActive) {
+ assert(fMvdHits);
+ assert(fMvdHitMatch);
+ assert(fMvdPoints);
+ for (Int_t iHit = 0; iHit < fMvdHits->GetEntriesFast(); iHit++) {
+ CbmMvdHit* hit = (CbmMvdHit*) fMvdHits->At(iHit);
+ assert(hit);
+ Int_t station = hit->GetStationNr();
+ assert(station >= 0 && station < fMvdNstations);
+ const CbmMatch* match = (const CbmMatch*) fMvdHitMatch->At(iHit);
+ assert(match);
+ if (match->GetNofLinks() <= 0) continue;
+ CbmLink link = match->GetMatchedLink();
+ CbmMvdPoint* pt = (CbmMvdPoint*) fMvdPoints->Get(link);
+ assert(pt);
+ link.SetIndex(pt->GetTrackID());
+ McTrackInfo& info = getMcTrackInfo(link);
+ info.fHitMap[station]++;
+ }
+ }
+
+ // pocess STS hits
+
+ for (Int_t iHit = 0; iHit < fStsHits->GetEntriesFast(); iHit++) {
+ CbmStsHit* hit = (CbmStsHit*) fStsHits->At(iHit);
+ assert(hit);
+
+ Int_t station = fStsSetup->GetStationNumber(hit->GetAddress());
+
+ // carefully check the hit match by looking at the strips from both sides
+
+ const CbmMatch* frontMatch = dynamic_cast<const CbmMatch*>(fStsClusterMatch->At(hit->GetFrontClusterId()));
+ assert(frontMatch);
+ if (frontMatch->GetNofLinks() <= 0) continue;
+
+ const CbmMatch* backMatch = dynamic_cast<const CbmMatch*>(fStsClusterMatch->At(hit->GetBackClusterId()));
+ assert(backMatch);
+ if (backMatch->GetNofLinks() <= 0) continue;
+
+ if (frontMatch->GetMatchedLink() == backMatch->GetMatchedLink()) {
+ CbmLink link = frontMatch->GetMatchedLink();
+ CbmStsPoint* pt = (CbmStsPoint*) fStsPoints->Get(link);
+ assert(pt);
+ link.SetIndex(pt->GetTrackID());
+ McTrackInfo& info = getMcTrackInfo(link);
+ info.fHitMap[fMvdNstations + station]++;
+ }
}
- LOG(debug) << GetName() << ": Filled hit map from " << nHits
- << " STS hits for " << fHitMap.size() << " MCTracks.";
+ LOG(debug) << GetName() << ": Filled hit map from " << fStsHits->GetEntriesFast() << " STS hits";
}
// -------------------------------------------------------------------------
// ------ Private method FillMatchMap ----------------------------------
-void CbmStsFindTracksQa::FillMatchMap(CbmEvent* event,
- Int_t& nRec,
- Int_t& nGhosts,
- Int_t& nClones) {
+void CbmStsFindTracksQa::FillMatchMap(Int_t& nRec, Int_t& nGhosts, Int_t& nClones)
+{
// Clear matching maps
- fMatchMap.clear();
- fQualiMap.clear();
+ for (auto it = fMcTrackInfoMap.begin(); it != fMcTrackInfoMap.end(); ++it) {
+ McTrackInfo& info = it->second;
+ info.fStsTrackMatch = -1;
+ info.fQuali = 0.;
+ info.fMatchedNHitsAll = 0;
+ info.fMatchedNHitsTrue = 0;
+ }
// Loop over StsTracks. Check matched MCtrack and fill maps.
- nGhosts = 0;
- nClones = 0;
- Int_t nTracks = (event ? event->GetNofData(ECbmDataType::kStsTrack)
- : fStsTracks->GetEntriesFast());
+ nGhosts = 0;
+ nClones = 0;
+ Int_t nTracks = fStsTracks->GetEntriesFast();
for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
// --- StsTrack
- Int_t trackIndex =
- (event ? event->GetIndex(ECbmDataType::kStsTrack, iTrack) : iTrack);
- CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(trackIndex);
+ CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(iTrack);
assert(stsTrack);
Int_t nHits = stsTrack->GetNofStsHits();
// --- TrackMatch
- CbmTrackMatchNew* match = (CbmTrackMatchNew*) fMatches->At(trackIndex);
+
+ assert(iTrack >= 0 && iTrack < fStsTrackMatches->GetEntriesFast());
+ CbmTrackMatchNew* match = (CbmTrackMatchNew*) fStsTrackMatches->At(iTrack);
assert(match);
Int_t nTrue = match->GetNofTrueHits();
- // --- Matched MCTrack
- Int_t mcTrackId = -1;
- if (nTrue > 0) mcTrackId = match->GetMatchedLink().GetIndex();
- if (mcTrackId < 0) {
+ // Check matching criterion (quota)
+ Double_t quali = Double_t(nTrue) / Double_t(nHits);
+
+ // Quality isn't good, it's a ghost
+
+ if (quali < fQuota) {
fhNhGhosts->Fill(nHits);
nGhosts++;
continue;
}
- // Check matching criterion (quota)
- Double_t quali = Double_t(nTrue) / Double_t(nHits);
- if (quali >= fQuota) {
-
- // No previous match for this MCTrack
- if (fMatchMap.find(mcTrackId) == fMatchMap.end()) {
- fMatchMap[mcTrackId] = iTrack;
- fQualiMap[mcTrackId] = quali;
- } //? no previous match
-
- // Previous match; take the better one
- else {
- if (fQualiMap[mcTrackId] < quali) {
- CbmStsTrack* oldTrack =
- (CbmStsTrack*) fStsTracks->At(fMatchMap[mcTrackId]);
- fhNhClones->Fill(Double_t(oldTrack->GetNofStsHits()));
- fMatchMap[mcTrackId] = iTrack;
- fQualiMap[mcTrackId] = quali;
- } //? new track matches better to MCTrack
-
- else
- fhNhClones->Fill(nHits);
- nClones++;
- } //? previous match found
-
- } //? true match ratio > quota
-
- // If not matched, it's a ghost
- else {
- fhNhGhosts->Fill(nHits);
- nGhosts++;
+ // Quality is good
+
+ // --- Matched MCTrack
+ assert(match->GetNofLinks() > 0);
+ const CbmLink& link = match->GetMatchedLink();
+ assert(link.GetIndex() >= 0);
+ McTrackInfo& info = getMcTrackInfo(link);
+
+ // previous match is better, this track is a clone
+ if ((quali < info.fQuali) || ((quali == info.fQuali) && (nTrue < info.fMatchedNHitsTrue))) {
+ fhNhClones->Fill(nHits);
+ nClones++;
+ continue;
}
+ // this track is better than the old one
+ if (info.fMatchedNHitsAll > 0) {
+ fhNhClones->Fill(info.fMatchedNHitsAll);
+ nClones++;
+ }
+ info.fStsTrackMatch = iTrack;
+ info.fQuali = quali;
+ info.fMatchedNHitsAll = nHits;
+ info.fMatchedNHitsTrue = nTrue;
+
} // Loop over StsTracks
+
nRec = nTracks;
LOG(debug) << GetName() << ": Filled match map for " << nRec
<< " STS tracks. Ghosts " << nGhosts << " Clones " << nClones;
diff --git a/reco/detectors/sts/CbmStsFindTracksQa.h b/reco/detectors/sts/CbmStsFindTracksQa.h
index db97586766..ea191d28db 100644
--- a/reco/detectors/sts/CbmStsFindTracksQa.h
+++ b/reco/detectors/sts/CbmStsFindTracksQa.h
@@ -10,19 +10,21 @@
** Quality check task for CbmStsFindTracks
**/
-
#ifndef CBMSTSTRACKFINDERQA_H
#define CBMSTSTRACKFINDERQA_H 1
+#include "CbmLink.h"
+
#include "FairTask.h"
+
#include "TStopwatch.h"
#include "TVector3.h"
class TH1F;
-class CbmEvent;
class CbmMCDataArray;
class CbmStsSetup;
-
+class CbmMCDataManager;
+class CbmTimeSlice;
class CbmStsFindTracksQa : public FairTask {
@@ -30,7 +32,6 @@ public:
/** Default constructor **/
CbmStsFindTracksQa(Int_t iVerbose = 1);
-
/** Standard constructor
*@param minHits Minimal number of StsHits for considered MCTracks
*@param quota True/all hits for track to be considered reconstructed
@@ -38,27 +39,21 @@ public:
**/
CbmStsFindTracksQa(Int_t minHits, Double_t quota, Int_t iVerbose);
-
/** Destructor **/
virtual ~CbmStsFindTracksQa();
-
/** Set parameter containers **/
virtual void SetParContainers();
-
/** Initialisation **/
virtual InitStatus Init();
-
/** Reinitialisation **/
virtual InitStatus ReInit();
-
/** Execution **/
virtual void Exec(Option_t* opt);
-
private:
/** Finish **/
virtual void Finish();
@@ -66,31 +61,27 @@ private:
/** Read the geometry parameters **/
InitStatus GetGeometry();
-
/** Get the target node from the geometry **/
void GetTargetPosition();
-
/** Create histograms **/
void CreateHistos();
-
/** Reset histograms and counters **/
void Reset();
+ // collect id's of MC events participating to the data
+ void CollectMcEvents(TClonesArray* Matches);
/** Fill a map from MCTrack index to number of corresponding StsHits **/
- void FillHitMap(CbmEvent* event);
-
+ void FillHitMap();
/** Fill a map from MCTrack index to matched StsTrack index
*@param nRec Number of reconstructed tracks (return)
*@param nGhosts Number of ghost tracks (return)
*@param nClones Number of clone tracks (return)
**/
- void
- FillMatchMap(CbmEvent* event, Int_t& nRec, Int_t& nGhosts, Int_t& nClones);
-
+ void FillMatchMap(Int_t& nRec, Int_t& nGhosts, Int_t& nClones);
/** Divide histograms (reco/all) with correct error for the efficiency
*@param histo1 reconstructed tracks
@@ -99,83 +90,130 @@ private:
**/
void DivideHistos(TH1* histo1, TH1* histo2, TH1* histo3);
+ struct McTrackInfo {
+ std::map<Int_t, Int_t> fHitMap = {}; /// Mvd / Sts station -> number of attached hits
+ Int_t fStsTrackMatch = -1; /// matched StsTrack index
+ Double_t fQuali = 0.; /// percentage of matched hits
+ Int_t fMatchedNHitsAll = 0; /// number of matched hits
+ Int_t fMatchedNHitsTrue = 0; /// number of matched hits
+ };
+
+ McTrackInfo& getMcTrackInfo(const CbmLink& link)
+ {
+ assert(fMcTrackInfoMap.find(link) != fMcTrackInfoMap.end());
+ return fMcTrackInfoMap[link];
+ }
+
+ std::map<CbmLink, McTrackInfo> fMcTrackInfoMap = {}; //! map track link -> track info
+
+ /// Setup
+ FairRootManager* fManager = nullptr; //!
+ CbmMCDataManager* fMcManager = nullptr; //!
+ CbmTimeSlice* fTimeSlice = nullptr; //!
+
+ /// MC tracks
+ CbmMCDataArray* fMCTracks = nullptr; //! MCtrack
+
+ /// Mvd
+ Bool_t fIsMvdActive = kTRUE; // is the Mvd module active
+ Int_t fMvdNstations = 0;
+ CbmMCDataArray* fMvdPoints = nullptr; //!
+ TClonesArray* fMvdCluster = nullptr; //!
+ TClonesArray* fMvdHits = nullptr; //!
+ TClonesArray* fMvdHitMatch = nullptr; //!
+
+ /// STS
+ CbmStsSetup* fStsSetup = nullptr; // STS setup interface
+ Int_t fStsNstations = 0; // Number of STS stations
+ CbmMCDataArray* fStsPoints = nullptr; //! StsPoints
+ TClonesArray* fStsHits = nullptr; //! StsHits
+ TClonesArray* fStsHitMatch = nullptr; //! StsHitMatch
+ TClonesArray* fStsClusterMatch = nullptr; //!StsClusterMatch
+ TClonesArray* fStsTracks = nullptr; //! StsTrack
+ TClonesArray* fStsTrackMatches = nullptr; //! StsTrackMatch
- /** Process one event
- ** @param event Pointer to event object
- **
- ** If a NULL pointer is given, the entire input array is processed
- ** (legacy mode)
- */
- void ProcessEvent(CbmEvent* event = NULL);
-
-
- /** Map from MCTrack index to number of attached StsHits in each station **/
- std::map<Int_t, std::map<Int_t, Int_t>> fHitMap;
+ /** Geometry parameters **/
+ TVector3 fTargetPos = {0., 0., 0.}; // Target centre position
+ /** Task parameters **/
- /** Map from MCTrack index to matched StsTrack index **/
- std::map<Int_t, Int_t> fMatchMap;
+ Int_t fMinStations = 4; // Minimal number of stations with hits for considered MCTrack
+ Double_t fQuota = 0.7; // True/all hits for track to be considered reconstructed
- /** Map from MCTrack index to percentage of matched hits **/
- std::map<Int_t, Double_t> fQualiMap;
+ TFolder fOutFolder = {"CbmStsFindTracksQa", "CbmStsFindTracksQa"}; /// output folder with histos and canvases
+ /** Histograms **/
- /** Pointers to data arrays **/
- TClonesArray* fEvents; //! Event
- CbmMCDataArray* fMCTracks; //! MCtrack
- CbmMCDataArray* fStsPoints; //! StsPoints
- TClonesArray* fStsHits; //! StsHits
- TClonesArray* fStsHitMatch; //! StsHitMatch
- TClonesArray* fStsTracks; //! StsTrack
- TClonesArray* fMatches; //! StsTrackMatch
+ // eff. vs. p, all
+ TH1F* fhMomAccAll = nullptr;
+ TH1F* fhMomRecAll = nullptr;
+ TH1F* fhMomEffAll = nullptr;
- /** Flag for legacy mode (event-by-event w/o event objects) **/
- Bool_t fLegacy;
+ // eff. vs. p, vertex
+ TH1F* fhMomAccPrim = nullptr;
+ TH1F* fhMomRecPrim = nullptr;
+ TH1F* fhMomEffPrim = nullptr;
+ // eff. vs. p, non-vertex
+ TH1F* fhMomAccSec = nullptr;
+ TH1F* fhMomRecSec = nullptr;
+ TH1F* fhMomEffSec = nullptr;
- /** Geometry parameters **/
- TVector3 fTargetPos; // Target centre position
- CbmStsSetup* fSetup; // STS setup interface
- Int_t fNStations; // Number of STS stations
+ // eff. vs. np, all
+ TH1F* fhNpAccAll = nullptr;
+ TH1F* fhNpRecAll = nullptr;
+ TH1F* fhNpEffAll = nullptr;
+ // eff. vs. np, vertex
+ TH1F* fhNpAccPrim = nullptr;
+ TH1F* fhNpRecPrim = nullptr;
+ TH1F* fhNpEffPrim = nullptr;
- /** Task parameters **/
- Int_t
- fMinStations; // Minimal number of stations with hits for considered MCTrack
- Double_t fQuota; // True/all hits for track to be considered reconstructed
+ // eff. vs. np, non-vertex
+ TH1F* fhNpAccSec = nullptr;
+ TH1F* fhNpRecSec = nullptr;
+ TH1F* fhNpEffSec = nullptr;
+ // eff. vs. z, non-vertex
+ TH1F* fhZAccSec = nullptr;
+ TH1F* fhZRecSec = nullptr;
+ TH1F* fhZEffSec = nullptr;
- /** Histograms **/
- TH1F *fhMomAccAll, *fhMomRecAll, *fhMomEffAll; // eff. vs. p, all
- TH1F *fhMomAccPrim, *fhMomRecPrim, *fhMomEffPrim; // eff. vs. p, vertex
- TH1F *fhMomAccSec, *fhMomRecSec, *fhMomEffSec; // eff. vs. p, non-vertex
- TH1F *fhNpAccAll, *fhNpRecAll, *fhNpEffAll; // eff. vs. np, all
- TH1F *fhNpAccPrim, *fhNpRecPrim, *fhNpEffPrim; // eff. vs. np, vertex
- TH1F *fhNpAccSec, *fhNpRecSec, *fhNpEffSec; // eff. vs. np, non-vertex
- TH1F *fhZAccSec, *fhZRecSec, *fhZEffSec; // eff. vs. z, non-vertex
- TH1F *fhNhClones, *fhNhGhosts; // # hits of clones and ghosts
+ // # hits of clones and ghosts
+ TH1F* fhNhClones = nullptr;
+ TH1F* fhNhGhosts = nullptr;
/** List of histograms **/
- TList* fHistoList;
+ TList* fHistoList = nullptr;
/** Counters **/
- Int_t fNAccAll, fNAccPrim, fNAccRef, fNAccSec;
- Int_t fNRecAll, fNRecPrim, fNRecRef, fNRecSec;
- Int_t fNGhosts, fNClones;
- Int_t fNEvents; /** Number of events with success **/
- Int_t fNEventsFailed; /** Number of events with failure **/
- Double_t fTime; /** Total real time used for good events **/
+ Int_t fNAll = 0;
+ Int_t fNAccAll = 0;
+ Int_t fNAccPrim = 0;
+ Int_t fNAccRef = 0;
+ Int_t fNAccRefLong = 0;
+ Int_t fNAccSec = 0;
+ Int_t fNRecAll = 0;
+ Int_t fNRecPrim = 0;
+ Int_t fNRecRef = 0;
+ Int_t fNRecRefLong = 0;
+ Int_t fNRecSec = 0;
+ Int_t fNGhosts = 0;
+ Int_t fNClones = 0;
+ Int_t fNEvents = 0; /** Number of events with success **/
+ Int_t fNEventsFailed = 0; /** Number of events with failure **/
+ Double_t fTime = 0.; /** Total real time used for good events **/
/** Timer **/
- TStopwatch fTimer;
+ TStopwatch fTimer = {};
CbmStsFindTracksQa(const CbmStsFindTracksQa&);
CbmStsFindTracksQa operator=(const CbmStsFindTracksQa&);
- ClassDef(CbmStsFindTracksQa, 2);
+ ClassDef(CbmStsFindTracksQa, 0);
};
diff --git a/sim/detectors/much/qa/CbmMuchDigitizerQa.cxx b/sim/detectors/much/qa/CbmMuchDigitizerQa.cxx
index 967038b0c4..ad0e463620 100644
--- a/sim/detectors/much/qa/CbmMuchDigitizerQa.cxx
+++ b/sim/detectors/much/qa/CbmMuchDigitizerQa.cxx
@@ -8,8 +8,11 @@
/// \date 21.10.2020
#include "CbmMuchDigitizerQa.h"
+
#include "CbmDigiManager.h"
#include "CbmLink.h"
+#include "CbmMCDataArray.h"
+#include "CbmMCDataManager.h"
#include "CbmMCTrack.h"
#include "CbmMatch.h"
#include "CbmMuchAddress.h"
@@ -23,6 +26,13 @@
#include "CbmMuchRecoDefs.h"
#include "CbmMuchStation.h"
#include "CbmQaCanvas.h"
+#include "CbmTimeSlice.h"
+
+#include <FairRootManager.h>
+#include <FairSink.h>
+#include <FairTask.h>
+#include <Logger.h>
+
#include "TClonesArray.h"
#include "TDatabasePDG.h"
#include "TF1.h"
@@ -32,21 +42,19 @@
#include "TParameter.h"
#include "TString.h"
#include "TStyle.h"
-#include <FairRootManager.h>
-#include <FairSink.h>
-#include <FairTask.h>
-#include <Logger.h>
#include <TAxis.h>
#include <TCanvas.h>
#include <TDirectory.h>
#include <TMath.h>
#include <TParticlePDG.h>
#include <TVector2.h>
+
#include <iostream>
#include <limits>
-#include <math.h>
#include <vector>
+#include <math.h>
+
using std::endl;
using std::vector;
@@ -84,8 +92,7 @@ void CbmMuchDigitizerQa::DeInit() {
fDigiMatches = nullptr;
fMCTracks = nullptr;
fOutFolder.Clear();
-
- SafeDelete(fPointInfos);
+ fMcPointInfoMap.clear();
for (uint i = 0; i < fvUsPadsFiredR.size(); i++) {
SafeDelete(fvUsPadsFiredR[i]);
@@ -151,10 +158,10 @@ void CbmMuchDigitizerQa::DeInit() {
// -------------------------------------------------------------------------
InitStatus CbmMuchDigitizerQa::Init() {
DeInit();
- fPointInfos = new TClonesArray("CbmMuchPointInfo", 10);
+ fMcPointInfoMap.clear();
TDirectory* oldDirectory = gDirectory;
- FairRootManager* fManager = FairRootManager::Instance();
+ fManager = FairRootManager::Instance();
if (!fManager) {
LOG(fatal) << "Can not find FairRootManager";
return kFATAL;
@@ -176,8 +183,21 @@ InitStatus CbmMuchDigitizerQa::Init() {
}
fDigiManager->Init();
- fMCTracks = (TClonesArray*) fManager->GetObject("MCTrack");
- fPoints = (TClonesArray*) fManager->GetObject("MuchPoint");
+ fMCTracks = nullptr;
+ fPoints = nullptr;
+
+ fMcManager = dynamic_cast<CbmMCDataManager*>(fManager->GetObject("MCDataManager"));
+
+ fTimeSlice = (CbmTimeSlice*) fManager->GetObject("TimeSlice.");
+ if (!fTimeSlice) { LOG(error) << GetName() << ": No time slice found"; }
+
+ if (fMcManager) {
+ // Get MCTrack array
+ fMCTracks = fMcManager->InitBranch("MCTrack");
+
+ // Get StsPoint array
+ fPoints = fMcManager->InitBranch("MuchPoint");
+ }
if (fMCTracks && fPoints
&& fDigiManager->IsMatchPresent(ECbmModuleId::kMuch)) {
@@ -635,12 +655,14 @@ void CbmMuchDigitizerQa::Exec(Option_t*) {
ProcessMCPoints();
FillChargePerPoint();
+
FillDigitizerPerformancePlots();
if (fVerbose > 1) {
PrintFrontLayerPoints();
PrintFrontLayerDigis();
}
+
gDirectory = oldDirectory;
return;
}
@@ -740,7 +762,7 @@ void CbmMuchDigitizerQa::OccupancyQa() {
// -------------------------------------------------------------------------
int CbmMuchDigitizerQa::ProcessMCPoints() {
- if (!fMCTracks || !fPoints) {
+ if (!fMCTracks || !fPoints || !fTimeSlice) {
LOG(debug)
<< " CbmMuchDigitizerQa::DigitizerQa(): Found no MC data. Skipping.";
return 0;
@@ -748,68 +770,75 @@ int CbmMuchDigitizerQa::ProcessMCPoints() {
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);
- if (!point) {
- LOG(error) << " Much MC point " << i << " out of "
- << fPoints->GetEntriesFast() << " doesn't exist";
- return -1;
- }
+ fMcPointInfoMap.clear();
- Int_t stId = CbmMuchAddress::GetStationIndex(point->GetDetectorID());
+ std::vector<CbmLink> events = fTimeSlice->GetMatch().GetLinks();
+ std::sort(events.begin(), events.end());
- if (stId < 0 || stId > fNstations) {
- LOG(error) << " Much MC point station " << stId << " is out of the range";
- return -1;
- }
+ for (uint iLink = 0; iLink < events.size(); iLink++) {
+ CbmLink link = events[iLink];
+ int nMuchPoints = fPoints->Size(link);
+ int nMcTracks = fMCTracks->Size(link);
+ for (Int_t ip = 0; ip < nMuchPoints; ip++) {
+ link.SetIndex(ip);
+ CbmMuchPoint* point = (CbmMuchPoint*) fPoints->Get(link);
+ if (!point) {
+ LOG(error) << " Much MC point " << ip << " out of " << nMuchPoints << " doesn't exist";
+ return -1;
+ }
- // Check if the point corresponds to a certain MC Track
- Int_t trackId = point->GetTrackID();
- if (trackId < 0 || trackId >= fMCTracks->GetEntriesFast()) {
- LOG(error) << " Much MC point track Id " << trackId
- << " is out of the range";
- return -1;
- }
+ Int_t stId = CbmMuchAddress::GetStationIndex(point->GetDetectorID());
- CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(trackId);
- if (!mcTrack) {
- LOG(error) << " MC track" << trackId << " is not found";
- return -1;
- }
+ if (stId < 0 || stId > fNstations) {
+ LOG(error) << " Much MC point station " << stId << " is out of the range";
+ return -1;
+ }
- Int_t pdgCode = mcTrack->GetPdgCode();
- TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pdgCode);
-
- point->PositionIn(vIn);
- point->PositionOut(vOut);
- point->Momentum(p);
- Double_t length = (vOut - vIn).Mag(); // Track length
- Double_t kine = 0; // Kinetic energy
- if (particle) {
- Double_t mass = particle->Mass();
- kine = sqrt(p.Mag2() + mass * mass) - mass;
- }
+ // Check if the point corresponds to a certain MC Track
+ Int_t trackId = point->GetTrackID();
+ if (trackId < 0 || trackId >= nMcTracks) {
+ LOG(error) << " Much MC point track Id " << trackId << " is out of the range";
+ return -1;
+ }
- // Get mother pdg code
- Int_t motherPdg = 0;
- Int_t motherId = mcTrack->GetMotherId();
- if (motherId >= fMCTracks->GetEntriesFast()) {
- LOG(error) << " MC track mother Id" << trackId << " is out of the range";
- return -1;
- }
- if (motherId >= 0) {
- CbmMCTrack* motherTrack = (CbmMCTrack*) fMCTracks->At(motherId);
- if (!motherTrack) {
+ CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->Get(link.GetFile(), link.GetEntry(), trackId);
+ if (!mcTrack) {
LOG(error) << " MC track" << trackId << " is not found";
return -1;
}
- motherPdg = motherTrack->GetPdgCode();
- }
- new ((*fPointInfos)[i])
- CbmMuchPointInfo(pdgCode, motherPdg, kine, length, stId);
- }
+
+ Int_t pdgCode = mcTrack->GetPdgCode();
+ TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pdgCode);
+
+ point->PositionIn(vIn);
+ point->PositionOut(vOut);
+ point->Momentum(p);
+ Double_t length = (vOut - vIn).Mag(); // Track length
+ Double_t kine = 0; // Kinetic energy
+ if (particle) {
+ Double_t mass = particle->Mass();
+ kine = sqrt(p.Mag2() + mass * mass) - mass;
+ }
+
+ // Get mother pdg code
+ Int_t motherPdg = 0;
+ Int_t motherId = mcTrack->GetMotherId();
+ if (motherId >= nMcTracks) {
+ LOG(error) << " MC track mother Id" << trackId << " is out of the range";
+ return -1;
+ }
+ if (motherId >= 0) {
+ CbmMCTrack* motherTrack = (CbmMCTrack*) fMCTracks->Get(link.GetFile(), link.GetEntry(), motherId);
+ if (!motherTrack) {
+ LOG(error) << " MC track" << trackId << " is not found";
+ return -1;
+ }
+ motherPdg = motherTrack->GetPdgCode();
+ }
+ fMcPointInfoMap.insert({link, CbmMuchPointInfo(pdgCode, motherPdg, kine, length, stId)});
+ } // points
+ } // events
return 0;
}
@@ -822,18 +851,21 @@ void CbmMuchDigitizerQa::FillChargePerPoint() {
return;
}
for (Int_t i = 0; i < fDigiManager->GetNofDigis(ECbmModuleId::kMuch); i++) {
- CbmMuchDigi* digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(i);
- CbmMatch* match =
- (CbmMatch*) fDigiManager->GetMatch(ECbmModuleId::kMuch, i);
+ CbmMuchDigi* digi = (CbmMuchDigi*) fDigiManager->Get<CbmMuchDigi>(i);
+ const CbmMatch* match = fDigiManager->GetMatch(ECbmModuleId::kMuch, i);
+ if (!match) {
+ LOG(error) << "digi has no match";
+ return;
+ }
const CbmMuchPad* pad = GetPad(digi->GetAddress());
Double_t 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);
- info->AddCharge(charge);
- info->AddArea(area);
+ const CbmLink& link = match->GetLink(pt);
+ Int_t charge = link.GetWeight();
+ CbmMuchPointInfo& info = getPointInfo(link);
+ info.AddCharge(charge);
+ info.AddArea(area);
}
}
}
@@ -847,16 +879,17 @@ void CbmMuchDigitizerQa::FillDigitizerPerformancePlots() {
return;
}
Bool_t verbose = (fVerbose > 2);
- for (Int_t i = 0; i < fPointInfos->GetEntriesFast(); i++) {
- CbmMuchPointInfo* info = (CbmMuchPointInfo*) fPointInfos->At(i);
+ int iPoint = 0;
+ for (auto it = fMcPointInfoMap.begin(); it != fMcPointInfoMap.end(); it++, ++iPoint) {
+ const CbmMuchPointInfo& info = it->second;
if (verbose) {
- printf("%i", i);
- info->Print();
+ printf("%i", iPoint);
+ info.Print();
}
- Double_t length = info->GetLength();
- Double_t kine = 1000 * (info->GetKine());
- Double_t charge = info->GetCharge();
- Int_t pdg = info->GetPdgCode();
+ Double_t length = info.GetLength();
+ Double_t kine = 1000 * (info.GetKine());
+ Double_t charge = info.GetCharge();
+ Int_t pdg = info.GetPdgCode();
if (charge <= 0) continue;
if (pdg == 22 || // photons
@@ -871,13 +904,13 @@ void CbmMuchDigitizerQa::FillDigitizerPerformancePlots() {
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;
+ 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);
fhNpadsVsS->Fill(area, nPads);
fhTrackCharge->Fill(charge);
- fvTrackCharge[info->GetStationId()]->Fill(charge);
+ fvTrackCharge[info.GetStationId()]->Fill(charge);
fhTrackChargeLog->Fill(log_charge);
fhTrackChargeVsTrackEnergyLog->Fill(log_kine, charge);
fhTrackChargeVsTrackLength->Fill(length, charge);
@@ -1004,21 +1037,26 @@ void CbmMuchDigitizerQa::DrawLengthCanvases() {
// -------------------------------------------------------------------------
void CbmMuchDigitizerQa::PrintFrontLayerPoints() {
- if (!fMCTracks || !fPoints) { 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());
+ if (!fMCTracks || !fTimeSlice) { return; }
+
+ const CbmMatch& match = fTimeSlice->GetMatch();
+ for (int iLink = 0; iLink < match.GetNofLinks(); iLink++) {
+ CbmLink link = match.GetLink(iLink);
+ int nMuchPoints = fPoints->Size(link);
+ for (Int_t ip = 0; ip < nMuchPoints; ip++) {
+ link.SetIndex(ip);
+ CbmMuchPoint* point = (CbmMuchPoint*) fPoints->Get(link);
+ if (!point) {
+ LOG(error) << " Much MC point " << ip << " out of " << nMuchPoints << " doesn't exist";
+ return;
+ }
+ 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", ip, point->GetXIn(), point->GetYIn(),
+ point->GetXOut(), point->GetYOut(), point->GetZIn());
+ }
}
}
diff --git a/sim/detectors/much/qa/CbmMuchDigitizerQa.h b/sim/detectors/much/qa/CbmMuchDigitizerQa.h
index 1ef0a78cd1..75b80d3f86 100644
--- a/sim/detectors/much/qa/CbmMuchDigitizerQa.h
+++ b/sim/detectors/much/qa/CbmMuchDigitizerQa.h
@@ -9,12 +9,17 @@
#ifndef CbmMuchDigitizerQa_H
#define CbmMuchDigitizerQa_H
+#include "CbmLink.h"
+#include "CbmMuchPointInfo.h"
+
#include "FairTask.h"
+
+#include "TFolder.h"
#include "TParameter.h"
-#include <Rtypes.h>
-#include <RtypesCore.h>
-#include <TFolder.h>
+
+#include <map>
#include <vector>
+
class TBuffer;
class TClass;
class TClonesArray;
@@ -28,6 +33,9 @@ class TH1F;
class TH2F;
class TVector2;
class CbmMuchPad;
+class CbmMCDataArray;
+class CbmMCDataManager;
+class CbmTimeSlice;
/// QA for the MUCH detector after a "digitization" step of the simulation.
/// The class reimplements corresponding QA checks from old CbmMuchHitFinderQa class
@@ -85,17 +93,29 @@ private:
TFolder* histFolder; /// folder wich contains histogramms
+ // setup
+ FairRootManager* fManager = nullptr;
+ CbmMCDataManager* fMcManager = nullptr;
+ CbmTimeSlice* fTimeSlice = nullptr;
+
// geometry
CbmMuchGeoScheme* fGeoScheme = nullptr;
Int_t fNstations = 0;
CbmDigiManager* fDigiManager = nullptr;
// containers
- TClonesArray* fPoints = nullptr;
+ CbmMCDataArray* fPoints = nullptr;
+ CbmMCDataArray* fMCTracks = nullptr;
TClonesArray* fDigis = nullptr;
TClonesArray* fDigiMatches = nullptr;
- TClonesArray* fMCTracks = nullptr;
- TClonesArray* fPointInfos = nullptr; /// temporary additional information
+
+ std::map<CbmLink, CbmMuchPointInfo> fMcPointInfoMap = {}; //! map point link -> point info
+
+ CbmMuchPointInfo& getPointInfo(const CbmLink& link)
+ {
+ assert(fMcPointInfoMap.find(link) != fMcPointInfoMap.end());
+ return fMcPointInfoMap[link];
+ }
TFolder fOutFolder; /// output folder with histos and canvases
TParameter<int> fNevents; /// number of processed events
diff --git a/sim/detectors/much/qa/CbmMuchTransportQa.cxx b/sim/detectors/much/qa/CbmMuchTransportQa.cxx
index 032bdd84f8..6a17b4ba83 100644
--- a/sim/detectors/much/qa/CbmMuchTransportQa.cxx
+++ b/sim/detectors/much/qa/CbmMuchTransportQa.cxx
@@ -8,6 +8,9 @@
/// \date 21.10.2020
#include "CbmMuchTransportQa.h"
+
+#include "CbmMCDataArray.h"
+#include "CbmMCDataManager.h"
#include "CbmMCTrack.h"
#include "CbmMuchAddress.h"
#include "CbmMuchGeoScheme.h"
@@ -15,22 +18,26 @@
#include "CbmMuchStation.h"
#include "CbmQaCanvas.h"
#include "CbmQaPie.h"
+#include "CbmTimeSlice.h"
+
+#include <FairRootManager.h>
+#include <FairSink.h>
+#include <FairTask.h>
+#include <Logger.h>
+
#include "TClonesArray.h"
#include "TDatabasePDG.h"
#include "TH1.h"
#include "TH2.h"
#include "TLegend.h"
#include "TStyle.h"
-#include <FairRootManager.h>
-#include <FairSink.h>
-#include <FairTask.h>
-#include <Logger.h>
#include <TAxis.h>
#include <TDirectory.h>
#include <TMath.h>
#include <TParameter.h>
#include <TString.h>
#include <TVector3.h>
+
#include <vector>
#define BINS_STA fNstations, 0, fNstations
@@ -118,16 +125,29 @@ InitStatus CbmMuchTransportQa::Init() {
DeInit();
TDirectory* oldDirectory = gDirectory;
- FairRootManager* manager = FairRootManager::Instance();
- fMcTracks = (TClonesArray*) manager->GetObject("MCTrack");
- fPoints = (TClonesArray*) manager->GetObject("MuchPoint");
- fNstations = CbmMuchGeoScheme::Instance()->GetNStations();
- histFolder = fOutFolder.AddFolder("hist", "Histogramms");
+ fManager = FairRootManager::Instance();
- if (!manager) {
+ if (!fManager) {
LOG(error) << "No FairRootManager found";
return kERROR;
}
+
+ fMcManager = dynamic_cast<CbmMCDataManager*>(fManager->GetObject("MCDataManager"));
+
+ fTimeSlice = (CbmTimeSlice*) fManager->GetObject("TimeSlice.");
+ if (!fTimeSlice) { LOG(error) << GetName() << ": No time slice found"; }
+
+ if (fMcManager) {
+ // Get MCTrack array
+ fMcTracks = fMcManager->InitBranch("MCTrack");
+
+ // Get StsPoint array
+ fPoints = fMcManager->InitBranch("MuchPoint");
+ }
+
+ fNstations = CbmMuchGeoScheme::Instance()->GetNStations();
+ histFolder = fOutFolder.AddFolder("hist", "Histogramms");
+
if (!fMcTracks) {
LOG(error) << "No MC tracks found";
return kERROR;
@@ -357,45 +377,53 @@ void CbmMuchTransportQa::Exec(Option_t*) {
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);
-
- for (Int_t i = 0; i < fPoints->GetEntriesFast(); i++) {
-
- CbmMuchPoint* point = (CbmMuchPoint*) fPoints->At(i);
- Int_t stId = CbmMuchAddress::GetStationIndex(point->GetDetectorID());
- UInt_t stMask = (1 << stId);
- Int_t trackId = point->GetTrackID();
- if (!point) {
- LOG(fatal) << "Much point " << i << " doesn't exist";
- break;
- } // Check if the point corresponds to a certain MC Track
- if (trackId < 0 || trackId >= fMcTracks->GetEntriesFast()) {
- LOG(fatal) << "Much point " << i << ": trackId " << trackId
- << " doesn't belong to [0," << fMcTracks->GetEntriesFast() - 1
- << "]";
- break;
- }
-
- CbmMCTrack* mcTrack = (CbmMCTrack*) fMcTracks->At(trackId);
- if (!mcTrack) {
- LOG(fatal) << "MC track " << trackId << " doesn't exist";
- break;
- }
-
- Int_t motherId = mcTrack->GetMotherId();
- Int_t pdgCode = mcTrack->GetPdgCode();
- if (pdgCode == 22 || // photons
- pdgCode == 2112) // neutrons
- {
- continue;
- }
- if (!(trackStaCross[trackId] & stMask)) {
- FillCountingHistos(stId, motherId, pdgCode);
+ const CbmMatch& sliceMatch = fTimeSlice->GetMatch();
+
+ for (int iLink = 0; iLink < sliceMatch.GetNofLinks(); iLink++) {
+ const CbmLink& eventLink = sliceMatch.GetLink(iLink);
+ int nMuchPoints = fPoints->Size(eventLink);
+ int nMcTracks = fMcTracks->Size(eventLink);
+
+ // bitmask tells which stations were crossed by mc track
+ std::vector<UInt_t> trackStaCross(nMcTracks, 0);
+
+ for (Int_t ip = 0; ip < nMuchPoints; ip++) {
+ CbmLink pointLink = eventLink;
+ pointLink.SetIndex(ip);
+ CbmMuchPoint* point = (CbmMuchPoint*) fPoints->Get(pointLink);
+ Int_t stId = CbmMuchAddress::GetStationIndex(point->GetDetectorID());
+ UInt_t stMask = (1 << stId);
+ Int_t trackId = point->GetTrackID();
+ if (!point) {
+ LOG(fatal) << "Much point " << ip << " doesn't exist";
+ break;
+ } // Check if the point corresponds to a certain MC Track
+ if (trackId < 0 || trackId >= nMcTracks) {
+ LOG(fatal) << "Much point " << ip << ": trackId " << trackId << " doesn't belong to [0," << nMcTracks - 1
+ << "]";
+ break;
+ }
+ CbmLink trackLink = eventLink;
+ trackLink.SetIndex(trackId);
+ CbmMCTrack* mcTrack = (CbmMCTrack*) fMcTracks->Get(trackLink);
+ if (!mcTrack) {
+ LOG(fatal) << "MC track " << trackId << " doesn't exist";
+ break;
+ }
+
+ Int_t motherId = mcTrack->GetMotherId();
+ Int_t pdgCode = mcTrack->GetPdgCode();
+ if (pdgCode == 22 || // photons
+ pdgCode == 2112) // neutrons
+ {
+ continue;
+ }
+
+ if (!(trackStaCross[trackId] & stMask)) { FillCountingHistos(stId, motherId, pdgCode); }
+ trackStaCross[trackId] |= stMask;
+ Fill2dSpatialDistributionHistos(point, stId);
}
- trackStaCross[trackId] |= stMask;
- Fill2dSpatialDistributionHistos(point, stId);
}
}
diff --git a/sim/detectors/much/qa/CbmMuchTransportQa.h b/sim/detectors/much/qa/CbmMuchTransportQa.h
index 73cc74e190..6dcd28c9a1 100644
--- a/sim/detectors/much/qa/CbmMuchTransportQa.h
+++ b/sim/detectors/much/qa/CbmMuchTransportQa.h
@@ -16,6 +16,10 @@
#include <RtypesCore.h>
#include <TFolder.h>
#include <vector>
+
+class CbmMCDataArray;
+class CbmMCDataManager;
+class CbmTimeSlice;
class CbmMuchPoint;
class CbmQaCanvas;
class TBuffer;
@@ -73,9 +77,15 @@ private:
/// geometry
Int_t fNstations = 0;
+ // setup
+ FairRootManager* fManager = nullptr;
+ CbmMCDataManager* fMcManager = nullptr;
+ CbmTimeSlice* fTimeSlice = nullptr;
+
/// containers
- TClonesArray* fPoints = nullptr;
- TClonesArray* fMcTracks = nullptr;
+ CbmMCDataArray* fPoints = nullptr;
+ CbmMCDataArray* fMcTracks = nullptr;
+
///
TFolder* histFolder; /// subfolder for histograms
--
GitLab