HBTAN/AliHBTWeightasCorrFctn.cxx

   1 #include "AliHBTWeightasCorrFctn.h"
   2 #include <TH1.h>
   3 #include <Riostream.h>
   4
   5 ///////////////////////////////////////////////////////
   6 //                                                   //
   7 // AliHBTWeightasCorrFctn.h                             //
   8 //                                                   //
   9 // Class for calculating 3D Weightas correlation        //
  10 // functions                                         //
  11 //                                                   //
  12 ///////////////////////////////////////////////////////
  13
  14 ClassImp(AliHBTWeightasCorrFctn)
  15
  16      AliHBTWeightasCorrFctn::AliHBTWeightasCorrFctn(const char* name, const char* title):
  17  AliHBTTwoPairFctn1D(name,title),
  18
  19      fNum(0x0),
  20      fDen(0x0),
  21      fRat(0x0)
  22
  23 {
  24 //ctor
  25 }
  26
  27 /******************************************************************/
  28 AliHBTWeightasCorrFctn::AliHBTWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
  29  AliHBTTwoPairFctn1D(name,title,nbins,maxXval,minXval),
  30
  31
  32      fNum(new TObjArray()),
  33      fDen(new TObjArray()),
  34      fRat(new TObjArray())
  35 {
  36      SetParams(nbins,maxXval, minXval);
  37 }
  38
  39 /******************************************************************/
  40 AliHBTWeightasCorrFctn::AliHBTWeightasCorrFctn(const AliHBTWeightasCorrFctn& in):
  41  AliHBTTwoPairFctn1D(in),
  42
  43
  44
  45      fNum((in.fNum)?(TObjArray*)in.fNum->Clone():0x0),
  46      fDen((in.fDen)?(TObjArray*)in.fDen->Clone():0x0),
  47      fRat((in.fRat)?(TObjArray*)in.fRat->Clone():0x0)
  48  {
  49 //ctor
  50 }
  51
  52 /******************************************************************/
  53
  54 AliHBTWeightasCorrFctn::~AliHBTWeightasCorrFctn()
  55 {
  56  //dtor
  57
  58      delete fNum;
  59      delete fDen;
  60      delete fRat;
  61
  62 }
  63
  64 /******************************************************************/
  65 void AliHBTWeightasCorrFctn::Write()
  66 {
  67 //out
  68      Int_t i;
  69 //     Int_t n=GetNumberOfIntervals();
  70      Double_t scale;
  71
  72      for(i=0;i<fNumberOfIntervals;i++){
  73           TH1D *num = ((TH1D*)fNum->At(i));
  74           TH1D *den = ((TH1D*)fDen->At(i));
  75           TH1D &rat = *((TH1D*)fRat->At(i));
  76           scale = Scale(num,den);
  77           Info("Write():","Scale in interval %d = %lf",i,scale);
  78           rat.Divide(num,den,scale);
  79
  80
  81           num->Write();
  82           den->Write();
  83           rat.Write();
  84      }
  85
  86 }
  87
  88 //-------------------------------------
  89 void AliHBTWeightasCorrFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
  90 {
  91     //Fills the numerator using pair from the same event
  92      trackpair = CheckPair(trackpair);
  93      if(partpair == 0x0) return;
  94      if(trackpair == 0x0) return;
  95
  96      Double_t weight = 1.0;
  97
  98      int n = fNumberOfIntervals;
  99      Double_t rplane=0.;   //reaction plane angle - 2 B determined
 100      Double_t phi=(trackpair->Particle1()->Phi()+trackpair->Particle2()->Phi())/2.-rplane; //deltaphi bo nie mam nic innego pod reka
 101
 102      phi=phi*360/(2*TMath::Pi());
 103      Double_t q=GetValue(trackpair, partpair);
 104      Int_t ntv;
 105      ntv =  (int)(phi*n/(360.));
 106
 107      TH1D *num = ((TH1D*)fNum->At(ntv));
 108
 109
 110      if ( trackpair && partpair)
 111      {
 112           if ( ( trackpair->Particle1()->GetPdgCode() == partpair->Particle1()->GetPdgCode()) &&
 113                ( trackpair->Particle2()->GetPdgCode() == partpair->Particle2()->GetPdgCode())    )
 114           {
 115                weight=partpair->GetWeight();
 116 //             Info("ProcessSameEvent","weight=%lf",weight);
 117           }
 118           num->Fill(q,weight);
 119      }
 120 }
 121
 122 /****************************************************************/
 123 void AliHBTWeightasCorrFctn::Init()
 124 {
 125      BuildHistos();
 126 }
 127 /****************************************************************/
 128
 129
 130 void AliHBTWeightasCorrFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
 131 {
 132      //process particles from different events
 133      Double_t rplane=0.;   //reaction plane angle - 2 B determined
 134      Double_t phi=(trackpair->Particle1()->Phi()+trackpair->Particle2()->Phi())/2.-rplane; //deltaphi bo nie mam nic innego pod reka
 135      phi=phi*360/(2*TMath::Pi());
 136      trackpair = CheckPair(trackpair);
 137
 138
 139      int n = fNumberOfIntervals;
 140
 141      Int_t ntv;
 142      ntv =  (int)(phi*n/(360.));
 143
 144      TH1D &den = *((TH1D*)fDen->At(ntv));
 145      if ( trackpair && partpair)
 146      {
 147      Double_t qout=GetValue(trackpair, partpair);
 148              den.Fill(qout);
 149      }
 150 }
 151
 152
 153 /******************************************************************/
 154
 155
 156 void AliHBTWeightasCorrFctn::SetParams(Int_t nbins, Float_t maxXval, Float_t minXval)
 157 {
 158 //set histogram parameters
 159      fnbins=nbins;
 160      fmaxXval= maxXval;
 161      fminXval=minXval;
 162 }
 163 TH1* AliHBTWeightasCorrFctn::GetResult()
 164 {
 165 //just for compatibility
 166      TH1D *den = ((TH1D*)fDen->UncheckedAt(1));
 167      return den;
 168  }
 169
 170
 171
 172
 173 ClassImp(AliHBTQOutWeightasCorrFctn)
 174
 175      AliHBTQOutWeightasCorrFctn::AliHBTQOutWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
 176 AliHBTWeightasCorrFctn(name,title,nbins,maxXval,minXval)
 177
 178 {
 179 //ct0r
 180 }
 181
 182 void AliHBTQOutWeightasCorrFctn::BuildHistos()
 183 {
 184 //buils histograms (allocates memory etc)
 185      Int_t i;
 186      int n=GetNumberOfIntervals();
 187
 188      int nbins=Getnbins();
 189
 190      double max = GetmaxXval();
 191      double min = GetminXval();
 192      char buff[10];
 193
 194
 195
 196      TH1D *num;
 197      TH1D *den;
 198      TH1D *rat;
 199
 200      TString nameNum = "WeightNumOut";
 201      TString nameDen = "WeightDenOut";
 202      TString nameRat = "WeightRatOut";
 203
 204      for(i=0;i<n;i++){
 205
 206           sprintf(buff,"%d",i);
 207
 208           nameNum +=TString(buff);
 209
 210           nameDen +=TString(buff);
 211           nameRat +=TString(buff);
 212
 213
 214           num = new TH1D(nameNum.Data(),nameNum.Data(),nbins,min,max);
 215           den = new TH1D(nameDen.Data(),nameDen.Data(),nbins,min,max);
 216           rat = new TH1D(nameRat.Data(),nameRat.Data(),nbins,min,max);
 217
 218           num->Sumw2();
 219           den->Sumw2();
 220           rat->Sumw2();
 221
 222           num->Reset();
 223           den->Reset();
 224           rat->Reset();
 225
 226           fNum->Add(num);
 227           fDen->Add(den);
 228           fRat->Add(rat);
 229
 230           nameNum = TString("WeightNumOut");
 231           nameDen = TString("WeightDenOut");
 232           nameRat = TString("WeightRatOut");
 233
 234      }
 235
 236
 237  }
 238
 239 ClassImp(AliHBTQSideWeightasCorrFctn)
 240
 241      AliHBTQSideWeightasCorrFctn::AliHBTQSideWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
 242 AliHBTWeightasCorrFctn(name,title,nbins,maxXval,minXval)
 243
 244 {
 245 //ct0r
 246 }
 247
 248 void AliHBTQSideWeightasCorrFctn::BuildHistos()
 249 {
 250 //builds histograms
 251      Int_t i;
 252      int n=GetNumberOfIntervals();
 253      int nbins=Getnbins();
 254
 255      double max = GetmaxXval();
 256      double min = GetminXval();
 257      char buff[10];
 258
 259
 260
 261      TH1D *num;
 262      TH1D *den;
 263      TH1D *rat;
 264
 265
 266      TString nameNum = "WeightNumSide";
 267      TString nameDen = "WeightDenSide";
 268      TString nameRat = "WeightRatSide";
 269
 270      for(i=0;i<n;i++){
 271
 272           sprintf(buff,"%d",i);
 273
 274           nameNum +=TString(buff);
 275
 276           nameDen +=TString(buff);
 277           nameRat +=TString(buff);
 278
 279
 280           num = new TH1D(nameNum.Data(),nameNum.Data(),nbins,min,max);
 281           den = new TH1D(nameDen.Data(),nameDen.Data(),nbins,min,max);
 282           rat = new TH1D(nameRat.Data(),nameRat.Data(),nbins,min,max);
 283
 284           num->Sumw2();
 285           den->Sumw2();
 286           rat->Sumw2();
 287
 288           num->Reset();
 289           den->Reset();
 290           rat->Reset();
 291
 292           fNum->Add(num);
 293           fDen->Add(den);
 294           fRat->Add(rat);
 295
 296           nameNum = TString("WeightNumSide");
 297           nameDen = TString("WeightDenSide");
 298           nameRat = TString("WeightRatSide");
 299
 300      }
 301
 302
 303  }
 304
 305 ClassImp(AliHBTQLongWeightasCorrFctn)
 306
 307      AliHBTQLongWeightasCorrFctn::AliHBTQLongWeightasCorrFctn(const char* name, const char* title, Int_t nbins, Float_t maxXval, Float_t minXval):
 308 AliHBTWeightasCorrFctn(name,title,nbins,maxXval,minXval)
 309
 310 {
 311 //ct0r
 312 }
 313
 314 void AliHBTQLongWeightasCorrFctn::BuildHistos()
 315 {
 316 //builds histograms
 317      Int_t i;
 318      int n=GetNumberOfIntervals();
 319      int nbins=Getnbins();
 320
 321      double max = GetmaxXval();
 322      double min = GetminXval();
 323      char buff[10];
 324
 325
 326
 327      TH1D *num;
 328      TH1D *den;
 329      TH1D *rat;
 330
 331      TString nameNum = "WeightNumLong";
 332      TString nameDen = "WeightDenLong";
 333      TString nameRat = "WeightRatLong";
 334
 335      for(i=0;i<n;i++){
 336
 337           sprintf(buff,"%d",i);
 338
 339           nameNum +=TString(buff);
 340
 341           nameDen +=TString(buff);
 342           nameRat +=TString(buff);
 343
 344
 345           num = new TH1D(nameNum.Data(),nameNum.Data(),nbins,min,max);
 346           den = new TH1D(nameDen.Data(),nameDen.Data(),nbins,min,max);
 347           rat = new TH1D(nameRat.Data(),nameRat.Data(),nbins,min,max);
 348
 349           num->Sumw2();
 350           den->Sumw2();
 351           rat->Sumw2();
 352
 353           num->Reset();
 354           den->Reset();
 355           rat->Reset();
 356
 357           fNum->Add(num);
 358           fDen->Add(den);
 359           fRat->Add(rat);
 360
 361           nameNum = TString("WeightNumLong");
 362           nameDen = TString("WeightDenLong");
 363           nameRat = TString("WeightRatLong");
 364
 365      }
 366
 367
 368  }
 369
 370 /********************************************************************/
 371 ClassImp(AliHBTasWeightQOSLCorrFctn)
 372      AliHBTasWeightQOSLCorrFctn::AliHBTasWeightQOSLCorrFctn(const char* name, const char* title):
 373  AliHBTTwoPairFctn3D(name,title),
 374
 375     fNum(0x0),
 376      fDen(0x0),
 377      fRat(0x0)
 378 //      fNumberOfIntervals(1)
 379
 380
 381 {
 382 //ctor
 383 }
 384
 385 /******************************************************************/
 386 AliHBTasWeightQOSLCorrFctn::AliHBTasWeightQOSLCorrFctn(const char* name, const char* title, Int_t nXbins, Double_t maxXval, Double_t minXval, Int_t nYbins, Double_t maxYval, Double_t minYval,Int_t nZbins, Double_t maxZval, Double_t minZval):
 387  AliHBTTwoPairFctn3D(name,title,nXbins,maxXval,minXval,nYbins,maxYval,minYval,nZbins,maxZval,minZval),
 388
 389
 390      fNum(new TObjArray()),
 391      fDen(new TObjArray()),
 392      fRat(new TObjArray())
 393 //      fNumberOfIntervals(1)
 394
 395 {
 396      SetParams(nXbins,maxXval,minXval,nYbins,maxYval,minYval,nZbins,maxZval,minZval);
 397 }
 398
 399 /******************************************************************/
 400 AliHBTasWeightQOSLCorrFctn::AliHBTasWeightQOSLCorrFctn(const AliHBTasWeightQOSLCorrFctn& in):
 401  AliHBTTwoPairFctn3D(in),
 402
 403
 404
 405      fNum((in.fNum)?(TObjArray*)in.fNum->Clone():0x0),
 406      fDen((in.fDen)?(TObjArray*)in.fDen->Clone():0x0),
 407      fRat((in.fRat)?(TObjArray*)in.fRat->Clone():0x0)
 408  {
 409 //ctor
 410 }
 411
 412 /******************************************************************/
 413
 414 AliHBTasWeightQOSLCorrFctn::~AliHBTasWeightQOSLCorrFctn()
 415 {
 416  //dtor
 417
 418      delete fNum;
 419      delete fDen;
 420      delete fRat;
 421
 422 }
 423
 424 /******************************************************************/
 425 void AliHBTasWeightQOSLCorrFctn::Write()
 426 {
 427 //out
 428      Int_t i;
 429 //     Int_t n=GetNumberOfIntervals();
 430     // Double_t scale;
 431
 432      for(i=0;i<fNumberOfIntervals;i++){
 433 Info("Write()","%d",i);
 434              TH3D *num = ((TH3D*)fNum->UncheckedAt(i));
 435           TH3D *den = ((TH3D*)fDen->UncheckedAt(i));
 436           TH3D &rat = *((TH3D*)fRat->At(i));
 437          // scale = Scale(num,den);
 438 //        Info("Write():","Scale in interval %d = %lf",i,scale);
 439           rat.Divide(num,den);//,scale);
 440
 441
 442           num->Write();
 443           den->Write();
 444           rat.Write();
 445      }
 446
 447 }
 448
 449 //-------------------------------------
 450 void AliHBTasWeightQOSLCorrFctn::ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
 451 {
 452     //Fills the numerator using pair from the same event
 453      trackpair = CheckPair(trackpair);
 454      if(partpair == 0x0) return;
 455      if(trackpair == 0x0) return;
 456
 457      Double_t weight = 1.0;
 458
 459      int n = fNumberOfIntervals;
 460      Double_t rplane=0.;   //reaction plane angle - 2 B determined
 461      Double_t phi=(trackpair->Particle1()->Phi()+trackpair->Particle2()->Phi())/2.-rplane; //deltaphi bo nie mam nic innego pod reka
 462
 463      phi=phi*360/(2*TMath::Pi());
 464      Double_t out = TMath::Abs(trackpair->GetQOutLCMS());
 465      Double_t side = TMath::Abs(trackpair->GetQSideLCMS());
 466      Double_t lon = TMath::Abs(trackpair->GetQLongLCMS());
 467
 468
 469      Int_t ntv;
 470      ntv =  (int)(phi*n/(360.));
 471
 472      TH3D *num = ((TH3D*)fNum->At(ntv));
 473
 474
 475      if ( trackpair && partpair)
 476      {
 477           if ( ( trackpair->Particle1()->GetPdgCode() == partpair->Particle1()->GetPdgCode()) &&
 478                ( trackpair->Particle2()->GetPdgCode() == partpair->Particle2()->GetPdgCode())    )
 479           {
 480                weight=partpair->GetWeight();
 481 //             Info("ProcessSameEvent","weight=%lf",weight);
 482           }
 483           num->Fill(out,side,lon,weight);
 484      }
 485 }
 486
 487 /****************************************************************/
 488 void AliHBTasWeightQOSLCorrFctn::Init()
 489 {
 490      BuildHistos();
 491 }
 492 /****************************************************************/
 493
 494
 495 void AliHBTasWeightQOSLCorrFctn::ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
 496 {
 497 //fill denominator
 498      Double_t rplane=0.;   //reaction plane angle - 2 B determined
 499      Double_t phi=(trackpair->Particle1()->Phi()+trackpair->Particle2()->Phi())/2.-rplane; //deltaphi bo nie mam nic innego pod reka
 500      phi=phi*360/(2*TMath::Pi());
 501
 502      Double_t out = TMath::Abs(trackpair->GetQOutLCMS());
 503      Double_t side = TMath::Abs(trackpair->GetQSideLCMS());
 504      Double_t lon = TMath::Abs(trackpair->GetQLongLCMS());
 505
 506      int n = fNumberOfIntervals;
 507
 508      Int_t ntv;
 509      ntv =  (int)(phi*n/(360.));
 510
 511      TH3D &den = *((TH3D*)fDen->At(ntv));
 512      if ( trackpair && partpair)
 513      {
 514           den.Fill(out,side,lon);
 515      }
 516 }
 517
 518
 519 /******************************************************************/
 520
 521
 522 void AliHBTasWeightQOSLCorrFctn::SetParams(Int_t nXbins, Float_t maxXval, Float_t minXval,Int_t nYbins, Float_t maxYval, Float_t minYval,Int_t nZbins, Float_t maxZval, Float_t minZval)
 523 {
 524 //set histogram parametera
 525      fnXbins=nXbins;
 526      fmaxXval= maxXval;
 527      fminXval=minXval;
 528      fnYbins=nYbins;
 529      fmaxYval= maxYval;
 530      fminYval=minYval;
 531
 532      fnZbins=nZbins;
 533      fmaxZval= maxZval;
 534      fminZval=minZval;
 535
 536
 537 }
 538 TH1* AliHBTasWeightQOSLCorrFctn::GetResult()
 539 {
 540  //just for compatibility
 541      TH1D *den = ((TH1D*)fDen->UncheckedAt(1));
 542      return den;
 543  }
 544
 545 void AliHBTasWeightQOSLCorrFctn::BuildHistos()
 546 {
 547     //builds histograms
 548      Int_t i;
 549      int n=GetNumberOfIntervals();
 550      int nXbins=GetnXbins();
 551      double maxX = GetmaxXval();
 552      double minX = GetminXval();
 553
 554      int nYbins=GetnYbins();
 555      double maxY = GetmaxXval();
 556      double minY = GetminXval();
 557      int nZbins=GetnZbins();
 558      double maxZ = GetmaxZval();
 559      double minZ = GetminZval();
 560      char buff[10];
 561
 562      TH3D *num;
 563      TH3D *den;
 564      TH3D *rat;
 565
 566      TString nameNum = "OSLWeightNum";
 567      TString nameDen = "OSLWeightDen";
 568      TString nameRat = "OSLWeightRat";
 569
 570      for(i=0;i<n;i++){
 571
 572           sprintf(buff,"%d",i);
 573
 574           nameNum +=TString(buff);
 575
 576           nameDen +=TString(buff);
 577           nameRat +=TString(buff);
 578
 579
 580           num = new TH3D(nameNum.Data(),nameNum.Data(),nXbins,minX,maxX,nYbins,minY,maxY,nZbins,minZ,maxZ);
 581           den = new TH3D(nameDen.Data(),nameDen.Data(),nXbins,minX,maxX,nYbins,minY,maxY,nZbins,minZ,maxZ);
 582           rat = new TH3D(nameRat.Data(),nameRat.Data(),nXbins,minX,maxX,nYbins,minY,maxY,nZbins,minZ,maxZ);
 583
 584           num->Sumw2();
 585           den->Sumw2();
 586           rat->Sumw2();
 587
 588           num->Reset();
 589           den->Reset();
 590           rat->Reset();
 591
 592           fNum->Add(num);
 593           fDen->Add(den);
 594           fRat->Add(rat);
 595
 596           nameNum = TString("OSLWeightNumLong");
 597           nameDen = TString("OSLWeightDenLong");
 598           nameRat = TString("OSLWeightRatLong");
 599
 600      }
 601 }
 602 Double_t AliHBTasWeightQOSLCorrFctn::Scale(TH3D *num, TH3D *den)
 603 {
 604 //scale histograms method, the default one didn't like me so I've done this one :)
 605 if (AliVAODParticle::GetDebug()>0) Info("Scale","Enetered Scale()");
 606   if(!num)
 607      {
 608          Error("Scale","No numerator");
 609          return 0.0;
 610      }
 611   if(!den)
 612      {
 613          Error("Scale","No denominator");
 614          return 0.0;
 615      }
 616
 617    if( (fNBinsToScaleX < 1) || (fNBinsToScaleY < 1) || (fNBinsToScaleZ < 1))
 618     {
 619          Error("Scale","Number of bins for scaling is smaller thnan 1");
 620          return 0.0;
 621                                             }
 622     UInt_t nbinsX = num->GetNbinsX();
 623   if (fNBinsToScaleX > nbinsX)
 624     {
 625         Error("Scale","Number of X bins for scaling is bigger thnan number of bins in histograms");
 626      return 0.0;
 627     }
 628
 629     UInt_t nbinsY = num->GetNbinsX();
 630    if (fNBinsToScaleY > nbinsY)
 631         {
 632          Error("Scale","Number of Y bins for scaling is bigger thnan number of bins in histograms");
 633         return 0.0;
 634         }
 635
 636    UInt_t nbinsZ = num->GetNbinsZ();
 637   if (fNBinsToScaleZ > nbinsZ)
 638     {
 639       Error("Scale","Number of Z bins for scaling is bigger thnan number of bins in histograms");
 640    return 0.0;
 641      }
 642
 643   if (AliVAODParticle::GetDebug()>0) Info("Scale","No errors detected");
 644     Int_t offsetX = nbinsX - fNBinsToScaleX - 1; //bin that we start loop over bins in axis X
 645     Int_t offsetY = nbinsY - fNBinsToScaleY - 1; //bin that we start loop over bins in axis Y
 646     Int_t offsetZ = nbinsZ - fNBinsToScaleZ - 1; //bin that we start loop over bins in axis Z
 647     Double_t densum = 0.0;
 648     Double_t numsum = 0.0;
 649
 650    for (UInt_t k = offsetZ; k<nbinsZ; k++)
 651      for (UInt_t j = offsetY; j<nbinsY; j++)
 652        for (UInt_t i = offsetX; i<nbinsX; i++)
 653          {
 654             if ( num->GetBinContent(i,j,k) > 0.0 )
 655               {
 656               densum += den->GetBinContent(i,j,k);
 657               numsum += num->GetBinContent(i,j,k);
 658                }
 659          }
 660      if(AliVAODParticle::GetDebug() > 0)
 661      Info("Scale","numsum=%f densum=%f fNBinsToScaleX=%d fNBinsToScaleY=%d fNBinsToScaleZ=%d",
 662               numsum,densum,fNBinsToScaleX,fNBinsToScaleY,fNBinsToScaleZ);
 663    if (numsum == 0) return 0.0;
 664      Double_t ret = densum/numsum;
 665   if(AliVAODParticle::GetDebug() > 0) Info("Scale","returning %f",ret);
 666   return ret;
 667 }
 668