CocoaModel/src/DeviationsFromFileSensor2D.cc

0001 //   COCOA class implementation file
0002 //Id: DeviationsFromFileSensor2D.cc
0003 //CAT: Model
0004 //
0005 //   History: v1.0
0006 //   Pedro Arce
0007
0008 #include "Alignment/CocoaModel/interface/DeviationsFromFileSensor2D.h"
0009 #include "Alignment/CocoaUtilities/interface/ALIFileIn.h"
0010 #include "Alignment/CocoaUtilities/interface/ALIUtils.h"
0011 #include <cassert>
0012 #include <cstdlib>
0013 #include <cmath>  // include floating-point std::abs functions
0014 #include <memory>
0015
0016 enum directions { xdir = 0, ydir = 1 };
0017
0018 ALIbool DeviationsFromFileSensor2D::theApply = true;
0019
0020 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0021 void DeviationsFromFileSensor2D::readFile(ALIFileIn& ifdevi) {
0022   verbose = ALIUtils::debug;
0023   //  verbose = 4;
0024   if (verbose >= 3)
0025     std::cout << "DeviationsFromFileSensor2D::readFile " << this << ifdevi.name() << std::endl;
0026
0027   theScanSenseX = 0;
0028   theScanSenseY = 0;
0029
0030   ALIuint nl = 1;
0031
0032   ALIdouble oldposX = 0, oldposY = 0;
0033   vd vcol;
0034   std::vector<ALIstring> wl;
0035   /*  //------ first line with dimension factors //now with global option
0036   ifdevi.getWordsInLine( wl );
0037   if( wl[0] != "DIMFACTOR:" || wl.size() != 3) {
0038     std::cerr << "Error reading sensor2D deviation file " << ifdevi.name() << std::endl
0039      << " first line has to be DIMFACTOR: 'posDimFactor' 'angDimFactor' " << std::endl;
0040     ALIUtils::dumpVS( wl, " ");
0041     exit(1);
0042   }
0043   ALIdouble posDimFactor = atof(wl[1].c_str());
0044   ALIdouble angDimFactor = atof(wl[2].c_str());
0045   */
0046
0047   for (;;) {
0048     ifdevi.getWordsInLine(wl);
0049     if (ifdevi.eof())
0050       break;
0051
0052     DeviationSensor2D* dev = new DeviationSensor2D();
0053     dev->fillData(wl);
0054
0055     if (verbose >= 5) {
0056       ALIUtils::dumpVS(wl, "deviation sensor2D", std::cout);
0057     }
0058
0059     //--- line 2 of data
0060     if (nl == 2) {
0061       //--------- get if scan is first in Y or X
0062       firstScanDir = ydir;
0063       if (verbose >= 3)
0064         std::cout << "firstScanDir " << firstScanDir << " " << dev->posX() << " " << oldposX << " " << dev->posY()
0065                   << " " << oldposY << std::endl;
0066       if (std::abs(dev->posX() - oldposX) > std::abs(dev->posY() - oldposY)) {
0067         std::cerr << "!!!EXITING first scan direction has to be Y for the moment " << ifdevi.name() << std::endl;
0068         firstScanDir = xdir;
0069         exit(1);
0070       }
0071       //-------- get sense of first scan direction
0072       if (firstScanDir == ydir) {
0073         if (dev->posY() > oldposY) {
0074           theScanSenseY = +1;
0075         } else {
0076           theScanSenseY = -1;
0077         }
0078         if (verbose >= 3)
0079           std::cout << " theScanSenseY " << theScanSenseY << std::endl;
0080       } else {
0081         if (dev->posX() > oldposX) {
0082           theScanSenseX = +1;
0083         } else {
0084           theScanSenseX = -1;
0085         }
0086         if (verbose >= 3)
0087           std::cout << " theScanSenseX " << theScanSenseX << std::endl;
0088       }
0089     }
0090
0091     //-      std::cout << "firstScanDir " << firstScanDir << " " <<  dev->posX() <<  " " << oldposX << " " <<  dev->posY() <<  " " << oldposY << std::endl;
0092
0093     //------- check if change of row: clear current std::vector of column values
0094     if ((firstScanDir == ydir && (dev->posY() - oldposY) * theScanSenseY < 0) ||
0095         (firstScanDir == xdir && (dev->posX() - oldposX) * theScanSenseX < 0)) {
0096       theDeviations.push_back(vcol);
0097       vcol.clear();
0098
0099       //-      std::cout << " theDevi size " << theDeviations.size() << " " << ifdevi.name()  << std::endl;
0100       //-------- get sense of second scan direction
0101       if (theScanSenseY == 0) {
0102         if (dev->posY() > oldposY) {
0103           theScanSenseY = +1;
0104         } else {
0105           theScanSenseY = -1;
0106         }
0107       }
0108       if (theScanSenseX == 0) {
0109         if (dev->posX() > oldposX) {
0110           theScanSenseX = +1;
0111         } else {
0112           theScanSenseX = -1;
0113         }
0114       }
0115       if (verbose >= 3)
0116         std::cout << " theScanSenseX " << theScanSenseX << " theScanSenseY " << theScanSenseY << std::endl;
0117     }
0118
0119     oldposX = dev->posX();
0120     oldposY = dev->posY();
0121
0122     //--- fill deviation std::vectors
0123     vcol.push_back(dev);
0124     nl++;
0125   }
0126   theDeviations.push_back(vcol);
0127
0128   //----- Calculate std::vector size
0129   ALIdouble dvsiz = (sqrt(ALIdouble(nl - 1)));
0130   //----- Check that it is a square of points (<=> vsiz is integer)
0131   ALIuint vsiz = ALIuint(dvsiz);
0132   if (vsiz != dvsiz) {
0133     if (ALIUtils::debug >= 0)
0134       std::cerr << "!!WARNING: error reading deviation file: number of X points <> number of Y points : Number of "
0135                    "points in X "
0136                 << dvsiz << " nl " << nl - 1 << " file " << ifdevi.name() << std::endl;
0137     //    exit(1);
0138   }
0139   theNPoints = vsiz;
0140
0141   if (verbose >= 4) {
0142     std::cout << " Filling deviation from file: " << ifdevi.name() << " theNPoints " << theNPoints << std::endl;
0143   }
0144 }
0145
0146 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0147 std::pair<ALIdouble, ALIdouble> DeviationsFromFileSensor2D::getDevis(ALIdouble intersX, ALIdouble intersY) {
0148   //intersX += 10000;
0149   //intersY += 10000;
0150   if (verbose >= 4)
0151     std::cout << " entering getdevis " << intersX << " " << intersY << " " << this << std::endl;
0152   vvd::iterator vvdite;
0153   vd::iterator vdite;
0154
0155   // assume scan is in Y first, else revers intersection coordinates
0156   /*  if( !firstScanDirY ) {
0157     ALIdouble tt = intersY;
0158     intersY = intersX;
0159     intersX = tt;
0160   }
0161 */
0162
0163   //---------- look which point in the deviation matrices correspond to intersX,Y
0164   //----- Look for each column, between which rows intersY is
0165   //assume first dir is Y
0166   auto yrows = std::make_unique<unsigned int[]>(theNPoints);
0167
0168   unsigned int ii = 0;
0169   ALIbool insideMatrix = false;
0170   for (vvdite = theDeviations.begin(); vvdite != (theDeviations.end() - 1); ++vvdite) {
0171     for (vdite = (*vvdite).begin(); vdite != ((*vvdite).end() - 1); ++vdite) {
0172       if (verbose >= 5)
0173         std::cout << " check posy " << (*(vdite))->posY() << " " << (*(vdite + 1))->posY() << " " << (*(vdite))
0174                   << std::endl;
0175       // if posy is between this point and previous one
0176
0177       //-     std::cout << "intersy" << intersY << " " <<  (*(vdite))->posY() << std::endl;
0178
0179       if ((intersY - (*(vdite))->posY()) * theScanSenseY > 0 &&
0180           (intersY - (*(vdite + 1))->posY()) * theScanSenseY < 0) {
0181         //-std::cout << " ii " << ii << std::endl;
0182         yrows[ii] = vdite - (*vvdite).begin();
0183         if (verbose >= 3)
0184           std::cout << intersY << " DeviationsFromFileSensor2D yrows " << ii << " " << yrows[ii] << " : "
0185                     << (*(vdite))->posY() << std::endl;
0186         insideMatrix = true;
0187         break;
0188       }
0189     }
0190     ii++;
0191   }
0192   assert(ii == theNPoints);
0193   if (insideMatrix == 0) {
0194     std::cerr << "!!EXITING intersection in Y outside matrix of deviations from file " << intersY << std::endl;
0195     exit(1);
0196   }
0197   insideMatrix = false;
0198
0199   vd thePoints;
0200   thePoints.clear();
0201   //----- For each row in 'yrows' look between which columns intersX is
0202   unsigned int rn;
0203   DeviationSensor2D *dev1, *dev2;
0204   for (ii = 0; ii < theNPoints - 1; ii++) {
0205     rn = yrows[ii];
0206     //-    std::cout << ii << " rn " << rn << std::endl;
0207     dev1 = (*(theDeviations.begin() + ii))[rn];  // column ii, row yrows[ii]
0208     rn = yrows[ii + 1];
0209     dev2 = (*(theDeviations.begin() + ii + 1))[rn];  // column ii+1, row yrows[ii+1]
0210     if ((intersX - dev1->posX()) * theScanSenseX > 0 && (intersX - dev2->posX()) * theScanSenseX < 0) {
0211       thePoints.push_back(dev1);
0212       thePoints.push_back(dev2);
0213       insideMatrix = true;
0214       if (verbose >= 3)
0215         std::cout << " column up " << ii << " " << dev1->posX() << " " << dev2->posX() << " : " << intersX << std::endl;
0216     }
0217
0218     rn = yrows[ii] + 1;
0219     if (rn == theNPoints)
0220       rn = theNPoints - 1;
0221     dev1 = (*(theDeviations.begin() + ii))[rn];  // column ii, row yrows[ii]+1
0222     rn = yrows[ii + 1] + 1;
0223     if (rn == theNPoints)
0224       rn = theNPoints - 1;
0225     dev2 = (*(theDeviations.begin() + ii + 1))[rn];  // column ii+1, row yrows[ii+1]+1
0226     if ((intersX - dev1->posX()) * theScanSenseX > 0 && (intersX - dev2->posX()) * theScanSenseX < 0) {
0227       thePoints.push_back(dev1);
0228       thePoints.push_back(dev2);
0229       if (verbose >= 3)
0230         std::cout << " column down " << ii << " " << dev1->posX() << " " << dev2->posX() << " : " << intersX
0231                   << std::endl;
0232     }
0233
0234     if (thePoints.size() == 4)
0235       break;
0236   }
0237
0238   if (insideMatrix == 0) {
0239     std::cerr << "!!EXITING intersection in X outside matrix of deviations from file " << intersX << std::endl;
0240     exit(1);
0241   }
0242
0243   //----------- If loop finished and not 4 points, point is outside scan bounds
0244
0245   //----- calculate deviation in x and y interpolating between four points
0246   ALIdouble dist, disttot = 0, deviX = 0, deviY = 0;
0247
0248   if (verbose >= 4)
0249     std::cout << " thepoints size " << thePoints.size() << std::endl;
0250
0251   for (ii = 0; ii < 4; ii++) {
0252     dist = sqrt(pow(thePoints[ii]->posX() - intersX, 2) + pow(thePoints[ii]->posY() - intersY, 2));
0253     disttot += 1. / dist;
0254     deviX += thePoints[ii]->devX() / dist;
0255     deviY += thePoints[ii]->devY() / dist;
0256     if (verbose >= 4) {
0257       //t      std::cout << ii << " point " << *thePoints[ii] << std::endl;
0258       std::cout << ii << " distances: " << dist << " " << deviX << " " << deviY << " devX " << thePoints[ii]->devX()
0259                 << " devY " << thePoints[ii]->devY() << std::endl;
0260     }
0261   }
0262   deviX /= disttot;
0263   deviY /= disttot;
0264
0265   // add offset
0266   deviX += theOffsetX;
0267   deviY += theOffsetY;
0268   if (verbose >= 4) {
0269     std::cout << " devisX/Y: " << deviX << " " << deviY << " intersX/Y " << intersX << " " << intersY << std::endl;
0270   }
0271
0272   //change sign!!!!!?!?!?!?!?!
0273   deviX *= -1;
0274   deviY *= -1;
0275   return std::pair<ALIdouble, ALIdouble>(deviX * 1.e-6, deviY * 1e-6);  // matrix is in microrad
0276 }