B2 analysis code
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / Nuclei / B2 / macros / RatioMult.C
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 // particle ratios as a function of multiplicity
17 // author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
18
19 #include <TROOT.h>
20 #include <TFile.h>
21 #include <TString.h>
22 #include <TGraphErrors.h>
23 #include <TF1.h>
24 #include <TCanvas.h>
25 #include <TMath.h>
26 #include <TStyle.h>
27
28 #include "B2.h"
29 #include "Config.h"
30
31 void GetRatio(Double_t x, Double_t y, Double_t errx, Double_t erry, Double_t& r, Double_t& rerr);
32 void Draw(TGraph* gr, Int_t marker, Int_t color, const TString& xtitle, const TString& ytitle);
33
34 void RatioMult(const TString& pSpectra = "~/alice/output/Proton-lhc10d-Mult-Spectra.root",
35                const TString& ptag     = "lhc10d",
36                const TString& dSpectra = "~/alice/output/Deuteron-lhc10d-Mult-Spectra.root",
37                const TString& dtag     = "lhc10d")
38 {
39 //
40 // particle ratios as a function of multiplicity (from fitting distributions)
41 //
42         using namespace B2mult;
43         
44         const Int_t kNpart = 2;
45         
46         const TString kProton[kNpart] = { "Proton", "AntiProton" };
47         const TString kDeuteron[kNpart] = { "Deuteron", "AntiDeuteron" };
48         
49         // open files
50         
51         TFile* fproton = new TFile(pSpectra.Data());
52         if(fproton->IsZombie()) exit(1);
53         
54         TFile* fdeuteron = new TFile(dSpectra.Data());
55         if(fdeuteron->IsZombie()) exit(1);
56         
57         // particle ratios
58         TGraphErrors*  grRatio[kNpart]    = { new TGraphErrors(), new TGraphErrors() };
59         TGraphErrors*  grMixRatio[kNpart] = { new TGraphErrors(), new TGraphErrors() };
60         
61         // model parameters
62         TGraphErrors*  grProtondNdy[kNpart]    = { new TGraphErrors(), new TGraphErrors() };
63         TGraphErrors*  grProtonN[kNpart]       = { new TGraphErrors(), new TGraphErrors() };
64         TGraphErrors*  grProtonC[kNpart]       = { new TGraphErrors(), new TGraphErrors() };
65         TGraphErrors*  grProtonChi2Ndf[kNpart] = { new TGraphErrors(), new TGraphErrors() };
66         //TGraphErrors*  grProtonProb[kNpart]    = { new TGraphErrors(), new TGraphErrors() };
67         
68         TGraphErrors*  grDeuterondNdy[kNpart]    = { new TGraphErrors(), new TGraphErrors() };
69         TGraphErrors*  grDeuteronN[kNpart]       = { new TGraphErrors(), new TGraphErrors() };
70         TGraphErrors*  grDeuteronC[kNpart]       = { new TGraphErrors(), new TGraphErrors() };
71         TGraphErrors*  grDeuteronChi2Ndf[kNpart] = { new TGraphErrors(), new TGraphErrors() };
72         //TGraphErrors*  grDeuteronProb[kNpart]    = { new TGraphErrors(), new TGraphErrors() };
73         
74         for(Int_t i=0; i<kNmult; ++i)
75         {
76                 TF1* fncProton[kNpart];
77                 TF1* fncDeuteron[kNpart];
78                 
79                 for(Int_t j=0; j<kNpart; ++j)
80                 {
81                         fncProton[j] = (TF1*)FindObj(fproton, ptag + "-" + kMultClass[i], kProton[j] + "_Fit_InvDiffYield_Pt");
82                         fncDeuteron[j] = (TF1*)FindObj(fdeuteron, dtag + "-" + kMultClass[i], kDeuteron[j] + "_Fit_InvDiffYield_Pt");
83                 }
84                 
85                 // integrated ratios
86                 
87                 Double_t dNdy[kNpart] = {fncProton[0]->GetParameter(0), fncDeuteron[0]->GetParameter(0)};
88                 Double_t dNdyErr[kNpart] = {fncProton[0]->GetParError(0), fncDeuteron[0]->GetParError(0)};
89                 
90                 Double_t dNdyBar[kNpart] = {fncProton[1]->GetParameter(0), fncDeuteron[1]->GetParameter(0)};
91                 Double_t dNdyBarErr[kNpart] = {fncProton[1]->GetParError(0), fncDeuteron[1]->GetParError(0)};
92                 
93                 // ratios
94                 Double_t ratio[kNpart], ratioErr[kNpart];
95                 for(Int_t j=0; j<kNpart; ++j)
96                 {
97                         GetRatio(dNdyBar[j], dNdy[j], dNdyBarErr[j], dNdyErr[j], ratio[j], ratioErr[j]);
98                         grRatio[j]->SetPoint(i, kKNOmult[i], ratio[j]);
99                         grRatio[j]->SetPointError(i, 0, ratioErr[j]);
100                 }
101                 
102                 // mixed ratios
103                 Double_t mixRatio[kNpart], mixRatioErr[kNpart];
104                 
105                 GetRatio(dNdy[1], dNdy[0], dNdyErr[1], dNdyErr[0], mixRatio[0], mixRatioErr[0]);
106                 GetRatio(dNdyBar[1], dNdyBar[0], dNdyBarErr[1], dNdyBarErr[0], mixRatio[1], mixRatioErr[1]);
107                 
108                 for(Int_t j=0; j<kNpart; ++j)
109                 {
110                         grMixRatio[j]->SetPoint(i, kKNOmult[i], mixRatio[j]);
111                         grMixRatio[j]->SetPointError(i, 0, mixRatioErr[j]);
112                 }
113                 
114                 // model parameters
115                 for(Int_t j=0; j<kNpart; ++j)
116                 {
117                         grProtondNdy[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetParameter(0));
118                         grProtondNdy[j]->SetPointError(i, 0, fncProton[j]->GetParError(0));
119                         
120                         grProtonN[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetParameter(1));
121                         grProtonN[j]->SetPointError(i, 0, fncProton[j]->GetParError(1));
122                         
123                         grProtonC[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetParameter(2));
124                         grProtonC[j]->SetPointError(i, 0, fncProton[j]->GetParError(2));
125                         
126                         grProtonChi2Ndf[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetChisquare()/fncProton[j]->GetNDF());
127                         //grProtonProb[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetProb());
128                         
129                         // deuteron
130                         grDeuterondNdy[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetParameter(0));
131                         grDeuterondNdy[j]->SetPointError(i, 0, fncDeuteron[j]->GetParError(0));
132                         
133                         grDeuteronN[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetParameter(1));
134                         grDeuteronN[j]->SetPointError(i, 0, fncDeuteron[j]->GetParError(1));
135                         
136                         grDeuteronC[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetParameter(2));
137                         grDeuteronC[j]->SetPointError(i, 0, fncDeuteron[j]->GetParError(2));
138                         
139                         grDeuteronChi2Ndf[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetChisquare()/fncDeuteron[j]->GetNDF());
140                         //grDeuteronProb[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetProb());
141                 }
142         }
143         
144         delete fproton;
145         delete fdeuteron;
146         
147         // draw
148         
149         TStyle* st = new TStyle();
150         
151         st->SetPadTickX(1);
152         st->SetPadTickY(1);
153         st->SetPadGridX(1);
154         st->SetPadGridY(1);
155         st->SetCanvasColor(0);
156         st->SetFrameBorderMode(0);
157         st->SetFrameFillColor(0);
158         st->SetLabelFont(62,"XYZ");
159         st->SetTitleFont(62,"XYZ");
160         
161         st->cd();
162         gROOT->ForceStyle();
163         
164         const Int_t kNCol = 4;
165         
166         TCanvas* c0 = new TCanvas("c0", "proton model parameters");
167         c0->Divide(kNCol,2);
168         
169         for(Int_t j=0; j<kNpart; ++j)
170         {
171                 c0->cd(kNCol*j+1);
172                 Draw(grProtondNdy[j], kFullCircle, kBlue, "z", "dN/dy");
173                 
174                 c0->cd(kNCol*j+2);
175                 Draw(grProtonN[j], kFullCircle, kBlue, "z", "n");
176                 
177                 c0->cd(kNCol*j+3);
178                 Draw(grProtonC[j], kFullCircle, kBlue, "z", "C");
179                 
180                 c0->cd(kNCol*j+4);
181                 Draw(grProtonChi2Ndf[j], kFullCircle, kBlue, "z", "#chi^{2}/ndf");
182                 
183         /*      c0->cd(kNCol*j+5);
184                 Draw(grProtonProb[j], kFullCircle, kBlue, "z", "Prob");
185         */
186         }
187         
188         TCanvas* c1 = new TCanvas("c1", "deuteron model parameters");
189         c1->Divide(kNCol,2);
190         
191         for(Int_t j=0; j<kNpart; ++j)
192         {
193                 c1->cd(kNCol*j+1);
194                 Draw(grDeuterondNdy[j], kFullCircle, kBlue, "z", "dN/dy");
195                 
196                 c1->cd(kNCol*j+2);
197                 Draw(grDeuteronN[j], kFullCircle, kBlue, "z", "n");
198                 
199                 c1->cd(kNCol*j+3);
200                 Draw(grDeuteronC[j], kFullCircle, kBlue, "z", "C");
201                 
202                 c1->cd(kNCol*j+4);
203                 Draw(grDeuteronChi2Ndf[j], kFullCircle, kBlue, "z", "#chi^{2}/ndf");
204                 
205         /*      c1->cd(kNCol*j+5);
206                 Draw(grDeuteronProb[j], kFullCircle, kBlue, "z", "Prob");
207         */
208         }
209         
210         TCanvas* c2 = new TCanvas("c2","particle ratios");
211         c2->Divide(2,2);
212         
213         c2->cd(1);
214         Draw(grRatio[0], kFullCircle, kRed, "z", "#bar{p}/p");
215         
216         c2->cd(2);
217         Draw(grRatio[1], kFullCircle, kRed, "z", "#bar{d}/d");
218         
219         c2->cd(3);
220         Draw(grMixRatio[0], kFullCircle, kRed, "z", "d/p");
221         
222         c2->cd(4);
223         Draw(grMixRatio[1], kFullCircle, kRed, "z", "#bar{d}/#bar{p}");
224 }
225
226 void GetRatio(Double_t x, Double_t y, Double_t errx, Double_t erry, Double_t& r, Double_t& rerr)
227 {
228 //
229 // get the ratio of x/y, its error and save in r and rerr
230 //
231         if(x == 0 || y == 0)
232         {
233                 r=0;
234                 rerr=1e+16;
235                 
236                 return;
237         }
238         
239         r = x/y;
240         rerr = r*TMath::Sqrt(TMath::Power(errx/x,2)+TMath::Power(erry/y,2));
241 }
242
243 void Draw(TGraph* gr, Int_t marker, Int_t color, const TString& xtitle, const TString& ytitle)
244 {
245 //
246 // draw a graph in current canvas
247 //
248         gr->SetMarkerStyle(marker);
249         gr->SetMarkerColor(color);
250         gr->SetLineColor(color);
251         gr->GetXaxis()->SetTitle(xtitle.Data());
252         gr->GetYaxis()->SetTitle(ytitle.Data());
253         gr->Draw("AP");
254 }