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