00001 void getSpectrumDAU(const char *pid)
00002 {
00003 gStyle->SetOptStat(0);
00004 gStyle->SetOptFit(0);
00005
00006 gStyle->SetErrorX(0);
00007
00008 Bool_t subtractNEUTRONS=kTRUE;
00009
00010
00011
00012
00013
00014 Double_t B[]={0.3,0.23,0.072,0.061};
00015
00016 Double_t T[]={0.205,0.215,0.179,0.173};
00017
00018 Double_t n[]={11.00,12.55,10.87,10.49};
00019
00020 TF1 *f_antineutron=new TF1("f_antineutron","[0]/pow((1.+(sqrt(x*x+1.) - 1.)/([1]*[2])),[2])",0.,15.);
00021 f_antineutron->SetParameters(B[1],T[1],n[1]);
00022
00023
00024 ifstream pQCDphotons("./datapoints/pQCD_Werner/rhic_cteq6_gamma_inv_sc1.dat");
00025 Float_t ppx[100];
00026 Float_t ppy[100];
00027 Int_t iii=0;
00028 cout<<"pqcd photons:"<<endl;
00029 while(iii<28){
00030 if(!pQCDphotons.good()) break;
00031 Float_t dummy=0.;
00032 pQCDphotons>>ppx[iii]>>dummy>>dummy>>ppy[iii];
00033 ppy[iii]*=1.e-09*7.5/42.;
00034 iii++;
00035 }
00036 TGraph *g_dirgamma=new TGraph(iii,ppx,ppy);
00037
00038 ifstream pQCDphotons05("./datapoints/pQCD_Werner/rhic_cteq6_gamma_inv_sc05.dat");
00039 Float_t ppx05[100];
00040 Float_t ppy05[100];
00041 Int_t iii05=0;
00042 while(iii05<28){
00043 if(!pQCDphotons05.good()) break;
00044 Float_t dummy=0.;
00045 pQCDphotons05>>ppx05[iii05]>>dummy>>dummy>>ppy05[iii05];
00046 ppy05[iii05]*=1.e-09*7.5/42.;
00047 iii05++;
00048 }
00049 TGraph *g_dirgamma05=new TGraph(iii05,ppx05,ppy05);
00050
00051 ifstream pQCDphotons2("./datapoints/pQCD_Werner/rhic_cteq6_gamma_inv_sc2.dat");
00052 Float_t ppx2[100];
00053 Float_t ppy2[100];
00054 Int_t iii2=0;
00055 while(iii2<28){
00056 if(!pQCDphotons2.good()) break;
00057 Float_t dummy=0.;
00058 pQCDphotons2>>ppx2[iii2]>>dummy>>dummy>>ppy2[iii2];
00059 ppy2[iii2]*=1.e-09*7.5/42.;
00060 iii2++;
00061 }
00062 TGraph *g_dirgamma2=new TGraph(iii2,ppx2,ppy2);
00063
00064
00065 ifstream phenix("./datapoints/phenix.dat");
00066 Float_t phex[100];
00067 Float_t phey[100];
00068 Float_t ephex[100];
00069 Float_t ephey[100];
00070 Int_t iphe=0;
00071 while(iphe<18){
00072 if(!phenix.good()) break;
00073 phenix>>phex[iphe]>>phey[iphe]>>ephey[iphe];
00074 ephex[iphe]=0.;
00075 iphe++;
00076 }
00077 TGraphErrors *phenix_dau=new TGraphErrors(iphe,phex,phey,ephex,ephey);
00078 phenix_dau->SetMarkerStyle(24);
00079 phenix_dau->SetName("phenix_dau");
00080
00081
00082 ifstream pQCDpions("./datapoints/pQCD_Werner/klaus_pi0inv_200_kkp_1.dat");
00083 Float_t pionx[100];
00084 Float_t piony[100];
00085 Int_t ipion=0;
00086 while(ipion<28){
00087 if(!pQCDpions.good()) break;
00088 pQCDpions>>pionx[ipion]>>piony[ipion];
00089 ipion++;
00090 }
00091 TGraphErrors *kkp=new TGraphErrors(ipion,pionx,piony);
00092 kkp->SetLineColor(54);
00093 kkp->SetName("kkp");
00094
00095
00096 ifstream pQCDpions05("./datapoints/pQCD_Werner/klaus_pi0inv_200_kkp_05.dat");
00097 Float_t pionx05[100];
00098 Float_t piony05[100];
00099 Int_t ipion05=0;
00100 while(ipion05<28){
00101 if(!pQCDpions05.good()) break;
00102 pQCDpions05>>pionx05[ipion05]>>piony05[ipion05];
00103 ipion05++;
00104 }
00105 TGraphErrors *kkp05=new TGraphErrors(ipion05,pionx05,piony05);
00106 kkp05->SetLineStyle(2);
00107 kkp05->SetLineColor(54);
00108 kkp05->SetName("kkp05");
00109
00110
00111 ifstream pQCDpions2("./datapoints/pQCD_Werner/klaus_pi0inv_200_kkp_2.dat");
00112 Float_t pionx2[100];
00113 Float_t piony2[100];
00114 Int_t ipion2=0;
00115 while(ipion2<28){
00116 if(!pQCDpions2.good()) break;
00117 pQCDpions2>>pionx2[ipion2]>>piony2[ipion2];
00118 ipion2++;
00119 }
00120 TGraphErrors *kkp2=new TGraphErrors(ipion2,pionx2,piony2);
00121 kkp2->SetLineStyle(2);
00122 kkp2->SetLineColor(54);
00123 kkp2->SetName("kkp2");
00124
00125 TFile *f_decaybg=new TFile("~/MyDecay/gammaDecayDAUSum.root","OPEN");
00126 TH1F *h_decaybg=(TH1F*)f_decaybg->Get("gamma");
00127 TH1F *h_decaypion=(TH1F*)f_decaybg->Get("gamma_pion");
00128 TF1 *fit_decay=new TF1("fit_decay","[0]/pow(x,[1])+[2]",.3,15.);
00129 TF1 *fit_piondecay=new TF1(*fit_decay);
00130 fit_decay->SetParameters(1.,1.,.5);
00131 fit_piondecay->SetParameters(1.,1.,.5);
00132 h_decaybg->Fit(fit_decay,"R0");
00133 h_decaypion->Fit(fit_piondecay,"R0Q");
00134
00135 TCanvas *c_test=new TCanvas();
00136 h_decaybg->Draw();
00137 fit_decay->Draw("same");
00138 c_test->SaveAs("test.eps");
00139
00140
00141 for(Int_t i=0;i<iii;i++){
00142 ppy[i]=ppy[i]/piony[i];
00143 ppy[i]=ppy[i]/fit_decay->Eval(ppx[i]);
00144 ppy[i]+=1.;
00145 ppy05[i]=ppy05[i]/piony05[i];
00146 ppy05[i]=ppy05[i]/fit_decay->Eval(ppx05[i]);
00147 ppy05[i]+=1.;
00148 ppy2[i]=ppy2[i]/piony2[i];
00149 ppy2[i]=ppy2[i]/fit_decay->Eval(ppx2[i]);
00150 ppy2[i]+=1.;
00151 }
00152 TGraphErrors *g_photonpqcd=new TGraphErrors(iii,ppx,ppy);
00153 g_photonpqcd->SetName("g_photonpqcd");
00154 TGraphErrors *g_photonpqcd05=new TGraphErrors(iii05,ppx05,ppy05);
00155 g_photonpqcd05->SetName("g_photonpqcd05");
00156 TGraphErrors *g_photonpqcd2=new TGraphErrors(iii2,ppx2,ppy2);
00157 g_photonpqcd2->SetName("g_photonpqcd2");
00158
00159
00160 TString dir("/star/u/russcher/gamma/analysis/output/dAu/");
00161 dir.Append(pid);
00162 dir.Append("/");
00163 dir.Append(pid);
00164 TString psout=dir;
00165 TString psout_eta=dir;
00166 TString eFile=dir;
00167 TString eFileGamma=dir;
00168 TString pi0File=dir;
00169 TString nbarFile=dir;
00170
00171 eFile.Append("pion_eff.root");
00172 eFileGamma.Append("gamma_eff.root");
00173 pi0File.Append("pi0_dAu.root");
00174 nbarFile.Append("antineutron_eff.root");
00175 psout_eta.Append("/dev/null");
00176 psout.Append("invmassplots.ps");
00177
00178
00179 gSystem->Load("$HOME/MyEvent/MyEvent");
00180 gSystem->Load("$HOME/gamma/analysis/lib/AnaCuts");
00181 gSystem->Load("$HOME/gamma/analysis/lib/EventMixer");
00182 gSystem->Load("$HOME/gamma/analysis/lib/Pi0Analysis");
00183
00184
00185 AnaCuts *cuts=new AnaCuts("dAu");
00186
00187
00188 TFile f(pi0File.Data(),"OPEN");
00189 TH2F *h_mb=new TH2F(*h_minvMB);
00190 TH2F *h_ht1=new TH2F(*h_minvHT1);
00191 TH2F *h_ht2=new TH2F(*h_minvHT2);
00192 TH1F *h_ev=new TH1F(*h_events);
00193
00194
00195 if(0){
00196 TFile *file_nbar=new TFile(nbarFile,"OPEN");
00197 TH1F *h_nbarEffMB=new TH1F(*h_effMB);
00198 h_nbarEffMB->Sumw2();
00199 TH1F *h_nbarEffHT1=new TH1F(*h_effHT1);
00200 h_nbarEffHT1->Sumw2();
00201 TH1F *h_nbarEffHT2=new TH1F(*h_effHT2);
00202 h_nbarEffHT2->Sumw2();
00203 }
00204
00205 int trigger=0;
00206 Int_t numberOfMB=0;
00207 Int_t numberOfHT1=0;
00208 Int_t numberOfHT2=0;
00209 for(Int_t i=1;i<=h_ev->GetNbinsX();i++)
00210 {
00211 trigger=(Int_t)h_ev->GetBinCenter(i);
00212 if(trigger&1) numberOfMB+=(Int_t)h_ev->GetBinContent(i);
00213 if(trigger&2) numberOfHT1+=(Int_t)h_ev->GetBinContent(i);
00214 if(trigger&4) numberOfHT2+=(Int_t)h_ev->GetBinContent(i);
00215 }
00216
00217 cout<<"number of mb: "<<numberOfMB<<endl;
00218 numberOfMB/=0.93;
00219 cout<<"nmb after 93% vertex eff.: "<<numberOfMB<<endl;
00220
00221 Float_t psMB=383.;
00222 Float_t psHT1=9.65;
00223 Float_t psHT2=1.;
00224
00225
00226 TFile g(eFile.Data(),"OPEN");
00227 TH1F *h_emb=new TH1F(*h_effMB);
00228 TH1F *h_eht1=new TH1F(*h_effHT1);
00229 TH1F *h_eht2=new TH1F(*h_effHT2);
00230
00231
00232
00233 TF1 *pion_cpv_corrMB=new TF1("pion_cpv_corrMB","1./(1.-0.01*(0.4+0.05*x))",0.,15.);
00234 TF1 *pion_cpv_corrHT1=new TF1("pion_cpv_corrHT1","1./(1.-0.01*(0.4+0.07*x))",0.,15.);
00235 TF1 *pion_cpv_corrHT2=new TF1("pion_cpv_corrHT2","1./(1.-0.01*(0.5+0.06*x))",0.,15.);
00236
00237 TF1 *gamma_cpv_corrMB=new TF1("gamma_cpv_corrMB","1./(1.-0.01*(4.1+0.4*x))",0.,15.);
00238 TF1 *gamma_cpv_corrHT1=new TF1("gamma_cpv_corrHT1","1./(1.-0.01*(3.4+0.5*x))",0.,15.);
00239 TF1 *gamma_cpv_corrHT2=new TF1("gamma_cpv_corrHT2","1./(1.-0.01*(4.9+0.3*x))",0.,15.);
00240 TF1 *gamma_cont_corrMB=new TF1("gamma_cont_corrMB","0.98",0.,15.);
00241 TF1 *gamma_cont_corrHT1=new TF1("gamma_cont_corrHT1","0.98",0.,15.);
00242 TF1 *gamma_cont_corrHT2=new TF1("gamma_cont_corrHT2","0.965",0.,15.);
00243
00244
00245
00246 TF1 *pion_conv_corrMB=new TF1("pion_conv_corrMB","1.15",0.,15.);
00247 TF1 *pion_conv_corrHT1=new TF1("pion_conv_corrHT1","1.15",0.,15.);
00248 TF1 *pion_conv_corrHT2=new TF1("pion_conv_corrHT2","1.15",0.,15.);
00249
00250 TF1 *gamma_conv_corrMB=new TF1("gamma_conv_corrMB","1.08",0.,15.);
00251 TF1 *gamma_conv_corrHT1=new TF1("gamma_conv_corrHT1","1.08",0.,15.);
00252 TF1 *gamma_conv_corrHT2=new TF1("gamma_conv_corrHT2","1.08",0.,15.);
00253
00254
00255 TFile binf("~/BinWidth/bincorrectionsDAU.root","OPEN");
00256 TH1F *h_binmb=new TH1F(*h4mb);
00257 TH1F *h_binht1=new TH1F(*h4ht1);
00258 TH1F *h_binht2=new TH1F(*h4ht2);
00259 h_binmb->Sumw2();
00260 h_binht1->Sumw2();
00261 h_binht2->Sumw2();
00262 for(Int_t i=1;i<=h_binmb->GetNbinsX();i++) h_binmb->SetBinError(i,0);
00263 for(Int_t i=1;i<=h_binht1->GetNbinsX();i++) h_binht1->SetBinError(i,0);
00264 for(Int_t i=1;i<=h_binht2->GetNbinsX();i++) h_binht2->SetBinError(i,0);
00265
00266
00267 Pi0Analysis *pi0=new Pi0Analysis(psout.Data(),psout_eta.Data(),"dAu");
00268 pi0->init("/dev/null");
00269 TH1F *pionYieldMB=new TH1F(*pi0->getYield(h_mb,"mb"));
00270 TH1F *pionYieldHT1=new TH1F(*pi0->getYield(h_ht1,"ht1"));
00271 TH1F *pionYieldHT2=new TH1F(*pi0->getYield(h_ht2,"ht2"));
00272 pi0->storeCanvases((dir+"canvases.root").Data());
00273
00274 cout<<"***************************************"<<endl;
00275 cout<<"got yield, dividing by rapidity bite!!!"<<endl;
00276 float dy_gamma=cuts->rapidityMaxCUT - cuts->rapidityMinCUT;
00277 float dy_pion=cuts->rapPionMaxCUT - cuts->rapPionMinCUT;
00278 cout<<"***************************************"<<endl;
00279 cout<<endl;
00280 cout<<" pion bite is "<<dy_pion<<endl;
00281 cout<<" gamma bite is "<<dy_gamma<<endl;
00282 cout<<endl;
00283 cout<<"***************************************"<<endl;
00284
00285 pionYieldMB->Scale(1./dy_pion);
00286 pionYieldHT1->Scale(1./dy_pion);
00287 pionYieldHT2->Scale(1./dy_pion);
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317 pionYieldMB->SetMarkerStyle(8);
00318 pionYieldMB->SetMarkerSize(1.0);
00319 pionYieldHT1->SetMarkerStyle(8);
00320 pionYieldHT1->SetMarkerSize(1.0);
00321 pionYieldHT1->SetMarkerColor(4);
00322 pionYieldHT2->SetMarkerStyle(8);
00323 pionYieldHT2->SetMarkerSize(1.0);
00324 pionYieldHT2->SetMarkerColor(2);
00325
00326 TF1 *scale=new TF1("scale","x",0.,15.);
00327
00328 pionYieldMB->SetNameTitle("pionYieldMB","corrected yield MB");
00329 pionYieldMB->Divide(h_emb);
00330 pionYieldMB->Scale(psMB/(psMB*numberOfMB*2.*TMath::Pi()));
00331 pionYieldMB->Divide(scale);
00332 pionYieldMB->Multiply(h_binmb);
00333 pionYieldMB->Multiply(pion_cpv_corrMB);
00334 pionYieldMB->Multiply(pion_conv_corrMB);
00335
00336 pionYieldHT1->SetNameTitle("pionYieldHT1","corrected yield HT1");
00337 pionYieldHT1->Divide(h_eht1);
00338 pionYieldHT1->Scale(psHT1/(psMB*numberOfMB*2.*TMath::Pi()));
00339 pionYieldHT1->Divide(scale);
00340 pionYieldHT1->Multiply(h_binht1);
00341 pionYieldHT1->Multiply(pion_cpv_corrHT1);
00342 pionYieldHT1->Multiply(pion_conv_corrHT1);
00343
00344 pionYieldHT2->SetNameTitle("pionYieldHT2","corrected yield HT2");
00345 pionYieldHT2->Divide(h_eht2);
00346 pionYieldHT2->Scale(psHT2/(psMB*numberOfMB*2.*TMath::Pi()));
00347 pionYieldHT2->Divide(scale);
00348 pionYieldHT2->Multiply(h_binht2);
00349 pionYieldHT2->Multiply(pion_cpv_corrHT2);
00350 pionYieldHT2->Multiply(pion_conv_corrHT2);
00351
00352
00353 TH1F *pionYieldMBratio=new TH1F(*pionYieldMB);
00354 TH1F *pionYieldHT1ratio=new TH1F(*pionYieldHT1);
00355 TH1F *pionYieldHT2ratio=new TH1F(*pionYieldHT2);
00356
00357 TH1F *pionXsMB=new TH1F(*pionYieldMB);
00358 TH1F *pionXsHT1=new TH1F(*pionYieldHT1);
00359 TH1F *pionXsHT2=new TH1F(*pionYieldHT2);
00360
00361 TH1F *pionXsMBnoErr=new TH1F(*pionXsMB);
00362 TH1F *pionXsHT1noErr=new TH1F(*pionXsHT1);
00363 TH1F *pionXsHT2noErr=new TH1F(*pionXsHT2);
00364 for(int i=1;i<=pionXsMBnoErr->GetNbinsX();i++){
00365 pionXsMBnoErr->SetBinError(i,0.);
00366 }
00367 for(int i=1;i<=pionXsHT1noErr->GetNbinsX();i++){
00368 pionXsHT1noErr->SetBinError(i,0.);
00369 }
00370 for(int i=1;i<=pionXsHT2noErr->GetNbinsX();i++){
00371 pionXsHT2noErr->SetBinError(i,0.);
00372 }
00373
00374 TGraphErrors *g_pionXsMB=new TGraphErrors(pionXsMB);
00375 g_pionXsMB->SetName("g_pionXsMB");
00376 removeThesePoints(g_pionXsMB,1);
00377
00378 TGraphErrors *g_pionXsHT1=new TGraphErrors(pionXsHT1);
00379 g_pionXsHT1->SetName("g_pionXsHT1");
00380 removeThesePoints(g_pionXsHT1,2);
00381
00382 TGraphErrors *g_pionXsHT2=new TGraphErrors(pionXsHT2);
00383 g_pionXsHT2->SetName("g_pionXsHT2");
00384 removeThesePoints(g_pionXsHT2,3);
00385
00386 if(0){
00387 cout<<endl<<"yield: x y ex ey"<<endl;
00388 cout<<"minbias"<<endl;
00389 printPoints(g_pionXsMB);
00390 cout<<endl<<"hightower-1"<<endl;
00391 printPoints(g_pionXsHT1);
00392 cout<<endl<<"hightower-2"<<endl;
00393 printPoints(g_pionXsHT2);
00394 cout<<endl;
00395 }
00396
00397
00398 TMultiGraph *m_pions_fit=new TMultiGraph();
00399 m_pions_fit->SetName("m_pions_fit");
00400 m_pions_fit->SetMinimum(1.0e-10);
00401 m_pions_fit->SetMaximum(.1);
00402
00403
00404 m_pions_fit->Add(g_pionXsMB);
00405 m_pions_fit->Add(g_pionXsHT1);
00406 m_pions_fit->Add(g_pionXsHT2);
00407
00408 TF1 *fitQCD=new TF1("fitQCD",sumpqcd,1.,15.,6);
00409 fitQCD->SetParameters(100.,-9.,1.,-8.5,1.,6.);
00410 fitQCD->FixParameter(4,2.);
00411 m_pions_fit->Fit(fitQCD,"R");
00412
00413 bool inclPhenix=true;
00414 bool inclSasha=false;
00415
00416 TCanvas *compare=new TCanvas("compare","compare;p_{T}:xsec (mb)",600,750);
00417 compare->cd();
00418
00419 TPad *padt=new TPad("padt","",0.0,0.3,1.,1.0);
00420 padt->SetBottomMargin(0.001);
00421 padt->SetLeftMargin(0.15);
00422 TPad *padb=new TPad("padb","",0.0,0.0,1.,0.3);
00423 padb->SetTopMargin(0.001);
00424 padb->SetBottomMargin(0.25);
00425 padb->SetLeftMargin(0.15);
00426
00427 padt->Draw();
00428 padt->cd();
00429 gPad->SetLogy();
00430
00431 TMultiGraph *m_pions=new TMultiGraph();
00432 m_pions->SetName("m_pions");
00433
00434 if(inclSasha){
00435
00436 }
00437 if(inclPhenix){
00438 m_pions->Add(phenix_dau);
00439 }
00440
00441 g_pionXsMB->Print();
00442 cout<<endl<<endl;
00443 g_pionXsHT1->Print();
00444 cout<<endl<<endl;
00445 g_pionXsHT2->Print();
00446 cout<<endl;
00447
00448
00449 m_pions->Add(g_pionXsMB);
00450 m_pions->Add(g_pionXsHT1);
00451 m_pions->Add(g_pionXsHT2);
00452
00453 m_pions->SetMinimum(5.0e-10);
00454 m_pions->SetMaximum(.099);
00455 m_pions->SetTitle(";p_{T} (GeV/c);invariant yield");
00456 m_pions->Draw("ap");
00457
00458
00459
00460 TLegend *leg=new TLegend(.5,.5,.85,.85);
00461
00462 if(inclPhenix) leg->AddEntry(phenix_dau,"PHENIX d+Au","p");
00463 if(inclSasha) leg->AddEntry(sasha_dau,"sasha's d+Au","p");
00464 leg->AddEntry(g_pionXsMB,"d+Au minimum bias","p");
00465 leg->AddEntry(g_pionXsHT1,"hightower 1","p");
00466 leg->AddEntry(g_pionXsHT2,"hightower 2","p");
00467
00468 leg->SetFillColor(0);
00469 leg->Draw();
00470
00471 compare->cd();
00472 padb->Draw();
00473 padb->cd();
00474
00475 TGraphErrors *phenix_dau2=new TGraphErrors(*phenix_dau);
00476 divideGraphWithFunction(phenix_dau2,fitQCD);
00477 TGraphErrors *g_pionXsMBcopy=new TGraphErrors(*g_pionXsMB);
00478 divideGraphWithFunction(g_pionXsMBcopy,fitQCD);
00479 TGraphErrors *g_pionXsHT1copy=new TGraphErrors(*g_pionXsHT1);
00480 divideGraphWithFunction(g_pionXsHT1copy,fitQCD);
00481 TGraphErrors *g_pionXsHT2copy=new TGraphErrors(*g_pionXsHT2);
00482 divideGraphWithFunction(g_pionXsHT2copy,fitQCD);
00483
00484 TMultiGraph *m_pions2=new TMultiGraph();
00485
00486 if(inclPhenix) m_pions2->Add(phenix_dau2);
00487 m_pions2->Add(g_pionXsMBcopy);
00488 m_pions2->Add(g_pionXsHT1copy);
00489 m_pions2->Add(g_pionXsHT2copy);
00490 m_pions2->SetMinimum(0.01);
00491 m_pions2->SetMaximum(1.99);
00492 m_pions2->Draw("ap");
00493
00494 compare->SaveAs((dir+"pionyield_dau.eps").Data());
00495 compare->SaveAs((dir+"pionyield_dau.root").Data());
00496
00497
00498
00499
00500
00501
00502
00503 TFile gg(eFile.Data(),"OPEN");
00504 TH1F *effGammaMB=new TH1F(*h_effDaughtersMB);
00505 TH1F *effGammaHT1=new TH1F(*h_effDaughtersHT1);
00506 TH1F *effGammaHT2=new TH1F(*h_effDaughtersHT2);
00507
00508
00509 TFile gg_single(eFileGamma.Data(),"OPEN");
00510 TH1F *effGammaSingleMB=new TH1F(*h_effMB);
00511 TH1F *effGammaSingleHT1=new TH1F(*h_effHT1);
00512 TH1F *effGammaSingleHT2=new TH1F(*h_effHT2);
00513
00514
00515 TFile ff(pi0File.Data(),"OPEN");
00516 TH1F *gammaYieldMB=new TH1F(*h_gammaMB);
00517 TH1F *gammaYieldHT1=new TH1F(*h_gammaHT1);
00518 TH1F *gammaYieldHT2=new TH1F(*h_gammaHT2);
00519
00520
00521 gammaYieldMB->Scale(1./dy_gamma);
00522 gammaYieldHT1->Scale(1./dy_gamma);
00523 gammaYieldHT2->Scale(1./dy_gamma);
00524
00525
00526 for(Int_t i=1;i<=gammaYieldMB->GetNbinsX();i++){
00527 gammaYieldMB->SetBinContent(i,gammaYieldMB->GetBinContent(i)/gammaYieldMB->GetXaxis()->GetBinWidth(i));
00528 gammaYieldMB->SetBinError(i,gammaYieldMB->GetBinError(i)/gammaYieldMB->GetXaxis()->GetBinWidth(i));
00529 }
00530 gammaYieldMB->Scale(psMB/(psMB*numberOfMB*2.*TMath::Pi()));
00531 gammaYieldMB->Divide(scale);
00532 gammaYieldMB->Multiply(h_binmb);
00533 gammaYieldMB->Divide(effGammaMB);
00534 gammaYieldMB->Multiply(gamma_cpv_corrMB);
00535 gammaYieldMB->Multiply(gamma_cont_corrMB);
00536 gammaYieldMB->Multiply(gamma_conv_corrMB);
00537
00538 gammaYieldMB->SetMarkerStyle(8);
00539 gammaYieldMB->SetMarkerSize(1.);
00540 TGraphErrors *g_inclPhotonsMB=new TGraphErrors(gammaYieldMB);
00541
00542
00543 for(Int_t i=1;i<=gammaYieldHT1->GetNbinsX();i++){
00544 gammaYieldHT1->SetBinContent(i,gammaYieldHT1->GetBinContent(i)/gammaYieldHT1->GetXaxis()->GetBinWidth(i));
00545 gammaYieldHT1->SetBinError(i,gammaYieldHT1->GetBinError(i)/gammaYieldHT1->GetXaxis()->GetBinWidth(i));
00546 }
00547 gammaYieldHT1->Scale(psHT1/(psMB*numberOfMB*2.*TMath::Pi()));
00548 gammaYieldHT1->Divide(scale);
00549 gammaYieldHT1->Multiply(h_binht1);
00550 gammaYieldHT1->Divide(effGammaHT1);
00551 gammaYieldHT1->Multiply(gamma_cpv_corrHT1);
00552 gammaYieldHT1->Multiply(gamma_cont_corrHT1);
00553 gammaYieldHT1->Multiply(gamma_conv_corrHT1);
00554
00555 gammaYieldHT1->SetMarkerStyle(8);
00556 gammaYieldHT1->SetMarkerSize(1.);
00557 gammaYieldHT1->SetMarkerColor(4);
00558 TGraphErrors *g_inclPhotonsHT1=new TGraphErrors(gammaYieldHT1);
00559
00560
00561 for(Int_t i=1;i<=gammaYieldHT2->GetNbinsX();i++){
00562 gammaYieldHT2->SetBinContent(i,gammaYieldHT2->GetBinContent(i)/gammaYieldHT2->GetXaxis()->GetBinWidth(i));
00563 gammaYieldHT2->SetBinError(i,gammaYieldHT2->GetBinError(i)/gammaYieldHT2->GetXaxis()->GetBinWidth(i));
00564 }
00565 gammaYieldHT2->Scale(psHT2/(psMB*numberOfMB*2.*TMath::Pi()));
00566 gammaYieldHT2->Divide(scale);
00567 gammaYieldHT2->Multiply(h_binht2);
00568 gammaYieldHT2->Divide(effGammaHT2);
00569 gammaYieldHT2->Multiply(gamma_cpv_corrHT2);
00570 gammaYieldHT2->Multiply(gamma_cont_corrHT2);
00571 gammaYieldHT2->Multiply(gamma_conv_corrHT2);
00572
00573 gammaYieldHT2->SetMarkerStyle(8);
00574 gammaYieldHT2->SetMarkerSize(1.);
00575 gammaYieldHT2->SetMarkerColor(2);
00576 TGraphErrors *g_inclPhotonsHT2=new TGraphErrors(gammaYieldHT2);
00577
00578
00579 removeThesePoints(g_inclPhotonsMB,1);
00580 removeThesePoints(g_inclPhotonsHT1,2);
00581 removeThesePoints(g_inclPhotonsHT2,3);
00582
00583
00584 TMultiGraph *m_incl=new TMultiGraph();
00585 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}}");
00586 m_incl->Add(g_inclPhotonsMB);
00587 m_incl->Add(g_inclPhotonsHT1);
00588 m_incl->Add(g_inclPhotonsHT2);
00589
00590 m_incl->Fit(fitQCD,"R0");
00591
00592 m_incl->SetMinimum(1.e-11);
00593 m_incl->SetMaximum(10.);
00594 TCanvas *c_incl=new TCanvas("c_incl","c_incl",600,400);
00595 gPad->SetLogy();
00596 m_incl->Draw("ap");
00597 c_incl->SaveAs((dir+"inclPhotonYield.eps").Data());
00598 c_incl->SaveAs((dir+"inclPhotonYield.root").Data());
00599
00600
00601
00602
00603 TH1F *gammaYieldMBratio=new TH1F(*gammaYieldMB);
00604 gammaYieldMBratio->SetName("gammaYieldMBratio");
00605 TH1F *gammaYieldHT1ratio=new TH1F(*gammaYieldHT1);
00606 gammaYieldHT1ratio->SetName("gammaYieldHT1ratio");
00607 TH1F *gammaYieldHT2ratio=new TH1F(*gammaYieldHT2);
00608 gammaYieldHT2ratio->SetName("gammaYieldHT2ratio");
00609
00610 gammaYieldMBratio->Divide(pionYieldMBratio);
00611 gammaYieldHT1ratio->Divide(pionYieldHT1ratio);
00612 gammaYieldHT2ratio->Divide(pionYieldHT2ratio);
00613
00614
00615 getRatio(gammaYieldMBratio,effGammaMB,effGammaSingleMB,fit_piondecay);
00616 getRatio(gammaYieldHT1ratio,effGammaHT1,effGammaSingleHT1,fit_piondecay);
00617 getRatio(gammaYieldHT2ratio,effGammaHT2,effGammaSingleHT2,fit_piondecay);
00618
00619 TH1F *gammaYieldMBratio_incl=new TH1F(*gammaYieldMBratio);
00620 TH1F *gammaYieldHT1ratio_incl=new TH1F(*gammaYieldHT1ratio);
00621 TH1F *gammaYieldHT2ratio_incl=new TH1F(*gammaYieldHT2ratio);
00622
00623
00624 TH1F *gammaYieldMBratioNoErr=new TH1F(*gammaYieldMBratio);
00625 TH1F *gammaYieldHT1ratioNoErr=new TH1F(*gammaYieldHT1ratio);
00626 TH1F *gammaYieldHT2ratioNoErr=new TH1F(*gammaYieldHT2ratio);
00627 for(int i=1;i<=gammaYieldMBratioNoErr->GetNbinsX();i++){
00628 gammaYieldMBratioNoErr->SetBinError(i,0.);
00629 }
00630 for(int i=1;i<=gammaYieldHT1ratioNoErr->GetNbinsX();i++){
00631 gammaYieldHT1ratioNoErr->SetBinError(i,0.);
00632 }
00633 for(int i=1;i<=gammaYieldHT2ratioNoErr->GetNbinsX();i++){
00634 gammaYieldHT2ratioNoErr->SetBinError(i,0.);
00635 }
00636
00637 TGraphErrors *g_ratioMB=new TGraphErrors(gammaYieldMBratio);
00638 g_ratioMB->SetName("g_ratioMB");
00639 TGraphErrors *g_ratioHT1=new TGraphErrors(gammaYieldHT1ratio);
00640 g_ratioHT1->SetName("g_ratioHT1");
00641 TGraphErrors *g_ratioHT2=new TGraphErrors(gammaYieldHT2ratio);
00642 g_ratioHT2->SetName("g_ratioHT2");
00643
00644
00645 removeThesePoints(g_ratioMB,1);
00646 removeThesePoints(g_ratioHT1,2);
00647 removeThesePoints(g_ratioHT2,3);
00648
00649 TCanvas *c_ratio=new TCanvas("c_ratio","c_ratio",400,200);
00650
00651 TMultiGraph *m_ratio=new TMultiGraph("m_ratio","d+Au 2003;p_{T};#gamma/#pi^{0}");
00652 m_ratio->Add(g_ratioMB);
00653 m_ratio->Add(g_ratioHT1);
00654 m_ratio->Add(g_ratioHT2);
00655
00656 m_ratio->Draw("ap");
00657 m_ratio->SetMinimum(.0);
00658 m_ratio->SetMaximum(2.);
00659
00660 TLegend *leg3=new TLegend(.35,.65,.65,.85);
00661 leg3->AddEntry(g_ratioMB,"minimum bias","p");
00662 leg3->AddEntry(g_ratioHT1,"hightower 1","p");
00663 leg3->AddEntry(g_ratioHT2,"hightower 2","p");
00664 leg3->AddEntry(fit_decay,"decay background (total)","l");
00665 leg3->AddEntry(fit_piondecay,"decay background (#pi^{0})","l");
00666 leg3->SetFillColor(0);
00667 leg3->Draw("same");
00668
00669 fit_decay->SetLineColor(13);
00670 fit_decay->SetLineWidth(1);
00671 fit_decay->SetLineColor(1);
00672 fit_decay->Draw("same");
00673 fit_piondecay->SetLineColor(13);
00674 fit_piondecay->SetLineWidth(1);
00675 fit_piondecay->SetLineStyle(2);
00676 fit_piondecay->SetLineColor(1);
00677 fit_piondecay->Draw("same");
00678
00679
00680 c_ratio->SaveAs((dir+"gammaOverPion.eps").Data());
00681 c_ratio->SaveAs((dir+"gammaOverPion.root").Data());
00682
00683
00684
00685 gammaYieldMBratio_incl->Multiply(pionYieldMBratio);
00686 gammaYieldHT1ratio_incl->Multiply(pionYieldHT1ratio);
00687 gammaYieldHT2ratio_incl->Multiply(pionYieldHT2ratio);
00688
00689 TGraphErrors *g_incl_corrMB=new TGraphErrors(gammaYieldMBratio_incl);
00690 TGraphErrors *g_incl_corrHT1=new TGraphErrors(gammaYieldHT1ratio_incl);
00691 TGraphErrors *g_incl_corrHT2=new TGraphErrors(gammaYieldHT2ratio_incl);
00692
00693 TCanvas *c_incl_corr=new TCanvas("c_incl_corr","c_incl_corr",400,300);
00694 gPad->SetLogy();
00695 TMultiGraph *m_incl_corr=new TMultiGraph();
00696 m_incl_corr->Add(g_incl_corrMB);
00697 m_incl_corr->Add(g_incl_corrHT1);
00698 m_incl_corr->Add(g_incl_corrHT2);
00699
00700 m_incl_corr->SetMinimum(1.e-11);
00701 m_incl_corr->SetMaximum(1.);
00702
00703 m_incl_corr->Draw("apX");
00704 c_incl_corr->SaveAs((dir+"inclPhotonYieldCorr.eps").Data());
00705 c_incl_corr->SaveAs((dir+"inclPhotonYieldCorr.root").Data());
00706
00707 TCanvas *c_doubleratio=new TCanvas("c_doubleratio","c_doubleratio",400,300);
00708 gStyle->SetOptStat(0);
00709 c_doubleratio->cd(1);
00710
00711 TH1F *gammaYieldMBdoubleratio=new TH1F(*gammaYieldMBratio);
00712 TH1F *gammaYieldHT1doubleratio=new TH1F(*gammaYieldHT1ratio);
00713 TH1F *gammaYieldHT2doubleratio=new TH1F(*gammaYieldHT2ratio);
00714
00715 gammaYieldMBdoubleratio->Divide(fit_decay);
00716 gammaYieldHT1doubleratio->Divide(fit_decay);
00717 gammaYieldHT2doubleratio->Divide(fit_decay);
00718
00719 TGraphErrors *g_doubleRatioMB=new TGraphErrors(gammaYieldMBdoubleratio);
00720 g_doubleRatioMB->SetName("g_doubleRatioMB");
00721 g_doubleRatioMB->SetMarkerStyle(8);
00722 TGraphErrors *g_doubleRatioHT1=new TGraphErrors(gammaYieldHT1doubleratio);
00723 g_doubleRatioHT1->SetName("g_doubleRatioHT1");
00724 g_doubleRatioHT1->SetMarkerStyle(8);
00725 TGraphErrors *g_doubleRatioHT2=new TGraphErrors(gammaYieldHT2doubleratio);
00726 g_doubleRatioHT2->SetName("g_doubleRatioHT2");
00727 g_doubleRatioHT2->SetMarkerStyle(8);
00728
00729 removeThesePoints(g_doubleRatioMB,1);
00730 removeThesePoints(g_doubleRatioHT1,2);
00731 removeThesePoints(g_doubleRatioHT2,3);
00732
00733 TMultiGraph *m_doubleratio=new TMultiGraph();
00734 m_doubleratio->SetMinimum(.5);
00735 m_doubleratio->SetMaximum(2.75);
00736
00737 cout<<endl;
00738 g_doubleRatioHT1->Print();
00739 cout<<endl<<endl;
00740 g_doubleRatioHT2->Print();
00741 cout<<endl;
00742
00743
00744 m_doubleratio->Add(g_doubleRatioHT1,"p");
00745 m_doubleratio->Add(g_doubleRatioHT2,"p");
00746
00747 g_photonpqcd->SetLineWidth(2);
00748 g_photonpqcd->SetLineColor(2);
00749 g_photonpqcd05->SetLineWidth(2);
00750 g_photonpqcd05->SetLineColor(2);
00751 g_photonpqcd05->SetLineStyle(2);
00752 g_photonpqcd2->SetLineWidth(2);
00753 g_photonpqcd2->SetLineColor(2);
00754 g_photonpqcd2->SetLineStyle(2);
00755
00756 m_doubleratio->Add(g_photonpqcd,"c");
00757 m_doubleratio->Add(g_photonpqcd05,"c");
00758 m_doubleratio->Add(g_photonpqcd2,"c");
00759
00760
00761
00762
00763 m_doubleratio->Draw("a");
00764
00765 m_doubleratio->GetXaxis()->SetTitle("p_{T} (GeV/c)");
00766 m_doubleratio->GetYaxis()->SetTitle("1 + #gamma_{dir}/#gamma_{incl}");
00767 m_doubleratio->GetXaxis()->SetRangeUser(2.,16.);
00768
00769
00770 TLegend *leg5=new TLegend(.15,.6,.6,.8);
00771 leg5->AddEntry(g_doubleRatioHT1,"hightower-1","p");
00772 leg5->AddEntry(g_doubleRatioHT2,"hightower-2","p");
00773 leg5->AddEntry(g_photonpqcd,"NLO (CTEQ6+KKP) #mu=p_{T}","l");
00774 leg5->AddEntry(g_photonpqcd05,"#mu=2p_{T}, #mu=p_{T}/2","l");
00775 leg5->SetFillColor(0);
00776 leg5->Draw("same");
00777
00778 c_doubleratio->cd(0);
00779 c_doubleratio->SaveAs((dir+"gammaDoubleRatio.eps").Data());
00780 c_doubleratio->SaveAs((dir+"gammaDoubleRatio.root").Data());
00781
00782
00783 TCanvas *c_dirphoton=new TCanvas("c_dirphoton","c_dirphoton",400,300);
00784 gStyle->SetOptStat(0);
00785 c_dirphoton->cd(1);
00786
00787 TH1F *dirphotonYieldMB=new TH1F(*gammaYieldMBdoubleratio);
00788 TH1F *dirphotonYieldHT1=new TH1F(*gammaYieldHT1doubleratio);
00789 TH1F *dirphotonYieldHT2=new TH1F(*gammaYieldHT2doubleratio);
00790
00791 TH1F *dirphotonYieldMBnoErr=new TH1F(*dirphotonYieldMB);
00792 for(int i=1;i<=dirphotonYieldMBnoErr->GetNbinsX();i++){
00793 dirphotonYieldMBnoErr->SetBinError(i,0.);
00794 }
00795 TH1F *dirphotonYieldHT1noErr=new TH1F(*dirphotonYieldHT1);
00796 for(int i=1;i<=dirphotonYieldHT1noErr->GetNbinsX();i++){
00797 dirphotonYieldHT1noErr->SetBinError(i,0.);
00798 }
00799 TH1F *dirphotonYieldHT2noErr=new TH1F(*dirphotonYieldHT2);
00800 for(int i=1;i<=dirphotonYieldHT1noErr->GetNbinsX();i++){
00801 dirphotonYieldHT2noErr->SetBinError(i,0.);
00802 }
00803
00804 TF1 *f_unity=new TF1("f_unity","1.",0.,15.);
00805 dirphotonYieldMB->Add(f_unity,-1.);
00806 dirphotonYieldHT1->Add(f_unity,-1.);
00807 dirphotonYieldHT2->Add(f_unity,-1.);
00808
00809
00810 dirphotonYieldMB->Divide(dirphotonYieldMBnoErr);
00811 dirphotonYieldHT1->Divide(dirphotonYieldHT1noErr);
00812 dirphotonYieldHT2->Divide(dirphotonYieldHT2noErr);
00813 dirphotonYieldMB->Multiply(gammaYieldMBratioNoErr);
00814 dirphotonYieldHT1->Multiply(gammaYieldHT1ratioNoErr);
00815 dirphotonYieldHT2->Multiply(gammaYieldHT2ratioNoErr);
00816 dirphotonYieldMB->Multiply(pionXsMBnoErr);
00817 dirphotonYieldHT1->Multiply(pionXsHT1noErr);
00818 dirphotonYieldHT2->Multiply(pionXsHT2noErr);
00819
00820
00821 TGraphErrors *g_dirphotonMB=new TGraphErrors(dirphotonYieldMB);
00822 g_dirphotonMB->SetName("g_dirphotonMB");
00823 g_dirphotonMB->SetMarkerStyle(8);
00824 TGraphErrors *g_dirphotonHT1=new TGraphErrors(dirphotonYieldHT1);
00825 g_dirphotonHT1->SetName("g_dirphotonHT1");
00826 g_dirphotonHT1->SetMarkerStyle(8);
00827 TGraphErrors *g_dirphotonHT2=new TGraphErrors(dirphotonYieldHT2);
00828 g_dirphotonHT2->SetName("g_dirphotonHT2");
00829 g_dirphotonHT2->SetMarkerStyle(8);
00830
00831
00832 removeThesePoints(g_dirphotonMB,1);
00833 removeThesePoints(g_dirphotonHT1,2);
00834 removeThesePoints(g_dirphotonHT2,3);
00835
00836 gPad->SetLogy();
00837
00838 TMultiGraph *m_dirphoton=new TMultiGraph();
00839 m_dirphoton->SetName("m_dirphoton");
00840 m_dirphoton->SetMinimum(1.0e-11);
00841 m_dirphoton->SetMaximum(0.1);
00842
00843 m_dirphoton->Add(g_dirgamma,"c");
00844 m_dirphoton->Add(g_dirgamma05,"c");
00845 m_dirphoton->Add(g_dirgamma2,"c");
00846
00847 cout<<"direct photons:"<<endl;
00848 g_dirphotonHT1->Print();
00849 cout<<endl;
00850 g_dirphotonHT2->Print();
00851 cout<<endl;
00852
00853 m_dirphoton->Add(g_dirphotonHT1,"p");
00854 m_dirphoton->Add(g_dirphotonHT2,"p");
00855
00856 m_dirphoton->Draw("a");
00857
00858 m_dirphoton->GetXaxis()->SetTitle("p_{T} (GeV/c)");
00859 m_dirphoton->GetYaxis()->SetTitle("1 - R^{-1}");
00860 m_dirphoton->GetXaxis()->SetRangeUser(2.,16.);
00861
00862
00863 TLegend *leg5=new TLegend(.15,.6,.6,.8);
00864 leg5->AddEntry(g_dirphotonHT1,"hightower-1","p");
00865 leg5->AddEntry(g_dirphotonHT2,"hightower-2","p");
00866 leg5->SetFillColor(0);
00867 leg5->Draw("same");
00868
00869 c_dirphoton->cd(0);
00870 c_dirphoton->SaveAs((dir+"gammaDirPhoton.eps").Data());
00871 c_dirphoton->SaveAs((dir+"gammaDirPhoton.root").Data());
00872
00873
00874
00875
00876
00877 return;
00878 }
