CMS 3D CMS Logo
MagneticFieldGrid.cc
Go to the documentation of this file.
1 // include header for MagneticFieldGrid (regular + extension for some trapezoids)
2 #include "MagneticFieldGrid.h"
3 #include <cassert>
4 
5 using namespace std;
6 
7 void MagneticFieldGrid::load(const string& name){
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 }
76 
77 int MagneticFieldGrid::gridType(){
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 }
90 
91 void MagneticFieldGrid::interpolateAtPoint(double X1, double X2, double X3, float &Bx, float &By, float &Bz){
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 }
143 
144 // FIXME: Signature should be:
145 //
146 // void MagneticFieldGrid::putCoordGetInd(double *pos, int *index)
147 //
148 void MagneticFieldGrid::putCoordGetInd(double X1, double X2, double X3, int &Index1, int &Index2, int &Index3){
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 }
217 
218 void MagneticFieldGrid::putIndicesGetB(int Index1, int Index2, int Index3, float &Bx, float &By, float &Bz){
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 }
226 
227 // FIXME: same as above.
228 void MagneticFieldGrid::putIndGetCoord(int Index1, int Index2, int Index3, double &X1, double &X2, double &X3){
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 }
291 
292 int MagneticFieldGrid::lineNumber(int Index1, int Index2, int Index3){
293  return Index1*NumberOfPoints[1]*NumberOfPoints[2] + Index2*NumberOfPoints[2] + Index3;
294 }
295 
296 void MagneticFieldGrid::BVector::putB3(float Bx, float By, float Bz){
297  B3[0] = Bx;
298  B3[1] = By;
299  B3[2] = Bz;
300  return;
301 }
302 
303 float MagneticFieldGrid::BVector::bx(){ return B3[0]; }
304 
305 float MagneticFieldGrid::BVector::by(){ return B3[1]; }
306 
307 float MagneticFieldGrid::BVector::bz(){ return B3[2]; }
MagneticFieldGrid::load
void load(const std::string &name)
load grid binary file
Definition: MagneticFieldGrid.cc:7
type
type
Definition: HCALResponse.h:21
mps_fire.i
i
Definition: mps_fire.py:338
MagneticFieldGrid::putCoordGetInd
void putCoordGetInd(double X1, double X2, double X3, int &Index1, int &Index2, int &Index3)
Definition: MagneticFieldGrid.cc:148
VectorFieldInterpolation::defineCellPoint100
void defineCellPoint100(double X1, double X2, double X3, double F1, double F2, double F3)
Definition: VectorFieldInterpolation.cc:16
MagneticFieldGrid::BVector::putB3
void putB3(float Bx, float By, float Bz)
Definition: MagneticFieldGrid.cc:296
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
MagneticFieldGrid::gridType
int gridType()
returns value of GridType (and eventually prints the type + short description)
Definition: MagneticFieldGrid.cc:77
binary_ifstream::close
void close()
Definition: binary_ifstream.cc:32
RecoTauValidation_cfi.text
text
Definition: RecoTauValidation_cfi.py:338
VectorFieldInterpolation::defineCellPoint110
void defineCellPoint110(double X1, double X2, double X3, double F1, double F2, double F3)
Definition: VectorFieldInterpolation.cc:38
MagneticFieldGrid::putIndGetCoord
void putIndGetCoord(int Index1, int Index2, int Index3, double &X1, double &X2, double &X3)
Definition: MagneticFieldGrid.cc:228
VectorFieldInterpolation::defineCellPoint101
void defineCellPoint101(double X1, double X2, double X3, double F1, double F2, double F3)
Definition: VectorFieldInterpolation.cc:60
VectorFieldInterpolation::putSCoordGetVField
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)
Definition: VectorFieldInterpolation.cc:93
VectorFieldInterpolation::defineCellPoint010
void defineCellPoint010(double X1, double X2, double X3, double F1, double F2, double F3)
Definition: VectorFieldInterpolation.cc:27
MagneticFieldGrid::lineNumber
int lineNumber(int Index1, int Index2, int Index3)
Definition: MagneticFieldGrid.cc:292
VectorFieldInterpolation::defineCellPoint111
void defineCellPoint111(double X1, double X2, double X3, double F1, double F2, double F3)
Definition: VectorFieldInterpolation.cc:82
VectorFieldInterpolation::defineCellPoint000
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)
Definition: VectorFieldInterpolation.cc:5
PFRecoTauDiscriminationByIsolation_cfi.offset
offset
Definition: PFRecoTauDiscriminationByIsolation_cfi.py:62
MagneticFieldGrid::putIndicesGetB
void putIndicesGetB(int Index1, int Index2, int Index3, float &Bx, float &By, float &Bz)
Definition: MagneticFieldGrid.cc:218
binary_ifstream
Definition: binary_ifstream.h:8
VectorFieldInterpolation::defineCellPoint001
void defineCellPoint001(double X1, double X2, double X3, double F1, double F2, double F3)
Definition: VectorFieldInterpolation.cc:49
gather_cfg.cout
cout
Definition: gather_cfg.py:144
VectorFieldInterpolation::defineCellPoint011
void defineCellPoint011(double X1, double X2, double X3, double F1, double F2, double F3)
Definition: VectorFieldInterpolation.cc:71
MagneticFieldGrid::interpolateAtPoint
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 ...
Definition: MagneticFieldGrid.cc:91
Generated for CMSSW Reference Manual by doxygen 1.8.11