StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Calall.C
1 //this macro is designed to calculate the single spin asymmetries from yellow and blue beams, and double spin asymmetry from pi0 study, based on the generated/fitted/normalized pi0 yield from FitPt2Mass.C in our fitting philosopy.
2 //author: Weihong He
3 
4 #include <iostream>
5 #include <string>
6 #include <fstream>
7  using namespace std;
8 //X:yellow Y:blue
9 //background_LL is calculated from BackgroundLL.C because we want to use another method: side-band analysis to background study, although you can also calculate background all here by simply replacing the all array to back array.
10 double epsilon_ll[7],ep_error[7],ptmean[7],epsilon_real[7],epsilon_back[7],ds[7],dr[7],dw[7],db[7],epsilon_X[7],epsilon_Y[7],dx[7],dy[7],ep_error_rb[7];
11 Calall(){
12  //x*[*] save epsilon from yellow beam, y*[*] save epsilon from blue beam, all*[*] save all values.
13  double x1[7],x2[7],x3[7],y1[7],y2[7],y3[7],all1[7],all2[7],all3[7];
14  //corresponding errors
15  double ex1[7],ex2[7],ex3[7],ey1[7],ey2[7],ey3[7],eall1[7],eall2[7],eall3[7];
16  double pt1[7],pt2[7],pt3[7];
17  gStyle->SetPadGridX(0);
18  gStyle->SetPadGridY(0);
19  gStyle->SetCanvasColor(0);
20  gStyle->SetOptStat(0);
21  c0=new TCanvas("epsilony","epsilony",700,500);
22  c0->Divide(1,1);
23  c1=new TCanvas("epsilonx","epsilonx",700,500);
24  c1->Divide(1,1);
25  c2=new TCanvas("all","all",700,500);
26  c2->Divide(1,1);
27  TH1F* h1=new TH1F("epsilon_LB","#epsilon_{LB} vs pT",60,4.,14.);
28  TH1F* h2=new TH1F("epsilon_LY","#epsilon_{LY} vs pT",60,4.,14.);
29  //theory curve histograms
30  TH1F* h5=new TH1F("A_LL","#vec{p}+#vec{p}->#pi^{0}X, #sqrt{s}=200Gev, 1.0#leq#eta#leq2.0;pT[Gev/c];A_{LL}",22,1.0,12.0);
31  TH1F* h6=new TH1F("A_LL","#vec{p}+#vec{p}->#pi^{0}X, #sqrt{s}=200Gev, 1.089#leq#eta#leq2.0,P=60%;[Gev/c];A_{LL}",22,1.0,12.0);
32  TH1F* h7=new TH1F("A_LL","#vec{p}+#vec{p}->#pi^{0}X, #sqrt{s}=200Gev, 1.089#leq#eta#leq2.0,P=60%;[Gev/c];A_{LL}",22,1.0,12.0);
33  TH1F* h8=new TH1F("A_LL","#vec{p}+#vec{p}->#pi^{0}X, #sqrt{s}=200Gev, 1.089#leq#eta#leq2.0,P=60%;[Gev/c];A_{LL}",22,1.0,12.0);
34 
35 
36  gStyle->SetOptStat(0);
37 
38  //dataset 1
39  ifstream yield1("yield3335557.txt");
40  doCalall(yield1);
41  for(int j=0;j<7;j++){
42  x1[j]=epsilon_X[j];
43  ex1[j]=dx[j];
44  y1[j]=epsilon_Y[j];
45  ey1[j]=dy[j];
46  all1[j]=epsilon_ll[j];
47  eall1[j]=ep_error_rb[j];
48  pt1[j]=ptmean[j];
49 
50  }
51 
52  //dataset 2
53  ifstream yield2("yield4444446.txt");
54  doCalall(yield2);
55  for(int j=0;j<7;j++){
56  x2[j]=epsilon_X[j];
57  ex2[j]=dx[j];
58  y2[j]=epsilon_Y[j];
59  ey2[j]=dy[j];
60  all2[j]=epsilon_ll[j];
61  eall2[j]=ep_error_rb[j];
62  pt2[j]=ptmean[j];
63 
64  }
65 
66  //dataset3
67  ifstream yield3("yield5556668.txt");
68  doCalall(yield3);
69  for(int j=0;j<7;j++){
70  x3[j]=epsilon_X[j];
71  ex3[j]=dx[j];
72  y3[j]=epsilon_Y[j];
73  ey3[j]=dy[j];
74  all3[j]=epsilon_ll[j];
75  eall3[j]=ep_error_rb[j];
76  pt3[j]=ptmean[j];
77 
78  }
79  double big[7],small[7];
80  for(int k=0;k<7;k++){
81  if(all1[k]>all2[k]) big[k]=all1[k];
82  else big[k]=all2[k];
83  if(all3[k]>big[k]) big[k]=all3[k];
84 
85  if(all1[k]<all2[k]) small[k]=all1[k];
86  else small[k]=all2[k];
87  if(all3[k]<small[k]) small[k]=all3[k];
88  }
89  double sys_error[7];
90  for(int m=0;m<7;m++){
91  sys_error[m]=big[m]-small[m];
92  //cout<<"syserror="<<sys_error[m]<<endl;
93  }
94  //calculate All
95 
96  const int nBins = 7;
97  float bins[nBins+1] = {4,5,6,7,8,9,10,12};
98 #if 1
99  TH1F* h4 = new TH1F("A_LL","#vec{p}+#vec{p}->#pi^{0}X, #sqrt{s}=200Gev, 1.0#leq#eta#leq2.0;pT[Gev/c];A_{LL};",nBins,bins);
100  h4->SetMarkerStyle(8);
101 #endif
102  TH1F* hSys = new TH1F("A_LL","#vec{p}+#vec{p}->#pi^{0}X, #sqrt{s}=200Gev, 1.0#leq#eta#leq2.0;pT[Gev/c];A_{LL};",nBins,bins);
103  TH1F* hB = new TH1F("A_LL","#vec{p}+#vec{p}->#pi^{0}X, #sqrt{s}=200Gev, 1.0#leq#eta#leq2.0;pT[Gev/c];A_{LL};",nBins,bins);
104  float baseLine=-0.17;
105  float sysErr[7] = {0.00239,0.00521,0.00485,0.00405,0.012,0.0108,0.0104};
106 
107  //point 1
108  h4->Fill(pt3[0],all3[0]);
109  int nbin0=h4->FindBin(pt3[0]);
110  h4->SetBinError(nbin0,eall3[0]);//+sys_error[0]);//+0.0005315);
111  h1->Fill(pt3[0],y3[0]);
112  int nbbin0=h1->FindBin(pt3[0]);
113  h1->SetBinError(nbbin0,ey3[0]);
114  h2->Fill(pt3[0],x3[0]);
115  int nbinn0=h2->FindBin(pt3[0]);
116  h2->SetBinError(nbinn0,ex3[0]);
117  cout<<"point1 error="<<eall3[0]<<" all="<<all3[0]<<endl;
118  hSys->SetBinContent(nbin0,baseLine + sysErr[0]);
119  hB->SetBinContent(nbin0,baseLine);
120  //hSys->SetBinError(nbin0,sysErr[0]/2.);
121  //point 2
122  h4->Fill(pt2[1],all2[1]);
123  int nbin1=h4->FindBin(pt2[1]);
124  h4->SetBinError(nbin1,eall2[1]);//+sys_error[1]);//+0.0002291);
125  h1->Fill(pt2[1],y2[1]);
126  int nbbin1=h1->FindBin(pt2[1]);
127  h1->SetBinError(nbbin1,ey2[1]);
128  h2->Fill(pt2[1],x2[1]);
129  int nbinn1=h2->FindBin(pt2[1]);
130  h2->SetBinError(nbinn1,ex2[1]);
131  cout<<"point2 error="<<eall2[1]<<" all="<<all2[1]<<endl;
132  hSys->SetBinContent(nbin1,baseLine + sysErr[1]);
133  //hSys->SetBinError(nbin1,sysErr[1]/2.);
134  hB->SetBinContent(nbin1,baseLine);
135  //point 3
136  h4->Fill(pt1[2],all1[2]);
137  int nbin2=h4->FindBin(pt1[2]);
138  h4->SetBinError(nbin2,eall1[2]);//+sys_error[2]);//+0.0005046);
139  h1->Fill(pt1[2],y1[2]);
140  int nbbin2=h1->FindBin(pt1[2]);
141  h1->SetBinError(nbbin2,ey1[2]);
142  h2->Fill(pt1[2],x1[2]);
143  int nbinn2=h2->FindBin(pt1[2]);
144  h2->SetBinError(nbinn2,ex1[2]);
145  cout<<"point3 error="<<eall1[2]<<" all="<<all1[2]<<endl;
146  hSys->SetBinContent(nbin2,baseLine + sysErr[2]);
147  //hSys->SetBinError(nbin2,sysErr[2]/2.);
148  hB->SetBinContent(nbin2,baseLine);
149  //point 4
150  h4->Fill(pt2[3],all2[3]);
151  int nbin3=h4->FindBin(pt2[3]);
152  h4->SetBinError(nbin3,eall2[3]);//+sys_error[3]);//+0.0004908);
153  h1->Fill(pt2[3],y2[3]);
154  int nbbin3=h1->FindBin(pt2[3]);
155  h1->SetBinError(nbbin3,ey2[3]);
156  h2->Fill(pt2[3],x2[3]);
157  int nbinn3=h2->FindBin(pt2[3]);
158  h2->SetBinError(nbinn3,ex2[3]);
159  cout<<"point4 error="<<eall2[3]<<" all="<<all2[3]<<endl;
160  hSys->SetBinContent(nbin3,baseLine + sysErr[3]);
161  //hSys->SetBinError(nbin3,sysErr[3]/2.);
162  hB->SetBinContent(nbin3,baseLine);
163  //point 5
164  h4->Fill(pt2[4],all2[4]);
165  int nbin4=h4->FindBin(pt2[4]);
166  h4->SetBinError(nbin4,eall2[4]);//+sys_error[4]);//+0.0002051);
167  h1->Fill(pt2[4],y2[4]);
168  int nbbin4=h1->FindBin(pt2[4]);
169  h1->SetBinError(nbbin4,ey2[4]);
170  h2->Fill(pt2[4],x2[4]);
171  int nbinn4=h2->FindBin(pt2[4]);
172  h2->SetBinError(nbinn4,ex2[4]);
173  cout<<"point5 error="<<eall2[4]<<" all="<<all2[4]<<endl;
174  hSys->SetBinContent(nbin4,baseLine + sysErr[4]);
175  //hSys->SetBinError(nbin4,sysErr[4]/2.);
176  hB->SetBinContent(nbin4,baseLine);
177  //point 6
178  h4->Fill(pt1[5],all1[5]);
179  int nbin5=h4->FindBin(pt1[5]);
180  h4->SetBinError(nbin5,eall1[5]);//+sys_error[5]);//+0.00055);
181  h1->Fill(pt1[5],y1[5]);
182  int nbbin5=h1->FindBin(pt1[5]);
183  h1->SetBinError(nbbin5,ey1[5]);
184  h2->Fill(pt1[5],x1[5]);
185  int nbinn5=h2->FindBin(pt1[5]);
186  h2->SetBinError(nbbin5,ex1[5]);
187  cout<<"point6 error="<<eall1[5]<<" all="<<all1[5]<<endl;
188  hSys->SetBinContent(nbin5,baseLine + sysErr[5]);
189  //hSys->SetBinError(nbin5,sysErr[5]/2.);
190  hB->SetBinContent(nbin5,baseLine);
191  //point 7
192  h4->Fill(pt2[6],all2[6]);
193  int nbin6=h4->FindBin(pt2[6]);
194  h4->SetBinError(nbin6,eall2[6]);//+sys_error[6]);//+0.001953);
195  h1->Fill(pt2[6],y2[6]);
196  int nbbin6=h1->FindBin(pt2[6]);
197  h1->SetBinError(nbbin6,ey2[6]);
198  h2->Fill(pt2[6],x2[6]);
199  int nbinn6=h2->FindBin(pt2[6]);
200  h2->SetBinError(nbinn6,ex2[6]);
201  cout<<"point7 error="<<eall2[6]<<" all="<<all2[6]<<endl;
202  hSys->SetBinContent(nbin6,baseLine + sysErr[6]);
203  //hSys->SetBinError(nbin6,sysErr[6]/2.);
204  hB->SetBinContent(nbin6,baseLine);
205  //open theory curve
206  ifstream unp("theorycurve/pion-unp-cteq6-rap1to2.dat");
207  ifstream g0("theorycurve/pion-pol-g0-rap1to2.dat");
208  ifstream plusg("theorycurve/pion-pol-max-rap1to2.dat");
209  ifstream minusg("theorycurve/pion-pol-maxminus-rap1to2.dat");
210  ifstream stdg("theorycurve/pion-pol-std-rap1to2.dat");
211  double tx1=0.0,tx2=0.0,tx3=0.0,tx4=0.0,tx5=0.0;
212  double td1=0.0,td2=0.0,td3=0.0,td4=0.0,td5=0.0;
213  double tb1=0.0,tb2=0.0,tb3=0.0,tb4=0.0,tb5=0.0;
214  double tc1=0.0,tc2=0.0,tc3=0.0,tc4=0.0,tc5=0.0;
215  double ty1=0.0,ty2=0.0,ty3=0.0,ty4=0.0,ty5=0.0;
216  int ccount=0;
217  double theoryd[18],theoryb[18],theoryc[18],theoryg[18];;
218  while(1){
219  unp>>ty1>>ty2>>ty3>>ty4>>ty5;
220  g0>>tx1>>tx2>>tx3>>tx4>>tx5;
221  plusg>>td1>>td2>>td3>>td4>>td5;
222  minusg>>tb1>>tb2>>tb3>>tb4>>tb5;
223  stdg>>tc1>>tc2>>tc3>>tc4>>tc5;
224  theoryd[ccount]=td5/ty5;
225  theoryb[ccount]=tb5/ty5;
226  theoryc[ccount]=tc5/ty5;
227  theoryg[ccount]=tx5/ty5;
228  //cout<<"x3="<<d3<<" theoryd="<<theoryg[ccount]<<endl;
229  h5->Fill(td3,theoryd[ccount]);
230  h6->Fill(tb3,theoryb[ccount]);
231  h7->Fill(tc3,theoryc[ccount]);
232  h8->Fill(tx3,theoryg[ccount]);
233  ccount++;
234  if(ccount>17) break;
235  }
236  h5->Fill(8.0,0.03);
237  h5->SetAxisRange(0,25,"X");
238  h5->Fill(9.0,0.0305);
239  h5->Fill(10.0,0.031);
240  h5->Fill(10.5,0.0312);
241  h5->Fill(11.5,0.031);
242  h6->Fill(8.0,-0.001);
243  h6->Fill(9.0,-0.0035);
244  h6->Fill(10.0,-0.0065);
245  h6->Fill(10.5,-0.0075);
246  h6->Fill(11.5,-0.01);
247  h7->Fill(8.0,0.0072);
248  h7->Fill(9.0,0.0079);
249  h7->Fill(10.0,0.0087);
250  h7->Fill(10.5,0.0091);
251  h7->Fill(11.5,0.0099);
252  h8->Fill(8.0,0.00152);
253  h8->Fill(9.0,0.00183);
254  h8->Fill(10.0,0.00204);
255  h8->Fill(10.5,0.0022);
256  h8->Fill(11.5,0.0026);
257  //plotting epsilon_ll, h1 for blue beam, h2 for yellow beam
258 #if 1
259  c0->cd(1);
260  h1->Draw();
261  h1->Fit("pol0","R","",4.0,12.0);
262  h1->SetMinimum(-0.05);
263  h1->SetMaximum(0.05);
264  TLegend* leg1 = new TLegend(0.15,0.65,0.5,0.85);
265  leg1->SetHeader("STAR 2006 Preliminary #pi^{0}");
266  leg1->SetBorderSize(0.0);
267  leg1->SetFillColor(0.0);
268  leg1->Draw();
269  c0->Print("epsilony.gif");
270  c1->cd(1);
271  h2->Draw();
272  h2->Fit("pol0","R","",4.0,12.0);
273  h2->SetMinimum(-0.05);
274  h2->SetMaximum(0.05);
275  TLegend* leg2 = new TLegend(0.15,0.65,0.5,0.85);
276  leg2->SetHeader("STAR 2006 Preliminary #pi^{0}");
277  leg2->SetBorderSize(0.0);
278  leg2->SetFillColor(0.0);
279  leg2->Draw();
280  c1->Print("epsilonx.gif");
281 #endif
282 
283 #if 1
284  //A_LL plotting
285  c2->cd(1);
286  h5->SetLineColor(kRed);
287  h5->Draw("C");
288  h5->SetMinimum(-0.18);
289  h5->SetMaximum(0.15);
290  hB->SetFillColor(11);
291  hB->Draw("same");
292  hSys->SetFillColor(10);
293  hSys->Draw("same");
294  h6->SetLineColor(kGreen);
295  h6->Draw("Csame");
296  h7->Draw("Csame");
297  h8->SetLineColor(kBlue);
298  h8->Draw("Csame");
299 
300  h4->Draw("E1same,p");
301 
302  TLegend* leg = new TLegend(0.15,0.6,0.4,0.8);
303  leg->AddEntry(h5,"#Delta{G}=G","L");
304  leg->AddEntry(h6,"#Delta{G}=-G","L");
305  leg->AddEntry(h7,"#Delta{G}=std","L");
306  leg->AddEntry(h8,"#Delta{G}=0","L");
307  leg->SetHeader("Vogelsang Prediction(GRSV)");
308  leg->SetBorderSize(0.0);
309  leg->SetFillColor(0.0);
310  leg->Draw();
311  TLegend* leg3 = new TLegend(0.6,0.75,0.9,0.85);
312  leg3->SetHeader("STAR 2006 Preliminary #pi^{0}");
313  leg3->SetBorderSize(0.0);
314  leg3->SetFillColor(0.0);
315  leg3->Draw();
316  c2->Print("all.gif");
317 #endif
318 
319 }
320 
321 void doCalall(ifstream& yield){
322  double epsilon_raw[7],weight[7],eps_lx_raw[7],eps_ly_raw[7],eps_lx_back[7],eps_ly_back[7];
323  string s1,s2,s3,s4,s5;
324  yield>>s1>>s2>>s3>>s4>>s5;
325  double totsum=0,realsum=0,backsum=0;
326  double psquare=0.3;
327  //cout<<s1<<endl;
328  for(int i=0;i<7;i++){
329  int flag=0;
330  int spin=0;
331  double totpi=0,realpi=0,backpi=0;
332  double pt=0.0;
333  double Nuu_sig=0,Nuu_real=0,Nuu_back=0,Nud_sig=0,Nud_real=0,Nud_back=0,Ndu_sig=0,Ndu_real=0,Ndu_back=0,Ndd_sig=0,Ndd_real=0,Ndd_back=0;
334  double sumtot=0,sumreal=0,sumback=0;
335  double sumpt=0.0;
336  double xx1=0.0,xx2=0.0,yy1=0.0,yy2=0.0;
337  double des=0.0,deb=0.0,der=0.0,dew=0.0,eepsilon_raw=0.0,eepsilon_real=0.0,eeps_lx_raw=0.0,eeps_ly_raw=0.0,eeps_lx_back=0.0,eeps_ly_back=0.0;
338  while(1)
339  {
340  yield>>spin>>pt>>totpi>>realpi>>backpi;
341  //cout<<"spin="<<spin<<" pt="<<pt<<" totpi="<<totpi<<" realpi="<<realpi<<" backpi="<<backpi<<endl;
342  sumtot+=totpi;
343  sumreal+=realpi;
344  sumback+=backpi;
345  sumpt+=pt*totpi;
346  totsum+=totpi;
347  realsum+=realpi;
348  backsum+=backpi;
349  //cout<<"sumpt="<<sumpt<<endl;
350  if(spin==0)
351  {
352  Nuu_sig+=totpi;
353  Nuu_real+=realpi;
354  Nuu_back+=backpi;
355  }
356  if(spin==1)
357  {
358  Nud_sig+=totpi;
359  Nud_real+=realpi;
360  Nud_back+=backpi;
361  }
362  if(spin==2)
363  {
364  Ndu_sig+=totpi;
365  Ndu_real+=realpi;
366  Ndu_back+=backpi;
367  }
368  if(spin==3)
369  {
370  Ndd_sig+=totpi;
371  Ndd_real+=realpi;
372  Ndd_back+=backpi;
373  }
374  flag++;
375  if(flag==4) break;
376  //if(flag==8) break;
377  }//end of while
378  xx1=Ndd_real+Ndu_real;
379  xx2=Nud_real+Nuu_real;
380  yy1=Ndd_real+Nud_real;
381  yy2=Ndu_real+Nuu_real;
382  epsilon_X[i]=-(xx1-xx2)/(xx1+xx2);
383  epsilon_Y[i]=-(yy1-yy2)/(yy1+yy2);
384  ptmean[i]=sumpt/sumtot;
385 
386  eepsilon_raw=(Nuu_sig-Nud_sig-Ndu_sig+Ndd_sig)/sumtot;
387  eepsilon_real=(Nuu_real-Nud_real-Ndu_real+Ndd_real)/sumreal;
388  eeps_lx_raw=(Ndd_sig-Nud_sig-Nuu_sig+Ndu_sig)/sumback;
389  eeps_ly_raw=(Ndd_sig+Nud_sig-Nuu_sig-Ndu_sig)/sumback;
390  eeps_lx_back=(Ndd_back-Nud_back-Nuu_back+Ndu_back)/sumback;
391  eeps_ly_back=(Ndd_back+Nud_back-Nuu_back-Ndu_back)/sumback;
392  epsilon_raw[i]=1.0/psquare*eepsilon_raw;
393  //all calculation from pure pi0 yield after background subtraction.
394  epsilon_real[i]=1.0/psquare*eepsilon_real;
395  eps_lx_raw[i]=1.0/psquare*eeps_lx_raw;
396  eps_ly_raw[i]=1.0/psquare*eeps_ly_raw;
397  eps_lx_back[i]=1.0/psquare*eeps_lx_back;
398  eps_ly_back[i]=1.0/psquare*eeps_ly_back;
399  weight[i]=sumback/sumtot;
400 
401  epsilon_back[i]=1.0/psquare*(Nuu_back-Nud_back-Ndu_back+Ndd_back)/sumback;
402  //all calculation
403  epsilon_ll[i]=(epsilon_raw[i]-weight[i]*epsilon_back[i])/(1.0-weight[i]);
404  //stat. error calculation
405  des=sqrt(sumtot)/sumtot;
406  der=sqrt(sumreal)/sumreal;
407  deb=sqrt(sumback)/sumback;
408  dew=sqrt(sumback)/sumtot;
409  ds[i]=1.0/psquare*des;
410  dr[i]=1.0/psquare*der;
411  db[i]=1.0/psquare*deb;
412  dw[i]=1.0/psquare*dew;
413 
414 
415  ep_error[i]=sqrt(pow(ds[i]/(1-weight[i]),2)+pow((weight[i]*db[i])/(1-weight[i]),2)+pow(dw[i]*(epsilon_raw[i]-epsilon_back[i])/pow(1-weight[i],2),2));
416  dx[i]=sqrt(pow(des/(1-weight[i]),2)+pow((weight[i]*deb)/(1-weight[i]),2)+pow(dew*(eeps_lx_raw-eeps_lx_back)/pow(1-weight[i],2),2));
417  dy[i]=sqrt(pow(des/(1-weight[i]),2)+pow((weight[i]*deb)/(1-weight[i]),2)+pow(dew*(eeps_ly_raw-eeps_ly_back)/pow(1-weight[i],2),2));
418  //the second way to calculate error of all
419  ep_error_rb[i]=sqrt(pow(dr[i],2)+pow(db[i],2));
420 
421 
422  }//end of for looping pt bins
423 
424 
425 }
Generated by   doxygen 1.8.5