WMuNuSelector.cc

Go to the documentation of this file.

//                                                                                                                      //
//                                    WMuNuSelector based on WMuNuCandidates                                            //
//                                                                                                                      //
//                                                                                                                      //
//    Filter of WMuNuCandidates for Analysis                                                                            //
//    --> From a WMuNuCandidate collection                                                                              //
//    --> Pre-Selection of events based in event cuts (trigger, Z rejection, ttbar rejection)                           //
//    --> The Ws are selected from the highest pt muon in the event (applying the standard WMuNu Selection cuts)        //
//                                                                                                                      //
//    --> Be careful: if this Selector is used as a filter for further analysis you still have to make sure that        //
//                    the W Candidate you use for your modules is the first one in the collection!!                     //
//                                                                                                                      //
//    Optionally, plots selection variables sequentially after cuts,                                                    //
//    and 2D histograms for background determination.                                                                   //
//                                                                                                                      //
//    For basic plots before & after cuts (without Candidate formalism), use WMuNuValidator.cc                          //
//                                                                                                                      //

#include "FWCore/Framework/interface/EDFilter.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/Event.h"
#include "DataFormats/Common/interface/TriggerResults.h"
#include "DataFormats/Common/interface/View.h"
#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/JetReco/interface/Jet.h"
#include "DataFormats/BeamSpot/interface/BeamSpot.h"
#include "AnalysisDataFormats/EWK/interface/WMuNuCandidate.h"
#include "TH1D.h"
#include "TH2D.h"

class WMuNuSelector : public edm::EDFilter {
public:
  WMuNuSelector (const edm::ParameterSet &);
  bool filter(edm::Event&, const edm::EventSetup&) override;
  void beginJob() override;
  void endJob() override;
  void init_histograms();
private:
  bool plotHistograms_;
  edm::EDGetTokenT<edm::TriggerResults> trigToken_;
  edm::EDGetTokenT<edm::View<reco::Muon> > muonToken_;
  edm::EDGetTokenT<edm::View<reco::Jet> > jetToken_;
  edm::EDGetTokenT<reco::BeamSpot> beamSpotToken_;
  edm::EDGetTokenT<reco::WMuNuCandidateCollection> WMuNuCollectionToken_;
  const std::string muonTrig_;
  double ptThrForZ1_;
  double ptThrForZ2_;
  double eJetMin_;
  int nJetMax_;
  double ptCut_;
  double etaCut_;
  bool isRelativeIso_;
  bool isCombinedIso_;
  double isoCut03_;
  double mtMin_;
  double mtMax_;
  double metMin_;
  double metMax_;
  double acopCut_;

  double dxyCut_;
  double normalizedChi2Cut_;
  int trackerHitsCut_;
  bool isAlsoTrackerMuon_;

  int selectByCharge_;

  double nall;
  double ntrig, npresel;
  double nsel;
  double ncharge;
  double nkin, nid,nacop,niso,nmass;



  std::map<std::string,TH1D*> h1_;
  std::map<std::string,TH2D*> h2_;

};
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/Common/interface/Handle.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"

#include "DataFormats/TrackReco/interface/Track.h"

#include "DataFormats/MuonReco/interface/MuonSelectors.h"

#include "DataFormats/METReco/interface/MET.h"

#include "DataFormats/GeometryVector/interface/Phi.h"

#include "FWCore/Common/interface/TriggerNames.h"


using namespace edm;
using namespace std;
using namespace reco;

WMuNuSelector::WMuNuSelector( const ParameterSet & cfg ) :
      // Fast selection (no histograms)
      plotHistograms_(cfg.getUntrackedParameter<bool> ("plotHistograms", true)),

      // Input collections
      trigToken_(consumes<TriggerResults>(cfg.getUntrackedParameter<edm::InputTag> ("TrigTag", edm::InputTag("TriggerResults::HLT")))),
      muonToken_(consumes<View<Muon> >(cfg.getUntrackedParameter<edm::InputTag> ("MuonTag", edm::InputTag("muons")))),
      jetToken_(consumes<View<Jet> >(cfg.getUntrackedParameter<edm::InputTag> ("JetTag", edm::InputTag("sisCone5CaloJets")))),
      beamSpotToken_(consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"))),
      WMuNuCollectionToken_(consumes<reco::WMuNuCandidateCollection>(cfg.getUntrackedParameter<edm::InputTag> ("WMuNuCollectionTag", edm::InputTag("WMuNus")))),


      // Preselection cuts
      muonTrig_(cfg.getUntrackedParameter<std::string>("MuonTrig", "HLT_Mu9")),
      ptThrForZ1_(cfg.getUntrackedParameter<double>("PtThrForZ1", 20.)),
      ptThrForZ2_(cfg.getUntrackedParameter<double>("PtThrForZ2", 10.)),
      eJetMin_(cfg.getUntrackedParameter<double>("EJetMin", 999999.)),
      nJetMax_(cfg.getUntrackedParameter<int>("NJetMax", 999999)),


      // Main cuts
      ptCut_(cfg.getUntrackedParameter<double>("PtCut", 25.)),
      etaCut_(cfg.getUntrackedParameter<double>("EtaCut", 2.1)),
      isRelativeIso_(cfg.getUntrackedParameter<bool>("IsRelativeIso", true)),
      isCombinedIso_(cfg.getUntrackedParameter<bool>("IsCombinedIso", false)),
      isoCut03_(cfg.getUntrackedParameter<double>("IsoCut03", 0.1)),
      mtMin_(cfg.getUntrackedParameter<double>("MtMin", 50.)),
      mtMax_(cfg.getUntrackedParameter<double>("MtMax", 200.)),
      metMin_(cfg.getUntrackedParameter<double>("MetMin", -999999.)),
      metMax_(cfg.getUntrackedParameter<double>("MetMax", 999999.)),
      acopCut_(cfg.getUntrackedParameter<double>("AcopCut", 2.)),

      // Muon quality cuts
      dxyCut_(cfg.getUntrackedParameter<double>("DxyCut", 0.2)),
      normalizedChi2Cut_(cfg.getUntrackedParameter<double>("NormalizedChi2Cut", 10.)),
      trackerHitsCut_(cfg.getUntrackedParameter<int>("TrackerHitsCut", 11)),
      isAlsoTrackerMuon_(cfg.getUntrackedParameter<bool>("IsAlsoTrackerMuon", true)),

      // W+/W- Selection
      selectByCharge_(cfg.getUntrackedParameter<int>("SelectByCharge", 0))
{
}

void WMuNuSelector::beginJob() {
      nall = 0;
      ntrig=0;
      npresel=0;
      ncharge = 0;
      nkin=0;
      nid=0;
      nacop=0;
      niso=0;
      nsel = 0;

   if(plotHistograms_){
     edm::Service<TFileService> fs;
     h1_["hNWCand_sel"]                  =fs->make<TH1D>("NWCand_sel","Nb. of WCandidates passing pre-selection (ordered by pt)",10,0.,10.);
     h1_["hPtMu_sel"]                    =fs->make<TH1D>("ptMu_sel","Pt mu",100,0.,100.);
     h1_["hEtaMu_sel"]                   =fs->make<TH1D>("etaMu_sel","Eta mu",50,-2.5,2.5);
     h1_["hd0_sel"]                      =fs->make<TH1D>("d0_sel","Impact parameter",1000,-1.,1.);
     h1_["hNHits_sel"]                   =fs->make<TH1D>("NumberOfValidHits_sel","Number of Hits in Silicon",100,0.,100.);
     h1_["hNormChi2_sel"]                =fs->make<TH1D>("NormChi2_sel","Chi2/ndof of global track",1000,0.,50.);
     h1_["hTracker_sel"]                 =fs->make<TH1D>("isTrackerMuon_sel","is Tracker Muon?",2,0.,2.);
     h1_["hMET_sel"]                     =fs->make<TH1D>("MET_sel","Missing Transverse Energy (GeV)", 300,0.,300.);
     h1_["hTMass_sel"]                   =fs->make<TH1D>("TMass_sel","Rec. Transverse Mass (GeV)",300,0.,300.);
     h1_["hAcop_sel"]                    =fs->make<TH1D>("Acop_sel","Mu-MET acoplanarity",50,0.,M_PI);
     h1_["hPtSum_sel"]                   =fs->make<TH1D>("ptSum_sel","Track Isolation, Sum pT (GeV)",200,0.,100.);
     h1_["hPtSumN_sel"]                  =fs->make<TH1D>("ptSumN_sel","Track Isolation, Sum pT/pT",1000,0.,10.);
     h1_["hCal_sel"]                     =fs->make<TH1D>("Cal_sel","Calorimetric isolation, HCAL+ECAL (GeV)",200,0.,100.);
     h1_["hIsoTotN_sel"]                 =fs->make<TH1D>("isoTotN_sel","(Sum pT + Cal)/pT",1000,0.,10.);
     h1_["hIsoTot_sel"]                  =fs->make<TH1D>("isoTot_sel","(Sum pT + Cal)",200,0.,100.);
     h2_["hTMass_PtSum_inclusive"]       =fs->make<TH2D>("TMass_PtSum_inclusive","Rec. Transverse Mass (GeV) vs Sum pT (GeV)",200,0.,100.,300,0.,300.);
     h2_["hTMass_PtSumNorm_inclusive"]   =fs->make<TH2D>("TMass_PtSumNorm_inclusive","Rec. Transverse Mass (GeV) vs Sum Pt / Pt", 1000,0,10,300,0,300);
     h2_["hTMass_TotIsoNorm_inclusive"]=fs->make<TH2D>("TMass_TotIsoNorm_inclusive","Rec. Transverse Mass (GeV) vs (Sum Pt + Cal)/Pt", 10000,0,10,200,0,200);
     h2_["hMET_PtSum_inclusive"]         =fs->make<TH2D>("MET_PtSum_inclusive","Missing Transverse Energy (GeV) vs Sum Pt (GeV)",200,0.,100.,300,0.,300.);
     h2_["hMET_PtSumNorm_inclusive"]     =fs->make<TH2D>("MET_PtSumNorm_inclusive","Missing Transverse Energy (GeV) vs Sum Pt/Pt",1000,0,10,300,0,300);
     h2_["hMET_TotIsoNorm_inclusive"]=fs->make<TH2D>("MET_TotIsoNorm_inclusive","Missing Transverse Energy (GeV) vs (SumPt + Cal)/Pt",10000,0,10,200,0,200);
   }
}


void WMuNuSelector::endJob() {
      double all = nall;
      double epresel = npresel/all;
      double etrig = ntrig/all;
      double ekin = nkin/all;
      double eid = nid/all;
      double eacop = nacop/all;
      double eiso = niso/all;
      double esel = nsel/all;

      LogVerbatim("") << "\n>>>>>> W SELECTION SUMMARY BEGIN >>>>>>>>>>>>>>>";
      LogVerbatim("") << "Total number of events analyzed: " << nall << " [events]";
      LogVerbatim("") << "Total number of events triggered: " << ntrig << " [events]";
      LogVerbatim("") << "Total number of events pre-selected: " << npresel << " [events]";
      LogVerbatim("") << "Total number of events after kinematic cuts: " << nkin << " [events]";
      LogVerbatim("") << "Total number of events after Muon ID cuts: " << nid << " [events]";
      LogVerbatim("") << "Total number of events after Acop cut: " << nacop << " [events]";
      LogVerbatim("") << "Total number of events after iso cut: " << niso << " [events]";
      LogVerbatim("") << "Total number of events selected: " << nsel << " [events]";
      LogVerbatim("") << "Efficiencies:";
      LogVerbatim("") << "Trigger Efficiency:                   " << "(" << setprecision(4) << etrig*100. <<" +/- "<< setprecision(2) << sqrt(etrig*(1-etrig)/all)*100. << ")%";
      LogVerbatim("") << "Pre-Selection Efficiency:             " << "(" << setprecision(4) << epresel*100. <<" +/- "<< setprecision(2) << sqrt(epresel*(1-epresel)/all)*100. << ")%";
      LogVerbatim("") << "Pt, Eta Selection Efficiency:         " << "(" << setprecision(4) << ekin*100. <<" +/- "<< setprecision(2) << sqrt(ekin*(1-ekin)/all)*100. << ")%";
      LogVerbatim("") << "MuonID Efficiency:                    " << "(" << setprecision(4) << eid*100. <<" +/- "<< setprecision(2) << sqrt(eid*(1-eid)/all)*100. << ")%";
      LogVerbatim("") << "Acop Efficiency:                      " << "(" << setprecision(4) << eacop*100. <<" +/- "<< setprecision(2) << sqrt(eacop*(1-eacop)/all)*100. << ")%";
      LogVerbatim("") << "Iso Efficiency:                       " << "(" << setprecision(4) << eiso*100. <<" +/- "<< setprecision(2) << sqrt(eiso*(1-eiso)/all)*100. << ")%";
      LogVerbatim("") << "Selection Efficiency:             " << "(" << setprecision(4) << esel*100. <<" +/- "<< setprecision(2) << sqrt(esel*(1-esel)/nall)*100. << ")%";

     if ( fabs(selectByCharge_)==1 ){
      esel = nsel/ncharge;
      LogVerbatim("") << "\n>>>>>> W+(-) SELECTION >>>>>>>>>>>>>>>";
      LogVerbatim("") << "Total number of W+(-) events pre-selected: " << ncharge << " [events]";
      LogVerbatim("") << "Total number of events selected: " << nsel << " [events]";
      LogVerbatim("") << "Selection Efficiency only W+(-): " << "(" << setprecision(4) << esel*100. <<" +/- "<< setprecision(2) << sqrt(esel*(1-esel)/ncharge)*100. << ")%";
     }
      LogVerbatim("") << ">>>>>> W SELECTION SUMMARY END   >>>>>>>>>>>>>>>\n";
}

bool WMuNuSelector::filter (Event & ev, const EventSetup &) {
      nall++;

      // Repeat Pre-Selection Cuts just in case...
      // Muon collection
      Handle<View<Muon> > muonCollection;
      if (!ev.getByToken(muonToken_, muonCollection)) {
            LogError("") << ">>> Muon collection does not exist !!!";
            return false;
      }
      unsigned int muonCollectionSize = muonCollection->size();

      // Trigger
      Handle<TriggerResults> triggerResults;
      if (!ev.getByToken(trigToken_, triggerResults)) {
            LogError("") << ">>> TRIGGER collection does not exist !!!";
            return false;
      }
      const edm::TriggerNames & triggerNames = ev.triggerNames(*triggerResults);
      bool trigger_fired = false;
      int itrig1 = triggerNames.triggerIndex(muonTrig_);
      if (triggerResults->accept(itrig1)) trigger_fired = true;
      LogTrace("") << ">>> Trigger bit: " << trigger_fired << " (" << muonTrig_ << ")";

      // Loop to reject/control Z->mumu is done separately
      unsigned int nmuonsForZ1 = 0;
      unsigned int nmuonsForZ2 = 0;
      for (unsigned int i=0; i<muonCollectionSize; i++) {
            const Muon& mu = muonCollection->at(i);
            if (!mu.isGlobalMuon()) continue;
            double pt = mu.pt();
            if (pt>ptThrForZ1_) nmuonsForZ1++;
            if (pt>ptThrForZ2_) nmuonsForZ2++;
      }
      LogTrace("") << "> Z rejection: muons above " << ptThrForZ1_ << " [GeV]: " << nmuonsForZ1;
      LogTrace("") << "> Z rejection: muons above " << ptThrForZ2_ << " [GeV]: " << nmuonsForZ2;

      // Jet collection
      Handle<View<Jet> > jetCollection;
      if (!ev.getByToken(jetToken_, jetCollection)) {
            LogError("") << ">>> JET collection does not exist !!!";
            return false;
      }
      unsigned int jetCollectionSize = jetCollection->size();
      int njets = 0;
      for (unsigned int i=0; i<jetCollectionSize; i++) {
            const Jet& jet = jetCollection->at(i);
            if (jet.et()>eJetMin_) njets++;
      }
      LogTrace("") << ">>> Total number of jets: " << jetCollectionSize;
      LogTrace("") << ">>> Number of jets above " << eJetMin_ << " [GeV]: " << njets;


      // Beam spot
      Handle<reco::BeamSpot> beamSpotHandle;
      if (!ev.getByToken(beamSpotToken_, beamSpotHandle)) {
            LogTrace("") << ">>> No beam spot found !!!";
            return false;
      }



      // Get WMuNu candidates from file:

      Handle<reco::WMuNuCandidateCollection> WMuNuCollection;
      if (!ev.getByToken(WMuNuCollectionToken_,WMuNuCollection) ) {
            LogTrace("") << ">>> WMuNu not found !!!";
            return false;
      }

      if(WMuNuCollection->empty()) {LogTrace("")<<"No WMuNu Candidates in the Event!"; return false;}
      if(WMuNuCollection->size() > 1) {LogTrace("")<<"This event contains more than one W Candidate";}

      // W->mu nu selection criteria

      LogTrace("") << "> WMuNu Candidate with: ";
      const WMuNuCandidate& WMuNu = WMuNuCollection->at(0);
      // WMuNuCandidates are ordered by Pt!
      // The Inclusive Selection WMuNu Candidate is the first one

      const reco::Muon & mu = WMuNu.getMuon();
      const reco::MET  & met =WMuNu.getNeutrino();
            if(plotHistograms_){
            h1_["hNWCand_sel"]->Fill(WMuNuCollection->size());
            }


      // Preselection cuts:

      if (!trigger_fired) {LogTrace("")<<"Event did not fire the Trigger"; return false;}
      ntrig++;

      if (nmuonsForZ1>=1 && nmuonsForZ2>=2) {LogTrace("")<<"Z Candidate!!"; return false;}
      if (njets>nJetMax_) {LogTrace("")<<"NJets > threshold";  return false;}

      npresel++;

      // Select Ws by charge:

      if (selectByCharge_*WMuNu.charge()==-1){ ncharge++;}


      // W->mu nu selection criteria

            if (!mu.isGlobalMuon()) return false;

            reco::TrackRef gm = mu.globalTrack();
            //reco::TrackRef tk = mu.innerTrack();

            // Pt,eta cuts
            double pt = mu.pt();
            double eta = mu.eta();
            LogTrace("") << "\t... Muon pt, eta: " << pt << " [GeV], " << eta;
                  if(plotHistograms_){ h1_["hPtMu_sel"]->Fill(pt);}
            if (pt<ptCut_) return false;
                  if(plotHistograms_){ h1_["hEtaMu_sel"]->Fill(eta);}
            if (fabs(eta)>etaCut_) return false;

            nkin++;


            // d0, chi2, nhits quality cuts
            double dxy = gm->dxy(beamSpotHandle->position());
            double normalizedChi2 = gm->normalizedChi2(); LogTrace("")<<"Im here"<<endl;
            double trackerHits = gm->hitPattern().numberOfValidTrackerHits();
            LogTrace("") << "\t... Muon dxy, normalizedChi2, trackerHits, isTrackerMuon?: " << dxy << " [cm], "<<normalizedChi2 << ", "<<trackerHits << ", " << mu.isTrackerMuon();

                  if(plotHistograms_){ h1_["hd0_sel"]->Fill(dxy);}
            if (!muon::isGoodMuon(mu,muon::GlobalMuonPromptTight) ) return false;
                  if(plotHistograms_){ h1_["hNormChi2_sel"]->Fill(normalizedChi2);}
            if (normalizedChi2>normalizedChi2Cut_) return false;
                  if(plotHistograms_){ h1_["hNHits_sel"]->Fill(trackerHits);}
            if (trackerHits<trackerHitsCut_) return false;
                  if(plotHistograms_){ h1_["hTracker_sel"]->Fill(mu.isTrackerMuon());}
            if (!mu.isTrackerMuon()) return false;

            nid++;

            // Acoplanarity cuts
            double acop = WMuNu.acop();
            LogTrace("") << "\t... acoplanarity: " << acop;

            // Isolation cuts
            double SumPt = mu.isolationR03().sumPt; double isovar=SumPt;
            double Cal   = mu.isolationR03().emEt + mu.isolationR03().hadEt; if(isCombinedIso_)isovar+=Cal;
                  if(plotHistograms_){
                  h1_["hPtSum_sel"]->Fill(SumPt);
                  h1_["hPtSumN_sel"]->Fill(SumPt/pt);
                  h1_["hCal_sel"]->Fill(Cal);
                  h1_["hIsoTot_sel"]->Fill( (SumPt+Cal));
                  h1_["hIsoTotN_sel"]->Fill( (SumPt+Cal) / pt );
                  }

            if (isRelativeIso_) isovar /= pt;
            bool iso = (isovar<=isoCut03_);
            LogTrace("") << "\t... isolation value" << isovar <<", isolated? " <<iso ;

            double met_et = met.pt();
            LogTrace("") << "\t... Met  pt: "<<WMuNu.getNeutrino().pt()<<"[GeV]";


            double massT = WMuNu.massT();
            double w_et = WMuNu.eT();

            LogTrace("") << "\t... W mass, W_et, W_px, W_py: " << massT << ", " << w_et << ", " << WMuNu.px() << ", " << WMuNu.py() << " [GeV]";


            // Plot 2D Histograms before final cuts
                  if(plotHistograms_ && acop<acopCut_){
                  h2_["hTMass_PtSum_inclusive"]->Fill(SumPt,massT);
                  h2_["hTMass_PtSumNorm_inclusive"]->Fill(SumPt/pt,massT);
                  h2_["hTMass_TotIsoNorm_inclusive"]->Fill((SumPt+Cal)/pt,massT);
                  h2_["hMET_PtSum_inclusive"]->Fill(SumPt,met_et);
                  h2_["hMET_PtSumNorm_inclusive"]->Fill(SumPt/pt,met_et);
                  h2_["hMET_TotIsoNorm_inclusive"]->Fill((SumPt+Cal)/pt,met_et);
                  }

            if (!iso) return false;

            niso++;

             if(plotHistograms_){ h1_["hAcop_sel"]->Fill(acop);}
            if (acop>=acopCut_) return false;

           nacop++;

            if(plotHistograms_){
                  h1_["hMET_sel"]->Fill(met_et);
                  h1_["hTMass_sel"]->Fill(massT);
            }


            if (massT<=mtMin_ || massT>=mtMax_)  return false;
            if (met_et<=metMin_ || met_et>=metMax_) return false;

           LogTrace("") << ">>>> Event ACCEPTED";
            nsel++;


            // (To be continued ;-) )

            return true;
}

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

DEFINE_FWK_MODULE(WMuNuSelector);