HBTAN/AliHBTFunction.cxx

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