CMS 3D CMS Logo

List of all members | Classes | Public Member Functions | Private Attributes
MagneticFieldGrid Class Reference

#include <MagneticFieldGrid.h>

Classes

class  BVector
 
class  HeaderType3
 

Public Member Functions

int gridType ()
 returns value of GridType (and eventually prints the type + short description) More...
 
void interpolateAtPoint (double X1, double X2, double X3, float &Bx, float &By, float &Bz)
 interpolates the magnetic field at input coordinate point and returns field values More...
 
int lineNumber (int Index1, int Index2, int Index3)
 
void load (const std::string &name)
 load grid binary file More...
 
 MagneticFieldGrid ()
 
void putCoordGetInd (double X1, double X2, double X3, int &Index1, int &Index2, int &Index3)
 
void putIndGetCoord (int Index1, int Index2, int Index3, double &X1, double &X2, double &X3)
 
void putIndicesGetB (int Index1, int Index2, int Index3, float &Bx, float &By, float &Bz)
 
 ~MagneticFieldGrid ()
 

Private Attributes

double BasicDistance0 [3]
 
double BasicDistance1 [3][3]
 
double BasicDistance2 [3][3]
 
bool EasyCoordinate [3]
 
std::vector< BVectorFieldValues
 
int GridType
 
int NumberOfPoints [3]
 
double ReferencePoint [3]
 
double RParAsFunOfPhi [4]
 

Detailed Description

load magnetic field grid from binary file remark: units are either (cm,cm,cm) or (cm,rad,cm) and Tesla for the magnetic field

additional functions either translate indices <-> coordinates, transfer data, or activate the interpolation between grid points

Author
: Volke.nosp@m.r.Dr.nosp@m.ollin.nosp@m.ger@.nosp@m.cern..nosp@m.ch

Modifications:

Definition at line 32 of file MagneticFieldGrid.h.

Constructor & Destructor Documentation

MagneticFieldGrid::MagneticFieldGrid ( )
inline

Definition at line 35 of file MagneticFieldGrid.h.

References mps_fire::i.

35  {
36  GridType = 0;
37  for (int i=0;i<3; ++i) {NumberOfPoints[i] = 0;};
38  for (int i=0;i<3; ++i) {ReferencePoint[i] = 0.;};
39  for (int i=0;i<3; ++i) {BasicDistance0[i] = 0.;};
40  for (int i=0;i<3; ++i) {for (int j=0;j<3; ++j) {BasicDistance1[i][j] = 0.;};};
41  for (int i=0;i<3; ++i) {for (int j=0;j<3; ++j) {BasicDistance2[i][j] = 0.;};};
42  for (int i=0;i<4; ++i) {RParAsFunOfPhi[i] = 0.;};
43  for (int i=0;i<3; ++i) {EasyCoordinate[i] = false;};
44  }
double BasicDistance1[3][3]
double BasicDistance2[3][3]
MagneticFieldGrid::~MagneticFieldGrid ( )
inline

Definition at line 46 of file MagneticFieldGrid.h.

46 {}

Member Function Documentation

int MagneticFieldGrid::gridType ( )

returns value of GridType (and eventually prints the type + short description)

Definition at line 77 of file MagneticFieldGrid.cc.

References gather_cfg::cout, and RecoTauValidation_cfi::text.

Referenced by GlobalGridWrapper::valueInTesla().

77  {
78  int type = GridType;
79  bool text = false;
80  if (text){
81  if (type == 0) cout << " grid type = " << type << " --> not determined" << endl;
82  if (type == 1) cout << " grid type = " << type << " --> (x,y,z) cube" << endl;
83  if (type == 2) cout << " grid type = " << type << " --> (x,y,z) trapezoid" << endl;
84  if (type == 3) cout << " grid type = " << type << " --> (r,phi,z) cube" << endl;
85  if (type == 4) cout << " grid type = " << type << " --> (r,phi,z) trapezoid" << endl;
86  if (type == 5) cout << " grid type = " << type << " --> (r,phi,z) 1/sin(phi)" << endl;
87  }
88  return type;
89 }
type
Definition: HCALResponse.h:21
void MagneticFieldGrid::interpolateAtPoint ( double  X1,
double  X2,
double  X3,
float &  Bx,
float &  By,
float &  Bz 
)

interpolates the magnetic field at input coordinate point and returns field values

Definition at line 91 of file MagneticFieldGrid.cc.

References VectorFieldInterpolation::defineCellPoint000(), VectorFieldInterpolation::defineCellPoint001(), VectorFieldInterpolation::defineCellPoint010(), VectorFieldInterpolation::defineCellPoint011(), VectorFieldInterpolation::defineCellPoint100(), VectorFieldInterpolation::defineCellPoint101(), VectorFieldInterpolation::defineCellPoint110(), VectorFieldInterpolation::defineCellPoint111(), objects.autophobj::float, mps_fire::i, diffTreeTool::index, hpstanc_transforms::max, and VectorFieldInterpolation::putSCoordGetVField().

Referenced by GlobalGridWrapper::valueInTesla().

91  {
92  double dB[3] = {0.,0.,0.};
93  // define interpolation object
94  VectorFieldInterpolation MagInterpol;
95  // calculate indices for "CellPoint000"
96  int index[3];
97  putCoordGetInd(X1,X2,X3,index[0],index[1],index[2]);
98  int index0[3] = {0,0,0};
99  int index1[3] = {0,0,0};
100  for (int i=0; i<3; ++i){
101  if (NumberOfPoints[i] > 1){
102  index0[i] = max(0,index[i]);
103  if (index0[i] > NumberOfPoints[i]-2) index0[i] = NumberOfPoints[i]-2;
104  index1[i] = max(1,index[i]+1);
105  if (index1[i] > NumberOfPoints[i]-1) index1[i] = NumberOfPoints[i]-1;
106  }
107  }
108  double tmpX[3];
109  float tmpB[3];
110  // define the corners of interpolation volume
111  // FIXME: should not unpack the arrays to then repack them as first thing.
112  putIndicesGetB(index0[0],index0[1],index0[2],tmpB[0],tmpB[1],tmpB[2]);
113  putIndGetCoord(index0[0],index0[1],index0[2],tmpX[0],tmpX[1],tmpX[2]);
114  MagInterpol.defineCellPoint000(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
115  putIndicesGetB(index1[0],index0[1],index0[2],tmpB[0],tmpB[1],tmpB[2]);
116  putIndGetCoord(index1[0],index0[1],index0[2],tmpX[0],tmpX[1],tmpX[2]);
117  MagInterpol.defineCellPoint100(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
118  putIndicesGetB(index0[0],index1[1],index0[2],tmpB[0],tmpB[1],tmpB[2]);
119  putIndGetCoord(index0[0],index1[1],index0[2],tmpX[0],tmpX[1],tmpX[2]);
120  MagInterpol.defineCellPoint010(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
121  putIndicesGetB(index1[0],index1[1],index0[2],tmpB[0],tmpB[1],tmpB[2]);
122  putIndGetCoord(index1[0],index1[1],index0[2],tmpX[0],tmpX[1],tmpX[2]);
123  MagInterpol.defineCellPoint110(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
124  putIndicesGetB(index0[0],index0[1],index1[2],tmpB[0],tmpB[1],tmpB[2]);
125  putIndGetCoord(index0[0],index0[1],index1[2],tmpX[0],tmpX[1],tmpX[2]);
126  MagInterpol.defineCellPoint001(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
127  putIndicesGetB(index1[0],index0[1],index1[2],tmpB[0],tmpB[1],tmpB[2]);
128  putIndGetCoord(index1[0],index0[1],index1[2],tmpX[0],tmpX[1],tmpX[2]);
129  MagInterpol.defineCellPoint101(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
130  putIndicesGetB(index0[0],index1[1],index1[2],tmpB[0],tmpB[1],tmpB[2]);
131  putIndGetCoord(index0[0],index1[1],index1[2],tmpX[0],tmpX[1],tmpX[2]);
132  MagInterpol.defineCellPoint011(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
133  putIndicesGetB(index1[0],index1[1],index1[2],tmpB[0],tmpB[1],tmpB[2]);
134  putIndGetCoord(index1[0],index1[1],index1[2],tmpX[0],tmpX[1],tmpX[2]);
135  MagInterpol.defineCellPoint111(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
136  // interpolate
137  MagInterpol.putSCoordGetVField(X1,X2,X3,dB[0],dB[1],dB[2]);
138  Bx = float(dB[0]);
139  By = float(dB[1]);
140  Bz = float(dB[2]);
141  return;
142 }
void putCoordGetInd(double X1, double X2, double X3, int &Index1, int &Index2, int &Index3)
void defineCellPoint100(double X1, double X2, double X3, double F1, double F2, double F3)
void defineCellPoint110(double X1, double X2, double X3, double F1, double F2, double F3)
void putIndGetCoord(int Index1, int Index2, int Index3, double &X1, double &X2, double &X3)
void defineCellPoint101(double X1, double X2, double X3, double F1, double F2, double F3)
void putSCoordGetVField(double X1, double X2, double X3, double &F1, double &F2, double &F3)
receive the interpolated field (out) at any point in space (in)
void defineCellPoint010(double X1, double X2, double X3, double F1, double F2, double F3)
void defineCellPoint111(double X1, double X2, double X3, double F1, double F2, double F3)
void defineCellPoint000(double X1, double X2, double X3, double F1, double F2, double F3)
provide the interpolation algorithm with 8 points, where the field is known (in)
void putIndicesGetB(int Index1, int Index2, int Index3, float &Bx, float &By, float &Bz)
void defineCellPoint001(double X1, double X2, double X3, double F1, double F2, double F3)
void defineCellPoint011(double X1, double X2, double X3, double F1, double F2, double F3)
int MagneticFieldGrid::lineNumber ( int  Index1,
int  Index2,
int  Index3 
)

Definition at line 292 of file MagneticFieldGrid.cc.

292  {
293  return Index1*NumberOfPoints[1]*NumberOfPoints[2] + Index2*NumberOfPoints[2] + Index3;
294 }
void MagneticFieldGrid::load ( const std::string &  name)

load grid binary file

Definition at line 7 of file MagneticFieldGrid.cc.

References binary_ifstream::close(), gather_cfg::cout, and MagneticFieldGrid::BVector::putB3().

Referenced by GlobalGridWrapper::GlobalGridWrapper().

7  {
8  binary_ifstream inFile(name);
9  inFile >> GridType;
10  // reading the header
11  switch (GridType){
12  case 1:
13  inFile >> NumberOfPoints[0] >> NumberOfPoints[1] >> NumberOfPoints[2];
14  inFile >> ReferencePoint[0] >> ReferencePoint[1] >> ReferencePoint[2];
15  inFile >> BasicDistance0[0] >> BasicDistance0[1] >> BasicDistance0[2];
16  break;
17  case 2:
18  inFile >> NumberOfPoints[0] >> NumberOfPoints[1] >> NumberOfPoints[2];
19  inFile >> ReferencePoint[0] >> ReferencePoint[1] >> ReferencePoint[2];
20  inFile >> BasicDistance0[0] >> BasicDistance0[1] >> BasicDistance0[2];
21  inFile >> BasicDistance1[0][0] >> BasicDistance1[1][0] >> BasicDistance1[2][0];
22  inFile >> BasicDistance1[0][1] >> BasicDistance1[1][1] >> BasicDistance1[2][1];
23  inFile >> BasicDistance1[0][2] >> BasicDistance1[1][2] >> BasicDistance1[2][2];
24  inFile >> BasicDistance2[0][0] >> BasicDistance2[1][0] >> BasicDistance2[2][0];
25  inFile >> BasicDistance2[0][1] >> BasicDistance2[1][1] >> BasicDistance2[2][1];
26  inFile >> BasicDistance2[0][2] >> BasicDistance2[1][2] >> BasicDistance2[2][2];
27  inFile >> EasyCoordinate[0] >> EasyCoordinate[1] >> EasyCoordinate[2];
28  break;
29  case 3:
30  inFile >> NumberOfPoints[0] >> NumberOfPoints[1] >> NumberOfPoints[2];
31  inFile >> ReferencePoint[0] >> ReferencePoint[1] >> ReferencePoint[2];
32  inFile >> BasicDistance0[0] >> BasicDistance0[1] >> BasicDistance0[2];
33  break;
34  case 4:
35  inFile >> NumberOfPoints[0] >> NumberOfPoints[1] >> NumberOfPoints[2];
36  inFile >> ReferencePoint[0] >> ReferencePoint[1] >> ReferencePoint[2];
37  inFile >> BasicDistance0[0] >> BasicDistance0[1] >> BasicDistance0[2];
38  inFile >> BasicDistance1[0][0] >> BasicDistance1[1][0] >> BasicDistance1[2][0];
39  inFile >> BasicDistance1[0][1] >> BasicDistance1[1][1] >> BasicDistance1[2][1];
40  inFile >> BasicDistance1[0][2] >> BasicDistance1[1][2] >> BasicDistance1[2][2];
41  inFile >> BasicDistance2[0][0] >> BasicDistance2[1][0] >> BasicDistance2[2][0];
42  inFile >> BasicDistance2[0][1] >> BasicDistance2[1][1] >> BasicDistance2[2][1];
43  inFile >> BasicDistance2[0][2] >> BasicDistance2[1][2] >> BasicDistance2[2][2];
44  inFile >> EasyCoordinate[0] >> EasyCoordinate[1] >> EasyCoordinate[2];
45  break;
46  case 5:
47  inFile >> NumberOfPoints[0] >> NumberOfPoints[1] >> NumberOfPoints[2];
48  inFile >> ReferencePoint[0] >> ReferencePoint[1] >> ReferencePoint[2];
49  inFile >> BasicDistance0[0] >> BasicDistance0[1] >> BasicDistance0[2];
50  inFile >> RParAsFunOfPhi[0] >> RParAsFunOfPhi[1] >> RParAsFunOfPhi[2] >> RParAsFunOfPhi[3];
51  break;
52  default:
53  assert(0); //this is a bug
54  }
55  //reading the field
56  float Bx, By, Bz;
57  BVector FieldEntry;
58  int nLines = NumberOfPoints[0]*NumberOfPoints[1]*NumberOfPoints[2];
59  FieldValues.reserve(nLines);
60 
61  for (int iLine=0; iLine<nLines; ++iLine){
62  inFile >> Bx >> By >> Bz;
63  FieldEntry.putB3(Bx,By,Bz);
64  FieldValues.push_back(FieldEntry);
65  }
66  // check completeness and close file
67  string lastEntry;
68  inFile >> lastEntry;
69  inFile.close();
70  if (lastEntry != "complete"){
71  GridType = 0;
72  cout << "error during file reading: file is not complete" << endl;
73  }
74  return;
75 }
std::vector< BVector > FieldValues
double BasicDistance1[3][3]
double BasicDistance2[3][3]
void MagneticFieldGrid::putCoordGetInd ( double  X1,
double  X2,
double  X3,
int &  Index1,
int &  Index2,
int &  Index3 
)

Definition at line 148 of file MagneticFieldGrid.cc.

References mps_fire::i, diffTreeTool::index, createfilelist::int, PFRecoTauDiscriminationByIsolation_cfi::offset, and funct::sin().

148  {
149  double pnt[3] = {X1,X2,X3};
150  int index[3];
151  switch (GridType){
152  case 1:{
153  for (int i=0; i<3; ++i){
154  index[i] = int((pnt[i]-ReferencePoint[i])/BasicDistance0[i]);
155  }
156  break;}
157  case 2:{
158  // FIXME: Should use else!
159  for (int i=0; i<3; ++i){
160  if (EasyCoordinate[i]){
161  index[i] = int((pnt[i]-ReferencePoint[i])/BasicDistance0[i]);
162  } else index[i] = 0;//computed below
163  }
164  for (int i=0; i<3; ++i){
165  if (!EasyCoordinate[i]){
166  double stepSize = BasicDistance0[i];
167  double offset = 0.0;
168  for (int j=0; j<3; ++j){
169  stepSize += BasicDistance1[i][j]*index[j];
170  offset += BasicDistance2[i][j]*index[j];
171  }
172  index[i] = int((pnt[i]-(ReferencePoint[i] + offset))/stepSize);
173  }
174  }
175  break;}
176  case 3:{
177  for (int i=0; i<3; ++i){
178  index[i] = int((pnt[i]-ReferencePoint[i])/BasicDistance0[i]);
179  }
180  break;}
181  case 4:{
182  // FIXME: should use else!
183  for (int i=0; i<3; ++i){
184  if (EasyCoordinate[i]){
185  index[i] = int((pnt[i]-ReferencePoint[i])/BasicDistance0[i]);
186  } else index[i] = 0;//computed below
187  }
188  for (int i=0; i<3; ++i){
189  if (!EasyCoordinate[i]){
190  double stepSize = BasicDistance0[i];
191  double offset = 0.0;
192  for (int j=0; j<3; ++j){
193  stepSize += BasicDistance1[i][j]*index[j];
194  offset += BasicDistance2[i][j]*index[j];
195  }
196  index[i] = int((pnt[i]-(ReferencePoint[i] + offset))/stepSize);
197  }
198  }
199  break;}
200  case 5:{
201  double sinPhi = sin(pnt[1]);
202  double stepSize = RParAsFunOfPhi[0] + RParAsFunOfPhi[1]/sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3]/sinPhi;
203  stepSize = stepSize/(NumberOfPoints[0]-1);
204  double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3]/sinPhi;
205  index[0] = int((pnt[0]-startingPoint)/stepSize);
206  index[1] = int((pnt[1]-ReferencePoint[1])/BasicDistance0[1]);
207  index[2] = int((pnt[2]-ReferencePoint[2])/BasicDistance0[2]);
208  break;}
209  default:
210  assert(0); //shouldn't be here
211  }
212  Index1 = index[0];
213  Index2 = index[1];
214  Index3 = index[2];
215  return;
216 }
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
double BasicDistance1[3][3]
double BasicDistance2[3][3]
void MagneticFieldGrid::putIndGetCoord ( int  Index1,
int  Index2,
int  Index3,
double &  X1,
double &  X2,
double &  X3 
)

Definition at line 228 of file MagneticFieldGrid.cc.

References mps_fire::i, diffTreeTool::index, PFRecoTauDiscriminationByIsolation_cfi::offset, and funct::sin().

228  {
229  int index[3] = {Index1, Index2, Index3};
230  double pnt[3];
231  switch (GridType){
232  case 1:{
233  for (int i=0; i<3; ++i){
234  pnt[i] = ReferencePoint[i] + BasicDistance0[i]*index[i];
235  }
236  break;}
237  case 2:{
238  for (int i=0; i<3; ++i){
239  if (EasyCoordinate[i]){
240  pnt[i] = ReferencePoint[i] + BasicDistance0[i]*index[i];
241  }
242  else {
243  double stepSize = BasicDistance0[i];
244  double offset = 0.0;
245  for (int j=0; j<3; ++j){
246  stepSize += BasicDistance1[i][j]*index[j];
247  offset += BasicDistance2[i][j]*index[j];
248  }
249  pnt[i] = ReferencePoint[i] + offset + stepSize*index[i];
250  }
251  }
252  break;}
253  case 3:{
254  for (int i=0; i<3; ++i){
255  pnt[i] = ReferencePoint[i] + BasicDistance0[i]*index[i];
256  }
257  break;}
258  case 4:{
259  for (int i=0; i<3; ++i){
260  if (EasyCoordinate[i]){
261  pnt[i] = ReferencePoint[i] + BasicDistance0[i]*index[i];
262  }
263  else {
264  double stepSize = BasicDistance0[i];
265  double offset = 0.0;
266  for (int j=0; j<3; ++j){
267  stepSize += BasicDistance1[i][j]*index[j];
268  offset += BasicDistance2[i][j]*index[j];
269  }
270  pnt[i] = ReferencePoint[i] + offset + stepSize*index[i];
271  }
272  }
273  break;}
274  case 5:{
275  pnt[2] = ReferencePoint[2] + BasicDistance0[2]*index[2];
276  pnt[1] = ReferencePoint[1] + BasicDistance0[1]*index[1];
277  double sinPhi = sin(pnt[1]);
278  double stepSize = RParAsFunOfPhi[0] + RParAsFunOfPhi[1]/sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3]/sinPhi;
279  stepSize = stepSize/(NumberOfPoints[0]-1);
280  double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3]/sinPhi;
281  pnt[0] = startingPoint + stepSize*index[0];
282  break;}
283  default:
284  assert(0);//bug if make it here
285  }
286  X1 = pnt[0];
287  X2 = pnt[1];
288  X3 = pnt[2];
289  return;
290 }
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
double BasicDistance1[3][3]
double BasicDistance2[3][3]
void MagneticFieldGrid::putIndicesGetB ( int  Index1,
int  Index2,
int  Index3,
float &  Bx,
float &  By,
float &  Bz 
)

Definition at line 218 of file MagneticFieldGrid.cc.

References MagneticFieldGrid::BVector::bx(), MagneticFieldGrid::BVector::by(), and MagneticFieldGrid::BVector::bz().

218  {
219  BVector FieldEntry;
220  FieldEntry = FieldValues.operator[](lineNumber(Index1, Index2, Index3));
221  Bx = FieldEntry.bx();
222  By = FieldEntry.by();
223  Bz = FieldEntry.bz();
224  return;
225 }
std::vector< BVector > FieldValues
int lineNumber(int Index1, int Index2, int Index3)

Member Data Documentation

double MagneticFieldGrid::BasicDistance0[3]
private

Definition at line 84 of file MagneticFieldGrid.h.

double MagneticFieldGrid::BasicDistance1[3][3]
private

Definition at line 85 of file MagneticFieldGrid.h.

double MagneticFieldGrid::BasicDistance2[3][3]
private

Definition at line 86 of file MagneticFieldGrid.h.

bool MagneticFieldGrid::EasyCoordinate[3]
private

Definition at line 88 of file MagneticFieldGrid.h.

std::vector<BVector> MagneticFieldGrid::FieldValues
private

Definition at line 90 of file MagneticFieldGrid.h.

int MagneticFieldGrid::GridType
private

Definition at line 80 of file MagneticFieldGrid.h.

int MagneticFieldGrid::NumberOfPoints[3]
private

Definition at line 82 of file MagneticFieldGrid.h.

double MagneticFieldGrid::ReferencePoint[3]
private

Definition at line 83 of file MagneticFieldGrid.h.

double MagneticFieldGrid::RParAsFunOfPhi[4]
private

Definition at line 87 of file MagneticFieldGrid.h.