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