CMS 3D CMS Logo

Public Member Functions | Protected Member Functions | Protected Attributes

pftools::LinearCalibrator Class Reference

#include <LinearCalibrator.h>

Inheritance diagram for pftools::LinearCalibrator:
pftools::Calibrator

List of all members.

Public Member Functions

LinearCalibratorclone () const
LinearCalibratorcreate () const
 LinearCalibrator ()
virtual ~LinearCalibrator ()

Protected Member Functions

virtual std::map
< DetectorElementPtr, double > 
getCalibrationCoefficientsCore () throw (PFToolsException&)
virtual TMatrixD & getHessian (const TMatrixD &eij, TMatrixD &hess, const TVectorD &truthE) const
virtual TVectorD & getProjections (const TMatrixD &eij, TVectorD &proj, const TVectorD &truthE) const
virtual void initEijMatrix (TMatrixD &eij, TVectorD &truthE)
 LinearCalibrator (const LinearCalibrator &lc)
virtual void populateDetElIndex ()

Protected Attributes

std::map< DetectorElementPtr,
unsigned > 
myDetElIndex

Detailed Description

Definition at line 24 of file LinearCalibrator.h.


Constructor & Destructor Documentation

LinearCalibrator::LinearCalibrator ( )

Definition at line 42 of file LinearCalibrator.cc.

Referenced by clone(), and create().

                                   {

}
LinearCalibrator::~LinearCalibrator ( ) [virtual]

Definition at line 46 of file LinearCalibrator.cc.

                                    {
}
LinearCalibrator::LinearCalibrator ( const LinearCalibrator lc) [protected]

Member Function Documentation

LinearCalibrator * LinearCalibrator::clone ( void  ) const [virtual]

Implements pftools::Calibrator.

Definition at line 48 of file LinearCalibrator.cc.

References LinearCalibrator().

                                                {
        return new LinearCalibrator(*this);
}
LinearCalibrator * LinearCalibrator::create ( void  ) const [virtual]

Implements pftools::Calibrator.

Definition at line 57 of file LinearCalibrator.cc.

References LinearCalibrator().

                                                 {
        return new LinearCalibrator();
}
std::map< DetectorElementPtr, double > LinearCalibrator::getCalibrationCoefficientsCore ( ) throw (PFToolsException&) [protected, virtual]

Reimplemented from pftools::Calibrator.

Definition at line 61 of file LinearCalibrator.cc.

References gather_cfg::cout, ExpressReco_HICollisions_FallBack::e, getHessian(), getProjections(), pftools::Calibrator::hasParticles(), initEijMatrix(), myDetElIndex, and convertSQLiteXML::ok.

                                         {
//      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;
}
TMatrixD & LinearCalibrator::getHessian ( const TMatrixD &  eij,
TMatrixD &  hess,
const TVectorD &  truthE 
) const [protected, virtual]

Definition at line 179 of file LinearCalibrator.cc.

References i, j, n, and funct::pow().

Referenced by getCalibrationCoefficientsCore().

                                              {
        //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;
}
TVectorD & LinearCalibrator::getProjections ( const TMatrixD &  eij,
TVectorD &  proj,
const TVectorD &  truthE 
) const [protected, virtual]

Definition at line 164 of file LinearCalibrator.cc.

References i, and j.

Referenced by getCalibrationCoefficientsCore().

                                              {
        //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;
}
void LinearCalibrator::initEijMatrix ( TMatrixD &  eij,
TVectorD &  truthE 
) [protected, virtual]

Definition at line 131 of file LinearCalibrator.cc.

References getHLTprescales::index, pftools::Calibrator::myDetectorElements, myDetElIndex, pftools::Calibrator::myParticleDeposits, L1TEmulatorMonitor_cff::p, and populateDetElIndex().

Referenced by getCalibrationCoefficientsCore().

                                                                    {
        //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;
        }

}
void LinearCalibrator::populateDetElIndex ( ) [protected, virtual]

Definition at line 197 of file LinearCalibrator.cc.

References getHLTprescales::index, pftools::Calibrator::myDetectorElements, and myDetElIndex.

Referenced by initEijMatrix().

                                          {
        //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;
                }
        }

}

Member Data Documentation