Skip to content

Commit 100710c

Browse files
gvolpe79gvolpe79
authored
Speed (#11448)
* Code optimization * Code optimization * Minors * Minors * Further optimization
1 parent 11cec9f commit 100710c

7 files changed

Lines changed: 138 additions & 146 deletions

File tree

‎DataFormats/Detectors/HMPID/include/DataFormatsHMP/Cluster.h‎

Lines changed: 16 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
1+
// Copyright 2020-2022 CERN and copyright holders of ALICE O2.
22
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
33
// All rights not expressly granted are reserved.
44
//
@@ -95,21 +95,21 @@ class Cluster
9595

9696
// public:
9797
protected:
98-
int mCh; // chamber number
99-
int mSi; // size of the formed cluster from which this cluster deduced
100-
int mSt; // flag to mark the quality of the cluster
101-
int mBox; // box contaning this cluster
102-
int mNlocMax; // number of local maxima in formed cluster
103-
int mMaxQpad; // abs pad number of a pad with the highest charge
104-
double mMaxQ; // that max charge value
105-
double mQRaw; // QDC value of the raw cluster
106-
double mQ; // QDC value of the actual cluster
107-
double mErrQ; // error on Q
108-
double mXX; // local x postion, [cm]
109-
double mErrX; // error on x postion, [cm]
110-
double mYY; // local y postion, [cm]
111-
double mErrY; // error on y postion, [cm]
112-
double mChi2; // some estimator of the fit quality
98+
int mCh; // chamber number
99+
int mSi; // size of the formed cluster from which this cluster deduced
100+
int mSt; // flag to mark the quality of the cluster
101+
int mBox; // box contaning this cluster
102+
int mNlocMax; // number of local maxima in formed cluster
103+
int mMaxQpad; // abs pad number of a pad with the highest charge
104+
double mMaxQ; // that max charge value
105+
double mQRaw; // QDC value of the raw cluster
106+
double mQ; // QDC value of the actual cluster
107+
double mErrQ; // error on Q
108+
double mXX; // local x postion, [cm]
109+
double mErrX; // error on x postion, [cm]
110+
double mYY; // local y postion, [cm]
111+
double mErrY; // error on y postion, [cm]
112+
double mChi2; // some estimator of the fit quality
113113
std::vector<const o2::hmpid::Digit*> mDigs; //! list of digits forming this cluster
114114

115115
public:

‎DataFormats/Detectors/HMPID/include/DataFormatsHMP/Digit.h‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
1+
// Copyright 2020-2022 CERN and copyright holders of ALICE O2.
22
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
33
// All rights not expressly granted are reserved.
44
//
@@ -9,7 +9,6 @@
99
// granted to it by virtue of its status as an Intergovernmental Organization
1010
// or submit itself to any jurisdiction.
1111

12-
///
1312
/// \file Digit.h
1413
/// \author Antonio Franco - INFN Bari
1514
/// \version 1.0

‎DataFormats/Detectors/HMPID/src/Cluster.cxx‎

Lines changed: 42 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
1+
// Copyright 2020-2022 CERN and copyright holders of ALICE O2.
22
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
33
// All rights not expressly granted are reserved.
44
//
@@ -167,7 +167,8 @@ void Cluster::corrSin()
167167
int py;
168168
o2::hmpid::Param::lors2Pad(mXX, mYY, pc, px, py); // tmp digit to get it center
169169
double x = mXX - mParam->lorsX(pc, px); // diff between cluster x and center of the pad contaning this cluster
170-
mXX += 3.31267e-2 * sin(2 * M_PI / 0.8 * x) - 2.66575e-3 * sin(4 * M_PI / 0.8 * x) + 2.80553e-3 * sin(6 * M_PI / 0.8 * x) + 0.0070;
170+
double xpi8on10 = M_PI / 0.8 * x;
171+
mXX += 3.31267e-2 * sin(2.0 * xpi8on10) - 2.66575e-3 * sin(4.0 * xpi8on10) + 2.80553e-3 * sin(6.0 * xpi8on10) + 0.0070;
171172
return;
172173
}
173174

@@ -197,54 +198,58 @@ void Cluster::fitFunc(int& iNpars, double* deriv, double& chi2, double* par, int
197198
for (int i = 0; i < nPads; i++) { // loop on all pads of the cluster
198199
double dQpadMath = 0;
199200
for (int j = 0; j < iNshape; j++) { // Mathiesons loop as all of them may contribute to this pad
200-
double fracMathi = o2::hmpid::Digit::intMathieson(par[3 * j], par[3 * j + 1], pClu->dig(i)->getPadID());
201-
dQpadMath += par[3 * j + 2] * fracMathi; // par[3*j+2] is charge par[3*j] is x par[3*j+1] is y of current Mathieson
201+
int baseOff = 3 * j;
202+
int baseOff1 = baseOff + 1;
203+
int baseOff2 = baseOff + 2;
204+
double fracMathi = o2::hmpid::Digit::intMathieson(par[baseOff], par[baseOff1], pClu->dig(i)->getPadID());
205+
dQpadMath += par[baseOff2] * fracMathi; // par[3*j+2] is charge par[3*j] is x par[3*j+1] is y of current Mathieson
202206
}
203207
if (dQpadMath > 0 && pClu->dig(i)->mQ > 0) {
204208
chi2 += std::pow((pClu->dig(i)->mQ - dQpadMath), 2.0) / pClu->dig(i)->mQ; // chi2 function to be minimized
205209
}
206210
}
207211
//---calculate gradients...
208212
if (iflag == 2) {
209-
float** derivPart;
210-
derivPart = new float*[iNpars];
211-
for (int j = 0; j < iNpars; j++) {
212-
deriv[j] = 0;
213-
derivPart[j] = new float[nPads];
214-
for (int i = 0; i < nPads; i++) {
215-
derivPart[j][i] = 0;
216-
}
217-
}
218-
for (int i = 0; i < nPads; i++) { // loop on all pads of the cluster
213+
std::vector<std::vector<float>> derivPart(iNpars, std::vector<float>(nPads, 0.0f));
214+
double dHalfPadX = o2::hmpid::Param::sizeHalfPadX();
215+
double dHalfPadY = o2::hmpid::Param::sizeHalfPadY();
216+
217+
for (int i = 0; i < nPads; i++) { // loop on all pads of the cluster
218+
int iPadId = pClu->dig(i)->getPadID();
219+
double lx = o2::hmpid::Digit::lorsX(iPadId);
220+
double ly = o2::hmpid::Digit::lorsY(iPadId);
219221
for (int j = 0; j < iNshape; j++) { // Mathiesons loop as all of them may contribute to this pad
220-
double fracMathi = o2::hmpid::Digit::intMathieson(par[3 * j], par[3 * j + 1], pClu->dig(i)->getPadID());
221-
double lx = o2::hmpid::Digit::lorsX(pClu->dig(i)->getPadID());
222-
double ly = o2::hmpid::Digit::lorsY(pClu->dig(i)->getPadID());
223-
derivPart[3 * j][i] += par[3 * j + 2] * (o2::hmpid::Digit::mathiesonX(par[3 * j] - lx - 0.5 * o2::hmpid::Param::sizePadX()) - o2::hmpid::Digit::mathiesonX(par[3 * j] - lx + 0.5 * o2::hmpid::Param::sizePadX())) *
224-
o2::hmpid::Digit::intPartMathiY(par[3 * j + 1], pClu->dig(i)->getPadID());
225-
derivPart[3 * j + 1][i] += par[3 * j + 2] * (o2::hmpid::Digit::mathiesonY(par[3 * j + 1] - ly - 0.5 * o2::hmpid::Param::sizePadY()) - o2::hmpid::Digit::mathiesonY(par[3 * j + 1] - ly + 0.5 * o2::hmpid::Param::sizePadY())) *
226-
o2::hmpid::Digit::intPartMathiX(par[3 * j], pClu->dig(i)->getPadID());
227-
derivPart[3 * j + 2][i] += fracMathi;
222+
int baseOff = 3 * j;
223+
int baseOff1 = baseOff + 1;
224+
int baseOff2 = baseOff + 2;
225+
double fracMathi = o2::hmpid::Digit::intMathieson(par[baseOff], par[baseOff1], iPadId);
226+
derivPart[baseOff][i] += par[baseOff2] * (o2::hmpid::Digit::mathiesonX(par[baseOff] - lx - dHalfPadX) - o2::hmpid::Digit::mathiesonX(par[baseOff] - lx + dHalfPadX)) *
227+
o2::hmpid::Digit::intPartMathiY(par[baseOff1], iPadId);
228+
derivPart[baseOff1][i] += par[baseOff2] * (o2::hmpid::Digit::mathiesonY(par[baseOff1] - ly - dHalfPadY) - o2::hmpid::Digit::mathiesonY(par[baseOff1] - ly + dHalfPadY)) *
229+
o2::hmpid::Digit::intPartMathiX(par[baseOff], iPadId);
230+
derivPart[baseOff2][i] += fracMathi;
228231
}
229232
}
230233
// loop on all pads of the cluster
231-
for (int i = 0; i < nPads; i++) { // loop on all pads of the cluster
232-
double dQpadMath = 0; // pad charge collector
234+
for (int i = 0; i < nPads; i++) { // loop on all pads of the cluster
235+
int iPadId = pClu->dig(i)->getPadID();
236+
double dPadmQ = pClu->dig(i)->mQ;
237+
double dQpadMath = 0.0; // pad charge collector
238+
double twoOverMq = 2.0 / dPadmQ;
233239
for (int j = 0; j < iNshape; j++) { // Mathiesons loop as all of them may contribute to this pad
234-
float fracMathi = o2::hmpid::Digit::intMathieson(par[3 * j], par[3 * j + 1], pClu->dig(i)->getPadID());
235-
dQpadMath += par[3 * j + 2] * fracMathi;
236-
if (dQpadMath > 0 && pClu->dig(i)->mQ > 0) {
237-
deriv[3 * j] += 2 / pClu->dig(i)->mQ * (pClu->dig(i)->mQ - dQpadMath) * derivPart[3 * j][i];
238-
deriv[3 * j + 1] += 2 / pClu->dig(i)->mQ * (pClu->dig(i)->mQ - dQpadMath) * derivPart[3 * j + 1][i];
239-
deriv[3 * j + 2] += 2 / pClu->dig(i)->mQ * (pClu->dig(i)->mQ - dQpadMath) * derivPart[3 * j + 2][i];
240+
int baseOff = 3 * j;
241+
int baseOff1 = baseOff + 1;
242+
int baseOff2 = baseOff + 2;
243+
double fracMathi = o2::hmpid::Digit::intMathieson(par[baseOff], par[baseOff1], iPadId);
244+
dQpadMath += par[baseOff2] * fracMathi;
245+
if (dQpadMath > 0 && dPadmQ > 0) {
246+
double appoggio = twoOverMq * (dPadmQ - dQpadMath);
247+
deriv[baseOff] += appoggio * derivPart[baseOff][i];
248+
deriv[baseOff1] += appoggio * derivPart[baseOff1][i];
249+
deriv[baseOff2] += appoggio * derivPart[baseOff2][i];
240250
}
241251
}
242252
}
243-
// delete array...
244-
for (int i = 0; i < iNpars; i++) {
245-
delete[] derivPart[i];
246-
}
247-
delete[] derivPart;
248253
} //---gradient calculations ended
249254
// fit ended. Final calculations
250255
} // FitFunction()
@@ -321,6 +326,7 @@ int Cluster::solve(std::vector<o2::hmpid::Cluster>* pCluLst, float* pSigmaCut, b
321326
// Arguments: pCluLst - cluster list pointer where to add new cluster(s)
322327
// isTryUnfold - flag to switch on/off unfolding
323328
// Returns: number of local maxima of original cluster
329+
324330
const int kMaxLocMax = 6; // max allowed number of loc max for fitting
325331
coG(); // First calculate CoG for the given cluster
326332
int iCluCnt = pCluLst->size(); // get current number of clusters already stored in the list by previous operations
@@ -497,6 +503,7 @@ void Cluster::digAdd(const Digit* pDig)
497503
// Adds a given digit to the list of digits belonging to this cluster, cluster is not owner of digits
498504
// Arguments: pDig - pointer to digit to be added
499505
// Returns: none
506+
500507
if (mDigs.size() == 0) { // create list of digits in the first invocation
501508
mSi = 0;
502509
// std::vector<o2::hmpid::Digit*> fDigs;

‎DataFormats/Detectors/HMPID/src/Digit.cxx‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
1+
// Copyright 2020-2022 CERN and copyright holders of ALICE O2.
22
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
33
// All rights not expressly granted are reserved.
44
//
@@ -325,8 +325,8 @@ Double_t Digit::qdcTot(Double_t e, Double_t time, Int_t pc, Int_t px, Int_t py,
325325
/// @return : a charge fraction [0-1] imposed into the pad
326326
double Digit::intPartMathiX(double x, int pad)
327327
{
328-
double shift1 = -lorsX(pad) + 0.5 * o2::hmpid::Param::sizePadX();
329-
double shift2 = -lorsX(pad) - 0.5 * o2::hmpid::Param::sizePadX();
328+
double shift1 = -lorsX(pad) + o2::hmpid::Param::sizeHalfPadX();
329+
double shift2 = -lorsX(pad) - o2::hmpid::Param::sizeHalfPadX();
330330

331331
double ux1 = o2::hmpid::Param::sqrtK3x() * tanh(o2::hmpid::Param::k2x() * (x + shift1) / o2::hmpid::Param::pitchAnodeCathode());
332332
double ux2 = o2::hmpid::Param::sqrtK3x() * tanh(o2::hmpid::Param::k2x() * (x + shift2) / o2::hmpid::Param::pitchAnodeCathode());
@@ -343,8 +343,8 @@ double Digit::intPartMathiX(double x, int pad)
343343
/// @return : a charge fraction [0-1] imposed into the pad
344344
double Digit::intPartMathiY(double y, int pad)
345345
{
346-
double shift1 = -lorsY(pad) + 0.5 * o2::hmpid::Param::sizePadY();
347-
double shift2 = -lorsY(pad) - 0.5 * o2::hmpid::Param::sizePadY();
346+
double shift1 = -lorsY(pad) + o2::hmpid::Param::sizeHalfPadY();
347+
double shift2 = -lorsY(pad) - o2::hmpid::Param::sizeHalfPadY();
348348

349349
double uy1 = o2::hmpid::Param::sqrtK3y() * tanh(o2::hmpid::Param::k2y() * (y + shift1) / o2::hmpid::Param::pitchAnodeCathode());
350350
double uy2 = o2::hmpid::Param::sqrtK3y() * tanh(o2::hmpid::Param::k2y() * (y + shift2) / o2::hmpid::Param::pitchAnodeCathode());

‎Detectors/HMPID/base/include/HMPIDBase/Param.h‎

Lines changed: 15 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
1+
// Copyright 2020-2022 CERN and copyright holders of ALICE O2.
22
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
33
// All rights not expressly granted are reserved.
44
//
@@ -16,25 +16,11 @@
1616
#include <TMath.h>
1717
#include <TNamed.h> //base class
1818
#include <TGeoManager.h> //Instance()
19-
20-
// ef : to use XYZVector
21-
#include <Math/Vector3D.h> //fields
22-
#include "Math/Vector3D.h" //fields
23-
24-
#include <TVector3.h> //Lors2Mars() Mars2Lors()
25-
26-
/* are these necessary?:
27-
#include <Math/GenVector/Rotation3D.h> //ef
28-
#include "Math/GenVector/RotationX.h" //ef
29-
#include "Math/GenVector/RotationY.h" //ef
30-
#include "Math/GenVector/RotationZ.h" //ef
31-
*/
19+
#include <TVector3.h> //Lors2Mars() Mars2Lors()
3220

3321
class TGeoVolume;
3422
class TGeoHMatrix;
3523

36-
using XYZVector = ROOT::Math::XYZVector;
37-
3824
namespace o2
3925
{
4026
namespace hmpid
@@ -58,8 +44,7 @@ class Param
5844

5945
void print(Option_t* opt = "") const; // print current parametrization
6046

61-
static Param* instance(); // pointer to Param singleton
62-
47+
static Param* instance(); // pointer to Param singleton
6348
static Param* instanceNoGeo(); // pointer to Param singleton without geometry.root for MOOD, displays, ...
6449
// geo info
6550
enum EChamberData { kMinCh = 0,
@@ -82,8 +67,10 @@ class Param
8267
kPadSigmaMasked = 20 }; // One can go up to 5 sigma cut, overflow is protected in AliHMPIDCalib
8368

8469
static float r2d() { return 57.2957795; }
85-
static float sizePadX() { return fgCellX; } // pad size x, [cm]
86-
static float sizePadY() { return fgCellY; } // pad size y, [cm]
70+
static float sizePadX() { return fgCellX; } // pad size x, [cm]
71+
static float sizePadY() { return fgCellY; } // pad size y, [cm]
72+
static float sizeHalfPadX() { return fgHalfCellX; } // half pad size x, [cm]
73+
static float sizeHalfPadY() { return fgHalfCellY; } // half pad size y, [cm]
8774

8875
static float sizePcX() { return fgPcX; } // PC size x
8976
static float sizePcY() { return fgPcY; } // PC size y
@@ -115,33 +102,8 @@ class Param
115102

116103
bool getInstType() const { return fgInstanceType; } // return if the instance is from geom or ideal
117104

118-
// is the point in dead area?
119-
static bool isInDead(float x, float y)
120-
{ // ef : moved function-definition from cxx
121-
// Check is the current point is outside of sensitive area or in dead zones
122-
// Arguments: x,y -position
123-
// Returns: 1 if not in sensitive zone
124-
for (Int_t iPc = 0; iPc < 6; iPc++) {
125-
if (x >= fgkMinPcX[iPc] && x <= fgkMaxPcX[iPc] && y >= fgkMinPcY[iPc] && y <= fgkMaxPcY[iPc]) {
126-
return kFALSE; // in current pc
127-
}
128-
}
129-
return kTRUE;
130-
}
131-
132-
static bool isDeadPad(Int_t padx, Int_t pady, Int_t ch) // ef : moved function-definition from cxx
133-
{
134-
135-
// Check is the current pad is active or not
136-
// Arguments: padx,pady pad integer coord
137-
// Returns: kTRUE if dead, kFALSE if active
138-
139-
if (fgMapPad[padx - 1][pady - 1][ch]) {
140-
return kFALSE; // current pad active
141-
}
142-
143-
return kTRUE;
144-
}
105+
static bool isInDead(float x, float y); // is the point in dead area?
106+
static bool isDeadPad(Int_t padx, Int_t pady, Int_t ch); // is a dead pad?
145107

146108
inline void setChStatus(Int_t ch, bool status = kTRUE);
147109
inline void setSectStatus(Int_t ch, Int_t sect, bool status);
@@ -228,17 +190,12 @@ class Param
228190
double l[3] = {x - mX, y - mY, z};
229191
mM[c]->LocalToMaster(l, m);
230192
}
231-
232-
// template <typename T = double>
233193
TVector3 lors2Mars(Int_t c, double x, double y, Int_t pl = kPc) const
234194
{
235195
double m[3];
236196
lors2Mars(c, x, y, m, pl);
237-
238-
return TVector3(m[0], m[1], m[2]); // TVector3(m);
239-
197+
return TVector3(m);
240198
} // MRS->LRS
241-
242199
void mars2Lors(Int_t c, double* m, double& x, double& y) const
243200
{
244201
double l[3];
@@ -255,13 +212,12 @@ class Param
255212
ph = TMath::ATan2(l[1], l[0]);
256213
}
257214
void lors2MarsVec(Int_t c, double* m, double* l) const { mM[c]->LocalToMasterVect(m, l); } // LRS->MRS
258-
259-
TVector3 norm(Int_t c) const // TVector3
215+
TVector3 norm(Int_t c) const
260216
{
261217
double n[3];
262218
norm(c, n);
263-
return TVector3(n[0], n[1], n[2]); // TVector3(n);
264-
} // norm
219+
return TVector3(n);
220+
} // norm
265221
void norm(Int_t c, double* n) const
266222
{
267223
double l[3] = {0, 0, 1};
@@ -327,8 +283,8 @@ class Param
327283
static Int_t fgThreshold; // sigma Cut
328284
static bool fgInstanceType; // kTRUE if from geomatry kFALSE if from ideal geometry
329285

330-
static float fgCellX, fgCellY, fgPcX, fgPcY, fgAllX, fgAllY; // definition of HMPID geometric parameters
331-
Param(bool noGeo); // default ctor is protected to enforce it to be singleton
286+
static float fgCellX, fgCellY, fgHalfCellX, fgHalfCellY, fgPcX, fgPcY, fgAllX, fgAllY; // definition of HMPID geometric parameters
287+
Param(bool noGeo); // default ctor is protected to enforce it to be singleton
332288

333289
static Param* fgInstance; // static pointer to instance of Param singleton
334290

0 commit comments

Comments
 (0)