B2 analysis code
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / Nuclei / B2 / AliLnSpectra.cxx
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 // invariant differential yields and cross sections
17 // author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
18
19 #include <Riostream.h>
20 #include <TFile.h>
21 #include <TH1D.h>
22 #include <TString.h>
23 #include <TMath.h>
24 #include <TGraphErrors.h>
25 #include <TF1.h>
26
27 #include "AliLnSpectra.h"
28 #include "B2.h"
29
30 ClassImp(AliLnSpectra)
31
32 AliLnSpectra::AliLnSpectra(const TString& particle, const TString& ptFilename, const TString& tag, const TString& outputFilename, const TString& otag, const Double_t xsec[3])
33 : TObject()
34 , fParticle(particle)
35 , fPtFilename(ptFilename)
36 , fTag(tag)
37 , fOutputFilename(outputFilename)
38 , fOutputTag(otag)
39 , fYMin(-0.5)
40 , fYMax(0.5)
41 , fNormToInel(1)
42 , fIsOnlyGen(0)
43 , fSysErr(1)
44 {
45 //
46 // constructor
47 //
48         for(Int_t i=0; i<3; ++i) fInelXsec[i] = xsec[i];
49 }
50
51 AliLnSpectra::~AliLnSpectra()
52 {
53 //
54 // destructor
55 //
56 }
57
58 Int_t AliLnSpectra::Exec()
59 {
60 //
61 // (invariant) differential yield projection in pt for |y| < 0.5
62 //
63         using namespace std;
64         
65         TFile* finput = new TFile(fPtFilename.Data(), "read");
66         if(finput->IsZombie()) exit(1);
67         
68         // pt and number of events
69         
70         TH1D* hPt = (TH1D*)FindObj(finput, fParticle + "_Pt");
71         TH1D* hStats = (TH1D*)FindObj(finput, fParticle + "_Stats");
72         
73         Double_t nEvent = hStats->Integral(3,3);
74         if(fNormToInel) nEvent = hStats->Integral(4,4);
75         
76         // ouputfile
77         
78         TFile* foutput = new TFile(fOutputFilename.Data(),"recreate");
79         
80         foutput->mkdir(fOutputTag.Data());
81         foutput->cd(fOutputTag.Data());
82         
83         // differential yield
84         
85         TGraphErrors* grDYieldPt = this->GetDiffYieldPt(hPt, nEvent, fParticle + "_DiffYield_Pt");
86         grDYieldPt->Write();
87         
88         TGraphErrors* grSysErrDYieldPt = this->AddSystError(grDYieldPt, fSysErr, fParticle + "_SysErr_DiffYield_Pt");
89         grSysErrDYieldPt->Write();
90         
91         TF1* fncTsallis0 = this->TsallisDiffYield(GetMass(fParticle), fParticle + "_Fit_DiffYield_Pt");
92         fncTsallis0->SetParLimits(0, 0, 1);
93         fncTsallis0->SetParLimits(1, 4, 50);
94         fncTsallis0->SetParLimits(2, 0.01, 10);
95         
96         grSysErrDYieldPt->Fit(fncTsallis0,"RNQ");
97         fncTsallis0->Write();
98         
99         // invariant differential yield
100         
101         TGraphErrors* grInvDYieldPt = this->GetInvDiffYieldPt(grDYieldPt, fParticle + "_InvDiffYield_Pt");
102         grInvDYieldPt->Write();
103         
104         TGraphErrors* grSysErrInvDYieldPt = this->AddSystError(grInvDYieldPt, fSysErr, fParticle + "_SysErr_InvDiffYield_Pt");
105         grSysErrInvDYieldPt->Write();
106         
107         TF1* fncTsallis1 = this->Tsallis(GetMass(fParticle), fParticle + "_Fit_InvDiffYield_Pt");
108         fncTsallis1->SetParLimits(0, 0, 1);
109         fncTsallis1->SetParLimits(1, 4, 50);
110         fncTsallis1->SetParLimits(2, 0.01, 10);
111         
112         grSysErrInvDYieldPt->Fit(fncTsallis1,"RNQ");
113         fncTsallis1->Write();
114         
115         // invariant differential cross section
116         
117         TGraphErrors* grInvXsectPt = GetInvDiffXsectionPt(grInvDYieldPt, fInelXsec, fParticle + "_InvDiffXSection_Pt");
118         grInvXsectPt->Write();
119         
120         TGraphErrors* grSysErrInvXsectPt = GetInvDiffXsectionPt(grSysErrInvDYieldPt, fInelXsec, fParticle + "_SysErr_InvDiffXSection_Pt");
121         grSysErrInvXsectPt->Write();
122         
123         TF1* fncTsallis2 = this->Tsallis(GetMass(fParticle), fInelXsec[0], fParticle + "_Fit_InvDiffXSection_Pt");
124         fncTsallis2->SetParLimits(0, 0, 1);
125         fncTsallis2->SetParLimits(1, 4, 50);
126         fncTsallis2->SetParLimits(2, 0.01, 10);
127         
128         grSysErrInvXsectPt->Fit(fncTsallis2,"RNQ");
129         fncTsallis2->Write();
130         
131         // clean
132         
133         delete fncTsallis0;
134         delete fncTsallis1;
135         delete fncTsallis2;
136         
137         delete grDYieldPt;
138         delete grInvDYieldPt;
139         delete grSysErrDYieldPt;
140         delete grSysErrInvDYieldPt;
141         delete grInvXsectPt;
142         
143         delete foutput;
144         delete finput;
145         
146         return 0;
147 }
148
149 TGraphErrors* AliLnSpectra::GetDiffYieldPt(const TH1D* hPt, Int_t nEvent, const TString& name) const
150 {
151 //
152 // projection of the differential yield in pt
153 // only statistical error of N is taken into account
154 //
155         TGraphErrors* grDYieldPt = new TGraphErrors();
156         grDYieldPt->SetName(name.Data());
157         
158         Double_t dy = fYMax - fYMin;
159         
160         TAxis* xaxis = hPt->GetXaxis();
161         
162         for(Int_t i=1, j=0; i < xaxis->GetNbins()+1; ++i)
163         {
164                 Double_t pt  = xaxis->GetBinCenter(i);
165                 Double_t dpt = xaxis->GetBinWidth(i);
166                 
167                 Double_t n = hPt->GetBinContent(i);
168                 Double_t eN = hPt->GetBinError(i);
169                 
170                 if(n==0) continue;
171                 
172                 Double_t yield = n/(nEvent*dpt*dy);
173                 Double_t yError = yield*eN/n;
174                 
175                 grDYieldPt->SetPoint(j,pt,yield);
176                 grDYieldPt->SetPointError(j++,0,yError);
177         }
178         
179         return grDYieldPt;
180 }
181
182 TGraphErrors* AliLnSpectra::GetInvDiffYieldPt(const TGraphErrors* grDYieldPt, const TString& name) const
183 {
184 //
185 // projection of the invariant differential yield in pt
186 // only statistical error of N is taken into account
187 //
188         TGraphErrors* grInvDYieldPt = new TGraphErrors();
189         grInvDYieldPt->SetName(name.Data());
190         
191         for(Int_t i=0; i < grDYieldPt->GetN(); ++i)
192         {
193                 Double_t pt, y;
194                 grDYieldPt->GetPoint(i,pt,y);
195                 
196                 Double_t errpt = grDYieldPt->GetErrorX(i);
197                 Double_t erry = grDYieldPt->GetErrorY(i);
198                 
199                 y = y/(2.*TMath::Pi()*pt);
200                 erry = erry / (2.*TMath::Pi()*pt);
201                 
202                 grInvDYieldPt->SetPoint(i, pt, y);
203                 grInvDYieldPt->SetPointError(i, errpt, erry);
204         }
205         
206         return grInvDYieldPt;
207 }
208
209 TGraphErrors* AliLnSpectra::GetInvDiffYieldPt(const TH1D* hPt, Int_t nEvent, const TString& name) const
210 {
211 //
212 // projection of the invariant differential yield in pt
213 // only statistical error of N is taken into account
214 //
215         TGraphErrors* grDYieldPt = this->GetDiffYieldPt(hPt, nEvent, "_tmp_diff_yield_");
216         TGraphErrors* grInvDYieldPt = this->GetInvDiffYieldPt(grDYieldPt, name);
217         
218         delete grDYieldPt;
219         return grInvDYieldPt;
220 }
221
222 TGraphErrors* AliLnSpectra::GetInvDiffXsectionPt(const TGraphErrors* grInvDYieldPt, const Double_t* sigma, const TString& name) const
223 {
224 //
225 // invariant differential cross section, only stat. errors
226 // (multiply each value for the total cross section)
227 //
228         Double_t xsec = sigma[0];
229         Double_t errXSec = sigma[1];
230         
231         TGraphErrors* grInvDXsecPt = new TGraphErrors();
232         grInvDXsecPt->SetName(name.Data());
233         
234         for(Int_t i=0; i < grInvDYieldPt->GetN(); ++i)
235         {
236                 Double_t pt, y;
237                 grInvDYieldPt->GetPoint(i,pt,y);
238                 
239                 Double_t errpt = grInvDYieldPt->GetErrorX(i);
240                 Double_t erry = grInvDYieldPt->GetErrorY(i);
241                 
242                 Double_t invxsec = xsec*y;
243                 Double_t err = invxsec*TMath::Sqrt(TMath::Power(errXSec/xsec,2) + TMath::Power(erry/y,2));
244                 
245                 grInvDXsecPt->SetPoint(i, pt, invxsec);
246                 grInvDXsecPt->SetPointError(i, errpt, err);
247         }
248         
249         return grInvDXsecPt;
250 }
251
252 TGraphErrors* AliLnSpectra::AddSystError(const TGraphErrors* gr, Double_t percent, const TString& name) const
253 {
254 //
255 // set the error of h as the given percent of its value
256 //
257         TGraphErrors* grSyst = new TGraphErrors();
258         grSyst->SetName(name.Data());
259         
260         for(Int_t i=0; i < gr->GetN(); ++i)
261         {
262                 Double_t x, y;
263                 gr->GetPoint(i,x,y);
264                 Double_t err = percent*y;
265                 
266                 grSyst->SetPoint(i, x, y);
267                 grSyst->SetPointError(i, 0.025, err);
268         }
269         
270         return grSyst;
271 }
272
273 TF1* AliLnSpectra::Tsallis(Double_t m0, const TString& name, Double_t xmin, Double_t xmax) const
274 {
275 //
276 // Tsallis distribution
277 // Phys. Rev. C 83, 064903 (2011)
278 // Phys. Rev. C 75, 064901 (2007)
279 //
280         TF1* fnc = new TF1(name.Data(), Form("[0]*([1]-1)*([1]-2)*TMath::Power(1+(sqrt(x*x+%f*%f)-%f)/([1]*[2]),-[1])/(2*TMath::Pi()*[1]*[2]*([1]*[2]+%f*([1]-2)))", m0, m0, m0, m0), xmin, xmax);
281         fnc->SetParNames("dN/dy","n","C");
282         fnc->SetParameters(0.1, 7, 0.2);
283         
284         return fnc;
285 }
286
287 TF1* AliLnSpectra::Tsallis(Double_t m0, Double_t xsect, const TString& name, Double_t xmin, Double_t xmax) const
288 {
289 //
290 // Tsallis distribution to fit to invariant cross section
291 // Phys. Rev. C 83, 064903 (2011)
292 // Phys. Rev. C 75, 064901 (2007)
293 //
294         TF1* fnc = new TF1(name.Data(), Form("%f*[0]*([1]-1)*([1]-2)*TMath::Power(1+(sqrt(x*x+%f*%f)-%f)/([1]*[2]),-[1])/(2*TMath::Pi()*[1]*[2]*([1]*[2]+%f*([1]-2)))",xsect, m0, m0, m0, m0), xmin, xmax);
295         fnc->SetParNames("dN/dy","n","C");
296         fnc->SetParameters(0.1, 7, 0.2);
297         
298         return fnc;
299 }
300
301 TF1* AliLnSpectra::TsallisDiffYield(Double_t m0, const TString& name, Double_t xmin, Double_t xmax) const
302 {
303 //
304 // Tsallis distribution to fit differential yield
305 // Phys. Rev. C 83, 064903 (2011)
306 // Phys. Rev. C 75, 064901 (2007)
307 //
308         TF1* fnc = new TF1(name.Data(), Form("x*[0]*([1]-1)*([1]-2)*TMath::Power(1+(sqrt(x*x+%f*%f)-%f)/([1]*[2]),-[1])/([1]*[2]*([1]*[2]+%f*([1]-2)))", m0, m0, m0, m0), xmin, xmax);
309         fnc->SetParNames("dN/dy","n","C");
310         fnc->SetParameters(0.1, 7, 0.2);
311         
312         return fnc;
313 }