00001 #include <TMath.h>
00002 #include <TGraphErrors.h>
00003 #include <TH1F.h>
00004 #include <TH2F.h>
00005 #include <TStyle.h>
00006 #include <TF1.h>
00007 #include <TCanvas.h>
00008 #include <TFile.h>
00009 #include <TPad.h>
00010 #include <TLegend.h>
00011 #include <TMultiGraph.h>
00012 #include <TColor.h>
00013 #include <assert.h>
00014
00015 #include <fstream>
00016 using namespace std;
00017
00018 #include <StEmcPool/StPhotonCommon/PhotonAnalysisSettings.h>
00019
00020 #include "macros/systematics/systematics_rda.h"
00021 #include "macros/systematics/systematics_dau.h"
00022 #include "macros/systematics/systematics_pp.h"
00023
00024 #include "AnaCuts.h"
00025 #include "Pi0Analysis.h"
00026
00027 #include "PhotonAnalysisUtil.h"
00028
00029 ClassImp(PhotonAnalysisUtil)
00030
00031 double pqcdlowpt(double *x,double *par)
00032 {
00033 double ret=par[0]*TMath::Power(1.+x[0],par[1]);
00034 return ret;
00035 }
00036 double pqcdhighpt(double *x,double *par)
00037 {
00038 double ret=par[0]*TMath::Power(1.+x[0],par[1]);
00039 return ret;
00040 }
00041 double sumpqcd(double *x,double *par)
00042 {
00043
00044
00045 double SW=1. - 1./(1.+TMath::Exp(TMath::Abs(par[4])*(x[0]-par[5])));
00046 double ret=SW*pqcdhighpt(x,&par[2]);
00047 ret+=(1.-SW)*pqcdlowpt(x,par);
00048
00049 return ret;
00050 }
00051 void printPoints(TGraphErrors *g)
00052 {
00053 for(Int_t i=0;i<g->GetN();i++){
00054 double x,y,ey;
00055 g->GetPoint(i,x,y);
00056 ey=g->GetErrorY(i);
00057 cout<<x<<" "<<y<<" "<<ey<<endl;
00058 }
00059 }
00060 void divideGraphWithFunction(TGraphErrors *graph,TF1 *func)
00061 {
00062 for(int i=0;i<graph->GetN();i++){
00063 double x=0.;
00064 double y=0.;
00065 graph->GetPoint(i,x,y);
00066 double ey=graph->GetErrorY(i);
00067 graph->SetPoint(i,x,y/func->Eval(x));
00068 graph->SetPointError(i,0.,ey/func->Eval(x));
00069 }
00070 }
00071 void divideGraphWithGraph(TGraphErrors *graph,TGraph *denom)
00072 {
00073 for(int i=0;i<graph->GetN();i++){
00074 double x=0.;
00075 double y=0.;
00076 graph->GetPoint(i,x,y);
00077 double ey=graph->GetErrorY(i);
00078 graph->SetPoint(i,x,y/denom->Eval(x));
00079 graph->SetPointError(i,0.,ey/denom->Eval(x));
00080 }
00081 }
00082 void getRatio(TH1F* h_ratio,TH1F* h_eff,TH1F* h_effSingle,TF1 *f_piondecay)
00083 {
00084
00085 TH1F *h_eff_copy=new TH1F(*h_eff);
00086 for(int i=1;i<=h_eff->GetNbinsX();i++){
00087 h_eff->SetBinError(i,0.);
00088 }
00089 TH1F *h_effSingle_copy=new TH1F(*h_effSingle);
00090 for(int i=1;i<=h_effSingle->GetNbinsX();i++){
00091 h_effSingle_copy->SetBinError(i,0.);
00092 }
00093
00094 h_ratio->Add(f_piondecay,-1.);
00095 for(int i=1;i<=h_ratio->GetNbinsX();i++){
00096
00097
00098 }
00099 h_eff_copy->Divide(h_effSingle_copy);
00100 h_ratio->Multiply(h_eff_copy);
00101 h_ratio->Add(f_piondecay);
00102 }
00103 void removeThesePoints(TGraphErrors *g,int trig)
00104 {
00105 if(trig==1){
00106 g->RemovePoint(0);
00107 g->RemovePoint(0);
00108
00109 g->RemovePoint(g->GetN()-1);
00110 }
00111 if(trig==2){
00112 g->RemovePoint(0);
00113 g->RemovePoint(0);
00114 g->RemovePoint(0);
00115 g->RemovePoint(0);
00116 g->RemovePoint(g->GetN()-1);
00117 g->RemovePoint(g->GetN()-1);
00118 }
00119 if(trig==3){
00120 g->RemovePoint(0);
00121 g->RemovePoint(0);
00122 g->RemovePoint(0);
00123 g->RemovePoint(0);
00124 g->RemovePoint(0);
00125 }
00126 }
00127
00128 void getPhotonSpectrum(const PhotonAnalysisSettings &settings) {
00129 gStyle->SetOptStat(0);
00130 gStyle->SetOptFit(0);
00131
00132 gStyle->SetErrorX(0);
00133
00134
00135
00136
00137
00138
00139
00140 double B[]={0.3,0.23,0.072,0.061};
00141
00142 double T[]={0.205,0.215,0.179,0.173};
00143
00144 double n[]={11.00,12.55,10.87,10.49};
00145
00146 double BR=35.8/63.9;
00147 double c_feeddown=0.8+0.2*BR;
00148 double CPVloss=1;
00149 TF1 *f_nbar=new TF1("f_nbar","[3]*[4]*[0]/TMath::Power((1.+(sqrt(x*x+1.) - 1.)/([1]*[2])),[2])",0.,15.);
00150 assert(f_nbar);
00151 if (settings.name == "pp05") f_nbar->SetParameters(B[3],T[3],n[3],c_feeddown,CPVloss);
00152 if (settings.name == "dAu") f_nbar->SetParameters(B[1],T[1],n[1],c_feeddown,CPVloss);
00153
00154
00155 ifstream pQCDphotons(TString(settings.input_datapoints_dir + "/pQCD_Werner/rhic_cteq6_gamma_inv_sc1.dat").Data());
00156 float ppx[100];
00157 float ppy[100];
00158 Int_t iii=0;
00159 cout<<"pqcd photons:"<<endl;
00160 while(iii<28){
00161 if(!pQCDphotons.good()) break;
00162 float dummy=0.;
00163 pQCDphotons>>ppx[iii]>>dummy>>dummy>>ppy[iii];
00164 if (settings.name == "pp05") ppy[iii]*=1.e-09;
00165 if (settings.name == "dAu") ppy[iii]*=1.e-09*2.21e3*7.5/42.0;
00166 iii++;
00167 }
00168 TGraph *g_dirgamma=new TGraph(iii,ppx,ppy);
00169 assert(g_dirgamma);
00170
00171 ifstream pQCDphotons05(TString(settings.input_datapoints_dir + "/pQCD_Werner/rhic_cteq6_gamma_inv_sc05.dat").Data());
00172 float ppx05[100];
00173 float ppy05[100];
00174 Int_t iii05=0;
00175 while(iii05<28){
00176 if(!pQCDphotons05.good()) break;
00177 float dummy=0.;
00178 pQCDphotons05>>ppx05[iii05]>>dummy>>dummy>>ppy05[iii05];
00179 if (settings.name == "pp05") ppy05[iii05]*=1.e-09;
00180 if (settings.name == "dAu") ppy05[iii05]*=1.e-09*2.21e3*7.5/42.0;
00181 iii05++;
00182 }
00183 TGraph *g_dirgamma05=new TGraph(iii05,ppx05,ppy05);
00184 assert(g_dirgamma05);
00185
00186 ifstream pQCDphotons2(TString(settings.input_datapoints_dir + "/pQCD_Werner/rhic_cteq6_gamma_inv_sc2.dat").Data());
00187 float ppx2[100];
00188 float ppy2[100];
00189 Int_t iii2=0;
00190 while(iii2<28){
00191 if(!pQCDphotons2.good()) break;
00192 float dummy=0.;
00193 pQCDphotons2>>ppx2[iii2]>>dummy>>dummy>>ppy2[iii2];
00194 if (settings.name == "pp05") ppy2[iii2]*=1.e-09;
00195 if (settings.name == "dAu") ppy2[iii2]*=1.e-09*2.21e3*7.5/42.0;
00196 iii2++;
00197 }
00198 TGraph *g_dirgamma2=new TGraph(iii2,ppx2,ppy2);
00199 assert(g_dirgamma2);
00200
00201 TString phenixFile;
00202 if (settings.name == "pp05") phenixFile = settings.input_datapoints_dir + "/phenix_xsec_pp.dat";
00203 if (settings.name == "dAu") phenixFile = settings.input_datapoints_dir + "/phenix.dat";
00204 ifstream phenix(phenixFile.Data());
00205 float phex[100];
00206 float phey[100];
00207 float ephey[100];
00208 Int_t iphe=0;
00209 while(iphe<16){
00210 if(!phenix.good()) break;
00211 if (settings.name == "pp05") phenix>>phex[iphe]>>phey[iphe];
00212 if (settings.name == "dAu") phenix>>phex[iphe]>>phey[iphe]>>ephey[iphe];
00213
00214 iphe++;
00215 }
00216 TGraphErrors *phenix_pp=new TGraphErrors(iphe,phex,phey);
00217 assert(phenix_pp);
00218 phenix_pp->SetMarkerStyle(24);
00219 phenix_pp->SetLineWidth(6);
00220 phenix_pp->SetName("phenix");
00221
00222
00223 ifstream frank(TString(settings.input_datapoints_dir + "/frank_pp05_new.dat").Data());
00224 float frx[100];
00225 float fry[100];
00226 float frex[100];
00227 float frey[100];
00228 Int_t ifr=0;
00229 while(ifr<12){
00230 if(!frank.good()) break;
00231 frank>>frx[ifr]>>fry[ifr]>>frey[ifr];
00232 frex[ifr]=0.;
00233 ifr++;
00234 }
00235 TGraphErrors *frank_pp=new TGraphErrors(ifr,frx,fry,frex,frey);
00236 assert(frank_pp);
00237 frank_pp->SetMarkerStyle(25);
00238 frank_pp->SetMarkerColor(4);
00239 frank_pp->SetName("frank_pp");
00240
00241
00242 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 };
00243 float xe_pp2005_sasha[]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
00244 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};
00245 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};
00246
00247 TGraphErrors *sasha_pp05=new TGraphErrors(13,x_pp2005_sasha,y_pp2005_sasha,xe_pp2005_sasha,ye_pp2005_sasha);
00248 assert(sasha_pp05);
00249 sasha_pp05->SetMarkerStyle(8);
00250 sasha_pp05->SetMarkerColor(TColor::GetColor(8,80,8));
00251 sasha_pp05->SetName("sasha_pp05");
00252
00253
00254 ifstream pQCDpions(TString(settings.input_datapoints_dir + "/pQCD_Werner/klaus_pi0inv_200_kkp_1.dat").Data());
00255 float pionx[100];
00256 float piony[100];
00257 Int_t ipion=0;
00258 while(ipion<28){
00259 if(!pQCDpions.good()) break;
00260 pQCDpions>>pionx[ipion]>>piony[ipion];
00261 if (settings.name == "pp05") piony[iii]*=1.e-09;
00262 if (settings.name == "dAu") piony[iii]*=1.e-09*2.21e3*7.5/42.0;
00263 ipion++;
00264 }
00265 TGraphErrors *kkp=new TGraphErrors(ipion,pionx,piony);
00266 assert(kkp);
00267 kkp->SetLineColor(54);
00268 kkp->SetName("kkp");
00269
00270
00271 ifstream pQCDpions05(TString(settings.input_datapoints_dir + "/pQCD_Werner/klaus_pi0inv_200_kkp_05.dat").Data());
00272 float pionx05[100];
00273 float piony05[100];
00274 Int_t ipion05=0;
00275 while(ipion05<28){
00276 if(!pQCDpions05.good()) break;
00277 pQCDpions05>>pionx05[ipion05]>>piony05[ipion05];
00278 if (settings.name == "pp05") piony05[iii05]*=1.e-09;
00279 if (settings.name == "dAu") piony05[iii05]*=1.e-09*2.21e3*7.5/42.0;
00280 ipion05++;
00281 }
00282 TGraphErrors *kkp05=new TGraphErrors(ipion05,pionx05,piony05);
00283 assert(kkp05);
00284 kkp05->SetLineStyle(2);
00285 kkp05->SetLineColor(54);
00286 kkp05->SetName("kkp05");
00287
00288
00289 ifstream pQCDpions2(TString(settings.input_datapoints_dir + "/pQCD_Werner/klaus_pi0inv_200_kkp_2.dat").Data());
00290 float pionx2[100];
00291 float piony2[100];
00292 Int_t ipion2=0;
00293 while(ipion2<28){
00294 if(!pQCDpions2.good()) break;
00295 pQCDpions2>>pionx2[ipion2]>>piony2[ipion2];
00296 if (settings.name == "pp05") piony2[iii2]*=1.e-09;
00297 if (settings.name == "dAu") piony2[iii2]*=1.e-09*2.21e3*7.5/42.0;
00298 ipion2++;
00299 }
00300 TGraphErrors *kkp2=new TGraphErrors(ipion2,pionx2,piony2);
00301 assert(kkp2);
00302 kkp2->SetLineStyle(2);
00303 kkp2->SetLineColor(54);
00304 kkp2->SetName("kkp2");
00305
00306 TFile f_decaybg(settings.input_decaybackground_file,"READ");
00307 const TH1F *h_decaybg_read=(const TH1F*)f_decaybg.Get("gamma");
00308 const TH1F *h_decaypion_read=(const TH1F*)f_decaybg.Get("gamma_pion");
00309 assert(h_decaybg_read);
00310 assert(h_decaypion_read);
00311 TH1F *h_decaybg = new TH1F(*h_decaybg_read);
00312 TH1F *h_decaypion = new TH1F(*h_decaypion_read);
00313 assert(h_decaybg);
00314 assert(h_decaypion);
00315 f_decaybg.Close();
00316 h_decaybg_read = 0;
00317 h_decaypion_read = 0;
00318
00319 TF1 *fit_decay=new TF1("fit_decay","[0]/TMath::Power(x,[1])+[2]",.3,15.);
00320 TF1 *fit_piondecay=new TF1(*fit_decay);
00321 assert(fit_decay);
00322 assert(fit_piondecay);
00323 fit_decay->SetParameters(1.,1.,.5);
00324 fit_piondecay->SetParameters(1.,1.,.5);
00325 h_decaybg->Fit(fit_decay,"R0Q");
00326 h_decaypion->Fit(fit_piondecay,"R0Q");
00327
00328
00329 for(Int_t i=0;i<iii;i++){
00330 ppy[i]=ppy[i]/piony[i];
00331 ppy[i]=ppy[i]/fit_decay->Eval(ppx[i]);
00332 ppy[i]+=1.;
00333 ppy05[i]=ppy05[i]/piony05[i];
00334 ppy05[i]=ppy05[i]/fit_decay->Eval(ppx05[i]);
00335 ppy05[i]+=1.;
00336 ppy2[i]=ppy2[i]/piony2[i];
00337 ppy2[i]=ppy2[i]/fit_decay->Eval(ppx2[i]);
00338 ppy2[i]+=1.;
00339 }
00340 TGraphErrors *g_photonpqcd=new TGraphErrors(iii,ppx,ppy);
00341 assert(g_photonpqcd);
00342 g_photonpqcd->SetName("g_photonpqcd");
00343 TGraphErrors *g_photonpqcd05=new TGraphErrors(iii05,ppx05,ppy05);
00344 assert(g_photonpqcd05);
00345 g_photonpqcd05->SetName("g_photonpqcd05");
00346 TGraphErrors *g_photonpqcd2=new TGraphErrors(iii2,ppx2,ppy2);
00347 assert(g_photonpqcd2);
00348 g_photonpqcd2->SetName("g_photonpqcd2");
00349
00350 AnaCuts *cuts=new AnaCuts(settings.name.Data());
00351 assert(cuts);
00352
00353
00354 TFile f(settings.input_pion_file.Data(),"READ");
00355 const TH2F *h_minvMB = (const TH2F*)f.Get("h_minvMB");
00356 const TH2F *h_minvHT1 = (const TH2F*)f.Get("h_minvHT1");
00357 const TH2F *h_minvHT2 = (const TH2F*)f.Get("h_minvHT2");
00358 const TH1F *h_events = (const TH1F*)f.Get("h_events");
00359 assert(h_minvMB);
00360 assert(h_minvHT1);
00361 assert(h_minvHT2);
00362 assert(h_events);
00363 TH2F *h_mb=new TH2F(*h_minvMB);
00364 TH2F *h_ht1=new TH2F(*h_minvHT1);
00365 TH2F *h_ht2=new TH2F(*h_minvHT2);
00366 TH1F *h_ev=new TH1F(*h_events);
00367 assert(h_mb);
00368 assert(h_ht1);
00369 assert(h_ht2);
00370 assert(h_ev);
00371 f.Close();
00372 h_minvMB = 0;
00373 h_minvHT1 = 0;
00374 h_minvHT2 = 0;
00375
00376 TFile file_nbar(settings.input_nbareff_file.Data(),"READ");
00377 const TH1F *h_effnbarMB = (const TH1F*)file_nbar.Get("h_effMB");
00378 const TH1F *h_effnbarHT1 = (const TH1F*)file_nbar.Get("h_effHT1");
00379 const TH1F *h_effnbarHT2 = (const TH1F*)file_nbar.Get("h_effHT2");
00380 assert(h_effnbarMB);
00381 assert(h_effnbarHT1);
00382 assert(h_effnbarHT2);
00383 TH1F *h_nbarEffMB=new TH1F(*h_effnbarMB);
00384 TH1F *h_nbarEffHT1=new TH1F(*h_effnbarHT1);
00385 TH1F *h_nbarEffHT2=new TH1F(*h_effnbarHT2);
00386 assert(h_nbarEffMB);
00387 assert(h_nbarEffHT1);
00388 assert(h_nbarEffHT2);
00389 h_nbarEffMB->Sumw2();
00390 h_nbarEffHT1->Sumw2();
00391 h_nbarEffHT2->Sumw2();
00392 file_nbar.Close();
00393 h_effnbarMB = 0;
00394 h_effnbarHT1 = 0;
00395 h_effnbarHT2 = 0;
00396
00397
00398 int trigger=0;
00399 Float_t numberOfMB=0;
00400 Float_t numberOfHT1=0;
00401 Float_t numberOfHT2=0;
00402 for(Int_t i=1;i<=h_ev->GetNbinsX();i++)
00403 {
00404 trigger=(Int_t)h_ev->GetBinCenter(i);
00405 if(trigger&1) numberOfMB+=(Int_t)h_ev->GetBinContent(i);
00406 if(trigger&2) numberOfHT1+=(Int_t)h_ev->GetBinContent(i);
00407 if(trigger&4) numberOfHT2+=(Int_t)h_ev->GetBinContent(i);
00408 }
00409
00410 cout<<"number of mb: "<<numberOfMB<<endl;
00411 if (settings.name == "dAu") {
00412 numberOfMB/=0.93;
00413 cout<<"nmb after 63% vertex eff.: "<<numberOfMB<<endl;
00414 }
00415
00416 Float_t psMB=1;
00417 Float_t psHT1=1;
00418 Float_t psHT2=1;
00419
00420 if (settings.name == "pp05") {
00421 psMB=32742;
00422 psHT1=3.89;
00423 psHT2=1.;
00424 }
00425 if (settings.name == "dAu") {
00426 psMB=383.;
00427 psHT1=9.65;
00428 psHT2=1.;
00429 }
00430
00431
00432 float nMBwithHT=1;
00433 if (settings.name == "pp05") {
00434 nMBwithHT = 328871;
00435 }
00436 if (settings.name == "dAu") {
00437 nMBwithHT = numberOfMB;
00438 }
00439
00440 TFile g(settings.input_pioneff_file.Data(),"READ");
00441 const TH1F *h_effpionMB = (const TH1F*)g.Get("h_effMB");
00442 const TH1F *h_effpionHT1 = (const TH1F*)g.Get("h_effHT1");
00443 const TH1F *h_effpionHT2 = (const TH1F*)g.Get("h_effHT2");
00444 assert(h_effpionMB);
00445 assert(h_effpionHT1);
00446 assert(h_effpionHT2);
00447 TH1F *h_emb=new TH1F(*h_effpionMB);
00448 TH1F *h_eht1=new TH1F(*h_effpionHT1);
00449 TH1F *h_eht2=new TH1F(*h_effpionHT2);
00450 assert(h_emb);
00451 assert(h_eht1);
00452 assert(h_eht2);
00453 g.Close();
00454 h_effpionMB = 0;
00455 h_effpionHT1 = 0;
00456 h_effpionHT2 = 0;
00457
00458
00459 TFile binf(settings.input_binwidth_file,"READ");
00460 const TH1F *h4mb = (const TH1F*)binf.Get("h4mb");
00461 const TH1F *h4ht1 = (const TH1F*)binf.Get("h4ht1");
00462 const TH1F *h4ht2 = (const TH1F*)binf.Get("h4ht2");
00463 assert(h4mb);
00464 assert(h4ht1);
00465 assert(h4ht2);
00466 TH1F *h_binmb=new TH1F(*h4mb);
00467 TH1F *h_binht1=new TH1F(*h4ht1);
00468 TH1F *h_binht2=new TH1F(*h4ht2);
00469 assert(h_binmb);
00470 assert(h_binht1);
00471 assert(h_binht2);
00472 binf.Close();
00473 h4mb = 0;
00474 h4ht1 = 0;
00475 h4ht2 = 0;
00476
00477 h_binmb->Sumw2();
00478 h_binht1->Sumw2();
00479 h_binht2->Sumw2();
00480 for(Int_t i=1;i<=h_binmb->GetNbinsX();i++) h_binmb->SetBinError(i,0);
00481 for(Int_t i=1;i<=h_binht1->GetNbinsX();i++) h_binht1->SetBinError(i,0);
00482 for(Int_t i=1;i<=h_binht2->GetNbinsX();i++) h_binht2->SetBinError(i,0);
00483
00484
00485
00486 TF1 *pion_cpv_corrMB=0;
00487 TF1 *pion_cpv_corrHT1=0;
00488 TF1 *pion_cpv_corrHT2=0;
00489
00490 TF1 *gamma_cpv_corrMB=0;
00491 TF1 *gamma_cpv_corrHT1=0;
00492 TF1 *gamma_cpv_corrHT2=0;
00493 TF1 *gamma_cont_corrMB=0;
00494 TF1 *gamma_cont_corrHT1=0;
00495 TF1 *gamma_cont_corrHT2=0;
00496
00497
00498 TF1 *pion_conv_corrMB=0;
00499 TF1 *pion_conv_corrHT1=0;
00500 TF1 *pion_conv_corrHT2=0;
00501
00502 TF1 *gamma_conv_corrMB=0;
00503 TF1 *gamma_conv_corrHT1=0;
00504 TF1 *gamma_conv_corrHT2=0;
00505 if (settings.name == "pp05") {
00506 pion_cpv_corrMB = new TF1("pion_cpv_corrMB","1./(1.-0.01*(0.3+0.0*x))",0.,15.);
00507 pion_cpv_corrHT1 = new TF1("pion_cpv_corrHT1","1./(1.-0.01*(-0.1+0.16*x))",0.,15.);
00508 pion_cpv_corrHT2 = new TF1("pion_cpv_corrHT2","1./(1.-0.01*(-0.2+0.18*x))",0.,15.);
00509 gamma_cpv_corrMB = new TF1("gamma_cpv_corrMB","1./(1.-0.01*(2.8+0.0*x))",0.,15.);
00510 gamma_cpv_corrHT1 = new TF1("gamma_cpv_corrHT1","1./(1.-0.01*(0.2+1.1*x))",0.,15.);
00511 gamma_cpv_corrHT2 = new TF1("gamma_cpv_corrHT2","1./(1.-0.01*(0.4+1.1*x))",0.,15.);
00512 gamma_cont_corrMB = new TF1("gamma_cont_corrMB","0.985",0.,15.);
00513 gamma_cont_corrHT1 = new TF1("gamma_cont_corrHT1","0.98",0.,15.);
00514 gamma_cont_corrHT2 = new TF1("gamma_cont_corrHT2","0.96",0.,15.);
00515 pion_conv_corrMB = new TF1("pion_conv_corrMB","1.1",0.,15.);
00516 pion_conv_corrHT1 = new TF1("pion_conv_corrHT1","1.1",0.,15.);
00517 pion_conv_corrHT2 = new TF1("pion_conv_corrHT2","1.1",0.,15.);
00518 gamma_conv_corrMB = new TF1("gamma_conv_corrMB","1.05",0.,15.);
00519 gamma_conv_corrHT1 = new TF1("gamma_conv_corrHT1","1.05",0.,15.);
00520 gamma_conv_corrHT2 = new TF1("gamma_conv_corrHT2","1.05",0.,15.);
00521 }
00522 if (settings.name == "dAu") {
00523 pion_cpv_corrMB = new TF1("pion_cpv_corrMB","1./(1.-0.01*(0.4+0.05*x))",0.,15.);
00524 pion_cpv_corrHT1 = new TF1("pion_cpv_corrHT1","1./(1.-0.01*(0.4+0.07*x))",0.,15.);
00525 pion_cpv_corrHT2 = new TF1("pion_cpv_corrHT2","1./(1.-0.01*(0.5+0.06*x))",0.,15.);
00526 gamma_cpv_corrMB = new TF1("gamma_cpv_corrMB","1./(1.-0.01*(4.1+0.4*x))",0.,15.);
00527 gamma_cpv_corrHT1 = new TF1("gamma_cpv_corrHT1","1./(1.-0.01*(3.4+0.5*x))",0.,15.);
00528 gamma_cpv_corrHT2 = new TF1("gamma_cpv_corrHT2","1./(1.-0.01*(4.9+0.3*x))",0.,15.);
00529 gamma_cont_corrMB = new TF1("gamma_cont_corrMB","0.98",0.,15.);
00530 gamma_cont_corrHT1 = new TF1("gamma_cont_corrHT1","0.98",0.,15.);
00531 gamma_cont_corrHT2 = new TF1("gamma_cont_corrHT2","0.965",0.,15.);
00532 pion_conv_corrMB = new TF1("pion_conv_corrMB","1.15",0.,15.);
00533 pion_conv_corrHT1 = new TF1("pion_conv_corrHT1","1.15",0.,15.);
00534 pion_conv_corrHT2 = new TF1("pion_conv_corrHT2","1.15",0.,15.);
00535 gamma_conv_corrMB = new TF1("gamma_conv_corrMB","1.08",0.,15.);
00536 gamma_conv_corrHT1 = new TF1("gamma_conv_corrHT1","1.08",0.,15.);
00537 gamma_conv_corrHT2 = new TF1("gamma_conv_corrHT2","1.08",0.,15.);
00538 }
00539 assert(pion_cpv_corrMB);
00540 assert(pion_cpv_corrHT1);
00541 assert(pion_cpv_corrHT2);
00542 assert(gamma_cpv_corrMB);
00543 assert(gamma_cpv_corrHT1);
00544 assert(gamma_cpv_corrHT2);
00545 assert(gamma_cont_corrMB);
00546 assert(gamma_cont_corrHT1);
00547 assert(gamma_cont_corrHT2);
00548 assert(pion_conv_corrMB);
00549 assert(pion_conv_corrHT1);
00550 assert(pion_conv_corrHT2);
00551 assert(gamma_conv_corrMB);
00552 assert(gamma_conv_corrHT1);
00553 assert(gamma_conv_corrHT2);
00554
00555
00556 Pi0Analysis *pi0=new Pi0Analysis(settings.output_invmassplots_file.Data(),settings.output_invmassplotseta_file.Data(),settings.name.Data());
00557 assert(pi0);
00558 pi0->init(settings.output_pionhistograms_file.Data());
00559 TH1F *pionYieldMB=new TH1F(*pi0->getYield(h_mb,"mb"));
00560 TH1F *pionYieldHT1=new TH1F(*pi0->getYield(h_ht1,"ht1"));
00561 TH1F *pionYieldHT2=new TH1F(*pi0->getYield(h_ht2,"ht2"));
00562 assert(pionYieldMB);
00563 assert(pionYieldHT1);
00564 assert(pionYieldHT2);
00565 pi0->storeCanvases(settings.output_pioncanvases_file.Data());
00566
00567
00568 cout<<"***************************************"<<endl;
00569 cout<<"got yield, dividing by rapidity bite!!!"<<endl;
00570 float dy_gamma=cuts->rapidityMaxCUT - cuts->rapidityMinCUT;
00571 float dy_pion=cuts->rapPionMaxCUT - cuts->rapPionMinCUT;
00572 cout<<"***************************************"<<endl;
00573 cout<<endl;
00574 cout<<" pion bite is "<<dy_pion<<endl;
00575 cout<<" gamma bite is "<<dy_gamma<<endl;
00576 cout<<endl;
00577 cout<<"***************************************"<<endl;
00578
00579 pionYieldMB->Scale(1./dy_pion);
00580 pionYieldHT1->Scale(1./dy_pion);
00581 pionYieldHT2->Scale(1./dy_pion);
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611 pionYieldMB->SetMarkerStyle(8);
00612 pionYieldMB->SetMarkerSize(1.0);
00613 pionYieldHT1->SetMarkerStyle(8);
00614 pionYieldHT1->SetMarkerSize(1.0);
00615 pionYieldHT1->SetMarkerColor(4);
00616 pionYieldHT2->SetMarkerStyle(8);
00617 pionYieldHT2->SetMarkerSize(1.0);
00618 pionYieldHT2->SetMarkerColor(2);
00619
00620 TF1 *scale=new TF1("scale","x",0.,15.);
00621 assert(scale);
00622
00623 pionYieldMB->SetNameTitle("pionYieldMB","corrected yield MB");
00624 pionYieldMB->Divide(h_emb);
00625 pionYieldMB->Scale(psMB/(psMB*numberOfMB*2.*TMath::Pi()));
00626 pionYieldMB->Divide(scale);
00627 pionYieldMB->Multiply(h_binmb);
00628 pionYieldMB->Multiply(pion_cpv_corrMB);
00629 pionYieldMB->Multiply(pion_conv_corrMB);
00630
00631 pionYieldHT1->SetNameTitle("pionYieldHT1","corrected yield HT1");
00632 pionYieldHT1->Divide(h_eht1);
00633 pionYieldHT1->Scale(psHT1/(psMB*nMBwithHT*2.*TMath::Pi()));
00634 pionYieldHT1->Divide(scale);
00635 pionYieldHT1->Multiply(h_binht1);
00636 pionYieldHT1->Multiply(pion_cpv_corrHT1);
00637 pionYieldHT1->Multiply(pion_conv_corrHT1);
00638
00639 pionYieldHT2->SetNameTitle("pionYieldHT2","corrected yield HT2");
00640 pionYieldHT2->Divide(h_eht2);
00641 pionYieldHT2->Scale(psHT2/(psMB*nMBwithHT*2.*TMath::Pi()));
00642 pionYieldHT2->Divide(scale);
00643 pionYieldHT2->Multiply(h_binht2);
00644 pionYieldHT2->Multiply(pion_cpv_corrHT2);
00645 pionYieldHT2->Multiply(pion_conv_corrHT2);
00646
00647
00648 TH1F *pionYieldMBratio=new TH1F(*pionYieldMB);
00649 TH1F *pionYieldHT1ratio=new TH1F(*pionYieldHT1);
00650 TH1F *pionYieldHT2ratio=new TH1F(*pionYieldHT2);
00651 assert(pionYieldMBratio);
00652 assert(pionYieldHT1ratio);
00653 assert(pionYieldHT2ratio);
00654
00655 TH1F *pionXsMB=new TH1F(*pionYieldMB);
00656 TH1F *pionXsHT1=new TH1F(*pionYieldHT1);
00657 TH1F *pionXsHT2=new TH1F(*pionYieldHT2);
00658 assert(pionXsMB);
00659 assert(pionXsHT1);
00660 assert(pionXsHT2);
00661 if (settings.name=="pp05") {
00662 pionXsMB->Scale((26.1/0.85));
00663 pionXsHT1->Scale((26.1/0.85));
00664 pionXsHT2->Scale((26.1/0.85));
00665 }
00666 if (settings.name=="dAu") {
00667 pionXsMB->Scale(2.21e3);
00668 pionXsHT1->Scale(2.21e3);
00669 pionXsHT2->Scale(2.21e3);
00670 }
00671
00672
00673 TH1F *pionXsMBnoErr=new TH1F(*pionXsMB);
00674 TH1F *pionXsHT1noErr=new TH1F(*pionXsHT1);
00675 TH1F *pionXsHT2noErr=new TH1F(*pionXsHT2);
00676 assert(pionXsMBnoErr);
00677 assert(pionXsHT1noErr);
00678 assert(pionXsHT2noErr);
00679 for(int i=1;i<=pionXsMBnoErr->GetNbinsX();i++){
00680 pionXsMBnoErr->SetBinError(i,0.);
00681 }
00682 for(int i=1;i<=pionXsHT1noErr->GetNbinsX();i++){
00683 pionXsHT1noErr->SetBinError(i,0.);
00684 }
00685 for(int i=1;i<=pionXsHT2noErr->GetNbinsX();i++){
00686 pionXsHT2noErr->SetBinError(i,0.);
00687 }
00688
00689 TGraphErrors *g_pionXsMB=new TGraphErrors(pionXsMB);
00690 assert(g_pionXsMB);
00691 g_pionXsMB->SetName("g_pionXsMB");
00692 removeThesePoints(g_pionXsMB,1);
00693
00694 TGraphErrors *g_pionXsHT1=new TGraphErrors(pionXsHT1);
00695 assert(g_pionXsHT1);
00696 g_pionXsHT1->SetName("g_pionXsHT1");
00697 removeThesePoints(g_pionXsHT1,2);
00698
00699 TGraphErrors *g_pionXsHT2=new TGraphErrors(pionXsHT2);
00700 assert(g_pionXsHT2);
00701 g_pionXsHT2->SetName("g_pionXsHT2");
00702 removeThesePoints(g_pionXsHT2,3);
00703
00704 if(1){
00705 cout<<endl<<"xsec: x y ex ey"<<endl;
00706 cout<<"minbias"<<endl;
00707 printPoints(g_pionXsMB);
00708 cout<<endl<<"hightower-1"<<endl;
00709 printPoints(g_pionXsHT1);
00710 cout<<endl<<"hightower-2"<<endl;
00711 printPoints(g_pionXsHT2);
00712 cout<<endl;
00713 }
00714
00715
00716 TMultiGraph *m_pions_fit=new TMultiGraph();
00717 assert(m_pions_fit);
00718 m_pions_fit->SetName("m_pions_fit");
00719 m_pions_fit->SetMinimum(5.0e-9);
00720 m_pions_fit->SetMaximum(0.99);
00721
00722
00723 m_pions_fit->Add(g_pionXsMB);
00724 m_pions_fit->Add(g_pionXsHT1);
00725 m_pions_fit->Add(g_pionXsHT2);
00726
00727 TF1 *fitQCD=new TF1("fitQCD",sumpqcd,1.,15.,6);
00728 assert(fitQCD);
00729 if (settings.name=="pp05") fitQCD->SetParameters(600.,-8.2,4.,-8.5,2.,2.);
00730 if (settings.name=="dAu") fitQCD->SetParameters(100.,-9.,1.,-8.5,1.,6.);
00731 fitQCD->FixParameter(4,2.);
00732 cout << "111" << endl;
00733
00734 cout << "222" << endl;
00735
00736 bool inclPhenix=true;
00737 bool inclFrank=true;
00738 bool inclSasha=true;
00739 bool inclPqcd=true;
00740
00741 TCanvas *compare=new TCanvas("compare","compare;p_{T}:xsec (mb)",600,750);
00742 assert(compare);
00743 compare->Divide(1, 2, 0, 0);
00744
00745 compare->cd(1);
00746 gPad->SetLogy();
00747
00748 TMultiGraph *m_pions=new TMultiGraph();
00749 assert(m_pions);
00750 m_pions->SetName("m_pions");
00751
00752 if(inclPqcd){
00753 m_pions->Add(kkp,"c");
00754 m_pions->Add(kkp05,"c");
00755 m_pions->Add(kkp2,"c");
00756 }
00757 if(inclSasha){
00758 m_pions->Add(sasha_pp05);
00759 }
00760 if(inclFrank){
00761 m_pions->Add(frank_pp);
00762 }
00763 if(inclPhenix){
00764 m_pions->Add(phenix_pp);
00765 }
00766
00767
00768
00769
00770
00771
00772
00773
00774 m_pions->Add(g_pionXsMB);
00775 m_pions->Add(g_pionXsHT1);
00776 m_pions->Add(g_pionXsHT2);
00777
00778 m_pions->SetMinimum(1.0e-9);
00779 m_pions->SetMaximum(1.);
00780
00781 m_pions->Draw("ap");
00782
00783
00784
00785 m_pions->GetXaxis()->SetLabelSize(0.);
00786 m_pions->GetXaxis()->SetTitle("p_{T} [GeV/c]");
00787 m_pions->GetYaxis()->SetTitle("Ed^{3}#sigma/dp^{3} [mb GeV^{-2} c^{3}]");
00788
00789 TLegend *leg=new TLegend(.5,.5,.85,.85);
00790 assert(leg);
00791
00792 if(inclPhenix) leg->AddEntry(phenix_pp,"PHENIX p+p","p");
00793 if(inclFrank) leg->AddEntry(frank_pp,"STAR preliminary (upd.)","p");
00794 if(inclSasha) leg->AddEntry(sasha_pp05,"O.Grebenyuk p+p","p");
00795 leg->AddEntry(g_pionXsMB,"p+p minimum bias","p");
00796 leg->AddEntry(g_pionXsHT1,"hightower 1","p");
00797 leg->AddEntry(g_pionXsHT2,"hightower 2","p");
00798 if(inclPqcd){
00799 leg->AddEntry(kkp,"kkp + CTEQ6m, #mu=p_{T}","l");
00800 leg->AddEntry(kkp2,"#mu=2p_{T},p_{T}/2","l");
00801 leg->Draw("same");
00802 }
00803
00804 leg->SetFillColor(0);
00805 leg->Draw();
00806
00807
00808
00809
00810
00811
00812 compare->cd(2);
00813
00814 TGraphErrors *sasha_pp05_overPqcd=new TGraphErrors(*sasha_pp05);
00815 assert(sasha_pp05_overPqcd);
00816 divideGraphWithGraph(sasha_pp05_overPqcd,kkp);
00817 TGraphErrors *phenix_pp05_overPqcd=new TGraphErrors(*phenix_pp);
00818 assert(phenix_pp05_overPqcd);
00819 divideGraphWithGraph(phenix_pp05_overPqcd,kkp);
00820 TGraphErrors *g_pionXsMB_overPqcd=new TGraphErrors(*g_pionXsMB);
00821 assert(g_pionXsMB_overPqcd);
00822 divideGraphWithGraph(g_pionXsMB_overPqcd,kkp);
00823 TGraphErrors *g_pionXsHT1_overPqcd=new TGraphErrors(*g_pionXsHT1);
00824 assert(g_pionXsHT1_overPqcd);
00825 divideGraphWithGraph(g_pionXsHT1_overPqcd,kkp);
00826 TGraphErrors *g_pionXsHT2_overPqcd=new TGraphErrors(*g_pionXsHT2);
00827 assert(g_pionXsHT2_overPqcd);
00828 divideGraphWithGraph(g_pionXsHT2_overPqcd,kkp);
00829 TGraphErrors *frank_pp05_overPqcd=new TGraphErrors(*frank_pp);
00830 assert(frank_pp05_overPqcd);
00831 divideGraphWithGraph(frank_pp05_overPqcd,kkp);
00832
00833 TGraphErrors *kkp05_ratio=new TGraphErrors(*kkp05);
00834 assert(kkp05_ratio);
00835 divideGraphWithGraph(kkp05_ratio,kkp);
00836 TGraphErrors *kkp_ratio=new TGraphErrors(*kkp);
00837 assert(kkp_ratio);
00838 divideGraphWithGraph(kkp_ratio,kkp);
00839 TGraphErrors *kkp2_ratio=new TGraphErrors(*kkp2);
00840 assert(kkp2_ratio);
00841 divideGraphWithGraph(kkp2_ratio,kkp);
00842
00843
00844 TGraphErrors *g_pionXsMB_sys=new TGraphErrors(*g_pionXsMB_overPqcd);
00845 assert(g_pionXsMB_sys);
00846 set_sys_pp_pion(g_pionXsMB_sys);
00847 TGraphErrors *g_pionXsHT1_sys=new TGraphErrors(*g_pionXsHT1_overPqcd);
00848 assert(g_pionXsHT1_sys);
00849 set_sys_pp_pion(g_pionXsHT1_sys);
00850 TGraphErrors *g_pionXsHT2_sys=new TGraphErrors(*g_pionXsHT2_overPqcd);
00851 assert(g_pionXsHT2_sys);
00852 set_sys_pp_pion(g_pionXsHT2_sys);
00853
00854 TMultiGraph *m_pions_over_pqcd=new TMultiGraph();
00855 assert(m_pions_over_pqcd);
00856
00857 m_pions_over_pqcd->SetMinimum(0.0);
00858 m_pions_over_pqcd->SetMaximum(2.5);
00859
00860 if(inclPqcd){
00861 m_pions_over_pqcd->Add(kkp05_ratio,"c");
00862 m_pions_over_pqcd->Add(kkp_ratio,"c");
00863 m_pions_over_pqcd->Add(kkp2_ratio,"c");
00864 }
00865 m_pions_over_pqcd->Add(g_pionXsMB_overPqcd);
00866 m_pions_over_pqcd->Add(g_pionXsHT1_overPqcd);
00867 m_pions_over_pqcd->Add(g_pionXsHT2_overPqcd);
00868
00869
00870
00871 if(inclPhenix) m_pions_over_pqcd->Add(phenix_pp05_overPqcd);
00872 if(inclFrank) m_pions_over_pqcd->Add(frank_pp05_overPqcd);
00873 if(inclSasha) m_pions_over_pqcd->Add(sasha_pp05_overPqcd);
00874
00875 m_pions_over_pqcd->Draw("ap");
00876 m_pions_over_pqcd->GetXaxis()->SetTitle("p_{T} [GeV/c]");
00877 m_pions_over_pqcd->GetYaxis()->SetTitle("Data / theory");
00878
00879 compare->SaveAs(settings.output_pionxsec_file.Data());
00880
00881 TMultiGraph *m_pions_over_fit=new TMultiGraph();
00882 assert(m_pions_over_fit);
00883 m_pions_over_fit->SetMinimum(0.01);
00884 m_pions_over_fit->SetMaximum(1.99);
00885
00886 TGraphErrors *g_pionXsMB_overFit=new TGraphErrors(*g_pionXsMB);
00887 assert(g_pionXsMB_overFit);
00888 divideGraphWithFunction(g_pionXsMB_overFit,fitQCD);
00889 TGraphErrors *g_pionXsHT1_overFit=new TGraphErrors(*g_pionXsHT1);
00890 assert(g_pionXsHT1_overFit);
00891 divideGraphWithFunction(g_pionXsHT1_overFit,fitQCD);
00892 TGraphErrors *g_pionXsHT2_overFit=new TGraphErrors(*g_pionXsHT2);
00893 assert(g_pionXsHT2_overFit);
00894 divideGraphWithFunction(g_pionXsHT2_overFit,fitQCD);
00895
00896 m_pions_over_fit->Add(g_pionXsMB_overFit);
00897 m_pions_over_fit->Add(g_pionXsHT1_overFit);
00898 m_pions_over_fit->Add(g_pionXsHT2_overFit);
00899
00900 m_pions_over_fit->Draw("ap");
00901
00902 compare->SaveAs(settings.output_pionxsecoverfit_file.Data());
00903
00904
00905
00906 TCanvas *compare2=new TCanvas("compare2","compare2;p_{T};yield divided by fit",600,300);
00907 assert(compare2);
00908 compare2->cd();
00909
00910
00911 TGraphErrors *frank_pp2=new TGraphErrors(*frank_pp);
00912 assert(frank_pp2);
00913 divideGraphWithFunction(frank_pp2,fitQCD);
00914 TGraphErrors *sasha_pp2=new TGraphErrors(*sasha_pp05);
00915 assert(sasha_pp2);
00916 divideGraphWithFunction(sasha_pp2,fitQCD);
00917 TGraphErrors *phenix_pp2=new TGraphErrors(*phenix_pp);
00918 assert(phenix_pp2);
00919 divideGraphWithFunction(phenix_pp2,fitQCD);
00920 TGraphErrors *g_pionXsMBcopy=new TGraphErrors(*g_pionXsMB);
00921 assert(g_pionXsMBcopy);
00922 divideGraphWithFunction(g_pionXsMBcopy,fitQCD);
00923 TGraphErrors *g_pionXsHT1copy=new TGraphErrors(*g_pionXsHT1);
00924 assert(g_pionXsHT1copy);
00925 divideGraphWithFunction(g_pionXsHT1copy,fitQCD);
00926 TGraphErrors *g_pionXsHT2copy=new TGraphErrors(*g_pionXsHT2);
00927 assert(g_pionXsHT2copy);
00928 divideGraphWithFunction(g_pionXsHT2copy,fitQCD);
00929
00930 inclPhenix=true;
00931 inclFrank=false;
00932 inclSasha=false;
00933 inclPqcd=false;
00934
00935 TMultiGraph *m_pions2=new TMultiGraph();
00936 assert(m_pions2);
00937 m_pions2->SetName("m_pions2");
00938 if(inclSasha) m_pions2->Add(sasha_pp2);
00939 if(inclFrank) m_pions2->Add(frank_pp2);
00940 if(inclPhenix) m_pions2->Add(phenix_pp2);
00941
00942 m_pions2->Add(g_pionXsMBcopy);
00943 m_pions2->Add(g_pionXsHT1copy);
00944 m_pions2->Add(g_pionXsHT2copy);
00945
00946 m_pions2->SetMinimum(0.000001);
00947 m_pions2->SetMaximum(3.);
00948 m_pions2->Draw("ap");
00949
00950 TLegend *legg=new TLegend(.25,.55,.65,.85);
00951 assert(legg);
00952 legg->AddEntry(g_pionXsMBcopy,"minimum bias","p");
00953 legg->AddEntry(g_pionXsHT1copy,"hightower 1","p");
00954 legg->AddEntry(g_pionXsHT2copy,"hightower 2","p");
00955 if(inclFrank) legg->AddEntry(frank_pp2,"Frank's p+p (upd.)","p");
00956 if(inclSasha) legg->AddEntry(sasha_pp2,"Sasha's p+p","p");
00957 if(inclPhenix) legg->AddEntry(phenix_pp,"PHENIX p+p","p");
00958 legg->Draw("same");
00959 legg->SetFillColor(0);
00960
00961 compare2->cd(0);
00962 compare2->SaveAs(settings.output_pionxsecratio_file.Data());
00963
00964
00965
00966
00967
00968
00969
00970 TFile gg(settings.input_pioneff_file.Data(),"READ");
00971 const TH1F *h_effDaughtersMB = (const TH1F*)gg.Get("h_effDaughtersMB");
00972 const TH1F *h_effDaughtersHT1 = (const TH1F*)gg.Get("h_effDaughtersHT1");
00973 const TH1F *h_effDaughtersHT2 = (const TH1F*)gg.Get("h_effDaughtersHT2");
00974 assert(h_effDaughtersMB);
00975 assert(h_effDaughtersHT1);
00976 assert(h_effDaughtersHT2);
00977 TH1F *effGammaMB=new TH1F(*h_effDaughtersMB);
00978 TH1F *effGammaHT1=new TH1F(*h_effDaughtersHT1);
00979 TH1F *effGammaHT2=new TH1F(*h_effDaughtersHT2);
00980 assert(effGammaMB);
00981 assert(effGammaHT1);
00982 assert(effGammaHT2);
00983 gg.Close();
00984 h_effDaughtersMB = 0;
00985 h_effDaughtersHT1 = 0;
00986 h_effDaughtersHT2 = 0;
00987
00988
00989 TFile gg_single(settings.input_gammaeff_file.Data(),"READ");
00990 const TH1F *h_effgammaMB = (const TH1F*)gg_single.Get("h_effMB");
00991 const TH1F *h_effgammaHT1 = (const TH1F*)gg_single.Get("h_effHT1");
00992 const TH1F *h_effgammaHT2 = (const TH1F*)gg_single.Get("h_effHT2");
00993 assert(h_effgammaMB);
00994 assert(h_effgammaHT1);
00995 assert(h_effgammaHT2);
00996 TH1F *effGammaSingleMB=new TH1F(*h_effgammaMB);
00997 TH1F *effGammaSingleHT1=new TH1F(*h_effgammaHT1);
00998 TH1F *effGammaSingleHT2=new TH1F(*h_effgammaHT2);
00999 assert(effGammaSingleMB);
01000 assert(effGammaSingleHT1);
01001 assert(effGammaSingleHT2);
01002 gg_single.Close();
01003 h_effgammaMB = 0;
01004 h_effgammaHT1 = 0;
01005 h_effgammaHT2 = 0;
01006
01007
01008 TFile ff(settings.input_pion_file.Data(),"READ");
01009 const TH1F *h_gammaMB = (const TH1F*)ff.Get("h_gammaMB");
01010 const TH1F *h_gammaHT1 = (const TH1F*)ff.Get("h_gammaHT1");
01011 const TH1F *h_gammaHT2 = (const TH1F*)ff.Get("h_gammaHT2");
01012 assert(h_gammaMB);
01013 assert(h_gammaHT1);
01014 assert(h_gammaHT2);
01015 TH1F *gammaYieldMB=new TH1F(*h_gammaMB);
01016 TH1F *gammaYieldHT1=new TH1F(*h_gammaHT1);
01017 TH1F *gammaYieldHT2=new TH1F(*h_gammaHT2);
01018 assert(gammaYieldMB);
01019 assert(gammaYieldHT1);
01020 assert(gammaYieldHT2);
01021 ff.Close();
01022 h_gammaMB = 0;
01023 h_gammaHT1 = 0;
01024 h_gammaHT2 = 0;
01025
01026
01027 gammaYieldMB->Scale(1./dy_gamma);
01028 gammaYieldHT1->Scale(1./dy_gamma);
01029 gammaYieldHT2->Scale(1./dy_gamma);
01030
01031
01032 for(Int_t i=1;i<=gammaYieldMB->GetNbinsX();i++){
01033 gammaYieldMB->SetBinContent(i,gammaYieldMB->GetBinContent(i)/gammaYieldMB->GetXaxis()->GetBinWidth(i));
01034 gammaYieldMB->SetBinError(i,gammaYieldMB->GetBinError(i)/gammaYieldMB->GetXaxis()->GetBinWidth(i));
01035 }
01036 gammaYieldMB->Scale(psMB/(psMB*numberOfMB*2.*TMath::Pi()));
01037 gammaYieldMB->Divide(scale);
01038 gammaYieldMB->Multiply(h_binmb);
01039 gammaYieldMB->Divide(effGammaMB);
01040 gammaYieldMB->Multiply(gamma_cpv_corrMB);
01041 gammaYieldMB->Multiply(gamma_cont_corrMB);
01042 gammaYieldMB->Multiply(gamma_conv_corrMB);
01043
01044 gammaYieldMB->SetMarkerStyle(8);
01045 gammaYieldMB->SetMarkerSize(1.);
01046 TGraphErrors *g_inclPhotonsMB=new TGraphErrors(gammaYieldMB);
01047 assert(g_inclPhotonsMB);
01048
01049 for(Int_t i=1;i<=gammaYieldHT1->GetNbinsX();i++){
01050 gammaYieldHT1->SetBinContent(i,gammaYieldHT1->GetBinContent(i)/gammaYieldHT1->GetXaxis()->GetBinWidth(i));
01051 gammaYieldHT1->SetBinError(i,gammaYieldHT1->GetBinError(i)/gammaYieldHT1->GetXaxis()->GetBinWidth(i));
01052 }
01053 gammaYieldHT1->Scale(psHT1/(psMB*nMBwithHT*2.*TMath::Pi()));
01054 gammaYieldHT1->Divide(scale);
01055 gammaYieldHT1->Multiply(h_binht1);
01056 gammaYieldHT1->Divide(effGammaHT1);
01057 gammaYieldHT1->Multiply(gamma_cpv_corrHT1);
01058 gammaYieldHT1->Multiply(gamma_cont_corrHT1);
01059 gammaYieldHT1->Multiply(gamma_conv_corrHT1);
01060
01061 gammaYieldHT1->SetMarkerStyle(8);
01062 gammaYieldHT1->SetMarkerSize(1.);
01063 gammaYieldHT1->SetMarkerColor(4);
01064 TGraphErrors *g_inclPhotonsHT1=new TGraphErrors(gammaYieldHT1);
01065 assert(g_inclPhotonsHT1);
01066
01067 for(Int_t i=1;i<=gammaYieldHT2->GetNbinsX();i++){
01068 gammaYieldHT2->SetBinContent(i,gammaYieldHT2->GetBinContent(i)/gammaYieldHT2->GetXaxis()->GetBinWidth(i));
01069 gammaYieldHT2->SetBinError(i,gammaYieldHT2->GetBinError(i)/gammaYieldHT2->GetXaxis()->GetBinWidth(i));
01070 }
01071 gammaYieldHT2->Scale(psHT2/(psMB*nMBwithHT*2.*TMath::Pi()));
01072 gammaYieldHT2->Divide(scale);
01073 gammaYieldHT2->Multiply(h_binht2);
01074 gammaYieldHT2->Divide(effGammaHT2);
01075 gammaYieldHT2->Multiply(gamma_cpv_corrHT2);
01076 gammaYieldHT2->Multiply(gamma_cont_corrHT2);
01077 gammaYieldHT2->Multiply(gamma_conv_corrHT2);
01078
01079 gammaYieldHT2->SetMarkerStyle(8);
01080 gammaYieldHT2->SetMarkerSize(1.);
01081 gammaYieldHT2->SetMarkerColor(2);
01082 TGraphErrors *g_inclPhotonsHT2=new TGraphErrors(gammaYieldHT2);
01083 assert(g_inclPhotonsHT2);
01084
01085 removeThesePoints(g_inclPhotonsMB,1);
01086 removeThesePoints(g_inclPhotonsHT1,2);
01087 removeThesePoints(g_inclPhotonsHT2,2);
01088
01089 TMultiGraph *m_incl=new TMultiGraph();
01090 assert(m_incl);
01091 m_incl->SetName("m_incl");
01092 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}}");
01093 m_incl->Add(g_inclPhotonsMB);
01094 m_incl->Add(g_inclPhotonsHT1);
01095 m_incl->Add(g_inclPhotonsHT2);
01096
01097
01098
01099 m_incl->SetMinimum(1.e-11);
01100 m_incl->SetMaximum(10.);
01101 TCanvas *c_incl=new TCanvas("c_incl","c_incl",600,400);
01102 assert(c_incl);
01103 gPad->SetLogy();
01104 m_incl->Draw("ap");
01105 c_incl->SaveAs(settings.output_inclphotonyield_file.Data());
01106
01107
01108
01109
01110 TH1F *gammaYieldMBratio=new TH1F(*gammaYieldMB);
01111 TH1F *gammaYieldHT1ratio=new TH1F(*gammaYieldHT1);
01112 TH1F *gammaYieldHT2ratio=new TH1F(*gammaYieldHT2);
01113 assert(gammaYieldMBratio);
01114 assert(gammaYieldHT1ratio);
01115 assert(gammaYieldHT2ratio);
01116 gammaYieldMBratio->SetName("gammaYieldMBratio");
01117 gammaYieldHT1ratio->SetName("gammaYieldHT1ratio");
01118 gammaYieldHT2ratio->SetName("gammaYieldHT2ratio");
01119
01120 gammaYieldMBratio->Divide(pionYieldMBratio);
01121 gammaYieldHT1ratio->Divide(pionYieldHT1ratio);
01122 gammaYieldHT2ratio->Divide(pionYieldHT2ratio);
01123
01124
01125
01126 getRatio(gammaYieldMBratio,effGammaMB,effGammaSingleMB,fit_piondecay);
01127 getRatio(gammaYieldHT1ratio,effGammaHT1,effGammaSingleHT1,fit_piondecay);
01128 getRatio(gammaYieldHT2ratio,effGammaHT2,effGammaSingleHT2,fit_piondecay);
01129
01130 TH1F *gammaYieldMBratio_incl=new TH1F(*gammaYieldMBratio);
01131 TH1F *gammaYieldHT1ratio_incl=new TH1F(*gammaYieldHT1ratio);
01132 TH1F *gammaYieldHT2ratio_incl=new TH1F(*gammaYieldHT2ratio);
01133 assert(gammaYieldMBratio_incl);
01134 assert(gammaYieldHT1ratio_incl);
01135 assert(gammaYieldHT2ratio_incl);
01136
01137 TH1F *gammaYieldMBratioNoErr=new TH1F(*gammaYieldMBratio);
01138 TH1F *gammaYieldHT1ratioNoErr=new TH1F(*gammaYieldHT1ratio);
01139 TH1F *gammaYieldHT2ratioNoErr=new TH1F(*gammaYieldHT2ratio);
01140 assert(gammaYieldMBratioNoErr);
01141 assert(gammaYieldHT1ratioNoErr);
01142 assert(gammaYieldHT2ratioNoErr);
01143 for(int i=1;i<=gammaYieldMBratioNoErr->GetNbinsX();i++)gammaYieldMBratioNoErr->SetBinError(i,0.);
01144 for(int i=1;i<=gammaYieldHT1ratioNoErr->GetNbinsX();i++)gammaYieldHT1ratioNoErr->SetBinError(i,0.);
01145 for(int i=1;i<=gammaYieldHT2ratioNoErr->GetNbinsX();i++)gammaYieldHT2ratioNoErr->SetBinError(i,0.);
01146
01147 TGraphErrors *g_ratioMB=new TGraphErrors(gammaYieldMBratio);
01148 TGraphErrors *g_ratioHT1=new TGraphErrors(gammaYieldHT1ratio);
01149 TGraphErrors *g_ratioHT2=new TGraphErrors(gammaYieldHT2ratio);
01150 assert(g_ratioMB);
01151 assert(g_ratioHT1);
01152 assert(g_ratioHT2);
01153 g_ratioMB->SetName("g_ratioMB");
01154 g_ratioHT1->SetName("g_ratioHT1");
01155 g_ratioHT2->SetName("g_ratioHT2");
01156
01157
01158 removeThesePoints(g_ratioMB,1);
01159 removeThesePoints(g_ratioHT1,2);
01160 removeThesePoints(g_ratioHT2,3);
01161
01162
01163 TCanvas *c_ratio=new TCanvas("c_ratio","c_ratio",400,200);
01164 assert(c_ratio);
01165
01166 TMultiGraph *m_ratio=new TMultiGraph("m_ratio","p+p 2005;p_{T};#gamma/#pi^{0}");
01167 assert(m_ratio);
01168 m_ratio->Add(g_ratioMB);
01169 m_ratio->Add(g_ratioHT1);
01170 m_ratio->Add(g_ratioHT2);
01171
01172 m_ratio->Draw("ap");
01173 m_ratio->SetMinimum(.001);
01174 m_ratio->SetMaximum(1.5);
01175
01176 TLegend *leg3=new TLegend(.35,.65,.65,.85);
01177 assert(leg3);
01178 leg3->AddEntry(g_ratioMB,"minimum bias","p");
01179 leg3->AddEntry(g_ratioHT1,"hightower 1","p");
01180 leg3->AddEntry(g_ratioHT2,"hightower 2","p");
01181 leg3->AddEntry(fit_decay,"decay background (total)","l");
01182 leg3->AddEntry(fit_piondecay,"decay background (#pi^{0})","l");
01183 leg3->SetFillColor(0);
01184 leg3->Draw("same");
01185
01186 fit_decay->SetLineColor(13);
01187 fit_decay->SetLineWidth(1);
01188 fit_decay->SetLineColor(1);
01189 fit_decay->Draw("same");
01190 fit_piondecay->SetLineColor(13);
01191 fit_piondecay->SetLineWidth(1);
01192 fit_piondecay->SetLineStyle(2);
01193 fit_piondecay->SetLineColor(1);
01194 fit_piondecay->Draw("same");
01195
01196 c_ratio->SaveAs(settings.output_gammaoverpion_file.Data());
01197
01198
01199 gammaYieldMBratio_incl->Multiply(pionYieldMBratio);
01200 gammaYieldHT1ratio_incl->Multiply(pionYieldHT1ratio);
01201 gammaYieldHT2ratio_incl->Multiply(pionYieldHT2ratio);
01202
01203
01204
01205
01206 TGraphErrors *g_incl_corrMB=new TGraphErrors(gammaYieldMBratio_incl);
01207 TGraphErrors *g_incl_corrHT1=new TGraphErrors(gammaYieldHT1ratio_incl);
01208 TGraphErrors *g_incl_corrHT2=new TGraphErrors(gammaYieldHT2ratio_incl);
01209 assert(g_incl_corrMB);
01210 assert(g_incl_corrHT1);
01211 assert(g_incl_corrHT2);
01212
01213 TCanvas *c_incl_corr=new TCanvas("c_incl_corr","c_incl_corr",400,300);
01214 assert(c_incl_corr);
01215 gPad->SetLogy();
01216 TMultiGraph *m_incl_corr=new TMultiGraph();
01217 assert(m_incl_corr);
01218 m_incl_corr->Add(g_incl_corrMB);
01219 m_incl_corr->Add(g_incl_corrHT1);
01220 m_incl_corr->Add(g_incl_corrHT2);
01221
01222 m_incl_corr->SetMinimum(1.e-11);
01223 m_incl_corr->SetMaximum(1.);
01224
01225 m_incl_corr->Draw("apX");
01226 c_incl_corr->SaveAs(settings.output_inclphotonyieldcorr_file.Data());
01227
01228 TCanvas *c_doubleratio=new TCanvas("c_doubleratio","c_doubleratio",400,300);
01229 assert(c_doubleratio);
01230 gStyle->SetOptStat(0);
01231 c_doubleratio->cd(1);
01232
01233 TH1F *gammaYieldMBdoubleratio=new TH1F(*gammaYieldMBratio);
01234 TH1F *gammaYieldHT1doubleratio=new TH1F(*gammaYieldHT1ratio);
01235 TH1F *gammaYieldHT2doubleratio=new TH1F(*gammaYieldHT2ratio);
01236 assert(gammaYieldMBdoubleratio);
01237 assert(gammaYieldHT1doubleratio);
01238 assert(gammaYieldHT2doubleratio);
01239
01240 gammaYieldMBdoubleratio->Divide(fit_decay);
01241 gammaYieldHT1doubleratio->Divide(fit_decay);
01242 gammaYieldHT2doubleratio->Divide(fit_decay);
01243
01244 TGraphErrors *g_doubleRatioMB=new TGraphErrors(gammaYieldMBdoubleratio);
01245 assert(g_doubleRatioMB);
01246 g_doubleRatioMB->SetName("g_doubleRatioMB");
01247 g_doubleRatioMB->SetMarkerStyle(8);
01248 TGraphErrors *g_doubleRatioHT1=new TGraphErrors(gammaYieldHT1doubleratio);
01249 assert(g_doubleRatioHT1);
01250 g_doubleRatioHT1->SetName("g_doubleRatioHT1");
01251 g_doubleRatioHT1->SetMarkerStyle(8);
01252 TGraphErrors *g_doubleRatioHT2=new TGraphErrors(gammaYieldHT2doubleratio);
01253 assert(g_doubleRatioHT2);
01254 g_doubleRatioHT2->SetName("g_doubleRatioHT2");
01255 g_doubleRatioHT2->SetMarkerStyle(8);
01256
01257 removeThesePoints(g_doubleRatioMB,1);
01258 removeThesePoints(g_doubleRatioHT1,2);
01259 removeThesePoints(g_doubleRatioHT2,3);
01260
01261 TMultiGraph *m_doubleratio=new TMultiGraph();
01262 assert(m_doubleratio);
01263 m_doubleratio->SetName("m_doubleratio");
01264 m_doubleratio->SetMinimum(.5);
01265 m_doubleratio->SetMaximum(2.75);
01266
01267 cout<<endl;
01268 g_doubleRatioHT1->Print();
01269 cout<<endl<<endl;
01270 g_doubleRatioHT2->Print();
01271 cout<<endl;
01272
01273
01274 m_doubleratio->Add(g_doubleRatioHT1,"p");
01275 m_doubleratio->Add(g_doubleRatioHT2,"p");
01276
01277 g_photonpqcd->SetLineWidth(2);
01278 g_photonpqcd->SetLineColor(2);
01279 g_photonpqcd05->SetLineWidth(2);
01280 g_photonpqcd05->SetLineColor(2);
01281 g_photonpqcd05->SetLineStyle(2);
01282 g_photonpqcd2->SetLineWidth(2);
01283 g_photonpqcd2->SetLineColor(2);
01284 g_photonpqcd2->SetLineStyle(2);
01285
01286 m_doubleratio->Add(g_photonpqcd,"c");
01287 m_doubleratio->Add(g_photonpqcd05,"c");
01288 m_doubleratio->Add(g_photonpqcd2,"c");
01289
01290
01291 TF1 *fitGamma2=new TF1("fitGamma2","1.+[0]*TMath::Power(x,[1])",2.,15.);
01292 assert(fitGamma2);
01293 g_photonpqcd->Fit(fitGamma2,"R0");
01294
01295 m_doubleratio->Draw("a");
01296
01297 m_doubleratio->GetXaxis()->SetTitle("p_{T} (GeV/c)");
01298 m_doubleratio->GetYaxis()->SetTitle("1 + #gamma_{dir}/#gamma_{incl}");
01299 m_doubleratio->GetXaxis()->SetRangeUser(2.,16.);
01300
01301
01302 TLegend *leg5=new TLegend(.15,.6,.6,.8);
01303 assert(leg5);
01304 leg5->AddEntry(g_doubleRatioHT1,"hightower-1","p");
01305 leg5->AddEntry(g_doubleRatioHT2,"hightower-2","p");
01306 leg5->AddEntry(g_photonpqcd,"NLO (CTEQ6+KKP) #mu=p_{T}","l");
01307 leg5->AddEntry(g_photonpqcd05,"#mu=2p_{T}, #mu=p_{T}/2","l");
01308 leg5->SetFillColor(0);
01309 leg5->Draw("same");
01310
01311 c_doubleratio->cd(0);
01312 c_doubleratio->SaveAs(settings.output_gammadoubleratio_file.Data());
01313
01314 {
01315
01316 TH1F *h_nbarContMB = new TH1F(*h_nbarEffMB);
01317 TH1F *h_nbarContHT1 = new TH1F(*h_nbarEffHT1);
01318 TH1F *h_nbarContHT2 = new TH1F(*h_nbarEffHT2);
01319 assert(h_nbarContMB);
01320 assert(h_nbarContHT1);
01321 assert(h_nbarContHT2);
01322 h_nbarContMB->Multiply(f_nbar);
01323 h_nbarContHT1->Multiply(f_nbar);
01324 h_nbarContHT2->Multiply(f_nbar);
01325
01326 h_nbarContMB->Divide(gammaYieldMB);
01327 h_nbarContHT1->Divide(gammaYieldHT1);
01328 h_nbarContHT2->Divide(gammaYieldHT2);
01329
01330 TH1F *effGammaMBnoerr = new TH1F(*effGammaMB);
01331 TH1F *effGammaHT1noerr = new TH1F(*effGammaHT1);
01332 TH1F *effGammaHT2noerr = new TH1F(*effGammaHT2);
01333 assert(effGammaMBnoerr);
01334 assert(effGammaHT1noerr);
01335 assert(effGammaHT2noerr);
01336 for (Int_t i = 1;i <= effGammaMBnoerr->GetXaxis()->GetNbins();i++) effGammaMBnoerr->SetBinError(i, 0);
01337 for (Int_t i = 1;i <= effGammaHT1noerr->GetXaxis()->GetNbins();i++) effGammaHT1noerr->SetBinError(i, 0);
01338 for (Int_t i = 1;i <= effGammaHT2noerr->GetXaxis()->GetNbins();i++) effGammaHT2noerr->SetBinError(i, 0);
01339 h_nbarContMB->Divide(effGammaMBnoerr);
01340 h_nbarContHT1->Divide(effGammaHT1noerr);
01341 h_nbarContHT2->Divide(effGammaHT2noerr);
01342
01343 TCanvas *c_nbar_cont = new TCanvas("c_nbar_cont", "Antineutron contamination");
01344 assert(c_nbar_cont);
01345 TH1F *h_nbar_cont = new TH1F("h_nbar_cont", "Antineutron contamination C_{0};p_{T} [GeV/c];C_{0}", 1000, 0, 10);
01346 assert(h_nbar_cont);
01347 h_nbar_cont->Draw();
01348 h_nbar_cont->GetYaxis()->SetRangeUser(0, 1.5);
01349
01350 h_nbarContMB->SetMarkerStyle(kOpenCircle);
01351 h_nbarContMB->SetMarkerSize(1.0);
01352 h_nbarContMB->SetMarkerColor(kBlack);
01353 h_nbarContMB->SetLineStyle(kSolid);
01354 h_nbarContMB->SetLineWidth(1);
01355 h_nbarContMB->SetLineColor(kBlack);
01356
01357 h_nbarContHT1->SetMarkerStyle(kOpenSquare);
01358 h_nbarContHT1->SetMarkerSize(1.0);
01359 h_nbarContHT1->SetMarkerColor(kBlue);
01360 h_nbarContHT1->SetLineStyle(kSolid);
01361 h_nbarContHT1->SetLineWidth(1);
01362 h_nbarContHT1->SetLineColor(kBlue);
01363
01364 h_nbarContHT2->SetMarkerStyle(kOpenTriangleUp);
01365 h_nbarContHT2->SetMarkerSize(1.0);
01366 h_nbarContHT2->SetMarkerColor(kRed);
01367 h_nbarContHT2->SetLineStyle(kSolid);
01368 h_nbarContHT2->SetLineWidth(1);
01369 h_nbarContHT2->SetLineColor(kRed);
01370
01371 h_nbarContMB->Draw("SAME");
01372 h_nbarContHT1->Draw("SAME");
01373 h_nbarContHT2->Draw("SAME");
01374
01375
01376 Float_t maxCont = h_nbarContMB->GetMaximum();
01377 TH1F *h_nbarContHT2_extreme = new TH1F(*h_nbarContHT2);
01378 assert(h_nbarContHT2_extreme);
01379 h_nbarContHT2_extreme->Scale(1.0/maxCont);
01380 h_nbarContHT2_extreme->SetLineStyle(kDashed);
01381 h_nbarContHT2_extreme->SetLineWidth(1);
01382 h_nbarContHT2_extreme->SetLineColor(kRed);
01383 h_nbarContHT2_extreme->Draw("SAME HIST C");
01384
01385
01386
01387
01388 TH1F *h_corr = new TH1F(*gammaYieldMBdoubleratio);
01389 assert(h_corr);
01390 for (Int_t i = 1;i < h_corr->GetXaxis()->GetNbins();i++) {
01391 h_corr->SetBinContent(i, 1.0);
01392 h_corr->SetBinError(i, 0.0);
01393 }
01394 h_corr->Divide(gammaYieldMBdoubleratio);
01395
01396 for (Int_t i = 1;i < h_corr->GetXaxis()->GetNbins();i++) {
01397 h_corr->SetBinContent(i, 1.0 - h_corr->GetBinContent(i));
01398 }
01399 h_corr->Divide(h_nbarContMB);
01400 h_corr->GetXaxis()->SetRangeUser(1.0, 4.0);
01401 TF1 *f_corr = new TF1("f_corr", "[0]");
01402 f_corr->SetRange(1.0, 4.0);
01403 h_corr->Fit(f_corr, "RQN");
01404 Float_t corrCoef = f_corr->GetParameter(0);
01405 TH1F *h_nbarContHT2_corr = new TH1F(*h_nbarContHT2);
01406 assert(h_nbarContHT2_corr);
01407 h_nbarContHT2_corr->Scale(corrCoef);
01408 h_nbarContHT2_corr->SetLineStyle(kSolid);
01409 h_nbarContHT2_corr->SetLineWidth(1);
01410 h_nbarContHT2_corr->SetLineColor(17);
01411 h_nbarContHT2_corr->SetFillColor(17);
01412 h_nbarContHT2_corr->SetFillStyle(1001);
01413 h_nbarContHT2_corr->Draw("SAME HIST AC");
01414
01415 TLegend *leg = new TLegend(0.5, 0.5, 0.8, 0.8);
01416 assert(leg);
01417 leg->AddEntry(h_nbarContMB, "MinBias", "P");
01418 leg->AddEntry(h_nbarContHT1, "HighTower-1", "P");
01419 leg->AddEntry(h_nbarContHT2, "HighTower-2", "P");
01420 leg->AddEntry(h_nbarContHT2_extreme, "extreme case", "L");
01421 leg->AddEntry(h_nbarContHT2_corr, "corrected", "F");
01422 leg->Draw();
01423
01424 h_nbar_cont->Draw("SAME");
01425 c_nbar_cont->SaveAs(settings.output_nbarcont_file.Data());
01426 }
01427
01428 TCanvas *c_dirphoton=new TCanvas("c_dirphoton","c_dirphoton",400,300);
01429 assert(c_dirphoton);
01430 gStyle->SetOptStat(0);
01431 c_dirphoton->cd(1);
01432
01433 TH1F *dirphotonYieldMB=new TH1F(*gammaYieldMBdoubleratio);
01434 TH1F *dirphotonYieldHT1=new TH1F(*gammaYieldHT1doubleratio);
01435 TH1F *dirphotonYieldHT2=new TH1F(*gammaYieldHT2doubleratio);
01436 assert(dirphotonYieldMB);
01437 assert(dirphotonYieldHT1);
01438 assert(dirphotonYieldHT2);
01439
01440 TH1F *dirphotonYieldMBnoErr=new TH1F(*dirphotonYieldMB);
01441 assert(dirphotonYieldMBnoErr);
01442 for(int i=1;i<=dirphotonYieldMBnoErr->GetNbinsX();i++){
01443 dirphotonYieldMBnoErr->SetBinError(i,0.);
01444 }
01445 TH1F *dirphotonYieldHT1noErr=new TH1F(*dirphotonYieldHT1);
01446 assert(dirphotonYieldHT1noErr);
01447 for(int i=1;i<=dirphotonYieldHT1noErr->GetNbinsX();i++){
01448 dirphotonYieldHT1noErr->SetBinError(i,0.);
01449 }
01450 TH1F *dirphotonYieldHT2noErr=new TH1F(*dirphotonYieldHT2);
01451 assert(dirphotonYieldHT2noErr);
01452 for(int i=1;i<=dirphotonYieldHT1noErr->GetNbinsX();i++){
01453 dirphotonYieldHT2noErr->SetBinError(i,0.);
01454 }
01455
01456 TF1 *f_unity=new TF1("f_unity","1.",0.,15.);
01457 assert(f_unity);
01458 dirphotonYieldMB->Add(f_unity,-1.);
01459 dirphotonYieldHT1->Add(f_unity,-1.);
01460 dirphotonYieldHT2->Add(f_unity,-1.);
01461
01462
01463 dirphotonYieldMB->Divide(dirphotonYieldMBnoErr);
01464 dirphotonYieldHT1->Divide(dirphotonYieldHT1noErr);
01465 dirphotonYieldHT2->Divide(dirphotonYieldHT2noErr);
01466 dirphotonYieldMB->Multiply(gammaYieldMBratioNoErr);
01467 dirphotonYieldHT1->Multiply(gammaYieldHT1ratioNoErr);
01468 dirphotonYieldHT2->Multiply(gammaYieldHT2ratioNoErr);
01469 dirphotonYieldMB->Multiply(pionXsMBnoErr);
01470 dirphotonYieldHT1->Multiply(pionXsHT1noErr);
01471 dirphotonYieldHT2->Multiply(pionXsHT2noErr);
01472
01473
01474 TGraphErrors *g_dirphotonMB=new TGraphErrors(dirphotonYieldMB);
01475 assert(g_dirphotonMB);
01476 g_dirphotonMB->SetName("g_dirphotonMB");
01477 g_dirphotonMB->SetMarkerStyle(8);
01478 TGraphErrors *g_dirphotonHT1=new TGraphErrors(dirphotonYieldHT1);
01479 assert(g_dirphotonHT1);
01480 g_dirphotonHT1->SetName("g_dirphotonHT1");
01481 g_dirphotonHT1->SetMarkerStyle(8);
01482 TGraphErrors *g_dirphotonHT2=new TGraphErrors(dirphotonYieldHT2);
01483 assert(g_dirphotonHT2);
01484 g_dirphotonHT2->SetName("g_dirphotonHT2");
01485 g_dirphotonHT2->SetMarkerStyle(8);
01486
01487
01488 removeThesePoints(g_dirphotonMB,1);
01489 removeThesePoints(g_dirphotonHT1,2);
01490 removeThesePoints(g_dirphotonHT2,3);
01491
01492 gPad->SetLogy();
01493
01494 TMultiGraph *m_dirphoton=new TMultiGraph();
01495 assert(m_dirphoton);
01496 m_dirphoton->SetName("m_dirphoton");
01497 m_dirphoton->SetMinimum(1.0e-11);
01498 m_dirphoton->SetMaximum(0.1);
01499
01500 m_dirphoton->Add(g_dirgamma,"c");
01501 m_dirphoton->Add(g_dirgamma05,"c");
01502 m_dirphoton->Add(g_dirgamma2,"c");
01503
01504 cout<<"direct photons:"<<endl;
01505 g_dirphotonHT1->Print();
01506 cout<<endl;
01507 g_dirphotonHT2->Print();
01508 cout<<endl;
01509
01510 m_dirphoton->Add(g_dirphotonHT1,"p");
01511 m_dirphoton->Add(g_dirphotonHT2,"p");
01512
01513 m_dirphoton->Draw("a");
01514
01515 m_dirphoton->GetXaxis()->SetTitle("p_{T} (GeV/c)");
01516 m_dirphoton->GetYaxis()->SetTitle("1 - R^{-1}");
01517 m_dirphoton->GetXaxis()->SetRangeUser(2.,16.);
01518
01519
01520 TLegend *leght=new TLegend(.15,.6,.6,.8);
01521 assert(leght);
01522 leght->AddEntry(g_dirphotonHT1,"hightower-1","p");
01523 leght->AddEntry(g_dirphotonHT2,"hightower-2","p");
01524 leght->SetFillColor(0);
01525 leght->Draw("same");
01526
01527 c_dirphoton->cd(0);
01528 c_dirphoton->SaveAs(settings.output_gammadirphoton_file.Data());
01529 }
01530
01531
