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.3 2009/10/12 22:45: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"
#include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
#include "DataFormats/Math/interface/deltaR.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){
  //
  //_____ stage 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){
    //
    //_____ stage 1
    //
    if ( is_trigger("HLT_Ele15_LW_L1R", iEvent) ){
    }
    else break;
    if ( has_electron("selectedLayer1Electrons", iEvent) ){
      c_1->count();
    }
    else break;
    //
    //_____ stage 2
    //
    cout << "CutFlow: passed stage 1: ";
    if ( is_electron("selectedLayer1Electrons",
                     30.0,
                     2.4,
                     0.02,
                     0.1,
                     iEvent) ){
      c_2->count();
      cout << endl;
    }
    else{
      cout << endl;
      break;
    }      
    //
    //_____ stage 3
    //
    cout << "CutFlow: passed stage 2: ";
    if ( !is_muon("selectedLayer1Muons",
                  30.0,
                  2.1,
                  11,
                  0.02,
                  10,
                  4000000.0,
                  6000000.0,
                  0.05,
                  iEvent) ){
      c_3->count();
      cout << endl;
    }
    else{
      cout << endl;
      break;
    }
    //
    //_____ stage 4
    //
    cout << "CutFlow: passed stage 3: ";
    int n_jets = count_jets("selectedLayer1Jets",
                            30.0,
                            2.4,
                            "selectedLayer1Electrons",
                            30.0,
                            2.4,
                            0.02,
                            0.3,
                            iEvent);
    if ( n_jets >= 0){
      c_4->count();
      cout << endl;
    }
    else{
      cout << endl;
      break;
    }
    //
    //_____ stage 5
    //
    cout << "CutFlow: passed stage 4: ";
    if ( n_jets >= 4 ){
      c_5->count();
      cout << endl;
    }
    else{
      cout << endl;
      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: " << c_1->getCount() << endl;
  cout << "Pass stage 2: " << c_2->getCount() << endl;
  cout << "Pass stage 3: " << c_3->getCount() << endl;
  cout << "Pass stage 4: " << c_4->getCount() << endl;
  cout << "Pass stage 5: " << 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 (mu+jets) or stage 2 (e+jets): ";
  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;
      }
    }
  }
  cout << "Number of good muons: " << passed_muons.getCount();
  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;
  double _d0 = 0.0;
  double _d0_pat = 0.0;
  double _trackiso = 10000.0;
  double _ecaloiso = 10000.0;
  double _ecaloiso2 = 10000.0;
  double _pat_ecaloiso = 10000.0;
  double _hcaloiso1 = 10000.0;
  double _hcaloiso2 = 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 robust-tight
    if( !(e->electronID("eidRobustTight")) ) continue;
    _pt = e->pt();
    _abseta = fabs(e->eta());
    _d0_pat = e->dB();
    _d0 = -(e->gsfTrack()->dxy(_bs));
    _trackiso = e->dr03IsolationVariables().tkSumPt;
    _ecaloiso = e->dr04IsolationVariables().ecalRecHitSumEt;
    _ecaloiso2 = e->dr03IsolationVariables().ecalRecHitSumEt;
    _pat_ecaloiso = e->isolation(pat::ECalIso);
    _hcaloiso1 = e->dr04IsolationVariables().hcalDepth1TowerSumEt;
    _hcaloiso2 = e->dr04IsolationVariables().hcalDepth2TowerSumEt;
    _caloiso = _ecaloiso + _hcaloiso1 + _hcaloiso2;
    //_caloiso = _pat_ecaloiso +_hcaloiso1 + _hcaloiso2;
    _reliso = (_trackiso + _caloiso)/_pt;
    if (_pt > mPt &&
        _abseta < mEta &&
        fabs(_d0) < mD0 &&
        fabs(_reliso) < mRelIso
        ){
      passed_electrons.count();
      electron_index = (int)(e - electrons->begin());
      if (passed_electrons.getCount()==1){
        cout << "pt=" << _pt << ", abseta=";
        cout << _abseta << ", d0=";
        cout << _d0 << ", trackiso=";
        cout << _trackiso << ", ecaloiso04=";
        cout << _ecaloiso << ", ecaloiso03=";
        cout << _ecaloiso2 << ", pat_ecaloiso=";
        cout << _pat_ecaloiso << ", hcaloiso1=";
        cout << _hcaloiso1 << ", hcaloiso2=";
        cout << _hcaloiso2 << ", caloiso04=";
        cout << _caloiso << ", reliso=";
        cout << _reliso;
        //cout << endl;
      }
    }
  }
  //cout << "Number of good electrons: " << passed_electrons.getCount();
  result = passed_electrons.getCount()==1;
  return result;
}
  

int CutFlow::count_jets(std::string jLabel,
                        double jPt,
                        double jEta,
                        std::string lLabel,
                        double lPt,
                        double lEta,
                        double lD0,
                        double lDR,
                        const edm::Event & iEvent){
  MeanCounter _njets("");
  Handle< vector< pat::Jet > >     jets;
  iEvent . getByLabel( jLabel, jets );
  Handle< vector< pat::Electron > >     electrons;
  iEvent . getByLabel( lLabel, 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 jets _______________________________________________
  //
  double j_pt = 0.0;
  double j_eta = 0.0;
  for ( vector<pat::Jet>::const_iterator jet = jets -> begin(); jet != jets -> end(); ++jet){
    j_pt = jet->pt();
    j_eta = jet->eta();
    //
    //_____ loop over electrons _______________________________________________
    //
    double e_pt = 0.0;
    double e_abseta = 0.0;
    double e_d0 = 0.0;
    double e_dr = 0.0;
    bool jet_removed = false;
    for ( vector<pat::Electron>::const_iterator e = electrons -> begin(); e != electrons -> end(); e++){
      //_____ check that the electron is robust-tight
      if( !(e->electronID("eidRobustTight")) ) continue;
      e_pt = e->pt();
      e_abseta = fabs(e->eta());
      e_d0 = -(e->gsfTrack()->dxy(_bs));
      e_dr = reco::deltaR(*jet, *e);
      if ( e_pt > lPt && e_abseta < lEta && fabs(e_d0)<lD0 && e_dr<lDR){
        jet_removed = true;
      }
    }
    if (j_pt>jPt && fabs(j_eta)<jEta && !jet_removed) _njets.count();
  }
  cout << "Number of jets after cleaning: " << _njets.getCount();
  return _njets.getCount();
}

//define this as a plug-in
DEFINE_FWK_MODULE(CutFlow);
Revision:	1.4
Committed:	Tue Oct 27 22:39:44 2009 UTC (15 years, 6 months ago) by kukartse
Content type:	text/plain
Branch:	MAIN
CVS Tags:	gak110409
Changes since 1.3:	+179 -78 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.4	// $Id: CutFlow.cc,v 1.3 2009/10/12 22:45:48 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.4	#include "DataFormats/Math/interface/deltaR.h"
31	kukartse	1.1
32			using namespace edm;
33
34			CutFlow::CutFlow(const edm::ParameterSet& iConfig){
35			c_total = new MeanCounter("Processing event number: ", 0, 1);
36			c_1 = new MeanCounter("CutFlow stage 1: ", 0, 1);
37			c_2 = new MeanCounter("CutFlow stage 2: ", 0, 1);
38			c_3 = new MeanCounter("CutFlow stage 3: ", 0, 1);
39			c_4 = new MeanCounter("CutFlow stage 4: ", 0, 1);
40			c_5 = new MeanCounter("CutFlow stage 5: ", 0, 1);
41			c_1->setPrintCount(false);
42			c_2->setPrintCount(false);
43			c_3->setPrintCount(false);
44			c_4->setPrintCount(false);
45			c_5->setPrintCount(false);
46	kukartse	1.2
47			muon_index = -1;
48			electron_index = -1;
49	kukartse	1.1	}
50
51
52			CutFlow::~CutFlow()
53			{
54			delete c_total;
55			delete c_1;
56			delete c_2;
57			delete c_3;
58			delete c_4;
59			delete c_5;
60			}
61
62
63			void
64			CutFlow::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
65			{
66			c_total->count();
67	kukartse	1.3	//
68			//_____ mu+jets cut flow
69			//
70			/*
71	kukartse	1.1	while(1){
72	kukartse	1.4	//
73			//_____ stage 1
74			//
75			if ( is_trigger("HLT_Mu9", iEvent) ){
76	kukartse	1.2	}
77			else break;
78			if ( has_global_muon("selectedLayer1Muons", iEvent) ){
79			c_1->count();
80			}
81	kukartse	1.1	else break;
82	kukartse	1.4
83	kukartse	1.1	if ( is_muon("selectedLayer1Muons",
84			20.0,
85			2.1,
86			11,
87			0.02,
88			10,
89			4.0,
90			6.0,
91			0.05,
92			iEvent) ){
93			c_2->count();
94			}
95			else break;
96	kukartse	1.2	if ( is_jets("selectedLayer1Jets",
97			30.0,
98			2.4,
99			iEvent) ){
100			c_3->count();
101			}
102			else break;
103			if ( no_loose_muon("selectedLayer1Muons",
104			2.5,
105			10.0,
106			0.2,
107			iEvent) ){
108			c_4->count();
109			}
110			else break;
111			if ( no_loose_electron("selectedLayer1Electrons",
112			2.5,
113			15.0,
114			0.2,
115			iEvent) ){
116			c_5->count();
117			}
118			else break;
119			cout << "Cut Flow: passed stage 5 " << endl;
120	kukartse	1.1	break;
121			}
122	kukartse	1.3	*/
123			//
124			//_____ e+jets cut flow
125			//
126			while(1){
127	kukartse	1.4	//
128			//_____ stage 1
129			//
130	kukartse	1.3	if ( is_trigger("HLT_Ele15_LW_L1R", iEvent) ){
131			}
132			else break;
133			if ( has_electron("selectedLayer1Electrons", iEvent) ){
134			c_1->count();
135			}
136			else break;
137	kukartse	1.4	//
138			//_____ stage 2
139			//
140			cout << "CutFlow: passed stage 1: ";
141	kukartse	1.3	if ( is_electron("selectedLayer1Electrons",
142	kukartse	1.4	30.0,
143			2.4,
144	kukartse	1.3	0.02,
145			0.1,
146			iEvent) ){
147			c_2->count();
148	kukartse	1.4	cout << endl;
149	kukartse	1.3	}
150	kukartse	1.4	else{
151			cout << endl;
152			break;
153			}
154			//
155			//_____ stage 3
156			//
157			cout << "CutFlow: passed stage 2: ";
158			if ( !is_muon("selectedLayer1Muons",
159			30.0,
160			2.1,
161			11,
162			0.02,
163			10,
164			4000000.0,
165			6000000.0,
166			0.05,
167			iEvent) ){
168	kukartse	1.3	c_3->count();
169	kukartse	1.4	cout << endl;
170	kukartse	1.3	}
171	kukartse	1.4	else{
172			cout << endl;
173			break;
174			}
175			//
176			//_____ stage 4
177			//
178			cout << "CutFlow: passed stage 3: ";
179			int n_jets = count_jets("selectedLayer1Jets",
180			30.0,
181			2.4,
182			"selectedLayer1Electrons",
183			30.0,
184			2.4,
185			0.02,
186			0.3,
187			iEvent);
188			if ( n_jets >= 0){
189	kukartse	1.3	c_4->count();
190	kukartse	1.4	cout << endl;
191	kukartse	1.3	}
192	kukartse	1.4	else{
193			cout << endl;
194			break;
195			}
196			//
197			//_____ stage 5
198			//
199			cout << "CutFlow: passed stage 4: ";
200			if ( n_jets >= 4 ){
201	kukartse	1.3	c_5->count();
202	kukartse	1.4	cout << endl;
203			}
204			else{
205			cout << endl;
206			break;
207	kukartse	1.3	}
208			cout << "Cut Flow: passed stage 5 " << endl;
209			break;
210			}
211	kukartse	1.1	}
212
213
214			void
215			CutFlow::beginJob()
216			{
217			}
218
219
220			void
221			CutFlow::endJob() {
222			cout << "========= Cut Flow Report ==========================================" << endl;
223			cout << "Total events: " << c_total->getCount() << endl;
224	kukartse	1.4	cout << "Pass stage 1: " << c_1->getCount() << endl;
225			cout << "Pass stage 2: " << c_2->getCount() << endl;
226			cout << "Pass stage 3: " << c_3->getCount() << endl;
227			cout << "Pass stage 4: " << c_4->getCount() << endl;
228			cout << "Pass stage 5: " << c_5->getCount() << endl;
229	kukartse	1.1	cout << "================================================================" << endl;
230			}
231
232	kukartse	1.2	bool CutFlow::is_trigger(std::string mTriggerName, const edm::Event& iEvent){
233	kukartse	1.1	bool result = false;
234			Handle<std::vector<pat::TriggerPath> > trigs;
235			iEvent.getByLabel("patTrigger", trigs);
236			for (std::vector<pat::TriggerPath>::const_iterator trig = trigs->begin();
237			trig!=trigs->end();
238			++trig){
239			//cout << "Trigger name and number: " << trig->name() << ", " << trig->index() << endl;
240	kukartse	1.2	if( trig->name().find(mTriggerName)!=std::string::npos ){
241	kukartse	1.1	result = trig->wasAccept();
242			}
243			}
244			return result;
245			}
246
247	kukartse	1.2	bool CutFlow::has_global_muon(std::string mLabel,
248			const edm::Event& iEvent){
249			bool result = false;
250	kukartse	1.4	//cout << "Cut Flow: passed trigger ";
251	kukartse	1.2	Handle< vector< pat::Muon > > muons;
252			iEvent . getByLabel( mLabel, muons );
253			//
254			//_____ loop over muons _______________________________________________
255			//
256			MeanCounter passed_muons("");
257			for ( vector<pat::Muon>::const_iterator mu = muons -> begin(); mu != muons -> end(); mu++){
258			if ( mu->isGlobalMuon() ){
259			passed_muons.count();
260			}
261			}
262			result = passed_muons.getCount()>0;
263	kukartse	1.4	//cout << endl;
264	kukartse	1.2	return result;
265			}
266
267
268	kukartse	1.1	bool CutFlow::is_muon(std::string mLabel,
269			double mPt,
270			double mEta,
271			int mNHits,
272			double mD0,
273			double mChi2Ndof,
274			double mEmVeto,
275			double mHadVeto,
276			double mRelIso,
277			const edm::Event& iEvent){
278			bool result = false;
279	kukartse	1.4	//cout << "Cut Flow: passed stage 1 (mu+jets) or stage 2 (e+jets): ";
280	kukartse	1.1	Handle< vector< pat::Muon > > muons;
281			iEvent . getByLabel( "selectedLayer1Muons", muons );
282			Handle< reco::BeamSpot > beamSpotHandle;
283			iEvent . getByLabel( "offlineBeamSpot", beamSpotHandle);
284			//
285			//_____ Beam spot _____________________________________________________
286			//
287			reco::BeamSpot beamSpot;
288			if ( beamSpotHandle.isValid() ){
289			beamSpot = *beamSpotHandle;
290			}
291			else{
292			edm::LogInfo("TtLJetsAnalyzer")
293			<< "No beam spot available from EventSetup \n";
294			}
295			const reco::BeamSpot::Point _bs(beamSpot.x0(), beamSpot.y0(),beamSpot.z0());
296			//
297			//_____ loop over muons _______________________________________________
298			//
299			double _pt = 0.0;
300			double _abseta = 0.0;
301			int _nhits = 0;
302			double _d0 = 0.0;
303			double _chi2ndof = 0.0;
304			double _emveto = 10000.0;
305			double _hadveto = 10000.0;
306			double _trackiso = 10000.0;
307			double _caloiso = 10000.0;
308			double _reliso = 10000.0;
309	kukartse	1.2	MeanCounter passed_muons("");
310	kukartse	1.1	for ( vector<pat::Muon>::const_iterator mu = muons -> begin(); mu != muons -> end(); mu++){
311	kukartse	1.2	//_____ check that the muon is global
312			if( !(mu->isGlobalMuon()) ) continue;
313	kukartse	1.1	_pt = mu->pt();
314			_abseta = fabs(mu->eta());
315	kukartse	1.2	//_d0 = mu->dB(); // global track
316			_d0 = -(mu->innerTrack()->dxy(_bs));
317	kukartse	1.1	_chi2ndof = mu->normChi2();
318			_nhits = mu->numberOfValidHits();
319			reco::TrackRef _track = mu->track();
320			//_____ check that there is a track
321			if(_track.isNull()) continue;
322			//_nhits = _track->numberOfValidHits();
323			//_d0 = -(_track->dxy(_bs));
324			//_chi2ndof = _track->normalizedChi2();
325			//_chi2ndof = (_track->chi2())/(_track->ndof());
326			//_____ check that isolation is valid
327			if ( !(mu->isIsolationValid()) ) continue;
328			_emveto = mu->isolationR03().emVetoEt;
329			_hadveto = mu->isolationR03().hadVetoEt;
330			_trackiso = mu->isolationR03().sumPt;
331			_caloiso = mu->isolationR03().emEt + mu->isolationR03().hadEt;
332			_reliso = (_trackiso + _caloiso)/_pt;
333			if (_pt > mPt &&
334			_abseta < mEta &&
335			_nhits >= mNHits &&
336	kukartse	1.2	fabs(_d0) < mD0 &&
337	kukartse	1.1	_chi2ndof < mChi2Ndof &&
338	kukartse	1.2	fabs(_emveto) < mEmVeto &&
339			fabs(_hadveto) < mHadVeto &&
340			fabs(_reliso) < mRelIso
341	kukartse	1.1	){
342	kukartse	1.2	passed_muons.count();
343			muon_index = (int)(mu - muons->begin());
344			if (passed_muons.getCount()==1){
345			cout << _pt << ", ";
346			cout << _abseta << ", ";
347			cout << _nhits << ", ";
348			cout << _d0 << ", ";
349			cout << _chi2ndof << ", ";
350			cout << _emveto << ", ";
351			cout << _hadveto << ", ";
352			cout << _trackiso << ", ";
353			cout << _caloiso << ", ";
354	kukartse	1.4	cout << _reliso;
355			//cout << endl;
356	kukartse	1.2	}
357			}
358			}
359	kukartse	1.4	cout << "Number of good muons: " << passed_muons.getCount();
360	kukartse	1.2	result = passed_muons.getCount()==1;
361			return result;
362			}
363
364
365			bool CutFlow::is_jets(std::string mLabel,
366			double mPt,
367			double mEta,
368			const edm::Event& iEvent){
369			bool result = false;
370
371	kukartse	1.4	//cout << "Cut Flow: passed stage 2: ";
372	kukartse	1.2
373			Handle< vector< pat::Jet > > jets;
374			iEvent . getByLabel( mLabel, jets );
375			MeanCounter c_jets("");
376			MeanCounter passed_jets("");
377			for ( vector<pat::Jet>::const_iterator jet = jets -> begin(); jet != jets -> end(); jet++ ){
378			c_jets.count();
379			double _pt = jet->pt();
380			double _eta = jet->eta();
381			if (
382			( fabs( _eta ) < mEta ) &&
383			( _pt > mPt )
384			)
385			{
386			passed_jets.count();
387			}
388			if (c_jets.getCount()<=4){
389			cout << _pt << ", ";
390	kukartse	1.4	cout << _eta;
391	kukartse	1.1	}
392			}
393	kukartse	1.4	//cout << endl;
394	kukartse	1.2	result = passed_jets.getCount() >=4;
395	kukartse	1.1	return result;
396			}
397	kukartse	1.2
398
399			bool CutFlow::no_loose_muon(std::string mLabel,
400			double mEta,
401			double mPt,
402			double mRelIso,
403			const edm::Event& iEvent){
404			bool result = false;
405	kukartse	1.4	//cout << "Cut Flow: passed stage 3: ";
406	kukartse	1.2	Handle< vector< pat::Muon > > muons;
407			iEvent . getByLabel( "selectedLayer1Muons", muons );
408			//
409			//_____ loop over muons _______________________________________________
410			//
411			double _pt = 0.0;
412			double _abseta = 0.0;
413			double _trackiso = 10000.0;
414			double _caloiso = 10000.0;
415			double _reliso = 10000.0;
416			MeanCounter passed_muons("");
417			for ( vector<pat::Muon>::const_iterator mu = muons -> begin(); mu != muons -> end(); mu++){
418			//_____ check that the muon is global
419			if( !(mu->isGlobalMuon()) ) continue;
420			_pt = mu->pt();
421			_abseta = fabs(mu->eta());
422			//_____ check that isolation is valid
423			if ( !(mu->isIsolationValid()) ) continue;
424			_trackiso = mu->isolationR03().sumPt;
425			_caloiso = mu->isolationR03().emEt + mu->isolationR03().hadEt;
426			_reliso = (_trackiso + _caloiso)/_pt;
427			int _index = (int)(mu - muons->begin());
428			if (_pt > mPt &&
429			_abseta < mEta &&
430			_reliso < mRelIso &&
431			_index != muon_index
432			){
433			passed_muons.count();
434			}
435			}
436			result = passed_muons.getCount()==0;
437	kukartse	1.4	//cout << endl;
438	kukartse	1.2	return result;
439			}
440
441
442			bool CutFlow::no_loose_electron(std::string mLabel,
443			double mEta,
444			double mPt,
445			double mRelIso,
446			const edm::Event& iEvent){
447			bool result = false;
448	kukartse	1.4	//cout << "Cut Flow: passed stage 4: ";
449	kukartse	1.2	Handle< vector< pat::Electron > > electrons;
450			iEvent . getByLabel( "selectedLayer1Electrons", electrons );
451			//
452			//_____ loop over electrons _______________________________________________
453			//
454			double _pt = 0.0;
455			double _abseta = 0.0;
456			double _trackiso = 10000.0;
457			double _caloiso = 10000.0;
458			double _reliso = 10000.0;
459			MeanCounter passed_electrons("");
460			for ( vector<pat::Electron>::const_iterator e = electrons -> begin(); e != electrons -> end(); e++){
461			_pt = e->pt();
462			_abseta = fabs(e->eta());
463			_trackiso = e->dr03IsolationVariables().tkSumPt;
464			_caloiso = e->dr03IsolationVariables().ecalRecHitSumEt + e->dr03IsolationVariables().hcalDepth1TowerSumEt + e->dr03IsolationVariables().hcalDepth2TowerSumEt;
465			_reliso = (_trackiso + _caloiso)/_pt;
466			int _index = (int)(e - electrons->begin());
467			if (_pt > mPt &&
468			_abseta < mEta &&
469			_reliso < mRelIso &&
470			_index != electron_index
471			){
472			passed_electrons.count();
473			}
474			}
475			result = passed_electrons.getCount()==0;
476	kukartse	1.4	//cout << endl;
477	kukartse	1.3	return result;
478			}
479
480
481			bool CutFlow::has_electron(std::string mLabel,
482			const edm::Event& iEvent){
483			bool result = false;
484			Handle< vector< pat::Electron > > electrons;
485			iEvent . getByLabel( mLabel, electrons );
486			result = electrons->size()>0;
487	kukartse	1.2	return result;
488			}
489
490	kukartse	1.1
491	kukartse	1.3	bool CutFlow::is_electron(std::string mLabel,
492			double mPt,
493			double mEta,
494			double mD0,
495			double mRelIso,
496			const edm::Event& iEvent){
497			bool result = false;
498	kukartse	1.4	//cout << "Cut Flow: passed stage 1: ";
499	kukartse	1.3	Handle< vector< pat::Electron > > electrons;
500			iEvent . getByLabel( "selectedLayer1Electrons", electrons );
501			Handle< reco::BeamSpot > beamSpotHandle;
502			iEvent . getByLabel( "offlineBeamSpot", beamSpotHandle);
503			//
504			//_____ Beam spot _____________________________________________________
505			//
506			reco::BeamSpot beamSpot;
507			if ( beamSpotHandle.isValid() ){
508			beamSpot = *beamSpotHandle;
509			}
510			else{
511			edm::LogInfo("TtLJetsAnalyzer")
512			<< "No beam spot available from EventSetup \n";
513			}
514			const reco::BeamSpot::Point _bs(beamSpot.x0(), beamSpot.y0(),beamSpot.z0());
515			//
516			//_____ loop over electrons _______________________________________________
517			//
518			double _pt = 0.0;
519			double _abseta = 0.0;
520			double _d0 = 0.0;
521	kukartse	1.4	double _d0_pat = 0.0;
522	kukartse	1.3	double _trackiso = 10000.0;
523	kukartse	1.4	double _ecaloiso = 10000.0;
524			double _ecaloiso2 = 10000.0;
525			double _pat_ecaloiso = 10000.0;
526			double _hcaloiso1 = 10000.0;
527			double _hcaloiso2 = 10000.0;
528	kukartse	1.3	double _caloiso = 10000.0;
529			double _reliso = 10000.0;
530			MeanCounter passed_electrons("");
531			for ( vector<pat::Electron>::const_iterator e = electrons -> begin(); e != electrons -> end(); e++){
532	kukartse	1.4	//_____ check that the electron is robust-tight
533			if( !(e->electronID("eidRobustTight")) ) continue;
534	kukartse	1.3	_pt = e->pt();
535			_abseta = fabs(e->eta());
536	kukartse	1.4	_d0_pat = e->dB();
537			_d0 = -(e->gsfTrack()->dxy(_bs));
538			_trackiso = e->dr03IsolationVariables().tkSumPt;
539			_ecaloiso = e->dr04IsolationVariables().ecalRecHitSumEt;
540			_ecaloiso2 = e->dr03IsolationVariables().ecalRecHitSumEt;
541			_pat_ecaloiso = e->isolation(pat::ECalIso);
542			_hcaloiso1 = e->dr04IsolationVariables().hcalDepth1TowerSumEt;
543			_hcaloiso2 = e->dr04IsolationVariables().hcalDepth2TowerSumEt;
544			_caloiso = _ecaloiso + _hcaloiso1 + _hcaloiso2;
545			//_caloiso = _pat_ecaloiso +_hcaloiso1 + _hcaloiso2;
546	kukartse	1.3	_reliso = (_trackiso + _caloiso)/_pt;
547			if (_pt > mPt &&
548			_abseta < mEta &&
549			fabs(_d0) < mD0 &&
550			fabs(_reliso) < mRelIso
551			){
552			passed_electrons.count();
553			electron_index = (int)(e - electrons->begin());
554			if (passed_electrons.getCount()==1){
555	kukartse	1.4	cout << "pt=" << _pt << ", abseta=";
556			cout << _abseta << ", d0=";
557			cout << _d0 << ", trackiso=";
558			cout << _trackiso << ", ecaloiso04=";
559			cout << _ecaloiso << ", ecaloiso03=";
560			cout << _ecaloiso2 << ", pat_ecaloiso=";
561			cout << _pat_ecaloiso << ", hcaloiso1=";
562			cout << _hcaloiso1 << ", hcaloiso2=";
563			cout << _hcaloiso2 << ", caloiso04=";
564			cout << _caloiso << ", reliso=";
565			cout << _reliso;
566			//cout << endl;
567	kukartse	1.3	}
568			}
569			}
570	kukartse	1.4	//cout << "Number of good electrons: " << passed_electrons.getCount();
571	kukartse	1.3	result = passed_electrons.getCount()==1;
572			return result;
573			}
574
575
576	kukartse	1.4	int CutFlow::count_jets(std::string jLabel,
577			double jPt,
578			double jEta,
579			std::string lLabel,
580			double lPt,
581			double lEta,
582			double lD0,
583			double lDR,
584			const edm::Event & iEvent){
585			MeanCounter _njets("");
586			Handle< vector< pat::Jet > > jets;
587			iEvent . getByLabel( jLabel, jets );
588			Handle< vector< pat::Electron > > electrons;
589			iEvent . getByLabel( lLabel, electrons );
590			Handle< reco::BeamSpot > beamSpotHandle;
591			iEvent . getByLabel( "offlineBeamSpot", beamSpotHandle);
592			//
593			//_____ Beam spot _____________________________________________________
594			//
595			reco::BeamSpot beamSpot;
596			if ( beamSpotHandle.isValid() ){
597			beamSpot = *beamSpotHandle;
598			}
599			else{
600			edm::LogInfo("TtLJetsAnalyzer")
601			<< "No beam spot available from EventSetup \n";
602			}
603			const reco::BeamSpot::Point _bs(beamSpot.x0(), beamSpot.y0(),beamSpot.z0());
604			//
605			//_____ loop over jets _______________________________________________
606			//
607			double j_pt = 0.0;
608			double j_eta = 0.0;
609			for ( vector<pat::Jet>::const_iterator jet = jets -> begin(); jet != jets -> end(); ++jet){
610			j_pt = jet->pt();
611			j_eta = jet->eta();
612			//
613			//_____ loop over electrons _______________________________________________
614			//
615			double e_pt = 0.0;
616			double e_abseta = 0.0;
617			double e_d0 = 0.0;
618			double e_dr = 0.0;
619			bool jet_removed = false;
620			for ( vector<pat::Electron>::const_iterator e = electrons -> begin(); e != electrons -> end(); e++){
621			//_____ check that the electron is robust-tight
622			if( !(e->electronID("eidRobustTight")) ) continue;
623			e_pt = e->pt();
624			e_abseta = fabs(e->eta());
625			e_d0 = -(e->gsfTrack()->dxy(_bs));
626			e_dr = reco::deltaR(jet, e);
627			if ( e_pt > lPt && e_abseta < lEta && fabs(e_d0)<lD0 && e_dr<lDR){
628			jet_removed = true;
629			}
630			}
631			if (j_pt>jPt && fabs(j_eta)<jEta && !jet_removed) _njets.count();
632			}
633			cout << "Number of jets after cleaning: " << _njets.getCount();
634			return _njets.getCount();
635			}
636	kukartse	1.3
637	kukartse	1.1	//define this as a plug-in
638			DEFINE_FWK_MODULE(CutFlow);