Skip to content

Commit 8cc7f07

Browse files
f3schf3sch
committed
ITS3: fix CheckClustersITS3 macro
Signed-off-by: Felix Schlepper <felix.schlepper@cern.ch>
1 parent 0318b9f commit 8cc7f07

1 file changed

Lines changed: 62 additions & 116 deletions

File tree

‎Detectors/Upgrades/ITS3/macros/test/CheckClustersITS3.C‎

Lines changed: 62 additions & 116 deletions
Original file line numberDiff line numberDiff line change
@@ -36,14 +36,15 @@
3636
#include "MathUtils/Utils.h"
3737
#include "SimulationDataFormat/MCCompLabel.h"
3838
#include "SimulationDataFormat/MCTruthContainer.h"
39+
#include "Steer/MCKinematicsReader.h"
3940
#include "DetectorsCommonDataFormats/DetectorNameConf.h"
4041
#endif
4142

4243
void CheckClustersITS3(const std::string& clusfile = "o2clus_its.root",
4344
const std::string& hitfile = "o2sim_HitsIT3.root",
4445
const std::string& inputGeom = "",
45-
std::string dictfile = "./ccdb/IT3/Calib/ClusterDictionary/snapshot.root",
46-
bool batch = false)
46+
std::string dictfile = "../ccdb/IT3/Calib/ClusterDictionary/snapshot.root",
47+
bool batch = true)
4748
{
4849
gROOT->SetBatch(batch);
4950

@@ -57,17 +58,15 @@ void CheckClustersITS3(const std::string& clusfile = "o2clus_its.root",
5758
using ROFRec = o2::itsmft::ROFRecord;
5859
using MC2ROF = o2::itsmft::MC2ROFRecord;
5960
using HitVec = std::vector<Hit>;
60-
using MC2HITS_map = std::unordered_map<uint64_t, int>; // maps (track_ID<<32 + chip_ID) to entry in the hit vector
61+
using MC2HITS_map = std::unordered_map<uint64_t, std::vector<int>>; // maps (track_ID<<32 + chip_ID) to entry in the hit vector
6162
std::array<MosaixSegmentation, 3> mMosaixSegmentations{0, 1, 2};
6263

63-
std::vector<HitVec*> hitVecPool;
64-
std::vector<MC2HITS_map> mc2hitVec;
65-
64+
ULong_t cClustersIB{0}, cClustersOB{0};
6665
ULong_t cPattValidIB{0}, cPattInvalidIB{0}, cLabelInvalidIB{0}, cNoMCIB{0};
6766
ULong_t cPattValidOB{0}, cPattInvalidOB{0}, cLabelInvalidOB{0}, cNoMCOB{0};
6867

6968
TFile fout("CheckClusters.root", "recreate");
70-
TNtuple nt("ntc", "cluster ntuple", "ev:lab:hlx:hlz:hgx:hgz:tx:tz:cgx:cgy:cgz:clx:cly:clz:dx:dy:dz:ex:ez:patid:rof:npx:id:eta:row:col:lay");
69+
TNtuple nt("ntc", "cluster ntuple", "ev:lab:hlx:hlz:hgx:hgz:tx:tz:cgx:cgy:cgz:clx:cly:clz:dx:dy:dz:ex:ez:patid:rof:npx:id:eta:row:col:lay:prim");
7170

7271
// Geometry
7372
o2::base::GeometryManager::loadGeometry(inputGeom);
@@ -80,6 +79,8 @@ void CheckClustersITS3(const std::string& clusfile = "o2clus_its.root",
8079
auto* hitTree = dynamic_cast<TTree*>(fileH.Get("o2sim"));
8180
std::vector<o2::itsmft::Hit>* hitArray = nullptr;
8281
hitTree->SetBranchAddress("IT3Hit", &hitArray);
82+
std::vector<HitVec*> hitVecPool;
83+
std::vector<MC2HITS_map> mc2hitVec;
8384
mc2hitVec.resize(hitTree->GetEntries());
8485
hitVecPool.resize(hitTree->GetEntries(), nullptr);
8586

@@ -106,65 +107,41 @@ void CheckClustersITS3(const std::string& clusfile = "o2clus_its.root",
106107
}
107108
dict.print();
108109

110+
//
111+
o2::steer::MCKinematicsReader reader("collisioncontext.root");
112+
109113
// ROFrecords
110114
std::vector<ROFRec> rofRecVec, *rofRecVecP = &rofRecVec;
111115
clusTree->SetBranchAddress("ITSClustersROF", &rofRecVecP);
112116

113117
// Cluster MC labels
114118
o2::dataformats::MCTruthContainer<o2::MCCompLabel>* clusLabArr = nullptr;
115119
std::vector<MC2ROF> mc2rofVec, *mc2rofVecP = &mc2rofVec;
116-
if ((hitTree != nullptr) && (clusTree->GetBranch("ITSClusterMCTruth") != nullptr)) {
117-
clusTree->SetBranchAddress("ITSClusterMCTruth", &clusLabArr);
118-
clusTree->SetBranchAddress("ITSClustersMC2ROF", &mc2rofVecP);
119-
}
120-
120+
clusTree->SetBranchAddress("ITSClusterMCTruth", &clusLabArr);
121121
clusTree->GetEntry(0);
122+
auto pattIt = patternsPtr->cbegin();
122123
unsigned int nROFRec = (int)rofRecVec.size();
123-
std::vector<int> mcEvMin(nROFRec, hitTree->GetEntries());
124-
std::vector<int> mcEvMax(nROFRec, -1);
125124

126-
// >> build min and max MC events used by each ROF
127-
for (int imc = mc2rofVec.size(); imc--;) {
128-
const auto& mc2rof = mc2rofVec[imc];
129-
// printf("MCRecord: ");
130-
// mc2rof.print();
131-
if (mc2rof.rofRecordID < 0) {
132-
continue; // this MC event did not contribute to any ROF
133-
}
134-
for (unsigned int irfd = mc2rof.maxROF - mc2rof.minROF + 1; irfd--;) {
135-
unsigned int irof = mc2rof.rofRecordID + irfd;
136-
if (irof >= nROFRec) {
137-
LOG(error) << "ROF=" << irof << " from MC2ROF record is >= N ROFs=" << nROFRec;
125+
// >> load all MC hits upfront
126+
for (int im = 0; im < (int)hitTree->GetEntries(); im++) {
127+
if (hitVecPool[im] == nullptr) {
128+
hitTree->SetBranchAddress("IT3Hit", &hitVecPool[im]);
129+
if (hitTree->GetEntry(im) <= 0 || hitVecPool[im] == nullptr) {
130+
LOGP(error, "Cannot read IT3Hit entry {} from {}", im, hitfile);
131+
return;
138132
}
139-
if (mcEvMin[irof] > imc) {
140-
mcEvMin[irof] = imc;
141-
}
142-
if (mcEvMax[irof] < imc) {
143-
mcEvMax[irof] = imc;
133+
auto& mc2hit = mc2hitVec[im];
134+
const auto* hv = hitVecPool[im];
135+
for (int ih = (int)hv->size(); ih--;) {
136+
const auto& hit = (*hv)[ih];
137+
uint64_t key = (uint64_t(hit.GetTrackID()) << 32) + hit.GetDetectorID();
138+
mc2hit[key].push_back(ih);
144139
}
145140
}
146141
}
147142

148-
LOGP(info, "Building min and max MC events used by each ROF");
149-
auto pattIt = patternsPtr->cbegin();
150143
for (unsigned int irof = 0; irof < nROFRec; irof++) {
151144
const auto& rofRec = rofRecVec[irof];
152-
// >> read and map MC events contributing to this ROF
153-
for (int im = mcEvMin[irof]; im <= mcEvMax[irof]; im++) {
154-
if (hitVecPool[im] == nullptr) {
155-
hitTree->SetBranchAddress("IT3Hit", &hitVecPool[im]);
156-
hitTree->GetEntry(im);
157-
auto& mc2hit = mc2hitVec[im];
158-
const auto* hitArray = hitVecPool[im];
159-
for (int ih = hitArray->size(); ih--;) {
160-
const auto& hit = (*hitArray)[ih];
161-
uint64_t key = (uint64_t(hit.GetTrackID()) << 32) + hit.GetDetectorID();
162-
mc2hit.emplace(key, ih);
163-
}
164-
}
165-
}
166-
167-
LOGP(debug, "Caching MC events contributing to this ROF {}", irof);
168145
for (int icl = 0; icl < rofRec.getNEntries(); icl++) {
169146
int clEntry = rofRec.getFirstEntry() + icl; // entry of icl-th cluster of this ROF in the vector of clusters
170147
const auto& cluster = (*clusArr)[clEntry];
@@ -175,7 +152,8 @@ void CheckClustersITS3(const std::string& clusfile = "o2clus_its.root",
175152
o2::math_utils::Point3D<float> locC;
176153
auto chipID = cluster.getSensorID();
177154
auto isIB = o2::its3::constants::detID::isDetITS3(chipID);
178-
auto layer = o2::its3::constants::detID::getDetID2Layer(chipID);
155+
(isIB) ? ++cClustersIB : ++cClustersOB;
156+
auto layer = gman->getLayer(chipID);
179157
auto clusterSize{-1};
180158
if (pattID == o2::itsmft::CompCluster::InvalidPatternID || dict.isGroup(pattID, isIB)) {
181159
o2::itsmft::ClusterPattern patt(pattIt);
@@ -199,31 +177,53 @@ void CheckClustersITS3(const std::string& clusfile = "o2clus_its.root",
199177
(isIB) ? ++cLabelInvalidIB : ++cLabelInvalidOB;
200178
continue;
201179
}
180+
auto track = reader.getTrack(lab);
181+
if (!track) {
182+
continue;
183+
}
184+
bool isPrimary = track->isPrimary();
202185

203186
// get MC info
204-
int trID = lab.getTrackID();
187+
const int trID = lab.getTrackID();
188+
const int evID = lab.getEventID();
189+
if (evID < 0 || evID >= (int)mc2hitVec.size()) {
190+
continue;
191+
}
205192
const auto& mc2hit = mc2hitVec[lab.getEventID()];
206-
const auto* hitArray = hitVecPool[lab.getEventID()];
207193
uint64_t key = (uint64_t(trID) << 32) + chipID;
208194
auto hitEntry = mc2hit.find(key);
209195
if (hitEntry == mc2hit.end()) {
210196
LOG(debug) << "Failed to find MC hit entry for Tr" << trID << " chipID" << chipID;
211197
(isIB) ? ++cNoMCIB : ++cNoMCOB;
212198
continue;
213199
}
214-
const auto& hit = (*hitArray)[hitEntry->second];
200+
201+
const o2::itsmft::Hit* hit = nullptr;
202+
int nCand{0};
203+
for (const auto ih : hitEntry->second) {
204+
const auto& candHit = (*hitVecPool[evID])[ih];
205+
const int lay = gman->getLayer(candHit.GetDetectorID());
206+
if (layer == lay) {
207+
hit = &candHit;
208+
++nCand;
209+
}
210+
}
211+
if (hit == nullptr || nCand > 1) {
212+
continue;
213+
}
214+
215215
//
216216
float dx = 0, dz = 0;
217217
int ievH = lab.getEventID();
218218
o2::math_utils::Point3D<float> locH, locHsta;
219219

220-
auto gloH = hit.GetPos();
221-
const auto& gloHsta = hit.GetPosStart();
222-
gloH.SetXYZ(0.5f * (gloH.X() + gloHsta.X()), 0.5f * (gloH.Y() + gloHsta.Y()), 0.5f * (gloH.Z() + gloHsta.Z()));
220+
auto gloH = hit->GetPos();
221+
const auto& gloHsta = hit->GetPosStart();
222+
gloH.SetXYZ((gloH.X() + gloHsta.X()) * 0.5f, 0.5f * (gloH.Y() + gloHsta.Y()), 0.5f * (gloH.Z() + gloHsta.Z()));
223223

224224
// mean local position of the hit
225-
locH = gman->getMatrixL2G(chipID) ^ (hit.GetPos()); // inverse conversion from global to local
226-
locHsta = gman->getMatrixL2G(chipID) ^ (hit.GetPosStart());
225+
locH = gman->getMatrixL2G(chipID) ^ (hit->GetPos()); // inverse conversion from global to local
226+
locHsta = gman->getMatrixL2G(chipID) ^ (hit->GetPosStart());
227227

228228
float x0, y0, z0, dltx, dlty, dltz, r;
229229
if (isIB) {
@@ -250,83 +250,29 @@ void CheckClustersITS3(const std::string& clusfile = "o2clus_its.root",
250250
dltz = locH.Z() - z0;
251251
r = (0.5 * (Segmentation::SensorLayerThickness - Segmentation::SensorLayerThicknessEff) - y0) / dlty;
252252
}
253-
locH.SetXYZ(x0 + r * dltx, y0 + r * dlty, z0 + r * dltz);
253+
locH.SetXYZ(x0 + (r * dltx), y0 + (r * dlty), z0 + (r * dltz));
254254

255255
float theta = std::acos(gloC.Z() / gloC.Rho());
256256
float eta = -std::log(std::tan(theta / 2));
257257

258-
std::array<float, 27> data = {(float)lab.getEventID(), (float)trID,
258+
std::array<float, 28> data = {(float)lab.getEventID(), (float)trID,
259259
locH.X(), locH.Z(),
260260
gloH.X(), gloH.Z(),
261261
dltx / dlty, dltz / dlty,
262262
gloC.X(), gloC.Y(), gloC.Z(),
263263
locC.X(), locC.Y(), locC.Z(),
264264
locC.X() - locH.X(), locC.Y() - locH.Y(), locC.Z() - locH.Z(),
265265
errX, errZ, (float)pattID,
266-
(float)rofRec.getROFrame(), (float)npix, (float)chipID, eta, (float)cluster.getRow(), (float)cluster.getCol(), (float)layer};
266+
(float)rofRec.getROFrame(), (float)npix, (float)chipID, eta, (float)cluster.getRow(), (float)cluster.getCol(), (float)layer, (float)isPrimary};
267267
nt.Fill(data.data());
268268
}
269269
}
270270

271+
LOGP(info, "IB {} clusters and OB {} clusters", cClustersIB, cClustersOB);
271272
LOGP(info, "IB {} valid PatternIDs and {} ({:.1f}%) invalid ones", cPattValidIB, cPattInvalidIB, ((float)cPattInvalidIB / (float)(cPattInvalidIB + cPattValidIB)) * 100);
272273
LOGP(info, "IB {} invalid Labels and {} with No MC Hit information ", cLabelInvalidIB, cNoMCIB);
273274
LOGP(info, "OB {} valid PatternIDs and {} ({:.1f}%) invalid ones", cPattValidOB, cPattInvalidOB, ((float)cPattInvalidOB / (float)(cPattInvalidOB + cPattValidOB)) * 100);
274275
LOGP(info, "OB {} invalid Labels and {} with No MC Hit information ", cLabelInvalidOB, cNoMCOB);
275-
276-
auto canvCgXCgY = new TCanvas("canvCgXCgY", "", 1600, 1600);
277-
canvCgXCgY->Divide(2, 2);
278-
canvCgXCgY->cd(1);
279-
nt.Draw("cgy:cgx>>h_cgy_vs_cgx_IB(1000, -5, 5, 1000, -5, 5)", "id < 3456", "colz");
280-
canvCgXCgY->cd(2);
281-
nt.Draw("cgy:cgz>>h_cgy_vs_cgz_IB(1000, -15, 15, 1000, -5, 5)", "id < 3456", "colz");
282-
canvCgXCgY->cd(3);
283-
nt.Draw("cgy:cgx>>h_cgy_vs_cgx_OB(1000, -50, 50, 1000, -50, 50)", "id >= 3456", "colz");
284-
canvCgXCgY->cd(4);
285-
nt.Draw("cgy:cgz>>h_cgy_vs_cgz_OB(1000, -100, 100, 1000, -50, 50)", "id >= 3456", "colz");
286-
canvCgXCgY->SaveAs("it3clusters_y_vs_x_vs_z.png");
287-
288-
auto canvdXdZ = new TCanvas("canvdXdZ", "", 1600, 800);
289-
canvdXdZ->Divide(2, 2);
290-
canvdXdZ->cd(1)->SetLogz();
291-
nt.Draw("dx:dz>>h_dx_vs_dz_IB(1000, -0.01, 0.01, 1000, -0.01, 0.01)", "id < 3456", "colz");
292-
canvdXdZ->cd(2)->SetLogz();
293-
nt.Draw("dx:dz>>h_dx_vs_dz_OB(1000, -0.01, 0.01, 1000, -0.01, 0.01)", "id >= 3456", "colz");
294-
canvdXdZ->cd(3)->SetLogz();
295-
nt.Draw("dx:dz>>h_dx_vs_dz_IB_z(1000, -0.01, 0.01, 1000, -0.01, 0.01)", "id < 3456 && abs(cgz) < 2", "colz");
296-
canvdXdZ->cd(4)->SetLogz();
297-
nt.Draw("dx:dz>>h_dx_vs_dz_OB_z(1000, -0.01, 0.01, 1000, -0.01, 0.01)", "id >= 3456 && abs(cgz) < 2", "colz");
298-
canvdXdZ->SaveAs("it3clusters_dx_vs_dz.png");
299-
300-
auto canvCHXZ = new TCanvas("canvCHXZ", "", 1600, 1600);
301-
canvCHXZ->Divide(2, 2);
302-
canvCHXZ->cd(1);
303-
nt.Draw("(cgx-hgx)*10000:eta>>h_chx_IB(101,-1.4,1.4,101,-50,50)", "id<3456", "prof");
304-
canvCHXZ->cd(2);
305-
nt.Draw("(cgx-hgx)*10000:eta>>h_chx_OB(101,-1.4,1.4,101,-50,50)", "id>=3456", "prof");
306-
canvCHXZ->cd(3);
307-
nt.Draw("(cgz-hgz)*10000:eta>>h_chz_IB(101,-1.4,1.4,101,-50,50)", "id<3456", "prof");
308-
canvCHXZ->cd(4);
309-
nt.Draw("(cgz-hgz)*10000:eta>>h_chz_OB(101,-1.4,1.4,101,-50,50)", "id>=3456", "prof");
310-
canvCgXCgY->SaveAs("it3clusters_xz_eta.png");
311-
312-
auto c1 = new TCanvas("p1", "pullX");
313-
c1->cd();
314-
c1->SetLogy();
315-
nt.Draw("dx/ex", "abs(dx/ex)<10&&patid<10");
316-
auto c2 = new TCanvas("p2", "pullZ");
317-
c2->cd();
318-
c2->SetLogy();
319-
nt.Draw("dz/ez", "abs(dz/ez)<10&&patid<10");
320-
321-
auto d1 = new TCanvas("d1", "deltaX");
322-
d1->cd();
323-
d1->SetLogy();
324-
nt.Draw("dx", "abs(dx)<5");
325-
auto d2 = new TCanvas("d2", "deltaZ");
326-
d2->cd();
327-
d2->SetLogy();
328-
nt.Draw("dz", "abs(dz)<5");
329-
330276
fout.cd();
331277
nt.Write();
332278
}

0 commit comments

Comments
 (0)