New reader for ITS tracking V2
[u/mrichter/AliRoot.git] / HBTAN / AliHBTFunction.cxx
CommitLineData
1b446896 1#include "AliHBTFunction.h"
2/******************************************************************/
3/*
4Piotr Krzysztof Skowronski
5Piotr.Skowronski@cern.ch
6Base classes for HBT functions
7
8 function
9 / \
10 / \
11 / \
12 / \
13 / \
14 / \
15 / \
16 two part four part
17 / | \ / | \
18 / | \ / | \
19 1D 2D 3D 1D 2D 3D
20
21 four particle functions are intendent to be resolution functions:
22 it is mecessary to have simulated particle pair corresponding to given
23 recontructed track pair in order to calculate function simualted value
24 and recontructed value to be further histogrammed
25
26*/
27/******************************************************************/
28/******************************************************************/
29
976183fd 30#include <iostream.h>
1b446896 31ClassImp( AliHBTFunction )
32
33AliHBTFunction::AliHBTFunction()
34{
35
36 fPairCut = new AliHBTEmptyPairCut(); //dummy cut
37}
38
39void AliHBTFunction::
40Write()
41 {
42 if (GetNumerator()) GetNumerator()->Write();
43 if (GetDenominator()) GetDenominator()->Write();
44 TH1* res = GetResult();
976183fd 45 if (res) res->Write();
1b446896 46 }
47/******************************************************************/
48
49TH1* AliHBTFunction::
50GetRatio(Double_t normfactor)
51 {
976183fd 52 //if (gDebug>0)
53 cout<<"Mormfactor is "<<normfactor<<" for "<<fName<<endl;
54
55 if (normfactor == 0.0)
56 {
57 Error("GetRatio","Scaling Factor is 0. Null poiner returned");
58 return 0x0;
59 }
1b446896 60 TString str = fName + " ratio";
61 TH1 *result = (TH1*)GetNumerator()->Clone(str.Data());
62
63 result->SetTitle(str.Data());
976183fd 64 //result->Sumw2();
1b446896 65
66 result->Divide(GetNumerator(),GetDenominator(),normfactor);
67
68 return result;
69
70 }
71/******************************************************************/
72void AliHBTFunction::SetPairCut(AliHBTPairCut* cut)
73{
74//Sets new Pair Cut. Old one is deleted
75//Note that it is created new object instead of simple pointer set
76//I do not want to have pointer
77//to object created somewhere else
78//because in that case I could not believe that
79//it would always exist (sb could delete it)
80//so we have always own copy
81
82 if(!cut)
83 {
84 Error("AliHBTFunction::SetPairCut","argument is NULL");
85 return;
86 }
87 delete fPairCut;
88 fPairCut = (AliHBTPairCut*)cut->Clone();
89
90}
91
92/******************************************************************/
93
94void AliHBTFunction::
95Rename(const Char_t * name)
96 {
97 //renames the function and histograms
98 SetName(name);
99 SetTitle(name);
100
101 TString numstr = fName + " Numerator"; //title and name of the
102 //numerator histogram
103 TString denstr = fName + " Denominator";//title and name of the
104 //denominator histogram
105
106 GetNumerator()->SetName(numstr.Data());
107 GetNumerator()->SetTitle(numstr.Data());
108
109 GetDenominator()->SetName(denstr.Data());
110 GetDenominator()->SetTitle(denstr.Data());
111
112 }
113
114void AliHBTFunction::
115Rename(const Char_t * name, const Char_t * title)
116 {
117 //renames and retitle the function and histograms
118
119 SetName(name);
120 SetTitle(title);
121
122 TString numstrn = fName + " Numerator"; //name of the
123 //numerator histogram
124
125 TString numstrt = fTitle + " Numerator"; //title of the
126 //numerator histogram
127
128 TString denstrn = fName + " Denominator";//name of the
129 //denominator histogram
130
131 TString denstrt = fTitle + " Denominator";//title of the
132 //denominator histogram
133
134
135 GetNumerator()->SetName(numstrn.Data());
136 GetNumerator()->SetTitle(numstrt.Data());
137
138 GetDenominator()->SetName(denstrn.Data());
139 GetDenominator()->SetTitle(denstrt.Data());
140
141
142 }
143
144/******************************************************************/
145/******************************************************************/
146/******************************************************************/
147
148ClassImp( AliHBTTwoPartFctn )
149
150/******************************************************************/
151/******************************************************************/
152/******************************************************************/
153
154ClassImp( AliHBTFourPartFctn)
155
156/******************************************************************/
157/******************************************************************/
158/******************************************************************/
159
160ClassImp( AliHBTTwoPartFctn1D )
161
162AliHBTTwoPartFctn1D::
163AliHBTTwoPartFctn1D(Int_t nbins, Double_t maxXval, Double_t minXval)
164 {
165 //Constructor of Two Part One Dimentional Function
166 // nbins: number of bins in histograms - default 100
167 // maxXval and minXval: range of histgram(s) default 0 - 0.15 (GeV)
168
169
170 TString numstr = fName + " Numerator"; //title and name of the
171 //numerator histogram
172 TString denstr = fName + " Denominator";//title and name of the
173 //denominator histogram
174
175 fNumerator = new TH1D(numstr.Data(),numstr.Data(),nbins,minXval,maxXval);
176 fDenominator = new TH1D(denstr.Data(),denstr.Data(),nbins,minXval,maxXval);
177
178 fNumerator->Sumw2();
179 fDenominator->Sumw2();
180
976183fd 181 fNBinsToScale = 30;
182
1b446896 183 }
184/******************************************************************/
185AliHBTTwoPartFctn1D::~AliHBTTwoPartFctn1D()
186{
187 delete fNumerator;
188 delete fDenominator;
189}
976183fd 190/******************************************************************/
191
192void AliHBTTwoPartFctn1D::ProcessSameEventParticles(AliHBTPair* pair)
193{
194 //Fills the numerator
195 pair = CheckPair(pair);
196 if(pair) fNumerator->Fill(GetValue(pair));
197}
198/******************************************************************/
199void AliHBTTwoPartFctn1D::ProcessDiffEventParticles(AliHBTPair* pair)
200 {
201 //fills denumerator
202 pair = CheckPair(pair);
203 if(pair) fDenominator->Fill(GetValue(pair));
204
205 }
206/******************************************************************/
207Double_t AliHBTTwoPartFctn1D::Scale()
208{
209 if(!fNumerator)
210 {
211 Error("Scale","No numerator");
212 return 0.0;
213 }
214 if(!fDenominator)
215 {
216 Error("Scale","No denominator");
217 return 0.0;
218 }
219
220 if(fNBinsToScale < 1)
221 {
222 return 0.0;
223 Error("Scale","Number of bins for scaling is smaller thnan 1");
224 }
225 Int_t nbins = fNumerator->GetNbinsX();
226 if (fNBinsToScale > nbins)
227 {
228 Error("Scale","Number of bins for scaling is bigger thnan number of bins in histograms");
229 return 0.0;
230 }
231 Double_t ratios[fNBinsToScale];
232
233 Int_t offset = nbins - fNBinsToScale - 1;
234 Int_t i;
235 for ( i = offset; i< nbins; i++)
236 {
237 if ( fNumerator->GetBinContent(i) == 0.0 )
238 {
239 ratios[i - offset] = -1.0; //since we play with histograms negative is impossible
240 //so it is good flag
241 }
242 else
243 {
244 ratios[i - offset] = fDenominator->GetBinContent(i)/fNumerator->GetBinContent(i);
245 }
246 }
247
248 Double_t sum = 0;
249 Int_t skipped = 0;
250 for (i = 0; i<fNBinsToScale; i++)
251 {
252 if (ratios[i] == -1.0) skipped++;
253 else sum += ratios[i];
254 }
255 cout<<"sum="<<sum<<" fNBinsToScale="<<fNBinsToScale<<" skipped="<<skipped<<endl;
256
257 return sum/(Double_t)(fNBinsToScale - skipped);
258}
259
1b446896 260/******************************************************************/
261/******************************************************************/
262/******************************************************************/
263
264ClassImp( AliHBTTwoPartFctn2D )
265
266AliHBTTwoPartFctn2D::
267AliHBTTwoPartFctn2D(Int_t nXbins, Double_t maxXval, Double_t minXval ,
268 Int_t nYbins, Double_t maxYval, Double_t minYval)
269
270{
271 TString numstr = fName + " Numerator"; //title and name of the
272 //numerator histogram
273 TString denstr = fName + " Denominator";//title and name of the
274 //denominator histogram
275
276 fNumerator = new TH2D(numstr.Data(),numstr.Data(),
277 nXbins,minXval,maxXval,
278 nYbins,minYval,maxYval);
279
280 fDenominator = new TH2D(denstr.Data(),denstr.Data(),
281 nXbins,minXval,maxXval,
282 nYbins,minYval,maxYval);
283
284 fNumerator->Sumw2();
285 fDenominator->Sumw2();
286
287}
288AliHBTTwoPartFctn2D::~AliHBTTwoPartFctn2D()
289{
290 delete fNumerator;
291 delete fDenominator;
292}
293void AliHBTTwoPartFctn2D::ProcessSameEventParticles(AliHBTPair* pair)
294{
295 pair = CheckPair(pair);
296 if(pair)
297 {
298 Double_t x,y;
299 GetValues(pair,x,y);
300 fNumerator->Fill(y,x);
301 }
302}
303
304void AliHBTTwoPartFctn2D::ProcessDiffEventParticles(AliHBTPair* pair)
305{
306 pair = CheckPair(pair);
307 if(pair)
308 {
309 Double_t x,y;
310 GetValues(pair,x,y);
311 fDenominator->Fill(y,x);
312 }
313
314}
315
316
317/******************************************************************/
318/******************************************************************/
319/******************************************************************/
320
321ClassImp( AliHBTTwoPartFctn3D)
322
323AliHBTTwoPartFctn3D::
324AliHBTTwoPartFctn3D(Int_t nXbins, Double_t maxXval, Double_t minXval,
325 Int_t nYbins, Double_t maxYval, Double_t minYval,
326 Int_t nZbins, Double_t maxZval, Double_t minZval)
327
328{
329 TString numstr = fName + " Numerator"; //title and name of the
330 //numerator histogram
331 TString denstr = fName + " Denominator";//title and name of the
332 //denominator histogram
333
334 fNumerator = new TH3D(numstr.Data(),numstr.Data(),
335 nXbins,minXval,maxXval,
336 nYbins,minYval,maxYval,
337 nZbins,minZval,maxZval);
338
339 fDenominator = new TH3D(denstr.Data(),denstr.Data(),
340 nXbins,minXval,maxXval,
341 nYbins,minYval,maxYval,
342 nZbins,minZval,maxZval);
343
344 fNumerator->Sumw2();
345 fDenominator->Sumw2();
346
347}
348
349
350AliHBTTwoPartFctn3D::~AliHBTTwoPartFctn3D()
351{
352 delete fNumerator;
353 delete fDenominator;
354}
355
356
357/******************************************************************/
358/******************************************************************/
359/******************************************************************/
360ClassImp( AliHBTFourPartFctn2D)
361
362
363AliHBTFourPartFctn2D::
364AliHBTFourPartFctn2D(Int_t nXbins, Double_t maxXval, Double_t minXval ,
365 Int_t nYbins, Double_t maxYval, Double_t minYval)
366
367{
368 TString numstr = fName + " Numerator"; //title and name of the
369 //numerator histogram
370 TString denstr = fName + " Denominator";//title and name of the
371 //denominator histogram
372
373 fNumerator = new TH2D(numstr.Data(),numstr.Data(),
374 nXbins,minXval,maxXval,
375 nYbins,minYval,maxYval);
376
377 fDenominator = new TH2D(denstr.Data(),denstr.Data(),
378 nXbins,minXval,maxXval,
379 nYbins,minYval,maxYval);
380
381 fNumerator->Sumw2();
382 fDenominator->Sumw2();
383
384}
385AliHBTFourPartFctn2D::~AliHBTFourPartFctn2D()
386{
387 delete fNumerator;
388 delete fDenominator;
389}
390void AliHBTFourPartFctn2D::
391ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
392{
393 partpair = CheckPair(partpair);
394 trackpair = CheckPair(trackpair);
395 if( partpair && trackpair)
396 {
397 Double_t x,y;
398 GetValues(trackpair,partpair,x,y);
399 fNumerator->Fill(y,x);
400 }
401}
402
403void AliHBTFourPartFctn2D::
404ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
405{
406 partpair = CheckPair(partpair);
407 trackpair = CheckPair(trackpair);
408 if( partpair && trackpair)
409 {
410 Double_t x,y;
411 GetValues(trackpair,partpair,x,y);
412 fDenominator->Fill(y,x);
413 }
414
415}
416