Tools/src/GenericMinL3Algorithm.cc

0001
0002 /** \file GenericMinL3Algorithm.cc
0003  *
0004  *
0005  * \author R.Ofierzynski, CERN
0006  */
0007
0008 #include "Calibration/Tools/interface/GenericMinL3Algorithm.h"
0009
0010 GenericMinL3Algorithm::GenericMinL3Algorithm(bool normalise)
0011     : normaliseFlag(normalise)
0012 // if normalize=true: for all events sum the e25ovTRP. then take the mean value and rescale all TrP
0013 {}
0014
0015 GenericMinL3Algorithm::~GenericMinL3Algorithm() {}
0016
0017 std::vector<float> GenericMinL3Algorithm::iterate(const std::vector<std::vector<float> >& eventMatrix,
0018                                                   const std::vector<float>& energyVector,
0019                                                   const int nIter) {
0020   std::vector<float> solution;
0021   std::vector<std::vector<float> > myEventMatrix(eventMatrix);
0022   int Nevents = eventMatrix.size();       // Number of events to calibrate with
0023   int Nchannels = eventMatrix[0].size();  // Number of channel coefficients
0024   std::vector<float> theCalibVector(Nchannels, 1.);
0025
0026   // Iterate the correction
0027   for (int iter = 1; iter <= nIter; iter++) {
0028     // make one iteration
0029     solution = iterate(myEventMatrix, energyVector);
0030
0031     if (solution.empty())
0032       return solution;
0033     // R.O.: or throw an exception, what's the standard CMS way ?
0034
0035     // re-calibrate eventMatrix with solution
0036     for (int i = 0; i < Nchannels; i++) {
0037       for (int ievent = 0; ievent < Nevents; ievent++) {
0038         myEventMatrix[ievent][i] *= solution[i];
0039       }
0040       // save solution into theCalibVector
0041       theCalibVector[i] *= solution[i];
0042     }
0043
0044   }  // end iterate the correction
0045
0046   return theCalibVector;
0047 }
0048
0049 std::vector<float> GenericMinL3Algorithm::iterate(const std::vector<std::vector<float> >& eventMatrix,
0050                                                   const std::vector<float>& energyVector) {
0051   std::vector<float> solution;
0052   std::vector<float> myEnergyVector(energyVector);
0053
0054   int Nevents = eventMatrix.size();       // Number of events to calibrate with
0055   int Nchannels = eventMatrix[0].size();  // Number of channel coefficients
0056
0057   // Sanity check
0058   if (Nevents != int(myEnergyVector.size())) {
0059     std::cout << "GenericMinL3Algorithm::iterate(): Error: bad matrix dimensions. Dropping out." << std::endl;
0060     return solution;  // empty vector !
0061   }
0062
0063   // initialize the solution vector with 1.
0064   solution.assign(Nchannels, 1.);
0065
0066   int ievent, i, j;
0067
0068   // if normalization flag is set, normalize energies
0069   float sumOverEnergy;
0070   if (normaliseFlag) {
0071     float scale = 0.;
0072
0073     std::cout << "GenericMinL3Algorithm::iterate(): Normalising event data" << std::endl;
0074
0075     for (i = 0; i < Nevents; i++) {
0076       sumOverEnergy = 0.;
0077       for (j = 0; j < Nchannels; j++) {
0078         sumOverEnergy += eventMatrix[i][j];
0079       }
0080       sumOverEnergy /= myEnergyVector[i];
0081       scale += sumOverEnergy;
0082     }
0083     scale /= Nevents;
0084     std::cout << "  Normalisation = " << scale << std::endl;
0085
0086     for (i = 0; i < Nevents; i++) {
0087       myEnergyVector[i] *= scale;
0088     }
0089   }  // end normalize energies
0090
0091   // This is where the real work goes on...
0092   float sum25, invsum25;
0093   float w;                                  // weight for event
0094   std::vector<float> wsum(Nchannels, 0.);   // sum of weights for a crystal
0095   std::vector<float> Ewsum(Nchannels, 0.);  // sum of products of weight*Etrue/E25
0096
0097   // Loop over events
0098   for (ievent = 0; ievent < Nevents; ievent++) {
0099     // Loop over the 5x5 to find the sum
0100     sum25 = 0.;
0101
0102     for (i = 0; i < Nchannels; i++) {
0103       sum25 += eventMatrix[ievent][i];
0104     }  //*calibs[i];
0105
0106     if (sum25 != 0.) {
0107       invsum25 = 1 / sum25;
0108       // Loop over the 5x5 again and calculate the weights for each xtal
0109       for (i = 0; i < Nchannels; i++) {
0110         w = eventMatrix[ievent][i] * invsum25;  // * calibs[i]
0111         wsum[i] += w;
0112         Ewsum[i] += (w * myEnergyVector[ievent] * invsum25);
0113       }
0114     } else {
0115       std::cout << " Debug: dropping null event: " << ievent << std::endl;
0116     }
0117
0118   }  // end Loop over events
0119
0120   // Apply correction factors to all channels not in the margin
0121   for (i = 0; i < Nchannels; i++) {
0122     if (wsum[i] != 0.) {
0123       solution[i] *= Ewsum[i] / wsum[i];
0124     } else {
0125       std::cout << "warning - no event data for crystal index " << i << std::endl;
0126     }
0127   }
0128
0129   return solution;
0130 }