From 5812bf3628f37c9a86d6e88a91f64948f7b0e1c9 Mon Sep 17 00:00:00 2001 From: "se.gorbunov" <se.gorbunov@gsi.de> Date: Sun, 21 Aug 2022 21:46:35 +0000 Subject: [PATCH] L1: qa task for trd2d tracking --- macro/L1/run_reco_L1global.C | 9 + reco/qa/CMakeLists.txt | 62 ++- reco/qa/CbmTrackingTrdQa.cxx | 843 +++++++++++++++++++++++++++++++++++ reco/qa/CbmTrackingTrdQa.h | 231 ++++++++++ reco/qa/RecoQaLinkDef.h | 1 + 5 files changed, 1141 insertions(+), 5 deletions(-) create mode 100644 reco/qa/CbmTrackingTrdQa.cxx create mode 100644 reco/qa/CbmTrackingTrdQa.h diff --git a/macro/L1/run_reco_L1global.C b/macro/L1/run_reco_L1global.C index 05cfa152d7..af0f6e5fa2 100644 --- a/macro/L1/run_reco_L1global.C +++ b/macro/L1/run_reco_L1global.C @@ -393,6 +393,15 @@ void run_reco_L1global(TString input = "", Int_t nTimeSlices = -1, Int_t firstTi FairTask* globalFindTracks = new CbmL1GlobalFindTracksEvents(globalTrackFinder); run->AddTask(globalFindTracks); std::cout << "-I- " << myName << ": Added task " << globalFindTracks->GetName() << std::endl; + + if (debugWithMC) { + CbmMatchRecoToMC* match2 = new CbmMatchRecoToMC(); + run->AddTask(match2); + } + + CbmTrackingTrdQa* trdTrackerQa = new CbmTrackingTrdQa; + run->AddTask(trdTrackerQa); + std::cout << "-I- " << myName << ": Added task " << trdTrackerQa->GetName() << std::endl; } // ------------------------------------------------------------------------ diff --git a/reco/qa/CMakeLists.txt b/reco/qa/CMakeLists.txt index 8fb6c18950..d64eef170d 100644 --- a/reco/qa/CMakeLists.txt +++ b/reco/qa/CMakeLists.txt @@ -1,17 +1,54 @@ -Set(INCLUDE_DIRECTORIES -${CMAKE_CURRENT_SOURCE_DIR} +# ---- Specify include directories ----------------------- +set(INCLUDE_DIRECTORIES + +# This directory +${CBMROOT_SOURCE_DIR}/reco/qa + +# external +${CBMROOT_SOURCE_DIR}/external/xpu/src + +# Reco +${CBMROOT_SOURCE_DIR}/reco/base + +#QA +${CBMROOT_SOURCE_DIR}/core/qa + +# Base ${CBMBASE_DIR} -${CBMDATA_DIR} -${CBMDATA_DIR}/mvd +# Data +${CBMDATA_DIR} +${CBMDATA_DIR}/base +${CBMDATA_DIR}/global ${CBMDATA_DIR}/sts +${CBMDATA_DIR}/raw +${CBMDATA_DIR}/mvd ${CBMDATA_DIR}/much ${CBMDATA_DIR}/trd ${CBMDATA_DIR}/tof + +# MVD +${CBMROOT_SOURCE_DIR}/mvd + +# STS +${CBMDETECTORBASE_DIR}/sts +${CBMDETECTORBASE_DIR}/trd +${CBMROOT_SOURCE_DIR}/reco/detectors/sts +${CBMROOT_SOURCE_DIR}/reco/detectors/sts/unpack +${CBMROOT_SOURCE_DIR}/reco/detectors/sts/qa +${CBMROOT_SOURCE_DIR}/reco/detectors/sts/experimental + ) +set(SYSTEM_INCLUDE_DIRECTORIES +${BASE_INCLUDE_DIRECTORIES} +${FLES_IPC_INCLUDE_DIRECTORY} # for fles infos for unpacker +) +# ---- End of include directories ------------------------ + + Set(SYSTEM_INCLUDE_DIRECTORIES ${BASE_INCLUDE_DIRECTORIES} ${ROOT_INCLUDE_DIR} @@ -32,6 +69,7 @@ link_directories(${LINK_DIRECTORIES}) set(SRCS CbmRecoQa.cxx + CbmTrackingTrdQa.cxx ) set(LINKDEF RecoQaLinkDef.h) @@ -39,7 +77,21 @@ set(LINKDEF RecoQaLinkDef.h) Set(LIBRARY_NAME CbmRecoQa) Set(DEPENDENCIES - CbmBase CbmData Base Core Hist + CbmBase CbmData Base Core Hist + + Base + CbmBase + CbmData + CbmSimSteer +# CbmGeoSetup + CbmTrdBase + CbmMvd + CbmStsBase + CbmRecoBase + CbmRecoSts + CbmQaBase + Core + Hist ) GENERATE_LIBRARY() diff --git a/reco/qa/CbmTrackingTrdQa.cxx b/reco/qa/CbmTrackingTrdQa.cxx new file mode 100644 index 0000000000..b2c5c65459 --- /dev/null +++ b/reco/qa/CbmTrackingTrdQa.cxx @@ -0,0 +1,843 @@ +/* Copyright (C) 2006-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergey Gorbunov, Denis Bertini [committer], Volker Friese, Florian Uhlig */ + +// ------------------------------------------------------------------------- +// ----- CbmTrackingTrdQa source file ----- +// ----- Created 11/01/06 by V. Friese ----- +// ------------------------------------------------------------------------- + +// Includes class header +#include "CbmTrackingTrdQa.h" + +// Includes from C++ +#include <cassert> +#include <iomanip> + +// Includes from ROOT +#include "TClonesArray.h" +#include "TGeoManager.h" +#include "TH1F.h" +#include "TH2F.h" + +// Includes from FairRoot +#include "FairEventHeader.h" +#include "FairRun.h" +#include <Logger.h> + +// Includes from CbmRoot +#include "CbmGlobalTrack.h" +#include "CbmMCDataArray.h" +#include "CbmMCDataManager.h" +#include "CbmMCTrack.h" +#include "CbmTimeSlice.h" +#include "CbmTrackMatchNew.h" +#include "CbmTrdHit.h" +#include "CbmTrdPoint.h" +#include "CbmTrdTrack.h" +#include "CbmTrdTrackingInterface.h" + +#include "FairRunAna.h" + +using std::fixed; +using std::right; +using std::setprecision; +using std::setw; + + +// ----- Default constructor ------------------------------------------- +CbmTrackingTrdQa::CbmTrackingTrdQa(Int_t iVerbose) : FairTask("CbmTrackingTrdQa", iVerbose) {} + +// ------------------------------------------------------------------------- + + +// ----- Standard constructor ------------------------------------------ +CbmTrackingTrdQa::CbmTrackingTrdQa(Int_t minStations, Double_t quota, Int_t iVerbose) + : FairTask("CbmTrackingTrdQa", iVerbose) + , fMinStations(minStations) + , fQuota(quota) +{ +} +// ------------------------------------------------------------------------- + + +// ----- Destructor ---------------------------------------------------- +CbmTrackingTrdQa::~CbmTrackingTrdQa() +{ + fHistoList->Delete(); + delete fHistoList; +} +// ------------------------------------------------------------------------- + + +// ----- Task execution ------------------------------------------------ +void CbmTrackingTrdQa::Exec(Option_t* /*opt*/) +{ + + LOG(debug) << GetName() << ": Process event "; + + // Timer + fTimer.Start(); + + // Eventwise counters + // Int_t nMCTracks = 0; + Int_t nTracks = 0; + Int_t nGhosts = 0; + Int_t nClones = 0; + Int_t nAll = 0; + Int_t nAcc = 0; + Int_t nRecAll = 0; + Int_t nPrim = 0; + Int_t nRecPrim = 0; + Int_t nFast = 0; + Int_t nRecFast = 0; + Int_t nFastLong = 0; + Int_t nRecFastLong = 0; + Int_t nSec = 0; + Int_t nRecSec = 0; + TVector3 momentum; + TVector3 vertex; + + // check consistency + //assert(fGlobalTracks->GetEntriesFast() == fGlobalTrackMatches->GetEntriesFast()); + assert(fTrdTracks->GetEntriesFast() == fTrdTrackMatches->GetEntriesFast()); + + std::vector<CbmLink> events = fTimeSlice->GetMatch().GetLinks(); + std::sort(events.begin(), events.end()); + + { + fMcTrackInfoMap.clear(); + McTrackInfo info; + info.fIsAccepted = false; + info.fNtimesReconstructed = 0; + info.fIsPrimary = false; + info.fIsFast = false; + info.fIsLong = false; + + std::cout << "MC events in slice : " << events.size() << std::endl; + + for (uint iLink = 0; iLink < events.size(); iLink++) { + CbmLink link = events[iLink]; + Int_t nMCTracks = fMCTracks->Size(link); + std::cout << "MC event " << iLink << " n mc tracks " << nMCTracks << std::endl; + for (Int_t iTr = 0; iTr < nMCTracks; iTr++) { + link.SetIndex(iTr); + fMcTrackInfoMap.insert({link, info}); + } + } + } + + // Fill hit and track maps + FillHitMap(); + FillTrackMatchMap(nTracks, nGhosts, nClones); + + int nMcTracks = fMcTrackInfoMap.size(); + + // Loop over MCTracks + int iMcTrack = 0; + for (auto itTrack = fMcTrackInfoMap.begin(); itTrack != fMcTrackInfoMap.end(); ++itTrack, ++iMcTrack) { + const CbmLink& link = itTrack->first; + McTrackInfo& info = itTrack->second; + CbmMCTrack* mcTrack = dynamic_cast<CbmMCTrack*>(fMCTracks->Get(link)); + assert(mcTrack); + + nAll++; + + // Check origin of MCTrack + // TODO: Track origin should rather be compared to MC event vertex + // But that is not available from MCDataManager + mcTrack->GetStartVertex(vertex); + if (TMath::Abs(vertex.Z() - fTargetPos.Z()) < 1.) { info.fIsPrimary = true; } + + // Get momentum + mcTrack->GetMomentum(momentum); + info.fPt = momentum.Pt(); + info.fIsFast = (info.fPt > 0.1); + + // Continue only for reconstructible tracks + + Int_t nStations = info.fHitMap.size(); + + int nContStations = 0; // Number of continious stations + { + int istaprev = -1; + int len = 0; + for (auto itSta = info.fHitMap.begin(); itSta != info.fHitMap.end(); itSta++) { + if (len == 0 || itSta->first == istaprev + 1) { len++; } + else { + len = 1; + } + if (nContStations < len) { nContStations = len; } + istaprev = itSta->first; + } + } + + info.fIsLong = (nContStations >= fTrdNstations); + + if (nStations < fMinStations) continue; // Too few stations + if (nContStations < fMinStations) continue; // Too few stations + + info.fIsAccepted = true; + nAcc++; + + if (info.fIsPrimary) { nPrim++; } + else { + nSec++; + } + + if (info.fIsFast) { nFast++; } + + Bool_t isFastLong = (info.fIsFast && info.fIsLong); + if (isFastLong) { nFastLong++; } + + // Fill histograms for reconstructible tracks + + fhPtAccAll->Fill(info.fPt); + fhNpAccAll->Fill(Double_t(nStations)); + if (info.fIsPrimary) { + fhPtAccPrim->Fill(info.fPt); + fhNpAccPrim->Fill(Double_t(nStations)); + } + else { + fhPtAccSec->Fill(info.fPt); + fhNpAccSec->Fill(Double_t(nStations)); + fhZAccSec->Fill(vertex.Z()); + } + + // Get matched GlobalTrack + Int_t globalTrackId = info.fGlobalTrackMatch; + Int_t trdTrackId = info.fTrdTrackMatch; + Double_t quali = info.fQuali; + // Bool_t isRec = kFALSE; + if (globalTrackId < 0) continue; + if (trdTrackId < 0) continue; + CbmTrdTrack* trdTrack = (CbmTrdTrack*) fTrdTracks->At(trdTrackId); + assert(trdTrack); + assert(quali >= fQuota); + CbmTrackMatchNew* match = (CbmTrackMatchNew*) fTrdTrackMatches->At(trdTrackId); + assert(match); + Int_t nTrue = match->GetNofTrueHits(); + Int_t nWrong = match->GetNofWrongHits(); + //Int_t nFake = match->GetNofFakeHits(); + Int_t nFake = 0; + Int_t nAllHits = trdTrack->GetNofHits(); + assert(nTrue + nWrong + nFake == nAllHits); + // Verbose output + LOG(debug1) << GetName() << ": MCTrack " << iMcTrack << ", stations " << nStations << ", hits " << nAllHits + << ", true hits " << nTrue; + + // Fill histograms for reconstructed tracks + nRecAll++; + fhPtRecAll->Fill(info.fPt); + fhNpRecAll->Fill(Double_t(nAllHits)); + if (info.fIsPrimary) { + nRecPrim++; + fhPtRecPrim->Fill(info.fPt); + fhNpRecPrim->Fill(Double_t(nAllHits)); + } + else { + nRecSec++; + fhPtRecSec->Fill(info.fPt); + fhNpRecSec->Fill(Double_t(nAllHits)); + fhZRecSec->Fill(vertex.Z()); + } + if (info.fIsFast) nRecFast++; + if (isFastLong) nRecFastLong++; + + } // Loop over MCTracks + + // Loop over MC points + + for (uint iLink = 0; iLink < events.size(); iLink++) { + CbmLink link = events[iLink]; + Int_t nMcPoints = fTrdPoints->Size(link); + std::cout << "MC event " << iLink << " n mc points " << nMcPoints << std::endl; + for (Int_t ip = 0; ip < nMcPoints; ip++) { + link.SetIndex(ip); + CbmTrdPoint* pt = (CbmTrdPoint*) fTrdPoints->Get(link); + assert(pt); + Int_t station = CbmTrdTrackingInterface::Instance()->GetTrackingStationIndex(pt->GetDetectorID()); + link.SetIndex(pt->GetTrackID()); + McTrackInfo& info = getMcTrackInfo(link); + if (info.fIsAccepted) { + fhStationEffXY[station].Fill(pt->GetX(), pt->GetY(), (info.fNtimesReconstructed > 0) ? 100. : 0.); + } + } + } + + // Calculate efficiencies + Double_t effAll = Double_t(nRecAll) / Double_t(nAcc); + Double_t effPrim = Double_t(nRecPrim) / Double_t(nPrim); + Double_t effFast = Double_t(nRecFast) / Double_t(nFast); + Double_t effFastLong = Double_t(nRecFastLong) / Double_t(nFastLong); + Double_t effSec = Double_t(nRecSec) / Double_t(nSec); + + fTimer.Stop(); + + + // Event summary + LOG(info) << "+ " << setw(20) << GetName() << ": Event " << setw(6) << right << fNEvents << ", real time " << fixed + << setprecision(6) << fTimer.RealTime() << " s, MC tracks: all " << nMcTracks << ", acc. " << nAcc + << ", rec. " << nRecAll << ", eff. " << setprecision(2) << 100. * effAll << " %"; + if (fair::Logger::Logging(fair::Severity::debug)) { + LOG(debug) << "---------- CbmTrackingTrdQa : Event summary ------------"; + LOG(debug) << "MCTracks : " << nAll << ", reconstructible: " << nAcc << ", reconstructed: " << nRecAll; + LOG(debug) << "Vertex : reconstructible: " << nPrim << ", reconstructed: " << nRecPrim << ", efficiency " + << effPrim * 100. << "%"; + LOG(debug) << "Fast : reconstructible: " << nFast << ", reconstructed: " << nRecFast << ", efficiency " + << effFast * 100. << "%"; + LOG(debug) << "Fast long : reconstructible: " << nFastLong << ", reconstructed: " << nRecFastLong << ", efficiency " + << effFastLong * 100. << "%"; + LOG(debug) << "Non-vertex : reconstructible: " << nSec << ", reconstructed: " << nRecSec << ", efficiency " + << effSec * 100. << "%"; + LOG(debug) << "TrdTracks " << nTracks << ", ghosts " << nGhosts << ", clones " << nClones; + LOG(debug) << "-----------------------------------------------------------\n"; + } + + + // Increase counters + fNAll += nAll; + fNAccAll += nAcc; + fNAccPrim += nPrim; + fNAccFast += nFast; + fNAccFastLong += nFastLong; + fNAccSec += nSec; + fNRecAll += nRecAll; + fNRecPrim += nRecPrim; + fNRecFast += nRecFast; + fNRecFastLong += nRecFastLong; + fNRecSec += nRecSec; + fNGhosts += nGhosts; + fNClones += nClones; + fNEvents++; + fTime += fTimer.RealTime(); +} +// ------------------------------------------------------------------------- + + +// ----- Public method SetParContainers -------------------------------- +void CbmTrackingTrdQa::SetParContainers() {} +// ------------------------------------------------------------------------- + + +// ----- Public method Init -------------------------------------------- +InitStatus CbmTrackingTrdQa::Init() +{ + + LOG(info) << "\n\n===================================================="; + LOG(info) << GetName() << ": Initialising..."; + + fManager = FairRootManager::Instance(); + assert(fManager); + + fMcManager = dynamic_cast<CbmMCDataManager*>(fManager->GetObject("MCDataManager")); + + assert(fMcManager); + + fTimeSlice = static_cast<CbmTimeSlice*>(fManager->GetObject("TimeSlice.")); + + if (fTimeSlice == nullptr) { LOG(fatal) << "CbmTrackingTrdQa: No time slice object"; } + + if (fMcManager) { + fMCTracks = fMcManager->InitBranch("MCTrack"); + fTrdPoints = fMcManager->InitBranch("TrdPoint"); + } + + assert(fMCTracks); + assert(fTrdPoints); + + // Get the geometry + InitStatus geoStatus = GetGeometry(); + if (geoStatus != kSUCCESS) { + LOG(error) << GetName() << "::Init: Error in reading geometry!"; + return geoStatus; + } + + + // TRD + + // Get TrdHit array + fTrdHits = (TClonesArray*) fManager->GetObject("TrdHit"); + assert(fTrdHits); + + // Get TrdHitMatch array + fTrdHitMatch = (TClonesArray*) fManager->GetObject("TrdHitMatch"); + assert(fTrdHitMatch); + + // Get GlobalTrack array + fGlobalTracks = (TClonesArray*) fManager->GetObject("GlobalTrack"); + assert(fGlobalTracks); + + // Get GlobalTrackMatch array + fGlobalTrackMatches = (TClonesArray*) fManager->GetObject("GlobalTrackMatch"); + //assert(fGlobalTrackMatches); + + // Get TrdTrack array + fTrdTracks = (TClonesArray*) fManager->GetObject("TrdTrack"); + assert(fTrdTracks); + + // Get TrdTrackMatch array + fTrdTrackMatches = (TClonesArray*) fManager->GetObject("TrdTrackMatch"); + assert(fTrdTrackMatches); + + + // Create histograms + CreateHistos(); + Reset(); + + // Output + LOG(info) << " Number of Trd stations : " << fTrdNstations; + LOG(info) << " Target position ( " << fTargetPos.X() << ", " << fTargetPos.Y() << ", " << fTargetPos.Z() << ") cm"; + LOG(info) << " Minimum number of Trd stations : " << fMinStations; + LOG(info) << " Matching quota : " << fQuota; + LOG(info) << "===================================================="; + + return geoStatus; +} +// ------------------------------------------------------------------------- + + +// ----- Public method ReInit ------------------------------------------ +InitStatus CbmTrackingTrdQa::ReInit() +{ + + LOG(info) << "\n\n===================================================="; + LOG(info) << GetName() << ": Re-initialising..."; + + // Get the geometry of target and Trd + InitStatus geoStatus = GetGeometry(); + if (geoStatus != kSUCCESS) { + LOG(error) << GetName() << "::Init: Error in reading geometry!"; + return geoStatus; + } + + // --- Screen log + LOG(info) << " Number of TRD stations : " << fTrdNstations; + LOG(info) << " Target position ( " << fTargetPos.X() << ", " << fTargetPos.Y() << ", " << fTargetPos.Z() << ") cm"; + LOG(info) << " Minimum number of TRD stations : " << fMinStations; + LOG(info) << " Matching quota : " << fQuota; + LOG(info) << "===================================================="; + + return geoStatus; +} +// ------------------------------------------------------------------------- + + +// ----- Private method Finish ----------------------------------------- +void CbmTrackingTrdQa::Finish() +{ + + // Divide histograms for efficiency calculation + DivideHistos(fhPtRecAll, fhPtAccAll, fhPtEffAll); + DivideHistos(fhPtRecPrim, fhPtAccPrim, fhPtEffPrim); + DivideHistos(fhPtRecSec, fhPtAccSec, fhPtEffSec); + DivideHistos(fhNpRecAll, fhNpAccAll, fhNpEffAll); + DivideHistos(fhNpRecPrim, fhNpAccPrim, fhNpEffPrim); + DivideHistos(fhNpRecSec, fhNpAccSec, fhNpEffSec); + DivideHistos(fhZRecSec, fhZAccSec, fhZEffSec); + + // Normalise histos for clones and ghosts to one event + if (fNEvents) { + fhNhClones->Scale(1. / Double_t(fNEvents)); + fhNhGhosts->Scale(1. / Double_t(fNEvents)); + } + + // Calculate integrated efficiencies and rates + Double_t effAll = Double_t(fNRecAll) / Double_t(fNAccAll); + Double_t effPrim = Double_t(fNRecPrim) / Double_t(fNAccPrim); + Double_t effFast = Double_t(fNRecFast) / Double_t(fNAccFast); + Double_t effFastLong = Double_t(fNRecFastLong) / Double_t(fNAccFastLong); + Double_t effSec = Double_t(fNRecSec) / Double_t(fNAccSec); + Double_t rateGhosts = Double_t(fNGhosts) / Double_t(fNRecAll); + Double_t rateClones = Double_t(fNClones) / Double_t(fNRecAll); + + // Run summary to screen + std::cout << std::endl; + LOG(info) << "====================================="; + LOG(info) << fName << ": Run summary "; + LOG(info) << "Events processed : " << fNEvents << setprecision(2); + LOG(info) << "Eff. all tracks : " << effAll * 100 << " % (" << fNRecAll << "/" << fNAccAll << ")"; + LOG(info) << "Eff. vertex tracks : " << effPrim * 100 << " % (" << fNRecPrim << "/" << fNAccPrim << ")"; + LOG(info) << "Eff. fast tracks : " << effFast * 100 << " % (" << fNRecFast << "/" << fNAccFast << ")"; + LOG(info) << "Eff. fast long tracks : " << effFastLong * 100 << " % (" << fNRecFastLong << "/" << fNAccFastLong + << ")"; + LOG(info) << "Eff. secondary tracks : " << effSec * 100 << " % (" << fNRecSec << "/" << fNAccSec << ")"; + LOG(info) << "Ghost rate : " << rateGhosts * 100 << " % (" << fNGhosts << "/" << fNRecAll << ")"; + LOG(info) << "Clone rate : " << rateClones * 100 << " % (" << fNClones << "/" << fNRecAll << ")"; + LOG(info) << "mc tracks/event " << fNAll / fNEvents << " accepted " << fNRecAll / fNEvents; + LOG(info) << "Time per event : " << setprecision(6) << fTime / Double_t(fNEvents) << " s"; + + + LOG(info) << "====================================="; + + // Write histos to output + /* + gDirectory->mkdir("CbmTrackingTrdQa"); + gDirectory->cd("CbmTrackingTrdQa"); + TIter next(fHistoList); + while (TH1* histo = ((TH1*) next())) + histo->Write(); + gDirectory->cd(".."); +*/ + FairSink* sink = FairRootManager::Instance()->GetSink(); + sink->WriteObject(&fOutFolder, nullptr); +} +// ------------------------------------------------------------------------- + + +// ----- Private method GetGeometry ------------------------------------ +InitStatus CbmTrackingTrdQa::GetGeometry() +{ + // Get target geometry + GetTargetPosition(); + + // Get TRD setup + auto trdInterface = CbmTrdTrackingInterface::Instance(); + assert(trdInterface); + fTrdNstations = trdInterface->GetNtrackingStations(); + return kSUCCESS; +} +// ------------------------------------------------------------------------- + + +// ----- Get target node ----------------------------------------------- +void CbmTrackingTrdQa::GetTargetPosition() +{ + + TGeoNode* target = NULL; + + gGeoManager->CdTop(); + TGeoNode* cave = gGeoManager->GetCurrentNode(); + for (Int_t iNode1 = 0; iNode1 < cave->GetNdaughters(); iNode1++) { + TString name = cave->GetDaughter(iNode1)->GetName(); + if (name.Contains("pipe", TString::kIgnoreCase)) { + LOG(debug) << "Found pipe node " << name; + gGeoManager->CdDown(iNode1); + break; + } + } + for (Int_t iNode2 = 0; iNode2 < gGeoManager->GetCurrentNode()->GetNdaughters(); iNode2++) { + TString name = gGeoManager->GetCurrentNode()->GetDaughter(iNode2)->GetName(); + if (name.Contains("pipevac1", TString::kIgnoreCase)) { + LOG(debug) << "Found vacuum node " << name; + gGeoManager->CdDown(iNode2); + break; + } + } + for (Int_t iNode3 = 0; iNode3 < gGeoManager->GetCurrentNode()->GetNdaughters(); iNode3++) { + TString name = gGeoManager->GetCurrentNode()->GetDaughter(iNode3)->GetName(); + if (name.Contains("target", TString::kIgnoreCase)) { + LOG(debug) << "Found target node " << name; + gGeoManager->CdDown(iNode3); + target = gGeoManager->GetCurrentNode(); + break; + } + } + if (!target) { + fTargetPos[0] = 0.; + fTargetPos[1] = 0.; + fTargetPos[2] = 0.; + } + else { + TGeoHMatrix* glbMatrix = gGeoManager->GetCurrentMatrix(); + Double_t* pos = glbMatrix->GetTranslation(); + fTargetPos[0] = pos[0]; + fTargetPos[1] = pos[1]; + fTargetPos[2] = pos[2]; + } + + gGeoManager->CdTop(); +} +// ------------------------------------------------------------------------- + + +// ----- Private method CreateHistos ----------------------------------- +void CbmTrackingTrdQa::CreateHistos() +{ + + fOutFolder.Clear(); + + // Histogram list + fHistoList = new TList(); + + // Momentum distributions + Double_t minPt = 0.; + Double_t maxPt = 10.; + Int_t nBinsPt = 40; + + fhStationEffXY.clear(); + + for (int i = 0; i < fTrdNstations; i++) { + //auto trdInterface = CbmTrdTrackingInterface::Instance(); + double dx = 150; //trdInterface->GetXmax(i); + double dy = 150; //trdInterface->GetYmax(i); + fhStationEffXY.emplace_back(Form("fhStationEffXY%i", i), Form("Efficiency XY: Station %i;X [cm];Y [cm]", i), 50, + -dx, dx, 50, -dy, dy); + fhStationEffXY[i].SetDirectory(0); + fhStationEffXY[i].SetOptStat(10); + fhStationEffXY[i].GetYaxis()->SetTitleOffset(1.4); + } + + for (int i = 0; i < fTrdNstations; i++) { + fHistoList->Add(&fhStationEffXY[i]); + } + + fhPtAccAll = new TH1F("hPtAccAll", "all reconstructable tracks", nBinsPt, minPt, maxPt); + fhPtRecAll = new TH1F("hPtRecAll", "all reconstructed tracks", nBinsPt, minPt, maxPt); + fhPtEffAll = new TH1F("hPtEffAll", "efficiency all tracks", nBinsPt, minPt, maxPt); + fhPtAccPrim = new TH1F("hPtAccPrim", "reconstructable vertex tracks", nBinsPt, minPt, maxPt); + fhPtRecPrim = new TH1F("hPtRecPrim", "reconstructed vertex tracks", nBinsPt, minPt, maxPt); + fhPtEffPrim = new TH1F("hPtEffPrim", "efficiency vertex tracks", nBinsPt, minPt, maxPt); + fhPtAccSec = new TH1F("hPtAccSec", "reconstructable non-vertex tracks", nBinsPt, minPt, maxPt); + fhPtRecSec = new TH1F("hPtRecSec", "reconstructed non-vertex tracks", nBinsPt, minPt, maxPt); + fhPtEffSec = new TH1F("hPtEffSec", "efficiency non-vertex tracks", nBinsPt, minPt, maxPt); + fHistoList->Add(fhPtAccAll); + fHistoList->Add(fhPtRecAll); + fHistoList->Add(fhPtEffAll); + fHistoList->Add(fhPtAccPrim); + fHistoList->Add(fhPtRecPrim); + fHistoList->Add(fhPtEffPrim); + fHistoList->Add(fhPtAccSec); + fHistoList->Add(fhPtRecSec); + fHistoList->Add(fhPtEffSec); + + // Number-of-points distributions + Double_t minNp = -0.5; + Double_t maxNp = 15.5; + Int_t nBinsNp = 16; + fhNpAccAll = new TH1F("hNpAccAll", "all reconstructable tracks", nBinsNp, minNp, maxNp); + fhNpRecAll = new TH1F("hNpRecAll", "all reconstructed tracks", nBinsNp, minNp, maxNp); + fhNpEffAll = new TH1F("hNpEffAll", "efficiency all tracks", nBinsNp, minNp, maxNp); + fhNpAccPrim = new TH1F("hNpAccPrim", "reconstructable vertex tracks", nBinsNp, minNp, maxNp); + fhNpRecPrim = new TH1F("hNpRecPrim", "reconstructed vertex tracks", nBinsNp, minNp, maxNp); + fhNpEffPrim = new TH1F("hNpEffPrim", "efficiency vertex tracks", nBinsNp, minNp, maxNp); + fhNpAccSec = new TH1F("hNpAccSec", "reconstructable non-vertex tracks", nBinsNp, minNp, maxNp); + fhNpRecSec = new TH1F("hNpRecSec", "reconstructed non-vertex tracks", nBinsNp, minNp, maxNp); + fhNpEffSec = new TH1F("hNpEffSec", "efficiency non-vertex tracks", nBinsNp, minNp, maxNp); + fHistoList->Add(fhNpAccAll); + fHistoList->Add(fhNpRecAll); + fHistoList->Add(fhNpEffAll); + fHistoList->Add(fhNpAccPrim); + fHistoList->Add(fhNpRecPrim); + fHistoList->Add(fhNpEffPrim); + fHistoList->Add(fhNpAccSec); + fHistoList->Add(fhNpRecSec); + fHistoList->Add(fhNpEffSec); + + // z(vertex) distributions + Double_t minZ = 0.; + Double_t maxZ = 50.; + Int_t nBinsZ = 50; + fhZAccSec = new TH1F("hZAccSec", "reconstructable non-vertex tracks", nBinsZ, minZ, maxZ); + fhZRecSec = new TH1F("hZRecSecl", "reconstructed non-vertex tracks", nBinsZ, minZ, maxZ); + fhZEffSec = new TH1F("hZEffRec", "efficiency non-vertex tracks", nBinsZ, minZ, maxZ); + fHistoList->Add(fhZAccSec); + fHistoList->Add(fhZRecSec); + fHistoList->Add(fhZEffSec); + + // Number-of-hit distributions + fhNhClones = new TH1F("hNhClones", "number of hits for clones", nBinsNp, minNp, maxNp); + fhNhGhosts = new TH1F("hNhGhosts", "number of hits for ghosts", nBinsNp, minNp, maxNp); + fHistoList->Add(fhNhClones); + fHistoList->Add(fhNhGhosts); + + TIter next(fHistoList); + while (TH1* histo = ((TH1*) next())) { + fOutFolder.Add(histo); + } +} +// ------------------------------------------------------------------------- + + +// ----- Private method Reset ------------------------------------------ +void CbmTrackingTrdQa::Reset() +{ + + TIter next(fHistoList); + while (TH1* histo = ((TH1*) next())) + histo->Reset(); + + fNAccAll = fNAccPrim = fNAccFast = fNAccFastLong = fNAccSec = 0; + fNRecAll = fNRecPrim = fNRecFast = fNRecFastLong = fNRecSec = 0; + fNGhosts = fNClones = fNEvents = 0; +} +// ------------------------------------------------------------------------- + + +// ----- Private method FillHitMap ------------------------------------- +void CbmTrackingTrdQa::FillHitMap() +{ + + // --- Fill hit map ( mcTrack -> ( station -> number of hits ) ) + + // pocess Trd hits + + for (Int_t iHit = 0; iHit < fTrdHits->GetEntriesFast(); iHit++) { + CbmTrdHit* hit = (CbmTrdHit*) fTrdHits->At(iHit); + assert(hit); + + if ((int) hit->GetClassType() != 1) { + // skip TRD-1D hit + continue; + } + + Int_t station = CbmTrdTrackingInterface::Instance()->GetTrackingStationIndex(hit); + + // carefully check the hit match by looking at the strips from both sides + + const CbmMatch* match = dynamic_cast<const CbmMatch*>(fTrdHitMatch->At(iHit)); + assert(match); + if (match->GetNofLinks() <= 0) continue; + + CbmLink link = match->GetMatchedLink(); + + CbmTrdPoint* pt = (CbmTrdPoint*) fTrdPoints->Get(link); + assert(pt); + assert(CbmTrdTrackingInterface::Instance()->GetTrackingStationIndex(pt->GetDetectorID()) == station); + + link.SetIndex(pt->GetTrackID()); + McTrackInfo& info = getMcTrackInfo(link); + info.fHitMap[station]++; + } + LOG(debug) << GetName() << ": Filled hit map from " << fTrdHits->GetEntriesFast() << " Trd hits"; +} +// ------------------------------------------------------------------------- + + +// ------ Private method FillMatchMap ---------------------------------- +void CbmTrackingTrdQa::FillTrackMatchMap(Int_t& nRec, Int_t& nGhosts, Int_t& nClones) +{ + + // Clear matching maps + for (auto it = fMcTrackInfoMap.begin(); it != fMcTrackInfoMap.end(); ++it) { + McTrackInfo& info = it->second; + info.fGlobalTrackMatch = -1; + info.fQuali = 0.; + info.fMatchedNHitsAll = 0; + info.fMatchedNHitsTrue = 0; + info.fNtimesReconstructed = 0; + } + + // Loop over GlobalTracks. Check matched MCtrack and fill maps. + nGhosts = 0; + nClones = 0; + nRec = 0; + + for (Int_t iGlobalTrack = 0; iGlobalTrack < fGlobalTracks->GetEntriesFast(); iGlobalTrack++) { + + // --- GlobalTrack + CbmGlobalTrack* globalTrack = (CbmGlobalTrack*) fGlobalTracks->At(iGlobalTrack); + assert(globalTrack); + + // --- TrackMatch + + assert(iGlobalTrack >= 0 && iGlobalTrack < fTrdTrackMatches->GetEntriesFast()); + //CbmTrackMatchNew* globalMatch = (CbmTrackMatchNew*) fGlobalTrackMatches->At(iGlobalTrack); + //assert(globalMatch); + + + int iTrdTrack = globalTrack->GetTrdTrackIndex(); + if (iTrdTrack < 0) continue; + CbmTrdTrack* trdTrack = (CbmTrdTrack*) fTrdTracks->At(iTrdTrack); + assert(trdTrack); + nRec++; + + // --- TrackMatch + + assert(iTrdTrack >= 0 && iTrdTrack < fTrdTrackMatches->GetEntriesFast()); + CbmTrackMatchNew* trdMatch = (CbmTrackMatchNew*) fTrdTrackMatches->At(iTrdTrack); + assert(trdMatch); + + Int_t nHits = trdTrack->GetNofHits(); + Int_t nTrue = trdMatch->GetNofTrueHits(); + + // Check matching criterion (quota) + Double_t quali = Double_t(nTrue) / Double_t(nHits); + + // Quality isn't good, it's a ghost + + if (quali < fQuota) { + fhNhGhosts->Fill(nHits); + nGhosts++; + continue; + } + + // Quality is good + + // --- Matched MCTrack + assert(trdMatch->GetNofLinks() > 0); + CbmLink link = trdMatch->GetMatchedLink(); + assert(link.GetIndex() >= 0); + link.SetFile(0); //?? + + McTrackInfo& info = getMcTrackInfo(link); + + // previous match is better, this track is a clone + if ((quali < info.fQuali) || ((quali == info.fQuali) && (nTrue < info.fMatchedNHitsTrue))) { + fhNhClones->Fill(nHits); + nClones++; + continue; + } + + // this track is better than the old one + if (info.fMatchedNHitsAll > 0) { + fhNhClones->Fill(info.fMatchedNHitsAll); + nClones++; + } + info.fGlobalTrackMatch = iGlobalTrack; + info.fTrdTrackMatch = iTrdTrack; + info.fQuali = quali; + info.fMatchedNHitsAll = nHits; + info.fMatchedNHitsTrue = nTrue; + info.fNtimesReconstructed++; + } // Loop over GlobalTracks + + LOG(debug) << GetName() << ": Filled match map for " << nRec << " Trd tracks. Ghosts " << nGhosts << " Clones " + << nClones; +} +// ------------------------------------------------------------------------- + + +// ----- Private method DivideHistos ----------------------------------- +void CbmTrackingTrdQa::DivideHistos(TH1* histo1, TH1* histo2, TH1* histo3) +{ + + if (!histo1 || !histo2 || !histo3) { + LOG(fatal) << GetName() << "::DivideHistos: " + << "NULL histogram pointer"; + } + + Int_t nBins = histo1->GetNbinsX(); + if (histo2->GetNbinsX() != nBins || histo3->GetNbinsX() != nBins) { + LOG(error) << GetName() << "::DivideHistos: " + << "Different bin numbers in histos"; + LOG(error) << histo1->GetName() << " " << histo1->GetNbinsX(); + LOG(error) << histo2->GetName() << " " << histo2->GetNbinsX(); + LOG(error) << histo3->GetName() << " " << histo3->GetNbinsX(); + return; + } + + Double_t c1, c2, c3, ce; + for (Int_t iBin = 0; iBin < nBins; iBin++) { + c1 = histo1->GetBinContent(iBin); + c2 = histo2->GetBinContent(iBin); + if (c2 != 0.) { + c3 = c1 / c2; + Double_t c4 = (c3 * (1. - c3) / c2); + if (c4 >= 0.) { ce = TMath::Sqrt(c3 * (1. - c3) / c2); } + else { + ce = 0; + } + } + else { + c3 = 0.; + ce = 0.; + } + histo3->SetBinContent(iBin, c3); + histo3->SetBinError(iBin, ce); + } +} +// ------------------------------------------------------------------------- + + +ClassImp(CbmTrackingTrdQa) diff --git a/reco/qa/CbmTrackingTrdQa.h b/reco/qa/CbmTrackingTrdQa.h new file mode 100644 index 0000000000..1497c378d9 --- /dev/null +++ b/reco/qa/CbmTrackingTrdQa.h @@ -0,0 +1,231 @@ +/* Copyright (C) 2006-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergey Gorbunov [committer] */ + +// ------------------------------------------------------------------------- +// ----- CbmTrackingTrdQa header file ----- +// ----- Created 15/08/22 ----- +// ------------------------------------------------------------------------- + + +/** CbmTrackingTrdQa.h + *@author S. Gorbunov + ** + ** Quality check task for Trd tracks + **/ + +#ifndef CbmTrackingTrdQa_H +#define CbmTrackingTrdQa_H 1 + +#include "CbmLink.h" +#include "CbmQaHist.h" + +#include <FairTask.h> + +#include "TH2F.h" +#include "TProfile2D.h" +#include <TFolder.h> +#include <TStopwatch.h> +#include <TVector3.h> + +#include "Riostream.h" + +class TH1F; +class CbmMCDataArray; +class CbmMCDataManager; +class CbmTimeSlice; + +class CbmTrackingTrdQa : public FairTask { + +public: + /** Default constructor **/ + CbmTrackingTrdQa(Int_t iVerbose = 1); + + /** Standard constructor + *@param minHits Minimal number of TrdHits for considered MCTracks + *@param quota True/all hits for track to be considered reconstructed + *@param iVerbose Verbosity level + **/ + CbmTrackingTrdQa(Int_t minHits, Double_t quota, Int_t iVerbose); + + /** Destructor **/ + virtual ~CbmTrackingTrdQa(); + + /** Set parameter containers **/ + virtual void SetParContainers(); + + /** Initialisation **/ + virtual InitStatus Init(); + + /** Reinitialisation **/ + virtual InitStatus ReInit(); + + /** Execution **/ + virtual void Exec(Option_t* opt); + +private: + /** Finish **/ + virtual void Finish(); + + /** Read the geometry parameters **/ + InitStatus GetGeometry(); + + /** Get the target node from the geometry **/ + void GetTargetPosition(); + + /** Create histograms **/ + void CreateHistos(); + + /** Reset histograms and counters **/ + void Reset(); + + // collect id's of MC events participating to the data + void CollectMcEvents(TClonesArray* Matches); + + /** Fill a map from MCTrack index to number of corresponding TrdHits **/ + void FillHitMap(); + + /** Fill a map from MCTrack index to matched GlobalTrack index + *@param nRec Number of reconstructed tracks (return) + *@param nGhosts Number of ghost tracks (return) + *@param nClones Number of clone tracks (return) + **/ + void FillTrackMatchMap(Int_t& nRec, Int_t& nGhosts, Int_t& nClones); + + /** Divide histograms (reco/all) with correct error for the efficiency + *@param histo1 reconstructed tracks + *@param histo2 all tracks (normalisation) + *@param histo3 efficiency + **/ + void DivideHistos(TH1* histo1, TH1* histo2, TH1* histo3); + + struct McTrackInfo { + std::map<Int_t, Int_t> fHitMap = {}; /// Trd station -> number of attached hits + Int_t fGlobalTrackMatch = -1; /// matched GlobalTrack index + Int_t fTrdTrackMatch = -1; /// matched TrdTrack index + Double_t fQuali = 0.; /// percentage of matched hits + Int_t fMatchedNHitsAll = 0; /// number of matched hits + Int_t fMatchedNHitsTrue = 0; /// number of matched hits + Bool_t fIsAccepted {0}; + Int_t fNtimesReconstructed {0}; + Bool_t fIsPrimary {0}; + Bool_t fIsFast {0}; + Bool_t fIsLong {0}; + Double_t fPt {0.}; + }; + + McTrackInfo& getMcTrackInfo(const CbmLink& link) + { + assert(fMcTrackInfoMap.find(link) != fMcTrackInfoMap.end()); + return fMcTrackInfoMap[link]; + } + + std::map<CbmLink, McTrackInfo> fMcTrackInfoMap = {}; //! map track link -> track info + + /// Setup + FairRootManager* fManager = nullptr; //! + CbmMCDataManager* fMcManager = nullptr; //! + CbmTimeSlice* fTimeSlice = nullptr; //! + + /// MC tracks + CbmMCDataArray* fMCTracks = nullptr; //! MCtrack + + /// Trd + Int_t fTrdNstations = 0; // Number of Trd stations + CbmMCDataArray* fTrdPoints = nullptr; //! TrdPoints + TClonesArray* fTrdHits = nullptr; //! TrdHits + TClonesArray* fTrdHitMatch = nullptr; //! TrdHitMatch + TClonesArray* fGlobalTracks = nullptr; //! GlobalTrack + TClonesArray* fGlobalTrackMatches = nullptr; //! GlobalTrackMatch + TClonesArray* fTrdTracks = nullptr; //! TrdTrack + TClonesArray* fTrdTrackMatches = nullptr; //! TrdTrackMatch + + /** Geometry parameters **/ + TVector3 fTargetPos = {0., 0., 0.}; // Target centre position + + /** Task parameters **/ + + Int_t fMinStations = 4; // Minimal number of stations with hits for considered MCTrack + + Double_t fQuota = 0.7; // True/all hits for track to be considered reconstructed + + TFolder fOutFolder = {"TrackingTrdQa", "TrackingTrdQa"}; /// output folder with histos and canvases + + /** Histograms **/ + + // eff. vs. XY + std::vector<CbmQaHist<TProfile2D>> fhStationEffXY; + + // eff. vs. pt, all + TH1F* fhPtAccAll = nullptr; + TH1F* fhPtRecAll = nullptr; + TH1F* fhPtEffAll = nullptr; + + // eff. vs. pt, vertex + TH1F* fhPtAccPrim = nullptr; + TH1F* fhPtRecPrim = nullptr; + TH1F* fhPtEffPrim = nullptr; + + // eff. vs. pt, non-vertex + TH1F* fhPtAccSec = nullptr; + TH1F* fhPtRecSec = nullptr; + TH1F* fhPtEffSec = nullptr; + + // eff. vs. np, all + TH1F* fhNpAccAll = nullptr; + TH1F* fhNpRecAll = nullptr; + TH1F* fhNpEffAll = nullptr; + + // eff. vs. np, vertex + TH1F* fhNpAccPrim = nullptr; + TH1F* fhNpRecPrim = nullptr; + TH1F* fhNpEffPrim = nullptr; + + // eff. vs. np, non-vertex + TH1F* fhNpAccSec = nullptr; + TH1F* fhNpRecSec = nullptr; + TH1F* fhNpEffSec = nullptr; + + // eff. vs. z, non-vertex + TH1F* fhZAccSec = nullptr; + TH1F* fhZRecSec = nullptr; + TH1F* fhZEffSec = nullptr; + + // # hits of clones and ghosts + TH1F* fhNhClones = nullptr; + TH1F* fhNhGhosts = nullptr; + + + /** List of histograms **/ + TList* fHistoList = nullptr; + + + /** Counters **/ + Int_t fNAll = 0; + Int_t fNAccAll = 0; + Int_t fNAccPrim = 0; + Int_t fNAccFast = 0; + Int_t fNAccFastLong = 0; + Int_t fNAccSec = 0; + Int_t fNRecAll = 0; + Int_t fNRecPrim = 0; + Int_t fNRecFast = 0; + Int_t fNRecFastLong = 0; + Int_t fNRecSec = 0; + Int_t fNGhosts = 0; + Int_t fNClones = 0; + Int_t fNEvents = 0; /** Number of events with success **/ + Int_t fNEventsFailed = 0; /** Number of events with failure **/ + Double_t fTime = 0.; /** Total real time used for good events **/ + + /** Timer **/ + TStopwatch fTimer = {}; + + CbmTrackingTrdQa(const CbmTrackingTrdQa&); + CbmTrackingTrdQa operator=(const CbmTrackingTrdQa&); + + ClassDef(CbmTrackingTrdQa, 0); +}; + + +#endif diff --git a/reco/qa/RecoQaLinkDef.h b/reco/qa/RecoQaLinkDef.h index 9975b3bee7..9568f48aca 100644 --- a/reco/qa/RecoQaLinkDef.h +++ b/reco/qa/RecoQaLinkDef.h @@ -11,5 +11,6 @@ #pragma link off all functions; #pragma link C++ class CbmRecoQa + ; +#pragma link C++ class CbmTrackingTrdQa + ; #endif -- GitLab