ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/StudyModule.C
Revision: 1.5
Committed: Fri Mar 2 08:36:07 2012 UTC (13 years, 2 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.4: +16 -4 lines
Log Message:
Added NP/FP signal distribution functions

File Contents

# User Rev Content
1 buchmann 1.1 #include <iostream>
2     #include <vector>
3     #include <sys/stat.h>
4    
5     #include <TCut.h>
6     #include <TROOT.h>
7     #include <TCanvas.h>
8     #include <TMath.h>
9     #include <TColor.h>
10     #include <TPaveText.h>
11     #include <TRandom.h>
12     #include <TH1.h>
13     #include <TH2.h>
14     #include <TF1.h>
15     #include <TSQLResult.h>
16     #include <TProfile.h>
17     #include <TPaveStats.h>
18    
19 buchmann 1.4 #include "GeneratorLevelStudyModule.C"
20    
21 buchmann 1.1 //#include "TTbar_stuff.C"
22     using namespace std;
23    
24     using namespace PlottingSetup;
25    
26    
27     void do_experimental_pred_obs_calculation(float cut ,string mcjzb,string datajzb, int mcordata) {
28     dout << "Crunching the numbers for JZB>" << cut << endl;
29     string xlabel="JZB [GeV] -- for algoritm internal use only!";
30     TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSSF&&cutnJets&&basiccut,mcordata,luminosity);
31     TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSOF&&cutnJets&&basiccut,mcordata,luminosity);
32     TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSSF&&cutnJets&&basiccut,mcordata,luminosity);
33     TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSOF&&cutnJets&&basiccut,mcordata,luminosity);
34    
35     TH1F *SBOSSFP;
36     TH1F *SBOSOFP;
37     TH1F *SBOSSFN;
38     TH1F *SBOSOFN;
39    
40     flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
41     if(PlottingSetup::RestrictToMassPeak) {
42     SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,1,cut,14000, xlabel, "events",cutOSSF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
43     SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,1,cut,14000, xlabel, "events",cutOSOF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
44     SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,1,cut,14000, xlabel, "events",cutOSSF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
45     SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,1,cut,14000, xlabel, "events",cutOSOF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
46     }
47    
48    
49     flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
50     if(PlottingSetup::RestrictToMassPeak) {
51     dout << " Observed : " << ZOSSFP->Integral() << endl;
52     dout << " Predicted: " << ZOSSFN->Integral() << " + (1/3)*(" << ZOSOFP->Integral() << "-" << ZOSOFN->Integral()<<") + (1/3)*(" << SBOSSFP->Integral() << "-" << SBOSSFN->Integral()<<") + (1/3)*(" << SBOSOFP->Integral() << "-" << SBOSOFN->Integral()<<")" << endl;
53     dout << " P(ZJets ) \t " << ZOSSFN->Integral() << endl;
54     dout << " P(e&mu;]) \t " << ZOSOFP->Integral() << "-" << ZOSOFN->Integral() << " = " << ZOSOFP->Integral()-ZOSOFN->Integral()<<endl;
55     dout << " P(ossf,sb]) \t " << SBOSSFP->Integral() << "-" << SBOSSFN->Integral()<<" = "<<SBOSSFP->Integral()-SBOSSFN->Integral()<<endl;
56     dout << " P(osof,sb]) \t " << SBOSOFP->Integral() << "-" << SBOSOFN->Integral()<<" = "<<SBOSOFP->Integral()-SBOSOFN->Integral()<<endl;
57     } else {
58     dout << " Observed : " << ZOSSFP->Integral() << endl;
59     dout << " Predicted: " << ZOSSFN->Integral() << " + (" << ZOSOFP->Integral() << "-" << ZOSOFN->Integral()<<")" << endl;
60     dout << " P(ZJets ) \t " << ZOSSFN->Integral() << endl;
61     dout << " P(e&mu;]) \t " << ZOSOFP->Integral() << "-" << ZOSOFN->Integral() << " = " << ZOSOFP->Integral()-ZOSOFN->Integral()<<endl;
62     }
63    
64    
65     delete ZOSSFP;
66     delete ZOSOFP;
67     delete ZOSSFN;
68     delete ZOSOFN;
69    
70     delete SBOSSFP;
71     delete SBOSOFP;
72     delete SBOSSFN;
73     delete SBOSOFN;
74     }
75    
76     void look_at_sidebands(string mcjzb, string datajzb, bool includejetcut, float cutat=0) {
77    
78     flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak -- this funciton is meaningless for the offpeak case
79     if(!PlottingSetup::RestrictToMassPeak) return;
80     dout << "Looking at sidebands ... " << endl;
81     int mcordata=data;//data // you can perform this study for mc or data ...
82    
83     TCut specialjetcut;
84     if(includejetcut) specialjetcut=cutnJets&&basiccut;
85     else specialjetcut="mll>0";
86     dout << "The datajzb variable is defined as " << datajzb << endl;
87     stringstream addcut;
88     addcut<<"(pt>"<<cutat<<")";
89     TCut additionalcut=addcut.str().c_str();
90    
91     int nbins=75; float min=51;float max=201;string xlabel="mll [GeV]";
92     TCanvas *c1 = new TCanvas("c1","c1");
93     c1->SetLogy(1);
94     TH1F *datahistoOSSF = allsamples.Draw("datahistoOSSF","mll",nbins,min,max, xlabel, "events",cutOSSF&&specialjetcut&&additionalcut,data,luminosity);
95     THStack mcstackOSSF = allsamples.DrawStack("mcstackOSSF","mll",nbins,min,max, xlabel, "events",cutOSSF&&specialjetcut&&additionalcut,mc,luminosity);
96    
97     datahistoOSSF->SetMinimum(1);
98     datahistoOSSF->Draw();
99     mcstackOSSF.Draw("same");
100     datahistoOSSF->Draw("same");
101     TLegend *kinleg = allsamples.allbglegend();
102     kinleg->AddEntry(datahistoOSSF,"OSSF (data)","p");
103     kinleg->Draw();
104     if(includejetcut) CompleteSave(c1,"sidebands/"+any2string(cutat)+"/OSSF");
105     else CompleteSave(c1,"sidebands/"+any2string(cutat)+"/OSSF_nojetcut");
106    
107     TH1F *datahistoOSOF = allsamples.Draw("datahistoOSOF","mll",nbins,min,max, xlabel, "events",cutOSOF&&specialjetcut&&additionalcut,data,luminosity);
108     THStack mcstackOSOF = allsamples.DrawStack("mcstackOSOF","mll",nbins,min,max, xlabel, "events",cutOSOF&&specialjetcut&&additionalcut,mc,luminosity);
109     // datahistoOSOF->SetMinimum(0.4);
110     datahistoOSOF->Draw();
111     mcstackOSOF.Draw("same");
112     datahistoOSOF->Draw("same");
113     TLegend *kinleg2 = allsamples.allbglegend();
114     kinleg2->AddEntry(datahistoOSOF,"OSOF (data)","p");
115     kinleg2->Draw();
116     if(includejetcut) CompleteSave(c1,"sidebands/"+any2string(cutat)+"/OSOF");
117     else CompleteSave(c1,"sidebands/"+any2string(cutat)+"/OSOF_nojetcut");
118    
119     TH1F *rawmlleemmData = allsamples.Draw("rawmlleemmData","mll",200,0,200, "mll [GeV]", "events", cutOSSF&&specialjetcut&&sidebandcut&&additionalcut,data, luminosity);
120     TH1F *rawmllemData = allsamples.Draw("rawmllemData" ,"mll",200,0,200, "mll [GeV]", "events", cutmass&&cutOSOF&&specialjetcut&&additionalcut,data, luminosity);
121     dout << "Number of events in peak for OSOF: " << rawmllemData->GetEntries() << endl;
122     dout << "Number of events in SB for OSSF: " << rawmlleemmData->GetEntries() << endl;
123    
124     TH1F *SFttbarZpeak = allsamples.Draw("SFttbarZpeak",mcjzb,100,-200,400, "JZB [GeV]", "events",cutmass&&cutOSSF&&specialjetcut&&additionalcut,mc,luminosity,allsamples.FindSample("TTJets"));
125     TH1F *OFttbarZpeak = allsamples.Draw("OFttbarZpeak",mcjzb,100,-200,400, "JZB [GeV]", "events",cutmass&&cutOSOF&&specialjetcut&&additionalcut,mc,luminosity,allsamples.FindSample("TTJets"));
126     TH1F *SFttbarsideb = allsamples.Draw("SFttbarsideb",mcjzb,100,-200,400, "JZB [GeV]", "events",cutOSSF&&specialjetcut&&sidebandcut&&additionalcut,mc,luminosity,allsamples.FindSample("TTJets"));
127    
128     SFttbarZpeak->SetLineColor(kBlack);
129     OFttbarZpeak->SetLineColor(kBlue);
130     OFttbarZpeak->SetMarkerColor(kBlue);
131     SFttbarsideb->SetLineColor(kPink);
132     SFttbarsideb->SetMarkerColor(kPink);
133    
134     SFttbarZpeak->Draw("histo");
135     OFttbarZpeak->Draw("histo,same");
136     SFttbarsideb->Draw("histo,same");
137    
138     TLegend *leg3 = new TLegend(0.6,0.8,0.89,0.89);
139     leg3->AddEntry(SFttbarZpeak,"SF ttbar Z peak","l");
140     leg3->AddEntry(OFttbarZpeak,"OF ttbar Z peak","l");
141     leg3->AddEntry(SFttbarsideb,"SF ttbar SB","l");
142     leg3->SetFillColor(kWhite);
143     leg3->SetLineColor(kWhite);
144     leg3->SetBorderSize(0);
145    
146     leg3->Draw();
147     if(includejetcut) CompleteSave(c1,"sidebands/"+any2string(cutat)+"/ttbar_comparison");
148     else CompleteSave(c1,"sidebands/"+any2string(cutat)+"/ttbar_comparison_nojetcut");
149    
150    
151     c1->SetLogy(0);
152    
153     SFttbarsideb->SetFillColor(TColor::GetColor("#F5A9A9"));
154     SFttbarsideb->SetLineColor(TColor::GetColor("#F5A9A9"));
155     SFttbarsideb->SetFillStyle(3004);
156     OFttbarZpeak->SetFillColor(TColor::GetColor("#819FF7"));
157     OFttbarZpeak->SetLineColor(TColor::GetColor("#819FF7"));
158     OFttbarZpeak->SetFillColor(kBlue);
159     OFttbarZpeak->SetFillStyle(3005);
160    
161     OFttbarZpeak->Rebin(2);
162     SFttbarZpeak->Rebin(2);
163     SFttbarsideb->Rebin(2);
164     OFttbarZpeak->Divide(SFttbarZpeak);
165     SFttbarsideb->Divide(SFttbarZpeak);
166     OFttbarZpeak->GetYaxis()->SetRangeUser(0,5);
167     OFttbarZpeak->GetYaxis()->SetTitle("ratio");
168     TF1 *centralfitO = new TF1("centralfitO","pol1",-40,120);
169     TF1 *centralfit1 = new TF1("centralfit1","pol1",-200,400);
170     // TF1 *centralfitS = new TF1("centralfitS","pol1",-40,120);
171     SFttbarsideb->Fit(centralfitO,"R");
172     //OFttbarZpeak->Fit(centralfitO,"R");
173     centralfit1->SetParameters(centralfitO->GetParameters());
174     // SFttbarZpeak->Fit(centralfitS,"R");
175     OFttbarZpeak->Draw("e5");
176     SFttbarsideb->Draw("e5,same");
177     OFttbarZpeak->Draw("same");
178     SFttbarsideb->Draw("same");
179     centralfit1->SetLineColor(kOrange);
180     // centralfitS->SetLineColor(kOrange);
181     centralfit1->SetLineWidth(2);
182     // centralfitS->SetLineWidth(2);
183     centralfit1->Draw("same");
184     // centralfitS->Draw("same");
185     TLine *oneline = new TLine(-200,1,400,1);
186     oneline->SetLineColor(kBlue);
187     TLine *point5 = new TLine(-200,0.5,400,0.5);
188     point5->SetLineStyle(2);
189     point5->SetLineColor(kGreen);
190     TLine *op5 = new TLine(-200,1.5,400,1.5);
191     op5->SetLineStyle(2);
192     op5->SetLineColor(kGreen);
193     TLine *point7 = new TLine(-200,0.7,400,0.7);
194     point7->SetLineStyle(2);
195     point7->SetLineColor(kBlack);
196     TLine *op7 = new TLine(-200,1.3,400,1.3);
197     op7->SetLineStyle(2);
198     op7->SetLineColor(kBlack);
199     oneline->Draw("same");
200     point5->Draw("same");
201     point7->Draw("same");
202     op5->Draw("same");
203     op7->Draw("same");
204     TLegend *leg4 = new TLegend(0.6,0.65,0.89,0.89);
205     leg4->AddEntry(OFttbarZpeak,"OF ttbar Z peak / truth","l");
206     leg4->AddEntry(SFttbarsideb,"SF ttbar SB / truth","l");
207     leg4->AddEntry(centralfit1,"Fit to [-40,120] GeV region (OF)","l");
208     leg4->AddEntry(point5,"50% systematic envelope","l");
209     leg4->AddEntry(point7,"30% systematic envelope","l");
210     // leg4->AddEntry(centralfitS,"Fit to [-40,120] GeV region (SF)","l");
211     leg4->SetFillColor(kWhite);
212     leg4->SetLineColor(kWhite);
213     leg4->SetBorderSize(0);
214     leg4->Draw("same");
215     if(includejetcut) CompleteSave(c1,"sidebands/"+any2string(cutat)+"/ttbar_comparison_ratio");
216     else CompleteSave(c1,"sidebands/"+any2string(cutat)+"/ttbar_comparison_ratio_nojetcut");
217     dout << "Moving on to predicted / observed yields! " << endl;
218     dout << "Sideband definition: " << (const char*) sidebandcut << endl;
219     /*
220     do_experimental_pred_obs_calculation(50,mcjzb,datajzb,mcordata);
221     do_experimental_pred_obs_calculation(75,mcjzb,datajzb,mcordata);
222     do_experimental_pred_obs_calculation(100,mcjzb,datajzb,mcordata);
223     do_experimental_pred_obs_calculation(125,mcjzb,datajzb,mcordata);
224     do_experimental_pred_obs_calculation(150,mcjzb,datajzb,mcordata);
225     */
226    
227     delete rawmlleemmData;
228     delete rawmllemData;
229     delete SFttbarZpeak;
230     delete OFttbarZpeak;
231     delete SFttbarsideb;
232     delete datahistoOSOF;
233     delete datahistoOSSF;
234    
235    
236     }
237    
238     void look_at_sidebands(string mcjzb, string datajzb) {
239     // for (int i=0;i<100;i+=10) {
240     int i=0;
241     {
242     look_at_sidebands(mcjzb,datajzb, true,i);
243     look_at_sidebands(mcjzb,datajzb, false,i);
244     }
245     }
246    
247     void find_sideband_definition() {
248     flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
249     if(!PlottingSetup::RestrictToMassPeak) return; // this function is meaningless for the offpeak analysis
250    
251     TH1F *mllttbar = allsamples.Draw("mllttbar","mll",145,55,200, "mll [GeV]", "events",cutOSSF&&cutnJets&&!cutmass,mc,luminosity,allsamples.FindSample("TTJets"));
252     TH1F *mllttbarz = allsamples.Draw("mllttbarz","mll",1,50,200, "mll [GeV]", "events",cutOSSF&&cutnJets&&cutmass,mc,luminosity,allsamples.FindSample("TTJets"));
253     float leftstop=0;
254     float rightstop=0;
255     int pos90=mllttbar->FindBin(90);
256     float leftsum=0; float rightsum=0;
257     float peaksum=mllttbarz->Integral();
258     for(int i=0;i<mllttbar->GetNbinsX()&&!(leftstop&&rightstop);i++) {
259     if(pos90-i<1) leftstop=mllttbar->GetBinLowEdge(1);
260     if(mllttbar->GetBinLowEdge(pos90+i)+mllttbar->GetBinWidth(pos90+i)>190) rightstop=190;
261     if(!leftstop) leftsum+=mllttbar->GetBinContent(pos90-i);
262     if(!rightstop) rightsum+=mllttbar->GetBinContent(pos90+i);
263     if(leftsum+rightsum>peaksum) {
264     if(!leftstop) leftstop=mllttbar->GetBinLowEdge(pos90-i);
265     if(!rightstop) rightstop=mllttbar->GetBinLowEdge(pos90+i)+mllttbar->GetBinWidth(pos90+i);
266     dout << "Found the boundaries! on the left: " << leftstop << " and on the right " << rightstop << endl;
267     dout << "Total sum : " << leftsum+rightsum << " which supposedly corresponds ~ to " << peaksum << endl;
268     }
269     }
270     TH1F *mllttbart = allsamples.Draw("mllttbart","mll",1,55,155, "mll [GeV]", "events",cutOSSF&&cutnJets&&!cutmass,mc,luminosity,allsamples.FindSample("TTJets"));
271     dout << mllttbart->Integral() << endl;
272    
273    
274     }
275    
276     void calculate_upper_limits(string mcjzb, string datajzb) {
277     write_warning("calculate_upper_limits","Upper limit calculation temporarily deactivated");
278     // write_warning("calculate_upper_limits","Calculation of SUSY upper limits has been temporarily suspended in favor of top discovery");
279     // rediscover_the_top(mcjzb,datajzb);
280     /*
281     TCanvas *c3 = new TCanvas("c3","c3");
282     c3->SetLogy(1);
283     vector<float> binning;
284     //binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
285     binning.push_back(50);
286     binning.push_back(100);
287     binning.push_back(150);
288     binning.push_back(200);
289     binning.push_back(500);
290     TH1F *datapredictiona = allsamples.Draw("datapredictiona", "-"+datajzb, binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity);
291     TH1F *datapredictionb = allsamples.Draw("datapredictionb", "-"+datajzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity);
292     TH1F *datapredictionc = allsamples.Draw("datapredictionc", datajzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity);
293     TH1F *dataprediction = (TH1F*)datapredictiona->Clone();
294     dataprediction->Add(datapredictionb,-1);
295     dataprediction->Add(datapredictionc);
296     TH1F *puresignal = allsamples.Draw("puresignal", mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
297     TH1F *signalpred = allsamples.Draw("signalpred", "-"+mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
298     TH1F *signalpredlo = allsamples.Draw("signalpredlo", "-"+mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
299     TH1F *signalpredro = allsamples.Draw("signalpredro", mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
300     TH1F *puredata = allsamples.Draw("puredata", datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
301     signalpred->Add(signalpredlo,-1);
302     signalpred->Add(signalpredro);
303     puresignal->Add(signalpred,-1);//subtracting signal contamination
304     ofstream myfile;
305     myfile.open ("ShapeFit_log.txt");
306     establish_upper_limits(puredata,dataprediction,puresignal,"LM4",myfile);
307     myfile.close();
308     */
309     }
310    
311     TH1F *runcheckhisto(string cut) {
312     string histoname=GetNumericHistoName();
313     TH1F *histo = new TH1F(histoname.c_str(),histoname.c_str(),100,163000,168000);
314     (allsamples.collection)[0].events->Draw(("runNum>>"+histoname).c_str(),cut.c_str());
315     return histo;
316     }
317    
318     void run_check() {
319     gROOT->SetStyle("Plain");
320     TCanvas *c1 = new TCanvas("runnum","runnum",800,1000);
321     c1->Divide(2,4);
322     flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
323     if(PlottingSetup::RestrictToMassPeak) c1->Divide(2,2); // there are only four regions ...
324    
325     c1->cd(1);
326     TH1F *ossfp = runcheckhisto((const char*)(cutmass&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)"));
327     ossfp->Draw();
328     TText *t1 = write_title("OSSF,P");t1->Draw();
329    
330     c1->cd(2);
331     TH1F *ossfn = runcheckhisto((const char*)(cutmass&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)<-100)"));
332     ossfn->Draw();
333     TText *t2 = write_title("OSSF,N");t2->Draw();
334    
335     c1->cd(3);
336     TH1F *osofp = runcheckhisto((const char*)(cutmass&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)"));
337     osofp->Draw();
338     TText *t3 = write_title("OSOF,P");t3->Draw();
339     c1->cd(4);
340     TH1F *osofn = runcheckhisto((const char*)(cutmass&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)<-100)"));
341     osofn->Draw();
342     TText *t4 = write_title("OSOF,N");t4->Draw();
343    
344     flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
345     if(PlottingSetup::RestrictToMassPeak) {
346     c1->cd(5);
347     TH1F *sbofp = runcheckhisto((const char*)(sidebandcut&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)"));
348     sbofp->Draw();
349     TText *t5 = write_title("SB,OSOF,P");t5->Draw();
350     c1->cd(6);
351     TH1F *sbofn = runcheckhisto((const char*)(cutOSOF&&cutnJets&&basiccut&&sidebandcut&&"((jzb[1]+0.06*pt-2.84727)<-100)"));
352     sbofn->Draw();
353     TText *t6 = write_title("SB,OSOF,N");t6->Draw();
354    
355     c1->cd(7);
356     TH1F *sbsfp = runcheckhisto((const char*)(sidebandcut&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)"));
357     sbsfp->Draw();
358     TText *t7 = write_title("SB,OSSF,P");t7->Draw();
359     c1->cd(8);
360     TH1F *sbsfn = runcheckhisto((const char*)(sidebandcut&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)<-100)"));
361     sbsfn->Draw();
362     TText *t8 = write_title("SB,OSSF,N");t8->Draw();
363     }
364    
365     c1->SaveAs("runNumber.png");
366     }
367    
368     TH1F *give_boson_pred(TCut bcut,string mcjzb) {
369     int nbins=50;
370     TH1F *jzbn = allsamples.Draw("jzbn","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events", bcut&&cutOSSF&&cutnJets&&cutmass,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
371     TH1F *jzbno = allsamples.Draw("jzbno","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",bcut&&cutOSOF&&cutnJets&&cutmass,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
372     TH1F *jzbpo = allsamples.Draw("jzbp",mcjzb,nbins,0,350, "JZB [GeV]", "events", bcut&&cutOSOF&&cutnJets&&cutmass,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
373    
374     //Sidebands
375     TH1F *jzbnos;
376     TH1F *jzbpos;
377     TH1F *jzbnss;
378     TH1F *jzbpss;
379    
380     flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
381     if(PlottingSetup::RestrictToMassPeak) {
382     jzbnos = allsamples.Draw("jzbnos","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",bcut&&cutOSOF&&cutnJets&&sidebandcut,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
383     jzbpos = allsamples.Draw("jzbpos",mcjzb,nbins,0,350, "JZB [GeV]", "events", bcut&&cutOSOF&&cutnJets&&sidebandcut,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
384     jzbnss = allsamples.Draw("jzbnss","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",bcut&&cutOSSF&&cutnJets&&sidebandcut,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
385     jzbpss = allsamples.Draw("jzbpss",mcjzb,nbins,0,350, "JZB [GeV]", "events", bcut&&cutOSSF&&cutnJets&&sidebandcut,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
386     }
387    
388    
389     TH1F *pred = (TH1F*)jzbn->Clone("pred");
390     flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
391     if(PlottingSetup::RestrictToMassPeak) {
392     pred->Add(jzbno,-1.0/3);
393     pred->Add(jzbpo,1.0/3);
394     pred->Add(jzbnos,-1.0/3);
395     pred->Add(jzbpos,1.0/3);
396     pred->Add(jzbnss,-1.0/3);
397     pred->Add(jzbpss,1.0/3);
398     } else {
399     pred->Add(jzbno,-1.0);
400     pred->Add(jzbpo,1.0);
401     }
402     pred->SetLineColor(kRed);
403     pred->SetMinimum(0);
404     delete jzbn;
405     delete jzbpo;
406     delete jzbno;
407     delete jzbpos;
408     delete jzbnos;
409     delete jzbpss;
410     delete jzbnss;
411    
412     return pred;
413     }
414    
415    
416     void show_dibosons(string datajzb, string mcjzb) {
417     TCut WW("(abs(genMID1)==24&&abs(genMID2)==24)||(abs(genGMID1)==24&&abs(genGMID2)==24)");
418     TCut ZZ("(abs(genMID1)==23&&abs(genMID2)==23)||(abs(genGMID1)==23&&abs(genGMID2)==23)");
419     TCut WZ("((abs(genMID1)==23&&abs(genMID2)==24)||(abs(genGMID1)==23&&abs(genGMID2)==24))||((abs(genMID1)==24&&abs(genMID2)==23)||(abs(genGMID1)==24&&abs(genGMID2)==23))");
420    
421     TCanvas *dibs = new TCanvas("dibs","dibs",900,900);
422     dibs->Divide(2,2);
423    
424     dibs->cd(1);
425     TH1F *wwjzbp = allsamples.Draw("wwjzbp",mcjzb,70,0,350, "JZB [GeV]", "events", WW&&cutOSSF&&cutnJets&&cutmass,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
426     TH1F *wwpred = (TH1F*) give_boson_pred(WW,mcjzb);
427     wwpred->Draw("histo");
428     wwjzbp->Draw("histo,same");
429     TLegend *leg = make_legend("WW");
430     leg->SetFillColor(kWhite);
431     leg->SetLineColor(kWhite);
432     leg->SetHeader("WW (MC)");
433     leg->AddEntry(wwjzbp,"Observed","l");
434     leg->AddEntry(wwpred,"Predicted","l");
435     leg->Draw("same");
436    
437     dibs->cd(2);
438     TH1F *wzjzbp = allsamples.Draw("wzjzbp",mcjzb,70,0,350, "JZB [GeV]", "events",WZ&&cutOSSF&&cutnJets&&cutmass,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
439     TH1F *wzpred = (TH1F*) give_boson_pred(WZ,mcjzb);
440     wzpred->Draw("histo");
441     wzjzbp->Draw("same,histo");
442     TLegend *leg2 = (TLegend*)leg->Clone("leg2");
443     leg2->SetHeader("WZ (MC)");
444     leg2->Draw("same");
445     DrawPrelim();
446    
447     dibs->cd(3);
448     TH1F *zzjzbp = allsamples.Draw("zzjzbp",mcjzb,70,0,350, "JZB [GeV]", "events",ZZ&&cutOSSF&&cutnJets&&cutmass,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
449     TH1F *zzpred = (TH1F*) give_boson_pred(ZZ,mcjzb);
450     zzpred->Draw("histo");
451     zzjzbp->Draw("same,histo");
452     TLegend *leg3 = (TLegend*)leg->Clone("leg2");
453     leg3->SetHeader("ZZ (MC)");
454     leg3->Draw("same");
455     leg3->Draw("same");
456     DrawPrelim();
457    
458     dibs->cd(4);
459     TH1F *alljzbp = allsamples.Draw("alljzbp",mcjzb,70,0,350, "JZB [GeV]", "events",(WW||WZ||ZZ)&&cutOSSF&&cutnJets&&cutmass,mc,luminosity,allsamples.FindSample("VVJetsTo4L_TuneD6T_7TeV"));
460     TH1F *allpred = (TH1F*) give_boson_pred((WW||WZ||ZZ),mcjzb);
461     allpred->Draw("histo");
462     alljzbp->Draw("same,histo");
463     TLegend *leg4 = (TLegend*)leg->Clone("leg2");
464     leg4->SetHeader("All dibosons (MC)");
465     leg4->Draw("same");
466     DrawPrelim();
467    
468     CompleteSave(dibs,"Studies/Dibosons");
469    
470     }
471    
472     void show_rare_samples(string datajzb, string mcjzb) {
473     TCanvas *rares = new TCanvas("rares","Rare Samples",900,900);
474    
475     int nbins=50;
476     TH1F *jzbp = raresample.Draw("jzbp", mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSSF&&cutnJets&&cutmass,mc,luminosity);
477    
478     TH1F *jzbn = raresample.Draw("jzbn", "-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSSF&&cutnJets&&cutmass,mc,luminosity);
479     TH1F *jzbno = raresample.Draw("jzbno", "-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSOF&&cutnJets&&cutmass,mc,luminosity);
480     TH1F *jzbpo = raresample.Draw("jzbp", mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSOF&&cutnJets&&cutmass,mc,luminosity);
481    
482     //Sidebands
483     TH1F *jzbnos = raresample.Draw("jzbnos","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSOF&&cutnJets&&sidebandcut,mc,luminosity);
484     TH1F *jzbpos = raresample.Draw("jzbpos", mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSOF&&cutnJets&&sidebandcut,mc,luminosity);
485     TH1F *jzbnss = raresample.Draw("jzbnss","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSSF&&cutnJets&&sidebandcut,mc,luminosity);
486     TH1F *jzbpss = raresample.Draw("jzbpss", mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSSF&&cutnJets&&sidebandcut,mc,luminosity);
487    
488     TH1F *pred = (TH1F*)jzbn->Clone("pred");
489     pred->Add(jzbno,-1.0/3);
490     pred->Add(jzbpo,1.0/3);
491     pred->Add(jzbnos,-1.0/3);
492     pred->Add(jzbpos,1.0/3);
493     pred->Add(jzbnss,-1.0/3);
494     pred->Add(jzbpss,1.0/3);
495     pred->SetLineColor(kRed);
496     pred->SetMinimum(0);
497     delete jzbn;
498     delete jzbpo;
499     delete jzbno;
500     delete jzbpos;
501     delete jzbnos;
502     delete jzbpss;
503     delete jzbnss;
504    
505     pred->Draw("histo");
506     pred->GetYaxis()->SetTitleOffset(1.3);
507     jzbp->Draw("histo,same");
508     TLegend *leg = make_legend("WW");
509     leg->SetFillColor(kWhite);
510     leg->SetLineColor(kWhite);
511     leg->SetHeader("Rare Samples (MC)");
512     leg->AddEntry(jzbp,"Observed","l");
513     leg->AddEntry(pred,"Predicted","l");
514     leg->Draw("same");
515    
516    
517     CompleteSave(rares,"Studies/Rare_Samples");
518     delete jzbp;
519     delete pred;
520     delete leg;
521     delete rares;
522     }
523    
524    
525     class signature {
526     public:
527     int runNum;
528     int eventNum;
529     int lumi;
530     };
531    
532     vector<signature> get_list_of_events(string cut) {
533     float jzb;
534     int runNum,lumi,eventNum;
535     (allsamples.collection)[0].events->SetBranchAddress("jzb",&jzb);
536     (allsamples.collection)[0].events->SetBranchAddress("eventNum",&eventNum);
537     (allsamples.collection)[0].events->SetBranchAddress("lumi",&lumi);
538     (allsamples.collection)[0].events->SetBranchAddress("runNum",&runNum);
539    
540     TTreeFormula *select = new TTreeFormula("select", cut.c_str()&&essentialcut, (allsamples.collection)[0].events);
541     vector<signature> allevents;
542     for (Int_t entry = 0 ; entry < (allsamples.collection)[0].events->GetEntries() ; entry++) {
543     (allsamples.collection)[0].events->LoadTree(entry);
544     if (select->EvalInstance()) {
545     (allsamples.collection)[0].events->GetEntry(entry);
546     signature newevent;
547     newevent.runNum=runNum;
548     newevent.eventNum=eventNum;
549     newevent.lumi=lumi;
550     allevents.push_back(newevent);
551     }
552     }
553     cout << "Done looping!" << endl;
554     return allevents;
555     }
556    
557     void make_double_plot(string variable, int nbins, float min, float max, float ymax, bool logscale,string xlabel, string filename,TCut observed, TCut predicted, bool is_data, bool noscale = false) {
558     TCut ibasiccut=basiccut;
559    
560     //Step 2: Refine the cut
561     TCanvas *ckin = new TCanvas("ckin","Kinematic Plots (in the making)",600,600);
562     ckin->SetLogy(logscale);
563    
564     TH1F *datahistoa = allsamples.Draw("datahistoa",variable,nbins,min,max, xlabel, "events",observed,is_data,luminosity);
565     TH1F *datahistob = allsamples.Draw("datahistob",variable,nbins,min,max, xlabel, "events",predicted,is_data,luminosity);
566    
567     datahistoa->SetLineColor(kBlue);
568     datahistob->SetLineColor(kRed);
569    
570     TLegend *kinleg = new TLegend(0.6,0.7,0.8,0.89);
571     kinleg->SetFillColor(kWhite);
572     kinleg->SetLineColor(kWhite);
573     kinleg->SetBorderSize(0);
574     kinleg->AddEntry(datahistoa,"Observed "+TString(is_data?"Data":"MC"),"l");
575     kinleg->AddEntry(datahistob,"Predicted "+TString(is_data?"Data":"MC"),"l");
576    
577     datahistoa->SetMaximum(ymax);
578     datahistoa->Draw("histo");
579     datahistob->SetLineStyle(2);
580     if ( !noscale ) datahistob->Scale(0.3);
581     datahistob->Draw("histo,sames");
582     TVirtualPad::Pad()->Update();
583    
584     TPaveStats *sa = (TPaveStats*)datahistoa->GetListOfFunctions()->FindObject("stats");
585     TPaveStats *sb = (TPaveStats*)datahistob->GetListOfFunctions()->FindObject("stats");
586     if ( sa && sb ) {
587     sa->SetTextColor(datahistoa->GetLineColor());
588     sb->SetTextColor(datahistob->GetLineColor());
589     sb->SetY1NDC(sb->GetY1NDC()-0.25);
590     sb->SetY2NDC(sb->GetY2NDC()-0.25);
591     TVirtualPad::Pad()->Update();
592     }
593     kinleg->Draw();
594     TText* write_cut = write_title(variable);
595     write_cut->Draw();
596     CompleteSave(ckin,"special_kin/"+filename);
597     datahistoa->Delete();
598     datahistob->Delete();
599     delete ckin;
600     }
601    
602     TH1F *give_lm0_pred(TCut bcut,string mcjzb) {
603     int nbins=50;
604     TH1F *jzbn = signalsamples.Draw("jzbn","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events", cutOSSF&&cutnJets&&cutmass,mc,luminosity,signalsamples.FindSample("LM0"));
605     TH1F *jzbno = signalsamples.Draw("jzbno","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSOF&&cutnJets&&cutmass,mc,luminosity,signalsamples.FindSample("LM0"));
606     TH1F *jzbpo = signalsamples.Draw("jzbp",mcjzb,nbins,0,350, "JZB [GeV]", "events", cutOSOF&&cutnJets&&cutmass,mc,luminosity,signalsamples.FindSample("LM0"));
607    
608     //Sidebands
609     flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
610     TH1F *jzbnos;
611     TH1F *jzbpos;
612     TH1F *jzbnss;
613     TH1F *jzbpss;
614    
615     if(PlottingSetup::RestrictToMassPeak) {
616     jzbnos = signalsamples.Draw("jzbnos","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSOF&&cutnJets&&sidebandcut,mc,luminosity,signalsamples.FindSample("LM0"));
617     jzbpos = signalsamples.Draw("jzbpos",mcjzb,nbins,0,350, "JZB [GeV]", "events", cutOSOF&&cutnJets&&sidebandcut,mc,luminosity,signalsamples.FindSample("LM0"));
618     jzbnss = signalsamples.Draw("jzbnss","-"+mcjzb,nbins,0,350, "JZB [GeV]", "events",cutOSSF&&cutnJets&&sidebandcut,mc,luminosity,signalsamples.FindSample("LM0"));
619     jzbpss = signalsamples.Draw("jzbpss",mcjzb,nbins,0,350, "JZB [GeV]", "events", cutOSSF&&cutnJets&&sidebandcut,mc,luminosity,signalsamples.FindSample("LM0"));
620     }
621    
622     TH1F *pred = (TH1F*)jzbn->Clone("pred");
623     flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
624     if(PlottingSetup::RestrictToMassPeak) {
625     pred->Add(jzbno,-1.0/3);
626     pred->Add(jzbpo,1.0/3);
627     pred->Add(jzbnos,-1.0/3);
628     pred->Add(jzbpos,1.0/3);
629     pred->Add(jzbnss,-1.0/3);
630     pred->Add(jzbpss,1.0/3);
631     } else {
632     pred->Add(jzbno,-1.0);
633     pred->Add(jzbpo,1.0);
634     }
635    
636     pred->SetLineColor(kRed);
637     pred->SetMinimum(0);
638    
639     delete jzbn;
640     delete jzbpo;
641     delete jzbno;
642     delete jzbpos;
643     delete jzbnos;
644     delete jzbpss;
645     delete jzbnss;
646     cout << "pred contains " << pred->Integral() << "entries" << endl;
647     return pred;
648     }
649    
650    
651     void lm0_illustration() {
652     TCanvas *can = new TCanvas("can","Signal Background Comparison Canvas");
653     can->SetLogy(1);
654    
655     int sbg_nbins=130;
656     float sbg_min=-500; //-110;
657     float sbg_max=800; //jzbHigh;
658    
659     float simulatedlumi=1000;//in pb please
660    
661     TH1F *JZBplotZJETs = allsamples.Draw("JZBplotZJETs",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,simulatedlumi,allsamples.FindSample("DYJetsToLL"));
662     TH1F *JZBplotLM0 = signalsamples.Draw("JZBplotLM0",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,simulatedlumi,signalsamples.FindSample("LM0"));
663     TH1F *JZBplotTtbar = allsamples.Draw("JZBplotTtbar",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,simulatedlumi,allsamples.FindSample("TTJets"));
664    
665     JZBplotTtbar->SetLineColor(allsamples.GetColor("TTJet"));
666     JZBplotZJETs->SetFillColor(allsamples.GetColor("DY"));
667     JZBplotZJETs->SetLineColor(kBlack);
668     JZBplotLM0->SetLineStyle(2);
669     JZBplotZJETs->SetMaximum(JZBplotZJETs->GetMaximum()*5);
670     JZBplotZJETs->SetMinimum(1);
671    
672     JZBplotTtbar->SetMaximum(JZBplotZJETs->GetMaximum());
673     JZBplotTtbar->SetMinimum(0.01);
674     JZBplotTtbar->SetFillColor(allsamples.GetColor("TTJets"));
675     JZBplotTtbar->DrawClone("histo");
676     JZBplotZJETs->Draw("histo,same");
677     JZBplotTtbar->SetFillColor(0);
678     JZBplotTtbar->DrawClone("histo,same");
679     JZBplotTtbar->SetFillColor(allsamples.GetColor("TTJets"));
680     JZBplotLM0->Draw("histo,same");
681    
682    
683     TLegend *signal_bg_comparison_leg2 = make_legend("",0.55,0.75,false);
684     signal_bg_comparison_leg2->AddEntry(JZBplotZJETs,"Z+Jets","f");
685     signal_bg_comparison_leg2->AddEntry(JZBplotTtbar,"t#bar{t}","f");
686     signal_bg_comparison_leg2->AddEntry(JZBplotLM0,"LM0","f");
687     signal_bg_comparison_leg2->Draw();
688     TText* title = write_title("CMS MC simulation, #sqrt{s}= 7 TeV @ L=1 fb^{-1}");
689     title->Draw();
690     CompleteSave(can,"LM0/jzb_bg_vs_signal_distribution__LM0");
691     can->SetLogy(0);
692     TH1F *lm0jzbp = signalsamples.Draw("lm0jzb",jzbvariablemc,70,0,350, "JZB [GeV]", "events",cutOSSF&&cutnJets&&cutmass,mc,luminosity,signalsamples.FindSample("LM0"));
693     TH1F *lm0pred = (TH1F*) give_lm0_pred("mll>0",jzbvariablemc);
694     lm0pred->Draw("histo");
695     lm0jzbp->Draw("histo,same");
696     TLegend *leg = make_legend("LM0");
697     leg->SetFillColor(kWhite);
698     leg->SetLineColor(kWhite);
699     leg->SetHeader("LM0 (MC)");
700     leg->AddEntry(lm0jzbp,"Observed","l");
701     leg->AddEntry(lm0pred,"Predicted","l");
702     leg->Draw("same");
703     CompleteSave(can,"LM0/LeftRight");
704    
705     delete lm0jzbp;
706     delete lm0pred;
707     }
708    
709     void kinematic_dist_of_pred_and_obs() {//former plot_list
710     gStyle->SetOptStat("oueMri");
711     TCut observed=(cutmass&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)");
712    
713     TCut predicted = (cutmass&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)")
714     || (sidebandcut&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)")
715     || (sidebandcut&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)");
716     TCut predictedMC = (cutmass&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.04*pt-1.82559)>100)")
717     || (sidebandcut&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.04*pt-1.82559)>100)")
718     || (sidebandcut&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.04*pt-1.82559)>100)");
719    
720     flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
721     if(!PlottingSetup::RestrictToMassPeak) {
722     predicted = TCut((cutmass&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)"));
723     predictedMC = TCut((cutmass&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.04*pt-1.82559)>100)"));
724     }
725     // TCut predicted=((cutmass&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)<-100)")
726     // ||(cutmass&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)")
727     // ||(cutmass&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)<-100)")
728     // ||(sidebandcut&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)")
729     // ||(sidebandcut&&cutOSOF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)<-100)")
730     // ||(sidebandcut&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)>100)")
731     // ||(sidebandcut&&cutOSSF&&cutnJets&&basiccut&&"((jzb[1]+0.06*pt-2.84727)<-100)"));
732    
733     bool doPF=false;
734     bool dolog=true;
735     bool nolog=false;
736    
737     // Mll: do not scale
738     make_double_plot("mll",20,50,150,11,nolog,"m_{ll} [GeV]","mll",observed,predicted,data,true);
739     make_double_plot("mll",20,50,150,11,nolog,"m_{ll} [GeV]","mll_MC",observed,predictedMC,mc,true);
740     make_double_plot("met[4]",20,0,400,11,nolog,"pfMET [GeV]","pfMET",observed,predicted,data);
741     make_double_plot("met[4]",20,0,400,11,nolog,"pfMET [GeV]","pfMET_MC",observed,predictedMC,mc);
742     make_double_plot("jetpt[0]",10,0,400,11,nolog,"leading jet p_{T} [GeV]","pfJetGoodPt_0",observed,predicted,data);
743     make_double_plot("jetpt[0]",10,0,400,11,nolog,"leading jet p_{T} [GeV]","pfJetGoodPt_0_MC",observed,predictedMC,mc);
744     make_double_plot("jeteta[0]",10,-5,5,11,nolog,"leading jet #eta","pfJetGoodEta_0",observed,predicted,data);
745     make_double_plot("jeteta[0]",10,-5,5,11,nolog,"leading jet #eta","pfJetGoodEta_0_MC",observed,predictedMC,mc);
746     make_double_plot("pt",10,0,300,11,nolog,"Z p_{T} [GeV]","Zpt",observed,predicted,data);
747     make_double_plot("pt",10,0,300,11,nolog,"Z p_{T} [GeV]","Zpt_MC",observed,predictedMC,mc);
748     make_double_plot("pt1",10,0,200,15,nolog,"p_{T} [GeV]","pt1",observed,predicted,data);
749     make_double_plot("pt1",10,0,200,15,nolog,"p_{T} [GeV]","pt1_MC",observed,predictedMC,mc);
750     make_double_plot("pt2",10,0,200,25,nolog,"p_{T} [GeV]","pt2",observed,predicted,data);
751     make_double_plot("pt2",10,0,200,25,nolog,"p_{T} [GeV]","pt2_MC",observed,predictedMC,mc);
752     make_double_plot("eta1",10,-5,5,11,nolog,"#eta_{1,l}","eta_1",observed,predicted,data);
753     make_double_plot("eta1",10,-5,5,11,nolog,"#eta_{1,l}","eta_1_MC",observed,predictedMC,mc);
754     make_double_plot("eta2",10,-5,5,11,nolog,"#eta_{2,l}","eta_2",observed,predicted,data);
755     make_double_plot("eta2",10,-5,5,11,nolog,"#eta_{2,l}","eta_2_MC",observed,predictedMC,mc);
756     make_double_plot("phi1-phi2",10,-6.0,6.0,11,nolog,"#phi_{1}-#phi_{2}","dphi",observed,predicted,data);
757     make_double_plot("phi1-phi2",10,-6.0,6.0,11,nolog,"#phi_{1}-#phi_{2}","dphi_MC",observed,predictedMC,mc);
758     make_double_plot("pfJetGoodNum",8,0.5,8.5,20,nolog,"nJets","nJets",observed,predicted,data);
759     make_double_plot("pfJetGoodNum",8,0.5,8.5,20,nolog,"nJets","nJets_MC",observed,predictedMC,mc);
760    
761     }
762 buchmann 1.2
763    
764 buchmann 1.4 void generator_study_plots() {
765    
766     //uncomment whichever one you want to see :-)
767    
768 buchmann 1.5
769     GenLevelStudy::X_vs_generation_lm4();
770 buchmann 1.4 /*
771     GenLevelStudy::compare_sms_lm4();
772     GenLevelStudy::MomentumFraction();
773     GenLevelStudy::AngleMETsumLSP();
774     GenLevelStudy::RatioMETsumLSP();
775     GenLevelStudy::AngleLSPLSP();
776     GenLevelStudy::AngleLSPZ();
777     GenLevelStudy::WidthIllustration();
778     GenLevelStudy::ZDecayIllustration();
779     GenLevelStudy::DeltaLSPmomentum();
780     GenLevelStudy::pStarIllustration(0.5);
781     GenLevelStudy::pStarIllustration(0.25);
782     GenLevelStudy::pStarIllustration(0.75);
783 buchmann 1.5 GenLevelStudy::pStarDistributions();
784 buchmann 1.4 GenLevelStudy::DrawJetBand();
785 buchmann 1.5 GenLevelStudy::DrawOnly100to150inJetBand();
786     GenLevelStudy::ImpactOfGluinoChi2MassDifference();
787     */
788 buchmann 1.4 }
789    
790    
791 buchmann 1.2 void jzb_negative_generator_study() {
792     write_warning(__FUNCTION__,"We are going to use a t5zz scan file, and \033[1;31m WON'T \033[1;35m cut on MassGlu/MassLSP in order to improve statistics. This is ok for small studies, for a real study you'll need to look at points individually ...");
793    
794    
795     scansample.AddSample("/scratch/buchmann/ntuples/GeneratorInformationInJZB___JZBplusSamples_TestingSMS_v7/SMS-T5zz_x-05_Mgluino-150to1200_mLSP-50to1150_7TeV-Pythia6Z__Summer11-PU_START42_V11_FastSim-v2__AODSIM___inclindex_v2.root","SMST5zz",1,1,false,false,0,kRed);
796     TCanvas *jcan = new TCanvas("jcan","jcan");
797     scansample.collection[scansample.collection.size()-1].events->Draw("(LSP1pt/LSP1Mopt):pureGeneratorJZB","genNjets>2","PROF");
798     TH1F *h = new TH1F("h","h",100,-500,500);
799     h->SetLineColor(kBlack);
800     TProfile *p = (TProfile*)jcan->GetPrimitive("htemp");
801     p->GetXaxis()->SetTitle("Generator JZB");
802     p->GetXaxis()->CenterTitle();
803     p->GetYaxis()->SetTitle("( LSP p_{T} ) / ( LSP mother p_{T} )");
804     p->GetYaxis()->CenterTitle();
805     p->SetLineColor(kBlue);
806 buchmann 1.3 TLegend* leg = make_legend("", 0.40, 0.75, false);
807 buchmann 1.2 leg->AddEntry(p,"LSP1 pt / LSP1 Mother pt","l");
808     leg->AddEntry(h,"Z pt / Z Mother pt","l");
809     leg->Draw();
810     TText *title = write_title("JZB as a function of the first LSP's momentum transfer");
811     title->Draw();
812     scansample.collection[scansample.collection.size()-1].events->Draw("(genZPt/LSP1Mopt):pureGeneratorJZB","genNjets>2","PROF,same");
813    
814     CompleteSave(jcan,"NegativeJZBStudy/LSPpt_LSPMopt");
815    
816     scansample.collection[scansample.collection.size()-1].events->Draw("angleLSPLSP:pureGeneratorJZB","genNjets>2","PROF");
817     TProfile *p1 = (TProfile*)jcan->GetPrimitive("htemp");
818     p1->GetXaxis()->SetTitle("Generator JZB");
819     p1->GetXaxis()->CenterTitle();
820     p1->GetYaxis()->SetTitle("#angle(LSP1,LSP2)");
821     p1->GetYaxis()->CenterTitle();
822     TText *title1 = write_title("JZB as a function of the angle between the two LSPs");
823     title1->Draw();
824     CompleteSave(jcan,"NegativeJZBStudy/AngleLSPLSP");
825    
826     TH1F *jzbdistributionvsz[5];
827     THStack zstack;
828     jcan->SetLogy(1);
829 buchmann 1.3 TLegend* leg2 = make_legend("", 0.15, 0.75, false);
830     leg2->SetX2(0.4);
831     TLegend* leg3 = make_legend("", 0.15, 0.75, false);
832     leg3->SetX2(0.4);
833 buchmann 1.2 for(int z=0;z<5;z++) {
834     stringstream specialcut;
835     if(z<4) specialcut << "genNjets>2&&(LSP1pt/LSP1Mopt)>" << z*0.2 << "&&(LSP1pt/LSP1Mopt)<" << (z+1)*0.2;
836     else specialcut << "genNjets>2&&(LSP1pt/LSP1Mopt)>" << z*0.2;
837     stringstream histtitle;
838     histtitle << "splitup_" << z;
839     stringstream ntitle;
840     if(z<4) ntitle << z*0.2 << " < z < " << (z+1)*0.2;
841     else ntitle << z*0.2 << " < z";
842     jzbdistributionvsz[z] = scansample.Draw(histtitle.str().c_str(),"pureGeneratorJZB",100,-500,500, "generator JZB (GeV)", "events",specialcut.str().c_str(),mc,1.0);
843     jzbdistributionvsz[z]->SetLineColor(z+1);
844     jzbdistributionvsz[z]->SetMarkerSize(0);
845     zstack.Add(jzbdistributionvsz[z]);
846 buchmann 1.3 leg2->AddEntry(jzbdistributionvsz[z],ntitle.str().c_str(),"f");
847     leg3->AddEntry(jzbdistributionvsz[z],ntitle.str().c_str(),"l");
848 buchmann 1.2 }
849    
850     jzbdistributionvsz[0]->GetYaxis()->SetRangeUser(2,800);
851     jzbdistributionvsz[0]->DrawNormalized("histo");
852     for(int z=0;z<5;z++) {
853     jzbdistributionvsz[z]->DrawNormalized("histo,same");
854     }
855    
856     // zstack.Draw("nostack,histo");
857 buchmann 1.3 leg3->Draw("same");
858 buchmann 1.2 CompleteSave(jcan,"NegativeJZBStudy/StackedAccordingToMomentumFractionIndividual");
859    
860     for(int z=0;z<5;z++) {
861     jzbdistributionvsz[z]->SetFillColor(z+1);
862     }
863    
864     zstack.Draw("histo");
865 buchmann 1.3 leg2->Draw("same");
866 buchmann 1.2 CompleteSave(jcan,"NegativeJZBStudy/StackedAccordingToMomentumFractionStacked");
867    
868     //varangle vasysyn
869    
870     // scansample.collection[scansample.collection.size()-1].events->Draw("(LSPPromptnessLevel[0]/4.0)*angleLSPLSP/(LSP1pt/LSP1Mopt):pureGeneratorJZB","genNjets>2","PROF");
871    
872    
873    
874     }
875    
876     string ReplaceCharacter(string originalstring,string replacethis,string replacewiththis)
877     {
878     int pos = originalstring.find(replacethis);
879     if(pos == -1) return originalstring;
880     originalstring.replace(pos, replacewiththis.length(), replacewiththis);
881     return originalstring;
882     }
883     string removefunnystring(string name) {
884     name=ReplaceCharacter(name,"[","_");
885     name=ReplaceCharacter(name,"]","_");
886     name=ReplaceCharacter(name,"{","_");
887     name=ReplaceCharacter(name,"}","_");
888     name=ReplaceCharacter(name,".","_");
889     name=ReplaceCharacter(name,",","_");
890     name=ReplaceCharacter(name,";","_");
891     name=ReplaceCharacter(name,":","_");
892     name=ReplaceCharacter(name,"'","_");
893     name=ReplaceCharacter(name,"$","_");
894     name=ReplaceCharacter(name,"@","_");
895     return name;
896     }
897    
898     void compare_lm4_sms_variable(TTree *eventsLM4, TTree *eventsSMS, string variable, int nbins, float xmin, float xmax, TCut cut, string saveas, bool dology=false) {
899     TCanvas *can = new TCanvas("can","can");
900     can->SetLogy(dology);
901     TH1F *hlm4 = new TH1F("hlm4","hlm4",nbins,xmin,xmax);
902     TH1F *hsms = new TH1F("hsms","hsms",nbins,xmin,xmax);
903     eventsLM4->Draw((variable+">>hlm4").c_str(),cut,"goff");
904     eventsSMS->Draw((variable+">>hsms").c_str(),cut,"goff");
905     hlm4->SetLineColor(kBlue);
906     hsms->SetLineColor(kRed);
907    
908     if(hlm4->Integral()>0) hlm4->Scale(1.0/hlm4->Integral());
909     else write_warning(__FUNCTION__,"Watch out, LM4 histo is empty!");
910     if(hsms->Integral()>0) hsms->Scale(1.0/hsms->Integral());
911     else write_warning(__FUNCTION__,"Watch out, SMS histo is empty!");
912    
913     float min=get_nonzero_minimum(hlm4);
914     float max=hlm4->GetMaximum();
915     if(get_nonzero_minimum(hsms)<min) min=get_nonzero_minimum(hsms);
916     if(hsms->GetMaximum()>max) max=hsms->GetMaximum();
917     if(dology) max*=5;
918     else max*=2;
919    
920     hlm4->GetYaxis()->SetRangeUser(min,max);
921     hlm4->GetXaxis()->SetTitle(variable.c_str());
922     hlm4->GetXaxis()->CenterTitle();
923     hlm4->Draw("histo");
924     hsms->Draw("histo,same");
925     TLegend *leg = make_legend("",0.2,0.98,false);
926     leg->SetY2(1.0);
927     leg->SetNColumns(2);
928     leg->AddEntry(hlm4,"LM4","l");
929     leg->AddEntry(hsms,"\"LM4\" SMS","l");
930     leg->Draw();
931     stringstream saveas2;
932     saveas2 << "ComparingLM4_SMS/" << removefunnystring(saveas);
933     CompleteSave(can,saveas2.str());
934     delete can;
935     delete hlm4;
936     delete hsms;
937     }
938    
939    
940     void compare_LM4_and_SMS() {
941     TFile *f1 = new TFile("/shome/lbaeni/jzb/LM4_SMS/SMS_LM4_JZB.root");
942     TTree *LM4events = (TTree*)f1->Get("events");
943 buchmann 1.4 TFile *f2 = new TFile("/scratch/buchmann/ntuples/GeneratorInformationInJZB___JZBplusSamples_TestingSMS_v5/LM4_SUSY_sftsht_7TeV-pythia6__Summer11-PU_S4_START42_V11-v2__withIndex.root");
944 buchmann 1.2 TTree *SMSevents = (TTree*)f2->Get("events");
945    
946     compare_lm4_sms_variable(LM4events, SMSevents, "mll",100,50,150,cutOSSF&&cutnJets&&basiccut,"mll",true);
947     compare_lm4_sms_variable(LM4events, SMSevents, "jzb[1]+0.04*pt-2.32212",100,-300,700,cutmass&&cutOSSF&&cutnJets&&basiccut,"jzb",true);
948 buchmann 1.3 compare_lm4_sms_variable(LM4events, SMSevents, "pureGeneratorJZB",100,-300.0,700.0,cutOSSF&&basiccut,"pureGeneratorJZB",true);
949 buchmann 1.2 compare_lm4_sms_variable(LM4events, SMSevents, "pfJetGoodNum",10,-0.5,9.5,cutmass&&cutOSSF&&basiccut,"pfJetGoodNum",true);
950     compare_lm4_sms_variable(LM4events, SMSevents, "pt",100,15.0,200.0,cutOSSF&&basiccut,"pt",false);
951     compare_lm4_sms_variable(LM4events, SMSevents, "pt1",100,15.0,100.0,cutOSSF&&basiccut,"pt1",false);
952     compare_lm4_sms_variable(LM4events, SMSevents, "pt2",100,15.0,100.0,cutOSSF&&basiccut,"pt2",false);
953     compare_lm4_sms_variable(LM4events, SMSevents, "met[4]",100,0.0,600.0,cutOSSF&&basiccut,"met",false);
954 buchmann 1.3 compare_lm4_sms_variable(LM4events, SMSevents, "genMET",100,0.0,600.0,cutOSSF&&basiccut,"genMET",false);
955     compare_lm4_sms_variable(LM4events, SMSevents, "genNjets",10,-0.5,9.5,cutOSSF&&basiccut,"genNjets",true);
956 buchmann 1.2 }
957 buchmann 1.5
958    
959     void compare_onpeak_offpeak_signal_distributions(string mcjzb, string datajzb, int sampleindex) {
960     cout << (signalsamples.collection)[sampleindex].filename << endl;
961     }
962    
963     void compare_onpeak_offpeak_signal_distributions(string mcjzb, string datajzb) {
964     for(int i=0;i<signalsamples.collection.size();i++) compare_onpeak_offpeak_signal_distributions(mcjzb,datajzb,i);
965     }
966