ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/algomez/FourTop/TreeAnalyzer/test/Analyzer.C
Revision: 1.2
Committed: Thu Dec 8 20:18:14 2011 UTC (13 years, 5 months ago) by algomez
Content type: text/plain
Branch: MAIN
Changes since 1.1: +219 -171 lines
Log Message:
adding some new plots

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     TString weightfilename = "/uscms/home/algomez/work/CMSSW_4_2_4/src/Yumiceva/TreeAnalyzer/test/Weight3D.root";
155     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.2 hmuons["N_muons"] = new TH1F("N_muons"+hname,"Number of Muons",6,-0.5,5.5);
164 algomez 1.1 hmuons["pt_1jet"] = new TH1F("muon_pt_1jet"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
165     hmuons["pt_2jet"] = new TH1F("muon_pt_2jet"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
166 algomez 1.2 hmuons["pt_fin"] = new TH1F("muon_pt_fin"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
167 algomez 1.1 hmuons["eta"] = new TH1F("muon_eta"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
168     hmuons["eta_1jet"] = new TH1F("muon_eta_1jet"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
169     hmuons["eta_2jet"] = new TH1F("muon_eta_2jet"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
170 algomez 1.2 hmuons["eta_fin"] = new TH1F("muon_eta_fin"+hname,"#eta^{#mu}", 20, -2.1, 2.1);
171 algomez 1.1 hmuons["phi"] = new TH1F("muon_phi"+hname,"#phi^{#mu}", 30, -3.15, 3.15);
172     hmuons["phi_1jet"] = new TH1F("muon_phi_1jet"+hname,"#phi^{#mu}", 30, -3.15, 3.15);
173     hmuons["phi_2jet"] = new TH1F("muon_phi_2jet"+hname,"#phi^{#mu}", 30, -3.15, 3.15);
174 algomez 1.2 hmuons["phi_fin"] = new TH1F("muon_phi_fin"+hname,"#phi^{#mu}", 30, -3.15, 3.15);
175 algomez 1.1 hmuons["reliso"] = new TH1F("muon_reliso"+hname,"Relative Isolation", 40, 0, 0.2);
176     hmuons["reliso_1jet"] = new TH1F("muon_reliso_1jet"+hname,"Relative Isolation", 40, 0, 0.2);
177     hmuons["reliso_2jet"] = new TH1F("muon_reliso_2jet"+hname,"Relative Isolation", 40, 0, 0.2);
178 algomez 1.2 hmuons["reliso_fin"] = new TH1F("muon_reliso_fin"+hname,"Relative Isolation", 40, 0, 0.2);
179 algomez 1.1 hmuons["deltaR_cut0"] = new TH1F("deltaR_cut0"+hname,"#DeltaR(#mu,jet)",80, 0, 4);
180     hmuons["deltaR"] = new TH1F("deltaR"+hname,"#DeltaR(#mu,jet)", 80, 0, 4);
181     hmuons["d0_cut1"] = new TH1F("d0_cut1"+hname,"#mu Impact Parameter [cm]",20,-0.1,0.1);
182     hmuons["pt_cut1"] = new TH1F("muon_pt_cut1"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
183     hmuons["pt_cut2"] = new TH1F("muon_pt_cut2"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
184     hmuons["pt_cut3"] = new TH1F("muon_pt_cut3"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
185     hmuons["pt"] = new TH1F("muon_pt"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 500);
186 algomez 1.2 //hmuons["pt"] = new TH1F("muon_pt"+hname,"p_{T}^{#mu} [GeV/c]", 50, 0, 800); // for 1 TeV
187 algomez 1.1 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     hMET["MET"] = new TH1F("MET"+hname,"Missing Transverse Energy [GeV]", 50, 0, 300);
224 algomez 1.2 //hMET["MET"] = new TH1F("MET"+hname,"Missing Transverse Energy [GeV]", 50, 0, 600); // for 1 TeV
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     hMET["Ht"] = new TH1F("Ht"+hname,"H_{T} [GeV]", 50, 0, 1000);
231     hMET["Htlep"] = new TH1F("Htlep"+hname,"H_{T,lep} [GeV]", 50, 0, 1000);
232     hMET["PzNu"] = new TH1F("PzNu"+hname,"p_{z} #nu [GeV]", 40, -300,300);
233     hMET["EtaNu"] = new TH1F("EtaNu"+hname,"#eta",50,-2.2,2.2);
234     hMET["LepWmass"] = new TH1F("LepWmass"+hname,"W#rightarrow#mu#nu Mass [GeV/c^{2}]",20, 0, 150);
235     hMET["LepWmass_topcut"] = new TH1F("LepWmass_topcut"+hname,"W#rightarrow#mu#nu Mass [GeV/c^{2}]",20, 0, 150);
236     hMET["LepWmassNoPt"]=new TH1F("LepWmassNoPt"+hname,"W#rightarrow#mu#nu Mass [GeV/c^{2}]",20, 0, 150);
237     hMET["deltaPhi"] = new TH1F("deltaPhi"+hname,"#Delta #phi(#mu,MET)",50, -3.15, 3.15);
238    
239     hM["WMt"] = new TH1F("Mt"+hname,"M_{T}(W) [GeV/c^{2}]", 50, 0, 300);
240     hM["WMt_2jet"] = new TH1F("Mt_2jet"+hname,"M_{T}(W) [GeV/c^{2}]", 50, 0, 300);
241 algomez 1.2 hM["WMt_fin"] = new TH1F("Mt_fin"+hname,"M_{T}(W) [GeV/c^{2}]", 50, 0, 300);
242 algomez 1.1 hM["dijet"] = new TH1F("dijet"+hname,"(jj) mass [GeV/c^{2}]", 50, 0, 1000);
243     hM["trijet"] = new TH1F("trijet"+hname,"(jjj) mass [GeV/c^{2}]", 50, 0, 1000);
244     hM["top_1btag"] = new TH1F("top_1btag"+hname,"top mass [GeV/c^{2}]",50,0,500);
245     hM["top_2btag"] = new TH1F("top_2btag"+hname,"top mass [GeV/c^{2}]",50,0,500);
246     hM["Wprime"] = new TH1F("Wprime"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
247     hM["Wprime_0btag"] = new TH1F("Wprime_0btag"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
248     hM["Wprime_0btag_light"] = new TH1F("Wprime_0btag_light"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
249     hM["Wprime_0btag_bb"] = new TH1F("Wprime_0btag_bb"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
250     hM["Wprime_0btag_cc"] = new TH1F("Wprime_0btag_cc"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
251     hM["Wprime_1btag"] = new TH1F("Wprime_1btag"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
252     hM["Wprime_1btag_systUp"] = new TH1F("Wprime_1btag_systUp"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
253     hM["Wprime_1btag_systDown"] = new TH1F("Wprime_1btag_systDown"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
254     if (fSample == "WJets")
255     {
256     hM["Wprime_0btag_light"] = new TH1F("Wprime_0btag_light"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
257     hM["Wprime_0btag_bb"] = new TH1F("Wprime_0btag_bb"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
258     hM["Wprime_0btag_cc"] = new TH1F("Wprime_0btag_cc"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
259     hM["Wprime_1btag_light"] = new TH1F("Wprime_1btag_light"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
260     hM["Wprime_1btag_bb"] = new TH1F("Wprime_1btag_bb"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
261     hM["Wprime_1btag_cc"] = new TH1F("Wprime_1btag_cc"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
262     }
263     hM["Wprime_1onlybtag"] = new TH1F("Wprime_1onlybtag"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
264     hM["Wprime_2btag"] = new TH1F("Wprime_2btag"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
265     hM["Wprime_2btag_systUp"] = new TH1F("Wprime_2btag_systUp"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
266     hM["Wprime_2btag_systDown"] = new TH1F("Wprime_2btag_systDown"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
267     if ( fSample.Contains("Wprime") )
268     {
269     hM["Wprime_0btag_MChad"] = new TH1F("Wprime_0btag_MChad"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
270     hM["Wprime_0btag_MCsemimu"] = new TH1F("Wprime_0btag_MCsemimu"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
271     hM["Wprime_0btag_MCsemie"] = new TH1F("Wprime_0btag_MCsemie"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
272     hM["Wprime_1btag_MChad"] = new TH1F("Wprime_1btag_MChad"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
273     hM["Wprime_1btag_MCsemimu"] = new TH1F("Wprime_1btag_MCsemimu"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
274     hM["Wprime_1btag_MCsemie"] = new TH1F("Wprime_1btag_MCsemie"+hname,"W' invariant mass [GeV/c^{2}]", 70, 100, 3000);
275     }
276    
277     hjets["pt"] = new TH1F("jet_pt"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
278 algomez 1.2 hjets["pt_b_mc"] = new TH1F("jet_pt_b_mc"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
279     hjets["pt_c_mc"] = new TH1F("jet_pt_c_mc"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
280     hjets["pt_l_mc"] = new TH1F("jet_pt_l_mc"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
281 algomez 1.1 hjets["pt_btag"] = new TH1F("jet_pt_btag"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
282     hjets["pt_btag_b"] = new TH1F("jet_pt_btag_b"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
283     hjets["pt_btag_c"] = new TH1F("jet_pt_btag_c"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
284     hjets["pt_btag_l"] = new TH1F("jet_pt_btag_l"+hname,"jet p_{T} [GeV/c]", 60, 30, 800);
285     hjets["1st_pt"] = new TH1F("jet1_pt"+hname,"Leading jet p_{T} [GeV/c]",60, 30, 800);
286 algomez 1.2 //hjets["1st_pt"] = new TH1F("jet1_pt"+hname,"Leading jet p_{T} [GeV/c]",60, 30, 1000); // for 1 TeV
287 algomez 1.1 hjets["2nd_pt"] = new TH1F("jet2_pt"+hname,"2nd jet p_{T} [GeV/c]",60, 30, 800);
288 algomez 1.2 //hjets["2nd_pt"] = new TH1F("jet2_pt"+hname,"2nd jet p_{T} [GeV/c]",60, 30, 1000); // for 1 TeV
289 algomez 1.1 hjets["3rd_pt"] = new TH1F("jet3_pt"+hname,"3rd jet p_{T} [GeV/c]",60, 30, 800);
290     hjets["4th_pt"] = new TH1F("jet4_pt"+hname,"4th jet p_{T} [GeV/c]",60, 30, 800);
291 algomez 1.2 hjets["1st_pt_fin"] = new TH1F("jet1_pt_fin"+hname,"Leading jet p_{T} [GeV/c]",60, 30, 800);
292     hjets["2nd_pt_fin"] = new TH1F("jet2_pt_fin"+hname,"2nd jet p_{T} [GeV/c]",60, 30, 800);
293 algomez 1.1 hjets["eta"] = new TH1F("jet_eta"+hname,"jet #eta",50, -2.4, 2.4);
294     hjets["1st_eta"] = new TH1F("jet1_eta"+hname,"Leading jet #eta",50, -2.4, 2.4);
295     hjets["2nd_eta"] = new TH1F("jet2_eta"+hname,"2nd jet #eta",50, -2.4, 2.4);
296 algomez 1.2 hjets["1st_eta_fin"] = new TH1F("jet1_eta_fin"+hname,"Leading jet #eta",50, -2.4, 2.4);
297     hjets["2nd_eta_fin"] = new TH1F("jet2_eta_fin"+hname,"2nd jet #eta",50, -2.4, 2.4);
298 algomez 1.1 hjets["phi"] = new TH1F("jet_phi"+hname,"jet #phi",50, 0, 3.15);
299 algomez 1.2 hjets["Njets"] = new TH1F("Njets"+hname,"jet multiplicity",13,0.5,12.5);
300 algomez 1.1 hjets["Njets_1tag"] = new TH1F("Njets_1tag"+hname,"jet multiplicity",6,0.5,6.5);
301     hjets["Njets_2tag"] = new TH1F("Njets_2tag"+hname,"jet multiplicity",6,0.5,6.5);
302     hjets["Nbtags_TCHPM"] = new TH1F("Nbjets_TCHPM_N2j"+hname,"Tagged b-jets",3,-0.5,2.5);
303 algomez 1.2 hjets["Nbtags_CSVM"] = new TH1F("Nbjets_CSVM_N2j"+hname,"Tagged b-jets",3,-0.5,2.5);
304 algomez 1.1 hjets["Dijet_deltaPhi"] = new TH1F("Dijet_deltaPhi"+hname,"#Delta #phi(j,j)",30,-3.15,3.15);
305     hjets["tb_deltaPhi"] = new TH1F("tb_deltaPhi"+hname,"#Delta #phi(t,b)",30,0.,3.15);
306 algomez 1.2 hjets["tb_deltaEta"] = new TH1F("tb_deltaEta"+hname,"#Delta #eta(t,b)",30,-5,5);
307 algomez 1.1 h2_pt_Wprime = new TH2F("toppt_vs_Wprime"+hname,"top p_{T} vs W' mass",60,0,1500,70, 100, 3000);
308 algomez 1.2 hjets["pt_Wprime"] = new TH1F("pt_Wprime"+hname,"W' p_{T} [GeV]",50,0,200);
309     hjets["pt_top"] = new TH1F("pt_top"+hname,"top p_{T} [GeV]",60,0,1500);
310     hjets["pt_b"] = new TH1F("pt_b"+hname,"b-jet p_{T} [GeV]",60,0,1500);
311     hjets["gen_deltaR_mub"] = new TH1F("gen_deltaR_mub"+hname,"#Delta R(#mu,b)",40,0,4);
312 algomez 1.1
313     map<string,TH1* > allhistos = hmuons;
314     allhistos.insert( helectrons.begin(), helectrons.end() );
315     allhistos.insert( hMET.begin(), hMET.end() );
316     allhistos.insert( hM.begin(), hM.end() );
317     allhistos.insert( hjets.begin(), hjets.end() );
318    
319     for ( map<string,TH1* >::const_iterator imap=allhistos.begin(); imap!=allhistos.end(); ++imap )
320     {
321     TH1 *temp = imap->second;
322     temp->Sumw2();
323     temp->SetXTitle( temp->GetTitle() );
324     }
325    
326     // cut flow
327     if (fChannel==1)
328     { //muon +jets
329     fCutLabels.push_back("Processed");
330     fCutLabels.push_back("OneIsoMu");
331     fCutLabels.push_back("LooseMuVeto");
332     fCutLabels.push_back("ElectronVeto");
333     fCutLabels.push_back("MET");
334     fCutLabels.push_back("1Jet");
335     fCutLabels.push_back("2Jet");
336     fCutLabels.push_back("3Jet");
337     fCutLabels.push_back("4Jet");
338     fCutLabels.push_back("2Jet1b");
339     fCutLabels.push_back("2Jet2b");
340     fCutLabels.push_back("MaxJets");
341     fCutLabels.push_back("phi");
342     fCutLabels.push_back("topmass");
343     }
344     else
345     { //electron+jets
346     fCutLabels.push_back("Processed");
347     fCutLabels.push_back("OneIsoElectron");
348     fCutLabels.push_back("LooseMuVeto");
349     fCutLabels.push_back("ZVeto");
350     fCutLabels.push_back("ConversionVeto");
351     fCutLabels.push_back("MET");
352     fCutLabels.push_back("1Jet");
353     fCutLabels.push_back("2Jet");
354     fCutLabels.push_back("3Jet");
355     fCutLabels.push_back("4Jet");
356     fCutLabels.push_back("2Jet1b");
357     fCutLabels.push_back("2Jet2b");
358     }
359     hcutflow = new TH1D("cutflow","cut flow", fCutLabels.size(), 0.5, fCutLabels.size() +0.5 );
360    
361     for ( vector<string>::const_iterator ivec= fCutLabels.begin(); ivec!=fCutLabels.end(); ++ivec)
362     {
363     cutmap[ *ivec ] = 0;
364     }
365    
366     double pu_weights[35] = {
367     0.0255506,
368     0.251923,
369     0.549605,
370     0.924918,
371     1.25977,
372     1.48142,
373     1.57923,
374     1.57799,
375     1.52152,
376     1.4414,
377     1.35889,
378     1.2841,
379     1.21927,
380     1.16125,
381     1.11244,
382     1.06446,
383     1.01666,
384     0.966989,
385     0.913378,
386     0.85774,
387     0.799258,
388     0.734225,
389     0.670242,
390     0.607641,
391     0.54542,
392     0.484084,
393     0.427491,
394     0.369787,
395     0.321454,
396     0.280706,
397     0.238499,
398     0.198961,
399     0.166742,
400     0.146428,
401     0.224425
402     };
403    
404     fpu_weights_vec.assign( pu_weights, pu_weights + 35 );
405    
406     // For JEC uncertainties
407 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");
408 algomez 1.1
409     }
410    
411     Bool_t Analyzer::Process(Long64_t entry)
412     {
413     // The Process() function is called for each entry in the tree (or possibly
414     // keyed object in the case of PROOF) to be processed. The entry argument
415     // specifies which entry in the currently loaded tree is to be processed.
416     // It can be passed to either Analyzer::GetEntry() or TBranch::GetEntry()
417     // to read either all or the required parts of the data. When processing
418     // keyed objects with PROOF, the object is already loaded and is available
419     // via the fObject pointer.
420     //
421     // This function should contain the "body" of the analysis. It can contain
422     // simple or elaborate selection criteria, run algorithms on the data
423     // of the event and typically fill histograms.
424     //
425     // The processing can be stopped by calling Abort().
426     //
427     // Use fStatus to set the return value of TTree::Process().
428     //
429     // The return value is currently not used.
430    
431     //TString option = GetOption();
432    
433     //if ( entry % 100 == 0 )
434     //cout<< "process entry " << entry << endl;
435    
436     //TString sEntry = Form("%f", float(entry) );
437     // Info("Process",
438     //"entry # %s", Form("%f", float(entry) ) );
439    
440     //fChain->GetTree()->GetEntry(entry);
441     fChain->GetEntry(entry);
442    
443     //if (entry>10) return kTRUE;
444    
445     // event info
446     //cout << "run: " << ntuple->run << " lumi: " << ntuple->lumi << endl;
447    
448     // get collections
449     vector< TopVertexEvent > primaryVertices = ntuple->vertices;
450     vector< TopMuonEvent > muons = ntuple->muons;
451 algomez 1.2 vector< TopElectronEvent > electrons = ntuple->PFelectrons; // use PF electrons (gsf collection is called "electrons")
452 algomez 1.1 vector< TopJetEvent > jets = ntuple->PFjets;
453    
454 algomez 1.2 // USE PF Isolation
455     fMuSelector.UsePFIsolation(true);
456     fEleSelector.UsePFIsolation(true);
457    
458 algomez 1.1 size_t total_pvs = primaryVertices.size();
459     size_t total_muons = muons.size();
460     size_t total_electrons = electrons.size();
461     size_t total_jets = jets.size();
462    
463     float PVz = 0.;
464     TLorentzVector p4muon, p4ele, p4lepton, p4MET;
465     TLorentzVector p4Nu, p4OtherNu;
466     TLorentzVector p4QCDmuon;
467    
468     vector< TLorentzVector > p4jets;
469    
470 algomez 1.2 ////////////////////
471     // GENERATOR
472     ///////////////////
473     if (fIsMC && fSample.Contains("Wp") )
474     {
475     TLorentzVector p4genLepton;
476     TLorentzVector p4genNu;
477     TLorentzVector p4genb;
478     if (ntuple->gen.bLep_pt>0)
479     {
480     p4genLepton.SetPtEtaPhiE(ntuple->gen.mu_pt,ntuple->gen.mu_eta,ntuple->gen.mu_phi,ntuple->gen.mu_e);
481     p4genNu.SetPtEtaPhiE(ntuple->gen.nu_pt,ntuple->gen.nu_eta,ntuple->gen.nu_phi,ntuple->gen.nu_e);
482     p4genb.SetPtEtaPhiE(ntuple->gen.bLep_pt,ntuple->gen.bLep_eta,ntuple->gen.bLep_phi,ntuple->gen.bLep_e);
483    
484     hjets["gen_deltaR_mub"]->Fill( p4genLepton.DeltaR( p4genb ) );
485     }
486     }
487 algomez 1.1 ////////////////////////////////////
488     // PRIMARY VERTICES
489     ///////////////////////////////////
490    
491     for ( size_t ipv=0; ipv < total_pvs; ++ipv) {
492    
493     if (ipv==0) PVz = primaryVertices[ipv].vz;
494    
495     }
496    
497     hPVs["N"]->Fill( total_pvs );
498     // calculate PU weight
499     double PUweight = 1.;
500    
501     if (fIsMC) {
502 algomez 1.2
503     Int_t mc_npvminus1 = TMath::Min(ntuple->gen.Bx_minus1,34);
504     Int_t mc_npv0 = TMath::Min(ntuple->gen.Bx_0,34);
505     Int_t mc_npvplus1 = TMath::Min(ntuple->gen.Bx_plus1,34);
506    
507     PUweight = f3Dweight->GetBinContent(mc_npvminus1,mc_npv0,mc_npvplus1);
508    
509     //int iibin = 0;
510     //for ( vector<double>::iterator ivec = fpu_weights_vec.begin(); ivec != fpu_weights_vec.end(); ++ivec )
511     // {
512     // int mc_npvs = ntuple->gen.Bx_0; // in-time pile up
513     ////int mc_npvs = (int)total_pvs;
514     //if ( mc_npvs >= iibin+1 ) PUweight = *ivec; // use the last weight for last bin
515     //if ( ( iibin <= mc_npvs ) && ( mc_npvs < iibin + 1 ) ) PUweight = *ivec;
516     //iibin++;
517     //}
518 algomez 1.1 }
519    
520     hPVs["Nreweight"]->Fill( total_pvs, PUweight );
521    
522 algomez 1.2 /////////////
523     // HLT scale factor for MC
524     ////////////
525     double SF_hlt = 1.;
526     if (fIsMC) SF_hlt = 0.966;
527    
528     // LETS INCLUDE THE TRIGGER SF INTO THE PU WEIGHTS
529     PUweight = PUweight*SF_hlt;
530    
531 algomez 1.1 cutmap["Processed"] += PUweight;
532    
533     if (fVerbose) cout << "done pv" << endl;
534     //////////////////////////////////
535     // MUONS
536     //////////////////////////////////
537 algomez 1.2 int ngoodIDmuons = 0;
538 algomez 1.1 int nloosemuons = 0;
539     int ntightmuons = 0;
540     int nqcdmuons = 0;
541    
542     double RelIso = -1.;
543     double deltaR = -1.;
544     double QCDRelIso = -1.;
545     double QCDdeltaR = -1;
546    
547 algomez 1.2
548 algomez 1.1 for ( size_t imu=0; imu < total_muons; ++imu) {
549    
550     TopMuonEvent muon = muons[imu];
551 algomez 1.2 hmuons["N_muons"]->Fill( total_muons );
552 algomez 1.1
553     h1test->Fill( muon.pt );
554     hmuons["pt_cut1"]->Fill( muon.pt, PUweight );
555 algomez 1.2
556     if ( fMuSelector.MuonID( muon, PVz ) ) ngoodIDmuons++;
557 algomez 1.1
558     // select only good muons
559 algomez 1.2
560 algomez 1.1 if ( fMuSelector.MuonLoose( muon ) )
561     {
562    
563     nloosemuons++;
564    
565     if ( fMuSelector.MuonTight( muon, PVz) ) hmuons["pt_cut2"]->Fill( muon.pt, PUweight );
566     if ( fMuSelector.MuonTightDeltaR( muon, PVz, jets) ) {
567     ntightmuons++;
568     deltaR = fMuSelector.GetDeltaR();
569     }
570     p4muon.SetPtEtaPhiE( muon.pt, muon.eta, muon.phi, muon.e );
571 algomez 1.2 RelIso = muon.pfreliso; //muon.reliso03;
572 algomez 1.1
573     }
574     // check muon in QCD control region
575     if ( fMuSelector.MuonRelax02IsoQCD( muon, PVz, jets ) )
576     {
577     nqcdmuons++;
578     // keep the leading muon for selection
579     if (nqcdmuons==1) {
580     p4QCDmuon.SetPtEtaPhiE( muon.pt, muon.eta, muon.phi, muon.e );
581 algomez 1.2 QCDRelIso = muon.pfreliso; //muon.reliso03;
582 algomez 1.1 QCDdeltaR = fMuSelector.GetDeltaR();
583     }
584     }
585    
586     }
587    
588     if (fVerbose) cout << "done muons" << endl;
589     //////////////////////////////////
590     // ELECTRONS
591     //////////////////////////////////
592     int nlooseelectrons = 0;
593     int ntightelectrons = 0;
594     bool IsConversion = false;
595    
596     for ( size_t iele=0; iele < total_electrons; ++iele) {
597    
598     TopElectronEvent electron = electrons[iele];
599    
600     if ( fEleSelector.ElectronLoose(electron) ) nlooseelectrons++;
601    
602     if ( fEleSelector.ElectronTight(electron, PVz ) )
603     {
604    
605     if (ntightelectrons == 0)
606     {
607     IsConversion = electron.IsConversion;
608     p4ele.SetPtEtaPhiE( electron.pt, electron.eta, electron.phi, electron.e );
609     helectrons["pt_cut2"]->Fill( p4ele.Pt(), PUweight );
610     helectrons["eta_cut2"]->Fill( p4ele.Eta(), PUweight );
611     helectrons["phi_cut2"]->Fill( p4ele.Phi(), PUweight );
612    
613     }
614     ntightelectrons++;
615     }
616     }
617    
618     if (fVerbose) cout << "done electron" << endl;
619    
620     if ( fChannel == 1 ) // muon+jets
621     {
622    
623     if (fdoQCD2SideBand) {
624    
625     if (nqcdmuons == 0) return kTRUE;
626     cutmap["OneIsoMu"] += PUweight;
627    
628     p4lepton = p4QCDmuon;
629     RelIso = QCDRelIso;
630     deltaR = QCDdeltaR;
631     } else {
632    
633 algomez 1.2 if ( ngoodIDmuons > 0 ) hmuons["Ngood"]->Fill( total_pvs, PUweight);
634     if ( ntightmuons > 0 ) hmuons["Niso"]->Fill( total_pvs, PUweight);
635    
636 algomez 1.1 if ( ntightmuons != 1 ) return kTRUE;
637     cutmap["OneIsoMu"] += PUweight;
638 algomez 1.2
639 algomez 1.1 if ( nloosemuons > 1 ) return kTRUE;
640     cutmap["LooseMuVeto"] += PUweight;
641    
642     if ( nlooseelectrons > 0 ) return kTRUE;
643     cutmap["ElectronVeto"] += PUweight;
644    
645     p4lepton = p4muon;
646     if (fVerbose) cout << "got a good lepton" << endl;
647     }
648    
649     hmuons["pt"]->Fill( p4lepton.Pt(), PUweight );
650     hmuons["eta"]->Fill( p4lepton.Eta(), PUweight );
651     hmuons["phi"]->Fill( p4lepton.Phi(), PUweight );
652     hmuons["reliso"]->Fill( RelIso, PUweight );
653     hmuons["deltaR"]->Fill( deltaR, PUweight );
654    
655     }
656     else // electron+jets
657     {
658     // pending ...
659     }
660    
661     if (fVerbose) cout << "done lepton selection " << endl;
662    
663     /////////////////////////////////
664     // MET
665     /////////////////////////////////
666    
667     p4MET.SetPtEtaPhiE( ntuple->PFMET,
668     0.,
669     ntuple->PFMETphi,
670     ntuple->PFMET );
671    
672     //temporal check using genMET
673     //p4MET.SetPtEtaPhiE( ntuple->gen.MET,
674     // 0.,
675     // ntuple->gen.METPhi,
676     // ntuple->gen.MET );
677    
678     if (fdoQCD1SideBand && p4MET.Et() > 20.) return kTRUE;
679     else if ( p4MET.Et() <= 20. && fdoQCD2SideBand==false ) return kTRUE;
680    
681     if (fVerbose) cout << "pass MET cut" << endl;
682    
683     //cutmap["MET"] += PUweight;
684     hMET["MET"]->Fill( p4MET.Pt(), PUweight );
685     hMET["phi"]->Fill( p4MET.Phi(), PUweight );
686     hMET["Ht"]->Fill( ntuple->PFHt, PUweight );
687     hMET["Htlep"]->Fill( ntuple->PFHt + p4lepton.Pt(), PUweight );
688    
689     double Wpt = p4lepton.Pt() + p4MET.Pt();
690     double Wpx = p4lepton.Px() + p4MET.Px();
691     double Wpy = p4lepton.Py() + p4MET.Py();
692     double WMt = sqrt(Wpt*Wpt-Wpx*Wpx-Wpy*Wpy);
693    
694 algomez 1.2 if ( WMt <= 40. ) return kTRUE;
695 algomez 1.1 cutmap["MET"] += PUweight;
696    
697     if (fVerbose) cout << "pass W MT cut " << endl;
698    
699     /////////////////////////////////
700     // estimate Pz of neutrino
701     ////////////////////////////////
702    
703     fzCalculator.SetMET(p4MET);
704     fzCalculator.SetLepton(p4lepton);
705     if (fChannel==2)
706     {
707     fzCalculator.SetLeptonType("electron");
708     }
709    
710     double pzNu = fzCalculator.Calculate();
711     p4Nu = TLorentzVector();
712     p4OtherNu = TLorentzVector();
713    
714     p4Nu.SetPxPyPzE(p4MET.Px(), p4MET.Py(), pzNu, sqrt(p4MET.Px()*p4MET.Px()+p4MET.Py()*p4MET.Py()+pzNu*pzNu));
715     //print "pzNu = " +str(pzNu)
716     //print "p4Nu =("+str(p4Nu.Px())+","+str(p4Nu.Py())+","+str(p4Nu.Pz())+","+str(p4Nu.E())
717     double pzOtherNu = fzCalculator.getOther();
718     p4OtherNu.SetPxPyPzE( p4MET.Px(), p4MET.Py(),pzOtherNu,sqrt(p4MET.Px()*p4MET.Px()+p4MET.Py()*p4MET.Py()+pzOtherNu*pzOtherNu));
719    
720     double WmassNoPt = (p4Nu+p4lepton).M();
721    
722     if ( fzCalculator.IsComplex() )
723     {
724     double ptNu1 = fzCalculator.getPtneutrino(1);
725     double ptNu2 = fzCalculator.getPtneutrino(2);
726     TLorentzVector p4Nu1tmp;
727     TLorentzVector p4Nu2tmp;
728    
729     p4Nu1tmp.SetPxPyPzE( ptNu1*p4MET.Px()/p4MET.Pt(), ptNu1*p4MET.Py()/p4MET.Pt(), pzNu, sqrt(ptNu1*ptNu1+pzNu*pzNu));
730     p4Nu2tmp.SetPxPyPzE( ptNu2*p4MET.Px()/p4MET.Pt(), ptNu2*p4MET.Py()/p4MET.Pt(), pzNu, sqrt(ptNu2*ptNu2+pzNu*pzNu));
731    
732     TLorentzVector Wtmp;
733     Wtmp = p4lepton + p4Nu1tmp;
734     double Wm1 = 0;
735     double Wm2 = 0;
736     Wm1 = Wtmp.M();
737     Wtmp = p4lepton + p4Nu2tmp;
738     Wm2 = Wtmp.M();
739     if ( fabs( Wm1 - 80.) < fabs( Wm2 - 80.) ) p4Nu = p4Nu1tmp;
740     else p4Nu = p4Nu2tmp;
741    
742     p4OtherNu = p4Nu; // since we chose the real part, the two solutions are the same.
743     }
744    
745    
746     hMET["PzNu"]->Fill(pzNu, PUweight ); //change this to 2d with two sol and as a function of jets
747    
748     TLorentzVector p4LepW = p4lepton + p4Nu;
749     TLorentzVector p4OtherLepW = p4lepton + p4OtherNu;
750    
751     //hMET["LepWmass"]->Fill(p4LepW.M(), PUweight );
752     //if ( fzCalculator.IsComplex() ) hMET["LepWmassComplex"]->Fill( p4LepW.M(), PUweight );
753    
754    
755     /////////////////////////////////
756     // JETS
757     ////////////////////////////////
758    
759     //JetCombinatorics myCombi = JetCombinatorics();
760    
761     int njets = 0;
762     map< string, vector<float> > bdisc;
763     map< string, vector<bool> > isTagb;
764     vector<int> listflavor;
765    
766     for ( size_t ijet=0; ijet < total_jets; ++ijet)
767     {
768    
769     TopJetEvent jet = jets[ijet];
770     double SF_JEC = 1.;
771     if (fdoJECunc)
772     {
773     fJECunc->setJetEta( jet.eta);
774     fJECunc->setJetPt( jet.pt);
775     double jec_unc = fJECunc->getUncertainty(true);
776     if (fVerbose) cout << "JEC uncertainty is " << jec_unc << endl;
777     if (fdoJECup) SF_JEC = 1.+jec_unc;
778     else SF_JEC = 1.-jec_unc;
779     }
780    
781     if ( SF_JEC*jet.pt > 35. && fabs(jet.eta) < 2.4 && SF_JEC*jets[0].pt > 120. )
782     {
783     if (fVerbose) cout << " jet pt " << SF_JEC*jet.pt << endl;
784    
785     //hjets["pt"]->Fill( jet.pt, PUweight );
786     //hjets["eta"]->Fill( jet.eta, PUweight );
787     //hjets["phi"]->Fill( jet.phi, PUweight );
788    
789     TLorentzVector tmpjet;
790     tmpjet.SetPtEtaPhiE(SF_JEC*jet.pt, jet.eta, jet.phi, SF_JEC*jet.e);
791     p4jets.push_back( tmpjet);
792     listflavor.push_back( jet.mc.flavor );
793    
794     if (fVerbose) {
795     cout << "done storing njets " << njets << endl;
796     cout << " bdisc " << jet.btag_TCHP << endl;
797 algomez 1.2 cout << " bdisc " << jet.btag_CSV << endl;
798 algomez 1.1 }
799     // store discriminators
800     bdisc["TCHP"].push_back( jet.btag_TCHP );
801 algomez 1.2 bdisc["CSV"].push_back( jet.btag_CSV );
802 algomez 1.1 if (fVerbose) cout << "store bdisc" << endl;
803     // TCHPL cut at 1.19
804     // check TCHPM cut at 1.93
805     if ( jet.btag_TCHP > 1.93 ) isTagb["TCHPM"].push_back(true);
806     else isTagb["TCHPM"].push_back(false);
807     if (fVerbose) cout << "done tchpl" << endl;
808     // check SSVHEM cut at 1.74
809 algomez 1.2 if ( jet.btag_CSV > 0.679) isTagb["CSVM"].push_back(true);
810     else isTagb["CSVM"].push_back(false);
811     if (fVerbose) cout << "done csvm" << endl;
812 algomez 1.1 // CSVM cut at 0.679
813    
814     njets++;
815     }
816     }
817    
818     if (njets>0) hM["WMt"]->Fill( WMt, PUweight );
819    
820     hjets["Njets"]->Fill( njets, PUweight );
821    
822     if (fVerbose) cout << "done jets" << endl;
823    
824     if (njets > 0 ) cutmap["1Jet"] += PUweight;
825     if (njets > 1 ) cutmap["2Jet"] += PUweight;
826     if (njets > 2 ) cutmap["3Jet"] += PUweight;
827     if (njets > 3 ) cutmap["4Jet"] += PUweight;
828    
829     //if (njets == 1)
830     // {
831     // hjets["1st_pt"]->Fill( p4jets[0].Pt(), PUweight );
832     // }
833 algomez 1.2 /*
834     if (njets==2)
835     {
836     hjets["2nd_pt"]->Fill( p4jets[1].Pt(), PUweight );
837     }
838    
839     if (njets==3)
840     {
841     hjets["3rd_pt"]->Fill( p4jets[2].Pt(), PUweight );
842     }
843    
844     if (njets >= 4)
845     {
846     hjets["4th_pt"]->Fill( p4jets[3].Pt(), PUweight );
847     }
848     */
849 algomez 1.1
850     if (njets > 1 )
851     {
852     // count partons
853     int number_of_b = 0;
854     int number_of_c = 0;
855     int number_of_l = 0;
856     int number_of_b_highpt = 0;
857     int number_of_c_highpt = 0;
858     float SFb_0tag = 1.;
859     float SFb_only1tag = 1.;
860     float SFb_1tag = 1.;
861     float SFb_2tag = 1.;
862     //float SFb_0tag_syst[2] = {1.}; // for systematics
863     float SFb_1tag_syst[2] = {1.};
864     float SFb_2tag_syst[2] = {1.};
865     //float SFb_1tag_systhighpt[2] = {1.};
866    
867     for ( size_t kk=0; kk < p4jets.size(); ++kk)
868     {
869     hjets["pt"]->Fill( p4jets[kk].Pt(), PUweight );
870 algomez 1.2 if ( abs(listflavor[kk])==5 && p4jets[kk].Pt()<=240 ) { number_of_b++; hjets["pt_b_mc"]->Fill( p4jets[kk].Pt(), PUweight );}
871     if ( abs(listflavor[kk])==4 && p4jets[kk].Pt()<=240 ) { number_of_c++; hjets["pt_c_mc"]->Fill( p4jets[kk].Pt(), PUweight );}
872 algomez 1.1 if ( abs(listflavor[kk])==1 || abs(listflavor[kk])==2 || abs(listflavor[kk])==3 || abs(listflavor[kk])==21 )
873 algomez 1.2 { number_of_l++; hjets["pt_l_mc"]->Fill( p4jets[kk].Pt(), PUweight );}
874     if ( abs(listflavor[kk])==5 && p4jets[kk].Pt()>240 ) { number_of_b_highpt++; hjets["pt_b_mc"]->Fill( p4jets[kk].Pt(), PUweight );}
875     if ( abs(listflavor[kk])==4 && p4jets[kk].Pt()>240 ) { number_of_c_highpt++; hjets["pt_c_mc"]->Fill( p4jets[kk].Pt(), PUweight );}
876 algomez 1.1
877     if ( isTagb["TCHPM"][kk] )
878     {
879     hjets["pt_btag"]->Fill( p4jets[kk].Pt(), PUweight );
880     if ( abs(listflavor[kk])==5 ) hjets["pt_btag_b"]->Fill( p4jets[kk].Pt(), PUweight );
881     if ( abs(listflavor[kk])==4 ) hjets["pt_btag_c"]->Fill( p4jets[kk].Pt(), PUweight );
882     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 );
883     }
884     }
885    
886     hPVs["Nreweight_2jet"]->Fill( total_pvs, PUweight );
887     hjets["1st_pt"]->Fill( p4jets[0].Pt(), PUweight );
888     hjets["2nd_pt"]->Fill( p4jets[1].Pt(), PUweight );
889     hjets["1st_eta"]->Fill( p4jets[0].Eta(), PUweight );
890     hjets["2nd_eta"]->Fill( p4jets[1].Eta(), PUweight );
891    
892     hmuons["pt_2jet"]->Fill( p4lepton.Pt(), PUweight );
893     hmuons["reliso_2jet"]->Fill( RelIso, PUweight );
894     hMET["MET_2jet"]->Fill( p4MET.Pt(), PUweight );
895 algomez 1.2
896 algomez 1.1 if (fIsMC) {
897     hMET["genMET_2jet"]->Fill( ntuple->gen.MET, PUweight );
898     hMET["deltaMET_2jet"]->Fill( p4MET.Pt() - ntuple->gen.MET, PUweight );
899     }
900    
901     hM["WMt_2jet"]->Fill( WMt, PUweight );
902    
903     TLorentzVector p4Dijet = p4jets[0] + p4jets[1];
904     double Dijet_deltaPhi = p4jets[0].DeltaPhi(p4jets[1]);
905     TLorentzVector p4Top = p4jets[2] + p4LepW;
906    
907     TLorentzVector p4Wprime = p4LepW + p4Dijet;
908     hM["Wprime"]->Fill( p4Wprime.M(), PUweight );
909    
910     if (fVerbose) cout << "Wprime mass calculated" << endl;
911    
912     // count number of b-tags
913     int Nbtags_TCHPM = 0;
914 algomez 1.2 int Nbtags_CSVM = 0;
915 algomez 1.1 for ( size_t itag=0; itag< isTagb["TCHPM"].size(); ++itag )
916     {
917     if ( isTagb["TCHPM"][itag] ) Nbtags_TCHPM++;
918 algomez 1.2 if ( isTagb["CSVM"][itag] ) Nbtags_CSVM++;
919 algomez 1.1 }
920    
921     // compute b-tag event weight
922     if ( fIsMC )
923     {
924    
925     // zeto tag
926     BTagWeight b0(0,0); // number of tags
927     BTagWeight::JetInfo bj(0.63,0.91); // mean MC eff and mean SF. For TCHPM=0.91\pm0.09, CSVM=0.96\pm0.096
928     BTagWeight::JetInfo cj(0.15,0.91);
929     double light_mceff = 0.017;
930     if ( 100 < p4jets[0].Pt() && p4jets[0].Pt() <= 200 ) light_mceff = 0.04;
931     if ( 200 < p4jets[0].Pt() && p4jets[0].Pt() <= 300 ) light_mceff = 0.08;
932     if ( 300 < p4jets[0].Pt() && p4jets[0].Pt() <= 400 ) light_mceff = 0.12;
933     if ( 400 < p4jets[0].Pt() ) light_mceff = 0.14;
934     BTagWeight::JetInfo lj(light_mceff,1.22); //for TCHPM=1.22, CSVM=1.08 \pm 0.13
935     // b-tag systematic UP 9% for b, 18% for c
936     BTagWeight::JetInfo bjUP(0.63,0.99);
937     BTagWeight::JetInfo cjUP(0.15,1.07);
938     // b-tag systemacit DOWN 9% for b, 18% for c
939     BTagWeight::JetInfo bjDOWN(0.63,0.83);
940     BTagWeight::JetInfo cjDOWN(0.15,0.75);
941     // for high pt jets > 240 UP 50% for b and c
942     BTagWeight::JetInfo bjUPhighpt(0.63,1.36);
943     BTagWeight::JetInfo cjUPhighpt(0.15,1.36);
944     // for high pt jets > 240 UP 50% for b and c
945     BTagWeight::JetInfo bjDOWNhighpt(0.63,0.46);
946     BTagWeight::JetInfo cjDOWNhighpt(0.15,0.46);
947    
948     vector<BTagWeight::JetInfo> j;
949     for(int i=0;i<number_of_b;i++)j.push_back(bj);
950     for(int i=0;i<number_of_b_highpt;i++)j.push_back(bj);
951     for(int i=0;i<number_of_c;i++)j.push_back(cj);
952     for(int i=0;i<number_of_c_highpt;i++)j.push_back(cj);
953     for(int i=0;i<number_of_l;i++)j.push_back(lj);
954    
955     if (Nbtags_TCHPM==0) {
956     SFb_0tag = b0.weight(j,0);
957     hjets["Nbtags_TCHPM"]->Fill( Nbtags_TCHPM, PUweight*SFb_0tag ); // fill bin 0
958     }
959     // only one tag
960     BTagWeight b11(1,1); // number of tags
961     if (Nbtags_TCHPM==1) {
962     SFb_only1tag = b11.weight(j,1);
963     hjets["Nbtags_TCHPM"]->Fill( Nbtags_TCHPM, PUweight*SFb_only1tag ); // fill bin 1
964     }
965     // at least one tag
966     BTagWeight b1(1,Nbtags_TCHPM); // number of tags
967     if (Nbtags_TCHPM>=1) {
968     SFb_1tag = b1.weight(j,1);
969    
970     vector<BTagWeight::JetInfo> jj;
971     for(int i=0;i<number_of_b;i++)jj.push_back(bjUP);
972     for(int i=0;i<number_of_b_highpt;i++)jj.push_back(bjUPhighpt);
973     for(int i=0;i<number_of_c;i++)jj.push_back(cjUP);
974     for(int i=0;i<number_of_c_highpt;i++)jj.push_back(cjUPhighpt);
975     for(int i=0;i<number_of_l;i++)jj.push_back(lj);
976     SFb_1tag_syst[0] = b1.weight(jj,1); //UP
977     vector<BTagWeight::JetInfo> jk;
978     for(int i=0;i<number_of_b;i++)jk.push_back(bjDOWN);
979     for(int i=0;i<number_of_b_highpt;i++)jk.push_back(bjDOWNhighpt);
980     for(int i=0;i<number_of_c;i++)jk.push_back(cjDOWN);
981     for(int i=0;i<number_of_c_highpt;i++)jk.push_back(cjDOWNhighpt);
982     for(int i=0;i<number_of_l;i++)jk.push_back(lj);
983     SFb_1tag_syst[1] = b1.weight(jk,1);//DOWN
984    
985     }
986     // at least two tags
987     BTagWeight b2(2,Nbtags_TCHPM); // number of tags
988     if (Nbtags_TCHPM>=2) {
989     SFb_2tag = b2.weight(j,2);
990     hjets["Nbtags_TCHPM"]->Fill( Nbtags_TCHPM, PUweight*SFb_2tag ); // fill bin >=2
991    
992     vector<BTagWeight::JetInfo> jj;
993     for(int i=0;i<number_of_b;i++)jj.push_back(bjUP);
994     for(int i=0;i<number_of_b_highpt;i++)jj.push_back(bjUPhighpt);
995     for(int i=0;i<number_of_c;i++)jj.push_back(cjUP);
996     for(int i=0;i<number_of_c_highpt;i++)jj.push_back(cjUPhighpt);
997     for(int i=0;i<number_of_l;i++)jj.push_back(lj);
998     SFb_2tag_syst[0] = b2.weight(jj,2); //UP
999    
1000     vector<BTagWeight::JetInfo> jk;
1001     for(int i=0;i<number_of_b;i++)jk.push_back(bjDOWN);
1002     for(int i=0;i<number_of_b_highpt;i++)jk.push_back(bjDOWNhighpt);
1003     for(int i=0;i<number_of_c;i++)jk.push_back(cjDOWN);
1004     for(int i=0;i<number_of_c_highpt;i++)jk.push_back(cjDOWNhighpt);
1005     for(int i=0;i<number_of_l;i++)jk.push_back(lj);
1006     SFb_2tag_syst[1] = b2.weight(jk,2);//DOWN
1007     }
1008     }
1009     else hjets["Nbtags_TCHPM"]->Fill( Nbtags_TCHPM );
1010    
1011     //hjets["Nbtags_TCHPM"]->Fill( Nbtags_TCHPM, PUweight*SFb );
1012 algomez 1.2 //hjets["Nbtags_CSVM"]->Fill( Nbtags_CSVM, PUweight*SFb );
1013 algomez 1.1
1014 algomez 1.2 if ( Nbtags_TCHPM == 0 && njets < 5)
1015 algomez 1.1 {
1016 algomez 1.2 // compute top mass
1017     int top_bjet_index = 0;
1018     double deltaM = 99999.;
1019     //double top_mass = 0.;
1020     //TLorentzVector p4Top;
1021     TLorentzVector bestp4Top;
1022     for ( size_t itejet = 0; itejet < p4jets.size(); ++itejet )
1023     {
1024    
1025     for ( int isolw =0; isolw<2; ++ isolw)
1026     {
1027    
1028     p4Top = p4jets[itejet] + p4LepW;
1029     if (isolw==1) p4Top = p4jets[itejet] + p4OtherLepW;
1030    
1031     if ( fabs( p4Top.M() - 172.) < deltaM ) {
1032     top_bjet_index = itejet;
1033     //top_mass = p4Top.M();
1034     deltaM = fabs( p4Top.M() - 172.);
1035     bestp4Top = p4Top;
1036     }
1037    
1038     }
1039     }
1040     p4Top = bestp4Top;
1041     bool passcutWlep = false;
1042     if ( p4LepW.M() < 90 ) passcutWlep = true;
1043    
1044     if (passcutWlep) {
1045    
1046     if ( top_bjet_index == 0 ) {
1047     p4Wprime = p4Top + p4jets[1];
1048 algomez 1.1 }
1049 algomez 1.2 else {
1050     p4Wprime = p4Top + p4jets[0];
1051 algomez 1.1 }
1052 algomez 1.2
1053     bool passcut = true;
1054     if (fdoMtopCut && ( p4Top.M() <= 130 || p4Top.M() >= 210 ) ) passcut = false;
1055    
1056     if (passcut) {
1057 algomez 1.1
1058 algomez 1.2 hM["Wprime_0btag"]->Fill( p4Wprime.M(), PUweight*SFb_0tag );
1059     if (fSample=="WJets")
1060     {
1061     int FH = ntuple->flavorHistory;
1062     if ( FH == 1 || FH == 2 || FH == 5 || FH == 7 || FH == 9 ) hM["Wprime_0btag_bb"]->Fill( p4Wprime.M(), PUweight*SFb_0tag );
1063     else if ( FH == 3 || FH == 4 || FH == 6 || FH == 8 || FH == 10) hM["Wprime_0btag_cc"]->Fill( p4Wprime.M(), PUweight*SFb_0tag );
1064     else if ( FH == 11 ) hM["Wprime_0btag_light"]->Fill( p4Wprime.M(), PUweight*SFb_0tag );
1065     }
1066     if (fSample.Contains("Wprime"))
1067     {
1068     if ( ntuple->gen.LeptonicChannel == 11 ) hM["Wprime_0btag_MCsemie"]->Fill( p4Wprime.M(), PUweight*SFb_0tag );
1069     else if ( ntuple->gen.LeptonicChannel == 13 ) hM["Wprime_0btag_MCsemimu"]->Fill( p4Wprime.M(), PUweight*SFb_0tag );
1070     else hM["Wprime_0btag_MChad"]->Fill( p4Wprime.M(), PUweight*SFb_0tag );
1071     }
1072     } // top mass cut
1073     }// Wlep cut
1074     }// 0 tag
1075 algomez 1.1
1076     // check if the two leading jets have a b jet
1077     if ( Nbtags_TCHPM >= 1 && (isTagb["TCHPM"][0] || isTagb["TCHPM"][1]) )
1078     {
1079    
1080     //cutmap["2Jet1b"] += PUweight*SFb_1tag;
1081     hjets["Njets_1tag"]->Fill( njets, PUweight*SFb_1tag );
1082     if (Nbtags_TCHPM > 1 ) hjets["Njets_2tag"]->Fill( njets, PUweight*SFb_2tag );
1083    
1084     if ( njets < 5 ) {
1085    
1086     //cutmap["2Jet1b"] += PUweight*SFb_1tag;
1087     // calculate dijet mass closest to W mass
1088     double the_dijet_mass = 0;
1089     double sigma2 = 13.*13.;
1090     double tmpchi2 = 999999.;
1091     for ( size_t itejet = 0; itejet < p4jets.size(); ++itejet )
1092     {
1093     for ( size_t jtejet = itejet+1; jtejet < p4jets.size(); ++jtejet )
1094     {
1095     TLorentzVector tmpvv = p4jets[itejet] + p4jets[jtejet];
1096     double tmpmass = tmpvv.M();
1097     if ( (tmpmass - 83.8 )*(tmpmass - 83.8)/sigma2 < tmpchi2 ) { tmpchi2 = (tmpmass - 83.8 )*(tmpmass - 83.8)/sigma2; the_dijet_mass = tmpmass; }
1098     }
1099     }
1100     hM["dijet"]->Fill( the_dijet_mass, PUweight*SFb_1tag );
1101     hjets["Dijet_deltaPhi"]->Fill( Dijet_deltaPhi, PUweight*SFb_1tag );
1102    
1103     // compute jjj mass for events with >=3 jets
1104     if ( njets > 2 ) {
1105     TLorentzVector p4Trijet = p4jets[0] + p4jets[1] + p4jets[2];
1106     hM["trijet"]->Fill( p4Trijet.M(), PUweight*SFb_1tag );
1107     }
1108     // compute top mass
1109     int top_bjet_index = 0;
1110     double deltaM = 99999.;
1111     //double top_mass = 0.;
1112     //TLorentzVector p4Top;
1113     TLorentzVector bestp4Top;
1114     for ( size_t itejet = 0; itejet < p4jets.size(); ++itejet )
1115     {
1116    
1117     for ( int isolw =0; isolw<2; ++ isolw)
1118     {
1119    
1120     p4Top = p4jets[itejet] + p4LepW;
1121     if (isolw==1) p4Top = p4jets[itejet] + p4OtherLepW;
1122    
1123     if ( fabs( p4Top.M() - 172.) < deltaM ) {
1124     top_bjet_index = itejet;
1125     //top_mass = p4Top.M();
1126     deltaM = fabs( p4Top.M() - 172.);
1127     bestp4Top = p4Top;
1128     }
1129    
1130     }
1131     }
1132     p4Top = bestp4Top;
1133     hM["top_1btag"]->Fill( p4Top.M(), PUweight*SFb_1tag );
1134     hMET["LepWmass"]->Fill(p4LepW.M(), PUweight*SFb_1tag );
1135     hMET["LepWmassNoPt"]->Fill(WmassNoPt, PUweight*SFb_1tag );
1136    
1137     bool passcutWlep = false;
1138 algomez 1.2 if ( p4LepW.M() < 90 ) passcutWlep = true;
1139 algomez 1.1
1140     if (passcutWlep) {
1141    
1142     if ( p4Top.M() > 130 && p4Top.M() < 210 ) hMET["LepWmass_topcut"]->Fill(p4LepW.M(), PUweight*SFb_1tag );
1143    
1144     double tb_deltaPhi = 0.;
1145     double tb_deltaR = 0.;
1146 algomez 1.2 double tb_deltaEta = -999;
1147     double pt_wprime = 0;
1148     double pt_top = 0;
1149     double pt_b = 0;
1150 algomez 1.1 if ( top_bjet_index == 0 ) {
1151     p4Wprime = p4Top + p4jets[1];
1152     tb_deltaPhi = p4Top.DeltaPhi( p4jets[1] );
1153     tb_deltaR = p4Top.DeltaR( p4jets[1] );
1154 algomez 1.2 tb_deltaEta = p4Top.Eta() - p4jets[1].Eta();
1155     pt_wprime = p4Wprime.Pt();
1156     pt_top = p4Top.Pt();
1157     pt_b = p4jets[1].Pt();
1158 algomez 1.1 }
1159     else {
1160     p4Wprime = p4Top + p4jets[0];
1161     tb_deltaPhi = fabs( p4Top.DeltaPhi( p4jets[0] ) );
1162     tb_deltaR = p4Top.DeltaR( p4jets[0] );
1163 algomez 1.2 tb_deltaEta = p4Top.Eta() - p4jets[0].Eta();
1164     pt_wprime= p4Wprime.Pt();
1165     pt_top = p4Top.Pt();
1166     pt_b = p4jets[0].Pt();
1167 algomez 1.1 }
1168    
1169     bool passcut = true;
1170     if (fdoMtopCut && ( p4Top.M() <= 130 || p4Top.M() >= 210 ) ) passcut = false;
1171    
1172     if (passcut) {
1173     cutmap["2Jet1b"] += PUweight*SFb_1tag;
1174     h2_pt_Wprime->Fill( p4Top.Pt(), p4Wprime.M(), PUweight*SFb_1tag );
1175     hM["Wprime_1btag"]->Fill( p4Wprime.M(), PUweight*SFb_1tag );
1176 algomez 1.2 // FINAL Plots
1177 algomez 1.1 hMET["deltaPhi"]->Fill( p4lepton.DeltaPhi( p4MET ), PUweight*SFb_1tag );
1178 algomez 1.2 hjets["tb_deltaPhi"]->Fill( tb_deltaPhi, PUweight*SFb_1tag );
1179     hjets["tb_deltaEta"]->Fill( tb_deltaEta, PUweight*SFb_1tag );
1180     //hjets["tb_deltaR"]->Fill( tb_deltaR, PUweight*SFb_1tag );
1181     hjets["pt_Wprime"]->Fill( pt_wprime, PUweight*SFb_1tag );
1182     hjets["pt_top"]->Fill( pt_top, PUweight*SFb_1tag );
1183     hjets["pt_b"]->Fill( pt_b, PUweight*SFb_1tag );
1184     hjets["1st_pt_fin"]->Fill( p4jets[0].Pt(), PUweight*SFb_1tag );
1185     hjets["2nd_pt_fin"]->Fill( p4jets[1].Pt(), PUweight*SFb_1tag );
1186     hjets["1st_eta_fin"]->Fill( p4jets[0].Eta(), PUweight*SFb_1tag );
1187     hjets["2nd_eta_fin"]->Fill( p4jets[1].Eta(), PUweight*SFb_1tag );
1188     hmuons["pt_fin"]->Fill( p4lepton.Pt(), PUweight*SFb_1tag );
1189     hmuons["reliso_fin"]->Fill( RelIso, PUweight*SFb_1tag );
1190     hMET["MET_fin"]->Fill( p4MET.Pt(), PUweight*SFb_1tag );
1191     hM["WMt_fin"]->Fill( WMt, PUweight*SFb_1tag );
1192 algomez 1.1
1193 algomez 1.2 if ( p4Wprime.M() > 1500.)
1194 algomez 1.1 {
1195     TString outstring = "run: ";
1196     outstring += TString(Form("%i",ntuple->run)) +" lumi: "+ TString(Form("%i",ntuple->lumi)) + " event: " + TString(Form("%i",ntuple->event));
1197     Info("Process",outstring);
1198     }
1199    
1200     if (fIsMC)
1201     {
1202     hM["Wprime_1btag_systUp"]->Fill( p4Wprime.M(), PUweight*SFb_1tag_syst[0] );
1203     hM["Wprime_1btag_systDown"]->Fill( p4Wprime.M(), PUweight*SFb_1tag_syst[1] );
1204     }
1205    
1206     if (fSample.Contains("Wprime"))
1207     {
1208     if ( ntuple->gen.LeptonicChannel == 11 ) hM["Wprime_1btag_MCsemie"]->Fill( p4Wprime.M(), PUweight*SFb_1tag );
1209     else if ( ntuple->gen.LeptonicChannel == 13 ) hM["Wprime_1btag_MCsemimu"]->Fill( p4Wprime.M(), PUweight*SFb_1tag );
1210     else hM["Wprime_1btag_MChad"]->Fill( p4Wprime.M(), PUweight*SFb_1tag );
1211     }
1212    
1213     if (fSample=="WJets")
1214     {
1215     int FH = ntuple->flavorHistory;
1216     if ( FH == 1 || FH == 2 || FH == 5 || FH == 7 || FH == 9 ) hM["Wprime_1btag_bb"]->Fill( p4Wprime.M(), PUweight*SFb_1tag );
1217     else if ( FH == 3 || FH == 4 || FH == 6 || FH == 8 || FH == 10) hM["Wprime_1btag_cc"]->Fill( p4Wprime.M(), PUweight*SFb_1tag );
1218     else if ( FH == 11 ) hM["Wprime_1btag_light"]->Fill( p4Wprime.M(), PUweight*SFb_1tag );
1219     }
1220    
1221     // check two b-tags
1222     if ( Nbtags_TCHPM > 1 )
1223     {
1224     cutmap["2Jet2b"] += PUweight*SFb_2tag;
1225     hM["Wprime_2btag"]->Fill( p4Wprime.M(), PUweight*SFb_2tag );
1226    
1227     if (fIsMC)
1228     {
1229     hM["Wprime_2btag_systUp"]->Fill( p4Wprime.M(), PUweight*SFb_2tag_syst[0] );
1230     hM["Wprime_2btag_systDown"]->Fill( p4Wprime.M(), PUweight*SFb_2tag_syst[1] );
1231     }
1232     }
1233    
1234     }//passcut mtop cut if requested
1235     }//passcutWlep cut
1236     }// njets < 5
1237    
1238     }
1239    
1240     if ( Nbtags_TCHPM == 1 ) hM["Wprime_1onlybtag"]->Fill( p4Wprime.M(), PUweight*SFb_only1tag );
1241    
1242     }
1243    
1244    
1245     if (fVerbose) cout << "done analysis" << endl;
1246    
1247     return kTRUE;
1248     }
1249    
1250     void Analyzer::SlaveTerminate()
1251     {
1252     // The SlaveTerminate() function is called after all entries or objects
1253     // have been processed. When running with PROOF SlaveTerminate() is called
1254     // on each slave server.
1255    
1256     // fill cutflow histogram
1257    
1258     int ibin = 1;
1259     for ( vector<string>::const_iterator ivec= fCutLabels.begin(); ivec != fCutLabels.end(); ++ivec )
1260     // for ( map<string, int >::const_iterator imap=cutmap.begin(); imap!=cutmap.end(); ++imap )
1261     {
1262     hcutflow->SetBinContent( ibin, cutmap[ *ivec ] );
1263     ibin++;
1264     }
1265     // Write the ntuple to the file
1266     if (fFile) {
1267     Bool_t cleanup = kFALSE;
1268     TDirectory *savedir = gDirectory;
1269     if ( h1test->GetEntries() > 0) {
1270     fFile->cd();
1271     h1test->Write();
1272     hcutflow->Write();
1273     h2_pt_Wprime->Write();
1274     fFile->mkdir("muons");
1275     fFile->cd("muons");
1276     for ( map<string,TH1* >::const_iterator imap=hmuons.begin(); imap!=hmuons.end(); ++imap )
1277     {
1278     TH1 *temp = imap->second;
1279     if ( temp->GetEntries() > 0 )
1280     temp->Write();
1281     //else cout << "Warning: empty histogram " << temp->GetName() << " will not be written to file." << endl;
1282     }
1283     fFile->cd();
1284     fFile->mkdir("PVs");
1285     fFile->cd("PVs");
1286     for ( map<string,TH1* >::const_iterator imap=hPVs.begin(); imap!=hPVs.end(); ++imap )
1287     {
1288     TH1 *temp = imap->second;
1289     if ( temp->GetEntries() > 0 )
1290     temp->Write();
1291     //else cout << "Warning: empty histogram " << temp->GetName() << " will not be written to file." << endl;
1292     }
1293     fFile->cd();
1294     fFile->mkdir("jets");
1295     fFile->cd("jets");
1296     for ( map<string,TH1* >::const_iterator imap=hjets.begin(); imap!=hjets.end(); ++imap )
1297     {
1298     TH1 *temp = imap->second;
1299     if ( temp->GetEntries() > 0 )
1300     temp->Write();
1301     }
1302     fFile->cd();
1303     fFile->mkdir("mass");
1304     fFile->cd("mass");
1305     for ( map<string,TH1* >::const_iterator imap=hM.begin(); imap!=hM.end(); ++imap )
1306     {
1307     TH1 *temp = imap->second;
1308     if ( temp->GetEntries() > 0 )
1309     temp->Write();
1310     }
1311     fFile->cd();
1312     fFile->mkdir("MET");
1313     fFile->cd("MET");
1314     for ( map<string,TH1* >::const_iterator imap=hMET.begin(); imap!=hMET.end(); ++imap )
1315     {
1316     TH1 *temp = imap->second;
1317     if ( temp->GetEntries() > 0 )
1318     temp->Write();
1319     }
1320    
1321     fFile->cd();
1322    
1323     fProofFile->Print();
1324     fOutput->Add(fProofFile);
1325     } else {
1326     cleanup = kTRUE;
1327     }
1328    
1329     h2_pt_Wprime->SetDirectory(0);
1330    
1331     h1test->SetDirectory(0);
1332     hcutflow->SetDirectory(0);
1333     gDirectory = savedir;
1334     fFile->Close();
1335     // Cleanup, if needed
1336     if (cleanup) {
1337     TUrl uf(*(fFile->GetEndpointUrl()));
1338     SafeDelete(fFile);
1339     gSystem->Unlink(uf.GetFile());
1340     SafeDelete(fProofFile);
1341     }
1342     }
1343    
1344     }
1345    
1346     void Analyzer::Terminate()
1347     {
1348     // The Terminate() function is the last function to be called during
1349     // a query. It always runs on the client, it can be used to present
1350     // the results graphically or save the results to file.
1351    
1352     Info("Terminate","Analyzer done.");
1353     }