/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoParticleFlow/PFClusterTools/src/LinearCalibrator.cc

Go to the documentation of this file.

#include "RecoParticleFlow/PFClusterTools/interface/LinearCalibrator.h"
#include <cassert>
#include <cmath>
#include "TVector.h"
#include "TDecompLU.h"
#include "TDecompSVD.h"

#include "TDecompBK.h"

#include <iostream>

using namespace pftools;

/* These utility functions allow you to print matrices and vectors to a stream,
 * and copy and paste the output directly into Octave or Matlab
 */
void printMat(std::ostream& s, const TMatrixD& input) {
        s << "\t[";
        for (int i(0); i < input.GetNrows(); i++) {
                for (int j(0); j < input.GetNcols(); j++) {
                        s << input[i][j]<< ", \t";
                }
                if (i != (input.GetNrows() - 1)) {
                        s << ";\n";
                        s << "\t";
                }
        }
        s << "\t]\n";
}

void printVec(std::ostream& s, const TVectorD& input) {
        s << "\t[";
        for (int i(0); i < input.GetNrows(); i++) {
                s << input[i];
                if (i != (input.GetNrows() - 1)) {
                        s << ",\n\t";
                }
        }
        s << "\t]\n";
}

LinearCalibrator::LinearCalibrator() {

}

LinearCalibrator::~LinearCalibrator() {
}
LinearCalibrator* LinearCalibrator::clone() const {
        return new LinearCalibrator(*this);
}

LinearCalibrator::LinearCalibrator(const LinearCalibrator& lc) {
        myDetectorElements = lc.myDetectorElements;
        myParticleDeposits = lc.myParticleDeposits;
}

LinearCalibrator* LinearCalibrator::create() const {
        return new LinearCalibrator();
}

std::map<DetectorElementPtr, double> LinearCalibrator::getCalibrationCoefficientsCore()
                throw(PFToolsException&) {
//      std::cout << __PRETTY_FUNCTION__
//                      << ": determining linear calibration coefficients...\n";
        if (!hasParticles()) {
                //I have no particles to calibrate to - throw exception.
                PFToolsException me("Calibrator has no particles for calibration!");
                throw me;
        }
        //std::cout << "\tGetting eij matrix...\n";
        TMatrixD eij;
        TVectorD truthE;
        initEijMatrix(eij, truthE);

//      std::cout << "\tEij matrix:\n";
//      printMat(std::cout, eij);

        //std::cout << "\tGetting projections...\n";
        TVectorD proj;
        TMatrixD hess;

        getProjections(eij, proj, truthE);
        //std::cout << "\tProjections:\n";
        //printVec(std::cout, proj);
        getHessian(eij, hess, truthE);

        TDecompLU lu;
        lu.SetMatrix(hess);
        //std::cout << "\tHessian:\n";
        //printMat(std::cout, hess);

        lu.SetTol(1e-25);
        TMatrixD hessInv = lu.Invert();
        //std::cout <<"\tInverse Hessian:\n";
        //printMat(std::cout, hessInv);
        TVectorD calibsSolved(eij.GetNcols());

        bool ok(true);
        calibsSolved = lu.Solve(proj, ok);
        if (ok) {
                //std::cout << "\tLU reports ok.\n";
        }
        else {
                std::cout << "\tWARNING: LU reports NOT ok.\n";
                //This usually happens when you've asked the calibrator to solve for the 'a' term, without including
                //a dummy 'a' term in each particle deposit.
                //Make sure you do
                //Deposition dOffset(offset, eta, phi, 1.0);
                //particle->addRecDeposition(dOffset);
                //particle->addTruthDeposition(dOffset);
                std::cout
                                << "\tDid you forget to add a dummy offset deposition to each particle?\n";
                std::cout << "\tThrowing an exception!"<< std::endl;
                PFToolsException
                                me("TDecompLU did not converge successfully when finding calibrations. Did you forget to add a dummy offset deposition to each particle?");
                throw me;
        }

        //std::cout << "\tCalibrations: \n";
        //printVec(std::cout, calibsSolved);

        std::map<DetectorElementPtr, double> answers;
        for (std::map<DetectorElementPtr, unsigned>::iterator it =
                        myDetElIndex.begin(); it != myDetElIndex.end(); ++it) {
                DetectorElementPtr de = (*it).first;
                answers[de] = calibsSolved[(*it).second];
        }
        return answers;
}

void LinearCalibrator::initEijMatrix(TMatrixD& eij, TVectorD& truthE) {
        //std::cout << __PRETTY_FUNCTION__ << "\n";
        //std::cout << "\tGetting detector element indices...\n";
        populateDetElIndex();
        eij.Clear();
        eij.Zero();

        truthE.ResizeTo(myParticleDeposits.size());
        truthE.Zero();

        //First determine number of calibration constants.

        eij.ResizeTo(myParticleDeposits.size(), myDetElIndex.size());

        //loop over all particle deposits
        unsigned index(0);
        for (std::vector<ParticleDepositPtr>::const_iterator cit =
                        myParticleDeposits.begin(); cit != myParticleDeposits.end(); ++cit) {
                ParticleDepositPtr p = *cit;
                //get each of the relevant detector elements

                for (std::vector<DetectorElementPtr>::const_iterator detElIt =
                                myDetectorElements.begin(); detElIt != myDetectorElements.end(); ++detElIt) {
                        DetectorElementPtr de = *detElIt;
                        eij[index][myDetElIndex[de]] = p->getRecEnergy(de);
                        //truthE[p->getId()] += p->getTruthEnergy(de);
                }
                truthE[index] = p->getTruthEnergy();
                ++index;
        }

}

TVectorD& LinearCalibrator::getProjections(const TMatrixD& eij, TVectorD& proj,
                const TVectorD& truthE) const {
        //std::cout << __PRETTY_FUNCTION__ << "\n";
        proj.ResizeTo(eij.GetNcols());
        proj.Zero();

        for (int j(0); j < eij.GetNcols(); ++j) {
                for (int i(0); i < eij.GetNrows(); ++i) {
                        proj[j] += eij[i][j] / truthE[i];
                }
        }

        return proj;
}

TMatrixD& LinearCalibrator::getHessian(const TMatrixD& eij, TMatrixD& hess,
                const TVectorD& truthE) const {
        //std::cout << __PRETTY_FUNCTION__ << "\n";
        unsigned nCalibs(eij.GetNcols());
        hess.ResizeTo(nCalibs, nCalibs);
        hess.Zero();

        for (unsigned i(0); i < nCalibs; ++i) {
                for (unsigned j(0); j < nCalibs; ++j) {
                        for (int n(0); n < eij.GetNrows(); ++n) {
                                hess[i][j] += eij[n][i] * eij[n][j]/ pow(truthE[n], 2.0);
                        }
                }
        }

        return hess;
}

void LinearCalibrator::populateDetElIndex() {
        //reserve index = 0 for the constant term, if we're told to compute it
        unsigned index(0);

        myDetElIndex.clear();
        //loop over known detector elements, and assign a unique row/column index to each
        for (std::vector<DetectorElementPtr>::const_iterator cit =
                        myDetectorElements.begin(); cit != myDetectorElements.end(); ++cit) {
                DetectorElementPtr de = *cit;
                //std::cout << "\tGot element: "<< *de;
                //check we don't have duplicate detector elements
                if (myDetElIndex.count(de) == 0) {
                        myDetElIndex[de] = index;
                        ++index;
                }
        }

}