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