pftools::LinearCalibrator Class Reference

#include <LinearCalibrator.h>

Inheritance diagram for pftools::LinearCalibrator:

Public Member Functions
LinearCalibrator *	clone () const
LinearCalibrator *	create () const
	LinearCalibrator ()
virtual	~LinearCalibrator ()
Public Member Functions inherited from pftools::Calibrator
void	addDetectorElement (DetectorElementPtr const de)
void	addParticleDeposit (ParticleDepositPtr pd)
	Calibrator ()
std::map< DetectorElementPtr, double >	getCalibrationCoefficients () throw ( PFToolsException&)
std::vector< ParticleDepositPtr >	getParticles ()
int	hasParticles () const
virtual	~Calibrator ()

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
Protected Attributes inherited from pftools::Calibrator
std::vector< DetectorElementPtr >	myDetectorElements
std::vector< ParticleDepositPtr >	myParticleDeposits

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

Definition at line 52 of file LinearCalibrator.cc.

References pftools::Calibrator::myDetectorElements, and pftools::Calibrator::myParticleDeposits.

                                                             {
    myDetectorElements = lc.myDetectorElements;
    myParticleDeposits = lc.myParticleDeposits;
}

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 ( ) 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 & )

protectedvirtual

Reimplemented from pftools::Calibrator.

Definition at line 61 of file LinearCalibrator.cc.

References gather_cfg::cout, alignCSCRings::e, getHessian(), getProjections(), pftools::Calibrator::hasParticles(), initEijMatrix(), hlt_dqm_clientPB-live_cfg::me, 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

protectedvirtual

Definition at line 179 of file LinearCalibrator.cc.

References i, j, gen::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

protectedvirtual

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  )

protectedvirtual

Definition at line 131 of file LinearCalibrator.cc.

References cmsHarvester::index, pftools::Calibrator::myDetectorElements, myDetElIndex, pftools::Calibrator::myParticleDeposits, AlCaHLTBitMon_ParallelJobs::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 ( )

protectedvirtual

Definition at line 197 of file LinearCalibrator.cc.

References cmsHarvester::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

std::map<DetectorElementPtr, unsigned> pftools::LinearCalibrator::myDetElIndex

protected

Definition at line 64 of file LinearCalibrator.h.

Referenced by getCalibrationCoefficientsCore(), initEijMatrix(), and populateDetElIndex().