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, 7 months ago) by kukartse
Content type:	text/plain
Branch:	MAIN
Log Message:	* empty log message *
#	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.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);