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