]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/EBYE/LRC/AliLRCAnalysis.cxx
CTP raw data does not contain CDH, so skip the warning about the wrong data size...
[u/mrichter/AliRoot.git] / PWG2 / EBYE / LRC / AliLRCAnalysis.cxx
CommitLineData
d03e5de4 1/**************************************************************************
2 * Author: Panos Christakoglou. *
3 * Contributors are mentioned in the code where appropriate. *
4 * *
5 * Permission to use, copy, modify and distribute this software and its *
6 * documentation strictly for non-commercial purposes is hereby granted *
7 * without fee, provided that the above copyright notice appears in all *
8 * copies and that both the copyright notice and this permission notice *
9 * appear in the supporting documentation. The authors make no claims *
10 * about the suitability of this software for any purpose. It is *
11 * provided "as is" without express or implied warranty. *
12 **************************************************************************/
13//-------------------------------------------------------------------------
14// Description:
15// This class include into LRC library for Long-Range Correlation analysis
16// it is base class for NN, PtN, PtPt
17// implements base methods for thees classes
18// Origin: Petr Naumenko, SPbSU-CERN, Petr.Naoumenko@cern.ch
19//-------------------------------------------------------------------------
20
21/* $Id$ */
22
23//-------------------------------------------------------------------------
24// LRC library for Long-Range Correlation analysis
25//
26// Origin: Petr Naumenko, SPbSU-CERN, Petr.Naoumenko@cern.ch
27//-------------------------------------------------------------------------
28
29#include "AliLRCAnalysis.h"
30
31ClassImp(AliLRCAnalysis)
32
33/******************************************************
34 * AliLRCAnalysis class
35 ******************************************************/
36
37AliLRCAnalysis::AliLRCAnalysis(): fPrAbs(new TH1D()), fPrRel(new TH1D()), fPrf(new TH1D()), fPrb(new TH1D()), fileHist(new TFile()), fdptb(.0), fEntries(0), fSx((char*)" "), fSy((char*)" "){
38//Empty constructor
39}
40
41AliLRCAnalysis::AliLRCAnalysis(const AliLRCAnalysis& a):fPrAbs(a.fPrAbs), fPrRel(a.fPrRel), fPrf(a.fPrf), fPrb(a.fPrb), fileHist(a.fileHist), fdptb(a.fdptb), fEntries(a.fEntries), fSx(a.fSx), fSy(a.fSy){
42//Constructor
43}
44
45AliLRCAnalysis& AliLRCAnalysis::operator= (const AliLRCAnalysis& a){
46//Operator =
47 if(this!=&a){
48 fPrAbs = a.fPrAbs;
49 fPrRel = a.fPrRel;
50 fPrf = a.fPrf;
51 fPrb = a.fPrb;
52 fileHist = a.fileHist;
53 fSx = a.fSx;
54 fSy = a.fSy;
55 fdptb = a.fdptb;
56 fEntries = a.fEntries;
57 }
58 return *this;
59}
60
61
62AliLRCAnalysis::~AliLRCAnalysis() {
63//Destructor
64 delete fPrAbs;
65 delete fPrRel;
66 delete fileHist;
67}
68
69double AliLRCAnalysis::HI2(TH1D *h, double a, double b, double xmin, double xmax) const {
70//hi square calculation of approximation 1d hist with ax+b between xmin and xmax
71 double trueN = 0;
72 double fhi2=0;
73 double num;
74 double fN=h->GetNbinsX();
75 double f = h->GetXaxis()->GetXmin();
76 double fdf = (h->GetXaxis()->GetXmax()-f)/fN;
77 int fNmin = int((xmin-f)/fdf)+1;
78 int fNmax = int((xmax-f)/fdf)+1;
79 for(int i=fNmin; i<=fNmax; i++) {
80 double fw = h->GetBinError(i);
81 if(fw!=0){
82 num = b*(i*fdf-fdf/2+f)+a-h->GetBinContent(i);
83 fhi2 = fhi2 + (num*num) / (fw*fw);
84 trueN++;
85 }
86 }
87 return fhi2/(trueN-2);
88}
89
90double AliLRCAnalysis::HI2(TH1D *h, double a, double b) const {
91//hi square calculation of approximation 1d hist with ax+b
92 double trueN = 0;
93 double fhi2=0;
94 double num;
95 int fN=h->GetNbinsX();
96 double f = h->GetXaxis()->GetXmin();
97 double fdf = (h->GetXaxis()->GetXmax()-f)/fN;
98 int fNmin = 1;
99 int fNmax = fN;
100 for(int i=fNmin; i<=fNmax; i++) {
101 double fw = h->GetBinError(i);
102 if(fw!=0){
103 num = b*(i*fdf-fdf/2+f)+a-h->GetBinContent(i);
104 fhi2 = fhi2 + (num*num) / (fw*fw);
105 trueN++;
106 }
107 }
108 return fhi2/(trueN-2);
109}
110
111
112void AliLRCAnalysis::CreateHist(char *name, char *nameAbs, char *nameRel, char *atitleF, char *atitleB, char *rtitleF, char *rtitleB, TH2D* sourceHist) {
113//Create absolute and relation var histogramm
114 TProfile* profX = (TProfile*) sourceHist->ProfileX(name, 1, sourceHist->GetNbinsY());
115 fEntries = (int) sourceHist->GetEntries();
116 fPrAbs = new TH1D(nameAbs, profX->GetTitle(), profX->GetXaxis()->GetNbins(), profX->GetXaxis()->GetXmin(), profX->GetXaxis()->GetXmax());
117 fPrAbs->SetOption("E");
118 fPrRel = new TH1D(nameRel, profX->GetTitle(), profX->GetXaxis()->GetNbins(), profX->GetXaxis()->GetXmin()/sourceHist->ProjectionX()->GetMean(), profX->GetXaxis()->GetXmax()/sourceHist->ProjectionX()->GetMean());
119 fPrRel->SetOption("E");
120 fPrAbs->GetXaxis()->SetTitle(atitleF);
121 fPrAbs->GetYaxis()->SetTitle(atitleB);
122 fPrRel->GetXaxis()->SetTitle(rtitleF);
123 fPrRel->GetYaxis()->SetTitle(rtitleB);
124 fPrf = (TH1D*) sourceHist->ProjectionX();
125 fPrb = (TH1D*) sourceHist->ProjectionY();
126 fPrf->GetXaxis()->SetTitle(atitleF);
127 fPrf->GetYaxis()->SetTitle("Tracks");
128 fPrb->GetXaxis()->SetTitle(atitleB);
129 fPrb->GetYaxis()->SetTitle("Tracks");
130 fSx = atitleF;
131 fSy = atitleB;
132}
133
134void AliLRCAnalysis::DrawAbs() {
135// Draw abs var histrogramm
136 double mnf;
137 double y1, y2, x1, x2;
138 Int_t i, n;
139 char str[50];
140 mnf = fPrf->GetMean();
141 fPrAbs->SetStats(0);
142 fPrAbs->Fit("pol1", "", "", mnf-2*sqrt(mnf), mnf+2*sqrt(mnf));
143 //fPrAbs->Fit("pol1");
144 TF1 *fit1 = fPrAbs->GetFunction("pol1");
145 y1=fPrAbs->GetBinContent(1)-fPrAbs->GetBinError(1);
146 y2=fPrAbs->GetBinContent(1)+fPrAbs->GetBinError(1);
147
148 n=fPrAbs->GetNbinsX();
149 for(i=2; i<=n; i++){
150 if(fPrAbs->GetBinContent(i)-fPrAbs->GetBinError(i)<y1)
151 y1=fPrAbs->GetBinContent(i)-fPrAbs->GetBinError(i);
152 if(fPrAbs->GetBinContent(i)+fPrAbs->GetBinError(i)>y2)
153 y2=fPrAbs->GetBinContent(i)+fPrAbs->GetBinError(i);
154 }
155 x1 = fPrAbs->GetXaxis()->GetXmin();
156 x2 = fPrAbs->GetXaxis()->GetXmax();
157 TPaveText *pt1 = new TPaveText(x1+(x2-x1)/4, y1+(y2-y1)/2, x2-(x2-x1)/4, y2);
158 pt1->SetTextSize(0.03);
159
160 sprintf(str, "Entries = %i", fEntries);
161 pt1->AddText(str);
162 sprintf(str, "a = %f #pm %f", fit1->GetParameter(0), fit1->GetParError(0));
163 pt1->AddText(str);
164 sprintf(str, "b = %f #pm %f", fit1->GetParameter(1), fit1->GetParError(1));
165 pt1->AddText(str);
166 sprintf(str, "#hat{#chi}^{2} = #chi^{2}/(n-2) = %f", HI2(fPrAbs, fit1->GetParameter(0), fit1->GetParameter(1), mnf-2*sqrt(mnf), mnf+2*sqrt(mnf)));
167 pt1->AddText(str);
168 sprintf(str, "<%s> = %f " , fSx, fPrf->GetMean());
169 pt1->AddText(str);
170 sprintf(str, "<%s> = %f", fSy, fPrb->GetMean());
171 pt1->AddText(str);
172
173 sprintf(str, "<<%s>> = %f " , fSx, fPrf->GetRMS());
174 pt1->AddText(str);
175 sprintf(str, "<<%s>> = %f", fSy, fPrb->GetRMS());
176 pt1->AddText(str);
177 if(fdptb){
178 sprintf(str, "d%s = %f", fSy, fdptb);
179 pt1->AddText(str);
180 }
181
182
183 pt1->SetTextAlign(12);
184 pt1->SetTextFont(42);
185 pt1->SetFillColor(0);
186
187 pt1->Draw();
188}
189
190void AliLRCAnalysis::DrawRel() {
191// Draw rel var histogramm
192double mnf;
193double y1, y2, x1, x2;
194Int_t i, n;
195char str[50];
196
197mnf = fPrf->GetMean();
198
199fPrRel->SetStats(0);
200AliLRCFit *fit1 = new AliLRCFit(fPrRel,1-2/sqrt(mnf), 1+2/sqrt(mnf));
201TF1 *f1 = new TF1("f1", "[0] + [1]*x", 0, fPrRel->GetXaxis()->GetXmax());
202f1->SetParameter(0,fit1->Geta());
203f1->SetParameter(1,fit1->Getb());
204fPrRel->Fit("f1", "", "", 1-2/sqrt(mnf), 1+2/sqrt(mnf));
205 y1=fPrRel->GetBinContent(1)-fPrRel->GetBinError(1);
206 y2=fPrRel->GetBinContent(1)+fPrRel->GetBinError(1);
207 n=fPrRel->GetNbinsX();
208 for(i=2; i<=n; i++) {
209 if(fPrRel->GetBinContent(i)-fPrRel->GetBinError(i)<y1)
210 y1=fPrRel->GetBinContent(i)-fPrRel->GetBinError(i);
211 if(fPrRel->GetBinContent(i)+fPrRel->GetBinError(i)>y2)
212 y2=fPrRel->GetBinContent(i)+fPrRel->GetBinError(i);
213 }
214 x1 = fPrRel->GetXaxis()->GetXmin();
215 x2 = fPrRel->GetXaxis()->GetXmax();
216 TPaveText *pt1 = new TPaveText(x1+(x2-x1)/4, y1+(y2-y1)/2, x2-(x2-x1)/4, y2);
217 pt1->SetTextSize(0.03);
218 sprintf(str, "Entries = %i", fEntries);
219 pt1->AddText(str);
220 sprintf(str, "a = %f #pm %f", fit1->Geta(), fit1->Getda());
221 pt1->AddText(str);
222 sprintf(str, "b = %f #pm %f", fit1->Getb(), fit1->Getdb());
223 pt1->AddText(str);
224 sprintf(str, "#hat{#chi}^{2} = #chi^{2}/(n-2) = %f", fit1->Gethi2());
225 pt1->AddText(str);
226 sprintf(str, "<%s> = %f " , fSx, fPrf->GetMean());
227 pt1->AddText(str);
228 sprintf(str, "<%s> = %f", fSy, fPrb->GetMean());
229 pt1->AddText(str);
230 sprintf(str, "<<%s>> = %f " , fSx, fPrf->GetRMS());
231 pt1->AddText(str);
232 sprintf(str, "<<%s>> = %f", fSy, fPrb->GetRMS());
233 pt1->AddText(str);
234 if(fdptb){
235 sprintf(str, "d%s = %f", fSy, fdptb);
236 pt1->AddText(str);
237 }
238 pt1->SetTextAlign(12);
239 pt1->SetTextFont(42);
240 pt1->SetFillColor(0);
241 pt1->Draw();
242}
243
244void AliLRCAnalysis::SetGraphics() const {
245// Set root graph style
246 gStyle->SetCanvasColor(10);
247 gStyle->SetFrameFillColor(10);
248 gStyle->SetStatColor(0);
249 gStyle->SetPadColor(0);
250}
251
252double AliLRCAnalysis::Integral(TH2D* source, Int_t nbin) const {
253// calculate the integrall for x bin and y bins of 2d histogramm
254 double sum = 0;
255 for(Int_t i = 1; i<=source->GetNbinsY(); i++) {
256 sum += source->GetYaxis()->GetBinCenter(i)*source->GetBinContent(nbin, i);
257 }
258 return sum;
259}
260
261void AliLRCAnalysis::SetErrors(TH2D* source, const char *name){
262// Calculate errors for NN
263 TProfile* profX = (TProfile*) source->ProfileX(name, 1, source->GetNbinsY());
264 for(int i = 0; i < profX->GetNbinsX(); i++)
265 {
266 fPrAbs->SetBinContent(i, profX->GetBinContent(i));
267 if(fPrf->GetBinContent(i)!=0)
268 {
269 fPrAbs->SetBinError(i,sqrt(profX->GetBinContent(i)/fPrf->GetBinContent(i)));
270 }
271 fPrRel->SetBinContent(i, fPrAbs->GetBinContent(i)/fPrb->GetMean());
272 fPrRel->SetBinError(i,fPrAbs->GetBinError(i)/fPrb->GetMean());
273 }
274}
275
276void AliLRCAnalysis::SetErrors(TH2D* source, const char *name, double ptd, TH2D* nb){
277//Calculate arrors for ptn and ptpt
278 TProfile* profX = (TProfile*) source->ProfileX(name, 1, source->GetNbinsY());
279 fdptb = ptd;
280 double pt;
281 for(int i = 0; i < profX->GetNbinsX(); i++)
282 {
283 fPrAbs->SetBinContent(i, profX->GetBinContent(i));
284 if(fPrf->GetBinContent(i)!=0)
285 {
286 pt = profX->GetBinContent(i);
287 fPrAbs->SetBinError(i,ptd*sqrt(Integral(nb,i))/fPrf->GetBinContent(i));
288 }
289 fPrRel->SetBinContent(i, fPrAbs->GetBinContent(i)/fPrb->GetMean());
290 fPrRel->SetBinError(i,fPrAbs->GetBinError(i)/fPrb->GetMean());
291 }
292
293}
294
295