00001
00002
00003
00004
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00021
00022 #include "FWCore/Framework/interface/EDFilter.h"
00023 #include "FWCore/Utilities/interface/InputTag.h"
00024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00025 #include "FWCore/Framework/interface/Event.h"
00026 #include "TH1D.h"
00027 #include "TH2D.h"
00028
00029 class WMuNuSelector : public edm::EDFilter {
00030 public:
00031 WMuNuSelector (const edm::ParameterSet &);
00032 virtual bool filter(edm::Event&, const edm::EventSetup&);
00033 virtual void beginJob();
00034 virtual void endJob();
00035 void init_histograms();
00036 private:
00037 bool plotHistograms_;
00038 edm::InputTag trigTag_;
00039 edm::InputTag muonTag_;
00040 edm::InputTag jetTag_;
00041 edm::InputTag WMuNuCollectionTag_;
00042 const std::string muonTrig_;
00043 double ptThrForZ1_;
00044 double ptThrForZ2_;
00045 double eJetMin_;
00046 int nJetMax_;
00047 double ptCut_;
00048 double etaCut_;
00049 bool isRelativeIso_;
00050 bool isCombinedIso_;
00051 double isoCut03_;
00052 double mtMin_;
00053 double mtMax_;
00054 double metMin_;
00055 double metMax_;
00056 double acopCut_;
00057
00058 double dxyCut_;
00059 double normalizedChi2Cut_;
00060 int trackerHitsCut_;
00061 bool isAlsoTrackerMuon_;
00062
00063 int selectByCharge_;
00064
00065 double nall;
00066 double ntrig, npresel;
00067 double nsel;
00068 double ncharge;
00069 double nkin, nid,nacop,niso,nmass;
00070
00071
00072
00073 std::map<std::string,TH1D*> h1_;
00074 std::map<std::string,TH2D*> h2_;
00075
00076 };
00077 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00078 #include "DataFormats/Common/interface/Handle.h"
00079
00080 #include "FWCore/ServiceRegistry/interface/Service.h"
00081 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00082
00083 #include "DataFormats/TrackReco/interface/Track.h"
00084 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00085
00086 #include "DataFormats/MuonReco/interface/Muon.h"
00087 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
00088
00089 #include "DataFormats/METReco/interface/MET.h"
00090 #include "DataFormats/JetReco/interface/Jet.h"
00091
00092 #include "DataFormats/GeometryVector/interface/Phi.h"
00093
00094 #include "FWCore/Common/interface/TriggerNames.h"
00095 #include "DataFormats/Common/interface/TriggerResults.h"
00096
00097 #include "DataFormats/Common/interface/View.h"
00098
00099 #include "AnalysisDataFormats/EWK/interface/WMuNuCandidate.h"
00100
00101
00102 using namespace edm;
00103 using namespace std;
00104 using namespace reco;
00105
00106 WMuNuSelector::WMuNuSelector( const ParameterSet & cfg ) :
00107
00108 plotHistograms_(cfg.getUntrackedParameter<bool> ("plotHistograms", true)),
00109
00110
00111 trigTag_(cfg.getUntrackedParameter<edm::InputTag> ("TrigTag", edm::InputTag("TriggerResults::HLT"))),
00112 muonTag_(cfg.getUntrackedParameter<edm::InputTag> ("MuonTag", edm::InputTag("muons"))),
00113 jetTag_(cfg.getUntrackedParameter<edm::InputTag> ("JetTag", edm::InputTag("sisCone5CaloJets"))),
00114 WMuNuCollectionTag_(cfg.getUntrackedParameter<edm::InputTag> ("WMuNuCollectionTag", edm::InputTag("WMuNus"))),
00115
00116
00117
00118 muonTrig_(cfg.getUntrackedParameter<std::string>("MuonTrig", "HLT_Mu9")),
00119 ptThrForZ1_(cfg.getUntrackedParameter<double>("PtThrForZ1", 20.)),
00120 ptThrForZ2_(cfg.getUntrackedParameter<double>("PtThrForZ2", 10.)),
00121 eJetMin_(cfg.getUntrackedParameter<double>("EJetMin", 999999.)),
00122 nJetMax_(cfg.getUntrackedParameter<int>("NJetMax", 999999)),
00123
00124
00125
00126 ptCut_(cfg.getUntrackedParameter<double>("PtCut", 25.)),
00127 etaCut_(cfg.getUntrackedParameter<double>("EtaCut", 2.1)),
00128 isRelativeIso_(cfg.getUntrackedParameter<bool>("IsRelativeIso", true)),
00129 isCombinedIso_(cfg.getUntrackedParameter<bool>("IsCombinedIso", false)),
00130 isoCut03_(cfg.getUntrackedParameter<double>("IsoCut03", 0.1)),
00131 mtMin_(cfg.getUntrackedParameter<double>("MtMin", 50.)),
00132 mtMax_(cfg.getUntrackedParameter<double>("MtMax", 200.)),
00133 metMin_(cfg.getUntrackedParameter<double>("MetMin", -999999.)),
00134 metMax_(cfg.getUntrackedParameter<double>("MetMax", 999999.)),
00135 acopCut_(cfg.getUntrackedParameter<double>("AcopCut", 2.)),
00136
00137
00138 dxyCut_(cfg.getUntrackedParameter<double>("DxyCut", 0.2)),
00139 normalizedChi2Cut_(cfg.getUntrackedParameter<double>("NormalizedChi2Cut", 10.)),
00140 trackerHitsCut_(cfg.getUntrackedParameter<int>("TrackerHitsCut", 11)),
00141 isAlsoTrackerMuon_(cfg.getUntrackedParameter<bool>("IsAlsoTrackerMuon", true)),
00142
00143
00144 selectByCharge_(cfg.getUntrackedParameter<int>("SelectByCharge", 0))
00145 {
00146 }
00147
00148 void WMuNuSelector::beginJob() {
00149 nall = 0;
00150 ntrig=0;
00151 npresel=0;
00152 ncharge = 0;
00153 nkin=0;
00154 nid=0;
00155 nacop=0;
00156 niso=0;
00157 nsel = 0;
00158
00159 if(plotHistograms_){
00160 edm::Service<TFileService> fs;
00161 h1_["hNWCand_sel"] =fs->make<TH1D>("NWCand_sel","Nb. of WCandidates passing pre-selection (ordered by pt)",10,0.,10.);
00162 h1_["hPtMu_sel"] =fs->make<TH1D>("ptMu_sel","Pt mu",100,0.,100.);
00163 h1_["hEtaMu_sel"] =fs->make<TH1D>("etaMu_sel","Eta mu",50,-2.5,2.5);
00164 h1_["hd0_sel"] =fs->make<TH1D>("d0_sel","Impact parameter",1000,-1.,1.);
00165 h1_["hNHits_sel"] =fs->make<TH1D>("NumberOfValidHits_sel","Number of Hits in Silicon",100,0.,100.);
00166 h1_["hNormChi2_sel"] =fs->make<TH1D>("NormChi2_sel","Chi2/ndof of global track",1000,0.,50.);
00167 h1_["hTracker_sel"] =fs->make<TH1D>("isTrackerMuon_sel","is Tracker Muon?",2,0.,2.);
00168 h1_["hMET_sel"] =fs->make<TH1D>("MET_sel","Missing Transverse Energy (GeV)", 300,0.,300.);
00169 h1_["hTMass_sel"] =fs->make<TH1D>("TMass_sel","Rec. Transverse Mass (GeV)",300,0.,300.);
00170 h1_["hAcop_sel"] =fs->make<TH1D>("Acop_sel","Mu-MET acoplanarity",50,0.,M_PI);
00171 h1_["hPtSum_sel"] =fs->make<TH1D>("ptSum_sel","Track Isolation, Sum pT (GeV)",200,0.,100.);
00172 h1_["hPtSumN_sel"] =fs->make<TH1D>("ptSumN_sel","Track Isolation, Sum pT/pT",1000,0.,10.);
00173 h1_["hCal_sel"] =fs->make<TH1D>("Cal_sel","Calorimetric isolation, HCAL+ECAL (GeV)",200,0.,100.);
00174 h1_["hIsoTotN_sel"] =fs->make<TH1D>("isoTotN_sel","(Sum pT + Cal)/pT",1000,0.,10.);
00175 h1_["hIsoTot_sel"] =fs->make<TH1D>("isoTot_sel","(Sum pT + Cal)",200,0.,100.);
00176 h2_["hTMass_PtSum_inclusive"] =fs->make<TH2D>("TMass_PtSum_inclusive","Rec. Transverse Mass (GeV) vs Sum pT (GeV)",200,0.,100.,300,0.,300.);
00177 h2_["hTMass_PtSumNorm_inclusive"] =fs->make<TH2D>("TMass_PtSumNorm_inclusive","Rec. Transverse Mass (GeV) vs Sum Pt / Pt", 1000,0,10,300,0,300);
00178 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);
00179 h2_["hMET_PtSum_inclusive"] =fs->make<TH2D>("MET_PtSum_inclusive","Missing Transverse Energy (GeV) vs Sum Pt (GeV)",200,0.,100.,300,0.,300.);
00180 h2_["hMET_PtSumNorm_inclusive"] =fs->make<TH2D>("MET_PtSumNorm_inclusive","Missing Transverse Energy (GeV) vs Sum Pt/Pt",1000,0,10,300,0,300);
00181 h2_["hMET_TotIsoNorm_inclusive"]=fs->make<TH2D>("MET_TotIsoNorm_inclusive","Missing Transverse Energy (GeV) vs (SumPt + Cal)/Pt",10000,0,10,200,0,200);
00182 }
00183 }
00184
00185
00186 void WMuNuSelector::endJob() {
00187 double all = nall;
00188 double epresel = npresel/all;
00189 double etrig = ntrig/all;
00190 double ekin = nkin/all;
00191 double eid = nid/all;
00192 double eacop = nacop/all;
00193 double eiso = niso/all;
00194 double esel = nsel/all;
00195
00196 LogVerbatim("") << "\n>>>>>> W SELECTION SUMMARY BEGIN >>>>>>>>>>>>>>>";
00197 LogVerbatim("") << "Total number of events analyzed: " << nall << " [events]";
00198 LogVerbatim("") << "Total number of events triggered: " << ntrig << " [events]";
00199 LogVerbatim("") << "Total number of events pre-selected: " << npresel << " [events]";
00200 LogVerbatim("") << "Total number of events after kinematic cuts: " << nkin << " [events]";
00201 LogVerbatim("") << "Total number of events after Muon ID cuts: " << nid << " [events]";
00202 LogVerbatim("") << "Total number of events after Acop cut: " << nacop << " [events]";
00203 LogVerbatim("") << "Total number of events after iso cut: " << niso << " [events]";
00204 LogVerbatim("") << "Total number of events selected: " << nsel << " [events]";
00205 LogVerbatim("") << "Efficiencies:";
00206 LogVerbatim("") << "Trigger Efficiency: " << "(" << setprecision(4) << etrig*100. <<" +/- "<< setprecision(2) << sqrt(etrig*(1-etrig)/all)*100. << ")%";
00207 LogVerbatim("") << "Pre-Selection Efficiency: " << "(" << setprecision(4) << epresel*100. <<" +/- "<< setprecision(2) << sqrt(epresel*(1-epresel)/all)*100. << ")%";
00208 LogVerbatim("") << "Pt, Eta Selection Efficiency: " << "(" << setprecision(4) << ekin*100. <<" +/- "<< setprecision(2) << sqrt(ekin*(1-ekin)/all)*100. << ")%";
00209 LogVerbatim("") << "MuonID Efficiency: " << "(" << setprecision(4) << eid*100. <<" +/- "<< setprecision(2) << sqrt(eid*(1-eid)/all)*100. << ")%";
00210 LogVerbatim("") << "Acop Efficiency: " << "(" << setprecision(4) << eacop*100. <<" +/- "<< setprecision(2) << sqrt(eacop*(1-eacop)/all)*100. << ")%";
00211 LogVerbatim("") << "Iso Efficiency: " << "(" << setprecision(4) << eiso*100. <<" +/- "<< setprecision(2) << sqrt(eiso*(1-eiso)/all)*100. << ")%";
00212 LogVerbatim("") << "Selection Efficiency: " << "(" << setprecision(4) << esel*100. <<" +/- "<< setprecision(2) << sqrt(esel*(1-esel)/nall)*100. << ")%";
00213
00214 if ( fabs(selectByCharge_)==1 ){
00215 esel = nsel/ncharge;
00216 LogVerbatim("") << "\n>>>>>> W+(-) SELECTION >>>>>>>>>>>>>>>";
00217 LogVerbatim("") << "Total number of W+(-) events pre-selected: " << ncharge << " [events]";
00218 LogVerbatim("") << "Total number of events selected: " << nsel << " [events]";
00219 LogVerbatim("") << "Selection Efficiency only W+(-): " << "(" << setprecision(4) << esel*100. <<" +/- "<< setprecision(2) << sqrt(esel*(1-esel)/ncharge)*100. << ")%";
00220 }
00221 LogVerbatim("") << ">>>>>> W SELECTION SUMMARY END >>>>>>>>>>>>>>>\n";
00222 }
00223
00224 bool WMuNuSelector::filter (Event & ev, const EventSetup &) {
00225 nall++;
00226
00227
00228
00229 Handle<View<Muon> > muonCollection;
00230 if (!ev.getByLabel(muonTag_, muonCollection)) {
00231 LogError("") << ">>> Muon collection does not exist !!!";
00232 return 0;
00233 }
00234 unsigned int muonCollectionSize = muonCollection->size();
00235
00236
00237 Handle<TriggerResults> triggerResults;
00238 if (!ev.getByLabel(trigTag_, triggerResults)) {
00239 LogError("") << ">>> TRIGGER collection does not exist !!!";
00240 return 0;
00241 }
00242 const edm::TriggerNames & triggerNames = ev.triggerNames(*triggerResults);
00243 bool trigger_fired = false;
00244 int itrig1 = triggerNames.triggerIndex(muonTrig_);
00245 if (triggerResults->accept(itrig1)) trigger_fired = true;
00246 LogTrace("") << ">>> Trigger bit: " << trigger_fired << " (" << muonTrig_ << ")";
00247
00248
00249 unsigned int nmuonsForZ1 = 0;
00250 unsigned int nmuonsForZ2 = 0;
00251 for (unsigned int i=0; i<muonCollectionSize; i++) {
00252 const Muon& mu = muonCollection->at(i);
00253 if (!mu.isGlobalMuon()) continue;
00254 double pt = mu.pt();
00255 if (pt>ptThrForZ1_) nmuonsForZ1++;
00256 if (pt>ptThrForZ2_) nmuonsForZ2++;
00257 }
00258 LogTrace("") << "> Z rejection: muons above " << ptThrForZ1_ << " [GeV]: " << nmuonsForZ1;
00259 LogTrace("") << "> Z rejection: muons above " << ptThrForZ2_ << " [GeV]: " << nmuonsForZ2;
00260
00261
00262 Handle<View<Jet> > jetCollection;
00263 if (!ev.getByLabel(jetTag_, jetCollection)) {
00264 LogError("") << ">>> JET collection does not exist !!!";
00265 return 0;
00266 }
00267 unsigned int jetCollectionSize = jetCollection->size();
00268 int njets = 0;
00269 for (unsigned int i=0; i<jetCollectionSize; i++) {
00270 const Jet& jet = jetCollection->at(i);
00271 if (jet.et()>eJetMin_) njets++;
00272 }
00273 LogTrace("") << ">>> Total number of jets: " << jetCollectionSize;
00274 LogTrace("") << ">>> Number of jets above " << eJetMin_ << " [GeV]: " << njets;
00275
00276
00277
00278 Handle<reco::BeamSpot> beamSpotHandle;
00279 if (!ev.getByLabel(InputTag("offlineBeamSpot"), beamSpotHandle)) {
00280 LogTrace("") << ">>> No beam spot found !!!";
00281 return false;
00282 }
00283
00284
00285
00286
00287
00288 Handle<reco::WMuNuCandidateCollection> WMuNuCollection;
00289 if (!ev.getByLabel(WMuNuCollectionTag_,WMuNuCollection) ) {
00290 LogTrace("") << ">>> WMuNu not found !!!";
00291 return false;
00292 }
00293
00294 if(WMuNuCollection->size() < 1) {LogTrace("")<<"No WMuNu Candidates in the Event!"; return 0;}
00295 if(WMuNuCollection->size() > 1) {LogTrace("")<<"This event contains more than one W Candidate";}
00296
00297
00298
00299 LogTrace("") << "> WMuNu Candidate with: ";
00300 const WMuNuCandidate& WMuNu = WMuNuCollection->at(0);
00301
00302
00303
00304 const reco::Muon & mu = WMuNu.getMuon();
00305 const reco::MET & met =WMuNu.getNeutrino();
00306 if(plotHistograms_){
00307 h1_["hNWCand_sel"]->Fill(WMuNuCollection->size());
00308 }
00309
00310
00311
00312
00313 if (!trigger_fired) {LogTrace("")<<"Event did not fire the Trigger"; return 0;}
00314 ntrig++;
00315
00316 if (nmuonsForZ1>=1 && nmuonsForZ2>=2) {LogTrace("")<<"Z Candidate!!"; return 0;}
00317 if (njets>nJetMax_) {LogTrace("")<<"NJets > threshold"; return 0;}
00318
00319 npresel++;
00320
00321
00322
00323 if (selectByCharge_*WMuNu.charge()==-1){ ncharge++;}
00324
00325
00326
00327
00328 if (!mu.isGlobalMuon()) return 0;
00329
00330 reco::TrackRef gm = mu.globalTrack();
00331
00332
00333
00334 double pt = mu.pt();
00335 double eta = mu.eta();
00336 LogTrace("") << "\t... Muon pt, eta: " << pt << " [GeV], " << eta;
00337 if(plotHistograms_){ h1_["hPtMu_sel"]->Fill(pt);}
00338 if (pt<ptCut_) return 0;
00339 if(plotHistograms_){ h1_["hEtaMu_sel"]->Fill(eta);}
00340 if (fabs(eta)>etaCut_) return 0;
00341
00342 nkin++;
00343
00344
00345
00346 double dxy = gm->dxy(beamSpotHandle->position());
00347 double normalizedChi2 = gm->normalizedChi2(); LogTrace("")<<"Im here"<<endl;
00348 double trackerHits = gm->hitPattern().numberOfValidTrackerHits();
00349 LogTrace("") << "\t... Muon dxy, normalizedChi2, trackerHits, isTrackerMuon?: " << dxy << " [cm], "<<normalizedChi2 << ", "<<trackerHits << ", " << mu.isTrackerMuon();
00350
00351 if(plotHistograms_){ h1_["hd0_sel"]->Fill(dxy);}
00352 if (!muon::isGoodMuon(mu,muon::GlobalMuonPromptTight) ) return 0;
00353 if(plotHistograms_){ h1_["hNormChi2_sel"]->Fill(normalizedChi2);}
00354 if (normalizedChi2>normalizedChi2Cut_) return 0;
00355 if(plotHistograms_){ h1_["hNHits_sel"]->Fill(trackerHits);}
00356 if (trackerHits<trackerHitsCut_) return 0;
00357 if(plotHistograms_){ h1_["hTracker_sel"]->Fill(mu.isTrackerMuon());}
00358 if (!mu.isTrackerMuon()) return 0;
00359
00360 nid++;
00361
00362
00363 double acop = WMuNu.acop();
00364 LogTrace("") << "\t... acoplanarity: " << acop;
00365
00366
00367 double SumPt = mu.isolationR03().sumPt; double isovar=SumPt;
00368 double Cal = mu.isolationR03().emEt + mu.isolationR03().hadEt; if(isCombinedIso_)isovar+=Cal;
00369 if(plotHistograms_){
00370 h1_["hPtSum_sel"]->Fill(SumPt);
00371 h1_["hPtSumN_sel"]->Fill(SumPt/pt);
00372 h1_["hCal_sel"]->Fill(Cal);
00373 h1_["hIsoTot_sel"]->Fill( (SumPt+Cal));
00374 h1_["hIsoTotN_sel"]->Fill( (SumPt+Cal) / pt );
00375 }
00376
00377 if (isRelativeIso_) isovar /= pt;
00378 bool iso = (isovar<=isoCut03_);
00379 LogTrace("") << "\t... isolation value" << isovar <<", isolated? " <<iso ;
00380
00381 double met_et = met.pt();
00382 LogTrace("") << "\t... Met pt: "<<WMuNu.getNeutrino().pt()<<"[GeV]";
00383
00384
00385 double massT = WMuNu.massT();
00386 double w_et = WMuNu.eT();
00387
00388 LogTrace("") << "\t... W mass, W_et, W_px, W_py: " << massT << ", " << w_et << ", " << WMuNu.px() << ", " << WMuNu.py() << " [GeV]";
00389
00390
00391
00392 if(plotHistograms_ && acop<acopCut_){
00393 h2_["hTMass_PtSum_inclusive"]->Fill(SumPt,massT);
00394 h2_["hTMass_PtSumNorm_inclusive"]->Fill(SumPt/pt,massT);
00395 h2_["hTMass_TotIsoNorm_inclusive"]->Fill((SumPt+Cal)/pt,massT);
00396 h2_["hMET_PtSum_inclusive"]->Fill(SumPt,met_et);
00397 h2_["hMET_PtSumNorm_inclusive"]->Fill(SumPt/pt,met_et);
00398 h2_["hMET_TotIsoNorm_inclusive"]->Fill((SumPt+Cal)/pt,met_et);
00399 }
00400
00401 if (!iso) return 0;
00402
00403 niso++;
00404
00405 if(plotHistograms_){ h1_["hAcop_sel"]->Fill(acop);}
00406 if (acop>=acopCut_) return 0;
00407
00408 nacop++;
00409
00410 if(plotHistograms_){
00411 h1_["hMET_sel"]->Fill(met_et);
00412 h1_["hTMass_sel"]->Fill(massT);
00413 }
00414
00415
00416 if (massT<=mtMin_ || massT>=mtMax_) return 0;
00417 if (met_et<=metMin_ || met_et>=metMax_) return 0;
00418
00419 LogTrace("") << ">>>> Event ACCEPTED";
00420 nsel++;
00421
00422
00423
00424
00425 return true;
00426 }
00427
00428 #include "FWCore/Framework/interface/MakerMacros.h"
00429
00430 DEFINE_FWK_MODULE(WMuNuSelector);