Merge pull request #13391 from ishvetso/PhotonIsolationCITK_80X_recipe

Isolation CITK 80X (rebased)

PhysicsTools/IsolationAlgos/plugins/CITKPFIsolationSumProducerForPUPPI.cc

@@ -0,0 +1,182 @@
+#include "FWCore/Framework/interface/global/EDProducer.h"
+#include "FWCore/Framework/interface/ConsumesCollector.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+#include "DataFormats/TrackReco/interface/TrackFwd.h"
+#include "DataFormats/MuonReco/interface/MuonFwd.h"
+#include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
+#include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
+#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
+#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
+#include "PhysicsTools/IsolationAlgos/interface/EventDependentAbsVeto.h"
+#include "DataFormats/Candidate/interface/CandidateFwd.h"
+#include "DataFormats/Candidate/interface/Candidate.h"
+#include "PhysicsTools/IsolationAlgos/interface/CITKIsolationConeDefinitionBase.h"
+#include "DataFormats/Common/interface/OwnVector.h"
+
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+
+#include <string>
+#include <unordered_map>
+
+//module to compute isolation sum weighted with PUPPI weights
+namespace citk {
+  class PFIsolationSumProducerForPUPPI : public edm::stream::EDProducer<> {
+
+  public:  
+    PFIsolationSumProducerForPUPPI(const edm::ParameterSet&);
+
+    virtual ~PFIsolationSumProducerForPUPPI() {}
+
+    virtual void beginLuminosityBlock(const edm::LuminosityBlock&,
+			      const edm::EventSetup&) override final;
+
+    virtual void produce(edm::Event&, const edm::EventSetup&) override final;
+
+  private:  
+    // datamembers
+    static constexpr unsigned kNPFTypes = 8;
+    typedef std::unordered_map<std::string,int> TypeMap;
+    typedef std::vector<std::unique_ptr<IsolationConeDefinitionBase> > IsoTypes;
+    typedef edm::View<reco::Candidate> CandView;
+    const TypeMap _typeMap;
+    edm::EDGetTokenT<CandView> _to_isolate, _isolate_with;
+    edm::EDGetTokenT<edm::ValueMap<float> > puppiValueMapToken_;//for puppiValueMap
+    edm::Handle<edm::ValueMap<float>> puppiValueMap;//puppiValueMap
+    // indexed by pf candidate type
+    std::array<IsoTypes,kNPFTypes> _isolation_types; 
+    std::array<std::vector<std::string>,kNPFTypes> _product_names;
+    bool useValueMapForPUPPI = true;
+  };
+}
+
+typedef citk::PFIsolationSumProducerForPUPPI CITKPFIsolationSumProducerForPUPPI;
+
+DEFINE_FWK_MODULE(CITKPFIsolationSumProducerForPUPPI);
+
+namespace citk {
+  PFIsolationSumProducerForPUPPI::PFIsolationSumProducerForPUPPI(const edm::ParameterSet& c) :
+    _typeMap( { {"h+",1},
+	        {"h0",5},
+		{"gamma",4},
+		{"electron",2},
+		{"muon",3},
+		{"HFh",6},
+		{"HFgamma",7} } ){
+    _to_isolate = 
+      consumes<CandView>(c.getParameter<edm::InputTag>("srcToIsolate"));
+    _isolate_with = 
+      consumes<CandView>(c.getParameter<edm::InputTag>("srcForIsolationCone"));
+      if (c.getParameter<edm::InputTag>("puppiValueMap").label().size() != 0) {
+        puppiValueMapToken_ = mayConsume<edm::ValueMap<float>>(c.getParameter<edm::InputTag>("puppiValueMap")); //getting token for puppiValueMap
+        useValueMapForPUPPI = true;
+      }
+      else useValueMapForPUPPI = false;
+    const std::vector<edm::ParameterSet>& isoDefs = 
+      c.getParameterSetVector("isolationConeDefinitions");
+    for( const auto& isodef : isoDefs ) {
+      const std::string& name = 
+	isodef.getParameter<std::string>("isolationAlgo");
+      const float coneSize = isodef.getParameter<double>("coneSize");
+      char buf[50];
+      std::sprintf(buf,"DR%.2f",coneSize);
+      std::string coneName(buf);
+      auto decimal = coneName.find('.');
+      if( decimal != std::string::npos ) coneName.erase(decimal,1);
+      const std::string& isotype = 
+	isodef.getParameter<std::string>("isolateAgainst");
+      IsolationConeDefinitionBase* theisolator =
+	CITKIsolationConeDefinitionFactory::get()->create(name,isodef);
+      theisolator->setConsumes(consumesCollector());
+      const auto thetype = _typeMap.find(isotype);
+      if( thetype == _typeMap.end() ) {
+	throw cms::Exception("InvalidIsolationType")
+	  << "Isolation type: " << isotype << " is not available in the "
+	  << "list of allowed isolations!.";
+      }
+      _isolation_types[thetype->second].emplace_back(theisolator);
+      const std::string dash("-");
+      std::string pname = isotype+dash+coneName+dash+theisolator->additionalCode();
+      _product_names[thetype->second].emplace_back(pname);
+      produces<edm::ValueMap<float> >(pname);
+    }
+  }
+
+  void  PFIsolationSumProducerForPUPPI::
+  beginLuminosityBlock(const edm::LuminosityBlock&,
+		       const edm::EventSetup& es) {
+    for( const auto& isolators_for_type : _isolation_types ) {
+      for( const auto& isolator : isolators_for_type ) {
+	isolator->getEventSetupInfo(es);
+      }
+    }
+  }
+
+  void  PFIsolationSumProducerForPUPPI::
+  produce(edm::Event& ev, const edm::EventSetup& es) {
+    typedef std::auto_ptr<edm::ValueMap<float> >  product_type;
+    typedef std::vector<float> product_values;
+    edm::Handle<CandView> to_isolate;
+    edm::Handle<CandView> isolate_with;
+    ev.getByToken(_to_isolate,to_isolate);
+    ev.getByToken(_isolate_with,isolate_with);
+    if(useValueMapForPUPPI)ev.getByToken(puppiValueMapToken_, puppiValueMap);
+
+    // the list of value vectors indexed as "to_isolate"
+    std::array<std::vector<product_values>,kNPFTypes> the_values;    
+    // get extra event info and setup value cache
+    unsigned i = 0;
+    for( const auto& isolators_for_type : _isolation_types ) {
+      the_values[i++].resize(isolators_for_type.size());
+      for( const auto& isolator : isolators_for_type ) {
+	isolator->getEventInfo(ev);
+      }
+    }
+    reco::PFCandidate helper; // to translate pdg id to type    
+    // loop over the candidates we are isolating and fill the values
+    for( size_t c = 0; c < to_isolate->size(); ++c ) {
+      auto cand_to_isolate = to_isolate->ptrAt(c);
+      std::array<std::vector<float>,kNPFTypes> cand_values;      
+      unsigned k = 0;
+      for( const auto& isolators_for_type : _isolation_types ) {
+	cand_values[k].resize(isolators_for_type.size());
+	for( auto& value : cand_values[k] ) value = 0.0;
+	++k;
+      }
+      for( size_t ic = 0; ic < isolate_with->size(); ++ic ) {
+        auto isocand = isolate_with->ptrAt(ic);
+        edm::Ptr<pat::PackedCandidate> aspackedCandidate(isocand);
+        auto isotype = helper.translatePdgIdToType(isocand->pdgId());	
+	      const auto& isolations = _isolation_types[isotype];	
+    	   for( unsigned i = 0; i < isolations.size(); ++ i  ) {
+    	  if( isolations[i]->isInIsolationCone(cand_to_isolate,isocand) ) {
+          double puppiWeight = 0.;
+    	    if (!useValueMapForPUPPI) puppiWeight = aspackedCandidate -> puppiWeight(); // if miniAOD, take puppiWeight directly from the object
+          else  puppiWeight = (*puppiValueMap)[isocand]; // if AOD, take puppiWeight from the valueMap
+          if (puppiWeight > 0.)cand_values[isotype][i] += (isocand->pt())*puppiWeight; // this is basically the main change to Lindsey's code: scale pt with puppiWeight for candidates with puppiWeight > 0.
+    	  }
+    	}
+    }
+      // add this candidate to isolation value list
+      for( unsigned i = 0; i < kNPFTypes; ++i ) {
+	for( unsigned j = 0; j < cand_values[i].size(); ++j ) {
+	  the_values[i][j].push_back(cand_values[i][j]);
+	}
+      }
+    }
+    // fill and put all products
+    for( unsigned i = 0; i < kNPFTypes; ++ i ) {
+      for( unsigned j = 0; j < the_values[i].size(); ++j ) {
+	product_type the_product( new edm::ValueMap<float> );
+	edm::ValueMap<float>::Filler fillerprod(*the_product);
+	fillerprod.insert(to_isolate, 
+			  the_values[i][j].begin(),
+			  the_values[i][j].end());
+	fillerprod.fill();
+	ev.put(the_product,_product_names[i][j]);
+      }
+    }
+  }
+}

RecoEgamma/EgammaIsolationAlgos/plugins/isolation_cones/ElectronPFIsolationWithMapBasedVeto.cc

@@ -0,0 +1,169 @@
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+
+#include "DataFormats/Candidate/interface/Candidate.h"
+#include "DataFormats/Candidate/interface/CandidateFwd.h"
+
+#include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
+#include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
+
+#include "DataFormats/PatCandidates/interface/Electron.h"
+
+#include "DataFormats/PatCandidates/interface/Photon.h"
+
+#include "DataFormats/Common/interface/ValueMap.h"
+#include "FWCore/Framework/interface/ConsumesCollector.h"
+
+#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
+#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
+
+#include "PhysicsTools/IsolationAlgos/interface/IsoDepositVetoFactory.h"
+#include "DataFormats/RecoCandidate/interface/IsoDepositVetos.h"
+
+#include "PhysicsTools/IsolationAlgos/interface/CITKIsolationConeDefinitionBase.h"
+#include "DataFormats/Math/interface/deltaR.h"
+
+#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
+
+
+#include <unordered_map>
+namespace reco {
+  typedef edm::Ptr<reco::GsfElectron> GsfElectronPtr;
+}
+
+namespace pat {
+  typedef edm::Ptr<pat::PackedCandidate> PackedCandidatePtr;
+  typedef edm::Ptr<pat::Electron> patElectronPtr;
+}
+namespace{
+	// This template function finds whether theCandidate is in thefootprint
+	// // collection. It is templated to be able to handle both reco and pat
+	// // photons (from AOD and miniAOD, respectively).
+	//
+	template <class T, class U>
+	bool isInFootprint(const T& thefootprint, const U& theCandidate) 
+	{
+	    for ( auto itr = thefootprint.begin(); itr != thefootprint.end(); ++itr ) 
+    	    { 
+      	      if( itr->key() == theCandidate.key() ) return true;
+   	    }
+ 	    return false;
+	}
+        //This function is needed because pfNoPileUpCandidates have changed keys, 
+        ////and thus the solution is to use sourceCandidatePtr(0) 
+        // // This function *shouldn't be used for packedCandidate*
+       template <class T, class U>
+        bool isInFootprintAlternative(const T& thefootprint, const U& theCandidate)
+        {
+            for ( auto itr = thefootprint.begin(); itr != thefootprint.end(); ++itr )
+            {
+              if( itr->key() == theCandidate->sourceCandidatePtr(0).key() ) return true;
+            }
+            return false;
+        }
+
+}
+class ElectronPFIsolationWithMapBasedVeto : public citk::IsolationConeDefinitionBase {
+public:
+  ElectronPFIsolationWithMapBasedVeto(const edm::ParameterSet& c) :
+    citk::IsolationConeDefinitionBase(c),
+    _isolateAgainst(c.getParameter<std::string>("isolateAgainst")), //isolate against either h+, h0 or gamma
+    _miniAODVertexCodes(c.getParameter<std::vector<unsigned> >("miniAODVertexCodes")){} //quality flags to be used for association with the vertex configurable, the vertex can be chosen
+
+  ElectronPFIsolationWithMapBasedVeto(const ElectronPFIsolationWithMapBasedVeto&) = delete;
+  ElectronPFIsolationWithMapBasedVeto& operator=(const ElectronPFIsolationWithMapBasedVeto&) =delete;
+
+  bool isInIsolationCone(const reco::CandidatePtr& photon,
+			 const reco::CandidatePtr& other) const override final;
+
+
+  // this object is needed for reco case
+ edm::Handle< edm::ValueMap<std::vector<reco::PFCandidateRef > > > particleBasedIsolationMap;			 
+ edm::EDGetTokenT<edm::ValueMap<std::vector<reco::PFCandidateRef > > > particleBasedIsolationToken_;
+
+ virtual void getEventInfo(const edm::Event& iEvent)
+  {
+      iEvent.getByToken(particleBasedIsolationToken_, particleBasedIsolationMap); 
+  };			 
+
+
+  //As far as I understand now, the object particleBasedIsolationMap should be fixed, so we don't configure the name
+  void setConsumes(edm::ConsumesCollector iC)
+  {
+      particleBasedIsolationToken_ = iC.mayConsume<edm::ValueMap<std::vector<reco::PFCandidateRef > > >(edm::InputTag("particleBasedIsolation", "gedGsfElectrons"));
+  }
+
+  //! Destructor
+  virtual ~ElectronPFIsolationWithMapBasedVeto(){};
+
+
+private:    
+  const std::string _isolateAgainst, _vertexCollection;
+  const std::vector<unsigned> _miniAODVertexCodes;
+
+
+};
+
+DEFINE_EDM_PLUGIN(CITKIsolationConeDefinitionFactory,
+		  ElectronPFIsolationWithMapBasedVeto,
+		  "ElectronPFIsolationWithMapBasedVeto");
+
+
+//This function defines whether particular PFCandidate is inside of isolation cone of photon or not by checking deltaR and whether footprint removal for this candidate should be done. Additionally, for miniAOD charged hadrons from the PV are considered. *** For AOD this should be done by the corresponding sequence beforehand!!! ***
+
+bool ElectronPFIsolationWithMapBasedVeto::isInIsolationCone(const reco::CandidatePtr& electron,  const reco::CandidatePtr& pfCandidate  ) const {
+
+  //convert the electron and candidate objects to the corresponding pat or reco objects. What is used depends on what is user running on: miniAOD or AOD
+  pat::patElectronPtr aspat_electronptr(electron);
+
+  pat::PackedCandidatePtr aspacked(pfCandidate);
+  reco::PFCandidatePtr aspf(pfCandidate);
+
+
+  bool inFootprint = false;
+  bool result = true;
+  const float deltar2 = reco::deltaR2(*electron,*pfCandidate); //calculate deltaR2 distance between PFCandidate and photon
+
+  // dealing here with patObjects: miniAOD case
+  if ( aspacked.get() )    
+  {
+    inFootprint = isInFootprint(aspat_electronptr ->associatedPackedPFCandidates(),aspacked);
+
+    //checking if the charged candidates come from the appropriate vertex
+    if( aspacked->charge() != 0 ) 
+    {
+      bool is_vertex_allowed = false;
+      for( const unsigned vtxtype : _miniAODVertexCodes ) 
+      {
+      	if( vtxtype == aspacked->fromPV() ) {
+      	  is_vertex_allowed = true;
+      	  break;
+      	}
+      }      
+
+     result &= (is_vertex_allowed);
+    }
+     //return true if the candidate is inside the cone and not in the footprint
+    result &= deltar2 < _coneSize2 && (!inFootprint);
+   }
+
+  // dealing here with recoObjects: AOD case
+  else if ( aspf.get())
+  {
+      inFootprint = isInFootprintAlternative((*particleBasedIsolationMap)[electron], pfCandidate); 
+      //return true if the candidate is inside the cone and not in the footprint
+      result &= deltar2 < _coneSize2 && (!inFootprint);
+  }
+
+  // throw exception if it is not a patObject or recoObject
+  else {
+    throw cms::Exception("InvalidIsolationInput")
+      << "The supplied candidate to be used as isolation "
+      << "was neither a reco::Photon nor a pat::Photon!";
+  }
+
+  return result;
+}