ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbbAnalysis/VHbbDataFormats/bin/Ntupler.cc
Revision: 1.11
Committed: Thu Sep 8 13:02:03 2011 UTC (13 years, 8 months ago) by arizzi
Content type: text/plain
Branch: MAIN
Changes since 1.10: +8 -10 lines
Log Message:
config file for t and p

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