Tools/interface/MinL3AlgoUniv.h

0001 #ifndef MinL3AlgoUniv_H
0002 #define MinL3AlgoUniv_H
0003
0004 /** \class MinL3AlgoUniv
0005  *  Implementation of the L3 Collaboration algorithm to solve a system Ax = B
0006  *  by minimization of |Ax-B| using an iterative linear approach
0007  *  This class should be universal, i.e. working with DetIds or whatever else
0008  *    will be invented to identify Subdetector parts
0009  *  The bookkeeping of the cluster size and its elements has to be done by the user.
0010  *
0011  * \author R.Ofierzynski, CERN
0012  */
0013
0014 #include <vector>
0015 #include <iostream>
0016 #include <map>
0017 #include <cmath>
0018
0019 template <class IDdet>
0020 class MinL3AlgoUniv {
0021 public:
0022   typedef std::map<IDdet, float> IDmap;
0023   typedef typename IDmap::value_type IDmapvalue;
0024   typedef typename IDmap::iterator iter_IDmap;
0025   typedef typename IDmap::const_iterator citer_IDmap;
0026
0027   /// Default constructor
0028   /// kweight_ = event weight
0029   MinL3AlgoUniv(float kweight_ = 0.);
0030
0031   /// Destructor
0032   ~MinL3AlgoUniv();
0033
0034   /// method doing the full calibration running nIter number of times,
0035   ///   recalibrating the event matrix after each iteration with the new solution
0036   /// returns the vector of calibration coefficients built from all iteration solutions
0037   /// >> also to be used also as recipe on how to use the calibration methods <<
0038   /// >> one-by-one with a re-selection of the events in between the iterations<<
0039   IDmap iterate(const std::vector<std::vector<float> >& eventMatrix,
0040                 const std::vector<std::vector<IDdet> >& idMatrix,
0041                 const std::vector<float>& energyVector,
0042                 const int& nIter,
0043                 const bool& normalizeFlag = false);
0044
0045   /// add event to the calculation of the calibration vector
0046   void addEvent(const std::vector<float>& myCluster, const std::vector<IDdet>& idCluster, const float& energy);
0047
0048   /// recalibrate before next iteration: give previous solution vector as argument
0049   std::vector<float> recalibrateEvent(const std::vector<float>& myCluster,
0050                                       const std::vector<IDdet>& idCluster,
0051                                       const IDmap& newCalibration);
0052
0053   /// get the solution at the end of the calibration as a map between
0054   /// DetIds and calibration constant
0055   IDmap getSolution(const bool resetsolution = true);
0056
0057   /// reset for new iteration
0058   void resetSolution();
0059
0060 private:
0061   float kweight;
0062   int countEvents;
0063   IDmap wsum;
0064   IDmap Ewsum;
0065 };
0066
0067 template <class IDdet>
0068 MinL3AlgoUniv<IDdet>::MinL3AlgoUniv(float kweight_) : kweight(kweight_), countEvents(0) {
0069   resetSolution();
0070 }
0071
0072 template <class IDdet>
0073 MinL3AlgoUniv<IDdet>::~MinL3AlgoUniv() {}
0074
0075 template <class IDdet>
0076 typename MinL3AlgoUniv<IDdet>::IDmap MinL3AlgoUniv<IDdet>::iterate(const std::vector<std::vector<float> >& eventMatrix,
0077                                                                    const std::vector<std::vector<IDdet> >& idMatrix,
0078                                                                    const std::vector<float>& energyVector,
0079                                                                    const int& nIter,
0080                                                                    const bool& normalizeFlag) {
0081   int Nevents = eventMatrix.size();  // Number of events to calibrate with
0082
0083   IDmap totalSolution;
0084   IDmap iterSolution;
0085   std::vector<std::vector<float> > myEventMatrix(eventMatrix);
0086   std::vector<float> myEnergyVector(energyVector);
0087
0088   int i;
0089
0090   // Iterate the correction
0091   for (int iter = 1; iter <= nIter; iter++) {
0092     // if normalization flag is set, normalize energies
0093     float sumOverEnergy;
0094     if (normalizeFlag) {
0095       float scale = 0.;
0096
0097       for (i = 0; i < Nevents; i++) {
0098         sumOverEnergy = 0.;
0099         for (unsigned j = 0; j < myEventMatrix[i].size(); j++) {
0100           sumOverEnergy += myEventMatrix[i][j];
0101         }
0102         sumOverEnergy /= myEnergyVector[i];
0103         scale += sumOverEnergy;
0104       }
0105       scale /= Nevents;
0106
0107       for (i = 0; i < Nevents; i++) {
0108         myEnergyVector[i] *= scale;
0109       }
0110     }  // end normalize energies
0111
0112     // now the real work starts:
0113     for (int iEvt = 0; iEvt < Nevents; iEvt++) {
0114       addEvent(myEventMatrix[iEvt], idMatrix[iEvt], myEnergyVector[iEvt]);
0115     }
0116     iterSolution = getSolution();
0117     if (iterSolution.empty())
0118       return iterSolution;
0119
0120     // re-calibrate eventMatrix with solution
0121     for (int ievent = 0; ievent < Nevents; ievent++) {
0122       myEventMatrix[ievent] = recalibrateEvent(myEventMatrix[ievent], idMatrix[ievent], iterSolution);
0123     }
0124
0125     // save solution into theCalibVector
0126     for (iter_IDmap i = iterSolution.begin(); i != iterSolution.end(); i++) {
0127       iter_IDmap itotal = totalSolution.find(i->first);
0128       if (itotal == totalSolution.end()) {
0129         totalSolution.insert(IDmapvalue(i->first, i->second));
0130       } else {
0131         itotal->second *= i->second;
0132       }
0133     }
0134
0135     //      resetSolution(); // reset for new iteration, now: getSolution does it automatically if not vetoed
0136   }  // end iterate correction
0137
0138   return totalSolution;
0139 }
0140
0141 template <class IDdet>
0142 void MinL3AlgoUniv<IDdet>::addEvent(const std::vector<float>& myCluster,
0143                                     const std::vector<IDdet>& idCluster,
0144                                     const float& energy) {
0145   countEvents++;
0146
0147   float w, invsumXmatrix;
0148   float eventw;
0149   // Loop over the crystal matrix to find the sum
0150   float sumXmatrix = 0.;
0151
0152   for (unsigned i = 0; i < myCluster.size(); i++) {
0153     sumXmatrix += myCluster[i];
0154   }
0155
0156   // event weighting
0157   eventw = 1 - fabs(1 - sumXmatrix / energy);
0158   eventw = pow(eventw, kweight);
0159
0160   if (sumXmatrix != 0.) {
0161     invsumXmatrix = 1 / sumXmatrix;
0162     // Loop over the crystal matrix (3x3,5x5,7x7) again and calculate the weights for each xtal
0163     for (unsigned i = 0; i < myCluster.size(); i++) {
0164       w = myCluster[i] * invsumXmatrix;
0165
0166       // include the weights into wsum, Ewsum
0167       iter_IDmap iwsum = wsum.find(idCluster[i]);
0168       if (iwsum == wsum.end())
0169         wsum.insert(IDmapvalue(idCluster[i], w * eventw));
0170       else
0171         iwsum->second += w * eventw;
0172
0173       iter_IDmap iEwsum = Ewsum.find(idCluster[i]);
0174       if (iEwsum == Ewsum.end())
0175         Ewsum.insert(IDmapvalue(idCluster[i], (w * eventw * energy * invsumXmatrix)));
0176       else
0177         iEwsum->second += (w * eventw * energy * invsumXmatrix);
0178     }
0179   }
0180   //  else {std::cout << " Debug: dropping null event: " << countEvents << std::endl;}
0181 }
0182
0183 template <class IDdet>
0184 typename MinL3AlgoUniv<IDdet>::IDmap MinL3AlgoUniv<IDdet>::getSolution(const bool resetsolution) {
0185   IDmap solution;
0186
0187   for (iter_IDmap i = wsum.begin(); i != wsum.end(); i++) {
0188     iter_IDmap iEwsum = Ewsum.find(i->first);
0189     float myValue = 1;
0190     if (i->second != 0)
0191       myValue = iEwsum->second / i->second;
0192
0193     solution.insert(IDmapvalue(i->first, myValue));
0194   }
0195
0196   if (resetsolution)
0197     resetSolution();
0198
0199   return solution;
0200 }
0201
0202 template <class IDdet>
0203 void MinL3AlgoUniv<IDdet>::resetSolution() {
0204   wsum.clear();
0205   Ewsum.clear();
0206 }
0207
0208 template <class IDdet>
0209 std::vector<float> MinL3AlgoUniv<IDdet>::recalibrateEvent(const std::vector<float>& myCluster,
0210                                                           const std::vector<IDdet>& idCluster,
0211                                                           const IDmap& newCalibration) {
0212   std::vector<float> newCluster(myCluster);
0213
0214   for (unsigned i = 0; i < myCluster.size(); i++) {
0215     citer_IDmap icalib = newCalibration.find(idCluster[i]);
0216     if (icalib != newCalibration.end()) {
0217       newCluster[i] *= icalib->second;
0218     } else {
0219       std::cout << "No calibration available for this element." << std::endl;
0220     }
0221   }
0222
0223   return newCluster;
0224 }
0225
0226 #endif  // MinL3AlgoUniv_H