Utils/src/CutFlow.cc

// -*- C++ -*-
//
// Package:    CutFlow
// Class:      CutFlow
// 
/**\class CutFlow CutFlow.h LJMet/Utils/interface/CutFlow.h

 Description: times various aspects of data processing

 Implementation:
     Keep this class as simple and fast as possible
*/
//
// Original Author:  "Gennadiy Kukartsev"
//         Created:  Tue Oct  6 13:39:12 CDT 2009
// $Id: CutFlow.cc,v 1.1 2009/10/08 22:51:48 kukartse Exp $
//
//

#include "LJMet/Utils/interface/CutFlow.h"
#include "FWCore/Framework/interface/TriggerNames.h"
#include "DataFormats/Common/interface/TriggerResults.h"
#include "DataFormats/PatCandidates/interface/Jet.h"
#include "DataFormats/PatCandidates/interface/Electron.h"
#include "DataFormats/PatCandidates/interface/Muon.h"
#include "DataFormats/PatCandidates/interface/MET.h"
#include "DataFormats/PatCandidates/interface/TriggerPath.h"
#include "DataFormats/BeamSpot/interface/BeamSpot.h"


using namespace edm;

CutFlow::CutFlow(const edm::ParameterSet& iConfig){
  c_total = new MeanCounter("Processing event number: ", 0, 1);
  c_1 = new MeanCounter("CutFlow stage 1: ", 0, 1);
  c_2 = new MeanCounter("CutFlow stage 2: ", 0, 1);
  c_3 = new MeanCounter("CutFlow stage 3: ", 0, 1);
  c_4 = new MeanCounter("CutFlow stage 4: ", 0, 1);
  c_5 = new MeanCounter("CutFlow stage 5: ", 0, 1);
  c_1->setPrintCount(false);
  c_2->setPrintCount(false);
  c_3->setPrintCount(false);
  c_4->setPrintCount(false);
  c_5->setPrintCount(false);
}


CutFlow::~CutFlow()
{
  delete c_total;
  delete c_1;
  delete c_2;
  delete c_3;
  delete c_4;
  delete c_5;
}


void
CutFlow::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
  c_total->count();
  while(1){
    if ( is_trigger("HLT_Mu9", iEvent) ) c_1->count();
    else break;
    if ( is_muon("selectedLayer1Muons",
                 20.0,
                 2.1,
                 11,
                 0.02,
                 10,
                 4.0,
                 6.0,
                 0.05,
                 iEvent) ){
      c_2->count();
    }
    else break;
    cout << "Cut Flow: passed muon" << endl;
    break;
  }
}


void 
CutFlow::beginJob()
{
}


void 
CutFlow::endJob() {
  cout << "========= Cut Flow Report ==========================================" << endl;
  cout << "Total events: " << c_total->getCount() << endl;
  cout << "Pass stage 1 (trigger+global): " << c_1->getCount() << endl;
  cout << "Pass stage 2 (muon kinematics): " << c_2->getCount() << endl;
  cout << "Pass stage 3 (jets): " << c_3->getCount() << endl;
  cout << "Pass stage 4 (no loose muon): " << c_4->getCount() << endl;
  cout << "Pass stage 5 (no loose electron): " << c_5->getCount() << endl;
  cout << "================================================================" << endl;
}

bool CutFlow::is_trigger(std::string mLabel, const edm::Event& iEvent){
  bool result = false;
  Handle<std::vector<pat::TriggerPath> > trigs;
  iEvent.getByLabel("patTrigger", trigs);
  for (std::vector<pat::TriggerPath>::const_iterator trig = trigs->begin();
       trig!=trigs->end();
       ++trig){
    //cout << "Trigger name and number: " << trig->name() << ",  " << trig->index() << endl;
    if( trig->name().find(mLabel)!=std::string::npos ){
      result = trig->wasAccept();
    }
  }
  return result;
}

bool CutFlow::is_muon(std::string mLabel,
                      double mPt,
                      double mEta, 
                      int mNHits,
                      double mD0,
                      double mChi2Ndof,
                      double mEmVeto,
                      double mHadVeto,
                      double mRelIso,
                      const edm::Event& iEvent){
  bool result = false;
  Handle< vector< pat::Muon > >     muons;
  iEvent . getByLabel( "selectedLayer1Muons", muons );
  Handle< reco::BeamSpot >          beamSpotHandle;
  iEvent . getByLabel( "offlineBeamSpot", beamSpotHandle);
  //
  //_____ Beam spot _____________________________________________________
  //
  reco::BeamSpot beamSpot;
  if ( beamSpotHandle.isValid() ){
    beamSpot = *beamSpotHandle;
  }
  else{
    edm::LogInfo("TtLJetsAnalyzer")
      << "No beam spot available from EventSetup \n";
  }
  const reco::BeamSpot::Point _bs(beamSpot.x0(), beamSpot.y0(),beamSpot.z0());
  //
  //_____ loop over muons _______________________________________________
  //
  double _pt = 0.0;
  double _abseta = 0.0;
  int _nhits = 0;
  double _d0 = 0.0;
  double _chi2ndof = 0.0;
  double _emveto = 10000.0;
  double _hadveto = 10000.0;
  double _trackiso = 10000.0;
  double _caloiso = 10000.0;
  double _reliso = 10000.0;
  for ( vector<pat::Muon>::const_iterator mu = muons -> begin(); mu != muons -> end(); mu++){
    _pt = mu->pt();
    _abseta = fabs(mu->eta());
    _d0 = mu->dB();
    _chi2ndof = mu->normChi2();
    _nhits = mu->numberOfValidHits();
    reco::TrackRef _track = mu->track();
    //_____ check that there is a track
    if(_track.isNull()) continue;
    //_nhits = _track->numberOfValidHits();
    //_d0 = -(_track->dxy(_bs));
    //_chi2ndof = _track->normalizedChi2();
    //_chi2ndof = (_track->chi2())/(_track->ndof());
    //_____ check that isolation is valid
    if ( !(mu->isIsolationValid()) ) continue;
    _emveto = mu->isolationR03().emVetoEt;
    _hadveto = mu->isolationR03().hadVetoEt;
    _trackiso = mu->isolationR03().sumPt;
    _caloiso = mu->isolationR03().emEt + mu->isolationR03().hadEt;
    _reliso = (_trackiso + _caloiso)/_pt;
    if (_pt > mPt &&
        _abseta < mEta &&
        _nhits >= mNHits &&
        _d0 < mD0 &&
        _chi2ndof < mChi2Ndof &&
        _emveto < mEmVeto &&
        _hadveto < mHadVeto &&
        _reliso < mRelIso
        ){
      result = true;
      break;
    }
  }
  cout << "Cut Flow: passed trigger: " << _pt << ", ";
  cout << _abseta << ", ";
  cout << _nhits << ", ";
  cout << _d0 << ", ";
  cout << _chi2ndof << ", ";
  cout << _emveto << ", ";
  cout << _hadveto << ", ";
  cout << _trackiso << ", ";
  cout << _caloiso << ", ";
  cout << _reliso << ", ";
  cout << endl;
  return result;
}
  

//define this as a plug-in
DEFINE_FWK_MODULE(CutFlow);
Revision:	1.1
Committed:	Fri Oct 9 19:56:45 2009 UTC (15 years, 6 months ago) by kukartse
Content type:	text/plain
Branch:	MAIN
Log Message:	* empty log message *
#	User	Rev	Content
1	kukartse	1.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.1 2009/10/08 22:51:48 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
30
31			using namespace edm;
32
33			CutFlow::CutFlow(const edm::ParameterSet& iConfig){
34			c_total = new MeanCounter("Processing event number: ", 0, 1);
35			c_1 = new MeanCounter("CutFlow stage 1: ", 0, 1);
36			c_2 = new MeanCounter("CutFlow stage 2: ", 0, 1);
37			c_3 = new MeanCounter("CutFlow stage 3: ", 0, 1);
38			c_4 = new MeanCounter("CutFlow stage 4: ", 0, 1);
39			c_5 = new MeanCounter("CutFlow stage 5: ", 0, 1);
40			c_1->setPrintCount(false);
41			c_2->setPrintCount(false);
42			c_3->setPrintCount(false);
43			c_4->setPrintCount(false);
44			c_5->setPrintCount(false);
45			}
46
47
48			CutFlow::~CutFlow()
49			{
50			delete c_total;
51			delete c_1;
52			delete c_2;
53			delete c_3;
54			delete c_4;
55			delete c_5;
56			}
57
58
59			void
60			CutFlow::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
61			{
62			c_total->count();
63			while(1){
64			if ( is_trigger("HLT_Mu9", iEvent) ) c_1->count();
65			else break;
66			if ( is_muon("selectedLayer1Muons",
67			20.0,
68			2.1,
69			11,
70			0.02,
71			10,
72			4.0,
73			6.0,
74			0.05,
75			iEvent) ){
76			c_2->count();
77			}
78			else break;
79			cout << "Cut Flow: passed muon" << endl;
80			break;
81			}
82			}
83
84
85			void
86			CutFlow::beginJob()
87			{
88			}
89
90
91			void
92			CutFlow::endJob() {
93			cout << "========= Cut Flow Report ==========================================" << endl;
94			cout << "Total events: " << c_total->getCount() << endl;
95			cout << "Pass stage 1 (trigger+global): " << c_1->getCount() << endl;
96			cout << "Pass stage 2 (muon kinematics): " << c_2->getCount() << endl;
97			cout << "Pass stage 3 (jets): " << c_3->getCount() << endl;
98			cout << "Pass stage 4 (no loose muon): " << c_4->getCount() << endl;
99			cout << "Pass stage 5 (no loose electron): " << c_5->getCount() << endl;
100			cout << "================================================================" << endl;
101			}
102
103			bool CutFlow::is_trigger(std::string mLabel, const edm::Event& iEvent){
104			bool result = false;
105			Handle<std::vector<pat::TriggerPath> > trigs;
106			iEvent.getByLabel("patTrigger", trigs);
107			for (std::vector<pat::TriggerPath>::const_iterator trig = trigs->begin();
108			trig!=trigs->end();
109			++trig){
110			//cout << "Trigger name and number: " << trig->name() << ", " << trig->index() << endl;
111			if( trig->name().find(mLabel)!=std::string::npos ){
112			result = trig->wasAccept();
113			}
114			}
115			return result;
116			}
117
118			bool CutFlow::is_muon(std::string mLabel,
119			double mPt,
120			double mEta,
121			int mNHits,
122			double mD0,
123			double mChi2Ndof,
124			double mEmVeto,
125			double mHadVeto,
126			double mRelIso,
127			const edm::Event& iEvent){
128			bool result = false;
129			Handle< vector< pat::Muon > > muons;
130			iEvent . getByLabel( "selectedLayer1Muons", muons );
131			Handle< reco::BeamSpot > beamSpotHandle;
132			iEvent . getByLabel( "offlineBeamSpot", beamSpotHandle);
133			//
134			//_____ Beam spot _____________________________________________________
135			//
136			reco::BeamSpot beamSpot;
137			if ( beamSpotHandle.isValid() ){
138			beamSpot = *beamSpotHandle;
139			}
140			else{
141			edm::LogInfo("TtLJetsAnalyzer")
142			<< "No beam spot available from EventSetup \n";
143			}
144			const reco::BeamSpot::Point _bs(beamSpot.x0(), beamSpot.y0(),beamSpot.z0());
145			//
146			//_____ loop over muons _______________________________________________
147			//
148			double _pt = 0.0;
149			double _abseta = 0.0;
150			int _nhits = 0;
151			double _d0 = 0.0;
152			double _chi2ndof = 0.0;
153			double _emveto = 10000.0;
154			double _hadveto = 10000.0;
155			double _trackiso = 10000.0;
156			double _caloiso = 10000.0;
157			double _reliso = 10000.0;
158			for ( vector<pat::Muon>::const_iterator mu = muons -> begin(); mu != muons -> end(); mu++){
159			_pt = mu->pt();
160			_abseta = fabs(mu->eta());
161			_d0 = mu->dB();
162			_chi2ndof = mu->normChi2();
163			_nhits = mu->numberOfValidHits();
164			reco::TrackRef _track = mu->track();
165			//_____ check that there is a track
166			if(_track.isNull()) continue;
167			//_nhits = _track->numberOfValidHits();
168			//_d0 = -(_track->dxy(_bs));
169			//_chi2ndof = _track->normalizedChi2();
170			//_chi2ndof = (_track->chi2())/(_track->ndof());
171			//_____ check that isolation is valid
172			if ( !(mu->isIsolationValid()) ) continue;
173			_emveto = mu->isolationR03().emVetoEt;
174			_hadveto = mu->isolationR03().hadVetoEt;
175			_trackiso = mu->isolationR03().sumPt;
176			_caloiso = mu->isolationR03().emEt + mu->isolationR03().hadEt;
177			_reliso = (_trackiso + _caloiso)/_pt;
178			if (_pt > mPt &&
179			_abseta < mEta &&
180			_nhits >= mNHits &&
181			_d0 < mD0 &&
182			_chi2ndof < mChi2Ndof &&
183			_emveto < mEmVeto &&
184			_hadveto < mHadVeto &&
185			_reliso < mRelIso
186			){
187			result = true;
188			break;
189			}
190			}
191			cout << "Cut Flow: passed trigger: " << _pt << ", ";
192			cout << _abseta << ", ";
193			cout << _nhits << ", ";
194			cout << _d0 << ", ";
195			cout << _chi2ndof << ", ";
196			cout << _emveto << ", ";
197			cout << _hadveto << ", ";
198			cout << _trackiso << ", ";
199			cout << _caloiso << ", ";
200			cout << _reliso << ", ";
201			cout << endl;
202			return result;
203			}
204
205
206			//define this as a plug-in
207			DEFINE_FWK_MODULE(CutFlow);