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