Skip to content

Commit dae1154

Browse files
spucillospucillo
authored
Merge branch 'AliceO2Group:master' into master
2 parents 87fa998 + 5235c64 commit dae1154

File tree

465 files changed

+47877
-19757
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

465 files changed

+47877
-19757
lines changed

‎.github/workflows/labeler.yml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@ jobs:
1616
steps:
1717
- name: Label the PR
1818
id: labeler
19-
uses: actions/labeler@v5
19+
uses: actions/labeler@v6
2020
with:
2121
repo-token: ${{ secrets.GITHUB_TOKEN }}
2222
sync-labels: true

‎.github/workflows/mega-linter.yml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -38,7 +38,7 @@ jobs:
3838
id: ml
3939
# You can override MegaLinter flavor used to have faster performances
4040
# More info at https://megalinter.io/flavors/
41-
uses: oxsecurity/megalinter@v8.8.0
41+
uses: oxsecurity/megalinter@v9.0.1
4242
env:
4343
# All available variables are described in documentation:
4444
# https://megalinter.io/configuration/

‎.github/workflows/stale.yml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@ jobs:
1313
stale:
1414
runs-on: ubuntu-latest
1515
steps:
16-
- uses: actions/stale@v9
16+
- uses: actions/stale@v10
1717
with:
1818
repo-token: ${{ secrets.GITHUB_TOKEN }}
1919
stale-pr-message: 'This PR has not been updated in the last 30 days. Is it still needed? Unless further action is taken, it will be closed in 5 days.'

‎ALICE3/Core/DelphesO2TrackSmearer.cxx‎

Lines changed: 135 additions & 67 deletions
Original file line numberDiff line numberDiff line change
@@ -10,14 +10,14 @@
1010
// or submit itself to any jurisdiction.
1111

1212
///
13-
/// @file DelphesO2TrackSmearer.cxx
14-
/// @brief Porting to O2Physics of DelphesO2 code.
13+
/// \file DelphesO2TrackSmearer.cxx
14+
/// \author Roberto Preghenella
15+
/// \brief Porting to O2Physics of DelphesO2 code.
1516
/// Minimal changes have been made to the original code for adaptation purposes, formatting and commented parts have been considered.
1617
/// Relevant sources:
1718
/// DelphesO2/src/lutCovm.hh https://github.com/AliceO2Group/DelphesO2/blob/master/src/lutCovm.hh
1819
/// DelphesO2/src/TrackSmearer.cc https://github.com/AliceO2Group/DelphesO2/blob/master/src/TrackSmearer.cc
1920
/// DelphesO2/src/TrackSmearer.hh https://github.com/AliceO2Group/DelphesO2/blob/master/src/TrackSmearer.hh
20-
/// @author: Roberto Preghenella
2121
/// @email: preghenella@bo.infn.it
2222
///
2323

@@ -36,6 +36,12 @@
3636

3737
#include "ALICE3/Core/DelphesO2TrackSmearer.h"
3838

39+
#include <CommonConstants/PhysicsConstants.h>
40+
#include <Framework/Logger.h>
41+
42+
#include <map>
43+
#include <string>
44+
3945
namespace o2
4046
{
4147
namespace delphes
@@ -45,35 +51,73 @@ namespace delphes
4551

4652
bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload)
4753
{
48-
auto ipdg = getIndexPDG(pdg);
54+
if (!filename || filename[0] == '\0') {
55+
LOG(info) << " --- No LUT file provided for PDG " << pdg << ". Skipping load.";
56+
return false;
57+
}
58+
const auto ipdg = getIndexPDG(pdg);
59+
LOGF(info, "Will load %s lut file ..: '%s'", getParticleName(pdg), filename);
4960
if (mLUTHeader[ipdg] && !forceReload) {
50-
std::cout << " --- LUT table for PDG " << pdg << " has been already loaded with index " << ipdg << std::endl;
61+
LOG(info) << " --- LUT table for PDG " << pdg << " has been already loaded with index " << ipdg << std::endl;
5162
return false;
5263
}
64+
if (strncmp(filename, "ccdb:", 5) == 0) { // Check if filename starts with "ccdb:"
65+
LOG(info) << " --- LUT file source identified as CCDB.";
66+
std::string path = std::string(filename).substr(5); // Remove "ccdb:" prefix
67+
const std::string outPath = "/tmp/LUTs/";
68+
filename = Form("%s/%s/snapshot.root", outPath.c_str(), path.c_str());
69+
std::ifstream checkFile(filename); // Check if file already exists
70+
if (!checkFile.is_open()) { // File does not exist, retrieve from CCDB
71+
LOG(info) << " --- CCDB source detected for PDG " << pdg << ": " << path;
72+
if (!mCcdbManager) {
73+
LOG(fatal) << " --- CCDB manager not set. Please set it before loading LUT from CCDB.";
74+
}
75+
std::map<std::string, std::string> metadata;
76+
mCcdbManager->getCCDBAccessor().retrieveBlob(path, outPath, metadata, 1);
77+
// Add CCDB handling logic here if needed
78+
LOG(info) << " --- Now retrieving LUT file from CCDB to: " << filename;
79+
} else { // File exists, proceed to load
80+
LOG(info) << " --- LUT file already exists: " << filename << ". Skipping download.";
81+
checkFile.close();
82+
}
83+
return loadTable(pdg, filename, forceReload);
84+
}
85+
5386
mLUTHeader[ipdg] = new lutHeader_t;
5487

5588
std::ifstream lutFile(filename, std::ifstream::binary);
5689
if (!lutFile.is_open()) {
57-
std::cout << " --- cannot open covariance matrix file for PDG " << pdg << ": " << filename << std::endl;
90+
LOG(info) << " --- cannot open covariance matrix file for PDG " << pdg << ": " << filename << std::endl;
5891
delete mLUTHeader[ipdg];
5992
mLUTHeader[ipdg] = nullptr;
6093
return false;
6194
}
6295
lutFile.read(reinterpret_cast<char*>(mLUTHeader[ipdg]), sizeof(lutHeader_t));
6396
if (lutFile.gcount() != sizeof(lutHeader_t)) {
64-
std::cout << " --- troubles reading covariance matrix header for PDG " << pdg << ": " << filename << std::endl;
97+
LOG(info) << " --- troubles reading covariance matrix header for PDG " << pdg << ": " << filename << std::endl;
6598
delete mLUTHeader[ipdg];
6699
mLUTHeader[ipdg] = nullptr;
67100
return false;
68101
}
69102
if (mLUTHeader[ipdg]->version != LUTCOVM_VERSION) {
70-
std::cout << " --- LUT header version mismatch: expected/detected = " << LUTCOVM_VERSION << "/" << mLUTHeader[ipdg]->version << std::endl;
103+
LOG(info) << " --- LUT header version mismatch: expected/detected = " << LUTCOVM_VERSION << "/" << mLUTHeader[ipdg]->version << std::endl;
71104
delete mLUTHeader[ipdg];
72105
mLUTHeader[ipdg] = nullptr;
73106
return false;
74107
}
75-
if (mLUTHeader[ipdg]->pdg != pdg) {
76-
std::cout << " --- LUT header PDG mismatch: expected/detected = " << pdg << "/" << mLUTHeader[ipdg]->pdg << std::endl;
108+
bool specialPdgCase = false;
109+
switch (pdg) { // Handle special cases
110+
case o2::constants::physics::kAlpha: // Special case: Allow Alpha particles to use He3 LUT
111+
specialPdgCase = (mLUTHeader[ipdg]->pdg == o2::constants::physics::kHelium3);
112+
if (specialPdgCase)
113+
LOG(info)
114+
<< " --- Alpha particles (PDG " << pdg << ") will use He3 LUT data (PDG " << mLUTHeader[ipdg]->pdg << ")" << std::endl;
115+
break;
116+
default:
117+
break;
118+
}
119+
if (mLUTHeader[ipdg]->pdg != pdg && !specialPdgCase) {
120+
LOG(info) << " --- LUT header PDG mismatch: expected/detected = " << pdg << "/" << mLUTHeader[ipdg]->pdg << std::endl;
77121
delete mLUTHeader[ipdg];
78122
mLUTHeader[ipdg] = nullptr;
79123
return false;
@@ -93,14 +137,14 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload)
93137
mLUTEntry[ipdg][inch][irad][ieta][ipt] = new lutEntry_t;
94138
lutFile.read(reinterpret_cast<char*>(mLUTEntry[ipdg][inch][irad][ieta][ipt]), sizeof(lutEntry_t));
95139
if (lutFile.gcount() != sizeof(lutEntry_t)) {
96-
std::cout << " --- troubles reading covariance matrix entry for PDG " << pdg << ": " << filename << std::endl;
140+
LOG(info) << " --- troubles reading covariance matrix entry for PDG " << pdg << ": " << filename << std::endl;
97141
return false;
98142
}
99143
}
100144
}
101145
}
102146
}
103-
std::cout << " --- read covariance matrix table for PDG " << pdg << ": " << filename << std::endl;
147+
LOG(info) << " --- read covariance matrix table for PDG " << pdg << ": " << filename << std::endl;
104148
mLUTHeader[ipdg]->print();
105149

106150
lutFile.close();
@@ -123,43 +167,58 @@ lutEntry_t*
123167
// Interpolate if requested
124168
auto fraction = mLUTHeader[ipdg]->nchmap.fracPositionWithinBin(nch);
125169
if (mInterpolateEfficiency) {
126-
if (fraction > 0.5) {
127-
if (mWhatEfficiency == 1) {
128-
if (inch < mLUTHeader[ipdg]->nchmap.nbins - 1) {
129-
interpolatedEff = (1.5f - fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff + (-0.5f + fraction) * mLUTEntry[ipdg][inch + 1][irad][ieta][ipt]->eff;
130-
} else {
131-
interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff;
132-
}
133-
}
134-
if (mWhatEfficiency == 2) {
135-
if (inch < mLUTHeader[ipdg]->nchmap.nbins - 1) {
136-
interpolatedEff = (1.5f - fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2 + (-0.5f + fraction) * mLUTEntry[ipdg][inch + 1][irad][ieta][ipt]->eff2;
137-
} else {
138-
interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2;
139-
}
170+
static constexpr float kFractionThreshold = 0.5f;
171+
if (fraction > kFractionThreshold) {
172+
switch (mWhatEfficiency) {
173+
case 1:
174+
if (inch < mLUTHeader[ipdg]->nchmap.nbins - 1) {
175+
interpolatedEff = (1.5f - fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff + (-0.5f + fraction) * mLUTEntry[ipdg][inch + 1][irad][ieta][ipt]->eff;
176+
} else {
177+
interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff;
178+
}
179+
break;
180+
case 2:
181+
if (inch < mLUTHeader[ipdg]->nchmap.nbins - 1) {
182+
interpolatedEff = (1.5f - fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2 + (-0.5f + fraction) * mLUTEntry[ipdg][inch + 1][irad][ieta][ipt]->eff2;
183+
} else {
184+
interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2;
185+
}
186+
break;
187+
default:
188+
LOG(fatal) << " --- getLUTEntry: unknown efficiency type " << mWhatEfficiency;
140189
}
141190
} else {
142-
float comparisonValue = mLUTHeader[ipdg]->nchmap.log ? log10(nch) : nch;
143-
if (mWhatEfficiency == 1) {
144-
if (inch > 0 && comparisonValue < mLUTHeader[ipdg]->nchmap.max) {
145-
interpolatedEff = (0.5f + fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff + (0.5f - fraction) * mLUTEntry[ipdg][inch - 1][irad][ieta][ipt]->eff;
146-
} else {
147-
interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff;
148-
}
149-
}
150-
if (mWhatEfficiency == 2) {
151-
if (inch > 0 && comparisonValue < mLUTHeader[ipdg]->nchmap.max) {
152-
interpolatedEff = (0.5f + fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2 + (0.5f - fraction) * mLUTEntry[ipdg][inch - 1][irad][ieta][ipt]->eff2;
153-
} else {
154-
interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2;
155-
}
191+
float comparisonValue = mLUTHeader[ipdg]->nchmap.log ? std::log10(nch) : nch;
192+
switch (mWhatEfficiency) {
193+
case 1:
194+
if (inch > 0 && comparisonValue < mLUTHeader[ipdg]->nchmap.max) {
195+
interpolatedEff = (0.5f + fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff + (0.5f - fraction) * mLUTEntry[ipdg][inch - 1][irad][ieta][ipt]->eff;
196+
} else {
197+
interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff;
198+
}
199+
break;
200+
case 2:
201+
if (inch > 0 && comparisonValue < mLUTHeader[ipdg]->nchmap.max) {
202+
interpolatedEff = (0.5f + fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2 + (0.5f - fraction) * mLUTEntry[ipdg][inch - 1][irad][ieta][ipt]->eff2;
203+
} else {
204+
interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2;
205+
}
206+
break;
207+
default:
208+
LOG(fatal) << " --- getLUTEntry: unknown efficiency type " << mWhatEfficiency;
156209
}
157210
}
158211
} else {
159-
if (mWhatEfficiency == 1)
160-
interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff;
161-
if (mWhatEfficiency == 2)
162-
interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2;
212+
switch (mWhatEfficiency) {
213+
case 1:
214+
interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff;
215+
break;
216+
case 2:
217+
interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2;
218+
break;
219+
default:
220+
LOG(fatal) << " --- getLUTEntry: unknown efficiency type " << mWhatEfficiency;
221+
}
163222
}
164223
return mLUTEntry[ipdg][inch][irad][ieta][ipt];
165224
} //;
@@ -172,10 +231,14 @@ bool TrackSmearer::smearTrack(O2Track& o2track, lutEntry_t* lutEntry, float inte
172231
// generate efficiency
173232
if (mUseEfficiency) {
174233
auto eff = 0.;
175-
if (mWhatEfficiency == 1)
176-
eff = lutEntry->eff;
177-
if (mWhatEfficiency == 2)
178-
eff = lutEntry->eff2;
234+
switch (mWhatEfficiency) {
235+
case 1:
236+
eff = lutEntry->eff;
237+
break;
238+
case 2:
239+
eff = lutEntry->eff2;
240+
break;
241+
}
179242
if (mInterpolateEfficiency)
180243
eff = interpolatedEff;
181244
if (gRandom->Uniform() > eff)
@@ -187,26 +250,28 @@ bool TrackSmearer::smearTrack(O2Track& o2track, lutEntry_t* lutEntry, float inte
187250
return false;
188251

189252
// transform params vector and smear
190-
double params_[5];
191-
for (int i = 0; i < 5; ++i) {
253+
static constexpr int kParSize = 5;
254+
double params[kParSize];
255+
for (int i = 0; i < kParSize; ++i) {
192256
double val = 0.;
193-
for (int j = 0; j < 5; ++j)
257+
for (int j = 0; j < kParSize; ++j)
194258
val += lutEntry->eigvec[j][i] * o2track.getParam(j);
195-
params_[i] = gRandom->Gaus(val, sqrt(lutEntry->eigval[i]));
259+
params[i] = gRandom->Gaus(val, std::sqrt(lutEntry->eigval[i]));
196260
}
197261
// transform back params vector
198-
for (int i = 0; i < 5; ++i) {
262+
for (int i = 0; i < kParSize; ++i) {
199263
double val = 0.;
200-
for (int j = 0; j < 5; ++j)
201-
val += lutEntry->eiginv[j][i] * params_[j];
264+
for (int j = 0; j < kParSize; ++j)
265+
val += lutEntry->eiginv[j][i] * params[j];
202266
o2track.setParam(val, i);
203267
}
204268
// should make a sanity check that par[2] sin(phi) is in [-1, 1]
205-
if (fabs(o2track.getParam(2)) > 1.) {
206-
std::cout << " --- smearTrack failed sin(phi) sanity check: " << o2track.getParam(2) << std::endl;
269+
if (std::fabs(o2track.getParam(2)) > 1.) {
270+
LOG(info) << " --- smearTrack failed sin(phi) sanity check: " << o2track.getParam(2) << std::endl;
207271
}
208272
// set covariance matrix
209-
for (int i = 0; i < 15; ++i)
273+
static constexpr int kCovMatSize = 15;
274+
for (int i = 0; i < kCovMatSize; ++i)
210275
o2track.setCov(lutEntry->covm[i], i);
211276
return isReconstructed;
212277
}
@@ -217,8 +282,11 @@ bool TrackSmearer::smearTrack(O2Track& o2track, int pdg, float nch)
217282
{
218283

219284
auto pt = o2track.getPt();
220-
if (abs(pdg) == 1000020030) {
221-
pt *= 2.f;
285+
switch (pdg) {
286+
case o2::constants::physics::kHelium3:
287+
case -o2::constants::physics::kHelium3:
288+
pt *= 2.f;
289+
break;
222290
}
223291
auto eta = o2track.getEta();
224292
float interpolatedEff = 0.0f;
@@ -234,7 +302,7 @@ double TrackSmearer::getPtRes(int pdg, float nch, float eta, float pt)
234302
{
235303
float dummy = 0.0f;
236304
auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy);
237-
auto val = sqrt(lutEntry->covm[14]) * lutEntry->pt;
305+
auto val = std::sqrt(lutEntry->covm[14]) * lutEntry->pt;
238306
return val;
239307
}
240308

@@ -244,9 +312,9 @@ double TrackSmearer::getEtaRes(int pdg, float nch, float eta, float pt)
244312
{
245313
float dummy = 0.0f;
246314
auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy);
247-
auto sigmatgl = sqrt(lutEntry->covm[9]); // sigmatgl2
248-
auto etaRes = fabs(sin(2.0 * atan(exp(-eta)))) * sigmatgl; // propagate tgl to eta uncertainty
249-
etaRes /= lutEntry->eta; // relative uncertainty
315+
auto sigmatgl = std::sqrt(lutEntry->covm[9]); // sigmatgl2
316+
auto etaRes = std::fabs(std::sin(2.0 * std::atan(std::exp(-eta)))) * sigmatgl; // propagate tgl to eta uncertainty
317+
etaRes /= lutEntry->eta; // relative uncertainty
250318
return etaRes;
251319
}
252320
/*****************************************************************/
@@ -255,7 +323,7 @@ double TrackSmearer::getAbsPtRes(int pdg, float nch, float eta, float pt)
255323
{
256324
float dummy = 0.0f;
257325
auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy);
258-
auto val = sqrt(lutEntry->covm[14]) * pow(lutEntry->pt, 2);
326+
auto val = std::sqrt(lutEntry->covm[14]) * lutEntry->pt * lutEntry->pt;
259327
return val;
260328
}
261329

@@ -265,8 +333,8 @@ double TrackSmearer::getAbsEtaRes(int pdg, float nch, float eta, float pt)
265333
{
266334
float dummy = 0.0f;
267335
auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy);
268-
auto sigmatgl = sqrt(lutEntry->covm[9]); // sigmatgl2
269-
auto etaRes = fabs(sin(2.0 * atan(exp(-eta)))) * sigmatgl; // propagate tgl to eta uncertainty
336+
auto sigmatgl = std::sqrt(lutEntry->covm[9]); // sigmatgl2
337+
auto etaRes = std::fabs(std::sin(2.0 * std::atan(std::exp(-eta)))) * sigmatgl; // propagate tgl to eta uncertainty
270338
return etaRes;
271339
}
272340
/*****************************************************************/

0 commit comments

Comments
 (0)