ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/algomez/FourTop/TreeAnalyzer/test/Analyzer.C
Revision: 1.11
Committed: Mon Jan 23 18:25:00 2012 UTC (13 years, 3 months ago) by algomez
Content type: text/plain
Branch: MAIN
Changes since 1.10: +57 -15 lines
Log Message:
include different B-tagging operating points for CSV

File Contents

# User Rev Content
1 algomez 1.1
2     // The class definition in Analyzer.h has been generated automatically
3     // by the ROOT utility TTree::MakeSelector(). This class is derived
4     // from the ROOT class TSelector. For more information on the TSelector
5     // framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.
6    
7     // The following methods are defined in this file:
8     // Begin(): called every time a loop on the tree starts,
9     // a convenient place to create your histograms.
10     // SlaveBegin(): called after Begin(), when on PROOF called only on the
11     // slave servers.
12     // Process(): called for each event, in this function you decide what
13     // to read and fill your histograms.
14     // SlaveTerminate: called at the end of the loop on the tree, when on PROOF
15     // called only on the slave servers.
16     // Terminate(): called at the end of the loop on the tree,
17     // a convenient place to draw/fit your histograms.
18     //
19     // To use this file, try the following session on your Tree T:
20     //
21     // Root > T->Process("Analyzer.C")
22     // Root > T->Process("Analyzer.C","some options")
23     // Root > T->Process("Analyzer.C+")
24     //
25    
26     #include "Analyzer.h"
27 algomez 1.9 //#include "BTagWeight.h"
28     #include "BTagSFUtil_lite.h"
29 algomez 1.1 //#include "Yumiceva/TreeAnalyzer/interface/JetCombinatorics.h"
30    
31     #include <TStyle.h>
32     #include <TSystem.h>
33    
34     #include <iostream>
35     #include <vector>
36     #include <stdio.h>
37     #include <stdlib.h>
38 algomez 1.8 #define NMAX 100
39 algomez 1.1
40     void Analyzer::ParseInput()
41     {
42    
43     if (fMyOpt.Contains("muon"))
44     {
45     fChannel = 1;
46     }
47     if (fMyOpt.Contains("electron"))
48     {
49     fChannel = 2;
50     }
51     if (fMyOpt.Contains("verbose"))
52     {
53     fVerbose = true;
54     }
55     if (fMyOpt.Contains("JECUP")) { fdoJECunc = true; fdoJECup = true; }
56     if (fMyOpt.Contains("JECDOWN")) { fdoJECunc = true; fdoJECup = false; }
57 algomez 1.2 if (fMyOpt.Contains("PUUP")) { fpuhistogram = "WHistUp";}
58     if (fMyOpt.Contains("PUDOWN")) { fpuhistogram = "WHistDown";}
59 algomez 1.8 if (fMyOpt.Contains("QCD1")) fdoQCD1SideBand = true;//anti-isolation
60     if (fMyOpt.Contains("QCD2")) fdoQCD2SideBand = true;//MET < 20 GeV
61 algomez 1.1 if (fMyOpt.Contains("mtop")) fdoMtopCut = true;
62     if (fMyOpt.Contains("outdir")) {
63     TString tmp = fMyOpt;
64     tmp = tmp.Remove(0,fMyOpt.Index("outdir")+7);
65     fOutdir = tmp+"/";
66     Info("Begin","output files will be written to directory: %s", fOutdir.Data());
67     }
68     if (fMyOpt.Contains("sample"))
69     {
70     TString tmp = fMyOpt;
71     tmp = tmp.Remove(0,fMyOpt.Index("sample")+7);
72     fSample = tmp;
73     if (fdoJECunc && fdoJECup==true) fSample += "_JECUP";
74     if (fdoJECunc && fdoJECup==false) fSample += "_JECDOWN";
75 algomez 1.2 if (fMyOpt.Contains("PUUP")) fSample += "_PUUP";
76     if (fMyOpt.Contains("PUDOWN")) fSample += "_PUDOWN";
77 algomez 1.1
78     Info("Begin","Histogram names will have suffix: %s", fSample.Data());
79    
80     if ( fSample.Contains("data") )
81     {
82     fIsMC = false;
83     Info("Begin","This sample is treated as DATA");
84     if (fdoQCD1SideBand) fSample = "dataQCD1";
85     if (fdoQCD2SideBand) fSample = "dataQCD2";
86     }
87 algomez 1.2 else
88     {
89     Info("Begin","This sample is treated as MC");
90     Info("Begin","we will use pileup set: %s", fpuhistogram.Data());
91     }
92    
93 algomez 1.1
94     }
95     }
96    
97     void Analyzer::Begin(TTree * /*tree*/)
98     {
99     // The Begin() function is called at the start of the query.
100     // When running with PROOF Begin() is only called on the client.
101     // The tree argument is deprecated (on PROOF 0 is passed).
102    
103     TString option = GetOption();
104     fMyOpt = option;
105     ParseInput();
106    
107     Info("Begin", "starting with process option: %s", option.Data());
108    
109     }
110    
111     void Analyzer::SlaveBegin(TTree * tree)
112     {
113    
114     // The SlaveBegin() function is called after the Begin() function.
115     // When running with PROOF SlaveBegin() is called on each slave server.
116     // The tree argument is deprecated (on PROOF 0 is passed).
117    
118     TString option = GetOption();
119     fMyOpt = option;
120     ParseInput();
121    
122     Info("SlaveBegin",
123     "starting with process option: %s (tree: %p)", option.Data(), tree);
124    
125     //initialize the Tree branch addresses
126     Init(tree);
127    
128     // We may be creating a dataset or a merge file: check it
129     TNamed *nm = dynamic_cast<TNamed *>(fInput->FindObject("SimpleNtuple.root"));
130     if (nm) {
131     // Just create the object
132     UInt_t opt = TProofOutputFile::kRegister | TProofOutputFile::kOverwrite | TProofOutputFile::kVerify;
133 algomez 1.8 fProofFile = new TProofOutputFile("SimpleNtuple.root", TProofOutputFile::kDataset, opt, nm->GetTitle());
134 algomez 1.1 } else {
135     // For the ntuple, we use the automatic file merging facility
136     // Check if an output URL has been given
137     TNamed *out = (TNamed *) fInput->FindObject("PROOF_OUTPUTFILE_LOCATION");
138     Info("SlaveBegin", "PROOF_OUTPUTFILE_LOCATION: %s", (out ? out->GetTitle() : "undef"));
139     TString tmpfilename = "results";
140     if ( fSample != "" ) tmpfilename += "_"+fSample+".root";
141     else tmpfilename = "results.root";
142     fProofFile = new TProofOutputFile(tmpfilename, (out ? out->GetTitle() : "M"));
143     out = (TNamed *) fInput->FindObject("PROOF_OUTPUTFILE");
144     if (out) fProofFile->SetOutputFileName(fOutdir + out->GetTitle());
145     }
146    
147     // Open the file
148     //TDirectory *savedir = gDirectory;
149     if (!(fFile = fProofFile->OpenFile("RECREATE"))) {
150     Warning("SlaveBegin", "problems opening file: %s/%s",
151     fProofFile->GetDir(), fProofFile->GetFileName());
152     }
153    
154 algomez 1.2 // Get PU weights
155 algomez 1.8 TString weightfilename = "/uscms/home/yumiceva/work/CMSSW_4_2_4/src/Yumiceva/TreeAnalyzer/test/Weight3Dfinebin4p7.root"; //Weight3D.root";
156 algomez 1.2 fweightfile = new TFile(weightfilename,"read");
157     f3Dweight = (TH1D*) fweightfile->Get(fpuhistogram);
158    
159 algomez 1.1 //create histograms
160     h1test = new TH1F("h1test","muon p_{T}",100,10.,400);
161     //fHist = new HistoManager(string(fSample));
162     TString hname = "_"+fSample;
163    
164 algomez 1.6 hmuons["N_cut0"] = new TH1F("N_muons_cut0"+hname,"Number of Muons",6,-0.5,5.5);
165     hmuons["N"] = new TH1F("N_muons"+hname,"Number of Muons",6,-0.5,5.5);
166 algomez 1.10 hmuons["N_1btag"] = new TH1F("N_muons_1btag"+hname,"Number of Muons",6,-0.5,5.5);
167 algomez 1.6 hmuons["Nelectrons_cut0"] = new TH1F("Nelectrons_cut0"+hname,"Number of Loose Electrons",6,-0.5,5.5);
168     hmuons["Nelectrons"] = new TH1F("Nelectrons"+hname,"Number of Loose Electrons",6,-0.5,5.5);
169 algomez 1.10 hmuons["Nelectrons_1btag"] = new TH1F("Nelectrons_1btag"+hname,"Number of Loose Electrons",6,-0.5,5.5);
170 algomez 1.1 hmuons["pt_1jet"] = new TH1F("muon_pt_1jet"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
171     hmuons["pt_2jet"] = new TH1F("muon_pt_2jet"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
172 algomez 1.2 hmuons["pt_fin"] = new TH1F("muon_pt_fin"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
173 algomez 1.1 hmuons["eta"] = new TH1F("muon_eta"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
174     hmuons["eta_1jet"] = new TH1F("muon_eta_1jet"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
175     hmuons["eta_2jet"] = new TH1F("muon_eta_2jet"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
176 algomez 1.2 hmuons["eta_fin"] = new TH1F("muon_eta_fin"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
177 algomez 1.1 hmuons["phi"] = new TH1F("muon_phi"+hname,"#phi^{#mu}", 30, -3.15, 3.15);
178     hmuons["phi_1jet"] = new TH1F("muon_phi_1jet"+hname,"#phi^{#mu}", 30, -3.15, 3.15);
179     hmuons["phi_2jet"] = new TH1F("muon_phi_2jet"+hname,"#phi^{#mu}", 30, -3.15, 3.15);
180 algomez 1.2 hmuons["phi_fin"] = new TH1F("muon_phi_fin"+hname,"#phi^{#mu}", 30, -3.15, 3.15);
181 algomez 1.1 hmuons["reliso"] = new TH1F("muon_reliso"+hname,"Relative Isolation", 40, 0, 0.2);
182 algomez 1.2 hmuons["reliso_fin"] = new TH1F("muon_reliso_fin"+hname,"Relative Isolation", 40, 0, 0.2);
183 algomez 1.1 hmuons["deltaR_cut0"] = new TH1F("deltaR_cut0"+hname,"#DeltaR(#mu,jet)",80, 0, 4);
184     hmuons["deltaR"] = new TH1F("deltaR"+hname,"#DeltaR(#mu,jet)", 80, 0, 4);
185 algomez 1.10 hmuons["deltaR_1btag"] = new TH1F("deltaR_1btag"+hname,"#DeltaR(#mu,jet)", 80, 0, 4);
186 algomez 1.1 hmuons["d0_cut1"] = new TH1F("d0_cut1"+hname,"#mu Impact Parameter [cm]",20,-0.1,0.1);
187     hmuons["pt_cut1"] = new TH1F("muon_pt_cut1"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
188     hmuons["pt_cut2"] = new TH1F("muon_pt_cut2"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
189     hmuons["pt_cut3"] = new TH1F("muon_pt_cut3"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
190     hmuons["pt"] = new TH1F("muon_pt"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
191 algomez 1.10 hmuons["pt_1btag"] = new TH1F("muon_pt_1btag"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
192 algomez 1.1 hmuons["dz"] = new TH1F("dz"+hname,"|z(#mu) - z_{PV}| [cm]", 25, 0, 1.);
193 algomez 1.2 hmuons["Niso"] = new TH1F("Niso"+hname,"Number of Primary Vertices", 25,-0.5,24.5);
194     hmuons["Ngood"] = new TH1F("Ngood"+hname,"Number of Primary Vertices",25,-0.5,24.5);
195 algomez 1.1
196     hPVs["N"] = new TH1F("NPV"+hname,"Number of PVs",25,-0.5,24.5);
197     hPVs["Nreweight"] = new TH1F("NPVreweight"+hname,"Number of PVs",25,-0.5,24.5);
198     hPVs["Nreweight_2jet"] = new TH1F("NPVreweight_2jet"+hname,"Number of PVs",25,-0.5,24.5);
199 algomez 1.10 hPVs["Nreweight_2jet_1btag"] = new TH1F("NPVreweight_2jet_1btag"+hname,"Number of PVs",25,-0.5,24.5);
200 algomez 1.1
201     helectrons["pt"] = new TH1F("electron_pt"+hname,"p_{T}^{#mu} [GeV/c]", 50, 20, 500);
202     helectrons["pt_cut2"] = new TH1F("electron_pt_cut2"+hname,"p_{T}^{#mu} [GeV/c]", 50, 20, 500);
203     helectrons["pt_1jet"] = new TH1F("electron_pt_1jet"+hname,"p_{T}^{#mu} [GeV/c]", 50, 20, 500);
204     helectrons["pt_2jet"] = new TH1F("electron_pt_2jet"+hname,"p_{T}^{#mu} [GeV/c]", 50, 20, 500);
205     helectrons["pt_3jet"] = new TH1F("electron_pt_3jet"+hname,"p_{T}^{#mu} [GeV/c]", 50, 20, 500);
206     helectrons["pt_4jet"] = new TH1F("electron_pt_4jet"+hname,"p_{T}^{#mu} [GeV/c]", 50, 20, 500);
207     helectrons["eta"] = new TH1F("electron_eta"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
208     helectrons["eta_cut2"] = new TH1F("electron_eta_cut2"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
209     helectrons["eta_1jet"] = new TH1F("electron_eta_1jet"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
210     helectrons["eta_2jet"] = new TH1F("electron_eta_2jet"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
211     helectrons["eta_3jet"] = new TH1F("electron_eta_3jet"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
212     helectrons["eta_4jet"] = new TH1F("electron_eta_4jet"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
213     helectrons["phi"] = new TH1F("electron_phi"+hname,"#phi^{#mu}", 20, 0, 3.15);
214     helectrons["phi_cut2"] = new TH1F("electron_phi_cut2"+hname,"#phi^{#mu}", 20, 0, 3.15);
215     helectrons["phi_1jet"] = new TH1F("electron_phi_1jet"+hname,"#phi^{#mu}", 20, 0, 3.15);
216     helectrons["phi_2jet"] = new TH1F("electron_phi_2jet"+hname,"#phi^{#mu}", 20, 0, 3.15);
217     helectrons["phi_3jet"] = new TH1F("electron_phi_3jet"+hname,"#phi^{#mu}", 20, 0, 3.15);
218     helectrons["phi_4jet"] = new TH1F("electron_phi_4jet"+hname,"#phi^{#mu}", 20, 0, 3.15);
219     helectrons["reliso"] = new TH1F("electron_reliso"+hname,"Relative Isolation", 40, 0, 0.2);
220     helectrons["reliso_1jet"] = new TH1F("electron_reliso_1jet"+hname,"Relative Isolation", 40, 0, 0.2);
221     helectrons["reliso_2jet"] = new TH1F("electron_reliso_2jet"+hname,"Relative Isolation", 40, 0, 0.2);
222     helectrons["reliso_3jet"] = new TH1F("electron_reliso_3jet"+hname,"Relative Isolation", 40, 0, 0.2);
223     helectrons["reliso_4jet"] = new TH1F("electron_reliso_4jet"+hname,"Relative Isolation", 40, 0, 0.2);
224     helectrons["deltaR_cut0"] = new TH1F("electron_deltaR_cut0"+hname,"#DeltaR(#mu,jet)",80, 0, 4);
225     helectrons["deltaR"] = new TH1F("electron_deltaR"+hname,"#DeltaR(#mu,jet)",80, 0, 4);
226     helectrons["d0_cut1"] = new TH1F("electron_d0_cut1"+hname,"#mu Impact Parameter [cm]",20,-0.1,0.1);
227     helectrons["dz"] = new TH1F("electron_dz"+hname,"|z(#mu) - z_{PV}| [cm]", 25, 0, 1.);
228    
229 algomez 1.6 hMET["MET"] = new TH1F("MET"+hname,"Missing Transverse Energy [GeV]", 50, 0, 500);
230 algomez 1.10 hMET["MET_1btag"] = new TH1F("MET_1btag"+hname,"Missing Transverse Energy [GeV]", 50, 0, 500);
231 algomez 1.6 hMET["MET_cut0"] = new TH1F("MET_cut0"+hname,"Missing Transverse Energy [GeV]", 50, 0, 300);
232 algomez 1.1 hMET["MET_2jet"] = new TH1F("MET_2jet"+hname,"Missing Transverse Energy [GeV]", 50, 0, 300);
233 algomez 1.2 hMET["MET_fin"] = new TH1F("MET_fin"+hname,"Missing Transverse Energy [GeV]", 50, 0, 300);
234 algomez 1.1 hMET["genMET_2jet"] = new TH1F("genMET_2jet"+hname,"Missing Transverse Energy [GeV]", 50, 0, 300);
235     hMET["deltaMET_2jet"] = new TH1F("deltaMET_2jet"+hname,"Missing Transverse Energy [GeV]", 50, -200, 200);
236     hMET["phi"] = new TH1F("MET_phi"+hname,"#phi Missing Transverse Energy [GeV]", 20, 0, 3.15);
237 algomez 1.6 hMET["Ht0"] = new TH1F("Ht0"+hname,"H_{T} [GeV]", 50, 0, 3000);
238     hMET["Ht"] = new TH1F("Ht"+hname,"H_{T} [GeV]", 50, 0, 3000);
239 algomez 1.10 hMET["Ht_1btag"] = new TH1F("Ht_1btag"+hname,"H_{T} [GeV]", 50, 0, 3000);
240 algomez 1.8 hMET["Htlep"] = new TH1F("Htlep"+hname,"H_{T,lep} [GeV]", 100, 0, 1000);
241 algomez 1.10 hMET["Stlep"] = new TH1F("Stlep"+hname,"S_{t}^{lep} [GeV]", 100, 0, 1000);
242     hMET["Stlep_1btag"] = new TH1F("Stlep_1btag"+hname,"S_{t}^{lep} [GeV]", 100, 0, 1000);
243     hMET["Stjet"] = new TH1F("Stjet"+hname,"S_{t}^{jet} [GeV]", 100, 0, 3000);
244     hMET["Stjet_1btag"] = new TH1F("Stjet_1btag"+hname,"S_{t}^{jet} [GeV]", 100, 0, 3000);
245 algomez 1.1 hMET["PzNu"] = new TH1F("PzNu"+hname,"p_{z} #nu [GeV]", 40, -300,300);
246     hMET["EtaNu"] = new TH1F("EtaNu"+hname,"#eta",50,-2.2,2.2);
247     hMET["LepWmass"] = new TH1F("LepWmass"+hname,"W#rightarrow#mu#nu Mass [GeV/c^{2}]",20, 0, 150);
248     hMET["LepWmass_topcut"] = new TH1F("LepWmass_topcut"+hname,"W#rightarrow#mu#nu Mass [GeV/c^{2}]",20, 0, 150);
249     hMET["LepWmassNoPt"]=new TH1F("LepWmassNoPt"+hname,"W#rightarrow#mu#nu Mass [GeV/c^{2}]",20, 0, 150);
250     hMET["deltaPhi"] = new TH1F("deltaPhi"+hname,"#Delta #phi(#mu,MET)",50, -3.15, 3.15);
251    
252 algomez 1.6 hM["WMt"] = new TH1F("Mt"+hname,"M_{T}(W) [GeV/c^{2}]", 50, 0, 500); // Transverse Mass sqrt(Wpt*Wpt - Wpx*Wpx - Wpy*Wpy)
253 algomez 1.1 hM["WMt_2jet"] = new TH1F("Mt_2jet"+hname,"M_{T}(W) [GeV/c^{2}]", 50, 0, 300);
254 algomez 1.10 hM["WMt_2jet_1btag"] = new TH1F("Mt_2jet_1btag"+hname,"M_{T}(W) [GeV/c^{2}]", 50, 0, 300);
255 algomez 1.2 hM["WMt_fin"] = new TH1F("Mt_fin"+hname,"M_{T}(W) [GeV/c^{2}]", 50, 0, 300);
256 algomez 1.1 hM["dijet"] = new TH1F("dijet"+hname,"(jj) mass [GeV/c^{2}]", 50, 0, 1000);
257     hM["trijet"] = new TH1F("trijet"+hname,"(jjj) mass [GeV/c^{2}]", 50, 0, 1000);
258     hM["top_1btag"] = new TH1F("top_1btag"+hname,"top mass [GeV/c^{2}]",50,0,500);
259     hM["top_2btag"] = new TH1F("top_2btag"+hname,"top mass [GeV/c^{2}]",50,0,500);
260 algomez 1.6
261 algomez 1.1 hjets["pt"] = new TH1F("jet_pt"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
262 algomez 1.2 hjets["pt_b_mc"] = new TH1F("jet_pt_b_mc"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
263     hjets["pt_c_mc"] = new TH1F("jet_pt_c_mc"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
264     hjets["pt_l_mc"] = new TH1F("jet_pt_l_mc"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
265 algomez 1.1 hjets["pt_btag"] = new TH1F("jet_pt_btag"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
266     hjets["pt_btag_b"] = new TH1F("jet_pt_btag_b"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
267     hjets["pt_btag_c"] = new TH1F("jet_pt_btag_c"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
268     hjets["pt_btag_l"] = new TH1F("jet_pt_btag_l"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
269 algomez 1.6 hjets["1st_pt_cut0"] = new TH1F("jet1_pt_cut0"+hname,"Leading jet p_{T} [GeV/c]",60, 30, 1000);
270     hjets["1st_pt"] = new TH1F("jet1_pt"+hname,"Leading jet p_{T} [GeV/c]",60, 30, 1000);
271 algomez 1.10 hjets["1st_pt_1btag"] = new TH1F("jet1_pt_1btag"+hname,"Leading jet p_{T} [GeV/c]",60, 30, 1000);
272 algomez 1.6 hjets["2nd_pt_cut0"] = new TH1F("jet2_pt_cut0"+hname,"2nd jet p_{T} [GeV/c]",60, 30, 1000);
273     hjets["2nd_pt"] = new TH1F("jet2_pt"+hname,"2nd jet p_{T} [GeV/c]",60, 30, 1000);
274 algomez 1.10 hjets["2nd_pt_1btag"] = new TH1F("jet2_pt_1btag"+hname,"2nd jet p_{T} [GeV/c]",60, 30, 1000);
275 algomez 1.6 hjets["3rd_pt_cut0"] = new TH1F("jet3_pt_cut0"+hname,"3rd jet p_{T} [GeV/c]",60, 30, 800);
276 algomez 1.1 hjets["3rd_pt"] = new TH1F("jet3_pt"+hname,"3rd jet p_{T} [GeV/c]",60, 30, 800);
277 algomez 1.10 hjets["3rd_pt_1btag"] = new TH1F("jet3_pt_1btag"+hname,"3rd jet p_{T} [GeV/c]",60, 30, 800);
278 algomez 1.6 hjets["4th_pt_cut0"] = new TH1F("jet4_pt_cut0"+hname,"4th jet p_{T} [GeV/c]",60, 30, 800);
279 algomez 1.1 hjets["4th_pt"] = new TH1F("jet4_pt"+hname,"4th jet p_{T} [GeV/c]",60, 30, 800);
280 algomez 1.10 hjets["4th_pt_1btag"] = new TH1F("jet4_pt_1btag"+hname,"4th jet p_{T} [GeV/c]",60, 30, 800);
281 algomez 1.6 hjets["5th_pt_cut0"] = new TH1F("jet5_pt_cut0"+hname,"5th jet p_{T} [GeV/c]",60, 30, 800);
282     hjets["5th_pt"] = new TH1F("jet5_pt"+hname,"5th jet p_{T} [GeV/c]",60, 30, 800);
283 algomez 1.10 hjets["5th_pt_1btag"] = new TH1F("jet5_pt_1btag"+hname,"5th jet p_{T} [GeV/c]",60, 30, 800);
284 algomez 1.6 hjets["6th_pt_cut0"] = new TH1F("jet6_pt_cut0"+hname,"6th jet p_{T} [GeV/c]",60, 30, 800);
285     hjets["6th_pt"] = new TH1F("jet6_pt"+hname,"6th jet p_{T} [GeV/c]",60, 30, 800);
286 algomez 1.10 hjets["6th_pt_1btag"] = new TH1F("jet6_pt_1btag"+hname,"6th jet p_{T} [GeV/c]",60, 30, 800);
287 algomez 1.6 hjets["7th_pt_cut0"] = new TH1F("jet7_pt_cut0"+hname,"7th jet p_{T} [GeV/c]",60, 30, 800);
288     hjets["7th_pt"] = new TH1F("jet7_pt"+hname,"7th jet p_{T} [GeV/c]",60, 30, 800);
289 algomez 1.10 hjets["7th_pt_1btag"] = new TH1F("jet7_pt_1btag"+hname,"7th jet p_{T} [GeV/c]",60, 30, 800);
290 algomez 1.2 hjets["1st_pt_fin"] = new TH1F("jet1_pt_fin"+hname,"Leading jet p_{T} [GeV/c]",60, 30, 800);
291     hjets["2nd_pt_fin"] = new TH1F("jet2_pt_fin"+hname,"2nd jet p_{T} [GeV/c]",60, 30, 800);
292 algomez 1.1 hjets["eta"] = new TH1F("jet_eta"+hname,"jet #eta",50, -2.4, 2.4);
293 algomez 1.6 hjets["1st_eta_cut0"] = new TH1F("jet1_eta_cut0"+hname,"Leading jet #eta",50, -2.4, 2.4);
294 algomez 1.1 hjets["1st_eta"] = new TH1F("jet1_eta"+hname,"Leading jet #eta",50, -2.4, 2.4);
295 algomez 1.10 hjets["1st_eta_1btag"] = new TH1F("jet1_eta_1btag"+hname,"Leading jet #eta",50, -2.4, 2.4);
296 algomez 1.6 hjets["2nd_eta_cut0"] = new TH1F("jet2_eta_cut0"+hname,"2nd jet #eta",50, -2.4, 2.4);
297 algomez 1.1 hjets["2nd_eta"] = new TH1F("jet2_eta"+hname,"2nd jet #eta",50, -2.4, 2.4);
298 algomez 1.10 hjets["2nd_eta_1btag"] = new TH1F("jet2_eta_1btag"+hname,"2nd jet #eta",50, -2.4, 2.4);
299 algomez 1.2 hjets["1st_eta_fin"] = new TH1F("jet1_eta_fin"+hname,"Leading jet #eta",50, -2.4, 2.4);
300     hjets["2nd_eta_fin"] = new TH1F("jet2_eta_fin"+hname,"2nd jet #eta",50, -2.4, 2.4);
301 algomez 1.1 hjets["phi"] = new TH1F("jet_phi"+hname,"jet #phi",50, 0, 3.15);
302 algomez 1.3 hjets["Njets"] = new TH1F("Njets"+hname,"jet multiplicity",13,-0.5,12.5);
303 algomez 1.6 hjets["Njets_cut0"] = new TH1F("Njets_cut0"+hname,"jet multiplicity",13,-0.5,12.5);
304 algomez 1.4 hjets["Njets_cut1"] = new TH1F("Njets_cut1"+hname,"jet multiplicity",13,-0.5,12.5);
305 algomez 1.10 hjets["Njets_1btag"] = new TH1F("Njets_1tag"+hname,"jet multiplicity",13,-0.5,12.5);
306     hjets["Njets_2btag"] = new TH1F("Njets_2tag"+hname,"jet multiplicity",6,0.5,6.5);
307     hjets["Nbtags_TCHPM"] = new TH1F("Nbjets_TCHPM"+hname,"Tagged b-jets",3,-0.5,2.5);
308 algomez 1.11 hjets["Nbtags_CSVL"] = new TH1F("Nbjets_CSVL"+hname,"Tagged b-jets",8,-0.5,7.5);
309 algomez 1.10 hjets["Nbtags_CSVM"] = new TH1F("Nbjets_CSVM"+hname,"Tagged b-jets",8,-0.5,7.5);
310 algomez 1.11 hjets["Nbtags_CSVT"] = new TH1F("Nbjets_CSVT"+hname,"Tagged b-jets",8,-0.5,7.5);
311 algomez 1.1 hjets["Dijet_deltaPhi"] = new TH1F("Dijet_deltaPhi"+hname,"#Delta #phi(j,j)",30,-3.15,3.15);
312 algomez 1.6 hjets["Dijet_deltaR_cut0"] = new TH1F("Dijet_deltaR_cut0"+hname,"#DeltaR(j,j)",80,0.,4.);
313     hjets["Dijet_deltaR"] = new TH1F("Dijet_deltaR"+hname,"#DeltaR(j,j)",80,0.,4.);
314 algomez 1.10 hjets["Dijet_deltaR_1btag"] = new TH1F("Dijet_deltaR_1btag"+hname,"#DeltaR(j,j)",80,0.,4.);
315 algomez 1.1 hjets["tb_deltaPhi"] = new TH1F("tb_deltaPhi"+hname,"#Delta #phi(t,b)",30,0.,3.15);
316 algomez 1.2 hjets["tb_deltaEta"] = new TH1F("tb_deltaEta"+hname,"#Delta #eta(t,b)",30,-5,5);
317     hjets["pt_top"] = new TH1F("pt_top"+hname,"top p_{T} [GeV]",60,0,1500);
318     hjets["pt_b"] = new TH1F("pt_b"+hname,"b-jet p_{T} [GeV]",60,0,1500);
319     hjets["gen_deltaR_mub"] = new TH1F("gen_deltaR_mub"+hname,"#Delta R(#mu,b)",40,0,4);
320 algomez 1.8 hjets["1st_bdisc"] = new TH1F("1st_bdisc"+hname,"b discriminator Leading jet", 50, -0.01, 1.01);
321 algomez 1.10 hjets["1st_bdisc_1btag"] = new TH1F("1st_bdisc_1btag"+hname,"b discriminator Leading jet", 50, -0.01, 1.01);
322 algomez 1.8 hjets["2nd_bdisc"] = new TH1F("2nd_bdisc"+hname,"b discriminator 2nd jet", 50, -0.01, 1.01);
323 algomez 1.10 hjets["2nd_bdisc_1btag"] = new TH1F("2nd_bdisc_1btag"+hname,"b discriminator 2nd jet", 50, -0.01, 1.01);
324 algomez 1.8 hjets["3rd_bdisc"] = new TH1F("3rd_bdisc"+hname,"b discriminator 3rd jet", 50, -0.01, 1.01);
325 algomez 1.10 hjets["3rd_bdisc_1btag"] = new TH1F("3rd_bdisc_1btag"+hname,"b discriminator 3rd jet", 50, -0.01, 1.01);
326 algomez 1.8 hjets["4th_bdisc"] = new TH1F("4th_bdisc"+hname,"b discriminator 4th jet", 50, -0.01, 1.01);
327 algomez 1.10 hjets["4th_bdisc_1btag"] = new TH1F("4th_bdisc_1btag"+hname,"b discriminator 4th jet", 50, -0.01, 1.01);
328 algomez 1.8 hjets["5th_bdisc"] = new TH1F("5th_bdisc"+hname,"b discriminator 5th jet", 50, -0.01, 1.01);
329     hjets["6th_bdisc"] = new TH1F("6th_bdisc"+hname,"b discriminator 6th jet", 50, -0.01, 1.01);
330     hjets["7th_bdisc"] = new TH1F("7th_bdisc"+hname,"b discriminator 7th jet", 50, -0.01, 1.01);
331 algomez 1.10 hjets["prod_bdisc"] = new TH1F("prod_bdisc"+hname,"discriminators product", 50, -0.01, 1.01);
332     hjets["prod_bdisc_1btag"] = new TH1F("prod_bdisc_1btag"+hname,"discriminators product", 50, -0.01, 1.01);
333 algomez 1.8
334    
335 algomez 1.1
336     map<string,TH1* > allhistos = hmuons;
337     allhistos.insert( helectrons.begin(), helectrons.end() );
338     allhistos.insert( hMET.begin(), hMET.end() );
339     allhistos.insert( hM.begin(), hM.end() );
340     allhistos.insert( hjets.begin(), hjets.end() );
341    
342     for ( map<string,TH1* >::const_iterator imap=allhistos.begin(); imap!=allhistos.end(); ++imap )
343     {
344     TH1 *temp = imap->second;
345     temp->Sumw2();
346     temp->SetXTitle( temp->GetTitle() );
347     }
348    
349     // cut flow
350     if (fChannel==1)
351     { //muon +jets
352     fCutLabels.push_back("Processed");
353     fCutLabels.push_back("OneIsoMu");
354     fCutLabels.push_back("LooseMuVeto");
355     fCutLabels.push_back("ElectronVeto");
356     fCutLabels.push_back("MET");
357     fCutLabels.push_back("1Jet");
358     fCutLabels.push_back("2Jet");
359     fCutLabels.push_back("3Jet");
360     fCutLabels.push_back("4Jet");
361     fCutLabels.push_back("2Jet1b");
362     fCutLabels.push_back("2Jet2b");
363     fCutLabels.push_back("MaxJets");
364     fCutLabels.push_back("phi");
365     fCutLabels.push_back("topmass");
366     }
367     else
368     { //electron+jets
369     fCutLabels.push_back("Processed");
370     fCutLabels.push_back("OneIsoElectron");
371     fCutLabels.push_back("LooseMuVeto");
372     fCutLabels.push_back("ZVeto");
373     fCutLabels.push_back("ConversionVeto");
374     fCutLabels.push_back("MET");
375     fCutLabels.push_back("1Jet");
376     fCutLabels.push_back("2Jet");
377     fCutLabels.push_back("3Jet");
378     fCutLabels.push_back("4Jet");
379     fCutLabels.push_back("2Jet1b");
380     fCutLabels.push_back("2Jet2b");
381     }
382     hcutflow = new TH1D("cutflow","cut flow", fCutLabels.size(), 0.5, fCutLabels.size() +0.5 );
383    
384     for ( vector<string>::const_iterator ivec= fCutLabels.begin(); ivec!=fCutLabels.end(); ++ivec)
385     {
386     cutmap[ *ivec ] = 0;
387     }
388    
389     double pu_weights[35] = {
390     0.0255506,
391     0.251923,
392     0.549605,
393     0.924918,
394     1.25977,
395     1.48142,
396     1.57923,
397     1.57799,
398     1.52152,
399     1.4414,
400     1.35889,
401     1.2841,
402     1.21927,
403     1.16125,
404     1.11244,
405     1.06446,
406     1.01666,
407     0.966989,
408     0.913378,
409     0.85774,
410     0.799258,
411     0.734225,
412     0.670242,
413     0.607641,
414     0.54542,
415     0.484084,
416     0.427491,
417     0.369787,
418     0.321454,
419     0.280706,
420     0.238499,
421     0.198961,
422     0.166742,
423     0.146428,
424     0.224425
425     };
426    
427     fpu_weights_vec.assign( pu_weights, pu_weights + 35 );
428    
429     // For JEC uncertainties
430 algomez 1.8 if (fdoJECunc) fJECunc = new JetCorrectionUncertainty("/uscms/home/yumiceva/work/CMSSW_4_2_4/src/Yumiceva/TreeAnalyzer/test/GR_R_42_V19_AK5PF_Uncertainty.txt");
431    
432 algomez 1.9 // LOT for b-tagging SF
433     //bSF_table.LoadTable("/uscms/home/yumiceva/work/CMSSW_4_2_4/src/Yumiceva/TreeAnalyzer/test/Table_CSVM_beff_SF.txt");
434     lSF_table.LoadTable("/uscms/home/yumiceva/work/CMSSW_4_2_4/src/Yumiceva/TreeAnalyzer/test/Table_CSVM_lmistag_SF.txt");
435     leff_table.LoadTable("/uscms/home/yumiceva/work/CMSSW_4_2_4/src/Yumiceva/TreeAnalyzer/test/Table_CSVM_lmistag_Eff.txt");
436    
437     //leff_tablebtag efficiency file and 2D histogram
438     //TString btagefffilename = "/uscms/home/weizou/work/NtupleMaker/CMSSW_4_2_4/src/Yumiceva/TreeAnalyzer/test/BtagEff.root";
439     //btagefffile = new TFile(btagefffilename,"read");
440     //ttbar
441     //f2Dttbarbtag = (TH2D*) btagefffile->Get("ttbar_csv_btageff");
442     //f2Dttbarctag = (TH2D*) btagefffile->Get("ttbar_csv_ctageff");
443     //f2Dttbarlighttag = (TH2D*) btagefffile->Get("ttbar_csv_lighttageff");
444     //wprime
445     //f2Dwprimebtag = (TH2D*) btagefffile->Get("wprime_csv_btageff");
446     //f2Dwprimectag = (TH2D*) btagefffile->Get("wprime_csv_ctageff");
447     //f2Dwprimelighttag = (TH2D*) btagefffile->Get("wprime_csv_lighttageff");
448    
449    
450 algomez 1.8 //------- Store information in a Tree
451     MyStoreTree = new StoreTreeVariable();
452    
453     if(fChannel == 1){
454     MyStoreTree->SetElectronFalse();
455     //MyStoreTree->SetJetFalse();
456     MyStoreTree->SetVertexFalse();
457     MyStoreTree->SetTriggerFalse();
458 algomez 1.10 //MyStoreTree->SetMetFalse();
459 algomez 1.8 MyStoreTree->SetMuonFalse();
460     }
461    
462     if(fChannel == 2){
463    
464     MyStoreTree->SetElectronFalse();
465     //MyStoreTree->SetJetFalse();
466     MyStoreTree->SetVertexFalse();
467     MyStoreTree->SetTriggerFalse();
468 algomez 1.10 //MyStoreTree->SetMetFalse();
469 algomez 1.8 MyStoreTree->SetMuonFalse();
470     }
471    
472     MyStoreTree->InitialAll();
473     //Get the Store Tree
474     MyStoreTree->GetStoreTree()->SetDirectory(fFile);
475     MyStoreTree->GetStoreTree()->AutoSave();
476 algomez 1.9 //////////////////////////////////////////
477 algomez 1.1
478     }
479    
480     Bool_t Analyzer::Process(Long64_t entry)
481     {
482     // The Process() function is called for each entry in the tree (or possibly
483     // keyed object in the case of PROOF) to be processed. The entry argument
484     // specifies which entry in the currently loaded tree is to be processed.
485     // It can be passed to either Analyzer::GetEntry() or TBranch::GetEntry()
486     // to read either all or the required parts of the data. When processing
487     // keyed objects with PROOF, the object is already loaded and is available
488     // via the fObject pointer.
489     //
490     // This function should contain the "body" of the analysis. It can contain
491     // simple or elaborate selection criteria, run algorithms on the data
492     // of the event and typically fill histograms.
493     //
494     // The processing can be stopped by calling Abort().
495     //
496     // Use fStatus to set the return value of TTree::Process().
497     //
498     // The return value is currently not used.
499    
500     //TString option = GetOption();
501    
502 algomez 1.8 if ( entry % 100 == 0 )
503     cout<< "process entry " << entry << endl;
504 algomez 1.1
505     //TString sEntry = Form("%f", float(entry) );
506     // Info("Process",
507     //"entry # %s", Form("%f", float(entry) ) );
508    
509     //fChain->GetTree()->GetEntry(entry);
510     fChain->GetEntry(entry);
511    
512     //if (entry>10) return kTRUE;
513    
514     // event info
515     //cout << "run: " << ntuple->run << " lumi: " << ntuple->lumi << endl;
516    
517     // get collections
518     vector< TopVertexEvent > primaryVertices = ntuple->vertices;
519     vector< TopMuonEvent > muons = ntuple->muons;
520 algomez 1.2 vector< TopElectronEvent > electrons = ntuple->PFelectrons; // use PF electrons (gsf collection is called "electrons")
521 algomez 1.1 vector< TopJetEvent > jets = ntuple->PFjets;
522    
523 algomez 1.2 // USE PF Isolation
524     fMuSelector.UsePFIsolation(true);
525     fEleSelector.UsePFIsolation(true);
526    
527 algomez 1.1 size_t total_pvs = primaryVertices.size();
528     size_t total_muons = muons.size();
529     size_t total_electrons = electrons.size();
530     size_t total_jets = jets.size();
531    
532     float PVz = 0.;
533     TLorentzVector p4muon, p4ele, p4lepton, p4MET;
534     TLorentzVector p4Nu, p4OtherNu;
535     TLorentzVector p4QCDmuon;
536    
537     vector< TLorentzVector > p4jets;
538 algomez 1.6 vector< TLorentzVector > p4Othermuon; // leading muon
539 algomez 1.1
540 algomez 1.2 ////////////////////
541     // GENERATOR
542     ///////////////////
543 algomez 1.6 if (fIsMC && fSample.Contains("4Top") )
544 algomez 1.2 {
545     TLorentzVector p4genLepton;
546     TLorentzVector p4genNu;
547     TLorentzVector p4genb;
548     if (ntuple->gen.bLep_pt>0)
549     {
550     p4genLepton.SetPtEtaPhiE(ntuple->gen.mu_pt,ntuple->gen.mu_eta,ntuple->gen.mu_phi,ntuple->gen.mu_e);
551     p4genNu.SetPtEtaPhiE(ntuple->gen.nu_pt,ntuple->gen.nu_eta,ntuple->gen.nu_phi,ntuple->gen.nu_e);
552     p4genb.SetPtEtaPhiE(ntuple->gen.bLep_pt,ntuple->gen.bLep_eta,ntuple->gen.bLep_phi,ntuple->gen.bLep_e);
553    
554     hjets["gen_deltaR_mub"]->Fill( p4genLepton.DeltaR( p4genb ) );
555     }
556     }
557 algomez 1.1 ////////////////////////////////////
558     // PRIMARY VERTICES
559     ///////////////////////////////////
560    
561 algomez 1.6 for ( size_t ipv=0; ipv < total_pvs; ++ipv)
562     {
563     if (ipv==0) PVz = primaryVertices[ipv].vz;
564     }
565 algomez 1.1
566 algomez 1.6 //hPVs["N"]->Fill( total_pvs );
567 algomez 1.1
568     // calculate PU weight
569     double PUweight = 1.;
570    
571     if (fIsMC) {
572 algomez 1.2
573     Int_t mc_npvminus1 = TMath::Min(ntuple->gen.Bx_minus1,34);
574     Int_t mc_npv0 = TMath::Min(ntuple->gen.Bx_0,34);
575     Int_t mc_npvplus1 = TMath::Min(ntuple->gen.Bx_plus1,34);
576    
577     PUweight = f3Dweight->GetBinContent(mc_npvminus1,mc_npv0,mc_npvplus1);
578    
579     //int iibin = 0;
580     //for ( vector<double>::iterator ivec = fpu_weights_vec.begin(); ivec != fpu_weights_vec.end(); ++ivec )
581     // {
582     // int mc_npvs = ntuple->gen.Bx_0; // in-time pile up
583     ////int mc_npvs = (int)total_pvs;
584     //if ( mc_npvs >= iibin+1 ) PUweight = *ivec; // use the last weight for last bin
585     //if ( ( iibin <= mc_npvs ) && ( mc_npvs < iibin + 1 ) ) PUweight = *ivec;
586     //iibin++;
587     //}
588 algomez 1.1 }
589    
590 algomez 1.5 // For 4 tops
591     if (fIsMC && fSample.Contains("4Top") ) {
592     PUweight = 1;
593     }
594    
595 algomez 1.1 hPVs["Nreweight"]->Fill( total_pvs, PUweight );
596    
597 algomez 1.2 /////////////
598     // HLT scale factor for MC
599     ////////////
600     double SF_hlt = 1.;
601 algomez 1.9 if (fIsMC) SF_hlt = 0.97; //0.966;
602     double SF_iso = 1.;
603     if (fIsMC) SF_iso = 0.996;
604     PUweight = PUweight*SF_hlt*SF_iso; // LETS INCLUDE THE TRIGGER SF INTO THE PU WEIGHTS
605    
606     //TEMP
607     //PUweight = 1.;
608 algomez 1.2
609 algomez 1.1 cutmap["Processed"] += PUweight;
610    
611     if (fVerbose) cout << "done pv" << endl;
612 algomez 1.6
613 algomez 1.1 //////////////////////////////////
614     // MUONS
615     //////////////////////////////////
616 algomez 1.2 int ngoodIDmuons = 0;
617 algomez 1.1 int nloosemuons = 0;
618     int ntightmuons = 0;
619     int nqcdmuons = 0;
620 algomez 1.9 int nmuonID[8] = {0};
621 algomez 1.1
622     double RelIso = -1.;
623     double deltaR = -1.;
624     double QCDRelIso = -1.;
625     double QCDdeltaR = -1;
626    
627 algomez 1.2
628 algomez 1.1 for ( size_t imu=0; imu < total_muons; ++imu) {
629    
630 algomez 1.6 TopMuonEvent muon = muons[imu];
631 algomez 1.9
632 algomez 1.6 h1test->Fill( muon.pt );
633     hmuons["pt_cut1"]->Fill( muon.pt, PUweight );
634 algomez 1.9
635     if (muon.pt > 32. && fabs(muon.eta)<2.1 && muon.IsTrackerMuon==1 ) {
636     nmuonID[0] += 1;
637     if (fabs(muon.d0)<0.02 ) {
638     nmuonID[1] += 1;
639     if (muon.trackerhits>=11 ) {
640     nmuonID[2] += 1;
641     if (muon.normchi2<10 ) {
642     nmuonID[3] +=1;
643     if (muon.muonhits>0 ) {
644     nmuonID[4] += 1;
645     if (muon.pixelhits >= 1 ) {
646     nmuonID[5] += 1;
647     if (muon.muonstations > 1 ) {
648     nmuonID[6] +=1;
649     if (fabs(muon.vz - PVz) < 1.) {
650     nmuonID[7] +=1;
651     }
652     }
653     }
654     }
655     }
656     }
657     }
658     }
659 algomez 1.6
660     if ( fMuSelector.MuonID( muon, PVz ) ) ngoodIDmuons++;
661 algomez 1.9
662 algomez 1.6 // select only good muons
663 algomez 1.9
664 algomez 1.6 if ( fMuSelector.MuonLoose( muon ) ) {
665 algomez 1.9
666 algomez 1.6 nloosemuons++;
667 algomez 1.9
668     if ( fMuSelector.MuonTight( muon, PVz) ) hmuons["pt_cut2"]->Fill( muon.pt, PUweight );
669 algomez 1.6 //hmuons["N_tisomuons"]->Fill( nloosemuons );
670     //hmuons["charge_tiso"]->Fill( muon.charge, PUweight );
671     if ( fMuSelector.MuonTightDeltaR( muon, PVz, jets) ) {
672     ntightmuons++;
673     deltaR = fMuSelector.GetDeltaR();
674     }
675 algomez 1.1
676 algomez 1.6 p4muon.SetPtEtaPhiE( muon.pt, muon.eta, muon.phi, muon.e );
677     RelIso = muon.pfreliso; //muon.reliso03;
678     p4Othermuon.push_back( p4muon ); // for leading muon
679 algomez 1.9 }
680 algomez 1.4
681 algomez 1.6 // check muon in QCD control region
682     if ( fMuSelector.MuonRelax02IsoQCD( muon, PVz, jets ) ) {
683     nqcdmuons++;
684 algomez 1.9
685 algomez 1.6 // keep the leading muon for selection
686     if (nqcdmuons==1) {
687     p4QCDmuon.SetPtEtaPhiE( muon.pt, muon.eta, muon.phi, muon.e );
688     QCDRelIso = muon.pfreliso; //muon.reliso03;
689     QCDdeltaR = fMuSelector.GetDeltaR();
690     }
691 algomez 1.1 }
692 algomez 1.6 }
693 algomez 1.1
694 algomez 1.6 if (fVerbose) cout << "done muons" << endl;
695 algomez 1.1
696     //////////////////////////////////
697     // ELECTRONS
698     //////////////////////////////////
699     int nlooseelectrons = 0;
700     int ntightelectrons = 0;
701     bool IsConversion = false;
702    
703     for ( size_t iele=0; iele < total_electrons; ++iele) {
704    
705 algomez 1.6 TopElectronEvent electron = electrons[iele];
706 algomez 1.1
707 algomez 1.6 if ( fEleSelector.ElectronLoose(electron) ) nlooseelectrons++;
708 algomez 1.1
709 algomez 1.6 if ( fEleSelector.ElectronTight(electron, PVz ) ) {
710 algomez 1.1
711 algomez 1.6 if (ntightelectrons == 0) {
712     IsConversion = electron.IsConversion;
713     p4ele.SetPtEtaPhiE( electron.pt, electron.eta, electron.phi, electron.e );
714     helectrons["pt_cut2"]->Fill( p4ele.Pt(), PUweight );
715     helectrons["eta_cut2"]->Fill( p4ele.Eta(), PUweight );
716     helectrons["phi_cut2"]->Fill( p4ele.Phi(), PUweight );
717     }
718     ntightelectrons++;
719     }
720 algomez 1.1 }
721     if (fVerbose) cout << "done electron" << endl;
722    
723 algomez 1.6 /////////////////////////////////////
724     // MUON/ELECTRON + JETS
725     /////////////////////////////////////
726 algomez 1.1
727 algomez 1.6 if ( fChannel == 1 ) {
728 algomez 1.1
729 algomez 1.6 if (fdoQCD2SideBand) {
730 algomez 1.1
731 algomez 1.6 if (nqcdmuons == 0) return kTRUE;
732     cutmap["OneIsoMu"] += PUweight;
733 algomez 1.1
734 algomez 1.6 p4lepton = p4QCDmuon;
735     RelIso = QCDRelIso;
736     deltaR = QCDdeltaR;
737     }
738     else {
739 algomez 1.9 /*
740     if (nmuonID[0] == 1) cutmap["MuID1"] += PUweight;
741     if (nmuonID[1] == 1) cutmap["MuID2"] +=PUweight;
742     if (nmuonID[2] == 1) cutmap["MuID3"] +=PUweight;
743     if (nmuonID[3] == 1) cutmap["MuID4"] +=PUweight;
744     if (nmuonID[4] == 1) cutmap["MuID5"] +=PUweight;
745     if (nmuonID[5] == 1) cutmap["MuID6"] +=PUweight;
746     if (nmuonID[6] == 1) cutmap["MuID7"] +=PUweight;
747     if (nmuonID[7] == 1) cutmap["MuID8"] +=PUweight;
748     if (ngoodIDmuons == 1) cutmap["MuID"] += PUweight;*/ // TEMP
749 algomez 1.2
750 algomez 1.6 if ( ngoodIDmuons > 0 ) hmuons["Ngood"]->Fill( total_pvs, PUweight);
751     if ( ntightmuons > 0 ) hmuons["Niso"]->Fill( total_pvs, PUweight);
752     if ( ntightmuons != 1 ) return kTRUE;
753     cutmap["OneIsoMu"] += PUweight;
754 algomez 1.2
755 algomez 1.6 if ( nloosemuons > 1 ) return kTRUE;
756     cutmap["LooseMuVeto"] += PUweight;
757 algomez 1.1
758 algomez 1.6 if ( nlooseelectrons > 0 ) return kTRUE;
759     cutmap["ElectronVeto"] += PUweight;
760 algomez 1.1
761 algomez 1.6 p4lepton = p4muon;
762     if (fVerbose) cout << "got a good lepton" << endl;
763     }
764    
765     hmuons["pt"]->Fill( p4lepton.Pt(), PUweight );
766     hmuons["eta"]->Fill( p4lepton.Eta(), PUweight );
767     hmuons["phi"]->Fill( p4lepton.Phi(), PUweight );
768     hmuons["reliso"]->Fill( RelIso, PUweight );
769     //hmuons["deltaR"]->Fill( deltaR, PUweight );
770 algomez 1.1
771 algomez 1.6 }
772 algomez 1.1 else // electron+jets
773     {
774     // pending ...
775     }
776    
777     if (fVerbose) cout << "done lepton selection " << endl;
778    
779     /////////////////////////////////
780 algomez 1.6 // MET + Ht
781 algomez 1.1 /////////////////////////////////
782    
783     p4MET.SetPtEtaPhiE( ntuple->PFMET,
784     0.,
785     ntuple->PFMETphi,
786     ntuple->PFMET );
787    
788     //temporal check using genMET
789     //p4MET.SetPtEtaPhiE( ntuple->gen.MET,
790     // 0.,
791     // ntuple->gen.METPhi,
792     // ntuple->gen.MET );
793    
794     if (fdoQCD1SideBand && p4MET.Et() > 20.) return kTRUE;
795     else if ( p4MET.Et() <= 20. && fdoQCD2SideBand==false ) return kTRUE;
796 algomez 1.6 if (fVerbose) cout << "pass MET cut" << endl;
797 algomez 1.1
798    
799     //cutmap["MET"] += PUweight;
800 algomez 1.6 //hMET["MET"]->Fill( p4MET.Pt(), PUweight );
801     //hMET["phi"]->Fill( p4MET.Phi(), PUweight );
802     //hMET["Ht0"]->Fill( ntuple->PFHt, PUweight );
803     //hMET["Htlep"]->Fill( ntuple->PFHt + p4lepton.Pt(), PUweight );
804    
805 algomez 1.1
806     double Wpt = p4lepton.Pt() + p4MET.Pt();
807     double Wpx = p4lepton.Px() + p4MET.Px();
808     double Wpy = p4lepton.Py() + p4MET.Py();
809     double WMt = sqrt(Wpt*Wpt-Wpx*Wpx-Wpy*Wpy);
810 algomez 1.3
811     //if ( WMt <= 40. ) return kTRUE;
812 algomez 1.1 cutmap["MET"] += PUweight;
813    
814     if (fVerbose) cout << "pass W MT cut " << endl;
815    
816     /////////////////////////////////
817     // estimate Pz of neutrino
818     ////////////////////////////////
819    
820     fzCalculator.SetMET(p4MET);
821     fzCalculator.SetLepton(p4lepton);
822 algomez 1.6 if (fChannel==2) {
823     fzCalculator.SetLeptonType("electron");
824     }
825 algomez 1.9
826     double pzNu = fzCalculator.Calculate();
827     double pzOtherNu = fzCalculator.getOther();
828     if (fabs(pzNu) > fabs(pzOtherNu) ) {
829     double tmppzNu = pzNu;
830     double tmppzOtherNu = pzOtherNu;
831     pzNu = tmppzOtherNu;
832     pzOtherNu = tmppzNu;
833     }
834 algomez 1.1
835     p4Nu = TLorentzVector();
836     p4OtherNu = TLorentzVector();
837    
838     p4Nu.SetPxPyPzE(p4MET.Px(), p4MET.Py(), pzNu, sqrt(p4MET.Px()*p4MET.Px()+p4MET.Py()*p4MET.Py()+pzNu*pzNu));
839     //print "pzNu = " +str(pzNu)
840     //print "p4Nu =("+str(p4Nu.Px())+","+str(p4Nu.Py())+","+str(p4Nu.Pz())+","+str(p4Nu.E())
841     p4OtherNu.SetPxPyPzE( p4MET.Px(), p4MET.Py(),pzOtherNu,sqrt(p4MET.Px()*p4MET.Px()+p4MET.Py()*p4MET.Py()+pzOtherNu*pzOtherNu));
842    
843 algomez 1.6 //double WmassNoPt = (p4Nu+p4lepton).M();
844 algomez 1.1
845 algomez 1.6 if ( fzCalculator.IsComplex() ) {
846     double ptNu1 = fzCalculator.getPtneutrino(1);
847     double ptNu2 = fzCalculator.getPtneutrino(2);
848     TLorentzVector p4Nu1tmp;
849     TLorentzVector p4Nu2tmp;
850 algomez 1.1
851 algomez 1.6 p4Nu1tmp.SetPxPyPzE( ptNu1*p4MET.Px()/p4MET.Pt(), ptNu1*p4MET.Py()/p4MET.Pt(), pzNu, sqrt(ptNu1*ptNu1+pzNu*pzNu));
852     p4Nu2tmp.SetPxPyPzE( ptNu2*p4MET.Px()/p4MET.Pt(), ptNu2*p4MET.Py()/p4MET.Pt(), pzNu, sqrt(ptNu2*ptNu2+pzNu*pzNu));
853 algomez 1.1
854 algomez 1.6 TLorentzVector Wtmp;
855     Wtmp = p4lepton + p4Nu1tmp;
856     double Wm1 = 0;
857     double Wm2 = 0;
858     Wm1 = Wtmp.M();
859     Wtmp = p4lepton + p4Nu2tmp;
860     Wm2 = Wtmp.M();
861     if ( fabs( Wm1 - 80.) < fabs( Wm2 - 80.) ) p4Nu = p4Nu1tmp;
862     else p4Nu = p4Nu2tmp;
863 algomez 1.1
864 algomez 1.6 p4OtherNu = p4Nu; // since we chose the real part, the two solutions are the same.
865     }
866 algomez 1.1
867    
868     hMET["PzNu"]->Fill(pzNu, PUweight ); //change this to 2d with two sol and as a function of jets
869    
870     TLorentzVector p4LepW = p4lepton + p4Nu;
871     TLorentzVector p4OtherLepW = p4lepton + p4OtherNu;
872    
873     //hMET["LepWmass"]->Fill(p4LepW.M(), PUweight );
874     //if ( fzCalculator.IsComplex() ) hMET["LepWmassComplex"]->Fill( p4LepW.M(), PUweight );
875    
876    
877     /////////////////////////////////
878     // JETS
879     ////////////////////////////////
880    
881 algomez 1.6 //JetCombinatorics myCombi = JetCombinatoric();
882 algomez 1.1
883     int njets = 0;
884 algomez 1.11 //double prod_bdisc =0;
885 algomez 1.8 MyStoreTree->GetJetVariable()->numjets = 0;
886 algomez 1.11 MyStoreTree->GetJetVariable()->numBjets_csvl = 0;
887     MyStoreTree->GetJetVariable()->numBjets_csvm = 0;
888     MyStoreTree->GetJetVariable()->numBjets_csvt = 0;
889 algomez 1.1 map< string, vector<float> > bdisc;
890     map< string, vector<bool> > isTagb;
891 algomez 1.9 map< string, vector<bool> > isTagbUp;
892     map< string, vector<bool> > isTagbDown;
893 algomez 1.1 vector<int> listflavor;
894 algomez 1.8 vector<float> bdiscriminator;
895    
896 algomez 1.1
897 algomez 1.6 for ( size_t ijet=0; ijet < total_jets; ++ijet) {
898 algomez 1.1
899 algomez 1.6 TopJetEvent jet = jets[ijet];
900     double SF_JEC = 1.;
901     if (fdoJECunc){
902     fJECunc->setJetEta( jet.eta);
903     fJECunc->setJetPt( jet.pt);
904     double jec_unc = fJECunc->getUncertainty(true);
905     if (fVerbose) cout << "JEC uncertainty is " << jec_unc << endl;
906     if (fdoJECup) SF_JEC = 1.+jec_unc;
907     else SF_JEC = 1.-jec_unc;
908 algomez 1.1 }
909    
910 algomez 1.7 if ( SF_JEC*jet.pt > 40. && fabs(jet.eta) < 2.4 && SF_JEC*jets[0].pt > 100. && SF_JEC*jets[1].pt > 60. ) { //first cut
911 algomez 1.6
912 algomez 1.8 //if (fVerbose) cout << " jet pt " << SF_JEC*jet.pt << endl;
913 algomez 1.6
914     //hjets["pt"]->Fill( jet.pt, PUweight );
915     //hjets["eta"]->Fill( jet.eta, PUweight );
916     //hjets["phi"]->Fill( jet.phi, PUweight );
917    
918     TLorentzVector tmpjet;
919     tmpjet.SetPtEtaPhiE(SF_JEC*jet.pt, jet.eta, jet.phi, SF_JEC*jet.e);
920     p4jets.push_back( tmpjet);
921     listflavor.push_back( jet.mc.flavor );
922 algomez 1.9 if (jet.btag_CSV > 0) bdiscriminator.push_back( jet.btag_CSV ); // for bdiscriminator
923 algomez 1.8 else bdiscriminator.push_back( 0 );
924 algomez 1.6
925     if (fVerbose) {
926     cout << "done storing njets " << njets << endl;
927 algomez 1.8 //cout << " bdisc " << jet.btag_TCHP << endl;
928 algomez 1.6 cout << " bdisc " << jet.btag_CSV << endl;
929     }
930    
931 algomez 1.7 // store discriminators
932 algomez 1.8 //bdisc["TCHP"].push_back( jet.btag_TCHP );
933 algomez 1.6 bdisc["CSV"].push_back( jet.btag_CSV );
934     if (fVerbose) cout << "store bdisc" << endl;
935     // TCHPL cut at 1.19
936     // check TCHPM cut at 1.93
937 algomez 1.7 //if ( jet.btag_TCHP > 1.93 ) isTagb["TCHPM"].push_back(true);
938     //else isTagb["TCHPM"].push_back(false);
939     //if (fVerbose) cout << "done tchpl" << endl;
940 algomez 1.6 // check SSVHEM cut at 1.74
941 algomez 1.11
942     if ( jet.btag_CSV > 0.244) {
943     isTagb["CSVL"].push_back(true);
944     //isTagbUp["CSVL"].push_back(true);
945     //isTagbDown["CSVL"].push_back(true);
946     } else {
947     isTagb["CSVL"].push_back(false);
948     //isTagbUp["CSVL"].push_back(false);
949     //isTagbDown["CSVL"].push_back(false);
950     }
951     // CSVM cut at 0.244 LOOSE
952 algomez 1.9 if ( jet.btag_CSV > 0.679) {
953     isTagb["CSVM"].push_back(true);
954     isTagbUp["CSVM"].push_back(true);
955     isTagbDown["CSVM"].push_back(true);
956     } else {
957     isTagb["CSVM"].push_back(false);
958     isTagbUp["CSVM"].push_back(false);
959     isTagbDown["CSVM"].push_back(false);
960     }
961 algomez 1.7 // CSVM cut at 0.679 MEDIUM
962 algomez 1.11
963     if ( jet.btag_CSV > 0.898) {
964     isTagb["CSVT"].push_back(true);
965     //isTagbUp["CSVT"].push_back(true);
966     //isTagbDown["CSVT"].push_back(true);
967     } else {
968     isTagb["CSVT"].push_back(false);
969     //isTagbUp["CSVT"].push_back(false);
970     //isTagbDown["CSVT"].push_back(false);
971     }
972     // CSVM cut at 0.898 TIGHT
973    
974     if (fVerbose) cout << "done csv" << endl;
975 algomez 1.7
976 algomez 1.6 njets++;
977     }
978     }
979 algomez 1.1
980 algomez 1.10 // if (njets>0) {
981     // /////////// plots without cuts
982     // hM["WMt"]->Fill( WMt, PUweight );
983     // hjets["Njets"]->Fill( njets, PUweight );
984     // hMET["Ht0"]->Fill( ntuple->PFHt, PUweight );
985     // hjets["Njets_cut0"]->Fill( njets, PUweight );
986     // hjets["1st_pt_cut0"]->Fill( p4jets[0].Pt(), PUweight );
987     // hjets["1st_eta_cut0"]->Fill( p4jets[0].Eta(), PUweight );
988     // hjets["2nd_pt_cut0"]->Fill( p4jets[1].Pt(), PUweight );
989     // hjets["2nd_eta_cut0"]->Fill( p4jets[1].Eta(), PUweight );
990     // hjets["3rd_pt_cut0"]->Fill( p4jets[2].Pt(), PUweight );
991     // hjets["4th_pt_cut0"]->Fill( p4jets[3].Pt(), PUweight );
992     // hjets["5th_pt_cut0"]->Fill( p4jets[4].Pt(), PUweight );
993     // hjets["6th_pt_cut0"]->Fill( p4jets[5].Pt(), PUweight );
994     // hjets["7th_pt_cut0"]->Fill( p4jets[6].Pt(), PUweight );
995     // hmuons["N_cut0"]->Fill( total_muons, PUweight );
996     // hmuons["Nelectrons_cut0"]->Fill( nlooseelectrons, PUweight );
997     // hMET["MET_cut0"]->Fill( p4MET.Pt(), PUweight );
998     // hmuons["deltaR_cut0"]->Fill( deltaR, PUweight );
999     // //hjets["Dijet_deltaR_cut0"]->Fill( deltaRjj, PUweight );
1000     // /////////////////////////////////////////////////////////////////////
1001     // }
1002 algomez 1.1
1003     if (fVerbose) cout << "done jets" << endl;
1004    
1005     if (njets > 0 ) cutmap["1Jet"] += PUweight;
1006     if (njets > 1 ) cutmap["2Jet"] += PUweight;
1007     if (njets > 2 ) cutmap["3Jet"] += PUweight;
1008     if (njets > 3 ) cutmap["4Jet"] += PUweight;
1009    
1010 algomez 1.10 if (njets >= 4) {
1011 algomez 1.1
1012 algomez 1.7 // count partons
1013     int number_of_b = 0;
1014     int number_of_c = 0;
1015     int number_of_l = 0;
1016    
1017 algomez 1.6 double Ht = 0;
1018 algomez 1.8 double Stlep = 0;
1019     double Stjet = 0;
1020 algomez 1.6 double deltaRjj = 999.;
1021 algomez 1.7
1022 algomez 1.8
1023 algomez 1.6 for ( size_t kk=0; kk < p4jets.size(); ++kk) {
1024 algomez 1.7 // Ht calculation
1025 algomez 1.6 Ht += p4jets[kk].Pt();
1026    
1027 algomez 1.8 // St Calculation
1028     Stlep = p4MET.Pt() + p4lepton.Pt() ;
1029     Stjet = p4MET.Pt() + p4lepton.Pt() + Ht ;
1030    
1031    
1032 algomez 1.7 // deltaR(jet,jet)
1033 algomez 1.6 for ( vector< TopJetEvent>::iterator ijet=jets.begin(); ijet != jets.end(); ++ijet) {
1034    
1035     TopJetEvent jet = *ijet;
1036     TLorentzVector tmpp4Jet;
1037     tmpp4Jet.SetPtEtaPhiE(jet.pt, jet.eta, jet.phi, jet.e );
1038     double tmpdeltaR = p4jets[kk].DeltaR(tmpp4Jet);
1039     if ( tmpdeltaR < 0.3 ) continue;
1040     if ( tmpdeltaR < deltaRjj ) deltaRjj = tmpdeltaR;
1041     }
1042 algomez 1.7
1043     // b-tag
1044     if ( abs(listflavor[kk])==5 && p4jets[kk].Pt()<=240 ) { number_of_b++; hjets["pt_b_mc"]->Fill( p4jets[kk].Pt(), PUweight );}
1045     if ( abs(listflavor[kk])==4 && p4jets[kk].Pt()<=240 ) { number_of_c++; hjets["pt_c_mc"]->Fill( p4jets[kk].Pt(), PUweight );}
1046     if ( abs(listflavor[kk])==1 || abs(listflavor[kk])==2 || abs(listflavor[kk])==3 || abs(listflavor[kk])==21 ) { number_of_l++; hjets["pt_l_mc"]->Fill( p4jets[kk].Pt(), PUweight );}
1047 algomez 1.9 //if ( abs(listflavor[kk])==5 && p4jets[kk].Pt()>240 ) { number_of_b_highpt++; hjets["pt_b_mc"]->Fill( p4jets[kk].Pt(), PUweight );}
1048     //if ( abs(listflavor[kk])==4 && p4jets[kk].Pt()>240 ) { number_of_c_highpt++; hjets["pt_c_mc"]->Fill( p4jets[kk].Pt(), PUweight );}
1049 algomez 1.8
1050 algomez 1.9
1051     hjets["pt"]->Fill( p4jets[kk].Pt(), PUweight );
1052     //float b_mc_eff = f2Dttbarbtag->GetBinContent(f2Dttbarbtag->FindBin(p4jets[kk].Pt(),fabs(p4jets[kk].Eta())));
1053     //float c_mc_eff = f2Dttbarctag->GetBinContent(f2Dttbarctag->FindBin(p4jets[kk].Pt(),fabs(p4jets[kk].Eta())));
1054     //float l_mc_eff = f2Dttbarlighttag->GetBinContent(f2Dttbarlighttag->FindBin(p4jets[kk].Pt(),fabs(p4jets[kk].Eta())));
1055     float b_data_eff = 0.705;
1056     float c_data_eff = b_data_eff/5.;
1057    
1058     float b_SF = 0.972; //bSF_table.GetValue( p4jets[kk].Pt(), fabs(p4jets[kk].Eta()) );
1059     float l_SF = lSF_table.GetValue( p4jets[kk].Pt(), fabs(p4jets[kk].Eta()) );
1060     float l_data_eff = leff_table.GetValue( p4jets[kk].Pt(), fabs(p4jets[kk].Eta()) );
1061    
1062     double seed = abs(static_cast<int>(sin(p4jets[kk].Phi()*1000000)*100000));
1063     BTagSFUtil btsfutil = BTagSFUtil( seed );
1064    
1065     if ( fIsMC && abs(listflavor[kk])==5 ) {
1066     number_of_b++; hjets["pt_b_mc"]->Fill( p4jets[kk].Pt(), PUweight );
1067     //Info("Process",TString("btag SF= ")+TString(Form("%f",b_SF)));
1068    
1069     bool tmp_IsTag = isTagb["CSVM"][kk];
1070     isTagb["CSVM"][kk] = btsfutil.modifyBTagsWithSF( tmp_IsTag, abs(listflavor[kk]), b_SF, b_data_eff, l_SF, l_data_eff);
1071    
1072     if (p4jets[kk].Pt()<=240 ) {
1073     tmp_IsTag = isTagbUp["CSVM"][kk];
1074     isTagbUp["CSVM"][kk] = btsfutil.modifyBTagsWithSF( tmp_IsTag , abs(listflavor[kk]), b_SF*1.04, b_data_eff, l_SF, l_data_eff);
1075     tmp_IsTag = isTagbDown["CSVM"][kk];
1076     isTagbDown["CSVM"][kk] = btsfutil.modifyBTagsWithSF( tmp_IsTag, abs(listflavor[kk]), b_SF*0.96, b_data_eff, l_SF, l_data_eff);
1077     } else {
1078     tmp_IsTag = isTagbUp["CSVM"][kk];
1079     isTagbUp["CSVM"][kk] = btsfutil.modifyBTagsWithSF( tmp_IsTag , abs(listflavor[kk]), b_SF*1.15, b_data_eff, l_SF, l_data_eff);
1080     tmp_IsTag = isTagbDown["CSVM"][kk];
1081     isTagbDown["CSVM"][kk] = btsfutil.modifyBTagsWithSF( tmp_IsTag, abs(listflavor[kk]), b_SF*0.85, b_data_eff, l_SF, l_data_eff);
1082     }
1083     }
1084    
1085     if ( fIsMC && abs(listflavor[kk])==4 ) {
1086     number_of_c++; hjets["pt_c_mc"]->Fill( p4jets[kk].Pt(), PUweight );
1087     bool tmp_IsTag = isTagb["CSVM"][kk];
1088     isTagb["CSVM"][kk] = btsfutil.modifyBTagsWithSF( tmp_IsTag, abs(listflavor[kk]), b_SF, c_data_eff, l_SF, l_data_eff);
1089     if (p4jets[kk].Pt()<=240 ) {
1090     tmp_IsTag = isTagbUp["CSVM"][kk];
1091     isTagbUp["CSVM"][kk] = btsfutil.modifyBTagsWithSF( tmp_IsTag , abs(listflavor[kk]), b_SF*1.08, c_data_eff, l_SF, l_data_eff);
1092     tmp_IsTag = isTagbDown["CSVM"][kk];
1093     isTagbDown["CSVM"][kk] = btsfutil.modifyBTagsWithSF( tmp_IsTag, abs(listflavor[kk]), b_SF*0.92, c_data_eff, l_SF, l_data_eff);
1094     } else {
1095     tmp_IsTag = isTagbUp["CSVM"][kk];
1096     isTagbUp["CSVM"][kk] = btsfutil.modifyBTagsWithSF( tmp_IsTag , abs(listflavor[kk]), b_SF*1.15, c_data_eff, l_SF, l_data_eff);
1097     tmp_IsTag = isTagbDown["CSVM"][kk];
1098     isTagbDown["CSVM"][kk] = btsfutil.modifyBTagsWithSF( tmp_IsTag, abs(listflavor[kk]), b_SF*0.85, c_data_eff, l_SF, l_data_eff);
1099     }
1100     }
1101    
1102     if ( fIsMC && (abs(listflavor[kk])==1 || abs(listflavor[kk])==2 || abs(listflavor[kk])==3 || abs(listflavor[kk])==21 )) {
1103     number_of_l++; hjets["pt_l_mc"]->Fill( p4jets[kk].Pt(), PUweight );
1104     bool tmp_IsTag = isTagb["CSVM"][kk];
1105     isTagb["CSVM"][kk] = btsfutil.modifyBTagsWithSF( tmp_IsTag, abs(listflavor[kk]), b_SF, b_data_eff, l_SF, l_data_eff);
1106     if (p4jets[kk].Pt()<=240 ) {
1107     tmp_IsTag = isTagbUp["CSVM"][kk];
1108     isTagbUp["CSVM"][kk] = btsfutil.modifyBTagsWithSF( tmp_IsTag , abs(listflavor[kk]), b_SF, b_data_eff, l_SF*1.13, l_data_eff);
1109     tmp_IsTag = isTagbDown["CSVM"][kk];
1110     isTagbDown["CSVM"][kk] = btsfutil.modifyBTagsWithSF( tmp_IsTag, abs(listflavor[kk]), b_SF, b_data_eff, l_SF*1.13, l_data_eff);
1111     } else {
1112     tmp_IsTag = isTagbUp["CSVM"][kk];
1113     isTagbUp["CSVM"][kk] = btsfutil.modifyBTagsWithSF( tmp_IsTag , abs(listflavor[kk]), b_SF, b_data_eff, l_SF*1.20, l_data_eff);
1114     tmp_IsTag = isTagbDown["CSVM"][kk];
1115     isTagbDown["CSVM"][kk] = btsfutil.modifyBTagsWithSF( tmp_IsTag, abs(listflavor[kk]), b_SF, b_data_eff, l_SF*1.20, l_data_eff);
1116     }
1117     }
1118    
1119     if ( isTagb["CSVM"][kk] ) {
1120     hjets["pt_btag"]->Fill( p4jets[kk].Pt(), PUweight );
1121     if ( abs(listflavor[kk])==5 ) hjets["pt_btag_b"]->Fill( p4jets[kk].Pt(), PUweight );
1122     if ( abs(listflavor[kk])==4 ) hjets["pt_btag_c"]->Fill( p4jets[kk].Pt(), PUweight );
1123     if ( abs(listflavor[kk])==1 || abs(listflavor[kk])==2 || abs(listflavor[kk])==3 || abs(listflavor[kk])==21 ) hjets["pt_btag_l"]->Fill( p4jets[kk].Pt(), PUweight );
1124     }
1125 algomez 1.1 }
1126    
1127 algomez 1.9 // W+jets h.f. corrections
1128     int FH = 2; // 0=Wbb, 1=Wcc, 2=Wqq
1129     if ( number_of_b > 0 ) FH = 0;
1130     else if ( number_of_b == 0 && number_of_c > 0 ) FH = 1;
1131     else FH = 2;
1132     float SF_W = 1.;
1133     if (fSample=="WJets") {
1134     // from TOP-11-003
1135     if ( FH == 0 ) SF_W = 1.21;
1136     if ( FH == 1 ) SF_W = 1.66;
1137     }
1138    
1139     // Split W+jets if requested
1140     if (fSample=="Wbb" && FH != 0 ) return kTRUE;
1141     if (fSample=="Wcc" && FH != 1 ) return kTRUE;
1142     if (fSample=="Wqq" && FH != 2 ) return kTRUE;
1143    
1144 algomez 1.7 // count number of b-tags
1145 algomez 1.11 //int Nbtags_TCHPM = 0;
1146     int Nbtags_CSVL = 0;
1147     int Nbtags_CSVT = 0;
1148 algomez 1.7 int Nbtags_CSVM = 0;
1149 algomez 1.10 int NbtagsUp_CSVM = 0;
1150     int NbtagsDown_CSVM = 0;
1151 algomez 1.11 //float SFb_0tag = 1.;
1152     //float SFb_only1tag = 1.;
1153 algomez 1.10 float SFb_1tag = 1.;
1154 algomez 1.11 //float SFb_2tag = 1.;
1155     //float SFb_0tag_syst[2] = {1.}; // for systematics
1156     //float SFb_1tag_syst[2] = {1.};
1157     //float SFb_2tag_syst[2] = {1.};
1158 algomez 1.10
1159 algomez 1.7 for ( size_t itag=0; itag< isTagb["CSVM"].size(); ++itag ) {
1160 algomez 1.11 //if ( isTagb["TCHPM"][itag] ) Nbtags_TCHPM++;
1161     if ( isTagb["CSVL"][itag] ) Nbtags_CSVL++;
1162 algomez 1.7 if ( isTagb["CSVM"][itag] ) Nbtags_CSVM++;
1163 algomez 1.10 if ( isTagbUp["CSVM"][itag] ) NbtagsUp_CSVM++;
1164     if ( isTagbDown["CSVM"][itag] ) NbtagsDown_CSVM++;
1165 algomez 1.11 if ( isTagb["CSVT"][itag] ) Nbtags_CSVT++;
1166 algomez 1.7 }
1167 algomez 1.10
1168     // store b tags
1169 algomez 1.11 //hjets["Nbtags_CSVL"]->Fill( Nbtags_CSVL, PUweight*SF_W );
1170     //hjets["Nbtags_CSVM"]->Fill( Nbtags_CSVM, PUweight*SF_W );
1171     //hjets["Nbtags_CSVT"]->Fill( Nbtags_CSVT, PUweight*SF_W );
1172 algomez 1.7
1173     // compute b-tag event weight
1174     if ( fIsMC ) {
1175 algomez 1.9 hMET["genMET_2jet"]->Fill( ntuple->gen.MET, PUweight*SF_W );
1176     hMET["deltaMET_2jet"]->Fill( p4MET.Pt() - ntuple->gen.MET, PUweight*SF_W );
1177     }
1178     /*
1179 algomez 1.7 // zeto tag
1180     BTagWeight b0(0,0); // number of tags
1181     //BTagWeight::JetInfo bj(0.63,0.91); // mean MC eff and mean SF. For TCHPM=0.91\pm0.09, CSVM=0.96\pm0.096
1182     //BTagWeight::JetInfo cj(0.15,0.91);
1183     BTagWeight::JetInfo bj(0.63,0.96);
1184     BTagWeight::JetInfo cj(0.15,0.96);
1185     double light_mceff = 0.017; //CHECK
1186     if ( 100 < p4jets[0].Pt() && p4jets[0].Pt() <= 200 ) light_mceff = 0.04;
1187     if ( 200 < p4jets[0].Pt() && p4jets[0].Pt() <= 300 ) light_mceff = 0.08;
1188     if ( 300 < p4jets[0].Pt() && p4jets[0].Pt() <= 400 ) light_mceff = 0.12;
1189     if ( 400 < p4jets[0].Pt() ) light_mceff = 0.14;
1190    
1191     //BTagWeight::JetInfo lj(light_mceff,1.22); //for TCHPM=1.22, CSVM=1.08 \pm 0.13
1192     BTagWeight::JetInfo lj(light_mceff,1.08);
1193    
1194     // b-tag systematic UP 9% for b, 18% for c
1195     // for CSVM 5% for b, and 10% for c
1196     //BTagWeight::JetInfo bjUP(0.63,0.99);
1197     //BTagWeight::JetInfo cjUP(0.15,1.07);
1198     BTagWeight::JetInfo bjUP(0.63,1.008);
1199     BTagWeight::JetInfo cjUP(0.15,1.056);
1200    
1201     // b-tag systemacit DOWN 9% for b, 18% for c
1202     //BTagWeight::JetInfo bjDOWN(0.63,0.83);
1203     //BTagWeight::JetInfo cjDOWN(0.15,0.75);
1204     BTagWeight::JetInfo bjDOWN(0.63,0.912);
1205     BTagWeight::JetInfo cjDOWN(0.15,0.864);
1206    
1207     // for high pt jets > 240 UP 50% for b and c
1208     // for CSVM
1209     //BTagWeight::JetInfo bjUPhighpt(0.63,1.36);
1210     //BTagWeight::JetInfo cjUPhighpt(0.15,1.36);
1211     BTagWeight::JetInfo bjUPhighpt(0.63,1.104);
1212     BTagWeight::JetInfo cjUPhighpt(0.15,1.104);
1213     // for high pt jets > 240 DOWN 50% for b and c
1214     //BTagWeight::JetInfo bjDOWNhighpt(0.63,0.46);
1215     //BTagWeight::JetInfo cjDOWNhighpt(0.15,0.46);
1216     BTagWeight::JetInfo bjDOWNhighpt(0.63,0.816);
1217     BTagWeight::JetInfo cjDOWNhighpt(0.15,0.816);
1218    
1219     vector<BTagWeight::JetInfo> j;
1220     for(int i=0;i<number_of_b;i++)j.push_back(bj);
1221     for(int i=0;i<number_of_b_highpt;i++)j.push_back(bj);
1222     for(int i=0;i<number_of_c;i++)j.push_back(cj);
1223     for(int i=0;i<number_of_c_highpt;i++)j.push_back(cj);
1224     for(int i=0;i<number_of_l;i++)j.push_back(lj);
1225    
1226     // changed to CSVM from TCHPM
1227     if (Nbtags_CSVM==0) {
1228     SFb_0tag = b0.weight(j,0);
1229     hjets["Nbtags_CSVM"]->Fill( Nbtags_CSVM, PUweight*SFb_0tag ); // fill bin 0
1230     }
1231    
1232     // only one tag
1233     BTagWeight b11(1,1); // number of tags
1234     if (Nbtags_CSVM==1) {
1235     SFb_only1tag = b11.weight(j,1);
1236     hjets["Nbtags_CSVM"]->Fill( Nbtags_CSVM, PUweight*SFb_only1tag ); // fill bin 1
1237     }
1238    
1239     // at least one tag
1240     BTagWeight b1(1,Nbtags_CSVM); // number of tags
1241     if (Nbtags_CSVM>=1) {
1242     SFb_1tag = b1.weight(j,1);
1243    
1244     // UP
1245     vector<BTagWeight::JetInfo> jj;
1246     for(int i=0;i<number_of_b;i++)jj.push_back(bjUP);
1247     for(int i=0;i<number_of_b_highpt;i++)jj.push_back(bjUPhighpt);
1248     for(int i=0;i<number_of_c;i++)jj.push_back(cjUP);
1249     for(int i=0;i<number_of_c_highpt;i++)jj.push_back(cjUPhighpt);
1250     for(int i=0;i<number_of_l;i++)jj.push_back(lj);
1251     SFb_1tag_syst[0] = b1.weight(jj,1);
1252    
1253     // DOWN
1254     vector<BTagWeight::JetInfo> jk;
1255     for(int i=0;i<number_of_b;i++)jk.push_back(bjDOWN);
1256     for(int i=0;i<number_of_b_highpt;i++)jk.push_back(bjDOWNhighpt);
1257     for(int i=0;i<number_of_c;i++)jk.push_back(cjDOWN);
1258     for(int i=0;i<number_of_c_highpt;i++)jk.push_back(cjDOWNhighpt);
1259     for(int i=0;i<number_of_l;i++)jk.push_back(lj);
1260     SFb_1tag_syst[1] = b1.weight(jk,1);
1261 algomez 1.8
1262 algomez 1.7 }
1263    
1264     // at least two tags
1265     BTagWeight b2(2,Nbtags_CSVM); // number of tags
1266     if (Nbtags_CSVM>=2) {
1267     SFb_2tag = b2.weight(j,2);
1268     hjets["Nbtags_CSVM"]->Fill( Nbtags_CSVM, PUweight*SFb_2tag ); // fill bin >=2
1269    
1270     // UP
1271     vector<BTagWeight::JetInfo> jj;
1272     for(int i=0;i<number_of_b;i++)jj.push_back(bjUP);
1273     for(int i=0;i<number_of_b_highpt;i++)jj.push_back(bjUPhighpt);
1274     for(int i=0;i<number_of_c;i++)jj.push_back(cjUP);
1275     for(int i=0;i<number_of_c_highpt;i++)jj.push_back(cjUPhighpt);
1276     for(int i=0;i<number_of_l;i++)jj.push_back(lj);
1277     SFb_2tag_syst[0] = b2.weight(jj,2);
1278    
1279     // DOWN
1280     vector<BTagWeight::JetInfo> jk;
1281     for(int i=0;i<number_of_b;i++)jk.push_back(bjDOWN);
1282     for(int i=0;i<number_of_b_highpt;i++)jk.push_back(bjDOWNhighpt);
1283     for(int i=0;i<number_of_c;i++)jk.push_back(cjDOWN);
1284     for(int i=0;i<number_of_c_highpt;i++)jk.push_back(cjDOWNhighpt);
1285     for(int i=0;i<number_of_l;i++)jk.push_back(lj);
1286     SFb_2tag_syst[1] = b2.weight(jk,2);
1287 algomez 1.9 }*/
1288     //else hjets["Nbtags_CSVM"]->Fill( Nbtags_CSVM );
1289 algomez 1.7
1290     //hjets["Nbtags_TCHPM"]->Fill( Nbtags_TCHPM, PUweight*SFb );
1291     //hjets["Nbtags_CSVM"]->Fill( Nbtags_CSVM, PUweight*SFb );
1292 algomez 1.6
1293 algomez 1.9
1294 algomez 1.7 // Cuts
1295 algomez 1.6 bool passcut = true;
1296     if ( Ht <= 300. ) passcut = false;
1297    
1298     if (passcut) {
1299 algomez 1.7 hPVs["Nreweight_2jet"]->Fill( total_pvs, PUweight );
1300     hMET["Ht"]->Fill( Ht, PUweight );
1301 algomez 1.10 hMET["MET"]->Fill( p4MET.Pt(), PUweight );
1302     hMET["Stlep"]->Fill( Stlep , PUweight );
1303     hMET["Stjet"]->Fill( Stjet , PUweight );
1304 algomez 1.6 hjets["Njets_cut1"]->Fill( njets, PUweight );
1305     hjets["1st_pt"]->Fill( p4jets[0].Pt(), PUweight );
1306     hjets["1st_eta"]->Fill( p4jets[0].Eta(), PUweight );
1307     hjets["2nd_pt"]->Fill( p4jets[1].Pt(), PUweight );
1308     hjets["2nd_eta"]->Fill( p4jets[1].Eta(), PUweight );
1309     hjets["3rd_pt"]->Fill( p4jets[2].Pt(), PUweight );
1310     hjets["4th_pt"]->Fill( p4jets[3].Pt(), PUweight );
1311     hjets["5th_pt"]->Fill( p4jets[4].Pt(), PUweight );
1312     hjets["6th_pt"]->Fill( p4jets[5].Pt(), PUweight );
1313     hjets["7th_pt"]->Fill( p4jets[6].Pt(), PUweight );
1314 algomez 1.10 hjets["Dijet_deltaR"]->Fill( deltaRjj, PUweight );
1315 algomez 1.6 hmuons["N"]->Fill( total_muons, PUweight );
1316     hmuons["Nelectrons"]->Fill( nlooseelectrons, PUweight );
1317 algomez 1.10 hmuons["pt"]->Fill( p4lepton.Pt(), PUweight );
1318 algomez 1.6 hmuons["deltaR"]->Fill( deltaR, PUweight );
1319 algomez 1.7 hM["WMt_2jet"]->Fill( WMt, PUweight );
1320    
1321 algomez 1.11 // B-tagging
1322     hjets["Nbtags_CSVL"]->Fill( Nbtags_CSVL, PUweight );
1323     hjets["Nbtags_CSVM"]->Fill( Nbtags_CSVM, PUweight );
1324     hjets["Nbtags_CSVT"]->Fill( Nbtags_CSVT, PUweight );
1325    
1326 algomez 1.10 // plot bdiscriminator
1327     if ( bdiscriminator[0] >= 0 ) { MyStoreTree->GetJetVariable()->bdisc_1st = bdiscriminator[0]; hjets["1st_bdisc"]->Fill( bdiscriminator[0], PUweight );}
1328     if ( bdiscriminator[1] >= 0 ) { MyStoreTree->GetJetVariable()->bdisc_2nd = bdiscriminator[1]; hjets["2nd_bdisc"]->Fill( bdiscriminator[1], PUweight );}
1329     if ( bdiscriminator[2] >= 0 ) { MyStoreTree->GetJetVariable()->bdisc_3rd = bdiscriminator[2]; hjets["3rd_bdisc"]->Fill( bdiscriminator[2], PUweight );}
1330     if ( bdiscriminator[3] >= 0 ) { MyStoreTree->GetJetVariable()->bdisc_4th = bdiscriminator[3]; hjets["4th_bdisc"]->Fill( bdiscriminator[3], PUweight );}
1331 algomez 1.11 //prod_bdisc = bdiscriminator[0]*bdiscriminator[1]*bdiscriminator[2]*bdiscriminator[3];
1332     //hjets["prod_bdisc"]->Fill( prod_bdisc, PUweight );
1333    
1334 algomez 1.10 MyStoreTree->GetGeneralVariable()->PUWeight = PUweight;
1335     MyStoreTree->GetMetVariable()->Ht = Ht;
1336     MyStoreTree->GetMetVariable()->Stlep = Stlep;
1337     MyStoreTree->GetMetVariable()->Stjet = Stjet;
1338     MyStoreTree->GetJetVariable()->numjets= njets;
1339 algomez 1.11 MyStoreTree->GetJetVariable()->numBjets_csvl= Nbtags_CSVL;
1340     MyStoreTree->GetJetVariable()->numBjets_csvm= Nbtags_CSVM;
1341     MyStoreTree->GetJetVariable()->numBjets_csvt= Nbtags_CSVT;
1342 algomez 1.10 //////////////////////////////////////////////////////////////////////////////////
1343    
1344 algomez 1.7 if ( Nbtags_CSVM >= 1 ) {
1345 algomez 1.10 hPVs["Nreweight_2jet_1btag"]->Fill( total_pvs, PUweight );
1346     hMET["Ht_1btag"]->Fill( Ht, PUweight );
1347     hMET["MET_1btag"]->Fill( p4MET.Pt(), PUweight );
1348     hMET["Stlep_1btag"]->Fill( Stlep , PUweight );
1349     hMET["Stjet_1btag"]->Fill( Stjet , PUweight );
1350     hjets["Njets_1btag"]->Fill(njets, PUweight*SFb_1tag );
1351     hjets["1st_pt_1btag"]->Fill( p4jets[0].Pt(), PUweight );
1352     hjets["1st_eta_1btag"]->Fill( p4jets[0].Eta(), PUweight );
1353     hjets["2nd_pt_1btag"]->Fill( p4jets[1].Pt(), PUweight );
1354     hjets["2nd_eta_1btag"]->Fill( p4jets[1].Eta(), PUweight );
1355     hjets["3rd_pt_1btag"]->Fill( p4jets[2].Pt(), PUweight );
1356     hjets["4th_pt_1btag"]->Fill( p4jets[3].Pt(), PUweight );
1357     hjets["5th_pt_1btag"]->Fill( p4jets[4].Pt(), PUweight );
1358     hjets["6th_pt_1btag"]->Fill( p4jets[5].Pt(), PUweight );
1359     hjets["7th_pt_1btag"]->Fill( p4jets[6].Pt(), PUweight );
1360     hjets["Dijet_deltaR_1btag"]->Fill( deltaRjj, PUweight );
1361     hmuons["N_1btag"]->Fill( total_muons, PUweight );
1362     hmuons["Nelectrons_1btag"]->Fill( nlooseelectrons, PUweight );
1363     hmuons["pt_1btag"]->Fill( p4lepton.Pt(), PUweight );
1364     hmuons["deltaR_1btag"]->Fill( deltaR, PUweight );
1365     hM["WMt_2jet_1btag"]->Fill( WMt, PUweight );
1366    
1367     // plot bdiscriminator
1368     //if ( bdiscriminator[0] >= 0 ) hjets["1st_bdisc"]->Fill( bdiscriminator[0], PUweight );
1369     hjets["1st_bdisc_1btag"]->Fill( bdiscriminator[0], PUweight );
1370     if ( bdiscriminator[1] >= 0 ) hjets["2nd_bdisc_1btag"]->Fill( bdiscriminator[1], PUweight );
1371     if ( bdiscriminator[2] >= 0 ) hjets["3rd_bdisc_1btag"]->Fill( bdiscriminator[2], PUweight );
1372     if ( bdiscriminator[3] >= 0 ) hjets["4th_bdisc_1btag"]->Fill( bdiscriminator[3], PUweight );
1373 algomez 1.11 //hjets["prod_bdisc_1btag"]->Fill( prod_bdisc, PUweight );
1374 algomez 1.10 //////////////////////////////////////////////////////////////////////////////////
1375 algomez 1.7 }
1376 algomez 1.4 }
1377 algomez 1.6 }
1378 algomez 1.1
1379 algomez 1.9 MyStoreTree->GetGeneralVariable()->Run = ntuple->run;
1380     MyStoreTree->GetGeneralVariable()->Lumi = ntuple->lumi;
1381     MyStoreTree->GetGeneralVariable()->Event = ntuple->event;
1382 algomez 1.8
1383     MyStoreTree->GetStoreTree()->Fill();
1384 algomez 1.1 if (fVerbose) cout << "done analysis" << endl;
1385     return kTRUE;
1386     }
1387    
1388     void Analyzer::SlaveTerminate()
1389     {
1390     // The SlaveTerminate() function is called after all entries or objects
1391     // have been processed. When running with PROOF SlaveTerminate() is called
1392     // on each slave server.
1393    
1394     // fill cutflow histogram
1395    
1396     int ibin = 1;
1397     for ( vector<string>::const_iterator ivec= fCutLabels.begin(); ivec != fCutLabels.end(); ++ivec )
1398     // for ( map<string, int >::const_iterator imap=cutmap.begin(); imap!=cutmap.end(); ++imap )
1399     {
1400     hcutflow->SetBinContent( ibin, cutmap[ *ivec ] );
1401     ibin++;
1402     }
1403 algomez 1.8
1404 algomez 1.1 // Write the ntuple to the file
1405     if (fFile) {
1406 algomez 1.8 Bool_t cleanup = kFALSE;
1407     TDirectory *savedir = gDirectory;
1408     if ( h1test->GetEntries() > 0) {
1409     fFile->cd();
1410     h1test->Write();
1411     hcutflow->Write();
1412     MyStoreTree->GetStoreTree()->Write();
1413     //h2_pt_Wprime->Write();
1414     fFile->mkdir("muons");
1415     fFile->cd("muons");
1416     for ( map<string,TH1* >::const_iterator imap=hmuons.begin(); imap!=hmuons.end(); ++imap ) {
1417     TH1 *temp = imap->second;
1418     if ( temp->GetEntries() > 0 )
1419     temp->Write();
1420     //else cout << "Warning: empty histogram " << temp->GetName() << " will not be written to file." << endl;
1421     }
1422     fFile->cd();
1423     fFile->mkdir("PVs");
1424     fFile->cd("PVs");
1425     for ( map<string,TH1* >::const_iterator imap=hPVs.begin(); imap!=hPVs.end(); ++imap ) {
1426     TH1 *temp = imap->second;
1427     if ( temp->GetEntries() > 0 )
1428     temp->Write();
1429     //else cout << "Warning: empty histogram " << temp->GetName() << " will not be written to file." << endl;
1430     }
1431     fFile->cd();
1432     fFile->mkdir("jets");
1433     fFile->cd("jets");
1434     for ( map<string,TH1* >::const_iterator imap=hjets.begin(); imap!=hjets.end(); ++imap ) {
1435     TH1 *temp = imap->second;
1436     if ( temp->GetEntries() > 0 )
1437     temp->Write();
1438     }
1439     fFile->cd();
1440     fFile->mkdir("mass");
1441     fFile->cd("mass");
1442     for ( map<string,TH1* >::const_iterator imap=hM.begin(); imap!=hM.end(); ++imap ) {
1443     TH1 *temp = imap->second;
1444     if ( temp->GetEntries() > 0 )
1445     temp->Write();
1446     }
1447     fFile->cd();
1448     fFile->mkdir("MET");
1449     fFile->cd("MET");
1450     for ( map<string,TH1* >::const_iterator imap=hMET.begin(); imap!=hMET.end(); ++imap ) {
1451     TH1 *temp = imap->second;
1452     if ( temp->GetEntries() > 0 )
1453     temp->Write();
1454     }
1455     fFile->cd();
1456     // fFile->mkdir("TMVA");
1457     // fFile->cd("TMVA");
1458     // fBdisc->Write();
1459     // fBdisc->Print();
1460     // TMVA->Write();
1461     //
1462     fFile->cd();
1463    
1464     fProofFile->Print();
1465     fOutput->Add(fProofFile);
1466     } else {
1467     cleanup = kTRUE;
1468     }
1469 algomez 1.1
1470    
1471     h1test->SetDirectory(0);
1472     hcutflow->SetDirectory(0);
1473     gDirectory = savedir;
1474     fFile->Close();
1475     // Cleanup, if needed
1476     if (cleanup) {
1477 algomez 1.8 TUrl uf(*(fFile->GetEndpointUrl()));
1478     SafeDelete(fFile);
1479     gSystem->Unlink(uf.GetFile());
1480     SafeDelete(fProofFile);
1481 algomez 1.1 }
1482     }
1483    
1484     }
1485    
1486     void Analyzer::Terminate()
1487     {
1488     // The Terminate() function is the last function to be called during
1489     // a query. It always runs on the client, it can be used to present
1490     // the results graphically or save the results to file.
1491    
1492     Info("Terminate","Analyzer done.");
1493     }