ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/ZbTools.C
Revision: 1.15
Committed: Tue Nov 27 07:25:55 2012 UTC (12 years, 5 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.14: +135 -52 lines
Log Message: 
Added new classes to store Z+b results for different scenarios (reference, inclusive, efficiency/mistag uncertainty up&down, medium WP, alpha up and down

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 <TLegendEntry.h>
18    
19     using namespace std;
20    
21     using namespace PlottingSetup;
22    
23     namespace ZbData {
24     vector<float> data_over_mc;
25    
26 buchmann 1.10 TCut ZplusBsel("pt1>20&&pt2>20&&mll>60&&mll<120&&id1==id2");
27 buchmann 1.9 TCut LeadingB("Zb3010_bTagProbCSVBP[0]>0.898");
28     TCut EtaB("abs(Zb3010_pfJetGoodEta[0])<1.3");
29 buchmann 1.10 TCut PhiZcut("abs(Zb3010_pfJetDphiZ[0])>2.7");
30 buchmann 1.1
31     }
32    
33 buchmann 1.15 class ZbPointResult {
34     public:
35     ZbPointResult(string var, float ptlow, float pthigh, float alphalow, float alphahigh, Value Data, Value MC);
36     float PtLow;
37     float PtHigh;
38     float AlphaLow;
39     float AlphaHigh;
40     string variable;
41     Value Data;
42     Value MC;
43     Value Ratio;
44     };
45    
46     ZbPointResult::ZbPointResult(string var, float ptlow, float pthigh, float alphalow, float alphahigh, Value data, Value mc) :
47     variable(var), PtLow(ptlow), PtHigh(pthigh), AlphaLow(alphalow), AlphaHigh(alphahigh), Data(data), MC(mc) {
48     Ratio = data/mc;
49     }
50    
51     std::ostream &operator<<(std::ostream &ostr, ZbPointResult &b)
52     {
53     ostr << b.PtLow << ";" << b.PtHigh << ";" << b.AlphaLow << ";" << b.AlphaHigh << ";" << b.Data.getValue() << ";" << b.Data.getError() << ";" << b.MC.getValue() << ";" << b.MC.getError() << ";" << b.Ratio.getValue() << ";" << b.Ratio.getError();
54     return ostr;
55     }
56    
57    
58     class ZbResultContainer {
59     public:
60     string name;
61     vector<ZbPointResult> MPF_Rabs_location;
62     Value MPF_Result_FaceValue;
63     Value Rabs_Result_FaceValue;
64     Value MPF_Result_Extrapolated;
65     Value Rabs_Result_Extrapolated;
66    
67     ZbResultContainer(string name);
68     ZbResultContainer(const ZbResultContainer &old);
69     };
70    
71    
72     ZbResultContainer::ZbResultContainer(string _name) {
73     name=_name;
74     cout<< "Result Container for \"" << name << "\" has been initialized" << endl;
75    
76     }
77    
78     std::ostream &operator<<(std::ostream &ostr, ZbResultContainer &b)
79     {
80     ostr << "Results for " << b.name << endl;
81     ostr << " MPF and Rabs locations: " << endl;
82     for(int i=0;i<b.MPF_Rabs_location.size();i++) ostr << " " << b.MPF_Rabs_location[i] << endl;
83     ostr << " Results: " << endl;
84     ostr << " Face Value: " << endl;
85     ostr << " MPF " << b.MPF_Result_FaceValue << endl;
86     ostr << " Rabs " << b.Rabs_Result_FaceValue << endl;
87     ostr << " Extrapolated: " << endl;
88     ostr << " MPF " << b.MPF_Result_Extrapolated << endl;
89     ostr << " Rabs " << b.Rabs_Result_Extrapolated << endl;
90    
91     return ostr;
92     }
93    
94    
95    
96 buchmann 1.1 using namespace ZbData;
97    
98 buchmann 1.7
99     //*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/* DELETE EVERYTHING BETWEEN THIS LINE /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
100     void SpecialCutFlow(string identifier);
101    
102     //*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/* AND THIS LINE /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
103    
104    
105 buchmann 1.1 void print_yield(TCut cut) {
106     TH1F *data = allsamples.Draw("data", "mll",1,0,10000000, "m_{ll} [GeV]", "events", cut,data,luminosity);
107 buchmann 1.7 dout << " Yields for " << (const char*) cut << endl;
108     dout << " data: " << data->Integral() << endl;
109 buchmann 1.9 THStack mcm = allsamples.DrawStack("mc", "mll",1,0,10000000, "m_{ll} [GeV]", "events", cut,mc,luminosity);
110 buchmann 1.1 int nhists = mcm.GetHists()->GetEntries();
111 buchmann 1.7 dout << "Going to loop over " << nhists << " histograms " << endl;
112 buchmann 1.1 TFile *test = new TFile("test.root","RECREATE");
113     mcm.Write();
114     test->Close();
115     for (int i = 0; i < nhists; i++) {
116     TH1* hist = static_cast<TH1*>(mcm.GetHists()->At(i));
117     cout << " " << hist->GetName() << ": " << hist->Integral() << endl;
118     }
119 buchmann 1.9
120 buchmann 1.1 cout << endl << endl;
121 buchmann 1.9
122 buchmann 1.1 delete data;
123     }
124    
125 buchmann 1.2 void draw_kin_variable(string variable, TCut cut, int nbins, float min, float max, string xlabel, string saveas, bool dologscale, string speciallabel="") {
126 buchmann 1.1
127     TCanvas *ckin = new TCanvas("ckin","kin variable canvas");
128     ckin->SetLogy(dologscale);
129     TH1F *datahisto = allsamples.Draw("datahisto",variable,nbins,min,max, xlabel, "events",cut,data,luminosity);
130     datahisto->SetMarkerSize(DataMarkerSize);
131     THStack mcstack = allsamples.DrawStack("mcstack",variable,nbins,min,max, xlabel, "events",cut,mc,luminosity);
132 buchmann 1.2 if(dologscale) datahisto->SetMaximum(3*datahisto->GetMaximum());
133     else datahisto->SetMaximum(1.5*datahisto->GetMaximum());
134     TText *speciall;
135     if(speciallabel!="") {
136     speciall = new TText(0.2,0.8,speciallabel.c_str());
137     speciall->Draw();
138     }
139 buchmann 1.1 datahisto->Draw("e1");
140     ckin->Update();
141     mcstack.Draw("same");
142     datahisto->Draw("same,e1");
143     TLegend *kinleg = allsamples.allbglegend();
144     kinleg->Draw();
145    
146     TPad *kinpad = new TPad("kinpad","kinpad",0,0,1,1);
147     kinpad->cd();
148     kinpad->SetLogy(dologscale);
149     datahisto->Draw("e1");
150     mcstack.Draw("same");
151     datahisto->Draw("same,e1");
152     datahisto->Draw("same,axis");
153     kinleg->Draw();
154     DrawPrelim();
155     saveas="kin/"+saveas;
156     save_with_ratio(datahisto,mcstack,kinpad->cd(),saveas);
157    
158    
159    
160    
161     ZbData::data_over_mc.push_back(datahisto->Integral()/CollapseStack(mcstack)->Integral());
162    
163     datahisto->Delete();
164     delete ckin;
165     }
166    
167     void print_all_b_yields() {
168     cout << "Basic selection with a b jet" << endl;
169 buchmann 1.9 print_yield(ZplusBsel&&"Zb3010_pfJetGoodNumBtag>0");
170 buchmann 1.2 cout << "Leading jet is a b-jet " << endl;
171 buchmann 1.1 print_yield(ZplusBsel&&LeadingB);
172     cout << "|#eta_{b}|<1.3 " << endl;
173     print_yield(ZplusBsel&&LeadingB&&EtaB);
174     cout << "|#delta#phi(b<Z)|>2.7" << endl;
175     print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
176 buchmann 1.2 cout << "10<ptZ<1000 GeV" << endl;
177 buchmann 1.1 print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000");
178 buchmann 1.13 cout << "#alpha < 0.35" << endl;
179     print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.35");
180 buchmann 1.1 cout << "#alpha < 0.3" << endl;
181 buchmann 1.9 print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.3");
182 buchmann 1.1 cout << "#alpha < 0.2" << endl;
183 buchmann 1.9 print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.2");
184 buchmann 1.1 cout << "#alpha < 0.15" << endl;
185 buchmann 1.9 print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.15");
186 buchmann 1.1 cout << "#alpha < 0.1" << endl;
187 buchmann 1.9 print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.1");
188     // cout << "#alpha < 0.05" << endl;
189     // print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.05");
190 buchmann 1.1 }
191    
192 buchmann 1.15 void TEST_DrawEvilCutFlow() {
193 buchmann 1.7 stringstream MegaCut;
194    
195 buchmann 1.9 MegaCut << " 1*(" << (const char*) (ZplusBsel&&TCut("Zb3010_pfJetGoodNumBtag>0")) << ")";
196 buchmann 1.7 MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB) << ")";
197     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB) << ")";
198     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut) << ")";
199     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000")) << ")";
200 buchmann 1.9 MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.3")) << ")";
201     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.2")) << ")";
202     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.15")) << ")";
203     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.1")) << ")";
204     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.05")) << ")";
205 buchmann 1.7 MegaCut << " ";
206    
207     TCanvas *can = new TCanvas("can","can");
208     can->SetLogy(1);
209     TCut basecut("mll>6||mll<10");
210    
211    
212    
213     TLegend *leg = allsamples.allbglegend();
214    
215    
216     TH1F *data = allsamples.Draw("data", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,data,luminosity);
217     THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
218    
219     float runningsum=0;
220     for(int i=data->GetNbinsX();i>0;i--) {
221     runningsum+=data->GetBinContent(i);
222     data->SetBinContent(i,runningsum);
223     }
224    
225     float minimum=data->GetMaximum();
226    
227     TH1F* h;
228     TIter nextOF(mcm.GetHists());
229     while ( h = (TH1F*)nextOF() ) {
230     float runningsum=0;
231     minimum=data->GetMaximum();//done deliberately at each step so get the minimum of the last (!) contribution
232     for(int i=h->GetNbinsX();i>0;i--) {
233     runningsum+=h->GetBinContent(i);
234     h->SetBinContent(i,runningsum);
235     if(runningsum<minimum&&runningsum>0) minimum=runningsum;
236     }
237     }
238    
239    
240    
241     data->SetMinimum(0.05*minimum);
242     can->cd(1)->SetBottomMargin(0.25);
243     data->GetXaxis()->SetBinLabel(1,"Pre-selection");
244     data->GetXaxis()->SetBinLabel(2,"Contains b jet");
245     data->GetXaxis()->SetBinLabel(3,"Leading jet is b-jet");
246     data->GetXaxis()->SetBinLabel(4,"|#eta_{b}|<1.3 ");
247     data->GetXaxis()->SetBinLabel(5,"|#delta#phi(b,Z)|>2.7");
248     data->GetXaxis()->SetBinLabel(6,"10<p_{T}^{Z}<1000 GeV");
249     data->GetXaxis()->SetBinLabel(7,"#alpha < 0.3");
250     data->GetXaxis()->SetBinLabel(8,"#alpha < 0.2");
251     data->GetXaxis()->SetBinLabel(9,"#alpha < 0.15");
252     data->GetXaxis()->SetBinLabel(10,"#alpha < 0.1");
253     data->GetXaxis()->SetBinLabel(11,"#alpha < 0.05");
254     data->GetXaxis()->LabelsOption("v");
255    
256     data->Draw("e1");
257     mcm.Draw("same");
258     data->Draw("same,e1,axis");
259     data->Draw("same,e1");
260     // data->Draw("same,TEXT");
261    
262     leg->Draw();
263    
264     DrawPrelim();
265    
266     CompleteSave(can,"CutFlow");
267     }
268    
269 buchmann 1.1 void draw_Zb_kin_vars() {
270 buchmann 1.6 draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",25,0,200,"Z p_{T}","Official/Zpt",1);
271 buchmann 1.9 draw_kin_variable("Zb3010_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt",1);
272     draw_kin_variable("Zb3010_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt",1);
273     draw_kin_variable("Zb3010_pfJetGoodPt[1]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,2,"p_{T}^{J2} / p_{T}^{Z}","Official/SubLeadingJetPt_Over_ZPt",1);
274     draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<50", 40,0,2,"#alpha","Official/alpha_pt_30_to_50",1);
275     draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>50&&pt<75", 40,0,2,"#alpha","Official/alpha_50_to_75",1);
276     draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>75&&pt<125", 40,0,2,"#alpha","Official/alpha_pt_75_to_125",1);
277     draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>125&&pt<1000",40,0,2,"#alpha","Official/alpha_pt_125_to_1000",1);
278     draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",40,0,2,"#alpha","Official/alpha_pt_30_to_1000",1);
279 buchmann 1.6
280 buchmann 1.10 draw_kin_variable("Zb3010_pfBJetDphiZ[0]",ZplusBsel&&LeadingB&&"pt>10&&pt<1000&&pfBJetDphiZ[0]>0",30,0,3.2,"#delta#phi (Z,b lead)","DeltaPhiZBlead",0);
281 buchmann 1.6
282     /*
283     draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,200,"p_{T}^{l1}","Lep/pt1",1);
284     draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id1==0",20,0,200,"p_{T}^{e1}","Lep/pt1_e",1);
285     draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id1==1",20,0,200,"p_{T}^{#mu1}","Lep/pt1_m",1);
286     draw_kin_variable("pt2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,200,"p_{T}^{l2}","Lep/pt2",1);
287     draw_kin_variable("pt2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id2==0",20,0,200,"p_{T}^{e2}","Lep/pt2_e",1);
288     draw_kin_variable("pt2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id2==1",20,0,200,"p_{T}^{#mu2}","Lep/pt2_m",1);
289     draw_kin_variable("eta1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,-3.1,3.1,"#eta^{l1}","Lep/eta1",1);
290     draw_kin_variable("eta1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id1==0",20,-3.1,3.1,"#eta^{e1}","Lep/eta1_e",1);
291     draw_kin_variable("eta1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id1==1",20,-3.1,3.1,"#eta^{m1}","Lep/eta1_m",1);
292     draw_kin_variable("eta2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,-3.1,3.1,"#eta^{l1}","Lep/eta2",1);
293     draw_kin_variable("eta2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id2==0",20,-3.1,3.1,"#eta^{e2}","Lep/eta2_e",1);
294     draw_kin_variable("eta2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id2==1",20,-3.1,3.1,"#eta^{#mu2}","Lep/eta2_m",1);
295     draw_kin_variable("numVtx",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",30,0,30,"N(Vtx)","Lep/nVtx",1);
296     */
297 buchmann 1.1 }
298    
299     void draw_mpf_vars() {
300     ZbData::data_over_mc.clear();
301 buchmann 1.9 draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.3",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p3",0,"#alpha<0.3, 10 GeV < p_{T}^{Z} < 1000 GeV");
302     draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.2",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p2",0,"#alpha<0.2, 10 GeV < p_{T}^{Z} < 1000 GeV");
303     draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.15",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p15",0,"#alpha<0.15, 10 GeV < p_{T}^{Z} < 1000 GeV");
304     draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.1",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p1",0,"#alpha<0.1, 10 GeV < p_{T}^{Z} < 1000 GeV");
305     draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.05",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p05",0,"#alpha<0.05, 10 GeV < p_{T}^{Z} < 1000 GeV");
306    
307     draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.3",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p3",0,"#alpha<0.3, 10 GeV < p_{T}^{Z} < 1000 GeV");
308     draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.2",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p2",0,"#alpha<0.2, 10 GeV < p_{T}^{Z} < 1000 GeV");
309     draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.15",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p15",0,"#alpha<0.15, 10 GeV < p_{T}^{Z} < 1000 GeV");
310     draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.1",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p1",0,"#alpha<0.1, 10 GeV < p_{T}^{Z} < 1000 GeV");
311     draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.05",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz__mpf_alpha_smaller_0p05",0,"#alpha<0.05, 10 GeV < p_{T}^{Z} < 1000 GeV");
312 buchmann 1.1 }
313    
314 buchmann 1.12 /*double GetMedian(TH1F *histo) {
315 buchmann 1.11 int numBins = histo->GetXaxis()->GetNbins();
316 buchmann 1.12 Double_t x[numBins];
317     Double_t y[numBins];
318 buchmann 1.11 for (int i = 1; i <= numBins; i++) {
319     x[i] = histo->GetBinCenter(i);
320     y[i] = histo->GetBinContent(i);
321     }
322 buchmann 1.12 return TMath::Median(numBins, &x, &y);
323     }*/
324 buchmann 1.11
325 buchmann 1.15 Value get_Zb_data_over_mc(ZbResultContainer &res, string variable, string variable2, TCut cut, string saveas, float pt1, float pt2, float a1, float a2) {
326 buchmann 1.2 // write_warning(__FUNCTION__,"Debugging this function, therefore always returning 3.1415+/-1.010101"); return Value(3.1415,1.010101);
327 buchmann 1.12 TCanvas *cn = new TCanvas("cn","cn");
328     string varname="MPF";
329     if(Contains(variable,"pfJetGood")) varname="R_{abs}";
330     TH1F *hdata = allsamples.Draw("hdata",variable,1200,0,30, varname, "events", cut,data,luminosity);
331     TH1F *hmc = allsamples.Draw("hmc" , variable,1200,0,30, varname, "events", cut, mc, luminosity);
332 buchmann 1.14 hmc->SetLineColor(kBlue);
333    
334 buchmann 1.10 float a=hdata->GetMean();
335 buchmann 1.2 float b=hmc->GetMean();
336     float da=hdata->GetMeanError();
337 buchmann 1.10 float db=hmc->GetMeanError();
338 buchmann 1.2 float factor = a / b;
339     float error = TMath::Sqrt( (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
340 buchmann 1.12
341 buchmann 1.14 TF1 *gaus_data = new TF1("gaus_data","gaus",0.7,1.3);
342     gaus_data->SetParameter(0,hdata->GetMaximum());
343     gaus_data->SetParameter(2,hdata->GetMean());
344     gaus_data->SetParameter(1,hdata->GetRMS());
345     hdata->Fit(gaus_data);
346     TF1 *gaus_mc = new TF1("gaus_mc","gaus",0.7,1.3);
347     gaus_mc->SetParameter(0,hmc->GetMaximum());
348     gaus_mc->SetParameter(1,hmc->GetMean());
349     gaus_mc->SetParameter(2,hmc->GetRMS());
350     hmc->Fit(gaus_mc);
351    
352    
353 buchmann 1.12 cn->cd();
354     hdata->GetYaxis()->SetTitle("events / 0.1");
355 buchmann 1.14 hdata->SetMarkerColor(kRed);
356     /* hdata->Rebin(4);
357 buchmann 1.12 hmc->Rebin(4);
358 buchmann 1.14
359     gaus_mc->SetParameter(0,gaus_mc->GetParameter(0)*4);//rebinning
360     gaus_data->SetParameter(0,gaus_data->GetParameter(0)*4);//rebinning
361     */
362     gaus_mc->SetLineColor(kBlue);
363     gaus_data->SetLineColor(kRed);
364    
365 buchmann 1.12 hdata->GetXaxis()->SetRangeUser(0,2);
366     hdata->Draw("e1");
367     hmc->Draw("histo,same");
368     hdata->Draw("e1,same");
369 buchmann 1.14 gaus_data->Draw("same");
370     gaus_mc->Draw("same");
371 buchmann 1.12 stringstream summary;
372 buchmann 1.14 summary << "#splitline{" << varname << "}{#splitline{data: " << std::setprecision(4) << a << " #pm " << std::setprecision(4) << da << "}{";
373     summary << "#splitline{mc: " << std::setprecision(4) << b << " #pm " << std::setprecision(4) << db << "}{#splitline{ratio: " << factor;
374     summary << " #pm " << error << "}{#splitline{}{#splitline{data fit: " << std::setprecision(4) << gaus_data->GetParameter(1) << "+/-";
375     summary << std::setprecision(4) << gaus_data->GetParError(1) << "}{#splitline{";
376     summary << "mc fit:" << std::setprecision(4) << gaus_mc->GetParameter(1) << "+/-" << std::setprecision(4) << gaus_mc->GetParError(1) << "}{ratio: ";
377     summary << gaus_data->GetParameter(1)/gaus_mc->GetParameter(1) << "+/-" << gaus_data->GetParameter(1)/gaus_mc->GetParameter(1) * sqrt((gaus_data->GetParError(1) * gaus_data->GetParError(1))/(gaus_data->GetParameter(1)*gaus_data->GetParameter(1)) + (gaus_mc->GetParError(1) * gaus_mc->GetParError(1))/(gaus_mc->GetParameter(1)*gaus_mc->GetParameter(1)) ) << "}}}}}}}";
378 buchmann 1.12 TText *infobox = write_title(summary.str());
379     infobox->SetX(0.75);
380     infobox->SetTextSize(0.03);
381     infobox->SetY(0.75);
382     infobox->Draw();
383    
384    
385     DrawPrelim();
386     CompleteSave(cn,"ResponsePlots/"+saveas);
387    
388 buchmann 1.15 res.MPF_Rabs_location.push_back(ZbPointResult(variable,pt1,pt2,a1,a2,Value(a,da),Value(b,db)));
389 buchmann 1.12 cout << "Data;" << a << ";"<<da<<";MC"<<b<<";"<<db<<endl;
390     delete cn;
391 buchmann 1.2 delete hdata;
392     delete hmc;
393 buchmann 1.9
394 buchmann 1.2 return Value(factor,error);
395     }
396    
397 buchmann 1.15 ZbResultContainer new_data_mc_agreement_2d(string AlphaVariable, string ContainerName, TCut BTagCut, TCut BTagWgt) {
398     ZbResultContainer results(ContainerName);
399 buchmann 1.2
400 buchmann 1.15 cutWeight=TCut("(weight*(weight<1000)*(is_data+(!is_data)))")*BTagWgt;
401     LeadingB=BTagCut;
402 buchmann 1.13
403 buchmann 1.3 gStyle->SetOptFit(0);
404 buchmann 1.2 TCut basecut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
405 buchmann 1.15 // const int nalphacuts=6;
406     // float alphacuts[nalphacuts] = {0.1,0.15,0.2,0.25,0.3,0.35};
407    
408     // const int nptcuts=5;
409     // float ptcuts[nptcuts]={30,50,75,125,1000};
410    
411     write_info(__FUNCTION__,"JUST TESTING - PLEASE REMOVE THIS LINE");const int nalphacuts=2;float alphacuts[nalphacuts] = {0.1,0.3};const int nptcuts=3;float ptcuts[nptcuts]={30,75,1000};
412 buchmann 1.5
413 buchmann 1.2
414     float MPF_Results[nalphacuts][nptcuts];
415     float MPF_Errors[nalphacuts][nptcuts];
416     float RABS_Results[nalphacuts][nptcuts];
417     float RABS_Errors[nalphacuts][nptcuts];
418    
419    
420 buchmann 1.9 for(int ia=0;ia<nalphacuts;ia++) {
421 buchmann 1.2 for(int ipt=0;ipt<nptcuts-1;ipt++) {
422     stringstream specialcut;
423 buchmann 1.15 float alow,ahigh;
424     if(ia>0) {
425     specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << alphacuts[ia] << " && " << AlphaVariable << ">" << alphacuts[ia-1] << "))";
426     alow=alphacuts[ia-1];
427     ahigh=alphacuts[ia];
428     } else {
429     specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << alphacuts[ia] << "))";
430     alow=0;
431     ahigh=alphacuts[ia];
432     }
433 buchmann 1.9 cout << specialcut.str() << endl;
434 buchmann 1.15 Value MPF_data_over_mc = get_Zb_data_over_mc(results,"mpf","",TCut(basecut && specialcut.str().c_str()),"MPF___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__alpha_"+any2string(alphacuts[ia]), ptcuts[ipt], ptcuts[ipt+1],alow,ahigh);
435     Value RABS_data_over_mc = get_Zb_data_over_mc(results,"Zb3010_pfJetGoodPt[0]/pt","",TCut(basecut && specialcut.str().c_str()),"Rabs___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__alpha_"+any2string(alphacuts[ia]), ptcuts[ipt], ptcuts[ipt+1],alow,ahigh);
436 buchmann 1.2
437     MPF_Results[ia][ipt]=MPF_data_over_mc.getValue();
438     MPF_Errors[ia][ipt]=MPF_data_over_mc.getError();
439    
440     RABS_Results[ia][ipt]=RABS_data_over_mc.getValue();
441     RABS_Errors[ia][ipt]=RABS_data_over_mc.getError();
442    
443 buchmann 1.9 cout << "alpha cut at " << alphacuts[ia+1] << " and pt range " << ptcuts[ipt] << " , " << ptcuts[ipt+1] << " \t : \t" << specialcut.str() << " \t \t --> " << MPF_data_over_mc << " (MPF) " << RABS_data_over_mc << " (RABS)" << endl;
444 buchmann 1.2 }
445     }
446    
447     float MPF_ExtraPolatedResults[nptcuts];
448     float RABS_ExtraPolatedResults[nptcuts];
449    
450     TCanvas *can = new TCanvas("can");
451 buchmann 1.14 can->cd(1)->SetLeftMargin(0.15);
452 buchmann 1.2
453     TGraphErrors *MPF_FinalGraph = new TGraphErrors();
454     TGraphErrors *RABS_FinalGraph = new TGraphErrors();
455 buchmann 1.13 TGraphErrors *MPF_FaceValueAtPoint3 = new TGraphErrors();
456     TGraphErrors *RABS_FaceValueAtPoint3 = new TGraphErrors();
457 buchmann 1.2
458     for(int ipt=0;ipt<nptcuts-1;ipt++) {
459    
460     stringstream filename;
461     filename << "Extrapolation/Zplusb_data_over_mc___" << ptcuts[ipt] << "_to_" << ptcuts[ipt+1];
462     cout << "Going to create histo for " << filename.str() << endl;
463     TGraphErrors *mpf_gr =new TGraphErrors();
464     TGraphErrors *rabs_gr =new TGraphErrors();
465 buchmann 1.13
466     int rabs_counter=0;
467     int mpf_counter=0;
468    
469 buchmann 1.2 for(int ia=0;ia<nalphacuts;ia++) {
470 buchmann 1.9 cout << "Setting a point for an alpha cut at " << alphacuts[ia] << " with value " << MPF_Results[ia][ipt] << " (MPF) and " << RABS_Results[ia][ipt] << " (RABS)" << endl;
471 buchmann 1.13 if(MPF_Results[ia][ipt]==MPF_Results[ia][ipt] && MPF_Errors[ia][ipt]==MPF_Errors[ia][ipt]) { // remove nan's
472     mpf_gr->SetPoint(mpf_counter,alphacuts[ia],MPF_Results[ia][ipt]);
473     mpf_gr->SetPointError(mpf_counter,0,MPF_Errors[ia][ipt]);
474     mpf_counter++;
475     } else {
476     dout << "Detected BS for MPF method (bin ignored) : " << endl;
477     dout << " alpha: " << alphacuts[ia] << " , pt range: " << ptcuts[ipt] << " < pt < " << ptcuts[ipt+1] << endl;
478     cout << " value (MPF) : " << MPF_Results[ia][ipt] << " +/- " << MPF_Errors[ia][ipt] << endl;
479     }
480     if(RABS_Results[ia][ipt]==RABS_Results[ia][ipt] && RABS_Errors[ia][ipt]==RABS_Errors[ia][ipt]) {
481     rabs_gr->SetPoint(rabs_counter,alphacuts[ia],RABS_Results[ia][ipt]);
482     rabs_gr->SetPointError(rabs_counter,0,RABS_Errors[ia][ipt]);
483     rabs_counter++;
484     } else {
485     dout << "Detected BS for RABS method (bin ignored) : " << endl;
486     dout << " alpha: " << alphacuts[ia] << " , pt range: " << ptcuts[ipt] << " < pt < " << ptcuts[ipt+1] << endl;
487     cout << " value (RABS) : " << RABS_Results[ia][ipt] << " +/- " << RABS_Errors[ia][ipt] << endl;
488     }
489     }
490    
491     stringstream specialcut;
492 buchmann 1.15 specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<0.3))";
493     Value MPF_data_over_mc = get_Zb_data_over_mc(results,"mpf","",TCut(basecut && specialcut.str().c_str()),"MPF___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__INCLUSIVEalpha_0.3",ptcuts[ipt],ptcuts[ipt+1],0,0.3);
494     Value RABS_data_over_mc = get_Zb_data_over_mc(results,"Zb3010_pfJetGoodPt[0]/pt","",TCut(basecut && specialcut.str().c_str()),"Rabs___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__INCLUSIVEalpha_0.3",ptcuts[ipt],ptcuts[ipt+1],0,0.3);
495 buchmann 1.14 cout << "About to add a point in pt (" << (ptcuts[ipt]+ptcuts[ipt+1])/2 << " ) to face value : " << MPF_data_over_mc << " (MPF) , " << RABS_data_over_mc << "(RABS)" << endl;
496 buchmann 1.13 MPF_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,MPF_data_over_mc.getValue());
497     MPF_FaceValueAtPoint3->SetPointError(ipt,0,MPF_data_over_mc.getError());
498     RABS_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,RABS_data_over_mc.getValue());
499     RABS_FaceValueAtPoint3->SetPointError(ipt,0,RABS_data_over_mc.getError());
500 buchmann 1.2
501 buchmann 1.13 can->cd();
502 buchmann 1.2 mpf_gr->SetMarkerStyle(21);
503     mpf_gr->SetMarkerColor(kBlue);
504     mpf_gr->Draw("AP");
505 buchmann 1.13 mpf_gr->Fit("pol1","");
506 buchmann 1.3 mpf_gr->GetXaxis()->SetTitle("#alpha");
507     mpf_gr->GetXaxis()->CenterTitle();
508 buchmann 1.9 mpf_gr->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
509 buchmann 1.3 mpf_gr->GetYaxis()->CenterTitle();
510 buchmann 1.2 TF1 *mpf_pol=(TF1*)mpf_gr->GetFunction("pol1");
511 buchmann 1.14 mpf_pol->SetLineWidth(0);
512 buchmann 1.13 TF1 *mpf_pol_complete = new TF1("mpf_pol_complete","[0]+[1]*x");
513     mpf_pol_complete->SetParameters(mpf_pol->GetParameters());
514     mpf_pol_complete->SetLineWidth(1);
515     mpf_pol_complete->SetLineColor(kBlue);
516     float min,max;
517     FindMinMax(mpf_gr,min,max);
518     TH1F *histo = new TH1F("histo","histo",1,0,0.4);
519 buchmann 1.14 histo->GetXaxis()->SetTitle("#alpha");
520     histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
521     histo->GetXaxis()->CenterTitle();
522     histo->GetYaxis()->CenterTitle();
523 buchmann 1.13 histo->SetStats(0);
524     histo->GetXaxis()->SetRangeUser(0,0.4);
525     histo->GetYaxis()->SetRangeUser(min,max);
526     mpf_gr->GetYaxis()->SetRangeUser(min,max);
527    
528     TPolyLine *mpf_fit_uncert = GetFitUncertaintyShape(mpf_gr, "pol1", min, max,0.0,0.4);
529     mpf_pol->SetLineColor(TColor::GetColor("#0101DF"));
530 buchmann 1.14 mpf_fit_uncert->SetFillColor(TColor::GetColor("#5353FC"));
531     histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
532 buchmann 1.13 histo->Draw();
533 buchmann 1.14 mpf_fit_uncert->Draw("F");
534     mpf_fit_uncert->Draw("");
535 buchmann 1.13 mpf_gr->Draw("P");
536     mpf_pol_complete->Draw("same");
537    
538 buchmann 1.2 float mpf_a=mpf_pol->GetParameter(0);
539     float mpf_b=mpf_pol->GetParameter(1);
540     MPF_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),mpf_a);
541     MPF_FinalGraph->SetPointError(ipt,0,mpf_pol->GetParError(0));
542    
543     MPF_ExtraPolatedResults[ipt]=mpf_a;
544    
545     stringstream mpf_info;
546 buchmann 1.14 mpf_info << "#splitline{#splitline{#splitline{" << ptcuts[ipt] << " GeV < p_{T}^{Z} < " << ptcuts[ipt+1] << " GeV }{ (MPF) }}{Extrapolated value: " << std::setprecision(4) << mpf_a << "}}{#chi^{2}/NDF : " << mpf_pol->GetChisquare() << " / " << mpf_pol->GetNDF() << "}";
547 buchmann 1.3
548 buchmann 1.2 TText *mpf_mark = write_title(mpf_info.str());
549     mpf_mark->SetX(0.75);
550     mpf_mark->SetTextSize(0.03);
551     mpf_mark->SetY(0.75);
552     mpf_mark->Draw();
553    
554     string filenamebkp=filename.str();
555     filename << "__MPF";
556 buchmann 1.3 DrawPrelim();
557 buchmann 1.15 CompleteSave(can,ContainerName+"/"+filename.str());
558 buchmann 1.2 cout << "MPF : " << mpf_a << " + " << mpf_b << " * alpha " << endl;
559    
560     filename.str("");
561    
562     rabs_gr->SetMarkerStyle(21);
563 buchmann 1.13 rabs_gr->SetMarkerColor(kRed);
564 buchmann 1.2 rabs_gr->Draw("AP");
565 buchmann 1.13 rabs_gr->Fit("pol1","");
566 buchmann 1.3 rabs_gr->GetXaxis()->SetTitle("#alpha");
567     rabs_gr->GetXaxis()->CenterTitle();
568 buchmann 1.14 rabs_gr->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
569 buchmann 1.3 rabs_gr->GetYaxis()->CenterTitle();
570 buchmann 1.13 TF1 *rabs_pol=(TF1*)rabs_gr->GetFunction("pol1");
571 buchmann 1.14 rabs_pol->SetLineWidth(0);
572 buchmann 1.13 TF1 *rabs_pol_complete = new TF1("rabs_pol_complete","[0]+[1]*x");
573     rabs_pol_complete->SetParameters(rabs_pol->GetParameters());
574     rabs_pol_complete->SetLineColor(kRed);
575     rabs_pol_complete->SetLineWidth(1);
576     FindMinMax(rabs_gr,min,max);
577     rabs_gr->GetYaxis()->SetRangeUser(min,max);
578     histo->GetYaxis()->SetRangeUser(min,max);
579     TPolyLine *rabs_fit_uncert = GetFitUncertaintyShape(rabs_gr, "pol1", min, max,0.0,0.4);
580 buchmann 1.14 rabs_pol->SetLineColor(TColor::GetColor("#FA5858"));
581     rabs_pol->SetFillColor(TColor::GetColor("#FA5858"));
582     rabs_fit_uncert->SetLineColor(TColor::GetColor("#FA5858"));
583     rabs_fit_uncert->SetFillColor(TColor::GetColor("#FA5858"));
584     histo->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
585 buchmann 1.13 histo->Draw();
586 buchmann 1.14 rabs_fit_uncert->Draw("F");
587     rabs_fit_uncert->Draw("");
588 buchmann 1.13 rabs_gr->Draw("P");
589     rabs_pol_complete->Draw("same");
590    
591 buchmann 1.3
592 buchmann 1.2 float rabs_a=rabs_pol->GetParameter(0);
593     float rabs_b=rabs_pol->GetParameter(1);
594     RABS_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),rabs_a);
595     RABS_FinalGraph->SetPointError(ipt,0,rabs_pol->GetParError(0));
596 buchmann 1.3
597 buchmann 1.2 cout << "!+!+!+!+!+!+!+!+!+ Just added a point to the final plot : " << rabs_a << " +/- " << rabs_pol->GetParError(0) << endl;
598    
599     stringstream rabs_info;
600 buchmann 1.14 rabs_info << "#splitline{#splitline{#splitline{" << ptcuts[ipt] << " GeV < p_{T}^{Z} < " << ptcuts[ipt+1] << " GeV }{ (R_{abs}) }}{Extrapolated value: " << std::setprecision(4) << rabs_a << "}}{#chi^{2}/NDF : " << rabs_pol->GetChisquare() << " / " << rabs_pol->GetNDF() << "}";
601 buchmann 1.2 TText *rabs_mark = write_title(rabs_info.str());
602     rabs_mark->SetX(0.75);
603     rabs_mark->SetTextSize(0.03);
604     rabs_mark->SetY(0.75);
605     rabs_mark->Draw();
606 buchmann 1.3
607    
608 buchmann 1.2 RABS_ExtraPolatedResults[ipt]=rabs_a;
609     filename << filenamebkp << "__RABS";
610 buchmann 1.3 DrawPrelim();
611 buchmann 1.15 CompleteSave(can,ContainerName+"/"+filename.str());
612 buchmann 1.2 cout << "RABS : " << rabs_a << " + " << rabs_b << " * alpha " << endl;
613     }
614    
615     MPF_FinalGraph->SetMarkerStyle(21);
616 buchmann 1.13 MPF_FinalGraph->SetMarkerColor(kBlue);
617     MPF_FinalGraph->Fit("pol0",""); // quiet, use minos
618 buchmann 1.3
619 buchmann 1.2 TF1 *mpf_pol0=(TF1*)MPF_FinalGraph->GetFunction("pol0");
620 buchmann 1.13 TF1 *mpf_pol0_complete = new TF1("mpf_pol0_complete","[0]+[1]*x");
621     mpf_pol0_complete->SetLineWidth(1);
622     mpf_pol0_complete->SetLineColor(kBlue);
623     mpf_pol0->SetLineColor(kBlue);
624 buchmann 1.2 MPF_FinalGraph->Draw("AP");
625 buchmann 1.3 MPF_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
626     MPF_FinalGraph->GetXaxis()->CenterTitle();
627 buchmann 1.9 MPF_FinalGraph->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
628 buchmann 1.3 MPF_FinalGraph->GetYaxis()->CenterTitle();
629     MPF_FinalGraph->Draw("AP");
630 buchmann 1.13 mpf_pol0_complete->Draw("same");
631 buchmann 1.3
632 buchmann 1.2 stringstream mpf_result;
633 buchmann 1.3 mpf_result << "#splitline{C_{abs}= " << std::setprecision(5) << 1/mpf_pol0->GetParameter(0) << " #pm " << std::setprecision(2) << mpf_pol0->GetParError(0)/(mpf_pol0->GetParameter(0)*mpf_pol0->GetParameter(0)) << "}{ (MPF)}";
634 buchmann 1.2 TText *rmpf_final_mark = write_title(mpf_result.str());
635     rmpf_final_mark->SetX(0.75);
636     rmpf_final_mark->SetY(0.75);
637     rmpf_final_mark->SetTextSize(0.03);
638     rmpf_final_mark->Draw();
639 buchmann 1.3 DrawPrelim();
640 buchmann 1.5 float MPFResult=1/mpf_pol0->GetParameter(0);
641     float MPFResultError=mpf_pol0->GetParError(0)/(mpf_pol0->GetParameter(0)*mpf_pol0->GetParameter(0));
642 buchmann 1.2
643     stringstream filename2;
644     filename2 << "Extrapolation/FINAL_RESULT_MPF";
645 buchmann 1.15 CompleteSave(can,ContainerName+"/"+filename2.str());
646 buchmann 1.2
647    
648     RABS_FinalGraph->SetMarkerStyle(21);
649     RABS_FinalGraph->SetMarkerStyle(21);
650     RABS_FinalGraph->SetMarkerColor(kRed);
651 buchmann 1.3 RABS_FinalGraph->Fit("pol0","QE"); // quiet, use minos
652     RABS_FinalGraph->Draw("AP");
653     RABS_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
654     RABS_FinalGraph->GetXaxis()->CenterTitle();
655 buchmann 1.9 RABS_FinalGraph->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
656 buchmann 1.3 RABS_FinalGraph->GetYaxis()->CenterTitle();
657 buchmann 1.2 RABS_FinalGraph->Draw("AP");
658     TF1 *rabs_pol0=(TF1*)RABS_FinalGraph->GetFunction("pol0");
659     stringstream rabs_result;
660 buchmann 1.3 rabs_result << "#splitline{C_{abs}= " << std::setprecision(5) << 1/rabs_pol0->GetParameter(0) << " #pm " << std::setprecision(2) << rabs_pol0->GetParError(0)/(rabs_pol0->GetParameter(0)*rabs_pol0->GetParameter(0)) << "}{ (R_{abs})}";;
661 buchmann 1.2 TText *rabs_final_mark = write_title(rabs_result.str());
662     rabs_final_mark->SetX(0.75);
663     rabs_final_mark->SetY(0.75);
664     rabs_final_mark->SetTextSize(0.03);
665     rabs_final_mark->Draw();
666 buchmann 1.3 DrawPrelim();
667 buchmann 1.2
668 buchmann 1.5 float RabsResult=1/rabs_pol0->GetParameter(0);
669     float RabsResultError=rabs_pol0->GetParError(0)/(rabs_pol0->GetParameter(0)*rabs_pol0->GetParameter(0));
670    
671     filename2.str("");
672     filename2 << "Extrapolation/FINAL_RESULT_RABS";
673 buchmann 1.15 CompleteSave(can,ContainerName+"/"+filename2.str());
674 buchmann 1.5
675 buchmann 1.9 can->cd();
676 buchmann 1.13 MPF_FaceValueAtPoint3->SetMarkerStyle(21);
677     MPF_FaceValueAtPoint3->SetMarkerColor(kBlue);
678     MPF_FaceValueAtPoint3->Fit("pol0","QE"); // quiet, use minos
679     MPF_FaceValueAtPoint3->Draw("AP");
680     MPF_FaceValueAtPoint3->GetXaxis()->SetTitle("p_{T}^{Z}");
681     MPF_FaceValueAtPoint3->GetXaxis()->CenterTitle();
682     MPF_FaceValueAtPoint3->GetYaxis()->SetTitle("<MPF>__{data}/<MPF>_{mc} for #alpha<0.3");
683     MPF_FaceValueAtPoint3->GetYaxis()->CenterTitle();
684     MPF_FaceValueAtPoint3->Draw("AP");
685     TF1 *faceval_pol0=(TF1*)MPF_FaceValueAtPoint3->GetFunction("pol0");
686     faceval_pol0->SetLineColor(kBlue);
687     faceval_pol0->SetLineWidth(1);
688     faceval_pol0->Draw("same");
689     TF1 *faceval_pol0_complete = new TF1("faceval_pol0_complete","[0]+[1]*x");
690     faceval_pol0_complete->SetParameters(faceval_pol0->GetParameters());
691 buchmann 1.5 stringstream faceval_result;
692     faceval_result << "#splitline{C_{abs}= " << std::setprecision(5) << 1/faceval_pol0->GetParameter(0) << " #pm " << std::setprecision(2) << faceval_pol0->GetParError(0)/(faceval_pol0->GetParameter(0)*faceval_pol0->GetParameter(0)) << "}{ (MPF at #alpha<0.3)}";;
693     TText *faceval_final_mark = write_title(faceval_result.str());
694     faceval_final_mark->SetX(0.75);
695     faceval_final_mark->SetY(0.75);
696     faceval_final_mark->SetTextSize(0.03);
697     faceval_final_mark->Draw();
698     DrawPrelim();
699    
700     float FaceValResult=1/faceval_pol0->GetParameter(0);
701     float FaceValResultError=faceval_pol0->GetParError(0)/(faceval_pol0->GetParameter(0)*faceval_pol0->GetParameter(0));
702    
703 buchmann 1.2 filename2.str("");
704 buchmann 1.9 filename2 << "Extrapolation/FINAL_RESULT_MPF_FaceValp3";
705 buchmann 1.15 CompleteSave(can,ContainerName+"/"+filename2.str());
706 buchmann 1.2
707 buchmann 1.13 can->cd();
708     RABS_FaceValueAtPoint3->SetMarkerStyle(21);
709     RABS_FaceValueAtPoint3->SetMarkerColor(kRed);
710     RABS_FaceValueAtPoint3->Fit("pol0","QE"); // quiet, use minos
711     RABS_FaceValueAtPoint3->Draw("AP");
712     RABS_FaceValueAtPoint3->GetXaxis()->SetTitle("p_{T}^{Z}");
713     RABS_FaceValueAtPoint3->GetXaxis()->CenterTitle();
714     RABS_FaceValueAtPoint3->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc} for #alpha<0.3");
715     RABS_FaceValueAtPoint3->GetYaxis()->CenterTitle();
716     RABS_FaceValueAtPoint3->Draw("AP");
717     TF1 *faceval_pol0RABS=(TF1*)RABS_FaceValueAtPoint3->GetFunction("pol0");
718     TF1 *faceval_pol0RABS_complete = new TF1("faceval_pol0_complete","[0]+[1]*x");
719     faceval_pol0RABS_complete->SetParameters(faceval_pol0RABS->GetParameters());
720     faceval_pol0RABS->SetLineWidth(1);
721     faceval_pol0RABS->SetLineColor(kRed);
722     faceval_pol0RABS_complete->Draw("same");
723     stringstream RABSfaceval_result;
724     RABSfaceval_result << "#splitline{C_{abs}= " << std::setprecision(5) << 1/faceval_pol0RABS->GetParameter(0) << " #pm " << std::setprecision(2) << faceval_pol0RABS->GetParError(0)/(faceval_pol0RABS->GetParameter(0)*faceval_pol0RABS->GetParameter(0)) << "}{ (R_{abs} at #alpha<0.3)}";;
725     TText *RABSfaceval_final_mark = write_title(RABSfaceval_result.str());
726     RABSfaceval_final_mark->SetX(0.75);
727     RABSfaceval_final_mark->SetY(0.75);
728     RABSfaceval_final_mark->SetTextSize(0.03);
729     RABSfaceval_final_mark->Draw();
730     DrawPrelim();
731    
732     float RABSFaceValResult=1/faceval_pol0RABS->GetParameter(0);
733     float RABSFaceValResultError=faceval_pol0RABS->GetParError(0)/(faceval_pol0RABS->GetParameter(0)*faceval_pol0RABS->GetParameter(0));
734    
735     filename2.str("");
736     filename2 << "Extrapolation/FINAL_RESULT_RABS_FaceValp3";
737 buchmann 1.15 CompleteSave(can,ContainerName+"/"+filename2.str());
738 buchmann 1.13
739 buchmann 1.5 cout << "FINAL RESULTS: " << endl;
740     cout << "MPF : " << MPFResult << " +/- " << MPFResultError << endl;
741     cout << "Rabs: " << RabsResult << " +/- " << RabsResultError << endl;
742 buchmann 1.13 cout << "Face value at 0.3 :" << FaceValResult << " +/- " << FaceValResultError << " (MPF) " << endl;
743     cout << "Face value at 0.3 :" << RABSFaceValResult << " +/- " << RABSFaceValResultError << " (RABS) " << endl;
744 buchmann 1.2
745     delete can;
746    
747 buchmann 1.15 write_info(__FUNCTION__,"Need to fill result capsule");
748     results.MPF_Result_FaceValue=Value(FaceValResult,FaceValResultError);
749     results.Rabs_Result_FaceValue=Value(RABSFaceValResult,RABSFaceValResultError);
750     results.MPF_Result_Extrapolated=Value(MPFResult,MPFResultError);
751     results.Rabs_Result_Extrapolated=Value(RabsResult,RabsResultError);
752     cout << results << endl;
753     return results;
754 buchmann 1.2 }
755    
756    
757 buchmann 1.15 void TEST_DoPUStudy(string identifier) {
758 buchmann 1.9
759     float numVtxcuts[7]={0,10,15,20};
760    
761     TCanvas *can = new TCanvas("can","can");
762     TH1F *malpha[8];
763     TH1F *dalpha[8];
764     cout << identifier << ";";
765     cout << endl;
766     for(int i=1;i<5;i++) {
767     stringstream sSpecialCut;
768     if(i<4) {
769     sSpecialCut << "numVtx>" << numVtxcuts[i-1] << " && numVtx<" << numVtxcuts[i];
770     cout << numVtxcuts[i-1] << " < nVtx < " << numVtxcuts[i] << ";";
771     } else {
772     sSpecialCut << "numVtx>" << numVtxcuts[i-1];
773     cout << numVtxcuts[i-1] << ";";
774     }
775    
776     TCut SpecialCut(sSpecialCut.str().c_str());
777    
778    
779     malpha[i] = allsamples.Draw("malpha"+any2string(i),"mpf",1000,0,10,"MPF","events",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&SpecialCut,mc,luminosity);
780     dalpha[i] = allsamples.Draw("dalpha"+any2string(i),"mpf",1000,0,10,"MPF","events",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&SpecialCut,data,luminosity);
781    
782    
783     float a=dalpha[i]->GetMean();
784     float b=malpha[i]->GetMean();
785     float da=dalpha[i]->GetMeanError();
786     float db=malpha[i]->GetMeanError();
787     float factor = a / b;
788     float error = TMath::Sqrt( (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
789    
790     cout << malpha[i]->GetMean() << ";" << malpha[i]->GetMeanError() << ";" << dalpha[i]->GetMean() << ";" << dalpha[i]->GetMeanError() << ";" << malpha[i]->Integral()<<";"<<dalpha[i]->Integral() << ";"<<factor << ";" << error << endl;
791     malpha[i]->SetLineColor(i+1);
792     if(i==0) malpha[i]->Draw("histo");
793     else malpha[i]->Draw("histo,same");
794     }
795    
796     CompleteSave(can,"AlphaOverview_"+identifier);
797     delete can;
798    
799     /* for(int i=1;i<8;i++) {
800     delete malpha[i];
801     delete dalpha[i];
802     }
803     */
804    
805     }
806    
807 buchmann 1.15 void TEST_ScenarioComparison() {
808 buchmann 1.13 //very dumb way to do this.
809 buchmann 1.9 TGraphAsymmErrors *gr[5];
810     TF1 *fit[5];
811     bool dofit=false;
812    
813     TCanvas *can = new TCanvas("can","can",500,500);
814    
815     TLegend *leg = new TLegend(0.2,0.2,0.5,0.45);
816     leg->SetBorderSize(0);
817     leg->SetFillColor(kWhite);
818    
819     for(int i=0;i<5;i++) {
820     gr[i] = new TGraphAsymmErrors();
821     double x[4] = {5,12.5,17.5,20};
822     double dxl[4] = {5,2.5,2.5,2.5};
823     double dxh[4] = {5,2.5,2.5,10};
824    
825     string name="";
826     string color="";
827    
828     if(i==0) {
829     double y[4] = {1.008,1.013,1.003,1.007};
830     double dyl[4] = {0.011,0.009,0.013,0.020};
831     double dyh[4] = {0.011,0.009,0.013,0.020};
832     name="30/30";
833     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
834     color="#a01d99";
835     }
836     if(i==1) {
837     double y[4] = {1.008+0.001,1.013+0.001,1.003+0.001,1.007+0.001};
838     double dyl[4] = {0.011,0.009,0.013,0.020};
839     double dyh[4] = {0.011,0.009,0.013,0.020};
840     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
841 buchmann 1.13 name="30/15";
842 buchmann 1.9 color="#2464bc";
843     }
844     if(i==2) {
845     double y[4] = {1.014,1.026,1.014,1.059};
846     double dyl[4] = {0.011,0.010,0.014,0.023};
847     double dyh[4] = {0.011,0.010,0.014,0.023};
848     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
849     name="15/15";
850     color="#f49230";
851     }
852     if(i==3) {
853     double y[4] = {1.015,1.020,1.008,1.050};
854     double dyl[4] = {0.011,0.010,0.014,0.022};
855     double dyh[4] = {0.011,0.010,0.014,0.022};
856     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
857     name="20/15";
858     color="#c72f3c";
859     }
860     if(i==4) {
861     double y[4] = {1.015+0.001,1.020+0.001,1.008+0.001,1.050+0.001};
862     double dyl[4] = {0.011,0.010,0.014,0.022};
863     double dyh[4] = {0.011,0.010,0.014,0.022};
864     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
865     name="20/20";
866     color="#7cb433";
867     }
868    
869    
870     gr[i]->SetTitle();
871     gr[i]->GetXaxis()->SetTitle("N(Vtx)");
872     gr[i]->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
873     gr[i]->GetXaxis()->CenterTitle();
874     gr[i]->GetYaxis()->CenterTitle();
875     gr[i]->SetName(name.c_str());
876     gr[i]->SetMarkerSize(0.00001);
877     gr[i]->GetYaxis()->SetRangeUser(0.90,1.10);
878     gr[i]->SetLineColor(TColor::GetColor(color.c_str()));
879     gr[i]->SetMarkerColor(TColor::GetColor(color.c_str()));
880     if(dofit) {
881     gr[i]->Fit("pol1");
882     fit[i] = (TF1*)gr[i]->GetFunction("pol1");
883     fit[i]->SetLineColor(TColor::GetColor(color.c_str()));
884     }
885     leg->AddEntry(gr[i],name.c_str(),"l");
886     if(i==0) gr[i]->Draw("AP*");
887     else gr[i]->Draw("P");
888     }
889     leg->Draw();
890     DrawPrelim();
891    
892     CompleteSave(can,"ScenarioComparison");
893     delete can;
894    
895     }
896    
897 buchmann 1.15 void TEST_ScenarioComparisonInclusive() {
898 buchmann 1.13 //dumbest way ever to do this. but ok we only need to do it once.
899 buchmann 1.9 TGraphAsymmErrors *gr[5];
900     TF1 *fit[5];
901     bool dofit=true;
902    
903     TCanvas *can = new TCanvas("can","can",500,500);
904    
905     TLegend *leg = new TLegend(0.2,0.2,0.5,0.45);
906     leg->SetBorderSize(0);
907     leg->SetFillColor(kWhite);
908    
909     for(int i=0;i<5;i++) {
910     gr[i] = new TGraphAsymmErrors();
911     double x[4] = {5,12.5,17.5,20};
912     double dxl[4] = {5,2.5,2.5,2.5};
913     double dxh[4] = {5,2.5,2.5,10};
914    
915     string name="";
916     string color="";
917    
918     if(i==0) {
919     double y[4] = {1.01103,1.01496,1.01486,1.02785};
920     double dyl[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
921     double dyh[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
922     name="30/30";
923     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
924     color="#a01d99";
925     }
926     if(i==1) {
927     double y[4] = {1.01103+0.001,1.01496+0.001,1.01486+0.001,1.02785+0.001};
928     double dyl[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
929     double dyh[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
930     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
931 buchmann 1.13 name="30/15";
932 buchmann 1.9 color="#2464bc";
933     }
934     if(i==2) {
935     double y[4] = {1.01575,1.02879,1.0325,1.05456};
936     double dyl[4] = {0.00288281,0.00278536,0.00367451,0.00590421};
937     double dyh[4] = {0.00288281,0.00278536,0.00367451,0.00590421};
938     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
939     name="15/15";
940     color="#f49230";
941     }
942     if(i==3) {
943     double y[4] = {1.01221,1.01922,1.02396,1.04677};
944     double dyl[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
945     double dyh[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
946     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
947     name="20/15";
948     color="#c72f3c";
949     }
950     if(i==4) {
951     double y[4] = {1.01221+0.001,1.01922+0.001,1.02396+0.001,1.04677+0.001};
952     double dyl[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
953     double dyh[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
954     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
955     name="20/20";
956     color="#7cb433";
957     }
958    
959    
960     gr[i]->SetTitle();
961     gr[i]->GetXaxis()->SetTitle("N(Vtx)");
962     gr[i]->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
963     gr[i]->GetXaxis()->CenterTitle();
964     gr[i]->GetYaxis()->CenterTitle();
965     gr[i]->SetName(name.c_str());
966     gr[i]->SetMarkerSize(0.00001);
967     gr[i]->GetYaxis()->SetRangeUser(0.90,1.10);
968     gr[i]->SetLineColor(TColor::GetColor(color.c_str()));
969     gr[i]->SetMarkerColor(TColor::GetColor(color.c_str()));
970     if(dofit) {
971     gr[i]->Fit("pol1");
972     fit[i] = (TF1*)gr[i]->GetFunction("pol1");
973     fit[i]->SetLineColor(TColor::GetColor(color.c_str()));
974     }
975     leg->AddEntry(gr[i],name.c_str(),"l");
976     if(i==0) gr[i]->Draw("AP*");
977     else gr[i]->Draw("P");
978     }
979     leg->Draw();
980     DrawPrelim();
981    
982     CompleteSave(can,"ScenarioComparisonInclusive");
983     delete can;
984    
985     }
986 buchmann 1.2
987 buchmann 1.15 void TEST_compare_selection(string identifier) {
988 buchmann 1.7 bool recognized_scenario=false;
989    
990     cout << "Running with identifier " << identifier << endl;
991 buchmann 1.9 if(identifier=="Zb3010") {
992     LeadingB=TCut ("Zb3010_bTagProbCSVBP[0]>0.898");
993     EtaB=TCut("abs(Zb3010_pfJetGoodEta[0])<1.3");
994     PhiZcut=TCut("abs(Zb3010_pfJetDphiZ[0])>2.7");
995     recognized_scenario=true;
996     }
997    
998     if(identifier=="Zb3010") {
999     LeadingB=TCut ("Zb3010_bTagProbCSVBP[0]>0.898");
1000     EtaB=TCut("abs(Zb3010_pfJetGoodEta[0])<1.3");
1001     PhiZcut=TCut("abs(Zb3010_pfJetDphiZ[0])>2.7");
1002 buchmann 1.7 recognized_scenario=true;
1003     }
1004    
1005 buchmann 1.9 if(identifier=="Zb40") {
1006     LeadingB=TCut ("Zb40_bTagProbCSVBP[0]>0.898");
1007     EtaB=TCut("abs(Zb40_pfJetGoodEta[0])<1.3");
1008     PhiZcut=TCut("abs(Zb40_pfJetDphiZ[0])>2.7");
1009     recognized_scenario=true;
1010     }
1011    
1012     if(identifier=="Zb301030") {
1013     LeadingB=TCut ("Zb301030_bTagProbCSVBP[0]>0.898");
1014     EtaB=TCut("abs(Zb301030_pfJetGoodEta[0])<1.3");
1015     PhiZcut=TCut("abs(Zb301030_pfJetDphiZ[0])>2.7");
1016 buchmann 1.7 recognized_scenario=true;
1017     }
1018    
1019     if(identifier=="Zb1530") {
1020 buchmann 1.9 LeadingB=TCut ("Zb1530_bTagProbCSVBP[0]>0.898");
1021 buchmann 1.7 EtaB=TCut("abs(Zb1530_pfJetGoodEta[0])<1.3");
1022     PhiZcut=TCut("abs(Zb1530_pfJetDphiZ[0])>2.7");
1023     recognized_scenario=true;
1024     }
1025    
1026 buchmann 1.9 if(identifier=="Zb301010") {
1027     LeadingB=TCut ("Zb301010_bTagProbCSVBP[0]>0.898");
1028     EtaB=TCut("abs(Zb301010_pfJetGoodEta[0])<1.3");
1029     PhiZcut=TCut("abs(Zb301010_pfJetDphiZ[0])>2.7");
1030     recognized_scenario=true;
1031     }
1032    
1033     if(identifier=="Zb1510") {
1034     LeadingB=TCut ("Zb1510_bTagProbCSVBP[0]>0.898");
1035     EtaB=TCut("abs(Zb1510_pfJetGoodEta[0])<1.3");
1036     PhiZcut=TCut("abs(Zb1510_pfJetDphiZ[0])>2.7");
1037     recognized_scenario=true;
1038     }
1039    
1040     if(identifier=="Zb301010") {
1041     LeadingB=TCut ("Zb301010_bTagProbCSVBP[0]>0.898");
1042     EtaB=TCut("abs(Zb301010_pfJetGoodEta[0])<1.3");
1043     PhiZcut=TCut("abs(Zb301010_pfJetDphiZ[0])>2.7");
1044     recognized_scenario=true;
1045     }
1046    
1047     if(identifier=="ZbMikko") {
1048     LeadingB=TCut ("ZbMikko_bTagProbCSVBP[0]>0.898 && ZbMikko__IsClean");
1049     EtaB=TCut("abs(ZbMikko_pfJetGoodEta[0])<1.3 && ZbMikko__IsClean");
1050     PhiZcut=TCut("abs(ZbMikko_pfJetDphiZ[0])>2.7 && ZbMikko__IsClean");
1051 buchmann 1.7 recognized_scenario=true;
1052     }
1053    
1054 buchmann 1.9 if(identifier=="Zb3010_SecEta3") {
1055     LeadingB=TCut ("Zb3010_SecEta3_bTagProbCSVBP[0]>0.898");
1056     EtaB=TCut("abs(Zb3010_SecEta3_pfJetGoodEta[0])<1.3");
1057     PhiZcut=TCut("abs(Zb3010_SecEta3_pfJetDphiZ[0])>2.7");
1058 buchmann 1.7 recognized_scenario=true;
1059     }
1060    
1061 buchmann 1.9 if(identifier=="Zb3010_SecEta5") {
1062     LeadingB=TCut ("Zb3010_SecEta5_bTagProbCSVBP[0]>0.898");
1063     EtaB=TCut("abs(Zb3010_SecEta5_pfJetGoodEta[0])<1.3");
1064     PhiZcut=TCut("abs(Zb3010_SecEta5_pfJetDphiZ[0])>2.7");
1065 buchmann 1.7 recognized_scenario=true;
1066     }
1067    
1068 buchmann 1.9 if(identifier=="Zb3010_p5Clean") {
1069     LeadingB=TCut ("Zb3010_p5Clean_bTagProbCSVBP[0]>0.898 && Zb3010_p5Clean_IsClean");
1070     EtaB=TCut("abs(Zb3010_p5Clean_pfJetGoodEta[0])<1.3 && Zb3010_p5Clean_IsClean");
1071     PhiZcut=TCut("abs(Zb3010_p5Clean_pfJetDphiZ[0])>2.7 && Zb3010_p5Clean_IsClean");
1072 buchmann 1.7 recognized_scenario=true;
1073     }
1074    
1075 buchmann 1.9 if(identifier=="Zb3010CHS") {
1076     LeadingB=TCut ("Zb3010CHS_bTagProbCSVBP[0]>0.898");
1077     EtaB=TCut("abs(Zb3010CHS_pfJetGoodEta[0])<1.3");
1078     PhiZcut=TCut("abs(Zb3010CHS_pfJetDphiZ[0])>2.7");
1079 buchmann 1.7 recognized_scenario=true;
1080     }
1081    
1082     assert(recognized_scenario);
1083    
1084     cout << "Cuts: " << endl;
1085     cout << " " << (const char*) LeadingB << endl;
1086     cout << " " << (const char*) EtaB << endl;
1087     cout << " " << (const char*) PhiZcut << endl;
1088    
1089     cout << endl << endl;
1090    
1091 buchmann 1.9 // ScenarioComparison();
1092     // ScenarioComparisonInclusive();
1093    
1094    
1095     // DoPUStudy(identifier);
1096    
1097     // SpecialCutFlow(identifier);
1098     // draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",17000,30,200,"Z p_{T}","Official/"+identifier+"/Zpt__"+identifier,1);
1099     // draw_kin_variable(identifier+"_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",40,0,200,"Leading Jet Pt (B)","Official/"+identifier+"/LeadingJetPt__"+identifier,1);
1100     // draw_kin_variable(identifier+"_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",40,0,200,"Sub-Leading Jet Pt [GeV]","Official/"+identifier+"/Jet2Pt__"+identifier,1);
1101     // draw_kin_variable(identifier+"_pfJetGoodPt[1]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,2,"p_{T}^{J2} / p_{T}^{Z}","Official/"+identifier+"/SubLeadingJetPt_Over_ZPt__"+identifier,1);
1102     draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<50"),20,0,2,"#alpha","Official/"+identifier+"/alpha_pt_30_to_50__"+identifier,1);
1103     draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>50&&pt<75"),20,0,2,"#alpha","Official/"+identifier+"/alpha_50_to_75__"+identifier,1);
1104     draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>75&&pt<125"),20,0,2,"#alpha","Official/"+identifier+"/alpha_pt_75_to_125__"+identifier,1);
1105     draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>125&&pt<1000"),20,0,2,"#alpha","Official/"+identifier+"/alpha_pt_125_to_1000__"+identifier,1);
1106     draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut,20,0,2,"#alpha","Official/"+identifier+"/alpha_full__"+identifier,1);
1107     // draw_kin_variable(identifier+"_pfBJetDphiZ[0]",ZplusBsel&&LeadingB&&"pt>10&&pt<1000&&pfBJetDphiZ[0]>0",30,0,3.2,"#delta#phi (Z,b lead)","Official/"+identifier+"/DeltaPhiZBlead__"+identifier,0);
1108 buchmann 1.7
1109 buchmann 1.9 // draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.3").c_str()),20,0,2,"Z+b MPF","Official/"+identifier+"/mpf_alpha_smaller_0p3___"+identifier,0,"#alpha<0.3, 10 GeV < p_{T}^{Z} < 1000 GeV");
1110     // draw_kin_variable(identifier+"_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.3").c_str()),20,0,2,"p_{T} b-jet / p_{T} Z","Official/"+identifier+"/ptb_over_ptz___alpha_smaller_0p3______"+identifier,0,"#alpha<0.3, 10 GeV < p_{T}^{Z} < 1000 GeV");
1111 buchmann 1.7
1112    
1113 buchmann 1.9 // */
1114 buchmann 1.7
1115     }
1116    
1117 buchmann 1.14 void GetNumberEventsInsideOutsideAlphaWindow(TCut specialcut, string title) {
1118 buchmann 1.13
1119 buchmann 1.14 // TCut cut = ZplusBsel&&LeadingB&&PhiZcut&&specialcut&&TCut("Zb3010_alpha<0.3");
1120     TCut cut = specialcut;
1121     TCut in = EtaB;
1122     TCut out = TCut("abs(Zb3010_pfJetGoodEta[0])>1.3");
1123 buchmann 1.13
1124     TH1F *data_IN = allsamples.Draw("data_IN", "mll",1,0,150, "nothing [GeV]", "events", cut&&in, data,luminosity);
1125     TH1F *data_OUT = allsamples.Draw("data_OUT", "mll",1,0,150, "nothing [GeV]", "events", cut&&out,data,luminosity);
1126     TH1F *mc_IN = allsamples.Draw("mc_IN", "mll",1,0,150, "nothing [GeV]", "events", cut&&in, mc ,luminosity);
1127     TH1F *mc_OUT = allsamples.Draw("mc_OUT", "mll",1,0,150, "nothing [GeV]", "events", cut&&out,mc ,luminosity);
1128    
1129 buchmann 1.14 dout << title << " : " << endl;
1130 buchmann 1.13 dout << " Data: In " << data_IN->Integral() << " vs out " << data_OUT->Integral() << endl;
1131     dout << " MC : In " << mc_IN->Integral() << " vs out " << mc_OUT->Integral() << endl;
1132 buchmann 1.14 dout << " ratio in : " << data_IN->Integral()/mc_IN->Integral() << " +/- " << (data_IN->Integral()/mc_IN->Integral()) * sqrt(1.0/mc_IN->Integral()+ 1.0/data_IN->Integral()) << endl;
1133     dout << " ratio out : " << data_OUT->Integral()/mc_OUT->Integral() << " +/- " << (data_OUT->Integral()/mc_OUT->Integral()) * sqrt(1.0/mc_OUT->Integral()+ 1.0/data_OUT->Integral()) << endl;
1134    
1135     delete data_IN;
1136     delete data_OUT;
1137     delete mc_IN;
1138     delete mc_OUT;
1139 buchmann 1.13 }
1140    
1141 buchmann 1.15 void TEST_GetNumberEventsInsideOutsideAlphaWindow() {
1142 buchmann 1.14 TCanvas *ca = new TCanvas("ca","ca");
1143     // GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==0"), "ee");
1144     // GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==1"), "mm");
1145     // GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2"), "ee/mm");
1146    
1147     cout << "BASE SELECTION" << endl;
1148     GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==0"), "ee");
1149     GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==1"), "mm");
1150     GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2"), "ee/mm");
1151    
1152     cout << "BASE SELECTION: Z+b" << endl;
1153     GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1154     GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1155     GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2"), "ee/mm");
1156    
1157     cout << "Leading B" << endl;
1158     GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1159     GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1160     GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2"), "ee/mm");
1161    
1162     cout << "PhiZcut" << endl;
1163     GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1164     GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1165     GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2"), "ee/mm");
1166    
1167     cout << "alpha cut" << endl;
1168     GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0")&&TCut("Zb3010_alpha<0.3"), "ee");
1169     GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1")&&TCut("Zb3010_alpha<0.3"), "mm");
1170     GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2")&&TCut("Zb3010_alpha<0.3"), "ee/mm");
1171    
1172     delete ca;
1173    
1174    
1175     }
1176    
1177    
1178    
1179 buchmann 1.13
1180 buchmann 1.15 void do_basic_ZB_analysis(bool DoCompleteAnalysis) {
1181 buchmann 1.12
1182     //https://twiki.cern.ch/twiki/bin/viewauth/CMS/BtagPOG#2011_Data_and_MC
1183     //https://twiki.cern.ch/twiki/pub/CMS/BtagPOG/SFb-ttbar_payload.txt
1184 buchmann 1.9
1185 buchmann 1.14 cout << "The lepton requirement is " << (const char*) leptoncut << endl;
1186 buchmann 1.12 bool doquick=true;
1187 buchmann 1.9
1188 buchmann 1.2 TCanvas *zbcanvas = new TCanvas("zbcanvas","zbcanvas");
1189 buchmann 1.1
1190 buchmann 1.9 // compare_selection("Zb3010_p5Clean");
1191     // compare_selection("Zb3010CHS");
1192     // compare_selection("Zb3010");
1193     // compare_selection("Zb301030");
1194     // compare_selection("Zb1530");
1195     // compare_selection("Zb301010");
1196     // compare_selection("Zb3010_SecEta3");
1197     // compare_selection("Zb3010_SecEta5");
1198     // compare_selection("Zb1510");
1199     // compare_selection("Zb301010");
1200     // compare_selection("ZbMikko");
1201     // compare_selection("Zb3010");
1202     // compare_selection("Zb40");
1203 buchmann 1.8
1204 buchmann 1.11 if(!doquick) {
1205     print_all_b_yields();
1206     draw_mpf_vars();
1207     // DrawEvilCutFlow();
1208    
1209     draw_Zb_kin_vars();
1210    
1211     }
1212 buchmann 1.7
1213 buchmann 1.14 // GetNumberEventsInsideOutsideAlphaWindow();
1214 buchmann 1.15
1215     dout <<"THIS IS WHERE IT GETS REALLY INTERESTING!!!!" << endl;
1216     if(DoCompleteAnalysis) {
1217     /* need to run the following scenarios:
1218     * - normal
1219     * - inclusive
1220     * - medium WP
1221     * - alpha up
1222     * - alpha down
1223     */
1224    
1225     ZbResultContainer inclusive = new_data_mc_agreement_2d("Zb3010_alpha","inclusive",TCut("Zb3010_bTagProbCSVBP[0]>-10000"),TCut("1.0"));
1226     ZbResultContainer Reference = new_data_mc_agreement_2d("Zb3010_alpha","reference",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgt"));
1227     ZbResultContainer effmistagUp = new_data_mc_agreement_2d("Zb3010_alpha","effmistagUp",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtUp"));
1228     ZbResultContainer effmistagDown = new_data_mc_agreement_2d("Zb3010_alpha","effmistagDown",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtDown"));
1229     ZbResultContainer medium = new_data_mc_agreement_2d("Zb3010_alpha","medium",TCut("Zb3010_bTagProbCSVBP[0]>0.679"),TCut("ZbCHS3010_BTagWgtMedium"));
1230     ZbResultContainer AlphaUp = new_data_mc_agreement_2d("Zb3010_alphaUP","AlphaUp",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgt"));
1231     ZbResultContainer AlphaDown = new_data_mc_agreement_2d("Zb3010_alphaDOWN","AlphaDown",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgt"));
1232    
1233     LeadingB=TCut("Zb3010_bTagProbCSVBP[0]>0.898");
1234     } else {
1235     ZbResultContainer Reference = new_data_mc_agreement_2d("Zb3010_alpha","reference",TCut("Zb3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgt"));
1236     }
1237    
1238 buchmann 1.1
1239 buchmann 1.2 delete zbcanvas;
1240 buchmann 1.1 }
1241 buchmann 1.7
1242    
1243    
1244    
1245     //*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/* DELETE EVERYTHING BELOW THIS LINE /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
1246    
1247    
1248     const char* concatenate(string bla1, string bla2) {
1249     stringstream bla;
1250     bla << bla1 << bla2;
1251     return bla.str().c_str();
1252     }
1253    
1254     void SpecialCutFlow(string identifier) {
1255     stringstream MegaCut;
1256    
1257    
1258    
1259     MegaCut << " 1*(" << (const char*) (ZplusBsel&&TCut(concatenate(identifier,"_pfJetGoodNumBtag>0"))) << ")";
1260     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB) << ")";
1261     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB) << ")";
1262     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut) << ")";
1263 buchmann 1.9 MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000").c_str())) << ")";
1264     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.3").c_str())) << ")";
1265     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.2").c_str())) << ")";
1266     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.15").c_str())) << ")";
1267     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.1").c_str())) << ")";
1268     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.05").c_str())) << ")";
1269 buchmann 1.7 MegaCut << " ";
1270    
1271     TCanvas *can = new TCanvas("can","can");
1272     can->SetLogy(1);
1273 buchmann 1.9 TCut basecut(ZplusBsel);
1274 buchmann 1.7
1275    
1276    
1277     TLegend *leg = allsamples.allbglegend();
1278    
1279    
1280     TH1F *data = allsamples.Draw("data", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,data,luminosity);
1281     THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
1282    
1283     float runningsum=0;
1284     for(int i=data->GetNbinsX();i>0;i--) {
1285     runningsum+=data->GetBinContent(i);
1286     data->SetBinContent(i,runningsum);
1287     }
1288    
1289     float minimum=data->GetMaximum();
1290    
1291     TH1F* h;
1292     TIter nextOF(mcm.GetHists());
1293    
1294     float mcsum=0;
1295     float zb=0;
1296     while ( h = (TH1F*)nextOF() ) {
1297     float runningsum=0;
1298     minimum=data->GetMaximum();//done deliberately at each step so get the minimum of the last (!) contribution
1299     for(int i=h->GetNbinsX();i>0;i--) {
1300     runningsum+=h->GetBinContent(i);
1301     h->SetBinContent(i,runningsum);
1302     if(runningsum<minimum&&runningsum>0) minimum=runningsum;
1303     }
1304     if(Contains(h->GetName(),"Z_b_")) {
1305     zb=h->GetBinContent(h->GetNbinsX());
1306     }
1307     mcsum+=h->GetBinContent(h->GetNbinsX());
1308     }
1309    
1310    
1311    
1312     data->SetMinimum(0.05*minimum);
1313     can->cd(1)->SetBottomMargin(0.25);
1314     data->GetXaxis()->SetBinLabel(1,"Pre-selection");
1315     data->GetXaxis()->SetBinLabel(2,"Contains b jet");
1316     data->GetXaxis()->SetBinLabel(3,"Leading jet is b-jet");
1317     data->GetXaxis()->SetBinLabel(4,"|#eta_{b}|<1.3 ");
1318     data->GetXaxis()->SetBinLabel(5,"|#delta#phi(b,Z)|>2.7");
1319     data->GetXaxis()->SetBinLabel(6,"10<p_{T}^{Z}<1000 GeV");
1320     data->GetXaxis()->SetBinLabel(7,"#alpha < 0.3");
1321     data->GetXaxis()->SetBinLabel(8,"#alpha < 0.2");
1322     data->GetXaxis()->SetBinLabel(9,"#alpha < 0.15");
1323     data->GetXaxis()->SetBinLabel(10,"#alpha < 0.1");
1324     data->GetXaxis()->SetBinLabel(11,"#alpha < 0.05");
1325     data->GetXaxis()->LabelsOption("v");
1326    
1327     stringstream purityinfo;
1328     purityinfo << "Purity: " << std::setprecision(3) << 100*zb/mcsum << " %";
1329     stringstream neventsinfo;
1330     neventsinfo << "Nevents: " << data->GetBinContent(data->GetNbinsX());
1331     TH1F *crap = new TH1F("crap","",1,0,1);
1332     crap->SetLineColor(kWhite);
1333     leg->AddEntry(crap,purityinfo.str().c_str(),"l");
1334     leg->AddEntry(crap,neventsinfo.str().c_str(),"l");
1335    
1336     data->Draw("e1");
1337     mcm.Draw("same");
1338     data->Draw("same,e1,axis");
1339     data->Draw("same,e1");
1340     // data->Draw("same,TEXT");
1341    
1342     leg->Draw();
1343    
1344     DrawPrelim();
1345    
1346     CompleteSave(can,"CutFlow___"+identifier);
1347    
1348     delete crap;
1349     delete data;
1350    
1351    
1352     }