ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/algomez/FourTop/TreeAnalyzer/test/Analyzer.C
Revision: 1.6
Committed: Fri Dec 16 15:46:26 2011 UTC (13 years, 4 months ago) by algomez
Content type: text/plain
Branch: MAIN
Changes since 1.5: +249 -661 lines
Log Message:
first cuts

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     #include "BTagWeight.h"
28     //#include "Yumiceva/TreeAnalyzer/interface/JetCombinatorics.h"
29    
30     #include <TStyle.h>
31     #include <TSystem.h>
32    
33     #include <iostream>
34     #include <vector>
35     #include <stdio.h>
36     #include <stdlib.h>
37    
38     void Analyzer::ParseInput()
39     {
40    
41     if (fMyOpt.Contains("muon"))
42     {
43     fChannel = 1;
44     }
45     if (fMyOpt.Contains("electron"))
46     {
47     fChannel = 2;
48     }
49     if (fMyOpt.Contains("verbose"))
50     {
51     fVerbose = true;
52     }
53     if (fMyOpt.Contains("JECUP")) { fdoJECunc = true; fdoJECup = true; }
54     if (fMyOpt.Contains("JECDOWN")) { fdoJECunc = true; fdoJECup = false; }
55 algomez 1.2 if (fMyOpt.Contains("PUUP")) { fpuhistogram = "WHistUp";}
56     if (fMyOpt.Contains("PUDOWN")) { fpuhistogram = "WHistDown";}
57 algomez 1.1 if (fMyOpt.Contains("QCD1")) fdoQCD1SideBand = true;
58     if (fMyOpt.Contains("QCD2")) fdoQCD2SideBand = true;
59     if (fMyOpt.Contains("mtop")) fdoMtopCut = true;
60     if (fMyOpt.Contains("outdir")) {
61     TString tmp = fMyOpt;
62     tmp = tmp.Remove(0,fMyOpt.Index("outdir")+7);
63     fOutdir = tmp+"/";
64     Info("Begin","output files will be written to directory: %s", fOutdir.Data());
65     }
66     if (fMyOpt.Contains("sample"))
67     {
68     TString tmp = fMyOpt;
69     tmp = tmp.Remove(0,fMyOpt.Index("sample")+7);
70     fSample = tmp;
71     if (fdoJECunc && fdoJECup==true) fSample += "_JECUP";
72     if (fdoJECunc && fdoJECup==false) fSample += "_JECDOWN";
73 algomez 1.2 if (fMyOpt.Contains("PUUP")) fSample += "_PUUP";
74     if (fMyOpt.Contains("PUDOWN")) fSample += "_PUDOWN";
75 algomez 1.1
76     Info("Begin","Histogram names will have suffix: %s", fSample.Data());
77    
78     if ( fSample.Contains("data") )
79     {
80     fIsMC = false;
81     Info("Begin","This sample is treated as DATA");
82     if (fdoQCD1SideBand) fSample = "dataQCD1";
83     if (fdoQCD2SideBand) fSample = "dataQCD2";
84     }
85 algomez 1.2 else
86     {
87     Info("Begin","This sample is treated as MC");
88     Info("Begin","we will use pileup set: %s", fpuhistogram.Data());
89     }
90    
91 algomez 1.1
92     }
93     }
94    
95     void Analyzer::Begin(TTree * /*tree*/)
96     {
97     // The Begin() function is called at the start of the query.
98     // When running with PROOF Begin() is only called on the client.
99     // The tree argument is deprecated (on PROOF 0 is passed).
100    
101     TString option = GetOption();
102     fMyOpt = option;
103     ParseInput();
104    
105     Info("Begin", "starting with process option: %s", option.Data());
106    
107     }
108    
109     void Analyzer::SlaveBegin(TTree * tree)
110     {
111    
112     // The SlaveBegin() function is called after the Begin() function.
113     // When running with PROOF SlaveBegin() is called on each slave server.
114     // The tree argument is deprecated (on PROOF 0 is passed).
115    
116     TString option = GetOption();
117     fMyOpt = option;
118     ParseInput();
119    
120     Info("SlaveBegin",
121     "starting with process option: %s (tree: %p)", option.Data(), tree);
122    
123     //initialize the Tree branch addresses
124     Init(tree);
125    
126     // We may be creating a dataset or a merge file: check it
127     TNamed *nm = dynamic_cast<TNamed *>(fInput->FindObject("SimpleNtuple.root"));
128     if (nm) {
129     // Just create the object
130     UInt_t opt = TProofOutputFile::kRegister | TProofOutputFile::kOverwrite | TProofOutputFile::kVerify;
131     fProofFile = new TProofOutputFile("SimpleNtuple.root",
132     TProofOutputFile::kDataset, opt, nm->GetTitle());
133     } else {
134     // For the ntuple, we use the automatic file merging facility
135     // Check if an output URL has been given
136     TNamed *out = (TNamed *) fInput->FindObject("PROOF_OUTPUTFILE_LOCATION");
137     Info("SlaveBegin", "PROOF_OUTPUTFILE_LOCATION: %s", (out ? out->GetTitle() : "undef"));
138     TString tmpfilename = "results";
139     if ( fSample != "" ) tmpfilename += "_"+fSample+".root";
140     else tmpfilename = "results.root";
141     fProofFile = new TProofOutputFile(tmpfilename, (out ? out->GetTitle() : "M"));
142     out = (TNamed *) fInput->FindObject("PROOF_OUTPUTFILE");
143     if (out) fProofFile->SetOutputFileName(fOutdir + out->GetTitle());
144     }
145    
146     // Open the file
147     //TDirectory *savedir = gDirectory;
148     if (!(fFile = fProofFile->OpenFile("RECREATE"))) {
149     Warning("SlaveBegin", "problems opening file: %s/%s",
150     fProofFile->GetDir(), fProofFile->GetFileName());
151     }
152    
153 algomez 1.2 // Get PU weights
154 algomez 1.6 TString weightfilename = "/uscms/home/algomez/work/CMSSW_4_2_4/src/Yumiceva/TreeAnalyzer/test/Weight3Dfinebin4p7.root"; //Weight3D.root";
155 algomez 1.2 fweightfile = new TFile(weightfilename,"read");
156     f3Dweight = (TH1D*) fweightfile->Get(fpuhistogram);
157    
158 algomez 1.1 //create histograms
159     h1test = new TH1F("h1test","muon p_{T}",100,10.,400);
160     //fHist = new HistoManager(string(fSample));
161     TString hname = "_"+fSample;
162    
163 algomez 1.6 hmuons["N_cut0"] = new TH1F("N_muons_cut0"+hname,"Number of Muons",6,-0.5,5.5);
164     hmuons["N"] = new TH1F("N_muons"+hname,"Number of Muons",6,-0.5,5.5);
165     hmuons["Nelectrons_cut0"] = new TH1F("Nelectrons_cut0"+hname,"Number of Loose Electrons",6,-0.5,5.5);
166     hmuons["Nelectrons"] = new TH1F("Nelectrons"+hname,"Number of Loose Electrons",6,-0.5,5.5);
167 algomez 1.1 hmuons["pt_1jet"] = new TH1F("muon_pt_1jet"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
168     hmuons["pt_2jet"] = new TH1F("muon_pt_2jet"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
169 algomez 1.2 hmuons["pt_fin"] = new TH1F("muon_pt_fin"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
170 algomez 1.1 hmuons["eta"] = new TH1F("muon_eta"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
171     hmuons["eta_1jet"] = new TH1F("muon_eta_1jet"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
172     hmuons["eta_2jet"] = new TH1F("muon_eta_2jet"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
173 algomez 1.2 hmuons["eta_fin"] = new TH1F("muon_eta_fin"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
174 algomez 1.1 hmuons["phi"] = new TH1F("muon_phi"+hname,"#phi^{#mu}", 30, -3.15, 3.15);
175     hmuons["phi_1jet"] = new TH1F("muon_phi_1jet"+hname,"#phi^{#mu}", 30, -3.15, 3.15);
176     hmuons["phi_2jet"] = new TH1F("muon_phi_2jet"+hname,"#phi^{#mu}", 30, -3.15, 3.15);
177 algomez 1.2 hmuons["phi_fin"] = new TH1F("muon_phi_fin"+hname,"#phi^{#mu}", 30, -3.15, 3.15);
178 algomez 1.1 hmuons["reliso"] = new TH1F("muon_reliso"+hname,"Relative Isolation", 40, 0, 0.2);
179 algomez 1.2 hmuons["reliso_fin"] = new TH1F("muon_reliso_fin"+hname,"Relative Isolation", 40, 0, 0.2);
180 algomez 1.1 hmuons["deltaR_cut0"] = new TH1F("deltaR_cut0"+hname,"#DeltaR(#mu,jet)",80, 0, 4);
181     hmuons["deltaR"] = new TH1F("deltaR"+hname,"#DeltaR(#mu,jet)", 80, 0, 4);
182     hmuons["d0_cut1"] = new TH1F("d0_cut1"+hname,"#mu Impact Parameter [cm]",20,-0.1,0.1);
183     hmuons["pt_cut1"] = new TH1F("muon_pt_cut1"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
184     hmuons["pt_cut2"] = new TH1F("muon_pt_cut2"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
185     hmuons["pt_cut3"] = new TH1F("muon_pt_cut3"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
186     hmuons["pt"] = new TH1F("muon_pt"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
187     hmuons["dz"] = new TH1F("dz"+hname,"|z(#mu) - z_{PV}| [cm]", 25, 0, 1.);
188 algomez 1.2 hmuons["Niso"] = new TH1F("Niso"+hname,"Number of Primary Vertices", 25,-0.5,24.5);
189     hmuons["Ngood"] = new TH1F("Ngood"+hname,"Number of Primary Vertices",25,-0.5,24.5);
190 algomez 1.1
191     hPVs["N"] = new TH1F("NPV"+hname,"Number of PVs",25,-0.5,24.5);
192     hPVs["Nreweight"] = new TH1F("NPVreweight"+hname,"Number of PVs",25,-0.5,24.5);
193     hPVs["Nreweight_2jet"] = new TH1F("NPVreweight_2jet"+hname,"Number of PVs",25,-0.5,24.5);
194    
195     helectrons["pt"] = new TH1F("electron_pt"+hname,"p_{T}^{#mu} [GeV/c]", 50, 20, 500);
196     helectrons["pt_cut2"] = new TH1F("electron_pt_cut2"+hname,"p_{T}^{#mu} [GeV/c]", 50, 20, 500);
197     helectrons["pt_1jet"] = new TH1F("electron_pt_1jet"+hname,"p_{T}^{#mu} [GeV/c]", 50, 20, 500);
198     helectrons["pt_2jet"] = new TH1F("electron_pt_2jet"+hname,"p_{T}^{#mu} [GeV/c]", 50, 20, 500);
199     helectrons["pt_3jet"] = new TH1F("electron_pt_3jet"+hname,"p_{T}^{#mu} [GeV/c]", 50, 20, 500);
200     helectrons["pt_4jet"] = new TH1F("electron_pt_4jet"+hname,"p_{T}^{#mu} [GeV/c]", 50, 20, 500);
201     helectrons["eta"] = new TH1F("electron_eta"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
202     helectrons["eta_cut2"] = new TH1F("electron_eta_cut2"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
203     helectrons["eta_1jet"] = new TH1F("electron_eta_1jet"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
204     helectrons["eta_2jet"] = new TH1F("electron_eta_2jet"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
205     helectrons["eta_3jet"] = new TH1F("electron_eta_3jet"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
206     helectrons["eta_4jet"] = new TH1F("electron_eta_4jet"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
207     helectrons["phi"] = new TH1F("electron_phi"+hname,"#phi^{#mu}", 20, 0, 3.15);
208     helectrons["phi_cut2"] = new TH1F("electron_phi_cut2"+hname,"#phi^{#mu}", 20, 0, 3.15);
209     helectrons["phi_1jet"] = new TH1F("electron_phi_1jet"+hname,"#phi^{#mu}", 20, 0, 3.15);
210     helectrons["phi_2jet"] = new TH1F("electron_phi_2jet"+hname,"#phi^{#mu}", 20, 0, 3.15);
211     helectrons["phi_3jet"] = new TH1F("electron_phi_3jet"+hname,"#phi^{#mu}", 20, 0, 3.15);
212     helectrons["phi_4jet"] = new TH1F("electron_phi_4jet"+hname,"#phi^{#mu}", 20, 0, 3.15);
213     helectrons["reliso"] = new TH1F("electron_reliso"+hname,"Relative Isolation", 40, 0, 0.2);
214     helectrons["reliso_1jet"] = new TH1F("electron_reliso_1jet"+hname,"Relative Isolation", 40, 0, 0.2);
215     helectrons["reliso_2jet"] = new TH1F("electron_reliso_2jet"+hname,"Relative Isolation", 40, 0, 0.2);
216     helectrons["reliso_3jet"] = new TH1F("electron_reliso_3jet"+hname,"Relative Isolation", 40, 0, 0.2);
217     helectrons["reliso_4jet"] = new TH1F("electron_reliso_4jet"+hname,"Relative Isolation", 40, 0, 0.2);
218     helectrons["deltaR_cut0"] = new TH1F("electron_deltaR_cut0"+hname,"#DeltaR(#mu,jet)",80, 0, 4);
219     helectrons["deltaR"] = new TH1F("electron_deltaR"+hname,"#DeltaR(#mu,jet)",80, 0, 4);
220     helectrons["d0_cut1"] = new TH1F("electron_d0_cut1"+hname,"#mu Impact Parameter [cm]",20,-0.1,0.1);
221     helectrons["dz"] = new TH1F("electron_dz"+hname,"|z(#mu) - z_{PV}| [cm]", 25, 0, 1.);
222    
223 algomez 1.6 hMET["MET"] = new TH1F("MET"+hname,"Missing Transverse Energy [GeV]", 50, 0, 500);
224     hMET["MET_cut0"] = new TH1F("MET_cut0"+hname,"Missing Transverse Energy [GeV]", 50, 0, 300);
225 algomez 1.1 hMET["MET_2jet"] = new TH1F("MET_2jet"+hname,"Missing Transverse Energy [GeV]", 50, 0, 300);
226 algomez 1.2 hMET["MET_fin"] = new TH1F("MET_fin"+hname,"Missing Transverse Energy [GeV]", 50, 0, 300);
227 algomez 1.1 hMET["genMET_2jet"] = new TH1F("genMET_2jet"+hname,"Missing Transverse Energy [GeV]", 50, 0, 300);
228     hMET["deltaMET_2jet"] = new TH1F("deltaMET_2jet"+hname,"Missing Transverse Energy [GeV]", 50, -200, 200);
229     hMET["phi"] = new TH1F("MET_phi"+hname,"#phi Missing Transverse Energy [GeV]", 20, 0, 3.15);
230 algomez 1.6 hMET["Ht0"] = new TH1F("Ht0"+hname,"H_{T} [GeV]", 50, 0, 3000);
231     hMET["Ht"] = new TH1F("Ht"+hname,"H_{T} [GeV]", 50, 0, 3000);
232 algomez 1.1 hMET["Htlep"] = new TH1F("Htlep"+hname,"H_{T,lep} [GeV]", 50, 0, 1000);
233     hMET["PzNu"] = new TH1F("PzNu"+hname,"p_{z} #nu [GeV]", 40, -300,300);
234     hMET["EtaNu"] = new TH1F("EtaNu"+hname,"#eta",50,-2.2,2.2);
235     hMET["LepWmass"] = new TH1F("LepWmass"+hname,"W#rightarrow#mu#nu Mass [GeV/c^{2}]",20, 0, 150);
236     hMET["LepWmass_topcut"] = new TH1F("LepWmass_topcut"+hname,"W#rightarrow#mu#nu Mass [GeV/c^{2}]",20, 0, 150);
237     hMET["LepWmassNoPt"]=new TH1F("LepWmassNoPt"+hname,"W#rightarrow#mu#nu Mass [GeV/c^{2}]",20, 0, 150);
238     hMET["deltaPhi"] = new TH1F("deltaPhi"+hname,"#Delta #phi(#mu,MET)",50, -3.15, 3.15);
239    
240 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)
241 algomez 1.1 hM["WMt_2jet"] = new TH1F("Mt_2jet"+hname,"M_{T}(W) [GeV/c^{2}]", 50, 0, 300);
242 algomez 1.4 hM["WMt_5jet"] = new TH1F("Mt_5jet"+hname,"M_{T}(W) [GeV/c^{2}]", 50, 0, 300); // MT for njets >= 5
243 algomez 1.2 hM["WMt_fin"] = new TH1F("Mt_fin"+hname,"M_{T}(W) [GeV/c^{2}]", 50, 0, 300);
244 algomez 1.1 hM["dijet"] = new TH1F("dijet"+hname,"(jj) mass [GeV/c^{2}]", 50, 0, 1000);
245     hM["trijet"] = new TH1F("trijet"+hname,"(jjj) mass [GeV/c^{2}]", 50, 0, 1000);
246     hM["top_1btag"] = new TH1F("top_1btag"+hname,"top mass [GeV/c^{2}]",50,0,500);
247     hM["top_2btag"] = new TH1F("top_2btag"+hname,"top mass [GeV/c^{2}]",50,0,500);
248 algomez 1.6
249 algomez 1.1 hjets["pt"] = new TH1F("jet_pt"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
250 algomez 1.2 hjets["pt_b_mc"] = new TH1F("jet_pt_b_mc"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
251     hjets["pt_c_mc"] = new TH1F("jet_pt_c_mc"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
252     hjets["pt_l_mc"] = new TH1F("jet_pt_l_mc"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
253 algomez 1.1 hjets["pt_btag"] = new TH1F("jet_pt_btag"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
254     hjets["pt_btag_b"] = new TH1F("jet_pt_btag_b"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
255     hjets["pt_btag_c"] = new TH1F("jet_pt_btag_c"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
256     hjets["pt_btag_l"] = new TH1F("jet_pt_btag_l"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
257 algomez 1.6 hjets["1st_pt_cut0"] = new TH1F("jet1_pt_cut0"+hname,"Leading jet p_{T} [GeV/c]",60, 30, 1000);
258     hjets["1st_pt"] = new TH1F("jet1_pt"+hname,"Leading jet p_{T} [GeV/c]",60, 30, 1000);
259     hjets["2nd_pt_cut0"] = new TH1F("jet2_pt_cut0"+hname,"2nd jet p_{T} [GeV/c]",60, 30, 1000);
260     hjets["2nd_pt"] = new TH1F("jet2_pt"+hname,"2nd jet p_{T} [GeV/c]",60, 30, 1000);
261     hjets["3rd_pt_cut0"] = new TH1F("jet3_pt_cut0"+hname,"3rd jet p_{T} [GeV/c]",60, 30, 800);
262 algomez 1.1 hjets["3rd_pt"] = new TH1F("jet3_pt"+hname,"3rd jet p_{T} [GeV/c]",60, 30, 800);
263 algomez 1.6 hjets["4th_pt_cut0"] = new TH1F("jet4_pt_cut0"+hname,"4th jet p_{T} [GeV/c]",60, 30, 800);
264 algomez 1.1 hjets["4th_pt"] = new TH1F("jet4_pt"+hname,"4th jet p_{T} [GeV/c]",60, 30, 800);
265 algomez 1.6 hjets["5th_pt_cut0"] = new TH1F("jet5_pt_cut0"+hname,"5th jet p_{T} [GeV/c]",60, 30, 800);
266     hjets["5th_pt"] = new TH1F("jet5_pt"+hname,"5th jet p_{T} [GeV/c]",60, 30, 800);
267     hjets["6th_pt_cut0"] = new TH1F("jet6_pt_cut0"+hname,"6th jet p_{T} [GeV/c]",60, 30, 800);
268     hjets["6th_pt"] = new TH1F("jet6_pt"+hname,"6th jet p_{T} [GeV/c]",60, 30, 800);
269     hjets["7th_pt_cut0"] = new TH1F("jet7_pt_cut0"+hname,"7th jet p_{T} [GeV/c]",60, 30, 800);
270     hjets["7th_pt"] = new TH1F("jet7_pt"+hname,"7th jet p_{T} [GeV/c]",60, 30, 800);
271 algomez 1.2 hjets["1st_pt_fin"] = new TH1F("jet1_pt_fin"+hname,"Leading jet p_{T} [GeV/c]",60, 30, 800);
272     hjets["2nd_pt_fin"] = new TH1F("jet2_pt_fin"+hname,"2nd jet p_{T} [GeV/c]",60, 30, 800);
273 algomez 1.1 hjets["eta"] = new TH1F("jet_eta"+hname,"jet #eta",50, -2.4, 2.4);
274 algomez 1.6 hjets["1st_eta_cut0"] = new TH1F("jet1_eta_cut0"+hname,"Leading jet #eta",50, -2.4, 2.4);
275 algomez 1.1 hjets["1st_eta"] = new TH1F("jet1_eta"+hname,"Leading jet #eta",50, -2.4, 2.4);
276 algomez 1.6 hjets["2nd_eta_cut0"] = new TH1F("jet2_eta_cut0"+hname,"2nd jet #eta",50, -2.4, 2.4);
277 algomez 1.1 hjets["2nd_eta"] = new TH1F("jet2_eta"+hname,"2nd jet #eta",50, -2.4, 2.4);
278 algomez 1.2 hjets["1st_eta_fin"] = new TH1F("jet1_eta_fin"+hname,"Leading jet #eta",50, -2.4, 2.4);
279     hjets["2nd_eta_fin"] = new TH1F("jet2_eta_fin"+hname,"2nd jet #eta",50, -2.4, 2.4);
280 algomez 1.1 hjets["phi"] = new TH1F("jet_phi"+hname,"jet #phi",50, 0, 3.15);
281 algomez 1.3 hjets["Njets"] = new TH1F("Njets"+hname,"jet multiplicity",13,-0.5,12.5);
282 algomez 1.6 hjets["Njets_cut0"] = new TH1F("Njets_cut0"+hname,"jet multiplicity",13,-0.5,12.5);
283 algomez 1.4 hjets["Njets_cut1"] = new TH1F("Njets_cut1"+hname,"jet multiplicity",13,-0.5,12.5);
284 algomez 1.1 hjets["Njets_1tag"] = new TH1F("Njets_1tag"+hname,"jet multiplicity",6,0.5,6.5);
285     hjets["Njets_2tag"] = new TH1F("Njets_2tag"+hname,"jet multiplicity",6,0.5,6.5);
286     hjets["Nbtags_TCHPM"] = new TH1F("Nbjets_TCHPM_N2j"+hname,"Tagged b-jets",3,-0.5,2.5);
287 algomez 1.2 hjets["Nbtags_CSVM"] = new TH1F("Nbjets_CSVM_N2j"+hname,"Tagged b-jets",3,-0.5,2.5);
288 algomez 1.1 hjets["Dijet_deltaPhi"] = new TH1F("Dijet_deltaPhi"+hname,"#Delta #phi(j,j)",30,-3.15,3.15);
289 algomez 1.6 hjets["Dijet_deltaR_cut0"] = new TH1F("Dijet_deltaR_cut0"+hname,"#DeltaR(j,j)",80,0.,4.);
290     hjets["Dijet_deltaR"] = new TH1F("Dijet_deltaR"+hname,"#DeltaR(j,j)",80,0.,4.);
291 algomez 1.1 hjets["tb_deltaPhi"] = new TH1F("tb_deltaPhi"+hname,"#Delta #phi(t,b)",30,0.,3.15);
292 algomez 1.2 hjets["tb_deltaEta"] = new TH1F("tb_deltaEta"+hname,"#Delta #eta(t,b)",30,-5,5);
293     hjets["pt_top"] = new TH1F("pt_top"+hname,"top p_{T} [GeV]",60,0,1500);
294     hjets["pt_b"] = new TH1F("pt_b"+hname,"b-jet p_{T} [GeV]",60,0,1500);
295     hjets["gen_deltaR_mub"] = new TH1F("gen_deltaR_mub"+hname,"#Delta R(#mu,b)",40,0,4);
296 algomez 1.1
297     map<string,TH1* > allhistos = hmuons;
298     allhistos.insert( helectrons.begin(), helectrons.end() );
299     allhistos.insert( hMET.begin(), hMET.end() );
300     allhistos.insert( hM.begin(), hM.end() );
301     allhistos.insert( hjets.begin(), hjets.end() );
302    
303     for ( map<string,TH1* >::const_iterator imap=allhistos.begin(); imap!=allhistos.end(); ++imap )
304     {
305     TH1 *temp = imap->second;
306     temp->Sumw2();
307     temp->SetXTitle( temp->GetTitle() );
308     }
309    
310     // cut flow
311     if (fChannel==1)
312     { //muon +jets
313     fCutLabels.push_back("Processed");
314     fCutLabels.push_back("OneIsoMu");
315     fCutLabels.push_back("LooseMuVeto");
316     fCutLabels.push_back("ElectronVeto");
317     fCutLabels.push_back("MET");
318     fCutLabels.push_back("1Jet");
319     fCutLabels.push_back("2Jet");
320     fCutLabels.push_back("3Jet");
321     fCutLabels.push_back("4Jet");
322     fCutLabels.push_back("2Jet1b");
323     fCutLabels.push_back("2Jet2b");
324     fCutLabels.push_back("MaxJets");
325     fCutLabels.push_back("phi");
326     fCutLabels.push_back("topmass");
327     }
328     else
329     { //electron+jets
330     fCutLabels.push_back("Processed");
331     fCutLabels.push_back("OneIsoElectron");
332     fCutLabels.push_back("LooseMuVeto");
333     fCutLabels.push_back("ZVeto");
334     fCutLabels.push_back("ConversionVeto");
335     fCutLabels.push_back("MET");
336     fCutLabels.push_back("1Jet");
337     fCutLabels.push_back("2Jet");
338     fCutLabels.push_back("3Jet");
339     fCutLabels.push_back("4Jet");
340     fCutLabels.push_back("2Jet1b");
341     fCutLabels.push_back("2Jet2b");
342     }
343     hcutflow = new TH1D("cutflow","cut flow", fCutLabels.size(), 0.5, fCutLabels.size() +0.5 );
344    
345     for ( vector<string>::const_iterator ivec= fCutLabels.begin(); ivec!=fCutLabels.end(); ++ivec)
346     {
347     cutmap[ *ivec ] = 0;
348     }
349    
350     double pu_weights[35] = {
351     0.0255506,
352     0.251923,
353     0.549605,
354     0.924918,
355     1.25977,
356     1.48142,
357     1.57923,
358     1.57799,
359     1.52152,
360     1.4414,
361     1.35889,
362     1.2841,
363     1.21927,
364     1.16125,
365     1.11244,
366     1.06446,
367     1.01666,
368     0.966989,
369     0.913378,
370     0.85774,
371     0.799258,
372     0.734225,
373     0.670242,
374     0.607641,
375     0.54542,
376     0.484084,
377     0.427491,
378     0.369787,
379     0.321454,
380     0.280706,
381     0.238499,
382     0.198961,
383     0.166742,
384     0.146428,
385     0.224425
386     };
387    
388     fpu_weights_vec.assign( pu_weights, pu_weights + 35 );
389    
390     // For JEC uncertainties
391 algomez 1.2 if (fdoJECunc) fJECunc = new JetCorrectionUncertainty("/uscms/home/algomez/work/CMSSW_4_2_4/src/Yumiceva/TreeAnalyzer/test/GR_R_42_V19_AK5PF_Uncertainty.txt");
392 algomez 1.1
393     }
394    
395     Bool_t Analyzer::Process(Long64_t entry)
396     {
397     // The Process() function is called for each entry in the tree (or possibly
398     // keyed object in the case of PROOF) to be processed. The entry argument
399     // specifies which entry in the currently loaded tree is to be processed.
400     // It can be passed to either Analyzer::GetEntry() or TBranch::GetEntry()
401     // to read either all or the required parts of the data. When processing
402     // keyed objects with PROOF, the object is already loaded and is available
403     // via the fObject pointer.
404     //
405     // This function should contain the "body" of the analysis. It can contain
406     // simple or elaborate selection criteria, run algorithms on the data
407     // of the event and typically fill histograms.
408     //
409     // The processing can be stopped by calling Abort().
410     //
411     // Use fStatus to set the return value of TTree::Process().
412     //
413     // The return value is currently not used.
414    
415     //TString option = GetOption();
416    
417     //if ( entry % 100 == 0 )
418     //cout<< "process entry " << entry << endl;
419    
420     //TString sEntry = Form("%f", float(entry) );
421     // Info("Process",
422     //"entry # %s", Form("%f", float(entry) ) );
423    
424     //fChain->GetTree()->GetEntry(entry);
425     fChain->GetEntry(entry);
426    
427     //if (entry>10) return kTRUE;
428    
429     // event info
430     //cout << "run: " << ntuple->run << " lumi: " << ntuple->lumi << endl;
431    
432     // get collections
433     vector< TopVertexEvent > primaryVertices = ntuple->vertices;
434     vector< TopMuonEvent > muons = ntuple->muons;
435 algomez 1.2 vector< TopElectronEvent > electrons = ntuple->PFelectrons; // use PF electrons (gsf collection is called "electrons")
436 algomez 1.1 vector< TopJetEvent > jets = ntuple->PFjets;
437    
438 algomez 1.2 // USE PF Isolation
439     fMuSelector.UsePFIsolation(true);
440     fEleSelector.UsePFIsolation(true);
441    
442 algomez 1.1 size_t total_pvs = primaryVertices.size();
443     size_t total_muons = muons.size();
444     size_t total_electrons = electrons.size();
445     size_t total_jets = jets.size();
446    
447     float PVz = 0.;
448     TLorentzVector p4muon, p4ele, p4lepton, p4MET;
449     TLorentzVector p4Nu, p4OtherNu;
450     TLorentzVector p4QCDmuon;
451    
452     vector< TLorentzVector > p4jets;
453 algomez 1.6 vector< TLorentzVector > p4Othermuon; // leading muon
454 algomez 1.1
455 algomez 1.2 ////////////////////
456     // GENERATOR
457     ///////////////////
458 algomez 1.6 if (fIsMC && fSample.Contains("4Top") )
459 algomez 1.2 {
460     TLorentzVector p4genLepton;
461     TLorentzVector p4genNu;
462     TLorentzVector p4genb;
463     if (ntuple->gen.bLep_pt>0)
464     {
465     p4genLepton.SetPtEtaPhiE(ntuple->gen.mu_pt,ntuple->gen.mu_eta,ntuple->gen.mu_phi,ntuple->gen.mu_e);
466     p4genNu.SetPtEtaPhiE(ntuple->gen.nu_pt,ntuple->gen.nu_eta,ntuple->gen.nu_phi,ntuple->gen.nu_e);
467     p4genb.SetPtEtaPhiE(ntuple->gen.bLep_pt,ntuple->gen.bLep_eta,ntuple->gen.bLep_phi,ntuple->gen.bLep_e);
468    
469     hjets["gen_deltaR_mub"]->Fill( p4genLepton.DeltaR( p4genb ) );
470     }
471     }
472 algomez 1.1 ////////////////////////////////////
473     // PRIMARY VERTICES
474     ///////////////////////////////////
475    
476 algomez 1.6 for ( size_t ipv=0; ipv < total_pvs; ++ipv)
477     {
478     if (ipv==0) PVz = primaryVertices[ipv].vz;
479     }
480 algomez 1.1
481 algomez 1.6 //hPVs["N"]->Fill( total_pvs );
482 algomez 1.1
483     // calculate PU weight
484     double PUweight = 1.;
485    
486     if (fIsMC) {
487 algomez 1.2
488     Int_t mc_npvminus1 = TMath::Min(ntuple->gen.Bx_minus1,34);
489     Int_t mc_npv0 = TMath::Min(ntuple->gen.Bx_0,34);
490     Int_t mc_npvplus1 = TMath::Min(ntuple->gen.Bx_plus1,34);
491    
492     PUweight = f3Dweight->GetBinContent(mc_npvminus1,mc_npv0,mc_npvplus1);
493    
494     //int iibin = 0;
495     //for ( vector<double>::iterator ivec = fpu_weights_vec.begin(); ivec != fpu_weights_vec.end(); ++ivec )
496     // {
497     // int mc_npvs = ntuple->gen.Bx_0; // in-time pile up
498     ////int mc_npvs = (int)total_pvs;
499     //if ( mc_npvs >= iibin+1 ) PUweight = *ivec; // use the last weight for last bin
500     //if ( ( iibin <= mc_npvs ) && ( mc_npvs < iibin + 1 ) ) PUweight = *ivec;
501     //iibin++;
502     //}
503 algomez 1.1 }
504    
505 algomez 1.5 // For 4 tops
506     if (fIsMC && fSample.Contains("4Top") ) {
507     PUweight = 1;
508     }
509    
510 algomez 1.1 hPVs["Nreweight"]->Fill( total_pvs, PUweight );
511    
512 algomez 1.2 /////////////
513     // HLT scale factor for MC
514     ////////////
515     double SF_hlt = 1.;
516     if (fIsMC) SF_hlt = 0.966;
517    
518 algomez 1.6 PUweight = PUweight*SF_hlt; // LETS INCLUDE THE TRIGGER SF INTO THE PU WEIGHTS
519 algomez 1.2
520 algomez 1.1 cutmap["Processed"] += PUweight;
521    
522     if (fVerbose) cout << "done pv" << endl;
523 algomez 1.6
524 algomez 1.1 //////////////////////////////////
525     // MUONS
526     //////////////////////////////////
527 algomez 1.2 int ngoodIDmuons = 0;
528 algomez 1.1 int nloosemuons = 0;
529     int ntightmuons = 0;
530     int nqcdmuons = 0;
531    
532     double RelIso = -1.;
533     double deltaR = -1.;
534     double QCDRelIso = -1.;
535     double QCDdeltaR = -1;
536    
537 algomez 1.2
538 algomez 1.1 for ( size_t imu=0; imu < total_muons; ++imu) {
539    
540 algomez 1.6 TopMuonEvent muon = muons[imu];
541 algomez 1.2
542 algomez 1.6 h1test->Fill( muon.pt );
543     //hmuons["N_muons"]->Fill( total_muons );
544     hmuons["pt_cut1"]->Fill( muon.pt, PUweight );
545    
546     if ( fMuSelector.MuonID( muon, PVz ) ) ngoodIDmuons++;
547    
548     // select only good muons
549     if ( fMuSelector.MuonLoose( muon ) ) {
550     nloosemuons++;
551    
552     if ( fMuSelector.MuonTight( muon, PVz) ) {
553     //hmuons["N_tisomuons"]->Fill( nloosemuons );
554     //hmuons["charge_tiso"]->Fill( muon.charge, PUweight );
555     hmuons["pt_cut2"]->Fill( muon.pt, PUweight );
556     }
557    
558     if ( fMuSelector.MuonTightDeltaR( muon, PVz, jets) ) {
559     ntightmuons++;
560     deltaR = fMuSelector.GetDeltaR();
561     }
562 algomez 1.1
563 algomez 1.6 p4muon.SetPtEtaPhiE( muon.pt, muon.eta, muon.phi, muon.e );
564     RelIso = muon.pfreliso; //muon.reliso03;
565 algomez 1.1
566 algomez 1.6 p4Othermuon.push_back( p4muon ); // for leading muon
567 algomez 1.1 }
568 algomez 1.4
569 algomez 1.6 // check muon in QCD control region
570     if ( fMuSelector.MuonRelax02IsoQCD( muon, PVz, jets ) ) {
571     nqcdmuons++;
572     // keep the leading muon for selection
573     if (nqcdmuons==1) {
574     p4QCDmuon.SetPtEtaPhiE( muon.pt, muon.eta, muon.phi, muon.e );
575     QCDRelIso = muon.pfreliso; //muon.reliso03;
576     QCDdeltaR = fMuSelector.GetDeltaR();
577     }
578 algomez 1.1 }
579 algomez 1.6 }
580 algomez 1.1
581 algomez 1.6 if (fVerbose) cout << "done muons" << endl;
582 algomez 1.1
583     //////////////////////////////////
584     // ELECTRONS
585     //////////////////////////////////
586     int nlooseelectrons = 0;
587     int ntightelectrons = 0;
588     bool IsConversion = false;
589    
590     for ( size_t iele=0; iele < total_electrons; ++iele) {
591    
592 algomez 1.6 TopElectronEvent electron = electrons[iele];
593 algomez 1.1
594 algomez 1.6 if ( fEleSelector.ElectronLoose(electron) ) nlooseelectrons++;
595 algomez 1.1
596 algomez 1.6 if ( fEleSelector.ElectronTight(electron, PVz ) ) {
597 algomez 1.1
598 algomez 1.6 if (ntightelectrons == 0) {
599     IsConversion = electron.IsConversion;
600     p4ele.SetPtEtaPhiE( electron.pt, electron.eta, electron.phi, electron.e );
601     helectrons["pt_cut2"]->Fill( p4ele.Pt(), PUweight );
602     helectrons["eta_cut2"]->Fill( p4ele.Eta(), PUweight );
603     helectrons["phi_cut2"]->Fill( p4ele.Phi(), PUweight );
604     }
605     ntightelectrons++;
606     }
607 algomez 1.1 }
608     if (fVerbose) cout << "done electron" << endl;
609    
610 algomez 1.6 /////////////////////////////////////
611     // MUON/ELECTRON + JETS
612     /////////////////////////////////////
613 algomez 1.1
614 algomez 1.6 if ( fChannel == 1 ) {
615 algomez 1.1
616 algomez 1.6 if (fdoQCD2SideBand) {
617 algomez 1.1
618 algomez 1.6 if (nqcdmuons == 0) return kTRUE;
619     cutmap["OneIsoMu"] += PUweight;
620 algomez 1.1
621 algomez 1.6 p4lepton = p4QCDmuon;
622     RelIso = QCDRelIso;
623     deltaR = QCDdeltaR;
624     }
625     else {
626 algomez 1.2
627 algomez 1.6 if ( ngoodIDmuons > 0 ) hmuons["Ngood"]->Fill( total_pvs, PUweight);
628     if ( ntightmuons > 0 ) hmuons["Niso"]->Fill( total_pvs, PUweight);
629     if ( ntightmuons != 1 ) return kTRUE;
630     cutmap["OneIsoMu"] += PUweight;
631 algomez 1.2
632 algomez 1.6 if ( nloosemuons > 1 ) return kTRUE;
633     cutmap["LooseMuVeto"] += PUweight;
634 algomez 1.1
635 algomez 1.6 if ( nlooseelectrons > 0 ) return kTRUE;
636     cutmap["ElectronVeto"] += PUweight;
637 algomez 1.1
638 algomez 1.6 p4lepton = p4muon;
639     if (fVerbose) cout << "got a good lepton" << endl;
640     }
641    
642     hmuons["pt"]->Fill( p4lepton.Pt(), PUweight );
643     hmuons["eta"]->Fill( p4lepton.Eta(), PUweight );
644     hmuons["phi"]->Fill( p4lepton.Phi(), PUweight );
645     hmuons["reliso"]->Fill( RelIso, PUweight );
646     //hmuons["deltaR"]->Fill( deltaR, PUweight );
647 algomez 1.1
648 algomez 1.6 }
649 algomez 1.1 else // electron+jets
650     {
651     // pending ...
652     }
653    
654     if (fVerbose) cout << "done lepton selection " << endl;
655    
656     /////////////////////////////////
657 algomez 1.6 // MET + Ht
658 algomez 1.1 /////////////////////////////////
659    
660     p4MET.SetPtEtaPhiE( ntuple->PFMET,
661     0.,
662     ntuple->PFMETphi,
663     ntuple->PFMET );
664    
665     //temporal check using genMET
666     //p4MET.SetPtEtaPhiE( ntuple->gen.MET,
667     // 0.,
668     // ntuple->gen.METPhi,
669     // ntuple->gen.MET );
670    
671     if (fdoQCD1SideBand && p4MET.Et() > 20.) return kTRUE;
672     else if ( p4MET.Et() <= 20. && fdoQCD2SideBand==false ) return kTRUE;
673 algomez 1.6 if (fVerbose) cout << "pass MET cut" << endl;
674 algomez 1.1
675    
676     //cutmap["MET"] += PUweight;
677 algomez 1.6 //hMET["MET"]->Fill( p4MET.Pt(), PUweight );
678     //hMET["phi"]->Fill( p4MET.Phi(), PUweight );
679     //hMET["Ht0"]->Fill( ntuple->PFHt, PUweight );
680     //hMET["Htlep"]->Fill( ntuple->PFHt + p4lepton.Pt(), PUweight );
681    
682 algomez 1.1
683     double Wpt = p4lepton.Pt() + p4MET.Pt();
684     double Wpx = p4lepton.Px() + p4MET.Px();
685     double Wpy = p4lepton.Py() + p4MET.Py();
686     double WMt = sqrt(Wpt*Wpt-Wpx*Wpx-Wpy*Wpy);
687 algomez 1.3
688     //if ( WMt <= 40. ) return kTRUE;
689 algomez 1.1 cutmap["MET"] += PUweight;
690    
691     if (fVerbose) cout << "pass W MT cut " << endl;
692    
693     /////////////////////////////////
694     // estimate Pz of neutrino
695     ////////////////////////////////
696    
697     fzCalculator.SetMET(p4MET);
698     fzCalculator.SetLepton(p4lepton);
699 algomez 1.6 if (fChannel==2) {
700     fzCalculator.SetLeptonType("electron");
701     }
702 algomez 1.1
703     double pzNu = fzCalculator.Calculate();
704     p4Nu = TLorentzVector();
705     p4OtherNu = TLorentzVector();
706    
707     p4Nu.SetPxPyPzE(p4MET.Px(), p4MET.Py(), pzNu, sqrt(p4MET.Px()*p4MET.Px()+p4MET.Py()*p4MET.Py()+pzNu*pzNu));
708     //print "pzNu = " +str(pzNu)
709     //print "p4Nu =("+str(p4Nu.Px())+","+str(p4Nu.Py())+","+str(p4Nu.Pz())+","+str(p4Nu.E())
710     double pzOtherNu = fzCalculator.getOther();
711     p4OtherNu.SetPxPyPzE( p4MET.Px(), p4MET.Py(),pzOtherNu,sqrt(p4MET.Px()*p4MET.Px()+p4MET.Py()*p4MET.Py()+pzOtherNu*pzOtherNu));
712    
713 algomez 1.6 //double WmassNoPt = (p4Nu+p4lepton).M();
714 algomez 1.1
715 algomez 1.6 if ( fzCalculator.IsComplex() ) {
716     double ptNu1 = fzCalculator.getPtneutrino(1);
717     double ptNu2 = fzCalculator.getPtneutrino(2);
718     TLorentzVector p4Nu1tmp;
719     TLorentzVector p4Nu2tmp;
720 algomez 1.1
721 algomez 1.6 p4Nu1tmp.SetPxPyPzE( ptNu1*p4MET.Px()/p4MET.Pt(), ptNu1*p4MET.Py()/p4MET.Pt(), pzNu, sqrt(ptNu1*ptNu1+pzNu*pzNu));
722     p4Nu2tmp.SetPxPyPzE( ptNu2*p4MET.Px()/p4MET.Pt(), ptNu2*p4MET.Py()/p4MET.Pt(), pzNu, sqrt(ptNu2*ptNu2+pzNu*pzNu));
723 algomez 1.1
724 algomez 1.6 TLorentzVector Wtmp;
725     Wtmp = p4lepton + p4Nu1tmp;
726     double Wm1 = 0;
727     double Wm2 = 0;
728     Wm1 = Wtmp.M();
729     Wtmp = p4lepton + p4Nu2tmp;
730     Wm2 = Wtmp.M();
731     if ( fabs( Wm1 - 80.) < fabs( Wm2 - 80.) ) p4Nu = p4Nu1tmp;
732     else p4Nu = p4Nu2tmp;
733 algomez 1.1
734 algomez 1.6 p4OtherNu = p4Nu; // since we chose the real part, the two solutions are the same.
735     }
736 algomez 1.1
737    
738     hMET["PzNu"]->Fill(pzNu, PUweight ); //change this to 2d with two sol and as a function of jets
739    
740     TLorentzVector p4LepW = p4lepton + p4Nu;
741     TLorentzVector p4OtherLepW = p4lepton + p4OtherNu;
742    
743     //hMET["LepWmass"]->Fill(p4LepW.M(), PUweight );
744     //if ( fzCalculator.IsComplex() ) hMET["LepWmassComplex"]->Fill( p4LepW.M(), PUweight );
745    
746    
747     /////////////////////////////////
748     // JETS
749     ////////////////////////////////
750    
751 algomez 1.6 //JetCombinatorics myCombi = JetCombinatoric();
752 algomez 1.1
753     int njets = 0;
754     map< string, vector<float> > bdisc;
755     map< string, vector<bool> > isTagb;
756     vector<int> listflavor;
757    
758 algomez 1.6 for ( size_t ijet=0; ijet < total_jets; ++ijet) {
759 algomez 1.1
760 algomez 1.6 TopJetEvent jet = jets[ijet];
761     double SF_JEC = 1.;
762     if (fdoJECunc){
763     fJECunc->setJetEta( jet.eta);
764     fJECunc->setJetPt( jet.pt);
765     double jec_unc = fJECunc->getUncertainty(true);
766     if (fVerbose) cout << "JEC uncertainty is " << jec_unc << endl;
767     if (fdoJECup) SF_JEC = 1.+jec_unc;
768     else SF_JEC = 1.-jec_unc;
769 algomez 1.1 }
770    
771    
772 algomez 1.6 if ( SF_JEC*jet.pt > 35. && fabs(jet.eta) < 2.4 && SF_JEC*jets[0].pt > 120. ) {
773    
774     if (fVerbose) cout << " jet pt " << SF_JEC*jet.pt << endl;
775    
776     //hjets["pt"]->Fill( jet.pt, PUweight );
777     //hjets["eta"]->Fill( jet.eta, PUweight );
778     //hjets["phi"]->Fill( jet.phi, PUweight );
779    
780     TLorentzVector tmpjet;
781     tmpjet.SetPtEtaPhiE(SF_JEC*jet.pt, jet.eta, jet.phi, SF_JEC*jet.e);
782     p4jets.push_back( tmpjet);
783     listflavor.push_back( jet.mc.flavor );
784    
785     if (fVerbose) {
786     cout << "done storing njets " << njets << endl;
787     cout << " bdisc " << jet.btag_TCHP << endl;
788     cout << " bdisc " << jet.btag_CSV << endl;
789     }
790    
791     /*/ store discriminators
792     bdisc["TCHP"].push_back( jet.btag_TCHP );
793     bdisc["CSV"].push_back( jet.btag_CSV );
794     if (fVerbose) cout << "store bdisc" << endl;
795     // TCHPL cut at 1.19
796     // check TCHPM cut at 1.93
797     if ( jet.btag_TCHP > 1.93 ) isTagb["TCHPM"].push_back(true);
798     else isTagb["TCHPM"].push_back(false);
799     if (fVerbose) cout << "done tchpl" << endl;
800     // check SSVHEM cut at 1.74
801     if ( jet.btag_CSV > 0.679) isTagb["CSVM"].push_back(true);
802     else isTagb["CSVM"].push_back(false);
803     if (fVerbose) cout << "done csvm" << endl;
804     // CSVM cut at 0.679
805     */
806     njets++;
807     }
808     }
809 algomez 1.1
810 algomez 1.6 if (njets>0) {
811     hM["WMt"]->Fill( WMt, PUweight );
812     hjets["Njets"]->Fill( njets, PUweight );
813     }
814 algomez 1.1
815     if (fVerbose) cout << "done jets" << endl;
816    
817     if (njets > 0 ) cutmap["1Jet"] += PUweight;
818     if (njets > 1 ) cutmap["2Jet"] += PUweight;
819     if (njets > 2 ) cutmap["3Jet"] += PUweight;
820     if (njets > 3 ) cutmap["4Jet"] += PUweight;
821    
822 algomez 1.6 if (njets >= 2) {
823 algomez 1.1
824 algomez 1.6 // Ht Calculation
825     double Ht = 0;
826     double deltaRjj = 999.;
827     for ( size_t kk=0; kk < p4jets.size(); ++kk) {
828     Ht += p4jets[kk].Pt();
829    
830     for ( vector< TopJetEvent>::iterator ijet=jets.begin(); ijet != jets.end(); ++ijet) {
831    
832     TopJetEvent jet = *ijet;
833     TLorentzVector tmpp4Jet;
834     tmpp4Jet.SetPtEtaPhiE(jet.pt, jet.eta, jet.phi, jet.e );
835     double tmpdeltaR = p4jets[kk].DeltaR(tmpp4Jet);
836     if ( tmpdeltaR < 0.3 ) continue;
837     if ( tmpdeltaR < deltaRjj ) deltaRjj = tmpdeltaR;
838     }
839 algomez 1.1 }
840    
841 algomez 1.4
842 algomez 1.6 hMET["Ht"]->Fill( Ht, PUweight );
843     hMET["Ht0"]->Fill( ntuple->PFHt, PUweight );
844     hjets["Njets_cut0"]->Fill( njets, PUweight );
845     hjets["1st_pt_cut0"]->Fill( p4jets[0].Pt(), PUweight );
846     hjets["1st_eta_cut0"]->Fill( p4jets[0].Eta(), PUweight );
847     hjets["2nd_pt_cut0"]->Fill( p4jets[1].Pt(), PUweight );
848     hjets["2nd_eta_cut0"]->Fill( p4jets[1].Eta(), PUweight );
849     hjets["3rd_pt_cut0"]->Fill( p4jets[2].Pt(), PUweight );
850     hjets["4th_pt_cut0"]->Fill( p4jets[3].Pt(), PUweight );
851     hjets["5th_pt_cut0"]->Fill( p4jets[4].Pt(), PUweight );
852     hjets["6th_pt_cut0"]->Fill( p4jets[5].Pt(), PUweight );
853     hjets["7th_pt_cut0"]->Fill( p4jets[6].Pt(), PUweight );
854     hmuons["N_cut0"]->Fill( total_muons, PUweight );
855     hmuons["Nelectrons_cut0"]->Fill( nlooseelectrons, PUweight );
856     hMET["MET_cut0"]->Fill( p4MET.Pt(), PUweight );
857     hmuons["deltaR_cut0"]->Fill( deltaR, PUweight );
858     hjets["Dijet_deltaR_cut0"]->Fill( deltaRjj, PUweight );
859    
860     bool passcut = true;
861     if ( Ht <= 300. ) passcut = false;
862    
863     if (passcut) {
864     hjets["Njets_cut1"]->Fill( njets, PUweight );
865     hjets["1st_pt"]->Fill( p4jets[0].Pt(), PUweight );
866     hjets["1st_eta"]->Fill( p4jets[0].Eta(), PUweight );
867     hjets["2nd_pt"]->Fill( p4jets[1].Pt(), PUweight );
868     hjets["2nd_eta"]->Fill( p4jets[1].Eta(), PUweight );
869     hjets["3rd_pt"]->Fill( p4jets[2].Pt(), PUweight );
870     hjets["4th_pt"]->Fill( p4jets[3].Pt(), PUweight );
871     hjets["5th_pt"]->Fill( p4jets[4].Pt(), PUweight );
872     hjets["6th_pt"]->Fill( p4jets[5].Pt(), PUweight );
873     hjets["7th_pt"]->Fill( p4jets[6].Pt(), PUweight );
874     hmuons["N"]->Fill( total_muons, PUweight );
875     hmuons["Nelectrons"]->Fill( nlooseelectrons, PUweight );
876     hMET["MET"]->Fill( p4MET.Pt(), PUweight );
877     hmuons["deltaR"]->Fill( deltaR, PUweight );
878     hjets["Dijet_deltaR"]->Fill( deltaRjj, PUweight );
879 algomez 1.4 }
880 algomez 1.6 }
881 algomez 1.1
882     if (fVerbose) cout << "done analysis" << endl;
883     return kTRUE;
884     }
885    
886     void Analyzer::SlaveTerminate()
887     {
888     // The SlaveTerminate() function is called after all entries or objects
889     // have been processed. When running with PROOF SlaveTerminate() is called
890     // on each slave server.
891    
892     // fill cutflow histogram
893    
894     int ibin = 1;
895     for ( vector<string>::const_iterator ivec= fCutLabels.begin(); ivec != fCutLabels.end(); ++ivec )
896     // for ( map<string, int >::const_iterator imap=cutmap.begin(); imap!=cutmap.end(); ++imap )
897     {
898     hcutflow->SetBinContent( ibin, cutmap[ *ivec ] );
899     ibin++;
900     }
901     // Write the ntuple to the file
902     if (fFile) {
903     Bool_t cleanup = kFALSE;
904     TDirectory *savedir = gDirectory;
905     if ( h1test->GetEntries() > 0) {
906     fFile->cd();
907     h1test->Write();
908     hcutflow->Write();
909 algomez 1.3 //h2_pt_Wprime->Write();
910 algomez 1.1 fFile->mkdir("muons");
911     fFile->cd("muons");
912     for ( map<string,TH1* >::const_iterator imap=hmuons.begin(); imap!=hmuons.end(); ++imap )
913     {
914     TH1 *temp = imap->second;
915     if ( temp->GetEntries() > 0 )
916     temp->Write();
917     //else cout << "Warning: empty histogram " << temp->GetName() << " will not be written to file." << endl;
918     }
919     fFile->cd();
920     fFile->mkdir("PVs");
921     fFile->cd("PVs");
922     for ( map<string,TH1* >::const_iterator imap=hPVs.begin(); imap!=hPVs.end(); ++imap )
923     {
924     TH1 *temp = imap->second;
925     if ( temp->GetEntries() > 0 )
926     temp->Write();
927     //else cout << "Warning: empty histogram " << temp->GetName() << " will not be written to file." << endl;
928     }
929     fFile->cd();
930     fFile->mkdir("jets");
931     fFile->cd("jets");
932     for ( map<string,TH1* >::const_iterator imap=hjets.begin(); imap!=hjets.end(); ++imap )
933     {
934     TH1 *temp = imap->second;
935     if ( temp->GetEntries() > 0 )
936     temp->Write();
937     }
938     fFile->cd();
939     fFile->mkdir("mass");
940     fFile->cd("mass");
941     for ( map<string,TH1* >::const_iterator imap=hM.begin(); imap!=hM.end(); ++imap )
942     {
943     TH1 *temp = imap->second;
944     if ( temp->GetEntries() > 0 )
945     temp->Write();
946     }
947     fFile->cd();
948     fFile->mkdir("MET");
949     fFile->cd("MET");
950     for ( map<string,TH1* >::const_iterator imap=hMET.begin(); imap!=hMET.end(); ++imap )
951     {
952     TH1 *temp = imap->second;
953     if ( temp->GetEntries() > 0 )
954     temp->Write();
955     }
956    
957     fFile->cd();
958    
959     fProofFile->Print();
960     fOutput->Add(fProofFile);
961     } else {
962     cleanup = kTRUE;
963     }
964    
965    
966     h1test->SetDirectory(0);
967     hcutflow->SetDirectory(0);
968     gDirectory = savedir;
969     fFile->Close();
970     // Cleanup, if needed
971     if (cleanup) {
972     TUrl uf(*(fFile->GetEndpointUrl()));
973     SafeDelete(fFile);
974     gSystem->Unlink(uf.GetFile());
975     SafeDelete(fProofFile);
976     }
977     }
978    
979     }
980    
981     void Analyzer::Terminate()
982     {
983     // The Terminate() function is the last function to be called during
984     // a query. It always runs on the client, it can be used to present
985     // the results graphically or save the results to file.
986    
987     Info("Terminate","Analyzer done.");
988     }