HBTAN/AliHBTFunction.cxx

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