Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-05-26 00:02:12

0001 /** \file HouseholderDecomposition.cc
0002  *
0003  *
0004  * \author Lorenzo Agostino, R.Ofierzynski, CERN
0005  */
0006 
0007 #include "Calibration/Tools/interface/HouseholderDecomposition.h"
0008 #include <cfloat>
0009 #include <cmath>
0010 #include <cstdlib>
0011 
0012 HouseholderDecomposition::HouseholderDecomposition(int squareMode_, int mineta_, int maxeta_, int minphi_, int maxphi_)
0013     : squareMode(squareMode_), countEvents(0), mineta(mineta_), maxeta(maxeta_), minphi(minphi_), maxphi(maxphi_) {
0014   Neta = maxeta - mineta + 1;
0015   if (mineta * maxeta < 0)
0016     Neta--;  // there's no eta index = 0
0017   Nphi = maxphi - minphi + 1;
0018   if (Nphi < 0)
0019     Nphi += 360;
0020 
0021   Nchannels = Neta * Nphi;  // no. of channels, get it from edges of the region
0022 
0023   Nxtals = squareMode * squareMode;  // no. of xtals in one event
0024 
0025   sigmaReplacement = 0.00001;  // the sum of columns is replaced by this value in case it is zero (e.g. dead crystal)
0026 }
0027 
0028 HouseholderDecomposition::~HouseholderDecomposition() {}
0029 
0030 std::vector<float> HouseholderDecomposition::runRegional(const std::vector<std::vector<float> >& eventMatrix,
0031                                                          const std::vector<int>& VmaxCeta,
0032                                                          const std::vector<int>& VmaxCphi,
0033                                                          const std::vector<float>& energyVector,
0034                                                          const int& nIter,
0035                                                          const int& regLength) {
0036   // make regions
0037   makeRegions(regLength);
0038 
0039   Nevents = eventMatrix.size();  // Number of events to calibrate with
0040 
0041   std::vector<float> totalSolution(Nchannels, 1.);
0042   std::vector<float> iterSolution(Nchannels, 1.);
0043   std::vector<std::vector<float> > myEventMatrix(eventMatrix);
0044 
0045   // loop over nIter
0046   for (int iter = 1; iter <= nIter; iter++) {
0047     // loop over regions
0048     for (unsigned int ireg = 0; ireg < regMinEta.size(); ireg++) {
0049       std::vector<float> regIterSolution, regEnergyVector;
0050       std::vector<int> regVmaxCeta, regVmaxCphi;
0051       std::vector<std::vector<float> > regEventMatrix;
0052 
0053       // initialize new instance with regional min,max indices
0054       HouseholderDecomposition regionalHH(
0055           squareMode, regMinEta[ireg], regMaxEta[ireg], regMinPhi[ireg], regMaxPhi[ireg]);
0056 
0057       // copy all events in region into new eventmatrix, energyvector, VmaxCeta, VmaxCphi
0058       for (unsigned int ia = 0; ia < VmaxCeta.size(); ia++) {
0059         if ((VmaxCeta[ia] >= regMinEta[ireg]) && (VmaxCeta[ia] <= regMaxEta[ireg]) &&
0060             (VmaxCphi[ia] >= regMinPhi[ireg]) && (VmaxCphi[ia] <= regMaxPhi[ireg])) {
0061           // save event, calculate new eventmatrix(truncated) and energy
0062           regVmaxCeta.push_back(VmaxCeta[ia]);
0063           regVmaxCphi.push_back(VmaxCphi[ia]);
0064 
0065           std::vector<float> regEvent = myEventMatrix[ia];
0066           float regEnergy = energyVector[ia];
0067           for (int i2 = 0; i2 < Nxtals; i2++) {
0068             int iFullReg = regionalHH.indexSqr2Reg(i2, VmaxCeta[ia], VmaxCphi[ia]);
0069             if (iFullReg < 0)  // crystal outside
0070             {
0071               regEnergy -= regEvent[i2];
0072               regEvent[i2] = 0.;
0073             }
0074           }
0075           regEventMatrix.push_back(regEvent);
0076           regEnergyVector.push_back(regEnergy);
0077         }
0078       }
0079 
0080       // calibrate
0081       //      std::cout << "HouseholderDecomposition::runRegional - Starting calibration of region " << ireg << ": eta "
0082       //           << regMinEta[ireg] << " to " << regMaxEta[ireg] << ", phi " << regMinPhi[ireg] << " to " << regMaxPhi[ireg] << std::endl;
0083       regIterSolution = regionalHH.iterate(regEventMatrix, regVmaxCeta, regVmaxCphi, regEnergyVector);
0084       //      std::cout << "HouseholderDecomposition::runRegional - calibration of region finished. " << std::endl;
0085 
0086       // save solution into global iterSolution
0087       // don't forget to delete the ones that are on the border !
0088       for (unsigned int i1 = 0; i1 < regIterSolution.size(); i1++) {
0089         int regFrame = regLength / 2;
0090         int currRegPhiRange = regMaxPhi[ireg] - regMinPhi[ireg] + 1;
0091         int currRegEta = i1 / currRegPhiRange + regMinEta[ireg];
0092         int currRegPhi = i1 % currRegPhiRange + regMinPhi[ireg];
0093         int newindex = -100;
0094         // if crystal well inside:
0095         if ((currRegEta >= (regMinEta[ireg] + regFrame * (!(regMinEta[ireg] == mineta)))) &&
0096             (currRegEta <= (regMaxEta[ireg] - regFrame * (!(regMaxEta[ireg] == maxeta)))) &&
0097             (currRegPhi >= (regMinPhi[ireg] + regFrame * (!(regMinPhi[ireg] == minphi)))) &&
0098             (currRegPhi <= (regMaxPhi[ireg] - regFrame * (!(regMaxPhi[ireg] == maxphi))))) {
0099           newindex = (currRegEta - mineta) * Nphi + currRegPhi - minphi;
0100           iterSolution[newindex] = regIterSolution[i1];
0101         }
0102       }
0103     }  // end loop over regions
0104 
0105     if (iterSolution.empty())
0106       return iterSolution;
0107 
0108     // re-calibrate eventMatrix with solution
0109     for (int ievent = 0; ievent < Nevents; ievent++) {
0110       myEventMatrix[ievent] = recalibrateEvent(myEventMatrix[ievent], VmaxCeta[ievent], VmaxCphi[ievent], iterSolution);
0111     }
0112 
0113     // save solution into theCalibVector
0114     for (int i = 0; i < Nchannels; i++) {
0115       totalSolution[i] *= iterSolution[i];
0116     }
0117 
0118   }  // end loop over nIter
0119 
0120   return totalSolution;
0121 }
0122 
0123 std::vector<float> HouseholderDecomposition::iterate(const std::vector<std::vector<float> >& eventMatrix,
0124                                                      const std::vector<int>& VmaxCeta,
0125                                                      const std::vector<int>& VmaxCphi,
0126                                                      const std::vector<float>& energyVector,
0127                                                      const int& nIter,
0128                                                      const bool& normalizeFlag) {
0129   Nevents = eventMatrix.size();  // Number of events to calibrate with
0130 
0131   std::vector<float> totalSolution(Nchannels, 1.);
0132   std::vector<float> iterSolution;
0133   std::vector<std::vector<float> > myEventMatrix(eventMatrix);
0134   std::vector<float> myEnergyVector(energyVector);
0135 
0136   int i, j;
0137 
0138   // Iterate the correction
0139   for (int iter = 1; iter <= nIter; iter++) {
0140     // if normalization flag is set, normalize energies
0141     float sumOverEnergy;
0142     if (normalizeFlag) {
0143       float scale = 0.;
0144 
0145       for (i = 0; i < Nevents; i++) {
0146         sumOverEnergy = 0.;
0147         for (j = 0; j < Nxtals; j++) {
0148           sumOverEnergy += myEventMatrix[i][j];
0149         }
0150         sumOverEnergy /= myEnergyVector[i];
0151         scale += sumOverEnergy;
0152       }
0153       scale /= Nevents;
0154 
0155       for (i = 0; i < Nevents; i++) {
0156         myEnergyVector[i] *= scale;
0157       }
0158     }  // end normalize energies
0159 
0160     // now the real work starts:
0161     iterSolution = iterate(myEventMatrix, VmaxCeta, VmaxCphi, myEnergyVector);
0162 
0163     if (iterSolution.empty())
0164       return iterSolution;
0165 
0166     // re-calibrate eventMatrix with solution
0167     for (int ievent = 0; ievent < Nevents; ievent++) {
0168       myEventMatrix[ievent] = recalibrateEvent(myEventMatrix[ievent], VmaxCeta[ievent], VmaxCphi[ievent], iterSolution);
0169     }
0170 
0171     for (int i = 0; i < Nchannels; i++) {
0172       // save solution into theCalibVector
0173       totalSolution[i] *= iterSolution[i];
0174     }
0175 
0176   }  // end iterate correction
0177 
0178   return totalSolution;
0179 }
0180 
0181 std::vector<float> HouseholderDecomposition::iterate(const std::vector<std::vector<float> >& eventMatrix,
0182                                                      const std::vector<int>& VmaxCeta,
0183                                                      const std::vector<int>& VmaxCphi,
0184                                                      const std::vector<float>& energyVectorOrig) {
0185   std::vector<float> solution;
0186 
0187   Nevents = eventMatrix.size();  // Number of events to calibrate with
0188 
0189   if (Nchannels > Nevents) {
0190     std::cout << "Householder::runIter(): more channels to calibrate than events available. " << std::endl;
0191     std::cout << "  Nchannels=" << Nchannels << std::endl;
0192     std::cout << "  Nevents=" << Nevents << std::endl;
0193     std::cout << " ******************    ERROR   *********************" << std::endl;
0194     return solution;  // empty vector
0195   }
0196 
0197   // input: eventMatrixOrig - the unzipped matrix
0198   eventMatrixOrig = unzipMatrix(eventMatrix, VmaxCeta, VmaxCphi);
0199 
0200   if (eventMatrixOrig.size() != energyVectorOrig.size()) {
0201     std::cout << "Householder::runIter(): matrix dimensions non-conformant. " << std::endl;
0202     std::cout << "  energyVectorOrig.size()=" << energyVectorOrig.size() << std::endl;
0203     std::cout << "  eventMatrixOrig.size()=" << eventMatrixOrig.size() << std::endl;
0204     std::cout << " ******************    ERROR   *********************" << std::endl;
0205     return solution;  // empty vector
0206   }
0207 
0208   int i, j;
0209   eventMatrixProc = eventMatrixOrig;
0210   energyVectorProc = energyVectorOrig;  // copy energyVectorOrig vector
0211   std::vector<float> e(Nchannels);
0212   alpha.assign(Nchannels, 0.);
0213   pivot.assign(Nchannels, 0);
0214 
0215   //--------------------
0216   bool decomposeSuccess = decompose();
0217 
0218   if (!decomposeSuccess) {
0219     std::cout << "Householder::runIter(): Failed: Singular condition in decomposition." << std::endl;
0220     std::cout << "***************** PROBLEM in DECOMPOSITION *************************" << std::endl;
0221     return solution;  // empty vector
0222   }
0223 
0224   /* DBL_EPSILON: Difference between 1.0 and the minimum float greater than 1.0 */
0225   float mydbleps = 2.22045e-16;  //DBL_EPSILON;
0226   float etasqr = mydbleps * mydbleps;
0227   //  std::cout << "LOOK at DBL_EPSILON :" << mydbleps <<std::endl;
0228 
0229   //--------------------
0230   // apply transformations to rhs - find solution vector
0231   solution.assign(Nchannels, 0.);
0232   solve(solution);
0233 
0234   // compute residual vector energyVectorProc
0235   for (i = 0; i < Nevents; i++) {
0236     energyVectorProc[i] = energyVectorOrig[i];
0237     for (j = 0; j < Nchannels; j++) {
0238       energyVectorProc[i] -= eventMatrixOrig[i][j] * solution[j];
0239     }
0240   }
0241 
0242   //--------------------
0243   // compute first correction vector e
0244   solve(e);
0245 
0246   float normy0 = 0.;
0247   float norme1 = 0.;
0248   float norme0;
0249 
0250   for (i = 0; i < Nchannels; i++) {
0251     normy0 += solution[i] * solution[i];
0252     norme1 += e[i] * e[i];
0253   }
0254 
0255   //  std::cout << "Householder::runIter(): applying first correction";
0256   //  std::cout << " normy0 = " << normy0;
0257   //  std::cout << " norme1 = " << norme1 << std::endl;
0258 
0259   // not attempt at obtaining the solution is made unless the norm of the first
0260   // correction  is significantly smaller than the norm of the initial solution
0261   if (norme1 > (0.0625 * normy0)) {
0262     //      std::cout << "Householder::runIter(): first correction is too large. Failed." << std::endl;
0263   }
0264 
0265   // improve the solution
0266   for (i = 0; i < Nchannels; i++) {
0267     solution[i] += e[i];
0268   }
0269 
0270   //  std::cout << "Householder::runIter(): improving solution...." << std::endl;
0271 
0272   //--------------------
0273   // only continue iteration if the correction was significant
0274   while (norme1 > (etasqr * normy0)) {
0275     //      std::cout << "Householder::runIter(): norme1 = " << norme1 << std::endl;
0276 
0277     for (i = 0; i < Nevents; i++) {
0278       energyVectorProc[i] = energyVectorOrig[i];
0279       for (j = 0; j < Nchannels; j++) {
0280         energyVectorProc[i] -= eventMatrixOrig[i][j] * solution[j];
0281       }
0282     }
0283 
0284     // compute next correction vector
0285     solve(e);
0286 
0287     norme0 = norme1;
0288     norme1 = 0.;
0289 
0290     for (i = 0; i < Nchannels; i++) {
0291       norme1 += e[i] * e[i];
0292     }
0293 
0294     // terminate iteration if the norm of the new correction failed to decrease
0295     // significantly compared to the norm of the previous correction
0296     if (norme1 > (0.0625 * norme0))
0297       break;
0298 
0299     // apply correction vector
0300     for (i = 0; i < Nchannels; i++) {
0301       solution[i] += e[i];
0302     }
0303   }
0304 
0305   //clean up
0306   eventMatrixOrig.clear();
0307   eventMatrixProc.clear();
0308   energyVectorProc.clear();
0309   alpha.clear();
0310   pivot.clear();
0311 
0312   return solution;
0313 }
0314 
0315 bool HouseholderDecomposition::decompose() {
0316   int i, j, jbar, k;
0317   float beta, sigma, alphak, eventMatrixkk;
0318   std::vector<float> y(Nchannels);
0319   std::vector<float> sum(Nchannels);
0320 
0321   //  std::cout << "Householder::decompose() started" << std::endl;
0322 
0323   for (j = 0; j < Nchannels; j++) {
0324     // jth column sum: squared sum for each crystal
0325     sum[j] = 0.;
0326     for (i = 0; i < Nevents; i++)
0327       sum[j] += eventMatrixProc[i][j] * eventMatrixProc[i][j];
0328 
0329     // bookkeeping vector
0330     pivot[j] = j;
0331   }
0332 
0333   for (k = 0; k < Nchannels; k++) {
0334     // kth Householder transformation
0335     sigma = sum[k];
0336     jbar = k;
0337 
0338     // go through all following columns
0339     // find the largest sumSquared in the following columns
0340     for (j = k + 1; j < Nchannels; j++) {
0341       if (sum[j] > sigma) {
0342         sigma = sum[j];
0343         jbar = j;
0344       }
0345     }
0346 
0347     if (jbar != k) {
0348       // column interchange:
0349       //     interchange within: bookkeeping vector, squaredSum, eventMatrixProc
0350 
0351       i = pivot[k];
0352       pivot[k] = pivot[jbar];
0353       pivot[jbar] = i;
0354 
0355       sum[jbar] = sum[k];
0356       sum[k] = sigma;
0357 
0358       for (i = 0; i < Nevents; i++) {
0359         sigma = eventMatrixProc[i][k];
0360         eventMatrixProc[i][k] = eventMatrixProc[i][jbar];
0361         eventMatrixProc[i][jbar] = sigma;
0362       }
0363     }  // end column interchange
0364 
0365     // now store in sigma the squared sum of the readoutEnergies for this column(crystal)
0366     sigma = 0.;
0367     for (i = k; i < Nevents; i++) {
0368       sigma += eventMatrixProc[i][k] * eventMatrixProc[i][k];
0369     }
0370 
0371     // found a zero-column, bail out
0372     if (sigma == 0.) {
0373       //      std::cout << "Householder::decompose() failed" << std::endl;
0374       //      return false;
0375       // rof 14.12.2006: workaround to avoid failure of algorithm because of dead crystals:
0376       sigma = sigmaReplacement;
0377       //      std::cout << "Householder::decompose - found zero column " << jbar << ", replacing sum of column elements by " << sigma << std::endl;
0378     }
0379 
0380     // the following paragraph acts only on the diagonal element:
0381     // if element=0, then calculate alpha & beta
0382 
0383     // take the diagonal element
0384     eventMatrixkk = eventMatrixProc[k][k];
0385 
0386     if (eventMatrixkk < 0.)
0387       alpha[k] = sqrt(sigma);
0388     else
0389       alpha[k] = sqrt(sigma) * (-1.);
0390 
0391     alphak = alpha[k];
0392 
0393     beta = 1 / (sigma - eventMatrixkk * alphak);
0394     // replace it
0395     eventMatrixProc[k][k] = eventMatrixkk - alphak;
0396 
0397     for (j = k + 1; j < Nchannels; j++) {
0398       y[j] = 0.;
0399 
0400       for (i = k; i < Nevents; i++) {
0401         y[j] += eventMatrixProc[i][k] * eventMatrixProc[i][j];
0402       }
0403 
0404       y[j] *= beta;
0405     }
0406 
0407     for (j = k + 1; j < Nchannels; j++) {
0408       for (i = k; i < Nevents; i++) {
0409         eventMatrixProc[i][j] -= eventMatrixProc[i][k] * y[j];
0410         sum[j] -= eventMatrixProc[k][j] * eventMatrixProc[k][j];
0411       }
0412     }
0413   }  // end of kth householder transformation
0414 
0415   //  std::cout << "Householder::decompose() finished" << std::endl;
0416 
0417   return true;
0418 }
0419 
0420 void HouseholderDecomposition::solve(std::vector<float>& y) {
0421   std::vector<float> z(Nchannels, 0.);
0422 
0423   float gamma;
0424   int i, j;
0425 
0426   //  std::cout << "Householder::solve() begin" << std::endl;
0427 
0428   for (j = 0; j < Nchannels; j++) {
0429     // apply jth transformation to the right hand side
0430     gamma = 0.;
0431     for (i = j; i < Nevents; i++) {
0432       gamma += eventMatrixProc[i][j] * energyVectorProc[i];
0433     }
0434     gamma /= (alpha[j] * eventMatrixProc[j][j]);
0435 
0436     for (i = j; i < Nevents; i++) {
0437       energyVectorProc[i] += gamma * eventMatrixProc[i][j];
0438     }
0439   }
0440 
0441   z[Nchannels - 1] = energyVectorProc[Nchannels - 1] / alpha[Nchannels - 1];
0442 
0443   for (i = Nchannels - 2; i >= 0; i--) {
0444     z[i] = energyVectorProc[i];
0445     for (j = i + 1; j < Nchannels; j++) {
0446       z[i] -= eventMatrixProc[i][j] * z[j];
0447     }
0448     z[i] /= alpha[i];
0449   }
0450 
0451   for (i = 0; i < Nchannels; i++) {
0452     y[pivot[i]] = z[i];
0453   }
0454 
0455   //  std::cout << "Householder::solve() finished." << std::endl;
0456 }
0457 
0458 std::vector<float> HouseholderDecomposition::recalibrateEvent(const std::vector<float>& eventSquare,
0459                                                               const int& maxCeta,
0460                                                               const int& maxCphi,
0461                                                               const std::vector<float>& recalibrateVector) {
0462   std::vector<float> newEventSquare(eventSquare);
0463   int iFull;
0464 
0465   for (int i = 0; i < Nxtals; i++) {
0466     iFull = indexSqr2Reg(i, maxCeta, maxCphi);
0467     if (iFull >= 0)
0468       newEventSquare[i] *= recalibrateVector[iFull];
0469   }
0470   return newEventSquare;
0471 }
0472 
0473 int HouseholderDecomposition::indexSqr2Reg(const int& sqrIndex, const int& maxCeta, const int& maxCphi) {
0474   int regionIndex;
0475 
0476   // get the current eta, phi indices
0477   int curr_eta = maxCeta - squareMode / 2 + sqrIndex % squareMode;
0478   if (curr_eta * maxCeta <= 0) {
0479     if (maxCeta > 0)
0480       curr_eta--;
0481     else
0482       curr_eta++;
0483   }  // JUMP over 0
0484 
0485   int curr_phi = maxCphi - squareMode / 2 + sqrIndex / squareMode;
0486   if (curr_phi < 1)
0487     curr_phi += 360;
0488   if (curr_phi > 360)
0489     curr_phi -= 360;
0490 
0491   bool negPhiDirection = (maxphi < minphi);
0492   int iFullphi;
0493 
0494   regionIndex = -1;
0495 
0496   if (curr_eta >= mineta && curr_eta <= maxeta)
0497     if ((!negPhiDirection && (curr_phi >= minphi && curr_phi <= maxphi)) ||
0498         (negPhiDirection && !(curr_phi >= minphi && curr_phi <= maxphi))) {
0499       iFullphi = curr_phi - minphi;
0500       if (iFullphi < 0)
0501         iFullphi += 360;
0502       regionIndex = (curr_eta - mineta) * (maxphi - minphi + 1 + 360 * negPhiDirection) + iFullphi;
0503     }
0504 
0505   return regionIndex;
0506 }
0507 
0508 std::vector<std::vector<float> > HouseholderDecomposition::unzipMatrix(
0509     const std::vector<std::vector<float> >& eventMatrix,
0510     const std::vector<int>& VmaxCeta,
0511     const std::vector<int>& VmaxCphi) {
0512   std::vector<std::vector<float> > fullMatrix;
0513 
0514   int iFull;
0515 
0516   for (int i = 0; i < Nevents; i++) {
0517     std::vector<float> foo(Nchannels, 0.);
0518     for (int k = 0; k < Nxtals; k++) {
0519       iFull = indexSqr2Reg(k, VmaxCeta[i], VmaxCphi[i]);
0520       if (iFull >= 0)
0521         foo[iFull] = eventMatrix[i][k];
0522     }
0523     fullMatrix.push_back(foo);
0524   }
0525 
0526   return fullMatrix;
0527 }
0528 
0529 void HouseholderDecomposition::makeRegions(const int& regLength) {
0530   //  int regFrame = regLength/2;
0531   int regFrame = squareMode / 2;
0532   // first eta:
0533   int remEta = Neta % regLength;
0534   if (remEta > regLength / 2)
0535     remEta -= regLength;
0536   int numSubRegEta = Neta / regLength + (remEta < 0) * 1;
0537 
0538   int addtoEta = remEta / numSubRegEta;
0539   int addtomoreEta =
0540       remEta % numSubRegEta;  // add "addtomore" number of times (addto+1), remaining times add just (addto)
0541 
0542   // then phi:
0543   int remPhi = Nphi % regLength;
0544   if (remPhi > regLength / 2)
0545     remPhi -= regLength;
0546   int numSubRegPhi = Nphi / regLength + (remPhi < 0) * 1;
0547 
0548   int addtoPhi = remPhi / numSubRegPhi;
0549   int addtomorePhi =
0550       remPhi % numSubRegPhi;  // add "addtomore" number of times (addto+-1), remaining times add just (addto)
0551 
0552   // now put it all together
0553   int startIndexEta = mineta;
0554   int startIndexPhi;
0555   int endIndexEta;
0556   int endIndexPhi;
0557   for (int i = 0; i < numSubRegEta; i++) {
0558     int addedLengthEta = regLength + addtoEta + addtomoreEta / abs(addtomoreEta) * (i < abs(addtomoreEta));
0559     endIndexEta = startIndexEta + addedLengthEta - 1;
0560 
0561     startIndexPhi = minphi;
0562     for (int j = 0; j < numSubRegPhi; j++) {
0563       int addedLengthPhi = regLength + addtoPhi + addtomorePhi / abs(addtomorePhi) * (j < abs(addtomorePhi));
0564       endIndexPhi = startIndexPhi + addedLengthPhi - 1;
0565 
0566       regMinPhi.push_back(startIndexPhi - regFrame * (j != 0));
0567       regMaxPhi.push_back(endIndexPhi + regFrame * (j != (numSubRegPhi - 1)));
0568       regMinEta.push_back(startIndexEta - regFrame * (i != 0));
0569       regMaxEta.push_back(endIndexEta + regFrame * (i != (numSubRegEta - 1)));
0570 
0571       startIndexPhi = endIndexPhi + 1;
0572     }
0573     startIndexEta = endIndexEta + 1;
0574   }
0575 
0576   //  // print it all
0577   //  std::cout << "Householder::makeRegions created the following regions for calibration:" << std::endl;
0578   //  for (int i=0; i<regMinEta.size(); i++)
0579   //    std::cout << "Region " << i << ": eta = " << regMinEta[i] << " to " << regMaxEta[i] << ", phi = " << regMinPhi[i] << " to " << regMaxPhi[i] << std::endl;
0580 }