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 *
#	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.2 2009/10/10 05:26:18 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
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	muon_index = -1;
47	electron_index = -1;
48	}
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	//
67	//_____ mu+jets cut flow
68	//
69	/*
70	while(1){
71	if ( is_trigger("HLT_Mu9", iEvent) ){
72	}
73	else break;
74	if ( has_global_muon("selectedLayer1Muons", iEvent) ){
75	c_1->count();
76	}
77	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	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	break;
116	}
117	*/
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	}
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	bool CutFlow::is_trigger(std::string mTriggerName, const edm::Event& iEvent){
186	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	if( trig->name().find(mTriggerName)!=std::string::npos ){
194	result = trig->wasAccept();
195	}
196	}
197	return result;
198	}
199
200	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	cout << endl;
217	return result;
218	}
219
220
221	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	cout << "Cut Flow: passed stage 1: ";
233	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	MeanCounter passed_muons("");
263	for ( vector<pat::Muon>::const_iterator mu = muons -> begin(); mu != muons -> end(); mu++){
264	//_____ check that the muon is global
265	if( !(mu->isGlobalMuon()) ) continue;
266	_pt = mu->pt();
267	_abseta = fabs(mu->eta());
268	//_d0 = mu->dB(); // global track
269	_d0 = -(mu->innerTrack()->dxy(_bs));
270	_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	fabs(_d0) < mD0 &&
290	_chi2ndof < mChi2Ndof &&
291	fabs(_emveto) < mEmVeto &&
292	fabs(_hadveto) < mHadVeto &&
293	fabs(_reliso) < mRelIso
294	){
295	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	}
344	}
345	cout << endl;
346	result = passed_jets.getCount() >=4;
347	return result;
348	}
349
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	cout << endl;
390	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	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	return result;
440	}
441
442
443	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	//define this as a plug-in
537	DEFINE_FWK_MODULE(CutFlow);