CMS 3D CMS Logo

HouseholderDecomposition.cc
Go to the documentation of this file.
1 
8 #include <cfloat>
9 #include <cmath>
10 #include <cstdlib>
11 
12 HouseholderDecomposition::HouseholderDecomposition(int squareMode_, int mineta_, int maxeta_, int minphi_, int maxphi_)
13  :squareMode(squareMode_), countEvents(0),mineta(mineta_), maxeta(maxeta_), minphi(minphi_), maxphi(maxphi_)
14 {
15  Neta = maxeta - mineta + 1;
16  if (mineta * maxeta < 0) Neta--; // there's no eta index = 0
17  Nphi = maxphi - minphi + 1;
18  if (Nphi <0) Nphi += 360;
19 
20  Nchannels = Neta * Nphi; // no. of channels, get it from edges of the region
21 
22  Nxtals = squareMode*squareMode; // no. of xtals in one event
23 
24  sigmaReplacement = 0.00001; // the sum of columns is replaced by this value in case it is zero (e.g. dead crystal)
25 
26 }
27 
29 {
30 }
31 
32 
33 std::vector<float> HouseholderDecomposition::runRegional(const std::vector<std::vector<float> >& eventMatrix, const std::vector<int>& VmaxCeta, const std::vector<int>& VmaxCphi, const std::vector<float>& energyVector, const int& nIter, const int& regLength)
34 {
35  // make regions
36  makeRegions(regLength);
37 
38  Nevents = eventMatrix.size(); // Number of events to calibrate with
39 
40  std::vector<float> totalSolution(Nchannels,1.);
41  std::vector<float> iterSolution(Nchannels,1.);
42  std::vector<std::vector<float> > myEventMatrix(eventMatrix);
43  std::vector<float> myEnergyVector(energyVector);
44 
45  // loop over nIter
46  for (int iter=1;iter<=nIter;iter++)
47  {
48  // loop over regions
49  for (unsigned int ireg=0; ireg<regMinEta.size(); ireg++)
50  {
51  std::vector<float> regIterSolution, regEnergyVector;
52  std::vector<int> regVmaxCeta, regVmaxCphi;
53  std::vector<std::vector<float> > regEventMatrix;
54 
55  // initialize new instance with regional min,max indices
56  HouseholderDecomposition regionalHH(squareMode,regMinEta[ireg],regMaxEta[ireg],regMinPhi[ireg],regMaxPhi[ireg]);
57 
58  // copy all events in region into new eventmatrix, energyvector, VmaxCeta, VmaxCphi
59  for (unsigned int ia=0; ia<VmaxCeta.size(); ia++)
60  {
61  if ((VmaxCeta[ia] >= regMinEta[ireg]) && (VmaxCeta[ia] <= regMaxEta[ireg])
62  && (VmaxCphi[ia] >= regMinPhi[ireg]) && (VmaxCphi[ia] <= regMaxPhi[ireg]))
63  {
64  // save event, calculate new eventmatrix(truncated) and energy
65  regVmaxCeta.push_back(VmaxCeta[ia]);
66  regVmaxCphi.push_back(VmaxCphi[ia]);
67 
68  std::vector<float> regEvent = myEventMatrix[ia];
69  float regEnergy = energyVector[ia];
70  for (int i2=0; i2<Nxtals; i2++)
71  {
72  int iFullReg = regionalHH.indexSqr2Reg(i2, VmaxCeta[ia], VmaxCphi[ia]);
73  if (iFullReg <0) // crystal outside
74  {
75  regEnergy -= regEvent[i2];
76  regEvent[i2] = 0.;
77  }
78  }
79  regEventMatrix.push_back(regEvent);
80  regEnergyVector.push_back(regEnergy);
81  }
82  }
83 
84  // calibrate
85  // std::cout << "HouseholderDecomposition::runRegional - Starting calibration of region " << ireg << ": eta "
86  // << regMinEta[ireg] << " to " << regMaxEta[ireg] << ", phi " << regMinPhi[ireg] << " to " << regMaxPhi[ireg] << std::endl;
87  regIterSolution = regionalHH.iterate(regEventMatrix, regVmaxCeta, regVmaxCphi, regEnergyVector);
88  // std::cout << "HouseholderDecomposition::runRegional - calibration of region finished. " << std::endl;
89 
90  // save solution into global iterSolution
91  // don't forget to delete the ones that are on the border !
92  for (unsigned int i1=0; i1<regIterSolution.size(); i1++)
93  {
94  int regFrame = regLength/2;
95  int currRegPhiRange = regMaxPhi[ireg] - regMinPhi[ireg] + 1;
96  int currRegEta = i1 / currRegPhiRange + regMinEta[ireg];
97  int currRegPhi = i1 % currRegPhiRange + regMinPhi[ireg];
98  int newindex = -100;
99  // if crystal well inside:
100  if ( (currRegEta >= (regMinEta[ireg]+regFrame*(!(regMinEta[ireg]==mineta))) ) &&
101  (currRegEta <= (regMaxEta[ireg]-regFrame*(!(regMaxEta[ireg]==maxeta))) ) &&
102  (currRegPhi >= (regMinPhi[ireg]+regFrame*(!(regMinPhi[ireg]==minphi))) ) &&
103  (currRegPhi <= (regMaxPhi[ireg]-regFrame*(!(regMaxPhi[ireg]==maxphi))) ) )
104  {
105  newindex = (currRegEta-mineta)*Nphi + currRegPhi-minphi;
106  iterSolution[newindex] = regIterSolution[i1];
107  }
108  }
109  } // end loop over regions
110 
111  if (iterSolution.empty()) return iterSolution;
112 
113  // re-calibrate eventMatrix with solution
114  for (int ievent = 0; ievent<Nevents; ievent++)
115  {
116  myEventMatrix[ievent] = recalibrateEvent(myEventMatrix[ievent], VmaxCeta[ievent], VmaxCphi[ievent], iterSolution);
117  }
118 
119  // save solution into theCalibVector
120  for (int i=0; i<Nchannels; i++)
121  {
122  totalSolution[i] *= iterSolution[i];
123  }
124 
125  } // end loop over nIter
126 
127  return totalSolution;
128 
129 }
130 
131 
132 std::vector<float> HouseholderDecomposition::iterate(const std::vector<std::vector<float> >& eventMatrix, const std::vector<int>& VmaxCeta, const std::vector<int>& VmaxCphi, const std::vector<float>& energyVector, const int& nIter, const bool& normalizeFlag)
133 {
134  Nevents = eventMatrix.size(); // Number of events to calibrate with
135 
136  std::vector<float> totalSolution(Nchannels,1.);
137  std::vector<float> iterSolution;
138  std::vector<std::vector<float> > myEventMatrix(eventMatrix);
139  std::vector<float> myEnergyVector(energyVector);
140 
141  int i, j;
142 
143  // Iterate the correction
144  for (int iter=1;iter<=nIter;iter++)
145  {
146 
147  // if normalization flag is set, normalize energies
148  float sumOverEnergy;
149  if (normalizeFlag)
150  {
151  float scale = 0.;
152 
153  for (i=0; i<Nevents; i++)
154  {
155  sumOverEnergy = 0.;
156  for (j=0; j<Nxtals; j++) {sumOverEnergy += myEventMatrix[i][j];}
157  sumOverEnergy /= myEnergyVector[i];
158  scale += sumOverEnergy;
159  }
160  scale /= Nevents;
161 
162  for (i=0; i<Nevents; i++) {myEnergyVector[i] *= scale;}
163  } // end normalize energies
164 
165 
166 
167  // now the real work starts:
168  iterSolution = iterate(myEventMatrix,VmaxCeta,VmaxCphi,myEnergyVector);
169 
170  if (iterSolution.empty()) return iterSolution;
171 
172  // re-calibrate eventMatrix with solution
173  for (int ievent = 0; ievent<Nevents; ievent++)
174  {
175  myEventMatrix[ievent] = recalibrateEvent(myEventMatrix[ievent], VmaxCeta[ievent], VmaxCphi[ievent], iterSolution);
176  }
177 
178  for (int i=0; i<Nchannels; i++)
179  {
180  // save solution into theCalibVector
181  totalSolution[i] *= iterSolution[i];
182  }
183 
184  } // end iterate correction
185 
186  return totalSolution;
187 }
188 
189 
190 
191 std::vector<float> HouseholderDecomposition::iterate(const std::vector<std::vector<float> >& eventMatrix, const std::vector<int>& VmaxCeta, const std::vector<int>& VmaxCphi, const std::vector<float>& energyVectorOrig)
192 {
193  std::vector<float> solution;
194 
195  Nevents=eventMatrix.size(); // Number of events to calibrate with
196 
197  if (Nchannels > Nevents)
198  {
199  std::cout << "Householder::runIter(): more channels to calibrate than events available. " << std::endl;
200  std::cout << " Nchannels=" << Nchannels << std::endl;
201  std::cout << " Nevents=" << Nevents << std::endl;
202  std::cout << " ****************** ERROR *********************" << std::endl;
203  return solution; // empty vector
204  }
205 
206  // input: eventMatrixOrig - the unzipped matrix
207  eventMatrixOrig = unzipMatrix(eventMatrix,VmaxCeta,VmaxCphi);
208 
209  if (eventMatrixOrig.size() != energyVectorOrig.size())
210  {
211  std::cout << "Householder::runIter(): matrix dimensions non-conformant. " << std::endl;
212  std::cout << " energyVectorOrig.size()=" << energyVectorOrig.size() << std::endl;
213  std::cout << " eventMatrixOrig.size()=" << eventMatrixOrig.size() << std::endl;
214  std::cout << " ****************** ERROR *********************" << std::endl;
215  return solution; // empty vector
216  }
217 
218 
219  int i,j;
221  energyVectorProc = energyVectorOrig; // copy energyVectorOrig vector
222  std::vector<float> e(Nchannels);
223  alpha.assign(Nchannels,0.);
224  pivot.assign(Nchannels,0);
225 
226 
227  //--------------------
228  bool decomposeSuccess = decompose();
229 
230  if( !decomposeSuccess )
231  {
232  std::cout << "Householder::runIter(): Failed: Singular condition in decomposition."<< std::endl;
233  std::cout << "***************** PROBLEM in DECOMPOSITION *************************"<<std::endl;
234  return solution; // empty vector
235  }
236 
237  /* DBL_EPSILON: Difference between 1.0 and the minimum float greater than 1.0 */
238  float mydbleps = 2.22045e-16; //DBL_EPSILON;
239  float etasqr = mydbleps*mydbleps;
240  // std::cout << "LOOK at DBL_EPSILON :" << mydbleps <<std::endl;
241 
242 
243  //--------------------
244  // apply transformations to rhs - find solution vector
245  solution.assign(Nchannels,0.);
246  solve(solution);
247 
248 
249  // compute residual vector energyVectorProc
250  for (i=0;i<Nevents;i++)
251  {
252  energyVectorProc[i] = energyVectorOrig[i];
253  for (j=0; j<Nchannels; j++)
254  {
255  energyVectorProc[i]-=eventMatrixOrig[i][j]*solution[j];
256  }
257  }
258 
259 
260  //--------------------
261  // compute first correction vector e
262  solve(e);
263 
264  float normy0=0.;
265  float norme1=0.;
266  float norme0;
267 
268 
269 
270 
271  for (i=0;i<Nchannels;i++)
272  {
273  normy0 += solution[i] * solution[i];
274  norme1 += e[i] * e[i];
275  }
276 
277 // std::cout << "Householder::runIter(): applying first correction";
278 // std::cout << " normy0 = " << normy0;
279 // std::cout << " norme1 = " << norme1 << std::endl;
280 
281  // not attempt at obtaining the solution is made unless the norm of the first
282  // correction is significantly smaller than the norm of the initial solution
283  if (norme1>(0.0625*normy0))
284  {
285  // std::cout << "Householder::runIter(): first correction is too large. Failed." << std::endl;
286  }
287 
288  // improve the solution
289  for (i=0; i<Nchannels; i++)
290  {
291  solution[i]+=e[i];
292  }
293 
294  // std::cout << "Householder::runIter(): improving solution...." << std::endl;
295 
296  //--------------------
297  // only continue iteration if the correction was significant
298  while (norme1>(etasqr*normy0))
299  {
300  // std::cout << "Householder::runIter(): norme1 = " << norme1 << std::endl;
301 
302  for (i=0; i<Nevents; i++)
303  {
304  energyVectorProc[i] = energyVectorOrig[i];
305  for (j=0; j<Nchannels; j++)
306  {
307  energyVectorProc[i] -= eventMatrixOrig[i][j] * solution[j];
308  }
309  }
310 
311  // compute next correction vector
312  solve(e);
313 
314  norme0 = norme1;
315  norme1 = 0.;
316 
317  for (i=0;i<Nchannels;i++)
318  {
319  norme1+=e[i]*e[i];
320  }
321 
322  // terminate iteration if the norm of the new correction failed to decrease
323  // significantly compared to the norm of the previous correction
324  if (norme1>(0.0625*norme0)) break;
325 
326  // apply correction vector
327  for (i=0;i<Nchannels;i++)
328  {
329  solution[i]+=e[i];
330  }
331  }
332 
333  //clean up
334  eventMatrixOrig.clear();
335  eventMatrixProc.clear();
336  energyVectorProc.clear();
337  alpha.clear();
338  pivot.clear();
339 
340 
341  return solution;
342 }
343 
344 
346 {
347  int i,j,jbar,k;
348  float beta,sigma,alphak,eventMatrixkk;
349  std::vector<float> y(Nchannels);
350  std::vector<float> sum(Nchannels);
351 
352  // std::cout << "Householder::decompose() started" << std::endl;
353 
354  for (j=0;j<Nchannels;j++)
355  {
356  // jth column sum: squared sum for each crystal
357  sum[j]=0.;
358  for (i=0;i<Nevents;i++)
359  sum[j]+=eventMatrixProc[i][j]*eventMatrixProc[i][j];
360 
361  // bookkeeping vector
362  pivot[j] = j;
363  }
364 
365  for (k=0;k<Nchannels;k++)
366  {
367  // kth Householder transformation
368  sigma = sum[k];
369  jbar = k;
370 
371  // go through all following columns
372  // find the largest sumSquared in the following columns
373  for (j=k+1;j<Nchannels;j++)
374  {
375  if (sum[j] > sigma)
376  {
377  sigma = sum[j];
378  jbar=j;
379  }
380  }
381 
382  if (jbar != k)
383  {
384  // column interchange:
385  // interchange within: bookkeeping vector, squaredSum, eventMatrixProc
386 
387  i = pivot[k];
388  pivot[k]=pivot[jbar];
389  pivot[jbar]=i;
390 
391  sum[jbar]=sum[k];
392  sum[k]=sigma;
393 
394  for (i=0;i<Nevents;i++)
395  {
396  sigma=eventMatrixProc[i][k];
398  eventMatrixProc[i][jbar]=sigma;
399  }
400  } // end column interchange
401 
402  // now store in sigma the squared sum of the readoutEnergies for this column(crystal)
403  sigma=0.;
404  for (i=k;i<Nevents;i++)
405  {
406  sigma+=eventMatrixProc[i][k]*eventMatrixProc[i][k];
407  }
408 
409  // found a zero-column, bail out
410  if (sigma == 0.)
411  {
412 // std::cout << "Householder::decompose() failed" << std::endl;
413 // return false;
414  // rof 14.12.2006: workaround to avoid failure of algorithm because of dead crystals:
415  sigma = sigmaReplacement;
416  // std::cout << "Householder::decompose - found zero column " << jbar << ", replacing sum of column elements by " << sigma << std::endl;
417  }
418 
419 
420  // the following paragraph acts only on the diagonal element:
421  // if element=0, then calculate alpha & beta
422 
423  // take the diagonal element
424  eventMatrixkk = eventMatrixProc[k][k];
425 
426  if (eventMatrixkk < 0.)
427  alpha[k] = sqrt(sigma);
428  else
429  alpha[k] = sqrt(sigma) * (-1.);
430 
431  alphak = alpha[k];
432 
433  beta = 1 / (sigma - eventMatrixkk * alphak);
434  // replace it
435  eventMatrixProc[k][k] = eventMatrixkk - alphak;
436 
437  for (j=k+1; j<Nchannels; j++)
438  {
439  y[j] = 0.;
440 
441  for (i=k; i<Nevents; i++)
442  {
443  y[j] += eventMatrixProc[i][k] * eventMatrixProc[i][j];
444  }
445 
446  y[j] *= beta;
447  }
448 
449  for (j=k+1; j<Nchannels; j++)
450  {
451  for (i=k; i<Nevents; i++)
452  {
453  eventMatrixProc[i][j] -= eventMatrixProc[i][k] * y[j];
454  sum[j] -= eventMatrixProc[k][j] * eventMatrixProc[k][j];
455  }
456  }
457  } // end of kth householder transformation
458 
459  // std::cout << "Householder::decompose() finished" << std::endl;
460 
461  return true;
462 }
463 
464 
465 void HouseholderDecomposition::solve(std::vector<float> &y)
466 {
467  std::vector<float> z(Nchannels,0.);
468 
469  float gamma;
470  int i,j;
471 
472  // std::cout << "Householder::solve() begin" << std::endl;
473 
474  for (j=0; j<Nchannels; j++)
475  {
476  // apply jth transformation to the right hand side
477  gamma = 0.;
478  for (i=j; i<Nevents; i++)
479  {
480  gamma += eventMatrixProc[i][j] * energyVectorProc[i];
481  }
482  gamma /= (alpha[j] * eventMatrixProc[j][j]);
483 
484  for (i=j; i<Nevents; i++)
485  {
486  energyVectorProc[i] += gamma * eventMatrixProc[i][j];
487  }
488  }
489 
490  z[Nchannels-1] = energyVectorProc[Nchannels-1] / alpha[Nchannels-1];
491 
492  for (i=Nchannels-2; i>=0; i--)
493  {
494  z[i] = energyVectorProc[i];
495  for (j=i+1; j<Nchannels; j++)
496  {
497  z[i] -= eventMatrixProc[i][j]*z[j];
498  }
499  z[i] /= alpha[i];
500  }
501 
502  for (i=0; i<Nchannels; i++)
503  {
504  y[pivot[i]] = z[i];
505  }
506 
507  // std::cout << "Householder::solve() finished." << std::endl;
508 
509 }
510 
511 
512 std::vector<float> HouseholderDecomposition::recalibrateEvent(const std::vector<float>& eventSquare, const int& maxCeta, const int& maxCphi, const std::vector<float>& recalibrateVector)
513 {
514  std::vector<float> newEventSquare(eventSquare);
515  int iFull;
516 
517  for (int i=0; i<Nxtals; i++)
518  {
519  iFull = indexSqr2Reg(i, maxCeta, maxCphi);
520  if (iFull >=0)
521  newEventSquare[i] *= recalibrateVector[iFull];
522  }
523  return newEventSquare;
524 }
525 
526 
527 int HouseholderDecomposition::indexSqr2Reg(const int& sqrIndex, const int& maxCeta, const int& maxCphi)
528 {
529  int regionIndex;
530 
531  // get the current eta, phi indices
532  int curr_eta = maxCeta - squareMode/2 + sqrIndex%squareMode;
533  if (curr_eta * maxCeta <= 0) {if (maxCeta > 0) curr_eta--; else curr_eta++; } // JUMP over 0
534 
535  int curr_phi = maxCphi - squareMode/2 + sqrIndex/squareMode;
536  if (curr_phi < 1) curr_phi += 360;
537  if (curr_phi > 360) curr_phi -= 360;
538 
539  bool negPhiDirection = (maxphi < minphi);
540  int iFullphi;
541 
542  regionIndex = -1;
543 
544  if (curr_eta >= mineta && curr_eta <= maxeta)
545  if ( (!negPhiDirection && (curr_phi >= minphi && curr_phi <= maxphi)) ||
546  (negPhiDirection && !(curr_phi >= minphi && curr_phi <= maxphi)) )
547  {
548  iFullphi = curr_phi - minphi;
549  if (iFullphi < 0) iFullphi += 360;
550  regionIndex = (curr_eta - mineta) * (maxphi - minphi + 1 + 360*negPhiDirection) + iFullphi;
551  }
552 
553  return regionIndex;
554 }
555 
556 std::vector<std::vector<float> > HouseholderDecomposition::unzipMatrix(const std::vector< std::vector<float> >& eventMatrix, const std::vector<int>& VmaxCeta, const std::vector<int>& VmaxCphi)
557 {
558  std::vector< std::vector<float> > fullMatrix;
559 
560  int iFull;
561 
562  for (int i=0; i<Nevents; i++)
563  {
564  std::vector<float> foo(Nchannels,0.);
565  for (int k=0; k<Nxtals; k++)
566  {
567  iFull = indexSqr2Reg(k, VmaxCeta[i], VmaxCphi[i]);
568  if (iFull >=0)
569  foo[iFull] = eventMatrix[i][k];
570  }
571  fullMatrix.push_back(foo);
572  }
573 
574  return fullMatrix;
575 }
576 
577 void HouseholderDecomposition::makeRegions(const int& regLength)
578 {
579  // int regFrame = regLength/2;
580  int regFrame = squareMode/2;
581  // first eta:
582  int remEta = Neta % regLength;
583  if (remEta > regLength/2) remEta -= regLength;
584  int numSubRegEta = Neta / regLength + (remEta < 0)*1;
585 
586  int addtoEta = remEta / numSubRegEta;
587  int addtomoreEta = remEta % numSubRegEta; // add "addtomore" number of times (addto+1), remaining times add just (addto)
588 
589  // then phi:
590  int remPhi = Nphi % regLength;
591  if (remPhi > regLength/2) remPhi -= regLength;
592  int numSubRegPhi = Nphi / regLength + (remPhi < 0)*1;
593 
594  int addtoPhi = remPhi / numSubRegPhi;
595  int addtomorePhi = remPhi % numSubRegPhi; // add "addtomore" number of times (addto+-1), remaining times add just (addto)
596 
597  // now put it all together
598  int addedLengthEta = 0;
599  int addedLengthPhi = 0;
600  int startIndexEta = mineta;
601  int startIndexPhi;
602  int endIndexEta = 0;
603  int endIndexPhi;
604  for (int i=0; i < numSubRegEta; i++)
605  {
606  addedLengthEta = regLength + addtoEta + addtomoreEta/abs(addtomoreEta)*(i<abs(addtomoreEta));
607  endIndexEta = startIndexEta + addedLengthEta - 1;
608 
609  startIndexPhi = minphi;
610  endIndexPhi = 0;
611  for (int j=0; j < numSubRegPhi; j++)
612  {
613  addedLengthPhi = regLength + addtoPhi + addtomorePhi/abs(addtomorePhi)*(j<abs(addtomorePhi));
614  endIndexPhi = startIndexPhi + addedLengthPhi - 1;
615 
616  regMinPhi.push_back(startIndexPhi - regFrame*(j!=0) );
617  regMaxPhi.push_back(endIndexPhi + regFrame*(j!=(numSubRegPhi-1)) );
618  regMinEta.push_back(startIndexEta - regFrame*(i!=0) );
619  regMaxEta.push_back(endIndexEta + regFrame*(i!=(numSubRegEta-1)) );
620 
621  startIndexPhi = endIndexPhi + 1;
622  }
623  startIndexEta = endIndexEta + 1;
624  }
625 
626 // // print it all
627 // std::cout << "Householder::makeRegions created the following regions for calibration:" << std::endl;
628 // for (int i=0; i<regMinEta.size(); i++)
629 // std::cout << "Region " << i << ": eta = " << regMinEta[i] << " to " << regMaxEta[i] << ", phi = " << regMinPhi[i] << " to " << regMaxPhi[i] << std::endl;
630 
631 }
const double beta
void solve(std::vector< float > &y)
std::vector< float > recalibrateEvent(const std::vector< float > &eventSquare, const int &maxCeta, const int &maxCphi, const std::vector< float > &recalibrateVector)
Recalibrate before next iteration: give previous solution vector as argument.
void makeRegions(const int &regLength)
std::vector< float > iterate(const std::vector< std::vector< float > > &eventMatrix, const std::vector< int > &VmaxCeta, const std::vector< int > &VmaxCphi, const std::vector< float > &energyVector, const int &nIter, const bool &normalizeFlag=false)
T sqrt(T t)
Definition: SSEVec.h:18
int indexSqr2Reg(const int &sqrIndex, const int &maxCeta, const int &maxCphi)
Method to translate from square indices to region indices.
std::vector< std::vector< float > > eventMatrixProc
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< std::vector< float > > eventMatrixOrig
int k[5][pyjets_maxn]
std::vector< float > energyVectorProc
std::vector< float > runRegional(const std::vector< std::vector< float > > &eventMatrix, const std::vector< int > &VmaxCeta, const std::vector< int > &VmaxCphi, const std::vector< float > &energyVector, const int &nIter, const int &regLength=5)
std::vector< std::vector< float > > unzipMatrix(const std::vector< std::vector< float > > &eventMatrix, const std::vector< int > &VmaxCeta, const std::vector< int > &VmaxCphi)
Unzips the skimmed matrix into a full matrix.
HouseholderDecomposition(int squareMode_=5, int mineta_=1, int maxeta_=85, int minphi_=1, int maxphi_=20)
Default constructor.