Fixed new fast pv3

HLTrigger/JetMET/BuildFile.xml

@@ -14,4 +14,6 @@
 <use   name="RecoMET/METAlgorithms"/>
 <use   name="HLTrigger/HLTcore"/>
 <use   name="RecoJets/JetProducers"/>
+<use   name="TrackingTools/IPTools"/>
+<use   name="TrackingTools/TransientTrack"/>
 <flags   EDM_PLUGIN="1"/>

HLTrigger/JetMET/src/PixelJetPuId.cc

@@ -0,0 +1,217 @@
+// -*- C++ -*-
+//
+// Package:    PixelJetPuId
+// Class:      PixelJetPuId
+// 
+/**\class PixelJetPuId PixelJetPuId.cc RecoBTag/PixelJetPuId/src/PixelJetPuId.cc
+
+ Description:
+ The PixelJetPuId module select all the pixel tracks compatible with a jet.
+ If the sum of the tracks momentum is under a threshold the jet is tagged as "PUjets".
+
+ Implementation:
+     [Notes on implementation]
+*/
+//
+// Original Author:  Silvio DONATO
+//         Created:  Wed Dec 18 10:05:40 CET 2013
+//
+//
+
+
+// system include files
+#include <memory>
+
+// user include files
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/EDProducer.h"
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+#include "FWCore/Utilities/interface/InputTag.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+
+#include "DataFormats/TrackReco/interface/Track.h"
+#include "DataFormats/JetReco/interface/CaloJet.h"
+#include "DataFormats/JetReco/interface/CaloJetCollection.h"
+#include "DataFormats/VertexReco/interface/VertexFwd.h"
+#include "DataFormats/VertexReco/interface/Vertex.h"
+
+#include "TrackingTools/IPTools/interface/IPTools.h"
+#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
+#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
+#include "TrackingTools/Records/interface/TransientTrackRecord.h"
+
+#include "DataFormats/Math/interface/deltaR.h"
+
+//
+// class declaration
+//
+
+class PixelJetPuId : public edm::EDProducer {
+   public:
+      PixelJetPuId(const edm::ParameterSet&);
+      virtual ~PixelJetPuId();
+
+      static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+   private:
+      virtual void produce(edm::Event&, const edm::EventSetup&);
+
+
+      // ----------member data ---------------------------
+     edm::InputTag m_primaryVertex;
+     edm::InputTag m_tracks;
+     edm::InputTag m_jets;
+     edm::EDGetTokenT<std::vector<reco::Track> > tracksToken;
+     edm::EDGetTokenT<edm::View<reco::CaloJet> > jetsToken;
+     edm::EDGetTokenT<reco::VertexCollection> primaryVertexToken;
+
+     double m_MinTrackPt; 
+     double m_MaxTrackChi2; 
+     double m_MaxTrackDistanceToJet; 
+
+     bool   m_fwjets;
+     double m_mineta_fwjets;
+     double m_minet_fwjets;
+
+     double m_MinGoodJetTrackPt;
+     double m_MinGoodJetTrackPtRatio; 
+};
+
+
+//
+// constructors and destructor
+//
+PixelJetPuId::PixelJetPuId(const edm::ParameterSet& iConfig)
+{
+  //InputTag
+  m_tracks           = iConfig.getParameter<edm::InputTag>("tracks");
+  tracksToken        = consumes<std::vector<reco::Track> >(m_tracks);
+  m_jets             = iConfig.getParameter<edm::InputTag>("jets");
+  jetsToken          = consumes<edm::View<reco::CaloJet> >(m_jets);
+  m_primaryVertex    = iConfig.getParameter<edm::InputTag>("primaryVertex");
+  primaryVertexToken = consumes<reco::VertexCollection>(m_primaryVertex);
+
+  //Tracks Selection
+  m_MinTrackPt             = iConfig.getParameter<double>("MinTrackPt");
+  m_MaxTrackDistanceToJet  = iConfig.getParameter<double>("MaxTrackDistanceToJet");
+  m_MaxTrackChi2           = iConfig.getParameter<double>("MaxTrackChi2");
+
+  //A jet is defined as a signal jet if Sum(trackPt) > minPt or Sum(comp.trackPt)/CaloJetPt > minPtRatio
+  m_MinGoodJetTrackPt      = iConfig.getParameter<double>("MinGoodJetTrackPt");
+  m_MinGoodJetTrackPtRatio = iConfig.getParameter<double>("MinGoodJetTrackPtRatio");
+
+  m_fwjets                 = iConfig.getParameter<bool>("UseForwardJetsAsNoPU");
+  m_mineta_fwjets          = iConfig.getParameter<double>("MinEtaForwardJets");
+  m_minet_fwjets           = iConfig.getParameter<double>("MinEtForwardJets");
+
+  produces<std::vector<reco::CaloJet> >(); 
+  produces<std::vector<reco::CaloJet> >("PUjets"); 
+}
+
+
+PixelJetPuId::~PixelJetPuId() {}
+
+
+void
+PixelJetPuId::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<edm::InputTag> ("jets",edm::InputTag("hltCaloJetL1FastJetCorrected"));
+  desc.add<edm::InputTag> ("tracks",edm::InputTag("hltPixelTracksNoPU"));
+  desc.add<edm::InputTag> ("primaryVertex",edm::InputTag("hltFastPVPixelVertices"));
+  desc.add<double>("MinGoodJetTrackPtRatio",0.045);
+  desc.add<double>("MinGoodJetTrackPt",1.8);
+  desc.add<double>("MaxTrackDistanceToJet",0.04);
+  desc.add<double>("MinTrackPt",0.6);
+  desc.add<double>("MaxTrackChi2",20.);
+  desc.add<bool>("UseForwardJetsAsNoPU",true);
+  desc.add<double>("MinEtaForwardJets",2.4);
+  desc.add<double>("MinEtForwardJets",40.);
+  descriptions.add("pixelJetPuId",desc);
+}
+
+//
+// member functions
+//
+
+// ------------ method called on each new Event  ------------
+void PixelJetPuId::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
+{
+  using namespace edm;
+  std::auto_ptr<std::vector<reco::CaloJet> > pOut(new std::vector<reco::CaloJet> );
+  std::auto_ptr<std::vector<reco::CaloJet> > pOut_PUjets(new std::vector<reco::CaloJet> );
+
+   //get tracks
+  Handle<std::vector<reco::Track> > tracks;
+  iEvent.getByToken(tracksToken, tracks);
+
+  //get jets
+  Handle<edm::View<reco::CaloJet> > jets;
+  iEvent.getByToken(jetsToken, jets);
+
+  //get primary vertices
+  Handle<reco::VertexCollection> primaryVertex;
+  iEvent.getByToken(primaryVertexToken, primaryVertex);
+
+  //get Transient Track Builder
+  edm::ESHandle<TransientTrackBuilder> builder;
+  iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder", builder);
+
+  //loop on trackIPTagInfos
+  if(primaryVertex->size()>0)
+    {
+      const reco::Vertex* pv = &*primaryVertex->begin();
+      //loop on jets
+      for(edm::View<reco::CaloJet>::const_iterator itJet = jets->begin();  itJet != jets->end(); itJet++ ) {
+
+	math::XYZVector jetMomentum = itJet->momentum();
+	GlobalVector direction(jetMomentum.x(), jetMomentum.y(), jetMomentum.z());
+
+	math::XYZVector trMomentum;
+
+	//loop on tracks
+	if(fabs(itJet->eta())>m_mineta_fwjets)
+	  {
+	    if((m_fwjets) && (itJet->et()>m_minet_fwjets))
+	      pOut->push_back(*itJet);// fill forward jet as signal jet
+	  }
+	else 
+	  {
+	    for(std::vector<reco::Track>::const_iterator itTrack = tracks->begin(); itTrack != tracks->end(); ++itTrack) 
+	      {
+		double deltaR2 = reco::deltaR2( jetMomentum, *itTrack );
+		if(deltaR2<0.25)
+		  {
+		    reco::TransientTrack transientTrack = builder->build(*itTrack);
+		    float jetTrackDistance = -((IPTools::jetTrackDistance(transientTrack, direction, *pv)).second).value();
+
+		    //select the tracks compabible with the jet
+		    if(( itTrack->pt() > m_MinTrackPt) && ( itTrack->normalizedChi2() < m_MaxTrackChi2) && (jetTrackDistance<m_MaxTrackDistanceToJet))
+		      {
+			trMomentum += itTrack->momentum(); //calculate the Sum(trackPt)
+		      }
+		  }
+	      }
+	    //if Sum(comp.trackPt)/CaloJetPt > minPtRatio or Sum(trackPt) > minPt  the jet is a signal jet
+	    if(trMomentum.rho()/jetMomentum.rho() > m_MinGoodJetTrackPtRatio || trMomentum.rho() > m_MinGoodJetTrackPt ) 
+	      {
+		pOut->push_back(*itJet);        // fill it as signal jet
+	      }
+	    else//else it is a PUjet
+	      {
+		pOut_PUjets->push_back(*itJet); // fill it as PUjets
+	      }
+	  }
+      }
+    }
+  iEvent.put(pOut);
+  iEvent.put(pOut_PUjets,"PUjets");
+
+}
+
+//define this as a plug-in
+DEFINE_FWK_MODULE(PixelJetPuId);

RecoPixelVertexing/PixelVertexFinding/BuildFile.xml

@@ -1,18 +1,19 @@
-<flags   EDM_PLUGIN="1"/>
-<use   name="DataFormats/VertexReco"/>
-<use   name="FWCore/Framework"/>
-<use   name="FWCore/ParameterSet"/>
-<use   name="FWCore/PluginManager"/>
-<use   name="Geometry/Records"/>
-<use   name="Geometry/TrackerGeometryBuilder"/>
-<use   name="DataFormats/GeometryCommonDetAlgo"/>
-<use   name="CommonTools/Clustering1D"/>
-<use   name="DataFormats/TrackReco"/>
-<use   name="RecoPixelVertexing/PixelTriplets"/>
-<use   name="RecoPixelVertexing/PixelTrackFitting"/>
+<flags EDM_PLUGIN="1"/>
+<use name="DataFormats/VertexReco"/>
+<use name="FWCore/Framework"/>
+<use name="FWCore/ParameterSet"/>
+<use name="FWCore/PluginManager"/>
+<use name="Geometry/Records"/>
+<use name="Geometry/TrackerGeometryBuilder"/>
+<use name="DataFormats/GeometryCommonDetAlgo"/>
+<use name="CommonTools/Clustering1D"/>
+<use name="DataFormats/TrackReco"/>
+<use name="RecoPixelVertexing/PixelTriplets"/>
+<use name="RecoPixelVertexing/PixelTrackFitting"/>
 
 <use name="DataFormats/SiPixelCluster"/>
 <use name="DataFormats/JetReco"/>
 <use name="RecoLocalTracker/Records"/>
 <use name="RecoLocalTracker/ClusterParameterEstimator"/>
 <use name="DataFormats/TrackerRecHit2D"/>
+

RecoPixelVertexing/PixelVertexFinding/interface/FindPeakFastPV.h

@@ -0,0 +1,48 @@
+#ifndef RecoPixelVertexing_FindPeakFastPV_h
+#define RecoPixelVertexing_FindPeakFastPV_h
+/** \class FindPeakFastPV FindPeakFastPV.h RecoPixelVertexing/PixelVertexFinding/FindPeakFastPV.h 
+ * Given *zProjections* and *zWeights* find the peak of width *m_zClusterWidth*. 
+ * Use only values with *zWeights*>*m_weightCut*. 
+ * Look near *oldVertex* within *m_zClusterSearchArea*.
+ *
+ *  \author Silvio Donato (SNS)
+ */
+
+
+float FindPeakFastPV(const std::vector<float> &zProjections, const std::vector<float> &zWeights, const float oldVertex, const float m_zClusterWidth, const float m_zClusterSearchArea, const float m_weightCut) {
+float centerWMax = oldVertex;
+if( m_zClusterWidth > 0 && m_zClusterSearchArea >0 ) 
+{
+  std::vector<float>::const_iterator itCenter = zProjections.begin();
+  std::vector<float>::const_iterator itLeftSide = zProjections.begin();
+  std::vector<float>::const_iterator itRightSide = zProjections.begin();
+  const float zClusterWidth = m_zClusterWidth*0.5; //take half zCluster width 
+  const float zLowerBound = oldVertex - zClusterWidth - m_zClusterSearchArea;
+  const float zUpperBound = oldVertex + zClusterWidth + m_zClusterSearchArea;
+  float maxW=0;
+  for(;itCenter!=zProjections.end(); itCenter++)  {
+    //loop only on the zProjections within oldVertex +/- (zClusterWidth + m_zClusterSearchArea)
+    if( (*itCenter < zLowerBound) || (*itCenter > zUpperBound) ) continue; 
+
+    while(itLeftSide != zProjections.end() && (*itCenter - *itLeftSide) > zClusterWidth  ) itLeftSide++;
+    while(itRightSide != zProjections.end() && (*itRightSide - *itCenter) < zClusterWidth  ) itRightSide++;
+    float nWeighted = 0;
+    float centerW= 0;
+
+    for(std::vector<float>::const_iterator ii = itLeftSide; ii != itRightSide; ii++) {
+      //loop inside the peak and calculate its weight
+      if(zWeights[ii-zProjections.begin()]<m_weightCut) continue;
+      nWeighted+=zWeights[ii-zProjections.begin()]; 
+      centerW+=(*ii)*zWeights[ii-zProjections.begin()]; 
+    }
+    centerW/=nWeighted; //calculate the weighted peak center
+    if(nWeighted > maxW) {
+      maxW=nWeighted;
+      centerWMax=centerW;	       
+    }
+  }
+}
+return centerWMax;
+}
+
+#endif