ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbbAnalysis/VHbbDataFormats/bin/Ntupler.cc
Revision: 1.49
Committed: Tue Oct 11 18:08:47 2011 UTC (13 years, 7 months ago) by arizzi
Content type: text/plain
Branch: MAIN
Changes since 1.48: +19 -2 lines
Log Message: 
add W decay mode (to be tested) and fix the COuntWithPU that was badly wrong Fill(PUweight) instead of Fill(1,PUweight). anyhow as the PU has been fixed the standard Count histo can be used in existing ntuples

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     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; particleIso[i]=-99; dxy[i]=-99; dz[i]=-99; type[i]=-99;
192     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     setSpecific(i,j);
209     }
210     template <class Input> void setSpecific(const Input & i, int j)
211     {
212     }
213 arizzi 1.4
214    
215    
216 arizzi 1.1
217 bortigno 1.30 float mass[MAXL]; //MT in case of W
218     float pt[MAXL];
219     float eta[MAXL];
220     float phi[MAXL];
221     float aodCombRelIso[MAXL];
222     float pfCombRelIso[MAXL];
223     float photonIso[MAXL];
224     float neutralHadIso[MAXL];
225     float chargedHadIso[MAXL];
226     float particleIso[MAXL];
227     float dxy[MAXL];
228     float dz[MAXL];
229     int type[MAXL];
230     float id80[MAXL];
231     float id95[MAXL];
232 arizzi 1.48 float vbtf[MAXL];
233 bortigno 1.30 };
234 arizzi 1.1
235 bortigno 1.30 template <> void LeptonInfo::setSpecific<VHbbEvent::ElectronInfo>(const VHbbEvent::ElectronInfo & i, int j){
236 arizzi 1.36 id80[j]=i.id80;
237     id95[j]=i.id95;
238 bortigno 1.30 }
239     template <> void LeptonInfo::setSpecific<VHbbEvent::MuonInfo>(const VHbbEvent::MuonInfo & i, int j){
240     dxy[j]=i.ipDb;
241     dz[j]=i.zPVPt;
242 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 &&
243     (i.pfChaIso+i.pfPhoIso+i.pfNeuIso)/i.p4.Pt()<.15 && fabs(i.p4.Eta())<2.4 && i.p4.Pt()>20 ) ;
244    
245 bortigno 1.30 }
246    
247 arizzi 1.4
248 bortigno 1.30 typedef struct
249     {
250     float et;
251     float sumet;
252     float sig;
253     float phi;
254     } METInfo;
255 arizzi 1.1
256 bortigno 1.30 typedef struct
257     {
258     float mht;
259     float ht;
260     float sig;
261     float phi;
262     } MHTInfo;
263 nmohr 1.5
264 bortigno 1.30 typedef struct
265     {
266     float mass;
267     float pt;
268     float wMass;
269     } TopInfo;
270 arizzi 1.1
271 bortigno 1.30 typedef struct
272     {
273     int run;
274     int lumi;
275     int event;
276     int json;
277     } EventInfo;
278 arizzi 1.1
279 bortigno 1.30 typedef struct
280     {
281     void set(const VHbbEvent::SimpleJet & j, int i)
282 arizzi 1.1 {
283 bortigno 1.30 pt[i]=j.p4.Pt();
284     eta[i]=j.p4.Eta();
285     phi[i]=j.p4.Phi();
286     csv[i]=j.csv;
287     numTracksSV[i] = j.vtxMass;
288     vtxMass[i]= j.vtxMass;
289     vtx3dL[i] = j.vtx3dL;
290     vtx3deL[i] = j.vtx3deL;
291     chf[i]=j.chargedHadronEFraction;
292     nhf[i] =j.neutralHadronEFraction;
293     cef[i] =j.chargedEmEFraction;
294     nef[i] =j.neutralEmEFraction;
295     nconstituents[i] =j.nConstituents;
296     nch[i]=j.ntracks;
297    
298     flavour[i]=j.flavour;
299     if(j.bestMCp4.Pt() > 0)
300     {
301     genPt[i]=j.bestMCp4.Pt();
302     genEta[i]=j.bestMCp4.Eta();
303     genPhi[i]=j.bestMCp4.Phi();
304     }
305     JECUnc[i]=j.jecunc;
306 arizzi 1.36 id[i]=jetId(i);
307     }
308     bool jetId(int i)
309     {
310     if(nhf[i] > 0.99) return false;
311     if(nef[i] > 0.99) return false;
312 arizzi 1.37 if(nconstituents[i] <= 1) return false;
313 arizzi 1.47 if(fabs(eta[i])<2.5) {
314 arizzi 1.36 if(cef[i] > 0.99) return false;
315 arizzi 1.37 if(chf[i] == 0) return false;
316 arizzi 1.36 if(nch[i]== 0) return false;
317 arizzi 1.37 }
318 arizzi 1.36 return true;
319 bortigno 1.30 }
320     void reset()
321     {
322     for(int i=0;i<MAXJ;i++) {
323     pt[i]=-99; eta[i]=-99; phi[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;
324 arizzi 1.1 }
325 bortigno 1.30 }
326     float pt[MAXJ];
327     float eta[MAXJ];
328     float phi[MAXJ];
329     float csv[MAXJ];
330     float cosTheta[MAXJ];
331     int numTracksSV[MAXJ];
332     float chf[MAXJ];
333     float nhf[MAXJ];
334     float cef[MAXJ];
335     float nef[MAXJ];
336     float nch[MAXJ];
337     float nconstituents[MAXJ];
338     float flavour[MAXJ];
339     float genPt[MAXJ];
340     float genEta[MAXJ];
341     float genPhi[MAXJ];
342     float JECUnc[MAXJ];
343     float vtxMass[MAXJ];
344     float vtx3dL [MAXJ];
345     float vtx3deL[MAXJ];
346 arizzi 1.36 bool id[MAXJ];
347 bortigno 1.30 } JetInfo;
348 arizzi 1.1
349     int main(int argc, char* argv[])
350     {
351     gROOT->Reset();
352    
353     TTree *_outTree;
354 bortigno 1.30 IVFInfo IVF;
355     //FIXME
356     // SimBHadronInfo SimBs;
357 arizzi 1.39 float rho,rho25;
358     int nPVs;
359 arizzi 1.9 METInfo MET;
360 arizzi 1.36 METInfo METnoPU;
361 arizzi 1.9 MHTInfo MHT;
362     TopInfo top;
363     EventInfo EVENT;
364 bortigno 1.30 // JetInfo jet1,jet2, addJet1, addJet2;
365 arizzi 1.9 // lepton1,lepton2;
366     JetInfo hJets, aJets;
367     LeptonInfo vLeptons, aLeptons;
368 arizzi 1.1 int naJets=0, nhJets=0;
369 arizzi 1.33 HiggsInfo H,SVH; //add here the fatjet higgs
370 arizzi 1.9 TrackInfo V;
371     int nvlep=0,nalep=0;
372 arizzi 1.1
373 arizzi 1.41 float HVdPhi,HVMass,HMETdPhi,VMt,deltaPullAngle,deltaPullAngleAK7,gendrcc,gendrbb, genZpt, genWpt, weightTrig, weightTrigMay,weightTrigV4, weightTrigMET, minDeltaPhijetMET, jetPt_minDeltaPhijetMET , PUweight;
374 bortigno 1.30 float weightEleRecoAndId,weightEleTrigJetMETPart, weightEleTrigElePart;
375 arizzi 1.24
376 bortigno 1.30 //FIXME
377     // float SimBs_dr = -99, SimBs_dEta= -99, SimBs_dPhi = -99, SimBs_Hmass = -99, SimBs_Hpt = -99, SimBs_Heta = -99, SimBs_Hphi = -99;
378 arizzi 1.49 int WplusMode,WminusMode;
379 bortigno 1.30 int Vtype,nSvs,nSimBs,numJets,numBJets,eventFlav;
380     // bool isMET80_CJ80, ispfMHT150, isMET80_2CJ20,isMET65_2CJ20, isJETID,isIsoMu17;
381     bool triggerFlags[500],hbhe;
382 arizzi 1.44 float btag1T2CSF=1.,btag2TSF=1.,btag1TSF=1.,btagA0CSF=1., btag2CSF=1.;
383 arizzi 1.1 // ----------------------------------------------------------------------
384     // First Part:
385     //
386     // * enable the AutoLibraryLoader
387     // * book the histograms of interest
388     // * open the input file
389     // ----------------------------------------------------------------------
390    
391     // load framework libraries
392 bortigno 1.30 gSystem->Load("libFWCoreFWLite");
393 arizzi 1.1 gSystem->Load("libDataFormatsFWLite");
394     AutoLibraryLoader::enable();
395    
396     // parse arguments
397     if ( argc < 2 ) {
398     return 0;
399     }
400    
401 arizzi 1.4 std::vector<VHbbCandidate> * candZlocal = new std::vector<VHbbCandidate>;
402     std::vector<VHbbCandidate> * candWlocal = new std::vector<VHbbCandidate>;
403    
404 arizzi 1.1 // get the python configuration
405     PythonProcessDesc builder(argv[1]);
406     const edm::ParameterSet& in = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("fwliteInput" );
407     const edm::ParameterSet& out = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("fwliteOutput");
408     const edm::ParameterSet& ana = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("Analyzer");
409 bortigno 1.30
410 nmohr 1.7 std::vector<edm::LuminosityBlockRange> jsonVector;
411     if ( in.exists("lumisToProcess") )
412 bortigno 1.30 {
413     std::vector<edm::LuminosityBlockRange> const & lumisTemp =
414     in.getUntrackedParameter<std::vector<edm::LuminosityBlockRange> > ("lumisToProcess");
415 nmohr 1.7 jsonVector.resize( lumisTemp.size() );
416     copy( lumisTemp.begin(), lumisTemp.end(), jsonVector.begin() );
417 bortigno 1.30 }
418 arizzi 1.1
419     // now get each parameter
420     int maxEvents_( in.getParameter<int>("maxEvents") );
421 nmohr 1.29 int skipEvents_( in.getParameter<int>("skipEvents") );
422 arizzi 1.1 unsigned int outputEvery_( in.getParameter<unsigned int>("outputEvery") );
423     std::string outputFile_( out.getParameter<std::string>("fileName" ) );
424     std::vector<std::string> triggers( ana.getParameter<std::vector<std::string> >("triggers") );
425 arizzi 1.6 double btagThr = ana.getParameter<double>("bJetCountThreshold" );
426 arizzi 1.4 bool fromCandidate = ana.getParameter<bool>("readFromCandidates");
427 arizzi 1.17 bool useHighestPtHiggsZ = ana.getParameter<bool>("useHighestPtHiggsZ");
428     bool useHighestPtHiggsW = ana.getParameter<bool>("useHighestPtHiggsW");
429 nmohr 1.31 HbbCandidateFinderAlgo * algoZ = new HbbCandidateFinderAlgo(ana.getParameter<bool>("verbose"), ana.getParameter<double>("jetPtThresholdZ"),useHighestPtHiggsZ);
430     HbbCandidateFinderAlgo * algoW = new HbbCandidateFinderAlgo(ana.getParameter<bool>("verbose"), ana.getParameter<double>("jetPtThresholdW"),useHighestPtHiggsW );
431 arizzi 1.4
432 arizzi 1.16 TriggerWeight triggerWeight(ana);
433 arizzi 1.17 BTagWeight btag(2); // 2 operating points "Custom" = 0.5 and "Tight = 0.898"
434     BTagSampleEfficiency btagEff( ana.getParameter<std::string>("btagEffFileName" ).c_str() );
435 arizzi 1.1
436     std::vector<std::string> inputFiles_( in.getParameter<std::vector<std::string> >("fileNames") );
437 bortigno 1.30 // std::string inputFile( in.getParameter<std::string> ("fileName") );
438 arizzi 1.1
439 arizzi 1.3 std::string PUmcfileName_ = in.getParameter<std::string> ("PUmcfileName") ;
440     std::string PUdatafileName_ = in.getParameter<std::string> ("PUdatafileName") ;
441 arizzi 1.25
442 arizzi 1.1 bool isMC_( ana.getParameter<bool>("isMC") );
443     TriggerReader trigger(isMC_);
444 arizzi 1.18 TriggerReader patFilters(false);
445 arizzi 1.25
446     edm::LumiReWeighting lumiWeights;
447     if(isMC_)
448 bortigno 1.30 {
449 arizzi 1.34 lumiWeights = edm::LumiReWeighting(PUmcfileName_,PUdatafileName_ , "pileup", "pileup");
450 bortigno 1.30 }
451 arizzi 1.1
452 bortigno 1.30 // TFile *_outPUFile = new TFile((outputFile_+"_PU").c_str(), "recreate");
453     TFile *_outFile = new TFile(outputFile_.c_str(), "recreate");
454     TH1F * count = new TH1F("Count","Count", 1,0,2 );
455     TH1F * countWithPU = new TH1F("CountWithPU","CountWithPU", 1,0,2 );
456 arizzi 1.41 TH1F * pu = new TH1F("pileup","",51,-0.5,50.5);
457 arizzi 1.1 _outTree = new TTree("tree", "myTree");
458    
459 arizzi 1.33 _outTree->Branch("H" , &H , "mass/F:pt/F:eta:phi/F:dR/F:dPhi/F:dEta/F");
460 arizzi 1.1 _outTree->Branch("V" , &V , "mass/F:pt/F:eta:phi/F");
461     _outTree->Branch("nhJets" , &nhJets , "nhJets/I");
462     _outTree->Branch("naJets" , &naJets , "naJets/I");
463 arizzi 1.2
464     _outTree->Branch("hJet_pt",hJets.pt ,"pt[nhJets]/F");
465     _outTree->Branch("hJet_eta",hJets.eta ,"eta[nhJets]/F");
466     _outTree->Branch("hJet_phi",hJets.phi ,"phi[nhJets]/F");
467     _outTree->Branch("hJet_csv",hJets.csv ,"csv[nhJets]/F");
468     _outTree->Branch("hJet_cosTheta",hJets.cosTheta ,"cosTheta[nhJets]/F");
469     _outTree->Branch("hJet_numTracksSV",hJets.numTracksSV ,"numTracksSV[nhJets]/I");
470     _outTree->Branch("hJet_chf",hJets.chf ,"chf[nhJets]/F");
471     _outTree->Branch("hJet_nhf",hJets.nhf ,"nhf[nhJets]/F");
472     _outTree->Branch("hJet_cef",hJets.cef ,"cef[nhJets]/F");
473     _outTree->Branch("hJet_nef",hJets.nef ,"nef[nhJets]/F");
474     _outTree->Branch("hJet_nch",hJets.nch ,"nch[nhJets]/F");
475     _outTree->Branch("hJet_nconstituents",hJets.nconstituents ,"nconstituents[nhJets]");
476     _outTree->Branch("hJet_flavour",hJets.flavour ,"flavour[nhJets]/F");
477     _outTree->Branch("hJet_genPt",hJets.genPt ,"genPt[nhJets]/F");
478     _outTree->Branch("hJet_genEta",hJets.genEta ,"genEta[nhJets]/F");
479     _outTree->Branch("hJet_genPhi",hJets.genPhi ,"genPhi[nhJets]/F");
480     _outTree->Branch("hJet_JECUnc",hJets.JECUnc ,"JECUnc[nhJets]/F");
481 arizzi 1.36 _outTree->Branch("hJet_vtxMass",hJets.vtxMass ,"vtxMass[nhJets]/F");
482     _outTree->Branch("hJet_vtx3dL",hJets.vtx3dL ,"vtx3dL[nhJets]/F");
483     _outTree->Branch("hJet_vtx3deL",hJets.vtx3deL ,"vtx3deL[nhJets]/F");
484     _outTree->Branch("hJet_id",hJets.id ,"id[nhJets]/b");
485 arizzi 1.2
486     _outTree->Branch("aJet_pt",aJets.pt ,"pt[naJets]/F");
487     _outTree->Branch("aJet_eta",aJets.eta ,"eta[naJets]/F");
488     _outTree->Branch("aJet_phi",aJets.phi ,"phi[naJets]/F");
489     _outTree->Branch("aJet_csv",aJets.csv ,"csv[naJets]/F");
490     _outTree->Branch("aJet_cosTheta",aJets.cosTheta ,"cosTheta[naJets]/F");
491     _outTree->Branch("aJet_numTracksSV",aJets.numTracksSV ,"numTracksSV[naJets]/I");
492     _outTree->Branch("aJet_chf",aJets.chf ,"chf[naJets]/F");
493     _outTree->Branch("aJet_nhf",aJets.nhf ,"nhf[naJets]/F");
494     _outTree->Branch("aJet_cef",aJets.cef ,"cef[naJets]/F");
495     _outTree->Branch("aJet_nef",aJets.nef ,"nef[naJets]/F");
496     _outTree->Branch("aJet_nch",aJets.nch ,"nch[naJets]/F");
497     _outTree->Branch("aJet_nconstituents",aJets.nconstituents ,"nconstituents[naJets]");
498     _outTree->Branch("aJet_flavour",aJets.flavour ,"flavour[naJets]/F");
499     _outTree->Branch("aJet_genPt",aJets.genPt ,"genPt[naJets]/F");
500     _outTree->Branch("aJet_genEta",aJets.genEta ,"genEta[naJets]/F");
501     _outTree->Branch("aJet_genPhi",aJets.genPhi ,"genPhi[naJets]/F");
502     _outTree->Branch("aJet_JECUnc",aJets.JECUnc ,"JECUnc[naJets]/F");
503 arizzi 1.24 _outTree->Branch("aJet_vtxMass",aJets.vtxMass ,"vtxMass[naJets]/F");
504     _outTree->Branch("aJet_vtx3dL",aJets.vtx3dL ,"vtx3dL[naJets]/F");
505     _outTree->Branch("aJet_vtx3deL",aJets.vtx3deL ,"vtx3deL[naJets]/F");
506 arizzi 1.36 _outTree->Branch("aJet_id",aJets.id ,"id[naJets]/b");
507 arizzi 1.1
508     _outTree->Branch("numJets" , &numJets , "numJets/I" );
509     _outTree->Branch("numBJets" , &numBJets , "numBJets/I" );
510     _outTree->Branch("deltaPullAngle", &deltaPullAngle , "deltaPullAngle/F");
511     _outTree->Branch("gendrcc" , &gendrcc , "gendrcc/F");
512     _outTree->Branch("gendrbb" , &gendrbb , "gendrbb/F");
513     _outTree->Branch("genZpt" , &genZpt , "genZpt/F");
514     _outTree->Branch("genWpt" , &genWpt , "genWpt/F");
515     _outTree->Branch("weightTrig" , &weightTrig , "weightTrig/F");
516 arizzi 1.38 _outTree->Branch("weightTrigMay" , &weightTrigMay , "weightTrigMay/F");
517     _outTree->Branch("weightTrigV4" , &weightTrigV4 , "weightTrigV4/F");
518 arizzi 1.22 _outTree->Branch("weightTrigMET" , &weightTrigMET , "weightTrigMET/F");
519 arizzi 1.24 _outTree->Branch("weightEleRecoAndId" , &weightEleRecoAndId , "weightEleRecoAndId/F");
520     _outTree->Branch("weightEleTrigJetMETPart" , &weightEleTrigJetMETPart , "weightEleTrigJetMETPart/F");
521     _outTree->Branch("weightEleTrigElePart" , &weightEleTrigElePart , "weightEleTrigElePart/F");
522    
523    
524 arizzi 1.1 _outTree->Branch("deltaPullAngleAK7", &deltaPullAngleAK7 , "deltaPullAngleAK7/F");
525     _outTree->Branch("PUweight", &PUweight , "PUweight/F");
526 arizzi 1.3 _outTree->Branch("eventFlav", &eventFlav , "eventFlav/I");
527 arizzi 1.1
528    
529    
530    
531     _outTree->Branch("Vtype" , &Vtype , "Vtype/I" );
532     _outTree->Branch("HVdPhi" , &HVdPhi , "HVdPhi/F" );
533 arizzi 1.41 _outTree->Branch("HVMass" , &HVMass , "HVMass/F" );
534 arizzi 1.19 _outTree->Branch("HMETdPhi" , &HMETdPhi , "HMETdPhi/F" );
535 arizzi 1.1 _outTree->Branch("VMt" , &VMt , "VMt/F" );
536    
537 arizzi 1.9 _outTree->Branch("nvlep" , &nvlep , "nvlep/I");
538     _outTree->Branch("nalep" , &nalep , "nalep/I");
539 arizzi 1.3
540 arizzi 1.9 _outTree->Branch("vLepton_mass",vLeptons.mass ,"mass[nvlep]/F");
541     _outTree->Branch("vLepton_pt",vLeptons.pt ,"pt[nvlep]/F");
542     _outTree->Branch("vLepton_eta",vLeptons.eta ,"eta[nvlep]");
543     _outTree->Branch("vLepton_phi",vLeptons.phi ,"phi[nvlep]/F");
544     _outTree->Branch("vLepton_aodCombRelIso",vLeptons.aodCombRelIso ,"aodCombRelIso[nvlep]/F");
545     _outTree->Branch("vLepton_pfCombRelIso",vLeptons.pfCombRelIso ,"pfCombRelIso[nvlep]/F");
546     _outTree->Branch("vLepton_photonIso",vLeptons.photonIso ,"photonIso[nvlep]/F");
547     _outTree->Branch("vLepton_neutralHadIso",vLeptons.neutralHadIso ,"neutralHadIso[nvlep]/F");
548     _outTree->Branch("vLepton_chargedHadIso",vLeptons.chargedHadIso ,"chargedHadIso[nvlep]/F");
549     _outTree->Branch("vLepton_particleIso",vLeptons.particleIso ,"particleIso[nvlep]/F");
550     _outTree->Branch("vLepton_dxy",vLeptons.dxy ,"dxy[nvlep]/F");
551     _outTree->Branch("vLepton_dz",vLeptons.dz ,"dz[nvlep]/F");
552     _outTree->Branch("vLepton_type",vLeptons.type ,"type[nvlep]/I");
553 arizzi 1.27 _outTree->Branch("vLepton_id80",vLeptons.id80 ,"id80[nvlep]/F");
554     _outTree->Branch("vLepton_id95",vLeptons.id95 ,"id95[nvlep]/F");
555 arizzi 1.48 _outTree->Branch("vLepton_vbtf",vLeptons.vbtf ,"vbtf[nvlep]/F");
556 arizzi 1.9
557     _outTree->Branch("aLepton_mass",aLeptons.mass ,"mass[nalep]/F");
558     _outTree->Branch("aLepton_pt",aLeptons.pt ,"pt[nalep]/F");
559     _outTree->Branch("aLepton_eta",aLeptons.eta ,"eta[nalep]");
560     _outTree->Branch("aLepton_phi",aLeptons.phi ,"phi[nalep]/F");
561     _outTree->Branch("aLepton_aodCombRelIso",aLeptons.aodCombRelIso ,"aodCombRelIso[nalep]/F");
562     _outTree->Branch("aLepton_pfCombRelIso",aLeptons.pfCombRelIso ,"pfCombRelIso[nalep]/F");
563     _outTree->Branch("aLepton_photonIso",aLeptons.photonIso ,"photonIso[nalep]/F");
564     _outTree->Branch("aLepton_neutralHadIso",aLeptons.neutralHadIso ,"neutralHadIso[nalep]/F");
565     _outTree->Branch("aLepton_chargedHadIso",aLeptons.chargedHadIso ,"chargedHadIso[nalep]/F");
566     _outTree->Branch("aLepton_particleIso",aLeptons.particleIso ,"particleIso[nalep]/F");
567     _outTree->Branch("aLepton_dxy",aLeptons.dxy ,"dxy[nalep]/F");
568     _outTree->Branch("aLepton_dz",aLeptons.dz ,"dz[nalep]/F");
569     _outTree->Branch("aLepton_type",aLeptons.type ,"type[nalep]/I");
570 arizzi 1.48 _outTree->Branch("aLepton_id80",aLeptons.id80 ,"id80[nalep]/F");
571     _outTree->Branch("aLepton_id95",aLeptons.id95 ,"id95[nalep]/F");
572     _outTree->Branch("aLepton_vbtf",aLeptons.vbtf ,"vbtf[nalep]/F");
573 arizzi 1.9
574 arizzi 1.8 _outTree->Branch("top" , &top , "mass/F:pt/F:wMass/F");
575 arizzi 1.49 _outTree->Branch("WplusMode" , &WplusMode , "WplusMode/I");
576     _outTree->Branch("WminusMode" , &WminusMode , "WminusMode/I");
577 arizzi 1.1
578 bortigno 1.30 //IVF
579     _outTree->Branch("nSvs",&nSvs ,"nSvs/I");
580 arizzi 1.33 _outTree->Branch("Sv_massBCand", &IVF.massBcand,"massBcand[nSvs]/F");
581     _outTree->Branch("Sv_massSv", &IVF.massSv,"massSv[nSvs]/F");
582     _outTree->Branch("Sv_pt", &IVF.pt,"pt[nSvs]/F");
583     _outTree->Branch("Sv_eta", &IVF.eta,"eta[nSvs]/F");
584     _outTree->Branch("Sv_phi", &IVF.phi,"phi[nSvs]/F");
585     _outTree->Branch("Sv_dist3D", &IVF.dist3D,"dist3D[nSvs]/F");
586     _outTree->Branch("Sv_dist2D", &IVF.dist2D,"dist2D[nSvs]/F");
587     _outTree->Branch("Sv_distSim2D", &IVF.distSig2D,"distSig2D[nSvs]/F");
588     _outTree->Branch("Sv_distSig3D", &IVF.distSig3D,"distSig3D[nSvs]/F");
589     _outTree->Branch("Sv_dist3D_norm", &IVF.dist3D_norm,"dist3D_norm[nSvs]/F");
590 bortigno 1.30 //IVF higgs candidate
591 arizzi 1.33 _outTree->Branch("SVH" , &SVH , "mass/F:pt/F:eta:phi/F:dR/F:dPhi/F:dEta/F");
592 bortigno 1.30
593    
594     //FIXME : need to update the EDMNtuple with BHadron infos
595     // //SimBHadron
596     // _outTree->Branch("nSimBs",&nSimBs ,"nSimBs/I");
597     // _outTree->Branch("SimBs_mass", &SimBs.mass,"mass[nSvs]/F");
598     // _outTree->Branch("SimBs_pt", &SimBs.pt,"pt[nSvs]/F");
599     // _outTree->Branch("SimBs_eta", &SimBs.eta,"eta[nSvs]/F");
600     // _outTree->Branch("SimBs_phi", &SimBs.phi,"phi[nSvs]/F");
601     // _outTree->Branch("SimBs_vtx_x", &SimBs.vtx_x,"vtx_x[nSvs]/F");
602     // _outTree->Branch("SimBs_vtx_y", &SimBs.vtx_y,"vtx_y[nSvs]/F");
603     // _outTree->Branch("SimBs_vtx_z", &SimBs.vtx_z,"vtx_z[nSvs]/F");
604     // _outTree->Branch("SimBs_pdgId", &SimBs.pdgId,"pdgId[nSvs]/F");
605     // _outTree->Branch("SimBs_status", &SimBs.status,"status[nSvs]/F");
606     // //SimBHadron Higgs Candidate
607     // _outTree->Branch("SimBs_dr", &SimBs_dr,"SimBs_dr/F");
608     // _outTree->Branch("SimBs_dPhi", &SimBs_dPhi,"SimBs_dPhi/F");
609     // _outTree->Branch("SimBs_dEta", &SimBs_dEta,"SimBs_dEta/F");
610     // _outTree->Branch("SimBs_Hmass", &SimBs_Hmass,"SimBs_Hmass/F");
611     // _outTree->Branch("SimBs_Hpt", &SimBs_Hpt,"SimBs_Hpt/F");
612     // _outTree->Branch("SimBs_Heta", &SimBs_Heta,"SimBs_Heta/F");
613     // _outTree->Branch("SimBs_Hphi", &SimBs_Hphi,"SimBs_Hphi/F");
614    
615    
616 arizzi 1.36 _outTree->Branch("rho" , &rho , "rho/F");
617     _outTree->Branch("rho25" , &rho25 , "rho25/F");
618     _outTree->Branch("nPVs" , &nPVs , "nPVs/I");
619     _outTree->Branch("METnoPU" , &METnoPU , "et/F:sumet:sig/F:phi/F");
620 arizzi 1.1 _outTree->Branch("MET" , &MET , "et/F:sumet:sig/F:phi/F");
621     _outTree->Branch("MHT" , &MHT , "mht/F:ht:sig/F:phi/F");
622     _outTree->Branch("minDeltaPhijetMET" , &minDeltaPhijetMET , "minDeltaPhijetMET/F");
623     _outTree->Branch("jetPt_minDeltaPhijetMET" , &jetPt_minDeltaPhijetMET , "jetPt_minDeltaPhijetMET/F");
624    
625 bortigno 1.30 std::stringstream s;
626     s << "triggerFlags[" << triggers.size() << "]/b";
627     _outTree->Branch("triggerFlags", triggerFlags, s.str().c_str());
628 nmohr 1.7
629 arizzi 1.17
630 bortigno 1.30 _outTree->Branch("EVENT" , &EVENT , "run/I:lumi/I:event/I:json/I");
631     _outTree->Branch("hbhe" , &hbhe , "hbhe/b");
632     _outTree->Branch("btag1TSF" , &btag1TSF , "btag1TSF/F");
633     _outTree->Branch("btag2TSF" , &btag2TSF , "btag2TSF/F");
634     _outTree->Branch("btag1T2CSF" , &btag1T2CSF , "btag1T2CSF/F");
635 arizzi 1.44 _outTree->Branch("btag2CSF" , &btag2CSF , "btag2CSF/F");
636     _outTree->Branch("btagA0CSF" , &btagA0CSF , "btagA0CSF/F");
637 arizzi 1.1
638 bortigno 1.30 int ievt=0;
639     int totalcount=0;
640 arizzi 1.1
641 bortigno 1.30 // TFile* inFile = new TFile(inputFile.c_str(), "read");
642 arizzi 1.1 for(unsigned int iFile=0; iFile<inputFiles_.size(); ++iFile) {
643     std::cout << iFile << std::endl;
644     TFile* inFile = TFile::Open(inputFiles_[iFile].c_str());
645     if(inFile==0) continue;
646    
647 bortigno 1.30 // loop the events
648 arizzi 1.1
649     fwlite::Event ev(inFile);
650 nmohr 1.29 for(ev.toBegin(); !ev.atEnd() ; ++ev, ++ievt)
651 arizzi 1.1 {
652 nmohr 1.29 if (ievt <= skipEvents_) continue;
653     if (maxEvents_ >= 0){
654     if (ievt > maxEvents_ + skipEvents_) break;
655     };
656 arizzi 1.16 count->Fill(1.);
657 arizzi 1.1
658 arizzi 1.25 fwlite::Handle< VHbbEventAuxInfo > vhbbAuxHandle;
659     vhbbAuxHandle.getByLabel(ev,"HbbAnalyzerNew");
660     const VHbbEventAuxInfo & aux = *vhbbAuxHandle.product();
661     PUweight=1.;
662 arizzi 1.1 if(isMC_){
663 arizzi 1.41
664 arizzi 1.1 // PU weights
665 bortigno 1.30 std::map<int, unsigned int>::const_iterator puit = aux.puInfo.pus.find(0);
666 arizzi 1.42 PUweight = lumiWeights.weight( puit->second );
667 arizzi 1.41 pu->Fill(puit->second);
668 bortigno 1.30 }
669 arizzi 1.49 countWithPU->Fill(1,PUweight);
670 nmohr 1.7
671 bortigno 1.30 //Write event info
672     EVENT.run = ev.id().run();
673     EVENT.lumi = ev.id().luminosityBlock();
674     EVENT.event = ev.id().event();
675     EVENT.json = jsonContainsEvent (jsonVector, ev);
676    
677     //FIXME : need to update EDM ntuple with BHadron infos
678     // // simBHadrons
679     // fwlite::Handle<SimBHadronCollection> SBHC;
680     // SBHC.getByLabel(ev, "bhadrons");
681     // const SimBHadronCollection sbhc = *(SBHC.product());
682    
683     const std::vector<VHbbCandidate> * candZ ;
684     const std::vector<VHbbCandidate> * candW ;
685     const VHbbEvent * iEvent =0;
686     if(fromCandidate)
687     {
688     fwlite::Handle< std::vector<VHbbCandidate> > vhbbCandHandleZ;
689     vhbbCandHandleZ.getByLabel(ev,"hbbBestCSVPt20Candidates");
690     candZ = vhbbCandHandleZ.product();
691    
692     fwlite::Handle< std::vector<VHbbCandidate> > vhbbCandHandle;
693     vhbbCandHandle.getByLabel(ev,"hbbHighestPtHiggsPt30Candidates");
694     candW = vhbbCandHandle.product();
695     }
696     else
697     {
698     candZlocal->clear();
699     candWlocal->clear();
700     fwlite::Handle< VHbbEvent > vhbbHandle;
701     vhbbHandle.getByLabel(ev,"HbbAnalyzerNew");
702     iEvent = vhbbHandle.product();
703     algoZ->run(vhbbHandle.product(),*candZlocal);
704     algoW->run(vhbbHandle.product(),*candWlocal);
705     candZ= candZlocal;
706     candW= candWlocal;
707    
708     /* for(size_t m=0;m<iEvent->muInfo.size();m++)
709     {
710    
711     if( fabs(iEvent->muInfo[m].p4.Pt()-28.118684) < 0.0001 ||
712     fabs(iEvent->muInfo[m].p4.Pt()-34.853199) < 0.0001 )
713     {
714     std::cout << "FOUND " << iEvent->muInfo[m].p4.Pt() << " " << EVENT.event << " " << candW->size() << " " << candZ->size() << std::endl;
715     }
716     }
717     */
718 arizzi 1.4
719 bortigno 1.30 }
720 arizzi 1.1
721 arizzi 1.3 const std::vector<VHbbCandidate> * cand = candZ;
722 arizzi 1.1
723    
724    
725 bortigno 1.30 /* fwlite::Handle< VHbbEvent > vhbbHandle;
726     vhbbHandle.getByLabel(ev,"HbbAnalyzerNew");
727     const VHbbEvent iEvent = *vhbbHandle.product();
728     */
729 arizzi 1.1
730 bortigno 1.30 // std::clog << "Filling tree "<< std::endl;
731     bool isW=false;
732 arizzi 1.21
733 arizzi 1.17
734 bortigno 1.30 if(cand->size() == 0 or cand->at(0).H.jets.size() < 2) continue;
735     if(cand->size() > 1 )
736 arizzi 1.3 {
737 bortigno 1.30 std::cout << "MULTIPLE CANDIDATES: " << cand->size() << std::endl;
738 arizzi 1.3 }
739 bortigno 1.30 if(cand->at(0).candidateType == VHbbCandidate::Wmun || cand->at(0).candidateType == VHbbCandidate::Wen ) { cand=candW; isW=true; }
740     if(cand->size() == 0)
741 arizzi 1.3 {
742 bortigno 1.30 // std::cout << "W event loss due to tigther cuts" << std::endl;
743 arizzi 1.3 continue;
744     }
745 arizzi 1.35
746     // secondary vtxs
747     fwlite::Handle<std::vector<reco::Vertex> > SVC;
748     SVC.getByLabel(ev,"selectedVertices");
749     const std::vector<reco::Vertex> svc = *(SVC.product());
750    
751 bortigno 1.30 const VHbbCandidate & vhCand = cand->at(0);
752     patFilters.setEvent(&ev,"VH");
753     hbhe = patFilters.accept("hbhe");
754 arizzi 1.18
755 bortigno 1.30 trigger.setEvent(&ev);
756     for(size_t j=0;j < triggers.size();j++)
757 arizzi 1.4 triggerFlags[j]=trigger.accept(triggers[j]);
758 arizzi 1.39
759 bortigno 1.30 eventFlav=0;
760     if(aux.mcBbar.size() > 0 || aux.mcB.size() > 0) eventFlav=5;
761     else if(aux.mcC.size() > 0) eventFlav=4;
762 arizzi 1.8
763 arizzi 1.3
764 bortigno 1.30 H.mass = vhCand.H.p4.M();
765     H.pt = vhCand.H.p4.Pt();
766     H.eta = vhCand.H.p4.Eta();
767     H.phi = vhCand.H.p4.Phi();
768     V.mass = vhCand.V.p4.M();
769 arizzi 1.36 if(isW) V.mass = vhCand.Mt();
770 bortigno 1.30 V.pt = vhCand.V.p4.Pt();
771     V.eta = vhCand.V.p4.Eta();
772     V.phi = vhCand.V.p4.Phi();
773     nhJets=2;
774     hJets.set(vhCand.H.jets[0],0);
775     hJets.set(vhCand.H.jets[1],1);
776     aJets.reset();
777     naJets=vhCand.additionalJets.size();
778     numBJets=0;
779     if(vhCand.H.jets[0].csv> btagThr) numBJets++;
780     if(vhCand.H.jets[1].csv> btagThr) numBJets++;
781     for( int j=0; j < naJets && j < MAXJ; j++ )
782 arizzi 1.6 {
783 bortigno 1.30 aJets.set(vhCand.additionalJets[j],j);
784     if(vhCand.additionalJets[j].csv> btagThr) numBJets++;
785 arizzi 1.6 }
786 bortigno 1.30 numJets = vhCand.additionalJets.size()+2;
787 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());
788     H.dPhi = deltaPhi(vhCand.H.jets[0].p4.Phi(),vhCand.H.jets[1].p4.Phi());
789     H.dEta= TMath::Abs( vhCand.H.jets[0].p4.Eta() - vhCand.H.jets[1].p4.Eta() );
790 bortigno 1.30 HVdPhi = fabs( deltaPhi(vhCand.H.p4.Phi(),vhCand.V.p4.Phi()) ) ;
791 arizzi 1.41 HVMass = (vhCand.H.p4 + vhCand.V.p4).M() ;
792 bortigno 1.30 HMETdPhi = fabs( deltaPhi(vhCand.H.p4.Phi(),vhCand.V.mets.at(0).p4.Phi()) ) ;
793     VMt = vhCand.Mt() ;
794     deltaPullAngle = vhCand.H.deltaTheta;
795    
796     hJets.cosTheta[0]= vhCand.H.helicities[0];
797     hJets.cosTheta[1]= vhCand.H.helicities[1];
798 arizzi 1.36 // METInfo calomet; METInfo tcmet; METInfo pfmet; METInfo mht; METInfo metNoPU
799 arizzi 1.24
800 bortigno 1.30 MET.et = vhCand.V.mets.at(0).p4.Pt();
801     MET.phi = vhCand.V.mets.at(0).p4.Phi();
802     MET.sumet = vhCand.V.mets.at(0).sumEt;
803     MET.sig = vhCand.V.mets.at(0).metSig;
804 arizzi 1.24
805 arizzi 1.36 METnoPU.et = iEvent->metNoPU.p4.Pt();
806     METnoPU.phi = iEvent->metNoPU.p4.Phi();
807     METnoPU.sumet = iEvent->metNoPU.sumEt;
808     METnoPU.sig = iEvent->metNoPU.metSig;
809    
810     rho = aux.puInfo.rho;
811     rho25 = aux.puInfo.rho25;
812     nPVs=aux.pvInfo.nVertices;
813    
814 arizzi 1.24 if(!fromCandidate) {
815     MHT.mht = iEvent->mht.p4.Pt();
816     MHT.phi = iEvent->mht.p4.Phi();
817     MHT.ht = iEvent->mht.sumEt;
818     MHT.sig = iEvent->mht.metSig;
819 bortigno 1.30 }
820    
821     Vtype = vhCand.candidateType;
822 arizzi 1.24
823    
824 bortigno 1.30 //Secondary Vertices
825     IVF.reset();
826     nSvs = svc.size();
827     const TVector3 recoPv = aux.pvInfo.firstPVInPT2;
828     const math::XYZPoint myPv(recoPv);
829     //look here for Matrix filling info http://project-mathlibs.web.cern.ch/project-mathlibs/sw/html/SMatrixDoc.html
830     std::vector<double> fillMatrix(6);
831     for (int i = 0; i<6; ++i) fillMatrix[i] = 0.;
832     fillMatrix[0] = TMath::Power(0.002,2);
833     fillMatrix[2] = TMath::Power(0.002,2);
834     fillMatrix[5] = TMath::Power(0.002,2);
835     const ROOT::Math::SMatrix<double, 3, 3, ROOT::Math::MatRepSym<double, 3> > myFakeMatrixError(fillMatrix.begin(),fillMatrix.end());
836     const reco::Vertex recoVtxPv(myPv, myFakeMatrixError);
837     for( int j=0; j < nSvs && j < MAXB; ++j ) {
838     const GlobalVector flightDir = flightDirection(recoPv,svc[j]);
839     reco::SecondaryVertex recoSv(recoVtxPv, svc[j], flightDir ,true);
840     IVF.set( recoSv, recoPv ,j);
841     }
842     if(nSvs > 2){
843     TLorentzVector BCands_H1, BCands_H2, BCands_H;
844     BCands_H1.SetPtEtaPhiM(IVF.pt[0], IVF.eta[0], IVF.phi[0], IVF.massBcand[0]);
845     BCands_H2.SetPtEtaPhiM(IVF.pt[1], IVF.eta[1], IVF.phi[1], IVF.massBcand[1]);
846     BCands_H = BCands_H1 + BCands_H2;
847 arizzi 1.33 SVH.dR = deltaR(IVF.eta[0], IVF.phi[0], IVF.eta[1], IVF.phi[1] );
848     SVH.dPhi = deltaPhi(IVF.phi[0], IVF.phi[1] );
849     SVH.dEta = TMath::Abs(IVF.eta[0] - IVF.eta[1] );
850     SVH.mass = BCands_H.M();
851     SVH.pt = BCands_H.Pt();
852     SVH.eta = BCands_H.Eta();
853     SVH.phi = BCands_H.Phi();
854 bortigno 1.30 }
855    
856     //FIXME : need to update EDM ntuple with simBhadron info
857     // //SimBhadron
858     // SimBs.reset();
859     // nSimBs = sbhc.size();
860     // for( int j=0; j < nSimBs && j < MAXB; ++j )
861     // SimBs.set( sbhc.at(j), j);
862     // if(nSimBs > 2){
863     // TLorentzVector SimBs_H1, SimBs_H2, SimBs_H;
864     // SimBs_H1.SetPtEtaPhiM(SimBs.pt[0], SimBs.eta[0], SimBs.phi[0], SimBs.mass[0]);
865     // SimBs_H2.SetPtEtaPhiM(SimBs.pt[1], SimBs.eta[1], SimBs.phi[1], SimBs.mass[1]);
866     // SimBs_H = SimBs_H1 + SimBs_H2;
867     // SimBs_dr = deltaR(SimBs.eta[0], SimBs.phi[0], SimBs.eta[1], SimBs.phi[1] );
868     // SimBs_dPhi = deltaPhi(SimBs.phi[0], SimBs.phi[1] );
869     // SimBs_dEta = TMath::Abs(SimBs.eta[0] - SimBs.eta[1] );
870     // SimBs_Hmass = SimBs_H.M();
871     // SimBs_Hpt = SimBs_H.Pt();
872     // SimBs_Heta = SimBs_H.Eta();
873     // SimBs_Hphi = SimBs_H.Phi();
874     // }
875    
876    
877     //Loop on jets
878     double maxBtag=-99999;
879     minDeltaPhijetMET = 999;
880     TLorentzVector bJet;
881     std::vector<std::vector<BTagWeight::JetInfo> > btagJetInfos;
882     std::vector<float> jet30eta;
883     std::vector<float> jet30pt;
884     if(fromCandidate)
885     {
886     //Loop on Higgs Jets
887     for(unsigned int j=0; j < vhCand.H.jets.size(); j++ ){
888     if (vhCand.H.jets[j].csv > maxBtag) { bJet=vhCand.H.jets[j].p4 ; maxBtag =vhCand.H.jets[j].csv; }
889 arizzi 1.39 if (fabs(deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.H.jets[j].p4.Phi())) < minDeltaPhijetMET)
890 bortigno 1.30 {
891 arizzi 1.39 minDeltaPhijetMET=fabs(deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.H.jets[j].p4.Phi()));
892 arizzi 1.24 jetPt_minDeltaPhijetMET=vhCand.H.jets[j].p4.Pt();
893 bortigno 1.30 }
894     btagJetInfos.push_back(btagEff.jetInfo(vhCand.H.jets[j]));
895     }
896     //Loop on Additional Jets
897     for(unsigned int j=0; j < vhCand.additionalJets.size(); j++ ){
898     if (vhCand.additionalJets[j].csv > maxBtag) { bJet=vhCand.additionalJets[j].p4 ; maxBtag =vhCand.additionalJets[j].csv; }
899 arizzi 1.45 /* if (fabs(deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.additionalJets[j].p4.Phi())) < minDeltaPhijetMET)
900 bortigno 1.30 {
901 arizzi 1.39 minDeltaPhijetMET=fabs(deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.additionalJets[j].p4.Phi()));
902 arizzi 1.24 jetPt_minDeltaPhijetMET=vhCand.additionalJets[j].p4.Pt();
903 arizzi 1.45 }*/
904 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
905     {
906     if(vhCand.additionalJets[j].p4.Pt() > 20)
907     btagJetInfos.push_back(btagEff.jetInfo(vhCand.additionalJets[j]));
908     }
909     }
910     }
911     else
912     {
913 arizzi 1.45 //Loop on all jets
914 bortigno 1.30 for(unsigned int j=0; j < iEvent->simpleJets2.size(); j++ ){
915     if (iEvent->simpleJets2[j].csv > maxBtag) { bJet=iEvent->simpleJets2[j].p4 ; maxBtag =iEvent->simpleJets2[j].csv; }
916 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)
917 arizzi 1.46 {
918     minDeltaPhijetMET=fabs(deltaPhi( vhCand.V.mets.at(0).p4.Phi(), iEvent->simpleJets2[j].p4.Phi()));
919     jetPt_minDeltaPhijetMET=iEvent->simpleJets2[j].p4.Pt();
920     }
921 bortigno 1.30 if(iEvent->simpleJets2[j].p4.Pt() > 30)
922     {
923     jet30eta.push_back(iEvent->simpleJets2[j].p4.Eta());
924     jet30pt.push_back(iEvent->simpleJets2[j].p4.Pt());
925     }
926 arizzi 1.24
927 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)
928     // 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
929     if( ( isW && ! useHighestPtHiggsW ) || ( ! isW && ! useHighestPtHiggsZ ) )
930     {
931     if(iEvent->simpleJets2[j].p4.Pt() > 20)
932     btagJetInfos.push_back(btagEff.jetInfo(iEvent->simpleJets2[j]));
933     }
934     }
935    
936 arizzi 1.45 //Loop on Higgs jets
937    
938     for(unsigned int j=0; j < vhCand.H.jets.size(); j++ ) {
939    
940 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
941     // by a criteria (pt of the dijet) that is not btag SF dependent
942     if(!( ( isW && ! useHighestPtHiggsW ) || ( ! isW && ! useHighestPtHiggsZ ) )) {
943 arizzi 1.45 btagJetInfos.push_back(btagEff.jetInfo(vhCand.H.jets[j]));
944 bortigno 1.30 }
945 arizzi 1.45 }
946 arizzi 1.24
947 bortigno 1.30 }
948     vLeptons.reset();
949     weightTrig = 1.; // better to default to 1
950 arizzi 1.38 weightTrigMay = -1.;
951     weightTrigV4 = -1.;
952 bortigno 1.30 TLorentzVector leptonForTop;
953     size_t firstAddMu=0;
954     size_t firstAddEle=0;
955     if(Vtype == VHbbCandidate::Zmumu ){
956     vLeptons.set(vhCand.V.muons[0],0,13);
957     vLeptons.set(vhCand.V.muons[1],1,13);
958     float cweightID = triggerWeight.scaleMuID(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleMuID(vLeptons.pt[1],vLeptons.eta[1]) ;
959     float weightTrig1 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[0],vLeptons.eta[0]);
960     float weightTrig2 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[1],vLeptons.eta[1]);
961     float cweightTrig = weightTrig1 + weightTrig2 - weightTrig1*weightTrig2;
962     weightTrig = cweightID * cweightTrig;
963     nvlep=2;
964     firstAddMu=2;
965     }
966     if( Vtype == VHbbCandidate::Zee ){
967     vLeptons.set(vhCand.V.electrons[0],0,11);
968     vLeptons.set(vhCand.V.electrons[1],1,11);
969     nvlep=2;
970     firstAddEle=2;
971     std::vector<float> pt,eta;
972     pt.push_back(vLeptons.pt[0]); eta.push_back(vLeptons.eta[0]);
973     pt.push_back(vLeptons.pt[1]); eta.push_back(vLeptons.eta[1]);
974     weightEleRecoAndId=triggerWeight.scaleID95Ele(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleRecoEle(vLeptons.pt[0],vLeptons.eta[0]) *
975     triggerWeight.scaleID95Ele(vLeptons.pt[1],vLeptons.eta[1]) * triggerWeight.scaleRecoEle(vLeptons.pt[1],vLeptons.eta[1]);
976     weightEleTrigElePart = triggerWeight.scaleDoubleEle17Ele8(pt,eta);
977     weightTrig = weightEleTrigElePart * weightEleRecoAndId;
978 arizzi 1.24
979    
980    
981 bortigno 1.30 }
982     if(Vtype == VHbbCandidate::Wmun ){
983     leptonForTop=vhCand.V.muons[0].p4;
984     vLeptons.set(vhCand.V.muons[0],0,13);
985     float cweightID = triggerWeight.scaleMuID(vLeptons.pt[0],vLeptons.eta[0]);
986     float weightTrig1 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[0],vLeptons.eta[0]);
987     float cweightTrig = weightTrig1;
988     weightTrig = cweightID * cweightTrig;
989     nvlep=1;
990     firstAddMu=1;
991     }
992     if( Vtype == VHbbCandidate::Wen ){
993     leptonForTop=vhCand.V.electrons[0].p4;
994     vLeptons.set(vhCand.V.electrons[0],0,11);
995     nvlep=1;
996     firstAddEle=1;
997 arizzi 1.38 weightTrigMay = triggerWeight.scaleSingleEleMay(vLeptons.pt[0],vLeptons.eta[0]);
998     weightTrigV4 = triggerWeight.scaleSingleEleV4(vLeptons.pt[0],vLeptons.eta[0]);
999 bortigno 1.30 weightEleRecoAndId=triggerWeight.scaleID80Ele(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleRecoEle(vLeptons.pt[0],vLeptons.eta[0]);
1000     weightEleTrigJetMETPart=triggerWeight.scaleJet30Jet25(jet30eta,jet30pt)*triggerWeight.scalePFMHTEle(MET.et);
1001 arizzi 1.38 weightEleTrigElePart= weightTrigV4; //this is for debugging only, checking only the V4 part
1002 bortigno 1.30
1003 arizzi 1.38 weightTrigMay*=weightEleRecoAndId;
1004     weightTrigV4*=weightEleRecoAndId;
1005     weightTrigV4*=weightEleTrigJetMETPart;
1006     weightTrig = weightTrigMay * 0.187 + weightTrigV4 * (1.-0.187); //FIXME: use proper lumi if we reload 2.fb
1007 bortigno 1.30 }
1008     if( Vtype == VHbbCandidate::Znn ){
1009     nvlep=0;
1010     float weightTrig1 = triggerWeight.scaleMetHLT(vhCand.V.mets.at(0).p4.Pt());
1011     weightTrig = weightTrig1;
1012     weightTrigMET = weightTrig1;
1013     }
1014 arizzi 1.38
1015     if(weightTrigMay < 0) weightTrigMay=weightTrig;
1016     if(weightTrigV4 < 0) weightTrigV4=weightTrig;
1017 arizzi 1.40 if(!isMC_)
1018     {
1019     weightTrig = 1.;
1020     weightTrigMay = 1.;
1021     weightTrigV4 = 1.;
1022    
1023     }
1024 bortigno 1.30 aLeptons.reset();
1025     nalep=0;
1026     if(fromCandidate)
1027 arizzi 1.21 {
1028     for(size_t j=firstAddMu;j< vhCand.V.muons.size();j++) aLeptons.set(vhCand.V.muons[j],nalep++,13);
1029     for(size_t j=firstAddEle;j< vhCand.V.electrons.size();j++) aLeptons.set(vhCand.V.electrons[j],nalep++,11);
1030     }
1031 bortigno 1.30 else
1032 arizzi 1.21 {
1033 bortigno 1.30 for(size_t j=0;j< iEvent->muInfo.size();j++)
1034     {
1035 arizzi 1.21 if((j!= vhCand.V.firstLepton && j!= vhCand.V.secondLepton) || ((Vtype != VHbbCandidate::Wmun ) && (Vtype != VHbbCandidate::Zmumu )) )
1036 bortigno 1.30 aLeptons.set(iEvent->muInfo[j],nalep++,13);
1037     }
1038     for(size_t j=0;j< iEvent->eleInfo.size();j++)
1039     {
1040 arizzi 1.21 if((j!= vhCand.V.firstLepton && j!= vhCand.V.secondLepton) || ((Vtype != VHbbCandidate::Wen ) && (Vtype != VHbbCandidate::Zee )))
1041 bortigno 1.30 aLeptons.set(iEvent->eleInfo[j],nalep++,11);
1042     }
1043 arizzi 1.8
1044 arizzi 1.21 }
1045 arizzi 1.8
1046 arizzi 1.24
1047 bortigno 1.30 if(isMC_)
1048     {
1049     //std::cout << "BTAGSF " << btagJetInfos.size() << " " << btag.weight<BTag1Tight2CustomFilter>(btagJetInfos) << std::endl;
1050     if ( btagJetInfos.size()< 10)
1051     {
1052     btag1T2CSF = btag.weight<BTag1Tight2CustomFilter>(btagJetInfos);
1053     btag2TSF = btag.weight<BTag2TightFilter>(btagJetInfos);
1054     btag1TSF = btag.weight<BTag1TightFilter>(btagJetInfos);
1055 arizzi 1.44 btagA0CSF = btag.weight<BTagAntiMax0CustomFilter>(btagJetInfos);
1056     btag2CSF = btag.weight<BTag2CustomFilter>(btagJetInfos);
1057 bortigno 1.30 }
1058     else
1059     {
1060     std::cout << "WARNING: combinatorics for " << btagJetInfos.size() << " jets is too high (>=10). use SF=1 " << std::endl;
1061     //TODO: revert to random throw for this cases
1062     btag1T2CSF = 1.;
1063     btag2TSF = 1.;
1064     btag1TSF = 1.;
1065     }
1066     }
1067 arizzi 1.17
1068 bortigno 1.30 if(maxBtag > -99999)
1069 arizzi 1.8 {
1070 bortigno 1.30 TopHypo topQuark = TopMassReco::topMass(leptonForTop,bJet,vhCand.V.mets.at(0).p4);
1071     top.mass = topQuark.p4.M();
1072     top.pt = topQuark.p4.Pt();
1073     top.wMass = topQuark.p4W.M();
1074 arizzi 1.8 } else {
1075 bortigno 1.30 top.mass = -99;
1076     top.pt = -99;
1077     top.wMass = -99;
1078     }
1079 arizzi 1.8
1080    
1081    
1082 bortigno 1.30 //FIXME: too much warnings... figure out why
1083     // gendrcc=aux.genCCDeltaR();
1084     // gendrbb=aux.genBBDeltaR();
1085     genZpt=aux.mcZ.size() > 0 ? aux.mcZ[0].p4.Pt():-99;
1086     genWpt=aux.mcW.size() > 0 ? aux.mcW[0].p4.Pt():-99;
1087 arizzi 1.49 WminusMode=-99;
1088     WplusMode=-99;
1089     for(unsigned int j=0; j<= aux.mcW.size();j++)
1090     {
1091     for(unsigned int k=0;k< aux.mcW[j].dauid.size();k++)
1092     {
1093     int idd=abs(aux.mcW[j].dauid[k]);
1094     if(idd==11 || idd==13 || idd==15|| (idd<=5 && idd >=1))
1095     {
1096     if(WminusMode==-99 && aux.mcW[j].charge ==-1) WminusMode = idd;
1097     if(WplusMode==-99 && aux.mcW[j].charge ==-1) WplusMode = idd;
1098     }
1099     }
1100 arizzi 1.24
1101 arizzi 1.49
1102     }
1103 arizzi 1.25
1104 arizzi 1.32 /// Compute pull angle from AK7
1105 arizzi 1.25 if(!fromCandidate){
1106 arizzi 1.32 std::vector<VHbbEvent::SimpleJet> ak7wrt1(iEvent->simpleJets3);
1107     std::vector<VHbbEvent::SimpleJet> ak7wrt2(iEvent->simpleJets3);
1108     if(ak7wrt1.size() > 1){
1109 arizzi 1.24 CompareDeltaR deltaRComparatorJ1(vhCand.H.jets[0].p4);
1110     CompareDeltaR deltaRComparatorJ2(vhCand.H.jets[1].p4);
1111 arizzi 1.32 std::sort( ak7wrt1.begin(),ak7wrt1.end(),deltaRComparatorJ1 );
1112     std::sort( ak7wrt2.begin(),ak7wrt2.end(),deltaRComparatorJ2 );
1113 arizzi 1.24 std::vector<VHbbEvent::SimpleJet> ak7_matched;
1114 arizzi 1.32 // if the matched are different save them
1115     if(ak7wrt1[0].p4.DeltaR(ak7wrt2[0].p4) > 0.1) {
1116     ak7_matched.push_back(ak7wrt1[0]);
1117     ak7_matched.push_back(ak7wrt2[0]);
1118     }
1119     // else look at the second best
1120     else{
1121     // ak7wrt1 is best
1122     if( ak7wrt1[1].p4.DeltaR(vhCand.H.jets[0].p4) < ak7wrt2[1].p4.DeltaR(vhCand.H.jets[1].p4))
1123     {
1124     ak7_matched.push_back(ak7wrt1[1]);
1125     ak7_matched.push_back(ak7wrt2[0]);
1126     }
1127     else
1128     {
1129     ak7_matched.push_back(ak7wrt1[0]);
1130     ak7_matched.push_back(ak7wrt2[1]);
1131     }
1132     }
1133 arizzi 1.24 CompareJetPt ptComparator;
1134     std::sort( ak7_matched.begin(),ak7_matched.end(),ptComparator );
1135     if(ak7_matched[0].p4.DeltaR(vhCand.H.jets[0].p4) < 0.5
1136 arizzi 1.32 and ak7_matched[1].p4.DeltaR(vhCand.H.jets[1].p4) < 0.5)
1137     {
1138     deltaPullAngleAK7 = VHbbCandidateTools::getDeltaTheta(ak7_matched[0],ak7_matched[1]);
1139     }
1140 arizzi 1.24 }
1141     }
1142 arizzi 1.3
1143 bortigno 1.30 _outTree->Fill();
1144 arizzi 1.4
1145 bortigno 1.30 }// closed event loop
1146 arizzi 1.1
1147 bortigno 1.30 std::cout << "closing the file: " << inputFiles_[iFile] << std::endl;
1148     inFile->Close();
1149     // close input file
1150     } // loop on files
1151 arizzi 1.1
1152    
1153 bortigno 1.30 std::cout << "Events: " << ievt <<std::endl;
1154     std::cout << "TotalCount: " << totalcount <<std::endl;
1155 arizzi 1.1
1156    
1157    
1158 bortigno 1.30 _outFile->cd();
1159 arizzi 1.1
1160 bortigno 1.30 _outTree->Write();
1161     _outFile->Write();
1162     _outFile->Close();
1163     return 0;
1164 arizzi 1.1 }
1165    
1166