Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:21:35

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) {
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         if (cutType == "sig")
0052           discriminator = (fabs(DeltaZ) / DeltaZ_Error) <= cutSigmaZ;
0053         else
0054           discriminator = fabs(DeltaZ) < cutDeltaZ;
0055 
0056         if (discriminator) {
0057           Pt_jets_X += track->px();
0058           Pt_jets_Y += track->py();
0059         }
0060       }
0061     }
0062   }
0063   double Var = -1;
0064 
0065   if (Algo == 1)
0066     Var = Track_Pt(Pt_jets_X, Pt_jets_Y) / jet_et;
0067   else if (Algo == 2) {
0068     if (Track_Pt(Pt_jets_X_tot, Pt_jets_Y_tot) != 0)
0069       Var = Track_Pt(Pt_jets_X, Pt_jets_Y) / Track_Pt(Pt_jets_X_tot, Pt_jets_Y_tot);
0070     else
0071       std::cout << "[Jets] JetVertexAssociation: Warning! problems for  Algo = 2: possible division by zero .."
0072                 << std::endl;
0073   } else {
0074     Var = Track_Pt(Pt_jets_X, Pt_jets_Y) / jet_et;
0075     std::cout << "[Jets] JetVertexAssociation: Warning! Algo = " << Algo << " not found; using Algo = 1" << std::endl;
0076   }
0077 
0078   //  cout<<"Var = "<<Var<<endl;
0079 
0080   if (Var >= threshold)
0081     parameter = std::pair<double, bool>(Var, true);
0082   else
0083     parameter = std::pair<double, bool>(Var, false);
0084 
0085   return parameter;
0086 }
0087 
0088 double JetVertexMain::DeltaR(double eta1, double eta2, double phi1, double phi2) {
0089   double dphi = fabs(phi1 - phi2);
0090   if (dphi > M_PI)
0091     dphi = 2 * M_PI - dphi;
0092   double deta = fabs(eta1 - eta2);
0093   return sqrt(dphi * dphi + deta * deta);
0094 }
0095 
0096 double JetVertexMain::Track_Pt(double px, double py) { return sqrt(px * px + py * py); }