30266ff99d8102b438ed2716beaa80afffc84652
[u/mrichter/AliRoot.git] / HBTAN / AliHBTFunction.cxx
1 #include "AliHBTFunction.h"
2
3 /* $Id: */
4
5 //--------------------------------------------------------------------
6 //AliHBTFunction
7 //Author: Piotr Krzysztof Skowronski
8 //Piotr.Skowronski@cern.ch
9 //Base classes for HBT functions
10 /*
11  OnePairFctn       Function                  TwoPairFctn
12    | \    \        /  |   \                      /|\   
13    |  \    \      /   |    \                    / | \    
14    |   \    \   1D   2D    3D                  /  |  \    
15    |    \    \  / \   |\   | \________________/__ |   \             
16    |     \    \/   \  | \__|_____________    /   \|    \             
17    |      \    \    \_|____|__           \  /     |     \            
18    |       \ /  \___  |    |  \           \/      |\     \              
19    |        /       \ |    |   \          /\      | \     \             
20    |       / \       \|    |    \        /  \     |  \     \              
21    |      /   \      /\    |     \      /    \    |   \     |              
22    |     /     \    /  \   |      \    /      \   |    \    |               
23    |    /       \  /    \  |       \  /        \  |     \   |
24    |   /         \/      \ |        \/          \ |      \  |               
25  OnePair1D   OnePair2D  OnePair3D  TwoPair1D  TwoPair2D TwoPair3D
26
27
28  four particle functions are intendent to be resolution functions:
29  it is mecessary to have simulated particle pair corresponding to given
30  recontructed track pair in order to calculate function simualted value 
31  and recontructed value, to be further histogrammed
32 */
33 //--------------------------------------------------------------------- 
34
35 #include <TH2.h>
36 #include <TH3.h>
37
38 /******************************************************************/
39 /******************************************************************/
40
41 ClassImp( AliHBTFunction )
42
43 AliHBTFunction::AliHBTFunction():
44   fPairCut(new AliHBTEmptyPairCut())   //dummy cut
45 {
46 //Default constructor
47 }
48 /******************************************************************/
49 AliHBTFunction::AliHBTFunction(const char* name,const char* title):
50   TNamed(name,title),
51   fPairCut(new AliHBTEmptyPairCut()) //dummy cut  
52 {
53 //Constructor  
54 }
55 /******************************************************************/
56
57 AliHBTFunction::~AliHBTFunction()
58 {
59 //destructor  
60   delete fPairCut;
61 }
62 /******************************************************************/
63
64 void AliHBTFunction::WriteFunction()
65 {
66 //writes result of the function to file
67    if (GetNumerator()) GetNumerator()->Write();
68    if (GetDenominator()) GetDenominator()->Write();
69    TH1* res = GetResult();
70    if (res) res->Write();
71 }
72 /******************************************************************/
73
74 TH1* AliHBTFunction::GetRatio(Double_t normfactor)
75  {
76  //returns ratio of numerator and denominator
77  //
78    if (gDebug>0) Info("GetRatio","Norm. Factor is %f for %s",normfactor,GetName());
79    
80    if (normfactor == 0.0)
81     {
82       Error("GetRatio","Scaling Factor is 0. Null poiner returned");
83       return 0x0;
84     }
85    TString str = fName + " ratio";
86    TH1 *result = (TH1*)GetNumerator()->Clone(str.Data());
87    
88    result->SetTitle(str.Data());
89    
90    result->Divide(GetNumerator(),GetDenominator(),normfactor);
91    
92    return result;
93    
94  }
95 /******************************************************************/
96 void AliHBTFunction::SetPairCut(AliHBTPairCut* cut)
97 {
98 //Sets new Pair Cut. Old one is deleted
99 //Note that it is created new object instead of simple pointer set
100 //I do not want to have pointer 
101 //to object created somewhere else
102 //because in that case I could not believe that 
103 //it would always exist (sb could delete it)
104 //so we have always own copy
105
106  if(!cut) 
107    {
108      Error("AliHBTFunction::SetPairCut","argument is NULL");
109      return;
110    }
111  delete fPairCut;
112  fPairCut = (AliHBTPairCut*)cut->Clone();
113  
114 }
115
116 /******************************************************************/
117
118 void AliHBTFunction::Rename(const Char_t * name)
119  {
120  //renames the function and histograms
121   SetName(name);
122   SetTitle(name);
123   
124   TString numstr = fName + " Numerator";  //title and name of the 
125                                            //numerator histogram
126   TString denstr = fName + " Denominator";//title and name of the 
127                                            //denominator histogram
128  
129   GetNumerator()->SetName(numstr.Data());
130   GetNumerator()->SetTitle(numstr.Data());
131   
132   GetDenominator()->SetName(denstr.Data());
133   GetDenominator()->SetTitle(denstr.Data());
134   
135  }
136
137 void AliHBTFunction::Rename(const Char_t * name, const Char_t * title)
138  {
139  //renames and retitle the function and histograms
140  
141   SetName(name);
142   SetTitle(title);
143   
144   TString numstrn = fName + " Numerator";  //name of the 
145                                            //numerator histogram
146
147   TString numstrt = fTitle + " Numerator";  //title of the 
148                                            //numerator histogram
149                    
150   TString denstrn = fName + " Denominator";//name of the 
151                                            //denominator histogram
152
153   TString denstrt = fTitle + " Denominator";//title of the 
154                                            //denominator histogram
155                    
156  
157   GetNumerator()->SetName(numstrn.Data());
158   GetNumerator()->SetTitle(numstrt.Data());
159   
160   GetDenominator()->SetName(denstrn.Data());
161   GetDenominator()->SetTitle(denstrt.Data());
162
163
164  }
165 /******************************************************************/
166
167 void AliHBTFunction::InitFunction()
168 {
169 //Iniotializes fctn.: Resets histograms
170 //In case histograms are not created in ctor, builds with default parameters
171   if ( !(GetNumerator()&&GetDenominator()) ) BuildHistos();
172   GetNumerator()->Reset();
173   GetDenominator()->Reset();
174 }
175
176 /******************************************************************/
177 /******************************************************************/
178 /******************************************************************/
179
180 ClassImp( AliHBTOnePairFctn )
181
182 /******************************************************************/
183 /******************************************************************/
184 /******************************************************************/
185
186 ClassImp( AliHBTTwoPairFctn)
187
188 /******************************************************************/
189 /******************************************************************/
190 /******************************************************************/
191
192 ClassImp( AliHBTFunction1D )
193
194 const Int_t AliHBTFunction1D::fgkDefaultNBins = 100;
195 const Float_t AliHBTFunction1D::fgkDefaultMin = 0.0;
196 const Float_t AliHBTFunction1D::fgkDefaultMax = 0.15;
197 const UInt_t AliHBTFunction1D::fgkDefaultNBinsToScale = 30;
198
199
200 AliHBTFunction1D::AliHBTFunction1D():
201  fNumerator(0x0),
202  fDenominator(0x0),
203  fNBinsToScale(fgkDefaultNBinsToScale)
204 {//default constructor
205 }
206 /******************************************************************/
207
208 AliHBTFunction1D::AliHBTFunction1D(Int_t nbins, Float_t maxXval, Float_t minXval):
209  fNumerator(0x0),
210  fDenominator(0x0),
211  fNBinsToScale(fgkDefaultNBinsToScale)
212 {
213  //Constructor of Two Part One Dimentional Function 
214  // nbins: number of bins in histograms - default 100
215  // maxXval and minXval: range of histgram(s) default 0 - 0.15 (GeV)
216  BuildHistos(nbins,maxXval,minXval);
217 }
218 /******************************************************************/
219 AliHBTFunction1D::AliHBTFunction1D(const Char_t *name, const Char_t *title):
220  AliHBTFunction(name,title),
221  fNumerator(0x0),
222  fDenominator(0x0),
223  fNBinsToScale(fgkDefaultNBinsToScale)
224 {//constructor
225 }
226 /******************************************************************/
227 AliHBTFunction1D::AliHBTFunction1D(const Char_t *name, const Char_t *title,
228                                    Int_t nbins, Float_t maxXval, Float_t minXval):
229  AliHBTFunction(name,title),
230  fNumerator(0x0),
231  fDenominator(0x0),
232  fNBinsToScale(fgkDefaultNBinsToScale)
233 {
234 //constructor
235   BuildHistos(nbins,maxXval,minXval);
236 }
237 /******************************************************************/
238
239 AliHBTFunction1D::~AliHBTFunction1D()
240 {
241 //destructor
242   delete fNumerator;
243   delete fDenominator;
244 }
245 /******************************************************************/
246 void AliHBTFunction1D::BuildHistos()
247 {
248 //builds histograms with default settings
249  BuildHistos(fgkDefaultNBins,fgkDefaultMax,fgkDefaultMin);
250 }
251
252 /******************************************************************/
253
254 void AliHBTFunction1D::BuildHistos(Int_t nbins, Float_t max, Float_t min)
255 {
256 //builds numarator and denominator hitograms
257   TString numstr = fName + " Numerator";  //title and name of the 
258                                           //numerator histogram
259   TString denstr = fName + " Denominator";//title and name of the 
260                                           //denominator histogram
261   fNumerator   = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
262   fDenominator = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
263    
264   fNumerator->Sumw2();
265   fDenominator->Sumw2();
266 }
267 /******************************************************************/
268
269 Double_t AliHBTFunction1D::Scale()
270 {
271  //Calculates the factor that should be used to scale 
272  //quatience of fNumerator and fDenominator to 1 at tail
273  return Scale(fNumerator,fDenominator);
274 }
275 /******************************************************************/
276
277 Double_t AliHBTFunction1D::Scale(TH1D* num,TH1D* den)
278 {
279  //Calculates the factor that should be used to scale 
280  //quatience of num and den to 1 at tail
281  
282   if (gDebug>0) Info("Scale","Enetered Scale()");
283   if(!num) 
284    {
285      Error("Scale","No numerator");
286      return 0.0;
287    }
288   if(!den) 
289    {
290      Error("Scale","No denominator");
291      return 0.0;
292    }
293   
294   if(fNBinsToScale < 1) 
295    {
296     return 0.0;
297     Error("Scale","Number of bins for scaling is smaller thnan 1");
298    }
299   UInt_t nbins = num->GetNbinsX();
300   if (fNBinsToScale > nbins) 
301    {
302     Error("Scale","Number of bins for scaling is bigger thnan number of bins in histograms");
303     return 0.0;
304    }
305   if (gDebug>0) Info("Scale","No errors detected");
306
307   Double_t ratio;
308   Double_t sum = 0;
309   Int_t n = 0;
310   
311   Int_t offset = nbins - fNBinsToScale - 1; 
312
313   for (UInt_t i = offset; i< nbins; i++)
314    {
315     if ( num->GetBinContent(i) > 0.0 )
316      {
317        ratio = den->GetBinContent(i)/num->GetBinContent(i);
318        sum += ratio;
319        n++;
320      }
321    }
322   
323   if(gDebug > 0) Info("Scale","sum=%f fNBinsToScale=%d n=%d",sum,fNBinsToScale,n);
324   
325   if (n == 0) return 0.0;
326   Double_t ret = sum/((Double_t)n);
327
328   if(gDebug > 0) Info("Scale","returning %f",ret);
329   return ret;
330
331
332 /******************************************************************/
333 /******************************************************************/
334 /******************************************************************/
335
336 //____________________
337 ///////////////////////////////////////////////////////
338 //                                                   //
339 // AliHBTFunction2D                                  //
340 //                                                   //
341 // Base Calss for 2-dimensinal Functions             //
342 //                                                   //
343 // Piotr.Skowronski@cern.ch                          //
344 // http://alisoft.cern.ch/people/skowron/analyzer    //
345 //                                                   //
346 ///////////////////////////////////////////////////////
347
348 ClassImp( AliHBTFunction2D )
349
350 const Int_t AliHBTFunction2D::fgkDefaultNBinsX = 200;//default number of Bins in X axis in histograms
351 const Float_t AliHBTFunction2D::fgkDefaultMinX = 0.0;//Default min value of X axis in histograms
352 const Float_t AliHBTFunction2D::fgkDefaultMaxX = 1.5;//Default max value of X axis inhistograms
353
354 const Int_t AliHBTFunction2D::fgkDefaultNBinsY = 200;//default number of Bins in histograms
355 const Float_t AliHBTFunction2D::fgkDefaultMinY = -0.15;//Default min value of histograms
356 const Float_t AliHBTFunction2D::fgkDefaultMaxY =  0.15;//Default max value of histograms
357
358 const UInt_t AliHBTFunction2D::fgkDefaultNBinsToScaleX = 30;//Default number of X bins used for scaling to tale
359 const UInt_t AliHBTFunction2D::fgkDefaultNBinsToScaleY = 30;//Default number of bins used for scaling to tale
360
361 /******************************************************************/
362 AliHBTFunction2D::AliHBTFunction2D():
363  fNumerator(0x0),
364  fDenominator(0x0),
365  fNBinsToScaleX(fgkDefaultNBinsToScaleX), 
366  fNBinsToScaleY(fgkDefaultNBinsToScaleY)
367 {//default constructor
368 }
369 /******************************************************************/
370 AliHBTFunction2D::AliHBTFunction2D(const Char_t *name, const Char_t *title):
371  AliHBTFunction(name,title),
372  fNumerator(0x0),
373  fDenominator(0x0),
374  fNBinsToScaleX(fgkDefaultNBinsToScaleX), 
375  fNBinsToScaleY(fgkDefaultNBinsToScaleY)
376 {//constructor
377 }
378 /******************************************************************/
379
380 AliHBTFunction2D::AliHBTFunction2D(Int_t nXbins, Double_t maxXval, Double_t minXval,
381                                    Int_t nYbins, Double_t maxYval, Double_t minYval):
382  fNumerator(0x0),
383  fDenominator(0x0),
384  fNBinsToScaleX(fgkDefaultNBinsToScaleX), 
385  fNBinsToScaleY(fgkDefaultNBinsToScaleY)
386 {
387   BuildHistos(nXbins,maxXval,minXval,nYbins,maxYval,minYval);
388 }         
389 /******************************************************************/
390
391 AliHBTFunction2D::AliHBTFunction2D(const Char_t *name, const Char_t *title,
392                                    Int_t nXbins, Double_t maxXval, Double_t minXval,
393                                    Int_t nYbins, Double_t maxYval, Double_t minYval):
394  AliHBTFunction(name,title),
395  fNumerator(0x0),
396  fDenominator(0x0),
397  fNBinsToScaleX(fgkDefaultNBinsToScaleX), 
398  fNBinsToScaleY(fgkDefaultNBinsToScaleY)
399 {
400   BuildHistos(nXbins,maxXval,minXval,nYbins,maxYval,minYval);
401 }         
402 /******************************************************************/
403
404 AliHBTFunction2D::~AliHBTFunction2D()
405 {
406   delete fNumerator;
407   delete fDenominator;
408 }
409 /******************************************************************/
410
411 void AliHBTFunction2D::BuildHistos()
412 {
413 //Creates default histograms
414   BuildHistos(fgkDefaultNBinsX,fgkDefaultMaxX,fgkDefaultMinX,
415               fgkDefaultNBinsY,fgkDefaultMaxY,fgkDefaultMinY);
416 }
417 /******************************************************************/
418
419 void AliHBTFunction2D::BuildHistos(Int_t nxbins, Float_t xmax, Float_t xmin,
420                                    Int_t nybins, Float_t ymax, Float_t ymin)
421 {
422   //Builds numerator and denominator histograms (2d-case)
423  TString numstr = fName + " Numerator";  //title and name of the 
424                                            //numerator histogram
425  TString denstr = fName + " Denominator";//title and name of the 
426                                            //denominator histogram
427          
428  fNumerator   = new TH2D(numstr.Data(),numstr.Data(),
429                          nxbins,xmin,xmax,nybins,ymin,ymax);
430                
431  fDenominator = new TH2D(denstr.Data(),denstr.Data(),
432                          nxbins,xmin,xmax,nybins,ymin,ymax);
433  
434  fNumerator->Sumw2();
435  fDenominator->Sumw2();
436 }
437 /******************************************************************/
438
439 void AliHBTFunction2D::SetNumberOfBinsToScale(UInt_t xn, UInt_t yn)
440 {
441 //defines area used for scaling factor calculation
442   fNBinsToScaleX = xn;
443   fNBinsToScaleY = yn;
444 }
445 /******************************************************************/
446
447 Double_t AliHBTFunction2D::Scale()
448 {
449   // Calculates the factor that should be used to scale 
450   // quatience of fNumerator and fDenominator to 1 at 
451   // given region
452   if (gDebug>0) Info("Scale","Enetered Scale()");
453   if(!fNumerator) 
454    {
455      Error("Scale","No numerator");
456      return 0.0;
457    }
458   if(!fDenominator) 
459    {
460      Error("Scale","No denominator");
461      return 0.0;
462    }
463   
464   if( (fNBinsToScaleX < 1) || (fNBinsToScaleY < 1) ) 
465    {
466     return 0.0;
467     Error("Scale","Number of bins for scaling is smaller thnan 1");
468    }
469   UInt_t nbinsX = fNumerator->GetNbinsX();
470   if (fNBinsToScaleX > nbinsX) 
471    {
472     Error("Scale","Number of X bins for scaling is bigger thnan number of bins in histograms");
473     return 0.0;
474    }
475    
476   UInt_t nbinsY = fNumerator->GetNbinsX();
477   if (fNBinsToScaleY > nbinsY) 
478    {
479     Error("Scale","Number of Y bins for scaling is bigger thnan number of bins in histograms");
480     return 0.0;
481    }
482
483   if (gDebug>0) Info("Scale","No errors detected");
484
485   Int_t offsetX = nbinsX - fNBinsToScaleX - 1; //bin that we start loop over bins in axis X
486   Int_t offsetY = nbinsY - fNBinsToScaleY - 1; //bin that we start loop over bins in axis X
487
488   Double_t ratio;
489   Double_t sum = 0;
490   Int_t n = 0;
491   
492   for (UInt_t j = offsetY; j< nbinsY; j++)
493     for (UInt_t i = offsetX; i< nbinsX; i++)
494      {
495       if ( fNumerator->GetBinContent(i,j) > 0.0 )
496        {
497          ratio = fDenominator->GetBinContent(i,j)/fNumerator->GetBinContent(i,j);
498          sum += ratio;
499          n++;
500        }
501      }
502   
503   if(gDebug > 0) Info("Scale","sum=%f fNBinsToScaleX=%d fNBinsToScaleY=%d n=%d",sum,fNBinsToScaleX,fNBinsToScaleY,n);
504   
505   if (n == 0) return 0.0;
506   Double_t ret = sum/((Double_t)n);
507
508   if(gDebug > 0) Info("Scale","returning %f",ret);
509   return ret;
510
511
512 /******************************************************************/
513 /******************************************************************/
514 /******************************************************************/
515
516 //____________________
517 ///////////////////////////////////////////////////////
518 //                                                   //
519 // AliHBTFunction3D                                  //
520 //                                                   //
521 // Base Calss for 3-dimensinal Functions             //
522 //                                                   //
523 // Piotr.Skowronski@cern.ch                          //
524 // http://alisoft.cern.ch/people/skowron/analyzer    //
525 //                                                   //
526 ///////////////////////////////////////////////////////
527
528 ClassImp( AliHBTFunction3D)
529
530 const Int_t AliHBTFunction3D::fgkDefaultNBinsX = 200;//default number of Bins in X axis in histograms
531 const Float_t AliHBTFunction3D::fgkDefaultMinX = 0.0;//Default min value of X axis in histograms
532 const Float_t AliHBTFunction3D::fgkDefaultMaxX = 1.5;//Default max value of X axis inhistograms
533
534 const Int_t AliHBTFunction3D::fgkDefaultNBinsY = 200;//default number of Bins in histograms
535 const Float_t AliHBTFunction3D::fgkDefaultMinY = -0.15;//Default min value of histograms
536 const Float_t AliHBTFunction3D::fgkDefaultMaxY =  0.15;//Default max value of histograms
537
538 const Int_t AliHBTFunction3D::fgkDefaultNBinsZ = 200;//default number of Bins in histograms
539 const Float_t AliHBTFunction3D::fgkDefaultMinZ = -0.15;//Default min value of histograms
540 const Float_t AliHBTFunction3D::fgkDefaultMaxZ =  0.15;//Default max value of histograms
541
542 const UInt_t AliHBTFunction3D::fgkDefaultNBinsToScaleX = 30;//Default number of X bins used for scaling to tale
543 const UInt_t AliHBTFunction3D::fgkDefaultNBinsToScaleY = 30;//Default number of bins used for scaling to tale
544 const UInt_t AliHBTFunction3D::fgkDefaultNBinsToScaleZ = 30;//Default number of bins used for scaling to tale
545
546 AliHBTFunction3D::AliHBTFunction3D():
547  fNumerator(0x0),
548  fDenominator(0x0),
549  fNBinsToScaleX(fgkDefaultNBinsToScaleX), 
550  fNBinsToScaleY(fgkDefaultNBinsToScaleY),
551  fNBinsToScaleZ(fgkDefaultNBinsToScaleZ)
552 {
553  //constructor
554 }
555 /******************************************************************/
556
557 AliHBTFunction3D::AliHBTFunction3D(const Char_t *name, const Char_t *title):
558  AliHBTFunction(name,title),
559  fNumerator(0x0),
560  fDenominator(0x0),
561  fNBinsToScaleX(fgkDefaultNBinsToScaleX), 
562  fNBinsToScaleY(fgkDefaultNBinsToScaleY),
563  fNBinsToScaleZ(fgkDefaultNBinsToScaleZ)
564 {
565  //constructor
566 }  
567 /******************************************************************/
568
569 AliHBTFunction3D::AliHBTFunction3D(Int_t nXbins, Double_t maxXval, Double_t minXval, 
570                   Int_t nYbins, Double_t maxYval, Double_t minYval, 
571                   Int_t nZbins, Double_t maxZval, Double_t minZval):
572  fNumerator(0x0),
573  fDenominator(0x0),
574  fNBinsToScaleX(fgkDefaultNBinsToScaleX), 
575  fNBinsToScaleY(fgkDefaultNBinsToScaleY),
576  fNBinsToScaleZ(fgkDefaultNBinsToScaleZ)
577 {
578  //constructor
579  BuildHistos( nXbins,maxXval,minXval,nYbins,maxYval,minYval,nZbins,maxZval,minZval);
580 }         
581 /******************************************************************/
582
583 AliHBTFunction3D::AliHBTFunction3D(const Char_t *name, const Char_t *title,
584                  Int_t nXbins, Double_t maxXval, Double_t minXval, 
585                  Int_t nYbins, Double_t maxYval, Double_t minYval, 
586                  Int_t nZbins, Double_t maxZval, Double_t minZval):
587  AliHBTFunction(name,title),
588  fNumerator(0x0),
589  fDenominator(0x0),
590  fNBinsToScaleX(fgkDefaultNBinsToScaleX), 
591  fNBinsToScaleY(fgkDefaultNBinsToScaleY),
592  fNBinsToScaleZ(fgkDefaultNBinsToScaleZ)
593 {
594  //constructor
595  BuildHistos( nXbins,maxXval,minXval,nYbins,maxYval,minYval,nZbins,maxZval,minZval);
596 }         
597 /******************************************************************/
598
599
600 AliHBTFunction3D::~AliHBTFunction3D()
601 {
602   delete fNumerator;
603   delete fDenominator;
604 }
605 /******************************************************************/
606
607 void AliHBTFunction3D::BuildHistos()
608 {
609 //Creates default histograms
610   BuildHistos(fgkDefaultNBinsX,fgkDefaultMaxX,fgkDefaultMinX,
611               fgkDefaultNBinsY,fgkDefaultMaxY,fgkDefaultMinY,
612               fgkDefaultNBinsZ,fgkDefaultMaxZ,fgkDefaultMinZ);
613 }
614 /******************************************************************/
615
616 void AliHBTFunction3D::BuildHistos(Int_t nxbins, Float_t xmax, Float_t xmin,
617                                    Int_t nybins, Float_t ymax, Float_t ymin,
618                        Int_t nzbins, Float_t zmax, Float_t zmin)
619 {
620   //Builds numerator and denominator histograms (3d-case)
621    TString numstr = fName + " Numerator";  //title and name of the 
622                                            //numerator histogram
623    TString denstr = fName + " Denominator";//title and name of the 
624                                            //denominator histogram
625          
626    fNumerator   = new TH3D(numstr.Data(),numstr.Data(),
627                            nxbins,xmin,xmin,nybins,ymin,ymin,nzbins,zmin,zmin);
628                
629    fDenominator = new TH3D(denstr.Data(),denstr.Data(),
630                            nxbins,xmin,xmin,nybins,ymin,ymin,nzbins,zmin,zmin);
631    
632    fNumerator->Sumw2();
633    fDenominator->Sumw2();
634 }
635
636 /******************************************************************/
637
638 Double_t AliHBTFunction3D::Scale()
639 {
640   // Calculates the factor that should be used to scale 
641   // quatience of fNumerator and fDenominator to 1 at 
642   // given volume
643   if (gDebug>0) Info("Scale","Enetered Scale()");
644   if(!fNumerator) 
645    {
646      Error("Scale","No numerator");
647      return 0.0;
648    }
649   if(!fDenominator) 
650    {
651      Error("Scale","No denominator");
652      return 0.0;
653    }
654   
655   if( (fNBinsToScaleX < 1) || (fNBinsToScaleY < 1) || (fNBinsToScaleZ < 1)) 
656    {
657     return 0.0;
658     Error("Scale","Number of bins for scaling is smaller thnan 1");
659    }
660   UInt_t nbinsX = fNumerator->GetNbinsX();
661   if (fNBinsToScaleX > nbinsX) 
662    {
663     Error("Scale","Number of X bins for scaling is bigger thnan number of bins in histograms");
664     return 0.0;
665    }
666    
667   UInt_t nbinsY = fNumerator->GetNbinsX();
668   if (fNBinsToScaleY > nbinsY) 
669    {
670     Error("Scale","Number of Y bins for scaling is bigger thnan number of bins in histograms");
671     return 0.0;
672    }
673
674   UInt_t nbinsZ = fNumerator->GetNbinsZ();
675   if (fNBinsToScaleZ > nbinsZ) 
676    {
677     Error("Scale","Number of Z bins for scaling is bigger thnan number of bins in histograms");
678     return 0.0;
679    }
680
681   if (gDebug>0) Info("Scale","No errors detected");
682
683   Int_t offsetX = nbinsX - fNBinsToScaleX - 1; //bin that we start loop over bins in axis X
684   Int_t offsetY = nbinsY - fNBinsToScaleY - 1; //bin that we start loop over bins in axis Y
685   Int_t offsetZ = nbinsZ - fNBinsToScaleZ - 1; //bin that we start loop over bins in axis Z
686
687   Double_t ratio;
688   Double_t sum = 0;
689   Int_t n = 0;
690   
691   for (UInt_t k = offsetZ; k<nbinsZ; k++)
692     for (UInt_t j = offsetY; j<nbinsY; j++)
693       for (UInt_t i = offsetX; i<nbinsX; i++)
694        {
695         if ( fNumerator->GetBinContent(i,j,k) > 0.0 )
696          {
697            ratio = fDenominator->GetBinContent(i,j,k)/fNumerator->GetBinContent(i,j,k);
698            sum += ratio;
699            n++;
700          }
701        }
702   
703   if(gDebug > 0) 
704     Info("Scale","sum=%f fNBinsToScaleX=%d fNBinsToScaleY=%d fNBinsToScaleZ=%d n=%d",
705           sum,fNBinsToScaleX,fNBinsToScaleY,fNBinsToScaleZ,n);
706   
707   if (n == 0) return 0.0;
708   Double_t ret = sum/((Double_t)n);
709
710   if(gDebug > 0) Info("Scale","returning %f",ret);
711   return ret;
712
713 /******************************************************************/
714
715 void AliHBTFunction3D::SetNumberOfBinsToScale(UInt_t xn, UInt_t yn,UInt_t zn)
716 {
717 //sets up the volume to be used for scaling to tail
718  fNBinsToScaleX = xn;
719  fNBinsToScaleY = yn;
720  fNBinsToScaleZ = zn;
721 }
722
723 /******************************************************************/
724 /******************************************************************/
725 /******************************************************************/
726 /******************************************************************/
727 /******************************************************************/
728 /******************************************************************/
729 /******************************************************************/
730 /******************************************************************/
731 /******************************************************************/
732
733 //____________________
734 ///////////////////////////////////////////////////////
735 //                                                   //
736 // AliHBTOnePairFctn1D                               //
737 //                                                   //
738 // Base Calss for 1-dimensinal Functions that need   //
739 // one pair to fill function                         //
740 //                                                   //
741 // Piotr.Skowronski@cern.ch                          //
742 // http://alisoft.cern.ch/people/skowron/analyzer    //
743 //                                                   //
744 ///////////////////////////////////////////////////////
745
746 ClassImp( AliHBTOnePairFctn1D )
747 /******************************************************************/
748
749 AliHBTOnePairFctn1D::AliHBTOnePairFctn1D(Int_t nbins, Float_t maxXval, Float_t minXval):
750  AliHBTFunction1D(nbins,maxXval,minXval)
751 {
752 }
753 /******************************************************************/
754
755 AliHBTOnePairFctn1D::AliHBTOnePairFctn1D(const Char_t *name, const Char_t *title):
756  AliHBTFunction1D(name,title)
757 {
758 }
759 /******************************************************************/
760
761 AliHBTOnePairFctn1D::AliHBTOnePairFctn1D(const Char_t *name, const Char_t *title,
762                                          Int_t nbins, Float_t maxXval, Float_t minXval):
763  AliHBTFunction1D(name,title,nbins,maxXval,minXval)
764 {
765 }
766 /******************************************************************/
767
768 void AliHBTOnePairFctn1D::ProcessSameEventParticles(AliHBTPair* pair)
769 {
770  //Fills the numerator using pair from the same event
771    pair = CheckPair(pair);
772    if(pair) fNumerator->Fill(GetValue(pair));
773 }
774 /******************************************************************/
775 void AliHBTOnePairFctn1D::ProcessDiffEventParticles(AliHBTPair* pair)
776  {
777   //Fills the denumerator using mixed pairs
778    pair = CheckPair(pair);
779    if(pair) fDenominator->Fill(GetValue(pair));
780   }
781
782 /******************************************************************/
783 /******************************************************************/
784 /******************************************************************/
785
786 //____________________
787 ///////////////////////////////////////////////////////
788 //                                                   //
789 // AliHBTOnePairFctn2D                               //
790 //                                                   //
791 // Base Calss for 2-dimensinal Functions that need   //
792 // one pair to fill function                         //
793 //                                                   //
794 // Piotr.Skowronski@cern.ch                          //
795 // http://alisoft.cern.ch/people/skowron/analyzer    //
796 //                                                   //
797 ///////////////////////////////////////////////////////
798
799 ClassImp( AliHBTOnePairFctn2D )
800 /******************************************************************/
801
802 AliHBTOnePairFctn2D::AliHBTOnePairFctn2D(const Char_t *name, const Char_t *title):
803  AliHBTFunction2D(name,title)
804 {
805 }
806 /******************************************************************/
807
808 AliHBTOnePairFctn2D::AliHBTOnePairFctn2D(Int_t nXbins, Double_t maxXval, Double_t minXval, 
809                       Int_t nYbins, Double_t maxYval, Double_t minYval):
810  AliHBTFunction2D(nXbins, maxXval, minXval, nYbins, maxYval, minYval)
811 {
812 }
813 /******************************************************************/
814
815 AliHBTOnePairFctn2D::AliHBTOnePairFctn2D(const Char_t *name, const Char_t *title,
816                       Int_t nXbins, Double_t maxXval, Double_t minXval, 
817                       Int_t nYbins, Double_t maxYval, Double_t minYval):
818  AliHBTFunction2D(name,title, nXbins, maxXval, minXval, nYbins, maxYval, minYval)
819 {
820 }
821 /******************************************************************/
822
823 void AliHBTOnePairFctn2D::ProcessSameEventParticles(AliHBTPair* pair)
824 {
825   // Fills the numerator using pairs from the same event
826   pair = CheckPair(pair);
827   if(pair) 
828    { 
829      Double_t x,y;
830      GetValues(pair,x,y);
831      fNumerator->Fill(x,y);
832    }
833 }
834 /******************************************************************/
835
836 void AliHBTOnePairFctn2D::ProcessDiffEventParticles(AliHBTPair* pair)
837 {
838   // Fills the denumerator using mixed pairs
839   pair = CheckPair(pair);
840   if(pair) 
841    { 
842      Double_t x,y;
843      GetValues(pair,x,y);
844      fDenominator->Fill(x,y);
845    }
846 }
847 /******************************************************************/
848 /******************************************************************/
849 /******************************************************************/
850 /******************************************************************/
851 //____________________
852 ///////////////////////////////////////////////////////
853 //                                                   //
854 // AliHBTOnePairFctn3D                               //
855 //                                                   //
856 // Base Calss for 3-dimensinal Functions that need   //
857 // one pair to fill function                         //
858 //                                                   //
859 // Piotr.Skowronski@cern.ch                          //
860 // http://alisoft.cern.ch/people/skowron/analyzer    //
861 //                                                   //
862 ///////////////////////////////////////////////////////
863 ClassImp( AliHBTOnePairFctn3D)
864
865 /******************************************************************/
866
867 AliHBTOnePairFctn3D::AliHBTOnePairFctn3D(const Char_t *name, const Char_t *title):
868  AliHBTFunction3D(name,title)
869 {
870 }
871 /******************************************************************/
872
873 AliHBTOnePairFctn3D::AliHBTOnePairFctn3D(Int_t nXbins, Double_t maxXval, Double_t minXval, 
874                       Int_t nYbins, Double_t maxYval, Double_t minYval, 
875                       Int_t nZbins, Double_t maxZval, Double_t minZval):
876  AliHBTFunction3D(nXbins, maxXval, minXval, nYbins, maxYval, minYval, nZbins, maxZval, minZval)
877 {
878 }
879 /******************************************************************/
880
881 AliHBTOnePairFctn3D::AliHBTOnePairFctn3D(const Char_t *name, const Char_t *title,
882                       Int_t nXbins, Double_t maxXval, Double_t minXval, 
883                       Int_t nYbins, Double_t maxYval, Double_t minYval, 
884                       Int_t nZbins, Double_t maxZval, Double_t minZval):
885  AliHBTFunction3D(name,title,nXbins, maxXval, minXval, nYbins, maxYval, minYval, nZbins, maxZval, minZval)
886 {
887 }
888 /******************************************************************/
889
890 void AliHBTOnePairFctn3D::ProcessSameEventParticles(AliHBTPair* pair)
891 {
892 //Reacts on pair coming from same event (real pairs)
893 //and fills numerator histogram
894   pair  = CheckPair(pair);
895   if( pair ) 
896    { 
897      Double_t x,y,z;
898      GetValues(pair,x,y,z);
899      fNumerator->Fill(x,y,z);
900    }
901 }
902 /******************************************************************/
903
904 void AliHBTOnePairFctn3D::ProcessDiffEventParticles(AliHBTPair* pair)
905 {
906 //Reacts on pair coming from different events (mixed pairs)
907 //and fills denominator histogram
908   pair  = CheckPair(pair);
909   if( pair ) 
910    { 
911      Double_t x,y,z;
912      GetValues(pair,x,y,z);
913      fDenominator->Fill(x,y,z);
914    }
915 }
916
917 /******************************************************************/
918 /******************************************************************/
919 /******************************************************************/
920
921 //____________________
922 ///////////////////////////////////////////////////////
923 //                                                   //
924 // AliHBTTwoPairFctn1D                               //
925 //                                                   //
926 // Base Calss for 1-dimensinal Functions that need   //
927 // two pair (simulated and reconstructed)            //
928 // to fill function                                  //
929 //                                                   //
930 // Piotr.Skowronski@cern.ch                          //
931 // http://alisoft.cern.ch/people/skowron/analyzer    //
932 //                                                   //
933 ///////////////////////////////////////////////////////
934 ClassImp(AliHBTTwoPairFctn1D)
935 /******************************************************************/
936
937 AliHBTTwoPairFctn1D::AliHBTTwoPairFctn1D(Int_t nbins, Float_t maxXval, Float_t minXval):
938  AliHBTFunction1D(nbins,maxXval,minXval)
939 {
940 }
941 /******************************************************************/
942
943 AliHBTTwoPairFctn1D::AliHBTTwoPairFctn1D(const Char_t *name, const Char_t *title):
944  AliHBTFunction1D(name,title)
945 {
946 }
947 /******************************************************************/
948
949 AliHBTTwoPairFctn1D::AliHBTTwoPairFctn1D(const Char_t *name, const Char_t *title,
950                                          Int_t nbins, Float_t maxXval, Float_t minXval):
951  AliHBTFunction1D(name,title,nbins,maxXval,minXval)
952 {
953 }
954 /******************************************************************/
955
956 void AliHBTTwoPairFctn1D::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
957 {
958   // Fills the numerator using pairs from the same event
959   partpair  = CheckPair(partpair);
960   if( partpair ) 
961    { 
962      Double_t x = GetValue(trackpair,partpair);
963      fNumerator->Fill(x);
964    }
965 }
966 /******************************************************************/
967
968 void AliHBTTwoPairFctn1D::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
969 {
970   // Fills the denumerator usin mixed pairs
971   partpair  = CheckPair(partpair);
972   if( partpair )
973    { 
974      Double_t x = GetValue(trackpair,partpair);
975      fDenominator->Fill(x);
976    }
977 }
978 /******************************************************************/
979 /******************************************************************/
980 /******************************************************************/
981
982 //____________________
983 ///////////////////////////////////////////////////////
984 //                                                   //
985 // AliHBTTwoPairFctn2D                               //
986 //                                                   //
987 // Base Calss for 2-dimensinal Functions that need   //
988 // two pair (simulated and reconstructed)            //
989 // to fill function                                  //
990 //                                                   //
991 // Piotr.Skowronski@cern.ch                          //
992 // http://alisoft.cern.ch/people/skowron/analyzer    //
993 //                                                   //
994 ///////////////////////////////////////////////////////
995
996 ClassImp(AliHBTTwoPairFctn2D)
997
998 /******************************************************************/
999
1000 AliHBTTwoPairFctn2D::AliHBTTwoPairFctn2D(const Char_t *name, const Char_t *title):
1001  AliHBTFunction2D(name,title)
1002 {
1003 }
1004 /******************************************************************/
1005
1006 AliHBTTwoPairFctn2D::AliHBTTwoPairFctn2D(Int_t nXbins, Double_t maxXval, Double_t minXval, 
1007                       Int_t nYbins, Double_t maxYval, Double_t minYval):
1008  AliHBTFunction2D(nXbins, maxXval, minXval, nYbins, maxYval, minYval)
1009 {
1010 }
1011 /******************************************************************/
1012
1013 AliHBTTwoPairFctn2D::AliHBTTwoPairFctn2D(const Char_t *name, const Char_t *title,
1014                       Int_t nXbins, Double_t maxXval, Double_t minXval, 
1015                       Int_t nYbins, Double_t maxYval, Double_t minYval):
1016  AliHBTFunction2D(name,title,nXbins, maxXval, minXval, nYbins, maxYval, minYval)
1017 {
1018 }
1019 /******************************************************************/
1020
1021 void AliHBTTwoPairFctn2D::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
1022 {
1023 //processes pair of particles coming from a same events (real pair)
1024   partpair  = CheckPair(partpair);  //check cuts
1025   if( partpair ) 
1026    { 
1027      Double_t x,y;
1028      GetValues(trackpair,partpair,x,y);
1029      fNumerator->Fill(x,y);
1030    }
1031 }
1032 /******************************************************************/
1033
1034 void AliHBTTwoPairFctn2D::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
1035 {
1036 //processes pair of particles coming from a different events (mixed pair)
1037   partpair  = CheckPair(partpair);
1038   if( partpair ) 
1039    { 
1040      Double_t x,y;
1041      GetValues(trackpair,partpair,x,y);
1042      fDenominator->Fill(x,y);
1043    }
1044 }
1045
1046 /******************************************************************/
1047 /******************************************************************/
1048 /******************************************************************/
1049
1050 //____________________
1051 ///////////////////////////////////////////////////////
1052 //                                                   //
1053 // AliHBTTwoPairFctn3D                               //
1054 //                                                   //
1055 // Base Calss for 3-dimensinal Functions that need   //
1056 // two pair (simulated and reconstructed)            //
1057 // to fill function                                  //
1058 //                                                   //
1059 // Piotr.Skowronski@cern.ch                          //
1060 // http://alisoft.cern.ch/people/skowron/analyzer    //
1061 //                                                   //
1062 ///////////////////////////////////////////////////////
1063
1064 ClassImp(AliHBTTwoPairFctn3D)
1065
1066 /******************************************************************/
1067
1068 AliHBTTwoPairFctn3D::AliHBTTwoPairFctn3D(const Char_t *name, const Char_t *title):
1069  AliHBTFunction3D(name,title)
1070 {
1071 }
1072 /******************************************************************/
1073
1074 AliHBTTwoPairFctn3D::AliHBTTwoPairFctn3D(Int_t nXbins, Double_t maxXval, Double_t minXval, 
1075                       Int_t nYbins, Double_t maxYval, Double_t minYval, 
1076                       Int_t nZbins, Double_t maxZval, Double_t minZval):
1077  AliHBTFunction3D(nXbins, maxXval, minXval, nYbins, maxYval, minYval, nZbins, maxZval, minZval)
1078 {
1079 }
1080 /******************************************************************/
1081
1082 AliHBTTwoPairFctn3D::AliHBTTwoPairFctn3D(const Char_t *name, const Char_t *title,
1083                       Int_t nXbins, Double_t maxXval, Double_t minXval, 
1084                       Int_t nYbins, Double_t maxYval, Double_t minYval, 
1085                       Int_t nZbins, Double_t maxZval, Double_t minZval):
1086  AliHBTFunction3D(name,title,nXbins, maxXval, minXval, nYbins, maxYval, minYval, nZbins, maxZval, minZval)
1087 {
1088 }
1089 /******************************************************************/
1090
1091 void AliHBTTwoPairFctn3D::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
1092 {
1093   // Fills th numerator using pairs from the same event
1094   partpair  = CheckPair(partpair);
1095   if( partpair ) 
1096    { 
1097      Double_t x,y,z;
1098      GetValues(trackpair,partpair,x,y,z);
1099      fNumerator->Fill(x,y,z);
1100    }
1101 }
1102 /******************************************************************/
1103
1104 void AliHBTTwoPairFctn3D::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
1105 {
1106   // Fills the denumerator using mixed pairs
1107   partpair  = CheckPair(partpair);
1108   if( partpair ) 
1109    { 
1110      Double_t x,y,z;
1111      GetValues(trackpair,partpair,x,y,z);
1112      fDenominator->Fill(x,y,z);
1113    }
1114
1115 }
1116
1117
1118 /******************************************************************/
1119 /******************************************************************/
1120 /******************************************************************/
1121