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