/data/refman/pasoursint/CMSSW_4_4_5_patch3/src/DataFormats/ParticleFlowReco/src/CaloWindow.cc

Go to the documentation of this file.

#include "DataFormats/ParticleFlowReco/interface/CaloWindow.h"
#include "DataFormats/Math/interface/deltaR.h"
#include <TMath.h>

using namespace pftools;
//using namespace edm;
using namespace std;

//Let's do CaloRing first
CaloRing::CaloRing() :
        panes_(1) {
        myPanes_.push_back(0);
}

CaloRing::CaloRing(unsigned nPanes) :
        panes_(nPanes) {
        for (unsigned j(0); j < nPanes; ++j)
                myPanes_.push_back(0);
}

double CaloRing::getEnergy(unsigned pane) const {
        if (pane < panes_)
                return myPanes_[pane];
        else
                return 0;
}

bool CaloRing::setEnergy(unsigned pane, double energy) {
        if (pane < panes_) {
                myPanes_[pane] = energy;
                return true;
        }
        return false;
}

bool CaloRing::addEnergy(unsigned pane, double energy) {
        if (pane < panes_) {
                myPanes_[pane] += energy;
                return true;
        }
        return false;
}

double CaloRing::totalE() const {
        double ans(0.0);
        for (std::vector<double>::const_iterator cit = myPanes_.begin(); cit
                        != myPanes_.end(); ++cit) {
                ans += *cit;
        }
        return ans;
}

void CaloRing::reset() {
        for (std::vector<double>::iterator cit = myPanes_.begin(); cit
                        != myPanes_.end(); ++cit)
                *cit = 0.0;

}

void CaloRing::printEnergies(std::ostream& s, double range) {
        for (std::vector<double>::iterator cit = myPanes_.begin(); cit
                        != myPanes_.end(); ++cit)
                s << (*cit)/range << "\t";
}

std::ostream& pftools::operator<<(std::ostream& s,
                const pftools::CaloRing& caloRing) {
        for (std::vector<double>::const_iterator it = caloRing.myPanes_.begin(); it
                        != caloRing.myPanes_.end(); ++it) {
                s << *it << "\t";
        }
        s << " => ring E = " << caloRing.totalE();
        return s;
}

std::ostream& pftools::operator<<(std::ostream& s,
                const pftools::CaloWindow& caloWindow) {

        s << "CaloWindow at (" << caloWindow.baryEta_ << ", "
                        << caloWindow.baryPhi_ << "):\n";
        double totalE(0.0);
        for (std::map<unsigned, pftools::CaloRing>::const_iterator cit =
                        caloWindow.energies_.begin(); cit != caloWindow.energies_.end(); ++cit) {
                unsigned ring = (*cit).first;
                const CaloRing& cr = (*cit).second;
                s << "Ring " << ring << ":\t" << cr << "\n";
                totalE += cr.totalE();
        }

        s << "\tTotal E = " << totalE << std::endl;
        return s;

}

vector<double> CaloWindow::stream(double normalisation) const {
        vector<double> stream;
        for (map<unsigned, pftools::CaloRing>::const_iterator cit =
                        energies_.begin(); cit != energies_.end(); ++cit) {
                CaloRing c = (*cit).second;
                std::vector<double> ringE = c.getEnergies();
                stream.insert(stream.end(), ringE.begin(), ringE.end());
        }
        if(normalisation != 1.0){
                for(vector<double>::iterator i = stream.begin(); i != stream.end(); ++i) {
                        double& entry = *i;
                        entry /= normalisation;
                }
        }
        return stream;

}

void CaloWindow::printEnergies(std::ostream& s, double range) {

        for (std::map<unsigned, pftools::CaloRing>::const_iterator cit =
                        energies_.begin(); cit != energies_.end(); ++cit) {
                CaloRing c = (*cit).second;
                c.printEnergies(s, range);
                s << "\t";
        }

}

CaloWindow::CaloWindow() :
        baryEta_(0.0), baryPhi_(0.0), nRings_(1), deltaR_(0.1), nPanes_(1),  axis_(0.0) {

}

CaloWindow::CaloWindow(double eta, double phi, unsigned nRings, double deltaR, unsigned nPanes, double axis) {

        init(eta, phi, nRings, deltaR, nPanes, axis);
}

void CaloWindow::reset() {
        energies_.clear();
}

void CaloWindow::init(double eta, double phi, unsigned nRings, double deltaR, unsigned nPanes, double axis) {
        baryEta_ = eta;
        baryPhi_ = phi;
        nRings_ = nRings;
        deltaR_ = deltaR;
        nPanes_ = nPanes;
        axis_ = axis;

        energies_.clear();
        CaloRing c(1);
        energies_[0] = c;
        for (unsigned j(1); j < nRings_; ++j) {

                CaloRing r(nPanes);
                energies_[j] = r;
        }
}

CaloWindow::~CaloWindow() {

}

std::pair<unsigned, unsigned> CaloWindow::relativePosition(double eta,
                double phi) const {
        //How far is this hit from the barycentre in deltaR?
        double dEta = eta - baryEta_;
        double dPhi = deltaPhi(phi, baryPhi_);

        double dR = deltaR(eta, phi, baryEta_, baryPhi_);

        unsigned ring = static_cast<unsigned> (floor(dR / deltaR_));

        if (ring == 0) {
                //LogDebug("CaloWindow") << "Relative position: adding to central ring\n";
                return std::pair<unsigned, unsigned>(0, 0);
        }

        if (ring >= nRings_) {
                return std::pair<unsigned, unsigned>(ring, 0);
        }

        double dTheta = 0;

        if (dEta > 0) {
                if (dPhi > 0)
                        dTheta = TMath::ATan(dPhi / dEta);
                else
                        dTheta = 2 * TMath::Pi() + TMath::ATan(dPhi / dEta);
        } else {
                if (dPhi > 0)
                        dTheta = TMath::Pi() + TMath::ATan(dPhi / dEta);
                else
                        dTheta = TMath::Pi() + TMath::ATan(dPhi / dEta);
        }
        //Rotation. Rotate theta into axis of caloWindow
        // and also by half a window pane.
        //TODO: bug check!!
        //double dThetaOrig(dTheta);
        double paneOn2 = TMath::Pi() / (*energies_.find(ring)).second.getNPanes();
        dTheta =  dTheta - axis_ - paneOn2;
        //Original theta between 0 and 2 Pi, but transform above might move us beyond this, so...
        //double dThetaCopy(dTheta);
        //std::cout << "dTheta " << dThetaOrig << ", axis " << axis_ << ", paneOn2 " << paneOn2 << ", thetaPrime " << dTheta << "\n";
        if(dTheta > 2 *TMath::Pi()) {
                //std::cout << "To infinity and beyond...\n";
                while(dTheta > 2 * TMath::Pi()) {
                        dTheta -= 2 * TMath::Pi();
                        //std::cout << "dTheta1 " << dTheta << "\n";
                }
        }
        else if (dTheta < 0) {
                //std::cout << "To infinity and beyond 2... dTheta = " << dTheta << "\n";
                while(dTheta < 0) {
                        dTheta += 2 * TMath::Pi();
                        //std::cout << "dTheta2 " << dTheta << "\n";
                }
        }

//      std::cout << "\tdTheta is " << dTheta << " rad \n";
        unsigned division = static_cast<unsigned> (floor((*energies_.find(ring)).second.getNPanes() * dTheta
                        / (TMath::Pi() * 2)));
        //LogDebug("CaloWindow") << "Ring is " << ring << ", pane is "
        //              << division << "\n";

        return std::pair<unsigned, unsigned>(ring, division);

}

bool CaloWindow::addHit(double eta, double phi, double energy) {

        std::pair<unsigned, unsigned> position = relativePosition(eta, phi);
        if (position.first >= nRings_) {
/*              double dEta = eta - baryEta_;
                double dPhi = deltaPhi(phi, baryPhi_);
                double dR = deltaR(eta, phi, baryEta_, baryPhi_);
                LogDebug("CaloWindow")
                                << "Hit is outside my range - it would be in ring "
                                << position.first << ", with dR = " << dR << std::endl;*/
                return false;
        }
//      std::cout << "Adding hit to ring " << position.first << " in position "
//                      << position.second << " with energy " << energy << "\n";

        CaloRing& c = energies_[position.first];
        c.addEnergy(position.second, energy);

        return true;
}
std::map<unsigned, double> CaloWindow::getRingEnergySummations() const {
        std::map<unsigned, double> answer;
        for (std::map<unsigned, CaloRing>::const_iterator cit = energies_.begin(); cit
                        != energies_.end(); ++cit) {
                std::pair<unsigned, CaloRing> pair = *cit;
                answer[pair.first] = pair.second.totalE();
        }

        return answer;
}

void TestCaloWindow::doTest() {

        CaloWindow cw(0.0, 0.0, 3, 0.1, 4);
        std::cout << cw << std::endl;
        cw.addHit(0, 0.05, 1.0);
        cw.addHit(0.22, 0.05, 2.0);
        cw.addHit(0, 0.8, 4.0);
        cw.addHit(0.2, 0, 1.0);
        cw.addHit(0.15, 0.15, 2.0);
        cw.addHit(0.0, 0.2, 3.0);
        cw.addHit(-0.15, 0.15, 4.0);
        cw.addHit(-0.2, 0, 5.0);
        cw.addHit(-0.15, -0.15, 6.0);
        cw.addHit(-0.0, -0.2, 7.0);
        cw.addHit(0.15, -0.15, 8.0);

        std::cout << cw << std::endl;
        std::cout << "bye bye" << std::endl;

}