ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbbAnalysis/VHbbDataFormats/bin/Ntupler.cc
Revision: 1.95
Committed: Wed Jun 6 16:03:18 2012 UTC (12 years, 11 months ago) by arizzi
Content type: text/plain
Branch: MAIN
CVS Tags: Step2ForV32_v2, Step2ForV32_v1
Changes since 1.94: +4 -1 lines
Log Message:
add electron presel variable

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