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