JetVertexAssociation/src/JetVertexMain.cc

0001
0002 #include "JetMETCorrections/JetVertexAssociation/interface/JetVertexMain.h"
0003 #include "DataFormats/JetReco/interface/CaloJet.h"
0004 #include "DataFormats/TrackReco/interface/Track.h"
0005 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include <cmath>
0008 #include <string>
0009 using namespace reco;
0010 using namespace edm;
0011
0012 JetVertexMain::JetVertexMain(const ParameterSet& parameters) {
0013   cutSigmaZ = parameters.getParameter<double>("JV_sigmaZ");
0014   cutDeltaZ = parameters.getParameter<double>("JV_deltaZ");
0015   threshold = parameters.getParameter<double>("JV_alpha_threshold");
0016   cone_size = parameters.getParameter<double>("JV_cone_size");
0017   Algo = parameters.getParameter<int>("JV_type_Algo");
0018   cutType = parameters.getParameter<std::string>("JV_cutType");
0019 }
0020
0021 std::pair<double, bool> JetVertexMain::Main(const reco::CaloJet& jet,
0022                                             edm::Handle<TrackCollection> tracks,
0023                                             double signal_vert_Z,
0024                                             double signal_vert_z_error) const {
0025   std::pair<double, bool> parameter;
0026
0027   double jet_et = jet.et();
0028   double jet_phi = jet.phi();
0029   double jet_eta = jet.eta();
0030
0031   //  cout<<"JET: "<<jet_et<<endl;
0032   double Pt_jets_X = 0.;
0033   double Pt_jets_Y = 0.;
0034   double Pt_jets_X_tot = 0.;
0035   double Pt_jets_Y_tot = 0.;
0036
0037   TrackCollection::const_iterator track = tracks->begin();
0038
0039   if (!tracks->empty()) {
0040     for (; track != tracks->end(); track++) {
0041       double Vertex_Z = track->vz();
0042       double Vertex_Z_Error = track->dzError();
0043       double track_eta = track->eta();
0044       double track_phi = track->phi();
0045
0046       if (DeltaR(track_eta, jet_eta, track_phi, jet_phi) < cone_size) {
0047         double DeltaZ = Vertex_Z - signal_vert_Z;
0048         double DeltaZ_Error = sqrt((Vertex_Z_Error * Vertex_Z_Error) + (signal_vert_z_error * signal_vert_z_error));
0049         Pt_jets_X_tot += track->px();
0050         Pt_jets_Y_tot += track->py();
0051         bool discriminator;
0052         if (cutType == "sig")
0053           discriminator = (fabs(DeltaZ) / DeltaZ_Error) <= cutSigmaZ;
0054         else
0055           discriminator = fabs(DeltaZ) < cutDeltaZ;
0056
0057         if (discriminator) {
0058           Pt_jets_X += track->px();
0059           Pt_jets_Y += track->py();
0060         }
0061       }
0062     }
0063   }
0064   double Var = -1;
0065
0066   if (Algo == 1)
0067     Var = Track_Pt(Pt_jets_X, Pt_jets_Y) / jet_et;
0068   else if (Algo == 2) {
0069     if (Track_Pt(Pt_jets_X_tot, Pt_jets_Y_tot) != 0)
0070       Var = Track_Pt(Pt_jets_X, Pt_jets_Y) / Track_Pt(Pt_jets_X_tot, Pt_jets_Y_tot);
0071     else
0072       std::cout << "[Jets] JetVertexAssociation: Warning! problems for  Algo = 2: possible division by zero .."
0073                 << std::endl;
0074   } else {
0075     Var = Track_Pt(Pt_jets_X, Pt_jets_Y) / jet_et;
0076     std::cout << "[Jets] JetVertexAssociation: Warning! Algo = " << Algo << " not found; using Algo = 1" << std::endl;
0077   }
0078
0079   //  cout<<"Var = "<<Var<<endl;
0080
0081   if (Var >= threshold)
0082     parameter = std::pair<double, bool>(Var, true);
0083   else
0084     parameter = std::pair<double, bool>(Var, false);
0085
0086   return parameter;
0087 }
0088
0089 double JetVertexMain::DeltaR(double eta1, double eta2, double phi1, double phi2) const {
0090   double dphi = fabs(phi1 - phi2);
0091   if (dphi > M_PI)
0092     dphi = 2 * M_PI - dphi;
0093   double deta = fabs(eta1 - eta2);
0094   return sqrt(dphi * dphi + deta * deta);
0095 }
0096
0097 double JetVertexMain::Track_Pt(double px, double py) const { return sqrt(px * px + py * py); }