Scaling to tail implemented
[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     two part          four part 
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  
36  fPairCut = new AliHBTEmptyPairCut(); //dummy cut
37 }
38
39 void AliHBTFunction::
40 Write()
41  {
42    if (GetNumerator()) GetNumerator()->Write();
43    if (GetDenominator()) GetDenominator()->Write();
44    TH1* res = GetResult();
45    if (res) res->Write();
46  }
47 /******************************************************************/
48
49 TH1* AliHBTFunction::
50 GetRatio(Double_t normfactor)
51  {
52    //if (gDebug>0) 
53    cout<<"Mormfactor is "<<normfactor<<" for "<<fName<<endl;
54    
55    if (normfactor == 0.0)
56     {
57       Error("GetRatio","Scaling Factor is 0. Null poiner returned");
58       return 0x0;
59     }
60    TString str = fName + " ratio";
61    TH1 *result = (TH1*)GetNumerator()->Clone(str.Data());
62    
63    result->SetTitle(str.Data());
64    //result->Sumw2();
65    
66    result->Divide(GetNumerator(),GetDenominator(),normfactor);
67    
68    return result;
69    
70  }
71 /******************************************************************/
72 void AliHBTFunction::SetPairCut(AliHBTPairCut* cut)
73 {
74 //Sets new Pair Cut. Old one is deleted
75 //Note that it is created new object instead of simple pointer set
76 //I do not want to have pointer 
77 //to object created somewhere else
78 //because in that case I could not believe that 
79 //it would always exist (sb could delete it)
80 //so we have always own copy
81
82  if(!cut) 
83    {
84      Error("AliHBTFunction::SetPairCut","argument is NULL");
85      return;
86    }
87  delete fPairCut;
88  fPairCut = (AliHBTPairCut*)cut->Clone();
89  
90 }
91
92 /******************************************************************/
93
94 void AliHBTFunction::
95 Rename(const Char_t * name)
96  {
97  //renames the function and histograms
98   SetName(name);
99   SetTitle(name);
100   
101   TString numstr = fName + " Numerator";  //title and name of the 
102                                            //numerator histogram
103   TString denstr = fName + " Denominator";//title and name of the 
104                                            //denominator histogram
105  
106   GetNumerator()->SetName(numstr.Data());
107   GetNumerator()->SetTitle(numstr.Data());
108   
109   GetDenominator()->SetName(denstr.Data());
110   GetDenominator()->SetTitle(denstr.Data());
111   
112  }
113
114 void AliHBTFunction::
115 Rename(const Char_t * name, const Char_t * title)
116  {
117  //renames and retitle the function and histograms
118  
119   SetName(name);
120   SetTitle(title);
121   
122   TString numstrn = fName + " Numerator";  //name of the 
123                                            //numerator histogram
124
125   TString numstrt = fTitle + " Numerator";  //title of the 
126                                            //numerator histogram
127                    
128   TString denstrn = fName + " Denominator";//name of the 
129                                            //denominator histogram
130
131   TString denstrt = fTitle + " Denominator";//title of the 
132                                            //denominator histogram
133                    
134  
135   GetNumerator()->SetName(numstrn.Data());
136   GetNumerator()->SetTitle(numstrt.Data());
137   
138   GetDenominator()->SetName(denstrn.Data());
139   GetDenominator()->SetTitle(denstrt.Data());
140
141
142  }
143
144 /******************************************************************/
145 /******************************************************************/
146 /******************************************************************/
147
148 ClassImp( AliHBTTwoPartFctn )
149
150 /******************************************************************/
151 /******************************************************************/
152 /******************************************************************/
153
154 ClassImp( AliHBTFourPartFctn)
155
156 /******************************************************************/
157 /******************************************************************/
158 /******************************************************************/
159
160 ClassImp( AliHBTTwoPartFctn1D )
161
162 AliHBTTwoPartFctn1D::
163 AliHBTTwoPartFctn1D(Int_t nbins, Double_t maxXval, Double_t minXval)
164  {
165  //Constructor of Two Part One Dimentional Function 
166  // nbins: number of bins in histograms - default 100
167  // maxXval and minXval: range of histgram(s) default 0 - 0.15 (GeV)
168  
169  
170    TString numstr = fName + " Numerator";  //title and name of the 
171                                            //numerator histogram
172    TString denstr = fName + " Denominator";//title and name of the 
173                                            //denominator histogram
174          
175    fNumerator   = new TH1D(numstr.Data(),numstr.Data(),nbins,minXval,maxXval);
176    fDenominator = new TH1D(denstr.Data(),denstr.Data(),nbins,minXval,maxXval);
177    
178    fNumerator->Sumw2();
179    fDenominator->Sumw2();
180    
181    fNBinsToScale = 30;
182    
183  }
184 /******************************************************************/
185 AliHBTTwoPartFctn1D::~AliHBTTwoPartFctn1D()
186 {
187   delete fNumerator;
188   delete fDenominator;
189 }
190 /******************************************************************/
191
192 void AliHBTTwoPartFctn1D::ProcessSameEventParticles(AliHBTPair* pair)
193 {
194  //Fills the numerator
195    pair = CheckPair(pair);
196    if(pair) fNumerator->Fill(GetValue(pair));
197 }
198 /******************************************************************/
199 void AliHBTTwoPartFctn1D::ProcessDiffEventParticles(AliHBTPair* pair)
200  {
201   //fills denumerator
202    pair = CheckPair(pair);
203    if(pair) fDenominator->Fill(GetValue(pair));
204
205   }
206 /******************************************************************/
207 Double_t AliHBTTwoPartFctn1D::Scale()
208 {
209   if(!fNumerator) 
210    {
211      Error("Scale","No numerator");
212      return 0.0;
213    }
214   if(!fDenominator) 
215    {
216      Error("Scale","No denominator");
217      return 0.0;
218    }
219   
220   if(fNBinsToScale < 1) 
221    {
222     return 0.0;
223     Error("Scale","Number of bins for scaling is smaller thnan 1");
224    }
225   Int_t nbins = fNumerator->GetNbinsX();
226   if (fNBinsToScale > nbins) 
227    {
228     Error("Scale","Number of bins for scaling is bigger thnan number of bins in histograms");
229     return 0.0;
230    }
231   Double_t ratios[fNBinsToScale];
232
233   Int_t offset = nbins - fNBinsToScale - 1; 
234   Int_t i;
235   for ( i = offset; i< nbins; i++)
236    {
237     if ( fNumerator->GetBinContent(i) == 0.0 )
238      {
239        ratios[i - offset] = -1.0; //since we play with histograms negative is impossible 
240                                   //so it is good flag
241      }
242     else
243      {
244        ratios[i - offset] = fDenominator->GetBinContent(i)/fNumerator->GetBinContent(i);
245      }
246    }
247  
248   Double_t sum = 0;
249   Int_t skipped = 0;
250   for (i = 0; i<fNBinsToScale; i++)
251    {
252     if (ratios[i] == -1.0) skipped++;
253     else sum += ratios[i];
254    }
255   cout<<"sum="<<sum<<" fNBinsToScale="<<fNBinsToScale<<" skipped="<<skipped<<endl;
256
257   return sum/(Double_t)(fNBinsToScale - skipped);
258
259
260 /******************************************************************/
261 /******************************************************************/
262 /******************************************************************/
263
264 ClassImp( AliHBTTwoPartFctn2D )
265
266 AliHBTTwoPartFctn2D::
267 AliHBTTwoPartFctn2D(Int_t nXbins, Double_t maxXval, Double_t minXval , 
268                     Int_t nYbins, Double_t maxYval, Double_t minYval)
269
270 {
271    TString numstr = fName + " Numerator";  //title and name of the 
272                                            //numerator histogram
273    TString denstr = fName + " Denominator";//title and name of the 
274                                            //denominator histogram
275          
276    fNumerator   = new TH2D(numstr.Data(),numstr.Data(),
277                            nXbins,minXval,maxXval,
278                nYbins,minYval,maxYval);
279                
280    fDenominator = new TH2D(denstr.Data(),denstr.Data(),
281                            nXbins,minXval,maxXval,
282                nYbins,minYval,maxYval);
283    
284    fNumerator->Sumw2();
285    fDenominator->Sumw2();
286
287 }         
288 AliHBTTwoPartFctn2D::~AliHBTTwoPartFctn2D()
289 {
290   delete fNumerator;
291   delete fDenominator;
292 }
293 void AliHBTTwoPartFctn2D::ProcessSameEventParticles(AliHBTPair* pair)
294 {
295   pair = CheckPair(pair);
296   if(pair) 
297    { 
298      Double_t x,y;
299      GetValues(pair,x,y);
300      fNumerator->Fill(y,x);
301    }
302 }
303
304 void AliHBTTwoPartFctn2D::ProcessDiffEventParticles(AliHBTPair* pair)
305 {
306   pair = CheckPair(pair);
307   if(pair) 
308    { 
309      Double_t x,y;
310      GetValues(pair,x,y);
311      fDenominator->Fill(y,x);
312    }
313
314 }
315
316
317 /******************************************************************/
318 /******************************************************************/
319 /******************************************************************/
320
321 ClassImp( AliHBTTwoPartFctn3D)
322
323 AliHBTTwoPartFctn3D::
324 AliHBTTwoPartFctn3D(Int_t nXbins, Double_t maxXval, Double_t minXval, 
325                     Int_t nYbins, Double_t maxYval, Double_t minYval, 
326                     Int_t nZbins, Double_t maxZval, Double_t minZval)
327
328 {
329    TString numstr = fName + " Numerator";  //title and name of the 
330                                            //numerator histogram
331    TString denstr = fName + " Denominator";//title and name of the 
332                                            //denominator histogram
333          
334    fNumerator   = new TH3D(numstr.Data(),numstr.Data(),
335                            nXbins,minXval,maxXval,
336                nYbins,minYval,maxYval,
337                nZbins,minZval,maxZval);
338                
339    fDenominator = new TH3D(denstr.Data(),denstr.Data(),
340                            nXbins,minXval,maxXval,
341                nYbins,minYval,maxYval,
342                nZbins,minZval,maxZval);
343    
344    fNumerator->Sumw2();
345    fDenominator->Sumw2();
346
347 }         
348
349
350 AliHBTTwoPartFctn3D::~AliHBTTwoPartFctn3D()
351 {
352   delete fNumerator;
353   delete fDenominator;
354 }
355
356
357 /******************************************************************/
358 /******************************************************************/
359 /******************************************************************/
360 ClassImp( AliHBTFourPartFctn2D)
361
362
363 AliHBTFourPartFctn2D::
364 AliHBTFourPartFctn2D(Int_t nXbins, Double_t maxXval, Double_t minXval , 
365                     Int_t nYbins, Double_t maxYval, Double_t minYval)
366
367 {
368    TString numstr = fName + " Numerator";  //title and name of the 
369                                            //numerator histogram
370    TString denstr = fName + " Denominator";//title and name of the 
371                                            //denominator histogram
372          
373    fNumerator   = new TH2D(numstr.Data(),numstr.Data(),
374                            nXbins,minXval,maxXval,
375                nYbins,minYval,maxYval);
376                
377    fDenominator = new TH2D(denstr.Data(),denstr.Data(),
378                            nXbins,minXval,maxXval,
379                nYbins,minYval,maxYval);
380    
381    fNumerator->Sumw2();
382    fDenominator->Sumw2();
383
384 }         
385 AliHBTFourPartFctn2D::~AliHBTFourPartFctn2D()
386 {
387   delete fNumerator;
388   delete fDenominator;
389 }
390 void AliHBTFourPartFctn2D::
391 ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
392 {
393   partpair  = CheckPair(partpair);
394   trackpair = CheckPair(trackpair);
395   if( partpair && trackpair) 
396    { 
397      Double_t x,y;
398      GetValues(trackpair,partpair,x,y);
399      fNumerator->Fill(y,x);
400    }
401 }
402
403 void AliHBTFourPartFctn2D::
404 ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
405 {
406   partpair  = CheckPair(partpair);
407   trackpair = CheckPair(trackpair);
408   if( partpair && trackpair)
409    { 
410      Double_t x,y;
411      GetValues(trackpair,partpair,x,y);
412      fDenominator->Fill(y,x);
413    }
414
415 }
416