ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbbAnalysis/VHbbDataFormats/bin/Ntupler.cc
Revision: 1.23
Committed: Wed Sep 14 15:12:56 2011 UTC (13 years, 8 months ago) by arizzi
Content type: text/plain
Branch: MAIN
Changes since 1.22: +2 -2 lines
Log Message: 
fix in BTag SF.. not use the same 3 times

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.4
46 nmohr 1.7 bool jsonContainsEvent (const std::vector< edm::LuminosityBlockRange > &jsonVec,
47     const edm::EventBase &event)
48     {
49     // if the jsonVec is empty, then no JSON file was provided so all
50     // events should pass
51     if (jsonVec.empty())
52     {
53     return true;
54     }
55     bool (* funcPtr) (edm::LuminosityBlockRange const &,
56     edm::LuminosityBlockID const &) = &edm::contains;
57     edm::LuminosityBlockID lumiID (event.id().run(),
58     event.id().luminosityBlock());
59     std::vector< edm::LuminosityBlockRange >::const_iterator iter =
60     std::find_if (jsonVec.begin(), jsonVec.end(),
61     boost::bind(funcPtr, _1, lumiID) );
62     return jsonVec.end() != iter;
63    
64     }
65    
66 arizzi 1.1
67    
68     typedef struct
69     {
70     float mass; //MT in case of W
71     float pt;
72     float eta;
73     float phi;
74 arizzi 1.9 } TrackInfo;
75 arizzi 1.1
76    
77 arizzi 1.9 struct LeptonInfo
78 arizzi 1.1 {
79     void reset()
80     {
81     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; }
82     }
83    
84     template <class Input> void set(const Input & i, int j,int t)
85     {
86     type[j]=t;
87     pt[j]=i.p4.Pt();
88     mass[j]=i.p4.M();
89     eta[j]=i.p4.Eta();
90     phi[j]=i.p4.Phi();
91     aodCombRelIso[j]=(i.hIso+i.eIso+i.tIso)/i.p4.Pt();
92     pfCombRelIso[j]=(i.pfChaIso+i.pfPhoIso+i.pfNeuIso)/i.p4.Pt();
93     photonIso[j]=i.pfPhoIso;
94     neutralHadIso[j]=i.pfNeuIso;
95     chargedHadIso[j]=i.pfChaIso;
96 arizzi 1.4 setID(i,j);
97 arizzi 1.1 }
98 arizzi 1.4 template <class Input> void setID(const Input & i, int j)
99     {
100     id[j]=-99;
101     }
102    
103    
104    
105 arizzi 1.1
106     float mass[MAXL]; //MT in case of W
107     float pt[MAXL];
108     float eta[MAXL];
109     float phi[MAXL];
110     float aodCombRelIso[MAXL];
111     float pfCombRelIso[MAXL];
112     float photonIso[MAXL];
113     float neutralHadIso[MAXL];
114     float chargedHadIso[MAXL];
115     float particleIso[MAXL];
116     float dxy[MAXL];
117     float dz[MAXL];
118     int type[MAXL];
119 arizzi 1.4 float id[MAXL];
120 arizzi 1.1 };
121    
122 arizzi 1.9 template <> void LeptonInfo::setID<VHbbEvent::ElectronInfo>(const VHbbEvent::ElectronInfo & i, int j){
123 arizzi 1.4 id[j]=i.id80r;
124     }
125    
126 arizzi 1.1 typedef struct
127     {
128     float et;
129     float sumet;
130     float sig;
131     float phi;
132 arizzi 1.9 } METInfo;
133 arizzi 1.1
134     typedef struct
135     {
136     float mht;
137     float ht;
138     float sig;
139     float phi;
140 arizzi 1.9 } MHTInfo;
141 nmohr 1.5
142     typedef struct
143     {
144     float mass;
145     float pt;
146     float wMass;
147 arizzi 1.9 } TopInfo;
148 arizzi 1.1
149 nmohr 1.7 typedef struct
150 arizzi 1.1 {
151     int run;
152     int lumi;
153     int event;
154 nmohr 1.7 int json;
155 arizzi 1.9 } EventInfo;
156 arizzi 1.1
157     typedef struct
158     {
159     void set(const VHbbEvent::SimpleJet & j, int i)
160     {
161     pt[i]=j.p4.Pt();
162     eta[i]=j.p4.Eta();
163     phi[i]=j.p4.Phi();
164     csv[i]=j.csv;
165 arizzi 1.3 //FIXME: numTracksSV (NEED EDM FIX)
166     //FIXME: chf; float nhf; float cef; float nef; float nch; nconstituents; (NEED EDM FIX)
167    
168 arizzi 1.1 flavour[i]=j.flavour;
169 arizzi 1.3 if(j.bestMCp4.Pt() > 0)
170     {
171     genPt[i]=j.bestMCp4.Pt();
172     genEta[i]=j.bestMCp4.Eta();
173     genPhi[i]=j.bestMCp4.Phi();
174     }
175     //FIXME JECUnc (NEED EDM FIX)
176    
177 arizzi 1.1 }
178     void reset()
179     {
180     for(int i=0;i<MAXJ;i++) {
181     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;
182     }
183     }
184     float pt[MAXJ];
185     float eta[MAXJ];
186     float phi[MAXJ];
187     float csv[MAXJ];
188     float cosTheta[MAXJ];
189     int numTracksSV[MAXJ];
190     float chf[MAXJ];
191     float nhf[MAXJ];
192     float cef[MAXJ];
193     float nef[MAXJ];
194     float nch[MAXJ];
195     float nconstituents[MAXJ];
196     float flavour[MAXJ];
197     float genPt[MAXJ];
198     float genEta[MAXJ];
199     float genPhi[MAXJ];
200     float JECUnc[MAXJ];
201    
202 arizzi 1.9 } JetInfo;
203 arizzi 1.1
204     int main(int argc, char* argv[])
205     {
206     gROOT->Reset();
207    
208     TTree *_outTree;
209 arizzi 1.9 METInfo MET;
210     MHTInfo MHT;
211     TopInfo top;
212     EventInfo EVENT;
213     // JetInfo jet1,jet2, addJet1, addJet2;
214     // lepton1,lepton2;
215     JetInfo hJets, aJets;
216     LeptonInfo vLeptons, aLeptons;
217 arizzi 1.1 int naJets=0, nhJets=0;
218 arizzi 1.9 TrackInfo H;
219     TrackInfo V;
220     int nvlep=0,nalep=0;
221 arizzi 1.1
222 arizzi 1.22 float jjdr,jjdPhi,jjdEta,HVdPhi,HMETdPhi,VMt,deltaPullAngle,deltaPullAngleAK7,gendrcc,gendrbb, genZpt, genWpt, weightTrig, weightTrigMET, minDeltaPhijetMET, jetPt_minDeltaPhijetMET , PUweight;
223 arizzi 1.3 int nofLeptons15,nofLeptons20, Vtype,numJets,numBJets,eventFlav;
224 arizzi 1.1 // bool isMET80_CJ80, ispfMHT150, isMET80_2CJ20,isMET65_2CJ20, isJETID,isIsoMu17;
225 arizzi 1.18 bool triggerFlags[500],hbhe;
226 arizzi 1.17 float btag1T2CSF=1.,btag2TSF=1.,btag1TSF=1.;
227 arizzi 1.1 // ----------------------------------------------------------------------
228     // First Part:
229     //
230     // * enable the AutoLibraryLoader
231     // * book the histograms of interest
232     // * open the input file
233     // ----------------------------------------------------------------------
234    
235     // load framework libraries
236     gSystem->Load( "libFWCoreFWLite" );
237     gSystem->Load("libDataFormatsFWLite");
238     AutoLibraryLoader::enable();
239    
240     // parse arguments
241     if ( argc < 2 ) {
242     return 0;
243     }
244    
245 arizzi 1.4 std::vector<VHbbCandidate> * candZlocal = new std::vector<VHbbCandidate>;
246     std::vector<VHbbCandidate> * candWlocal = new std::vector<VHbbCandidate>;
247    
248 arizzi 1.1 // get the python configuration
249     PythonProcessDesc builder(argv[1]);
250     const edm::ParameterSet& in = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("fwliteInput" );
251     const edm::ParameterSet& out = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("fwliteOutput");
252     const edm::ParameterSet& ana = builder.processDesc()->getProcessPSet()->getParameter<edm::ParameterSet>("Analyzer");
253 nmohr 1.7 std::vector<edm::LuminosityBlockRange> jsonVector;
254     if ( in.exists("lumisToProcess") )
255     {
256     std::vector<edm::LuminosityBlockRange> const & lumisTemp =
257     in.getUntrackedParameter<std::vector<edm::LuminosityBlockRange> > ("lumisToProcess");
258     jsonVector.resize( lumisTemp.size() );
259     copy( lumisTemp.begin(), lumisTemp.end(), jsonVector.begin() );
260     }
261 arizzi 1.1
262     // now get each parameter
263     int maxEvents_( in.getParameter<int>("maxEvents") );
264     unsigned int outputEvery_( in.getParameter<unsigned int>("outputEvery") );
265     std::string outputFile_( out.getParameter<std::string>("fileName" ) );
266    
267    
268     std::vector<std::string> triggers( ana.getParameter<std::vector<std::string> >("triggers") );
269 arizzi 1.6 double btagThr = ana.getParameter<double>("bJetCountThreshold" );
270 arizzi 1.4 bool fromCandidate = ana.getParameter<bool>("readFromCandidates");
271 arizzi 1.17
272     bool useHighestPtHiggsZ = ana.getParameter<bool>("useHighestPtHiggsZ");
273     bool useHighestPtHiggsW = ana.getParameter<bool>("useHighestPtHiggsW");
274     HbbCandidateFinderAlgo * algoZ = new HbbCandidateFinderAlgo(ana.getParameter<bool>("verbose"), ana.getParameter<double>("jetPtThresholdZ"),useHighestPtHiggsZ );
275     HbbCandidateFinderAlgo * algoW = new HbbCandidateFinderAlgo(ana.getParameter<bool>("verbose"), ana.getParameter<double>("jetPtThresholdW"),useHighestPtHiggsW);
276 arizzi 1.4
277 arizzi 1.16 TriggerWeight triggerWeight(ana);
278 arizzi 1.17 BTagWeight btag(2); // 2 operating points "Custom" = 0.5 and "Tight = 0.898"
279     BTagSampleEfficiency btagEff( ana.getParameter<std::string>("btagEffFileName" ).c_str() );
280 arizzi 1.1
281     std::vector<std::string> inputFiles_( in.getParameter<std::vector<std::string> >("fileNames") );
282     // std::string inputFile( in.getParameter<std::string> ("fileName") );
283    
284    
285 arizzi 1.3 std::string PUmcfileName_ = in.getParameter<std::string> ("PUmcfileName") ;
286     std::string PUdatafileName_ = in.getParameter<std::string> ("PUdatafileName") ;
287 arizzi 1.1 bool isMC_( ana.getParameter<bool>("isMC") );
288     TriggerReader trigger(isMC_);
289 arizzi 1.18 TriggerReader patFilters(false);
290 arizzi 1.1
291 arizzi 1.6 // TFile *_outPUFile = new TFile((outputFile_+"_PU").c_str(), "recreate");
292     // TH1F * pu = new TH1F("pileup","",-0.5,24.5,25);
293 arizzi 1.1 TFile *_outFile = new TFile(outputFile_.c_str(), "recreate");
294 arizzi 1.16 TH1F * count = new TH1F("Count","Count", 1,0,2 );
295     TH1F * countWithPU = new TH1F("CountWithPU","CountWithPU", 1,0,2 );
296 arizzi 1.1 _outTree = new TTree("tree", "myTree");
297    
298     _outTree->Branch("H" , &H , "mass/F:pt/F:eta:phi/F");
299     _outTree->Branch("V" , &V , "mass/F:pt/F:eta:phi/F");
300     _outTree->Branch("nhJets" , &nhJets , "nhJets/I");
301     _outTree->Branch("naJets" , &naJets , "naJets/I");
302 arizzi 1.2
303     _outTree->Branch("hJet_pt",hJets.pt ,"pt[nhJets]/F");
304     _outTree->Branch("hJet_eta",hJets.eta ,"eta[nhJets]/F");
305     _outTree->Branch("hJet_phi",hJets.phi ,"phi[nhJets]/F");
306     _outTree->Branch("hJet_csv",hJets.csv ,"csv[nhJets]/F");
307     _outTree->Branch("hJet_cosTheta",hJets.cosTheta ,"cosTheta[nhJets]/F");
308     _outTree->Branch("hJet_numTracksSV",hJets.numTracksSV ,"numTracksSV[nhJets]/I");
309     _outTree->Branch("hJet_chf",hJets.chf ,"chf[nhJets]/F");
310     _outTree->Branch("hJet_nhf",hJets.nhf ,"nhf[nhJets]/F");
311     _outTree->Branch("hJet_cef",hJets.cef ,"cef[nhJets]/F");
312     _outTree->Branch("hJet_nef",hJets.nef ,"nef[nhJets]/F");
313     _outTree->Branch("hJet_nch",hJets.nch ,"nch[nhJets]/F");
314     _outTree->Branch("hJet_nconstituents",hJets.nconstituents ,"nconstituents[nhJets]");
315     _outTree->Branch("hJet_flavour",hJets.flavour ,"flavour[nhJets]/F");
316     _outTree->Branch("hJet_genPt",hJets.genPt ,"genPt[nhJets]/F");
317     _outTree->Branch("hJet_genEta",hJets.genEta ,"genEta[nhJets]/F");
318     _outTree->Branch("hJet_genPhi",hJets.genPhi ,"genPhi[nhJets]/F");
319     _outTree->Branch("hJet_JECUnc",hJets.JECUnc ,"JECUnc[nhJets]/F");
320    
321     _outTree->Branch("aJet_pt",aJets.pt ,"pt[naJets]/F");
322     _outTree->Branch("aJet_eta",aJets.eta ,"eta[naJets]/F");
323     _outTree->Branch("aJet_phi",aJets.phi ,"phi[naJets]/F");
324     _outTree->Branch("aJet_csv",aJets.csv ,"csv[naJets]/F");
325     _outTree->Branch("aJet_cosTheta",aJets.cosTheta ,"cosTheta[naJets]/F");
326     _outTree->Branch("aJet_numTracksSV",aJets.numTracksSV ,"numTracksSV[naJets]/I");
327     _outTree->Branch("aJet_chf",aJets.chf ,"chf[naJets]/F");
328     _outTree->Branch("aJet_nhf",aJets.nhf ,"nhf[naJets]/F");
329     _outTree->Branch("aJet_cef",aJets.cef ,"cef[naJets]/F");
330     _outTree->Branch("aJet_nef",aJets.nef ,"nef[naJets]/F");
331     _outTree->Branch("aJet_nch",aJets.nch ,"nch[naJets]/F");
332     _outTree->Branch("aJet_nconstituents",aJets.nconstituents ,"nconstituents[naJets]");
333     _outTree->Branch("aJet_flavour",aJets.flavour ,"flavour[naJets]/F");
334     _outTree->Branch("aJet_genPt",aJets.genPt ,"genPt[naJets]/F");
335     _outTree->Branch("aJet_genEta",aJets.genEta ,"genEta[naJets]/F");
336     _outTree->Branch("aJet_genPhi",aJets.genPhi ,"genPhi[naJets]/F");
337     _outTree->Branch("aJet_JECUnc",aJets.JECUnc ,"JECUnc[naJets]/F");
338    
339 arizzi 1.1
340     _outTree->Branch("jjdr" , &jjdr , "jjdr/F" );
341     _outTree->Branch("jjdPhi" , &jjdPhi , "jjdPhi/F" );
342 arizzi 1.6 _outTree->Branch("jjdEta" , &jjdEta , "jjdEta/F" );
343 arizzi 1.1 _outTree->Branch("numJets" , &numJets , "numJets/I" );
344     _outTree->Branch("numBJets" , &numBJets , "numBJets/I" );
345     _outTree->Branch("nofLeptons15" , &nofLeptons15 , "nofLeptons15/I" );
346     _outTree->Branch("nofLeptons20" , &nofLeptons20 , "nofLeptons20/I" );
347     _outTree->Branch("deltaPullAngle", &deltaPullAngle , "deltaPullAngle/F");
348     _outTree->Branch("gendrcc" , &gendrcc , "gendrcc/F");
349     _outTree->Branch("gendrbb" , &gendrbb , "gendrbb/F");
350     _outTree->Branch("genZpt" , &genZpt , "genZpt/F");
351     _outTree->Branch("genWpt" , &genWpt , "genWpt/F");
352     _outTree->Branch("weightTrig" , &weightTrig , "weightTrig/F");
353 arizzi 1.22 _outTree->Branch("weightTrigMET" , &weightTrigMET , "weightTrigMET/F");
354 arizzi 1.1 _outTree->Branch("deltaPullAngleAK7", &deltaPullAngleAK7 , "deltaPullAngleAK7/F");
355     _outTree->Branch("PUweight", &PUweight , "PUweight/F");
356 arizzi 1.3 _outTree->Branch("eventFlav", &eventFlav , "eventFlav/I");
357 arizzi 1.1
358    
359    
360    
361     _outTree->Branch("Vtype" , &Vtype , "Vtype/I" );
362     _outTree->Branch("HVdPhi" , &HVdPhi , "HVdPhi/F" );
363 arizzi 1.19 _outTree->Branch("HMETdPhi" , &HMETdPhi , "HMETdPhi/F" );
364 arizzi 1.1 _outTree->Branch("VMt" , &VMt , "VMt/F" );
365    
366 arizzi 1.9 _outTree->Branch("nvlep" , &nvlep , "nvlep/I");
367     _outTree->Branch("nalep" , &nalep , "nalep/I");
368 arizzi 1.3
369 arizzi 1.9 _outTree->Branch("vLepton_mass",vLeptons.mass ,"mass[nvlep]/F");
370     _outTree->Branch("vLepton_pt",vLeptons.pt ,"pt[nvlep]/F");
371     _outTree->Branch("vLepton_eta",vLeptons.eta ,"eta[nvlep]");
372     _outTree->Branch("vLepton_phi",vLeptons.phi ,"phi[nvlep]/F");
373     _outTree->Branch("vLepton_aodCombRelIso",vLeptons.aodCombRelIso ,"aodCombRelIso[nvlep]/F");
374     _outTree->Branch("vLepton_pfCombRelIso",vLeptons.pfCombRelIso ,"pfCombRelIso[nvlep]/F");
375     _outTree->Branch("vLepton_photonIso",vLeptons.photonIso ,"photonIso[nvlep]/F");
376     _outTree->Branch("vLepton_neutralHadIso",vLeptons.neutralHadIso ,"neutralHadIso[nvlep]/F");
377     _outTree->Branch("vLepton_chargedHadIso",vLeptons.chargedHadIso ,"chargedHadIso[nvlep]/F");
378     _outTree->Branch("vLepton_particleIso",vLeptons.particleIso ,"particleIso[nvlep]/F");
379     _outTree->Branch("vLepton_dxy",vLeptons.dxy ,"dxy[nvlep]/F");
380     _outTree->Branch("vLepton_dz",vLeptons.dz ,"dz[nvlep]/F");
381     _outTree->Branch("vLepton_type",vLeptons.type ,"type[nvlep]/I");
382     _outTree->Branch("vLepton_id",vLeptons.id ,"id[nvlep]/F");
383    
384     _outTree->Branch("aLepton_mass",aLeptons.mass ,"mass[nalep]/F");
385     _outTree->Branch("aLepton_pt",aLeptons.pt ,"pt[nalep]/F");
386     _outTree->Branch("aLepton_eta",aLeptons.eta ,"eta[nalep]");
387     _outTree->Branch("aLepton_phi",aLeptons.phi ,"phi[nalep]/F");
388     _outTree->Branch("aLepton_aodCombRelIso",aLeptons.aodCombRelIso ,"aodCombRelIso[nalep]/F");
389     _outTree->Branch("aLepton_pfCombRelIso",aLeptons.pfCombRelIso ,"pfCombRelIso[nalep]/F");
390     _outTree->Branch("aLepton_photonIso",aLeptons.photonIso ,"photonIso[nalep]/F");
391     _outTree->Branch("aLepton_neutralHadIso",aLeptons.neutralHadIso ,"neutralHadIso[nalep]/F");
392     _outTree->Branch("aLepton_chargedHadIso",aLeptons.chargedHadIso ,"chargedHadIso[nalep]/F");
393     _outTree->Branch("aLepton_particleIso",aLeptons.particleIso ,"particleIso[nalep]/F");
394     _outTree->Branch("aLepton_dxy",aLeptons.dxy ,"dxy[nalep]/F");
395     _outTree->Branch("aLepton_dz",aLeptons.dz ,"dz[nalep]/F");
396     _outTree->Branch("aLepton_type",aLeptons.type ,"type[nalep]/I");
397     _outTree->Branch("aLepton_id",aLeptons.id ,"id[nalep]/F");
398    
399 arizzi 1.8 _outTree->Branch("top" , &top , "mass/F:pt/F:wMass/F");
400 arizzi 1.1
401     _outTree->Branch("MET" , &MET , "et/F:sumet:sig/F:phi/F");
402     _outTree->Branch("MHT" , &MHT , "mht/F:ht:sig/F:phi/F");
403     _outTree->Branch("minDeltaPhijetMET" , &minDeltaPhijetMET , "minDeltaPhijetMET/F");
404     _outTree->Branch("jetPt_minDeltaPhijetMET" , &jetPt_minDeltaPhijetMET , "jetPt_minDeltaPhijetMET/F");
405    
406     std::stringstream s;
407     s << "triggerFlags[" << triggers.size() << "]/b";
408     _outTree->Branch("triggerFlags", triggerFlags, s.str().c_str());
409 nmohr 1.7
410 arizzi 1.17
411 nmohr 1.7 _outTree->Branch("EVENT" , &EVENT , "run/I:lumi/I:event/I:json/I");
412 arizzi 1.18 _outTree->Branch("hbhe" , &hbhe , "hbhe/b");
413 arizzi 1.17 _outTree->Branch("btag1TSF" , &btag1TSF , "btag1TSF/F");
414     _outTree->Branch("btag2TSF" , &btag2TSF , "btag2TSF/F");
415     _outTree->Branch("btag1T2CSF" , &btag1T2CSF , "btag1T2CSF/F");
416 arizzi 1.1
417     int ievt=0;
418     int totalcount=0;
419    
420     // TFile* inFile = new TFile(inputFile.c_str(), "read");
421     for(unsigned int iFile=0; iFile<inputFiles_.size(); ++iFile) {
422     std::cout << iFile << std::endl;
423     TFile* inFile = TFile::Open(inputFiles_[iFile].c_str());
424     if(inFile==0) continue;
425    
426     // loop the events
427    
428     fwlite::Event ev(inFile);
429     for(ev.toBegin(); !ev.atEnd(); ++ev, ++ievt)
430     {
431 arizzi 1.16 count->Fill(1.);
432 arizzi 1.1
433     if(isMC_){
434     // PU weights
435 arizzi 1.3
436 arizzi 1.9 // edm::LumiReWeighting LumiWeights_ = edm::LumiReWeighting(PUmcfileName_,PUdatafileName_ , "pileup", "pileup");
437     // double avg=0;
438 arizzi 1.3 //FIXME: PU (NEED EDM FIX)
439 arizzi 1.1 // if( PUintimeSizes.isValid() && PUouttime1minusSizes.isValid() && PUouttime1plusSizes.isValid()){
440     // avg = (double)( *PUintimeSizes );
441     // }
442 arizzi 1.3 PUweight = 1.0; // FIXME: LumiWeights_.weight3BX( avg /3.); (NEED EDM FIX)
443 arizzi 1.1 }
444 arizzi 1.16 countWithPU->Fill(PUweight);
445 nmohr 1.7
446     //Write event info
447     EVENT.run = ev.id().run();
448     EVENT.lumi = ev.id().luminosityBlock();
449     EVENT.event = ev.id().event();
450     EVENT.json = jsonContainsEvent (jsonVector, ev);
451 arizzi 1.3
452 arizzi 1.4
453     const std::vector<VHbbCandidate> * candZ ;
454     const std::vector<VHbbCandidate> * candW ;
455 arizzi 1.13 const VHbbEvent * iEvent =0;
456 arizzi 1.4 if(fromCandidate)
457     {
458 arizzi 1.3 fwlite::Handle< std::vector<VHbbCandidate> > vhbbCandHandleZ;
459     vhbbCandHandleZ.getByLabel(ev,"hbbBestCSVPt20Candidates");
460 arizzi 1.4 candZ = vhbbCandHandleZ.product();
461 arizzi 1.3
462     fwlite::Handle< std::vector<VHbbCandidate> > vhbbCandHandle;
463     vhbbCandHandle.getByLabel(ev,"hbbHighestPtHiggsPt30Candidates");
464 arizzi 1.4 candW = vhbbCandHandle.product();
465     }
466     else
467     {
468     candZlocal->clear();
469     candWlocal->clear();
470     fwlite::Handle< VHbbEvent > vhbbHandle;
471     vhbbHandle.getByLabel(ev,"HbbAnalyzerNew");
472 arizzi 1.13 iEvent = vhbbHandle.product();
473 arizzi 1.4 algoZ->run(vhbbHandle.product(),*candZlocal);
474     algoW->run(vhbbHandle.product(),*candWlocal);
475     candZ= candZlocal;
476     candW= candWlocal;
477    
478    
479 arizzi 1.1
480 arizzi 1.4 }
481 arizzi 1.1
482 arizzi 1.3 const std::vector<VHbbCandidate> * cand = candZ;
483 arizzi 1.1
484    
485     fwlite::Handle< VHbbEventAuxInfo > vhbbAuxHandle;
486     vhbbAuxHandle.getByLabel(ev,"HbbAnalyzerNew");
487     const VHbbEventAuxInfo & aux = *vhbbAuxHandle.product();
488    
489     /* fwlite::Handle< VHbbEvent > vhbbHandle;
490     vhbbHandle.getByLabel(ev,"HbbAnalyzerNew");
491     const VHbbEvent iEvent = *vhbbHandle.product();
492     */
493    
494     // std::clog << "Filling tree "<< std::endl;
495 arizzi 1.17 bool isW=false;
496 arizzi 1.21
497 arizzi 1.17
498 arizzi 1.3 if(cand->size() == 0 or cand->at(0).H.jets.size() < 2) continue;
499     if(cand->size() > 1 )
500     {
501     std::cout << "MULTIPLE CANDIDATES: " << cand->size() << std::endl;
502     }
503 arizzi 1.17 if(cand->at(0).candidateType == VHbbCandidate::Wmun || cand->at(0).candidateType == VHbbCandidate::Wen ) { cand=candW; isW=true; }
504 arizzi 1.3 if(cand->size() == 0)
505     {
506     // std::cout << "W event loss due to tigther cuts" << std::endl;
507     continue;
508     }
509     const VHbbCandidate & vhCand = cand->at(0);
510 arizzi 1.18 patFilters.setEvent(&ev,"VH");
511     hbhe = patFilters.accept("hbhe");
512    
513 arizzi 1.4 trigger.setEvent(&ev);
514     for(size_t j=0;j < triggers.size();j++)
515     triggerFlags[j]=trigger.accept(triggers[j]);
516 arizzi 1.3
517     eventFlav=0;
518     if(aux.mcBbar.size() > 0 || aux.mcB.size() > 0) eventFlav=5;
519     else if(aux.mcC.size() > 0) eventFlav=4;
520 arizzi 1.8
521 arizzi 1.3
522 arizzi 1.1 H.mass = vhCand.H.p4.M();
523     H.pt = vhCand.H.p4.Pt();
524     H.eta = vhCand.H.p4.Eta();
525     H.phi = vhCand.H.p4.Phi();
526     V.mass = vhCand.V.p4.M();
527     V.pt = vhCand.V.p4.Pt();
528     V.eta = vhCand.V.p4.Eta();
529     V.phi = vhCand.V.p4.Phi();
530     nhJets=2;
531     hJets.set(vhCand.H.jets[0],0);
532     hJets.set(vhCand.H.jets[1],1);
533     aJets.reset();
534     naJets=vhCand.additionalJets.size();
535 arizzi 1.6 numBJets=0;
536     if(vhCand.H.jets[0].csv> btagThr) numBJets++;
537     if(vhCand.H.jets[1].csv> btagThr) numBJets++;
538     for( int j=0; j < naJets && j < MAXJ; j++ )
539     {
540     aJets.set(vhCand.additionalJets[j],j);
541     if(vhCand.additionalJets[j].csv> btagThr) numBJets++;
542     }
543 arizzi 1.1 numJets = vhCand.additionalJets.size()+2;
544     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());
545     jjdPhi = deltaPhi(vhCand.H.jets[0].p4.Phi(),vhCand.H.jets[1].p4.Phi());
546 arizzi 1.6 jjdEta= TMath::Abs( vhCand.H.jets[0].p4.Eta() - vhCand.H.jets[1].p4.Eta() );
547 arizzi 1.19 HVdPhi = fabs( deltaPhi(vhCand.H.p4.Phi(),vhCand.V.p4.Phi()) ) ;
548 arizzi 1.20 HMETdPhi = fabs( deltaPhi(vhCand.H.p4.Phi(),vhCand.V.mets.at(0).p4.Phi()) ) ;
549 nmohr 1.14 VMt = vhCand.Mt() ;
550 arizzi 1.1 deltaPullAngle = vhCand.H.deltaTheta;
551 arizzi 1.16
552 arizzi 1.1
553 arizzi 1.3 hJets.cosTheta[0]= vhCand.H.helicities[0];
554     hJets.cosTheta[1]= vhCand.H.helicities[1];
555 arizzi 1.1
556     MET.et = vhCand.V.mets.at(0).p4.Pt();
557     MET.phi = vhCand.V.mets.at(0).p4.Phi();
558     MET.sumet = vhCand.V.mets.at(0).sumEt;
559     MET.sig = vhCand.V.mets.at(0).metSig;
560 arizzi 1.3 //FIXME add MHT _outTree->Branch("MHT" , &MHT , "mht/F:ht:sig/F:phi/F"); (NEED EDM FIX)
561 arizzi 1.1 Vtype = vhCand.candidateType;
562 arizzi 1.9 vLeptons.reset();
563 arizzi 1.15 weightTrig = 1.; // better to default to 1
564 arizzi 1.8 TLorentzVector leptonForTop;
565 arizzi 1.9 size_t firstAddMu=0;
566     size_t firstAddEle=0;
567 arizzi 1.1 if(Vtype == VHbbCandidate::Zmumu ){
568 arizzi 1.13 vLeptons.set(vhCand.V.muons[0],0,13);
569     vLeptons.set(vhCand.V.muons[1],1,13);
570 arizzi 1.16 float cweightID = triggerWeight.scaleMuID(vLeptons.pt[0],vLeptons.eta[0]) * triggerWeight.scaleMuID(vLeptons.pt[1],vLeptons.eta[1]) ;
571     float weightTrig1 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[0],vLeptons.eta[0]);
572     float weightTrig2 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[1],vLeptons.eta[1]);
573 arizzi 1.1 float cweightTrig = weightTrig1 + weightTrig2 - weightTrig1*weightTrig2;
574     weightTrig = cweightID * cweightTrig;
575 arizzi 1.9 nvlep=2;
576     firstAddMu=2;
577 arizzi 1.1 }
578     if( Vtype == VHbbCandidate::Zee ){
579 arizzi 1.9 vLeptons.set(vhCand.V.electrons[0],0,11);
580     vLeptons.set(vhCand.V.electrons[1],1,11);
581     nvlep=2;
582     firstAddEle=2;
583 arizzi 1.3 //FIXME: trigger weights for electrons
584 arizzi 1.1 }
585     if(Vtype == VHbbCandidate::Wmun ){
586 arizzi 1.8 leptonForTop=vhCand.V.muons[0].p4;
587 arizzi 1.13 vLeptons.set(vhCand.V.muons[0],0,13);
588 arizzi 1.16 float cweightID = triggerWeight.scaleMuID(vLeptons.pt[0],vLeptons.eta[0]);
589     float weightTrig1 = triggerWeight.scaleMuIsoHLT(vLeptons.pt[0],vLeptons.eta[0]);
590 arizzi 1.1 float cweightTrig = weightTrig1;
591     weightTrig = cweightID * cweightTrig;
592 arizzi 1.9 nvlep=1;
593     firstAddMu=1;
594 arizzi 1.1 }
595     if( Vtype == VHbbCandidate::Wen ){
596 arizzi 1.8 leptonForTop=vhCand.V.electrons[0].p4;
597 arizzi 1.9 vLeptons.set(vhCand.V.electrons[0],0,11);
598     nvlep=1;
599     firstAddEle=1;
600 arizzi 1.1 }
601 arizzi 1.3 if( Vtype == VHbbCandidate::Znn ){
602 arizzi 1.9 nvlep=0;
603 arizzi 1.22 float weightTrig1 = triggerWeight.scaleMetHLT(vhCand.V.mets.at(0).p4.Pt());
604     weightTrig = weightTrig1;
605     weightTrigMET = weightTrig1;
606 arizzi 1.3 }
607 arizzi 1.9
608     aLeptons.reset();
609     nalep=0;
610 arizzi 1.21 if(fromCandidate)
611     {
612     for(size_t j=firstAddMu;j< vhCand.V.muons.size();j++) aLeptons.set(vhCand.V.muons[j],nalep++,13);
613     for(size_t j=firstAddEle;j< vhCand.V.electrons.size();j++) aLeptons.set(vhCand.V.electrons[j],nalep++,11);
614     }
615     else
616     {
617     for(size_t j=0;j< iEvent->muInfo.size();j++)
618     {
619     if((j!= vhCand.V.firstLepton && j!= vhCand.V.secondLepton) || ((Vtype != VHbbCandidate::Wmun ) && (Vtype != VHbbCandidate::Zmumu )) )
620     aLeptons.set(iEvent->muInfo[j],nalep++,13);
621     }
622     for(size_t j=0;j< iEvent->eleInfo.size();j++)
623     {
624     if((j!= vhCand.V.firstLepton && j!= vhCand.V.secondLepton) || ((Vtype != VHbbCandidate::Wen ) && (Vtype != VHbbCandidate::Zee )))
625     aLeptons.set(iEvent->eleInfo[j],nalep++,11);
626     }
627 arizzi 1.8
628 arizzi 1.21 }
629 arizzi 1.8
630 arizzi 1.17 //Loop on jets
631    
632 arizzi 1.8 double maxBtag=-99999;
633 arizzi 1.16 minDeltaPhijetMET = 999;
634 arizzi 1.8 TLorentzVector bJet;
635 arizzi 1.17 std::vector<std::vector<BTagWeight::JetInfo> > btagJetInfos;
636 arizzi 1.8 for(unsigned int j=0; j < vhCand.H.jets.size(); j++ ){
637     if (vhCand.H.jets[j].csv > maxBtag) { bJet=vhCand.H.jets[j].p4 ; maxBtag =vhCand.H.jets[j].csv; }
638 arizzi 1.16 if (deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.H.jets[j].p4.Phi()) < minDeltaPhijetMET)
639     {
640     minDeltaPhijetMET=deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.H.jets[j].p4.Phi());
641     jetPt_minDeltaPhijetMET=vhCand.H.jets[j].p4.Pt();
642     }
643 arizzi 1.17 btagJetInfos.push_back(btagEff.jetInfo(vhCand.H.jets[j]));
644 arizzi 1.8 }
645     for(unsigned int j=0; j < vhCand.additionalJets.size(); j++ ){
646     if (vhCand.additionalJets[j].csv > maxBtag) { bJet=vhCand.additionalJets[j].p4 ; maxBtag =vhCand.additionalJets[j].csv; }
647 arizzi 1.16 if (deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.additionalJets[j].p4.Phi()) < minDeltaPhijetMET)
648     {
649     minDeltaPhijetMET=deltaPhi( vhCand.V.mets.at(0).p4.Phi(), vhCand.additionalJets[j].p4.Phi());
650     jetPt_minDeltaPhijetMET=vhCand.additionalJets[j].p4.Pt();
651     }
652 arizzi 1.17 if( ( isW && ! useHighestPtHiggsW ) || ( ! isW && ! useHighestPtHiggsZ ) ) // btag SF computed using only H-jets if best-H made with dijetPt rather than best CSV
653     {
654 arizzi 1.22 if(vhCand.additionalJets[j].p4.Pt() > 20)
655 arizzi 1.17 btagJetInfos.push_back(btagEff.jetInfo(vhCand.additionalJets[j]));
656     }
657     }
658     if(isMC_)
659     {
660     //std::cout << "BTAGSF " << btagJetInfos.size() << " " << btag.weight<BTag1Tight2CustomFilter>(btagJetInfos) << std::endl;
661     btag1T2CSF = btag.weight<BTag1Tight2CustomFilter>(btagJetInfos);
662 arizzi 1.23 btag2TSF = btag.weight<BTag2TightFilter>(btagJetInfos);
663     btag1TSF = btag.weight<BTag1TightFilter>(btagJetInfos);;
664 arizzi 1.16 }
665 arizzi 1.17
666 arizzi 1.8 if(maxBtag > -99999)
667     {
668     TopHypo topQuark = TopMassReco::topMass(leptonForTop,bJet,vhCand.V.mets.at(0).p4);
669     top.mass = topQuark.p4.M();
670     top.pt = topQuark.p4.Pt();
671     top.wMass = topQuark.p4W.M();
672     } else {
673     top.mass = -99;
674     top.pt = -99;
675     top.wMass = -99;
676     }
677    
678    
679    
680 arizzi 1.1 //FIXME _outTree->Branch("nofLeptons15" , &nofLeptons15 , "nofLeptons15/I" );
681     nofLeptons20= vhCand.additionalLeptons();
682     // if(aux.mcC.size() >=2)
683     // std::cout << "C Must not be zero and it is ... " << aux.mcC[1].p4.Pt() << std::endl;
684     // if(aux.mcB.size() >=1)
685     // std::cout << "B Must not be zero and it is ... " << aux.mcB[0].p4.Pt() << std::endl;
686    
687 arizzi 1.3 // FIXME gendrcc=aux.genCCDeltaR(); (NEED EDM FIX)
688    
689     // FIXME gendrbb=aux.genBBDeltaR(); (NEED EDM FIX)
690 arizzi 1.1 genZpt=aux.mcZ.size() > 0 ? aux.mcZ[0].p4.Pt():-99;
691     genWpt=aux.mcW.size() > 0 ? aux.mcW[0].p4.Pt():-99;
692 arizzi 1.3
693 arizzi 1.1 //FIXME: _outTree->Branch("deltaPullAngleAK7", &deltaPullAngleAK7 , "deltaPullAngleAK7/F");
694    
695    
696    
697 arizzi 1.17
698 arizzi 1.1 _outTree->Fill();
699 arizzi 1.4
700 arizzi 1.1 }// closed event loop
701    
702     std::cout << "closing the file: " << inputFiles_[iFile] << std::endl;
703     inFile->Close();
704     // close input file
705     } // loop on files
706    
707    
708     std::cout << "Events: " << ievt <<std::endl;
709     std::cout << "TotalCount: " << totalcount <<std::endl;
710    
711    
712    
713     _outFile->cd();
714    
715     _outTree->Write();
716     _outFile->Write();
717     _outFile->Close();
718     return 0;
719     }
720    
721