Skip to content

Commit a365645

Browse files
njacazionjacazio
authored
Add writer for the LUT from AODs (#79)
- Using the binning of the input histograms
1 parent 952ab74 commit a365645

File tree

2 files changed

+321
-0
lines changed

2 files changed

+321
-0
lines changed

‎examples/scripts/createO2tables.py‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -64,6 +64,9 @@ def main(configuration_file,
6464
arguments = locals() # List of arguments to put into the log
6565
parser = configparser.RawConfigParser()
6666
parser.read(configuration_file)
67+
if config_entry not in parser.keys():
68+
fatal_msg(f"Did not find configuration entry '{config_entry}' in config file",
69+
configuration_file + "\n\t Available entries:\n\t\t" + "\n\t\t".join(list(parser.keys())))
6770

6871
run_cmd("./clean.sh > /dev/null 2>&1", check_status=False)
6972
# Dictionary of fetched options

‎src/lutWrite_aod.cc‎

Lines changed: 318 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,318 @@
1+
/// @author: Roberto Preghenella
2+
/// @email: [email protected]
3+
4+
/// @author: Nicolo' Jacazio
5+
/// @email: [email protected]
6+
7+
#include "TCanvas.h"
8+
#include "TDatabasePDG.h"
9+
#include "TFile.h"
10+
#include "TMatrixD.h"
11+
#include "TMatrixDSymEigen.h"
12+
#include "TProfile2D.h"
13+
#include "TProfile3D.h"
14+
#include "TVectorD.h"
15+
#include "fwdRes/fwdRes.C"
16+
#include "lutCovm.hh"
17+
#include <Riostream.h>
18+
19+
void diagonalise(lutEntry_t& lutEntry);
20+
21+
bool fwdSolve(float* covm, float pt = 0.1, float eta = 0.0,
22+
float mass = 0.13957000)
23+
{
24+
if (fwdRes(covm, pt, eta, mass) < 0)
25+
return false;
26+
return true;
27+
}
28+
29+
void lutWrite_aod(const char* filename = "/tmp/lutCovm.pi.aod.dat",
30+
int pdg = 211,
31+
float field = 0.2, int layer = 0, int what = 0,
32+
int efftype = 0,
33+
const char* infilename = "/tmp/AnalysisResults_LUT.root",
34+
float minPt = 0.f,
35+
float maxPt = 80.f,
36+
float minEta = -4.f,
37+
float maxEta = 4.f)
38+
{
39+
40+
std::map<int, std::string> partname{ { 11, "electron" }, { 13, "muon" }, { 211, "pion" }, { 321, "kaon" }, { 2212, "proton" } };
41+
const std::string dn = "alice3-lutmaker-" + partname[pdg];
42+
43+
// Get the input from the analysis results
44+
TFile f(infilename);
45+
if (!f.IsOpen()) {
46+
Printf("Did not find %s", infilename);
47+
return;
48+
}
49+
// f.ls();
50+
TDirectory* d = nullptr;
51+
f.GetObject(dn.c_str(), d);
52+
if (!d) {
53+
Printf("Did not find %s", dn.c_str());
54+
f.ls();
55+
return;
56+
}
57+
// d->ls();
58+
std::map<std::string, TH1F*> h{ { "eta", nullptr }, { "pt", nullptr } };
59+
std::map<std::string, TProfile2D*> m{ { "CovMat_cYY", nullptr },
60+
{ "CovMat_cZY", nullptr },
61+
{ "CovMat_cZZ", nullptr },
62+
{ "CovMat_cSnpY", nullptr },
63+
{ "CovMat_cSnpZ", nullptr },
64+
{ "CovMat_cSnpSnp", nullptr },
65+
{ "CovMat_cTglY", nullptr },
66+
{ "CovMat_cTglZ", nullptr },
67+
{ "CovMat_cTglSnp", nullptr },
68+
{ "CovMat_cTglTgl", nullptr },
69+
{ "CovMat_c1PtY", nullptr },
70+
{ "CovMat_c1PtZ", nullptr },
71+
{ "CovMat_c1PtSnp", nullptr },
72+
{ "CovMat_c1PtTgl", nullptr },
73+
{ "CovMat_c1Pt21Pt2", nullptr },
74+
{ "Efficiency", nullptr } };
75+
76+
struct binning {
77+
int n = 0;
78+
float min = 0;
79+
float max = 0;
80+
bool log = false;
81+
};
82+
83+
binning histo_eta_bins;
84+
binning histo_pt_bins;
85+
for (auto const& i : h) {
86+
87+
auto setBinning = [&h, &i](binning& b) {
88+
const auto j = h[i.first];
89+
if (b.n == 0) {
90+
Printf("Setting bin for %s", i.first.c_str());
91+
b.n = j->GetXaxis()->GetNbins();
92+
b.min = j->GetXaxis()->GetBinLowEdge(1);
93+
b.max = j->GetXaxis()->GetBinUpEdge(b.n);
94+
}
95+
if (std::abs(j->GetXaxis()->GetBinWidth(1) - j->GetXaxis()->GetBinWidth(j->GetNbinsX())) > 0.0001f) {
96+
b.log = true;
97+
}
98+
};
99+
100+
h[i.first] = ((TH1F*)d->Get(i.first.c_str()));
101+
h[i.first]->SetDirectory(0);
102+
if (i.first == "eta") {
103+
setBinning(histo_eta_bins);
104+
} else if (i.first == "pt") {
105+
setBinning(histo_pt_bins);
106+
}
107+
}
108+
109+
for (auto const& i : m) {
110+
auto checkBinning = [&m, &i, histo_eta_bins, histo_pt_bins]() {
111+
const auto j = m[i.first];
112+
const char* n = i.first.c_str();
113+
// X
114+
const TAxis* x = j->GetXaxis();
115+
if (histo_pt_bins.n != x->GetNbins()) {
116+
Printf("Different number of bins on X for %s: %i vs %i", n, histo_pt_bins.n, x->GetNbins());
117+
return false;
118+
}
119+
if (std::abs(histo_pt_bins.min - x->GetBinLowEdge(1)) > 0.0001f) {
120+
Printf("Different starting on X for %s: %f vs %f, diff. is %f", n, histo_pt_bins.min, x->GetBinLowEdge(1), histo_pt_bins.min - x->GetBinLowEdge(1));
121+
return false;
122+
}
123+
if (std::abs(histo_pt_bins.max - x->GetBinUpEdge(x->GetNbins())) > 0.0001f) {
124+
Printf("Different ending on X for %s: %f vs %f, diff. is %f", n, histo_pt_bins.max, x->GetBinUpEdge(x->GetNbins()), histo_pt_bins.max - x->GetBinUpEdge(x->GetNbins()));
125+
return false;
126+
}
127+
// Y
128+
const TAxis* y = j->GetYaxis();
129+
if (histo_eta_bins.n != y->GetNbins()) {
130+
Printf("Different number of bins on Y for %s: %i vs %i", n, histo_eta_bins.n, y->GetNbins());
131+
return false;
132+
}
133+
if (std::abs(histo_eta_bins.min - y->GetBinLowEdge(1)) > 0.0001f) {
134+
Printf("Different starting on Y for %s: %f vs %f, diff. is %f", n, histo_eta_bins.min, y->GetBinLowEdge(1), histo_eta_bins.min - y->GetBinLowEdge(1));
135+
return false;
136+
}
137+
if (std::abs(histo_eta_bins.max - y->GetBinUpEdge(y->GetNbins())) > 0.0001f) {
138+
Printf("Different ending on Y for %s: %f vs %f, diff. is %f", n, histo_eta_bins.max, y->GetBinUpEdge(y->GetNbins()), histo_eta_bins.max - y->GetBinUpEdge(y->GetNbins()));
139+
return false;
140+
}
141+
return true;
142+
};
143+
// m[i.first] = (TProfile3D*)d->Get(i.first.c_str());
144+
// m[i.first] = ((TProfile3D*)d->Get(i.first.c_str()))->Project3DProfile("yx");
145+
m[i.first] = ((TProfile2D*)d->Get(i.first.c_str()));
146+
m[i.first]->SetDirectory(0);
147+
if (!checkBinning()) {
148+
Printf("Something went wrong, stopping");
149+
return;
150+
}
151+
}
152+
153+
f.Close();
154+
155+
// output file
156+
ofstream lutFile(filename, std::ofstream::binary);
157+
158+
// write header
159+
lutHeader_t lutHeader;
160+
// pid
161+
lutHeader.pdg = pdg;
162+
lutHeader.mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
163+
lutHeader.field = field;
164+
// nch
165+
lutHeader.nchmap.log = true;
166+
lutHeader.nchmap.nbins = 1;
167+
lutHeader.nchmap.min = 0.;
168+
lutHeader.nchmap.max = 4.;
169+
// radius
170+
lutHeader.radmap.log = false;
171+
lutHeader.radmap.nbins = 1;
172+
lutHeader.radmap.min = 0.;
173+
lutHeader.radmap.max = 100.;
174+
// eta
175+
lutHeader.etamap.log = false;
176+
lutHeader.etamap.nbins = histo_eta_bins.n;
177+
lutHeader.etamap.min = histo_eta_bins.min;
178+
lutHeader.etamap.max = histo_eta_bins.max;
179+
Printf("LUT eta: %i, [%f, %f]", lutHeader.etamap.nbins, lutHeader.etamap.min, lutHeader.etamap.max);
180+
// pt
181+
lutHeader.ptmap.log = histo_pt_bins.log;
182+
lutHeader.ptmap.nbins = histo_pt_bins.n;
183+
lutHeader.ptmap.min = histo_pt_bins.log ? std::log10(histo_pt_bins.min) : histo_pt_bins.min;
184+
lutHeader.ptmap.max = histo_pt_bins.log ? std::log10(histo_pt_bins.max) : histo_pt_bins.max;
185+
Printf("LUT pt: %i, [%f, %f]%s", lutHeader.ptmap.nbins, lutHeader.ptmap.min, lutHeader.ptmap.max, lutHeader.ptmap.log ? " LOG AXIS!" : "");
186+
lutFile.write(reinterpret_cast<char*>(&lutHeader), sizeof(lutHeader));
187+
// entries
188+
const int nnch = lutHeader.nchmap.nbins;
189+
const int nrad = lutHeader.radmap.nbins;
190+
const int neta = lutHeader.etamap.nbins;
191+
const int npt = lutHeader.ptmap.nbins;
192+
lutEntry_t lutEntry;
193+
194+
auto resetCovM = [&lutEntry]() {
195+
lutEntry.valid = false;
196+
for (int i = 0; i < 15; ++i) {
197+
lutEntry.covm[i] = 0.;
198+
}
199+
};
200+
201+
TH1F* hptcalls = (TH1F*)h["pt"]->Clone("hptcalls");
202+
hptcalls->Reset();
203+
204+
TH1F* hetacalls = (TH1F*)h["eta"]->Clone("hetacalls");
205+
hetacalls->Reset();
206+
207+
// write entries
208+
for (int inch = 0; inch < nnch; ++inch) {
209+
for (int irad = 0; irad < nrad; ++irad) {
210+
for (int ieta = 0; ieta < neta; ++ieta) {
211+
lutEntry.eta = lutHeader.etamap.eval(ieta);
212+
hetacalls->Fill(lutEntry.eta);
213+
const int bin_eta = h["eta"]->FindBin(lutEntry.eta);
214+
if (ieta == 0 && bin_eta != 1) {
215+
Printf("First eta is not the first bin");
216+
return;
217+
}
218+
if (lutEntry.eta < minEta || lutEntry.eta > maxEta) {
219+
continue;
220+
}
221+
for (int ipt = 0; ipt < npt; ++ipt) {
222+
lutEntry.pt = lutHeader.ptmap.eval(ipt);
223+
hptcalls->Fill(lutEntry.pt);
224+
const int bin_pt = h["pt"]->FindBin(lutEntry.pt);
225+
if (ipt == 0 && bin_pt != 1) {
226+
Printf("First pt is not the first bin");
227+
return;
228+
}
229+
if (lutEntry.pt < minPt || lutEntry.pt > maxPt) {
230+
continue;
231+
}
232+
if (bin_eta <= 0 || bin_eta > h["eta"]->GetNbinsX()) {
233+
resetCovM();
234+
} else if (h["eta"]->GetBinContent(bin_eta) <= 0.f) {
235+
resetCovM();
236+
} else if (bin_pt <= 0 || bin_pt > h["pt"]->GetNbinsX()) {
237+
resetCovM();
238+
} else if (h["pt"]->GetBinContent(bin_pt) <= 0.f) {
239+
resetCovM();
240+
} else {
241+
if (fabs(lutEntry.eta) < .3) { // Barrel
242+
// const int bin = m["Efficiency"]->FindBin(lutEntry.pt, lutEntry.eta, 3.14);
243+
const int bin = m["Efficiency"]->FindBin(lutEntry.pt, lutEntry.eta);
244+
lutEntry.eff = m["Efficiency"]->GetBinContent(bin);
245+
lutEntry.covm[0] = m["CovMat_cYY"]->GetBinContent(bin);
246+
lutEntry.covm[1] = m["CovMat_cZY"]->GetBinContent(bin);
247+
lutEntry.covm[2] = m["CovMat_cZZ"]->GetBinContent(bin);
248+
lutEntry.covm[3] = m["CovMat_cSnpY"]->GetBinContent(bin);
249+
lutEntry.covm[4] = m["CovMat_cSnpZ"]->GetBinContent(bin);
250+
lutEntry.covm[5] = m["CovMat_cSnpSnp"]->GetBinContent(bin);
251+
lutEntry.covm[6] = m["CovMat_cTglY"]->GetBinContent(bin);
252+
lutEntry.covm[7] = m["CovMat_cTglZ"]->GetBinContent(bin);
253+
lutEntry.covm[8] = m["CovMat_cTglSnp"]->GetBinContent(bin);
254+
lutEntry.covm[9] = m["CovMat_cTglTgl"]->GetBinContent(bin);
255+
lutEntry.covm[10] = m["CovMat_c1PtY"]->GetBinContent(bin);
256+
lutEntry.covm[11] = m["CovMat_c1PtZ"]->GetBinContent(bin);
257+
lutEntry.covm[12] = m["CovMat_c1PtSnp"]->GetBinContent(bin);
258+
lutEntry.covm[13] = m["CovMat_c1PtTgl"]->GetBinContent(bin);
259+
lutEntry.covm[14] = m["CovMat_c1Pt21Pt2"]->GetBinContent(bin);
260+
261+
lutEntry.valid = true;
262+
} else {
263+
lutEntry.eff = 1.;
264+
resetCovM();
265+
// // printf(" --- fwdSolve: pt = %f, eta = %f, mass = %f, field=%f \n", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field);
266+
// if (!fwdSolve(lutEntry.covm, lutEntry.pt, lutEntry.eta, lutHeader.mass)) {
267+
// // printf(" --- fwdSolve: error \n");
268+
// }
269+
}
270+
}
271+
diagonalise(lutEntry);
272+
if (lutEntry.valid) {
273+
Printf("Writing valid entry at pT %f and eta %f:", lutEntry.pt, lutEntry.eta);
274+
lutEntry.print();
275+
}
276+
lutFile.write(reinterpret_cast<char*>(&lutEntry), sizeof(lutEntry_t));
277+
}
278+
}
279+
}
280+
}
281+
282+
lutFile.close();
283+
TCanvas* can = new TCanvas();
284+
can->Divide(2);
285+
can->cd(1);
286+
hptcalls->Scale(1. / hetacalls->GetNbinsX());
287+
hptcalls->Draw("HIST");
288+
can->cd(2);
289+
hetacalls->Draw("HIST");
290+
}
291+
292+
void diagonalise(lutEntry_t& lutEntry)
293+
{
294+
// Printf(" --- diagonalise: pt = %f, eta = %f", lutEntry.pt, lutEntry.eta);
295+
TMatrixDSym m(5);
296+
double fcovm[5][5];
297+
for (int i = 0, k = 0; i < 5; ++i)
298+
for (int j = 0; j < i + 1; ++j, ++k) {
299+
fcovm[i][j] = lutEntry.covm[k];
300+
fcovm[j][i] = lutEntry.covm[k];
301+
}
302+
m.SetMatrixArray((double*)fcovm);
303+
TMatrixDSymEigen eigen(m);
304+
// eigenvalues vector
305+
TVectorD eigenVal = eigen.GetEigenValues();
306+
for (int i = 0; i < 5; ++i)
307+
lutEntry.eigval[i] = eigenVal[i];
308+
// eigenvectors matrix
309+
TMatrixD eigenVec = eigen.GetEigenVectors();
310+
for (int i = 0; i < 5; ++i)
311+
for (int j = 0; j < 5; ++j)
312+
lutEntry.eigvec[i][j] = eigenVec[i][j];
313+
// inverse eigenvectors matrix
314+
eigenVec.Invert();
315+
for (int i = 0; i < 5; ++i)
316+
for (int j = 0; j < 5; ++j)
317+
lutEntry.eiginv[i][j] = eigenVec[i][j];
318+
}

0 commit comments

Comments
 (0)