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.2 2009/10/10 05:26:18 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();
  //
  //_____ mu+jets cut flow
  //
  /*
  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;
  }
  */
  //
  //_____ e+jets cut flow
  //
  while(1){
    if ( is_trigger("HLT_Ele15_LW_L1R", iEvent) ){
    }
    else break;
    if ( has_electron("selectedLayer1Electrons", iEvent) ){
      c_1->count();
    }
    else break;
    if ( is_electron("selectedLayer1Electrons",
                     20.0,
                     2.1,
                     0.02,
                     0.1,
                     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;
  cout << endl;
  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();
    }
  }
  result = passed_muons.getCount()==0;
  cout << endl;
  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());
    if (_pt > mPt &&
        _abseta < mEta &&
        _reliso < mRelIso &&
        _index != electron_index
        ){
      passed_electrons.count();
    }
  }
  result = passed_electrons.getCount()==0;
  cout << endl;
  return result;
}


bool CutFlow::has_electron(std::string mLabel,
                           const edm::Event& iEvent){
  bool result = false;
  Handle< vector< pat::Electron > >     electrons;
  iEvent . getByLabel( mLabel, electrons );
  result = electrons->size()>0;
  return result;
}


bool CutFlow::is_electron(std::string mLabel,
                      double mPt,
                      double mEta, 
                      double mD0,
                      double mRelIso,
                      const edm::Event& iEvent){
  bool result = false;
  cout << "Cut Flow: passed stage 1: ";
  Handle< vector< pat::Electron > >     electrons;
  iEvent . getByLabel( "selectedLayer1Electrons", electrons );
  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 electrons _______________________________________________
  //
  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_electrons("");
  for ( vector<pat::Electron>::const_iterator e = electrons -> begin(); e != electrons -> end(); e++){
    //_____ check that the electron is global
    if( !(e->isGlobalElectron()) ) continue;
    _pt = e->pt();
    _abseta = fabs(e->eta());
    //_d0 = e->dB(); // global track
    _d0 = -(e->innerTrack()->dxy(_bs));
    _chi2ndof = e->normChi2();
    _nhits = e->numberOfValidHits();
    reco::TrackRef _track = e->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 ( !(e->isIsolationValid()) ) continue;
    _emveto = e->isolationR03().emVetoEt;
    _hadveto = e->isolationR03().hadVetoEt;
    _trackiso = e->isolationR03().sumPt;
    _caloiso = e->isolationR03().emEt + e->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_electrons.count();
      electron_index = (int)(e - electrons->begin());
      if (passed_electrons.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_electrons.getCount()==1;
  return result;
}
  

//define this as a plug-in
DEFINE_FWK_MODULE(CutFlow);
Revision:	1.3
Committed:	Mon Oct 12 22:45:48 2009 UTC (15 years, 6 months ago) by kukartse
Content type:	text/plain
Branch:	MAIN
Changes since 1.2:	+158 -15 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.3	// $Id: CutFlow.cc,v 1.2 2009/10/10 05:26:18 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	kukartse	1.3	//
67			//_____ mu+jets cut flow
68			//
69			/*
70	kukartse	1.1	while(1){
71	kukartse	1.2	if ( is_trigger("HLT_Mu9", iEvent) ){
72			}
73			else break;
74			if ( has_global_muon("selectedLayer1Muons", iEvent) ){
75			c_1->count();
76			}
77	kukartse	1.1	else break;
78			if ( is_muon("selectedLayer1Muons",
79			20.0,
80			2.1,
81			11,
82			0.02,
83			10,
84			4.0,
85			6.0,
86			0.05,
87			iEvent) ){
88			c_2->count();
89			}
90			else break;
91	kukartse	1.2	if ( is_jets("selectedLayer1Jets",
92			30.0,
93			2.4,
94			iEvent) ){
95			c_3->count();
96			}
97			else break;
98			if ( no_loose_muon("selectedLayer1Muons",
99			2.5,
100			10.0,
101			0.2,
102			iEvent) ){
103			c_4->count();
104			}
105			else break;
106			if ( no_loose_electron("selectedLayer1Electrons",
107			2.5,
108			15.0,
109			0.2,
110			iEvent) ){
111			c_5->count();
112			}
113			else break;
114			cout << "Cut Flow: passed stage 5 " << endl;
115	kukartse	1.1	break;
116			}
117	kukartse	1.3	*/
118			//
119			//_____ e+jets cut flow
120			//
121			while(1){
122			if ( is_trigger("HLT_Ele15_LW_L1R", iEvent) ){
123			}
124			else break;
125			if ( has_electron("selectedLayer1Electrons", iEvent) ){
126			c_1->count();
127			}
128			else break;
129			if ( is_electron("selectedLayer1Electrons",
130			20.0,
131			2.1,
132			0.02,
133			0.1,
134			iEvent) ){
135			c_2->count();
136			}
137			else break;
138			if ( is_jets("selectedLayer1Jets",
139			30.0,
140			2.4,
141			iEvent) ){
142			c_3->count();
143			}
144			else break;
145			if ( no_loose_muon("selectedLayer1Muons",
146			2.5,
147			10.0,
148			0.2,
149			iEvent) ){
150			c_4->count();
151			}
152			else break;
153			if ( no_loose_electron("selectedLayer1Electrons",
154			2.5,
155			15.0,
156			0.2,
157			iEvent) ){
158			c_5->count();
159			}
160			else break;
161			cout << "Cut Flow: passed stage 5 " << endl;
162			break;
163			}
164	kukartse	1.1	}
165
166
167			void
168			CutFlow::beginJob()
169			{
170			}
171
172
173			void
174			CutFlow::endJob() {
175			cout << "========= Cut Flow Report ==========================================" << endl;
176			cout << "Total events: " << c_total->getCount() << endl;
177			cout << "Pass stage 1 (trigger+global): " << c_1->getCount() << endl;
178			cout << "Pass stage 2 (muon kinematics): " << c_2->getCount() << endl;
179			cout << "Pass stage 3 (jets): " << c_3->getCount() << endl;
180			cout << "Pass stage 4 (no loose muon): " << c_4->getCount() << endl;
181			cout << "Pass stage 5 (no loose electron): " << c_5->getCount() << endl;
182			cout << "================================================================" << endl;
183			}
184
185	kukartse	1.2	bool CutFlow::is_trigger(std::string mTriggerName, const edm::Event& iEvent){
186	kukartse	1.1	bool result = false;
187			Handle<std::vector<pat::TriggerPath> > trigs;
188			iEvent.getByLabel("patTrigger", trigs);
189			for (std::vector<pat::TriggerPath>::const_iterator trig = trigs->begin();
190			trig!=trigs->end();
191			++trig){
192			//cout << "Trigger name and number: " << trig->name() << ", " << trig->index() << endl;
193	kukartse	1.2	if( trig->name().find(mTriggerName)!=std::string::npos ){
194	kukartse	1.1	result = trig->wasAccept();
195			}
196			}
197			return result;
198			}
199
200	kukartse	1.2	bool CutFlow::has_global_muon(std::string mLabel,
201			const edm::Event& iEvent){
202			bool result = false;
203			cout << "Cut Flow: passed trigger ";
204			Handle< vector< pat::Muon > > muons;
205			iEvent . getByLabel( mLabel, muons );
206			//
207			//_____ loop over muons _______________________________________________
208			//
209			MeanCounter passed_muons("");
210			for ( vector<pat::Muon>::const_iterator mu = muons -> begin(); mu != muons -> end(); mu++){
211			if ( mu->isGlobalMuon() ){
212			passed_muons.count();
213			}
214			}
215			result = passed_muons.getCount()>0;
216	kukartse	1.3	cout << endl;
217	kukartse	1.2	return result;
218			}
219
220
221	kukartse	1.1	bool CutFlow::is_muon(std::string mLabel,
222			double mPt,
223			double mEta,
224			int mNHits,
225			double mD0,
226			double mChi2Ndof,
227			double mEmVeto,
228			double mHadVeto,
229			double mRelIso,
230			const edm::Event& iEvent){
231			bool result = false;
232	kukartse	1.2	cout << "Cut Flow: passed stage 1: ";
233	kukartse	1.1	Handle< vector< pat::Muon > > muons;
234			iEvent . getByLabel( "selectedLayer1Muons", muons );
235			Handle< reco::BeamSpot > beamSpotHandle;
236			iEvent . getByLabel( "offlineBeamSpot", beamSpotHandle);
237			//
238			//_____ Beam spot _____________________________________________________
239			//
240			reco::BeamSpot beamSpot;
241			if ( beamSpotHandle.isValid() ){
242			beamSpot = *beamSpotHandle;
243			}
244			else{
245			edm::LogInfo("TtLJetsAnalyzer")
246			<< "No beam spot available from EventSetup \n";
247			}
248			const reco::BeamSpot::Point _bs(beamSpot.x0(), beamSpot.y0(),beamSpot.z0());
249			//
250			//_____ loop over muons _______________________________________________
251			//
252			double _pt = 0.0;
253			double _abseta = 0.0;
254			int _nhits = 0;
255			double _d0 = 0.0;
256			double _chi2ndof = 0.0;
257			double _emveto = 10000.0;
258			double _hadveto = 10000.0;
259			double _trackiso = 10000.0;
260			double _caloiso = 10000.0;
261			double _reliso = 10000.0;
262	kukartse	1.2	MeanCounter passed_muons("");
263	kukartse	1.1	for ( vector<pat::Muon>::const_iterator mu = muons -> begin(); mu != muons -> end(); mu++){
264	kukartse	1.2	//_____ check that the muon is global
265			if( !(mu->isGlobalMuon()) ) continue;
266	kukartse	1.1	_pt = mu->pt();
267			_abseta = fabs(mu->eta());
268	kukartse	1.2	//_d0 = mu->dB(); // global track
269			_d0 = -(mu->innerTrack()->dxy(_bs));
270	kukartse	1.1	_chi2ndof = mu->normChi2();
271			_nhits = mu->numberOfValidHits();
272			reco::TrackRef _track = mu->track();
273			//_____ check that there is a track
274			if(_track.isNull()) continue;
275			//_nhits = _track->numberOfValidHits();
276			//_d0 = -(_track->dxy(_bs));
277			//_chi2ndof = _track->normalizedChi2();
278			//_chi2ndof = (_track->chi2())/(_track->ndof());
279			//_____ check that isolation is valid
280			if ( !(mu->isIsolationValid()) ) continue;
281			_emveto = mu->isolationR03().emVetoEt;
282			_hadveto = mu->isolationR03().hadVetoEt;
283			_trackiso = mu->isolationR03().sumPt;
284			_caloiso = mu->isolationR03().emEt + mu->isolationR03().hadEt;
285			_reliso = (_trackiso + _caloiso)/_pt;
286			if (_pt > mPt &&
287			_abseta < mEta &&
288			_nhits >= mNHits &&
289	kukartse	1.2	fabs(_d0) < mD0 &&
290	kukartse	1.1	_chi2ndof < mChi2Ndof &&
291	kukartse	1.2	fabs(_emveto) < mEmVeto &&
292			fabs(_hadveto) < mHadVeto &&
293			fabs(_reliso) < mRelIso
294	kukartse	1.1	){
295	kukartse	1.2	passed_muons.count();
296			muon_index = (int)(mu - muons->begin());
297			if (passed_muons.getCount()==1){
298			cout << _pt << ", ";
299			cout << _abseta << ", ";
300			cout << _nhits << ", ";
301			cout << _d0 << ", ";
302			cout << _chi2ndof << ", ";
303			cout << _emveto << ", ";
304			cout << _hadveto << ", ";
305			cout << _trackiso << ", ";
306			cout << _caloiso << ", ";
307			cout << _reliso << ", ";
308			cout << endl;
309			}
310			}
311			}
312			result = passed_muons.getCount()==1;
313			return result;
314			}
315
316
317			bool CutFlow::is_jets(std::string mLabel,
318			double mPt,
319			double mEta,
320			const edm::Event& iEvent){
321			bool result = false;
322
323			cout << "Cut Flow: passed stage 2: ";
324
325			Handle< vector< pat::Jet > > jets;
326			iEvent . getByLabel( mLabel, jets );
327			MeanCounter c_jets("");
328			MeanCounter passed_jets("");
329			for ( vector<pat::Jet>::const_iterator jet = jets -> begin(); jet != jets -> end(); jet++ ){
330			c_jets.count();
331			double _pt = jet->pt();
332			double _eta = jet->eta();
333			if (
334			( fabs( _eta ) < mEta ) &&
335			( _pt > mPt )
336			)
337			{
338			passed_jets.count();
339			}
340			if (c_jets.getCount()<=4){
341			cout << _pt << ", ";
342			cout << _eta << ", ";
343	kukartse	1.1	}
344			}
345			cout << endl;
346	kukartse	1.2	result = passed_jets.getCount() >=4;
347	kukartse	1.1	return result;
348			}
349	kukartse	1.2
350
351			bool CutFlow::no_loose_muon(std::string mLabel,
352			double mEta,
353			double mPt,
354			double mRelIso,
355			const edm::Event& iEvent){
356			bool result = false;
357			cout << "Cut Flow: passed stage 3: ";
358			Handle< vector< pat::Muon > > muons;
359			iEvent . getByLabel( "selectedLayer1Muons", muons );
360			//
361			//_____ loop over muons _______________________________________________
362			//
363			double _pt = 0.0;
364			double _abseta = 0.0;
365			double _trackiso = 10000.0;
366			double _caloiso = 10000.0;
367			double _reliso = 10000.0;
368			MeanCounter passed_muons("");
369			for ( vector<pat::Muon>::const_iterator mu = muons -> begin(); mu != muons -> end(); mu++){
370			//_____ check that the muon is global
371			if( !(mu->isGlobalMuon()) ) continue;
372			_pt = mu->pt();
373			_abseta = fabs(mu->eta());
374			//_____ check that isolation is valid
375			if ( !(mu->isIsolationValid()) ) continue;
376			_trackiso = mu->isolationR03().sumPt;
377			_caloiso = mu->isolationR03().emEt + mu->isolationR03().hadEt;
378			_reliso = (_trackiso + _caloiso)/_pt;
379			int _index = (int)(mu - muons->begin());
380			if (_pt > mPt &&
381			_abseta < mEta &&
382			_reliso < mRelIso &&
383			_index != muon_index
384			){
385			passed_muons.count();
386			}
387			}
388			result = passed_muons.getCount()==0;
389	kukartse	1.3	cout << endl;
390	kukartse	1.2	return result;
391			}
392
393
394			bool CutFlow::no_loose_electron(std::string mLabel,
395			double mEta,
396			double mPt,
397			double mRelIso,
398			const edm::Event& iEvent){
399			bool result = false;
400			cout << "Cut Flow: passed stage 4: ";
401			Handle< vector< pat::Electron > > electrons;
402			iEvent . getByLabel( "selectedLayer1Electrons", electrons );
403			//
404			//_____ loop over electrons _______________________________________________
405			//
406			double _pt = 0.0;
407			double _abseta = 0.0;
408			double _trackiso = 10000.0;
409			double _caloiso = 10000.0;
410			double _reliso = 10000.0;
411			MeanCounter passed_electrons("");
412			for ( vector<pat::Electron>::const_iterator e = electrons -> begin(); e != electrons -> end(); e++){
413			_pt = e->pt();
414			_abseta = fabs(e->eta());
415			_trackiso = e->dr03IsolationVariables().tkSumPt;
416			_caloiso = e->dr03IsolationVariables().ecalRecHitSumEt + e->dr03IsolationVariables().hcalDepth1TowerSumEt + e->dr03IsolationVariables().hcalDepth2TowerSumEt;
417			_reliso = (_trackiso + _caloiso)/_pt;
418			int _index = (int)(e - electrons->begin());
419			if (_pt > mPt &&
420			_abseta < mEta &&
421			_reliso < mRelIso &&
422			_index != electron_index
423			){
424			passed_electrons.count();
425			}
426			}
427			result = passed_electrons.getCount()==0;
428	kukartse	1.3	cout << endl;
429			return result;
430			}
431
432
433			bool CutFlow::has_electron(std::string mLabel,
434			const edm::Event& iEvent){
435			bool result = false;
436			Handle< vector< pat::Electron > > electrons;
437			iEvent . getByLabel( mLabel, electrons );
438			result = electrons->size()>0;
439	kukartse	1.2	return result;
440			}
441
442	kukartse	1.1
443	kukartse	1.3	bool CutFlow::is_electron(std::string mLabel,
444			double mPt,
445			double mEta,
446			double mD0,
447			double mRelIso,
448			const edm::Event& iEvent){
449			bool result = false;
450			cout << "Cut Flow: passed stage 1: ";
451			Handle< vector< pat::Electron > > electrons;
452			iEvent . getByLabel( "selectedLayer1Electrons", electrons );
453			Handle< reco::BeamSpot > beamSpotHandle;
454			iEvent . getByLabel( "offlineBeamSpot", beamSpotHandle);
455			//
456			//_____ Beam spot _____________________________________________________
457			//
458			reco::BeamSpot beamSpot;
459			if ( beamSpotHandle.isValid() ){
460			beamSpot = *beamSpotHandle;
461			}
462			else{
463			edm::LogInfo("TtLJetsAnalyzer")
464			<< "No beam spot available from EventSetup \n";
465			}
466			const reco::BeamSpot::Point _bs(beamSpot.x0(), beamSpot.y0(),beamSpot.z0());
467			//
468			//_____ loop over electrons _______________________________________________
469			//
470			double _pt = 0.0;
471			double _abseta = 0.0;
472			int _nhits = 0;
473			double _d0 = 0.0;
474			double _chi2ndof = 0.0;
475			double _emveto = 10000.0;
476			double _hadveto = 10000.0;
477			double _trackiso = 10000.0;
478			double _caloiso = 10000.0;
479			double _reliso = 10000.0;
480			MeanCounter passed_electrons("");
481			for ( vector<pat::Electron>::const_iterator e = electrons -> begin(); e != electrons -> end(); e++){
482			//_____ check that the electron is global
483			if( !(e->isGlobalElectron()) ) continue;
484			_pt = e->pt();
485			_abseta = fabs(e->eta());
486			//_d0 = e->dB(); // global track
487			_d0 = -(e->innerTrack()->dxy(_bs));
488			_chi2ndof = e->normChi2();
489			_nhits = e->numberOfValidHits();
490			reco::TrackRef _track = e->track();
491			//_____ check that there is a track
492			if(_track.isNull()) continue;
493			//_nhits = _track->numberOfValidHits();
494			//_d0 = -(_track->dxy(_bs));
495			//_chi2ndof = _track->normalizedChi2();
496			//_chi2ndof = (_track->chi2())/(_track->ndof());
497			//_____ check that isolation is valid
498			if ( !(e->isIsolationValid()) ) continue;
499			_emveto = e->isolationR03().emVetoEt;
500			_hadveto = e->isolationR03().hadVetoEt;
501			_trackiso = e->isolationR03().sumPt;
502			_caloiso = e->isolationR03().emEt + e->isolationR03().hadEt;
503			_reliso = (_trackiso + _caloiso)/_pt;
504			if (_pt > mPt &&
505			_abseta < mEta &&
506			_nhits >= mNHits &&
507			fabs(_d0) < mD0 &&
508			_chi2ndof < mChi2Ndof &&
509			fabs(_emveto) < mEmVeto &&
510			fabs(_hadveto) < mHadVeto &&
511			fabs(_reliso) < mRelIso
512			){
513			passed_electrons.count();
514			electron_index = (int)(e - electrons->begin());
515			if (passed_electrons.getCount()==1){
516			cout << _pt << ", ";
517			cout << _abseta << ", ";
518			cout << _nhits << ", ";
519			cout << _d0 << ", ";
520			cout << _chi2ndof << ", ";
521			cout << _emveto << ", ";
522			cout << _hadveto << ", ";
523			cout << _trackiso << ", ";
524			cout << _caloiso << ", ";
525			cout << _reliso << ", ";
526			cout << endl;
527			}
528			}
529			}
530			result = passed_electrons.getCount()==1;
531			return result;
532			}
533
534
535
536	kukartse	1.1	//define this as a plug-in
537			DEFINE_FWK_MODULE(CutFlow);