ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbbAnalysis/VHbbDataFormats/bin/Ntupler.cc
Revision: 1.111.2.3
Committed: Wed Nov 28 03:11:54 2012 UTC (12 years, 5 months ago) by ntran
Content type: text/plain
Branch: hbbsubstructDevPostHCP
CVS Tags: hbbsubstructDev_8
Changes since 1.111.2.2: +4 -2 lines
Log Message:
add qjets

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