Tools/src/MinL3Algorithm.cc

0001
0002 /** \file MinL3Algorithm.cc
0003  *
0004  * \author R.Ofierzynski, CERN
0005  */
0006
0007 #include "Calibration/Tools/interface/MinL3Algorithm.h"
0008 #include <cmath>
0009
0010 MinL3Algorithm::MinL3Algorithm(float kweight_, int squareMode_, int mineta_, int maxeta_, int minphi_, int maxphi_)
0011     : kweight(kweight_),
0012       squareMode(squareMode_),
0013       mineta(mineta_),
0014       maxeta(maxeta_),
0015       minphi(minphi_),
0016       maxphi(maxphi_),
0017       countEvents(0) {
0018   int Neta = maxeta - mineta + 1;
0019   if (mineta * maxeta < 0)
0020     Neta--;  // there's no eta index = 0
0021   int Nphi = maxphi - minphi + 1;
0022   if (Nphi < 0)
0023     Nphi += 360;
0024
0025   Nchannels = Neta * Nphi;  // no. of channels, get it from edges of the region
0026
0027   Nxtals = squareMode * squareMode;  // no. of xtals in one event
0028
0029   wsum.assign(Nchannels, 0.);
0030   Ewsum.assign(Nchannels, 0.);
0031 }
0032
0033 MinL3Algorithm::~MinL3Algorithm() {}
0034
0035 std::vector<float> MinL3Algorithm::iterate(const std::vector<std::vector<float> >& eventMatrix,
0036                                            const std::vector<int>& VmaxCeta,
0037                                            const std::vector<int>& VmaxCphi,
0038                                            const std::vector<float>& energyVector,
0039                                            const int& nIter,
0040                                            const bool& normalizeFlag) {
0041   int Nevents = eventMatrix.size();  // Number of events to calibrate with
0042
0043   std::vector<float> totalSolution(Nchannels, 1.);
0044   std::vector<float> iterSolution;
0045   std::vector<std::vector<float> > myEventMatrix(eventMatrix);
0046   std::vector<float> myEnergyVector(energyVector);
0047
0048   int i, j;
0049
0050   // Iterate the correction
0051   for (int iter = 1; iter <= nIter; iter++) {
0052     // if normalization flag is set, normalize energies
0053     float sumOverEnergy;
0054     if (normalizeFlag) {
0055       float scale = 0.;
0056
0057       for (i = 0; i < Nevents; i++) {
0058         sumOverEnergy = 0.;
0059         for (j = 0; j < Nxtals; j++) {
0060           sumOverEnergy += myEventMatrix[i][j];
0061         }
0062         sumOverEnergy /= myEnergyVector[i];
0063         scale += sumOverEnergy;
0064       }
0065       scale /= Nevents;
0066
0067       for (i = 0; i < Nevents; i++) {
0068         myEnergyVector[i] *= scale;
0069       }
0070     }  // end normalize energies
0071
0072     // now the real work starts:
0073     for (int iEvt = 0; iEvt < Nevents; iEvt++) {
0074       addEvent(myEventMatrix[iEvt], VmaxCeta[iEvt], VmaxCphi[iEvt], myEnergyVector[iEvt]);
0075     }
0076     iterSolution = getSolution();
0077     if (iterSolution.empty())
0078       return iterSolution;
0079
0080     // re-calibrate eventMatrix with solution
0081     for (int ievent = 0; ievent < Nevents; ievent++) {
0082       myEventMatrix[ievent] = recalibrateEvent(myEventMatrix[ievent], VmaxCeta[ievent], VmaxCphi[ievent], iterSolution);
0083     }
0084
0085     for (int i = 0; i < Nchannels; i++) {
0086       // save solution into theCalibVector
0087       totalSolution[i] *= iterSolution[i];
0088     }
0089     //      resetSolution(); // reset for new iteration, now: getSolution does it automatically if not vetoed
0090   }  // end iterate correction
0091
0092   return totalSolution;
0093 }
0094
0095 void MinL3Algorithm::addEvent(const std::vector<float>& eventSquare,
0096                               const int& maxCeta,
0097                               const int& maxCphi,
0098                               const float& energy) {
0099   countEvents++;
0100
0101   float w, invsumXmatrix;
0102   float eventw;
0103   int iFull, i;
0104   // Loop over the crystal matrix to find the sum
0105   float sumXmatrix = 0.;
0106
0107   for (i = 0; i < Nxtals; i++) {
0108     sumXmatrix += eventSquare[i];
0109   }
0110
0111   // event weighting
0112   eventw = 1 - fabs(1 - sumXmatrix / energy);
0113   eventw = pow(eventw, kweight);
0114
0115   if (sumXmatrix != 0.) {
0116     invsumXmatrix = 1 / sumXmatrix;
0117     // Loop over the crystal matrix (3x3,5x5,7x7) again and calculate the weights for each xtal
0118     for (i = 0; i < Nxtals; i++) {
0119       w = eventSquare[i] * invsumXmatrix;
0120
0121       iFull = indexSqr2Reg(i, maxCeta, maxCphi);
0122       if (iFull >= 0) {
0123         // with event weighting:
0124         wsum[iFull] += w * eventw;
0125         Ewsum[iFull] += (w * eventw * energy * invsumXmatrix);
0126         //        wsum[iFull] += w;
0127         //        Ewsum[iFull] += (w * energy * invsumXmatrix);
0128       }
0129     }
0130   }
0131   //  else {std::cout << " Debug: dropping null event: " << countEvents << std::endl;}
0132 }
0133
0134 std::vector<float> MinL3Algorithm::getSolution(bool resetsolution) {
0135   std::vector<float> solution(Nchannels, 1.);
0136
0137   for (int i = 0; i < Nchannels; i++) {
0138     if (wsum[i] != 0.) {
0139       solution[i] *= Ewsum[i] / wsum[i];
0140     }
0141     //      else
0142     //  { std::cout << "warning - no event data for crystal index (reduced region) " << i << std::endl; }
0143   }
0144
0145   if (resetsolution)
0146     resetSolution();
0147
0148   return solution;
0149 }
0150
0151 void MinL3Algorithm::resetSolution() {
0152   wsum.assign(Nchannels, 0.);
0153   Ewsum.assign(Nchannels, 0.);
0154 }
0155
0156 std::vector<float> MinL3Algorithm::recalibrateEvent(const std::vector<float>& eventSquare,
0157                                                     const int& maxCeta,
0158                                                     const int& maxCphi,
0159                                                     const std::vector<float>& recalibrateVector) {
0160   std::vector<float> newEventSquare(eventSquare);
0161   int iFull;
0162
0163   for (int i = 0; i < Nxtals; i++) {
0164     iFull = indexSqr2Reg(i, maxCeta, maxCphi);
0165     if (iFull >= 0)
0166       newEventSquare[i] *= recalibrateVector[iFull];
0167   }
0168   return newEventSquare;
0169 }
0170
0171 int MinL3Algorithm::indexSqr2Reg(const int& sqrIndex, const int& maxCeta, const int& maxCphi) {
0172   int regionIndex;
0173
0174   // get the current eta, phi indices
0175   int curr_eta = maxCeta - squareMode / 2 + sqrIndex % squareMode;
0176   if (curr_eta * maxCeta <= 0) {
0177     if (maxCeta > 0)
0178       curr_eta--;
0179     else
0180       curr_eta++;
0181   }  // JUMP over 0
0182
0183   int curr_phi = maxCphi - squareMode / 2 + sqrIndex / squareMode;
0184   if (curr_phi < 1)
0185     curr_phi += 360;
0186   if (curr_phi > 360)
0187     curr_phi -= 360;
0188
0189   bool negPhiDirection = (maxphi < minphi);
0190   int iFullphi;
0191
0192   regionIndex = -1;
0193
0194   if (curr_eta >= mineta && curr_eta <= maxeta)
0195     if ((!negPhiDirection && (curr_phi >= minphi && curr_phi <= maxphi)) ||
0196         (negPhiDirection && !(curr_phi >= minphi && curr_phi <= maxphi))) {
0197       iFullphi = curr_phi - minphi;
0198       if (iFullphi < 0)
0199         iFullphi += 360;
0200       regionIndex = (curr_eta - mineta) * (maxphi - minphi + 1 + 360 * negPhiDirection) + iFullphi;
0201     }
0202
0203   return regionIndex;
0204 }