ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbbAnalysis/VHbbDataFormats/bin/Ntupler.cc
Revision: 1.4
Committed: Wed Sep 7 16:40:36 2011 UTC (13 years, 8 months ago) by arizzi
Content type: text/plain
Branch: MAIN
Changes since 1.3: +86 -33 lines
Log Message:
options to read from the VHEvent

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