Skip to content

Commit fc56502

Browse files
njacazionjacazio
committed
Add writer for the LUT from AODs
- Using the binning of the input histograms
1 parent 952ab74 commit fc56502

File tree

2 files changed

+282
-0
lines changed

2 files changed

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

0 commit comments

Comments
 (0)