MuIsolatorResultProducer.h

Go to the documentation of this file.

#ifndef MuonIsolationProducers_MuIsolatorResultProducer_H
#define MuonIsolationProducers_MuIsolatorResultProducer_H

#include "FWCore/Framework/interface/EDProducer.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "RecoMuon/MuonIsolation/interface/MuIsoBaseIsolator.h"

#include "FWCore/Framework/interface/Event.h"

#include "DataFormats/Common/interface/RefToBase.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "DataFormats/BeamSpot/interface/BeamSpot.h"
#include "DataFormats/Common/interface/ValueMap.h"


#include <string>

namespace edm { class EventSetup; }

struct muisorhelper {
  typedef muonisolation::MuIsoBaseIsolator Isolator;
  typedef Isolator::Result Result;
  typedef Isolator::ResultType ResultType;
  typedef std::vector<Result> Results;
  typedef Isolator::DepositContainer DepositContainer;
  
  template<typename BT>
  class CandMap {
  public:
    typedef typename edm::RefToBase<BT> key_type;
    typedef typename edm::Handle<edm::View<BT> > handle_type;
    typedef DepositContainer value_type;
    typedef std::pair<key_type, value_type> pair_type;
    typedef typename std::vector<pair_type > map_type;
    typedef typename map_type::iterator iterator;
    
    map_type& get() { return cMap_;}
    const map_type& get() const {return cMap_;}
    const handle_type& handle() const { return handle_;}
    void setHandle(const handle_type& rhs) { handle_ = rhs;}
  private:
    map_type cMap_;
    handle_type handle_;
  };
  
};

template <typename BT= reco::Candidate>
  class MuIsolatorResultProducer : public edm::EDProducer {
  
  public:
  MuIsolatorResultProducer(const edm::ParameterSet&);
  
  ~MuIsolatorResultProducer() override;
  
  void produce(edm::Event&, const edm::EventSetup&) override;
  
  private:
  typedef muisorhelper::Isolator Isolator;
  typedef muisorhelper::Result Result;
  typedef muisorhelper::ResultType ResultType;
  typedef muisorhelper::Results Results;
  typedef muisorhelper::DepositContainer DepositContainer;
  
  typedef muisorhelper::CandMap<BT> CandMap;
  
  struct DepositConf { 
    edm::InputTag tag; 
    double weight; 
    double threshold;
  };
  
  struct VetoCuts { 
    bool selectAll; 
    double muAbsEtaMax; 
    double muPtMin;
    double muAbsZMax;
    double muD0Max;
  };
  
  
  
  void callWhatProduces();
  
  unsigned int initAssociation(edm::Event& event, CandMap& candMapT) const;
  
  void initVetos(std::vector<reco::IsoDeposit::Vetos*>& vetos, CandMap& candMap) const;
  
  template <typename RT>
  void writeOutImpl(edm::Event& event, const CandMap& candMapT, const Results& results) const;
  
  void writeOut(edm::Event& event, const CandMap& candMap, const Results& results) const;
  
  edm::ParameterSet theConfig;
  std::vector<DepositConf> theDepositConfs;
  
  bool theRemoveOtherVetos;
  VetoCuts theVetoCuts;
  
  Isolator* theIsolator;
  ResultType theResultType;
  
  
  std::string theBeamlineOption;
  edm::InputTag theBeamSpotLabel;
  reco::TrackBase::Point theBeam;
};


template<typename BT>
template<typename RT> inline
void MuIsolatorResultProducer<BT>::writeOutImpl(edm::Event& event, const CandMap& candMapT, 
                        const Results& results) const {
  std::vector<RT> resV(results.size());   
  for (unsigned int i = 0; i< resV.size(); ++i) resV[i] = results[i].val<RT>(); 
  auto outMap = std::make_unique<edm::ValueMap<RT>>(); 
  typename edm::ValueMap<RT>::Filler filler(*outMap); 

  if (!candMapT.get().empty()){
    filler.insert(candMapT.handle(), resV.begin(), resV.end()); 
    filler.fill(); 
  }
  
  event.put(std::move(outMap)); 
}


template<typename BT>
inline void MuIsolatorResultProducer<BT>::writeOut(edm::Event& event, 
                           const CandMap& candMapT, 
                           const Results& results) const {
  
  std::string metname = "RecoMuon|MuonIsolationProducers";
  LogDebug(metname)<<"Before calling writeOutMap  with result type "<<theIsolator->resultType();
  
  if (theResultType == Isolator::ISOL_INT_TYPE) writeOutImpl<int>(event, candMapT, results);
  if (theResultType == Isolator::ISOL_FLOAT_TYPE) writeOutImpl<float>(event, candMapT, results);
  if (theResultType == Isolator::ISOL_BOOL_TYPE) writeOutImpl<bool>(event, candMapT, results);
}


template<typename BT>
inline void MuIsolatorResultProducer<BT>::callWhatProduces() {
  if (theResultType == Isolator::ISOL_FLOAT_TYPE) produces<edm::ValueMap<float> >();
  if (theResultType == Isolator::ISOL_INT_TYPE  ) produces<edm::ValueMap<int> >();
  if (theResultType == Isolator::ISOL_BOOL_TYPE ) produces<edm::ValueMap<bool> >();      
}

// Framework
#include "FWCore/Framework/interface/EDProducer.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "DataFormats/Common/interface/Handle.h"

#include "FWCore/Framework/interface/ESHandle.h"

#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
#include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
#include "DataFormats/Common/interface/RefToBase.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"

#include "RecoMuon/MuonIsolation/interface/Range.h"
#include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h"

#include "RecoMuon/MuonIsolation/interface/IsolatorByDeposit.h"
#include "RecoMuon/MuonIsolation/interface/IsolatorByDepositCount.h"
#include "RecoMuon/MuonIsolation/interface/IsolatorByNominalEfficiency.h"

#include <string>


template<typename BT>
MuIsolatorResultProducer<BT>::MuIsolatorResultProducer(const edm::ParameterSet& par) :
  theConfig(par),
  theRemoveOtherVetos(par.getParameter<bool>("RemoveOtherVetos")),
  theIsolator(nullptr),
  theBeam(0,0,0)
{
  LogDebug("RecoMuon|MuonIsolation")<<" MuIsolatorResultProducer CTOR";
  
  std::vector<edm::ParameterSet> depositInputs = 
    par.getParameter<std::vector<edm::ParameterSet> >("InputMuIsoDeposits");    

  std::vector<double> dWeights( depositInputs.size());
  std::vector<double> dThresholds( depositInputs.size());

  for (unsigned int iDep = 0; iDep < depositInputs.size(); ++iDep){
    DepositConf dConf;
    dConf.tag = depositInputs[iDep].getParameter<edm::InputTag>("DepositTag");
    dConf.weight = depositInputs[iDep].getParameter<double>("DepositWeight");
    dConf.threshold = depositInputs[iDep].getParameter<double>("DepositThreshold");
    
    dWeights[iDep] = dConf.weight;
    dThresholds[iDep] = dConf.threshold;

    theDepositConfs.push_back(dConf);
  }

  edm::ParameterSet isoPset = par.getParameter<edm::ParameterSet>("IsolatorPSet");
  std::string isolatorType = isoPset.getParameter<std::string>("ComponentName");
  if ( isolatorType == "IsolatorByDeposit"){    
    std::string coneSizeType = isoPset.getParameter<std::string>("ConeSizeType");
    if (coneSizeType == "FixedConeSize"){
      float coneSize(isoPset.getParameter<double>("coneSize"));

      theIsolator = new muonisolation::IsolatorByDeposit(coneSize, dWeights, dThresholds);

      //      theIsolator = new IsolatorByDeposit(isoPset);
    } else if (coneSizeType == "CutsConeSize"){
//       Cuts cuts(isoPset.getParameter<edm::ParameterSet>("CutsPSet"));
      
//       theIsolator = new IsolatorByDeposit(coneSize, dWeights, dThresholds);
    }
  } else if ( isolatorType == "IsolatorByNominalEfficiency"){
    theIsolator = new muonisolation::IsolatorByNominalEfficiency("noname", std::vector<std::string>(1,"8:0.97"), dWeights);
  } else if ( isolatorType == "IsolatorByDepositCount"){    
    std::string coneSizeType = isoPset.getParameter<std::string>("ConeSizeType");
    if (coneSizeType == "FixedConeSize"){
      float coneSize(isoPset.getParameter<double>("coneSize"));
      
      theIsolator = new muonisolation::IsolatorByDepositCount(coneSize, dThresholds);
      
      //      theIsolator = new IsolatorByDeposit(isoPset);
    } else if (coneSizeType == "CutsConeSize"){
      //       Cuts cuts(isoPset.getParameter<edm::ParameterSet>("CutsPSet"));
      
      //       theIsolator = new IsolatorByDeposit(coneSize, dWeights, dThresholds);
    }
  }
  
  if (theIsolator == nullptr ){
    edm::LogError("MuonIsolationProducers")<<"Failed to initialize an isolator";
  }
  theResultType = theIsolator->resultType();

  callWhatProduces();

  if (theRemoveOtherVetos){
    edm::ParameterSet vetoPSet = par.getParameter<edm::ParameterSet>("VetoPSet");
    theVetoCuts.selectAll = vetoPSet.getParameter<bool>("SelectAll");

    if (! theVetoCuts.selectAll){
      theVetoCuts.muAbsEtaMax = vetoPSet.getParameter<double>("MuAbsEtaMax");
      theVetoCuts.muPtMin     = vetoPSet.getParameter<double>("MuPtMin");
      theVetoCuts.muAbsZMax   = vetoPSet.getParameter<double>("MuAbsZMax");
      theVetoCuts.muD0Max      = vetoPSet.getParameter<double>("MuD0Max");    
      theBeamlineOption = par.getParameter<std::string>("BeamlineOption");
      theBeamSpotLabel = par.getParameter<edm::InputTag>("BeamSpotLabel");
    }
  }
}

template<typename BT>
MuIsolatorResultProducer<BT>::~MuIsolatorResultProducer(){
  if (theIsolator) delete theIsolator;
  LogDebug("RecoMuon|MuIsolatorResultProducer")<<" MuIsolatorResultProducer DTOR";
}


template<typename BT>
void MuIsolatorResultProducer<BT>::produce(edm::Event& event, const edm::EventSetup& eventSetup){
  
  std::string metname = "RecoMuon|MuonIsolationProducers";
  LogDebug(metname)<<" Muon Deposit producing..."
           <<" BEGINING OF EVENT " <<"================================";

  theBeam = reco::TrackBase::Point(0,0, 0); 
 
  if(theRemoveOtherVetos && ! theVetoCuts.selectAll){
    if (theBeamlineOption == "BeamSpotFromEvent"){ 
      //pick beamSpot 
      reco::BeamSpot beamSpot; 
      edm::Handle<reco::BeamSpot> beamSpotH; 
      
      event.getByLabel(theBeamSpotLabel,beamSpotH); 
      
      if (beamSpotH.isValid()){ 
    theBeam = beamSpotH->position(); 
    LogTrace(metname)<<"Extracted beam point at "<<theBeam<<std::endl; 
      } 
    }
  } 

  CandMap candMapT;
  
  unsigned int colSize = initAssociation(event, candMapT);

  Results results(colSize);

  std::vector<reco::IsoDeposit::Vetos*> vetoDeps(theDepositConfs.size(), nullptr);

  if (colSize != 0){
    if (theRemoveOtherVetos){

      initVetos(vetoDeps, candMapT);
    }

    for (unsigned int muI=0; muI < colSize; ++muI){
      results[muI] = theIsolator->result(candMapT.get()[muI].second, *(candMapT.get()[muI].first));
      
      if (results[muI].typeF()!= theIsolator->resultType()){
    edm::LogError("MuonIsolationProducers")<<"Failed to get result from the isolator";
      }
    }
    
  }

  LogDebug(metname)<<"Ready to write out results of size "<<results.size();
  writeOut(event, candMapT, results);

  for(unsigned int iDep = 0; iDep< vetoDeps.size(); ++iDep){
    if (vetoDeps[iDep]){
      delete vetoDeps[iDep];
      vetoDeps[iDep] = nullptr;
    }
  }
}

template<typename BT>
unsigned int
MuIsolatorResultProducer<BT>::initAssociation(edm::Event& event, CandMap& candMapT) const {
  std::string metname = "RecoMuon|MuonIsolationProducers";
  
  typedef reco::IsoDepositMap::container CT;

  for (unsigned int iMap = 0; iMap < theDepositConfs.size(); ++iMap){
    edm::Handle<reco::IsoDepositMap> depH;
    event.getByLabel(theDepositConfs[iMap].tag, depH);
    LogDebug(metname) <<"Got Deposits of size "<<depH->size();
    if (depH->empty()) continue;

    typename CandMap::handle_type keyH;
    event.get(depH->begin().id(), keyH);
    candMapT.setHandle(keyH);
    typename CT::const_iterator depHCI = depH->begin().begin();
    typename CT::const_iterator depEnd = depH->begin().end();
    unsigned int keyI=0;
    for (; depHCI != depEnd; ++depHCI, ++keyI){

      typename CandMap::key_type muPtr(keyH->refAt(keyI));
      if (iMap == 0) candMapT.get().push_back(typename CandMap::pair_type(muPtr, DepositContainer(theDepositConfs.size())));
      typename CandMap::iterator muI = candMapT.get().begin();
      for (; muI != candMapT.get().end(); ++muI){
    if (muI->first == muPtr) break;
      }
      if (muI->first != muPtr){
    edm::LogError("MuonIsolationProducers")<<"Failed to align muon map";
      }
      muI->second[iMap].dep = &*depHCI; 
    }
  }

  LogDebug(metname)<<"Picked and aligned nDeps = "<<candMapT.get().size();
  return candMapT.get().size();
}

template <typename BT >
void MuIsolatorResultProducer<BT>::initVetos(std::vector<reco::IsoDeposit::Vetos*>& vetos, CandMap& candMapT) const {
  std::string metname = "RecoMuon|MuonIsolationProducers";
  

  if (theRemoveOtherVetos){
    LogDebug(metname)<<"Start checking for vetos based on input/expected vetos.size of "<<vetos.size()
             <<" passed at "<<&vetos
             <<" and an input map.size of "<<candMapT.get().size();

    unsigned int muI = 0;
    for (; muI < candMapT.get().size(); ++muI) {
      typename CandMap::key_type mu = candMapT.get()[muI].first;
      double d0 = ( (mu->vx() - theBeam.x() )* mu->py() - (mu->vy() - theBeam.y())* mu->px() ) / mu->pt();
      LogDebug(metname)<<"Muon at index "<<muI;
      if (theVetoCuts.selectAll 
      || (fabs(mu->eta()) < theVetoCuts.muAbsEtaMax
          && mu->pt() > theVetoCuts.muPtMin
          && fabs(mu->vz()) < theVetoCuts.muAbsZMax
          && fabs(d0) < theVetoCuts.muD0Max
          )
      ){
    LogDebug(metname)<<"muon passes the cuts";
    for (unsigned int iDep =0; iDep < candMapT.get()[muI].second.size(); ++iDep){
      if (vetos[iDep] == nullptr) vetos[iDep] = new reco::IsoDeposit::Vetos();

      vetos[iDep]->push_back(candMapT.get()[muI].second[iDep].dep->veto());
    }
      }
    }

    LogDebug(metname)<<"Assigning vetos";
    muI = 0;
    for (; muI < candMapT.get().size(); ++muI) {
      for(unsigned int iDep =0; iDep < candMapT.get()[muI].second.size(); ++iDep){
    candMapT.get()[muI].second[iDep].vetos = vetos[iDep];
      }
    }
    LogDebug(metname)<<"Done with vetos";
  }
}


  
#endif