Skip to content

Commit 46648d6

Browse files
Lucas-NLshahor02
Lucas-NL
authored and
shahor02
committed
Adapted CheckTopologies for MFT
1 parent 3dd2eb9 commit 46648d6

3 files changed

Lines changed: 312 additions & 0 deletions

File tree

‎Detectors/ITSMFT/MFT/macros/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,3 +9,4 @@
99
# submit itself to any jurisdiction.
1010

1111
add_subdirectory(mapping)
12+
add_subdirectory(test)

‎Detectors/ITSMFT/MFT/macros/test/CMakeLists.txt‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
# Copyright CERN and copyright holders of ALICE O2. This software is distributed
2+
# under the terms of the GNU General Public License v3 (GPL Version 3), copied
3+
# verbatim in the file "COPYING".
4+
#
5+
# See http://alice-o2.web.cern.ch/license for full licensing information.
6+
#
7+
# In applying this license CERN does not waive the privileges and immunities
8+
# granted to it by virtue of its status as an Intergovernmental Organization or
9+
# submit itself to any jurisdiction.
10+
11+
o2_add_test_root_macro(CheckTopologies.C
12+
PUBLIC_LINK_LIBRARIES O2::MathUtils
13+
O2::MFTBase
14+
O2::ITSMFTReconstruction
15+
O2::ITSMFTSimulation
16+
O2::DataFormatsITSMFT
17+
O2::SimulationDataFormat
18+
LABELS mft)
19+

‎Detectors/ITSMFT/MFT/macros/test/CheckTopologies.C‎

Lines changed: 292 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,292 @@
1+
/// \file CheckTopologies.C
2+
/// Macros to test the generation of a dictionary of topologies. Three dictionaries are generated: one with signal-cluster only, one with noise-clusters only and one with all the clusters.
3+
4+
#if !defined(__CLING__) || defined(__ROOTCLING__)
5+
#include <TAxis.h>
6+
#include <TCanvas.h>
7+
#include <TSystem.h>
8+
#include <TFile.h>
9+
#include <TH1F.h>
10+
#include <TH2F.h>
11+
#include <TNtuple.h>
12+
#include <TString.h>
13+
#include <TStyle.h>
14+
#include <TTree.h>
15+
#include <TStopwatch.h>
16+
#include <fstream>
17+
#include <string>
18+
19+
#include "MathUtils/Utils.h"
20+
#include "MFTBase/GeometryTGeo.h"
21+
#include "ITSMFTReconstruction/BuildTopologyDictionary.h"
22+
#include "DataFormatsITSMFT/CompCluster.h"
23+
#include "DataFormatsITSMFT/ClusterTopology.h"
24+
#include "DataFormatsITSMFT/ROFRecord.h"
25+
#include "ITSMFTSimulation/Hit.h"
26+
#include "MathUtils/Cartesian3D.h"
27+
#include "SimulationDataFormat/MCCompLabel.h"
28+
#include "SimulationDataFormat/MCTruthContainer.h"
29+
#include "DetectorsCommonDataFormats/NameConf.h"
30+
#include "SimulationDataFormat/DigitizationContext.h"
31+
#include "Framework/Logger.h"
32+
#include <unordered_map>
33+
#endif
34+
35+
void CheckTopologies(std::string clusfile = "mftclusters.root",
36+
std::string hitfile = "o2sim_HitsMFT.root",
37+
std::string collContextfile = "collisioncontext.root",
38+
std::string inputGeom = "")
39+
{
40+
const int QEDSourceID = 99; // Clusters from this MC source correspond to QED electrons
41+
42+
using namespace o2::base;
43+
using namespace o2::mft;
44+
45+
using o2::itsmft::BuildTopologyDictionary;
46+
using o2::itsmft::ClusterTopology;
47+
using o2::itsmft::CompCluster;
48+
using o2::itsmft::CompClusterExt;
49+
using o2::itsmft::Hit;
50+
using ROFRec = o2::itsmft::ROFRecord;
51+
using MC2ROF = o2::itsmft::MC2ROFRecord;
52+
using HitVec = std::vector<Hit>;
53+
using MC2HITS_map = std::unordered_map<uint64_t, int>; // maps (track_ID<<16 + chip_ID) to entry in the hit vector
54+
std::unordered_map<int, int> hadronicMCMap; // mapping from MC event entry to hadronic event ID
55+
std::vector<HitVec*> hitVecPool;
56+
std::vector<MC2HITS_map> mc2hitVec;
57+
const o2::steer::DigitizationContext* digContext = nullptr;
58+
TStopwatch sw;
59+
sw.Start();
60+
61+
// Geometry
62+
o2::base::GeometryManager::loadGeometry(inputGeom);
63+
auto gman = o2::mft::GeometryTGeo::Instance();
64+
gman->fillMatrixCache(o2::utils::bit2Mask(o2::TransformType::T2L, o2::TransformType::T2GRot,
65+
o2::TransformType::L2G)); // request cached transforms
66+
67+
// Hits
68+
TFile* fileH = nullptr;
69+
TTree* hitTree = nullptr;
70+
71+
if (!hitfile.empty() && !collContextfile.empty() && !gSystem->AccessPathName(hitfile.c_str()) && !gSystem->AccessPathName(collContextfile.c_str())) {
72+
fileH = TFile::Open(hitfile.data());
73+
hitTree = (TTree*)fileH->Get("o2sim");
74+
mc2hitVec.resize(hitTree->GetEntries());
75+
hitVecPool.resize(hitTree->GetEntries(), nullptr);
76+
digContext = o2::steer::DigitizationContext::loadFromFile(collContextfile);
77+
78+
auto& intGlo = digContext->getEventParts(digContext->isQEDProvided());
79+
int hadrID = -1, nGlo = intGlo.size(), nHadro = 0;
80+
for (int iglo = 0; iglo < nGlo; iglo++) {
81+
const auto& parts = intGlo[iglo];
82+
bool found = false;
83+
for (auto& part : parts) {
84+
if (part.sourceID == 0) { // we use underlying background
85+
hadronicMCMap[iglo] = part.entryID;
86+
found = true;
87+
nHadro++;
88+
break;
89+
}
90+
}
91+
if (!found) {
92+
hadronicMCMap[iglo] = -1;
93+
}
94+
}
95+
if (nHadro < hitTree->GetEntries()) {
96+
LOG(FATAL) << "N=" << nHadro << " hadronic events < "
97+
<< " N=" << hitTree->GetEntries() << " Hit enties.";
98+
}
99+
}
100+
101+
// Clusters
102+
TFile* fileCl = TFile::Open(clusfile.data());
103+
TTree* clusTree = (TTree*)fileCl->Get("o2sim");
104+
std::vector<CompClusterExt>* clusArr = nullptr;
105+
clusTree->SetBranchAddress("MFTClusterComp", &clusArr);
106+
std::vector<unsigned char>* patternsPtr = nullptr;
107+
auto pattBranch = clusTree->GetBranch("MFTClusterPatt");
108+
if (pattBranch) {
109+
pattBranch->SetAddress(&patternsPtr);
110+
}
111+
112+
// ROFrecords
113+
std::vector<ROFRec> rofRecVec, *rofRecVecP = &rofRecVec;
114+
clusTree->SetBranchAddress("MFTClustersROF", &rofRecVecP);
115+
116+
// Cluster MC labels
117+
o2::dataformats::MCTruthContainer<o2::MCCompLabel>* clusLabArr = nullptr;
118+
std::vector<MC2ROF> mc2rofVec, *mc2rofVecP = &mc2rofVec;
119+
if (hitTree && clusTree->GetBranch("MFTClusterMCTruth")) {
120+
clusTree->SetBranchAddress("MFTClusterMCTruth", &clusLabArr);
121+
clusTree->SetBranchAddress("MFTClustersMC2ROF", &mc2rofVecP);
122+
}
123+
clusTree->GetEntry(0);
124+
125+
// Topologies dictionaries: 1) all clusters 2) signal clusters only 3) noise clusters only
126+
BuildTopologyDictionary completeDictionary;
127+
BuildTopologyDictionary signalDictionary;
128+
BuildTopologyDictionary noiseDictionary;
129+
130+
int nROFRec = (int)rofRecVec.size();
131+
std::vector<int> mcEvMin, mcEvMax;
132+
133+
if (clusLabArr) { // >> build min and max MC events used by each ROF
134+
mcEvMin.resize(nROFRec, hitTree->GetEntries());
135+
mcEvMax.resize(nROFRec, -1);
136+
for (int imc = mc2rofVec.size(); imc--;) {
137+
int hadrID = hadronicMCMap[imc];
138+
if (hadrID < 0) {
139+
continue;
140+
}
141+
const auto& mc2rof = mc2rofVec[imc];
142+
if (mc2rof.rofRecordID < 0) {
143+
continue; // this MC event did not contribute to any ROF
144+
}
145+
for (int irfd = mc2rof.maxROF - mc2rof.minROF + 1; irfd--;) {
146+
int irof = mc2rof.rofRecordID + irfd;
147+
if (mcEvMin[irof] > hadrID) {
148+
mcEvMin[irof] = hadrID;
149+
}
150+
if (mcEvMax[irof] < hadrID) {
151+
mcEvMax[irof] = hadrID;
152+
}
153+
}
154+
}
155+
} // << build min and max MC events used by each ROF
156+
157+
auto pattIdx = patternsPtr->cbegin();
158+
for (int irof = 0; irof < nROFRec; irof++) {
159+
const auto& rofRec = rofRecVec[irof];
160+
161+
rofRec.print();
162+
163+
if (clusLabArr) { // >> read and map MC events contributing to this ROF
164+
for (int im = mcEvMin[irof]; im <= mcEvMax[irof]; im++) {
165+
if (!hitVecPool[im]) {
166+
hitTree->SetBranchAddress("MFTHit", &hitVecPool[im]);
167+
hitTree->GetEntry(im);
168+
auto& mc2hit = mc2hitVec[im];
169+
const auto* hitArray = hitVecPool[im];
170+
for (int ih = hitArray->size(); ih--;) {
171+
const auto& hit = (*hitArray)[ih];
172+
uint64_t key = (uint64_t(hit.GetTrackID()) << 32) + hit.GetDetectorID();
173+
mc2hit.emplace(key, ih);
174+
}
175+
}
176+
}
177+
} // << cache MC events contributing to this ROF
178+
179+
for (int icl = 0; icl < rofRec.getNEntries(); icl++) {
180+
int clEntry = rofRec.getFirstEntry() + icl; // entry of icl-th cluster of this ROF in the vector of clusters
181+
// do we read MC data?
182+
183+
const auto& cluster = (*clusArr)[clEntry];
184+
185+
if (cluster.getPatternID() != CompCluster::InvalidPatternID) {
186+
LOG(WARNING) << "Encountered patternID = " << cluster.getPatternID() << " != " << CompCluster::InvalidPatternID;
187+
LOG(WARNING) << "Clusters have already been generated with a dictionary! Quitting";
188+
return;
189+
}
190+
191+
ClusterTopology topology;
192+
o2::itsmft::ClusterPattern pattern(pattIdx);
193+
topology.setPattern(pattern);
194+
195+
float dX = BuildTopologyDictionary::IgnoreVal, dZ = BuildTopologyDictionary::IgnoreVal;
196+
if (clusLabArr) {
197+
const auto& lab = (clusLabArr->getLabels(clEntry))[0];
198+
auto srcID = lab.getSourceID();
199+
if (lab.isValid() && srcID != QEDSourceID) { // use MC truth info only for non-QED and non-noise clusters
200+
auto trID = lab.getTrackID();
201+
const auto& mc2hit = mc2hitVec[lab.getEventID()];
202+
const auto* hitArray = hitVecPool[lab.getEventID()];
203+
Int_t chipID = cluster.getSensorID();
204+
uint64_t key = (uint64_t(trID) << 32) + chipID;
205+
auto hitEntry = mc2hit.find(key);
206+
if (hitEntry != mc2hit.end()) {
207+
const auto& hit = (*hitArray)[hitEntry->second];
208+
auto locH = gman->getMatrixL2G(chipID) ^ (hit.GetPos()); // inverse conversion from global to local
209+
auto locHsta = gman->getMatrixL2G(chipID) ^ (hit.GetPosStart());
210+
locH.SetXYZ(0.5 * (locH.X() + locHsta.X()), 0.5 * (locH.Y() + locHsta.Y()), 0.5 * (locH.Z() + locHsta.Z()));
211+
const auto locC = o2::itsmft::TopologyDictionary::getClusterCoordinates(cluster, pattern);
212+
dX = locH.X() - locC.X();
213+
dZ = locH.Z() - locC.Z();
214+
} else {
215+
printf("Failed to find MC hit entry for Tr:%d chipID:%d\n", trID, chipID);
216+
lab.print();
217+
}
218+
signalDictionary.accountTopology(topology, dX, dZ);
219+
} else {
220+
noiseDictionary.accountTopology(topology, dX, dZ);
221+
}
222+
}
223+
completeDictionary.accountTopology(topology, dX, dZ);
224+
}
225+
// clean MC cache for events which are not needed anymore
226+
if (clusLabArr) {
227+
int irfNext = irof;
228+
int limMC = irfNext == nROFRec ? hitVecPool.size() : mcEvMin[irfNext]; // can clean events up to this
229+
for (int imc = mcEvMin[irof]; imc < limMC; imc++) {
230+
delete hitVecPool[imc];
231+
hitVecPool[imc] = nullptr;
232+
mc2hitVec[imc].clear();
233+
}
234+
}
235+
}
236+
auto dID = o2::detectors::DetID::MFT;
237+
238+
completeDictionary.setThreshold(0.0001);
239+
completeDictionary.groupRareTopologies();
240+
completeDictionary.printDictionaryBinary(o2::base::NameConf::getDictionaryFileName(dID, "", ".bin"));
241+
completeDictionary.printDictionary(o2::base::NameConf::getDictionaryFileName(dID, "", ".txt"));
242+
completeDictionary.saveDictionaryRoot(o2::base::NameConf::getDictionaryFileName(dID, "", ".root"));
243+
244+
TFile histogramOutput("histograms.root", "recreate");
245+
TCanvas* cComplete = new TCanvas("cComplete", "Distribution of all the topologies");
246+
cComplete->cd();
247+
cComplete->SetLogy();
248+
TH1F* hComplete = nullptr;
249+
o2::itsmft::TopologyDictionary::getTopologyDistribution(completeDictionary.getDictionary(), hComplete, "hComplete");
250+
hComplete->SetDirectory(0);
251+
hComplete->Draw("hist");
252+
hComplete->Write();
253+
cComplete->Write();
254+
255+
TCanvas* cNoise = nullptr;
256+
TCanvas* cSignal = nullptr;
257+
TH1F* hNoise = nullptr;
258+
TH1F* hSignal = nullptr;
259+
260+
if (clusLabArr) {
261+
noiseDictionary.setThreshold(0.0001);
262+
noiseDictionary.groupRareTopologies();
263+
noiseDictionary.printDictionaryBinary(o2::base::NameConf::getDictionaryFileName(dID, "noise", ".bin"));
264+
noiseDictionary.printDictionary(o2::base::NameConf::getDictionaryFileName(dID, "noise", ".txt"));
265+
noiseDictionary.saveDictionaryRoot(o2::base::NameConf::getDictionaryFileName(dID, "noise", ".root"));
266+
signalDictionary.setThreshold(0.0001);
267+
signalDictionary.groupRareTopologies();
268+
signalDictionary.printDictionaryBinary(o2::base::NameConf::getDictionaryFileName(dID, "signal", ".bin"));
269+
signalDictionary.printDictionary(o2::base::NameConf::getDictionaryFileName(dID, "signal", ".txt"));
270+
signalDictionary.saveDictionaryRoot(o2::base::NameConf::getDictionaryFileName(dID, "signal", ".root"));
271+
cNoise = new TCanvas("cNoise", "Distribution of noise topologies");
272+
cNoise->cd();
273+
cNoise->SetLogy();
274+
o2::itsmft::TopologyDictionary::getTopologyDistribution(noiseDictionary.getDictionary(), hNoise, "hNoise");
275+
hNoise->SetDirectory(0);
276+
hNoise->Draw("hist");
277+
histogramOutput.cd();
278+
hNoise->Write();
279+
cNoise->Write();
280+
cSignal = new TCanvas("cSignal", "cSignal");
281+
cSignal->cd();
282+
cSignal->SetLogy();
283+
o2::itsmft::TopologyDictionary::getTopologyDistribution(signalDictionary.getDictionary(), hSignal, "hSignal");
284+
hSignal->SetDirectory(0);
285+
hSignal->Draw("hist");
286+
histogramOutput.cd();
287+
hSignal->Write();
288+
cSignal->Write();
289+
sw.Stop();
290+
sw.Print();
291+
}
292+
}

0 commit comments

Comments
 (0)