JetCharge/interface/JetCharge.h

0001 #ifndef RecoJets_JetCharge_JetCharge_H
0002 #define RecoJets_JetCharge_JetCharge_H
0003
0004 #include "DataFormats/Math/interface/LorentzVector.h"
0005 #include "DataFormats/Math/interface/Vector.h"
0006 //#include "DataFormats/Candidate/interface/CandidateFwd.h"
0007 #include "DataFormats/TrackReco/interface/Track.h"
0008 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0009 #include "DataFormats/Candidate/interface/Candidate.h"
0010 #include "DataFormats/Candidate/interface/CompositeCandidate.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include <Math/VectorUtil.h>
0013 #include <Math/Rotation3D.h>
0014 #include <Math/RotationZ.h>
0015 #include <Math/RotationY.h>
0016 #include <cmath>
0017
0018 class JetCharge {
0019 public:
0020   enum Variable { Pt, RelPt, RelEta, DeltaR, Unit };
0021   typedef reco::Particle::LorentzVector LorentzVector;
0022   typedef reco::Particle::Vector Vector;
0023
0024   JetCharge(Variable var, double exponent = 1.0) : var_(var), exp_(exponent) {}
0025   JetCharge(const edm::ParameterSet &iCfg);
0026
0027   double charge(const LorentzVector &lv, const reco::TrackCollection &vec) const;
0028   double charge(const LorentzVector &lv, const reco::TrackRefVector &vec) const;
0029   double charge(const LorentzVector &lv, const reco::CandidateCollection &vec) const;
0030   double charge(const reco::Candidate &parent) const;
0031   //double charge(const LorentzVector &lv, const reco::CandidateRefVector &vec) const ;
0032
0033 private:
0034   // functions
0035   template <typename T, typename C>
0036   double chargeFromRef(const LorentzVector &lv, const C &vec) const;
0037   template <typename T, typename C>
0038   double chargeFromVal(const LorentzVector &lv, const C &vec) const;
0039   template <typename T, typename IT>
0040   double chargeFromValIterator(const LorentzVector &lv, const IT &begin, const IT &end) const;
0041
0042   template <class T>
0043   double getWeight(const LorentzVector &lv, const T &obj) const;
0044
0045   // data members
0046   Variable var_;
0047   double exp_;
0048 };
0049 template <typename T, typename IT>
0050 double JetCharge::chargeFromValIterator(const LorentzVector &lv, const IT &begin, const IT &end) const {
0051   double num = 0.0, den = 0.0;
0052   for (IT it = begin; it != end; ++it) {
0053     const T &obj = *it;
0054     double w = getWeight(lv, obj);
0055     den += w;
0056     num += w * obj.charge();
0057   }
0058   return (den > 0.0 ? num / den : 0.0);
0059 }
0060
0061 template <typename T, typename C>
0062 double JetCharge::chargeFromVal(const LorentzVector &lv, const C &vec) const {
0063   typedef typename C::const_iterator IT;
0064   return JetCharge::chargeFromValIterator<T, IT>(lv, vec.begin(), vec.end());
0065 }
0066
0067 template <typename T, typename C>
0068 double JetCharge::chargeFromRef(const LorentzVector &lv, const C &vec) const {
0069   typedef typename C::const_iterator IT;
0070   double num = 0.0, den = 0.0;
0071   for (IT it = vec.begin(), end = vec.end(); it != end; ++it) {
0072     const T &obj = *it;
0073     double w = getWeight(lv, *obj);
0074     den += w;
0075     num += w * obj->charge();
0076   }
0077   return (den > 0.0 ? num / den : 0.0);
0078 }
0079
0080 template <class T>
0081 double JetCharge::getWeight(const LorentzVector &lv, const T &obj) const {
0082   double ret;
0083   switch (var_) {
0084     case Pt:
0085       ret = obj.pt();
0086       break;
0087     case DeltaR:
0088       ret = ROOT::Math::VectorUtil::DeltaR(lv.Vect(), obj.momentum());
0089       break;
0090     case RelPt:
0091     case RelEta:
0092       ret = lv.Vect().Dot(obj.momentum()) / (lv.P() * obj.p());    // cos(theta)
0093       ret = (var_ == RelPt ? std::sqrt(1 - ret * ret) * obj.p() :  // p * sin(theta) = pt
0094                  -0.5 * std::log((1 - ret) / (1 + ret)));          // = - log tan theta/2 = eta
0095       break;
0096     case Unit:
0097     default:
0098       ret = 1.0;
0099   }
0100   return (exp_ == 1.0 ? ret : (ret > 0 ? std::pow(ret, exp_) : -std::pow(std::abs(ret), exp_)));
0101 }
0102
0103 #endif