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