Reference Manual

Go to the documentation of this file.
00001 // include header for MagneticFieldGrid (regular + extension for some trapezoids)
00002 #include "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   FieldValues.reserve(nLines);
00057 
00058   for (int iLine=0; iLine<nLines; ++iLine){
00059     inFile >> Bx >> By >> Bz;
00060     FieldEntry.putB3(Bx,By,Bz);
00061     FieldValues.push_back(FieldEntry);
00062   }
00063   // check completeness and close file
00064   string lastEntry;
00065   inFile >> lastEntry;
00066   inFile.close();
00067   if (lastEntry != "complete"){
00068     GridType = 0;
00069     cout << "error during file reading: file is not complete" << endl;
00070   }
00071   return;
00072 }
00073 
00074 int MagneticFieldGrid::gridType(){
00075   int type = GridType;
00076   bool text = false;
00077   if (text){
00078     if (type == 0) cout << "  grid type = " << type << "  -->  not determined" << endl;
00079     if (type == 1) cout << "  grid type = " << type << "  -->  (x,y,z) cube" << endl;
00080     if (type == 2) cout << "  grid type = " << type << "  -->  (x,y,z) trapezoid" << endl;
00081     if (type == 3) cout << "  grid type = " << type << "  -->  (r,phi,z) cube" << endl;
00082     if (type == 4) cout << "  grid type = " << type << "  -->  (r,phi,z) trapezoid" << endl;
00083     if (type == 5) cout << "  grid type = " << type << "  -->  (r,phi,z) 1/sin(phi)" << endl;
00084   }
00085   return type;
00086 }
00087 
00088 void MagneticFieldGrid::interpolateAtPoint(double X1, double X2, double X3, float &Bx, float &By, float &Bz){
00089   double dB[3] = {0.,0.,0.};
00090   // define interpolation object
00091   VectorFieldInterpolation MagInterpol;
00092   // calculate indices for "CellPoint000"
00093   int index[3];
00094   putCoordGetInd(X1,X2,X3,index[0],index[1],index[2]);
00095   int index0[3] = {0,0,0};
00096   int index1[3] = {0,0,0};
00097   for (int i=0; i<3; ++i){
00098     if (NumberOfPoints[i] > 1){
00099                                            index0[i] = max(0,index[i]);
00100       if (index0[i] > NumberOfPoints[i]-2) index0[i] = NumberOfPoints[i]-2;
00101                                            index1[i] = max(1,index[i]+1);
00102       if (index1[i] > NumberOfPoints[i]-1) index1[i] = NumberOfPoints[i]-1;
00103     }
00104   }
00105   double tmpX[3];
00106   float  tmpB[3];
00107   // define the corners of interpolation volume
00108   // FIXME: should not unpack the arrays to then repack them as first thing.
00109   putIndicesGetB(index0[0],index0[1],index0[2],tmpB[0],tmpB[1],tmpB[2]);
00110   putIndGetCoord(index0[0],index0[1],index0[2],tmpX[0],tmpX[1],tmpX[2]);
00111   MagInterpol.defineCellPoint000(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
00112   putIndicesGetB(index1[0],index0[1],index0[2],tmpB[0],tmpB[1],tmpB[2]);
00113   putIndGetCoord(index1[0],index0[1],index0[2],tmpX[0],tmpX[1],tmpX[2]);
00114   MagInterpol.defineCellPoint100(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
00115   putIndicesGetB(index0[0],index1[1],index0[2],tmpB[0],tmpB[1],tmpB[2]);
00116   putIndGetCoord(index0[0],index1[1],index0[2],tmpX[0],tmpX[1],tmpX[2]);
00117   MagInterpol.defineCellPoint010(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
00118   putIndicesGetB(index1[0],index1[1],index0[2],tmpB[0],tmpB[1],tmpB[2]);
00119   putIndGetCoord(index1[0],index1[1],index0[2],tmpX[0],tmpX[1],tmpX[2]);
00120   MagInterpol.defineCellPoint110(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
00121   putIndicesGetB(index0[0],index0[1],index1[2],tmpB[0],tmpB[1],tmpB[2]);
00122   putIndGetCoord(index0[0],index0[1],index1[2],tmpX[0],tmpX[1],tmpX[2]);
00123   MagInterpol.defineCellPoint001(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
00124   putIndicesGetB(index1[0],index0[1],index1[2],tmpB[0],tmpB[1],tmpB[2]);
00125   putIndGetCoord(index1[0],index0[1],index1[2],tmpX[0],tmpX[1],tmpX[2]);
00126   MagInterpol.defineCellPoint101(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
00127   putIndicesGetB(index0[0],index1[1],index1[2],tmpB[0],tmpB[1],tmpB[2]);
00128   putIndGetCoord(index0[0],index1[1],index1[2],tmpX[0],tmpX[1],tmpX[2]);
00129   MagInterpol.defineCellPoint011(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
00130   putIndicesGetB(index1[0],index1[1],index1[2],tmpB[0],tmpB[1],tmpB[2]);
00131   putIndGetCoord(index1[0],index1[1],index1[2],tmpX[0],tmpX[1],tmpX[2]);
00132   MagInterpol.defineCellPoint111(tmpX[0],tmpX[1],tmpX[2],double(tmpB[0]),double(tmpB[1]),double(tmpB[2]));
00133   // interpolate
00134   MagInterpol.putSCoordGetVField(X1,X2,X3,dB[0],dB[1],dB[2]);
00135   Bx = float(dB[0]);
00136   By = float(dB[1]);
00137   Bz = float(dB[2]);
00138   return;
00139 }
00140 
00141 // FIXME: Signature should be:
00142 //
00143 //     void MagneticFieldGrid::putCoordGetInd(double *pos, int *index) 
00144 //
00145 void MagneticFieldGrid::putCoordGetInd(double X1, double X2, double X3, int &Index1, int &Index2, int &Index3){
00146   double pnt[3] = {X1,X2,X3};
00147   int index[3];
00148   switch (GridType){
00149   case 1:
00150     for (int i=0; i<3; ++i){
00151       index[i] = int((pnt[i]-ReferencePoint[i])/BasicDistance0[i]);
00152     }
00153     break;
00154   case 2:
00155     // FIXME: Should use else!
00156     for (int i=0; i<3; ++i){
00157       if (EasyCoordinate[i]){
00158         index[i] = int((pnt[i]-ReferencePoint[i])/BasicDistance0[i]);
00159       }
00160     }
00161     for (int i=0; i<3; ++i){
00162       if (!EasyCoordinate[i]){
00163         double stepSize = BasicDistance0[i];
00164         double offset   = 0.0;
00165         for (int j=0; j<3; ++j){
00166           stepSize += BasicDistance1[i][j]*index[j];
00167           offset   += BasicDistance2[i][j]*index[j];
00168         }
00169         index[i] = int((pnt[i]-(ReferencePoint[i] + offset))/stepSize);
00170       }
00171     }
00172     break;
00173   case 3:
00174     for (int i=0; i<3; ++i){
00175       index[i] = int((pnt[i]-ReferencePoint[i])/BasicDistance0[i]);
00176     }
00177     break;
00178   case 4:
00179     // FIXME: should use else!
00180     for (int i=0; i<3; ++i){
00181       if (EasyCoordinate[i]){
00182         index[i] = int((pnt[i]-ReferencePoint[i])/BasicDistance0[i]);
00183       }
00184     }
00185     for (int i=0; i<3; ++i){
00186       if (!EasyCoordinate[i]){
00187         double stepSize = BasicDistance0[i];
00188         double offset   = 0.0;
00189         for (int j=0; j<3; ++j){
00190           stepSize += BasicDistance1[i][j]*index[j];
00191           offset   += BasicDistance2[i][j]*index[j];
00192         }
00193         index[i] = int((pnt[i]-(ReferencePoint[i] + offset))/stepSize);
00194       }
00195     }
00196     break;
00197   case 5:
00198     double sinPhi = sin(pnt[1]);
00199     double stepSize = RParAsFunOfPhi[0] + RParAsFunOfPhi[1]/sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3]/sinPhi;
00200     stepSize =  stepSize/(NumberOfPoints[0]-1);
00201     double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3]/sinPhi;
00202     index[0] = int((pnt[0]-startingPoint)/stepSize);
00203     index[1] = int((pnt[1]-ReferencePoint[1])/BasicDistance0[1]);
00204     index[2] = int((pnt[2]-ReferencePoint[2])/BasicDistance0[2]);
00205     break;
00206   }
00207   Index1 = index[0];
00208   Index2 = index[1];
00209   Index3 = index[2];
00210   return;
00211 }
00212 
00213 void MagneticFieldGrid::putIndicesGetB(int Index1, int Index2, int Index3, float &Bx, float &By, float &Bz){
00214   BVector FieldEntry;
00215   FieldEntry = FieldValues.operator[](lineNumber(Index1, Index2, Index3));
00216   Bx = FieldEntry.bx();
00217   By = FieldEntry.by();
00218   Bz = FieldEntry.bz();
00219   return;
00220 }
00221 
00222 // FIXME: same as above.
00223 void MagneticFieldGrid::putIndGetCoord(int Index1, int Index2, int Index3, double &X1, double &X2, double &X3){
00224   int index[3] = {Index1, Index2, Index3};
00225   double pnt[3];
00226   switch (GridType){
00227   case 1:
00228     for (int i=0; i<3; ++i){
00229       pnt[i] = ReferencePoint[i] + BasicDistance0[i]*index[i];
00230     }
00231     break;
00232   case 2:
00233     for (int i=0; i<3; ++i){
00234       if (EasyCoordinate[i]){
00235         pnt[i] = ReferencePoint[i] + BasicDistance0[i]*index[i];
00236       }
00237       else {
00238         double stepSize = BasicDistance0[i];
00239         double offset   = 0.0;
00240         for (int j=0; j<3; ++j){
00241           stepSize += BasicDistance1[i][j]*index[j];
00242           offset   += BasicDistance2[i][j]*index[j];
00243         }
00244         pnt[i] = ReferencePoint[i] + offset + stepSize*index[i];
00245       }
00246     }
00247     break;
00248   case 3:
00249     for (int i=0; i<3; ++i){
00250       pnt[i] = ReferencePoint[i] + BasicDistance0[i]*index[i];
00251     }
00252     break;
00253   case 4:
00254     for (int i=0; i<3; ++i){
00255       if (EasyCoordinate[i]){
00256         pnt[i] = ReferencePoint[i] + BasicDistance0[i]*index[i];
00257       }
00258       else {
00259         double stepSize = BasicDistance0[i];
00260         double offset   = 0.0;
00261         for (int j=0; j<3; ++j){
00262           stepSize += BasicDistance1[i][j]*index[j];
00263           offset   += BasicDistance2[i][j]*index[j];
00264         }
00265         pnt[i] = ReferencePoint[i] + offset + stepSize*index[i];
00266       }
00267     }
00268     break;
00269   case 5:
00270     pnt[2] = ReferencePoint[2] + BasicDistance0[2]*index[2];
00271     pnt[1] = ReferencePoint[1] + BasicDistance0[1]*index[1];
00272     double sinPhi = sin(pnt[1]);
00273     double stepSize = RParAsFunOfPhi[0] + RParAsFunOfPhi[1]/sinPhi - RParAsFunOfPhi[2] - RParAsFunOfPhi[3]/sinPhi;
00274     stepSize =  stepSize/(NumberOfPoints[0]-1);
00275     double startingPoint = RParAsFunOfPhi[2] + RParAsFunOfPhi[3]/sinPhi;
00276     pnt[0] = startingPoint + stepSize*index[0];
00277     break;
00278   }
00279   X1 = pnt[0];
00280   X2 = pnt[1];
00281   X3 = pnt[2];
00282   return;
00283 }
00284 
00285 int MagneticFieldGrid::lineNumber(int Index1, int Index2, int Index3){
00286   return Index1*NumberOfPoints[1]*NumberOfPoints[2] + Index2*NumberOfPoints[2] + Index3;
00287 }
00288 
00289 void MagneticFieldGrid::BVector::putB3(float Bx, float By, float Bz){
00290   B3[0] = Bx;
00291   B3[1] = By;
00292   B3[2] = Bz;
00293   return;
00294 }
00295 
00296 float MagneticFieldGrid::BVector::bx(){  return B3[0]; }
00297 
00298 float MagneticFieldGrid::BVector::by(){  return B3[1]; }
00299 
00300 float MagneticFieldGrid::BVector::bz(){  return B3[2]; }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/MagneticField/Interpolation/src/MagneticFieldGrid.cc
