ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/algomez/FourTop/TreeAnalyzer/test/Analyzer.C
Revision: 1.3
Committed: Fri Dec 9 18:38:04 2011 UTC (13 years, 5 months ago) by algomez
Content type: text/plain
Branch: MAIN
Changes since 1.2: +17 -16 lines
Log Message:
getting rid of all the Wprime analysis

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 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)
240 algomez 1.1 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.3 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 algomez 1.3
690 algomez 1.1 double Wpt = p4lepton.Pt() + p4MET.Pt();
691     double Wpx = p4lepton.Px() + p4MET.Px();
692     double Wpy = p4lepton.Py() + p4MET.Py();
693     double WMt = sqrt(Wpt*Wpt-Wpx*Wpx-Wpy*Wpy);
694 algomez 1.3
695     //if ( WMt <= 40. ) return kTRUE;
696 algomez 1.1 cutmap["MET"] += PUweight;
697    
698     if (fVerbose) cout << "pass W MT cut " << endl;
699    
700     /////////////////////////////////
701     // estimate Pz of neutrino
702     ////////////////////////////////
703    
704     fzCalculator.SetMET(p4MET);
705     fzCalculator.SetLepton(p4lepton);
706     if (fChannel==2)
707     {
708     fzCalculator.SetLeptonType("electron");
709     }
710    
711     double pzNu = fzCalculator.Calculate();
712     p4Nu = TLorentzVector();
713     p4OtherNu = TLorentzVector();
714    
715     p4Nu.SetPxPyPzE(p4MET.Px(), p4MET.Py(), pzNu, sqrt(p4MET.Px()*p4MET.Px()+p4MET.Py()*p4MET.Py()+pzNu*pzNu));
716     //print "pzNu = " +str(pzNu)
717     //print "p4Nu =("+str(p4Nu.Px())+","+str(p4Nu.Py())+","+str(p4Nu.Pz())+","+str(p4Nu.E())
718     double pzOtherNu = fzCalculator.getOther();
719     p4OtherNu.SetPxPyPzE( p4MET.Px(), p4MET.Py(),pzOtherNu,sqrt(p4MET.Px()*p4MET.Px()+p4MET.Py()*p4MET.Py()+pzOtherNu*pzOtherNu));
720    
721     double WmassNoPt = (p4Nu+p4lepton).M();
722    
723     if ( fzCalculator.IsComplex() )
724     {
725     double ptNu1 = fzCalculator.getPtneutrino(1);
726     double ptNu2 = fzCalculator.getPtneutrino(2);
727     TLorentzVector p4Nu1tmp;
728     TLorentzVector p4Nu2tmp;
729    
730     p4Nu1tmp.SetPxPyPzE( ptNu1*p4MET.Px()/p4MET.Pt(), ptNu1*p4MET.Py()/p4MET.Pt(), pzNu, sqrt(ptNu1*ptNu1+pzNu*pzNu));
731     p4Nu2tmp.SetPxPyPzE( ptNu2*p4MET.Px()/p4MET.Pt(), ptNu2*p4MET.Py()/p4MET.Pt(), pzNu, sqrt(ptNu2*ptNu2+pzNu*pzNu));
732    
733     TLorentzVector Wtmp;
734     Wtmp = p4lepton + p4Nu1tmp;
735     double Wm1 = 0;
736     double Wm2 = 0;
737     Wm1 = Wtmp.M();
738     Wtmp = p4lepton + p4Nu2tmp;
739     Wm2 = Wtmp.M();
740     if ( fabs( Wm1 - 80.) < fabs( Wm2 - 80.) ) p4Nu = p4Nu1tmp;
741     else p4Nu = p4Nu2tmp;
742    
743     p4OtherNu = p4Nu; // since we chose the real part, the two solutions are the same.
744     }
745    
746    
747     hMET["PzNu"]->Fill(pzNu, PUweight ); //change this to 2d with two sol and as a function of jets
748    
749     TLorentzVector p4LepW = p4lepton + p4Nu;
750     TLorentzVector p4OtherLepW = p4lepton + p4OtherNu;
751    
752     //hMET["LepWmass"]->Fill(p4LepW.M(), PUweight );
753     //if ( fzCalculator.IsComplex() ) hMET["LepWmassComplex"]->Fill( p4LepW.M(), PUweight );
754    
755    
756     /////////////////////////////////
757     // JETS
758     ////////////////////////////////
759    
760     //JetCombinatorics myCombi = JetCombinatorics();
761    
762     int njets = 0;
763     map< string, vector<float> > bdisc;
764     map< string, vector<bool> > isTagb;
765     vector<int> listflavor;
766    
767     for ( size_t ijet=0; ijet < total_jets; ++ijet)
768     {
769    
770     TopJetEvent jet = jets[ijet];
771     double SF_JEC = 1.;
772     if (fdoJECunc)
773     {
774     fJECunc->setJetEta( jet.eta);
775     fJECunc->setJetPt( jet.pt);
776     double jec_unc = fJECunc->getUncertainty(true);
777     if (fVerbose) cout << "JEC uncertainty is " << jec_unc << endl;
778     if (fdoJECup) SF_JEC = 1.+jec_unc;
779     else SF_JEC = 1.-jec_unc;
780     }
781    
782     if ( SF_JEC*jet.pt > 35. && fabs(jet.eta) < 2.4 && SF_JEC*jets[0].pt > 120. )
783     {
784     if (fVerbose) cout << " jet pt " << SF_JEC*jet.pt << endl;
785    
786     //hjets["pt"]->Fill( jet.pt, PUweight );
787     //hjets["eta"]->Fill( jet.eta, PUweight );
788     //hjets["phi"]->Fill( jet.phi, PUweight );
789    
790     TLorentzVector tmpjet;
791     tmpjet.SetPtEtaPhiE(SF_JEC*jet.pt, jet.eta, jet.phi, SF_JEC*jet.e);
792     p4jets.push_back( tmpjet);
793     listflavor.push_back( jet.mc.flavor );
794    
795     if (fVerbose) {
796     cout << "done storing njets " << njets << endl;
797     cout << " bdisc " << jet.btag_TCHP << endl;
798 algomez 1.2 cout << " bdisc " << jet.btag_CSV << endl;
799 algomez 1.1 }
800     // store discriminators
801     bdisc["TCHP"].push_back( jet.btag_TCHP );
802 algomez 1.2 bdisc["CSV"].push_back( jet.btag_CSV );
803 algomez 1.1 if (fVerbose) cout << "store bdisc" << endl;
804     // TCHPL cut at 1.19
805     // check TCHPM cut at 1.93
806     if ( jet.btag_TCHP > 1.93 ) isTagb["TCHPM"].push_back(true);
807     else isTagb["TCHPM"].push_back(false);
808     if (fVerbose) cout << "done tchpl" << endl;
809     // check SSVHEM cut at 1.74
810 algomez 1.2 if ( jet.btag_CSV > 0.679) isTagb["CSVM"].push_back(true);
811     else isTagb["CSVM"].push_back(false);
812     if (fVerbose) cout << "done csvm" << endl;
813 algomez 1.1 // CSVM cut at 0.679
814    
815     njets++;
816     }
817     }
818    
819 algomez 1.3 if (njets>0) hM["WMt"]->Fill( WMt, PUweight );
820 algomez 1.1
821     hjets["Njets"]->Fill( njets, PUweight );
822    
823     if (fVerbose) cout << "done jets" << endl;
824    
825     if (njets > 0 ) cutmap["1Jet"] += PUweight;
826     if (njets > 1 ) cutmap["2Jet"] += PUweight;
827     if (njets > 2 ) cutmap["3Jet"] += PUweight;
828     if (njets > 3 ) cutmap["4Jet"] += PUweight;
829    
830     //if (njets == 1)
831     // {
832     // hjets["1st_pt"]->Fill( p4jets[0].Pt(), PUweight );
833     // }
834 algomez 1.2 /*
835     if (njets==2)
836     {
837     hjets["2nd_pt"]->Fill( p4jets[1].Pt(), PUweight );
838     }
839    
840     if (njets==3)
841     {
842     hjets["3rd_pt"]->Fill( p4jets[2].Pt(), PUweight );
843     }
844    
845     if (njets >= 4)
846     {
847     hjets["4th_pt"]->Fill( p4jets[3].Pt(), PUweight );
848     }
849     */
850 algomez 1.1
851     if (njets > 1 )
852     {
853     // count partons
854     int number_of_b = 0;
855     int number_of_c = 0;
856     int number_of_l = 0;
857     int number_of_b_highpt = 0;
858     int number_of_c_highpt = 0;
859     float SFb_0tag = 1.;
860     float SFb_only1tag = 1.;
861     float SFb_1tag = 1.;
862     float SFb_2tag = 1.;
863     //float SFb_0tag_syst[2] = {1.}; // for systematics
864     float SFb_1tag_syst[2] = {1.};
865     float SFb_2tag_syst[2] = {1.};
866     //float SFb_1tag_systhighpt[2] = {1.};
867    
868     for ( size_t kk=0; kk < p4jets.size(); ++kk)
869     {
870     hjets["pt"]->Fill( p4jets[kk].Pt(), PUweight );
871 algomez 1.2 if ( abs(listflavor[kk])==5 && p4jets[kk].Pt()<=240 ) { number_of_b++; hjets["pt_b_mc"]->Fill( p4jets[kk].Pt(), PUweight );}
872     if ( abs(listflavor[kk])==4 && p4jets[kk].Pt()<=240 ) { number_of_c++; hjets["pt_c_mc"]->Fill( p4jets[kk].Pt(), PUweight );}
873 algomez 1.1 if ( abs(listflavor[kk])==1 || abs(listflavor[kk])==2 || abs(listflavor[kk])==3 || abs(listflavor[kk])==21 )
874 algomez 1.2 { number_of_l++; hjets["pt_l_mc"]->Fill( p4jets[kk].Pt(), PUweight );}
875     if ( abs(listflavor[kk])==5 && p4jets[kk].Pt()>240 ) { number_of_b_highpt++; hjets["pt_b_mc"]->Fill( p4jets[kk].Pt(), PUweight );}
876     if ( abs(listflavor[kk])==4 && p4jets[kk].Pt()>240 ) { number_of_c_highpt++; hjets["pt_c_mc"]->Fill( p4jets[kk].Pt(), PUweight );}
877 algomez 1.1
878     if ( isTagb["TCHPM"][kk] )
879     {
880     hjets["pt_btag"]->Fill( p4jets[kk].Pt(), PUweight );
881     if ( abs(listflavor[kk])==5 ) hjets["pt_btag_b"]->Fill( p4jets[kk].Pt(), PUweight );
882     if ( abs(listflavor[kk])==4 ) hjets["pt_btag_c"]->Fill( p4jets[kk].Pt(), PUweight );
883     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 );
884     }
885     }
886    
887     hPVs["Nreweight_2jet"]->Fill( total_pvs, PUweight );
888     hjets["1st_pt"]->Fill( p4jets[0].Pt(), PUweight );
889     hjets["2nd_pt"]->Fill( p4jets[1].Pt(), PUweight );
890     hjets["1st_eta"]->Fill( p4jets[0].Eta(), PUweight );
891     hjets["2nd_eta"]->Fill( p4jets[1].Eta(), PUweight );
892    
893     hmuons["pt_2jet"]->Fill( p4lepton.Pt(), PUweight );
894     hmuons["reliso_2jet"]->Fill( RelIso, PUweight );
895     hMET["MET_2jet"]->Fill( p4MET.Pt(), PUweight );
896 algomez 1.2
897 algomez 1.1 if (fIsMC) {
898     hMET["genMET_2jet"]->Fill( ntuple->gen.MET, PUweight );
899     hMET["deltaMET_2jet"]->Fill( p4MET.Pt() - ntuple->gen.MET, PUweight );
900     }
901    
902     hM["WMt_2jet"]->Fill( WMt, PUweight );
903    
904     TLorentzVector p4Dijet = p4jets[0] + p4jets[1];
905     double Dijet_deltaPhi = p4jets[0].DeltaPhi(p4jets[1]);
906     TLorentzVector p4Top = p4jets[2] + p4LepW;
907    
908 algomez 1.3 //TLorentzVector p4Wprime = p4LepW + p4Dijet;
909     //hM["Wprime"]->Fill( p4Wprime.M(), PUweight );
910 algomez 1.1
911 algomez 1.3 //if (fVerbose) cout << "Wprime mass calculated" << endl;
912 algomez 1.1
913     // count number of b-tags
914     int Nbtags_TCHPM = 0;
915 algomez 1.2 int Nbtags_CSVM = 0;
916 algomez 1.1 for ( size_t itag=0; itag< isTagb["TCHPM"].size(); ++itag )
917     {
918     if ( isTagb["TCHPM"][itag] ) Nbtags_TCHPM++;
919 algomez 1.2 if ( isTagb["CSVM"][itag] ) Nbtags_CSVM++;
920 algomez 1.1 }
921    
922     // compute b-tag event weight
923     if ( fIsMC )
924     {
925    
926     // zeto tag
927     BTagWeight b0(0,0); // number of tags
928     BTagWeight::JetInfo bj(0.63,0.91); // mean MC eff and mean SF. For TCHPM=0.91\pm0.09, CSVM=0.96\pm0.096
929     BTagWeight::JetInfo cj(0.15,0.91);
930     double light_mceff = 0.017;
931     if ( 100 < p4jets[0].Pt() && p4jets[0].Pt() <= 200 ) light_mceff = 0.04;
932     if ( 200 < p4jets[0].Pt() && p4jets[0].Pt() <= 300 ) light_mceff = 0.08;
933     if ( 300 < p4jets[0].Pt() && p4jets[0].Pt() <= 400 ) light_mceff = 0.12;
934     if ( 400 < p4jets[0].Pt() ) light_mceff = 0.14;
935     BTagWeight::JetInfo lj(light_mceff,1.22); //for TCHPM=1.22, CSVM=1.08 \pm 0.13
936     // b-tag systematic UP 9% for b, 18% for c
937     BTagWeight::JetInfo bjUP(0.63,0.99);
938     BTagWeight::JetInfo cjUP(0.15,1.07);
939     // b-tag systemacit DOWN 9% for b, 18% for c
940     BTagWeight::JetInfo bjDOWN(0.63,0.83);
941     BTagWeight::JetInfo cjDOWN(0.15,0.75);
942     // for high pt jets > 240 UP 50% for b and c
943     BTagWeight::JetInfo bjUPhighpt(0.63,1.36);
944     BTagWeight::JetInfo cjUPhighpt(0.15,1.36);
945     // for high pt jets > 240 UP 50% for b and c
946     BTagWeight::JetInfo bjDOWNhighpt(0.63,0.46);
947     BTagWeight::JetInfo cjDOWNhighpt(0.15,0.46);
948    
949     vector<BTagWeight::JetInfo> j;
950     for(int i=0;i<number_of_b;i++)j.push_back(bj);
951     for(int i=0;i<number_of_b_highpt;i++)j.push_back(bj);
952     for(int i=0;i<number_of_c;i++)j.push_back(cj);
953     for(int i=0;i<number_of_c_highpt;i++)j.push_back(cj);
954     for(int i=0;i<number_of_l;i++)j.push_back(lj);
955    
956     if (Nbtags_TCHPM==0) {
957     SFb_0tag = b0.weight(j,0);
958     hjets["Nbtags_TCHPM"]->Fill( Nbtags_TCHPM, PUweight*SFb_0tag ); // fill bin 0
959     }
960     // only one tag
961     BTagWeight b11(1,1); // number of tags
962     if (Nbtags_TCHPM==1) {
963     SFb_only1tag = b11.weight(j,1);
964     hjets["Nbtags_TCHPM"]->Fill( Nbtags_TCHPM, PUweight*SFb_only1tag ); // fill bin 1
965     }
966     // at least one tag
967     BTagWeight b1(1,Nbtags_TCHPM); // number of tags
968     if (Nbtags_TCHPM>=1) {
969     SFb_1tag = b1.weight(j,1);
970    
971     vector<BTagWeight::JetInfo> jj;
972     for(int i=0;i<number_of_b;i++)jj.push_back(bjUP);
973     for(int i=0;i<number_of_b_highpt;i++)jj.push_back(bjUPhighpt);
974     for(int i=0;i<number_of_c;i++)jj.push_back(cjUP);
975     for(int i=0;i<number_of_c_highpt;i++)jj.push_back(cjUPhighpt);
976     for(int i=0;i<number_of_l;i++)jj.push_back(lj);
977     SFb_1tag_syst[0] = b1.weight(jj,1); //UP
978     vector<BTagWeight::JetInfo> jk;
979     for(int i=0;i<number_of_b;i++)jk.push_back(bjDOWN);
980     for(int i=0;i<number_of_b_highpt;i++)jk.push_back(bjDOWNhighpt);
981     for(int i=0;i<number_of_c;i++)jk.push_back(cjDOWN);
982     for(int i=0;i<number_of_c_highpt;i++)jk.push_back(cjDOWNhighpt);
983     for(int i=0;i<number_of_l;i++)jk.push_back(lj);
984     SFb_1tag_syst[1] = b1.weight(jk,1);//DOWN
985    
986     }
987     // at least two tags
988     BTagWeight b2(2,Nbtags_TCHPM); // number of tags
989     if (Nbtags_TCHPM>=2) {
990     SFb_2tag = b2.weight(j,2);
991     hjets["Nbtags_TCHPM"]->Fill( Nbtags_TCHPM, PUweight*SFb_2tag ); // fill bin >=2
992    
993     vector<BTagWeight::JetInfo> jj;
994     for(int i=0;i<number_of_b;i++)jj.push_back(bjUP);
995     for(int i=0;i<number_of_b_highpt;i++)jj.push_back(bjUPhighpt);
996     for(int i=0;i<number_of_c;i++)jj.push_back(cjUP);
997     for(int i=0;i<number_of_c_highpt;i++)jj.push_back(cjUPhighpt);
998     for(int i=0;i<number_of_l;i++)jj.push_back(lj);
999     SFb_2tag_syst[0] = b2.weight(jj,2); //UP
1000    
1001     vector<BTagWeight::JetInfo> jk;
1002     for(int i=0;i<number_of_b;i++)jk.push_back(bjDOWN);
1003     for(int i=0;i<number_of_b_highpt;i++)jk.push_back(bjDOWNhighpt);
1004     for(int i=0;i<number_of_c;i++)jk.push_back(cjDOWN);
1005     for(int i=0;i<number_of_c_highpt;i++)jk.push_back(cjDOWNhighpt);
1006     for(int i=0;i<number_of_l;i++)jk.push_back(lj);
1007     SFb_2tag_syst[1] = b2.weight(jk,2);//DOWN
1008     }
1009     }
1010     else hjets["Nbtags_TCHPM"]->Fill( Nbtags_TCHPM );
1011    
1012     //hjets["Nbtags_TCHPM"]->Fill( Nbtags_TCHPM, PUweight*SFb );
1013 algomez 1.2 //hjets["Nbtags_CSVM"]->Fill( Nbtags_CSVM, PUweight*SFb );
1014 algomez 1.1
1015 algomez 1.2 if ( Nbtags_TCHPM == 0 && njets < 5)
1016 algomez 1.1 {
1017 algomez 1.2 // compute top mass
1018     int top_bjet_index = 0;
1019     double deltaM = 99999.;
1020     //double top_mass = 0.;
1021     //TLorentzVector p4Top;
1022     TLorentzVector bestp4Top;
1023     for ( size_t itejet = 0; itejet < p4jets.size(); ++itejet )
1024     {
1025    
1026     for ( int isolw =0; isolw<2; ++ isolw)
1027     {
1028    
1029     p4Top = p4jets[itejet] + p4LepW;
1030     if (isolw==1) p4Top = p4jets[itejet] + p4OtherLepW;
1031    
1032     if ( fabs( p4Top.M() - 172.) < deltaM ) {
1033     top_bjet_index = itejet;
1034     //top_mass = p4Top.M();
1035     deltaM = fabs( p4Top.M() - 172.);
1036     bestp4Top = p4Top;
1037     }
1038    
1039     }
1040     }
1041     p4Top = bestp4Top;
1042 algomez 1.3 /*bool passcutWlep = false;
1043 algomez 1.2 if ( p4LepW.M() < 90 ) passcutWlep = true;
1044    
1045     if (passcutWlep) {
1046    
1047     if ( top_bjet_index == 0 ) {
1048     p4Wprime = p4Top + p4jets[1];
1049 algomez 1.1 }
1050 algomez 1.2 else {
1051     p4Wprime = p4Top + p4jets[0];
1052 algomez 1.1 }
1053 algomez 1.2
1054     bool passcut = true;
1055     if (fdoMtopCut && ( p4Top.M() <= 130 || p4Top.M() >= 210 ) ) passcut = false;
1056    
1057     if (passcut) {
1058 algomez 1.1
1059 algomez 1.2 hM["Wprime_0btag"]->Fill( p4Wprime.M(), PUweight*SFb_0tag );
1060     if (fSample=="WJets")
1061     {
1062     int FH = ntuple->flavorHistory;
1063     if ( FH == 1 || FH == 2 || FH == 5 || FH == 7 || FH == 9 ) hM["Wprime_0btag_bb"]->Fill( p4Wprime.M(), PUweight*SFb_0tag );
1064     else if ( FH == 3 || FH == 4 || FH == 6 || FH == 8 || FH == 10) hM["Wprime_0btag_cc"]->Fill( p4Wprime.M(), PUweight*SFb_0tag );
1065     else if ( FH == 11 ) hM["Wprime_0btag_light"]->Fill( p4Wprime.M(), PUweight*SFb_0tag );
1066     }
1067     if (fSample.Contains("Wprime"))
1068     {
1069     if ( ntuple->gen.LeptonicChannel == 11 ) hM["Wprime_0btag_MCsemie"]->Fill( p4Wprime.M(), PUweight*SFb_0tag );
1070     else if ( ntuple->gen.LeptonicChannel == 13 ) hM["Wprime_0btag_MCsemimu"]->Fill( p4Wprime.M(), PUweight*SFb_0tag );
1071     else hM["Wprime_0btag_MChad"]->Fill( p4Wprime.M(), PUweight*SFb_0tag );
1072     }
1073     } // top mass cut
1074 algomez 1.3 }// Wlep cut */
1075     }// 0 tag
1076 algomez 1.1
1077     // check if the two leading jets have a b jet
1078     if ( Nbtags_TCHPM >= 1 && (isTagb["TCHPM"][0] || isTagb["TCHPM"][1]) )
1079     {
1080    
1081     //cutmap["2Jet1b"] += PUweight*SFb_1tag;
1082     hjets["Njets_1tag"]->Fill( njets, PUweight*SFb_1tag );
1083     if (Nbtags_TCHPM > 1 ) hjets["Njets_2tag"]->Fill( njets, PUweight*SFb_2tag );
1084    
1085     if ( njets < 5 ) {
1086    
1087     //cutmap["2Jet1b"] += PUweight*SFb_1tag;
1088     // calculate dijet mass closest to W mass
1089     double the_dijet_mass = 0;
1090     double sigma2 = 13.*13.;
1091     double tmpchi2 = 999999.;
1092     for ( size_t itejet = 0; itejet < p4jets.size(); ++itejet )
1093     {
1094     for ( size_t jtejet = itejet+1; jtejet < p4jets.size(); ++jtejet )
1095     {
1096     TLorentzVector tmpvv = p4jets[itejet] + p4jets[jtejet];
1097     double tmpmass = tmpvv.M();
1098     if ( (tmpmass - 83.8 )*(tmpmass - 83.8)/sigma2 < tmpchi2 ) { tmpchi2 = (tmpmass - 83.8 )*(tmpmass - 83.8)/sigma2; the_dijet_mass = tmpmass; }
1099     }
1100     }
1101     hM["dijet"]->Fill( the_dijet_mass, PUweight*SFb_1tag );
1102     hjets["Dijet_deltaPhi"]->Fill( Dijet_deltaPhi, PUweight*SFb_1tag );
1103    
1104     // compute jjj mass for events with >=3 jets
1105     if ( njets > 2 ) {
1106     TLorentzVector p4Trijet = p4jets[0] + p4jets[1] + p4jets[2];
1107     hM["trijet"]->Fill( p4Trijet.M(), PUweight*SFb_1tag );
1108     }
1109     // compute top mass
1110     int top_bjet_index = 0;
1111     double deltaM = 99999.;
1112     //double top_mass = 0.;
1113     //TLorentzVector p4Top;
1114     TLorentzVector bestp4Top;
1115     for ( size_t itejet = 0; itejet < p4jets.size(); ++itejet )
1116     {
1117    
1118     for ( int isolw =0; isolw<2; ++ isolw)
1119     {
1120    
1121     p4Top = p4jets[itejet] + p4LepW;
1122     if (isolw==1) p4Top = p4jets[itejet] + p4OtherLepW;
1123    
1124     if ( fabs( p4Top.M() - 172.) < deltaM ) {
1125     top_bjet_index = itejet;
1126     //top_mass = p4Top.M();
1127     deltaM = fabs( p4Top.M() - 172.);
1128     bestp4Top = p4Top;
1129     }
1130    
1131     }
1132     }
1133     p4Top = bestp4Top;
1134     hM["top_1btag"]->Fill( p4Top.M(), PUweight*SFb_1tag );
1135     hMET["LepWmass"]->Fill(p4LepW.M(), PUweight*SFb_1tag );
1136     hMET["LepWmassNoPt"]->Fill(WmassNoPt, PUweight*SFb_1tag );
1137    
1138 algomez 1.3 /*bool passcutWlep = false;
1139 algomez 1.2 if ( p4LepW.M() < 90 ) passcutWlep = true;
1140 algomez 1.1
1141     if (passcutWlep) {
1142    
1143     if ( p4Top.M() > 130 && p4Top.M() < 210 ) hMET["LepWmass_topcut"]->Fill(p4LepW.M(), PUweight*SFb_1tag );
1144    
1145     double tb_deltaPhi = 0.;
1146     double tb_deltaR = 0.;
1147 algomez 1.2 double tb_deltaEta = -999;
1148     double pt_wprime = 0;
1149     double pt_top = 0;
1150     double pt_b = 0;
1151 algomez 1.1 if ( top_bjet_index == 0 ) {
1152     p4Wprime = p4Top + p4jets[1];
1153     tb_deltaPhi = p4Top.DeltaPhi( p4jets[1] );
1154     tb_deltaR = p4Top.DeltaR( p4jets[1] );
1155 algomez 1.2 tb_deltaEta = p4Top.Eta() - p4jets[1].Eta();
1156     pt_wprime = p4Wprime.Pt();
1157     pt_top = p4Top.Pt();
1158     pt_b = p4jets[1].Pt();
1159 algomez 1.1 }
1160     else {
1161     p4Wprime = p4Top + p4jets[0];
1162     tb_deltaPhi = fabs( p4Top.DeltaPhi( p4jets[0] ) );
1163     tb_deltaR = p4Top.DeltaR( p4jets[0] );
1164 algomez 1.2 tb_deltaEta = p4Top.Eta() - p4jets[0].Eta();
1165     pt_wprime= p4Wprime.Pt();
1166     pt_top = p4Top.Pt();
1167     pt_b = p4jets[0].Pt();
1168 algomez 1.1 }
1169    
1170     bool passcut = true;
1171     if (fdoMtopCut && ( p4Top.M() <= 130 || p4Top.M() >= 210 ) ) passcut = false;
1172    
1173     if (passcut) {
1174     cutmap["2Jet1b"] += PUweight*SFb_1tag;
1175     h2_pt_Wprime->Fill( p4Top.Pt(), p4Wprime.M(), PUweight*SFb_1tag );
1176     hM["Wprime_1btag"]->Fill( p4Wprime.M(), PUweight*SFb_1tag );
1177 algomez 1.2 // FINAL Plots
1178 algomez 1.1 hMET["deltaPhi"]->Fill( p4lepton.DeltaPhi( p4MET ), PUweight*SFb_1tag );
1179 algomez 1.2 hjets["tb_deltaPhi"]->Fill( tb_deltaPhi, PUweight*SFb_1tag );
1180     hjets["tb_deltaEta"]->Fill( tb_deltaEta, PUweight*SFb_1tag );
1181     //hjets["tb_deltaR"]->Fill( tb_deltaR, PUweight*SFb_1tag );
1182     hjets["pt_Wprime"]->Fill( pt_wprime, PUweight*SFb_1tag );
1183     hjets["pt_top"]->Fill( pt_top, PUweight*SFb_1tag );
1184     hjets["pt_b"]->Fill( pt_b, PUweight*SFb_1tag );
1185     hjets["1st_pt_fin"]->Fill( p4jets[0].Pt(), PUweight*SFb_1tag );
1186     hjets["2nd_pt_fin"]->Fill( p4jets[1].Pt(), PUweight*SFb_1tag );
1187     hjets["1st_eta_fin"]->Fill( p4jets[0].Eta(), PUweight*SFb_1tag );
1188     hjets["2nd_eta_fin"]->Fill( p4jets[1].Eta(), PUweight*SFb_1tag );
1189     hmuons["pt_fin"]->Fill( p4lepton.Pt(), PUweight*SFb_1tag );
1190     hmuons["reliso_fin"]->Fill( RelIso, PUweight*SFb_1tag );
1191     hMET["MET_fin"]->Fill( p4MET.Pt(), PUweight*SFb_1tag );
1192     hM["WMt_fin"]->Fill( WMt, PUweight*SFb_1tag );
1193 algomez 1.1
1194 algomez 1.2 if ( p4Wprime.M() > 1500.)
1195 algomez 1.1 {
1196     TString outstring = "run: ";
1197     outstring += TString(Form("%i",ntuple->run)) +" lumi: "+ TString(Form("%i",ntuple->lumi)) + " event: " + TString(Form("%i",ntuple->event));
1198     Info("Process",outstring);
1199     }
1200    
1201     if (fIsMC)
1202     {
1203     hM["Wprime_1btag_systUp"]->Fill( p4Wprime.M(), PUweight*SFb_1tag_syst[0] );
1204     hM["Wprime_1btag_systDown"]->Fill( p4Wprime.M(), PUweight*SFb_1tag_syst[1] );
1205     }
1206    
1207     if (fSample.Contains("Wprime"))
1208     {
1209     if ( ntuple->gen.LeptonicChannel == 11 ) hM["Wprime_1btag_MCsemie"]->Fill( p4Wprime.M(), PUweight*SFb_1tag );
1210     else if ( ntuple->gen.LeptonicChannel == 13 ) hM["Wprime_1btag_MCsemimu"]->Fill( p4Wprime.M(), PUweight*SFb_1tag );
1211     else hM["Wprime_1btag_MChad"]->Fill( p4Wprime.M(), PUweight*SFb_1tag );
1212     }
1213    
1214     if (fSample=="WJets")
1215     {
1216     int FH = ntuple->flavorHistory;
1217     if ( FH == 1 || FH == 2 || FH == 5 || FH == 7 || FH == 9 ) hM["Wprime_1btag_bb"]->Fill( p4Wprime.M(), PUweight*SFb_1tag );
1218     else if ( FH == 3 || FH == 4 || FH == 6 || FH == 8 || FH == 10) hM["Wprime_1btag_cc"]->Fill( p4Wprime.M(), PUweight*SFb_1tag );
1219     else if ( FH == 11 ) hM["Wprime_1btag_light"]->Fill( p4Wprime.M(), PUweight*SFb_1tag );
1220     }
1221    
1222     // check two b-tags
1223     if ( Nbtags_TCHPM > 1 )
1224     {
1225     cutmap["2Jet2b"] += PUweight*SFb_2tag;
1226     hM["Wprime_2btag"]->Fill( p4Wprime.M(), PUweight*SFb_2tag );
1227    
1228     if (fIsMC)
1229     {
1230     hM["Wprime_2btag_systUp"]->Fill( p4Wprime.M(), PUweight*SFb_2tag_syst[0] );
1231     hM["Wprime_2btag_systDown"]->Fill( p4Wprime.M(), PUweight*SFb_2tag_syst[1] );
1232     }
1233     }
1234    
1235     }//passcut mtop cut if requested
1236 algomez 1.3 }//passcutWlep cut */
1237 algomez 1.1 }// njets < 5
1238    
1239     }
1240    
1241 algomez 1.3 // if ( Nbtags_TCHPM == 1 ) hM["Wprime_1onlybtag"]->Fill( p4Wprime.M(), PUweight*SFb_only1tag );
1242 algomez 1.1
1243     }
1244    
1245    
1246     if (fVerbose) cout << "done analysis" << endl;
1247    
1248     return kTRUE;
1249     }
1250    
1251     void Analyzer::SlaveTerminate()
1252     {
1253     // The SlaveTerminate() function is called after all entries or objects
1254     // have been processed. When running with PROOF SlaveTerminate() is called
1255     // on each slave server.
1256    
1257     // fill cutflow histogram
1258    
1259     int ibin = 1;
1260     for ( vector<string>::const_iterator ivec= fCutLabels.begin(); ivec != fCutLabels.end(); ++ivec )
1261     // for ( map<string, int >::const_iterator imap=cutmap.begin(); imap!=cutmap.end(); ++imap )
1262     {
1263     hcutflow->SetBinContent( ibin, cutmap[ *ivec ] );
1264     ibin++;
1265     }
1266     // Write the ntuple to the file
1267     if (fFile) {
1268     Bool_t cleanup = kFALSE;
1269     TDirectory *savedir = gDirectory;
1270     if ( h1test->GetEntries() > 0) {
1271     fFile->cd();
1272     h1test->Write();
1273     hcutflow->Write();
1274 algomez 1.3 //h2_pt_Wprime->Write();
1275 algomez 1.1 fFile->mkdir("muons");
1276     fFile->cd("muons");
1277     for ( map<string,TH1* >::const_iterator imap=hmuons.begin(); imap!=hmuons.end(); ++imap )
1278     {
1279     TH1 *temp = imap->second;
1280     if ( temp->GetEntries() > 0 )
1281     temp->Write();
1282     //else cout << "Warning: empty histogram " << temp->GetName() << " will not be written to file." << endl;
1283     }
1284     fFile->cd();
1285     fFile->mkdir("PVs");
1286     fFile->cd("PVs");
1287     for ( map<string,TH1* >::const_iterator imap=hPVs.begin(); imap!=hPVs.end(); ++imap )
1288     {
1289     TH1 *temp = imap->second;
1290     if ( temp->GetEntries() > 0 )
1291     temp->Write();
1292     //else cout << "Warning: empty histogram " << temp->GetName() << " will not be written to file." << endl;
1293     }
1294     fFile->cd();
1295     fFile->mkdir("jets");
1296     fFile->cd("jets");
1297     for ( map<string,TH1* >::const_iterator imap=hjets.begin(); imap!=hjets.end(); ++imap )
1298     {
1299     TH1 *temp = imap->second;
1300     if ( temp->GetEntries() > 0 )
1301     temp->Write();
1302     }
1303     fFile->cd();
1304     fFile->mkdir("mass");
1305     fFile->cd("mass");
1306     for ( map<string,TH1* >::const_iterator imap=hM.begin(); imap!=hM.end(); ++imap )
1307     {
1308     TH1 *temp = imap->second;
1309     if ( temp->GetEntries() > 0 )
1310     temp->Write();
1311     }
1312     fFile->cd();
1313     fFile->mkdir("MET");
1314     fFile->cd("MET");
1315     for ( map<string,TH1* >::const_iterator imap=hMET.begin(); imap!=hMET.end(); ++imap )
1316     {
1317     TH1 *temp = imap->second;
1318     if ( temp->GetEntries() > 0 )
1319     temp->Write();
1320     }
1321    
1322     fFile->cd();
1323    
1324     fProofFile->Print();
1325     fOutput->Add(fProofFile);
1326     } else {
1327     cleanup = kTRUE;
1328     }
1329    
1330 algomez 1.3 //h2_pt_Wprime->SetDirectory(0);
1331 algomez 1.1
1332     h1test->SetDirectory(0);
1333     hcutflow->SetDirectory(0);
1334     gDirectory = savedir;
1335     fFile->Close();
1336     // Cleanup, if needed
1337     if (cleanup) {
1338     TUrl uf(*(fFile->GetEndpointUrl()));
1339     SafeDelete(fFile);
1340     gSystem->Unlink(uf.GetFile());
1341     SafeDelete(fProofFile);
1342     }
1343     }
1344    
1345     }
1346    
1347     void Analyzer::Terminate()
1348     {
1349     // The Terminate() function is the last function to be called during
1350     // a query. It always runs on the client, it can be used to present
1351     // the results graphically or save the results to file.
1352    
1353     Info("Terminate","Analyzer done.");
1354     }