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)
82  cout << " grid type = " << type << " --> not determined" << endl;
83  if (type == 1)
84  cout << " grid type = " << type << " --> (x,y,z) cube" << endl;
85  if (type == 2)
86  cout << " grid type = " << type << " --> (x,y,z) trapezoid" << endl;
87  if (type == 3)
88  cout << " grid type = " << type << " --> (r,phi,z) cube" << endl;
89  if (type == 4)
90  cout << " grid type = " << type << " --> (r,phi,z) trapezoid" << endl;
91  if (type == 5)
92  cout << " grid type = " << type << " --> (r,phi,z) 1/sin(phi)" << endl;
93  }
94  return type;
95 }
96 
97 void MagneticFieldGrid::interpolateAtPoint(double X1, double X2, double X3, float &Bx, float &By, float &Bz) {
98  double dB[3] = {0., 0., 0.};
99  // define interpolation object
100  VectorFieldInterpolation MagInterpol;
101  // calculate indices for "CellPoint000"
102  int index[3];
103  putCoordGetInd(X1, X2, X3, index[0], index[1], index[2]);
104  int index0[3] = {0, 0, 0};
105  int index1[3] = {0, 0, 0};
106  for (int i = 0; i < 3; ++i) {
107  if (NumberOfPoints[i] > 1) {
108  index0[i] = max(0, index[i]);
109  if (index0[i] > NumberOfPoints[i] - 2)
110  index0[i] = NumberOfPoints[i] - 2;
111  index1[i] = max(1, index[i] + 1);
112  if (index1[i] > NumberOfPoints[i] - 1)
113  index1[i] = NumberOfPoints[i] - 1;
114  }
115  }
116  double tmpX[3];
117  float tmpB[3];
118  // define the corners of interpolation volume
119  // FIXME: should not unpack the arrays to then repack them as first thing.
120  putIndicesGetB(index0[0], index0[1], index0[2], tmpB[0], tmpB[1], tmpB[2]);
121  putIndGetCoord(index0[0], index0[1], index0[2], tmpX[0], tmpX[1], tmpX[2]);
122  MagInterpol.defineCellPoint000(tmpX[0], tmpX[1], tmpX[2], double(tmpB[0]), double(tmpB[1]), double(tmpB[2]));
123  putIndicesGetB(index1[0], index0[1], index0[2], tmpB[0], tmpB[1], tmpB[2]);
124  putIndGetCoord(index1[0], index0[1], index0[2], tmpX[0], tmpX[1], tmpX[2]);
125  MagInterpol.defineCellPoint100(tmpX[0], tmpX[1], tmpX[2], double(tmpB[0]), double(tmpB[1]), double(tmpB[2]));
126  putIndicesGetB(index0[0], index1[1], index0[2], tmpB[0], tmpB[1], tmpB[2]);
127  putIndGetCoord(index0[0], index1[1], index0[2], tmpX[0], tmpX[1], tmpX[2]);
128  MagInterpol.defineCellPoint010(tmpX[0], tmpX[1], tmpX[2], double(tmpB[0]), double(tmpB[1]), double(tmpB[2]));
129  putIndicesGetB(index1[0], index1[1], index0[2], tmpB[0], tmpB[1], tmpB[2]);
130  putIndGetCoord(index1[0], index1[1], index0[2], tmpX[0], tmpX[1], tmpX[2]);
131  MagInterpol.defineCellPoint110(tmpX[0], tmpX[1], tmpX[2], double(tmpB[0]), double(tmpB[1]), double(tmpB[2]));
132  putIndicesGetB(index0[0], index0[1], index1[2], tmpB[0], tmpB[1], tmpB[2]);
133  putIndGetCoord(index0[0], index0[1], index1[2], tmpX[0], tmpX[1], tmpX[2]);
134  MagInterpol.defineCellPoint001(tmpX[0], tmpX[1], tmpX[2], double(tmpB[0]), double(tmpB[1]), double(tmpB[2]));
135  putIndicesGetB(index1[0], index0[1], index1[2], tmpB[0], tmpB[1], tmpB[2]);
136  putIndGetCoord(index1[0], index0[1], index1[2], tmpX[0], tmpX[1], tmpX[2]);
137  MagInterpol.defineCellPoint101(tmpX[0], tmpX[1], tmpX[2], double(tmpB[0]), double(tmpB[1]), double(tmpB[2]));
138  putIndicesGetB(index0[0], index1[1], index1[2], tmpB[0], tmpB[1], tmpB[2]);
139  putIndGetCoord(index0[0], index1[1], index1[2], tmpX[0], tmpX[1], tmpX[2]);
140  MagInterpol.defineCellPoint011(tmpX[0], tmpX[1], tmpX[2], double(tmpB[0]), double(tmpB[1]), double(tmpB[2]));
141  putIndicesGetB(index1[0], index1[1], index1[2], tmpB[0], tmpB[1], tmpB[2]);
142  putIndGetCoord(index1[0], index1[1], index1[2], tmpX[0], tmpX[1], tmpX[2]);
143  MagInterpol.defineCellPoint111(tmpX[0], tmpX[1], tmpX[2], double(tmpB[0]), double(tmpB[1]), double(tmpB[2]));
144  // interpolate
145  MagInterpol.putSCoordGetVField(X1, X2, X3, dB[0], dB[1], dB[2]);
146  Bx = float(dB[0]);
147  By = float(dB[1]);
148  Bz = float(dB[2]);
149  return;
150 }
151 
152 // FIXME: Signature should be:
153 //
154 // void MagneticFieldGrid::putCoordGetInd(double *pos, int *index)
155 //
156 void MagneticFieldGrid::putCoordGetInd(double X1, double X2, double X3, int &Index1, int &Index2, int &Index3) {
157  double pnt[3] = {X1, X2, X3};
158  int index[3];
159  switch (GridType) {
160  case 1: {
161  for (int i = 0; i < 3; ++i) {
162  index[i] = int((pnt[i] - ReferencePoint[i]) / BasicDistance0[i]);
163  }
164  break;
165  }
166  case 2: {
167  // FIXME: Should use else!
168  for (int i = 0; i < 3; ++i) {
169  if (EasyCoordinate[i]) {
170  index[i] = int((pnt[i] - ReferencePoint[i]) / BasicDistance0[i]);
171  } else
172  index[i] = 0; //computed below
173  }
174  for (int i = 0; i < 3; ++i) {
175  if (!EasyCoordinate[i]) {
176  double stepSize = BasicDistance0[i];
177  double offset = 0.0;
178  for (int j = 0; j < 3; ++j) {
179  stepSize += BasicDistance1[i][j] * index[j];
180  offset += BasicDistance2[i][j] * index[j];
181  }
182  index[i] = int((pnt[i] - (ReferencePoint[i] + offset)) / stepSize);
183  }
184  }
185  break;
186  }
187  case 3: {
188  for (int i = 0; i < 3; ++i) {
189  index[i] = int((pnt[i] - ReferencePoint[i]) / BasicDistance0[i]);
190  }
191  break;
192  }
193  case 4: {
194  // FIXME: should use else!
195  for (int i = 0; i < 3; ++i) {
196  if (EasyCoordinate[i]) {
197  index[i] = int((pnt[i] - ReferencePoint[i]) / BasicDistance0[i]);
198  } else
199  index[i] = 0; //computed below
200  }
201  for (int i = 0; i < 3; ++i) {
202  if (!EasyCoordinate[i]) {
203  double stepSize = BasicDistance0[i];
204  double offset = 0.0;
205  for (int j = 0; j < 3; ++j) {
206  stepSize += BasicDistance1[i][j] * index[j];
207  offset += BasicDistance2[i][j] * index[j];
208  }
209  index[i] = int((pnt[i] - (ReferencePoint[i] + offset)) / stepSize);
210  }
211  }
212  break;
213  }
214  case 5: {
215  double sinPhi = sin(pnt[1]);
216  double stepSize = RParAsFunOfPhi[0] + RParAsFunOfPhi[1] / sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3] / sinPhi;
217  stepSize = stepSize / (NumberOfPoints[0] - 1);
218  double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3] / sinPhi;
219  index[0] = int((pnt[0] - startingPoint) / stepSize);
220  index[1] = int((pnt[1] - ReferencePoint[1]) / BasicDistance0[1]);
221  index[2] = int((pnt[2] - ReferencePoint[2]) / BasicDistance0[2]);
222  break;
223  }
224  default:
225  assert(0); //shouldn't be here
226  }
227  Index1 = index[0];
228  Index2 = index[1];
229  Index3 = index[2];
230  return;
231 }
232 
233 void MagneticFieldGrid::putIndicesGetB(int Index1, int Index2, int Index3, float &Bx, float &By, float &Bz) {
234  BVector FieldEntry;
235  FieldEntry = FieldValues.operator[](lineNumber(Index1, Index2, Index3));
236  Bx = FieldEntry.bx();
237  By = FieldEntry.by();
238  Bz = FieldEntry.bz();
239  return;
240 }
241 
242 // FIXME: same as above.
243 void MagneticFieldGrid::putIndGetCoord(int Index1, int Index2, int Index3, double &X1, double &X2, double &X3) {
244  int index[3] = {Index1, Index2, Index3};
245  double pnt[3];
246  switch (GridType) {
247  case 1: {
248  for (int i = 0; i < 3; ++i) {
249  pnt[i] = ReferencePoint[i] + BasicDistance0[i] * index[i];
250  }
251  break;
252  }
253  case 2: {
254  for (int i = 0; i < 3; ++i) {
255  if (EasyCoordinate[i]) {
256  pnt[i] = ReferencePoint[i] + BasicDistance0[i] * index[i];
257  } else {
258  double stepSize = BasicDistance0[i];
259  double offset = 0.0;
260  for (int j = 0; j < 3; ++j) {
261  stepSize += BasicDistance1[i][j] * index[j];
262  offset += BasicDistance2[i][j] * index[j];
263  }
264  pnt[i] = ReferencePoint[i] + offset + stepSize * index[i];
265  }
266  }
267  break;
268  }
269  case 3: {
270  for (int i = 0; i < 3; ++i) {
271  pnt[i] = ReferencePoint[i] + BasicDistance0[i] * index[i];
272  }
273  break;
274  }
275  case 4: {
276  for (int i = 0; i < 3; ++i) {
277  if (EasyCoordinate[i]) {
278  pnt[i] = ReferencePoint[i] + BasicDistance0[i] * index[i];
279  } else {
280  double stepSize = BasicDistance0[i];
281  double offset = 0.0;
282  for (int j = 0; j < 3; ++j) {
283  stepSize += BasicDistance1[i][j] * index[j];
284  offset += BasicDistance2[i][j] * index[j];
285  }
286  pnt[i] = ReferencePoint[i] + offset + stepSize * index[i];
287  }
288  }
289  break;
290  }
291  case 5: {
292  pnt[2] = ReferencePoint[2] + BasicDistance0[2] * index[2];
293  pnt[1] = ReferencePoint[1] + BasicDistance0[1] * index[1];
294  double sinPhi = sin(pnt[1]);
295  double stepSize = RParAsFunOfPhi[0] + RParAsFunOfPhi[1] / sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3] / sinPhi;
296  stepSize = stepSize / (NumberOfPoints[0] - 1);
297  double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3] / sinPhi;
298  pnt[0] = startingPoint + stepSize * index[0];
299  break;
300  }
301  default:
302  assert(0); //bug if make it here
303  }
304  X1 = pnt[0];
305  X2 = pnt[1];
306  X3 = pnt[2];
307  return;
308 }
309 
310 int MagneticFieldGrid::lineNumber(int Index1, int Index2, int Index3) {
311  return Index1 * NumberOfPoints[1] * NumberOfPoints[2] + Index2 * NumberOfPoints[2] + Index3;
312 }
313 
314 void MagneticFieldGrid::BVector::putB3(float Bx, float By, float Bz) {
315  B3[0] = Bx;
316  B3[1] = By;
317  B3[2] = Bz;
318  return;
319 }
320 
321 float MagneticFieldGrid::BVector::bx() { return B3[0]; }
322 
323 float MagneticFieldGrid::BVector::by() { return B3[1]; }
324 
325 float MagneticFieldGrid::BVector::bz() { return B3[2]; }
MagneticFieldGrid::gridType
int gridType()
returns value of GridType (and eventually prints the type + short description)
Definition: MagneticFieldGrid.cc:77
MagneticFieldGrid::putIndicesGetB
void putIndicesGetB(int Index1, int Index2, int Index3, float &Bx, float &By, float &Bz)
Definition: MagneticFieldGrid.cc:233
VectorFieldInterpolation
Definition: VectorFieldInterpolation.h:54
mps_fire.i
i
Definition: mps_fire.py:428
dqmMemoryStats.float
float
Definition: dqmMemoryStats.py:127
VectorFieldInterpolation::defineCellPoint011
void defineCellPoint011(double X1, double X2, double X3, double F1, double F2, double F3)
Definition: VectorFieldInterpolation.cc:64
MagneticFieldGrid::BVector::bz
float bz()
Definition: MagneticFieldGrid.cc:325
gather_cfg.cout
cout
Definition: gather_cfg.py:144
cms::cuda::assert
assert(be >=bs)
submitDQMOfflineCAF.nLines
nLines
Definition: submitDQMOfflineCAF.py:676
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
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::defineCellPoint100
void defineCellPoint100(double X1, double X2, double X3, double F1, double F2, double F3)
Definition: VectorFieldInterpolation.cc:14
VectorFieldInterpolation::defineCellPoint111
void defineCellPoint111(double X1, double X2, double X3, double F1, double F2, double F3)
Definition: VectorFieldInterpolation.cc:74
MagneticFieldGrid::lineNumber
int lineNumber(int Index1, int Index2, int Index3)
Definition: MagneticFieldGrid.cc:310
MagneticFieldGrid::putCoordGetInd
void putCoordGetInd(double X1, double X2, double X3, int &Index1, int &Index2, int &Index3)
Definition: MagneticFieldGrid.cc:156
Phase1L1TJetProducer_cfi.sinPhi
sinPhi
Definition: Phase1L1TJetProducer_cfi.py:39
binary_ifstream
Definition: binary_ifstream.h:8
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
VectorFieldInterpolation::defineCellPoint110
void defineCellPoint110(double X1, double X2, double X3, double F1, double F2, double F3)
Definition: VectorFieldInterpolation.cc:34
MagneticFieldGrid.h
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
type
type
Definition: SiPixelVCal_PayloadInspector.cc:39
gainCalibHelper::gainCalibPI::type
type
Definition: SiPixelGainCalibHelper.h:40
MagneticFieldGrid::BVector::putB3
void putB3(float Bx, float By, float Bz)
Definition: MagneticFieldGrid.cc:314
createfilelist.int
int
Definition: createfilelist.py:10
MagneticFieldGrid::BVector::bx
float bx()
Definition: MagneticFieldGrid.cc:321
binary_ifstream::close
void close()
Definition: binary_ifstream.cc:21
MagneticFieldGrid::BVector
Definition: MagneticFieldGrid.h:80
VectorFieldInterpolation::defineCellPoint101
void defineCellPoint101(double X1, double X2, double X3, double F1, double F2, double F3)
Definition: VectorFieldInterpolation.cc:54
MagneticFieldGrid::putIndGetCoord
void putIndGetCoord(int Index1, int Index2, int Index3, double &X1, double &X2, double &X3)
Definition: MagneticFieldGrid.cc:243
VectorFieldInterpolation::defineCellPoint010
void defineCellPoint010(double X1, double X2, double X3, double F1, double F2, double F3)
Definition: VectorFieldInterpolation.cc:24
std
Definition: JetResolutionObject.h:76
MagneticFieldGrid::load
void load(const std::string &name)
load grid binary file
Definition: MagneticFieldGrid.cc:7
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:97
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
AlignmentPI::index
index
Definition: AlignmentPayloadInspectorHelper.h:46
MagneticFieldGrid::BVector::by
float by()
Definition: MagneticFieldGrid.cc:323
runonSM.text
text
Definition: runonSM.py:43
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
hltrates_dqm_sourceclient-live_cfg.offset
offset
Definition: hltrates_dqm_sourceclient-live_cfg.py:82
VectorFieldInterpolation::defineCellPoint001
void defineCellPoint001(double X1, double X2, double X3, double F1, double F2, double F3)
Definition: VectorFieldInterpolation.cc:44
Generated for CMSSW Reference Manual by doxygen 1.8.16