StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
makeBkgdFiles.C
1 // ********************************************************
2 // This code makes the necessary root histos for the
3 // background subtraction and systematic uncertainty calculations
4 // Set charge == 1 for postive charges
5 // Set charge == -1 for negative charges
6 // Set charge == 0 for the sum of both
7 //
8 // Set two_or_four == 2 for 2 GeV wide bin histograms
9 // Set two_or_four == 4 for 4 GeV wide bin hisrograms
10 // ********************************************************
11 
12 void makeBkgdFiles(int charge, int two_or_four) {
13 
14  // ******************************************************
15  // Read in the data set and all the histograms needed
16  // to do the actual calculation
17  // ******************************************************
18 
19  gStyle->SetOptDate(0);
20  TFile *f1 = new TFile("run9setP1234.wana.hist.root");
21 
22  // hybrid2
23  //TFile *f1 = new TFile("run9setP1234.wana.hist.root");
24  //TFile *f1 = new TFile("noZTag/run9setP1234.wana.hist.root");
25 
26  // get the signal and missing endcap backgrounds
27  if (charge == 1) {
28  TH1F *signal = (TH1F*)f1->Get("pos_muclustpTbal_wE");
29  TH1F *signal_wo_eemc = (TH1F*)f1->Get("pos_muclustpTbal_noE");
30  } else if (charge == -1) {
31  TH1F *signal = (TH1F*)f1->Get("neg_muclustpTbal_wE");
32  TH1F *signal_wo_eemc = (TH1F*)f1->Get("neg_muclustpTbal_noE");
33  } else if (charge == 0) {
34  TH1F *signal = (TH1F*)f1->Get("muclustPtBal");
35  TH1F *signal_wo_eemc = (TH1F*)f1->Get("muclustPtBalnoE");
36  }
37 
38  // ******************************************************
39  // Read in all the MC data sets for background
40  // subtractions and studies
41  // ******************************************************
42 
43  TFile *MC_fs[6];
44 
45  /*MC_fs[1] = new TFile("rcf10010.wana.hist.root"); // W+ -> e++nu
46  MC_fs[0] = new TFile("rcf10011.wana.hist.root"); // W- -> e-+nu
47  MC_fs[2] = new TFile("rcf10012.wana.hist.root"); // W -> tau+nu
48  MC_fs[3] = new TFile("rcf10014.wana.hist.root"); // W -> any
49  MC_fs[4] = new TFile("rcf10015.wana.hist.root"); // Z -> any
50 */
51 
52 
53  /*MC_fs[0] = new TFile("rck10010.wana.hist.root"); // W+ -> e++nu
54  MC_fs[1] = new TFile("rck10011.wana.hist.root"); // W- -> e-+nu
55  MC_fs[2] = new TFile("rck10012.wana.hist.root"); // W -> tau+nu
56  MC_fs[3] = new TFile("rck10017.wana.hist.root"); // W -> any
57  MC_fs[4] = new TFile("rck10018.wana.hist.root"); // Z -> any
58 */
59 
60  MC_fs[0] = new TFile("rcn10010.wana.hist.root"); // W+ -> e++nu
61  MC_fs[1] = new TFile("rcn10011.wana.hist.root"); // W- -> e-+nu
62  MC_fs[2] = new TFile("rcn10012.wana.hist.root"); // W -> tau+nu
63  MC_fs[3] = new TFile("rcn10017.wana.hist.root"); // W -> any
64  MC_fs[4] = new TFile("rcn10018.wana.hist.root"); // Z -> any
65 
66 
67 /*
68  MC_fs[0] = new TFile("noZTag/rck10010.wana.hist.root"); // W+ -> e++nu
69  MC_fs[1] = new TFile("noZTag/rck10011.wana.hist.root"); // W- -> e-+nu
70  MC_fs[2] = new TFile("noZTag/rck10012.wana.hist.root"); // W -> tau+nu
71  MC_fs[3] = new TFile("noZTag/rck10017.wana.hist.root"); // W -> any
72  MC_fs[4] = new TFile("noZTag/rck10018.wana.hist.root"); // Z -> any
73 */
74 
75  // Because of a screw up between the RFF and FF in the simulation
76  // and reconstruction the charge sign of the MC needs to be inverted
77  TH1F *MC_dists_raw[5][3];
78  for (int i=0; i<5; i++) {
79  if (charge == 1) {
80  MC_dists_raw[i][0] = (TH1F*)MC_fs[i]->Get("pos_muclustpTbal_wE");
81  MC_dists_raw[i][1] = (TH1F*)MC_fs[i]->Get("pos_muclustpTbal_noE");
82  //MC_dists_raw[i][2] = (TH1F*)MC_fs[i]->Get("pos_muclustpTbal_bckgrd");
83  } else if (charge == -1) {
84  MC_dists_raw[i][0] = (TH1F*)MC_fs[i]->Get("neg_muclustpTbal_wE");
85  MC_dists_raw[i][1] = (TH1F*)MC_fs[i]->Get("neg_muclustpTbal_noE");
87  } else if (charge == 0) {
88  MC_dists_raw[i][0] = (TH1F*)MC_fs[i]->Get("muclustPtBal");
89  MC_dists_raw[i][1] = (TH1F*)MC_fs[i]->Get("muclustPtBalnoE");
90  MC_dists_raw[i][2] = (TH1F*)MC_fs[i]->Get("muclustPtBal_bckgrd");
91  }
92  MC_dists_raw[i][2] = (TH1F*)MC_fs[i]->Get("muclustPtBal_bckgrd");
93 
94  }
95 
96  // Old lumi estimates
97  //float lumi[5] = {116.5,318.5,76.9,1220.,45.6}; // old MC
98  //float lumi[5] = {30.,63.,61.,14.7,33.};
99 
100  // w/ Z and W xsec scaled to RHICBOS from PYTHiA
101  //float lumi[5] = {30./(124./98.7),63./(41./32.9),61./((124.+41.)/(98.7+32.9)),14.7/((124+41)/(98.7+32.9)),33./(10./8.23)};
102 
103  // w/ Z and W xsec scaled to RHICBOS from PYTHiA - new MC
104  float lumi[5] = {30.4/(124./98.7),31.4/(41./32.9),30.3/((124.+41.)/(98.7+32.9)),30.8/((124+41)/(98.7+32.9)),33.4/(10./8.23)};
105 
106  float lumi_fact[6];
107  for (int i=0; i<5; i++) {lumi_fact[i] = 12.0/lumi[i];}
108 
109  for (int i=0; i<5; i++) {
110  for (int j=0; j<3; j++) {
111  MC_dists_raw[i][j]->Scale(lumi_fact[i]);
112  }
113  }
114 
115 
116  // Repack the histograms to mesh with the odd staggered binning that
117  // is being used
118  char str[200];
119  TH1F *MC_dists_repack[5][3];
120  for (int i=0; i<5; i++) {
121  sprintf(str,"mcclustPtBal_%d",i);
122  MC_dists_repack[i][0] = new TH1F(str,str,49,1.,99.);
123  sprintf(str,"mcclustPtBalnoE_%d",i);
124  MC_dists_repack[i][1] = new TH1F(str,str,49,1.,99.);
125  sprintf(str,"mcclustPtBal_bckgrd_%d",i);
126  MC_dists_repack[i][2] = new TH1F(str,str,49,1.,99.);
127  }
128 
129  for (int i=0; i<5; i++) {
130  for (int j=0; j<3; j++) {
131  for (int k=1; k<=49; k++) {
132  MC_dists_repack[i][j]->SetBinContent(k,MC_dists_raw[i][j]->GetBinContent(2*k)+
133  MC_dists_raw[i][j]->GetBinContent(2*k+1));
134  }
135  }
136  }
137 
138  // Make the EEMC background the real one of interest now
139  for (int i=0; i<5; i++) {
140  MC_dists_repack[i][1]->Add(MC_dists_repack[i][0],-1);
141  MC_dists_raw[i][1]->Add(MC_dists_raw[i][0],-1.);
142  }
143 
144  // **********************************************
145  // Do some otherstuff
146  // **********************************************
147 
148  TH1F *eemc_bkgd = signal_wo_eemc->Clone();
149  eemc_bkgd->Add(signal,-1.);
150  TH1F *signal_final = signal->Clone();
151  signal_final->Add(eemc_bkgd,-1);
152 
153  TH1F *eemc_bkgd2 = new TH1F("eemc_bkgd2","eemc_bkgd2",49,1,99);
154  TH1F *zsig_bkgd2 = new TH1F("zsig_bkgd2","zsig_bkgd2",49,1,99);
155  TH1F *zeemc_bkgd2 = new TH1F("zeemc_bgkd2","zeemc_bkgd2",49,1,99);
156  TH1F *zback_bkgd2 = new TH1F("zback_bkgd2","zback_bkgd2",49,1,99);
157 
158  TH1F *wanysig_bkgd2 = new TH1F("wanysig_bkgd2","wanysig_bkgd2",49,1,99);
159  TH1F *wanyeemc_bkgd2 = new TH1F("wanyeemc_bgkd2","wanyeemc_bkgd2",49,1,99);
160  TH1F *wanyback_bkgd2 = new TH1F("wanyback_bkgd2","wanyback_bkgd2",49,1,99);
161 
162  TH1F *zsig = MC_dists_raw[4][0]->Clone();
163  TH1F *zeemc = MC_dists_raw[4][1]->Clone();
164  TH1F *zback = MC_dists_raw[4][2]->Clone();
165 
166  signal_final->Add(zsig,-1.);
167 
168  // First get the "nominal" background shape
169  TH1F *bkgd_shape1 = (TH1F*)f1->Get("muclustPtBal_bckgrd");
170 
171  TH1F *bkgd_shape_nom = new TH1F("bkgd_shape","bkgd_shape",49,1,99);
172  for (int i=1; i<=49; i++) {
173  bkgd_shape_nom->SetBinContent(i,bkgd_shape1->GetBinContent(2*i)+
174  bkgd_shape1->GetBinContent(2*i+1));
175  }
176  TH1F *bkgd_shape_nom2 = bkgd_shape_nom->Clone();
177 
178  TH1F *signal_final2 = new TH1F("signal_final2","signal_final2",49,1,99);
179  signal_final2->SetLineColor(2);
180  signal_final2->SetLineWidth(2.*signal_final2->GetLineWidth());
181  TH1F *signal2 = new TH1F("signal2","signal2",49,1,99);
182  for (int i=1; i<=49; i++) {
183  signal_final2->SetBinContent(i,signal_final->GetBinContent(2*i)+
184  signal_final->GetBinContent(2*i+1));
185  signal2->SetBinContent(i,signal->GetBinContent(2*i)+
186  signal->GetBinContent(2*i+1));
187  eemc_bkgd2->SetBinContent(i,eemc_bkgd->GetBinContent(2*i)+
188  eemc_bkgd->GetBinContent(2*i+1));
189  zsig_bkgd2->SetBinContent(i,zsig->GetBinContent(2*i)+
190  zsig->GetBinContent(2*i+1));
191  zeemc_bkgd2->SetBinContent(i,zeemc->GetBinContent(2*i)+
192  zeemc->GetBinContent(2*i+1));
193  zback_bkgd2->SetBinContent(i,zback->GetBinContent(2*i)+
194  zback->GetBinContent(2*i+1));
195  wanysig_bkgd2->SetBinContent(i,MC_dists_raw[3][0]->GetBinContent(2*i)+
196  MC_dists_raw[3][0]->GetBinContent(2*i+1));
197  wanyeemc_bkgd2->SetBinContent(i,MC_dists_raw[3][1]->GetBinContent(2*i)+
198  MC_dists_raw[3][1]->GetBinContent(2*i+1));
199  wanyback_bkgd2->SetBinContent(i,MC_dists_raw[3][2]->GetBinContent(2*i)+
200  MC_dists_raw[3][2]->GetBinContent(2*i+1));
201  }
202 
203  TCanvas *can2 = new TCanvas("can2","can2",0,0,600,400);
204  signal2->Draw();
205  signal_final2->Draw("same");
206  //can2->Print("signal_w_eemc.eps");
207  //can2->Print("signal_w_eemc.png");
208 
209 
210  // **********************************************
211  // plot all the business
212  // **********************************************
213  eemc_bkgd->SetLineColor(8);
214  eemc_bkgd->SetLineWidth(2.*eemc_bkgd->GetLineWidth());
215  bkgd_shape_nom2->SetLineColor(4);
216  bkgd_shape_nom2->SetLineWidth(2.*bkgd_shape_nom2->GetLineWidth());
217 
218  TCanvas *can8 = new TCanvas("can8","can8",0,0,3000,1800);
219  can8->Divide(5,3);
220  for (int i=0; i<5; i++) {
221  for (int j=0; j<3; j++) {
222  can8->cd(5*(j)+i+1);
223  //cout << 5*(j)+i+1 << endl;
224  gPad->SetGridx(0);
225  gPad->SetGridy(0);
226  if (j == 0) {
227  signal_final2->Draw();
228  } else if (j == 1) {
229  gPad->SetLogy();
230  eemc_bkgd->Draw();
231  } else if (j == 2) {
232  gPad->SetLogy();
233  bkgd_shape_nom2->Draw();
234  }
235  MC_dists_repack[i][j]->Draw("same");
236  }
237  }
238  //can8->Print("phys_bkgds.eps");
239  //can8->Print("phys_bkgds.png");
240 
241  for (int i=0; i<5; i++) {
242  float signal_sum = 0.;
243  float eemc_sum = 0.;
244  float back_sum = 0.;
245  for (int j=13; j<=49; j++) {
246  signal_sum += MC_dists_repack[i][0]->GetBinContent(j);
247  eemc_sum += MC_dists_repack[i][1]->GetBinContent(j);
248  back_sum += MC_dists_repack[i][2]->GetBinContent(j);
249  }
250  }
251 
252  //Tau "background" needs no further scaling
253  float taufrac=1.5;
254  // calculate a fractional scale factor for the charge separated
255  // background due to tau->e
256  //float tauplus_xsec = 9.8e-8;
257  //float tauminus_xsec = 3.1e-8;
258  //float taufrac;
259  //if (charge == 1) {
260  // taufrac = tauplus_xsec/(tauplus_xsec+tauminus_xsec);
261  //} else if (charge == -1) {
262  // taufrac = tauminus_xsec/(tauplus_xsec+tauminus_xsec);
263  //} else if (charge == 0) {
264  // taufrac = 1.;
265  //}
266 
267 
268  // ******************************************************
269  // Do the iterative normalization of the W signal to
270  // to the background for the nominal shape
271  // ******************************************************
272 
273  // subtract off the tau background
274  signal_final2->Add(MC_dists_repack[2][0],-1.*taufrac);
275 
276  // Now lets try a new background function
277  // Fit the peak with a line from 23<ET<39
278  TF1 *func1 = new TF1("func1","[0]+[1]*x",23,39);
279  func1->SetParameter(0,0.);
280  func1->SetParameter(1,0.);
281 
282  TCanvas *can4 = new TCanvas("can4","can4",0,0,600,400);
283  signal_final2->Draw();
284 
285  // Calculate the background for the nominal signal
286  // by doing the 20 iterations over the signal peak
287  float signal_in_norm[50];
288  TH1F *bkgd_shape_unnorm[20];
289  TH1F *signal_for_new[20];
290  for (int i=0; i<20; i++) {
291  bkgd_shape_unnorm[i] = (TH1F*)bkgd_shape_nom->Clone();
292  signal_for_new[i] = (TH1F*)signal_final2->Clone();
293 
294  // subtract off the Z signal from the background
295  bkgd_shape_unnorm[i]->Add(zback_bkgd2,-1.);
296 
297  // calculate the W signal in the normalization bins
298  signal_in_norm[8] = func1->Integral(15,17);
299  signal_in_norm[9] = func1->Integral(17,19);
300  signal_in_norm[10] = func1->Integral(19,21);
301  //cout << "Integrals = " << signal_in_norm[8] << " " << signal_in_norm[9] << " " << signal_in_norm[10] << endl;
302  for (int j=8; j<=10; j++) {
303  if (signal_in_norm[j] < 0) {signal_in_norm[j] = 0.;}
304  }
305 
306  // calculate the normalization factor using the
307  // signal subtraction normalization bins
308  float normt = 0, normb = 0.;
309  for (int k=8; k<10; k++) {
310  if (bkgd_shape_unnorm[i]->GetBinContent(k) > 0) {
311  normt += signal_final2->GetBinContent(k)-signal_in_norm[k];
312  normb += bkgd_shape_unnorm[i]->GetBinContent(k);
313  }
314  }
315  if (normb > 0 && normt > 0) {
316  float norm = normt/normb;
317  bkgd_shape_unnorm[i]->Scale(norm);
318  bkgd_shape_unnorm[i]->Draw("same");
319  //cout << "norm " << i << " = " << norm << endl;
320  }
321 
322  // With the new signal estimate, calculate the normalization
323  // factor
324  for (int j=1; j<=49; j++) {
325  if (bkgd_shape_unnorm[i]->GetBinContent(j) < 0) {bkgd_shape_unnorm[i]->SetBinContent(j,0.);}
326  }
327  signal_for_new[i]->Add(bkgd_shape_unnorm[i],-1.);
328  signal_for_new[i]->Fit(func1,"RQ");
329  }
330  //can4->Print("new_norm.eps");
331  //can4->Print("new_norm.png");
332 
333  TH1F *signal_in_norm_region = new TH1F("signal_in_norm_region","signal_in_norm_region",49,1.,99.);
334 
335  signal_in_norm_region->SetBinContent(8,signal_in_norm[8]);
336  signal_in_norm_region->SetBinContent(9,signal_in_norm[9]);
337 
338  TH1F *new_bkgd = new TH1F("new_bkgd","new_bkgd",49,1.,99.);
339  new_bkgd = (TH1F*)bkgd_shape_unnorm[19]->Clone();
340  new_bkgd->SetName("new_bkgd");
341 
342  TCanvas *can5 = new TCanvas("can5","can5",0,0,600,400);
343  signal_final2->Draw();
344  new_bkgd->Draw("same");
345  signal_in_norm_region->Draw("same");
346  //can5->Print("signal_new.eps");
347  //can5->Print("signal_new.png");
348 
349  // ******************************************************
350  // Calculate all the 1200 shapes for the background
351  // systematic shape study
352  // ******************************************************
353 
354  // get the various background shapes made by
355  // using different cuts for the awayside backgrounds
356  TH2F *bkgd_hists_from_file[21];
357  TH2F *bkgd_hists_from_file2[21];
358  char str[200];
359  for (int i=0; i<=20; i++) {
360 
361 // if (charge == 1) {
362 // sprintf(str,"pos_failAwaySide_Awayside_pt_bin_%d",i);
363 // bkgd_hists_from_file[i] = (TH2F*)f1->Get(str);
364 // } else if (charge == -1) {
365 // sprintf(str,"neg_failAwaySide_Awayside_pt_bin_%d",i);
366 // bkgd_hists_from_file[i] = (TH2F*)f1->Get(str);
367 // } else if (charge == 0) {
368  sprintf(str,"pos_failAwaySide_Awayside_pt_bin_%d",i);
369  bkgd_hists_from_file[i] = (TH2F*)f1->Get(str);
370  sprintf(str,"neg_failAwaySide_Awayside_pt_bin_%d",i);
371  bkgd_hists_from_file2[i] = (TH2F*)f1->Get(str);
372  bkgd_hists_from_file[i]->Add(bkgd_hists_from_file2[i]);
373 // }
374  }
375 
376  // Now do the rebinning
377  TH1F *bkgd_hists1[21][21];
378  TH1F *bkgd_hists2[21][21];
379  TH1F *bkgd_hists3[21][21];
380  int jval=0;
381  for (int i=0; i<=20; i++) {
382  int j=jval;
383 // for (int j=0; j<20; j++) {
384  sprintf(str,"bkgd_hist1_%d_%d",i,j);
385  bkgd_hists1[i][j] = new TH1F(str,str,49,1,99);
386  sprintf(str,"bkgd_hist2_%d_%d",i,j);
387  bkgd_hists2[i][j] = new TH1F(str,str,49,1,99);
388  sprintf(str,"bkgd_hist3_%d_%d",i,j);
389  bkgd_hists3[i][j] = new TH1F(str,str,49,1,99);
390  for (int k=1; k<=49; k++) {
391  bkgd_hists1[i][j]->SetBinContent(k,bkgd_hists_from_file[0]->GetBinContent(2*k,i+1)+bkgd_hists_from_file[0]->GetBinContent(2*k+1,i+1));
392  bkgd_hists2[i][j]->SetBinContent(k,bkgd_hists_from_file[0]->GetBinContent(2*k,i+1)+bkgd_hists_from_file[0]->GetBinContent(2*k+1,i+1));
393  bkgd_hists3[i][j]->SetBinContent(k,bkgd_hists_from_file[0]->GetBinContent(2*k,i+1)+bkgd_hists_from_file[0]->GetBinContent(2*k+1,i+1));
394  }
395  // }
396  }
397 
398 
399 
400  // initiaize the iterative fit functions
401  TF1 *func2 = new TF1("func2","[0]+[1]*x",23,39);
402  func2->SetParameter(0,0.);
403  func2->SetParameter(1,0.);
404  TF1 *func3 = new TF1("func3","[0]+[1]*x",23,39);
405  func3->SetParameter(0,0.);
406  func3->SetParameter(1,0.);
407 
408  // Now loop over the the 1200 possibilities (400 loops and 3 normalization
409  // regions)
410  float final_sig_in_norm[21][21][3];
411  float final_chisquare[21][21];
412  float signal_in_norm1[50];
413  float signal_in_norm2[50];
414  float signal_in_norm3[50];
415  TH1F *new_bkgd_hists1[21][21];
416  TH1F *new_bkgd_hists2[21][21];
417  TH1F *new_bkgd_hists3[21][21];
418  TH1F *bkgd_shape_unnorm1[20];
419  TH1F *bkgd_shape_unnorm2[20];
420  TH1F *bkgd_shape_unnorm3[20];
421  TH1F *signal_for_new1[20];
422  TH1F *signal_for_new2[20];
423  TH1F *signal_for_new3[20];
424  for (int i=0; i<=20; i++) {
425  int j=jval;
426 // for (int j=0; j<20; j++) {
427 
428  func1->SetParameter(0,0.);
429  func1->SetParameter(1,0.);
430  func2->SetParameter(0,0.);
431  func2->SetParameter(1,0.);
432  func3->SetParameter(0,0.);
433  func3->SetParameter(1,0.);
434 
435  for (int l=0; l<20; l++) {
436  bkgd_shape_unnorm1[l] = (TH1F*)bkgd_hists1[i][j]->Clone();
437  bkgd_shape_unnorm2[l] = (TH1F*)bkgd_hists2[i][j]->Clone();
438  bkgd_shape_unnorm3[l] = (TH1F*)bkgd_hists3[i][j]->Clone();
439  signal_for_new1[l] = (TH1F*)signal_final2->Clone();
440  signal_for_new2[l] = (TH1F*)signal_final2->Clone();
441  signal_for_new3[l] = (TH1F*)signal_final2->Clone();
442 
443  // subtract off the Z contamination
444  bkgd_shape_unnorm1[l]->Add(zback_bkgd2,-1.);
445  bkgd_shape_unnorm2[l]->Add(zback_bkgd2,-1.);
446  bkgd_shape_unnorm3[l]->Add(zback_bkgd2,-1.);
447 
448  // calculate the signal in the normalization bins
449  signal_in_norm1[8] = func1->Integral(15,17);
450  signal_in_norm1[9] = func1->Integral(17,19);
451  signal_in_norm1[10] = func1->Integral(19,21);
452  signal_in_norm2[8] = func2->Integral(15,17);
453  signal_in_norm2[9] = func2->Integral(17,19);
454  signal_in_norm2[10] = func2->Integral(19,21);
455  signal_in_norm3[8] = func3->Integral(15,17);
456  signal_in_norm3[9] = func3->Integral(17,19);
457  signal_in_norm3[10] = func3->Integral(19,21);
458 
459  for (int m=8; m<=10; m++) {
460  if (signal_in_norm1[m] < 0) {signal_in_norm1[m] = 0.;}
461  if (signal_in_norm2[m] < 0) {signal_in_norm2[m] = 0.;}
462  if (signal_in_norm3[m] < 0) {signal_in_norm3[m] = 0.;}
463  }
464 
465  // calculate the normalization factor for 1 bin
466  float normt = 0, normb = 0.;
467  for (int k=8; k<=8; k++) {
468  if (bkgd_shape_unnorm1[l]->GetBinContent(k) > 0) {
469  normt += signal_final2->GetBinContent(k)-signal_in_norm1[k];
470  normb += bkgd_shape_unnorm1[l]->GetBinContent(k);
471  }
472  }
473  if (normb > 0 && normt > 0) {
474  float norm = normt/normb;
475  bkgd_shape_unnorm1[l]->Scale(norm);
476  }
477  for (int m=1; m<=49; m++) {
478  if (bkgd_shape_unnorm1[l]->GetBinContent(m) < 0) {bkgd_shape_unnorm1[l]->SetBinContent(m,0.);}
479  }
480  signal_for_new1[l]->Add(bkgd_shape_unnorm1[l],-1.);
481  signal_for_new1[l]->Fit(func1,"RQ");
482 
483  // calculate the normalization factor for 2 bins
484  normt = 0.; normb = 0.;
485  for (int k=8; k<=9; k++) {
486  if (bkgd_shape_unnorm2[l]->GetBinContent(k) > 0) {
487  normt += signal_final2->GetBinContent(k)-signal_in_norm2[k];
488  normb += bkgd_shape_unnorm2[l]->GetBinContent(k);
489  }
490  }
491  if (normb > 0 && normt > 0) {
492  float norm = normt/normb;
493  bkgd_shape_unnorm2[l]->Scale(norm);
494  //cout << "norm " << i << " " << l << " = " << norm << endl;
495  }
496 
497  for (int m=1; m<=49; m++) {
498  if (bkgd_shape_unnorm2[l]->GetBinContent(m) < 0) {bkgd_shape_unnorm2[l]->SetBinContent(m,0.);}
499  }
500  signal_for_new2[l]->Add(bkgd_shape_unnorm2[l],-1.);
501  signal_for_new2[l]->Fit(func2,"RQ");
502 
503  // calculate the normalization factor for 3 bins
504  normt = 0.; normb = 0.;
505  for (int k=8; k<=10; k++) {
506  if (bkgd_shape_unnorm3[l]->GetBinContent(k) > 0) {
507  normt += signal_final2->GetBinContent(k)-signal_in_norm3[k];
508  normb += bkgd_shape_unnorm3[l]->GetBinContent(k);
509  }
510  }
511  if (normb > 0 && normt > 0) {
512  float norm = normt/normb;
513  bkgd_shape_unnorm3[l]->Scale(norm);
514  }
515  for (int m=1; m<=49; m++) {
516  if (bkgd_shape_unnorm3[l]->GetBinContent(m) < 0) {bkgd_shape_unnorm3[l]->SetBinContent(m,0.);}
517  }
518  signal_for_new3[l]->Add(bkgd_shape_unnorm3[l],-1.);
519  signal_for_new3[l]->Fit(func3,"RQ");
520  } // end of for loop over l
521 
522  // Save the last iteration as the background histogram
523  new_bkgd_hists1[i][j] = (TH1F*)bkgd_shape_unnorm1[19]->Clone();
524  new_bkgd_hists2[i][j] = (TH1F*)bkgd_shape_unnorm2[19]->Clone();
525  new_bkgd_hists3[i][j] = (TH1F*)bkgd_shape_unnorm3[19]->Clone();
526 
527  //}
528  }
529 
530  // plot all the 1200 histograms
531  gStyle->SetTitleBorderSize(0);
532  TCanvas *can6 = new TCanvas("can6","can6",0,0,600,400);
533  signal_final2->SetStats(kFALSE);
534  if (charge == 1) {
535  signal_final2->SetTitle("W+ Background Shapes");
536  } else if (charge == -1) {
537  signal_final2->SetTitle("W- Background Shapes");
538  }
539  signal_final2->GetXaxis()->SetRangeUser(0.,70.);
540  signal_final2->GetXaxis()->SetTitle("2x2 Cluster E_{T} (GeV)");
541  signal_final2->Draw();
542  for (int i=0; i<=20; i++) {
543  int j=jval;
544  //for (int j=0; j<20; j++) {
545  new_bkgd_hists1[i][j]->Draw("same");
546  new_bkgd_hists2[i][j]->Draw("same");
547  new_bkgd_hists3[i][j]->Draw("same");
548  //}
549  }
550  new_bkgd->SetLineColor(4);
551  //new_bkgd->SetLineWidth(4.*new_bkgd->GetLineWidth());
552  new_bkgd->Draw("same");
553  if (charge == 1) {
554  //can6->Print("Wplus_bkgd_shapes.eps");
555  //can6->Print("Wplus_bkgd_shapes.png");
556  } else if (charge == -1) {
557  //can6->Print("Wminus_bkgd_shapes.eps");
558  //can6->Print("Wminus_bkgd_shapes.png");
559  }
560 
561  TH1F *chi2s = new TH1F("chi2s","chi2s",50,0.,10.);
562  for (int i=0; i<=20; i++) {
563  int j=jval;
564  //for (int j=0; j<20; j++) {
565  chi2s->Fill(final_chisquare[i][j]);
566  //}
567  }
568 
569  TCanvas *can7 = new TCanvas("can7","can7",0,0,600,400);
570  chi2s->Draw();
571  //can7->Print("chi2s.eps");
572  //can7->Print("chi2s.png");
573 
574  // ************************************************
575  // Now calculate the all background numbers and their
576  // systematic uncertainties (and spit them out to histograms)
577  // ************************************************
578 
579  // First get the simple numbers (backgrounds and their
580  // statistical uncertainties)
581  TH1F *tauhist = MC_dists_repack[2][0]->Clone();
582  //TH1F *zsig = MC_dists_repack[3][0]->Clone();
583  //TH1F *zeemc = MC_dists_repack[3][1]->Clone();
584  //TH1F *zback = MC_dists_repack[3][2]->Clone();
585  tauhist->Scale(taufrac);
586  float tau_norm = lumi_fact[2];
587  float Z_norm = lumi_fact[4];
588  float bkgd_sum = 0.;
589  float signal_sum = 0.;
590  float raw_sum = 0.;
591  float QCD_sum = 0., tau_sum = 0., eemc_sum = 0.;
592  float QCD_raw_sum = 0.;
593  float Wany_bkgd_sum = 0.;
594  float Wany_eemc_sum = 0.;
595  float zsig_sum = 0., zeemc_sum = 0.,zback_sum = 0.;
596  float wanysig_sum = 0., wanyeemc_sum = 0.,wanyback_sum = 0.;
597  for (int i=13; i<=26; i++) {
598  //cout << i << " " << new_bkgd->GetBinCenter(i) << " " << new_bkgd->GetBinContent(i) << " " << tauhist->GetBinContent(i) << " " << eemc_bkgd2->GetBinContent(i) << endl;
599  bkgd_sum += new_bkgd->GetBinContent(i);
600  //bkgd_sum += tauhist->GetBinContent(i);
601  bkgd_sum += eemc_bkgd2->GetBinContent(i);
602  QCD_raw_sum += bkgd_shape_nom->GetBinContent(i);
603  QCD_sum += new_bkgd->GetBinContent(i);
604  tau_sum += tauhist->GetBinContent(i);
605  eemc_sum += eemc_bkgd2->GetBinContent(i);
606  signal_sum += signal_final2->GetBinContent(i);
607  raw_sum += signal2->GetBinContent(i);
608  zsig_sum += zsig_bkgd2->GetBinContent(i);
609  zeemc_sum += zeemc_bkgd2->GetBinContent(i);
610  zback_sum += zback_bkgd2->GetBinContent(i);
611  wanysig_sum += wanysig_bkgd2->GetBinContent(i);
612  wanyeemc_sum += wanyeemc_bkgd2->GetBinContent(i);
613  wanyback_sum += wanyback_bkgd2->GetBinContent(i);
614 
615  }
616  cout << "The total background for ET>25 is " << bkgd_sum+zsig_sum << endl;
617  cout << "QCD = " << QCD_sum << ", tau = " << tau_sum << ", eemc = " << eemc_sum << ", and Z = " << zsig_sum << endl;
618  cout << "Raw = " << raw_sum << endl;
619  cout << "W Signal (w/o tau) = " << signal_sum-QCD_sum << endl;
620  cout << "Z in sig = " << zsig_sum << endl;
621  cout << "Z in eemc = " << zeemc_sum << endl;
622  cout << "Z in back = " << zback_sum << endl;
623  cout << "Wany in sig = " << wanysig_sum << endl;
624  cout << "Wany in eemc = " << wanyeemc_sum << endl;
625  cout << "Wany in back = " << wanyback_sum << endl;
626  cout << "QCD raw in back = " << QCD_raw_sum << endl;
627  cout << "The QCD stat unc. is " << norm*sqrt(QCD_sum/norm) << endl;
628  float tau_stat = tau_norm*taufrac*sqrt(tau_sum/(tau_norm*taufrac));
629  cout << "The tau stat unc. is " << tau_stat << endl;
630  float eemc_stat = sqrt(eemc_sum);
631  cout << "The eemc stat unc. is " << eemc_stat << endl;
632  float Z_stat = Z_norm*sqrt(zsig_sum/Z_norm);
633  cout << "The Z stat unc. is " << Z_stat << endl;
634  float tot_stat = sqrt(tau_stat*tau_stat+eemc_stat*eemc_stat+Z_stat*Z_stat);
635  cout << "tau+eemc+Z stat unc. is " << sqrt(tau_stat*tau_stat+eemc_stat*eemc_stat+Z_stat*Z_stat) << endl;
636 
637  cout << "f_tau = " << tau_sum/raw_sum << endl;
638  cout << "f_QCD = " << QCD_sum/raw_sum << endl;
639  cout << "f_EEMC = " << eemc_sum/raw_sum << endl;
640  cout << "f_Z = " << zsig_sum/raw_sum << endl;
641 
642  // Set up some histograms to hold all the errors that
643  // are calculated
644  TH1F *raw_stat_err2 = new TH1F("raw_stat_err2","raw_stat_err2",49,1.,99.);
645  TH1F *QCD_stat_err2 = new TH1F("QCD_stat_err2","QCD_stat_err2",49,1.,99.);
646  TH1F *eemc_stat_err2 = new TH1F("eemc_stat_err2","eemc_stat_err2",49,1.,99.);
647  TH1F *tau_stat_err2 = new TH1F("tau_stat_err2","tau_stat_err2",49,1.,99.);
648  TH1F *QCD_syst_high_err = new TH1F("QCD_syst_high_err","QCD_syst_high_err",49,1.,99.);
649  TH1F *QCD_syst_low_err = new TH1F("QCD_syst_low_err","QCD_syst_low_err",49,1.,99.);
650  TH1F *zsig_stat_err2 = new TH1F("zsig_stat_err2","zsig_stat_err2",49,1.,99.);
651  TH1F *zback_stat_err2 = new TH1F("zback_stat_err2","zback_stat_err2",49,1.,99.);
652 
653 
654  for (int i=1; i<=26; i++) {
655  raw_stat_err2->SetBinContent(i,signal2->GetBinContent(i));
656  QCD_stat_err2->SetBinContent(i,fabs(norm*new_bkgd->GetBinContent(i)));
657  eemc_stat_err2->SetBinContent(i,eemc_bkgd2->GetBinContent(i));
658  tau_stat_err2->SetBinContent(i,tau_norm*taufrac*tauhist->GetBinContent(i));
659  zsig_stat_err2->SetBinContent(i,Z_norm*zsig_bkgd2->GetBinContent(i));
660  zback_stat_err2->SetBinContent(i,norm*norm*Z_norm*zback_bkgd2->GetBinContent(i));
661  //cout << "Error " << i << " " << raw_stat_err2->GetBinCenter(i) << " " << raw_stat_err2->GetBinContent(i) << " " << QCD_stat_err2->GetBinContent(i) << " " << eemc_stat_err2->GetBinContent(i) << " " << tau_stat_err2->GetBinContent(i) << endl;
662  }
663 
664  // Now go through all the 1200 background shapes and find the
665  // high and low in each ET bin to give the maximum extent uncertainty
666  // for each ET bin (and also sum the low and high to give the
667  // overall number used in the xsec analysis)
668 
669  TH1F *low_bkgd = new_bkgd->Clone();
670  low_bkgd->SetName("low_bkgd");
671  TH1F *high_bkgd = new_bkgd->Clone();
672  high_bkgd->SetName("high_bkgd");
673 
674  float low_sum = 0.;
675  float high_sum = 0.;
676  for (int i=7; i<=49; i++) {
677  float high = 0.;
678  float low = 10000.;
679  for (int j=0; j<=20; j++) {
680  int k=jval;
681  //for (int k=0; k<20; k++) {
682  if (new_bkgd_hists1[j][k]->GetBinContent(i) < low) {
683  if (new_bkgd_hists1[j][k]->GetBinContent(i) >= 0) {
684  low = new_bkgd_hists1[j][k]->GetBinContent(i);
685  }
686  }
687  if (new_bkgd_hists1[j][k]->GetBinContent(i) > high) {
688  high = new_bkgd_hists1[j][k]->GetBinContent(i);
689  }
690 
691  if (new_bkgd_hists2[j][k]->GetBinContent(i) < low) {
692  if (new_bkgd_hists2[j][k]->GetBinContent(i) >= 0) {
693  low = new_bkgd_hists2[j][k]->GetBinContent(i);
694  }
695  }
696  if (new_bkgd_hists2[j][k]->GetBinContent(i) > high) {
697  high = new_bkgd_hists2[j][k]->GetBinContent(i);
698  }
699 
700  if (new_bkgd_hists3[j][k]->GetBinContent(i) < low) {
701  if (new_bkgd_hists3[j][k]->GetBinContent(i) >= 0) {
702  low = new_bkgd_hists3[j][k]->GetBinContent(i);
703  }
704  }
705  if (new_bkgd_hists3[j][k]->GetBinContent(i) > high) {
706  high = new_bkgd_hists3[j][k]->GetBinContent(i);
707  }
708  //cout << i << " low = " << low << " high = " << high << endl;
709  //} // end of k-loop
710  } // end of j-loop
711 
712  // calculate the sum
713  low_bkgd->SetBinContent(i,0.);
714  if ((low != 10000) && (new_bkgd->GetBinContent(i)-low > 0)) {
715  if ((i >= 13)&&(i<=26)) {low_sum += low;}
716  low_bkgd->SetBinContent(i,low);
717  }
718  if ((i >= 13)&&(i<=26)) {high_sum += high;}
719  high_bkgd->SetBinContent(i,high);
720  cout << i << " " << low_sum << " " << high_sum << endl;
721  //cout << QCD_syst_low_err->GetBinCenter(i) << " low = " << low << " high = " << high << " nom = " << new_bkgd->GetBinContent(i) << endl;
722  // set the bin-by-bin error too
723  if ((low != 10000) && (new_bkgd->GetBinContent(i)-low > 0)) {
724  QCD_syst_low_err->SetBinContent(i,new_bkgd->GetBinContent(i)-low);
725  } else {
726  QCD_syst_low_err->SetBinContent(i,0.);
727  }
728  QCD_syst_high_err->SetBinContent(i,high-new_bkgd->GetBinContent(i));
729 
730  } // end of i=loop
731 
732  cout << "QCD shape sys. unc. calc************" << endl;
733  cout << "The QCD low sum = " << low_sum << endl;
734  cout << "The QCD high sum = " << high_sum << endl;
735 
736  cout << "The QCD low error = " << QCD_sum-low_sum << endl;
737  cout << "The QCD high error = " << high_sum-QCD_sum << endl;
738 
739  cout << "Total low error = " << sqrt((QCD_sum-low_sum)*(QCD_sum-low_sum)+tot_stat*tot_stat) << endl;
740  cout << "Total high error = " << sqrt((high_sum-QCD_sum)*(high_sum-QCD_sum)+tot_stat*tot_stat) << endl;
741 
742  // ******************************************************
743  // Write out the histograms of interest to a
744  // file
745  // ******************************************************
746 
747  if (charge == 1) {
748  TFile *f2 = new TFile("bkgd_histos_pos_final.root","recreate");
749  } else if (charge == -1) {
750  TFile *f2 = new TFile("bkgd_histos_neg_final.root","recreate");
751  } else if (charge == 0) {
752  TFile *f2 = new TFile("bkgd_histos_sum.root","recreate");
753  }
754 
755  f2->cd();
756 
757  if (two_or_four == 2) {
758 
759  TH1F *signal_final3 = new TH1F("signal_final3","signal_final3",49,1.,99.);
760  for (int i=1; i<=49; i++) {
761  signal_final3->SetBinContent(i,signal_final2->GetBinContent(i));
762  }
763  signal_final3->Add(new_bkgd,-1.);
764 
765  tauhist->Write();
766  zsig->Write();
767  signal2->Write();
768  signal_final2->Write();
769  signal_final3->Write();
770  eemc_bkgd2->Write();
771  new_bkgd->Write();
772  low_bkgd->Write();
773  high_bkgd->Write();
774  signal_in_norm_region->Write();
775 
776  raw_stat_err2->Write();
777  QCD_stat_err2->Write();
778  eemc_stat_err2->Write();
779  tau_stat_err2->Write();
780  QCD_syst_high_err->Write();
781  QCD_syst_low_err->Write();
782  zsig_stat_err2->Write();
783  zback_stat_err2->Write();
784 
785  for (int i=1; i<=49; i++) {
786  cout.setf(ios::fixed);
787  cout.precision(2);
788  cout << " " << signal2->GetBinCenter(i) << " & " << signal2->GetBinContent(i) << " & " << signal_final3->GetBinContent(i) << " & " << QCD_stat_err2->GetBinContent(i) << " & " << eemc_stat_err2->GetBinContent(i) << " & " << tau_stat_err2->GetBinContent(i) << " & " << zsig_stat_err2->GetBinContent(i) << " & " << zback_stat_err2->GetBinContent(i) << " & " << QCD_syst_high_err->GetBinContent(i) << " & " << QCD_syst_low_err->GetBinContent(i) << " \\\\" << endl;
789 
790  }
791 
792  f2->Close();
793 
794  } else if (two_or_four == 4) {
795 
796  TH1F *signal_final3 = new TH1F("signal_final3","signal_final3",49,1.,99.);
797  for (int i=1; i<=49; i++) {
798  signal_final3->SetBinContent(i,signal_final2->GetBinContent(i));
799  }
800  signal_final3->Add(new_bkgd,-1.);
801 
802  TH1F *tauhist_repack = new TH1F("tauhist_r","tauhist_r",24,3.,99.);
803  TH1F *signal2_repack = new TH1F("signal2_r","signal2_r",24,3.,99.);
804  TH1F *signal_final2_repack = new TH1F("signal_final2_r","signal_final2_r",24,3.,99.);
805  TH1F *signal_final3_repack = new TH1F("signal_final3_r","signal_final3_r",24,3.,99.);
806  TH1F *eemc_bkgd2_repack = new TH1F("eemc_bkgd2_r","eemc_bkgd2_r",24,3.,99.);
807  TH1F *new_bkgd_repack = new TH1F("new_bkgd_r","new_bkgd_r",24,3.,99.);
808  TH1F *low_bkgd_repack = new TH1F("low_bkgd_r","low_bkgd_r",24,3.,99.);
809  TH1F *high_bkgd_repack = new TH1F("high_bkgd_r","high_bkgd_r",24,3.,99.);
810  TH1F *signal_in_norm_region_repack = new TH1F("signal_in_norm_region_r","signal_in_norm_region_r",24,3.,99.);
811 
812  TH1F *raw_stat_err2_repack = new TH1F("raw_stat_err2_r","raw_stat_err2_r",24,3.,99.);
813  TH1F *QCD_stat_err2_repack = new TH1F("QCD_stat_err2_r","QCD_stat_err2_r",24,3.,99.);
814  TH1F *eemc_stat_err2_repack = new TH1F("eemc_stat_err2_r","eemc_stat_err2_r",24,3.,99.);
815  TH1F *tau_stat_err2_repack = new TH1F("tau_stat_err2_r","tau_stat_err2_r",24,3.,99.);
816  TH1F *QCD_syst_high_err_repack = new TH1F("QCD_syst_high_err_r","QCD_syst_high_err_r",24,3.,99.);
817  TH1F *QCD_syst_low_err_repack = new TH1F("QCD_syst_low_err_r","QCD_syst_low_err_r",24,3.,99.);
818  TH1F *zsig_stat_err2_repack = new TH1F("zsig_stat_err2_r","zsig_stat_err2_r",24,3.,99.);
819  TH1F *zback_stat_err2_repack = new TH1F("zback_stat_err2_r","zback_stat_err2_r",24,3.,99.);
820 
821  for (int i=1; i<=24; i++) {
822  tauhist_repack->SetBinContent(i,tauhist->GetBinContent(2*i)+
823  tauhist->GetBinContent(2*i+1));;
824  signal2_repack->SetBinContent(i,signal2->GetBinContent(2*i)+
825  signal2->GetBinContent(2*i+1));
826  signal_final2_repack->SetBinContent(i,signal_final2->GetBinContent(2*i)+
827  signal_final2->GetBinContent(2*i+1));
828  signal_final3_repack->SetBinContent(i,signal_final3->GetBinContent(2*i)+
829  signal_final3->GetBinContent(2*i+1));
830  eemc_bkgd2_repack->SetBinContent(i,eemc_bkgd2->GetBinContent(2*i)+
831  eemc_bkgd2->GetBinContent(2*i+1));
832  new_bkgd_repack->SetBinContent(i,new_bkgd->GetBinContent(2*i)+
833  new_bkgd->GetBinContent(2*i+1));
834  low_bkgd_repack->SetBinContent(i,low_bkgd->GetBinContent(2*i)+
835  low_bkgd->GetBinContent(2*i+1));
836  high_bkgd_repack->SetBinContent(i,high_bkgd->GetBinContent(2*i)+
837  high_bkgd->GetBinContent(2*i+1));
838  signal_in_norm_region_repack->SetBinContent(i,signal_in_norm_region->GetBinContent(2*i)+signal_in_norm_region->GetBinContent(2*i+1));
839 
840  raw_stat_err2_repack->SetBinContent(i,raw_stat_err2->GetBinContent(2*i)+
841  raw_stat_err2->GetBinContent(2*i+1));
842  QCD_stat_err2_repack->SetBinContent(i,QCD_stat_err2->GetBinContent(2*i)+
843  QCD_stat_err2->GetBinContent(2*i+1));
844  eemc_stat_err2_repack->SetBinContent(i,eemc_stat_err2->GetBinContent(2*i)+
845  eemc_stat_err2->GetBinContent(2*i+1));
846  tau_stat_err2_repack->SetBinContent(i,tau_stat_err2->GetBinContent(2*i)+
847  tau_stat_err2->GetBinContent(2*i+1));
848  QCD_syst_high_err_repack->SetBinContent(i,QCD_syst_high_err->GetBinContent(2*i)+QCD_syst_high_err->GetBinContent(2*i+1));
849  QCD_syst_low_err_repack->SetBinContent(i,QCD_syst_low_err->GetBinContent(2*i)+QCD_syst_low_err->GetBinContent(2*i+1));
850  zsig_stat_err2_repack->SetBinContent(i,zsig_stat_err2->GetBinContent(2*i)+
851  zback_stat_err2->GetBinContent(2*i+1));
852  zback_stat_err2_repack->SetBinContent(i,zback_stat_err2->GetBinContent(2*i)+
853  zback_stat_err2->GetBinContent(2*i+1));
854  }
855 
856  tauhist_repack->Write();
857  signal2_repack->Write();
858  signal_final2_repack->Write();
859  signal_final3_repack->Write();
860  eemc_bkgd2_repack->Write();
861  new_bkgd_repack->Write();
862  low_bkgd_repack->Write();
863  high_bkgd_repack->Write();
864  signal_in_norm_region_repack->Write();
865 
866  raw_stat_err2_repack->Write();
867  QCD_stat_err2_repack->Write();
868  eemc_stat_err2_repack->Write();
869  tau_stat_err2_repack->Write();
870  QCD_syst_high_err_repack->Write();
871  QCD_syst_low_err_repack->Write();
872  zsig_stat_err2_repack->Write();
873  zback_stat_err2_repack->Write();
874 
875  for (int i=1; i<=24; i++) {
876  cout.setf(ios::fixed);
877  cout.precision(2);
878  cout << " " << signal2_repack->GetBinCenter(i) << " & " << signal2_repack->GetBinContent(i) << " & " << signal_final3_repack->GetBinContent(i) << " & " << QCD_stat_err2_repack->GetBinContent(i) << " & " << eemc_stat_err2_repack->GetBinContent(i) << " & " << tau_stat_err2_repack->GetBinContent(i) << " & " << zsig_stat_err2_repack->GetBinContent(i) << " & " << zback_stat_err2_repack->GetBinContent(i) << " & " << QCD_syst_high_err_repack->GetBinContent(i) << " & " << QCD_syst_low_err_repack->GetBinContent(i) << " \\\\" << endl;
879  }
880 
881  f2->Close();
882 
883  }
884 
885 
886 } // end of macro
Generated by   doxygen 1.8.5