StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
getSpectrumPP05.C
1 
2 #include "systematics/systematics_pp.h"
3 
4 double pqcdlowpt(double *x,double *par)
5 {
6  double ret=par[0]*pow(1.+x[0],par[1]);
7  return ret;
8 }
9 double pqcdhighpt(double *x,double *par)
10 {
11  double ret=par[0]*pow(1.+x[0],par[1]);
12  return ret;
13 }
14 double sumpqcd(double *x,double *par)
15 {
16  //transition (SW=1 is old)
17  //params: 2x low, 2x hi, 2x trans
18  double SW=1. - 1./(1.+exp(fabs(par[4])*(x[0]-par[5])));
19  double ret=SW*pqcdhighpt(x,&par[2]);
20  ret+=(1.-SW)*pqcdlowpt(x,par);
21  //ret=pqcdlowpt(x,par);
22  return ret;
23 }
24 void printPoints(TGraphErrors *g)
25 {
26  for(Int_t i=0;i<g->GetN();i++){
27  double x,y,ey;
28  g->GetPoint(i,x,y);
29  ey=g->GetErrorY(i);
30  cout<<x<<" "<<y<<" "<<ey<<endl;
31  }
32 }
33 void divideGraphWithFunction(TGraphErrors *graph,TF1 *func)
34 {
35  for(int i=0;i<graph->GetN();i++){
36  double x=0.;
37  double y=0.;
38  graph->GetPoint(i,x,y);
39  double ey=graph->GetErrorY(i);
40  graph->SetPoint(i,x,y/func->Eval(x));
41  graph->SetPointError(i,0.,ey/func->Eval(x));
42  }
43 }
44 void divideGraphWithGraph(TGraphErrors *graph,TGraph *denom)
45 {
46  for(int i=0;i<graph->GetN();i++){
47  double x=0.;
48  double y=0.;
49  graph->GetPoint(i,x,y);
50  double ey=graph->GetErrorY(i);
51  graph->SetPoint(i,x,y/denom->Eval(x));
52  graph->SetPointError(i,0.,ey/denom->Eval(x));
53  }
54 }
55 void getRatio(TH1F* h_ratio,TH1F* h_eff,TH1F* h_effSingle,TF1 *f_piondecay)
56 {
57  //copy, don't want to mess with the pointers:
58  TH1F *h_eff_copy=new TH1F(*h_eff);
59  for(int i=1;i<=h_eff->GetNbinsX();i++){
60  h_eff->SetBinError(i,0.);
61  }
62  TH1F *h_effSingle_copy=new TH1F(*h_effSingle);
63  for(int i=1;i<=h_effSingle->GetNbinsX();i++){
64  h_effSingle_copy->SetBinError(i,0.);
65  }
66  //correct for single photon fraction:
67  h_ratio->Add(f_piondecay,-1.);
68  for(int i=1;i<=h_ratio->GetNbinsX();i++){
69  //cout<<i<<endl;
70  //if(h_ratio->GetBinContent(i)<0.) cout<<"negative"<<endl;
71  }
72  h_eff_copy->Divide(h_effSingle_copy);
73  h_ratio->Multiply(h_eff_copy);
74  h_ratio->Add(f_piondecay);
75 }
76 void removeThesePoints(TGraphErrors *g,int trig)
77 {
78  if(trig==1){
79  g->RemovePoint(0);
80  g->RemovePoint(0);
81  //additional for pp:
82  g->RemovePoint(g->GetN()-1);
83  }
84  if(trig==2){
85  g->RemovePoint(0);
86  g->RemovePoint(0);
87  g->RemovePoint(0);
88  g->RemovePoint(0);
89  g->RemovePoint(g->GetN()-1);
90  g->RemovePoint(g->GetN()-1);
91  }
92  if(trig==3){
93  g->RemovePoint(0);
94  g->RemovePoint(0);
95  g->RemovePoint(0);
96  g->RemovePoint(0);
97  g->RemovePoint(0);
98  }
99 }
100 void getSpectrumPP05(const char *pid)
101 {
102  gStyle->SetOptStat(0);
103  gStyle->SetOptFit(0);
104 
105  gStyle->SetErrorX(0);
106 
107 
108 
109  bool subtractNEUTRONS=kTRUE;
110 
111  //neutron data:
112  //hadron data for Levy function, nucl-ex/0601033:
113  //--> d^2N/(2*pi*pT*N*dpT*dy) = B/((1+((mT - m0)/nT))^n)
114  // {p-dAu; pbar-dAu; p-pp; pbar-pp} and m0 = m_neutron = 1.0 GeV.
115  double B[]={0.3,0.23,0.072,0.061};//->[0]
116  //double eB[]={0.01,0.01,0.005,0.005};
117  double T[]={0.205,0.215,0.179,0.173};//->[1]
118  //double eT[]={0.004,0.005,0.006,0.006};
119  double n[]={11.00,12.55,10.87,10.49};//->[2]
120  //double en[]={0.29,0.41,0.43,0.40};
121  double BR=35.8/63.9;//->p + pi- (63.9%) : ->n + pi0 (35.8%)
122  double c_feeddown=0.8+0.2*BR;
123  double CPVloss=0.98;//for minbias
124  TF1 *f_nbar=new TF1("f_nbar","[3]*[4]*[0]/pow((1.+(sqrt(x*x+1.) - 1.)/([1]*[2])),[2])",0.,15.);
125  f_nbar->SetParameters(B[3],T[3],n[3],c_feeddown,CPVloss);
126 
127  //get direct gammas pQCD:
128  ifstream pQCDphotons("./datapoints/pQCD_Werner/rhic_cteq6_gamma_inv_sc1.dat");
129  float ppx[100];
130  float ppy[100];
131  Int_t iii=0;
132  cout<<"pqcd photons:"<<endl;
133  while(iii<28){
134  if(!pQCDphotons.good()) break;
135  float dummy=0.;
136  pQCDphotons>>ppx[iii]>>dummy>>dummy>>ppy[iii];
137  ppy[iii]*=1.e-09;//convert to mb
138  iii++;
139  }
140  TGraph *g_dirgamma=new TGraph(iii,ppx,ppy);
141  //get direct gammas pQCD: scale 0.5*pT
142  ifstream pQCDphotons05("./datapoints/pQCD_Werner/rhic_cteq6_gamma_inv_sc05.dat");
143  float ppx05[100];
144  float ppy05[100];
145  Int_t iii05=0;
146  while(iii05<28){
147  if(!pQCDphotons05.good()) break;
148  float dummy=0.;
149  pQCDphotons05>>ppx05[iii05]>>dummy>>dummy>>ppy05[iii05];
150  ppy05[iii05]*=1.e-09;//convert to mb
151  iii05++;
152  }
153  TGraph *g_dirgamma05=new TGraph(iii05,ppx05,ppy05);
154 
155  //get direct gammas pQCD: scale 2*pT
156  ifstream pQCDphotons2("./datapoints/pQCD_Werner/rhic_cteq6_gamma_inv_sc2.dat");
157  float ppx2[100];
158  float ppy2[100];
159  Int_t iii2=0;
160  while(iii2<28){
161  if(!pQCDphotons2.good()) break;
162  float dummy=0.;
163  pQCDphotons2>>ppx2[iii2]>>dummy>>dummy>>ppy2[iii2];
164  ppy2[iii2]*=1.e-09;//convert to mb
165  iii2++;
166  }
167  TGraph *g_dirgamma2=new TGraph(iii2,ppx2,ppy2);
168  //get phenix pions in pp
169  ifstream phenix("./datapoints/phenix_xsec_pp.dat");
170  float phex[100];
171  float phey[100];
172  Int_t iphe=0;
173  while(iphe<16){
174  if(!phenix.good()) break;
175  phenix>>phex[iphe]>>phey[iphe];
176  //is in mb already!!
177  iphe++;
178  }
179  TGraphErrors *phenix_pp=new TGraphErrors(iphe,phex,phey);
180  phenix_pp->SetMarkerStyle(24);
181  phenix_pp->SetLineWidth(6);
182  phenix_pp->SetName("phenix_pp");
183 
184  //frank's spin2006 pions:
185  ifstream frank("./datapoints/frank_pp05_new.dat");
186  float frx[100];
187  float fry[100];
188  float frex[100];
189  float frey[100];
190  Int_t ifr=0;
191  while(ifr<12){
192  if(!frank.good()) break;
193  frank>>frx[ifr]>>fry[ifr]>>frey[ifr];
194  frex[ifr]=0.;
195  ifr++;
196  }
197  TGraphErrors *frank_pp=new TGraphErrors(ifr,frx,fry,frex,frey);
198  frank_pp->SetMarkerStyle(25);
199  frank_pp->SetMarkerColor(4);
200  frank_pp->SetName("frank_pp");
201 
202  //sasha's pions:
203  float x_pp2005_sasha[]={1.20358,1.70592,2.21238,2.71916,3.4032,4.4224,5.43571,6.44468,7.45056,8.83794,10.8659,12.8831,15.2692 };
204  float xe_pp2005_sasha[]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
205  float y_pp2005_sasha[]={0.212971,0.0337464,0.00811345,0.00248195,0.000467495,6.37683e-05,1.4487e-05,3.67539e-06,1.26723e-06,3.3676e-07,8.01941e-08,1.93813e-08,5.56059e-09};
206  float ye_pp2005_sasha[]={0.00463771,0.000949529,0.000323778,0.000139023,3.68181e-05,7.61527e-07,2.28359e-07,9.19064e-08,4.93116e-08,8.16844e-09,3.66273e-09,1.75057e-09,7.36722e-10};
207 
208  TGraphErrors *sasha_pp05=new TGraphErrors(13,x_pp2005_sasha,y_pp2005_sasha,xe_pp2005_sasha,ye_pp2005_sasha);
209  sasha_pp05->SetMarkerStyle(8);
210  sasha_pp05->SetMarkerColor(TColor::GetColor(8,80,8));
211  sasha_pp05->SetName("sasha_pp05");
212 
213  //get pions KKP scale pT
214  ifstream pQCDpions("./datapoints/pQCD_Werner/klaus_pi0inv_200_kkp_1.dat");
215  float pionx[100];
216  float piony[100];
217  Int_t ipion=0;
218  while(ipion<28){
219  if(!pQCDpions.good()) break;
220  pQCDpions>>pionx[ipion]>>piony[ipion];
221  piony[ipion]*=1.e-09;//convert to mb
222  ipion++;
223  }
224  TGraphErrors *kkp=new TGraphErrors(ipion,pionx,piony);
225  kkp->SetLineColor(54);
226  kkp->SetName("kkp");
227 
228  //get pions KKP scale 0.5*pT
229  ifstream pQCDpions05("./datapoints/pQCD_Werner/klaus_pi0inv_200_kkp_05.dat");
230  float pionx05[100];
231  float piony05[100];
232  Int_t ipion05=0;
233  while(ipion05<28){
234  if(!pQCDpions05.good()) break;
235  pQCDpions05>>pionx05[ipion05]>>piony05[ipion05];
236  piony05[ipion05]*=1.e-09;//convert to mb
237  ipion05++;
238  }
239  TGraphErrors *kkp05=new TGraphErrors(ipion05,pionx05,piony05);
240  kkp05->SetLineStyle(2);
241  kkp05->SetLineColor(54);
242  kkp05->SetName("kkp05");
243 
244  //get pions KKP scale 2*pT
245  ifstream pQCDpions2("./datapoints/pQCD_Werner/klaus_pi0inv_200_kkp_2.dat");
246  float pionx2[100];
247  float piony2[100];
248  Int_t ipion2=0;
249  while(ipion2<28){
250  if(!pQCDpions2.good()) break;
251  pQCDpions2>>pionx2[ipion2]>>piony2[ipion2];
252  piony2[ipion2]*=1.e-09;//convert to mb
253  ipion2++;
254  }
255  TGraphErrors *kkp2=new TGraphErrors(ipion2,pionx2,piony2);
256  kkp2->SetLineStyle(2);
257  kkp2->SetLineColor(54);
258  kkp2->SetName("kkp2");
259 
260  TFile *f_decaybg=new TFile("~/MyDecay/gammaDecayPPSum.root","OPEN");
261  TH1F *h_decaybg=(TH1F*)f_decaybg->Get("gamma");
262  TH1F *h_decaypion=(TH1F*)f_decaybg->Get("gamma_pion");
263  TF1 *fit_decay=new TF1("fit_decay","[0]/pow(x,[1])+[2]",.3,15.);
264  TF1 *fit_piondecay=new TF1(*fit_decay);
265  fit_decay->SetParameters(1.,1.,.5);
266  fit_piondecay->SetParameters(1.,1.,.5);
267  h_decaybg->Fit(fit_decay,"R0Q");
268  h_decaypion->Fit(fit_piondecay,"R0Q");
269 
270  //take ratio gamma_direct/pion and divide by gamma/pi and then +1:
271  for(Int_t i=0;i<iii;i++){
272  ppy[i]=ppy[i]/piony[i];
273  ppy[i]=ppy[i]/fit_decay->Eval(ppx[i]);
274  ppy[i]+=1.;
275  ppy05[i]=ppy05[i]/piony05[i];
276  ppy05[i]=ppy05[i]/fit_decay->Eval(ppx05[i]);
277  ppy05[i]+=1.;
278  ppy2[i]=ppy2[i]/piony2[i];
279  ppy2[i]=ppy2[i]/fit_decay->Eval(ppx2[i]);
280  ppy2[i]+=1.;
281  }
282  TGraphErrors *g_photonpqcd=new TGraphErrors(iii,ppx,ppy);
283  g_photonpqcd->SetName("g_photonpqcd");
284  TGraphErrors *g_photonpqcd05=new TGraphErrors(iii05,ppx05,ppy05);
285  g_photonpqcd05->SetName("g_photonpqcd05");
286  TGraphErrors *g_photonpqcd2=new TGraphErrors(iii2,ppx2,ppy2);
287  g_photonpqcd2->SetName("g_photonpqcd2");
288 
289  //set outputfiles
290  TString dir("/star/u/russcher/gamma/analysis/output/pp05/");
291  dir.Append(pid);
292  dir.Append("/");
293  dir.Append(pid);
294  TString psout=dir;
295  TString psout_eta=dir;
296  TString psout2=dir;
297  TString psout3=dir;
298  TString eFile=dir;
299  TString eFileGamma=dir;
300  TString pi0File=dir;
301  TString nbarFile=dir;
302  TString nbarOut=dir;
303 
304  eFile.Append("pion_eff.root");
305  eFileGamma.Append("gamma_eff.root");
306  pi0File.Append("pi0_pp05.root");
307  nbarFile.Append("antineutron_eff.root");
308  nbarOut.Append("nbar_contam.ps");
309  psout2.Append("correctedspec.ps");
310  psout_eta.Append("/dev/null");
311  psout.Append("invmassplots.ps");
312 
313  //load *.so
314  gSystem->Load("$HOME/MyEvent/MyEvent");
315  gSystem->Load("$HOME/gamma/analysis/lib/AnaCuts");
316  gSystem->Load("$HOME/gamma/analysis/lib/EventMixer");
317  gSystem->Load("$HOME/gamma/analysis/lib/Pi0Analysis");
318 
319 
320  AnaCuts *cuts=new AnaCuts("pp05");
321 
322  //get inv mass hists
323  TFile f(pi0File.Data(),"OPEN");
324  TH2F *h_mb=new TH2F(*h_minvMB);
325  TH2F *h_ht1=new TH2F(*h_minvHT1);
326  TH2F *h_ht2=new TH2F(*h_minvHT2);
327  TH1F *h_ev=new TH1F(*h_events);
328 
329  //get anti-neutron
330  TFile *file_nbar=new TFile(nbarFile,"OPEN");
331  TH1F *h_nbarEffMB=new TH1F(*h_effMB);
332  h_nbarEffMB->Sumw2();
333  TH1F *h_nbarEffHT1=new TH1F(*h_effHT1);
334  h_nbarEffHT1->Sumw2();
335  TH1F *h_nbarEffHT2=new TH1F(*h_effHT2);
336  h_nbarEffHT2->Sumw2();
337 
338  //get prescales
339  int trigger=0;
340  Int_t numberOfMB=0;
341  Int_t numberOfHT1=0;
342  Int_t numberOfHT2=0;
343  for(Int_t i=1;i<=h_ev->GetNbinsX();i++)
344  {
345  trigger=(Int_t)h_ev->GetBinCenter(i);
346  if(trigger&1) numberOfMB+=(Int_t)h_ev->GetBinContent(i);
347  if(trigger&2) numberOfHT1+=(Int_t)h_ev->GetBinContent(i);
348  if(trigger&4) numberOfHT2+=(Int_t)h_ev->GetBinContent(i);
349  }
350 
351  cout<<"number of mb: "<<numberOfMB<<endl;
352  //numberOfMB/=0.63;
353  //cout<<"nmb after 62% vertex eff.: "<<numberOfMB<<endl;
354 
355  float psMB=32742;
356  float psHT1=3.89;
357  float psHT2=1.;
358 
359  //mb events for hightower:
360  float nMBwithHT=328871;
361 
362  //get efficiencies+acceptance
363  TFile g(eFile.Data(),"OPEN");
364  TH1F *h_emb=new TH1F(*h_effMB);
365  TH1F *h_eht1=new TH1F(*h_effHT1);
366  TH1F *h_eht2=new TH1F(*h_effHT2);
367 
368  //bin corrections
369  TFile binf("~/BinWidth/bincorrectionsPP.root","OPEN");
370  TH1F *h_binmb=new TH1F(*h4mb);
371  TH1F *h_binht1=new TH1F(*h4ht1);
372  TH1F *h_binht2=new TH1F(*h4ht2);
373  h_binmb->Sumw2();
374  h_binht1->Sumw2();
375  h_binht2->Sumw2();
376  for(Int_t i=1;i<=h_binmb->GetNbinsX();i++) h_binmb->SetBinError(i,0);
377  for(Int_t i=1;i<=h_binht1->GetNbinsX();i++) h_binht1->SetBinError(i,0);
378  for(Int_t i=1;i<=h_binht2->GetNbinsX();i++) h_binht2->SetBinError(i,0);
379 
380  //corrections, all multiplicative, in case of
381  //pions:
382  TF1 *pion_cpv_corrMB=new TF1("pion_cpv_corrMB","1./(1.-0.01*(0.3+0.0*x))",0.,15.);
383  TF1 *pion_cpv_corrHT1=new TF1("pion_cpv_corrHT1","1./(1.-0.01*(-0.1+0.16*x))",0.,15.);
384  TF1 *pion_cpv_corrHT2=new TF1("pion_cpv_corrHT2","1./(1.-0.01*(-0.2+0.18*x))",0.,15.);
385  //photons:
386  TF1 *gamma_cpv_corrMB=new TF1("gamma_cpv_corrMB","1./(1.-0.01*(2.8+0.0*x))",0.,15.);
387  TF1 *gamma_cpv_corrHT1=new TF1("gamma_cpv_corrHT1","1./(1.-0.01*(0.2+1.1*x))",0.,15.);
388  TF1 *gamma_cpv_corrHT2=new TF1("gamma_cpv_corrHT2","1./(1.-0.01*(0.4+1.1*x))",0.,15.);
389  TF1 *gamma_cont_corrMB=new TF1("gamma_cont_corrMB","0.985",0.,15.);
390  TF1 *gamma_cont_corrHT1=new TF1("gamma_cont_corrHT1","0.98",0.,15.);
391  TF1 *gamma_cont_corrHT2=new TF1("gamma_cont_corrHT2","0.96",0.,15.);
392 
393  //missing material
394  //pions:
395  TF1 *pion_conv_corrMB=new TF1("pion_conv_corrMB","1.1",0.,15.);
396  TF1 *pion_conv_corrHT1=new TF1("pion_conv_corrHT1","1.1",0.,15.);
397  TF1 *pion_conv_corrHT2=new TF1("pion_conv_corrHT2","1.1",0.,15.);
398  //photons:
399  TF1 *gamma_conv_corrMB=new TF1("gamma_conv_corrMB","1.05",0.,15.);
400  TF1 *gamma_conv_corrHT1=new TF1("gamma_conv_corrHT1","1.05",0.,15.);
401  TF1 *gamma_conv_corrHT2=new TF1("gamma_conv_corrHT2","1.05",0.,15.);
402 
403 
404  //get yield
405  Pi0Analysis *pi0=new Pi0Analysis(psout.Data(),psout_eta.Data(),"pp05");
406  pi0->init("/dev/null");
407  TH1F *pionYieldMB=new TH1F(*pi0->getYield(h_mb,"mb"));
408  TH1F *pionYieldHT1=new TH1F(*pi0->getYield(h_ht1,"ht1"));
409  TH1F *pionYieldHT2=new TH1F(*pi0->getYield(h_ht2,"ht2"));
410  pi0->storeCanvases((dir+"canvases.root").Data());
411 
412 
413  cout<<"***************************************"<<endl;
414  cout<<"got yield, dividing by rapidity bite!!!"<<endl;
415  float dy_gamma=cuts->rapidityMaxCUT - cuts->rapidityMinCUT;
416  float dy_pion=cuts->rapPionMaxCUT - cuts->rapPionMinCUT;
417  cout<<"***************************************"<<endl;
418  cout<<endl;
419  cout<<" pion bite is "<<dy_pion<<endl;
420  cout<<" gamma bite is "<<dy_gamma<<endl;
421  cout<<endl;
422  cout<<"***************************************"<<endl;
423 
424  pionYieldMB->Scale(1./dy_pion);
425  pionYieldHT1->Scale(1./dy_pion);
426  pionYieldHT2->Scale(1./dy_pion);
427 
428  //set yield to zero
429  /*
430  for(Int_t i=0;i<pionYieldHT2->GetNbinsX();i++)
431  {
432  if(i<1)
433  {
434  pionYieldMB->SetBinContent(i+1,0);
435  pionYieldMB->SetBinError(i+1,0);
436  }
437  if(i<4)
438  {
439  pionYieldHT1->SetBinContent(i+1,0);
440  pionYieldHT1->SetBinError(i+1,0);
441  }
442  if(i>12)
443  {
444  pionYieldHT1->SetBinContent(i+1,0);
445  pionYieldHT1->SetBinError(i+1,0);
446  }
447  if(i<6)
448  {
449  pionYieldHT2->SetBinContent(i+1,0);
450  pionYieldHT2->SetBinError(i+1,0);
451  }
452  }
453  */
454 
455  //set colors:
456  pionYieldMB->SetMarkerStyle(8);
457  pionYieldMB->SetMarkerSize(1.0);
458  pionYieldHT1->SetMarkerStyle(8);
459  pionYieldHT1->SetMarkerSize(1.0);
460  pionYieldHT1->SetMarkerColor(4);
461  pionYieldHT2->SetMarkerStyle(8);
462  pionYieldHT2->SetMarkerSize(1.0);
463  pionYieldHT2->SetMarkerColor(2);
464 
465  TF1 *scale=new TF1("scale","x",0.,15.);
466 
467  pionYieldMB->SetNameTitle("pionYieldMB","corrected yield MB");
468  pionYieldMB->Divide(h_emb);
469  pionYieldMB->Scale(psMB/(psMB*numberOfMB*2.*TMath::Pi()));
470  pionYieldMB->Divide(scale);
471  pionYieldMB->Multiply(h_binmb);
472  pionYieldMB->Multiply(pion_cpv_corrMB);
473  pionYieldMB->Multiply(pion_conv_corrMB);
474 
475  pionYieldHT1->SetNameTitle("pionYieldHT1","corrected yield HT1");
476  pionYieldHT1->Divide(h_eht1);
477  pionYieldHT1->Scale(psHT1/(psMB*nMBwithHT*2.*TMath::Pi()));
478  pionYieldHT1->Divide(scale);
479  pionYieldHT1->Multiply(h_binht1);
480  pionYieldHT1->Multiply(pion_cpv_corrHT1);
481  pionYieldHT1->Multiply(pion_conv_corrHT1);
482 
483  pionYieldHT2->SetNameTitle("pionYieldHT2","corrected yield HT2");
484  pionYieldHT2->Divide(h_eht2);
485  pionYieldHT2->Scale(psHT2/(psMB*nMBwithHT*2.*TMath::Pi()));
486  pionYieldHT2->Divide(scale);
487  pionYieldHT2->Multiply(h_binht2);
488  pionYieldHT2->Multiply(pion_cpv_corrHT2);
489  pionYieldHT2->Multiply(pion_conv_corrHT2);
490 
491  //create pion yield for double ratio:
492  TH1F *pionYieldMBratio=new TH1F(*pionYieldMB);
493  TH1F *pionYieldHT1ratio=new TH1F(*pionYieldHT1);
494  TH1F *pionYieldHT2ratio=new TH1F(*pionYieldHT2);
495 
496 
497  TH1F *pionXsMB=new TH1F(*pionYieldMB);
498  TH1F *pionXsHT1=new TH1F(*pionYieldHT1);
499  TH1F *pionXsHT2=new TH1F(*pionYieldHT2);
500  pionXsMB->Scale((26.1/0.85));//times xs_bbc/bbc_eff
501  pionXsHT1->Scale((26.1/0.85));
502  pionXsHT2->Scale((26.1/0.85));
503 
504 
505  TH1F *pionXsMBnoErr=new TH1F(*pionXsMB);
506  TH1F *pionXsHT1noErr=new TH1F(*pionXsHT1);
507  TH1F *pionXsHT2noErr=new TH1F(*pionXsHT2);
508  for(int i=1;i<=pionXsMBnoErr->GetNbinsX();i++){
509  pionXsMBnoErr->SetBinError(i,0.);
510  }
511  for(int i=1;i<=pionXsHT1noErr->GetNbinsX();i++){
512  pionXsHT1noErr->SetBinError(i,0.);
513  }
514  for(int i=1;i<=pionXsHT2noErr->GetNbinsX();i++){
515  pionXsHT2noErr->SetBinError(i,0.);
516  }
517 
518  TGraphErrors *g_pionXsMB=new TGraphErrors(pionXsMB);
519  g_pionXsMB->SetName("g_pionXsMB");
520  removeThesePoints(g_pionXsMB,1);
521 
522  TGraphErrors *g_pionXsHT1=new TGraphErrors(pionXsHT1);
523  g_pionXsHT1->SetName("g_pionXsHT1");
524  removeThesePoints(g_pionXsHT1,2);
525 
526  TGraphErrors *g_pionXsHT2=new TGraphErrors(pionXsHT2);
527  g_pionXsHT2->SetName("g_pionXsHT2");
528  removeThesePoints(g_pionXsHT2,3);
529 
530  if(0){
531  cout<<endl<<"xsec: x y ex ey"<<endl;
532  cout<<"minbias"<<endl;
533  printPoints(g_pionXsMB);
534  cout<<endl<<"hightower-1"<<endl;
535  printPoints(g_pionXsHT1);
536  cout<<endl<<"hightower-2"<<endl;
537  printPoints(g_pionXsHT2);
538  cout<<endl;
539  }
540 
541 
542  TMultiGraph *m_pions_fit=new TMultiGraph();
543  m_pions_fit->SetName("m_pions_fit");
544  m_pions_fit->SetMinimum(5.0e-9);
545  m_pions_fit->SetMaximum(0.99);
546 
547 
548  m_pions_fit->Add(g_pionXsMB);
549  m_pions_fit->Add(g_pionXsHT1);
550  m_pions_fit->Add(g_pionXsHT2);
551 
552  TF1 *fitQCD=new TF1("fitQCD",sumpqcd,1.,15.,6);
553  fitQCD->SetParameters(600.,-8.2,4.,-8.5,2.,2.);
554  fitQCD->FixParameter(4,2.);
555  m_pions_fit->Fit(fitQCD,"R");
556 
557  bool inclPhenix=false;
558  bool inclFrank=false;
559  bool inclSasha=false;
560  bool inclPqcd=false;
561 
562  TCanvas *compare=new TCanvas("compare","compare;p_{T}:xsec (mb)",600,750);
563  compare->cd();
564 
565  TPad *padt=new TPad("padt","",0.0,0.3,1.,1.0);
566  padt->SetBottomMargin(0.001);
567  padt->SetLeftMargin(0.15);
568  TPad *padb=new TPad("padb","",0.0,0.0,1.,0.3);
569  padb->SetTopMargin(0.001);
570  padb->SetBottomMargin(0.25);
571  padb->SetLeftMargin(0.15);
572 
573  padt->Draw();
574  padt->cd();
575  gPad->SetLogy();
576 
577  TMultiGraph *m_pions=new TMultiGraph();
578  m_pions->SetName("m_pions");
579 
580  if(inclPqcd){
581  m_pions->Add(kkp,"c");
582  m_pions->Add(kkp05,"c");
583  m_pions->Add(kkp2,"c");
584  }
585  if(inclSasha){
586  m_pions->Add(sasha_pp05);
587  }
588  if(inclFrank){
589  m_pions->Add(frank_pp);
590  }
591  if(inclPhenix){
592  m_pions->Add(phenix_pp);
593  }
594 
595  //g_pionXsMB->Print();
596  //cout<<endl<<endl;
597  //g_pionXsHT1->Print();
598  //cout<<endl<<endl;
599  //g_pionXsHT2->Print();
600  //cout<<endl;
601 
602  m_pions->Add(g_pionXsMB);
603  m_pions->Add(g_pionXsHT1);
604  m_pions->Add(g_pionXsHT2);
605 
606  m_pions->SetMinimum(1.0e-9);
607  m_pions->SetMaximum(1.);
608 
609  m_pions->Draw("ap");
610 
611  //fitQCD->Draw("same");
612 
613  m_pions->GetXaxis()->SetLabelSize(0.);
614  m_pions->GetYaxis()->SetTitle("Ed^{3}#sigma/dp^{3} (mb GeV^{-2}c^{2})");
615 
616  TLegend *leg=new TLegend(.5,.5,.85,.85);
617 
618  if(inclPhenix) leg->AddEntry(phenix_pp,"PHENIX p+p","p");
619  if(inclFrank) leg->AddEntry(frank_pp,"STAR preliminary (upd.)","p");
620  if(inclSasha) leg->AddEntry(sasha_pp05,"O.Grebenyuk p+p","p");
621  leg->AddEntry(g_pionXsMB,"p+p minimum bias","p");
622  leg->AddEntry(g_pionXsHT1,"hightower 1","p");
623  leg->AddEntry(g_pionXsHT2,"hightower 2","p");
624  if(inclPqcd){
625  leg->AddEntry(kkp,"kkp + CTEQ6m, #mu=p_{T}","l");
626  leg->AddEntry(kkp2,"#mu=2p_{T},p_{T}/2","l");
627  leg->Draw("same");
628  }
629 
630  leg->SetFillColor(0);
631  leg->Draw();
632 
633  compare->cd();
634  padb->Draw();
635  padb->cd();
636 
637  TGraphErrors *sasha_pp05_overPqcd=new TGraphErrors(*sasha_pp05);
638  divideGraphWithGraph(sasha_pp05_overPqcd,kkp);
639  TGraphErrors *phenix_pp05_overPqcd=new TGraphErrors(*phenix_pp);
640  divideGraphWithGraph(phenix_pp05_overPqcd,kkp);
641  TGraphErrors *g_pionXsMB_overPqcd=new TGraphErrors(*g_pionXsMB);
642  divideGraphWithGraph(g_pionXsMB_overPqcd,kkp);
643  TGraphErrors *g_pionXsHT1_overPqcd=new TGraphErrors(*g_pionXsHT1);
644  divideGraphWithGraph(g_pionXsHT1_overPqcd,kkp);
645  TGraphErrors *g_pionXsHT2_overPqcd=new TGraphErrors(*g_pionXsHT2);
646  divideGraphWithGraph(g_pionXsHT2_overPqcd,kkp);
647  TGraphErrors *frank_pp05_overPqcd=new TGraphErrors(*frank_pp);
648  divideGraphWithGraph(frank_pp05_overPqcd,kkp);
649 
650  TGraphErrors *kkp05_ratio=new TGraphErrors(*kkp05);
651  divideGraphWithGraph(kkp05_ratio,kkp);
652  TGraphErrors *kkp_ratio=new TGraphErrors(*kkp);
653  divideGraphWithGraph(kkp_ratio,kkp);
654  TGraphErrors *kkp2_ratio=new TGraphErrors(*kkp2);
655  divideGraphWithGraph(kkp2_ratio,kkp);
656 
657  //systematic errors:
658  TGraphErrors *g_pionXsMB_sys=new TGraphErrors(*g_pionXsMB_overPqcd);
659  set_sys_pp_pion(g_pionXsMB_sys);
660  TGraphErrors *g_pionXsHT1_sys=new TGraphErrors(*g_pionXsHT1_overPqcd);
661  set_sys_pp_pion(g_pionXsHT1_sys);
662  TGraphErrors *g_pionXsHT2_sys=new TGraphErrors(*g_pionXsHT2_overPqcd);
663  set_sys_pp_pion(g_pionXsHT2_sys);
664 
665  TMultiGraph *m_pions_over_pqcd=new TMultiGraph();
666 
667  m_pions_over_pqcd->SetMinimum(0.0);
668  m_pions_over_pqcd->SetMaximum(2.5);
669 
670  if(inclPqcd){
671  m_pions_over_pqcd->Add(kkp05_ratio,"c");
672  m_pions_over_pqcd->Add(kkp_ratio,"c");
673  m_pions_over_pqcd->Add(kkp2_ratio,"c");
674  }
675  m_pions_over_pqcd->Add(g_pionXsMB_overPqcd);
676  m_pions_over_pqcd->Add(g_pionXsHT1_overPqcd);
677  m_pions_over_pqcd->Add(g_pionXsHT2_overPqcd);
678  //m_pions_over_pqcd->Add(g_pionXsMB_sys,"c");
679  //m_pions_over_pqcd->Add(g_pionXsHT1_sys,"c");
680  //m_pions_over_pqcd->Add(g_pionXsHT2_sys,"c");
681  if(inclPhenix) m_pions_over_pqcd->Add(phenix_pp05_overPqcd);
682  if(inclFrank) m_pions_over_pqcd->Add(frank_pp05_overPqcd);
683  if(inclSasha) m_pions_over_pqcd->Add(sasha_pp05_overPqcd);
684 
685  //m_pions_over_pqcd->Draw("ap");
686  //m_pions_over_pqcd->GetXaxis()->SetTitle("p_{T} (GeV/c)");
687 
688  compare->SaveAs((dir+"pionxsec_pp.eps").Data());
689  compare->SaveAs((dir+"pionxsec_pp.root").Data());
690 
691 
692  TMultiGraph *m_pions_over_fit=new TMultiGraph();
693  m_pions_over_fit->SetMinimum(0.01);
694  m_pions_over_fit->SetMaximum(1.99);
695 
696  TGraphErrors *g_pionXsMB_overFit=new TGraphErrors(*g_pionXsMB);
697  divideGraphWithFunction(g_pionXsMB_overFit,fitQCD);
698  TGraphErrors *g_pionXsHT1_overFit=new TGraphErrors(*g_pionXsHT1);
699  divideGraphWithFunction(g_pionXsHT1_overFit,fitQCD);
700  TGraphErrors *g_pionXsHT2_overFit=new TGraphErrors(*g_pionXsHT2);
701  divideGraphWithFunction(g_pionXsHT2_overFit,fitQCD);
702 
703  m_pions_over_fit->Add(g_pionXsMB_overFit);
704  m_pions_over_fit->Add(g_pionXsHT1_overFit);
705  m_pions_over_fit->Add(g_pionXsHT2_overFit);
706 
707  m_pions_over_fit->Draw("ap");
708 
709  compare->SaveAs((dir+"pionxsec_pp_overfit.eps").Data());
710  compare->SaveAs((dir+"pionxsec_pp_overfit.root").Data());
711 
712 
713 
714  TCanvas *compare2=new TCanvas("compare2","compare2;p_{T};yield divided by fit",600,300);
715  compare2->cd();
716 
717  //divide by fit:
718  TGraphErrors *frank_pp2=new TGraphErrors(*frank_pp);
719  divideGraphWithFunction(frank_pp2,fitQCD);
720  TGraphErrors *sasha_pp2=new TGraphErrors(*sasha_pp05);
721  divideGraphWithFunction(sasha_pp2,fitQCD);
722  TGraphErrors *phenix_pp2=new TGraphErrors(*phenix_pp);
723  divideGraphWithFunction(phenix_pp2,fitQCD);
724  TGraphErrors *g_pionXsMBcopy=new TGraphErrors(*g_pionXsMB);
725  divideGraphWithFunction(g_pionXsMBcopy,fitQCD);
726  TGraphErrors *g_pionXsHT1copy=new TGraphErrors(*g_pionXsHT1);
727  divideGraphWithFunction(g_pionXsHT1copy,fitQCD);
728  TGraphErrors *g_pionXsHT2copy=new TGraphErrors(*g_pionXsHT2);
729  divideGraphWithFunction(g_pionXsHT2copy,fitQCD);
730 
731  inclPhenix=true;
732  inclFrank=false;
733  inclSasha=false;
734  inclPqcd=false;
735 
736  TMultiGraph *m_pions2=new TMultiGraph();
737  m_pions2->SetName("m_pions2");
738  if(inclSasha) m_pions2->Add(sasha_pp2);
739  if(inclFrank) m_pions2->Add(frank_pp2);
740  if(inclPhenix) m_pions2->Add(phenix_pp2);
741 
742  m_pions2->Add(g_pionXsMBcopy);
743  m_pions2->Add(g_pionXsHT1copy);
744  m_pions2->Add(g_pionXsHT2copy);
745 
746  m_pions2->SetMinimum(0.000001);
747  m_pions2->SetMaximum(3.);
748  m_pions2->Draw("ap");
749 
750  TLegend *legg=new TLegend(.25,.55,.65,.85);
751  legg->AddEntry(g_pionXsMBcopy,"minimum bias","p");
752  legg->AddEntry(g_pionXsHT1copy,"hightower 1","p");
753  legg->AddEntry(g_pionXsHT2copy,"hightower 2","p");
754  if(inclFrank) legg->AddEntry(frank_pp2,"Frank's p+p (upd.)","p");
755  if(inclSasha) legg->AddEntry(sasha_pp2,"Sasha's p+p","p");
756  if(inclPhenix) legg->AddEntry(phenix_pp,"PHENIX p+p","p");
757  legg->Draw("same");
758  legg->SetFillColor(0);
759 
760  compare2->cd(0);
761  compare2->SaveAs((dir+"pionxsec_pp_ratio.eps").Data());
762  compare2->SaveAs((dir+"pionxsec_pp_ratio.root").Data());
763 
764 
765  //********************************************
766  // Get double ratio:
767  //********************************************
768 
769  //pion decay photon eff:
770  TFile gg(eFile.Data(),"OPEN");
771  TH1F *effGammaMB=new TH1F(*h_effDaughtersMB);
772  TH1F *effGammaHT1=new TH1F(*h_effDaughtersHT1);
773  TH1F *effGammaHT2=new TH1F(*h_effDaughtersHT2);
774 
775  //single photon eff:
776  TFile gg_single(eFileGamma.Data(),"OPEN");
777  TH1F *effGammaSingleMB=new TH1F(*h_effMB);
778  TH1F *effGammaSingleHT1=new TH1F(*h_effHT1);
779  TH1F *effGammaSingleHT2=new TH1F(*h_effHT2);
780 
781  //raw neutral clusters:
782  TFile ff(pi0File.Data(),"OPEN");
783  TH1F *gammaYieldMB=new TH1F(*h_gammaMB);
784  TH1F *gammaYieldHT1=new TH1F(*h_gammaHT1);
785  TH1F *gammaYieldHT2=new TH1F(*h_gammaHT2);
786 
787  //divide rap. bite:
788  gammaYieldMB->Scale(1./dy_gamma);
789  gammaYieldHT1->Scale(1./dy_gamma);
790  gammaYieldHT2->Scale(1./dy_gamma);
791 
792 
793  for(Int_t i=1;i<=gammaYieldMB->GetNbinsX();i++){
794  gammaYieldMB->SetBinContent(i,gammaYieldMB->GetBinContent(i)/gammaYieldMB->GetXaxis()->GetBinWidth(i));
795  gammaYieldMB->SetBinError(i,gammaYieldMB->GetBinError(i)/gammaYieldMB->GetXaxis()->GetBinWidth(i));
796  }
797  gammaYieldMB->Scale(psMB/(psMB*numberOfMB*2.*TMath::Pi()));
798  gammaYieldMB->Divide(scale);// ../pT
799  gammaYieldMB->Multiply(h_binmb);
800  gammaYieldMB->Divide(effGammaMB);
801  gammaYieldMB->Multiply(gamma_cpv_corrMB);
802  gammaYieldMB->Multiply(gamma_cont_corrMB);
803  gammaYieldMB->Multiply(gamma_conv_corrMB);
804 
805  gammaYieldMB->SetMarkerStyle(8);
806  gammaYieldMB->SetMarkerSize(1.);
807  TGraphErrors *g_inclPhotonsMB=new TGraphErrors(gammaYieldMB);
808 
809 
810  for(Int_t i=1;i<=gammaYieldHT1->GetNbinsX();i++){
811  gammaYieldHT1->SetBinContent(i,gammaYieldHT1->GetBinContent(i)/gammaYieldHT1->GetXaxis()->GetBinWidth(i));
812  gammaYieldHT1->SetBinError(i,gammaYieldHT1->GetBinError(i)/gammaYieldHT1->GetXaxis()->GetBinWidth(i));
813  }
814  gammaYieldHT1->Scale(psHT1/(psMB*nMBwithHT*2.*TMath::Pi()));
815  gammaYieldHT1->Divide(scale);// ../pT
816  gammaYieldHT1->Multiply(h_binht1);
817  gammaYieldHT1->Divide(effGammaHT1);
818  gammaYieldHT1->Multiply(gamma_cpv_corrHT1);
819  gammaYieldHT1->Multiply(gamma_cont_corrHT1);
820  gammaYieldHT1->Multiply(gamma_conv_corrHT1);
821 
822  gammaYieldHT1->SetMarkerStyle(8);
823  gammaYieldHT1->SetMarkerSize(1.);
824  gammaYieldHT1->SetMarkerColor(4);
825  TGraphErrors *g_inclPhotonsHT1=new TGraphErrors(gammaYieldHT1);
826 
827  for(Int_t i=1;i<=gammaYieldHT2->GetNbinsX();i++){
828  gammaYieldHT2->SetBinContent(i,gammaYieldHT2->GetBinContent(i)/gammaYieldHT2->GetXaxis()->GetBinWidth(i));
829  gammaYieldHT2->SetBinError(i,gammaYieldHT2->GetBinError(i)/gammaYieldHT2->GetXaxis()->GetBinWidth(i));
830  }
831  gammaYieldHT2->Scale(psHT2/(psMB*nMBwithHT*2.*TMath::Pi()));
832  gammaYieldHT2->Divide(scale);// ../pT
833  gammaYieldHT2->Multiply(h_binht2);
834  gammaYieldHT2->Divide(effGammaHT2);
835  gammaYieldHT2->Multiply(gamma_cpv_corrHT2);
836  gammaYieldHT2->Multiply(gamma_cont_corrHT2);
837  gammaYieldHT2->Multiply(gamma_conv_corrHT2);
838 
839  gammaYieldHT2->SetMarkerStyle(8);
840  gammaYieldHT2->SetMarkerSize(1.);
841  gammaYieldHT2->SetMarkerColor(2);
842  TGraphErrors *g_inclPhotonsHT2=new TGraphErrors(gammaYieldHT2);
843 
844  removeThesePoints(g_inclPhotonsMB,1);
845  removeThesePoints(g_inclPhotonsHT1,2);
846  removeThesePoints(g_inclPhotonsHT2,2);
847 
848  TMultiGraph *m_incl=new TMultiGraph();
849  m_incl->SetName("m_incl");
850  m_incl->SetTitle("inclusive photon invariant yield 0<y<1;p_{T} (GeV/c);#frac{1}{2#piNp_{T}} #frac{d^{2}N}{dydp_{T}}");
851  m_incl->Add(g_inclPhotonsMB);
852  m_incl->Add(g_inclPhotonsHT1);
853  m_incl->Add(g_inclPhotonsHT2);
854 
855  m_incl->Fit(fitQCD,"R0");
856 
857  m_incl->SetMinimum(1.e-11);
858  m_incl->SetMaximum(10.);
859  TCanvas *c_incl=new TCanvas("c_incl","c_incl",600,400);
860  gPad->SetLogy();
861  m_incl->Draw("ap");
862  c_incl->SaveAs((dir+"inclPhotonYield.eps").Data());
863  c_incl->SaveAs((dir+"inclPhotonYield.root").Data());
864 
865 
866 
867  //get ratio:
868  TH1F *gammaYieldMBratio=new TH1F(*gammaYieldMB);
869  gammaYieldMBratio->SetName("gammaYieldMBratio");
870  TH1F *gammaYieldHT1ratio=new TH1F(*gammaYieldHT1);
871  gammaYieldHT1ratio->SetName("gammaYieldHT1ratio");
872  TH1F *gammaYieldHT2ratio=new TH1F(*gammaYieldHT2);
873  gammaYieldHT2ratio->SetName("gammaYieldHT2ratio");
874 
875  gammaYieldMBratio->Divide(pionYieldMBratio);
876  gammaYieldHT1ratio->Divide(pionYieldHT1ratio);
877  gammaYieldHT2ratio->Divide(pionYieldHT2ratio);
878 
879 
880  //correct gamma over pion ratio, using two efficiencies:
881  getRatio(gammaYieldMBratio,effGammaMB,effGammaSingleMB,fit_piondecay);
882  getRatio(gammaYieldHT1ratio,effGammaHT1,effGammaSingleHT1,fit_piondecay);
883  getRatio(gammaYieldHT2ratio,effGammaHT2,effGammaSingleHT2,fit_piondecay);
884 
885  TH1F *gammaYieldMBratio_incl=new TH1F(*gammaYieldMBratio);
886  TH1F *gammaYieldHT1ratio_incl=new TH1F(*gammaYieldHT1ratio);
887  TH1F *gammaYieldHT2ratio_incl=new TH1F(*gammaYieldHT2ratio);
888 
889  TH1F *gammaYieldMBratioNoErr=new TH1F(*gammaYieldMBratio);
890  TH1F *gammaYieldHT1ratioNoErr=new TH1F(*gammaYieldHT1ratio);
891  TH1F *gammaYieldHT2ratioNoErr=new TH1F(*gammaYieldHT2ratio);
892  for(int i=1;i<=gammaYieldMBratioNoErr->GetNbinsX();i++){
893  gammaYieldMBratioNoErr->SetBinError(i,0.);
894  }
895  for(int i=1;i<=gammaYieldHT1ratioNoErr->GetNbinsX();i++){
896  gammaYieldHT1ratioNoErr->SetBinError(i,0.);
897  }
898  for(int i=1;i<=gammaYieldHT2ratioNoErr->GetNbinsX();i++){
899  gammaYieldHT2ratioNoErr->SetBinError(i,0.);
900  }
901 
902  TGraphErrors *g_ratioMB=new TGraphErrors(gammaYieldMBratio);
903  g_ratioMB->SetName("g_ratioMB");
904  TGraphErrors *g_ratioHT1=new TGraphErrors(gammaYieldHT1ratio);
905  g_ratioHT1->SetName("g_ratioHT1");
906  TGraphErrors *g_ratioHT2=new TGraphErrors(gammaYieldHT2ratio);
907  g_ratioHT2->SetName("g_ratioHT2");
908 
909 
910  removeThesePoints(g_ratioMB,1);
911  removeThesePoints(g_ratioHT1,2);
912  removeThesePoints(g_ratioHT2,3);
913 
914 
915  TCanvas *c_ratio=new TCanvas("c_ratio","c_ratio",400,200);
916 
917  TMultiGraph *m_ratio=new TMultiGraph("m_ratio","p+p 2005;p_{T};#gamma/#pi^{0}");
918  m_ratio->Add(g_ratioMB);
919  m_ratio->Add(g_ratioHT1);
920  m_ratio->Add(g_ratioHT2);
921 
922  m_ratio->Draw("ap");
923  m_ratio->SetMinimum(.001);
924  m_ratio->SetMaximum(1.5);
925 
926  TLegend *leg3=new TLegend(.35,.65,.65,.85);
927  leg3->AddEntry(g_ratioMB,"minimum bias","p");
928  leg3->AddEntry(g_ratioHT1,"hightower 1","p");
929  leg3->AddEntry(g_ratioHT2,"hightower 2","p");
930  leg3->AddEntry(fit_decay,"decay background (total)","l");
931  leg3->AddEntry(fit_piondecay,"decay background (#pi^{0})","l");
932  leg3->SetFillColor(0);
933  leg3->Draw("same");
934 
935  fit_decay->SetLineColor(13);
936  fit_decay->SetLineWidth(1);
937  fit_decay->SetLineColor(1);
938  fit_decay->Draw("same");
939  fit_piondecay->SetLineColor(13);
940  fit_piondecay->SetLineWidth(1);
941  fit_piondecay->SetLineStyle(2);
942  fit_piondecay->SetLineColor(1);
943  fit_piondecay->Draw("same");
944 
945  c_ratio->SaveAs((dir+"gammaOverPion.eps").Data());
946  c_ratio->SaveAs((dir+"gammaOverPion.root").Data());
947 
948 
949  //create fully corrected incl. photons:
950  gammaYieldMBratio_incl->Multiply(pionYieldMBratio);
951  gammaYieldHT1ratio_incl->Multiply(pionYieldHT1ratio);
952  gammaYieldHT2ratio_incl->Multiply(pionYieldHT2ratio);
953  //gammaYieldMBratio_incl->Scale(7.5);
954  //gammaYieldHT1ratio_incl->Scale(7.5);
955  //gammaYieldHT2ratio_incl->Scale(7.5);
956 
957  TGraphErrors *g_incl_corrMB=new TGraphErrors(gammaYieldMBratio_incl);
958  TGraphErrors *g_incl_corrHT1=new TGraphErrors(gammaYieldHT1ratio_incl);
959  TGraphErrors *g_incl_corrHT2=new TGraphErrors(gammaYieldHT2ratio_incl);
960 
961  TCanvas *c_incl_corr=new TCanvas("c_incl_corr","c_incl_corr",400,300);
962  gPad->SetLogy();
963  TMultiGraph *m_incl_corr=new TMultiGraph();
964  m_incl_corr->Add(g_incl_corrMB);
965  m_incl_corr->Add(g_incl_corrHT1);
966  m_incl_corr->Add(g_incl_corrHT2);
967 
968  m_incl_corr->SetMinimum(1.e-11);
969  m_incl_corr->SetMaximum(1.);
970 
971  m_incl_corr->Draw("apX");
972  c_incl_corr->SaveAs((dir+"inclPhotonYieldCorr.eps").Data());
973  c_incl_corr->SaveAs((dir+"inclPhotonYieldCorr.root").Data());
974 
975 
976 
977  TCanvas *c_doubleratio=new TCanvas("c_doubleratio","c_doubleratio",400,300);
978  gStyle->SetOptStat(0);
979  c_doubleratio->cd(1);
980 
981  TH1F *gammaYieldMBdoubleratio=new TH1F(*gammaYieldMBratio);
982  TH1F *gammaYieldHT1doubleratio=new TH1F(*gammaYieldHT1ratio);
983  TH1F *gammaYieldHT2doubleratio=new TH1F(*gammaYieldHT2ratio);
984 
985  gammaYieldMBdoubleratio->Divide(fit_decay);
986  gammaYieldHT1doubleratio->Divide(fit_decay);
987  gammaYieldHT2doubleratio->Divide(fit_decay);
988 
989  TGraphErrors *g_doubleRatioMB=new TGraphErrors(gammaYieldMBdoubleratio);
990  g_doubleRatioMB->SetName("g_doubleRatioMB");
991  g_doubleRatioMB->SetMarkerStyle(8);
992  TGraphErrors *g_doubleRatioHT1=new TGraphErrors(gammaYieldHT1doubleratio);
993  g_doubleRatioHT1->SetName("g_doubleRatioHT1");
994  g_doubleRatioHT1->SetMarkerStyle(8);
995  TGraphErrors *g_doubleRatioHT2=new TGraphErrors(gammaYieldHT2doubleratio);
996  g_doubleRatioHT2->SetName("g_doubleRatioHT2");
997  g_doubleRatioHT2->SetMarkerStyle(8);
998 
999  removeThesePoints(g_doubleRatioMB,1);
1000  removeThesePoints(g_doubleRatioHT1,2);
1001  removeThesePoints(g_doubleRatioHT2,3);
1002 
1003  TMultiGraph *m_doubleratio=new TMultiGraph();
1004  m_doubleratio->SetName("m_doubleratio");
1005  m_doubleratio->SetMinimum(.5);
1006  m_doubleratio->SetMaximum(2.75);
1007 
1008  cout<<endl;
1009  g_doubleRatioHT1->Print();
1010  cout<<endl<<endl;
1011  g_doubleRatioHT2->Print();
1012  cout<<endl;
1013 
1014  //m_doubleratio->Add(g_doubleRatioMB,"p");
1015  m_doubleratio->Add(g_doubleRatioHT1,"p");
1016  m_doubleratio->Add(g_doubleRatioHT2,"p");
1017 
1018  g_photonpqcd->SetLineWidth(2);
1019  g_photonpqcd->SetLineColor(2);
1020  g_photonpqcd05->SetLineWidth(2);
1021  g_photonpqcd05->SetLineColor(2);
1022  g_photonpqcd05->SetLineStyle(2);
1023  g_photonpqcd2->SetLineWidth(2);
1024  g_photonpqcd2->SetLineColor(2);
1025  g_photonpqcd2->SetLineStyle(2);
1026 
1027  m_doubleratio->Add(g_photonpqcd,"c");
1028  m_doubleratio->Add(g_photonpqcd05,"c");
1029  m_doubleratio->Add(g_photonpqcd2,"c");
1030 
1031  //appropriate fit to photon pqcd result
1032  TF1 *fitGamma2=new TF1("fitGamma2","1.+[0]*pow(x,[1])",2.,15.);
1033  g_photonpqcd->Fit(fitGamma2,"R0");
1034 
1035  m_doubleratio->Draw("a");
1036 
1037  m_doubleratio->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1038  m_doubleratio->GetYaxis()->SetTitle("1 + #gamma_{dir}/#gamma_{incl}");
1039  m_doubleratio->GetXaxis()->SetRangeUser(2.,16.);
1040 
1041 
1042  TLegend *leg5=new TLegend(.15,.6,.6,.8);
1043  leg5->AddEntry(g_doubleRatioHT1,"hightower-1","p");
1044  leg5->AddEntry(g_doubleRatioHT2,"hightower-2","p");
1045  leg5->AddEntry(g_photonpqcd,"NLO (CTEQ6+KKP) #mu=p_{T}","l");
1046  leg5->AddEntry(g_photonpqcd05,"#mu=2p_{T}, #mu=p_{T}/2","l");
1047  leg5->SetFillColor(0);
1048  leg5->Draw("same");
1049 
1050  c_doubleratio->cd(0);
1051  c_doubleratio->SaveAs((dir+"gammaDoubleRatio.eps").Data());
1052  c_doubleratio->SaveAs((dir+"gammaDoubleRatio.root").Data());
1053 
1054 
1055  TCanvas *c_dirphoton=new TCanvas("c_dirphoton","c_dirphoton",400,300);
1056  gStyle->SetOptStat(0);
1057  c_dirphoton->cd(1);
1058 
1059  TH1F *dirphotonYieldMB=new TH1F(*gammaYieldMBdoubleratio);
1060  TH1F *dirphotonYieldHT1=new TH1F(*gammaYieldHT1doubleratio);
1061  TH1F *dirphotonYieldHT2=new TH1F(*gammaYieldHT2doubleratio);
1062 
1063  TH1F *dirphotonYieldMBnoErr=new TH1F(*dirphotonYieldMB);
1064  for(int i=1;i<=dirphotonYieldMBnoErr->GetNbinsX();i++){
1065  dirphotonYieldMBnoErr->SetBinError(i,0.);
1066  }
1067  TH1F *dirphotonYieldHT1noErr=new TH1F(*dirphotonYieldHT1);
1068  for(int i=1;i<=dirphotonYieldHT1noErr->GetNbinsX();i++){
1069  dirphotonYieldHT1noErr->SetBinError(i,0.);
1070  }
1071  TH1F *dirphotonYieldHT2noErr=new TH1F(*dirphotonYieldHT2);
1072  for(int i=1;i<=dirphotonYieldHT1noErr->GetNbinsX();i++){
1073  dirphotonYieldHT2noErr->SetBinError(i,0.);
1074  }
1075 
1076  TF1 *f_unity=new TF1("f_unity","1.",0.,15.);
1077  dirphotonYieldMB->Add(f_unity,-1.);
1078  dirphotonYieldHT1->Add(f_unity,-1.);
1079  dirphotonYieldHT2->Add(f_unity,-1.);
1080 
1081 
1082  dirphotonYieldMB->Divide(dirphotonYieldMBnoErr);
1083  dirphotonYieldHT1->Divide(dirphotonYieldHT1noErr);
1084  dirphotonYieldHT2->Divide(dirphotonYieldHT2noErr);
1085  dirphotonYieldMB->Multiply(gammaYieldMBratioNoErr);
1086  dirphotonYieldHT1->Multiply(gammaYieldHT1ratioNoErr);
1087  dirphotonYieldHT2->Multiply(gammaYieldHT2ratioNoErr);
1088  dirphotonYieldMB->Multiply(pionXsMBnoErr);
1089  dirphotonYieldHT1->Multiply(pionXsHT1noErr);
1090  dirphotonYieldHT2->Multiply(pionXsHT2noErr);
1091 
1092 
1093  TGraphErrors *g_dirphotonMB=new TGraphErrors(dirphotonYieldMB);
1094  g_dirphotonMB->SetName("g_dirphotonMB");
1095  g_dirphotonMB->SetMarkerStyle(8);
1096  TGraphErrors *g_dirphotonHT1=new TGraphErrors(dirphotonYieldHT1);
1097  g_dirphotonHT1->SetName("g_dirphotonHT1");
1098  g_dirphotonHT1->SetMarkerStyle(8);
1099  TGraphErrors *g_dirphotonHT2=new TGraphErrors(dirphotonYieldHT2);
1100  g_dirphotonHT2->SetName("g_dirphotonHT2");
1101  g_dirphotonHT2->SetMarkerStyle(8);
1102 
1103 
1104  removeThesePoints(g_dirphotonMB,1);
1105  removeThesePoints(g_dirphotonHT1,2);
1106  removeThesePoints(g_dirphotonHT2,3);
1107 
1108  gPad->SetLogy();
1109 
1110  TMultiGraph *m_dirphoton=new TMultiGraph();
1111  m_dirphoton->SetName("m_dirphoton");
1112  m_dirphoton->SetMinimum(1.0e-11);
1113  m_dirphoton->SetMaximum(0.1);
1114 
1115  m_dirphoton->Add(g_dirgamma,"c");
1116  m_dirphoton->Add(g_dirgamma05,"c");
1117  m_dirphoton->Add(g_dirgamma2,"c");
1118 
1119  cout<<"direct photons:"<<endl;
1120  g_dirphotonHT1->Print();
1121  cout<<endl;
1122  g_dirphotonHT2->Print();
1123  cout<<endl;
1124 
1125  m_dirphoton->Add(g_dirphotonHT1,"p");
1126  m_dirphoton->Add(g_dirphotonHT2,"p");
1127 
1128  m_dirphoton->Draw("a");
1129 
1130  m_dirphoton->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1131  m_dirphoton->GetYaxis()->SetTitle("1 - R^{-1}");
1132  m_dirphoton->GetXaxis()->SetRangeUser(2.,16.);
1133 
1134 
1135  TLegend *leg5=new TLegend(.15,.6,.6,.8);
1136  leg5->AddEntry(g_dirphotonHT1,"hightower-1","p");
1137  leg5->AddEntry(g_dirphotonHT2,"hightower-2","p");
1138  leg5->SetFillColor(0);
1139  leg5->Draw("same");
1140 
1141  c_dirphoton->cd(0);
1142  c_dirphoton->SaveAs((dir+"gammaDirPhoton.eps").Data());
1143  c_dirphoton->SaveAs((dir+"gammaDirPhoton.root").Data());
1144 
1145 
1146 
1147 
1148 
1149  return;
1150 }
Generated by   doxygen 1.8.5