ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/LJMet/Utils/src/CutFlow.cc
Revision: 1.5
Committed: Tue Nov 10 22:29:35 2009 UTC (15 years, 5 months ago) by kukartse
Content type: text/plain
Branch: MAIN
CVS Tags: V00-00-02, V00-00-01_BeamSplash09, HEAD
Changes since 1.4: +2 -2 lines
Log Message:
analyzer for checking RECO experimental data

File Contents

# Content
1 // -*- C++ -*-
2 //
3 // Package: CutFlow
4 // Class: CutFlow
5 //
6 /**\class CutFlow CutFlow.h LJMet/Utils/interface/CutFlow.h
7
8 Description: times various aspects of data processing
9
10 Implementation:
11 Keep this class as simple and fast as possible
12 */
13 //
14 // Original Author: "Gennadiy Kukartsev"
15 // Created: Tue Oct 6 13:39:12 CDT 2009
16 // $Id: CutFlow.cc,v 1.4 2009/10/27 22:39:44 kukartse Exp $
17 //
18 //
19
20 #include "LJMet/Utils/interface/CutFlow.h"
21 #include "FWCore/Framework/interface/TriggerNames.h"
22 #include "DataFormats/Common/interface/TriggerResults.h"
23 #include "DataFormats/PatCandidates/interface/Jet.h"
24 #include "DataFormats/PatCandidates/interface/Electron.h"
25 #include "DataFormats/PatCandidates/interface/Muon.h"
26 #include "DataFormats/PatCandidates/interface/MET.h"
27 #include "DataFormats/PatCandidates/interface/TriggerPath.h"
28 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
29 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
30 #include "DataFormats/Math/interface/deltaR.h"
31
32 using namespace edm;
33
34 CutFlow::CutFlow(const edm::ParameterSet& iConfig){
35 c_total = new MeanCounter("Processing event number: ", 0, 1);
36 c_1 = new MeanCounter("CutFlow stage 1: ", 0, 1);
37 c_2 = new MeanCounter("CutFlow stage 2: ", 0, 1);
38 c_3 = new MeanCounter("CutFlow stage 3: ", 0, 1);
39 c_4 = new MeanCounter("CutFlow stage 4: ", 0, 1);
40 c_5 = new MeanCounter("CutFlow stage 5: ", 0, 1);
41 c_1->setPrintCount(false);
42 c_2->setPrintCount(false);
43 c_3->setPrintCount(false);
44 c_4->setPrintCount(false);
45 c_5->setPrintCount(false);
46
47 muon_index = -1;
48 electron_index = -1;
49 }
50
51
52 CutFlow::~CutFlow()
53 {
54 delete c_total;
55 delete c_1;
56 delete c_2;
57 delete c_3;
58 delete c_4;
59 delete c_5;
60 }
61
62
63 void
64 CutFlow::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
65 {
66 c_total->count();
67 //
68 //_____ mu+jets cut flow
69 //
70 /*
71 while(1){
72 //
73 //_____ stage 1
74 //
75 if ( is_trigger("HLT_Mu9", iEvent) ){
76 }
77 else break;
78 if ( has_global_muon("selectedLayer1Muons", iEvent) ){
79 c_1->count();
80 }
81 else break;
82
83 if ( is_muon("selectedLayer1Muons",
84 20.0,
85 2.1,
86 11,
87 0.02,
88 10,
89 4.0,
90 6.0,
91 0.05,
92 iEvent) ){
93 c_2->count();
94 }
95 else break;
96 if ( is_jets("selectedLayer1Jets",
97 30.0,
98 2.4,
99 iEvent) ){
100 c_3->count();
101 }
102 else break;
103 if ( no_loose_muon("selectedLayer1Muons",
104 2.5,
105 10.0,
106 0.2,
107 iEvent) ){
108 c_4->count();
109 }
110 else break;
111 if ( no_loose_electron("selectedLayer1Electrons",
112 2.5,
113 15.0,
114 0.2,
115 iEvent) ){
116 c_5->count();
117 }
118 else break;
119 cout << "Cut Flow: passed stage 5 " << endl;
120 break;
121 }
122 */
123 //
124 //_____ e+jets cut flow
125 //
126 while(1){
127 //
128 //_____ stage 1
129 //
130 if ( is_trigger("HLT_Ele15_LW_L1R", iEvent) ){
131 }
132 else break;
133 if ( has_electron("selectedLayer1Electrons", iEvent) ){
134 c_1->count();
135 }
136 else break;
137 //
138 //_____ stage 2
139 //
140 cout << "CutFlow: passed stage 1: ";
141 if ( is_electron("selectedLayer1Electrons",
142 30.0,
143 2.4,
144 0.02,
145 0.1,
146 iEvent) ){
147 c_2->count();
148 cout << endl;
149 }
150 else{
151 cout << endl;
152 break;
153 }
154 //
155 //_____ stage 3
156 //
157 cout << "CutFlow: passed stage 2: ";
158 if ( !is_muon("selectedLayer1Muons",
159 30.0,
160 2.1,
161 11,
162 0.02,
163 10,
164 4000000.0,
165 6000000.0,
166 0.05,
167 iEvent) ){
168 c_3->count();
169 cout << endl;
170 }
171 else{
172 cout << endl;
173 break;
174 }
175 //
176 //_____ stage 4
177 //
178 cout << "CutFlow: passed stage 3: ";
179 int n_jets = count_jets("selectedLayer1Jets",
180 30.0,
181 2.4,
182 "selectedLayer1Electrons",
183 30.0,
184 2.4,
185 0.02,
186 0.3,
187 iEvent);
188 if ( n_jets >= 0){
189 c_4->count();
190 cout << endl;
191 }
192 else{
193 cout << endl;
194 break;
195 }
196 //
197 //_____ stage 5
198 //
199 cout << "CutFlow: passed stage 4: ";
200 if ( n_jets >= 4 ){
201 c_5->count();
202 cout << endl;
203 }
204 else{
205 cout << endl;
206 break;
207 }
208 cout << "Cut Flow: passed stage 5 " << endl;
209 break;
210 }
211 }
212
213
214 void
215 CutFlow::beginJob()
216 {
217 }
218
219
220 void
221 CutFlow::endJob() {
222 cout << "========= Cut Flow Report ==========================================" << endl;
223 cout << "Total events: " << c_total->getCount() << endl;
224 cout << "Pass stage 1: " << c_1->getCount() << endl;
225 cout << "Pass stage 2: " << c_2->getCount() << endl;
226 cout << "Pass stage 3: " << c_3->getCount() << endl;
227 cout << "Pass stage 4: " << c_4->getCount() << endl;
228 cout << "Pass stage 5: " << c_5->getCount() << endl;
229 cout << "================================================================" << endl;
230 }
231
232 bool CutFlow::is_trigger(std::string mTriggerName, const edm::Event& iEvent){
233 bool result = false;
234 Handle<std::vector<pat::TriggerPath> > trigs;
235 iEvent.getByLabel("patTrigger", trigs);
236 for (std::vector<pat::TriggerPath>::const_iterator trig = trigs->begin();
237 trig!=trigs->end();
238 ++trig){
239 //cout << "Trigger name and number: " << trig->name() << ", " << trig->index() << endl;
240 if( trig->name().find(mTriggerName)!=std::string::npos ){
241 result = trig->wasAccept();
242 }
243 }
244 return result;
245 }
246
247 bool CutFlow::has_global_muon(std::string mLabel,
248 const edm::Event& iEvent){
249 bool result = false;
250 //cout << "Cut Flow: passed trigger ";
251 Handle< vector< pat::Muon > > muons;
252 iEvent . getByLabel( mLabel, muons );
253 //
254 //_____ loop over muons _______________________________________________
255 //
256 MeanCounter passed_muons("");
257 for ( vector<pat::Muon>::const_iterator mu = muons -> begin(); mu != muons -> end(); mu++){
258 if ( mu->isGlobalMuon() ){
259 passed_muons.count();
260 }
261 }
262 result = passed_muons.getCount()>0;
263 //cout << endl;
264 return result;
265 }
266
267
268 bool CutFlow::is_muon(std::string mLabel,
269 double mPt,
270 double mEta,
271 int mNHits,
272 double mD0,
273 double mChi2Ndof,
274 double mEmVeto,
275 double mHadVeto,
276 double mRelIso,
277 const edm::Event& iEvent){
278 bool result = false;
279 //cout << "Cut Flow: passed stage 1 (mu+jets) or stage 2 (e+jets): ";
280 Handle< vector< pat::Muon > > muons;
281 iEvent . getByLabel( "selectedLayer1Muons", muons );
282 Handle< reco::BeamSpot > beamSpotHandle;
283 iEvent . getByLabel( "offlineBeamSpot", beamSpotHandle);
284 //
285 //_____ Beam spot _____________________________________________________
286 //
287 reco::BeamSpot beamSpot;
288 if ( beamSpotHandle.isValid() ){
289 beamSpot = *beamSpotHandle;
290 }
291 else{
292 edm::LogInfo("TtLJetsAnalyzer")
293 << "No beam spot available from EventSetup \n";
294 }
295 const reco::BeamSpot::Point _bs(beamSpot.x0(), beamSpot.y0(),beamSpot.z0());
296 //
297 //_____ loop over muons _______________________________________________
298 //
299 double _pt = 0.0;
300 double _abseta = 0.0;
301 int _nhits = 0;
302 double _d0 = 0.0;
303 double _chi2ndof = 0.0;
304 double _emveto = 10000.0;
305 double _hadveto = 10000.0;
306 double _trackiso = 10000.0;
307 double _caloiso = 10000.0;
308 double _reliso = 10000.0;
309 MeanCounter passed_muons("");
310 for ( vector<pat::Muon>::const_iterator mu = muons -> begin(); mu != muons -> end(); mu++){
311 //_____ check that the muon is global
312 if( !(mu->isGlobalMuon()) ) continue;
313 _pt = mu->pt();
314 _abseta = fabs(mu->eta());
315 //_d0 = mu->dB(); // global track
316 _d0 = -(mu->innerTrack()->dxy(_bs));
317 _chi2ndof = mu->normChi2();
318 _nhits = mu->numberOfValidHits();
319 reco::TrackRef _track = mu->track();
320 //_____ check that there is a track
321 if(_track.isNull()) continue;
322 //_nhits = _track->numberOfValidHits();
323 //_d0 = -(_track->dxy(_bs));
324 //_chi2ndof = _track->normalizedChi2();
325 //_chi2ndof = (_track->chi2())/(_track->ndof());
326 //_____ check that isolation is valid
327 if ( !(mu->isIsolationValid()) ) continue;
328 _emveto = mu->isolationR03().emVetoEt;
329 _hadveto = mu->isolationR03().hadVetoEt;
330 _trackiso = mu->isolationR03().sumPt;
331 _caloiso = mu->isolationR03().emEt + mu->isolationR03().hadEt;
332 _reliso = (_trackiso + _caloiso)/_pt;
333 if (_pt > mPt &&
334 _abseta < mEta &&
335 _nhits >= mNHits &&
336 fabs(_d0) < mD0 &&
337 _chi2ndof < mChi2Ndof &&
338 fabs(_emveto) < mEmVeto &&
339 fabs(_hadveto) < mHadVeto &&
340 fabs(_reliso) < mRelIso
341 ){
342 passed_muons.count();
343 muon_index = (int)(mu - muons->begin());
344 if (passed_muons.getCount()==1){
345 cout << _pt << ", ";
346 cout << _abseta << ", ";
347 cout << _nhits << ", ";
348 cout << _d0 << ", ";
349 cout << _chi2ndof << ", ";
350 cout << _emveto << ", ";
351 cout << _hadveto << ", ";
352 cout << _trackiso << ", ";
353 cout << _caloiso << ", ";
354 cout << _reliso;
355 //cout << endl;
356 }
357 }
358 }
359 cout << "Number of good muons: " << passed_muons.getCount();
360 result = passed_muons.getCount()==1;
361 return result;
362 }
363
364
365 bool CutFlow::is_jets(std::string mLabel,
366 double mPt,
367 double mEta,
368 const edm::Event& iEvent){
369 bool result = false;
370
371 //cout << "Cut Flow: passed stage 2: ";
372
373 Handle< vector< pat::Jet > > jets;
374 iEvent . getByLabel( mLabel, jets );
375 MeanCounter c_jets("");
376 MeanCounter passed_jets("");
377 for ( vector<pat::Jet>::const_iterator jet = jets -> begin(); jet != jets -> end(); jet++ ){
378 c_jets.count();
379 double _pt = jet->pt();
380 double _eta = jet->eta();
381 if (
382 ( fabs( _eta ) < mEta ) &&
383 ( _pt > mPt )
384 )
385 {
386 passed_jets.count();
387 }
388 if (c_jets.getCount()<=4){
389 cout << _pt << ", ";
390 cout << _eta;
391 }
392 }
393 //cout << endl;
394 result = passed_jets.getCount() >=4;
395 return result;
396 }
397
398
399 bool CutFlow::no_loose_muon(std::string mLabel,
400 double mEta,
401 double mPt,
402 double mRelIso,
403 const edm::Event& iEvent){
404 bool result = false;
405 //cout << "Cut Flow: passed stage 3: ";
406 Handle< vector< pat::Muon > > muons;
407 iEvent . getByLabel( "selectedLayer1Muons", muons );
408 //
409 //_____ loop over muons _______________________________________________
410 //
411 double _pt = 0.0;
412 double _abseta = 0.0;
413 double _trackiso = 10000.0;
414 double _caloiso = 10000.0;
415 double _reliso = 10000.0;
416 MeanCounter passed_muons("");
417 for ( vector<pat::Muon>::const_iterator mu = muons -> begin(); mu != muons -> end(); mu++){
418 //_____ check that the muon is global
419 if( !(mu->isGlobalMuon()) ) continue;
420 _pt = mu->pt();
421 _abseta = fabs(mu->eta());
422 //_____ check that isolation is valid
423 if ( !(mu->isIsolationValid()) ) continue;
424 _trackiso = mu->isolationR03().sumPt;
425 _caloiso = mu->isolationR03().emEt + mu->isolationR03().hadEt;
426 _reliso = (_trackiso + _caloiso)/_pt;
427 int _index = (int)(mu - muons->begin());
428 if (_pt > mPt &&
429 _abseta < mEta &&
430 _reliso < mRelIso &&
431 _index != muon_index
432 ){
433 passed_muons.count();
434 }
435 }
436 result = passed_muons.getCount()==0;
437 //cout << endl;
438 return result;
439 }
440
441
442 bool CutFlow::no_loose_electron(std::string mLabel,
443 double mEta,
444 double mPt,
445 double mRelIso,
446 const edm::Event& iEvent){
447 bool result = false;
448 //cout << "Cut Flow: passed stage 4: ";
449 Handle< vector< pat::Electron > > electrons;
450 iEvent . getByLabel( "selectedLayer1Electrons", electrons );
451 //
452 //_____ loop over electrons _______________________________________________
453 //
454 double _pt = 0.0;
455 double _abseta = 0.0;
456 double _trackiso = 10000.0;
457 double _caloiso = 10000.0;
458 double _reliso = 10000.0;
459 MeanCounter passed_electrons("");
460 for ( vector<pat::Electron>::const_iterator e = electrons -> begin(); e != electrons -> end(); e++){
461 _pt = e->pt();
462 _abseta = fabs(e->eta());
463 _trackiso = e->dr03IsolationVariables().tkSumPt;
464 _caloiso = e->dr03IsolationVariables().ecalRecHitSumEt + e->dr03IsolationVariables().hcalDepth1TowerSumEt + e->dr03IsolationVariables().hcalDepth2TowerSumEt;
465 _reliso = (_trackiso + _caloiso)/_pt;
466 int _index = (int)(e - electrons->begin());
467 if (_pt > mPt &&
468 _abseta < mEta &&
469 _reliso < mRelIso &&
470 _index != electron_index
471 ){
472 passed_electrons.count();
473 }
474 }
475 result = passed_electrons.getCount()==0;
476 //cout << endl;
477 return result;
478 }
479
480
481 bool CutFlow::has_electron(std::string mLabel,
482 const edm::Event& iEvent){
483 bool result = false;
484 Handle< vector< pat::Electron > > electrons;
485 iEvent . getByLabel( mLabel, electrons );
486 result = electrons->size()>0;
487 return result;
488 }
489
490
491 bool CutFlow::is_electron(std::string mLabel,
492 double mPt,
493 double mEta,
494 double mD0,
495 double mRelIso,
496 const edm::Event& iEvent){
497 bool result = false;
498 //cout << "Cut Flow: passed stage 1: ";
499 Handle< vector< pat::Electron > > electrons;
500 iEvent . getByLabel( "selectedLayer1Electrons", electrons );
501 Handle< reco::BeamSpot > beamSpotHandle;
502 iEvent . getByLabel( "offlineBeamSpot", beamSpotHandle);
503 //
504 //_____ Beam spot _____________________________________________________
505 //
506 reco::BeamSpot beamSpot;
507 if ( beamSpotHandle.isValid() ){
508 beamSpot = *beamSpotHandle;
509 }
510 else{
511 edm::LogInfo("TtLJetsAnalyzer")
512 << "No beam spot available from EventSetup \n";
513 }
514 const reco::BeamSpot::Point _bs(beamSpot.x0(), beamSpot.y0(),beamSpot.z0());
515 //
516 //_____ loop over electrons _______________________________________________
517 //
518 double _pt = 0.0;
519 double _abseta = 0.0;
520 double _d0 = 0.0;
521 double _d0_pat = 0.0;
522 double _trackiso = 10000.0;
523 double _ecaloiso = 10000.0;
524 double _ecaloiso2 = 10000.0;
525 double _pat_ecaloiso = 10000.0;
526 double _hcaloiso1 = 10000.0;
527 double _hcaloiso2 = 10000.0;
528 double _caloiso = 10000.0;
529 double _reliso = 10000.0;
530 MeanCounter passed_electrons("");
531 for ( vector<pat::Electron>::const_iterator e = electrons -> begin(); e != electrons -> end(); e++){
532 //_____ check that the electron is robust-tight
533 if( !(e->electronID("eidRobustTight")) ) continue;
534 _pt = e->pt();
535 _abseta = fabs(e->eta());
536 _d0_pat = e->dB();
537 _d0 = -(e->gsfTrack()->dxy(_bs));
538 _trackiso = e->dr03IsolationVariables().tkSumPt;
539 _ecaloiso = e->dr04IsolationVariables().ecalRecHitSumEt;
540 _ecaloiso2 = e->dr03IsolationVariables().ecalRecHitSumEt;
541 //_pat_ecaloiso = e->isolation(pat::ECalIso);
542 _hcaloiso1 = e->dr04IsolationVariables().hcalDepth1TowerSumEt;
543 _hcaloiso2 = e->dr04IsolationVariables().hcalDepth2TowerSumEt;
544 _caloiso = _ecaloiso + _hcaloiso1 + _hcaloiso2;
545 //_caloiso = _pat_ecaloiso +_hcaloiso1 + _hcaloiso2;
546 _reliso = (_trackiso + _caloiso)/_pt;
547 if (_pt > mPt &&
548 _abseta < mEta &&
549 fabs(_d0) < mD0 &&
550 fabs(_reliso) < mRelIso
551 ){
552 passed_electrons.count();
553 electron_index = (int)(e - electrons->begin());
554 if (passed_electrons.getCount()==1){
555 cout << "pt=" << _pt << ", abseta=";
556 cout << _abseta << ", d0=";
557 cout << _d0 << ", trackiso=";
558 cout << _trackiso << ", ecaloiso04=";
559 cout << _ecaloiso << ", ecaloiso03=";
560 cout << _ecaloiso2 << ", pat_ecaloiso=";
561 cout << _pat_ecaloiso << ", hcaloiso1=";
562 cout << _hcaloiso1 << ", hcaloiso2=";
563 cout << _hcaloiso2 << ", caloiso04=";
564 cout << _caloiso << ", reliso=";
565 cout << _reliso;
566 //cout << endl;
567 }
568 }
569 }
570 //cout << "Number of good electrons: " << passed_electrons.getCount();
571 result = passed_electrons.getCount()==1;
572 return result;
573 }
574
575
576 int CutFlow::count_jets(std::string jLabel,
577 double jPt,
578 double jEta,
579 std::string lLabel,
580 double lPt,
581 double lEta,
582 double lD0,
583 double lDR,
584 const edm::Event & iEvent){
585 MeanCounter _njets("");
586 Handle< vector< pat::Jet > > jets;
587 iEvent . getByLabel( jLabel, jets );
588 Handle< vector< pat::Electron > > electrons;
589 iEvent . getByLabel( lLabel, electrons );
590 Handle< reco::BeamSpot > beamSpotHandle;
591 iEvent . getByLabel( "offlineBeamSpot", beamSpotHandle);
592 //
593 //_____ Beam spot _____________________________________________________
594 //
595 reco::BeamSpot beamSpot;
596 if ( beamSpotHandle.isValid() ){
597 beamSpot = *beamSpotHandle;
598 }
599 else{
600 edm::LogInfo("TtLJetsAnalyzer")
601 << "No beam spot available from EventSetup \n";
602 }
603 const reco::BeamSpot::Point _bs(beamSpot.x0(), beamSpot.y0(),beamSpot.z0());
604 //
605 //_____ loop over jets _______________________________________________
606 //
607 double j_pt = 0.0;
608 double j_eta = 0.0;
609 for ( vector<pat::Jet>::const_iterator jet = jets -> begin(); jet != jets -> end(); ++jet){
610 j_pt = jet->pt();
611 j_eta = jet->eta();
612 //
613 //_____ loop over electrons _______________________________________________
614 //
615 double e_pt = 0.0;
616 double e_abseta = 0.0;
617 double e_d0 = 0.0;
618 double e_dr = 0.0;
619 bool jet_removed = false;
620 for ( vector<pat::Electron>::const_iterator e = electrons -> begin(); e != electrons -> end(); e++){
621 //_____ check that the electron is robust-tight
622 if( !(e->electronID("eidRobustTight")) ) continue;
623 e_pt = e->pt();
624 e_abseta = fabs(e->eta());
625 e_d0 = -(e->gsfTrack()->dxy(_bs));
626 e_dr = reco::deltaR(*jet, *e);
627 if ( e_pt > lPt && e_abseta < lEta && fabs(e_d0)<lD0 && e_dr<lDR){
628 jet_removed = true;
629 }
630 }
631 if (j_pt>jPt && fabs(j_eta)<jEta && !jet_removed) _njets.count();
632 }
633 cout << "Number of jets after cleaning: " << _njets.getCount();
634 return _njets.getCount();
635 }
636
637 //define this as a plug-in
638 DEFINE_FWK_MODULE(CutFlow);