Interpolation/src/SpecialCylindricalMFGrid.cc

0001 #include "SpecialCylindricalMFGrid.h"
0002 #include "MagneticField/Interpolation/interface/binary_ifstream.h"
0003 #include "LinearGridInterpolator3D.h"
0004
0005 #include <iostream>
0006
0007 using namespace std;
0008
0009 SpecialCylindricalMFGrid::SpecialCylindricalMFGrid(binary_ifstream& inFile,
0010                                                    const GloballyPositioned<float>& vol,
0011                                                    int gridType)
0012     : MFGrid3D(vol) {
0013   if (gridType == 5) {
0014     sector1 = false;
0015   } else if (gridType == 6) {
0016     sector1 = true;
0017   } else {
0018     cout << "ERROR wrong SpecialCylindricalMFGrid type " << gridType << endl;
0019     sector1 = false;
0020   }
0021
0022   int n1, n2, n3;
0023   inFile >> n1 >> n2 >> n3;
0024 #ifdef DEBUG_GRID
0025   cout << "n1 " << n1 << " n2 " << n2 << " n3 " << n3 << endl;
0026 #endif
0027   double xref, yref, zref;
0028   inFile >> xref >> yref >> zref;
0029   double stepx, stepy, stepz;
0030   inFile >> stepx >> stepy >> stepz;
0031
0032   double RParAsFunOfPhi[4];  // R = f(phi) or const. (0,2: const. par. ; 1,3: const./sin(phi));
0033   inFile >> RParAsFunOfPhi[0] >> RParAsFunOfPhi[1] >> RParAsFunOfPhi[2] >> RParAsFunOfPhi[3];
0034
0035   vector<BVector> fieldValues;
0036   float Bx, By, Bz;
0037   int nLines = n1 * n2 * n3;
0038   fieldValues.reserve(nLines);
0039   for (int iLine = 0; iLine < nLines; ++iLine) {
0040     inFile >> Bx >> By >> Bz;
0041     // This would be fine only if local r.f. has the axes oriented as the global r.f.
0042     // For this volume we know that the local and global r.f. have different axis
0043     // orientation, so we do not try to be clever.
0044     //    fieldValues.push_back(BVector(Bx,By,Bz));
0045
0046     // Preserve double precision!
0047     Vector3DBase<double, LocalTag> lB = frame().toLocal(Vector3DBase<double, GlobalTag>(Bx, By, Bz));
0048     fieldValues.push_back(BVector(lB.x(), lB.y(), lB.z()));
0049   }
0050   // check completeness
0051   string lastEntry;
0052   inFile >> lastEntry;
0053   if (lastEntry != "complete") {
0054     cout << "ERROR during file reading: file is not complete" << endl;
0055   }
0056
0057   GlobalPoint grefp(GlobalPoint::Cylindrical(xref, yref, zref));
0058
0059 #ifdef DEBUG_GRID
0060   LocalPoint lrefp = frame().toLocal(grefp);
0061   cout << "Grid reference point in grid system: " << xref << "," << yref << "," << zref << endl;
0062   cout << "Grid reference point in global x,y,z: " << grefp << endl;
0063   cout << "Grid reference point in local x,y,z: " << lrefp << endl;
0064   cout << "steps " << stepx << "," << stepy << "," << stepz << endl;
0065   cout << "RParAsFunOfPhi[0...4] = ";
0066   for (int i = 0; i < 4; ++i)
0067     cout << RParAsFunOfPhi[i] << " ";
0068   cout << endl;
0069 #endif
0070
0071   Grid1D gridX(0, n1 - 1, n1);  // offset and step size not constant
0072   Grid1D gridY(yref, yref + stepy * (n2 - 1), n2);
0073   Grid1D gridZ(grefp.z(), grefp.z() + stepz * (n3 - 1), n3);
0074
0075   grid_ = GridType(gridX, gridY, gridZ, fieldValues);
0076
0077   stepConstTerm_ = (RParAsFunOfPhi[0] - RParAsFunOfPhi[2]) / (n1 - 1);
0078   stepPhiTerm_ = (RParAsFunOfPhi[1] - RParAsFunOfPhi[3]) / (n1 - 1);
0079   startConstTerm_ = RParAsFunOfPhi[2];
0080   startPhiTerm_ = RParAsFunOfPhi[3];
0081
0082   // Activate/deactivate timers
0083   //   static SimpleConfigurable<bool> timerOn(false,"MFGrid:timing");
0084   //   (*TimingReport::current()).switchOn("MagneticFieldProvider::valueInTesla(SpecialCylindricalMFGrid)",timerOn);
0085 }
0086
0087 MFGrid::LocalVector SpecialCylindricalMFGrid::uncheckedValueInTesla(const LocalPoint& p) const {
0088   //   static TimingReport::Item & timer= (*TimingReport::current())["MagneticFieldProvider::valueInTesla(SpecialCylindricalMFGrid)"];
0089   //   TimeMe t(timer,false);
0090
0091   LinearGridInterpolator3D interpol(grid_);
0092   double a, b, c;
0093   toGridFrame(p, a, b, c);
0094   // the following holds if B values was not converted to local coords -- see ctor
0095   //   GlobalVector gv( interpol.interpolate( a, b, c)); // grid in global frame
0096   //   return frame().toLocal(gv);           // must return a local vector
0097   return LocalVector(interpol.interpolate(a, b, c));
0098 }
0099
0100 void SpecialCylindricalMFGrid::dump() const {}
0101
0102 MFGrid::LocalPoint SpecialCylindricalMFGrid::fromGridFrame(double a, double b, double c) const {
0103   double sinPhi;  // sin or cos depending on wether we are at phi=0 or phi=pi/2
0104   if (sector1) {
0105     sinPhi = cos(b);
0106   } else {
0107     sinPhi = sin(b);
0108   }
0109
0110   double R = a * stepSize(sinPhi) + startingPoint(sinPhi);
0111   // "OLD" convention of phi.
0112   //  GlobalPoint gp( GlobalPoint::Cylindrical(R, Geom::pi() - b, c));
0113   GlobalPoint gp(GlobalPoint::Cylindrical(R, b, c));
0114   return frame().toLocal(gp);
0115 }
0116
0117 void SpecialCylindricalMFGrid::toGridFrame(const LocalPoint& p, double& a, double& b, double& c) const {
0118   GlobalPoint gp = frame().toGlobal(p);
0119   double sinPhi;  // sin or cos depending on wether we are at phi=0 or phi=pi/2
0120   if (sector1) {
0121     sinPhi = cos(gp.phi());
0122   } else {
0123     sinPhi = sin(gp.phi());
0124   }
0125   a = (gp.perp() - startingPoint(sinPhi)) / stepSize(sinPhi);
0126   // FIXME: "OLD" convention of phi.
0127   // b = Geom::pi() - gp.phi();
0128   b = gp.phi();
0129   c = gp.z();
0130
0131 #ifdef DEBUG_GRID
0132   if (sector1) {
0133     cout << "toGridFrame: sinPhi ";
0134   } else {
0135     cout << "toGridFrame: cosPhi ";
0136   }
0137   cout << sinPhi << " LocalPoint " << p << " GlobalPoint " << gp << endl
0138        << " a " << a << " b " << b << " c " << c << endl;
0139
0140 #endif
0141 }