diff --git a/core/base/CbmTrackingDetectorInterfaceBase.cxx b/core/base/CbmTrackingDetectorInterfaceBase.cxx index 1cc11e8fc38f60fea3f8ef7df46cfca1f3162c46..8dc142c8c0fb27bc5f166db7c0d02a34e5c56ffd 100644 --- a/core/base/CbmTrackingDetectorInterfaceBase.cxx +++ b/core/base/CbmTrackingDetectorInterfaceBase.cxx @@ -45,38 +45,43 @@ bool CbmTrackingDetectorInterfaceBase::Check() const } } - // Size along X-axis - auto xMax = this->GetXmax(iSt); - if (xMax < std::numeric_limits<double>::epsilon() || std::isnan(xMax)) { - msg << prefix << " zero, negative or NaN X-size (" << xMax << " cm)\n"; - res = false && res; - } + { // Size along X-axis + auto xMax = this->GetXmax(iSt); + auto xMin = this->GetXmin(iSt); + if (!(xMax > xMin) || std::isnan(xMax) || std::isnan(xMin)) { + msg << prefix << " zero, negative or NaN X-size (xMin = " << xMin << ", xMax = " << xMax << " cm)\n"; + res = false && res; + } + } - // Size along Y-axis - auto yMax = this->GetYmax(iSt); - if (yMax < std::numeric_limits<double>::epsilon() || std::isnan(yMax)) { - msg << prefix << " zero, negative or NaN Y-size (" << yMax << " cm)\n"; - res = false && res; + { // Size along Y-axis + auto yMax = this->GetYmax(iSt); + auto yMin = this->GetYmin(iSt); + if (!(yMax > yMin) || std::isnan(yMax) || std::isnan(yMin)) { + msg << prefix << " zero, negative or NaN Y-size (yMin = " << yMin << ", yMax = " << yMax << " cm)\n"; + res = false && res; + } } } - // Position along beam axis - std::vector<double> zPositions(this->GetNtrackingStations()); - for (int iSt = 0; iSt < this->GetNtrackingStations(); ++iSt) { - zPositions[iSt] = this->GetZ(iSt); - } - std::set<double> zPositionSet(zPositions.begin(), zPositions.end()); - if (zPositions.size() != zPositionSet.size()) { - msg << "\t- Some of stations have the same z position component:\n"; - for (size_t iSt = 0; iSt < zPositions.size(); ++iSt) { - msg << "\t\tstation " << iSt << ", z = " << zPositions[iSt] << " cm\n"; + { // Position along beam axis + std::vector<double> zPositions(this->GetNtrackingStations()); + for (int iSt = 0; iSt < this->GetNtrackingStations(); ++iSt) { + zPositions[iSt] = this->GetZ(iSt); + } + std::set<double> zPositionSet(zPositions.begin(), zPositions.end()); + if (zPositions.size() != zPositionSet.size()) { + msg << "\t- Some of stations have the same z position component:\n"; + for (size_t iSt = 0; iSt < zPositions.size(); ++iSt) { + msg << "\t\tstation " << iSt << ", z = " << zPositions[iSt] << " cm\n"; + } + res = false && res; } - res = false && res; } } if (!res) { - LOG(error) << msg.str() << "\033[4mErrors above mean that CA tracking cannot be used for the current version of " + LOG(error) << msg.str() << "\033[4mErrors above mean that the CA tracking cannot be used with the current version of " << this->GetDetectorName() << " setup. Please, check if the " << this->GetDetectorName() << " setup parameters and the corresponding tracking detector interface are initialized properly\033[0m"; } diff --git a/core/detectors/tof/CbmTofTrackingInterface.cxx b/core/detectors/tof/CbmTofTrackingInterface.cxx index 871f8d7c9e08997e4fe1eeb10f62e4594e702ad6..11707595da51c000fd326e394929e8f0186132fc 100644 --- a/core/detectors/tof/CbmTofTrackingInterface.cxx +++ b/core/detectors/tof/CbmTofTrackingInterface.cxx @@ -16,6 +16,7 @@ #include "FairDetector.h" #include "FairRunAna.h" #include <Logger.h> +#include <limits> ClassImp(CbmTofTrackingInterface) @@ -40,29 +41,80 @@ CbmTofTrackingInterface::~CbmTofTrackingInterface() InitStatus CbmTofTrackingInterface::Init() { // create digitization parameters from geometry file - CbmTofCreateDigiPar* tofDigiPar = new CbmTofCreateDigiPar("TOF Digi Producer", "TOF task"); + auto tofDigiPar = CbmTofCreateDigiPar("TOF Digi Producer", "TOF task"); LOG(info) << "Create DigiPar"; - tofDigiPar->Init(); - // ** ToF tracking station position z-components initialization ** + tofDigiPar.Init(); + + // ** ToF tracking station geometrical information initialization ** // Init ToF stations position z-components. For each ToF tracking station the position z-component is calculated // as an average of the components for each ToF module inside the tracking station. std::vector<int> nTofStationModules(this->GetNtrackingStations(), 0); // Number of ToF modules for a given tracking station + auto nStations = this->GetNtrackingStations(); + fTofStationXMin.clear(); + fTofStationXMax.clear(); + fTofStationYMin.clear(); + fTofStationYMax.clear(); fTofStationZ.clear(); - fTofStationZ.resize(this->GetNtrackingStations(), 0); - - for (int iSmType = 0; iSmType < fDigiBdfPar->GetNbSmTypes(); iSmType++) { - for (int iSm = 0; iSm < fDigiBdfPar->GetNbSm(iSmType); iSm++) { - for (int iRpc = 0; iRpc < fDigiBdfPar->GetNbRpc(iSmType); iRpc++) { - int iAddr = CbmTofAddress::GetUniqueAddress(iSm, iRpc, 0, 0, iSmType); - int station = fDigiBdfPar->GetTrackingStation(iSmType, iSm, iRpc); - auto* channelInfo = dynamic_cast<CbmTofCell*>(fDigiPar->GetCell(iAddr)); - if (nullptr == channelInfo) { break; } - if (station < 0) { continue; } - fTofStationZ[station] += channelInfo->GetZ(); - nTofStationModules[station] += 1; + fTofStationZMin.clear(); + fTofStationZMax.clear(); + fTofStationXMin.resize(nStations, std::numeric_limits<double>::max()); + fTofStationXMax.resize(nStations, std::numeric_limits<double>::lowest()); + fTofStationYMin.resize(nStations, std::numeric_limits<double>::max()); + fTofStationYMax.resize(nStations, std::numeric_limits<double>::lowest()); + fTofStationZ.resize(nStations, 0); + fTofStationZMin.resize(nStations, std::numeric_limits<double>::max()); + fTofStationZMax.resize(nStations, std::numeric_limits<double>::lowest()); + + for (int iSmType = 0; iSmType < fDigiBdfPar->GetNbSmTypes(); ++iSmType) { + for (int iSm = 0; iSm < fDigiBdfPar->GetNbSm(iSmType); ++iSm) { + for (int iRpc = 0; iRpc < fDigiBdfPar->GetNbRpc(iSmType); ++iRpc) { + for (int iCh = 0; iCh < fDigiBdfPar->GetNbChan(iSmType, iRpc); ++iCh) { + for (int iSide = 0; iSide < 1; ++iSide) { + auto address = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iCh, iSide, iSmType); + + int iStation = fDigiBdfPar->GetTrackingStation(iSmType, iSm, iRpc); // Local index of tracking station + + auto* pChannelInfo = dynamic_cast<CbmTofCell*>(fDigiPar->GetCell(address)); + if (nullptr == pChannelInfo) { + LOG(warn) << fName << ": CbmTofCell object is not defined for iSmType = " << iSmType << ", iSm = " << iSm + << ", iRpc = " << iRpc; + continue; + } + + // Tracking station sizes + auto chPosX = pChannelInfo->GetX(); + auto chPosY = pChannelInfo->GetY(); + auto chPosZ = pChannelInfo->GetZ(); + auto chSizeX = pChannelInfo->GetSizex(); + auto chSizeY = pChannelInfo->GetSizey(); + auto chLoBoarderX = chPosX - chSizeX; + auto chUpBoarderX = chPosX + chSizeX; + auto chLoBoarderY = chPosY - chSizeY; + auto chUpBoarderY = chPosY + chSizeY; + + LOG(info) << "TOF Cell: iSmType = " << iSmType << ", iSm = " << iSm << ", iRpc = " << iRpc + << ", iCh = " << iCh << ", iSide = " << iSide << ", addr = " << address << ", sta = " + << iStation << ", pos = " << chPosX << ',' << chPosY << ',' << chPosZ << ", xRange = " + << chLoBoarderX << "," << chUpBoarderX << ", yRange = " << chLoBoarderY << "," << chUpBoarderY; + + // Cuts on T0 and undefined station ID + if (5 == iSmType) { continue; } // Skip T0 + if (iStation < 0) { continue; } + + fTofStationZ[iStation] += chPosZ; + if (chPosZ > fTofStationZMax[iStation]) { fTofStationZMax[iStation] = chPosZ; } + if (chPosZ < fTofStationZMin[iStation]) { fTofStationZMin[iStation] = chPosZ; } + if (chUpBoarderX > fTofStationXMax[iStation]) { fTofStationXMax[iStation] = chUpBoarderX; } + if (chUpBoarderY > fTofStationYMax[iStation]) { fTofStationYMax[iStation] = chUpBoarderY; } + if (chLoBoarderX < fTofStationXMin[iStation]) { fTofStationXMin[iStation] = chLoBoarderX; } + if (chLoBoarderY < fTofStationYMin[iStation]) { fTofStationYMin[iStation] = chLoBoarderY; } + + nTofStationModules[iStation] += 1; + } + } } } } diff --git a/core/detectors/tof/CbmTofTrackingInterface.h b/core/detectors/tof/CbmTofTrackingInterface.h index 81833eab0ff3f9a422e3169f34eddbb872fb99ee..be2c2a37aa260ff26ff5e2cc169bc529d8b4b6d1 100644 --- a/core/detectors/tof/CbmTofTrackingInterface.h +++ b/core/detectors/tof/CbmTofTrackingInterface.h @@ -25,7 +25,6 @@ #include "TMath.h" -#include <iostream> #include <vector> class CbmTofDigiPar; @@ -80,29 +79,29 @@ public: return iSt; } - // TODO: SZh 30.08.2023: Provide automatic definintion of bounds + // TODO: SZh. 10.09.2023: Provide automatic selection /// @brief Gets upper bound of a station along the X-axis /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) /// @return Size of station along the X-axis [cm] - double GetXmax(int /*stationId*/) const { return 100.; } + double GetXmax(int stationId) const { return 100; } - // TODO: SZh 30.08.2023: Provide automatic definintion of bounds + // TODO: SZh. 10.09.2023: Provide automatic definition /// @brief Gets lower bound of a station along the X-axis /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) /// @return Size of station along the X-axis [cm] - double GetXmin(int /*stationId*/) const { return -100.; } + double GetXmin(int stationId) const { return -100; } - // TODO: SZh 30.08.2023: Provide automatic definintion of bounds + // TODO: SZh. 10.09.2023: Provide automatic definition /// @brief Gets max size of a station along the Y-axis /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) /// @return Size of station along the Y-axis [cm] - double GetYmax(int /*stationId*/) const { return 100.; } + double GetYmax(int stationId) const { return 100; } - // TODO: SZh 30.08.2023: Provide automatic definintion of bounds + // TODO: SZh. 10.09.2023: Provide automatic definition /// @brief Gets max size of a station along the Y-axis /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) /// @return Size of station along the Y-axis [cm] - double GetYmin(int /*stationId*/) const { return -100.; } + double GetYmin(int stationId) const { return -100; } /// @brief Gets reference z of the station /// @param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) @@ -112,12 +111,12 @@ public: /// Gets max z of the station /// \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 GetZ(stationId) + 5.; } + double GetZmax(int stationId) const { return fTofStationZMax[stationId] + 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 GetZ(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]) @@ -142,7 +141,13 @@ private: CbmTofDigiPar* fDigiPar {nullptr}; CbmTofDigiBdfPar* fDigiBdfPar {nullptr}; - std::vector<float> fTofStationZ {}; ///< ToF stations position z-component [cm] + std::vector<double> fTofStationXMin {}; ///< Lower bounds of TOF stations along x-axis [cm] + std::vector<double> fTofStationXMax {}; ///< Upper bounds of TOF stations along x-axis [cm] + std::vector<double> fTofStationYMin {}; ///< Lower bounds of TOF stations along y-axis [cm] + std::vector<double> fTofStationYMax {}; ///< Upper bounds of TOF stations along y-axis [cm] + std::vector<double> fTofStationZ {}; ///< Centers of TOF stations along z-axis [cm] + std::vector<double> fTofStationZMin {}; ///< Lower bounds of TOF stations along z-axis [cm] + std::vector<double> fTofStationZMax {}; ///< Upper bounds of TOF stations along z-axis [cm] ClassDef(CbmTofTrackingInterface, 0); }; diff --git a/core/qa/CbmQaIO.h b/core/qa/CbmQaIO.h index b99f475cd141b486c9d85d95809456c1da2ae6d5..ef89c77ea68f726fef5df4b34d3bfcce8d3684f5 100644 --- a/core/qa/CbmQaIO.h +++ b/core/qa/CbmQaIO.h @@ -252,6 +252,8 @@ T* CbmQaIO::MakeQaObject(TString sName, TString sTitle, Args... args) } } + //LOG(info) << "CbmQaIO: init ROOT object \"" << sName << '\"'; // Debug + const char* title = sTitle.Data(); if (bUseConfig) { if constexpr (std::is_base_of_v<TH1, T>) { diff --git a/macro/mcbm/mcbm_qa.C b/macro/mcbm/mcbm_qa.C index 52709d434d06653a34f011fe97ccbfde35765698..1a67d7523fdeaa69e0c3181fde4937955badbad1 100644 --- a/macro/mcbm/mcbm_qa.C +++ b/macro/mcbm/mcbm_qa.C @@ -212,28 +212,27 @@ void mcbm_qa(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_2020_03_test", run->AddTask(muchHitFinderQa); // CA Input QA - //auto* pCaInputMuch = new CbmCaInputQaMuch(verbose, bUseMC); - //pCaInputMuch->SetEfficiencyThrsh(0.5, 0, 100); - //run->AddTask(pCaInputMuch); + auto* pCaInputMuch = new CbmCaInputQaMuch(verbose, bUseMC); + pCaInputMuch->SetEfficiencyThrsh(0.5, 0, 100); + run->AddTask(pCaInputMuch); } // ------------------------------------------------------------------------ // ----- TRD QA ----------------------------------------------------------- if (bUseTrd) { // CA Input QA - //auto* pCaInputTrd = new CbmCaInputQaTrd(verbose, bUseMC); - //pCaInputTrd->SetEfficiencyThrsh(0.5, 0, 100); - ////pCaInputTrd->SetConfigName(qaConfig); - //run->AddTask(pCaInputTrd); + auto* pCaInputTrd = new CbmCaInputQaTrd(verbose, bUseMC); + pCaInputTrd->SetEfficiencyThrsh(0.5, 0, 100); + run->AddTask(pCaInputTrd); } // ------------------------------------------------------------------------ // ----- TOF QA ----------------------------------------------------------- if (bUseTof) { // CA Input QA - //auto* pCaInputTof = new CbmCaInputQaTof(verbose, bUseMC); - //pCaInputTof->SetEfficiencyThrsh(0.5, 0, 100); - //run->AddTask(pCaInputTof); + auto* pCaInputTof = new CbmCaInputQaTof(verbose, bUseMC); + pCaInputTof->SetEfficiencyThrsh(0.5, 0, 100); + run->AddTask(pCaInputTof); } // ------------------------------------------------------------------------ diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt index 0b74f0cc066cefc9c6186ef2c8a820efc2809676..4f849847864423aae9a64785775450b71a7ca249 100644 --- a/reco/L1/CMakeLists.txt +++ b/reco/L1/CMakeLists.txt @@ -81,10 +81,11 @@ set(SRCS qa/CbmTrackerInputQaTrd.cxx qa/CbmTrackerInputQaTof.cxx qa/CbmCaInputQaBase.cxx + qa/CbmCaInputQaMvd.cxx qa/CbmCaInputQaSts.cxx - #qa/CbmCaInputQaMuch.cxx - #qa/CbmCaInputQaTrd.cxx - #qa/CbmCaInputQaTof.cxx + qa/CbmCaInputQaMuch.cxx + qa/CbmCaInputQaTrd.cxx + qa/CbmCaInputQaTof.cxx qa/CbmCaInputQaSetup.cxx qa/CbmCaOutputQa.cxx qa/CbmCaTrackTypeQa.cxx diff --git a/reco/L1/L1LinkDef.h b/reco/L1/L1LinkDef.h index c1dc866aae060386d505547df0245713a05101c7..3e232991793e3649e68bd265abca7e4d882cf563 100644 --- a/reco/L1/L1LinkDef.h +++ b/reco/L1/L1LinkDef.h @@ -25,10 +25,11 @@ #pragma link C++ class CbmTrackerInputQaTrd + ; #pragma link C++ class CbmTrackerInputQaTof + ; #pragma link C++ class CbmCaInputQaBase + ; -//#pragma link C++ class CbmCaInputQaMuch + ; +#pragma link C++ class CbmCaInputQaMvd + ; +#pragma link C++ class CbmCaInputQaMuch + ; #pragma link C++ class CbmCaInputQaSts + ; -//#pragma link C++ class CbmCaInputQaTrd + ; -//#pragma link C++ class CbmCaInputQaTof + ; +#pragma link C++ class CbmCaInputQaTrd + ; +#pragma link C++ class CbmCaInputQaTof + ; #pragma link C++ enum cbm::ca::ETrackType; #pragma link C++ enum class L1DetectorID; #pragma link C++ class cbm::ca::OutputQa + ; diff --git a/reco/L1/qa/CbmCaHitQaData.cxx b/reco/L1/qa/CbmCaHitQaData.cxx new file mode 100644 index 0000000000000000000000000000000000000000..0c295c22ff0511b1a4cf5e66cea1f1b3055c8efd --- /dev/null +++ b/reco/L1/qa/CbmCaHitQaData.cxx @@ -0,0 +1,20 @@ +/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// @file CbmCaHitQaData.cxx +/// @date 01.09.2023 +/// @brief A helper class to store hit and MC-point parameter and calculate related quantities (implementation) +/// @author S.Zharko <s.zharko@gsi.de> + +#include "CbmCaHitQaData.h" + +using cbm::ca::HitQaData; + +// --------------------------------------------------------------------------------------------------------------------- +// +void HitQaData::Reset() +{ + this->operator=(HitQaData()); +} + diff --git a/reco/L1/qa/CbmCaHitQaData.h b/reco/L1/qa/CbmCaHitQaData.h new file mode 100644 index 0000000000000000000000000000000000000000..b453d0c69b4f35b776dd3559b5cbb19a4e8c3aec --- /dev/null +++ b/reco/L1/qa/CbmCaHitQaData.h @@ -0,0 +1,324 @@ +/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// @file CbmCaHitQaData.h +/// @date 01.09.2023 +/// @brief A helper class to store hit and MC-point parameter and calculate related quantities (header) +/// @author S.Zharko <s.zharko@gsi.de> + +#ifndef CbmCaHitQaData_h +#define CbmCaHitQaData_h 1 + +#include <cmath> +#include <limits> +#include <tuple> + +namespace cbm::ca +{ + /// @class HitQaData + /// @brief Contains necessary data to calculate hit residuals and pulls + class HitQaData { + public: + /// @brief Default constructor + HitQaData() = default; + + /// @brief Destructor + ~HitQaData() = default; + + /// @brief Copy constructor + HitQaData(const HitQaData& ) = default; + + /// @brief Move constructor + HitQaData(HitQaData&& ) = default; + + /// @brief Copy assignment operator + HitQaData& operator=(const HitQaData& ) = default; + + /// @brief Move assignment operator + HitQaData& operator=(HitQaData&& ) = default; + + /// @brief Gets hit u-coordinate error + /// @return Hit u-coordinate error [cm] + double GetHitDu() const + { + auto cU = cos(fPhiU); + auto sU = sin(fPhiU); + return sqrt(cU * cU * GetHitDx() * GetHitDx() + 2. * sU * cU * GetHitDxy() + sU * sU * GetHitDy() * GetHitDy()); + } + + /// @brief Gets hit v-coordinate error + /// @return Hit v-coordinate error [cm] + double GetHitDv() const + { + auto cV = cos(fPhiV); + auto sV = sin(fPhiV); + return sqrt(cV * cV * GetHitDx() * GetHitDx() + 2. * sV * cV * GetHitDxy() + sV * sV * GetHitDy() * GetHitDy()); + } + + /// @brief Gets hit u- and v-coordinate covariance + /// @return Hit u- and v-coordinate covariance [cm2] + double GetHitDuv() const + { + auto cU = cos(fPhiU); + auto sU = sin(fPhiU); + auto cV = cos(fPhiV); + auto sV = sin(fPhiV); + return GetHitDx() * GetHitDx() * cU * cV + GetHitDxy() * (sU * cV + cU * sV) + GetHitDy() * GetHitDy() * sU * sV; + } + + /// @brief Gets hit x-coordinate error + /// @return Hit x-coordinate error [cm] + double GetHitDx() const { return fHitDx; } + + /// @brief Gets hit x- and y-coordinate covariance + /// @return Hit x- and y-coordinate covariance [cm2] + double GetHitDxy() const { return fHitDxy; } + + /// @brief Gets hit y-coordinate error + /// @return Hit y-coordinate error [cm] + double GetHitDy() const { return fHitDy; } + + /// @brief Gets hit index + /// @return Hit index + int GetHitIndex() const { return fHitID; } + + /// @brief Gets Pearson correlation coefficient for u- and v-coordinates + /// @return Pearson correlation coefficient for u- and v-coordinates + double GetHitRuv() const { return GetHitDuv() / GetHitDu() / GetHitDv(); } + + /// @brief Gets hit time + /// @return hit time [ns] + double GetHitTime() const { return fHitTime; } + + /// @brief Gets hit time error + /// @return Hit time error [ns] + double GetHitTimeError() const { return fHitTimeError; } + + /// @brief Gets hit u-coordinate + /// @return hit u-coordinate [cm] + double GetHitU() const { return GetHitX() * cos(fPhiU) + GetHitY() * sin(fPhiU); } + + /// @brief Gets hit v-coordinate + /// @return hit v-coordinate [cm] + double GetHitV() const { return GetHitX() * cos(fPhiV) + GetHitY() * sin(fPhiV); } + + /// @brief Gets hit x-coordinate + /// @return hit x-coordinate [cm] + double GetHitX() const { return fHitX; } + + /// @brief Gets hit y-coordinate + /// @return hit y-coordinate [cm] + double GetHitY() const { return fHitY; } + + /// @brief Gets hit z-coordinate + /// @return hit z-coordinate [cm] + double GetHitZ() const { return fHitZ; } + + /// @brief Gets Flag: if track has hits + /// @return Flag: if track has hits + bool GetIfTrackHasHits() const { return fbTrackHasHits; } + + /// @brief Gets Flag: if track selected + /// @return Flag: if track selected + bool GetIfTrackSelected() const { return fbTrackSelected; } + + /// @brief Gets residual of time + /// @return Residual of time [ns] + double GetResidualTime() const { return GetHitTime() - GetPointTime(); } + + /// @brief Gets residual of u-coordinate + /// @return Residual of u-coordinate [cm] + double GetResidualU() const { return GetHitU() - GetPointU(); } + + /// @brief Gets residual of v-coordinate + /// @return Residual of v-coordinate [cm] + double GetResidualV() const { return GetHitV() - GetPointV(); } + + /// @brief Gets residual of x-coordinate + /// @return Residual of x-coordinate [cm] + double GetResidualX() const { return GetHitX() - GetPointX(); } + + /// @brief Gets residual of y-coordinate + /// @return Residual of y-coordinate [cm] + double GetResidualY() const { return GetHitY() - GetPointY(); } + + /// @brief Gets front strips stereo angle + /// @return Front strips stereo angle [rad] + double GetPhiU() const { return fPhiU; } + + /// @brief Gets back strips stereo angle + /// @return Back strips stereo angle [rad] + double GetPhiV() const { return fPhiV; } + + /// @brief Gets point index + /// @return A tuple (pointID, eventID, fileID) + std::tuple<int, int, int> GetPointID() const + { + return std::make_tuple(fPointID, fMCEventID, fMCFileID); + } + + /// @brief Gets point time + /// @return Point time [ns] + double GetPointTime() const { return fPointTime; } + + /// @brief Gets point u-coordinate + /// @return Point u-coordinate [cm] + double GetPointU() const { return GetPointX() * cos(fPhiU) + GetPointY() * sin(fPhiU); } + + /// @brief Gets point v-coordinate + /// @return Point v-coordinate [cm] + double GetPointV() const { return GetPointX() * cos(fPhiV) + GetPointY() * sin(fPhiV); } + + /// @brief Gets point x-coordinate + /// @return Point x-coordinate [cm] + double GetPointX() const { return fPointX; } + + /// @brief Gets point y-coordinate + /// @return Point y-coordinate [cm] + double GetPointY() const { return fPointY; } + + /// @brief Gets point z-coordinate + /// @param pointZ Point z-coordinate [cm] + double GetPointZ() const { return fPointZ; } + + /// @brief Gets pull of time + /// @return Pull of time + double GetPullTime() const { return GetResidualTime() / GetHitTimeError(); } + + /// @brief Gets pull of u-coordinate + /// @return Pull of u-coordinate + double GetPullU() const { return GetResidualU() / GetHitDu(); } + + /// @brief Gets pull of v-coordinate + /// @return Pull of v-coordinate + double GetPullV() const { return GetResidualV() / GetHitDv(); } + + /// @brief Gets pull of x-coordinate + /// @return Pull of x-coordinate + double GetPullX() const { return GetResidualX() / GetHitDx(); } + + /// @brief Gets pull of y-coordinate + /// @return Pull of y-coordinate + double GetPullY() const { return GetResidualY() / GetHitDy(); } + + /// @brief Gets station local index + /// @return Station local index + int GetStationID() const { return fStationID; } + + /// @brief Resets data fields + void Reset() { this->operator=(std::move(HitQaData())); } + + /// @brief Sets hit x-coordinate error + /// @param hitDx Hit x-coordinate error [cm] + void SetHitDx(double hitDx) { fHitDx = hitDx; } + + /// @brief Sets hit x- and y-coordinate covariance + /// @param hitDxy Hit x- and y-coordinate covariance [cm2] + void SetHitDxy(double hitDxy) { fHitDxy = hitDxy; } + + /// @brief Sets hit y-coordinate error + /// @param hitDy Hit y-coordinate error [cm] + void SetHitDy(double hitDy) { fHitDy = hitDy; } + + /// @brief Sets hit index + /// @param hitID Hit index + void SetHitIndex(int hitID) { fHitID = hitID; } + + /// @brief Sets hit time + /// @param hitTime Hit time [ns] + void SetHitTime(double hitTime) { fHitTime = hitTime; } + + /// @brief Sets hit time error + /// @param hitTimeError Hit time error [ns] + void SetHitTimeError(double hitTimeError) { fHitTimeError = hitTimeError; } + + /// @brief Sets hit x-coordinate + /// @param hitX Hit x-coordinate [cm] + void SetHitX(double hitX) { fHitX = hitX; } + + /// @brief Sets hit y-coordinate + /// @param hitY Hit y-coordinate [cm] + void SetHitY(double hitY) { fHitY = hitY; } + + /// @brief Sets hit z-coordinate + /// @param hitZ Hit z-coordinate [cm] + void SetHitZ(double hitZ) { fHitZ = hitZ; } + + /// @brief Sets Flag: if track has hits + /// @param ifTrackHasHits Flag: if track has hits + void SetIfTrackHasHits(bool ifTrackHasHits) { fbTrackHasHits = ifTrackHasHits; } + + /// @brief Sets Flag: if track selected + /// @param ifTrackSelected Flag: if track selected + void SetIfTrackSelected(bool ifTrackSelected) { fbTrackSelected = ifTrackSelected; } + + /// @brief Sets front strips stereo angle + /// @param phiU Front strips stereo angle [rad] + void SetPhiU(double phiU) { fPhiU = phiU; } + + /// @brief Sets back strips stereo angle + /// @param phiV Back strips stereo angle [rad] + void SetPhiV(double phiV) { fPhiV = phiV; } + + /// @brief Sets point index + /// @param pointID Index of point + /// @param eventID Index of MC event + /// @param fileID Index of MC file + void SetPointID(int pointID, int eventID, int fileID) + { + fPointID = pointID; + fMCEventID = eventID; + fMCFileID = fileID; + } + + /// @brief Sets point time + /// @param pointTime Point time [ns] + void SetPointTime(double pointTime) { fPointTime = pointTime; } + + /// @brief Sets point x-coordinate + /// @param pointX Point x-coordinate [cm] + void SetPointX(double pointX) { fPointX = pointX; } + + /// @brief Sets point y-coordinate + /// @param pointY Point y-coordinate [cm] + void SetPointY(double pointY) { fPointY = pointY; } + + /// @brief Sets point z-coordinate + /// @param pointZ Point z-coordinate [cm] + void SetPointZ(double pointZ) { fPointZ = pointZ; } + + /// @brief Sets station local index + /// @return Station local index + void SetStationID(int iStLoc) { fStationID = iStLoc; } + + private: + static constexpr double kNAN = std::numeric_limits<double>::signaling_NaN(); + + double fPhiU = kNAN; ///< Stereo angle for front strips [rad] + double fPhiV = kNAN; ///< Stereo anele for back strips [rad] + double fHitX = kNAN; ///< Hit x-coordinate [cm] + double fHitY = kNAN; ///< Hit y-coordinate [cm] + double fHitZ = kNAN; ///< Hit z-coordinate [cm] + double fHitTime = kNAN; ///< Hit time [ns] + double fHitDx = kNAN; ///< Hit x-coordinate error [cm] + double fHitDy = kNAN; ///< Hit y-coordinate error [cm] + double fHitDxy = kNAN; ///< Hit x- and y-coordinate covariance [cm2] + double fHitTimeError = kNAN; ///< Hit time error [ns] + double fPointX = kNAN; ///< Point x-coordinate [cm] + double fPointY = kNAN; ///< Point y-coordinate [cm] + double fPointZ = kNAN; ///< Point z-coordinate [cm] + double fPointTime = kNAN; ///< Point time [ns] + int fStationID = -1; ///< Local index of tracking station + int fHitID = -1; ///< Index of hit + int fPointID = -1; ///< Index of MC point + int fMCEventID = -1; ///< Index of MC event + int fMCFileID = -1; ///< Index of MC file id + bool fbTrackSelected = false; ///< Flag: if track selected + bool fbTrackHasHits = false; ///< Flag: if track has hits + + }; + +} + +#endif // CbmCaHitQaData_h diff --git a/reco/L1/qa/CbmCaInputQaBase.cxx b/reco/L1/qa/CbmCaInputQaBase.cxx new file mode 100644 index 0000000000000000000000000000000000000000..552a16737318f56d30f7a1413f836070285b9bb2 --- /dev/null +++ b/reco/L1/qa/CbmCaInputQaBase.cxx @@ -0,0 +1,1423 @@ +/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// @file CbmCaInputQaBase.cxx +/// @date 30.08.2023 +/// @brief Class providing basic functionality for CA input QA-tasks (implementation) +/// @author S.Zharko <s.zharko@gsi.de> + +#include "CbmCaInputQaBase.h" +#include "CbmTrackingDetectorInterfaceBase.h" +#include "Logger.h" +#include "CbmAddress.h" +#include "CbmMCDataArray.h" +#include "CbmMCEventList.h" +#include "CbmMCTrack.h" +#include "CbmMatch.h" +#include "CbmQaCanvas.h" +#include "CbmQaTable.h" +#include "CbmQaUtil.h" +#include "CbmMvdHit.h" +#include "CbmMvdPoint.h" +#include "CbmStsCluster.h" +#include "CbmStsHit.h" +#include "CbmStsPoint.h" +#include "CbmMuchPixelHit.h" +#include "CbmMuchPoint.h" +#include "CbmTofAddress.h" +#include "CbmTrdHit.h" +#include "CbmTrdPoint.h" +#include "CbmTofHit.h" +#include "CbmTofPoint.h" + +#include "CbmTimeSlice.h" +#include "FairMCPoint.h" + +#include "FairRootManager.h" +#include "Logger.h" + +#include "TBox.h" +#include "TClonesArray.h" +#include "TEllipse.h" +#include "TF1.h" +#include "TFormula.h" +#include "TGraphAsymmErrors.h" +#include "TH1F.h" +#include "TH2F.h" +#include "TMath.h" +#include "TProfile.h" +#include "TProfile2D.h" +#include "TStyle.h" + +#include <algorithm> +#include <fstream> +#include <iomanip> +#include <numeric> +#include <tuple> + +#include "L1Constants.h" + +namespace phys = L1Constants::phys; // from L1Constants.h + + +// --------------------------------------------------------------------------------------------------------------------- +// +template <L1DetectorID DetID> +CbmCaInputQaBase<DetID>::CbmCaInputQaBase(const char* name, int verbose, bool isMCUsed) +: CbmQaTask(name, verbose, isMCUsed) +{} + +// --------------------------------------------------------------------------------------------------------------------- +// +template <L1DetectorID DetID> +bool CbmCaInputQaBase<DetID>::Check() +{ + bool res = true; + + int nSt = fpDetInterface->GetNtrackingStations(); + + + // ************************************************************** + // ** Basic checks, available both for real and simulated data ** + // ************************************************************** + + // Fill efficiency distributions + + for (int iSt = 0; iSt < nSt; ++iSt) { + TProfile2D* effxy = fvpe_reco_eff_vs_xy[iSt]; + for (int i = 1; i < effxy->GetNbinsX() - 1; i++) { + for (int j = 1; j < effxy->GetNbinsY() - 1; j++) { + int bin = effxy->GetBin(i, j); + if (effxy->GetBinEntries(bin) >= 1) { + fvph_reco_eff[iSt]->Fill(effxy->GetBinContent(bin)); + fvph_reco_eff[nSt]->Fill(effxy->GetBinContent(bin)); + } + } + } + } + + // ----- Checks for mismatches in the ordering of the stations + // + std::vector<double> vStationPos(nSt, 0.); + for (int iSt = 0; iSt < nSt; ++iSt) { + vStationPos[iSt] = fpDetInterface->GetZ(iSt); + } + + if (!std::is_sorted(vStationPos.cbegin(), vStationPos.cend(), [](int l, int r) { return l <= r; })) { + if (fVerbose > 0) { + LOG(error) << fName << ": stations are ordered improperly along the beam axis:"; + for (auto z : vStationPos) { + LOG(error) << "\t- " << z; + } + } + res = false; + } + + // ----- Checks for mismatch between station and hit z positions + // The purpose of this block is to be ensured, that hits belong to the correct tracking station. For each tracking + // station a unified position z-coordinate is defined, which generally differs from the corresponding positions of + // reconstructed hits. This is due to non-regular position along the beam axis for each detector sensor. Thus, the + // positions of the hits along z-axis are distributed relatively to the position of the station in some range. + // If hits goes out this range, it can signal about a mismatch between hits and geometry. For limits of the range + // one should select large enough values to cover the difference described above and in the same time small enough + // to avoid overlaps with the neighboring stations. + // For each station, a distribution of z_{hit} - z_{st} is integrated over the defined range and scaled by the + // total number of entries to the distribution. If this value is smaller then unity, then some of the hits belong to + // another station. + // + for (int iSt = 0; iSt < nSt; ++iSt) { + int nHits = fvph_hit_station_delta_z[iSt]->GetEntries(); + if (!nHits) { + LOG_IF(error, fVerbose > 0) << fName << ": station " << iSt << " does not have hits"; + res = false; + continue; + } + int iBinMin = fvph_hit_station_delta_z[iSt]->FindBin(-fMaxDiffZStHit); + int iBinMax = fvph_hit_station_delta_z[iSt]->FindBin(+fMaxDiffZStHit); + + if (fvph_hit_station_delta_z[iSt]->Integral(iBinMin, iBinMax) < nHits) { + LOG_IF(error, fVerbose > 0) << fName << ": station " << iSt << " has mismatches in hit z-positions"; + res = false; + } + } + + + // ******************************************************* + // ** Additional checks, if MC information is available ** + // ******************************************************* + + if (IsMCUsed()) { + // ----- Check efficiencies + // Error is raised, if any station has integrated efficiency lower then a selected threshold. + // Integrated efficiency is estimated with fitting the efficiency curve for hit distances (r) with a uniform + // distribution in the range from kEffRangeMin to kEffRangeMax, where the efficiency is assigned to be constant + // + // Fit efficiency curves + LOG(info) << "-- Hit efficiency integrated over hit distance from station center"; + + auto* pEffTable = MakeQaObject<CbmQaTable>("vs Station/eff_table", "Efficiency table", nSt, 2); + pEffTable->SetNamesOfCols({"Station ID", "Efficiency"}); + pEffTable->SetColWidth(20); + + for (int iSt = 0; iSt < nSt; ++iSt) { + auto eff = fvph_reco_eff[iSt]->GetMean(); + pEffTable->SetCell(iSt, 0, iSt); + pEffTable->SetCell(iSt, 1, eff); + res = CheckRange("Hit finder efficiency in station " + std::to_string(iSt), eff, fEffThrsh, 1.000); + } + + LOG(info) << '\n' << pEffTable->ToString(3); + + // ----- Checks for residuals + // Check hits for potential biases from central values + + auto* pResidualsTable = + MakeQaObject<CbmQaTable>("vs Station/residuals_mean", "Residual mean values in different stations", nSt, 4); + pResidualsTable->SetNamesOfCols({"Station ID", "Residual(x) [cm]", "Residual(y) [cm]", "Residual(t) [ns]"}); + pResidualsTable->SetColWidth(20); + + // Fit residuals + for (int iSt = 0; iSt <= nSt; ++iSt) { + + cbm::qa::util::SetLargeStats(fvph_res_x[iSt]); + cbm::qa::util::SetLargeStats(fvph_res_y[iSt]); + cbm::qa::util::SetLargeStats(fvph_res_t[iSt]); + + pResidualsTable->SetCell(iSt, 0, iSt); + pResidualsTable->SetCell(iSt, 1, fvph_res_x[iSt]->GetStdDev()); + pResidualsTable->SetCell(iSt, 2, fvph_res_y[iSt]->GetStdDev()); + pResidualsTable->SetCell(iSt, 3, fvph_res_t[iSt]->GetStdDev()); + } + + LOG(info) << '\n' << pResidualsTable->ToString(8); + + // ----- Checks for pulls + // + // For the hit pull is defined as a ration of the difference between hit and MC-point position or time component + // values and an error of the hit position (or time) component, calculated in the same z-positions. In ideal case, + // when the resolutions are well determined for detecting stations, the pull distributions should have a RMS equal + // to unity. Here we allow a RMS of the pull distributions to be varied in a predefined range. If the RMS runs out + // this range, QA task fails. + + auto* pPullsTable = + MakeQaObject<CbmQaTable>("vs Station/pulls_rms", "Pulls std. dev. values in different stations", nSt, 4); + pPullsTable->SetNamesOfCols({"Station ID", "Pull(x) sigm", "Pull(y) sigm", "Pull(t) sigm"}); + pPullsTable->SetColWidth(20); + + for (int iSt = 0; iSt < nSt + 1; ++iSt) { + + // Fit pull distributions for nicer representation. Fit results are not used in further checks. + + cbm::qa::util::SetLargeStats(fvph_pull_x[iSt]); + cbm::qa::util::SetLargeStats(fvph_pull_y[iSt]); + cbm::qa::util::SetLargeStats(fvph_pull_t[iSt]); + + // Check the pull quality + res = CheckRangePull(fvph_pull_x[iSt]) && res; + res = CheckRangePull(fvph_pull_y[iSt]) && res; + res = CheckRangePull(fvph_pull_t[iSt]) && res; + + pPullsTable->SetCell(iSt, 0, iSt); + pPullsTable->SetCell(iSt, 1, fvph_pull_x[iSt]->GetStdDev()); + pPullsTable->SetCell(iSt, 2, fvph_pull_y[iSt]->GetStdDev()); + pPullsTable->SetCell(iSt, 3, fvph_pull_t[iSt]->GetStdDev()); + } + + LOG(info) << '\n' << pPullsTable->ToString(3); + } // McUsed + + return res; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +template <L1DetectorID DetID> +void CbmCaInputQaBase<DetID>::DeInit() +{ + // Vectors with pointers to histograms + fvph_hit_xy.clear(); + fvph_hit_zx.clear(); + fvph_hit_zy.clear(); + + fvph_hit_station_delta_z.clear(); + + fvph_hit_dx.clear(); + fvph_hit_dy.clear(); + fvph_hit_du.clear(); + fvph_hit_dv.clear(); + fvph_hit_kuv.clear(); + fvph_hit_dt.clear(); + + fvph_n_points_per_hit.clear(); + + fvph_point_xy.clear(); + fvph_point_zx.clear(); + fvph_point_zy.clear(); + + fvph_point_hit_delta_z.clear(); + + fvph_res_x.clear(); + fvph_res_y.clear(); + fvph_res_u.clear(); + fvph_res_v.clear(); + fvph_res_t.clear(); + + fvph_pull_x.clear(); + fvph_pull_y.clear(); + fvph_pull_u.clear(); + fvph_pull_v.clear(); + fvph_pull_t.clear(); + + fvph_res_x_vs_x.clear(); + fvph_res_y_vs_y.clear(); + fvph_res_u_vs_u.clear(); + fvph_res_v_vs_v.clear(); + fvph_res_t_vs_t.clear(); + + fvph_pull_x_vs_x.clear(); + fvph_pull_y_vs_y.clear(); + fvph_pull_u_vs_u.clear(); + fvph_pull_v_vs_v.clear(); + fvph_pull_t_vs_t.clear(); + + fvpe_reco_eff_vs_xy.clear(); + fvph_reco_eff.clear(); +} + +// --------------------------------------------------------------------------------------------------------------------- +// +template <L1DetectorID DetID> +void CbmCaInputQaBase<DetID>::FillHistograms() +{ + int nSt = fpDetInterface->GetNtrackingStations(); + int nHits = fpHits->GetEntriesFast(); + int nMCevents = (IsMCUsed()) ? fpMCEventList->GetNofEvents() : -1; + + // TODO: SZh 06.09.2023: Probably, this approach can fail, if there are several input files are used. Thus I propose + // to use unordered_map with a CbmLink key type. + std::vector<std::vector<std::vector<int>>> + vNofHitsPerMcTrack; // Number of hits per MC track per station in different MC events + + if (IsMCUsed()) { + vNofHitsPerMcTrack.resize(nMCevents); + for (int iE = 0; iE < nMCevents; ++iE) { + int iFile = fpMCEventList->GetFileIdByIndex(iE); + int iEvent = fpMCEventList->GetEventIdByIndex(iE); + int nMcTracks = fpMCTracks->Size(iFile, iEvent); + vNofHitsPerMcTrack[iE].resize(nSt); + for (int iSt = 0; iSt < nSt; iSt++) { + vNofHitsPerMcTrack[iE][iSt].resize(nMcTracks, 0); + } + } + } + + for (int iH = 0; iH < nHits; ++iH) { + fHitQaData.Reset(); + 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 constexpr (L1DetectorID::kTof == DetID) { + auto address = pHit->GetAddress(); + if (0x00202806 == address || 0x00002806 == address) { continue; } // TEST + if (5 == CbmTofAddress::GetSmType(address)) { continue; } // Skip T0 hits + } + + // ************************* + // ** Reconstructed hit QA + + // 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; + + auto [stPhiU, stPhiV] = fpDetInterface->GetStereoAnglesSensor(pHit->GetAddress()); + + fHitQaData.SetPhiU(stPhiU); + fHitQaData.SetPhiV(stPhiV); + fHitQaData.SetHitX(pHit->GetX()); + fHitQaData.SetHitY(pHit->GetY()); + fHitQaData.SetHitZ(pHit->GetZ()); + fHitQaData.SetHitTime(pHit->GetTime()); + + fHitQaData.SetHitDx(pHit->GetDx()); + fHitQaData.SetHitDy(pHit->GetDy()); + fHitQaData.SetHitDxy(pHit->GetDxy()); + fHitQaData.SetHitTimeError(pHit->GetTimeError()); + fHitQaData.SetStationID(iSt); + + // Per station distributions + fvph_hit_xy[iSt]->Fill(fHitQaData.GetHitX(), fHitQaData.GetHitY()); + fvph_hit_zx[iSt]->Fill(fHitQaData.GetHitZ(), fHitQaData.GetHitX()); + fvph_hit_zy[iSt]->Fill(fHitQaData.GetHitZ(), fHitQaData.GetHitY()); + + fvph_hit_station_delta_z[iSt]->Fill(fHitQaData.GetHitZ() - fpDetInterface->GetZ(iSt)); + + fvph_hit_dx[iSt]->Fill(fHitQaData.GetHitDx()); + fvph_hit_dy[iSt]->Fill(fHitQaData.GetHitDy()); + fvph_hit_du[iSt]->Fill(fHitQaData.GetHitDu()); + fvph_hit_dv[iSt]->Fill(fHitQaData.GetHitDv()); + fvph_hit_kuv[iSt]->Fill(fHitQaData.GetHitRuv()); + fvph_hit_dt[iSt]->Fill(fHitQaData.GetHitTimeError()); + + // Sum over station distributions + fvph_hit_xy[nSt]->Fill(fHitQaData.GetHitX(), fHitQaData.GetHitY()); + fvph_hit_zx[nSt]->Fill(fHitQaData.GetHitZ(), fHitQaData.GetHitX()); + fvph_hit_zy[nSt]->Fill(fHitQaData.GetHitZ(), fHitQaData.GetHitY()); + + fvph_hit_station_delta_z[nSt]->Fill(fHitQaData.GetHitZ() - fpDetInterface->GetZ(iSt)); + + fvph_hit_dx[nSt]->Fill(fHitQaData.GetHitDx()); + fvph_hit_dy[nSt]->Fill(fHitQaData.GetHitDy()); + fvph_hit_du[nSt]->Fill(fHitQaData.GetHitDu()); + fvph_hit_dv[nSt]->Fill(fHitQaData.GetHitDv()); + fvph_hit_kuv[nSt]->Fill(fHitQaData.GetHitRuv()); + fvph_hit_dt[nSt]->Fill(fHitQaData.GetHitTimeError()); + + // ********************** + // ** MC information QA + + if (IsMCUsed()) { + const auto* pHitMatch = dynamic_cast<CbmMatch*>(fpHitMatches->At(iH)); + assert(pHitMatch); + + // Evaluate number of hits per one MC point + int nMCpoints = 0; // Number of MC points for this hit + for (int iLink = 0; iLink < pHitMatch->GetNofLinks(); ++iLink) { + const auto& link = pHitMatch->GetLink(iLink); + + int iP = link.GetIndex(); // Index of MC point + + // Skip noisy links + if (iP < 0) { continue; } + + ++nMCpoints; + + int iE = fpMCEventList->GetEventIndex(link); + + 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 << ')'; + + // 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; + + 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 << ')'; + + vNofHitsPerMcTrack[iE][iSt][iTr]++; + } + + fvph_n_points_per_hit[iSt]->Fill(nMCpoints); + fvph_n_points_per_hit[nSt]->Fill(nMCpoints); + + // + // Fill the following histograms exclusively for isolated hits with a one-to-one correspondence to the mc track + // The mc track must satisfy the cuts + // + + if (pHitMatch->GetNofLinks() != 1) { continue; } + + // The best link to in the match (probably, the cut on nMCpoints is meaningless) + assert(pHitMatch->GetNofLinks() > 0); // Should be always true due to the cut above + const auto& bestPointLink = pHitMatch->GetMatchedLink(); + + // Skip noisy links + if (bestPointLink.GetIndex() < 0) { continue; } + + // 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; + + // 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: "; + + double t0MC = fpMCEventList->GetEventTime(bestPointLink); + LOG_IF(fatal, t0MC < 0) << fName << ": MC time zero is lower then 0 ns: " << t0MC; + + // cut on the mc track quality + if (!IsTrackSelected(pMCTrack, pMCPoint)) { continue; } + + { // skip the case when the mc point participates to several hits + int iE = fpMCEventList->GetEventIndex(bestTrackLink); + if (vNofHitsPerMcTrack[iE][iSt][pMCPoint->GetTrackID()] != 1) { continue; } + } + + // ----- MC point properties + // + double mass = pMCTrack->GetMass(); + //double charge = pMCTrack->GetCharge(); + //double pdg = pMCTrack->GetPdgCode(); + + // Entrance position and time + // NOTE: SZh 04.09.2023: Methods GetX(), GetY(), GetZ() for MVD, STS, MUCH, TRD and TOF always return + // positions of track at entrance to the active volume. + double xMC = pMCPoint->FairMCPoint::GetX(); + double yMC = pMCPoint->FairMCPoint::GetY(); + double zMC = pMCPoint->FairMCPoint::GetZ(); + double tMC = pMCPoint->GetTime() + t0MC; + + // MC point entrance momenta + // NOTE: SZh 04.09.2023: Methods GetPx(), GetPy(), GetPz() for MVD, STS, MUCH, TRD and TOF always return + // the momentum components of track at entrance to the active volume. + double pxMC = pMCPoint->GetPx(); + double pyMC = pMCPoint->GetPy(); + double pzMC = pMCPoint->GetPz(); + double pMC = sqrt(pxMC * pxMC + pyMC * pyMC + pzMC * pzMC); + + // MC point exit momenta + // double pxMCo = pMCPoint->GetPxOut(); + // double pyMCo = pMCPoint->GetPyOut(); + // double pzMCo = pMCPoint->GetPzOut(); + // double pMCo = sqrt(pxMCo * pxMCo + pyMCo * pyMCo + pzMCo * pzMCo); + + // Position and time shifted to z-coordinate of the hit + // TODO: SZh 04.09.2023: Probably, the approximation + double shiftZ = fHitQaData.GetHitZ() - zMC; // Difference between hit and point z positions + double xMCs = xMC + shiftZ * pxMC / pzMC; + double yMCs = yMC + shiftZ * pyMC / pzMC; + double tMCs = tMC + shiftZ / (pzMC * phys::kSpeedOfLight) * sqrt(mass * mass + pMC * pMC); + + fHitQaData.SetPointTime(tMCs); + fHitQaData.SetPointX(xMCs); + fHitQaData.SetPointY(yMCs); + fHitQaData.SetPointZ(fHitQaData.GetHitZ()); + + double zRes = kNAN; + if constexpr (L1DetectorID::kTof == DetID) { + zRes = fHitQaData.GetHitZ() - pMCPoint->GetZ(); + } + else { + zRes = fHitQaData.GetHitZ() - 0.5 * (pMCPoint->GetZ() + pMCPoint->GetZOut()); + } + fvph_point_hit_delta_z[iSt]->Fill(zRes); + + double xRes = fHitQaData.GetResidualX(); + double yRes = fHitQaData.GetResidualY(); + double uRes = fHitQaData.GetResidualU(); + double vRes = fHitQaData.GetResidualV(); + double tRes = fHitQaData.GetResidualTime(); + + double xPull = fHitQaData.GetPullX(); + double yPull = fHitQaData.GetPullY(); + double uPull = fHitQaData.GetPullU(); + double vPull = fHitQaData.GetPullV(); + double tPull = fHitQaData.GetPullTime(); + + fvph_res_x[iSt]->Fill(xRes); + fvph_res_y[iSt]->Fill(yRes); + fvph_res_u[iSt]->Fill(uRes); + fvph_res_v[iSt]->Fill(vRes); + fvph_res_t[iSt]->Fill(tRes); + + fvph_pull_x[iSt]->Fill(xPull); + fvph_pull_y[iSt]->Fill(yPull); + fvph_pull_u[iSt]->Fill(uPull); + fvph_pull_v[iSt]->Fill(vPull); + fvph_pull_t[iSt]->Fill(tPull); + + fvph_res_x_vs_x[iSt]->Fill(fHitQaData.GetPointX(), xRes); + fvph_res_y_vs_y[iSt]->Fill(fHitQaData.GetPointY(), yRes); + fvph_res_u_vs_u[iSt]->Fill(fHitQaData.GetPointU(), uRes); + fvph_res_v_vs_v[iSt]->Fill(fHitQaData.GetPointV(), vRes); + fvph_res_t_vs_t[iSt]->Fill(fHitQaData.GetPointTime(), tRes); + + fvph_pull_x_vs_x[iSt]->Fill(fHitQaData.GetPointX(), xPull); + fvph_pull_y_vs_y[iSt]->Fill(fHitQaData.GetPointY(), yPull); + fvph_pull_u_vs_u[iSt]->Fill(fHitQaData.GetPointU(), uPull); + fvph_pull_v_vs_v[iSt]->Fill(fHitQaData.GetPointV(), vPull); + fvph_pull_t_vs_t[iSt]->Fill(fHitQaData.GetPointTime(), tPull); + + // fill summary histos + + fvph_point_hit_delta_z[nSt]->Fill(zRes); + + fvph_res_x[nSt]->Fill(xRes); + fvph_res_y[nSt]->Fill(yRes); + fvph_res_u[nSt]->Fill(uRes); + fvph_res_v[nSt]->Fill(vRes); + fvph_res_t[nSt]->Fill(tRes); + + fvph_pull_x[nSt]->Fill(xPull); + fvph_pull_y[nSt]->Fill(yPull); + fvph_pull_u[nSt]->Fill(uPull); + fvph_pull_v[nSt]->Fill(vPull); + fvph_pull_t[nSt]->Fill(tPull); + + fvph_res_x_vs_x[nSt]->Fill(fHitQaData.GetPointX(), xRes); + fvph_res_y_vs_y[nSt]->Fill(fHitQaData.GetPointY(), yRes); + fvph_res_u_vs_u[nSt]->Fill(fHitQaData.GetPointU(), uRes); + fvph_res_v_vs_v[nSt]->Fill(fHitQaData.GetPointV(), vRes); + fvph_res_t_vs_t[nSt]->Fill(fHitQaData.GetPointTime(), tRes); + + fvph_pull_x_vs_x[nSt]->Fill(fHitQaData.GetPointX(), xPull); + fvph_pull_y_vs_y[nSt]->Fill(fHitQaData.GetPointY(), yPull); + fvph_pull_u_vs_u[nSt]->Fill(fHitQaData.GetPointU(), uPull); + fvph_pull_v_vs_v[nSt]->Fill(fHitQaData.GetPointV(), vPull); + fvph_pull_t_vs_t[nSt]->Fill(fHitQaData.GetPointTime(), tPull); + + } + FillHistogramsPerHit(); + } // loop over hits: end + + // Fill hit reconstruction efficiencies + if (IsMCUsed()) { + for (int iE = 0; iE < nMCevents; ++iE) { + int iFile = fpMCEventList->GetFileIdByIndex(iE); + int iEvent = fpMCEventList->GetEventIdByIndex(iE); + int nPoints = fpMCPoints->Size(iFile, iEvent); + int nTracks = fpMCTracks->Size(iFile, iEvent); + + // If efficiency for the track at the station is already evaluated + std::vector<std::vector<bool>> vIsTrackProcessed(nSt); + for (int iSt = 0; iSt < nSt; iSt++) { + vIsTrackProcessed[iSt].resize(nTracks, 0); + } + + for (int iP = 0; iP < nPoints; ++iP) { + fHitQaData.Reset(); + fHitQaData.SetPointID(iP, iEvent, iFile); + + 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; + + int address = pMCPoint->GetDetectorID(); + int iSt = fpDetInterface->GetTrackingStationIndex(address); + LOG_IF(fatal, iSt < 0 || iSt >= nSt) + << fName << ": MC point for FEI = " << iFile << ", " << iEvent << ", " << iP << " and address " << address + << " has wrong station index: iSt = " << iSt; + + // NOTE: SZh 04.09.2023: Methods GetX(), GetY(), GetZ() for MVD, STS, MUCH, TRD and TOF always return + // positions of track at entrance to the active volume. + fHitQaData.SetPointX(pMCPoint->FairMCPoint::GetX()); + fHitQaData.SetPointY(pMCPoint->FairMCPoint::GetY()); + fHitQaData.SetPointZ(pMCPoint->FairMCPoint::GetZ()); + + fvph_point_xy[iSt]->Fill(fHitQaData.GetPointX(), fHitQaData.GetPointY()); + fvph_point_zx[iSt]->Fill(fHitQaData.GetPointZ(), fHitQaData.GetPointX()); + fvph_point_zy[iSt]->Fill(fHitQaData.GetPointZ(), fHitQaData.GetPointY()); + + fvph_point_xy[nSt]->Fill(fHitQaData.GetPointX(), fHitQaData.GetPointY()); + fvph_point_zx[nSt]->Fill(fHitQaData.GetPointZ(), fHitQaData.GetPointX()); + fvph_point_zy[nSt]->Fill(fHitQaData.GetPointZ(), fHitQaData.GetPointY()); + + 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 << ')'; + + 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; + + // check efficiency only once per mc track per station + if (vIsTrackProcessed[iSt][iTr]) { continue; } + + vIsTrackProcessed[iSt][iTr] = true; + + // cut on the mc track quality + if (!fHitQaData.GetIfTrackSelected()) { continue; } + + // Conditions under which point is accounted as reconstructed: point + bool ifTrackHasHits = fHitQaData.GetIfTrackSelected(); + + fvpe_reco_eff_vs_xy[iSt]->Fill(fHitQaData.GetPointX(), fHitQaData.GetPointY(), ifTrackHasHits); + fvpe_reco_eff_vs_xy[nSt]->Fill(fHitQaData.GetPointX(), fHitQaData.GetPointY(), ifTrackHasHits); + } // loop over MC-points: end + } // loop over MC-events: end + } // MC is used: end +} + +// --------------------------------------------------------------------------------------------------------------------- +// +template <L1DetectorID DetID> +InitStatus CbmCaInputQaBase<DetID>::InitDataBranches() +{ + LOG_IF(fatal, !fpDetInterface) << "\033[1;31m" << fName << ": tracking detector interface is undefined\033[0m"; + + // FairRootManager + auto* pFairRootManager = FairRootManager::Instance(); + LOG_IF(fatal, !pFairRootManager) << "\033[1;31m" << fName << ": FairRootManager instance is a null pointer\033[0m"; + + // Time-slice + fpTimeSlice = dynamic_cast<CbmTimeSlice*>(pFairRootManager->GetObject("TimeSlice.")); + LOG_IF(fatal, !fpTimeSlice) << "\033[1;31m" << fName << ": time-slice branch is not found\033[0m"; + + // ----- Reconstructed data-branches initialization + + // Hits container + fpHits = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject(cbm::ca::kDetHitBrName[DetID])); + LOG_IF(fatal, !fpHits) << "\033[1;31m" << fName << ": container of reconstructed hits is not found\033[0m"; + + // ----- MC-information branches initialization + if (IsMCUsed()) { + // MC manager + fpMCDataManager = dynamic_cast<CbmMCDataManager*>(pFairRootManager->GetObject("MCDataManager")); + LOG_IF(fatal, !fpMCDataManager) << "\033[1;31m" << fName << ": MC data manager branch is not found\033[0m"; + + // MC event list + fpMCEventList = dynamic_cast<CbmMCEventList*>(pFairRootManager->GetObject("MCEventList.")); + LOG_IF(fatal, !fpMCEventList) << "\033[1;31m" << fName << ": MC event list branch is not found\033[0m"; + + // MC tracks + fpMCTracks = fpMCDataManager->InitBranch("MCTrack"); + LOG_IF(fatal, !fpMCTracks) << "\033[1;31m" << fName << ": MC track branch is not found\033[0m"; + + // MC points + fpMCPoints = fpMCDataManager->InitBranch(cbm::ca::kDetPointBrName[DetID]); + LOG_IF(fatal, !fpMCTracks) << "\033[1;31m" << fName << ": MC point branch is not found\033[0m"; + + // Hit matches + const char* hitMatchName = Form("%sMatch", cbm::ca::kDetHitBrName[DetID]); + fpHitMatches = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject(hitMatchName)); + LOG_IF(fatal, !fpHitMatches) << "\033[1;31m]" << fName << ": hit match branch is not found\033[0m"; + } + + // get hit distribution ranges from the geometry + int nSt = fpDetInterface->GetNtrackingStations(); + + frXmin.resize(nSt + 1, 0.); + frXmax.resize(nSt + 1, 0.); + + frYmin.resize(nSt + 1, 0.); + frYmax.resize(nSt + 1, 0.); + + frZmin.resize(nSt + 1, 0.); + frZmax.resize(nSt + 1, 0.); + + for (int i = nSt; i >= 0; --i) { + { // read the ranges for station i + int j = (i == nSt ? 0 : i); + + frXmin[i] = fpDetInterface->GetXmin(j); + frXmax[i] = fpDetInterface->GetXmax(j); + + frYmin[i] = fpDetInterface->GetYmin(j); + frYmax[i] = fpDetInterface->GetYmax(j); + + frZmin[i] = fpDetInterface->GetZmin(j); + frZmax[i] = fpDetInterface->GetZmax(j); + } + + // update overall ranges + + frXmin[nSt] = std::min(frXmin[nSt], frXmin[i]); + frXmax[nSt] = std::max(frXmax[nSt], frXmax[i]); + + frYmin[nSt] = std::min(frYmin[nSt], frYmin[i]); + frYmax[nSt] = std::max(frYmax[nSt], frYmax[i]); + + frZmin[nSt] = std::min(frZmin[nSt], frZmin[i]); + frZmax[nSt] = std::max(frZmax[nSt], frZmax[i]); + } + + // add 5% margins for a better view + + for (int i = 0; i <= nSt; ++i) { + double dx = 0.05 * fabs(frXmax[i] - frXmin[i]); + frXmin[i] -= dx; + frXmax[i] += dx; + + double dy = 0.05 * fabs(frYmax[i] - frYmin[i]); + frYmin[i] -= dy; + frYmax[i] += dy; + + double dz = 0.05 * fabs(frZmax[i] - frZmin[i]); + frZmin[i] -= dz; + frZmax[i] += dz; + } + + + return kSUCCESS; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +template <L1DetectorID DetID> +InitStatus CbmCaInputQaBase<DetID>::InitHistograms() +{ + int nSt = fpDetInterface->GetNtrackingStations(); + + //SetStoringMode(EStoringMode::kSAMEDIR); + SetStoringMode(EStoringMode::kSUBDIR); + + // ----- Subdirectories + MakeQaDirectory("Summary"); + MakeQaDirectory("Summary/vs Station"); + MakeQaDirectory("Summary/vs N digi"); + MakeQaDirectory("All stations"); + + // ----- Histograms initialization + fvph_hit_xy.resize(nSt + 1, nullptr); + fvph_hit_zy.resize(nSt + 1, nullptr); + fvph_hit_zx.resize(nSt + 1, nullptr); + + fvph_hit_station_delta_z.resize(nSt + 1, nullptr); + + fvph_hit_dx.resize(nSt + 1, nullptr); + fvph_hit_dy.resize(nSt + 1, nullptr); + fvph_hit_dt.resize(nSt + 1, nullptr); + fvph_hit_dv.resize(nSt + 1, nullptr); + fvph_hit_kuv.resize(nSt + 1, nullptr); + fvph_hit_du.resize(nSt + 1, nullptr); + + std::string detName = fpDetInterface->GetDetectorName(); + + for (int iSt = 0; iSt <= nSt; ++iSt) { + TString sF = ""; // Subfolder + sF += (iSt == nSt) ? "All stations/" : Form("Station %d/", iSt); + TString nsuff = (iSt == nSt) ? "" : Form("_st%d", iSt); // Histogram name suffix + TString tsuff = (iSt == nSt) ? "" : Form(" in %s station %d", detName.c_str(), iSt); // Histogram title suffix + TString sN = ""; + TString sT = ""; + + // place directories in a prefered order + MakeQaDirectory(sF + "occup/"); + if (IsMCUsed()) { + MakeQaDirectory(sF + "res/"); + MakeQaDirectory(sF + "pull/"); + MakeQaDirectory(sF + "eff/"); + } + MakeQaDirectory(sF + "err/"); + + int nBinsZ = (iSt < nSt) ? fNbins : fNbinsZ; + + // Hit occupancy + sN = (TString) "hit_xy" + nsuff; + sT = (TString) "Hit occupancy in xy-Plane" + tsuff + ";x_{hit} [cm];y_{hit} [cm]"; + fvph_hit_xy[iSt] = + MakeQaObject<TH2F>(sF + "occup/" + sN, sT, fNbins, frXmin[iSt], frXmax[iSt], fNbins, frYmin[iSt], frYmax[iSt]); + + sN = (TString) "hit_zx" + nsuff; + sT = (TString) "Hit occupancy in xz-Plane" + tsuff + ";z_{hit} [cm];x_{hit} [cm]"; + fvph_hit_zx[iSt] = + MakeQaObject<TH2F>(sF + "occup/" + sN, sT, nBinsZ, frZmin[iSt], frZmax[iSt], fNbins, frXmin[iSt], frXmax[iSt]); + + sN = (TString) "hit_zy" + nsuff; + sT = (TString) "Hit occupancy in yz-plane" + tsuff + ";z_{hit} [cm];y_{hit} [cm]"; + fvph_hit_zy[iSt] = + MakeQaObject<TH2F>(sF + "occup/" + sN, sT, nBinsZ, frZmin[iSt], frZmax[iSt], fNbins, frYmin[iSt], frYmax[iSt]); + + // Hit errors + sN = (TString) "hit_dx" + nsuff; + sT = (TString) "Hit position error along x-axis" + tsuff + ";dx_{hit} [cm]"; + fvph_hit_dx[iSt] = MakeQaObject<TH1F>(sF + "err/" + sN, sT, fNbins, fRHitDx[0], fRHitDx[1]); + + sN = (TString) "hit_dy" + nsuff; + sT = (TString) "Hit position error along y-axis" + tsuff + ";dy_{hit} [cm]"; + fvph_hit_dy[iSt] = MakeQaObject<TH1F>(sF + "err/" + sN, sT, fNbins, fRHitDy[0], fRHitDy[1]); + + sN = (TString) "hit_du" + nsuff; + sT = (TString) "Hit position error across front strips" + tsuff + ";du_{hit} [cm]"; + fvph_hit_du[iSt] = MakeQaObject<TH1F>(sF + "err/" + sN, sT, fNbins, fRHitDu[0], fRHitDu[1]); + + sN = (TString) "hit_dv" + nsuff; + sT = (TString) "Hit position error across back strips" + tsuff + ";dv_{hit} [cm]"; + fvph_hit_dv[iSt] = MakeQaObject<TH1F>(sF + "err/" + sN, sT, fNbins, fRHitDv[0], fRHitDv[1]); + + sN = (TString) "hit_kuv" + nsuff; + sT = (TString) "Hit error correlation between front and back strips" + tsuff + ";kuv_{hit} [unitless]"; + fvph_hit_kuv[iSt] = MakeQaObject<TH1F>(sF + "err/" + sN, sT, fNbins / 2 * 2 + 1, -1.1, 1.1); + + sN = (TString) "hit_dt" + nsuff; + sT = (TString) "Hit time error" + tsuff + ";dt_{hit} [ns]"; + fvph_hit_dt[iSt] = MakeQaObject<TH1F>(sF + "err/" + sN, sT, fNbins, fRHitDt[0], fRHitDt[1]); + + sN = (TString) "hit_station_delta_z" + nsuff; + sT = (TString) "Different between hit and station z-positions" + tsuff + ";z_{hit} - z_{st} [cm]"; + fvph_hit_station_delta_z[iSt] = MakeQaObject<TH1F>(sF + sN, sT, fNbins, -5., 5.); + + } // loop over station index: end + + // ----- Initialize histograms, which are use MC-information + if (IsMCUsed()) { + // Resize histogram vectors + fvph_n_points_per_hit.resize(nSt + 1, nullptr); + fvph_point_xy.resize(nSt + 1, nullptr); + fvph_point_zx.resize(nSt + 1, nullptr); + fvph_point_zy.resize(nSt + 1, nullptr); + fvph_point_hit_delta_z.resize(nSt + 1, nullptr); + fvph_res_x.resize(nSt + 1, nullptr); + fvph_res_y.resize(nSt + 1, nullptr); + fvph_res_u.resize(nSt + 1, nullptr); + fvph_res_v.resize(nSt + 1, nullptr); + fvph_res_t.resize(nSt + 1, nullptr); + fvph_pull_x.resize(nSt + 1, nullptr); + fvph_pull_y.resize(nSt + 1, nullptr); + fvph_pull_u.resize(nSt + 1, nullptr); + fvph_pull_v.resize(nSt + 1, nullptr); + fvph_pull_t.resize(nSt + 1, nullptr); + fvph_res_x_vs_x.resize(nSt + 1, nullptr); + fvph_res_y_vs_y.resize(nSt + 1, nullptr); + fvph_res_u_vs_u.resize(nSt + 1, nullptr); + fvph_res_v_vs_v.resize(nSt + 1, nullptr); + fvph_res_t_vs_t.resize(nSt + 1, nullptr); + fvph_pull_x_vs_x.resize(nSt + 1, nullptr); + fvph_pull_y_vs_y.resize(nSt + 1, nullptr); + fvph_pull_u_vs_u.resize(nSt + 1, nullptr); + fvph_pull_v_vs_v.resize(nSt + 1, nullptr); + fvph_pull_t_vs_t.resize(nSt + 1, nullptr); + fvpe_reco_eff_vs_xy.resize(nSt + 1, nullptr); + fvph_reco_eff.resize(nSt + 1, nullptr); + + for (int iSt = 0; iSt <= nSt; ++iSt) { + TString sF = ""; // Subfolder + sF += (iSt == nSt) ? "All stations/" : Form("Station %d/", iSt); + TString nsuff = (iSt == nSt) ? "" : Form("_st%d", iSt); // Histogram name suffix + TString tsuff = (iSt == nSt) ? "" : Form(" in %s station %d", detName.c_str(), iSt); // Histogram title suffix + TString sN = ""; + TString sT = ""; + + sN = (TString) "n_points_per_hit" + nsuff; + sT = (TString) "Number of points per hit" + tsuff + ";N_{point}/hit"; + fvph_n_points_per_hit[iSt] = MakeQaObject<TH1F>(sF + sN, sT, 10, -0.5, 9.5); + + // Point occupancy + sN = (TString) "point_xy" + nsuff; + sT = (TString) "Point occupancy in XY plane" + tsuff + ";x_{MC} [cm];y_{MC} [cm]"; + fvph_point_xy[iSt] = + MakeQaObject<TH2F>(sF + "occup/" + sN, sT, fNbins, frXmin[iSt], frXmax[iSt], fNbins, frYmin[iSt], frYmax[iSt]); + + sN = (TString) "point_zx" + nsuff; + sT = (TString) "Point Occupancy in XZ plane" + tsuff + ";z_{MC} [cm];x_{MC} [cm]"; + fvph_point_zx[iSt] = + MakeQaObject<TH2F>(sF + "occup/" + sN, sT, fNbinsZ, frZmin[iSt], frZmax[iSt], fNbins, frXmin[iSt], frXmax[iSt]); + + sN = (TString) "point_zy" + nsuff; + sT = (TString) "Point Occupancy in YZ Plane" + tsuff + ";z_{MC} [cm];y_{MC} [cm]"; + fvph_point_zy[iSt] = + MakeQaObject<TH2F>(sF + "occup/" + sN, sT, fNbinsZ, frZmin[iSt], frZmax[iSt], fNbins, frYmin[iSt], frYmax[iSt]); + + // Difference between MC input z and hit z coordinates + sN = (TString) "point_hit_delta_z" + nsuff; + sT = (TString) "Distance between " + detName + " point and hit along z axis" + tsuff + ";z_{reco} - z_{MC} [cm]"; + fvph_point_hit_delta_z[iSt] = MakeQaObject<TH1F>(sF + sN, sT, fNbins, -0.05, 0.05); + + sN = (TString) "res_x" + nsuff; + sT = (TString) "Residuals for X" + tsuff + ";x_{reco} - x_{MC} [cm]"; + fvph_res_x[iSt] = MakeQaObject<TH1F>(sF + "res/" + sN, sT, fNbins, fRResX[0], fRResX[1]); + + sN = (TString) "res_y" + nsuff; + sT = (TString) "Residuals for Y" + tsuff + ";y_{reco} - y_{MC} [cm]"; + fvph_res_y[iSt] = MakeQaObject<TH1F>(sF + "res/" + sN, sT, fNbins, fRResY[0], fRResY[1]); + + sN = (TString) "res_u" + nsuff; + sT = (TString) "Residuals for Front strip coordinate" + tsuff + ";u_{reco} - u_{MC} [cm]"; + fvph_res_u[iSt] = MakeQaObject<TH1F>(sF + "res/" + sN, sT, fNbins, fRResU[0], fRResU[1]); + + sN = (TString) "res_v" + nsuff; + sT = (TString) "Residuals for Back strip coordinate" + tsuff + ";v_{reco} - v_{MC} [cm]"; + fvph_res_v[iSt] = MakeQaObject<TH1F>(sF + "res/" + sN, sT, fNbins, fRResV[0], fRResV[1]); + + sN = (TString) "res_t" + nsuff; + sT = (TString) "Residuals for Time" + tsuff + ";t_{reco} - t_{MC} [ns]"; + fvph_res_t[iSt] = MakeQaObject<TH1F>(sF + "res/" + sN, sT, fNbins, fRResT[0], fRResT[1]); + + sN = (TString) "pull_x" + nsuff; + sT = (TString) "Pulls for X" + tsuff + ";(x_{reco} - x_{MC}) / #sigma_{x}^{reco}"; + fvph_pull_x[iSt] = MakeQaObject<TH1F>(sF + "pull/" + sN, sT, kNbinsPull, kRPull[0], kRPull[1]); + + sN = (TString) "pull_y" + nsuff; + sT = (TString) "Pulls for Y" + tsuff + ";(y_{reco} - y_{MC}) / #sigma_{y}^{reco}"; + fvph_pull_y[iSt] = MakeQaObject<TH1F>(sF + "pull/" + sN, sT, kNbinsPull, kRPull[0], kRPull[1]); + + sN = (TString) "pull_u" + nsuff; + sT = (TString) "Pulls for Front strip coordinate" + tsuff + ";(u_{reco} - u_{MC}) / #sigma_{u}^{reco}"; + fvph_pull_u[iSt] = MakeQaObject<TH1F>(sF + "pull/" + sN, sT, kNbinsPull, kRPull[0], kRPull[1]); + + sN = (TString) "pull_v" + nsuff; + sT = (TString) "Pulls for Back strip coordinate" + tsuff + ";(v_{reco} - v_{MC}) / #sigma_{v}^{reco}"; + fvph_pull_v[iSt] = MakeQaObject<TH1F>(sF + "pull/" + sN, sT, kNbinsPull, kRPull[0], kRPull[1]); + + sN = (TString) "pull_t" + nsuff; + sT = (TString) "Pulls for Time" + tsuff + ";(t_{reco} - t_{MC}) / #sigma_{t}^{reco}"; + fvph_pull_t[iSt] = MakeQaObject<TH1F>(sF + "pull/" + sN, sT, kNbinsPull, kRPull[0], kRPull[1]); + + sN = (TString) "res_x_vs_x" + nsuff; + sT = (TString) "Residuals for X" + tsuff + ";x_{MC} [cm];x_{reco} - x_{MC} [cm]"; + fvph_res_x_vs_x[iSt] = + MakeQaObject<TH2F>(sF + "res/" + sN, sT, fNbins, frXmin[iSt], frXmax[iSt], fNbins, fRResX[0], fRResX[1]); + + sN = (TString) "res_y_vs_y" + nsuff; + sT = (TString) "Residuals for Y" + tsuff + ";y_{MC} [cm];y_{reco} - y_{MC} [cm]"; + fvph_res_y_vs_y[iSt] = + MakeQaObject<TH2F>(sF + "res/" + sN, sT, fNbins, frYmin[iSt], frYmax[iSt], fNbins, fRResY[0], fRResY[1]); + + sN = (TString) "res_u_vs_u" + nsuff; + sT = (TString) "Residuals for Front strip coordinate" + tsuff + ";u_{MC} [cm];u_{reco} - u_{MC} [cm]"; + fvph_res_u_vs_u[iSt] = + MakeQaObject<TH2F>(sF + "res/" + sN, sT, fNbins, frXmin[iSt], frXmax[iSt], fNbins, fRResU[0], fRResU[1]); + + sN = (TString) "res_v_vs_v" + nsuff; + sT = (TString) "Residuals for Back strip coordinate" + tsuff + ";v_{MC} [cm];v_{reco} - v_{MC} [cm]"; + fvph_res_v_vs_v[iSt] = + MakeQaObject<TH2F>(sF + "res/" + sN, sT, fNbins, frXmin[iSt], frXmax[iSt], fNbins, fRResV[0], fRResV[1]); + + sN = (TString) "res_t_vs_t" + nsuff; + sT = (TString) "Residuals for Time" + tsuff + ";t_{MC} [ns];t_{reco} - t_{MC} [ns]"; + fvph_res_t_vs_t[iSt] = MakeQaObject<TH2F>(sF + "res/" + sN, sT, fNbins, 0, 0, fNbins, fRResT[0], fRResT[1]); + + sN = (TString) "pull_x_vs_x" + nsuff; + sT = (TString) "Pulls for X" + tsuff + ";x_{MC} [cm];(x_{reco} - x_{MC}) / #sigma_{x}^{reco}"; + fvph_pull_x_vs_x[iSt] = + MakeQaObject<TH2F>(sF + "pull/" + sN, sT, fNbins, frXmin[iSt], frXmax[iSt], kNbinsPull, kRPull[0], kRPull[1]); + + sN = (TString) "pull_y_vs_y" + nsuff; + sT = (TString) "Pulls for Y" + tsuff + ";y_{MC} [cm];(y_{reco} - y_{MC}) / #sigma_{y}^{reco}"; + fvph_pull_y_vs_y[iSt] = + MakeQaObject<TH2F>(sF + "pull/" + sN, sT, fNbins, frYmin[iSt], frYmax[iSt], kNbinsPull, kRPull[0], kRPull[1]); + + sN = (TString) "pull_u_vs_u" + nsuff; + sT = + (TString) "Pulls for Front strip coordinate" + tsuff + ";u_{MC} [cm];(u_{reco} - u_{MC}) / #sigma_{u}^{reco}"; + fvph_pull_u_vs_u[iSt] = + MakeQaObject<TH2F>(sF + "pull/" + sN, sT, fNbins, frXmin[iSt], frXmax[iSt], kNbinsPull, kRPull[0], kRPull[1]); + + sN = (TString) "pull_v_vs_v" + nsuff; + sT = + (TString) "Pulls for Back strip coordinate" + tsuff + ";v_{MC} [cm];(v_{reco} - v_{MC}) / #sigma_{v}^{reco}"; + fvph_pull_v_vs_v[iSt] = + MakeQaObject<TH2F>(sF + "pull/" + sN, sT, fNbins, frXmin[iSt], frXmax[iSt], kNbinsPull, kRPull[0], kRPull[1]); + + sN = (TString) "pull_t_vs_t" + nsuff; + sT = (TString) "Pulls for Time" + tsuff + ";t_{MC} [ns];(t_{reco} - t_{MC}) / #sigma_{t}^{reco}"; + fvph_pull_t_vs_t[iSt] = MakeQaObject<TH2F>(sF + "pull/" + sN, sT, fNbins, 0, 0, fNbins, kRPull[0], kRPull[1]); + + sN = (TString) "reco_eff_vs_xy" + nsuff; + sT = (TString) "Hit rec. efficiency in XY" + tsuff + ";x_{MC} [cm];y_{MC} [cm];#epsilon"; + fvpe_reco_eff_vs_xy[iSt] = MakeQaObject<TProfile2D>(sF + "eff/" + sN, sT, fNbins, frXmin[iSt], frXmax[iSt], + fNbins, frYmin[iSt], frYmax[iSt]); + + sN = (TString) "reco_eff" + nsuff; + sT = (TString) "Hit rec. efficiency in XY bins" + tsuff + ";eff"; + fvph_reco_eff[iSt] = MakeQaObject<TH1F>(sF + "eff/" + sN, sT, 130, -0.005, 1.30 - 0.005); + } // stations + } + return kSUCCESS; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +// TODO: rename the method to Finish or Finalise or MakeSummary, since it does more than just canvases +// +template <L1DetectorID DetID> +InitStatus CbmCaInputQaBase<DetID>::InitCanvases() +{ + gStyle->SetOptFit(1); + int nSt = fpDetInterface->GetNtrackingStations(); + + // ******************** + // ** Drawing options + + // Contours + constexpr auto contColor = kOrange + 7; + constexpr auto contWidth = 2; // Line width + constexpr auto contStyle = 2; // Line style + constexpr auto contFill = 0; // Fill style + + // ********************** + // ** summary for all stations + + { + auto* canv = MakeQaObject<CbmQaCanvas>("occ_hit", "Hit Occupancy", 1600, 800); + canv->Divide2D(3); + canv->cd(1); + fvph_hit_xy[nSt]->DrawCopy("colz", ""); + canv->cd(2); + fvph_hit_zx[nSt]->DrawCopy("colz", ""); + canv->cd(3); + fvph_hit_zy[nSt]->DrawCopy("colz", ""); + } + + { + auto* canv = MakeQaObject<CbmQaCanvas>("occ_point", "Point Occupancy", 1600, 800); + canv->Divide2D(3); + + canv->cd(1); + fvph_point_xy[nSt]->DrawCopy("colz", ""); + canv->cd(2); + fvph_point_zx[nSt]->DrawCopy("colz", ""); + canv->cd(3); + fvph_point_zy[nSt]->DrawCopy("colz", ""); + } + + { + auto* canv = MakeQaObject<CbmQaCanvas>("residual", "Hit Residuals", 1600, 800); + canv->Divide2D(6); + canv->cd(1); + fvph_res_x[nSt]->DrawCopy("colz", ""); + canv->cd(2); + fvph_res_y[nSt]->DrawCopy("colz", ""); + canv->cd(3); + fvph_res_t[nSt]->DrawCopy("colz", ""); + canv->cd(4); + fvph_res_u[nSt]->DrawCopy("colz", ""); + canv->cd(5); + fvph_res_v[nSt]->DrawCopy("colz", ""); + } + + { + auto* canv = MakeQaObject<CbmQaCanvas>("pull", "Hit Pulls", 1600, 800); + canv->Divide2D(6); + canv->cd(1); + fvph_pull_x[nSt]->DrawCopy("colz", ""); + canv->cd(2); + fvph_pull_y[nSt]->DrawCopy("colz", ""); + canv->cd(3); + fvph_pull_t[nSt]->DrawCopy("colz", ""); + canv->cd(4); + fvph_pull_u[nSt]->DrawCopy("colz", ""); + canv->cd(5); + fvph_pull_v[nSt]->DrawCopy("colz", ""); + } + + { + auto* canv = MakeQaObject<CbmQaCanvas>("eff", "Hit Reconstruction Efficiency", 1600, 800); + canv->Divide2D(2); + canv->cd(1); + fvpe_reco_eff_vs_xy[nSt]->DrawCopy("colz", ""); + canv->cd(2); + fvph_reco_eff[nSt]->DrawCopy("colz", ""); + } + + { + auto* canv = MakeQaObject<CbmQaCanvas>("err", "Hit Errors", 1600, 800); + canv->Divide2D(6); + canv->cd(1); + fvph_hit_dx[nSt]->DrawCopy(); + canv->cd(2); + fvph_hit_dy[nSt]->DrawCopy(); + canv->cd(3); + fvph_hit_dt[nSt]->DrawCopy(); + canv->cd(4); + fvph_hit_du[nSt]->DrawCopy(); + canv->cd(5); + fvph_hit_dv[nSt]->DrawCopy(); + canv->cd(6); + fvph_hit_kuv[nSt]->DrawCopy(); + } + + { + auto* canv = MakeQaObject<CbmQaCanvas>("other", "Other histograms", 1600, 800); + canv->Divide2D(3); + canv->cd(1); + fvph_n_points_per_hit[nSt]->DrawCopy("colz", ""); + canv->cd(2); + fvph_hit_station_delta_z[nSt]->DrawCopy("colz", ""); + canv->cd(3); + fvph_point_hit_delta_z[nSt]->DrawCopy("colz", ""); + } + + + // ********************************* + // ** Hit occupancies + + { // ** Occupancy in XY plane + auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/hit_xy", "", 1600, 800); + canv->Divide2D(nSt); + for (int iSt = 0; iSt < nSt; ++iSt) { + canv->cd(iSt + 1); + fvph_hit_xy[iSt]->DrawCopy("colz", ""); + + // Square boarders of the active area of the station + double stXmin = fpDetInterface->GetXmin(iSt); + double stXmax = fpDetInterface->GetXmax(iSt); + double stYmin = fpDetInterface->GetYmin(iSt); + double stYmax = fpDetInterface->GetYmax(iSt); + + auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax); + pBox->SetLineWidth(contWidth); + pBox->SetLineStyle(contStyle); + pBox->SetLineColor(contColor); + pBox->SetFillStyle(contFill); + pBox->Draw("SAME"); + } + } + + { // ** Occupancy in XZ plane + auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/hit_zx", "", 1600, 800); + canv->cd(); + fvph_hit_zx[nSt]->DrawCopy("colz", ""); + for (int iSt = 0; iSt < nSt; ++iSt) { + // Station positions in detector IFS + double stZmin = fpDetInterface->GetZmin(iSt); + double stZmax = fpDetInterface->GetZmax(iSt); + double stHmin = fpDetInterface->GetXmin(iSt); + double stHmax = fpDetInterface->GetXmax(iSt); + + auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax); + pBox->SetLineWidth(contWidth); + pBox->SetLineStyle(contStyle); + pBox->SetLineColor(contColor); + pBox->SetFillStyle(contFill); + pBox->Draw("SAME"); + } + } + + { // ** Occupancy in YZ plane + auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/hit_zy", "", 1600, 800); + canv->cd(); + fvph_hit_zy[nSt]->DrawCopy("colz", ""); + for (int iSt = 0; iSt < nSt; ++iSt) { + // Station positions in detector IFS + double stZmin = fpDetInterface->GetZmin(iSt); + double stZmax = fpDetInterface->GetZmax(iSt); + double stHmin = fpDetInterface->GetYmin(iSt); + double stHmax = fpDetInterface->GetYmax(iSt); + + auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax); + pBox->SetLineWidth(contWidth); + pBox->SetLineStyle(contStyle); + pBox->SetLineColor(contColor); + pBox->SetFillStyle(contFill); + pBox->Draw("SAME"); + } + } + + // ************ + // ************ + + if (IsMCUsed()) { + + // ********************** + // ** Point occupancies + + { // ** Occupancy in XY plane + auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/point_xy", "", 1600, 800); + canv->Divide2D(nSt); + for (int iSt = 0; iSt < nSt; ++iSt) { + canv->cd(iSt + 1); + fvph_point_xy[iSt]->DrawCopy("colz", ""); + + // Square boarders of the active area of the station + double stXmin = fpDetInterface->GetXmin(iSt); + double stXmax = fpDetInterface->GetXmax(iSt); + double stYmin = fpDetInterface->GetYmin(iSt); + double stYmax = fpDetInterface->GetYmax(iSt); + + auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax); + pBox->SetLineWidth(contWidth); + pBox->SetLineStyle(contStyle); + pBox->SetLineColor(contColor); + pBox->SetFillStyle(contFill); + pBox->Draw("SAME"); + } + } + + { // ** Occupancy in XZ plane + auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/point_zx", "", 1600, 800); + canv->cd(); + fvph_point_zx[nSt]->DrawCopy("colz", ""); + for (int iSt = 0; iSt < nSt; ++iSt) { + // Station positions in detector IFS + double stZmin = fpDetInterface->GetZmin(iSt); + double stZmax = fpDetInterface->GetZmax(iSt); + double stHmin = fpDetInterface->GetXmin(iSt); + double stHmax = fpDetInterface->GetXmax(iSt); + + auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax); + pBox->SetLineWidth(contWidth); + pBox->SetLineStyle(contStyle); + pBox->SetLineColor(contColor); + pBox->SetFillStyle(contFill); + pBox->Draw("SAME"); + } + } + + { // ** Occupancy in YZ plane + auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/point_zy", "", 1600, 800); + canv->cd(); + fvph_point_zy[nSt]->DrawCopy("colz", ""); + for (int iSt = 0; iSt < nSt; ++iSt) { + // Station positions in detector IFS + double stZmin = fpDetInterface->GetZmin(iSt); + double stZmax = fpDetInterface->GetZmax(iSt); + double stHmin = fpDetInterface->GetYmin(iSt); + double stHmax = fpDetInterface->GetYmax(iSt); + + auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax); + pBox->SetLineWidth(contWidth); + pBox->SetLineStyle(contStyle); + pBox->SetLineColor(contColor); + pBox->SetFillStyle(contFill); + pBox->Draw("SAME"); + } + } + + // ************** + // ** Residuals + + { // x-coordinate + auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/res_x", "Residuals for x coordinate", 1600, 800); + canv->Divide2D(nSt); + for (int iSt = 0; iSt < nSt; ++iSt) { + canv->cd(iSt + 1); + fvph_res_x[iSt]->DrawCopy("", ""); + } + } + + { // y-coordinate + auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/res_y", "Residuals for y coordinate", 1600, 800); + canv->Divide2D(nSt); + for (int iSt = 0; iSt < nSt; ++iSt) { + canv->cd(iSt + 1); + fvph_res_y[iSt]->DrawCopy("", ""); + } + } + + { // u-coordinate + auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/res_u", "Residuals for u coordinate", 1600, 800); + canv->Divide2D(nSt); + for (int iSt = 0; iSt < nSt; ++iSt) { + canv->cd(iSt + 1); + fvph_res_u[iSt]->DrawCopy("", ""); + } + } + + { // v-coordinate + auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/res_v", "Residuals for v coordinate", 1600, 800); + canv->Divide2D(nSt); + for (int iSt = 0; iSt < nSt; ++iSt) { + canv->cd(iSt + 1); + fvph_res_v[iSt]->DrawCopy("", ""); + } + } + + { // t-coordinate + auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/res_t", "Residuals for time", 1600, 800); + canv->Divide2D(nSt); + for (int iSt = 0; iSt < nSt; ++iSt) { + canv->cd(iSt + 1); + fvph_res_t[iSt]->DrawCopy("", ""); + } + } + + + // ********** + // ** Pulls + + { // x-coordinate + auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/pull_x", "Pulls for x coordinate", 1600, 800); + canv->Divide2D(nSt); + for (int iSt = 0; iSt < nSt; ++iSt) { + canv->cd(iSt + 1); + fvph_pull_x[iSt]->DrawCopy("", ""); + } + } + + { // y-coordinate + auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/pull_y", "Pulls for y coordinate", 1600, 800); + canv->Divide2D(nSt); + for (int iSt = 0; iSt < nSt; ++iSt) { + canv->cd(iSt + 1); + fvph_pull_y[iSt]->DrawCopy("", ""); + } + } + + { // u-coordinate + auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/pull_u", "Pulls for u coordinate", 1600, 800); + canv->Divide2D(nSt); + for (int iSt = 0; iSt < nSt; ++iSt) { + canv->cd(iSt + 1); + fvph_pull_u[iSt]->DrawCopy("", ""); + } + } + + { // v-coordinate + auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/pull_v", "Pulls for v coordinate", 1600, 800); + canv->Divide2D(nSt); + for (int iSt = 0; iSt < nSt; ++iSt) { + canv->cd(iSt + 1); + fvph_pull_v[iSt]->DrawCopy("", ""); + } + } + + { // t-coordinate + auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/pull_t", "Pulls for time", 1600, 800); + canv->Divide2D(nSt); + for (int iSt = 0; iSt < nSt; ++iSt) { + canv->cd(iSt + 1); + fvph_pull_t[iSt]->DrawCopy("", ""); + } + } + + // ************************************ + // ** Hit reconstruction efficiencies + + { + auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/reco_eff", "Hit efficiencies in xy bins", 1600, 800); + canv->Divide2D(nSt); + for (int iSt = 0; iSt < nSt; ++iSt) { + canv->cd(iSt + 1); + //fvph_reco_eff[iSt]->SetMarkerStyle(20); + fvph_reco_eff[iSt]->DrawCopy(); + } + } + + { + auto* canv = MakeQaObject<CbmQaCanvas>("vs Station/reco_eff_vs_xy", "Hit efficiencies wrt x and y", 1600, 800); + canv->Divide2D(nSt); + for (int iSt = 0; iSt < nSt; ++iSt) { + canv->cd(iSt + 1); + fvpe_reco_eff_vs_xy[iSt]->DrawCopy("colz", ""); + } + } + } + + return kSUCCESS; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +template <L1DetectorID DetID> +bool CbmCaInputQaBase<DetID>::IsTrackSelected(const CbmMCTrack* track, const Point_t* point) const +{ + assert(point); + + if (fTrackSelectionPrimary && track->GetMotherId() >= 0) { return false; } + + double px = std::numeric_limits<double>::signaling_NaN(); + double py = std::numeric_limits<double>::signaling_NaN(); + double pz = std::numeric_limits<double>::signaling_NaN(); + if constexpr (L1DetectorID::kTof == DetID) { + px = point->GetPx(); + py = point->GetPy(); + pz = point->GetPz(); + } + else { + px = point->GetPxOut(); + py = point->GetPyOut(); + pz = point->GetPzOut(); + } + double p = sqrt(px * px + py * py + pz * pz); + + if (pz <= 0.) { return false; } + + if (p < fTrackSelectionMomentum) { return false; } + + if (TMath::ATan2(sqrt(px * px + py * py), pz) * TMath::RadToDeg() > fTrackSelectionAngle) { return false; } + + return true; +} + +template class CbmCaInputQaBase<L1DetectorID::kMvd>; +template class CbmCaInputQaBase<L1DetectorID::kSts>; +template class CbmCaInputQaBase<L1DetectorID::kMuch>; +template class CbmCaInputQaBase<L1DetectorID::kTrd>; +template class CbmCaInputQaBase<L1DetectorID::kTof>; + + diff --git a/reco/L1/qa/CbmCaInputQaBase.h b/reco/L1/qa/CbmCaInputQaBase.h index f775b093ebc485f04bac68a579569ab852965eba..ea0b4791f9a8f6e0e0ed547a9ac8d2f7f6ce1973 100644 --- a/reco/L1/qa/CbmCaInputQaBase.h +++ b/reco/L1/qa/CbmCaInputQaBase.h @@ -42,6 +42,7 @@ class TProfile2D; /// A QA-task class, which provides assurance of MuCh hits and geometry template <L1DetectorID DetID> class CbmCaInputQaBase : public CbmQaTask { +protected: using Point_t = cbm::ca::PointTypes_t::at<DetID>; ///< Point type for detector ID using Hit_t = cbm::ca::HitTypes_t::at<DetID>; ///< Hit type for detector ID @@ -110,16 +111,16 @@ protected: void FillHistograms(); /// @brief Fills histograms per hit - void FillHistogramsPerHit() {} + virtual void FillHistogramsPerHit() {} /// @brief Fills histograms per MC point - void FillHistogramsPerPoint() {} + virtual void FillHistogramsPerPoint() {} /// @brief De-initializes histograms void DeInit(); /// @brief Defines parameter for a derived class - void DefineParameters(); + virtual void DefineParameters() = 0; /// @brief Checks ranges for mean and standard deviation /// @return False, if variable exits the range diff --git a/reco/L1/qa/CbmCaInputQaMuch.cxx b/reco/L1/qa/CbmCaInputQaMuch.cxx index 1d92dff15f6f45cfe4b1d982d6f6cfae60b4d08f..4ca71307bd5453c939c63361aebe147ed2a8bfc1 100644 --- a/reco/L1/qa/CbmCaInputQaMuch.cxx +++ b/reco/L1/qa/CbmCaInputQaMuch.cxx @@ -8,221 +8,23 @@ /// \author S.Zharko <s.zharko@gsi.de> #include "CbmCaInputQaMuch.h" - -#include "CbmAddress.h" -#include "CbmMCDataArray.h" -#include "CbmMCEventList.h" -#include "CbmMCTrack.h" -#include "CbmMatch.h" -#include "CbmMuchPixelHit.h" -#include "CbmMuchPoint.h" #include "CbmMuchTrackingInterface.h" -#include "CbmQaCanvas.h" -#include "CbmQaEff.h" -#include "CbmTimeSlice.h" - -#include "FairRootManager.h" -#include "Logger.h" - -#include "TBox.h" -#include "TClonesArray.h" -#include "TEllipse.h" -#include "TF1.h" -#include "TFormula.h" -#include "TGraphAsymmErrors.h" -#include "TH1F.h" -#include "TH2F.h" -#include "TMath.h" -#include "TStyle.h" - -#include <algorithm> -#include <fstream> -#include <iomanip> -#include <numeric> - -#include "L1Constants.h" ClassImp(CbmCaInputQaMuch); -namespace phys = L1Constants::phys; // from L1Constants.h - // --------------------------------------------------------------------------------------------------------------------- // -CbmCaInputQaMuch::CbmCaInputQaMuch(int verbose, bool isMCUsed) : CbmQaTask("CbmCaInputQaMuch", verbose, isMCUsed) {} +CbmCaInputQaMuch::CbmCaInputQaMuch(int verbose, bool isMCUsed) : CbmCaInputQaBase("CbmCaInputQaMuch", verbose, isMCUsed) +{ + // Default parameters of task + DefineParameters(); +} // --------------------------------------------------------------------------------------------------------------------- // bool CbmCaInputQaMuch::Check() { - bool res = true; - - int nSt = fpDetInterface->GetNtrackingStations(); - - - // ************************************************************** - // ** Basic checks, available both for real and simulated data ** - // ************************************************************** - - // ----- Checks for mismatches in the ordering of the stations - // - std::vector<double> vStationPos(nSt, 0.); - for (int iSt = 0; iSt < nSt; ++iSt) { - vStationPos[iSt] = fpDetInterface->GetZref(iSt); - } - - if (!std::is_sorted(vStationPos.cbegin(), vStationPos.cend(), [](int l, int r) { return l <= r; })) { - if (fVerbose > 0) { - LOG(error) << fName << ": stations are ordered improperly along the beam axis:"; - for (auto zPos : vStationPos) { - LOG(error) << "\t- " << zPos; - } - } - res = false; - } - - // ----- Checks for mismatch between station and hit z positions - // The purpose of this block is to be ensured, that hits belong to the correct tracking station. For each tracking - // station a unified position z-coordinate is defined, which generally differs from the corresponding positions of - // reconstructed hits. This is due to non-regular position along the beam axis for each detector sensor. Thus, the - // positions of the hits along z-axis are distributed relatively to the position of the station in some range. - // If hits goes out this range, it can signal about a mismatch between hits and geometry. For limits of the range - // one should select large enough values to cover the difference described above and in the same time small enough - // to avoid overlaps with the neighboring stations. - // For each station, a distribution of z_{hit} - z_{st} is integrated over the defined range and scaled by the - // total number of entries to the distribution. If this value is smaller then unity, then some of the hits belong to - // another station. - // - for (int iSt = 0; iSt < nSt; ++iSt) { - int nHits = fvph_hit_station_delta_z[iSt]->GetEntries(); - if (!nHits) { - LOG_IF(error, fVerbose > 0) << fName << ": station " << iSt << " does not have hits"; - res = false; - continue; - } - int iBinMin = fvph_hit_station_delta_z[iSt]->FindBin(-fMaxDiffZStHit); - int iBinMax = fvph_hit_station_delta_z[iSt]->FindBin(+fMaxDiffZStHit); - - if (fvph_hit_station_delta_z[iSt]->Integral(iBinMin, iBinMax) < nHits) { - LOG_IF(error, fVerbose > 0) << fName << ": station " << iSt << " has mismatches in hit z-positions"; - res = false; - } - } - - - // ******************************************************* - // ** Additional checks, if MC information is available ** - // ******************************************************* - - if (IsMCUsed()) { - // ----- Check efficiencies - // Error is raised, if any station has integrated efficiency lower then a selected threshold. - // Integrated efficiency is estimated with fitting the efficiency curve for hit distances (r) with a uniform - // distribution in the range from kEffRangeMin to kEffRangeMax, where the efficiency is assigned to be constant - // - // Fit efficiency curves - LOG(info) << "-- Hit efficiency integrated over hit distance from station center"; - LOG(info) << "\tintegration range: [" << fEffRange[0] << ", " << fEffRange[1] << "] cm"; - - auto* pEffTable = MakeQaObject<CbmQaTable>("eff_table", "Efficiency table", nSt, 2); - pEffTable->SetNamesOfCols({"Station ID", "Efficiency"}); - pEffTable->SetColWidth(20); - - for (int iSt = 0; iSt < nSt; ++iSt) { - auto [eff, effEL, effEU] = fvpe_reco_eff_vs_r[iSt]->GetTotalEfficiency(); - pEffTable->SetCell(iSt, 0, iSt); - pEffTable->SetCell(iSt, 1, eff); - res = CheckRange("Hit finder efficiency in station " + std::to_string(iSt), eff, fEffThrsh, 1.000); - } - - LOG(info) << '\n' << pEffTable->ToString(3); - - // ----- Checks for residuals - // Check hits for potential biases from central values - - // Function to fit a residual distribution, returns a structure - auto FitResidualsAndGetMean = [&](TH1* pHist) { - auto fit = TF1("fitres", "gausn"); - double statMean = pHist->GetMean(); - double statSigm = pHist->GetStdDev(); - fit.SetParameters(100., statMean, statSigm); - pHist->Fit("fitres", "MQ"); - pHist->Fit("fitres", "MQ"); - pHist->Fit("fitres", "MQ"); - // NOTE: Several fit procedures are made to avoid empty fit results - std::array<double, 3> result; - result[0] = fit.GetParameter(1); - result[1] = -fit.GetParameter(2) * fResMeanThrsh; - result[2] = +fit.GetParameter(2) * fResMeanThrsh; - return result; - }; - - auto* pResidualsTable = - MakeQaObject<CbmQaTable>("residuals_mean", "Residual mean values in different stations", nSt, 4); - pResidualsTable->SetNamesOfCols({"Station ID", "Residual(x) [cm]", "Residual(y) [cm]", "Residual(t) [ns]"}); - pResidualsTable->SetColWidth(20); - - // Fill residuals fit results and check boundaries - for (int iSt = 0; iSt < nSt; ++iSt) { - auto fitX = FitResidualsAndGetMean(fvph_res_x[iSt]); - auto fitY = FitResidualsAndGetMean(fvph_res_y[iSt]); - auto fitT = FitResidualsAndGetMean(fvph_res_t[iSt]); - res = CheckRange("Mean of x residual [cm] in station " + std::to_string(iSt), fitX[0], fitX[1], fitX[2]) && res; - res = CheckRange("Mean of y residual [cm] in station " + std::to_string(iSt), fitY[0], fitY[1], fitY[2]) && res; - res = CheckRange("Mean of t residual [ns] in station " + std::to_string(iSt), fitT[0], fitT[1], fitT[2]) && res; - pResidualsTable->SetCell(iSt, 0, iSt); - pResidualsTable->SetCell(iSt, 1, fitX[0]); - pResidualsTable->SetCell(iSt, 2, fitX[1]); - pResidualsTable->SetCell(iSt, 3, fitX[2]); - } - - LOG(info) << '\n' << pResidualsTable->ToString(8); - - // ----- Checks for pulls - // For the hit pull is defined as a ration of the difference between hit and MC-point position or time component - // values and an error of the hit position (or time) component, calculated in the same z-positions. In ideal case, - // when the resolutions are well determined for detecting stations, the pull distributions should have a RMS equal - // to unity. Here we allow a RMS of the pull distributions to be varied in a predefined range. If the RMS runs out - // this range, QA task fails. - // TODO: Add checks for center - - // Fit pull distributions - - // ** Kaniadakis Gaussian distribution, gives smaller chi2 / NDF - //if (!gROOT->FindObject("Expk")) { - // new TFormula("Expk", "TMath::Power(TMath::Sqrt(1 + x[1] * x[1] * x[0] * x[0]) + x[0] * x[1], 1./x[1])"); - //} - //TF1* pPullFit = new TF1("pullFitGausn", "[0] * Expk(-[2] * (x - [1]) * (x - [1]), [3])", -10., 10.); - - auto FitPullsAndGetSigma = [&](TH1* pHist) { - auto fit = TF1("fitpull", "gausn(0)"); - fit.SetParameters(100, 0.001, 1.000); - pHist->Fit("fitpull", "Q"); - return fit.GetParameter(2); - }; - - auto* pPullsTable = MakeQaObject<CbmQaTable>("pulls_rms", "Pulls std. dev. values in different stations", nSt, 4); - pPullsTable->SetNamesOfCols({"Station ID", "Pull(x) sigm", "Pull(y) sigm", "Pull(z) sigm"}); - pPullsTable->SetColWidth(20); - - for (int iSt = 0; iSt < nSt; ++iSt) { - double rmsPullX = FitPullsAndGetSigma(fvph_pull_x[iSt]); - double rmsPullY = FitPullsAndGetSigma(fvph_pull_y[iSt]); - double rmsPullT = FitPullsAndGetSigma(fvph_pull_t[iSt]); - - double rmsPullHi = 1. + fPullWidthThrsh; - double rmsPullLo = 1. - fPullWidthThrsh; - - res = CheckRange("Rms for x pull in station " + std::to_string(iSt), rmsPullX, rmsPullLo, rmsPullHi) && res; - res = CheckRange("Rms for y pull in station " + std::to_string(iSt), rmsPullY, rmsPullLo, rmsPullHi) && res; - res = CheckRange("Rms for t pull in station " + std::to_string(iSt), rmsPullT, rmsPullLo, rmsPullHi) && res; - pPullsTable->SetCell(iSt, 0, iSt); - pPullsTable->SetCell(iSt, 1, rmsPullX); - pPullsTable->SetCell(iSt, 2, rmsPullY); - pPullsTable->SetCell(iSt, 3, rmsPullT); - } - - LOG(info) << '\n' << pPullsTable->ToString(3); - } + bool res = CbmCaInputQaBase::Check(); return res; } @@ -231,281 +33,44 @@ bool CbmCaInputQaMuch::Check() // void CbmCaInputQaMuch::DeInit() { - // Vectors with pointers to histograms - fvph_hit_ypos_vs_xpos.clear(); - fvph_hit_xpos_vs_zpos.clear(); - fvph_hit_ypos_vs_zpos.clear(); - - fvph_hit_station_delta_z.clear(); - - fvph_hit_dx.clear(); - fvph_hit_dy.clear(); - fvph_hit_dt.clear(); - - fvph_n_points_per_hit.clear(); - - fvph_point_ypos_vs_xpos.clear(); - fvph_point_xpos_vs_zpos.clear(); - fvph_point_ypos_vs_zpos.clear(); - - fvph_point_hit_delta_z.clear(); - - fvph_res_x.clear(); - fvph_res_y.clear(); - fvph_res_t.clear(); - - fvph_pull_x.clear(); - fvph_pull_y.clear(); - fvph_pull_t.clear(); - - fvph_res_x_vs_x.clear(); - fvph_res_y_vs_y.clear(); - fvph_res_t_vs_t.clear(); - - fvph_pull_x_vs_x.clear(); - fvph_pull_y_vs_y.clear(); - fvph_pull_t_vs_t.clear(); - - fvpe_reco_eff_vs_xy.clear(); - fvpe_reco_eff_vs_r.clear(); + CbmCaInputQaBase::DeInit(); } // --------------------------------------------------------------------------------------------------------------------- // -void CbmCaInputQaMuch::FillHistograms() +void CbmCaInputQaMuch::DefineParameters() { - int nSt = fpDetInterface->GetNtrackingStations(); - int nHits = fpHits->GetEntriesFast(); - int nMCevents = (IsMCUsed()) ? fpMCEventList->GetNofEvents() : -1; - - std::vector<std::vector<int>> vNofHitsPerPoint; // Number of hits per MC point in different MC events - - if (IsMCUsed()) { - vNofHitsPerPoint.resize(nMCevents); - for (int iE = 0; iE < nMCevents; ++iE) { - int iFile = fpMCEventList->GetFileIdByIndex(iE); - int iEvent = fpMCEventList->GetEventIdByIndex(iE); - int nPoints = fpMCPoints->Size(iFile, iEvent); - vNofHitsPerPoint[iE].resize(nPoints, 0); - } - } - - for (int iH = 0; iH < nHits; ++iH) { - const auto* pHit = dynamic_cast<const CbmMuchPixelHit*>(fpHits->At(iH)); - LOG_IF(fatal, !pHit) << fName << ": hit with iH = " << iH << " is not an CbmMuchPixelHit (dynamic cast failed)"; - - // ************************* - // ** Reconstructed hit QA - - // 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; - - // Hit position - double xHit = pHit->GetX(); - double yHit = pHit->GetY(); - double zHit = pHit->GetZ(); - - double dxHit = pHit->GetDx(); - double dyHit = pHit->GetDy(); - //double dxyHit = pHit->GetDxy(); - - // Hit time - double tHit = pHit->GetTime(); - double dtHit = pHit->GetTimeError(); - - fvph_hit_ypos_vs_xpos[iSt]->Fill(xHit, yHit); - fvph_hit_xpos_vs_zpos[iSt]->Fill(zHit, xHit); - fvph_hit_ypos_vs_zpos[iSt]->Fill(zHit, yHit); - - fvph_hit_station_delta_z[iSt]->Fill(zHit - fpDetInterface->GetZref(iSt)); - - fvph_hit_dx[iSt]->Fill(dxHit); - fvph_hit_dy[iSt]->Fill(dyHit); - fvph_hit_dt[iSt]->Fill(dtHit); - - fvph_hit_ypos_vs_xpos[nSt]->Fill(xHit, yHit); - fvph_hit_xpos_vs_zpos[nSt]->Fill(zHit, xHit); - fvph_hit_ypos_vs_zpos[nSt]->Fill(zHit, yHit); - fvph_hit_dx[nSt]->Fill(dxHit); - fvph_hit_dy[nSt]->Fill(dyHit); - fvph_hit_dt[nSt]->Fill(dtHit); - - - // ********************** - // ** MC information QA - - if (IsMCUsed()) { - CbmMatch* pHitMatch = dynamic_cast<CbmMatch*>(fpHitMatches->At(iH)); - LOG_IF(fatal, !pHitMatch) << fName << ": match object not found for hit ID " << iH; - - // Evaluate number of hits per one MC point - int nMCpoints = 0; // Number of MC points for this hit - for (int iLink = 0; iLink < pHitMatch->GetNofLinks(); ++iLink) { - const auto& link = pHitMatch->GetLink(iLink); - - int iP = link.GetIndex(); // Index of MC point - - // Skip noisy links - if (iP < 0) { continue; } - - ++nMCpoints; - - int iE = fpMCEventList->GetEventIndex(link); - - 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 << ')'; - - LOG_IF(fatal, iP >= (int) vNofHitsPerPoint[iE].size()) - << fName << ": index of MC point is out of range (hit id = " << iH << ", link id = " << iLink - << ", event id = " << iE << ", point id = " << iP << ')'; - vNofHitsPerPoint[iE][iP]++; - } - - fvph_n_points_per_hit[iSt]->Fill(nMCpoints); - fvph_n_points_per_hit[nSt]->Fill(nMCpoints); - - if (nMCpoints != 1) { continue; } - - // The best link to in the match (probably, the cut on nMCpoints is meaningless) - const auto& bestPointLink = pHitMatch->GetMatchedLink(); - - // Skip noisy links - if (bestPointLink.GetIndex() < 0) { continue; } - - // Point matched by the best link - const auto* pMCPoint = dynamic_cast<const CbmMuchPoint*>(fpMCPoints->Get(bestPointLink)); - LOG_IF(fatal, !pMCPoint) << fName << ": MC point object does not exist for hit " << iH; - - // 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: "; - - double t0MC = fpMCEventList->GetEventTime(bestPointLink); - LOG_IF(fatal, t0MC < 0) << fName << ": MC time zero is lower then 0 ns: " << t0MC; - - - // ----- MC point properties - // - double mass = pMCTrack->GetMass(); - //double charge = pMCTrack->GetCharge(); - //double pdg = pMCTrack->GetPdgCode(); - - // Entrance position and time - double xMC = pMCPoint->GetXIn(); - double yMC = pMCPoint->GetYIn(); - double zMC = pMCPoint->GetZIn(); - double tMC = pMCPoint->GetTime() + t0MC; - - // MC point entrance momenta - double pxMC = pMCPoint->GetPx(); - double pyMC = pMCPoint->GetPy(); - double pzMC = pMCPoint->GetPz(); - double pMC = sqrt(pxMC * pxMC + pyMC * pyMC + pzMC * pzMC); - - // MC point exit momenta - // double pxMCo = pMCPoint->GetPxOut(); - // double pyMCo = pMCPoint->GetPyOut(); - // double pzMCo = pMCPoint->GetPzOut(); - // double pMCo = sqrt(pxMCo * pxMCo + pyMCo * pyMCo + pzMCo * pzMCo); - - // Position and time shifted to z-coordinate of the hit - double shiftZ = zHit - zMC; // Difference between hit and point z positions - double xMCs = xMC + shiftZ * pxMC / pzMC; - double yMCs = yMC + shiftZ * pyMC / pzMC; - double tMCs = tMC + shiftZ / (pzMC * phys::kSpeedOfLight) * sqrt(mass * mass + pMC * pMC); - - // Residuals - double xRes = xHit - xMCs; - double yRes = yHit - yMCs; - double tRes = tHit - tMCs; - - // ------ Cuts - - if (std::fabs(pMCPoint->GetPzOut()) < fMinMomentum) { continue; } // CUT ON MINIMUM MOMENTUM - //if (pMCo < cuts::kMinP) { continue; } // Momentum threshold - - - fvph_point_ypos_vs_xpos[iSt]->Fill(xMC, yMC); - fvph_point_xpos_vs_zpos[iSt]->Fill(zMC, xMC); - fvph_point_ypos_vs_zpos[iSt]->Fill(zMC, yMC); - - fvph_point_hit_delta_z[iSt]->Fill(shiftZ); - - fvph_res_x[iSt]->Fill(xRes); - fvph_res_y[iSt]->Fill(yRes); - fvph_res_t[iSt]->Fill(tRes); - - fvph_pull_x[iSt]->Fill(xRes / dxHit); - fvph_pull_y[iSt]->Fill(yRes / dyHit); - fvph_pull_t[iSt]->Fill(tRes / dtHit); - - fvph_res_x_vs_x[iSt]->Fill(xHit, xRes); - fvph_res_y_vs_y[iSt]->Fill(yHit, yRes); - fvph_res_t_vs_t[iSt]->Fill(tHit, tRes); - - fvph_pull_x_vs_x[iSt]->Fill(xHit, xRes / dxHit); - fvph_pull_y_vs_y[iSt]->Fill(yHit, yRes / dyHit); - fvph_pull_t_vs_t[iSt]->Fill(tHit, tRes / dtHit); - - fvph_point_ypos_vs_xpos[nSt]->Fill(xMC, yMC); - fvph_point_xpos_vs_zpos[nSt]->Fill(zMC, xMC); - fvph_point_ypos_vs_zpos[nSt]->Fill(zMC, yMC); - - fvph_point_hit_delta_z[nSt]->Fill(shiftZ); - - fvph_res_x[nSt]->Fill(xRes); - fvph_res_y[nSt]->Fill(yRes); - fvph_res_t[nSt]->Fill(tRes); - - fvph_pull_x[nSt]->Fill(xRes / dxHit); - fvph_pull_y[nSt]->Fill(yRes / dyHit); - fvph_pull_t[nSt]->Fill(tRes / dtHit); - - fvph_res_x_vs_x[nSt]->Fill(xHit, xRes); - fvph_res_y_vs_y[nSt]->Fill(yHit, yRes); - fvph_res_t_vs_t[nSt]->Fill(tHit, tRes); - - fvph_pull_x_vs_x[nSt]->Fill(xHit, xRes / dxHit); - fvph_pull_y_vs_y[nSt]->Fill(yHit, yRes / dyHit); - fvph_pull_t_vs_t[nSt]->Fill(tHit, tRes / dtHit); - } - } // loop over hits: end - - // Fill hit reconstruction efficiencies - if (IsMCUsed()) { - for (int iE = 0; iE < nMCevents; ++iE) { - int iFile = fpMCEventList->GetFileIdByIndex(iE); - int iEvent = fpMCEventList->GetEventIdByIndex(iE); - int nPoints = fpMCPoints->Size(iFile, iEvent); - std::cout << "points:\n"; - for (int iP = 0; iP < nPoints; ++iP) { - const auto* pMCPoint = dynamic_cast<const CbmMuchPoint*>(fpMCPoints->Get(iFile, iEvent, iP)); - LOG_IF(fatal, !pMCPoint) << fName << ": MC point does not exist for iFile = " << iFile - << ", iEvent = " << iEvent << ", iP = " << iP; - int address = pMCPoint->GetDetectorID(); - int iSt = fpDetInterface->GetTrackingStationIndex(address); - LOG_IF(fatal, iSt < 0 || iSt >= nSt) - << fName << ": MC point for FEI = " << iFile << ", " << iEvent << ", " << iP << " and address " << address - << " has wrong station index: iSt = " << iSt; - - double xMC = pMCPoint->GetXIn(); - double yMC = pMCPoint->GetYIn(); - double rMC = sqrt(xMC * xMC + yMC * yMC); - - bool ifPointHasHits = (vNofHitsPerPoint[iE][iP] > 0); - fvpe_reco_eff_vs_xy[iSt]->Fill(ifPointHasHits, xMC, yMC); - fvpe_reco_eff_vs_xy[nSt]->Fill(ifPointHasHits, xMC, yMC); - - fvpe_reco_eff_vs_r[iSt]->Fill(ifPointHasHits, rMC); - fvpe_reco_eff_vs_r[nSt]->Fill(ifPointHasHits, rMC); + auto SetRange = [] (std::array<double, 2>& range, double min, double max) { range[0] = min; range[1] = max; }; + // Hit errors + SetRange(fRHitDx, 0.0000, 5.00); // [cm] + SetRange(fRHitDy, 0.0000, 5.00); // [cm] + SetRange(fRHitDu, 0.0000, 5.00); // [cm] + SetRange(fRHitDv, 0.0000, 5.00); // [cm] + SetRange(fRHitDt, 0.0000, 0.15); // [ns] + // Residuals + SetRange(fRResX, -2.00, 2.00); + SetRange(fRResY, -4.00, 4.00); + SetRange(fRResU, -2.00, 2.00); + SetRange(fRResV, -4.00, 4.00); + SetRange(fRResT, -5.00, 5.00); + // QA result selection criteria + SetRange(fEffRange, 10.0, 30.0); ///< Range for hit efficiency approximation + fResMeanThrsh = 0.50; ///< Maximum allowed deviation of residual mean from zero [sigma] + fPullMeanThrsh = 0.10; ///< Maximum allowed deviation of pull mean from zero + fPullWidthThrsh = 2.00; ///< Maximum allowed deviation of pull width from unity + fEffThrsh = 0.50; ///< Threshold for hit efficiency in the selected range + fMaxDiffZStHit = 1.00; ///< Maximum allowed difference between z-position of hit and station [cm] + fMinMomentum = 0.05; ///< Minimum momentum of particle [GeV/c] + // Track selection criteria + fTrackSelectionPrimary = true; ///< must the track come from the primary vertex + fTrackSelectionMomentum = 0.1; ///< track momentum GeV when it exits the tracking station + fTrackSelectionAngle = 60.; ///< track incl. angle [grad] when it exits the station +} - } // loop over MC-points: end - } // loop over MC-events: end - } // MC is used: end +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmCaInputQaMuch::FillHistogramsPerHit() +{ } // --------------------------------------------------------------------------------------------------------------------- @@ -515,45 +80,9 @@ InitStatus CbmCaInputQaMuch::InitDataBranches() // STS tracking detector interface fpDetInterface = CbmMuchTrackingInterface::Instance(); - LOG_IF(fatal, !fpDetInterface) << "\033[1;31m" << fName << ": tracking detector interface is undefined\033[0m"; - - // FairRootManager - auto* pFairRootManager = FairRootManager::Instance(); - LOG_IF(fatal, !pFairRootManager) << "\033[1;31m" << fName << ": FairRootManager instance is a null pointer\033[0m"; - - // Time-slice - fpTimeSlice = dynamic_cast<CbmTimeSlice*>(pFairRootManager->GetObject("TimeSlice.")); - LOG_IF(fatal, !fpTimeSlice) << "\033[1;31m" << fName << ": time-slice branch is not found\033[0m"; - - // ----- Reconstructed data-branches initialization - - // Hits container - fpHits = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject("MuchPixelHit")); - LOG_IF(fatal, !fpHits) << "\033[1;31m" << fName << ": container of reconstructed hits in MuCh is not found\033[0m"; - - // ----- MC-information branches initialization - if (IsMCUsed()) { - // MC manager - fpMCDataManager = dynamic_cast<CbmMCDataManager*>(pFairRootManager->GetObject("MCDataManager")); - LOG_IF(fatal, !fpMCDataManager) << "\033[1;31m" << fName << ": MC data manager branch is not found\033[0m"; - - // MC event list - fpMCEventList = dynamic_cast<CbmMCEventList*>(pFairRootManager->GetObject("MCEventList.")); - LOG_IF(fatal, !fpMCEventList) << "\033[1;31m" << fName << ": MC event list branch is not found\033[0m"; - - // MC tracks - fpMCTracks = fpMCDataManager->InitBranch("MCTrack"); - LOG_IF(fatal, !fpMCTracks) << "\033[1;31m" << fName << ": MC track branch is not found\033[0m"; - - // MC points - fpMCPoints = fpMCDataManager->InitBranch("MuchPoint"); - LOG_IF(fatal, !fpMCTracks) << "\033[1;31m" << fName << ": MC point branch is not found for MuCh\033[0m"; - - // Hit matches - fpHitMatches = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject("MuchPixelHitMatch")); - LOG_IF(fatal, !fpHitMatches) << "\033[1;31m]" << fName << ": hit match branch is not found for MuCh\033[0m"; - } - + auto baseInitStatus = CbmCaInputQaBase::InitDataBranches(); + if (kSUCCESS != baseInitStatus) { return baseInitStatus; } + return kSUCCESS; } @@ -561,239 +90,9 @@ InitStatus CbmCaInputQaMuch::InitDataBranches() // InitStatus CbmCaInputQaMuch::InitCanvases() { - gStyle->SetOptFit(1); - int nSt = fpDetInterface->GetNtrackingStations(); - - // ******************** - // ** Drawing options - - // Contours - constexpr auto contColor = kOrange + 7; - constexpr auto contWidth = 2; // Line width - constexpr auto contStyle = 2; // Line style - constexpr auto contFill = 0; // Fill style - - // ********************************* - // ** Hit occupancies - - // ** Occupancy in XY plane - auto* pc_hit_ypos_vs_xpos = MakeQaObject<CbmQaCanvas>("hit_ypos_vs_xpos", "", 1600, 800); - pc_hit_ypos_vs_xpos->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_hit_ypos_vs_xpos->cd(iSt + 1); - fvph_hit_ypos_vs_xpos[iSt]->DrawCopy("colz", ""); - - // Square boarders of station, which are used for the field approximation - double stXmin = +fpDetInterface->GetXmax(iSt); - double stXmax = -fpDetInterface->GetXmax(iSt); - double stYmin = -fpDetInterface->GetYmax(iSt); - double stYmax = +fpDetInterface->GetYmax(iSt); - - auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - - // ----- Occupancy projections - auto* pc_hit_xpos = MakeQaObject<CbmQaCanvas>("hit_xpos", "", 1400, 800); - pc_hit_xpos->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_hit_xpos->cd(iSt + 1); - auto* phProj = fvph_hit_ypos_vs_xpos[iSt]->ProjectionY(); - phProj->DrawCopy(); - } - - auto* pc_hit_ypos = MakeQaObject<CbmQaCanvas>("hit_ypos", "", 1400, 800); - pc_hit_ypos->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_hit_ypos->cd(iSt + 1); - auto* phProj = fvph_hit_ypos_vs_xpos[iSt]->ProjectionX(); - phProj->DrawCopy(); - } - - - // ** Occupancy in XZ plane - MakeQaObject<CbmQaCanvas>("hit_xpos_vs_zpos", "", 600, 400); - fvph_hit_xpos_vs_zpos[nSt]->DrawCopy("colz", ""); - for (int iSt = 0; iSt < nSt; ++iSt) { - // Station positions in detector IFS - double stZmin = fpDetInterface->GetZmin(iSt); - double stZmax = fpDetInterface->GetZmax(iSt); - double stHmin = -fpDetInterface->GetXmax(iSt); - double stHmax = +fpDetInterface->GetXmax(iSt); - - auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - - // ** Occupancy in YZ plane - MakeQaObject<CbmQaCanvas>("hit_ypos_vs_zpos", "", 600, 400); - fvph_hit_ypos_vs_zpos[nSt]->DrawCopy("colz", ""); - for (int iSt = 0; iSt < nSt; ++iSt) { - // Station positions in detector IFS - double stZmin = fpDetInterface->GetZmin(iSt); - double stZmax = fpDetInterface->GetZmax(iSt); - double stHmin = -fpDetInterface->GetYmax(iSt); - double stHmax = +fpDetInterface->GetYmax(iSt); - - auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - - // ************ - // ************ - - if (IsMCUsed()) { - - // ********************** - // ** Point occupancies - - // ** Occupancy in XY plane - auto* pc_point_ypos_vs_xpos = MakeQaObject<CbmQaCanvas>("point_ypos_vs_xpos", "", 1600, 800); - pc_point_ypos_vs_xpos->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_point_ypos_vs_xpos->cd(iSt + 1); - fvph_point_ypos_vs_xpos[iSt]->DrawCopy("colz", ""); - - // Square boarders of station, which are used for the field approximation - double stXmin = +fpDetInterface->GetXmax(iSt); - double stXmax = -fpDetInterface->GetXmax(iSt); - double stYmin = -fpDetInterface->GetYmax(iSt); - double stYmax = +fpDetInterface->GetYmax(iSt); - - auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - - // ** Occupancy in XZ plane - MakeQaObject<CbmQaCanvas>("point_xpos_vs_zpos", "", 600, 400); - fvph_point_xpos_vs_zpos[nSt]->SetStats(false); - fvph_point_xpos_vs_zpos[nSt]->DrawCopy("colz", ""); - for (int iSt = 0; iSt < nSt; ++iSt) { - // Station positions in detector IFS - double stZmin = fpDetInterface->GetZmin(iSt); - double stZmax = fpDetInterface->GetZmax(iSt); - double stHmin = -fpDetInterface->GetXmax(iSt); - double stHmax = +fpDetInterface->GetXmax(iSt); - - auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - - // ** Occupancy in YZ plane - MakeQaObject<CbmQaCanvas>("point_ypos_vs_zpos", "", 600, 400); - fvph_point_ypos_vs_zpos[nSt]->SetStats(false); - fvph_point_ypos_vs_zpos[nSt]->DrawCopy("colz", ""); - for (int iSt = 0; iSt < nSt; ++iSt) { - // Station positions in detector IFS - double stZmin = fpDetInterface->GetZmin(iSt); - double stZmax = fpDetInterface->GetZmax(iSt); - double stHmin = -fpDetInterface->GetYmax(iSt); - double stHmax = +fpDetInterface->GetYmax(iSt); - - auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - - // ************** - // ** Residuals - - // x-coordinate - auto* pc_res_x = MakeQaObject<CbmQaCanvas>("res_x", "Residuals for x coordinate", 1600, 800); - pc_res_x->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_res_x->cd(iSt + 1); - fvph_res_x[iSt]->DrawCopy("", ""); - } - - // y-coordinate - auto* pc_res_y = MakeQaObject<CbmQaCanvas>("res_y", "Residuals for y coordinate", 1600, 800); - pc_res_y->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_res_y->cd(iSt + 1); - fvph_res_y[iSt]->DrawCopy("", ""); - } - - // time - auto* pc_res_t = MakeQaObject<CbmQaCanvas>("res_t", "Residuals for time", 1600, 800); - pc_res_t->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_res_t->cd(iSt + 1); - fvph_res_t[iSt]->DrawCopy("", ""); - } - - - // ********** - // ** Pulls - - // x-coordinate - auto* pc_pull_x = MakeQaObject<CbmQaCanvas>("pull_x", "Pulls for x coordinate", 1600, 800); - pc_pull_x->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_pull_x->cd(iSt + 1); - fvph_pull_x[iSt]->DrawCopy("", ""); - } - - // y-coordinate - auto* pc_pull_y = MakeQaObject<CbmQaCanvas>("pull_y", "Pulls for y coordinate", 1600, 800); - pc_pull_y->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_pull_y->cd(iSt + 1); - fvph_pull_y[iSt]->DrawCopy("", ""); - } - - // time - auto* pc_pull_t = MakeQaObject<CbmQaCanvas>("pull_t", "Pulls for time", 1600, 800); - pc_pull_t->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_pull_t->cd(iSt + 1); - fvph_pull_t[iSt]->DrawCopy("", ""); - } - - - // ************************************ - // ** Hit reconstruction efficiencies - - auto* pc_reco_eff_vs_r = - MakeQaObject<CbmQaCanvas>("reco_eff_vs_r", "Hit efficiencies wrt distance from center", 1600, 800); - pc_reco_eff_vs_r->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_reco_eff_vs_r->cd(iSt + 1); - fvpe_reco_eff_vs_r[iSt]->SetMarkerStyle(20); - fvpe_reco_eff_vs_r[iSt]->DrawCopy("AP", ""); - } - - auto* pc_reco_eff_vs_xy = MakeQaObject<CbmQaCanvas>("reco_eff_vs_xy", "Hit efficiencies wrt x and y", 1600, 800); - pc_reco_eff_vs_xy->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_reco_eff_vs_xy->cd(iSt + 1); - fvpe_reco_eff_vs_xy[iSt]->DrawCopy("colz", ""); - } - } - + auto baseInitStatus = CbmCaInputQaBase::InitCanvases(); + if (kSUCCESS != baseInitStatus) { return baseInitStatus; } + return kSUCCESS; } @@ -801,172 +100,9 @@ InitStatus CbmCaInputQaMuch::InitCanvases() // InitStatus CbmCaInputQaMuch::InitHistograms() { - int nSt = fpDetInterface->GetNtrackingStations(); - - // ----- Histograms initialization - fvph_hit_ypos_vs_xpos.resize(nSt + 1, nullptr); - fvph_hit_ypos_vs_zpos.resize(nSt + 1, nullptr); - fvph_hit_xpos_vs_zpos.resize(nSt + 1, nullptr); - - fvph_hit_station_delta_z.resize(nSt); - - fvph_hit_dx.resize(nSt + 1, nullptr); - fvph_hit_dy.resize(nSt + 1, nullptr); - fvph_hit_dt.resize(nSt + 1, nullptr); - - for (int iSt = 0; iSt <= nSt; ++iSt) { - TString nsuff = (iSt == nSt) ? "" : Form("_st%d", iSt); // Histogram name suffix - TString tsuff = (iSt == nSt) ? "" : Form(" in MuCh station %d", iSt); // Histogram title suffix - TString sN = ""; - TString sT = ""; - - // Hit occupancy - sN = (TString) "hit_ypos_vs_xpos" + nsuff; - sT = (TString) "Hit occupancy in xy-Plane" + tsuff + ";x_{hit} [cm];y_{hit} [cm]"; - fvph_hit_ypos_vs_xpos[iSt] = MakeQaObject<TH2F>(sN, sT, kNbins, kRHitX[0], kRHitX[1], kNbins, kRHitY[0], kRHitY[1]); - - sN = (TString) "hit_xpos_vs_zpos" + nsuff; - sT = (TString) "Hit occupancy in xz-Plane" + tsuff + ";z_{hit} [cm];x_{hit} [cm]"; - fvph_hit_xpos_vs_zpos[iSt] = - MakeQaObject<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitX[0], kRHitX[1]); - - sN = (TString) "hit_ypos_vs_zpos" + nsuff; - sT = (TString) "Hit occupancy in yz-plane" + tsuff + ";z_{hit} [cm];y_{hit} [cm]"; - fvph_hit_ypos_vs_zpos[iSt] = - MakeQaObject<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitY[0], kRHitY[1]); - - // Hit errors - sN = (TString) "hit_dx" + nsuff; - sT = (TString) "Hit position error along x-axis" + tsuff + ";dx_{hit} [cm]"; - fvph_hit_dx[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, kRHitDx[0], kRHitDx[1]); - - sN = (TString) "hit_dy" + nsuff; - sT = (TString) "Hit position error along y-axis" + tsuff + ";dy_{hit} [cm]"; - fvph_hit_dy[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, kRHitDy[0], kRHitDy[1]); - - sN = (TString) "hit_dt" + nsuff; - sT = (TString) "Hit time error" + tsuff + ";dt_{hit} [ns]"; - fvph_hit_dt[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, kRHitDt[0], kRHitDt[1]); - - if (iSt == nSt) { continue; } // Following histograms are defined only for separate station - - sN = (TString) "hit_station_delta_z" + nsuff; - sT = (TString) "Different between hit and station z-positions" + tsuff + ";z_{hit} - z_{st} [cm]"; - fvph_hit_station_delta_z[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, -5., 5); + auto baseInitStatus = CbmCaInputQaBase::InitHistograms(); + if (kSUCCESS != baseInitStatus) { return baseInitStatus; } - } // loop over station index: end - - // ----- Initialize histograms, which are use MC-information - if (IsMCUsed()) { - // Resize histogram vectors - fvph_n_points_per_hit.resize(nSt + 1, nullptr); - fvph_point_ypos_vs_xpos.resize(nSt + 1, nullptr); - fvph_point_xpos_vs_zpos.resize(nSt + 1, nullptr); - fvph_point_ypos_vs_zpos.resize(nSt + 1, nullptr); - fvph_point_hit_delta_z.resize(nSt + 1, nullptr); - fvph_res_x.resize(nSt + 1, nullptr); - fvph_res_y.resize(nSt + 1, nullptr); - fvph_res_t.resize(nSt + 1, nullptr); - fvph_pull_x.resize(nSt + 1, nullptr); - fvph_pull_y.resize(nSt + 1, nullptr); - fvph_pull_t.resize(nSt + 1, nullptr); - fvph_res_x_vs_x.resize(nSt + 1, nullptr); - fvph_res_y_vs_y.resize(nSt + 1, nullptr); - fvph_res_t_vs_t.resize(nSt + 1, nullptr); - fvph_pull_x_vs_x.resize(nSt + 1, nullptr); - fvph_pull_y_vs_y.resize(nSt + 1, nullptr); - fvph_pull_t_vs_t.resize(nSt + 1, nullptr); - fvpe_reco_eff_vs_xy.resize(nSt + 1, nullptr); - fvpe_reco_eff_vs_r.resize(nSt + 1, nullptr); - - for (int iSt = 0; iSt <= nSt; ++iSt) { - TString nsuff = (iSt == nSt) ? "" : Form("_st%d", iSt); // Histogram name suffix - TString tsuff = (iSt == nSt) ? "" : Form(" in MuCh station %d", iSt); // Histogram title suffix - TString sN = ""; - TString sT = ""; - - sN = (TString) "n_points_per_hit" + nsuff; - sT = (TString) "Number of points per hit" + tsuff + ";N_{point}/hit"; - fvph_n_points_per_hit[iSt] = MakeQaObject<TH1F>(sN, sT, 10, -0.5, 9.5); - - // Point occupancy - sN = (TString) "point_ypos_vs_xpos" + nsuff; - sT = (TString) "Point occupancy in XY plane" + tsuff + ";x_{MC} [cm];y_{MC} [cm]"; - fvph_point_ypos_vs_xpos[iSt] = - MakeQaObject<TH2F>(sN, sT, kNbins, kRHitX[0], kRHitX[1], kNbins, kRHitY[0], kRHitY[1]); - - sN = (TString) "point_xpos_vs_zpos" + nsuff; - sT = (TString) "Point Occupancy in XZ plane" + tsuff + ";z_{MC} [cm];x_{MC} [cm]"; - fvph_point_xpos_vs_zpos[iSt] = - MakeQaObject<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitX[0], kRHitX[1]); - - sN = (TString) "point_ypos_vs_zpos" + nsuff; - sT = (TString) "Point Occupancy in YZ Plane" + tsuff + ";z_{MC} [cm];y_{MC} [cm]"; - fvph_point_ypos_vs_zpos[iSt] = - MakeQaObject<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitY[0], kRHitY[1]); - - // Difference between MC input z and hit z coordinates - sN = (TString) "point_hit_delta_z" + nsuff; - sT = (TString) "Distance between point and hit along z axis" + tsuff + ";z_{reco} - z_{MC} [cm]"; - fvph_point_hit_delta_z[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, -0.04, 0.04); - - sN = (TString) "res_x" + nsuff; - sT = (TString) "Residuals for x-coordinate" + tsuff + ";x_{reco} - x_{MC} [cm]"; - fvph_res_x[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, kRResX[0], kRResX[1]); - - sN = (TString) "res_y" + nsuff; - sT = (TString) "Residuals for y-coordinate" + tsuff + ";y_{reco} - y_{MC} [cm]"; - fvph_res_y[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, kRResY[0], kRResY[1]); - - sN = (TString) "res_t" + nsuff; - sT = (TString) "Residuals for time" + tsuff + ";t_{reco} - t_{MC} [ns]"; - fvph_res_t[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, kRResT[0], kRResT[1]); - - sN = (TString) "pull_x" + nsuff; - sT = (TString) "Pulls for x-coordinate" + tsuff + ";(x_{reco} - x_{MC}) / #sigma_{x}^{reco}"; - fvph_pull_x[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, kRPullX[0], kRPullX[1]); - - sN = (TString) "pull_y" + nsuff; - sT = (TString) "Pulls for y-coordinate" + tsuff + ";(y_{reco} - y_{MC}) / #sigma_{y}^{reco}"; - fvph_pull_y[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, kRPullY[0], kRPullY[1]); - - sN = (TString) "pull_t" + nsuff; - sT = (TString) "Pulls for time" + tsuff + ";(t_{reco} - t_{MC}) / #sigma_{t}^{reco}"; - fvph_pull_t[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, kRPullT[0], kRPullT[1]); - - sN = (TString) "res_x_vs_x" + nsuff; - sT = (TString) "Residuals for x-coordinate" + tsuff + ";x_{reco} [cm];x_{reco} - x_{MC} [cm]"; - fvph_res_x_vs_x[iSt] = MakeQaObject<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResX[0], kRResX[1]); - - sN = (TString) "res_y_vs_y" + nsuff; - sT = (TString) "Residuals for y-coordinate" + tsuff + ";y_{reco} [cm];y_{reco} - y_{MC} [cm]"; - fvph_res_y_vs_y[iSt] = MakeQaObject<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResY[0], kRResY[1]); - - sN = (TString) "res_t_vs_t" + nsuff; - sT = (TString) "Residuals for time" + tsuff + "t_{reco} [ns];t_{reco} - t_{MC} [ns]"; - fvph_res_t_vs_t[iSt] = MakeQaObject<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResT[0], kRResT[1]); - - sN = (TString) "pull_x_vs_x" + nsuff; - sT = (TString) "Pulls for x-coordinate" + tsuff + "x_{reco} [cm];(x_{reco} - x_{MC}) / #sigma_{x}^{reco}"; - fvph_pull_x_vs_x[iSt] = MakeQaObject<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullX[0], kRPullX[1]); - - sN = (TString) "pull_y_vs_y" + nsuff; - sT = (TString) "Pulls for y-coordinate" + tsuff + ";y_{reco} [cm];(y_{reco} - y_{MC}) / #sigma_{y}^{reco}"; - fvph_pull_y_vs_y[iSt] = MakeQaObject<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullY[0], kRPullY[1]); - - sN = (TString) "pull_t_vs_t" + nsuff; - sT = (TString) "Pulls for time" + tsuff + ";t_{reco} [ns];(t_{reco} - t_{MC}) / #sigma_{t}^{reco}"; - fvph_pull_t_vs_t[iSt] = MakeQaObject<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullT[0], kRPullT[1]); - - - sN = (TString) "reco_eff_vs_xy" + nsuff; - sT = (TString) "Hit rec. efficiency" + tsuff + ";x_{MC} [cm];y_{MC} [cm];#epsilon"; - fvpe_reco_eff_vs_xy[iSt] = MakeQaObject<CbmQaEff>(sN, sT, 100, -50, 50, 100, -50, 50); - - sN = (TString) "reco_eff_vs_r" + nsuff; - sT = (TString) "Hit rec. efficiency" + tsuff + ";r_{MC} [cm];#epsilon"; - fvpe_reco_eff_vs_r[iSt] = MakeQaObject<CbmQaEff>(sN, sT, 100, 0, 100); - } - } return kSUCCESS; } + diff --git a/reco/L1/qa/CbmCaInputQaMuch.h b/reco/L1/qa/CbmCaInputQaMuch.h index 48e65d3853cf7085a89afa26187c706caed4f120..c8f190a712bd7479d41cfeb6396674146b85a370 100644 --- a/reco/L1/qa/CbmCaInputQaMuch.h +++ b/reco/L1/qa/CbmCaInputQaMuch.h @@ -12,155 +12,48 @@ #define CbmCaInputQaMuch_h 1 #include "CbmCaInputQaBase.h" -#include "CbmMCDataManager.h" -#include "CbmQaTask.h" - -#include "TMath.h" - -#include <set> -#include <unordered_map> -#include <vector> - -class CbmMCEventList; -class CbmMCDataArray; -class CbmMCDataManager; -class CbmTimeSlice; -class CbmMatch; -class CbmMuchPixelHit; -class CbmMuchTrackingInterface; -class TClonesArray; -class TH1F; -class CbmQaEff; - -/// A QA-task class, which provides assurance of MuCh hits and geometry -class CbmCaInputQaMuch : public CbmQaTask, public CbmCaInputQaBase { + +/// A QA-task class, which provides assurance of TOF hits and geometry +/// +class CbmCaInputQaMuch : public CbmCaInputQaBase<L1DetectorID::kMuch> { public: - /// Constructor from parameters - /// \param verbose Verbose level - /// \param isMCUsed Flag, whether MC information is available for this task + /// @brief Constructor from parameters + /// @param verbose Verbose level + /// @param isMCUsed Flag, whether MC information is available for this task CbmCaInputQaMuch(int verbose, bool isMCUsed); - ClassDef(CbmCaInputQaMuch, 0); - protected: - // ******************************************** - // ** Virtual method override from CbmQaTask ** - // ******************************************** - - /// Checks results of the QA - /// \return Success flag + /// @brief Checks results of the QA + /// @return Success flag bool Check(); - /// Initializes data branches + /// @brief Initializes data branches InitStatus InitDataBranches(); - /// Initializes canvases + /// @brief Initializes canvases InitStatus InitCanvases(); - /// Initializes histograms + /// @brief Initializes histograms InitStatus InitHistograms(); - /// Fills histograms per event or time-slice - void FillHistograms(); + /// @brief Fills histograms per event or time-slice + void FillHistograms() { CbmCaInputQaBase::FillHistograms(); } - /// De-initializes histograms + /// @brief De-initializes histograms void DeInit(); -private: - // ********************* - // ** Private methods ** - // ********************* - - // ********************* - // ** Class variables ** - // ********************* - - - // ** Data branches ** - CbmMuchTrackingInterface* fpDetInterface = nullptr; ///< Instance of detector interface - - CbmTimeSlice* fpTimeSlice = nullptr; ///< Pointer to current time-slice - - TClonesArray* fpHits = nullptr; ///< Array of hits - - CbmMCDataManager* fpMCDataManager = nullptr; ///< MC data manager - CbmMCEventList* fpMCEventList = nullptr; ///< MC event list - CbmMCDataArray* fpMCTracks = nullptr; ///< Array of MC tracks - CbmMCDataArray* fpMCPoints = nullptr; ///< Array of MC points - - TClonesArray* fpHitMatches = nullptr; ///< Array of hit matches - - // ** Histograms ** - - static constexpr double kEffRangeMin = 10.; ///< Lower limit of hit distance for efficiency integration [cm] - static constexpr double kEffRangeMax = 30.; ///< Upper limit of hit distance for efficiency integration [cm] - - static constexpr double kRHitX[2] = {-100., 100}; ///< Range for hit x coordinate [cm] - static constexpr double kRHitY[2] = {-100., 100}; ///< Range for hit y coordinate [cm] - static constexpr double kRHitZ[2] = {145., 335}; ///< Range for hit z coordinate [cm] - - static constexpr int kNbins = 200; ///< General number of bins - static constexpr int kNbinsZ = 800; ///< Number of bins for z coordinate axis + /// @brief Defines parameters of the task + void DefineParameters(); - static constexpr double kRHitDx[2] = {-.05, .005}; ///< Range for hit x coordinate [cm] - static constexpr double kRHitDy[2] = {-.05, .005}; ///< Range for hit y coordinate [cm] - static constexpr double kRHitDt[2] = {-10., 10.}; ///< Range for hit time [ns] + /// @brief Fills histograms per hit + void FillHistogramsPerHit(); + + /// @brief Fills histograms per MC point + void FillHistogramsPerPoint() {} - static constexpr double kRResX[2] = {-.02, .02}; ///< Range for residual of x coordinate [cm] - static constexpr double kRResY[2] = {-.1, .1}; ///< Range for residual of y coordinate [cm] - static constexpr double kRResT[2] = {-15., 15.}; ///< Range for residual of time [ns] - - static constexpr double kRPullX[2] = {-10., 10.}; ///< Range for pull of x coordinate - static constexpr double kRPullY[2] = {-10., 10.}; ///< Range for pull of y coordinate - static constexpr double kRPullT[2] = {-5., 5.}; ///< Range for pull of time - - - // NOTE: the last element of each vector stands for integral distribution over all stations - - // hit occupancy - std::vector<TH2F*> fvph_hit_ypos_vs_xpos; ///< hit ypos vs xpos in different stations - std::vector<TH2F*> fvph_hit_xpos_vs_zpos; ///< hit xpos vs zpos in different stations - std::vector<TH2F*> fvph_hit_ypos_vs_zpos; ///< hit ypos vs zpos in different stations - - // difference between hit and station z positions - std::vector<TH1F*> fvph_hit_station_delta_z; ///< Difference between station and hit z positions [cm] - - // hit errors - std::vector<TH1F*> fvph_hit_dx; - std::vector<TH1F*> fvph_hit_dy; - std::vector<TH1F*> fvph_hit_dt; - - // MC points occupancy - std::vector<TH1F*> fvph_n_points_per_hit; ///< number of points per one hit in different stations - - - std::vector<TH2F*> fvph_point_ypos_vs_xpos; ///< point ypos vs xpos in different stations - std::vector<TH2F*> fvph_point_xpos_vs_zpos; ///< point xpos vs zpos in different stations - std::vector<TH2F*> fvph_point_ypos_vs_zpos; ///< point ypos vs zpos in different stations - - std::vector<TH1F*> fvph_point_hit_delta_z; ///< difference between z of hit and MC point in different stations - - // Residuals - std::vector<TH1F*> fvph_res_x; ///< residual for x coordinate in different stations - std::vector<TH1F*> fvph_res_y; ///< residual for y coordinate in different stations - std::vector<TH1F*> fvph_res_t; ///< residual for t coordinate in different stations - - std::vector<TH2F*> fvph_res_x_vs_x; ///< residual for x coordinate in different stations - std::vector<TH2F*> fvph_res_y_vs_y; ///< residual for y coordinate in different stations - std::vector<TH2F*> fvph_res_t_vs_t; ///< residual for t coordinate in different stations - - // Pulls - std::vector<TH1F*> fvph_pull_x; ///< pull for x coordinate in different stations - std::vector<TH1F*> fvph_pull_y; ///< pull for y coordinate in different stations - std::vector<TH1F*> fvph_pull_t; ///< pull for t coordinate in different stations - - std::vector<TH2F*> fvph_pull_x_vs_x; ///< pull for x coordinate in different stations - std::vector<TH2F*> fvph_pull_y_vs_y; ///< pull for y coordinate in different stations - std::vector<TH2F*> fvph_pull_t_vs_t; ///< pull for t coordinate in different stations +private: - // Hit efficiencies - std::vector<CbmQaEff*> fvpe_reco_eff_vs_xy; ///< Efficiency of hit reconstruction vs. x and y coordinates of MC - std::vector<CbmQaEff*> fvpe_reco_eff_vs_r; ///< Efficiency of hit reconstruction vs. distance from center + ClassDef(CbmCaInputQaMuch, 0); }; #endif // CbmCaInputQaMuch_h diff --git a/reco/L1/qa/CbmCaInputQaMvd.cxx b/reco/L1/qa/CbmCaInputQaMvd.cxx new file mode 100644 index 0000000000000000000000000000000000000000..92050240f0d128c73a2885c6fc645af0a73ab8d8 --- /dev/null +++ b/reco/L1/qa/CbmCaInputQaMvd.cxx @@ -0,0 +1,107 @@ +/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// @file CbmCaInputQaMvd.cxx +/// @date 10.09.2023 +/// @brief QA-task for CA tracking input from MVD (implementation) +/// @author S.Zharko <s.zharko@gsi.de> + +#include "CbmCaInputQaMvd.h" +#include "CbmMvdTrackingInterface.h" + +ClassImp(CbmCaInputQaMvd); + +// --------------------------------------------------------------------------------------------------------------------- +// +CbmCaInputQaMvd::CbmCaInputQaMvd(int verbose, bool isMCUsed) : CbmCaInputQaBase("CbmCaInputQaMvd", verbose, isMCUsed) +{ + // Default parameters of task + DefineParameters(); +} + +// --------------------------------------------------------------------------------------------------------------------- +// +bool CbmCaInputQaMvd::Check() +{ + bool res = CbmCaInputQaBase::Check(); + + return res; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmCaInputQaMvd::DeInit() +{ + CbmCaInputQaBase::DeInit(); +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmCaInputQaMvd::DefineParameters() +{ + auto SetRange = [] (std::array<double, 2>& range, double min, double max) { range[0] = min; range[1] = max; }; + // Hit errors + SetRange(fRHitDx, 0.0000, 5.00); // [cm] + SetRange(fRHitDy, 0.0000, 5.00); // [cm] + SetRange(fRHitDu, 0.0000, 5.00); // [cm] + SetRange(fRHitDv, 0.0000, 5.00); // [cm] + SetRange(fRHitDt, 0.0000, 0.15); // [ns] + // Residuals + SetRange(fRResX, -2.00, 2.00); + SetRange(fRResY, -4.00, 4.00); + SetRange(fRResU, -2.00, 2.00); + SetRange(fRResV, -4.00, 4.00); + SetRange(fRResT, -0.50, 0.50); + // QA result selection criteria + SetRange(fEffRange, 10.0, 30.0); ///< Range for hit efficiency approximation + fResMeanThrsh = 0.50; ///< Maximum allowed deviation of residual mean from zero [sigma] + fPullMeanThrsh = 0.10; ///< Maximum allowed deviation of pull mean from zero + fPullWidthThrsh = 2.00; ///< Maximum allowed deviation of pull width from unity + fEffThrsh = 0.50; ///< Threshold for hit efficiency in the selected range + fMaxDiffZStHit = 1.00; ///< Maximum allowed difference between z-position of hit and station [cm] + fMinMomentum = 0.05; ///< Minimum momentum of particle [GeV/c] + // Track selection criteria + fTrackSelectionPrimary = true; ///< must the track come from the primary vertex + fTrackSelectionMomentum = 0.1; ///< track momentum GeV when it exits the tracking station + fTrackSelectionAngle = 60.; ///< track incl. angle [grad] when it exits the station +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmCaInputQaMvd::FillHistogramsPerHit() +{ +} + +// --------------------------------------------------------------------------------------------------------------------- +// +InitStatus CbmCaInputQaMvd::InitDataBranches() +{ + // STS tracking detector interface + fpDetInterface = CbmMvdTrackingInterface::Instance(); + + auto baseInitStatus = CbmCaInputQaBase::InitDataBranches(); + if (kSUCCESS != baseInitStatus) { return baseInitStatus; } + + return kSUCCESS; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +InitStatus CbmCaInputQaMvd::InitCanvases() +{ + auto baseInitStatus = CbmCaInputQaBase::InitCanvases(); + if (kSUCCESS != baseInitStatus) { return baseInitStatus; } + + return kSUCCESS; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +InitStatus CbmCaInputQaMvd::InitHistograms() +{ + auto baseInitStatus = CbmCaInputQaBase::InitHistograms(); + if (kSUCCESS != baseInitStatus) { return baseInitStatus; } + + return kSUCCESS; +} diff --git a/reco/L1/qa/CbmCaInputQaMvd.h b/reco/L1/qa/CbmCaInputQaMvd.h new file mode 100644 index 0000000000000000000000000000000000000000..7715156b76b8a73d18c9b4595af62a6a78c632f6 --- /dev/null +++ b/reco/L1/qa/CbmCaInputQaMvd.h @@ -0,0 +1,59 @@ +/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// \file CbmCaInputQaMvd.h +/// \date 10.09.2023 +/// \brief QA-task for CA tracking input from MVD (header) +/// \author S.Zharko <s.zharko@gsi.de> + + +#ifndef CbmCaInputQaMvd_h +#define CbmCaInputQaMvd_h 1 + +#include "CbmCaInputQaBase.h" + +/// A QA-task class, which provides assurance of TOF hits and geometry +/// +class CbmCaInputQaMvd : public CbmCaInputQaBase<L1DetectorID::kTrd> { +public: + /// @brief Constructor from parameters + /// @param verbose Verbose level + /// @param isMCUsed Flag, whether MC information is available for this task + CbmCaInputQaMvd(int verbose, bool isMCUsed); + +protected: + /// @brief Checks results of the QA + /// @return Success flag + bool Check(); + + /// @brief Initializes data branches + InitStatus InitDataBranches(); + + /// @brief Initializes canvases + InitStatus InitCanvases(); + + /// @brief Initializes histograms + InitStatus InitHistograms(); + + /// @brief Fills histograms per event or time-slice + void FillHistograms() { CbmCaInputQaBase::FillHistograms(); } + + /// @brief De-initializes histograms + void DeInit(); + + /// @brief Defines parameters of the task + void DefineParameters(); + + /// @brief Fills histograms per hit + void FillHistogramsPerHit(); + + /// @brief Fills histograms per MC point + void FillHistogramsPerPoint() {} + +private: + + ClassDef(CbmCaInputQaMvd, 0); +}; + +#endif // CbmCaInputQaMvd_h diff --git a/reco/L1/qa/CbmCaInputQaSts.cxx b/reco/L1/qa/CbmCaInputQaSts.cxx index d98c2148b7473d67d30be7ce34859cf2afc7ac13..e512ec3d6dda119d01c230d0b15b407d2723e4f1 100644 --- a/reco/L1/qa/CbmCaInputQaSts.cxx +++ b/reco/L1/qa/CbmCaInputQaSts.cxx @@ -53,13 +53,10 @@ namespace phys = L1Constants::phys; // from L1Constants.h // --------------------------------------------------------------------------------------------------------------------- // -CbmCaInputQaSts::CbmCaInputQaSts(int verbose, bool isMCUsed) -: CbmCaInputQaBase("CbmCaInputQaSts", verbose, isMCUsed) +CbmCaInputQaSts::CbmCaInputQaSts(int verbose, bool isMCUsed): CbmCaInputQaBase("CbmCaInputQaSts", verbose, isMCUsed) { - // TODO: remove an enlarged threshold when the pulls are fixed - SetPullMeanDeviationThrsh(0.1); - SetPullWidthDeviationThrsh(2.); - SetEfficiencyThrsh(0.5, 0, 100); + // Default parameters of task + DefineParameters(); } // --------------------------------------------------------------------------------------------------------------------- @@ -112,8 +109,8 @@ void CbmCaInputQaSts::DefineParameters() // QA result selection criteria SetRange(fEffRange, 0.0, 100.0); ///< Range for hit efficiency approximation fResMeanThrsh = 0.50; ///< Maximum allowed deviation of residual mean from zero [sigma] - fPullMeanThrsh = 0.05; ///< Maximum allowed deviation of pull mean from zero - fPullWidthThrsh = 0.10; ///< Maximum allowed deviation of pull width from unity + fPullMeanThrsh = 0.10; ///< Maximum allowed deviation of pull mean from zero + fPullWidthThrsh = 2.00; ///< Maximum allowed deviation of pull width from unity fEffThrsh = 0.50; ///< Threshold for hit efficiency in the selected range fMaxDiffZStHit = 1.00; ///< Maximum allowed difference between z-position of hit and station [cm] fMinMomentum = 0.05; ///< Minimum momentum of particle [GeV/c] diff --git a/reco/L1/qa/CbmCaInputQaSts.h b/reco/L1/qa/CbmCaInputQaSts.h index 6fef1d4d6c26b2b1afbd229629137554b2d7ff5e..404c86c1f5839cb5f444c9b6ca90c9599b1643c6 100644 --- a/reco/L1/qa/CbmCaInputQaSts.h +++ b/reco/L1/qa/CbmCaInputQaSts.h @@ -45,10 +45,6 @@ public: CbmCaInputQaSts(int verbose, bool isMCUsed); protected: - // ******************************************** - // ** Virtual method override from CbmQaTask ** - // ******************************************** - /// Checks results of the QA and returns some flag bool Check(); @@ -70,21 +66,12 @@ protected: /// @brief Defines parameters of the task void DefineParameters(); - - /// Checks ranges for mean and standard deviation - /// \return False, if variable exits the range - bool CheckRangePull(TH1* h) - { - return CbmQaTask::CheckRange(h, fPullMeanThrsh, 1. - fPullWidthThrsh, 1. + fPullWidthThrsh); - } - /// @brief Fills histograms per hit void FillHistogramsPerHit(); /// @brief Fills histograms per MC point void FillHistogramsPerPoint(); - private: // ----- Data branches TClonesArray* fpClusters = nullptr; ///< Array of hit clusters diff --git a/reco/L1/qa/CbmCaInputQaTof.cxx b/reco/L1/qa/CbmCaInputQaTof.cxx index 4feb827090073de8ea62d87657ba92c3352d98fa..81704750adddbd6284fe235ec0da9941d0cdca46 100644 --- a/reco/L1/qa/CbmCaInputQaTof.cxx +++ b/reco/L1/qa/CbmCaInputQaTof.cxx @@ -10,28 +10,18 @@ #include "CbmCaInputQaTof.h" #include "CbmAddress.h" -#include "CbmLink.h" -#include "CbmMCDataArray.h" -#include "CbmMCEventList.h" -#include "CbmMCTrack.h" -#include "CbmMatch.h" -#include "CbmQaCanvas.h" -#include "CbmQaEff.h" -#include "CbmTimeSlice.h" #include "CbmTofHit.h" -#include "CbmTofInteraction.h" #include "CbmTofPoint.h" #include "CbmTofTrackingInterface.h" +#include "FairRunAna.h" +#include "FairRuntimeDb.h" #include "FairRootManager.h" #include "Logger.h" #include "TBox.h" #include "TClonesArray.h" -#include "TEllipse.h" #include "TF1.h" -#include "TFormula.h" -#include "TGraphAsymmErrors.h" #include "TH1F.h" #include "TH2F.h" #include "TMath.h" @@ -41,7 +31,6 @@ #include <fstream> #include <iomanip> #include <numeric> -#include <set> // TMP!!!! #include "CaToolsLinkKey.h" #include "L1Constants.h" @@ -52,181 +41,17 @@ namespace phys = L1Constants::phys; // from L1Constants.h // --------------------------------------------------------------------------------------------------------------------- // -CbmCaInputQaTof::CbmCaInputQaTof(int verbose, bool isMCUsed) : CbmQaTask("CbmCaInputQaTof", verbose, isMCUsed) {} +CbmCaInputQaTof::CbmCaInputQaTof(int verbose, bool isMCUsed) : CbmCaInputQaBase("CbmCaInputQaTof", verbose, isMCUsed) +{ + // Default parameters of task + DefineParameters(); +} // --------------------------------------------------------------------------------------------------------------------- // bool CbmCaInputQaTof::Check() { - bool res = true; - - int nSt = fpDetInterface->GetNtrackingStations(); - - - // ************************************************************** - // ** Basic checks, available both for real and simulated data ** - // ************************************************************** - - // ----- Checks for mismatches in the ordering of the stations - // - std::vector<double> vStationPos(nSt, 0.); - for (int iSt = 0; iSt < nSt; ++iSt) { - vStationPos[iSt] = fpDetInterface->GetZref(iSt); - } - - if (!std::is_sorted(vStationPos.cbegin(), vStationPos.cend(), [](int l, int r) { return l <= r; })) { - if (fVerbose > 0) { - LOG(error) << fName << ": stations are ordered improperly along the beam axis:"; - for (auto zPos : vStationPos) { - LOG(error) << "\t- " << zPos; - } - } - res = false; - } - - // ----- Checks for mismatch between station and hit z positions - // The purpose of this block is to be ensured, that hits belong to the correct tracking station. For each tracking - // station a unified position z-coordinate is defined, which generally differs from the corresponding positions of - // reconstructed hits. This is due to non-regular position along the beam axis for each detector sensor. Thus, the - // positions of the hits along z-axis are distributed relatively to the position of the station in some range. - // If hits goes out this range, it can signal about a mismatch between hits and geometry. For limits of the range - // one should select large enough values to cover the difference described above and in the same time small enough - // to avoid overlaps with the neighboring stations. - // For each station, a distribution of z_{hit} - z_{st} is integrated over the defined range and scaled by the - // total number of entries to the distribution. If this value is smaller then unity, then some of the hits belong to - // another station. - // - for (int iSt = 0; iSt < nSt; ++iSt) { - int nHits = fvph_hit_station_delta_z[iSt]->GetEntries(); - if (!nHits) { - LOG_IF(error, fVerbose > 0) << fName << ": station " << iSt << " does not have hits"; - res = false; - continue; - } - int iBinMin = fvph_hit_station_delta_z[iSt]->FindBin(-fMaxDiffZStHit); - int iBinMax = fvph_hit_station_delta_z[iSt]->FindBin(+fMaxDiffZStHit); - - if (fvph_hit_station_delta_z[iSt]->Integral(iBinMin, iBinMax) < nHits) { - LOG_IF(error, fVerbose > 0) << fName << ": station " << iSt << " has mismatches in hit z-positions"; - res = false; - } - } - - - // ******************************************************* - // ** Additional checks, if MC information is available ** - // ******************************************************* - - if (IsMCUsed()) { - // ----- Check efficiencies - // Error is raised, if any station has integrated efficiency lower then a selected threshold. - // Integrated efficiency is estimated with fitting the efficiency curve for hit distances (r) with a uniform - // distribution in the range from kEffRangeMin to kEffRangeMax, where the efficiency is assigned to be constant - // - // Fit efficiency curves - LOG(info) << "-- Hit efficiency integrated over hit distance from station center"; - LOG(info) << "\tintegration range: [" << fEffRange[0] << ", " << fEffRange[1] << "] cm"; - - auto* pEffTable = MakeQaObject<CbmQaTable>("eff_table", "Efficiency table", nSt, 2); - pEffTable->SetNamesOfCols({"Station ID", "Efficiency"}); - pEffTable->SetColWidth(20); - - for (int iSt = 0; iSt < nSt; ++iSt) { - auto [eff, effEL, effEU] = fvpe_reco_eff_vs_r[iSt]->GetTotalEfficiency(); - pEffTable->SetCell(iSt, 0, iSt); - pEffTable->SetCell(iSt, 1, eff); - res = CheckRange("Hit finder efficiency in station " + std::to_string(iSt), eff, fEffThrsh, 1.000); - } - - LOG(info) << '\n' << pEffTable->ToString(3); - - // ----- Checks for residuals - // Check hits for potential biases from central values - - // Function to fit a residual distribution, returns a structure - auto FitResidualsAndGetMean = [&](TH1* pHist) { - auto fit = TF1("fitres", "gausn"); - double statMean = pHist->GetMean(); - double statSigm = pHist->GetStdDev(); - fit.SetParameters(100., statMean, statSigm); - pHist->Fit("fitres", "MQ"); - pHist->Fit("fitres", "MQ"); - pHist->Fit("fitres", "MQ"); - // NOTE: Several fit procedures are made to avoid empty fit results - std::array<double, 3> result; - result[0] = fit.GetParameter(1); - result[1] = -fit.GetParameter(2) * fResMeanThrsh; - result[2] = +fit.GetParameter(2) * fResMeanThrsh; - return result; - }; - - auto* pResidualsTable = - MakeQaObject<CbmQaTable>("residuals_mean", "Residual mean values in different stations", nSt, 4); - pResidualsTable->SetNamesOfCols({"Station ID", "Residual(x) [cm]", "Residual(y) [cm]", "Residual(t) [ns]"}); - pResidualsTable->SetColWidth(20); - - // Fill residuals fit results and check boundaries - for (int iSt = 0; iSt < nSt; ++iSt) { - auto fitX = FitResidualsAndGetMean(fvph_res_x[iSt]); - auto fitY = FitResidualsAndGetMean(fvph_res_y[iSt]); - auto fitT = FitResidualsAndGetMean(fvph_res_t[iSt]); - res = CheckRange("Mean of x residual [cm] in station " + std::to_string(iSt), fitX[0], fitX[1], fitX[2]) && res; - res = CheckRange("Mean of y residual [cm] in station " + std::to_string(iSt), fitY[0], fitY[1], fitY[2]) && res; - res = CheckRange("Mean of t residual [ns] in station " + std::to_string(iSt), fitT[0], fitT[1], fitT[2]) && res; - pResidualsTable->SetCell(iSt, 0, iSt); - pResidualsTable->SetCell(iSt, 1, fitX[0]); - pResidualsTable->SetCell(iSt, 2, fitX[1]); - pResidualsTable->SetCell(iSt, 3, fitX[2]); - } - - LOG(info) << '\n' << pResidualsTable->ToString(8); - - // ----- Checks for pulls - // For the hit pull is defined as a ration of the difference between hit and MC-point position or time component - // values and an error of the hit position (or time) component, calculated in the same z-positions. In ideal case, - // when the resolutions are well determined for detecting stations, the pull distributions should have a RMS equal - // to unity. Here we allow a RMS of the pull distributions to be varied in a predefined range. If the RMS runs out - // this range, QA task fails. - // TODO: Add checks for center - - // Fit pull distributions - - // ** Kaniadakis Gaussian distribution, gives smaller chi2 / NDF - //if (!gROOT->FindObject("Expk")) { - // new TFormula("Expk", "TMath::Power(TMath::Sqrt(1 + x[1] * x[1] * x[0] * x[0]) + x[0] * x[1], 1./x[1])"); - //} - //TF1* pPullFit = new TF1("pullFitGausn", "[0] * Expk(-[2] * (x - [1]) * (x - [1]), [3])", -10., 10.); - - auto FitPullsAndGetSigma = [&](TH1* pHist) { - auto fit = TF1("fitpull", "gausn(0)"); - fit.SetParameters(100, 0.001, 1.000); - pHist->Fit("fitpull", "Q"); - return fit.GetParameter(2); - }; - - auto* pPullsTable = MakeQaObject<CbmQaTable>("pulls_rms", "Pulls std. dev. values in different stations", nSt, 4); - pPullsTable->SetNamesOfCols({"Station ID", "Pull(x) sigm", "Pull(y) sigm", "Pull(z) sigm"}); - pPullsTable->SetColWidth(20); - - for (int iSt = 0; iSt < nSt; ++iSt) { - double rmsPullX = FitPullsAndGetSigma(fvph_pull_x[iSt]); - double rmsPullY = FitPullsAndGetSigma(fvph_pull_y[iSt]); - double rmsPullT = FitPullsAndGetSigma(fvph_pull_t[iSt]); - - double rmsPullHi = 1. + fPullWidthThrsh; - double rmsPullLo = 1. - fPullWidthThrsh; - - res = CheckRange("Rms for x pull in station " + std::to_string(iSt), rmsPullX, rmsPullLo, rmsPullHi) && res; - res = CheckRange("Rms for y pull in station " + std::to_string(iSt), rmsPullY, rmsPullLo, rmsPullHi) && res; - res = CheckRange("Rms for t pull in station " + std::to_string(iSt), rmsPullT, rmsPullLo, rmsPullHi) && res; - pPullsTable->SetCell(iSt, 0, iSt); - pPullsTable->SetCell(iSt, 1, rmsPullX); - pPullsTable->SetCell(iSt, 2, rmsPullY); - pPullsTable->SetCell(iSt, 3, rmsPullT); - } - - LOG(info) << '\n' << pPullsTable->ToString(3); - } + bool res = CbmCaInputQaBase::Check(); return res; } @@ -235,370 +60,69 @@ bool CbmCaInputQaTof::Check() // void CbmCaInputQaTof::DeInit() { - // Vectors with pointers to histograms - fvph_hit_ypos_vs_xpos.clear(); - fvph_hit_xpos_vs_zpos.clear(); - fvph_hit_ypos_vs_zpos.clear(); - - fvph_hit_station_delta_z.clear(); - - fvph_hit_dx.clear(); - fvph_hit_dy.clear(); - fvph_hit_dt.clear(); - - fvph_n_points_per_hit.clear(); - - fvph_point_ypos_vs_xpos.clear(); - fvph_point_xpos_vs_zpos.clear(); - fvph_point_ypos_vs_zpos.clear(); - - fvph_point_hit_delta_z.clear(); - - fvph_res_x.clear(); - fvph_res_y.clear(); - fvph_res_t.clear(); - - fvph_pull_x.clear(); - fvph_pull_y.clear(); - fvph_pull_t.clear(); - - fvph_res_x_vs_x.clear(); - fvph_res_y_vs_y.clear(); - fvph_res_t_vs_t.clear(); - - fvph_pull_x_vs_x.clear(); - fvph_pull_y_vs_y.clear(); - fvph_pull_t_vs_t.clear(); - - fvpe_reco_eff_vs_xy.clear(); - fvpe_reco_eff_vs_r.clear(); + CbmCaInputQaBase::DeInit(); + fvph_hit_xy_vs_cell.clear(); + fvph_hit_zx_vs_cell.clear(); + fvph_hit_zy_vs_cell.clear(); } // --------------------------------------------------------------------------------------------------------------------- // -void CbmCaInputQaTof::FillHistograms() +void CbmCaInputQaTof::DefineParameters() { - int nSt = fpDetInterface->GetNtrackingStations(); - int nHits = fpHits->GetEntriesFast(); - - // ----- Fill TOF interactions - std::shared_ptr<TClonesArray> pMCInteractions = nullptr; // Array of MC interactions - std::shared_ptr<TClonesArray> pHitInterMatches = nullptr; // Array of hit matches to MC interactions - if (IsMCUsed()) { - pMCInteractions = std::make_shared<TClonesArray>("CbmTofInteraction"); - pHitInterMatches = std::make_shared<TClonesArray>("CbmMatch", nHits); - FillInteractions(pMCInteractions, pHitInterMatches); - } // IsMCUsed - - std::vector<int> vNofHitsPerInteraction; // Number of hits per MC point in TS - - if (IsMCUsed()) { vNofHitsPerInteraction.resize(pMCInteractions->GetEntriesFast()); } - - for (int iH = 0; iH < nHits; ++iH) { - const auto* pHit = dynamic_cast<const CbmTofHit*>(fpHits->At(iH)); - LOG_IF(fatal, !pHit) << fName << ": hit with iH = " << iH << " is not an CbmTofHit (dynamic cast failed)"; - - // FIXME: Is this cut still needed?? - int32_t address = pHit->GetAddress(); - if (0x00202806 == address || 0x00002806 == address) { continue; } - - if (5 == CbmTofAddress::GetSmType(pHit->GetAddress())) { continue; } + auto SetRange = [] (std::array<double, 2>& range, double min, double max) { range[0] = min; range[1] = max; }; + // Hit errors + SetRange(fRHitDx, 0.0000, 5.00); // [cm] + SetRange(fRHitDy, 0.0000, 5.00); // [cm] + SetRange(fRHitDu, 0.0000, 5.00); // [cm] + SetRange(fRHitDv, 0.0000, 5.00); // [cm] + SetRange(fRHitDt, 0.0000, 0.15); // [ns] + // Residuals + SetRange(fRResX, -2.00, 2.00); + SetRange(fRResY, -4.00, 4.00); + SetRange(fRResU, -2.00, 2.00); + SetRange(fRResV, -4.00, 4.00); + SetRange(fRResT, -0.50, 0.50); + // QA result selection criteria + SetRange(fEffRange, 10.0, 30.0); ///< Range for hit efficiency approximation + fResMeanThrsh = 0.50; ///< Maximum allowed deviation of residual mean from zero [sigma] + fPullMeanThrsh = 0.10; ///< Maximum allowed deviation of pull mean from zero + fPullWidthThrsh = 2.00; ///< Maximum allowed deviation of pull width from unity + fEffThrsh = 0.50; ///< Threshold for hit efficiency in the selected range + fMaxDiffZStHit = 1.00; ///< Maximum allowed difference between z-position of hit and station [cm] + fMinMomentum = 0.05; ///< Minimum momentum of particle [GeV/c] + // Track selection criteria + fTrackSelectionPrimary = true; ///< must the track come from the primary vertex + fTrackSelectionMomentum = 0.1; ///< track momentum GeV when it exits the tracking station + fTrackSelectionAngle = 60.; ///< track incl. angle [grad] when it exits the station - // ************************* - // ** Reconstructed hit QA - - // 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; - - // Hit position - double xHit = pHit->GetX(); - double yHit = pHit->GetY(); - double zHit = pHit->GetZ(); - - double dxHit = pHit->GetDx(); - double dyHit = pHit->GetDy(); - //double dxyHit = pHit->GetDxy(); - - // Hit time - double tHit = pHit->GetTime(); - double dtHit = pHit->GetTimeError(); - - fvph_hit_ypos_vs_xpos[iSt]->Fill(xHit, yHit); - fvph_hit_xpos_vs_zpos[iSt]->Fill(zHit, xHit); - fvph_hit_ypos_vs_zpos[iSt]->Fill(zHit, yHit); - - fvph_hit_station_delta_z[iSt]->Fill(zHit - fpDetInterface->GetZref(iSt)); - - fvph_hit_dx[iSt]->Fill(dxHit); - fvph_hit_dy[iSt]->Fill(dyHit); - fvph_hit_dt[iSt]->Fill(dtHit); - - fvph_hit_ypos_vs_xpos[nSt]->Fill(xHit, yHit); - fvph_hit_xpos_vs_zpos[nSt]->Fill(zHit, xHit); - fvph_hit_ypos_vs_zpos[nSt]->Fill(zHit, yHit); - fvph_hit_dx[nSt]->Fill(dxHit); - fvph_hit_dy[nSt]->Fill(dyHit); - fvph_hit_dt[nSt]->Fill(dtHit); - - - // ********************** - // ** MC information QA - - if (IsMCUsed()) { - CbmMatch* pHitMatch = dynamic_cast<CbmMatch*>(pHitInterMatches->At(iH)); - LOG_IF(fatal, !pHitMatch) << fName << ": match object not found for hit ID " << iH; - - // NOTE: Here we treat interactions simply as MC points - // Evaluate number of hits per one MC point - int nMCpoints = 0; // Number of MC points for this hit - int nLinks = pHitMatch->GetNofLinks(); - for (int iLink = 0; iLink < nLinks; ++iLink) { - const auto& link = pHitMatch->GetLink(iLink); - - int iP = link.GetIndex(); // Index of MC interaction - - // Skip noisy links - if (iP < 0) { continue; } - - ++nMCpoints; - - LOG_IF(fatal, iP >= (int) vNofHitsPerInteraction.size()) - << fName << ": index of MC point is out of range (hit id = " << iH << ", link id = " << iLink - << ", point id = " << iP << ')'; - vNofHitsPerInteraction[iP]++; - } - - fvph_n_points_per_hit[iSt]->Fill(nMCpoints); - fvph_n_points_per_hit[nSt]->Fill(nMCpoints); - - if (nMCpoints != 1) { continue; } // ?? - - // The best link to in the match (probably, the cut on nMCpoints is meaningless) - const auto& bestPointLink = pHitMatch->GetMatchedLink(); - - // Skip noisy links - if (bestPointLink.GetIndex() < 0) { continue; } - - // Point matched by the best link - const auto* pMCPoint = dynamic_cast<const CbmTofInteraction*>(pMCInteractions->At(bestPointLink.GetIndex())); - LOG_IF(fatal, !pMCPoint) << fName << ": MC point object does not exist for hit " << iH; - - // 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: "; - - double t0MC = fpMCEventList->GetEventTime(bestPointLink); - LOG_IF(fatal, t0MC < 0) << fName << ": MC time zero is lower then 0 ns: " << t0MC; - - // ----- MC point properties - // - double mass = pMCTrack->GetMass(); - //double charge = pMCTrack->GetCharge(); - //double pdg = pMCTrack->GetPdgCode(); - - // Entrance position and time - double xMC = pMCPoint->GetX(); - double yMC = pMCPoint->GetY(); - double zMC = pMCPoint->GetZ(); - double tMC = pMCPoint->GetTime() + t0MC; - - // MC point entrance momenta - double pxMC = pMCPoint->GetPx(); - double pyMC = pMCPoint->GetPy(); - double pzMC = pMCPoint->GetPz(); - double pMC = sqrt(pxMC * pxMC + pyMC * pyMC + pzMC * pzMC); - - // Position and time shifted to z-coordinate of the hit - double shiftZ = zHit - zMC; // Difference between hit and point z positions - double xMCs = xMC + shiftZ * pxMC / pzMC; - double yMCs = yMC + shiftZ * pyMC / pzMC; - double tMCs = tMC + shiftZ / (pzMC * phys::kSpeedOfLight) * sqrt(mass * mass + pMC * pMC); - - // Residuals - double xRes = xHit - xMCs; - double yRes = yHit - yMCs; - double tRes = tHit - tMCs; - - // ------ Cuts - - if (std::fabs(pMCPoint->GetPz()) < fMinMomentum) { continue; } // CUT ON MINIMUM MOMENTUM - //if (pMCo < cuts::kMinP) { continue; } // Momentum threshold - - - fvph_point_ypos_vs_xpos[iSt]->Fill(xMC, yMC); - fvph_point_xpos_vs_zpos[iSt]->Fill(zMC, xMC); - fvph_point_ypos_vs_zpos[iSt]->Fill(zMC, yMC); - - fvph_point_hit_delta_z[iSt]->Fill(shiftZ); - - fvph_res_x[iSt]->Fill(xRes); - fvph_res_y[iSt]->Fill(yRes); - fvph_res_t[iSt]->Fill(tRes); - - fvph_pull_x[iSt]->Fill(xRes / dxHit); - fvph_pull_y[iSt]->Fill(yRes / dyHit); - fvph_pull_t[iSt]->Fill(tRes / dtHit); - - fvph_res_x_vs_x[iSt]->Fill(xHit, xRes); - fvph_res_y_vs_y[iSt]->Fill(yHit, yRes); - fvph_res_t_vs_t[iSt]->Fill(tHit, tRes); - - fvph_pull_x_vs_x[iSt]->Fill(xHit, xRes / dxHit); - fvph_pull_y_vs_y[iSt]->Fill(yHit, yRes / dyHit); - fvph_pull_t_vs_t[iSt]->Fill(tHit, tRes / dtHit); - - fvph_point_ypos_vs_xpos[nSt]->Fill(xMC, yMC); - fvph_point_xpos_vs_zpos[nSt]->Fill(zMC, xMC); - fvph_point_ypos_vs_zpos[nSt]->Fill(zMC, yMC); - - fvph_point_hit_delta_z[nSt]->Fill(shiftZ); - - fvph_res_x[nSt]->Fill(xRes); - fvph_res_y[nSt]->Fill(yRes); - fvph_res_t[nSt]->Fill(tRes); - - fvph_pull_x[nSt]->Fill(xRes / dxHit); - fvph_pull_y[nSt]->Fill(yRes / dyHit); - fvph_pull_t[nSt]->Fill(tRes / dtHit); - - fvph_res_x_vs_x[nSt]->Fill(xHit, xRes); - fvph_res_y_vs_y[nSt]->Fill(yHit, yRes); - fvph_res_t_vs_t[nSt]->Fill(tHit, tRes); - - fvph_pull_x_vs_x[nSt]->Fill(xHit, xRes / dxHit); - fvph_pull_y_vs_y[nSt]->Fill(yHit, yRes / dyHit); - fvph_pull_t_vs_t[nSt]->Fill(tHit, tRes / dtHit); - } - } // loop over hits: end - - // Fill hit reconstruction efficiencies - if (IsMCUsed()) { - int nPoints = pMCInteractions->GetEntriesFast(); - - for (int iP = 0; iP < nPoints; ++iP) { - const auto* pMCPoint = dynamic_cast<const CbmTofInteraction*>(pMCInteractions->At(iP)); - LOG_IF(fatal, !pMCPoint) << fName << ": MC point does not exist for interaction " << iP; - //int address = pMCPoint->GetDetectorID(); - int iSt = GetStationID(pMCPoint); - LOG_IF(fatal, iSt < 0 || iSt >= nSt) - << fName << ": MC point for index = " << iP << " has wrong station index: iSt = " << iSt; - - double xMC = pMCPoint->GetX(); - double yMC = pMCPoint->GetY(); - double rMC = sqrt(xMC * xMC + yMC * yMC); - - // Conditions under which point is accounted as reconstructed: point - bool ifPointHasHits = (vNofHitsPerInteraction[iP] > 0); - - fvpe_reco_eff_vs_xy[iSt]->Fill(ifPointHasHits, xMC, yMC); - fvpe_reco_eff_vs_xy[nSt]->Fill(ifPointHasHits, xMC, yMC); - - fvpe_reco_eff_vs_r[iSt]->Fill(ifPointHasHits, rMC); - fvpe_reco_eff_vs_r[nSt]->Fill(ifPointHasHits, rMC); - } // iP - } // MC is used: end - - // Clear TClonesArray objects - if (IsMCUsed()) { - pHitInterMatches->Clear(); - pMCInteractions->Clear(); - } // IsMCUsed } // --------------------------------------------------------------------------------------------------------------------- // -void CbmCaInputQaTof::FillInteractions(std::shared_ptr<TClonesArray>& pvInters, - std::shared_ptr<TClonesArray>& pvMatches) +void CbmCaInputQaTof::FillHistogramsPerHit() { - std::unordered_map<ca::tools::LinkKey, int> mPointToInter; // Map of point (defined by a link key) to interaction - - int nHits = fpHits->GetEntriesFast(); - int nMCevents = (IsMCUsed()) ? fpMCEventList->GetNofEvents() : -1; - - // ----- Fill array of interactions - int iInter = -1; // Index of interaction (global over all MC events) - for (int iE = 0; iE < nMCevents; ++iE) { - int iFile = fpMCEventList->GetFileIdByIndex(iE); - int iEvent = fpMCEventList->GetEventIdByIndex(iE); - int nMCPoints = fpMCPoints->Size(iFile, iEvent); - - int iDprev = -1; // detector ID for previous point - int iTprev = -1; // track ID for previous point - for (int iP = 0; iP < nMCPoints; ++iP) { - const auto* pMCPoint = dynamic_cast<const CbmTofPoint*>(fpMCPoints->Get(iFile, iEvent, iP)); - LOG_IF(fatal, !pMCPoint) << fName << ": MC point is not defined for link (f,e,i) = " << iFile << ' ' << iEvent - << ' ' << iP; - - int iDcurr = pMCPoint->GetDetectorID(); // Current detector ID - int iTcurr = pMCPoint->GetTrackID(); // Current track ID - if (iDcurr != iDprev || iTcurr != iTprev) { - iInter++; - iDprev = iDcurr; - iTprev = iTcurr; - new ((*pvInters)[iInter]) CbmTofInteraction(); // Add empty interaction object - } - // Update current intercation with point, and updtate the map - auto* pCurrInter = static_cast<CbmTofInteraction*>(pvInters->At(iInter)); - pCurrInter->AddPoint(pMCPoint); - mPointToInter[ca::tools::LinkKey(iP, iEvent, iFile)] = iInter; - } // iP - } // iE - - // ----- Match interactions with hits - for (int iH = 0; iH < nHits; ++iH) { - const auto* pHit = dynamic_cast<const CbmTofHit*>(fpHits->At(iH)); - LOG_IF(fatal, !pHit) << fName << ": Hit is not defined for iH = " << iH; - - // FIXME: Is this cut still needed?? - int32_t address = pHit->GetAddress(); - if (0x00202806 == address || 0x00002806 == address) { continue; } - - auto* pHitInterMatch = new ((*pvMatches)[iH]) CbmMatch(); - - // Collect create interaction links from point links and fill the match - auto* pHitPointMatch = dynamic_cast<CbmMatch*>(fpHitMatches->At(iH)); - int nPointLinks = pHitPointMatch->GetNofLinks(); - for (int iL = 0; iL < nPointLinks; ++iL) { - const auto& link = pHitPointMatch->GetLink(iL); - int iFile = link.GetFile(); - int iEvent = link.GetEntry(); - int iP = link.GetIndex(); - double weight = link.GetWeight(); - int iInteraction = mPointToInter[ca::tools::LinkKey(iP, iEvent, iFile)]; - - // Create a TOF interaction link with a proper index - pHitInterMatch->AddLink(weight, iInteraction, iEvent, iFile); - - // Update interaction with information from point matched to the hit - // Since the averaged interaction in general may not be on the MC-trajectory of the particle, - // we decided to shift the interaction to the matched point. Here we assume, that there is only one - // matched point from the interaction. - auto* pPoint = static_cast<CbmTofPoint*>(fpMCPoints->Get(iFile, iEvent, iP)); - auto* pInter = static_cast<CbmTofInteraction*>(pvInters->At(iInteraction)); - pInter->SetFromPoint(pPoint); - } - } // iH -} - -// --------------------------------------------------------------------------------------------------------------------- -// -int CbmCaInputQaTof::GetStationID(const CbmTofPoint* pPoint) const -{ - double dist = 1000.; // NOTE: arbitrary large number - int iStSelect = -1; - // We select the station, which center is closest to the MC point - for (int iSt = 0; iSt < fpDetInterface->GetNtrackingStations(); ++iSt) { - if (std::fabs(pPoint->GetZ() - fpDetInterface->GetZref(iSt)) < dist) { - dist = std::fabs(pPoint->GetZ() - fpDetInterface->GetZref(iSt)); - iStSelect = iSt; - } + 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"; + + { // Check, if the hit is created on the one of the defined RPCs. If not, save to address into map + auto address = pHit->GetAddress(); + + auto iSmType = CbmTofAddress::GetSmType(address); + auto iSm = CbmTofAddress::GetSmId(address); + auto iRpc = CbmTofAddress::GetRpcId(address); + auto iCh = CbmTofAddress::GetChannelId(address); + double xHit = pHit->GetX(); + double yHit = pHit->GetY(); + double zHit = pHit->GetZ(); + LOG(info) << "Filling " << iSmType << ' ' << iSm << ' ' << iRpc << ' ' << iCh << " |-> " << xHit << ' ' << yHit << ' ' + << zHit; + fvph_hit_xy_vs_cell[iSmType][iSm][iRpc][iCh]->Fill(xHit, yHit); + fvph_hit_zx_vs_cell[iSmType][iSm][iRpc][iCh]->Fill(zHit, xHit); + fvph_hit_zy_vs_cell[iSmType][iSm][iRpc][iCh]->Fill(zHit, yHit); } - return iStSelect; } // --------------------------------------------------------------------------------------------------------------------- @@ -608,44 +132,12 @@ InitStatus CbmCaInputQaTof::InitDataBranches() // STS tracking detector interface fpDetInterface = CbmTofTrackingInterface::Instance(); - LOG_IF(fatal, !fpDetInterface) << "\033[1;31m" << fName << ": tracking detector interface is undefined\033[0m"; - - // FairRootManager - auto* pFairRootManager = FairRootManager::Instance(); - LOG_IF(fatal, !pFairRootManager) << "\033[1;31m" << fName << ": FairRootManager instance is a null pointer\033[0m"; - - // Time-slice - fpTimeSlice = dynamic_cast<CbmTimeSlice*>(pFairRootManager->GetObject("TimeSlice.")); - LOG_IF(fatal, !fpTimeSlice) << "\033[1;31m" << fName << ": time-slice branch is not found\033[0m"; - - // ----- Reconstructed data-branches initialization - - // Hits container - fpHits = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject("TofHit")); - LOG_IF(fatal, !fpHits) << "\033[1;31m" << fName << ": container of reconstructed hits in TOF is not found\033[0m"; - - // ----- MC-information branches initialization - if (IsMCUsed()) { - // MC manager - fpMCDataManager = dynamic_cast<CbmMCDataManager*>(pFairRootManager->GetObject("MCDataManager")); - LOG_IF(fatal, !fpMCDataManager) << "\033[1;31m" << fName << ": MC data manager branch is not found\033[0m"; - - // MC event list - fpMCEventList = dynamic_cast<CbmMCEventList*>(pFairRootManager->GetObject("MCEventList.")); - LOG_IF(fatal, !fpMCEventList) << "\033[1;31m" << fName << ": MC event list branch is not found\033[0m"; - - // MC tracks - fpMCTracks = fpMCDataManager->InitBranch("MCTrack"); - LOG_IF(fatal, !fpMCTracks) << "\033[1;31m" << fName << ": MC track branch is not found\033[0m"; - - // MC points - fpMCPoints = fpMCDataManager->InitBranch("TofPoint"); - LOG_IF(fatal, !fpMCTracks) << "\033[1;31m" << fName << ": MC point branch is not found for TOF\033[0m"; - - // Hit matches - fpHitMatches = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject("TofHitMatch")); - LOG_IF(fatal, !fpHitMatches) << "\033[1;31m]" << fName << ": hit match branch is not found for TOF\033[0m"; - } + auto baseInitStatus = CbmCaInputQaBase::InitDataBranches(); + if (kSUCCESS != baseInitStatus) { return baseInitStatus; } + + auto* pRuntimeDb = FairRunAna::Instance()->GetRuntimeDb(); + fDigiPar = dynamic_cast<CbmTofDigiPar*>(pRuntimeDb->getContainer("CbmTofDigiPar")); + fDigiBdfPar = dynamic_cast<CbmTofDigiBdfPar*>(pRuntimeDb->getContainer("CbmTofDigiBdfPar")); return kSUCCESS; } @@ -654,236 +146,56 @@ InitStatus CbmCaInputQaTof::InitDataBranches() // InitStatus CbmCaInputQaTof::InitCanvases() { - gStyle->SetOptFit(1); - int nSt = fpDetInterface->GetNtrackingStations(); - - // ******************** - // ** Drawing options - - // Contours - constexpr auto contColor = kOrange + 7; - constexpr auto contWidth = 2; // Line width - constexpr auto contStyle = 2; // Line style - constexpr auto contFill = 0; // Fill style - - // ********************************* - // ** Hit occupancies - - // ** Occupancy in XY plane - auto* pc_hit_ypos_vs_xpos = MakeQaObject<CbmQaCanvas>("hit_ypos_vs_xpos", "", 1600, 800); - pc_hit_ypos_vs_xpos->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_hit_ypos_vs_xpos->cd(iSt + 1); - fvph_hit_ypos_vs_xpos[iSt]->DrawCopy("colz", ""); - - // Square boarders of station, which are used for the field approximation - double stXmin = +fpDetInterface->GetXmax(iSt); - double stXmax = -fpDetInterface->GetXmax(iSt); - double stYmin = -fpDetInterface->GetYmax(iSt); - double stYmax = +fpDetInterface->GetYmax(iSt); - - auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - - // ----- Occupancy projections - auto* pc_hit_xpos = MakeQaObject<CbmQaCanvas>("hit_xpos", "", 1400, 800); - pc_hit_xpos->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_hit_xpos->cd(iSt + 1); - auto* phProj = fvph_hit_ypos_vs_xpos[iSt]->ProjectionY(); - phProj->DrawCopy(); - } - - auto* pc_hit_ypos = MakeQaObject<CbmQaCanvas>("hit_ypos", "", 1400, 800); - pc_hit_ypos->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_hit_ypos->cd(iSt + 1); - auto* phProj = fvph_hit_ypos_vs_xpos[iSt]->ProjectionX(); - phProj->DrawCopy(); - } - - - // ** Occupancy in XZ plane - MakeQaObject<CbmQaCanvas>("hit_xpos_vs_zpos", "", 600, 400); - fvph_hit_xpos_vs_zpos[nSt]->DrawCopy("colz", ""); - for (int iSt = 0; iSt < nSt; ++iSt) { - // Station positions in detector IFS - double stZmin = fpDetInterface->GetZmin(iSt); - double stZmax = fpDetInterface->GetZmax(iSt); - double stHmin = -fpDetInterface->GetXmax(iSt); - double stHmax = +fpDetInterface->GetXmax(iSt); - - auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - - // ** Occupancy in YZ plane - MakeQaObject<CbmQaCanvas>("hit_ypos_vs_zpos", "", 600, 400); - fvph_hit_ypos_vs_zpos[nSt]->DrawCopy("colz", ""); - for (int iSt = 0; iSt < nSt; ++iSt) { - // Station positions in detector IFS - double stZmin = fpDetInterface->GetZmin(iSt); - double stZmax = fpDetInterface->GetZmax(iSt); - double stHmin = -fpDetInterface->GetYmax(iSt); - double stHmax = +fpDetInterface->GetYmax(iSt); - - auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - - // ************ - // ************ - - if (IsMCUsed()) { - - // ********************** - // ** Point occupancies - - // ** Occupancy in XY plane - auto* pc_point_ypos_vs_xpos = MakeQaObject<CbmQaCanvas>("point_ypos_vs_xpos", "", 1600, 800); - pc_point_ypos_vs_xpos->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_point_ypos_vs_xpos->cd(iSt + 1); - fvph_point_ypos_vs_xpos[iSt]->DrawCopy("colz", ""); - - // Square boarders of station, which are used for the field approximation - double stXmin = +fpDetInterface->GetXmax(iSt); - double stXmax = -fpDetInterface->GetXmax(iSt); - double stYmin = -fpDetInterface->GetYmax(iSt); - double stYmax = +fpDetInterface->GetYmax(iSt); - - auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - - // ** Occupancy in XZ plane - MakeQaObject<CbmQaCanvas>("point_xpos_vs_zpos", "", 600, 400); - fvph_point_xpos_vs_zpos[nSt]->SetStats(false); - fvph_point_xpos_vs_zpos[nSt]->DrawCopy("colz", ""); - for (int iSt = 0; iSt < nSt; ++iSt) { - // Station positions in detector IFS - double stZmin = fpDetInterface->GetZmin(iSt); - double stZmax = fpDetInterface->GetZmax(iSt); - double stHmin = -fpDetInterface->GetXmax(iSt); - double stHmax = +fpDetInterface->GetXmax(iSt); - - auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); + auto baseInitStatus = CbmCaInputQaBase::InitCanvases(); + if (kSUCCESS != baseInitStatus) { return baseInitStatus; } + + // Hit occupancy vs TOF cell + { + auto DrawBox = [&] (double xLo, double yLo, double xUp, double yUp) { + auto* pBox = new TBox(xLo, yLo, xUp, yUp); + pBox->SetLineWidth(kOrange + 7); + pBox->SetLineStyle(2); + pBox->SetLineColor(2); + pBox->SetFillStyle(0); pBox->Draw("SAME"); - } - - // ** Occupancy in YZ plane - MakeQaObject<CbmQaCanvas>("point_ypos_vs_zpos", "", 600, 400); - fvph_point_ypos_vs_zpos[nSt]->SetStats(false); - fvph_point_ypos_vs_zpos[nSt]->DrawCopy("colz", ""); - for (int iSt = 0; iSt < nSt; ++iSt) { - // Station positions in detector IFS - double stZmin = fpDetInterface->GetZmin(iSt); - double stZmax = fpDetInterface->GetZmax(iSt); - double stHmin = -fpDetInterface->GetYmax(iSt); - double stHmax = +fpDetInterface->GetYmax(iSt); - - auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - - // ************** - // ** Residuals - - // x-coordinate - auto* pc_res_x = MakeQaObject<CbmQaCanvas>("res_x", "Residuals for x coordinate", 1600, 800); - pc_res_x->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_res_x->cd(iSt + 1); - fvph_res_x[iSt]->DrawCopy("", ""); - } - - // y-coordinate - auto* pc_res_y = MakeQaObject<CbmQaCanvas>("res_y", "Residuals for y coordinate", 1600, 800); - pc_res_y->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_res_y->cd(iSt + 1); - fvph_res_y[iSt]->DrawCopy("", ""); - } - - // time - auto* pc_res_t = MakeQaObject<CbmQaCanvas>("res_t", "Residuals for time", 1600, 800); - pc_res_t->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_res_t->cd(iSt + 1); - fvph_res_t[iSt]->DrawCopy("", ""); - } - - - // ********** - // ** Pulls - - // x-coordinate - auto* pc_pull_x = MakeQaObject<CbmQaCanvas>("pull_x", "Pulls for x coordinate", 1600, 800); - pc_pull_x->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_pull_x->cd(iSt + 1); - fvph_pull_x[iSt]->DrawCopy("", ""); - } - - // y-coordinate - auto* pc_pull_y = MakeQaObject<CbmQaCanvas>("pull_y", "Pulls for y coordinate", 1600, 800); - pc_pull_y->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_pull_y->cd(iSt + 1); - fvph_pull_y[iSt]->DrawCopy("", ""); - } - - // time - auto* pc_pull_t = MakeQaObject<CbmQaCanvas>("pull_t", "Pulls for time", 1600, 800); - pc_pull_t->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_pull_t->cd(iSt + 1); - fvph_pull_t[iSt]->DrawCopy("", ""); - } - - - // ************************************ - // ** Hit reconstruction efficiencies + }; - auto* pc_reco_eff_vs_r = - MakeQaObject<CbmQaCanvas>("reco_eff_vs_r", "Hit efficiencies wrt distance from center", 1600, 800); - pc_reco_eff_vs_r->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_reco_eff_vs_r->cd(iSt + 1); - fvpe_reco_eff_vs_r[iSt]->SetMarkerStyle(20); - fvpe_reco_eff_vs_r[iSt]->DrawCopy("AP", ""); - } - auto* pc_reco_eff_vs_xy = MakeQaObject<CbmQaCanvas>("reco_eff_vs_xy", "Hit efficiencies wrt x and y", 1600, 800); - pc_reco_eff_vs_xy->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_reco_eff_vs_xy->cd(iSt + 1); - fvpe_reco_eff_vs_xy[iSt]->DrawCopy("colz", ""); + for (int iSmType = 0; iSmType < fDigiBdfPar->GetNbSmTypes(); ++iSmType) { + if (iSmType == 5) { continue; } // skip T0 + for (int iSm = 0; iSm < fDigiBdfPar->GetNbSm(iSmType); ++iSm) { + for (int iRpc = 0; iRpc < fDigiBdfPar->GetNbRpc(iSmType); ++iRpc) { + for (int iCh = 0; iCh < fDigiBdfPar->GetNbChan(iSmType, iRpc); ++iCh) { + const char* dir = "occup_cell/"; + TString name = Form("%s/occup_xy_smt%d_sm%d_rpc%d_ch%d", dir, iSmType, iSm, iRpc, iCh); + + auto address = CbmTofAddress::GetUniqueAddress(iSm, iRpc, iCh, /*side = */ 0, iSmType); + const auto* pCell = dynamic_cast<const CbmTofCell*>(fDigiPar->GetCell(address)); + auto xLo = pCell->GetX() - 0.5 * pCell->GetSizex(); + auto xUp = pCell->GetX() + 0.5 * pCell->GetSizex(); + auto yLo = pCell->GetY() - 0.5 * pCell->GetSizey(); + auto yUp = pCell->GetY() + 0.5 * pCell->GetSizey(); + auto zLo = pCell->GetZ() - 1.0; + auto zUp = pCell->GetZ() + 1.0; + + auto* canv = MakeQaObject<CbmQaCanvas>(name, "", 1500, 800); + canv->Divide(3,1); + { + canv->cd(1); + fvph_hit_xy_vs_cell[iSmType][iSm][iRpc][iCh]->DrawCopy("colz", ""); + DrawBox(xLo, yLo, xUp, yUp); + + canv->cd(2); + fvph_hit_zx_vs_cell[iSmType][iSm][iRpc][iCh]->DrawCopy("colz", ""); + DrawBox(zLo, xLo, zUp, xUp); + + canv->cd(3); + fvph_hit_zy_vs_cell[iSmType][iSm][iRpc][iCh]->DrawCopy("colz", ""); + DrawBox(zLo, yLo, zUp, yUp); + } + } + } + } } } @@ -894,172 +206,58 @@ InitStatus CbmCaInputQaTof::InitCanvases() // InitStatus CbmCaInputQaTof::InitHistograms() { - int nSt = fpDetInterface->GetNtrackingStations(); - - // ----- Histograms initialization - fvph_hit_ypos_vs_xpos.resize(nSt + 1, nullptr); - fvph_hit_ypos_vs_zpos.resize(nSt + 1, nullptr); - fvph_hit_xpos_vs_zpos.resize(nSt + 1, nullptr); - - fvph_hit_station_delta_z.resize(nSt); - - fvph_hit_dx.resize(nSt + 1, nullptr); - fvph_hit_dy.resize(nSt + 1, nullptr); - fvph_hit_dt.resize(nSt + 1, nullptr); - - for (int iSt = 0; iSt <= nSt; ++iSt) { - TString nsuff = (iSt == nSt) ? "" : Form("_st%d", iSt); // Histogram name suffix - TString tsuff = (iSt == nSt) ? "" : Form(" in TOF station %d", iSt); // Histogram title suffix - TString sN = ""; - TString sT = ""; - - // Hit occupancy - sN = (TString) "hit_ypos_vs_xpos" + nsuff; - sT = (TString) "Hit occupancy in xy-Plane" + tsuff + ";x_{hit} [cm];y_{hit} [cm]"; - fvph_hit_ypos_vs_xpos[iSt] = MakeQaObject<TH2F>(sN, sT, kNbins, kRHitX[0], kRHitX[1], kNbins, kRHitY[0], kRHitY[1]); - - sN = (TString) "hit_xpos_vs_zpos" + nsuff; - sT = (TString) "Hit occupancy in xz-Plane" + tsuff + ";z_{hit} [cm];x_{hit} [cm]"; - fvph_hit_xpos_vs_zpos[iSt] = - MakeQaObject<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitX[0], kRHitX[1]); - - sN = (TString) "hit_ypos_vs_zpos" + nsuff; - sT = (TString) "Hit occupancy in yz-plane" + tsuff + ";z_{hit} [cm];y_{hit} [cm]"; - fvph_hit_ypos_vs_zpos[iSt] = - MakeQaObject<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitY[0], kRHitY[1]); - - // Hit errors - sN = (TString) "hit_dx" + nsuff; - sT = (TString) "Hit position error along x-axis" + tsuff + ";dx_{hit} [cm]"; - fvph_hit_dx[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, kRHitDx[0], kRHitDx[1]); - - sN = (TString) "hit_dy" + nsuff; - sT = (TString) "Hit position error along y-axis" + tsuff + ";dy_{hit} [cm]"; - fvph_hit_dy[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, kRHitDy[0], kRHitDy[1]); - - sN = (TString) "hit_dt" + nsuff; - sT = (TString) "Hit time error" + tsuff + ";dt_{hit} [ns]"; - fvph_hit_dt[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, kRHitDt[0], kRHitDt[1]); - - if (iSt == nSt) { continue; } // Following histograms are defined only for separate station - - sN = (TString) "hit_station_delta_z" + nsuff; - sT = (TString) "Different between hit and station z-positions" + tsuff + ";z_{hit} - z_{st} [cm]"; - fvph_hit_station_delta_z[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, -5., 5); - - } // loop over station index: end - - // ----- Initialize histograms, which are use MC-information - if (IsMCUsed()) { - // Resize histogram vectors - fvph_n_points_per_hit.resize(nSt + 1, nullptr); - fvph_point_ypos_vs_xpos.resize(nSt + 1, nullptr); - fvph_point_xpos_vs_zpos.resize(nSt + 1, nullptr); - fvph_point_ypos_vs_zpos.resize(nSt + 1, nullptr); - fvph_point_hit_delta_z.resize(nSt + 1, nullptr); - fvph_res_x.resize(nSt + 1, nullptr); - fvph_res_y.resize(nSt + 1, nullptr); - fvph_res_t.resize(nSt + 1, nullptr); - fvph_pull_x.resize(nSt + 1, nullptr); - fvph_pull_y.resize(nSt + 1, nullptr); - fvph_pull_t.resize(nSt + 1, nullptr); - fvph_res_x_vs_x.resize(nSt + 1, nullptr); - fvph_res_y_vs_y.resize(nSt + 1, nullptr); - fvph_res_t_vs_t.resize(nSt + 1, nullptr); - fvph_pull_x_vs_x.resize(nSt + 1, nullptr); - fvph_pull_y_vs_y.resize(nSt + 1, nullptr); - fvph_pull_t_vs_t.resize(nSt + 1, nullptr); - fvpe_reco_eff_vs_xy.resize(nSt + 1, nullptr); - fvpe_reco_eff_vs_r.resize(nSt + 1, nullptr); - - for (int iSt = 0; iSt <= nSt; ++iSt) { - TString nsuff = (iSt == nSt) ? "" : Form("_st%d", iSt); // Histogram name suffix - TString tsuff = (iSt == nSt) ? "" : Form(" in TOF station %d", iSt); // Histogram title suffix - TString sN = ""; - TString sT = ""; - - sN = (TString) "n_points_per_hit" + nsuff; - sT = (TString) "Number of points per hit" + tsuff + ";N_{point}/hit"; - fvph_n_points_per_hit[iSt] = MakeQaObject<TH1F>(sN, sT, 10, -0.5, 9.5); - - // Point occupancy - sN = (TString) "point_ypos_vs_xpos" + nsuff; - sT = (TString) "Point occupancy in XY plane" + tsuff + ";x_{MC} [cm];y_{MC} [cm]"; - fvph_point_ypos_vs_xpos[iSt] = - MakeQaObject<TH2F>(sN, sT, kNbins, kRHitX[0], kRHitX[1], kNbins, kRHitY[0], kRHitY[1]); - - sN = (TString) "point_xpos_vs_zpos" + nsuff; - sT = (TString) "Point Occupancy in XZ plane" + tsuff + ";z_{MC} [cm];x_{MC} [cm]"; - fvph_point_xpos_vs_zpos[iSt] = - MakeQaObject<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitX[0], kRHitX[1]); - - sN = (TString) "point_ypos_vs_zpos" + nsuff; - sT = (TString) "Point Occupancy in YZ Plane" + tsuff + ";z_{MC} [cm];y_{MC} [cm]"; - fvph_point_ypos_vs_zpos[iSt] = - MakeQaObject<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitY[0], kRHitY[1]); - - // Difference between MC input z and hit z coordinates - sN = (TString) "point_hit_delta_z" + nsuff; - sT = (TString) "Distance between point and hit along z axis" + tsuff + ";z_{reco} - z_{MC} [cm]"; - fvph_point_hit_delta_z[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, -0.04, 0.04); - - sN = (TString) "res_x" + nsuff; - sT = (TString) "Residuals for x-coordinate" + tsuff + ";x_{reco} - x_{MC} [cm]"; - fvph_res_x[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, kRResX[0], kRResX[1]); - - sN = (TString) "res_y" + nsuff; - sT = (TString) "Residuals for y-coordinate" + tsuff + ";y_{reco} - y_{MC} [cm]"; - fvph_res_y[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, kRResY[0], kRResY[1]); - - sN = (TString) "res_t" + nsuff; - sT = (TString) "Residuals for time" + tsuff + ";t_{reco} - t_{MC} [ns]"; - fvph_res_t[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, kRResT[0], kRResT[1]); - - sN = (TString) "pull_x" + nsuff; - sT = (TString) "Pulls for x-coordinate" + tsuff + ";(x_{reco} - x_{MC}) / #sigma_{x}^{reco}"; - fvph_pull_x[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, kRPullX[0], kRPullX[1]); - - sN = (TString) "pull_y" + nsuff; - sT = (TString) "Pulls for y-coordinate" + tsuff + ";(y_{reco} - y_{MC}) / #sigma_{y}^{reco}"; - fvph_pull_y[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, kRPullY[0], kRPullY[1]); - - sN = (TString) "pull_t" + nsuff; - sT = (TString) "Pulls for time" + tsuff + ";(t_{reco} - t_{MC}) / #sigma_{t}^{reco}"; - fvph_pull_t[iSt] = MakeQaObject<TH1F>(sN, sT, kNbins, kRPullT[0], kRPullT[1]); - - sN = (TString) "res_x_vs_x" + nsuff; - sT = (TString) "Residuals for x-coordinate" + tsuff + ";x_{reco} [cm];x_{reco} - x_{MC} [cm]"; - fvph_res_x_vs_x[iSt] = MakeQaObject<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResX[0], kRResX[1]); - - sN = (TString) "res_y_vs_y" + nsuff; - sT = (TString) "Residuals for y-coordinate" + tsuff + ";y_{reco} [cm];y_{reco} - y_{MC} [cm]"; - fvph_res_y_vs_y[iSt] = MakeQaObject<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResY[0], kRResY[1]); - - sN = (TString) "res_t_vs_t" + nsuff; - sT = (TString) "Residuals for time" + tsuff + "t_{reco} [ns];t_{reco} - t_{MC} [ns]"; - fvph_res_t_vs_t[iSt] = MakeQaObject<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResT[0], kRResT[1]); - - sN = (TString) "pull_x_vs_x" + nsuff; - sT = (TString) "Pulls for x-coordinate" + tsuff + "x_{reco} [cm];(x_{reco} - x_{MC}) / #sigma_{x}^{reco}"; - fvph_pull_x_vs_x[iSt] = MakeQaObject<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullX[0], kRPullX[1]); - - sN = (TString) "pull_y_vs_y" + nsuff; - sT = (TString) "Pulls for y-coordinate" + tsuff + ";y_{reco} [cm];(y_{reco} - y_{MC}) / #sigma_{y}^{reco}"; - fvph_pull_y_vs_y[iSt] = MakeQaObject<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullY[0], kRPullY[1]); - - sN = (TString) "pull_t_vs_t" + nsuff; - sT = (TString) "Pulls for time" + tsuff + ";t_{reco} [ns];(t_{reco} - t_{MC}) / #sigma_{t}^{reco}"; - fvph_pull_t_vs_t[iSt] = MakeQaObject<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullT[0], kRPullT[1]); - - - sN = (TString) "reco_eff_vs_xy" + nsuff; - sT = (TString) "Hit rec. efficiency" + tsuff + ";x_{MC} [cm];y_{MC} [cm];#epsilon"; - fvpe_reco_eff_vs_xy[iSt] = MakeQaObject<CbmQaEff>(sN, sT, 100, -50, 50, 100, -50, 50); - - sN = (TString) "reco_eff_vs_r" + nsuff; - sT = (TString) "Hit rec. efficiency" + tsuff + ";r_{MC} [cm];#epsilon"; - fvpe_reco_eff_vs_r[iSt] = MakeQaObject<CbmQaEff>(sN, sT, 100, 0, 100); + auto baseInitStatus = CbmCaInputQaBase::InitHistograms(); + if (kSUCCESS != baseInitStatus) { return baseInitStatus; } + + // Hit occupancy vs. TOF cell + MakeQaDirectory("occup_cell/"); + int nSmTypes = fDigiBdfPar->GetNbSmTypes(); + fvph_hit_xy_vs_cell.resize(nSmTypes); + fvph_hit_zx_vs_cell.resize(nSmTypes); + fvph_hit_zy_vs_cell.resize(nSmTypes); + for (int iSmType = 0; iSmType < nSmTypes; ++iSmType) { + if (iSmType == 5) { continue; } // skip T0 + MakeQaDirectory(Form("occup_cell/sm_type_%d", iSmType)); + int nSm = fDigiBdfPar->GetNbSm(iSmType); + int nRpc = fDigiBdfPar->GetNbRpc(iSmType); + + fvph_hit_xy_vs_cell[iSmType].resize(nSm); + fvph_hit_zx_vs_cell[iSmType].resize(nSm); + fvph_hit_zy_vs_cell[iSmType].resize(nSm); + for (int iSm = 0; iSm < fDigiBdfPar->GetNbSm(iSmType); ++iSm) { + MakeQaDirectory(Form("occup_cell/sm_type_%d/sm_%d/", iSmType, iSm)); + fvph_hit_xy_vs_cell[iSmType][iSm].resize(nRpc); + fvph_hit_zx_vs_cell[iSmType][iSm].resize(nRpc); + fvph_hit_zy_vs_cell[iSmType][iSm].resize(nRpc); + for (int iRpc = 0; iRpc < nRpc; ++iRpc) { + int nCh = fDigiBdfPar->GetNbChan(iSmType, iRpc); + fvph_hit_xy_vs_cell[iSmType][iSm][iRpc].resize(nCh, nullptr); + fvph_hit_zx_vs_cell[iSmType][iSm][iRpc].resize(nCh, nullptr); + fvph_hit_zy_vs_cell[iSmType][iSm][iRpc].resize(nCh, nullptr); + for (int iCh = 0; iCh < nCh; ++iCh) { + const char* dir = Form("occup_cell/sm_type_%d/sm_%d/rpc%d_ch%d", iSmType, iSm, iRpc, iCh); + TString name = Form("%s/occup_xy_smt%d_sm%d_rpc%d_ch%d", dir, iSmType, iSm, iRpc, iCh); + TString title = Form("Hit Occupancy in xy-Plane for iSmType = %d, iSm = %d, iRpc = %d, iCh = %d", + iSmType, iSm, iRpc, iCh); + title += ";x_{hit} [cm];y_{hit} [cm]"; + fvph_hit_xy_vs_cell[iSmType][iSm][iRpc][iCh] + = MakeQaObject<TH2F>(name, title, fNbXo, fLoXo, fUpXo, fNbYo, fLoYo, fUpYo); + name = Form("%s/occup_zx_smt%d_sm%d_rpc%d_ch%d", dir, iSmType, iSm, iRpc, iCh); + title = Form("Hit Occupancy in zx-Plane for iSmType = %d, iSm = %d, iRpc = %d, iCh = %d", + iSmType, iSm, iRpc, iCh); + title += ";z_{hit} [cm];x_{hit} [cm]"; + fvph_hit_zx_vs_cell[iSmType][iSm][iRpc][iCh] + = MakeQaObject<TH2F>(name, title, fNbinsZ, frZmin.back(), frZmax.back(), fNbXo, fLoXo, fUpXo); + name = Form("%s/occup_zy_smt%d_sm%d_rpc%d_ch%d", dir, iSmType, iSm, iRpc, iCh); + title = Form("Hit Occupancy in zy-Plane for iSmType = %d, iSm = %d, iRpc = %d", iSmType, iSm, iRpc); + title += ";z_{hit} [cm];y_{hit} [cm]"; + fvph_hit_zy_vs_cell[iSmType][iSm][iRpc][iCh] + = MakeQaObject<TH2F>(name, title, fNbinsZ, frZmin.back(), frZmax.back(), fNbYo, fLoYo, fUpYo); + } + } } } + + return kSUCCESS; } diff --git a/reco/L1/qa/CbmCaInputQaTof.h b/reco/L1/qa/CbmCaInputQaTof.h index 02062816a7586b642e77dfc395bb64e762c3d0e9..58b440a149636cae6526fdb3a84a1944ff1e23d1 100644 --- a/reco/L1/qa/CbmCaInputQaTof.h +++ b/reco/L1/qa/CbmCaInputQaTof.h @@ -14,12 +14,16 @@ #include "CbmCaInputQaBase.h" #include "CbmMCDataManager.h" #include "CbmQaTask.h" +#include "CbmTofCell.h" #include "TMath.h" -#include <set> +#include <tuple> #include <unordered_map> #include <vector> +#include <string> +#include <sstream> +#include <iomanip> class CbmMCEventList; class CbmMCDataArray; @@ -32,147 +36,68 @@ class CbmTofTrackingInterface; class TClonesArray; class TH1F; class CbmQaEff; +class CbmTofDigiPar; +class CbmTofDigiBdfPar; /// A QA-task class, which provides assurance of TOF hits and geometry /// -class CbmCaInputQaTof : public CbmQaTask, public CbmCaInputQaBase { +class CbmCaInputQaTof : public CbmCaInputQaBase<L1DetectorID::kTof> { public: - /// Constructor from parameters - /// \param verbose Verbose level - /// \param isMCUsed Flag, whether MC information is available for this task + /// @brief Constructor from parameters + /// @param verbose Verbose level + /// @param isMCUsed Flag, whether MC information is available for this task CbmCaInputQaTof(int verbose, bool isMCUsed); - ClassDef(CbmCaInputQaTof, 0); - protected: - // ******************************************** - // ** Virtual method override from CbmQaTask ** - // ******************************************** - - /// Checks results of the QA - /// \return Success flag + /// @brief Checks results of the QA + /// @return Success flag bool Check(); - /// Initializes data branches + /// @brief Initializes data branches InitStatus InitDataBranches(); - /// Initializes canvases + /// @brief Initializes canvases InitStatus InitCanvases(); - /// Initializes histograms + /// @brief Initializes histograms InitStatus InitHistograms(); - /// Fills histograms per event or time-slice - void FillHistograms(); + /// @brief Fills histograms per event or time-slice + void FillHistograms() { CbmCaInputQaBase::FillHistograms(); } - /// De-initializes histograms + /// @brief De-initializes histograms void DeInit(); -private: - // ********************* - // ** Private methods ** - // ********************* - - /// Gets index of station by pointer to MC point - /// \param pPoint Pointer to TOF MC point - int GetStationID(const CbmTofPoint* pPoint) const; - - /// Fill vector of interactions and corresponding match objects - /// \param pvInters Pointer to interactions array - /// \param pvMatches Pointer to hit->interaction matches - void FillInteractions(std::shared_ptr<TClonesArray>& pvInters, std::shared_ptr<TClonesArray>& pvMatches); - - // ********************* - // ** Class variables ** - // ********************* - - - // ** Data branches ** - CbmTofTrackingInterface* fpDetInterface = nullptr; ///< Instance of detector interface - - CbmTimeSlice* fpTimeSlice = nullptr; ///< Pointer to current time-slice - - TClonesArray* fpHits = nullptr; ///< Array of hits - - CbmMCDataManager* fpMCDataManager = nullptr; ///< MC data manager - CbmMCEventList* fpMCEventList = nullptr; ///< MC event list - CbmMCDataArray* fpMCTracks = nullptr; ///< Array of MC tracks - CbmMCDataArray* fpMCPoints = nullptr; ///< Array of MC points - - TClonesArray* fpHitMatches = nullptr; ///< Array of hit matches - - // ** Histograms ** - - static constexpr double kEffRangeMin = 10.; ///< Lower limit of hit distance for efficiency integration [cm] - static constexpr double kEffRangeMax = 30.; ///< Upper limit of hit distance for efficiency integration [cm] - - static constexpr double kRHitX[2] = {-200., 200.}; ///< Range for hit x coordinate [cm] - static constexpr double kRHitY[2] = {-200., 200.}; ///< Range for hit y coordinate [cm] - static constexpr double kRHitZ[2] = {750., 850.}; ///< Range for hit z coordinate [cm] - - static constexpr int kNbins = 200; ///< General number of bins - static constexpr int kNbinsZ = 800; ///< Number of bins for z coordinate axis + /// @brief Defines parameters of the task + void DefineParameters(); - static constexpr double kRHitDx[2] = {-.05, .005}; ///< Range for hit x coordinate [cm] - static constexpr double kRHitDy[2] = {-.05, .005}; ///< Range for hit y coordinate [cm] - static constexpr double kRHitDt[2] = {-10., 10.}; ///< Range for hit time [ns] + /// @brief Fills histograms per hit + void FillHistogramsPerHit(); + + /// @brief Fills histograms per MC point + void FillHistogramsPerPoint() {} - static constexpr double kRResX[2] = {-2., 2.}; ///< Range for residual of x coordinate [cm] - static constexpr double kRResY[2] = {-4., 4.}; ///< Range for residual of y coordinate [cm] - static constexpr double kRResT[2] = {-.5, .5}; ///< Range for residual of time [ns] - - static constexpr double kRPullX[2] = {-10., 10.}; ///< Range for pull of x coordinate - static constexpr double kRPullY[2] = {-10., 10.}; ///< Range for pull of y coordinate - static constexpr double kRPullT[2] = {-5., 5.}; ///< Range for pull of time - - - // NOTE: the last element of each vector stands for integral distribution over all stations - - // hit occupancy - std::vector<TH2F*> fvph_hit_ypos_vs_xpos; ///< hit ypos vs xpos in different stations - std::vector<TH2F*> fvph_hit_xpos_vs_zpos; ///< hit xpos vs zpos in different stations - std::vector<TH2F*> fvph_hit_ypos_vs_zpos; ///< hit ypos vs zpos in different stations - - // difference between hit and station z positions - std::vector<TH1F*> fvph_hit_station_delta_z; ///< Difference between station and hit z positions [cm] - - // hit errors - std::vector<TH1F*> fvph_hit_dx; - std::vector<TH1F*> fvph_hit_dy; - std::vector<TH1F*> fvph_hit_dt; - - // MC points occupancy - std::vector<TH1F*> fvph_n_points_per_hit; ///< number of points per one hit in different stations - - - std::vector<TH2F*> fvph_point_ypos_vs_xpos; ///< point ypos vs xpos in different stations - std::vector<TH2F*> fvph_point_xpos_vs_zpos; ///< point xpos vs zpos in different stations - std::vector<TH2F*> fvph_point_ypos_vs_zpos; ///< point ypos vs zpos in different stations - - std::vector<TH1F*> fvph_point_hit_delta_z; ///< difference between z of hit and MC point in different stations - - // Residuals - std::vector<TH1F*> fvph_res_x; ///< residual for x coordinate in different stations - std::vector<TH1F*> fvph_res_y; ///< residual for y coordinate in different stations - std::vector<TH1F*> fvph_res_t; ///< residual for t coordinate in different stations - - std::vector<TH2F*> fvph_res_x_vs_x; ///< residual for x coordinate in different stations - std::vector<TH2F*> fvph_res_y_vs_y; ///< residual for y coordinate in different stations - std::vector<TH2F*> fvph_res_t_vs_t; ///< residual for t coordinate in different stations - - // Pulls - std::vector<TH1F*> fvph_pull_x; ///< pull for x coordinate in different stations - std::vector<TH1F*> fvph_pull_y; ///< pull for y coordinate in different stations - std::vector<TH1F*> fvph_pull_t; ///< pull for t coordinate in different stations - - std::vector<TH2F*> fvph_pull_x_vs_x; ///< pull for x coordinate in different stations - std::vector<TH2F*> fvph_pull_y_vs_y; ///< pull for y coordinate in different stations - std::vector<TH2F*> fvph_pull_t_vs_t; ///< pull for t coordinate in different stations +private: + /// @brief Fills channel info map + void FillChannelInfoMap(); + + CbmTofDigiPar* fDigiPar = nullptr; + CbmTofDigiBdfPar* fDigiBdfPar = nullptr; + + int fNbXo = 400; ///< Number of bins in occupancy + int fNbYo = 400; ///< Number of bins in occupancy + int fUpXo = +100; ///< X max in occupancy [cm] + int fLoXo = -100; ///< X min in occupancy [cm] + int fUpYo = +100; ///< Y max in occupancy [cm] + int fLoYo = -100; ///< Y min in occupancy [cm] + template <typename T> + using Vector4D_t = std::vector<std::vector<std::vector<std::vector<T>>>>; + Vector4D_t<TH2F*> fvph_hit_xy_vs_cell; ///< Hit occupancy vs. TOF cell ix XY-plane + Vector4D_t<TH2F*> fvph_hit_zx_vs_cell; ///< Hit occupancy vs. TOF cell in ZX-plane + Vector4D_t<TH2F*> fvph_hit_zy_vs_cell; ///< Hit occupancy vs. TOF cell in ZY-plane - // Hit efficiencies - std::vector<CbmQaEff*> fvpe_reco_eff_vs_xy; ///< Efficiency of hit reconstruction vs. x and y coordinates of MC - std::vector<CbmQaEff*> fvpe_reco_eff_vs_r; ///< Efficiency of hit reconstruction vs. distance from center + ClassDef(CbmCaInputQaTof, 0); }; #endif // CbmCaInputQaTof_h diff --git a/reco/L1/qa/CbmCaInputQaTrd.cxx b/reco/L1/qa/CbmCaInputQaTrd.cxx index c9efb402e48a18f30ad82015f2b7326d22054471..22d13ada4f7b1bfbcf65b75947b5421cf397522c 100644 --- a/reco/L1/qa/CbmCaInputQaTrd.cxx +++ b/reco/L1/qa/CbmCaInputQaTrd.cxx @@ -2,227 +2,29 @@ SPDX-License-Identifier: GPL-3.0-only Authors: Sergei Zharko [committer] */ -/// \file CbmCaInputQaTrd.cxx -/// \date 13.01.2023 -/// \brief QA-task for CA tracking input from TRD detector (implementation) -/// \author S.Zharko <s.zharko@gsi.de> +/// @file CbmCaInputQaTrd.cxx +/// @date 13.01.2023 +/// @brief QA-task for CA tracking input from TRD detector (implementation) +/// @author S.Zharko <s.zharko@gsi.de> #include "CbmCaInputQaTrd.h" - -#include "CbmAddress.h" -#include "CbmMCDataArray.h" -#include "CbmMCEventList.h" -#include "CbmMCTrack.h" -#include "CbmMatch.h" -#include "CbmQaCanvas.h" -#include "CbmQaEff.h" -#include "CbmTimeSlice.h" -#include "CbmTrdHit.h" -#include "CbmTrdPoint.h" #include "CbmTrdTrackingInterface.h" -#include "FairRootManager.h" -#include "Logger.h" - -#include "TBox.h" -#include "TClonesArray.h" -#include "TEllipse.h" -#include "TF1.h" -#include "TFormula.h" -#include "TGraphAsymmErrors.h" -#include "TH1F.h" -#include "TH2F.h" -#include "TMath.h" -#include "TStyle.h" - -#include <algorithm> -#include <fstream> -#include <iomanip> -#include <numeric> - -#include "L1Constants.h" - ClassImp(CbmCaInputQaTrd); -namespace phys = L1Constants::phys; // from L1Constants.h - // --------------------------------------------------------------------------------------------------------------------- // -CbmCaInputQaTrd::CbmCaInputQaTrd(int verbose, bool isMCUsed) : CbmQaTask("CbmCaInputQaTrd", verbose, isMCUsed) {} +CbmCaInputQaTrd::CbmCaInputQaTrd(int verbose, bool isMCUsed) : CbmCaInputQaBase("CbmCaInputQaTrd", verbose, isMCUsed) +{ + // Default parameters of task + DefineParameters(); +} // --------------------------------------------------------------------------------------------------------------------- // bool CbmCaInputQaTrd::Check() { - bool res = true; - - int nSt = fpDetInterface->GetNtrackingStations(); - - - // ************************************************************** - // ** Basic checks, available both for real and simulated data ** - // ************************************************************** - - // ----- Checks for mismatches in the ordering of the stations - // - std::vector<double> vStationPos(nSt, 0.); - for (int iSt = 0; iSt < nSt; ++iSt) { - vStationPos[iSt] = fpDetInterface->GetZref(iSt); - } - - if (!std::is_sorted(vStationPos.cbegin(), vStationPos.cend(), [](int l, int r) { return l <= r; })) { - if (fVerbose > 0) { - LOG(error) << fName << ": stations are ordered improperly along the beam axis:"; - for (auto zPos : vStationPos) { - LOG(error) << "\t- " << zPos; - } - } - res = false; - } - - // ----- Checks for mismatch between station and hit z positions - // The purpose of this block is to be ensured, that hits belong to the correct tracking station. For each tracking - // station a unified position z-coordinate is defined, which generally differs from the corresponding positions of - // reconstructed hits. This is due to non-regular position along the beam axis for each detector sensor. Thus, the - // positions of the hits along z-axis are distributed relatively to the position of the station in some range. - // If hits goes out this range, it can signal about a mismatch between hits and geometry. For limits of the range - // one should select large enough values to cover the difference described above and in the same time small enough - // to avoid overlaps with the neighboring stations. - // For each station, a distribution of z_{hit} - z_{st} is integrated over the defined range and scaled by the - // total number of entries to the distribution. If this value is smaller then unity, then some of the hits belong to - // another station. - // - for (int iSt = 0; iSt < nSt; ++iSt) { - int nHits = fvph_hit_station_delta_z[iSt]->GetEntries(); - if (!nHits) { - LOG_IF(error, fVerbose > 0) << fName << ": station " << iSt << " does not have hits"; - res = false; - continue; - } - int iBinMin = fvph_hit_station_delta_z[iSt]->FindBin(-fMaxDiffZStHit); - int iBinMax = fvph_hit_station_delta_z[iSt]->FindBin(+fMaxDiffZStHit); - - if (fvph_hit_station_delta_z[iSt]->Integral(iBinMin, iBinMax) < nHits) { - LOG_IF(error, fVerbose > 0) << fName << ": station " << iSt << " has mismatches in hit z-positions"; - res = false; - } - } - - - // ******************************************************* - // ** Additional checks, if MC information is available ** - // ******************************************************* - - if (IsMCUsed()) { - // ----- Check efficiencies - // Error is raised, if any station has integrated efficiency lower then a selected threshold. - // Integrated efficiency is estimated with fitting the efficiency curve for hit distances (r) with a uniform - // distribution in the range from kEffRangeMin to kEffRangeMax, where the efficiency is assigned to be constant - // - // Fit efficiency curves - LOG(info) << "-- Hit efficiency integrated over hit distance from station center"; - LOG(info) << "\tintegration range: [" << fEffRange[0] << ", " << fEffRange[1] << "] cm"; - - auto* pEffTable = MakeQaObject<CbmQaTable>("eff_table", "Efficiency table", nSt, 2); - pEffTable->SetNamesOfCols({"Station ID", "Efficiency"}); - pEffTable->SetColWidth(20); - - for (int iSt = 0; iSt < nSt; ++iSt) { - auto eff = fvpe_reco_eff_vs_r[iSt]->GetMean(2); - pEffTable->SetCell(iSt, 0, iSt); - pEffTable->SetCell(iSt, 1, eff); - res = CheckRange("Hit finder efficiency in station " + std::to_string(iSt), eff, fEffThrsh, 1.000); - } - - LOG(info) << '\n' << pEffTable->ToString(3); - - // ----- Checks for residuals - // Check hits for potential biases from central values - - // Function to fit a residual distribution, returns a structure - auto FitResidualsAndGetMean = [&](TH1* pHist) { - auto fit = TF1("fitres", "gausn"); - double statMean = pHist->GetMean(); - double statSigm = pHist->GetStdDev(); - fit.SetParameters(100., statMean, statSigm); - pHist->Fit("fitres", "MQ"); - pHist->Fit("fitres", "MQ"); - pHist->Fit("fitres", "MQ"); - // NOTE: Several fit procedures are made to avoid empty fit results - std::array<double, 3> result; - result[0] = fit.GetParameter(1); - result[1] = -fit.GetParameter(2) * fResMeanThrsh; - result[2] = +fit.GetParameter(2) * fResMeanThrsh; - return result; - }; - - auto* pResidualsTable = - MakeQaObject<CbmQaTable>("residuals_mean", "Residual mean values in different stations", nSt, 4); - pResidualsTable->SetNamesOfCols({"Station ID", "Residual(x) [cm]", "Residual(y) [cm]", "Residual(t) [ns]"}); - pResidualsTable->SetColWidth(20); - - // Fill residuals fit results and check boundaries - for (int iSt = 0; iSt < nSt; ++iSt) { - auto fitX = FitResidualsAndGetMean(fvph_res_x[iSt]); - auto fitY = FitResidualsAndGetMean(fvph_res_y[iSt]); - auto fitT = FitResidualsAndGetMean(fvph_res_t[iSt]); - res = CheckRange("Mean of x residual [cm] in station " + std::to_string(iSt), fitX[0], fitX[1], fitX[2]) && res; - res = CheckRange("Mean of y residual [cm] in station " + std::to_string(iSt), fitY[0], fitY[1], fitY[2]) && res; - res = CheckRange("Mean of t residual [ns] in station " + std::to_string(iSt), fitT[0], fitT[1], fitT[2]) && res; - pResidualsTable->SetCell(iSt, 0, iSt); - pResidualsTable->SetCell(iSt, 1, fitX[0]); - pResidualsTable->SetCell(iSt, 2, fitX[1]); - pResidualsTable->SetCell(iSt, 3, fitX[2]); - } - - LOG(info) << '\n' << pResidualsTable->ToString(8); - - // ----- Checks for pulls - // For the hit pull is defined as a ration of the difference between hit and MC-point position or time component - // values and an error of the hit position (or time) component, calculated in the same z-positions. In ideal case, - // when the resolutions are well determined for detecting stations, the pull distributions should have a RMS equal - // to unity. Here we allow a RMS of the pull distributions to be varied in a predefined range. If the RMS runs out - // this range, QA task fails. - // TODO: Add checks for center - - // Fit pull distributions - - // ** Kaniadakis Gaussian distribution, gives smaller chi2 / NDF - //if (!gROOT->FindObject("Expk")) { - // new TFormula("Expk", "TMath::Power(TMath::Sqrt(1 + x[1] * x[1] * x[0] * x[0]) + x[0] * x[1], 1./x[1])"); - //} - //TF1* pPullFit = new TF1("pullFitGausn", "[0] * Expk(-[2] * (x - [1]) * (x - [1]), [3])", -10., 10.); - - auto FitPullsAndGetSigma = [&](TH1* pHist) { - auto fit = TF1("fitpull", "gausn(0)"); - fit.SetParameters(100, 0.001, 1.000); - pHist->Fit("fitpull", "Q"); - return fit.GetParameter(2); - }; - - auto* pPullsTable = MakeQaObject<CbmQaTable>("pulls_rms", "Pulls std. dev. values in different stations", nSt, 4); - pPullsTable->SetNamesOfCols({"Station ID", "Pull(x) sigm", "Pull(y) sigm", "Pull(z) sigm"}); - pPullsTable->SetColWidth(20); - - for (int iSt = 0; iSt < nSt; ++iSt) { - double rmsPullX = FitPullsAndGetSigma(fvph_pull_x[iSt]); - double rmsPullY = FitPullsAndGetSigma(fvph_pull_y[iSt]); - double rmsPullT = FitPullsAndGetSigma(fvph_pull_t[iSt]); - - double rmsPullHi = 1. + fPullWidthThrsh; - double rmsPullLo = 1. - fPullWidthThrsh; - - res = CheckRange("Rms for x pull in station " + std::to_string(iSt), rmsPullX, rmsPullLo, rmsPullHi) && res; - res = CheckRange("Rms for y pull in station " + std::to_string(iSt), rmsPullY, rmsPullLo, rmsPullHi) && res; - res = CheckRange("Rms for t pull in station " + std::to_string(iSt), rmsPullT, rmsPullLo, rmsPullHi) && res; - pPullsTable->SetCell(iSt, 0, iSt); - pPullsTable->SetCell(iSt, 1, rmsPullX); - pPullsTable->SetCell(iSt, 2, rmsPullY); - pPullsTable->SetCell(iSt, 3, rmsPullT); - } - - LOG(info) << '\n' << pPullsTable->ToString(3); - } + bool res = CbmCaInputQaBase::Check(); return res; } @@ -231,250 +33,44 @@ bool CbmCaInputQaTrd::Check() // void CbmCaInputQaTrd::DeInit() { - // Vectors with pointers to histograms - fvph_hit_ypos_vs_xpos.clear(); - fvph_hit_xpos_vs_zpos.clear(); - fvph_hit_ypos_vs_zpos.clear(); - - fvph_hit_station_delta_z.clear(); - - fvph_hit_dx.clear(); - fvph_hit_dy.clear(); - fvph_hit_dt.clear(); - - fvph_n_points_per_hit.clear(); - - fvph_point_ypos_vs_xpos.clear(); - fvph_point_xpos_vs_zpos.clear(); - fvph_point_ypos_vs_zpos.clear(); - - fvph_point_hit_delta_z.clear(); - - fvph_res_x.clear(); - fvph_res_y.clear(); - fvph_res_t.clear(); - - fvph_pull_x.clear(); - fvph_pull_y.clear(); - fvph_pull_t.clear(); - - fvph_res_x_vs_x.clear(); - fvph_res_y_vs_y.clear(); - fvph_res_t_vs_t.clear(); - - fvph_pull_x_vs_x.clear(); - fvph_pull_y_vs_y.clear(); - fvph_pull_t_vs_t.clear(); - - fvpe_reco_eff_vs_xy.clear(); - fvpe_reco_eff_vs_r.clear(); + CbmCaInputQaBase::DeInit(); } // --------------------------------------------------------------------------------------------------------------------- // -void CbmCaInputQaTrd::FillHistograms() +void CbmCaInputQaTrd::DefineParameters() { - int nSt = fpDetInterface->GetNtrackingStations(); - int nHits = fpHits->GetEntriesFast(); - int nMCevents = (IsMCUsed()) ? fpMCEventList->GetNofEvents() : -1; - - std::vector<std::vector<int>> vNofHitsPerPoint; // Number of hits per MC point in different MC events - - if (IsMCUsed()) { - vNofHitsPerPoint.resize(nMCevents); - for (int iE = 0; iE < nMCevents; ++iE) { - int iFile = fpMCEventList->GetFileIdByIndex(iE); - int iEvent = fpMCEventList->GetEventIdByIndex(iE); - int nPoints = fpMCPoints->Size(iFile, iEvent); - vNofHitsPerPoint[iE].resize(nPoints, 0); - } - } - - for (int iH = 0; iH < nHits; ++iH) { - const auto* pHit = dynamic_cast<const CbmTrdHit*>(fpHits->At(iH)); - LOG_IF(fatal, !pHit) << fName << ": hit with iH = " << iH << " is not an CbmTrdHit (dynamic cast failed)"; - - // ************************* - // ** Reconstructed hit QA - - // 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; - - //int hitType = pHit->GetClassType(); // TRD-1D, TRD-2D - - // Hit position - double xHit = pHit->GetX(); - double yHit = pHit->GetY(); - double zHit = pHit->GetZ(); - - double dxHit = pHit->GetDx(); - double dyHit = pHit->GetDy(); - - // Hit time - double tHit = pHit->GetTime(); - double dtHit = pHit->GetTimeError(); - - fvph_hit_ypos_vs_xpos[iSt]->Fill(xHit, yHit); - fvph_hit_xpos_vs_zpos[iSt]->Fill(zHit, xHit); - fvph_hit_ypos_vs_zpos[iSt]->Fill(zHit, yHit); - - fvph_hit_station_delta_z[iSt]->Fill(zHit - fpDetInterface->GetZref(iSt)); - - fvph_hit_dx[iSt]->Fill(dxHit); - fvph_hit_dy[iSt]->Fill(dyHit); - fvph_hit_dt[iSt]->Fill(dtHit); - - // ********************** - // ** MC information QA - - if (IsMCUsed()) { - CbmMatch* pHitMatch = dynamic_cast<CbmMatch*>(fpHitMatches->At(iH)); - LOG_IF(fatal, !pHitMatch) << fName << ": match object not found for hit ID " << iH; - - // Evaluate number of hits per one MC point - int nMCpoints = 0; // Number of MC points for this hit - - for (int iLink = 0; iLink < pHitMatch->GetNofLinks(); ++iLink) { - const auto& link = pHitMatch->GetLink(iLink); - - int iP = link.GetIndex(); // Index of MC point - - // Skip noisy links - if (iP < 0) { continue; } - - ++nMCpoints; - - int iE = fpMCEventList->GetEventIndex(link); - - 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 << ')'; - - LOG_IF(fatal, iP >= (int) vNofHitsPerPoint[iE].size()) - << fName << ": index of MC point is out of range (hit id = " << iH << ", link id = " << iLink - << ", event id = " << iE << ", point id = " << iP << ')'; - vNofHitsPerPoint[iE][iP]++; - } - - fvph_n_points_per_hit[iSt]->Fill(nMCpoints); - - if (nMCpoints != 1) { continue; } // ?? - - // The best link to in the match (probably, the cut on nMCpoints is meaningless) - const auto& bestPointLink = pHitMatch->GetMatchedLink(); - - // Skip noisy links - if (bestPointLink.GetIndex() < 0) { continue; } - - // Point matched by the best link - const auto* pMCPoint = dynamic_cast<const CbmTrdPoint*>(fpMCPoints->Get(bestPointLink)); - LOG_IF(fatal, !pMCPoint) << fName << ": MC point object does not exist for hit " << iH; - - // 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: "; - - double t0MC = fpMCEventList->GetEventTime(bestPointLink); - LOG_IF(fatal, t0MC < 0) << fName << ": MC time zero is lower then 0 ns: " << t0MC; - - // ----- MC point properties - // - double mass = pMCTrack->GetMass(); - //double charge = pMCTrack->GetCharge(); - //double pdg = pMCTrack->GetPdgCode(); - - // Entrance position and time - double xMC = pMCPoint->GetX(); - double yMC = pMCPoint->GetY(); - double zMC = pMCPoint->GetZ(); - double tMC = pMCPoint->GetTime() + t0MC; - - // MC point entrance momenta - double pxMC = pMCPoint->GetPx(); - double pyMC = pMCPoint->GetPy(); - double pzMC = pMCPoint->GetPz(); - double pMC = sqrt(pxMC * pxMC + pyMC * pyMC + pzMC * pzMC); - - // MC point exit momenta - // double pxMCo = pMCPoint->GetPxOut(); - // double pyMCo = pMCPoint->GetPyOut(); - // double pzMCo = pMCPoint->GetPzOut(); - // double pMCo = sqrt(pxMCo * pxMCo + pyMCo * pyMCo + pzMCo * pzMCo); - - // Position and time shifted to z-coordinate of the hit - double shiftZ = zHit - zMC; // Difference between hit and point z positions - double xMCs = xMC + shiftZ * pxMC / pzMC; - double yMCs = yMC + shiftZ * pyMC / pzMC; - double tMCs = tMC + shiftZ / (pzMC * phys::kSpeedOfLight) * sqrt(mass * mass + pMC * pMC); - - // Residuals - double xRes = xHit - xMCs; - double yRes = yHit - yMCs; - double tRes = tHit - tMCs; - - // ------ Cuts - - if (std::fabs(pMCPoint->GetPzOut()) < fMinMomentum) { continue; } // CUT ON MINIMUM MOMENTUM - //if (pMCo < cuts::kMinP) { continue; } // Momentum threshold - - - fvph_point_ypos_vs_xpos[iSt]->Fill(xMC, yMC); - fvph_point_xpos_vs_zpos[iSt]->Fill(zMC, xMC); - fvph_point_ypos_vs_zpos[iSt]->Fill(zMC, yMC); - - fvph_point_hit_delta_z[iSt]->Fill(shiftZ); - - fvph_res_x[iSt]->Fill(xRes); - fvph_res_y[iSt]->Fill(yRes); - fvph_res_t[iSt]->Fill(tRes); - - fvph_pull_x[iSt]->Fill(xRes / dxHit); - fvph_pull_y[iSt]->Fill(yRes / dyHit); - fvph_pull_t[iSt]->Fill(tRes / dtHit); - - fvph_res_x_vs_x[iSt]->Fill(xHit, xRes); - fvph_res_y_vs_y[iSt]->Fill(yHit, yRes); - fvph_res_t_vs_t[iSt]->Fill(tHit, tRes); - - fvph_pull_x_vs_x[iSt]->Fill(xHit, xRes / dxHit); - fvph_pull_y_vs_y[iSt]->Fill(yHit, yRes / dyHit); - fvph_pull_t_vs_t[iSt]->Fill(tHit, tRes / dtHit); - } - } // loop over hits: end - - // Fill hit reconstruction efficiencies - if (IsMCUsed()) { - for (int iE = 0; iE < nMCevents; ++iE) { - int iFile = fpMCEventList->GetFileIdByIndex(iE); - int iEvent = fpMCEventList->GetEventIdByIndex(iE); - int nPoints = fpMCPoints->Size(iFile, iEvent); - - for (int iP = 0; iP < nPoints; ++iP) { - const auto* pMCPoint = dynamic_cast<const CbmTrdPoint*>(fpMCPoints->Get(iFile, iEvent, iP)); - LOG_IF(fatal, !pMCPoint) << fName << ": MC point does not exist for iFile = " << iFile - << ", iEvent = " << iEvent << ", iP = " << iP; - int address = pMCPoint->GetDetectorID(); - int iSt = fpDetInterface->GetTrackingStationIndex(address); - LOG_IF(fatal, iSt < 0 || iSt >= nSt) - << fName << ": MC point for FEI = " << iFile << ", " << iEvent << ", " << iP << " and address " << address - << " has wrong station index: iSt = " << iSt; - - double xMC = pMCPoint->GetXIn(); - double yMC = pMCPoint->GetYIn(); - double rMC = sqrt(xMC * xMC + yMC * yMC); - - // Conditions under which point is accounted as reconstructed: point - bool ifPointHasHits = (vNofHitsPerPoint[iE][iP] > 0); - - fvpe_reco_eff_vs_xy[iSt]->Fill(xMC, yMC, ifPointHasHits); - fvpe_reco_eff_vs_r[iSt]->Fill(rMC, ifPointHasHits); + auto SetRange = [] (std::array<double, 2>& range, double min, double max) { range[0] = min; range[1] = max; }; + // Hit errors + SetRange(fRHitDx, 0.0000, 5.00); // [cm] + SetRange(fRHitDy, 0.0000, 5.00); // [cm] + SetRange(fRHitDu, 0.0000, 5.00); // [cm] + SetRange(fRHitDv, 0.0000, 5.00); // [cm] + SetRange(fRHitDt, 0.0000, 10.00); // [ns] + // Residuals + SetRange(fRResX, -2.00, 2.00); + SetRange(fRResY, -4.00, 4.00); + SetRange(fRResU, -2.00, 2.00); + SetRange(fRResV, -4.00, 4.00); + SetRange(fRResT, -0.50, 0.50); + // QA result selection criteria + SetRange(fEffRange, 10.0, 30.0); ///< Range for hit efficiency approximation + fResMeanThrsh = 0.50; ///< Maximum allowed deviation of residual mean from zero [sigma] + fPullMeanThrsh = 0.10; ///< Maximum allowed deviation of pull mean from zero + fPullWidthThrsh = 2.00; ///< Maximum allowed deviation of pull width from unity + fEffThrsh = 0.50; ///< Threshold for hit efficiency in the selected range + fMaxDiffZStHit = 1.00; ///< Maximum allowed difference between z-position of hit and station [cm] + fMinMomentum = 0.05; ///< Minimum momentum of particle [GeV/c] + // Track selection criteria + fTrackSelectionPrimary = true; ///< must the track come from the primary vertex + fTrackSelectionMomentum = 0.1; ///< track momentum GeV when it exits the tracking station + fTrackSelectionAngle = 60.; ///< track incl. angle [grad] when it exits the station +} - } // loop over MC-points: end - } // loop over MC-events: end - } // MC is used: end +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmCaInputQaTrd::FillHistogramsPerHit() +{ } // --------------------------------------------------------------------------------------------------------------------- @@ -484,45 +80,9 @@ InitStatus CbmCaInputQaTrd::InitDataBranches() // STS tracking detector interface fpDetInterface = CbmTrdTrackingInterface::Instance(); - LOG_IF(fatal, !fpDetInterface) << "\033[1;31m" << fName << ": tracking detector interface is undefined\033[0m"; - - // FairRootManager - auto* pFairRootManager = FairRootManager::Instance(); - LOG_IF(fatal, !pFairRootManager) << "\033[1;31m" << fName << ": FairRootManager instance is a null pointer\033[0m"; - - // Time-slice - fpTimeSlice = dynamic_cast<CbmTimeSlice*>(pFairRootManager->GetObject("TimeSlice.")); - LOG_IF(fatal, !fpTimeSlice) << "\033[1;31m" << fName << ": time-slice branch is not found\033[0m"; - - // ----- Reconstructed data-branches initialization - - // Hits container - fpHits = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject("TrdHit")); - LOG_IF(fatal, !fpHits) << "\033[1;31m" << fName << ": container of reconstructed hits in TRD is not found\033[0m"; - - // ----- MC-information branches initialization - if (IsMCUsed()) { - // MC manager - fpMCDataManager = dynamic_cast<CbmMCDataManager*>(pFairRootManager->GetObject("MCDataManager")); - LOG_IF(fatal, !fpMCDataManager) << "\033[1;31m" << fName << ": MC data manager branch is not found\033[0m"; - - // MC event list - fpMCEventList = dynamic_cast<CbmMCEventList*>(pFairRootManager->GetObject("MCEventList.")); - LOG_IF(fatal, !fpMCEventList) << "\033[1;31m" << fName << ": MC event list branch is not found\033[0m"; - - // MC tracks - fpMCTracks = fpMCDataManager->InitBranch("MCTrack"); - LOG_IF(fatal, !fpMCTracks) << "\033[1;31m" << fName << ": MC track branch is not found\033[0m"; - - // MC points - fpMCPoints = fpMCDataManager->InitBranch("TrdPoint"); - LOG_IF(fatal, !fpMCTracks) << "\033[1;31m" << fName << ": MC point branch is not found for TRD\033[0m"; - - // Hit matches - fpHitMatches = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject("TrdHitMatch")); - LOG_IF(fatal, !fpHitMatches) << "\033[1;31m]" << fName << ": hit match branch is not found for TRD\033[0m"; - } - + auto baseInitStatus = CbmCaInputQaBase::InitDataBranches(); + if (kSUCCESS != baseInitStatus) { return baseInitStatus; } + return kSUCCESS; } @@ -530,164 +90,9 @@ InitStatus CbmCaInputQaTrd::InitDataBranches() // InitStatus CbmCaInputQaTrd::InitCanvases() { - gStyle->SetOptFit(1); - int nSt = fpDetInterface->GetNtrackingStations(); - - // ******************** - // ** Drawing options - - // Contours - constexpr auto contColor = kOrange + 7; - constexpr auto contWidth = 2; // Line width - constexpr auto contStyle = 2; // Line style - constexpr auto contFill = 0; // Fill style - - // ********************************* - // ** Hit occupancies - - // ** Occupancy in XY plane - auto* pc_hit_ypos_vs_xpos = MakeQaObject<CbmQaCanvas>("hit_ypos_vs_xpos", "", 1600, 800); - pc_hit_ypos_vs_xpos->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_hit_ypos_vs_xpos->cd(iSt + 1); - fvph_hit_ypos_vs_xpos[iSt]->DrawCopy("colz", ""); - - - // Square boarders of station, which are used for the field approximation - double stXmin = +fpDetInterface->GetXmax(iSt); - double stXmax = -fpDetInterface->GetXmax(iSt); - double stYmin = -fpDetInterface->GetYmax(iSt); - double stYmax = +fpDetInterface->GetYmax(iSt); - - auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - - // ----- Occupancy projections - auto* pc_hit_xpos = MakeQaObject<CbmQaCanvas>("hit_xpos", "", 1400, 800); - pc_hit_xpos->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_hit_xpos->cd(iSt + 1); - auto* phProj = fvph_hit_ypos_vs_xpos[iSt]->ProjectionY(); - phProj->DrawCopy(); - } - - auto* pc_hit_ypos = MakeQaObject<CbmQaCanvas>("hit_ypos", "", 1400, 800); - pc_hit_ypos->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_hit_ypos->cd(iSt + 1); - auto* phProj = fvph_hit_ypos_vs_xpos[iSt]->ProjectionX(); - phProj->DrawCopy(); - } - - // ************ - // ************ - - if (IsMCUsed()) { - - // ********************** - // ** Point occupancies - - // ** Occupancy in XY plane - auto* pc_point_ypos_vs_xpos = MakeQaObject<CbmQaCanvas>("point_ypos_vs_xpos", "", 1600, 800); - pc_point_ypos_vs_xpos->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_point_ypos_vs_xpos->cd(iSt + 1); - fvph_point_ypos_vs_xpos[iSt]->DrawCopy("colz", ""); - - // Square boarders of station, which are used for the field approximation - double stXmin = +fpDetInterface->GetXmax(iSt); - double stXmax = -fpDetInterface->GetXmax(iSt); - double stYmin = -fpDetInterface->GetYmax(iSt); - double stYmax = +fpDetInterface->GetYmax(iSt); - - auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax); - pBox->SetLineWidth(contWidth); - pBox->SetLineStyle(contStyle); - pBox->SetLineColor(contColor); - pBox->SetFillStyle(contFill); - pBox->Draw("SAME"); - } - - // ************** - // ** Residuals - - // x-coordinate - auto* pc_res_x = MakeQaObject<CbmQaCanvas>("res_x", "Residuals for x coordinate", 1600, 800); - pc_res_x->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_res_x->cd(iSt + 1); - fvph_res_x[iSt]->DrawCopy("", ""); - } - - // y-coordinate - auto* pc_res_y = MakeQaObject<CbmQaCanvas>("res_y", "Residuals for y coordinate", 1600, 800); - pc_res_y->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_res_y->cd(iSt + 1); - fvph_res_y[iSt]->DrawCopy("", ""); - } - - // time - auto* pc_res_t = MakeQaObject<CbmQaCanvas>("res_t", "Residuals for time", 1600, 800); - pc_res_t->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_res_t->cd(iSt + 1); - fvph_res_t[iSt]->DrawCopy("", ""); - } - - - // ********** - // ** Pulls - - // x-coordinate - auto* pc_pull_x = MakeQaObject<CbmQaCanvas>("pull_x", "Pulls for x coordinate", 1600, 800); - pc_pull_x->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_pull_x->cd(iSt + 1); - fvph_pull_x[iSt]->DrawCopy("", ""); - } - - // y-coordinate - auto* pc_pull_y = MakeQaObject<CbmQaCanvas>("pull_y", "Pulls for y coordinate", 1600, 800); - pc_pull_y->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_pull_y->cd(iSt + 1); - fvph_pull_y[iSt]->DrawCopy("", ""); - } - - // time - auto* pc_pull_t = MakeQaObject<CbmQaCanvas>("pull_t", "Pulls for time", 1600, 800); - pc_pull_t->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_pull_t->cd(iSt + 1); - fvph_pull_t[iSt]->DrawCopy("", ""); - } - - - // ************************************ - // ** Hit reconstruction efficiencies - - auto* pc_reco_eff_vs_r = - MakeQaObject<CbmQaCanvas>("reco_eff_vs_r", "Hit efficiencies wrt distance from center", 1600, 800); - pc_reco_eff_vs_r->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_reco_eff_vs_r->cd(iSt + 1); - fvpe_reco_eff_vs_r[iSt]->DrawCopy("AP", ""); - } - - auto* pc_reco_eff_vs_xy = MakeQaObject<CbmQaCanvas>("reco_eff_vs_xy", "Hit efficiencies wrt x and y", 1600, 800); - pc_reco_eff_vs_xy->Divide2D(nSt); - for (int iSt = 0; iSt < nSt; ++iSt) { - pc_reco_eff_vs_xy->cd(iSt + 1); - fvpe_reco_eff_vs_xy[iSt]->DrawCopy("colz", ""); - } - } - + auto baseInitStatus = CbmCaInputQaBase::InitCanvases(); + if (kSUCCESS != baseInitStatus) { return baseInitStatus; } + return kSUCCESS; } @@ -695,168 +100,8 @@ InitStatus CbmCaInputQaTrd::InitCanvases() // InitStatus CbmCaInputQaTrd::InitHistograms() { - int nSt = fpDetInterface->GetNtrackingStations(); - - // ----- Histograms initialization - fvph_hit_ypos_vs_xpos.resize(nSt, nullptr); - fvph_hit_ypos_vs_zpos.resize(nSt, nullptr); - fvph_hit_xpos_vs_zpos.resize(nSt, nullptr); - - fvph_hit_station_delta_z.resize(nSt); - - fvph_hit_dx.resize(nSt, nullptr); - fvph_hit_dy.resize(nSt, nullptr); - fvph_hit_dt.resize(nSt, nullptr); - - if (IsMCUsed()) { - fvph_n_points_per_hit.resize(nSt, nullptr); - fvph_point_ypos_vs_xpos.resize(nSt, nullptr); - fvph_point_xpos_vs_zpos.resize(nSt, nullptr); - fvph_point_ypos_vs_zpos.resize(nSt, nullptr); - fvph_point_hit_delta_z.resize(nSt, nullptr); - fvph_res_x.resize(nSt, nullptr); - fvph_res_y.resize(nSt, nullptr); - fvph_res_t.resize(nSt, nullptr); - fvph_pull_x.resize(nSt, nullptr); - fvph_pull_y.resize(nSt, nullptr); - fvph_pull_t.resize(nSt, nullptr); - fvph_res_x_vs_x.resize(nSt, nullptr); - fvph_res_y_vs_y.resize(nSt, nullptr); - fvph_res_t_vs_t.resize(nSt, nullptr); - fvph_pull_x_vs_x.resize(nSt, nullptr); - fvph_pull_y_vs_y.resize(nSt, nullptr); - fvph_pull_t_vs_t.resize(nSt, nullptr); - fvpe_reco_eff_vs_xy.resize(nSt, nullptr); - fvpe_reco_eff_vs_r.resize(nSt, nullptr); - } - - for (int iSt = 0; iSt < nSt; ++iSt) { - TString nsuff = Form("_st%d", iSt); - TString tsuff = Form(" in TRD station %d", iSt); - TString npref = Form("station%d/", iSt); - TString xyBinningTag = Form("xy_station%d", iSt); - - TString sN = ""; - TString sT = ""; - - double xMax = fpDetInterface->GetXmax(iSt); - double yMax = fpDetInterface->GetYmax(iSt); - double xRes = 0.05; // [cm] - double yRes = 0.20; // [cm] - int nBinsX = static_cast<int>(2 * xMax / xRes); - int nBinsY = static_cast<int>(2 * yMax / yRes); - - // Hit occupancy - sN = npref + "hit_ypos_vs_xpos" + nsuff + ";" + xyBinningTag; - sT = (TString) "Hit occupancy in xy-Plane" + tsuff + ";x_{hit} [cm];y_{hit} [cm]"; - fvph_hit_ypos_vs_xpos[iSt] = MakeQaObject<TH2F>(sN, sT, nBinsX, -xMax, xMax, nBinsY, -yMax, yMax); - - sN = npref + "hit_xpos_vs_zpos" + nsuff; - sT = (TString) "Hit occupancy in xz-Plane" + tsuff + ";z_{hit} [cm];x_{hit} [cm]"; - fvph_hit_xpos_vs_zpos[iSt] = MakeQaObject<TH2F>(sN, sT, fNbinsZ, fvRHitZ[0], fvRHitZ[1], nBinsX, -xMax, xMax); - - sN = npref + "hit_ypos_vs_zpos" + nsuff; - sT = (TString) "Hit occupancy in yz-plane" + tsuff + ";z_{hit} [cm];y_{hit} [cm]"; - fvph_hit_ypos_vs_zpos[iSt] = MakeQaObject<TH2F>(sN, sT, fNbinsZ, fvRHitZ[0], fvRHitZ[1], nBinsY, -yMax, yMax); - - // Hit errors - sN = npref + "hit_dx" + nsuff; - sT = (TString) "Hit position error along x-axis" + tsuff + ";dx_{hit} [cm]"; - fvph_hit_dx[iSt] = MakeQaObject<TH1F>(sN, sT, fNbins, fvRHitDx[0], fvRHitDx[1]); - - sN = npref + "hit_dy" + nsuff; - sT = (TString) "Hit position error along y-axis" + tsuff + ";dy_{hit} [cm]"; - fvph_hit_dy[iSt] = MakeQaObject<TH1F>(sN, sT, fNbins, fvRHitDy[0], fvRHitDy[1]); - - sN = npref + "hit_dt" + nsuff; - sT = (TString) "Hit time error" + tsuff + ";dt_{hit} [ns]"; - fvph_hit_dt[iSt] = MakeQaObject<TH1F>(sN, sT, fNbins, fvRHitDt[0], fvRHitDt[1]); - - sN = npref + "hit_station_delta_z" + nsuff; - sT = (TString) "Different between hit and station z-positions" + tsuff + ";z_{hit} - z_{st} [cm]"; - fvph_hit_station_delta_z[iSt] = MakeQaObject<TH1F>(sN, sT, fNbins, -5., 5); - - // ----- Initialize histograms, which are use MC-information - if (IsMCUsed()) { - sN = npref + "n_points_per_hit" + nsuff; - sT = (TString) "Number of points per hit" + tsuff + ";N_{point}/hit"; - fvph_n_points_per_hit[iSt] = MakeQaObject<TH1F>(sN, sT, 10, -0.5, 9.5); - - // Point occupancy - sN = npref + "point_ypos_vs_xpos" + nsuff + ";" + xyBinningTag; - sT = (TString) "Point occupancy in XY plane" + tsuff + ";x_{MC} [cm];y_{MC} [cm]"; - fvph_point_ypos_vs_xpos[iSt] = MakeQaObject<TH2F>(sN, sT, nBinsX, -xMax, xMax, nBinsY, -yMax, yMax); - - sN = npref + "point_xpos_vs_zpos" + nsuff; - sT = (TString) "Point Occupancy in XZ plane" + tsuff + ";z_{MC} [cm];x_{MC} [cm]"; - fvph_point_xpos_vs_zpos[iSt] = MakeQaObject<TH2F>(sN, sT, fNbinsZ, fvRHitZ[0], fvRHitZ[1], nBinsX, -xMax, xMax); - - sN = npref + "point_ypos_vs_zpos" + nsuff; - sT = (TString) "Point Occupancy in YZ Plane" + tsuff + ";z_{MC} [cm];y_{MC} [cm]"; - fvph_point_ypos_vs_zpos[iSt] = MakeQaObject<TH2F>(sN, sT, fNbinsZ, fvRHitZ[0], fvRHitZ[1], nBinsY, -yMax, yMax); - - // Difference between MC input z and hit z coordinates - sN = npref + "point_hit_delta_z" + nsuff; - sT = (TString) "Distance between point and hit along z axis" + tsuff + ";z_{reco} - z_{MC} [cm]"; - fvph_point_hit_delta_z[iSt] = MakeQaObject<TH1F>(sN, sT, fNbins, -0.04, 0.04); - - sN = npref + "res_x" + nsuff; - sT = (TString) "Residuals for x-coordinate" + tsuff + ";x_{reco} - x_{MC} [cm]"; - fvph_res_x[iSt] = MakeQaObject<TH1F>(sN, sT, fNbins, fvRResX[0], fvRResX[1]); - - sN = npref + "res_y" + nsuff; - sT = (TString) "Residuals for y-coordinate" + tsuff + ";y_{reco} - y_{MC} [cm]"; - fvph_res_y[iSt] = MakeQaObject<TH1F>(sN, sT, fNbins, fvRResY[0], fvRResY[1]); - - sN = npref + "res_t" + nsuff; - sT = (TString) "Residuals for time" + tsuff + ";t_{reco} - t_{MC} [ns]"; - fvph_res_t[iSt] = MakeQaObject<TH1F>(sN, sT, fNbins, fvRResT[0], fvRResT[1]); - - sN = npref + "pull_x" + nsuff; - sT = (TString) "Pulls for x-coordinate" + tsuff + ";(x_{reco} - x_{MC}) / #sigma_{x}^{reco}"; - fvph_pull_x[iSt] = MakeQaObject<TH1F>(sN, sT, fNbins, fvRPullX[0], fvRPullX[1]); - - sN = npref + "pull_y" + nsuff; - sT = (TString) "Pulls for y-coordinate" + tsuff + ";(y_{reco} - y_{MC}) / #sigma_{y}^{reco}"; - fvph_pull_y[iSt] = MakeQaObject<TH1F>(sN, sT, fNbins, fvRPullY[0], fvRPullY[1]); - - sN = npref + "pull_t" + nsuff; - sT = (TString) "Pulls for time" + tsuff + ";(t_{reco} - t_{MC}) / #sigma_{t}^{reco}"; - fvph_pull_t[iSt] = MakeQaObject<TH1F>(sN, sT, fNbins, fvRPullT[0], fvRPullT[1]); - - sN = npref + "res_x_vs_x" + nsuff; - sT = (TString) "Residuals for x-coordinate" + tsuff + ";x_{reco} [cm];x_{reco} - x_{MC} [cm]"; - fvph_res_x_vs_x[iSt] = MakeQaObject<TH2F>(sN, sT, fNbins, 0, 0, fNbins, fvRResX[0], fvRResX[1]); - - sN = npref + "res_y_vs_y" + nsuff; - sT = (TString) "Residuals for y-coordinate" + tsuff + ";y_{reco} [cm];y_{reco} - y_{MC} [cm]"; - fvph_res_y_vs_y[iSt] = MakeQaObject<TH2F>(sN, sT, fNbins, 0, 0, fNbins, fvRResY[0], fvRResY[1]); - - sN = npref + "res_t_vs_t" + nsuff; - sT = (TString) "Residuals for time" + tsuff + "t_{reco} [ns];t_{reco} - t_{MC} [ns]"; - fvph_res_t_vs_t[iSt] = MakeQaObject<TH2F>(sN, sT, fNbins, 0, 0, fNbins, fvRResT[0], fvRResT[1]); - - sN = npref + "pull_x_vs_x" + nsuff; - sT = (TString) "Pulls for x-coordinate" + tsuff + "x_{reco} [cm];(x_{reco} - x_{MC}) / #sigma_{x}^{reco}"; - fvph_pull_x_vs_x[iSt] = MakeQaObject<TH2F>(sN, sT, fNbins, 0, 0, fNbins, fvRPullX[0], fvRPullX[1]); - - sN = npref + "pull_y_vs_y" + nsuff; - sT = (TString) "Pulls for y-coordinate" + tsuff + ";y_{reco} [cm];(y_{reco} - y_{MC}) / #sigma_{y}^{reco}"; - fvph_pull_y_vs_y[iSt] = MakeQaObject<TH2F>(sN, sT, fNbins, 0, 0, fNbins, fvRPullY[0], fvRPullY[1]); - - sN = npref + "pull_t_vs_t" + nsuff; - sT = (TString) "Pulls for time" + tsuff + ";t_{reco} [ns];(t_{reco} - t_{MC}) / #sigma_{t}^{reco}"; - fvph_pull_t_vs_t[iSt] = MakeQaObject<TH2F>(sN, sT, fNbins, 0, 0, fNbins, fvRPullT[0], fvRPullT[1]); - - sN = npref + "reco_eff_vs_xy" + nsuff + ";" + xyBinningTag; - sT = (TString) "Hit rec. efficiency" + tsuff + ";x_{MC} [cm];y_{MC} [cm];#epsilon"; - fvpe_reco_eff_vs_xy[iSt] = - MakeQaObject<TProfile2D>(sN, sT, fNbins, fvRHitX[0], fvRHitX[1], fNbins, fvRHitY[0], fvRHitY[1]); + auto baseInitStatus = CbmCaInputQaBase::InitHistograms(); + if (kSUCCESS != baseInitStatus) { return baseInitStatus; } - sN = npref + "reco_eff_vs_r" + nsuff; - sT = (TString) "Hit rec. efficiency" + tsuff + ";r_{MC} [cm];#epsilon"; - fvpe_reco_eff_vs_r[iSt] = MakeQaObject<TProfile>(sN, sT, 100, 0, 100, 0., 1.); - } - } return kSUCCESS; } diff --git a/reco/L1/qa/CbmCaInputQaTrd.h b/reco/L1/qa/CbmCaInputQaTrd.h index 70bd9e7f481ac400d18a8c480d349d028ba8e4fb..54b2303368ff3ba85311e7d9856d17d8c61e4f4b 100644 --- a/reco/L1/qa/CbmCaInputQaTrd.h +++ b/reco/L1/qa/CbmCaInputQaTrd.h @@ -12,155 +12,48 @@ #define CbmCaInputQaTrd_h 1 #include "CbmCaInputQaBase.h" -#include "CbmMCDataManager.h" -#include "CbmQaTask.h" - -#include "TMath.h" - -#include <set> -#include <unordered_map> -#include <vector> - -class CbmMCEventList; -class CbmMCDataArray; -class CbmMCDataManager; -class CbmTimeSlice; -class CbmMatch; -class CbmTrdHit; -class CbmTrdTrackingInterface; -class TClonesArray; -class TH1F; -class TH2F; -class TProfile; -class TProfile2D; - -/// A QA-task class, which provides assurance of TRD hits and geometry -class CbmCaInputQaTrd : public CbmQaTask, public CbmCaInputQaBase { + +/// A QA-task class, which provides assurance of TOF hits and geometry +/// +class CbmCaInputQaTrd : public CbmCaInputQaBase<L1DetectorID::kTrd> { public: - /// Constructor from parameters - /// \param verbose Verbose level - /// \param isMCUsed Flag, whether MC information is available for this task + /// @brief Constructor from parameters + /// @param verbose Verbose level + /// @param isMCUsed Flag, whether MC information is available for this task CbmCaInputQaTrd(int verbose, bool isMCUsed); - ClassDef(CbmCaInputQaTrd, 0); - protected: - // ******************************************** - // ** Virtual method override from CbmQaTask ** - // ******************************************** - - /// Checks results of the QA - /// \return Success flag + /// @brief Checks results of the QA + /// @return Success flag bool Check(); - /// Initializes data branches + /// @brief Initializes data branches InitStatus InitDataBranches(); - /// Initializes canvases + /// @brief Initializes canvases InitStatus InitCanvases(); - /// Initializes histograms + /// @brief Initializes histograms InitStatus InitHistograms(); - /// Fills histograms per event or time-slice - void FillHistograms(); + /// @brief Fills histograms per event or time-slice + void FillHistograms() { CbmCaInputQaBase::FillHistograms(); } - /// De-initializes histograms + /// @brief De-initializes histograms void DeInit(); -private: - // ********************* - // ** Private methods ** - // ********************* - - // ********************* - // ** Class variables ** - // ********************* - - - // ** Data branches ** - CbmTrdTrackingInterface* fpDetInterface = nullptr; ///< Instance of detector interface - - CbmTimeSlice* fpTimeSlice = nullptr; ///< Pointer to current time-slice - - TClonesArray* fpHits = nullptr; ///< Array of hits - - CbmMCDataManager* fpMCDataManager = nullptr; ///< MC data manager - CbmMCEventList* fpMCEventList = nullptr; ///< MC event list - CbmMCDataArray* fpMCTracks = nullptr; ///< Array of MC tracks - CbmMCDataArray* fpMCPoints = nullptr; ///< Array of MC points - - TClonesArray* fpHitMatches = nullptr; ///< Array of hit matches - - std::vector<char> fvTrdStationType; ///< Station type: TRD-1D and TRD-2D - - // ** Histograms ** - double fEffRangeMin = 35.; ///< Lower limit of hit distance for efficiency integration [cm] - double fEffRangeMax = 70.; ///< Upper limit of hit distance for efficiency integration [cm] - - std::array<double, 2> fvRHitX = {-100., 100}; ///< Range for hit x coordinate [cm] - std::array<double, 2> fvRHitY = {-100., 100}; ///< Range for hit y coordinate [cm] - std::array<double, 2> fvRHitZ = {30., 300.}; ///< Range for hit z coordinate [cm] + /// @brief Defines parameters of the task + void DefineParameters(); - int fNbins = 100; ///< General number of bins - int fNbinsZ = 800; ///< Number of bins for z coordinate axis + /// @brief Fills histograms per hit + void FillHistogramsPerHit(); + + /// @brief Fills histograms per MC point + void FillHistogramsPerPoint() {} - std::array<double, 2> fvRHitDx = {-.05, .005}; ///< Range for hit x coordinate [cm] - std::array<double, 2> fvRHitDy = {-.05, .005}; ///< Range for hit y coordinate [cm] - std::array<double, 2> fvRHitDt = {-10., 10.}; ///< Range for hit time [ns] - - std::array<double, 2> fvRResX = {-2., 2.}; ///< Range for residual of x coordinate [cm] - std::array<double, 2> fvRResY = {-2., 2.}; ///< Range for residual of y coordinate [cm] - std::array<double, 2> fvRResT = {-20., 20.}; ///< Range for residual of time [ns] - - std::array<double, 2> fvRPullX = {-10., 10.}; ///< Range for pull of x coordinate - std::array<double, 2> fvRPullY = {-10., 10.}; ///< Range for pull of y coordinate - std::array<double, 2> fvRPullT = {-5., 5.}; ///< Range for pull of time - - // hit occupancy - std::vector<TH2F*> fvph_hit_ypos_vs_xpos; ///< hit ypos vs xpos in different stations - std::vector<TH2F*> fvph_hit_xpos_vs_zpos; ///< hit xpos vs zpos in different stations - std::vector<TH2F*> fvph_hit_ypos_vs_zpos; ///< hit ypos vs zpos in different stations - - // difference between hit and station z positions - std::vector<TH1F*> fvph_hit_station_delta_z; ///< Difference between station and hit z positions [cm] - - // hit errors - std::vector<TH1F*> fvph_hit_dx; - std::vector<TH1F*> fvph_hit_dy; - std::vector<TH1F*> fvph_hit_dt; - - // MC points occupancy - std::vector<TH1F*> fvph_n_points_per_hit; ///< number of points per one hit in different stations - - - std::vector<TH2F*> fvph_point_ypos_vs_xpos; ///< point ypos vs xpos in different stations - std::vector<TH2F*> fvph_point_xpos_vs_zpos; ///< point xpos vs zpos in different stations - std::vector<TH2F*> fvph_point_ypos_vs_zpos; ///< point ypos vs zpos in different stations - - std::vector<TH1F*> fvph_point_hit_delta_z; ///< difference between z of hit and MC point in different stations - - // Residuals - std::vector<TH1F*> fvph_res_x; ///< residual for x coordinate in different stations - std::vector<TH1F*> fvph_res_y; ///< residual for y coordinate in different stations - std::vector<TH1F*> fvph_res_t; ///< residual for t coordinate in different stations - - std::vector<TH2F*> fvph_res_x_vs_x; ///< residual for x coordinate in different stations - std::vector<TH2F*> fvph_res_y_vs_y; ///< residual for y coordinate in different stations - std::vector<TH2F*> fvph_res_t_vs_t; ///< residual for t coordinate in different stations - - // Pulls - std::vector<TH1F*> fvph_pull_x; ///< pull for x coordinate in different stations - std::vector<TH1F*> fvph_pull_y; ///< pull for y coordinate in different stations - std::vector<TH1F*> fvph_pull_t; ///< pull for t coordinate in different stations - - std::vector<TH2F*> fvph_pull_x_vs_x; ///< pull for x coordinate in different stations - std::vector<TH2F*> fvph_pull_y_vs_y; ///< pull for y coordinate in different stations - std::vector<TH2F*> fvph_pull_t_vs_t; ///< pull for t coordinate in different stations +private: - // Hit efficiencies - std::vector<TProfile2D*> fvpe_reco_eff_vs_xy; ///< Efficiency of hit reconstruction vs. x and y coordinates of MC - std::vector<TProfile*> fvpe_reco_eff_vs_r; ///< Efficiency of hit reconstruction vs. distance from center + ClassDef(CbmCaInputQaTrd, 0); }; #endif // CbmCaInputQaTrd_h