Removing try/catch-es
[u/mrichter/AliRoot.git] / 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 }