diff --git a/core/detectors/tof/CbmTofTrackingInterface.cxx b/core/detectors/tof/CbmTofTrackingInterface.cxx
index 1f11a22e7dfd5d454a7136a8cc5d19eee443868c..a5d4ce763f12683687627d11b6ef5e19da867b61 100644
--- a/core/detectors/tof/CbmTofTrackingInterface.cxx
+++ b/core/detectors/tof/CbmTofTrackingInterface.cxx
@@ -10,7 +10,6 @@
***************************************************************************************************/
#include "CbmTofTrackingInterface.h"
-
#include "CbmTofCreateDigiPar.h"
#include "FairDetector.h"
@@ -98,6 +97,9 @@ InitStatus CbmTofTrackingInterface::Init()
auto chLoBoarderY = chPosY - chSizeY;
auto chUpBoarderY = chPosY + chSizeY;
+ LOG(info) << "iSmType=" << iSmType << ", iSm=" << iSm << ", iRpc=" << iRpc << ": trSta=" << iStation
+ << ", z=" << chPosZ;
+
// Cuts on Bmon and undefined station ID
if (5 == iSmType) { continue; } // Skip Bmon
if (iStation < 0) { continue; }
diff --git a/core/detectors/tof/CbmTofTrackingInterface.h b/core/detectors/tof/CbmTofTrackingInterface.h
index 68a8f8a5eb65ea68b6b0914fc7160d0bcd1f7769..c17316d0cd797a31638d2e7d08d7a66b0e1ee052 100644
--- a/core/detectors/tof/CbmTofTrackingInterface.h
+++ b/core/detectors/tof/CbmTofTrackingInterface.h
@@ -130,13 +130,15 @@ public:
/// FIXME: Replace GetZref(stationId) + 5 with fTofStationZMax[stationId] + 0.5
/// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
/// \return max Z of the station [cm]
- double GetZmax(int stationId) const { return GetZref(stationId) + 5.; }
+ //double GetZmax(int stationId) const { return GetZref(stationId) + 5.; }
+ double GetZmax(int stationId) const { return fTofStationZMax[stationId] + .5; }
/// FIXME: Replace GetZref(stationId) - 5 with fTofStationZMin[stationId] - 0.5
/// @brief Gets min z of the station
/// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
/// @return min Z of the station [cm]
- double GetZmin(int stationId) const { return GetZref(stationId) - 5.; }
+ //double GetZmin(int stationId) const { return GetZref(stationId) - 5.; }
+ double GetZmin(int stationId) const { return fTofStationZMin[stationId] - .5; }
/// @brief Check if station provides time measurements
/// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 8ea34bed88bb474712400cae81f7ff2fb0f2b2c3..87188873f21b30b6d14cea4dd87d9df385ecd50c 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -53,8 +53,6 @@ set(SRCS
catools/CaToolsWFExpression.cxx
catools/CaToolsMaterialHelper.cxx
- qa/CbmTrackerInputQaTrd.cxx
- qa/CbmTrackerInputQaTof.cxx
qa/CbmCaInputQaBase.cxx
qa/CbmCaInputQaMvd.cxx
qa/CbmCaInputQaSts.cxx
diff --git a/reco/L1/CbmCaMCModule.h b/reco/L1/CbmCaMCModule.h
index 526eb622ff1a3878c5831fe4a4f6165bbe206041..d2cb8ed62e068f20ed037c55c697ed86cdf4f136 100644
--- a/reco/L1/CbmCaMCModule.h
+++ b/reco/L1/CbmCaMCModule.h
@@ -384,6 +384,10 @@ namespace cbm::ca
// // ----- Get station index
int iStLoc = fvpDetInterface[DetID]->GetTrackingStationIndex(pExtPoint);
+ if (iStLoc < 0) {
+ return std::nullopt;
+ }
+
int stationID = fpParameters->GetStationIndexActive(iStLoc, DetID);
if (stationID == -1) {
return std::nullopt;
diff --git a/reco/L1/CbmCaTimeSliceReader.cxx b/reco/L1/CbmCaTimeSliceReader.cxx
index da259489940e2b6b1065d60c4c1e09bff3068ef8..f396a1f638dcf9d57904616eb79fe1519951db6e 100644
--- a/reco/L1/CbmCaTimeSliceReader.cxx
+++ b/reco/L1/CbmCaTimeSliceReader.cxx
@@ -377,6 +377,9 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector()
auto* pHit = static_cast<Hit_t*>(fvpBrHits[DetID]->At(iHext));
int iStGeom = pDetInterface->GetTrackingStationIndex(pHit);
+ if (iStGeom < 0) {
+ continue; // NOTE: sensors with iStGeom = -1 are ignored in tracking
+ }
// Fill out detector specific data
if constexpr (ca::EDetectorID::kSts == DetID) {
@@ -389,11 +392,6 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector()
if (5 == CbmTofAddress::GetSmType(pHit->GetAddress())) {
continue;
}
-
- // FIXME: Figure it out, if this cut is still needed (introduced a year ago for mCBM)
- if (ECbmCaTrackingMode::kMCBM == fTrackingMode && pHit->GetZ() > 400) {
- continue;
- }
}
int iStActive = fpParameters->GetStationIndexActive(iStGeom, DetID);
diff --git a/reco/L1/qa/CbmCaInputQaBase.cxx b/reco/L1/qa/CbmCaInputQaBase.cxx
index 80a15ad140e379f636dfa4c3b7fcd33af4a133a9..613c1c1d91a6c3e86d394ff460dcc9a6f5107afd 100644
--- a/reco/L1/qa/CbmCaInputQaBase.cxx
+++ b/reco/L1/qa/CbmCaInputQaBase.cxx
@@ -326,13 +326,13 @@ void CbmCaInputQaBase<DetID>::FillHistograms()
fHitQaData.SetHitIndex(iH);
const auto* pHit = dynamic_cast<const CbmPixelHit*>(fpHits->At(iH));
- LOG_IF(fatal, !pHit) << fName << ": hit with iH = " << iH << " is not an CbmStsHit (dynamic cast failed)";
+ if (!pHit) {
+ LOG(error) << fName << ": hit with iH = " << iH << " is not an CbmStsHit (dynamic cast failed)";
+ continue;
+ }
if constexpr (ca::EDetectorID::kTof == DetID) {
auto address = pHit->GetAddress();
- if (0x00202806 == address || 0x00002806 == address) {
- continue;
- } // TEST
if (5 == CbmTofAddress::GetSmType(address)) {
continue;
} // Skip Bmon hits
@@ -346,8 +346,14 @@ void CbmCaInputQaBase<DetID>::FillHistograms()
// Hit station geometry info
int iSt = fpDetInterface->GetTrackingStationIndex(pHit);
- LOG_IF(fatal, iSt < 0 || iSt >= nSt) << fName << ": index of station (" << iSt << ") is out of range "
- << "for hit with id = " << iH;
+ if (iSt < 0) {
+ continue;
+ }
+
+ if (iSt >= nSt) {
+ LOG(error) << fName << ": index of station (" << iSt << ") is out of range for hit with id = " << iH;
+ continue;
+ }
auto [stPhiU, stPhiV] = fpDetInterface->GetStereoAnglesSensor(pHit->GetAddress());
@@ -417,19 +423,26 @@ void CbmCaInputQaBase<DetID>::FillHistograms()
// TODO: debug
//if (iE < 0) continue;
- LOG_IF(fatal, iE < 0 || iE >= nMCevents)
- << fName << ": id of MC event is out of range (hit id = " << iH << ", link id = " << iLink
- << ", event id = " << iE << ", mc point ID = " << iP << ')';
+ if (iE < 0 || iE >= nMCevents) {
+ LOG(error) << fName << ": id of MC event is out of range (hit id = " << iH << ", link id = " << iLink
+ << ", event id = " << iE << ", mc point ID = " << iP << ')';
+ continue;
+ }
// matched point
const auto* pMCPoint = dynamic_cast<const Point_t*>(fpMCPoints->Get(link));
- LOG_IF(fatal, !pMCPoint) << fName << ": MC point object does not exist for hit " << iH;
+ if (!pMCPoint) {
+ LOG(error) << fName << ": MC point object does not exist for hit " << iH;
+ continue;
+ }
int iTr = pMCPoint->GetTrackID();
- LOG_IF(fatal, iTr >= (int) vNofHitsPerMcTrack[iE][iSt].size())
- << fName << ": index of MC track is out of range (hit id = " << iH << ", link id = " << iLink
- << ", event id = " << iE << ", track id = " << iTr << ')';
+ if (iTr >= static_cast<int>(vNofHitsPerMcTrack[iE][iSt].size())) {
+ LOG(error) << fName << ": index of MC track is out of range (hit id = " << iH << ", link id = " << iLink
+ << ", event id = " << iE << ", track id = " << iTr << ')';
+ continue;
+ }
vNofHitsPerMcTrack[iE][iSt][iTr]++;
}
@@ -458,16 +471,25 @@ void CbmCaInputQaBase<DetID>::FillHistograms()
// Point matched by the best link
const auto* pMCPoint = dynamic_cast<const Point_t*>(fpMCPoints->Get(bestPointLink));
fHitQaData.SetPointID(bestPointLink.GetIndex(), bestPointLink.GetEntry(), bestPointLink.GetFile());
- LOG_IF(fatal, !pMCPoint) << fName << ": MC point object does not exist for hit " << iH;
+ if (!pMCPoint) {
+ LOG(error) << fName << ": MC point object does not exist for hit " << iH;
+ continue;
+ }
// MC track for this point
CbmLink bestTrackLink = bestPointLink;
bestTrackLink.SetIndex(pMCPoint->GetTrackID());
const auto* pMCTrack = dynamic_cast<const CbmMCTrack*>(fpMCTracks->Get(bestTrackLink));
- LOG_IF(fatal, !pMCTrack) << fName << ": MC track object does not exist for hit " << iH << " and link: ";
+ if (!pMCTrack) {
+ LOG(error) << fName << ": MC track object does not exist for hit " << iH << " and link: ";
+ continue;
+ }
double t0MC = fpMCEventList->GetEventTime(bestPointLink);
- LOG_IF(fatal, t0MC < 0) << fName << ": MC time zero is lower then 0 ns: " << t0MC;
+ if (t0MC < 0) {
+ LOG(error) << fName << ": MC time zero is lower then 0 ns: " << t0MC;
+ continue;
+ }
// cut on the mc track quality
if (!IsTrackSelected(pMCTrack, pMCPoint)) {
@@ -617,12 +639,19 @@ void CbmCaInputQaBase<DetID>::FillHistograms()
fMonitor.IncrementCounter(EMonitorKey::kMcPoint);
const auto* pMCPoint = dynamic_cast<const Point_t*>(fpMCPoints->Get(iFile, iEvent, iP));
- LOG_IF(fatal, !pMCPoint) << fName << ": MC point does not exist for iFile = " << iFile
- << ", iEvent = " << iEvent << ", iP = " << iP;
+ if (!pMCPoint) {
+ LOG(error) << fName << ": MC point does not exist for iFile = " << iFile << ", iEvent = " << iEvent
+ << ", iP = " << iP;
+ continue;
+ }
int address = pMCPoint->GetDetectorID();
int iSt = fpDetInterface->GetTrackingStationIndex(pMCPoint);
- if (iSt < 0 || iSt >= nSt) {
+ if (iSt < 0) {
+ continue; // NOTE: legit, such sensors must be ignored in tracking
+ }
+
+ if (iSt >= nSt) {
LOG(error) << fName << ": MC point for FEI = " << iFile << ", " << iEvent << ", " << iP << " and address "
<< address << " has wrong station index: iSt = " << iSt;
fMonitor.IncrementCounter(EMonitorKey::kMcPointWrongStation);
@@ -645,16 +674,21 @@ void CbmCaInputQaBase<DetID>::FillHistograms()
int iTr = pMCPoint->GetTrackID();
- LOG_IF(fatal, iTr >= nTracks) << fName << ": index of MC track is out of range (point id = " << iP
- << ", event id = " << iE << ", track id = " << iTr << ')';
-
+ if (iTr >= nTracks) {
+ LOG(error) << fName << ": index of MC track is out of range (point id = " << iP << ", event id = " << iE
+ << ", track id = " << iTr << ')';
+ continue;
+ }
const auto* pMCTrack = dynamic_cast<const CbmMCTrack*>(fpMCTracks->Get(iFile, iEvent, iTr));
fHitQaData.SetIfTrackHasHits(vNofHitsPerMcTrack[iE][iSt][iTr] > 0);
fHitQaData.SetIfTrackSelected(IsTrackSelected(pMCTrack, pMCPoint));
- LOG_IF(fatal, pMCTrack == nullptr) << fName << ": null MC track pointer for file id = " << iFile
- << ", event id = " << iEvent << ", track id = " << iTr;
+ if (!pMCTrack) {
+ LOG(error) << fName << ": null MC track pointer for file id = " << iFile << ", event id = " << iEvent
+ << ", track id = " << iTr;
+ continue;
+ }
// check efficiency only once per mc track per station
if (vIsTrackProcessed[iSt][iTr]) {
diff --git a/reco/L1/qa/CbmCaInputQaSetup.cxx b/reco/L1/qa/CbmCaInputQaSetup.cxx
index 4078b7df27758bf9e3396ca288c98ebabcc9c8cf..4cec7abf9afafe4a3cc1c2cb66cb3a8f4d32cea2 100644
--- a/reco/L1/qa/CbmCaInputQaSetup.cxx
+++ b/reco/L1/qa/CbmCaInputQaSetup.cxx
@@ -70,6 +70,10 @@ void InputQaSetup::FillHistogramsDet()
}
int iStLoc = fvpDetInterface[DetID]->GetTrackingStationIndex(address);
+ if (iStLoc < 0) {
+ continue;
+ }
+
int iStGeo = fpParameters->GetStationIndexGeometry(iStLoc, DetID);
auto xHit = pHit->GetX();
auto yHit = pHit->GetY();
diff --git a/reco/L1/qa/CbmCaInputQaTof.cxx b/reco/L1/qa/CbmCaInputQaTof.cxx
index 226374c5c4a6ca68957431fa258350610a9dbb43..497e0d6f254c44f70e8c3af82150df549fd7d17f 100644
--- a/reco/L1/qa/CbmCaInputQaTof.cxx
+++ b/reco/L1/qa/CbmCaInputQaTof.cxx
@@ -112,7 +112,10 @@ void CbmCaInputQaTof::FillHistograms()
LOG(info) << fName << ": ===== Hit Sample";
for (int iHit = 0; iHit < fpHits->GetEntriesFast(); ++iHit) {
const auto* pHit = dynamic_cast<const Hit_t*>(fpHits->At(iHit));
- LOG_IF(fatal, !pHit) << fName << "::FillHistogramsPerHit: hit with index " << iHit << " not found";
+ if (!pHit) {
+ LOG(error) << fName << "::FillHistogramsPerHit: hit with index " << iHit << " not found";
+ continue;
+ }
int address = pHit->GetAddress();
int iHitSmType = CbmTofAddress::GetSmType(address);
@@ -170,7 +173,10 @@ void CbmCaInputQaTof::FillHistogramsPerHit()
{
auto iHit = GetHitQaData().GetHitIndex();
const auto* pHit = dynamic_cast<const Hit_t*>(fpHits->At(iHit));
- LOG_IF(fatal, !pHit) << fName << "::FillHistogramsPerHit: hit with index " << iHit << " not found";
+ if (!pHit) {
+ LOG(error) << fName << "::FillHistogramsPerHit: hit with index " << iHit << " not found";
+ return;
+ }
int address = pHit->GetAddress();
int iHitSmType = CbmTofAddress::GetSmType(address);
diff --git a/reco/L1/qa/CbmTrackerInputQaTof.cxx b/reco/L1/qa/CbmTrackerInputQaTof.cxx
deleted file mode 100644
index dd8a63c122a4d95ac2dd140faa4ac9ced28b7f93..0000000000000000000000000000000000000000
--- a/reco/L1/qa/CbmTrackerInputQaTof.cxx
+++ /dev/null
@@ -1,739 +0,0 @@
-/* Copyright (C) 2021 Facility for Antiproton and Ion Research in Europe, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov[committer]*/
-
-/// @file CbmTrackerInputQaTof.cxx
-/// @author Sergey Gorbunov
-/// @date 02.11.2021
-
-#include "CbmTrackerInputQaTof.h"
-
-#include "CbmDefs.h"
-#include "CbmLink.h"
-#include "CbmMCDataArray.h"
-#include "CbmMCDataManager.h"
-#include "CbmMCEventList.h"
-#include "CbmMCTrack.h"
-#include "CbmMatch.h"
-#include "CbmQaCanvas.h"
-#include "CbmTofCell.h"
-#include "CbmTofDigiPar.h" // in tof/TofParam
-#include "CbmTofTrackingInterface.h"
-//#include "CbmSetup.h"
-#include "CbmTimeSlice.h"
-//#include "CbmTofCluster.h"
-#include "CbmTofDigi.h"
-#include "CbmTofHit.h"
-//#include "CbmTofParModDigi.h" // for CbmTofModule
-//#include "CbmTofParModGeo.h"
-//#include "CbmTofParSetDigi.h" // for CbmTofParSetDigi
-//#include "CbmTofParSetGeo.h"
-#include "CbmQaUtil.h"
-#include "CbmTofPoint.h"
-
-#include <FairRootManager.h>
-#include <FairRunAna.h>
-#include <FairRuntimeDb.h>
-#include <FairSink.h>
-#include <FairTask.h>
-#include <Logger.h>
-
-#include <TClonesArray.h>
-#include <TDatabasePDG.h>
-#include <TGenericClassInfo.h>
-#include <TGeoManager.h>
-#include <TGeoNode.h>
-#include <TMathBase.h>
-#include <TObjArray.h>
-#include <TParameter.h>
-#include <TParticlePDG.h>
-#include <TString.h>
-#include <TStyle.h>
-
-#include <algorithm>
-#include <cstdio>
-#include <iostream>
-#include <map>
-#include <vector>
-
-ClassImp(CbmTrackerInputQaTof);
-
-// -------------------------------------------------------------------------
-
-CbmTrackerInputQaTof::CbmTrackerInputQaTof(const char* name, Int_t verbose) : FairTask(name, verbose)
-{
- // Constructor
-
- // Create a list of histogramms
-
- fHistList.clear();
- fHistList.reserve(20);
- fHistList.push_back(&fhResidualX);
- fHistList.push_back(&fhResidualY);
- fHistList.push_back(&fhResidualT);
- //fHistList.push_back(&fhResidualZ);
- fHistList.push_back(&fhPullX);
- fHistList.push_back(&fhPullY);
- fHistList.push_back(&fhPullT);
-
-
- // Keep the ownership on the histograms to avoid double deletion
- for (unsigned int i = 0; i < fHistList.size(); i++) {
- fHistList[i]->SetDirectory(nullptr);
- }
-}
-
-// -------------------------------------------------------------------------
-CbmTrackerInputQaTof::~CbmTrackerInputQaTof()
-{ /// Destructor
- DeInit();
-}
-
-
-// -------------------------------------------------------------------------
-void CbmTrackerInputQaTof::DeInit()
-{
- fIsTofInSetup = 0;
- fIsMcPresent = false;
- fNtrackingStations = 0;
-
- fTimeSlice = nullptr;
-
- fMcManager = nullptr;
- fMcTracks = nullptr;
- fMcPoints = nullptr;
-
- fClusters = nullptr;
- fHits = nullptr;
- fHitMatches = nullptr;
-
- fOutFolder.Clear();
-
- fHistFolder = nullptr;
-
- fNevents.SetVal(0);
-
- fhPointsPerHit.clear();
- fhHitsPerPoint.clear();
-}
-// -------------------------------------------------------------------------
-
-void CbmTrackerInputQaTof::SetParContainers()
-{
- fTofDigiPar = nullptr;
- fTofGeoPar = nullptr;
-
- // Get run and runtime database
- FairRunAna* ana = FairRunAna::Instance();
- if (!ana) {
- LOG(fatal) << "No FairRunAna instance";
- }
-
- FairRuntimeDb* rtdb = ana->GetRuntimeDb();
- if (!rtdb) {
- LOG(fatal) << "No FairRuntimeDb instance";
- }
-
- // fTofDigiPar = dynamic_cast<CbmTofDigiPar*>(rtdb->getContainer("CbmTofDigiPar"));
- // fTofGeoPar = dynamic_cast<CbmTofParSetGeo*>(rtdb->getContainer("CbmTofParSetGeo"));
-
- fTofDigiPar = (CbmTofDigiPar*) rtdb->getContainer("CbmTofDigiPar");
- // fTofGeoPar = rtdb->getContainer("CbmGeoTofPar");
-}
-
-// -------------------------------------------------------------------------
-InitStatus CbmTrackerInputQaTof::ReInit()
-{
- // Initialize and check the setup
-
- DeInit();
-
- // Check if Tof is present in the CBM setup
- // A straightforward way to do so is CbmSetup::Instance()->IsActive(ECbmModuleId::kTof),
- // but unfortunatelly CbmSetup class is unaccesible.
- // ( CbmSimSteer library requires libfreetype that has linking problems on MacOsX )
- // Therefore let's simply check if any Tof material is present in the geometry.
- //
- fIsTofInSetup = 0;
- {
- TObjArray* topNodes = gGeoManager->GetTopNode()->GetNodes();
- for (Int_t iTopNode = 0; iTopNode < topNodes->GetEntriesFast(); iTopNode++) {
- TGeoNode* topNode = static_cast<TGeoNode*>(topNodes->At(iTopNode));
- if (TString(topNode->GetName()).Contains("tof")) {
- fIsTofInSetup = 1;
- break;
- }
- }
- }
-
- if (!fIsTofInSetup) {
- LOG(info) << "Tof is not in the setup, do nothing";
- return kSUCCESS;
- }
-
- // check parameter containers
-
- if (!fTofDigiPar) {
- LOG(error) << "No CbmTofParSetDigi container in FairRuntimeDb";
- return kERROR;
- }
-
- // if (!fTofGeoPar) {
- // LOG(error) << "No CbmTofParSetGeo container in FairRuntimeDb";
- // return kERROR;
- // }
-
- FairRootManager* manager = FairRootManager::Instance();
-
- fTimeSlice = dynamic_cast<CbmTimeSlice*>(manager->GetObject("TimeSlice."));
- if (!fTimeSlice) {
- LOG(error) << "No time slice found";
- return kERROR;
- }
-
- fHits = dynamic_cast<TClonesArray*>(manager->GetObject("TofHit"));
-
- if (!fHits) {
- LOG(error) << "No Tof hit array found";
- return kERROR;
- }
-
- // fClusters = dynamic_cast<TClonesArray*>(manager->GetObject("TofCluster"));
- //
- // if (!fClusters) {
- // LOG(error) << "No Tof cluster array found";
- // return kERROR;
- // }
-
- fMcManager = dynamic_cast<CbmMCDataManager*>(manager->GetObject("MCDataManager"));
-
- fIsMcPresent = (fMcManager != nullptr);
-
- if (fIsMcPresent) {
-
- fMcEventList = dynamic_cast<CbmMCEventList*>(manager->GetObject("MCEventList."));
- fMcTracks = fMcManager->InitBranch("MCTrack");
- fMcPoints = fMcManager->InitBranch("TofPoint");
- fHitMatches = dynamic_cast<TClonesArray*>(manager->GetObject("TofHitMatch"));
-
- fHitDigiMatches = dynamic_cast<TClonesArray*>(manager->GetObject("TofDigiMatch"));
- fDigis = dynamic_cast<TClonesArray*>(manager->GetObject("TofCalDigi"));
- if (nullptr == fDigis) {
- LOG(info) << "No calibrated tof digi vector in the input files => trying original vector";
- fDigis = dynamic_cast<TClonesArray*>(manager->GetObject("TofDigi"));
- if (nullptr == fDigis) {
- LOG(info) << "No original tof digi vector in the input files! Ignore TOF!";
- }
- }
- fTofCalDigiMatch = dynamic_cast<TClonesArray*>(manager->GetObject("TofCalDigiMatch"));
-
-
- // we assume that when Tof is in the setup and the MC manager is present,
- // the Tof MC data must be present too
-
- if (!fMcEventList) {
- LOG(error) << ": No MCEventList data!";
- return kERROR;
- }
-
- if (!fMcTracks) {
- LOG(error) << "No MC tracks found";
- return kERROR;
- }
-
- if (!fMcPoints) {
- LOG(error) << "No MC points found";
- return kERROR;
- }
-
- if (!fHitMatches) {
- LOG(error) << "No Tof hit matches found";
- return kERROR;
- }
-
- // if (!fDigis) {
- // LOG(error) << "No Tof digis found";
- // return kERROR;
- // }
- }
-
- // init geometry setup and analyse it
-
- auto retVal = GeometryQa();
-
- if (fNtrackingStations <= 0) {
- LOG(error) << "can not count Tof tracking stations";
- return kERROR;
- }
-
- // initialise histograms
- fOutFolder.SetOwner(false);
- fHistFolder = fOutFolder.AddFolder("rawHist", "Raw Histograms");
- gStyle->SetOptStat(0);
-
- fNevents.SetVal(0);
- fHistFolder->Add(&fNevents);
-
- for (unsigned int i = 0; i < fHistList.size(); i++) {
- fHistList[i]->Reset();
- fHistFolder->Add(fHistList[i]);
- }
-
- fhPointsPerHit.clear();
- fhHitsPerPoint.clear();
- fhEfficiencyR.clear();
- fhEfficiencyXY.clear();
-
- fhPointsPerHit.reserve(fNtrackingStations);
- fhHitsPerPoint.reserve(fNtrackingStations);
- fhEfficiencyR.reserve(fNtrackingStations);
- fhEfficiencyXY.reserve(fNtrackingStations);
-
- for (Int_t i = 0; i < fNtrackingStations; i++) {
- fhPointsPerHit.emplace_back(Form("hMcPointsPerHit%i", i),
- Form("MC Points per Hit: Station %i;N mc points / hit", i), 10, -0.5, 9.5);
- fhPointsPerHit[i].SetDirectory(nullptr);
- fhPointsPerHit[i].SetLineWidth(2);
- fhPointsPerHit[i].SetOptStat(110);
- fHistFolder->Add(&fhPointsPerHit[i]);
- }
-
- for (Int_t i = 0; i < fNtrackingStations; i++) {
- fhHitsPerPoint.emplace_back(Form("hHitsPerMcPoint%i", i),
- Form("Hits per MC Point: Station %i; N hits / mc point", i), 10, -0.5, 9.5);
- fhHitsPerPoint[i].SetDirectory(nullptr);
- fhHitsPerPoint[i].SetLineWidth(2);
- fhHitsPerPoint[i].SetOptStat(110);
- fHistFolder->Add(&fhHitsPerPoint[i]);
- }
-
- for (Int_t i = 0; i < fNtrackingStations; i++) {
-
- double dx = 350;
- double dy = 350;
- double dr = sqrt(dx * dx + dy * dy);
-
- fhEfficiencyR.emplace_back(Form("hEfficiencyR%i", i), Form("Efficiency R: Station %i;R [cm]", i), 100, 0, dr);
- fhEfficiencyR[i].SetDirectory(nullptr);
- fhEfficiencyR[i].SetOptStat(1110);
- fHistFolder->Add(&fhEfficiencyR[i]);
-
- fhEfficiencyXY.emplace_back(Form("hEfficiencyXY%i", i), Form("Efficiency XY: Station %i;X [cm];Y [cm]", i), 50, -dx,
- dx, 50, -dy, dy);
- fhEfficiencyXY[i].SetDirectory(nullptr);
- fhEfficiencyXY[i].SetOptStat(10);
- fhEfficiencyXY[i].GetYaxis()->SetTitleOffset(1.4);
- fHistFolder->Add(&fhEfficiencyXY[i]);
- }
-
- fCanvResidual.Clear();
- fCanvResidual.Divide2D(6);
-
- fCanvPull.Clear();
- fCanvPull.Divide2D(6);
-
- fCanvEfficiencyR.Clear();
- fCanvEfficiencyR.Divide2D(fNtrackingStations);
-
- fCanvEfficiencyXY.Clear();
- fCanvEfficiencyXY.Divide2D(fNtrackingStations);
-
- fCanvPointsPerHit.Clear();
- fCanvPointsPerHit.Divide2D(fNtrackingStations);
-
- fCanvHitsPerPoint.Clear();
- fCanvHitsPerPoint.Divide2D(fNtrackingStations);
-
- fOutFolder.Add(&fCanvResidual);
- fOutFolder.Add(&fCanvPull);
- fOutFolder.Add(&fCanvEfficiencyR);
- fOutFolder.Add(&fCanvEfficiencyXY);
- fOutFolder.Add(&fCanvPointsPerHit);
- fOutFolder.Add(&fCanvHitsPerPoint);
-
- return retVal;
-}
-
-
-// -------------------------------------------------------------------------
-void CbmTrackerInputQaTof::Exec(Option_t*)
-{
- if (!fIsTofInSetup) {
- return;
- }
-
- // update number of processed events
- fNevents.SetVal(fNevents.GetVal() + 1);
- LOG(debug) << "Event: " << fNevents.GetVal();
- ResolutionQa();
-}
-
-
-// -------------------------------------------------------------------------
-void CbmTrackerInputQaTof::Finish()
-{
- FairSink* sink = FairRootManager::Instance()->GetSink();
- if (sink) {
- sink->WriteObject(&GetQa(), nullptr);
- }
-}
-
-int CbmTrackerInputQaTof::GetStationIndex(CbmTofPoint* p)
-{
- // if (p->GetZ() > 800) return -1;
- float dist = 1000;
- int iSta = -1;
- auto tofInterface = CbmTofTrackingInterface::Instance();
- for (int iSt = 0; iSt < fNtrackingStations; iSt++) {
- if (fabs(p->GetZ() - tofInterface->GetZref(iSt)) < dist) {
- dist = fabs(p->GetZ() - tofInterface->GetZref(iSt));
- iSta = iSt;
- }
- }
- return iSta;
-}
-
-
-// -------------------------------------------------------------------------
-InitStatus CbmTrackerInputQaTof::GeometryQa()
-{
- // check geometry of the Tof modules
-
- auto tofInterface = CbmTofTrackingInterface::Instance();
-
- fNtrackingStations = tofInterface->GetNtrackingStations();
-
- double lastZ = -1;
- for (int iStation = 0; iStation < fNtrackingStations; iStation++) {
-
- // tofInterface->GetZref(iStation);
- // tofInterface->GetTimeResolution(iStation);
- // tofInterface->IsTimeInfoProvided(iStation);
-
- double staZ = tofInterface->GetZref(iStation);
-
- // check that the stations are properly ordered in Z
- if (((iStation > 0) && (staZ <= lastZ)) || ((staZ != staZ))) {
- LOG(error) << "Tof trackig stations are not properly ordered in Z";
- return kERROR;
- }
- lastZ = staZ;
- }
-
- return kSUCCESS;
-}
-
-// -------------------------------------------------------------------------
-void CbmTrackerInputQaTof::ResolutionQa()
-{
-
- if (!fHitMatches || !fHits || !fMcEventList) return;
-
- Int_t nHits = fHits->GetEntriesFast();
- Int_t nHitMatches = fHitMatches->GetEntriesFast();
-
- if (nHits != nHitMatches) {
- LOG(fatal) << "CbmTrackerInputQaTof::ResolutionQa => N hits " << nHits << " not matching N matches " << nHitMatches;
- }
-
- auto tofInterface = CbmTofTrackingInterface::Instance();
-
- fNtrackingStations = tofInterface->GetNtrackingStations();
-
- int nMcEvents = fMcEventList->GetNofEvents();
-
- // TOF hit is usually created out of several MC points,
- // but is linked to only one of those points.
- // It should be taken into account.
-
- // N TOF hits per MC point
- std::vector<int> pointNhits[nMcEvents];
-
- // choose only one MC point per plane per MC track
- std::vector<int> selectedPoints[nMcEvents];
-
- //LOG(info) << "CbmTrackerInputQaTof:: timeslice nHits " << nHits << " nHitMatches " << nHitMatches;
- //LOG(info) << "CbmTrackerInputQaTof:: start processing " << nMcEvents << " events.. ";
-
- for (int iE = 0; iE < nMcEvents; iE++) {
- int fileId = fMcEventList->GetFileIdByIndex(iE);
- int eventId = fMcEventList->GetEventIdByIndex(iE);
-
- const int nMcTracks = fMcTracks->Size(fileId, eventId);
- const int nMcPoints = fMcPoints->Size(fileId, eventId);
-
- //LOG(info) << "CbmTrackerInputQaTof:: process Mc event " << iE << " fileId " << fileId << " eventId " << eventId << " nMcTracks "
- // << nMcTracks << " nMcPoints " << nMcPoints;
-
- pointNhits[iE].resize(nMcPoints, 0);
- selectedPoints[iE].clear();
- selectedPoints[iE].reserve(nMcPoints);
-
- std::vector<char> selectedMcPointPerStation[fNtrackingStations];
-
- for (Int_t iS = 0; iS < fNtrackingStations; iS++) {
- selectedMcPointPerStation[iS].resize(nMcTracks, 0);
- }
-
- // check which points of the current MC event are matched to hits,
- // fill N hits per mc point
- std::vector<char> isPointMatched(nMcPoints, 0);
- for (Int_t iHit = 0; iHit < nHits; iHit++) { // loop over hits in the time slice
-
- CbmMatch* pHitMatch = static_cast<CbmMatch*>(fHitMatches->At(iHit));
-
- for (int iLink = 0; iLink < pHitMatch->GetNofLinks(); iLink++) {
- int iFile = pHitMatch->GetLink(iLink).GetFile();
- int iEvent = pHitMatch->GetLink(iLink).GetEntry();
- int iMC = pHitMatch->GetLink(iLink).GetIndex();
- // check if the link belongs to the current Mc event iE
- if ((fileId != iFile) || (eventId != iEvent)) continue;
- assert(iMC >= 0 && iMC < nMcPoints);
- isPointMatched[iMC] = 1;
- pointNhits[iE][iMC]++;
- }
- }
-
- // select the matched points first
- for (Int_t iMC = 0; iMC < nMcPoints; iMC++) {
- if (isPointMatched[iMC] == 0) continue;
- CbmTofPoint* p = dynamic_cast<CbmTofPoint*>(fMcPoints->Get(fileId, eventId, iMC));
- int PointStation = GetStationIndex(p);
- if (PointStation < 0) continue;
- if (!selectedMcPointPerStation[PointStation][p->GetTrackID()]) {
- selectedPoints[iE].push_back(iMC);
- }
- selectedMcPointPerStation[PointStation][p->GetTrackID()] = 1;
- }
-
- // select the unmatched points if there is no corresponding matched point
- for (Int_t iMC = 0; iMC < nMcPoints; iMC++) {
- CbmTofPoint* p = dynamic_cast<CbmTofPoint*>(fMcPoints->Get(fileId, eventId, iMC));
- int PointStation = GetStationIndex(p);
- if (PointStation < 0) continue;
- if (!selectedMcPointPerStation[PointStation][p->GetTrackID()]) selectedPoints[iE].push_back(iMC);
- selectedMcPointPerStation[PointStation][p->GetTrackID()] = 1;
- }
- }
-
-
- // loop over hit in the time slice
- for (Int_t iHit = 0; iHit < nHits; iHit++) {
-
- CbmTofHit* hit = dynamic_cast<CbmTofHit*>(fHits->At(iHit));
- if (!hit) {
- LOG(fatal) << "Tof hit N " << iHit << " doesn't exist";
- }
- const int address = hit->GetAddress();
-
- //should ignore these hits for some reason
- if (0x00202806 == address || 0x00002806 == address) continue;
-
- int StationIndex = tofInterface->GetTrackingStationIndex(hit);
-
- if (StationIndex < 0 || StationIndex >= fNtrackingStations) {
- LOG(fatal) << "Tof hit layer id " << StationIndex << " is out of range";
- }
-
- CbmMatch* myHitMatch = dynamic_cast<CbmMatch*>(fHitMatches->At(iHit));
- if (!myHitMatch) {
- LOG(fatal) << "Tof hit match N " << iHit << " doesn't exist";
- }
-
- // take corresponding MC point
-
- const CbmLink& bestLink = myHitMatch->GetMatchedLink();
-
- // skip hits from the noise digis
- if (bestLink.GetIndex() < 0) {
- continue;
- }
-
- CbmTofPoint* p = dynamic_cast<CbmTofPoint*>(fMcPoints->Get(bestLink));
- if (!p) {
- LOG(error) << "CbmTrackerInputQaTof:: link points to a non-existing MC point";
- }
-
- if (p->GetDetectorID() != (int) hit->GetAddress()) {
- LOG(error) << "CbmTrackerInputQaTof:: mc point module address differs from the hit module address";
- }
-
- // how many MC points?
- int nHitPoints = 0;
- for (int i = 0; i < myHitMatch->GetNofLinks(); i++) {
- const CbmLink& link = myHitMatch->GetLink(i);
- if (link.GetIndex() >= 0) {
- int iE = fMcEventList->GetEventIndex(link);
- int indexPoint = link.GetIndex();
- bool mcPointValid = 0;
- int nPointsEv = selectedPoints[iE].size();
- for (int ip = 0; ip < nPointsEv; ip++) {
- if (selectedPoints[iE][ip] == indexPoint) mcPointValid = 1;
- }
- if (!mcPointValid) continue;
- nHitPoints++;
- }
- }
- fhPointsPerHit[StationIndex].Fill(nHitPoints);
-
- // check that the nominal station Z is not far from the active material
-
- // the cut of 1 cm is choosen arbitrary and can be changed
- //
- // const CbmTofParModGeo* pGeo = dynamic_cast<const CbmTofParModGeo*>(fTofGeoPar->GetModulePar(ModuleId));
- // if (!pGeo) {
- // LOG(fatal) << " No Geo params for station " << StationIndex << " module " << ModuleId << " found ";
- // return;
- // }
-
- // double staZ = tofInterface->GetZref(StationIndex); // module->GetZ(); //+ 410;
- // if ((staZ < p->GetZ() - 1.) || (staZ > p->GetZ() + 1.)) {
- // LOG(error) << "Tof station " << StationIndex << ": active material Z[" << p->GetZ() << " cm," << p->GetZ()
- // << " cm] is too far from the nominal station Z " << staZ << " cm";
- // return;
- // }
- // // the cut of 1 cm is choosen arbitrary and can be changed
- // if (fabs(hit->GetZ() - staZ) > 1.) {
- // LOG(error) << "Tof station " << StationIndex << ": hit Z " << hit->GetZ()
- // << " cm, is too far from the nominal station Z " << staZ << " cm";
- // return;
- // }
- // }
-
- // residual and pull
-
- if (nHitPoints != 1) continue; // only check residual for non-mixed hits
-
- double t0 = fMcEventList->GetEventTime(bestLink);
- if (t0 < 0) {
- LOG(error) << "MC event time doesn't exist for a Tof link";
- return;
- }
-
- double x = p->GetX(); // take the "In" point since the time is stored only for this point
- double y = p->GetY();
- double z = p->GetZ();
- double t = t0 + p->GetTime();
- double px = p->GetPx();
- double py = p->GetPy();
- double pz = p->GetPz();
-
- // skip very slow particles
- if (fabs(p->GetPz()) < 0.01) continue;
-
- // transport the particle from the MC point to the hit Z
- double dz = hit->GetZ() - z;
- x += dz * px / pz;
- y += dz * py / pz;
-
- fhResidualZ.Fill(dz);
-
- // get particle mass
- double mass = 0;
- {
- CbmLink mcTrackLink = bestLink;
- mcTrackLink.SetIndex(p->GetTrackID());
- CbmMCTrack* mcTrack = dynamic_cast<CbmMCTrack*>(fMcTracks->Get(mcTrackLink));
- if (!mcTrack) {
- LOG(error) << "MC track " << p->GetTrackID() << " doesn't exist";
- return;
- }
- Int_t pdg = mcTrack->GetPdgCode();
- if (pdg < 9999999 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) {
- mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
- }
- }
-
- constexpr double speedOfLight = 29.979246; // cm/ns
- TVector3 mom3;
- p->Momentum(mom3);
- t += dz / (pz * speedOfLight) * sqrt(mass * mass + mom3.Mag2());
-
- double dx = hit->GetX() - x;
- double dy = hit->GetY() - y;
- double dt = hit->GetTime() - t;
- double sx = hit->GetDx();
- double sy = hit->GetDy();
- double st = hit->GetTimeError();
-
- // residuals
- fhResidualX.Fill(dx);
- fhResidualY.Fill(dy);
- fhResidualT.Fill(dt);
-
- // pulls
- if ((sx < 1.e-5) || (sy < 1.e-5) || (st < 1.e-5)) {
- LOG(error) << "Tof hit errors are not properly set: errX " << hit->GetDx() << " errY " << hit->GetDy() << " errT "
- << st;
- // return;
- }
-
- fhPullX.Fill(dx / sx);
- fhPullY.Fill(dy / sy);
- fhPullT.Fill(dt / st);
-
- } // iHit
-
- // fill efficiency: N hits per MC point
-
- for (int iE = 0; iE < nMcEvents; iE++) {
- int fileId = fMcEventList->GetFileIdByIndex(iE);
- int eventId = fMcEventList->GetEventIdByIndex(iE);
- int nPointsEv = selectedPoints[iE].size();
- for (int ip = 0; ip < nPointsEv; ip++) {
-
- int iPoint = selectedPoints[iE][ip];
- CbmTofPoint* p = dynamic_cast<CbmTofPoint*>(fMcPoints->Get(fileId, eventId, iPoint));
- if (p == nullptr) {
- LOG(error) << "MC point doesn't exist";
- return;
- }
-
- int iSt = GetStationIndex(p);
- if (iSt < 0) continue;
-
- int StationIndex = iSt; //CbmTofTrackingInterface::Instance()->GetTrackingStationIndex(p);
- if (StationIndex < 0 || StationIndex >= fNtrackingStations) {
- LOG(fatal) << "Tof point layer id " << StationIndex << " is out of range, z of point = " << p->GetZ();
- return;
- }
- fhHitsPerPoint[StationIndex].Fill(pointNhits[iE][iPoint]);
- fhEfficiencyR[StationIndex].Fill(sqrt(p->GetX() * p->GetX() + p->GetY() * p->GetY()),
- (pointNhits[iE][iPoint] > 0));
- fhEfficiencyXY[StationIndex].Fill(p->GetX(), p->GetY(), (pointNhits[iE][iPoint] > 0));
- }
- }
-}
-
-
-// -------------------------------------------------------------------------
-TFolder& CbmTrackerInputQaTof::GetQa()
-{
- gStyle->SetPaperSize(20, 20);
-
- for (Int_t i = 0; i < fNtrackingStations; i++) {
- fCanvHitsPerPoint.cd(i + 1);
- fhHitsPerPoint[i].DrawCopy("", "");
- fCanvPointsPerHit.cd(i + 1);
- fhPointsPerHit[i].DrawCopy("", "");
-
- fCanvEfficiencyR.cd(i + 1);
- fhEfficiencyR[i].DrawCopy("colz", "");
-
- fCanvEfficiencyXY.cd(i + 1);
- fhEfficiencyXY[i].DrawCopy("colz", "");
- }
-
- for (UInt_t i = 0; i < fHistList.size(); i++) {
- cbm::qa::util::FitKaniadakisGaussian(fHistList[i]);
- }
-
- for (UInt_t i = 0; i < 3; i++) {
- fCanvResidual.cd(i + 1);
- fHistList[i]->DrawCopy("", "");
- fCanvPull.cd(i + 1);
- fHistList[i + 3]->DrawCopy("", "");
- }
-
- return fOutFolder;
-}
diff --git a/reco/L1/qa/CbmTrackerInputQaTof.h b/reco/L1/qa/CbmTrackerInputQaTof.h
deleted file mode 100644
index c6d18aae529badea2ebfe59a0e78d1935869e069..0000000000000000000000000000000000000000
--- a/reco/L1/qa/CbmTrackerInputQaTof.h
+++ /dev/null
@@ -1,162 +0,0 @@
-/* Copyright (C) 2021 Facility for Antiproton and Ion Research in Europe, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov[committer]*/
-
-/// @file CbmTrackerInputQaTof.h
-/// @author Sergey Gorbunov
-/// @date 02.11.2021
-
-
-#pragma once
-
-#include "CbmQaCanvas.h"
-#include "CbmQaHist.h"
-
-#include <FairTask.h>
-
-#include <TFolder.h>
-#include <TH1D.h>
-#include <TH1F.h>
-#include <TH1I.h>
-#include <TH2F.h>
-#include <TParameter.h>
-#include <TProfile.h>
-#include <TProfile2D.h>
-
-class CbmMCDataManager;
-class CbmMCEventList;
-class CbmMCDataArray;
-class CbmTimeSlice;
-class CbmTofDigiPar;
-class CbmTofParSetGeo;
-class CbmTofPoint;
-
-class TClonesArray;
-
-///
-/// CbmTrackerInputQaTof class represents all the CA tracker requirements for the Tof detector.
-/// When this QA test is passed, the tracker must work (at least that is the idea).
-///
-/// The class ensures that the tracker has the correct understanding of the Tof geometry and interfaces
-/// and checks the quality of all tracking-related aspects of the Tof data.
-///
-class CbmTrackerInputQaTof : public FairTask {
-
- public:
- /// Constructor
- CbmTrackerInputQaTof(const char* name = "TrackerInputQaTof", Int_t verbose = 1);
-
- /// Destructor
- ~CbmTrackerInputQaTof();
-
- /// FairTask: Intialisation at begin of run.
- InitStatus Init() override { return ReInit(); }
-
- /// FairTask: Reinitialisation.
- InitStatus ReInit() override;
-
- /// FairTask: Intialise parameter containers.
- void SetParContainers() override;
-
- /// FairTask: Action at end of run. For this task and all of the subtasks.
- void Finish() override;
-
- /// TTask: Clear all data structures created by a previous execution of a task.
- void Clear(Option_t* /*option*/ = "") override {}
-
- /// TTask: Process a timeslice
- void Exec(Option_t*) override;
-
- /// Prepare the Qa output and return it as a reference to an internal TFolder.
- /// The reference is non-const, because FairSink can not write const objects
- TFolder& GetQa();
-
- private:
- /// Check the geometry
- InitStatus GeometryQa();
-
- /// Analysis of hit uncertainty (pull) distributions
- void ResolutionQa();
-
- /// Free the memory etc.
- void DeInit();
-
- int GetStationIndex(CbmTofPoint* p);
-
-
- // Setup
-
- Bool_t fIsTofInSetup{false};
- Bool_t fIsMcPresent{false};
-
- Int_t fNtrackingStations{0};
-
- CbmTimeSlice* fTimeSlice{nullptr};
- CbmTofParSetGeo* fTofGeoPar{nullptr};
- CbmTofDigiPar* fTofDigiPar{nullptr};
-
- /// MC data
- CbmMCEventList* fMcEventList{nullptr}; // list of MC events connected to the current time slice
- CbmMCDataManager* fMcManager{nullptr};
- CbmMCDataArray* fMcTracks{nullptr};
- CbmMCDataArray* fMcPoints{nullptr};
-
- /// Data
- TClonesArray* fClusters{nullptr};
- TClonesArray* fHits{nullptr};
- TClonesArray* fHitMatches{nullptr};
-
- TClonesArray* fHitDigiMatches{nullptr};
- TClonesArray* fDigis{nullptr};
- TClonesArray* fTofCalDigiMatch{nullptr};
-
- /// Output
-
-
- TFolder fOutFolder{"TrackerInputQaTof", "TrackerInputQaTof"}; /// output folder with histos and canvases
- TFolder* fHistFolder{nullptr}; /// subfolder for histograms
-
- TParameter<int> fNevents{"nEvents", 0}; /// number of processed events
-
- /// Histogram for Residual Distribution
- CbmQaHist<TH1D> fhResidualX{"hResidualX", "Tof: Residual X;(X_{reco} - X_{MC})(cm)", 100, -10, 10};
- CbmQaHist<TH1D> fhResidualY{"hRresidualY", "Tof: Residual Y;(Y_{reco} - Y_{MC})(cm)", 100, -10, 10};
- CbmQaHist<TH1D> fhResidualT{"hRresidualT", "Tof: Residual T;(T_{reco} - T_{MC})(ns)", 100, -50, 50};
- CbmQaHist<TH1D> fhResidualZ{"hRresidualZ", "Tof: Residual Z;(Z_{reco} - Z_{MC})(cm)", 100, -100, 100};
-
- /// Histogram for PULL Distribution
- CbmQaHist<TH1D> fhPullX{"hPullX", "Tof: Pull X;(X_{reco} - X_{MC}) / #sigmaU_{reco}", 100, -5, 5};
- CbmQaHist<TH1D> fhPullY{"hPullX", "Tof: Pull Y;(Y_{reco} - Y_{MC}) / #sigmaV_{reco}", 100, -5, 5};
- CbmQaHist<TH1D> fhPullT{"hPullT", "Tof: Pull T;(T_{reco} - T_{MC}) / #sigmaT_{reco}", 100, -5, 5};
-
- /// List of the above histogramms
- std::vector<CbmQaHist<TH1D>*> fHistList;
-
- /// hits purity
- std::vector<CbmQaHist<TH1I>> fhPointsPerHit;
-
- /// hits efficiency
- std::vector<CbmQaHist<TH1I>> fhHitsPerPoint;
-
- /// hits efficiency
- std::vector<CbmQaHist<TProfile2D>> fhEfficiencyXY;
- std::vector<CbmQaHist<TProfile>> fhEfficiencyR;
-
- /// Canvaces: collection of histogramms
-
- CbmQaCanvas fCanvResidual{"cResidual", "Residual Distribution", 3 * 500, 2 * 500};
- CbmQaCanvas fCanvPull{"cPull", "Pull Distribution", 3 * 500, 2 * 500};
- CbmQaCanvas fCanvEfficiencyXY{"cEfficiencyXY", "Efficiency XY: % reconstructed McPoint", 2 * 500, 2 * 500};
- CbmQaCanvas fCanvEfficiencyR{"cEfficiencyR", "Efficiency R: % reconstructed McPoint", 2 * 500, 2 * 500};
- CbmQaCanvas fCanvHitsPerPoint{"cHitsPerMcPoint", "Efficiency: Hits Per McPoint", 2 * 500, 2 * 500};
- CbmQaCanvas fCanvPointsPerHit{"cMcPointsPerHit", "Purity: McPoints per Hit", 2 * 500, 2 * 500};
-
- private:
- /// Suppressed copy constructor
- CbmTrackerInputQaTof(const CbmTrackerInputQaTof&) = delete;
-
- /// Suppressed assignment operator
- CbmTrackerInputQaTof& operator=(const CbmTrackerInputQaTof&) = delete;
-
- ClassDefOverride(CbmTrackerInputQaTof, 0);
-};
diff --git a/reco/L1/qa/CbmTrackerInputQaTrd.cxx b/reco/L1/qa/CbmTrackerInputQaTrd.cxx
deleted file mode 100644
index be4464037d7ae91e7d28f880d1667981f9db5107..0000000000000000000000000000000000000000
--- a/reco/L1/qa/CbmTrackerInputQaTrd.cxx
+++ /dev/null
@@ -1,757 +0,0 @@
-/* Copyright (C) 2021 Facility for Antiproton and Ion Research in Europe, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov[committer]*/
-
-/// @file CbmTrackerInputQaTrd.cxx
-/// @author Sergey Gorbunov
-/// @date 02.11.2021
-
-#include "CbmTrackerInputQaTrd.h"
-
-#include "CbmDefs.h"
-#include "CbmDigiManager.h"
-#include "CbmLink.h"
-#include "CbmMCDataArray.h"
-#include "CbmMCDataManager.h"
-#include "CbmMCEventList.h"
-#include "CbmMCTrack.h"
-#include "CbmMatch.h"
-#include "CbmQaCanvas.h"
-//#include "CbmSetup.h"
-#include "CbmQaUtil.h"
-#include "CbmTimeSlice.h"
-#include "CbmTrdCluster.h"
-#include "CbmTrdDigi.h"
-#include "CbmTrdHit.h"
-#include "CbmTrdParModDigi.h" // for CbmTrdModule
-#include "CbmTrdParModGeo.h"
-#include "CbmTrdParSetDigi.h" // for CbmTrdParSetDigi
-#include "CbmTrdParSetGeo.h"
-#include "CbmTrdPoint.h"
-
-#include <FairRootManager.h>
-#include <FairRunAna.h>
-#include <FairRuntimeDb.h>
-#include <FairSink.h>
-#include <FairTask.h>
-#include <Logger.h>
-
-#include <TClonesArray.h>
-#include <TDatabasePDG.h>
-#include <TGenericClassInfo.h>
-#include <TGeoManager.h>
-#include <TGeoNode.h>
-#include <TMathBase.h>
-#include <TObjArray.h>
-#include <TParameter.h>
-#include <TParticlePDG.h>
-#include <TString.h>
-#include <TStyle.h>
-
-#include <algorithm>
-#include <cstdio>
-#include <iostream>
-#include <map>
-#include <vector>
-
-ClassImp(CbmTrackerInputQaTrd);
-
-// -------------------------------------------------------------------------
-
-CbmTrackerInputQaTrd::CbmTrackerInputQaTrd(const char* name, Int_t verbose) : FairTask(name, verbose)
-{
- // Constructor
-
- // Create a list of histogramms
-
- fHistList.clear();
- fHistList.reserve(20);
- fHistList.push_back(&fh1DresidualU);
- fHistList.push_back(&fh1DresidualV);
- fHistList.push_back(&fh1DresidualT);
- fHistList.push_back(&fh2DresidualX);
- fHistList.push_back(&fh2DresidualY);
- fHistList.push_back(&fh2DresidualT);
- fHistList.push_back(&fh1DpullU);
- fHistList.push_back(&fh1DpullV);
- fHistList.push_back(&fh1DpullT);
- fHistList.push_back(&fh2DpullX);
- fHistList.push_back(&fh2DpullY);
- fHistList.push_back(&fh2DpullT);
-
- // Keep the ownership on the histograms to avoid double deletion
- for (unsigned int i = 0; i < fHistList.size(); i++) {
- fHistList[i]->SetDirectory(nullptr);
- }
-}
-
-// -------------------------------------------------------------------------
-CbmTrackerInputQaTrd::~CbmTrackerInputQaTrd()
-{ /// Destructor
- DeInit();
-}
-
-
-// -------------------------------------------------------------------------
-void CbmTrackerInputQaTrd::DeInit()
-{
- fIsTrdInSetup = 0;
- fIsMcPresent = false;
- fNtrackingStations = 0;
-
- fTimeSlice = nullptr;
- fDigiManager = nullptr;
-
- fMcManager = nullptr;
- fMcTracks = nullptr;
- fMcPoints = nullptr;
-
- fClusters = nullptr;
- fHits = nullptr;
- fHitMatches = nullptr;
-
- fOutFolder.Clear();
-
- fHistFolder = nullptr;
-
- fNevents.SetVal(0);
-
- fhPointsPerHit.clear();
- fhHitsPerPoint.clear();
-}
-// -------------------------------------------------------------------------
-
-void CbmTrackerInputQaTrd::SetParContainers()
-{
- // FIXME: SZh 13.01.2023: With tracking detector interfaces these data bases are no longer needed
- // to be initialized in this class
- fTrdDigiPar = nullptr;
- fTrdGeoPar = nullptr;
-
- // Get run and runtime database
- FairRunAna* ana = FairRunAna::Instance();
- if (!ana) {
- LOG(fatal) << "No FairRunAna instance";
- }
-
- FairRuntimeDb* rtdb = ana->GetRuntimeDb();
- if (!rtdb) {
- LOG(fatal) << "No FairRuntimeDb instance";
- }
-
- fTrdDigiPar = dynamic_cast<CbmTrdParSetDigi*>(rtdb->getContainer("CbmTrdParSetDigi"));
- fTrdGeoPar = dynamic_cast<CbmTrdParSetGeo*>(rtdb->getContainer("CbmTrdParSetGeo"));
-}
-
-// -------------------------------------------------------------------------
-InitStatus CbmTrackerInputQaTrd::ReInit()
-{
- // Initialize and check the setup
-
- DeInit();
-
- // Check if TRD is present in the CBM setup
- // A straightforward way to do so is CbmSetup::Instance()->IsActive(ECbmModuleId::kTrd),
- // but unfortunatelly CbmSetup class is unaccesible.
- // ( CbmSimSteer library requires libfreetype that has linking problems on MacOsX )
- // Therefore let's simply check if any TRD material is present in the geometry.
- //
- fIsTrdInSetup = 0;
- {
- TObjArray* topNodes = gGeoManager->GetTopNode()->GetNodes();
- for (Int_t iTopNode = 0; iTopNode < topNodes->GetEntriesFast(); iTopNode++) {
- TGeoNode* topNode = static_cast<TGeoNode*>(topNodes->At(iTopNode));
- if (TString(topNode->GetName()).Contains("trd")) {
- fIsTrdInSetup = 1;
- break;
- }
- }
- }
-
- if (!fIsTrdInSetup) {
- LOG(info) << "TRD is not in the setup, do nothing";
- return kSUCCESS;
- }
-
- // check parameter containers
-
- if (!fTrdDigiPar) {
- LOG(error) << "No CbmTrdParSetDigi container in FairRuntimeDb";
- return kERROR;
- }
-
- if (!fTrdGeoPar) {
- LOG(error) << "No CbmTrdParSetGeo container in FairRuntimeDb";
- return kERROR;
- }
-
- FairRootManager* manager = FairRootManager::Instance();
- fDigiManager = CbmDigiManager::Instance();
- fDigiManager->Init();
-
- fTimeSlice = dynamic_cast<CbmTimeSlice*>(manager->GetObject("TimeSlice."));
- if (!fTimeSlice) {
- LOG(error) << "No time slice found";
- return kERROR;
- }
-
- fHits = dynamic_cast<TClonesArray*>(manager->GetObject("TrdHit"));
-
- if (!fHits) {
- LOG(error) << "No TRD hit array found";
- return kERROR;
- }
-
- fClusters = dynamic_cast<TClonesArray*>(manager->GetObject("TrdCluster"));
-
- if (!fClusters) {
- LOG(error) << "No TRD cluster array found";
- return kERROR;
- }
-
- fMcManager = dynamic_cast<CbmMCDataManager*>(manager->GetObject("MCDataManager"));
-
- fIsMcPresent = (fMcManager != nullptr);
-
- if (fIsMcPresent) {
-
- fMcEventList = dynamic_cast<CbmMCEventList*>(manager->GetObject("MCEventList."));
- fMcTracks = fMcManager->InitBranch("MCTrack");
- fMcPoints = fMcManager->InitBranch("TrdPoint");
- fHitMatches = dynamic_cast<TClonesArray*>(manager->GetObject("TrdHitMatch"));
-
- // we assume that when TRD is in the setup and the MC manager is present,
- // the TRD MC data must be present too
-
- if (!fMcEventList) {
- LOG(error) << ": No MCEventList data!";
- return kERROR;
- }
-
- if (!fMcTracks) {
- LOG(error) << "No MC tracks found";
- return kERROR;
- }
-
- if (!fMcPoints) {
- LOG(error) << "No MC points found";
- return kERROR;
- }
-
- if (!fHitMatches) {
- LOG(error) << "No TRD hit matches found";
- return kERROR;
- }
-
- if (!fDigiManager->IsMatchPresent(ECbmModuleId::kTrd)) {
- LOG(error) << "No TRD digi matches found";
- return kERROR;
- }
- }
-
- // count TRD stations
- // TODO: put the code below to some TRD geometry interface
-
- fNtrackingStations = 0;
- {
- Int_t layerCounter = 0;
- TObjArray* topNodes = gGeoManager->GetTopNode()->GetNodes();
- for (Int_t iTopNode = 0; iTopNode < topNodes->GetEntriesFast(); iTopNode++) {
- TGeoNode* topNode = static_cast<TGeoNode*>(topNodes->At(iTopNode));
- if (TString(topNode->GetName()).Contains("trd")) {
- TObjArray* layers = topNode->GetNodes();
- for (Int_t iLayer = 0; iLayer < layers->GetEntriesFast(); iLayer++) {
- TGeoNode* layer = static_cast<TGeoNode*>(layers->At(iLayer));
- if (TString(layer->GetName()).Contains("layer")) {
- layerCounter++;
- }
- }
- }
- }
- fNtrackingStations = layerCounter;
- }
-
- LOG(debug) << "Init(): fNtrackingStations = " << fNtrackingStations;
-
- if (fNtrackingStations <= 0) {
- LOG(error) << "can not count TRD tracking stations";
- return kERROR;
- }
-
- // initialise histograms
- fOutFolder.SetOwner(false);
- fHistFolder = fOutFolder.AddFolder("rawHist", "Raw Histograms");
- gStyle->SetOptStat(0);
-
- fNevents.SetVal(0);
- fHistFolder->Add(&fNevents);
-
- for (unsigned int i = 0; i < fHistList.size(); i++) {
- fHistList[i]->Reset();
- fHistFolder->Add(fHistList[i]);
- }
-
- fhPointsPerHit.clear();
- fhHitsPerPoint.clear();
- fhEfficiencyR.clear();
- fhEfficiencyXY.clear();
-
- fhPointsPerHit.reserve(fNtrackingStations);
- fhHitsPerPoint.reserve(fNtrackingStations);
- fhEfficiencyR.reserve(fNtrackingStations);
- fhEfficiencyXY.reserve(fNtrackingStations);
-
- for (Int_t i = 0; i < fNtrackingStations; i++) {
- fhPointsPerHit.emplace_back(Form("hMcPointsPerHit%i", i),
- Form("MC Points per Hit: Station %i;N mc points / hit", i), 10, -0.5, 9.5);
- fhPointsPerHit[i].SetDirectory(nullptr);
- fhPointsPerHit[i].SetLineWidth(2);
- fhPointsPerHit[i].SetOptStat(110);
- fHistFolder->Add(&fhPointsPerHit[i]);
- }
-
- for (Int_t i = 0; i < fNtrackingStations; i++) {
- fhHitsPerPoint.emplace_back(Form("hHitsPerMcPoint%i", i),
- Form("Hits per MC Point: Station %i; N hits / mc point", i), 10, -0.5, 9.5);
- fhHitsPerPoint[i].SetDirectory(nullptr);
- fhHitsPerPoint[i].SetLineWidth(2);
- fhHitsPerPoint[i].SetOptStat(110);
- fHistFolder->Add(&fhHitsPerPoint[i]);
- }
-
- for (Int_t i = 0; i < fNtrackingStations; i++) {
-
- double dx = 350;
- double dy = 350;
- double dr = sqrt(dx * dx + dy * dy);
-
- fhEfficiencyR.emplace_back(Form("hEfficiencyR%i", i), Form("Efficiency R: Station %i;R [cm]", i), 100, 0, dr);
- fhEfficiencyR[i].SetDirectory(nullptr);
- fhEfficiencyR[i].SetOptStat(1110);
- fHistFolder->Add(&fhEfficiencyR[i]);
-
- fhEfficiencyXY.emplace_back(Form("hEfficiencyXY%i", i), Form("Efficiency XY: Station %i;X [cm];Y [cm]", i), 50, -dx,
- dx, 50, -dy, dy);
- fhEfficiencyXY[i].SetDirectory(nullptr);
- fhEfficiencyXY[i].SetOptStat(10);
- fhEfficiencyXY[i].GetYaxis()->SetTitleOffset(1.4);
- fHistFolder->Add(&fhEfficiencyXY[i]);
- }
-
- fCanvResidual.Clear();
- fCanvResidual.Divide2D(6);
-
- fCanvPull.Clear();
- fCanvPull.Divide2D(6);
-
- fCanvEfficiencyR.Clear();
- fCanvEfficiencyR.Divide2D(fNtrackingStations);
-
- fCanvEfficiencyXY.Clear();
- fCanvEfficiencyXY.Divide2D(fNtrackingStations);
-
- fCanvPointsPerHit.Clear();
- fCanvPointsPerHit.Divide2D(fNtrackingStations);
-
- fCanvHitsPerPoint.Clear();
- fCanvHitsPerPoint.Divide2D(fNtrackingStations);
-
- fOutFolder.Add(&fCanvResidual);
- fOutFolder.Add(&fCanvPull);
- fOutFolder.Add(&fCanvEfficiencyR);
- fOutFolder.Add(&fCanvEfficiencyXY);
- fOutFolder.Add(&fCanvPointsPerHit);
- fOutFolder.Add(&fCanvHitsPerPoint);
-
- // analyse the geometry setup
- return GeometryQa();
-}
-
-
-// -------------------------------------------------------------------------
-void CbmTrackerInputQaTrd::Exec(Option_t*)
-{
- if (!fIsTrdInSetup) {
- return;
- }
-
- // update number of processed events
- fNevents.SetVal(fNevents.GetVal() + 1);
- LOG(debug) << "Event: " << fNevents.GetVal();
- ResolutionQa();
-}
-
-
-// -------------------------------------------------------------------------
-void CbmTrackerInputQaTrd::Finish()
-{
- FairSink* sink = FairRootManager::Instance()->GetSink();
- if (sink) {
- sink->WriteObject(&GetQa(), nullptr);
- }
-}
-
-
-// -------------------------------------------------------------------------
-InitStatus CbmTrackerInputQaTrd::GeometryQa()
-{
- // check geometry of the TRD modules
-
- double lastZ = -1;
- for (int iStation = 0; iStation < fNtrackingStations; iStation++) {
-
- int ModuleId = fTrdDigiPar->GetModuleId(iStation);
-
- const CbmTrdParModGeo* pGeo = dynamic_cast<const CbmTrdParModGeo*>(fTrdGeoPar->GetModulePar(ModuleId));
- if (!pGeo) {
- LOG(fatal) << " No Geo params for station " << iStation << " module " << ModuleId << " found ";
- return kERROR;
- }
-
- double staZ = pGeo->GetZ();
-
- // check that the stations are properly ordered in Z
- if ((iStation > 0) && (staZ <= lastZ)) {
- LOG(error) << "TRD trackig stations are not properly ordered in Z";
- return kERROR;
- }
- lastZ = staZ;
-
- //TODO: what are these 3 and 6 here in L1 code? Why are they hardcoed?
- // if (num == 0 || num == 2 || num == 4) geo.push_back(3);
- // if (num == 1 || num == 3) geo.push_back(6);
- }
-
- return kSUCCESS;
-}
-
-// -------------------------------------------------------------------------
-void CbmTrackerInputQaTrd::ResolutionQa()
-{
- if (!fDigiManager->IsMatchPresent(ECbmModuleId::kTrd)) return;
-
- if (!fIsMcPresent) return;
-
- Int_t nHits = fHits->GetEntriesFast();
- Int_t nClusters = fClusters->GetEntriesFast();
- Int_t nDigis = fDigiManager->GetNofDigis(ECbmModuleId::kTrd);
-
- int nMcEvents = fMcEventList->GetNofEvents();
-
- // Vector of integers parallel to mc points
-
- std::vector<std::vector<int>> pointNhits;
- pointNhits.resize(nMcEvents);
- for (int iE = 0; iE < nMcEvents; iE++) {
- int fileId = fMcEventList->GetFileIdByIndex(iE);
- int eventId = fMcEventList->GetEventIdByIndex(iE);
- int nPoints = fMcPoints->Size(fileId, eventId);
- pointNhits[iE].resize(nPoints, 0);
- }
-
- for (Int_t iHit = 0; iHit < nHits; iHit++) {
-
- CbmTrdHit* hit = dynamic_cast<CbmTrdHit*>(fHits->At(iHit));
- if (!hit) {
- LOG(error) << "TRD hit N " << iHit << " doesn't exist";
- return;
- }
-
- if ((int) hit->GetClassType() > 1) {
- // class type 0: trd-1D
- // class type 1: trd-2D
- // the return type is (currently) boolean, so this error should never happen
- LOG(error) << "Unknown detector class type " << hit->GetClassType() << " for TRD hit N " << iHit;
- return;
- }
-
- const int address = hit->GetAddress();
- int StationIndex = CbmTrdAddress::GetLayerId(address);
-
- if (StationIndex < 0 || StationIndex >= fNtrackingStations) {
- LOG(fatal) << "TRD hit layer id " << StationIndex << " is out of range";
- return;
- }
-
- Int_t clusterId = hit->GetRefId();
- if (clusterId < 0 || clusterId >= nClusters) {
- LOG(error) << "TRD cluster id " << clusterId << " is out of range";
- return;
- }
-
- CbmTrdCluster* cluster = dynamic_cast<CbmTrdCluster*>(fClusters->At(clusterId));
- if (!cluster) {
- LOG(error) << "TRD cluster N " << clusterId << " doesn't exist";
- return;
- }
-
- if (cluster->GetAddress() != address) {
- LOG(error) << "TRD hit address " << address << " differs from its cluster address " << cluster->GetAddress();
- return;
- }
-
- Int_t nClusterDigis = cluster->GetNofDigis();
-
- if (nClusterDigis <= 0) {
- LOG(error) << "hit match for TRD cluster " << clusterId << " has no digis ";
- return;
- }
-
- // custom finder of the digi matches
-
- CbmMatch myHitMatch;
- for (Int_t iDigi = 0; iDigi < nClusterDigis; iDigi++) {
- Int_t digiIdx = cluster->GetDigi(iDigi);
- if (digiIdx < 0 || digiIdx >= nDigis) {
- LOG(error) << "TRD cluster: digi index " << digiIdx << " out of range ";
- return;
- }
- const CbmTrdDigi* digi = CbmDigiManager::Instance()->Get<CbmTrdDigi>(digiIdx);
- if (!digi) {
- LOG(fatal) << "TRD digi " << iDigi << " not found";
- return;
- }
-
- if (digi->GetAddressModule() != address) {
- LOG(error) << "TRD hit address " << address << " differs from its digi address " << digi->GetAddressModule();
- return;
- }
- const CbmMatch* match = dynamic_cast<const CbmMatch*>(fDigiManager->GetMatch(ECbmModuleId::kTrd, digiIdx));
- if (!match) {
- LOG(fatal) << "TRD digi match " << digiIdx << " not found";
- return;
- }
- myHitMatch.AddLinks(*match);
- }
- if (myHitMatch.GetNofLinks() == 0) {
- continue;
- }
-
- const CbmLink& bestLink = myHitMatch.GetMatchedLink();
-
- { // check if the hit match is correct
- CbmMatch* hitMatch = dynamic_cast<CbmMatch*>(fHitMatches->At(iHit));
- if (hitMatch == nullptr) {
- LOG(error) << "hit match for TRD hit " << iHit << " doesn't exist";
- }
- else {
- const CbmLink& link = hitMatch->GetMatchedLink();
- if ((link != bestLink) && (link.GetWeight() != bestLink.GetWeight())) {
- LOG(error) << "hit match for TRD hit " << iHit << " doesn't correspond to digi matches";
- }
- }
- }
-
- // how many MC points? fill N hits per mc point
-
- int nHitPoints = 0;
- for (int i = 0; i < myHitMatch.GetNofLinks(); i++) {
- const CbmLink& link = myHitMatch.GetLink(i);
- if (link.GetIndex() >= 0) { // not a noise digi
- nHitPoints++;
- int iE = fMcEventList->GetEventIndex(link);
- if (iE < 0 || iE >= nMcEvents) {
- LOG(error) << "link points to a non-existing MC event";
- return;
- }
- if (link.GetIndex() >= (int) pointNhits[iE].size()) {
- LOG(error) << "link points to a non-existing MC point";
- return;
- }
- pointNhits[iE][link.GetIndex()]++;
- }
- }
-
- fhPointsPerHit[StationIndex].Fill(nHitPoints);
-
- // take corresponding MC point
-
- // skip hits from the noise digis
- if (bestLink.GetIndex() < 0) {
- continue;
- }
-
- CbmTrdPoint* p = dynamic_cast<CbmTrdPoint*>(fMcPoints->Get(bestLink));
- if (p == nullptr) {
- LOG(error) << "link points to a non-existing MC point";
- return;
- }
-
- if (p->GetModuleAddress() != (int) CbmTrdAddress::GetModuleAddress(address)) {
- LOG(error) << "mc point module address differs from the hit module address";
- return;
- }
-
- // check that the nominal station Z is not far from the active material
- {
- // the cut of 1 cm is choosen arbitrary and can be changed
-
- int ModuleId = fTrdDigiPar->GetModuleId(StationIndex);
-
- const CbmTrdParModGeo* pGeo = dynamic_cast<const CbmTrdParModGeo*>(fTrdGeoPar->GetModulePar(ModuleId));
- if (!pGeo) {
- LOG(fatal) << " No Geo params for station " << StationIndex << " module " << ModuleId << " found ";
- return;
- }
-
- double staZ = pGeo->GetZ(); // module->GetZ(); //+ 410;
- if ((staZ < p->GetZIn() - 1.) || (staZ > p->GetZOut() + 1.)) {
- LOG(error) << "TRD station " << StationIndex << ": active material Z[" << p->GetZIn() << " cm," << p->GetZOut()
- << " cm] is too far from the nominal station Z " << staZ << " cm";
- return;
- }
- // the cut of 1 cm is choosen arbitrary and can be changed
- if (fabs(hit->GetZ() - staZ) > 1.) {
- LOG(error) << "TRD station " << StationIndex << ": hit Z " << hit->GetZ()
- << " cm, is too far from the nominal station Z " << staZ << " cm";
- return;
- }
- }
-
- // residual and pull
-
- if (nHitPoints != 1) continue; // only check residual for non-mixed hits
-
- double t0 = fMcEventList->GetEventTime(bestLink);
- if (t0 < 0) {
- LOG(error) << "MC event time doesn't exist for a TRD link";
- return;
- }
-
- double x = p->GetX(); // take the "In" point since the time is stored only for this point
- double y = p->GetY();
- double z = p->GetZ();
- double t = t0 + p->GetTime();
- double px = p->GetPx();
- double py = p->GetPy();
- double pz = p->GetPz();
-
- // skip very slow particles
- if (fabs(p->GetPzOut()) < 0.01) continue;
-
- // transport the particle from the MC point to the hit Z
- double dz = hit->GetZ() - z;
- x += dz * px / pz;
- y += dz * py / pz;
-
- // get particle mass
- double mass = 0;
- {
- CbmLink mcTrackLink = bestLink;
- mcTrackLink.SetIndex(p->GetTrackID());
- CbmMCTrack* mcTrack = dynamic_cast<CbmMCTrack*>(fMcTracks->Get(mcTrackLink));
- if (!mcTrack) {
- LOG(error) << "MC track " << p->GetTrackID() << " doesn't exist";
- return;
- }
- Int_t pdg = mcTrack->GetPdgCode();
- if (pdg < 9999999 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) {
- mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
- }
- }
-
- constexpr double speedOfLight = 29.979246; // cm/ns
- TVector3 mom3;
- p->Momentum(mom3);
- t += dz / (pz * speedOfLight) * sqrt(mass * mass + mom3.Mag2());
-
- double du = hit->GetX() - x;
- double dv = hit->GetY() - y;
- double dt = hit->GetTime() - t;
- double su = hit->GetDx();
- double sv = hit->GetDy();
- double st = hit->GetTimeError();
-
- // residuals
- if (hit->GetClassType() == 0) {
- if (StationIndex % 2) {
- std::swap(du, dv);
- std::swap(su, sv);
- }
- fh1DresidualU.Fill(du);
- fh1DresidualV.Fill(dv);
- fh1DresidualT.Fill(hit->GetTime() - t);
- }
- else {
- fh2DresidualX.Fill(du);
- fh2DresidualY.Fill(dv);
- fh2DresidualT.Fill(hit->GetTime() - t);
- }
-
- // pulls
- if ((su < 1.e-5) || (sv < 1.e-5) || (st < 1.e-5)) {
- LOG(error) << "TRD hit errors are not properly set: errX " << hit->GetDx() << " errY " << hit->GetDy() << " errT "
- << st;
- return;
- }
-
- if (hit->GetClassType() == 0) {
- fh1DpullU.Fill(du / su);
- fh1DpullV.Fill(dv / sv);
- fh1DpullT.Fill(dt / st);
- }
- else {
- fh2DpullX.Fill(du / su);
- fh2DpullY.Fill(dv / sv);
- fh2DpullT.Fill(dt / st);
- }
-
- } // iHit
-
- // fill efficiency: N hits per MC point
-
- for (int iE = 0; iE < nMcEvents; iE++) {
- int fileId = fMcEventList->GetFileIdByIndex(iE);
- int eventId = fMcEventList->GetEventIdByIndex(iE);
- int nPoints = fMcPoints->Size(fileId, eventId);
- for (int ip = 0; ip < nPoints; ip++) {
- CbmTrdPoint* p = dynamic_cast<CbmTrdPoint*>(fMcPoints->Get(fileId, eventId, ip));
- if (p == nullptr) {
- LOG(error) << "MC point doesn't exist";
- return;
- }
- auto address = p->GetModuleAddress();
- int StationIndex = CbmTrdAddress::GetLayerId(address);
- if (StationIndex < 0 || StationIndex >= fNtrackingStations) {
- LOG(fatal) << "TRD hit layer id " << StationIndex << " for module " << address << " is out of range";
- return;
- }
- fhHitsPerPoint[StationIndex].Fill(pointNhits[iE][ip]);
- fhEfficiencyR[StationIndex].Fill(sqrt(p->GetX() * p->GetX() + p->GetY() * p->GetY()), (pointNhits[iE][ip] > 0));
- fhEfficiencyXY[StationIndex].Fill(p->GetX(), p->GetY(), (pointNhits[iE][ip] > 0));
- }
- }
-}
-
-
-// -------------------------------------------------------------------------
-TFolder& CbmTrackerInputQaTrd::GetQa()
-{
- gStyle->SetPaperSize(20, 20);
-
- for (Int_t i = 0; i < fNtrackingStations; i++) {
- fCanvHitsPerPoint.cd(i + 1);
- fhHitsPerPoint[i].DrawCopy("", "");
- fCanvPointsPerHit.cd(i + 1);
- fhPointsPerHit[i].DrawCopy("", "");
-
- fCanvEfficiencyR.cd(i + 1);
- fhEfficiencyR[i].DrawCopy("colz", "");
-
- fCanvEfficiencyXY.cd(i + 1);
- fhEfficiencyXY[i].DrawCopy("colz", "");
- }
-
- for (UInt_t i = 0; i < fHistList.size(); i++) {
- cbm::qa::util::FitKaniadakisGaussian(fHistList[i]);
- }
-
- for (UInt_t i = 0; i < 6; i++) {
- fCanvResidual.cd(i + 1);
- fHistList[i]->DrawCopy("", "");
- fCanvPull.cd(i + 1);
- fHistList[6 + i]->DrawCopy("", "");
- }
-
- return fOutFolder;
-}
diff --git a/reco/L1/qa/CbmTrackerInputQaTrd.h b/reco/L1/qa/CbmTrackerInputQaTrd.h
deleted file mode 100644
index 14fcb7f7248fc98f0f39e37a00d06d80ddcdf342..0000000000000000000000000000000000000000
--- a/reco/L1/qa/CbmTrackerInputQaTrd.h
+++ /dev/null
@@ -1,164 +0,0 @@
-/* Copyright (C) 2021 Facility for Antiproton and Ion Research in Europe, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov[committer]*/
-
-/// @file CbmTrackerInputQaTrd.h
-/// @author Sergey Gorbunov
-/// @date 02.11.2021
-
-
-#pragma once
-
-#include "CbmQaCanvas.h"
-#include "CbmQaHist.h"
-
-#include <FairTask.h>
-
-#include <TFolder.h>
-#include <TH1D.h>
-#include <TH1F.h>
-#include <TH1I.h>
-#include <TH2F.h>
-#include <TParameter.h>
-#include <TProfile.h>
-#include <TProfile2D.h>
-
-class CbmDigiManager;
-class CbmMCDataManager;
-class CbmMCEventList;
-class CbmMCDataArray;
-class CbmTimeSlice;
-class CbmTrdParSetGeo;
-class CbmTrdParSetDigi;
-
-class TClonesArray;
-
-///
-/// CbmTrackerInputQaTrd class represents all the CA tracker requirements for the TRD detector.
-/// When this QA test is passed, the tracker must work (at least that is the idea).
-///
-/// The class ensures that the tracker has the correct understanding of the TRD geometry and interfaces
-/// and checks the quality of all tracking-related aspects of the TRD data.
-///
-class CbmTrackerInputQaTrd : public FairTask {
-
- public:
- /// Constructor
- CbmTrackerInputQaTrd(const char* name = "TrackerInputQaTrd", Int_t verbose = 1);
-
- /// Destructor
- ~CbmTrackerInputQaTrd();
-
- /// FairTask: Intialisation at begin of run.
- InitStatus Init() override { return ReInit(); }
-
- /// FairTask: Reinitialisation.
- InitStatus ReInit() override;
-
- /// FairTask: Intialise parameter containers.
- void SetParContainers() override;
-
- /// FairTask: Action at end of run. For this task and all of the subtasks.
- void Finish() override;
-
- /// TTask: Clear all data structures created by a previous execution of a task.
- void Clear(Option_t* /*option*/ = "") override {}
-
- /// TTask: Process a timeslice
- void Exec(Option_t*) override;
-
- /// Prepare the Qa output and return it as a reference to an internal TFolder.
- /// The reference is non-const, because FairSink can not write const objects
- TFolder& GetQa();
-
- private:
- /// Check the geometry
- InitStatus GeometryQa();
-
- /// Analysis of hit uncertainty (pull) distributions
- void ResolutionQa();
-
- /// Free the memory etc.
- void DeInit();
-
-
- // Setup
-
- Bool_t fIsTrdInSetup{false};
- Bool_t fIsMcPresent{false};
-
- Int_t fNtrackingStations{0};
-
- CbmTimeSlice* fTimeSlice{nullptr};
- CbmTrdParSetGeo* fTrdGeoPar{nullptr};
- CbmTrdParSetDigi* fTrdDigiPar{nullptr};
- CbmDigiManager* fDigiManager{nullptr};
-
- /// MC data
- CbmMCEventList* fMcEventList{nullptr}; // list of MC events connected to the current time slice
- CbmMCDataManager* fMcManager{nullptr};
- CbmMCDataArray* fMcTracks{nullptr};
- CbmMCDataArray* fMcPoints{nullptr};
-
- /// Data
- TClonesArray* fClusters{nullptr};
- TClonesArray* fHits{nullptr};
- TClonesArray* fHitMatches{nullptr};
-
- /// Output
-
-
- TFolder fOutFolder{"TrackerInputQaTrd", "TrackerInputQaTrd"}; /// output folder with histos and canvases
- TFolder* fHistFolder{nullptr}; /// subfolder for histograms
-
- TParameter<int> fNevents{"nEvents", 0}; /// number of processed events
-
- /// Histogram for Residual Distribution
- CbmQaHist<TH1D> fh1DresidualU{"h1DresidualU", "Trd1D: Residual U;(U_{reco} - U_{MC})(cm)", 100, -.5, .5};
- CbmQaHist<TH1D> fh1DresidualV{"h1DresidualV", "Trd1D: Residual V;(V_{reco} - V_{MC})(cm)", 100, -10, 10};
- CbmQaHist<TH1D> fh1DresidualT{"h1DresidualT", "Trd1D: Residual T;(T_{reco} - T_{MC})(ns)", 100, -50, 50};
-
- CbmQaHist<TH1D> fh2DresidualX{"h2DresidualX", "Trd2D: Residual X;(X_{reco} - X_{MC})(cm)", 100, -5, 5};
- CbmQaHist<TH1D> fh2DresidualY{"h2DresidualY", "Trd2D: Residual Y;(Y_{reco} - Y_{MC})(cm)", 100, -5, 5};
- CbmQaHist<TH1D> fh2DresidualT{"h2DresidualT", "Trd2D: Residual T;(T_{reco} - T_{MC})(ns)", 100, -1000, 1000};
-
- /// Histogram for PULL Distribution
- CbmQaHist<TH1D> fh1DpullU{"h1DpullU", "Trd1D: Pull U;(U_{reco} - U_{MC}) / #sigmaU_{reco}", 100, -5, 5};
- CbmQaHist<TH1D> fh1DpullV{"h1DpullV", "Trd1D: Pull V;(V_{reco} - V_{MC}) / #sigmaV_{reco}", 100, -5, 5};
- CbmQaHist<TH1D> fh1DpullT{"h1DpullT", "Trd1D: Pull T;(T_{reco} - T_{MC}) / #sigmaT_{reco}", 100, -5, 5};
-
- CbmQaHist<TH1D> fh2DpullX{"h2DpullX", "Trd2D: Pull X;(X_{reco} - X_{MC}) / #sigmaX_{reco}", 100, -5, 5};
- CbmQaHist<TH1D> fh2DpullY{"h2DpullY", "Trd2D: Pull Y;(Y_{reco} - Y_{MC}) / #sigmaY_{reco}", 100, -5, 5};
- CbmQaHist<TH1D> fh2DpullT{"h2DpullT", "Trd2D: Pull T;(T_{reco} - T_{MC}) / #sigmaT_{reco}", 100, -5, 5};
-
- /// List of the above histogramms
- std::vector<CbmQaHist<TH1D>*> fHistList;
-
- /// hits purity
- std::vector<CbmQaHist<TH1I>> fhPointsPerHit;
-
- /// hits efficiency
- std::vector<CbmQaHist<TH1I>> fhHitsPerPoint;
-
- /// hits efficiency
- std::vector<CbmQaHist<TProfile2D>> fhEfficiencyXY;
- std::vector<CbmQaHist<TProfile>> fhEfficiencyR;
-
- /// Canvaces: collection of histogramms
-
- CbmQaCanvas fCanvResidual{"cResidual", "Residual Distribution", 3 * 500, 2 * 500};
- CbmQaCanvas fCanvPull{"cPull", "Pull Distribution", 3 * 500, 2 * 500};
- CbmQaCanvas fCanvEfficiencyXY{"cEfficiencyXY", "Efficiency XY: % reconstructed McPoint", 2 * 500, 2 * 500};
- CbmQaCanvas fCanvEfficiencyR{"cEfficiencyR", "Efficiency R: % reconstructed McPoint", 2 * 500, 2 * 500};
- CbmQaCanvas fCanvHitsPerPoint{"cHitsPerMcPoint", "Efficiency: Hits Per McPoint", 2 * 500, 2 * 500};
- CbmQaCanvas fCanvPointsPerHit{"cMcPointsPerHit", "Purity: McPoints per Hit", 2 * 500, 2 * 500};
-
- private:
- /// Suppressed copy constructor
- CbmTrackerInputQaTrd(const CbmTrackerInputQaTrd&) = delete;
-
- /// Suppressed assignment operator
- CbmTrackerInputQaTrd& operator=(const CbmTrackerInputQaTrd&) = delete;
-
- ClassDefOverride(CbmTrackerInputQaTrd, 0);
-};
diff --git a/reco/L1/utils/CbmCaIdealHitProducerDet.h b/reco/L1/utils/CbmCaIdealHitProducerDet.h
index eaff2a8243977e482a2eeb3a10f74ed46269b434..35c558cb6dee65d77c36489246cffa4064db4498 100644
--- a/reco/L1/utils/CbmCaIdealHitProducerDet.h
+++ b/reco/L1/utils/CbmCaIdealHitProducerDet.h
@@ -527,6 +527,9 @@ namespace cbm::ca
// Get setting
int iSt = fpDetInterface->GetTrackingStationIndex(pHit);
+ if (iSt < 0) {
+ continue;
+ }
const auto& setting = fvStationPars[iSt];
// Skip hits from stations, where new hits are to be created from points
@@ -621,6 +624,10 @@ namespace cbm::ca
}
auto* pPoint = static_cast<Point_t*>(fpBrPoints->Get(fileId, eventId, iP));
int iSt = fpDetInterface->GetTrackingStationIndex(pPoint->GetDetectorID());
+ if (iSt < 0) {
+ continue;
+ }
+
const auto& setting = fvStationPars[iSt];
if (setting.fMode != 2) {
continue;