From 41bc6d47f006be43c75a6810b928bc71ae9e3c4e Mon Sep 17 00:00:00 2001 From: Matteo Concas Date: Thu, 5 Jan 2023 14:51:33 +0100 Subject: [PATCH] Add CheckVertices macro Finalise diagnostic macro --- .../ITSMFT/ITS/macros/test/CMakeLists.txt | 7 + .../ITSMFT/ITS/macros/test/CheckVertices.C | 321 ++++++++++++++++++ 2 files changed, 328 insertions(+) create mode 100644 Detectors/ITSMFT/ITS/macros/test/CheckVertices.C diff --git a/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt b/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt index cfd86b1a055b1..550e0373e4f44 100644 --- a/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt +++ b/Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt @@ -73,6 +73,13 @@ o2_add_test_root_macro(CheckSquasher.C O2::DataFormatsITSMFT LABELS its) +o2_add_test_root_macro(CheckVertices.C + PUBLIC_LINK_LIBRARIES O2::SimulationDataFormat + O2::ITSBase + O2::DataFormatsITS + O2::DataFormatsITSMFT + LABELS its) + o2_add_test_root_macro(DisplayTrack.C PUBLIC_LINK_LIBRARIES O2::ITSBase O2::DataFormatsITSMFT diff --git a/Detectors/ITSMFT/ITS/macros/test/CheckVertices.C b/Detectors/ITSMFT/ITS/macros/test/CheckVertices.C new file mode 100644 index 0000000000000..9d8b1bd49a6cf --- /dev/null +++ b/Detectors/ITSMFT/ITS/macros/test/CheckVertices.C @@ -0,0 +1,321 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#include + +#include +#include + +// #include "CommonUtils/RootSerializableKeyValueStore.h" +#include "Framework/Logger.h" +#include "ITSBase/GeometryTGeo.h" +#include "ReconstructionDataFormats/Vertex.h" +#include "SimulationDataFormat/MCEventHeader.h" +#include "SimulationDataFormat/TrackReference.h" +#include "SimulationDataFormat/MCTrack.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#include "DataFormatsITSMFT/CompCluster.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "DataFormatsITS/TrackITS.h" +#endif + +o2::MCCompLabel getMainLabel(std::vector& labs); + +struct ParticleInfo { + int event; + int pdg; + float pt; + float eta; + float phi; + int mother; + int first; + unsigned short clusters = 0u; + unsigned char isReco = 0u; + unsigned char isFake = 0u; + bool isPrimary = 0u; + unsigned char storedStatus = 2; /// not stored = 2, fake = 1, good = 0 + bool canContribToVertex = false; + std::array rofs = {-1, -1, -1, -1, -1, -1, -1}; /// readout frames of corresponding clusters + o2::its::TrackITS track; + o2::MCCompLabel lab; +}; + +struct RofInfo { + void print(); + void uniqeff(); + int id = 0; + std::vector eventIds; // ID of events in rof + std::vector usedIds; // EvtID used to calculate actual efficiency + std::vector parts; // Particle usable for vertexing + std::vector> vertLabels; // Labels associated to contributors to vertex + std::unordered_map> simVerts; // Simulated vertices of events that can be spot in current rof + std::vector>> recoVerts; // Vertices found in current ROF + float recoeff; // Vertexing efficiency +}; + +void RofInfo::print() +{ + std::cout << "\n=================================== ROF " << id << " ============================================ \n"; + // Simulated vertices + for (auto& sV : simVerts) { + std::cout << "\tSimulated vertex for event: " << sV.first << " vertex:" + << " x= " << sV.second[0] + << " y= " << sV.second[1] + << " z= " << sV.second[2] + << std::endl; + std::cout << "\t\tPotentially contributing tracks:\n"; + for (auto& part : parts) { + if (part.lab.getEventID() == sV.first && part.canContribToVertex) { + std::cout << "\t\t\t" << part.lab << "\t" << part.pt << " [GeV]\t" << part.pdg << std::endl; + } + } + std::cout << std::endl; + } + + // Reconstructed vertices + for (size_t iV{0}; iV < recoVerts.size(); ++iV) { + std::cout << "\tReconstructed vertex for event: " << getMainLabel(vertLabels[iV]).getEventID() << ":" + << " x= " << recoVerts[iV].getX() + << " y= " << recoVerts[iV].getY() + << " z= " << recoVerts[iV].getZ() + << std::endl; + std::cout << "\t\tContributor labels:\n"; + for (auto& l : vertLabels[iV]) { + std::cout << "\t\t\t" << l << std::endl; + } + } + + // Efficiency + if (simVerts.size() || recoVerts.size()) { + std::cout << "\n\tEfficiency: " << recoeff * 100 << " %\n"; + } +} + +void RofInfo::uniqeff() +{ + auto c{0}; + int current{-42}; + std::sort(parts.begin(), parts.end(), [](ParticleInfo& lp, ParticleInfo& rp) { return lp.lab.getEventID() > rp.lab.getEventID(); }); // sorting at this point should be harmless. + for (auto& p : parts) { + if (p.lab.getEventID() != current) { + eventIds.push_back(p.lab.getEventID()); + current = p.lab.getEventID(); + } + } + + usedIds.resize(eventIds.size(), false); + for (size_t iV{0}; iV < vertLabels.size(); ++iV) { + auto label = getMainLabel(vertLabels[iV]); + for (size_t evId{0}; evId < eventIds.size(); ++evId) { + if (eventIds[evId] == label.getEventID() && !usedIds[evId]) { + usedIds[evId] = true; + ++c; + } + } + } + recoeff = (float)c / (float)eventIds.size(); +} + +#pragma link C++ class ParticleInfo + ; +#pragma link C++ class RofInfo + ; + +void CheckVertices(const int dumprof = -1, std::string path = "tf1/", std::string tracfile = "o2trac_its.root", std::string clusfile = "o2clus_its.root", std::string kinefile = "sgn_1_Kine.root") +{ + using Vertex = o2::dataformats::Vertex>; + using namespace o2::dataformats; + using namespace o2::itsmft; + using namespace o2::its; + + // Geometry + o2::base::GeometryManager::loadGeometry(path.data()); + auto gman = o2::its::GeometryTGeo::Instance(); + + // MC tracks and event header + TFile* file0 = TFile::Open((path + kinefile).data()); + TTree* mcTree = (TTree*)gFile->Get("o2sim"); + mcTree->SetBranchStatus("*", 0); // disable all branches + mcTree->SetBranchStatus("MCEventHeader*", 1); + mcTree->SetBranchStatus("MCTrack*", 1); + + std::vector* mcArr = nullptr; + mcTree->SetBranchAddress("MCTrack", &mcArr); + MCEventHeader* eventHeader = nullptr; + mcTree->SetBranchAddress("MCEventHeader.", &eventHeader); + + // Clusters + TFile::Open((path + clusfile).data()); + TTree* clusTree = (TTree*)gFile->Get("o2sim"); + std::vector* clusArr = nullptr; + clusTree->SetBranchAddress("ITSClusterComp", &clusArr); + std::vector* clusROFRecords = nullptr; + clusTree->SetBranchAddress("ITSClustersROF", &clusROFRecords); + + // Cluster MC labels + o2::dataformats::MCTruthContainer* clusLabArr = nullptr; + clusTree->SetBranchAddress("ITSClusterMCTruth", &clusLabArr); + + // Reconstructed vertices + TFile* recFile = TFile::Open((path + tracfile).data()); + TTree* recTree = (TTree*)recFile->Get("o2sim"); + + std::vector* recVerArr = nullptr; + recTree->SetBranchAddress("Vertices", &recVerArr); + std::vector* recVerROFArr = nullptr; + recTree->SetBranchAddress("VerticesROF", &recVerROFArr); + std::vector* recLabelsArr = nullptr; + recTree->SetBranchAddress("ITSVertexMCTruth", &recLabelsArr); + + // Process + // Fill MC info + auto nev{mcTree->GetEntriesFast()}; + std::vector> info(nev); + std::vector> simVerts; + for (auto n{0}; n < nev; ++n) { + mcTree->GetEvent(n); + info[n].resize(mcArr->size()); + // Event header + for (unsigned int mcI{0}; mcI < mcArr->size(); ++mcI) { + auto part = mcArr->at(mcI); + info[n][mcI].event = n; + info[n][mcI].pdg = part.GetPdgCode(); + info[n][mcI].pt = part.GetPt(); + info[n][mcI].phi = part.GetPhi(); + info[n][mcI].eta = part.GetEta(); + info[n][mcI].isPrimary = part.isPrimary(); + } + simVerts.push_back({eventHeader->GetX(), eventHeader->GetY(), eventHeader->GetZ()}); + } + + // Fill ROF info and complement MC info with cluster info + std::vector rofinfo; + for (int frame = 0; frame < clusTree->GetEntriesFast(); frame++) { // Cluster frames + if (!clusTree->GetEvent(frame)) + continue; + rofinfo.resize(clusROFRecords->size()); + for (size_t rof{0}; rof < clusROFRecords->size(); ++rof) { + for (int iClus{clusROFRecords->at(rof).getFirstEntry()}; iClus < clusROFRecords->at(rof).getFirstEntry() + clusROFRecords->at(rof).getNEntries(); ++iClus) { + auto lab = (clusLabArr->getLabels(iClus))[0]; + if (!lab.isValid() || lab.getSourceID() != 0 || !lab.isCorrect()) + continue; + + int trackID, evID, srcID; + bool fake; + lab.get(trackID, evID, srcID, fake); + if (evID < 0 || evID >= (int)info.size()) { + std::cout << "Cluster MC label eventID out of range" << std::endl; + continue; + } + if (trackID < 0 || trackID >= (int)info[evID].size()) { + std::cout << "Cluster MC label trackID out of range" << std::endl; + continue; + } + info[evID][trackID].lab = lab; // seems redundant but we are going to copy these info and loosing the nice evt/tr_id ordering + const CompClusterExt& c = (*clusArr)[iClus]; + auto layer = gman->getLayer(c.getSensorID()); + info[evID][trackID].clusters |= 1 << layer; + info[evID][trackID].rofs[layer] = rof; + } + } + } + + for (size_t evt{0}; evt < info.size(); ++evt) { + auto& evInfo = info[evt]; + int ntrackable{0}; + int nusable{0}; + for (auto& part : evInfo) { + if (part.clusters & (1 << 0) && part.clusters & (1 << 1) && part.clusters & (1 << 2)) { + ++ntrackable; + if (part.rofs[0] > -1 && part.rofs[0] == part.rofs[1] && part.rofs[1] == part.rofs[2]) { + ++nusable; + part.canContribToVertex = true; + rofinfo[part.rofs[0]].parts.push_back(part); + int trackID, evID, srcID; + bool fake; + part.lab.get(trackID, evID, srcID, fake); + rofinfo[part.rofs[0]].simVerts[evID] = simVerts[evID]; + } + } + } + } + + // Reco vertices processing + for (int frame = 0; frame < recTree->GetEntriesFast(); frame++) { // Vertices frames + if (!recTree->GetEvent(frame)) { + continue; + } + // loop on rof records + int contLabIdx{0}; + for (size_t iRecord{0}; iRecord < recVerROFArr->size(); ++iRecord) { + auto& rec = recVerROFArr->at(iRecord); + auto verStartIdx = rec.getFirstEntry(), verSize = rec.getNEntries(); + int totContrib{0}, nVerts{0}; + rofinfo[iRecord].id = iRecord; + rofinfo[iRecord].vertLabels.resize(verSize); + int vertCounter{0}; + for (int iVertex{verStartIdx}; iVertex < verStartIdx + verSize; ++iVertex, ++vertCounter) { + auto vert = recVerArr->at(iVertex); + rofinfo[iRecord].recoVerts.push_back(vert); + totContrib += vert.getNContributors(); + nVerts += 1; + for (int ic{0}; ic < vert.getNContributors(); ++ic, ++contLabIdx) { + rofinfo[iRecord].vertLabels[vertCounter].push_back(recLabelsArr->at(contLabIdx)); + // std::cout << "Pushed " << rofinfo[iRecord].vertLabels[vertCounter].back() << " at position " << rofinfo[iRecord].vertLabels[vertCounter].size() << std::endl; + } + } + } + } + // Epilog + LOGP(info, "ROF inspection summary"); + size_t nvt{0}, nevts{0}; + if (dumprof < 0) { + for (size_t iROF{0}; iROF < rofinfo.size(); ++iROF) { + auto& rof = rofinfo[iROF]; + nvt += rof.recoVerts.size(); + nevts += rof.simVerts.size(); + rof.uniqeff(); + rof.print(); + } + } else { + rofinfo[dumprof].print(); + rofinfo[dumprof].uniqeff(); + nvt += rofinfo[dumprof].recoVerts.size(); + nevts += rofinfo[dumprof].simVerts.size(); + } + LOGP(info, "Summary:"); + LOGP(info, "\nFound {} vertices in {} usable events out of {}", nvt, nevts, simVerts.size()); +} + +o2::MCCompLabel getMainLabel(std::vector& labs) +{ + o2::MCCompLabel lab; + size_t max_count = 0; + for (size_t i = 0; i < labs.size(); i++) { + size_t count = 1; + for (size_t j = i + 1; j < labs.size(); j++) + if (labs[i] == labs[j]) + count++; + if (count > max_count) + max_count = count; + } + + for (size_t i = 0; i < labs.size(); i++) { + size_t count = 1; + for (size_t j = i + 1; j < labs.size(); j++) + if (labs[i] == labs[j]) + count++; + if (count == max_count) + lab = labs[i]; + } + return lab; +} \ No newline at end of file