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