/ElectroWeakAnalysis/ZMuMu/plugins/ZGlobalVsSAIsolationAnalyzer.cc

  1#include "DataFormats/Common/interface/AssociationVector.h"
  2#include "FWCore/Framework/interface/EDAnalyzer.h"
  3#include "FWCore/Utilities/interface/EDMException.h"
  4#include "DataFormats/Candidate/interface/Candidate.h"
  5#include "DataFormats/Candidate/interface/Particle.h"
  6#include "DataFormats/Candidate/interface/CandidateFwd.h"
  7#include "FWCore/ParameterSet/interface/ParameterSet.h"
  8#include "FWCore/Framework/interface/Event.h"
  9#include "FWCore/Framework/interface/EventSetup.h"
 10#include "FWCore/Utilities/interface/InputTag.h"
 11#include "FWCore/ServiceRegistry/interface/Service.h"
 12#include "CommonTools/UtilAlgos/interface/TFileService.h"
 13#include "DataFormats/Math/interface/deltaR.h"
 14#include "DataFormats/PatCandidates/interface/Muon.h"
 15#include "DataFormats/TrackReco/interface/Track.h"
 16#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
 17#include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
 18#include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
 19#include "DataFormats/PatCandidates/interface/Isolation.h"
 20#include "DataFormats/Common/interface/ValueMap.h"
 21#include "DataFormats/PatCandidates/interface/GenericParticle.h"
 22#include "DataFormats/RecoCandidate/interface/IsoDepositVetos.h"
 23#include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h"
 24#include <vector>
 25#include <string>
 26#include <iostream>
 27
 28using namespace edm;
 29using namespace std;
 30using namespace reco;
 31using namespace isodeposit;
 32
 33class ZGlobalVsSAIsolationAnalyzer : public edm::EDAnalyzer {
 34public:
 35  typedef math::XYZVector Vector;
 36  ZGlobalVsSAIsolationAnalyzer(const edm::ParameterSet& cfg);
 37private:
 38  virtual void analyze(const edm::Event& event, const edm::EventSetup& setup) override;
 39  virtual void endJob() override;
 40  InputTag src_;
 41  double dRVeto;
 42  double dRTrk, dREcal, dRHcal;
 43  double ptThreshold, etEcalThreshold, etHcalThreshold;
 44  double alpha, beta;
 45  double isoCut_;
 46  unsigned long selGlobal_, selSA_, totGlobal_, totSA_;
 47  bool isolated(const Direction & dir, const pat::IsoDeposit * trkIsoDep, 
 48		const pat::IsoDeposit * ecalIsoDep, const pat::IsoDeposit * hcalIsoDep);
 49  void evaluate(const reco::Candidate* dau);
 50};
 51
 52ZGlobalVsSAIsolationAnalyzer::ZGlobalVsSAIsolationAnalyzer(const ParameterSet& cfg):
 53  src_(cfg.getParameter<InputTag>("src")),
 54  dRVeto(cfg.getParameter<double>("veto")),
 55  dRTrk(cfg.getParameter<double>("deltaRTrk")),
 56  dREcal(cfg.getParameter<double>("deltaREcal")),
 57  dRHcal(cfg.getParameter<double>("deltaRHcal")),
 58  ptThreshold(cfg.getParameter<double>("ptThreshold")),
 59  etEcalThreshold(cfg.getParameter<double>("etEcalThreshold")),
 60  etHcalThreshold(cfg.getParameter<double>("etHcalThreshold")),
 61  alpha(cfg.getParameter<double>("alpha")),
 62  beta(cfg.getParameter<double>("beta")),
 63  isoCut_(cfg.getParameter<double>("isoCut")),
 64  selGlobal_(0), selSA_(0), totGlobal_(0), totSA_(0) {
 65}
 66
 67bool ZGlobalVsSAIsolationAnalyzer::isolated(const Direction & dir, const pat::IsoDeposit * trkIsoDep, 
 68					    const pat::IsoDeposit * ecalIsoDep, const pat::IsoDeposit * hcalIsoDep) {
 69  IsoDeposit::AbsVetos vetoTrk, vetoEcal, vetoHcal;
 70  vetoTrk.push_back(new ConeVeto(dir, dRVeto));
 71  vetoTrk.push_back(new ThresholdVeto(ptThreshold));
 72  vetoEcal.push_back(new ConeVeto(dir, 0.));
 73  vetoEcal.push_back(new ThresholdVeto(etEcalThreshold));
 74  vetoHcal.push_back(new ConeVeto(dir, 0.));
 75  vetoHcal.push_back(new ThresholdVeto(etHcalThreshold));
 76  
 77  double trkIso = trkIsoDep->sumWithin(dir, dRTrk, vetoTrk);
 78  double ecalIso = ecalIsoDep->sumWithin(dir, dREcal, vetoEcal);
 79  double hcalIso = hcalIsoDep->sumWithin(dir, dRHcal, vetoHcal);
 80  double iso = alpha*((0.5*(1+beta)*ecalIso) + (0.5*(1-beta)*hcalIso)) + (1-alpha)*trkIso;
 81  return iso < isoCut_;      
 82}
 83
 84void ZGlobalVsSAIsolationAnalyzer::evaluate(const reco::Candidate* dau) {
 85  const pat::Muon * mu = dynamic_cast<const pat::Muon *>(&*dau->masterClone());
 86  if(mu == 0) throw Exception(errors::InvalidReference) << "Daughter is not a muon!\n";
 87  const pat::IsoDeposit * trkIsoDep = mu->isoDeposit(pat::TrackIso);
 88  const pat::IsoDeposit * ecalIsoDep = mu->isoDeposit(pat::EcalIso);
 89  const pat::IsoDeposit * hcalIsoDep = mu->isoDeposit(pat::HcalIso);
 90  // global muon
 91  {
 92    Direction dir = Direction(mu->eta(), mu->phi());
 93    if(isolated(dir, trkIsoDep, ecalIsoDep, hcalIsoDep)) selGlobal_++;
 94    totGlobal_++;
 95  }
 96  // stand-alone
 97  {
 98    TrackRef sa = dau->get<TrackRef,reco::StandAloneMuonTag>();
 99    Direction dir = Direction(sa->eta(), sa->phi());
100    if(isolated(dir, trkIsoDep, ecalIsoDep, hcalIsoDep)) selSA_++;
101    totSA_++;   
102  }
103}
104
105void ZGlobalVsSAIsolationAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup) {
106  Handle<CandidateView> dimuons;
107  event.getByLabel(src_, dimuons);
108  for(unsigned int i=0; i< dimuons->size(); ++ i) {
109    const Candidate & zmm = (* dimuons)[i];
110    evaluate(zmm.daughter(0));
111    evaluate(zmm.daughter(1));
112  }    
113}
114  
115void ZGlobalVsSAIsolationAnalyzer::endJob() {
116  cout << "Isolation efficiency report:" << endl;
117  double eff, err;
118  eff =  double(selGlobal_)/double(totGlobal_);
119  err = sqrt(eff*(1.-eff)/double(totGlobal_));
120  cout <<"Global: " << selGlobal_ << "/" << totGlobal_ << " = " <<eff <<"+/-" << err<< endl;
121  eff = double(selSA_)/double(totSA_);
122  err = sqrt(eff*(1.-eff)/double(totSA_));  
123  cout <<"St.Al.: " << selSA_ << "/" << totSA_ << " = " << eff <<"+/-" << err << endl;
124}
125
126#include "FWCore/Framework/interface/MakerMacros.h"
127
128DEFINE_FWK_MODULE(ZGlobalVsSAIsolationAnalyzer);