CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MagneticFieldGrid.cc
Go to the documentation of this file.
1 // include header for MagneticFieldGrid (regular + extension for some trapezoids)
2 #include "MagneticField/Interpolation/src/MagneticFieldGrid.h"
3 
4 using namespace std;
5 
6 void MagneticFieldGrid::load(const string& name){
7  binary_ifstream inFile(name);
8  inFile >> GridType;
9  // reading the header
10  switch (GridType){
11  case 1:
12  inFile >> NumberOfPoints[0] >> NumberOfPoints[1] >> NumberOfPoints[2];
13  inFile >> ReferencePoint[0] >> ReferencePoint[1] >> ReferencePoint[2];
14  inFile >> BasicDistance0[0] >> BasicDistance0[1] >> BasicDistance0[2];
15  break;
16  case 2:
17  inFile >> NumberOfPoints[0] >> NumberOfPoints[1] >> NumberOfPoints[2];
18  inFile >> ReferencePoint[0] >> ReferencePoint[1] >> ReferencePoint[2];
19  inFile >> BasicDistance0[0] >> BasicDistance0[1] >> BasicDistance0[2];
20  inFile >> BasicDistance1[0][0] >> BasicDistance1[1][0] >> BasicDistance1[2][0];
21  inFile >> BasicDistance1[0][1] >> BasicDistance1[1][1] >> BasicDistance1[2][1];
22  inFile >> BasicDistance1[0][2] >> BasicDistance1[1][2] >> BasicDistance1[2][2];
23  inFile >> BasicDistance2[0][0] >> BasicDistance2[1][0] >> BasicDistance2[2][0];
24  inFile >> BasicDistance2[0][1] >> BasicDistance2[1][1] >> BasicDistance2[2][1];
25  inFile >> BasicDistance2[0][2] >> BasicDistance2[1][2] >> BasicDistance2[2][2];
26  inFile >> EasyCoordinate[0] >> EasyCoordinate[1] >> EasyCoordinate[2];
27  break;
28  case 3:
29  inFile >> NumberOfPoints[0] >> NumberOfPoints[1] >> NumberOfPoints[2];
30  inFile >> ReferencePoint[0] >> ReferencePoint[1] >> ReferencePoint[2];
31  inFile >> BasicDistance0[0] >> BasicDistance0[1] >> BasicDistance0[2];
32  break;
33  case 4:
34  inFile >> NumberOfPoints[0] >> NumberOfPoints[1] >> NumberOfPoints[2];
35  inFile >> ReferencePoint[0] >> ReferencePoint[1] >> ReferencePoint[2];
36  inFile >> BasicDistance0[0] >> BasicDistance0[1] >> BasicDistance0[2];
37  inFile >> BasicDistance1[0][0] >> BasicDistance1[1][0] >> BasicDistance1[2][0];
38  inFile >> BasicDistance1[0][1] >> BasicDistance1[1][1] >> BasicDistance1[2][1];
39  inFile >> BasicDistance1[0][2] >> BasicDistance1[1][2] >> BasicDistance1[2][2];
40  inFile >> BasicDistance2[0][0] >> BasicDistance2[1][0] >> BasicDistance2[2][0];
41  inFile >> BasicDistance2[0][1] >> BasicDistance2[1][1] >> BasicDistance2[2][1];
42  inFile >> BasicDistance2[0][2] >> BasicDistance2[1][2] >> BasicDistance2[2][2];
43  inFile >> EasyCoordinate[0] >> EasyCoordinate[1] >> EasyCoordinate[2];
44  break;
45  case 5:
46  inFile >> NumberOfPoints[0] >> NumberOfPoints[1] >> NumberOfPoints[2];
47  inFile >> ReferencePoint[0] >> ReferencePoint[1] >> ReferencePoint[2];
48  inFile >> BasicDistance0[0] >> BasicDistance0[1] >> BasicDistance0[2];
49  inFile >> RParAsFunOfPhi[0] >> RParAsFunOfPhi[1] >> RParAsFunOfPhi[2] >> RParAsFunOfPhi[3];
50  break;
51  }
52  //reading the field
53  float Bx, By, Bz;
54  BVector FieldEntry;
55  int nLines = NumberOfPoints[0]*NumberOfPoints[1]*NumberOfPoints[2];
56  FieldValues.reserve(nLines);
57 
58  for (int iLine=0; iLine<nLines; ++iLine){
59  inFile >> Bx >> By >> Bz;
60  FieldEntry.putB3(Bx,By,Bz);
61  FieldValues.push_back(FieldEntry);
62  }
63  // check completeness and close file
64  string lastEntry;
65  inFile >> lastEntry;
66  inFile.close();
67  if (lastEntry != "complete"){
68  GridType = 0;
69  cout << "error during file reading: file is not complete" << endl;
70  }
71  return;
72 }
73 
74 int MagneticFieldGrid::gridType(){
75  int type = GridType;
76  bool text = false;
77  if (text){
78  if (type == 0) cout << " grid type = " << type << " --> not determined" << endl;
79  if (type == 1) cout << " grid type = " << type << " --> (x,y,z) cube" << endl;
80  if (type == 2) cout << " grid type = " << type << " --> (x,y,z) trapezoid" << endl;
81  if (type == 3) cout << " grid type = " << type << " --> (r,phi,z) cube" << endl;
82  if (type == 4) cout << " grid type = " << type << " --> (r,phi,z) trapezoid" << endl;
83  if (type == 5) cout << " grid type = " << type << " --> (r,phi,z) 1/sin(phi)" << endl;
84  }
85  return type;
86 }
87 
88 void MagneticFieldGrid::interpolateAtPoint(double X1, double X2, double X3, float &Bx, float &By, float &Bz){
89  double dB[3] = {0.,0.,0.};
90  // define interpolation object
91  VectorFieldInterpolation MagInterpol;
92  // calculate indices for "CellPoint000"
93  int index[3];
94  putCoordGetInd(X1,X2,X3,index[0],index[1],index[2]);
95  int index0[3] = {0,0,0};
96  int index1[3] = {0,0,0};
97  for (int i=0; i<3; ++i){
98  if (NumberOfPoints[i] > 1){
99  index0[i] = max(0,index[i]);
100  if (index0[i] > NumberOfPoints[i]-2) index0[i] = NumberOfPoints[i]-2;
101  index1[i] = max(1,index[i]+1);
102  if (index1[i] > NumberOfPoints[i]-1) index1[i] = NumberOfPoints[i]-1;
103  }
104  }
105  double tmpX[3];
106  float tmpB[3];
107  // define the corners of interpolation volume
108  // FIXME: should not unpack the arrays to then repack them as first thing.
109  putIndicesGetB(index0[0],index0[1],index0[2],tmpB[0],tmpB[1],tmpB[2]);
110  putIndGetCoord(index0[0],index0[1],index0[2],tmpX[0],tmpX[1],tmpX[2]);
111  MagInterpol.defineCellPoint000(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
112  putIndicesGetB(index1[0],index0[1],index0[2],tmpB[0],tmpB[1],tmpB[2]);
113  putIndGetCoord(index1[0],index0[1],index0[2],tmpX[0],tmpX[1],tmpX[2]);
114  MagInterpol.defineCellPoint100(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
115  putIndicesGetB(index0[0],index1[1],index0[2],tmpB[0],tmpB[1],tmpB[2]);
116  putIndGetCoord(index0[0],index1[1],index0[2],tmpX[0],tmpX[1],tmpX[2]);
117  MagInterpol.defineCellPoint010(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
118  putIndicesGetB(index1[0],index1[1],index0[2],tmpB[0],tmpB[1],tmpB[2]);
119  putIndGetCoord(index1[0],index1[1],index0[2],tmpX[0],tmpX[1],tmpX[2]);
120  MagInterpol.defineCellPoint110(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
121  putIndicesGetB(index0[0],index0[1],index1[2],tmpB[0],tmpB[1],tmpB[2]);
122  putIndGetCoord(index0[0],index0[1],index1[2],tmpX[0],tmpX[1],tmpX[2]);
123  MagInterpol.defineCellPoint001(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
124  putIndicesGetB(index1[0],index0[1],index1[2],tmpB[0],tmpB[1],tmpB[2]);
125  putIndGetCoord(index1[0],index0[1],index1[2],tmpX[0],tmpX[1],tmpX[2]);
126  MagInterpol.defineCellPoint101(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
127  putIndicesGetB(index0[0],index1[1],index1[2],tmpB[0],tmpB[1],tmpB[2]);
128  putIndGetCoord(index0[0],index1[1],index1[2],tmpX[0],tmpX[1],tmpX[2]);
129  MagInterpol.defineCellPoint011(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
130  putIndicesGetB(index1[0],index1[1],index1[2],tmpB[0],tmpB[1],tmpB[2]);
131  putIndGetCoord(index1[0],index1[1],index1[2],tmpX[0],tmpX[1],tmpX[2]);
132  MagInterpol.defineCellPoint111(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
133  // interpolate
134  MagInterpol.putSCoordGetVField(X1,X2,X3,dB[0],dB[1],dB[2]);
135  Bx = float(dB[0]);
136  By = float(dB[1]);
137  Bz = float(dB[2]);
138  return;
139 }
140 
141 // FIXME: Signature should be:
142 //
143 // void MagneticFieldGrid::putCoordGetInd(double *pos, int *index)
144 //
145 void MagneticFieldGrid::putCoordGetInd(double X1, double X2, double X3, int &Index1, int &Index2, int &Index3){
146  double pnt[3] = {X1,X2,X3};
147  int index[3];
148  switch (GridType){
149  case 1:
150  for (int i=0; i<3; ++i){
151  index[i] = int((pnt[i]-ReferencePoint[i])/BasicDistance0[i]);
152  }
153  break;
154  case 2:
155  // FIXME: Should use else!
156  for (int i=0; i<3; ++i){
157  if (EasyCoordinate[i]){
158  index[i] = int((pnt[i]-ReferencePoint[i])/BasicDistance0[i]);
159  }
160  }
161  for (int i=0; i<3; ++i){
162  if (!EasyCoordinate[i]){
163  double stepSize = BasicDistance0[i];
164  double offset = 0.0;
165  for (int j=0; j<3; ++j){
166  stepSize += BasicDistance1[i][j]*index[j];
167  offset += BasicDistance2[i][j]*index[j];
168  }
169  index[i] = int((pnt[i]-(ReferencePoint[i] + offset))/stepSize);
170  }
171  }
172  break;
173  case 3:
174  for (int i=0; i<3; ++i){
175  index[i] = int((pnt[i]-ReferencePoint[i])/BasicDistance0[i]);
176  }
177  break;
178  case 4:
179  // FIXME: should use else!
180  for (int i=0; i<3; ++i){
181  if (EasyCoordinate[i]){
182  index[i] = int((pnt[i]-ReferencePoint[i])/BasicDistance0[i]);
183  }
184  }
185  for (int i=0; i<3; ++i){
186  if (!EasyCoordinate[i]){
187  double stepSize = BasicDistance0[i];
188  double offset = 0.0;
189  for (int j=0; j<3; ++j){
190  stepSize += BasicDistance1[i][j]*index[j];
191  offset += BasicDistance2[i][j]*index[j];
192  }
193  index[i] = int((pnt[i]-(ReferencePoint[i] + offset))/stepSize);
194  }
195  }
196  break;
197  case 5:
198  double sinPhi = sin(pnt[1]);
199  double stepSize = RParAsFunOfPhi[0] + RParAsFunOfPhi[1]/sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3]/sinPhi;
200  stepSize = stepSize/(NumberOfPoints[0]-1);
201  double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3]/sinPhi;
202  index[0] = int((pnt[0]-startingPoint)/stepSize);
203  index[1] = int((pnt[1]-ReferencePoint[1])/BasicDistance0[1]);
204  index[2] = int((pnt[2]-ReferencePoint[2])/BasicDistance0[2]);
205  break;
206  }
207  Index1 = index[0];
208  Index2 = index[1];
209  Index3 = index[2];
210  return;
211 }
212 
213 void MagneticFieldGrid::putIndicesGetB(int Index1, int Index2, int Index3, float &Bx, float &By, float &Bz){
214  BVector FieldEntry;
215  FieldEntry = FieldValues.operator[](lineNumber(Index1, Index2, Index3));
216  Bx = FieldEntry.bx();
217  By = FieldEntry.by();
218  Bz = FieldEntry.bz();
219  return;
220 }
221 
222 // FIXME: same as above.
223 void MagneticFieldGrid::putIndGetCoord(int Index1, int Index2, int Index3, double &X1, double &X2, double &X3){
224  int index[3] = {Index1, Index2, Index3};
225  double pnt[3];
226  switch (GridType){
227  case 1:
228  for (int i=0; i<3; ++i){
229  pnt[i] = ReferencePoint[i] + BasicDistance0[i]*index[i];
230  }
231  break;
232  case 2:
233  for (int i=0; i<3; ++i){
234  if (EasyCoordinate[i]){
235  pnt[i] = ReferencePoint[i] + BasicDistance0[i]*index[i];
236  }
237  else {
238  double stepSize = BasicDistance0[i];
239  double offset = 0.0;
240  for (int j=0; j<3; ++j){
241  stepSize += BasicDistance1[i][j]*index[j];
242  offset += BasicDistance2[i][j]*index[j];
243  }
244  pnt[i] = ReferencePoint[i] + offset + stepSize*index[i];
245  }
246  }
247  break;
248  case 3:
249  for (int i=0; i<3; ++i){
250  pnt[i] = ReferencePoint[i] + BasicDistance0[i]*index[i];
251  }
252  break;
253  case 4:
254  for (int i=0; i<3; ++i){
255  if (EasyCoordinate[i]){
256  pnt[i] = ReferencePoint[i] + BasicDistance0[i]*index[i];
257  }
258  else {
259  double stepSize = BasicDistance0[i];
260  double offset = 0.0;
261  for (int j=0; j<3; ++j){
262  stepSize += BasicDistance1[i][j]*index[j];
263  offset += BasicDistance2[i][j]*index[j];
264  }
265  pnt[i] = ReferencePoint[i] + offset + stepSize*index[i];
266  }
267  }
268  break;
269  case 5:
270  pnt[2] = ReferencePoint[2] + BasicDistance0[2]*index[2];
271  pnt[1] = ReferencePoint[1] + BasicDistance0[1]*index[1];
272  double sinPhi = sin(pnt[1]);
273  double stepSize = RParAsFunOfPhi[0] + RParAsFunOfPhi[1]/sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3]/sinPhi;
274  stepSize = stepSize/(NumberOfPoints[0]-1);
275  double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3]/sinPhi;
276  pnt[0] = startingPoint + stepSize*index[0];
277  break;
278  }
279  X1 = pnt[0];
280  X2 = pnt[1];
281  X3 = pnt[2];
282  return;
283 }
284 
285 int MagneticFieldGrid::lineNumber(int Index1, int Index2, int Index3){
286  return Index1*NumberOfPoints[1]*NumberOfPoints[2] + Index2*NumberOfPoints[2] + Index3;
287 }
288 
289 void MagneticFieldGrid::BVector::putB3(float Bx, float By, float Bz){
290  B3[0] = Bx;
291  B3[1] = By;
292  B3[2] = Bz;
293  return;
294 }
295 
296 float MagneticFieldGrid::BVector::bx(){ return B3[0]; }
297 
298 float MagneticFieldGrid::BVector::by(){ return B3[1]; }
299 
300 float MagneticFieldGrid::BVector::bz(){ return B3[2]; }
MagneticFieldGrid::load
void load(const std::string &name)
load grid binary file
Definition: MagneticFieldGrid.cc:6
type
type
Definition: HCALResponse.h:22
i
int i
Definition: DBlmapReader.cc:9
MagneticFieldGrid::putCoordGetInd
void putCoordGetInd(double X1, double X2, double X3, int &Index1, int &Index2, int &Index3)
Definition: MagneticFieldGrid.cc:145
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:289
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:74
binary_ifstream::close
void close()
Definition: binary_ifstream.cc:32
VectorFieldInterpolation::defineCellPoint110
void defineCellPoint110(double X1, double X2, double X3, double F1, double F2, double F3)
Definition: VectorFieldInterpolation.cc:38
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
MagneticFieldGrid::putIndGetCoord
void putIndGetCoord(int Index1, int Index2, int Index3, double &X1, double &X2, double &X3)
Definition: MagneticFieldGrid.cc:223
VectorFieldInterpolation::defineCellPoint101
void defineCellPoint101(double X1, double X2, double X3, double F1, double F2, double F3)
Definition: VectorFieldInterpolation.cc:60
j
int j
Definition: DBlmapReader.cc:9
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
runonSM.text
tuple text
Definition: runonSM.py:42
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:285
VectorFieldInterpolation::defineCellPoint111
void defineCellPoint111(double X1, double X2, double X3, double F1, double F2, double F3)
Definition: VectorFieldInterpolation.cc:82
evf::evtn::offset
unsigned int offset(bool)
Definition: GlobalEventNumber.cc:60
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
MagneticFieldGrid::putIndicesGetB
void putIndicesGetB(int Index1, int Index2, int Index3, float &Bx, float &By, float &Bz)
Definition: MagneticFieldGrid.cc:213
binary_ifstream
Definition: binary_ifstream.h:7
VectorFieldInterpolation::defineCellPoint001
void defineCellPoint001(double X1, double X2, double X3, double F1, double F2, double F3)
Definition: VectorFieldInterpolation.cc:49
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
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:88
Generated for CMSSW Reference Manual by doxygen 1.8.5