diff --git a/core/base/CbmTrackingDetectorInterfaceBase.h b/core/base/CbmTrackingDetectorInterfaceBase.h index ceb1d6a88a6341bf9acdd5ccf32d5030471ab6b4..dc9580829e2d38d7454a12b82fa26969c47a5c89 100644 --- a/core/base/CbmTrackingDetectorInterfaceBase.h +++ b/core/base/CbmTrackingDetectorInterfaceBase.h @@ -12,7 +12,9 @@ #ifndef CbmTrackingDetectorInterfaceBase_h #define CbmTrackingDetectorInterfaceBase_h 1 -/// Abstract class, which should be inherited by every detecting subsystem interface class +class CbmPixelHit; + +/// Abstract class, which should be inherited by every detecting subsystem tracking interface class /// class CbmTrackingDetectorInterfaceBase { public: @@ -22,25 +24,10 @@ public: /// Gets actual number of stations, provided by the current geometry setup virtual int GetNtrackingStations() const = 0; - /// Gets time resolution for a station - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Time resolution [ns] - virtual double GetTimeResolution(int stationId) const = 0; - - /// Gets z component of the station position - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Z position of station [cm] - virtual double GetZ(int stationId) const = 0; - - /// Gets max size 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] - virtual double GetXmax(int stationId) const = 0; - - /// Gets max size of a station along the Y-axis + /// Gets station radiation length /// \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] - virtual double GetYmax(int stationId) const = 0; + /// \return Radiation length [cm] + virtual double GetRadLength(int stationId) const = 0; /// Gets size of inner radius of station /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) @@ -52,16 +39,6 @@ public: /// \return Size of station outer radius [cm] virtual double GetRmax(int stationId) const = 0; - /// Gets station thickness along the Z-axis - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Station thickness [cm] - virtual double GetThickness(int stationId) const = 0; - - /// Gets station radiation length - /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) - /// \return Radiation length [cm] - virtual double GetRadLength(int stationId) const = 0; - /// Gets front strips stereo angle /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) /// \return Absolute stereo angle for front strips [rad] @@ -82,10 +59,46 @@ public: /// \return Spatial resolution (RMS) for front strips [cm] virtual double GetStripsSpatialRmsBack(int stationId) const = 0; + /// Gets station thickness along the Z-axis + /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// \return Station thickness [cm] + virtual double GetThickness(int stationId) const = 0; + + /// Gets time resolution for a station + /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// \return Time resolution [ns] + virtual double GetTimeResolution(int stationId) const = 0; + + /// Gets a tracking station of a CbmPixelHit + /// \param hit A pointer to CbmPixelHit + /// \return Local index of the tracking station + virtual int GetTrackingStationIndex(const CbmPixelHit* hit) const = 0; + + /// Gets max size 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] + virtual double GetXmax(int stationId) const = 0; + + /// 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] + virtual double GetYmax(int stationId) const = 0; + + /// Gets z component of the station position + /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1]) + /// \return Z position of station [cm] + virtual double GetZ(int stationId) const = 0; + /// Check if station provides time measurements /// \param stationId Tracking station ID in the setup /// \return Flag: true - station provides time measurements, false - station does not provide time measurements virtual bool IsTimeInfoProvided(int stationId) const = 0; + + /// Technical flag: true - all downcasts are done with dynamic_cast followed by checks for nullptr (increases + // computing time, better for debug), false - all downcasts are done with static_cast without sanity checks + // (decreases computting time, but can cause bugs) + + static constexpr bool kUseDynamicCast {true}; }; #endif // CbmTrackingDetectorInterfaceBase_h diff --git a/core/detectors/much/CbmMuchTrackingInterface.h b/core/detectors/much/CbmMuchTrackingInterface.h index ea49fc092fa02e0e8f1b608266fb393a3ed08d1c..0b637738a6ca6f8f63749c3bbd6815b2747588e3 100644 --- a/core/detectors/much/CbmMuchTrackingInterface.h +++ b/core/detectors/much/CbmMuchTrackingInterface.h @@ -16,6 +16,7 @@ #include "CbmMuchModuleGem.h" #include "CbmMuchPad.h" #include "CbmMuchStation.h" +#include "CbmPixelHit.h" #include "CbmTrackingDetectorInterfaceBase.h" #include "FairTask.h" @@ -98,6 +99,15 @@ public: /// \return Time resolution [ns] double GetTimeResolution(int /*stationId*/) const { return 3.9; } + /// Gets a tracking station of a CbmPixelHit + /// \param hit A pointer to CbmPixelHit + /// \return Local index of the tracking station + int GetTrackingStationIndex(const CbmPixelHit* hit) const + { + return CbmMuchGeoScheme::GetStationIndex(hit->GetAddress()) * 3 + + CbmMuchGeoScheme::GetLayerIndex(hit->GetAddress()); + } + /// Gets max size 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] diff --git a/core/detectors/sts/CbmStsTrackingInterface.h b/core/detectors/sts/CbmStsTrackingInterface.h index a105e1f8d37c86343cc200335610ee2956558e28..d3ffb736194d759782c5bf96ead509f84b7548e7 100644 --- a/core/detectors/sts/CbmStsTrackingInterface.h +++ b/core/detectors/sts/CbmStsTrackingInterface.h @@ -12,6 +12,7 @@ #ifndef CbmStsTrackingInterface_h #define CbmStsTrackingInterface_h 1 +#include "CbmPixelHit.h" #include "CbmStsParSetModule.h" #include "CbmStsParSetSensor.h" #include "CbmStsParSetSensorCond.h" @@ -103,6 +104,14 @@ public: /// \return Time resolution [ns] double GetTimeResolution(int /*stationId*/) const { return 5.; } + /// Gets a tracking station of a CbmPixelHit + /// \param hit A pointer to CbmPixelHit + /// \return Local index of the tracking station + int GetTrackingStationIndex(const CbmPixelHit* hit) const + { + return CbmStsSetup::Instance()->GetStationNumber(hit->GetAddress()); + } + /// Gets max size 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] diff --git a/core/detectors/tof/CbmTofTrackingInterface.h b/core/detectors/tof/CbmTofTrackingInterface.h index 8a7e4eb480c390aa95a875430ad44d5b47946adb..2c3f2e222a6efdda3f5492e01d669546606fcf83 100644 --- a/core/detectors/tof/CbmTofTrackingInterface.h +++ b/core/detectors/tof/CbmTofTrackingInterface.h @@ -12,6 +12,8 @@ #ifndef CbmTofTrackingInterface_h #define CbmTofTrackingInterface_h 1 +#include "CbmPixelHit.h" +#include "CbmTofAddress.h" #include "CbmTofDigiBdfPar.h" #include "CbmTofHit.h" #include "CbmTrackingDetectorInterfaceBase.h" @@ -94,12 +96,16 @@ public: double GetTimeResolution(int /*stationId*/) const { return 0.075; } /// Gets a tracking station of a ToF hit - /// \param hit A pointer to ToF hit - __attribute__((always_inline)) int GetTrackingStation(CbmTofHit* hit) const + /// \param hit A pointer to ToF hit + /// \return Local index of the tracking station + int GetTrackingStationIndex(const CbmPixelHit* hit) const { - int iSt = fDigiBdfPar->GetTrackingStation(hit); + int iSt = fDigiBdfPar->GetTrackingStation(CbmTofAddress::GetSmType(hit->GetAddress()), + CbmTofAddress::GetSmId(hit->GetAddress()), + CbmTofAddress::GetRpcId(hit->GetAddress())); + if (fIfMissingHits) { - /// Recalculate station index for inconsistent hits + // Recalculate station index for inconsistent hits if (hit->GetX() > 20. && hit->GetZ() > 270. && 1 == iSt) { iSt = 2; } } return iSt; diff --git a/core/detectors/trd/CbmTrdTrackingInterface.h b/core/detectors/trd/CbmTrdTrackingInterface.h index 1291540a1151482583122cd1979581c667d5d186..f193db43381efbce19166fc2eb6aa831776cf59a 100644 --- a/core/detectors/trd/CbmTrdTrackingInterface.h +++ b/core/detectors/trd/CbmTrdTrackingInterface.h @@ -12,7 +12,9 @@ #ifndef CbmTrdTrackingInterface_h #define CbmTrdTrackingInterface_h 1 +#include "CbmPixelHit.h" #include "CbmTrackingDetectorInterfaceBase.h" +#include "CbmTrdAddress.h" #include "CbmTrdParModDigi.h" #include "CbmTrdParSetDigi.h" @@ -93,6 +95,11 @@ public: /// \return Time resolution [ns] double GetTimeResolution(int /*stationId*/) const { return 10.; } + /// Gets a tracking station of a CbmPixelHit + /// \param hit A pointer to CbmPixelHit + /// \return Local index of the tracking station + int GetTrackingStationIndex(const CbmPixelHit* hit) const { return CbmTrdAddress::GetLayerId(hit->GetAddress()); } + /// Gets max size 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] diff --git a/macro/run/run_qa.C b/macro/run/run_qa.C index f462df4b65dc9d11bd6617fbfda91b6c4c25fd8d..17e1d89a247e49c78c22860e148427a06d60c1f7 100644 --- a/macro/run/run_qa.C +++ b/macro/run/run_qa.C @@ -160,6 +160,8 @@ void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "d run->AddTask(mcManager); // ------------------------------------------------------------------------ + run->AddTask(new CbmTrackingDetectorInterfaceInit()); // Geometry interface initializer for tracker + // ----- Match reco to MC ------ CbmMatchRecoToMC* matchTask = new CbmMatchRecoToMC(); run->AddTask(matchTask); @@ -189,6 +191,7 @@ void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "d if (CbmSetup::Instance()->IsActive(ECbmModuleId::kSts)) { //run->AddTask(new CbmStsDigitizeQa()); //opens lots of windows run->AddTask(new CbmStsFindTracksQa()); + run->AddTask(new CbmTrackingInputQaSts()); } // ------------------------------------------------------------------------ diff --git a/mvd/CbmMvdTrackingInterface.h b/mvd/CbmMvdTrackingInterface.h index 6868e1a0396bd51ee4511f90142a9850fcced63f..b2df5b841841a6bd636ca1faab47ac7896748949 100644 --- a/mvd/CbmMvdTrackingInterface.h +++ b/mvd/CbmMvdTrackingInterface.h @@ -13,7 +13,9 @@ #define CbmMvdTrackingInterface_h 1 #include "CbmMvdDetector.h" +#include "CbmMvdHit.h" #include "CbmMvdStationPar.h" +#include "CbmPixelHit.h" #include "CbmTrackingDetectorInterfaceBase.h" #include "FairTask.h" @@ -100,6 +102,20 @@ public: /// \return Time resolution [ns] double GetTimeResolution(int /*stationId*/) const { return 1000.; } + /// Gets a tracking station of a CbmPixelHit + /// \param hit A pointer to CbmPixelHit + /// \return Local index of the tracking station + int GetTrackingStationIndex(const CbmPixelHit* hit) const + { + auto hitMvd = [&] { + if constexpr (kUseDynamicCast) { return dynamic_cast<const CbmMvdHit*>(hit); } + else { + return static_cast<const CbmMvdHit*>(hit); + } + }(); + return hitMvd->GetStationNr(); + } + /// Gets max size of a tracking 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] diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt index cc8a444983f5f144b9a11bd372c6e35436b08421..f21620c1f06bfafd718ee3afad74dea22957eb46 100644 --- a/reco/L1/CMakeLists.txt +++ b/reco/L1/CMakeLists.txt @@ -157,6 +157,8 @@ ParticleFinder/CbmL1PFFitter.cxx ParticleFinder/CbmL1PFMCParticle.cxx qa/CbmTrackerInputQaTrd.cxx +qa/CbmTrackingInputQaBase.cxx +qa/CbmTrackingInputQaSts.cxx ) set(HEADERS @@ -189,6 +191,8 @@ OffLineInterface/CbmL1GlobalFindTracksEvents.h L1Algo/L1Def.h L1Algo/L1Vector.h qa/CbmTrackerInputQaTrd.h +qa/CbmTrackingInputQaBase.h +qa/CbmTrackingInputQaSts.h ) diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx index efc2fa5dc693ba62f72bcbc97a7e9ae4d24dc372..e74aedc781882285781511a73a800fa1ec223be0 100644 --- a/reco/L1/CbmL1ReadEvent.cxx +++ b/reco/L1/CbmL1ReadEvent.cxx @@ -1014,9 +1014,9 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength, th.id = tmpHits.size(); - if (0x00202806 == mh->GetAddress() || 0x00002806 == mh->GetAddress()) continue; + if (0x00202806 == mh->GetAddress() || 0x00002806 == mh->GetAddress()) continue; // TODO: Why? (S.Zharko) - int sttof = CbmTofTrackingInterface::Instance()->GetTrackingStation(mh); + int sttof = CbmTofTrackingInterface::Instance()->GetTrackingStationIndex(mh); if (sttof < 0) continue; diff --git a/reco/L1/L1LinkDef.h b/reco/L1/L1LinkDef.h index aa220f0e37adf8e97793ad0b8fcfa29813bd9769..a4407477eccbcf4b506afd09aae87ad39f88d318 100644 --- a/reco/L1/L1LinkDef.h +++ b/reco/L1/L1LinkDef.h @@ -29,6 +29,9 @@ //#pragma link C++ class CbmL1SttHit+; //#pragma link C++ class CbmL1SttTrackFinder+; //#pragma link C++ class CbmL1SttTrack+; +#pragma link C++ class CbmTrackingInputQaBase + ; +#pragma link C++ class CbmTrackingInputQaSts + ; #pragma link C++ class CbmTrackerInputQaTrd + ; + #endif diff --git a/reco/L1/qa/CbmTrackingInputQaBase.cxx b/reco/L1/qa/CbmTrackingInputQaBase.cxx new file mode 100644 index 0000000000000000000000000000000000000000..d0ffaf24753298f0d1a21fe77de763bcb09c0350 --- /dev/null +++ b/reco/L1/qa/CbmTrackingInputQaBase.cxx @@ -0,0 +1,319 @@ +/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergey Gorbunov, Sergei Zharko [committer] */ + +/*************************************************************************************************** + * @file CbmTrackingInputQaBase.cxx + * @brief Base class for the tracking input QA (definition) + * @since 24.06.2022 + * @author S.Zharko <s.zharko@gsi.de> + ***************************************************************************************************/ + +#include "CbmTrackingInputQaBase.h" + +#include "CbmCluster.h" +#include "CbmDigiManager.h" +#include "CbmMCEventList.h" +#include "CbmMatch.h" +#include "CbmPixelHit.h" +#include "CbmStsAddress.h" // TODO: TMP for tests, must be removed!!!! (S.Zharko) +#include "CbmTimeSlice.h" + +#include "FairLogger.h" + +#include "TClonesArray.h" + +#include <iostream> + +// -------------------------------------------------------------------------------------------------------------------- +// +CbmTrackingInputQaBase::CbmTrackingInputQaBase(const char* name, int verbose) : FairTask(name, verbose) {} + +// -------------------------------------------------------------------------------------------------------------------- +// +CbmTrackingInputQaBase::~CbmTrackingInputQaBase() {} + +ClassImp(CbmTrackingInputQaBase); + +// +// -------------------------------------------------------------------------------------------------------------------- +// +void CbmTrackingInputQaBase::Clear(Option_t*) +{ + std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::Clear(Option_t*)\n"; +} +// +// -------------------------------------------------------------------------------------------------------------------- +// +void CbmTrackingInputQaBase::Exec(Option_t*) +{ + std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::Exec(Option_t*)\n"; + + // Run resolution Qa + ResolutionQa(); + + // Update number of events + fNevents.SetVal(fNevents.GetVal() + 1); +} +// +// -------------------------------------------------------------------------------------------------------------------- +// +void CbmTrackingInputQaBase::Finish() { std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::Finish()\n"; } +// +// -------------------------------------------------------------------------------------------------------------------- +// +InitStatus CbmTrackingInputQaBase::Init() +{ + std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::Init()\n"; + + // If detector is not present, fpDetectorInterface is nullptr + if (!fpDetectorInterface) { LOG(fatal) << "CbmTrackingInputQa: Tracking detector interface is undefined"; } + + // Input data initialization + InitInputData(); + + // Base histograms initialization + InitBaseHistograms(); + + // Register histograms in output file + // TODO: Add subfolders for each subsystem like + fpOutFolder = + std::make_shared<TFolder>(TString("TrackingInputQa_") + fDetName, TString("TrackingInputQa_") + fDetName); + fpOutFolder->SetOwner(false); + fpOutFolderHists = fpOutFolder->AddFolder("rawHist", "Raw histograms"); + gStyle->SetOptStat(0); + + fNevents.SetVal(0); + fpOutFolderHists->Add(&fNevents); + + // Register histograms in the output folder + for (auto* histo : fHistList) { + fpOutFolderHists->Add(histo); + } + + // Geometry QA + GeometryQa(); + + return kSUCCESS; +} + +// -------------------------------------------------------------------------------------------------------------------- +// +InitStatus CbmTrackingInputQaBase::InitInputData() +{ + std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::InitInputData()\n"; + + auto* manager = FairRootManager::Instance(); + fpDigiManager = CbmDigiManager::Instance(); + fpDigiManager->Init(); // NOTE: is initialized only once + + // Reconstructed data initialization + fpTimeSlice = dynamic_cast<CbmTimeSlice*>(manager->GetObject("TimeSlice.")); + if (!fpTimeSlice) { + LOG(error) << "Time slice was not found"; + return kERROR; + } + + fpHits = GetHitsInput(); + if (!fpHits) { + LOG(error) << "Hits input array was not found for " << fDetTitle; + return kERROR; + } + + fpClusters = GetClustersInput(); + if (!fpClusters) { + LOG(error) << "Cluster input array was not found for " << fDetTitle; + return kERROR; + } + + // MC data initialization + fpMcManager = dynamic_cast<CbmMCDataManager*>(manager->GetObject("MCDataManager")); + fIsMcPresent = bool(fpMcManager); + + if (fIsMcPresent) { + fpMcEventList = dynamic_cast<CbmMCEventList*>(manager->GetObject("MCEventList.")); + fpMcTracks = fpMcManager->InitBranch("MCTrack"); + fpMcPoints = GetMcPointsInput(); + fpHitMatches = GetHitMatchesInput(); + + if (!fpMcEventList) { + LOG(error) << "MCEventList data were not found"; + return kERROR; + } + + if (!fpMcTracks) { + LOG(error) << "MC tracks data were not found"; + return kERROR; + } + + if (!fpMcPoints) { + LOG(error) << "MC points data were not found for " << fDetTitle; + return kERROR; + } + + if (!fpHitMatches) { + LOG(error) << "Hits match data were not found for " << fDetTitle; + return kERROR; + } + + if (!fpDigiManager->IsMatchPresent(fDetId)) { + LOG(error) << "No digi matches were found for " << fDetTitle; + return kERROR; + } + } + + return kSUCCESS; +} + +// -------------------------------------------------------------------------------------------------------------------- +// +InitStatus CbmTrackingInputQaBase::InitBaseHistograms() +{ + std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::InitBaseHistograms()\n"; + + fHistResidualX.clear(); + fHistResidualY.clear(); + fHistResidualT.clear(); + + fHistPullX.clear(); + fHistPullY.clear(); + fHistPullT.clear(); + + fHistPointsPerHit.clear(); + fHistHitsPerPoint.clear(); + + int nTrackingStations = fpDetectorInterface->GetNtrackingStations(); + + fHistResidualX.reserve(nTrackingStations); + fHistResidualY.reserve(nTrackingStations); + fHistResidualT.reserve(nTrackingStations); + + fHistPullX.reserve(nTrackingStations); + fHistPullY.reserve(nTrackingStations); + fHistPullT.reserve(nTrackingStations); + + fHistPointsPerHit.reserve(nTrackingStations); + fHistHitsPerPoint.reserve(nTrackingStations); + + // Create histograms + TString histName = ""; + TString histTitle = ""; + int nBins = -1; + int rMin = 0.; + int rMax = 0.; + for (int iSt = 0; iSt < nTrackingStations; ++iSt) { + // Residuals + histName = fDetName + Form("_st%d_h1ResidualX", iSt); + histTitle = fDetTitle + Form(" Station %d: Residual X; (X_{reco} - X_{MC}) [cm]", iSt); + std::tie(nBins, rMin, rMax) = fRangeResidualX; + fHistResidualX.emplace_back(histName, histTitle, nBins, rMin, rMax); + + histName = fDetName + Form("_st%d_h1ResidualY", iSt); + histTitle = fDetTitle + Form(" Station %d: Residual Y; (Y_{reco} - Y_{MC}) [cm]", iSt); + std::tie(nBins, rMin, rMax) = fRangeResidualY; + fHistResidualY.emplace_back(histName, histTitle, nBins, rMin, rMax); + + histName = fDetName + Form("_st%d_h1ResidualT", iSt); + histTitle = fDetTitle + Form(" Station %d: Residual T; (T_{reco} - T_{MC}) [cm]", iSt); + std::tie(nBins, rMin, rMax) = fRangeResidualT; + fHistResidualT.emplace_back(histName, histTitle, nBins, rMin, rMax); + + // Pulls + histName = fDetName + Form("_st%d_h1PullX", iSt); + histTitle = fDetTitle + Form(" Station %d: Pull X; (X_{reco} - X_{MC}) / #sigmaX_{reco}", iSt); + std::tie(nBins, rMin, rMax) = fRangePullX; + fHistPullX.emplace_back(histName, histTitle, nBins, rMin, rMax); + + histName = fDetName + Form("_st%d_h1PullY", iSt); + histTitle = fDetTitle + Form(" Station %d: Pull Y; (Y_{reco} - Y_{MC}) / #sigmaY_{reco}", iSt); + std::tie(nBins, rMin, rMax) = fRangePullY; + fHistPullY.emplace_back(histName, histTitle, nBins, rMin, rMax); + + histName = fDetName + Form("_st%d_h1PullT", iSt); + histTitle = fDetTitle + Form(" Station %d: Pull T; (T_{reco} - T_{MC}) / #sigmaT_{reco}", iSt); + std::tie(nBins, rMin, rMax) = fRangePullT; + fHistPullT.emplace_back(histName, histTitle, nBins, rMin, rMax); + + // MC points per one hit + histName = fDetName + Form("_st%d_h1PointsPerHit", iSt); + histTitle = fDetTitle + Form(" Station %d: MC Points per Hit; N_{MC Points} per hit", iSt); + std::tie(nBins, rMin, rMax) = fRangePointsPerHit; + fHistPointsPerHit.emplace_back(histName, histTitle, nBins, rMin, rMax); + + // Hits per one mc point + histName = fDetName + Form("_st%d_h1HitsPerPoint", iSt); + histTitle = fDetTitle + Form(" Station %d: Hit per MC points; N_{hits} per MC point", iSt); + std::tie(nBins, rMin, rMax) = fRangeHitsPerPoint; + fHistHitsPerPoint.emplace_back(histName, histTitle, nBins, rMin, rMax); + } + + // Register histograms + for (int iSt = 0; iSt < nTrackingStations; ++iSt) { + RegisterHist(&fHistResidualX[iSt]); + RegisterHist(&fHistResidualY[iSt]); + RegisterHist(&fHistResidualT[iSt]); + + RegisterHist(&fHistPullX[iSt]); + RegisterHist(&fHistPullY[iSt]); + RegisterHist(&fHistPullT[iSt]); + + RegisterHist(&fHistPointsPerHit[iSt]); + RegisterHist(&fHistHitsPerPoint[iSt]); + } + + return kSUCCESS; +} + +// -------------------------------------------------------------------------------------------------------------------- +// +void CbmTrackingInputQaBase::Reset() { std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::Reset()\n"; } + +// -------------------------------------------------------------------------------------------------------------------- +// +void CbmTrackingInputQaBase::SetParContainers() +{ + std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::SetParContainers()\n"; +} + +// -------------------------------------------------------------------------------------------------------------------- +// +InitStatus CbmTrackingInputQaBase::GeometryQa() +{ + std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::GeometryQa()\n"; + return kSUCCESS; +} + +// -------------------------------------------------------------------------------------------------------------------- +// +void CbmTrackingInputQaBase::ResolutionQa() +{ + std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::ResolutionQa()\n"; + + // Resolution QA is impossible without MC information + if (!fIsMcPresent) { return; } + + int nHits = fpHits->GetEntriesFast(); // Hits in the event/time slice + int nClusters = fpClusters->GetEntriesFast(); // Clusters in the event/time slice + int nDigis = fpDigiManager->GetNofDigis(fDetId); // Number of digies + int nStations = fpDetectorInterface->GetNtrackingStations(); // Number of tracking stations + + // ********************************************* + // ** Loop over hits in this event/time slice ** + // ********************************************* + + for (int iHit = 0; iHit < nHits; ++iHit) { + /// TODO: put it into a function from iHit + + auto* hit = dynamic_cast<CbmPixelHit*>(fpHits->At(iHit)); + if (!hit) { LOG(error) << "There is an empty hit at position " << iHit << " (" << fDetTitle << ")"; } + + // Tracking station index + int iStation = fpDetectorInterface->GetTrackingStationIndex(hit); + if (iStation < 0 || iStation >= nStations) { + LOG(fatal) << "Wrong tracking station index was calculated by the tracking detector interface for " << fDetTitle + << ": station index = " << iStation << " (expected in the range [0, " << nStations - 1 << "])"; + return; + } + + } // Loop over hits in this event/time slice: END +} diff --git a/reco/L1/qa/CbmTrackingInputQaBase.h b/reco/L1/qa/CbmTrackingInputQaBase.h new file mode 100644 index 0000000000000000000000000000000000000000..76feaf7ba904b90471fe3eeef942be6e82e36803 --- /dev/null +++ b/reco/L1/qa/CbmTrackingInputQaBase.h @@ -0,0 +1,195 @@ +/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergey Gorbunov, Sergei Zharko [committer] */ + +/*************************************************************************************************** + * @file CbmTrackingInputQaBase.h + * @brief Base class for the tracking input QA (declaration) + * @since 24.06.2022 + * @author S.Zharko <s.zharko@gsi.de> + ***************************************************************************************************/ + +#ifndef CbmTrackingInputQaBase_h +#define CbmTrackingInputQaBase_h 1 + +#include "CbmMCDataManager.h" +#include "CbmQaCanvas.h" +#include "CbmQaHist.h" +#include "CbmSetup.h" +#include "CbmTrackingDetectorInterfaceBase.h" // Communicate via tracking detector interface + +#include <FairTask.h> + +#include "TFolder.h" +#include "TH1D.h" +#include "TH1I.h" +#include "TH2D.h" +#include "TParameter.h" +#include "TProfile.h" +#include "TString.h" +#include <tuple> + +#include <memory> +#include <vector> + +class CbmDigiManager; +class CbmMCEventList; +class CbmMCDataArray; +class CbmTimeSlice; +class TClonesArray; +class CbmCluster; +class CbmMatch; + +/// Class CbmTrackingInputQaBase is an abstract FairTask inherited class, providing interface for a +/// tracking detector QA including input geometry and data checks. One should inherit this class +/// for a particular detector subsystem. +/// +class CbmTrackingInputQaBase : public FairTask { +public: + /// Constructor from the task name and verbosity level + CbmTrackingInputQaBase(const char* name = "CbmTrackingInputQaBase", int verbose = 1); + + /// Destructor + virtual ~CbmTrackingInputQaBase(); + + /// Suppress copy and move + CbmTrackingInputQaBase(const CbmTrackingInputQaBase&) = delete; + CbmTrackingInputQaBase(CbmTrackingInputQaBase&&) = delete; + CbmTrackingInputQaBase& operator=(const CbmTrackingInputQaBase&) = delete; + CbmTrackingInputQaBase& operator=(CbmTrackingInputQaBase&&) = delete; + + /// TTask: Clears all data structures created by a previous execution of the task + /// \param option Specify action for the TTask::Exec (not defined yet) + void Clear(Option_t* option = ""); + + /// Ttask: Executes the task (processes a timeslice) + /// \param option Specify action for the TTask::Exec (not defined yet) + void Exec(Option_t* /*option*/); + + /// FairTask: Action at hte end of the run + void Finish(); + + //TFolder& GetQa(); + + /// FairTask: Task initialization an the beginning of the run + InitStatus Init(); + + /// FairTask: Task reinitialization + InitStatus ReInit() { return Init(); } + + /// FairTask: Set parameter containers (if any) + void SetParContainers(); + + /// Resets data fields + void Reset(); + + /// Gets pointer to hits array of a particular detector type + virtual TClonesArray* GetHitsInput() = 0; + + /// Gets pointer to clusters array of a particular detector type + virtual TClonesArray* GetClustersInput() = 0; + + /// Gets pointer to hit matches array for a particular detector type + virtual TClonesArray* GetHitMatchesInput() = 0; + + /// Gets pointer to MC points for a particular detector type + virtual CbmMCDataArray* GetMcPointsInput() = 0; + + ClassDef(CbmTrackingInputQaBase, 0); // Base class for tracking input QA + +protected: + /// Registers histogram + void RegisterHist(TH1* histo) + { + fHistList.push_back(histo); + histo->SetDirectory(0); + } + + /// Checks geometry + InitStatus GeometryQa(); // TODO: Final or not final? (S.Zharko) + + /// Checks resolution + /// A general scheme of resolution checks is provided in this method. One should specify particular operations in the + /// derived classes for each particular detector.... + void ResolutionQa(); // TODO: Final or not final? (S.Zharko) + + + // *********** + // ** Setup ** + // *********** + + + CbmTrackingDetectorInterfaceBase* fpDetectorInterface {nullptr}; ///< Pointer to current tracking detector I/F + CbmTimeSlice* fpTimeSlice {nullptr}; ///< Pointer to current time slice + CbmDigiManager* fpDigiManager {nullptr}; ///< Pointer to digi manager + + CbmMCEventList* fpMcEventList {nullptr}; + CbmMCDataManager* fpMcManager {nullptr}; + + // ********** + // ** Data ** + // ********** + + TClonesArray* fpClusters {nullptr}; + TClonesArray* fpHits {nullptr}; + TClonesArray* fpHitMatches {nullptr}; + + CbmMCDataArray* fpMcTracks {nullptr}; + CbmMCDataArray* fpMcPoints {nullptr}; + + // ******************************** + // ** List of base QA histograms ** + // ******************************** + + // NOTE: One can add extra detector specific histograms in the Init method of the derived class. In that case, the + // histograms should be registered via RegisterHisto method + + std::vector<CbmQaHist<TH1D>> fHistResidualX; ///< Residual distributions for X vs index of a station + std::vector<CbmQaHist<TH1D>> fHistResidualY; ///< Residual distributions for X vs index of a station + std::vector<CbmQaHist<TH1D>> fHistResidualT; ///< Residual distributions for T vs index of a station + + std::vector<CbmQaHist<TH1D>> fHistPullX; ///< Pulls distribution for X vs. index of a station + std::vector<CbmQaHist<TH1D>> fHistPullY; ///< Pulls distribution for Y vs. index of a station + std::vector<CbmQaHist<TH1D>> fHistPullT; ///< Pulls distribution for T vs. index of a station + + std::vector<CbmQaHist<TH1I>> fHistPointsPerHit; ///< Distribution of MC points per one hit vs. index of a station + std::vector<CbmQaHist<TH1I>> fHistHitsPerPoint; ///< Distribution of hits per one MC point vs. index of a station + + // Histogram properties, can be modified in the constructor of a derived class + std::tuple<int, double, double> fRangeResidualX {100, -5., 5.}; + std::tuple<int, double, double> fRangeResidualY {100, -10., 10.}; + std::tuple<int, double, double> fRangeResidualT {100, -50., 50.}; + + std::tuple<int, double, double> fRangePullX {100, -5., 5.}; + std::tuple<int, double, double> fRangePullY {100, -5., 5.}; + std::tuple<int, double, double> fRangePullT {100, -5., 5.}; + + std::tuple<int, double, double> fRangePointsPerHit {10, -0.5, 9.5}; + std::tuple<int, double, double> fRangeHitsPerPoint {10, -0.5, 9.5}; + + TParameter<long> fNevents {"nEvents", 0}; ///< Number of processed events + + // ************ + // ** Output ** + // ************ + + std::shared_ptr<TFolder> fpOutFolder {nullptr}; ///< Output folder for this QA task + TFolder* fpOutFolderHists {nullptr}; ///< Subfolder for raw histograms + + ECbmModuleId fDetId; ///< ID of the detector subsystem + TString fDetName {""}; ///< Name of the detector subsystem, used for i/o (e.g., mvd, sts, much, trd1d, trd2d, tof) + TString fDetTitle {""}; ///< Title of the detector subsystem, used for logs (e.g., MVD, STS, MuCh, TRD-1D, ...) + +private: + /// Initializes base histograms + InitStatus InitBaseHistograms(); + + /// Initializes input data arrays + InitStatus InitInputData(); + + std::vector<TH1*> fHistList; ///< List of the pointers to all histograms contained in the class + + bool fIsMcPresent {false}; +}; + +#endif // CbmTrackingInputQaBase_h diff --git a/reco/L1/qa/CbmTrackingInputQaSts.cxx b/reco/L1/qa/CbmTrackingInputQaSts.cxx new file mode 100644 index 0000000000000000000000000000000000000000..b02229e00f732c8b72312839b35ced78bf8af79c --- /dev/null +++ b/reco/L1/qa/CbmTrackingInputQaSts.cxx @@ -0,0 +1,69 @@ +/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergey Gorbunov, Sergei Zharko [committer] */ + +/*************************************************************************************************** + * @file CbmTrackingInputQaSts.cxx + * @brief Class for the tracking input QA for STS (definition) + * @since 27.06.2022 + * @author S.Zharko <s.zharko@gsi.de> + ***************************************************************************************************/ + +#include "CbmTrackingInputQaSts.h" + +#include "CbmDigiManager.h" +#include "CbmMCDataArray.h" +#include "CbmMCEventList.h" +#include "CbmMatch.h" +#include "CbmSetup.h" +#include "CbmStsCluster.h" +#include "CbmStsDigi.h" +#include "CbmStsHit.h" +#include "CbmStsTrackingInterface.h" + +ClassImp(CbmTrackingInputQaSts); + +// -------------------------------------------------------------------------------------------------------------------- +// +CbmTrackingInputQaSts::CbmTrackingInputQaSts(int verbose) : CbmTrackingInputQaBase("CbmTrackingInputQaSts", verbose) +{ + std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaSts::CbmTrackingInputQaSts()\n"; + + // Setup detector-specific parameters + fDetId = ECbmModuleId::kSts; + fDetName = "sts"; + fDetTitle = "STS"; + fpDetectorInterface = CbmStsTrackingInterface::Instance(); +} + +// -------------------------------------------------------------------------------------------------------------------- +// +CbmTrackingInputQaSts::~CbmTrackingInputQaSts() {} + +// -------------------------------------------------------------------------------------------------------------------- +// +InitStatus CbmTrackingInputQaSts::Init() +{ + std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaSts::Init()\n"; + // Check, if STS is in the setup + if (!CbmSetup::Instance()->IsActive(ECbmModuleId::kSts)) { + LOG(warn) << "STS is not active in the current setup"; + return kSUCCESS; + } + + // Initialize base histograms + CbmTrackingInputQaBase::Init(); + + // Additional STS-specific histograms can be initialized here + // ... + + return kSUCCESS; +} + +// -------------------------------------------------------------------------------------------------------------------- +// + + +// -------------------------------------------------------------------------------------------------------------------- + +// -------------------------------------------------------------------------------------------------------------------- diff --git a/reco/L1/qa/CbmTrackingInputQaSts.h b/reco/L1/qa/CbmTrackingInputQaSts.h new file mode 100644 index 0000000000000000000000000000000000000000..f866655a9982e917885c56d87d8038abae2c6aa7 --- /dev/null +++ b/reco/L1/qa/CbmTrackingInputQaSts.h @@ -0,0 +1,55 @@ +/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergey Gorbunov, Sergei Zharko [committer] */ + +/*************************************************************************************************** + * @file CbmTrackingInputQaSts.h + * @brief Class for the tracking input QA for STS (declaration) + * @since 27.06.2022 + * @author S.Zharko <s.zharko@gsi.de> + ***************************************************************************************************/ + +#ifndef CbmTrackingInputQaSts_h +#define CbmTrackingInputQaSts_h 1 + +#include "CbmTrackingInputQaBase.h" + +#include "FairRootManager.h" + +#include "TClonesArray.h" + +class CbmMCDataArray; +class CbmMatch; +class CbmCluster; + +class CbmTrackingInputQaSts : public CbmTrackingInputQaBase { +public: + CbmTrackingInputQaSts(int verbose = 1); + ~CbmTrackingInputQaSts(); + + InitStatus Init(); + + + /// Gets pointer to hits array of a particular detector type + TClonesArray* GetHitsInput() { return dynamic_cast<TClonesArray*>(FairRootManager::Instance()->GetObject("StsHit")); } + + /// Gets pointer to clusters array of a particular detector type + TClonesArray* GetClustersInput() + { + return dynamic_cast<TClonesArray*>(FairRootManager::Instance()->GetObject("StsCluster")); + } + + /// Gets pointer to hit matches array for a particular detector type + TClonesArray* GetHitMatchesInput() + { + return dynamic_cast<TClonesArray*>(FairRootManager::Instance()->GetObject("StsHitMatch")); + } + + /// Gets pointer to MC points for a particular detector type + CbmMCDataArray* GetMcPointsInput() { return fpMcManager->InitBranch("StsPoint"); } + + + ClassDef(CbmTrackingInputQaSts, 0); +}; + +#endif // CbmTrackingInputQaSts_h