Skip to content

Commit 27410de

Browse files
fprinofprino
authored
Simple macro for ITS tracks vs BC in ROF
1 parent c12208e commit 27410de

1 file changed

Lines changed: 387 additions & 0 deletions

File tree

‎Detectors/ITSMFT/ITS/macros/test/CheckITSTracksVsROF.C‎

Lines changed: 387 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,387 @@
1+
/// \file CheckITSTracksVsROF.C
2+
/// \brief Simple macro to check ITS tracks in AO2D files as a function of bunch crossing within a ROF
3+
4+
#if !defined(__CLING__) || defined(__ROOTCLING__)
5+
#include <TSystem.h>
6+
#include <TFile.h>
7+
#include <TTree.h>
8+
#include <TKey.h>
9+
#include <TH1F.h>
10+
#include <TH2F.h>
11+
#include <TProfile.h>
12+
#include <TCanvas.h>
13+
#include <CCDB/BasicCCDBManager.h>
14+
#include <CommonConstants/LHCConstants.h>
15+
#include "DataFormatsParameters/GRPLHCIFData.h"
16+
#include "DataFormatsFIT/Triggers.h"
17+
#include <DataFormatsFT0/Digit.h>
18+
#include <Framework/DataTypes.h>
19+
#include <DataFormatsParameters/AggregatedRunInfo.h>
20+
#endif
21+
22+
int32_t nBCsPerOrbit = o2::constants::lhc::LHCMaxBunches;
23+
int32_t nBCsPerITSROF = 198;
24+
int32_t offsetITSROF = 64; // from current analysis of delay of ITS ROF with respect to other detectors
25+
constexpr float minCollTimeFT0A = -5.f;
26+
constexpr float maxCollTimeFT0A = 5.f;
27+
constexpr float minCollTimeFT0C = -5.f;
28+
constexpr float maxCollTimeFT0C = 5.f;
29+
constexpr float minCollTimeFV0A = -5.f;
30+
constexpr float maxCollTimeFV0A = 5.f;
31+
constexpr float maxZv = 10.f;
32+
constexpr float maxTrackEta = 0.8f;
33+
34+
void CheckITSTracksVsROF(int maxFiles = 6, int maxDF = 99999)
35+
{
36+
37+
printf("nBCsPerOrbit = %d\n", nBCsPerOrbit);
38+
// define histograms
39+
TH1F* hTrNoColl = new TH1F("hTrNoColl", "", 2, -0.5, 1.5);
40+
hTrNoColl->GetXaxis()->SetBinLabel(1, "No collision assoc.");
41+
hTrNoColl->GetXaxis()->SetBinLabel(2, "With collision assoc.");
42+
TH1F* hTrTime = new TH1F("hTrTime", ";time;entries", 100, -10000., 10000.);
43+
TH1F* hCollTime = new TH1F("hCollTime", ";time;entries", 100, -20., 20.);
44+
TH1F* hVz = new TH1F("hVz", ";z_{v} (cm);entries", 100, -20., 20.);
45+
TH2F* hTrVsCollTime = new TH2F("hTrVsCollTime", ";coll time;track time", 100, -20., 20., 100, -10000., 10000.);
46+
int nBinsROF = 395;
47+
float minROF = -197.5;
48+
float maxROF = 197.5;
49+
TH2F* hCluTrackVsBCinROF = new TH2F("hCluTrackVsBC", ";bcInROF;N_{ITSclusters}", nBinsROF, minROF, maxROF, 8, -0.5, 7.5);
50+
TH1F* hNEventsVsBC = new TH1F("hNEventsVsBC", ";bc;N_{Collisions}", 3601, -0.5, 3600.5);
51+
TH1F* hNTracksVsBC = new TH1F("hNTracksVsBC", ";bc;N_{Tracks}", 3601, -0.5, 3600.5);
52+
TH1F* hNTracksPerEventVsBC = new TH1F("hNTracksPerEventVsBC", ";bc;N_{Tracks} per event", 3601, -0.5, 3600.5);
53+
TH1F* hNEventsVsBCinROF = new TH1F("hNEventsVsBCinROF", ";bcInROF;N_{Collisions}", nBinsROF, minROF, maxROF);
54+
TH1F* hNTracksVsBCinROF = new TH1F("hNTracksVsBCinROF", ";bcInROF;N_{Tracks}", nBinsROF, minROF, maxROF);
55+
TH1F* hNTracksPerEventVsBCinROF = new TH1F("hNTracksPerEventVsBCinROF", ";bcInROF;N_{Tracks} per event", nBinsROF, minROF, maxROF);
56+
TH2F* hNPVContribVsBCinROF = new TH2F("hNPVContribVsBC", ";bcInROF;N_{PVcontrib}", nBinsROF, minROF, maxROF, 151, -0.5, 150.5);
57+
58+
// process all AO2D files in subdirectories of the current directory
59+
gSystem->Exec("find ./ -name AO2D.root | sort > filelist.txt");
60+
int nFiles = 0;
61+
std::ifstream aodList("filelist.txt");
62+
std::string line;
63+
std::vector<std::string> filNames;
64+
while (std::getline(aodList, line)) {
65+
if (!line.empty()) {
66+
filNames.push_back(line);
67+
++nFiles;
68+
}
69+
}
70+
if (nFiles > maxFiles)
71+
nFiles = maxFiles;
72+
73+
for (int jf = 0; jf < nFiles; ++jf) {
74+
TFile* fil = TFile::Open(filNames[jf].c_str());
75+
if (!fil || fil->IsZombie()) {
76+
printf("Cannot open %s\n", filNames[jf].c_str());
77+
continue;
78+
}
79+
TList* l = (TList*)fil->GetListOfKeys();
80+
int nKeys = l->GetEntries();
81+
for (int j = 0; j < nKeys; ++j) {
82+
if (j >= maxDF)
83+
break;
84+
TKey* k = (TKey*)l->At(j);
85+
TString cname = k->GetClassName();
86+
if (cname == "TDirectoryFile") {
87+
TDirectoryFile* d = (TDirectoryFile*)fil->Get(k->GetName());
88+
printf("%s %s\n", fil->GetName(), d->GetName());
89+
90+
// Access TTree with BC info
91+
TTree* tb = (TTree*)d->Get("O2bc_001");
92+
if (!tb) {
93+
printf("TTree O2bc_001 missing\n");
94+
continue;
95+
}
96+
printf(" BC Tree entries = %lld\n", tb->GetEntries());
97+
ULong64_t gloBC;
98+
int nRun;
99+
tb->SetBranchAddress("fRunNumber", &nRun);
100+
tb->SetBranchAddress("fGlobalBC", &gloBC);
101+
102+
// Access TTree with collision info
103+
TTree* tc = (TTree*)d->Get("O2collision_001");
104+
if (!tc) {
105+
printf("TTree O2collision_001 missing\n");
106+
continue;
107+
}
108+
printf(" Collision Tree entries = %lld\n", tc->GetEntries());
109+
float xv, yv, zv, ctime;
110+
int iBC;
111+
ushort nPVcontrib;
112+
tc->SetBranchAddress("fPosX", &xv);
113+
tc->SetBranchAddress("fPosY", &yv);
114+
tc->SetBranchAddress("fPosZ", &zv);
115+
tc->SetBranchAddress("fCollisionTime", &ctime);
116+
tc->SetBranchAddress("fNumContrib", &nPVcontrib);
117+
tc->SetBranchAddress("fIndexBCs", &iBC);
118+
119+
// Access TTree with FT0 info
120+
TTree* tft = (TTree*)d->Get("O2ft0");
121+
if (!tft) {
122+
printf("TTree O2ft0 missing\n");
123+
continue;
124+
}
125+
printf(" FT0 Tree entries = %lld\n", tft->GetEntries());
126+
float timeft0a, timeft0c;
127+
int iBCft0;
128+
UChar_t trigMask;
129+
tft->SetBranchAddress("fIndexBCs", &iBCft0);
130+
tft->SetBranchAddress("fTimeA", &timeft0a);
131+
tft->SetBranchAddress("fTimeC", &timeft0c);
132+
tft->SetBranchAddress("fTriggerMask", &trigMask);
133+
134+
// Access TTree with FV0 info
135+
TTree* tfv = (TTree*)d->Get("O2fv0a");
136+
if (!tfv) {
137+
printf("TTree O2fv0a missing\n");
138+
continue;
139+
}
140+
printf(" FV0A Tree entries = %lld\n", tfv->GetEntries());
141+
float timefv0a;
142+
int iBCfv0;
143+
tfv->SetBranchAddress("fIndexBCs", &iBCfv0);
144+
tfv->SetBranchAddress("fTime", &timefv0a);
145+
146+
// Access TTree with track parameters at their innermost update
147+
TTree* tt = (TTree*)d->Get("O2track_iu");
148+
if (!tt) {
149+
printf("TTree O2track_iu missing\n");
150+
continue;
151+
}
152+
printf(" Track Tree entries = %lld\n", tt->GetEntries());
153+
// 5 parameters defining the track helix: y, z, sin(phi), tan(lambda), q/pt
154+
// x is a coordinate along the track
155+
// alpha: angle between track reference system and global reference system
156+
float x, alp, y, z, snp, tgl, qpt;
157+
int iColl;
158+
tt->SetBranchAddress("fIndexCollisions", &iColl);
159+
tt->SetBranchAddress("fX", &x);
160+
tt->SetBranchAddress("fAlpha", &alp);
161+
tt->SetBranchAddress("fY", &y);
162+
tt->SetBranchAddress("fZ", &z);
163+
tt->SetBranchAddress("fSnp", &snp);
164+
tt->SetBranchAddress("fTgl", &tgl);
165+
tt->SetBranchAddress("fSigned1Pt", &qpt);
166+
167+
// Access TTree with track extra information
168+
TTree* te = (TTree*)d->Get("O2trackextra_002");
169+
if (!te)
170+
te = (TTree*)d->Get("O2trackextra_001");
171+
if (!te) {
172+
printf("TTree O2trackextra_002 and O2trackextra_001 both missing\n");
173+
continue;
174+
}
175+
uint itsclusiz;
176+
uint8_t nTPCclu;
177+
uint trflag;
178+
float itschi2, tpcchi2, trtime;
179+
te->SetBranchAddress("fITSClusterSizes", &itsclusiz);
180+
te->SetBranchAddress("fITSChi2NCl", &itschi2);
181+
te->SetBranchAddress("fTPCNClsFindable", &nTPCclu);
182+
te->SetBranchAddress("fTPCChi2NCl", &tpcchi2);
183+
te->SetBranchAddress("fTrackTime", &trtime);
184+
te->SetBranchAddress("fFlags", &trflag);
185+
186+
// get run information
187+
if (tb->GetEntries() == 0) {
188+
printf("Empty BC tree\n");
189+
continue;
190+
}
191+
tb->GetEntry(0);
192+
auto runInfo = o2::parameters::AggregatedRunInfo::buildAggregatedRunInfo(o2::ccdb::BasicCCDBManager::instance(), nRun);
193+
// uint64_t sorTimestamp = runInfo.sor;
194+
// uint64_t eorTimestamp = runInfo.eor;
195+
int64_t bcSOR = runInfo.orbitSOR * nBCsPerOrbit;
196+
int64_t nBCsPerTF = runInfo.orbitsPerTF * nBCsPerOrbit;
197+
auto calcBcInROF = [&](ULong64_t bcGlobal) {
198+
int bc = (bcGlobal + nBCsPerOrbit) % nBCsPerITSROF;
199+
bc -= offsetITSROF;
200+
return bc;
201+
};
202+
203+
// fill maps for BC and FT0
204+
int nFT0 = tft->GetEntries();
205+
std::unordered_map<int, int> bcToFT0;
206+
bcToFT0.reserve(nFT0);
207+
for (int jft = 0; jft < nFT0; ++jft) {
208+
tft->GetEntry(jft);
209+
bcToFT0[iBCft0] = jft;
210+
}
211+
// fill maps for BC and FV0
212+
int nFV0 = tfv->GetEntries();
213+
std::unordered_map<int, int> bcToFV0;
214+
bcToFV0.reserve(nFV0);
215+
for (int jfv = 0; jfv < nFV0; ++jfv) {
216+
tfv->GetEntry(jfv);
217+
bcToFV0[iBCfv0] = jfv;
218+
}
219+
// tag good collisions
220+
int nColl = tc->GetEntries();
221+
std::vector<int> ft0indices(nColl, -1);
222+
std::vector<int> fv0indices(nColl, -1);
223+
std::vector<bool> isGoodColl(nColl, false);
224+
for (int i = 0; i < nColl; ++i) {
225+
ft0indices[i] = -1;
226+
tc->GetEntry(i);
227+
if (iBC < 0 || iBC >= tb->GetEntries()) {
228+
printf("ERROR: iBC out of range\n");
229+
continue;
230+
}
231+
tb->GetEntry(iBC);
232+
auto it = bcToFT0.find(iBC);
233+
if (it != bcToFT0.end())
234+
ft0indices[i] = it->second;
235+
auto iv = bcToFV0.find(iBC);
236+
if (iv != bcToFV0.end())
237+
fv0indices[i] = iv->second;
238+
if (std::abs(zv) > maxZv)
239+
continue;
240+
int jft0 = ft0indices[i];
241+
int jfv0 = fv0indices[i];
242+
if (jft0 >= 0 && jfv0 >= 0) {
243+
tft->GetEntry(jft0);
244+
tfv->GetEntry(jfv0);
245+
bool isTVX = TESTBIT(trigMask, o2::ft0::Triggers::bitVertex);
246+
int64_t bcInTF = (static_cast<int64_t>(gloBC) - bcSOR) % nBCsPerTF;
247+
bool isTFBorderOK = (bcInTF > 300 && bcInTF < nBCsPerTF - 4000);
248+
if (isTFBorderOK && isTVX &&
249+
timeft0a > minCollTimeFT0A && timeft0a < maxCollTimeFT0A &&
250+
timeft0c > minCollTimeFT0C && timeft0c < maxCollTimeFT0C &&
251+
timefv0a > minCollTimeFV0A && timefv0a < maxCollTimeFV0A) {
252+
isGoodColl[i] = true;
253+
int bcInITSROF = calcBcInROF(gloBC);
254+
hNEventsVsBC->Fill(gloBC % nBCsPerOrbit);
255+
hNEventsVsBCinROF->Fill(bcInITSROF);
256+
hCollTime->Fill(ctime);
257+
hVz->Fill(zv);
258+
}
259+
}
260+
}
261+
262+
// loop over tracks
263+
int nTracks = tt->GetEntries();
264+
for (int i = 0; i < nTracks; ++i) {
265+
// O2track_iu and O2trackextra_xxx entries are aligned 1-to-1
266+
tt->GetEntry(i);
267+
te->GetEntry(i);
268+
// compute the ITS cluster map, which can be used in track selecion
269+
int nITSclu = 0;
270+
uint itsCluMap = 0;
271+
for (int jLay = 0; jLay < 7; ++jLay) {
272+
if ((itsclusiz >> (jLay * 4)) & 0xf) {
273+
nITSclu++;
274+
itsCluMap |= (1 << jLay);
275+
}
276+
}
277+
// track selections
278+
if (!(trflag & o2::aod::track::PVContributor))
279+
continue;
280+
if (nITSclu == 0 || itschi2 < 0. || tpcchi2 < 0.)
281+
continue;
282+
// compute eta from track dip angle lambda via tgl = tan(lambda)
283+
float cl = 1. / std::sqrt(1. + tgl * tgl);
284+
float sl = cl * tgl;
285+
float eta = 0.5 * std::log((1 + sl) / (1 - sl));
286+
if (std::abs(eta) > maxTrackEta)
287+
continue;
288+
// collision selections
289+
if (iColl >= 0 && iColl < nColl && isGoodColl[iColl]) {
290+
tc->GetEntry(iColl);
291+
hTrTime->Fill(trtime);
292+
hTrVsCollTime->Fill(ctime, trtime);
293+
tb->GetEntry(iBC);
294+
int bcInITSROF = calcBcInROF(gloBC);
295+
hNTracksVsBC->Fill(gloBC % nBCsPerOrbit);
296+
hCluTrackVsBCinROF->Fill(bcInITSROF, nITSclu);
297+
hNTracksVsBCinROF->Fill(bcInITSROF);
298+
hNPVContribVsBCinROF->Fill(bcInITSROF, nPVcontrib);
299+
hTrNoColl->Fill(1);
300+
} else {
301+
hTrNoColl->Fill(0);
302+
}
303+
}
304+
}
305+
}
306+
}
307+
308+
// plotting
309+
TProfile* pctr = hNPVContribVsBCinROF->ProfileX();
310+
TProfile* pclu = hCluTrackVsBCinROF->ProfileX();
311+
TH1F* hFrac7 = new TH1F("hFrac7", ";bcInROF;Fraction of tracks with 7 clusters", hCluTrackVsBCinROF->GetXaxis()->GetNbins(), hCluTrackVsBCinROF->GetXaxis()->GetXmin(), hCluTrackVsBCinROF->GetXaxis()->GetXmax());
312+
for (int jx = 1; jx <= hCluTrackVsBCinROF->GetXaxis()->GetNbins(); jx++) {
313+
double tot = 0.;
314+
for (int jy = 1; jy <= hCluTrackVsBCinROF->GetYaxis()->GetNbins(); jy++) {
315+
double c = hCluTrackVsBCinROF->GetBinContent(jx, jy);
316+
tot += c;
317+
}
318+
double seven = hCluTrackVsBCinROF->GetBinContent(jx, 8);
319+
double f = 0.;
320+
double ef = 0.;
321+
if (tot > 0) {
322+
f = seven / tot;
323+
ef = std::sqrt(f * (1 - f) / tot);
324+
}
325+
hFrac7->SetBinContent(jx, f);
326+
hFrac7->SetBinError(jx, ef);
327+
}
328+
for (int jx = 1; jx <= hNTracksVsBC->GetXaxis()->GetNbins(); jx++) {
329+
double nt = hNTracksVsBC->GetBinContent(jx);
330+
double ne = hNEventsVsBC->GetBinContent(jx);
331+
if (ne > 0)
332+
hNTracksPerEventVsBC->SetBinContent(jx, nt / ne);
333+
}
334+
for (int jx = 1; jx <= hNTracksVsBCinROF->GetXaxis()->GetNbins(); jx++) {
335+
double nt = hNTracksVsBCinROF->GetBinContent(jx);
336+
double ne = hNEventsVsBCinROF->GetBinContent(jx);
337+
if (ne > 0)
338+
hNTracksPerEventVsBCinROF->SetBinContent(jx, nt / ne);
339+
}
340+
341+
hCluTrackVsBCinROF->SetStats(0);
342+
hFrac7->SetStats(0);
343+
hNPVContribVsBCinROF->SetStats(0);
344+
345+
TCanvas* ctr = new TCanvas("ctr", "", 1400, 500);
346+
ctr->Divide(3, 1);
347+
ctr->cd(1);
348+
hTrNoColl->Draw();
349+
ctr->cd(2);
350+
hCollTime->Draw();
351+
ctr->cd(3);
352+
hVz->Draw();
353+
354+
TCanvas* cbc = new TCanvas("cbc", "", 1400, 800);
355+
cbc->Divide(1, 3);
356+
cbc->cd(1);
357+
hNEventsVsBC->Draw();
358+
cbc->cd(2);
359+
hNTracksVsBC->Draw();
360+
cbc->cd(3);
361+
hNTracksPerEventVsBC->Draw();
362+
363+
TCanvas* crof = new TCanvas("crof", "", 1400, 800);
364+
crof->Divide(3, 2);
365+
crof->cd(1);
366+
hNEventsVsBCinROF->Draw();
367+
crof->cd(2);
368+
hNTracksVsBCinROF->Draw();
369+
crof->cd(3);
370+
hNTracksPerEventVsBCinROF->Draw();
371+
crof->cd(4);
372+
gPad->SetLogz();
373+
hCluTrackVsBCinROF->Draw("colz");
374+
pclu->SetLineWidth(2);
375+
pclu->Draw("same");
376+
crof->cd(5);
377+
hFrac7->SetMinimum(0.);
378+
hFrac7->SetMaximum(1.);
379+
hFrac7->SetLineWidth(2);
380+
hFrac7->Draw();
381+
crof->cd(6);
382+
gPad->SetLogz();
383+
hNPVContribVsBCinROF->Draw("colz");
384+
pctr->SetLineWidth(2);
385+
pctr->Draw("same");
386+
crof->SaveAs("ITSTracksVsBCinROF.png");
387+
}

0 commit comments

Comments
 (0)