EcalCalibAlgos/src/IMACalibBlock.cc

0001 /**
0002 */
0003
0004 #include "Calibration/EcalCalibAlgos/interface/IMACalibBlock.h"
0005 #include "Calibration/Tools/interface/BlockSolver.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/Utilities/interface/isFinite.h"
0008 #include "TH1F.h"
0009 #include "TFile.h"
0010 #include <cassert>
0011
0012 // -----------------------------------------------------
0013
0014 IMACalibBlock::IMACalibBlock(const int numberOfElements)
0015     : VEcalCalibBlock(numberOfElements), m_kaliVector(m_numberOfElements), m_kaliMatrix(evalX2Size()) {
0016   reset();
0017 }
0018
0019 // -----------------------------------------------------
0020
0021 IMACalibBlock::~IMACalibBlock() {}
0022
0023 // -----------------------------------------------------
0024
0025 void IMACalibBlock::Fill(std::map<int, double>::const_iterator MapBegin,
0026                          std::map<int, double>::const_iterator MapEnd,
0027                          double pTk,
0028                          double pSubtract,
0029                          double sigma) {
0030   //  std::cerr<<"\n\nfaccio il fill!\n";
0031   double inverror = 1. / sigma;
0032   //LP fist loop over energies
0033   for (std::map<int, double>::const_iterator itMap1 = MapBegin; itMap1 != MapEnd; ++itMap1) {
0034     //      std::cerr<<"numero "<<itMap1->first<<" vale "<<itMap1->second<<"\t";
0035     for (std::map<int, double>::const_iterator itMap2 = itMap1; itMap2 != MapEnd; ++itMap2) {
0036       //LP calculate the chi square value
0037       double dummy = itMap1->second * itMap2->second;
0038       dummy *= inverror;
0039       //LP fill the calib matrix
0040       m_kaliMatrix.at(itMap1->first * m_numberOfElements + itMap2->first) += dummy;
0041     }  //LP second loop over xtals
0042
0043     //LP calculate the vector value
0044     double dummy = pTk;
0045     dummy -= pSubtract;  //LP borders
0046     dummy *= itMap1->second;
0047     dummy *= inverror;
0048     //LP fill the calib vector
0049     m_kaliVector.at(itMap1->first) += dummy;
0050   }  //LP first loop over energies
0051   return;
0052 }
0053
0054 //------------------------------------------------------------
0055
0056 void IMACalibBlock::complete() {
0057   int bef;
0058   int aft;
0059   for (unsigned int i = 0; i < m_numberOfElements; ++i) {
0060     for (unsigned int j = i + 1; j < m_numberOfElements; ++j) {
0061       bef = (i * m_numberOfElements + j);
0062       aft = (j * m_numberOfElements + i);
0063       m_kaliMatrix.at(aft) = m_kaliMatrix.at(bef);
0064     }  //LP second loop over xtals
0065   }    //LP first loop over xtals
0066
0067   return;
0068 }
0069
0070 // ------------------------------------------------------------
0071
0072 int IMACalibBlock::solve(int usingBlockSolver, double min, double max) {
0073   complete();
0074   // TH1F vettore ("vettore","vettore",10,-0.1,9.9);
0075   // TH1F matrice ("matrice","matrice",100,-0.1,99.9);
0076   CLHEP::HepMatrix kaliMatrix(m_numberOfElements, m_numberOfElements);
0077   // for (std::vector<double>::iterator it = m_kaliVector.begin();
0078   //         it!= m_kaliVector.end();++it)
0079   //     vettore.Fill(it-m_kaliVector.begin(),*it);
0080   riempiMtr(m_kaliMatrix, kaliMatrix);
0081   // for (std::vector<double>::iterator it = m_kaliMatrix.begin();
0082   //         it!= m_kaliMatrix.end();++it)
0083   //     matrice.Fill(it-m_kaliMatrix.begin(),*it);
0084   //  TFile f ("fileInteressante.root","RECREATE");
0085   //  vettore.Write();
0086   // matrice.Write();
0087   CLHEP::HepVector kaliVector(m_numberOfElements);
0088   riempiVtr(m_kaliVector, kaliVector);
0089   //PG linear system solution
0090   CLHEP::HepVector result = CLHEP::solve(kaliMatrix, kaliVector);
0091   for (int i = 0; i < kaliVector.num_row(); ++i)
0092     if (result.normsq() < min * kaliVector.num_row() || result.normsq() > max * kaliVector.num_row()) {
0093       if (usingBlockSolver) {
0094         edm::LogWarning("IML") << "using  blocSlover " << std::endl;
0095         BlockSolver()(kaliMatrix, kaliVector, result);
0096       } else {
0097         edm::LogWarning("IML") << "coeff out of range " << std::endl;
0098         for (int i = 0; i < kaliVector.num_row(); ++i)
0099           result[i] = 1.;
0100       }
0101     }
0102   fillMap(result);
0103   return 0;
0104 }
0105
0106 //------------------------------------------------------------
0107
0108 inline int IMACalibBlock::evalX2Size() { return m_numberOfElements * m_numberOfElements; }
0109
0110 // ------------------------------------------------------------
0111
0112 void IMACalibBlock::riempiMtr(const std::vector<double>& piena, CLHEP::HepMatrix& vuota) {
0113   unsigned int max = m_numberOfElements;
0114
0115   assert(piena.size() == max * max);
0116   assert(vuota.num_row() == int(max));
0117   assert(vuota.num_col() == int(max));
0118   for (unsigned int i = 0; i < max; ++i)
0119     for (unsigned int j = 0; j < max; ++j)
0120       if (edm::isNotFinite(piena[i * max + j]))
0121         vuota[i][j] = 0.;
0122       else
0123         vuota[i][j] = piena[i * max + j];
0124
0125   return;
0126 }
0127
0128 // ------------------------------------------------------------
0129
0130 void IMACalibBlock::riempiVtr(const std::vector<double>& pieno, CLHEP::HepVector& vuoto) {
0131   int max = m_numberOfElements;
0132   assert(vuoto.num_row() == max);
0133   for (int i = 0; i < max; ++i)
0134     if (edm::isNotFinite(pieno[i]))
0135       vuoto[i] = 0.;
0136     else
0137       vuoto[i] = pieno[i];
0138
0139   return;
0140 }
0141
0142 // ------------------------------------------------------------
0143
0144 void IMACalibBlock::reset() {
0145   for (std::vector<double>::iterator vecIt = m_kaliVector.begin(); vecIt != m_kaliVector.end(); ++vecIt) {
0146     *vecIt = 0.;
0147   }
0148
0149   for (std::vector<double>::iterator vecIt = m_kaliMatrix.begin(); vecIt != m_kaliMatrix.end(); ++vecIt) {
0150     *vecIt = 0.;
0151   }
0152 }
0153
0154 // ------------------------------------------------------------
0155 //LP sta provando a fare i seguenti due metodi come locali. Speriamo di non fare stronzate.
0156
0157 void IMACalibBlock::fillMap(const CLHEP::HepVector& result) {
0158   for (unsigned int i = 0; i < m_numberOfElements; ++i) {
0159     m_coefficients[i] = result[i];
0160   }
0161
0162   return;
0163 }