JetAssociationAlgorithms/src/JetTracksAssociationDRVertexAssigned.cc

0001 // Associate jets with tracks by simple "dR" criteria

0002 // Fedor Ratnikov (UMd), Aug. 28, 2007

0003
0004 #include "RecoJets/JetAssociationAlgorithms/interface/JetTracksAssociationDRVertexAssigned.h"
0005
0006 #include "DataFormats/JetReco/interface/Jet.h"
0007 #include "DataFormats/TrackReco/interface/Track.h"
0008
0009 #include "DataFormats/Math/interface/deltaR.h"
0010 #include "DataFormats/Math/interface/Vector3D.h"
0011 #include "DataFormats/Math/interface/Point3D.h"
0012 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0013
0014 JetTracksAssociationDRVertexAssigned::JetTracksAssociationDRVertexAssigned(double fDr) : mDeltaR2Threshold(fDr * fDr) {}
0015
0016 void JetTracksAssociationDRVertexAssigned::produce(reco::JetTracksAssociation::Container* fAssociation,
0017                                                    const std::vector<edm::RefToBase<reco::Jet> >& fJets,
0018                                                    const std::vector<reco::TrackRef>& fTracks,
0019                                                    const reco::VertexCollection& vertices) const {
0020   // cache tracks kinematics

0021   std::vector<math::RhoEtaPhiVector> trackP3s;
0022   std::map<int, double> trackvert;
0023
0024   // std::cout<<" Number of vertices "<<vertices.size()<<std::endl;

0025
0026   trackP3s.reserve(fTracks.size());
0027   for (unsigned i = 0; i < fTracks.size(); ++i) {
0028     const reco::Track* track = &*(fTracks[i]);
0029     trackP3s.push_back(math::RhoEtaPhiVector(track->p(), track->eta(), track->phi()));
0030
0031     // OK: Look for the tracks not associated with vertices

0032
0033     const reco::TrackBaseRef ttr1(fTracks[i]);
0034
0035     int trackhasvert = -1;
0036     for (reco::VertexCollection::const_iterator iv = vertices.begin(); iv != vertices.end(); iv++) {
0037       std::vector<reco::TrackBaseRef>::const_iterator rr = find((*iv).tracks_begin(), (*iv).tracks_end(), ttr1);
0038       if (rr != (*iv).tracks_end()) {
0039         trackhasvert = 1;
0040         trackvert[i] = (*iv).position().z();
0041         // std::cout<<" Z "<<i<<" "<<trackhasvert<<" "<<(*iv).position().z()<<std::endl;

0042         break;
0043       }
0044     }  // all vertices

0045     if (trackhasvert < 0) {
0046       // Take impact parameter of the track as vertex position

0047       math::XYZPoint ppt(0., 0., 0.);
0048       trackvert[i] = track->dz(ppt);
0049       //    std::cout<<" Z "<<i<<" "<<trackhasvert<<" "<<track->dz(ppt)<<" "<<track->vz()<<" "<<track->vx()<<" "<<

0050       //    track->vy()<<" "<<track->pz()<<std::endl;

0051     }
0052     // OK

0053   }  // tracks

0054
0055   for (unsigned j = 0; j < fJets.size(); ++j) {
0056     reco::TrackRefVector assoTracks;
0057     const reco::Jet* jet = &*(fJets[j]);
0058     double jetEta = jet->eta();
0059     double jetPhi = jet->phi();
0060     double neweta = 0;
0061     for (unsigned t = 0; t < fTracks.size(); ++t) {
0062       std::map<int, double>::iterator cur = trackvert.find(t);
0063       if (cur != trackvert.end()) {
0064         neweta = jet->physicsEta((*cur).second, jetEta);
0065       } else {
0066         neweta = jetEta;
0067         //std::cout<<" Lost track - not in map "<<std::endl;

0068       }
0069
0070       //std::cout<<" Old eta-new eta "<<(*cur).second<<" "<<jetEta<<" "<<neweta<<" Track "<<t<<" "<<trackP3s[t].eta()<<std::endl;

0071
0072       double dR2 = deltaR2(neweta, jetPhi, trackP3s[t].eta(), trackP3s[t].phi());
0073       if (dR2 < mDeltaR2Threshold)
0074         assoTracks.push_back(fTracks[t]);
0075     }
0076
0077     reco::JetTracksAssociation::setValue(fAssociation, fJets[j], assoTracks);
0078   }
0079 }