Major JME developments for 73x

CommonTools/CandAlgos/plugins/CandViewRefMerger.cc

@@ -26,9 +26,8 @@ class CandViewRefMerger : public edm::EDProducer {
     for(std::vector<edm::EDGetTokenT<reco::CandidateView> >::const_iterator i = srcTokens_.begin(); i != srcTokens_.end(); ++i) {
       edm::Handle<reco::CandidateView> src;
       evt.getByToken(*i, src);
-      reco::CandidateBaseRefVector refs = src->refVector();
-      for(reco::CandidateBaseRefVector::const_iterator j = refs.begin(); j != refs.end(); ++j)
-	out->push_back(*j);
+      for(size_t j = 0; j < src->size(); ++j)
+	out->push_back(src->refAt(j));
     }
     evt.put(out);
   }

CommonTools/PileupModules/BuildFile.xml

@@ -0,0 +1,10 @@
+<use   name="DataFormats/Common"/>
+<use   name="DataFormats/ParticleFlowCandidate"/>
+<use   name="FWCore/ParameterSet"/>
+<use   name="root"/>
+<use   name="rootrflx"/>
+<use   name="rootcore"/>
+<use   name="fastjet"/>
+<export>
+  <lib   name="1"/>
+</export>

CommonTools/PileupModules/interface/PuppiAlgo.h

@@ -0,0 +1,52 @@
+#ifndef CommonTools_PileupModules_PuppiAlgo_h
+#define CommonTools_PileupModules_PuppiAlgo_h
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "fastjet/PseudoJet.hh"
+#include <vector>
+
+class PuppiAlgo{ 
+public:
+  PuppiAlgo(edm::ParameterSet &iConfig);
+  ~PuppiAlgo();
+  //Computing Mean and RMS
+  void   reset();
+  void   add(const fastjet::PseudoJet &iParticle,const double &iVal,const unsigned int iAlgo);
+  void   computeMedRMS(const unsigned int &iAlgo,const double &iPVFrac);
+  //Get the Weight
+  double compute(std::vector<double> const &iVals,double iChi2) const;
+  //Helpers
+  inline double ptMin() const { return fPtMin; }
+  inline double etaMin() const { return fEtaMin; }
+  inline double etaMax() const { return fEtaMax; }
+  inline int    numAlgos () const { return fNAlgos;}
+  inline int    algoId  ( unsigned int iAlgo) const { return fAlgoId.at(iAlgo); }
+  inline bool   isCharged  ( unsigned int iAlgo) const { return fCharged.at(iAlgo); }
+  inline double coneSize  ( unsigned int iAlgo) const { return fConeSize.at(iAlgo); }
+  inline double neutralPt  (int iNPV) const { return fNeutralPtMin + iNPV * fNeutralPtSlope; }
+
+private:  
+  unsigned int   fNAlgos;
+  float  fEtaMax;
+  float  fEtaMin;
+  float  fPtMin ;
+  double fNeutralPtMin;
+  double fNeutralPtSlope;
+  std::vector<float>  fPups;
+  std::vector<float>  fPupsPV;
+  std::vector<int>    fAlgoId;
+  std::vector<bool>   fCharged;
+  std::vector<bool>   fAdjust;
+  std::vector<int>    fCombId;
+  std::vector<double> fConeSize;
+  std::vector<double> fRMSPtMin;
+  std::vector<double> fRMSScaleFactor;
+  std::vector<double> fRMS;
+  std::vector<double> fMedian;
+  std::vector<double> fMean;
+  std::vector<int>    fNCount;
+};
+
+
+#endif

CommonTools/PileupModules/interface/PuppiContainer.h

@@ -0,0 +1,69 @@
+#ifndef COMMONTOOLS_PUPPI_PUPPICONTAINER_H_
+#define COMMONTOOLS_PUPPI_PUPPICONTAINER_H_
+
+#include "CommonTools/PileupModules/interface/PuppiAlgo.h"
+#include "CommonTools/PileupModules/interface/RecoObj.h"
+#include "fastjet/internal/base.hh"
+#include "fastjet/PseudoJet.hh"
+
+//FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
+
+
+//......................
+class PuppiContainer{
+public:
+
+
+  // Helper class designed to store Puppi information inside of fastjet pseudojets.
+  // In CMSSW we use the user_index to refer to the index of the input collection, 
+  // but Puppi uses it to decide between NHs, PV CHs, and PU CHs. Instead,
+  // make that a register. 
+  class PuppiUserInfo : public fastjet::PseudoJet::UserInfoBase {
+   public : 
+     PuppiUserInfo( int puppi_register = -1) : puppi_register_(puppi_register) {}
+     virtual ~PuppiUserInfo(){}
+
+     void set_puppi_register(int i) { puppi_register_ = i; }
+
+     inline int puppi_register() const { return puppi_register_; }
+
+   protected : 
+     int puppi_register_;     /// Used by puppi algorithm to decide neutrals vs PV vs PU
+  };
+
+
+
+
+    PuppiContainer(const edm::ParameterSet &iConfig);
+    ~PuppiContainer(); 
+    void initialize(const std::vector<RecoObj> &iRecoObjects);
+    std::vector<fastjet::PseudoJet> const & pfParticles() const { return fPFParticles; }    
+    std::vector<fastjet::PseudoJet> const & pvParticles() const { return fChargedPV; }        
+    std::vector<double> const & puppiWeights() ;
+    std::vector<fastjet::PseudoJet> const & puppiParticles() const { return fPupParticles;}
+
+protected:
+    double  goodVar      (fastjet::PseudoJet const &iPart,std::vector<fastjet::PseudoJet> const &iParts, int iOpt,double iRCone);
+    void    getRMSAvg    (int iOpt,std::vector<fastjet::PseudoJet> const &iConstits,std::vector<fastjet::PseudoJet> const &iParticles,std::vector<fastjet::PseudoJet> const &iChargeParticles);
+    double  getChi2FromdZ(double iDZ);
+    int     getPuppiId   ( float iPt, float iEta);
+    double  var_within_R (int iId, const std::vector<fastjet::PseudoJet> & particles, const fastjet::PseudoJet& centre, double R);  
+
+    std::vector<RecoObj>   fRecoParticles;
+    std::vector<fastjet::PseudoJet> fPFParticles;
+    std::vector<fastjet::PseudoJet> fChargedPV;
+    std::vector<fastjet::PseudoJet> fPupParticles;
+    std::vector<double>    fWeights;
+    std::vector<double>    fVals;
+    bool   fApplyCHS;
+    bool   fUseExp;
+    double fNeutralMinPt;
+    double fNeutralSlope;
+    double fPuppiWeightCut;
+    int    fNAlgos;
+    int    fNPV;
+    double fPVFrac;
+    std::vector<PuppiAlgo> fPuppiAlgo;
+};
+#endif
+

CommonTools/PileupModules/interface/RecoObj.h

@@ -0,0 +1,31 @@
+#ifndef CommonTools_PileupModules_PUPPI_RECOOBJ_HH
+#define CommonTools_PileupModules_PUPPI_RECOOBJ_HH
+
+class RecoObj 
+{
+public:
+      RecoObj():
+	pt(0), eta(0), phi(0), m(0),
+	id(0),pfType(-1),vtxId(-1),
+	trkChi2(0),vtxChi2(0),
+	time(0),depth(0),
+	expProb(0),expChi2PU(0),expChi2(0),
+	dZ(0),d0(0),charge(0)
+    {}
+    ~RecoObj(){}
+
+    float         pt, eta, phi, m;  // kinematics
+    int           id;
+    int           pfType;
+    int           vtxId;               // Vertex Id from Vertex Collection
+    float         trkChi2;             // Track Chi2
+    float         vtxChi2;             // Vertex Chi2
+    float         time,depth;    // Usefule Info
+    float         expProb;
+    float         expChi2PU;
+    float         expChi2;
+    float         dZ;
+    float         d0;
+    int           charge;
+};
+#endif

CommonTools/PileupModules/plugins/BuildFile.xml

@@ -0,0 +1,15 @@
+<export>
+</export>
+<library   name="CommonToolsPileupModules_plugins" file="*.cc">
+  <use   name="DataFormats/ParticleFlowCandidate"/>
+  <use   name="DataFormats/JetReco"/>
+  <use   name="FWCore/Framework"/>
+  <use   name="FWCore/MessageLogger"/>
+  <use   name="FWCore/ParameterSet"/>
+  <use   name="JetMETCorrections/Objects"/>
+  <use   name="CommonTools/Utils"/>
+  <use   name="CommonTools/ParticleFlow"/>
+  <use   name="CommonTools/PileupModules"/>
+  <use   name="fastjet"/>
+  <use   name="fastjet-contrib"/>
+</library>

CommonTools/PileupModules/plugins/PuppiProducer.cc

@@ -0,0 +1,202 @@
+// 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 "DataFormats/Common/interface/ValueMap.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "DataFormats/Common/interface/View.h"
+#include "DataFormats/Candidate/interface/Candidate.h"
+#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
+#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
+#include "DataFormats/VertexReco/interface/VertexFwd.h"
+#include "DataFormats/VertexReco/interface/Vertex.h"
+#include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
+#include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
+//Main File
+#include "fastjet/PseudoJet.hh"
+#include "CommonTools/PileupModules/plugins/PuppiProducer.h"
+
+
+
+// ------------------------------------------------------------------------------------------
+PuppiProducer::PuppiProducer(const edm::ParameterSet& iConfig) {
+  fUseDZ     = iConfig.getParameter<bool>("UseDeltaZCut");
+  fDZCut     = iConfig.getParameter<double>("DeltaZCut");
+  fPuppiContainer = std::unique_ptr<PuppiContainer> ( new PuppiContainer(iConfig) );
+
+  tokenPFCandidates_
+    = consumes<CandidateView>(iConfig.getParameter<edm::InputTag>("candName"));
+  tokenVertices_
+    = consumes<VertexCollection>(iConfig.getParameter<edm::InputTag>("vertexName"));
+
+
+  produces<edm::ValueMap<float> > ("PuppiWeights");
+  produces<edm::ValueMap<LorentzVector> > ("PuppiP4s");
+  produces<PFOutputCollection>();
+
+
+}
+// ------------------------------------------------------------------------------------------
+PuppiProducer::~PuppiProducer(){
+}
+// ------------------------------------------------------------------------------------------
+void PuppiProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+
+  // Get PFCandidate Collection
+  edm::Handle<CandidateView> hPFProduct;
+  iEvent.getByToken(tokenPFCandidates_,hPFProduct);
+  const CandidateView *pfCol = hPFProduct.product();
+
+  // Get vertex collection w/PV as the first entry?
+  edm::Handle<reco::VertexCollection> hVertexProduct;
+  iEvent.getByToken(tokenVertices_,hVertexProduct);
+  const reco::VertexCollection *pvCol = hVertexProduct.product();
+
+  //Fill the reco objects
+  fRecoObjCollection.clear();
+  for(CandidateView::const_iterator itPF = pfCol->begin(); itPF!=pfCol->end(); itPF++) {
+    RecoObj pReco;
+    pReco.pt  = itPF->pt();
+    pReco.eta = itPF->eta();
+    pReco.phi = itPF->phi();
+    pReco.m   = itPF->mass();
+    pReco.charge = itPF->charge(); 
+    const reco::Vertex *closestVtx = 0;
+    double pDZ    = -9999; 
+    double pD0    = -9999; 
+    int    pVtxId = -9999; 
+    bool lFirst = true;
+    const pat::PackedCandidate *lPack = dynamic_cast<const pat::PackedCandidate*>(&(*itPF));
+    if(lPack == 0 ) { 
+      const reco::PFCandidate *pPF = dynamic_cast<const reco::PFCandidate*>(&(*itPF));
+      for(reco::VertexCollection::const_iterator iV = pvCol->begin(); iV!=pvCol->end(); ++iV) {
+	if(lFirst) { 
+	  if      ( pPF->trackRef().isNonnull()    ) pDZ = pPF->trackRef()   ->dz(iV->position());
+	  else if ( pPF->gsfTrackRef().isNonnull() ) pDZ = pPF->gsfTrackRef()->dz(iV->position());
+	  if      ( pPF->trackRef().isNonnull()    ) pD0 = pPF->trackRef()   ->d0();
+	  else if ( pPF->gsfTrackRef().isNonnull() ) pD0 = pPF->gsfTrackRef()->d0();
+	  lFirst = false;
+	  if(pDZ > -9999) pVtxId = 0; 
+	}
+	if(iV->trackWeight(pPF->trackRef())>0) {
+	  closestVtx  = &(*iV);
+	  break;
+	}
+	pVtxId++;
+      }
+    } else if(lPack->vertexRef().isNonnull() )  {
+      pDZ        = lPack->dz(); 
+      pD0        = lPack->dxy(); 
+      closestVtx = &(*(lPack->vertexRef()));
+      pVtxId     = (lPack->fromPV() !=  (pat::PackedCandidate::PVUsedInFit));
+      if(lPack->fromPV() ==  (pat::PackedCandidate::PVLoose || pat::PackedCandidate::PVTight)) closestVtx = 0;
+    }
+    pReco.dZ      = pDZ;
+    pReco.d0      = pD0;
+
+    if(closestVtx == 0) pReco.vtxId = -1;
+    if(closestVtx != 0) pReco.vtxId = pVtxId;
+    //if(closestVtx != 0) pReco.vtxChi2 = closestVtx->trackWeight(itPF->trackRef());
+    //Set the id for Puppi Algo: 0 is neutral pfCandidate, id = 1 for particles coming from PV and id = 2 for charged particles from non-leading vertex
+    pReco.id       = 0; 
+
+    if(closestVtx != 0 && pVtxId == 0 && fabs(pReco.charge) > 0) pReco.id = 1;
+    if(closestVtx != 0 && pVtxId >  0 && fabs(pReco.charge) > 0) pReco.id = 2;
+    //Add a dZ cut if wanted (this helps)
+    if(fUseDZ && pDZ > -9999 && closestVtx == 0 && (fabs(pDZ) < fDZCut)) pReco.id = 1; 
+    if(fUseDZ && pDZ > -9999 && closestVtx == 0 && (fabs(pDZ) > fDZCut)) pReco.id = 2; 
+
+    //std::cout << "pVtxId = " << pVtxId << ", and charge = " << itPF->charge() << ", and closestVtx = " << closestVtx << ", and id = " << pReco.id << std::endl;
+
+    fRecoObjCollection.push_back(pReco);
+  }
+  fPuppiContainer->initialize(fRecoObjCollection);
+
+  //Compute the weights
+  const std::vector<double> lWeights = fPuppiContainer->puppiWeights();
+  //Fill it into the event
+  std::auto_ptr<edm::ValueMap<float> > lPupOut(new edm::ValueMap<float>());
+  edm::ValueMap<float>::Filler  lPupFiller(*lPupOut);
+  lPupFiller.insert(hPFProduct,lWeights.begin(),lWeights.end());
+  lPupFiller.fill();
+
+  //Fill a new PF Candidate Collection and write out the ValueMap of the new p4s.
+  // Since the size of the ValueMap must be equal to the input collection, we need
+  // to search the "puppi" particles to find a match for each input. If none is found,
+  // the input is set to have a four-vector of 0,0,0,0
+  const std::vector<fastjet::PseudoJet> lCandidates = fPuppiContainer->puppiParticles();
+  fPuppiCandidates.reset( new PFOutputCollection );
+  std::auto_ptr<edm::ValueMap<LorentzVector> > p4PupOut(new edm::ValueMap<LorentzVector>());
+  LorentzVectorCollection puppiP4s;
+  for ( auto i0 = hPFProduct->begin(),
+	  i0begin = hPFProduct->begin(),
+	  i0end = hPFProduct->end(); i0 != i0end; ++i0 ) {
+    //for(unsigned int i0 = 0; i0 < lCandidates.size(); i0++) {
+    //reco::PFCandidate pCand;
+    reco::PFCandidate const * tmpCandPF = dynamic_cast<reco::PFCandidate const *>( &*i0 );
+    if ( tmpCandPF == 0 ) {
+      throw cms::Exception("LogicError") << "Non-PFCandidates treated as PFCandidates in PuppiProducer." << std::endl;
+    }
+    auto id = (tmpCandPF != 0) ? tmpCandPF->translatePdgIdToType(i0->pdgId()) : reco::PFCandidate::X;
+    reco::PFCandidate pCand( i0->charge(),
+			     i0->p4(),
+			     id );
+    LorentzVector pVec = i0->p4();
+    int val = i0 - i0begin;
+
+    // Find the Puppi particle matched to the input collection using the "user_index" of the object. 
+    auto puppiMatched = find_if( lCandidates.begin(), lCandidates.end(), [&val]( fastjet::PseudoJet const & i ){ return i.user_index() == val; } );
+    if ( puppiMatched != lCandidates.end() ) {
+      pVec.SetPxPyPzE(puppiMatched->px(),puppiMatched->py(),puppiMatched->pz(),puppiMatched->E());
+    } else {
+      pVec.SetPxPyPzE( 0, 0, 0, 0);
+    }
+    pCand.setP4(pVec);
+    puppiP4s.push_back( pVec );
+    fPuppiCandidates->push_back(pCand);
+  }
+
+  //Compute the modified p4s
+  edm::ValueMap<LorentzVector>::Filler  p4PupFiller(*p4PupOut);
+  p4PupFiller.insert(hPFProduct,puppiP4s.begin(), puppiP4s.end() );
+  p4PupFiller.fill();
+
+  iEvent.put(lPupOut,"PuppiWeights");
+  iEvent.put(p4PupOut,"PuppiP4s");
+  iEvent.put(fPuppiCandidates);
+}
+
+// ------------------------------------------------------------------------------------------
+void PuppiProducer::beginJob() {
+}
+// ------------------------------------------------------------------------------------------
+void PuppiProducer::endJob() {
+}
+// ------------------------------------------------------------------------------------------
+void PuppiProducer::beginRun(edm::Run&, edm::EventSetup const&) {
+}
+// ------------------------------------------------------------------------------------------
+void PuppiProducer::endRun(edm::Run&, edm::EventSetup const&) {
+}
+// ------------------------------------------------------------------------------------------
+void PuppiProducer::beginLuminosityBlock(edm::LuminosityBlock&, edm::EventSetup const&) {
+}
+// ------------------------------------------------------------------------------------------
+void PuppiProducer::endLuminosityBlock(edm::LuminosityBlock&, edm::EventSetup const&) {
+}
+// ------------------------------------------------------------------------------------------
+void PuppiProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  //The following says we do not know what parameters are allowed so do no validation
+  // Please change this to state exactly what you do use, even if it is no parameters
+  edm::ParameterSetDescription desc;
+  desc.setUnknown();
+  descriptions.addDefault(desc);
+}
+//define this as a plug-in
+DEFINE_FWK_MODULE(PuppiProducer);