[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliCumulantsFunctions.cxx

//******************************* 
// flow analysis with cumulants *   
// author: Ante Bilandzic       * 
// email: anteb@nikhef.nl       *
//******************************* 

#define AliCumulantsFunctions_cxx

#include "Riostream.h"
#include "AliFlowCommonConstants.h"
#include "TChain.h"
#include "TFile.h"
#include "TList.h" //NEW
#include "TParticle.h"

#include "TProfile.h"
#include "TProfile2D.h" 
#include "TProfile3D.h"

#include "AliFlowEventSimple.h"
#include "AliFlowTrackSimple.h"
#include "AliFlowAnalysisWithCumulants.h"
#include "AliFlowCumuConstants.h"
#include "AliFlowCommonConstants.h"
#include "AliCumulantsFunctions.h"


ClassImp(AliCumulantsFunctions)

//________________________________________________________________________

AliCumulantsFunctions::AliCumulantsFunctions():  
  fIG(NULL),
  fDGRe(NULL),
  fDGIm(NULL),
  fifr(NULL),
  fdfr2(NULL), 
  fdfr4(NULL), 
  fdfr6(NULL),
  fdfr8(NULL),
  fAvMult(0)
 {
   //default constructor 
 }

AliCumulantsFunctions::~AliCumulantsFunctions(){
  //desctructor
}

AliCumulantsFunctions::AliCumulantsFunctions(TProfile2D *IntGen, TProfile3D *DiffGenRe, TProfile3D *DiffGenIm, TH1D *ifr, TH1D *dfr2, TH1D *dfr4, TH1D *dfr6, TH1D *dfr8, Double_t CvM):
  fIG(IntGen),
  fDGRe(DiffGenRe),
  fDGIm(DiffGenIm),
  fifr(ifr),
  fdfr2(dfr2), 
  fdfr4(dfr4), 
  fdfr6(dfr6),
  fdfr8(dfr8),
  fAvMult(CvM)
 {
   //custom constructor 
 }
  
//___________________________________________________________________________
void AliCumulantsFunctions::Calculate(){
 //calculate final flow estimates
 
     static const Int_t fgkQmax=AliFlowCumuConstants::kQmax;//needed for numerics
     static const Int_t fgkPmax=AliFlowCumuConstants::kPmax;//needed for numerics  
     static const Int_t fgkFlow=AliFlowCumuConstants::kFlow;//integrated flow coefficient to be calculated
     static const Int_t fgkMltpl=AliFlowCumuConstants::kMltpl;//the multiple in p=m*n (diff. flow) 
     static const Int_t fgknBins=100;//number of pt bins //to be improved
     Double_t fR0=AliFlowCumuConstants::fgR0;//needed for numerics
     Double_t fPtMax=AliFlowCommonConstants::GetPtMax();//maximum pt
     Double_t fPtMin=AliFlowCommonConstants::GetPtMin();//minimum pt
     Double_t fBinWidth=(fPtMax-fPtMin)/fgknBins;//width of pt bin (in GeV)   
     
  Double_t fBvG[fgkPmax][fgkQmax]={0.};
        
  for(Int_t p=0;p<fgkPmax;p++){
   for(Int_t q=0;q<fgkQmax;q++){ 
    fBvG[p][q]=fIG->GetBinContent(p+1,q+1);
   }
  } 
      
 
  /////////////////////////////////////////////////////////////////////////////      
  //////////////////gen. function for the cumulants////////////////////////////
  /////////////////////////////////////////////////////////////////////////////
  
  Double_t fC[fgkPmax][fgkQmax]={0.};
  for (Int_t p=0;p<fgkPmax;p++){
    for (Int_t q=0;q<fgkQmax;q++){
      fC[p][q]=1.*fAvMult*(pow(fBvG[p][q],(1./fAvMult))-1.); //to be improved
    }
  }
 
  /////////////////////////////////////////////////////////////////////////////
  ///////avaraging the gen. function for the cumulants over azimuth////////////
  /////////////////////////////////////////////////////////////////////////////
  
  Double_t fCAv[fgkPmax]={0.}; 
  for (Int_t p=0;p<fgkPmax;p++){
    Double_t fTempHere=0.; 
    for (Int_t q=0;q<fgkQmax;q++){
      fTempHere+=1.*fC[p][q];
    } 
    fCAv[p]=1.*fTempHere/fgkQmax;
  }
  
  /////////////////////////////////////////////////////////////////////////////
  //////////////////////////////////final results//////////////////////////////
  /////////////////////////////////////////////////////////////////////////////
  
  Double_t fCumulant[fgkPmax];//c_{iFlow}\{2(p+1)\}
  
  fCumulant[0]=(1./(fR0*fR0)) * (8.*fCAv[0] - 14.*fCAv[1] + (56./3.)*fCAv[2] - (35./2.)*fCAv[3] + 
			      (56./5.)*fCAv[4] - (14./3.)*fCAv[5] + (8./7.)*fCAv[6] - (1./8.)*fCAv[7]);

  fCumulant[1]=(1./pow(fR0,4.)) * ((-1924./35.)*fCAv[0] + (621./5.)*fCAv[1] - (8012./45.)*fCAv[2] + 
				 (691./4.)*fCAv[3] - (564./5.)*fCAv[4] + (2143./45.)*fCAv[5] - 
				 (412./35.)*fCAv[6] + (363./280.)*fCAv[7]);

  fCumulant[2]=(1./pow(fR0,6.)) * (349.*fCAv[0] - (18353./20.)*fCAv[1] + (7173./5.)*fCAv[2] - 
				 1457.*fCAv[3] + (4891./5.)*fCAv[4] - (1683./4.)*fCAv[5] + 
				 (527./5.)*fCAv[6] - (469./40.)*fCAv[7]);

  fCumulant[3]=(1./pow(fR0,8.)) * ((-10528./5.)*fCAv[0] + (30578./5.)*fCAv[1] - (51456./5.)*fCAv[2] + 
				 10993.*fCAv[3] - (38176./5.)*fCAv[4] + (16818./5.)*fCAv[5] - 
				 (4288./5.)*fCAv[6] + (967./10.)*fCAv[7]);

  fCumulant[4]=(1./pow(fR0,10.)) * (11500.*fCAv[0] - 35800.*fCAv[1] + 63900.*fCAv[2] - 71600.*fCAv[3] + 
				  51620.*fCAv[4] - 23400.*fCAv[5] + 6100.*fCAv[6] - 700.*fCAv[7]);

  fCumulant[5]=(1./pow(fR0,12.)) * (-52560.*fCAv[0] + 172080.*fCAv[1] - 321840.*fCAv[2] + 376200.*fCAv[3] - 
				  281520.*fCAv[4] + 131760.*fCAv[5] - 35280.*fCAv[6] + 4140.*fCAv[7]);

  fCumulant[6]=(1./pow(fR0,14.)) * (176400.*fCAv[0] - 599760.*fCAv[1] + 1164240.*fCAv[2] - 1411200.*fCAv[3] + 
				  1093680.*fCAv[4] - 529200.*fCAv[5] + 146160.*fCAv[6] - 17640.*fCAv[7]);

  fCumulant[7]=(1./pow(fR0,16.)) * (-322560*fCAv[0] + 1128960.*fCAv[1] - 2257920.*fCAv[2] + 2822400.*fCAv[3] - 
				  2257920.*fCAv[4] + 1128960.*fCAv[5] - 322560.*fCAv[6] + 40320.*fCAv[7]);
  
  
  cout<<""<<endl;
  cout<<"***************************"<<endl;
  cout<<"cumulants:"<<endl;
  
  cout<<" c_"<<fgkFlow<<"{2} = "<<fCumulant[0]<<endl; 
  cout<<" c_"<<fgkFlow<<"{4} = "<<fCumulant[1]<<endl;
  cout<<" c_"<<fgkFlow<<"{6} = "<<fCumulant[2]<<endl;
  cout<<" c_"<<fgkFlow<<"{8} = "<<fCumulant[3]<<endl; 
  cout<<"c_"<<fgkFlow<<"{10} = "<<fCumulant[4]<<endl; 
  cout<<"c_"<<fgkFlow<<"{12} = "<<fCumulant[5]<<endl;
  cout<<"c_"<<fgkFlow<<"{14} = "<<fCumulant[6]<<endl; 
  cout<<"c_"<<fgkFlow<<"{16} = "<<fCumulant[7]<<endl; 
  
  cout<<""<<endl;
  cout<<"integrated flow: "<<endl;
  
  
  Double_t fV2=0.,fV4=0.,fV6=0.,fV8=0.;
  Double_t fSdQ[4]={0.};
  Double_t fChiQ[4]={0.};
   if (fCumulant[0]>=0.){ 
    fV2=sqrt(fCumulant[0]);    
    //if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV2*fAvM,2.)>0.){
     //fChiQ[0]=fAvM*fV2/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV2*fAvM,2.),0.5);
     //fSdQ[0]=pow(((1./(2.*fAvM*nEvents))*((1.+1.*pow(fChiQ[0],2))/(1.*pow(fChiQ[0],2)))),0.5);
     cout<<" v_"<<fgkFlow<<"{2} = "<<100.*fV2<<"%, chi{2} = "<<fChiQ[0]<<", sd{2} = "<<100.*fSdQ[0]<<"%"<<endl;//to be improved (2->fgkFlow)
     //fCommonHistsRes2->FillIntegratedFlow(100.*fV2,100.*fSdQ[0]);
     //fCommonHistsRes2->FillChi(fChiQ[0]);
     fifr->SetBinContent(1,100.*fV2);
    //}
   //} else {
    //cout<<" v_"<<fgkFlow<<"{2} = Im"<<endl;  
   } else {
    //cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl; 
   }
  if (fCumulant[1]<=0.){
    fV4=pow(-fCumulant[1],(1./4.));
    //if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV4*fAvM,2.)>0.){
     //fChiQ[1]=fAvM*fV4/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV4*fAvM,2.),0.5);
     //fSdQ[1]=(1./(pow(2.*fAvM*nEvents,.5)))*pow((1.+2.*pow(fChiQ[1],2)+(1./4.)*pow(fChiQ[1],4.)+(1./4.)*pow(fChiQ[1],6.))/((1./4.)*pow(fChiQ[1],6.)),.5);
     cout<<" v_"<<fgkFlow<<"{4} = "<<100.*fV4<<"%, chi{4} = "<<fChiQ[1]<<", sd{4} = "<<100.*fSdQ[1]<<"%"<<endl;//to be improved (2->fgkFlow)
     //fCommonHistsRes4->FillInteg8ratedFlow(100.*fV4,100.*fSdQ[1]);
     //fCommonHistsRes4->FillChi(fChiQ[1]);
     fifr->SetBinContent(2,100.*fV4);
    //} else {
      //cout<<" v_"<<fgkFlow<<"{4} = Im"<<endl;
    //}
 // } else {
   // cout<<" v_"<<fgkFlow<<"{4} = Im"<<endl;  
  } 
  if (fCumulant[2]>=0.){
    fV6=pow((1./4.)*fCumulant[2],(1./6.));
    //if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV6*fAvM,2.)>0.){
     //fChiQ[2]=fAvM*fV6/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV6*fAvM,2.),0.5);
     //fSdQ[2]=(1./(pow(2.*fAvM*nEvents,.5)))*pow((3.+18.*pow(fChiQ[2],2)+9.*pow(fChiQ[2],4.)+28.*pow(fChiQ[2],6.)+12.*pow(fChiQ[2],8.)+24.*pow(fChiQ[2],10.))/(24.*pow(fChiQ[2],10.)),.5);
     cout<<" v_"<<fgkFlow<<"{6} = "<<100.*fV6<<"%, chi{6} = "<<fChiQ[2]<<", sd{6} = "<<100.*fSdQ[2]<<"%"<<endl;//to be improved (2->fgkFlow)
     //fCommonHistsRes6->FillIntegratedFlow(100.*fV6,100.*fSdQ[2]);
     //fCommonHistsRes6->FillChi(fChiQ[2]);
     fifr->SetBinContent(3,100.*fV6);
    //} else {
     // cout<<" v_"<<fgkFlow<<"{6} = Im"<<endl; 
   // }
  } else {
    //cout<<" v_"<<fgkFlow<<"{6} = Im"<<endl;  
  }
  if (fCumulant[3]<=0.){
    fV8=pow(-(1./33.)*fCumulant[3],(1./8.));
    //if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV8*fAvM,2.)>0.){
     //fChiQ[3]=fAvM*fV8/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV8*fAvM,2.),0.5);
     //fSdQ[3]=(1./(pow(2.*fAvM*nEvents,.5)))*pow((12.+96.*pow(fChiQ[3],2)+72.*pow(fChiQ[3],4.)+304.*pow(fChiQ[3],6.)+257.*pow(fChiQ[3],8.)+804.*pow(fChiQ[3],10.)+363.*pow(fChiQ[3],12.)+726.*pow(fChiQ[3],14.))/(726.*pow(fChiQ[3],14.)),.5);
      cout<<" v_"<<fgkFlow<<"{8} = "<<100.*fV8<<"%, chi{8} = "<<fChiQ[3]<<", sd{8} = "<<100.*fSdQ[3]<<"%"<<endl;//to be improved (2->fgkFlow)
      //fCommonHistsRes8->FillIntegratedFlow(100.*fV8,100.*fSdQ[3]);
      //fCommonHistsRes8->FillChi(fChiQ[3]);
      fifr->SetBinContent(4,100.*fV8);
     //} else {
       //cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl;a
     //}
  } else {
    //cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl;     
  }
  if (fCumulant[4]>=0.){//fHistProIntFlow
    cout<<"v_"<<fgkFlow<<"{10} = "<<100.*pow((1./456.)*fCumulant[4],(1./10.))<<"%"<<endl;//to be improved (2->fgkFlow)
    fifr->SetBinContent(5,100.*pow((1./456.)*fCumulant[4],(1./10.)));
  } else {
      cout<<"v_"<<fgkFlow<<"{10} = Im"<<endl;  //to be improved (2->fgkFlow)
  }
  if (fCumulant[5]<=0.){
    cout<<"v_"<<fgkFlow<<"{12} = "<<100.*pow(-(1./9460.)*fCumulant[5],(1./12.))<<"%"<<endl;//to be improved (2->fgkFlow)
    fifr->SetBinContent(6,100.*pow(-(1./9460.)*fCumulant[5],(1./12.)));
  } else {
    cout<<"v_"<<fgkFlow<<"{12} = Im"<<endl;  //to be improved (2->fgkFlow)
  }
  if (fCumulant[6]>=0.){
    cout<<"v_"<<fgkFlow<<"{14} = "<<100.*pow((1./274800.)*fCumulant[6],(1./14.))<<"%"<<endl;//to be improved (2->fgkFlow)
    fifr->SetBinContent(7,100.*pow((1./274800.)*fCumulant[6],(1./14.)));
  } else {
    cout<<"v_"<<fgkFlow<<"{14} = Im"<<endl;  //to be improved (2->fgkFlow)
  }
  if (fCumulant[7]<=0.){
    cout<<"v_"<<fgkFlow<<"{16} = "<<100.*pow(-(1./10643745.)*fCumulant[7],(1./16.))<<"%"<<endl;//to be improved (2->fgkFlow)
    fifr->SetBinContent(8,100.*pow(-(1./10643745.)*fCumulant[7],(1./16.)));
  } else {
    cout<<"v_"<<fgkFlow<<"{16} = Im"<<endl;  //to be improved (2->fgkFlow)
  }
  //cout<<"***************************"<<endl;
      
  //DIFFERENTIAL FLOW CALCULATIONS STARTS HERE!!!
  
  Double_t fX[fgknBins][fgkPmax][fgkQmax]={0.};//see the text bellow relation (11) in PG
  Double_t fY[fgknBins][fgkPmax][fgkQmax]={0.};
  
  for(Int_t b=0;b<fgknBins;b++){
    for(Int_t p=0;p<fgkPmax;p++){
      for(Int_t q=0;q<fgkQmax;q++){
	fX[b][p][q]=fDGRe->GetBinContent(b+1,p+1,q+1)/fBvG[p][q];
	fY[b][p][q]=fDGIm->GetBinContent(b+1,p+1,q+1)/fBvG[p][q];
      }
    }   
  } 
  
 
  Double_t fD[fgknBins][fgkPmax]={0.};//implementing relation (11) from PG
  
  for (Int_t b=0;b<fgknBins;b++){
    for (Int_t p=0;p<fgkPmax;p++){
      Double_t fTempHere3=0.; 
      for (Int_t q=0;q<fgkQmax;q++){
	fTempHere3+=cos(fgkMltpl*2.*q*TMath::Pi()/fgkQmax)*fX[b][p][q] + sin(fgkMltpl*2.*q*TMath::Pi()/fgkQmax)*fY[b][p][q];
      } 
      fD[b][p]=1.*(pow(fR0*pow(p+1.,.5),fgkMltpl)/fgkQmax)*fTempHere3;
      if(fD[b][p]){
       //cout<<"And this "<<b<<" "<<fD[b][p]<<endl;
      }  
    }
  } 
  
  Double_t fDiffCumulant2[fgknBins]={0.};//implementing relation (12) from PG
  Double_t fDiffCumulant4[fgknBins]={0.};
  Double_t fDiffCumulant6[fgknBins]={0.};
  Double_t fDiffCumulant8[fgknBins]={0.};
  
  for (Int_t b=0;b<fgknBins;b++){
    fDiffCumulant2[b]=(1./(fR0*fR0))*(4.*fD[b][0]-3.*fD[b][1]+(4./3.)*fD[b][2]-(1./4.)*fD[b][3]);
    fDiffCumulant4[b]=(1./pow(fR0,4.))*((-26./3.)*fD[b][0]+(19./2.)*fD[b][1]-(14./3.)*fD[b][2]+(11./12.)*fD[b][3]);
    fDiffCumulant6[b]=(1./pow(fR0,6.))*(18.*fD[b][0]-24.*fD[b][1]+14.*fD[b][2]-3.*fD[b][3]);
    fDiffCumulant8[b]=(1./pow(fR0,8.))*(-24.*fD[b][0]+36.*fD[b][1]-24.*fD[b][2]+6.*fD[b][3]);
  }
  
  Double_t fv2[fgknBins],fv4[fgknBins],fv6[fgknBins],fv8[fgknBins];
  //Double_t fAvPt[fgknBins];
  Double_t fSddiff2[fgknBins],fSddiff4[fgknBins];

  cout<<"number of pt bins: "<<fgknBins<<endl;
  cout<<"****************************************"<<endl;
  for (Int_t b=0;b<fgknBins;b++){ 
    //if(fBinNoOfParticles[b]==0)continue;
    //fAvPt[b]=fBinMeanPt[b]/fBinNoOfParticles[b];
    cout<<"pt bin: "<<b*fBinWidth<<"-"<<(b+1)*fBinWidth<<" GeV"<<endl;
    //cout<<"number of particles in this pt bin: "<<tempNo->GetBinContent(b+1,1,1)<<endl;
    //cout<<"mean pt in this bin: "<<fAvPt[b]<<" GeV"<<endl;
    if(fCumulant[0]>=0){
      fv2[b]=100.*fDiffCumulant2[b]/pow(fCumulant[0],.5);
      //if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV2*fAvM,2.)>0.){
       //fSddiff2[b]=pow((1./(2.*fBinNoOfParticles[b]))*((1.+pow(fChiQ[0],2.))/pow(fChiQ[0],2.)),0.5);
       cout<<"v'_2/2{2} = "<<fv2[b]<<"%, "<<" "<<"sd{2} = "<<100.*fSddiff2[b]<<"%"<<endl;
       //fDiffFlowResults2->SetBinContent(b+1,fv2[b],100.*fSddiff2[b]);
       fdfr2->SetBinContent(b+1,fv2[b]);
      //} else {
         //cout<<"v'_2/2{2} = Im"<<endl;
      //}
    }else{
      cout<<"v'_2/2{2} = Im"<<endl;
    } 
    if(fCumulant[1]<=0){
      fv4[b]=-100.*fDiffCumulant4[b]/pow(-fCumulant[1],.75);
      //if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV4*fAvM,2.)>0.){
       //fSddiff4[b]=pow((1./(2.*fBinNoOfParticles[b]))*((2.+6.*pow(fChiQ[1],2.)+pow(fChiQ[1],4.)+pow(fChiQ[1],6.))/pow(fChiQ[1],6.)),0.5);
       cout<<"v'_2/2{4} = "<<fv4[b]<<"%, "<<" "<<"sd{4} = "<<100.*fSddiff4[b]<<"%"<<endl;
       //fCommonHistsRes4->FillDifferentialFlow(b+1,fv4[b],100.*fSddiff4[b]);
       fdfr4->SetBinContent(b+1,fv4[b]);
      //} else {
       // cout<<"v'_2/2{4} = Im"<<endl;
      //} 
    }else{
      cout<<"v'_2/2{4} = Im"<<endl;
    }  
    if(fCumulant[2]>=0){
      cout<<"v'_2/2{6} = "<<100.*fDiffCumulant6[b]/(4.*pow((1./4.)*fCumulant[2],(5./6.)))<<"%"<<endl;
      fv6[b]=100.*fDiffCumulant6[b]/(4.*pow((1./4.)*fCumulant[2],(5./6.)));
      //fCommonHistsRes6->FillDifferentialFlow(b+1,fv6[b],0.);
      fdfr6->SetBinContent(b+1,fv6[b]);
    }else{
      cout<<"v'_2/2{6} = Im"<<endl;
    }     
    if(fCumulant[3]<=0){
      cout<<"v'_2/2{8} = "<<-100.*fDiffCumulant8[b]/(33.*pow(-(1./33.)*fCumulant[3],(7./8.)))<<"%"<<endl;
      fv8[b]=-100.*fDiffCumulant8[b]/(33.*pow(-(1./33.)*fCumulant[3],(7./8.))); 
      //fCommonHistsRes8->FillDifferentialFlow(b+1,fv8[b],0.);
      fdfr8->SetBinContent(b+1,fv8[b]);
    }else{
      cout<<"v'_2/2{8} = Im"<<endl;
    }       
    cout<<"****************************************"<<endl;
  }  
}
Commit	Line	Data
aaebd73d	1	//*******************************
	2	// flow analysis with cumulants *
	3	// author: Ante Bilandzic *
	4	// email: anteb@nikhef.nl *
	5	//*******************************
	6
	7	#define AliCumulantsFunctions_cxx
	8
	9	#include "Riostream.h"
	10	#include "AliFlowCommonConstants.h"
	11	#include "TChain.h"
	12	#include "TFile.h"
	13	#include "TList.h" //NEW
	14	#include "TParticle.h"
	15
	16	#include "TProfile.h"
	17	#include "TProfile2D.h"
	18	#include "TProfile3D.h"
	19
	20	#include "AliFlowEventSimple.h"
	21	#include "AliFlowTrackSimple.h"
	22	#include "AliFlowAnalysisWithCumulants.h"
	23	#include "AliFlowCumuConstants.h"
	24	#include "AliFlowCommonConstants.h"
	25	#include "AliCumulantsFunctions.h"
	26
	27
	28	ClassImp(AliCumulantsFunctions)
	29
	30	//________________________________________________________________________
	31
	32	AliCumulantsFunctions::AliCumulantsFunctions():
	33	fIG(NULL),
	34	fDGRe(NULL),
	35	fDGIm(NULL),
	36	fifr(NULL),
	37	fdfr2(NULL),
	38	fdfr4(NULL),
	39	fdfr6(NULL),
	40	fdfr8(NULL),
	41	fAvMult(0)
	42	{
	43	//default constructor
	44	}
	45
	46	AliCumulantsFunctions::~AliCumulantsFunctions(){
	47	//desctructor
	48	}
	49
	50	AliCumulantsFunctions::AliCumulantsFunctions(TProfile2D IntGen, TProfile3D DiffGenRe, TProfile3D DiffGenIm, TH1D ifr, TH1D dfr2, TH1D dfr4, TH1D dfr6, TH1D dfr8, Double_t CvM):
	51	fIG(IntGen),
	52	fDGRe(DiffGenRe),
	53	fDGIm(DiffGenIm),
	54	fifr(ifr),
	55	fdfr2(dfr2),
	56	fdfr4(dfr4),
	57	fdfr6(dfr6),
	58	fdfr8(dfr8),
	59	fAvMult(CvM)
	60	{
	61	//custom constructor
82956bee	62	}
aaebd73d	63
	64	//___________________________________________________________________________
	65	void AliCumulantsFunctions::Calculate(){
	66	//calculate final flow estimates
	67
	68	static const Int_t fgkQmax=AliFlowCumuConstants::kQmax;//needed for numerics
	69	static const Int_t fgkPmax=AliFlowCumuConstants::kPmax;//needed for numerics
	70	static const Int_t fgkFlow=AliFlowCumuConstants::kFlow;//integrated flow coefficient to be calculated
	71	static const Int_t fgkMltpl=AliFlowCumuConstants::kMltpl;//the multiple in p=m*n (diff. flow)
	72	static const Int_t fgknBins=100;//number of pt bins //to be improved
	73	Double_t fR0=AliFlowCumuConstants::fgR0;//needed for numerics
	74	Double_t fPtMax=AliFlowCommonConstants::GetPtMax();//maximum pt
	75	Double_t fPtMin=AliFlowCommonConstants::GetPtMin();//minimum pt
	76	Double_t fBinWidth=(fPtMax-fPtMin)/fgknBins;//width of pt bin (in GeV)
	77
	78	Double_t fBvG[fgkPmax][fgkQmax]={0.};
	79
	80	for(Int_t p=0;p<fgkPmax;p++){
	81	for(Int_t q=0;q<fgkQmax;q++){
	82	fBvG[p][q]=fIG->GetBinContent(p+1,q+1);
	83	}
	84	}
	85
	86
	87	/////////////////////////////////////////////////////////////////////////////
	88	//////////////////gen. function for the cumulants////////////////////////////
	89	/////////////////////////////////////////////////////////////////////////////
	90
	91	Double_t fC[fgkPmax][fgkQmax]={0.};
	92	for (Int_t p=0;p<fgkPmax;p++){
	93	for (Int_t q=0;q<fgkQmax;q++){
	94	fC[p][q]=1.fAvMult(pow(fBvG[p][q],(1./fAvMult))-1.); //to be improved
	95	}
	96	}
	97
	98	/////////////////////////////////////////////////////////////////////////////
	99	///////avaraging the gen. function for the cumulants over azimuth////////////
	100	/////////////////////////////////////////////////////////////////////////////
	101
	102	Double_t fCAv[fgkPmax]={0.};
	103	for (Int_t p=0;p<fgkPmax;p++){
	104	Double_t fTempHere=0.;
	105	for (Int_t q=0;q<fgkQmax;q++){
	106	fTempHere+=1.*fC[p][q];
	107	}
	108	fCAv[p]=1.*fTempHere/fgkQmax;
	109	}
	110
	111	/////////////////////////////////////////////////////////////////////////////
	112	//////////////////////////////////final results//////////////////////////////
	113	/////////////////////////////////////////////////////////////////////////////
	114
	115	Double_t fCumulant[fgkPmax];//c_{iFlow}\{2(p+1)\}
	116
	117	fCumulant[0]=(1./(fR0fR0)) (8.fCAv[0] - 14.fCAv[1] + (56./3.)fCAv[2] - (35./2.)fCAv[3] +
	118	(56./5.)fCAv[4] - (14./3.)fCAv[5] + (8./7.)fCAv[6] - (1./8.)fCAv[7]);
	119
	120	fCumulant[1]=(1./pow(fR0,4.)) * ((-1924./35.)fCAv[0] + (621./5.)fCAv[1] - (8012./45.)*fCAv[2] +
	121	(691./4.)fCAv[3] - (564./5.)fCAv[4] + (2143./45.)*fCAv[5] -
	122	(412./35.)fCAv[6] + (363./280.)fCAv[7]);
	123
	124	fCumulant[2]=(1./pow(fR0,6.)) * (349.fCAv[0] - (18353./20.)fCAv[1] + (7173./5.)*fCAv[2] -
	125	1457.fCAv[3] + (4891./5.)fCAv[4] - (1683./4.)*fCAv[5] +
	126	(527./5.)fCAv[6] - (469./40.)fCAv[7]);
127
128	fCumulant[3]=(1./pow(fR0,8.)) * ((-10528./5.)fCAv[0] + (30578./5.)fCAv[1] - (51456./5.)*fCAv[2] +
129	10993.fCAv[3] - (38176./5.)fCAv[4] + (16818./5.)*fCAv[5] -
130	(4288./5.)fCAv[6] + (967./10.)fCAv[7]);
131
132	fCumulant[4]=(1./pow(fR0,10.)) * (11500.fCAv[0] - 35800.fCAv[1] + 63900.fCAv[2] - 71600.fCAv[3] +
133	51620.fCAv[4] - 23400.fCAv[5] + 6100.fCAv[6] - 700.fCAv[7]);
134
135	fCumulant[5]=(1./pow(fR0,12.)) * (-52560.fCAv[0] + 172080.fCAv[1] - 321840.fCAv[2] + 376200.fCAv[3] -
136	281520.fCAv[4] + 131760.fCAv[5] - 35280.fCAv[6] + 4140.fCAv[7]);
137
138	fCumulant[6]=(1./pow(fR0,14.)) * (176400.fCAv[0] - 599760.fCAv[1] + 1164240.fCAv[2] - 1411200.fCAv[3] +
139	1093680.fCAv[4] - 529200.fCAv[5] + 146160.fCAv[6] - 17640.fCAv[7]);
140
141	fCumulant[7]=(1./pow(fR0,16.)) * (-322560fCAv[0] + 1128960.fCAv[1] - 2257920.fCAv[2] + 2822400.fCAv[3] -
142	2257920.fCAv[4] + 1128960.fCAv[5] - 322560.fCAv[6] + 40320.fCAv[7]);
143
144
145	cout<<""<<endl;
146	cout<<"***************************"<<endl;
147	cout<<"cumulants:"<<endl;
148
149	cout<<" c_"<<fgkFlow<<"{2} = "<<fCumulant[0]<<endl;
150	cout<<" c_"<<fgkFlow<<"{4} = "<<fCumulant[1]<<endl;
151	cout<<" c_"<<fgkFlow<<"{6} = "<<fCumulant[2]<<endl;
152	cout<<" c_"<<fgkFlow<<"{8} = "<<fCumulant[3]<<endl;
153	cout<<"c_"<<fgkFlow<<"{10} = "<<fCumulant[4]<<endl;
154	cout<<"c_"<<fgkFlow<<"{12} = "<<fCumulant[5]<<endl;
155	cout<<"c_"<<fgkFlow<<"{14} = "<<fCumulant[6]<<endl;
156	cout<<"c_"<<fgkFlow<<"{16} = "<<fCumulant[7]<<endl;
157
158	cout<<""<<endl;
159	cout<<"integrated flow: "<<endl;
160
161
162	Double_t fV2=0.,fV4=0.,fV6=0.,fV8=0.;
163	Double_t fSdQ[4]={0.};
164	Double_t fChiQ[4]={0.};
165	if (fCumulant[0]>=0.){
166	fV2=sqrt(fCumulant[0]);
167	//if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV2*fAvM,2.)>0.){
168	//fChiQ[0]=fAvMfV2/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV2fAvM,2.),0.5);
169	//fSdQ[0]=pow(((1./(2.fAvMnEvents))((1.+1.pow(fChiQ[0],2))/(1.*pow(fChiQ[0],2)))),0.5);
170	cout<<" v_"<<fgkFlow<<"{2} = "<<100.fV2<<"%, chi{2} = "<<fChiQ[0]<<", sd{2} = "<<100.fSdQ[0]<<"%"<<endl;//to be improved (2->fgkFlow)
171	//fCommonHistsRes2->FillIntegratedFlow(100.fV2,100.fSdQ[0]);
172	//fCommonHistsRes2->FillChi(fChiQ[0]);
173	fifr->SetBinContent(1,100.*fV2);
174	//}
175	//} else {
176	//cout<<" v_"<<fgkFlow<<"{2} = Im"<<endl;
177	} else {
178	//cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl;
179	}
180	if (fCumulant[1]<=0.){
181	fV4=pow(-fCumulant[1],(1./4.));
182	//if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV4*fAvM,2.)>0.){
183	//fChiQ[1]=fAvMfV4/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV4fAvM,2.),0.5);
184	//fSdQ[1]=(1./(pow(2.fAvMnEvents,.5)))pow((1.+2.pow(fChiQ[1],2)+(1./4.)pow(fChiQ[1],4.)+(1./4.)pow(fChiQ[1],6.))/((1./4.)*pow(fChiQ[1],6.)),.5);
185	cout<<" v_"<<fgkFlow<<"{4} = "<<100.fV4<<"%, chi{4} = "<<fChiQ[1]<<", sd{4} = "<<100.fSdQ[1]<<"%"<<endl;//to be improved (2->fgkFlow)
186	//fCommonHistsRes4->FillInteg8ratedFlow(100.fV4,100.fSdQ[1]);
187	//fCommonHistsRes4->FillChi(fChiQ[1]);
188	fifr->SetBinContent(2,100.*fV4);
189	//} else {
190	//cout<<" v_"<<fgkFlow<<"{4} = Im"<<endl;
191	//}
192	// } else {
193	// cout<<" v_"<<fgkFlow<<"{4} = Im"<<endl;
194	}
195	if (fCumulant[2]>=0.){
196	fV6=pow((1./4.)*fCumulant[2],(1./6.));
197	//if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV6*fAvM,2.)>0.){
198	//fChiQ[2]=fAvMfV6/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV6fAvM,2.),0.5);
199	//fSdQ[2]=(1./(pow(2.fAvMnEvents,.5)))pow((3.+18.pow(fChiQ[2],2)+9.pow(fChiQ[2],4.)+28.pow(fChiQ[2],6.)+12.pow(fChiQ[2],8.)+24.pow(fChiQ[2],10.))/(24.*pow(fChiQ[2],10.)),.5);
200	cout<<" v_"<<fgkFlow<<"{6} = "<<100.fV6<<"%, chi{6} = "<<fChiQ[2]<<", sd{6} = "<<100.fSdQ[2]<<"%"<<endl;//to be improved (2->fgkFlow)
201	//fCommonHistsRes6->FillIntegratedFlow(100.fV6,100.fSdQ[2]);
202	//fCommonHistsRes6->FillChi(fChiQ[2]);
203	fifr->SetBinContent(3,100.*fV6);
204	//} else {
205	// cout<<" v_"<<fgkFlow<<"{6} = Im"<<endl;
206	// }
207	} else {
208	//cout<<" v_"<<fgkFlow<<"{6} = Im"<<endl;
209	}
210	if (fCumulant[3]<=0.){
211	fV8=pow(-(1./33.)*fCumulant[3],(1./8.));
212	//if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV8*fAvM,2.)>0.){
213	//fChiQ[3]=fAvMfV8/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV8fAvM,2.),0.5);
214	//fSdQ[3]=(1./(pow(2.fAvMnEvents,.5)))pow((12.+96.pow(fChiQ[3],2)+72.pow(fChiQ[3],4.)+304.pow(fChiQ[3],6.)+257.pow(fChiQ[3],8.)+804.pow(fChiQ[3],10.)+363.pow(fChiQ[3],12.)+726.pow(fChiQ[3],14.))/(726.*pow(fChiQ[3],14.)),.5);
215	cout<<" v_"<<fgkFlow<<"{8} = "<<100.fV8<<"%, chi{8} = "<<fChiQ[3]<<", sd{8} = "<<100.fSdQ[3]<<"%"<<endl;//to be improved (2->fgkFlow)
216	//fCommonHistsRes8->FillIntegratedFlow(100.fV8,100.fSdQ[3]);
217	//fCommonHistsRes8->FillChi(fChiQ[3]);
218	fifr->SetBinContent(4,100.*fV8);
219	//} else {
220	//cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl;a
221	//}
222	} else {
223	//cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl;
224	}
225	if (fCumulant[4]>=0.){//fHistProIntFlow
226	cout<<"v_"<<fgkFlow<<"{10} = "<<100.pow((1./456.)fCumulant[4],(1./10.))<<"%"<<endl;//to be improved (2->fgkFlow)
227	fifr->SetBinContent(5,100.pow((1./456.)fCumulant[4],(1./10.)));
228	} else {
229	cout<<"v_"<<fgkFlow<<"{10} = Im"<<endl; //to be improved (2->fgkFlow)
230	}
231	if (fCumulant[5]<=0.){
232	cout<<"v_"<<fgkFlow<<"{12} = "<<100.pow(-(1./9460.)fCumulant[5],(1./12.))<<"%"<<endl;//to be improved (2->fgkFlow)
233	fifr->SetBinContent(6,100.pow(-(1./9460.)fCumulant[5],(1./12.)));
234	} else {
235	cout<<"v_"<<fgkFlow<<"{12} = Im"<<endl; //to be improved (2->fgkFlow)
236	}
237	if (fCumulant[6]>=0.){
238	cout<<"v_"<<fgkFlow<<"{14} = "<<100.pow((1./274800.)fCumulant[6],(1./14.))<<"%"<<endl;//to be improved (2->fgkFlow)
239	fifr->SetBinContent(7,100.pow((1./274800.)fCumulant[6],(1./14.)));
240	} else {
241	cout<<"v_"<<fgkFlow<<"{14} = Im"<<endl; //to be improved (2->fgkFlow)
242	}
243	if (fCumulant[7]<=0.){
244	cout<<"v_"<<fgkFlow<<"{16} = "<<100.pow(-(1./10643745.)fCumulant[7],(1./16.))<<"%"<<endl;//to be improved (2->fgkFlow)
245	fifr->SetBinContent(8,100.pow(-(1./10643745.)fCumulant[7],(1./16.)));
246	} else {
247	cout<<"v_"<<fgkFlow<<"{16} = Im"<<endl; //to be improved (2->fgkFlow)
248	}
249	//cout<<"***************************"<<endl;
250
251	//DIFFERENTIAL FLOW CALCULATIONS STARTS HERE!!!
252
253	Double_t fX[fgknBins][fgkPmax][fgkQmax]={0.};//see the text bellow relation (11) in PG
254	Double_t fY[fgknBins][fgkPmax][fgkQmax]={0.};
255
256	for(Int_t b=0;b<fgknBins;b++){
257	for(Int_t p=0;p<fgkPmax;p++){
258	for(Int_t q=0;q<fgkQmax;q++){
259	fX[b][p][q]=fDGRe->GetBinContent(b+1,p+1,q+1)/fBvG[p][q];
260	fY[b][p][q]=fDGIm->GetBinContent(b+1,p+1,q+1)/fBvG[p][q];
261	}
262	}
263	}
264
265
266	Double_t fD[fgknBins][fgkPmax]={0.};//implementing relation (11) from PG
267
268	for (Int_t b=0;b<fgknBins;b++){
269	for (Int_t p=0;p<fgkPmax;p++){
270	Double_t fTempHere3=0.;
271	for (Int_t q=0;q<fgkQmax;q++){
272	fTempHere3+=cos(fgkMltpl2.qTMath::Pi()/fgkQmax)fX[b][p][q] + sin(fgkMltpl2.qTMath::Pi()/fgkQmax)fY[b][p][q];
273	}
274	fD[b][p]=1.(pow(fR0pow(p+1.,.5),fgkMltpl)/fgkQmax)*fTempHere3;
275	if(fD[b][p]){
276	//cout<<"And this "<<b<<" "<<fD[b][p]<<endl;
277	}
278	}
279	}
280
281	Double_t fDiffCumulant2[fgknBins]={0.};//implementing relation (12) from PG
282	Double_t fDiffCumulant4[fgknBins]={0.};
283	Double_t fDiffCumulant6[fgknBins]={0.};
284	Double_t fDiffCumulant8[fgknBins]={0.};
285
286	for (Int_t b=0;b<fgknBins;b++){
287	fDiffCumulant2[b]=(1./(fR0fR0))(4.fD[b][0]-3.fD[b][1]+(4./3.)fD[b][2]-(1./4.)fD[b][3]);
288	fDiffCumulant4[b]=(1./pow(fR0,4.))((-26./3.)fD[b][0]+(19./2.)fD[b][1]-(14./3.)fD[b][2]+(11./12.)*fD[b][3]);
289	fDiffCumulant6[b]=(1./pow(fR0,6.))(18.fD[b][0]-24.fD[b][1]+14.fD[b][2]-3.*fD[b][3]);
290	fDiffCumulant8[b]=(1./pow(fR0,8.))(-24.fD[b][0]+36.fD[b][1]-24.fD[b][2]+6.*fD[b][3]);
291	}
292
293	Double_t fv2[fgknBins],fv4[fgknBins],fv6[fgknBins],fv8[fgknBins];
294	//Double_t fAvPt[fgknBins];
295	Double_t fSddiff2[fgknBins],fSddiff4[fgknBins];
296
297	cout<<"number of pt bins: "<<fgknBins<<endl;
298	cout<<"****************************************"<<endl;
299	for (Int_t b=0;b<fgknBins;b++){
300	//if(fBinNoOfParticles[b]==0)continue;
301	//fAvPt[b]=fBinMeanPt[b]/fBinNoOfParticles[b];
302	cout<<"pt bin: "<<bfBinWidth<<"-"<<(b+1)fBinWidth<<" GeV"<<endl;
303	//cout<<"number of particles in this pt bin: "<<tempNo->GetBinContent(b+1,1,1)<<endl;
304	//cout<<"mean pt in this bin: "<<fAvPt[b]<<" GeV"<<endl;
305	if(fCumulant[0]>=0){
306	fv2[b]=100.*fDiffCumulant2[b]/pow(fCumulant[0],.5);
307	//if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV2*fAvM,2.)>0.){
308	//fSddiff2[b]=pow((1./(2.fBinNoOfParticles[b]))((1.+pow(fChiQ[0],2.))/pow(fChiQ[0],2.)),0.5);
309	cout<<"v'_2/2{2} = "<<fv2[b]<<"%, "<<" "<<"sd{2} = "<<100.*fSddiff2[b]<<"%"<<endl;
310	//fDiffFlowResults2->SetBinContent(b+1,fv2[b],100.*fSddiff2[b]);
311	fdfr2->SetBinContent(b+1,fv2[b]);
312	//} else {
313	//cout<<"v'_2/2{2} = Im"<<endl;
314	//}
315	}else{
316	cout<<"v'_2/2{2} = Im"<<endl;
317	}
318	if(fCumulant[1]<=0){
319	fv4[b]=-100.*fDiffCumulant4[b]/pow(-fCumulant[1],.75);
320	//if (fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV4*fAvM,2.)>0.){
321	//fSddiff4[b]=pow((1./(2.fBinNoOfParticles[b]))((2.+6.*pow(fChiQ[1],2.)+pow(fChiQ[1],4.)+pow(fChiQ[1],6.))/pow(fChiQ[1],6.)),0.5);
322	cout<<"v'_2/2{4} = "<<fv4[b]<<"%, "<<" "<<"sd{4} = "<<100.*fSddiff4[b]<<"%"<<endl;
323	//fCommonHistsRes4->FillDifferentialFlow(b+1,fv4[b],100.*fSddiff4[b]);
324	fdfr4->SetBinContent(b+1,fv4[b]);
325	//} else {
326	// cout<<"v'_2/2{4} = Im"<<endl;
327	//}
328	}else{
329	cout<<"v'_2/2{4} = Im"<<endl;
330	}
331	if(fCumulant[2]>=0){
332	cout<<"v'_2/2{6} = "<<100.fDiffCumulant6[b]/(4.pow((1./4.)*fCumulant[2],(5./6.)))<<"%"<<endl;
333	fv6[b]=100.fDiffCumulant6[b]/(4.pow((1./4.)*fCumulant[2],(5./6.)));
334	//fCommonHistsRes6->FillDifferentialFlow(b+1,fv6[b],0.);
335	fdfr6->SetBinContent(b+1,fv6[b]);
336	}else{
337	cout<<"v'_2/2{6} = Im"<<endl;
338	}
339	if(fCumulant[3]<=0){
340	cout<<"v'_2/2{8} = "<<-100.fDiffCumulant8[b]/(33.pow(-(1./33.)*fCumulant[3],(7./8.)))<<"%"<<endl;
341	fv8[b]=-100.fDiffCumulant8[b]/(33.pow(-(1./33.)*fCumulant[3],(7./8.)));
342	//fCommonHistsRes8->FillDifferentialFlow(b+1,fv8[b],0.);
343	fdfr8->SetBinContent(b+1,fv8[b]);
344	}else{
345	cout<<"v'_2/2{8} = Im"<<endl;
346	}
347	cout<<"****************************************"<<endl;
348	}
349	}
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364