JetCharge.h

Go to the documentation of this file.

#ifndef RecoJets_JetCharge_JetCharge_H
#define RecoJets_JetCharge_JetCharge_H

#include "DataFormats/Math/interface/LorentzVector.h"
#include "DataFormats/Math/interface/Vector.h"
//#include "DataFormats/Candidate/interface/CandidateFwd.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/Candidate/interface/CompositeCandidate.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include <Math/VectorUtil.h>
#include <Math/Rotation3D.h>
#include <Math/RotationZ.h>
#include <Math/RotationY.h>
#include <cmath>

class JetCharge {
public:
    enum Variable { Pt, RelPt, RelEta, DeltaR, Unit };
    typedef reco::Particle::LorentzVector  LorentzVector;
    typedef reco::Particle::Vector         Vector;

    JetCharge(Variable var, double exponent=1.0) : var_(var), exp_(exponent) { }
    JetCharge(const edm::ParameterSet &iCfg) ;

    double charge(const LorentzVector &lv, const reco::TrackCollection &vec) const ;
    double charge(const LorentzVector &lv, const reco::TrackRefVector &vec) const ;
    double charge(const LorentzVector &lv, const reco::CandidateCollection &vec) const ;
    double charge(const reco::Candidate &parent) const ;
    //double charge(const LorentzVector &lv, const reco::CandidateRefVector &vec) const ;

private:
    // functions
        template<typename T, typename C>
    double chargeFromRef(const LorentzVector &lv, const C &vec) const ;
        template<typename T, typename C>
    double chargeFromVal(const LorentzVector &lv, const C &vec) const ;
        template<typename T, typename IT>
        double chargeFromValIterator(const LorentzVector &lv, const IT &begin, const IT &end) const ;

        template <class T>
    double getWeight(const LorentzVector &lv, const T& obj) const ;

    // data members
    Variable var_; double exp_;


};
template<typename T, typename IT>
double JetCharge::chargeFromValIterator(const LorentzVector &lv, const IT &begin, const IT &end) const {
    double num = 0.0, den = 0.0;
    for (IT it = begin; it != end ; ++it) {
        const T & obj = *it;
        double w = getWeight(lv, obj);
        den += w;
        num += w * obj.charge(); 
    }
    return (den > 0.0 ? num/den : 0.0);
}

template<typename T, typename C>
double JetCharge::chargeFromVal(const LorentzVector &lv, const C &vec) const {
    typedef typename C::const_iterator IT;
    return JetCharge::chargeFromValIterator<T,IT>(lv, vec.begin(), vec.end());
}

template<typename T, typename C>
double JetCharge::chargeFromRef(const LorentzVector &lv, const C &vec) const {
    typedef typename C::const_iterator IT;
    double num = 0.0, den = 0.0;
    for (IT it = vec.begin(), end = vec.end(); it != end ; ++it) {
        const T & obj = *it;
        double w = getWeight(lv, *obj);
        den += w;
        num += w * obj->charge();
    }
    return (den > 0.0 ? num/den : 0.0);
}

template <class T>
double JetCharge::getWeight(const LorentzVector &lv, const T& obj) const { 
    double ret;
    switch (var_) {
        case Pt: 
            ret = obj.pt(); 
            break;
        case DeltaR: 
            ret = ROOT::Math::VectorUtil::DeltaR(lv.Vect(), obj.momentum());
            break;
        case RelPt: 
        case RelEta: 
            ret =  lv.Vect().Dot(obj.momentum())/(lv.P() * obj.p()); // cos(theta)
            ret =  (var_ == RelPt ? 
                std::sqrt(1 - ret*ret) * obj.p() :    // p * sin(theta) = pt
            - 0.5 * std::log((1-ret)/(1+ret)));   // = - log tan theta/2 = eta
            break;
        case Unit:
        default:
            ret = 1.0;
    }
    return (exp_ == 1.0 ? ret :  (ret > 0 ? 
                                std::pow(ret,exp_) : 
                                - std::pow(std::abs(ret), exp_)));
}


#endif