Skip to content

Commit 0302436

Browse files
maciaccomaciacco
committed
add first version of QA macro to check digitization
1 parent 25f4dfe commit 0302436

3 files changed

Lines changed: 308 additions & 0 deletions

File tree

‎Detectors/Upgrades/ALICE3/IOTOF/macros/CMakeLists.txt‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -14,3 +14,14 @@ o2_add_test_root_macro(defineIOTOFGeo.C
1414

1515
o2_add_test_root_macro(drawTOFGeometry.C
1616
LABELS alice3)
17+
18+
o2_add_test_root_macro(CheckDigitsIOTOF.C
19+
PUBLIC_LINK_LIBRARIES O2::ITSMFTBase
20+
O2::ITSMFTSimulation
21+
O2::IOTOFBase
22+
O2::IOTOFSimulation
23+
O2::MathUtils
24+
O2::SimulationDataFormat
25+
O2::DetectorsBase
26+
O2::Steer
27+
LABELS iotof COMPILE_ONLY)

‎Detectors/Upgrades/ALICE3/IOTOF/macros/CheckDigitsIOTOF.C‎

Lines changed: 296 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,296 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
12+
/// \file CheckDigitsIOTOF.C
13+
/// \brief Simple macro to check TF3 digits
14+
15+
#if !defined(__CLING__) || defined(__ROOTCLING__)
16+
#include <TCanvas.h>
17+
#include <TFile.h>
18+
#include <TH2F.h>
19+
#include <TNtuple.h>
20+
#include <TString.h>
21+
#include <TTree.h>
22+
#include <TLine.h>
23+
#include <TStyle.h>
24+
25+
#include "IOTOFSimulation/Segmentation.h"
26+
#include "IOTOFBase/IOTOFBaseParameters.h"
27+
#include "IOTOFBase/GeometryTGeo.h"
28+
#include "DataFormatsIOTOF/Digit.h"
29+
#include "ITSMFTSimulation/Hit.h"
30+
#include "MathUtils/Utils.h"
31+
#include "SimulationDataFormat/ConstMCTruthContainer.h"
32+
#include "SimulationDataFormat/IOMCTruthContainerView.h"
33+
#include "SimulationDataFormat/MCCompLabel.h"
34+
#include "DetectorsBase/GeometryManager.h"
35+
#include "CCDB/BasicCCDBManager.h"
36+
37+
#include "DataFormatsITSMFT/ROFRecord.h"
38+
39+
#endif
40+
41+
#define ENABLE_UPGRADES
42+
43+
void addTLines(float pitch)
44+
{
45+
// Add grid lines at multiples of pitch on the current pad
46+
if (!gPad)
47+
return;
48+
49+
gPad->Update();
50+
51+
Double_t xmin = gPad->GetUxmin();
52+
Double_t xmax = gPad->GetUxmax();
53+
Double_t ymin = gPad->GetUymin();
54+
Double_t ymax = gPad->GetUymax();
55+
56+
// Calculate the first vertical line position (multiple of pitch)
57+
int nLinesX = 0;
58+
for (float x = xmin; x <= xmax && nLinesX < 1000; x += pitch, nLinesX++) {
59+
TLine* line = new TLine(x, ymin, x, ymax);
60+
line->SetLineStyle(2);
61+
line->SetLineColor(kGray);
62+
line->Draw("same");
63+
}
64+
65+
// Calculate the first horizontal line position (multiple of pitch)
66+
int nLinesY = 0;
67+
for (float y = ymin; y <= ymax && nLinesY < 1000; y += pitch, nLinesY++) {
68+
TLine* line = new TLine(xmin, y, xmax, y);
69+
line->SetLineStyle(2);
70+
line->SetLineColor(kGray);
71+
line->Draw("same");
72+
}
73+
74+
gPad->Modified();
75+
gPad->Update();
76+
}
77+
78+
void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", std::string hitfile = "o2sim_HitsTF3.root", std::string inputGeom = "o2sim_geometry.root")
79+
{
80+
gStyle->SetPalette(55);
81+
82+
using namespace o2::base;
83+
using namespace o2::iotof;
84+
85+
using o2::iotof::Digit;
86+
using o2::itsmft::Hit;
87+
88+
o2::conf::ConfigurableParam::updateFromString("IOTOFBase.segmentedInnerTOF=true;IOTOFBase.segmentedOuterTOF=true;IOTOFBase.enableForwardTOF=false;IOTOFBase.enableBackwardTOF=false");
89+
90+
auto seg = o2::iotof::Segmentation::Instance();
91+
92+
TFile* f = TFile::Open("CheckDigits.root", "recreate");
93+
94+
TNtuple* nt = new TNtuple("ntd", "digit ntuple", "id:x:y:z:rowD:colD:rowH:colH:xlH:zlH:xlcH:zlcH:dx:dz");
95+
TNtuple* nt2 = new TNtuple("ntd2", "digit ntuple", "id:z:dxH:dzH"); /// maximum number of elements in a tuple = 15: doing a new tuple to store more variables
96+
97+
auto& iotofPars = IOTOFBaseParam::Instance();
98+
99+
// Geometry
100+
o2::base::GeometryManager::loadGeometry(inputGeom);
101+
auto* gman = o2::iotof::GeometryTGeo::Instance();
102+
gman->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G));
103+
104+
// IOTOF response plane y = half sensor thickness
105+
float depthMax = iotofPars.sensorThickness; // fallback (no CCDB)
106+
const float yPlaneIOTOF = depthMax / 2.;
107+
108+
// Hits
109+
TFile* hitFile = TFile::Open(hitfile.data());
110+
TTree* hitTree = (TTree*)hitFile->Get("o2sim");
111+
int nevH = hitTree->GetEntries(); // hits are stored as one event per entry
112+
std::vector<std::vector<o2::itsmft::Hit>*> hitArray(nevH, nullptr);
113+
114+
std::vector<std::unordered_map<uint64_t, int>> mc2hitVec(nevH);
115+
116+
// Digits
117+
TFile* digFile = TFile::Open(digifile.data());
118+
TTree* digTree = (TTree*)digFile->Get("o2sim");
119+
120+
std::vector<o2::iotof::Digit>* digArr{nullptr};
121+
std::vector<o2::itsmft::ROFRecord>* rofRecordsArr{nullptr};
122+
o2::dataformats::IOMCTruthContainerView* plabelsArr{nullptr};
123+
124+
digTree->SetBranchAddress("TF3Digit", &digArr);
125+
digTree->SetBranchAddress("TF3DigitROF", &rofRecordsArr);
126+
digTree->SetBranchAddress("TF3DigitMCTruth", &plabelsArr);
127+
128+
digTree->GetEntry(0);
129+
130+
// Load all MC hit events upfront and build the hit lookup map.
131+
for (int im = 0; im < nevH; ++im) {
132+
hitTree->SetBranchAddress("TF3Hit", &hitArray[im]);
133+
hitTree->GetEntry(im);
134+
auto& mc2hit = mc2hitVec[im];
135+
for (int ih = hitArray[im]->size(); ih--;) {
136+
const auto& hit = (*hitArray[im])[ih];
137+
uint64_t key = (uint64_t(hit.GetTrackID()) << 32) + hit.GetDetectorID();
138+
mc2hit.emplace(key, ih);
139+
}
140+
}
141+
142+
auto& rofArr = *rofRecordsArr;
143+
const int nROFRec = (int)rofArr.size();
144+
145+
o2::dataformats::ConstMCTruthContainer<o2::MCCompLabel> labels;
146+
plabelsArr->copyandflatten(labels);
147+
148+
// LOOP on : ROFRecord array
149+
for (unsigned int iROF = 0; iROF < rofArr.size(); ++iROF) {
150+
151+
const unsigned int rofIndex = rofArr[iROF].getFirstEntry();
152+
const unsigned int rofNEntries = rofArr[iROF].getNEntries();
153+
154+
// LOOP on : digits array
155+
for (unsigned int iDigit = rofIndex; iDigit < rofIndex + rofNEntries; iDigit++) {
156+
if (iDigit % 1000 == 0) {
157+
std::cout << "Reading digit " << iDigit << " / " << digArr->size() << std::endl;
158+
}
159+
160+
Int_t ix = (*digArr)[iDigit].getRow(), iz = (*digArr)[iDigit].getColumn();
161+
Int_t iDetID = (*digArr)[iDigit].getChipIndex();
162+
Int_t subDetID = gman->getIOTOFLayer(iDetID);
163+
164+
Float_t x = 0.f, y = 0.f, z = 0.f;
165+
Float_t x_flat = 0.f, z_flat = 0.f;
166+
167+
// Float_t t = (*digArr)[iDigit].getTime();
168+
169+
if (subDetID >= 0) {
170+
seg->detectorToLocal(ix, iz, x, z, subDetID);
171+
}
172+
173+
o2::math_utils::Point3D<float> locD(x, y, z); // local Digit
174+
175+
Int_t chipID = (*digArr)[iDigit].getChipIndex();
176+
177+
auto lab = (labels.getLabels(iDigit))[0];
178+
179+
if (!lab.isValid()) { // not a noise
180+
continue;
181+
}
182+
183+
int trID = lab.getTrackID();
184+
185+
const auto gloD = gman->getMatrixL2G(chipID)(locD); // convert to global
186+
187+
std::unordered_map<uint64_t, int>* mc2hit = &mc2hitVec[lab.getEventID()];
188+
189+
// get MC info
190+
uint64_t key = (uint64_t(trID) << 32) + chipID;
191+
auto hitEntry = mc2hit->find(key);
192+
193+
if (hitEntry == mc2hit->end()) {
194+
LOG(error) << "Failed to find MC hit entry for Tr" << trID << " chipID" << chipID;
195+
continue;
196+
}
197+
198+
////// HITS
199+
Hit& hit = (*hitArray[lab.getEventID()])[hitEntry->second];
200+
201+
auto xyzLocE = gman->getMatrixL2G(chipID) ^ (hit.GetPos()); // inverse conversion from global to local
202+
auto xyzLocS = gman->getMatrixL2G(chipID) ^ (hit.GetPosStart());
203+
204+
// Hit local reference: use response plane interpolation
205+
o2::math_utils::Vector3D<float> locH; /// Hit reference (at response plane)
206+
o2::math_utils::Vector3D<float> locHS; /// Hit, start pos
207+
locHS.SetCoordinates(xyzLocS.X(), xyzLocS.Y(), xyzLocS.Z());
208+
o2::math_utils::Vector3D<float> locHE; /// Hit, end pos
209+
locHE.SetCoordinates(xyzLocE.X(), xyzLocE.Y(), xyzLocE.Z());
210+
o2::math_utils::Vector3D<float> locHF;
211+
212+
// IOTOF: Interpolate to response plane
213+
float x0 = locHS.X(), y0 = locHS.Y(), z0 = locHS.Z();
214+
float dltx = locHE.X() - x0, dlty = locHE.Y() - y0, dltz = locHE.Z() - z0;
215+
float r = (std::abs(dlty) > 1e-9f) ? (yPlaneIOTOF - y0) / dlty : 0.5f;
216+
locH.SetCoordinates(x0 + r * dltx, yPlaneIOTOF, z0 + r * dltz);
217+
218+
int row = 0, col = 0;
219+
float xlc = 0., zlc = 0.;
220+
221+
seg->localToDetector(locH.X(), locH.Z(), row, col, subDetID);
222+
seg->detectorToLocal(row, col, xlc, zlc, subDetID);
223+
nt->Fill(chipID, /// detector ID
224+
gloD.X(), gloD.Y(), gloD.Z(), /// global position retrieved from the digit: digit (row, col) ->local position -> global potision
225+
ix, iz, /// row and column of the digit
226+
row, col, /// row and col retrieved from the hit: hit global position -> hit local position -> detector position (row, col)
227+
locH.X(), locH.Z(), /// x and z of the hit in the local reference frame: hit global position -> hit local position
228+
xlc, zlc, /// x and z of the hit in the local frame: hit global position -> hit local position -> detector position (row, col) -> local position
229+
locHS.X() - locD.X(), locHS.Z() - locD.Z()); /// difference in x and z between the hit and the digit in the local frame
230+
nt2->Fill(chipID, gloD.Z(), locHS.X() - locHE.X(), locHS.Z() - locHE.Z()); /// differences between local hit start and hit end positions
231+
232+
} // end loop on digits array
233+
234+
} // end loop on ROFRecords
235+
236+
// digit maps in the xy and yz planes
237+
auto canvXY = new TCanvas("canvXY", "", 1600, 800);
238+
canvXY->Divide(2, 1);
239+
canvXY->cd(1);
240+
nt->Draw("y:x>>h_y_vs_x_IOTOF(1000, -100, 100, 1000, -100, 100)", "id >= 0 && id < 14352", "colz");
241+
canvXY->cd(2);
242+
nt->Draw("y:z>>h_y_vs_z_IOTOF(1000, -400, 400, 1000, -100, 100)", "id >= 0 && id < 14352", "colz");
243+
canvXY->SaveAs("tf3digits_y_vs_x_vs_z.pdf");
244+
245+
// z distributions
246+
auto canvZ = new TCanvas("canvZ", "", 800, 800);
247+
canvZ->cd();
248+
nt->Draw("z>>h_z_IOTOF(500, -70, 70)", "id >= 0 && id < 14352 ");
249+
canvZ->SaveAs("tf3digits_z.pdf");
250+
251+
// dz distributions (difference between local position of digits and hits in x and z)
252+
auto canvdZ = new TCanvas("canvdZ", "", 800, 800);
253+
canvdZ->cd();
254+
nt->Draw("dz>>h_dz_ML(500, -0.05, 0.05)", "id >= 0 && id < 14352 ");
255+
canvdZ->SaveAs("tf3digits_dz.pdf");
256+
canvdZ->SaveAs("tf3digits_dz.root");
257+
258+
// distributions of differences between local positions of digits and hits in x and z
259+
auto canvdXdZ = new TCanvas("canvdXdZ", "", 1600, 800);
260+
canvdXdZ->Divide(2, 1);
261+
canvdXdZ->cd(1);
262+
nt->Draw("dx:dz>>h_dx_vs_dz_ITOF(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 0 && id < 960", "colz");
263+
addTLines(0.01);
264+
auto h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_ITOF");
265+
Info("ITOF", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4);
266+
Info("ITOF", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4);
267+
canvdXdZ->cd(2);
268+
nt->Draw("dx:dz>>h_dx_vs_dz_OTOF(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 960 && id < 14352", "colz");
269+
addTLines(0.01);
270+
h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_OTOF");
271+
Info("OTOF", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4);
272+
Info("OTOF", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4);
273+
canvdXdZ->SaveAs("tf3digits_dx_vs_dz.pdf");
274+
canvdXdZ->SaveAs("tf3digits_dx_vs_dz.root");
275+
276+
// distribution of differences between hit start and hit end in local coordinates
277+
auto canvdXdZHit = new TCanvas("canvdXdZHit", "", 1600, 800);
278+
canvdXdZHit->Divide(2, 1);
279+
canvdXdZHit->cd(1);
280+
LOG(info) << "dxH, dzH";
281+
nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_ITOF(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 0 && id < 960", "colz");
282+
addTLines(0.01);
283+
h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_ITOF");
284+
Info("ITOF", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4);
285+
Info("ITOF", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4);
286+
canvdXdZHit->cd(2);
287+
nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_OTOF(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 960 && id < 14352", "colz");
288+
addTLines(0.01);
289+
h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_OTOF");
290+
Info("OTOF", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4);
291+
Info("OTOF", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4);
292+
canvdXdZHit->SaveAs("trkdigits_dxH_vs_dzH.pdf");
293+
294+
f->Write();
295+
f->Close();
296+
}

‎Detectors/Upgrades/ALICE3/IOTOF/simulation/src/IOTOFSimulationLinkDef.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,7 @@
2323
#pragma link C++ class o2::base::DetImpl < o2::iotof::Detector> + ;
2424

2525
#pragma link C++ class o2::iotof::Digitizer + ;
26+
#pragma link C++ class o2::iotof::Segmentation + ;
2627
#pragma link C++ class o2::iotof::DPLDigitizerParam + ;
2728
#pragma link C++ class o2::conf::ConfigurableParamHelper < o2::iotof::DPLDigitizerParam> + ;
2829

0 commit comments

Comments
 (0)