Skip to content

Commit ad2746d

Browse files
mpuccioMohammadAlTurany
mpuccio
authored and
MohammadAlTurany
committed
Latest track class from shahor02/trackO2 repository
1 parent fc6cc0a commit ad2746d

3 files changed

Lines changed: 467 additions & 483 deletions

File tree

‎Detectors/Base/include/DetectorsBase/Track.h‎

Lines changed: 96 additions & 110 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,8 @@
55
#ifndef ALICEO2_BASE_TRACK
66
#define ALICEO2_BASE_TRACK
77

8-
#include <stdio.h>
8+
#include <iostream>
9+
#include <algorithm>
910
#include <string.h>
1011
#include "DetectorsBase/Constants.h"
1112
#include "DetectorsBase/Utils.h"
@@ -16,39 +17,41 @@ namespace AliceO2 {
1617

1718
using namespace AliceO2::Base::Constants;
1819
using namespace AliceO2::Base::Utils;
20+
using namespace std;
1921

2022
// aliases for track elements
21-
enum {kX,kAlpha,
22-
kY,kZ,kSnp,kTgl,kQ2Pt,
23-
kSigY2,
23+
enum {kY,kZ,kSnp,kTgl,kQ2Pt};
24+
enum {kSigY2,
2425
kSigZY,kSigZ2,
2526
kSigSnpY,kSigSnpZ,kSigSnp2,
2627
kSigTglY,kSigTglZ,kSigTglSnp,kSigTgl2,
2728
kSigQ2PtY,kSigQ2PtZ,kSigQ2PtSnp,kSigQ2PtTgl,kSigQ2Pt2};
28-
enum {kNParams=5,kCovMatSize=15,kTrackPSize=kNParams+2,kTrackPCSize=kTrackPSize+kCovMatSize
29-
,kLabCovMatSize=21};
3029

31-
constexpr float
32-
kCY2max=100*100, // SigmaY<=100cm
33-
kCZ2max=100*100, // SigmaZ<=100cm
34-
kCSnp2max=1*1, // SigmaSin<=1
35-
kCTgl2max=1*1, // SigmaTan<=1
36-
kC1Pt2max=100*100; // Sigma1/Pt<=100 1/GeV
30+
constexpr int
31+
kNParams=5,
32+
kCovMatSize=15,
33+
kLabCovMatSize=21;
3734

35+
constexpr float
36+
kCY2max = 100*100, // SigmaY<=100cm
37+
kCZ2max = 100*100, // SigmaZ<=100cm
38+
kCSnp2max = 1*1, // SigmaSin<=1
39+
kCTgl2max = 1*1, // SigmaTan<=1
40+
kC1Pt2max = 100 * 100, // Sigma1/Pt<=100 1/GeV
41+
kCalcdEdxAuto = -999.f ; // value indicating request for dedx calculation
42+
43+
// helper function
44+
float BetheBlochSolid(float bg, float rho=2.33f,float kp1=0.20f,float kp2=3.00f,
45+
float meanI=173e-9f,float meanZA=0.49848f);
46+
void g3helx3(float qfield, float step,float vect[7]);
3847

39-
class TrackPar { // track parameterization, kinematics only
40-
public:
41-
TrackPar() {memset(mP,0,kTrackPSize*sizeof(float));}
42-
TrackPar(float x,float alpha, const float par[kNParams]);
43-
TrackPar(const float xyz[3],const float pxpypz[3],int sign, bool sectorAlpha=true);
4448

45-
float& operator[](int i) { return mP[i]; }
46-
float operator[](int i) const { return mP[i]; }
47-
operator float*() const { return (float*)mP; }
49+
class TrackParBase { // track parameterization, kinematics only. This base class cannot be instantiated
50+
public:
4851

49-
const float* GetParam() const { return &mP[kY]; }
50-
float GetX() const { return mP[kX]; }
51-
float GetAlpha() const { return mP[kAlpha]; }
52+
const float* GetParam() const { return mP; }
53+
float GetX() const { return mX; }
54+
float GetAlpha() const { return mAlpha; }
5255
float GetY() const { return mP[kY]; }
5356
float GetZ() const { return mP[kZ]; }
5457
float GetSnp() const { return mP[kSnp]; }
@@ -73,73 +76,50 @@ namespace AliceO2 {
7376
bool PropagateParamTo(float xk, const float b[3]);
7477
void InvertParam();
7578

76-
void Print();
77-
78-
static void g3helx3(float qfield, float step,float vect[7]);
79+
void PrintParam() const;
7980

8081
protected:
81-
float mP[kTrackPSize]; // x,alpha + 5 parameters
82+
// to keep this class non-virtual but derivable the c-tors and d-tor are protected
83+
TrackParBase() : mX{0.},mAlpha{0.},mP{0.f} {}
84+
TrackParBase(float x,float alpha, const float par[kNParams]);
85+
TrackParBase(const float xyz[3],const float pxpypz[3],int sign, bool sectorAlpha=true);
86+
TrackParBase(const TrackParBase&) = default;
87+
TrackParBase(TrackParBase&&) = default;
88+
TrackParBase& operator=(const TrackParBase& src) = default;
89+
~TrackParBase() = default;
90+
//
91+
float mX; /// X of track evaluation
92+
float mAlpha; /// track frame angle
93+
float mP[kNParams]; /// 5 parameters: Y,Z,sin(phi),tg(lambda),q/pT
8294
};
8395

84-
class TrackParCov { // track+error parameterization
96+
// rootcint does not swallow final keyword here
97+
class TrackParCov final : public TrackParBase { // track+error parameterization
8598
public:
86-
TrackParCov() { memset(mPC,0,kTrackPCSize*sizeof(float)); }
99+
TrackParCov() : TrackParBase{}, mC{0.f} { }
87100
TrackParCov(float x,float alpha, const float par[kNParams], const float cov[kCovMatSize]);
88-
TrackParCov(const float xyz[3],const float pxpypz[3],const float[kLabCovMatSize],
89-
int sign, bool sectorAlpha=true);
90-
91-
operator TrackPar*() { return reinterpret_cast<TrackPar*>(this); }
92-
operator TrackPar() { return *reinterpret_cast<TrackPar*>(this); }
93-
operator TrackPar&() { return *reinterpret_cast<TrackPar*>(this); }
94-
95-
float& operator[](int i) { return mPC[i]; }
96-
float operator[](int i) const { return mPC[i]; }
97-
operator float*() const { return (float*)mPC; }
98-
const float* GetParam() const { return &mPC[kY]; }
99-
const float* GetCov() const { return &mPC[kSigY2]; }
100-
101-
float GetX() const { return mPC[kX]; }
102-
float GetAlpha() const { return mPC[kAlpha]; }
103-
float GetY() const { return mPC[kY]; }
104-
float GetZ() const { return mPC[kZ]; }
105-
float GetSnp() const { return mPC[kSnp]; }
106-
float GetTgl() const { return mPC[kTgl]; }
107-
float GetQ2Pt() const { return mPC[kQ2Pt]; }
108-
109-
110-
float GetSigmaY2() const { return mPC[kSigY2]; }
111-
float GetSigmaZY() const { return mPC[kSigZY]; }
112-
float GetSigmaZ2() const { return mPC[kSigZ2]; }
113-
float GetSigmaSnpY() const { return mPC[kSigSnpY]; }
114-
float GetSigmaSnpZ() const { return mPC[kSigSnpZ]; }
115-
float GetSigmaSnp2() const { return mPC[kSigSnp2]; }
116-
float GetSigmaTglY() const { return mPC[kSigTglY]; }
117-
float GetSigmaTglZ() const { return mPC[kSigTglZ]; }
118-
float GetSigmaTglSnp() const { return mPC[kSigTglSnp]; }
119-
float GetSigmaTgl2() const { return mPC[kSigTgl2]; }
120-
float GetSigma1PtY() const { return mPC[kSigQ2PtY]; }
121-
float GetSigma1PtZ() const { return mPC[kSigQ2PtZ]; }
122-
float GetSigma1PtSnp() const { return mPC[kSigQ2PtSnp]; }
123-
float GetSigma1PtTgl() const { return mPC[kSigQ2PtTgl]; }
124-
float GetSigma1Pt2() const { return mPC[kSigQ2Pt2]; }
125-
126-
// derived getters
127-
float GetCurvature(float b) const { return mPC[kQ2Pt]*b*kB2C;}
128-
float GetSign() const { return mPC[kQ2Pt]>0 ? 1.f:-1.f;}
129-
float GetP() const { return Param()->GetP(); }
130-
float GetPt() const { return Param()->GetPt(); }
131-
float GetPhi() const { return Param()->GetPhi(); }
132-
float GetPhiPos() const { return Param()->GetPhiPos(); }
133-
void GetXYZ(float xyz[3]) const { Param()->GetXYZ(xyz); }
134-
bool GetPxPyPz(float pxyz[3]) const { return Param()->GetPxPyPz(pxyz); }
135-
bool GetPosDir(float posdirp[9]) const { return Param()->GetPosDir(posdirp); }
136-
137-
// parameters manipulation
138-
bool RotateParam(float alpha) { return Param()->RotateParam(alpha); }
139-
bool PropagateParamTo(float xk, float b) { return Param()->PropagateParamTo(xk,b); }
140-
bool PropagateParamTo(float xk, const float b[3]) {return Param()->PropagateParamTo(xk,b); }
141-
void InvertParam() { Param()->InvertParam(); }
142-
101+
TrackParCov(const float xyz[3],const float pxpypz[3],const float[kLabCovMatSize], int sign, bool sectorAlpha=true);
102+
103+
const float* GetCov() const { return mC; }
104+
float GetSigmaY2() const { return mC[kSigY2]; }
105+
float GetSigmaZY() const { return mC[kSigZY]; }
106+
float GetSigmaZ2() const { return mC[kSigZ2]; }
107+
float GetSigmaSnpY() const { return mC[kSigSnpY]; }
108+
float GetSigmaSnpZ() const { return mC[kSigSnpZ]; }
109+
float GetSigmaSnp2() const { return mC[kSigSnp2]; }
110+
float GetSigmaTglY() const { return mC[kSigTglY]; }
111+
float GetSigmaTglZ() const { return mC[kSigTglZ]; }
112+
float GetSigmaTglSnp() const { return mC[kSigTglSnp]; }
113+
float GetSigmaTgl2() const { return mC[kSigTgl2]; }
114+
float GetSigma1PtY() const { return mC[kSigQ2PtY]; }
115+
float GetSigma1PtZ() const { return mC[kSigQ2PtZ]; }
116+
float GetSigma1PtSnp() const { return mC[kSigQ2PtSnp]; }
117+
float GetSigma1PtTgl() const { return mC[kSigQ2PtTgl]; }
118+
float GetSigma1Pt2() const { return mC[kSigQ2Pt2]; }
119+
120+
void Print() const;
121+
122+
// parameters + covmat manipulation
143123
bool Rotate(float alpha);
144124
bool PropagateTo(float xk, float b);
145125
bool PropagateTo(float xk, const float b[3]);
@@ -152,69 +132,75 @@ namespace AliceO2 {
152132

153133
void ResetCovariance(float s2=0);
154134
void CheckCovariance();
155-
void Print();
156-
157-
protected:
158-
// internal cast to TrackPar
159-
const TrackPar* Param() const { return reinterpret_cast<const TrackPar*>(this); }
160-
TrackPar* Param() { return reinterpret_cast<TrackPar*>(this); }
161-
bool TrackPar2Momentum(float p[3], float alpha);
162135

163136
protected:
164-
float mPC[kTrackPCSize]; // x, alpha + 5 parameters + 15 errors
137+
float mC[kCovMatSize]; // x, alpha + 5 parameters + 15 errors
165138

166-
static const float kCalcdEdxAuto; // value indicating request for dedx calculation
167139
};
168140

141+
class TrackPar final : public TrackParBase { // track parameterization only
142+
public:
143+
TrackPar() {}
144+
TrackPar(float x,float alpha, const float par[kNParams]) : TrackParBase{x,alpha,par} {}
145+
TrackPar(const float xyz[3], const float pxpypz[3],int sign, bool sectorAlpha=true);
146+
TrackPar(const TrackParCov& src) : TrackParBase{static_cast<const TrackParBase&>(src)} {}
147+
//
148+
void Print() const {PrintParam();}
149+
};
169150

170151
//____________________________________________________________
171-
inline TrackPar::TrackPar(float x, float alpha, const float par[kNParams]) {
152+
inline TrackParBase::TrackParBase(float x, float alpha, const float par[kNParams]) : mX{x}, mAlpha{alpha} {
172153
// explicit constructor
173-
mP[kX] = x;
174-
mP[kAlpha] = alpha;
175-
memcpy(&mP[kY],par,kNParams*sizeof(float));
154+
std::copy(par, par + kNParams, mP);
176155
}
177156

178157
//_______________________________________________________
179-
inline void TrackPar::GetXYZ(float xyz[3]) const {
158+
inline void TrackParBase::GetXYZ(float xyz[3]) const {
180159
// track coordinates in lab frame
181-
xyz[0] = GetX();
160+
xyz[0] = GetX();
182161
xyz[1] = GetY();
183162
xyz[2] = GetZ();
184163
RotateZ(xyz,GetAlpha());
185164
}
186165

187166
//_______________________________________________________
188-
inline float TrackPar::GetPhiPos() const {
167+
inline float TrackParBase::GetPhiPos() const {
189168
// angle of track position
190169
float xy[2]={GetX(),GetY()};
191170
return atan2(xy[1],xy[0]);
192171
}
193172

194173
//____________________________________________________________
195-
inline float TrackPar::GetP() const {
174+
inline float TrackParBase::GetP() const {
196175
// return the track momentum
197-
float ptI = fabs(mP[kQ2Pt]);
198-
return (ptI>kAlmost0) ? sqrtf(1.f+ mP[kTgl]*mP[kTgl])/ptI : kVeryBig;
176+
float ptI = fabs(GetQ2Pt());
177+
return (ptI>kAlmost0) ? sqrtf(1.f+ GetTgl()*GetTgl())/ptI : kVeryBig;
199178
}
200179

201180
//____________________________________________________________
202-
inline float TrackPar::GetPt() const {
181+
inline float TrackParBase::GetPt() const {
203182
// return the track transverse momentum
204-
float ptI = fabs(mP[kQ2Pt]);
183+
float ptI = fabs(GetQ2Pt());
205184
return (ptI>kAlmost0) ? 1.f/ptI : kVeryBig;
206185
}
207186

208187
//============================================================
209188

210189
//____________________________________________________________
211-
inline TrackParCov::TrackParCov(float x, float alpha, const float par[kNParams], const float cov[kCovMatSize]) {
190+
inline TrackParCov::TrackParCov(float x, float alpha, const float par[kNParams], const float cov[kCovMatSize])
191+
: TrackParBase{x,alpha,par} {
192+
// explicit constructor
193+
std::copy(cov, cov + kCovMatSize, mC);
194+
}
195+
196+
//============================================================
197+
198+
//____________________________________________________________
199+
inline TrackPar::TrackPar(const float xyz[3], const float pxpypz[3],int sign, bool sectorAlpha)
200+
: TrackParBase{xyz,pxpypz,sign,sectorAlpha} {
212201
// explicit constructor
213-
mPC[kX] = x;
214-
mPC[kAlpha] = alpha;
215-
memcpy(&mPC[kY],par,kNParams*sizeof(float));
216-
memcpy(&mPC[kSigY2],cov,kCovMatSize*sizeof(float));
217202
}
203+
218204
}
219205
}
220206
}

‎Detectors/Base/include/DetectorsBase/Utils.h‎

Lines changed: 9 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -10,27 +10,27 @@
1010

1111
namespace AliceO2 {
1212
namespace Base {
13-
namespace Utils {
13+
namespace Utils {
1414

1515
using namespace AliceO2::Base::Constants;
1616

17-
inline void BringTo02Pi(float &phi) {
17+
inline void BringTo02Pi(float &phi) {
1818
// ensure angle in [0:2pi] for the input in [-pi:pi] or [0:pi]
19-
if (phi < 0) phi += k2PI;
19+
if (phi < 0) phi += k2PI;
2020
}
2121

22-
inline void BringTo02PiGen(float &phi) {
22+
inline void BringTo02PiGen(float &phi) {
2323
// ensure angle in [0:2pi] for the any input angle
2424
while(phi<0) {phi += k2PI;}
25-
while(phi>k2PI) {phi -= k2PI;}
25+
while(phi>k2PI) {phi -= k2PI;}
2626
}
2727

28-
inline void BringToPMPi(float &phi) {
28+
inline void BringToPMPi(float &phi) {
2929
// ensure angle in [-pi:pi] for the input in [-pi:pi] or [0:pi]
30-
if (phi > kPI) phi -= k2PI;
30+
if (phi > kPI) phi -= k2PI;
3131
}
3232

33-
inline void BringToPMPiGen(float &phi) {
33+
inline void BringToPMPiGen(float &phi) {
3434
// ensure angle in [-pi:pi] for any input angle
3535
while(phi<-kPI) {phi += k2PI;}
3636
while(phi> kPI) {phi -= k2PI;}
@@ -46,7 +46,7 @@ namespace AliceO2 {
4646
// transforms vector in tracking frame alpha to global frame
4747
float sn,cs, x=xy[0];
4848
sincosf(alpha,sn,cs);
49-
xy[0]=x*cs - xy[1]*sn;
49+
xy[0]=x*cs - xy[1]*sn;
5050
xy[1]=x*sn + xy[1]*cs;
5151
}
5252

@@ -67,41 +67,6 @@ namespace AliceO2 {
6767
return Sector2Angle(Angle2Sector(phi));
6868
}
6969

70-
inline float BetheBlochSolid(float bg, float rho=2.33f,float kp1=0.20f,float kp2=3.00f,
71-
float meanI=173e-9f,float meanZA=0.49848f) {
72-
//
73-
// This is the parameterization of the Bethe-Bloch formula inspired by Geant.
74-
//
75-
// bg - beta*gamma
76-
// rho - density [g/cm^3]
77-
// kp1 - density effect first junction point
78-
// kp2 - density effect second junction point
79-
// meanI - mean excitation energy [GeV]
80-
// meanZA - mean Z/A
81-
//
82-
// The default values for the kp* parameters are for silicon.
83-
// The returned value is in [GeV/(g/cm^2)].
84-
//
85-
constexpr float mK = 0.307075e-3f; // [GeV*cm^2/g]
86-
constexpr float me = 0.511e-3f; // [GeV/c^2]
87-
kp1 *= 2.303f;
88-
kp2 *= 2.303f;
89-
float bg2 = bg*bg;
90-
float maxT= 2.f*me*bg2; // neglecting the electron mass
91-
92-
//*** Density effect
93-
float d2=0.;
94-
const float x = log(bg);
95-
const float lhwI = log(28.816*1e-9*sqrtf(rho*meanZA)/meanI);
96-
if (x > kp2) d2 = lhwI + x - 0.5;
97-
else if (x > kp1) {
98-
double r=(kp2-x)/(kp2-kp1);
99-
d2 = lhwI + x - 0.5 + (0.5 - lhwI - kp1)*r*r*r;
100-
}
101-
return mK*meanZA*(1+bg2)/bg2*(0.5*log(2*me*bg2*maxT/(meanI*meanI)) - bg2/(1+bg2) - d2);
102-
}
103-
104-
10570
}
10671
}
10772
}

0 commit comments

Comments
 (0)