ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/HbbAnalysis/plugins/HbbTreeMaker.cc
Revision: 1.24
Committed: Mon Oct 25 11:03:45 2010 UTC (14 years, 6 months ago) by amagnan
Content type: text/plain
Branch: MAIN
Changes since 1.23: +17 -0 lines
Log Message:
add vertex rho and ak7 jet collection to tree

File Contents

# Content
1 #include "DataFormats/Common/interface/Handle.h"
2 #include "DataFormats/Common/interface/TriggerResults.h"
3 #include "DataFormats/Common/interface/HLTenums.h"
4 #include "DataFormats/Common/interface/ValueMap.h"
5
6 #include "DataFormats/Candidate/interface/Candidate.h"
7 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
8 #include "DataFormats/Scalers/interface/DcsStatus.h"
9 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
10 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
11
12 #include "DataFormats/MuonReco/interface/Muon.h"
13 #include "DataFormats/MuonReco/interface/MuonFwd.h"
14 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
15 #include "DataFormats/MuonReco/interface/MuonChamberMatch.h"
16 #include "DataFormats/MuonReco/interface/MuonIsolation.h"
17 #include "DataFormats/MuonReco/interface/MuonEnergy.h"
18 #include "DataFormats/MuonReco/interface/MuonTime.h"
19 #include "DataFormats/MuonReco/interface/MuonQuality.h"
20
21 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
22 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
23 #include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
24
25
26 #include "DataFormats/VertexReco/interface/Vertex.h"
27 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
28 #include "DataFormats/TrackReco/interface/Track.h"
29 #include "DataFormats/TrackReco/interface/TrackFwd.h"
30
31 #include "DataFormats/PatCandidates/interface/Lepton.h"
32 #include "DataFormats/PatCandidates/interface/Muon.h"
33 #include "DataFormats/PatCandidates/interface/Electron.h"
34 #include "DataFormats/PatCandidates/interface/Tau.h"
35 #include "DataFormats/PatCandidates/interface/Jet.h"
36 #include "DataFormats/PatCandidates/interface/MET.h"
37
38 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
39 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
40
41 #include "FWCore/ServiceRegistry/interface/Service.h"
42 #include "FWCore/Framework/interface/Run.h"
43 #include "FWCore/Framework/interface/ESHandle.h"
44 #include "CommonTools/UtilAlgos/interface/TFileService.h"
45
46 #include "PhysicsTools/SelectorUtils/interface/SimpleCutBasedElectronIDSelectionFunctor.h"
47 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
48 #include "MagneticField/Engine/interface/MagneticField.h"
49
50 #include "UserCode/HbbAnalysis/interface/Electron.hh"
51 #include "UserCode/HbbAnalysis/interface/Muon.hh"
52 #include "UserCode/HbbAnalysis/interface/Tau.hh"
53 #include "UserCode/HbbAnalysis/interface/Jet.hh"
54 #include "UserCode/HbbAnalysis/interface/Met.hh"
55 #include "UserCode/HbbAnalysis/interface/Trigger.hh"
56 #include "UserCode/HbbAnalysis/interface/Vertex.hh"
57 #include "UserCode/HbbAnalysis/interface/L1Object.hh"
58 #include "UserCode/HbbAnalysis/interface/HLTObject.hh"
59
60 #include "UserCode/HbbAnalysis/plugins/HbbTreeMaker.hh"
61
62 using namespace HbbAnalysis;
63
64 HbbTreeMaker::HbbTreeMaker(const edm::ParameterSet & pset):
65 debug_(pset.getParameter<int>("DEBUG")),
66 processData_(false),
67 flavour_(pset.getParameter<unsigned int>("JetFlavour")),
68 doGen_(pset.getParameter<bool>("DOGEN")),
69 genParticleSrc_(pset.getParameter<edm::InputTag>("GenParticles")),
70 trackSrc_(pset.getParameter<edm::InputTag>("Tracks")),
71 dcsSrc_(pset.getParameter<edm::InputTag>("DCSInfo")),
72 electronSrc_(pset.getParameter<edm::InputTag>("Electrons")),
73 muonSrc_(pset.getParameter<edm::InputTag>("Muons")),
74 caloTauSrc_(pset.getParameter<edm::InputTag>("CaloTaus")),
75 pfTauSrc_(pset.getParameter<edm::InputTag>("PFTaus")),
76 caloJetSrc_(pset.getParameter<edm::InputTag>("CaloJets")),
77 jptJetSrc_(pset.getParameter<edm::InputTag>("JPTJets")),
78 ak7JetSrc_(pset.getParameter<edm::InputTag>("AK7Jets")),
79 pfJetSrc_(pset.getParameter<edm::InputTag>("PFJets")),
80 caloMetSrc_(pset.getParameter<edm::InputTag>("CaloMET")),
81 tcMetSrc_(pset.getParameter<edm::InputTag>("TCMET")),
82 pfMetSrc_(pset.getParameter<edm::InputTag>("PFMET")),
83 //pairSrc_(pset.getParameter<edm::InputTag>("Pair")),
84 //mmPairSrc_(pset.getParameter<edm::InputTag>("MuMuPair")),
85 //etPairSrc_(pset.getParameter<edm::InputTag>("ETauPair")),
86 //mtPairSrc_(pset.getParameter<edm::InputTag>("MuTauPair")),
87 vertexSrc_(pset.getParameter<edm::InputTag>("Vertex")),
88 triggerSrc_(pset.getParameter<edm::InputTag>("Trigger")),
89 hltPaths_(pset.getParameter<std::vector<std::string> >("HLTPaths")),
90 l1CenJetSrc_(pset.getParameter<edm::InputTag>("L1CenJet")),
91 l1TauJetSrc_(pset.getParameter<edm::InputTag>("L1TauJet")),
92 l1FwdJetSrc_(pset.getParameter<edm::InputTag>("L1FwdJet")),
93 hltSummarySrc_(pset.getParameter<edm::InputTag>("HLTSummary")),
94 event_(0),
95 tree_(0),
96 hltTagsElec_(pset.getParameter<std::vector<edm::InputTag> >("HLTTagsElectrons")),
97 hltTagsMu_(pset.getParameter<std::vector<edm::InputTag> >("HLTTagsMuons")),
98 hltTagsJet_(pset.getParameter<std::vector<edm::InputTag> >("HLTTagsJets"))
99 //processVec_(pset.getParameter<std::vector<std::string> >("ProcessVec"))
100 {//constructor
101
102
103 }//constructor
104
105 HbbTreeMaker::~HbbTreeMaker()
106 {//destructor
107 }//destructor
108
109
110
111 void HbbTreeMaker::beginJob(){//beginJob
112
113
114 edm::Service<TFileService> lFileService;
115
116 TFileDirectory lDir = lFileService->mkdir("HbbAnalysis");
117
118 tree_ = lDir.make<TTree>("Tree","Tree");
119 tree_->Branch("HbbEvent","HbbAnalysis::HbbEvent",&event_);
120 event_ = new HbbEvent();
121
122 if (debug_) std::cout << "Initialising JetFlavour : " << std::endl;
123
124 lDir = lFileService->mkdir("JetFlavours");
125 jetFlav_.Initialise(lDir, debug_, flavour_);
126
127
128 }//beginJob
129
130 void HbbTreeMaker::endJob(){//endJob
131 if (!processData_) jetFlav_.printSummary();
132
133 //tree_->Write();
134 //delete tree_;
135 //delete event_;
136
137 }//endJob
138
139 void HbbTreeMaker::analyze(const edm::Event& aEvt, const edm::EventSetup& aEvtSetup){//analyze
140
141 //dumpContent(aEvt,"","MuTauPAT");
142 if (debug_) std::cout << "Processing event #" << aEvt.id().event() << std::endl;
143
144
145 static bool firstEvent = true;
146
147 event_->Clear();
148 event_->event(aEvt.id().event());
149 processData_ = aEvt.isRealData();
150
151 event_->run(aEvt.run());
152 event_->lumiBlock(aEvt.luminosityBlock());
153 event_->isRealData(aEvt.isRealData());
154 event_->bx(aEvt.bunchCrossing());
155
156 edm::Handle<reco::BeamSpot> lBeamSpot;
157 try {
158 aEvt.getByLabel(edm::InputTag("offlineBeamSpot"), lBeamSpot);
159 } catch(cms::Exception& e) {
160 std::cout<< "AMM: No beam spot found. Exception : " << e.what() << "." << std::endl;
161 }
162 HbbBeamSpot(lBeamSpot,event_->beamSpot());
163
164
165 if (!processData_) {
166 edm::Handle<edm::HepMCProduct> mcProd;
167 edm::Handle<GenRunInfoProduct> lRunProd;
168 try {
169 aEvt.getByLabel("generator", mcProd );
170 } catch(cms::Exception& e) {
171 std::cout << "AMM: Collection HepMCProduct not available! Exception : " << e.what() << ". " << std::endl;
172 }
173 try {
174 aEvt.getRun().getByLabel("generator",lRunProd);
175 } catch(cms::Exception& e) {
176 std::cout << "AMM: Collection GenInfoProduct not available! Exception : " << e.what() << ". " << std::endl;
177 }
178
179 HbbGenInfo(mcProd,lRunProd);
180 }
181
182 edm::Handle<reco::GenParticleCollection> lGenParticles;
183 try {
184 aEvt.getByLabel(genParticleSrc_,lGenParticles);
185 if (debug_) std::cout << "** ngenParticles = " << lGenParticles->size() << std::endl;
186 } catch(cms::Exception& e) {
187 if (!processData_) std::cout << "AMM: Collection genParticles not available! Exception : " << e.what() << ". " << std::endl;
188 }
189
190 if (doGen_ && !processData_) HbbParticles(lGenParticles,event_->particles());
191
192 unsigned int lNPartons = 0;
193 if (!processData_) lNPartons = jetFlav_.fillPartons(lGenParticles);
194
195 if (debug_ && !processData_) std::cout << "--- Number of partons = " << lNPartons << "." << std::endl;
196
197
198 edm::Handle<std::vector<reco::Vertex> > lRecoVertices;
199 try {
200 aEvt.getByLabel(vertexSrc_, lRecoVertices);
201 if (!lRecoVertices.isValid()){
202 edm::LogInfo("ERROR")<< "Error! Can't get vertex by label. ";
203 }
204 if (debug_) std::cout << "** vertexcollection = " << lRecoVertices->size() << " elements." << std::endl;
205 } catch(cms::Exception& e) {
206 std::cout << "AMM: Collection " << vertexSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
207 }
208
209 HbbVertices(lRecoVertices,event_->vertices());
210
211 //to select elecID
212 edm::Handle<reco::TrackCollection> lCtfTracks;
213 try {
214 aEvt.getByLabel(trackSrc_,lCtfTracks);
215 if (!lCtfTracks.isValid()){
216 edm::LogInfo("ERROR")<< "Error! Can't get tracks by label. ";
217 }
218 if (debug_) std::cout << "** trackcollection = " << lCtfTracks->size() << " elements." << std::endl;
219 } catch(cms::Exception& e) {
220 std::cout << "AMM: Collection " << trackSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
221 }
222
223 // need the magnetic field
224 double lBfield = 0;
225 // if isData then derive bfield using the
226 // magnet current from DcsStatus
227 // otherwise take it from the IdealMagneticFieldRecord
228 if (processData_) {
229 // scale factor = 3.801/18166.0 which are
230 // average values taken over a stable two
231 // week period
232 edm::Handle<DcsStatusCollection> dcsHandle;
233 try {
234 aEvt.getByLabel(dcsSrc_, dcsHandle);
235 if (!dcsHandle.isValid()){
236 edm::LogInfo("ERROR")<< "Error! Can't get tracks by label. ";
237 }
238 if (debug_) std::cout << "** DCScollection = " << dcsHandle->size() << " elements." << std::endl;
239 } catch(cms::Exception& e) {
240 std::cout << "AMM: Collection " << dcsSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
241 }
242 if (dcsHandle.isValid()){
243 float currentToBFieldScaleFactor = 2.09237036221512717e-04;
244 float current = (*dcsHandle)[0].magnetCurrent();
245 lBfield = current*currentToBFieldScaleFactor;
246 }
247 }
248 else {
249 edm::ESHandle<MagneticField> magneticField;
250 aEvtSetup.get<IdealMagneticFieldRecord>().get(magneticField);
251 lBfield = magneticField->inTesla(GlobalPoint(0.,0.,0.)).z();
252 }
253
254 edm::Handle<std::vector<pat::Electron> > lElectronCollection;
255
256 try {
257 aEvt.getByLabel(electronSrc_,lElectronCollection);
258 if (!lElectronCollection.isValid()){
259 edm::LogInfo("ERROR")<< "Error! Can't get electron by label. ";
260 }
261 if (debug_) std::cout << "** electroncollection = " << lElectronCollection->size() << " elements." << std::endl;
262 } catch(cms::Exception& e) {
263 std::cout << "AMM: Collection " << electronSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
264 }
265
266 HbbElectrons(lCtfTracks,lBfield,lElectronCollection,event_->electrons());
267
268 edm::Handle<std::vector<pat::Muon> > lMuonCollection;
269
270 try {
271 aEvt.getByLabel(muonSrc_,lMuonCollection);
272 if (!lMuonCollection.isValid()){
273 edm::LogInfo("ERROR")<< "Error! Can't get muon by label. ";
274 }
275 if (debug_) std::cout << "** muoncollection = " << lMuonCollection->size() << " elements." << std::endl;
276 } catch(cms::Exception& e) {
277 std::cout << "AMM: Collection " << muonSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
278 }
279
280 HbbMuons(lMuonCollection,lRecoVertices,lBeamSpot,event_->muons());
281
282
283 if (!( caloTauSrc_.label()=="" && caloTauSrc_.instance()=="" )) {
284 edm::Handle<std::vector<pat::Tau> > lTauCollection;
285
286 try {
287 aEvt.getByLabel(caloTauSrc_,lTauCollection);
288 if (!lTauCollection.isValid()){
289 edm::LogInfo("ERROR")<< "Error! Can't get caloTau by label. ";
290 }
291 if (debug_) std::cout << "** caloTaucollection = " << lTauCollection->size() << " elements." << std::endl;
292 } catch(cms::Exception& e) {
293 std::cout << "AMM: Collection " << caloTauSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
294 }
295
296 HbbTaus(lTauCollection,lRecoVertices,event_->caloTaus());
297 }
298
299 edm::Handle<std::vector<pat::Tau> > lPFTauCollection;
300
301 try {
302 aEvt.getByLabel(pfTauSrc_,lPFTauCollection);
303 if (!lPFTauCollection.isValid()){
304 edm::LogInfo("ERROR")<< "Error! Can't get pfTau by label. ";
305 }
306 if (debug_) std::cout << "** pfTaucollection = " << lPFTauCollection->size() << " elements." << std::endl;
307 } catch(cms::Exception& e) {
308 std::cout << "AMM: Collection " << pfTauSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
309 }
310
311 HbbTaus(lPFTauCollection,lRecoVertices,event_->pfTaus());
312
313 edm::Handle<std::vector<pat::Jet> > lCaloJetCollection;
314
315 try {
316 aEvt.getByLabel(caloJetSrc_,lCaloJetCollection);
317 if (!lCaloJetCollection.isValid()){
318 edm::LogInfo("ERROR")<< "Error! Can't get caloJets by label. ";
319 }
320 if (debug_) std::cout << "** caloJetcollection = " << lCaloJetCollection->size() << " elements." << std::endl;
321 } catch(cms::Exception& e) {
322 std::cout << "AMM: Collection " << caloJetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
323 }
324
325 //std::cout << "Processing calo jets:" << std::endl;
326 HbbJets(lCaloJetCollection,jetFlav_,lGenParticles,event_->caloJets());
327
328 edm::Handle<std::vector<pat::Jet> > lJptJetCollection;
329
330 try {
331 aEvt.getByLabel(jptJetSrc_,lJptJetCollection);
332 if (!lJptJetCollection.isValid()){
333 edm::LogInfo("ERROR")<< "Error! Can't get jptJets by label. ";
334 }
335 if (debug_) std::cout << "** jptJetcollection = " << lJptJetCollection->size() << " elements." << std::endl;
336 } catch(cms::Exception& e) {
337 std::cout << "AMM: Collection " << jptJetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
338 }
339
340 //std::cout << "Processing JPT jets:" << std::endl;
341 HbbJets(lJptJetCollection,jetFlav_,lGenParticles,event_->jptJets());
342
343 edm::Handle<std::vector<pat::Jet> > lAk7JetCollection;
344
345 try {
346 aEvt.getByLabel(ak7JetSrc_,lAk7JetCollection);
347 if (!lAkJetCollection.isValid()){
348 edm::LogInfo("ERROR")<< "Error! Can't get ak7Jets by label. ";
349 }
350 if (debug_) std::cout << "** ak7Jetcollection = " << lAk7JetCollection->size() << " elements." << std::endl;
351 } catch(cms::Exception& e) {
352 std::cout << "AMM: Collection " << ak7JetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
353 }
354
355 //std::cout << "Processing JPT jets:" << std::endl;
356 HbbJets(lAk7JetCollection,jetFlav_,lGenParticles,event_->ak7Jets());
357
358 edm::Handle<std::vector<pat::Jet> > lPfJetCollection;
359
360 try {
361 aEvt.getByLabel(pfJetSrc_,lPfJetCollection);
362 if (!lPfJetCollection.isValid()){
363 edm::LogInfo("ERROR")<< "Error! Can't get pfJets by label. ";
364 }
365 if (debug_) std::cout << "** pfJetcollection = " << lPfJetCollection->size() << " elements." << std::endl;
366 } catch(cms::Exception& e) {
367 std::cout << "AMM: Collection " << pfJetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
368 }
369
370 //std::cout << "Processing PF jets:" << std::endl;
371 HbbJets(lPfJetCollection,jetFlav_,lGenParticles,event_->pfJets());
372
373 edm::Handle<std::vector<pat::MET> > lCaloMetCol;
374
375 try {
376 aEvt.getByLabel(caloMetSrc_,lCaloMetCol);
377 if (!lCaloMetCol.isValid()){
378 edm::LogInfo("ERROR")<< "Error! Can't get caloMet by label. ";
379 }
380 if (debug_) std::cout << "** caloMet = " << lCaloMetCol->size() << " elements." << std::endl;
381 } catch(cms::Exception& e) {
382 std::cout << "AMM: Collection " << caloMetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
383 }
384
385 HbbMet(lCaloMetCol,event_->caloMet());
386
387 edm::Handle<std::vector<pat::MET> > lTcMetCol;
388
389 try {
390 aEvt.getByLabel(tcMetSrc_,lTcMetCol);
391 if (!lTcMetCol.isValid()){
392 edm::LogInfo("ERROR")<< "Error! Can't get tcMet by label. ";
393 }
394 if (debug_) std::cout << "** tcMet = " << lTcMetCol->size() << " elements." << std::endl;
395 } catch(cms::Exception& e) {
396 std::cout << "AMM: Collection " << tcMetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
397 }
398
399 HbbMet(lTcMetCol,event_->tcMet());
400
401 edm::Handle<std::vector<pat::MET> > lPfMetCol;
402
403 try {
404 aEvt.getByLabel(pfMetSrc_,lPfMetCol);
405 if (!lPfMetCol.isValid()){
406 edm::LogInfo("ERROR")<< "Error! Can't get pfMet by label. ";
407 }
408 if (debug_) std::cout << "** pfMet = " << lPfMetCol->size() << " elements." << std::endl;
409 } catch(cms::Exception& e) {
410 std::cout << "AMM: Collection " << pfMetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
411 }
412
413 HbbMet(lPfMetCol,event_->pfMet());
414
415
416
417 //triggers
418 edm::Handle<edm::TriggerResults> lTrigCol;
419 try {
420 aEvt.getByLabel(triggerSrc_,lTrigCol);
421 if (!lTrigCol.isValid()){
422 edm::LogInfo("ERROR")<< "Error! Can't get triggers by label. ";
423 }
424 if (debug_) std::cout << "** triggers = " << lTrigCol->size() << std::endl;
425 } catch(cms::Exception& e) {
426 std::cout << "AMM: Collection " << triggerSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
427 }
428 if (lTrigCol.isValid()) {
429 const edm::TriggerNames & lNames = aEvt.triggerNames(*lTrigCol);
430 HbbTrigger(lTrigCol,lNames,event_->triggers());
431 }
432
433 //L1 jet objects
434
435 edm::Handle<l1extra::L1JetParticleCollection> lCenJet;
436 edm::Handle<l1extra::L1JetParticleCollection> lFwdJet;
437 edm::Handle<l1extra::L1JetParticleCollection> lTauJet;
438 try {
439 aEvt.getByLabel(l1CenJetSrc_,lCenJet);
440 aEvt.getByLabel(l1TauJetSrc_,lTauJet);
441 aEvt.getByLabel(l1FwdJetSrc_,lFwdJet);
442 } catch(cms::Exception& e) {
443 std::cout << "AMM: L1 Collections " << l1CenJetSrc_ << ", " << l1TauJetSrc_ << ", " << l1FwdJetSrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
444 }
445 if (lCenJet.isValid()) {
446 HbbL1Objects(lCenJet,event_->l1Jets(),1);
447 }
448 if (lTauJet.isValid()) {
449 HbbL1Objects(lTauJet,event_->l1Jets(),2);
450 }
451 if (lFwdJet.isValid()) {
452 HbbL1Objects(lFwdJet,event_->l1Jets(),3);
453 }
454
455
456 edm::Handle<trigger::TriggerEvent> lHltSummary;
457 try {
458 aEvt.getByLabel(hltSummarySrc_,lHltSummary);
459 } catch(cms::Exception& e) {
460 std::cout << "AMM: collection " << hltSummarySrc_ << " not available! Exception : " << e.what() << ". " << std::endl;
461 }
462 if (lHltSummary.isValid()) {
463 HbbHLTObjects(lHltSummary,event_->hltObjects());
464 }
465
466
467
468 //event_->order();
469 tree_->Fill();
470 firstEvent = false;
471
472 }//analyze
473
474
475 void HbbTreeMaker::HbbBeamSpot(const edm::Handle<reco::BeamSpot> & aBeamSpot,
476 HbbAnalysis::BeamSpot & aBS)
477 {
478
479 if (aBeamSpot.isValid()) {
480 HbbAnalysis::BeamSpotVars lbs;
481 lbs.x0 = aBeamSpot->x0();
482 lbs.y0 = aBeamSpot->y0();
483 lbs.z0 = aBeamSpot->z0();
484 lbs.sigmaZ = aBeamSpot->sigmaZ();
485 lbs.BeamWidthX = aBeamSpot->BeamWidthX();
486 lbs.BeamWidthY = aBeamSpot->BeamWidthY();
487 lbs.x0Error = aBeamSpot->x0Error();
488 lbs.y0Error = aBeamSpot->y0Error();
489 lbs.z0Error = aBeamSpot->z0Error();
490 lbs.sigmaZ0Error = aBeamSpot->sigmaZ0Error();
491 lbs.BeamWidthXError = aBeamSpot->BeamWidthXError();
492 lbs.BeamWidthYError = aBeamSpot->BeamWidthYError();
493
494 aBS.bsVars(lbs);
495 }
496
497 }
498
499
500 void HbbTreeMaker::HbbElectrons(const edm::Handle<reco::TrackCollection> & aCtfTracks,
501 const double aBfield,
502 const edm::Handle<std::vector<pat::Electron> > & aCol,
503 std::vector<HbbAnalysis::Electron> & aVec)
504 {//HbbElectrons
505
506 if (aCol.isValid()){
507 if (aCol->size() > 0) {
508
509 unsigned int iEle = 0;
510 for (std::vector<pat::Electron>::const_iterator iter = aCol->begin();
511 iter != aCol->end();
512 ++iter,++iEle)
513 {
514 //if (debug_ > 1) std::cout << "**** Ele #" << iEle << ", pT,eta,phi = " << (*iter).pt() << " " << (*iter).eta() << " " << (*iter).phi() << std::endl;
515 // get the track collection
516
517
518
519 // define your function for WP95 using relative isolations
520 SimpleCutBasedElectronIDSelectionFunctor patSele95(SimpleCutBasedElectronIDSelectionFunctor::relIso95, aBfield, aCtfTracks);
521
522 HbbAnalysis::GenVars lGen;
523 if (!processData_ && (*iter).genLepton()){
524 lGen.valid = true;
525 lGen.E = (*iter).genLepton()->energy();
526 lGen.pT = (*iter).genLepton()->pt();
527 lGen.eta = (*iter).genLepton()->eta();
528 lGen.y = (*iter).genLepton()->y();
529 lGen.phi = (*iter).genLepton()->phi();
530 lGen.charge = (*iter).genLepton()->charge();
531 lGen.pdgId = (*iter).genLepton()->pdgId();
532 lGen.status = (*iter).genLepton()->status();
533 lGen.mass = (*iter).genLepton()->mass();
534 lGen.vx = (*iter).genLepton()->vx();
535 lGen.vy = (*iter).genLepton()->vy();
536 lGen.vz = (*iter).genLepton()->vz();
537 }
538 else {
539 lGen.valid = false;
540 lGen.E = 0;
541 lGen.pT = 0;
542 lGen.eta = 0;
543 lGen.y = 0;
544 lGen.phi = 0;
545 lGen.charge = 0;
546 lGen.pdgId = 0;
547 lGen.status = 0;
548 lGen.mass = 0;
549 lGen.vx = 0;
550 lGen.vy = 0;
551 lGen.vz = 0;
552 }
553
554 HbbAnalysis::BaseVars lReco;
555 lReco.E = (*iter).energy();
556 lReco.pT = (*iter).pt();
557 lReco.eta = (*iter).eta();
558 lReco.y = (*iter).y();
559 lReco.phi = (*iter).phi();
560 lReco.charge = (*iter).charge();
561 lReco.vx = (*iter).vx();
562 lReco.vy = (*iter).vy();
563 lReco.vz = (*iter).vz();
564
565 HbbAnalysis::SCVars lSC;
566 lSC.sigmaEtaEta = (*iter).scSigmaEtaEta();
567 lSC.sigmaIEtaIEta = (*iter).scSigmaIEtaIEta();
568 lSC.e1x5 = (*iter).scE1x5();
569 lSC.e2x5Max = (*iter).scE2x5Max();
570 lSC.e5x5 = (*iter).scE5x5();
571 lSC.eOverP = (*iter).eSuperClusterOverP();
572
573 HbbAnalysis::EleIsoVars lIso;
574 lIso.calo = (*iter).caloIso();
575 lIso.track = (*iter).trackIso();
576 lIso.ecal = (*iter).ecalIso();
577 lIso.hcal = (*iter).hcalIso();
578
579 HbbAnalysis::EleIDVars lID;
580 lID.idAndIso = patSele95(*iter);
581 lID.electronIDs = (*iter).electronIDs();
582 lID.hOverE = (*iter).hadronicOverEm();
583 lID.deltaPhiIn = (*iter).deltaPhiSuperClusterTrackAtVtx();
584 lID.deltaEtaIn = (*iter).deltaEtaSuperClusterTrackAtVtx();
585 lID.ecalDrivenSeed = (*iter).ecalDrivenSeed();
586 lID.trackerDrivenSeed = (*iter).trackerDrivenSeed();
587
588 HbbAnalysis::Electron lObj(lGen,lReco,lSC,lIso,lID);
589 aVec.push_back(lObj);
590 }
591 }
592 }
593
594 }//HbbElectrons
595
596
597 void HbbTreeMaker::HbbMuons(const edm::Handle<std::vector<pat::Muon> > & aCol,
598 const edm::Handle<std::vector<reco::Vertex> > & aRecoVertices,
599 const edm::Handle<reco::BeamSpot> & aBeamSpot,
600 std::vector<HbbAnalysis::Muon> & aVec)
601 {//HbbMuons
602
603 if (aCol.isValid()){
604 if (aCol->size() > 0) {
605
606 unsigned int iEle = 0;
607 for (std::vector<pat::Muon>::const_iterator iter = aCol->begin();
608 iter != aCol->end();
609 ++iter,++iEle)
610 {
611 const reco::Muon* recoMuon = dynamic_cast<const reco::Muon*>((*iter).originalObject());
612
613 //select only global and tracker muons
614 if (!recoMuon->isGlobalMuon() || !recoMuon->isTrackerMuon()) continue;
615
616 HbbAnalysis::GenVars lGen;
617 if (!processData_ && (*iter).genLepton()){
618 lGen.valid = true;
619 lGen.E = (*iter).genLepton()->energy();
620 lGen.pT = (*iter).genLepton()->pt();
621 lGen.eta = (*iter).genLepton()->eta();
622 lGen.y = (*iter).genLepton()->y();
623 lGen.phi = (*iter).genLepton()->phi();
624 lGen.charge = (*iter).genLepton()->charge();
625 lGen.pdgId = (*iter).genLepton()->pdgId();
626 lGen.status = (*iter).genLepton()->status();
627 lGen.mass = (*iter).genLepton()->mass();
628 lGen.vx = (*iter).genLepton()->vx();
629 lGen.vy = (*iter).genLepton()->vy();
630 lGen.vz = (*iter).genLepton()->vz();
631 }
632 else {
633 lGen.valid = false;
634 lGen.E = 0;
635 lGen.pT = 0;
636 lGen.eta = 0;
637 lGen.y = 0;
638 lGen.phi = 0;
639 lGen.charge = 0;
640 lGen.pdgId = 0;
641 lGen.status = 0;
642 lGen.mass = 0;
643 lGen.vx = 0;
644 lGen.vy = 0;
645 lGen.vz = 0;
646 }
647
648 HbbAnalysis::BaseVars lReco;
649 lReco.E = (*iter).energy();
650 lReco.pT = (*iter).pt();
651 lReco.eta = (*iter).eta();
652 lReco.y = (*iter).y();
653 lReco.phi = (*iter).phi();
654 lReco.charge = (*iter).charge();
655 lReco.vx = (*iter).vx();
656 lReco.vy = (*iter).vy();
657 lReco.vz = (*iter).vz();
658
659 HbbAnalysis::MuTrkVars lTrk;
660 reco::TrackRef glmuon = recoMuon->globalTrack();
661 reco::TrackRef tkmuon = recoMuon->innerTrack();
662
663 if ( glmuon.isAvailable() && glmuon.isNonnull() ) {
664 if (aBeamSpot.isValid()) {
665 lTrk.dxy = glmuon->dxy(aBeamSpot->position());
666 lTrk.dz = glmuon->dz(aBeamSpot->position());
667 }
668 else {
669 lTrk.dxy = 0;
670 lTrk.dz = 0;
671 }
672 lTrk.normalizedChi2 = glmuon->normalizedChi2();
673 lTrk.muonHits = glmuon->hitPattern().numberOfValidMuonHits();
674 lTrk.charge = glmuon->charge();
675 }
676 else {
677 lTrk.dxy = 0;
678 lTrk.dz = 0;
679 lTrk.normalizedChi2 = 0;
680 lTrk.muonHits = 0;
681 lTrk.charge = 0;
682 }
683 if ( tkmuon.isAvailable() && tkmuon.isNonnull() ) {
684 lTrk.trackerHits = tkmuon->hitPattern().numberOfValidTrackerHits();
685 lTrk.pixelHits = tkmuon->hitPattern().numberOfValidPixelHits();
686 }
687 else {
688 lTrk.trackerHits = 0;
689 lTrk.pixelHits = 0;
690 }
691
692 HbbAnalysis::MuIsoVars lIsoR03;
693 lIsoR03.sumPt = recoMuon->isolationR03().sumPt;
694 lIsoR03.emEt = recoMuon->isolationR03().emEt;
695 lIsoR03.hadEt = recoMuon->isolationR03().hadEt;
696 lIsoR03.nTracks = recoMuon->isolationR03().nTracks;
697 lIsoR03.nJets = recoMuon->isolationR03().nJets;
698
699 HbbAnalysis::MuIsoVars lIsoR05;
700 lIsoR05.sumPt = recoMuon->isolationR05().sumPt;
701 lIsoR05.emEt = recoMuon->isolationR05().emEt;
702 lIsoR05.hadEt = recoMuon->isolationR05().hadEt;
703 lIsoR05.nTracks = recoMuon->isolationR05().nTracks;
704 lIsoR05.nJets = recoMuon->isolationR05().nJets;
705
706 HbbAnalysis::MuIDVars lID;
707 //if (recoMuon->isGlobalMuon()) lID.type = 1;
708 //else if (recoMuon->isTrackerMuon()) lID.type = 2;
709 //else if (recoMuon->isStandAloneMuon()) lID.type = 3;
710 //else if (recoMuon->isCaloMuon()) lID.type = 4;
711 //else lID.type = 0;
712 lID.type = (recoMuon->isCaloMuon()<<3) | (recoMuon->isStandAloneMuon()<<2) | (recoMuon->isTrackerMuon()<<1) | recoMuon->isGlobalMuon();
713
714 if (muon::isGoodMuon(*recoMuon,muon::AllGlobalMuons)) lID.ids.push_back(1);
715 if (muon::isGoodMuon(*recoMuon,muon::AllStandAloneMuons)) lID.ids.push_back(2);
716 if (muon::isGoodMuon(*recoMuon,muon::AllTrackerMuons)) lID.ids.push_back(3);
717 if (muon::isGoodMuon(*recoMuon,muon::TrackerMuonArbitrated)) lID.ids.push_back(4);
718 if (muon::isGoodMuon(*recoMuon,muon::AllArbitrated)) lID.ids.push_back(5);
719 if (muon::isGoodMuon(*recoMuon,muon::GlobalMuonPromptTight)) lID.ids.push_back(6);
720 if (muon::isGoodMuon(*recoMuon,muon::TMLastStationLoose)) lID.ids.push_back(7);
721 if (muon::isGoodMuon(*recoMuon,muon::TMLastStationTight)) lID.ids.push_back(8);
722 if (muon::isGoodMuon(*recoMuon,muon::TM2DCompatibilityLoose)) lID.ids.push_back(9);
723 if (muon::isGoodMuon(*recoMuon,muon::TM2DCompatibilityTight)) lID.ids.push_back(10);
724 if (muon::isGoodMuon(*recoMuon,muon::TMOneStationLoose)) lID.ids.push_back(11);
725 if (muon::isGoodMuon(*recoMuon,muon::TMOneStationTight)) lID.ids.push_back(12);
726 if (muon::isGoodMuon(*recoMuon,muon::TMLastStationOptimizedLowPtLoose)) lID.ids.push_back(13);
727 if (muon::isGoodMuon(*recoMuon,muon::TMLastStationOptimizedLowPtTight)) lID.ids.push_back(14);
728 if (muon::isGoodMuon(*recoMuon,muon::GMTkChiCompatibility)) lID.ids.push_back(15);
729 if (muon::isGoodMuon(*recoMuon,muon::GMStaChiCompatibility)) lID.ids.push_back(16);
730 if (muon::isGoodMuon(*recoMuon,muon::GMTkKinkTight)) lID.ids.push_back(17);
731 if (muon::isGoodMuon(*recoMuon,muon::TMLastStationAngLoose)) lID.ids.push_back(18);
732 if (muon::isGoodMuon(*recoMuon,muon::TMLastStationAngTight)) lID.ids.push_back(19);
733 if (muon::isGoodMuon(*recoMuon,muon::TMOneStationAngLoose)) lID.ids.push_back(20);
734 if (muon::isGoodMuon(*recoMuon,muon::TMOneStationAngTight)) lID.ids.push_back(21);
735 if (muon::isGoodMuon(*recoMuon,muon::TMLastStationOptimizedBarrelLowPtLoose)) lID.ids.push_back(22);
736 if (muon::isGoodMuon(*recoMuon,muon::TMLastStationOptimizedBarrelLowPtTight)) lID.ids.push_back(23);
737
738 lID.caloCompat = recoMuon->caloCompatibility();
739 lID.segCompat = muon::segmentCompatibility(*recoMuon);
740 lID.nChambers = recoMuon->numberOfChambers();
741 lID.nMatchesLoose = recoMuon->numberOfMatches(reco::Muon::NoArbitration);
742 lID.nMatchesMedium = recoMuon->numberOfMatches(reco::Muon::SegmentArbitration);
743 lID.nMatchesTight = recoMuon->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
744
745 HbbAnalysis::Muon lObj(lGen,lReco,lTrk,lIsoR03,lIsoR05,lID);
746 aVec.push_back(lObj);
747 }
748 }
749 }
750
751 }//HbbMuons
752
753 void HbbTreeMaker::HbbTaus(const edm::Handle<std::vector<pat::Tau> > & aCol,
754 const edm::Handle<std::vector<reco::Vertex> > & aRecoVertices,
755 std::vector<HbbAnalysis::Tau> & aVec)
756 {//HbbTaus
757
758 if (aCol.isValid()){
759 if (aCol->size() > 0) {
760
761 unsigned int iEle = 0;
762 for (std::vector<pat::Tau>::const_iterator iter = aCol->begin();
763 iter != aCol->end();
764 ++iter,++iEle)
765 {
766 //if (debug_ > 1) std::cout << "**** Ele #" << iEle << ", pT,eta,phi = " << (*iter).pt() << " " << (*iter).eta() << " " << (*iter).phi() << std::endl;
767
768 HbbAnalysis::GenVars lGen;
769 if (!processData_ && (*iter).genLepton()){
770 lGen.valid = true;
771 lGen.E = (*iter).genLepton()->energy();
772 lGen.pT = (*iter).genLepton()->pt();
773 lGen.eta = (*iter).genLepton()->eta();
774 lGen.y = (*iter).genLepton()->y();
775 lGen.phi = (*iter).genLepton()->phi();
776 lGen.charge = (*iter).genLepton()->charge();
777 lGen.pdgId = (*iter).genLepton()->pdgId();
778 lGen.status = (*iter).genLepton()->status();
779 lGen.mass = (*iter).genLepton()->mass();
780 lGen.vx = (*iter).genLepton()->vx();
781 lGen.vy = (*iter).genLepton()->vy();
782 lGen.vz = (*iter).genLepton()->vz();
783 }
784 else {
785 lGen.valid = false;
786 lGen.E = 0;
787 lGen.pT = 0;
788 lGen.eta = 0;
789 lGen.y = 0;
790 lGen.phi = 0;
791 lGen.charge = 0;
792 lGen.pdgId = 0;
793 lGen.status = 0;
794 lGen.mass = 0;
795 lGen.vx = 0;
796 lGen.vy = 0;
797 lGen.vz = 0;
798 }
799
800
801 HbbAnalysis::GenVars lGenJet;
802 if (!processData_ && (*iter).genJet()){
803 lGenJet.valid = true;
804 lGenJet.E = (*iter).genJet()->energy();
805 lGenJet.pT = (*iter).genJet()->pt();
806 lGenJet.eta = (*iter).genJet()->eta();
807 lGenJet.y = (*iter).genJet()->y();
808 lGenJet.phi = (*iter).genJet()->phi();
809 lGenJet.charge = (*iter).genJet()->charge();
810 lGenJet.pdgId = (*iter).genJet()->pdgId();
811 lGenJet.status = (*iter).genJet()->status();
812 lGenJet.mass = (*iter).genJet()->mass();
813 lGenJet.vx = (*iter).genJet()->vx();
814 lGenJet.vy = (*iter).genJet()->vy();
815 lGenJet.vz = (*iter).genJet()->vz();
816 }
817 else {
818 lGenJet.valid = false;
819 lGenJet.E = 0;
820 lGenJet.pT = 0;
821 lGenJet.eta = 0;
822 lGenJet.y = 0;
823 lGenJet.phi = 0;
824 lGenJet.charge = 0;
825 lGenJet.pdgId = 0;
826 lGenJet.status = 0;
827 lGenJet.mass = 0;
828 lGenJet.vx = 0;
829 lGenJet.vy = 0;
830 lGenJet.vz = 0;
831 }
832
833 HbbAnalysis::BaseVars lReco;
834 lReco.E = (*iter).energy();
835 lReco.pT = (*iter).pt();
836 lReco.eta = (*iter).eta();
837 lReco.y = (*iter).y();
838 lReco.phi = (*iter).phi();
839 lReco.charge = (*iter).charge();
840 lReco.vx = (*iter).vx();
841 lReco.vy = (*iter).vy();
842 lReco.vz = (*iter).vz();
843
844 bool isCalo = (*iter).isCaloTau();
845 bool isPF = (*iter).isPFTau();
846
847 assert (isCalo == !isPF);
848
849 HbbAnalysis::TrkVars lLead;
850 if (isCalo) {
851 lLead.signedSipt = (*iter).leadTracksignedSipt();
852
853 reco::TrackRef lCaloTrk = (*iter).leadTrack();
854
855 if ( lCaloTrk.isAvailable() && lCaloTrk.isNonnull() ) {
856 lLead.pT = lCaloTrk->pt();
857 lLead.eta = lCaloTrk->eta();
858 lLead.phi = lCaloTrk->phi();
859
860 lLead.matchDist = reco::deltaR(lCaloTrk->momentum(), (*iter).p4());
861
862 if ( aRecoVertices->size() >= 1 ) {
863 const reco::Vertex& thePrimaryEventVertex = (*aRecoVertices->begin());
864 lLead.IPxy = lCaloTrk->dxy(thePrimaryEventVertex.position());
865 lLead.IPz = lCaloTrk->dz(thePrimaryEventVertex.position());
866 }
867 else {
868 lLead.IPxy = 0;
869 lLead.IPz = 0;
870 }
871 }
872 else {
873 lLead.pT = 0;
874 lLead.eta = 0;
875 lLead.phi = 0;
876 lLead.matchDist = 0;
877 lLead.IPxy = 0;
878 lLead.IPz = 0;
879 lLead.signedSipt = 0;
880 }
881 }//calo
882 else {
883 lLead.signedSipt = (*iter).leadPFChargedHadrCandsignedSipt();
884
885 reco::PFCandidateRef lHadrCand = (*iter).leadPFChargedHadrCand();
886
887 if ( lHadrCand.isAvailable() && lHadrCand.isNonnull() ) {
888 lLead.pT = lHadrCand->pt();
889 lLead.eta = lHadrCand->eta();
890 lLead.phi = lHadrCand->phi();
891
892 lLead.matchDist = reco::deltaR(lHadrCand->p4(), (*iter).p4());
893
894 reco::TrackRef lTrk = lHadrCand->trackRef();
895
896 if ( aRecoVertices->size() >= 1) {
897 const reco::Vertex& thePrimaryEventVertex = (*aRecoVertices->begin());
898 lLead.IPxy = lTrk->dxy(thePrimaryEventVertex.position());
899 lLead.IPz = lTrk->dz(thePrimaryEventVertex.position());
900 }
901 else {
902 lLead.IPxy = 0;
903 lLead.IPz = 0;
904 }
905 }
906 else {
907 lLead.pT = 0;
908 lLead.eta = 0;
909 lLead.phi = 0;
910 lLead.matchDist = 0;
911 lLead.IPxy = 0;
912 lLead.IPz = 0;
913 lLead.signedSipt = 0;
914 }
915
916 }//pf
917
918
919
920 const std::vector<std::pair<std::string,float> > & lIDs = (*iter).tauIDs();
921 bool lPrint = false;
922 //a few warning if additional discriminants are found:
923 if ((isPF && lIDs.size() != 16) ||
924 (isCalo && lIDs.size() != 3))
925 lPrint = true;
926
927 if (lPrint) {
928 std::cout << "!!!!!!! Discriminants changed, please update tree !!!!!!!" << std::endl;
929 std::cout << "--- isCaloTau = " << isCalo << ", isPFTau = " << isPF << std::endl;
930 std::cout << "------ ID names = " << std::endl;
931 }
932
933
934 HbbAnalysis::CaloTauIDVars lcaloId;
935 lcaloId.byIsolation = 0;
936 lcaloId.byLeadingTrackFinding = 0;
937 lcaloId.byLeadingTrackPtCut = 0;
938
939 HbbAnalysis::PFTauIDVars lpfId;
940 lpfId.byLeadingTrackFinding = 0;
941 lpfId.byLeadingTrackPtCut = 0;
942 lpfId.byTrackIsolation = 0;
943 lpfId.byECALIsolation = 0;
944 lpfId.byIsolation = 0;
945 lpfId.againstElectron = 0;
946 lpfId.againstMuon = 0;
947 lpfId.byIsolationUsingLeadingPion = 0;
948 lpfId.byTaNC = 0;
949 lpfId.byTaNCfrHalfPercent = 0;
950 lpfId.byTaNCfrOnePercent = 0;
951 lpfId.byTaNCfrQuarterPercent = 0;
952 lpfId.byTaNCfrTenthPercent = 0;
953 lpfId.ecalIsolationUsingLeadingPion = 0;
954 lpfId.leadingPionPtCut = 0;
955 lpfId.trackIsolationUsingLeadingPion = 0;
956
957
958
959 for (unsigned int id(0); id<lIDs.size(); ++id){
960
961 std::string lName = lIDs.at(id).first;
962 float lDiscri = lIDs.at(id).second;
963
964 if (lPrint) std::cout << "--------- " << lName << " = " << lDiscri << std::endl;
965 //pf
966
967 if (isPF) {
968 if (lName.find("leadingTrackFinding") != lName.npos) lpfId.byLeadingTrackFinding = lDiscri;
969 if (lName.find("leadingTrackPtCut") != lName.npos) lpfId.byLeadingTrackPtCut = lDiscri;
970 if (lName.find("trackIsolationUsingLeadingPion") != lName.npos) lpfId.trackIsolationUsingLeadingPion = lDiscri;
971 if (lName.find("trackIsolation") != lName.npos &&
972 lName.find("trackIsolationUsingLeadingPion") == lName.npos) lpfId.byTrackIsolation = lDiscri;
973 if (lName.find("ecalIsolationUsingLeadingPion") != lName.npos) lpfId.ecalIsolationUsingLeadingPion = lDiscri;
974 if (lName.find("ecalIsolation") != lName.npos &&
975 lName.find("ecalIsolationUsingLeadingPion") == lName.npos) lpfId.byECALIsolation = lDiscri;
976 if (lName.find("byIsolationUsingLeadingPion") != lName.npos) lpfId.byIsolationUsingLeadingPion = lDiscri;
977 if (lName.find("byIsolation") != lName.npos &&
978 lName.find("byIsolationUsingLeadingPion") == lName.npos) lpfId.byIsolation = lDiscri;
979 if (lName.find("againstElectron") != lName.npos) lpfId.againstElectron = lDiscri;
980 if (lName.find("againstMuon") != lName.npos) lpfId.againstMuon = lDiscri;
981 if (lName.find("byTaNCfrHalfPercent") != lName.npos) lpfId.byTaNCfrHalfPercent = lDiscri;
982 if (lName.find("byTaNCfrOnePercent") != lName.npos) lpfId.byTaNCfrOnePercent = lDiscri;
983 if (lName.find("byTaNCfrQuarterPercent") != lName.npos) lpfId.byTaNCfrQuarterPercent = lDiscri;
984 if (lName.find("byTaNCfrTenthPercent") != lName.npos) lpfId.byTaNCfrTenthPercent = lDiscri;
985 if (lName.find("byTaNC") != lName.npos &&
986 lName.find("byTaNCfr") == lName.npos) lpfId.byTaNC = lDiscri;
987 if (lName.find("leadingPionPtCut") != lName.npos) lpfId.leadingPionPtCut = lDiscri;
988 }
989 if (isCalo){
990 if (lName.find("byIsolation") != lName.npos) lcaloId.byIsolation = lDiscri;
991 if (lName.find("leadingTrackFinding") != lName.npos) lcaloId.byLeadingTrackFinding = lDiscri;
992 if (lName.find("leadingTrackPtCut") != lName.npos) lcaloId.byLeadingTrackPtCut = lDiscri;
993 }
994
995
996 }
997
998
999 if (isCalo){
1000
1001 HbbAnalysis::CaloTauIsoVars lcaloIso;
1002 lcaloIso.nIsoTracks = (*iter).isolationTracks().size();
1003 lcaloIso.nSigTracks = (*iter).signalTracks().size();
1004 lcaloIso.leadTrackHCAL3x3hitsEtSum = (*iter).leadTrackHCAL3x3hitsEtSum();
1005 lcaloIso.leadTrackHCAL3x3hottesthitDEta = (*iter).leadTrackHCAL3x3hottesthitDEta();
1006 lcaloIso.signalTracksInvariantMass = (*iter).signalTracksInvariantMass();
1007 lcaloIso.tracksInvariantMass = (*iter).TracksInvariantMass();
1008 lcaloIso.isolationTracksPtSum = (*iter).isolationTracksPtSum();
1009 lcaloIso.isolationECALhitsEtSum = (*iter).isolationECALhitsEtSum();
1010 lcaloIso.maximumHCALhitEt = (*iter).maximumHCALhitEt();
1011
1012 HbbAnalysis::Tau lObj(lGen,lGenJet,lReco,lLead,lcaloIso,lcaloId);
1013 aVec.push_back(lObj);
1014 }
1015 else {
1016
1017 HbbAnalysis::PFTauIsoVars lpfIso;
1018 lpfIso.nSigCands = (*iter).signalPFCands().size();
1019 lpfIso.nIsoCands = (*iter).isolationPFCands().size();
1020 lpfIso.maximumHCALPFClusterEt = (*iter).maximumHCALPFClusterEt();
1021 //lpfIso.emFraction = (*iter).emFraction();
1022 lpfIso.hcalTotOverPLead = (*iter).hcalTotOverPLead();
1023 lpfIso.hcalMaxOverPLead = (*iter).hcalMaxOverPLead();
1024 lpfIso.hcal3x3OverPLead = (*iter).hcal3x3OverPLead();
1025 lpfIso.ecalStripSumEOverPLead = (*iter).ecalStripSumEOverPLead();
1026 lpfIso.bremsRecoveryEOverPLead = (*iter).bremsRecoveryEOverPLead();
1027
1028 HbbAnalysis::PFTauCandVars lHadr;
1029 lHadr.nSigCands = (*iter).signalPFChargedHadrCands().size();
1030 lHadr.nIsoCands = (*iter).isolationPFChargedHadrCands().size();
1031 lHadr.isolationPtSum = (*iter).isolationPFChargedHadrCandsPtSum();
1032
1033 HbbAnalysis::PFTauCandVars lNeutr;
1034 lNeutr.nSigCands = (*iter).signalPFNeutrHadrCands().size();
1035 lNeutr.nIsoCands = (*iter).isolationPFNeutrHadrCands().size();
1036 lNeutr.isolationPtSum = (*iter).neutralHadronIso();
1037
1038 HbbAnalysis::PFTauCandVars lGamma;
1039 lGamma.nSigCands = (*iter).signalPFGammaCands().size();
1040 lGamma.nIsoCands = (*iter).isolationPFGammaCands().size();
1041 lGamma.isolationPtSum = (*iter).isolationPFGammaCandsEtSum();
1042
1043 HbbAnalysis::PFTauEleIDVars lEleId;
1044 if ( (*iter).electronPreIDTrack().isAvailable() && (*iter).electronPreIDTrack().isNonnull() ) {
1045 lEleId.pT = (*iter).electronPreIDTrack()->pt();
1046 lEleId.eta = (*iter).electronPreIDTrack()->eta();
1047 lEleId.phi = (*iter).electronPreIDTrack()->phi();
1048 }
1049 else {
1050 lEleId.pT = 0;
1051 lEleId.eta = 0;
1052 lEleId.phi = 0;
1053 }
1054
1055 lEleId.output = (*iter).electronPreIDOutput();
1056 lEleId.decision = (*iter).electronPreIDDecision();
1057
1058 HbbAnalysis::PFTauMuIDVars lMuId;
1059 lMuId.caloCompat = (*iter).caloComp();
1060 lMuId.segCompat = (*iter).segComp();
1061 lMuId.decision = (*iter).muonDecision();
1062
1063 HbbAnalysis::Tau lObj(lGen,lGenJet,lReco,lLead,lpfIso,lHadr,lNeutr,lGamma,lpfId,lEleId,lMuId);
1064 aVec.push_back(lObj);
1065 }
1066 }
1067 }
1068 }
1069
1070 }//HbbTaus
1071
1072 void HbbTreeMaker::HbbJets(const edm::Handle<std::vector<pat::Jet> > & aCol,
1073 const HbbAnalysis::JetFlavour & aJetFlav,
1074 const edm::Handle<reco::GenParticleCollection> & aGenParticles,
1075 std::vector<HbbAnalysis::Jet> & aVec)
1076 {//HbbJets
1077
1078 if (aCol.isValid()){//valid
1079 if (aCol->size() > 0) {//non empty
1080
1081 unsigned int iEle = 0;
1082 for (std::vector<pat::Jet>::const_iterator iter = aCol->begin();
1083 iter != aCol->end();
1084 ++iter,++iEle)
1085 {//loop on element
1086 //if (debug_ > 1) std::cout << "**** Ele #" << iEle << ", pT,eta,phi = " << (*iter).pt() << " " << (*iter).eta() << " " << (*iter).phi() << std::endl;
1087
1088 HbbAnalysis::GenVars lGen;
1089 if (!processData_ && (*iter).genParton()){
1090 lGen.valid = true;
1091 lGen.E = (*iter).genParton()->energy();
1092 lGen.pT = (*iter).genParton()->pt();
1093 lGen.eta = (*iter).genParton()->eta();
1094 lGen.y = (*iter).genParton()->y();
1095 lGen.phi = (*iter).genParton()->phi();
1096 lGen.charge = (*iter).genParton()->charge();
1097 lGen.pdgId = (*iter).genParton()->pdgId();
1098 lGen.status = (*iter).genParton()->status();
1099 lGen.mass = (*iter).genParton()->mass();
1100 lGen.vx = (*iter).genParton()->vx();
1101 lGen.vy = (*iter).genParton()->vy();
1102 lGen.vz = (*iter).genParton()->vz();
1103 }
1104 else {
1105 lGen.valid = false;
1106 lGen.E = 0;
1107 lGen.pT = 0;
1108 lGen.eta = 0;
1109 lGen.y = 0;
1110 lGen.phi = 0;
1111 lGen.charge = 0;
1112 lGen.pdgId = 0;
1113 lGen.status = 0;
1114 lGen.mass = 0;
1115 lGen.vx = 0;
1116 lGen.vy = 0;
1117 lGen.vz = 0;
1118 }
1119
1120
1121 HbbAnalysis::GenVars lGenJet;
1122 if (!processData_ && (*iter).genJet()){
1123 lGenJet.valid = true;
1124 lGenJet.E = (*iter).genJet()->energy();
1125 lGenJet.pT = (*iter).genJet()->pt();
1126 lGenJet.eta = (*iter).genJet()->eta();
1127 lGenJet.y = (*iter).genJet()->y();
1128 lGenJet.phi = (*iter).genJet()->phi();
1129 lGenJet.charge = (*iter).genJet()->charge();
1130 lGenJet.pdgId = (*iter).genJet()->pdgId();
1131 lGenJet.status = (*iter).genJet()->status();
1132 lGenJet.mass = (*iter).genJet()->mass();
1133 lGenJet.vx = (*iter).genJet()->vx();
1134 lGenJet.vy = (*iter).genJet()->vy();
1135 lGenJet.vz = (*iter).genJet()->vz();
1136 }
1137 else {
1138 lGenJet.valid = false;
1139 lGenJet.E = 0;
1140 lGenJet.pT = 0;
1141 lGenJet.eta = 0;
1142 lGenJet.y = 0;
1143 lGenJet.phi = 0;
1144 lGenJet.charge = 0;
1145 lGenJet.pdgId = 0;
1146 lGenJet.status = 0;
1147 lGenJet.mass = 0;
1148 lGenJet.vx = 0;
1149 lGenJet.vy = 0;
1150 lGenJet.vz = 0;
1151 }
1152
1153 HbbAnalysis::BaseVars lReco;
1154 lReco.E = (*iter).energy();
1155 lReco.pT = (*iter).pt();
1156 lReco.eta = (*iter).eta();
1157 //lReco.etaDet = (*iter).detectorEta((*iter).vz(),(*iter).eta());
1158 lReco.y = (*iter).y();
1159 lReco.phi = (*iter).phi();
1160 lReco.charge = (*iter).jetCharge();
1161 lReco.vx = (*iter).vx();
1162 lReco.vy = (*iter).vy();
1163 lReco.vz = (*iter).vz();
1164 //lReco.myEtaDet = HbbAnalysis::EtaDetector(lReco);
1165
1166 unsigned int lFlav = 0;
1167 if (!processData_ && aJetFlav.nPartons()) {
1168 lFlav = aJetFlav.partonMatchingGenJet((*iter),aGenParticles,0.4).second;
1169 }
1170
1171 HbbAnalysis::JetVars lCommon;
1172 lCommon.flavour = lFlav;
1173 lCommon.partonFlavour = (*iter).partonFlavour();
1174 lCommon.nAssociatedTracks = ((*iter).associatedTracks()).size();
1175 lCommon.rawpT = (*iter).pt();
1176 if ((*iter).hasCorrFactors() && (*iter).corrFactor("raw")>0)
1177 lCommon.rawpT = (*iter).pt()/((*iter).corrFactor("raw"));
1178
1179 lCommon.etaMean = (*iter).etaPhiStatistics().etaMean;
1180 lCommon.phiMean = (*iter).etaPhiStatistics().phiMean;
1181 lCommon.etaEtaMoment = (*iter).etaPhiStatistics().etaEtaMoment;
1182 lCommon.phiPhiMoment = (*iter).etaPhiStatistics().phiPhiMoment;
1183 lCommon.etaPhiMoment = (*iter).etaPhiStatistics().etaPhiMoment;
1184
1185 lCommon.l1Match = false;//TODO
1186 lCommon.hltMatch = false;//TODO
1187
1188 //lCommon.rawEta = (*iter).;
1189 //lCommon.rawPhi = (*iter).;
1190 // p_hasJetCorrFactors->Fill(aJet.hasJetCorrFactors());
1191
1192 HbbAnalysis::JetBtagVars lBtag;
1193 lBtag.cSV = (*iter).bDiscriminator("combinedSecondaryVertexBJetTags");
1194 lBtag.cSVMVA = (*iter).bDiscriminator("combinedSecondaryVertexMVABJetTags");
1195 lBtag.iPMVA = (*iter).bDiscriminator("impactParameterMVABJetTags");
1196 lBtag.bProba = (*iter).bDiscriminator("jetBProbabilityBJetTags");
1197 lBtag.probability = (*iter).bDiscriminator("jetProbabilityBJetTags");
1198 lBtag.sSVHP = (*iter).bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
1199 lBtag.sSVHE = (*iter).bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
1200 lBtag.softElectronByPt = ((*iter).bDiscriminator("softElectronByPtBJetTags") < -1000) ? -1000 : (*iter).bDiscriminator("softElectronByPtBJetTags");
1201 lBtag.softElectronByIP3d = ((*iter).bDiscriminator("softElectronByIP3dBJetTags") < -1000) ? -1000 : (*iter).bDiscriminator("softElectronByIP3dBJetTags");
1202 lBtag.softMuon = ((*iter).bDiscriminator("softMuonBJetTags") < -1000) ? -1000 : (*iter).bDiscriminator("softMuonBJetTags");
1203 lBtag.softMuonByPt = ((*iter).bDiscriminator("softMuonByPtBJetTags") < -1000) ? -1000 : (*iter).bDiscriminator("softMuonByPtBJetTags");
1204 lBtag.softMuonByIP3d = ((*iter).bDiscriminator("softMuonByIP3dBJetTags") < -1000) ? -1000 : (*iter).bDiscriminator("softMuonByIP3dBJetTags");
1205
1206 //std::cout << " -- New values of b-discri for jet : " << iEle << std::endl
1207 // << " ---- softElecs: " << lBtag.softElectronByPt << " " << lBtag.softElectronByIP3d << std::endl
1208 // << " ---- softMus: " << lBtag.softMuon << " " << lBtag.softMuonByPt << " " << lBtag.softMuonByIP3d << std::endl;
1209
1210
1211 lBtag.tCHE = (*iter).bDiscriminator("trackCountingHighEffBJetTags");
1212 lBtag.tCHP = (*iter).bDiscriminator("trackCountingHighPurBJetTags");
1213
1214
1215 bool isCalo = (*iter).isCaloJet();
1216 bool isPF = (*iter).isPFJet();
1217 bool isJPT = (*iter).isJPTJet();
1218
1219 //std::cout << " -- isCalo = " << isCalo << ", isPF = " << isPF << std::endl;
1220
1221 //assert (isCalo == !isPF);
1222 HbbAnalysis::JetIDVars lId;
1223 lId.fHPD = 0;
1224 lId.fRBX = 0;
1225 lId.n90Hits = 0;
1226 lId.hitsInN90 = 0;
1227 lId.restrictedEMF = 0;
1228 lId.nHCALTowers = 0;
1229 lId.nECALTowers = 0;
1230 if (isCalo || isJPT) {
1231 lId.fHPD = (*iter).jetID().fHPD;
1232 lId.fRBX = (*iter).jetID().fRBX;
1233 lId.n90Hits = (*iter).jetID().n90Hits;
1234 lId.hitsInN90 = (*iter).jetID().hitsInN90;
1235 lId.restrictedEMF = (*iter).jetID().restrictedEMF;
1236 lId.nHCALTowers = (*iter).jetID().nHCALTowers;
1237 lId.nECALTowers = (*iter).jetID().nECALTowers;
1238 }
1239 if (isCalo){
1240 HbbAnalysis::CaloJetVars lCalo;
1241 lCalo.maxEInEmTowers = (*iter).maxEInEmTowers();
1242 lCalo.maxEInHadTowers = (*iter).maxEInHadTowers();
1243 lCalo.energyFractionHadronic = (*iter).energyFractionHadronic();
1244 lCalo.emEnergyFraction = (*iter).emEnergyFraction();
1245 lCalo.hadEnergyInHB = (*iter).hadEnergyInHB();
1246 lCalo.hadEnergyInHO = (*iter).hadEnergyInHO();
1247 lCalo.hadEnergyInHE = (*iter).hadEnergyInHE();
1248 lCalo.hadEnergyInHF = (*iter).hadEnergyInHF();
1249 lCalo.emEnergyInEB = (*iter).emEnergyInEB();
1250 lCalo.emEnergyInEE = (*iter).emEnergyInEE();
1251 lCalo.emEnergyInHF = (*iter).emEnergyInHF();
1252 lCalo.towersArea = (*iter).towersArea();
1253 lCalo.n90 = (*iter).n90();
1254 lCalo.n60 = (*iter).n60();
1255
1256 HbbAnalysis::Jet lObj(lGen,lGenJet,lReco,lCommon,lCalo,lBtag,lId);
1257 aVec.push_back(lObj);
1258 }
1259 HbbAnalysis::JPTPFJetVars lJPTPF;
1260 lJPTPF.chargedHadronEnergy = 0;
1261 lJPTPF.chargedHadronEnergyFraction = 0;
1262 lJPTPF.neutralHadronEnergy = 0;
1263 lJPTPF.neutralHadronEnergyFraction = 0;
1264 lJPTPF.chargedEmEnergy = 0;
1265 lJPTPF.chargedEmEnergyFraction = 0;
1266 lJPTPF.neutralEmEnergy = 0;
1267 lJPTPF.neutralEmEnergyFraction = 0;
1268 lJPTPF.chargedMultiplicity = 0;
1269 lJPTPF.muonMultiplicity = 0;
1270
1271 if (isPF || isJPT){
1272 lJPTPF.chargedEmEnergy = (*iter).chargedEmEnergy();
1273 lJPTPF.neutralEmEnergy = (*iter).neutralEmEnergy();
1274 lJPTPF.muonMultiplicity = (*iter).muonMultiplicity();
1275 lJPTPF.chargedHadronEnergy = (*iter).chargedHadronEnergy();
1276 lJPTPF.neutralHadronEnergy = (*iter).neutralHadronEnergy();
1277 lJPTPF.chargedMultiplicity = (*iter).chargedMultiplicity();
1278
1279 lJPTPF.chargedHadronEnergyFraction = (*iter).chargedHadronEnergyFraction();
1280 lJPTPF.neutralHadronEnergyFraction = (*iter).neutralHadronEnergyFraction();
1281 lJPTPF.chargedEmEnergyFraction = (*iter).chargedEmEnergyFraction();
1282 lJPTPF.neutralEmEnergyFraction = (*iter).neutralEmEnergyFraction();
1283 }
1284 if (isJPT) {
1285 HbbAnalysis::JPTJetVars lJPT;
1286 // const reco::TrackRefVector& pionsInVertexInCalo () const{return jptSpecific().pionsInVertexInCalo; }
1287 // const reco::TrackRefVector& pionsInVertexOutCalo() const{return jptSpecific().pionsInVertexOutCalo;}
1288 // const reco::TrackRefVector& pionsOutVertexInCalo() const{return jptSpecific().pionsOutVertexInCalo;}
1289 // const reco::TrackRefVector& muonsInVertexInCalo () const{return jptSpecific().muonsInVertexInCalo; }
1290 // const reco::TrackRefVector& muonsInVertexOutCalo() const{return jptSpecific().muonsInVertexOutCalo;}
1291 // const reco::TrackRefVector& muonsOutVertexInCalo() const{return jptSpecific().muonsOutVertexInCalo;}
1292 // const reco::TrackRefVector& elecsInVertexInCalo () const{return jptSpecific().elecsInVertexInCalo; }
1293 // const reco::TrackRefVector& elecsInVertexOutCalo() const{return jptSpecific().elecsInVertexOutCalo;}
1294 // const reco::TrackRefVector& elecsOutVertexInCalo() const{return jptSpecific().elecsOutVertexInCalo;}
1295 lJPT.zspCorrection = (*iter).zspCorrection();
1296 lJPT.elecMultiplicity = (*iter).elecMultiplicity();
1297
1298
1299 HbbAnalysis::CaloJetVars lCalo;
1300 const reco::JPTJet * lJetRef = dynamic_cast<const reco::JPTJet *>((*iter).originalObject());
1301 edm::RefToBase<reco::Jet> jptjetRef = lJetRef->getCaloJetRef();
1302
1303 // Dynamic_cast to CaloJet object
1304 reco::CaloJet const * rawcalojet =
1305 dynamic_cast<reco::CaloJet const *>( &* jptjetRef);
1306
1307 lCalo.maxEInEmTowers = rawcalojet->maxEInEmTowers();
1308 lCalo.maxEInHadTowers = rawcalojet->maxEInHadTowers();
1309 lCalo.energyFractionHadronic = rawcalojet->energyFractionHadronic();
1310 lCalo.emEnergyFraction = rawcalojet->emEnergyFraction();
1311 lCalo.hadEnergyInHB = rawcalojet->hadEnergyInHB();
1312 lCalo.hadEnergyInHO = rawcalojet->hadEnergyInHO();
1313 lCalo.hadEnergyInHE = rawcalojet->hadEnergyInHE();
1314 lCalo.hadEnergyInHF = rawcalojet->hadEnergyInHF();
1315 lCalo.emEnergyInEB = rawcalojet->emEnergyInEB();
1316 lCalo.emEnergyInEE = rawcalojet->emEnergyInEE();
1317 lCalo.emEnergyInHF = rawcalojet->emEnergyInHF();
1318 lCalo.towersArea = rawcalojet->towersArea();
1319 lCalo.n90 = rawcalojet->n90();
1320 lCalo.n60 = rawcalojet->n60();
1321
1322 lCommon.etaMean = rawcalojet->etaPhiStatistics().etaMean;
1323 lCommon.phiMean = rawcalojet->etaPhiStatistics().phiMean;
1324 lCommon.etaEtaMoment = rawcalojet->etaPhiStatistics().etaEtaMoment;
1325 lCommon.phiPhiMoment = rawcalojet->etaPhiStatistics().phiPhiMoment;
1326 lCommon.etaPhiMoment = rawcalojet->etaPhiStatistics().etaPhiMoment;
1327
1328 HbbAnalysis::Jet lObj(lGen,lGenJet,lReco,lCommon,lCalo,lJPT,lJPTPF,lBtag,lId);
1329 aVec.push_back(lObj);
1330 }
1331 if (isPF) {
1332 HbbAnalysis::PFJetVars lPF;
1333 lPF.chargedMuEnergy = (*iter).chargedMuEnergy();
1334 lPF.chargedMuEnergyFraction = (*iter).chargedMuEnergyFraction();
1335 lPF.neutralMultiplicity = (*iter).neutralMultiplicity();
1336
1337 HbbAnalysis::Jet lObj(lGen,lGenJet,lReco,lCommon,lPF,lJPTPF,lBtag);
1338 aVec.push_back(lObj);
1339 }
1340
1341 }//loop on element
1342 }//non empty
1343 }//valid
1344
1345 }//HbbJets
1346
1347 void HbbTreeMaker::HbbMet(const edm::Handle<std::vector<pat::MET> > & aCol,
1348 HbbAnalysis::Met & aMet)
1349 {//HbbMet
1350 HbbAnalysis::MetVars lGen;
1351 HbbAnalysis::MetVars lReco;
1352
1353 assert(aCol->size() == 1);
1354 const pat::MET & lMet = *(aCol->begin());
1355
1356 lReco.mET = lMet.pt();
1357 lReco.mEx = lMet.px();
1358 lReco.mEy = lMet.py();
1359 lReco.sumET = lMet.sumEt();
1360 lReco.phi = lMet.phi();
1361 lReco.mEtSig = lMet.mEtSig();
1362
1363 const reco::GenMET *lGenMET = lMet.genMET();
1364
1365 if (!processData_ && lGenMET){
1366 lGen.mET = lGenMET->pt();
1367 lGen.mEx = lGenMET->px();
1368 lGen.mEy = lGenMET->py();
1369 lGen.sumET = lGenMET->sumEt();
1370 lGen.phi = lGenMET->phi();
1371 lGen.mEtSig = lGenMET->mEtSig();
1372 }
1373 else {
1374 lGen.mET = 0;
1375 lGen.mEx = 0;
1376 lGen.mEy = 0;
1377 lGen.sumET = 0;
1378 lGen.phi = 0;
1379 lGen.mEtSig = 0;
1380 }
1381
1382 aMet.genVars(lGen);
1383 aMet.recoVars(lReco);
1384
1385 }//HbbMet
1386
1387 void HbbTreeMaker::HbbTrigger(const edm::Handle<edm::TriggerResults> & aCol,
1388 const edm::TriggerNames & aNames,
1389 std::vector<HbbAnalysis::Trigger> & aVec)
1390 {//HbbTrigger
1391
1392 if (aCol.isValid()){
1393 //edm::TriggerNames lNames;
1394 //lNames.init(*aCol);
1395
1396 //for (unsigned int j=0; j<hltPaths_.size(); j++) {
1397 //bool valid=false;
1398 for (unsigned int i=0; i<aCol->size(); ++i){
1399 std::string trueName = aNames.triggerName(i);//.at(i);
1400
1401 HbbAnalysis::TriggerVars lTrigVars;
1402 lTrigVars.name = trueName;
1403 lTrigVars.index = i;
1404 lTrigVars.accept = aCol->accept(i);
1405
1406 HbbAnalysis::Trigger lObj(lTrigVars);
1407
1408 aVec.push_back(lObj);
1409 }
1410
1411 }
1412
1413 }
1414
1415 void HbbTreeMaker::HbbParticles(const edm::Handle<reco::GenParticleCollection> & aCol,
1416 std::vector<HbbAnalysis::GenParticle> & aVec)
1417 {//genparticles
1418
1419 //add genParticle inside vector
1420 for (unsigned int mccount = 0;mccount<aCol->size();++mccount)
1421 {//loop on particles
1422 const reco::Candidate & p = (*aCol)[mccount];
1423
1424 HbbAnalysis::MCVars lMC;
1425 lMC.index = mccount;
1426 lMC.E = p.energy();
1427 lMC.pT = p.pt();
1428 lMC.eta = p.eta();
1429 lMC.y = p.y();
1430 lMC.phi = p.phi();
1431 lMC.pdgId = p.pdgId();
1432 lMC.status = p.status();
1433 HbbAnalysis::GenParticle lObj(lMC);
1434 aVec.push_back(lObj);
1435
1436 }//loop on particles
1437
1438 //now need to get the list of parents....
1439
1440 for (unsigned int mccount = 0;mccount<aCol->size();++mccount)
1441 {//loop on particles
1442 const reco::Candidate & p = (*aCol)[mccount];
1443
1444 unsigned int nD = p.numberOfDaughters();
1445 //std::cout << mccount << " has " << nD << " daughter(s) : " ;
1446 for(unsigned int dau=0;dau<nD;++dau){//loop on daughters
1447 const reco::Candidate * pDau = p.daughter(dau);
1448 //std::cout << pDau << " ("<<std::flush;
1449 for (unsigned int mccount2 = 0;mccount2<aCol->size();++mccount2)
1450 {//loop on particles
1451 const reco::Candidate & p2 = (*aCol)[mccount2];
1452 //std::cout<<"DBG: mccount2 = "<<mccount2<<" gen size = "<<data_->gen().size()<<std::endl;
1453
1454 HbbAnalysis::GenParticle & part2 = aVec.at(mccount2);
1455 if (pDau == &p2) {
1456 part2.setParent(mccount);
1457 //std::cout << &p2 << ", index = " << mccount2 << "), "<<std::flush;
1458 break;
1459 }
1460 //else std::cout << std::endl << "****no match " << mccount2 << " " << &p2 << std::endl;
1461 //if (p2.pdgId() == pDau->pdgId() && fabs(p2.energy()-pDau->energy())<0.000001 && fabs(p2.eta()-pDau->eta())<0.000001 && fabs(p2.phi()-pDau->phi())<0.000001){
1462 //part2.parent(mccount);
1463 //}
1464 }//loop on particles
1465 }//loop on daughters
1466 //std::cout << std::endl;
1467
1468 }//loop on particles
1469
1470 if (debug_){
1471 for (unsigned int imc(0); imc < aVec.size(); ++imc){
1472 HbbAnalysis::GenParticle & p = aVec.at(imc);
1473 p.print();
1474 }
1475 }
1476
1477
1478 }//genparticles
1479
1480
1481
1482 void HbbTreeMaker::HbbVertices(const edm::Handle<std::vector<reco::Vertex> > & aCol,
1483 std::vector<HbbAnalysis::Vertex> & aVec)
1484 {
1485
1486 for (std::vector<reco::Vertex>::const_iterator iter = aCol->begin();
1487 iter != aCol->end();
1488 ++iter)
1489 {
1490
1491 if (!((*iter).isValid()) || (*iter).isFake()) continue;
1492
1493 HbbAnalysis::VertexVars lVtx;
1494
1495 if ((*iter).tracksSize() > 0) {
1496 for (std::vector<reco::TrackBaseRef>::const_iterator lTrk = (*iter).tracks_begin();
1497 lTrk != (*iter).tracks_end();
1498 ++lTrk) {
1499 lVtx.trackWeights.push_back((*iter).trackWeight(*lTrk));
1500 }
1501 }
1502
1503 if (lVtx.trackWeights.size() != (*iter).tracksSize()) {
1504 std::cout<< " -- Problem with tracks, size is not as expected ! "
1505 << std::endl
1506 << " --- Size of recoVertex tracks = " << (*iter).tracksSize()
1507 << std::endl
1508 << " --- Size of trackWeights = " << lVtx.trackWeights.size()
1509 << std::endl;
1510 }
1511
1512 lVtx.chi2 = (*iter).chi2();
1513 lVtx.ndof = (*iter).ndof();
1514 lVtx.x = (*iter).x();
1515 lVtx.y = (*iter).y();
1516 lVtx.z = (*iter).z();
1517 lVtx.rho = (*iter).position().Rho();
1518 lVtx.xError = (*iter).xError();
1519 lVtx.yError = (*iter).yError();
1520 lVtx.zError = (*iter).zError();
1521 lVtx.cov01 = (*iter).covariance(0,1);
1522 lVtx.cov02 = (*iter).covariance(0,2);
1523 lVtx.cov12 = (*iter).covariance(1,2);
1524
1525 HbbAnalysis::Vertex lObj(lVtx);
1526 aVec.push_back(lObj);
1527 }
1528
1529 }
1530
1531 void HbbTreeMaker::HbbL1Objects(const edm::Handle<l1extra::L1JetParticleCollection> & aCol,
1532 std::vector<HbbAnalysis::L1Object> & aVec,
1533 const unsigned int aType)
1534 {
1535
1536 for (l1extra::L1JetParticleCollection::const_iterator iter = aCol->begin();
1537 iter != aCol->end();
1538 ++iter)
1539 {
1540
1541 HbbAnalysis::L1Vars lVars;
1542 lVars.pT = (*iter).pt();
1543 lVars.ET = (*iter).et();
1544 lVars.eta = (*iter).eta();
1545 lVars.phi = (*iter).phi();
1546 lVars.bx = (*iter).bx();
1547 lVars.type = aType;
1548
1549 HbbAnalysis::L1Object lObj(lVars);
1550 aVec.push_back(lObj);
1551
1552 }
1553
1554 }
1555
1556
1557 void HbbTreeMaker::HbbHLTObjects(const edm::Handle<trigger::TriggerEvent> & aCol,
1558 std::vector<HbbAnalysis::HLTObject> & aVec)
1559 {
1560
1561
1562 //cout
1563 //const trigger::TriggerObjectCollection & toc(aCol->getObjects());
1564 //for (unsigned i=0; i<aCol->sizeFilters(); i++){
1565 //const edm::InputTag & lTag = aCol->filterTag(i);
1566 //std::cout << " -- ele " << i << ": " << lTag << ", key ids = " ;
1567 ////const int index = aCol->filterIndex(lTag);
1568 ////if (index >= aCol->sizeFilters()) {
1569 ////std::cout << " -- Filter " << aTag << " not found !" << std::endl;
1570 ////continue;
1571 ////}
1572 //const trigger::Keys & k = aCol->filterKeys(i);
1573 //for (trigger::Keys::const_iterator ki = k.begin(); ki !=k.end(); ++ki ) {
1574 //std::cout << toc[*ki].id() << ", " ;
1575 //}
1576 //std::cout << std::endl;
1577 //}
1578
1579 for (unsigned i(0); i<hltTagsElec_.size(); ++i){
1580 fillHLTVector(hltTagsElec_[i],aCol,aVec,1);
1581 }
1582 for (unsigned i(0); i<hltTagsMu_.size(); ++i){
1583 fillHLTVector(hltTagsMu_[i],aCol,aVec,2);
1584 }
1585 for (unsigned i(0); i<hltTagsJet_.size(); ++i){
1586 fillHLTVector(hltTagsJet_[i],aCol,aVec,3);
1587 }
1588 }
1589
1590 void HbbTreeMaker::fillHLTVector(const edm::InputTag & aTag,
1591 const edm::Handle<trigger::TriggerEvent> & aCol,
1592 std::vector<HbbAnalysis::HLTObject> & aVec,
1593 const unsigned int aType){
1594
1595 const trigger::TriggerObjectCollection & toc(aCol->getObjects());
1596 const int index = aCol->filterIndex(aTag);
1597 if (index >= aCol->sizeFilters()) {
1598 //std::cout << " -- Filter " << aTag << " not found !" << std::endl;
1599 return;
1600 }
1601
1602 const trigger::Keys & k = aCol->filterKeys(index);
1603
1604 for (trigger::Keys::const_iterator ki = k.begin(); ki !=k.end(); ++ki ) {
1605 HbbAnalysis::HLTVars lVars;
1606 lVars.pT = toc[*ki].pt();
1607 lVars.eta = toc[*ki].eta();
1608 lVars.phi = toc[*ki].phi();
1609 lVars.id = toc[*ki].id();
1610 lVars.mass = toc[*ki].mass();
1611 lVars.type = aType;
1612 HbbAnalysis::HLTObject lObj(lVars);
1613 aVec.push_back(lObj);
1614 }
1615 }
1616
1617 void HbbTreeMaker::HbbGenInfo(const edm::Handle<edm::HepMCProduct> & amcProd,
1618 const edm::Handle<GenRunInfoProduct> & aRunProd)
1619 {
1620 const HepMC::GenEvent * genEvt = amcProd->GetEvent();
1621 HbbAnalysis::GenInfo lInfo;
1622
1623 lInfo.processID(genEvt->signal_process_id());
1624 lInfo.ptHat(genEvt->event_scale());
1625
1626 lInfo.internalXS(aRunProd->internalXSec().value(),
1627 aRunProd->internalXSec().error());
1628
1629 lInfo.externalXSLO(aRunProd->externalXSecLO().value(),
1630 aRunProd->externalXSecLO().error());
1631
1632 lInfo.externalXSNLO(aRunProd->externalXSecNLO().value(),
1633 aRunProd->externalXSecNLO().error());
1634
1635 lInfo.filterEff(aRunProd->filterEfficiency());
1636
1637 event_->generatorInfo(lInfo);
1638
1639 }
1640
1641 #include "FWCore/Framework/interface/MakerMacros.h"
1642 DEFINE_FWK_MODULE(HbbTreeMaker);
1643
1644