Skip to content

Commit 575e725

Browse files
shahor02sawenzel
shahor02
authored and
sawenzel
committed
Ported fast field param from Shuto for full Barrel, |Z|<550
1 parent 734e5c8 commit 575e725

4 files changed

Lines changed: 7982 additions & 807 deletions

File tree

‎Common/Field/include/Field/MagFieldFast.h‎

Lines changed: 42 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -17,78 +17,75 @@
1717
#include <Rtypes.h>
1818
#include <string>
1919

20-
namespace o2 { namespace field {
21-
20+
namespace o2
21+
{
22+
namespace field
23+
{
2224
// Fast polynomial parametrization of Alice magnetic field, to be used for reconstruction.
23-
// Solenoid part fitted by Shuto Yamasaki from AliMagWrapCheb in the |Z|<260Interface and R<500 cm
25+
// Solenoid part fitted by Shuto Yamasaki from AliMagWrapCheb in the |Z|<260Interface and R<500 cm
2426
// Dipole part: to do
25-
class MagFieldFast
27+
class MagFieldFast
2628
{
27-
2829
public:
29-
enum {kNSolRRanges=5, kNSolZRanges=2, kNQuadrants=4, kNPolCoefs=20};
30-
enum EDim {kX,kY,kZ,kNDim};
31-
struct SolParam { float parBxyz[kNDim][kNPolCoefs];};
32-
33-
MagFieldFast(const std::string inpFName="");
34-
MagFieldFast(float factor, int nomField = 5, const std::string inpFmt="$(O2_ROOT)/share/Common/maps/sol%dk.txt");
30+
enum { kNSolRRanges = 5, kNSolZRanges = 22, kNQuadrants = 4, kNPolCoefs = 20 };
31+
enum EDim { kX, kY, kZ, kNDim };
32+
struct SolParam {
33+
float parBxyz[kNDim][kNPolCoefs];
34+
};
35+
36+
MagFieldFast(const std::string inpFName = "");
37+
MagFieldFast(float factor, int nomField = 5, const std::string inpFmt = "$(O2_ROOT)/share/Common/maps/sol%dk.txt");
3538
MagFieldFast(const MagFieldFast& src) = default;
3639
~MagFieldFast() = default;
37-
40+
3841
bool LoadData(const std::string inpFName);
3942

4043
bool Field(const double xyz[3], double bxyz[3]) const;
41-
bool Field(const float xyz[3], float bxyz[3]) const;
44+
bool Field(const float xyz[3], float bxyz[3]) const;
4245
bool GetBcomp(EDim comp, const double xyz[3], double& b) const;
43-
bool GetBcomp(EDim comp, const float xyz[3], float& b) const;
44-
45-
bool GetBx(const double xyz[3], double& bx) const {return GetBcomp(kX,xyz,bx);}
46-
bool GetBx(const float xyz[3], float& bx) const {return GetBcomp(kX,xyz,bx);}
47-
bool GetBy(const double xyz[3], double& by) const {return GetBcomp(kY,xyz,by);}
48-
bool GetBy(const float xyz[3], float& by) const {return GetBcomp(kY,xyz,by);}
49-
bool GetBz(const double xyz[3], double& bz) const {return GetBcomp(kZ,xyz,bz);}
50-
bool GetBz(const float xyz[3], float& bz) const {return GetBcomp(kZ,xyz,bz);}
51-
52-
void setFactorSol(float v=1.f) {mFactorSol = v;}
53-
float getFactorSol() const {return mFactorSol;}
54-
46+
bool GetBcomp(EDim comp, const float xyz[3], float& b) const;
47+
48+
bool GetBx(const double xyz[3], double& bx) const { return GetBcomp(kX, xyz, bx); }
49+
bool GetBx(const float xyz[3], float& bx) const { return GetBcomp(kX, xyz, bx); }
50+
bool GetBy(const double xyz[3], double& by) const { return GetBcomp(kY, xyz, by); }
51+
bool GetBy(const float xyz[3], float& by) const { return GetBcomp(kY, xyz, by); }
52+
bool GetBz(const double xyz[3], double& bz) const { return GetBcomp(kZ, xyz, bz); }
53+
bool GetBz(const float xyz[3], float& bz) const { return GetBcomp(kZ, xyz, bz); }
54+
void setFactorSol(float v = 1.f) { mFactorSol = v; }
55+
float getFactorSol() const { return mFactorSol; }
5556
protected:
57+
bool GetSegment(const float xyz[3], int& zSeg, int& rSeg, int& quadrant) const;
58+
static const float kSolR2Max[kNSolRRanges]; // Rmax2 of each range
59+
static const float kSolZMax; // max |Z| for solenoid parametrization
5660

57-
bool GetSegment(const float xyz[3], int& zSeg,int &rSeg, int &quadrant) const;
58-
static const float kSolR2Max[kNSolRRanges]; // Rmax2 of each range
59-
static const float kSolZMax; // max |Z| for solenoid parametrization
60-
61-
62-
int GetQuadrant(float x,float y) const
61+
int GetQuadrant(float x, float y) const
6362
{
6463
/// get point quadrant
65-
return y>0 ? (x>0 ? 0:1) : (x>0 ? 3:2);
64+
return y > 0 ? (x > 0 ? 0 : 1) : (x > 0 ? 3 : 2);
6665
}
67-
68-
float CalcPol(const float* cf, float x,float y, float z) const;
66+
67+
float CalcPol(const float* cf, float x, float y, float z) const;
6968

7069
float mFactorSol; // scaling factor
7170
SolParam mSolPar[kNSolRRanges][kNSolZRanges][kNQuadrants];
72-
73-
ClassDef(MagFieldFast,1)
71+
72+
ClassDef(MagFieldFast, 1)
7473
};
7574

76-
inline float MagFieldFast::CalcPol(const float* cf, float x,float y, float z) const
75+
inline float MagFieldFast::CalcPol(const float* cf, float x, float y, float z) const
7776
{
78-
7977
/** calculate polynomial
8078
* cf[0] + cf[1]*x + cf[2]*y + cf[3]*z + cf[4]*xx + cf[5]*xy + cf[6]*xz + cf[7]*yy + cf[8]*yz + cf[9]*zz +
81-
* cf[10]*xxx + cf[11]*xxy + cf[12]*xxz + cf[13]*xyy + cf[14]*xyz + cf[15]*xzz + cf[16]*yyy + cf[17]*yyz + cf[18]*yzz + cf[19]*zzz
79+
* cf[10]*xxx + cf[11]*xxy + cf[12]*xxz + cf[13]*xyy + cf[14]*xyz + cf[15]*xzz + cf[16]*yyy + cf[17]*yyz +
80+
*cf[18]*yzz + cf[19]*zzz
8281
**/
8382

84-
float val = cf[0] +
85-
x*(cf[1] + x*(cf[4] + x*cf[10] + y*cf[11] + z*cf[12]) + y*(cf[5]+z*cf[14]) ) +
86-
y*(cf[2] + y*(cf[7] + x*cf[13] + y*cf[16] + z*cf[17]) + z*(cf[8]) ) +
87-
z*(cf[3] + z*(cf[9] + x*cf[15] + y*cf[18] + z*cf[19]) + x*(cf[6]) );
88-
83+
float val = cf[0] + x * (cf[1] + x * (cf[4] + x * cf[10] + y * cf[11] + z * cf[12]) + y * (cf[5] + z * cf[14])) +
84+
y * (cf[2] + y * (cf[7] + x * cf[13] + y * cf[16] + z * cf[17]) + z * (cf[8])) +
85+
z * (cf[3] + z * (cf[9] + x * cf[15] + y * cf[18] + z * cf[19]) + x * (cf[6]));
86+
8987
return val;
9088
}
91-
9289
}
9390
}
9491

‎Common/Field/src/MagFieldFast.cxx‎

Lines changed: 75 additions & 68 deletions
Original file line numberDiff line numberDiff line change
@@ -23,17 +23,15 @@
2323
using namespace o2::field;
2424
using namespace std;
2525

26-
ClassImp(o2::field::MagFieldFast)
26+
ClassImp(o2::field::MagFieldFast);
2727

28+
const float MagFieldFast::kSolR2Max[MagFieldFast::kNSolRRanges] = { 80.f * 80.f, 250.f * 250.f, 400.f * 400.f,
29+
423.f * 423.f, 500.f * 500.f };
2830

29-
const float MagFieldFast::kSolR2Max[MagFieldFast::kNSolRRanges] =
30-
{80.f*80.f,250.f*250.f,400.f*400.f,423.f*423.f, 500.f*500.f};
31-
32-
const float MagFieldFast::kSolZMax = 260.0f;
31+
const float MagFieldFast::kSolZMax = 550.0f;
3332

3433
//_______________________________________________________________________
35-
MagFieldFast::MagFieldFast(const string inpFName) :
36-
mFactorSol(1.f)
34+
MagFieldFast::MagFieldFast(const string inpFName) : mFactorSol(1.f)
3735
{
3836
// c-tor
3937
if (!inpFName.empty() && !LoadData(inpFName)) {
@@ -42,15 +40,14 @@ mFactorSol(1.f)
4240
}
4341

4442
//_______________________________________________________________________
45-
MagFieldFast::MagFieldFast(float factor, int nomField, const string inpFmt) :
46-
mFactorSol(factor)
43+
MagFieldFast::MagFieldFast(float factor, int nomField, const string inpFmt) : mFactorSol(factor)
4744
{
4845
// c-tor
49-
if (nomField!=2 && nomField!=5) {
46+
if (nomField != 2 && nomField != 5) {
5047
LOG(FATAL) << "No parametrization for nominal field of " << nomField << " kG" << FairLogger::endl;
51-
}
48+
}
5249
TString pth;
53-
pth.Form(inpFmt.data(),nomField);
50+
pth.Form(inpFmt.data(), nomField);
5451
if (!LoadData(pth.Data())) {
5552
LOG(FATAL) << "Failed to initialize from " << pth.Data() << FairLogger::endl;
5653
}
@@ -61,62 +58,64 @@ bool MagFieldFast::LoadData(const string inpFName)
6158
{
6259
// load field from text file
6360

64-
std::ifstream in(gSystem->ExpandPathName(inpFName.data()),std::ifstream::in);
61+
std::ifstream in(gSystem->ExpandPathName(inpFName.data()), std::ifstream::in);
6562
if (in.fail()) {
6663
LOG(FATAL) << "Failed to open file " << inpFName << FairLogger::endl;
6764
return false;
6865
}
6966
std::string line;
70-
int valI,component = -1, nParams = 0, header[4] = {-1,-1,-1,-1}; // iR, iZ, iQuadrant, nVal
67+
int valI, component = -1, nParams = 0, header[4] = { -1, -1, -1, -1 }; // iR, iZ, iQuadrant, nVal
7168
SolParam* curParam = nullptr;
72-
73-
while (std::getline(in, line)) {
7469

75-
if (line.empty() || line[0]=='#') continue; // empy or comment
70+
while (std::getline(in, line)) {
71+
if (line.empty() || line[0] == '#')
72+
continue; // empy or comment
7673
std::stringstream ss(line);
77-
int cnt=0;
78-
79-
if (component<0) {
80-
while (cnt<4 && (ss>>header[cnt++]));
81-
if (cnt!=4) {
82-
LOG(FATAL) << "Wrong header " << line << FairLogger::endl;
83-
return false;
74+
int cnt = 0;
75+
76+
if (component < 0) {
77+
while (cnt < 4 && (ss >> header[cnt++]))
78+
;
79+
if (cnt != 4) {
80+
LOG(FATAL) << "Wrong header " << line << FairLogger::endl;
81+
return false;
8482
}
8583
curParam = &mSolPar[header[0]][header[1]][header[2]];
86-
}
87-
else {
88-
while (cnt<header[3] && (ss>>curParam->parBxyz[component][cnt++]));
89-
if (cnt!=header[3]) {
90-
LOG(FATAL) << "Wrong data (npar="<< cnt <<") for param " << header[0] << " " << header[1] << " "
91-
<< header[2] << " " << header[3] << " " << line << FairLogger::endl;
92-
return false;
84+
} else {
85+
while (cnt < header[3] && (ss >> curParam->parBxyz[component][cnt++]))
86+
;
87+
if (cnt != header[3]) {
88+
LOG(FATAL) << "Wrong data (npar=" << cnt << ") for param " << header[0] << " " << header[1] << " " << header[2]
89+
<< " " << header[3] << " " << line << FairLogger::endl;
90+
return false;
9391
}
94-
}
92+
}
9593
component++;
96-
if (component>2) {
94+
if (component > 2) {
9795
component = -1; // next header expected
9896
nParams++;
9997
}
10098
}
10199
//
102100
LOG(INFO) << "Loaded " << nParams << " params from " << inpFName << FairLogger::endl;
103-
if (nParams != kNSolRRanges*kNSolZRanges*kNQuadrants) {
104-
LOG(FATAL) << "Was expecting " << kNSolRRanges*kNSolZRanges*kNQuadrants << " params" << FairLogger::endl;
101+
if (nParams != kNSolRRanges * kNSolZRanges * kNQuadrants) {
102+
LOG(FATAL) << "Was expecting " << kNSolRRanges * kNSolZRanges * kNQuadrants << " params" << FairLogger::endl;
105103
}
106104
return true;
107105
}
108-
106+
109107
//_______________________________________________________________________
110108
bool MagFieldFast::Field(const double xyz[3], double bxyz[3]) const
111109
{
112110
// get field
113-
const float fxyz[3]={float(xyz[0]),float(xyz[1]),float(xyz[2])};
114-
int zSeg,rSeg,quadrant;
115-
if (!GetSegment(fxyz,zSeg,rSeg,quadrant)) return false;
116-
const SolParam *par = &mSolPar[rSeg][zSeg][quadrant];
117-
bxyz[kX] = CalcPol(par->parBxyz[kX],fxyz[kX],fxyz[kY],fxyz[kZ])*mFactorSol;
118-
bxyz[kY] = CalcPol(par->parBxyz[kY],fxyz[kX],fxyz[kY],fxyz[kZ])*mFactorSol;
119-
bxyz[kZ] = CalcPol(par->parBxyz[kZ],fxyz[kX],fxyz[kY],fxyz[kZ])*mFactorSol;
111+
const float fxyz[3] = { float(xyz[0]), float(xyz[1]), float(xyz[2]) };
112+
int zSeg, rSeg, quadrant;
113+
if (!GetSegment(fxyz, zSeg, rSeg, quadrant))
114+
return false;
115+
const SolParam* par = &mSolPar[rSeg][zSeg][quadrant];
116+
bxyz[kX] = CalcPol(par->parBxyz[kX], fxyz[kX], fxyz[kY], fxyz[kZ]) * mFactorSol;
117+
bxyz[kY] = CalcPol(par->parBxyz[kY], fxyz[kX], fxyz[kY], fxyz[kZ]) * mFactorSol;
118+
bxyz[kZ] = CalcPol(par->parBxyz[kZ], fxyz[kX], fxyz[kY], fxyz[kZ]) * mFactorSol;
120119
//
121120
return true;
122121
}
@@ -125,11 +124,12 @@ bool MagFieldFast::Field(const double xyz[3], double bxyz[3]) const
125124
bool MagFieldFast::GetBcomp(EDim comp, const double xyz[3], double& b) const
126125
{
127126
// get field
128-
const float fxyz[3]={float(xyz[0]),float(xyz[1]),float(xyz[2])};
129-
int zSeg,rSeg,quadrant;
130-
if (!GetSegment(fxyz,zSeg,rSeg,quadrant)) return false;
131-
const SolParam *par = &mSolPar[rSeg][zSeg][quadrant];
132-
b = CalcPol(par->parBxyz[comp],fxyz[kX],fxyz[kY],fxyz[kZ])*mFactorSol;
127+
const float fxyz[3] = { float(xyz[0]), float(xyz[1]), float(xyz[2]) };
128+
int zSeg, rSeg, quadrant;
129+
if (!GetSegment(fxyz, zSeg, rSeg, quadrant))
130+
return false;
131+
const SolParam* par = &mSolPar[rSeg][zSeg][quadrant];
132+
b = CalcPol(par->parBxyz[comp], fxyz[kX], fxyz[kY], fxyz[kZ]) * mFactorSol;
133133
//
134134
return true;
135135
}
@@ -138,10 +138,11 @@ bool MagFieldFast::GetBcomp(EDim comp, const double xyz[3], double& b) const
138138
bool MagFieldFast::GetBcomp(EDim comp, const float xyz[3], float& b) const
139139
{
140140
// get field
141-
int zSeg,rSeg,quadrant;
142-
if (!GetSegment(xyz,zSeg,rSeg,quadrant)) return false;
143-
const SolParam *par = &mSolPar[rSeg][zSeg][quadrant];
144-
b = CalcPol(par->parBxyz[comp],xyz[kX],xyz[kY],xyz[kZ])*mFactorSol;
141+
int zSeg, rSeg, quadrant;
142+
if (!GetSegment(xyz, zSeg, rSeg, quadrant))
143+
return false;
144+
const SolParam* par = &mSolPar[rSeg][zSeg][quadrant];
145+
b = CalcPol(par->parBxyz[comp], xyz[kX], xyz[kY], xyz[kZ]) * mFactorSol;
145146
//
146147
return true;
147148
}
@@ -150,33 +151,39 @@ bool MagFieldFast::GetBcomp(EDim comp, const float xyz[3], float& b) const
150151
bool MagFieldFast::Field(const float xyz[3], float bxyz[3]) const
151152
{
152153
// get field
153-
int zSeg,rSeg,quadrant;
154-
if (!GetSegment(xyz,zSeg,rSeg,quadrant)) return false;
155-
const SolParam *par = &mSolPar[rSeg][zSeg][quadrant];
156-
bxyz[kX] = CalcPol(par->parBxyz[kX],xyz[kX],xyz[kY],xyz[kZ])*mFactorSol;
157-
bxyz[kY] = CalcPol(par->parBxyz[kY],xyz[kX],xyz[kY],xyz[kZ])*mFactorSol;
158-
bxyz[kZ] = CalcPol(par->parBxyz[kZ],xyz[kX],xyz[kY],xyz[kZ])*mFactorSol;
154+
int zSeg, rSeg, quadrant;
155+
if (!GetSegment(xyz, zSeg, rSeg, quadrant))
156+
return false;
157+
const SolParam* par = &mSolPar[rSeg][zSeg][quadrant];
158+
bxyz[kX] = CalcPol(par->parBxyz[kX], xyz[kX], xyz[kY], xyz[kZ]) * mFactorSol;
159+
bxyz[kY] = CalcPol(par->parBxyz[kY], xyz[kX], xyz[kY], xyz[kZ]) * mFactorSol;
160+
bxyz[kZ] = CalcPol(par->parBxyz[kZ], xyz[kX], xyz[kY], xyz[kZ]) * mFactorSol;
159161
//
160162
return true;
161163
}
162164

163165
//_______________________________________________________________________
164-
bool MagFieldFast::GetSegment(const float xyz[3], int& zSeg,int &rSeg, int &quadrant) const
166+
bool MagFieldFast::GetSegment(const float xyz[3], int& zSeg, int& rSeg, int& quadrant) const
165167
{
166168
// get segment of point location
167169
const float &x = xyz[kX], &y = xyz[kY], &z = xyz[kZ];
170+
const float zGridSpaceInv = 1.f / (kSolZMax * 2 / kNSolZRanges);
168171
zSeg = -1;
169-
if (z<kSolZMax) {
170-
if (z>-kSolZMax) zSeg = z<0.f ? 0:1; // solenoid params
171-
else { // need to check dipole params
172+
if (z < kSolZMax) {
173+
if (z > -kSolZMax)
174+
zSeg = (z + kSolZMax) * zGridSpaceInv; // solenoid params
175+
else { // need to check dipole params
172176
return false;
173177
}
174-
}
175-
else return false;
178+
} else
179+
return false;
176180
// R segment
177-
float xx = x*x, yy = y*y, rr = xx+yy;
178-
for (rSeg=0;rSeg<kNSolRRanges;rSeg++) if (rr < kSolR2Max[rSeg]) break;
179-
if (rSeg==kNSolRRanges) return kFALSE;
180-
quadrant = GetQuadrant(x,y);
181+
float xx = x * x, yy = y * y, rr = xx + yy;
182+
for (rSeg = 0; rSeg < kNSolRRanges; rSeg++)
183+
if (rr < kSolR2Max[rSeg])
184+
break;
185+
if (rSeg == kNSolRRanges)
186+
return kFALSE;
187+
quadrant = GetQuadrant(x, y);
181188
return true;
182189
}

0 commit comments

Comments
 (0)