CMSSW Reference Manual

00001 // include header for MagneticFieldGrid (regular + extension for some trapezoids)
00002 #include "MagneticField/Interpolation/src/MagneticFieldGrid.h"
00003 
00004 using namespace std;
00005 
00006 void MagneticFieldGrid::load(const string& name){
00007   binary_ifstream inFile(name);
00008   inFile >> GridType;
00009   // reading the header
00010   switch (GridType){
00011   case 1:
00012     inFile >> NumberOfPoints[0]    >> NumberOfPoints[1]    >> NumberOfPoints[2];
00013     inFile >> ReferencePoint[0]    >> ReferencePoint[1]    >> ReferencePoint[2];
00014     inFile >> BasicDistance0[0]    >> BasicDistance0[1]    >> BasicDistance0[2];
00015     break;
00016   case 2:
00017     inFile >> NumberOfPoints[0]    >> NumberOfPoints[1]    >> NumberOfPoints[2];
00018     inFile >> ReferencePoint[0]    >> ReferencePoint[1]    >> ReferencePoint[2];
00019     inFile >> BasicDistance0[0]    >> BasicDistance0[1]    >> BasicDistance0[2];
00020     inFile >> BasicDistance1[0][0] >> BasicDistance1[1][0] >> BasicDistance1[2][0];
00021     inFile >> BasicDistance1[0][1] >> BasicDistance1[1][1] >> BasicDistance1[2][1];
00022     inFile >> BasicDistance1[0][2] >> BasicDistance1[1][2] >> BasicDistance1[2][2];
00023     inFile >> BasicDistance2[0][0] >> BasicDistance2[1][0] >> BasicDistance2[2][0];
00024     inFile >> BasicDistance2[0][1] >> BasicDistance2[1][1] >> BasicDistance2[2][1];
00025     inFile >> BasicDistance2[0][2] >> BasicDistance2[1][2] >> BasicDistance2[2][2];
00026     inFile >> EasyCoordinate[0]    >> EasyCoordinate[1]    >> EasyCoordinate[2];
00027     break;
00028   case 3:
00029     inFile >> NumberOfPoints[0]    >> NumberOfPoints[1]    >> NumberOfPoints[2];
00030     inFile >> ReferencePoint[0]    >> ReferencePoint[1]    >> ReferencePoint[2];
00031     inFile >> BasicDistance0[0]    >> BasicDistance0[1]    >> BasicDistance0[2];
00032     break;
00033   case 4:
00034     inFile >> NumberOfPoints[0]    >> NumberOfPoints[1]    >> NumberOfPoints[2];
00035     inFile >> ReferencePoint[0]    >> ReferencePoint[1]    >> ReferencePoint[2];
00036     inFile >> BasicDistance0[0]    >> BasicDistance0[1]    >> BasicDistance0[2];
00037     inFile >> BasicDistance1[0][0] >> BasicDistance1[1][0] >> BasicDistance1[2][0];
00038     inFile >> BasicDistance1[0][1] >> BasicDistance1[1][1] >> BasicDistance1[2][1];
00039     inFile >> BasicDistance1[0][2] >> BasicDistance1[1][2] >> BasicDistance1[2][2];
00040     inFile >> BasicDistance2[0][0] >> BasicDistance2[1][0] >> BasicDistance2[2][0];
00041     inFile >> BasicDistance2[0][1] >> BasicDistance2[1][1] >> BasicDistance2[2][1];
00042     inFile >> BasicDistance2[0][2] >> BasicDistance2[1][2] >> BasicDistance2[2][2];
00043     inFile >> EasyCoordinate[0]    >> EasyCoordinate[1]    >> EasyCoordinate[2];
00044     break;
00045   case 5:
00046     inFile >> NumberOfPoints[0]    >> NumberOfPoints[1]    >> NumberOfPoints[2];
00047     inFile >> ReferencePoint[0]    >> ReferencePoint[1]    >> ReferencePoint[2];
00048     inFile >> BasicDistance0[0]    >> BasicDistance0[1]    >> BasicDistance0[2];
00049     inFile >> RParAsFunOfPhi[0]    >> RParAsFunOfPhi[1]    >> RParAsFunOfPhi[2]    >> RParAsFunOfPhi[3];
00050     break;
00051   }
00052   //reading the field
00053   float Bx, By, Bz;
00054   BVector FieldEntry;
00055   int nLines = NumberOfPoints[0]*NumberOfPoints[1]*NumberOfPoints[2];
00056   for (int iLine=0; iLine<nLines; ++iLine){
00057     inFile >> Bx >> By >> Bz;
00058     FieldEntry.putB3(Bx,By,Bz);
00059     FieldValues.push_back(FieldEntry);
00060   }
00061   // check completeness and close file
00062   string lastEntry;
00063   inFile >> lastEntry;
00064   inFile.close();
00065   if (lastEntry != "complete"){
00066     GridType = 0;
00067     cout << "error during file reading: file is not complete" << endl;
00068   }
00069   return;
00070 }
00071 
00072 int MagneticFieldGrid::gridType(){
00073   int type = GridType;
00074   bool text = false;
00075   if (text){
00076     if (type == 0) cout << "  grid type = " << type << "  -->  not determined" << endl;
00077     if (type == 1) cout << "  grid type = " << type << "  -->  (x,y,z) cube" << endl;
00078     if (type == 2) cout << "  grid type = " << type << "  -->  (x,y,z) trapezoid" << endl;
00079     if (type == 3) cout << "  grid type = " << type << "  -->  (r,phi,z) cube" << endl;
00080     if (type == 4) cout << "  grid type = " << type << "  -->  (r,phi,z) trapezoid" << endl;
00081     if (type == 5) cout << "  grid type = " << type << "  -->  (r,phi,z) 1/sin(phi)" << endl;
00082   }
00083   return type;
00084 }
00085 
00086 void MagneticFieldGrid::interpolateAtPoint(double X1, double X2, double X3, float &Bx, float &By, float &Bz){
00087   double dB[3] = {0.,0.,0.};
00088   // define interpolation object
00089   VectorFieldInterpolation MagInterpol;
00090   // calculate indices for "CellPoint000"
00091   int index[3];
00092   putCoordGetInd(X1,X2,X3,index[0],index[1],index[2]);
00093   int index0[3] = {0,0,0};
00094   int index1[3] = {0,0,0};
00095   for (int i=0; i<3; ++i){
00096     if (NumberOfPoints[i] > 1){
00097                                            index0[i] = max(0,index[i]);
00098       if (index0[i] > NumberOfPoints[i]-2) index0[i] = NumberOfPoints[i]-2;
00099                                            index1[i] = max(1,index[i]+1);
00100       if (index1[i] > NumberOfPoints[i]-1) index1[i] = NumberOfPoints[i]-1;
00101     }
00102   }
00103   double tmpX[3];
00104   float  tmpB[3];
00105   // define the corners of interpolation volume
00106   putIndicesGetB(index0[0],index0[1],index0[2],tmpB[0],tmpB[1],tmpB[2]);
00107   putIndGetCoord(index0[0],index0[1],index0[2],tmpX[0],tmpX[1],tmpX[2]);
00108   MagInterpol.defineCellPoint000(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
00109   putIndicesGetB(index1[0],index0[1],index0[2],tmpB[0],tmpB[1],tmpB[2]);
00110   putIndGetCoord(index1[0],index0[1],index0[2],tmpX[0],tmpX[1],tmpX[2]);
00111   MagInterpol.defineCellPoint100(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
00112   putIndicesGetB(index0[0],index1[1],index0[2],tmpB[0],tmpB[1],tmpB[2]);
00113   putIndGetCoord(index0[0],index1[1],index0[2],tmpX[0],tmpX[1],tmpX[2]);
00114   MagInterpol.defineCellPoint010(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
00115   putIndicesGetB(index1[0],index1[1],index0[2],tmpB[0],tmpB[1],tmpB[2]);
00116   putIndGetCoord(index1[0],index1[1],index0[2],tmpX[0],tmpX[1],tmpX[2]);
00117   MagInterpol.defineCellPoint110(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
00118   putIndicesGetB(index0[0],index0[1],index1[2],tmpB[0],tmpB[1],tmpB[2]);
00119   putIndGetCoord(index0[0],index0[1],index1[2],tmpX[0],tmpX[1],tmpX[2]);
00120   MagInterpol.defineCellPoint001(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
00121   putIndicesGetB(index1[0],index0[1],index1[2],tmpB[0],tmpB[1],tmpB[2]);
00122   putIndGetCoord(index1[0],index0[1],index1[2],tmpX[0],tmpX[1],tmpX[2]);
00123   MagInterpol.defineCellPoint101(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
00124   putIndicesGetB(index0[0],index1[1],index1[2],tmpB[0],tmpB[1],tmpB[2]);
00125   putIndGetCoord(index0[0],index1[1],index1[2],tmpX[0],tmpX[1],tmpX[2]);
00126   MagInterpol.defineCellPoint011(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
00127   putIndicesGetB(index1[0],index1[1],index1[2],tmpB[0],tmpB[1],tmpB[2]);
00128   putIndGetCoord(index1[0],index1[1],index1[2],tmpX[0],tmpX[1],tmpX[2]);
00129   MagInterpol.defineCellPoint111(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
00130   // interpolate
00131   MagInterpol.putSCoordGetVField(X1,X2,X3,dB[0],dB[1],dB[2]);
00132   Bx = float(dB[0]);
00133   By = float(dB[1]);
00134   Bz = float(dB[2]);
00135   return;
00136 }
00137 
00138 void MagneticFieldGrid::putCoordGetInd(double X1, double X2, double X3, int &Index1, int &Index2, int &Index3){
00139   double pnt[3] = {X1,X2,X3};
00140   int index[3];
00141   switch (GridType){
00142   case 1:
00143     for (int i=0; i<3; ++i){
00144       index[i] = int((pnt[i]-ReferencePoint[i])/BasicDistance0[i]);
00145     }
00146     break;
00147   case 2:
00148     for (int i=0; i<3; ++i){
00149       if (EasyCoordinate[i]){
00150         index[i] = int((pnt[i]-ReferencePoint[i])/BasicDistance0[i]);
00151       }
00152     }
00153     for (int i=0; i<3; ++i){
00154       if (!EasyCoordinate[i]){
00155         double stepSize = BasicDistance0[i];
00156         double offset   = 0.0;
00157         for (int j=0; j<3; ++j){
00158           stepSize += BasicDistance1[i][j]*index[j];
00159           offset   += BasicDistance2[i][j]*index[j];
00160         }
00161         index[i] = int((pnt[i]-(ReferencePoint[i] + offset))/stepSize);
00162       }
00163     }
00164     break;
00165   case 3:
00166     for (int i=0; i<3; ++i){
00167       index[i] = int((pnt[i]-ReferencePoint[i])/BasicDistance0[i]);
00168     }
00169     break;
00170   case 4:
00171     for (int i=0; i<3; ++i){
00172       if (EasyCoordinate[i]){
00173         index[i] = int((pnt[i]-ReferencePoint[i])/BasicDistance0[i]);
00174       }
00175     }
00176     for (int i=0; i<3; ++i){
00177       if (!EasyCoordinate[i]){
00178         double stepSize = BasicDistance0[i];
00179         double offset   = 0.0;
00180         for (int j=0; j<3; ++j){
00181           stepSize += BasicDistance1[i][j]*index[j];
00182           offset   += BasicDistance2[i][j]*index[j];
00183         }
00184         index[i] = int((pnt[i]-(ReferencePoint[i] + offset))/stepSize);
00185       }
00186     }
00187     break;
00188   case 5:
00189     double sinPhi = sin(pnt[1]);
00190     double stepSize = RParAsFunOfPhi[0] + RParAsFunOfPhi[1]/sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3]/sinPhi;
00191     stepSize =  stepSize/(NumberOfPoints[0]-1);
00192     double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3]/sinPhi;
00193     index[0] = int((pnt[0]-startingPoint)/stepSize);
00194     index[1] = int((pnt[1]-ReferencePoint[1])/BasicDistance0[1]);
00195     index[2] = int((pnt[2]-ReferencePoint[2])/BasicDistance0[2]);
00196     break;
00197   }
00198   Index1 = index[0];
00199   Index2 = index[1];
00200   Index3 = index[2];
00201   return;
00202 }
00203 
00204 void MagneticFieldGrid::putIndicesGetB(int Index1, int Index2, int Index3, float &Bx, float &By, float &Bz){
00205   BVector FieldEntry;
00206   FieldEntry = FieldValues.operator[](lineNumber(Index1, Index2, Index3));
00207   Bx = FieldEntry.bx();
00208   By = FieldEntry.by();
00209   Bz = FieldEntry.bz();
00210   return;
00211 }
00212 
00213 void MagneticFieldGrid::putIndGetCoord(int Index1, int Index2, int Index3, double &X1, double &X2, double &X3){
00214   int index[3] = {Index1, Index2, Index3};
00215   double pnt[3];
00216   switch (GridType){
00217   case 1:
00218     for (int i=0; i<3; ++i){
00219       pnt[i] = ReferencePoint[i] + BasicDistance0[i]*index[i];
00220     }
00221     break;
00222   case 2:
00223     for (int i=0; i<3; ++i){
00224       if (EasyCoordinate[i]){
00225         pnt[i] = ReferencePoint[i] + BasicDistance0[i]*index[i];
00226       }
00227       else {
00228         double stepSize = BasicDistance0[i];
00229         double offset   = 0.0;
00230         for (int j=0; j<3; ++j){
00231           stepSize += BasicDistance1[i][j]*index[j];
00232           offset   += BasicDistance2[i][j]*index[j];
00233         }
00234         pnt[i] = ReferencePoint[i] + offset + stepSize*index[i];
00235       }
00236     }
00237     break;
00238   case 3:
00239     for (int i=0; i<3; ++i){
00240       pnt[i] = ReferencePoint[i] + BasicDistance0[i]*index[i];
00241     }
00242     break;
00243   case 4:
00244     for (int i=0; i<3; ++i){
00245       if (EasyCoordinate[i]){
00246         pnt[i] = ReferencePoint[i] + BasicDistance0[i]*index[i];
00247       }
00248       else {
00249         double stepSize = BasicDistance0[i];
00250         double offset   = 0.0;
00251         for (int j=0; j<3; ++j){
00252           stepSize += BasicDistance1[i][j]*index[j];
00253           offset   += BasicDistance2[i][j]*index[j];
00254         }
00255         pnt[i] = ReferencePoint[i] + offset + stepSize*index[i];
00256       }
00257     }
00258     break;
00259   case 5:
00260     pnt[2] = ReferencePoint[2] + BasicDistance0[2]*index[2];
00261     pnt[1] = ReferencePoint[1] + BasicDistance0[1]*index[1];
00262     double sinPhi = sin(pnt[1]);
00263     double stepSize = RParAsFunOfPhi[0] + RParAsFunOfPhi[1]/sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3]/sinPhi;
00264     stepSize =  stepSize/(NumberOfPoints[0]-1);
00265     double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3]/sinPhi;
00266     pnt[0] = startingPoint + stepSize*index[0];
00267     break;
00268   }
00269   X1 = pnt[0];
00270   X2 = pnt[1];
00271   X3 = pnt[2];
00272   return;
00273 }
00274 
00275 int MagneticFieldGrid::lineNumber(int Index1, int Index2, int Index3){
00276   return Index1*NumberOfPoints[1]*NumberOfPoints[2] + Index2*NumberOfPoints[2] + Index3;
00277 }
00278 
00279 void MagneticFieldGrid::BVector::putB3(float Bx, float By, float Bz){
00280   B3[0] = Bx;
00281   B3[1] = By;
00282   B3[2] = Bz;
00283   return;
00284 }
00285 
00286 float MagneticFieldGrid::BVector::bx(){  return B3[0]; }
00287 
00288 float MagneticFieldGrid::BVector::by(){  return B3[1]; }
00289 
00290 float MagneticFieldGrid::BVector::bz(){  return B3[2]; }
MagneticFieldGrid.cc
