IMACalibBlock.cc

Go to the documentation of this file.

 
#include "Calibration/EcalCalibAlgos/interface/IMACalibBlock.h"
#include "Calibration/Tools/interface/BlockSolver.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/isFinite.h"
#include "TH1F.h"
#include "TFile.h"
#include <cassert>
 
// -----------------------------------------------------
 
IMACalibBlock::IMACalibBlock(const int numberOfElements)
    : VEcalCalibBlock(numberOfElements), m_kaliVector(m_numberOfElements), m_kaliMatrix(evalX2Size()) {
  reset();
}
 
// -----------------------------------------------------
 
IMACalibBlock::~IMACalibBlock() {}
 
// -----------------------------------------------------
 
void IMACalibBlock::Fill(std::map<int, double>::const_iterator MapBegin,
                         std::map<int, double>::const_iterator MapEnd,
                         double pTk,
                         double pSubtract,
                         double sigma) {
  //  std::cerr<<"\n\nfaccio il fill!\n";
  double inverror = 1. / sigma;
  //LP fist loop over energies
  for (std::map<int, double>::const_iterator itMap1 = MapBegin; itMap1 != MapEnd; ++itMap1) {
    //      std::cerr<<"numero "<<itMap1->first<<" vale "<<itMap1->second<<"\t";
    for (std::map<int, double>::const_iterator itMap2 = itMap1; itMap2 != MapEnd; ++itMap2) {
      //LP calculate the chi square value
      double dummy = itMap1->second * itMap2->second;
      dummy *= inverror;
      //LP fill the calib matrix
      m_kaliMatrix.at(itMap1->first * m_numberOfElements + itMap2->first) += dummy;
    }  //LP second loop over xtals
 
    //LP calculate the vector value
    double dummy = pTk;
    dummy -= pSubtract;  //LP borders
    dummy *= itMap1->second;
    dummy *= inverror;
    //LP fill the calib vector
    m_kaliVector.at(itMap1->first) += dummy;
  }  //LP first loop over energies
  return;
}
 
//------------------------------------------------------------
 
void IMACalibBlock::complete() {
  int bef;
  int aft;
  for (unsigned int i = 0; i < m_numberOfElements; ++i) {
    for (unsigned int j = i + 1; j < m_numberOfElements; ++j) {
      bef = (i * m_numberOfElements + j);
      aft = (j * m_numberOfElements + i);
      m_kaliMatrix.at(aft) = m_kaliMatrix.at(bef);
    }  //LP second loop over xtals
  }    //LP first loop over xtals
 
  return;
}
 
// ------------------------------------------------------------
 
int IMACalibBlock::solve(int usingBlockSolver, double min, double max) {
  complete();
  // TH1F vettore ("vettore","vettore",10,-0.1,9.9);
  // TH1F matrice ("matrice","matrice",100,-0.1,99.9);
  CLHEP::HepMatrix kaliMatrix(m_numberOfElements, m_numberOfElements);
  // for (std::vector<double>::iterator it = m_kaliVector.begin();
  //         it!= m_kaliVector.end();++it)
  //     vettore.Fill(it-m_kaliVector.begin(),*it);
  riempiMtr(m_kaliMatrix, kaliMatrix);
  // for (std::vector<double>::iterator it = m_kaliMatrix.begin();
  //         it!= m_kaliMatrix.end();++it)
  //     matrice.Fill(it-m_kaliMatrix.begin(),*it);
  //  TFile f ("fileInteressante.root","RECREATE");
  //  vettore.Write();
  // matrice.Write();
  CLHEP::HepVector kaliVector(m_numberOfElements);
  riempiVtr(m_kaliVector, kaliVector);
  //PG linear system solution
  CLHEP::HepVector result = CLHEP::solve(kaliMatrix, kaliVector);
  for (int i = 0; i < kaliVector.num_row(); ++i)
    if (result.normsq() < min * kaliVector.num_row() || result.normsq() > max * kaliVector.num_row()) {
      if (usingBlockSolver) {
        edm::LogWarning("IML") << "using  blocSlover " << std::endl;
        BlockSolver()(kaliMatrix, kaliVector, result);
      } else {
        edm::LogWarning("IML") << "coeff out of range " << std::endl;
        for (int i = 0; i < kaliVector.num_row(); ++i)
          result[i] = 1.;
      }
    }
  fillMap(result);
  return 0;
}
 
//------------------------------------------------------------
 
inline int IMACalibBlock::evalX2Size() { return m_numberOfElements * m_numberOfElements; }
 
// ------------------------------------------------------------
 
void IMACalibBlock::riempiMtr(const std::vector<double>& piena, CLHEP::HepMatrix& vuota) {
  unsigned int max = m_numberOfElements;
 
  assert(piena.size() == max * max);
  assert(vuota.num_row() == int(max));
  assert(vuota.num_col() == int(max));
  for (unsigned int i = 0; i < max; ++i)
    for (unsigned int j = 0; j < max; ++j)
      if (edm::isNotFinite(piena[i * max + j]))
        vuota[i][j] = 0.;
      else
        vuota[i][j] = piena[i * max + j];
 
  return;
}
 
// ------------------------------------------------------------
 
void IMACalibBlock::riempiVtr(const std::vector<double>& pieno, CLHEP::HepVector& vuoto) {
  int max = m_numberOfElements;
  assert(vuoto.num_row() == max);
  for (int i = 0; i < max; ++i)
    if (edm::isNotFinite(pieno[i]))
      vuoto[i] = 0.;
    else
      vuoto[i] = pieno[i];
 
  return;
}
 
// ------------------------------------------------------------
 
void IMACalibBlock::reset() {
  for (std::vector<double>::iterator vecIt = m_kaliVector.begin(); vecIt != m_kaliVector.end(); ++vecIt) {
    *vecIt = 0.;
  }
 
  for (std::vector<double>::iterator vecIt = m_kaliMatrix.begin(); vecIt != m_kaliMatrix.end(); ++vecIt) {
    *vecIt = 0.;
  }
}
 
// ------------------------------------------------------------
//LP sta provando a fare i seguenti due metodi come locali. Speriamo di non fare stronzate.
 
void IMACalibBlock::fillMap(const CLHEP::HepVector& result) {
  for (unsigned int i = 0; i < m_numberOfElements; ++i) {
    m_coefficients[i] = result[i];
  }
 
  return;
}