StRoot: StSpinPool/StWalgoB2009/macros/plotMCDataComp

00001 // This macro takes as its input, the output of the
00002 // the makeBkgdFiles.C macro and the MC files for W+ and W-
00003 // and compares/plots the data and MC expectations
00004 
00005 void plotMCDataComp_4GeV() {
00006 
00007   TFile *f1 = new TFile("bkgd_histos_pos_output4.root");
00008   TFile *f2 = new TFile("bkgd_histos_neg_output4.root");
00009   TFile *f3 = new TFile("rcf10010_atMCscale.wana.hist.root");
00010   TFile *f4 = new TFile("rcf10011_atMCscale.wana.hist.root");
00011 
00012   TH1F *pos_data = (TH1F*)f1->Get("signal_final3");
00013   TH1F *neg_data = (TH1F*)f2->Get("signal_final3");
00014   TH1F *pos_mc   = (TH1F*)f3->Get("muclustPtBal");
00015   TH1F *neg_mc   = (TH1F*)f4->Get("muclustPtBal");
00016   TH1F *pos_data_repack = new TH1F("pos_data_repack","pos_data_repack",24,3.,99.);
00017   TH1F *neg_data_repack = new TH1F("neg_data_repack","neg_data_repack",24,3.,99.);
00018 
00019   // repack the data and MC into 4 GeV wide bins
00020   for (int i=1; i<=24; i++) {
00021     pos_data_repack->SetBinContent(i,pos_data->GetBinContent(2*i)+
00022                                      pos_data->GetBinContent(2*i+1));
00023     neg_data_repack->SetBinContent(i,neg_data->GetBinContent(2*i)+
00024                                      neg_data->GetBinContent(2*i+1));
00025   }
00026   
00027   TH1F *pos_mc_repack = new TH1F("pos_mc_repack","pos_mc_repack",24,3.,99.);
00028   TH1F *neg_mc_repack = new TH1F("neg_mc_repack","neg_mc_repack",24,3.,99.);
00029   for (int i=1; i<=24; i++) {
00030     cout << "Summing bins " << pos_mc->GetBinCenter(4*i) << ","
00031          << pos_mc->GetBinCenter(4*i+1) << ","   
00032          << pos_mc->GetBinCenter(4*i+2) << ","
00033          << pos_mc->GetBinCenter(4*i+3) << " into bin " << pos_mc_repack->GetBinCenter(i) << endl;
00034     pos_mc_repack->SetBinContent(i,pos_mc->GetBinContent(4*i)+
00035                                    pos_mc->GetBinContent(4*i+1)+
00036                                    pos_mc->GetBinContent(4*i+2)+
00037                                    pos_mc->GetBinContent(4*i+3));
00038     neg_mc_repack->SetBinContent(i,neg_mc->GetBinContent(4*i)+
00039                                    neg_mc->GetBinContent(4*i+1)+
00040                                    neg_mc->GetBinContent(4*i+2)+
00041                                    neg_mc->GetBinContent(4*i+3));
00042   } 
00043 
00044   // Find MC normalization factors to scale by
00045   // and scale the MC
00046   float pos_data_cnt = 0., pos_mc_cnt = 0.;
00047   float neg_data_cnt = 0., neg_mc_cnt = 0.;
00048   for (int i=4; i<=24; i++) {
00049     cout << i << " " << pos_data_repack->GetBinContent(i) << endl;
00050     pos_data_cnt += pos_data_repack->GetBinContent(i);
00051     neg_data_cnt += neg_data_repack->GetBinContent(i);
00052     pos_mc_cnt   += pos_mc_repack->GetBinContent(i);
00053     neg_mc_cnt   += neg_mc_repack->GetBinContent(i);
00054   }
00055   pos_mc_repack->Scale(pos_data_cnt/pos_mc_cnt);
00056   neg_mc_repack->Scale(neg_data_cnt/neg_mc_cnt);
00057   cout << "Pos Data = " << pos_data_cnt << " and MC = " << pos_mc_cnt << endl;
00058   cout << "Neg Data = " << neg_data_cnt << " and MC = " << neg_mc_cnt << endl;
00059   cout << "Norm ratio = " << pos_data_cnt/pos_mc_cnt << endl;
00060 
00061   // *******************************************************
00062   // Grab the statistical and QCD bkgd shape uncertainties
00063   // *******************************************************
00064  
00065   // Get all the necessary histograms for the calculation
00066   // of all the uncertainties
00067   TH1F *pos_raw_stat_err2 = (TH1F*)f1->Get("raw_stat_err2");
00068   TH1F *pos_QCD_stat_err2 = (TH1F*)f1->Get("QCD_stat_err2");
00069   TH1F *pos_eemc_stat_err2 = (TH1F*)f1->Get("eemc_stat_err2");
00070   TH1F *pos_tau_stat_err2 = (TH1F*)f1->Get("tau_stat_err2");
00071   TH1F *pos_QCD_syst_low_err = (TH1F*)f1->Get("QCD_syst_low_err");
00072   TH1F *pos_QCD_syst_high_err = (TH1F*)f1->Get("QCD_syst_high_err");
00073 
00074   TH1F *neg_raw_stat_err2 = (TH1F*)f2->Get("raw_stat_err2");
00075   TH1F *neg_QCD_stat_err2 = (TH1F*)f2->Get("QCD_stat_err2");
00076   TH1F *neg_eemc_stat_err2 = (TH1F*)f2->Get("eemc_stat_err2");
00077   TH1F *neg_tau_stat_err2 = (TH1F*)f2->Get("tau_stat_err2");
00078   TH1F *neg_QCD_syst_low_err = (TH1F*)f2->Get("QCD_syst_low_err");
00079   TH1F *neg_QCD_syst_high_err = (TH1F*)f2->Get("QCD_syst_high_err");
00080 
00081   float pos_high[24], pos_low[24];
00082   float neg_high[24], neg_low[24];
00083   float bin_cent[24], pos_bin_val[24], neg_bin_val[24];
00084   float ex[24];
00085   float err;
00086 
00087   // calculate the statistical uncertainty and the background
00088   // shape systematic uncertainty 
00089   for (int i=1; i<=24; i++) {
00090     ex[i-1] = 2.0; 
00091     err = 0.;
00092     err += pos_raw_stat_err2->GetBinContent(2*i)+
00093            pos_raw_stat_err2->GetBinContent(2*i+1);
00094     err += pos_QCD_stat_err2->GetBinContent(2*i)+
00095            pos_QCD_stat_err2->GetBinContent(2*i+1);
00096     err += pos_eemc_stat_err2->GetBinContent(2*i)+
00097            pos_eemc_stat_err2->GetBinContent(2*i+1);
00098     err += pos_tau_stat_err2->GetBinContent(2*i)+
00099            pos_tau_stat_err2->GetBinContent(2*i+1);
00100     pos_data_repack->SetBinError(i,sqrt(err));
00101     err = 0.;
00102     err += neg_raw_stat_err2->GetBinContent(2*i)+
00103            neg_raw_stat_err2->GetBinContent(2*i+1);
00104     err += neg_QCD_stat_err2->GetBinContent(2*i)+
00105            neg_QCD_stat_err2->GetBinContent(2*i+1);
00106     err += neg_eemc_stat_err2->GetBinContent(2*i)+
00107            neg_eemc_stat_err2->GetBinContent(2*i+1);
00108     err += neg_tau_stat_err2->GetBinContent(2*i)+
00109            neg_tau_stat_err2->GetBinContent(2*i+1);
00110     neg_data_repack->SetBinError(i,sqrt(err));
00111 
00112     bin_cent[i-1] = pos_data_repack->GetBinCenter(i);
00113     pos_bin_val[i-1] = pos_data_repack->GetBinContent(i);
00114     neg_bin_val[i-1] = neg_data_repack->GetBinContent(i);
00115 
00116     // The QCD background shape systematic errors are
00117     // put in their arrays
00118     pos_high[i-1] = pos_QCD_syst_high_err->GetBinContent(2*i)+
00119                     pos_QCD_syst_high_err->GetBinContent(2*i+1);
00120     neg_high[i-1] = neg_QCD_syst_high_err->GetBinContent(2*i)+
00121                     neg_QCD_syst_high_err->GetBinContent(2*i+1);
00122     pos_low[i-1]  = pos_QCD_syst_low_err->GetBinContent(2*i)+
00123                     pos_QCD_syst_low_err->GetBinContent(2*i+1);
00124     neg_low[i-1]  = neg_QCD_syst_low_err->GetBinContent(2*i)+
00125                     neg_QCD_syst_low_err->GetBinContent(2*i+1);
00126   }
00127 
00128   // *******************************************************
00129   // Calculate an energy scale systematic
00130   // *******************************************************
00131 
00132   TFile *f5 = new TFile("bkgd_histos_pos_energy_up.root");
00133   TFile *f6 = new TFile("bkgd_histos_pos_energy_down.root");
00134   TFile *f7 = new TFile("bkgd_histos_neg_energy_up.root");
00135   TFile *f8 = new TFile("bkgd_histos_neg_energy_down.root");
00136 
00137   TH1F *pos_data_en_up   = (TH1F*)f5->Get("signal_final3");
00138   TH1F *pos_data_en_down = (TH1F*)f6->Get("signal_final3");
00139   TH1F *neg_data_en_up   = (TH1F*)f7->Get("signal_final3");
00140   TH1F *neg_data_en_down = (TH1F*)f8->Get("signal_final3");
00141 
00142   TH1F *pos_data_en_up_repack = new TH1F("pos_data_en_up_repack","pos_data_en_up_repack",24,3.,99.);
00143   TH1F *pos_data_en_down_repack = new TH1F("pos_data_en_down_repack","pos_data_en_down_repack",24,3.,99.);
00144   TH1F *neg_data_en_up_repack = new TH1F("neg_data_en_up_repack","neg_data_en_up_repack",24,3.,99.);
00145   TH1F *neg_data_en_down_repack = new TH1F("neg_data_en_down_repack","neg_data_en_down_repack",24,3.,99.);
00146 
00147   for (int i=1; i<=24; i++) {
00148     pos_data_en_up_repack->SetBinContent(i,pos_data_en_up->GetBinContent(2*i)+
00149                                            pos_data_en_up->GetBinContent(2*i+1));
00150     pos_data_en_down_repack->SetBinContent(i,pos_data_en_down->GetBinContent(2*i)+
00151                                              pos_data_en_down->GetBinContent(2*i+1));
00152     neg_data_en_up_repack->SetBinContent(i,neg_data_en_up->GetBinContent(2*i)+
00153                                            neg_data_en_up->GetBinContent(2*i+1));
00154     neg_data_en_down_repack->SetBinContent(i,neg_data_en_down->GetBinContent(2*i)+
00155                                              neg_data_en_down->GetBinContent(2*i+1));
00156 
00157   }
00158 
00159   for (int i=1; i<=24; i++) {
00160     cout << "W+ -> " << pos_data_en_up_repack->GetBinContent(i) << " " << pos_data_repack->GetBinContent(i) << " " << pos_data_en_down_repack->GetBinContent(i) << endl;
00161   }
00162   for (int i=1; i<=24; i++) {
00163     cout << "W- -> " << neg_data_en_up_repack->GetBinContent(i) << " " << neg_data_repack->GetBinContent(i) << " " << neg_data_en_down_repack->GetBinContent(i) << endl;
00164   }
00165 
00166   // Now loop over all the bins of the three histograms (energy low, 
00167   // normal, energy high) and in each bin find the lowest and the
00168   // highest and then call that the systematic (normally it will be
00169   // the two energy scaled data sets) ... and do it for both charges
00170   float pos_data_ensys_low[24];
00171   float pos_data_ensys_high[24];
00172   float neg_data_ensys_low[24];
00173   float neg_data_ensys_high[24];
00174   for (int i=1; i<24; i++) {
00175 
00176     // for positive charge
00177     float low = 10000.;
00178     if (pos_data_en_up_repack->GetBinContent(i) < low) {
00179       low = pos_data_en_up_repack->GetBinContent(i);
00180     }
00181     if (pos_data_repack->GetBinContent(i) < low) {
00182       low = pos_data_repack->GetBinContent(i);
00183     }
00184     if (pos_data_en_down_repack->GetBinContent(i) < low) {
00185       low = pos_data_en_down_repack->GetBinContent(i);
00186     }
00187     pos_data_ensys_low[i-1] = fabs(low-pos_data_repack->GetBinContent(i));
00188 
00189     float high = 0.;
00190     if (pos_data_en_up_repack->GetBinContent(i) > high) {
00191       high = pos_data_en_up_repack->GetBinContent(i);
00192     }
00193     if (pos_data_repack->GetBinContent(i) > high) {
00194       high = pos_data_repack->GetBinContent(i);
00195     }
00196     if (pos_data_en_down_repack->GetBinContent(i) > high) {
00197       high = pos_data_en_down_repack->GetBinContent(i);
00198     }
00199     pos_data_ensys_high[i-1] = fabs(high-pos_data_repack->GetBinContent(i));
00200 
00201     // for negative charge
00202     low = 10000.;
00203     if (neg_data_en_up_repack->GetBinContent(i) < low) {
00204       low = neg_data_en_up_repack->GetBinContent(i);    }
00205     if (neg_data_repack->GetBinContent(i) < low) {
00206       low = neg_data_repack->GetBinContent(i);
00207     }
00208     if (neg_data_en_down_repack->GetBinContent(i) < low) {
00209       low = neg_data_en_down_repack->GetBinContent(i);
00210     }
00211     neg_data_ensys_low[i-1] = fabs(low-neg_data_repack->GetBinContent(i));
00212     high = 0.;
00213     if (neg_data_en_up_repack->GetBinContent(i) > high) {      high = neg_data_en_up_repack->GetBinContent(i);
00214     }
00215     if (neg_data_repack->GetBinContent(i) > high) {
00216       high = neg_data_repack->GetBinContent(i);
00217     }
00218     if (neg_data_en_down_repack->GetBinContent(i) > high) {
00219       high = neg_data_en_down_repack->GetBinContent(i);
00220     }
00221     neg_data_ensys_high[i-1] = fabs(high-neg_data_repack->GetBinContent(i));
00222 
00223     cout << i << " " << pos_data_ensys_low[i-1] << " " << pos_data_ensys_high[i-1] << 
00224 " " << neg_data_ensys_low[i-1] << " " << neg_data_ensys_high[i-1] << endl;
00225  
00226     // convention is switched (add to the QCD background systematic)
00227     // assume fully correlated (probably not true)
00228     pos_high[i-1] += pos_data_ensys_low[i-1];
00229     pos_low[i-1]  += pos_data_ensys_high[i-1];
00230     neg_high[i-1] += neg_data_ensys_low[i-1];
00231     neg_low[i-1]  += neg_data_ensys_high[i-1];
00232   }
00233 
00234   TGraphAsymmErrors *pos_syst_err_box = new TGraphAsymmErrors(24,bin_cent,pos_bin_val,ex,ex,pos_high,pos_low);
00235   TGraphAsymmErrors *neg_syst_err_box = new TGraphAsymmErrors(24,bin_cent,neg_bin_val,ex,ex,neg_high,neg_low);
00236 
00237   // ******************************************************
00238   // Calculate the MC to Data ratios now (even though they
00239   // are not plotted)
00240   // ******************************************************
00241 
00242   TH1F *pos_ratio = pos_data_repack->Clone();
00243   TH1F *neg_ratio = neg_data_repack->Clone();
00244   float pos_ratio_syst_high[24], pos_ratio_syst_low[24]; 
00245   float neg_ratio_syst_high[24], neg_ratio_syst_low[24]; 
00246   float pos_ratio_bin_val[24], neg_ratio_bin_val[24];
00247 
00248   for (int i=1; i<=24; i++) {
00249     if ((pos_mc_repack->GetBinContent(i) > 0)  && (pos_data_repack->GetBinContent(i) > 0)) {
00250       pos_ratio->SetBinContent(i,pos_data_repack->GetBinContent(i)/
00251                                  pos_mc_repack->GetBinContent(i));
00252       pos_ratio->SetBinError(i,pos_data_repack->GetBinError(i)/pos_mc_repack->GetBinContent(i));
00253       pos_ratio_bin_val[i-1] = pos_data_repack->GetBinContent(i)/
00254                                pos_mc_repack->GetBinContent(i);
00255       pos_ratio_syst_high[i-1] = (pos_data_repack->GetBinContent(i)+pos_high[i-1])/pos_mc_repack->GetBinContent(i);
00256       pos_ratio_syst_low[i-1] = (pos_data_repack->GetBinContent(i)-pos_low[i-1])/pos_mc_repack->GetBinContent(i);
00257     } else {
00258       pos_ratio->SetBinContent(i,0.);
00259       pos_ratio->SetBinError(i,0.);
00260       pos_ratio_syst_high[i-1] = 0.;
00261       pos_ratio_syst_low[i-1] = 0.;
00262       pos_ratio_bin_val[i-1] = 0.;
00263     }
00264     if ((neg_mc_repack->GetBinContent(i) > 0) && (neg_data_repack->GetBinContent(i) > 0)) {
00265       neg_ratio->SetBinContent(i,neg_data_repack->GetBinContent(i)/
00266                                  neg_mc_repack->GetBinContent(i));
00267       neg_ratio->SetBinError(i,neg_data_repack->GetBinError(i)/neg_mc_repack->GetBinContent(i));
00268       neg_ratio_bin_val[i-1] = neg_data_repack->GetBinContent(i)/
00269                              neg_mc_repack->GetBinContent(i);
00270       neg_ratio_syst_high[i-1] = (neg_data_repack->GetBinContent(i)+neg_high[i-1])/neg_mc_repack->GetBinContent(i);
00271       neg_ratio_syst_low[i-1] = (neg_data_repack->GetBinContent(i)-neg_low[i-1])/neg_mc_repack->GetBinContent(i);
00272     } else {
00273       neg_ratio->SetBinContent(i,0.);
00274       neg_ratio->SetBinError(i,0.);
00275       neg_ratio_syst_high[i-1] = 0.;
00276       neg_ratio_syst_low[i-1] = 0.;
00277       neg_ratio_bin_val[i-1] = 0.;
00278     }
00279   }
00280 
00281   for (int i=0; i<24; i++) {
00282     pos_ratio_syst_high[i] -= pos_ratio_bin_val[i];
00283     pos_ratio_syst_low[i] = pos_ratio_bin_val[i]-pos_ratio_syst_low[i];
00284     neg_ratio_syst_high[i] -= neg_ratio_bin_val[i];
00285     neg_ratio_syst_low[i] = neg_ratio_bin_val[i]-neg_ratio_syst_low[i];
00286   }
00287   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);
00288   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);
00289 
00290   // ******************************************************
00291   // Plots versus ET (no longer the ratios)
00292   // ******************************************************
00293   gStyle->SetTitleBorderSize(0);
00294   gStyle->SetOptDate(0); 
00295   TCanvas *can1 = new TCanvas("can1","can1",0,0,500,500);
00296 
00297   gPad->SetBorderMode(0);
00298   gPad->SetGridx(0);
00299   gPad->SetGridy(0);
00300   pos_mc_repack->GetYaxis()->SetRangeUser(-10.,180.);
00301   pos_mc_repack->GetXaxis()->SetRangeUser(0.,70.);
00302   pos_mc_repack->GetXaxis()->SetTitle("EMC Cluster E_{T} (GeV)");
00303   pos_mc_repack->GetYaxis()->SetTitleOffset(1.2);
00304   pos_mc_repack->GetYaxis()->SetTitle("Counts");
00305   pos_mc_repack->SetLineColor(2);
00306   pos_mc_repack->SetStats(kFALSE);
00307   pos_mc_repack->SetTitle("");
00308   pos_mc_repack->SetLineWidth(3.*pos_mc_repack->GetLineWidth());
00309   pos_data_repack->SetLineWidth(2.*pos_data_repack->GetLineWidth());
00310   pos_mc_repack->Draw();
00311   pos_syst_err_box->SetFillColor(14);
00312   pos_syst_err_box->Draw("e2same");
00313   pos_mc_repack->Draw("same");
00314   pos_data_repack->Draw("same");
00315   TLegend *leg1 = new TLegend(0.5,0.7,0.88,0.88);
00316   leg1->SetFillColor(0);
00317   leg1->SetBorderSize(0);
00318   leg1->AddEntry(pos_mc_repack,"MC Norm. to Signal","l");
00319   leg1->AddEntry(pos_data_repack,"Signal","le");
00320   leg1->AddEntry(pos_syst_err_box,"Syst. Unc.","f");
00321   leg1->Draw("same");
00322   TLatex *lat1 = new TLatex();
00323   lat1->SetNDC();
00324   lat1->DrawLatex(0.14,0.84,"p+p#rightarrow W^{+}#rightarrow e^{+}+#nu");
00325   lat1->DrawLatex(0.14,0.74,"Run9 STAR");
00326   lat1->DrawLatex(0.14,0.69,"Preliminary");
00327   lat1->DrawLatex(0.14,0.64,"#sqrt{s}=500 GeV");
00328   can1->Print("data_mc_ET_comp_output4_W+_4GeV.eps");
00329 
00330   TCanvas *can2 = new TCanvas("can2","can2",0,0,500,500);
00331   gPad->SetBorderMode(0);
00332   gPad->SetGridx(0);
00333   gPad->SetGridy(0);
00334   neg_mc_repack->GetXaxis()->SetRangeUser(0.,70.);
00335   neg_mc_repack->GetYaxis()->SetRangeUser(-10.,100.);
00336   neg_mc_repack->GetXaxis()->SetTitle("EMC Cluster E_{T} (GeV)");
00337   neg_mc_repack->GetYaxis()->SetTitleOffset(1.2);
00338   neg_mc_repack->GetYaxis()->SetTitle("Counts");
00339   neg_mc_repack->SetLineColor(2); 
00340   neg_mc_repack->SetStats(kFALSE);
00341   neg_mc_repack->SetTitle("");
00342   neg_mc_repack->SetLineWidth(3.*neg_mc_repack->GetLineWidth());
00343   neg_data_repack->SetLineWidth(2.*neg_data_repack->GetLineWidth());
00344   neg_mc_repack->Draw();
00345   neg_syst_err_box->SetFillColor(14);
00346   neg_syst_err_box->Draw("e2same");
00347   neg_mc_repack->Draw("same");
00348   neg_data_repack->Draw("same"); 
00349   leg1->Draw("same");
00350   lat1->DrawLatex(0.14,0.84,"p+p#rightarrow W^{-}#rightarrow e^{-}+#bar{#nu}");
00351   lat1->DrawLatex(0.14,0.74,"Run9 STAR");
00352   lat1->DrawLatex(0.14,0.69,"Preliminary");
00353   lat1->DrawLatex(0.14,0.64,"#sqrt{s}=500 GeV");
00354   can2->Print("data_mc_ET_comp_output4_W-_4GeV.eps");
00355 
00356 } // end of macro
StSpinPool/StWalgoB2009/macros/plotMCDataComp_4GeV.C
