Tools/interface/HouseholderDecomposition.h

0001 #ifndef HouseholderDecomposition_H
0002 #define HouseholderDecomposition_H
0003
0004 /** \class HouseholderDecomposition
0005  *  Implementation of the QR decomposition of a matrix using Householder transformation
0006  *
0007  *  13.03.2007: R.Ofierzynski
0008  *   - using a reduced matrix
0009  *   - implements Regional Householder Algorithm
0010  *
0011  * \author Lorenzo Agostino, R.Ofierzynski, CERN
0012  */
0013
0014 #include <vector>
0015 #include <iostream>
0016
0017 class HouseholderDecomposition {
0018 public:
0019   /// Default constructor
0020   HouseholderDecomposition(int squareMode_ = 5, int mineta_ = 1, int maxeta_ = 85, int minphi_ = 1, int maxphi_ = 20);
0021
0022   /// Destructor
0023   ~HouseholderDecomposition();
0024
0025   /// Run Regional HouseholderAlgorithm (fast version), that splits matrix into regional matrices and inverts them separately.
0026   /// Returns the final vector of calibration coefficients.
0027   /// input: eventMatrix - the skimmed event matrix,
0028   ///        VmaxCeta, VmaxCphi - vectors containing eta and phi indices of the maximum containment crystal for each event,
0029   ///        energyVector - the energy vector,
0030   ///        nIter - number of iterations to be performed,
0031   ///        regLength - default length of the region (in eta- and phi-indices), regLength=5 recommended
0032   /// Comment from author: if you use the same events, 2 iterations are recommended; the second iteration gives corrections of the order of 0.0001
0033   std::vector<float> runRegional(const std::vector<std::vector<float> >& eventMatrix,
0034                                  const std::vector<int>& VmaxCeta,
0035                                  const std::vector<int>& VmaxCphi,
0036                                  const std::vector<float>& energyVector,
0037                                  const int& nIter,
0038                                  const int& regLength = 5);
0039
0040   /// Run the Householder Algorithm several times (nIter). Returns the final vector of calibration coefficients.
0041   /// Comment from author: unless you do a new selection in between the iterations you don't need to run more than once;
0042   ///                      a second iteration on the same events does not improve the result in this case
0043   std::vector<float> iterate(const std::vector<std::vector<float> >& eventMatrix,
0044                              const std::vector<int>& VmaxCeta,
0045                              const std::vector<int>& VmaxCphi,
0046                              const std::vector<float>& energyVector,
0047                              const int& nIter,
0048                              const bool& normalizeFlag = false);
0049
0050   /// Run the Householder Algorithm. Returns the vector of calibration coefficients.
0051   std::vector<float> iterate(const std::vector<std::vector<float> >& eventMatrix,
0052                              const std::vector<int>& VmaxCeta,
0053                              const std::vector<int>& VmaxCphi,
0054                              const std::vector<float>& energyVectorOrig);
0055
0056   /// Recalibrate before next iteration: give previous solution vector as argument
0057   std::vector<float> recalibrateEvent(const std::vector<float>& eventSquare,
0058                                       const int& maxCeta,
0059                                       const int& maxCphi,
0060                                       const std::vector<float>& recalibrateVector);
0061
0062   /// Method to translate from square indices to region indices
0063   int indexSqr2Reg(const int& sqrIndex, const int& maxCeta, const int& maxCphi);
0064
0065 private:
0066   /// Make decomposition
0067   /// input: m=number of events, n=number of channels, qr=event matrix
0068   /// output: qr = transformed event matrix, alpha, pivot
0069   /// returns a boolean value, true if decomposition worked, false if it didn't
0070   bool decompose();
0071
0072   /// Apply transformations to rhs
0073   /// output: r = residual vector (energy vector), y = solution
0074   void solve(std::vector<float>& y);
0075
0076   /// Unzips the skimmed matrix into a full matrix
0077   std::vector<std::vector<float> > unzipMatrix(const std::vector<std::vector<float> >& eventMatrix,
0078                                                const std::vector<int>& VmaxCeta,
0079                                                const std::vector<int>& VmaxCphi);
0080
0081   /// Determines the regions used for splitting of the full matrix and calibrating separately
0082   /// used by the public runRegional method
0083   void makeRegions(const int& regLength);
0084
0085   int squareMode, countEvents;
0086   int mineta, maxeta, minphi, maxphi, Neta, Nphi;
0087   int Nchannels, Nxtals, Nevents;
0088   std::vector<std::vector<float> > eventMatrixOrig;
0089   std::vector<std::vector<float> > eventMatrixProc;
0090   std::vector<float> energyVectorProc;
0091   std::vector<float> alpha;
0092   std::vector<int> pivot;
0093
0094   std::vector<int> regMinPhi, regMaxPhi, regMinEta, regMaxEta;
0095   float sigmaReplacement;
0096 };
0097
0098 #endif  // HouseholderDecomposition_H