Skip to content

Commit bbe57b4

Browse files
wenyaCernwenyaCern
authored
Merge branch 'AliceO2Group:master' into master
2 parents 2098a14 + 907006c commit bbe57b4

Some content is hidden

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

46 files changed

+2453
-1486
lines changed

‎ALICE3/Core/DelphesO2TrackSmearer.cxx‎

Lines changed: 123 additions & 66 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,62 @@ 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
}
75108
if (mLUTHeader[ipdg]->pdg != pdg) {
76-
std::cout << " --- LUT header PDG mismatch: expected/detected = " << pdg << "/" << mLUTHeader[ipdg]->pdg << std::endl;
109+
LOG(info) << " --- LUT header PDG mismatch: expected/detected = " << pdg << "/" << mLUTHeader[ipdg]->pdg << std::endl;
77110
delete mLUTHeader[ipdg];
78111
mLUTHeader[ipdg] = nullptr;
79112
return false;
@@ -93,14 +126,14 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload)
93126
mLUTEntry[ipdg][inch][irad][ieta][ipt] = new lutEntry_t;
94127
lutFile.read(reinterpret_cast<char*>(mLUTEntry[ipdg][inch][irad][ieta][ipt]), sizeof(lutEntry_t));
95128
if (lutFile.gcount() != sizeof(lutEntry_t)) {
96-
std::cout << " --- troubles reading covariance matrix entry for PDG " << pdg << ": " << filename << std::endl;
129+
LOG(info) << " --- troubles reading covariance matrix entry for PDG " << pdg << ": " << filename << std::endl;
97130
return false;
98131
}
99132
}
100133
}
101134
}
102135
}
103-
std::cout << " --- read covariance matrix table for PDG " << pdg << ": " << filename << std::endl;
136+
LOG(info) << " --- read covariance matrix table for PDG " << pdg << ": " << filename << std::endl;
104137
mLUTHeader[ipdg]->print();
105138

106139
lutFile.close();
@@ -123,43 +156,58 @@ lutEntry_t*
123156
// Interpolate if requested
124157
auto fraction = mLUTHeader[ipdg]->nchmap.fracPositionWithinBin(nch);
125158
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-
}
159+
static constexpr float kFractionThreshold = 0.5f;
160+
if (fraction > kFractionThreshold) {
161+
switch (mWhatEfficiency) {
162+
case 1:
163+
if (inch < mLUTHeader[ipdg]->nchmap.nbins - 1) {
164+
interpolatedEff = (1.5f - fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff + (-0.5f + fraction) * mLUTEntry[ipdg][inch + 1][irad][ieta][ipt]->eff;
165+
} else {
166+
interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff;
167+
}
168+
break;
169+
case 2:
170+
if (inch < mLUTHeader[ipdg]->nchmap.nbins - 1) {
171+
interpolatedEff = (1.5f - fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2 + (-0.5f + fraction) * mLUTEntry[ipdg][inch + 1][irad][ieta][ipt]->eff2;
172+
} else {
173+
interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2;
174+
}
175+
break;
176+
default:
177+
LOG(fatal) << " --- getLUTEntry: unknown efficiency type " << mWhatEfficiency;
140178
}
141179
} 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-
}
180+
float comparisonValue = mLUTHeader[ipdg]->nchmap.log ? std::log10(nch) : nch;
181+
switch (mWhatEfficiency) {
182+
case 1:
183+
if (inch > 0 && comparisonValue < mLUTHeader[ipdg]->nchmap.max) {
184+
interpolatedEff = (0.5f + fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff + (0.5f - fraction) * mLUTEntry[ipdg][inch - 1][irad][ieta][ipt]->eff;
185+
} else {
186+
interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff;
187+
}
188+
break;
189+
case 2:
190+
if (inch > 0 && comparisonValue < mLUTHeader[ipdg]->nchmap.max) {
191+
interpolatedEff = (0.5f + fraction) * mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2 + (0.5f - fraction) * mLUTEntry[ipdg][inch - 1][irad][ieta][ipt]->eff2;
192+
} else {
193+
interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2;
194+
}
195+
break;
196+
default:
197+
LOG(fatal) << " --- getLUTEntry: unknown efficiency type " << mWhatEfficiency;
156198
}
157199
}
158200
} 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;
201+
switch (mWhatEfficiency) {
202+
case 1:
203+
interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff;
204+
break;
205+
case 2:
206+
interpolatedEff = mLUTEntry[ipdg][inch][irad][ieta][ipt]->eff2;
207+
break;
208+
default:
209+
LOG(fatal) << " --- getLUTEntry: unknown efficiency type " << mWhatEfficiency;
210+
}
163211
}
164212
return mLUTEntry[ipdg][inch][irad][ieta][ipt];
165213
} //;
@@ -172,10 +220,14 @@ bool TrackSmearer::smearTrack(O2Track& o2track, lutEntry_t* lutEntry, float inte
172220
// generate efficiency
173221
if (mUseEfficiency) {
174222
auto eff = 0.;
175-
if (mWhatEfficiency == 1)
176-
eff = lutEntry->eff;
177-
if (mWhatEfficiency == 2)
178-
eff = lutEntry->eff2;
223+
switch (mWhatEfficiency) {
224+
case 1:
225+
eff = lutEntry->eff;
226+
break;
227+
case 2:
228+
eff = lutEntry->eff2;
229+
break;
230+
}
179231
if (mInterpolateEfficiency)
180232
eff = interpolatedEff;
181233
if (gRandom->Uniform() > eff)
@@ -187,26 +239,28 @@ bool TrackSmearer::smearTrack(O2Track& o2track, lutEntry_t* lutEntry, float inte
187239
return false;
188240

189241
// transform params vector and smear
190-
double params_[5];
191-
for (int i = 0; i < 5; ++i) {
242+
static constexpr int kParSize = 5;
243+
double params[kParSize];
244+
for (int i = 0; i < kParSize; ++i) {
192245
double val = 0.;
193-
for (int j = 0; j < 5; ++j)
246+
for (int j = 0; j < kParSize; ++j)
194247
val += lutEntry->eigvec[j][i] * o2track.getParam(j);
195-
params_[i] = gRandom->Gaus(val, sqrt(lutEntry->eigval[i]));
248+
params[i] = gRandom->Gaus(val, std::sqrt(lutEntry->eigval[i]));
196249
}
197250
// transform back params vector
198-
for (int i = 0; i < 5; ++i) {
251+
for (int i = 0; i < kParSize; ++i) {
199252
double val = 0.;
200-
for (int j = 0; j < 5; ++j)
201-
val += lutEntry->eiginv[j][i] * params_[j];
253+
for (int j = 0; j < kParSize; ++j)
254+
val += lutEntry->eiginv[j][i] * params[j];
202255
o2track.setParam(val, i);
203256
}
204257
// 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;
258+
if (std::fabs(o2track.getParam(2)) > 1.) {
259+
LOG(info) << " --- smearTrack failed sin(phi) sanity check: " << o2track.getParam(2) << std::endl;
207260
}
208261
// set covariance matrix
209-
for (int i = 0; i < 15; ++i)
262+
static constexpr int kCovMatSize = 15;
263+
for (int i = 0; i < kCovMatSize; ++i)
210264
o2track.setCov(lutEntry->covm[i], i);
211265
return isReconstructed;
212266
}
@@ -217,8 +271,11 @@ bool TrackSmearer::smearTrack(O2Track& o2track, int pdg, float nch)
217271
{
218272

219273
auto pt = o2track.getPt();
220-
if (abs(pdg) == 1000020030) {
221-
pt *= 2.f;
274+
switch (pdg) {
275+
case o2::constants::physics::kHelium3:
276+
case -o2::constants::physics::kHelium3:
277+
pt *= 2.f;
278+
break;
222279
}
223280
auto eta = o2track.getEta();
224281
float interpolatedEff = 0.0f;
@@ -234,7 +291,7 @@ double TrackSmearer::getPtRes(int pdg, float nch, float eta, float pt)
234291
{
235292
float dummy = 0.0f;
236293
auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy);
237-
auto val = sqrt(lutEntry->covm[14]) * lutEntry->pt;
294+
auto val = std::sqrt(lutEntry->covm[14]) * lutEntry->pt;
238295
return val;
239296
}
240297

@@ -244,9 +301,9 @@ double TrackSmearer::getEtaRes(int pdg, float nch, float eta, float pt)
244301
{
245302
float dummy = 0.0f;
246303
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
304+
auto sigmatgl = std::sqrt(lutEntry->covm[9]); // sigmatgl2
305+
auto etaRes = std::fabs(std::sin(2.0 * std::atan(std::exp(-eta)))) * sigmatgl; // propagate tgl to eta uncertainty
306+
etaRes /= lutEntry->eta; // relative uncertainty
250307
return etaRes;
251308
}
252309
/*****************************************************************/
@@ -255,7 +312,7 @@ double TrackSmearer::getAbsPtRes(int pdg, float nch, float eta, float pt)
255312
{
256313
float dummy = 0.0f;
257314
auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy);
258-
auto val = sqrt(lutEntry->covm[14]) * pow(lutEntry->pt, 2);
315+
auto val = std::sqrt(lutEntry->covm[14]) * lutEntry->pt * lutEntry->pt;
259316
return val;
260317
}
261318

@@ -265,8 +322,8 @@ double TrackSmearer::getAbsEtaRes(int pdg, float nch, float eta, float pt)
265322
{
266323
float dummy = 0.0f;
267324
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
325+
auto sigmatgl = std::sqrt(lutEntry->covm[9]); // sigmatgl2
326+
auto etaRes = std::fabs(std::sin(2.0 * std::atan(std::exp(-eta)))) * sigmatgl; // propagate tgl to eta uncertainty
270327
return etaRes;
271328
}
272329
/*****************************************************************/

‎ALICE3/Core/DelphesO2TrackSmearer.h‎

Lines changed: 38 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -24,12 +24,14 @@
2424
#ifndef ALICE3_CORE_DELPHESO2TRACKSMEARER_H_
2525
#define ALICE3_CORE_DELPHESO2TRACKSMEARER_H_
2626

27-
#include <map>
28-
#include <iostream>
29-
#include <fstream>
27+
#include <CCDB/BasicCCDBManager.h>
28+
#include <ReconstructionDataFormats/Track.h>
3029

31-
#include "TRandom.h"
32-
#include "ReconstructionDataFormats/Track.h"
30+
#include <TRandom.h>
31+
32+
#include <fstream>
33+
#include <iostream>
34+
#include <map>
3335

3436
///////////////////////////////
3537
/// DelphesO2/src/lutCovm.hh //
@@ -85,7 +87,7 @@ struct map_t {
8587
if (bin > nbins - 1)
8688
return nbins - 1;
8789
return bin;
88-
} //;
90+
} //;
8991
void print() { printf("nbins = %d, min = %f, max = %f, log = %s \n", nbins, min, max, log ? "on" : "off"); } //;
9092
};
9193

@@ -214,10 +216,34 @@ class TrackSmearer
214216
return 7; // Helium3
215217
default:
216218
return 2; // Default: pion
217-
} //;
218-
} //;
219+
}
220+
}
219221

220-
void setdNdEta(float val) { mdNdEta = val; } //;
222+
const char* getParticleName(int pdg)
223+
{
224+
switch (abs(pdg)) {
225+
case 11:
226+
return "electron";
227+
case 13:
228+
return "muon";
229+
case 211:
230+
return "pion";
231+
case 321:
232+
return "kaon";
233+
case 2212:
234+
return "proton";
235+
case 1000010020:
236+
return "deuteron";
237+
case 1000010030:
238+
return "triton";
239+
case 1000020030:
240+
return "helium3";
241+
default:
242+
return "pion"; // Default: pion
243+
}
244+
}
245+
void setdNdEta(float val) { mdNdEta = val; } //;
246+
void setCcdbManager(o2::ccdb::BasicCCDBManager* mgr) { mCcdbManager = mgr; } //;
221247

222248
protected:
223249
static constexpr unsigned int nLUTs = 8; // Number of LUT available
@@ -228,6 +254,9 @@ class TrackSmearer
228254
bool mSkipUnreconstructed = true; // don't smear tracks that are not reco'ed
229255
int mWhatEfficiency = 1;
230256
float mdNdEta = 1600.;
257+
258+
private:
259+
o2::ccdb::BasicCCDBManager* mCcdbManager = nullptr;
231260
};
232261

233262
} // namespace delphes

0 commit comments

Comments
 (0)