PWG2/SPECTRA/Fit/PlotCombine7TeV.C

   1 class AliOADBContainer;
   2 class AliOADBPWG2Spectra;
   3 Int_t colors[]={ kBlack,kRed,kBlue+2,kGreen+2, kSpring, kTeal, kAzure};
   4 Int_t markers[]={21,22,23,24,25,26,27,28};
   5 Float_t min[]={0.1,0.2,0.3,-1.0,-1.0,-1.0,0.5,0.5,0.8,-1.0,-1.0,-1.0};
   6 Float_t max[]={0.5,0.5,0.5,-1.0,-1.0,-1.0,2.0,2.0,2.5,1.5,1.5,1.7};
   7 Float_t min2[]={0.1,0.2,0.3,-1.0,-1.0,-1.0,0.5,0.5,0.8,-1.0,-1.0,-1.0};
   8 Float_t max2[]={0.5,0.5,0.5,-1.0,-1.0,-1.0,2.0,2.0,2.5,0.1,0.1,0.1};
   9 TString mcmodels[]={"Phojet","Pythia109","Pythia320"};
  10 const char * fgkDetectorNames[] = {"ITS", "ITSTPC", "TPC", "TOF", "TOFTPC", "Dummy", "Dummy"};
  11 const char * fgkPidTypeNames[]  = {"GaussFit", "NSigma", "Bayes", "Kinks"};
  12 const char * partext[]={"#pi^{+}","K^{+}","p","#pi^{-}","K^{-}","#bar{p}"};
  13 void ALICEWorkInProgress(TCanvas *c,TString today="11/05/2010", TString label = "ALICE performance");
  14 TH1D* DrawHisto(AliOADBContainer* fOADBContainer, TCanvas* c1,Int_t ana,Int_t method,Int_t par,Int_t charge,Int_t color, Int_t marker,TLegend* leg, TString estimator ,Int_t bin);
  15 Double_t myLevyPt(Double_t *pt, Double_t *par);
  16 TH1D* CombineSpectra(TList* l,Int_t startpoint, Float_t min* , Float_t* max);
  17 void DrawandFit(TH1D* hist,TCanvas* c1,Int_t par,TCanvas* c1ratio,Float_t* yileds,Float_t* yiledserr);
  18 void GetMCratios(TString filename, TList* ktopi,TList* ptopi);
  19 void plotRatios(TCanvas* c1, TH1D* data , TList* mc,TLegend* leg ,Int_t opt=0);
  20 TH1D* Plotratiodatafit(TH1D* hist1, TF1* fun);
  21 TCanvas* Ktopi(Float_t value,Float_t err);
  22 void DrawMulRatios(TH1D** bin1,TH1D** bin2,TCanvas* c,TLegend* leg);
  23
  24
  25 void PlotCombine7TeV()
  26 {
  27         gSystem->Load("libCore.so");
  28         gSystem->Load("libGeom.so");
  29         gSystem->Load("libPhysics.so");
  30         gSystem->Load("libVMC");
  31         gSystem->Load("libTree");
  32         gSystem->Load("libProof");
  33         gSystem->Load("libMatrix");
  34         gSystem->Load("libSTEERBase");
  35         gSystem->Load("libESD");
  36         gSystem->Load("libAOD");
  37         gSystem->Load("libANALYSIS");
  38         gSystem->Load("libOADB");
  39         gSystem->Load("libANALYSISalice");
  40         gSystem->Load("libTENDER");
  41         gSystem->Load("libCORRFW");
  42         gSystem->Load("libPWG0base");
  43         gSystem->Load("libMinuit");
  44         gSystem->Load("libPWG2spectra");
  45
  46         gROOT->SetStyle("Plain");
  47         gStyle->SetPalette(1,0);
  48         gStyle->SetOptFit(0);
  49         gStyle->SetOptStat(0);
  50         Int_t det[]={AliOADBPWG2Spectra::kITSsa,AliOADBPWG2Spectra::kITSTPC,AliOADBPWG2Spectra::kTOF,AliOADBPWG2Spectra::kTOFTPC} ;
  51         Int_t method[]={AliOADBPWG2Spectra::kGaussFit,AliOADBPWG2Spectra::kGaussFit,AliOADBPWG2Spectra::kGaussFit,AliOADBPWG2Spectra::kNSigma };
  52
  53
  54          TString fileName = AliOADBPWG2Spectra::GetOADBPWG2SpectraFileName();
  55         TFile * f = new TFile (fileName);
  56         AliOADBContainer* fOADBContainer = (AliOADBContainer*) f->Get("Corrected");
  57         f->Close();
  58         if(!fOADBContainer)
  59                 return;
  60         AliOADBPWG2Spectra* fOADBSpectra = (AliOADBPWG2Spectra*) fOADBContainer->GetObject(116562);
  61         if(!fOADBSpectra)
  62                 return;
  63         fOADBSpectra->Print();
  64         TCanvas* c1= new TCanvas("pos","pos",1200,800);
  65         c1->cd()->SetLogy();
  66         TLegend* Leg1 = new TLegend(0.1,0.1,0.45,0.45,"","NDC");
  67         Leg1->SetNColumns(3);
  68         Leg1->SetTextSize(0.027);
  69         Leg1->SetFillStyle(kFALSE);
  70         Leg1->SetLineColor(kWhite);
  71         Leg1->SetBorderSize(0);
  72
  73
  74         TCanvas* c2= new TCanvas("neg","neg",1200,800);
  75         c2->cd()->SetLogy();
  76         TLegend* Leg2 = new TLegend(0.1,0.1,0.45,0.45,"","NDC");
  77         Leg2->SetNColumns(3);
  78         Leg2->SetTextSize(0.027);
  79         Leg2->SetFillStyle(kFALSE);
  80         Leg2->SetLineColor(kWhite);
  81         Leg2->SetBorderSize(0);
  82
  83
  84         TDatime date;
  85         TList* posparthist=new TList();
  86         TList* negparthist=new TList();
  87
  88         TList* posparthistbin1=new TList();
  89         TList* negparthistbin1=new TList();
  90
  91         TList* posparthistbin4=new TList();
  92         TList* negparthistbin4=new TList();
  93
  94         for(int i=0;i<4;i++)// ITS , ITSTPC, TOF TPCTOF
  95         {
  96                 for(int j=0;j<3;j++) //pion kaon proton
  97                 {
  98                         posparthist->Add(DrawHisto(fOADBSpectra,c1,det[i],method[i],j,0,colors[j],markers[i],Leg1));
  99                         negparthist->Add(DrawHisto(fOADBSpectra,c2,det[i],method[i],j,1,colors[j],markers[i],Leg2));
 100                         posparthistbin1->Add(DrawHisto(fOADBSpectra,0x0,det[i],method[i],j,0,colors[j],markers[i],0x0,"SPD2",1));
 101                         negparthistbin1->Add(DrawHisto(fOADBSpectra,0x0,det[i],method[i],j,1,colors[j],markers[i],0x0,"SPD2",1));
 102                         posparthistbin4->Add(DrawHisto(fOADBSpectra,0x0,det[i],method[i],j,0,colors[j],markers[i],0x0,"SPD2",4));
 103                         negparthistbin4->Add(DrawHisto(fOADBSpectra,0x0,det[i],method[i],j,1,colors[j],markers[i],0x0,"SPD2",4));
 104                 }
 105         }
 106         //negparthistbin4->ls();
 107         ALICEWorkInProgress(c1,Form("%d/%d/%d",date.GetDay(),date.GetMonth(),date.GetYear()),"ALICE preliminary");
 108         Leg1->Draw();
 109         ALICEWorkInProgress(c2,Form("%d/%d/%d",date.GetDay(),date.GetMonth(),date.GetYear()),"ALICE preliminary");
 110         Leg2->Draw();
 111
 112         c1->SaveAs("Pos.eps");
 113         c2->SaveAs("Neg.eps");
 114
 115
 116         TH1D* combinepos[3];
 117         TH1D* combineneg[3];
 118         TH1D* combineposbin1[3];
 119         TH1D* combinenegbin1[3];
 120         TH1D* combineposbin4[3];
 121         TH1D* combinenegbin4[3];
 122         TCanvas* c1fit= new TCanvas("posfit","posfit",1200,800);
 123         c1fit->cd()->SetLogy();
 124         TCanvas* c2fit= new TCanvas("negfit","negfit",1200,800);
 125         c2fit->cd()->SetLogy();
 126         TCanvas* c1fitratio= new TCanvas("posfitratio","posfitrato",1200,800);
 127         TCanvas* c2fitratio= new TCanvas("negfitratio","negfitratio",1200,800);
 128         Float_t yileds[3]={0.0,0.0,0.0};
 129         Float_t yiledserr[3]={0.0,0.0,0.0};
 130         for(int j=0;j<3;j++) //pion kaon proton
 131         {
 132                         combinepos[j]=CombineSpectra(posparthist,j, min,max);
 133                         combineposbin1[j]=CombineSpectra(posparthistbin1,j,min2,max2);
 134                         combineposbin4[j]=CombineSpectra(posparthistbin4,j,min2,max2);
 135                         DrawandFit(combinepos[j],c1fit,j,c1fitratio,yileds,yiledserr);
 136                         combinepos[j]->Sumw2();
 137                         combineposbin1[j]->Sumw2();
 138                         combineposbin4[j]->Sumw2();
 139                         combineposbin1[j]->Divide(combinepos[j]);
 140                         combineposbin4[j]->Divide(combinepos[j]);
 141                         combineneg[j]=CombineSpectra(negparthist,j, min,max);
 142                         combinenegbin1[j]=CombineSpectra(negparthistbin1,j,min2,max2);
 143                         combinenegbin4[j]=CombineSpectra(negparthistbin4,j,min2,max2);
 144                         DrawandFit(combineneg[j],c2fit,j,c2fitratio,yileds,yiledserr);
 145                         combineneg[j]->Sumw2();
 146                         combinenegbin1[j]->Sumw2();
 147                         combinenegbin4[j]->Sumw2();
 148                         combinenegbin1[j]->Divide(combinepos[j]);
 149                         combinenegbin4[j]->Divide(combinepos[j]);
 150                         combinepos[j]->Add(combineneg[j]);
 151
 152         }
 153         ALICEWorkInProgress(c1fit,Form("%d/%d/%d",date.GetDay(),date.GetMonth(),date.GetYear()),"ALICE preliminary");
 154         ALICEWorkInProgress(c2fit,Form("%d/%d/%d",date.GetDay(),date.GetMonth(),date.GetYear()),"ALICE preliminary");
 155         c1fit->SaveAs("posfit.eps");
 156         c2fit->SaveAs("negfit.eps");
 157         c1fitratio->SaveAs("posfitratio.eps");
 158         c2fitratio->SaveAs("negfitratio.eps");
 159         combinepos[1]->Divide(combinepos[0]);
 160         combinepos[2]->Divide(combinepos[0]);
 161
 162         TString mcfiles[3]={"/home/marek/Analysis/Spectra/7TeVLHC10b/MCmodels/AnalyseFastPhojet7000AccCut.root","/home/marek/Analysis/Spectra/7TeVLHC10b/MCmodels/AnalyseFastPythia7000AccCutTune109.root","/home/marek/Analysis/Spectra/7TeVLHC10b/MCmodels/AnalyseFastPythia7000AccCutTune320.root"};
 163         TList* ktopi=new TList();
 164         TList* ptopi=new TList();
 165         //ktopi->ls();
 166         for(int i=0;i<3;i++)
 167                 GetMCratios(mcfiles[i],ktopi,ptopi);
 168         //ktopi->ls();
 169         TCanvas* c1MC= new TCanvas("K/pi","K/pi",1200,800);
 170         c1MC->cd();
 171
 172         TLegend* Leg3 = new TLegend(0.1,0.35,0.45,0.75,"","NDC");
 173         Leg3->SetTextSize(0.027);
 174         Leg3->SetFillStyle(kFALSE);
 175         Leg3->SetLineColor(kWhite);
 176         Leg3->SetBorderSize(0);
 177         plotRatios(c1MC,combinepos[1],ktopi,Leg3,1);
 178         Leg3->Draw();
 179         TCanvas* c2MC= new TCanvas("p/pi","p/pi",1200,800);
 180         TLegend* Leg4 = new TLegend(0.1,0.25,0.45,0.65,"","NDC");
 181         Leg4->SetTextSize(0.027);
 182         Leg4->SetFillStyle(kFALSE);
 183         Leg4->SetLineColor(kWhite);
 184         Leg4->SetBorderSize(0);
 185         plotRatios(c2MC,combinepos[2],ptopi,Leg4,0);
 186         Leg4->Draw();
 187         ALICEWorkInProgress(c1MC,Form("%d/%d/%d",date.GetDay(),date.GetMonth(),date.GetYear()),"ALICE preliminary");
 188         ALICEWorkInProgress(c2MC,Form("%d/%d/%d",date.GetDay(),date.GetMonth(),date.GetYear()),"ALICE preliminary");
 189
 190         c1MC->SaveAs("spectraKtopiMC.eps");
 191         c2MC->SaveAs("spectraptopiMC.eps");
 192
 193         TCanvas* cK2pi=0x0;
 194         if(yileds[0]>0.0)
 195         {
 196                 cK2pi=Ktopi(yileds[1]/yileds[0],yileds[1]/yileds[0]*TMath::Sqrt((yiledserr[0]/yileds[0])*(yiledserr[0]/yileds[0])+(yiledserr[1]/yileds[1])*(yiledserr[1]/yileds[1])));
 197                 ALICEWorkInProgress(cK2pi,Form("%d/%d/%d",date.GetDay(),date.GetMonth(),date.GetYear()),"ALICE preliminary");
 198         }
 199         cK2pi->SaveAs("K2pi.eps");
 200         TCanvas* c1mul= new TCanvas("posmul","posmul",1200,800);
 201
 202         c1mul->cd();
 203         TLegend* Leg5 = new TLegend(0.2,0.2,0.45,0.45,"","NDC");
 204         Leg5->SetNColumns(2);
 205         Leg5->SetTextSize(0.027);
 206         Leg5->SetFillStyle(kFALSE);
 207         Leg5->SetLineColor(kWhite);
 208         Leg5->SetBorderSize(0);
 209         DrawMulRatios(combineposbin1,combineposbin4,c1mul,Leg5,0);
 210         Leg5->Draw();
 211         TCanvas* c2mul= new TCanvas("negmul","negmul",1200,800);
 212         TLegend* Leg6 = new TLegend(0.2,0.2,0.45,0.45,"","NDC");
 213         Leg6->SetNColumns(2);
 214         Leg6->SetTextSize(0.027);
 215         Leg6->SetFillStyle(kFALSE);
 216         Leg6->SetLineColor(kWhite);
 217         Leg6->SetBorderSize(0);
 218         DrawMulRatios(combinenegbin1,combinenegbin4,c2mul,Leg6,1);
 219         Leg6->Draw();
 220         c2mul->cd();
 221         ALICEWorkInProgress(c1mul,Form("%d/%d/%d",date.GetDay(),date.GetMonth(),date.GetYear()),"ALICE preliminary");
 222         ALICEWorkInProgress(c2mul,Form("%d/%d/%d",date.GetDay(),date.GetMonth(),date.GetYear()),"ALICE preliminary");
 223         c1mul->SaveAs("MulPos.eps");
 224         c2mul->SaveAs("MulNeg.eps");
 225 }
 226 void ALICEWorkInProgress(TCanvas *c,TString today, TString label)
 227 {
 228                 c->cd();
 229           TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo",0.72,0.72,0.89,0.89);
 230           //TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo",0.72,0.62,0.85,0.75);
 231           myPadLogo->SetFillColor(0);
 232           myPadLogo->SetBorderMode(0);
 233           myPadLogo->SetBorderSize(2);
 234           myPadLogo->SetFrameBorderMode(0);
 235           myPadLogo->SetLeftMargin(0.0);
 236           myPadLogo->SetTopMargin(0.0);
 237           myPadLogo->SetBottomMargin(0.0);
 238           myPadLogo->SetRightMargin(0.0);
 239           myPadLogo->SetFillStyle(0);
 240           myPadLogo->Draw();
 241           myPadLogo->cd();
 242           TASImage *myAliceLogo = new TASImage("alice_logo_trans.png");
 243           myAliceLogo->Draw();
 244           c->cd();
 245           TPaveText* t1=new TPaveText(0.418103, 0.837798, 0.656609, 0.888393,"NDC");
 246           t1->SetFillStyle(0);
 247           t1->SetBorderSize(0);
 248           t1->AddText(0.,0.,label);
 249           t1->SetTextColor(kRed);
 250           t1->SetTextFont(42);
 251           t1->SetTextSize(0.04);
 252           t1->Draw();
 253           TPaveText* t2=new TPaveText(0.418103, 0.80, 0.656609, 0.84,"NDC");
 254           t2->SetFillStyle(0);
 255           t2->SetBorderSize(0);
 256           t2->AddText(0.,0.,"pp at #sqrt{s} = 7 TeV");
 257           t2->SetTextColor(kRed);
 258           t2->SetTextFont(42);
 259           t2->SetTextSize(0.027);
 260           t2->Draw();
 261           TPaveText* t3=new TPaveText(0.418103, 0.76, 0.656609, 0.80,"NDC");
 262           t3->SetFillStyle(0);
 263           t3->SetBorderSize(0);
 264          // t3->AddText(0.,0.,"Statistical and systematic errors");
 265           t3->AddText(0.,0.,"Statistical  errors");
 266           t3->SetTextColor(kRed);
 267           t3->SetTextFont(42);
 268           t3->SetTextSize(0.027);
 269           t3->Draw();
 270            TPaveText* t2=new TPaveText(0.65,0.65,0.89,0.7,"NDC");
 271            t2->SetFillStyle(0);
 272            t2->SetBorderSize(0);
 273            t2->SetTextColor(kRed);
 274            t2->SetTextFont(52);
 275            t2->AddText(0.,0.,today.Data());
 276            t2->Draw();
 277 }
 278
 279 TH1D* DrawHisto(AliOADBPWG2Spectra* fOADBSpectra, TCanvas* c1, Int_t ana,Int_t method,Int_t par,Int_t charge,Int_t color, Int_t marker,TLegend* leg,TString estimator="MB" ,Int_t bin=-1)
 280 {
 281
 282    TH1D * h = 0;
 283     TString binname=estimator;
 284    if(estimator.CompareTo("MB")==0)
 285    {
 286                 cout<<fOADBSpectra->GetHistoName(ana,method,par,charge,"MB")<<endl;
 287                 h = fOADBSpectra->GetHisto(ana,method,par,charge,"MB");
 288    // Draw the selected histogram
 289         }
 290         else
 291         {
 292                 cout<<fOADBSpectra->GetHistoName(ana,method,par,charge,estimator.Data(),bin)<<endl;
 293                 h = fOADBSpectra->GetHisto(ana,method,par,charge,estimator.Data(),bin);
 294                 binname+=bin;
 295         }
 296    if(!h)
 297    {
 298      cout << "Cannot get pointer to histo" << endl;
 299      return 0x0;
 300    }
 301    h->GetYaxis()->SetRangeUser(0.001,10.0);
 302     h->GetXaxis()->SetRangeUser(0.0,3.0);
 303
 304     if(charge==0)
 305                 h->SetTitle(Form("Positive_%s",binname.Data()));
 306     else
 307                 h->SetTitle(Form("Negative_%s",binname.Data()));
 308     h->SetXTitle("p_{T} (GeV/c)");
 309      h->SetYTitle("1/N_{events} dN/dp_{T} |y|<0.5");
 310
 311    h->SetMarkerColor(color);
 312    h->SetMarkerStyle(marker);
 313
 314
 315
 316
 317  if(c1&&leg)
 318  {
 319         TString opt = "E1";
 320         c1->cd();
 321     c1->GetListOfPrimitives()->Print();
 322     if (c1->GetListOfPrimitives()->GetEntries()>0)
 323                 opt += "same";
 324         c1->Update();
 325         c1->Modified();
 326         c1->Update();
 327
 328         leg->AddEntry(h,Form("%s_{%s_%s}",partext[charge*3+par],fgkDetectorNames[ana],fgkPidTypeNames[method]),"p");
 329     h->Draw(opt.Data());
 330    // TCanvas::Update() draws the frame, after which it can be changed
 331 }
 332    return h;
 333 }
 334 Double_t myLevyPt(Double_t *pt, Double_t *par)
 335 {
 336
 337   Double_t pdNdy  = par[0];
 338   Double_t pTemp = par[1];
 339   Double_t pPower = par[2];
 340   Double_t pMass  = par[3];
 341   Double_t pBigCoef = ((pPower-1)*(pPower-2)) / (pPower*pTemp*(pPower*pTemp+pMass*(pPower-2)));
 342
 343   Double_t pInPower = 1.0 + (TMath::Sqrt(pt[0]*pt[0]+pMass*pMass)-pMass) / (pPower*pTemp);
 344   return pdNdy * pt[0] * pBigCoef * TMath::Power(pInPower,(-1.0)*pPower);
 345 }
 346 TH1D* CombineSpectra(TList* l,Int_t startpoint, Float_t* min, Float_t* max)
 347 {
 348         //l->ls();
 349 //      ((TH1D*)l->At(0))->Print("all");
 350         TH1D* hcombine=new TH1D(*((TH1D*)l->At(0)));
 351         hcombine->SetName(Form("%s_combine",((TH1D*)l->At(startpoint))->GetName()));
 352         //cout<<hcombine->GetName()<<endl;
 353         for (int i=1; i<=hcombine->GetNbinsX();i++)
 354         {
 355                 Float_t pt=hcombine->GetXaxis()->GetBinCenter(i);
 356                 Int_t use[4]={1,1,1,1};
 357                 Float_t values[4]={0.0,0.0,0.0,0.0};
 358                 Float_t errors[4]={-1.0,-1.0,-1.0,-1.0};
 359         //      cout<<pt<<endl;
 360                 for(int j=0;j<4;j++)
 361                 {
 362                         if(min[3*j+startpoint]>0.0&&min[3*j+startpoint]>pt)
 363                                 use[j]=0;
 364                         if(max[3*j+startpoint]>0.0&&max[3*j+startpoint]<pt)
 365                                 use[j]=0;
 366                         if(use[j])
 367                         {
 368                                 values[j]=((TH1D*)l->At(3*j+startpoint))->GetBinContent(i);
 369                                 if(values[j]>0.0)
 370                                         errors[j]=((TH1D*)l->At(3*j+startpoint))->GetBinError(i);
 371                                 else
 372                                 {
 373                                         values[j]=0.0;
 374                                         use[j]=0;
 375                                 }
 376                         }
 377                 }
 378                 if(use[0]+use[1]+use[2]+use[3]==0)
 379                 {
 380                                 hcombine->SetBinContent(i,0.0);
 381                                 hcombine->SetBinError(i,0.0);
 382                 }
 383                 else
 384                 {
 385                         hcombine->SetBinContent(i,(values[0]+values[1]+values[2]+values[3])/(use[0]+use[1]+use[2]+use[3]));
 386                         hcombine->SetBinError(i,TMath::MaxElement(4,errors));
 387                 }
 388                 //cout<<use[0]+use[1]+use[2]+use[3]<<endl;
 389         }
 390         return hcombine;
 391 }
 392 void DrawandFit(TH1D* hist,TCanvas* c1, Int_t par , TCanvas* cratio,Float_t* yileds,Float_t* yiledserr)
 393 {
 394         c1->cd();
 395         TF1 *f1 = new TF1(Form("%s_fit",hist->GetName()),myLevyPt,0.1,3.0,4);
 396         f1->SetParNames("dN/dy","T","n","mass");
 397         f1->FixParameter(3,AliPID::ParticleMass(par+2));
 398         f1->SetParameter(0,hist->GetMaximum()/2.0);
 399
 400         f1->SetParameter(1,0.1);
 401         f1->SetParLimits(1,0.01,10.0);
 402         f1->SetParameter(2,3.0);
 403         f1->SetParLimits(2,2.1,1000.0);
 404
 405         hist->SetMarkerColor(colors[par]);
 406     hist->SetMarkerStyle(markers[par]);
 407
 408
 409     hist->Fit(Form("%s_fit",hist->GetName()),"I0");
 410
 411     TString opt = "E1";
 412
 413      c1->GetListOfPrimitives()->Print();
 414      if (c1->GetListOfPrimitives()->GetEntries()>0)
 415                 opt += "same";
 416                 cout<<opt<<endl;
 417         hist->DrawCopy(opt.Data());
 418         f1->DrawCopy("Lsame");
 419         cratio->cd();
 420         Plotratiodatafit(hist,f1)->Draw(opt.Data());
 421         yileds[par]+=f1->GetParameter(0);
 422         yiledserr[par]=TMath::Sqrt(yiledserr[par]*yiledserr[par]+f1->GetParError(0)*f1->GetParError(0));
 423 }
 424 void GetMCratios(TString filename, TList* ktopi,TList* ptopi)
 425 {
 426                 TFile* f=TFile::Open(filename.Data());
 427                 if(!f)
 428                                         return;
 429                 TH1F* fpiplus=(TH1F*)f->Get("h1PrimariesPiPlusPtVar");
 430                 TH1F* fpiminus=(TH1F*)f->Get("h1PrimariesPiMinusPtVar");
 431                 TH1F* fKplus=(TH1F*)f->Get("h1PrimariesKPlusPtVar");
 432                 TH1F* fKminus=(TH1F*)f->Get("h1PrimariesKMinusPtVar");
 433                 TH1F* fpplus=(TH1F*)f->Get("h1PrimariesProtonPtVar");
 434                 TH1F* fpminus=(TH1F*)f->Get("h1PrimariesAntiProtonPtVar");
 435
 436                 fpiplus->Sumw2();
 437                 fpiminus->Sumw2();
 438                 fKplus->Sumw2();
 439                 fKminus->Sumw2();
 440                 fpplus->Sumw2();
 441                 fpminus->Sumw2();
 442
 443                 fpiplus->Add(fpiminus);
 444                 fKplus->Add(fKminus);
 445                 fpplus->Add(fpminus);
 446
 447                 fKplus->Divide(fpiplus);
 448                 fpplus->Divide(fpiplus);
 449
 450                 ktopi->Add(fKplus);
 451                 ptopi->Add(fpplus);
 452 }
 453 void plotRatios(TCanvas* c1, TH1D* data , TList* mc,TLegend* leg, Int_t opt)
 454 {
 455         c1->cd();
 456         data->SetMarkerStyle(markers[0]);
 457         data->SetMarkerColor(colors[0]);
 458         data->Draw("E1");
 459         leg->AddEntry(data,"data","p");
 460         if(opt==1)
 461                 data->SetTitle("K/#pi");
 462         else
 463                 data->SetTitle("p/#pi");
 464         data->GetYaxis()->SetTitle("");
 465         data->GetYaxis()->SetRangeUser(0.0,0.5);
 466         data->GetXaxis()->SetRangeUser(0.0,2.5);
 467
 468         for (int i=0;i<3;i++)
 469         {
 470                         TH1F* mchist=(TH1F*)mc->At(i);
 471                         leg->AddEntry(mchist,mcmodels[i].Data(),"p");
 472                         mchist->SetMarkerStyle(markers[i+1]);
 473                         mchist->SetMarkerColor(colors[i+1]);
 474                         mchist->Draw("E1same");
 475         }
 476 }
 477 TH1D* Plotratiodatafit(TH1D* hist1, TF1* fun)
 478 {
 479         TString name(hist1->GetName());
 480         TH1D* histratio=new TH1D(*hist1);
 481
 482         histratio->SetName(Form("ratio%s",name.Data()));
 483         //histratio->SetTitle("ratio");
 484         histratio->GetYaxis()->SetTitle("ratio data/fit");
 485         for (int i=1;i<histratio->GetNbinsX();i++)
 486         {
 487                 if(histratio->GetBinContent(i)>0.0)
 488                 {
 489                         histratio->SetBinContent(i,histratio->GetBinContent(i)/fun->Eval(histratio->GetXaxis()->GetBinCenter(i)));
 490                         histratio->SetBinError(i,histratio->GetBinError(i)/fun->Eval(histratio->GetXaxis()->GetBinCenter(i)));
 491                 }
 492         }
 493         histratio->GetYaxis()->SetRangeUser(0.8,1.2);
 494         return histratio;
 495 }
 496 TCanvas* Ktopi(Float_t value,Float_t err)
 497 {
 498 //=========Macro generated from canvas: c1/c1
 499 //=========  (Thu Oct 14 14:18:11 2010) by ROOT version5.26/00b
 500    TCanvas *c1 = new TCanvas("K2pi", "K2pi",0,44,1105,782);
 501    c1->Range(0.4639814,-0.024875,3.950796,0.223875);
 502    c1->SetFillColor(0);
 503    c1->SetBorderMode(0);
 504    c1->SetBorderSize(2);
 505    c1->SetLogx();
 506    c1->SetTickx(1);
 507    c1->SetTicky(1);
 508    c1->SetFrameBorderMode(0);
 509    c1->SetFrameBorderMode(0);
 510
 511    TH2F *__1 = new TH2F("__1","",4002,1.5,9000.5,199,0,0.199);
 512    __1->SetDirectory(0);
 513    __1->SetStats(0);
 514    __1->GetXaxis()->SetTitle("#sqrt{s} (GeV)");
 515    __1->GetXaxis()->SetRangeUser(6,9002);
 516    __1->GetXaxis()->CenterTitle(true);
 517    __1->GetXaxis()->SetTitleSize(0.05);
 518    __1->GetXaxis()->SetTitleOffset(0.87);
 519    __1->GetYaxis()->SetTitle("K/#pi");
 520    __1->GetYaxis()->CenterTitle(true);
 521    __1->GetYaxis()->SetTitleSize(0.06);
 522    __1->GetYaxis()->SetTitleOffset(0.67);
 523    __1->Draw("");
 524
 525    TGraphErrors *gre = new TGraphErrors(1);
 526    gre->SetName("Graph");
 527    gre->SetTitle("Graph");
 528    gre->SetFillColor(1);
 529    gre->SetMarkerColor(3);
 530    gre->SetMarkerStyle(29);
 531    gre->SetMarkerSize(3.1);
 532    gre->SetPoint(0,199.9004,0.1030843);
 533    gre->SetPointError(0,0,0.008);
 534
 535    TH1F *Graph1 = new TH1F("Graph1","Graph",100,199.9,201.1);
 536    Graph1->SetMinimum(0.0934);
 537    Graph1->SetMaximum(0.1126);
 538    Graph1->SetDirectory(0);
 539    Graph1->SetStats(0);
 540    gre->SetHistogram(Graph1);
 541
 542    gre->Draw("p");
 543
 544    gre = new TGraphErrors(1);
 545    gre->SetName("Graph");
 546    gre->SetTitle("Graph");
 547    gre->SetFillColor(1);
 548    gre->SetMarkerColor(2);
 549    gre->SetMarkerStyle(20);
 550    gre->SetMarkerSize(2.5);
 551    //gre->SetPoint(0,897.8542,0.1235148);
 552    //gre->SetPointError(0,0,0.008624);
 553    //gre->SetPoint(0,897.8542,0.122);  // new vale from MF: K/pi ratio = 0.122 +- 0.01  Oct. 14, 2010
 554    //    gre->SetPointError(0,0,0.01);
 555
 556    gre->SetPoint(0,897.8542,0.12294);  // new vale from MF: K/pi ratio = 0.12294 +- 0.012  Feb.. 21, 2011
 557    gre->SetPointError(0,0,0.012);
 558
 559    gre->SetPoint(1,7000,value);  // new vale from MF: K/pi ratio = 0.12294 +- 0.012  Feb.. 21, 2011
 560    gre->SetPointError(1,0,err);
 561
 562
 563    TH1F *Graph2 = new TH1F("Graph2","Graph",100,899.9,901.1);
 564    Graph2->SetMinimum(0.1128512);
 565    Graph2->SetMaximum(0.1335488);
 566    Graph2->SetDirectory(0);
 567    Graph2->SetStats(0);
 568    gre->SetHistogram(Graph2);
 569
 570    gre->Draw("p");
 571
 572    gre = new TGraphErrors(1);
 573    gre->SetName("Graph");
 574    gre->SetTitle("Graph");
 575    gre->SetFillColor(1);
 576    gre->SetMarkerColor(4);
 577    gre->SetMarkerStyle(22);
 578    gre->SetMarkerSize(2.2);
 579    gre->SetPoint(0,17.24717,0.08050324);
 580    gre->SetPointError(0,0,0.0034);
 581
 582    TH1F *Graph3 = new TH1F("Graph3","Graph",100,17.2,18.4);
 583    Graph3->SetMinimum(0.07634);
 584    Graph3->SetMaximum(0.0845);
 585    Graph3->SetDirectory(0);
 586    Graph3->SetStats(0);
 587    gre->SetHistogram(Graph3);
 588
 589    gre->Draw("p");
 590
 591    gre = new TGraphErrors(4);
 592    gre->SetName("Graph");
 593    gre->SetTitle("Graph");
 594    gre->SetFillColor(1);
 595    gre->SetMarkerStyle(28);
 596    gre->SetMarkerSize(2.5);
 597    gre->SetPoint(0,300,0.105);
 598    gre->SetPointError(0,0,0.02);
 599    gre->SetPoint(1,534.998,0.112045);
 600    gre->SetPointError(1,0,0.01);
 601    gre->SetPoint(2,994.3588,0.1041596);
 602    gre->SetPointError(2,0,0.008);
 603    gre->SetPoint(3,1795.007,0.1131203);
 604    gre->SetPointError(3,0,0.005);
 605
 606    TH1F *Graph4 = new TH1F("Graph4","Graph",100,150,1950);
 607    Graph4->SetMinimum(0.081);
 608    Graph4->SetMaximum(0.129);
 609    Graph4->SetDirectory(0);
 610    Graph4->SetStats(0);
 611    gre->SetHistogram(Graph4);
 612
 613    gre->Draw("p");
 614
 615    gre = new TGraphErrors(1);
 616    gre->SetName("Graph");
 617    gre->SetTitle("Graph");
 618    gre->SetFillColor(1);
 619    gre->SetMarkerStyle(24);
 620    gre->SetMarkerSize(2.6);
 621    gre->SetPoint(0,546.8309,0.09519885);
 622    gre->SetPointError(0,0,0.0114);
 623
 624    TH1F *Graph5 = new TH1F("Graph5","Graph",100,544.9,546.1);
 625    Graph5->SetMinimum(0.08132);
 626    Graph5->SetMaximum(0.10868);
 627    Graph5->SetDirectory(0);
 628    Graph5->SetStats(0);
 629    gre->SetHistogram(Graph5);
 630
 631    gre->Draw("p");
 632    c1->Modified();
 633    c1->cd();
 634    c1->SetSelected(c1);
 635   // c1->ToggleToolBar();
 636   return c1;
 637 }
 638 void DrawMulRatios(TH1D** bin1,TH1D** bin2,TCanvas* c,TLegend* leg,Int_t charge=0)
 639 {
 640         c->cd();
 641         for(int i=0;i<3;i++)
 642         {
 643                 bin1[i]->GetYaxis()->SetTitle("ratio over MB spectra");
 644                 bin1[i]->GetYaxis()->SetRangeUser(0.0,4.0);
 645                 bin1[i]->SetMarkerStyle(markers[i]);
 646                 bin1[i]->SetMarkerColor(colors[i]);
 647                 leg->AddEntry(bin1[i],Form("%s_{bin1}",partext[i+charge*3],"p"));
 648                 if(i==0)
 649                         bin1[i]->Draw("E1");
 650                 else
 651                         bin1[i]->Draw("E1same");
 652                 bin2[i]->GetYaxis()->SetTitle("ratio over MB spectra");
 653                 bin2[i]->GetYaxis()->SetRangeUser(0.0,4.0);
 654                 bin2[i]->SetMarkerStyle(markers[i+3]);
 655                 bin2[i]->SetMarkerColor(colors[i]);
 656                 bin2[i]->Draw("E1same") ;
 657                 leg->AddEntry(bin2[i],Form("%s_{bin4}",partext[i+charge*3],"p"));
 658         }
 659 }