/Users/jmonk/Physics/ForIA/src/Jet.cxx

Go to the documentation of this file.

#include "ForIA/Jet.hh"
#include "ForIA/JetKey.hh"
#include "ForIA/Track.hh"
#include "ForIA/Utils.hh"
#include "ForIA/Units.hh"
#include <stdexcept>
#include <gsl/gsl_cdf.h>
#include <cmath>

namespace ForIA{
  
  Jet::Jet(JetKeyConstPtr key): IFourMomentum(), m_key(key),
  m_gotpx(false), m_gotpy(false), m_gotpz(false), m_gotRapidity(false),
  m_gotLARQualityMean(false), m_gotChargedFraction(false),
  m_gotBadLooser(false),
  m_gotMatchedTracks(false){
    
  }
  
  Jet::Jet(JetKeyConstPtr &key,
           double E, double pT, double mass, double eta, double phi): IFourMomentum(),
  m_key(key),
  m_E(E), m_PT(pT), m_mass(mass), m_phi(phi), m_eta(eta),
  m_gotpx(false), m_gotpy(false), m_gotpz(false), m_gotRapidity(false),
  m_gotLARQualityMean(false), m_gotChargedFraction(false),
  m_gotBadLooser(false),
  m_gotMatchedTracks(false){
    
  }
  
  double Jet::E() const{
    return m_E;//*(1+getCorrection());
  }
  double Jet::ET() const{
    double et = E()*sin(theta());
    return et;
  }
  double Jet::PT() const{
    return m_PT;//*(1+getCorrection());
  }
  double Jet::mass() const{
    return m_mass;
  }
  double Jet::theta() const{
    double theta = 2*atan(exp(-1.0*m_eta));
    return theta;
  }
  double Jet::phi() const{
    return m_phi;
  }
  double Jet::eta() const{
    return m_eta;
  }
  
  double Jet::px()const{
    if(m_gotpx) return m_px;
    m_px = PT() * sin(phi());
    m_gotpx = true;
    return m_px;
  }

  double Jet::py()const{
    if(m_gotpy) return m_py;
    m_py = PT() * cos(phi());
    m_gotpy = true;
    return m_py;
  }
  
  double Jet::pz()const{
    if(m_gotpz) return m_pz;
    double tanp = tan(theta());
    
    if(fabs(tanp) < 1.e-12){
      m_pz = sqrt(E()*E() - mass()*mass());
    }else{
      m_pz = PT() / tanp;
    }
    
    m_gotpz = true;
    return m_pz;
  }
  
  double Jet::rParam()const{
    return m_key->rParam();
  }
  
  double Jet::rapidity() const{
    if(m_gotRapidity) return m_rapidity;
    double pz = PT() / tan(theta());
    m_rapidity = 0.5 * log((E() + pz)/ (E() - pz));
    m_gotRapidity = true;
    return m_rapidity;
  }
  
  double Jet::HECFraction()const{
    return m_hecFraction;
  }
  
  double Jet::HECQuality()const{
    return m_hecQuality;
  }
  
  double Jet::LArQuality()const{
    return m_larQuality;
  }
  
  double Jet::LArQualityMean()const{
    
    if(m_gotLARQualityMean) return m_LARQualityMean;
    
    // Would it be quicker to bit shift 16 times?
    static const double bit16 = 1./65535.;
    m_LARQualityMean = bit16 * m_averageLARQualityFraction;
    m_gotLARQualityMean = true;
    return m_LARQualityMean;
  }
  
  double Jet::negativeEnergy()const{
    return m_negEnergy;
  }
  
  double Jet::EMFraction()const{
    return m_emFrac;
  }
  
  double Jet::chargedFraction()const{
    if(m_gotChargedFraction) return m_chargedFraction;
    m_chargedFraction = m_trackPTSum / PT();
    m_gotChargedFraction=true;
    return m_chargedFraction;
  }
  
  double Jet::fracSamplingMax()const{
    return m_fracSamplingMax;
  }
  
  bool Jet::isBadLooser()const{
    if(m_gotBadLooser) return m_badLooser;
    
    double modEta = fabs(emScaleEta());
    
    bool hecSpike = ((HECFraction() > 0.5 && fabs(HECQuality()) > 0.5 && LArQualityMean() > 0.8) || 
                     fabs(negativeEnergy()) > 60.*GeV);
    
    bool emNoise = (EMFraction() > 0.95 && fabs(LArQuality()) > 0.8 &&
                    LArQualityMean() > 0.8 && modEta < 2.8);
    
    bool cosmic = ((EMFraction() < 0.05 && chargedFraction() < 0.05 && modEta < 2.) ||
                   (EMFraction() < 0.05 && modEta >= 2.) ||
                   (fracSamplingMax() > 0.99 && modEta < 2.));
    /*
    std::cout<<"eta, EM eta = "<<fabs(eta())<<", "<<modEta<<std::endl;
    
    std::cout<<"HECf, HECQ, LArQMean, negE = " << HECFraction()<<", "<<HECQuality()<<", "<<LArQualityMean()<<", "<<negativeEnergy()<<std::endl;
    std::cout<<"hecSpike = "<<hecSpike<<std::endl;
    
    std::cout<<"EMf, LArQ, LArQMean, eta = "<<EMFraction()<<", "<<LArQuality()<<", "<<LArQualityMean()<<", "<<emscale_eta()<<std::endl;
    std::cout<<"emNoise = "<<emNoise<<std::endl;
    
    std::cout<<"EMf, Chf, eta, FMax = "<<EMFraction()<<", "<<chargedFraction()<<", "<<emscale_eta()<<", "<<fracSamplingMax()<<std::endl;
    std::cout<<"cosmic = "<<cosmic<<std::endl;
    */
    return hecSpike || emNoise || cosmic;
  }
  
  bool Jet::isUgly()const{
    return m_tileGap3 > 0.5 || m_bchCorrCell > 0.5;
  }
  
  int Jet::isBadLoose() const{
    return m_isBadLoose;
  }
  int Jet::isBadMedium() const{
    return m_isBadMedium;
  }
  int Jet::isBadTight() const{
    return m_isBadTight;
  }
//  double Jet::emscale_E() const{
//    return m_emscale_E;
//  }
  double Jet::emScaleEta() const{
    return m_emScaleEta;
  }

  void Jet::scaleEnergy(double scale, bool percentual){
    if(percentual){
      m_E = m_E*(1.+scale);
    }else{
      m_E *= scale;
    }
    
    m_gotpx = false;
    m_gotpy = false;
    m_gotpz = false;
    
  return;
  }
  
  const TrackVector &Jet::matchedTracks()const{
    
    if(m_gotMatchedTracks) return m_matchedTracks;
    
    m_event->matchTracksToJets(*m_key);
    m_gotMatchedTracks = true;
    
    return m_matchedTracks;
  }
  
  
  std::ostream &operator << (std::ostream &out, const Jet &jet){
    out<<"Jet: { E: "<<jet.m_E;
    out<<" , PT: "<<jet.m_PT;
    out<<" , eta: "<<jet.m_eta;
    out<<" , phi: "<<jet.m_phi;
    out<<" , isUgly: "<<jet.isUgly();
    out<<" , isBadLoose: "<<jet.m_isBadLoose;
    out<<" , isBadMedium: "<<jet.m_isBadMedium;
    out<<" , isBadTight: "<<jet.m_isBadTight;
    out<<"}";
    return out;
  }

}