ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbbAnalysis/VHbbDataFormats/bin/Ntupler.cc
Revision: 1.8
Committed: Thu Sep 8 10:32:32 2011 UTC (13 years, 8 months ago) by arizzi
Content type: text/plain
Branch: MAIN
Changes since 1.7: +34 -12 lines
Log Message: 
interface/TopMassReco.h

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