Skip to content

Commit 9d15d37

Browse files
mhemmer-cernmhemmer-cern
committed
[EMCAL-688] Update ClusterFactory
- Remove TMath calls with std function calls - Fix wrong comments for `thetaToEta` and `etaToTheta` in the .cxx - Update `evalDispersion` function to calculate weights, eta and phi per cell only once in a single loop instead of having two loops
1 parent 2aa6d61 commit 9d15d37

2 files changed

Lines changed: 78 additions & 78 deletions

File tree

‎Detectors/EMCAL/base/include/EMCALBase/ClusterFactory.h‎

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -409,12 +409,14 @@ class ClusterFactory
409409
void evalTime(std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const;
410410

411411
///
412-
/// Converts Theta (Radians) to Eta (Radians)
413-
float thetaToEta(float arg) const;
412+
/// \brief Converts Theta (Radians) to Eta (Radians)
413+
/// \param theta theta
414+
float thetaToEta(float theta) const;
414415

415416
///
416-
/// Converts Eta (Radians) to Theta (Radians)
417-
float etaToTheta(float arg) const;
417+
/// \brief Converts Eta (Radians) to Theta (Radians)
418+
/// \param eta eta
419+
float etaToTheta(float eta) const;
418420

419421
private:
420422
o2::emcal::Geometry* mGeomPtr = nullptr;

‎Detectors/EMCAL/base/src/ClusterFactory.cxx‎

Lines changed: 72 additions & 74 deletions
Original file line numberDiff line numberDiff line change
@@ -21,9 +21,13 @@
2121
#include "EMCALBase/Geometry.h"
2222
// #include "MathUtils/Cartesian.h"
2323

24+
#include "CommonConstants/MathConstants.h"
25+
2426
#include <Rtypes.h>
2527

28+
#include <algorithm>
2629
#include <array>
30+
#include <cmath>
2731
#include <span>
2832

2933
using namespace o2::emcal;
@@ -141,81 +145,72 @@ o2::emcal::AnalysisCluster ClusterFactory<InputType>::buildCluster(int clusterIn
141145
}
142146

143147
///
144-
/// Calculates the dispersion of the shower at the origin of the cluster
145-
/// in cell units
148+
/// \brief Calculates the dispersion of the shower at the origin of the cluster in cell units
149+
/// \param inputsIndices span of the input cell Indices
150+
/// \param clusterAnalysis AnalysisCluster for which the elips axis is calculated
146151
//____________________________________________________________________________
147152
template <class InputType>
148153
void ClusterFactory<InputType>::evalDispersion(std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
149154
{
150155
double d = 0., wtot = 0.;
151-
int nstat = 0;
152156

153-
// Calculates the dispersion in cell units
157+
if (clusterAnalysis.E() <= 0) {
158+
clusterAnalysis.setDispersion(0.);
159+
return;
160+
}
161+
162+
struct CellWeight {
163+
double eta, phi, w;
164+
};
165+
std::vector<CellWeight> cellData;
166+
cellData.reserve(inputsIndices.size());
167+
154168
double etaMean = 0.0, phiMean = 0.0;
155169

156-
// Calculate mean values
157170
for (auto iInput : inputsIndices) {
171+
if (mInputsContainer[iInput].getEnergy() <= 0) {
172+
continue;
173+
}
158174

159-
if (clusterAnalysis.E() > 0 && mInputsContainer[iInput].getEnergy() > 0) {
160-
auto [nSupMod, nModule, nIphi, nIeta] = mGeomPtr->GetCellIndex(mInputsContainer[iInput].getTower());
161-
auto [iphi, ieta] = mGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod, nModule, nIphi, nIeta);
175+
auto [nSupMod, nModule, nIphi, nIeta] = mGeomPtr->GetCellIndex(mInputsContainer[iInput].getTower());
176+
auto [iphi, ieta] = mGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod, nModule, nIphi, nIeta);
162177

163-
// In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
164-
// C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
165-
if (mSharedCluster && nSupMod % 2) {
166-
ieta += EMCAL_COLS;
167-
}
178+
// In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
179+
// C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
180+
if (mSharedCluster && nSupMod % 2) {
181+
ieta += EMCAL_COLS;
182+
}
168183

169-
auto etai = static_cast<double>(ieta);
170-
auto phii = static_cast<double>(iphi);
171-
double w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E()));
184+
auto etai = static_cast<double>(ieta);
185+
auto phii = static_cast<double>(iphi);
186+
double w = std::max(0., static_cast<double>(mLogWeight + std::log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E())));
172187

173-
if (w > 0.0) {
174-
phiMean += phii * w;
175-
etaMean += etai * w;
176-
wtot += w;
177-
}
188+
if (w > 0.0) {
189+
cellData.push_back({etai, phii, w});
190+
phiMean += phii * w;
191+
etaMean += etai * w;
192+
wtot += w;
178193
}
179194
}
180195

181196
if (wtot > 0) {
182197
phiMean /= wtot;
183198
etaMean /= wtot;
184199
} else {
185-
LOG(error) << Form("Wrong weight %f\n", wtot);
200+
LOG(error) << "Wrong weight " << wtot;
186201
}
187202

188-
// Calculate dispersion
189-
for (auto iInput : inputsIndices) {
190-
191-
if (clusterAnalysis.E() > 0 && mInputsContainer[iInput].getEnergy() > 0) {
192-
auto [nSupMod, nModule, nIphi, nIeta] = mGeomPtr->GetCellIndex(mInputsContainer[iInput].getTower());
193-
auto [iphi, ieta] = mGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod, nModule, nIphi, nIeta);
194-
195-
// In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
196-
// C Side impair SM, nSupMod%2=1; A side pair SM, nSupMod%2=0
197-
if (mSharedCluster && nSupMod % 2) {
198-
ieta += EMCAL_COLS;
199-
}
200-
201-
auto etai = static_cast<double>(ieta);
202-
auto phii = static_cast<double>(iphi);
203-
double w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E()));
204-
205-
if (w > 0.0) {
206-
nstat++;
207-
d += w * ((etai - etaMean) * (etai - etaMean) + (phii - phiMean) * (phii - phiMean));
208-
}
209-
}
203+
for (const auto& c : cellData) {
204+
d += c.w * ((c.eta - etaMean) * (c.eta - etaMean) + (c.phi - phiMean) * (c.phi - phiMean));
210205
}
211206

212-
if (wtot > 0 && nstat > 1) {
207+
if (wtot > 0 && cellData.size() > 1) {
213208
d /= wtot;
214209
} else {
215210
d = 0.;
216211
}
217212

218-
clusterAnalysis.setDispersion(TMath::Sqrt(d));
213+
clusterAnalysis.setDispersion(std::sqrt(d));
219214
}
220215

221216
///
@@ -247,7 +242,7 @@ void ClusterFactory<InputType>::evalLocalPosition(std::span<const int> inputsInd
247242
}
248243

249244
if (mLogWeight > 0.0) {
250-
w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E()));
245+
w = std::max(0., static_cast<double>(mLogWeight + std::log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E())));
251246
} else {
252247
w = mInputsContainer[iInput].getEnergy(); // just energy
253248
}
@@ -266,7 +261,7 @@ void ClusterFactory<InputType>::evalLocalPosition(std::span<const int> inputsInd
266261
// cout << " wtot " << wtot << endl;
267262

268263
if (wtot > 0) {
269-
// xRMS = TMath::Sqrt(x2m - xMean*xMean);
264+
// xRMS = std::sqrt(x2m - xMean*xMean);
270265
for (int i = 0; i < 3; i++) {
271266
clXYZ[i] /= wtot;
272267

@@ -275,7 +270,7 @@ void ClusterFactory<InputType>::evalLocalPosition(std::span<const int> inputsInd
275270
clRmsXYZ[i] = clRmsXYZ[i] - clXYZ[i] * clXYZ[i];
276271

277272
if (clRmsXYZ[i] > 0.0) {
278-
clRmsXYZ[i] = TMath::Sqrt(clRmsXYZ[i]);
273+
clRmsXYZ[i] = std::sqrt(clRmsXYZ[i]);
279274
} else {
280275
clRmsXYZ[i] = 0;
281276
}
@@ -320,7 +315,7 @@ void ClusterFactory<InputType>::evalGlobalPosition(std::span<const int> inputsIn
320315
mGeomPtr->GetGlobal(lxyzi, xyzi, mGeomPtr->GetSuperModuleNumber(mInputsContainer[iInput].getTower()));
321316

322317
if (mLogWeight > 0.0) {
323-
w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E()));
318+
w = std::max(0., static_cast<double>(mLogWeight + std::log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E())));
324319
} else {
325320
w = mInputsContainer[iInput].getEnergy(); // just energy
326321
}
@@ -339,7 +334,7 @@ void ClusterFactory<InputType>::evalGlobalPosition(std::span<const int> inputsIn
339334
// cout << " wtot " << wtot << endl;
340335

341336
if (wtot > 0) {
342-
// xRMS = TMath::Sqrt(x2m - xMean*xMean);
337+
// xRMS = std::sqrt(x2m - xMean*xMean);
343338
for (i = 0; i < 3; i++) {
344339
clXYZ[i] /= wtot;
345340

@@ -348,7 +343,7 @@ void ClusterFactory<InputType>::evalGlobalPosition(std::span<const int> inputsIn
348343
clRmsXYZ[i] = clRmsXYZ[i] - clXYZ[i] * clXYZ[i];
349344

350345
if (clRmsXYZ[i] > 0.0) {
351-
clRmsXYZ[i] = TMath::Sqrt(clRmsXYZ[i]);
346+
clRmsXYZ[i] = std::sqrt(clRmsXYZ[i]);
352347
} else {
353348
clRmsXYZ[i] = 0;
354349
}
@@ -386,7 +381,7 @@ void ClusterFactory<InputType>::evalLocalPositionFit(double deff, double mLogWei
386381
}
387382

388383
if (mLogWeight > 0.0) {
389-
w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E()));
384+
w = std::max(0., static_cast<double>(mLogWeight + std::log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E())));
390385
} else {
391386
w = mInputsContainer[iInput].getEnergy(); // just energy
392387
}
@@ -405,7 +400,7 @@ void ClusterFactory<InputType>::evalLocalPositionFit(double deff, double mLogWei
405400
// cout << " wtot " << wtot << endl;
406401

407402
if (wtot > 0) {
408-
// xRMS = TMath::Sqrt(x2m - xMean*xMean);
403+
// xRMS = std::sqrt(x2m - xMean*xMean);
409404
for (i = 0; i < 3; i++) {
410405
clXYZ[i] /= wtot;
411406

@@ -414,7 +409,7 @@ void ClusterFactory<InputType>::evalLocalPositionFit(double deff, double mLogWei
414409
clRmsXYZ[i] = clRmsXYZ[i] - clXYZ[i] * clXYZ[i];
415410

416411
if (clRmsXYZ[i] > 0.0) {
417-
clRmsXYZ[i] = TMath::Sqrt(clRmsXYZ[i]);
412+
clRmsXYZ[i] = std::sqrt(clRmsXYZ[i]);
418413
} else {
419414
clRmsXYZ[i] = 0;
420415
}
@@ -458,8 +453,8 @@ void ClusterFactory<InputType>::getDeffW0(const double esum, double& deff, doubl
458453
e = esum < 0.5 ? 0.5 : esum;
459454
e = e > 100. ? 100. : e;
460455

461-
deff = kdp0 + kdp1 * TMath::Log(e);
462-
w0 = kwp0 / (1. + TMath::Exp(kwp1 * (e + kwp2)));
456+
deff = kdp0 + kdp1 * std::log(e);
457+
w0 = kwp0 / (1. + std::exp(kwp1 * (e + kwp2)));
463458
}
464459

465460
///
@@ -485,9 +480,9 @@ void ClusterFactory<InputType>::evalCoreEnergy(std::span<const int> inputsIndice
485480
for (auto iInput : inputsIndices) {
486481

487482
auto [eta, phi] = mGeomPtr->EtaPhiFromIndex(mInputsContainer[iInput].getTower());
488-
phi = phi * TMath::DegToRad();
483+
phi = phi * o2::constants::math::Deg2Rad;
489484

490-
double distance = TMath::Sqrt((eta - etaPoint) * (eta - etaPoint) + (phi - phiPoint) * (phi - phiPoint));
485+
double distance = std::sqrt((eta - etaPoint) * (eta - etaPoint) + (phi - phiPoint) * (phi - phiPoint));
491486

492487
if (distance < mCoreRadius) {
493488
coreEnergy += mInputsContainer[iInput].getEnergy();
@@ -555,8 +550,9 @@ void ClusterFactory<InputType>::evalNExMax(std::span<const int> inputsIndices, A
555550
}
556551

557552
///
558-
/// Calculates the axis of the shower ellipsoid in eta and phi
559-
/// in cell units
553+
/// \brief Calculates the axis of the shower ellipsoid in eta and phi in cell units
554+
/// \param inputsIndices span of the input cell Indices
555+
/// \param clusterAnalysis AnalysisCluster for which the elips axis is calculated
560556
//____________________________________________________________________________
561557
template <class InputType>
562558
void ClusterFactory<InputType>::evalElipsAxis(std::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
@@ -584,7 +580,7 @@ void ClusterFactory<InputType>::evalElipsAxis(std::span<const int> inputsIndices
584580
auto etai = static_cast<double>(ieta);
585581
auto phii = static_cast<double>(iphi);
586582

587-
double w = TMath::Max(0., mLogWeight + TMath::Log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E()));
583+
double w = std::max(0., static_cast<double>(mLogWeight + std::log(mInputsContainer[iInput].getEnergy() / clusterAnalysis.E())));
588584
// clusterAnalysis.E() summed amplitude of inputs, i.e. energy of cluster
589585
// Gives smaller value of lambda than log weight
590586
// w = mEnergyList[iInput] / clusterAnalysis.E(); // Nov 16, 2006 - try just energy
@@ -609,18 +605,18 @@ void ClusterFactory<InputType>::evalElipsAxis(std::span<const int> inputsIndices
609605
dxz /= wtot;
610606
dxz -= x * z;
611607

612-
lambda[0] = 0.5 * (dxx + dzz) + TMath::Sqrt(0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz);
608+
lambda[0] = 0.5 * (dxx + dzz) + std::sqrt(0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz);
613609

614610
if (lambda[0] > 0) {
615-
lambda[0] = TMath::Sqrt(lambda[0]);
611+
lambda[0] = std::sqrt(lambda[0]);
616612
} else {
617613
lambda[0] = 0;
618614
}
619615

620-
lambda[1] = 0.5 * (dxx + dzz) - TMath::Sqrt(0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz);
616+
lambda[1] = 0.5 * (dxx + dzz) - std::sqrt(0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz);
621617

622618
if (lambda[1] > 0) { // To avoid exception if numerical errors lead to negative lambda.
623-
lambda[1] = TMath::Sqrt(lambda[1]);
619+
lambda[1] = std::sqrt(lambda[1]);
624620
} else {
625621
lambda[1] = 0.;
626622
}
@@ -862,13 +858,13 @@ void ClusterFactory<InputType>::evalTime(std::span<const int> inputsIndices, Ana
862858
template <class InputType>
863859
double ClusterFactory<InputType>::tMaxInCm(const double e, const int key) const
864860
{
865-
const double ca = 4.82; // shower max parameter - first guess; ca=TMath::Log(1000./8.07)
861+
const double ca = 4.82; // shower max parameter - first guess; ca=std::log(1000./8.07)
866862
double tmax = 0.; // position of electromagnetic shower max in cm
867863

868864
const double x0 = 1.31; // radiation lenght (cm)
869865

870866
if (e > 0.1) {
871-
tmax = TMath::Log(e) + ca;
867+
tmax = std::log(e) + ca;
872868
if (key == 0) {
873869
tmax += 0.5;
874870
} else {
@@ -881,21 +877,23 @@ double ClusterFactory<InputType>::tMaxInCm(const double e, const int key) const
881877
}
882878

883879
///
884-
/// Converts Theta (Radians) to Eta (Radians)
880+
/// \brief Converts Eta (Radians) to Theta (Radians)
881+
/// \param eta eta
885882
//______________________________________________________________________________
886883
template <class InputType>
887-
float ClusterFactory<InputType>::etaToTheta(float arg) const
884+
float ClusterFactory<InputType>::etaToTheta(float eta) const
888885
{
889-
return (2. * TMath::ATan(TMath::Exp(-arg)));
886+
return (2.f * std::atan(std::exp(-eta)));
890887
}
891888

892889
///
893-
/// Converts Eta (Radians) to Theta (Radians)
890+
/// \brief Converts Theta (Radians) to Eta (Radians)
891+
/// \param theta theta
894892
//______________________________________________________________________________
895893
template <class InputType>
896-
float ClusterFactory<InputType>::thetaToEta(float arg) const
894+
float ClusterFactory<InputType>::thetaToEta(float theta) const
897895
{
898-
return (-1 * TMath::Log(TMath::Tan(0.5 * arg)));
896+
return (-1.f * std::log(std::tan(0.5f * theta)));
899897
}
900898

901899
template <class InputType>

0 commit comments

Comments
 (0)