PWGCF/FEMTOSCOPY/AliFemtoUser/AliFemtoCorrFctnDEtaDPhiCorrections.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 //                                                                            //
   3 // AliFemtoCorrFctnDEtaDPhiCorrections - A correlation function that analyzes //
   4 // two particle correlations with respect to the azimuthal angle (phi)        //
   5 // and pseudorapidity (eta) difference                                        //
   6 //                                                                            //
   7 // Authors: Adam Kisiel Adam.Kisiel@cern.ch                                   //
   8 //                                                                            //
   9 ////////////////////////////////////////////////////////////////////////////////
  10
  11 #include "AliFemtoCorrFctnDEtaDPhiCorrections.h"
  12 #include "AliFemtoModelHiddenInfo.h"
  13 //#include "AliFemtoHisto.hh"
  14 #include <cstdio>
  15 #include <TMath.h>
  16
  17 #ifdef __ROOT__
  18 ClassImp(AliFemtoCorrFctnDEtaDPhiCorrections)
  19 #endif
  20
  21 #define PIH 1.57079632679489656
  22 #define PIT 6.28318530717958623
  23
  24 //____________________________
  25 AliFemtoCorrFctnDEtaDPhiCorrections::AliFemtoCorrFctnDEtaDPhiCorrections(char* title, const int& aPhiBins=20, const int& aEtaBins=20):
  26   AliFemtoCorrFctn(),
  27   fDPhiDEtaNumerator(0),
  28   fDPhiDEtaDenominator(0),
  29   fDPhiNumerator(0),
  30   fDPhiDenominator(0),
  31   fDCosNumerator(0),
  32   fDCosDenominator(0),
  33   fDoPtAnalysis(0),
  34   fDPhiPtNumerator(0),
  35   fDPhiPtDenominator(0),
  36   fDCosPtNumerator(0),
  37   fDCosPtDenominator(0),
  38   fPhi(0),
  39   fEta(0),
  40   fYtYtNumerator(0),
  41   fYtYtDenominator(0),
  42   fIfCorrectionHist(kNone),
  43   fIfCorrection(0),
  44   fPtCorrectionsNum(0),
  45   fPtCorrectionsDen(0),
  46   fEtaCorrectionsNum(0),
  47   fEtaCorrectionsDen(0),
  48   fCorrFactorTab(0),
  49   fpTab(0),
  50   fPartType(kNoCorrection),
  51   fphiL(0),
  52   fphiT(0)
  53 {
  54
  55   fphiL = (-(int)(aPhiBins/4)+0.5)*2.*TMath::Pi()/aPhiBins;
  56   fphiT = 2*TMath::Pi()+(-(int)(aPhiBins/4)+0.5)*2.*TMath::Pi()/aPhiBins;
  57
  58   // set up numerator
  59   char tTitNumD[101] = "NumDPhiDEta";
  60   strncat(tTitNumD,title, 100);
  61   fDPhiDEtaNumerator = new TH2D(tTitNumD,title,aPhiBins,fphiL,fphiT,aEtaBins,-2.0,2.0);
  62   // set up denominator
  63   char tTitDenD[101] = "DenDPhiDEta";
  64   strncat(tTitDenD,title, 100);
  65   fDPhiDEtaDenominator = new TH2D(tTitDenD,title,aPhiBins,fphiL,fphiT,aEtaBins,-2.0,2.0);
  66
  67   // set up numerator
  68   char tTitNumDPhi[101] = "NumDPhi";
  69   strncat(tTitNumDPhi,title, 100);
  70   fDPhiNumerator = new TH1D(tTitNumDPhi,title,aPhiBins*2,fphiL, fphiT);
  71   // set up denominator
  72   char tTitDenDPhi[101] = "DenDPhi";
  73   strncat(tTitDenDPhi,title, 100);
  74   fDPhiDenominator = new TH1D(tTitDenDPhi,title,aPhiBins*2,fphiL, fphiT);
  75
  76   // set up numerator
  77   char tTitNumDCos[101] = "NumDCos";
  78   strncat(tTitNumDCos,title, 100);
  79   fDCosNumerator = new TH1D(tTitNumDCos,title,aPhiBins*2,-1.0,1.0);
  80   // set up denominator
  81   char tTitDenDCos[101] = "DenDCos";
  82   strncat(tTitDenDCos,title, 100);
  83   fDCosDenominator = new TH1D(tTitDenDCos,title,aPhiBins*2,-1.0,1.0);
  84
  85   char tTitPhi[101] = "Phi";
  86   strncat(tTitPhi,title, 100);
  87   fPhi = new TH1D(tTitPhi,title,90,-TMath::Pi(),TMath::Pi());
  88
  89   char tTitEta[101] = "Eta";
  90   strncat(tTitEta,title, 100);
  91   fEta = new TH1D(tTitEta,title,90,-1.2,1.2);
  92
  93   // set up numerator
  94   char tTitYtNum[101] = "NumYtYt";
  95   strncat(tTitYtNum,title, 100);
  96   fYtYtNumerator = new TH2D(tTitYtNum,title,aPhiBins,1,5,aEtaBins,1,5);
  97   // set up denominator
  98   char tTitYtYtDen[101] = "DenYtYt";
  99   strncat(tTitYtYtDen,title, 100);
 100   fYtYtDenominator = new TH2D(tTitYtYtDen,title,aPhiBins,1,5,aEtaBins,1,5);
 101
 102
 103   char tTitPtCorrectionsNum[101] = "NumpT1pT2EtaPhi";
 104   strncat(tTitPtCorrectionsNum,title, 100);
 105   char tTitPtCorrectionsDen[101] = "DenpT1pT2EtaPhi";
 106   strncat(tTitPtCorrectionsDen,title, 100);
 107
 108   Int_t nbins[4] = {20,20,aPhiBins,aEtaBins};
 109   Double_t xmin[4] = {0,0,-0.5*TMath::Pi(),-2.0};
 110   Double_t xmax[4] = {4,4,1.5*TMath::Pi(),2.0};
 111
 112
 113   fPtCorrectionsNum = new THnSparseF(tTitPtCorrectionsNum,title,4,nbins,xmin,xmax);
 114   fPtCorrectionsDen = new THnSparseF(tTitPtCorrectionsDen,title,4,nbins,xmin,xmax);
 115
 116   char tTitEtaCorrectionsNum[101] = "NumEta1Eta2EtaPhi";
 117   strncat(tTitEtaCorrectionsNum,title, 100);
 118   char tTitEtaCorrectionsDen[101] = "DenEta1Eta2EtaPhi";
 119   strncat(tTitEtaCorrectionsDen,title, 100);
 120
 121   Double_t xmineta[4] = {-1,1,-0.5*TMath::Pi(),-2.0};
 122   Double_t xmaxeta[4] = {-1,1,1.5*TMath::Pi(),2.0};
 123
 124   fEtaCorrectionsNum = new THnSparseF(tTitEtaCorrectionsNum,title,4,nbins,xmineta,xmaxeta);
 125   fEtaCorrectionsDen = new THnSparseF(tTitEtaCorrectionsDen,title,4,nbins,xmineta,xmaxeta);
 126
 127   // THnSparse(const char* name, const char* title, Int_t dim,
 128   //           const Int_t* nbins, const Double_t* xmin, const Double_t* xmax,
 129   //           Int_t chunksize);
 130
 131   // to enable error bar calculation...
 132   fDPhiDEtaNumerator->Sumw2();
 133   fDPhiDEtaDenominator->Sumw2();
 134   fDPhiNumerator->Sumw2();
 135   fDPhiDenominator->Sumw2();
 136   fDCosNumerator->Sumw2();
 137   fDCosDenominator->Sumw2();
 138   fPhi->Sumw2();
 139   fEta->Sumw2();
 140   fYtYtNumerator->Sumw2();
 141   fYtYtDenominator->Sumw2();
 142   fPtCorrectionsNum->Sumw2();
 143   fPtCorrectionsDen->Sumw2();
 144
 145
 146
 147 }
 148
 149 //____________________________
 150 AliFemtoCorrFctnDEtaDPhiCorrections::AliFemtoCorrFctnDEtaDPhiCorrections(const AliFemtoCorrFctnDEtaDPhiCorrections& aCorrFctn) :
 151   AliFemtoCorrFctn(),
 152   fDPhiDEtaNumerator(0),
 153   fDPhiDEtaDenominator(0),
 154   fDPhiNumerator(0),
 155   fDPhiDenominator(0),
 156   fDCosNumerator(0),
 157   fDCosDenominator(0),
 158   fDoPtAnalysis(0),
 159   fDPhiPtNumerator(0),
 160   fDPhiPtDenominator(0),
 161   fDCosPtNumerator(0),
 162   fDCosPtDenominator(0),
 163   fPhi(0),
 164   fEta(0),
 165   fYtYtNumerator(0),
 166   fYtYtDenominator(0),
 167   fIfCorrectionHist(kNone),
 168   fIfCorrection(0),
 169   fPtCorrectionsNum(0),
 170   fPtCorrectionsDen(0),
 171   fEtaCorrectionsNum(0),
 172   fEtaCorrectionsDen(0),
 173   fCorrFactorTab(0),
 174   fpTab(0),
 175   fPartType(kNoCorrection),
 176   fphiL(0),
 177   fphiT(0)
 178 {
 179   // copy constructor
 180   if (aCorrFctn.fDPhiDEtaNumerator)
 181     fDPhiDEtaNumerator = new TH2D(*aCorrFctn.fDPhiDEtaNumerator);
 182   else
 183     fDPhiDEtaNumerator = 0;
 184   if (aCorrFctn.fDPhiDEtaDenominator)
 185     fDPhiDEtaDenominator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
 186   else
 187     fDPhiDEtaDenominator = 0;
 188
 189   if (aCorrFctn.fDPhiNumerator)
 190     fDPhiNumerator = new TH1D(*aCorrFctn.fDPhiNumerator);
 191   else
 192     fDPhiNumerator = 0;
 193   if (aCorrFctn.fDPhiDenominator)
 194     fDPhiDenominator = new TH1D(*aCorrFctn.fDPhiDenominator);
 195   else
 196     fDPhiDenominator = 0;
 197
 198   if (aCorrFctn.fDCosNumerator)
 199     fDCosNumerator = new TH1D(*aCorrFctn.fDCosNumerator);
 200   else
 201     fDCosNumerator = 0;
 202   if (aCorrFctn.fDCosDenominator)
 203     fDCosDenominator = new TH1D(*aCorrFctn.fDCosDenominator);
 204   else
 205     fDCosDenominator = 0;
 206
 207   if (aCorrFctn.fDPhiPtNumerator)
 208     fDPhiPtNumerator = new TH2D(*aCorrFctn.fDPhiPtNumerator);
 209   else
 210     fDPhiPtNumerator = 0;
 211   if (aCorrFctn.fDPhiPtDenominator)
 212     fDPhiPtDenominator = new TH2D(*aCorrFctn.fDPhiPtDenominator);
 213   else
 214     fDPhiPtDenominator = 0;
 215
 216   if (aCorrFctn.fDCosPtNumerator)
 217     fDCosPtNumerator = new TH2D(*aCorrFctn.fDCosPtNumerator);
 218   else
 219     fDCosPtNumerator = 0;
 220   if (aCorrFctn.fDCosPtDenominator)
 221     fDCosPtDenominator = new TH2D(*aCorrFctn.fDCosPtDenominator);
 222   else
 223     fDCosPtDenominator = 0;
 224  if (aCorrFctn.fPhi)
 225     fPhi = new TH1D(*aCorrFctn.fPhi);
 226   else
 227     fPhi = 0;
 228  if (aCorrFctn.fEta)
 229     fEta = new TH1D(*aCorrFctn.fEta);
 230   else
 231     fEta = 0;
 232
 233  if (aCorrFctn.fYtYtNumerator)
 234    fYtYtNumerator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
 235  else
 236    fYtYtNumerator = 0;
 237
 238  if (aCorrFctn.fYtYtDenominator)
 239    fYtYtDenominator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
 240  else
 241     fYtYtDenominator = 0;
 242
 243   fphiL = aCorrFctn.fphiL;
 244   fphiT = aCorrFctn.fphiT;
 245
 246   fPartType = aCorrFctn.fPartType;
 247
 248 //  if (aCorrFctn.fPtCorrectionsNum)
 249 //    fPtCorrectionsNum = new THnSparseF(*aCorrFctn.fPtCorrectionsNum);
 250 //    else
 251 //    fPtCorrectionsNum = 0;
 252
 253 // if (aCorrFctn.fPtCorrectionsDen)
 254 //    fPtCorrectionsDen = new THnSparseF(*aCorrFctn.fPtCorrectionsDen);
 255 //  else
 256 //    fPtCorrectionsDen = 0;
 257
 258
 259
 260 }
 261 //____________________________
 262 AliFemtoCorrFctnDEtaDPhiCorrections::~AliFemtoCorrFctnDEtaDPhiCorrections(){
 263   // destructor
 264   delete fDPhiDEtaNumerator;
 265   delete fDPhiDEtaDenominator;
 266   delete fDPhiNumerator;
 267   delete fDPhiDenominator;
 268   delete fDCosNumerator;
 269   delete fDCosDenominator;
 270   if (fDoPtAnalysis) {
 271     delete fDPhiPtNumerator;
 272     delete fDPhiPtDenominator;
 273     delete fDCosPtNumerator;
 274     delete fDCosPtDenominator;
 275   }
 276   delete fPhi;
 277   delete fEta;
 278
 279   delete fYtYtNumerator;
 280   delete fYtYtDenominator;
 281
 282   delete fPtCorrectionsNum;
 283   delete fPtCorrectionsDen;
 284   delete fEtaCorrectionsNum;
 285   delete fEtaCorrectionsDen;
 286 }
 287 //_________________________
 288 AliFemtoCorrFctnDEtaDPhiCorrections& AliFemtoCorrFctnDEtaDPhiCorrections::operator=(const AliFemtoCorrFctnDEtaDPhiCorrections& aCorrFctn)
 289 {
 290   // assignment operator
 291   if (this == &aCorrFctn)
 292     return *this;
 293
 294   if (aCorrFctn.fDPhiDEtaNumerator)
 295     fDPhiDEtaNumerator = new TH2D(*aCorrFctn.fDPhiDEtaNumerator);
 296   else
 297     fDPhiDEtaNumerator = 0;
 298   if (aCorrFctn.fDPhiDEtaDenominator)
 299     fDPhiDEtaDenominator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
 300   else
 301     fDPhiDEtaDenominator = 0;
 302
 303   if (aCorrFctn.fDPhiNumerator)
 304     fDPhiNumerator = new TH1D(*aCorrFctn.fDPhiNumerator);
 305   else
 306     fDPhiNumerator = 0;
 307   if (aCorrFctn.fDPhiDenominator)
 308     fDPhiDenominator = new TH1D(*aCorrFctn.fDPhiDenominator);
 309   else
 310     fDPhiDenominator = 0;
 311
 312   if (aCorrFctn.fDCosNumerator)
 313     fDCosNumerator = new TH1D(*aCorrFctn.fDCosNumerator);
 314   else
 315     fDCosNumerator = 0;
 316   if (aCorrFctn.fDCosDenominator)
 317     fDCosDenominator = new TH1D(*aCorrFctn.fDCosDenominator);
 318   else
 319     fDCosDenominator = 0;
 320
 321   if (aCorrFctn.fDPhiPtNumerator)
 322     fDPhiPtNumerator = new TH2D(*aCorrFctn.fDPhiPtNumerator);
 323   else
 324     fDPhiPtNumerator = 0;
 325   if (aCorrFctn.fDPhiPtDenominator)
 326     fDPhiPtDenominator = new TH2D(*aCorrFctn.fDPhiPtDenominator);
 327   else
 328     fDPhiPtDenominator = 0;
 329
 330   if (aCorrFctn.fDCosPtNumerator)
 331     fDCosPtNumerator = new TH2D(*aCorrFctn.fDCosPtNumerator);
 332   else
 333     fDCosPtNumerator = 0;
 334   if (aCorrFctn.fDCosPtDenominator)
 335     fDCosPtDenominator = new TH2D(*aCorrFctn.fDCosPtDenominator);
 336   else
 337     fDCosPtDenominator = 0;
 338   if (aCorrFctn.fPhi)
 339     fPhi = new TH1D(*aCorrFctn.fPhi);
 340   else
 341     fPhi = 0;
 342   if (aCorrFctn.fEta)
 343     fEta = new TH1D(*aCorrFctn.fEta);
 344   else
 345     fEta = 0;
 346
 347  if (aCorrFctn.fYtYtNumerator)
 348    fYtYtNumerator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
 349  else
 350    fYtYtNumerator = 0;
 351
 352  if (aCorrFctn.fYtYtDenominator)
 353    fYtYtDenominator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
 354  else
 355    fYtYtDenominator = 0;
 356
 357  fIfCorrectionHist = kNone;
 358  fIfCorrection = 0;
 359
 360   fphiL = aCorrFctn.fphiL;
 361   fphiT = aCorrFctn.fphiT;
 362
 363   fPartType = aCorrFctn.fPartType;
 364
 365 // if (aCorrFctn.fPtCorrectionsNum)
 366 //    fPtCorrectionsNum = new THnSparseF(*aCorrFctn.fPtCorrectionsNum);
 367 //  else
 368 //    fPtCorrectionsNum = 0;
 369
 370 // if (aCorrFctn.fPtCorrectionsDen)
 371 //    fPtCorrectionsDen = new THnSparseF(*aCorrFctn.fPtCorrectionsDen);
 372 //  else
 373 //    fPtCorrectionsDen = 0;
 374
 375
 376
 377   return *this;
 378 }
 379 //_________________________
 380 void AliFemtoCorrFctnDEtaDPhiCorrections::Finish(){
 381   // here is where we should normalize, fit, etc...
 382   // we should NOT Draw() the histos (as I had done it below),
 383   // since we want to insulate ourselves from root at this level
 384   // of the code.  Do it instead at root command line with browser.
 385   //  mShareNumerator->Draw();
 386   //mShareDenominator->Draw();
 387   //mRatio->Draw();
 388
 389 }
 390
 391 //____________________________
 392 AliFemtoString AliFemtoCorrFctnDEtaDPhiCorrections::Report(){
 393   // create report
 394   string stemp = "TPC Ncls Correlation Function Report:\n";
 395   char ctemp[100];
 396   snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fDPhiDEtaNumerator->GetEntries());
 397   stemp += ctemp;
 398   snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDPhiDEtaDenominator->GetEntries());
 399   stemp += ctemp;
 400   //  stemp += mCoulombWeight->Report();
 401   AliFemtoString returnThis = stemp;
 402   return returnThis;
 403 }
 404 //____________________________
 405 void AliFemtoCorrFctnDEtaDPhiCorrections::AddRealPair( AliFemtoPair* pair){
 406   // add real (effect) pair
 407   if (fPairCut)
 408     if (!fPairCut->Pass(pair)) return;
 409
 410   /*double phi1 = pair->Track1()->Track()->P().Phi();
 411   double phi2 = pair->Track2()->Track()->P().Phi();
 412   double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
 413   double eta2 = pair->Track2()->Track()->P().PseudoRapidity();*/
 414
 415   double phi1 = pair->Track1()->FourMomentum().Phi();
 416   double phi2 = pair->Track2()->FourMomentum().Phi();
 417   double eta1 = pair->Track1()->FourMomentum().PseudoRapidity();
 418   double eta2 = pair->Track2()->FourMomentum().PseudoRapidity();
 419
 420   double dphi = phi1 - phi2;
 421   while (dphi<fphiL) dphi+=PIT;
 422   while (dphi>fphiT) dphi-=PIT;
 423
 424   double deta = eta1 - eta2;
 425
 426    double px1 = pair->Track1()->Track()->P().x();
 427    double py1 = pair->Track1()->Track()->P().y();
 428    //double pz1 = pair->Track1()->Track()->P().z();
 429
 430    double px2 = pair->Track2()->Track()->P().x();
 431    double py2 = pair->Track2()->Track()->P().y();
 432    //double pz2 = pair->Track2()->Track()->P().z();
 433
 434    double pt1 = TMath::Hypot(px1, py1);
 435    double pt2 = TMath::Hypot(px2, py2);
 436
 437    double corrweight;
 438    if (fIfCorrection) corrweight = CalculateCorrectionWeight(pt1, pt2);
 439 /*   double ptmin = pt1>pt2 ? pt2 : pt1;
 440
 441    double cosphi = (px1*px2 + py1*py2 + pz1*pz2)/
 442      sqrt((px1*px1 + py1*py1 + pz1*pz1)*(px2*px2 + py2*py2 + pz2*pz2));
 443 */
 444    if (fIfCorrection)
 445       fDPhiDEtaNumerator->Fill(dphi, deta, corrweight);
 446    else
 447       fDPhiDEtaNumerator->Fill(dphi, deta);
 448
 449   fDPhiNumerator->Fill(dphi);
 450 //   fDCosNumerator->Fill(cosphi);
 451
 452   if (fDoPtAnalysis) {
 453 //     fDPhiPtNumerator->Fill(dphi, ptmin);
 454 //     fDCosPtNumerator->Fill(cosphi, ptmin);
 455   }
 456
 457   fPhi->Fill(phi1);
 458   fEta->Fill(eta1);
 459
 460   double PionMass = 0.13956995;
 461   double yt1 = TMath::Log(sqrt(1+(pt1/PionMass)*(pt1/PionMass))+(pt1/PionMass));
 462   double yt2 = TMath::Log(sqrt(1+(pt2/PionMass)*(pt2/PionMass))+(pt2/PionMass));
 463   fYtYtNumerator->Fill(yt1,yt2);
 464
 465   if(fIfCorrectionHist)
 466     {
 467       if(fIfCorrectionHist == kPt){
 468         Double_t val[] = {pt1,pt2,dphi,deta};
 469         fPtCorrectionsNum->Fill(val);
 470       }
 471       if(fIfCorrectionHist == kEta){
 472         Double_t val[] = {eta1,eta2,dphi,deta};
 473         fEtaCorrectionsNum->Fill(val);
 474       }
 475
 476     }
 477
 478 }
 479 //____________________________
 480 void AliFemtoCorrFctnDEtaDPhiCorrections::AddMixedPair( AliFemtoPair* pair){
 481   // add mixed (background) pair
 482   if (fPairCut)
 483     if (!fPairCut->Pass(pair)) return;
 484
 485   /*double phi1 = pair->Track1()->Track()->P().Phi();
 486   double phi2 = pair->Track2()->Track()->P().Phi();
 487   double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
 488   double eta2 = pair->Track2()->Track()->P().PseudoRapidity();*/
 489
 490   double phi1 = pair->Track1()->FourMomentum().Phi();
 491   double phi2 = pair->Track2()->FourMomentum().Phi();
 492   double eta1 = pair->Track1()->FourMomentum().PseudoRapidity();
 493   double eta2 = pair->Track2()->FourMomentum().PseudoRapidity();
 494
 495   double dphi = phi1 - phi2;
 496   while (dphi<fphiL) dphi+=PIT;
 497   while (dphi>fphiT) dphi-=PIT;
 498
 499   double deta = eta1 - eta2;
 500
 501    double px1 = pair->Track1()->Track()->P().x();
 502    double py1 = pair->Track1()->Track()->P().y();
 503    //double pz1 = pair->Track1()->Track()->P().z();
 504
 505    double px2 = pair->Track2()->Track()->P().x();
 506    double py2 = pair->Track2()->Track()->P().y();
 507    //double pz2 = pair->Track2()->Track()->P().z();
 508
 509    double pt1 = TMath::Hypot(px1, py1);
 510    double pt2 = TMath::Hypot(px2, py2);
 511 //   double ptmin = pt1>pt2 ? pt2 : pt1;
 512
 513 //   double cosphi = (px1*px2 + py1*py2 + pz1*pz2)/
 514 //     sqrt((px1*px1 + py1*py1 + pz1*pz1)*(px2*px2 + py2*py2 + pz2*pz2));
 515
 516
 517    double corrweight=-999;
 518    if (fIfCorrection) corrweight = CalculateCorrectionWeight(pt1, pt2);
 519
 520    if(fIfCorrection)
 521       fDPhiDEtaDenominator->Fill(dphi, deta, corrweight);
 522    else
 523       fDPhiDEtaDenominator->Fill(dphi, deta);
 524
 525   fDPhiDenominator->Fill(dphi);
 526 //   fDCosDenominator->Fill(cosphi);
 527
 528   //if (fDoPtAnalysis) {
 529     //   fDPhiPtDenominator->Fill(dphi, ptmin);
 530     //   fDCosPtDenominator->Fill(cosphi, ptmin);
 531   //}
 532
 533   double PionMass = 0.13956995;
 534     double yt1 = TMath::Log(sqrt(1+(pt1/PionMass)*(pt1/PionMass))+(pt1/PionMass));
 535     double yt2 = TMath::Log(sqrt(1+(pt2/PionMass)*(pt2/PionMass))+(pt2/PionMass));
 536     fYtYtDenominator->Fill(yt1,yt2);
 537
 538   if(fIfCorrectionHist)
 539     {
 540       if(fIfCorrectionHist == kPt){
 541         Double_t val[] = {pt1,pt2,dphi,deta};
 542         fPtCorrectionsDen->Fill(val);
 543       }
 544       if(fIfCorrectionHist == kEta){
 545         Double_t val[] = {eta1,eta2,dphi,deta};
 546         fEtaCorrectionsDen->Fill(val);
 547       }
 548     }
 549
 550 }
 551
 552
 553 void AliFemtoCorrFctnDEtaDPhiCorrections::WriteHistos()
 554 {
 555   // Write out result histograms
 556   fDPhiDEtaNumerator->Write();
 557   fDPhiDEtaDenominator->Write();
 558   /*fDPhiNumerator->Write();
 559   fDPhiDenominator->Write();
 560   fDCosNumerator->Write();
 561   fDCosDenominator->Write();
 562   if (fDoPtAnalysis) {
 563     fDPhiPtNumerator->Write();
 564     fDPhiPtDenominator->Write();
 565     fDCosPtNumerator->Write();
 566     fDCosPtDenominator->Write();
 567     }*/
 568   fPhi->Write();
 569   fEta->Write();
 570
 571   if(fIfCorrectionHist){
 572     if(fIfCorrectionHist==kPt){
 573     fPtCorrectionsNum->Write();
 574     fPtCorrectionsDen->Write();}
 575     if(fIfCorrectionHist==kEta){
 576     fEtaCorrectionsNum->Write();
 577     fEtaCorrectionsDen->Write();}
 578   }
 579 }
 580
 581 TList* AliFemtoCorrFctnDEtaDPhiCorrections::GetOutputList()
 582 {
 583   // Prepare the list of objects to be written to the output
 584   TList *tOutputList = new TList();
 585
 586   tOutputList->Add(fDPhiDEtaNumerator);
 587   tOutputList->Add(fDPhiDEtaDenominator);
 588   /*tOutputList->Add(fDPhiNumerator);
 589   tOutputList->Add(fDPhiDenominator);
 590   tOutputList->Add(fDCosNumerator);
 591   tOutputList->Add(fDCosDenominator);
 592   if (fDoPtAnalysis) {
 593     tOutputList->Add(fDPhiPtNumerator);
 594     tOutputList->Add(fDPhiPtDenominator);
 595     tOutputList->Add(fDCosPtNumerator);
 596     tOutputList->Add(fDCosPtDenominator);
 597     }*/
 598   tOutputList->Add(fPhi);
 599   tOutputList->Add(fEta);
 600   tOutputList->Add(fYtYtNumerator);
 601   tOutputList->Add(fYtYtDenominator);
 602
 603   if(fIfCorrectionHist){
 604     if(fIfCorrection==kPt){
 605       tOutputList->Add(fPtCorrectionsNum);
 606       tOutputList->Add(fPtCorrectionsDen);
 607     }
 608     if(fIfCorrectionHist==kEta){
 609       tOutputList->Add(fEtaCorrectionsNum);
 610       tOutputList->Add(fEtaCorrectionsDen);
 611     }
 612   }
 613   return tOutputList;
 614
 615 }
 616
 617 void AliFemtoCorrFctnDEtaDPhiCorrections::SetDoPtAnalysis(int do2d)
 618 {
 619   fDoPtAnalysis = do2d;
 620
 621   int aPhiBins = fDPhiDEtaNumerator->GetNbinsX();
 622   const char *title = fDPhiDEtaNumerator->GetTitle();
 623
 624   // set up numerator
 625   char tTitNumDPhiPt[101] = "NumDPhiPt";
 626   strncat(tTitNumDPhiPt,title, 100);
 627   fDPhiPtNumerator = new TH2D(tTitNumDPhiPt,title,aPhiBins*2,-0.5*TMath::Pi(),3./2.*TMath::Pi(), 30, 0.0, 3.0);
 628   // set up denominator
 629   char tTitDenDPhiPt[101] = "DenDPhiPt";
 630   strncat(tTitDenDPhiPt,title, 100);
 631   fDPhiPtDenominator = new TH2D(tTitDenDPhiPt,title,aPhiBins*2,-0.5*TMath::Pi(),3./2.*TMath::Pi(), 30, 0.0, 3.0);
 632
 633   // set up numerator
 634   char tTitNumDCosPt[101] = "NumDCosPt";
 635   strncat(tTitNumDCosPt,title, 100);
 636   fDCosPtNumerator = new TH2D(tTitNumDCosPt,title,aPhiBins*2,-1.0,1.0, 30, 0.0, 3.0);
 637   // set up denominator
 638   char tTitDenDCosPt[101] = "DenDCosPt";
 639   strncat(tTitDenDCosPt,title, 100);
 640   fDCosPtDenominator = new TH2D(tTitDenDCosPt,title,aPhiBins*2,-1.0,1.0, 30, 0.0, 3.0);
 641
 642   fDPhiPtNumerator->Sumw2();
 643   fDPhiPtDenominator->Sumw2();
 644   fDCosPtNumerator->Sumw2();
 645   fDCosPtDenominator->Sumw2();
 646
 647 }
 648
 649 void AliFemtoCorrFctnDEtaDPhiCorrections::SetDoCorrections(bool doCorr)
 650 {
 651   fIfCorrection = doCorr;
 652 }
 653
 654 void AliFemtoCorrFctnDEtaDPhiCorrections::SetDoCorrectionsHist(CorrectionType doCorr)
 655 {
 656   fIfCorrectionHist = doCorr;
 657 }
 658
 659 void AliFemtoCorrFctnDEtaDPhiCorrections::LoadCorrectionTabFromFile(const char *pTtab, const char *corrTab)
 660 {
 661
 662   double val=-10000;
 663
 664   ifstream ifile1;
 665   ifile1.open(pTtab);
 666   if(ifile1)
 667     {
 668       int nrEntries1;
 669       ifile1>>nrEntries1;
 670       fpTab = new double[nrEntries1];
 671       int i=0;
 672       while(ifile1>>val)
 673         {
 674           fpTab[i] = val;
 675           i++;
 676         }
 677     }
 678   else
 679     {
 680       cout<<"No pT values file open!"<<endl;
 681     }
 682   ifile1.close();
 683
 684
 685   ifstream ifile2;
 686   ifile2.open(corrTab);
 687   if(ifile2)
 688     {
 689       int nrEntries2;
 690       ifile2>>nrEntries2;
 691       fCorrFactorTab = new double[nrEntries2];
 692       int i=0;
 693       while(ifile2>>val)
 694         {
 695           fCorrFactorTab[i] = val;
 696           cout<<"fCorrFactorTab: "<<fCorrFactorTab[i]<<endl;
 697           i++;
 698         }
 699     }
 700   else
 701     {
 702       cout<<"No corrections file open!"<<endl;
 703     }
 704   ifile2.close();
 705 }
 706
 707
 708 double AliFemtoCorrFctnDEtaDPhiCorrections::CalculateCorrectionWeight(double pT1, double pT2)
 709 {
 710
 711     double w1=0., w2=0.;
 712     if(pT1>0 && pT1<5 && pT2>0 && pT2<5)
 713     {
 714       if(pT1<pT2)
 715        {
 716           for (int piter = 0 ; piter<200 ; piter++)
 717           {
 718
 719              if(pT1>= fpTab[piter] && pT1< fpTab[piter+1])
 720                {
 721                  w1=fCorrFactorTab[piter];
 722                }
 723              if(pT2>= fpTab[piter] && pT2< fpTab[piter+1])
 724              {
 725                 w2=fCorrFactorTab[piter];
 726                 break;
 727              }
 728           }
 729        }
 730        else if(pT1>pT2)
 731        {
 732           for (int piter = 0 ; piter<200 ; piter++)
 733           {
 734              if(pT2>= fpTab[piter] && pT2< fpTab[piter+1])
 735                 w2=fCorrFactorTab[piter];
 736              if(pT1>= fpTab[piter] && pT1< fpTab[piter+1])
 737              {
 738                 w1=fCorrFactorTab[piter];
 739                 break;
 740              }
 741           }
 742        }
 743        else //pT1==pT2
 744        {
 745           for (int piter = 0 ; piter<200 ; piter++)
 746           {
 747              if(pT1>= fpTab[piter] && pT1< fpTab[piter+1])
 748              {
 749                 w1=fCorrFactorTab[piter];
 750                 w2=fCorrFactorTab[piter];
 751                 break;
 752              }
 753           }
 754        }
 755        return w1*w2;
 756     }
 757     else
 758        return 0;
 759    return 0;
 760 }