/data/refman/pasoursint/CMSSW_5_3_4/src/Calibration/Tools/src/MinL3Algorithm.cc

Go to the documentation of this file.

#include "Calibration/Tools/interface/MinL3Algorithm.h"
#include <math.h>

MinL3Algorithm::MinL3Algorithm(float kweight_, int squareMode_, int mineta_, int maxeta_, int minphi_, int maxphi_)
  :kweight(kweight_), squareMode(squareMode_), mineta(mineta_), maxeta(maxeta_), minphi(minphi_), maxphi(maxphi_), countEvents(0)
{
  int Neta = maxeta - mineta + 1;
  if (mineta * maxeta < 0) Neta--; // there's no eta index = 0
  int Nphi = maxphi - minphi + 1;
  if (Nphi <0) Nphi += 360;
  
  Nchannels = Neta * Nphi; // no. of channels, get it from edges of the region

  Nxtals = squareMode*squareMode; // no. of xtals in one event

  wsum.assign(Nchannels,0.);
  Ewsum.assign(Nchannels,0.);
}


MinL3Algorithm::~MinL3Algorithm()
{
}


std::vector<float> MinL3Algorithm::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)
{
  int Nevents = eventMatrix.size(); // Number of events to calibrate with

  std::vector<float> totalSolution(Nchannels,1.);
  std::vector<float> iterSolution;
  std::vector<std::vector<float> > myEventMatrix(eventMatrix);
  std::vector<float> myEnergyVector(energyVector);

  int i, j;

  // Iterate the correction
  for (int iter=1;iter<=nIter;iter++) 
    {

      // if normalization flag is set, normalize energies
      float sumOverEnergy;
      if (normalizeFlag)
        {
          float scale = 0.;
          
          for (i=0; i<Nevents; i++)
            {
              sumOverEnergy = 0.;
              for (j=0; j<Nxtals; j++) {sumOverEnergy += myEventMatrix[i][j];}
              sumOverEnergy /= myEnergyVector[i];
              scale += sumOverEnergy;
            }
          scale /= Nevents;
          
          for (i=0; i<Nevents; i++) {myEnergyVector[i] *= scale;}         
        } // end normalize energies

      // now the real work starts:
      for (int iEvt=0; iEvt < Nevents; iEvt++)
        {
          addEvent(myEventMatrix[iEvt], VmaxCeta[iEvt], VmaxCphi[iEvt], myEnergyVector[iEvt]);
        }
      iterSolution = getSolution();
      if (iterSolution.empty()) return iterSolution;

      // re-calibrate eventMatrix with solution
      for (int ievent = 0; ievent<Nevents; ievent++)
        {
          myEventMatrix[ievent] = recalibrateEvent(myEventMatrix[ievent], VmaxCeta[ievent], VmaxCphi[ievent], iterSolution);
        }

      for (int i=0; i<Nchannels; i++) 
        {
          // save solution into theCalibVector
          totalSolution[i] *= iterSolution[i];
        }
      //      resetSolution(); // reset for new iteration, now: getSolution does it automatically if not vetoed
    } // end iterate correction

  return totalSolution;
}


void MinL3Algorithm::addEvent(const std::vector<float>& eventSquare, const int& maxCeta, const int& maxCphi, const float& energy)
{
  countEvents++;

  float w, invsumXmatrix;
  float eventw;
  int iFull, i;
  // Loop over the crystal matrix to find the sum
  float sumXmatrix=0.;
      
  for (i=0; i<Nxtals; i++) { sumXmatrix+=eventSquare[i]; }
      
  // event weighting
  eventw = 1 - fabs(1 - sumXmatrix/energy);
  eventw = pow(eventw,kweight);
      
  if (sumXmatrix != 0.)
    {
      invsumXmatrix = 1/sumXmatrix;
      // Loop over the crystal matrix (3x3,5x5,7x7) again and calculate the weights for each xtal
      for (i=0; i<Nxtals; i++) 
        {               
          w = eventSquare[i] * invsumXmatrix;

          iFull = indexSqr2Reg(i, maxCeta, maxCphi);
          if (iFull >= 0)
            {
              // with event weighting:
              wsum[iFull] += w*eventw;
              Ewsum[iFull] += (w*eventw * energy * invsumXmatrix);
              //              wsum[iFull] += w;
              //              Ewsum[iFull] += (w * energy * invsumXmatrix);
            }
        }
    }
  //  else {std::cout << " Debug: dropping null event: " << countEvents << std::endl;}
}


std::vector<float> MinL3Algorithm::getSolution(bool resetsolution)
{
  std::vector<float> solution(Nchannels,1.);

  for (int i=0; i<Nchannels; i++) 
    {
      if (wsum[i] != 0.) 
        { solution[i]*=Ewsum[i]/wsum[i];}
      //      else 
      //        { std::cout << "warning - no event data for crystal index (reduced region) " << i << std::endl; }
    }
  
  if (resetsolution) resetSolution();

  return solution;
}


void MinL3Algorithm::resetSolution()
{
  wsum.assign(Nchannels,0.);
  Ewsum.assign(Nchannels,0.);
}


std::vector<float> MinL3Algorithm::recalibrateEvent(const std::vector<float>& eventSquare, const int& maxCeta, const int& maxCphi, const std::vector<float>& recalibrateVector)
{
  std::vector<float> newEventSquare(eventSquare);
  int iFull;

  for (int i=0; i<Nxtals; i++) 
    {
      iFull = indexSqr2Reg(i, maxCeta, maxCphi);
      if (iFull >=0)
        newEventSquare[i] *= recalibrateVector[iFull];
    }
  return newEventSquare;
}


int MinL3Algorithm::indexSqr2Reg(const int& sqrIndex, const int& maxCeta, const int& maxCphi)
{
  int regionIndex;

  // get the current eta, phi indices
  int curr_eta = maxCeta - squareMode/2 + sqrIndex%squareMode;
  if (curr_eta * maxCeta <= 0) {if (maxCeta > 0) curr_eta--; else curr_eta++; }  // JUMP over 0

  int curr_phi = maxCphi - squareMode/2 + sqrIndex/squareMode;
  if (curr_phi < 1) curr_phi += 360;
  if (curr_phi > 360) curr_phi -= 360;

  bool negPhiDirection = (maxphi < minphi);
  int iFullphi;

  regionIndex = -1;

  if (curr_eta >= mineta && curr_eta <= maxeta)
    if ( (!negPhiDirection && (curr_phi >= minphi && curr_phi <= maxphi)) ||
         (negPhiDirection && !(curr_phi >= minphi && curr_phi <= maxphi))      ) 
      {
        iFullphi = curr_phi - minphi;
        if (iFullphi < 0) iFullphi += 360;
        regionIndex = (curr_eta - mineta) * (maxphi - minphi + 1 + 360*negPhiDirection) + iFullphi;
      }

  return regionIndex;
}