ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/VHbbAnalysis/VHbbDataFormats/bin/Ntupler.cc
Revision: 1.3
Committed: Wed Sep 7 12:27:45 2011 UTC (13 years, 8 months ago) by arizzi
Content type: text/plain
Branch: MAIN
Changes since 1.2: +68 -54 lines
Log Message:
update of ntupler

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