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