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/09 19:56:45 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"
#include "DataFormats/EgammaCandidates/interface/GsfElectron.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);

  muon_index = -1;
  electron_index = -1;
}


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) ){
    }
    else break;
    if ( has_global_muon("selectedLayer1Muons", 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;
    if ( is_jets("selectedLayer1Jets",
                 30.0,
                 2.4,
                 iEvent) ){
      c_3->count();
    }
    else break;
    if ( no_loose_muon("selectedLayer1Muons",
                       2.5,
                       10.0,
                       0.2,
                       iEvent) ){
      c_4->count();
    }
    else break;
    if ( no_loose_electron("selectedLayer1Electrons",
                           2.5,
                           15.0,
                           0.2,
                           iEvent) ){
      c_5->count();
    }
    else break;
    cout << "Cut Flow: passed stage 5 " << 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 mTriggerName, 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(mTriggerName)!=std::string::npos ){
      result = trig->wasAccept();
    }
  }
  return result;
}

bool CutFlow::has_global_muon(std::string mLabel,
                              const edm::Event& iEvent){
  bool result = false;
  cout << "Cut Flow: passed trigger ";
  Handle< vector< pat::Muon > >     muons;
  iEvent . getByLabel( mLabel, muons );
  //
  //_____ loop over muons _______________________________________________
  //
  MeanCounter passed_muons("");
  for ( vector<pat::Muon>::const_iterator mu = muons -> begin(); mu != muons -> end(); mu++){
    if ( mu->isGlobalMuon() ){
      passed_muons.count();
    }
  }
  result = passed_muons.getCount()>0;
  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;
  cout << "Cut Flow: passed stage 1: ";
  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;
  MeanCounter passed_muons("");
  for ( vector<pat::Muon>::const_iterator mu = muons -> begin(); mu != muons -> end(); mu++){
    //_____ check that the muon is global
    if( !(mu->isGlobalMuon()) ) continue;
    _pt = mu->pt();
    _abseta = fabs(mu->eta());
    //_d0 = mu->dB(); // global track
    _d0 = -(mu->innerTrack()->dxy(_bs));
    _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 &&
        fabs(_d0) < mD0 &&
        _chi2ndof < mChi2Ndof &&
        fabs(_emveto) < mEmVeto &&
        fabs(_hadveto) < mHadVeto &&
        fabs(_reliso) < mRelIso
        ){
      passed_muons.count();
      muon_index = (int)(mu - muons->begin());
      if (passed_muons.getCount()==1){
        cout << _pt << ", ";
        cout << _abseta << ", ";
        cout << _nhits << ", ";
        cout << _d0 << ", ";
        cout << _chi2ndof << ", ";
        cout << _emveto << ", ";
        cout << _hadveto << ", ";
        cout << _trackiso << ", ";
        cout << _caloiso << ", ";
        cout << _reliso << ", ";
        cout << endl;
      }
    }
  }
  result = passed_muons.getCount()==1;
  return result;
}
  

bool CutFlow::is_jets(std::string mLabel,
                      double mPt,
                      double mEta,
                      const edm::Event& iEvent){
  bool result = false;

  cout << "Cut Flow: passed stage 2: ";

  Handle< vector< pat::Jet > >      jets;
  iEvent . getByLabel( mLabel, jets );
  MeanCounter c_jets("");
  MeanCounter passed_jets("");
  for ( vector<pat::Jet>::const_iterator jet = jets -> begin(); jet != jets -> end(); jet++ ){
    c_jets.count();
    double _pt = jet->pt();
    double _eta = jet->eta();
    if (
        ( fabs( _eta ) < mEta ) &&
        ( _pt > mPt )
        )
      {
        passed_jets.count();
      }
    if (c_jets.getCount()<=4){
      cout << _pt << ", ";
      cout << _eta << ", ";
    }
  }
  cout << endl;
  result = passed_jets.getCount() >=4;
  return result;
}


bool CutFlow::no_loose_muon(std::string mLabel,
                            double mEta,
                            double mPt,
                            double mRelIso,
                            const edm::Event& iEvent){
  bool result = false;
  cout << "Cut Flow: passed stage 3: ";
  Handle< vector< pat::Muon > >     muons;
  iEvent . getByLabel( "selectedLayer1Muons", muons );
  //
  //_____ loop over muons _______________________________________________
  //
  double _pt = 0.0;
  double _abseta = 0.0;
  double _trackiso = 10000.0;
  double _caloiso = 10000.0;
  double _reliso = 10000.0;
  MeanCounter passed_muons("");
  for ( vector<pat::Muon>::const_iterator mu = muons -> begin(); mu != muons -> end(); mu++){
    //_____ check that the muon is global
    if( !(mu->isGlobalMuon()) ) continue;
    _pt = mu->pt();
    _abseta = fabs(mu->eta());
    //_____ check that isolation is valid
    if ( !(mu->isIsolationValid()) ) continue;
    _trackiso = mu->isolationR03().sumPt;
    _caloiso = mu->isolationR03().emEt + mu->isolationR03().hadEt;
    _reliso = (_trackiso + _caloiso)/_pt;
    int _index = (int)(mu - muons->begin());
    if (_pt > mPt &&
        _abseta < mEta &&
        _reliso < mRelIso &&
        _index != muon_index
        ){
      passed_muons.count();
      cout << "DEBUG: found a loose muon!!!!!" << endl;
      cout << "DEBUG: index = " << _index << endl;
      cout << "DEBUG: abseta = " << _abseta << endl;
      cout << "DEBUG: pt = " << _pt << endl;
      cout << "DEBUG: reliso = " << _reliso << endl << endl;
    }
  }
  result = passed_muons.getCount()==0;
  return result;
}


bool CutFlow::no_loose_electron(std::string mLabel,
                                double mEta,
                                double mPt,
                                double mRelIso,
                                const edm::Event& iEvent){
  bool result = false;
  cout << "Cut Flow: passed stage 4: ";
  Handle< vector< pat::Electron > >     electrons;
  iEvent . getByLabel( "selectedLayer1Electrons", electrons );
  //
  //_____ loop over electrons _______________________________________________
  //
  double _pt = 0.0;
  double _abseta = 0.0;
  double _trackiso = 10000.0;
  double _caloiso = 10000.0;
  double _reliso = 10000.0;
  MeanCounter passed_electrons("");
  for ( vector<pat::Electron>::const_iterator e = electrons -> begin(); e != electrons -> end(); e++){
    _pt = e->pt();
    _abseta = fabs(e->eta());
    _trackiso = e->dr03IsolationVariables().tkSumPt;
    _caloiso = e->dr03IsolationVariables().ecalRecHitSumEt + e->dr03IsolationVariables().hcalDepth1TowerSumEt + e->dr03IsolationVariables().hcalDepth2TowerSumEt;
    _reliso = (_trackiso + _caloiso)/_pt;
    int _index = (int)(e - electrons->begin());
    //cout << "DEBUG: " << _index << endl;
    //cout << "DEBUG: " << _abseta << endl;
    //cout << "DEBUG: " << _pt << endl;
    //cout << "DEBUG: " << _reliso << endl << endl;
    if (_pt > mPt &&
        _abseta < mEta &&
        _reliso < mRelIso &&
        _index != electron_index
        ){
      passed_electrons.count();
      cout << "DEBUG: found a loose electron!!!!!" << endl;
      cout << "DEBUG: index = " << _index << endl;
      cout << "DEBUG: abseta = " << _abseta << endl;
      cout << "DEBUG: pt = " << _pt << endl;
      cout << "DEBUG: reliso = " << _reliso << endl << endl;
    }
  }
  result = passed_electrons.getCount()==0;
  return result;
}


//define this as a plug-in
DEFINE_FWK_MODULE(CutFlow);
Revision:	1.2
Committed:	Sat Oct 10 05:26:18 2009 UTC (15 years, 6 months ago) by kukartse
Content type:	text/plain
Branch:	MAIN
Changes since 1.1:	+211 -24 lines
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	kukartse	1.2	// $Id: CutFlow.cc,v 1.1 2009/10/09 19:56:45 kukartse Exp $
17	kukartse	1.1	//
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	kukartse	1.2	#include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
30	kukartse	1.1
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	kukartse	1.2
46			muon_index = -1;
47			electron_index = -1;
48	kukartse	1.1	}
49
50
51			CutFlow::~CutFlow()
52			{
53			delete c_total;
54			delete c_1;
55			delete c_2;
56			delete c_3;
57			delete c_4;
58			delete c_5;
59			}
60
61
62			void
63			CutFlow::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
64			{
65			c_total->count();
66			while(1){
67	kukartse	1.2	if ( is_trigger("HLT_Mu9", iEvent) ){
68			}
69			else break;
70			if ( has_global_muon("selectedLayer1Muons", iEvent) ){
71			c_1->count();
72			}
73	kukartse	1.1	else break;
74			if ( is_muon("selectedLayer1Muons",
75			20.0,
76			2.1,
77			11,
78			0.02,
79			10,
80			4.0,
81			6.0,
82			0.05,
83			iEvent) ){
84			c_2->count();
85			}
86			else break;
87	kukartse	1.2	if ( is_jets("selectedLayer1Jets",
88			30.0,
89			2.4,
90			iEvent) ){
91			c_3->count();
92			}
93			else break;
94			if ( no_loose_muon("selectedLayer1Muons",
95			2.5,
96			10.0,
97			0.2,
98			iEvent) ){
99			c_4->count();
100			}
101			else break;
102			if ( no_loose_electron("selectedLayer1Electrons",
103			2.5,
104			15.0,
105			0.2,
106			iEvent) ){
107			c_5->count();
108			}
109			else break;
110			cout << "Cut Flow: passed stage 5 " << endl;
111	kukartse	1.1	break;
112			}
113			}
114
115
116			void
117			CutFlow::beginJob()
118			{
119			}
120
121
122			void
123			CutFlow::endJob() {
124			cout << "========= Cut Flow Report ==========================================" << endl;
125			cout << "Total events: " << c_total->getCount() << endl;
126			cout << "Pass stage 1 (trigger+global): " << c_1->getCount() << endl;
127			cout << "Pass stage 2 (muon kinematics): " << c_2->getCount() << endl;
128			cout << "Pass stage 3 (jets): " << c_3->getCount() << endl;
129			cout << "Pass stage 4 (no loose muon): " << c_4->getCount() << endl;
130			cout << "Pass stage 5 (no loose electron): " << c_5->getCount() << endl;
131			cout << "================================================================" << endl;
132			}
133
134	kukartse	1.2	bool CutFlow::is_trigger(std::string mTriggerName, const edm::Event& iEvent){
135	kukartse	1.1	bool result = false;
136			Handle<std::vector<pat::TriggerPath> > trigs;
137			iEvent.getByLabel("patTrigger", trigs);
138			for (std::vector<pat::TriggerPath>::const_iterator trig = trigs->begin();
139			trig!=trigs->end();
140			++trig){
141			//cout << "Trigger name and number: " << trig->name() << ", " << trig->index() << endl;
142	kukartse	1.2	if( trig->name().find(mTriggerName)!=std::string::npos ){
143	kukartse	1.1	result = trig->wasAccept();
144			}
145			}
146			return result;
147			}
148
149	kukartse	1.2	bool CutFlow::has_global_muon(std::string mLabel,
150			const edm::Event& iEvent){
151			bool result = false;
152			cout << "Cut Flow: passed trigger ";
153			Handle< vector< pat::Muon > > muons;
154			iEvent . getByLabel( mLabel, muons );
155			//
156			//_____ loop over muons _______________________________________________
157			//
158			MeanCounter passed_muons("");
159			for ( vector<pat::Muon>::const_iterator mu = muons -> begin(); mu != muons -> end(); mu++){
160			if ( mu->isGlobalMuon() ){
161			passed_muons.count();
162			}
163			}
164			result = passed_muons.getCount()>0;
165			return result;
166			}
167
168
169	kukartse	1.1	bool CutFlow::is_muon(std::string mLabel,
170			double mPt,
171			double mEta,
172			int mNHits,
173			double mD0,
174			double mChi2Ndof,
175			double mEmVeto,
176			double mHadVeto,
177			double mRelIso,
178			const edm::Event& iEvent){
179			bool result = false;
180	kukartse	1.2	cout << "Cut Flow: passed stage 1: ";
181	kukartse	1.1	Handle< vector< pat::Muon > > muons;
182			iEvent . getByLabel( "selectedLayer1Muons", muons );
183			Handle< reco::BeamSpot > beamSpotHandle;
184			iEvent . getByLabel( "offlineBeamSpot", beamSpotHandle);
185			//
186			//_____ Beam spot _____________________________________________________
187			//
188			reco::BeamSpot beamSpot;
189			if ( beamSpotHandle.isValid() ){
190			beamSpot = *beamSpotHandle;
191			}
192			else{
193			edm::LogInfo("TtLJetsAnalyzer")
194			<< "No beam spot available from EventSetup \n";
195			}
196			const reco::BeamSpot::Point _bs(beamSpot.x0(), beamSpot.y0(),beamSpot.z0());
197			//
198			//_____ loop over muons _______________________________________________
199			//
200			double _pt = 0.0;
201			double _abseta = 0.0;
202			int _nhits = 0;
203			double _d0 = 0.0;
204			double _chi2ndof = 0.0;
205			double _emveto = 10000.0;
206			double _hadveto = 10000.0;
207			double _trackiso = 10000.0;
208			double _caloiso = 10000.0;
209			double _reliso = 10000.0;
210	kukartse	1.2	MeanCounter passed_muons("");
211	kukartse	1.1	for ( vector<pat::Muon>::const_iterator mu = muons -> begin(); mu != muons -> end(); mu++){
212	kukartse	1.2	//_____ check that the muon is global
213			if( !(mu->isGlobalMuon()) ) continue;
214	kukartse	1.1	_pt = mu->pt();
215			_abseta = fabs(mu->eta());
216	kukartse	1.2	//_d0 = mu->dB(); // global track
217			_d0 = -(mu->innerTrack()->dxy(_bs));
218	kukartse	1.1	_chi2ndof = mu->normChi2();
219			_nhits = mu->numberOfValidHits();
220			reco::TrackRef _track = mu->track();
221			//_____ check that there is a track
222			if(_track.isNull()) continue;
223			//_nhits = _track->numberOfValidHits();
224			//_d0 = -(_track->dxy(_bs));
225			//_chi2ndof = _track->normalizedChi2();
226			//_chi2ndof = (_track->chi2())/(_track->ndof());
227			//_____ check that isolation is valid
228			if ( !(mu->isIsolationValid()) ) continue;
229			_emveto = mu->isolationR03().emVetoEt;
230			_hadveto = mu->isolationR03().hadVetoEt;
231			_trackiso = mu->isolationR03().sumPt;
232			_caloiso = mu->isolationR03().emEt + mu->isolationR03().hadEt;
233			_reliso = (_trackiso + _caloiso)/_pt;
234			if (_pt > mPt &&
235			_abseta < mEta &&
236			_nhits >= mNHits &&
237	kukartse	1.2	fabs(_d0) < mD0 &&
238	kukartse	1.1	_chi2ndof < mChi2Ndof &&
239	kukartse	1.2	fabs(_emveto) < mEmVeto &&
240			fabs(_hadveto) < mHadVeto &&
241			fabs(_reliso) < mRelIso
242	kukartse	1.1	){
243	kukartse	1.2	passed_muons.count();
244			muon_index = (int)(mu - muons->begin());
245			if (passed_muons.getCount()==1){
246			cout << _pt << ", ";
247			cout << _abseta << ", ";
248			cout << _nhits << ", ";
249			cout << _d0 << ", ";
250			cout << _chi2ndof << ", ";
251			cout << _emveto << ", ";
252			cout << _hadveto << ", ";
253			cout << _trackiso << ", ";
254			cout << _caloiso << ", ";
255			cout << _reliso << ", ";
256			cout << endl;
257			}
258			}
259			}
260			result = passed_muons.getCount()==1;
261			return result;
262			}
263
264
265			bool CutFlow::is_jets(std::string mLabel,
266			double mPt,
267			double mEta,
268			const edm::Event& iEvent){
269			bool result = false;
270
271			cout << "Cut Flow: passed stage 2: ";
272
273			Handle< vector< pat::Jet > > jets;
274			iEvent . getByLabel( mLabel, jets );
275			MeanCounter c_jets("");
276			MeanCounter passed_jets("");
277			for ( vector<pat::Jet>::const_iterator jet = jets -> begin(); jet != jets -> end(); jet++ ){
278			c_jets.count();
279			double _pt = jet->pt();
280			double _eta = jet->eta();
281			if (
282			( fabs( _eta ) < mEta ) &&
283			( _pt > mPt )
284			)
285			{
286			passed_jets.count();
287			}
288			if (c_jets.getCount()<=4){
289			cout << _pt << ", ";
290			cout << _eta << ", ";
291	kukartse	1.1	}
292			}
293			cout << endl;
294	kukartse	1.2	result = passed_jets.getCount() >=4;
295	kukartse	1.1	return result;
296			}
297	kukartse	1.2
298
299			bool CutFlow::no_loose_muon(std::string mLabel,
300			double mEta,
301			double mPt,
302			double mRelIso,
303			const edm::Event& iEvent){
304			bool result = false;
305			cout << "Cut Flow: passed stage 3: ";
306			Handle< vector< pat::Muon > > muons;
307			iEvent . getByLabel( "selectedLayer1Muons", muons );
308			//
309			//_____ loop over muons _______________________________________________
310			//
311			double _pt = 0.0;
312			double _abseta = 0.0;
313			double _trackiso = 10000.0;
314			double _caloiso = 10000.0;
315			double _reliso = 10000.0;
316			MeanCounter passed_muons("");
317			for ( vector<pat::Muon>::const_iterator mu = muons -> begin(); mu != muons -> end(); mu++){
318			//_____ check that the muon is global
319			if( !(mu->isGlobalMuon()) ) continue;
320			_pt = mu->pt();
321			_abseta = fabs(mu->eta());
322			//_____ check that isolation is valid
323			if ( !(mu->isIsolationValid()) ) continue;
324			_trackiso = mu->isolationR03().sumPt;
325			_caloiso = mu->isolationR03().emEt + mu->isolationR03().hadEt;
326			_reliso = (_trackiso + _caloiso)/_pt;
327			int _index = (int)(mu - muons->begin());
328			if (_pt > mPt &&
329			_abseta < mEta &&
330			_reliso < mRelIso &&
331			_index != muon_index
332			){
333			passed_muons.count();
334			cout << "DEBUG: found a loose muon!!!!!" << endl;
335			cout << "DEBUG: index = " << _index << endl;
336			cout << "DEBUG: abseta = " << _abseta << endl;
337			cout << "DEBUG: pt = " << _pt << endl;
338			cout << "DEBUG: reliso = " << _reliso << endl << endl;
339			}
340			}
341			result = passed_muons.getCount()==0;
342			return result;
343			}
344
345
346			bool CutFlow::no_loose_electron(std::string mLabel,
347			double mEta,
348			double mPt,
349			double mRelIso,
350			const edm::Event& iEvent){
351			bool result = false;
352			cout << "Cut Flow: passed stage 4: ";
353			Handle< vector< pat::Electron > > electrons;
354			iEvent . getByLabel( "selectedLayer1Electrons", electrons );
355			//
356			//_____ loop over electrons _______________________________________________
357			//
358			double _pt = 0.0;
359			double _abseta = 0.0;
360			double _trackiso = 10000.0;
361			double _caloiso = 10000.0;
362			double _reliso = 10000.0;
363			MeanCounter passed_electrons("");
364			for ( vector<pat::Electron>::const_iterator e = electrons -> begin(); e != electrons -> end(); e++){
365			_pt = e->pt();
366			_abseta = fabs(e->eta());
367			_trackiso = e->dr03IsolationVariables().tkSumPt;
368			_caloiso = e->dr03IsolationVariables().ecalRecHitSumEt + e->dr03IsolationVariables().hcalDepth1TowerSumEt + e->dr03IsolationVariables().hcalDepth2TowerSumEt;
369			_reliso = (_trackiso + _caloiso)/_pt;
370			int _index = (int)(e - electrons->begin());
371			//cout << "DEBUG: " << _index << endl;
372			//cout << "DEBUG: " << _abseta << endl;
373			//cout << "DEBUG: " << _pt << endl;
374			//cout << "DEBUG: " << _reliso << endl << endl;
375			if (_pt > mPt &&
376			_abseta < mEta &&
377			_reliso < mRelIso &&
378			_index != electron_index
379			){
380			passed_electrons.count();
381			cout << "DEBUG: found a loose electron!!!!!" << endl;
382			cout << "DEBUG: index = " << _index << endl;
383			cout << "DEBUG: abseta = " << _abseta << endl;
384			cout << "DEBUG: pt = " << _pt << endl;
385			cout << "DEBUG: reliso = " << _reliso << endl << endl;
386			}
387			}
388			result = passed_electrons.getCount()==0;
389			return result;
390			}
391
392	kukartse	1.1
393			//define this as a plug-in
394			DEFINE_FWK_MODULE(CutFlow);