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