Skip to content

Commit fb8d91b

Browse files
aalkinaalkin
committed
add the rest of the methods
1 parent 63df43b commit fb8d91b

File tree

2 files changed

+231
-8
lines changed

2 files changed

+231
-8
lines changed

‎ALICE3/Core/FlatTrackSmearer.cxx‎

Lines changed: 230 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,8 @@
1515

1616
#include <Framework/RuntimeError.h>
1717

18+
#include <TRandom.h>
19+
1820
namespace o2::delphes
1921
{
2022
int TrackSmearer::getIndexPDG(int pdg) const
@@ -69,6 +71,15 @@ const char* TrackSmearer::getParticleName(int pdg) const
6971
}
7072
}
7173

74+
void TrackSmearer::setWhatEfficiency(int val)
75+
{
76+
// FIXME: this really should be an enum
77+
if (val > 2) {
78+
throw framework::runtime_error_f("getLUTEntry: unknown efficiency type %d", mWhatEfficiency);
79+
}
80+
mWhatEfficiency = val;
81+
}
82+
7283
bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload)
7384
{
7485
if (!filename || filename[0] == '\0') {
@@ -102,10 +113,9 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload)
102113
LOGF(error, "%s", framework::error_from_ref(ref).what);
103114
return false;
104115
}
105-
mHeaders[ipdg] = &mLUTData[ipdg].getHeaderRef();
106116

107117
LOGF(info, "Successfully read LUT for PDG %d: %s", pdg, localFilename.c_str());
108-
mHeaders[ipdg]->print();
118+
mLUTData[ipdg].getHeaderRef().print();
109119
return true;
110120
}
111121

@@ -126,10 +136,9 @@ bool TrackSmearer::adoptTable(int pdg, const uint8_t* buffer, size_t size, bool
126136
} catch (framework::RuntimeErrorRef ref) {
127137
LOGF(error, "%s", framework::error_from_ref(ref).what);
128138
}
129-
mHeaders[ipdg] = &mLUTData[ipdg].getHeaderRef();
130139

131140
LOGF(info, "Successfully adopted LUT for PDG %d", pdg);
132-
mHeaders[ipdg]->print();
141+
mLUTData[ipdg].getHeaderRef().print();
133142
return true;
134143
}
135144

@@ -150,10 +159,9 @@ bool TrackSmearer::viewTable(int pdg, const uint8_t* buffer, size_t size, bool f
150159
} catch (framework::RuntimeErrorRef ref) {
151160
LOGF(error, "%s", framework::error_from_ref(ref).what);
152161
}
153-
mHeaders[ipdg] = &mLUTData[ipdg].getHeaderRef();
154162

155163
LOGF(info, "Successfully adopted LUT for PDG %d", pdg);
156-
mHeaders[ipdg]->print();
164+
mLUTData[ipdg].getHeaderRef().print();
157165
return true;
158166
}
159167

@@ -181,5 +189,221 @@ bool TrackSmearer::checkSpecialCase(int pdg, lutHeader_t const& header)
181189
return specialPdgCase;
182190
}
183191

192+
const lutHeader_t* TrackSmearer::getLUTHeader(int pdg) const
193+
{
194+
const int ipdg = getIndexPDG(pdg);
195+
if (!mLUTData[ipdg].isLoaded()) {
196+
return nullptr;
197+
}
198+
return &mLUTData[ipdg].getHeaderRef();
199+
}
200+
201+
const lutEntry_t* TrackSmearer::getLUTEntry(const int pdg, const float nch, const float radius, const float eta, const float pt, float& interpolatedEff) const
202+
{
203+
const int ipdg = getIndexPDG(pdg);
204+
if (!mLUTData[ipdg].isLoaded()) {
205+
return nullptr;
206+
}
207+
208+
const auto& header = mLUTData[ipdg].getHeaderRef();
209+
210+
auto inch = header.nchmap.find(nch);
211+
auto irad = header.radmap.find(radius);
212+
auto ieta = header.etamap.find(eta);
213+
auto ipt = header.ptmap.find(pt);
214+
215+
// Interpolate efficiency if requested
216+
auto fraction = header.nchmap.fracPositionWithinBin(nch);
217+
if (mInterpolateEfficiency) {
218+
static constexpr float kFractionThreshold = 0.5f;
219+
if (fraction > kFractionThreshold) {
220+
switch (mWhatEfficiency) {
221+
case 1: {
222+
const auto* entry_curr = mLUTData[ipdg].getEntryRef(inch, irad, ieta, ipt);
223+
if (inch < header.nchmap.nbins - 1) {
224+
const auto* entry_next = mLUTData[ipdg].getEntryRef(inch + 1, irad, ieta, ipt);
225+
interpolatedEff = (1.5f - fraction) * entry_curr->eff + (-0.5f + fraction) * entry_next->eff;
226+
} else {
227+
interpolatedEff = entry_curr->eff;
228+
}
229+
break;
230+
}
231+
case 2: {
232+
const auto* entry_curr = mLUTData[ipdg].getEntryRef(inch, irad, ieta, ipt);
233+
if (inch < header.nchmap.nbins - 1) {
234+
const auto* entry_next = mLUTData[ipdg].getEntryRef(inch + 1, irad, ieta, ipt);
235+
interpolatedEff = (1.5f - fraction) * entry_curr->eff2 + (-0.5f + fraction) * entry_next->eff2;
236+
} else {
237+
interpolatedEff = entry_curr->eff2;
238+
}
239+
break;
240+
}
241+
}
242+
} else {
243+
float comparisonValue = header.nchmap.log ? std::log10(nch) : nch;
244+
switch (mWhatEfficiency) {
245+
case 1: {
246+
const auto* entry_curr = mLUTData[ipdg].getEntryRef(inch, irad, ieta, ipt);
247+
if (inch > 0 && comparisonValue < header.nchmap.max) {
248+
const auto* entry_prev = mLUTData[ipdg].getEntryRef(inch - 1, irad, ieta, ipt);
249+
interpolatedEff = (0.5f + fraction) * entry_curr->eff + (0.5f - fraction) * entry_prev->eff;
250+
} else {
251+
interpolatedEff = entry_curr->eff;
252+
}
253+
break;
254+
}
255+
case 2: {
256+
const auto* entry_curr = mLUTData[ipdg].getEntryRef(inch, irad, ieta, ipt);
257+
if (inch > 0 && comparisonValue < header.nchmap.max) {
258+
const auto* entry_prev = mLUTData[ipdg].getEntryRef(inch - 1, irad, ieta, ipt);
259+
interpolatedEff = (0.5f + fraction) * entry_curr->eff2 + (0.5f - fraction) * entry_prev->eff2;
260+
} else {
261+
interpolatedEff = entry_curr->eff2;
262+
}
263+
break;
264+
}
265+
}
266+
}
267+
} else {
268+
const auto* entry = mLUTData[ipdg].getEntryRef(inch, irad, ieta, ipt);
269+
if (entry) {
270+
switch (mWhatEfficiency) {
271+
case 1:
272+
interpolatedEff = entry->eff;
273+
break;
274+
case 2:
275+
interpolatedEff = entry->eff2;
276+
break;
277+
}
278+
}
279+
}
280+
281+
return mLUTData[ipdg].getEntryRef(inch, irad, ieta, ipt);
282+
}
283+
284+
bool TrackSmearer::smearTrack(O2Track& o2track, const lutEntry_t* lutEntry, float interpolatedEff)
285+
{
286+
bool isReconstructed = true;
287+
288+
// Generate efficiency
289+
if (mUseEfficiency) {
290+
auto eff = 0.f;
291+
switch (mWhatEfficiency) {
292+
case 1:
293+
eff = lutEntry->eff;
294+
break;
295+
case 2:
296+
eff = lutEntry->eff2;
297+
break;
298+
}
299+
if (mInterpolateEfficiency) {
300+
eff = interpolatedEff;
301+
}
302+
if (gRandom->Uniform() > eff) {//FIXME: use a fixed RNG instead of whatever ROOT has as a default
303+
isReconstructed = false;
304+
}
305+
}
306+
307+
// Return false already now in case not reco'ed
308+
if (!isReconstructed && mSkipUnreconstructed) {
309+
return false;
310+
}
311+
312+
// Transform params vector and smear
313+
static constexpr int kParSize = 5;
314+
double params[kParSize];
315+
for (int i = 0; i < kParSize; ++i) {
316+
double val = 0.;
317+
for (int j = 0; j < kParSize; ++j) {
318+
val += lutEntry->eigvec[j][i] * o2track.getParam(j);
319+
}
320+
params[i] = gRandom->Gaus(val, std::sqrt(lutEntry->eigval[i]));
321+
}
322+
323+
// Transform back params vector
324+
for (int i = 0; i < kParSize; ++i) {
325+
double val = 0.;
326+
for (int j = 0; j < kParSize; ++j) {
327+
val += lutEntry->eiginv[j][i] * params[j];
328+
}
329+
o2track.setParam(val, i);
330+
}
331+
332+
// Sanity check that par[2] sin(phi) is in [-1, 1]
333+
if (std::fabs(o2track.getParam(2)) > 1.) {
334+
LOGF(warn, "smearTrack failed sin(phi) sanity check: %f", o2track.getParam(2));
335+
}
336+
337+
// Set covariance matrix
338+
static constexpr int kCovMatSize = 15;
339+
for (int i = 0; i < kCovMatSize; ++i) {
340+
o2track.setCov(lutEntry->covm[i], i);
341+
}
342+
343+
return isReconstructed;
344+
}
345+
346+
bool TrackSmearer::smearTrack(O2Track& o2track, int pdg, float nch)
347+
{
348+
auto pt = o2track.getPt();
349+
switch (pdg) {
350+
case o2::constants::physics::kHelium3:
351+
case -o2::constants::physics::kHelium3:
352+
pt *= 2.f;
353+
break;
354+
}
355+
356+
auto eta = o2track.getEta();
357+
float interpolatedEff = 0.0f;
358+
const lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0.f, eta, pt, interpolatedEff);
359+
360+
if (!lutEntry || !lutEntry->valid) {
361+
return false;
362+
}
363+
364+
return smearTrack(o2track, lutEntry, interpolatedEff);
365+
}
366+
367+
double TrackSmearer::getPtRes(const int pdg, const float nch, const float eta, const float pt) const
368+
{
369+
float dummy = 0.0f;
370+
const lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0.f, eta, pt, dummy);
371+
auto val = std::sqrt(lutEntry->covm[14]) * lutEntry->pt;
372+
return val;
373+
}
374+
375+
double TrackSmearer::getEtaRes(const int pdg, const float nch, const float eta, const float pt) const
376+
{
377+
float dummy = 0.0f;
378+
const lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0.f, eta, pt, dummy);
379+
auto sigmatgl = std::sqrt(lutEntry->covm[9]); // sigmatgl2
380+
auto etaRes = std::fabs(std::sin(2.0 * std::atan(std::exp(-eta)))) * sigmatgl; // propagate tgl to eta uncertainty
381+
etaRes /= lutEntry->eta; // relative uncertainty
382+
return etaRes;
383+
}
384+
385+
double TrackSmearer::getAbsPtRes(const int pdg, const float nch, const float eta, const float pt) const
386+
{
387+
float dummy = 0.0f;
388+
const lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0.f, eta, pt, dummy);
389+
auto val = std::sqrt(lutEntry->covm[14]) * lutEntry->pt * lutEntry->pt;
390+
return val;
391+
}
392+
393+
double TrackSmearer::getAbsEtaRes(const int pdg, const float nch, const float eta, const float pt) const
394+
{
395+
float dummy = 0.0f;
396+
const lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0.f, eta, pt, dummy);
397+
auto sigmatgl = std::sqrt(lutEntry->covm[9]); // sigmatgl2
398+
auto etaRes = std::fabs(std::sin(2.0 * std::atan(std::exp(-eta)))) * sigmatgl; // propagate tgl to eta uncertainty
399+
return etaRes;
400+
}
401+
402+
double TrackSmearer::getEfficiency(const int pdg, const float nch, const float eta, const float pt) const
403+
{
404+
float efficiency = 0.0f;
405+
(void)getLUTEntry(pdg, nch, 0.f, eta, pt, efficiency);
406+
return efficiency;
407+
}
184408

185409
} // namespace o2::delphes

‎ALICE3/Core/FlatTrackSmearer.h‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -37,7 +37,7 @@ class TrackSmearer
3737
void useEfficiency(bool val) { mUseEfficiency = val; }
3838
void interpolateEfficiency(bool val) { mInterpolateEfficiency = val; }
3939
void skipUnreconstructed(bool val) { mSkipUnreconstructed = val; }
40-
void setWhatEfficiency(int val) { mWhatEfficiency = val; }
40+
void setWhatEfficiency(int val);
4141

4242
const lutHeader_t* getLUTHeader(int pdg) const;
4343
const lutEntry_t* getLUTEntry(int pdg, float nch, float radius, float eta, float pt, float& interpolatedEff) const;
@@ -60,7 +60,6 @@ class TrackSmearer
6060

6161
protected:
6262
static constexpr unsigned int nLUTs = 9; // Number of LUT available
63-
lutHeader_t const* mHeaders[nLUTs]; // header references for quick access
6463
FlatLutData mLUTData[nLUTs]; // Flat data storage
6564

6665
bool mUseEfficiency = true;

0 commit comments

Comments
 (0)