HbbAnalyzer/plugins/HbbCandidateFinder.cc

#include "VHbbAnalysis/HbbAnalyzer/interface/HbbCandidateFinder.h"

struct CompareJetPt {
  bool operator()( const VHbbEvent::SimpleJet& j1, const  VHbbEvent::SimpleJet& j2 ) const {
    return j1.fourMomentum.Pt() > j2.fourMomentum.Pt();
  }
};

  struct CompareBTag {
    bool operator()(const  VHbbEvent::SimpleJet& j1, const  VHbbEvent::SimpleJet& j2 ) const {
      return j1.csv > j2.csv;
    }
  };

HbbCandidateFinder::HbbCandidateFinder(const edm::ParameterSet& iConfig):   vhbbevent_(iConfig.getParameter<edm::InputTag>("VHbbEventLabel")), verbose_ (iConfig.getParameter<bool>("verbose")), jetPtThreshold(iConfig.getParameter<double>("jetPtThreshold")){
  produces<std::vector<VHbbCandidate > >();
}

HbbCandidateFinder::~HbbCandidateFinder(){}

void HbbCandidateFinder::beginJob(){}
void HbbCandidateFinder::endJob(){}


float HbbCandidateFinder::getDeltaTheta( VHbbEvent::SimpleJet * j1, VHbbEvent::SimpleJet * j2 ){return -1.;}


void HbbCandidateFinder::produce( edm::Event& iEvent, const edm::EventSetup& iEventSetup){
  
  std::auto_ptr<std::vector<VHbbCandidate> >  vHbbCandidates( new std::vector<VHbbCandidate>  );

  edm::Handle<VHbbEvent>  vHbbEvent;
  //  iEvent.getByLabel(vhbbevent_, vHbbEvent);
  iEvent.getByType(vHbbEvent);
  

  //
  // start searching for candidates
  //

  //  hbbCandidateFinderAlgo(vHbbCandidates, vHbbEvent-> result());
  // do nothing for a test
  
  run(vHbbEvent.product(),vHbbCandidates);


  if (verbose_)
    std::cout <<" Pushing VHbb candidates: "<<vHbbCandidates->size()<<std::endl;

  iEvent.put(vHbbCandidates);  

}

void HbbCandidateFinder::run (const VHbbEvent* event, std::auto_ptr<std::vector<VHbbCandidate> > & candidates){
  //
  // first find the jets
  //

  VHbbEvent::SimpleJet j1,j2;
  std::vector<VHbbEvent::SimpleJet> addJets;
  bool foundJets = findDiJets(event->simpleJets2,j1,j2,addJets) ;

  if (verbose_){
    std::cout <<" Found Dijets: "<<foundJets<< " Additional: "<<addJets.size()<< std::endl;
  }
  
  if (foundJets == false) return;
  //
  // search for a dilepton - just 
  //
  std::vector<VHbbEvent::MuonInfo> mu;
  findMuons(event->muInfo,mu);
  std::vector<VHbbEvent::ElectronInfo> ele;
  findElectrons(event->eleInfo,ele);
  
  if (verbose_){
    std::cout <<" Electrons: "<< ele.size()<<std::endl;
    std::cout <<" Muons    : "<< mu.size()<<std::endl;
  }
  if (ele.size()<1 && mu.size() < 1 ) return;

  //
  // fill!
  //
  VHbbCandidate tempMu;
  VHbbCandidate tempE;
  tempMu.H.jets.push_back(j1);
  tempMu.H.jets.push_back(j2);
  tempMu.H.fourMomentum = (j1).fourMomentum+(j2).fourMomentum;
  tempMu.additionalJets = addJets;
  tempE = tempMu;
  TLorentzVector pMu,pE;

  if (mu.size()){
    for (std::vector<VHbbEvent::MuonInfo>::iterator it = mu.begin(); it !=mu.end(); ++it){
      tempMu.V.muons.push_back(*it);
   }
  }
  if (ele.size()){
    for (std::vector<VHbbEvent::ElectronInfo>::iterator it = ele.begin(); it !=ele.end(); ++it){
      tempE.V.electrons.push_back(*it);
   }
  }

  if (tempMu.V.muons.size()){
   candidates->push_back(tempMu);
  }
  if (tempE.V.electrons.size()){
   candidates->push_back(tempE);
  }
}

bool HbbCandidateFinder::findDiJets (const std::vector<VHbbEvent::SimpleJet>& jets, VHbbEvent::SimpleJet& j1, VHbbEvent::SimpleJet& j2,std::vector<VHbbEvent::SimpleJet>& addJets){
  
 std::vector<VHbbEvent::SimpleJet> tempJets;
 
 for (unsigned int i=0 ; i< jets.size(); ++i){
   if (jets[i].fourMomentum.Pt()> jetPtThreshold)
     tempJets.push_back(jets[i]);
 }
 
 CompareBTag  bTagComparator;

 if (tempJets.size()<2) return false;
 
 std::sort(tempJets.begin(), tempJets.end(), bTagComparator);
 
 j1 = tempJets[0];
 j2 = tempJets[1];
 //
 // additional jets
 //
 if (tempJets.size()>2){
   for (unsigned int i=2 ; i< tempJets.size(); ++i){
     addJets.push_back(tempJets[i]);
   }
 }
  CompareJetPt ptComparator;

  std::sort(addJets.begin(), addJets.end(), ptComparator);
 return true;
}

void HbbCandidateFinder::findMuons(const std::vector<VHbbEvent::MuonInfo>& muons, std::vector<VHbbEvent::MuonInfo>& out){
  /* Use:
For both W -> mu nu and Z -> mu mu, we adopt the standard VBTF muon selection described in VbtfWmunuBaselineSelection. The explicit cuts are reproduced here:

    We use RECO Muons that are both Global and Tracker
    chi2/ndof < 10 for the global muon fit
    The track associated to the muon must have
        >= 1 pixel hits
        >= 10 pixel + strip hits
        >= 1 valid hit in the muon chambers
        >= 2 muon stations
        |dxy| < 0.2
        |eta| < 2.4 
    Relative combined isolation (R) is required to be < 0.15
        R = [Sum pT(trks) + Et(ECAL) + Et(HCAL)] / pT(mu) computed in a cone of radius 0.3 in eta-phi 
    pT(mu) > 20 GeV 
  */
  for (std::vector<VHbbEvent::MuonInfo>::const_iterator it = muons.begin(); it!= muons.end(); ++it){
    if (
        (*it).nPixelHits>= 1 &&
        (*it).nHits +(*it).nPixelHits >= 10 &&
        //
        // sta muon not saved????
        (*it).ipDb<.2 &&
        ((*it).hIso+(*it).eIso+(*it).tIso)/(*it).fourMomentum.Pt()<.15 &&
        (*it).fourMomentum.Eta()<2.4 &&
        (*it).fourMomentum.Pt()>20) {
      out.push_back(*it);
    }
  }
}


void HbbCandidateFinder::findElectrons(const std::vector<VHbbEvent::ElectronInfo>& electrons, std::vector<VHbbEvent::ElectronInfo>& out){
  /*
We adopt the standard cut-based selection from VBTF described in detail here.

    Z -> ee
        gsf electrons
        VBTF WP95
        |eta|<2.5, excluding the gap 1.44 < |eta| < 1.57
        pT(e) > 20 

    W -> e nu
        gsf electrons
        VBTF WP80
        |eta|<2.5, excluding the gap 1.44 < |eta| < 1.57
        pT(e) > 30 
  */

  for (std::vector<VHbbEvent::ElectronInfo>::const_iterator it = electrons.begin(); it!= electrons.end(); ++it){
    if (
        // fake
        //      (*it).id95>  &&
        std::abs((*it).fourMomentum.Eta()) < 2.5 &&
        !(      std::abs((*it).fourMomentum.Eta()) < 1.57 &&    std::abs((*it).fourMomentum.Eta()) > 1.44) &&
        (*it).fourMomentum.Pt()>30
        ){
      out.push_back(*it);
    }  
  }
}


//define this as a plug-in
DEFINE_FWK_MODULE(HbbCandidateFinder);
Revision:	1.3
Committed:	Thu Jun 16 12:16:02 2011 UTC (13 years, 11 months ago) by tboccali
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Jun16th2011
Changes since 1.2:	+146 -66 lines
Log Message:	candidate initial selection
#	User	Rev	Content
1	tboccali	1.1	#include "VHbbAnalysis/HbbAnalyzer/interface/HbbCandidateFinder.h"
2
3	tboccali	1.3	struct CompareJetPt {
4			bool operator()( const VHbbEvent::SimpleJet& j1, const VHbbEvent::SimpleJet& j2 ) const {
5			return j1.fourMomentum.Pt() > j2.fourMomentum.Pt();
6			}
7			};
8	tboccali	1.1
9	tboccali	1.3	struct CompareBTag {
10			bool operator()(const VHbbEvent::SimpleJet& j1, const VHbbEvent::SimpleJet& j2 ) const {
11			return j1.csv > j2.csv;
12			}
13			};
14	tboccali	1.1
15	tboccali	1.3	HbbCandidateFinder::HbbCandidateFinder(const edm::ParameterSet& iConfig): vhbbevent_(iConfig.getParameter<edm::InputTag>("VHbbEventLabel")), verbose_ (iConfig.getParameter<bool>("verbose")), jetPtThreshold(iConfig.getParameter<double>("jetPtThreshold")){
16	tboccali	1.1	produces<std::vector<VHbbCandidate > >();
17			}
18
19			HbbCandidateFinder::~HbbCandidateFinder(){}
20
21			void HbbCandidateFinder::beginJob(){}
22			void HbbCandidateFinder::endJob(){}
23
24
25			float HbbCandidateFinder::getDeltaTheta( VHbbEvent::SimpleJet * j1, VHbbEvent::SimpleJet * j2 ){return -1.;}
26
27
28
29			void HbbCandidateFinder::produce( edm::Event& iEvent, const edm::EventSetup& iEventSetup){
30
31			std::auto_ptr<std::vector<VHbbCandidate> > vHbbCandidates( new std::vector<VHbbCandidate> );
32
33			edm::Handle<VHbbEvent> vHbbEvent;
34	tboccali	1.3	// iEvent.getByLabel(vhbbevent_, vHbbEvent);
35			iEvent.getByType(vHbbEvent);
36	tboccali	1.1
37
38			//
39			// start searching for candidates
40			//
41
42			// hbbCandidateFinderAlgo(vHbbCandidates, vHbbEvent-> result());
43			// do nothing for a test
44
45	tboccali	1.3	run(vHbbEvent.product(),vHbbCandidates);
46
47	tboccali	1.2
48	tboccali	1.3	if (verbose_)
49			std::cout <<" Pushing VHbb candidates: "<<vHbbCandidates->size()<<std::endl;
50	tboccali	1.2
51	tboccali	1.1	iEvent.put(vHbbCandidates);
52
53			}
54
55	tboccali	1.2	void HbbCandidateFinder::run (const VHbbEvent* event, std::auto_ptr<std::vector<VHbbCandidate> > & candidates){
56			//
57			// first find the jets
58			//
59	tboccali	1.3
60			VHbbEvent::SimpleJet j1,j2;
61			std::vector<VHbbEvent::SimpleJet> addJets;
62			bool foundJets = findDiJets(event->simpleJets2,j1,j2,addJets) ;
63
64			if (verbose_){
65			std::cout <<" Found Dijets: "<<foundJets<< " Additional: "<<addJets.size()<< std::endl;
66			}
67	tboccali	1.2
68	tboccali	1.3	if (foundJets == false) return;
69	tboccali	1.2	//
70			// search for a dilepton - just
71			//
72	tboccali	1.3	std::vector<VHbbEvent::MuonInfo> mu;
73			findMuons(event->muInfo,mu);
74			std::vector<VHbbEvent::ElectronInfo> ele;
75			findElectrons(event->eleInfo,ele);
76
77			if (verbose_){
78			std::cout <<" Electrons: "<< ele.size()<<std::endl;
79			std::cout <<" Muons : "<< mu.size()<<std::endl;
80			}
81			if (ele.size()<1 && mu.size() < 1 ) return;
82	tboccali	1.2
83			//
84			// fill!
85			//
86	tboccali	1.3	VHbbCandidate tempMu;
87			VHbbCandidate tempE;
88			tempMu.H.jets.push_back(j1);
89			tempMu.H.jets.push_back(j2);
90			tempMu.H.fourMomentum = (j1).fourMomentum+(j2).fourMomentum;
91			tempMu.additionalJets = addJets;
92			tempE = tempMu;
93			TLorentzVector pMu,pE;
94
95			if (mu.size()){
96			for (std::vector<VHbbEvent::MuonInfo>::iterator it = mu.begin(); it !=mu.end(); ++it){
97			tempMu.V.muons.push_back(*it);
98			}
99			}
100			if (ele.size()){
101			for (std::vector<VHbbEvent::ElectronInfo>::iterator it = ele.begin(); it !=ele.end(); ++it){
102			tempE.V.electrons.push_back(*it);
103			}
104	tboccali	1.2	}
105
106	tboccali	1.3	if (tempMu.V.muons.size()){
107			candidates->push_back(tempMu);
108			}
109			if (tempE.V.electrons.size()){
110			candidates->push_back(tempE);
111			}
112	tboccali	1.2	}
113
114	tboccali	1.3	bool HbbCandidateFinder::findDiJets (const std::vector<VHbbEvent::SimpleJet>& jets, VHbbEvent::SimpleJet& j1, VHbbEvent::SimpleJet& j2,std::vector<VHbbEvent::SimpleJet>& addJets){
115	tboccali	1.2
116	tboccali	1.3	std::vector<VHbbEvent::SimpleJet> tempJets;
117
118			for (unsigned int i=0 ; i< jets.size(); ++i){
119			if (jets[i].fourMomentum.Pt()> jetPtThreshold)
120			tempJets.push_back(jets[i]);
121			}
122
123			CompareBTag bTagComparator;
124
125			if (tempJets.size()<2) return false;
126
127			std::sort(tempJets.begin(), tempJets.end(), bTagComparator);
128
129			j1 = tempJets[0];
130			j2 = tempJets[1];
131			//
132			// additional jets
133			//
134			if (tempJets.size()>2){
135			for (unsigned int i=2 ; i< tempJets.size(); ++i){
136			addJets.push_back(tempJets[i]);
137			}
138			}
139			CompareJetPt ptComparator;
140
141			std::sort(addJets.begin(), addJets.end(), ptComparator);
142			return true;
143			}
144
145			void HbbCandidateFinder::findMuons(const std::vector<VHbbEvent::MuonInfo>& muons, std::vector<VHbbEvent::MuonInfo>& out){
146			/* Use:
147			For both W -> mu nu and Z -> mu mu, we adopt the standard VBTF muon selection described in VbtfWmunuBaselineSelection. The explicit cuts are reproduced here:
148
149			We use RECO Muons that are both Global and Tracker
150			chi2/ndof < 10 for the global muon fit
151			The track associated to the muon must have
152			>= 1 pixel hits
153			>= 10 pixel + strip hits
154			>= 1 valid hit in the muon chambers
155			>= 2 muon stations
156			\|dxy\| < 0.2
157			\|eta\| < 2.4
158			Relative combined isolation (R) is required to be < 0.15
159			R = [Sum pT(trks) + Et(ECAL) + Et(HCAL)] / pT(mu) computed in a cone of radius 0.3 in eta-phi
160			pT(mu) > 20 GeV
161			*/
162			for (std::vector<VHbbEvent::MuonInfo>::const_iterator it = muons.begin(); it!= muons.end(); ++it){
163			if (
164			(*it).nPixelHits>= 1 &&
165			(it).nHits +(it).nPixelHits >= 10 &&
166			//
167			// sta muon not saved????
168			(*it).ipDb<.2 &&
169			((it).hIso+(it).eIso+(it).tIso)/(it).fourMomentum.Pt()<.15 &&
170			(*it).fourMomentum.Eta()<2.4 &&
171			(*it).fourMomentum.Pt()>20) {
172			out.push_back(*it);
173	tboccali	1.2	}
174			}
175			}
176
177
178	tboccali	1.3	void HbbCandidateFinder::findElectrons(const std::vector<VHbbEvent::ElectronInfo>& electrons, std::vector<VHbbEvent::ElectronInfo>& out){
179			/*
180			We adopt the standard cut-based selection from VBTF described in detail here.
181
182			Z -> ee
183			gsf electrons
184			VBTF WP95
185			\|eta\|<2.5, excluding the gap 1.44 < \|eta\| < 1.57
186			pT(e) > 20
187
188			W -> e nu
189			gsf electrons
190			VBTF WP80
191			\|eta\|<2.5, excluding the gap 1.44 < \|eta\| < 1.57
192			pT(e) > 30
193			*/
194
195			for (std::vector<VHbbEvent::ElectronInfo>::const_iterator it = electrons.begin(); it!= electrons.end(); ++it){
196			if (
197			// fake
198			// (*it).id95> &&
199			std::abs((*it).fourMomentum.Eta()) < 2.5 &&
200			!( std::abs((it).fourMomentum.Eta()) < 1.57 && std::abs((it).fourMomentum.Eta()) > 1.44) &&
201			(*it).fourMomentum.Pt()>30
202			){
203			out.push_back(*it);
204			}
205			}
206	tboccali	1.2	}
207
208	tboccali	1.1
209			//define this as a plug-in
210			DEFINE_FWK_MODULE(HbbCandidateFinder);