StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
plWJacob.C
1 int blueCol=4;
2 int yellCol=kYellow;
3 int stage=3;
4 float yMax=175;
5 float etMax=70;
6 
7 void plWJacob(int charge=0) {
8  gStyle->SetPalette(1,0);
9  gStyle->SetOptStat(1000000);
10 
11  TString iPath="/star/data05/scratch/stevens4/wAnalysisOut/apsXsec/data/";
12  TString core="run9setABCD";
13 
14  TFile *fd = TFile::Open(iPath+core+".wana.hist.root");
15 
16  assert(fd->IsOpen());
17 
18  if(charge==0){
19  //charge summed histos
20  hwE=(TH1F*)fd->Get("muclustPtBal"); assert(hwE);
21  hnE=(TH1F*)fd->Get("muclustPtBalnoE"); assert(hnE);
22  hQCD=(TH1F*)fd->Get("muclustPtBal_bckgrd"); assert(hQCD);
23  }
24  else if(charge==1){
25  //charge summed histos
26  hwE=(TH1F*)fd->Get("pubclustPtBalP"); assert(hwE);
27  hnE=(TH1F*)fd->Get("pubclustPtBalnoEP"); assert(hnE);
28  hQCD=(TH1F*)fd->Get("pubclustPtBal_bckgrdP"); assert(hQCD);
29  }
30  else if(charge==-1){
31  //charge summed histos
32  hwE=(TH1F*)fd->Get("pubclustPtBalN"); assert(hwE);
33  hnE=(TH1F*)fd->Get("pubclustPtBalnoEN"); assert(hnE);
34  hQCD=(TH1F*)fd->Get("pubclustPtBal_bckgrdN"); assert(hQCD);
35  }
36 
37  // plot raw spectra ..........
38  c=new TCanvas();
39  c->Divide(2,3);
40  c->cd(2);
41  hwE->Draw();
42  c->cd(4);
43  hnE->Draw(); // hnE->SetAxisRange(0.,etMax);
44 
45  //.. compute Intermediate results
46  TAxis* axX=hwE->GetXaxis();
47  float x1=axX->GetXmin()+1.0;
48  float x2=axX->GetXmax()+1.0;
49  int nbX=axX->GetNbins();
50 
51  // assumption: output histos will have 100 bins [1,101] GeV
52 
53  // those 3 have errors =sqrt(N)
54  hA=new TH1F("hA","W candidates (input)",nbX,x1,x2);
55  hAA=new TH1F("hAA","W candidates- 2nd Endcap",nbX,x1,x2);
56 
57  hB=new TH1F("hB","rejected by including Endcap",nbX,x1,x2);
58  hC=new TH1F("hC","rejected by inverting cut",nbX,x1,x2);
59 
60  // .......copy & subtract input spectra, subtract 2nd endcap
61  float sumA=0,sumB=0,sumC=0; // over ET<20
62  for(int i=1;i<=nbX-1;i++) {
63  float y1=hwE->GetBinContent(i+1);
64  float y2=hnE->GetBinContent(i+1);
65  float y3=hQCD->GetBinContent(i+1);
66  // printf("%d , %.1f, %.0f, %.0f, %.0f \n",i,hA->GetBinCenter(i),y1,y2,y3);
67  // compute yields & comulative variance
68  float yA=y1;
69  float yB=y2-y1; // Endcap background, no adding errors in quadrature
70  float yC=y3; // rejected by away ET
71  hA->SetBinContent(i,yA); // W candidates
72  hB->SetBinContent(i,yB); // endcapBCG
73  hC->SetBinContent(i,yC); // away BCG
74  }
75 
76  hC->Sumw2(); hAA->Sumw2();
77  // subtract 2nd endcap
78  hAA->Add(hA,1.); hAA->Add(hB,-1.);
79 
80  c->cd(1); hA->Draw();// hB->SetFillColor(8);
81 
82  c->cd(3); hB->Draw(); hB->SetFillColor(8);
83 
84 
85  c->cd(5); hC->Draw("he");
86  hC->SetLineColor(kRed);
87  hC->SetLineWidth(2);
88 
89  hA->Rebin(); hAA->Rebin(); hB->Rebin(); hC->Rebin();
90  nbX/=2;
91  hA->SetLineWidth(2.); hA->SetAxisRange(0.,etMax);
92  hB->SetAxisRange(0.,etMax);
93  hC->SetAxisRange(0.,etMax);
94 
95  float sumA=0, sumC=0; // over ET =[15,21] GeV, delET=6 GeV
96  for(int i=1;i<=nbX;i++) {// decides on backg norm. range
97  float y1=hAA->GetBinContent(i);
98  float y2=hC->GetBinContent(i);
99  // printf("%d %f %f \n",i,y1,y2);
100  if(i>=8 && i<=9) { sumA+=y1; sumC+=y2;}
101  }
102  float fact=sumA/sumC;
103  printf("sumA=%.1f sumC=%.1f fac=%.3f\n",sumA, sumC, fact);
104  hC->Scale(fact);
105 
106  //..........................................
107  //..........................................
108  // ..... "Preliminary" Jacobian peak
109  //..........................................
110  //..........................................
111  TString tt="W2009-prelimB"; tt+=stage;
112 
113  c=new TCanvas(tt,tt,600,600);
114  gStyle->SetOptStat(0);
115  c->SetFillColor(kWhite);
116  c->SetLeftMargin(0.15);
117  c->SetFrameBorderMode(0);
118 
119 
120 
121  // input W spectrum
122  hA3=(TH1F*) hA->Clone();
123  hA3->SetTitle("PT Balance Cut; EMC cluster #font[72]{E_{T}} (GeV) ;Counts");
124  hA3->GetYaxis()->SetTitleOffset(1.5);
125  hA3->SetAxisRange(7,69); hA3->SetMinimum(-39); hA3->SetMaximum(yMax);
126  hA3->Draw("h"); hA3->SetLineWidth(2);
127  ln=new TLine(8,0,70,0); ln->Draw();
128  printf("raw sum=%.0f\n",hA3->Integral());
129  hA3->SetTitle(iPath+core);
130 
131  // pure Ws
132  hDD=(TH1F*) hAA->Clone();
133  hDD->Add(hC,-1.); // subtracted scaled away ET
134  hDD->SetMarkerStyle(8);
135  hDD2=(TH1F*) hDD->Clone();
136  hDD2->SetMarkerStyle(8); hDD2->SetMarkerSize(1);
137  hDD->SetFillColor(yellCol); //hDD2->SetLineWidth(2);
138 
139  if( stage>2) { hDD->Draw("same h"); hDD2->Draw("same h"); }
140 
141  // away backg+ Endcap backg
142  hCB=(TH1F*) hC->Clone();
143  hCB->Add(hB,1.);
144  hCB->SetMarkerStyle(8); hCB->SetMarkerSize(1); hCB->SetMarkerColor(blueCol);
145  hCB->SetLineColor(blueCol); hCB->SetLineWidth(2);
146 
147  hCB->SetAxisRange(14.5,65);
148  if(stage>1) hCB->Draw("p e h same");
149 
150  // draw inut again w/ dashed style
151  hA3b=(TH1F*) hA3->Clone();
152  hA3b->SetLineStyle(2);
153  hA3b->Draw("same"); hA3b->SetAxisRange(15,18);
154 
155 
156  lg=new TLegend(0.32,0.64,0.66,0.80); // top
157  TDatime dt;lg->SetHeader(dt.AsString());
158  //lg->SetHeader(" cluster |#eta|_{ }<_{ }1");
159  if(charge==0)
160  lg->AddEntry(hA,"#font[52]{W} candidates (charge summed)","l");
161  else if(charge==1)
162  lg->AddEntry(hA,"#font[52]{W+} candidates ","l");
163  else
164  lg->AddEntry(hA,"#font[52]{W-} candidates ","l");
165  lg->AddEntry(hCB,"QCD backg. est.","lfep");
166  lg->AddEntry(hDD,"Backg. subtr. #font[52]{W}s","lfep");
167  lg->SetFillColor(kWhite);
168  lg->SetLineColor(kWhite);
169  lg->Draw();
170 
171  TLatex *lat1 = new TLatex(0.32,0.84,"STAR #font[72]{#vec{p}_{ }+_{ }#vec{p}} #sqrt{#font[72]{s}}_{ }=_{ }500 GeV");
172  lat1->SetNDC(); lat1->Draw("same");
173 
174 
175  tx=new TText(14,20,"software threshold"); tx->SetTextAngle(90);
176  tx->SetTextSize(0.03);//tx->SetTextColor(kMagenta);
177  tx->Draw();
178 
179  //....... Joe's additional variables
180  // Calcate the total relative error squared for the cross section
181  // which is roughly (S+2B)/S^2
182  // and for the AL which is roughly (S+B?)/S^2
183 
184  float tot_xsec_err_num = 0.;
185  float tot_xsec_err_den = 0.;
186  float tot_al_err_num = 0.;
187  float tot_al_err_den = 0.;
188  for (int i=11; i<=hA->GetNbinsX(); i++) {// defines ET range for x-section
189  tot_xsec_err_num += hA->GetBinError(i)*hA->GetBinError(i);
190  tot_xsec_err_num += hB->GetBinError(i)*hB->GetBinError(i);
191  tot_xsec_err_num += hC->GetBinError(i)*hC->GetBinError(i);
192  tot_xsec_err_den += hDD->GetBinContent(i);
193  // cout << i << " " << hA->GetBinContent(i) << endl;
194  }
195 
196  for (int i=13; i<=hA->GetNbinsX(); i++) { // defines ET range for AL
197  tot_al_err_num += hA->GetBinError(i)*hA->GetBinError(i);
198  tot_al_err_den += hDD->GetBinContent(i);
199  }
200 
201  float tot_xsec_rel_err = sqrt(tot_xsec_err_num)/tot_xsec_err_den;
202  float tot_al_err = sqrt(tot_al_err_num)/tot_al_err_den;
203  printf(" Joe measures: xsec_rel_err=%.3f AL_err=%.3f\n", tot_xsec_rel_err, tot_al_err);
204 
205  char text[100];
206  sprintf(text,"Error estim: xSec_rel=%.4f AL=%.3f", tot_xsec_rel_err, tot_al_err);
207  tx=new TText(21,-18,text); tx->Draw(); tx->SetTextSize(0.03);
208 
209  //...... end Joe' calculation
210 
211  return; // not print
212  for(int i=1;i<=nbX;i++) {
213  printf("%d , %.1f,",i,hA->GetBinCenter(i));
214  printf(" %.1f, %.1f, ", hA->GetBinContent(i), hA->GetBinError(i));
215  printf(" %.1f, %.1f, ", hCB->GetBinContent(i), hCB->GetBinError(i));
216  printf(" %.1f, %.1f\n ", hDD->GetBinContent(i), hDD->GetBinError(i));
217  }
218 
219 }
Generated by   doxygen 1.8.5