Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ITS: add macro to check seeds #10555

Merged
merged 1 commit into from
Jan 15, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions Detectors/ITSMFT/ITS/macros/test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
321 changes: 321 additions & 0 deletions Detectors/ITSMFT/ITS/macros/test/CheckVertices.C
Original file line number Diff line number Diff line change
@@ -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 <vector>
#include <unordered_map>

#include <TFile.h>
#include <TTree.h>

// #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<o2::MCCompLabel>& 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<int, 7> 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<int> eventIds; // ID of events in rof
std::vector<bool> usedIds; // EvtID used to calculate actual efficiency
std::vector<ParticleInfo> parts; // Particle usable for vertexing
std::vector<std::vector<o2::MCCompLabel>> vertLabels; // Labels associated to contributors to vertex
std::unordered_map<int, std::array<double, 3>> simVerts; // Simulated vertices of events that can be spot in current rof <evtId, pos[3]>
std::vector<o2::dataformats::Vertex<o2::dataformats::TimeStamp<int>>> 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<o2::dataformats::TimeStamp<int>>;
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<o2::MCTrack>* 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<CompClusterExt>* clusArr = nullptr;
clusTree->SetBranchAddress("ITSClusterComp", &clusArr);
std::vector<o2::itsmft::ROFRecord>* clusROFRecords = nullptr;
clusTree->SetBranchAddress("ITSClustersROF", &clusROFRecords);

// Cluster MC labels
o2::dataformats::MCTruthContainer<o2::MCCompLabel>* clusLabArr = nullptr;
clusTree->SetBranchAddress("ITSClusterMCTruth", &clusLabArr);

// Reconstructed vertices
TFile* recFile = TFile::Open((path + tracfile).data());
TTree* recTree = (TTree*)recFile->Get("o2sim");

std::vector<Vertex>* recVerArr = nullptr;
recTree->SetBranchAddress("Vertices", &recVerArr);
std::vector<ROFRecord>* recVerROFArr = nullptr;
recTree->SetBranchAddress("VerticesROF", &recVerROFArr);
std::vector<o2::MCCompLabel>* recLabelsArr = nullptr;
recTree->SetBranchAddress("ITSVertexMCTruth", &recLabelsArr);

// Process
// Fill MC info
auto nev{mcTree->GetEntriesFast()};
std::vector<std::vector<ParticleInfo>> info(nev);
std::vector<std::array<double, 3>> 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> 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<o2::MCCompLabel>& 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;
}