ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbbAnalysis/VHbbDataFormats/bin/Ntupler.cc
Revision: 1.60
Committed: Sun Oct 30 20:49:10 2011 UTC (13 years, 6 months ago) by arizzi
Content type: text/plain
Branch: MAIN
Changes since 1.59: +10 -4 lines
Log Message:
add new Zee triggers

File Contents

# User Rev Content
1 arizzi 1.1 #include <TH1F.h>
2     #include "PhysicsTools/Utilities/interface/LumiReWeighting.h"
3     #include <TH2F.h>
4     #include <TROOT.h>
5     #include <TFile.h>
6     #include <TTree.h>
7     #include <TSystem.h>
8     #include "DataFormats/FWLite/interface/Event.h"
9     #include "DataFormats/FWLite/interface/Handle.h"
10     #include "FWCore/FWLite/interface/AutoLibraryLoader.h"
11     #include "DataFormats/MuonReco/interface/Muon.h"
12     #include "DataFormats/PatCandidates/interface/Muon.h"
13     #include "PhysicsTools/FWLite/interface/TFileService.h"
14     #include "FWCore/ParameterSet/interface/ProcessDesc.h"
15     #include "FWCore/PythonParameterSet/interface/PythonProcessDesc.h"
16     #include "DataFormats/Math/interface/deltaR.h"
17     #include "DataFormats/Math/interface/deltaPhi.h"
18    
19     #include "DataFormats/FWLite/interface/LuminosityBlock.h"
20     #include "DataFormats/FWLite/interface/Run.h"
21     #include "DataFormats/Luminosity/interface/LumiSummary.h"
22    
23 arizzi 1.4 #include "VHbbAnalysis/VHbbDataFormats/interface/HbbCandidateFinderAlgo.h"
24     #include "VHbbAnalysis/VHbbDataFormats/src/HbbCandidateFinderAlgo.cc"
25 arizzi 1.1
26     #include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEvent.h"
27     #include "VHbbAnalysis/VHbbDataFormats/interface/VHbbEventAuxInfo.h"
28     #include "VHbbAnalysis/VHbbDataFormats/interface/VHbbCandidate.h"
29     #include "VHbbAnalysis/VHbbDataFormats/interface/TriggerReader.h"
30 nmohr 1.5 #include "VHbbAnalysis/VHbbDataFormats/interface/TopMassReco.h"
31 arizzi 1.1
32 bortigno 1.30 //for IVF
33     #include "RecoBTag/SecondaryVertex/interface/SecondaryVertex.h"
34     #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
35     #include "DataFormats/GeometryVector/interface/GlobalVector.h"
36     #include "DataFormats/GeometryVector/interface/VectorUtil.h"
37     #include <DataFormats/GeometrySurface/interface/Surface.h>
38     #include "Math/SMatrix.h"
39 arizzi 1.16
40     //Move class definition to Ntupler.h ?
41     //#include "VHbbAnalysis/VHbbDataFormats/interface/Ntupler.h"
42    
43 bortigno 1.30 //for SimBhad: FIXME
44     //#include "VHbbAnalysis/BAnalysis/interface/SimBHadron.h"
45     //btagging
46 arizzi 1.16 #include "VHbbAnalysis/VHbbDataFormats/interface/BTagWeight.h"
47 bortigno 1.30 //trigger
48 arizzi 1.16 #include "VHbbAnalysis/VHbbDataFormats/interface/TriggerWeight.h"
49    
50 arizzi 1.1 #include <sstream>
51     #include <string>
52    
53     #define MAXJ 30
54     #define MAXL 10
55 bortigno 1.30 #define MAXB 10
56 arizzi 1.1
57 arizzi 1.24 struct CompareDeltaR {
58 bortigno 1.30 CompareDeltaR(TLorentzVector p4dir_): p4dir(p4dir_) {}
59     bool operator()( const VHbbEvent::SimpleJet& j1, const VHbbEvent::SimpleJet& j2 ) const {
60     return j1.p4.DeltaR(p4dir) > j2.p4.DeltaR(p4dir);
61     }
62     TLorentzVector p4dir;
63 arizzi 1.24 };
64 arizzi 1.4
65 bortigno 1.30 const GlobalVector flightDirection(const TVector3 pv, const reco::Vertex &sv){
66     GlobalVector fdir(sv.position().X() - pv.X(),
67     sv.position().Y() - pv.Y(),
68     sv.position().Z() - pv.Z());
69     return fdir;
70     }
71    
72 nmohr 1.7 bool jsonContainsEvent (const std::vector< edm::LuminosityBlockRange > &jsonVec,
73     const edm::EventBase &event)
74     {
75 bortigno 1.30 // if the jsonVec is empty, then no JSON file was provided so all
76     // events should pass
77     if (jsonVec.empty())
78     {
79 nmohr 1.7 return true;
80 bortigno 1.30 }
81     bool (* funcPtr) (edm::LuminosityBlockRange const &,
82     edm::LuminosityBlockID const &) = &edm::contains;
83     edm::LuminosityBlockID lumiID (event.id().run(),
84     event.id().luminosityBlock());
85     std::vector< edm::LuminosityBlockRange >::const_iterator iter =
86     std::find_if (jsonVec.begin(), jsonVec.end(),
87     boost::bind(funcPtr, _1, lumiID) );
88     return jsonVec.end() != iter;
89 nmohr 1.7
90     }
91    
92 arizzi 1.1
93 bortigno 1.30 //FIXME : need to update EDM ntuple with SimBhadron infos
94     // typedef struct
95     // {
96     // void set(const SimBHadron & sbhc, int i){
97     // mass[i] = sbhc.mass();
98     // pt[i] = sbhc.pt();
99     // eta[i] = sbhc.eta();
100     // phi[i] = sbhc.phi();
101     // vtx_x[i] = sbhc.decPosition.x();
102     // vtx_y[i] = sbhc.decPosition.y();
103     // vtx_z[i] = sbhc.decPosition.z();
104     // pdgId[i] = sbhc.pdgId();
105     // status[i] = sbhc.status();
106     // };
107     // void reset(){
108     // for(int i=0; i < MAXB; ++i){
109     // 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;
110     // }
111     // };
112     // float mass[MAXB];
113     // float pt[MAXB];
114     // float eta[MAXB];
115     // float phi[MAXB];
116     // float vtx_x[MAXB];
117     // float vtx_y[MAXB];
118     // float vtx_z[MAXB];
119     // int pdgId[MAXB];
120     // int status[MAXB];
121     // // int quarkStatus[MAXB];
122     // // int brotherStatus[MAXB];
123     // // int otherId[MAXB];
124     // // bool etaOk[MAXB];
125     // // bool simOk[MAXB];
126     // // bool trackOk[MAXB];
127     // // bool cutOk[MAXB];
128     // // bool cutNewOk[MAXB];
129     // // bool mcMatchOk[MAXB];
130     // // bool matchOk[MAXB];
131     // } SimBHadronInfo;
132    
133    
134     typedef struct
135     {
136     void set( const reco::SecondaryVertex & recoSv, const TVector3 recoPv, int isv){
137     pt[isv] = recoSv.p4().Pt();
138     eta[isv] = flightDirection(recoPv,recoSv).eta();
139     phi[isv] = flightDirection(recoPv,recoSv).phi();
140     massBcand[isv] = recoSv.p4().M();
141     massSv[isv] = recoSv.p4().M();
142     dist3D[isv] = recoSv.dist3d().value();
143     distSig3D[isv] = recoSv.dist3d().significance();
144     dist2D[isv] = recoSv.dist2d().value();
145     distSig2D[isv] = recoSv.dist2d().significance();
146     dist3D_norm[isv] = recoSv.dist3d().value()/recoSv.p4().Gamma();
147     };
148     void reset(){
149     for(int i = 0; i < MAXB; ++i){
150     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;
151     }
152     };
153     float massBcand[MAXB];
154     float massSv[MAXB];
155     float pt[MAXB];
156     float eta[MAXB];
157     float phi[MAXB];
158     float dist3D[MAXB];
159     float distSig3D[MAXB];
160     float dist2D[MAXB];
161     float distSig2D[MAXB];
162     float dist3D_norm[MAXB];
163     } IVFInfo;
164 arizzi 1.1
165 bortigno 1.30
166     typedef struct
167     {
168 arizzi 1.33 float mass;
169     float pt;
170     float eta;
171     float phi;
172     float dR;
173     float dPhi;
174     float dEta;
175     } HiggsInfo;
176    
177     typedef struct
178     {
179 bortigno 1.30 float mass; //MT in case of W
180     float pt;
181     float eta;
182     float phi;
183     } TrackInfo;
184 arizzi 1.1
185    
186 arizzi 1.9 struct LeptonInfo
187 bortigno 1.30 {
188     void reset()
189 arizzi 1.1 {
190 arizzi 1.48 for(int i =0; i < MAXL;i++){
191 dlopes 1.55 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;
192 arizzi 1.48 id80[i]=-99; id95[i]=-99; vbtf[i]=-99;
193     }
194 bortigno 1.30 }
195 arizzi 1.1
196 bortigno 1.30 template <class Input> void set(const Input & i, int j,int t)
197     {
198     type[j]=t;
199     pt[j]=i.p4.Pt();
200     mass[j]=i.p4.M();
201     eta[j]=i.p4.Eta();
202     phi[j]=i.p4.Phi();
203     aodCombRelIso[j]=(i.hIso+i.eIso+i.tIso)/i.p4.Pt();
204     pfCombRelIso[j]=(i.pfChaIso+i.pfPhoIso+i.pfNeuIso)/i.p4.Pt();
205     photonIso[j]=i.pfPhoIso;
206     neutralHadIso[j]=i.pfNeuIso;
207     chargedHadIso[j]=i.pfChaIso;
208 dlopes 1.55 chargedPUIso[j]=i.pfChaPUIso;
209 bortigno 1.30 setSpecific(i,j);
210     }
211     template <class Input> void setSpecific(const Input & i, int j)
212     {
213     }
214 arizzi 1.4
215    
216    
217 arizzi 1.1
218 bortigno 1.30 float mass[MAXL]; //MT in case of W
219     float pt[MAXL];
220     float eta[MAXL];
221     float phi[MAXL];
222     float aodCombRelIso[MAXL];
223     float pfCombRelIso[MAXL];
224     float photonIso[MAXL];
225     float neutralHadIso[MAXL];
226     float chargedHadIso[MAXL];
227 dlopes 1.55 float chargedPUIso[MAXL];
228 bortigno 1.30 float particleIso[MAXL];
229     float dxy[MAXL];
230     float dz[MAXL];
231     int type[MAXL];
232     float id80[MAXL];
233     float id95[MAXL];
234 arizzi 1.48 float vbtf[MAXL];
235 bortigno 1.30 };
236 arizzi 1.1
237 bortigno 1.30 template <> void LeptonInfo::setSpecific<VHbbEvent::ElectronInfo>(const VHbbEvent::ElectronInfo & i, int j){
238 arizzi 1.36 id80[j]=i.id80;
239     id95[j]=i.id95;
240 bortigno 1.30 }
241     template <> void LeptonInfo::setSpecific<VHbbEvent::MuonInfo>(const VHbbEvent::MuonInfo & i, int j){
242     dxy[j]=i.ipDb;
243     dz[j]=i.zPVPt;
244 arizzi 1.48 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 &&
245     (i.pfChaIso+i.pfPhoIso+i.pfNeuIso)/i.p4.Pt()<.15 && fabs(i.p4.Eta())<2.4 && i.p4.Pt()>20 ) ;
246    
247 bortigno 1.30 }
248    
249 arizzi 1.4
250 bortigno 1.30 typedef struct
251     {
252     float et;
253     float sumet;
254     float sig;
255     float phi;
256     } METInfo;
257 arizzi 1.1
258 bortigno 1.30 typedef struct
259     {
260     float mht;
261     float ht;
262     float sig;
263     float phi;
264     } MHTInfo;
265 nmohr 1.5
266 bortigno 1.30 typedef struct
267     {
268     float mass;
269     float pt;
270     float wMass;
271     } TopInfo;
272 arizzi 1.1
273 bortigno 1.30 typedef struct
274     {
275     int run;
276     int lumi;
277     int event;
278     int json;
279     } EventInfo;
280 arizzi 1.1
281 bortigno 1.30 typedef struct
282     {
283     void set(const VHbbEvent::SimpleJet & j, int i)
284 arizzi 1.1 {
285 bortigno 1.30 pt[i]=j.p4.Pt();
286     eta[i]=j.p4.Eta();
287     phi[i]=j.p4.Phi();
288 arizzi 1.51 e[i]=j.p4.E();
289 bortigno 1.30 csv[i]=j.csv;
290     numTracksSV[i] = j.vtxMass;
291     vtxMass[i]= j.vtxMass;
292     vtx3dL[i] = j.vtx3dL;
293     vtx3deL[i] = j.vtx3deL;
294     chf[i]=j.chargedHadronEFraction;
295     nhf[i] =j.neutralHadronEFraction;
296     cef[i] =j.chargedEmEFraction;
297     nef[i] =j.neutralEmEFraction;
298     nconstituents[i] =j.nConstituents;
299     nch[i]=j.ntracks;
300    
301     flavour[i]=j.flavour;
302     if(j.bestMCp4.Pt() > 0)
303     {
304     genPt[i]=j.bestMCp4.Pt();
305     genEta[i]=j.bestMCp4.Eta();
306     genPhi[i]=j.bestMCp4.Phi();
307     }
308     JECUnc[i]=j.jecunc;
309 arizzi 1.36 id[i]=jetId(i);
310     }
311     bool jetId(int i)
312     {
313     if(nhf[i] > 0.99) return false;
314     if(nef[i] > 0.99) return false;
315 arizzi 1.37 if(nconstituents[i] <= 1) return false;
316 arizzi 1.47 if(fabs(eta[i])<2.5) {
317 arizzi 1.36 if(cef[i] > 0.99) return false;
318 arizzi 1.37 if(chf[i] == 0) return false;
319 arizzi 1.36 if(nch[i]== 0) return false;
320 arizzi 1.37 }
321 arizzi 1.36 return true;
322 bortigno 1.30 }
323     void reset()
324     {
325     for(int i=0;i<MAXJ;i++) {
326 arizzi 1.51 pt[i]=-99; eta[i]=-99; phi[i]=-99;e[i]=-99;csv[i]=-99; 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; genPt[i]=-99; genEta[i]=-99; genPhi[i]=-99; JECUnc[i]=-99;
327 arizzi 1.1 }
328 bortigno 1.30 }
329     float pt[MAXJ];
330     float eta[MAXJ];
331     float phi[MAXJ];
332 arizzi 1.51 float e[MAXJ];
333 bortigno 1.30 float csv[MAXJ];
334     float cosTheta[MAXJ];
335     int numTracksSV[MAXJ];
336     float chf[MAXJ];
337     float nhf[MAXJ];
338     float cef[MAXJ];
339     float nef[MAXJ];
340     float nch[MAXJ];
341     float nconstituents[MAXJ];
342     float flavour[MAXJ];
343     float genPt[MAXJ];
344     float genEta[MAXJ];
345     float genPhi[MAXJ];
346     float JECUnc[MAXJ];
347     float vtxMass[MAXJ];
348     float vtx3dL [MAXJ];
349     float vtx3deL[MAXJ];
350 arizzi 1.36 bool id[MAXJ];
351 bortigno 1.30 } JetInfo;
352 arizzi 1.1
353     int main(int argc, char* argv[])
354     {
355     gROOT->Reset();
356    
357     TTree *_outTree;
358 bortigno 1.30 IVFInfo IVF;
359     //FIXME
360     // SimBHadronInfo SimBs;
361 arizzi 1.39 float rho,rho25;
362     int nPVs;
363 arizzi 1.9 METInfo MET;
364 arizzi 1.52 METInfo fakeMET;
365 arizzi 1.36 METInfo METnoPU;
366 arizzi 1.58 METInfo METnoPUCh;
367 arizzi 1.9 MHTInfo MHT;
368     TopInfo top;
369     EventInfo EVENT;
370 bortigno 1.30 // JetInfo jet1,jet2, addJet1, addJet2;
371 arizzi 1.9 // lepton1,lepton2;
372     JetInfo hJets, aJets;
373     LeptonInfo vLeptons, aLeptons;
374 arizzi 1.1 int naJets=0, nhJets=0;
375 arizzi 1.33 HiggsInfo H,SVH; //add here the fatjet higgs
376 arizzi 1.9 TrackInfo V;
377     int nvlep=0,nalep=0;
378 arizzi 1.1
379 degrutto 1.54
380    
381     float HVdPhi,HVMass,HMETdPhi,VMt,deltaPullAngle,deltaPullAngleAK7,deltaPullAngle2,deltaPullAngle2AK7,gendrcc,gendrbb, genZpt, genWpt, weightTrig, weightTrigMay,weightTrigV4, weightTrigMET, weightTrigOrMu30, minDeltaPhijetMET, jetPt_minDeltaPhijetMET , PUweight, PUweight2011B;
382    
383    
384 arizzi 1.60 float weightEleRecoAndId,weightEleTrigJetMETPart, weightEleTrigElePart,weightEleTrigEleAugPart;
385 arizzi 1.24
386 bortigno 1.30 //FIXME
387     // float SimBs_dr = -99, SimBs_dEta= -99, SimBs_dPhi = -99, SimBs_Hmass = -99, SimBs_Hpt = -99, SimBs_Heta = -99, SimBs_Hphi = -99;
388 arizzi 1.49 int WplusMode,WminusMode;
389 bortigno 1.30 int Vtype,nSvs,nSimBs,numJets,numBJets,eventFlav;
390     // bool isMET80_CJ80, ispfMHT150, isMET80_2CJ20,isMET65_2CJ20, isJETID,isIsoMu17;
391 arizzi 1.56 bool triggerFlags[500],hbhe,ecalFlag;
392 arizzi 1.60 float btag1T2CSF=1.,btag2TSF=1.,btag1TSF=1.,btagA0CSF=1., btagA0TSF=1., btag2CSF=1., btag1TA1C=1.;
393 arizzi 1.1 // ----------------------------------------------------------------------
394     // First Part:
395     //
396     // * enable the AutoLibraryLoader
397     // * book the histograms of interest
398     // * open the input file
399     // ----------------------------------------------------------------------
400    
401     // load framework libraries
402 bortigno 1.30 gSystem->Load("libFWCoreFWLite");
403 arizzi 1.1 gSystem->Load("libDataFormatsFWLite");
404     AutoLibraryLoader::enable();
405    
406     // parse arguments
407     if ( argc < 2 ) {
408     return 0;
409     }
410    
411 arizzi 1.4 std::vector<VHbbCandidate> * candZlocal = new std::vector<VHbbCandidate>;
412     std::vector<VHbbCandidate> * candWlocal = new std::vector<VHbbCandidate>;
413    
414 arizzi 1.1 // get the python configuration
415     PythonProcessDesc builder(argv[1]);
416     const edm::ParameterSet& in = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("fwliteInput" );
417     const edm::ParameterSet& out = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("fwliteOutput");
418     const edm::ParameterSet& ana = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("Analyzer");
419 bortigno 1.30
420 nmohr 1.7 std::vector<edm::LuminosityBlockRange> jsonVector;
421     if ( in.exists("lumisToProcess") )
422 bortigno 1.30 {
423     std::vector<edm::LuminosityBlockRange> const & lumisTemp =
424     in.getUntrackedParameter<std::vector<edm::LuminosityBlockRange> > ("lumisToProcess");
425 nmohr 1.7 jsonVector.resize( lumisTemp.size() );
426     copy( lumisTemp.begin(), lumisTemp.end(), jsonVector.begin() );
427 bortigno 1.30 }
428 arizzi 1.1
429     // now get each parameter
430     int maxEvents_( in.getParameter<int>("maxEvents") );
431 nmohr 1.29 int skipEvents_( in.getParameter<int>("skipEvents") );
432 arizzi 1.1 unsigned int outputEvery_( in.getParameter<unsigned int>("outputEvery") );
433     std::string outputFile_( out.getParameter<std::string>("fileName" ) );
434     std::vector<std::string> triggers( ana.getParameter<std::vector<std::string> >("triggers") );
435 arizzi 1.6 double btagThr = ana.getParameter<double>("bJetCountThreshold" );
436 arizzi 1.4 bool fromCandidate = ana.getParameter<bool>("readFromCandidates");
437 arizzi 1.17 bool useHighestPtHiggsZ = ana.getParameter<bool>("useHighestPtHiggsZ");
438     bool useHighestPtHiggsW = ana.getParameter<bool>("useHighestPtHiggsW");
439 nmohr 1.31 HbbCandidateFinderAlgo * algoZ = new HbbCandidateFinderAlgo(ana.getParameter<bool>("verbose"), ana.getParameter<double>("jetPtThresholdZ"),useHighestPtHiggsZ);
440     HbbCandidateFinderAlgo * algoW = new HbbCandidateFinderAlgo(ana.getParameter<bool>("verbose"), ana.getParameter<double>("jetPtThresholdW"),useHighestPtHiggsW );
441 arizzi 1.4
442 arizzi 1.16 TriggerWeight triggerWeight(ana);
443 arizzi 1.17 BTagWeight btag(2); // 2 operating points "Custom" = 0.5 and "Tight = 0.898"
444     BTagSampleEfficiency btagEff( ana.getParameter<std::string>("btagEffFileName" ).c_str() );
445 arizzi 1.1
446     std::vector<std::string> inputFiles_( in.getParameter<std::vector<std::string> >("fileNames") );
447 bortigno 1.30 // std::string inputFile( in.getParameter<std::string> ("fileName") );
448 arizzi 1.1
449 arizzi 1.3 std::string PUmcfileName_ = in.getParameter<std::string> ("PUmcfileName") ;
450 degrutto 1.57 std::string PUmcfileName2011B_ = in.getParameter<std::string> ("PUmcfileName2011B") ;
451    
452 arizzi 1.3 std::string PUdatafileName_ = in.getParameter<std::string> ("PUdatafileName") ;
453 degrutto 1.54 std::string PUdatafileName2011B_ = in.getParameter<std::string> ("PUdatafileName2011B") ;
454    
455    
456    
457 arizzi 1.1 bool isMC_( ana.getParameter<bool>("isMC") );
458     TriggerReader trigger(isMC_);
459 arizzi 1.18 TriggerReader patFilters(false);
460 arizzi 1.25
461     edm::LumiReWeighting lumiWeights;
462 degrutto 1.54 edm::LumiReWeighting lumiWeights2011B;
463 arizzi 1.25 if(isMC_)
464 bortigno 1.30 {
465 arizzi 1.34 lumiWeights = edm::LumiReWeighting(PUmcfileName_,PUdatafileName_ , "pileup", "pileup");
466 degrutto 1.54
467 degrutto 1.57 lumiWeights2011B = edm::LumiReWeighting(PUmcfileName2011B_,PUdatafileName2011B_ , "pileup", "pileup");
468     //lumiWeights2011B.weight3D_init(); // generate the weights the fisrt time;
469 degrutto 1.54 lumiWeights2011B.weight3D_init("Weight3D.root");
470    
471    
472 bortigno 1.30 }
473 arizzi 1.1
474 bortigno 1.30 // TFile *_outPUFile = new TFile((outputFile_+"_PU").c_str(), "recreate");
475     TFile *_outFile = new TFile(outputFile_.c_str(), "recreate");
476     TH1F * count = new TH1F("Count","Count", 1,0,2 );
477     TH1F * countWithPU = new TH1F("CountWithPU","CountWithPU", 1,0,2 );
478 arizzi 1.58 TH1F * countWithPU2011B = new TH1F("CountWithPU2011B","CountWithPU2011B", 1,0,2 );
479    
480 arizzi 1.41 TH1F * pu = new TH1F("pileup","",51,-0.5,50.5);
481 arizzi 1.1 _outTree = new TTree("tree", "myTree");
482    
483 arizzi 1.33 _outTree->Branch("H" , &H , "mass/F:pt/F:eta:phi/F:dR/F:dPhi/F:dEta/F");
484 arizzi 1.1 _outTree->Branch("V" , &V , "mass/F:pt/F:eta:phi/F");
485     _outTree->Branch("nhJets" , &nhJets , "nhJets/I");
486     _outTree->Branch("naJets" , &naJets , "naJets/I");
487 arizzi 1.2
488     _outTree->Branch("hJet_pt",hJets.pt ,"pt[nhJets]/F");
489     _outTree->Branch("hJet_eta",hJets.eta ,"eta[nhJets]/F");
490     _outTree->Branch("hJet_phi",hJets.phi ,"phi[nhJets]/F");
491 arizzi 1.51 _outTree->Branch("hJet_e",hJets.e ,"e[nhJets]/F");
492 arizzi 1.2 _outTree->Branch("hJet_csv",hJets.csv ,"csv[nhJets]/F");
493     _outTree->Branch("hJet_cosTheta",hJets.cosTheta ,"cosTheta[nhJets]/F");
494     _outTree->Branch("hJet_numTracksSV",hJets.numTracksSV ,"numTracksSV[nhJets]/I");
495     _outTree->Branch("hJet_chf",hJets.chf ,"chf[nhJets]/F");
496     _outTree->Branch("hJet_nhf",hJets.nhf ,"nhf[nhJets]/F");
497     _outTree->Branch("hJet_cef",hJets.cef ,"cef[nhJets]/F");
498     _outTree->Branch("hJet_nef",hJets.nef ,"nef[nhJets]/F");
499     _outTree->Branch("hJet_nch",hJets.nch ,"nch[nhJets]/F");
500     _outTree->Branch("hJet_nconstituents",hJets.nconstituents ,"nconstituents[nhJets]");
501     _outTree->Branch("hJet_flavour",hJets.flavour ,"flavour[nhJets]/F");
502     _outTree->Branch("hJet_genPt",hJets.genPt ,"genPt[nhJets]/F");
503     _outTree->Branch("hJet_genEta",hJets.genEta ,"genEta[nhJets]/F");
504     _outTree->Branch("hJet_genPhi",hJets.genPhi ,"genPhi[nhJets]/F");
505     _outTree->Branch("hJet_JECUnc",hJets.JECUnc ,"JECUnc[nhJets]/F");
506 arizzi 1.36 _outTree->Branch("hJet_vtxMass",hJets.vtxMass ,"vtxMass[nhJets]/F");
507     _outTree->Branch("hJet_vtx3dL",hJets.vtx3dL ,"vtx3dL[nhJets]/F");
508     _outTree->Branch("hJet_vtx3deL",hJets.vtx3deL ,"vtx3deL[nhJets]/F");
509     _outTree->Branch("hJet_id",hJets.id ,"id[nhJets]/b");
510 arizzi 1.2
511     _outTree->Branch("aJet_pt",aJets.pt ,"pt[naJets]/F");
512     _outTree->Branch("aJet_eta",aJets.eta ,"eta[naJets]/F");
513     _outTree->Branch("aJet_phi",aJets.phi ,"phi[naJets]/F");
514 arizzi 1.51 _outTree->Branch("aJet_e",aJets.e ,"e[naJets]/F");
515 arizzi 1.2 _outTree->Branch("aJet_csv",aJets.csv ,"csv[naJets]/F");
516     _outTree->Branch("aJet_cosTheta",aJets.cosTheta ,"cosTheta[naJets]/F");
517     _outTree->Branch("aJet_numTracksSV",aJets.numTracksSV ,"numTracksSV[naJets]/I");
518     _outTree->Branch("aJet_chf",aJets.chf ,"chf[naJets]/F");
519     _outTree->Branch("aJet_nhf",aJets.nhf ,"nhf[naJets]/F");
520     _outTree->Branch("aJet_cef",aJets.cef ,"cef[naJets]/F");
521     _outTree->Branch("aJet_nef",aJets.nef ,"nef[naJets]/F");
522     _outTree->Branch("aJet_nch",aJets.nch ,"nch[naJets]/F");
523     _outTree->Branch("aJet_nconstituents",aJets.nconstituents ,"nconstituents[naJets]");
524     _outTree->Branch("aJet_flavour",aJets.flavour ,"flavour[naJets]/F");
525     _outTree->Branch("aJet_genPt",aJets.genPt ,"genPt[naJets]/F");
526     _outTree->Branch("aJet_genEta",aJets.genEta ,"genEta[naJets]/F");
527     _outTree->Branch("aJet_genPhi",aJets.genPhi ,"genPhi[naJets]/F");
528     _outTree->Branch("aJet_JECUnc",aJets.JECUnc ,"JECUnc[naJets]/F");
529 arizzi 1.24 _outTree->Branch("aJet_vtxMass",aJets.vtxMass ,"vtxMass[naJets]/F");
530     _outTree->Branch("aJet_vtx3dL",aJets.vtx3dL ,"vtx3dL[naJets]/F");
531     _outTree->Branch("aJet_vtx3deL",aJets.vtx3deL ,"vtx3deL[naJets]/F");
532 arizzi 1.36 _outTree->Branch("aJet_id",aJets.id ,"id[naJets]/b");
533 arizzi 1.1
534     _outTree->Branch("numJets" , &numJets , "numJets/I" );
535     _outTree->Branch("numBJets" , &numBJets , "numBJets/I" );
536     _outTree->Branch("deltaPullAngle", &deltaPullAngle , "deltaPullAngle/F");
537 arizzi 1.52 _outTree->Branch("deltaPullAngle2", &deltaPullAngle2 , "deltaPullAngle2/F");
538 arizzi 1.1 _outTree->Branch("gendrcc" , &gendrcc , "gendrcc/F");
539     _outTree->Branch("gendrbb" , &gendrbb , "gendrbb/F");
540     _outTree->Branch("genZpt" , &genZpt , "genZpt/F");
541     _outTree->Branch("genWpt" , &genWpt , "genWpt/F");
542     _outTree->Branch("weightTrig" , &weightTrig , "weightTrig/F");
543 arizzi 1.38 _outTree->Branch("weightTrigMay" , &weightTrigMay , "weightTrigMay/F");
544     _outTree->Branch("weightTrigV4" , &weightTrigV4 , "weightTrigV4/F");
545 arizzi 1.22 _outTree->Branch("weightTrigMET" , &weightTrigMET , "weightTrigMET/F");
546 arizzi 1.51 _outTree->Branch("weightTrigOrMu30" , &weightTrigOrMu30 , "weightTrigOrMu30/F");
547 arizzi 1.24 _outTree->Branch("weightEleRecoAndId" , &weightEleRecoAndId , "weightEleRecoAndId/F");
548     _outTree->Branch("weightEleTrigJetMETPart" , &weightEleTrigJetMETPart , "weightEleTrigJetMETPart/F");
549     _outTree->Branch("weightEleTrigElePart" , &weightEleTrigElePart , "weightEleTrigElePart/F");
550 arizzi 1.60 _outTree->Branch("weightEleTrigEleAugPart" , &weightEleTrigEleAugPart , "weightEleTrigEleAugPart/F");
551 arizzi 1.24
552    
553 arizzi 1.1 _outTree->Branch("deltaPullAngleAK7", &deltaPullAngleAK7 , "deltaPullAngleAK7/F");
554 arizzi 1.52 _outTree->Branch("deltaPullAngle2AK7", &deltaPullAngle2AK7 , "deltaPullAngle2AK7/F");
555 arizzi 1.1 _outTree->Branch("PUweight", &PUweight , "PUweight/F");
556 degrutto 1.54 _outTree->Branch("PUweight2011B", &PUweight2011B , "PUweight2011B/F");
557 arizzi 1.3 _outTree->Branch("eventFlav", &eventFlav , "eventFlav/I");
558 arizzi 1.1
559    
560    
561    
562     _outTree->Branch("Vtype" , &Vtype , "Vtype/I" );
563     _outTree->Branch("HVdPhi" , &HVdPhi , "HVdPhi/F" );
564 arizzi 1.41 _outTree->Branch("HVMass" , &HVMass , "HVMass/F" );
565 arizzi 1.19 _outTree->Branch("HMETdPhi" , &HMETdPhi , "HMETdPhi/F" );
566 arizzi 1.1 _outTree->Branch("VMt" , &VMt , "VMt/F" );
567    
568 arizzi 1.9 _outTree->Branch("nvlep" , &nvlep , "nvlep/I");
569     _outTree->Branch("nalep" , &nalep , "nalep/I");
570 arizzi 1.3
571 arizzi 1.9 _outTree->Branch("vLepton_mass",vLeptons.mass ,"mass[nvlep]/F");
572     _outTree->Branch("vLepton_pt",vLeptons.pt ,"pt[nvlep]/F");
573     _outTree->Branch("vLepton_eta",vLeptons.eta ,"eta[nvlep]");
574     _outTree->Branch("vLepton_phi",vLeptons.phi ,"phi[nvlep]/F");
575     _outTree->Branch("vLepton_aodCombRelIso",vLeptons.aodCombRelIso ,"aodCombRelIso[nvlep]/F");
576     _outTree->Branch("vLepton_pfCombRelIso",vLeptons.pfCombRelIso ,"pfCombRelIso[nvlep]/F");
577     _outTree->Branch("vLepton_photonIso",vLeptons.photonIso ,"photonIso[nvlep]/F");
578     _outTree->Branch("vLepton_neutralHadIso",vLeptons.neutralHadIso ,"neutralHadIso[nvlep]/F");
579     _outTree->Branch("vLepton_chargedHadIso",vLeptons.chargedHadIso ,"chargedHadIso[nvlep]/F");
580 dlopes 1.55 _outTree->Branch("vLepton_chargedPUIso",vLeptons.chargedPUIso ,"chargedPUIso[nvlep]/F");
581 arizzi 1.9 _outTree->Branch("vLepton_particleIso",vLeptons.particleIso ,"particleIso[nvlep]/F");
582     _outTree->Branch("vLepton_dxy",vLeptons.dxy ,"dxy[nvlep]/F");
583     _outTree->Branch("vLepton_dz",vLeptons.dz ,"dz[nvlep]/F");
584     _outTree->Branch("vLepton_type",vLeptons.type ,"type[nvlep]/I");
585 arizzi 1.27 _outTree->Branch("vLepton_id80",vLeptons.id80 ,"id80[nvlep]/F");
586     _outTree->Branch("vLepton_id95",vLeptons.id95 ,"id95[nvlep]/F");
587 arizzi 1.48 _outTree->Branch("vLepton_vbtf",vLeptons.vbtf ,"vbtf[nvlep]/F");
588 arizzi 1.9
589     _outTree->Branch("aLepton_mass",aLeptons.mass ,"mass[nalep]/F");
590     _outTree->Branch("aLepton_pt",aLeptons.pt ,"pt[nalep]/F");
591     _outTree->Branch("aLepton_eta",aLeptons.eta ,"eta[nalep]");
592     _outTree->Branch("aLepton_phi",aLeptons.phi ,"phi[nalep]/F");
593     _outTree->Branch("aLepton_aodCombRelIso",aLeptons.aodCombRelIso ,"aodCombRelIso[nalep]/F");
594     _outTree->Branch("aLepton_pfCombRelIso",aLeptons.pfCombRelIso ,"pfCombRelIso[nalep]/F");
595     _outTree->Branch("aLepton_photonIso",aLeptons.photonIso ,"photonIso[nalep]/F");
596     _outTree->Branch("aLepton_neutralHadIso",aLeptons.neutralHadIso ,"neutralHadIso[nalep]/F");
597     _outTree->Branch("aLepton_chargedHadIso",aLeptons.chargedHadIso ,"chargedHadIso[nalep]/F");
598 dlopes 1.55 _outTree->Branch("aLepton_chargedPUIso",aLeptons.chargedPUIso ,"chargedPUIso[nalep]/F");
599 arizzi 1.9 _outTree->Branch("aLepton_particleIso",aLeptons.particleIso ,"particleIso[nalep]/F");
600     _outTree->Branch("aLepton_dxy",aLeptons.dxy ,"dxy[nalep]/F");
601     _outTree->Branch("aLepton_dz",aLeptons.dz ,"dz[nalep]/F");
602     _outTree->Branch("aLepton_type",aLeptons.type ,"type[nalep]/I");
603 arizzi 1.48 _outTree->Branch("aLepton_id80",aLeptons.id80 ,"id80[nalep]/F");
604     _outTree->Branch("aLepton_id95",aLeptons.id95 ,"id95[nalep]/F");
605     _outTree->Branch("aLepton_vbtf",aLeptons.vbtf ,"vbtf[nalep]/F");
606 arizzi 1.9
607 arizzi 1.8 _outTree->Branch("top" , &top , "mass/F:pt/F:wMass/F");
608 arizzi 1.49 _outTree->Branch("WplusMode" , &WplusMode , "WplusMode/I");
609     _outTree->Branch("WminusMode" , &WminusMode , "WminusMode/I");
610 arizzi 1.1
611 bortigno 1.30 //IVF
612     _outTree->Branch("nSvs",&nSvs ,"nSvs/I");
613 arizzi 1.33 _outTree->Branch("Sv_massBCand", &IVF.massBcand,"massBcand[nSvs]/F");
614     _outTree->Branch("Sv_massSv", &IVF.massSv,"massSv[nSvs]/F");
615     _outTree->Branch("Sv_pt", &IVF.pt,"pt[nSvs]/F");
616     _outTree->Branch("Sv_eta", &IVF.eta,"eta[nSvs]/F");
617     _outTree->Branch("Sv_phi", &IVF.phi,"phi[nSvs]/F");
618     _outTree->Branch("Sv_dist3D", &IVF.dist3D,"dist3D[nSvs]/F");
619     _outTree->Branch("Sv_dist2D", &IVF.dist2D,"dist2D[nSvs]/F");
620     _outTree->Branch("Sv_distSim2D", &IVF.distSig2D,"distSig2D[nSvs]/F");
621     _outTree->Branch("Sv_distSig3D", &IVF.distSig3D,"distSig3D[nSvs]/F");
622     _outTree->Branch("Sv_dist3D_norm", &IVF.dist3D_norm,"dist3D_norm[nSvs]/F");
623 bortigno 1.30 //IVF higgs candidate
624 arizzi 1.33 _outTree->Branch("SVH" , &SVH , "mass/F:pt/F:eta:phi/F:dR/F:dPhi/F:dEta/F");
625 bortigno 1.30
626    
627     //FIXME : need to update the EDMNtuple with BHadron infos
628     // //SimBHadron
629     // _outTree->Branch("nSimBs",&nSimBs ,"nSimBs/I");
630     // _outTree->Branch("SimBs_mass", &SimBs.mass,"mass[nSvs]/F");
631     // _outTree->Branch("SimBs_pt", &SimBs.pt,"pt[nSvs]/F");
632     // _outTree->Branch("SimBs_eta", &SimBs.eta,"eta[nSvs]/F");
633     // _outTree->Branch("SimBs_phi", &SimBs.phi,"phi[nSvs]/F");
634     // _outTree->Branch("SimBs_vtx_x", &SimBs.vtx_x,"vtx_x[nSvs]/F");
635     // _outTree->Branch("SimBs_vtx_y", &SimBs.vtx_y,"vtx_y[nSvs]/F");
636     // _outTree->Branch("SimBs_vtx_z", &SimBs.vtx_z,"vtx_z[nSvs]/F");
637     // _outTree->Branch("SimBs_pdgId", &SimBs.pdgId,"pdgId[nSvs]/F");
638     // _outTree->Branch("SimBs_status", &SimBs.status,"status[nSvs]/F");
639     // //SimBHadron Higgs Candidate
640     // _outTree->Branch("SimBs_dr", &SimBs_dr,"SimBs_dr/F");
641     // _outTree->Branch("SimBs_dPhi", &SimBs_dPhi,"SimBs_dPhi/F");
642     // _outTree->Branch("SimBs_dEta", &SimBs_dEta,"SimBs_dEta/F");
643     // _outTree->Branch("SimBs_Hmass", &SimBs_Hmass,"SimBs_Hmass/F");
644     // _outTree->Branch("SimBs_Hpt", &SimBs_Hpt,"SimBs_Hpt/F");
645     // _outTree->Branch("SimBs_Heta", &SimBs_Heta,"SimBs_Heta/F");
646     // _outTree->Branch("SimBs_Hphi", &SimBs_Hphi,"SimBs_Hphi/F");
647    
648    
649 arizzi 1.36 _outTree->Branch("rho" , &rho , "rho/F");
650     _outTree->Branch("rho25" , &rho25 , "rho25/F");
651     _outTree->Branch("nPVs" , &nPVs , "nPVs/I");
652     _outTree->Branch("METnoPU" , &METnoPU , "et/F:sumet:sig/F:phi/F");
653 arizzi 1.58 _outTree->Branch("METnoPUCh" , &METnoPUCh , "et/F:sumet:sig/F:phi/F");
654 arizzi 1.1 _outTree->Branch("MET" , &MET , "et/F:sumet:sig/F:phi/F");
655 arizzi 1.52 _outTree->Branch("fakeMET" , &fakeMET , "et/F:sumet:sig/F:phi/F");
656 arizzi 1.1 _outTree->Branch("MHT" , &MHT , "mht/F:ht:sig/F:phi/F");
657     _outTree->Branch("minDeltaPhijetMET" , &minDeltaPhijetMET , "minDeltaPhijetMET/F");
658     _outTree->Branch("jetPt_minDeltaPhijetMET" , &jetPt_minDeltaPhijetMET , "jetPt_minDeltaPhijetMET/F");
659    
660 bortigno 1.30 std::stringstream s;
661     s << "triggerFlags[" << triggers.size() << "]/b";
662     _outTree->Branch("triggerFlags", triggerFlags, s.str().c_str());
663 nmohr 1.7
664 arizzi 1.17
665 bortigno 1.30 _outTree->Branch("EVENT" , &EVENT , "run/I:lumi/I:event/I:json/I");
666     _outTree->Branch("hbhe" , &hbhe , "hbhe/b");
667 arizzi 1.56 _outTree->Branch("ecalFlag" , &ecalFlag , "ecalFlag/b");
668 bortigno 1.30 _outTree->Branch("btag1TSF" , &btag1TSF , "btag1TSF/F");
669     _outTree->Branch("btag2TSF" , &btag2TSF , "btag2TSF/F");
670     _outTree->Branch("btag1T2CSF" , &btag1T2CSF , "btag1T2CSF/F");
671 arizzi 1.44 _outTree->Branch("btag2CSF" , &btag2CSF , "btag2CSF/F");
672     _outTree->Branch("btagA0CSF" , &btagA0CSF , "btagA0CSF/F");
673 arizzi 1.60 _outTree->Branch("btagA0TSF" , &btagA0TSF , "btagA0TSF/F");
674 arizzi 1.51 _outTree->Branch("btag1TA1C" , &btag1TA1C , "btag1TA1C/F");
675 arizzi 1.1
676 bortigno 1.30 int ievt=0;
677     int totalcount=0;
678 arizzi 1.1
679 bortigno 1.30 // TFile* inFile = new TFile(inputFile.c_str(), "read");
680 arizzi 1.1 for(unsigned int iFile=0; iFile<inputFiles_.size(); ++iFile) {
681     std::cout << iFile << std::endl;
682     TFile* inFile = TFile::Open(inputFiles_[iFile].c_str());
683     if(inFile==0) continue;
684    
685 bortigno 1.30 // loop the events
686 arizzi 1.1
687     fwlite::Event ev(inFile);
688 nmohr 1.29 for(ev.toBegin(); !ev.atEnd() ; ++ev, ++ievt)
689 arizzi 1.1 {
690 nmohr 1.29 if (ievt <= skipEvents_) continue;
691     if (maxEvents_ >= 0){
692     if (ievt > maxEvents_ + skipEvents_) break;
693     };
694 arizzi 1.16 count->Fill(1.);
695 arizzi 1.1
696 arizzi 1.25 fwlite::Handle< VHbbEventAuxInfo > vhbbAuxHandle;
697     vhbbAuxHandle.getByLabel(ev,"HbbAnalyzerNew");
698     const VHbbEventAuxInfo & aux = *vhbbAuxHandle.product();
699     PUweight=1.;
700 degrutto 1.54 PUweight2011B=1.;
701 arizzi 1.1 if(isMC_){
702 arizzi 1.41
703 degrutto 1.54 // PU weights // Run2011A
704 bortigno 1.30 std::map<int, unsigned int>::const_iterator puit = aux.puInfo.pus.find(0);
705 arizzi 1.42 PUweight = lumiWeights.weight( puit->second );
706 arizzi 1.41 pu->Fill(puit->second);
707 degrutto 1.54 // PU weight Run2011B
708     // PU weight Run2011B
709     std::map<int, unsigned int>::const_iterator puit0 = aux.puInfo.pus.find(0);
710     std::map<int, unsigned int>::const_iterator puitm1 = aux.puInfo.pus.find(-1);
711     std::map<int, unsigned int>::const_iterator puitp1 = aux.puInfo.pus.find(+1);
712    
713     PUweight2011B = lumiWeights2011B.weight3D( puitm1->second, puit0->second,puitp1->second);
714    
715 bortigno 1.30 }
716 arizzi 1.49 countWithPU->Fill(1,PUweight);
717 arizzi 1.58 countWithPU2011B->Fill(1,PUweight2011B);
718 nmohr 1.7
719 bortigno 1.30 //Write event info
720     EVENT.run = ev.id().run();
721     EVENT.lumi = ev.id().luminosityBlock();
722     EVENT.event = ev.id().event();
723     EVENT.json = jsonContainsEvent (jsonVector, ev);
724    
725     //FIXME : need to update EDM ntuple with BHadron infos
726     // // simBHadrons
727     // fwlite::Handle<SimBHadronCollection> SBHC;
728     // SBHC.getByLabel(ev, "bhadrons");
729     // const SimBHadronCollection sbhc = *(SBHC.product());
730    
731     const std::vector<VHbbCandidate> * candZ ;
732     const std::vector<VHbbCandidate> * candW ;
733     const VHbbEvent * iEvent =0;
734     if(fromCandidate)
735     {
736     fwlite::Handle< std::vector<VHbbCandidate> > vhbbCandHandleZ;
737     vhbbCandHandleZ.getByLabel(ev,"hbbBestCSVPt20Candidates");
738     candZ = vhbbCandHandleZ.product();
739    
740     fwlite::Handle< std::vector<VHbbCandidate> > vhbbCandHandle;
741     vhbbCandHandle.getByLabel(ev,"hbbHighestPtHiggsPt30Candidates");
742     candW = vhbbCandHandle.product();
743     }
744     else
745     {
746     candZlocal->clear();
747     candWlocal->clear();
748     fwlite::Handle< VHbbEvent > vhbbHandle;
749     vhbbHandle.getByLabel(ev,"HbbAnalyzerNew");
750     iEvent = vhbbHandle.product();
751     algoZ->run(vhbbHandle.product(),*candZlocal);
752     algoW->run(vhbbHandle.product(),*candWlocal);
753     candZ= candZlocal;
754     candW= candWlocal;
755    
756     /* for(size_t m=0;m<iEvent->muInfo.size();m++)
757     {
758    
759     if( fabs(iEvent->muInfo[m].p4.Pt()-28.118684) < 0.0001 ||
760     fabs(iEvent->muInfo[m].p4.Pt()-34.853199) < 0.0001 )
761     {
762     std::cout << "FOUND " << iEvent->muInfo[m].p4.Pt() << " " << EVENT.event << " " << candW->size() << " " << candZ->size() << std::endl;
763     }
764     }
765     */
766 arizzi 1.4
767 bortigno 1.30 }
768 arizzi 1.1
769 arizzi 1.3 const std::vector<VHbbCandidate> * cand = candZ;
770 arizzi 1.1
771    
772    
773 bortigno 1.30 /* fwlite::Handle< VHbbEvent > vhbbHandle;
774     vhbbHandle.getByLabel(ev,"HbbAnalyzerNew");
775     const VHbbEvent iEvent = *vhbbHandle.product();
776     */
777 arizzi 1.1
778 bortigno 1.30 // std::clog << "Filling tree "<< std::endl;
779     bool isW=false;
780 arizzi 1.21
781 arizzi 1.17
782 bortigno 1.30 if(cand->size() == 0 or cand->at(0).H.jets.size() < 2) continue;
783     if(cand->size() > 1 )
784 arizzi 1.3 {
785 bortigno 1.30 std::cout << "MULTIPLE CANDIDATES: " << cand->size() << std::endl;
786 arizzi 1.3 }
787 bortigno 1.30 if(cand->at(0).candidateType == VHbbCandidate::Wmun || cand->at(0).candidateType == VHbbCandidate::Wen ) { cand=candW; isW=true; }
788     if(cand->size() == 0)
789 arizzi 1.3 {
790 bortigno 1.30 // std::cout << "W event loss due to tigther cuts" << std::endl;
791 arizzi 1.3 continue;
792     }
793 arizzi 1.35
794 degrutto 1.54
795 arizzi 1.35 // secondary vtxs
796     fwlite::Handle<std::vector<reco::Vertex> > SVC;
797 arizzi 1.52 SVC.getByLabel(ev,"bcandidates");
798 arizzi 1.35 const std::vector<reco::Vertex> svc = *(SVC.product());
799    
800 bortigno 1.30 const VHbbCandidate & vhCand = cand->at(0);
801     patFilters.setEvent(&ev,"VH");
802     hbhe = patFilters.accept("hbhe");
803 arizzi 1.56 ecalFlag = patFilters.accept("ecalFilter");
804 arizzi 1.18
805 bortigno 1.30 trigger.setEvent(&ev);
806     for(size_t j=0;j < triggers.size();j++)
807 arizzi 1.4 triggerFlags[j]=trigger.accept(triggers[j]);
808 arizzi 1.39
809 bortigno 1.30 eventFlav=0;
810     if(aux.mcBbar.size() > 0 || aux.mcB.size() > 0) eventFlav=5;
811     else if(aux.mcC.size() > 0) eventFlav=4;
812 arizzi 1.8
813 arizzi 1.53 Vtype = vhCand.candidateType;
814 arizzi 1.3
815 bortigno 1.30 H.mass = vhCand.H.p4.M();
816     H.pt = vhCand.H.p4.Pt();
817     H.eta = vhCand.H.p4.Eta();
818     H.phi = vhCand.H.p4.Phi();
819     V.mass = vhCand.V.p4.M();
820 arizzi 1.36 if(isW) V.mass = vhCand.Mt();
821 bortigno 1.30 V.pt = vhCand.V.p4.Pt();
822     V.eta = vhCand.V.p4.Eta();
823     V.phi = vhCand.V.p4.Phi();
824     nhJets=2;
825     hJets.set(vhCand.H.jets[0],0);
826     hJets.set(vhCand.H.jets[1],1);
827     aJets.reset();
828     naJets=vhCand.additionalJets.size();
829     numBJets=0;
830     if(vhCand.H.jets[0].csv> btagThr) numBJets++;
831     if(vhCand.H.jets[1].csv> btagThr) numBJets++;
832     for( int j=0; j < naJets && j < MAXJ; j++ )
833 arizzi 1.6 {
834 bortigno 1.30 aJets.set(vhCand.additionalJets[j],j);
835     if(vhCand.additionalJets[j].csv> btagThr) numBJets++;
836 arizzi 1.6 }
837 bortigno 1.30 numJets = vhCand.additionalJets.size()+2;
838 arizzi 1.33 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());
839     H.dPhi = deltaPhi(vhCand.H.jets[0].p4.Phi(),vhCand.H.jets[1].p4.Phi());
840     H.dEta= TMath::Abs( vhCand.H.jets[0].p4.Eta() - vhCand.H.jets[1].p4.Eta() );
841 bortigno 1.30 HVdPhi = fabs( deltaPhi(vhCand.H.p4.Phi(),vhCand.V.p4.Phi()) ) ;
842 arizzi 1.41 HVMass = (vhCand.H.p4 + vhCand.V.p4).M() ;
843 bortigno 1.30 HMETdPhi = fabs( deltaPhi(vhCand.H.p4.Phi(),vhCand.V.mets.at(0).p4.Phi()) ) ;
844     VMt = vhCand.Mt() ;
845 arizzi 1.52 //eltaPullAngle = vhCand.H.deltaTheta;
846     deltaPullAngle = VHbbCandidateTools::getDeltaTheta(vhCand.H.jets[0],vhCand.H.jets[1]);
847     deltaPullAngle2 = VHbbCandidateTools::getDeltaTheta(vhCand.H.jets[1],vhCand.H.jets[0]);
848    
849 bortigno 1.30
850     hJets.cosTheta[0]= vhCand.H.helicities[0];
851     hJets.cosTheta[1]= vhCand.H.helicities[1];
852 arizzi 1.36 // METInfo calomet; METInfo tcmet; METInfo pfmet; METInfo mht; METInfo metNoPU
853 bortigno 1.30 MET.et = vhCand.V.mets.at(0).p4.Pt();
854     MET.phi = vhCand.V.mets.at(0).p4.Phi();
855     MET.sumet = vhCand.V.mets.at(0).sumEt;
856     MET.sig = vhCand.V.mets.at(0).metSig;
857 arizzi 1.24
858 arizzi 1.52
859     fakeMET.sumet = 0;
860     fakeMET.sig = 0;
861     fakeMET.et = 0;
862     fakeMET.phi = 0;
863     if( Vtype == VHbbCandidate::Zmumu) {
864     TVector3 mu1 = vhCand.V.muons[0].p4.Vect();
865     TVector3 mu2 = vhCand.V.muons[1].p4.Vect();
866     // Not needed with PFMET
867     // mu1.SetMag( mu1.Mag() - vhCand.V.muons[0].emEnergy - vhCand.V.muons[0].hadEnergy);
868     // mu2.SetMag( mu2.Mag() - vhCand.V.muons[1].emEnergy - vhCand.V.muons[1].hadEnergy);
869     TVector3 sum = vhCand.V.mets.at(0).p4.Vect() + mu1 + mu2;
870     fakeMET.et = sum.Pt();
871     fakeMET.phi = sum.Phi();
872     fakeMET.sumet = vhCand.V.mets.at(0).sumEt - mu1.Pt() - mu2.Pt();
873     }
874    
875 arizzi 1.36 METnoPU.et = iEvent->metNoPU.p4.Pt();
876     METnoPU.phi = iEvent->metNoPU.p4.Phi();
877     METnoPU.sumet = iEvent->metNoPU.sumEt;
878     METnoPU.sig = iEvent->metNoPU.metSig;
879 arizzi 1.59 METnoPUCh.et = iEvent->metCh.p4.Pt();
880     METnoPUCh.phi = iEvent->metCh.p4.Phi();
881     METnoPUCh.sumet = iEvent->metCh.sumEt;
882     METnoPUCh.sig = iEvent->metCh.metSig;
883 arizzi 1.36
884 arizzi 1.58 METnoPUCh.et = iEvent->metCh.p4.Pt();
885     METnoPUCh.phi = iEvent->metCh.p4.Phi();
886     METnoPUCh.sumet = iEvent->metCh.sumEt;
887     METnoPUCh.sig = iEvent->metCh.metSig;
888    
889 arizzi 1.36 rho = aux.puInfo.rho;
890     rho25 = aux.puInfo.rho25;
891     nPVs=aux.pvInfo.nVertices;
892    
893 arizzi 1.24 if(!fromCandidate) {
894     MHT.mht = iEvent->mht.p4.Pt();
895     MHT.phi = iEvent->mht.p4.Phi();
896     MHT.ht = iEvent->mht.sumEt;
897     MHT.sig = iEvent->mht.metSig;
898 bortigno 1.30 }
899    
900 arizzi 1.24
901    
902 bortigno 1.30 //Secondary Vertices
903     IVF.reset();
904     nSvs = svc.size();
905     const TVector3 recoPv = aux.pvInfo.firstPVInPT2;
906     const math::XYZPoint myPv(recoPv);
907     //look here for Matrix filling info http://project-mathlibs.web.cern.ch/project-mathlibs/sw/html/SMatrixDoc.html
908     std::vector<double> fillMatrix(6);
909     for (int i = 0; i<6; ++i) fillMatrix[i] = 0.;
910     fillMatrix[0] = TMath::Power(0.002,2);
911     fillMatrix[2] = TMath::Power(0.002,2);
912     fillMatrix[5] = TMath::Power(0.002,2);
913     const ROOT::Math::SMatrix<double, 3, 3, ROOT::Math::MatRepSym<double, 3> > myFakeMatrixError(fillMatrix.begin(),fillMatrix.end());
914     const reco::Vertex recoVtxPv(myPv, myFakeMatrixError);
915     for( int j=0; j < nSvs && j < MAXB; ++j ) {
916     const GlobalVector flightDir = flightDirection(recoPv,svc[j]);
917     reco::SecondaryVertex recoSv(recoVtxPv, svc[j], flightDir ,true);
918     IVF.set( recoSv, recoPv ,j);
919     }
920     if(nSvs > 2){
921     TLorentzVector BCands_H1, BCands_H2, BCands_H;
922     BCands_H1.SetPtEtaPhiM(IVF.pt[0], IVF.eta[0], IVF.phi[0], IVF.massBcand[0]);
923     BCands_H2.SetPtEtaPhiM(IVF.pt[1], IVF.eta[1], IVF.phi[1], IVF.massBcand[1]);
924     BCands_H = BCands_H1 + BCands_H2;
925 arizzi 1.33 SVH.dR = deltaR(IVF.eta[0], IVF.phi[0], IVF.eta[1], IVF.phi[1] );
926     SVH.dPhi = deltaPhi(IVF.phi[0], IVF.phi[1] );
927     SVH.dEta = TMath::Abs(IVF.eta[0] - IVF.eta[1] );
928     SVH.mass = BCands_H.M();
929     SVH.pt = BCands_H.Pt();
930     SVH.eta = BCands_H.Eta();
931     SVH.phi = BCands_H.Phi();
932 bortigno 1.30 }
933    
934     //FIXME : need to update EDM ntuple with simBhadron info
935     // //SimBhadron
936     // SimBs.reset();
937     // nSimBs = sbhc.size();
938     // for( int j=0; j < nSimBs && j < MAXB; ++j )
939     // SimBs.set( sbhc.at(j), j);
940     // if(nSimBs > 2){
941     // TLorentzVector SimBs_H1, SimBs_H2, SimBs_H;
942     // SimBs_H1.SetPtEtaPhiM(SimBs.pt[0], SimBs.eta[0], SimBs.phi[0], SimBs.mass[0]);
943     // SimBs_H2.SetPtEtaPhiM(SimBs.pt[1], SimBs.eta[1], SimBs.phi[1], SimBs.mass[1]);
944     // SimBs_H = SimBs_H1 + SimBs_H2;
945     // SimBs_dr = deltaR(SimBs.eta[0], SimBs.phi[0], SimBs.eta[1], SimBs.phi[1] );
946     // SimBs_dPhi = deltaPhi(SimBs.phi[0], SimBs.phi[1] );
947     // SimBs_dEta = TMath::Abs(SimBs.eta[0] - SimBs.eta[1] );
948     // SimBs_Hmass = SimBs_H.M();
949     // SimBs_Hpt = SimBs_H.Pt();
950     // SimBs_Heta = SimBs_H.Eta();
951     // SimBs_Hphi = SimBs_H.Phi();
952     // }
953    
954    
955     //Loop on jets
956     double maxBtag=-99999;
957     minDeltaPhijetMET = 999;
958     TLorentzVector bJet;
959     std::vector<std::vector<BTagWeight::JetInfo> > btagJetInfos;
960     std::vector<float> jet30eta;
961     std::vector<float> jet30pt;
962     if(fromCandidate)
963     {
964     //Loop on Higgs Jets
965     for(unsigned int j=0; j < vhCand.H.jets.size(); j++ ){
966     if (vhCand.H.jets[j].csv > maxBtag) { bJet=vhCand.H.jets[j].p4 ; maxBtag =vhCand.H.jets[j].csv; }
967 arizzi 1.39 if (fabs(deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.H.jets[j].p4.Phi())) < minDeltaPhijetMET)
968 bortigno 1.30 {
969 arizzi 1.39 minDeltaPhijetMET=fabs(deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.H.jets[j].p4.Phi()));
970 arizzi 1.24 jetPt_minDeltaPhijetMET=vhCand.H.jets[j].p4.Pt();
971 bortigno 1.30 }
972     btagJetInfos.push_back(btagEff.jetInfo(vhCand.H.jets[j]));
973     }
974     //Loop on Additional Jets
975     for(unsigned int j=0; j < vhCand.additionalJets.size(); j++ ){
976     if (vhCand.additionalJets[j].csv > maxBtag) { bJet=vhCand.additionalJets[j].p4 ; maxBtag =vhCand.additionalJets[j].csv; }
977 arizzi 1.45 /* if (fabs(deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.additionalJets[j].p4.Phi())) < minDeltaPhijetMET)
978 bortigno 1.30 {
979 arizzi 1.39 minDeltaPhijetMET=fabs(deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.additionalJets[j].p4.Phi()));
980 arizzi 1.24 jetPt_minDeltaPhijetMET=vhCand.additionalJets[j].p4.Pt();
981 arizzi 1.45 }*/
982 bortigno 1.30 if( ( isW && ! useHighestPtHiggsW ) || ( ! isW && ! useHighestPtHiggsZ ) ) // btag SF computed using only H-jets if best-H made with dijetPt rather than best CSV
983     {
984     if(vhCand.additionalJets[j].p4.Pt() > 20)
985     btagJetInfos.push_back(btagEff.jetInfo(vhCand.additionalJets[j]));
986     }
987     }
988     }
989     else
990     {
991 arizzi 1.45 //Loop on all jets
992 bortigno 1.30 for(unsigned int j=0; j < iEvent->simpleJets2.size(); j++ ){
993     if (iEvent->simpleJets2[j].csv > maxBtag) { bJet=iEvent->simpleJets2[j].p4 ; maxBtag =iEvent->simpleJets2[j].csv; }
994 arizzi 1.47 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)
995 arizzi 1.46 {
996     minDeltaPhijetMET=fabs(deltaPhi( vhCand.V.mets.at(0).p4.Phi(), iEvent->simpleJets2[j].p4.Phi()));
997     jetPt_minDeltaPhijetMET=iEvent->simpleJets2[j].p4.Pt();
998     }
999 bortigno 1.30 if(iEvent->simpleJets2[j].p4.Pt() > 30)
1000     {
1001     jet30eta.push_back(iEvent->simpleJets2[j].p4.Eta());
1002     jet30pt.push_back(iEvent->simpleJets2[j].p4.Pt());
1003     }
1004 arizzi 1.24
1005 bortigno 1.30 //For events made with highest CSV, all jets in the event should be taken into account for "tagging" SF (anti tagging is a mess)
1006     // 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
1007     if( ( isW && ! useHighestPtHiggsW ) || ( ! isW && ! useHighestPtHiggsZ ) )
1008     {
1009 arizzi 1.58 if(iEvent->simpleJets2[j].p4.Pt() > 20 && fabs(iEvent->simpleJets2[j].p4.Eta()) < 2.5)
1010 bortigno 1.30 btagJetInfos.push_back(btagEff.jetInfo(iEvent->simpleJets2[j]));
1011     }
1012     }
1013    
1014 arizzi 1.45 //Loop on Higgs jets
1015    
1016     for(unsigned int j=0; j < vhCand.H.jets.size(); j++ ) {
1017    
1018 bortigno 1.30 //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
1019     // by a criteria (pt of the dijet) that is not btag SF dependent
1020     if(!( ( isW && ! useHighestPtHiggsW ) || ( ! isW && ! useHighestPtHiggsZ ) )) {
1021 arizzi 1.45 btagJetInfos.push_back(btagEff.jetInfo(vhCand.H.jets[j]));
1022 bortigno 1.30 }
1023 arizzi 1.45 }
1024 arizzi 1.24
1025 bortigno 1.30 }
1026     vLeptons.reset();
1027     weightTrig = 1.; // better to default to 1
1028 arizzi 1.38 weightTrigMay = -1.;
1029     weightTrigV4 = -1.;
1030 arizzi 1.51 weightTrigOrMu30 = 1.;
1031 bortigno 1.30 TLorentzVector leptonForTop;
1032     size_t firstAddMu=0;
1033     size_t firstAddEle=0;
1034     if(Vtype == VHbbCandidate::Zmumu ){
1035     vLeptons.set(vhCand.V.muons[0],0,13);
1036     vLeptons.set(vhCand.V.muons[1],1,13);
1037     float cweightID = triggerWeight.scaleMuID(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleMuID(vLeptons.pt[1],vLeptons.eta[1]) ;
1038     float weightTrig1 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[0],vLeptons.eta[0]);
1039     float weightTrig2 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[1],vLeptons.eta[1]);
1040     float cweightTrig = weightTrig1 + weightTrig2 - weightTrig1*weightTrig2;
1041     weightTrig = cweightID * cweightTrig;
1042     nvlep=2;
1043     firstAddMu=2;
1044     }
1045     if( Vtype == VHbbCandidate::Zee ){
1046     vLeptons.set(vhCand.V.electrons[0],0,11);
1047     vLeptons.set(vhCand.V.electrons[1],1,11);
1048     nvlep=2;
1049     firstAddEle=2;
1050     std::vector<float> pt,eta;
1051     pt.push_back(vLeptons.pt[0]); eta.push_back(vLeptons.eta[0]);
1052     pt.push_back(vLeptons.pt[1]); eta.push_back(vLeptons.eta[1]);
1053     weightEleRecoAndId=triggerWeight.scaleID95Ele(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleRecoEle(vLeptons.pt[0],vLeptons.eta[0]) *
1054     triggerWeight.scaleID95Ele(vLeptons.pt[1],vLeptons.eta[1]) * triggerWeight.scaleRecoEle(vLeptons.pt[1],vLeptons.eta[1]);
1055     weightEleTrigElePart = triggerWeight.scaleDoubleEle17Ele8(pt,eta);
1056 arizzi 1.60 weightEleTrigEleAugPart = triggerWeight.scaleDoubleEle17Ele8Aug(pt,eta);
1057     weightTrig = (weightEleTrigElePart*1.14+weightEleTrigEleAugPart*0.98 )/2.12 * weightEleRecoAndId;
1058 arizzi 1.24
1059    
1060    
1061 bortigno 1.30 }
1062     if(Vtype == VHbbCandidate::Wmun ){
1063     leptonForTop=vhCand.V.muons[0].p4;
1064     vLeptons.set(vhCand.V.muons[0],0,13);
1065     float cweightID = triggerWeight.scaleMuID(vLeptons.pt[0],vLeptons.eta[0]);
1066     float weightTrig1 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[0],vLeptons.eta[0]);
1067     float cweightTrig = weightTrig1;
1068     weightTrig = cweightID * cweightTrig;
1069 arizzi 1.51 float weightTrig1OrMu30 = triggerWeight.scaleMuOr30IsoHLT(vLeptons.pt[0],vLeptons.eta[0]);
1070     weightTrigOrMu30 = cweightID*(cweightTrig+weightTrig1OrMu30);
1071 bortigno 1.30 nvlep=1;
1072     firstAddMu=1;
1073     }
1074     if( Vtype == VHbbCandidate::Wen ){
1075     leptonForTop=vhCand.V.electrons[0].p4;
1076     vLeptons.set(vhCand.V.electrons[0],0,11);
1077     nvlep=1;
1078     firstAddEle=1;
1079 arizzi 1.38 weightTrigMay = triggerWeight.scaleSingleEleMay(vLeptons.pt[0],vLeptons.eta[0]);
1080     weightTrigV4 = triggerWeight.scaleSingleEleV4(vLeptons.pt[0],vLeptons.eta[0]);
1081 bortigno 1.30 weightEleRecoAndId=triggerWeight.scaleID80Ele(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleRecoEle(vLeptons.pt[0],vLeptons.eta[0]);
1082 arizzi 1.59 weightEleTrigJetMETPart=triggerWeight.scaleJet30Jet25(jet30pt,jet30eta)*triggerWeight.scalePFMHTEle(MET.et);
1083 arizzi 1.38 weightEleTrigElePart= weightTrigV4; //this is for debugging only, checking only the V4 part
1084 bortigno 1.30
1085 arizzi 1.38 weightTrigMay*=weightEleRecoAndId;
1086     weightTrigV4*=weightEleRecoAndId;
1087     weightTrigV4*=weightEleTrigJetMETPart;
1088 arizzi 1.60 // weightTrig = weightTrigMay * 0.187 + weightTrigV4 * (1.-0.187); //FIXME: use proper lumi if we reload 2.fb
1089     weightTrig = (weightTrigMay * 0.215 + weightTrigV4 * 1.915)/ 2.13; //FIXME: use proper lumi if we reload 2.fb
1090 bortigno 1.30 }
1091     if( Vtype == VHbbCandidate::Znn ){
1092     nvlep=0;
1093     float weightTrig1 = triggerWeight.scaleMetHLT(vhCand.V.mets.at(0).p4.Pt());
1094     weightTrig = weightTrig1;
1095     weightTrigMET = weightTrig1;
1096     }
1097 arizzi 1.38
1098     if(weightTrigMay < 0) weightTrigMay=weightTrig;
1099     if(weightTrigV4 < 0) weightTrigV4=weightTrig;
1100 arizzi 1.40 if(!isMC_)
1101     {
1102     weightTrig = 1.;
1103     weightTrigMay = 1.;
1104     weightTrigV4 = 1.;
1105    
1106     }
1107 bortigno 1.30 aLeptons.reset();
1108     nalep=0;
1109     if(fromCandidate)
1110 arizzi 1.21 {
1111     for(size_t j=firstAddMu;j< vhCand.V.muons.size();j++) aLeptons.set(vhCand.V.muons[j],nalep++,13);
1112     for(size_t j=firstAddEle;j< vhCand.V.electrons.size();j++) aLeptons.set(vhCand.V.electrons[j],nalep++,11);
1113     }
1114 bortigno 1.30 else
1115 arizzi 1.21 {
1116 bortigno 1.30 for(size_t j=0;j< iEvent->muInfo.size();j++)
1117     {
1118 arizzi 1.21 if((j!= vhCand.V.firstLepton && j!= vhCand.V.secondLepton) || ((Vtype != VHbbCandidate::Wmun ) && (Vtype != VHbbCandidate::Zmumu )) )
1119 bortigno 1.30 aLeptons.set(iEvent->muInfo[j],nalep++,13);
1120     }
1121     for(size_t j=0;j< iEvent->eleInfo.size();j++)
1122     {
1123 arizzi 1.21 if((j!= vhCand.V.firstLepton && j!= vhCand.V.secondLepton) || ((Vtype != VHbbCandidate::Wen ) && (Vtype != VHbbCandidate::Zee )))
1124 bortigno 1.30 aLeptons.set(iEvent->eleInfo[j],nalep++,11);
1125     }
1126 arizzi 1.8
1127 arizzi 1.21 }
1128 arizzi 1.8
1129 arizzi 1.24
1130 bortigno 1.30 if(isMC_)
1131     {
1132     //std::cout << "BTAGSF " << btagJetInfos.size() << " " << btag.weight<BTag1Tight2CustomFilter>(btagJetInfos) << std::endl;
1133     if ( btagJetInfos.size()< 10)
1134     {
1135     btag1T2CSF = btag.weight<BTag1Tight2CustomFilter>(btagJetInfos);
1136     btag2TSF = btag.weight<BTag2TightFilter>(btagJetInfos);
1137     btag1TSF = btag.weight<BTag1TightFilter>(btagJetInfos);
1138 arizzi 1.44 btagA0CSF = btag.weight<BTagAntiMax0CustomFilter>(btagJetInfos);
1139 arizzi 1.60 btagA0TSF = btag.weight<BTagAntiMax0TightFilter>(btagJetInfos);
1140 arizzi 1.44 btag2CSF = btag.weight<BTag2CustomFilter>(btagJetInfos);
1141 arizzi 1.51 btag1TA1C = btag.weight<BTag1TightAndMax1CustomFilter>(btagJetInfos);
1142 bortigno 1.30 }
1143     else
1144     {
1145     std::cout << "WARNING: combinatorics for " << btagJetInfos.size() << " jets is too high (>=10). use SF=1 " << std::endl;
1146     //TODO: revert to random throw for this cases
1147     btag1T2CSF = 1.;
1148     btag2TSF = 1.;
1149     btag1TSF = 1.;
1150 arizzi 1.51 btagA0CSF = 1.;
1151 arizzi 1.60 btagA0TSF = 1.;
1152 arizzi 1.51 btag2CSF = 1.;
1153     btag1TA1C = 1.;
1154 bortigno 1.30 }
1155     }
1156 arizzi 1.17
1157 bortigno 1.30 if(maxBtag > -99999)
1158 arizzi 1.8 {
1159 bortigno 1.30 TopHypo topQuark = TopMassReco::topMass(leptonForTop,bJet,vhCand.V.mets.at(0).p4);
1160     top.mass = topQuark.p4.M();
1161     top.pt = topQuark.p4.Pt();
1162     top.wMass = topQuark.p4W.M();
1163 arizzi 1.8 } else {
1164 bortigno 1.30 top.mass = -99;
1165     top.pt = -99;
1166     top.wMass = -99;
1167     }
1168 arizzi 1.8
1169    
1170    
1171 bortigno 1.30 //FIXME: too much warnings... figure out why
1172     // gendrcc=aux.genCCDeltaR();
1173     // gendrbb=aux.genBBDeltaR();
1174     genZpt=aux.mcZ.size() > 0 ? aux.mcZ[0].p4.Pt():-99;
1175     genWpt=aux.mcW.size() > 0 ? aux.mcW[0].p4.Pt():-99;
1176 arizzi 1.49 WminusMode=-99;
1177     WplusMode=-99;
1178 arizzi 1.50 for(unsigned int j=0; j< aux.mcW.size();j++)
1179 arizzi 1.49 {
1180     for(unsigned int k=0;k< aux.mcW[j].dauid.size();k++)
1181     {
1182     int idd=abs(aux.mcW[j].dauid[k]);
1183     if(idd==11 || idd==13 || idd==15|| (idd<=5 && idd >=1))
1184     {
1185     if(WminusMode==-99 && aux.mcW[j].charge ==-1) WminusMode = idd;
1186     if(WplusMode==-99 && aux.mcW[j].charge ==-1) WplusMode = idd;
1187     }
1188     }
1189 arizzi 1.24
1190 arizzi 1.49
1191     }
1192 arizzi 1.25
1193 arizzi 1.32 /// Compute pull angle from AK7
1194 arizzi 1.25 if(!fromCandidate){
1195 arizzi 1.32 std::vector<VHbbEvent::SimpleJet> ak7wrt1(iEvent->simpleJets3);
1196     std::vector<VHbbEvent::SimpleJet> ak7wrt2(iEvent->simpleJets3);
1197     if(ak7wrt1.size() > 1){
1198 arizzi 1.24 CompareDeltaR deltaRComparatorJ1(vhCand.H.jets[0].p4);
1199     CompareDeltaR deltaRComparatorJ2(vhCand.H.jets[1].p4);
1200 arizzi 1.32 std::sort( ak7wrt1.begin(),ak7wrt1.end(),deltaRComparatorJ1 );
1201     std::sort( ak7wrt2.begin(),ak7wrt2.end(),deltaRComparatorJ2 );
1202 arizzi 1.24 std::vector<VHbbEvent::SimpleJet> ak7_matched;
1203 arizzi 1.32 // if the matched are different save them
1204     if(ak7wrt1[0].p4.DeltaR(ak7wrt2[0].p4) > 0.1) {
1205     ak7_matched.push_back(ak7wrt1[0]);
1206     ak7_matched.push_back(ak7wrt2[0]);
1207     }
1208     // else look at the second best
1209     else{
1210     // ak7wrt1 is best
1211     if( ak7wrt1[1].p4.DeltaR(vhCand.H.jets[0].p4) < ak7wrt2[1].p4.DeltaR(vhCand.H.jets[1].p4))
1212     {
1213     ak7_matched.push_back(ak7wrt1[1]);
1214     ak7_matched.push_back(ak7wrt2[0]);
1215     }
1216     else
1217     {
1218     ak7_matched.push_back(ak7wrt1[0]);
1219     ak7_matched.push_back(ak7wrt2[1]);
1220     }
1221     }
1222 arizzi 1.24 CompareJetPt ptComparator;
1223     std::sort( ak7_matched.begin(),ak7_matched.end(),ptComparator );
1224     if(ak7_matched[0].p4.DeltaR(vhCand.H.jets[0].p4) < 0.5
1225 arizzi 1.32 and ak7_matched[1].p4.DeltaR(vhCand.H.jets[1].p4) < 0.5)
1226     {
1227     deltaPullAngleAK7 = VHbbCandidateTools::getDeltaTheta(ak7_matched[0],ak7_matched[1]);
1228 arizzi 1.52 deltaPullAngle2AK7 = VHbbCandidateTools::getDeltaTheta(ak7_matched[1],ak7_matched[0]);
1229 arizzi 1.32 }
1230 arizzi 1.24 }
1231     }
1232 arizzi 1.3
1233 bortigno 1.30 _outTree->Fill();
1234 arizzi 1.4
1235 bortigno 1.30 }// closed event loop
1236 arizzi 1.1
1237 bortigno 1.30 std::cout << "closing the file: " << inputFiles_[iFile] << std::endl;
1238     inFile->Close();
1239     // close input file
1240     } // loop on files
1241 arizzi 1.1
1242    
1243 bortigno 1.30 std::cout << "Events: " << ievt <<std::endl;
1244     std::cout << "TotalCount: " << totalcount <<std::endl;
1245 arizzi 1.1
1246    
1247    
1248 bortigno 1.30 _outFile->cd();
1249 arizzi 1.1
1250 bortigno 1.30 _outTree->Write();
1251     _outFile->Write();
1252     _outFile->Close();
1253     return 0;
1254 arizzi 1.1 }
1255    
1256