StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
getSpectrumDAU.C
1 void getSpectrumDAU(const char *pid)
2 {
3  gStyle->SetOptStat(0);
4  gStyle->SetOptFit(0);
5 
6  gStyle->SetErrorX(0);
7 
8  Bool_t subtractNEUTRONS=kTRUE;
9 
10  //neutron data:
11  //hadron data for Levy function, nucl-ex/0601033:
12  //--> d^2N/(2*pi*pT*N*dpT*dy) = B/((1+((mT - m0)/nT))^n)
13  // {p-dAu; pbar-dAu; p-pp; pbar-pp} and m0 = m_neutron = 1.0 GeV.
14  Double_t B[]={0.3,0.23,0.072,0.061};//->[0]
15  //Double_t eB[]={0.01,0.01,0.005,0.005};
16  Double_t T[]={0.205,0.215,0.179,0.173};//->[1]
17  //Double_t eT[]={0.004,0.005,0.006,0.006};
18  Double_t n[]={11.00,12.55,10.87,10.49};//->[2]
19  //Double_t en[]={0.29,0.41,0.43,0.40};
20  TF1 *f_antineutron=new TF1("f_antineutron","[0]/pow((1.+(sqrt(x*x+1.) - 1.)/([1]*[2])),[2])",0.,15.);
21  f_antineutron->SetParameters(B[1],T[1],n[1]);
22 
23  //get direct gammas pQCD:
24  ifstream pQCDphotons("./datapoints/pQCD_Werner/rhic_cteq6_gamma_inv_sc1.dat");
25  Float_t ppx[100];
26  Float_t ppy[100];
27  Int_t iii=0;
28  cout<<"pqcd photons:"<<endl;
29  while(iii<28){
30  if(!pQCDphotons.good()) break;
31  Float_t dummy=0.;
32  pQCDphotons>>ppx[iii]>>dummy>>dummy>>ppy[iii];
33  ppy[iii]*=1.e-09*7.5/42.;//convert to mb * Nbin / sigma_inel
34  iii++;
35  }
36  TGraph *g_dirgamma=new TGraph(iii,ppx,ppy);
37  //get direct gammas pQCD: scale 0.5*pT
38  ifstream pQCDphotons05("./datapoints/pQCD_Werner/rhic_cteq6_gamma_inv_sc05.dat");
39  Float_t ppx05[100];
40  Float_t ppy05[100];
41  Int_t iii05=0;
42  while(iii05<28){
43  if(!pQCDphotons05.good()) break;
44  Float_t dummy=0.;
45  pQCDphotons05>>ppx05[iii05]>>dummy>>dummy>>ppy05[iii05];
46  ppy05[iii05]*=1.e-09*7.5/42.;//convert to mb
47  iii05++;
48  }
49  TGraph *g_dirgamma05=new TGraph(iii05,ppx05,ppy05);
50  //get direct gammas pQCD: scale 2*pT
51  ifstream pQCDphotons2("./datapoints/pQCD_Werner/rhic_cteq6_gamma_inv_sc2.dat");
52  Float_t ppx2[100];
53  Float_t ppy2[100];
54  Int_t iii2=0;
55  while(iii2<28){
56  if(!pQCDphotons2.good()) break;
57  Float_t dummy=0.;
58  pQCDphotons2>>ppx2[iii2]>>dummy>>dummy>>ppy2[iii2];
59  ppy2[iii2]*=1.e-09*7.5/42.;//convert to mb
60  iii2++;
61  }
62  TGraph *g_dirgamma2=new TGraph(iii2,ppx2,ppy2);
63 
64  //get phenix pions in dau
65  ifstream phenix("./datapoints/phenix.dat");
66  Float_t phex[100];
67  Float_t phey[100];
68  Float_t ephex[100];
69  Float_t ephey[100];
70  Int_t iphe=0;
71  while(iphe<18){
72  if(!phenix.good()) break;
73  phenix>>phex[iphe]>>phey[iphe]>>ephey[iphe];
74  ephex[iphe]=0.;
75  iphe++;
76  }
77  TGraphErrors *phenix_dau=new TGraphErrors(iphe,phex,phey,ephex,ephey);
78  phenix_dau->SetMarkerStyle(24);
79  phenix_dau->SetName("phenix_dau");
80 
81  //get pions KKP scale pT
82  ifstream pQCDpions("./datapoints/pQCD_Werner/klaus_pi0inv_200_kkp_1.dat");
83  Float_t pionx[100];
84  Float_t piony[100];
85  Int_t ipion=0;
86  while(ipion<28){
87  if(!pQCDpions.good()) break;
88  pQCDpions>>pionx[ipion]>>piony[ipion];
89  ipion++;
90  }
91  TGraphErrors *kkp=new TGraphErrors(ipion,pionx,piony);
92  kkp->SetLineColor(54);
93  kkp->SetName("kkp");
94 
95  //get pions KKP scale 0.5*pT
96  ifstream pQCDpions05("./datapoints/pQCD_Werner/klaus_pi0inv_200_kkp_05.dat");
97  Float_t pionx05[100];
98  Float_t piony05[100];
99  Int_t ipion05=0;
100  while(ipion05<28){
101  if(!pQCDpions05.good()) break;
102  pQCDpions05>>pionx05[ipion05]>>piony05[ipion05];
103  ipion05++;
104  }
105  TGraphErrors *kkp05=new TGraphErrors(ipion05,pionx05,piony05);
106  kkp05->SetLineStyle(2);
107  kkp05->SetLineColor(54);
108  kkp05->SetName("kkp05");
109 
110  //get pions KKP scale 2*pT
111  ifstream pQCDpions2("./datapoints/pQCD_Werner/klaus_pi0inv_200_kkp_2.dat");
112  Float_t pionx2[100];
113  Float_t piony2[100];
114  Int_t ipion2=0;
115  while(ipion2<28){
116  if(!pQCDpions2.good()) break;
117  pQCDpions2>>pionx2[ipion2]>>piony2[ipion2];
118  ipion2++;
119  }
120  TGraphErrors *kkp2=new TGraphErrors(ipion2,pionx2,piony2);
121  kkp2->SetLineStyle(2);
122  kkp2->SetLineColor(54);
123  kkp2->SetName("kkp2");
124 
125  TFile *f_decaybg=new TFile("~/MyDecay/gammaDecayDAUSum.root","OPEN");
126  TH1F *h_decaybg=(TH1F*)f_decaybg->Get("gamma");
127  TH1F *h_decaypion=(TH1F*)f_decaybg->Get("gamma_pion");
128  TF1 *fit_decay=new TF1("fit_decay","[0]/pow(x,[1])+[2]",.3,15.);
129  TF1 *fit_piondecay=new TF1(*fit_decay);
130  fit_decay->SetParameters(1.,1.,.5);
131  fit_piondecay->SetParameters(1.,1.,.5);
132  h_decaybg->Fit(fit_decay,"R0");
133  h_decaypion->Fit(fit_piondecay,"R0Q");
134 
135  TCanvas *c_test=new TCanvas();
136  h_decaybg->Draw();
137  fit_decay->Draw("same");
138  c_test->SaveAs("test.eps");
139 
140  //take ratio gamma_direct/pion and divide by gamma/pi and then +1:
141  for(Int_t i=0;i<iii;i++){
142  ppy[i]=ppy[i]/piony[i];
143  ppy[i]=ppy[i]/fit_decay->Eval(ppx[i]);
144  ppy[i]+=1.;
145  ppy05[i]=ppy05[i]/piony05[i];
146  ppy05[i]=ppy05[i]/fit_decay->Eval(ppx05[i]);
147  ppy05[i]+=1.;
148  ppy2[i]=ppy2[i]/piony2[i];
149  ppy2[i]=ppy2[i]/fit_decay->Eval(ppx2[i]);
150  ppy2[i]+=1.;
151  }
152  TGraphErrors *g_photonpqcd=new TGraphErrors(iii,ppx,ppy);
153  g_photonpqcd->SetName("g_photonpqcd");
154  TGraphErrors *g_photonpqcd05=new TGraphErrors(iii05,ppx05,ppy05);
155  g_photonpqcd05->SetName("g_photonpqcd05");
156  TGraphErrors *g_photonpqcd2=new TGraphErrors(iii2,ppx2,ppy2);
157  g_photonpqcd2->SetName("g_photonpqcd2");
158 
159  //set outputfiles
160  TString dir("/star/u/russcher/gamma/analysis/output/dAu/");
161  dir.Append(pid);
162  dir.Append("/");
163  dir.Append(pid);
164  TString psout=dir;
165  TString psout_eta=dir;
166  TString eFile=dir;
167  TString eFileGamma=dir;
168  TString pi0File=dir;
169  TString nbarFile=dir;
170 
171  eFile.Append("pion_eff.root");
172  eFileGamma.Append("gamma_eff.root");
173  pi0File.Append("pi0_dAu.root");
174  nbarFile.Append("antineutron_eff.root");
175  psout_eta.Append("/dev/null");
176  psout.Append("invmassplots.ps");
177 
178  //load *.so
179  gSystem->Load("$HOME/MyEvent/MyEvent");
180  gSystem->Load("$HOME/gamma/analysis/lib/AnaCuts");
181  gSystem->Load("$HOME/gamma/analysis/lib/EventMixer");
182  gSystem->Load("$HOME/gamma/analysis/lib/Pi0Analysis");
183 
184 
185  AnaCuts *cuts=new AnaCuts("dAu");
186 
187  //get inv mass hists
188  TFile f(pi0File.Data(),"OPEN");
189  TH2F *h_mb=new TH2F(*h_minvMB);
190  TH2F *h_ht1=new TH2F(*h_minvHT1);
191  TH2F *h_ht2=new TH2F(*h_minvHT2);
192  TH1F *h_ev=new TH1F(*h_events);
193 
194  //get anti-neutron
195  if(0){
196  TFile *file_nbar=new TFile(nbarFile,"OPEN");
197  TH1F *h_nbarEffMB=new TH1F(*h_effMB);
198  h_nbarEffMB->Sumw2();
199  TH1F *h_nbarEffHT1=new TH1F(*h_effHT1);
200  h_nbarEffHT1->Sumw2();
201  TH1F *h_nbarEffHT2=new TH1F(*h_effHT2);
202  h_nbarEffHT2->Sumw2();
203  }
204  //get prescales
205  int trigger=0;
206  Int_t numberOfMB=0;
207  Int_t numberOfHT1=0;
208  Int_t numberOfHT2=0;
209  for(Int_t i=1;i<=h_ev->GetNbinsX();i++)
210  {
211  trigger=(Int_t)h_ev->GetBinCenter(i);
212  if(trigger&1) numberOfMB+=(Int_t)h_ev->GetBinContent(i);
213  if(trigger&2) numberOfHT1+=(Int_t)h_ev->GetBinContent(i);
214  if(trigger&4) numberOfHT2+=(Int_t)h_ev->GetBinContent(i);
215  }
216 
217  cout<<"number of mb: "<<numberOfMB<<endl;
218  numberOfMB/=0.93;
219  cout<<"nmb after 93% vertex eff.: "<<numberOfMB<<endl;
220 
221  Float_t psMB=383.;
222  Float_t psHT1=9.65;
223  Float_t psHT2=1.;
224 
225  //get efficiencies+acceptance
226  TFile g(eFile.Data(),"OPEN");
227  TH1F *h_emb=new TH1F(*h_effMB);
228  TH1F *h_eht1=new TH1F(*h_effHT1);
229  TH1F *h_eht2=new TH1F(*h_effHT2);
230 
231  //corrections, all multiplicative, in case of
232  //pions:
233  TF1 *pion_cpv_corrMB=new TF1("pion_cpv_corrMB","1./(1.-0.01*(0.4+0.05*x))",0.,15.);
234  TF1 *pion_cpv_corrHT1=new TF1("pion_cpv_corrHT1","1./(1.-0.01*(0.4+0.07*x))",0.,15.);
235  TF1 *pion_cpv_corrHT2=new TF1("pion_cpv_corrHT2","1./(1.-0.01*(0.5+0.06*x))",0.,15.);
236  //photons:
237  TF1 *gamma_cpv_corrMB=new TF1("gamma_cpv_corrMB","1./(1.-0.01*(4.1+0.4*x))",0.,15.);
238  TF1 *gamma_cpv_corrHT1=new TF1("gamma_cpv_corrHT1","1./(1.-0.01*(3.4+0.5*x))",0.,15.);
239  TF1 *gamma_cpv_corrHT2=new TF1("gamma_cpv_corrHT2","1./(1.-0.01*(4.9+0.3*x))",0.,15.);
240  TF1 *gamma_cont_corrMB=new TF1("gamma_cont_corrMB","0.98",0.,15.);
241  TF1 *gamma_cont_corrHT1=new TF1("gamma_cont_corrHT1","0.98",0.,15.);
242  TF1 *gamma_cont_corrHT2=new TF1("gamma_cont_corrHT2","0.965",0.,15.);
243 
244  //missing material
245  //pions:
246  TF1 *pion_conv_corrMB=new TF1("pion_conv_corrMB","1.15",0.,15.);
247  TF1 *pion_conv_corrHT1=new TF1("pion_conv_corrHT1","1.15",0.,15.);
248  TF1 *pion_conv_corrHT2=new TF1("pion_conv_corrHT2","1.15",0.,15.);
249  //photons:
250  TF1 *gamma_conv_corrMB=new TF1("gamma_conv_corrMB","1.08",0.,15.);
251  TF1 *gamma_conv_corrHT1=new TF1("gamma_conv_corrHT1","1.08",0.,15.);
252  TF1 *gamma_conv_corrHT2=new TF1("gamma_conv_corrHT2","1.08",0.,15.);
253 
254  //bin corrections
255  TFile binf("~/BinWidth/bincorrectionsDAU.root","OPEN");
256  TH1F *h_binmb=new TH1F(*h4mb);
257  TH1F *h_binht1=new TH1F(*h4ht1);
258  TH1F *h_binht2=new TH1F(*h4ht2);
259  h_binmb->Sumw2();
260  h_binht1->Sumw2();
261  h_binht2->Sumw2();
262  for(Int_t i=1;i<=h_binmb->GetNbinsX();i++) h_binmb->SetBinError(i,0);
263  for(Int_t i=1;i<=h_binht1->GetNbinsX();i++) h_binht1->SetBinError(i,0);
264  for(Int_t i=1;i<=h_binht2->GetNbinsX();i++) h_binht2->SetBinError(i,0);
265 
266  //get yield
267  Pi0Analysis *pi0=new Pi0Analysis(psout.Data(),psout_eta.Data(),"dAu");
268  pi0->init("/dev/null");
269  TH1F *pionYieldMB=new TH1F(*pi0->getYield(h_mb,"mb"));
270  TH1F *pionYieldHT1=new TH1F(*pi0->getYield(h_ht1,"ht1"));
271  TH1F *pionYieldHT2=new TH1F(*pi0->getYield(h_ht2,"ht2"));
272  pi0->storeCanvases((dir+"canvases.root").Data());
273 
274  cout<<"***************************************"<<endl;
275  cout<<"got yield, dividing by rapidity bite!!!"<<endl;
276  float dy_gamma=cuts->rapidityMaxCUT - cuts->rapidityMinCUT;
277  float dy_pion=cuts->rapPionMaxCUT - cuts->rapPionMinCUT;
278  cout<<"***************************************"<<endl;
279  cout<<endl;
280  cout<<" pion bite is "<<dy_pion<<endl;
281  cout<<" gamma bite is "<<dy_gamma<<endl;
282  cout<<endl;
283  cout<<"***************************************"<<endl;
284 
285  pionYieldMB->Scale(1./dy_pion);
286  pionYieldHT1->Scale(1./dy_pion);
287  pionYieldHT2->Scale(1./dy_pion);
288 
289  //set yield to zero
290  /*
291  for(Int_t i=0;i<pionYieldHT2->GetNbinsX();i++)
292  {
293  if(i<1)
294  {
295  pionYieldMB->SetBinContent(i+1,0);
296  pionYieldMB->SetBinError(i+1,0);
297  }
298  if(i<4)
299  {
300  pionYieldHT1->SetBinContent(i+1,0);
301  pionYieldHT1->SetBinError(i+1,0);
302  }
303  if(i>10)
304  {
305  pionYieldHT1->SetBinContent(i+1,0);
306  pionYieldHT1->SetBinError(i+1,0);
307  }
308  if(i<6)
309  {
310  pionYieldHT2->SetBinContent(i+1,0);
311  pionYieldHT2->SetBinError(i+1,0);
312  }
313  }
314  */
315 
316  //set colors:
317  pionYieldMB->SetMarkerStyle(8);
318  pionYieldMB->SetMarkerSize(1.0);
319  pionYieldHT1->SetMarkerStyle(8);
320  pionYieldHT1->SetMarkerSize(1.0);
321  pionYieldHT1->SetMarkerColor(4);
322  pionYieldHT2->SetMarkerStyle(8);
323  pionYieldHT2->SetMarkerSize(1.0);
324  pionYieldHT2->SetMarkerColor(2);
325 
326  TF1 *scale=new TF1("scale","x",0.,15.);
327 
328  pionYieldMB->SetNameTitle("pionYieldMB","corrected yield MB");
329  pionYieldMB->Divide(h_emb);
330  pionYieldMB->Scale(psMB/(psMB*numberOfMB*2.*TMath::Pi()));
331  pionYieldMB->Divide(scale);
332  pionYieldMB->Multiply(h_binmb);
333  pionYieldMB->Multiply(pion_cpv_corrMB);
334  pionYieldMB->Multiply(pion_conv_corrMB);
335 
336  pionYieldHT1->SetNameTitle("pionYieldHT1","corrected yield HT1");
337  pionYieldHT1->Divide(h_eht1);
338  pionYieldHT1->Scale(psHT1/(psMB*numberOfMB*2.*TMath::Pi()));
339  pionYieldHT1->Divide(scale);
340  pionYieldHT1->Multiply(h_binht1);
341  pionYieldHT1->Multiply(pion_cpv_corrHT1);
342  pionYieldHT1->Multiply(pion_conv_corrHT1);
343 
344  pionYieldHT2->SetNameTitle("pionYieldHT2","corrected yield HT2");
345  pionYieldHT2->Divide(h_eht2);
346  pionYieldHT2->Scale(psHT2/(psMB*numberOfMB*2.*TMath::Pi()));
347  pionYieldHT2->Divide(scale);
348  pionYieldHT2->Multiply(h_binht2);
349  pionYieldHT2->Multiply(pion_cpv_corrHT2);
350  pionYieldHT2->Multiply(pion_conv_corrHT2);
351 
352  //create pion yield for double ratio:
353  TH1F *pionYieldMBratio=new TH1F(*pionYieldMB);
354  TH1F *pionYieldHT1ratio=new TH1F(*pionYieldHT1);
355  TH1F *pionYieldHT2ratio=new TH1F(*pionYieldHT2);
356 
357  TH1F *pionXsMB=new TH1F(*pionYieldMB);
358  TH1F *pionXsHT1=new TH1F(*pionYieldHT1);
359  TH1F *pionXsHT2=new TH1F(*pionYieldHT2);
360 
361  TH1F *pionXsMBnoErr=new TH1F(*pionXsMB);
362  TH1F *pionXsHT1noErr=new TH1F(*pionXsHT1);
363  TH1F *pionXsHT2noErr=new TH1F(*pionXsHT2);
364  for(int i=1;i<=pionXsMBnoErr->GetNbinsX();i++){
365  pionXsMBnoErr->SetBinError(i,0.);
366  }
367  for(int i=1;i<=pionXsHT1noErr->GetNbinsX();i++){
368  pionXsHT1noErr->SetBinError(i,0.);
369  }
370  for(int i=1;i<=pionXsHT2noErr->GetNbinsX();i++){
371  pionXsHT2noErr->SetBinError(i,0.);
372  }
373 
374  TGraphErrors *g_pionXsMB=new TGraphErrors(pionXsMB);
375  g_pionXsMB->SetName("g_pionXsMB");
376  removeThesePoints(g_pionXsMB,1);
377 
378  TGraphErrors *g_pionXsHT1=new TGraphErrors(pionXsHT1);
379  g_pionXsHT1->SetName("g_pionXsHT1");
380  removeThesePoints(g_pionXsHT1,2);
381 
382  TGraphErrors *g_pionXsHT2=new TGraphErrors(pionXsHT2);
383  g_pionXsHT2->SetName("g_pionXsHT2");
384  removeThesePoints(g_pionXsHT2,3);
385 
386  if(0){
387  cout<<endl<<"yield: x y ex ey"<<endl;
388  cout<<"minbias"<<endl;
389  printPoints(g_pionXsMB);
390  cout<<endl<<"hightower-1"<<endl;
391  printPoints(g_pionXsHT1);
392  cout<<endl<<"hightower-2"<<endl;
393  printPoints(g_pionXsHT2);
394  cout<<endl;
395  }
396 
397 
398  TMultiGraph *m_pions_fit=new TMultiGraph();
399  m_pions_fit->SetName("m_pions_fit");
400  m_pions_fit->SetMinimum(1.0e-10);
401  m_pions_fit->SetMaximum(.1);
402 
403 
404  m_pions_fit->Add(g_pionXsMB);
405  m_pions_fit->Add(g_pionXsHT1);
406  m_pions_fit->Add(g_pionXsHT2);
407 
408  TF1 *fitQCD=new TF1("fitQCD",sumpqcd,1.,15.,6);
409  fitQCD->SetParameters(100.,-9.,1.,-8.5,1.,6.);
410  fitQCD->FixParameter(4,2.);
411  m_pions_fit->Fit(fitQCD,"R");
412 
413  bool inclPhenix=true;
414  bool inclSasha=false;
415 
416  TCanvas *compare=new TCanvas("compare","compare;p_{T}:xsec (mb)",600,750);
417  compare->cd();
418 
419  TPad *padt=new TPad("padt","",0.0,0.3,1.,1.0);
420  padt->SetBottomMargin(0.001);
421  padt->SetLeftMargin(0.15);
422  TPad *padb=new TPad("padb","",0.0,0.0,1.,0.3);
423  padb->SetTopMargin(0.001);
424  padb->SetBottomMargin(0.25);
425  padb->SetLeftMargin(0.15);
426 
427  padt->Draw();
428  padt->cd();
429  gPad->SetLogy();
430 
431  TMultiGraph *m_pions=new TMultiGraph();
432  m_pions->SetName("m_pions");
433 
434  if(inclSasha){
435  //
436  }
437  if(inclPhenix){
438  m_pions->Add(phenix_dau);
439  }
440 
441  g_pionXsMB->Print();
442  cout<<endl<<endl;
443  g_pionXsHT1->Print();
444  cout<<endl<<endl;
445  g_pionXsHT2->Print();
446  cout<<endl;
447 
448 
449  m_pions->Add(g_pionXsMB);
450  m_pions->Add(g_pionXsHT1);
451  m_pions->Add(g_pionXsHT2);
452 
453  m_pions->SetMinimum(5.0e-10);
454  m_pions->SetMaximum(.099);
455  m_pions->SetTitle(";p_{T} (GeV/c);invariant yield");
456  m_pions->Draw("ap");
457 
458  //fitQCD->Draw("same");
459 
460  TLegend *leg=new TLegend(.5,.5,.85,.85);
461 
462  if(inclPhenix) leg->AddEntry(phenix_dau,"PHENIX d+Au","p");
463  if(inclSasha) leg->AddEntry(sasha_dau,"sasha's d+Au","p");
464  leg->AddEntry(g_pionXsMB,"d+Au minimum bias","p");
465  leg->AddEntry(g_pionXsHT1,"hightower 1","p");
466  leg->AddEntry(g_pionXsHT2,"hightower 2","p");
467 
468  leg->SetFillColor(0);
469  leg->Draw();
470 
471  compare->cd();
472  padb->Draw();
473  padb->cd();
474 
475  TGraphErrors *phenix_dau2=new TGraphErrors(*phenix_dau);
476  divideGraphWithFunction(phenix_dau2,fitQCD);
477  TGraphErrors *g_pionXsMBcopy=new TGraphErrors(*g_pionXsMB);
478  divideGraphWithFunction(g_pionXsMBcopy,fitQCD);
479  TGraphErrors *g_pionXsHT1copy=new TGraphErrors(*g_pionXsHT1);
480  divideGraphWithFunction(g_pionXsHT1copy,fitQCD);
481  TGraphErrors *g_pionXsHT2copy=new TGraphErrors(*g_pionXsHT2);
482  divideGraphWithFunction(g_pionXsHT2copy,fitQCD);
483 
484  TMultiGraph *m_pions2=new TMultiGraph();
485  //if(inclSasha) m_pions2->Add(sasha_dau2);
486  if(inclPhenix) m_pions2->Add(phenix_dau2);
487  m_pions2->Add(g_pionXsMBcopy);
488  m_pions2->Add(g_pionXsHT1copy);
489  m_pions2->Add(g_pionXsHT2copy);
490  m_pions2->SetMinimum(0.01);
491  m_pions2->SetMaximum(1.99);
492  m_pions2->Draw("ap");
493 
494  compare->SaveAs((dir+"pionyield_dau.eps").Data());
495  compare->SaveAs((dir+"pionyield_dau.root").Data());
496 
497 
498  //********************************************
499  // Get double ratio:
500  //********************************************
501 
502  //pion decay photon eff:
503  TFile gg(eFile.Data(),"OPEN");
504  TH1F *effGammaMB=new TH1F(*h_effDaughtersMB);
505  TH1F *effGammaHT1=new TH1F(*h_effDaughtersHT1);
506  TH1F *effGammaHT2=new TH1F(*h_effDaughtersHT2);
507 
508  //single photon eff:
509  TFile gg_single(eFileGamma.Data(),"OPEN");
510  TH1F *effGammaSingleMB=new TH1F(*h_effMB);
511  TH1F *effGammaSingleHT1=new TH1F(*h_effHT1);
512  TH1F *effGammaSingleHT2=new TH1F(*h_effHT2);
513 
514  //raw neutral clusters:
515  TFile ff(pi0File.Data(),"OPEN");
516  TH1F *gammaYieldMB=new TH1F(*h_gammaMB);
517  TH1F *gammaYieldHT1=new TH1F(*h_gammaHT1);
518  TH1F *gammaYieldHT2=new TH1F(*h_gammaHT2);
519 
520  //divide rap. bite:
521  gammaYieldMB->Scale(1./dy_gamma);
522  gammaYieldHT1->Scale(1./dy_gamma);
523  gammaYieldHT2->Scale(1./dy_gamma);
524 
525 
526  for(Int_t i=1;i<=gammaYieldMB->GetNbinsX();i++){
527  gammaYieldMB->SetBinContent(i,gammaYieldMB->GetBinContent(i)/gammaYieldMB->GetXaxis()->GetBinWidth(i));
528  gammaYieldMB->SetBinError(i,gammaYieldMB->GetBinError(i)/gammaYieldMB->GetXaxis()->GetBinWidth(i));
529  }
530  gammaYieldMB->Scale(psMB/(psMB*numberOfMB*2.*TMath::Pi()));
531  gammaYieldMB->Divide(scale);// ../pT
532  gammaYieldMB->Multiply(h_binmb);
533  gammaYieldMB->Divide(effGammaMB);
534  gammaYieldMB->Multiply(gamma_cpv_corrMB);
535  gammaYieldMB->Multiply(gamma_cont_corrMB);
536  gammaYieldMB->Multiply(gamma_conv_corrMB);
537 
538  gammaYieldMB->SetMarkerStyle(8);
539  gammaYieldMB->SetMarkerSize(1.);
540  TGraphErrors *g_inclPhotonsMB=new TGraphErrors(gammaYieldMB);
541 
542 
543  for(Int_t i=1;i<=gammaYieldHT1->GetNbinsX();i++){
544  gammaYieldHT1->SetBinContent(i,gammaYieldHT1->GetBinContent(i)/gammaYieldHT1->GetXaxis()->GetBinWidth(i));
545  gammaYieldHT1->SetBinError(i,gammaYieldHT1->GetBinError(i)/gammaYieldHT1->GetXaxis()->GetBinWidth(i));
546  }
547  gammaYieldHT1->Scale(psHT1/(psMB*numberOfMB*2.*TMath::Pi()));
548  gammaYieldHT1->Divide(scale);// ../pT
549  gammaYieldHT1->Multiply(h_binht1);
550  gammaYieldHT1->Divide(effGammaHT1);
551  gammaYieldHT1->Multiply(gamma_cpv_corrHT1);
552  gammaYieldHT1->Multiply(gamma_cont_corrHT1);
553  gammaYieldHT1->Multiply(gamma_conv_corrHT1);
554 
555  gammaYieldHT1->SetMarkerStyle(8);
556  gammaYieldHT1->SetMarkerSize(1.);
557  gammaYieldHT1->SetMarkerColor(4);
558  TGraphErrors *g_inclPhotonsHT1=new TGraphErrors(gammaYieldHT1);
559 
560 
561  for(Int_t i=1;i<=gammaYieldHT2->GetNbinsX();i++){
562  gammaYieldHT2->SetBinContent(i,gammaYieldHT2->GetBinContent(i)/gammaYieldHT2->GetXaxis()->GetBinWidth(i));
563  gammaYieldHT2->SetBinError(i,gammaYieldHT2->GetBinError(i)/gammaYieldHT2->GetXaxis()->GetBinWidth(i));
564  }
565  gammaYieldHT2->Scale(psHT2/(psMB*numberOfMB*2.*TMath::Pi()));
566  gammaYieldHT2->Divide(scale);// ../pT
567  gammaYieldHT2->Multiply(h_binht2);
568  gammaYieldHT2->Divide(effGammaHT2);
569  gammaYieldHT2->Multiply(gamma_cpv_corrHT2);
570  gammaYieldHT2->Multiply(gamma_cont_corrHT2);
571  gammaYieldHT2->Multiply(gamma_conv_corrHT2);
572 
573  gammaYieldHT2->SetMarkerStyle(8);
574  gammaYieldHT2->SetMarkerSize(1.);
575  gammaYieldHT2->SetMarkerColor(2);
576  TGraphErrors *g_inclPhotonsHT2=new TGraphErrors(gammaYieldHT2);
577 
578 
579  removeThesePoints(g_inclPhotonsMB,1);
580  removeThesePoints(g_inclPhotonsHT1,2);
581  removeThesePoints(g_inclPhotonsHT2,3);
582 
583 
584  TMultiGraph *m_incl=new TMultiGraph();
585  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}}");
586  m_incl->Add(g_inclPhotonsMB);
587  m_incl->Add(g_inclPhotonsHT1);
588  m_incl->Add(g_inclPhotonsHT2);
589 
590  m_incl->Fit(fitQCD,"R0");
591 
592  m_incl->SetMinimum(1.e-11);
593  m_incl->SetMaximum(10.);
594  TCanvas *c_incl=new TCanvas("c_incl","c_incl",600,400);
595  gPad->SetLogy();
596  m_incl->Draw("ap");
597  c_incl->SaveAs((dir+"inclPhotonYield.eps").Data());
598  c_incl->SaveAs((dir+"inclPhotonYield.root").Data());
599 
600 
601 
602  //get ratio:
603  TH1F *gammaYieldMBratio=new TH1F(*gammaYieldMB);
604  gammaYieldMBratio->SetName("gammaYieldMBratio");
605  TH1F *gammaYieldHT1ratio=new TH1F(*gammaYieldHT1);
606  gammaYieldHT1ratio->SetName("gammaYieldHT1ratio");
607  TH1F *gammaYieldHT2ratio=new TH1F(*gammaYieldHT2);
608  gammaYieldHT2ratio->SetName("gammaYieldHT2ratio");
609 
610  gammaYieldMBratio->Divide(pionYieldMBratio);
611  gammaYieldHT1ratio->Divide(pionYieldHT1ratio);
612  gammaYieldHT2ratio->Divide(pionYieldHT2ratio);
613 
614  //correct gamma over pion ratio, using two efficiencies:
615  getRatio(gammaYieldMBratio,effGammaMB,effGammaSingleMB,fit_piondecay);
616  getRatio(gammaYieldHT1ratio,effGammaHT1,effGammaSingleHT1,fit_piondecay);
617  getRatio(gammaYieldHT2ratio,effGammaHT2,effGammaSingleHT2,fit_piondecay);
618 
619  TH1F *gammaYieldMBratio_incl=new TH1F(*gammaYieldMBratio);
620  TH1F *gammaYieldHT1ratio_incl=new TH1F(*gammaYieldHT1ratio);
621  TH1F *gammaYieldHT2ratio_incl=new TH1F(*gammaYieldHT2ratio);
622 
623 
624  TH1F *gammaYieldMBratioNoErr=new TH1F(*gammaYieldMBratio);
625  TH1F *gammaYieldHT1ratioNoErr=new TH1F(*gammaYieldHT1ratio);
626  TH1F *gammaYieldHT2ratioNoErr=new TH1F(*gammaYieldHT2ratio);
627  for(int i=1;i<=gammaYieldMBratioNoErr->GetNbinsX();i++){
628  gammaYieldMBratioNoErr->SetBinError(i,0.);
629  }
630  for(int i=1;i<=gammaYieldHT1ratioNoErr->GetNbinsX();i++){
631  gammaYieldHT1ratioNoErr->SetBinError(i,0.);
632  }
633  for(int i=1;i<=gammaYieldHT2ratioNoErr->GetNbinsX();i++){
634  gammaYieldHT2ratioNoErr->SetBinError(i,0.);
635  }
636 
637  TGraphErrors *g_ratioMB=new TGraphErrors(gammaYieldMBratio);
638  g_ratioMB->SetName("g_ratioMB");
639  TGraphErrors *g_ratioHT1=new TGraphErrors(gammaYieldHT1ratio);
640  g_ratioHT1->SetName("g_ratioHT1");
641  TGraphErrors *g_ratioHT2=new TGraphErrors(gammaYieldHT2ratio);
642  g_ratioHT2->SetName("g_ratioHT2");
643 
644 
645  removeThesePoints(g_ratioMB,1);
646  removeThesePoints(g_ratioHT1,2);
647  removeThesePoints(g_ratioHT2,3);
648 
649  TCanvas *c_ratio=new TCanvas("c_ratio","c_ratio",400,200);
650 
651  TMultiGraph *m_ratio=new TMultiGraph("m_ratio","d+Au 2003;p_{T};#gamma/#pi^{0}");
652  m_ratio->Add(g_ratioMB);
653  m_ratio->Add(g_ratioHT1);
654  m_ratio->Add(g_ratioHT2);
655 
656  m_ratio->Draw("ap");
657  m_ratio->SetMinimum(.0);
658  m_ratio->SetMaximum(2.);
659 
660  TLegend *leg3=new TLegend(.35,.65,.65,.85);
661  leg3->AddEntry(g_ratioMB,"minimum bias","p");
662  leg3->AddEntry(g_ratioHT1,"hightower 1","p");
663  leg3->AddEntry(g_ratioHT2,"hightower 2","p");
664  leg3->AddEntry(fit_decay,"decay background (total)","l");
665  leg3->AddEntry(fit_piondecay,"decay background (#pi^{0})","l");
666  leg3->SetFillColor(0);
667  leg3->Draw("same");
668 
669  fit_decay->SetLineColor(13);
670  fit_decay->SetLineWidth(1);
671  fit_decay->SetLineColor(1);
672  fit_decay->Draw("same");
673  fit_piondecay->SetLineColor(13);
674  fit_piondecay->SetLineWidth(1);
675  fit_piondecay->SetLineStyle(2);
676  fit_piondecay->SetLineColor(1);
677  fit_piondecay->Draw("same");
678 
679 
680  c_ratio->SaveAs((dir+"gammaOverPion.eps").Data());
681  c_ratio->SaveAs((dir+"gammaOverPion.root").Data());
682 
683 
684  //create fully corrected incl. photons:
685  gammaYieldMBratio_incl->Multiply(pionYieldMBratio);
686  gammaYieldHT1ratio_incl->Multiply(pionYieldHT1ratio);
687  gammaYieldHT2ratio_incl->Multiply(pionYieldHT2ratio);
688 
689  TGraphErrors *g_incl_corrMB=new TGraphErrors(gammaYieldMBratio_incl);
690  TGraphErrors *g_incl_corrHT1=new TGraphErrors(gammaYieldHT1ratio_incl);
691  TGraphErrors *g_incl_corrHT2=new TGraphErrors(gammaYieldHT2ratio_incl);
692 
693  TCanvas *c_incl_corr=new TCanvas("c_incl_corr","c_incl_corr",400,300);
694  gPad->SetLogy();
695  TMultiGraph *m_incl_corr=new TMultiGraph();
696  m_incl_corr->Add(g_incl_corrMB);
697  m_incl_corr->Add(g_incl_corrHT1);
698  m_incl_corr->Add(g_incl_corrHT2);
699 
700  m_incl_corr->SetMinimum(1.e-11);
701  m_incl_corr->SetMaximum(1.);
702 
703  m_incl_corr->Draw("apX");
704  c_incl_corr->SaveAs((dir+"inclPhotonYieldCorr.eps").Data());
705  c_incl_corr->SaveAs((dir+"inclPhotonYieldCorr.root").Data());
706 
707  TCanvas *c_doubleratio=new TCanvas("c_doubleratio","c_doubleratio",400,300);
708  gStyle->SetOptStat(0);
709  c_doubleratio->cd(1);
710 
711  TH1F *gammaYieldMBdoubleratio=new TH1F(*gammaYieldMBratio);
712  TH1F *gammaYieldHT1doubleratio=new TH1F(*gammaYieldHT1ratio);
713  TH1F *gammaYieldHT2doubleratio=new TH1F(*gammaYieldHT2ratio);
714 
715  gammaYieldMBdoubleratio->Divide(fit_decay);
716  gammaYieldHT1doubleratio->Divide(fit_decay);
717  gammaYieldHT2doubleratio->Divide(fit_decay);
718 
719  TGraphErrors *g_doubleRatioMB=new TGraphErrors(gammaYieldMBdoubleratio);
720  g_doubleRatioMB->SetName("g_doubleRatioMB");
721  g_doubleRatioMB->SetMarkerStyle(8);
722  TGraphErrors *g_doubleRatioHT1=new TGraphErrors(gammaYieldHT1doubleratio);
723  g_doubleRatioHT1->SetName("g_doubleRatioHT1");
724  g_doubleRatioHT1->SetMarkerStyle(8);
725  TGraphErrors *g_doubleRatioHT2=new TGraphErrors(gammaYieldHT2doubleratio);
726  g_doubleRatioHT2->SetName("g_doubleRatioHT2");
727  g_doubleRatioHT2->SetMarkerStyle(8);
728 
729  removeThesePoints(g_doubleRatioMB,1);
730  removeThesePoints(g_doubleRatioHT1,2);
731  removeThesePoints(g_doubleRatioHT2,3);
732 
733  TMultiGraph *m_doubleratio=new TMultiGraph();
734  m_doubleratio->SetMinimum(.5);
735  m_doubleratio->SetMaximum(2.75);
736 
737  cout<<endl;
738  g_doubleRatioHT1->Print();
739  cout<<endl<<endl;
740  g_doubleRatioHT2->Print();
741  cout<<endl;
742 
743  //m_doubleratio->Add(g_doubleRatioMB,"p");
744  m_doubleratio->Add(g_doubleRatioHT1,"p");
745  m_doubleratio->Add(g_doubleRatioHT2,"p");
746 
747  g_photonpqcd->SetLineWidth(2);
748  g_photonpqcd->SetLineColor(2);
749  g_photonpqcd05->SetLineWidth(2);
750  g_photonpqcd05->SetLineColor(2);
751  g_photonpqcd05->SetLineStyle(2);
752  g_photonpqcd2->SetLineWidth(2);
753  g_photonpqcd2->SetLineColor(2);
754  g_photonpqcd2->SetLineStyle(2);
755 
756  m_doubleratio->Add(g_photonpqcd,"c");
757  m_doubleratio->Add(g_photonpqcd05,"c");
758  m_doubleratio->Add(g_photonpqcd2,"c");
759 
760  //appropriate fit to photon pqcd result
761  //TF1 *fitGamma2=new TF1("fitGamma2","1.+[0]*pow(x,[1])",2.,15.);
762 
763  m_doubleratio->Draw("a");
764 
765  m_doubleratio->GetXaxis()->SetTitle("p_{T} (GeV/c)");
766  m_doubleratio->GetYaxis()->SetTitle("1 + #gamma_{dir}/#gamma_{incl}");
767  m_doubleratio->GetXaxis()->SetRangeUser(2.,16.);
768 
769 
770  TLegend *leg5=new TLegend(.15,.6,.6,.8);
771  leg5->AddEntry(g_doubleRatioHT1,"hightower-1","p");
772  leg5->AddEntry(g_doubleRatioHT2,"hightower-2","p");
773  leg5->AddEntry(g_photonpqcd,"NLO (CTEQ6+KKP) #mu=p_{T}","l");
774  leg5->AddEntry(g_photonpqcd05,"#mu=2p_{T}, #mu=p_{T}/2","l");
775  leg5->SetFillColor(0);
776  leg5->Draw("same");
777 
778  c_doubleratio->cd(0);
779  c_doubleratio->SaveAs((dir+"gammaDoubleRatio.eps").Data());
780  c_doubleratio->SaveAs((dir+"gammaDoubleRatio.root").Data());
781 
782 
783  TCanvas *c_dirphoton=new TCanvas("c_dirphoton","c_dirphoton",400,300);
784  gStyle->SetOptStat(0);
785  c_dirphoton->cd(1);
786 
787  TH1F *dirphotonYieldMB=new TH1F(*gammaYieldMBdoubleratio);
788  TH1F *dirphotonYieldHT1=new TH1F(*gammaYieldHT1doubleratio);
789  TH1F *dirphotonYieldHT2=new TH1F(*gammaYieldHT2doubleratio);
790 
791  TH1F *dirphotonYieldMBnoErr=new TH1F(*dirphotonYieldMB);
792  for(int i=1;i<=dirphotonYieldMBnoErr->GetNbinsX();i++){
793  dirphotonYieldMBnoErr->SetBinError(i,0.);
794  }
795  TH1F *dirphotonYieldHT1noErr=new TH1F(*dirphotonYieldHT1);
796  for(int i=1;i<=dirphotonYieldHT1noErr->GetNbinsX();i++){
797  dirphotonYieldHT1noErr->SetBinError(i,0.);
798  }
799  TH1F *dirphotonYieldHT2noErr=new TH1F(*dirphotonYieldHT2);
800  for(int i=1;i<=dirphotonYieldHT1noErr->GetNbinsX();i++){
801  dirphotonYieldHT2noErr->SetBinError(i,0.);
802  }
803 
804  TF1 *f_unity=new TF1("f_unity","1.",0.,15.);
805  dirphotonYieldMB->Add(f_unity,-1.);
806  dirphotonYieldHT1->Add(f_unity,-1.);
807  dirphotonYieldHT2->Add(f_unity,-1.);
808 
809 
810  dirphotonYieldMB->Divide(dirphotonYieldMBnoErr);
811  dirphotonYieldHT1->Divide(dirphotonYieldHT1noErr);
812  dirphotonYieldHT2->Divide(dirphotonYieldHT2noErr);
813  dirphotonYieldMB->Multiply(gammaYieldMBratioNoErr);
814  dirphotonYieldHT1->Multiply(gammaYieldHT1ratioNoErr);
815  dirphotonYieldHT2->Multiply(gammaYieldHT2ratioNoErr);
816  dirphotonYieldMB->Multiply(pionXsMBnoErr);
817  dirphotonYieldHT1->Multiply(pionXsHT1noErr);
818  dirphotonYieldHT2->Multiply(pionXsHT2noErr);
819 
820 
821  TGraphErrors *g_dirphotonMB=new TGraphErrors(dirphotonYieldMB);
822  g_dirphotonMB->SetName("g_dirphotonMB");
823  g_dirphotonMB->SetMarkerStyle(8);
824  TGraphErrors *g_dirphotonHT1=new TGraphErrors(dirphotonYieldHT1);
825  g_dirphotonHT1->SetName("g_dirphotonHT1");
826  g_dirphotonHT1->SetMarkerStyle(8);
827  TGraphErrors *g_dirphotonHT2=new TGraphErrors(dirphotonYieldHT2);
828  g_dirphotonHT2->SetName("g_dirphotonHT2");
829  g_dirphotonHT2->SetMarkerStyle(8);
830 
831 
832  removeThesePoints(g_dirphotonMB,1);
833  removeThesePoints(g_dirphotonHT1,2);
834  removeThesePoints(g_dirphotonHT2,3);
835 
836  gPad->SetLogy();
837 
838  TMultiGraph *m_dirphoton=new TMultiGraph();
839  m_dirphoton->SetName("m_dirphoton");
840  m_dirphoton->SetMinimum(1.0e-11);
841  m_dirphoton->SetMaximum(0.1);
842 
843  m_dirphoton->Add(g_dirgamma,"c");
844  m_dirphoton->Add(g_dirgamma05,"c");
845  m_dirphoton->Add(g_dirgamma2,"c");
846 
847  cout<<"direct photons:"<<endl;
848  g_dirphotonHT1->Print();
849  cout<<endl;
850  g_dirphotonHT2->Print();
851  cout<<endl;
852 
853  m_dirphoton->Add(g_dirphotonHT1,"p");
854  m_dirphoton->Add(g_dirphotonHT2,"p");
855 
856  m_dirphoton->Draw("a");
857 
858  m_dirphoton->GetXaxis()->SetTitle("p_{T} (GeV/c)");
859  m_dirphoton->GetYaxis()->SetTitle("1 - R^{-1}");
860  m_dirphoton->GetXaxis()->SetRangeUser(2.,16.);
861 
862 
863  TLegend *leg5=new TLegend(.15,.6,.6,.8);
864  leg5->AddEntry(g_dirphotonHT1,"hightower-1","p");
865  leg5->AddEntry(g_dirphotonHT2,"hightower-2","p");
866  leg5->SetFillColor(0);
867  leg5->Draw("same");
868 
869  c_dirphoton->cd(0);
870  c_dirphoton->SaveAs((dir+"gammaDirPhoton.eps").Data());
871  c_dirphoton->SaveAs((dir+"gammaDirPhoton.root").Data());
872 
873 
874 
875 
876 
877  return;
878 }
Generated by   doxygen 1.8.5