ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/algomez/FourTop/TreeAnalyzer/test/Analyzer.C
Revision: 1.5
Committed: Mon Dec 12 18:49:31 2011 UTC (13 years, 4 months ago) by algomez
Content type: text/plain
Branch: MAIN
Changes since 1.4: +5 -0 lines
Log Message:
fix for 4tops samples

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