StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
plotMCDataComp_2GeV.C
1 // This macro takes as its input, the output of the
2 // the makeBkgdFiles.C macro and the MC files for W+ and W-
3 // and compares/plots the data and MC expectations
4 
5 void plotMCDataComp_2GeV() {
6 
7  TFile *f1 = new TFile("bkgd_histos_pos_output4.root");
8  TFile *f2 = new TFile("bkgd_histos_neg_output4.root");
9  TFile *f3 = new TFile("rcf10010_atMCscale.wana.hist.root");
10  TFile *f4 = new TFile("rcf10011_atMCscale.wana.hist.root");
11 
12  TH1F *pos_data = (TH1F*)f1->Get("signal_final3");
13  TH1F *neg_data = (TH1F*)f2->Get("signal_final3");
14  TH1F *pos_mc = (TH1F*)f3->Get("muclustPtBal");
15  TH1F *neg_mc = (TH1F*)f4->Get("muclustPtBal");
16  TH1F *pos_data_repack = new TH1F("pos_data_repack","pos_data_repack",49,1.,99.);
17  TH1F *neg_data_repack = new TH1F("neg_data_repack","neg_data_repack",49,1.,99.);
18 
19  // repack the data to 2 GeV wide bins
20  for (int i=1; i<=49; i++) {
21  pos_data_repack->SetBinContent(i,pos_data->GetBinContent(i));
22  neg_data_repack->SetBinContent(i,neg_data->GetBinContent(i));
23  }
24 
25  TH1F *pos_mc_repack = new TH1F("pos_mc_repack","pos_mc_repack",49,1.,99.);
26  TH1F *neg_mc_repack = new TH1F("neg_mc_repack","neg_mc_repack",49,1.,99.);
27  for (int i=1; i<=49; i++) {
28  pos_mc_repack->SetBinContent(i,pos_mc->GetBinContent(2*i-1)+
29  pos_mc->GetBinContent(2*i));
30  neg_mc_repack->SetBinContent(i,neg_mc->GetBinContent(2*i-1)+
31  neg_mc->GetBinContent(2*i));
32  }
33 
34  // Find MC normalization factors to scale by
35  // and scale the MC
36  float pos_data_cnt = 0., pos_mc_cnt = 0.;
37  float neg_data_cnt = 0., neg_mc_cnt = 0.;
38  for (int i=13; i<=49; i++) {
39  cout << i << " " << pos_data_repack->GetBinContent(i) << endl;
40  pos_data_cnt += pos_data_repack->GetBinContent(i);
41  neg_data_cnt += neg_data_repack->GetBinContent(i);
42  pos_mc_cnt += pos_mc_repack->GetBinContent(i);
43  neg_mc_cnt += neg_mc_repack->GetBinContent(i);
44  }
45  pos_mc_repack->Scale(pos_data_cnt/pos_mc_cnt);
46  neg_mc_repack->Scale(neg_data_cnt/neg_mc_cnt);
47  cout << "Pos Data = " << pos_data_cnt << " and MC = " << pos_mc_cnt << endl;
48  cout << "Neg Data = " << neg_data_cnt << " and MC = " << neg_mc_cnt << endl;
49  cout << "Norm ratio = " << pos_data_cnt/pos_mc_cnt << endl;
50 
51  // *******************************************************
52  // Grab the statistical and QCD bkgd shape uncertainties
53  // *******************************************************
54 
55  // Get all the necessary histograms for the calculation
56  // of all the uncertainties
57  TH1F *pos_raw_stat_err2 = (TH1F*)f1->Get("raw_stat_err2");
58  TH1F *pos_QCD_stat_err2 = (TH1F*)f1->Get("QCD_stat_err2");
59  TH1F *pos_eemc_stat_err2 = (TH1F*)f1->Get("eemc_stat_err2");
60  TH1F *pos_tau_stat_err2 = (TH1F*)f1->Get("tau_stat_err2");
61  TH1F *pos_QCD_syst_low_err = (TH1F*)f1->Get("QCD_syst_low_err");
62  TH1F *pos_QCD_syst_high_err = (TH1F*)f1->Get("QCD_syst_high_err");
63 
64  TH1F *neg_raw_stat_err2 = (TH1F*)f2->Get("raw_stat_err2");
65  TH1F *neg_QCD_stat_err2 = (TH1F*)f2->Get("QCD_stat_err2");
66  TH1F *neg_eemc_stat_err2 = (TH1F*)f2->Get("eemc_stat_err2");
67  TH1F *neg_tau_stat_err2 = (TH1F*)f2->Get("tau_stat_err2");
68  TH1F *neg_QCD_syst_low_err = (TH1F*)f2->Get("QCD_syst_low_err");
69  TH1F *neg_QCD_syst_high_err = (TH1F*)f2->Get("QCD_syst_high_err");
70 
71  float pos_high[49], pos_low[49];
72  float neg_high[49], neg_low[49];
73  float bin_cent[49], pos_bin_val[49], neg_bin_val[49];
74  float ex[49];
75  float err;
76 
77  // calculate the statistical uncertainty and the background
78  // shape systematic uncertainty
79  for (int i=1; i<=49; i++) {
80  ex[i-1] = 2.0;
81 
82  err = 0.;
83  err += pos_raw_stat_err2->GetBinContent(i);
84  err += pos_QCD_stat_err2->GetBinContent(i);
85  err += pos_eemc_stat_err2->GetBinContent(i);
86  err += pos_tau_stat_err2->GetBinContent(i);
87  pos_data_repack->SetBinError(i,sqrt(err));
88 
89  err = 0.;
90  err += neg_raw_stat_err2->GetBinContent(i);
91  err += neg_QCD_stat_err2->GetBinContent(i);
92  err += neg_eemc_stat_err2->GetBinContent(i);
93  err += neg_tau_stat_err2->GetBinContent(i);
94  neg_data_repack->SetBinError(i,sqrt(err));
95 
96  bin_cent[i-1] = pos_data_repack->GetBinCenter(i);
97  pos_bin_val[i-1] = pos_data_repack->GetBinContent(i);
98  neg_bin_val[i-1] = neg_data_repack->GetBinContent(i);
99 
100  // The QCD background shape systematic errors are
101  // put in their arrays
102  pos_high[i-1] = pos_QCD_syst_high_err->GetBinContent(i);
103  neg_high[i-1] = neg_QCD_syst_high_err->GetBinContent(i);
104  pos_low[i-1] = pos_QCD_syst_low_err->GetBinContent(i);
105  neg_low[i-1] = neg_QCD_syst_low_err->GetBinContent(i);
106  }
107 
108  // *******************************************************
109  // Calculate an energy scale systematic
110  // *******************************************************
111 
112  TFile *f5 = new TFile("bkgd_histos_pos_energy_up.root");
113  TFile *f6 = new TFile("bkgd_histos_pos_energy_down.root");
114  TFile *f7 = new TFile("bkgd_histos_neg_energy_up.root");
115  TFile *f8 = new TFile("bkgd_histos_neg_energy_down.root");
116 
117  TH1F *pos_data_en_up = (TH1F*)f5->Get("signal_final3");
118  TH1F *pos_data_en_down = (TH1F*)f6->Get("signal_final3");
119  TH1F *neg_data_en_up = (TH1F*)f7->Get("signal_final3");
120  TH1F *neg_data_en_down = (TH1F*)f8->Get("signal_final3");
121 
122  TH1F *pos_data_en_up_repack = new TH1F("pos_data_en_up_repack","pos_data_en_up_repack",49,1.,99.);
123  TH1F *pos_data_en_down_repack = new TH1F("pos_data_en_down_repack","pos_data_en_down_repack",49,1.,99.);
124  TH1F *neg_data_en_up_repack = new TH1F("neg_data_en_up_repack","neg_data_en_up_repack",49,1.,99.);
125  TH1F *neg_data_en_down_repack = new TH1F("neg_data_en_down_repack","neg_data_en_down_repack",49,1.,99.);
126 
127  for (int i=1; i<=49; i++) {
128  pos_data_en_up_repack->SetBinContent(i,pos_data_en_up->GetBinContent(i));
129  pos_data_en_down_repack->SetBinContent(i,pos_data_en_down->GetBinContent(i));
130  neg_data_en_up_repack->SetBinContent(i,neg_data_en_up->GetBinContent(i));
131  neg_data_en_down_repack->SetBinContent(i,neg_data_en_down->GetBinContent(i));
132  }
133 
134  for (int i=1; i<=49; i++) {
135  cout << "W+ -> " << pos_data_en_up_repack->GetBinContent(i) << " " << pos_data_repack->GetBinContent(i) << " " << pos_data_en_down_repack->GetBinContent(i) << endl;
136  }
137  for (int i=1; i<=49; i++) {
138  cout << "W- -> " << neg_data_en_up_repack->GetBinContent(i) << " " << neg_data_repack->GetBinContent(i) << " " << neg_data_en_down_repack->GetBinContent(i) << endl;
139  }
140 
141  // Now loop over all the bins of the three histograms (energy low,
142  // normal, energy high) and in each bin find the lowest and the
143  // highest and then call that the systematic (normally it will be
144  // the two energy scaled data sets) ... and do it for both charges
145  float pos_data_ensys_low[49];
146  float pos_data_ensys_high[49];
147  float neg_data_ensys_low[49];
148  float neg_data_ensys_high[49];
149  for (int i=1; i<=49; i++) {
150 
151  // for positive charge
152  float low = 10000.;
153  if (pos_data_en_up_repack->GetBinContent(i) < low) {
154  low = pos_data_en_up_repack->GetBinContent(i);
155  }
156  if (pos_data_repack->GetBinContent(i) < low) {
157  low = pos_data_repack->GetBinContent(i);
158  }
159  if (pos_data_en_down_repack->GetBinContent(i) < low) {
160  low = pos_data_en_down_repack->GetBinContent(i);
161  }
162  pos_data_ensys_low[i-1] = fabs(low-pos_data_repack->GetBinContent(i));
163 
164  float high = 0.;
165  if (pos_data_en_up_repack->GetBinContent(i) > high) {
166  high = pos_data_en_up_repack->GetBinContent(i);
167  }
168  if (pos_data_repack->GetBinContent(i) > high) {
169  high = pos_data_repack->GetBinContent(i);
170  }
171  if (pos_data_en_down_repack->GetBinContent(i) > high) {
172  high = pos_data_en_down_repack->GetBinContent(i);
173  }
174  pos_data_ensys_high[i-1] = fabs(high-pos_data_repack->GetBinContent(i));
175 
176  // for negative charge
177  low = 10000.;
178  if (neg_data_en_up_repack->GetBinContent(i) < low) {
179  low = neg_data_en_up_repack->GetBinContent(i); }
180  if (neg_data_repack->GetBinContent(i) < low) {
181  low = neg_data_repack->GetBinContent(i);
182  }
183  if (neg_data_en_down_repack->GetBinContent(i) < low) {
184  low = neg_data_en_down_repack->GetBinContent(i);
185  }
186  neg_data_ensys_low[i-1] = fabs(low-neg_data_repack->GetBinContent(i));
187  high = 0.;
188  if (neg_data_en_up_repack->GetBinContent(i) > high) { high = neg_data_en_up_repack->GetBinContent(i);
189  }
190  if (neg_data_repack->GetBinContent(i) > high) {
191  high = neg_data_repack->GetBinContent(i);
192  }
193  if (neg_data_en_down_repack->GetBinContent(i) > high) {
194  high = neg_data_en_down_repack->GetBinContent(i);
195  }
196  neg_data_ensys_high[i-1] = fabs(high-neg_data_repack->GetBinContent(i));
197 
198  cout << i << " " << pos_data_ensys_low[i-1] << " " << pos_data_ensys_high[i-1] <<
199 " " << neg_data_ensys_low[i-1] << " " << neg_data_ensys_high[i-1] << endl;
200 
201  // convention is switched (add to the QCD background systematic)
202  // assume fully correlated (probably not true)
203  pos_high[i-1] += pos_data_ensys_low[i-1];
204  pos_low[i-1] += pos_data_ensys_high[i-1];
205  neg_high[i-1] += neg_data_ensys_low[i-1];
206  neg_low[i-1] += neg_data_ensys_high[i-1];
207  }
208 
209  TGraphAsymmErrors *pos_syst_err_box = new TGraphAsymmErrors(49,bin_cent,pos_bin_val,ex,ex,pos_high,pos_low);
210  TGraphAsymmErrors *neg_syst_err_box = new TGraphAsymmErrors(49,bin_cent,neg_bin_val,ex,ex,neg_high,neg_low);
211 
212  // ******************************************************
213  // Calculate the MC to Data ratios now (even though they
214  // are not plotted)
215  // ******************************************************
216 
217  TH1F *pos_ratio = pos_data_repack->Clone();
218  TH1F *neg_ratio = neg_data_repack->Clone();
219  float pos_ratio_syst_high[49], pos_ratio_syst_low[49];
220  float neg_ratio_syst_high[49], neg_ratio_syst_low[49];
221  float pos_ratio_bin_val[49], neg_ratio_bin_val[49];
222 
223  for (int i=1; i<=49; i++) {
224  if ((pos_mc_repack->GetBinContent(i) > 0) && (pos_data_repack->GetBinContent(i) > 0)) {
225  pos_ratio->SetBinContent(i,pos_data_repack->GetBinContent(i)/
226  pos_mc_repack->GetBinContent(i));
227  pos_ratio->SetBinError(i,pos_data_repack->GetBinError(i)/pos_mc_repack->GetBinContent(i));
228  pos_ratio_bin_val[i-1] = pos_data_repack->GetBinContent(i)/
229  pos_mc_repack->GetBinContent(i);
230  pos_ratio_syst_high[i-1] = (pos_data_repack->GetBinContent(i)+pos_high[i-1])/pos_mc_repack->GetBinContent(i);
231  pos_ratio_syst_low[i-1] = (pos_data_repack->GetBinContent(i)-pos_low[i-1])/pos_mc_repack->GetBinContent(i);
232  } else {
233  pos_ratio->SetBinContent(i,0.);
234  pos_ratio->SetBinError(i,0.);
235  pos_ratio_syst_high[i-1] = 0.;
236  pos_ratio_syst_low[i-1] = 0.;
237  pos_ratio_bin_val[i-1] = 0.;
238  }
239 
240  if ((neg_mc_repack->GetBinContent(i) > 0) && (neg_data_repack->GetBinContent(i) > 0)) {
241  neg_ratio->SetBinContent(i,neg_data_repack->GetBinContent(i)/
242  neg_mc_repack->GetBinContent(i));
243  neg_ratio->SetBinError(i,neg_data_repack->GetBinError(i)/neg_mc_repack->GetBinContent(i));
244  neg_ratio_bin_val[i-1] = neg_data_repack->GetBinContent(i)/
245  neg_mc_repack->GetBinContent(i);
246  neg_ratio_syst_high[i-1] = (neg_data_repack->GetBinContent(i)+neg_high[i-1])/neg_mc_repack->GetBinContent(i);
247  neg_ratio_syst_low[i-1] = (neg_data_repack->GetBinContent(i)-neg_low[i-1])/neg_mc_repack->GetBinContent(i);
248  } else {
249  neg_ratio->SetBinContent(i,0.);
250  neg_ratio->SetBinError(i,0.);
251  neg_ratio_syst_high[i-1] = 0.;
252  neg_ratio_syst_low[i-1] = 0.;
253  neg_ratio_bin_val[i-1] = 0.;
254  }
255  }
256 
257  for (int i=0; i<49; i++) {
258  pos_ratio_syst_high[i] -= pos_ratio_bin_val[i];
259  pos_ratio_syst_low[i] = pos_ratio_bin_val[i]-pos_ratio_syst_low[i];
260  neg_ratio_syst_high[i] -= neg_ratio_bin_val[i];
261  neg_ratio_syst_low[i] = neg_ratio_bin_val[i]-neg_ratio_syst_low[i];
262  }
263  TGraphAsymmErrors *pos_ratio_syst_err_box = new TGraphAsymmErrors(49,bin_cent,pos_ratio_bin_val,ex,ex,pos_ratio_syst_high,pos_ratio_syst_low);
264  TGraphAsymmErrors *neg_ratio_syst_err_box = new TGraphAsymmErrors(49,bin_cent,neg_ratio_bin_val,ex,ex,neg_ratio_syst_high,neg_ratio_syst_low);
265 
266  // ******************************************************
267  // Plots versus ET (no longer the ratios)
268  // ******************************************************
269  gStyle->SetTitleBorderSize(0);
270  gStyle->SetOptDate(0);
271  TCanvas *can1 = new TCanvas("can1","can1",0,0,500,500);
272 
273  gPad->SetBorderMode(0);
274  gPad->SetGridx(0);
275  gPad->SetGridy(0);
276  pos_mc_repack->GetYaxis()->SetRangeUser(-10.,100.);
277  pos_mc_repack->GetXaxis()->SetRangeUser(0.,70.);
278  pos_mc_repack->SetLineColor(2);
279  pos_mc_repack->SetStats(kFALSE);
280  pos_mc_repack->GetYaxis()->SetTitle("Counts");
281  pos_mc_repack->GetYaxis()->SetTitleOffset(1.2);
282  pos_mc_repack->GetXaxis()->SetTitle("EMC Cluster E_{T} (GeV)");
283  pos_mc_repack->SetTitle("");
284  pos_mc_repack->SetLineWidth(3.*pos_mc_repack->GetLineWidth());
285  pos_data_repack->SetLineWidth(2.*pos_data_repack->GetLineWidth());
286  pos_mc_repack->Draw();
287  pos_syst_err_box->SetFillColor(14);
288  pos_syst_err_box->Draw("e2same");
289  pos_mc_repack->Draw("same");
290  pos_data_repack->Draw("same");
291  TLegend *leg1 = new TLegend(0.5,0.7,0.88,0.88);
292  leg1->SetFillColor(0);
293  leg1->SetBorderSize(0);
294  leg1->AddEntry(pos_mc_repack,"MC Norm. to Signal","l");
295  leg1->AddEntry(pos_data_repack,"Signal","le");
296  leg1->AddEntry(pos_syst_err_box,"Syst. Unc.","f");
297  leg1->Draw("same");
298  TLatex *lat1 = new TLatex();
299  lat1->SetNDC();
300  lat1->DrawLatex(0.14,0.84,"p+p#rightarrow W^{+}#rightarrow e^{+}+#nu");
301  lat1->DrawLatex(0.14,0.74,"Run9 STAR");
302  lat1->DrawLatex(0.14,0.69,"Preliminary");
303  lat1->DrawLatex(0.14,0.64,"#sqrt{s}=500 GeV");
304  can1->Print("data_mc_ET_comp_output4_W+_2GeV.eps");
305 
306  TCanvas *can2 = new TCanvas("can2","can2",0,0,500,500);
307  gPad->SetBorderMode(0);
308  gPad->SetGridx(0);
309  gPad->SetGridy(0);
310  neg_mc_repack->GetYaxis()->SetRangeUser(-5.,50.);
311  neg_mc_repack->GetXaxis()->SetRangeUser(0.,70.);
312  neg_mc_repack->SetLineColor(2);
313  neg_mc_repack->SetStats(kFALSE);
314  neg_mc_repack->GetYaxis()->SetTitle("Counts");
315  neg_mc_repack->GetYaxis()->SetTitleOffset(1.2);
316  neg_mc_repack->GetXaxis()->SetTitle("EMC Cluster E_{T} (GeV)");
317  neg_mc_repack->SetTitle("");
318  neg_mc_repack->SetLineWidth(3.*neg_mc_repack->GetLineWidth());
319  neg_data_repack->SetLineWidth(2.*neg_data_repack->GetLineWidth());
320  neg_mc_repack->Draw();
321  neg_syst_err_box->SetFillColor(14);
322  neg_syst_err_box->Draw("e2same");
323  neg_mc_repack->Draw("same");
324  neg_data_repack->Draw("same");
325  leg1->Draw("same");
326  lat1->DrawLatex(0.14,0.84,"p+p#rightarrow W^{-}#rightarrow e^{-}+#bar{#nu}");
327  lat1->DrawLatex(0.14,0.74,"Run9 STAR");
328  lat1->DrawLatex(0.14,0.69,"Preliminary");
329  lat1->DrawLatex(0.14,0.64,"#sqrt{s}=500 GeV");
330  can2->Print("data_mc_ET_comp_output4_W-_2GeV.eps");
331 
332 } // end of macro
Generated by   doxygen 1.8.5