ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbbAnalysis/VHbbDataFormats/bin/Ntupler.cc
Revision: 1.29
Committed: Mon Sep 26 08:13:21 2011 UTC (13 years, 7 months ago) by nmohr
Content type: text/plain
Branch: MAIN
Changes since 1.28: +6 -1 lines
Log Message:
Add maxEvents and skipEvents

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 arizzi 1.16
33     //Move class definition to Ntupler.h ?
34     //#include "VHbbAnalysis/VHbbDataFormats/interface/Ntupler.h"
35    
36     #include "VHbbAnalysis/VHbbDataFormats/interface/BTagWeight.h"
37     #include "VHbbAnalysis/VHbbDataFormats/interface/TriggerWeight.h"
38    
39 arizzi 1.1 #include <sstream>
40     #include <string>
41    
42     #define MAXJ 30
43     #define MAXL 10
44    
45 arizzi 1.24 struct CompareDeltaR {
46     CompareDeltaR(TLorentzVector p4dir_): p4dir(p4dir_) {}
47     bool operator()( const VHbbEvent::SimpleJet& j1, const VHbbEvent::SimpleJet& j2 ) const {
48     return j1.p4.DeltaR(p4dir) > j2.p4.DeltaR(p4dir);
49     }
50     TLorentzVector p4dir;
51     };
52 arizzi 1.4
53 nmohr 1.7 bool jsonContainsEvent (const std::vector< edm::LuminosityBlockRange > &jsonVec,
54     const edm::EventBase &event)
55     {
56     // if the jsonVec is empty, then no JSON file was provided so all
57     // events should pass
58     if (jsonVec.empty())
59     {
60     return true;
61     }
62     bool (* funcPtr) (edm::LuminosityBlockRange const &,
63     edm::LuminosityBlockID const &) = &edm::contains;
64     edm::LuminosityBlockID lumiID (event.id().run(),
65     event.id().luminosityBlock());
66     std::vector< edm::LuminosityBlockRange >::const_iterator iter =
67     std::find_if (jsonVec.begin(), jsonVec.end(),
68     boost::bind(funcPtr, _1, lumiID) );
69     return jsonVec.end() != iter;
70    
71     }
72    
73 arizzi 1.1
74    
75     typedef struct
76     {
77     float mass; //MT in case of W
78     float pt;
79     float eta;
80     float phi;
81 arizzi 1.9 } TrackInfo;
82 arizzi 1.1
83    
84 arizzi 1.9 struct LeptonInfo
85 arizzi 1.1 {
86     void reset()
87     {
88     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; }
89     }
90    
91     template <class Input> void set(const Input & i, int j,int t)
92     {
93     type[j]=t;
94     pt[j]=i.p4.Pt();
95     mass[j]=i.p4.M();
96     eta[j]=i.p4.Eta();
97     phi[j]=i.p4.Phi();
98     aodCombRelIso[j]=(i.hIso+i.eIso+i.tIso)/i.p4.Pt();
99     pfCombRelIso[j]=(i.pfChaIso+i.pfPhoIso+i.pfNeuIso)/i.p4.Pt();
100     photonIso[j]=i.pfPhoIso;
101     neutralHadIso[j]=i.pfNeuIso;
102     chargedHadIso[j]=i.pfChaIso;
103 arizzi 1.27 setSpecific(i,j);
104 arizzi 1.1 }
105 arizzi 1.27 template <class Input> void setSpecific(const Input & i, int j)
106 arizzi 1.4 {
107     }
108    
109    
110    
111 arizzi 1.1
112     float mass[MAXL]; //MT in case of W
113     float pt[MAXL];
114     float eta[MAXL];
115     float phi[MAXL];
116     float aodCombRelIso[MAXL];
117     float pfCombRelIso[MAXL];
118     float photonIso[MAXL];
119     float neutralHadIso[MAXL];
120     float chargedHadIso[MAXL];
121     float particleIso[MAXL];
122     float dxy[MAXL];
123     float dz[MAXL];
124     int type[MAXL];
125 arizzi 1.27 float id80[MAXL];
126     float id95[MAXL];
127 arizzi 1.1 };
128    
129 arizzi 1.27 template <> void LeptonInfo::setSpecific<VHbbEvent::ElectronInfo>(const VHbbEvent::ElectronInfo & i, int j){
130     id80[j]=i.id80r;
131     id95[j]=i.id95r;
132     }
133     template <> void LeptonInfo::setSpecific<VHbbEvent::MuonInfo>(const VHbbEvent::MuonInfo & i, int j){
134     dxy[j]=i.ipDb;
135     dz[j]=i.zPVPt;
136 arizzi 1.4 }
137    
138 arizzi 1.1 typedef struct
139     {
140     float et;
141     float sumet;
142     float sig;
143     float phi;
144 arizzi 1.9 } METInfo;
145 arizzi 1.1
146     typedef struct
147     {
148     float mht;
149     float ht;
150     float sig;
151     float phi;
152 arizzi 1.9 } MHTInfo;
153 nmohr 1.5
154     typedef struct
155     {
156     float mass;
157     float pt;
158     float wMass;
159 arizzi 1.9 } TopInfo;
160 arizzi 1.1
161 nmohr 1.7 typedef struct
162 arizzi 1.1 {
163     int run;
164     int lumi;
165     int event;
166 nmohr 1.7 int json;
167 arizzi 1.9 } EventInfo;
168 arizzi 1.1
169     typedef struct
170     {
171     void set(const VHbbEvent::SimpleJet & j, int i)
172     {
173     pt[i]=j.p4.Pt();
174     eta[i]=j.p4.Eta();
175     phi[i]=j.p4.Phi();
176     csv[i]=j.csv;
177 arizzi 1.24 numTracksSV[i] = j.vtxMass;
178     vtxMass[i]= j.vtxMass;
179     vtx3dL[i] = j.vtx3dL;
180     vtx3deL[i] = j.vtx3deL;
181     chf[i]=j.chargedHadronEFraction;
182     nhf[i] =j.neutralHadronEFraction;
183     cef[i] =j.chargedEmEFraction;
184     nef[i] =j.neutralEmEFraction;
185     nconstituents[i] =j.nConstituents;
186     nch[i]=j.ntracks;
187 arizzi 1.3
188 arizzi 1.1 flavour[i]=j.flavour;
189 arizzi 1.3 if(j.bestMCp4.Pt() > 0)
190     {
191     genPt[i]=j.bestMCp4.Pt();
192     genEta[i]=j.bestMCp4.Eta();
193     genPhi[i]=j.bestMCp4.Phi();
194     }
195 arizzi 1.24 JECUnc[i]=j.jecunc;
196 arizzi 1.3
197 arizzi 1.1 }
198     void reset()
199     {
200     for(int i=0;i<MAXJ;i++) {
201     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;
202     }
203     }
204     float pt[MAXJ];
205     float eta[MAXJ];
206     float phi[MAXJ];
207     float csv[MAXJ];
208     float cosTheta[MAXJ];
209     int numTracksSV[MAXJ];
210     float chf[MAXJ];
211     float nhf[MAXJ];
212     float cef[MAXJ];
213     float nef[MAXJ];
214     float nch[MAXJ];
215     float nconstituents[MAXJ];
216     float flavour[MAXJ];
217     float genPt[MAXJ];
218     float genEta[MAXJ];
219     float genPhi[MAXJ];
220     float JECUnc[MAXJ];
221 arizzi 1.24 float vtxMass[MAXJ];
222     float vtx3dL [MAXJ];
223     float vtx3deL[MAXJ];
224 arizzi 1.9 } JetInfo;
225 arizzi 1.1
226     int main(int argc, char* argv[])
227     {
228     gROOT->Reset();
229    
230     TTree *_outTree;
231 arizzi 1.9 METInfo MET;
232     MHTInfo MHT;
233     TopInfo top;
234     EventInfo EVENT;
235     // JetInfo jet1,jet2, addJet1, addJet2;
236     // lepton1,lepton2;
237     JetInfo hJets, aJets;
238     LeptonInfo vLeptons, aLeptons;
239 arizzi 1.1 int naJets=0, nhJets=0;
240 arizzi 1.9 TrackInfo H;
241     TrackInfo V;
242     int nvlep=0,nalep=0;
243 arizzi 1.1
244 arizzi 1.22 float jjdr,jjdPhi,jjdEta,HVdPhi,HMETdPhi,VMt,deltaPullAngle,deltaPullAngleAK7,gendrcc,gendrbb, genZpt, genWpt, weightTrig, weightTrigMET, minDeltaPhijetMET, jetPt_minDeltaPhijetMET , PUweight;
245 arizzi 1.24 float weightEleRecoAndId,weightEleTrigJetMETPart, weightEleTrigElePart;
246    
247     int Vtype,numJets,numBJets,eventFlav;
248 arizzi 1.1 // bool isMET80_CJ80, ispfMHT150, isMET80_2CJ20,isMET65_2CJ20, isJETID,isIsoMu17;
249 arizzi 1.18 bool triggerFlags[500],hbhe;
250 arizzi 1.17 float btag1T2CSF=1.,btag2TSF=1.,btag1TSF=1.;
251 arizzi 1.1 // ----------------------------------------------------------------------
252     // First Part:
253     //
254     // * enable the AutoLibraryLoader
255     // * book the histograms of interest
256     // * open the input file
257     // ----------------------------------------------------------------------
258    
259     // load framework libraries
260     gSystem->Load( "libFWCoreFWLite" );
261     gSystem->Load("libDataFormatsFWLite");
262     AutoLibraryLoader::enable();
263    
264     // parse arguments
265     if ( argc < 2 ) {
266     return 0;
267     }
268    
269 arizzi 1.4 std::vector<VHbbCandidate> * candZlocal = new std::vector<VHbbCandidate>;
270     std::vector<VHbbCandidate> * candWlocal = new std::vector<VHbbCandidate>;
271    
272 arizzi 1.1 // get the python configuration
273     PythonProcessDesc builder(argv[1]);
274     const edm::ParameterSet& in = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("fwliteInput" );
275     const edm::ParameterSet& out = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("fwliteOutput");
276     const edm::ParameterSet& ana = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("Analyzer");
277 nmohr 1.7 std::vector<edm::LuminosityBlockRange> jsonVector;
278     if ( in.exists("lumisToProcess") )
279     {
280     std::vector<edm::LuminosityBlockRange> const & lumisTemp =
281     in.getUntrackedParameter<std::vector<edm::LuminosityBlockRange> > ("lumisToProcess");
282     jsonVector.resize( lumisTemp.size() );
283     copy( lumisTemp.begin(), lumisTemp.end(), jsonVector.begin() );
284     }
285 arizzi 1.1
286     // now get each parameter
287     int maxEvents_( in.getParameter<int>("maxEvents") );
288 nmohr 1.29 int skipEvents_( in.getParameter<int>("skipEvents") );
289 arizzi 1.1 unsigned int outputEvery_( in.getParameter<unsigned int>("outputEvery") );
290     std::string outputFile_( out.getParameter<std::string>("fileName" ) );
291     std::vector<std::string> triggers( ana.getParameter<std::vector<std::string> >("triggers") );
292 arizzi 1.6 double btagThr = ana.getParameter<double>("bJetCountThreshold" );
293 arizzi 1.4 bool fromCandidate = ana.getParameter<bool>("readFromCandidates");
294 arizzi 1.17 bool useHighestPtHiggsZ = ana.getParameter<bool>("useHighestPtHiggsZ");
295     bool useHighestPtHiggsW = ana.getParameter<bool>("useHighestPtHiggsW");
296     HbbCandidateFinderAlgo * algoZ = new HbbCandidateFinderAlgo(ana.getParameter<bool>("verbose"), ana.getParameter<double>("jetPtThresholdZ"),useHighestPtHiggsZ );
297     HbbCandidateFinderAlgo * algoW = new HbbCandidateFinderAlgo(ana.getParameter<bool>("verbose"), ana.getParameter<double>("jetPtThresholdW"),useHighestPtHiggsW);
298 arizzi 1.4
299 arizzi 1.16 TriggerWeight triggerWeight(ana);
300 arizzi 1.17 BTagWeight btag(2); // 2 operating points "Custom" = 0.5 and "Tight = 0.898"
301     BTagSampleEfficiency btagEff( ana.getParameter<std::string>("btagEffFileName" ).c_str() );
302 arizzi 1.1
303     std::vector<std::string> inputFiles_( in.getParameter<std::vector<std::string> >("fileNames") );
304     // std::string inputFile( in.getParameter<std::string> ("fileName") );
305    
306    
307 arizzi 1.3 std::string PUmcfileName_ = in.getParameter<std::string> ("PUmcfileName") ;
308     std::string PUdatafileName_ = in.getParameter<std::string> ("PUdatafileName") ;
309 arizzi 1.25
310 arizzi 1.1 bool isMC_( ana.getParameter<bool>("isMC") );
311     TriggerReader trigger(isMC_);
312 arizzi 1.18 TriggerReader patFilters(false);
313 arizzi 1.25
314     edm::LumiReWeighting lumiWeights;
315     if(isMC_)
316     {
317 arizzi 1.26 lumiWeights = edm::LumiReWeighting(PUmcfileName_,PUdatafileName_ , "pileup", "pileup");
318 arizzi 1.25 }
319 arizzi 1.1
320 arizzi 1.6 // TFile *_outPUFile = new TFile((outputFile_+"_PU").c_str(), "recreate");
321     // TH1F * pu = new TH1F("pileup","",-0.5,24.5,25);
322 arizzi 1.1 TFile *_outFile = new TFile(outputFile_.c_str(), "recreate");
323 arizzi 1.16 TH1F * count = new TH1F("Count","Count", 1,0,2 );
324     TH1F * countWithPU = new TH1F("CountWithPU","CountWithPU", 1,0,2 );
325 arizzi 1.1 _outTree = new TTree("tree", "myTree");
326    
327     _outTree->Branch("H" , &H , "mass/F:pt/F:eta:phi/F");
328     _outTree->Branch("V" , &V , "mass/F:pt/F:eta:phi/F");
329     _outTree->Branch("nhJets" , &nhJets , "nhJets/I");
330     _outTree->Branch("naJets" , &naJets , "naJets/I");
331 arizzi 1.2
332     _outTree->Branch("hJet_pt",hJets.pt ,"pt[nhJets]/F");
333     _outTree->Branch("hJet_eta",hJets.eta ,"eta[nhJets]/F");
334     _outTree->Branch("hJet_phi",hJets.phi ,"phi[nhJets]/F");
335     _outTree->Branch("hJet_csv",hJets.csv ,"csv[nhJets]/F");
336     _outTree->Branch("hJet_cosTheta",hJets.cosTheta ,"cosTheta[nhJets]/F");
337     _outTree->Branch("hJet_numTracksSV",hJets.numTracksSV ,"numTracksSV[nhJets]/I");
338     _outTree->Branch("hJet_chf",hJets.chf ,"chf[nhJets]/F");
339     _outTree->Branch("hJet_nhf",hJets.nhf ,"nhf[nhJets]/F");
340     _outTree->Branch("hJet_cef",hJets.cef ,"cef[nhJets]/F");
341     _outTree->Branch("hJet_nef",hJets.nef ,"nef[nhJets]/F");
342     _outTree->Branch("hJet_nch",hJets.nch ,"nch[nhJets]/F");
343     _outTree->Branch("hJet_nconstituents",hJets.nconstituents ,"nconstituents[nhJets]");
344     _outTree->Branch("hJet_flavour",hJets.flavour ,"flavour[nhJets]/F");
345     _outTree->Branch("hJet_genPt",hJets.genPt ,"genPt[nhJets]/F");
346     _outTree->Branch("hJet_genEta",hJets.genEta ,"genEta[nhJets]/F");
347     _outTree->Branch("hJet_genPhi",hJets.genPhi ,"genPhi[nhJets]/F");
348     _outTree->Branch("hJet_JECUnc",hJets.JECUnc ,"JECUnc[nhJets]/F");
349 arizzi 1.24 _outTree->Branch("hJet_vtxMass",hJets.vtxMass ,"vtxMass[naJets]/F");
350     _outTree->Branch("hJet_vtx3dL",hJets.vtx3dL ,"vtx3dL[naJets]/F");
351     _outTree->Branch("hJet_vtx3deL",hJets.vtx3deL ,"vtx3deL[naJets]/F");
352 arizzi 1.2
353     _outTree->Branch("aJet_pt",aJets.pt ,"pt[naJets]/F");
354     _outTree->Branch("aJet_eta",aJets.eta ,"eta[naJets]/F");
355     _outTree->Branch("aJet_phi",aJets.phi ,"phi[naJets]/F");
356     _outTree->Branch("aJet_csv",aJets.csv ,"csv[naJets]/F");
357     _outTree->Branch("aJet_cosTheta",aJets.cosTheta ,"cosTheta[naJets]/F");
358     _outTree->Branch("aJet_numTracksSV",aJets.numTracksSV ,"numTracksSV[naJets]/I");
359     _outTree->Branch("aJet_chf",aJets.chf ,"chf[naJets]/F");
360     _outTree->Branch("aJet_nhf",aJets.nhf ,"nhf[naJets]/F");
361     _outTree->Branch("aJet_cef",aJets.cef ,"cef[naJets]/F");
362     _outTree->Branch("aJet_nef",aJets.nef ,"nef[naJets]/F");
363     _outTree->Branch("aJet_nch",aJets.nch ,"nch[naJets]/F");
364     _outTree->Branch("aJet_nconstituents",aJets.nconstituents ,"nconstituents[naJets]");
365     _outTree->Branch("aJet_flavour",aJets.flavour ,"flavour[naJets]/F");
366     _outTree->Branch("aJet_genPt",aJets.genPt ,"genPt[naJets]/F");
367     _outTree->Branch("aJet_genEta",aJets.genEta ,"genEta[naJets]/F");
368     _outTree->Branch("aJet_genPhi",aJets.genPhi ,"genPhi[naJets]/F");
369     _outTree->Branch("aJet_JECUnc",aJets.JECUnc ,"JECUnc[naJets]/F");
370 arizzi 1.24 _outTree->Branch("aJet_vtxMass",aJets.vtxMass ,"vtxMass[naJets]/F");
371     _outTree->Branch("aJet_vtx3dL",aJets.vtx3dL ,"vtx3dL[naJets]/F");
372     _outTree->Branch("aJet_vtx3deL",aJets.vtx3deL ,"vtx3deL[naJets]/F");
373 arizzi 1.1
374     _outTree->Branch("jjdr" , &jjdr , "jjdr/F" );
375     _outTree->Branch("jjdPhi" , &jjdPhi , "jjdPhi/F" );
376 arizzi 1.6 _outTree->Branch("jjdEta" , &jjdEta , "jjdEta/F" );
377 arizzi 1.1 _outTree->Branch("numJets" , &numJets , "numJets/I" );
378     _outTree->Branch("numBJets" , &numBJets , "numBJets/I" );
379     _outTree->Branch("deltaPullAngle", &deltaPullAngle , "deltaPullAngle/F");
380     _outTree->Branch("gendrcc" , &gendrcc , "gendrcc/F");
381     _outTree->Branch("gendrbb" , &gendrbb , "gendrbb/F");
382     _outTree->Branch("genZpt" , &genZpt , "genZpt/F");
383     _outTree->Branch("genWpt" , &genWpt , "genWpt/F");
384     _outTree->Branch("weightTrig" , &weightTrig , "weightTrig/F");
385 arizzi 1.22 _outTree->Branch("weightTrigMET" , &weightTrigMET , "weightTrigMET/F");
386 arizzi 1.24 _outTree->Branch("weightEleRecoAndId" , &weightEleRecoAndId , "weightEleRecoAndId/F");
387     _outTree->Branch("weightEleTrigJetMETPart" , &weightEleTrigJetMETPart , "weightEleTrigJetMETPart/F");
388     _outTree->Branch("weightEleTrigElePart" , &weightEleTrigElePart , "weightEleTrigElePart/F");
389    
390    
391 arizzi 1.1 _outTree->Branch("deltaPullAngleAK7", &deltaPullAngleAK7 , "deltaPullAngleAK7/F");
392     _outTree->Branch("PUweight", &PUweight , "PUweight/F");
393 arizzi 1.3 _outTree->Branch("eventFlav", &eventFlav , "eventFlav/I");
394 arizzi 1.1
395    
396    
397    
398     _outTree->Branch("Vtype" , &Vtype , "Vtype/I" );
399     _outTree->Branch("HVdPhi" , &HVdPhi , "HVdPhi/F" );
400 arizzi 1.19 _outTree->Branch("HMETdPhi" , &HMETdPhi , "HMETdPhi/F" );
401 arizzi 1.1 _outTree->Branch("VMt" , &VMt , "VMt/F" );
402    
403 arizzi 1.9 _outTree->Branch("nvlep" , &nvlep , "nvlep/I");
404     _outTree->Branch("nalep" , &nalep , "nalep/I");
405 arizzi 1.3
406 arizzi 1.9 _outTree->Branch("vLepton_mass",vLeptons.mass ,"mass[nvlep]/F");
407     _outTree->Branch("vLepton_pt",vLeptons.pt ,"pt[nvlep]/F");
408     _outTree->Branch("vLepton_eta",vLeptons.eta ,"eta[nvlep]");
409     _outTree->Branch("vLepton_phi",vLeptons.phi ,"phi[nvlep]/F");
410     _outTree->Branch("vLepton_aodCombRelIso",vLeptons.aodCombRelIso ,"aodCombRelIso[nvlep]/F");
411     _outTree->Branch("vLepton_pfCombRelIso",vLeptons.pfCombRelIso ,"pfCombRelIso[nvlep]/F");
412     _outTree->Branch("vLepton_photonIso",vLeptons.photonIso ,"photonIso[nvlep]/F");
413     _outTree->Branch("vLepton_neutralHadIso",vLeptons.neutralHadIso ,"neutralHadIso[nvlep]/F");
414     _outTree->Branch("vLepton_chargedHadIso",vLeptons.chargedHadIso ,"chargedHadIso[nvlep]/F");
415     _outTree->Branch("vLepton_particleIso",vLeptons.particleIso ,"particleIso[nvlep]/F");
416     _outTree->Branch("vLepton_dxy",vLeptons.dxy ,"dxy[nvlep]/F");
417     _outTree->Branch("vLepton_dz",vLeptons.dz ,"dz[nvlep]/F");
418     _outTree->Branch("vLepton_type",vLeptons.type ,"type[nvlep]/I");
419 arizzi 1.27 _outTree->Branch("vLepton_id80",vLeptons.id80 ,"id80[nvlep]/F");
420     _outTree->Branch("vLepton_id95",vLeptons.id95 ,"id95[nvlep]/F");
421 arizzi 1.9
422     _outTree->Branch("aLepton_mass",aLeptons.mass ,"mass[nalep]/F");
423     _outTree->Branch("aLepton_pt",aLeptons.pt ,"pt[nalep]/F");
424     _outTree->Branch("aLepton_eta",aLeptons.eta ,"eta[nalep]");
425     _outTree->Branch("aLepton_phi",aLeptons.phi ,"phi[nalep]/F");
426     _outTree->Branch("aLepton_aodCombRelIso",aLeptons.aodCombRelIso ,"aodCombRelIso[nalep]/F");
427     _outTree->Branch("aLepton_pfCombRelIso",aLeptons.pfCombRelIso ,"pfCombRelIso[nalep]/F");
428     _outTree->Branch("aLepton_photonIso",aLeptons.photonIso ,"photonIso[nalep]/F");
429     _outTree->Branch("aLepton_neutralHadIso",aLeptons.neutralHadIso ,"neutralHadIso[nalep]/F");
430     _outTree->Branch("aLepton_chargedHadIso",aLeptons.chargedHadIso ,"chargedHadIso[nalep]/F");
431     _outTree->Branch("aLepton_particleIso",aLeptons.particleIso ,"particleIso[nalep]/F");
432     _outTree->Branch("aLepton_dxy",aLeptons.dxy ,"dxy[nalep]/F");
433     _outTree->Branch("aLepton_dz",aLeptons.dz ,"dz[nalep]/F");
434     _outTree->Branch("aLepton_type",aLeptons.type ,"type[nalep]/I");
435 arizzi 1.27 _outTree->Branch("aLepton_id80",aLeptons.id80 ,"id[nalep]/F");
436     _outTree->Branch("aLepton_id95",aLeptons.id95 ,"id[nalep]/F");
437 arizzi 1.9
438 arizzi 1.8 _outTree->Branch("top" , &top , "mass/F:pt/F:wMass/F");
439 arizzi 1.1
440     _outTree->Branch("MET" , &MET , "et/F:sumet:sig/F:phi/F");
441     _outTree->Branch("MHT" , &MHT , "mht/F:ht:sig/F:phi/F");
442     _outTree->Branch("minDeltaPhijetMET" , &minDeltaPhijetMET , "minDeltaPhijetMET/F");
443     _outTree->Branch("jetPt_minDeltaPhijetMET" , &jetPt_minDeltaPhijetMET , "jetPt_minDeltaPhijetMET/F");
444    
445     std::stringstream s;
446     s << "triggerFlags[" << triggers.size() << "]/b";
447     _outTree->Branch("triggerFlags", triggerFlags, s.str().c_str());
448 nmohr 1.7
449 arizzi 1.17
450 nmohr 1.7 _outTree->Branch("EVENT" , &EVENT , "run/I:lumi/I:event/I:json/I");
451 arizzi 1.18 _outTree->Branch("hbhe" , &hbhe , "hbhe/b");
452 arizzi 1.17 _outTree->Branch("btag1TSF" , &btag1TSF , "btag1TSF/F");
453     _outTree->Branch("btag2TSF" , &btag2TSF , "btag2TSF/F");
454     _outTree->Branch("btag1T2CSF" , &btag1T2CSF , "btag1T2CSF/F");
455 arizzi 1.1
456     int ievt=0;
457     int totalcount=0;
458    
459     // TFile* inFile = new TFile(inputFile.c_str(), "read");
460     for(unsigned int iFile=0; iFile<inputFiles_.size(); ++iFile) {
461     std::cout << iFile << std::endl;
462     TFile* inFile = TFile::Open(inputFiles_[iFile].c_str());
463     if(inFile==0) continue;
464    
465     // loop the events
466    
467     fwlite::Event ev(inFile);
468 nmohr 1.29 for(ev.toBegin(); !ev.atEnd() ; ++ev, ++ievt)
469 arizzi 1.1 {
470 nmohr 1.29 if (ievt <= skipEvents_) continue;
471     if (maxEvents_ >= 0){
472     if (ievt > maxEvents_ + skipEvents_) break;
473     };
474 arizzi 1.16 count->Fill(1.);
475 arizzi 1.1
476 arizzi 1.25 fwlite::Handle< VHbbEventAuxInfo > vhbbAuxHandle;
477     vhbbAuxHandle.getByLabel(ev,"HbbAnalyzerNew");
478     const VHbbEventAuxInfo & aux = *vhbbAuxHandle.product();
479     PUweight=1.;
480 arizzi 1.1 if(isMC_){
481     // PU weights
482 arizzi 1.25 std::map<int, unsigned int>::const_iterator puit = aux.puInfo.pus.find(0);
483 tboccali 1.28 PUweight = lumiWeights.weight( puit->second /3. );
484 arizzi 1.1 }
485 arizzi 1.16 countWithPU->Fill(PUweight);
486 nmohr 1.7
487     //Write event info
488     EVENT.run = ev.id().run();
489     EVENT.lumi = ev.id().luminosityBlock();
490     EVENT.event = ev.id().event();
491     EVENT.json = jsonContainsEvent (jsonVector, ev);
492 arizzi 1.3
493 arizzi 1.4
494     const std::vector<VHbbCandidate> * candZ ;
495     const std::vector<VHbbCandidate> * candW ;
496 arizzi 1.13 const VHbbEvent * iEvent =0;
497 arizzi 1.4 if(fromCandidate)
498     {
499 arizzi 1.3 fwlite::Handle< std::vector<VHbbCandidate> > vhbbCandHandleZ;
500     vhbbCandHandleZ.getByLabel(ev,"hbbBestCSVPt20Candidates");
501 arizzi 1.4 candZ = vhbbCandHandleZ.product();
502 arizzi 1.3
503     fwlite::Handle< std::vector<VHbbCandidate> > vhbbCandHandle;
504     vhbbCandHandle.getByLabel(ev,"hbbHighestPtHiggsPt30Candidates");
505 arizzi 1.4 candW = vhbbCandHandle.product();
506     }
507     else
508     {
509     candZlocal->clear();
510     candWlocal->clear();
511     fwlite::Handle< VHbbEvent > vhbbHandle;
512     vhbbHandle.getByLabel(ev,"HbbAnalyzerNew");
513 arizzi 1.13 iEvent = vhbbHandle.product();
514 arizzi 1.4 algoZ->run(vhbbHandle.product(),*candZlocal);
515     algoW->run(vhbbHandle.product(),*candWlocal);
516     candZ= candZlocal;
517     candW= candWlocal;
518    
519 arizzi 1.27 /* for(size_t m=0;m<iEvent->muInfo.size();m++)
520     {
521 arizzi 1.4
522 arizzi 1.27 if( fabs(iEvent->muInfo[m].p4.Pt()-28.118684) < 0.0001 ||
523     fabs(iEvent->muInfo[m].p4.Pt()-34.853199) < 0.0001 )
524     {
525     std::cout << "FOUND " << iEvent->muInfo[m].p4.Pt() << " " << EVENT.event << " " << candW->size() << " " << candZ->size() << std::endl;
526     }
527     }
528     */
529 arizzi 1.1
530 arizzi 1.4 }
531 arizzi 1.1
532 arizzi 1.3 const std::vector<VHbbCandidate> * cand = candZ;
533 arizzi 1.1
534    
535    
536     /* fwlite::Handle< VHbbEvent > vhbbHandle;
537     vhbbHandle.getByLabel(ev,"HbbAnalyzerNew");
538     const VHbbEvent iEvent = *vhbbHandle.product();
539     */
540    
541     // std::clog << "Filling tree "<< std::endl;
542 arizzi 1.17 bool isW=false;
543 arizzi 1.21
544 arizzi 1.17
545 arizzi 1.3 if(cand->size() == 0 or cand->at(0).H.jets.size() < 2) continue;
546     if(cand->size() > 1 )
547     {
548     std::cout << "MULTIPLE CANDIDATES: " << cand->size() << std::endl;
549     }
550 arizzi 1.17 if(cand->at(0).candidateType == VHbbCandidate::Wmun || cand->at(0).candidateType == VHbbCandidate::Wen ) { cand=candW; isW=true; }
551 arizzi 1.3 if(cand->size() == 0)
552     {
553     // std::cout << "W event loss due to tigther cuts" << std::endl;
554     continue;
555     }
556     const VHbbCandidate & vhCand = cand->at(0);
557 arizzi 1.18 patFilters.setEvent(&ev,"VH");
558     hbhe = patFilters.accept("hbhe");
559    
560 arizzi 1.4 trigger.setEvent(&ev);
561     for(size_t j=0;j < triggers.size();j++)
562     triggerFlags[j]=trigger.accept(triggers[j]);
563 arizzi 1.3
564     eventFlav=0;
565     if(aux.mcBbar.size() > 0 || aux.mcB.size() > 0) eventFlav=5;
566     else if(aux.mcC.size() > 0) eventFlav=4;
567 arizzi 1.8
568 arizzi 1.3
569 arizzi 1.1 H.mass = vhCand.H.p4.M();
570     H.pt = vhCand.H.p4.Pt();
571     H.eta = vhCand.H.p4.Eta();
572     H.phi = vhCand.H.p4.Phi();
573     V.mass = vhCand.V.p4.M();
574     V.pt = vhCand.V.p4.Pt();
575     V.eta = vhCand.V.p4.Eta();
576     V.phi = vhCand.V.p4.Phi();
577     nhJets=2;
578     hJets.set(vhCand.H.jets[0],0);
579     hJets.set(vhCand.H.jets[1],1);
580     aJets.reset();
581     naJets=vhCand.additionalJets.size();
582 arizzi 1.6 numBJets=0;
583     if(vhCand.H.jets[0].csv> btagThr) numBJets++;
584     if(vhCand.H.jets[1].csv> btagThr) numBJets++;
585     for( int j=0; j < naJets && j < MAXJ; j++ )
586     {
587     aJets.set(vhCand.additionalJets[j],j);
588     if(vhCand.additionalJets[j].csv> btagThr) numBJets++;
589     }
590 arizzi 1.1 numJets = vhCand.additionalJets.size()+2;
591     jjdr = 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());
592     jjdPhi = deltaPhi(vhCand.H.jets[0].p4.Phi(),vhCand.H.jets[1].p4.Phi());
593 arizzi 1.6 jjdEta= TMath::Abs( vhCand.H.jets[0].p4.Eta() - vhCand.H.jets[1].p4.Eta() );
594 arizzi 1.19 HVdPhi = fabs( deltaPhi(vhCand.H.p4.Phi(),vhCand.V.p4.Phi()) ) ;
595 arizzi 1.20 HMETdPhi = fabs( deltaPhi(vhCand.H.p4.Phi(),vhCand.V.mets.at(0).p4.Phi()) ) ;
596 nmohr 1.14 VMt = vhCand.Mt() ;
597 arizzi 1.1 deltaPullAngle = vhCand.H.deltaTheta;
598 arizzi 1.16
599 arizzi 1.1
600 arizzi 1.3 hJets.cosTheta[0]= vhCand.H.helicities[0];
601     hJets.cosTheta[1]= vhCand.H.helicities[1];
602 arizzi 1.1
603     MET.et = vhCand.V.mets.at(0).p4.Pt();
604     MET.phi = vhCand.V.mets.at(0).p4.Phi();
605     MET.sumet = vhCand.V.mets.at(0).sumEt;
606     MET.sig = vhCand.V.mets.at(0).metSig;
607 arizzi 1.24
608    
609     if(!fromCandidate) {
610     MHT.mht = iEvent->mht.p4.Pt();
611     MHT.phi = iEvent->mht.p4.Phi();
612     MHT.ht = iEvent->mht.sumEt;
613     MHT.sig = iEvent->mht.metSig;
614    
615     }
616    
617    
618 arizzi 1.1 Vtype = vhCand.candidateType;
619 arizzi 1.24
620     //Loop on jets
621     double maxBtag=-99999;
622     minDeltaPhijetMET = 999;
623     TLorentzVector bJet;
624     std::vector<std::vector<BTagWeight::JetInfo> > btagJetInfos;
625     std::vector<float> jet30eta;
626     std::vector<float> jet30pt;
627     if(fromCandidate)
628     {
629     //Loop on Higgs Jets
630     for(unsigned int j=0; j < vhCand.H.jets.size(); j++ ){
631     if (vhCand.H.jets[j].csv > maxBtag) { bJet=vhCand.H.jets[j].p4 ; maxBtag =vhCand.H.jets[j].csv; }
632     if (deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.H.jets[j].p4.Phi()) < minDeltaPhijetMET)
633     {
634     minDeltaPhijetMET=deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.H.jets[j].p4.Phi());
635     jetPt_minDeltaPhijetMET=vhCand.H.jets[j].p4.Pt();
636     }
637     btagJetInfos.push_back(btagEff.jetInfo(vhCand.H.jets[j]));
638     }
639     //Loop on Additional Jets
640     for(unsigned int j=0; j < vhCand.additionalJets.size(); j++ ){
641     if (vhCand.additionalJets[j].csv > maxBtag) { bJet=vhCand.additionalJets[j].p4 ; maxBtag =vhCand.additionalJets[j].csv; }
642     if (deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.additionalJets[j].p4.Phi()) < minDeltaPhijetMET)
643     {
644     minDeltaPhijetMET=deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.additionalJets[j].p4.Phi());
645     jetPt_minDeltaPhijetMET=vhCand.additionalJets[j].p4.Pt();
646     }
647     if( ( isW && ! useHighestPtHiggsW ) || ( ! isW && ! useHighestPtHiggsZ ) ) // btag SF computed using only H-jets if best-H made with dijetPt rather than best CSV
648     {
649     if(vhCand.additionalJets[j].p4.Pt() > 20)
650     btagJetInfos.push_back(btagEff.jetInfo(vhCand.additionalJets[j]));
651     }
652     }
653     }
654     else
655     {
656     for(unsigned int j=0; j < iEvent->simpleJets2.size(); j++ ){
657     if (iEvent->simpleJets2[j].csv > maxBtag) { bJet=iEvent->simpleJets2[j].p4 ; maxBtag =iEvent->simpleJets2[j].csv; }
658     if (deltaPhi( vhCand.V.mets.at(0).p4.Phi(), iEvent->simpleJets2[j].p4.Phi()) < minDeltaPhijetMET)
659     {
660     minDeltaPhijetMET=deltaPhi( vhCand.V.mets.at(0).p4.Phi(), iEvent->simpleJets2[j].p4.Phi());
661     jetPt_minDeltaPhijetMET=iEvent->simpleJets2[j].p4.Pt();
662     }
663     if(iEvent->simpleJets2[j].p4.Pt() > 30)
664     {
665     jet30eta.push_back(iEvent->simpleJets2[j].p4.Eta());
666     jet30pt.push_back(iEvent->simpleJets2[j].p4.Pt());
667     }
668    
669     //For events made with highest CSV, all jets in the event should be taken into account for "tagging" SF (anti tagging is a mess)
670     // 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
671     if( ( isW && ! useHighestPtHiggsW ) || ( ! isW && ! useHighestPtHiggsZ ) )
672     {
673     if(iEvent->simpleJets2[j].p4.Pt() > 20)
674     btagJetInfos.push_back(btagEff.jetInfo(iEvent->simpleJets2[j]));
675     }
676     }
677    
678     //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
679     // by a criteria (pt of the dijet) that is not btag SF dependent
680     if(!( ( isW && ! useHighestPtHiggsW ) || ( ! isW && ! useHighestPtHiggsZ ) )) {
681     for(unsigned int j=0; j < vhCand.H.jets.size(); j++ ) btagJetInfos.push_back(btagEff.jetInfo(vhCand.H.jets[j]));
682     }
683    
684     }
685 arizzi 1.9 vLeptons.reset();
686 arizzi 1.15 weightTrig = 1.; // better to default to 1
687 arizzi 1.8 TLorentzVector leptonForTop;
688 arizzi 1.9 size_t firstAddMu=0;
689     size_t firstAddEle=0;
690 arizzi 1.1 if(Vtype == VHbbCandidate::Zmumu ){
691 arizzi 1.13 vLeptons.set(vhCand.V.muons[0],0,13);
692     vLeptons.set(vhCand.V.muons[1],1,13);
693 arizzi 1.16 float cweightID = triggerWeight.scaleMuID(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleMuID(vLeptons.pt[1],vLeptons.eta[1]) ;
694     float weightTrig1 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[0],vLeptons.eta[0]);
695     float weightTrig2 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[1],vLeptons.eta[1]);
696 arizzi 1.1 float cweightTrig = weightTrig1 + weightTrig2 - weightTrig1*weightTrig2;
697     weightTrig = cweightID * cweightTrig;
698 arizzi 1.9 nvlep=2;
699     firstAddMu=2;
700 arizzi 1.1 }
701     if( Vtype == VHbbCandidate::Zee ){
702 arizzi 1.9 vLeptons.set(vhCand.V.electrons[0],0,11);
703     vLeptons.set(vhCand.V.electrons[1],1,11);
704     nvlep=2;
705     firstAddEle=2;
706 arizzi 1.24 std::vector<float> pt,eta;
707     pt.push_back(vLeptons.pt[0]); eta.push_back(vLeptons.eta[0]);
708     pt.push_back(vLeptons.pt[1]); eta.push_back(vLeptons.eta[1]);
709     weightEleRecoAndId=triggerWeight.scaleID95Ele(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleRecoEle(vLeptons.pt[0],vLeptons.eta[0]) *
710     triggerWeight.scaleID95Ele(vLeptons.pt[1],vLeptons.eta[1]) * triggerWeight.scaleRecoEle(vLeptons.pt[1],vLeptons.eta[1]);
711     weightEleTrigElePart = triggerWeight.scaleDoubleEle17Ele8(pt,eta);
712     weightTrig = weightEleTrigElePart * weightEleRecoAndId;
713    
714    
715    
716 arizzi 1.1 }
717     if(Vtype == VHbbCandidate::Wmun ){
718 arizzi 1.8 leptonForTop=vhCand.V.muons[0].p4;
719 arizzi 1.13 vLeptons.set(vhCand.V.muons[0],0,13);
720 arizzi 1.16 float cweightID = triggerWeight.scaleMuID(vLeptons.pt[0],vLeptons.eta[0]);
721     float weightTrig1 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[0],vLeptons.eta[0]);
722 arizzi 1.1 float cweightTrig = weightTrig1;
723     weightTrig = cweightID * cweightTrig;
724 arizzi 1.9 nvlep=1;
725     firstAddMu=1;
726 arizzi 1.1 }
727     if( Vtype == VHbbCandidate::Wen ){
728 arizzi 1.8 leptonForTop=vhCand.V.electrons[0].p4;
729 arizzi 1.9 vLeptons.set(vhCand.V.electrons[0],0,11);
730     nvlep=1;
731     firstAddEle=1;
732 arizzi 1.24 float weightMay = triggerWeight.scaleSingleEleMay(vLeptons.pt[0],vLeptons.eta[0]);
733     float weightV4 = triggerWeight.scaleSingleEleV4(vLeptons.pt[0],vLeptons.eta[0]);
734     weightEleRecoAndId=triggerWeight.scaleID80Ele(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleRecoEle(vLeptons.pt[0],vLeptons.eta[0]);
735     weightEleTrigJetMETPart=triggerWeight.scaleJet30Jet25(jet30eta,jet30pt)*triggerWeight.scalePFMHTEle(MET.et);
736     weightEleTrigElePart= weightV4; //this is for debugging only, checking only the V4 part
737    
738     weightMay*=weightEleRecoAndId;
739     weightV4*=weightEleRecoAndId;
740     weightV4*=weightEleTrigJetMETPart;
741     weightTrig = weightMay * 0.187 + weightV4 * (1.-0.187); //FIXME: use proper lumi if we reload 2.fb
742 arizzi 1.1 }
743 arizzi 1.3 if( Vtype == VHbbCandidate::Znn ){
744 arizzi 1.9 nvlep=0;
745 arizzi 1.22 float weightTrig1 = triggerWeight.scaleMetHLT(vhCand.V.mets.at(0).p4.Pt());
746     weightTrig = weightTrig1;
747     weightTrigMET = weightTrig1;
748 arizzi 1.3 }
749 arizzi 1.9
750     aLeptons.reset();
751     nalep=0;
752 arizzi 1.21 if(fromCandidate)
753     {
754     for(size_t j=firstAddMu;j< vhCand.V.muons.size();j++) aLeptons.set(vhCand.V.muons[j],nalep++,13);
755     for(size_t j=firstAddEle;j< vhCand.V.electrons.size();j++) aLeptons.set(vhCand.V.electrons[j],nalep++,11);
756     }
757     else
758     {
759     for(size_t j=0;j< iEvent->muInfo.size();j++)
760     {
761     if((j!= vhCand.V.firstLepton && j!= vhCand.V.secondLepton) || ((Vtype != VHbbCandidate::Wmun ) && (Vtype != VHbbCandidate::Zmumu )) )
762     aLeptons.set(iEvent->muInfo[j],nalep++,13);
763     }
764     for(size_t j=0;j< iEvent->eleInfo.size();j++)
765     {
766     if((j!= vhCand.V.firstLepton && j!= vhCand.V.secondLepton) || ((Vtype != VHbbCandidate::Wen ) && (Vtype != VHbbCandidate::Zee )))
767     aLeptons.set(iEvent->eleInfo[j],nalep++,11);
768     }
769 arizzi 1.8
770 arizzi 1.21 }
771 arizzi 1.8
772 arizzi 1.24
773 arizzi 1.17 if(isMC_)
774     {
775     //std::cout << "BTAGSF " << btagJetInfos.size() << " " << btag.weight<BTag1Tight2CustomFilter>(btagJetInfos) << std::endl;
776 arizzi 1.24 if ( btagJetInfos.size()< 10)
777     {
778 arizzi 1.17 btag1T2CSF = btag.weight<BTag1Tight2CustomFilter>(btagJetInfos);
779 arizzi 1.23 btag2TSF = btag.weight<BTag2TightFilter>(btagJetInfos);
780 arizzi 1.24 btag1TSF = btag.weight<BTag1TightFilter>(btagJetInfos);
781     }
782     else
783     {
784     std::cout << "WARNING: combinatorics for " << btagJetInfos.size() << " jets is too high (>=10). use SF=1 " << std::endl;
785 arizzi 1.25 //TODO: revert to random throw for this cases
786 arizzi 1.24 btag1T2CSF = 1.;
787     btag2TSF = 1.;
788     btag1TSF = 1.;
789     }
790 arizzi 1.16 }
791 arizzi 1.17
792 arizzi 1.8 if(maxBtag > -99999)
793     {
794     TopHypo topQuark = TopMassReco::topMass(leptonForTop,bJet,vhCand.V.mets.at(0).p4);
795     top.mass = topQuark.p4.M();
796     top.pt = topQuark.p4.Pt();
797     top.wMass = topQuark.p4W.M();
798     } else {
799     top.mass = -99;
800     top.pt = -99;
801     top.wMass = -99;
802     }
803    
804    
805    
806 arizzi 1.26 //FIXME: too much warnings... figure out why
807     // gendrcc=aux.genCCDeltaR();
808     // gendrbb=aux.genBBDeltaR();
809 arizzi 1.24 genZpt=aux.mcZ.size() > 0 ? aux.mcZ[0].p4.Pt():-99;
810     genWpt=aux.mcW.size() > 0 ? aux.mcW[0].p4.Pt():-99;
811    
812 arizzi 1.25
813     /// Compute pull angle from AK7
814     if(!fromCandidate){
815     std::vector<VHbbEvent::SimpleJet> ak7 = iEvent->simpleJets3;
816 arizzi 1.24 if(ak7.size() > 1){
817     CompareDeltaR deltaRComparatorJ1(vhCand.H.jets[0].p4);
818     CompareDeltaR deltaRComparatorJ2(vhCand.H.jets[1].p4);
819     std::vector<VHbbEvent::SimpleJet> ak7_matched;
820    
821     std::sort( ak7.begin(),ak7.end(),deltaRComparatorJ1 );
822     ak7_matched.push_back(ak7[0]);
823    
824     if(ak7[1].p4.DeltaR(vhCand.H.jets[0].p4) > 0.5)
825     ak7.erase(ak7.begin());
826    
827     std::sort( ak7.begin(),ak7.end(),deltaRComparatorJ2 );
828 arizzi 1.25 if( abs(ak7[0].p4.Pt() - ak7_matched[1].p4.Pt()) > 0.1 )
829     ak7_matched.push_back(ak7[0]);
830     else
831     ak7_matched.push_back(ak7[1]);
832 arizzi 1.24
833     CompareJetPt ptComparator;
834     std::sort( ak7_matched.begin(),ak7_matched.end(),ptComparator );
835     if(ak7_matched[0].p4.DeltaR(vhCand.H.jets[0].p4) < 0.5
836     and ak7_matched[1].p4.DeltaR(vhCand.H.jets[1].p4) < 0.5 ){
837     deltaPullAngleAK7 = VHbbCandidateTools::getDeltaTheta(ak7_matched[0],ak7_matched[1]);
838     }
839     }
840     }
841 arizzi 1.3
842 arizzi 1.25
843 arizzi 1.1
844    
845 arizzi 1.17
846 arizzi 1.1 _outTree->Fill();
847 arizzi 1.4
848 arizzi 1.1 }// closed event loop
849    
850     std::cout << "closing the file: " << inputFiles_[iFile] << std::endl;
851     inFile->Close();
852     // close input file
853     } // loop on files
854    
855    
856     std::cout << "Events: " << ievt <<std::endl;
857     std::cout << "TotalCount: " << totalcount <<std::endl;
858    
859    
860    
861     _outFile->cd();
862    
863     _outTree->Write();
864     _outFile->Write();
865     _outFile->Close();
866     return 0;
867     }
868    
869