StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
PhotonAnalysisUtil.cxx
1 #include <TMath.h>
2 #include <TGraphErrors.h>
3 #include <TH1F.h>
4 #include <TH2F.h>
5 #include <TStyle.h>
6 #include <TF1.h>
7 #include <TCanvas.h>
8 #include <TFile.h>
9 #include <TPad.h>
10 #include <TLegend.h>
11 #include <TMultiGraph.h>
12 #include <TColor.h>
13 #include <assert.h>
14 
15 #include <fstream>
16 using namespace std;
17 
18 #include <StEmcPool/StPhotonCommon/PhotonAnalysisSettings.h>
19 
20 #include "macros/systematics/systematics_rda.h"
21 #include "macros/systematics/systematics_dau.h"
22 #include "macros/systematics/systematics_pp.h"
23 
24 #include "AnaCuts.h"
25 #include "Pi0Analysis.h"
26 
27 #include "PhotonAnalysisUtil.h"
28 
29 ClassImp(PhotonAnalysisUtil)
30 
31 double pqcdlowpt(double *x,double *par)
32 {
33  double ret=par[0]*TMath::Power(1.+x[0],par[1]);
34  return ret;
35 }
36 double pqcdhighpt(double *x,double *par)
37 {
38  double ret=par[0]*TMath::Power(1.+x[0],par[1]);
39  return ret;
40 }
41 double sumpqcd(double *x,double *par)
42 {
43  //transition (SW=1 is old)
44  //params: 2x low, 2x hi, 2x trans
45  double SW=1. - 1./(1.+TMath::Exp(TMath::Abs(par[4])*(x[0]-par[5])));
46  double ret=SW*pqcdhighpt(x,&par[2]);
47  ret+=(1.-SW)*pqcdlowpt(x,par);
48  //ret=pqcdlowpt(x,par);
49  return ret;
50 }
51 void printPoints(TGraphErrors *g)
52 {
53  for(Int_t i=0;i<g->GetN();i++){
54  double x,y,ey;
55  g->GetPoint(i,x,y);
56  ey=g->GetErrorY(i);
57  cout<<x<<" "<<y<<" "<<ey<<endl;
58  }
59 }
60 void divideGraphWithFunction(TGraphErrors *graph,TF1 *func)
61 {
62  for(int i=0;i<graph->GetN();i++){
63  double x=0.;
64  double y=0.;
65  graph->GetPoint(i,x,y);
66  double ey=graph->GetErrorY(i);
67  graph->SetPoint(i,x,y/func->Eval(x));
68  graph->SetPointError(i,0.,ey/func->Eval(x));
69  }
70 }
71 void divideGraphWithGraph(TGraphErrors *graph,TGraph *denom)
72 {
73  for(int i=0;i<graph->GetN();i++){
74  double x=0.;
75  double y=0.;
76  graph->GetPoint(i,x,y);
77  double ey=graph->GetErrorY(i);
78  graph->SetPoint(i,x,y/denom->Eval(x));
79  graph->SetPointError(i,0.,ey/denom->Eval(x));
80  }
81 }
82 void getRatio(TH1F* h_ratio,TH1F* h_eff,TH1F* h_effSingle,TF1 *f_piondecay)
83 {
84  //copy, don't want to mess with the pointers:
85  TH1F *h_eff_copy=new TH1F(*h_eff);
86  for(int i=1;i<=h_eff->GetNbinsX();i++){
87  h_eff->SetBinError(i,0.);
88  }
89  TH1F *h_effSingle_copy=new TH1F(*h_effSingle);
90  for(int i=1;i<=h_effSingle->GetNbinsX();i++){
91  h_effSingle_copy->SetBinError(i,0.);
92  }
93  //correct for single photon fraction:
94  h_ratio->Add(f_piondecay,-1.);
95  for(int i=1;i<=h_ratio->GetNbinsX();i++){
96  //cout<<i<<endl;
97  //if(h_ratio->GetBinContent(i)<0.) cout<<"negative"<<endl;
98  }
99  h_eff_copy->Divide(h_effSingle_copy);
100  h_ratio->Multiply(h_eff_copy);
101  h_ratio->Add(f_piondecay);
102 }
103 void removeThesePoints(TGraphErrors *g,int trig)
104 {
105  if(trig==1){
106  g->RemovePoint(0);
107  g->RemovePoint(0);
108  //additional for pp:
109  g->RemovePoint(g->GetN()-1);
110  }
111  if(trig==2){
112  g->RemovePoint(0);
113  g->RemovePoint(0);
114  g->RemovePoint(0);
115  g->RemovePoint(0);
116  g->RemovePoint(g->GetN()-1);
117  g->RemovePoint(g->GetN()-1);
118  }
119  if(trig==3){
120  g->RemovePoint(0);
121  g->RemovePoint(0);
122  g->RemovePoint(0);
123  g->RemovePoint(0);
124  g->RemovePoint(0);
125  }
126 }
127 
128 void getPhotonSpectrum(const PhotonAnalysisSettings &settings) {
129  gStyle->SetOptStat(0);
130  gStyle->SetOptFit(0);
131 
132  gStyle->SetErrorX(0);
133 
134  //bool subtractNEUTRONS=kTRUE;
135 
136  //neutron data:
137  //hadron data for Levy function, nucl-ex/0601033:
138  //--> d^2N/(2*pi*pT*N*dpT*dy) = B/((1+((mT - m0)/nT))^n)
139  // {p-dAu; pbar-dAu; p-pp; pbar-pp} and m0 = m_neutron = 1.0 GeV.
140  double B[]={0.3,0.23,0.072,0.061};//->[0]
141  //double eB[]={0.01,0.01,0.005,0.005};
142  double T[]={0.205,0.215,0.179,0.173};//->[1]
143  //double eT[]={0.004,0.005,0.006,0.006};
144  double n[]={11.00,12.55,10.87,10.49};//->[2]
145  //double en[]={0.29,0.41,0.43,0.40};
146  double BR=35.8/63.9;//->p + pi- (63.9%) : ->n + pi0 (35.8%)
147  double c_feeddown=0.8+0.2*BR;
148  double CPVloss=1;//0.98;//for minbias
149  TF1 *f_nbar=new TF1("f_nbar","[3]*[4]*[0]/TMath::Power((1.+(sqrt(x*x+1.) - 1.)/([1]*[2])),[2])",0.,15.);
150 assert(f_nbar);
151  if (settings.name == "pp05") f_nbar->SetParameters(B[3],T[3],n[3],c_feeddown,CPVloss);
152  if (settings.name == "dAu") f_nbar->SetParameters(B[1],T[1],n[1],c_feeddown,CPVloss);
153 
154  //get direct gammas pQCD:
155  ifstream pQCDphotons(TString(settings.input_datapoints_dir + "/pQCD_Werner/rhic_cteq6_gamma_inv_sc1.dat").Data());
156  float ppx[100];
157  float ppy[100];
158  Int_t iii=0;
159  cout<<"pqcd photons:"<<endl;
160  while(iii<28){
161  if(!pQCDphotons.good()) break;
162  float dummy=0.;
163  pQCDphotons>>ppx[iii]>>dummy>>dummy>>ppy[iii];
164  if (settings.name == "pp05") ppy[iii]*=1.e-09;//convert to mb
165  if (settings.name == "dAu") ppy[iii]*=1.e-09*2.21e3*7.5/42.0; //convert to mb * sigma_dAu * Nbin / sigma_inel
166  iii++;
167  }
168  TGraph *g_dirgamma=new TGraph(iii,ppx,ppy);
169 assert(g_dirgamma);
170  //get direct gammas pQCD: scale 0.5*pT
171  ifstream pQCDphotons05(TString(settings.input_datapoints_dir + "/pQCD_Werner/rhic_cteq6_gamma_inv_sc05.dat").Data());
172  float ppx05[100];
173  float ppy05[100];
174  Int_t iii05=0;
175  while(iii05<28){
176  if(!pQCDphotons05.good()) break;
177  float dummy=0.;
178  pQCDphotons05>>ppx05[iii05]>>dummy>>dummy>>ppy05[iii05];
179  if (settings.name == "pp05") ppy05[iii05]*=1.e-09;//convert to mb
180  if (settings.name == "dAu") ppy05[iii05]*=1.e-09*2.21e3*7.5/42.0; //convert to mb * sigma_dAu * Nbin / sigma_inel
181  iii05++;
182  }
183  TGraph *g_dirgamma05=new TGraph(iii05,ppx05,ppy05);
184 assert(g_dirgamma05);
185  //get direct gammas pQCD: scale 2*pT
186  ifstream pQCDphotons2(TString(settings.input_datapoints_dir + "/pQCD_Werner/rhic_cteq6_gamma_inv_sc2.dat").Data());
187  float ppx2[100];
188  float ppy2[100];
189  Int_t iii2=0;
190  while(iii2<28){
191  if(!pQCDphotons2.good()) break;
192  float dummy=0.;
193  pQCDphotons2>>ppx2[iii2]>>dummy>>dummy>>ppy2[iii2];
194  if (settings.name == "pp05") ppy2[iii2]*=1.e-09;//convert to mb
195  if (settings.name == "dAu") ppy2[iii2]*=1.e-09*2.21e3*7.5/42.0; //convert to mb * sigma_dAu * Nbin / sigma_inel
196  iii2++;
197  }
198  TGraph *g_dirgamma2=new TGraph(iii2,ppx2,ppy2);
199 assert(g_dirgamma2);
200  //get phenix pions in pp or dAu
201  TString phenixFile;
202  if (settings.name == "pp05") phenixFile = settings.input_datapoints_dir + "/phenix_xsec_pp.dat";
203  if (settings.name == "dAu") phenixFile = settings.input_datapoints_dir + "/phenix.dat";
204  ifstream phenix(phenixFile.Data());
205  float phex[100];
206  float phey[100];
207  float ephey[100];
208  Int_t iphe=0;
209  while(iphe<16){
210  if(!phenix.good()) break;
211  if (settings.name == "pp05") phenix>>phex[iphe]>>phey[iphe];
212  if (settings.name == "dAu") phenix>>phex[iphe]>>phey[iphe]>>ephey[iphe];
213  //is in mb already!!
214  iphe++;
215  }
216  TGraphErrors *phenix_pp=new TGraphErrors(iphe,phex,phey);
217 assert(phenix_pp);
218  phenix_pp->SetMarkerStyle(24);
219  phenix_pp->SetLineWidth(6);
220  phenix_pp->SetName("phenix");
221 
222  //frank's spin2006 pions:
223  ifstream frank(TString(settings.input_datapoints_dir + "/frank_pp05_new.dat").Data());
224  float frx[100];
225  float fry[100];
226  float frex[100];
227  float frey[100];
228  Int_t ifr=0;
229  while(ifr<12){
230  if(!frank.good()) break;
231  frank>>frx[ifr]>>fry[ifr]>>frey[ifr];
232  frex[ifr]=0.;
233  ifr++;
234  }
235  TGraphErrors *frank_pp=new TGraphErrors(ifr,frx,fry,frex,frey);
236 assert(frank_pp);
237  frank_pp->SetMarkerStyle(25);
238  frank_pp->SetMarkerColor(4);
239  frank_pp->SetName("frank_pp");
240 
241  //sasha's pions:
242  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 };
243  float xe_pp2005_sasha[]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
244  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};
245  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};
246 
247  TGraphErrors *sasha_pp05=new TGraphErrors(13,x_pp2005_sasha,y_pp2005_sasha,xe_pp2005_sasha,ye_pp2005_sasha);
248 assert(sasha_pp05);
249  sasha_pp05->SetMarkerStyle(8);
250  sasha_pp05->SetMarkerColor(TColor::GetColor(8,80,8));
251  sasha_pp05->SetName("sasha_pp05");
252 
253  //get pions KKP scale pT
254  ifstream pQCDpions(TString(settings.input_datapoints_dir + "/pQCD_Werner/klaus_pi0inv_200_kkp_1.dat").Data());
255  float pionx[100];
256  float piony[100];
257  Int_t ipion=0;
258  while(ipion<28){
259  if(!pQCDpions.good()) break;
260  pQCDpions>>pionx[ipion]>>piony[ipion];
261  if (settings.name == "pp05") piony[iii]*=1.e-09;//convert to mb
262  if (settings.name == "dAu") piony[iii]*=1.e-09*2.21e3*7.5/42.0; //convert to mb * sigma_dAu * Nbin / sigma_inel
263  ipion++;
264  }
265  TGraphErrors *kkp=new TGraphErrors(ipion,pionx,piony);
266 assert(kkp);
267  kkp->SetLineColor(54);
268  kkp->SetName("kkp");
269 
270  //get pions KKP scale 0.5*pT
271  ifstream pQCDpions05(TString(settings.input_datapoints_dir + "/pQCD_Werner/klaus_pi0inv_200_kkp_05.dat").Data());
272  float pionx05[100];
273  float piony05[100];
274  Int_t ipion05=0;
275  while(ipion05<28){
276  if(!pQCDpions05.good()) break;
277  pQCDpions05>>pionx05[ipion05]>>piony05[ipion05];
278  if (settings.name == "pp05") piony05[iii05]*=1.e-09;//convert to mb
279  if (settings.name == "dAu") piony05[iii05]*=1.e-09*2.21e3*7.5/42.0; //convert to mb * sigma_dAu * Nbin / sigma_inel
280  ipion05++;
281  }
282  TGraphErrors *kkp05=new TGraphErrors(ipion05,pionx05,piony05);
283 assert(kkp05);
284  kkp05->SetLineStyle(2);
285  kkp05->SetLineColor(54);
286  kkp05->SetName("kkp05");
287 
288  //get pions KKP scale 2*pT
289  ifstream pQCDpions2(TString(settings.input_datapoints_dir + "/pQCD_Werner/klaus_pi0inv_200_kkp_2.dat").Data());
290  float pionx2[100];
291  float piony2[100];
292  Int_t ipion2=0;
293  while(ipion2<28){
294  if(!pQCDpions2.good()) break;
295  pQCDpions2>>pionx2[ipion2]>>piony2[ipion2];
296  if (settings.name == "pp05") piony2[iii2]*=1.e-09;//convert to mb
297  if (settings.name == "dAu") piony2[iii2]*=1.e-09*2.21e3*7.5/42.0; //convert to mb * sigma_dAu * Nbin / sigma_inel
298  ipion2++;
299  }
300  TGraphErrors *kkp2=new TGraphErrors(ipion2,pionx2,piony2);
301 assert(kkp2);
302  kkp2->SetLineStyle(2);
303  kkp2->SetLineColor(54);
304  kkp2->SetName("kkp2");
305 
306  TFile f_decaybg(settings.input_decaybackground_file,"READ");
307  const TH1F *h_decaybg_read=(const TH1F*)f_decaybg.Get("gamma");
308  const TH1F *h_decaypion_read=(const TH1F*)f_decaybg.Get("gamma_pion");
309 assert(h_decaybg_read);
310 assert(h_decaypion_read);
311  TH1F *h_decaybg = new TH1F(*h_decaybg_read);
312  TH1F *h_decaypion = new TH1F(*h_decaypion_read);
313 assert(h_decaybg);
314 assert(h_decaypion);
315  f_decaybg.Close();
316  h_decaybg_read = 0;
317  h_decaypion_read = 0;
318 
319  TF1 *fit_decay=new TF1("fit_decay","[0]/TMath::Power(x,[1])+[2]",.3,15.);
320  TF1 *fit_piondecay=new TF1(*fit_decay);
321 assert(fit_decay);
322 assert(fit_piondecay);
323  fit_decay->SetParameters(1.,1.,.5);
324  fit_piondecay->SetParameters(1.,1.,.5);
325  h_decaybg->Fit(fit_decay,"R0Q");
326  h_decaypion->Fit(fit_piondecay,"R0Q");
327 
328  //take ratio gamma_direct/pion and divide by gamma/pi and then +1:
329  for(Int_t i=0;i<iii;i++){
330  ppy[i]=ppy[i]/piony[i];
331  ppy[i]=ppy[i]/fit_decay->Eval(ppx[i]);
332  ppy[i]+=1.;
333  ppy05[i]=ppy05[i]/piony05[i];
334  ppy05[i]=ppy05[i]/fit_decay->Eval(ppx05[i]);
335  ppy05[i]+=1.;
336  ppy2[i]=ppy2[i]/piony2[i];
337  ppy2[i]=ppy2[i]/fit_decay->Eval(ppx2[i]);
338  ppy2[i]+=1.;
339  }
340  TGraphErrors *g_photonpqcd=new TGraphErrors(iii,ppx,ppy);
341 assert(g_photonpqcd);
342  g_photonpqcd->SetName("g_photonpqcd");
343  TGraphErrors *g_photonpqcd05=new TGraphErrors(iii05,ppx05,ppy05);
344 assert(g_photonpqcd05);
345  g_photonpqcd05->SetName("g_photonpqcd05");
346  TGraphErrors *g_photonpqcd2=new TGraphErrors(iii2,ppx2,ppy2);
347 assert(g_photonpqcd2);
348  g_photonpqcd2->SetName("g_photonpqcd2");
349 
350  AnaCuts *cuts=new AnaCuts(settings.name.Data());
351 assert(cuts);
352 
353  //get inv mass hists
354  TFile f(settings.input_pion_file.Data(),"READ");
355  const TH2F *h_minvMB = (const TH2F*)f.Get("h_minvMB");
356  const TH2F *h_minvHT1 = (const TH2F*)f.Get("h_minvHT1");
357  const TH2F *h_minvHT2 = (const TH2F*)f.Get("h_minvHT2");
358  const TH1F *h_events = (const TH1F*)f.Get("h_events");
359 assert(h_minvMB);
360 assert(h_minvHT1);
361 assert(h_minvHT2);
362 assert(h_events);
363  TH2F *h_mb=new TH2F(*h_minvMB);
364  TH2F *h_ht1=new TH2F(*h_minvHT1);
365  TH2F *h_ht2=new TH2F(*h_minvHT2);
366  TH1F *h_ev=new TH1F(*h_events);
367 assert(h_mb);
368 assert(h_ht1);
369 assert(h_ht2);
370 assert(h_ev);
371  f.Close();
372  h_minvMB = 0;
373  h_minvHT1 = 0;
374  h_minvHT2 = 0;
375  //get anti-neutron
376  TFile file_nbar(settings.input_nbareff_file.Data(),"READ");
377  const TH1F *h_effnbarMB = (const TH1F*)file_nbar.Get("h_effMB");
378  const TH1F *h_effnbarHT1 = (const TH1F*)file_nbar.Get("h_effHT1");
379  const TH1F *h_effnbarHT2 = (const TH1F*)file_nbar.Get("h_effHT2");
380 assert(h_effnbarMB);
381 assert(h_effnbarHT1);
382 assert(h_effnbarHT2);
383  TH1F *h_nbarEffMB=new TH1F(*h_effnbarMB);
384  TH1F *h_nbarEffHT1=new TH1F(*h_effnbarHT1);
385  TH1F *h_nbarEffHT2=new TH1F(*h_effnbarHT2);
386 assert(h_nbarEffMB);
387 assert(h_nbarEffHT1);
388 assert(h_nbarEffHT2);
389  h_nbarEffMB->Sumw2();
390  h_nbarEffHT1->Sumw2();
391  h_nbarEffHT2->Sumw2();
392  file_nbar.Close();
393  h_effnbarMB = 0;
394  h_effnbarHT1 = 0;
395  h_effnbarHT2 = 0;
396 
397  //get prescales
398  int trigger=0;
399  Float_t numberOfMB=0;
400  Float_t numberOfHT1=0;
401  Float_t numberOfHT2=0;
402  for(Int_t i=1;i<=h_ev->GetNbinsX();i++)
403  {
404  trigger=(Int_t)h_ev->GetBinCenter(i);
405  if(trigger&1) numberOfMB+=(Int_t)h_ev->GetBinContent(i);
406  if(trigger&2) numberOfHT1+=(Int_t)h_ev->GetBinContent(i);
407  if(trigger&4) numberOfHT2+=(Int_t)h_ev->GetBinContent(i);
408  }
409 
410  cout<<"number of mb: "<<numberOfMB<<endl;
411  if (settings.name == "dAu") {
412  numberOfMB/=0.93;
413  cout<<"nmb after 63% vertex eff.: "<<numberOfMB<<endl;
414  }
415 
416  Float_t psMB=1;
417  Float_t psHT1=1;
418  Float_t psHT2=1;
419 
420  if (settings.name == "pp05") {
421  psMB=32742;
422  psHT1=3.89;
423  psHT2=1.;
424  }
425  if (settings.name == "dAu") {
426  psMB=383.;
427  psHT1=9.65;
428  psHT2=1.;
429  }
430 
431  //mb events for hightower:
432  float nMBwithHT=1;
433  if (settings.name == "pp05") {
434  nMBwithHT = 328871;
435  }
436  if (settings.name == "dAu") {
437  nMBwithHT = numberOfMB;
438  }
439  //get efficiencies+acceptance
440  TFile g(settings.input_pioneff_file.Data(),"READ");
441  const TH1F *h_effpionMB = (const TH1F*)g.Get("h_effMB");
442  const TH1F *h_effpionHT1 = (const TH1F*)g.Get("h_effHT1");
443  const TH1F *h_effpionHT2 = (const TH1F*)g.Get("h_effHT2");
444 assert(h_effpionMB);
445 assert(h_effpionHT1);
446 assert(h_effpionHT2);
447  TH1F *h_emb=new TH1F(*h_effpionMB);
448  TH1F *h_eht1=new TH1F(*h_effpionHT1);
449  TH1F *h_eht2=new TH1F(*h_effpionHT2);
450 assert(h_emb);
451 assert(h_eht1);
452 assert(h_eht2);
453  g.Close();
454  h_effpionMB = 0;
455  h_effpionHT1 = 0;
456  h_effpionHT2 = 0;
457 
458  //bin corrections
459  TFile binf(settings.input_binwidth_file,"READ");
460  const TH1F *h4mb = (const TH1F*)binf.Get("h4mb");
461  const TH1F *h4ht1 = (const TH1F*)binf.Get("h4ht1");
462  const TH1F *h4ht2 = (const TH1F*)binf.Get("h4ht2");
463 assert(h4mb);
464 assert(h4ht1);
465 assert(h4ht2);
466  TH1F *h_binmb=new TH1F(*h4mb);
467  TH1F *h_binht1=new TH1F(*h4ht1);
468  TH1F *h_binht2=new TH1F(*h4ht2);
469 assert(h_binmb);
470 assert(h_binht1);
471 assert(h_binht2);
472  binf.Close();
473  h4mb = 0;
474  h4ht1 = 0;
475  h4ht2 = 0;
476 
477  h_binmb->Sumw2();
478  h_binht1->Sumw2();
479  h_binht2->Sumw2();
480  for(Int_t i=1;i<=h_binmb->GetNbinsX();i++) h_binmb->SetBinError(i,0);
481  for(Int_t i=1;i<=h_binht1->GetNbinsX();i++) h_binht1->SetBinError(i,0);
482  for(Int_t i=1;i<=h_binht2->GetNbinsX();i++) h_binht2->SetBinError(i,0);
483 
484  //corrections, all multiplicative, in case of
485  //pions:
486  TF1 *pion_cpv_corrMB=0;
487  TF1 *pion_cpv_corrHT1=0;
488  TF1 *pion_cpv_corrHT2=0;
489  //photons:
490  TF1 *gamma_cpv_corrMB=0;
491  TF1 *gamma_cpv_corrHT1=0;
492  TF1 *gamma_cpv_corrHT2=0;
493  TF1 *gamma_cont_corrMB=0;
494  TF1 *gamma_cont_corrHT1=0;
495  TF1 *gamma_cont_corrHT2=0;
496  //missing material
497  //pions:
498  TF1 *pion_conv_corrMB=0;
499  TF1 *pion_conv_corrHT1=0;
500  TF1 *pion_conv_corrHT2=0;
501  //photons:
502  TF1 *gamma_conv_corrMB=0;
503  TF1 *gamma_conv_corrHT1=0;
504  TF1 *gamma_conv_corrHT2=0;
505  if (settings.name == "pp05") {
506  pion_cpv_corrMB = new TF1("pion_cpv_corrMB","1./(1.-0.01*(0.3+0.0*x))",0.,15.);
507  pion_cpv_corrHT1 = new TF1("pion_cpv_corrHT1","1./(1.-0.01*(-0.1+0.16*x))",0.,15.);
508  pion_cpv_corrHT2 = new TF1("pion_cpv_corrHT2","1./(1.-0.01*(-0.2+0.18*x))",0.,15.);
509  gamma_cpv_corrMB = new TF1("gamma_cpv_corrMB","1./(1.-0.01*(2.8+0.0*x))",0.,15.);
510  gamma_cpv_corrHT1 = new TF1("gamma_cpv_corrHT1","1./(1.-0.01*(0.2+1.1*x))",0.,15.);
511  gamma_cpv_corrHT2 = new TF1("gamma_cpv_corrHT2","1./(1.-0.01*(0.4+1.1*x))",0.,15.);
512  gamma_cont_corrMB = new TF1("gamma_cont_corrMB","0.985",0.,15.);
513  gamma_cont_corrHT1 = new TF1("gamma_cont_corrHT1","0.98",0.,15.);
514  gamma_cont_corrHT2 = new TF1("gamma_cont_corrHT2","0.96",0.,15.);
515  pion_conv_corrMB = new TF1("pion_conv_corrMB","1.1",0.,15.);
516  pion_conv_corrHT1 = new TF1("pion_conv_corrHT1","1.1",0.,15.);
517  pion_conv_corrHT2 = new TF1("pion_conv_corrHT2","1.1",0.,15.);
518  gamma_conv_corrMB = new TF1("gamma_conv_corrMB","1.05",0.,15.);
519  gamma_conv_corrHT1 = new TF1("gamma_conv_corrHT1","1.05",0.,15.);
520  gamma_conv_corrHT2 = new TF1("gamma_conv_corrHT2","1.05",0.,15.);
521  }
522  if (settings.name == "dAu") {
523  pion_cpv_corrMB = new TF1("pion_cpv_corrMB","1./(1.-0.01*(0.4+0.05*x))",0.,15.);
524  pion_cpv_corrHT1 = new TF1("pion_cpv_corrHT1","1./(1.-0.01*(0.4+0.07*x))",0.,15.);
525  pion_cpv_corrHT2 = new TF1("pion_cpv_corrHT2","1./(1.-0.01*(0.5+0.06*x))",0.,15.);
526  gamma_cpv_corrMB = new TF1("gamma_cpv_corrMB","1./(1.-0.01*(4.1+0.4*x))",0.,15.);
527  gamma_cpv_corrHT1 = new TF1("gamma_cpv_corrHT1","1./(1.-0.01*(3.4+0.5*x))",0.,15.);
528  gamma_cpv_corrHT2 = new TF1("gamma_cpv_corrHT2","1./(1.-0.01*(4.9+0.3*x))",0.,15.);
529  gamma_cont_corrMB = new TF1("gamma_cont_corrMB","0.98",0.,15.);
530  gamma_cont_corrHT1 = new TF1("gamma_cont_corrHT1","0.98",0.,15.);
531  gamma_cont_corrHT2 = new TF1("gamma_cont_corrHT2","0.965",0.,15.);
532  pion_conv_corrMB = new TF1("pion_conv_corrMB","1.15",0.,15.);
533  pion_conv_corrHT1 = new TF1("pion_conv_corrHT1","1.15",0.,15.);
534  pion_conv_corrHT2 = new TF1("pion_conv_corrHT2","1.15",0.,15.);
535  gamma_conv_corrMB = new TF1("gamma_conv_corrMB","1.08",0.,15.);
536  gamma_conv_corrHT1 = new TF1("gamma_conv_corrHT1","1.08",0.,15.);
537  gamma_conv_corrHT2 = new TF1("gamma_conv_corrHT2","1.08",0.,15.);
538  }
539 assert(pion_cpv_corrMB);
540 assert(pion_cpv_corrHT1);
541 assert(pion_cpv_corrHT2);
542 assert(gamma_cpv_corrMB);
543 assert(gamma_cpv_corrHT1);
544 assert(gamma_cpv_corrHT2);
545 assert(gamma_cont_corrMB);
546 assert(gamma_cont_corrHT1);
547 assert(gamma_cont_corrHT2);
548 assert(pion_conv_corrMB);
549 assert(pion_conv_corrHT1);
550 assert(pion_conv_corrHT2);
551 assert(gamma_conv_corrMB);
552 assert(gamma_conv_corrHT1);
553 assert(gamma_conv_corrHT2);
554 
555  //get yield
556  Pi0Analysis *pi0=new Pi0Analysis(settings.output_invmassplots_file.Data(),settings.output_invmassplotseta_file.Data(),settings.name.Data());
557 assert(pi0);
558  pi0->init(settings.output_pionhistograms_file.Data());
559  TH1F *pionYieldMB=new TH1F(*pi0->getYield(h_mb,"mb"));
560  TH1F *pionYieldHT1=new TH1F(*pi0->getYield(h_ht1,"ht1"));
561  TH1F *pionYieldHT2=new TH1F(*pi0->getYield(h_ht2,"ht2"));
562 assert(pionYieldMB);
563 assert(pionYieldHT1);
564 assert(pionYieldHT2);
565  pi0->storeCanvases(settings.output_pioncanvases_file.Data());
566 
567 
568  cout<<"***************************************"<<endl;
569  cout<<"got yield, dividing by rapidity bite!!!"<<endl;
570  float dy_gamma=cuts->rapidityMaxCUT - cuts->rapidityMinCUT;
571  float dy_pion=cuts->rapPionMaxCUT - cuts->rapPionMinCUT;
572  cout<<"***************************************"<<endl;
573  cout<<endl;
574  cout<<" pion bite is "<<dy_pion<<endl;
575  cout<<" gamma bite is "<<dy_gamma<<endl;
576  cout<<endl;
577  cout<<"***************************************"<<endl;
578 
579  pionYieldMB->Scale(1./dy_pion);
580  pionYieldHT1->Scale(1./dy_pion);
581  pionYieldHT2->Scale(1./dy_pion);
582 
583  //set yield to zero
584  /*
585  for(Int_t i=0;i<pionYieldHT2->GetNbinsX();i++)
586  {
587  if(i<1)
588  {
589  pionYieldMB->SetBinContent(i+1,0);
590  pionYieldMB->SetBinError(i+1,0);
591  }
592  if(i<4)
593  {
594  pionYieldHT1->SetBinContent(i+1,0);
595  pionYieldHT1->SetBinError(i+1,0);
596  }
597  if(i>12)
598  {
599  pionYieldHT1->SetBinContent(i+1,0);
600  pionYieldHT1->SetBinError(i+1,0);
601  }
602  if(i<6)
603  {
604  pionYieldHT2->SetBinContent(i+1,0);
605  pionYieldHT2->SetBinError(i+1,0);
606  }
607  }
608  */
609 
610  //set colors:
611  pionYieldMB->SetMarkerStyle(8);
612  pionYieldMB->SetMarkerSize(1.0);
613  pionYieldHT1->SetMarkerStyle(8);
614  pionYieldHT1->SetMarkerSize(1.0);
615  pionYieldHT1->SetMarkerColor(4);
616  pionYieldHT2->SetMarkerStyle(8);
617  pionYieldHT2->SetMarkerSize(1.0);
618  pionYieldHT2->SetMarkerColor(2);
619 
620  TF1 *scale=new TF1("scale","x",0.,15.);
621 assert(scale);
622 
623  pionYieldMB->SetNameTitle("pionYieldMB","corrected yield MB");
624  pionYieldMB->Divide(h_emb);
625  pionYieldMB->Scale(psMB/(psMB*numberOfMB*2.*TMath::Pi()));
626  pionYieldMB->Divide(scale);
627  pionYieldMB->Multiply(h_binmb);
628  pionYieldMB->Multiply(pion_cpv_corrMB);
629  pionYieldMB->Multiply(pion_conv_corrMB);
630 
631  pionYieldHT1->SetNameTitle("pionYieldHT1","corrected yield HT1");
632  pionYieldHT1->Divide(h_eht1);
633  pionYieldHT1->Scale(psHT1/(psMB*nMBwithHT*2.*TMath::Pi()));
634  pionYieldHT1->Divide(scale);
635  pionYieldHT1->Multiply(h_binht1);
636  pionYieldHT1->Multiply(pion_cpv_corrHT1);
637  pionYieldHT1->Multiply(pion_conv_corrHT1);
638 
639  pionYieldHT2->SetNameTitle("pionYieldHT2","corrected yield HT2");
640  pionYieldHT2->Divide(h_eht2);
641  pionYieldHT2->Scale(psHT2/(psMB*nMBwithHT*2.*TMath::Pi()));
642  pionYieldHT2->Divide(scale);
643  pionYieldHT2->Multiply(h_binht2);
644  pionYieldHT2->Multiply(pion_cpv_corrHT2);
645  pionYieldHT2->Multiply(pion_conv_corrHT2);
646 
647  //create pion yield for double ratio:
648  TH1F *pionYieldMBratio=new TH1F(*pionYieldMB);
649  TH1F *pionYieldHT1ratio=new TH1F(*pionYieldHT1);
650  TH1F *pionYieldHT2ratio=new TH1F(*pionYieldHT2);
651 assert(pionYieldMBratio);
652 assert(pionYieldHT1ratio);
653 assert(pionYieldHT2ratio);
654 
655  TH1F *pionXsMB=new TH1F(*pionYieldMB);
656  TH1F *pionXsHT1=new TH1F(*pionYieldHT1);
657  TH1F *pionXsHT2=new TH1F(*pionYieldHT2);
658 assert(pionXsMB);
659 assert(pionXsHT1);
660 assert(pionXsHT2);
661  if (settings.name=="pp05") {
662  pionXsMB->Scale((26.1/0.85));//times xs_bbc/bbc_eff
663  pionXsHT1->Scale((26.1/0.85));
664  pionXsHT2->Scale((26.1/0.85));
665  }
666  if (settings.name=="dAu") {
667  pionXsMB->Scale(2.21e3);//times cross section
668  pionXsHT1->Scale(2.21e3);
669  pionXsHT2->Scale(2.21e3);
670  }
671 
672 
673  TH1F *pionXsMBnoErr=new TH1F(*pionXsMB);
674  TH1F *pionXsHT1noErr=new TH1F(*pionXsHT1);
675  TH1F *pionXsHT2noErr=new TH1F(*pionXsHT2);
676 assert(pionXsMBnoErr);
677 assert(pionXsHT1noErr);
678 assert(pionXsHT2noErr);
679  for(int i=1;i<=pionXsMBnoErr->GetNbinsX();i++){
680  pionXsMBnoErr->SetBinError(i,0.);
681  }
682  for(int i=1;i<=pionXsHT1noErr->GetNbinsX();i++){
683  pionXsHT1noErr->SetBinError(i,0.);
684  }
685  for(int i=1;i<=pionXsHT2noErr->GetNbinsX();i++){
686  pionXsHT2noErr->SetBinError(i,0.);
687  }
688 
689  TGraphErrors *g_pionXsMB=new TGraphErrors(pionXsMB);
690 assert(g_pionXsMB);
691  g_pionXsMB->SetName("g_pionXsMB");
692  removeThesePoints(g_pionXsMB,1);
693 
694  TGraphErrors *g_pionXsHT1=new TGraphErrors(pionXsHT1);
695 assert(g_pionXsHT1);
696  g_pionXsHT1->SetName("g_pionXsHT1");
697  removeThesePoints(g_pionXsHT1,2);
698 
699  TGraphErrors *g_pionXsHT2=new TGraphErrors(pionXsHT2);
700 assert(g_pionXsHT2);
701  g_pionXsHT2->SetName("g_pionXsHT2");
702  removeThesePoints(g_pionXsHT2,3);
703 
704  if(1){
705  cout<<endl<<"xsec: x y ex ey"<<endl;
706  cout<<"minbias"<<endl;
707  printPoints(g_pionXsMB);
708  cout<<endl<<"hightower-1"<<endl;
709  printPoints(g_pionXsHT1);
710  cout<<endl<<"hightower-2"<<endl;
711  printPoints(g_pionXsHT2);
712  cout<<endl;
713  }
714 
715 
716  TMultiGraph *m_pions_fit=new TMultiGraph();
717 assert(m_pions_fit);
718  m_pions_fit->SetName("m_pions_fit");
719  m_pions_fit->SetMinimum(5.0e-9);
720  m_pions_fit->SetMaximum(0.99);
721 
722 
723  m_pions_fit->Add(g_pionXsMB);
724  m_pions_fit->Add(g_pionXsHT1);
725  m_pions_fit->Add(g_pionXsHT2);
726 
727  TF1 *fitQCD=new TF1("fitQCD",sumpqcd,1.,15.,6);
728 assert(fitQCD);
729  if (settings.name=="pp05") fitQCD->SetParameters(600.,-8.2,4.,-8.5,2.,2.);
730  if (settings.name=="dAu") fitQCD->SetParameters(100.,-9.,1.,-8.5,1.,6.);
731  fitQCD->FixParameter(4,2.);
732 cout << "111" << endl;
733 // m_pions_fit->Fit(fitQCD,"RV");
734 cout << "222" << endl;
735 
736  bool inclPhenix=true;
737  bool inclFrank=true;
738  bool inclSasha=true;
739  bool inclPqcd=true;
740 
741  TCanvas *compare=new TCanvas("compare","compare;p_{T}:xsec (mb)",600,750);
742 assert(compare);
743  compare->Divide(1, 2, 0, 0);
744 
745  compare->cd(1);
746  gPad->SetLogy();
747 
748  TMultiGraph *m_pions=new TMultiGraph();
749 assert(m_pions);
750  m_pions->SetName("m_pions");
751 
752  if(inclPqcd){
753  m_pions->Add(kkp,"c");
754  m_pions->Add(kkp05,"c");
755  m_pions->Add(kkp2,"c");
756  }
757  if(inclSasha){
758  m_pions->Add(sasha_pp05);
759  }
760  if(inclFrank){
761  m_pions->Add(frank_pp);
762  }
763  if(inclPhenix){
764  m_pions->Add(phenix_pp);
765  }
766 
767  //g_pionXsMB->Print();
768  //cout<<endl<<endl;
769  //g_pionXsHT1->Print();
770  //cout<<endl<<endl;
771  //g_pionXsHT2->Print();
772  //cout<<endl;
773 
774  m_pions->Add(g_pionXsMB);
775  m_pions->Add(g_pionXsHT1);
776  m_pions->Add(g_pionXsHT2);
777 
778  m_pions->SetMinimum(1.0e-9);
779  m_pions->SetMaximum(1.);
780 
781  m_pions->Draw("ap");
782 
783  //fitQCD->Draw("same");
784 
785  m_pions->GetXaxis()->SetLabelSize(0.);
786  m_pions->GetXaxis()->SetTitle("p_{T} [GeV/c]");
787  m_pions->GetYaxis()->SetTitle("Ed^{3}#sigma/dp^{3} [mb GeV^{-2} c^{3}]");
788 
789  TLegend *leg=new TLegend(.5,.5,.85,.85);
790 assert(leg);
791 
792  if(inclPhenix) leg->AddEntry(phenix_pp,"PHENIX p+p","p");
793  if(inclFrank) leg->AddEntry(frank_pp,"STAR preliminary (upd.)","p");
794  if(inclSasha) leg->AddEntry(sasha_pp05,"O.Grebenyuk p+p","p");
795  leg->AddEntry(g_pionXsMB,"p+p minimum bias","p");
796  leg->AddEntry(g_pionXsHT1,"hightower 1","p");
797  leg->AddEntry(g_pionXsHT2,"hightower 2","p");
798  if(inclPqcd){
799  leg->AddEntry(kkp,"kkp + CTEQ6m, #mu=p_{T}","l");
800  leg->AddEntry(kkp2,"#mu=2p_{T},p_{T}/2","l");
801  leg->Draw("same");
802  }
803 
804  leg->SetFillColor(0);
805  leg->Draw();
806 
807 /*
808  compare->cd();
809  padb->Draw();
810  padb->cd();
811 */
812  compare->cd(2);
813 
814  TGraphErrors *sasha_pp05_overPqcd=new TGraphErrors(*sasha_pp05);
815 assert(sasha_pp05_overPqcd);
816  divideGraphWithGraph(sasha_pp05_overPqcd,kkp);
817  TGraphErrors *phenix_pp05_overPqcd=new TGraphErrors(*phenix_pp);
818 assert(phenix_pp05_overPqcd);
819  divideGraphWithGraph(phenix_pp05_overPqcd,kkp);
820  TGraphErrors *g_pionXsMB_overPqcd=new TGraphErrors(*g_pionXsMB);
821 assert(g_pionXsMB_overPqcd);
822  divideGraphWithGraph(g_pionXsMB_overPqcd,kkp);
823  TGraphErrors *g_pionXsHT1_overPqcd=new TGraphErrors(*g_pionXsHT1);
824 assert(g_pionXsHT1_overPqcd);
825  divideGraphWithGraph(g_pionXsHT1_overPqcd,kkp);
826  TGraphErrors *g_pionXsHT2_overPqcd=new TGraphErrors(*g_pionXsHT2);
827 assert(g_pionXsHT2_overPqcd);
828  divideGraphWithGraph(g_pionXsHT2_overPqcd,kkp);
829  TGraphErrors *frank_pp05_overPqcd=new TGraphErrors(*frank_pp);
830 assert(frank_pp05_overPqcd);
831  divideGraphWithGraph(frank_pp05_overPqcd,kkp);
832 
833  TGraphErrors *kkp05_ratio=new TGraphErrors(*kkp05);
834 assert(kkp05_ratio);
835  divideGraphWithGraph(kkp05_ratio,kkp);
836  TGraphErrors *kkp_ratio=new TGraphErrors(*kkp);
837 assert(kkp_ratio);
838  divideGraphWithGraph(kkp_ratio,kkp);
839  TGraphErrors *kkp2_ratio=new TGraphErrors(*kkp2);
840 assert(kkp2_ratio);
841  divideGraphWithGraph(kkp2_ratio,kkp);
842 
843  //systematic errors:
844  TGraphErrors *g_pionXsMB_sys=new TGraphErrors(*g_pionXsMB_overPqcd);
845 assert(g_pionXsMB_sys);
846  set_sys_pp_pion(g_pionXsMB_sys);
847  TGraphErrors *g_pionXsHT1_sys=new TGraphErrors(*g_pionXsHT1_overPqcd);
848 assert(g_pionXsHT1_sys);
849  set_sys_pp_pion(g_pionXsHT1_sys);
850  TGraphErrors *g_pionXsHT2_sys=new TGraphErrors(*g_pionXsHT2_overPqcd);
851 assert(g_pionXsHT2_sys);
852  set_sys_pp_pion(g_pionXsHT2_sys);
853 
854  TMultiGraph *m_pions_over_pqcd=new TMultiGraph();
855 assert(m_pions_over_pqcd);
856 
857  m_pions_over_pqcd->SetMinimum(0.0);
858  m_pions_over_pqcd->SetMaximum(2.5);
859 
860  if(inclPqcd){
861  m_pions_over_pqcd->Add(kkp05_ratio,"c");
862  m_pions_over_pqcd->Add(kkp_ratio,"c");
863  m_pions_over_pqcd->Add(kkp2_ratio,"c");
864  }
865  m_pions_over_pqcd->Add(g_pionXsMB_overPqcd);
866  m_pions_over_pqcd->Add(g_pionXsHT1_overPqcd);
867  m_pions_over_pqcd->Add(g_pionXsHT2_overPqcd);
868  //m_pions_over_pqcd->Add(g_pionXsMB_sys,"c");
869  //m_pions_over_pqcd->Add(g_pionXsHT1_sys,"c");
870  //m_pions_over_pqcd->Add(g_pionXsHT2_sys,"c");
871  if(inclPhenix) m_pions_over_pqcd->Add(phenix_pp05_overPqcd);
872  if(inclFrank) m_pions_over_pqcd->Add(frank_pp05_overPqcd);
873  if(inclSasha) m_pions_over_pqcd->Add(sasha_pp05_overPqcd);
874 
875  m_pions_over_pqcd->Draw("ap");
876  m_pions_over_pqcd->GetXaxis()->SetTitle("p_{T} [GeV/c]");
877  m_pions_over_pqcd->GetYaxis()->SetTitle("Data / theory");
878 
879  compare->SaveAs(settings.output_pionxsec_file.Data());
880 
881  TMultiGraph *m_pions_over_fit=new TMultiGraph();
882 assert(m_pions_over_fit);
883  m_pions_over_fit->SetMinimum(0.01);
884  m_pions_over_fit->SetMaximum(1.99);
885 
886  TGraphErrors *g_pionXsMB_overFit=new TGraphErrors(*g_pionXsMB);
887 assert(g_pionXsMB_overFit);
888  divideGraphWithFunction(g_pionXsMB_overFit,fitQCD);
889  TGraphErrors *g_pionXsHT1_overFit=new TGraphErrors(*g_pionXsHT1);
890 assert(g_pionXsHT1_overFit);
891  divideGraphWithFunction(g_pionXsHT1_overFit,fitQCD);
892  TGraphErrors *g_pionXsHT2_overFit=new TGraphErrors(*g_pionXsHT2);
893 assert(g_pionXsHT2_overFit);
894  divideGraphWithFunction(g_pionXsHT2_overFit,fitQCD);
895 
896  m_pions_over_fit->Add(g_pionXsMB_overFit);
897  m_pions_over_fit->Add(g_pionXsHT1_overFit);
898  m_pions_over_fit->Add(g_pionXsHT2_overFit);
899 
900  m_pions_over_fit->Draw("ap");
901 
902  compare->SaveAs(settings.output_pionxsecoverfit_file.Data());
903 
904 
905 
906  TCanvas *compare2=new TCanvas("compare2","compare2;p_{T};yield divided by fit",600,300);
907 assert(compare2);
908  compare2->cd();
909 
910  //divide by fit:
911  TGraphErrors *frank_pp2=new TGraphErrors(*frank_pp);
912 assert(frank_pp2);
913  divideGraphWithFunction(frank_pp2,fitQCD);
914  TGraphErrors *sasha_pp2=new TGraphErrors(*sasha_pp05);
915 assert(sasha_pp2);
916  divideGraphWithFunction(sasha_pp2,fitQCD);
917  TGraphErrors *phenix_pp2=new TGraphErrors(*phenix_pp);
918 assert(phenix_pp2);
919  divideGraphWithFunction(phenix_pp2,fitQCD);
920  TGraphErrors *g_pionXsMBcopy=new TGraphErrors(*g_pionXsMB);
921 assert(g_pionXsMBcopy);
922  divideGraphWithFunction(g_pionXsMBcopy,fitQCD);
923  TGraphErrors *g_pionXsHT1copy=new TGraphErrors(*g_pionXsHT1);
924 assert(g_pionXsHT1copy);
925  divideGraphWithFunction(g_pionXsHT1copy,fitQCD);
926  TGraphErrors *g_pionXsHT2copy=new TGraphErrors(*g_pionXsHT2);
927 assert(g_pionXsHT2copy);
928  divideGraphWithFunction(g_pionXsHT2copy,fitQCD);
929 
930  inclPhenix=true;
931  inclFrank=false;
932  inclSasha=false;
933  inclPqcd=false;
934 
935  TMultiGraph *m_pions2=new TMultiGraph();
936 assert(m_pions2);
937  m_pions2->SetName("m_pions2");
938  if(inclSasha) m_pions2->Add(sasha_pp2);
939  if(inclFrank) m_pions2->Add(frank_pp2);
940  if(inclPhenix) m_pions2->Add(phenix_pp2);
941 
942  m_pions2->Add(g_pionXsMBcopy);
943  m_pions2->Add(g_pionXsHT1copy);
944  m_pions2->Add(g_pionXsHT2copy);
945 
946  m_pions2->SetMinimum(0.000001);
947  m_pions2->SetMaximum(3.);
948  m_pions2->Draw("ap");
949 
950  TLegend *legg=new TLegend(.25,.55,.65,.85);
951 assert(legg);
952  legg->AddEntry(g_pionXsMBcopy,"minimum bias","p");
953  legg->AddEntry(g_pionXsHT1copy,"hightower 1","p");
954  legg->AddEntry(g_pionXsHT2copy,"hightower 2","p");
955  if(inclFrank) legg->AddEntry(frank_pp2,"Frank's p+p (upd.)","p");
956  if(inclSasha) legg->AddEntry(sasha_pp2,"Sasha's p+p","p");
957  if(inclPhenix) legg->AddEntry(phenix_pp,"PHENIX p+p","p");
958  legg->Draw("same");
959  legg->SetFillColor(0);
960 
961  compare2->cd(0);
962  compare2->SaveAs(settings.output_pionxsecratio_file.Data());
963 
964 
965  //********************************************
966  // Get double ratio:
967  //********************************************
968 
969  //pion decay photon eff:
970  TFile gg(settings.input_pioneff_file.Data(),"READ");
971  const TH1F *h_effDaughtersMB = (const TH1F*)gg.Get("h_effDaughtersMB");
972  const TH1F *h_effDaughtersHT1 = (const TH1F*)gg.Get("h_effDaughtersHT1");
973  const TH1F *h_effDaughtersHT2 = (const TH1F*)gg.Get("h_effDaughtersHT2");
974 assert(h_effDaughtersMB);
975 assert(h_effDaughtersHT1);
976 assert(h_effDaughtersHT2);
977  TH1F *effGammaMB=new TH1F(*h_effDaughtersMB);
978  TH1F *effGammaHT1=new TH1F(*h_effDaughtersHT1);
979  TH1F *effGammaHT2=new TH1F(*h_effDaughtersHT2);
980 assert(effGammaMB);
981 assert(effGammaHT1);
982 assert(effGammaHT2);
983  gg.Close();
984  h_effDaughtersMB = 0;
985  h_effDaughtersHT1 = 0;
986  h_effDaughtersHT2 = 0;
987 
988  //single photon eff:
989  TFile gg_single(settings.input_gammaeff_file.Data(),"READ");
990  const TH1F *h_effgammaMB = (const TH1F*)gg_single.Get("h_effMB");
991  const TH1F *h_effgammaHT1 = (const TH1F*)gg_single.Get("h_effHT1");
992  const TH1F *h_effgammaHT2 = (const TH1F*)gg_single.Get("h_effHT2");
993 assert(h_effgammaMB);
994 assert(h_effgammaHT1);
995 assert(h_effgammaHT2);
996  TH1F *effGammaSingleMB=new TH1F(*h_effgammaMB);
997  TH1F *effGammaSingleHT1=new TH1F(*h_effgammaHT1);
998  TH1F *effGammaSingleHT2=new TH1F(*h_effgammaHT2);
999 assert(effGammaSingleMB);
1000 assert(effGammaSingleHT1);
1001 assert(effGammaSingleHT2);
1002  gg_single.Close();
1003  h_effgammaMB = 0;
1004  h_effgammaHT1 = 0;
1005  h_effgammaHT2 = 0;
1006 
1007  //raw neutral clusters:
1008  TFile ff(settings.input_pion_file.Data(),"READ");
1009  const TH1F *h_gammaMB = (const TH1F*)ff.Get("h_gammaMB");
1010  const TH1F *h_gammaHT1 = (const TH1F*)ff.Get("h_gammaHT1");
1011  const TH1F *h_gammaHT2 = (const TH1F*)ff.Get("h_gammaHT2");
1012 assert(h_gammaMB);
1013 assert(h_gammaHT1);
1014 assert(h_gammaHT2);
1015  TH1F *gammaYieldMB=new TH1F(*h_gammaMB);
1016  TH1F *gammaYieldHT1=new TH1F(*h_gammaHT1);
1017  TH1F *gammaYieldHT2=new TH1F(*h_gammaHT2);
1018 assert(gammaYieldMB);
1019 assert(gammaYieldHT1);
1020 assert(gammaYieldHT2);
1021  ff.Close();
1022  h_gammaMB = 0;
1023  h_gammaHT1 = 0;
1024  h_gammaHT2 = 0;
1025 
1026  //divide rap. bite:
1027  gammaYieldMB->Scale(1./dy_gamma);
1028  gammaYieldHT1->Scale(1./dy_gamma);
1029  gammaYieldHT2->Scale(1./dy_gamma);
1030 
1031 
1032  for(Int_t i=1;i<=gammaYieldMB->GetNbinsX();i++){
1033  gammaYieldMB->SetBinContent(i,gammaYieldMB->GetBinContent(i)/gammaYieldMB->GetXaxis()->GetBinWidth(i));
1034  gammaYieldMB->SetBinError(i,gammaYieldMB->GetBinError(i)/gammaYieldMB->GetXaxis()->GetBinWidth(i));
1035  }
1036  gammaYieldMB->Scale(psMB/(psMB*numberOfMB*2.*TMath::Pi()));
1037  gammaYieldMB->Divide(scale);// ../pT
1038  gammaYieldMB->Multiply(h_binmb);
1039  gammaYieldMB->Divide(effGammaMB);
1040  gammaYieldMB->Multiply(gamma_cpv_corrMB);
1041  gammaYieldMB->Multiply(gamma_cont_corrMB);
1042  gammaYieldMB->Multiply(gamma_conv_corrMB);
1043 
1044  gammaYieldMB->SetMarkerStyle(8);
1045  gammaYieldMB->SetMarkerSize(1.);
1046  TGraphErrors *g_inclPhotonsMB=new TGraphErrors(gammaYieldMB);
1047 assert(g_inclPhotonsMB);
1048 
1049  for(Int_t i=1;i<=gammaYieldHT1->GetNbinsX();i++){
1050  gammaYieldHT1->SetBinContent(i,gammaYieldHT1->GetBinContent(i)/gammaYieldHT1->GetXaxis()->GetBinWidth(i));
1051  gammaYieldHT1->SetBinError(i,gammaYieldHT1->GetBinError(i)/gammaYieldHT1->GetXaxis()->GetBinWidth(i));
1052  }
1053  gammaYieldHT1->Scale(psHT1/(psMB*nMBwithHT*2.*TMath::Pi()));
1054  gammaYieldHT1->Divide(scale);// ../pT
1055  gammaYieldHT1->Multiply(h_binht1);
1056  gammaYieldHT1->Divide(effGammaHT1);
1057  gammaYieldHT1->Multiply(gamma_cpv_corrHT1);
1058  gammaYieldHT1->Multiply(gamma_cont_corrHT1);
1059  gammaYieldHT1->Multiply(gamma_conv_corrHT1);
1060 
1061  gammaYieldHT1->SetMarkerStyle(8);
1062  gammaYieldHT1->SetMarkerSize(1.);
1063  gammaYieldHT1->SetMarkerColor(4);
1064  TGraphErrors *g_inclPhotonsHT1=new TGraphErrors(gammaYieldHT1);
1065 assert(g_inclPhotonsHT1);
1066 
1067  for(Int_t i=1;i<=gammaYieldHT2->GetNbinsX();i++){
1068  gammaYieldHT2->SetBinContent(i,gammaYieldHT2->GetBinContent(i)/gammaYieldHT2->GetXaxis()->GetBinWidth(i));
1069  gammaYieldHT2->SetBinError(i,gammaYieldHT2->GetBinError(i)/gammaYieldHT2->GetXaxis()->GetBinWidth(i));
1070  }
1071  gammaYieldHT2->Scale(psHT2/(psMB*nMBwithHT*2.*TMath::Pi()));
1072  gammaYieldHT2->Divide(scale);// ../pT
1073  gammaYieldHT2->Multiply(h_binht2);
1074  gammaYieldHT2->Divide(effGammaHT2);
1075  gammaYieldHT2->Multiply(gamma_cpv_corrHT2);
1076  gammaYieldHT2->Multiply(gamma_cont_corrHT2);
1077  gammaYieldHT2->Multiply(gamma_conv_corrHT2);
1078 
1079  gammaYieldHT2->SetMarkerStyle(8);
1080  gammaYieldHT2->SetMarkerSize(1.);
1081  gammaYieldHT2->SetMarkerColor(2);
1082  TGraphErrors *g_inclPhotonsHT2=new TGraphErrors(gammaYieldHT2);
1083 assert(g_inclPhotonsHT2);
1084 
1085  removeThesePoints(g_inclPhotonsMB,1);
1086  removeThesePoints(g_inclPhotonsHT1,2);
1087  removeThesePoints(g_inclPhotonsHT2,2);
1088 
1089  TMultiGraph *m_incl=new TMultiGraph();
1090 assert(m_incl);
1091  m_incl->SetName("m_incl");
1092  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}}");
1093  m_incl->Add(g_inclPhotonsMB);
1094  m_incl->Add(g_inclPhotonsHT1);
1095  m_incl->Add(g_inclPhotonsHT2);
1096 
1097  //m_incl->Fit(fitQCD,"R0");
1098 
1099  m_incl->SetMinimum(1.e-11);
1100  m_incl->SetMaximum(10.);
1101  TCanvas *c_incl=new TCanvas("c_incl","c_incl",600,400);
1102 assert(c_incl);
1103  gPad->SetLogy();
1104  m_incl->Draw("ap");
1105  c_incl->SaveAs(settings.output_inclphotonyield_file.Data());
1106 
1107 
1108 
1109  //get ratio:
1110  TH1F *gammaYieldMBratio=new TH1F(*gammaYieldMB);
1111  TH1F *gammaYieldHT1ratio=new TH1F(*gammaYieldHT1);
1112  TH1F *gammaYieldHT2ratio=new TH1F(*gammaYieldHT2);
1113 assert(gammaYieldMBratio);
1114 assert(gammaYieldHT1ratio);
1115 assert(gammaYieldHT2ratio);
1116  gammaYieldMBratio->SetName("gammaYieldMBratio");
1117  gammaYieldHT1ratio->SetName("gammaYieldHT1ratio");
1118  gammaYieldHT2ratio->SetName("gammaYieldHT2ratio");
1119 
1120  gammaYieldMBratio->Divide(pionYieldMBratio);
1121  gammaYieldHT1ratio->Divide(pionYieldHT1ratio);
1122  gammaYieldHT2ratio->Divide(pionYieldHT2ratio);
1123 
1124 
1125  //correct gamma over pion ratio, using two efficiencies:
1126  getRatio(gammaYieldMBratio,effGammaMB,effGammaSingleMB,fit_piondecay);
1127  getRatio(gammaYieldHT1ratio,effGammaHT1,effGammaSingleHT1,fit_piondecay);
1128  getRatio(gammaYieldHT2ratio,effGammaHT2,effGammaSingleHT2,fit_piondecay);
1129 
1130  TH1F *gammaYieldMBratio_incl=new TH1F(*gammaYieldMBratio);
1131  TH1F *gammaYieldHT1ratio_incl=new TH1F(*gammaYieldHT1ratio);
1132  TH1F *gammaYieldHT2ratio_incl=new TH1F(*gammaYieldHT2ratio);
1133 assert(gammaYieldMBratio_incl);
1134 assert(gammaYieldHT1ratio_incl);
1135 assert(gammaYieldHT2ratio_incl);
1136 
1137  TH1F *gammaYieldMBratioNoErr=new TH1F(*gammaYieldMBratio);
1138  TH1F *gammaYieldHT1ratioNoErr=new TH1F(*gammaYieldHT1ratio);
1139  TH1F *gammaYieldHT2ratioNoErr=new TH1F(*gammaYieldHT2ratio);
1140 assert(gammaYieldMBratioNoErr);
1141 assert(gammaYieldHT1ratioNoErr);
1142 assert(gammaYieldHT2ratioNoErr);
1143  for(int i=1;i<=gammaYieldMBratioNoErr->GetNbinsX();i++)gammaYieldMBratioNoErr->SetBinError(i,0.);
1144  for(int i=1;i<=gammaYieldHT1ratioNoErr->GetNbinsX();i++)gammaYieldHT1ratioNoErr->SetBinError(i,0.);
1145  for(int i=1;i<=gammaYieldHT2ratioNoErr->GetNbinsX();i++)gammaYieldHT2ratioNoErr->SetBinError(i,0.);
1146 
1147  TGraphErrors *g_ratioMB=new TGraphErrors(gammaYieldMBratio);
1148  TGraphErrors *g_ratioHT1=new TGraphErrors(gammaYieldHT1ratio);
1149  TGraphErrors *g_ratioHT2=new TGraphErrors(gammaYieldHT2ratio);
1150 assert(g_ratioMB);
1151 assert(g_ratioHT1);
1152 assert(g_ratioHT2);
1153  g_ratioMB->SetName("g_ratioMB");
1154  g_ratioHT1->SetName("g_ratioHT1");
1155  g_ratioHT2->SetName("g_ratioHT2");
1156 
1157 
1158  removeThesePoints(g_ratioMB,1);
1159  removeThesePoints(g_ratioHT1,2);
1160  removeThesePoints(g_ratioHT2,3);
1161 
1162 
1163  TCanvas *c_ratio=new TCanvas("c_ratio","c_ratio",400,200);
1164 assert(c_ratio);
1165 
1166  TMultiGraph *m_ratio=new TMultiGraph("m_ratio","p+p 2005;p_{T};#gamma/#pi^{0}");
1167 assert(m_ratio);
1168  m_ratio->Add(g_ratioMB);
1169  m_ratio->Add(g_ratioHT1);
1170  m_ratio->Add(g_ratioHT2);
1171 
1172  m_ratio->Draw("ap");
1173  m_ratio->SetMinimum(.001);
1174  m_ratio->SetMaximum(1.5);
1175 
1176  TLegend *leg3=new TLegend(.35,.65,.65,.85);
1177 assert(leg3);
1178  leg3->AddEntry(g_ratioMB,"minimum bias","p");
1179  leg3->AddEntry(g_ratioHT1,"hightower 1","p");
1180  leg3->AddEntry(g_ratioHT2,"hightower 2","p");
1181  leg3->AddEntry(fit_decay,"decay background (total)","l");
1182  leg3->AddEntry(fit_piondecay,"decay background (#pi^{0})","l");
1183  leg3->SetFillColor(0);
1184  leg3->Draw("same");
1185 
1186  fit_decay->SetLineColor(13);
1187  fit_decay->SetLineWidth(1);
1188  fit_decay->SetLineColor(1);
1189  fit_decay->Draw("same");
1190  fit_piondecay->SetLineColor(13);
1191  fit_piondecay->SetLineWidth(1);
1192  fit_piondecay->SetLineStyle(2);
1193  fit_piondecay->SetLineColor(1);
1194  fit_piondecay->Draw("same");
1195 
1196  c_ratio->SaveAs(settings.output_gammaoverpion_file.Data());
1197 
1198  //create fully corrected incl. photons:
1199  gammaYieldMBratio_incl->Multiply(pionYieldMBratio);
1200  gammaYieldHT1ratio_incl->Multiply(pionYieldHT1ratio);
1201  gammaYieldHT2ratio_incl->Multiply(pionYieldHT2ratio);
1202  //gammaYieldMBratio_incl->Scale(7.5);
1203  //gammaYieldHT1ratio_incl->Scale(7.5);
1204  //gammaYieldHT2ratio_incl->Scale(7.5);
1205 
1206  TGraphErrors *g_incl_corrMB=new TGraphErrors(gammaYieldMBratio_incl);
1207  TGraphErrors *g_incl_corrHT1=new TGraphErrors(gammaYieldHT1ratio_incl);
1208  TGraphErrors *g_incl_corrHT2=new TGraphErrors(gammaYieldHT2ratio_incl);
1209 assert(g_incl_corrMB);
1210 assert(g_incl_corrHT1);
1211 assert(g_incl_corrHT2);
1212 
1213  TCanvas *c_incl_corr=new TCanvas("c_incl_corr","c_incl_corr",400,300);
1214 assert(c_incl_corr);
1215  gPad->SetLogy();
1216  TMultiGraph *m_incl_corr=new TMultiGraph();
1217 assert(m_incl_corr);
1218  m_incl_corr->Add(g_incl_corrMB);
1219  m_incl_corr->Add(g_incl_corrHT1);
1220  m_incl_corr->Add(g_incl_corrHT2);
1221 
1222  m_incl_corr->SetMinimum(1.e-11);
1223  m_incl_corr->SetMaximum(1.);
1224 
1225  m_incl_corr->Draw("apX");
1226  c_incl_corr->SaveAs(settings.output_inclphotonyieldcorr_file.Data());
1227 
1228  TCanvas *c_doubleratio=new TCanvas("c_doubleratio","c_doubleratio",400,300);
1229 assert(c_doubleratio);
1230  gStyle->SetOptStat(0);
1231  c_doubleratio->cd(1);
1232 
1233  TH1F *gammaYieldMBdoubleratio=new TH1F(*gammaYieldMBratio);
1234  TH1F *gammaYieldHT1doubleratio=new TH1F(*gammaYieldHT1ratio);
1235  TH1F *gammaYieldHT2doubleratio=new TH1F(*gammaYieldHT2ratio);
1236 assert(gammaYieldMBdoubleratio);
1237 assert(gammaYieldHT1doubleratio);
1238 assert(gammaYieldHT2doubleratio);
1239 
1240  gammaYieldMBdoubleratio->Divide(fit_decay);
1241  gammaYieldHT1doubleratio->Divide(fit_decay);
1242  gammaYieldHT2doubleratio->Divide(fit_decay);
1243 
1244  TGraphErrors *g_doubleRatioMB=new TGraphErrors(gammaYieldMBdoubleratio);
1245 assert(g_doubleRatioMB);
1246  g_doubleRatioMB->SetName("g_doubleRatioMB");
1247  g_doubleRatioMB->SetMarkerStyle(8);
1248  TGraphErrors *g_doubleRatioHT1=new TGraphErrors(gammaYieldHT1doubleratio);
1249 assert(g_doubleRatioHT1);
1250  g_doubleRatioHT1->SetName("g_doubleRatioHT1");
1251  g_doubleRatioHT1->SetMarkerStyle(8);
1252  TGraphErrors *g_doubleRatioHT2=new TGraphErrors(gammaYieldHT2doubleratio);
1253 assert(g_doubleRatioHT2);
1254  g_doubleRatioHT2->SetName("g_doubleRatioHT2");
1255  g_doubleRatioHT2->SetMarkerStyle(8);
1256 
1257  removeThesePoints(g_doubleRatioMB,1);
1258  removeThesePoints(g_doubleRatioHT1,2);
1259  removeThesePoints(g_doubleRatioHT2,3);
1260 
1261  TMultiGraph *m_doubleratio=new TMultiGraph();
1262 assert(m_doubleratio);
1263  m_doubleratio->SetName("m_doubleratio");
1264  m_doubleratio->SetMinimum(.5);
1265  m_doubleratio->SetMaximum(2.75);
1266 
1267  cout<<endl;
1268  g_doubleRatioHT1->Print();
1269  cout<<endl<<endl;
1270  g_doubleRatioHT2->Print();
1271  cout<<endl;
1272 
1273  //m_doubleratio->Add(g_doubleRatioMB,"p");
1274  m_doubleratio->Add(g_doubleRatioHT1,"p");
1275  m_doubleratio->Add(g_doubleRatioHT2,"p");
1276 
1277  g_photonpqcd->SetLineWidth(2);
1278  g_photonpqcd->SetLineColor(2);
1279  g_photonpqcd05->SetLineWidth(2);
1280  g_photonpqcd05->SetLineColor(2);
1281  g_photonpqcd05->SetLineStyle(2);
1282  g_photonpqcd2->SetLineWidth(2);
1283  g_photonpqcd2->SetLineColor(2);
1284  g_photonpqcd2->SetLineStyle(2);
1285 
1286  m_doubleratio->Add(g_photonpqcd,"c");
1287  m_doubleratio->Add(g_photonpqcd05,"c");
1288  m_doubleratio->Add(g_photonpqcd2,"c");
1289 
1290  //appropriate fit to photon pqcd result
1291  TF1 *fitGamma2=new TF1("fitGamma2","1.+[0]*TMath::Power(x,[1])",2.,15.);
1292 assert(fitGamma2);
1293  g_photonpqcd->Fit(fitGamma2,"R0");
1294 
1295  m_doubleratio->Draw("a");
1296 
1297  m_doubleratio->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1298  m_doubleratio->GetYaxis()->SetTitle("1 + #gamma_{dir}/#gamma_{incl}");
1299  m_doubleratio->GetXaxis()->SetRangeUser(2.,16.);
1300 
1301 
1302  TLegend *leg5=new TLegend(.15,.6,.6,.8);
1303 assert(leg5);
1304  leg5->AddEntry(g_doubleRatioHT1,"hightower-1","p");
1305  leg5->AddEntry(g_doubleRatioHT2,"hightower-2","p");
1306  leg5->AddEntry(g_photonpqcd,"NLO (CTEQ6+KKP) #mu=p_{T}","l");
1307  leg5->AddEntry(g_photonpqcd05,"#mu=2p_{T}, #mu=p_{T}/2","l");
1308  leg5->SetFillColor(0);
1309  leg5->Draw("same");
1310 
1311  c_doubleratio->cd(0);
1312  c_doubleratio->SaveAs(settings.output_gammadoubleratio_file.Data());
1313 
1314  {
1315  // calculating antineutron contamination C0
1316  TH1F *h_nbarContMB = new TH1F(*h_nbarEffMB);
1317  TH1F *h_nbarContHT1 = new TH1F(*h_nbarEffHT1);
1318  TH1F *h_nbarContHT2 = new TH1F(*h_nbarEffHT2);
1319 assert(h_nbarContMB);
1320 assert(h_nbarContHT1);
1321 assert(h_nbarContHT2);
1322  h_nbarContMB->Multiply(f_nbar);
1323  h_nbarContHT1->Multiply(f_nbar);
1324  h_nbarContHT2->Multiply(f_nbar);
1325 
1326  h_nbarContMB->Divide(gammaYieldMB);
1327  h_nbarContHT1->Divide(gammaYieldHT1);
1328  h_nbarContHT2->Divide(gammaYieldHT2);
1329 
1330  TH1F *effGammaMBnoerr = new TH1F(*effGammaMB);
1331  TH1F *effGammaHT1noerr = new TH1F(*effGammaHT1);
1332  TH1F *effGammaHT2noerr = new TH1F(*effGammaHT2);
1333 assert(effGammaMBnoerr);
1334 assert(effGammaHT1noerr);
1335 assert(effGammaHT2noerr);
1336  for (Int_t i = 1;i <= effGammaMBnoerr->GetXaxis()->GetNbins();i++) effGammaMBnoerr->SetBinError(i, 0);
1337  for (Int_t i = 1;i <= effGammaHT1noerr->GetXaxis()->GetNbins();i++) effGammaHT1noerr->SetBinError(i, 0);
1338  for (Int_t i = 1;i <= effGammaHT2noerr->GetXaxis()->GetNbins();i++) effGammaHT2noerr->SetBinError(i, 0);
1339  h_nbarContMB->Divide(effGammaMBnoerr);
1340  h_nbarContHT1->Divide(effGammaHT1noerr);
1341  h_nbarContHT2->Divide(effGammaHT2noerr);
1342 
1343  TCanvas *c_nbar_cont = new TCanvas("c_nbar_cont", "Antineutron contamination");
1344 assert(c_nbar_cont);
1345  TH1F *h_nbar_cont = new TH1F("h_nbar_cont", "Antineutron contamination C_{0};p_{T} [GeV/c];C_{0}", 1000, 0, 10);
1346 assert(h_nbar_cont);
1347  h_nbar_cont->Draw();
1348  h_nbar_cont->GetYaxis()->SetRangeUser(0, 1.5);
1349 
1350  h_nbarContMB->SetMarkerStyle(kOpenCircle);
1351  h_nbarContMB->SetMarkerSize(1.0);
1352  h_nbarContMB->SetMarkerColor(kBlack);
1353  h_nbarContMB->SetLineStyle(kSolid);
1354  h_nbarContMB->SetLineWidth(1);
1355  h_nbarContMB->SetLineColor(kBlack);
1356 
1357  h_nbarContHT1->SetMarkerStyle(kOpenSquare);
1358  h_nbarContHT1->SetMarkerSize(1.0);
1359  h_nbarContHT1->SetMarkerColor(kBlue);
1360  h_nbarContHT1->SetLineStyle(kSolid);
1361  h_nbarContHT1->SetLineWidth(1);
1362  h_nbarContHT1->SetLineColor(kBlue);
1363 
1364  h_nbarContHT2->SetMarkerStyle(kOpenTriangleUp);
1365  h_nbarContHT2->SetMarkerSize(1.0);
1366  h_nbarContHT2->SetMarkerColor(kRed);
1367  h_nbarContHT2->SetLineStyle(kSolid);
1368  h_nbarContHT2->SetLineWidth(1);
1369  h_nbarContHT2->SetLineColor(kRed);
1370 
1371  h_nbarContMB->Draw("SAME");
1372  h_nbarContHT1->Draw("SAME");
1373  h_nbarContHT2->Draw("SAME");
1374 
1375  // extreme case - at maximum C0 reaches unity
1376  Float_t maxCont = h_nbarContMB->GetMaximum();
1377  TH1F *h_nbarContHT2_extreme = new TH1F(*h_nbarContHT2);
1378 assert(h_nbarContHT2_extreme);
1379  h_nbarContHT2_extreme->Scale(1.0/maxCont);
1380  h_nbarContHT2_extreme->SetLineStyle(kDashed);
1381  h_nbarContHT2_extreme->SetLineWidth(1);
1382  h_nbarContHT2_extreme->SetLineColor(kRed);
1383  h_nbarContHT2_extreme->Draw("SAME HIST C");
1384 
1385  // another case - at 1 < pT < 4 GeV/c direct photons vanish,
1386  // therefore an excess of MinBias double ratio over unity in that region
1387  // is attributed to antineutrons, so that C0 = 1 - (1/R)
1388  TH1F *h_corr = new TH1F(*gammaYieldMBdoubleratio);
1389 assert(h_corr);
1390  for (Int_t i = 1;i < h_corr->GetXaxis()->GetNbins();i++) {
1391  h_corr->SetBinContent(i, 1.0);
1392  h_corr->SetBinError(i, 0.0);
1393  }
1394  h_corr->Divide(gammaYieldMBdoubleratio);
1395  //h_corr->Scale(-1.0);
1396  for (Int_t i = 1;i < h_corr->GetXaxis()->GetNbins();i++) {
1397  h_corr->SetBinContent(i, 1.0 - h_corr->GetBinContent(i));
1398  }
1399  h_corr->Divide(h_nbarContMB);
1400  h_corr->GetXaxis()->SetRangeUser(1.0, 4.0);
1401  TF1 *f_corr = new TF1("f_corr", "[0]");
1402  f_corr->SetRange(1.0, 4.0);
1403  h_corr->Fit(f_corr, "RQN");
1404  Float_t corrCoef = f_corr->GetParameter(0);
1405  TH1F *h_nbarContHT2_corr = new TH1F(*h_nbarContHT2);
1406 assert(h_nbarContHT2_corr);
1407  h_nbarContHT2_corr->Scale(corrCoef);
1408  h_nbarContHT2_corr->SetLineStyle(kSolid);
1409  h_nbarContHT2_corr->SetLineWidth(1);
1410  h_nbarContHT2_corr->SetLineColor(17);
1411  h_nbarContHT2_corr->SetFillColor(17);
1412  h_nbarContHT2_corr->SetFillStyle(1001);
1413  h_nbarContHT2_corr->Draw("SAME HIST AC");
1414 
1415  TLegend *leg = new TLegend(0.5, 0.5, 0.8, 0.8);
1416 assert(leg);
1417  leg->AddEntry(h_nbarContMB, "MinBias", "P");
1418  leg->AddEntry(h_nbarContHT1, "HighTower-1", "P");
1419  leg->AddEntry(h_nbarContHT2, "HighTower-2", "P");
1420  leg->AddEntry(h_nbarContHT2_extreme, "extreme case", "L");
1421  leg->AddEntry(h_nbarContHT2_corr, "corrected", "F");
1422  leg->Draw();
1423 
1424  h_nbar_cont->Draw("SAME");
1425  c_nbar_cont->SaveAs(settings.output_nbarcont_file.Data());
1426  }
1427 
1428  TCanvas *c_dirphoton=new TCanvas("c_dirphoton","c_dirphoton",400,300);
1429 assert(c_dirphoton);
1430  gStyle->SetOptStat(0);
1431  c_dirphoton->cd(1);
1432 
1433  TH1F *dirphotonYieldMB=new TH1F(*gammaYieldMBdoubleratio);
1434  TH1F *dirphotonYieldHT1=new TH1F(*gammaYieldHT1doubleratio);
1435  TH1F *dirphotonYieldHT2=new TH1F(*gammaYieldHT2doubleratio);
1436 assert(dirphotonYieldMB);
1437 assert(dirphotonYieldHT1);
1438 assert(dirphotonYieldHT2);
1439 
1440  TH1F *dirphotonYieldMBnoErr=new TH1F(*dirphotonYieldMB);
1441 assert(dirphotonYieldMBnoErr);
1442  for(int i=1;i<=dirphotonYieldMBnoErr->GetNbinsX();i++){
1443  dirphotonYieldMBnoErr->SetBinError(i,0.);
1444  }
1445  TH1F *dirphotonYieldHT1noErr=new TH1F(*dirphotonYieldHT1);
1446 assert(dirphotonYieldHT1noErr);
1447  for(int i=1;i<=dirphotonYieldHT1noErr->GetNbinsX();i++){
1448  dirphotonYieldHT1noErr->SetBinError(i,0.);
1449  }
1450  TH1F *dirphotonYieldHT2noErr=new TH1F(*dirphotonYieldHT2);
1451 assert(dirphotonYieldHT2noErr);
1452  for(int i=1;i<=dirphotonYieldHT1noErr->GetNbinsX();i++){
1453  dirphotonYieldHT2noErr->SetBinError(i,0.);
1454  }
1455 
1456  TF1 *f_unity=new TF1("f_unity","1.",0.,15.);
1457 assert(f_unity);
1458  dirphotonYieldMB->Add(f_unity,-1.);
1459  dirphotonYieldHT1->Add(f_unity,-1.);
1460  dirphotonYieldHT2->Add(f_unity,-1.);
1461 
1462 
1463  dirphotonYieldMB->Divide(dirphotonYieldMBnoErr);
1464  dirphotonYieldHT1->Divide(dirphotonYieldHT1noErr);
1465  dirphotonYieldHT2->Divide(dirphotonYieldHT2noErr);
1466  dirphotonYieldMB->Multiply(gammaYieldMBratioNoErr);
1467  dirphotonYieldHT1->Multiply(gammaYieldHT1ratioNoErr);
1468  dirphotonYieldHT2->Multiply(gammaYieldHT2ratioNoErr);
1469  dirphotonYieldMB->Multiply(pionXsMBnoErr);
1470  dirphotonYieldHT1->Multiply(pionXsHT1noErr);
1471  dirphotonYieldHT2->Multiply(pionXsHT2noErr);
1472 
1473 
1474  TGraphErrors *g_dirphotonMB=new TGraphErrors(dirphotonYieldMB);
1475 assert(g_dirphotonMB);
1476  g_dirphotonMB->SetName("g_dirphotonMB");
1477  g_dirphotonMB->SetMarkerStyle(8);
1478  TGraphErrors *g_dirphotonHT1=new TGraphErrors(dirphotonYieldHT1);
1479 assert(g_dirphotonHT1);
1480  g_dirphotonHT1->SetName("g_dirphotonHT1");
1481  g_dirphotonHT1->SetMarkerStyle(8);
1482  TGraphErrors *g_dirphotonHT2=new TGraphErrors(dirphotonYieldHT2);
1483 assert(g_dirphotonHT2);
1484  g_dirphotonHT2->SetName("g_dirphotonHT2");
1485  g_dirphotonHT2->SetMarkerStyle(8);
1486 
1487 
1488  removeThesePoints(g_dirphotonMB,1);
1489  removeThesePoints(g_dirphotonHT1,2);
1490  removeThesePoints(g_dirphotonHT2,3);
1491 
1492  gPad->SetLogy();
1493 
1494  TMultiGraph *m_dirphoton=new TMultiGraph();
1495 assert(m_dirphoton);
1496  m_dirphoton->SetName("m_dirphoton");
1497  m_dirphoton->SetMinimum(1.0e-11);
1498  m_dirphoton->SetMaximum(0.1);
1499 
1500  m_dirphoton->Add(g_dirgamma,"c");
1501  m_dirphoton->Add(g_dirgamma05,"c");
1502  m_dirphoton->Add(g_dirgamma2,"c");
1503 
1504  cout<<"direct photons:"<<endl;
1505  g_dirphotonHT1->Print();
1506  cout<<endl;
1507  g_dirphotonHT2->Print();
1508  cout<<endl;
1509 
1510  m_dirphoton->Add(g_dirphotonHT1,"p");
1511  m_dirphoton->Add(g_dirphotonHT2,"p");
1512 
1513  m_dirphoton->Draw("a");
1514 
1515  m_dirphoton->GetXaxis()->SetTitle("p_{T} (GeV/c)");
1516  m_dirphoton->GetYaxis()->SetTitle("1 - R^{-1}");
1517  m_dirphoton->GetXaxis()->SetRangeUser(2.,16.);
1518 
1519 
1520  TLegend *leght=new TLegend(.15,.6,.6,.8);
1521 assert(leght);
1522  leght->AddEntry(g_dirphotonHT1,"hightower-1","p");
1523  leght->AddEntry(g_dirphotonHT2,"hightower-2","p");
1524  leght->SetFillColor(0);
1525  leght->Draw("same");
1526 
1527  c_dirphoton->cd(0);
1528  c_dirphoton->SaveAs(settings.output_gammadirphoton_file.Data());
1529 }
1530 
1531 
Generated by   doxygen 1.8.5