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
3636
3737#include " ALICE3/Core/DelphesO2TrackSmearer.h"
3838
39+ #include < CommonConstants/PhysicsConstants.h>
3940#include < Framework/Logger.h>
4041
42+ #include < map>
43+ #include < string>
44+
4145namespace o2
4246{
4347namespace delphes
@@ -54,7 +58,7 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload)
5458 const auto ipdg = getIndexPDG (pdg);
5559 LOGF (info, " Will load %s lut file ..: '%s'" , getParticleName (pdg), filename);
5660 if (mLUTHeader [ipdg] && !forceReload) {
57- 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;
5862 return false ;
5963 }
6064 if (strncmp (filename, " ccdb:" , 5 ) == 0 ) { // Check if filename starts with "ccdb:"
@@ -83,26 +87,26 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload)
8387
8488 std::ifstream lutFile (filename, std::ifstream::binary);
8589 if (!lutFile.is_open ()) {
86- 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;
8791 delete mLUTHeader [ipdg];
8892 mLUTHeader [ipdg] = nullptr ;
8993 return false ;
9094 }
9195 lutFile.read (reinterpret_cast <char *>(mLUTHeader [ipdg]), sizeof (lutHeader_t));
9296 if (lutFile.gcount () != sizeof (lutHeader_t)) {
93- 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;
9498 delete mLUTHeader [ipdg];
9599 mLUTHeader [ipdg] = nullptr ;
96100 return false ;
97101 }
98102 if (mLUTHeader [ipdg]->version != LUTCOVM_VERSION) {
99- 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;
100104 delete mLUTHeader [ipdg];
101105 mLUTHeader [ipdg] = nullptr ;
102106 return false ;
103107 }
104108 if (mLUTHeader [ipdg]->pdg != pdg) {
105- 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;
106110 delete mLUTHeader [ipdg];
107111 mLUTHeader [ipdg] = nullptr ;
108112 return false ;
@@ -122,14 +126,14 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload)
122126 mLUTEntry [ipdg][inch][irad][ieta][ipt] = new lutEntry_t;
123127 lutFile.read (reinterpret_cast <char *>(mLUTEntry [ipdg][inch][irad][ieta][ipt]), sizeof (lutEntry_t));
124128 if (lutFile.gcount () != sizeof (lutEntry_t)) {
125- 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;
126130 return false ;
127131 }
128132 }
129133 }
130134 }
131135 }
132- 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;
133137 mLUTHeader [ipdg]->print ();
134138
135139 lutFile.close ();
@@ -152,43 +156,58 @@ lutEntry_t*
152156 // Interpolate if requested
153157 auto fraction = mLUTHeader [ipdg]->nchmap .fracPositionWithinBin (nch);
154158 if (mInterpolateEfficiency ) {
155- if (fraction > 0.5 ) {
156- if (mWhatEfficiency == 1 ) {
157- if (inch < mLUTHeader [ipdg]->nchmap .nbins - 1 ) {
158- interpolatedEff = (1 .5f - fraction) * mLUTEntry [ipdg][inch][irad][ieta][ipt]->eff + (-0 .5f + fraction) * mLUTEntry [ipdg][inch + 1 ][irad][ieta][ipt]->eff ;
159- } else {
160- interpolatedEff = mLUTEntry [ipdg][inch][irad][ieta][ipt]->eff ;
161- }
162- }
163- if (mWhatEfficiency == 2 ) {
164- if (inch < mLUTHeader [ipdg]->nchmap .nbins - 1 ) {
165- interpolatedEff = (1 .5f - fraction) * mLUTEntry [ipdg][inch][irad][ieta][ipt]->eff2 + (-0 .5f + fraction) * mLUTEntry [ipdg][inch + 1 ][irad][ieta][ipt]->eff2 ;
166- } else {
167- interpolatedEff = mLUTEntry [ipdg][inch][irad][ieta][ipt]->eff2 ;
168- }
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 ;
169178 }
170179 } else {
171- float comparisonValue = mLUTHeader [ipdg]->nchmap .log ? log10 (nch) : nch;
172- if (mWhatEfficiency == 1 ) {
173- if (inch > 0 && comparisonValue < mLUTHeader [ipdg]->nchmap .max ) {
174- interpolatedEff = (0 .5f + fraction) * mLUTEntry [ipdg][inch][irad][ieta][ipt]->eff + (0 .5f - fraction) * mLUTEntry [ipdg][inch - 1 ][irad][ieta][ipt]->eff ;
175- } else {
176- interpolatedEff = mLUTEntry [ipdg][inch][irad][ieta][ipt]->eff ;
177- }
178- }
179- if (mWhatEfficiency == 2 ) {
180- if (inch > 0 && comparisonValue < mLUTHeader [ipdg]->nchmap .max ) {
181- interpolatedEff = (0 .5f + fraction) * mLUTEntry [ipdg][inch][irad][ieta][ipt]->eff2 + (0 .5f - fraction) * mLUTEntry [ipdg][inch - 1 ][irad][ieta][ipt]->eff2 ;
182- } else {
183- interpolatedEff = mLUTEntry [ipdg][inch][irad][ieta][ipt]->eff2 ;
184- }
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 ;
185198 }
186199 }
187200 } else {
188- if (mWhatEfficiency == 1 )
189- interpolatedEff = mLUTEntry [ipdg][inch][irad][ieta][ipt]->eff ;
190- if (mWhatEfficiency == 2 )
191- 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+ }
192211 }
193212 return mLUTEntry [ipdg][inch][irad][ieta][ipt];
194213} // ;
@@ -201,10 +220,14 @@ bool TrackSmearer::smearTrack(O2Track& o2track, lutEntry_t* lutEntry, float inte
201220 // generate efficiency
202221 if (mUseEfficiency ) {
203222 auto eff = 0 .;
204- if (mWhatEfficiency == 1 )
205- eff = lutEntry->eff ;
206- if (mWhatEfficiency == 2 )
207- 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+ }
208231 if (mInterpolateEfficiency )
209232 eff = interpolatedEff;
210233 if (gRandom ->Uniform () > eff)
@@ -216,26 +239,28 @@ bool TrackSmearer::smearTrack(O2Track& o2track, lutEntry_t* lutEntry, float inte
216239 return false ;
217240
218241 // transform params vector and smear
219- double params_[5 ];
220- 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) {
221245 double val = 0 .;
222- for (int j = 0 ; j < 5 ; ++j)
246+ for (int j = 0 ; j < kParSize ; ++j)
223247 val += lutEntry->eigvec [j][i] * o2track.getParam (j);
224- params_ [i] = gRandom ->Gaus (val, sqrt (lutEntry->eigval [i]));
248+ params [i] = gRandom ->Gaus (val, std:: sqrt (lutEntry->eigval [i]));
225249 }
226250 // transform back params vector
227- for (int i = 0 ; i < 5 ; ++i) {
251+ for (int i = 0 ; i < kParSize ; ++i) {
228252 double val = 0 .;
229- for (int j = 0 ; j < 5 ; ++j)
230- val += lutEntry->eiginv [j][i] * params_ [j];
253+ for (int j = 0 ; j < kParSize ; ++j)
254+ val += lutEntry->eiginv [j][i] * params [j];
231255 o2track.setParam (val, i);
232256 }
233257 // should make a sanity check that par[2] sin(phi) is in [-1, 1]
234- if (fabs (o2track.getParam (2 )) > 1 .) {
235- 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;
236260 }
237261 // set covariance matrix
238- for (int i = 0 ; i < 15 ; ++i)
262+ static constexpr int kCovMatSize = 15 ;
263+ for (int i = 0 ; i < kCovMatSize ; ++i)
239264 o2track.setCov (lutEntry->covm [i], i);
240265 return isReconstructed;
241266}
@@ -246,8 +271,11 @@ bool TrackSmearer::smearTrack(O2Track& o2track, int pdg, float nch)
246271{
247272
248273 auto pt = o2track.getPt ();
249- if (abs (pdg) == 1000020030 ) {
250- pt *= 2 .f ;
274+ switch (pdg) {
275+ case o2::constants::physics::kHelium3 :
276+ case -o2::constants::physics::kHelium3 :
277+ pt *= 2 .f ;
278+ break ;
251279 }
252280 auto eta = o2track.getEta ();
253281 float interpolatedEff = 0 .0f ;
@@ -263,7 +291,7 @@ double TrackSmearer::getPtRes(int pdg, float nch, float eta, float pt)
263291{
264292 float dummy = 0 .0f ;
265293 auto lutEntry = getLUTEntry (pdg, nch, 0 ., eta, pt, dummy);
266- auto val = sqrt (lutEntry->covm [14 ]) * lutEntry->pt ;
294+ auto val = std:: sqrt (lutEntry->covm [14 ]) * lutEntry->pt ;
267295 return val;
268296}
269297
@@ -273,9 +301,9 @@ double TrackSmearer::getEtaRes(int pdg, float nch, float eta, float pt)
273301{
274302 float dummy = 0 .0f ;
275303 auto lutEntry = getLUTEntry (pdg, nch, 0 ., eta, pt, dummy);
276- auto sigmatgl = sqrt (lutEntry->covm [9 ]); // sigmatgl2
277- auto etaRes = fabs (sin (2.0 * atan (exp (-eta)))) * sigmatgl; // propagate tgl to eta uncertainty
278- 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
279307 return etaRes;
280308}
281309/* ****************************************************************/
@@ -284,7 +312,7 @@ double TrackSmearer::getAbsPtRes(int pdg, float nch, float eta, float pt)
284312{
285313 float dummy = 0 .0f ;
286314 auto lutEntry = getLUTEntry (pdg, nch, 0 ., eta, pt, dummy);
287- auto val = sqrt (lutEntry->covm [14 ]) * pow ( lutEntry->pt , 2 ) ;
315+ auto val = std:: sqrt (lutEntry->covm [14 ]) * lutEntry->pt * lutEntry-> pt ;
288316 return val;
289317}
290318
@@ -294,8 +322,8 @@ double TrackSmearer::getAbsEtaRes(int pdg, float nch, float eta, float pt)
294322{
295323 float dummy = 0 .0f ;
296324 auto lutEntry = getLUTEntry (pdg, nch, 0 ., eta, pt, dummy);
297- auto sigmatgl = sqrt (lutEntry->covm [9 ]); // sigmatgl2
298- 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
299327 return etaRes;
300328}
301329/* ****************************************************************/
0 commit comments