Skip to content

Speed #11448

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
shahor02 merged 5 commits into from
Jun 9, 2023
Merged

Speed #11448

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
6bbba9a
Code optimization
gvolpe79 Jun 7, 2023
d6133bc
Code optimization
gvolpe79 Jun 7, 2023
0d3fc7a
Minors
gvolpe79 Jun 7, 2023
27cb4c3
Minors
gvolpe79 Jun 7, 2023
9459322
Further optimization
gvolpe79 Jun 7, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 16 additions & 16 deletions DataFormats/Detectors/HMPID/include/DataFormatsHMP/Cluster.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// Copyright 2020-2022 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
Expand Down Expand Up @@ -95,21 +95,21 @@ class Cluster

// public:
protected:
int mCh; // chamber number
int mSi; // size of the formed cluster from which this cluster deduced
int mSt; // flag to mark the quality of the cluster
int mBox; // box contaning this cluster
int mNlocMax; // number of local maxima in formed cluster
int mMaxQpad; // abs pad number of a pad with the highest charge
double mMaxQ; // that max charge value
double mQRaw; // QDC value of the raw cluster
double mQ; // QDC value of the actual cluster
double mErrQ; // error on Q
double mXX; // local x postion, [cm]
double mErrX; // error on x postion, [cm]
double mYY; // local y postion, [cm]
double mErrY; // error on y postion, [cm]
double mChi2; // some estimator of the fit quality
int mCh; // chamber number
int mSi; // size of the formed cluster from which this cluster deduced
int mSt; // flag to mark the quality of the cluster
int mBox; // box contaning this cluster
int mNlocMax; // number of local maxima in formed cluster
int mMaxQpad; // abs pad number of a pad with the highest charge
double mMaxQ; // that max charge value
double mQRaw; // QDC value of the raw cluster
double mQ; // QDC value of the actual cluster
double mErrQ; // error on Q
double mXX; // local x postion, [cm]
double mErrX; // error on x postion, [cm]
double mYY; // local y postion, [cm]
double mErrY; // error on y postion, [cm]
double mChi2; // some estimator of the fit quality
std::vector<const o2::hmpid::Digit*> mDigs; //! list of digits forming this cluster

public:
Expand Down
3 changes: 1 addition & 2 deletions DataFormats/Detectors/HMPID/include/DataFormatsHMP/Digit.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// Copyright 2020-2022 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
Expand All @@ -9,7 +9,6 @@
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.

///
/// \file Digit.h
/// \author Antonio Franco - INFN Bari
/// \version 1.0
Expand Down
77 changes: 42 additions & 35 deletions DataFormats/Detectors/HMPID/src/Cluster.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// Copyright 2020-2022 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
Expand Down Expand Up @@ -167,7 +167,8 @@ void Cluster::corrSin()
int py;
o2::hmpid::Param::lors2Pad(mXX, mYY, pc, px, py); // tmp digit to get it center
double x = mXX - mParam->lorsX(pc, px); // diff between cluster x and center of the pad contaning this cluster
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;
double xpi8on10 = M_PI / 0.8 * x;
mXX += 3.31267e-2 * sin(2.0 * xpi8on10) - 2.66575e-3 * sin(4.0 * xpi8on10) + 2.80553e-3 * sin(6.0 * xpi8on10) + 0.0070;
return;
}

Expand Down Expand Up @@ -197,54 +198,58 @@ void Cluster::fitFunc(int& iNpars, double* deriv, double& chi2, double* par, int
for (int i = 0; i < nPads; i++) { // loop on all pads of the cluster
double dQpadMath = 0;
for (int j = 0; j < iNshape; j++) { // Mathiesons loop as all of them may contribute to this pad
double fracMathi = o2::hmpid::Digit::intMathieson(par[3 * j], par[3 * j + 1], pClu->dig(i)->getPadID());
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
int baseOff = 3 * j;
int baseOff1 = baseOff + 1;
int baseOff2 = baseOff + 2;
double fracMathi = o2::hmpid::Digit::intMathieson(par[baseOff], par[baseOff1], pClu->dig(i)->getPadID());
dQpadMath += par[baseOff2] * fracMathi; // par[3*j+2] is charge par[3*j] is x par[3*j+1] is y of current Mathieson
}
if (dQpadMath > 0 && pClu->dig(i)->mQ > 0) {
chi2 += std::pow((pClu->dig(i)->mQ - dQpadMath), 2.0) / pClu->dig(i)->mQ; // chi2 function to be minimized
}
}
//---calculate gradients...
if (iflag == 2) {
float** derivPart;
derivPart = new float*[iNpars];
for (int j = 0; j < iNpars; j++) {
deriv[j] = 0;
derivPart[j] = new float[nPads];
for (int i = 0; i < nPads; i++) {
derivPart[j][i] = 0;
}
}
for (int i = 0; i < nPads; i++) { // loop on all pads of the cluster
std::vector<std::vector<float>> derivPart(iNpars, std::vector<float>(nPads, 0.0f));
double dHalfPadX = o2::hmpid::Param::sizeHalfPadX();
double dHalfPadY = o2::hmpid::Param::sizeHalfPadY();

for (int i = 0; i < nPads; i++) { // loop on all pads of the cluster
int iPadId = pClu->dig(i)->getPadID();
double lx = o2::hmpid::Digit::lorsX(iPadId);
double ly = o2::hmpid::Digit::lorsY(iPadId);
for (int j = 0; j < iNshape; j++) { // Mathiesons loop as all of them may contribute to this pad
double fracMathi = o2::hmpid::Digit::intMathieson(par[3 * j], par[3 * j + 1], pClu->dig(i)->getPadID());
double lx = o2::hmpid::Digit::lorsX(pClu->dig(i)->getPadID());
double ly = o2::hmpid::Digit::lorsY(pClu->dig(i)->getPadID());
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())) *
o2::hmpid::Digit::intPartMathiY(par[3 * j + 1], pClu->dig(i)->getPadID());
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())) *
o2::hmpid::Digit::intPartMathiX(par[3 * j], pClu->dig(i)->getPadID());
derivPart[3 * j + 2][i] += fracMathi;
int baseOff = 3 * j;
int baseOff1 = baseOff + 1;
int baseOff2 = baseOff + 2;
double fracMathi = o2::hmpid::Digit::intMathieson(par[baseOff], par[baseOff1], iPadId);
derivPart[baseOff][i] += par[baseOff2] * (o2::hmpid::Digit::mathiesonX(par[baseOff] - lx - dHalfPadX) - o2::hmpid::Digit::mathiesonX(par[baseOff] - lx + dHalfPadX)) *
o2::hmpid::Digit::intPartMathiY(par[baseOff1], iPadId);
derivPart[baseOff1][i] += par[baseOff2] * (o2::hmpid::Digit::mathiesonY(par[baseOff1] - ly - dHalfPadY) - o2::hmpid::Digit::mathiesonY(par[baseOff1] - ly + dHalfPadY)) *
o2::hmpid::Digit::intPartMathiX(par[baseOff], iPadId);
derivPart[baseOff2][i] += fracMathi;
}
}
// loop on all pads of the cluster
for (int i = 0; i < nPads; i++) { // loop on all pads of the cluster
double dQpadMath = 0; // pad charge collector
for (int i = 0; i < nPads; i++) { // loop on all pads of the cluster
int iPadId = pClu->dig(i)->getPadID();
double dPadmQ = pClu->dig(i)->mQ;
double dQpadMath = 0.0; // pad charge collector
double twoOverMq = 2.0 / dPadmQ;
for (int j = 0; j < iNshape; j++) { // Mathiesons loop as all of them may contribute to this pad
float fracMathi = o2::hmpid::Digit::intMathieson(par[3 * j], par[3 * j + 1], pClu->dig(i)->getPadID());
dQpadMath += par[3 * j + 2] * fracMathi;
if (dQpadMath > 0 && pClu->dig(i)->mQ > 0) {
deriv[3 * j] += 2 / pClu->dig(i)->mQ * (pClu->dig(i)->mQ - dQpadMath) * derivPart[3 * j][i];
deriv[3 * j + 1] += 2 / pClu->dig(i)->mQ * (pClu->dig(i)->mQ - dQpadMath) * derivPart[3 * j + 1][i];
deriv[3 * j + 2] += 2 / pClu->dig(i)->mQ * (pClu->dig(i)->mQ - dQpadMath) * derivPart[3 * j + 2][i];
int baseOff = 3 * j;
int baseOff1 = baseOff + 1;
int baseOff2 = baseOff + 2;
double fracMathi = o2::hmpid::Digit::intMathieson(par[baseOff], par[baseOff1], iPadId);
dQpadMath += par[baseOff2] * fracMathi;
if (dQpadMath > 0 && dPadmQ > 0) {
double appoggio = twoOverMq * (dPadmQ - dQpadMath);
deriv[baseOff] += appoggio * derivPart[baseOff][i];
deriv[baseOff1] += appoggio * derivPart[baseOff1][i];
deriv[baseOff2] += appoggio * derivPart[baseOff2][i];
}
}
}
// delete array...
for (int i = 0; i < iNpars; i++) {
delete[] derivPart[i];
}
delete[] derivPart;
} //---gradient calculations ended
// fit ended. Final calculations
} // FitFunction()
Expand Down Expand Up @@ -321,6 +326,7 @@ int Cluster::solve(std::vector<o2::hmpid::Cluster>* pCluLst, float* pSigmaCut, b
// Arguments: pCluLst - cluster list pointer where to add new cluster(s)
// isTryUnfold - flag to switch on/off unfolding
// Returns: number of local maxima of original cluster

const int kMaxLocMax = 6; // max allowed number of loc max for fitting
coG(); // First calculate CoG for the given cluster
int iCluCnt = pCluLst->size(); // get current number of clusters already stored in the list by previous operations
Expand Down Expand Up @@ -497,6 +503,7 @@ void Cluster::digAdd(const Digit* pDig)
// Adds a given digit to the list of digits belonging to this cluster, cluster is not owner of digits
// Arguments: pDig - pointer to digit to be added
// Returns: none

if (mDigs.size() == 0) { // create list of digits in the first invocation
mSi = 0;
// std::vector<o2::hmpid::Digit*> fDigs;
Expand Down
10 changes: 5 additions & 5 deletions DataFormats/Detectors/HMPID/src/Digit.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// Copyright 2020-2022 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
Expand Down Expand Up @@ -325,8 +325,8 @@ Double_t Digit::qdcTot(Double_t e, Double_t time, Int_t pc, Int_t px, Int_t py,
/// @return : a charge fraction [0-1] imposed into the pad
double Digit::intPartMathiX(double x, int pad)
{
double shift1 = -lorsX(pad) + 0.5 * o2::hmpid::Param::sizePadX();
double shift2 = -lorsX(pad) - 0.5 * o2::hmpid::Param::sizePadX();
double shift1 = -lorsX(pad) + o2::hmpid::Param::sizeHalfPadX();
double shift2 = -lorsX(pad) - o2::hmpid::Param::sizeHalfPadX();

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

double uy1 = o2::hmpid::Param::sqrtK3y() * tanh(o2::hmpid::Param::k2y() * (y + shift1) / o2::hmpid::Param::pitchAnodeCathode());
double uy2 = o2::hmpid::Param::sqrtK3y() * tanh(o2::hmpid::Param::k2y() * (y + shift2) / o2::hmpid::Param::pitchAnodeCathode());
Expand Down
74 changes: 15 additions & 59 deletions Detectors/HMPID/base/include/HMPIDBase/Param.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// Copyright 2020-2022 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
Expand All @@ -16,25 +16,11 @@
#include <TMath.h>
#include <TNamed.h> //base class
#include <TGeoManager.h> //Instance()

// ef : to use XYZVector
#include <Math/Vector3D.h> //fields
#include "Math/Vector3D.h" //fields

#include <TVector3.h> //Lors2Mars() Mars2Lors()

/* are these necessary?:
#include <Math/GenVector/Rotation3D.h> //ef
#include "Math/GenVector/RotationX.h" //ef
#include "Math/GenVector/RotationY.h" //ef
#include "Math/GenVector/RotationZ.h" //ef
*/
#include <TVector3.h> //Lors2Mars() Mars2Lors()

class TGeoVolume;
class TGeoHMatrix;

using XYZVector = ROOT::Math::XYZVector;

namespace o2
{
namespace hmpid
Expand All @@ -58,8 +44,7 @@ class Param

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

static Param* instance(); // pointer to Param singleton

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

static float r2d() { return 57.2957795; }
static float sizePadX() { return fgCellX; } // pad size x, [cm]
static float sizePadY() { return fgCellY; } // pad size y, [cm]
static float sizePadX() { return fgCellX; } // pad size x, [cm]
static float sizePadY() { return fgCellY; } // pad size y, [cm]
static float sizeHalfPadX() { return fgHalfCellX; } // half pad size x, [cm]
static float sizeHalfPadY() { return fgHalfCellY; } // half pad size y, [cm]

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

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

// is the point in dead area?
static bool isInDead(float x, float y)
{ // ef : moved function-definition from cxx
// Check is the current point is outside of sensitive area or in dead zones
// Arguments: x,y -position
// Returns: 1 if not in sensitive zone
for (Int_t iPc = 0; iPc < 6; iPc++) {
if (x >= fgkMinPcX[iPc] && x <= fgkMaxPcX[iPc] && y >= fgkMinPcY[iPc] && y <= fgkMaxPcY[iPc]) {
return kFALSE; // in current pc
}
}
return kTRUE;
}

static bool isDeadPad(Int_t padx, Int_t pady, Int_t ch) // ef : moved function-definition from cxx
{

// Check is the current pad is active or not
// Arguments: padx,pady pad integer coord
// Returns: kTRUE if dead, kFALSE if active

if (fgMapPad[padx - 1][pady - 1][ch]) {
return kFALSE; // current pad active
}

return kTRUE;
}
static bool isInDead(float x, float y); // is the point in dead area?
static bool isDeadPad(Int_t padx, Int_t pady, Int_t ch); // is a dead pad?

inline void setChStatus(Int_t ch, bool status = kTRUE);
inline void setSectStatus(Int_t ch, Int_t sect, bool status);
Expand Down Expand Up @@ -228,17 +190,12 @@ class Param
double l[3] = {x - mX, y - mY, z};
mM[c]->LocalToMaster(l, m);
}

// template <typename T = double>
TVector3 lors2Mars(Int_t c, double x, double y, Int_t pl = kPc) const
{
double m[3];
lors2Mars(c, x, y, m, pl);

return TVector3(m[0], m[1], m[2]); // TVector3(m);

return TVector3(m);
} // MRS->LRS

void mars2Lors(Int_t c, double* m, double& x, double& y) const
{
double l[3];
Expand All @@ -255,13 +212,12 @@ class Param
ph = TMath::ATan2(l[1], l[0]);
}
void lors2MarsVec(Int_t c, double* m, double* l) const { mM[c]->LocalToMasterVect(m, l); } // LRS->MRS

TVector3 norm(Int_t c) const // TVector3
TVector3 norm(Int_t c) const
{
double n[3];
norm(c, n);
return TVector3(n[0], n[1], n[2]); // TVector3(n);
} // norm
return TVector3(n);
} // norm
void norm(Int_t c, double* n) const
{
double l[3] = {0, 0, 1};
Expand Down Expand Up @@ -327,8 +283,8 @@ class Param
static Int_t fgThreshold; // sigma Cut
static bool fgInstanceType; // kTRUE if from geomatry kFALSE if from ideal geometry

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

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

Expand Down
Loading