StRoot: StEmcPool/StPhotonAnalysis/macros/getSpectrumPP05.C Source File

00001 
00002 #include "systematics/systematics_pp.h"
00003 
00004 double pqcdlowpt(double *x,double *par)
00005 {
00006   double ret=par[0]*pow(1.+x[0],par[1]);
00007   return ret;
00008 }
00009 double pqcdhighpt(double *x,double *par)
00010 {
00011   double ret=par[0]*pow(1.+x[0],par[1]);
00012   return ret;
00013 }
00014 double sumpqcd(double *x,double *par)
00015 {
00016   //transition (SW=1 is old)
00017   //params: 2x low, 2x hi, 2x trans
00018   double SW=1. - 1./(1.+exp(fabs(par[4])*(x[0]-par[5])));
00019   double ret=SW*pqcdhighpt(x,&par[2]);
00020   ret+=(1.-SW)*pqcdlowpt(x,par);
00021   //ret=pqcdlowpt(x,par);
00022   return ret;
00023 }
00024 void printPoints(TGraphErrors *g)
00025 {
00026   for(Int_t i=0;i<g->GetN();i++){
00027     double x,y,ey;
00028     g->GetPoint(i,x,y);
00029     ey=g->GetErrorY(i);
00030     cout<<x<<"   "<<y<<"   "<<ey<<endl;
00031   }
00032 }
00033 void divideGraphWithFunction(TGraphErrors *graph,TF1 *func)
00034 {
00035   for(int i=0;i<graph->GetN();i++){
00036     double x=0.;
00037     double y=0.;
00038     graph->GetPoint(i,x,y);
00039     double ey=graph->GetErrorY(i);
00040     graph->SetPoint(i,x,y/func->Eval(x));
00041     graph->SetPointError(i,0.,ey/func->Eval(x));
00042   }
00043 }
00044 void divideGraphWithGraph(TGraphErrors *graph,TGraph *denom)
00045 {
00046   for(int i=0;i<graph->GetN();i++){
00047     double x=0.;
00048     double y=0.;
00049     graph->GetPoint(i,x,y);
00050     double ey=graph->GetErrorY(i);
00051     graph->SetPoint(i,x,y/denom->Eval(x));
00052     graph->SetPointError(i,0.,ey/denom->Eval(x));
00053   }
00054 }
00055 void getRatio(TH1F* h_ratio,TH1F* h_eff,TH1F* h_effSingle,TF1 *f_piondecay)
00056 {
00057   //copy, don't want to mess with the pointers:
00058   TH1F *h_eff_copy=new TH1F(*h_eff);
00059   for(int i=1;i<=h_eff->GetNbinsX();i++){
00060     h_eff->SetBinError(i,0.);
00061   }
00062   TH1F *h_effSingle_copy=new TH1F(*h_effSingle);
00063   for(int i=1;i<=h_effSingle->GetNbinsX();i++){
00064     h_effSingle_copy->SetBinError(i,0.);
00065   }
00066   //correct for single photon fraction:
00067   h_ratio->Add(f_piondecay,-1.);
00068   for(int i=1;i<=h_ratio->GetNbinsX();i++){
00069     //cout<<i<<endl;
00070     //if(h_ratio->GetBinContent(i)<0.) cout<<"negative"<<endl;
00071   }
00072   h_eff_copy->Divide(h_effSingle_copy);
00073   h_ratio->Multiply(h_eff_copy);
00074   h_ratio->Add(f_piondecay);
00075 }
00076 void removeThesePoints(TGraphErrors *g,int trig)
00077 {
00078   if(trig==1){
00079     g->RemovePoint(0);
00080     g->RemovePoint(0);
00081     //additional for pp:
00082     g->RemovePoint(g->GetN()-1);
00083   }
00084   if(trig==2){
00085     g->RemovePoint(0);
00086     g->RemovePoint(0);
00087     g->RemovePoint(0);
00088     g->RemovePoint(0);
00089     g->RemovePoint(g->GetN()-1);
00090     g->RemovePoint(g->GetN()-1);
00091   }
00092   if(trig==3){
00093     g->RemovePoint(0);
00094     g->RemovePoint(0);
00095     g->RemovePoint(0);
00096     g->RemovePoint(0);
00097     g->RemovePoint(0);
00098   }
00099 }
00100 void getSpectrumPP05(const char *pid)
00101 {
00102   gStyle->SetOptStat(0);
00103   gStyle->SetOptFit(0);
00104 
00105   gStyle->SetErrorX(0);
00106 
00107 
00108 
00109   bool subtractNEUTRONS=kTRUE;
00110 
00111   //neutron data:
00112   //hadron data for Levy function, nucl-ex/0601033:
00113   //-->  d^2N/(2*pi*pT*N*dpT*dy) = B/((1+((mT - m0)/nT))^n)
00114   // {p-dAu; pbar-dAu; p-pp; pbar-pp} and m0 = m_neutron = 1.0 GeV.
00115   double B[]={0.3,0.23,0.072,0.061};//->[0]
00116   //double eB[]={0.01,0.01,0.005,0.005};
00117   double T[]={0.205,0.215,0.179,0.173};//->[1]
00118   //double eT[]={0.004,0.005,0.006,0.006};
00119   double n[]={11.00,12.55,10.87,10.49};//->[2]
00120   //double en[]={0.29,0.41,0.43,0.40};
00121   double BR=35.8/63.9;//->p + pi- (63.9%) : ->n + pi0 (35.8%)
00122   double c_feeddown=0.8+0.2*BR;
00123   double CPVloss=0.98;//for minbias
00124   TF1 *f_nbar=new TF1("f_nbar","[3]*[4]*[0]/pow((1.+(sqrt(x*x+1.) - 1.)/([1]*[2])),[2])",0.,15.);
00125   f_nbar->SetParameters(B[3],T[3],n[3],c_feeddown,CPVloss);
00126 
00127   //get direct gammas pQCD:
00128   ifstream pQCDphotons("./datapoints/pQCD_Werner/rhic_cteq6_gamma_inv_sc1.dat");
00129   float ppx[100];
00130   float ppy[100];
00131   Int_t iii=0;
00132   cout<<"pqcd photons:"<<endl;
00133   while(iii<28){
00134     if(!pQCDphotons.good()) break;
00135     float dummy=0.;
00136     pQCDphotons>>ppx[iii]>>dummy>>dummy>>ppy[iii];
00137     ppy[iii]*=1.e-09;//convert to mb
00138     iii++;
00139   }
00140   TGraph *g_dirgamma=new TGraph(iii,ppx,ppy);
00141   //get direct gammas pQCD: scale 0.5*pT
00142   ifstream pQCDphotons05("./datapoints/pQCD_Werner/rhic_cteq6_gamma_inv_sc05.dat");
00143   float ppx05[100];
00144   float ppy05[100];
00145   Int_t iii05=0;
00146   while(iii05<28){
00147     if(!pQCDphotons05.good()) break;
00148     float dummy=0.;
00149     pQCDphotons05>>ppx05[iii05]>>dummy>>dummy>>ppy05[iii05];
00150     ppy05[iii05]*=1.e-09;//convert to mb
00151     iii05++;
00152   }
00153   TGraph *g_dirgamma05=new TGraph(iii05,ppx05,ppy05);
00154 
00155   //get direct gammas pQCD: scale 2*pT
00156   ifstream pQCDphotons2("./datapoints/pQCD_Werner/rhic_cteq6_gamma_inv_sc2.dat");
00157   float ppx2[100];
00158   float ppy2[100];
00159   Int_t iii2=0;
00160   while(iii2<28){
00161     if(!pQCDphotons2.good()) break;
00162     float dummy=0.;
00163     pQCDphotons2>>ppx2[iii2]>>dummy>>dummy>>ppy2[iii2];
00164     ppy2[iii2]*=1.e-09;//convert to mb
00165     iii2++;
00166   }
00167   TGraph *g_dirgamma2=new TGraph(iii2,ppx2,ppy2);
00168   //get phenix pions in pp
00169   ifstream phenix("./datapoints/phenix_xsec_pp.dat");
00170   float phex[100];
00171   float phey[100];
00172   Int_t iphe=0;
00173   while(iphe<16){
00174     if(!phenix.good()) break;
00175     phenix>>phex[iphe]>>phey[iphe];
00176     //is in mb already!!
00177     iphe++;
00178   }
00179   TGraphErrors *phenix_pp=new TGraphErrors(iphe,phex,phey);
00180   phenix_pp->SetMarkerStyle(24);
00181   phenix_pp->SetLineWidth(6);
00182   phenix_pp->SetName("phenix_pp");
00183 
00184   //frank's spin2006 pions:
00185   ifstream frank("./datapoints/frank_pp05_new.dat");
00186   float frx[100];
00187   float fry[100];
00188   float frex[100];
00189   float frey[100];
00190   Int_t ifr=0;
00191   while(ifr<12){
00192     if(!frank.good()) break;
00193     frank>>frx[ifr]>>fry[ifr]>>frey[ifr];
00194     frex[ifr]=0.;
00195     ifr++;
00196   }
00197   TGraphErrors *frank_pp=new TGraphErrors(ifr,frx,fry,frex,frey);
00198   frank_pp->SetMarkerStyle(25);
00199   frank_pp->SetMarkerColor(4);
00200   frank_pp->SetName("frank_pp");
00201 
00202   //sasha's pions:
00203   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 };
00204   float xe_pp2005_sasha[]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
00205   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};
00206   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};
00207 
00208   TGraphErrors *sasha_pp05=new TGraphErrors(13,x_pp2005_sasha,y_pp2005_sasha,xe_pp2005_sasha,ye_pp2005_sasha);
00209   sasha_pp05->SetMarkerStyle(8);
00210   sasha_pp05->SetMarkerColor(TColor::GetColor(8,80,8));
00211   sasha_pp05->SetName("sasha_pp05");
00212 
00213   //get pions KKP  scale pT
00214   ifstream pQCDpions("./datapoints/pQCD_Werner/klaus_pi0inv_200_kkp_1.dat");
00215   float pionx[100];
00216   float piony[100];
00217   Int_t ipion=0;
00218   while(ipion<28){
00219     if(!pQCDpions.good()) break;
00220     pQCDpions>>pionx[ipion]>>piony[ipion];
00221     piony[ipion]*=1.e-09;//convert to mb
00222     ipion++;
00223   }
00224   TGraphErrors *kkp=new TGraphErrors(ipion,pionx,piony);
00225   kkp->SetLineColor(54);
00226   kkp->SetName("kkp");
00227 
00228   //get pions KKP  scale 0.5*pT
00229   ifstream pQCDpions05("./datapoints/pQCD_Werner/klaus_pi0inv_200_kkp_05.dat");
00230   float pionx05[100];
00231   float piony05[100];
00232   Int_t ipion05=0;
00233   while(ipion05<28){
00234     if(!pQCDpions05.good()) break;
00235     pQCDpions05>>pionx05[ipion05]>>piony05[ipion05];
00236     piony05[ipion05]*=1.e-09;//convert to mb
00237     ipion05++;
00238   }
00239   TGraphErrors *kkp05=new TGraphErrors(ipion05,pionx05,piony05);
00240   kkp05->SetLineStyle(2);
00241   kkp05->SetLineColor(54);
00242   kkp05->SetName("kkp05");
00243 
00244   //get pions KKP  scale 2*pT
00245   ifstream pQCDpions2("./datapoints/pQCD_Werner/klaus_pi0inv_200_kkp_2.dat");
00246   float pionx2[100];
00247   float piony2[100];
00248   Int_t ipion2=0;
00249   while(ipion2<28){
00250     if(!pQCDpions2.good()) break;
00251     pQCDpions2>>pionx2[ipion2]>>piony2[ipion2];
00252     piony2[ipion2]*=1.e-09;//convert to mb
00253     ipion2++;
00254   }
00255   TGraphErrors *kkp2=new TGraphErrors(ipion2,pionx2,piony2);
00256   kkp2->SetLineStyle(2);
00257   kkp2->SetLineColor(54);
00258   kkp2->SetName("kkp2");
00259 
00260   TFile *f_decaybg=new TFile("~/MyDecay/gammaDecayPPSum.root","OPEN");
00261   TH1F *h_decaybg=(TH1F*)f_decaybg->Get("gamma");
00262   TH1F *h_decaypion=(TH1F*)f_decaybg->Get("gamma_pion");
00263   TF1 *fit_decay=new TF1("fit_decay","[0]/pow(x,[1])+[2]",.3,15.);
00264   TF1 *fit_piondecay=new TF1(*fit_decay);
00265   fit_decay->SetParameters(1.,1.,.5);
00266   fit_piondecay->SetParameters(1.,1.,.5);
00267   h_decaybg->Fit(fit_decay,"R0Q");
00268   h_decaypion->Fit(fit_piondecay,"R0Q");
00269 
00270   //take ratio gamma_direct/pion and divide by gamma/pi and then +1:
00271   for(Int_t i=0;i<iii;i++){
00272     ppy[i]=ppy[i]/piony[i];
00273     ppy[i]=ppy[i]/fit_decay->Eval(ppx[i]);
00274     ppy[i]+=1.;
00275     ppy05[i]=ppy05[i]/piony05[i];
00276     ppy05[i]=ppy05[i]/fit_decay->Eval(ppx05[i]);
00277     ppy05[i]+=1.;
00278     ppy2[i]=ppy2[i]/piony2[i];
00279     ppy2[i]=ppy2[i]/fit_decay->Eval(ppx2[i]);
00280     ppy2[i]+=1.;
00281   }
00282   TGraphErrors *g_photonpqcd=new TGraphErrors(iii,ppx,ppy);
00283   g_photonpqcd->SetName("g_photonpqcd");
00284   TGraphErrors *g_photonpqcd05=new TGraphErrors(iii05,ppx05,ppy05);
00285   g_photonpqcd05->SetName("g_photonpqcd05");
00286   TGraphErrors *g_photonpqcd2=new TGraphErrors(iii2,ppx2,ppy2);
00287   g_photonpqcd2->SetName("g_photonpqcd2");
00288 
00289   //set outputfiles
00290   TString dir("/star/u/russcher/gamma/analysis/output/pp05/");
00291   dir.Append(pid);
00292   dir.Append("/");
00293   dir.Append(pid);
00294   TString psout=dir;
00295   TString psout_eta=dir;
00296   TString psout2=dir;
00297   TString psout3=dir;
00298   TString eFile=dir;
00299   TString eFileGamma=dir;
00300   TString pi0File=dir;
00301   TString nbarFile=dir;
00302   TString nbarOut=dir;
00303   
00304   eFile.Append("pion_eff.root");
00305   eFileGamma.Append("gamma_eff.root");
00306   pi0File.Append("pi0_pp05.root");
00307   nbarFile.Append("antineutron_eff.root");
00308   nbarOut.Append("nbar_contam.ps");
00309   psout2.Append("correctedspec.ps");
00310   psout_eta.Append("/dev/null");
00311   psout.Append("invmassplots.ps");
00312 
00313   //load *.so
00314   gSystem->Load("$HOME/MyEvent/MyEvent");
00315   gSystem->Load("$HOME/gamma/analysis/lib/AnaCuts");
00316   gSystem->Load("$HOME/gamma/analysis/lib/EventMixer");
00317   gSystem->Load("$HOME/gamma/analysis/lib/Pi0Analysis");
00318 
00319 
00320   AnaCuts *cuts=new AnaCuts("pp05");
00321 
00322   //get inv mass hists
00323   TFile f(pi0File.Data(),"OPEN");
00324   TH2F *h_mb=new TH2F(*h_minvMB);
00325   TH2F *h_ht1=new TH2F(*h_minvHT1);
00326   TH2F *h_ht2=new TH2F(*h_minvHT2);
00327   TH1F *h_ev=new TH1F(*h_events);
00328 
00329   //get anti-neutron
00330   TFile *file_nbar=new TFile(nbarFile,"OPEN");
00331   TH1F *h_nbarEffMB=new TH1F(*h_effMB);
00332   h_nbarEffMB->Sumw2();
00333   TH1F *h_nbarEffHT1=new TH1F(*h_effHT1);
00334   h_nbarEffHT1->Sumw2();
00335   TH1F *h_nbarEffHT2=new TH1F(*h_effHT2);
00336   h_nbarEffHT2->Sumw2();
00337 
00338   //get prescales
00339   int trigger=0;
00340   Int_t numberOfMB=0;
00341   Int_t numberOfHT1=0;
00342   Int_t numberOfHT2=0;
00343   for(Int_t i=1;i<=h_ev->GetNbinsX();i++)
00344     {
00345       trigger=(Int_t)h_ev->GetBinCenter(i);
00346       if(trigger&1) numberOfMB+=(Int_t)h_ev->GetBinContent(i);
00347       if(trigger&2) numberOfHT1+=(Int_t)h_ev->GetBinContent(i);
00348       if(trigger&4) numberOfHT2+=(Int_t)h_ev->GetBinContent(i);
00349     }
00350 
00351   cout<<"number of mb: "<<numberOfMB<<endl;
00352   //numberOfMB/=0.63;
00353   //cout<<"nmb after 62% vertex eff.: "<<numberOfMB<<endl;
00354 
00355   float psMB=32742;
00356   float psHT1=3.89;
00357   float psHT2=1.;
00358 
00359   //mb events for hightower:
00360   float nMBwithHT=328871;
00361 
00362   //get efficiencies+acceptance
00363   TFile g(eFile.Data(),"OPEN");
00364   TH1F *h_emb=new TH1F(*h_effMB);
00365   TH1F *h_eht1=new TH1F(*h_effHT1);
00366   TH1F *h_eht2=new TH1F(*h_effHT2);
00367 
00368   //bin corrections
00369   TFile binf("~/BinWidth/bincorrectionsPP.root","OPEN"); 
00370   TH1F *h_binmb=new TH1F(*h4mb);
00371   TH1F *h_binht1=new TH1F(*h4ht1);
00372   TH1F *h_binht2=new TH1F(*h4ht2);
00373   h_binmb->Sumw2();
00374   h_binht1->Sumw2();
00375   h_binht2->Sumw2();
00376   for(Int_t i=1;i<=h_binmb->GetNbinsX();i++) h_binmb->SetBinError(i,0);
00377   for(Int_t i=1;i<=h_binht1->GetNbinsX();i++) h_binht1->SetBinError(i,0);
00378   for(Int_t i=1;i<=h_binht2->GetNbinsX();i++) h_binht2->SetBinError(i,0);
00379 
00380   //corrections, all multiplicative, in case of
00381   //pions:
00382   TF1 *pion_cpv_corrMB=new TF1("pion_cpv_corrMB","1./(1.-0.01*(0.3+0.0*x))",0.,15.);
00383   TF1 *pion_cpv_corrHT1=new TF1("pion_cpv_corrHT1","1./(1.-0.01*(-0.1+0.16*x))",0.,15.);
00384   TF1 *pion_cpv_corrHT2=new TF1("pion_cpv_corrHT2","1./(1.-0.01*(-0.2+0.18*x))",0.,15.);
00385   //photons:
00386   TF1 *gamma_cpv_corrMB=new TF1("gamma_cpv_corrMB","1./(1.-0.01*(2.8+0.0*x))",0.,15.);
00387   TF1 *gamma_cpv_corrHT1=new TF1("gamma_cpv_corrHT1","1./(1.-0.01*(0.2+1.1*x))",0.,15.);
00388   TF1 *gamma_cpv_corrHT2=new TF1("gamma_cpv_corrHT2","1./(1.-0.01*(0.4+1.1*x))",0.,15.);
00389   TF1 *gamma_cont_corrMB=new TF1("gamma_cont_corrMB","0.985",0.,15.);
00390   TF1 *gamma_cont_corrHT1=new TF1("gamma_cont_corrHT1","0.98",0.,15.);
00391   TF1 *gamma_cont_corrHT2=new TF1("gamma_cont_corrHT2","0.96",0.,15.);
00392 
00393   //missing material
00394   //pions:
00395   TF1 *pion_conv_corrMB=new TF1("pion_conv_corrMB","1.1",0.,15.);
00396   TF1 *pion_conv_corrHT1=new TF1("pion_conv_corrHT1","1.1",0.,15.);
00397   TF1 *pion_conv_corrHT2=new TF1("pion_conv_corrHT2","1.1",0.,15.);
00398   //photons:
00399   TF1 *gamma_conv_corrMB=new TF1("gamma_conv_corrMB","1.05",0.,15.);
00400   TF1 *gamma_conv_corrHT1=new TF1("gamma_conv_corrHT1","1.05",0.,15.);
00401   TF1 *gamma_conv_corrHT2=new TF1("gamma_conv_corrHT2","1.05",0.,15.);
00402 
00403 
00404   //get yield
00405   Pi0Analysis *pi0=new Pi0Analysis(psout.Data(),psout_eta.Data(),"pp05");
00406   pi0->init("/dev/null");
00407   TH1F *pionYieldMB=new TH1F(*pi0->getYield(h_mb,"mb"));
00408   TH1F *pionYieldHT1=new TH1F(*pi0->getYield(h_ht1,"ht1"));
00409   TH1F *pionYieldHT2=new TH1F(*pi0->getYield(h_ht2,"ht2"));
00410   pi0->storeCanvases((dir+"canvases.root").Data());
00411 
00412 
00413   cout<<"***************************************"<<endl;
00414   cout<<"got yield, dividing by rapidity bite!!!"<<endl;
00415   float dy_gamma=cuts->rapidityMaxCUT - cuts->rapidityMinCUT;
00416   float dy_pion=cuts->rapPionMaxCUT - cuts->rapPionMinCUT;
00417   cout<<"***************************************"<<endl;
00418   cout<<endl;
00419   cout<<"    pion bite is "<<dy_pion<<endl;
00420   cout<<"    gamma bite is "<<dy_gamma<<endl;
00421   cout<<endl;
00422   cout<<"***************************************"<<endl;
00423 
00424   pionYieldMB->Scale(1./dy_pion);
00425   pionYieldHT1->Scale(1./dy_pion);
00426   pionYieldHT2->Scale(1./dy_pion);
00427 
00428   //set yield to zero
00429   /*
00430   for(Int_t i=0;i<pionYieldHT2->GetNbinsX();i++)
00431     {
00432       if(i<1)
00433         {
00434           pionYieldMB->SetBinContent(i+1,0);
00435           pionYieldMB->SetBinError(i+1,0);
00436         }
00437       if(i<4)
00438         {
00439           pionYieldHT1->SetBinContent(i+1,0);
00440           pionYieldHT1->SetBinError(i+1,0);
00441         }
00442       if(i>12)
00443         {
00444           pionYieldHT1->SetBinContent(i+1,0);
00445           pionYieldHT1->SetBinError(i+1,0);
00446         }
00447       if(i<6) 
00448         {
00449           pionYieldHT2->SetBinContent(i+1,0);
00450           pionYieldHT2->SetBinError(i+1,0);
00451         }
00452     }
00453   */
00454 
00455   //set colors:
00456   pionYieldMB->SetMarkerStyle(8);
00457   pionYieldMB->SetMarkerSize(1.0);
00458   pionYieldHT1->SetMarkerStyle(8);
00459   pionYieldHT1->SetMarkerSize(1.0); 
00460   pionYieldHT1->SetMarkerColor(4);
00461   pionYieldHT2->SetMarkerStyle(8);
00462   pionYieldHT2->SetMarkerSize(1.0);
00463   pionYieldHT2->SetMarkerColor(2);
00464   
00465   TF1 *scale=new TF1("scale","x",0.,15.);
00466 
00467   pionYieldMB->SetNameTitle("pionYieldMB","corrected yield MB");
00468   pionYieldMB->Divide(h_emb);
00469   pionYieldMB->Scale(psMB/(psMB*numberOfMB*2.*TMath::Pi()));
00470   pionYieldMB->Divide(scale);
00471   pionYieldMB->Multiply(h_binmb);
00472   pionYieldMB->Multiply(pion_cpv_corrMB);
00473   pionYieldMB->Multiply(pion_conv_corrMB);
00474 
00475   pionYieldHT1->SetNameTitle("pionYieldHT1","corrected yield HT1");
00476   pionYieldHT1->Divide(h_eht1);
00477   pionYieldHT1->Scale(psHT1/(psMB*nMBwithHT*2.*TMath::Pi()));
00478   pionYieldHT1->Divide(scale);
00479   pionYieldHT1->Multiply(h_binht1);
00480   pionYieldHT1->Multiply(pion_cpv_corrHT1);
00481   pionYieldHT1->Multiply(pion_conv_corrHT1);
00482 
00483   pionYieldHT2->SetNameTitle("pionYieldHT2","corrected yield HT2");
00484   pionYieldHT2->Divide(h_eht2);
00485   pionYieldHT2->Scale(psHT2/(psMB*nMBwithHT*2.*TMath::Pi()));
00486   pionYieldHT2->Divide(scale);
00487   pionYieldHT2->Multiply(h_binht2);
00488   pionYieldHT2->Multiply(pion_cpv_corrHT2);
00489   pionYieldHT2->Multiply(pion_conv_corrHT2);
00490 
00491   //create pion yield for double ratio:
00492   TH1F *pionYieldMBratio=new TH1F(*pionYieldMB);
00493   TH1F *pionYieldHT1ratio=new TH1F(*pionYieldHT1);
00494   TH1F *pionYieldHT2ratio=new TH1F(*pionYieldHT2);
00495   
00496 
00497   TH1F *pionXsMB=new TH1F(*pionYieldMB);
00498   TH1F *pionXsHT1=new TH1F(*pionYieldHT1);
00499   TH1F *pionXsHT2=new TH1F(*pionYieldHT2);
00500   pionXsMB->Scale((26.1/0.85));//times xs_bbc/bbc_eff
00501   pionXsHT1->Scale((26.1/0.85));
00502   pionXsHT2->Scale((26.1/0.85));
00503 
00504 
00505   TH1F *pionXsMBnoErr=new TH1F(*pionXsMB);
00506   TH1F *pionXsHT1noErr=new TH1F(*pionXsHT1);
00507   TH1F *pionXsHT2noErr=new TH1F(*pionXsHT2);
00508   for(int i=1;i<=pionXsMBnoErr->GetNbinsX();i++){
00509     pionXsMBnoErr->SetBinError(i,0.);
00510   }
00511   for(int i=1;i<=pionXsHT1noErr->GetNbinsX();i++){
00512     pionXsHT1noErr->SetBinError(i,0.);
00513   }
00514   for(int i=1;i<=pionXsHT2noErr->GetNbinsX();i++){
00515     pionXsHT2noErr->SetBinError(i,0.);
00516   }
00517 
00518   TGraphErrors *g_pionXsMB=new TGraphErrors(pionXsMB);
00519   g_pionXsMB->SetName("g_pionXsMB");
00520   removeThesePoints(g_pionXsMB,1);
00521 
00522   TGraphErrors *g_pionXsHT1=new TGraphErrors(pionXsHT1);
00523   g_pionXsHT1->SetName("g_pionXsHT1");
00524   removeThesePoints(g_pionXsHT1,2);
00525 
00526   TGraphErrors *g_pionXsHT2=new TGraphErrors(pionXsHT2);
00527   g_pionXsHT2->SetName("g_pionXsHT2");
00528   removeThesePoints(g_pionXsHT2,3);
00529 
00530   if(0){
00531     cout<<endl<<"xsec: x  y  ex  ey"<<endl;
00532     cout<<"minbias"<<endl;
00533     printPoints(g_pionXsMB);
00534     cout<<endl<<"hightower-1"<<endl;
00535     printPoints(g_pionXsHT1);
00536     cout<<endl<<"hightower-2"<<endl;
00537     printPoints(g_pionXsHT2);
00538     cout<<endl;
00539   }
00540   
00541 
00542   TMultiGraph *m_pions_fit=new TMultiGraph();
00543   m_pions_fit->SetName("m_pions_fit");
00544   m_pions_fit->SetMinimum(5.0e-9);
00545   m_pions_fit->SetMaximum(0.99);
00546 
00547 
00548   m_pions_fit->Add(g_pionXsMB);
00549   m_pions_fit->Add(g_pionXsHT1);
00550   m_pions_fit->Add(g_pionXsHT2);
00551 
00552   TF1 *fitQCD=new TF1("fitQCD",sumpqcd,1.,15.,6);
00553   fitQCD->SetParameters(600.,-8.2,4.,-8.5,2.,2.);
00554   fitQCD->FixParameter(4,2.);
00555   m_pions_fit->Fit(fitQCD,"R");
00556   
00557   bool inclPhenix=false;
00558   bool inclFrank=false;
00559   bool inclSasha=false;
00560   bool inclPqcd=false;
00561 
00562   TCanvas *compare=new TCanvas("compare","compare;p_{T}:xsec (mb)",600,750);
00563   compare->cd();
00564 
00565   TPad *padt=new TPad("padt","",0.0,0.3,1.,1.0);
00566   padt->SetBottomMargin(0.001);
00567   padt->SetLeftMargin(0.15);
00568   TPad *padb=new TPad("padb","",0.0,0.0,1.,0.3);
00569   padb->SetTopMargin(0.001);
00570   padb->SetBottomMargin(0.25);
00571   padb->SetLeftMargin(0.15);
00572 
00573   padt->Draw();
00574   padt->cd();
00575   gPad->SetLogy();
00576 
00577   TMultiGraph *m_pions=new TMultiGraph();
00578   m_pions->SetName("m_pions");
00579 
00580   if(inclPqcd){
00581     m_pions->Add(kkp,"c");
00582     m_pions->Add(kkp05,"c");
00583     m_pions->Add(kkp2,"c");
00584   }
00585   if(inclSasha){
00586     m_pions->Add(sasha_pp05);
00587   }
00588   if(inclFrank){
00589     m_pions->Add(frank_pp);
00590   }
00591   if(inclPhenix){
00592     m_pions->Add(phenix_pp);
00593   }
00594 
00595   //g_pionXsMB->Print();
00596   //cout<<endl<<endl;
00597   //g_pionXsHT1->Print();
00598   //cout<<endl<<endl;
00599   //g_pionXsHT2->Print();
00600   //cout<<endl;
00601 
00602   m_pions->Add(g_pionXsMB);
00603   m_pions->Add(g_pionXsHT1);
00604   m_pions->Add(g_pionXsHT2);
00605 
00606   m_pions->SetMinimum(1.0e-9);
00607   m_pions->SetMaximum(1.);
00608 
00609   m_pions->Draw("ap");
00610 
00611   //fitQCD->Draw("same");
00612 
00613   m_pions->GetXaxis()->SetLabelSize(0.);
00614   m_pions->GetYaxis()->SetTitle("Ed^{3}#sigma/dp^{3} (mb GeV^{-2}c^{2})");
00615 
00616   TLegend *leg=new TLegend(.5,.5,.85,.85);
00617 
00618   if(inclPhenix) leg->AddEntry(phenix_pp,"PHENIX p+p","p");
00619   if(inclFrank) leg->AddEntry(frank_pp,"STAR preliminary (upd.)","p");
00620   if(inclSasha) leg->AddEntry(sasha_pp05,"O.Grebenyuk p+p","p");
00621   leg->AddEntry(g_pionXsMB,"p+p minimum bias","p");
00622   leg->AddEntry(g_pionXsHT1,"hightower 1","p");
00623   leg->AddEntry(g_pionXsHT2,"hightower 2","p");
00624   if(inclPqcd){
00625     leg->AddEntry(kkp,"kkp + CTEQ6m, #mu=p_{T}","l");
00626     leg->AddEntry(kkp2,"#mu=2p_{T},p_{T}/2","l");
00627     leg->Draw("same");
00628   }
00629 
00630   leg->SetFillColor(0);
00631   leg->Draw();
00632 
00633   compare->cd();
00634   padb->Draw();
00635   padb->cd();
00636 
00637   TGraphErrors *sasha_pp05_overPqcd=new TGraphErrors(*sasha_pp05);
00638   divideGraphWithGraph(sasha_pp05_overPqcd,kkp);
00639   TGraphErrors *phenix_pp05_overPqcd=new TGraphErrors(*phenix_pp);
00640   divideGraphWithGraph(phenix_pp05_overPqcd,kkp);
00641   TGraphErrors *g_pionXsMB_overPqcd=new TGraphErrors(*g_pionXsMB);
00642   divideGraphWithGraph(g_pionXsMB_overPqcd,kkp);
00643   TGraphErrors *g_pionXsHT1_overPqcd=new TGraphErrors(*g_pionXsHT1);
00644   divideGraphWithGraph(g_pionXsHT1_overPqcd,kkp);
00645   TGraphErrors *g_pionXsHT2_overPqcd=new TGraphErrors(*g_pionXsHT2);
00646   divideGraphWithGraph(g_pionXsHT2_overPqcd,kkp);
00647   TGraphErrors *frank_pp05_overPqcd=new TGraphErrors(*frank_pp);
00648   divideGraphWithGraph(frank_pp05_overPqcd,kkp);
00649 
00650   TGraphErrors *kkp05_ratio=new TGraphErrors(*kkp05);
00651   divideGraphWithGraph(kkp05_ratio,kkp);
00652   TGraphErrors *kkp_ratio=new TGraphErrors(*kkp);
00653   divideGraphWithGraph(kkp_ratio,kkp);
00654   TGraphErrors *kkp2_ratio=new TGraphErrors(*kkp2);
00655   divideGraphWithGraph(kkp2_ratio,kkp);
00656 
00657   //systematic errors:
00658   TGraphErrors *g_pionXsMB_sys=new TGraphErrors(*g_pionXsMB_overPqcd);
00659   set_sys_pp_pion(g_pionXsMB_sys);
00660   TGraphErrors *g_pionXsHT1_sys=new TGraphErrors(*g_pionXsHT1_overPqcd);
00661   set_sys_pp_pion(g_pionXsHT1_sys);
00662   TGraphErrors *g_pionXsHT2_sys=new TGraphErrors(*g_pionXsHT2_overPqcd);
00663   set_sys_pp_pion(g_pionXsHT2_sys);
00664   
00665   TMultiGraph *m_pions_over_pqcd=new TMultiGraph();
00666 
00667   m_pions_over_pqcd->SetMinimum(0.0);
00668   m_pions_over_pqcd->SetMaximum(2.5);
00669 
00670   if(inclPqcd){
00671     m_pions_over_pqcd->Add(kkp05_ratio,"c");
00672     m_pions_over_pqcd->Add(kkp_ratio,"c");
00673     m_pions_over_pqcd->Add(kkp2_ratio,"c");
00674   }
00675   m_pions_over_pqcd->Add(g_pionXsMB_overPqcd);
00676   m_pions_over_pqcd->Add(g_pionXsHT1_overPqcd);
00677   m_pions_over_pqcd->Add(g_pionXsHT2_overPqcd);
00678   //m_pions_over_pqcd->Add(g_pionXsMB_sys,"c");
00679   //m_pions_over_pqcd->Add(g_pionXsHT1_sys,"c");
00680   //m_pions_over_pqcd->Add(g_pionXsHT2_sys,"c");
00681   if(inclPhenix) m_pions_over_pqcd->Add(phenix_pp05_overPqcd);
00682   if(inclFrank) m_pions_over_pqcd->Add(frank_pp05_overPqcd);
00683   if(inclSasha) m_pions_over_pqcd->Add(sasha_pp05_overPqcd);
00684 
00685   //m_pions_over_pqcd->Draw("ap");
00686   //m_pions_over_pqcd->GetXaxis()->SetTitle("p_{T} (GeV/c)");
00687 
00688   compare->SaveAs((dir+"pionxsec_pp.eps").Data());
00689   compare->SaveAs((dir+"pionxsec_pp.root").Data());
00690 
00691 
00692   TMultiGraph *m_pions_over_fit=new TMultiGraph();
00693   m_pions_over_fit->SetMinimum(0.01);
00694   m_pions_over_fit->SetMaximum(1.99);
00695 
00696   TGraphErrors *g_pionXsMB_overFit=new TGraphErrors(*g_pionXsMB);
00697   divideGraphWithFunction(g_pionXsMB_overFit,fitQCD);
00698   TGraphErrors *g_pionXsHT1_overFit=new TGraphErrors(*g_pionXsHT1);
00699   divideGraphWithFunction(g_pionXsHT1_overFit,fitQCD);
00700   TGraphErrors *g_pionXsHT2_overFit=new TGraphErrors(*g_pionXsHT2);
00701   divideGraphWithFunction(g_pionXsHT2_overFit,fitQCD);
00702 
00703   m_pions_over_fit->Add(g_pionXsMB_overFit);
00704   m_pions_over_fit->Add(g_pionXsHT1_overFit);
00705   m_pions_over_fit->Add(g_pionXsHT2_overFit);
00706 
00707   m_pions_over_fit->Draw("ap");
00708 
00709   compare->SaveAs((dir+"pionxsec_pp_overfit.eps").Data());
00710   compare->SaveAs((dir+"pionxsec_pp_overfit.root").Data());
00711 
00712 
00713 
00714   TCanvas *compare2=new TCanvas("compare2","compare2;p_{T};yield divided by fit",600,300);
00715   compare2->cd();
00716 
00717   //divide by fit:
00718   TGraphErrors *frank_pp2=new TGraphErrors(*frank_pp);  
00719   divideGraphWithFunction(frank_pp2,fitQCD);
00720   TGraphErrors *sasha_pp2=new TGraphErrors(*sasha_pp05);
00721   divideGraphWithFunction(sasha_pp2,fitQCD);
00722   TGraphErrors *phenix_pp2=new TGraphErrors(*phenix_pp);
00723   divideGraphWithFunction(phenix_pp2,fitQCD);
00724   TGraphErrors *g_pionXsMBcopy=new TGraphErrors(*g_pionXsMB);
00725   divideGraphWithFunction(g_pionXsMBcopy,fitQCD);
00726   TGraphErrors *g_pionXsHT1copy=new TGraphErrors(*g_pionXsHT1);
00727   divideGraphWithFunction(g_pionXsHT1copy,fitQCD);
00728   TGraphErrors *g_pionXsHT2copy=new TGraphErrors(*g_pionXsHT2);
00729   divideGraphWithFunction(g_pionXsHT2copy,fitQCD);
00730 
00731   inclPhenix=true;
00732   inclFrank=false;
00733   inclSasha=false;
00734   inclPqcd=false;
00735 
00736   TMultiGraph *m_pions2=new TMultiGraph();
00737   m_pions2->SetName("m_pions2");
00738   if(inclSasha) m_pions2->Add(sasha_pp2);
00739   if(inclFrank) m_pions2->Add(frank_pp2);
00740   if(inclPhenix) m_pions2->Add(phenix_pp2);
00741   
00742   m_pions2->Add(g_pionXsMBcopy);
00743   m_pions2->Add(g_pionXsHT1copy);
00744   m_pions2->Add(g_pionXsHT2copy);
00745 
00746   m_pions2->SetMinimum(0.000001);
00747   m_pions2->SetMaximum(3.);
00748   m_pions2->Draw("ap");
00749 
00750   TLegend *legg=new TLegend(.25,.55,.65,.85);
00751   legg->AddEntry(g_pionXsMBcopy,"minimum bias","p");
00752   legg->AddEntry(g_pionXsHT1copy,"hightower 1","p");
00753   legg->AddEntry(g_pionXsHT2copy,"hightower 2","p");
00754   if(inclFrank) legg->AddEntry(frank_pp2,"Frank's p+p (upd.)","p");
00755   if(inclSasha) legg->AddEntry(sasha_pp2,"Sasha's p+p","p");
00756   if(inclPhenix) legg->AddEntry(phenix_pp,"PHENIX p+p","p");
00757   legg->Draw("same");
00758   legg->SetFillColor(0);
00759 
00760   compare2->cd(0);
00761   compare2->SaveAs((dir+"pionxsec_pp_ratio.eps").Data());
00762   compare2->SaveAs((dir+"pionxsec_pp_ratio.root").Data());
00763 
00764 
00765   //********************************************
00766   //   Get double ratio:
00767   //********************************************
00768 
00769   //pion decay photon eff:
00770   TFile gg(eFile.Data(),"OPEN");  
00771   TH1F *effGammaMB=new TH1F(*h_effDaughtersMB);
00772   TH1F *effGammaHT1=new TH1F(*h_effDaughtersHT1);
00773   TH1F *effGammaHT2=new TH1F(*h_effDaughtersHT2);
00774   
00775   //single photon eff:
00776   TFile gg_single(eFileGamma.Data(),"OPEN");
00777   TH1F *effGammaSingleMB=new TH1F(*h_effMB);
00778   TH1F *effGammaSingleHT1=new TH1F(*h_effHT1);
00779   TH1F *effGammaSingleHT2=new TH1F(*h_effHT2);
00780 
00781   //raw neutral clusters:
00782   TFile ff(pi0File.Data(),"OPEN");
00783   TH1F *gammaYieldMB=new TH1F(*h_gammaMB);
00784   TH1F *gammaYieldHT1=new TH1F(*h_gammaHT1);
00785   TH1F *gammaYieldHT2=new TH1F(*h_gammaHT2);
00786 
00787   //divide rap. bite:
00788   gammaYieldMB->Scale(1./dy_gamma);
00789   gammaYieldHT1->Scale(1./dy_gamma);
00790   gammaYieldHT2->Scale(1./dy_gamma);
00791 
00792 
00793   for(Int_t i=1;i<=gammaYieldMB->GetNbinsX();i++){
00794     gammaYieldMB->SetBinContent(i,gammaYieldMB->GetBinContent(i)/gammaYieldMB->GetXaxis()->GetBinWidth(i));
00795     gammaYieldMB->SetBinError(i,gammaYieldMB->GetBinError(i)/gammaYieldMB->GetXaxis()->GetBinWidth(i));
00796   }
00797   gammaYieldMB->Scale(psMB/(psMB*numberOfMB*2.*TMath::Pi()));
00798   gammaYieldMB->Divide(scale);// ../pT
00799   gammaYieldMB->Multiply(h_binmb);
00800   gammaYieldMB->Divide(effGammaMB);
00801   gammaYieldMB->Multiply(gamma_cpv_corrMB);
00802   gammaYieldMB->Multiply(gamma_cont_corrMB);
00803   gammaYieldMB->Multiply(gamma_conv_corrMB);
00804 
00805   gammaYieldMB->SetMarkerStyle(8);
00806   gammaYieldMB->SetMarkerSize(1.);
00807   TGraphErrors *g_inclPhotonsMB=new TGraphErrors(gammaYieldMB);
00808 
00809 
00810   for(Int_t i=1;i<=gammaYieldHT1->GetNbinsX();i++){
00811     gammaYieldHT1->SetBinContent(i,gammaYieldHT1->GetBinContent(i)/gammaYieldHT1->GetXaxis()->GetBinWidth(i));
00812     gammaYieldHT1->SetBinError(i,gammaYieldHT1->GetBinError(i)/gammaYieldHT1->GetXaxis()->GetBinWidth(i));
00813   }
00814   gammaYieldHT1->Scale(psHT1/(psMB*nMBwithHT*2.*TMath::Pi()));
00815   gammaYieldHT1->Divide(scale);// ../pT
00816   gammaYieldHT1->Multiply(h_binht1);
00817   gammaYieldHT1->Divide(effGammaHT1);
00818   gammaYieldHT1->Multiply(gamma_cpv_corrHT1);
00819   gammaYieldHT1->Multiply(gamma_cont_corrHT1);
00820   gammaYieldHT1->Multiply(gamma_conv_corrHT1);
00821 
00822   gammaYieldHT1->SetMarkerStyle(8);
00823   gammaYieldHT1->SetMarkerSize(1.);
00824   gammaYieldHT1->SetMarkerColor(4);
00825   TGraphErrors *g_inclPhotonsHT1=new TGraphErrors(gammaYieldHT1);
00826 
00827   for(Int_t i=1;i<=gammaYieldHT2->GetNbinsX();i++){
00828     gammaYieldHT2->SetBinContent(i,gammaYieldHT2->GetBinContent(i)/gammaYieldHT2->GetXaxis()->GetBinWidth(i));
00829     gammaYieldHT2->SetBinError(i,gammaYieldHT2->GetBinError(i)/gammaYieldHT2->GetXaxis()->GetBinWidth(i));
00830   }
00831   gammaYieldHT2->Scale(psHT2/(psMB*nMBwithHT*2.*TMath::Pi()));
00832   gammaYieldHT2->Divide(scale);// ../pT
00833   gammaYieldHT2->Multiply(h_binht2);
00834   gammaYieldHT2->Divide(effGammaHT2);
00835   gammaYieldHT2->Multiply(gamma_cpv_corrHT2);
00836   gammaYieldHT2->Multiply(gamma_cont_corrHT2);
00837   gammaYieldHT2->Multiply(gamma_conv_corrHT2);
00838 
00839   gammaYieldHT2->SetMarkerStyle(8);
00840   gammaYieldHT2->SetMarkerSize(1.);
00841   gammaYieldHT2->SetMarkerColor(2);
00842   TGraphErrors *g_inclPhotonsHT2=new TGraphErrors(gammaYieldHT2);
00843 
00844   removeThesePoints(g_inclPhotonsMB,1);
00845   removeThesePoints(g_inclPhotonsHT1,2);
00846   removeThesePoints(g_inclPhotonsHT2,2);
00847 
00848   TMultiGraph *m_incl=new TMultiGraph();
00849   m_incl->SetName("m_incl");
00850   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}}");
00851   m_incl->Add(g_inclPhotonsMB);
00852   m_incl->Add(g_inclPhotonsHT1);
00853   m_incl->Add(g_inclPhotonsHT2);
00854 
00855   m_incl->Fit(fitQCD,"R0");
00856 
00857   m_incl->SetMinimum(1.e-11);
00858   m_incl->SetMaximum(10.);
00859   TCanvas *c_incl=new TCanvas("c_incl","c_incl",600,400);
00860   gPad->SetLogy();
00861   m_incl->Draw("ap");
00862   c_incl->SaveAs((dir+"inclPhotonYield.eps").Data());
00863   c_incl->SaveAs((dir+"inclPhotonYield.root").Data());
00864 
00865 
00866 
00867   //get ratio:
00868   TH1F *gammaYieldMBratio=new TH1F(*gammaYieldMB);
00869   gammaYieldMBratio->SetName("gammaYieldMBratio");
00870   TH1F *gammaYieldHT1ratio=new TH1F(*gammaYieldHT1);
00871   gammaYieldHT1ratio->SetName("gammaYieldHT1ratio");
00872   TH1F *gammaYieldHT2ratio=new TH1F(*gammaYieldHT2);
00873   gammaYieldHT2ratio->SetName("gammaYieldHT2ratio");
00874 
00875   gammaYieldMBratio->Divide(pionYieldMBratio);
00876   gammaYieldHT1ratio->Divide(pionYieldHT1ratio);
00877   gammaYieldHT2ratio->Divide(pionYieldHT2ratio);
00878 
00879 
00880   //correct gamma over pion ratio, using two efficiencies:
00881   getRatio(gammaYieldMBratio,effGammaMB,effGammaSingleMB,fit_piondecay);
00882   getRatio(gammaYieldHT1ratio,effGammaHT1,effGammaSingleHT1,fit_piondecay);
00883   getRatio(gammaYieldHT2ratio,effGammaHT2,effGammaSingleHT2,fit_piondecay);
00884 
00885   TH1F *gammaYieldMBratio_incl=new TH1F(*gammaYieldMBratio);
00886   TH1F *gammaYieldHT1ratio_incl=new TH1F(*gammaYieldHT1ratio);
00887   TH1F *gammaYieldHT2ratio_incl=new TH1F(*gammaYieldHT2ratio);
00888   
00889   TH1F *gammaYieldMBratioNoErr=new TH1F(*gammaYieldMBratio);
00890   TH1F *gammaYieldHT1ratioNoErr=new TH1F(*gammaYieldHT1ratio);
00891   TH1F *gammaYieldHT2ratioNoErr=new TH1F(*gammaYieldHT2ratio);
00892   for(int i=1;i<=gammaYieldMBratioNoErr->GetNbinsX();i++){
00893     gammaYieldMBratioNoErr->SetBinError(i,0.);
00894   }
00895   for(int i=1;i<=gammaYieldHT1ratioNoErr->GetNbinsX();i++){
00896     gammaYieldHT1ratioNoErr->SetBinError(i,0.);
00897   }
00898   for(int i=1;i<=gammaYieldHT2ratioNoErr->GetNbinsX();i++){
00899     gammaYieldHT2ratioNoErr->SetBinError(i,0.);
00900   }
00901 
00902   TGraphErrors *g_ratioMB=new TGraphErrors(gammaYieldMBratio);
00903   g_ratioMB->SetName("g_ratioMB");
00904   TGraphErrors *g_ratioHT1=new TGraphErrors(gammaYieldHT1ratio);
00905   g_ratioHT1->SetName("g_ratioHT1");
00906   TGraphErrors *g_ratioHT2=new TGraphErrors(gammaYieldHT2ratio);
00907   g_ratioHT2->SetName("g_ratioHT2");
00908 
00909 
00910   removeThesePoints(g_ratioMB,1);
00911   removeThesePoints(g_ratioHT1,2);
00912   removeThesePoints(g_ratioHT2,3);
00913 
00914 
00915   TCanvas *c_ratio=new TCanvas("c_ratio","c_ratio",400,200);
00916 
00917   TMultiGraph *m_ratio=new TMultiGraph("m_ratio","p+p 2005;p_{T};#gamma/#pi^{0}");
00918   m_ratio->Add(g_ratioMB);
00919   m_ratio->Add(g_ratioHT1);
00920   m_ratio->Add(g_ratioHT2);
00921 
00922   m_ratio->Draw("ap");
00923   m_ratio->SetMinimum(.001);
00924   m_ratio->SetMaximum(1.5);
00925 
00926   TLegend *leg3=new TLegend(.35,.65,.65,.85);
00927   leg3->AddEntry(g_ratioMB,"minimum bias","p");
00928   leg3->AddEntry(g_ratioHT1,"hightower 1","p");
00929   leg3->AddEntry(g_ratioHT2,"hightower 2","p");
00930   leg3->AddEntry(fit_decay,"decay background (total)","l");
00931   leg3->AddEntry(fit_piondecay,"decay background (#pi^{0})","l");
00932   leg3->SetFillColor(0);
00933   leg3->Draw("same");
00934 
00935   fit_decay->SetLineColor(13);
00936   fit_decay->SetLineWidth(1);
00937   fit_decay->SetLineColor(1);
00938   fit_decay->Draw("same");
00939   fit_piondecay->SetLineColor(13);
00940   fit_piondecay->SetLineWidth(1);
00941   fit_piondecay->SetLineStyle(2);
00942   fit_piondecay->SetLineColor(1);
00943   fit_piondecay->Draw("same");
00944 
00945   c_ratio->SaveAs((dir+"gammaOverPion.eps").Data());
00946   c_ratio->SaveAs((dir+"gammaOverPion.root").Data());
00947 
00948 
00949   //create fully corrected incl. photons:
00950   gammaYieldMBratio_incl->Multiply(pionYieldMBratio);
00951   gammaYieldHT1ratio_incl->Multiply(pionYieldHT1ratio);
00952   gammaYieldHT2ratio_incl->Multiply(pionYieldHT2ratio);
00953   //gammaYieldMBratio_incl->Scale(7.5);
00954   //gammaYieldHT1ratio_incl->Scale(7.5);
00955   //gammaYieldHT2ratio_incl->Scale(7.5);
00956   
00957   TGraphErrors *g_incl_corrMB=new TGraphErrors(gammaYieldMBratio_incl);
00958   TGraphErrors *g_incl_corrHT1=new TGraphErrors(gammaYieldHT1ratio_incl);
00959   TGraphErrors *g_incl_corrHT2=new TGraphErrors(gammaYieldHT2ratio_incl);
00960 
00961   TCanvas *c_incl_corr=new TCanvas("c_incl_corr","c_incl_corr",400,300);
00962   gPad->SetLogy();
00963   TMultiGraph *m_incl_corr=new TMultiGraph();
00964   m_incl_corr->Add(g_incl_corrMB);
00965   m_incl_corr->Add(g_incl_corrHT1);
00966   m_incl_corr->Add(g_incl_corrHT2);
00967 
00968   m_incl_corr->SetMinimum(1.e-11);
00969   m_incl_corr->SetMaximum(1.);
00970 
00971   m_incl_corr->Draw("apX");
00972   c_incl_corr->SaveAs((dir+"inclPhotonYieldCorr.eps").Data());
00973   c_incl_corr->SaveAs((dir+"inclPhotonYieldCorr.root").Data());
00974 
00975 
00976 
00977   TCanvas *c_doubleratio=new TCanvas("c_doubleratio","c_doubleratio",400,300);
00978   gStyle->SetOptStat(0);
00979   c_doubleratio->cd(1);
00980 
00981   TH1F *gammaYieldMBdoubleratio=new TH1F(*gammaYieldMBratio);
00982   TH1F *gammaYieldHT1doubleratio=new TH1F(*gammaYieldHT1ratio);
00983   TH1F *gammaYieldHT2doubleratio=new TH1F(*gammaYieldHT2ratio);
00984 
00985   gammaYieldMBdoubleratio->Divide(fit_decay);
00986   gammaYieldHT1doubleratio->Divide(fit_decay);
00987   gammaYieldHT2doubleratio->Divide(fit_decay);
00988 
00989   TGraphErrors *g_doubleRatioMB=new TGraphErrors(gammaYieldMBdoubleratio);
00990   g_doubleRatioMB->SetName("g_doubleRatioMB");
00991   g_doubleRatioMB->SetMarkerStyle(8);
00992   TGraphErrors *g_doubleRatioHT1=new TGraphErrors(gammaYieldHT1doubleratio);
00993   g_doubleRatioHT1->SetName("g_doubleRatioHT1");
00994   g_doubleRatioHT1->SetMarkerStyle(8);
00995   TGraphErrors *g_doubleRatioHT2=new TGraphErrors(gammaYieldHT2doubleratio);
00996   g_doubleRatioHT2->SetName("g_doubleRatioHT2");
00997   g_doubleRatioHT2->SetMarkerStyle(8);
00998 
00999   removeThesePoints(g_doubleRatioMB,1);
01000   removeThesePoints(g_doubleRatioHT1,2);
01001   removeThesePoints(g_doubleRatioHT2,3);
01002 
01003   TMultiGraph *m_doubleratio=new TMultiGraph();
01004   m_doubleratio->SetName("m_doubleratio");
01005   m_doubleratio->SetMinimum(.5);
01006   m_doubleratio->SetMaximum(2.75);
01007 
01008   cout<<endl;
01009   g_doubleRatioHT1->Print();
01010   cout<<endl<<endl;
01011   g_doubleRatioHT2->Print();
01012   cout<<endl;
01013 
01014   //m_doubleratio->Add(g_doubleRatioMB,"p");
01015   m_doubleratio->Add(g_doubleRatioHT1,"p");
01016   m_doubleratio->Add(g_doubleRatioHT2,"p");
01017 
01018   g_photonpqcd->SetLineWidth(2);
01019   g_photonpqcd->SetLineColor(2);
01020   g_photonpqcd05->SetLineWidth(2);
01021   g_photonpqcd05->SetLineColor(2);
01022   g_photonpqcd05->SetLineStyle(2);
01023   g_photonpqcd2->SetLineWidth(2);
01024   g_photonpqcd2->SetLineColor(2);
01025   g_photonpqcd2->SetLineStyle(2);
01026 
01027   m_doubleratio->Add(g_photonpqcd,"c");
01028   m_doubleratio->Add(g_photonpqcd05,"c");
01029   m_doubleratio->Add(g_photonpqcd2,"c");
01030 
01031   //appropriate fit to photon pqcd result
01032   TF1 *fitGamma2=new TF1("fitGamma2","1.+[0]*pow(x,[1])",2.,15.);
01033   g_photonpqcd->Fit(fitGamma2,"R0");
01034 
01035   m_doubleratio->Draw("a");
01036   
01037   m_doubleratio->GetXaxis()->SetTitle("p_{T} (GeV/c)");
01038   m_doubleratio->GetYaxis()->SetTitle("1 + #gamma_{dir}/#gamma_{incl}");
01039   m_doubleratio->GetXaxis()->SetRangeUser(2.,16.);
01040   
01041 
01042   TLegend *leg5=new TLegend(.15,.6,.6,.8);
01043   leg5->AddEntry(g_doubleRatioHT1,"hightower-1","p");
01044   leg5->AddEntry(g_doubleRatioHT2,"hightower-2","p");
01045   leg5->AddEntry(g_photonpqcd,"NLO (CTEQ6+KKP) #mu=p_{T}","l");
01046   leg5->AddEntry(g_photonpqcd05,"#mu=2p_{T}, #mu=p_{T}/2","l");
01047   leg5->SetFillColor(0);
01048   leg5->Draw("same");
01049 
01050   c_doubleratio->cd(0);  
01051   c_doubleratio->SaveAs((dir+"gammaDoubleRatio.eps").Data());
01052   c_doubleratio->SaveAs((dir+"gammaDoubleRatio.root").Data());
01053 
01054 
01055   TCanvas *c_dirphoton=new TCanvas("c_dirphoton","c_dirphoton",400,300);
01056   gStyle->SetOptStat(0);
01057   c_dirphoton->cd(1);
01058 
01059   TH1F *dirphotonYieldMB=new TH1F(*gammaYieldMBdoubleratio);
01060   TH1F *dirphotonYieldHT1=new TH1F(*gammaYieldHT1doubleratio);
01061   TH1F *dirphotonYieldHT2=new TH1F(*gammaYieldHT2doubleratio);
01062 
01063   TH1F *dirphotonYieldMBnoErr=new TH1F(*dirphotonYieldMB);
01064   for(int i=1;i<=dirphotonYieldMBnoErr->GetNbinsX();i++){
01065     dirphotonYieldMBnoErr->SetBinError(i,0.);
01066   }
01067   TH1F *dirphotonYieldHT1noErr=new TH1F(*dirphotonYieldHT1);
01068   for(int i=1;i<=dirphotonYieldHT1noErr->GetNbinsX();i++){
01069     dirphotonYieldHT1noErr->SetBinError(i,0.);
01070   }
01071   TH1F *dirphotonYieldHT2noErr=new TH1F(*dirphotonYieldHT2);
01072   for(int i=1;i<=dirphotonYieldHT1noErr->GetNbinsX();i++){
01073     dirphotonYieldHT2noErr->SetBinError(i,0.);
01074   }
01075 
01076   TF1 *f_unity=new TF1("f_unity","1.",0.,15.);
01077   dirphotonYieldMB->Add(f_unity,-1.);
01078   dirphotonYieldHT1->Add(f_unity,-1.);
01079   dirphotonYieldHT2->Add(f_unity,-1.);
01080 
01081 
01082   dirphotonYieldMB->Divide(dirphotonYieldMBnoErr);
01083   dirphotonYieldHT1->Divide(dirphotonYieldHT1noErr);
01084   dirphotonYieldHT2->Divide(dirphotonYieldHT2noErr);
01085   dirphotonYieldMB->Multiply(gammaYieldMBratioNoErr);
01086   dirphotonYieldHT1->Multiply(gammaYieldHT1ratioNoErr);
01087   dirphotonYieldHT2->Multiply(gammaYieldHT2ratioNoErr);
01088   dirphotonYieldMB->Multiply(pionXsMBnoErr);
01089   dirphotonYieldHT1->Multiply(pionXsHT1noErr);
01090   dirphotonYieldHT2->Multiply(pionXsHT2noErr);
01091 
01092 
01093   TGraphErrors *g_dirphotonMB=new TGraphErrors(dirphotonYieldMB);
01094   g_dirphotonMB->SetName("g_dirphotonMB");
01095   g_dirphotonMB->SetMarkerStyle(8);
01096   TGraphErrors *g_dirphotonHT1=new TGraphErrors(dirphotonYieldHT1);
01097   g_dirphotonHT1->SetName("g_dirphotonHT1");
01098   g_dirphotonHT1->SetMarkerStyle(8);
01099   TGraphErrors *g_dirphotonHT2=new TGraphErrors(dirphotonYieldHT2);
01100   g_dirphotonHT2->SetName("g_dirphotonHT2");
01101   g_dirphotonHT2->SetMarkerStyle(8);
01102 
01103 
01104   removeThesePoints(g_dirphotonMB,1);
01105   removeThesePoints(g_dirphotonHT1,2);
01106   removeThesePoints(g_dirphotonHT2,3);
01107 
01108   gPad->SetLogy();
01109 
01110   TMultiGraph *m_dirphoton=new TMultiGraph();
01111   m_dirphoton->SetName("m_dirphoton");
01112   m_dirphoton->SetMinimum(1.0e-11);
01113   m_dirphoton->SetMaximum(0.1);
01114 
01115   m_dirphoton->Add(g_dirgamma,"c");
01116   m_dirphoton->Add(g_dirgamma05,"c");
01117   m_dirphoton->Add(g_dirgamma2,"c");
01118 
01119   cout<<"direct photons:"<<endl;
01120   g_dirphotonHT1->Print();
01121   cout<<endl;
01122   g_dirphotonHT2->Print();
01123   cout<<endl;
01124 
01125   m_dirphoton->Add(g_dirphotonHT1,"p");
01126   m_dirphoton->Add(g_dirphotonHT2,"p");
01127 
01128   m_dirphoton->Draw("a");
01129   
01130   m_dirphoton->GetXaxis()->SetTitle("p_{T} (GeV/c)");
01131   m_dirphoton->GetYaxis()->SetTitle("1 - R^{-1}");
01132   m_dirphoton->GetXaxis()->SetRangeUser(2.,16.);
01133   
01134 
01135   TLegend *leg5=new TLegend(.15,.6,.6,.8);
01136   leg5->AddEntry(g_dirphotonHT1,"hightower-1","p");
01137   leg5->AddEntry(g_dirphotonHT2,"hightower-2","p");
01138   leg5->SetFillColor(0);
01139   leg5->Draw("same");
01140 
01141   c_dirphoton->cd(0);  
01142   c_dirphoton->SaveAs((dir+"gammaDirPhoton.eps").Data());
01143   c_dirphoton->SaveAs((dir+"gammaDirPhoton.root").Data());
01144 
01145 
01146 
01147 
01148 
01149   return;
01150 }