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