/data/refman/pasoursint/CMSSW_6_1_2_SLHC4_patch1/src/DataFormats/ParticleFlowReco/src/CaloEllipse.cc

Go to the documentation of this file.

/*
 * CaloEllipse.cc
 *
 *  Created on: 24-Mar-2009
 *      Author: jamie
 */

#include "DataFormats/ParticleFlowReco/interface/CaloEllipse.h"

#include <cmath>


using namespace pftools;

CaloEllipse::CaloEllipse() {


}

CaloEllipse::~CaloEllipse() {

}
double CaloEllipse::getTheta() const {
        if (dataPoints_.size() < 2) {
                return 0.0;
        }

        PointVector meanAdj;

        Point mean = getPosition();

        double sum_dxdx(0.0);
        double sum_dydy(0.0);
        double sum_dxdy(0.0);
        PointCit it = dataPoints_.begin();
        for (; it != dataPoints_.end(); ++it) {
                const Point& p = *it;
                Point q(p.first - mean.first, p.second - mean.second);
                meanAdj.push_back(q);
                sum_dxdx += q.first * q.first;
                sum_dydy += q.second * q.second;
                sum_dxdy += q.first * q.second;
        }
        double theta = 0.5 * atan(2.0 * sum_dxdy / (sum_dydy - sum_dxdx));
        return theta;
}

Point CaloEllipse::getMajorMinorAxes(double sigma) const {

        if (dataPoints_.size() < 2) {
                return Point(0, 0);
        }

        Point mean = getPosition();

        PointCit it = dataPoints_.begin();
        double sum_xx(0.0);
        double sum_yy(0.0);
        double theta = getTheta();

        for (; it != dataPoints_.end(); ++it) {
                const Point& p = *it;
                double X = cos(theta) * (p.first - mean.first) - sin(theta) * (p.second
                                - mean.second);
                double Y = sin(theta) * (p.first - mean.first) + cos(theta) * (p.second
                                - mean.second);
                sum_xx += X * X;
                sum_yy += Y * Y;
        }

        double a = sigma * sqrt(sum_xx / dataPoints_.size());
        double b = sigma * sqrt(sum_yy / dataPoints_.size());

        double major, minor;

        if(a > b) {
                major = a;
                minor = b;
        } else {
                major = b;
                minor = a;
        }

        return Point(major, minor);
}

double CaloEllipse::getEccentricity() const {
        Point p = getMajorMinorAxes();
        double a = p.first;
        double b = p.second;
        if(a == 0)
                return 0;
        double ecc = sqrt((a * a - b * b) / (a * a));
        return ecc;
}

Point CaloEllipse::getPosition() const {

        if (dataPoints_.size() < 1) {
                return Point(0, 0);
        }

        double x_tot(0.0);
        double y_tot(0.0);
        PointCit it = dataPoints_.begin();
        for (; it != dataPoints_.end(); ++it) {
                const Point& p = *it;
                x_tot += p.first;
                y_tot += p.second;

        }

        return Point(x_tot / dataPoints_.size(), y_tot / dataPoints_.size());
}

void CaloEllipse::resetCaches() {
        cachedTheta_ = 0.0;
        cachedMinor_ = 0.0;
        cachedMajor_ = 0.0;

}

void CaloEllipse::makeCaches() {
        cachedTheta_ = getTheta();
        Point axes = getMajorMinorAxes();
        cachedMajor_ = axes.first;
        cachedMinor_ = axes.second;
}

void CaloEllipse::reset() {
        dataPoints_.clear();
        resetCaches();
}

std::ostream& pftools::operator<<(std::ostream& s, const pftools::CaloEllipse& em) {
        s << "CaloEllipse at position = " << toString(em.getPosition()) << ", theta = "
                        << em.getTheta() << ", major/minor axes = " << toString(
                        em.getMajorMinorAxes()) << ", eccentricity = " << em.getEccentricity() << "\n";

        return s;
}