[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / AliFemtoUser / AliFemtoCorrFctnDEtaDPhi.cxx

////////////////////////////////////////////////////////////////////////////////
//                                                                            //
// AliFemtoCorrFctnDEtaDPhi - A correlation function that analyzes            //
// two particle correlations with respect to the azimuthal angle (phi)        //
// and pseudorapidity (eta) difference                                        //
//                                                                            //
// Authors: Adam Kisiel Adam.Kisiel@cern.ch                                   //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////

#include "AliFemtoCorrFctnDEtaDPhi.h"
#include "AliFemtoModelHiddenInfo.h"
//#include "AliFemtoHisto.hh"
#include <cstdio>
#include <TMath.h>

#ifdef __ROOT__ 
ClassImp(AliFemtoCorrFctnDEtaDPhi)
#endif
  
#define PIH 1.57079632679489656
#define PIT 6.28318530717958623

//____________________________
AliFemtoCorrFctnDEtaDPhi::AliFemtoCorrFctnDEtaDPhi(char* title, const int& aPhiBins=20, const int& aEtaBins=20):
  AliFemtoCorrFctn(),
  fDPhiDEtaNumerator(0),
  fDPhiDEtaDenominator(0),
  fDPhiNumerator(0),
  fDPhiDenominator(0),
  fDCosNumerator(0),
  fDCosDenominator(0),
  fDoPtAnalysis(0),
  fDPhiPtNumerator(0),
  fDPhiPtDenominator(0),
  fDCosPtNumerator(0),
  fDCosPtDenominator(0),
  fPhi(0),
  fEta(0),
  fYtYtNumerator(0),
  fYtYtDenominator(0),
  fIfCorrectionHist(kNone),
  fPtCorrectionsNum(0),
  fPtCorrectionsDen(0),
  fEtaCorrectionsNum(0),
  fEtaCorrectionsDen(0),
  fphiL(0),
  fphiT(0)
{

  fphiL = (-(int)(aPhiBins/4)+0.5)*2.*TMath::Pi()/aPhiBins;
  fphiT = 2*TMath::Pi()+(-(int)(aPhiBins/4)+0.5)*2.*TMath::Pi()/aPhiBins;

  // set up numerator
  char tTitNumD[101] = "NumDPhiDEta";
  strncat(tTitNumD,title, 100);
  fDPhiDEtaNumerator = new TH2D(tTitNumD,title,aPhiBins,fphiL,fphiT,aEtaBins,-2.0,2.0);
  // set up denominator
  char tTitDenD[101] = "DenDPhiDEta";
  strncat(tTitDenD,title, 100);
  fDPhiDEtaDenominator = new TH2D(tTitDenD,title,aPhiBins,fphiL,fphiT,aEtaBins,-2.0,2.0);

  // set up numerator
  char tTitNumDPhi[101] = "NumDPhi";
  strncat(tTitNumDPhi,title, 100);
  fDPhiNumerator = new TH1D(tTitNumDPhi,title,aPhiBins*2,-0.5*TMath::Pi(),1.5*TMath::Pi());
  // set up denominator
  char tTitDenDPhi[101] = "DenDPhi";
  strncat(tTitDenDPhi,title, 100);
  fDPhiDenominator = new TH1D(tTitDenDPhi,title,aPhiBins*2,-0.5*TMath::Pi(),1.5*TMath::Pi());

  // set up numerator
  char tTitNumDCos[101] = "NumDCos";
  strncat(tTitNumDCos,title, 100);
  fDCosNumerator = new TH1D(tTitNumDCos,title,aPhiBins*2,-1.0,1.0);
  // set up denominator
  char tTitDenDCos[101] = "DenDCos";
  strncat(tTitDenDCos,title, 100);
  fDCosDenominator = new TH1D(tTitDenDCos,title,aPhiBins*2,-1.0,1.0);

  char tTitPhi[101] = "Phi";
  strncat(tTitPhi,title, 100);
  fPhi = new TH1D(tTitPhi,title,90,-TMath::Pi(),TMath::Pi());

  char tTitEta[101] = "Eta";
  strncat(tTitEta,title, 100);
  fEta = new TH1D(tTitEta,title,90,-1.2,1.2);
  
  // set up numerator
  char tTitYtNum[101] = "NumYtYt";
  strncat(tTitYtNum,title, 100);
  fYtYtNumerator = new TH2D(tTitYtNum,title,aPhiBins,1,5,aEtaBins,1,5);
  // set up denominator
  char tTitYtYtDen[101] = "DenYtYt";
  strncat(tTitYtYtDen,title, 100);
  fYtYtDenominator = new TH2D(tTitYtYtDen,title,aPhiBins,1,5,aEtaBins,1,5);


  char tTitPtCorrectionsNum[101] = "NumpT1pT2EtaPhi";
  strncat(tTitPtCorrectionsNum,title, 100);
  char tTitPtCorrectionsDen[101] = "DenpT1pT2EtaPhi";
  strncat(tTitPtCorrectionsDen,title, 100);

  Int_t nbins[4] = {20,20,aPhiBins,aEtaBins};
  Double_t xmin[4] = {0,0,-0.5*TMath::Pi(),-2.0};
  Double_t xmax[4] = {4,4,1.5*TMath::Pi(),2.0};


  fPtCorrectionsNum = new THnSparseF(tTitPtCorrectionsNum,title,4,nbins,xmin,xmax);
  fPtCorrectionsDen = new THnSparseF(tTitPtCorrectionsDen,title,4,nbins,xmin,xmax);

  char tTitEtaCorrectionsNum[101] = "NumEta1Eta2EtaPhi";
  strncat(tTitEtaCorrectionsNum,title, 100);
  char tTitEtaCorrectionsDen[101] = "DenEta1Eta2EtaPhi";
  strncat(tTitEtaCorrectionsDen,title, 100);

  Double_t xmineta[4] = {-1,1,-0.5*TMath::Pi(),-2.0};
  Double_t xmaxeta[4] = {-1,1,1.5*TMath::Pi(),2.0};

  fEtaCorrectionsNum = new THnSparseF(tTitEtaCorrectionsNum,title,4,nbins,xmineta,xmaxeta);
  fEtaCorrectionsDen = new THnSparseF(tTitEtaCorrectionsDen,title,4,nbins,xmineta,xmaxeta);

  // THnSparse(const char* name, const char* title, Int_t dim,
  //           const Int_t* nbins, const Double_t* xmin, const Double_t* xmax,
  //           Int_t chunksize);

  // to enable error bar calculation...
  fDPhiDEtaNumerator->Sumw2();
  fDPhiDEtaDenominator->Sumw2();
  fDPhiNumerator->Sumw2();
  fDPhiDenominator->Sumw2();
  fDCosNumerator->Sumw2();
  fDCosDenominator->Sumw2();
  fPhi->Sumw2();
  fEta->Sumw2();
  fYtYtNumerator->Sumw2();
  fYtYtDenominator->Sumw2();
  fPtCorrectionsNum->Sumw2();
  fPtCorrectionsDen->Sumw2();
}

//____________________________
AliFemtoCorrFctnDEtaDPhi::AliFemtoCorrFctnDEtaDPhi(const AliFemtoCorrFctnDEtaDPhi& aCorrFctn) :
  AliFemtoCorrFctn(),
  fDPhiDEtaNumerator(0),
  fDPhiDEtaDenominator(0),
  fDPhiNumerator(0),
  fDPhiDenominator(0),
  fDCosNumerator(0),
  fDCosDenominator(0),
  fDoPtAnalysis(0),
  fDPhiPtNumerator(0),
  fDPhiPtDenominator(0),
  fDCosPtNumerator(0),
  fDCosPtDenominator(0),
  fPhi(0),
  fEta(0),
  fYtYtNumerator(0),
  fYtYtDenominator(0),
  fIfCorrectionHist(kNone),
  fPtCorrectionsNum(0),
  fPtCorrectionsDen(0),
  fEtaCorrectionsNum(0),
  fEtaCorrectionsDen(0),
  fphiL(0),
  fphiT(0)
{
  // copy constructor
  if (aCorrFctn.fDPhiDEtaNumerator)
    fDPhiDEtaNumerator = new TH2D(*aCorrFctn.fDPhiDEtaNumerator);
  else
    fDPhiDEtaNumerator = 0;
  if (aCorrFctn.fDPhiDEtaDenominator)
    fDPhiDEtaDenominator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
  else
    fDPhiDEtaDenominator = 0;

  if (aCorrFctn.fDPhiNumerator)
    fDPhiNumerator = new TH1D(*aCorrFctn.fDPhiNumerator);
  else
    fDPhiNumerator = 0;
  if (aCorrFctn.fDPhiDenominator)
    fDPhiDenominator = new TH1D(*aCorrFctn.fDPhiDenominator);
  else
    fDPhiDenominator = 0;

  if (aCorrFctn.fDCosNumerator)
    fDCosNumerator = new TH1D(*aCorrFctn.fDCosNumerator);
  else
    fDCosNumerator = 0;
  if (aCorrFctn.fDCosDenominator)
    fDCosDenominator = new TH1D(*aCorrFctn.fDCosDenominator);
  else
    fDCosDenominator = 0;

  if (aCorrFctn.fDPhiPtNumerator)
    fDPhiPtNumerator = new TH2D(*aCorrFctn.fDPhiPtNumerator);
  else
    fDPhiPtNumerator = 0;
  if (aCorrFctn.fDPhiPtDenominator)
    fDPhiPtDenominator = new TH2D(*aCorrFctn.fDPhiPtDenominator);
  else
    fDPhiPtDenominator = 0;

  if (aCorrFctn.fDCosPtNumerator)
    fDCosPtNumerator = new TH2D(*aCorrFctn.fDCosPtNumerator);
  else
    fDCosPtNumerator = 0;
  if (aCorrFctn.fDCosPtDenominator)
    fDCosPtDenominator = new TH2D(*aCorrFctn.fDCosPtDenominator);
  else
    fDCosPtDenominator = 0;
 if (aCorrFctn.fPhi)
    fPhi = new TH1D(*aCorrFctn.fPhi);
  else
    fPhi = 0;
 if (aCorrFctn.fEta)
    fEta = new TH1D(*aCorrFctn.fEta);
  else
    fEta = 0;

 if (aCorrFctn.fYtYtNumerator)
   fYtYtNumerator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
 else 
   fYtYtNumerator = 0;

 if (aCorrFctn.fYtYtDenominator)
   fYtYtDenominator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
 else 
   fYtYtDenominator = 0;

  fphiL = aCorrFctn.fphiL;
  fphiT = aCorrFctn.fphiT;

//  if (aCorrFctn.fPtCorrectionsNum)
//    fPtCorrectionsNum = new THnSparseF(*aCorrFctn.fPtCorrectionsNum);
//    else 
//    fPtCorrectionsNum = 0;

// if (aCorrFctn.fPtCorrectionsDen)
//    fPtCorrectionsDen = new THnSparseF(*aCorrFctn.fPtCorrectionsDen);
//  else 
//    fPtCorrectionsDen = 0;


}
//____________________________
AliFemtoCorrFctnDEtaDPhi::~AliFemtoCorrFctnDEtaDPhi(){
  // destructor
  delete fDPhiDEtaNumerator;
  delete fDPhiDEtaDenominator;
  delete fDPhiNumerator;
  delete fDPhiDenominator;
  delete fDCosNumerator;
  delete fDCosDenominator;
  if (fDoPtAnalysis) {
    delete fDPhiPtNumerator;
    delete fDPhiPtDenominator;
    delete fDCosPtNumerator;
    delete fDCosPtDenominator;
  }
  delete fPhi;
  delete fEta;

  delete fYtYtNumerator;
  delete fYtYtDenominator;

  delete fPtCorrectionsNum;
  delete fPtCorrectionsDen;
  delete fEtaCorrectionsNum;
  delete fEtaCorrectionsDen;
}
//_________________________
AliFemtoCorrFctnDEtaDPhi& AliFemtoCorrFctnDEtaDPhi::operator=(const AliFemtoCorrFctnDEtaDPhi& aCorrFctn)
{
  // assignment operator
  if (this == &aCorrFctn)
    return *this;

  if (aCorrFctn.fDPhiDEtaNumerator)
    fDPhiDEtaNumerator = new TH2D(*aCorrFctn.fDPhiDEtaNumerator);
  else
    fDPhiDEtaNumerator = 0;
  if (aCorrFctn.fDPhiDEtaDenominator)
    fDPhiDEtaDenominator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
  else
    fDPhiDEtaDenominator = 0;

  if (aCorrFctn.fDPhiNumerator)
    fDPhiNumerator = new TH1D(*aCorrFctn.fDPhiNumerator);
  else
    fDPhiNumerator = 0;
  if (aCorrFctn.fDPhiDenominator)
    fDPhiDenominator = new TH1D(*aCorrFctn.fDPhiDenominator);
  else
    fDPhiDenominator = 0;

  if (aCorrFctn.fDCosNumerator)
    fDCosNumerator = new TH1D(*aCorrFctn.fDCosNumerator);
  else
    fDCosNumerator = 0;
  if (aCorrFctn.fDCosDenominator)
    fDCosDenominator = new TH1D(*aCorrFctn.fDCosDenominator);
  else
    fDCosDenominator = 0;

  if (aCorrFctn.fDPhiPtNumerator)
    fDPhiPtNumerator = new TH2D(*aCorrFctn.fDPhiPtNumerator);
  else
    fDPhiPtNumerator = 0;
  if (aCorrFctn.fDPhiPtDenominator)
    fDPhiPtDenominator = new TH2D(*aCorrFctn.fDPhiPtDenominator);
  else
    fDPhiPtDenominator = 0;

  if (aCorrFctn.fDCosPtNumerator)
    fDCosPtNumerator = new TH2D(*aCorrFctn.fDCosPtNumerator);
  else
    fDCosPtNumerator = 0;
  if (aCorrFctn.fDCosPtDenominator)
    fDCosPtDenominator = new TH2D(*aCorrFctn.fDCosPtDenominator);
  else
    fDCosPtDenominator = 0;
  if (aCorrFctn.fPhi)
    fPhi = new TH1D(*aCorrFctn.fPhi);
  else
    fPhi = 0;
  if (aCorrFctn.fEta)
    fEta = new TH1D(*aCorrFctn.fEta);
  else
    fEta = 0;

 if (aCorrFctn.fYtYtNumerator)
   fYtYtNumerator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
 else 
   fYtYtNumerator = 0;

 if (aCorrFctn.fYtYtDenominator)
   fYtYtDenominator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
 else 
   fYtYtDenominator = 0;

 fIfCorrectionHist = kNone;

  fphiL = aCorrFctn.fphiL;
  fphiT = aCorrFctn.fphiT;

// if (aCorrFctn.fPtCorrectionsNum)
//    fPtCorrectionsNum = new THnSparseF(*aCorrFctn.fPtCorrectionsNum);
//  else 
//    fPtCorrectionsNum = 0;

// if (aCorrFctn.fPtCorrectionsDen)
//    fPtCorrectionsDen = new THnSparseF(*aCorrFctn.fPtCorrectionsDen);
//  else 
//    fPtCorrectionsDen = 0;


  return *this;
}
//_________________________
void AliFemtoCorrFctnDEtaDPhi::Finish(){
  // here is where we should normalize, fit, etc...
  // we should NOT Draw() the histos (as I had done it below),
  // since we want to insulate ourselves from root at this level
  // of the code.  Do it instead at root command line with browser.
  //  mShareNumerator->Draw();
  //mShareDenominator->Draw();
  //mRatio->Draw();

}

//____________________________
AliFemtoString AliFemtoCorrFctnDEtaDPhi::Report(){
  // create report
  string stemp = "TPC Ncls Correlation Function Report:\n";
  char ctemp[100];
  snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fDPhiDEtaNumerator->GetEntries());
  stemp += ctemp;
  snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDPhiDEtaDenominator->GetEntries());
  stemp += ctemp;
  //  stemp += mCoulombWeight->Report();
  AliFemtoString returnThis = stemp;
  return returnThis;
}
//____________________________
void AliFemtoCorrFctnDEtaDPhi::AddRealPair( AliFemtoPair* pair){
  // add real (effect) pair
  if (fPairCut)
    if (!fPairCut->Pass(pair)) return;

  /*double phi1 = pair->Track1()->Track()->P().Phi();
  double phi2 = pair->Track2()->Track()->P().Phi();
  double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
  double eta2 = pair->Track2()->Track()->P().PseudoRapidity();*/

  double phi1 = pair->Track1()->FourMomentum().Phi();
  double phi2 = pair->Track2()->FourMomentum().Phi();
  double eta1 = pair->Track1()->FourMomentum().PseudoRapidity();
  double eta2 = pair->Track2()->FourMomentum().PseudoRapidity();

  double dphi = phi1 - phi2;
  while (dphi<fphiL) dphi+=PIT;
  while (dphi>fphiT) dphi-=PIT;

  double deta = eta1 - eta2;

   double px1 = pair->Track1()->Track()->P().x();
   double py1 = pair->Track1()->Track()->P().y();
   //double pz1 = pair->Track1()->Track()->P().z();

   double px2 = pair->Track2()->Track()->P().x();
   double py2 = pair->Track2()->Track()->P().y();
   //double pz2 = pair->Track2()->Track()->P().z();

   double pt1 = TMath::Hypot(px1, py1);
   double pt2 = TMath::Hypot(px2, py2);
//   double ptmin = pt1>pt2 ? pt2 : pt1;

//   double cosphi = (px1*px2 + py1*py2 + pz1*pz2)/
//     sqrt((px1*px1 + py1*py1 + pz1*pz1)*(px2*px2 + py2*py2 + pz2*pz2));

  fDPhiDEtaNumerator->Fill(dphi, deta);

  fDPhiNumerator->Fill(dphi);
//   fDCosNumerator->Fill(cosphi);

  if (fDoPtAnalysis) {
//     fDPhiPtNumerator->Fill(dphi, ptmin);
//     fDCosPtNumerator->Fill(cosphi, ptmin);
  }

  fPhi->Fill(phi1);
  fEta->Fill(eta1);

  double PionMass = 0.13956995;
  double yt1 = TMath::Log(sqrt(1+(pt1/PionMass)*(pt1/PionMass))+(pt1/PionMass));
  double yt2 = TMath::Log(sqrt(1+(pt2/PionMass)*(pt2/PionMass))+(pt2/PionMass));
  fYtYtNumerator->Fill(yt1,yt2);

  if(fIfCorrectionHist)
    {
      if(fIfCorrectionHist == kPt){
	Double_t val[] = {pt1,pt2,dphi,deta};
	fPtCorrectionsNum->Fill(val);
      }
      if(fIfCorrectionHist == kEta){
	Double_t val[] = {eta1,eta2,dphi,deta};
	fEtaCorrectionsNum->Fill(val);
      }

    }

}
//____________________________
void AliFemtoCorrFctnDEtaDPhi::AddMixedPair( AliFemtoPair* pair){
  // add mixed (background) pair
  if (fPairCut)
    if (!fPairCut->Pass(pair)) return;

  /*double phi1 = pair->Track1()->Track()->P().Phi();
  double phi2 = pair->Track2()->Track()->P().Phi();
  double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
  double eta2 = pair->Track2()->Track()->P().PseudoRapidity();*/

  double phi1 = pair->Track1()->FourMomentum().Phi();
  double phi2 = pair->Track2()->FourMomentum().Phi();
  double eta1 = pair->Track1()->FourMomentum().PseudoRapidity();
  double eta2 = pair->Track2()->FourMomentum().PseudoRapidity();

  double dphi = phi1 - phi2;
  while (dphi<fphiL) dphi+=PIT;
  while (dphi>fphiT) dphi-=PIT;

  double deta = eta1 - eta2;

   double px1 = pair->Track1()->Track()->P().x();
   double py1 = pair->Track1()->Track()->P().y();
   //double pz1 = pair->Track1()->Track()->P().z();

   double px2 = pair->Track2()->Track()->P().x();
   double py2 = pair->Track2()->Track()->P().y();
   //double pz2 = pair->Track2()->Track()->P().z();

   double pt1 = TMath::Hypot(px1, py1);
   double pt2 = TMath::Hypot(px2, py2);
//   double ptmin = pt1>pt2 ? pt2 : pt1;

//   double cosphi = (px1*px2 + py1*py2 + pz1*pz2)/
//     sqrt((px1*px1 + py1*py1 + pz1*pz1)*(px2*px2 + py2*py2 + pz2*pz2));

  fDPhiDEtaDenominator->Fill(dphi, deta);

  fDPhiDenominator->Fill(dphi);
//   fDCosDenominator->Fill(cosphi);

  //if (fDoPtAnalysis) {
    //   fDPhiPtDenominator->Fill(dphi, ptmin);
    //   fDCosPtDenominator->Fill(cosphi, ptmin);
  //}

  double PionMass = 0.13956995;
    double yt1 = TMath::Log(sqrt(1+(pt1/PionMass)*(pt1/PionMass))+(pt1/PionMass));
    double yt2 = TMath::Log(sqrt(1+(pt2/PionMass)*(pt2/PionMass))+(pt2/PionMass));
    fYtYtDenominator->Fill(yt1,yt2);

  if(fIfCorrectionHist)
    {
      if(fIfCorrectionHist == kPt){
	Double_t val[] = {pt1,pt2,dphi,deta};
	fPtCorrectionsDen->Fill(val);
      }
      if(fIfCorrectionHist == kEta){
	Double_t val[] = {eta1,eta2,dphi,deta};
	fEtaCorrectionsDen->Fill(val);
      }
    }

}


void AliFemtoCorrFctnDEtaDPhi::WriteHistos()
{
  // Write out result histograms
  fDPhiDEtaNumerator->Write();
  fDPhiDEtaDenominator->Write();
  /*fDPhiNumerator->Write();
  fDPhiDenominator->Write();
  fDCosNumerator->Write();
  fDCosDenominator->Write();*/
  if (fDoPtAnalysis) {
    fDPhiPtNumerator->Write();
    fDPhiPtDenominator->Write();
    fDCosPtNumerator->Write();
    fDCosPtDenominator->Write();
  }
  fPhi->Write();
  fEta->Write();
  
  if(fIfCorrectionHist){
    if(fIfCorrectionHist==kPt){
    fPtCorrectionsNum->Write();
    fPtCorrectionsDen->Write();}
    if(fIfCorrectionHist==kEta){
    fEtaCorrectionsNum->Write();
    fEtaCorrectionsDen->Write();}
  }
}

TList* AliFemtoCorrFctnDEtaDPhi::GetOutputList()
{
  // Prepare the list of objects to be written to the output
  TList *tOutputList = new TList();

  tOutputList->Add(fDPhiDEtaNumerator);
  tOutputList->Add(fDPhiDEtaDenominator);
  /*tOutputList->Add(fDPhiNumerator);
  tOutputList->Add(fDPhiDenominator);
  tOutputList->Add(fDCosNumerator);
  tOutputList->Add(fDCosDenominator);
  if (fDoPtAnalysis) {
    tOutputList->Add(fDPhiPtNumerator);
    tOutputList->Add(fDPhiPtDenominator);
    tOutputList->Add(fDCosPtNumerator);
    tOutputList->Add(fDCosPtDenominator);
    }*/
  tOutputList->Add(fPhi);
  tOutputList->Add(fEta);
  tOutputList->Add(fYtYtNumerator);
  tOutputList->Add(fYtYtDenominator);

  if(fIfCorrectionHist){
    if(fIfCorrectionHist==kPt){
      tOutputList->Add(fPtCorrectionsNum);
      tOutputList->Add(fPtCorrectionsDen);
    }
    if(fIfCorrectionHist==kEta){
      tOutputList->Add(fEtaCorrectionsNum);
      tOutputList->Add(fEtaCorrectionsDen);
    }
  }
  return tOutputList;

}

void AliFemtoCorrFctnDEtaDPhi::SetDoPtAnalysis(int do2d)
{
  fDoPtAnalysis = do2d;
  
  int aPhiBins = fDPhiDEtaNumerator->GetNbinsX();
  const char *title = fDPhiDEtaNumerator->GetTitle();

  // set up numerator
  char tTitNumDPhiPt[101] = "NumDPhiPt";
  strncat(tTitNumDPhiPt,title, 100);
  fDPhiPtNumerator = new TH2D(tTitNumDPhiPt,title,aPhiBins*2,-0.5*TMath::Pi(),3./2.*TMath::Pi(), 30, 0.0, 3.0);
  // set up denominator
  char tTitDenDPhiPt[101] = "DenDPhiPt";
  strncat(tTitDenDPhiPt,title, 100);
  fDPhiPtDenominator = new TH2D(tTitDenDPhiPt,title,aPhiBins*2,-0.5*TMath::Pi(),3./2.*TMath::Pi(), 30, 0.0, 3.0);

  // set up numerator
  char tTitNumDCosPt[101] = "NumDCosPt";
  strncat(tTitNumDCosPt,title, 100);
  fDCosPtNumerator = new TH2D(tTitNumDCosPt,title,aPhiBins*2,-1.0,1.0, 30, 0.0, 3.0);
  // set up denominator
  char tTitDenDCosPt[101] = "DenDCosPt";
  strncat(tTitDenDCosPt,title, 100);
  fDCosPtDenominator = new TH2D(tTitDenDCosPt,title,aPhiBins*2,-1.0,1.0, 30, 0.0, 3.0);

  fDPhiPtNumerator->Sumw2();
  fDPhiPtDenominator->Sumw2();
  fDCosPtNumerator->Sumw2();
  fDCosPtDenominator->Sumw2();
  
}

void AliFemtoCorrFctnDEtaDPhi::SetDo4DCorrectionHist(CorrectionType doCorr)
{
  fIfCorrectionHist = doCorr;
}
Commit	Line	Data
76ce4b5b	1	////////////////////////////////////////////////////////////////////////////////
	2	// //
	3	// AliFemtoCorrFctnDEtaDPhi - 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 "AliFemtoCorrFctnDEtaDPhi.h"
	12	#include "AliFemtoModelHiddenInfo.h"
	13	//#include "AliFemtoHisto.hh"
	14	#include <cstdio>
	15	#include <TMath.h>
	16
	17	#ifdef __ROOT__
	18	ClassImp(AliFemtoCorrFctnDEtaDPhi)
	19	#endif
	20
	21	#define PIH 1.57079632679489656
	22	#define PIT 6.28318530717958623
76ce4b5b	23
	24	//____________________________
	25	AliFemtoCorrFctnDEtaDPhi::AliFemtoCorrFctnDEtaDPhi(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),
973a91f8	37	fDCosPtDenominator(0),
973a91f8	38	fPhi(0),
3b97e69c	39	fEta(0),
3b97e69c	40	fYtYtNumerator(0),
d381642d	41	fYtYtDenominator(0),
41fa2a11	42	fIfCorrectionHist(kNone),
d381642d	43	fPtCorrectionsNum(0),
	44	fPtCorrectionsDen(0),
	45	fEtaCorrectionsNum(0),
0cf72019	46	fEtaCorrectionsDen(0),
	47	fphiL(0),
	48	fphiT(0)
76ce4b5b	49	{
4d019bab	50
	51	fphiL = (-(int)(aPhiBins/4)+0.5)2.TMath::Pi()/aPhiBins;
	52	fphiT = 2TMath::Pi()+(-(int)(aPhiBins/4)+0.5)2.*TMath::Pi()/aPhiBins;
	53
76ce4b5b	54	// set up numerator
	55	char tTitNumD[101] = "NumDPhiDEta";
	56	strncat(tTitNumD,title, 100);
4d019bab	57	fDPhiDEtaNumerator = new TH2D(tTitNumD,title,aPhiBins,fphiL,fphiT,aEtaBins,-2.0,2.0);
76ce4b5b	58	// set up denominator
	59	char tTitDenD[101] = "DenDPhiDEta";
	60	strncat(tTitDenD,title, 100);
4d019bab	61	fDPhiDEtaDenominator = new TH2D(tTitDenD,title,aPhiBins,fphiL,fphiT,aEtaBins,-2.0,2.0);
76ce4b5b	62
	63	// set up numerator
	64	char tTitNumDPhi[101] = "NumDPhi";
	65	strncat(tTitNumDPhi,title, 100);
	66	fDPhiNumerator = new TH1D(tTitNumDPhi,title,aPhiBins2,-0.5TMath::Pi(),1.5*TMath::Pi());
	67	// set up denominator
	68	char tTitDenDPhi[101] = "DenDPhi";
	69	strncat(tTitDenDPhi,title, 100);
	70	fDPhiDenominator = new TH1D(tTitDenDPhi,title,aPhiBins2,-0.5TMath::Pi(),1.5*TMath::Pi());
	71
	72	// set up numerator
	73	char tTitNumDCos[101] = "NumDCos";
	74	strncat(tTitNumDCos,title, 100);
	75	fDCosNumerator = new TH1D(tTitNumDCos,title,aPhiBins*2,-1.0,1.0);
	76	// set up denominator
	77	char tTitDenDCos[101] = "DenDCos";
	78	strncat(tTitDenDCos,title, 100);
	79	fDCosDenominator = new TH1D(tTitDenDCos,title,aPhiBins*2,-1.0,1.0);
	80
973a91f8	81	char tTitPhi[101] = "Phi";
	82	strncat(tTitPhi,title, 100);
	83	fPhi = new TH1D(tTitPhi,title,90,-TMath::Pi(),TMath::Pi());
	84
	85	char tTitEta[101] = "Eta";
	86	strncat(tTitEta,title, 100);
	87	fEta = new TH1D(tTitEta,title,90,-1.2,1.2);
	88
3b97e69c	89	// set up numerator
	90	char tTitYtNum[101] = "NumYtYt";
	91	strncat(tTitYtNum,title, 100);
	92	fYtYtNumerator = new TH2D(tTitYtNum,title,aPhiBins,1,5,aEtaBins,1,5);
	93	// set up denominator
	94	char tTitYtYtDen[101] = "DenYtYt";
	95	strncat(tTitYtYtDen,title, 100);
	96	fYtYtDenominator = new TH2D(tTitYtYtDen,title,aPhiBins,1,5,aEtaBins,1,5);
	97
	98
d381642d	99	char tTitPtCorrectionsNum[101] = "NumpT1pT2EtaPhi";
	100	strncat(tTitPtCorrectionsNum,title, 100);
	101	char tTitPtCorrectionsDen[101] = "DenpT1pT2EtaPhi";
	102	strncat(tTitPtCorrectionsDen,title, 100);
	103
	104	Int_t nbins[4] = {20,20,aPhiBins,aEtaBins};
	105	Double_t xmin[4] = {0,0,-0.5*TMath::Pi(),-2.0};
	106	Double_t xmax[4] = {4,4,1.5*TMath::Pi(),2.0};
	107
	108
	109	fPtCorrectionsNum = new THnSparseF(tTitPtCorrectionsNum,title,4,nbins,xmin,xmax);
	110	fPtCorrectionsDen = new THnSparseF(tTitPtCorrectionsDen,title,4,nbins,xmin,xmax);
	111
	112	char tTitEtaCorrectionsNum[101] = "NumEta1Eta2EtaPhi";
	113	strncat(tTitEtaCorrectionsNum,title, 100);
	114	char tTitEtaCorrectionsDen[101] = "DenEta1Eta2EtaPhi";
	115	strncat(tTitEtaCorrectionsDen,title, 100);
	116
	117	Double_t xmineta[4] = {-1,1,-0.5*TMath::Pi(),-2.0};
	118	Double_t xmaxeta[4] = {-1,1,1.5*TMath::Pi(),2.0};
	119
	120	fEtaCorrectionsNum = new THnSparseF(tTitEtaCorrectionsNum,title,4,nbins,xmineta,xmaxeta);
	121	fEtaCorrectionsDen = new THnSparseF(tTitEtaCorrectionsDen,title,4,nbins,xmineta,xmaxeta);
	122
	123	// THnSparse(const char* name, const char* title, Int_t dim,
	124	// const Int_t* nbins, const Double_t* xmin, const Double_t* xmax,
	125	// Int_t chunksize);
973a91f8	126
76ce4b5b	127	// to enable error bar calculation...
	128	fDPhiDEtaNumerator->Sumw2();
	129	fDPhiDEtaDenominator->Sumw2();
	130	fDPhiNumerator->Sumw2();
	131	fDPhiDenominator->Sumw2();
	132	fDCosNumerator->Sumw2();
	133	fDCosDenominator->Sumw2();
973a91f8	134	fPhi->Sumw2();
973a91f8	135	fEta->Sumw2();
3b97e69c	136	fYtYtNumerator->Sumw2();
3b97e69c	137	fYtYtDenominator->Sumw2();
d381642d	138	fPtCorrectionsNum->Sumw2();
d381642d	139	fPtCorrectionsDen->Sumw2();
76ce4b5b	140	}
	141
	142	//____________________________
	143	AliFemtoCorrFctnDEtaDPhi::AliFemtoCorrFctnDEtaDPhi(const AliFemtoCorrFctnDEtaDPhi& aCorrFctn) :
	144	AliFemtoCorrFctn(),
	145	fDPhiDEtaNumerator(0),
	146	fDPhiDEtaDenominator(0),
	147	fDPhiNumerator(0),
	148	fDPhiDenominator(0),
	149	fDCosNumerator(0),
	150	fDCosDenominator(0),
	151	fDoPtAnalysis(0),
	152	fDPhiPtNumerator(0),
	153	fDPhiPtDenominator(0),
	154	fDCosPtNumerator(0),
973a91f8	155	fDCosPtDenominator(0),
973a91f8	156	fPhi(0),
3b97e69c	157	fEta(0),
3b97e69c	158	fYtYtNumerator(0),
d381642d	159	fYtYtDenominator(0),
41fa2a11	160	fIfCorrectionHist(kNone),
d381642d	161	fPtCorrectionsNum(0),
	162	fPtCorrectionsDen(0),
	163	fEtaCorrectionsNum(0),
0cf72019	164	fEtaCorrectionsDen(0),
	165	fphiL(0),
	166	fphiT(0)
76ce4b5b	167	{
	168	// copy constructor
	169	if (aCorrFctn.fDPhiDEtaNumerator)
	170	fDPhiDEtaNumerator = new TH2D(*aCorrFctn.fDPhiDEtaNumerator);
	171	else
	172	fDPhiDEtaNumerator = 0;
	173	if (aCorrFctn.fDPhiDEtaDenominator)
	174	fDPhiDEtaDenominator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
	175	else
	176	fDPhiDEtaDenominator = 0;
	177
	178	if (aCorrFctn.fDPhiNumerator)
	179	fDPhiNumerator = new TH1D(*aCorrFctn.fDPhiNumerator);
	180	else
	181	fDPhiNumerator = 0;
	182	if (aCorrFctn.fDPhiDenominator)
	183	fDPhiDenominator = new TH1D(*aCorrFctn.fDPhiDenominator);
	184	else
	185	fDPhiDenominator = 0;
	186
	187	if (aCorrFctn.fDCosNumerator)
	188	fDCosNumerator = new TH1D(*aCorrFctn.fDCosNumerator);
	189	else
	190	fDCosNumerator = 0;
	191	if (aCorrFctn.fDCosDenominator)
	192	fDCosDenominator = new TH1D(*aCorrFctn.fDCosDenominator);
	193	else
	194	fDCosDenominator = 0;
	195
	196	if (aCorrFctn.fDPhiPtNumerator)
	197	fDPhiPtNumerator = new TH2D(*aCorrFctn.fDPhiPtNumerator);
	198	else
	199	fDPhiPtNumerator = 0;
	200	if (aCorrFctn.fDPhiPtDenominator)
	201	fDPhiPtDenominator = new TH2D(*aCorrFctn.fDPhiPtDenominator);
	202	else
	203	fDPhiPtDenominator = 0;
	204
	205	if (aCorrFctn.fDCosPtNumerator)
	206	fDCosPtNumerator = new TH2D(*aCorrFctn.fDCosPtNumerator);
	207	else
	208	fDCosPtNumerator = 0;
	209	if (aCorrFctn.fDCosPtDenominator)
	210	fDCosPtDenominator = new TH2D(*aCorrFctn.fDCosPtDenominator);
	211	else
	212	fDCosPtDenominator = 0;
973a91f8	213	if (aCorrFctn.fPhi)
	214	fPhi = new TH1D(*aCorrFctn.fPhi);
	215	else
	216	fPhi = 0;
	217	if (aCorrFctn.fEta)
	218	fEta = new TH1D(*aCorrFctn.fEta);
	219	else
	220	fEta = 0;
76ce4b5b	221
3b97e69c	222	if (aCorrFctn.fYtYtNumerator)
	223	fYtYtNumerator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
	224	else
	225	fYtYtNumerator = 0;
	226
	227	if (aCorrFctn.fYtYtDenominator)
	228	fYtYtDenominator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
	229	else
	230	fYtYtDenominator = 0;
d381642d	231
4d019bab	232	fphiL = aCorrFctn.fphiL;
	233	fphiT = aCorrFctn.fphiT;
	234
d381642d	235	// if (aCorrFctn.fPtCorrectionsNum)
	236	// fPtCorrectionsNum = new THnSparseF(*aCorrFctn.fPtCorrectionsNum);
	237	// else
	238	// fPtCorrectionsNum = 0;
	239
	240	// if (aCorrFctn.fPtCorrectionsDen)
	241	// fPtCorrectionsDen = new THnSparseF(*aCorrFctn.fPtCorrectionsDen);
	242	// else
	243	// fPtCorrectionsDen = 0;
	244
	245
	246
76ce4b5b	247	}
	248	//____________________________
	249	AliFemtoCorrFctnDEtaDPhi::~AliFemtoCorrFctnDEtaDPhi(){
	250	// destructor
	251	delete fDPhiDEtaNumerator;
	252	delete fDPhiDEtaDenominator;
	253	delete fDPhiNumerator;
	254	delete fDPhiDenominator;
	255	delete fDCosNumerator;
	256	delete fDCosDenominator;
	257	if (fDoPtAnalysis) {
	258	delete fDPhiPtNumerator;
	259	delete fDPhiPtDenominator;
	260	delete fDCosPtNumerator;
	261	delete fDCosPtDenominator;
	262	}
973a91f8	263	delete fPhi;
973a91f8	264	delete fEta;
3b97e69c	265
	266	delete fYtYtNumerator;
	267	delete fYtYtDenominator;
d381642d	268
	269	delete fPtCorrectionsNum;
	270	delete fPtCorrectionsDen;
	271	delete fEtaCorrectionsNum;
	272	delete fEtaCorrectionsDen;
76ce4b5b	273	}
	274	//_________________________
	275	AliFemtoCorrFctnDEtaDPhi& AliFemtoCorrFctnDEtaDPhi::operator=(const AliFemtoCorrFctnDEtaDPhi& aCorrFctn)
	276	{
	277	// assignment operator
	278	if (this == &aCorrFctn)
	279	return *this;
	280
	281	if (aCorrFctn.fDPhiDEtaNumerator)
	282	fDPhiDEtaNumerator = new TH2D(*aCorrFctn.fDPhiDEtaNumerator);
	283	else
	284	fDPhiDEtaNumerator = 0;
	285	if (aCorrFctn.fDPhiDEtaDenominator)
	286	fDPhiDEtaDenominator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
	287	else
	288	fDPhiDEtaDenominator = 0;
	289
	290	if (aCorrFctn.fDPhiNumerator)
	291	fDPhiNumerator = new TH1D(*aCorrFctn.fDPhiNumerator);
	292	else
	293	fDPhiNumerator = 0;
	294	if (aCorrFctn.fDPhiDenominator)
	295	fDPhiDenominator = new TH1D(*aCorrFctn.fDPhiDenominator);
	296	else
	297	fDPhiDenominator = 0;
	298
	299	if (aCorrFctn.fDCosNumerator)
	300	fDCosNumerator = new TH1D(*aCorrFctn.fDCosNumerator);
	301	else
	302	fDCosNumerator = 0;
	303	if (aCorrFctn.fDCosDenominator)
	304	fDCosDenominator = new TH1D(*aCorrFctn.fDCosDenominator);
	305	else
	306	fDCosDenominator = 0;
	307
	308	if (aCorrFctn.fDPhiPtNumerator)
	309	fDPhiPtNumerator = new TH2D(*aCorrFctn.fDPhiPtNumerator);
	310	else
	311	fDPhiPtNumerator = 0;
	312	if (aCorrFctn.fDPhiPtDenominator)
	313	fDPhiPtDenominator = new TH2D(*aCorrFctn.fDPhiPtDenominator);
	314	else
	315	fDPhiPtDenominator = 0;
	316
	317	if (aCorrFctn.fDCosPtNumerator)
	318	fDCosPtNumerator = new TH2D(*aCorrFctn.fDCosPtNumerator);
	319	else
	320	fDCosPtNumerator = 0;
	321	if (aCorrFctn.fDCosPtDenominator)
	322	fDCosPtDenominator = new TH2D(*aCorrFctn.fDCosPtDenominator);
	323	else
	324	fDCosPtDenominator = 0;
973a91f8	325	if (aCorrFctn.fPhi)
	326	fPhi = new TH1D(*aCorrFctn.fPhi);
	327	else
	328	fPhi = 0;
	329	if (aCorrFctn.fEta)
	330	fEta = new TH1D(*aCorrFctn.fEta);
	331	else
	332	fEta = 0;
76ce4b5b	333
3b97e69c	334	if (aCorrFctn.fYtYtNumerator)
	335	fYtYtNumerator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
	336	else
	337	fYtYtNumerator = 0;
	338
	339	if (aCorrFctn.fYtYtDenominator)
	340	fYtYtDenominator = new TH2D(*aCorrFctn.fDPhiDEtaDenominator);
	341	else
	342	fYtYtDenominator = 0;
	343
41fa2a11	344	fIfCorrectionHist = kNone;
d381642d	345
4d019bab	346	fphiL = aCorrFctn.fphiL;
	347	fphiT = aCorrFctn.fphiT;
	348
d381642d	349	// if (aCorrFctn.fPtCorrectionsNum)
	350	// fPtCorrectionsNum = new THnSparseF(*aCorrFctn.fPtCorrectionsNum);
	351	// else
	352	// fPtCorrectionsNum = 0;
	353
	354	// if (aCorrFctn.fPtCorrectionsDen)
	355	// fPtCorrectionsDen = new THnSparseF(*aCorrFctn.fPtCorrectionsDen);
	356	// else
	357	// fPtCorrectionsDen = 0;
	358
	359
	360
76ce4b5b	361	return *this;
	362	}
	363	//_________________________
	364	void AliFemtoCorrFctnDEtaDPhi::Finish(){
	365	// here is where we should normalize, fit, etc...
	366	// we should NOT Draw() the histos (as I had done it below),
	367	// since we want to insulate ourselves from root at this level
	368	// of the code. Do it instead at root command line with browser.
	369	// mShareNumerator->Draw();
	370	//mShareDenominator->Draw();
	371	//mRatio->Draw();
	372
	373	}
	374
	375	//____________________________
	376	AliFemtoString AliFemtoCorrFctnDEtaDPhi::Report(){
	377	// create report
	378	string stemp = "TPC Ncls Correlation Function Report:\n";
	379	char ctemp[100];
	380	snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fDPhiDEtaNumerator->GetEntries());
	381	stemp += ctemp;
	382	snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDPhiDEtaDenominator->GetEntries());
	383	stemp += ctemp;
	384	// stemp += mCoulombWeight->Report();
	385	AliFemtoString returnThis = stemp;
	386	return returnThis;
	387	}
	388	//____________________________
	389	void AliFemtoCorrFctnDEtaDPhi::AddRealPair( AliFemtoPair* pair){
	390	// add real (effect) pair
	391	if (fPairCut)
	392	if (!fPairCut->Pass(pair)) return;
	393
973a91f8	394	/*double phi1 = pair->Track1()->Track()->P().Phi();
76ce4b5b	395	double phi2 = pair->Track2()->Track()->P().Phi();
76ce4b5b	396	double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
973a91f8	397	double eta2 = pair->Track2()->Track()->P().PseudoRapidity();*/
	398
	399	double phi1 = pair->Track1()->FourMomentum().Phi();
	400	double phi2 = pair->Track2()->FourMomentum().Phi();
	401	double eta1 = pair->Track1()->FourMomentum().PseudoRapidity();
	402	double eta2 = pair->Track2()->FourMomentum().PseudoRapidity();
76ce4b5b	403
76ce4b5b	404	double dphi = phi1 - phi2;
4d019bab	405	while (dphi<fphiL) dphi+=PIT;
4d019bab	406	while (dphi>fphiT) dphi-=PIT;
76ce4b5b	407
	408	double deta = eta1 - eta2;
	409
d381642d	410	double px1 = pair->Track1()->Track()->P().x();
	411	double py1 = pair->Track1()->Track()->P().y();
	412	//double pz1 = pair->Track1()->Track()->P().z();
76ce4b5b	413
d381642d	414	double px2 = pair->Track2()->Track()->P().x();
	415	double py2 = pair->Track2()->Track()->P().y();
	416	//double pz2 = pair->Track2()->Track()->P().z();
76ce4b5b	417
d381642d	418	double pt1 = TMath::Hypot(px1, py1);
d381642d	419	double pt2 = TMath::Hypot(px2, py2);
76ce4b5b	420	// double ptmin = pt1>pt2 ? pt2 : pt1;
	421
	422	// double cosphi = (px1px2 + py1py2 + pz1*pz2)/
	423	// sqrt((px1px1 + py1py1 + pz1pz1)(px2px2 + py2py2 + pz2*pz2));
	424
	425	fDPhiDEtaNumerator->Fill(dphi, deta);
	426
	427	fDPhiNumerator->Fill(dphi);
	428	// fDCosNumerator->Fill(cosphi);
	429
	430	if (fDoPtAnalysis) {
	431	// fDPhiPtNumerator->Fill(dphi, ptmin);
	432	// fDCosPtNumerator->Fill(cosphi, ptmin);
	433	}
	434
973a91f8	435	fPhi->Fill(phi1);
	436	fEta->Fill(eta1);
	437
3b97e69c	438	double PionMass = 0.13956995;
d381642d	439	double yt1 = TMath::Log(sqrt(1+(pt1/PionMass)*(pt1/PionMass))+(pt1/PionMass));
d381642d	440	double yt2 = TMath::Log(sqrt(1+(pt2/PionMass)*(pt2/PionMass))+(pt2/PionMass));
3b97e69c	441	fYtYtNumerator->Fill(yt1,yt2);
3b97e69c	442
41fa2a11	443	if(fIfCorrectionHist)
d381642d	444	{
41fa2a11	445	if(fIfCorrectionHist == kPt){
d381642d	446	Double_t val[] = {pt1,pt2,dphi,deta};
	447	fPtCorrectionsNum->Fill(val);
	448	}
41fa2a11	449	if(fIfCorrectionHist == kEta){
d381642d	450	Double_t val[] = {eta1,eta2,dphi,deta};
	451	fEtaCorrectionsNum->Fill(val);
	452	}
	453
	454	}
	455
76ce4b5b	456	}
	457	//____________________________
	458	void AliFemtoCorrFctnDEtaDPhi::AddMixedPair( AliFemtoPair* pair){
	459	// add mixed (background) pair
	460	if (fPairCut)
	461	if (!fPairCut->Pass(pair)) return;
	462
973a91f8	463	/*double phi1 = pair->Track1()->Track()->P().Phi();
76ce4b5b	464	double phi2 = pair->Track2()->Track()->P().Phi();
76ce4b5b	465	double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
973a91f8	466	double eta2 = pair->Track2()->Track()->P().PseudoRapidity();*/
	467
	468	double phi1 = pair->Track1()->FourMomentum().Phi();
	469	double phi2 = pair->Track2()->FourMomentum().Phi();
	470	double eta1 = pair->Track1()->FourMomentum().PseudoRapidity();
	471	double eta2 = pair->Track2()->FourMomentum().PseudoRapidity();
76ce4b5b	472
76ce4b5b	473	double dphi = phi1 - phi2;
4d019bab	474	while (dphi<fphiL) dphi+=PIT;
4d019bab	475	while (dphi>fphiT) dphi-=PIT;
76ce4b5b	476
	477	double deta = eta1 - eta2;
	478
d381642d	479	double px1 = pair->Track1()->Track()->P().x();
	480	double py1 = pair->Track1()->Track()->P().y();
	481	//double pz1 = pair->Track1()->Track()->P().z();
76ce4b5b	482
d381642d	483	double px2 = pair->Track2()->Track()->P().x();
	484	double py2 = pair->Track2()->Track()->P().y();
	485	//double pz2 = pair->Track2()->Track()->P().z();
76ce4b5b	486
d381642d	487	double pt1 = TMath::Hypot(px1, py1);
d381642d	488	double pt2 = TMath::Hypot(px2, py2);
76ce4b5b	489	// double ptmin = pt1>pt2 ? pt2 : pt1;
	490
	491	// double cosphi = (px1px2 + py1py2 + pz1*pz2)/
	492	// sqrt((px1px1 + py1py1 + pz1pz1)(px2px2 + py2py2 + pz2*pz2));
	493
	494	fDPhiDEtaDenominator->Fill(dphi, deta);
	495
	496	fDPhiDenominator->Fill(dphi);
	497	// fDCosDenominator->Fill(cosphi);
	498
3b97e69c	499	//if (fDoPtAnalysis) {
76ce4b5b	500	// fDPhiPtDenominator->Fill(dphi, ptmin);
76ce4b5b	501	// fDCosPtDenominator->Fill(cosphi, ptmin);
3b97e69c	502	//}
	503
	504	double PionMass = 0.13956995;
d381642d	505	double yt1 = TMath::Log(sqrt(1+(pt1/PionMass)*(pt1/PionMass))+(pt1/PionMass));
	506	double yt2 = TMath::Log(sqrt(1+(pt2/PionMass)*(pt2/PionMass))+(pt2/PionMass));
	507	fYtYtDenominator->Fill(yt1,yt2);
	508
41fa2a11	509	if(fIfCorrectionHist)
d381642d	510	{
41fa2a11	511	if(fIfCorrectionHist == kPt){
d381642d	512	Double_t val[] = {pt1,pt2,dphi,deta};
	513	fPtCorrectionsDen->Fill(val);
	514	}
41fa2a11	515	if(fIfCorrectionHist == kEta){
d381642d	516	Double_t val[] = {eta1,eta2,dphi,deta};
	517	fEtaCorrectionsDen->Fill(val);
	518	}
	519	}
	520
76ce4b5b	521	}
	522
	523
	524	void AliFemtoCorrFctnDEtaDPhi::WriteHistos()
	525	{
	526	// Write out result histograms
	527	fDPhiDEtaNumerator->Write();
	528	fDPhiDEtaDenominator->Write();
973a91f8	529	/*fDPhiNumerator->Write();
76ce4b5b	530	fDPhiDenominator->Write();
76ce4b5b	531	fDCosNumerator->Write();
41fa2a11	532	fDCosDenominator->Write();*/
76ce4b5b	533	if (fDoPtAnalysis) {
	534	fDPhiPtNumerator->Write();
	535	fDPhiPtDenominator->Write();
	536	fDCosPtNumerator->Write();
	537	fDCosPtDenominator->Write();
41fa2a11	538	}
973a91f8	539	fPhi->Write();
973a91f8	540	fEta->Write();
3b97e69c	541
41fa2a11	542	if(fIfCorrectionHist){
41fa2a11	543	if(fIfCorrectionHist==kPt){
d381642d	544	fPtCorrectionsNum->Write();
d381642d	545	fPtCorrectionsDen->Write();}
41fa2a11	546	if(fIfCorrectionHist==kEta){
d381642d	547	fEtaCorrectionsNum->Write();
	548	fEtaCorrectionsDen->Write();}
	549	}
76ce4b5b	550	}
	551
	552	TList* AliFemtoCorrFctnDEtaDPhi::GetOutputList()
	553	{
	554	// Prepare the list of objects to be written to the output
	555	TList *tOutputList = new TList();
	556
	557	tOutputList->Add(fDPhiDEtaNumerator);
	558	tOutputList->Add(fDPhiDEtaDenominator);
973a91f8	559	/*tOutputList->Add(fDPhiNumerator);
76ce4b5b	560	tOutputList->Add(fDPhiDenominator);
	561	tOutputList->Add(fDCosNumerator);
	562	tOutputList->Add(fDCosDenominator);
	563	if (fDoPtAnalysis) {
	564	tOutputList->Add(fDPhiPtNumerator);
	565	tOutputList->Add(fDPhiPtDenominator);
	566	tOutputList->Add(fDCosPtNumerator);
	567	tOutputList->Add(fDCosPtDenominator);
973a91f8	568	}*/
	569	tOutputList->Add(fPhi);
	570	tOutputList->Add(fEta);
3b97e69c	571	tOutputList->Add(fYtYtNumerator);
3b97e69c	572	tOutputList->Add(fYtYtDenominator);
d381642d	573
41fa2a11	574	if(fIfCorrectionHist){
41fa2a11	575	if(fIfCorrectionHist==kPt){
d381642d	576	tOutputList->Add(fPtCorrectionsNum);
	577	tOutputList->Add(fPtCorrectionsDen);
	578	}
41fa2a11	579	if(fIfCorrectionHist==kEta){
d381642d	580	tOutputList->Add(fEtaCorrectionsNum);
	581	tOutputList->Add(fEtaCorrectionsDen);
	582	}
	583	}
76ce4b5b	584	return tOutputList;
	585
	586	}
	587
	588	void AliFemtoCorrFctnDEtaDPhi::SetDoPtAnalysis(int do2d)
	589	{
	590	fDoPtAnalysis = do2d;
	591
	592	int aPhiBins = fDPhiDEtaNumerator->GetNbinsX();
	593	const char *title = fDPhiDEtaNumerator->GetTitle();
	594
	595	// set up numerator
	596	char tTitNumDPhiPt[101] = "NumDPhiPt";
	597	strncat(tTitNumDPhiPt,title, 100);
4d019bab	598	fDPhiPtNumerator = new TH2D(tTitNumDPhiPt,title,aPhiBins2,-0.5TMath::Pi(),3./2.*TMath::Pi(), 30, 0.0, 3.0);
76ce4b5b	599	// set up denominator
	600	char tTitDenDPhiPt[101] = "DenDPhiPt";
	601	strncat(tTitDenDPhiPt,title, 100);
4d019bab	602	fDPhiPtDenominator = new TH2D(tTitDenDPhiPt,title,aPhiBins2,-0.5TMath::Pi(),3./2.*TMath::Pi(), 30, 0.0, 3.0);
76ce4b5b	603
	604	// set up numerator
	605	char tTitNumDCosPt[101] = "NumDCosPt";
	606	strncat(tTitNumDCosPt,title, 100);
	607	fDCosPtNumerator = new TH2D(tTitNumDCosPt,title,aPhiBins*2,-1.0,1.0, 30, 0.0, 3.0);
	608	// set up denominator
	609	char tTitDenDCosPt[101] = "DenDCosPt";
	610	strncat(tTitDenDCosPt,title, 100);
	611	fDCosPtDenominator = new TH2D(tTitDenDCosPt,title,aPhiBins*2,-1.0,1.0, 30, 0.0, 3.0);
	612
	613	fDPhiPtNumerator->Sumw2();
	614	fDPhiPtDenominator->Sumw2();
	615	fDCosPtNumerator->Sumw2();
	616	fDCosPtDenominator->Sumw2();
	617
	618	}
d381642d	619
41fa2a11	620	void AliFemtoCorrFctnDEtaDPhi::SetDo4DCorrectionHist(CorrectionType doCorr)
d381642d	621	{
41fa2a11	622	fIfCorrectionHist = doCorr;
d381642d	623	}