CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
WMuNuSelector.cc
Go to the documentation of this file.
1 // //
3 // WMuNuSelector based on WMuNuCandidates //
4 // //
6 // //
7 // Filter of WMuNuCandidates for Analysis //
8 // --> From a WMuNuCandidate collection //
9 // --> Pre-Selection of events based in event cuts (trigger, Z rejection, ttbar rejection) //
10 // --> The Ws are selected from the highest pt muon in the event (applying the standard WMuNu Selection cuts) //
11 // //
12 // --> Be careful: if this Selector is used as a filter for further analysis you still have to make sure that //
13 // the W Candidate you use for your modules is the first one in the collection!! //
14 // //
15 // Optionally, plots selection variables sequentially after cuts, //
16 // and 2D histograms for background determination. //
17 // //
18 // For basic plots before & after cuts (without Candidate formalism), use WMuNuValidator.cc //
19 // //
21 
32 #include "TH1D.h"
33 #include "TH2D.h"
34 
35 class WMuNuSelector : public edm::EDFilter {
36 public:
38  virtual bool filter(edm::Event&, const edm::EventSetup&) override;
39  virtual void beginJob() override;
40  virtual void endJob() override;
41  void init_histograms();
42 private:
50  double ptThrForZ1_;
51  double ptThrForZ2_;
52  double eJetMin_;
53  int nJetMax_;
54  double ptCut_;
55  double etaCut_;
58  double isoCut03_;
59  double mtMin_;
60  double mtMax_;
61  double metMin_;
62  double metMax_;
63  double acopCut_;
64 
65  double dxyCut_;
69 
71 
72  double nall;
73  double ntrig, npresel;
74  double nsel;
75  double ncharge;
76  double nkin, nid,nacop,niso,nmass;
77 
78 
79 
80  std::map<std::string,TH1D*> h1_;
81  std::map<std::string,TH2D*> h2_;
82 
83 };
86 
89 
91 
93 
95 
97 
99 
100 
101 using namespace edm;
102 using namespace std;
103 using namespace reco;
104 
106  // Fast selection (no histograms)
107  plotHistograms_(cfg.getUntrackedParameter<bool> ("plotHistograms", true)),
108 
109  // Input collections
110  trigToken_(consumes<TriggerResults>(cfg.getUntrackedParameter<edm::InputTag> ("TrigTag", edm::InputTag("TriggerResults::HLT")))),
111  muonToken_(consumes<View<Muon> >(cfg.getUntrackedParameter<edm::InputTag> ("MuonTag", edm::InputTag("muons")))),
112  jetToken_(consumes<View<Jet> >(cfg.getUntrackedParameter<edm::InputTag> ("JetTag", edm::InputTag("sisCone5CaloJets")))),
113  beamSpotToken_(consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"))),
114  WMuNuCollectionToken_(consumes<reco::WMuNuCandidateCollection>(cfg.getUntrackedParameter<edm::InputTag> ("WMuNuCollectionTag", edm::InputTag("WMuNus")))),
115 
116 
117  // Preselection cuts
118  muonTrig_(cfg.getUntrackedParameter<std::string>("MuonTrig", "HLT_Mu9")),
119  ptThrForZ1_(cfg.getUntrackedParameter<double>("PtThrForZ1", 20.)),
120  ptThrForZ2_(cfg.getUntrackedParameter<double>("PtThrForZ2", 10.)),
121  eJetMin_(cfg.getUntrackedParameter<double>("EJetMin", 999999.)),
122  nJetMax_(cfg.getUntrackedParameter<int>("NJetMax", 999999)),
123 
124 
125  // Main cuts
126  ptCut_(cfg.getUntrackedParameter<double>("PtCut", 25.)),
127  etaCut_(cfg.getUntrackedParameter<double>("EtaCut", 2.1)),
128  isRelativeIso_(cfg.getUntrackedParameter<bool>("IsRelativeIso", true)),
129  isCombinedIso_(cfg.getUntrackedParameter<bool>("IsCombinedIso", false)),
130  isoCut03_(cfg.getUntrackedParameter<double>("IsoCut03", 0.1)),
131  mtMin_(cfg.getUntrackedParameter<double>("MtMin", 50.)),
132  mtMax_(cfg.getUntrackedParameter<double>("MtMax", 200.)),
133  metMin_(cfg.getUntrackedParameter<double>("MetMin", -999999.)),
134  metMax_(cfg.getUntrackedParameter<double>("MetMax", 999999.)),
135  acopCut_(cfg.getUntrackedParameter<double>("AcopCut", 2.)),
136 
137  // Muon quality cuts
138  dxyCut_(cfg.getUntrackedParameter<double>("DxyCut", 0.2)),
139  normalizedChi2Cut_(cfg.getUntrackedParameter<double>("NormalizedChi2Cut", 10.)),
140  trackerHitsCut_(cfg.getUntrackedParameter<int>("TrackerHitsCut", 11)),
141  isAlsoTrackerMuon_(cfg.getUntrackedParameter<bool>("IsAlsoTrackerMuon", true)),
142 
143  // W+/W- Selection
144  selectByCharge_(cfg.getUntrackedParameter<int>("SelectByCharge", 0))
145 {
146 }
147 
149  nall = 0;
150  ntrig=0;
151  npresel=0;
152  ncharge = 0;
153  nkin=0;
154  nid=0;
155  nacop=0;
156  niso=0;
157  nsel = 0;
158 
159  if(plotHistograms_){
161  h1_["hNWCand_sel"] =fs->make<TH1D>("NWCand_sel","Nb. of WCandidates passing pre-selection (ordered by pt)",10,0.,10.);
162  h1_["hPtMu_sel"] =fs->make<TH1D>("ptMu_sel","Pt mu",100,0.,100.);
163  h1_["hEtaMu_sel"] =fs->make<TH1D>("etaMu_sel","Eta mu",50,-2.5,2.5);
164  h1_["hd0_sel"] =fs->make<TH1D>("d0_sel","Impact parameter",1000,-1.,1.);
165  h1_["hNHits_sel"] =fs->make<TH1D>("NumberOfValidHits_sel","Number of Hits in Silicon",100,0.,100.);
166  h1_["hNormChi2_sel"] =fs->make<TH1D>("NormChi2_sel","Chi2/ndof of global track",1000,0.,50.);
167  h1_["hTracker_sel"] =fs->make<TH1D>("isTrackerMuon_sel","is Tracker Muon?",2,0.,2.);
168  h1_["hMET_sel"] =fs->make<TH1D>("MET_sel","Missing Transverse Energy (GeV)", 300,0.,300.);
169  h1_["hTMass_sel"] =fs->make<TH1D>("TMass_sel","Rec. Transverse Mass (GeV)",300,0.,300.);
170  h1_["hAcop_sel"] =fs->make<TH1D>("Acop_sel","Mu-MET acoplanarity",50,0.,M_PI);
171  h1_["hPtSum_sel"] =fs->make<TH1D>("ptSum_sel","Track Isolation, Sum pT (GeV)",200,0.,100.);
172  h1_["hPtSumN_sel"] =fs->make<TH1D>("ptSumN_sel","Track Isolation, Sum pT/pT",1000,0.,10.);
173  h1_["hCal_sel"] =fs->make<TH1D>("Cal_sel","Calorimetric isolation, HCAL+ECAL (GeV)",200,0.,100.);
174  h1_["hIsoTotN_sel"] =fs->make<TH1D>("isoTotN_sel","(Sum pT + Cal)/pT",1000,0.,10.);
175  h1_["hIsoTot_sel"] =fs->make<TH1D>("isoTot_sel","(Sum pT + Cal)",200,0.,100.);
176  h2_["hTMass_PtSum_inclusive"] =fs->make<TH2D>("TMass_PtSum_inclusive","Rec. Transverse Mass (GeV) vs Sum pT (GeV)",200,0.,100.,300,0.,300.);
177  h2_["hTMass_PtSumNorm_inclusive"] =fs->make<TH2D>("TMass_PtSumNorm_inclusive","Rec. Transverse Mass (GeV) vs Sum Pt / Pt", 1000,0,10,300,0,300);
178  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);
179  h2_["hMET_PtSum_inclusive"] =fs->make<TH2D>("MET_PtSum_inclusive","Missing Transverse Energy (GeV) vs Sum Pt (GeV)",200,0.,100.,300,0.,300.);
180  h2_["hMET_PtSumNorm_inclusive"] =fs->make<TH2D>("MET_PtSumNorm_inclusive","Missing Transverse Energy (GeV) vs Sum Pt/Pt",1000,0,10,300,0,300);
181  h2_["hMET_TotIsoNorm_inclusive"]=fs->make<TH2D>("MET_TotIsoNorm_inclusive","Missing Transverse Energy (GeV) vs (SumPt + Cal)/Pt",10000,0,10,200,0,200);
182  }
183 }
184 
185 
187  double all = nall;
188  double epresel = npresel/all;
189  double etrig = ntrig/all;
190  double ekin = nkin/all;
191  double eid = nid/all;
192  double eacop = nacop/all;
193  double eiso = niso/all;
194  double esel = nsel/all;
195 
196  LogVerbatim("") << "\n>>>>>> W SELECTION SUMMARY BEGIN >>>>>>>>>>>>>>>";
197  LogVerbatim("") << "Total number of events analyzed: " << nall << " [events]";
198  LogVerbatim("") << "Total number of events triggered: " << ntrig << " [events]";
199  LogVerbatim("") << "Total number of events pre-selected: " << npresel << " [events]";
200  LogVerbatim("") << "Total number of events after kinematic cuts: " << nkin << " [events]";
201  LogVerbatim("") << "Total number of events after Muon ID cuts: " << nid << " [events]";
202  LogVerbatim("") << "Total number of events after Acop cut: " << nacop << " [events]";
203  LogVerbatim("") << "Total number of events after iso cut: " << niso << " [events]";
204  LogVerbatim("") << "Total number of events selected: " << nsel << " [events]";
205  LogVerbatim("") << "Efficiencies:";
206  LogVerbatim("") << "Trigger Efficiency: " << "(" << setprecision(4) << etrig*100. <<" +/- "<< setprecision(2) << sqrt(etrig*(1-etrig)/all)*100. << ")%";
207  LogVerbatim("") << "Pre-Selection Efficiency: " << "(" << setprecision(4) << epresel*100. <<" +/- "<< setprecision(2) << sqrt(epresel*(1-epresel)/all)*100. << ")%";
208  LogVerbatim("") << "Pt, Eta Selection Efficiency: " << "(" << setprecision(4) << ekin*100. <<" +/- "<< setprecision(2) << sqrt(ekin*(1-ekin)/all)*100. << ")%";
209  LogVerbatim("") << "MuonID Efficiency: " << "(" << setprecision(4) << eid*100. <<" +/- "<< setprecision(2) << sqrt(eid*(1-eid)/all)*100. << ")%";
210  LogVerbatim("") << "Acop Efficiency: " << "(" << setprecision(4) << eacop*100. <<" +/- "<< setprecision(2) << sqrt(eacop*(1-eacop)/all)*100. << ")%";
211  LogVerbatim("") << "Iso Efficiency: " << "(" << setprecision(4) << eiso*100. <<" +/- "<< setprecision(2) << sqrt(eiso*(1-eiso)/all)*100. << ")%";
212  LogVerbatim("") << "Selection Efficiency: " << "(" << setprecision(4) << esel*100. <<" +/- "<< setprecision(2) << sqrt(esel*(1-esel)/nall)*100. << ")%";
213 
214  if ( fabs(selectByCharge_)==1 ){
215  esel = nsel/ncharge;
216  LogVerbatim("") << "\n>>>>>> W+(-) SELECTION >>>>>>>>>>>>>>>";
217  LogVerbatim("") << "Total number of W+(-) events pre-selected: " << ncharge << " [events]";
218  LogVerbatim("") << "Total number of events selected: " << nsel << " [events]";
219  LogVerbatim("") << "Selection Efficiency only W+(-): " << "(" << setprecision(4) << esel*100. <<" +/- "<< setprecision(2) << sqrt(esel*(1-esel)/ncharge)*100. << ")%";
220  }
221  LogVerbatim("") << ">>>>>> W SELECTION SUMMARY END >>>>>>>>>>>>>>>\n";
222 }
223 
225  nall++;
226 
227  // Repeat Pre-Selection Cuts just in case...
228  // Muon collection
229  Handle<View<Muon> > muonCollection;
230  if (!ev.getByToken(muonToken_, muonCollection)) {
231  LogError("") << ">>> Muon collection does not exist !!!";
232  return 0;
233  }
234  unsigned int muonCollectionSize = muonCollection->size();
235 
236  // Trigger
238  if (!ev.getByToken(trigToken_, triggerResults)) {
239  LogError("") << ">>> TRIGGER collection does not exist !!!";
240  return 0;
241  }
242  const edm::TriggerNames & triggerNames = ev.triggerNames(*triggerResults);
243  bool trigger_fired = false;
244  int itrig1 = triggerNames.triggerIndex(muonTrig_);
245  if (triggerResults->accept(itrig1)) trigger_fired = true;
246  LogTrace("") << ">>> Trigger bit: " << trigger_fired << " (" << muonTrig_ << ")";
247 
248  // Loop to reject/control Z->mumu is done separately
249  unsigned int nmuonsForZ1 = 0;
250  unsigned int nmuonsForZ2 = 0;
251  for (unsigned int i=0; i<muonCollectionSize; i++) {
252  const Muon& mu = muonCollection->at(i);
253  if (!mu.isGlobalMuon()) continue;
254  double pt = mu.pt();
255  if (pt>ptThrForZ1_) nmuonsForZ1++;
256  if (pt>ptThrForZ2_) nmuonsForZ2++;
257  }
258  LogTrace("") << "> Z rejection: muons above " << ptThrForZ1_ << " [GeV]: " << nmuonsForZ1;
259  LogTrace("") << "> Z rejection: muons above " << ptThrForZ2_ << " [GeV]: " << nmuonsForZ2;
260 
261  // Jet collection
262  Handle<View<Jet> > jetCollection;
263  if (!ev.getByToken(jetToken_, jetCollection)) {
264  LogError("") << ">>> JET collection does not exist !!!";
265  return 0;
266  }
267  unsigned int jetCollectionSize = jetCollection->size();
268  int njets = 0;
269  for (unsigned int i=0; i<jetCollectionSize; i++) {
270  const Jet& jet = jetCollection->at(i);
271  if (jet.et()>eJetMin_) njets++;
272  }
273  LogTrace("") << ">>> Total number of jets: " << jetCollectionSize;
274  LogTrace("") << ">>> Number of jets above " << eJetMin_ << " [GeV]: " << njets;
275 
276 
277  // Beam spot
278  Handle<reco::BeamSpot> beamSpotHandle;
279  if (!ev.getByToken(beamSpotToken_, beamSpotHandle)) {
280  LogTrace("") << ">>> No beam spot found !!!";
281  return false;
282  }
283 
284 
285 
286  // Get WMuNu candidates from file:
287 
289  if (!ev.getByToken(WMuNuCollectionToken_,WMuNuCollection) ) {
290  LogTrace("") << ">>> WMuNu not found !!!";
291  return false;
292  }
293 
294  if(WMuNuCollection->size() < 1) {LogTrace("")<<"No WMuNu Candidates in the Event!"; return 0;}
295  if(WMuNuCollection->size() > 1) {LogTrace("")<<"This event contains more than one W Candidate";}
296 
297  // W->mu nu selection criteria
298 
299  LogTrace("") << "> WMuNu Candidate with: ";
300  const WMuNuCandidate& WMuNu = WMuNuCollection->at(0);
301  // WMuNuCandidates are ordered by Pt!
302  // The Inclusive Selection WMuNu Candidate is the first one
303 
304  const reco::Muon & mu = WMuNu.getMuon();
305  const reco::MET & met =WMuNu.getNeutrino();
306  if(plotHistograms_){
307  h1_["hNWCand_sel"]->Fill(WMuNuCollection->size());
308  }
309 
310 
311  // Preselection cuts:
312 
313  if (!trigger_fired) {LogTrace("")<<"Event did not fire the Trigger"; return 0;}
314  ntrig++;
315 
316  if (nmuonsForZ1>=1 && nmuonsForZ2>=2) {LogTrace("")<<"Z Candidate!!"; return 0;}
317  if (njets>nJetMax_) {LogTrace("")<<"NJets > threshold"; return 0;}
318 
319  npresel++;
320 
321  // Select Ws by charge:
322 
323  if (selectByCharge_*WMuNu.charge()==-1){ ncharge++;}
324 
325 
326  // W->mu nu selection criteria
327 
328  if (!mu.isGlobalMuon()) return 0;
329 
330  reco::TrackRef gm = mu.globalTrack();
331  //reco::TrackRef tk = mu.innerTrack();
332 
333  // Pt,eta cuts
334  double pt = mu.pt();
335  double eta = mu.eta();
336  LogTrace("") << "\t... Muon pt, eta: " << pt << " [GeV], " << eta;
337  if(plotHistograms_){ h1_["hPtMu_sel"]->Fill(pt);}
338  if (pt<ptCut_) return 0;
339  if(plotHistograms_){ h1_["hEtaMu_sel"]->Fill(eta);}
340  if (fabs(eta)>etaCut_) return 0;
341 
342  nkin++;
343 
344 
345  // d0, chi2, nhits quality cuts
346  double dxy = gm->dxy(beamSpotHandle->position());
347  double normalizedChi2 = gm->normalizedChi2(); LogTrace("")<<"Im here"<<endl;
348  double trackerHits = gm->hitPattern().numberOfValidTrackerHits();
349  LogTrace("") << "\t... Muon dxy, normalizedChi2, trackerHits, isTrackerMuon?: " << dxy << " [cm], "<<normalizedChi2 << ", "<<trackerHits << ", " << mu.isTrackerMuon();
350 
351  if(plotHistograms_){ h1_["hd0_sel"]->Fill(dxy);}
352  if (!muon::isGoodMuon(mu,muon::GlobalMuonPromptTight) ) return 0;
353  if(plotHistograms_){ h1_["hNormChi2_sel"]->Fill(normalizedChi2);}
354  if (normalizedChi2>normalizedChi2Cut_) return 0;
355  if(plotHistograms_){ h1_["hNHits_sel"]->Fill(trackerHits);}
356  if (trackerHits<trackerHitsCut_) return 0;
357  if(plotHistograms_){ h1_["hTracker_sel"]->Fill(mu.isTrackerMuon());}
358  if (!mu.isTrackerMuon()) return 0;
359 
360  nid++;
361 
362  // Acoplanarity cuts
363  double acop = WMuNu.acop();
364  LogTrace("") << "\t... acoplanarity: " << acop;
365 
366  // Isolation cuts
367  double SumPt = mu.isolationR03().sumPt; double isovar=SumPt;
368  double Cal = mu.isolationR03().emEt + mu.isolationR03().hadEt; if(isCombinedIso_)isovar+=Cal;
369  if(plotHistograms_){
370  h1_["hPtSum_sel"]->Fill(SumPt);
371  h1_["hPtSumN_sel"]->Fill(SumPt/pt);
372  h1_["hCal_sel"]->Fill(Cal);
373  h1_["hIsoTot_sel"]->Fill( (SumPt+Cal));
374  h1_["hIsoTotN_sel"]->Fill( (SumPt+Cal) / pt );
375  }
376 
377  if (isRelativeIso_) isovar /= pt;
378  bool iso = (isovar<=isoCut03_);
379  LogTrace("") << "\t... isolation value" << isovar <<", isolated? " <<iso ;
380 
381  double met_et = met.pt();
382  LogTrace("") << "\t... Met pt: "<<WMuNu.getNeutrino().pt()<<"[GeV]";
383 
384 
385  double massT = WMuNu.massT();
386  double w_et = WMuNu.eT();
387 
388  LogTrace("") << "\t... W mass, W_et, W_px, W_py: " << massT << ", " << w_et << ", " << WMuNu.px() << ", " << WMuNu.py() << " [GeV]";
389 
390 
391  // Plot 2D Histograms before final cuts
392  if(plotHistograms_ && acop<acopCut_){
393  h2_["hTMass_PtSum_inclusive"]->Fill(SumPt,massT);
394  h2_["hTMass_PtSumNorm_inclusive"]->Fill(SumPt/pt,massT);
395  h2_["hTMass_TotIsoNorm_inclusive"]->Fill((SumPt+Cal)/pt,massT);
396  h2_["hMET_PtSum_inclusive"]->Fill(SumPt,met_et);
397  h2_["hMET_PtSumNorm_inclusive"]->Fill(SumPt/pt,met_et);
398  h2_["hMET_TotIsoNorm_inclusive"]->Fill((SumPt+Cal)/pt,met_et);
399  }
400 
401  if (!iso) return 0;
402 
403  niso++;
404 
405  if(plotHistograms_){ h1_["hAcop_sel"]->Fill(acop);}
406  if (acop>=acopCut_) return 0;
407 
408  nacop++;
409 
410  if(plotHistograms_){
411  h1_["hMET_sel"]->Fill(met_et);
412  h1_["hTMass_sel"]->Fill(massT);
413  }
414 
415 
416  if (massT<=mtMin_ || massT>=mtMax_) return 0;
417  if (met_et<=metMin_ || met_et>=metMax_) return 0;
418 
419  LogTrace("") << ">>>> Event ACCEPTED";
420  nsel++;
421 
422 
423  // (To be continued ;-) )
424 
425  return true;
426 }
427 
429 
float hadEt
hcal sum-Et
Definition: MuonIsolation.h:9
int i
Definition: DBlmapReader.cc:9
const std::string muonTrig_
virtual edm::TriggerNames const & triggerNames(edm::TriggerResults const &triggerResults) const
Definition: Event.cc:208
tuple cfg
Definition: looper.py:259
edm::EDGetTokenT< edm::View< reco::Muon > > muonToken_
float sumPt
sum-pt of tracks
Definition: MuonIsolation.h:7
virtual double et() const
transverse energy
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:449
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
edm::EDGetTokenT< edm::TriggerResults > trigToken_
bool isTrackerMuon() const
Definition: Muon.h:219
Base class for all types of Jets.
Definition: Jet.h:20
bool isAlsoTrackerMuon_
bool isGlobalMuon() const
Definition: Muon.h:218
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
bool ev
double ptThrForZ1_
T eta() const
edm::EDGetTokenT< edm::View< reco::Jet > > jetToken_
virtual double eta() const
momentum pseudorapidity
virtual double pt() const
transverse momentum
edm::EDGetTokenT< reco::BeamSpot > beamSpotToken_
unsigned int triggerIndex(std::string const &name) const
Definition: TriggerNames.cc:32
Definition: MET.h:42
virtual void beginJob() override
T sqrt(T t)
Definition: SSEVec.h:48
double eT() const
float emEt
ecal sum-Et
Definition: MuonIsolation.h:8
virtual int charge() const
electric charge
Definition: LeafCandidate.h:91
double acop() const
void init_histograms()
const int mu
Definition: Constants.h:22
static std::string const triggerResults
Definition: EdmProvDump.cc:40
virtual void endJob() override
std::map< std::string, TH1D * > h1_
#define LogTrace(id)
#define M_PI
bool isGoodMuon(const reco::Muon &muon, SelectionType type, reco::Muon::ArbitrationType arbitrationType=reco::Muon::SegmentAndTrackArbitration)
main GoodMuon wrapper call
std::vector< reco::WMuNuCandidate > WMuNuCandidateCollection
virtual double px() const
x coordinate of momentum vector
WMuNuSelector(const edm::ParameterSet &)
double massT() const
edm::EDGetTokenT< reco::WMuNuCandidateCollection > WMuNuCollectionToken_
std::map< std::string, TH2D * > h2_
double ptThrForZ2_
volatile std::atomic< bool > shutdown_flag false
double normalizedChi2Cut_
virtual bool filter(edm::Event &, const edm::EventSetup &) override
const reco::Muon & getMuon() const
const reco::MET & getNeutrino() const
virtual double py() const
y coordinate of momentum vector
const MuonIsolation & isolationR03() const
Definition: Muon.h:158
virtual TrackRef globalTrack() const
reference to Track reconstructed in both tracked and muon detector
Definition: Muon.h:54