ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbbAnalysis/VHbbDataFormats/bin/Ntupler.cc
Revision: 1.38
Committed: Mon Oct 3 17:13:58 2011 UTC (13 years, 7 months ago) by arizzi
Content type: text/plain
Branch: MAIN
CVS Tags: AR_Oct5Ntuple
Changes since 1.37: +16 -9 lines
Log Message: 
as proposed by david, add two additional weightTrig one for May10 and one for Prompt V4 and later periods to allow later combination with proper relative weights. For channels where it doesnt matter I copy the weightTrig value in both V4 and May10 ones

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