ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbbAnalysis/VHbbDataFormats/bin/Ntupler7TeV.cc
Revision: 1.2
Committed: Thu May 24 15:58:24 2012 UTC (12 years, 11 months ago) by bortigno
Content type: text/plain
Branch: MAIN
CVS Tags: Step2ForV32_v1, Step2ForV32_v0, Step2ForV31_v0
Changes since 1.1: +64 -1 lines
Log Message: 
PDFweight added for 7TeV->8TeV reweighting. It works only with LHE based geneators. For the others generators and data PDFweight=1.

File Contents

# Content
1 #include <TH1F.h>
2 #include <TH3F.h>
3 #include "PhysicsTools/Utilities/interface/LumiReWeighting.h"
4 #include "PhysicsTools/Utilities/interface/Lumi3DReWeighting.h"
5 #include <TH2F.h>
6 #include <TROOT.h>
7 #include <TFile.h>
8 #include <TTree.h>
9 #include <TSystem.h>
10 #include "DataFormats/FWLite/interface/Event.h"
11 #include "DataFormats/FWLite/interface/Handle.h"
12 #include "FWCore/FWLite/interface/AutoLibraryLoader.h"
13 #include "DataFormats/MuonReco/interface/Muon.h"
14 #include "DataFormats/PatCandidates/interface/Muon.h"
15 #include "PhysicsTools/FWLite/interface/TFileService.h"
16 #include "FWCore/ParameterSet/interface/ProcessDesc.h"
17 #include "FWCore/PythonParameterSet/interface/PythonProcessDesc.h"
18 #include "DataFormats/Math/interface/deltaR.h"
19 #include "DataFormats/Math/interface/deltaPhi.h"
20
21 #include "DataFormats/FWLite/interface/LuminosityBlock.h"
22 #include "DataFormats/FWLite/interface/Run.h"
23 #include "DataFormats/Luminosity/interface/LumiSummary.h"
24
25 #include "VHbbAnalysis/VHbbDataFormats/interface/HbbCandidateFinderAlgo.h"
26 #include "VHbbAnalysis/VHbbDataFormats/src/HbbCandidateFinderAlgo.cc"
27
28 #include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEvent.h"
29 #include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEventAuxInfo.h"
30 #include "VHbbAnalysis/VHbbDataFormats/interface/VHbbCandidate.h"
31 #include "VHbbAnalysis/VHbbDataFormats/interface/TriggerReader.h"
32 #include "VHbbAnalysis/VHbbDataFormats/interface/TopMassReco.h"
33
34 //for IVF
35 #include "RecoBTag/SecondaryVertex/interface/SecondaryVertex.h"
36 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
37 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
38 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
39 #include <DataFormats/GeometrySurface/interface/Surface.h>
40 #include "Math/SMatrix.h"
41
42 //for LHE info
43 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
44 #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
45 //pdf reweighting
46
47
48 //Move class definition to Ntupler.h ?
49 //#include "VHbbAnalysis/VHbbDataFormats/interface/Ntupler.h"
50
51 #include "ZSV/BAnalysis/interface/SimBHadron.h"
52 //btagging
53 #include "VHbbAnalysis/VHbbDataFormats/interface/BTagWeight.h"
54 //trigger
55 #include "VHbbAnalysis/VHbbDataFormats/interface/TriggerWeight.h"
56
57 #include <sstream>
58 #include <string>
59
60 #define MAXJ 30
61 #define MAXL 10
62 #define MAXB 10
63 #define MAXT 60
64 #define nMetUnc 24
65 //eleEnDown/Up, muEn, tauEn, JES, JER, Unclustered for type1 MET [0-11] and than type1p2 MET [12-23]
66
67 struct CompareDeltaR {
68 CompareDeltaR(TLorentzVector p4dir_): p4dir(p4dir_) {}
69 bool operator()( const VHbbEvent::SimpleJet& j1, const VHbbEvent::SimpleJet& j2 ) const {
70 return j1.p4.DeltaR(p4dir) > j2.p4.DeltaR(p4dir);
71 }
72 TLorentzVector p4dir;
73 };
74
75 const GlobalVector flightDirection(const TVector3 pv, const reco::Vertex &sv){
76 GlobalVector fdir(sv.position().X() - pv.X(),
77 sv.position().Y() - pv.Y(),
78 sv.position().Z() - pv.Z());
79 return fdir;
80 }
81
82 bool jsonContainsEvent (const std::vector< edm::LuminosityBlockRange > &jsonVec,
83 const edm::EventBase &event)
84 {
85 // if the jsonVec is empty, then no JSON file was provided so all
86 // events should pass
87 if (jsonVec.empty())
88 {
89 return true;
90 }
91 bool (* funcPtr) (edm::LuminosityBlockRange const &,
92 edm::LuminosityBlockID const &) = &edm::contains;
93 edm::LuminosityBlockID lumiID (event.id().run(),
94 event.id().luminosityBlock());
95 std::vector< edm::LuminosityBlockRange >::const_iterator iter =
96 std::find_if (jsonVec.begin(), jsonVec.end(),
97 boost::bind(funcPtr, _1, lumiID) );
98 return jsonVec.end() != iter;
99
100 }
101
102
103 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
104 namespace LHAPDF {
105 void initPDFSet(int nset, int setid, int member=0);
106 int numberPDF(int nset);
107 void usePDFMember(int nset, int member);
108 double xfx(int nset, double x, double Q, int fl);
109 double getXmin(int nset, int member);
110 double getXmax(int nset, int member);
111 double getQ2min(int nset, int member);
112 double getQ2max(int nset, int member);
113 void extrapolate(bool extrapolate=true);
114 }
115
116
117 float resolutionBias(float eta)
118 {
119 // return 0;//Nominal!
120 if(eta< 1.1) return 0.05;
121 if(eta< 2.5) return 0.10;
122 if(eta< 5) return 0.30;
123 return 0;
124 }
125
126
127 typedef struct
128 {
129 void set(const SimBHadron & sbhc, int i){
130 mass[i] = sbhc.mass();
131 pt[i] = sbhc.pt();
132 eta[i] = sbhc.eta();
133 phi[i] = sbhc.phi();
134 vtx_x[i] = sbhc.decPosition.x();
135 vtx_y[i] = sbhc.decPosition.y();
136 vtx_z[i] = sbhc.decPosition.z();
137 pdgId[i] = sbhc.pdgId();
138 status[i] = sbhc.status();
139 };
140 void reset(){
141 for(int i=0; i < MAXB; ++i){
142 mass[i] = -99; pt[i] = -99; eta[i] = -99; phi[i] = -99; vtx_x[i] = -99; vtx_y[i] = -99; vtx_z[i] = -99; pdgId[i] = -99; status[i] = -99;
143 }
144 };
145 float mass[MAXB];
146 float pt[MAXB];
147 float eta[MAXB];
148 float phi[MAXB];
149 float vtx_x[MAXB];
150 float vtx_y[MAXB];
151 float vtx_z[MAXB];
152 int pdgId[MAXB];
153 int status[MAXB];
154 // int quarkStatus[MAXB];
155 // int brotherStatus[MAXB];
156 // int otherId[MAXB];
157 // bool etaOk[MAXB];
158 // bool simOk[MAXB];
159 // bool trackOk[MAXB];
160 // bool cutOk[MAXB];
161 // bool cutNewOk[MAXB];
162 // bool mcMatchOk[MAXB];
163 // bool matchOk[MAXB];
164 } SimBHadronInfo;
165
166
167 typedef struct
168 {
169 void set( const reco::SecondaryVertex & recoSv, const TVector3 recoPv, int isv){
170 pt[isv] = recoSv.p4().Pt();
171 eta[isv] = flightDirection(recoPv,recoSv).eta();
172 phi[isv] = flightDirection(recoPv,recoSv).phi();
173 massBcand[isv] = recoSv.p4().M();
174 massSv[isv] = recoSv.p4().M();
175 dist3D[isv] = recoSv.dist3d().value();
176 distSig3D[isv] = recoSv.dist3d().significance();
177 dist2D[isv] = recoSv.dist2d().value();
178 distSig2D[isv] = recoSv.dist2d().significance();
179 dist3D_norm[isv] = recoSv.dist3d().value()/recoSv.p4().Gamma();
180 };
181 void reset(){
182 for(int i = 0; i < MAXB; ++i){
183 massBcand[i] = -99; massSv[i]= -99; pt[i] = -99; eta[i] = -99; phi[i] = -99; dist3D[i] = -99; distSig3D[i] = -99; dist2D[i] = -99; distSig2D[i] = -99; dist3D_norm[i] = -99;
184 }
185 };
186 float massBcand[MAXB];
187 float massSv[MAXB];
188 float pt[MAXB];
189 float eta[MAXB];
190 float phi[MAXB];
191 float dist3D[MAXB];
192 float distSig3D[MAXB];
193 float dist2D[MAXB];
194 float distSig2D[MAXB];
195 float dist3D_norm[MAXB];
196 } IVFInfo;
197
198
199 typedef struct
200 {
201 int HiggsFlag;
202 float mass;
203 float pt;
204 float eta;
205 float phi;
206 float dR;
207 float dPhi;
208 float dEta;
209 } HiggsInfo;
210
211 typedef struct
212 {
213 int FatHiggsFlag;
214 float mass;
215 float pt;
216 float eta;
217 float phi;
218 float filteredmass;
219 float filteredpt;
220 float filteredeta;
221 float filteredphi;
222 // float dR;
223 // float dPhi;
224 // float dEta;
225 } FatHiggsInfo;
226
227
228
229 typedef struct
230 {
231 float mass;
232 float pt;
233 float eta;
234 float phi;
235 float status;
236 float charge;
237 float momid;
238 } genParticleInfo;
239
240
241
242 typedef struct
243
244 {
245 float bmass;
246 float bpt;
247 float beta;
248 float bphi;
249 float bstatus;
250 float wdau1mass;
251 float wdau1pt;
252 float wdau1eta;
253 float wdau1phi;
254 float wdau1id;
255 float wdau2mass;
256 float wdau2pt;
257 float wdau2eta;
258 float wdau2phi;
259 float wdau2id;
260
261
262 } genTopInfo;
263
264
265
266
267 typedef struct
268 {
269 bool HiggsCSVtkSharing;
270 bool HiggsIPtkSharing;
271 bool HiggsSVtkSharing;
272 bool FatHiggsCSVtkSharing;
273 bool FatHiggsIPtkSharing;
274 bool FatHiggsSVtkSharing;
275 } TrackSharingInfo;
276
277 typedef struct
278 {
279 float mass; //MT in case of W
280 float pt;
281 float eta;
282 float phi;
283 } TrackInfo;
284
285
286 struct LeptonInfo
287 {
288 void reset()
289 {
290 for(int i =0; i < MAXL;i++){
291 mass[i]=-99; pt[i]=-99; eta[i]=-99; phi[i]=-99; aodCombRelIso[i]=-99; pfCombRelIso[i]=-99; photonIso[i]=-99; neutralHadIso[i]=-99; chargedHadIso[i]=-99; chargedPUIso[i]=-99; particleIso[i]=-99; dxy[i]=-99; dz[i]=-99; type[i]=-99; genPt[i]=-99; genEta[i]=-99; genPhi[i]=-99;
292 id80[i]=-99; id95[i]=-99; vbtf[i]=-99; id80NoIso[i]=-99;
293 charge[i]=-99;
294 }
295 }
296
297 template <class Input> void set(const Input & i, int j,int t)
298 {
299 type[j]=t;
300 pt[j]=i.p4.Pt();
301 mass[j]=i.p4.M();
302 eta[j]=i.p4.Eta();
303 phi[j]=i.p4.Phi();
304 aodCombRelIso[j]=(i.hIso+i.eIso+i.tIso)/i.p4.Pt();
305 pfCombRelIso[j]=(i.pfChaIso+i.pfPhoIso+i.pfNeuIso)/i.p4.Pt();
306 photonIso[j]=i.pfPhoIso;
307 neutralHadIso[j]=i.pfNeuIso;
308 chargedHadIso[j]=i.pfChaIso;
309 chargedPUIso[j]=i.pfChaPUIso;
310 charge[j]=i.charge;
311 if(i.mcFourMomentum.Pt() > 0)
312 {
313 genPt[j]=i.mcFourMomentum.Pt();
314 genEta[j]=i.mcFourMomentum.Eta();
315 genPhi[j]=i.mcFourMomentum.Phi();
316 }
317 setSpecific(i,j);
318 }
319 template <class Input> void setSpecific(const Input & i, int j)
320 {
321 }
322
323
324
325
326 float mass[MAXL]; //MT in case of W
327 float pt[MAXL];
328 float eta[MAXL];
329 float phi[MAXL];
330 float aodCombRelIso[MAXL];
331 float pfCombRelIso[MAXL];
332 float photonIso[MAXL];
333 float neutralHadIso[MAXL];
334 float chargedHadIso[MAXL];
335 float chargedPUIso[MAXL];
336 float particleIso[MAXL];
337 float genPt[MAXL];
338 float genEta[MAXL];
339 float genPhi[MAXL];
340 float dxy[MAXL];
341 float dz[MAXL];
342 int type[MAXL];
343 float id80[MAXL];
344 float id95[MAXL];
345 float vbtf[MAXL];
346 float id80NoIso[MAXL];
347 float charge[MAXL];
348
349 };
350
351 template <> void LeptonInfo::setSpecific<VHbbEvent::ElectronInfo>(const VHbbEvent::ElectronInfo & i, int j){
352 id80[j]=i.id80;
353 id95[j]=i.id95;
354 id80NoIso[j]=(i.innerHits ==0 && !(fabs(i.convDist)<0.02 && fabs(i.convDcot)<0.02) &&
355 ((i.isEB && i.sihih<0.01 && fabs(i.Dphi)<0.06 && fabs(i.Deta)<0.004) || (i.isEE && i.sihih<0.03 && fabs(i.Dphi)<0.03 && fabs(i.Deta)<0.007)));
356 }
357 template <> void LeptonInfo::setSpecific<VHbbEvent::MuonInfo>(const VHbbEvent::MuonInfo & i, int j){
358 dxy[j]=i.ipDb;
359 dz[j]=i.zPVPt;
360 vbtf[j]=( i.globChi2<10 && i.nPixelHits>= 1 && i.globNHits != 0 && i.nHits > 10 && i.cat & 0x1 && i.cat & 0x2 && i.nMatches >=2 && i.ipDb<.2 &&
361 (i.pfChaIso+i.pfPhoIso+i.pfNeuIso)/i.p4.Pt()<.15 && fabs(i.p4.Eta())<2.4 && i.p4.Pt()>20 ) ;
362
363 }
364
365
366 typedef struct
367 {
368 float et;
369 float sumet;
370 float sig;
371 float phi;
372 } METInfo;
373
374
375 typedef struct
376 {
377
378 float et[nMetUnc]; float phi[nMetUnc]; float sumet[nMetUnc];
379 template <class Input> void set(const Input & i, int j)
380 {
381 et[j]=i.p4.Pt();
382 phi[j]=i.p4.Phi();
383 sumet[j]=i.sumEt;
384
385 }
386 } METUncInfo ;
387
388
389 typedef struct
390 {
391 float mht;
392 float ht;
393 float sig;
394 float phi;
395 } MHTInfo;
396
397 typedef struct
398 {
399 float mass;
400 float pt;
401 float wMass;
402 } TopInfo;
403
404 typedef struct
405 {
406 int run;
407 int lumi;
408 int event;
409 int json;
410 } EventInfo;
411
412 typedef struct
413 {
414 void set(const VHbbEvent::SimpleJet & j, int i)
415 {
416 pt[i]=j.p4.Pt();
417 eta[i]=j.p4.Eta();
418 phi[i]=j.p4.Phi();
419 e[i]=j.p4.E();
420 csv[i]=j.csv;
421 csvivf[i]=j.csvivf;
422 cmva[i]=j.cmva;
423 numTracksSV[i] = j.vtxNTracks;
424 vtxMass[i]= j.vtxMass;
425 vtxPt[i]= j.vtxP4.Pt();
426 vtxEta[i]= j.vtxP4.Eta();
427 vtxPhi[i]= j.vtxP4.Phi();
428 vtxE[i]= j.vtxP4.E();
429 vtx3dL[i] = j.vtx3dL;
430 vtx3deL[i] = j.vtx3deL;
431 chf[i]=j.chargedHadronEFraction;
432 nhf[i] =j.neutralHadronEFraction;
433 cef[i] =j.chargedEmEFraction;
434 nef[i] =j.neutralEmEFraction;
435 nconstituents[i] =j.nConstituents;
436 nch[i]=j.ntracks;
437 SF_CSVL[i]=j.SF_CSVL;
438 SF_CSVM[i]=j.SF_CSVM;
439 SF_CSVT[i]=j.SF_CSVT;
440 SF_CSVLerr[i]=j.SF_CSVLerr;
441 SF_CSVMerr[i]=j.SF_CSVMerr;
442 SF_CSVTerr[i]=j.SF_CSVTerr;
443
444 flavour[i]=j.flavour;
445 isSemiLeptMCtruth[i]=j.isSemiLeptMCtruth;
446 isSemiLept[i]=j.isSemiLept;
447 SoftLeptpdgId[i] = j.SoftLeptpdgId;
448 SoftLeptIdlooseMu[i] = j.SoftLeptIdlooseMu;
449 SoftLeptId95[i] = j.SoftLeptId95;
450 SoftLeptPt[i] = j.SoftLeptPt;
451 SoftLeptdR[i] = j.SoftLeptdR;
452 SoftLeptptRel[i] = j.SoftLeptptRel;
453 SoftLeptRelCombIso[i] = j.SoftLeptRelCombIso;
454 if(j.bestMCp4.Pt() > 0)
455 {
456 genPt[i]=j.bestMCp4.Pt();
457 genEta[i]=j.bestMCp4.Eta();
458 genPhi[i]=j.bestMCp4.Phi();
459 }
460 JECUnc[i]=j.jecunc;
461 id[i]=jetId(i);
462 ptRaw[i]=j.ptRaw;
463 ptLeadTrack[i]=j.ptLeadTrack;
464
465 }
466 bool jetId(int i)
467 {
468 if(nhf[i] > 0.99) return false;
469 if(nef[i] > 0.99) return false;
470 if(nconstituents[i] <= 1) return false;
471 if(fabs(eta[i])<2.5) {
472 if(cef[i] > 0.99) return false;
473 if(chf[i] == 0) return false;
474 if(nch[i]== 0) return false;
475 }
476 return true;
477 }
478 void reset()
479 {
480 for(int i=0;i<MAXJ;i++) {
481 pt[i]=-99; eta[i]=-99; phi[i]=-99;e[i]=-99;csv[i]=-99;
482 csvivf[i]=-99; cmva[i]=-99;
483 cosTheta[i]=-99; numTracksSV[i]=-99; chf[i]=-99; nhf[i]=-99; cef[i]=-99; nef[i]=-99; nch[i]=-99; nconstituents[i]=-99; flavour[i]=-99; isSemiLeptMCtruth[i]=-99; isSemiLept[i]=-99;
484 SoftLeptpdgId[i] = -99; SoftLeptIdlooseMu[i] = -99; SoftLeptId95[i] = -99; SoftLeptPt[i] = -99; SoftLeptdR[i] = -99; SoftLeptptRel[i] = -99; SoftLeptRelCombIso[i] = -99;
485 genPt[i]=-99; genEta[i]=-99; genPhi[i]=-99; JECUnc[i]=-99; ptRaw[i]=-99.; ptLeadTrack[i]=-99.;
486 }
487 }
488 float pt[MAXJ];
489 float eta[MAXJ];
490 float phi[MAXJ];
491 float e[MAXJ];
492 float csv[MAXJ];
493 float csvivf[MAXJ];
494 float cmva[MAXJ];
495 float cosTheta[MAXJ];
496 int numTracksSV[MAXJ];
497 float chf[MAXJ];
498 float nhf[MAXJ];
499 float cef[MAXJ];
500 float nef[MAXJ];
501 float nch[MAXJ];
502 float nconstituents[MAXJ];
503 float flavour[MAXJ];
504 int isSemiLept[MAXJ];
505 int isSemiLeptMCtruth[MAXJ];
506 int SoftLeptpdgId[MAXJ] ;
507 int SoftLeptIdlooseMu[MAXJ] ;
508 int SoftLeptId95[MAXJ] ;
509 float SoftLeptPt[MAXJ] ;
510 float SoftLeptdR[MAXJ] ;
511 float SoftLeptptRel[MAXJ] ;
512 float SoftLeptRelCombIso[MAXJ];
513 float genPt[MAXJ];
514 float genEta[MAXJ];
515 float genPhi[MAXJ];
516 float JECUnc[MAXJ];
517 float vtxMass[MAXJ];
518 float vtxPt[MAXJ];
519 float vtxEta[MAXJ];
520 float vtxPhi[MAXJ];
521 float vtxE[MAXJ];
522 float vtx3dL [MAXJ];
523 float vtx3deL[MAXJ];
524 bool id[MAXJ];
525 float SF_CSVL[MAXJ];
526 float SF_CSVM[MAXJ];
527 float SF_CSVT[MAXJ];
528 float SF_CSVLerr[MAXJ];
529 float SF_CSVMerr[MAXJ];
530 float SF_CSVTerr[MAXJ];
531 float ptRaw[MAXJ];
532 float ptLeadTrack[MAXJ];
533 } JetInfo;
534
535 int main(int argc, char* argv[])
536 {
537 gROOT->Reset();
538
539
540
541 TTree *_outTree;
542 IVFInfo IVF;
543 SimBHadronInfo SimBs;
544 float rho,rho25;
545 int nPVs;
546 METInfo MET;
547 METInfo fakeMET;
548 METInfo METnoPU;
549 METInfo METnoPUCh;
550 METInfo METtype1corr;
551 METInfo METtype1p2corr;
552 METInfo METnoPUtype1corr;
553 METInfo METnoPUtype1p2corr;
554
555 MHTInfo MHT;
556
557 METUncInfo metUnc;
558
559
560 TopInfo top;
561 EventInfo EVENT;
562 // JetInfo jet1,jet2, addJet1, addJet2;
563 // lepton1,lepton2;
564 JetInfo hJets, aJets, fathFilterJets, aJetsFat;
565 LeptonInfo vLeptons, aLeptons;
566 int naJets=0, nhJets=0, nfathFilterJets=0, naJetsFat=0;
567 HiggsInfo H,SVH,SimBsH;
568 FatHiggsInfo FatH;
569 genParticleInfo genZ, genZstar, genWstar, genW, genH, genB, genBbar; //add here the fatjet higgs
570 genTopInfo genTop, genTbar;
571 TrackInfo V;
572 int nvlep=0,nalep=0;
573 float lheV_pt=0; //for the Madgraph sample stitching
574 float PDFweight=1.; // for pdf reweighting (only madgraph)
575 TrackSharingInfo TkSharing; // track sharing info;
576
577 float HVdPhi,HVMass,HMETdPhi,VMt,deltaPullAngle,deltaPullAngleAK7,deltaPullAngle2,deltaPullAngle2AK7,gendrcc,gendrbb, genZpt, genWpt, genHpt, weightTrig, weightTrigMay,weightTrigV4, weightTrigMET, weightTrigOrMu30, minDeltaPhijetMET, jetPt_minDeltaPhijetMET , PUweight, PUweight2011B;
578 float PU0,PUp1,PUm1;
579
580 float weightEleRecoAndId,weightEleTrigJetMETPart, weightEleTrigElePart,weightEleTrigEleAugPart;
581 float weightTrigMET80, weightTrigMET100, weightTrig2CJet20 , weightTrigMET150 , weightTrigMET802CJet, weightTrigMET1002CJet, weightTrigMETLP ;
582
583 int WplusMode,WminusMode;
584 int Vtype,nSvs=0,nSimBs=0,numJets,numBJets,eventFlav;
585 // bool isMET80_CJ80, ispfMHT150, isMET80_2CJ20,isMET65_2CJ20, isJETID,isIsoMu17;
586 bool triggerFlags[500],hbhe,ecalFlag,totalKinematics, cschaloFlag, hcallaserFlag, trackingfailureFlag ;
587
588
589 float btag1T2CSF=1.,btag2TSF=1.,btag1TSF=1.,btagA0CSF=1., btagA0TSF=1., btag2CSF=1., btag1TA1C=1.;
590 // ----------------------------------------------------------------------
591 // First Part:
592 //
593 // * enable the AutoLibraryLoader
594 // * book the histograms of interest
595 // * open the input file
596 // ----------------------------------------------------------------------
597
598 // load framework libraries
599 gSystem->Load("libFWCoreFWLite");
600 gSystem->Load("libDataFormatsFWLite");
601 AutoLibraryLoader::enable();
602
603 //init the PDF set. Need to check for each sample beacuse pdfset and pdfmember info missing in AODSIM. (for Madgraph always (1,10042,0) : mctqe6l)
604 // LHAPDF::initPDFSet(1,_pdfset, _pdfmember);
605 LHAPDF::initPDFSet(1, 10042, 0);
606
607 // parse arguments
608 if ( argc < 2 ) {
609 return 0;
610 }
611
612 std::vector<VHbbCandidate> * candZlocal = new std::vector<VHbbCandidate>;
613 std::vector<VHbbCandidate> * candWlocal = new std::vector<VHbbCandidate>;
614
615 // get the python configuration
616 PythonProcessDesc builder(argv[1]);
617 const edm::ParameterSet& in = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("fwliteInput" );
618 const edm::ParameterSet& out = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("fwliteOutput");
619 const edm::ParameterSet& ana = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("Analyzer");
620
621 std::vector<edm::LuminosityBlockRange> jsonVector;
622 if ( in.exists("lumisToProcess") )
623 {
624 std::vector<edm::LuminosityBlockRange> const & lumisTemp =
625 in.getUntrackedParameter<std::vector<edm::LuminosityBlockRange> > ("lumisToProcess");
626 jsonVector.resize( lumisTemp.size() );
627 copy( lumisTemp.begin(), lumisTemp.end(), jsonVector.begin() );
628 }
629
630 // now get each parameter
631 int maxEvents_( in.getParameter<int>("maxEvents") );
632 int skipEvents_( in.getParameter<int>("skipEvents") );
633 int runMin_( in.getParameter<int>("runMin") );
634 int runMax_( in.getParameter<int>("runMax") );
635 unsigned int outputEvery_( in.getParameter<unsigned int>("outputEvery") );
636 std::string outputFile_( out.getParameter<std::string>("fileName" ) );
637 std::vector<std::string> triggers( ana.getParameter<std::vector<std::string> >("triggers") );
638 double btagThr = ana.getParameter<double>("bJetCountThreshold" );
639 bool fromCandidate = ana.getParameter<bool>("readFromCandidates");
640 bool useHighestPtHiggsZ = ana.getParameter<bool>("useHighestPtHiggsZ");
641 bool useHighestPtHiggsW = ana.getParameter<bool>("useHighestPtHiggsW");
642 HbbCandidateFinderAlgo * algoZ = new HbbCandidateFinderAlgo(ana.getParameter<bool>("verbose"), ana.getParameter<double>("jetPtThresholdZ"),useHighestPtHiggsZ);
643 HbbCandidateFinderAlgo * algoW = new HbbCandidateFinderAlgo(ana.getParameter<bool>("verbose"), ana.getParameter<double>("jetPtThresholdW"),useHighestPtHiggsW );
644 HbbCandidateFinderAlgo * algoRecoverLowPt = new HbbCandidateFinderAlgo(ana.getParameter<bool>("verbose"), 15, true);
645
646 TriggerWeight triggerWeight(ana);
647 BTagWeight btag(2); // 2 operating points "Custom" = 0.5 and "Tight = 0.898"
648 BTagSampleEfficiency btagEff( ana.getParameter<std::string>("btagEffFileName" ).c_str() );
649
650 std::vector<std::string> inputFiles_( in.getParameter<std::vector<std::string> >("fileNames") );
651 // std::string inputFile( in.getParameter<std::string> ("fileName") );
652
653 std::string PUmcfileName_ = in.getParameter<std::string> ("PUmcfileName") ;
654 std::string PUmcfileName2011B_ = in.getParameter<std::string> ("PUmcfileName2011B") ;
655
656 std::string PUdatafileName_ = in.getParameter<std::string> ("PUdatafileName") ;
657 std::string PUdatafileName2011B_ = in.getParameter<std::string> ("PUdatafileName2011B") ;
658 std::string Weight3DfileName_ = in.getParameter<std::string> ("Weight3DfileName") ;
659
660
661
662 bool isMC_( ana.getParameter<bool>("isMC") );
663 TriggerReader trigger(isMC_);
664 TriggerReader patFilters(false);
665
666 edm::LumiReWeighting lumiWeights;
667 edm::Lumi3DReWeighting lumiWeights2011B;
668 if(isMC_)
669 {
670 lumiWeights = edm::LumiReWeighting(PUmcfileName_,PUdatafileName_ , "pileup", "pileup");
671
672 lumiWeights2011B = edm::Lumi3DReWeighting(PUmcfileName2011B_,PUdatafileName2011B_ , "pileup", "pileup");
673 if(Weight3DfileName_!="")
674 { lumiWeights2011B.weight3D_init(Weight3DfileName_.c_str()); }
675 else
676 {
677 lumiWeights2011B.weight3D_init(1.0); // generate the weights the fisrt time;
678 }
679
680 }
681
682 // TFile *_outPUFile = new TFile((outputFile_+"_PU").c_str(), "recreate");
683 TFile *_outFile = new TFile(outputFile_.c_str(), "recreate");
684 TH1F * count = new TH1F("Count","Count", 1,0,2 );
685 TH1F * countWithPU = new TH1F("CountWithPU","CountWithPU", 1,0,2 );
686 TH1F * countWithPU2011B = new TH1F("CountWithPU2011B","CountWithPU2011B", 1,0,2 );
687 TH3F * input3DPU = new TH3F("Input3DPU","Input3DPU", 36,-0.5,35.5,36,-0.5,35.5, 36,-0.5,35.5 );
688
689 TH1F * pu = new TH1F("pileup","",51,-0.5,50.5);
690 _outTree = new TTree("tree", "myTree");
691
692 _outTree->Branch("H" , &H , "HiggsFlag/I:mass/F:pt/F:eta:phi/F:dR/F:dPhi/F:dEta/F");
693 _outTree->Branch("V" , &V , "mass/F:pt/F:eta:phi/F");
694 _outTree->Branch("FatH" , &FatH , "FatHiggsFlag/I:mass/F:pt/F:eta:phi/F:filteredmass/F:filteredpt/F:filteredeta/F:filteredphi/F");
695 _outTree->Branch("lheV_pt" , &lheV_pt , "lheV_pt/F");
696 _outTree->Branch("PDFweight" , &PDFweight , "PDFweight/F");
697 _outTree->Branch("genZ" , &genZ , "mass/F:pt/F:eta:phi/F:status/F:charge:momid/F");
698 _outTree->Branch("genZstar" , &genZstar , "mass/F:pt/F:eta:phi/F:status/F:charge:momid/F");
699 _outTree->Branch("genW" , &genW , "mass/F:pt/F:eta:phi/F:status/F:charge:momid/F");
700 _outTree->Branch("genWstar" , &genWstar , "mass/F:pt/F:eta:phi/F:status/F:charge:momid/F");
701 _outTree->Branch("genH" , &genH , "mass/F:pt/F:eta:phi/F:status/F:charge:momid/F");
702 _outTree->Branch("genB" , &genB , "mass/F:pt/F:eta:phi/F:status/F:charge:momid/F");
703 _outTree->Branch("genBbar" , &genBbar , "mass/F:pt/F:eta:phi/F:status/F:charge:momid/F");
704
705 _outTree->Branch("genTop" , &genTop , "bmass/F:bpt/F:beta:bphi/F:bstatus/F:wdau1mass/F:wdau1pt/F:wdau1eta:wdau1phi/F:wdau1id/F:wdau2mass/F:wdau2pt/F:wdau2eta:wdau2phi/F:wdau2id/F");
706 _outTree->Branch("genTbar" , &genTbar , "bmass/F:bpt/F:beta:bphi/F:bstatus/F:wdau1mass/F:wdau1pt/F:wdau1eta:wdau1phi/F:wdau1id/F:wdau2mass/F:wdau2pt/F:wdau2eta:wdau2phi/F:wdau2id/F");
707
708 _outTree->Branch("TkSharing", &TkSharing, "HiggsCSVtkSharing/b:HiggsIPtkSharing:HiggsSVtkSharing:FatHiggsCSVtkSharing:FatHiggsIPtkSharing:FatHiggsSVtkSharing");
709
710 _outTree->Branch("nhJets" , &nhJets , "nhJets/I");
711 _outTree->Branch("nfathFilterJets", &nfathFilterJets, "nfathFilterJets/I");
712 _outTree->Branch("naJets" , &naJets , "naJets/I");
713
714
715
716 _outTree->Branch("hJet_pt",hJets.pt ,"pt[nhJets]/F");
717 _outTree->Branch("hJet_eta",hJets.eta ,"eta[nhJets]/F");
718 _outTree->Branch("hJet_phi",hJets.phi ,"phi[nhJets]/F");
719 _outTree->Branch("hJet_e",hJets.e ,"e[nhJets]/F");
720 _outTree->Branch("hJet_csv",hJets.csv ,"csv[nhJets]/F");
721 _outTree->Branch("hJet_csvivf",hJets.csvivf ,"csvivf[nhJets]/F");
722 _outTree->Branch("hJet_cmva",hJets.cmva ,"cmva[nhJets]/F");
723 _outTree->Branch("hJet_cosTheta",hJets.cosTheta ,"cosTheta[nhJets]/F");
724 _outTree->Branch("hJet_numTracksSV",hJets.numTracksSV ,"numTracksSV[nhJets]/I");
725 _outTree->Branch("hJet_chf",hJets.chf ,"chf[nhJets]/F");
726 _outTree->Branch("hJet_nhf",hJets.nhf ,"nhf[nhJets]/F");
727 _outTree->Branch("hJet_cef",hJets.cef ,"cef[nhJets]/F");
728 _outTree->Branch("hJet_nef",hJets.nef ,"nef[nhJets]/F");
729 _outTree->Branch("hJet_nch",hJets.nch ,"nch[nhJets]/F");
730 _outTree->Branch("hJet_nconstituents",hJets.nconstituents ,"nconstituents[nhJets]");
731 _outTree->Branch("hJet_flavour",hJets.flavour ,"flavour[nhJets]/F");
732 _outTree->Branch("hJet_isSemiLept",hJets.isSemiLept ,"isSemiLept[nhJets]/I");
733 _outTree->Branch("hJet_isSemiLeptMCtruth",hJets.isSemiLeptMCtruth ,"isSemiLeptMCtruth[nhJets]/I");
734 _outTree->Branch("hJet_SoftLeptpdgId", hJets.SoftLeptpdgId , "SoftLeptpdgId[nhJets]/I");
735 _outTree->Branch("hJet_SoftLeptIdlooseMu", hJets.SoftLeptIdlooseMu , "SoftLeptIdlooseMu[nhJets]/I");
736 _outTree->Branch("hJet_SoftLeptId95", hJets.SoftLeptId95 , "SoftLeptId95[nhJets]/I");
737 _outTree->Branch("hJet_SoftLeptPt", hJets.SoftLeptPt , "SoftLeptPt[nhJets]/F");
738 _outTree->Branch("hJet_SoftLeptdR", hJets.SoftLeptdR , "SoftLeptdR[nhJets]/F");
739 _outTree->Branch("hJet_SoftLeptptRel", hJets.SoftLeptptRel , "SoftLeptptRel[nhJets]/F");
740 _outTree->Branch("hJet_SoftLeptRelCombIso", hJets.SoftLeptRelCombIso , "SoftLeptRelCombIso[nhJets]/F");
741 _outTree->Branch("hJet_genPt",hJets.genPt ,"genPt[nhJets]/F");
742 _outTree->Branch("hJet_genEta",hJets.genEta ,"genEta[nhJets]/F");
743 _outTree->Branch("hJet_genPhi",hJets.genPhi ,"genPhi[nhJets]/F");
744 _outTree->Branch("hJet_JECUnc",hJets.JECUnc ,"JECUnc[nhJets]/F");
745 _outTree->Branch("hJet_vtxMass",hJets.vtxMass ,"vtxMass[nhJets]/F");
746 _outTree->Branch("hJet_vtxPt",hJets.vtxPt ,"vtxPt[nhJets]/F");
747 _outTree->Branch("hJet_vtxEta",hJets.vtxEta ,"vtxEta[nhJets]/F");
748 _outTree->Branch("hJet_vtxPhi",hJets.vtxPhi ,"vtxPhi[nhJets]/F");
749 _outTree->Branch("hJet_vtxE",hJets.vtxE ,"vtxE[nhJets]/F");
750 _outTree->Branch("hJet_vtx3dL",hJets.vtx3dL ,"vtx3dL[nhJets]/F");
751 _outTree->Branch("hJet_vtx3deL",hJets.vtx3deL ,"vtx3deL[nhJets]/F");
752 _outTree->Branch("hJet_id",hJets.id ,"id[nhJets]/b");
753 _outTree->Branch("hJet_SF_CSVL",hJets.SF_CSVL ,"SF_CSVL[nhJets]/b");
754 _outTree->Branch("hJet_SF_CSVM",hJets.SF_CSVM ,"SF_CSVM[nhJets]/b");
755 _outTree->Branch("hJet_SF_CSVT",hJets.SF_CSVT ,"SF_CSVT[nhJets]/b");
756 _outTree->Branch("hJet_SF_CSVLerr",hJets.SF_CSVLerr ,"SF_CSVLerr[nhJets]/b");
757 _outTree->Branch("hJet_SF_CSVMerr",hJets.SF_CSVMerr ,"SF_CSVMerr[nhJets]/b");
758 _outTree->Branch("hJet_SF_CSVTerr",hJets.SF_CSVTerr ,"SF_CSVTerr[nhJets]/b");
759 _outTree->Branch("hJet_ptRaw",hJets.ptRaw ,"ptRaw[nhJets]/F");
760 _outTree->Branch("hJet_ptLeadTrack",hJets.ptLeadTrack ,"ptLeadTrack[nhJets]/F");
761
762 _outTree->Branch("fathFilterJets_pt",fathFilterJets.pt ,"pt[nfathFilterJets]/F");
763 _outTree->Branch("fathFilterJets_eta",fathFilterJets.eta ,"eta[nfathFilterJets]/F");
764 _outTree->Branch("fathFilterJets_phi",fathFilterJets.phi ,"phi[nfathFilterJets]/F");
765 _outTree->Branch("fathFilterJets_e",fathFilterJets.e ,"e[nfathFilterJets]/F");
766 _outTree->Branch("fathFilterJets_csv",fathFilterJets.csv ,"csv[nfathFilterJets]/F");
767 _outTree->Branch("fathFilterJets_chf",fathFilterJets.chf ,"chf[nfathFilterJets]/F");
768 _outTree->Branch("fathFilterJets_ptRaw",fathFilterJets.ptRaw ,"ptRaw[nfathFilterJets]/F");
769 _outTree->Branch("fathFilterJets_ptLeadTrack",fathFilterJets.ptLeadTrack ,"ptLeadTrack[nfathFilterJets]/F");
770 _outTree->Branch("fathFilterJets_flavour",fathFilterJets.flavour ,"flavour[nfathFilterJets]/F");
771 _outTree->Branch("fathFilterJets_isSemiLept",fathFilterJets.isSemiLept ,"isSemiLept[nfathFilterJets]/F");
772 _outTree->Branch("fathFilterJets_isSemiLeptMCtruth",fathFilterJets.isSemiLeptMCtruth ,"isSemiLeptMCtruth[nfathFilterJets]/F");
773 _outTree->Branch("fathFilterJets_genPt",fathFilterJets.genPt ,"genPt[nfathFilterJets]/F");
774 _outTree->Branch("fathFilterJets_genEta",fathFilterJets.genEta ,"genEta[nfathFilterJets]/F");
775 _outTree->Branch("fathFilterJets_genPhi",fathFilterJets.genPhi ,"genPhi[nfathFilterJets]/F");
776 _outTree->Branch("fathFilterJets_vtxMass",fathFilterJets.vtxMass ,"vtxMass[nfathFilterJets]/F");
777 _outTree->Branch("fathFilterJets_vtx3dL",fathFilterJets.vtx3dL ,"vtx3dL[nfathFilterJets]/F");
778 _outTree->Branch("fathFilterJets_vtx3deL",fathFilterJets.vtx3deL ,"vtx3deL[nfathFilterJets]/F");
779 _outTree->Branch("fathFilterJets_vtxPt",fathFilterJets.vtxPt ,"vtxPt[nfathFilterJets]/F");
780 _outTree->Branch("fathFilterJets_vtxEta",fathFilterJets.vtxEta ,"vtxEta[nfathFilterJets]/F");
781 _outTree->Branch("fathFilterJets_vtxPhi",fathFilterJets.vtxPhi ,"vtxPhi[nfathFilterJets]/F");
782 _outTree->Branch("fathFilterJets_vtxE",fathFilterJets.vtxE ,"vtxE[nfathFilterJets]/F");
783
784
785 // _outTree->Branch("fathFilterJets_pt",fathFilterJets.pt ,"pt[nfathFilterJets]/F");
786 // _outTree->Branch("fathFilterJets_eta",fathFilterJets.eta ,"eta[nfathFilterJets]/F");
787 // _outTree->Branch("fathFilterJets_phi",fathFilterJets.phi ,"phi[nfathFilterJets]/F");
788 // _outTree->Branch("fathFilterJets_e",fathFilterJets.e ,"e[nfathFilterJets]/F");
789 // _outTree->Branch("fathFilterJets_csv",fathFilterJets.csv ,"csv[nfathFilterJets]/F");
790 _outTree->Branch("fathFilterJets_csvivf",fathFilterJets.csvivf ,"csvivf[nfathFilterJets]/F");
791 _outTree->Branch("fathFilterJets_cmva",fathFilterJets.cmva ,"cmva[nfathFilterJets]/F");
792 // _outTree->Branch("fathFilterJets_flavour",fathFilterJets.flavour ,"flavour[nfathFilterJets]/F");
793
794
795 _outTree->Branch("aJet_pt",aJets.pt ,"pt[naJets]/F");
796 _outTree->Branch("aJet_eta",aJets.eta ,"eta[naJets]/F");
797 _outTree->Branch("aJet_phi",aJets.phi ,"phi[naJets]/F");
798 _outTree->Branch("aJet_e",aJets.e ,"e[naJets]/F");
799 _outTree->Branch("aJet_csv",aJets.csv ,"csv[naJets]/F");
800 _outTree->Branch("aJet_csvivf",aJets.csvivf ,"csvivf[naJets]/F");
801 _outTree->Branch("aJet_cmva",aJets.cmva ,"cmva[naJets]/F");
802 _outTree->Branch("aJet_cosTheta",aJets.cosTheta ,"cosTheta[naJets]/F");
803 _outTree->Branch("aJet_numTracksSV",aJets.numTracksSV ,"numTracksSV[naJets]/I");
804 _outTree->Branch("aJet_chf",aJets.chf ,"chf[naJets]/F");
805 _outTree->Branch("aJet_nhf",aJets.nhf ,"nhf[naJets]/F");
806 _outTree->Branch("aJet_cef",aJets.cef ,"cef[naJets]/F");
807 _outTree->Branch("aJet_nef",aJets.nef ,"nef[naJets]/F");
808 _outTree->Branch("aJet_nch",aJets.nch ,"nch[naJets]/F");
809 _outTree->Branch("aJet_nconstituents",aJets.nconstituents ,"nconstituents[naJets]");
810 _outTree->Branch("aJet_flavour",aJets.flavour ,"flavour[naJets]/F");
811 _outTree->Branch("aJet_isSemiLept",aJets.isSemiLept ,"isSemiLept[naJets]/I");
812 _outTree->Branch("aJet_isSemiLeptMCtruth",aJets.isSemiLeptMCtruth ,"isSemiLeptMCtruth[naJets]/I");
813 _outTree->Branch("aJet_SoftLeptpdgId",aJets.SoftLeptpdgId , "SoftLeptpdgId[naJets]/I");
814 _outTree->Branch("aJet_SoftLeptIdlooseMu", aJets.SoftLeptIdlooseMu , "SoftLeptIdlooseMu[naJets]/I");
815 _outTree->Branch("aJet_SoftLeptId95", aJets.SoftLeptId95 , "SoftLeptId95[naJets]/I");
816 _outTree->Branch("aJet_SoftLeptPt", aJets.SoftLeptPt , "SoftLeptPt[naJets]/F");
817 _outTree->Branch("aJet_SoftLeptdR", aJets.SoftLeptdR , "SoftLeptdR[naJets]/F");
818 _outTree->Branch("aJet_SoftLeptptRel", aJets.SoftLeptptRel , "SoftLeptptRel[naJets]/F");
819 _outTree->Branch("aJet_SoftLeptRelCombIso", aJets.SoftLeptRelCombIso , "SoftLeptRelCombIso[naJets]/F");
820
821 _outTree->Branch("aJet_genPt",aJets.genPt ,"genPt[naJets]/F");
822 _outTree->Branch("aJet_genEta",aJets.genEta ,"genEta[naJets]/F");
823 _outTree->Branch("aJet_genPhi",aJets.genPhi ,"genPhi[naJets]/F");
824 _outTree->Branch("aJet_JECUnc",aJets.JECUnc ,"JECUnc[naJets]/F");
825 _outTree->Branch("aJet_vtxMass",aJets.vtxMass ,"vtxMass[naJets]/F");
826 _outTree->Branch("aJet_vtx3dL",aJets.vtx3dL ,"vtx3dL[naJets]/F");
827 _outTree->Branch("aJet_vtx3deL",aJets.vtx3deL ,"vtx3deL[naJets]/F");
828 _outTree->Branch("aJet_id",aJets.id ,"id[naJets]/b");
829 _outTree->Branch("aJet_SF_CSVL",aJets.SF_CSVL ,"SF_CSVL[naJets]/b");
830 _outTree->Branch("aJet_SF_CSVM",aJets.SF_CSVM ,"SF_CSVM[naJets]/b");
831 _outTree->Branch("aJet_SF_CSVT",aJets.SF_CSVT ,"SF_CSVT[naJets]/b");
832 _outTree->Branch("aJet_SF_CSVLerr",aJets.SF_CSVLerr ,"SF_CSVLerr[naJets]/b");
833 _outTree->Branch("aJet_SF_CSVMerr",aJets.SF_CSVMerr ,"SF_CSVMerr[naJets]/b");
834 _outTree->Branch("aJet_SF_CSVTerr",aJets.SF_CSVTerr ,"SF_CSVTerr[naJets]/b");
835 _outTree->Branch("aJet_ptRaw",aJets.ptRaw ,"ptRaw[naJets]/F");
836
837 _outTree->Branch("naJetsFat" , &naJetsFat , "naJetsFat/I");
838 _outTree->Branch("aJetFat_pt",aJetsFat.pt ,"pt[naJetsFat]/F");
839 _outTree->Branch("aJetFat_eta",aJetsFat.eta ,"eta[naJetsFat]/F");
840 _outTree->Branch("aJetFat_phi",aJetsFat.phi ,"phi[naJetsFat]/F");
841 _outTree->Branch("aJetFat_e",aJetsFat.e ,"e[naJetsFat]/F");
842 _outTree->Branch("aJetFat_csv",aJetsFat.csv ,"csv[naJetsFat]/F");
843
844
845 _outTree->Branch("numJets" , &numJets , "numJets/I" );
846 _outTree->Branch("numBJets" , &numBJets , "numBJets/I" );
847 _outTree->Branch("deltaPullAngle", &deltaPullAngle , "deltaPullAngle/F");
848 _outTree->Branch("deltaPullAngle2", &deltaPullAngle2 , "deltaPullAngle2/F");
849 _outTree->Branch("gendrcc" , &gendrcc , "gendrcc/F");
850 _outTree->Branch("gendrbb" , &gendrbb , "gendrbb/F");
851 _outTree->Branch("genZpt" , &genZpt , "genZpt/F");
852 _outTree->Branch("genWpt" , &genWpt , "genWpt/F");
853 _outTree->Branch("genHpt" , &genHpt , "genHpt/F");
854 _outTree->Branch("weightTrig" , &weightTrig , "weightTrig/F");
855 _outTree->Branch("weightTrigMay" , &weightTrigMay , "weightTrigMay/F");
856 _outTree->Branch("weightTrigV4" , &weightTrigV4 , "weightTrigV4/F");
857 _outTree->Branch("weightTrigMET" , &weightTrigMET , "weightTrigMET/F");
858 _outTree->Branch("weightTrigOrMu30" , &weightTrigOrMu30 , "weightTrigOrMu30/F");
859 _outTree->Branch("weightEleRecoAndId" , &weightEleRecoAndId , "weightEleRecoAndId/F");
860 _outTree->Branch("weightEleTrigJetMETPart" , &weightEleTrigJetMETPart , "weightEleTrigJetMETPart/F");
861 _outTree->Branch("weightEleTrigElePart" , &weightEleTrigElePart , "weightEleTrigElePart/F");
862 _outTree->Branch("weightEleTrigEleAugPart" , &weightEleTrigEleAugPart , "weightEleTrigEleAugPart/F");
863
864 _outTree->Branch("weightTrigMET80" , &weightTrigMET80 , "weightTrigMET80/F");
865 _outTree->Branch("weightTrigMET100" , &weightTrigMET100 , "weightTrigMET100/F");
866 _outTree->Branch("weightTrig2CJet20" , &weightTrig2CJet20 , "weightTrig2CJet20/F");
867 _outTree->Branch("weightTrigMET150" , &weightTrigMET150 , "weightTrigMET150/F");
868 _outTree->Branch("weightTrigMET802CJet" , &weightTrigMET802CJet , "weightTrigMET802CJet/F");
869 _outTree->Branch("weightTrigMET1002CJet" , &weightTrigMET1002CJet , "weightTrigMET1002CJet/F");
870 _outTree->Branch("weightTrigMETLP" , &weightTrigMETLP , "weightTrigMETLP/F");
871
872
873 _outTree->Branch("deltaPullAngleAK7", &deltaPullAngleAK7 , "deltaPullAngleAK7/F");
874 _outTree->Branch("deltaPullAngle2AK7", &deltaPullAngle2AK7 , "deltaPullAngle2AK7/F");
875 _outTree->Branch("PU0", &PU0 , "PU0/F");
876 _outTree->Branch("PUm1", &PUm1 , "PUm1/F");
877 _outTree->Branch("PUp1", &PUp1 , "PUp1/F");
878 _outTree->Branch("PUweight", &PUweight , "PUweight/F");
879 _outTree->Branch("PUweight2011B", &PUweight2011B , "PUweight2011B/F");
880 _outTree->Branch("eventFlav", &eventFlav , "eventFlav/I");
881
882
883
884
885 _outTree->Branch("Vtype" , &Vtype , "Vtype/I" );
886 _outTree->Branch("HVdPhi" , &HVdPhi , "HVdPhi/F" );
887 _outTree->Branch("HVMass" , &HVMass , "HVMass/F" );
888 _outTree->Branch("HMETdPhi" , &HMETdPhi , "HMETdPhi/F" );
889 _outTree->Branch("VMt" , &VMt , "VMt/F" );
890
891 _outTree->Branch("nvlep" , &nvlep , "nvlep/I");
892 _outTree->Branch("nalep" , &nalep , "nalep/I");
893
894 _outTree->Branch("vLepton_mass",vLeptons.mass ,"mass[nvlep]/F");
895 _outTree->Branch("vLepton_pt",vLeptons.pt ,"pt[nvlep]/F");
896 _outTree->Branch("vLepton_eta",vLeptons.eta ,"eta[nvlep]");
897 _outTree->Branch("vLepton_phi",vLeptons.phi ,"phi[nvlep]/F");
898 _outTree->Branch("vLepton_aodCombRelIso",vLeptons.aodCombRelIso ,"aodCombRelIso[nvlep]/F");
899 _outTree->Branch("vLepton_pfCombRelIso",vLeptons.pfCombRelIso ,"pfCombRelIso[nvlep]/F");
900 _outTree->Branch("vLepton_photonIso",vLeptons.photonIso ,"photonIso[nvlep]/F");
901 _outTree->Branch("vLepton_neutralHadIso",vLeptons.neutralHadIso ,"neutralHadIso[nvlep]/F");
902 _outTree->Branch("vLepton_chargedHadIso",vLeptons.chargedHadIso ,"chargedHadIso[nvlep]/F");
903 _outTree->Branch("vLepton_chargedPUIso",vLeptons.chargedPUIso ,"chargedPUIso[nvlep]/F");
904 _outTree->Branch("vLepton_particleIso",vLeptons.particleIso ,"particleIso[nvlep]/F");
905 _outTree->Branch("vLepton_dxy",vLeptons.dxy ,"dxy[nvlep]/F");
906 _outTree->Branch("vLepton_dz",vLeptons.dz ,"dz[nvlep]/F");
907 _outTree->Branch("vLepton_type",vLeptons.type ,"type[nvlep]/I");
908 _outTree->Branch("vLepton_id80",vLeptons.id80 ,"id80[nvlep]/F");
909 _outTree->Branch("vLepton_id95",vLeptons.id95 ,"id95[nvlep]/F");
910 _outTree->Branch("vLepton_vbtf",vLeptons.vbtf ,"vbtf[nvlep]/F");
911 _outTree->Branch("vLepton_id80NoIso",vLeptons.id80NoIso ,"id80NoIso[nvlep]/F");
912 _outTree->Branch("vLepton_genPt",vLeptons.genPt ,"genPt[nvlep]/F");
913 _outTree->Branch("vLepton_genEta",vLeptons.genEta ,"genEta[nvlep]");
914 _outTree->Branch("vLepton_genPhi",vLeptons.genPhi ,"genPhi[nvlep]/F");
915 _outTree->Branch("vLepton_charge",vLeptons.charge ,"charge[nvlep]/F");
916
917 _outTree->Branch("aLepton_mass",aLeptons.mass ,"mass[nalep]/F");
918 _outTree->Branch("aLepton_pt",aLeptons.pt ,"pt[nalep]/F");
919 _outTree->Branch("aLepton_eta",aLeptons.eta ,"eta[nalep]");
920 _outTree->Branch("aLepton_phi",aLeptons.phi ,"phi[nalep]/F");
921 _outTree->Branch("aLepton_aodCombRelIso",aLeptons.aodCombRelIso ,"aodCombRelIso[nalep]/F");
922 _outTree->Branch("aLepton_pfCombRelIso",aLeptons.pfCombRelIso ,"pfCombRelIso[nalep]/F");
923 _outTree->Branch("aLepton_photonIso",aLeptons.photonIso ,"photonIso[nalep]/F");
924 _outTree->Branch("aLepton_neutralHadIso",aLeptons.neutralHadIso ,"neutralHadIso[nalep]/F");
925 _outTree->Branch("aLepton_chargedHadIso",aLeptons.chargedHadIso ,"chargedHadIso[nalep]/F");
926 _outTree->Branch("aLepton_chargedPUIso",aLeptons.chargedPUIso ,"chargedPUIso[nalep]/F");
927 _outTree->Branch("aLepton_particleIso",aLeptons.particleIso ,"particleIso[nalep]/F");
928 _outTree->Branch("aLepton_dxy",aLeptons.dxy ,"dxy[nalep]/F");
929 _outTree->Branch("aLepton_dz",aLeptons.dz ,"dz[nalep]/F");
930 _outTree->Branch("aLepton_type",aLeptons.type ,"type[nalep]/I");
931 _outTree->Branch("aLepton_id80",aLeptons.id80 ,"id80[nalep]/F");
932 _outTree->Branch("aLepton_id95",aLeptons.id95 ,"id95[nalep]/F");
933 _outTree->Branch("aLepton_vbtf",aLeptons.vbtf ,"vbtf[nalep]/F");
934 _outTree->Branch("aLepton_id80NoIso",aLeptons.id80NoIso ,"id80NoIso[nalep]/F");
935 _outTree->Branch("aLepton_genPt",aLeptons.genPt ,"genPt[nalep]/F");
936 _outTree->Branch("aLepton_genEta",aLeptons.genEta ,"genEta[nalep]");
937 _outTree->Branch("aLepton_genPhi",aLeptons.genPhi ,"genPhi[nalep]/F");
938 _outTree->Branch("aLepton_charge",aLeptons.charge ,"charge[nalep]/F");
939
940 _outTree->Branch("top" , &top , "mass/F:pt/F:wMass/F");
941 _outTree->Branch("WplusMode" , &WplusMode , "WplusMode/I");
942 _outTree->Branch("WminusMode" , &WminusMode , "WminusMode/I");
943
944 //IVF
945 _outTree->Branch("nSvs",&nSvs ,"nSvs/I");
946 _outTree->Branch("Sv_massBCand", &IVF.massBcand,"massBcand[nSvs]/F");
947 _outTree->Branch("Sv_massSv", &IVF.massSv,"massSv[nSvs]/F");
948 _outTree->Branch("Sv_pt", &IVF.pt,"pt[nSvs]/F");
949 _outTree->Branch("Sv_eta", &IVF.eta,"eta[nSvs]/F");
950 _outTree->Branch("Sv_phi", &IVF.phi,"phi[nSvs]/F");
951 _outTree->Branch("Sv_dist3D", &IVF.dist3D,"dist3D[nSvs]/F");
952 _outTree->Branch("Sv_dist2D", &IVF.dist2D,"dist2D[nSvs]/F");
953 _outTree->Branch("Sv_distSim2D", &IVF.distSig2D,"distSig2D[nSvs]/F");
954 _outTree->Branch("Sv_distSig3D", &IVF.distSig3D,"distSig3D[nSvs]/F");
955 _outTree->Branch("Sv_dist3D_norm", &IVF.dist3D_norm,"dist3D_norm[nSvs]/F");
956 //IVF higgs candidate
957 _outTree->Branch("SVH" , &SVH , "mass/F:pt/F:eta:phi/F:dR/F:dPhi/F:dEta/F");
958
959
960
961 // //SimBHadron
962 _outTree->Branch("nSimBs",&nSimBs ,"nSimBs/I");
963 _outTree->Branch("SimBs_mass", &SimBs.mass,"mass[nSimBs]/F");
964 _outTree->Branch("SimBs_pt", &SimBs.pt,"pt[nSimBs]/F");
965 _outTree->Branch("SimBs_eta", &SimBs.eta,"eta[nSimBs]/F");
966 _outTree->Branch("SimBs_phi", &SimBs.phi,"phi[nSimBs]/F");
967 _outTree->Branch("SimBs_vtx_x", &SimBs.vtx_x,"vtx_x[nSimBs]/F");
968 _outTree->Branch("SimBs_vtx_y", &SimBs.vtx_y,"vtx_y[nSimBs]/F");
969 _outTree->Branch("SimBs_vtx_z", &SimBs.vtx_z,"vtx_z[nSimBs]/F");
970 _outTree->Branch("SimBs_pdgId", &SimBs.pdgId,"pdgId[nSimBs]/F");
971 _outTree->Branch("SimBs_status", &SimBs.status,"status[nSimBs]/F");
972 //SimBHadron Higgs Candidate
973 _outTree->Branch("SimBsH" , &SimBsH , "mass/F:pt/F:eta:phi/F:dR/F:dPhi/F:dEta/F");
974
975 _outTree->Branch("rho" , &rho , "rho/F");
976 _outTree->Branch("rho25" , &rho25 , "rho25/F");
977 _outTree->Branch("nPVs" , &nPVs , "nPVs/I");
978 _outTree->Branch("METnoPU" , &METnoPU , "et/F:sumet:sig/F:phi/F");
979 _outTree->Branch("METnoPUCh" , &METnoPUCh , "et/F:sumet:sig/F:phi/F");
980 _outTree->Branch("MET" , &MET , "et/F:sumet:sig/F:phi/F");
981 _outTree->Branch("METtype1corr" , &METtype1corr , "et/F:sumet:sig/F:phi/F");
982 _outTree->Branch("METtype1p2corr" , &METtype1p2corr , "et/F:sumet:sig/F:phi/F");
983 _outTree->Branch("METnoPUtype1corr" , &METnoPUtype1corr , "et/F:sumet:sig/F:phi/F");
984 _outTree->Branch("METnoPUtype1p2corr" , &METnoPUtype1p2corr , "et/F:sumet:sig/F:phi/F");
985
986 _outTree->Branch("metUnc_et",&metUnc.et ,"et[24]/F");
987 _outTree->Branch("metUnc_phi",&metUnc.phi ,"phi[24]/F");
988 _outTree->Branch("metUnc_sumet",&metUnc.sumet ,"sumet[24]/F");
989
990
991
992 _outTree->Branch("fakeMET" , &fakeMET , "et/F:sumet:sig/F:phi/F");
993 _outTree->Branch("MHT" , &MHT , "mht/F:ht:sig/F:phi/F");
994 _outTree->Branch("minDeltaPhijetMET" , &minDeltaPhijetMET , "minDeltaPhijetMET/F");
995 _outTree->Branch("jetPt_minDeltaPhijetMET" , &jetPt_minDeltaPhijetMET , "jetPt_minDeltaPhijetMET/F");
996
997 std::stringstream s;
998 s << "triggerFlags[" << triggers.size() << "]/b";
999 _outTree->Branch("triggerFlags", triggerFlags, s.str().c_str());
1000
1001 _outTree->Branch("EVENT" , &EVENT , "run/I:lumi/I:event/I:json/I");
1002 _outTree->Branch("hbhe" , &hbhe , "hbhe/b");
1003 _outTree->Branch("totalKinematics" , &totalKinematics , "totalKinematics/b");
1004 _outTree->Branch("ecalFlag" , &ecalFlag , "ecalFlag/b");
1005 _outTree->Branch("cschaloFlag" , &cschaloFlag , "cschaloFlag/b");
1006 _outTree->Branch("hcallaserFlag" , &hcallaserFlag , "hcallaserFlag/b");
1007 _outTree->Branch("trackingfailureFlag" , &trackingfailureFlag , "trackingfailureFlag/b");
1008 _outTree->Branch("btag1TSF" , &btag1TSF , "btag1TSF/F");
1009 _outTree->Branch("btag2TSF" , &btag2TSF , "btag2TSF/F");
1010 _outTree->Branch("btag1T2CSF" , &btag1T2CSF , "btag1T2CSF/F");
1011 _outTree->Branch("btag2CSF" , &btag2CSF , "btag2CSF/F");
1012 _outTree->Branch("btagA0CSF" , &btagA0CSF , "btagA0CSF/F");
1013 _outTree->Branch("btagA0TSF" , &btagA0TSF , "btagA0TSF/F");
1014 _outTree->Branch("btag1TA1C" , &btag1TA1C , "btag1TA1C/F");
1015
1016 int ievt=0;
1017 int totalcount=0;
1018
1019 // TFile* inFile = new TFile(inputFile.c_str(), "read");
1020 for(unsigned int iFile=0; iFile<inputFiles_.size(); ++iFile) {
1021 std::cout << iFile << std::endl;
1022 TFile* inFile = TFile::Open(inputFiles_[iFile].c_str());
1023 if(inFile==0) continue;
1024
1025 // loop the events
1026
1027 fwlite::Run run(inFile);
1028 fwlite::Event ev(inFile);
1029 for(ev.toBegin(); !ev.atEnd() ; ++ev, ++ievt)
1030 {
1031 if (ievt <= skipEvents_) continue;
1032 if (maxEvents_ >= 0){
1033 if (ievt > maxEvents_ + skipEvents_) break;
1034 };
1035
1036 fwlite::Handle< VHbbEventAuxInfo > vhbbAuxHandle;
1037 vhbbAuxHandle.getByLabel(ev,"HbbAnalyzerNew");
1038 const VHbbEventAuxInfo & aux = *vhbbAuxHandle.product();
1039
1040 if(EVENT.run < runMin_ && runMin_ > 0) continue;
1041 if(EVENT.run > runMax_ && runMax_ > 0) continue;
1042
1043 count->Fill(1.);
1044
1045
1046
1047
1048
1049
1050 PUweight=1.;
1051 PUweight2011B=1.;
1052 if(isMC_){
1053
1054 // PU weights // Run2011A
1055 std::map<int, unsigned int>::const_iterator puit = aux.puInfo.pus.find(0);
1056 int npu =puit->second ;
1057 PUweight = lumiWeights.weight( npu );
1058 pu->Fill(puit->second);
1059 // PU weight Run2011B
1060 // PU weight Run2011B
1061 std::map<int, unsigned int>::const_iterator puit0 = aux.puInfo.pus.find(0);
1062 std::map<int, unsigned int>::const_iterator puitm1 = aux.puInfo.pus.find(-1);
1063 std::map<int, unsigned int>::const_iterator puitp1 = aux.puInfo.pus.find(+1);
1064 PU0=puit0->second;
1065 PUp1=puitp1->second;
1066 PUm1=puitm1->second;
1067 input3DPU->Fill(PUm1,PU0,PUp1);
1068 PUweight2011B = lumiWeights2011B.weight3D( puitm1->second, puit0->second,puitp1->second);
1069
1070 }
1071 countWithPU->Fill(1,PUweight);
1072 countWithPU2011B->Fill(1,PUweight2011B);
1073
1074
1075 //LHE Infos
1076 fwlite::Handle<LHEEventProduct> evt;
1077
1078 // std::cout << "Label for lhe = " << evt.getBranchNameFor(ev,"source") << std::endl;
1079 if( !((evt.getBranchNameFor(ev,"source")).empty()) ){
1080 evt.getByLabel(ev,"source");
1081
1082 //std::cout << "LHEEventProduct found!" << std::endl;
1083 bool lCheck=false;
1084 bool lbarCheck=false;
1085 bool vlCheck=false;
1086 bool vlbarCheck=false;
1087 int idl, idlbar;
1088 TLorentzVector l,lbar,vl,vlbar,V_tlv;
1089 const lhef::HEPEUP hepeup_ = evt->hepeup();
1090 double _origECMS = 7000; // energy in GeV
1091 double _newECMS = 8000; // energy in GeV
1092 float Q = hepeup_.SCALUP;
1093 int id1 = hepeup_.IDUP[0];
1094 double x1 = fabs(hepeup_.PUP[0][2]/(_origECMS/2));
1095 double x1prime = fabs(hepeup_.PUP[0][2]/(_newECMS/2));
1096 int id2 = hepeup_.IDUP[1];
1097 double x2 = fabs(hepeup_.PUP[1][2]/(_origECMS/2));
1098 double x2prime = fabs(hepeup_.PUP[1][2]/(_newECMS/2));
1099 // std::cout << "*******LHECOMWeightProducer*******\n" <<
1100 // " Q : " << Q << "\n" <<
1101 // " id1: " << id1 << "\n" <<
1102 // " x1 : " << x1 << "\n" <<
1103 // " x1': " << x1prime << "\n" <<
1104 // " id2: " << id2 << "\n" <<
1105 // " x2 : " << x2 << "\n" <<
1106 // " x2': " << x2prime << std::endl;
1107 //gluon is 0 in the LHAPDF numberin
1108 if (id1 == 21)
1109 id1 = 0;
1110 if (id2 == 21)
1111 id2 = 0;
1112 int _pdfmember = 0;
1113 LHAPDF::usePDFMember(1,_pdfmember);
1114 double oldpdf1 = LHAPDF::xfx(1, x1, Q, id1)/x1;
1115 double oldpdf2 = LHAPDF::xfx(1, x2, Q, id2)/x2;
1116 double newpdf1 = LHAPDF::xfx(1, x1prime, Q, id1)/x1prime;
1117 double newpdf2 = LHAPDF::xfx(1, x2prime, Q, id2)/x2prime;
1118 PDFweight = (newpdf1/oldpdf1)*(newpdf2/oldpdf2);
1119
1120 // std::cout << " xfx1 : " << oldpdf1 << "\n" <<
1121 // " xfx2 : " << oldpdf2 << "\n" <<
1122 // " xfx1': " << newpdf1 << "\n" <<
1123 // " xfx2': " << newpdf2 << "\n" <<
1124 // " weight:" << (newpdf1/oldpdf1)*(newpdf2/oldpdf2) << std::endl;
1125
1126 const std::vector<lhef::HEPEUP::FiveVector> pup_ = hepeup_.PUP; // px, py, pz, E, M
1127 for(unsigned int i=0; i<pup_.size(); ++i){
1128 int id=hepeup_.IDUP[i]; //pdgId
1129
1130 if(id==11){ l.SetPxPyPzE(hepeup_.PUP[i][0],hepeup_.PUP[i][1],hepeup_.PUP[i][2],hepeup_.PUP[i][3]); lCheck=true;}
1131 if(id==-11){ lbar.SetPxPyPzE(hepeup_.PUP[i][0],hepeup_.PUP[i][1],hepeup_.PUP[i][2],hepeup_.PUP[i][3]); lbarCheck=true;}
1132 if(id==12){ vl.SetPxPyPzE(hepeup_.PUP[i][0],hepeup_.PUP[i][1],hepeup_.PUP[i][2],hepeup_.PUP[i][3]); vlCheck=true;}
1133 if(id==-12){ vlbar.SetPxPyPzE(hepeup_.PUP[i][0],hepeup_.PUP[i][1],hepeup_.PUP[i][2],hepeup_.PUP[i][3]); vlbarCheck=true;}
1134
1135 if(id==13){ l.SetPxPyPzE(hepeup_.PUP[i][0],hepeup_.PUP[i][1],hepeup_.PUP[i][2],hepeup_.PUP[i][3]); lCheck=true;}
1136 if(id==-13){ lbar.SetPxPyPzE(hepeup_.PUP[i][0],hepeup_.PUP[i][1],hepeup_.PUP[i][2],hepeup_.PUP[i][3]); lbarCheck=true;}
1137 if(id==14){ vl.SetPxPyPzE(hepeup_.PUP[i][0],hepeup_.PUP[i][1],hepeup_.PUP[i][2],hepeup_.PUP[i][3]); vlCheck=true;}
1138 if(id==-14){ vlbar.SetPxPyPzE(hepeup_.PUP[i][0],hepeup_.PUP[i][1],hepeup_.PUP[i][2],hepeup_.PUP[i][3]); vlbarCheck=true;}
1139
1140 if(id==15){ l.SetPxPyPzE(hepeup_.PUP[i][0],hepeup_.PUP[i][1],hepeup_.PUP[i][2],hepeup_.PUP[i][3]); lCheck=true;}
1141 if(id==-15){ lbar.SetPxPyPzE(hepeup_.PUP[i][0],hepeup_.PUP[i][1],hepeup_.PUP[i][2],hepeup_.PUP[i][3]); lbarCheck=true;}
1142 if(id==16){ vl.SetPxPyPzE(hepeup_.PUP[i][0],hepeup_.PUP[i][1],hepeup_.PUP[i][2],hepeup_.PUP[i][3]); vlCheck=true;}
1143 if(id==-16){ vlbar.SetPxPyPzE(hepeup_.PUP[i][0],hepeup_.PUP[i][1],hepeup_.PUP[i][2],hepeup_.PUP[i][3]); vlbarCheck=true;}
1144
1145 }
1146 if( lCheck && lbarCheck ) V_tlv = l + lbar; // ZtoLL
1147 if( vlCheck && vlbarCheck ) V_tlv = vl + vlbar; // ZtoNuNu
1148 if( lCheck && vlbarCheck ) V_tlv = l + vlbar; // WToLNu
1149 if( lbarCheck && vlCheck ) V_tlv = lbar + vl; // WToLNu
1150 lheV_pt = V_tlv.Pt();
1151 }
1152
1153 //std::cout << "lhe V pt = " << lheV_pt << std::endl;
1154
1155 //Write event info
1156 EVENT.run = ev.id().run();
1157 EVENT.lumi = ev.id().luminosityBlock();
1158 EVENT.event = ev.id().event();
1159 EVENT.json = jsonContainsEvent (jsonVector, ev);
1160
1161 // simBHadrons
1162 const SimBHadronCollection *sbhc;
1163 if(isMC_){
1164 fwlite::Handle<SimBHadronCollection> SBHC;
1165 SBHC.getByLabel(ev, "bhadrons");
1166 sbhc = SBHC.product();
1167 }
1168
1169 const std::vector<VHbbCandidate> * candZ ;
1170 const std::vector<VHbbCandidate> * candW ;
1171 VHbbEvent modifiedEvent;;
1172 const VHbbEvent * iEvent =0;
1173 if(fromCandidate)
1174 {
1175 fwlite::Handle< std::vector<VHbbCandidate> > vhbbCandHandleZ;
1176 vhbbCandHandleZ.getByLabel(ev,"hbbBestCSVPt20Candidates");
1177 candZ = vhbbCandHandleZ.product();
1178
1179 fwlite::Handle< std::vector<VHbbCandidate> > vhbbCandHandle;
1180 vhbbCandHandle.getByLabel(ev,"hbbHighestPtHiggsPt30Candidates");
1181 candW = vhbbCandHandle.product();
1182 }
1183 else
1184 {
1185 candZlocal->clear();
1186 candWlocal->clear();
1187 fwlite::Handle< VHbbEvent > vhbbHandle;
1188 vhbbHandle.getByLabel(ev,"HbbAnalyzerNew");
1189 modifiedEvent = *vhbbHandle.product();
1190 if(isMC_)
1191 {
1192 iEvent= &modifiedEvent;
1193
1194 for(size_t j=0; j< modifiedEvent.simpleJets2.size() ; j++)
1195 {
1196 TLorentzVector & p4 = modifiedEvent.simpleJets2[j].p4;
1197 TLorentzVector & mcp4 = modifiedEvent.simpleJets2[j].bestMCp4;
1198 if ((fabs(p4.Pt() - mcp4.Pt())/ p4.Pt())<0.5) { //Limit the effect to the core
1199 float cor = (p4.Pt()+resolutionBias(fabs(p4.Eta()))*(p4.Pt()-mcp4.Pt()))/p4.Pt();
1200 p4.SetPtEtaPhiE(p4.Pt()*cor,p4.Eta(), p4.Phi(), p4.E()*cor);
1201 }
1202 }
1203 } else
1204 {
1205 iEvent = vhbbHandle.product();
1206 }
1207
1208 algoZ->run(iEvent,*candZlocal);
1209 algoW->run(iEvent,*candWlocal);
1210
1211 if(candZlocal->size() == 0 or candZlocal->at(0).H.jets.size() < 2) //recover low pt
1212 {
1213 candZlocal->clear();
1214 candWlocal->clear();
1215 algoRecoverLowPt->run(iEvent,*candZlocal);
1216 algoRecoverLowPt->run(iEvent,*candWlocal);
1217 }
1218
1219 candZ= candZlocal;
1220 candW= candWlocal;
1221 /* for(size_t m=0;m<iEvent->muInfo.size();m++)
1222 {
1223
1224 if( fabs(iEvent->muInfo[m].p4.Pt()-28.118684) < 0.0001 ||
1225 fabs(iEvent->muInfo[m].p4.Pt()-34.853199) < 0.0001 )
1226 {
1227 std::cout << "FOUND " << iEvent->muInfo[m].p4.Pt() << " " << EVENT.event << " " << candW->size() << " " << candZ->size() << std::endl;
1228 }
1229 }
1230 */
1231
1232 }
1233
1234 const std::vector<VHbbCandidate> * cand = candZ;
1235
1236
1237
1238 /* fwlite::Handle< VHbbEvent > vhbbHandle;
1239 vhbbHandle.getByLabel(ev,"HbbAnalyzerNew");
1240 const VHbbEvent iEvent = *vhbbHandle.product();
1241 */
1242
1243 // std::clog << "Filling tree "<< std::endl;
1244 bool isW=false;
1245
1246 // to check how much we gain with jets subtraction
1247
1248
1249
1250 genHpt=aux.mcH.size() > 0 ? aux.mcH[0].p4.Pt():-99;
1251
1252 if(cand->size() == 0 or cand->at(0).H.jets.size() < 2) continue;
1253 //if(cand->size() == 0 ) continue;
1254 //std::cout << "cand->size() " << cand->size() << std::endl;
1255 //std::cout << "cand->at(0).H.jets.size() " << cand->at(0).H.jets.size() << std::endl;
1256 if(cand->size() > 1 )
1257 {
1258 std::cout << "MULTIPLE CANDIDATES: " << cand->size() << std::endl;
1259 }
1260 if(cand->at(0).candidateType == VHbbCandidate::Wmun || cand->at(0).candidateType == VHbbCandidate::Wen ) { cand=candW; isW=true; }
1261 if(cand->size() == 0)
1262 {
1263 // std::cout << "W event loss due to tigther cuts" << std::endl;
1264 continue;
1265 }
1266
1267
1268 // secondary vtxs
1269 fwlite::Handle<std::vector<reco::Vertex> > SVC;
1270 SVC.getByLabel(ev,"bcandidates");
1271 const std::vector<reco::Vertex> svc = *(SVC.product());
1272
1273 const VHbbCandidate & vhCand = cand->at(0);
1274 patFilters.setEvent(&ev,"VH");
1275 hbhe = patFilters.accept("hbhe");
1276 ecalFlag = patFilters.accept("ecalFilter");
1277 totalKinematics = patFilters.accept("totalKinematics");
1278 cschaloFlag = patFilters.accept("cschaloFilter");
1279 hcallaserFlag = patFilters.accept("hcallaserFilter");
1280 trackingfailureFlag = patFilters.accept("trackingfailureFilter");
1281
1282 trigger.setEvent(&ev);
1283 for(size_t j=0;j < triggers.size();j++)
1284 triggerFlags[j]=trigger.accept(triggers[j]);
1285
1286 eventFlav=0;
1287
1288
1289
1290 if(aux.mcBbar.size() > 0 || aux.mcB.size() > 0) eventFlav=5;
1291 else if(aux.mcC.size() > 0) eventFlav=4;
1292
1293
1294 Vtype = vhCand.candidateType;
1295
1296 if(vhCand.H.HiggsFlag) H.HiggsFlag=1; else H.HiggsFlag=0;
1297
1298 if(vhCand.H.HiggsFlag){
1299 H.mass = vhCand.H.p4.M();
1300 H.pt = vhCand.H.p4.Pt();
1301
1302
1303 H.eta = vhCand.H.p4.Eta();
1304 H.phi = vhCand.H.p4.Phi();
1305 }
1306
1307 if(vhCand.FatH.FatHiggsFlag) FatH.FatHiggsFlag =1; else FatH.FatHiggsFlag=0;
1308 fathFilterJets.reset();
1309 aJetsFat.reset();
1310 if(vhCand.FatH.FatHiggsFlag){
1311 FatH.mass= vhCand.FatH.p4.M();
1312 FatH.pt = vhCand.FatH.p4.Pt();
1313 FatH.eta = vhCand.FatH.p4.Eta();
1314 FatH.phi = vhCand.FatH.p4.Phi();
1315
1316 // if(vhCand.FatH.FatHiggsFlag) vhCand.FatH.subjetsSize;
1317 nfathFilterJets=vhCand.FatH.subjetsSize;
1318 for( int j=0; j < nfathFilterJets; j++ ){
1319 fathFilterJets.set(vhCand.FatH.jets[j],j);
1320 }
1321
1322 if(nfathFilterJets==2){
1323 FatH.filteredmass=(vhCand.FatH.jets[0].p4+vhCand.FatH.jets[1].p4).M();
1324 FatH.filteredpt=(vhCand.FatH.jets[0].p4+vhCand.FatH.jets[1].p4).Pt();
1325 FatH.filteredeta=(vhCand.FatH.jets[0].p4+vhCand.FatH.jets[1].p4).Eta();
1326 FatH.filteredphi=(vhCand.FatH.jets[0].p4+vhCand.FatH.jets[1].p4).Phi();
1327 }
1328 else if(nfathFilterJets==3){
1329 FatH.filteredmass=(vhCand.FatH.jets[0].p4+vhCand.FatH.jets[1].p4+vhCand.FatH.jets[2].p4).M();
1330 FatH.filteredpt=(vhCand.FatH.jets[0].p4+vhCand.FatH.jets[1].p4+vhCand.FatH.jets[2].p4).Pt();
1331 FatH.filteredeta=(vhCand.FatH.jets[0].p4+vhCand.FatH.jets[1].p4+vhCand.FatH.jets[2].p4).Eta();
1332 FatH.filteredphi=(vhCand.FatH.jets[0].p4+vhCand.FatH.jets[1].p4+vhCand.FatH.jets[2].p4).Phi();
1333 }
1334
1335 naJetsFat=vhCand.additionalJetsFat.size();
1336 for( int j=0; j < naJetsFat && j < MAXJ; j++ )
1337 {
1338 aJetsFat.set(vhCand.additionalJetsFat[j],j);
1339 }
1340
1341
1342 } // FatHiggsFlag
1343
1344 hJets.reset();
1345 aJets.reset();
1346 if(vhCand.H.HiggsFlag){
1347
1348 nhJets=2;
1349 hJets.set(vhCand.H.jets[0],0);
1350 hJets.set(vhCand.H.jets[1],1);
1351
1352
1353 aJets.reset();
1354
1355 naJets=vhCand.additionalJets.size();
1356 numBJets=0;
1357 if(vhCand.H.jets[0].csv> btagThr) numBJets++;
1358 if(vhCand.H.jets[1].csv> btagThr) numBJets++;
1359 for( int j=0; j < naJets && j < MAXJ; j++ )
1360 {
1361 aJets.set(vhCand.additionalJets[j],j);
1362 if(vhCand.additionalJets[j].csv> btagThr) numBJets++;
1363 }
1364 numJets = vhCand.additionalJets.size()+2;
1365 H.dR = deltaR(vhCand.H.jets[0].p4.Eta(),vhCand.H.jets[0].p4.Phi(),vhCand.H.jets[1].p4.Eta(),vhCand.H.jets[1].p4.Phi());
1366 H.dPhi = deltaPhi(vhCand.H.jets[0].p4.Phi(),vhCand.H.jets[1].p4.Phi());
1367 H.dEta= TMath::Abs( vhCand.H.jets[0].p4.Eta() - vhCand.H.jets[1].p4.Eta() );
1368 HVdPhi = fabs( deltaPhi(vhCand.H.p4.Phi(),vhCand.V.p4.Phi()) ) ;
1369 HVMass = (vhCand.H.p4 + vhCand.V.p4).M() ;
1370 HMETdPhi = fabs( deltaPhi(vhCand.H.p4.Phi(),vhCand.V.mets.at(0).p4.Phi()) ) ;
1371 //eltaPullAngle = vhCand.H.deltaTheta;
1372
1373 deltaPullAngle = VHbbCandidateTools::getDeltaTheta(vhCand.H.jets[0],vhCand.H.jets[1]);
1374 deltaPullAngle2 = VHbbCandidateTools::getDeltaTheta(vhCand.H.jets[1],vhCand.H.jets[0]);
1375 hJets.cosTheta[0]= vhCand.H.helicities[0];
1376 hJets.cosTheta[1]= vhCand.H.helicities[1];
1377 } // Higgs Flag
1378
1379 V.mass = vhCand.V.p4.M();
1380 if(isW) V.mass = vhCand.Mt();
1381 V.pt = vhCand.V.p4.Pt();
1382 V.eta = vhCand.V.p4.Eta();
1383 V.phi = vhCand.V.p4.Phi();
1384 VMt = vhCand.Mt() ;
1385
1386
1387 // METInfo calomet; METInfo tcmet; METInfo pfmet; METInfo mht; METInfo metNoPU
1388 MET.et = vhCand.V.mets.at(0).p4.Pt();
1389 MET.phi = vhCand.V.mets.at(0).p4.Phi();
1390 MET.sumet = vhCand.V.mets.at(0).sumEt;
1391 MET.sig = vhCand.V.mets.at(0).metSig;
1392
1393
1394 fakeMET.sumet = 0;
1395 fakeMET.sig = 0;
1396 fakeMET.et = 0;
1397 fakeMET.phi = 0;
1398 if( Vtype == VHbbCandidate::Zmumu) {
1399 TVector3 mu1 = vhCand.V.muons[0].p4.Vect();
1400 TVector3 mu2 = vhCand.V.muons[1].p4.Vect();
1401 // Not needed with PFMET
1402 // mu1.SetMag( mu1.Mag() - vhCand.V.muons[0].emEnergy - vhCand.V.muons[0].hadEnergy);
1403 // mu2.SetMag( mu2.Mag() - vhCand.V.muons[1].emEnergy - vhCand.V.muons[1].hadEnergy);
1404 TVector3 sum = vhCand.V.mets.at(0).p4.Vect() + mu1 + mu2;
1405 fakeMET.et = sum.Pt();
1406 fakeMET.phi = sum.Phi();
1407 fakeMET.sumet = vhCand.V.mets.at(0).sumEt - mu1.Pt() - mu2.Pt();
1408 }
1409
1410
1411
1412 METnoPU.et = iEvent->metNoPU.p4.Pt();
1413 METnoPU.phi = iEvent->metNoPU.p4.Phi();
1414 METnoPU.sumet = iEvent->metNoPU.sumEt;
1415 METnoPU.sig = iEvent->metNoPU.metSig;
1416 METnoPUCh.et = iEvent->metCh.p4.Pt();
1417 METnoPUCh.phi = iEvent->metCh.p4.Phi();
1418 METnoPUCh.sumet = iEvent->metCh.sumEt;
1419 METnoPUCh.sig = iEvent->metCh.metSig;
1420
1421 METnoPUCh.et = iEvent->metCh.p4.Pt();
1422 METnoPUCh.phi = iEvent->metCh.p4.Phi();
1423 METnoPUCh.sumet = iEvent->metCh.sumEt;
1424 METnoPUCh.sig = iEvent->metCh.metSig;
1425
1426
1427 METtype1corr.et = iEvent->pfmetType1corr.p4.Pt();
1428 METtype1corr.phi = iEvent->pfmetType1corr.p4.Phi();
1429 METtype1corr.sumet = iEvent->pfmetType1corr.sumEt;
1430 METtype1corr.sig = iEvent->pfmetType1corr.metSig;
1431
1432
1433 METtype1p2corr.et = iEvent->pfmetType1p2corr.p4.Pt();
1434 METtype1p2corr.phi = iEvent->pfmetType1p2corr.p4.Phi();
1435 METtype1p2corr.sumet = iEvent->pfmetType1p2corr.sumEt;
1436 METtype1p2corr.sig = iEvent->pfmetType1p2corr.metSig;
1437
1438
1439 METnoPUtype1corr.et = iEvent->pfmetNoPUType1corr.p4.Pt();
1440 METnoPUtype1corr.phi = iEvent->pfmetNoPUType1corr.p4.Phi();
1441 METnoPUtype1corr.sumet = iEvent->pfmetNoPUType1corr.sumEt;
1442 METnoPUtype1corr.sig = iEvent->pfmetNoPUType1corr.metSig;
1443
1444
1445 METnoPUtype1p2corr.et = iEvent->pfmetNoPUType1p2corr.p4.Pt();
1446 METnoPUtype1p2corr.phi = iEvent->pfmetNoPUType1p2corr.p4.Phi();
1447 METnoPUtype1p2corr.sumet = iEvent->pfmetNoPUType1p2corr.sumEt;
1448 METnoPUtype1p2corr.sig = iEvent->pfmetNoPUType1p2corr.metSig;
1449
1450 // std::cout << " iEvent->metUncInfo.size() " << iEvent->metUncInfo.size() << std::endl;
1451 for(size_t m=0;m<iEvent->metUncInfo.size();m++)
1452 {
1453 metUnc.set(iEvent->metUncInfo[m], m );
1454 // std::cout << "metUncInfo[" << m <<" ].et = " << metUnc.et[m] << std::endl;
1455 }
1456
1457
1458 rho = aux.puInfo.rho;
1459 rho25 = aux.puInfo.rho25;
1460 nPVs=aux.pvInfo.nVertices;
1461
1462 if(!fromCandidate) {
1463 MHT.mht = iEvent->mht.p4.Pt();
1464 MHT.phi = iEvent->mht.p4.Phi();
1465 MHT.ht = iEvent->mht.sumEt;
1466 MHT.sig = iEvent->mht.metSig;
1467 }
1468
1469
1470 /////////
1471 // track sharing flags:
1472 ////////
1473 TkSharing.HiggsCSVtkSharing = TkSharing.HiggsIPtkSharing = TkSharing.HiggsSVtkSharing = TkSharing.FatHiggsCSVtkSharing = TkSharing.FatHiggsIPtkSharing = TkSharing.FatHiggsSVtkSharing = false;
1474
1475 // csv tracks
1476 if(vhCand.H.HiggsFlag){
1477
1478 if (vhCand.H.jets[0].csvNTracks > 0 && vhCand.H.jets[1].csvNTracks > 0){
1479 for (int t=0;t!=vhCand.H.jets[0].csvNTracks;t++){
1480 for (int ti=0;ti!=vhCand.H.jets[1].csvNTracks;ti++){
1481 if ((int)vhCand.H.jets[0].csvTrackIds[t] == (int)vhCand.H.jets[1].csvTrackIds[ti]){
1482 TkSharing.HiggsCSVtkSharing = true;
1483 }// same trackID
1484 }// loop tracks in second hjet
1485 }// loop tracks in first hjet
1486 }// if tracks in jet
1487
1488 // ip tracks
1489 if (vhCand.H.jets[0].btagNTracks > 0 && vhCand.H.jets[1].btagNTracks > 0){
1490 for (int t=0;t!=vhCand.H.jets[0].btagNTracks;t++){
1491 for (int ti=0;ti!=vhCand.H.jets[1].btagNTracks;ti++){
1492 if ((int)vhCand.H.jets[0].btagTrackIds[t] == (int)vhCand.H.jets[1].btagTrackIds[ti]){
1493 TkSharing.HiggsIPtkSharing = true;
1494 }// same trackID
1495 }// loop tracks in second hjet
1496 }// loop tracks in first hjet
1497 }// if tracks in jet
1498
1499 // sv tracks
1500 if (vhCand.H.jets[0].vtxNTracks > 0 && vhCand.H.jets[1].vtxNTracks > 0){
1501 for (int t=0;t!=vhCand.H.jets[0].vtxNTracks;t++){
1502 for (int ti=0;ti!=vhCand.H.jets[1].vtxNTracks;ti++){
1503 if ((int)vhCand.H.jets[0].vtxTrackIds[t] == (int)vhCand.H.jets[1].vtxTrackIds[ti]){
1504 TkSharing.HiggsSVtkSharing = true;
1505 }// same trackID
1506 }// loop tracks in second hjet
1507 }// loop tracks in first hjet
1508 }// if tracks in jet
1509
1510 } // Di-jet Higgs Flag
1511
1512 // tracksharing for Filtered jets:
1513 if(vhCand.FatH.FatHiggsFlag && nfathFilterJets > 1){
1514
1515 // csv tracks
1516 if (vhCand.FatH.jets[0].csvNTracks > 0 && vhCand.FatH.jets[1].csvNTracks > 0){
1517 for (int t=0;t!=vhCand.FatH.jets[0].csvNTracks;t++){
1518 for (int ti=0;ti!=vhCand.FatH.jets[1].csvNTracks;ti++){
1519 if ((int)vhCand.FatH.jets[0].csvTrackIds[t] == (int)vhCand.FatH.jets[1].csvTrackIds[ti]){
1520 TkSharing.FatHiggsCSVtkSharing = true;
1521 }// same trackID
1522 }// loop tracks in second hjet
1523 }// loop tracks in first hjet
1524 }// if tracks in jet
1525
1526 // ip tracks
1527 if (vhCand.FatH.jets[0].btagNTracks > 0 && vhCand.FatH.jets[1].btagNTracks > 0){
1528 for (int t=0;t!=vhCand.FatH.jets[0].btagNTracks;t++){
1529 for (int ti=0;ti!=vhCand.FatH.jets[1].btagNTracks;ti++){
1530 if ((int)vhCand.FatH.jets[0].btagTrackIds[t] == (int)vhCand.FatH.jets[1].btagTrackIds[ti]){
1531 TkSharing.FatHiggsIPtkSharing = true;
1532 }// same trackID
1533 }// loop tracks in second hjet
1534 }// loop tracks in first hjet
1535 }// if tracks in jet
1536
1537 // sv tracks
1538 if (vhCand.FatH.jets[0].vtxNTracks > 0 && vhCand.FatH.jets[1].vtxNTracks > 0){
1539 for (int t=0;t!=vhCand.FatH.jets[0].vtxNTracks;t++){
1540 for (int ti=0;ti!=vhCand.FatH.jets[1].vtxNTracks;ti++){
1541 if ((int)vhCand.FatH.jets[0].vtxTrackIds[t] == (int)vhCand.FatH.jets[1].vtxTrackIds[ti]){
1542 TkSharing.FatHiggsSVtkSharing = true;
1543 }// same trackID
1544 }// loop tracks in second hjet
1545 }// loop tracks in first hjet
1546 }// if tracks in jet
1547
1548 }// fatH
1549
1550 ////////
1551 ////////
1552
1553 //Secondary Vertices
1554 IVF.reset();
1555 nSvs = svc.size();
1556 const TVector3 recoPv = aux.pvInfo.firstPVInPT2;
1557 const math::XYZPoint myPv(recoPv);
1558
1559 //FAKE ERROR MATRIX
1560 // //look here for Matrix filling info http://project-mathlibs.web.cern.ch/project-mathlibs/sw/html/SMatrixDoc.html
1561 // std::vector<double> fillMatrix(6);
1562 // for (int i = 0; i<6; ++i) fillMatrix[i] = 0.;
1563 // fillMatrix[0] = TMath::Power(0.002,2);
1564 // fillMatrix[2] = TMath::Power(0.002,2);
1565 // fillMatrix[5] = TMath::Power(0.002,2);
1566 // const ROOT::Math::SMatrix<double, 3, 3, ROOT::Math::MatRepSym<double, 3> > myFakeMatrixError(fillMatrix.begin(),fillMatrix.end());
1567 // const reco::Vertex recoVtxPv(myPv, myFakeMatrixError);
1568
1569 // REAL ERROR MATRIX
1570 const reco::Vertex recoVtxPv(myPv, aux.pvInfo.efirstPVInPT2);
1571 for( int j=0; j < nSvs && j < MAXB; ++j ) {
1572 const GlobalVector flightDir = flightDirection(recoPv,svc[j]);
1573 reco::SecondaryVertex recoSv(recoVtxPv, svc[j], flightDir ,true);
1574 IVF.set( recoSv, recoPv ,j);
1575 }
1576 if(nSvs > 1){
1577 TLorentzVector BCands_H1, BCands_H2, BCands_H;
1578 BCands_H1.SetPtEtaPhiM(IVF.pt[0], IVF.eta[0], IVF.phi[0], IVF.massBcand[0]);
1579 BCands_H2.SetPtEtaPhiM(IVF.pt[1], IVF.eta[1], IVF.phi[1], IVF.massBcand[1]);
1580 BCands_H = BCands_H1 + BCands_H2;
1581 SVH.dR = deltaR(IVF.eta[0], IVF.phi[0], IVF.eta[1], IVF.phi[1] );
1582 SVH.dPhi = deltaPhi(IVF.phi[0], IVF.phi[1] );
1583 SVH.dEta = TMath::Abs(IVF.eta[0] - IVF.eta[1] );
1584 SVH.mass = BCands_H.M();
1585 SVH.pt = BCands_H.Pt();
1586 SVH.eta = BCands_H.Eta();
1587 SVH.phi = BCands_H.Phi();
1588 }
1589
1590 //SimBhadron
1591 SimBs.reset();
1592 if(isMC_){
1593 nSimBs = sbhc->size();
1594 for( int j=0; j < nSimBs && j < MAXB; ++j )
1595 SimBs.set( sbhc->at(j), j);
1596 if(nSimBs > 1){
1597 TLorentzVector SimBs_H1, SimBs_H2, SimBs_H;
1598 SimBs_H1.SetPtEtaPhiM(SimBs.pt[0], SimBs.eta[0], SimBs.phi[0], SimBs.mass[0]);
1599 SimBs_H2.SetPtEtaPhiM(SimBs.pt[1], SimBs.eta[1], SimBs.phi[1], SimBs.mass[1]);
1600 SimBs_H = SimBs_H1 + SimBs_H2;
1601 SimBsH.dR = deltaR(SimBs.eta[0], SimBs.phi[0], SimBs.eta[1], SimBs.phi[1] );
1602 SimBsH.dPhi = deltaPhi(SimBs.phi[0], SimBs.phi[1] );
1603 SimBsH.dEta = TMath::Abs(SimBs.eta[0] - SimBs.eta[1] );
1604 SimBsH.mass = SimBs_H.M();
1605 SimBsH.pt = SimBs_H.Pt();
1606 SimBsH.eta = SimBs_H.Eta();
1607 SimBsH.phi = SimBs_H.Phi();
1608 }
1609 }
1610
1611 //Loop on jets
1612 double maxBtag=-99999;
1613 minDeltaPhijetMET = 999;
1614 TLorentzVector bJet;
1615 std::vector<std::vector<BTagWeight::JetInfo> > btagJetInfos;
1616 std::vector<float> jet10eta;
1617 std::vector<float> jet10pt;
1618 std::vector<float> jet30eta;
1619 std::vector<float> jet30pt;
1620 if(fromCandidate)
1621 {
1622 //Loop on Higgs Jets
1623 for(unsigned int j=0; j < vhCand.H.jets.size(); j++ ){
1624 if (vhCand.H.jets[j].csv > maxBtag) { bJet=vhCand.H.jets[j].p4 ; maxBtag =vhCand.H.jets[j].csv; }
1625 if (fabs(deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.H.jets[j].p4.Phi())) < minDeltaPhijetMET)
1626 {
1627 minDeltaPhijetMET=fabs(deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.H.jets[j].p4.Phi()));
1628 jetPt_minDeltaPhijetMET=vhCand.H.jets[j].p4.Pt();
1629 }
1630 btagJetInfos.push_back(btagEff.jetInfo(vhCand.H.jets[j]));
1631 }
1632 //Loop on Additional Jets
1633 for(unsigned int j=0; j < vhCand.additionalJets.size(); j++ ){
1634 if (vhCand.additionalJets[j].csv > maxBtag) { bJet=vhCand.additionalJets[j].p4 ; maxBtag =vhCand.additionalJets[j].csv; }
1635 /* if (fabs(deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.additionalJets[j].p4.Phi())) < minDeltaPhijetMET)
1636 {
1637 minDeltaPhijetMET=fabs(deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.additionalJets[j].p4.Phi()));
1638 jetPt_minDeltaPhijetMET=vhCand.additionalJets[j].p4.Pt();
1639 }*/
1640 if( ( isW && ! useHighestPtHiggsW ) || ( ! isW && ! useHighestPtHiggsZ ) ) // btag SF computed using only H-jets if best-H made with dijetPt rather than best CSV
1641 {
1642 if(vhCand.additionalJets[j].p4.Pt() > 20)
1643 btagJetInfos.push_back(btagEff.jetInfo(vhCand.additionalJets[j]));
1644 }
1645 }
1646 }
1647 else
1648 {
1649 //Loop on all jets
1650 for(unsigned int j=0; j < iEvent->simpleJets2.size(); j++ ){
1651 if (iEvent->simpleJets2[j].csv > maxBtag) { bJet=iEvent->simpleJets2[j].p4 ; maxBtag =iEvent->simpleJets2[j].csv; }
1652 if ( iEvent->simpleJets2[j].p4.Pt() > 20 && fabs(iEvent->simpleJets2[j].p4.Eta()) < 2.5&& fabs(deltaPhi( vhCand.V.mets.at(0).p4.Phi(), iEvent->simpleJets2[j].p4.Phi())) < minDeltaPhijetMET)
1653 {
1654 minDeltaPhijetMET=fabs(deltaPhi( vhCand.V.mets.at(0).p4.Phi(), iEvent->simpleJets2[j].p4.Phi()));
1655 jetPt_minDeltaPhijetMET=iEvent->simpleJets2[j].p4.Pt();
1656 }
1657 if(iEvent->simpleJets2[j].p4.Pt() > 10)
1658 {
1659 jet10eta.push_back(iEvent->simpleJets2[j].p4.Eta());
1660 jet10pt.push_back(iEvent->simpleJets2[j].p4.Pt());
1661 }
1662 if(iEvent->simpleJets2[j].p4.Pt() > 30)
1663 {
1664 jet30eta.push_back(iEvent->simpleJets2[j].p4.Eta());
1665 jet30pt.push_back(iEvent->simpleJets2[j].p4.Pt());
1666 }
1667
1668 //For events made with highest CSV, all jets in the event should be taken into account for "tagging" SF (anti tagging is a mess)
1669 // because for example a light jet not used for the Higgs can have in reality a higher CSV due to SF > 1 and become a higgs jet
1670 if( ( isW && ! useHighestPtHiggsW ) || ( ! isW && ! useHighestPtHiggsZ ) )
1671 {
1672 if(iEvent->simpleJets2[j].p4.Pt() > 20 && fabs(iEvent->simpleJets2[j].p4.Eta()) < 2.5)
1673 btagJetInfos.push_back(btagEff.jetInfo(iEvent->simpleJets2[j]));
1674 }
1675 }
1676
1677 //Loop on Higgs jets
1678
1679 if(vhCand.H.HiggsFlag){
1680 for(unsigned int j=0; j < vhCand.H.jets.size(); j++ ) {
1681
1682 //if we use the highest pt pair, only the two higgs jet should be used to compute the SF because the other jets are excluded
1683 // by a criteria (pt of the dijet) that is not btag SF dependent
1684 if(!( ( isW && ! useHighestPtHiggsW ) || ( ! isW && ! useHighestPtHiggsZ ) )) {
1685 btagJetInfos.push_back(btagEff.jetInfo(vhCand.H.jets[j]));
1686 }
1687 }
1688 }// HiggsFlag
1689
1690 }
1691 vLeptons.reset();
1692 weightTrig = 1.; // better to default to 1
1693 weightTrigMay = -1.;
1694 weightTrigV4 = -1.;
1695 weightTrigOrMu30 = 1.;
1696 TLorentzVector leptonForTop;
1697 size_t firstAddMu=0;
1698 size_t firstAddEle=0;
1699 if(Vtype == VHbbCandidate::Zmumu ){
1700 vLeptons.set(vhCand.V.muons[0],0,13);
1701 vLeptons.set(vhCand.V.muons[1],1,13);
1702 float cweightID = triggerWeight.scaleMuID(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleMuID(vLeptons.pt[1],vLeptons.eta[1]) ;
1703 float weightTrig1 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[0],vLeptons.eta[0]);
1704 float weightTrig2 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[1],vLeptons.eta[1]);
1705 float cweightTrig = weightTrig1 + weightTrig2 - weightTrig1*weightTrig2;
1706 weightTrig = cweightID * cweightTrig;
1707 nvlep=2;
1708 firstAddMu=2;
1709 }
1710 if( Vtype == VHbbCandidate::Zee ){
1711 vLeptons.set(vhCand.V.electrons[0],0,11);
1712 vLeptons.set(vhCand.V.electrons[1],1,11);
1713 nvlep=2;
1714 firstAddEle=2;
1715 std::vector<float> pt,eta;
1716 pt.push_back(vLeptons.pt[0]); eta.push_back(vLeptons.eta[0]);
1717 pt.push_back(vLeptons.pt[1]); eta.push_back(vLeptons.eta[1]);
1718 weightEleRecoAndId=triggerWeight.scaleID95Ele(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleRecoEle(vLeptons.pt[0],vLeptons.eta[0]) *
1719 triggerWeight.scaleID95Ele(vLeptons.pt[1],vLeptons.eta[1]) * triggerWeight.scaleRecoEle(vLeptons.pt[1],vLeptons.eta[1]);
1720 weightEleTrigElePart = triggerWeight.scaleDoubleEle17Ele8(pt,eta);
1721 weightEleTrigEleAugPart = triggerWeight.scaleDoubleEle17Ele8Aug(pt,eta);
1722 weightTrig = (weightEleTrigElePart*1.14+weightEleTrigEleAugPart*0.98 )/2.12 * weightEleRecoAndId;
1723
1724
1725
1726 }
1727 if(Vtype == VHbbCandidate::Wmun ){
1728 leptonForTop=vhCand.V.muons[0].p4;
1729 vLeptons.set(vhCand.V.muons[0],0,13);
1730 float cweightID = triggerWeight.scaleMuID(vLeptons.pt[0],vLeptons.eta[0]);
1731 float weightTrig1 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[0],vLeptons.eta[0]);
1732 float cweightTrig = weightTrig1;
1733 weightTrig = cweightID * cweightTrig;
1734 float weightTrig1OrMu30 = triggerWeight.scaleMuOr30IsoHLT(vLeptons.pt[0],vLeptons.eta[0]);
1735 weightTrigOrMu30 = cweightID*weightTrig1OrMu30;
1736 nvlep=1;
1737 firstAddMu=1;
1738 }
1739 if( Vtype == VHbbCandidate::Wen ){
1740 leptonForTop=vhCand.V.electrons[0].p4;
1741 vLeptons.set(vhCand.V.electrons[0],0,11);
1742 nvlep=1;
1743 firstAddEle=1;
1744 weightTrigMay = triggerWeight.scaleSingleEleMay(vLeptons.pt[0],vLeptons.eta[0]);
1745 weightTrigV4 = triggerWeight.scaleSingleEleV4(vLeptons.pt[0],vLeptons.eta[0]);
1746 weightEleRecoAndId=triggerWeight.scaleID80Ele(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleRecoEle(vLeptons.pt[0],vLeptons.eta[0]);
1747 weightEleTrigJetMETPart=triggerWeight.scaleJet30Jet25(jet30pt,jet30eta)*triggerWeight.scalePFMHTEle(MET.et);
1748 weightEleTrigElePart= weightTrigV4; //this is for debugging only, checking only the V4 part
1749
1750 weightTrigMay*=weightEleRecoAndId;
1751 weightTrigV4*=weightEleRecoAndId;
1752 weightTrigV4*=weightEleTrigJetMETPart;
1753 // weightTrig = weightTrigMay * 0.187 + weightTrigV4 * (1.-0.187); //FIXME: use proper lumi if we reload 2.fb
1754 weightTrig = (weightTrigMay * 0.215 + weightTrigV4 * 1.915)/ 2.13; //FIXME: use proper lumi if we reload 2.fb
1755 }
1756
1757 if(isMC_)
1758 {
1759 weightTrigMET80 = triggerWeight.scaleMET80(MET.et);
1760 weightTrigMET100 = triggerWeight.scaleMET80(MET.et);
1761 weightTrig2CJet20 = triggerWeight.scale2CentralJet( jet10pt, jet10eta);
1762 weightTrigMET150 = triggerWeight.scaleMET150(MET.et);
1763 weightTrigMET802CJet= weightTrigMET80 * weightTrig2CJet20;
1764 weightTrigMET1002CJet= weightTrigMET100 * weightTrig2CJet20;
1765 }
1766 if( Vtype == VHbbCandidate::Znn ){
1767 nvlep=0;
1768 float weightTrig1 = triggerWeight.scaleMetHLT(vhCand.V.mets.at(0).p4.Pt());
1769 weightTrigMETLP = weightTrig1;
1770 weightTrig = weightTrigMET150 + weightTrigMET802CJet - weightTrigMET802CJet*weightTrigMET150;
1771 }
1772
1773 if(weightTrigMay < 0) weightTrigMay=weightTrig;
1774 if(weightTrigV4 < 0) weightTrigV4=weightTrig;
1775 if(!isMC_)
1776 {
1777 weightTrig = 1.;
1778 weightTrigMay = 1.;
1779 weightTrigV4 = 1.;
1780 weightEleRecoAndId= 1.;
1781 weightEleTrigJetMETPart= 1.;
1782 weightEleTrigElePart= 1.;
1783 weightEleTrigEleAugPart=1.;
1784 weightTrigMET80= 1.;
1785 weightTrigMET100= 1.;
1786 weightTrig2CJet20= 1.;
1787 weightTrigMET150= 1.;
1788 weightTrigMET802CJet= 1.;
1789 weightTrigMET1002CJet= 1.;
1790 weightTrigMETLP = 1.;
1791
1792 }
1793 aLeptons.reset();
1794 nalep=0;
1795 if(fromCandidate)
1796 {
1797 for(size_t j=firstAddMu;j< vhCand.V.muons.size();j++) aLeptons.set(vhCand.V.muons[j],nalep++,13);
1798 for(size_t j=firstAddEle;j< vhCand.V.electrons.size();j++) aLeptons.set(vhCand.V.electrons[j],nalep++,11);
1799 }
1800 else
1801 {
1802 for(size_t j=0;j< iEvent->muInfo.size();j++)
1803 {
1804 if((j!= vhCand.V.firstLepton && j!= vhCand.V.secondLepton) || ((Vtype != VHbbCandidate::Wmun ) && (Vtype != VHbbCandidate::Zmumu )) )
1805 aLeptons.set(iEvent->muInfo[j],nalep++,13);
1806 }
1807 for(size_t j=0;j< iEvent->eleInfo.size();j++)
1808 {
1809 if((j!= vhCand.V.firstLepton && j!= vhCand.V.secondLepton) || ((Vtype != VHbbCandidate::Wen ) && (Vtype != VHbbCandidate::Zee )))
1810 aLeptons.set(iEvent->eleInfo[j],nalep++,11);
1811 }
1812
1813 }
1814
1815
1816 if(isMC_)
1817 {
1818 //std::cout << "BTAGSF " << btagJetInfos.size() << " " << btag.weight<BTag1Tight2CustomFilter>(btagJetInfos) << std::endl;
1819 if ( btagJetInfos.size()< 10)
1820 {
1821 btag1T2CSF = btag.weight<BTag1Tight2CustomFilter>(btagJetInfos);
1822 btag2TSF = btag.weight<BTag2TightFilter>(btagJetInfos);
1823 btag1TSF = btag.weight<BTag1TightFilter>(btagJetInfos);
1824 btagA0CSF = btag.weight<BTagAntiMax0CustomFilter>(btagJetInfos);
1825 btagA0TSF = btag.weight<BTagAntiMax0TightFilter>(btagJetInfos);
1826 btag2CSF = btag.weight<BTag2CustomFilter>(btagJetInfos);
1827 btag1TA1C = btag.weight<BTag1TightAndMax1CustomFilter>(btagJetInfos);
1828 }
1829 else
1830 {
1831 std::cout << "WARNING: combinatorics for " << btagJetInfos.size() << " jets is too high (>=10). use SF=1 " << std::endl;
1832 //TODO: revert to random throw for this cases
1833 btag1T2CSF = 1.;
1834 btag2TSF = 1.;
1835 btag1TSF = 1.;
1836 btagA0CSF = 1.;
1837 btagA0TSF = 1.;
1838 btag2CSF = 1.;
1839 btag1TA1C = 1.;
1840 }
1841 }
1842
1843 if(maxBtag > -99999)
1844 {
1845 TopHypo topQuark = TopMassReco::topMass(leptonForTop,bJet,vhCand.V.mets.at(0).p4);
1846 top.mass = topQuark.p4.M();
1847 top.pt = topQuark.p4.Pt();
1848 top.wMass = topQuark.p4W.M();
1849 } else {
1850 top.mass = -99;
1851 top.pt = -99;
1852 top.wMass = -99;
1853 }
1854
1855
1856
1857 //FIXME: too much warnings... figure out why
1858 // gendrcc=aux.genCCDeltaR();
1859 // gendrbb=aux.genBBDeltaR();
1860
1861
1862
1863 genZpt=aux.mcZ.size() > 0 ? aux.mcZ[0].p4.Pt():-99;
1864 genWpt=aux.mcW.size() > 0 ? aux.mcW[0].p4.Pt():-99;
1865
1866 // Z* is status=3 and Nmother=2 (q and qbar)
1867 // Z is status=2 and Ndau=2 (mup and mum)
1868 for (unsigned int i=0; i<aux.mcZ.size(); i++){
1869 if (aux.mcZ[i].status==3) {
1870 genZstar.mass = aux.mcZ[i].p4.M();
1871 genZstar.pt = aux.mcZ[i].p4.Pt();
1872 genZstar.eta = aux.mcZ[i].p4.Eta();
1873 genZstar.phi = aux.mcZ[i].p4.Phi();
1874 genZstar.status = aux.mcZ[i].status;
1875 genZstar.charge = aux.mcZ[i].charge;
1876 if (aux.mcZ[i].momid!=-99) genZstar.momid = aux.mcZ[i].momid ;
1877 }
1878
1879
1880
1881 if (aux.mcZ[i].dauid.size()>1) {
1882 if ( abs(aux.mcZ[i].dauid[0])==13 || abs(aux.mcZ[i].dauid[0])==11 ) {
1883 genZ.mass = aux.mcZ[i].p4.M();
1884 genZ.pt = aux.mcZ[i].p4.Pt();
1885 genZ.eta = aux.mcZ[i].p4.Eta();
1886 genZ.phi = aux.mcZ[i].p4.Phi();
1887 genZ.status = aux.mcZ[i].status;
1888 genZ.charge = aux.mcZ[i].charge;
1889 if ( aux.mcZ[i].momid!=-99) genZ.momid = aux.mcZ[i].momid;
1890 }
1891 }
1892 }
1893
1894
1895
1896 for (unsigned int i=0; i<aux.mcW.size(); i++){
1897 if ( aux.mcW[i].momid==6 && aux.mcW[i].dauid.size()>1 ){
1898 genTop.wdau1mass= aux.mcW[i].dauFourMomentum[0].M();
1899 genTop.wdau1pt= aux.mcW[i].dauFourMomentum[0].Pt();
1900 genTop.wdau1eta= aux.mcW[i].dauFourMomentum[0].Eta();
1901 genTop.wdau1phi= aux.mcW[i].dauFourMomentum[0].Phi();
1902 genTop.wdau1id= aux.mcW[i].dauid[0];
1903
1904 genTop.wdau2mass= aux.mcW[i].dauFourMomentum[1].M();
1905 genTop.wdau2pt= aux.mcW[i].dauFourMomentum[1].Pt();
1906 genTop.wdau2eta= aux.mcW[i].dauFourMomentum[1].Eta();
1907 genTop.wdau2phi= aux.mcW[i].dauFourMomentum[1].Phi();
1908 genTop.wdau2id= aux.mcW[i].dauid[1];
1909
1910 }
1911
1912
1913 if ( aux.mcW[i].momid==-6 && aux.mcW[i].dauid.size()>1 ){
1914 genTbar.wdau1mass= aux.mcW[i].dauFourMomentum[0].M();
1915 genTbar.wdau1pt= aux.mcW[i].dauFourMomentum[0].Pt();
1916 genTbar.wdau1eta= aux.mcW[i].dauFourMomentum[0].Eta();
1917 genTbar.wdau1phi= aux.mcW[i].dauFourMomentum[0].Phi();
1918 genTbar.wdau1id= aux.mcW[i].dauid[0];
1919
1920 genTbar.wdau2mass= aux.mcW[i].dauFourMomentum[1].M();
1921 genTbar.wdau2pt= aux.mcW[i].dauFourMomentum[1].Pt();
1922 genTbar.wdau2eta= aux.mcW[i].dauFourMomentum[1].Eta();
1923 genTbar.wdau2phi= aux.mcW[i].dauFourMomentum[1].Phi();
1924 genTbar.wdau2id= aux.mcW[i].dauid[1];
1925
1926 }
1927
1928 if (aux.mcW[i].status==3 ) {
1929 genWstar.mass = aux.mcW[i].p4.M();
1930 genWstar.pt = aux.mcW[i].p4.Pt();
1931 genWstar.eta = aux.mcW[i].p4.Eta();
1932 genWstar.phi = aux.mcW[i].p4.Phi();
1933 genWstar.status = aux.mcW[i].status;
1934 genWstar.charge = aux.mcW[i].charge;
1935 if ( aux.mcW[i].momid!=-99) genWstar.momid = aux.mcW[i].momid;
1936 }
1937 if (aux.mcW[i].dauid.size()>1 && (abs(aux.mcW[i].dauid[0])==13 || abs(aux.mcW[i].dauid[0])==11 )) {
1938 genW.mass = aux.mcW[i].p4.M();
1939 genW.pt = aux.mcW[i].p4.Pt();
1940 genW.eta = aux.mcW[i].p4.Eta();
1941 genW.phi = aux.mcW[i].p4.Phi();
1942 genW.status = aux.mcW[i].status;
1943 genW.charge = aux.mcW[i].charge;
1944 if (aux.mcW[i].momid!=-99) genW.momid = aux.mcW[i].momid;
1945 }
1946 }
1947 // b coming from Higgs
1948 for (unsigned int i=0; i<aux.mcB.size(); i++){
1949 if (abs(aux.mcB[i].momid)!=5) {
1950 genB.mass = aux.mcB[i].p4.M();
1951 genB.pt = aux.mcB[i].p4.Pt();
1952 genB.eta = aux.mcB[i].p4.Eta();
1953 genB.phi = aux.mcB[i].p4.Phi();
1954 genB.status = aux.mcB[i].status;
1955 genB.charge = aux.mcB[i].charge;
1956 if (aux.mcB[i].momid!=-99) genB.momid = aux.mcB[i].momid;
1957 }
1958
1959 if ( aux.mcB[i].momid==6 ){
1960 genTop.bmass = aux.mcB[i].p4.M();
1961 genTop.bpt = aux.mcB[i].p4.Pt();
1962 genTop.beta = aux.mcB[i].p4.Eta();
1963 genTop.bphi = aux.mcB[i].p4.Phi();
1964 genTop.bstatus = aux.mcB[i].status;
1965
1966 }
1967 }
1968
1969 for (unsigned int i=0; i<aux.mcBbar.size(); i++){
1970 if (abs(aux.mcBbar[i].momid)!=5 ) {
1971 genBbar.mass = aux.mcBbar[i].p4.M();
1972 genBbar.pt = aux.mcBbar[i].p4.Pt();
1973 genBbar.eta = aux.mcBbar[i].p4.Eta();
1974 genBbar.phi = aux.mcBbar[i].p4.Phi();
1975 genBbar.status = aux.mcBbar[i].status;
1976 if (aux.mcBbar[i].momid!=-99) genBbar.momid = aux.mcBbar[i].momid;
1977 }
1978 if ( aux.mcBbar[i].momid==-6 ){
1979 genTbar.bmass = aux.mcBbar[i].p4.M();
1980 genTbar.bpt = aux.mcBbar[i].p4.Pt();
1981 genTbar.beta = aux.mcBbar[i].p4.Eta();
1982 genTbar.bphi = aux.mcBbar[i].p4.Phi();
1983 genTbar.bstatus = aux.mcBbar[i].status;
1984
1985
1986 }
1987
1988 }
1989
1990
1991
1992
1993 if (aux.mcH.size()>0) {
1994 genH.mass = aux.mcH[0].p4.M();
1995 genH.pt = aux.mcH[0].p4.Pt();
1996 genH.eta = aux.mcH[0].p4.Eta();
1997 genH.phi = aux.mcH[0].p4.Phi();
1998 genH.status = aux.mcH[0].status;
1999 genH.charge = aux.mcH[0].charge;
2000 if (aux.mcH[0].momid!=-99) genH.momid = aux.mcH[0].momid;
2001 }
2002
2003
2004
2005
2006 WminusMode=-99;
2007 WplusMode=-99;
2008 for(unsigned int j=0; j< aux.mcW.size();j++)
2009 {
2010 for(unsigned int k=0;k< aux.mcW[j].dauid.size();k++)
2011 {
2012 int idd=abs(aux.mcW[j].dauid[k]);
2013 if(idd==11 || idd==13 || idd==15|| (idd<=5 && idd >=1))
2014 {
2015 if(WminusMode==-99 && aux.mcW[j].charge ==-1) WminusMode = idd;
2016 if(WplusMode==-99 && aux.mcW[j].charge ==+1) WplusMode = idd;
2017 }
2018 }
2019 /*
2020 /// now check if a semileptonic W is also in a bjets....
2021 if ( ( (WminusMode==11 || WminusMode==13 || WminusMode==15 ) || (WplusMode==11 || WplusMode==13 || WplusMode==15 )) && deltaR(vhCand.H.jets[0].p4.Eta(),vhCand.H.jets[0].p4.Phi(), aux.mcW[j].p4.Eta(), aux.mcW[j].p4.Phi())<0.3 ) hJets.isSemiLeptMCtruth[0]=1;
2022 if ( ( (WminusMode==11 || WminusMode==13 || WminusMode==15 ) || (WplusMode==11 || WplusMode==13 || WplusMode==15 )) && deltaR(vhCand.H.jets[1].p4.Eta(),vhCand.H.jets[1].p4.Phi(), aux.mcW[j].p4.Eta(), aux.mcW[j].p4.Phi())<0.3 ) hJets.isSemiLeptMCtruth[1]=1;
2023
2024 for( int j=0; j < naJets && j < MAXJ; j++ )
2025 {
2026 if ((idd==11 || idd==13 || idd==15 ) && deltaR(vhCand.additionalJets[j].p4.Eta(),vhCand.additionalJets[j].p4.Phi(), aux.mcW[j].p4.Eta(), aux.mcW[j].p4.Phi()) <0.3) aJets.isSemiLept[j]=1;
2027
2028 }
2029 */
2030
2031 }
2032 /// Compute pull angle from AK7
2033 if(vhCand.H.HiggsFlag){
2034 if(!fromCandidate){
2035 std::vector<VHbbEvent::SimpleJet> ak7wrt1(iEvent->simpleJets3);
2036 std::vector<VHbbEvent::SimpleJet> ak7wrt2(iEvent->simpleJets3);
2037 if(ak7wrt1.size() > 1){
2038 CompareDeltaR deltaRComparatorJ1(vhCand.H.jets[0].p4);
2039 CompareDeltaR deltaRComparatorJ2(vhCand.H.jets[1].p4);
2040 std::sort( ak7wrt1.begin(),ak7wrt1.end(),deltaRComparatorJ1 );
2041 std::sort( ak7wrt2.begin(),ak7wrt2.end(),deltaRComparatorJ2 );
2042 std::vector<VHbbEvent::SimpleJet> ak7_matched;
2043 // if the matched are different save them
2044 if(ak7wrt1[0].p4.DeltaR(ak7wrt2[0].p4) > 0.1) {
2045 ak7_matched.push_back(ak7wrt1[0]);
2046 ak7_matched.push_back(ak7wrt2[0]);
2047 }
2048 // else look at the second best
2049 else{
2050 // ak7wrt1 is best
2051 if( ak7wrt1[1].p4.DeltaR(vhCand.H.jets[0].p4) < ak7wrt2[1].p4.DeltaR(vhCand.H.jets[1].p4))
2052 {
2053 ak7_matched.push_back(ak7wrt1[1]);
2054 ak7_matched.push_back(ak7wrt2[0]);
2055 }
2056 else
2057 {
2058 ak7_matched.push_back(ak7wrt1[0]);
2059 ak7_matched.push_back(ak7wrt2[1]);
2060 }
2061 }
2062 CompareJetPt ptComparator;
2063 std::sort( ak7_matched.begin(),ak7_matched.end(),ptComparator );
2064 if(ak7_matched[0].p4.DeltaR(vhCand.H.jets[0].p4) < 0.5
2065 and ak7_matched[1].p4.DeltaR(vhCand.H.jets[1].p4) < 0.5)
2066 {
2067 deltaPullAngleAK7 = VHbbCandidateTools::getDeltaTheta(ak7_matched[0],ak7_matched[1]);
2068 deltaPullAngle2AK7 = VHbbCandidateTools::getDeltaTheta(ak7_matched[1],ak7_matched[0]);
2069 }
2070 }
2071 }
2072 }//HiggsFlag
2073
2074 _outTree->Fill();
2075
2076 }// closed event loop
2077
2078 std::cout << "closing the file: " << inputFiles_[iFile] << std::endl;
2079 inFile->Close();
2080 // close input file
2081 } // loop on files
2082
2083
2084 std::cout << "Events: " << ievt <<std::endl;
2085 std::cout << "TotalCount: " << totalcount <<std::endl;
2086
2087
2088
2089 _outFile->cd();
2090
2091 _outTree->Write();
2092 _outFile->Write();
2093 _outFile->Close();
2094 return 0;
2095 }
2096
2097