From d613093881abc977b7cf8e18a9475826197a3705 Mon Sep 17 00:00:00 2001 From: snelling Date: Thu, 10 Jun 2010 20:17:50 +0000 Subject: [PATCH] cleanup of GFC cumulants (some well deserved attention) --- .../AliFlowCommon/AliCumulantsFunctions.cxx | 2152 ---------- .../AliFlowCommon/AliCumulantsFunctions.h | 112 - .../AliFlowAnalysisWithCumulants.cxx | 3658 ++++++++++++----- .../AliFlowAnalysisWithCumulants.h | 462 ++- .../AliFlowCommon/AliFlowCumuConstants.cxx | 61 - .../FLOW/AliFlowCommon/AliFlowCumuConstants.h | 111 - PWG2/PWG2flowCommonLinkDef.h | 2 - PWG2/libPWG2flowCommon.pkg | 2 - 8 files changed, 2904 insertions(+), 3656 deletions(-) delete mode 100644 PWG2/FLOW/AliFlowCommon/AliCumulantsFunctions.cxx delete mode 100644 PWG2/FLOW/AliFlowCommon/AliCumulantsFunctions.h delete mode 100644 PWG2/FLOW/AliFlowCommon/AliFlowCumuConstants.cxx delete mode 100644 PWG2/FLOW/AliFlowCommon/AliFlowCumuConstants.h diff --git a/PWG2/FLOW/AliFlowCommon/AliCumulantsFunctions.cxx b/PWG2/FLOW/AliFlowCommon/AliCumulantsFunctions.cxx deleted file mode 100644 index d68be0cc214..00000000000 --- a/PWG2/FLOW/AliFlowCommon/AliCumulantsFunctions.cxx +++ /dev/null @@ -1,2152 +0,0 @@ -/************************************************************************* -* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. * -* * -* Author: The ALICE Off-line Project. * -* Contributors are mentioned in the code where appropriate. * -* * -* Permission to use, copy, modify and distribute this software and its * -* documentation strictly for non-commercial purposes is hereby granted * -* without fee, provided that the above copyright notice appears in all * -* copies and that both the copyright notice and this permission notice * -* appear in the supporting documentation. The authors make no claims * -* about the suitability of this software for any purpose. It is * -* provided "as is" without express or implied warranty. * -**************************************************************************/ - -/********************************** - * functions and equations needed * - * for calculation of cumulants * - * and final flow estimates * - * * - * author: Ante Bilandzic * - * (anteb@nikhef.nl) * - *********************************/ - -#define AliCumulantsFunctions_cxx - -#include "Riostream.h" -#include "TChain.h" -#include "TFile.h" -#include "TList.h" -#include "TParticle.h" -#include -#include - -#include "TProfile.h" -#include "TProfile2D.h" -#include "TProfile3D.h" -#include "TH1.h" -#include "TH1D.h" - -#include "AliFlowEventSimple.h" -#include "AliFlowTrackSimple.h" -#include "AliFlowAnalysisWithCumulants.h" -#include "AliFlowCumuConstants.h" -#include "AliFlowCommonConstants.h" -#include "AliFlowCommonHist.h" -#include "AliFlowCommonHistResults.h" -#include "AliCumulantsFunctions.h" - -ClassImp(AliCumulantsFunctions) - -//================================================================================================================_ - -AliCumulantsFunctions::AliCumulantsFunctions(): - fIntGenFun(NULL), - fIntGenFun4(NULL),//only for other system of Eq. - fIntGenFun6(NULL),//only for other system of Eq. - fIntGenFun8(NULL),//only for other system of Eq. - fIntGenFun16(NULL),//only for other system of Eq. - fAvMult4(NULL),//only for other system of Eq. - fAvMult6(NULL),//only for other system of Eq. - fAvMult8(NULL),//only for other system of Eq. - fAvMult16(NULL),//only for other system of Eq. - fDiffPtRPGenFunRe(NULL), - fDiffPtRPGenFunIm(NULL), - fPtBinRPNoOfParticles(NULL), - fDiffEtaRPGenFunRe(NULL), - fDiffEtaRPGenFunIm(NULL), - fEtaBinRPNoOfParticles(NULL), - fDiffPtPOIGenFunRe(NULL), - fDiffPtPOIGenFunIm(NULL), - fPtBinPOINoOfParticles(NULL), - fDiffEtaPOIGenFunRe(NULL), - fDiffEtaPOIGenFunIm(NULL), - fEtaBinPOINoOfParticles(NULL), - fifr(NULL), - fdfr2(NULL), - fdfr4(NULL), - fdfr6(NULL), - fdfr8(NULL), - fAvMult(NULL), - fQVector(NULL), - fAverageOfSquaredWeight(NULL), - fchr2nd(NULL), - fchr4th(NULL), - fchr6th(NULL), - fchr8th(NULL), - fch(NULL)//common control histograms - /* - fdRe0(NULL), - fdRe1(NULL), - fdRe2(NULL), - fdRe3(NULL), - fdRe4(NULL), - fdRe5(NULL), - fdRe6(NULL), - fdRe7(NULL), - fdIm0(NULL), - fdIm1(NULL), - fdIm2(NULL), - fdIm3(NULL), - fdIm4(NULL), - fdIm5(NULL), - fdIm6(NULL), - fdIm7(NULL) - */ -{ - //default constructor -} - -AliCumulantsFunctions::~AliCumulantsFunctions() -{ - //destructor -} - -AliCumulantsFunctions::AliCumulantsFunctions(TProfile2D *intGenFun, TProfile2D *intGenFun4, TProfile2D *intGenFun6, TProfile2D *intGenFun8, TProfile2D *intGenFun16, TProfile *avMult4, TProfile *avMult6, TProfile *avMult8, TProfile *avMult16, TProfile3D *diffPtRPGenFunRe, TProfile3D *diffPtRPGenFunIm, TProfile *ptBinRPNoOfParticles, TProfile3D *diffEtaRPGenFunRe, TProfile3D *diffEtaRPGenFunIm, TProfile *etaBinRPNoOfParticles, TProfile3D *diffPtPOIGenFunRe, TProfile3D *diffPtPOIGenFunIm, TProfile *ptBinPOINoOfParticles, TProfile3D *diffEtaPOIGenFunRe, TProfile3D *diffEtaPOIGenFunIm, TProfile *etaBinPOINoOfParticles, TH1D *ifr, TH1D *dfr2, TH1D *dfr4, TH1D *dfr6, TH1D *dfr8, TProfile *avMult, TProfile *qVector, TProfile *averageOfSquaredWeight, AliFlowCommonHistResults *chr2nd, AliFlowCommonHistResults *chr4th, AliFlowCommonHistResults *chr6th, AliFlowCommonHistResults *chr8th, AliFlowCommonHist *ch): - fIntGenFun(intGenFun), - fIntGenFun4(intGenFun4),//only for other system of Eq. - fIntGenFun6(intGenFun6),//only for other system of Eq. - fIntGenFun8(intGenFun8),//only for other system of Eq. - fIntGenFun16(intGenFun16),//only for other system of Eq. - fAvMult4(avMult4),//only for other system of Eq. - fAvMult6(avMult6),//only for other system of Eq. - fAvMult8(avMult8),//only for other system of Eq. - fAvMult16(avMult16),//only for other system of Eq. - fDiffPtRPGenFunRe(diffPtRPGenFunRe), - fDiffPtRPGenFunIm(diffPtRPGenFunIm), - fPtBinRPNoOfParticles(ptBinRPNoOfParticles), - fDiffEtaRPGenFunRe(diffEtaRPGenFunRe), - fDiffEtaRPGenFunIm(diffEtaRPGenFunIm), - fEtaBinRPNoOfParticles(etaBinRPNoOfParticles), - fDiffPtPOIGenFunRe(diffPtPOIGenFunRe), - fDiffPtPOIGenFunIm(diffPtPOIGenFunIm), - fPtBinPOINoOfParticles(ptBinPOINoOfParticles), - fDiffEtaPOIGenFunRe(diffEtaPOIGenFunRe), - fDiffEtaPOIGenFunIm(diffEtaPOIGenFunIm), - fEtaBinPOINoOfParticles(etaBinPOINoOfParticles), - fifr(ifr), - fdfr2(dfr2), - fdfr4(dfr4), - fdfr6(dfr6), - fdfr8(dfr8), - fAvMult(avMult), - fQVector(qVector), - fAverageOfSquaredWeight(averageOfSquaredWeight), - fchr2nd(chr2nd), - fchr4th(chr4th), - fchr6th(chr6th), - fchr8th(chr8th), - fch(ch)//common control histograms - /* - fdRe0(dRe0), - fdRe1(dRe1), - fdRe2(dRe2), - fdRe3(dRe3), - fdRe4(dRe4), - fdRe5(dRe5), - fdRe6(dRe6), - fdRe7(dRe7), - fdIm0(dIm0), - fdIm1(dIm1), - fdIm2(dIm2), - fdIm3(dIm3), - fdIm4(dIm4), - fdIm5(dIm5), - fdIm6(dIm6), - fdIm7(dIm7) - */ -{ - //custom constructor -} - -//================================================================================================================ - -void AliCumulantsFunctions::Calculate() -{ - //calculate cumulants and final integrated and differential flow estimates and store the results into output histograms - const Int_t cQmax=AliFlowCumuConstants::GetMaster()->GetQmax(); //needed for numerics - const Int_t cPmax=AliFlowCumuConstants::GetMaster()->GetPmax(); //needed for numerics - const Int_t cQmax4=AliFlowCumuConstants::GetMaster()->GetQmax4(); //needed for numerics - const Int_t cPmax4=AliFlowCumuConstants::GetMaster()->GetPmax4(); //needed for numerics - const Int_t cQmax6=AliFlowCumuConstants::GetMaster()->GetQmax6(); //needed for numerics - const Int_t cPmax6=AliFlowCumuConstants::GetMaster()->GetPmax6(); //needed for numerics - const Int_t cQmax8=AliFlowCumuConstants::GetMaster()->GetQmax8(); //needed for numerics - const Int_t cPmax8=AliFlowCumuConstants::GetMaster()->GetPmax8(); //needed for numerics - const Int_t cQmax16=AliFlowCumuConstants::GetMaster()->GetQmax16(); //needed for numerics - const Int_t cPmax16=AliFlowCumuConstants::GetMaster()->GetPmax16(); //needed for numerics - - const Int_t cFlow=AliFlowCumuConstants::GetMaster()->GetFlow(); //integrated flow coefficient to be calculated - const Int_t cMltpl=AliFlowCumuConstants::GetMaster()->GetMltpl(); //the multiple in p=m*n (diff. flow) - const Int_t cnBinsPt=100; //number of pt bins //to be improved - const Int_t cnBinsEta=80; //number of eta bins //to be improved - - Double_t fR0=AliFlowCumuConstants::GetMaster()->GetR0(); //needed for numerics - //Double_t fPtMax=AliFlowCommonConstants::GetMaster()->GetPtMax(); //maximum pt - //Double_t fPtMin=AliFlowCommonConstants::GetMaster()->GetPtMin(); //minimum pt - //Double_t fBinWidthPt=(fPtMax-fPtMin)/cnBinsPt; //width of pt bin (in GeV) - - Bool_t fOtherEquations=AliFlowCumuConstants::GetMaster()->GetOtherEquations(); //numerical equations for cumulants solved up to different highest order - - //avarage selected multiplicity - Double_t dAvM=0.; - if(fAvMult) - { - dAvM=fAvMult->GetBinContent(1); - } - - //number of events - Int_t nEvents=0; - if(fAvMult) - { - nEvents=(Int_t)(fAvMult->GetBinEntries(1)); - } - - // - Double_t dAvQx=0.,dAvQy=0.,dAvQ2x=0.,dAvQ2y=0.; - if(fQVector) - { - dAvQx = fQVector->GetBinContent(1); // - dAvQy = fQVector->GetBinContent(2); // - dAvQ2x = fQVector->GetBinContent(3); //<(Q_x)^2> - dAvQ2y = fQVector->GetBinContent(4); //<(Q_y)^2> - } - - // - Double_t dAvw2 = 1.; - if(nEvents>0) - { - dAvw2 = fAverageOfSquaredWeight->GetBinContent(1); - if(dAvw2 == 0) - { - cout< is empty. Nothing will be calculated !!!!"< - TMatrixD dAvG(cPmax, cQmax); - Bool_t someAvGEntryIsNegative=kFALSE; - if(fIntGenFun) - { - for(Int_t p=0;pGetBinContent(p+1,q+1); - if(dAvG(p,q)<0.) - { - someAvGEntryIsNegative=kTRUE; - } - } - } - } - - ///////////////////////////////////////////////////////////////////////////// - //////////////////gen. function for the cumulants//////////////////////////// - ///////////////////////////////////////////////////////////////////////////// - - TMatrixD dC(cPmax, cQmax);//C[p][q] - if(dAvM>0 && someAvGEntryIsNegative==kFALSE) - { - for(Int_t p=0;p - Double_t tempHere=0.; - for(Int_t p=0;p0 && dAvM>0 && dAvw2>0 && cumulant[0]>=0.) - { - if((dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(cumulant[0],(1./2.))*(dAvM/dAvw2),2.)>0.)) - { - chiQ[0]=(dAvM*dV2)/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2*dAvM/dAvw2,2.),0.5); // to be improved, analogously for higher orders - } - if(chiQ[0]) - { - sdQ[0]=pow(((1./(2.*dAvM*nEvents))*((1.+2.*pow(chiQ[0],2))/(2.*pow(chiQ[0],2)))),0.5); - } - cout<<" v_"<0 && dAvM>0 && dAvw2>0 && cumulant[1]<=0.) - { - if((dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(-cumulant[1],(1./4.))*(dAvM/dAvw2),2.)>0.)) - { - chiQ[1]=(dAvM*dV4)/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4*dAvM/dAvw2,2.),0.5); - } - if(chiQ[1]) - { - sdQ[1]=(1./(pow(2.*dAvM*nEvents,0.5)))*pow((1.+4.*pow(chiQ[1],2)+1.*pow(chiQ[1],4.)+2.*pow(chiQ[1],6.))/(2.*pow(chiQ[1],6.)),0.5); - } - cout<<" v_"<0 && dAvM>0 && dAvw2>0 && cumulant[2]>=0.) - { - if((dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow((1./4.)*cumulant[2],(1./6.))*(dAvM/dAvw2),2.)>0.)) - { - chiQ[2]=(dAvM*dV6)/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV6*dAvM/dAvw2,2.),0.5); - } - if(chiQ[2]) - { - sdQ[2]=(1./(pow(2.*dAvM*nEvents,0.5)))*pow((3.+18.*pow(chiQ[2],2)+9.*pow(chiQ[2],4.)+28.*pow(chiQ[2],6.)+12.*pow(chiQ[2],8.)+24.*pow(chiQ[2],10.))/(24.*pow(chiQ[2],10.)),0.5); - } - cout<<" v_"<0 && dAvM>0 && dAvw2>0 && cumulant[3]<=0.) - { - if((dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(-(1./33.)*cumulant[3],(1./8.))*(dAvM/dAvw2),2.)>0.)) - { - chiQ[3]=(dAvM*dV8)/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV8*dAvM/dAvw2,2.),0.5); - } - if(chiQ[3]) - { - sdQ[3]=(1./(pow(2.*dAvM*nEvents,0.5)))*pow((12.+96.*pow(chiQ[3],2.)+72.*pow(chiQ[3],4.)+304.*pow(chiQ[3],6.)+257.*pow(chiQ[3],8.)+804.*pow(chiQ[3],10.)+363.*pow(chiQ[3],12.)+726.*pow(chiQ[3],14.))/(726.*pow(chiQ[3],14.)),0.5); - } - cout<<" v_"<=0.) - { - cout<<"v_"<GetBinEntries(b+1); - - for(Int_t p=0;pGetBinContent(b+1,p+1,q+1)/dAvG(p,q); - } - if(fDiffPtRPGenFunIm) - { - ptYRP[index3d(b,p,q,cnBinsPt,cPmax)]=fDiffPtRPGenFunIm->GetBinContent(b+1,p+1,q+1)/dAvG(p,q); - } - } - } - } - } - - TVectorD etaXRP(cnBinsEta*cPmax*cQmax); - TVectorD etaYRP(cnBinsEta*cPmax*cQmax); - TVectorD etaBinRPNoOfParticles(cnBinsEta); - - //3D profiles (for eta) - for(Int_t b=0;bGetBinEntries(b); - for(Int_t p=0;pGetBinContent(b+1,p+1,q+1)/dAvG(p,q); - } - if(fDiffEtaRPGenFunIm) - { - etaYRP[index3d(b,p,q,cnBinsEta,cPmax)]=fDiffEtaRPGenFunIm->GetBinContent(b+1,p+1,q+1)/dAvG(p,q); - } - } - } - } - } - - - //------------------------------------------------------------------------------------------------------------------------------- - //final results for differential flow (in pt): - TMatrixD ptDRP(cnBinsPt, cPmax); - Double_t tempSumForptDRP=0.; - for (Int_t b=0;b0) - { - v2ptRP[b]=ptRPDiffCumulant2[b]/pow(cumulant[0],0.5); - if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2*dAvM,2.)>0. && ptBinRPNoOfParticles[b]>0.) - { - if(chiQ[0]>0) - { - sdRPDiff2pt[b]=pow((1./(2.*ptBinRPNoOfParticles[b]))*((1.+pow(chiQ[0],2.))/pow(chiQ[0],2.)),0.5); - } - //cout<<"v'_2/2{2} = "<SetBinContent(b+1,v2ptRP[b]); - fdfr2->SetBinError(b+1,sdRPDiff2pt[b]); - //common histogram (to be removed): - fchr2nd->FillDifferentialFlow(b+1,v2ptRP[b],sdRPDiff2pt[b]); - // Fill common result histogram: - if(TMath::Abs(v2ptRP[b])>1.e-44) fchr2nd->FillDifferentialFlowPtRP(b+1,v2ptRP[b],sdRPDiff2pt[b]); - - } else { - //cout<<"v'_2/2{2} = Im"<0.&&ptBinRPNoOfParticles[b]>0.) // to be improved - { - if(chiQ[1]>0) - { - sdRPDiff4pt[b]=pow((1./(2.*ptBinRPNoOfParticles[b]))*((2.+6.*pow(chiQ[1],2.)+pow(chiQ[1],4.)+pow(chiQ[1],6.))/pow(chiQ[1],6.)),0.5); - } - //cout<<"v'_2/2{4} = "<SetBinContent(b+1,v4ptRP[b]); - fdfr4->SetBinError(b+1,sdRPDiff4pt[b]); - //common histogram (to be removed): - fchr4th->FillDifferentialFlow(b+1,v4ptRP[b],sdRPDiff4pt[b]); - // Fill common result histogram: - if(TMath::Abs(v4ptRP[b])>1.e-44) fchr4th->FillDifferentialFlowPtRP(b+1,v4ptRP[b],sdRPDiff4pt[b]); - - } else { - //cout<<"v'_2/2{4} = Im"<0){ - //cout<<"v'_2/2{6} = "<<100.*ptRPDiffCumulant6[b]/(4.*pow((1./4.)*cumulant[2],(5./6.)))<<"%"<SetBinContent(b+1,v6ptRP[b]); - //common histogram (to be removed): - fchr6th->FillDifferentialFlow(b+1,v6ptRP[b],0.); - // Fill common result histogram: - if(TMath::Abs(v6ptRP[b])>1.e-44) fchr6th->FillDifferentialFlowPtRP(b+1,v6ptRP[b],0.); - - }else{ - //cout<<"v'_2/2{6} = Im"<SetBinContent(b+1,v8ptRP[b]); - //common histogram (to be removed): - fchr8th->FillDifferentialFlow(b+1,v8ptRP[b],0.); - // Fill common result histogram: - if(TMath::Abs(v8ptRP[b])>1.e-44) fchr8th->FillDifferentialFlowPtRP(b+1,v8ptRP[b],0.); - - }else{ - //cout<<"v'_2/2{8} = Im"<GetHistPtRP()) - { - dSumOfYieldRP+=(fch->GetHistPtRP())->GetBinContent(b); - if(fchr2nd->GetHistDiffFlowPtRP()) - { - dV2RP+=((fchr2nd->GetHistDiffFlowPtRP())->GetBinContent(b))*(fch->GetHistPtRP())->GetBinContent(b); - dV2RPError+=pow((fch->GetHistPtRP())->GetBinContent(b),2.)*pow((fchr2nd->GetHistDiffFlowPtRP())->GetBinError(b),2.); - } - if(fchr4th->GetHistDiffFlowPtRP()) - { - dV4RP+=((fchr4th->GetHistDiffFlowPtRP())->GetBinContent(b))*(fch->GetHistPtRP())->GetBinContent(b); - dV4RPError+=pow((fch->GetHistPtRP())->GetBinContent(b),2.)*pow((fchr4th->GetHistDiffFlowPtRP())->GetBinError(b),2.); - } - if(fchr6th->GetHistDiffFlowPtRP()) - { - dV6RP+=((fchr6th->GetHistDiffFlowPtRP())->GetBinContent(b))*(fch->GetHistPtRP())->GetBinContent(b); - dV6RPError+=pow((fch->GetHistPtRP())->GetBinContent(b),2.)*pow((fchr6th->GetHistDiffFlowPtRP())->GetBinError(b),2.); - } - if(fchr8th->GetHistDiffFlowPtRP()) - { - dV8RP+=((fchr8th->GetHistDiffFlowPtRP())->GetBinContent(b))*(fch->GetHistPtRP())->GetBinContent(b); - dV8RPError+=pow((fch->GetHistPtRP())->GetBinContent(b),2.)*pow((fchr8th->GetHistDiffFlowPtRP())->GetBinError(b),2.); - } - } - } - - if(dSumOfYieldRP) - { - dV2RP/=dSumOfYieldRP; - dV2RPError/=(dSumOfYieldRP*dSumOfYieldRP); - dV4RP/=dSumOfYieldRP; - dV4RPError/=(dSumOfYieldRP*dSumOfYieldRP); - dV6RP/=dSumOfYieldRP; - dV6RPError/=(dSumOfYieldRP*dSumOfYieldRP); - dV8RP/=dSumOfYieldRP; - dV8RPError/=(dSumOfYieldRP*dSumOfYieldRP); - fchr2nd->FillIntegratedFlowRP(dV2RP,pow(dV2RPError,0.5)); - fchr4th->FillIntegratedFlowRP(dV4RP,pow(dV4RPError,0.5)); - fchr6th->FillIntegratedFlowRP(dV6RP,pow(dV6RPError,0.5));//to be improved (errors needed) - fchr8th->FillIntegratedFlowRP(dV8RP,pow(dV8RPError,0.5));//to be improved (errors needed) - } - - cout<0) - { - v2etaRP[b]=etaRPDiffCumulant2[b]/pow(cumulant[0],0.5); - if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2*dAvM,2.)>0. && etaBinRPNoOfParticles[b]>0.) // to be improved - { - if(chiQ[0]>0) - { - sdRPDiff2eta[b]=pow((1./(2.*etaBinRPNoOfParticles[b]))*((1.+pow(chiQ[0],2.))/pow(chiQ[0],2.)),0.5); - } - //cout<<"v'_2/2{2} = "<SetBinContent(b+1,v2etaRP[b]); - //fdfr2->SetBinError(b+1,sdDiff2eta[b]); - //common histogram: - //fchr2nd->FillDifferentialFlow(b+1,v2etaRP[b],sdDiff2eta[b]) - // Fill common result histogram: - if(TMath::Abs(v2etaRP[b])>1.e-44) fchr2nd->FillDifferentialFlowEtaRP(b+1,v2etaRP[b],sdRPDiff2eta[b]); - - } else { - //cout<<"v'_2/2{2} = Im"<0.&&etaBinRPNoOfParticles[b]>0.) // to be improved - { - if(chiQ[1]) - { - sdRPDiff4eta[b]=pow((1./(2.*etaBinRPNoOfParticles[b]))*((2.+6.*pow(chiQ[1],2.)+pow(chiQ[1],4.)+pow(chiQ[1],6.))/pow(chiQ[1],6.)),0.5); - } - //cout<<"v'_2/2{4} = "<SetBinContent(b+1,v4eta[b]); - //fdfr4->SetBinError(b+1,sdDiff4eta[b]); - //common histogram: - //fchr4th->FillDifferentialFlow(b+1,v4eta[b],sdDiff4eta[b]); - // Fill common result histogram: - if(TMath::Abs(v4etaRP[b])>1.e-44) fchr4th->FillDifferentialFlowEtaRP(b+1,v4etaRP[b],sdRPDiff4eta[b]); - - } else { - //cout<<"v'_2/2{4} = Im"<0){ - //cout<<"v'_2/2{6} = "<<100.*etaRPDiffCumulant6[b]/(4.*pow((1./4.)*cumulant[2],(5./6.)))<<"%"<SetBinContent(b+1,v6eta[b]); - //common histogram: - //fchr6th->FillDifferentialFlow(b+1,v6eta[b],0.); - // Fill common result histogram: - if(TMath::Abs(v6etaRP[b])>1.e-44) fchr6th->FillDifferentialFlowEtaRP(b+1,v6etaRP[b],0.); - - }else{ - //cout<<"v'_2/2{6} = Im"<SetBinContent(b+1,v8eta[b]); - //common histogram: - //fchr8th->FillDifferentialFlow(b+1,v8eta[b],0.); - // Fill common result histogram: - if(TMath::Abs(v8etaRP[b])>1.e-44) fchr8th->FillDifferentialFlowEtaRP(b+1,v8etaRP[b],0.); - - }else{ - //cout<<"v'_2/2{8} = Im"<GetBinEntries(b+1); - for(Int_t p=0;pGetBinContent(b+1,p+1,q+1)/dAvG(p,q); - } - if(fDiffPtPOIGenFunIm) - { - ptY[index3d(b,p,q,cnBinsPt,cPmax)]=fDiffPtPOIGenFunIm->GetBinContent(b+1,p+1,q+1)/dAvG(p,q); - } - } - } - } - } - - TVectorD etaX(cnBinsEta*cPmax*cQmax); - TVectorD etaY(cnBinsEta*cPmax*cQmax); - TVectorD etaBinPOINoOfParticles(cnBinsEta); - - //3D profiles (for eta) - for(Int_t b=0;bGetBinEntries(b+1); - for(Int_t p=0;pGetBinContent(b+1,p+1,q+1)/dAvG(p,q); - } - if(fDiffEtaPOIGenFunIm) - { - etaY[index3d(b,p,q,cnBinsEta,cPmax)]=fDiffEtaPOIGenFunIm->GetBinContent(b+1,p+1,q+1)/dAvG(p,q); - } - } - } - } - } - - /* - if(dAvM){ - for(Int_t b=0;bGetBinContent(b+1,q+1)/AvG[0][q]; - Y[b][0][q]=fdIm0->GetBinContent(b+1,q+1)/AvG[0][q]; - //-------------------------------------------------- - X[b][1][q]=fdRe1->GetBinContent(b+1,q+1)/AvG[1][q]; - Y[b][1][q]=fdIm1->GetBinContent(b+1,q+1)/AvG[1][q]; - //-------------------------------------------------- - X[b][2][q]=fdRe2->GetBinContent(b+1,q+1)/AvG[2][q]; - Y[b][2][q]=fdIm2->GetBinContent(b+1,q+1)/AvG[2][q]; - //-------------------------------------------------- - X[b][3][q]=fdRe3->GetBinContent(b+1,q+1)/AvG[3][q]; - Y[b][3][q]=fdIm3->GetBinContent(b+1,q+1)/AvG[3][q]; - //-------------------------------------------------- - X[b][4][q]=fdRe4->GetBinContent(b+1,q+1)/AvG[4][q]; - Y[b][4][q]=fdIm4->GetBinContent(b+1,q+1)/AvG[4][q]; - //-------------------------------------------------- - X[b][5][q]=fdRe5->GetBinContent(b+1,q+1)/AvG[5][q]; - Y[b][5][q]=fdIm5->GetBinContent(b+1,q+1)/AvG[5][q]; - //-------------------------------------------------- - X[b][6][q]=fdRe6->GetBinContent(b+1,q+1)/AvG[6][q]; - Y[b][6][q]=fdIm6->GetBinContent(b+1,q+1)/AvG[6][q]; - //-------------------------------------------------- - X[b][7][q]=fdRe7->GetBinContent(b+1,q+1)/AvG[7][q]; - Y[b][7][q]=fdIm7->GetBinContent(b+1,q+1)/AvG[7][q]; - } - //} - } - } - */ - - //------------------------------------------------------------------------------------------------------------------------------- - //final results for differential flow (in pt): - TMatrixD ptD(cnBinsPt,cPmax); - Double_t tempSumForPtD=0.; - for (Int_t b=0;b0) - { - v2pt[b]=ptDiffCumulant2[b]/pow(cumulant[0],0.5); - if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2*dAvM,2.)>0.&&ptBinPOINoOfParticles[b]>0.) - { - if(chiQ[0]>0) - { - sdDiff2pt[b]=pow((1./(2.*ptBinPOINoOfParticles[b]))*((1.+pow(chiQ[0],2.))/pow(chiQ[0],2.)),0.5); - } - //cout<<"v'_2/2{2} = "<SetBinContent(b+1,v2pt[b]); - fdfr2->SetBinError(b+1,sdDiff2pt[b]); - //common histogram: - fchr2nd->FillDifferentialFlow(b+1,v2pt[b],sdDiff2pt[b]); - // Fill common result histogram: - if(TMath::Abs(v2pt[b])>1.e-44) fchr2nd->FillDifferentialFlowPtPOI(b+1,v2pt[b],sdDiff2pt[b]); - - } else { - //cout<<"v'_2/2{2} = Im"<0.&&ptBinPOINoOfParticles[b]>0.) - { - if(chiQ[1]>0) - { - sdDiff4pt[b]=pow((1./(2.*ptBinPOINoOfParticles[b]))*((2.+6.*pow(chiQ[1],2.)+pow(chiQ[1],4.)+pow(chiQ[1],6.))/pow(chiQ[1],6.)),0.5); - } - //cout<<"v'_2/2{4} = "<SetBinContent(b+1,v4pt[b]); - fdfr4->SetBinError(b+1,sdDiff4pt[b]); - //common histogram: - fchr4th->FillDifferentialFlow(b+1,v4pt[b],sdDiff4pt[b]); - // Fill common result histogram: - if(TMath::Abs(v4pt[b])>1.e-44) fchr4th->FillDifferentialFlowPtPOI(b+1,v4pt[b],sdDiff4pt[b]); - - } else { - //cout<<"v'_2/2{4} = Im"<0){ - //cout<<"v'_2/2{6} = "<<100.*ptDiffCumulant6[b]/(4.*pow((1./4.)*cumulant[2],(5./6.)))<<"%"<SetBinContent(b+1,v6pt[b]); - //common histogram: - fchr6th->FillDifferentialFlow(b+1,v6pt[b],0.); - // Fill common result histogram: - if(TMath::Abs(v6pt[b])>1.e-44) fchr6th->FillDifferentialFlowPtPOI(b+1,v6pt[b],0.); - - }else{ - //cout<<"v'_2/2{6} = Im"<SetBinContent(b+1,v8pt[b]); - //common histogram: - fchr8th->FillDifferentialFlow(b+1,v8pt[b],0.); - // Fill common result histogram: - if(TMath::Abs(v8pt[b])>1.e-44) fchr8th->FillDifferentialFlowPtPOI(b+1,v8pt[b],0.); - - }else{ - //cout<<"v'_2/2{8} = Im"<GetHistPtPOI()) - { - dSumOfYieldPOI+=(fch->GetHistPtPOI())->GetBinContent(b); - if(fchr2nd->GetHistDiffFlowPtPOI()) - { - dV2POI+=((fchr2nd->GetHistDiffFlowPtPOI())->GetBinContent(b))*(fch->GetHistPtPOI())->GetBinContent(b); - dV2POIError+=pow((fch->GetHistPtPOI())->GetBinContent(b),2.)*pow((fchr2nd->GetHistDiffFlowPtPOI())->GetBinError(b),2.); - } - if(fchr4th->GetHistDiffFlowPtPOI()) - { - dV4POI+=((fchr4th->GetHistDiffFlowPtPOI())->GetBinContent(b))*(fch->GetHistPtPOI())->GetBinContent(b); - dV4POIError+=pow((fch->GetHistPtPOI())->GetBinContent(b),2.)*pow((fchr4th->GetHistDiffFlowPtPOI())->GetBinError(b),2.); - } - if(fchr6th->GetHistDiffFlowPtPOI()) - { - dV6POI+=((fchr6th->GetHistDiffFlowPtPOI())->GetBinContent(b))*(fch->GetHistPtPOI())->GetBinContent(b); - dV6POIError+=pow((fch->GetHistPtPOI())->GetBinContent(b),2.)*pow((fchr6th->GetHistDiffFlowPtPOI())->GetBinError(b),2.); - } - if(fchr8th->GetHistDiffFlowPtPOI()) - { - dV8POI+=((fchr8th->GetHistDiffFlowPtPOI())->GetBinContent(b))*(fch->GetHistPtPOI())->GetBinContent(b); - dV8POIError+=pow((fch->GetHistPtPOI())->GetBinContent(b),2.)*pow((fchr8th->GetHistDiffFlowPtPOI())->GetBinError(b),2.); - } - } - } - - if(dSumOfYieldPOI) - { - dV2POI/=dSumOfYieldPOI; - dV2POIError/=(dSumOfYieldPOI*dSumOfYieldPOI); - dV4POI/=dSumOfYieldPOI; - dV4POIError/=(dSumOfYieldPOI*dSumOfYieldPOI); - dV6POI/=dSumOfYieldPOI; - dV6POIError/=(dSumOfYieldPOI*dSumOfYieldPOI); - dV8POI/=dSumOfYieldPOI; - dV8POIError/=(dSumOfYieldPOI*dSumOfYieldPOI); - fchr2nd->FillIntegratedFlowPOI(dV2POI,pow(dV2POIError,0.5)); - fchr4th->FillIntegratedFlowPOI(dV4POI,pow(dV4POIError,0.5)); - fchr6th->FillIntegratedFlowPOI(dV6POI,pow(dV6POIError,0.5));//to be improved (errors needed) - fchr8th->FillIntegratedFlowPOI(dV8POI,pow(dV8POIError,0.5));//to be improved (errors needed) - } - - cout<0) - { - v2eta[b]=etaDiffCumulant2[b]/pow(cumulant[0],0.5); - if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2*dAvM,2.)>0.&&etaBinPOINoOfParticles[b]>0.) // to be improved - { - if(chiQ[0]>0) - { - sdDiff2eta[b]=pow((1./(2.*etaBinPOINoOfParticles[b]))*((1.+pow(chiQ[0],2.))/pow(chiQ[0],2.)),0.5); - } - //cout<<"v'_2/2{2} = "<SetBinContent(b+1,v2eta[b]); - fdfr2->SetBinError(b+1,sdDiff2eta[b]); - //common histogram: - //fchr2nd->FillDifferentialFlow(b+1,v2eta[b],sdDiff2eta[b]); - // Fill common result histogram: - if(TMath::Abs(v2eta[b])>1.e-44) fchr2nd->FillDifferentialFlowEtaPOI(b+1,v2eta[b],sdDiff2eta[b]); - - } else { - //cout<<"v'_2/2{2} = Im"<0.&&etaBinPOINoOfParticles[b]>0.) // to be improved - { - if(chiQ[1]>0) - { - sdDiff4eta[b]=pow((1./(2.*etaBinPOINoOfParticles[b]))*((2.+6.*pow(chiQ[1],2.)+pow(chiQ[1],4.)+pow(chiQ[1],6.))/pow(chiQ[1],6.)),0.5); - } - //cout<<"v'_2/2{4} = "<SetBinContent(b+1,v4eta[b]); - fdfr4->SetBinError(b+1,sdDiff4eta[b]); - //common histogram: - //fchr4th->FillDifferentialFlow(b+1,v4eta[b],sdDiff4eta[b]); - // Fill common result histogram: - if(TMath::Abs(v4eta[b])>1.e-44) fchr4th->FillDifferentialFlowEtaPOI(b+1,v4eta[b],sdDiff4eta[b]); - - } else { - //cout<<"v'_2/2{4} = Im"<0){ - //cout<<"v'_2/2{6} = "<<100.*etaDiffCumulant6[b]/(4.*pow((1./4.)*cumulant[2],(5./6.)))<<"%"<SetBinContent(b+1,v6eta[b]); - //common histogram: - //fchr6th->FillDifferentialFlow(b+1,v6eta[b],0.); - // Fill common result histogram: - if(TMath::Abs(v6eta[b])>1.e-44) fchr6th->FillDifferentialFlowEtaPOI(b+1,v6eta[b],0.); - - }else{ - //cout<<"v'_2/2{6} = Im"<SetBinContent(b+1,v8eta[b]); - //common histogram: - //fchr8th->FillDifferentialFlow(b+1,v8eta[b],0.); - // Fill common result histogram: - if(TMath::Abs(v8eta[b])>1.e-44) fchr8th->FillDifferentialFlowEtaPOI(b+1,v8eta[b],0.); - - }else{ - //cout<<"v'_2/2{8} = Im"<GetBinContent(1); - } - - //number of events - Int_t nEvents4=0; - if(fAvMult4) - { - nEvents4=(Int_t)(fAvMult4->GetBinEntries(1)); - } - - TMatrixD dAvG4(cPmax4,cQmax4); - Bool_t someAvGEntryIsNegative4=kFALSE; - for(Int_t p=0;pGetBinContent(p+1,q+1); - if(dAvG4(p,q)<0.) - { - someAvGEntryIsNegative4=kTRUE; - } - } - } - ///////////////////////////////////////////////////////////////////////////// - //////////////////gen. function for the cumulants//////////////////////////// - ///////////////////////////////////////////////////////////////////////////// - - TMatrixD dC4(cPmax4,cQmax4);//C[p][q] - if(dAvM4>0 && someAvGEntryIsNegative4==kFALSE) - { - for (Int_t p=0;p - for (Int_t p=0;p=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(cumulant4[0],(1./2.))*dAvM4,2.)>0.)) - { - chiQo4[0]=dAvM4*dV2o4/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2o4*dAvM4,2.),0.5); - if(chiQo4[0]) - { - sdQo4[0]=pow(((1./(2.*dAvM4*nEvents4))*((1.+2.*pow(chiQo4[0],2))/(2.*pow(chiQo4[0],2)))),0.5); - } - cout<<" v_"<0.)) - { - chiQo4[1]=dAvM4*dV4o4/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4o4*dAvM4,2.),0.5); - if(chiQo4[1]) - { - sdQo4[1]=(1./(pow(2.*dAvM4*nEvents4,0.5)))*pow((1.+4.*pow(chiQo4[1],2)+1.*pow(chiQo4[1],4.)+2.*pow(chiQo4[1],6.))/(2.*pow(chiQo4[1],6.)),0.5); - } - cout<<" v_"<GetBinContent(1); - } - - //number of events - Int_t nEvents6=0; - if(fAvMult6) - { - nEvents6=(Int_t)(fAvMult6->GetBinEntries(1)); - } - - TMatrixD dAvG6(cPmax6,cQmax6); - Bool_t someAvGEntryIsNegative6=kFALSE; - for(Int_t p=0;pGetBinContent(p+1,q+1); - if(dAvG6(p,q)<0.) - { - someAvGEntryIsNegative6=kTRUE; - } - } - } - ///////////////////////////////////////////////////////////////////////////// - //////////////////gen. function for the cumulants//////////////////////////// - ///////////////////////////////////////////////////////////////////////////// - - TMatrixD dC6(cPmax6,cQmax6);//C[p][q] - if(dAvM6>0 && someAvGEntryIsNegative6==kFALSE) - { - for (Int_t p=0;p - Double_t tempHere6=0.; - for (Int_t p=0;p=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(cumulant6[0],(1./2.))*dAvM6,2.)>0.)) - { - chiQo6[0]=dAvM6*dV2o6/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2o6*dAvM6,2.),0.5); - if(chiQo6[0]) - { - sdQo6[0]=pow(((1./(2.*dAvM6*nEvents6))*((1.+2.*pow(chiQo6[0],2))/(2.*pow(chiQo6[0],2)))),0.5); - } - cout<<" v_"<0.)) - { - chiQo6[1]=dAvM6*dV4o6/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4o6*dAvM6,2.),0.5); - if(chiQo6[1]) - { - sdQo6[1]=(1./(pow(2.*dAvM6*nEvents6,0.5)))*pow((1.+4.*pow(chiQo6[1],2)+1.*pow(chiQo6[1],4.)+2.*pow(chiQo6[1],6.))/(2.*pow(chiQo6[1],6.)),0.5); - } - cout<<" v_"<=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow((1./4.)*cumulant6[2],(1./6.))*dAvM6,2.)>0.)) - { - chiQo6[2]=dAvM6*dV6/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV6o6*dAvM6,2.),0.5); - if(chiQo6[2]) - { - sdQo6[2]=(1./(pow(2.*dAvM6*nEvents6,0.5)))*pow((3.+18.*pow(chiQo6[2],2.)+9.*pow(chiQo6[2],4.)+28.*pow(chiQo6[2],6.)+12.*pow(chiQo6[2],8.)+24.*pow(chiQo6[2],10.))/(24.*pow(chiQo6[2],10.)),0.5); - } - cout<<" v_"<GetBinContent(1); - } - - //number of events - Int_t nEvents8=0; - if(fAvMult8) - { - nEvents8=(Int_t)(fAvMult8->GetBinEntries(1)); - } - - TMatrixD dAvG8(cPmax8,cQmax8); - Bool_t someAvGEntryIsNegative8=kFALSE; - for(Int_t p=0;pGetBinContent(p+1,q+1); - if(dAvG8(p,q)<0.) - { - someAvGEntryIsNegative8=kTRUE; - } - } - } - ///////////////////////////////////////////////////////////////////////////// - //////////////////gen. function for the cumulants//////////////////////////// - ///////////////////////////////////////////////////////////////////////////// - - TMatrixD dC8(cPmax8,cQmax8);//C[p][q] - if(dAvM8>0 && someAvGEntryIsNegative8==kFALSE) - { - for (Int_t p=0;p - Double_t tempHere8=0.; - for (Int_t p=0;p=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(cumulant8[0],(1./2.))*dAvM8,2.)>0.)) - { - chiQo8[0]=dAvM8*dV2o8/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2o8*dAvM8,2.),0.5); - if(chiQo8[0]) - { - sdQo8[0]=pow(((1./(2.*dAvM8*nEvents8))*((1.+2.*pow(chiQo8[0],2.))/(2.*pow(chiQo8[0],2)))),0.5); - } - cout<<" v_"<0.)) - { - chiQo8[1]=dAvM8*dV4o8/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4o8*dAvM8,2.),0.5); - if(chiQo8[1]) - { - sdQo8[1]=(1./(pow(2.*dAvM8*nEvents8,0.5)))*pow((1.+4.*pow(chiQo8[1],2)+1.*pow(chiQo8[1],4.)+2.*pow(chiQo8[1],6.))/(2.*pow(chiQo8[1],6.)),0.5); - } - cout<<" v_"<=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow((1./4.)*cumulant8[2],(1./6.))*dAvM8,2.)>0.)) - { - chiQo8[2]=dAvM8*dV6o8/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV6o8*dAvM8,2.),0.5); - if(chiQo8[2]) - { - sdQo8[2]=(1./(pow(2.*dAvM8*nEvents8,0.5)))*pow((3.+18.*pow(chiQo8[2],2)+9.*pow(chiQo8[2],4.)+28.*pow(chiQo8[2],6.)+12.*pow(chiQo8[2],8.)+24.*pow(chiQo8[2],10.))/(24.*pow(chiQo8[2],10.)),0.5); - } - cout<<" v_"<0.)) - { - chiQo8[3]=dAvM8*dV8o8/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV8o8*dAvM8,2.),0.5); - if(chiQo8[3]) - { - sdQo8[3]=(1./(pow(2.*dAvM8*nEvents8,0.5)))*pow((12.+96.*pow(chiQo8[3],2)+72.*pow(chiQo8[3],4.)+304.*pow(chiQo8[3],6.)+257.*pow(chiQo8[3],8.)+804.*pow(chiQo8[3],10.)+363.*pow(chiQo8[3],12.)+726.*pow(chiQo8[3],14.))/(726.*pow(chiQo8[3],14.)),0.5); - } - cout<<" v_"<GetBinContent(1); - } - - //number of events - Int_t nEvents16=0; - if(fAvMult16) - { - nEvents16=(Int_t)(fAvMult16->GetBinEntries(1)); - } - - TMatrixD dAvG16(cPmax16,cQmax16); - Bool_t someAvGEntryIsNegative16=kFALSE; - for(Int_t p=0;pGetBinContent(p+1,q+1); - if(dAvG16(p,q)<0.) - { - someAvGEntryIsNegative16=kTRUE; - } - } - } - - ///////////////////////////////////////////////////////////////////////////// - //////////////////gen. function for the cumulants//////////////////////////// - ///////////////////////////////////////////////////////////////////////////// - - TMatrixD dC16(cPmax16,cQmax16);//C16[p][q] - if(dAvM16>0 && someAvGEntryIsNegative16==kFALSE) - { - for(Int_t p=0;p - Double_t tempHere16=0.; - for (Int_t p=0;p=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(cumulant16[0],(1./2.))*dAvM16,2.)>0.)) - { - chiQo16[0]=dAvM16*dV2o16/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2o16*dAvM16,2.),0.5); - if(chiQo16[0]) - { - sdQo16[0]=pow(((1./(2.*dAvM16*nEvents16))*((1.+2.*pow(chiQo16[0],2))/(2.*pow(chiQo16[0],2)))),0.5); - } - cout<<" v_"<0.)) - { - chiQo16[1]=dAvM16*dV4o16/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4o16*dAvM16,2.),0.5); - if(chiQo16[1]) - { - sdQo16[1]=(1./(pow(2.*dAvM16*nEvents16,0.5)))*pow((1.+4.*pow(chiQo16[1],2.)+1.*pow(chiQo16[1],4.)+2.*pow(chiQo16[1],6.))/(2.*pow(chiQo16[1],6.)),0.5); - } - cout<<" v_"<=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow((1./4.)*cumulant16[2],(1./6.))*dAvM16,2.)>0.)) - { - chiQo16[2]=dAvM16*dV6o16/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV6o16*dAvM16,2.),0.5); - if(chiQo16[2]) - { - sdQo16[2]=(1./(pow(2.*dAvM16*nEvents16,0.5)))*pow((3.+18.*pow(chiQo16[2],2)+9.*pow(chiQo16[2],4.)+28.*pow(chiQo16[2],6.)+12.*pow(chiQo16[2],8.)+24.*pow(chiQo16[2],10.))/(24.*pow(chiQo16[2],10.)),0.5); - } - cout<<" v_"<0.)) - { - chiQo16[3]=dAvM16*dV8o16/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV8o16*dAvM16,2.),0.5); - if(chiQo16[3]) - { - sdQo16[3]=(1./(pow(2.*dAvM16*nEvents16,0.5)))*pow((12.+96.*pow(chiQo16[3],2)+72.*pow(chiQo16[3],4.)+304.*pow(chiQo16[3],6.)+257.*pow(chiQo16[3],8.)+804.*pow(chiQo16[3],10.)+363.*pow(chiQo16[3],12.)+726.*pow(chiQo16[3],14.))/(726.*pow(chiQo16[3],14.)),0.5); - } - cout<<" v_"<=0.) - { - cout<<"v_"<=0.) - { - cout<<"v_"<, 2nd bin: , 3rd bin: <(Q_x)^2>, 4th bin: <(Q_y)^2>) - - TProfile *fAverageOfSquaredWeight; // - - AliFlowCommonHistResults *fchr2nd; //final results for 2nd order int. and diff. flow stored in the common histograms - AliFlowCommonHistResults *fchr4th; //final results for 4th order int. and diff. flow stored in the common histograms - AliFlowCommonHistResults *fchr6th; //final results for 6th order int. and diff. flow stored in the common histograms - AliFlowCommonHistResults *fchr8th; //final results for 8th order int. and diff. flow stored in the common histograms - - AliFlowCommonHist *fch; //common control histogram - - /* - TProfile2D *fdRe0,*fdRe1,*fdRe2,*fdRe3,*fdRe4,*fdRe5,*fdRe6,*fdRe7;//differential flow - TProfile2D *fdIm0,*fdIm1,*fdIm2,*fdIm3,*fdIm4,*fdIm5,*fdIm6,*fdIm7;//differential flow - */ - - Int_t index3d(Int_t i, Int_t j, Int_t k, Int_t maxi, Int_t maxj) {return (i*maxj+j+maxi*maxj*k);} //for indexing 3d arrays - - ClassDef(AliCumulantsFunctions, 0); -}; - -//================================================================================================================ - -#endif - - - - - diff --git a/PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithCumulants.cxx b/PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithCumulants.cxx index 5ec79ed751f..35598533489 100644 --- a/PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithCumulants.cxx +++ b/PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithCumulants.cxx @@ -13,12 +13,17 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/******************************** - * flow analysis with cumulants * - * * - * author: Ante Bilandzic * - * (anteb@nikhef.nl) * - *******************************/ +/* $Id$ */ + +/************************************************* + * Flow analysis with cumulants. In this class * + * cumulants are calculated by making use of the * + * formalism of generating functions proposed by * + * Ollitrault et al. * + * * + * Author: Ante Bilandzic * + * (abilandzic@gmail.com) * + *************************************************/ #define AliFlowAnalysisWithCumulants_cxx @@ -29,9 +34,10 @@ #include "TProfile.h" #include "TProfile2D.h" #include "TProfile3D.h" -#include "TH1.h" #include "TH1F.h" #include "TH1D.h" +#include "TMatrixD.h" +#include "TVectorD.h" #include "AliFlowCommonConstants.h" #include "AliFlowCommonHist.h" @@ -39,127 +45,119 @@ #include "AliFlowEventSimple.h" #include "AliFlowTrackSimple.h" #include "AliFlowAnalysisWithCumulants.h" -#include "AliFlowCumuConstants.h" -#include "AliCumulantsFunctions.h" #include "AliFlowVector.h" - //================================================================================================================ ClassImp(AliFlowAnalysisWithCumulants) -AliFlowAnalysisWithCumulants::AliFlowAnalysisWithCumulants(): - fTrack(NULL), +AliFlowAnalysisWithCumulants::AliFlowAnalysisWithCumulants(): fHistList(NULL), - fWeightsList(NULL), - fR0(0), - fPtMax(0), - fPtMin(0), - fBinWidthPt(0), - fgknBinsPt(0), - fEtaMax(0), - fEtaMin(0), - fBinWidthEta(0), - fgknBinsEta(0), - fAvQx(0), - fAvQy(0), - fAvQ2x(0), - fAvQ2y(0), - fAvMultIntFlowGFC(NULL), - fQVectorComponentsGFC(NULL), - fIntFlowResultsGFC(NULL), - fDiffFlowResults2ndOrderGFC(NULL), - fDiffFlowResults4thOrderGFC(NULL), - fDiffFlowResults6thOrderGFC(NULL), - fDiffFlowResults8thOrderGFC(NULL), + fHistListName(NULL), + fAnalysisSettings(NULL), + fCommonHists(NULL), fCommonHistsResults2nd(NULL), fCommonHistsResults4th(NULL), fCommonHistsResults6th(NULL), fCommonHistsResults8th(NULL), - fIntFlowGenFun(NULL), - fIntFlowGenFun4(NULL),//(only for other system of Eq.) - fIntFlowGenFun6(NULL),//(only for other system of Eq.) - fIntFlowGenFun8(NULL),//(only for other system of Eq.) - fIntFlowGenFun16(NULL),//(only for other system of Eq.) - fAvMultIntFlow4GFC(NULL),//(only for other system of Eq.) - fAvMultIntFlow6GFC(NULL),//(only for other system of Eq.) - fAvMultIntFlow8GFC(NULL),//(only for other system of Eq.) - fAvMultIntFlow16GFC(NULL),//(only for other system of Eq.) - fDiffFlowPtRPGenFunRe(NULL), - fDiffFlowPtRPGenFunIm(NULL), - fPtBinRPNoOfParticles(NULL), - fDiffFlowEtaRPGenFunRe(NULL), - fDiffFlowEtaRPGenFunIm(NULL), - fEtaBinRPNoOfParticles(NULL), - fDiffFlowPtPOIGenFunRe(NULL), - fDiffFlowPtPOIGenFunIm(NULL), - fPtBinPOINoOfParticles(NULL), - fDiffFlowEtaPOIGenFunRe(NULL), - fDiffFlowEtaPOIGenFunIm(NULL), - fEtaBinPOINoOfParticles(NULL), - /* - fDiffFlowGenFunRe0(NULL), - fDiffFlowGenFunRe1(NULL), - fDiffFlowGenFunRe2(NULL), - fDiffFlowGenFunRe3(NULL), - fDiffFlowGenFunRe4(NULL), - fDiffFlowGenFunRe5(NULL), - fDiffFlowGenFunRe6(NULL), - fDiffFlowGenFunRe7(NULL), - fDiffFlowGenFunIm0(NULL), - fDiffFlowGenFunIm1(NULL), - fDiffFlowGenFunIm2(NULL), - fDiffFlowGenFunIm3(NULL), - fDiffFlowGenFunIm4(NULL), - fDiffFlowGenFunIm5(NULL), - fDiffFlowGenFunIm6(NULL), - fDiffFlowGenFunIm7(NULL), - */ - fCommonHists(NULL), - fOtherEquations(kFALSE), + fnBinsPhi(0), + fPhiMin(0), + fPhiMax(0), + fPhiBinWidth(0), + fnBinsPt(0), + fPtMin(0), + fPtMax(0), + fPtBinWidth(0), + fnBinsEta(0), + fEtaMin(0), + fEtaMax(0), + fEtaBinWidth(0), + fHarmonic(2), + fMultiple(1), + fR0(2.2), + fWeightsList(NULL), fUsePhiWeights(kFALSE), fUsePtWeights(kFALSE), fUseEtaWeights(kFALSE), - fAverageOfSquaredWeight(NULL) -{ - //constructor - fHistList = new TList(); - fWeightsList = new TList(); - fWeightsList->SetName("Weights"); - fWeightsList->SetOwner(kTRUE); - fR0=AliFlowCumuConstants::GetMaster()->GetR0(); - //Pt: - fPtMax=AliFlowCommonConstants::GetMaster()->GetPtMax(); - fPtMin=AliFlowCommonConstants::GetMaster()->GetPtMin(); - fgknBinsPt=AliFlowCommonConstants::GetMaster()->GetNbinsPt(); - if(fgknBinsPt) - { - fBinWidthPt=(fPtMax-fPtMin)/fgknBinsPt; - } - //Eta: - fEtaMax=AliFlowCommonConstants::GetMaster()->GetEtaMax(); - fEtaMin=AliFlowCommonConstants::GetMaster()->GetEtaMin(); - fgknBinsEta=AliFlowCommonConstants::GetMaster()->GetNbinsEta(); - if(fgknBinsEta) + fUseParticleWeights(NULL), + fPhiWeights(NULL), + fPtWeights(NULL), + fEtaWeights(NULL), + fMultiplicityWeight(NULL), + fReferenceFlowList(NULL), + fReferenceFlowProfiles(NULL), + fReferenceFlowResults(NULL), + fReferenceFlowFlags(NULL), + fnBinsMult(10000), + fMinMult(0.), + fMaxMult(10000.), + fGEBE(NULL), + fReferenceFlowGenFun(NULL), + fQvectorComponents(NULL), + fAverageOfSquaredWeight(NULL), + fAvM(0.), + fnEvts(0), + fReferenceFlowCumulants(NULL), + fReferenceFlow(NULL), + fChi(NULL), + fDiffFlowList(NULL), + fDiffFlowProfiles(NULL), + fDiffFlowResults(NULL), + fDiffFlowFlags(NULL), + fTuningList(NULL), + fTuningProfiles(NULL), + fTuningResults(NULL), + fTuningFlags(NULL), + fTuneParameters(kFALSE) { - fBinWidthEta=(fEtaMax-fEtaMin)/fgknBinsEta; - } + // Constructor. + + // Base list to hold all output objects: + fHistList = new TList(); + fHistListName = new TString("cobjGFC"); + fHistList->SetName(fHistListName->Data()); + fHistList->SetOwner(kTRUE); + + // List to hold histograms with phi, pt and eta weights: + fWeightsList = new TList(); + fWeightsList->SetName("Weights"); + fWeightsList->SetOwner(kTRUE); + fHistList->Add(fWeightsList); + + // Multiplicity weight: + fMultiplicityWeight = new TString("unit"); + + // Initialize all arrays: + this->InitializeArrays(); - fOtherEquations=AliFlowCumuConstants::GetMaster()->GetOtherEquations(); -} +} // end of AliFlowAnalysisWithCumulants::AliFlowAnalysisWithCumulants() + +//================================================================================================================ AliFlowAnalysisWithCumulants::~AliFlowAnalysisWithCumulants() { - //desctructor + // Desctructor. + delete fHistList; delete fWeightsList; -} + +} // end of AliFlowAnalysisWithCumulants::~AliFlowAnalysisWithCumulants() //================================================================================================================ void AliFlowAnalysisWithCumulants::Init() { - //various output histograms + // Initialize and book all objects. + + // a) Cross check if the user settings make sense before starting; + // b) Access all common constants; + // c) Book and fill weights histograms; + // d) Book and nest all lists in the base list fHistList; + // e) Book and fill profile holding analysis settings; + // f) Book common control and results histograms; + // g) Store flags for reference flow; + // h) Store flags for differential flow; + // i) Book all objects needed for tuning. //save old value and prevent histograms from being added to directory //to avoid name clashes in case multiple analaysis objects are used @@ -167,1023 +165,2701 @@ void AliFlowAnalysisWithCumulants::Init() Bool_t oldHistAddStatus = TH1::AddDirectoryStatus(); TH1::AddDirectory(kFALSE); - //average multiplicity - fAvMultIntFlowGFC = new TProfile("fAvMultIntFlowGFC","Average Weighted Multiplicity",1,0,1,"s"); - fAvMultIntFlowGFC->SetXTitle(""); - fAvMultIntFlowGFC->SetYTitle(""); - fAvMultIntFlowGFC->SetLabelSize(0.06); - fAvMultIntFlowGFC->SetMarkerStyle(25); - fAvMultIntFlowGFC->SetLabelOffset(0.01); - (fAvMultIntFlowGFC->GetXaxis())->SetBinLabel(1,"Average Weighted Multiplicity"); - fHistList->Add(fAvMultIntFlowGFC); - - //averages of Q-vector components (1st bin: , 2nd bin: , 3rd bin: <(Q_x)^2>, 4th bin: <(Q_y)^2>) - fQVectorComponentsGFC = new TProfile("fQVectorComponentsGFC","Average of Q-vector components",4,0.,4.); - fQVectorComponentsGFC->SetXTitle(""); - fQVectorComponentsGFC->SetYTitle(""); - fQVectorComponentsGFC->SetLabelSize(0.06); - fQVectorComponentsGFC->SetMarkerStyle(25); - (fQVectorComponentsGFC->GetXaxis())->SetBinLabel(1,""); - (fQVectorComponentsGFC->GetXaxis())->SetBinLabel(2,""); - (fQVectorComponentsGFC->GetXaxis())->SetBinLabel(3,""); - (fQVectorComponentsGFC->GetXaxis())->SetBinLabel(4,""); - fHistList->Add(fQVectorComponentsGFC); - - //final results for integrated flow (v_n{2}, v_n{4},..., v_n{16}) from cumulants (by default n=2) - fIntFlowResultsGFC = new TH1D("fIntFlowResultsGFC","Integrated Flow From Cumulants (Generating Function)",4,0,4); - fIntFlowResultsGFC->SetXTitle(""); - fIntFlowResultsGFC->SetYTitle(""); - fIntFlowResultsGFC->SetLabelSize(0.06); - //fIntFlowResultsGFC->SetTickLength(1); - fIntFlowResultsGFC->SetMarkerStyle(25); - (fIntFlowResultsGFC->GetXaxis())->SetBinLabel(1,"v_{n}{2}"); - (fIntFlowResultsGFC->GetXaxis())->SetBinLabel(2,"v_{n}{4}"); - (fIntFlowResultsGFC->GetXaxis())->SetBinLabel(3,"v_{n}{6}"); - (fIntFlowResultsGFC->GetXaxis())->SetBinLabel(4,"v_{n}{8}"); - fHistList->Add(fIntFlowResultsGFC); - - //final results for differential flow v_p/n{2} (by default p=n=2) - fDiffFlowResults2ndOrderGFC = new TH1D("fDiffFlowResults2ndOrderGFC","v'_2/2{2}",fgknBinsPt,fPtMin,fPtMax); - fDiffFlowResults2ndOrderGFC->SetXTitle("pt [GeV]"); - fDiffFlowResults2ndOrderGFC->SetYTitle(""); - fHistList->Add(fDiffFlowResults2ndOrderGFC); - - //final results for differential flow v_p/n{4} (by default p=n=2) - fDiffFlowResults4thOrderGFC = new TH1D("fDiffFlowResults4thOrderGFC","v'_2/2{4}",fgknBinsPt,fPtMin,fPtMax); - fDiffFlowResults4thOrderGFC->SetXTitle("pt [GeV]"); - fDiffFlowResults4thOrderGFC->SetYTitle(""); - fHistList->Add(fDiffFlowResults4thOrderGFC); - - //final results for differential flow v_p/n{6} (by default p=n=2) - fDiffFlowResults6thOrderGFC = new TH1D("fDiffFlowResults6thOrderGFC","v'_2/2{6}",fgknBinsPt,fPtMin,fPtMax); - fDiffFlowResults6thOrderGFC->SetXTitle("pt [GeV]"); - fDiffFlowResults6thOrderGFC->SetYTitle(""); - fHistList->Add(fDiffFlowResults6thOrderGFC); - - //final results for differential flow v_p/n{8} (by default p=n=2) - fDiffFlowResults8thOrderGFC = new TH1D("fDiffFlowResults8thOrderGFC","v'_2/2{8}",fgknBinsPt,fPtMin,fPtMax); - fDiffFlowResults8thOrderGFC->SetXTitle("pt [GeV]"); - fDiffFlowResults8thOrderGFC->SetYTitle(""); - fHistList->Add(fDiffFlowResults8thOrderGFC); - - //avarage of the generating function for integrated flow - fIntFlowGenFun = new TProfile2D("fIntFlowGenFun","",fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax); - fIntFlowGenFun->SetXTitle("p"); - fIntFlowGenFun->SetYTitle("q"); - fHistList->Add(fIntFlowGenFun); - - if(fOtherEquations) - { - //avarage of the generating function for integrated flow (only for other system of Eq. - up to 4th order) - fIntFlowGenFun4 = new TProfile2D("fIntFlowGenFun4","",fgkPmax4,0.,(Double_t)fgkPmax4,fgkQmax4,0.,(Double_t)fgkQmax4); - fIntFlowGenFun4->SetXTitle("p4"); - fIntFlowGenFun4->SetYTitle("q4"); - fHistList->Add(fIntFlowGenFun4); - - //avarage of the generating function for integrated flow (only for other system of Eq. - up to 6th order) - fIntFlowGenFun6 = new TProfile2D("fIntFlowGenFun6","",fgkPmax6,0.,(Double_t)fgkPmax6,fgkQmax6,0.,(Double_t)fgkQmax6); - fIntFlowGenFun6->SetXTitle("p6"); - fIntFlowGenFun6->SetYTitle("q6"); - fHistList->Add(fIntFlowGenFun6); - - //avarage of the generating function for integrated flow (only for other system of Eq. - up to 8th order) - fIntFlowGenFun8 = new TProfile2D("fIntFlowGenFun8","",fgkPmax8,0.,(Double_t)fgkPmax8,fgkQmax8,0.,(Double_t)fgkQmax8); - fIntFlowGenFun8->SetXTitle("p8"); - fIntFlowGenFun8->SetYTitle("q8"); - fHistList->Add(fIntFlowGenFun8); - - //avarage of the generating function for integrated flow (only for other system of Eq. - up to 16th order) - fIntFlowGenFun16 = new TProfile2D("fIntFlowGenFun16","",fgkPmax16,0.,(Double_t)fgkPmax16,fgkQmax16,0.,(Double_t)fgkQmax16); - fIntFlowGenFun16->SetXTitle("p16"); - fIntFlowGenFun16->SetYTitle("q16"); - fHistList->Add(fIntFlowGenFun16); - - //average multiplicity (only for other system of Eq. - up to 4th order) - fAvMultIntFlow4GFC = new TProfile("fAvMultIntFlow4GFC","Average Multiplicity",1,0,1,"s"); - fAvMultIntFlow4GFC->SetXTitle(""); - fAvMultIntFlow4GFC->SetYTitle(""); - fAvMultIntFlow4GFC->SetLabelSize(0.06); - fAvMultIntFlow4GFC->SetMarkerStyle(25); - fAvMultIntFlow4GFC->SetLabelOffset(0.01); - (fAvMultIntFlow4GFC->GetXaxis())->SetBinLabel(1,"Average Multiplicity"); - fHistList->Add(fAvMultIntFlow4GFC); - - //average multiplicity (only for other system of Eq. - up to 6th order) - fAvMultIntFlow6GFC = new TProfile("fAvMultIntFlow6GFC","Average Multiplicity",1,0,1,"s"); - fAvMultIntFlow6GFC->SetXTitle(""); - fAvMultIntFlow6GFC->SetYTitle(""); - fAvMultIntFlow6GFC->SetLabelSize(0.06); - fAvMultIntFlow6GFC->SetMarkerStyle(25); - fAvMultIntFlow6GFC->SetLabelOffset(0.01); - (fAvMultIntFlow6GFC->GetXaxis())->SetBinLabel(1,"Average Multiplicity"); - fHistList->Add(fAvMultIntFlow6GFC); - - //average multiplicity (only for other system of Eq. - up to 8th order) - fAvMultIntFlow8GFC = new TProfile("fAvMultIntFlow8GFC","Average Multiplicity",1,0,1,"s"); - fAvMultIntFlow8GFC->SetXTitle(""); - fAvMultIntFlow8GFC->SetYTitle(""); - fAvMultIntFlow8GFC->SetLabelSize(0.06); - fAvMultIntFlow8GFC->SetMarkerStyle(25); - fAvMultIntFlow8GFC->SetLabelOffset(0.01); - (fAvMultIntFlow8GFC->GetXaxis())->SetBinLabel(1,"Average Multiplicity"); - fHistList->Add(fAvMultIntFlow8GFC); - - //average multiplicity (only for other system of Eq. - up to 16th order) - fAvMultIntFlow16GFC = new TProfile("fAvMultIntFlow16GFC","Average Multiplicity",1,0,1,"s"); - fAvMultIntFlow16GFC->SetXTitle(""); - fAvMultIntFlow16GFC->SetYTitle(""); - fAvMultIntFlow16GFC->SetLabelSize(0.06); - fAvMultIntFlow16GFC->SetMarkerStyle(25); - fAvMultIntFlow16GFC->SetLabelOffset(0.01); - (fAvMultIntFlow16GFC->GetXaxis())->SetBinLabel(1,"Average Multiplicity"); - fHistList->Add(fAvMultIntFlow16GFC); - } - - //avarage of the real part of generating function for differential flow in Pt - fDiffFlowPtRPGenFunRe = new TProfile3D("fDiffFlowPtRPGenFunRe","",fgknBinsPt,(Double_t)(fPtMin/fBinWidthPt),(Double_t)(fPtMax/fBinWidthPt),fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowPtRPGenFunRe->SetXTitle("b"); - fDiffFlowPtRPGenFunRe->SetYTitle("p"); - fDiffFlowPtRPGenFunRe->SetZTitle("q"); - fDiffFlowPtRPGenFunRe->SetTitleOffset(1.44,"X"); - fDiffFlowPtRPGenFunRe->SetTitleOffset(1.44,"Y"); - fHistList->Add(fDiffFlowPtRPGenFunRe); - - //avarage of the imaginary part of generating function for differential flow in Pt - fDiffFlowPtRPGenFunIm = new TProfile3D("fDiffFlowPtRPGenFunIm","",fgknBinsPt,(Double_t)(fPtMin/fBinWidthPt),(Double_t)(fPtMax/fBinWidthPt),fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowPtRPGenFunIm->SetXTitle("b"); - fDiffFlowPtRPGenFunIm->SetYTitle("p"); - fDiffFlowPtRPGenFunIm->SetZTitle("q"); - fDiffFlowPtRPGenFunIm->SetTitleOffset(1.44,"X"); - fDiffFlowPtRPGenFunIm->SetTitleOffset(1.44,"Y"); - fHistList->Add(fDiffFlowPtRPGenFunIm); - - //number of particles per pt bin - fPtBinRPNoOfParticles = new TProfile("fPtBinRPNoOfParticles","Number of particles per #p_{t} bin",fgknBinsPt,fPtMin,fPtMax); - fPtBinRPNoOfParticles->SetXTitle("pt [GeV]"); - fPtBinRPNoOfParticles->SetYTitle(""); - fHistList->Add(fPtBinRPNoOfParticles); - - //avarage of the real part of generating function for differential flow in Eta - fDiffFlowEtaRPGenFunRe = new TProfile3D("fDiffFlowEtaRPGenFunRe","",fgknBinsEta,(Double_t)(fEtaMin/fBinWidthEta),(Double_t)(fEtaMax/fBinWidthEta),fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowEtaRPGenFunRe->SetXTitle("b"); - fDiffFlowEtaRPGenFunRe->SetYTitle("p"); - fDiffFlowEtaRPGenFunRe->SetZTitle("q"); - fDiffFlowEtaRPGenFunRe->SetTitleOffset(1.44,"X"); - fDiffFlowEtaRPGenFunRe->SetTitleOffset(1.44,"Y"); - fHistList->Add(fDiffFlowEtaRPGenFunRe); - - //avarage of the imaginary part of generating function for differential flow in Eta - fDiffFlowEtaRPGenFunIm = new TProfile3D("fDiffFlowEtaRPGenFunIm","",fgknBinsEta,(Double_t)(fEtaMin/fBinWidthEta),(Double_t)(fEtaMax/fBinWidthEta),fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowEtaRPGenFunIm->SetXTitle("b"); - fDiffFlowEtaRPGenFunIm->SetYTitle("p"); - fDiffFlowEtaRPGenFunIm->SetZTitle("q"); - fDiffFlowEtaRPGenFunIm->SetTitleOffset(1.44,"X"); - fDiffFlowEtaRPGenFunIm->SetTitleOffset(1.44,"Y"); - fHistList->Add(fDiffFlowEtaRPGenFunIm); - - //number of particles per eta bin - fEtaBinRPNoOfParticles = new TProfile("fEtaBinRPNoOfParticles","Number of particles per #eta bin",fgknBinsEta,fEtaMin,fEtaMax); - fEtaBinRPNoOfParticles->SetXTitle("#eta"); - fEtaBinRPNoOfParticles->SetYTitle(""); - fHistList->Add(fEtaBinRPNoOfParticles); - - //avarage of the real part of generating function for differential flow in Pt - fDiffFlowPtPOIGenFunRe = new TProfile3D("fDiffFlowPtPOIGenFunRe","",fgknBinsPt,(Double_t)(fPtMin/fBinWidthPt),(Double_t)(fPtMax/fBinWidthPt),fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowPtPOIGenFunRe->SetXTitle("b"); - fDiffFlowPtPOIGenFunRe->SetYTitle("p"); - fDiffFlowPtPOIGenFunRe->SetZTitle("q"); - fDiffFlowPtPOIGenFunRe->SetTitleOffset(1.44,"X"); - fDiffFlowPtPOIGenFunRe->SetTitleOffset(1.44,"Y"); - fHistList->Add(fDiffFlowPtPOIGenFunRe); - - //avarage of the imaginary part of generating function for differential flow in Pt - fDiffFlowPtPOIGenFunIm = new TProfile3D("fDiffFlowPtPOIGenFunIm","",fgknBinsPt,(Double_t)(fPtMin/fBinWidthPt),(Double_t)(fPtMax/fBinWidthPt),fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowPtPOIGenFunIm->SetXTitle("b"); - fDiffFlowPtPOIGenFunIm->SetYTitle("p"); - fDiffFlowPtPOIGenFunIm->SetZTitle("q"); - fDiffFlowPtPOIGenFunIm->SetTitleOffset(1.44,"X"); - fDiffFlowPtPOIGenFunIm->SetTitleOffset(1.44,"Y"); - fHistList->Add(fDiffFlowPtPOIGenFunIm); - - //number of particles per pt bin - fPtBinPOINoOfParticles = new TProfile("fPtBinPOINoOfParticles","Number of particles per #p_{t} bin",fgknBinsPt,fPtMin,fPtMax); - fPtBinPOINoOfParticles->SetXTitle("pt [GeV]"); - fPtBinPOINoOfParticles->SetYTitle(""); - fHistList->Add(fPtBinPOINoOfParticles); - - //avarage of the real part of generating function for differential flow in Eta - fDiffFlowEtaPOIGenFunRe = new TProfile3D("fDiffFlowEtaPOIGenFunRe","",fgknBinsEta,(Double_t)(fEtaMin/fBinWidthEta),(Double_t)(fEtaMax/fBinWidthEta),fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowEtaPOIGenFunRe->SetXTitle("b"); - fDiffFlowEtaPOIGenFunRe->SetYTitle("p"); - fDiffFlowEtaPOIGenFunRe->SetZTitle("q"); - fDiffFlowEtaPOIGenFunRe->SetTitleOffset(1.44,"X"); - fDiffFlowEtaPOIGenFunRe->SetTitleOffset(1.44,"Y"); - fHistList->Add(fDiffFlowEtaPOIGenFunRe); - - //avarage of the imaginary part of generating function for differential flow in Eta - fDiffFlowEtaPOIGenFunIm = new TProfile3D("fDiffFlowEtaPOIGenFunIm","",fgknBinsEta,(Double_t)(fEtaMin/fBinWidthEta),(Double_t)(fEtaMax/fBinWidthEta),fgkPmax,0.,(Double_t)fgkPmax,fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowEtaPOIGenFunIm->SetXTitle("b"); - fDiffFlowEtaPOIGenFunIm->SetYTitle("p"); - fDiffFlowEtaPOIGenFunIm->SetZTitle("q"); - fDiffFlowEtaPOIGenFunIm->SetTitleOffset(1.44,"X"); - fDiffFlowEtaPOIGenFunIm->SetTitleOffset(1.44,"Y"); - fHistList->Add(fDiffFlowEtaPOIGenFunIm); - - //number of particles per eta bin - fEtaBinPOINoOfParticles = new TProfile("fEtaBinPOINoOfParticles","Number of particles per #eta bin",fgknBinsEta,fEtaMin,fEtaMax); - fEtaBinPOINoOfParticles->SetXTitle("#eta"); - fEtaBinPOINoOfParticles->SetYTitle(""); - fHistList->Add(fEtaBinPOINoOfParticles); - - /* - fDiffFlowGenFunRe0 = new TProfile2D("fDiffFlowGenFunRe0","Re()",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowGenFunRe0->SetXTitle("b"); - fDiffFlowGenFunRe0->SetYTitle("q"); - fHistList->Add(fDiffFlowGenFunRe0); - - fDiffFlowGenFunIm0 = new TProfile2D("fDiffFlowGcout<<"HEY M1"<)",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowGenFunIm0->SetXTitle("b"); - fDiffFlowGenFunIm0->SetYTitle("q"); - fHistList->Add(fDiffFlowGenFunIm0); - - fDiffFlowGenFunRe1 = new TProfile2D("fDiffFlowGenFunRe1","Re()",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowGenFunRe1->SetXTitle("b"); - fDiffFlowGenFunRe1->SetYTitle("q"); - fHistList->Add(fDiffFlowGenFunRe1); - - fDiffFlowGenFunIm1 = new TProfile2D("fDiffFlowGenFunIm1","Im()",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowGenFunIm1->SetXTitle("b"); - fDiffFlowGenFunIm1->SetYTitle("q"); - fHistList->Add(fDiffFlowGenFunIm1); - - fDiffFlowGenFunRe2 = new TProfile2D("fDiffFlowGenFunRe2","Re()",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowGenFunRe2->SetXTitle("b"); - fDiffFlowGenFunRe2->SetYTitle("q"); - fHistList->Add(fDiffFlowGenFunRe2); - - fDiffFlowGenFunIm2 = new TProfile2D("fDiffFlowGenFunIm2","Im()",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowGenFunIm2->SetXTitle("b"); - fDiffFlowGenFunIm2->SetYTitle("q"); - fHistList->Add(fDiffFlowGenFunIm2); - - fDiffFlowGenFunRe3 = new TProfile2D("fDiffFlowGenFunRe3","Re()",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowGenFunRe3->SetXTitle("b"); - fDiffFlowGenFunRe3->SetYTitle("q"); - fHistList->Add(fDiffFlowGenFunRe3); - - fDiffFlowGenFunIm3 = new TProfile2D("fDiffFlowGenFunIm3","Im()",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowGenFunIm3->SetXTitle("b"); - fDiffFlowGenFunIm3->SetYTitle("q"); - fHistList->Add(fDiffFlowGenFunIm3); - - fDiffFlowGenFunRe4 = new TProfile2D("fDiffFlowGenFunRe4","Re()",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowGenFunRe4->SetXTitle("b"); - fDiffFlowGenFunRe4->SetYTitle("q"); - fHistList->Add(fDiffFlowGenFunRe4); - - fDiffFlowGenFunIm4 = new TProfile2D("fDiffFlowGenFunIm4","Im()",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowGenFunIm4->SetXTitle("b"); - fDiffFlowGenFunIm4->SetYTitle("q"); - fHistList->Add(fDiffFlowGenFunIm4); - - fDiffFlowGenFunRe5 = new TProfile2D("fDiffFlowGenFunRe5","Re()",fgkQmax,0.,(Double_t)fgkQmax,fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth)); - fDiffFlowGenFunRe5->SetXTitle("b"); - fDiffFlowGenFunRe5->SetYTitle("q"); - fHistList->Add(fDiffFlowGenFunRe5); - - fDiffFlowGenFunIm5 = new TProfile2D("fDiffFlowGenFunIm5","Im()",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowGenFunIm5->SetXTitle("b"); - fDiffFlowGenFunIm5->SetYTitle("q"); - fHistList->Add(fDiffFlowGenFunIm5); - - fDiffFlowGenFunRe6 = new TProfile2D("fDiffFlowGenFunRe6","Re()",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowGenFunRe6->SetXTitle("b"); - fDiffFlowGenFunRe6->SetYTitle("q"); - fHistList->Add(fDiffFlowGenFunRe6); - - fDiffFlowGenFunIm6 = new TProfile2D("fDiffFlowGenFunIm6","Im()",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowGenFunIm6->SetXTitle("b"); - fDiffFlowGenFunIm6->SetYTitle("q"); - fHistList->Add(fDiffFlowGenFunIm6); - - fDiffFlowGenFunRe7 = new TProfile2D("fDiffFlowGenFunRe7","Re()",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowGenFunRe7->SetXTitle("b"); - fDiffFlowGenFunRe7->SetYTitle("q"); - fHistList->Add(fDiffFlowGenFunRe7); - - fDiffFlowGenFunIm7 = new TProfile2D("fDiffFlowGenFunIm7","Im()",fgknBins,(Double_t)(fPtMin/fBinWidth),(Double_t)(fPtMax/fBinWidth),fgkQmax,0.,(Double_t)fgkQmax); - fDiffFlowGenFunIm7->SetXTitle("b"); - fDiffFlowGenFunIm7->SetYTitle("q"); - fHistList->Add(fDiffFlowGenFunIm7); - */ - - //common control histograms - fCommonHists = new AliFlowCommonHist("AliFlowCommonHistGFC"); - fHistList->Add(fCommonHists); - - //common histograms for final results (2nd order) - fCommonHistsResults2nd = new AliFlowCommonHistResults("AliFlowCommonHistResults2ndOrderGFC"); - fHistList->Add(fCommonHistsResults2nd); + this->CrossCheckSettings(); + this->AccessConstants(); + this->BookAndFillWeightsHistograms(); + this->BookAndNestAllLists(); + this->BookProfileHoldingSettings(); + this->BookCommonHistograms(); + this->BookEverythingForReferenceFlow(); + this->BookEverythingForDiffFlow(); + this->StoreReferenceFlowFlags(); + this->StoreDiffFlowFlags(); + if(fTuneParameters){this->BookEverythingForTuning();} - //common histograms for final results (4th order) - fCommonHistsResults4th = new AliFlowCommonHistResults("AliFlowCommonHistResults4thOrderGFC"); - fHistList->Add(fCommonHistsResults4th); + (fCommonHists->GetHarmonic())->Fill(0.5,fHarmonic); // to be improved (moved somewhere else?) - //common histograms for final results (6th order) - fCommonHistsResults6th = new AliFlowCommonHistResults("AliFlowCommonHistResults6thOrderGFC"); - fHistList->Add(fCommonHistsResults6th); - - //common histograms for final results (8th order) - fCommonHistsResults8th = new AliFlowCommonHistResults("AliFlowCommonHistResults8thOrderGFC"); - fHistList->Add(fCommonHistsResults8th); - - // - fAverageOfSquaredWeight = new TProfile("fAverageOfSquaredWeight","",1,0,1); - fAverageOfSquaredWeight->SetLabelSize(0.06); - fAverageOfSquaredWeight->SetMarkerStyle(25); - fAverageOfSquaredWeight->SetLabelOffset(0.01); - (fAverageOfSquaredWeight->GetXaxis())->SetBinLabel(1,""); - fHistList->Add(fAverageOfSquaredWeight); - - // add list fWeightsList with weights to the main list - fHistList->Add(fWeightsList); - TH1::AddDirectory(oldHistAddStatus); -}//end of Init() + +} // end of void AliFlowAnalysisWithCumulants::Init() //================================================================================================================ void AliFlowAnalysisWithCumulants::Make(AliFlowEventSimple* anEvent) { - //running over data: - Int_t nPrim = anEvent->NumberOfTracks(); //total multiplicity + // Running over data only in this method. - Int_t nEventNSelTracksRP = anEvent->GetEventNSelTracksRP(); //selected multiplicity (particles used for int. flow) + // a) Check all pointers used in this method; + // b) If tuning enabled, fill generating functions for different values of tuning parameters; + // c) For default values of tuning parameters (r0 = 2.2 and cutoff at 10th order): + // c1) Fill common control histograms; + // c2) Fill generating function for reference flow; + // c3) Fill profile holding average values of various Q-vector components; + // c4) Fill generating function for differential flow. - Int_t n = 2; // int flow harmonic (to be improved) + this->CheckPointersUsedInMake(); + if(fTuneParameters) {this->FillGeneratingFunctionsForDifferentTuningParameters(anEvent);} + Int_t nRP = anEvent->GetEventNSelTracksRP(); // number of RPs (i.e. number of particles used to determine the reaction plane) + if(nRP<10) {return;} // generating function formalism make sense only for nRPs >= 10 for default settings + fCommonHists->FillControlHistograms(anEvent); + this->FillGeneratingFunctionForReferenceFlow(anEvent); + this->FillQvectorComponents(anEvent); + this->FillGeneratingFunctionForDiffFlow(anEvent); + +} // end of void AliFlowAnalysisWithCumulants::Make() + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::Finish() +{ + // Calculate the final results. - //--------------------------------------------------------------------------------------------------------- - // weights: - Bool_t useWeights = fUsePhiWeights||fUsePtWeights||fUseEtaWeights; + // a) Check all pointers used in this method; + // b) Access settings for analysis with Generating Function Cumulants; + // c) From relevant common control histogram get average multiplicity of RPs and number of events; + // d) Calculate cumulants for reference flow; + // e) Calculate from isotropic cumulants reference flow; + // f) Calculate error for reference flow estimates; + // g) Store the final results for reference flow in common hist results; + // h) Print on the screen the final results for reference flow; + // i) Calculate cumulants for differential flow; + // j) Calculate differential flow for RPs/POIs vs pt/eta from cumulants; + // k) Calculate integrated flow of RPs and POIs; + // l) Print on the screen the final results for integrated flow of RPs and POIs; + // m) If tuning enabled, calculate results for different tuning parameters. - TH1F *phiWeights = NULL; // histogram with phi weights - TH1D *ptWeights = NULL; // histogram with pt weights - TH1D *etaWeights = NULL; // histogram with eta weights + this->CheckPointersUsedInFinish(); + this->AccessSettings(); + this->GetAvMultAndNoOfEvts(); + this->CalculateCumulantsForReferenceFlow(); + this->CalculateReferenceFlow(); + this->CalculateReferenceFlowError(); + this->FillCommonHistResultsForReferenceFlow(); + if(fPrintFinalResults[0]){this->PrintFinalResults("RF");} + this->CalculateCumulantsForDiffFlow("RP","pt"); + this->CalculateCumulantsForDiffFlow("RP","eta"); + this->CalculateCumulantsForDiffFlow("POI","pt"); + this->CalculateCumulantsForDiffFlow("POI","eta"); + this->CalculateDifferentialFlow("RP","pt"); + this->CalculateDifferentialFlow("RP","eta"); + this->CalculateDifferentialFlow("POI","pt"); + this->CalculateDifferentialFlow("POI","eta"); + this->CalculateDifferentialFlowErrors("RP","pt"); + this->CalculateDifferentialFlowErrors("RP","eta"); + this->CalculateDifferentialFlowErrors("POI","pt"); + this->CalculateDifferentialFlowErrors("POI","eta"); + this->FillCommonHistResultsForDifferentialFlow("RP"); + this->FillCommonHistResultsForDifferentialFlow("POI"); + this->CalculateIntegratedFlow("RP"); + this->CalculateIntegratedFlow("POI"); + if(fPrintFinalResults[1]){this->PrintFinalResults("RP");} + if(fPrintFinalResults[2]){this->PrintFinalResults("POI");} + if(fTuneParameters){this->FinalizeTuning();} + +} // end of void AliFlowAnalysisWithCumulants::Finish() + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::FinalizeTuning() +{ + // Finalize results with tuned inerpolating parameters. + + for(Int_t r=0;r<10;r++) + { + if(TMath::Abs(fTuningR0[r])<1.e-10) continue; // protection against division by r0 bellow + for(Int_t pq=0;pq<5;pq++) + { + Int_t pMax = fTuningGenFun[r][pq]->GetXaxis()->GetNbins(); + Int_t qMax = fTuningGenFun[r][pq]->GetYaxis()->GetNbins(); + // + TMatrixD dAvG(pMax,qMax); + dAvG.Zero(); + Bool_t someAvGEntryIsNegative = kFALSE; + for(Int_t p=0;pGetBinContent(fTuningGenFun[r][pq]->GetBin(p+1,q+1)); + if(dAvG(p,q)<0.) + { + someAvGEntryIsNegative = kTRUE; + cout< is negative !!!! GFC results are meaningless for r0 = %f, pq = %i.",p,q,fTuningR0[r],pq)<0. && !someAvGEntryIsNegative) + { + for(Int_t p=0;p (Remark: here <> stands for average over azimuth): + TVectorD dAvC(pMax); + dAvC.Zero(); + for(Int_t p=0;p=0.) {flow[0] = pow(cumulant[0],1./2.);} + if(cumulant[1]<=0.) {flow[1] = pow(-1.*cumulant[1],1./4.);} + } + else if(pMax==3) + { + cumulant[0] = (1./(fTuningR0[r]*fTuningR0[r]))*(3.*dAvC[0]-(3./2.)*dAvC[1]+(1./3.)*dAvC[2]); + cumulant[1] = (2./pow(fTuningR0[r],4.))*((-5.)*dAvC[0]+4.*dAvC[1]-1.*dAvC[2]); + cumulant[2] = (6./pow(fTuningR0[r],6.))*(3.*dAvC[0]-3.*dAvC[1]+1.*dAvC[2]); + if(cumulant[0]>=0.) {flow[0] = pow(cumulant[0],1./2.);} + if(cumulant[1]<=0.) {flow[1] = pow(-1.*cumulant[1],1./4.);} + if(cumulant[2]>=0.) {flow[2] = pow((1./4.)*cumulant[2],1./6.);} + } + else if(pMax==4) + { + cumulant[0] = (1./(fTuningR0[r]*fTuningR0[r]))*(4.*dAvC[0]-3.*dAvC[1]+(4./3.)*dAvC[2]-(1./4.)*dAvC[3]); + cumulant[1] = (1./pow(fTuningR0[r],4.))*((-52./3.)*dAvC[0]+19.*dAvC[1]-(28./3.)*dAvC[2]+(11./6.)*dAvC[3]); + cumulant[2] = (3./pow(fTuningR0[r],6.))*(18.*dAvC[0]-24.*dAvC[1]+14.*dAvC[2]-3.*dAvC[3]); + cumulant[3] = (24./pow(fTuningR0[r],8.))*((-4.)*dAvC[0]+6.*dAvC[1]-4.*dAvC[2]+1.*dAvC[3]); + if(cumulant[0]>=0.) {flow[0] = pow(cumulant[0],1./2.);} + if(cumulant[1]<=0.) {flow[1] = pow(-1.*cumulant[1],1./4.);} + if(cumulant[2]>=0.) {flow[2] = pow((1./4.)*cumulant[2],1./6.);} + if(cumulant[3]<=0.) {flow[3] = pow((-1./33.)*cumulant[3],1./8.);} + } + else if(pMax==5) + { + cumulant[0] = (-1./(60*fTuningR0[r]*fTuningR0[r]))*((-300.)*dAvC[0]+300.*dAvC[1]-200.*dAvC[2]+75.*dAvC[3]-12.*dAvC[4]); + cumulant[1] = (-1./(6.*pow(fTuningR0[r],4.)))*(154.*dAvC[0]-214.*dAvC[1]+156.*dAvC[2]-61.*dAvC[3]+10.*dAvC[4]); + cumulant[2] = (3./(2.*pow(fTuningR0[r],6.)))*(71.*dAvC[0]-118.*dAvC[1]+98.*dAvC[2]-41.*dAvC[3]+7.*dAvC[4]); + cumulant[3] = (-24./pow(fTuningR0[r],8.))*(14.*dAvC[0]-26.*dAvC[1]+24.*dAvC[2]-11.*dAvC[3]+2.*dAvC[4]); + cumulant[4] = (120./pow(fTuningR0[r],10.))*(5.*dAvC[0]-10.*dAvC[1]+10.*dAvC[2]-5.*dAvC[3]+1.*dAvC[4]); + if(cumulant[0]>=0.) {flow[0] = pow(cumulant[0],1./2.);} + if(cumulant[1]<=0.) {flow[1] = pow(-1.*cumulant[1],1./4.);} + if(cumulant[2]>=0.) {flow[2] = pow((1./4.)*cumulant[2],1./6.);} + if(cumulant[3]<=0.) {flow[3] = pow((-1./33.)*cumulant[3],1./8.);} + if(cumulant[4]>=0.) {flow[4] = pow((1./456.)*cumulant[4],1./10.);} + } + else if(pMax==8) + { + cumulant[0] = (1./(fTuningR0[r]*fTuningR0[r]))*(8.*dAvC[0]-14.*dAvC[1]+(56./3.)*dAvC[2]-(35./2.)*dAvC[3] + + (56./5.)*dAvC[4]-(14./3.)*dAvC[5]+(8./7.)*dAvC[6]-(1./8.)*dAvC[7]); + cumulant[1] = (1./pow(fTuningR0[r],4.))*((-1924./35.)*dAvC[0]+(621./5.)*dAvC[1]-(8012./45.)*dAvC[2] + + (691./4.)*dAvC[3]-(564./5.)*dAvC[4]+(2143./45.)*dAvC[5]-(412./35.)*dAvC[6]+(363./280.)*dAvC[7]); + cumulant[2] = (1./pow(fTuningR0[r],6.))*(349.*dAvC[0]-(18353./20.)*dAvC[1]+(7173./5.)*dAvC[2] + - 1457.*dAvC[3]+(4891./5.)*dAvC[4]-(1683./4.)*dAvC[5]+(527./5.)*dAvC[6]-(469./40.)*dAvC[7]); + cumulant[3] = (1./pow(fTuningR0[r],8.))*((-10528./5.)*dAvC[0]+(30578./5.)*dAvC[1]-(51456./5.)*dAvC[2] + + 10993.*dAvC[3]-(38176./5.)*dAvC[4]+(16818./5.)*dAvC[5]-(4288./5.)*dAvC[6]+(967./10.)*dAvC[7]); + cumulant[4] = (1./pow(fTuningR0[r],10.))*(11500.*dAvC[0]-35800.*dAvC[1]+63900.*dAvC[2]-71600.*dAvC[3] + + 51620.*dAvC[4]-23400.*dAvC[5]+6100.*dAvC[6]-700.*dAvC[7]); + cumulant[5] = (1./pow(fTuningR0[r],12.))*(-52560.*dAvC[0]+172080.*dAvC[1]-321840.*dAvC[2]+376200.*dAvC[3] + - 281520.*dAvC[4]+131760.*dAvC[5]-35280.*dAvC[6]+4140.*dAvC[7]); + cumulant[6] = (1./pow(fTuningR0[r],14.))*(176400.*dAvC[0]-599760.*dAvC[1]+1164240.*dAvC[2]-1411200.*dAvC[3] + + 1093680.*dAvC[4]-529200.*dAvC[5]+146160.*dAvC[6]-17640.*dAvC[7]); + cumulant[7] = (1./pow(fTuningR0[r],16.))*(-322560*dAvC[0]+1128960.*dAvC[1]-2257920.*dAvC[2]+2822400.*dAvC[3] + - 2257920.*dAvC[4]+1128960.*dAvC[5]-322560.*dAvC[6]+40320.*dAvC[7]); + if(cumulant[0]>=0.) {flow[0] = pow(cumulant[0],1./2.);} + if(cumulant[1]<=0.) {flow[1] = pow(-1.*cumulant[1],1./4.);} + if(cumulant[2]>=0.) {flow[2] = pow((1./4.)*cumulant[2],1./6.);} + if(cumulant[3]<=0.) {flow[3] = pow((-1./33.)*cumulant[3],1./8.);} + if(cumulant[4]>=0.) {flow[4] = pow((1./456.)*cumulant[4],1./10.);} + if(cumulant[5]<=0.) {flow[5] = pow((-1./9460.)*cumulant[5],1./12.);} + if(cumulant[6]>=0.) {flow[6] = pow((1./274800.)*cumulant[6],1./14.);} + if(cumulant[7]<=0.) {flow[7] = pow((-1./10643745.)*cumulant[7],1./16.);} + } + // Store cumulants and reference flow: + for(Int_t co=0;coSetBinContent(co+1,cumulant[co]); + fTuningFlow[r][pq]->SetBinContent(co+1,flow[co]); + } + } // end of for(Int_t pq=0;pq<5;pq++) + } // end of for(Int_t r=0;r<10;r++) + +} // end of void AliFlowAnalysisWithCumulants::FinalizeTuning() +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::FillGeneratingFunctionsForDifferentTuningParameters(AliFlowEventSimple *anEvent) +{ + // Fill generating function for reference flow evaluated for different tuning parameters. + + Int_t pMax[5] = {2,3,4,5,8}; + Int_t qMax[5] = {5,7,9,11,17}; + + // Particle variables and weights: + Double_t dPhi = 0.; // azimuthal angle in the laboratory frame + Double_t dPt = 0.; // transverse momentum + Double_t dEta = 0.; // pseudorapidity Double_t wPhi = 1.; // phi weight Double_t wPt = 1.; // pt weight Double_t wEta = 1.; // eta weight - - if(useWeights) + + Int_t nPrim = anEvent->NumberOfTracks(); // nPrim = total number of primary tracks, i.e. nPrim = nRP + nPOI + rest, where: + // nRP = # of particles used to determine the reaction plane; + // nPOI = # of particles of interest for a detailed flow analysis; + // rest = # of particles which are not niether RPs nor POIs. + + Int_t nRP = anEvent->GetEventNSelTracksRP(); // nRP = # of particles used to determine the reaction plane; + + Double_t tuningGenFunEBE[10][5][8][17] = {{{{0.}}}}; + for(Int_t r=0;r<10;r++) { - if(!fWeightsList) - { - cout<<" WARNING: fWeightsList is NULL pointer in GFC::Make() !!!!"<(fWeightsList->FindObject("phi_weights")); - if(!phiWeights) - { - cout<<" WARNING: couldn't access the histogram with phi weights in GFC::Make() !!!!"<(fWeightsList->FindObject("pt_weights")); - if(!ptWeights) + for(Int_t p=0;pGetTrack(i); + if(aftsTrack && aftsTrack->InRPSelection()) { - etaWeights = dynamic_cast(fWeightsList->FindObject("eta_weights")); - if(!etaWeights) + // Access particle variables and weights: + dPhi = aftsTrack->Phi(); + dPt = aftsTrack->Pt(); + dEta = aftsTrack->Eta(); + if(fUsePhiWeights && fnBinsPhi) // determine phi weight for this particle: { - cout<<" WARNING: couldn't access the histogram with eta weights in GFC::Make() !!!!"<GetBinContent(1+(Int_t)(TMath::Floor(dPhi*fnBinsPhi/TMath::TwoPi()))); } - } - } - //--------------------------------------------------------------------------------------------------------- + if(fUsePtWeights && fnBinsPt) // determine pt weight for this particle: + { + wPt = fPtWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fPtBinWidth))); + } + if(fUseEtaWeights && fEtaBinWidth) // determine eta weight for this particle: + { + wEta = fEtaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fEtaBinWidth))); + } + // Fill the generating functions: + for(Int_t r=0;r<10;r++) // 10 different values for interpolating parameter r0 + { + if(TMath::Abs(fTuningR0[r])<1.e-10) continue; + for(Int_t pq=0;pq<5;pq++) // 5 different values for set (pMax,qMax) + { + for(Int_t p=0;pInRPSelection()) + } // end of for(Int_t i=0;i9) //generating function formalism applied here make sense only for selected multiplicity >= 10 - { - //fill the common control histograms - fCommonHists->FillControlHistograms(anEvent); + // Store G[p][q]: + for(Int_t r=0;r<10;r++) + { + for(Int_t pq=0;pq<5;pq++) + { + for(Int_t p=0;pFill((Double_t)p,(Double_t)q,tuningGenFunEBE[r][pq][p][q],1.);} + } + } + } + } - //initializing the generating function G[p][q] for integrated flow - Double_t genfunG[fgkPmax][fgkQmax]; +} // end of void AliFlowAnalysisWithCumulants::FillGeneratingFunctionsForDifferentTuningParameters(AliFlowEventSimple *anEvent) + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::FillGeneratingFunctionForReferenceFlow(AliFlowEventSimple *anEvent) +{ + // Fill generating function for reference flow for current event. + + // Particle variables and weights: + Double_t dPhi = 0.; // azimuthal angle in the laboratory frame + Double_t dPt = 0.; // transverse momentum + Double_t dEta = 0.; // pseudorapidity + Double_t wPhi = 1.; // phi weight + Double_t wPt = 1.; // pt weight + Double_t wEta = 1.; // eta weight + + Int_t nPrim = anEvent->NumberOfTracks(); // nPrim = total number of primary tracks, i.e. nPrim = nRP + nPOI + rest, where: + // nRP = # of particles used to determine the reaction plane; + // nPOI = # of particles of interest for a detailed flow analysis; + // rest = # of particles which are not niether RPs nor POIs. - for(Int_t p=0;pGetEventNSelTracksRP(); // nRP = # of particles used to determine the reaction plane; + + // Initializing the generating function G[p][q] for reference flow for current event: + Int_t pMax = fGEBE->GetNrows(); + Int_t qMax = fGEBE->GetNcols(); + for(Int_t p=0;pGetTrack(i); - if(fTrack && fTrack->InRPSelection()) + AliFlowTrackSimple *aftsTrack = anEvent->GetTrack(i); + if(aftsTrack && aftsTrack->InRPSelection()) { - nSelTracksRP++; - // get azimuthal angle, momentum and pseudorapidity of a particle: - dPhi = fTrack->Phi(); - dPt = fTrack->Pt(); - dEta = fTrack->Eta(); - // phi weights: - if(fUsePhiWeights) - { - nBinsPhi = phiWeights->GetNbinsX(); - if(nBinsPhi) - { - wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi()))); - } - } - // pt weights: - if(fUsePtWeights) - { - if(fBinWidthPt) - { - wPt = ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fBinWidthPt))); - } - } - // eta weights: - if(fUseEtaWeights) - { - if(fBinWidthEta) - { - wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fBinWidthEta))); - } + crossCheckRP++; + // Access particle variables and weights: + dPhi = aftsTrack->Phi(); + dPt = aftsTrack->Pt(); + dEta = aftsTrack->Eta(); + if(fUsePhiWeights && fnBinsPhi) // determine phi weight for this particle: + { + wPhi = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*fnBinsPhi/TMath::TwoPi()))); } - // evaluate the generating function: - for(Int_t p=0;pGetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fPtBinWidth))); + } + if(fUseEtaWeights && fEtaBinWidth) // determine eta weight for this particle: + { + wEta = fEtaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fEtaBinWidth))); + } + // Fill the generating function: + for(Int_t p=0;p - fAverageOfSquaredWeight->Fill(0.,pow(wPhi*wPt*wEta,2.),1); - } // end of if(fTrack && fTrack->InRPSelection()) + // Fill the profile to calculate <>: + fAverageOfSquaredWeight->Fill(0.5,pow(wPhi*wPt*wEta,2.),1.); + } // end of if(aftsTrack && aftsTrack->InRPSelection()) } // end of for(Int_t i=0;i - for(Int_t p=0;pGetEventNSelTracksRP())) + { + cout<: + // Determine first the event weight for G[p][q]: + // (to be improved - this can be implemented much better, this shall be executed only once out of Make(), eventWeight should be a data member) + Double_t eventWeight = 0.; + if(!strcmp(fMultiplicityWeight->Data(),"unit")) + { + eventWeight = 1.; + } else if(!strcmp(fMultiplicityWeight->Data(),"multiplicity")) + { + eventWeight = anEvent->GetEventNSelTracksRP(); + } + // Store G[p][q] weighted appropriately: + for(Int_t p=0;pFill((Double_t)p,(Double_t)q,genfunG[p][q],1); + fReferenceFlowGenFun->Fill((Double_t)p,(Double_t)q,(*fGEBE)(p,q),eventWeight); } } - //storing the selected multiplicity (if fAvMultIntFlow is not filled here then you had wrongly selected the particles used for integrated flow) - if(nSelTracksRP==nEventNSelTracksRP) +} // end of void AliFlowAnalysisWithCumulants::FillGeneratingFunctionForReferenceFlow(AliFlowEventSimple* anEvent) + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::FillQvectorComponents(AliFlowEventSimple* anEvent) +{ + // Fill components of Q-vector for current event (needed for error calculation). + + // Remark: Components are stored in profile fQvectorComponents whose binning is organized as follows: + // 1st bin: Q_x + // 2nd bin: Q_y + // 3rd bin: (Q_x)^2 + // 4th bin: (Q_y)^2 + + AliFlowVector afv; + afv.Set(0.,0.); + afv.SetMult(0); + + Int_t n = 2; // to be removed + + if(anEvent) { - fAvMultIntFlowGFC->Fill(0.,nSelTracksRP,1); - } else - { - cout<<"WARNING: nSelTracksRP != nEventNSelTracksRP in GFC::Make(). Something is wrong with RP flagging !!!!"<GetQ(1*n,fWeightsList,fUsePhiWeights,fUsePtWeights,fUseEtaWeights); // get the Q-vector for this event + fQvectorComponents->Fill(0.5,afv.X(),1.); // in the 1st bin fill Q_x + fQvectorComponents->Fill(1.5,afv.Y(),1.); // in the 2nd bin fill Q_y + fQvectorComponents->Fill(2.5,pow(afv.X(),2.),1.); // in the 3rd bin fill (Q_x)^2 + fQvectorComponents->Fill(3.5,pow(afv.Y(),2.),1.); // in the 4th bin fill (Q_y)^2 + } + +} // end of void AliFlowAnalysisWithCumulants::FillQvectorComponents(AliFlowEventSimple* anEvent) + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::FillGeneratingFunctionForDiffFlow(AliFlowEventSimple* anEvent) +{ + // Fill generating function for differential flow for the current event. + + // Remark 0: Generating function D[b][p][q] is a complex number => real and imaginary part are calculated separately + // (b denotes pt or eta bin); + // Remark 1: Note that bellow G[p][q] is needed, the value of generating function for reference flow for the CURRENT event. + // This values is obtained in method FillGeneratingFunctionForReferenceFlow() as TMatrixD fGEBE; + // Remark 2: Results for D[b][p][q] are stored in 3D profiles fDiffFlowGenFun[0=Re,1=Im][0=RP,1=POI][0=pt,1=eta] in order to + // automatically get and at the end of the day. + + // Particle variables and weights: + Double_t dPhi = 0.; // azimuthal angle in the laboratory frame + Double_t dPt = 0.; // transverse momentum + Double_t dEta = 0.; // pseudorapidity + Double_t wPhi = 1.; // phi weight + Double_t wPt = 1.; // pt weight + Double_t wEta = 1.; // eta weight + + // pMax and qMax: + Int_t pMax = fGEBE->GetNrows(); + Int_t qMax = fGEBE->GetNcols(); + + Int_t nPrim = anEvent->NumberOfTracks(); // nPrim = total number of primary tracks, i.e. nPrim = nRP + nPOI + rest, where: + // nRP = # of particles used to determine the reaction plane; + // nPOI = # of particles of interest for a detailed flow analysis; + // rest = # of particles which are not niether RPs nor POIs. - // calculating Q-vector of event (needed for errors) - AliFlowVector fQVector; - fQVector.Set(0.,0.); - fQVector.SetMult(0); - fQVector=anEvent->GetQ(1*n,fWeightsList,fUsePhiWeights,fUsePtWeights,fUseEtaWeights); // get the Q vector for this event - fQVectorComponentsGFC->Fill(0.,fQVector.X(),1); // in the 1st bin fill Q_x - fQVectorComponentsGFC->Fill(1.,fQVector.Y(),1); // in the 2nd bin fill Q_y - fQVectorComponentsGFC->Fill(2.,pow(fQVector.X(),2.),1); // in the 3rd bin fill (Q_x)^2 - fQVectorComponentsGFC->Fill(3.,pow(fQVector.Y(),2.),1); // in the 4th bin fill (Q_y)^2 - - //3D profiles for differential flow in pt and eta - //second loop over data: evaluating the generating function D[b][p][q] for differential flow - //remark 0: D[b][p][q] is a complex number => real and imaginary part are calculated separately - //remark 1: note that bellow G[p][q] is needed, the value of generating function for integrated flow for the CURRENT event - //remark 2: results are stored in 3D profiles in order to automatically get and + Int_t nRP = anEvent->GetEventNSelTracksRP(); // nRP = # of particles used to determine the reaction plane + + // Start the second loop over event in order to evaluate the generating function D[b][p][q] for differential flow: for(Int_t i=0;iGetTrack(i); - if(fTrack) + AliFlowTrackSimple *aftsTrack = anEvent->GetTrack(i); + if(aftsTrack) { - if(fTrack->InRPSelection() && fTrack->InPOISelection()) + if(!(aftsTrack->InRPSelection() || aftsTrack->InPOISelection())) continue; + // Differential flow of POIs: + if(aftsTrack->InPOISelection()) { - // get azimuthal angle, momentum and pseudorapidity of a particle: - dPhi = fTrack->Phi(); - dPt = fTrack->Pt(); - dEta = fTrack->Eta(); - - fPtBinPOINoOfParticles->Fill(dPt,dPt,1.); - fEtaBinPOINoOfParticles->Fill(dEta,dEta,1.); - - // phi weights: - if(fUsePhiWeights) - { - nBinsPhi = phiWeights->GetNbinsX(); - if(nBinsPhi) - { - wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi()))); - } - } - // pt weights: - if(fUsePtWeights) - { - if(fBinWidthPt) - { - wPt = ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fBinWidthPt))); - } - } - // eta weights: - if(fUseEtaWeights) - { - if(fBinWidthEta) - { - wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fBinWidthEta))); - } - } + // Get azimuthal angle, momentum and pseudorapidity of a particle: + dPhi = aftsTrack->Phi(); + dPt = aftsTrack->Pt(); + dEta = aftsTrack->Eta(); + Double_t ptEta[2] = {dPt,dEta}; - for(Int_t p=0;pFill(dPt/fBinWidthPt,(Double_t)p,(Double_t)q,genfunG[p][q]*cos(fgkMltpl*fgkFlow*dPhi)/(1.+wPhi*wPt*wEta*(2.*fR0*sqrt(p+1.)/nSelTracksRP)*cos(fgkFlow*dPhi-2.*q*TMath::Pi()/fgkQmax)),1.); - //imaginary part (Pt) - fDiffFlowPtPOIGenFunIm->Fill(dPt/fBinWidthPt,(Double_t)p,(Double_t)q,genfunG[p][q]*sin(fgkMltpl*fgkFlow*dPhi)/(1.+wPhi*wPt*wEta*(2.*fR0*sqrt(p+1.)/nSelTracksRP)*cos(fgkFlow*dPhi-2.*q*TMath::Pi()/fgkQmax)),1.); - //real part (Eta) - fDiffFlowEtaPOIGenFunRe->Fill(dEta/fBinWidthEta,(Double_t)p,(Double_t)q,genfunG[p][q]*cos(fgkMltpl*fgkFlow*dPhi)/(1.+wPhi*wPt*wEta*(2.*fR0*sqrt(p+1.)/nSelTracksRP)*cos(fgkFlow*dPhi-2.*q*TMath::Pi()/fgkQmax)),1.); - //imaginary part (Eta) - fDiffFlowEtaPOIGenFunIm->Fill(dEta/fBinWidthEta,(Double_t)p,(Double_t)q,genfunG[p][q]*sin(fgkMltpl*fgkFlow*dPhi)/(1.+wPhi*wPt*wEta*(2.*fR0*sqrt(p+1.)/nSelTracksRP)*cos(fgkFlow*dPhi-2.*q*TMath::Pi()/fgkQmax)),1.); - } + // Count number of POIs in pt/eta bin: + for(Int_t pe=0;pe<2;pe++) + { + fNoOfParticlesInBin[1][pe]->Fill(ptEta[pe],ptEta[pe],1.); } - } - else if(fTrack->InPOISelection() && !(fTrack->InRPSelection())) - { - // get azimuthal angle, momentum and pseudorapidity of a particle: - dPhi = fTrack->Phi(); - dPt = fTrack->Pt(); - dEta = fTrack->Eta(); - - fPtBinPOINoOfParticles->Fill(dPt,dPt,1.); - fEtaBinPOINoOfParticles->Fill(dEta,dEta,1.); - for(Int_t p=0;pInRPSelection())) // particle was flagged only as POI { - for(Int_t q=0;qFill(dPt/fBinWidthPt,(Double_t)p,(Double_t)q,genfunG[p][q]*cos(fgkMltpl*fgkFlow*dPhi),1.); - //imaginary part (Pt) - fDiffFlowPtPOIGenFunIm->Fill(dPt/fBinWidthPt,(Double_t)p,(Double_t)q,genfunG[p][q]*sin(fgkMltpl*fgkFlow*dPhi),1.); - //real part (Eta) - fDiffFlowEtaPOIGenFunRe->Fill(dEta/fBinWidthEta,(Double_t)p,(Double_t)q,genfunG[p][q]*cos(fgkMltpl*fgkFlow*dPhi),1.); - //imaginary part (Eta) - fDiffFlowEtaPOIGenFunIm->Fill(dEta/fBinWidthEta,(Double_t)p,(Double_t)q,genfunG[p][q]*sin(fgkMltpl*fgkFlow*dPhi),1.); - } - } - } - //RP: - if(fTrack->InRPSelection()) - { - // get azimuthal angle, momentum and pseudorapidity of a particle: - dPhi = fTrack->Phi(); - dPt = fTrack->Pt(); - dEta = fTrack->Eta(); - - fPtBinRPNoOfParticles->Fill(dPt,dPt,1.); - fEtaBinRPNoOfParticles->Fill(dEta,dEta,1.); - - // phi weights: - if(fUsePhiWeights) + for(Int_t q=0;qFill(ptEta[pe],(Double_t)p,(Double_t)q, // to be improved - hardwired weight 1. in the line bellow + (*fGEBE)(p,q)*cos(fMultiple*fHarmonic*dPhi),1.); + } + else if(ri==1) // Imaginary part (to be improved - this can be implemented better) + { + fDiffFlowGenFun[ri][1][pe]->Fill(ptEta[pe],(Double_t)p,(Double_t)q, // to be improved - hardwired weight 1. in the line bellow + (*fGEBE)(p,q)*sin(fMultiple*fHarmonic*dPhi),1.); + } + } // end of for(Int_t pe=0;pe<2;pe++) + } // end of for(Int_t ri=0;ri<2;ri++) + } // end of for(Int_t q=0;qInRPSelection())) // particle was flagged only as POI + else if(aftsTrack->InRPSelection()) // particle was flagged both as RP and POI { - nBinsPhi = phiWeights->GetNbinsX(); - if(nBinsPhi) + // If particle weights were used, get them: + if(fUsePhiWeights && fnBinsPhi) // determine phi weight for this particle: { - wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi()))); + wPhi = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*fnBinsPhi/TMath::TwoPi()))); } - } - // pt weights: - if(fUsePtWeights) - { - if(fBinWidthPt) + if(fUsePtWeights && fnBinsPt) // determine pt weight for this particle: { - wPt = ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fBinWidthPt))); - } - } - // eta weights: - if(fUseEtaWeights) - { - if(fBinWidthEta) + wPt = fPtWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fPtBinWidth))); + } + if(fUseEtaWeights && fEtaBinWidth) // determine eta weight for this particle: { - wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fBinWidthEta))); - } + wEta = fEtaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fEtaBinWidth))); + } + // Fill generating function: + for(Int_t p=0;pFill(ptEta[pe],(Double_t)p,(Double_t)q, // to be improved - hardwired weight 1. in the line bellow + (*fGEBE)(p,q)*cos(fMultiple*fHarmonic*dPhi)/(1.+wPhi*wPt*wEta*(2.*fR0*sqrt(p+1.)/nRP)*cos(fHarmonic*dPhi-2.*q*TMath::Pi()/qMax)),1.); + } + else if(ri==1) // Imaginary part (to be improved - this can be implemented better) + { + fDiffFlowGenFun[ri][1][pe]->Fill(ptEta[pe],(Double_t)p,(Double_t)q, // to be improved - hardwired weight 1. in the line bellow + (*fGEBE)(p,q)*sin(fMultiple*fHarmonic*dPhi)/(1.+wPhi*wPt*wEta*(2.*fR0*sqrt(p+1.)/nRP)*cos(fHarmonic*dPhi-2.*q*TMath::Pi()/qMax)),1.); + } + } // end of for(Int_t pe=0;pe<2;pe++) + } // end of for(Int_t ri=0;ri<2;ri++) + } // end of for(Int_t q=0;qInRPSelection()) // particle was flagged both as RP and POI + } // end of if(aftsTrack->InPOISelection()) + // Differential flow of RPs: + if(aftsTrack->InRPSelection()) + { + // Get azimuthal angle, momentum and pseudorapidity of a particle: + dPhi = aftsTrack->Phi(); + dPt = aftsTrack->Pt(); + dEta = aftsTrack->Eta(); + Double_t ptEta[2] = {dPt,dEta}; + + // Count number of RPs in pt/eta bin: + for(Int_t pe=0;pe<2;pe++) + { + fNoOfParticlesInBin[0][pe]->Fill(ptEta[pe],ptEta[pe],1.); } - for(Int_t p=0;pFill(dPt/fBinWidthPt,(Double_t)p,(Double_t)q,genfunG[p][q]*cos(fgkMltpl*fgkFlow*dPhi)/(1.+wPhi*wPt*wEta*(2.*fR0*sqrt(p+1.)/nSelTracksRP)*cos(fgkFlow*dPhi-2.*q*TMath::Pi()/fgkQmax)),1.); - //imaginary part (Pt) - fDiffFlowPtRPGenFunIm->Fill(dPt/fBinWidthPt,(Double_t)p,(Double_t)q,genfunG[p][q]*sin(fgkMltpl*fgkFlow*dPhi)/(1.+wPhi*wPt*wEta*(2.*fR0*sqrt(p+1.)/nSelTracksRP)*cos(fgkFlow*dPhi-2.*q*TMath::Pi()/fgkQmax)),1.); - //real part (Eta) - fDiffFlowEtaRPGenFunRe->Fill(dEta/fBinWidthEta,(Double_t)p,(Double_t)q,genfunG[p][q]*cos(fgkMltpl*fgkFlow*dPhi)/(1.+wPhi*wPt*wEta*(2.*fR0*sqrt(p+1.)/nSelTracksRP)*cos(fgkFlow*dPhi-2.*q*TMath::Pi()/fgkQmax)),1.); - //imaginary part (Eta) - fDiffFlowEtaRPGenFunIm->Fill(dEta/fBinWidthEta,(Double_t)p,(Double_t)q,genfunG[p][q]*sin(fgkMltpl*fgkFlow*dPhi)/(1.+wPhi*wPt*wEta*(2.*fR0*sqrt(p+1.)/nSelTracksRP)*cos(fgkFlow*dPhi-2.*q*TMath::Pi()/fgkQmax)),1.); - } + wPhi = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*fnBinsPhi/TMath::TwoPi()))); } - }//end of if(fTrack->InRPSelection()) - }//end of if(fTrack) - }//ending the second loop over data - - - - /* - //sixteen 2D profiles for differential flow - for(Int_t i=0;iGetTrack(i); - if (fTrack && fTrack->InPOISelection()){ - //for(Int_t p=0;pFill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[0][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(0.+1.)/nSelTracksRP)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.); - //imaginary part - fDiffFlowGenFunIm0->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[0][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(0.+1.)/nSelTracksRP)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.); - //----------------------------------------------------------------------- - //real part - fDiffFlowGenFunRe1->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[1][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(1.+1.)/nSelTracksRP)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.); - //imaginary part - fDiffFlowGenFunIm1->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[1][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(1.+1.)/nSelTracksRP)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.); - //----------------------------------------------------------------------- - //real part - fDiffFlowGenFunRe2->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[2][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(2.+1.)/nSelTracksRP)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.); - //imaginary part - fDiffFlowGenFunIm2->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[2][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(2.+1.)/nSelTracksRP)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.); - //----------------------------------------------------------------------- - //real part - fDiffFlowGenFunRe3->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[3][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(3.+1.)/nSelTracksRP)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.); - //imaginary part - fDiffFlowGenFunIm3->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[3][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(3.+1.)/nSelTracksRP)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.); - //----------------------------------------------------------------------- - //real part - fDiffFlowGenFunRe4->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[4][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(4.+1.)/nSelTracksRP)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.); - //imaginary part - fDiffFlowGenFunIm4->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[4][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(4.+1.)/nSelTracksRP)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.); - //----------------------------------------------------------------------- - //real part - fDiffFlowGenFunRe5->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[5][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(5.+1.)/nSelTracksRP)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.); - //imaginary part - fDiffFlowGenFunIm5->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[5][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(5.+1.)/nSelTracksRP)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.); - //----------------------------------------------------------------------- - //real part - fDiffFlowGenFunRe6->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[6][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(6.+1.)/nSelTracksRP)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.); - //imaginary part - fDiffFlowGenFunIm6->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[6][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(6.+1.)/nSelTracksRP)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.); - //----------------------------------------------------------------------- - //real part - fDiffFlowGenFunRe7->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[7][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(7.+1.)/nSelTracksRP)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.); - //imaginary part - fDiffFlowGenFunIm7->Fill(fTrack->Pt()/fBinWidth,(Double_t)q,genfunG[7][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/(1.+(2.*fR0*sqrt(7.+1.)/nSelTracksRP)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax)),1.); - } - //} - } - }//ending the second loop over data - */ - - }//end of if(nEventNSelTracksRP>9) - - - - - - - - - - - - - - //off the record: numerical equations for cumulants solved up to different highest order - if(fOtherEquations) - { - //running over data - Int_t nPrimOE = anEvent->NumberOfTracks();//total multiplicity - - Int_t nEventNSelTracksRPOE = anEvent->GetEventNSelTracksRP(); - - Double_t genfunG4[fgkPmax4][fgkQmax4]; - Double_t genfunG6[fgkPmax6][fgkQmax6]; - Double_t genfunG8[fgkPmax8][fgkQmax8]; - Double_t genfunG16[fgkPmax16][fgkQmax16]; - for(Int_t p=0;p15) fAvMultIntFlow16GFC->Fill(0.,nEventNSelTracksRPOE,1); - if(nEventNSelTracksRPOE>7) fAvMultIntFlow8GFC->Fill(0.,nEventNSelTracksRPOE,1); - if(nEventNSelTracksRPOE>5) fAvMultIntFlow6GFC->Fill(0.,nEventNSelTracksRPOE,1); - if(nEventNSelTracksRPOE>3) fAvMultIntFlow4GFC->Fill(0.,nEventNSelTracksRPOE,1); - - //first loop over data: evaluating the generating function G[p][q] for integrated flow - for(Int_t i=0;iGetTrack(i); - if(fTrack && fTrack->InRPSelection()) - { - for(Int_t p=0;pGetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fPtBinWidth))); + } + if(fUseEtaWeights && fEtaBinWidth) // determine eta weight for this particle: + { + wEta = fEtaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fEtaBinWidth))); + } + // Fill generating function: + for(Int_t p=0;p15) + for(Int_t ri=0;ri<2;ri++) { - genfunG16[p][q]*=(1.+(2.*fR0*sqrt(p+1.)/nEventNSelTracksRPOE)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax16)); - } - if(p7) - { - genfunG8[p][q]*=(1.+(2.*fR0*sqrt(p+1.)/nEventNSelTracksRPOE)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax8)); - } - if(p5) + if(ri==0) // Real part (to be improved - this can be implemented better) { - genfunG6[p][q]*=(1.+(2.*fR0*sqrt(p+1.)/nEventNSelTracksRPOE)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax6)); - } - if(pFill(ptEta[pe],(Double_t)p,(Double_t)q, // to be improved - hardwired weight 1. in the line bellow + (*fGEBE)(p,q)*cos(fMultiple*fHarmonic*dPhi)/(1.+wPhi*wPt*wEta*(2.*fR0*sqrt(p+1.)/nRP)*cos(fHarmonic*dPhi-2.*q*TMath::Pi()/qMax)),1.); + } + else if(ri==1) // Imaginary part (to be improved - this can be implemented better) { - if(nEventNSelTracksRPOE>3) - { - genfunG4[p][q]*=(1.+(2.*fR0*sqrt(p+1.)/nEventNSelTracksRPOE)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax4)); - } + fDiffFlowGenFun[ri][0][pe]->Fill(ptEta[pe],(Double_t)p,(Double_t)q, // to be improved - hardwired weight 1. in the line bellow + (*fGEBE)(p,q)*sin(fMultiple*fHarmonic*dPhi)/(1.+wPhi*wPt*wEta*(2.*fR0*sqrt(p+1.)/nRP)*cos(fHarmonic*dPhi-2.*q*TMath::Pi()/qMax)),1.); } - } + } // end of for(Int_t pe=0;pe<2;pe++) + } // end of for(Int_t ri=0;ri<2;ri++) + } // end of for(Int_t q=0;qInRPSelection()) + } // end of if(aftsTrack) + } // end of for(Int_t i=0;iSetHistList(outputListHistos); + if(!fHistList) + { + cout<GetPointersForBaseHistograms(); + this->GetPointersForCommonControlHistograms(); + this->GetPointersForCommonResultsHistograms(); + this->GetPointersForParticleWeightsHistograms(); + this->GetPointersForReferenceFlowObjects(); + this->GetPointersForDiffFlowObjects(); + } else + { + cout<(fHistList->FindObject(analysisSettingsName.Data())); + if(analysisSettings) + { + this->SetAnalysisSettings(analysisSettings); + } else + { + cout<(fHistList->FindObject(commonHistsName.Data())); + if(commonHist) + { + this->SetCommonHists(commonHist); + } else + { + cout< + (fHistList->FindObject(commonHistResults2ndOrderName.Data())); + if(commonHistRes2nd) + { + this->SetCommonHistsResults2nd(commonHistRes2nd); + } else + { + cout< + (fHistList->FindObject(commonHistResults4thOrderName.Data())); + if(commonHistRes4th) + { + this->SetCommonHistsResults4th(commonHistRes4th); + } else + { + cout< + (fHistList->FindObject(commonHistResults6thOrderName.Data())); + if(commonHistRes6th) + { + this->SetCommonHistsResults6th(commonHistRes6th); + } else + { + cout< + (fHistList->FindObject(commonHistResults8thOrderName.Data())); + if(commonHistRes8th) + { + this->SetCommonHistsResults8th(commonHistRes8th); + } else + { + cout<(fHistList->FindObject("Weights")); + if(weightsList) this->SetWeightsList(weightsList); + TString fUseParticleWeightsName = "fUseParticleWeights"; + TProfile *useParticleWeights = dynamic_cast(weightsList->FindObject(fUseParticleWeightsName.Data())); + if(useParticleWeights) + { + this->SetUseParticleWeights(useParticleWeights); + fUsePhiWeights = (Int_t)fUseParticleWeights->GetBinContent(1); + fUsePtWeights = (Int_t)fUseParticleWeights->GetBinContent(2); + fUseEtaWeights = (Int_t)fUseParticleWeights->GetBinContent(3); + } + +} // end of void AliFlowAnalysisWithCumulants::GetPointersForParticleWeightsHistograms(); + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::GetPointersForReferenceFlowObjects() +{ + // Get pointers for all objects relevant for calculation of reference flow. + + // a) Get pointers to all lists relevant for reference flow; + // b) Get pointer to profile holding flags; + // c) Get pointers to all objects in the list fReferenceFlowProfiles; + // d) Get pointers to all objects in the list fReferenceFlowResults. + + // a) Get pointers to all lists relevant for reference flow: + TList *referenceFlowList = dynamic_cast(fHistList->FindObject("Reference Flow")); + if(!referenceFlowList) + { + cout<(referenceFlowList->FindObject("Profiles")); + if(!referenceFlowProfiles) + { + cout<(referenceFlowList->FindObject("Results")); + if(!referenceFlowResults) + { + cout<(referenceFlowList->FindObject(referenceFlowFlagsName.Data())); + if(referenceFlowFlags) + { + this->SetReferenceFlowFlags(referenceFlowFlags); + } else + { + cout<(referenceFlowProfiles->FindObject(referenceFlowGenFunName.Data())); + if(referenceFlowGenFun) + { + this->SetReferenceFlowGenFun(referenceFlowGenFun); + } else + { + cout<(referenceFlowProfiles->FindObject(qvectorComponentsName.Data())); + if(qvectorComponents) + { + this->SetQvectorComponents(qvectorComponents); + } else + { + cout<>, where w = wPhi*wPt*wEta: + TString averageOfSquaredWeightName = "fAverageOfSquaredWeight"; + TProfile *averageOfSquaredWeight = dynamic_cast(referenceFlowProfiles->FindObject(averageOfSquaredWeightName.Data())); + if(averageOfSquaredWeight) + { + this->SetAverageOfSquaredWeight(averageOfSquaredWeight); + } else + { + cout<(referenceFlowResults->FindObject(referenceFlowCumulantsName.Data())); + if(referenceFlowCumulants) + { + this->SetReferenceFlowCumulants(referenceFlowCumulants); + } else + { + cout<(referenceFlowResults->FindObject(referenceFlowName.Data())); + if(referenceFlow) + { + this->SetReferenceFlow(referenceFlow); + } else + { + cout<(referenceFlowResults->FindObject(chiName.Data())); + if(chi) + { + this->SetChi(chi); + } else + { + cout<(fHistList->FindObject("Differential Flow")); // to be improved (hardwired name) + if(!diffFlowList) + { + cout<(diffFlowList->FindObject("Profiles")); // to be improved (hardwired name) + if(!diffFlowProfiles) + { + cout<(diffFlowList->FindObject("Results")); // to be improved (hardwired name) + if(!diffFlowResults) + { + cout<(diffFlowList->FindObject(diffFlowFlagsName.Data())); + if(diffFlowFlags) + { + this->SetDiffFlowFlags(diffFlowFlags); + } else + { + cout< // to be improved - harwired name fDiffFlowGenFun in the line bellow + (diffFlowProfiles->FindObject(Form("fDiffFlowGenFun (%s, %s, %s)",reIm[ri].Data(),rpPoi[rp].Data(),ptEta[pe].Data()))); + if(diffFlowGenFun[ri][rp][pe]) + { + this->SetDiffFlowGenFun(diffFlowGenFun[ri][rp][pe],ri,rp,pe); + } else + { + cout< // to be improved - harwired name fNoOfParticlesInBin in the line bellow + (diffFlowProfiles->FindObject(Form("fNoOfParticlesInBin (%s, %s)",rpPoi[rp].Data(),ptEta[pe].Data()))); + if(noOfParticlesInBin[rp][pe]) + { + this->SetNoOfParticlesInBin(noOfParticlesInBin[rp][pe],rp,pe); + } else + { + cout< // to be improved - hardwired name fDiffFlowCumulants in the line bellow + (diffFlowResults->FindObject(Form("fDiffFlowCumulants (%s, %s, %s)",rpPoi[rp].Data(),ptEta[pe].Data(),order[co].Data()))); + if(diffFlowCumulants[rp][pe][co]) + { + this->SetDiffFlowCumulants(diffFlowCumulants[rp][pe][co],rp,pe,co); + } else + { + cout<InRPSelection()) - }//ending the loop over data + } + } + } + // Differential flow per pt/eta bin for RPs/POIs: + TH1D *diffFlow[2][2][4] = {{{NULL}}}; + for(Int_t rp=0;rp<2;rp++) + { + for(Int_t pe=0;pe<2;pe++) + { + for(Int_t co=0;co<4;co++) + { + diffFlow[rp][pe][co] = dynamic_cast // to be improved - hardwired name fDiffFlow in the line bellow + (diffFlowResults->FindObject(Form("fDiffFlow (%s, %s, %s)",rpPoi[rp].Data(),ptEta[pe].Data(),order[co].Data()))); + if(diffFlow[rp][pe][co]) + { + this->SetDiffFlow(diffFlow[rp][pe][co],rp,pe,co); + } else + { + cout<GetHistPtPOI())->Clone(); + } else if(rpPoi == "RP") + { + yieldPt = (TH1F*)(fCommonHists->GetHistPtRP())->Clone(); + } + + TH1D *flow2ndPt = (TH1D*)fDiffFlow[rp][0][0]->Clone(); + TH1D *flow4thPt = (TH1D*)fDiffFlow[rp][0][1]->Clone(); + TH1D *flow6thPt = (TH1D*)fDiffFlow[rp][0][2]->Clone(); + TH1D *flow8thPt = (TH1D*)fDiffFlow[rp][0][3]->Clone(); + Double_t dvn2nd = 0., dvn4th = 0., dvn6th = 0., dvn8th = 0.; // differential flow + Double_t dErrvn2nd = 0., dErrvn4th = 0., dErrvn6th = 0., dErrvn8th = 0.; // error on differential flow + Double_t dVn2nd = 0., dVn4th = 0., dVn6th = 0., dVn8th = 0.; // integrated flow + Double_t dErrVn2nd = 0., dErrVn4th = 0., dErrVn6th = 0., dErrVn8th = 0.; // error on integrated flow + Double_t dYield = 0.; // pt yield + Double_t dSum2nd = 0., dSum4th = 0., dSum6th = 0., dSum8th = 0.; // needed for normalizing integrated flow + fnBinsPt = flow2ndPt->GetXaxis()->GetNbins(); + // looping over pt bins: + for(Int_t p=1;p<=fnBinsPt;p++) + { + dvn2nd = flow2ndPt->GetBinContent(p); + dvn4th = flow4thPt->GetBinContent(p); + dvn6th = flow6thPt->GetBinContent(p); + dvn8th = flow8thPt->GetBinContent(p); + + dErrvn2nd = flow2ndPt->GetBinError(p); + dErrvn4th = flow4thPt->GetBinError(p); + dErrvn6th = flow6thPt->GetBinError(p); + dErrvn8th = flow8thPt->GetBinError(p); + + dYield = yieldPt->GetBinContent(p); + + dVn2nd += dvn2nd*dYield; + dVn4th += dvn4th*dYield; + dVn6th += dvn6th*dYield; + dVn8th += dvn8th*dYield; + + dSum2nd += dYield; + dSum4th += dYield; + dSum6th += dYield; + dSum8th += dYield; + + dErrVn2nd += dYield*dYield*dErrvn2nd*dErrvn2nd; + dErrVn4th += dYield*dYield*dErrvn4th*dErrvn4th; + dErrVn6th += dYield*dYield*dErrvn6th*dErrvn6th; + dErrVn8th += dYield*dYield*dErrvn8th*dErrvn8th; + } // end of for(Int_t p=1;p<=fnBinsPt;p++) + + // normalizing the results for integrated flow: + if(dSum2nd) + { + dVn2nd /= dSum2nd; + dErrVn2nd /= (dSum2nd*dSum2nd); + dErrVn2nd = TMath::Sqrt(dErrVn2nd); + } + if(dSum4th) + { + dVn4th /= dSum4th; + dErrVn4th /= (dSum4th*dSum4th); + dErrVn4th = TMath::Sqrt(dErrVn4th); + } + if(dSum6th) + { + dVn6th /= dSum6th; + dErrVn6th /= (dSum6th*dSum6th); + dErrVn6th = TMath::Sqrt(dErrVn6th); + } + if(dSum8th) + { + dVn8th /= dSum8th; + dErrVn8th /= (dSum8th*dSum8th); + dErrVn8th = TMath::Sqrt(dErrVn8th); + } + + // storing the results for integrated flow in common hist results: + if(rpPoi == "POI") + { + fCommonHistsResults2nd->FillIntegratedFlowPOI(dVn2nd,dErrVn2nd); + fCommonHistsResults4th->FillIntegratedFlowPOI(dVn4th,dErrVn4th); + fCommonHistsResults6th->FillIntegratedFlowPOI(dVn6th,dErrVn6th); + fCommonHistsResults8th->FillIntegratedFlowPOI(dVn8th,dErrVn8th); + } + else if(rpPoi == "RP") + { + fCommonHistsResults2nd->FillIntegratedFlowRP(dVn2nd,dErrVn2nd); + fCommonHistsResults4th->FillIntegratedFlowRP(dVn4th,dErrVn4th); + fCommonHistsResults6th->FillIntegratedFlowRP(dVn6th,dErrVn6th); + fCommonHistsResults8th->FillIntegratedFlowRP(dVn8th,dErrVn8th); + } + + delete flow2ndPt; + delete flow4thPt; + delete flow6thPt; + delete flow8thPt; + delete yieldPt; + +} // end of void AliFlowAnalysisWithCumulants::CalculateIntegratedFlow(TString rpPoi) + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::FillCommonHistResultsForDifferentialFlow(TString rpPoi) +{ + // Fill common result histograms for differential flow. + // (to be improved - this method can be implemented much better) + + Int_t rp = -1; + + if(rpPoi == "RP") + { + rp = 0; + } else if(rpPoi == "POI") + { + rp = 1; + } + + // pt: + for(Int_t p=1;p<=fnBinsPt;p++) + { + // Result: + Double_t v2 = fDiffFlow[rp][0][0]->GetBinContent(p); + Double_t v4 = fDiffFlow[rp][0][1]->GetBinContent(p); + Double_t v6 = fDiffFlow[rp][0][2]->GetBinContent(p); + Double_t v8 = fDiffFlow[rp][0][3]->GetBinContent(p); + // Error: + Double_t v2Error = fDiffFlow[rp][0][0]->GetBinError(p); + Double_t v4Error = fDiffFlow[rp][0][1]->GetBinError(p); + Double_t v6Error = fDiffFlow[rp][0][2]->GetBinError(p); + Double_t v8Error = fDiffFlow[rp][0][3]->GetBinError(p); + // Fill common hist results: + if(rpPoi == "RP") + { + fCommonHistsResults2nd->FillDifferentialFlowPtRP(p,v2,v2Error); + fCommonHistsResults4th->FillDifferentialFlowPtRP(p,v4,v4Error); + fCommonHistsResults6th->FillDifferentialFlowPtRP(p,v6,v6Error); + fCommonHistsResults8th->FillDifferentialFlowPtRP(p,v8,v8Error); + } else if(rpPoi == "POI") + { + fCommonHistsResults2nd->FillDifferentialFlowPtPOI(p,v2,v2Error); + fCommonHistsResults4th->FillDifferentialFlowPtPOI(p,v4,v4Error); + fCommonHistsResults6th->FillDifferentialFlowPtPOI(p,v6,v6Error); + fCommonHistsResults8th->FillDifferentialFlowPtPOI(p,v8,v8Error); + } + } // end of for(Int_t p=1;p<=fnBinsPt;p++) + + // eta: + for(Int_t e=1;e<=fnBinsEta;e++) + { + // Results: + Double_t v2 = fDiffFlow[rp][1][0]->GetBinContent(e); + Double_t v4 = fDiffFlow[rp][1][1]->GetBinContent(e); + Double_t v6 = fDiffFlow[rp][1][2]->GetBinContent(e); + Double_t v8 = fDiffFlow[rp][1][3]->GetBinContent(e); + // Errors: + Double_t v2Error = fDiffFlow[rp][1][0]->GetBinError(e); + Double_t v4Error = fDiffFlow[rp][1][1]->GetBinError(e); + Double_t v6Error = fDiffFlow[rp][1][2]->GetBinError(e); + Double_t v8Error = fDiffFlow[rp][1][3]->GetBinError(e); + // Fill common hist results: + if(rpPoi == "RP") + { + fCommonHistsResults2nd->FillDifferentialFlowEtaRP(e,v2,v2Error); + fCommonHistsResults4th->FillDifferentialFlowEtaRP(e,v4,v4Error); + fCommonHistsResults6th->FillDifferentialFlowEtaRP(e,v6,v6Error); + fCommonHistsResults8th->FillDifferentialFlowEtaRP(e,v8,v8Error); + } else if(rpPoi == "POI") + { + fCommonHistsResults2nd->FillDifferentialFlowEtaPOI(e,v2,v2Error); + fCommonHistsResults4th->FillDifferentialFlowEtaPOI(e,v4,v4Error); + fCommonHistsResults6th->FillDifferentialFlowEtaPOI(e,v6,v6Error); + fCommonHistsResults8th->FillDifferentialFlowEtaPOI(e,v8,v8Error); + } + } // end of for(Int_t e=1;e<=fnBinsEta;e++) + +} // end of void AliFlowAnalysisWithCumulants::FillCommonHistResultsForDifferentialFlow(TString rpPoi) + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::CalculateDifferentialFlow(TString rpPoi, TString ptEta) +{ + // Calculate differential flow for RPs/POIs vs pt/eta from cumulants. + + Int_t rp = -1; // RP or POI + Int_t pe = -1; // pt or eta + + if(rpPoi == "RP") + { + rp = 0; + } else if(rpPoi == "POI") + { + rp = 1; + } + if(ptEta == "pt") + { + pe = 0; + } else if(ptEta == "eta") + { + pe = 1; + } + + // Reference cumulants: + Double_t gfc2 = fReferenceFlowCumulants->GetBinContent(1); // reference 2nd order cumulant + Double_t gfc4 = fReferenceFlowCumulants->GetBinContent(2); // reference 4th order cumulant + Double_t gfc6 = fReferenceFlowCumulants->GetBinContent(3); // reference 6th order cumulant + Double_t gfc8 = fReferenceFlowCumulants->GetBinContent(4); // reference 8th order cumulant + + Int_t nBins = fDiffFlowCumulants[rp][pe][0]->GetXaxis()->GetNbins(); + + for(Int_t b=1;b<=nBins;b++) + { + // Differential cumulants: + Double_t gfd2 = fDiffFlowCumulants[rp][pe][0]->GetBinContent(b); // differential 2nd order cumulant + Double_t gfd4 = fDiffFlowCumulants[rp][pe][1]->GetBinContent(b); // differential 4th order cumulant + Double_t gfd6 = fDiffFlowCumulants[rp][pe][2]->GetBinContent(b); // differential 6th order cumulant + Double_t gfd8 = fDiffFlowCumulants[rp][pe][3]->GetBinContent(b); // differential 8th order cumulant + // Differential flow: + Double_t v2 = 0.; // v'{2,GFC} + Double_t v4 = 0.; // v'{4,GFC} + Double_t v6 = 0.; // v'{6,GFC} + Double_t v8 = 0.; // v'{8,GFC} + // 2nd order: + if(gfc2>0.) + { + v2 = gfd2/pow(gfc2,0.5); + fDiffFlow[rp][pe][0]->SetBinContent(b,v2); + } + // 4th order: + if(gfc4<0.) + { + v4 = -gfd4/pow(-gfc4,.75); + fDiffFlow[rp][pe][1]->SetBinContent(b,v4); + } + // 6th order: + if(gfc6>0.) + { + v6 = gfd6/(4.*pow((1./4.)*gfc6,(5./6.))); + fDiffFlow[rp][pe][2]->SetBinContent(b,v6); + } + // 8th order: + if(gfc8<0.) + { + v8 = -gfd8/(33.*pow(-(1./33.)*gfc8,(7./8.))); + fDiffFlow[rp][pe][3]->SetBinContent(b,v8); + } + } // end of for(Int_t b=1;b<=nBins;b++) + +} // end of void AliFlowAnalysisWithCumulants::CalculateDifferentialFlow(TString rpPoi,TString ptEta) + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::CalculateDifferentialFlowErrors(TString rpPoi,TString ptEta) +{ + // Calculate errors of differential flow. + + Int_t rp = -1; // RP or POI + Int_t pe = -1; // pt or eta + + if(rpPoi == "RP") + { + rp = 0; + } else if(rpPoi == "POI") + { + rp = 1; + } + if(ptEta == "pt") + { + pe = 0; + } else if(ptEta == "eta") + { + pe = 1; + } + + // Resolution chi: + Double_t chi2 = fChi->GetBinContent(1); + Double_t chi4 = fChi->GetBinContent(2); + //Double_t chi6 = fChi->GetBinContent(3); + //Double_t chi8 = fChi->GetBinContent(4); + + Int_t nBins = fNoOfParticlesInBin[rp][pe]->GetXaxis()->GetNbins(); + for(Int_t b=1;b<=nBins;b++) + { + Int_t nParticles = (Int_t)fNoOfParticlesInBin[rp][pe]->GetBinEntries(b); + // Error of 2nd order estimate: + if(chi2>0. && nParticles>0.) + { + Double_t v2Error = pow((1./(2.*nParticles))*((1.+pow(chi2,2.))/pow(chi2,2.)),0.5); + fDiffFlow[rp][pe][0]->SetBinError(b,v2Error); + } + // Error of 4th order estimate: + if(chi4>0. && nParticles>0.) + { + Double_t v4Error = pow((1./(2.*nParticles))*((2.+6.*pow(chi4,2.)+pow(chi4,4.)+pow(chi4,6.))/pow(chi4,6.)),0.5); + fDiffFlow[rp][pe][1]->SetBinError(b,v4Error); + } + // Error of 6th order estimate: + //if(chi6>0. && nParticles>0.) + //{ + // Double_t v6Error = ... // to be improved - yet to be calculated + fDiffFlow[rp][pe][2]->SetBinError(b,0.); + //} + // Error of 8th order estimate: + //if(chi8>0. && nParticles>0.) + //{ + // Double_t v8Error = ... // to be improved - yet to be calculated + fDiffFlow[rp][pe][3]->SetBinError(b,0.); + //} + } // end of for(Int_t b=1;b<=nBins;b++) + +} // end of void AliFlowAnalysisWithCumulants::CalculateDifferentialFlowErrors(TString rpPoi,TString ptEta) + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::CalculateCumulantsForDiffFlow(TString rpPoi,TString ptEta) +{ + // Calculate cumulants for differential flow. + + Int_t rp = -1; // RP or POI + Int_t pe = -1; // pt or eta + + if(rpPoi == "RP") + { + rp = 0; + } else if(rpPoi == "POI") + { + rp = 1; + } + if(ptEta == "pt") + { + pe = 0; + } else if(ptEta == "eta") + { + pe = 1; + } + + // [nBins][pMax][qMax]: + Int_t nBins = fDiffFlowGenFun[0][rp][pe]->GetXaxis()->GetNbins(); + Int_t pMax = fDiffFlowGenFun[0][rp][pe]->GetYaxis()->GetNbins(); + Int_t qMax = fDiffFlowGenFun[0][rp][pe]->GetZaxis()->GetNbins(); + // + TMatrixD dAvG(pMax,qMax); + dAvG.Zero(); + for(Int_t p=0;pGetBinContent(fReferenceFlowGenFun->GetBin(p+1,q+1)); + } + } + // Loop over pt/eta bins and calculate differential cumulants: + for(Int_t b=0;bGetBinEntries(b+1); + for(Int_t p=0;p1.e-44) + { + Double_t X = fDiffFlowGenFun[0][rp][pe]->GetBinContent(fDiffFlowGenFun[0][rp][pe]->GetBin(b+1,p+1,q+1))/dAvG(p,q); // see Ollitrault's Practical guide (Eq. 11) + Double_t Y = fDiffFlowGenFun[1][rp][pe]->GetBinContent(fDiffFlowGenFun[0][rp][pe]->GetBin(b+1,p+1,q+1))/dAvG(p,q); // see Ollitrault's Practical guide (Eq. 11) + tempSum += cos(fMultiple*2.*q*TMath::Pi()/qMax)*X + + sin(fMultiple*2.*q*TMath::Pi()/qMax)*Y; + } + } + D[p] = (pow(fR0*pow(p+1.0,0.5),fMultiple)/qMax)*tempSum; + } + gfc[0] = (1./(fR0*fR0))*(5.*D[0]-5.*D[1]+(10./3.)*D[2]-(5./4.)*D[3]+(1./5.)*D[4]); + gfc[1] = (1./pow(fR0,4.))*((-77./6.)*D[0]+(107./6.)*D[1]-(13./1.)*D[2]+(61./12.)*D[3]-(5./6.)*D[4]); + gfc[2] = (1./pow(fR0,6.))*((71./2.)*D[0]-59.*D[1]+49.*D[2]-(41./2.)*D[3]+(7./2.)*D[4]); + gfc[3] = (1./pow(fR0,8.))*(-84.*D[0]+156.*D[1]-144.*D[2]+66.*D[3]-12.*D[4]); + // gfc[4] = (1./pow(fR0,10.))*(120.*D[0]-240.*D[1]+240.*D[2]-120.*D[3]+24.*D[4]); // 10th order cumulant (to be improved - where to store it?) + // Store cumulants: + for(Int_t co=0;co<4;co++) + { + fDiffFlowCumulants[rp][pe][co]->SetBinContent(b+1,gfc[co]); + } + } + +} // end of void AliFlowAnalysisWithCumulants::CalculateCumulantsForDiffFlow(TString rpPoi, TString ptEta) + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::PrintFinalResults(TString type) +{ + // Printing on the screen the final results for reference flow and for integrated flow of RPs and POIs. + + Int_t n = fHarmonic; + + Double_t dVn[4] = {0.}; // array to hold Vn{2}, Vn{4}, Vn{6} and Vn{8} + Double_t dVnErr[4] = {0.}; // array to hold errors of Vn{2}, Vn{4}, Vn{6} and Vn{8} + + if(type == "RF") + { + dVn[0] = (fCommonHistsResults2nd->GetHistIntFlow())->GetBinContent(1); + dVnErr[0] = (fCommonHistsResults2nd->GetHistIntFlow())->GetBinError(1); + dVn[1] = (fCommonHistsResults4th->GetHistIntFlow())->GetBinContent(1); + dVnErr[1] = (fCommonHistsResults4th->GetHistIntFlow())->GetBinError(1); + dVn[2] = (fCommonHistsResults6th->GetHistIntFlow())->GetBinContent(1); + dVnErr[2] = (fCommonHistsResults6th->GetHistIntFlow())->GetBinError(1); + dVn[3] = (fCommonHistsResults8th->GetHistIntFlow())->GetBinContent(1); + dVnErr[3] = (fCommonHistsResults8th->GetHistIntFlow())->GetBinError(1); + } else if(type == "RP") + { + dVn[0] = (fCommonHistsResults2nd->GetHistIntFlowRP())->GetBinContent(1); + dVnErr[0] = (fCommonHistsResults2nd->GetHistIntFlowRP())->GetBinError(1); + dVn[1] = (fCommonHistsResults4th->GetHistIntFlowRP())->GetBinContent(1); + dVnErr[1] = (fCommonHistsResults4th->GetHistIntFlowRP())->GetBinError(1); + dVn[2] = (fCommonHistsResults6th->GetHistIntFlowRP())->GetBinContent(1); + dVnErr[2] = (fCommonHistsResults6th->GetHistIntFlowRP())->GetBinError(1); + dVn[3] = (fCommonHistsResults8th->GetHistIntFlowRP())->GetBinContent(1); + dVnErr[3] = (fCommonHistsResults8th->GetHistIntFlowRP())->GetBinError(1); + } else if(type == "POI") + { + dVn[0] = (fCommonHistsResults2nd->GetHistIntFlowPOI())->GetBinContent(1); + dVnErr[0] = (fCommonHistsResults2nd->GetHistIntFlowPOI())->GetBinError(1); + dVn[1] = (fCommonHistsResults4th->GetHistIntFlowPOI())->GetBinContent(1); + dVnErr[1] = (fCommonHistsResults4th->GetHistIntFlowPOI())->GetBinError(1); + dVn[2] = (fCommonHistsResults6th->GetHistIntFlowPOI())->GetBinContent(1); + dVnErr[2] = (fCommonHistsResults6th->GetHistIntFlowPOI())->GetBinError(1); + dVn[3] = (fCommonHistsResults8th->GetHistIntFlowPOI())->GetBinContent(1); + dVnErr[3] = (fCommonHistsResults8th->GetHistIntFlowPOI())->GetBinError(1); + } else + { + cout< = "<<(Double_t)fCommonHists->GetHistMultRP()->GetMean()< = "<<(Double_t)fCommonHists->GetHistMultRP()->GetMean()< = "<<(Double_t)fCommonHists->GetHistMultPOI()->GetMean()<GetBinContent(1); + Double_t v4 = fReferenceFlow->GetBinContent(2); + Double_t v6 = fReferenceFlow->GetBinContent(3); + Double_t v8 = fReferenceFlow->GetBinContent(4); + // Errors: + Double_t v2Error = fReferenceFlow->GetBinError(1); + Double_t v4Error = fReferenceFlow->GetBinError(2); + Double_t v6Error = fReferenceFlow->GetBinError(3); + Double_t v8Error = fReferenceFlow->GetBinError(4); + // Fill results end errors in common hist results: + fCommonHistsResults2nd->FillIntegratedFlow(v2,v2Error); + fCommonHistsResults4th->FillIntegratedFlow(v4,v4Error); + fCommonHistsResults6th->FillIntegratedFlow(v6,v6Error); + fCommonHistsResults8th->FillIntegratedFlow(v8,v8Error); + // Chi: + Double_t chi2 = fChi->GetBinContent(1); + Double_t chi4 = fChi->GetBinContent(2); + Double_t chi6 = fChi->GetBinContent(3); + Double_t chi8 = fChi->GetBinContent(4); + // Fill resolution chi in common hist results: + fCommonHistsResults2nd->FillChi(chi2); + fCommonHistsResults4th->FillChi(chi4); + fCommonHistsResults6th->FillChi(chi6); + fCommonHistsResults8th->FillChi(chi8); + +} // end of AliFlowAnalysisWithCumulants::FillCommonHistResultsForReferenceFlow() + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::CalculateReferenceFlowError() +{ + // Calculate error of reference flow harmonics. + + // Generating Function Cumulants: + Double_t gfc2 = fReferenceFlowCumulants->GetBinContent(1); // GFC{2} + Double_t gfc4 = fReferenceFlowCumulants->GetBinContent(2); // GFC{4} + Double_t gfc6 = fReferenceFlowCumulants->GetBinContent(3); // GFC{6} + Double_t gfc8 = fReferenceFlowCumulants->GetBinContent(4); // GFC{8} + // Reference flow estimates: + Double_t v2 = fReferenceFlow->GetBinContent(1); // v{2,GFC} + Double_t v4 = fReferenceFlow->GetBinContent(2); // v{4,GFC} + Double_t v6 = fReferenceFlow->GetBinContent(3); // v{6,GFC} + Double_t v8 = fReferenceFlow->GetBinContent(4); // v{8,GFC} + // Statistical errors of reference flow estimates: + Double_t v2Error = 0.; // statistical error of v{2,GFC} + Double_t v4Error = 0.; // statistical error of v{4,GFC} + Double_t v6Error = 0.; // statistical error of v{6,GFC} + Double_t v8Error = 0.; // statistical error of v{8,GFC} + // Chi: + Double_t chi2 = 0.; + Double_t chi4 = 0.; + Double_t chi6 = 0.; + Double_t chi8 = 0.; + // : + Double_t dAvQx = fQvectorComponents->GetBinContent(1); // + Double_t dAvQy = fQvectorComponents->GetBinContent(2); // + Double_t dAvQ2x = fQvectorComponents->GetBinContent(3); // <(Q_x)^2> + Double_t dAvQ2y = fQvectorComponents->GetBinContent(4); // <(Q_y)^2> + // : + Double_t dAvw2 = 1.; + if(fnEvts>0) + { + dAvw2 = fAverageOfSquaredWeight->GetBinContent(1); + if(TMath::Abs(dAvw2)<1.e-44) + { + cout<0. && fAvM>0. && dAvw2>0. && gfc2>=0.) + { + if((dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(gfc2,(1./2.))*(fAvM/dAvw2),2.)>0.)) + { + chi2 = (fAvM*v2)/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(v2*fAvM/dAvw2,2.),0.5); + } + if(TMath::Abs(chi2)>1.e-44) + { + v2Error = pow(((1./(2.*fAvM*fnEvts))*((1.+2.*pow(chi2,2))/(2.*pow(chi2,2)))),0.5); + } + } + // Calculating statistical error of v{4,GFC}: + if(fnEvts>0 && fAvM>0 && dAvw2>0 && gfc4<=0.) + { + if((dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(-gfc4,(1./4.))*(fAvM/dAvw2),2.)>0.)) + { + chi4 = (fAvM*v4)/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(v4*fAvM/dAvw2,2.),0.5); + } + if(TMath::Abs(chi4)>1.e-44) + { + v4Error = (1./(pow(2.*fAvM*fnEvts,0.5)))*pow((1.+4.*pow(chi4,2)+1.*pow(chi4,4.)+2.*pow(chi4,6.))/(2.*pow(chi4,6.)),0.5); + } + } + // Calculating statistical error of v{6,GFC}: + if(fnEvts>0 && fAvM>0 && dAvw2>0 && gfc6>=0.) + { + if((dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow((1./4.)*gfc6,(1./6.))*(fAvM/dAvw2),2.)>0.)) + { + chi6 = (fAvM*v6)/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(v6*fAvM/dAvw2,2.),0.5); + } + if(TMath::Abs(chi6)>1.e-44) + { + v6Error = (1./(pow(2.*fAvM*fnEvts,0.5)))*pow((3.+18.*pow(chi6,2)+9.*pow(chi6,4.)+28.*pow(chi6,6.) + +12.*pow(chi6,8.)+24.*pow(chi6,10.))/(24.*pow(chi6,10.)),0.5); + } + } + // Calculating statistical error of v{8,GFC}: + if(fnEvts>0 && fAvM>0 && dAvw2>0 && gfc8<=0.) + { + if((dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(-(1./33.)*gfc8,(1./8.))*(fAvM/dAvw2),2.)>0.)) + { + chi8=(fAvM*v8)/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(v8*fAvM/dAvw2,2.),0.5); + } + if(TMath::Abs(chi8)>1.e-44) + { + v8Error = (1./(pow(2.*fAvM*fnEvts,0.5)))*pow((12.+96.*pow(chi8,2.)+72.*pow(chi8,4.)+304.*pow(chi8,6.) + +257.*pow(chi8,8.)+804.*pow(chi8,10.)+363.*pow(chi8,12.)+726.*pow(chi8,14.))/(726.*pow(chi8,14.)),0.5); + } + } + + // Store errors for reference flow: + fReferenceFlow->SetBinError(1,v2Error); + fReferenceFlow->SetBinError(2,v4Error); + fReferenceFlow->SetBinError(3,v6Error); + fReferenceFlow->SetBinError(4,v8Error); + // Store resolution chi: + fChi->SetBinContent(1,chi2); + fChi->SetBinContent(2,chi4); + fChi->SetBinContent(3,chi6); + fChi->SetBinContent(4,chi8); + +} // end of void AliFlowAnalysisWithCumulants::CalculateReferenceFlowError() + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::CalculateReferenceFlow() +{ + // Calculate from isotropic cumulants reference flow. + + // Generating Function Cumulants: + Double_t gfc2 = fReferenceFlowCumulants->GetBinContent(1); // GFC{2} + Double_t gfc4 = fReferenceFlowCumulants->GetBinContent(2); // GFC{4} + Double_t gfc6 = fReferenceFlowCumulants->GetBinContent(3); // GFC{6} + Double_t gfc8 = fReferenceFlowCumulants->GetBinContent(4); // GFC{8} + // Double_t gfc10 = fReferenceFlowCumulants->GetBinContent(5); // GFC{10} overflow + // Reference flow estimates: + Double_t v2 = 0.; // v{2,GFC} + Double_t v4 = 0.; // v{4,GFC} + Double_t v6 = 0.; // v{6,GFC} + Double_t v8 = 0.; // v{8,GFC} + // Double_t v10 = 0.; // v{10,GFC} overflow + // Calculate reference flow estimates from Q-cumulants: + if(gfc2>=0.) v2 = pow(gfc2,1./2.); + if(gfc4<=0.) v4 = pow(-1.*gfc4,1./4.); + if(gfc6>=0.) v6 = pow((1./4.)*gfc6,1./6.); + if(gfc8<=0.) v8 = pow((-1./33.)*gfc8,1./8.); + // if(gfc10>=0.) v10 = pow((1./456.)*gfc10,1./10.); + // Store results for reference flow: + fReferenceFlow->SetBinContent(1,v2); + fReferenceFlow->SetBinContent(2,v4); + fReferenceFlow->SetBinContent(3,v6); + fReferenceFlow->SetBinContent(4,v8); + // fReferenceFlow->SetBinContent(5,v10); // overflow + +} // end of void AliFlowAnalysisWithCumulants::CalculateReferenceFlow() + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::CalculateCumulantsForReferenceFlow() +{ + // Calculate cumulants for reference flow. + + Int_t pMax = fReferenceFlowGenFun->GetXaxis()->GetNbins(); + Int_t qMax = fReferenceFlowGenFun->GetYaxis()->GetNbins(); + + // + TMatrixD dAvG(pMax,qMax); + dAvG.Zero(); + Bool_t someAvGEntryIsNegative = kFALSE; + for(Int_t p=0;pGetBinContent(fReferenceFlowGenFun->GetBin(p+1,q+1)); + if(dAvG(p,q)<0.) + { + someAvGEntryIsNegative = kTRUE; + cout< is negative !!!! GFC results are meaningless.",p,q)<0. && !someAvGEntryIsNegative) + { + for(Int_t p=0;p (Remark: here <> stands for average over azimuth): + TVectorD dAvC(pMax); + dAvC.Zero(); + for(Int_t p=0;pSetBinContent(co+1,cumulant[co]); + } + +} // end of void AliFlowAnalysisWithCumulants::CalculateCumulantsForReferenceFlow() + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::GetAvMultAndNoOfEvts() +{ + // From relevant common control histogram get average multiplicity of RPs and number of events. + + fAvM = (Double_t)fCommonHists->GetHistMultRP()->GetMean(); + fnEvts = (Int_t)fCommonHists->GetHistMultRP()->GetEntries(); + +} // end of void AliFlowAnalysisWithCumulants::GetAvMultAndNoOfEvts() + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::InitializeArrays() +{ + // Initialize all arrays. + + for(Int_t ri=0;ri<2;ri++) + { + for(Int_t rp=0;rp<2;rp++) + { + for(Int_t pe=0;pe<2;pe++) + { + fDiffFlowGenFun[ri][rp][pe] = NULL; + } + } + } + for(Int_t rp=0;rp<2;rp++) + { + for(Int_t pe=0;pe<2;pe++) + { + fNoOfParticlesInBin[rp][pe] = NULL; + } + } + for(Int_t rp=0;rp<2;rp++) + { + for(Int_t pe=0;pe<2;pe++) + { + for(Int_t co=0;co<4;co++) + { + fDiffFlowCumulants[rp][pe][co] = NULL; + fDiffFlow[rp][pe][co] = NULL; + } + } + } + for(Int_t i=0;i<3;i++) + { + fPrintFinalResults[i] = kTRUE; + } + for(Int_t r=0;r<10;r++) + { + fTuningR0[r] = 0.; + for(Int_t pq=0;pq<5;pq++) + { + fTuningGenFun[r][pq] = NULL; + fTuningCumulants[r][pq] = NULL; + fTuningFlow[r][pq] = NULL; + } + } + +} // end of void AliFlowAnalysisWithCumulants::InitializeArrays() + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::CrossCheckSettings() +{ + // Cross-check the user settings before starting. + + // a) Cross check if the choice for multiplicity weight make sense. + + // a) Cross check if the choice for multiplicity weight make sense: + if(strcmp(fMultiplicityWeight->Data(),"unit") && + strcmp(fMultiplicityWeight->Data(),"multiplicity")) + { + cout<Data()<<"\"."< - for(Int_t p=0;pGetNbinsPhi(); + fPhiMin = AliFlowCommonConstants::GetMaster()->GetPhiMin(); + fPhiMax = AliFlowCommonConstants::GetMaster()->GetPhiMax(); + if(fnBinsPhi) fPhiBinWidth = (fPhiMax-fPhiMin)/fnBinsPhi; + fnBinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt(); + fPtMin = AliFlowCommonConstants::GetMaster()->GetPtMin(); + fPtMax = AliFlowCommonConstants::GetMaster()->GetPtMax(); + if(fnBinsPt) fPtBinWidth = (fPtMax-fPtMin)/fnBinsPt; + fnBinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta(); + fEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin(); + fEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax(); + if(fnBinsEta) fEtaBinWidth = (fEtaMax-fEtaMin)/fnBinsEta; + +} // end of void AliFlowAnalysisWithCumulants::AccessConstants() + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::BookAndFillWeightsHistograms() +{ + // book and fill histograms which hold phi, pt and eta weights + + if(!fWeightsList) + { + cout<<"WARNING (GFC): fWeightsList is NULL in AFAWGFC::BAFWH() !!!!"<SetLabelSize(0.06); + (fUseParticleWeights->GetXaxis())->SetBinLabel(1,"w_{#phi}"); + (fUseParticleWeights->GetXaxis())->SetBinLabel(2,"w_{p_{T}}"); + (fUseParticleWeights->GetXaxis())->SetBinLabel(3,"w_{#eta}"); + fUseParticleWeights->Fill(0.5,(Int_t)fUsePhiWeights); + fUseParticleWeights->Fill(1.5,(Int_t)fUsePtWeights); + fUseParticleWeights->Fill(2.5,(Int_t)fUseEtaWeights); + fWeightsList->Add(fUseParticleWeights); + if(fUsePhiWeights) { - for(Int_t q=0;qFindObject("phi_weights")) { - if(nEventNSelTracksRPOE>15) fIntFlowGenFun16->Fill((Double_t)p,(Double_t)q,genfunG16[p][q],1); - if(p(fWeightsList->FindObject("phi_weights")); + if(TMath::Abs(fPhiWeights->GetBinWidth(1)-fPhiBinWidth)>pow(10.,-6.)) { - if(nEventNSelTracksRPOE>7) fIntFlowGenFun8->Fill((Double_t)p,(Double_t)q,genfunG8[p][q],1); - if(p5) fIntFlowGenFun6->Fill((Double_t)p,(Double_t)q,genfunG6[p][q],1); - if(p3) fIntFlowGenFun4->Fill((Double_t)p,(Double_t)q,genfunG4[p][q],1); - } + cout<FindObject(\"phi_weights\") is NULL in AFAWGFC::BAFWH() !!!!"<FindObject("pt_weights")) + { + fPtWeights = dynamic_cast(fWeightsList->FindObject("pt_weights")); + if(TMath::Abs(fPtWeights->GetBinWidth(1)-fPtBinWidth)>pow(10.,-6.)) + { + cout<FindObject(\"pt_weights\") is NULL in AFAWGFC::BAFWH() !!!!"<FindObject("eta_weights")) + { + fEtaWeights = dynamic_cast(fWeightsList->FindObject("eta_weights")); + if(TMath::Abs(fEtaWeights->GetBinWidth(1)-fEtaBinWidth)>pow(10.,-6.)) + { + cout<FindObject(\"eta_weights\") is NULL in AFAWGFC::BAFWH() !!!!"<SetTickLength(-0.01,"Y"); + fReferenceFlowFlags->SetMarkerStyle(25); + fReferenceFlowFlags->SetLabelSize(0.05); + fReferenceFlowFlags->SetLabelOffset(0.02,"Y"); + fReferenceFlowFlags->GetXaxis()->SetBinLabel(1,"Particle weights"); + fReferenceFlowFlags->GetXaxis()->SetBinLabel(2,"Event weights"); + fReferenceFlowList->Add(fReferenceFlowFlags); + + // c) Book all event-by-event quantities: + fGEBE = new TMatrixD(pMax,qMax); + + // d) Book all profiles: + // Average of the generating function for reference flow : + fReferenceFlowGenFun = new TProfile2D("fReferenceFlowGenFun","#LTG[p][q]#GT",pMax,0.,(Double_t)pMax,qMax,0.,(Double_t)qMax); + fReferenceFlowGenFun->SetXTitle("p"); + fReferenceFlowGenFun->SetYTitle("q"); + fReferenceFlowProfiles->Add(fReferenceFlowGenFun); + // Averages of Q-vector components: + fQvectorComponents = new TProfile("fQvectorComponents","Averages of Q-vector components",4,0.,4.); + fQvectorComponents->SetLabelSize(0.06); + fQvectorComponents->SetMarkerStyle(25); + fQvectorComponents->GetXaxis()->SetBinLabel(1,"#LTQ_{x}#GT"); // Q_{x} + fQvectorComponents->GetXaxis()->SetBinLabel(2,"#LTQ_{y}#GT"); // Q_{y} + fQvectorComponents->GetXaxis()->SetBinLabel(3,"#LTQ_{x}^{2}#GT"); // Q_{x}^{2} + fQvectorComponents->GetXaxis()->SetBinLabel(4,"#LTQ_{y}^{2}#GT"); // Q_{y}^{2} + fReferenceFlowProfiles->Add(fQvectorComponents); + // <>, where w = wPhi*wPt*wEta: + fAverageOfSquaredWeight = new TProfile("fAverageOfSquaredWeight","#LT#LTw^{2}#GT#GT",1,0,1); + fAverageOfSquaredWeight->SetLabelSize(0.06); + fAverageOfSquaredWeight->SetMarkerStyle(25); + fAverageOfSquaredWeight->SetLabelOffset(0.01); + fAverageOfSquaredWeight->GetXaxis()->SetBinLabel(1,"#LT#LTw^{2}#GT#GT"); + fReferenceFlowProfiles->Add(fAverageOfSquaredWeight); + + // e) Book all histograms: + // Final results for isotropic cumulants for reference flow: + TString referenceFlowCumulantsName = "fReferenceFlowCumulants"; + fReferenceFlowCumulants = new TH1D(referenceFlowCumulantsName.Data(),"Isotropic Generating Function Cumulants for reference flow",4,0,4); // to be improved (hw 4) + fReferenceFlowCumulants->SetLabelSize(0.05); + fReferenceFlowCumulants->SetMarkerStyle(25); + fReferenceFlowCumulants->GetXaxis()->SetBinLabel(1,"GFC{2}"); + fReferenceFlowCumulants->GetXaxis()->SetBinLabel(2,"GFC{4}"); + fReferenceFlowCumulants->GetXaxis()->SetBinLabel(3,"GFC{6}"); + fReferenceFlowCumulants->GetXaxis()->SetBinLabel(4,"GFC{8}"); + fReferenceFlowResults->Add(fReferenceFlowCumulants); + // Final results for reference flow: + fReferenceFlow = new TH1D("fReferenceFlow","Reference flow",4,0,4); // to be improved (hardwired 4) + fReferenceFlow->SetLabelSize(0.06); + fReferenceFlow->SetMarkerStyle(25); + fReferenceFlow->GetXaxis()->SetBinLabel(1,"v_{n}{2,GFC}"); + fReferenceFlow->GetXaxis()->SetBinLabel(2,"v_{n}{4,GFC}"); + fReferenceFlow->GetXaxis()->SetBinLabel(3,"v_{n}{6,GFC}"); + fReferenceFlow->GetXaxis()->SetBinLabel(4,"v_{n}{8,GFC}"); + fReferenceFlowResults->Add(fReferenceFlow); + // Final results for resolution: + fChi = new TH1D("fChi","Resolution",4,0,4); // to be improved (hardwired 4) + fChi->SetLabelSize(0.06); + fChi->SetMarkerStyle(25); + fChi->GetXaxis()->SetBinLabel(1,"#chi_{2}"); + fChi->GetXaxis()->SetBinLabel(2,"#chi_{4}"); + fChi->GetXaxis()->SetBinLabel(3,"#chi_{6}"); + fChi->GetXaxis()->SetBinLabel(4,"#chi_{8}"); + fReferenceFlowResults->Add(fChi); + +} // end of void AliFlowAnalysisWithCumulants::BookEverythingForReferenceFlow() + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::BookEverythingForTuning() +{ + // Book all objects relevant for tuning. + + // a) Define pMax's and qMax's: + // b) Book profile to hold all tuning parameters and flags; + // c) Book all profiles; + // d) Book all histograms. + + // a) Define pMax's and qMax's: + Int_t pMax[5] = {2,3,4,5,8}; + Int_t qMax[5] = {5,7,9,11,17}; + + // b) Book profile to hold all tuning parameters and flags: + TString tuningFlagsName = "fTuningFlags"; + fTuningFlags = new TProfile(tuningFlagsName.Data(),"Tuning parameters",10,0,10); + // fTuningFlags->SetTickLength(-0.01,"Y"); + fTuningFlags->SetMarkerStyle(25); + fTuningFlags->SetLabelSize(0.05); + fTuningFlags->SetLabelOffset(0.02,"X"); + for(Int_t r=1;r<=10;r++) + { + fTuningFlags->GetXaxis()->SetBinLabel(r,Form("r_{0,%d}",r-1)); + fTuningFlags->Fill(r-0.5,fTuningR0[r-1],1.); + } + fTuningList->Add(fTuningFlags); + + // c) Book all profiles: + // Average of the generating function for reference flow for different tuning parameters: + for(Int_t r=0;r<10;r++) + { + for(Int_t pq=0;pq<5;pq++) + { + fTuningGenFun[r][pq] = new TProfile2D(Form("fTuningGenFun (r_{0,%i}, pq set %i)",r,pq), + Form("#LTG[p][q]#GT for r_{0} = %f, p_{max} = %i, q_{max} = %i",fTuningR0[r],pMax[pq],qMax[pq]), + pMax[pq],0.,(Double_t)pMax[pq],qMax[pq],0.,(Double_t)qMax[pq]); + fTuningGenFun[r][pq]->SetXTitle("p"); + fTuningGenFun[r][pq]->SetYTitle("q"); + fTuningProfiles->Add(fTuningGenFun[r][pq]); } } -}//end of if(fOtherEquations) - -}//end of Make() + + // d) Book all histograms: + // Final results for isotropic cumulants for reference flow for different tuning parameters: + for(Int_t r=0;r<10;r++) + { + for(Int_t pq=0;pq<5;pq++) + { + fTuningCumulants[r][pq] = new TH1D(Form("fTuningCumulants (r_{0,%i}, pq set %i)",r,pq), + Form("GFC for r_{0} = %f, p_{max} = %i, q_{max} = %i",fTuningR0[r],pMax[pq],qMax[pq]), + pMax[pq],0,pMax[pq]); + // fTuningCumulants[r][pq]->SetLabelSize(0.05); + fTuningCumulants[r][pq]->SetMarkerStyle(25); + for(Int_t b=1;b<=pMax[pq];b++) + { + fTuningCumulants[r][pq]->GetXaxis()->SetBinLabel(b,Form("GFC{%i}",2*b)); + } + fTuningResults->Add(fTuningCumulants[r][pq]); + } + } + // Final results for reference flow for different tuning parameters: + for(Int_t r=0;r<10;r++) + { + for(Int_t pq=0;pq<5;pq++) + { + fTuningFlow[r][pq] = new TH1D(Form("fTuningFlow (r_{0,%i}, pq set %i)",r,pq), + Form("Reference flow for r_{0} = %f, p_{max} = %i, q_{max} = %i",fTuningR0[r],pMax[pq],qMax[pq]), + pMax[pq],0,pMax[pq]); + // fTuningFlow[r][pq]->SetLabelSize(0.06); + fTuningFlow[r][pq]->SetMarkerStyle(25); + for(Int_t b=1;b<=pMax[pq];b++) + { + fTuningFlow[r][pq]->GetXaxis()->SetBinLabel(b,Form("v{%i,GFC}",2*b)); + } + fTuningResults->Add(fTuningFlow[r][pq]); + } + } + +} // end of void AliFlowAnalysisWithCumulants::BookEverythingForTuning() //================================================================================================================ -void AliFlowAnalysisWithCumulants::GetOutputHistograms(TList *outputListHistos) +void AliFlowAnalysisWithCumulants::BookEverythingForDiffFlow() { - // get the pointers to all output histograms before calling Finish() - if (outputListHistos) - { - //Get the common histograms from the output list - //histograms to store the final results - TH1D *intFlowResults = dynamic_cast(outputListHistos->FindObject("fIntFlowResultsGFC")); - TH1D *diffFlowResults2 = dynamic_cast(outputListHistos->FindObject("fDiffFlowResults2ndOrderGFC")); - TH1D *diffFlowResults4 = dynamic_cast(outputListHistos->FindObject("fDiffFlowResults4thOrderGFC")); - TH1D *diffFlowResults6 = dynamic_cast(outputListHistos->FindObject("fDiffFlowResults6thOrderGFC")); - TH1D *diffFlowResults8 = dynamic_cast(outputListHistos->FindObject("fDiffFlowResults8thOrderGFC")); - - //common histograms to store the final results the integrated and differential flow - AliFlowCommonHistResults *commonHistRes2nd = dynamic_cast(outputListHistos->FindObject("AliFlowCommonHistResults2ndOrderGFC")); - AliFlowCommonHistResults *commonHistRes4th = dynamic_cast(outputListHistos->FindObject("AliFlowCommonHistResults4thOrderGFC")); - AliFlowCommonHistResults *commonHistRes6th = dynamic_cast(outputListHistos->FindObject("AliFlowCommonHistResults6thOrderGFC")); - AliFlowCommonHistResults *commonHistRes8th = dynamic_cast(outputListHistos->FindObject("AliFlowCommonHistResults8thOrderGFC")); - - //common control histogram - AliFlowCommonHist *commonHists = dynamic_cast(outputListHistos->FindObject("AliFlowCommonHistGFC")); - - //profiles with average values of generating functions for int. and diff. flow - TProfile2D *intFlowGenFun = dynamic_cast(outputListHistos->FindObject("fIntFlowGenFun")); - - TProfile2D *intFlowGenFun4 = dynamic_cast(outputListHistos->FindObject("fIntFlowGenFun4")); //only for other system of Eq. - TProfile2D *intFlowGenFun6 = dynamic_cast(outputListHistos->FindObject("fIntFlowGenFun6")); //only for other system of Eq. - TProfile2D *intFlowGenFun8 = dynamic_cast(outputListHistos->FindObject("fIntFlowGenFun8")); //only for other system of Eq. - TProfile2D *intFlowGenFun16 = dynamic_cast(outputListHistos->FindObject("fIntFlowGenFun16")); //only for other system of Eq. - - //RP, Pt: - TProfile3D *diffFlowPtRPGenFunRe = dynamic_cast(outputListHistos->FindObject("fDiffFlowPtRPGenFunRe")); - TProfile3D *diffFlowPtRPGenFunIm = dynamic_cast(outputListHistos->FindObject("fDiffFlowPtRPGenFunIm")); - TProfile *ptBinRPNoOfParticles = dynamic_cast(outputListHistos->FindObject("fPtBinRPNoOfParticles")); - - //RP, Eta: - TProfile3D *diffFlowEtaRPGenFunRe = dynamic_cast(outputListHistos->FindObject("fDiffFlowEtaRPGenFunRe")); - TProfile3D *diffFlowEtaRPGenFunIm = dynamic_cast(outputListHistos->FindObject("fDiffFlowEtaRPGenFunIm")); - TProfile *etaBinRPNoOfParticles = dynamic_cast(outputListHistos->FindObject("fEtaBinRPNoOfParticles")); - - //POI, Pt: - TProfile3D *diffFlowPtPOIGenFunRe = dynamic_cast(outputListHistos->FindObject("fDiffFlowPtPOIGenFunRe")); - TProfile3D *diffFlowPtPOIGenFunIm = dynamic_cast(outputListHistos->FindObject("fDiffFlowPtPOIGenFunIm")); - TProfile *ptBinPOINoOfParticles = dynamic_cast(outputListHistos->FindObject("fPtBinPOINoOfParticles")); - - //POI, Eta: - TProfile3D *diffFlowEtaPOIGenFunRe = dynamic_cast(outputListHistos->FindObject("fDiffFlowEtaPOIGenFunRe")); - TProfile3D *diffFlowEtaPOIGenFunIm = dynamic_cast(outputListHistos->FindObject("fDiffFlowEtaPOIGenFunIm")); - TProfile *etaBinPOINoOfParticles = dynamic_cast(outputListHistos->FindObject("fEtaBinPOINoOfParticles")); - - //average selected multiplicity (for int. flow) - TProfile *avMult = dynamic_cast(outputListHistos->FindObject("fAvMultIntFlowGFC")); - - TProfile *avMult4 = dynamic_cast(outputListHistos->FindObject("fAvMultIntFlow4GFCA")); //only for other system of Eq. - TProfile *avMult6 = dynamic_cast(outputListHistos->FindObject("fAvMultIntFlow6GFCA")); //only for other system of Eq. - TProfile *avMult8 = dynamic_cast(outputListHistos->FindObject("fAvMultIntFlow8GFCA")); //only for other system of Eq. - TProfile *avMult16 = dynamic_cast(outputListHistos->FindObject("fAvMultIntFlow16GFCA")); //only for other system of Eq. - - //average values of Q-vector components (1st bin: , 2nd bin: , 3rd bin: <(Q_x)^2>, 4th bin: <(Q_y)^2>) - TProfile *qVectorComponents = dynamic_cast(outputListHistos->FindObject("fQVectorComponentsGFC")); - - // - TProfile *averageOfSquaredWeight = dynamic_cast(outputListHistos->FindObject("fAverageOfSquaredWeight")); - - /* - TProfile2D *diffFlowPtGenFunRe0 = dynamic_cast(outputListHistos->FindObject("fdiffFlowPtGenFunRe0")); - TProfile2D *diffFlowPtGenFunRe1 = dynamic_cast(outputListHistos->FindObject("fdiffFlowPtGenFunRe1")); - TProfile2D *diffFlowPtGenFunRe2 = dynamic_cast(outputListHistos->FindObject("fdiffFlowPtGenFunRe2")); - TProfile2D *diffFlowPtGenFunRe3 = dynamic_cast(outputListHistos->FindObject("fdiffFlowPtGenFunRe3")); - TProfile2D *diffFlowPtGenFunRe4 = dynamic_cast(outputListHistos->FindObject("fdiffFlowPtGenFunRe4")); - TProfile2D *diffFlowPtGenFunRe5 = dynamic_cast(outputListHistos->FindObject("fdiffFlowPtGenFunRe5")); - TProfile2D *diffFlowPtGenFunRe6 = dynamic_cast(outputListHistos->FindObject("fdiffFlowPtGenFunRe6")); - TProfile2D *diffFlowPtGenFunRe7 = dynamic_cast(outputListHistos->FindObject("fdiffFlowPtGenFunRe7")); - TProfile2D *diffFlowPtGenFunIm0 = dynamic_cast(outputListHistos->FindObject("fdiffFlowPtGenFunIm0")); - TProfile2D *diffFlowPtGenFunIm1 = dynamic_cast(outputListHistos->FindObject("fdiffFlowPtGenFunIm1")); - TProfile2D *diffFlowPtGenFunIm2 = dynamic_cast(outputListHistos->FindObject("fdiffFlowPtGenFunIm2")); - TProfile2D *diffFlowPtGenFunIm3 = dynamic_cast(outputListHistos->FindObject("fdiffFlowPtGenFunIm3")); - TProfile2D *diffFlowPtGenFunIm4 = dynamic_cast(outputListHistos->FindObject("fdiffFlowPtGenFunIm4")); - TProfile2D *diffFlowPtGenFunIm5 = dynamic_cast(outputListHistos->FindObject("fdiffFlowPtGenFunIm5")); - TProfile2D *diffFlowPtGenFunIm6 = dynamic_cast(outputListHistos->FindObject("fdiffFlowPtGenFunIm6")); - TProfile2D *diffFlowPtGenFunIm7 = dynamic_cast(outputListHistos->FindObject("fdiffFlowPtGenFunIm7")); - */ - - //profile with avarage selected multiplicity for int. flow - //TProfile *avMult = dynamic_cast(outputListHistos->FindObject("fAvMultIntFlow")); - - //profile with avarage values of Q-vector components (1st bin: , 2nd bin: , 3rd bin: <(Q_x)^2>, 4th bin: <(Q_y)^2>) - //TProfile *QVectorComponents = dynamic_cast(outputListHistos->FindObject("fQVectorComponents")); - - //q-distribution - //TH1D *qDist = dynamic_cast(outputListHistos->FindObject("fQDist")); - - //AliCumulantsFunctions finalResults(intFlowGenFun,NULL,NULL, intFlowResults,diffFlowResults2,diffFlowResults4,diffFlowResults6,diffFlowResults8,avMult,QVectorComponents,qDist,diffFlowPtGenFunRe0,diffFlowPtGenFunRe1,diffFlowPtGenFunRe2, diffFlowPtGenFunRe3,diffFlowPtGenFunRe4,diffFlowPtGenFunRe5,diffFlowPtGenFunRe6,diffFlowPtGenFunRe7,diffFlowPtGenFunIm0,diffFlowPtGenFunIm1, diffFlowPtGenFunIm2,diffFlowPtGenFunIm3,diffFlowPtGenFunIm4,diffFlowPtGenFunIm5,diffFlowPtGenFunIm6,diffFlowPtGenFunIm7); - - //AliCumulantsFunctions finalResults(intFlowGenFun,diffFlowPtGenFunRe,diffFlowPtGenFunIm, intFlowResults,diffFlowResults2,diffFlowResults4,diffFlowResults6,diffFlowResults8,avMult,QVectorComponents,qDist); - - //finalResults.Calculate(); - - //--------------------------------------------------- - - this->SetIntFlowResults(intFlowResults); - this->SetDiffFlowResults2nd(diffFlowResults2); - this->SetDiffFlowResults4th(diffFlowResults4); - this->SetDiffFlowResults6th(diffFlowResults6); - this->SetDiffFlowResults8th(diffFlowResults8); - - this->SetCommonHistsResults2nd(commonHistRes2nd); - this->SetCommonHistsResults4th(commonHistRes4th); - this->SetCommonHistsResults6th(commonHistRes6th); - this->SetCommonHistsResults8th(commonHistRes8th); - - this->SetCommonHists(commonHists); - - this->SetIntFlowGenFun(intFlowGenFun); - - this->SetIntFlowGenFun4(intFlowGenFun4); //only for other system of Eq. - this->SetIntFlowGenFun6(intFlowGenFun6); //only for other system of Eq. - this->SetIntFlowGenFun8(intFlowGenFun8); //only for other system of Eq. - this->SetIntFlowGenFun16(intFlowGenFun16); //only for other system of Eq. - - this->SetDiffFlowPtRPGenFunRe(diffFlowPtRPGenFunRe); - this->SetDiffFlowPtRPGenFunIm(diffFlowPtRPGenFunIm); - this->SetNumberOfParticlesPerPtBinRP(ptBinRPNoOfParticles); - - this->SetDiffFlowEtaRPGenFunRe(diffFlowEtaRPGenFunRe); - this->SetDiffFlowEtaRPGenFunIm(diffFlowEtaRPGenFunIm); - this->SetNumberOfParticlesPerEtaBinRP(etaBinRPNoOfParticles); - - this->SetDiffFlowPtPOIGenFunRe(diffFlowPtPOIGenFunRe); - this->SetDiffFlowPtPOIGenFunIm(diffFlowPtPOIGenFunIm); - this->SetNumberOfParticlesPerPtBinPOI(ptBinPOINoOfParticles); - - this->SetDiffFlowEtaPOIGenFunRe(diffFlowEtaPOIGenFunRe); - this->SetDiffFlowEtaPOIGenFunIm(diffFlowEtaPOIGenFunIm); - this->SetNumberOfParticlesPerEtaBinPOI(etaBinPOINoOfParticles); - - this->SetAverageMultiplicity(avMult); - - this->SetAverageMultiplicity4(avMult4); //only for other system of Eq. - this->SetAverageMultiplicity6(avMult6); //only for other system of Eq. - this->SetAverageMultiplicity8(avMult8); //only for other system of Eq. - this->SetAverageMultiplicity16(avMult16); //only for other system of Eq. + // Book all objects relevant for calculation of differential flow. + + // a) Define static constants for array's boundaries; + // b) Define local variables and local flags for booking; + // c) Book profile to hold all flags for differential flow; + // d) Book all event-by-event quantities; + // e) Book all profiles; + // f) Book all histograms. + + // a) Define static constants for array's boundaries: + static const Int_t pMax = 5; + static const Int_t qMax = 11; + + // b) Define local variables and local flags for booking: + Int_t nBinsPtEta[2] = {fnBinsPt,fnBinsEta}; + Double_t minPtEta[2] = {fPtMin,fEtaMin}; + Double_t maxPtEta[2] = {fPtMax,fEtaMax}; + TString reIm[2] = {"Re","Im"}; + TString rpPoi[2] = {"RP","POI"}; + TString ptEta[2] = {"p_{t}","#eta"}; + TString order[4] = {"2nd order","4th order","6th order","8th order"}; + + // c) Book profile to hold all flags for differential flow: + TString diffFlowFlagsName = "fDiffFlowFlags"; + fDiffFlowFlags = new TProfile(diffFlowFlagsName.Data(),"Flags for Differential Flow",1,0,1); + fDiffFlowFlags->SetTickLength(-0.01,"Y"); + fDiffFlowFlags->SetMarkerStyle(25); + fDiffFlowFlags->SetLabelSize(0.05); + fDiffFlowFlags->SetLabelOffset(0.02,"Y"); + fDiffFlowFlags->GetXaxis()->SetBinLabel(1,"..."); + fDiffFlowList->Add(fDiffFlowFlags); + + // d) Book all event-by-event quantities: + // ... (to be improved - perhaps not needed) + + // e) Book all profiles: + // Generating functions for differential flow: + for(Int_t ri=0;ri<2;ri++) + { + for(Int_t rp=0;rp<2;rp++) + { + for(Int_t pe=0;pe<2;pe++) + { + fDiffFlowGenFun[ri][rp][pe] = new TProfile3D(Form("fDiffFlowGenFun (%s, %s, %s)",reIm[ri].Data(),rpPoi[rp].Data(),ptEta[pe].Data()), + Form("#LT%s[D[%s-bin][p][q]]#GT for %ss",reIm[ri].Data(),ptEta[pe].Data(),rpPoi[rp].Data()), + nBinsPtEta[pe],minPtEta[pe],maxPtEta[pe],pMax,0.,(Double_t)pMax,qMax,0.,(Double_t)qMax); + fDiffFlowGenFun[ri][rp][pe]->SetXTitle(ptEta[pe].Data()); + fDiffFlowGenFun[ri][rp][pe]->SetYTitle("p"); + fDiffFlowGenFun[ri][rp][pe]->SetZTitle("q"); + fDiffFlowGenFun[ri][rp][pe]->SetTitleOffset(1.44,"X"); + fDiffFlowGenFun[ri][rp][pe]->SetTitleOffset(1.44,"Y"); + fDiffFlowProfiles->Add(fDiffFlowGenFun[ri][rp][pe]); + // to be improved - alternative // nBinsPtEta[pe],(Double_t)(fPtMin/fPtBinWidth),(Double_t)(fPtMax/fPtBinWidth),pMax,0.,(Double_t)pMax,qMax,0.,(Double_t)qMax); + } + } + } + // Number of particles in pt/eta bin for RPs/POIs: + for(Int_t rp=0;rp<2;rp++) + { + for(Int_t pe=0;pe<2;pe++) + { + fNoOfParticlesInBin[rp][pe] = new TProfile(Form("fNoOfParticlesInBin (%s, %s)",rpPoi[rp].Data(),ptEta[pe].Data()), + Form("Number of %ss per %s bin",rpPoi[rp].Data(),ptEta[pe].Data()), + nBinsPtEta[pe],minPtEta[pe],maxPtEta[pe]); + fNoOfParticlesInBin[rp][pe]->SetXTitle(ptEta[pe].Data()); + fDiffFlowProfiles->Add(fNoOfParticlesInBin[rp][pe]); + } + } + // Differential cumulants per pt/eta bin for RPs/POIs: + for(Int_t rp=0;rp<2;rp++) + { + for(Int_t pe=0;pe<2;pe++) + { + for(Int_t co=0;co<4;co++) + { + fDiffFlowCumulants[rp][pe][co] = new TH1D(Form("fDiffFlowCumulants (%s, %s, %s)",rpPoi[rp].Data(),ptEta[pe].Data(),order[co].Data()), + Form("Differential %s cumulant for %ss vs %s",order[co].Data(),rpPoi[rp].Data(),ptEta[pe].Data()), + nBinsPtEta[pe],minPtEta[pe],maxPtEta[pe]); + fDiffFlowCumulants[rp][pe][co]->SetXTitle(ptEta[pe].Data()); + fDiffFlowResults->Add(fDiffFlowCumulants[rp][pe][co]); + } + } + } + // Differential flow per pt/eta bin for RPs/POIs: + for(Int_t rp=0;rp<2;rp++) + { + for(Int_t pe=0;pe<2;pe++) + { + for(Int_t co=0;co<4;co++) + { + fDiffFlow[rp][pe][co] = new TH1D(Form("fDiffFlow (%s, %s, %s)",rpPoi[rp].Data(),ptEta[pe].Data(),order[co].Data()), + Form("Differential flow from %s cumulant for %ss vs %s",order[co].Data(),rpPoi[rp].Data(),ptEta[pe].Data()), + nBinsPtEta[pe],minPtEta[pe],maxPtEta[pe]); + fDiffFlow[rp][pe][co]->SetXTitle(ptEta[pe].Data()); + fDiffFlowResults->Add(fDiffFlow[rp][pe][co]); + } + } + } + +}// end of void AliFlowAnalysisWithCumulants::BookEverythingForDiffFlow() + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::StoreReferenceFlowFlags() +{ + // Store all flags for reference flow in profile fReferenceFlowFlags. + + if(!fReferenceFlowFlags) + { + cout<Fill(0.5,(Double_t)fUsePhiWeights||fUsePtWeights||fUseEtaWeights); + // Which event weight was used to weight generating function event-by-event: + if(strcmp(fMultiplicityWeight->Data(),"unit")) + { + fReferenceFlowFlags->Fill(1.5,0.); // 0 = "unit" (default) + } else if(strcmp(fMultiplicityWeight->Data(),"multiplicity")) + { + fReferenceFlowFlags->Fill(1.5,1.); // 1 = "multiplicity" + } + +} // end of void AliFlowAnalysisWithCumulants::StoreReferenceFlowFlags() + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::StoreDiffFlowFlags() +{ + // Store all flags for differential flow in profile fDiffFlowFlags. + + if(!fDiffFlowFlags) + { + cout<Fill(0.5,(Double_t) ... ); + +} // end of void AliFlowAnalysisWithCumulants::StoreDiffFlowFlags() + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::BookAndNestAllLists() +{ + // Book and nest all list in base list fHistList. + + // a) Book and nest lists for reference flow; + // b) Book and nest lists for differential flow; + // c) Book and nest lists for tuning. - this->SetQVectorComponents(qVectorComponents); + // a) Book and nest all lists for reference flow: + fReferenceFlowList = new TList(); + fReferenceFlowList->SetName("Reference Flow"); + fReferenceFlowList->SetOwner(kTRUE); + fHistList->Add(fReferenceFlowList); + fReferenceFlowProfiles = new TList(); + fReferenceFlowProfiles->SetName("Profiles"); + fReferenceFlowProfiles->SetOwner(kTRUE); + fReferenceFlowList->Add(fReferenceFlowProfiles); + fReferenceFlowResults = new TList(); + fReferenceFlowResults->SetName("Results"); + fReferenceFlowResults->SetOwner(kTRUE); + fReferenceFlowList->Add(fReferenceFlowResults); + // b) Book and nest lists for differential flow: + fDiffFlowList = new TList(); + fDiffFlowList->SetName("Differential Flow"); + fDiffFlowList->SetOwner(kTRUE); + fHistList->Add(fDiffFlowList); + fDiffFlowProfiles = new TList(); + fDiffFlowProfiles->SetName("Profiles"); + fDiffFlowProfiles->SetOwner(kTRUE); + fDiffFlowList->Add(fDiffFlowProfiles); + fDiffFlowResults = new TList(); + fDiffFlowResults->SetName("Results"); + fDiffFlowResults->SetOwner(kTRUE); + fDiffFlowList->Add(fDiffFlowResults); + // c) Book and nest lists for tuning: + if(fTuneParameters) + { + fTuningList = new TList(); + fTuningList->SetName("Tuning"); + fTuningList->SetOwner(kTRUE); + fHistList->Add(fTuningList); + fTuningProfiles = new TList(); + fTuningProfiles->SetName("Profiles"); + fTuningProfiles->SetOwner(kTRUE); + fTuningList->Add(fTuningProfiles); + fTuningResults = new TList(); + fTuningResults->SetName("Results"); + fTuningResults->SetOwner(kTRUE); + fTuningList->Add(fTuningResults); + } - this->SetAverageOfSquaredWeight(averageOfSquaredWeight); +} // end of void AliFlowAnalysisWithCumulants::BookAndNestAllLists() + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::BookProfileHoldingSettings() +{ + // Book profile to hold all analysis settings. + + TString analysisSettingsName = "fAnalysisSettings"; + fAnalysisSettings = new TProfile(analysisSettingsName.Data(),"Settings for analysis with Generating Function Cumulants",6,0.,6.); + fAnalysisSettings->GetXaxis()->SetLabelSize(0.035); + fAnalysisSettings->GetXaxis()->SetBinLabel(1,"Harmonic"); + fAnalysisSettings->Fill(0.5,fHarmonic); + fAnalysisSettings->GetXaxis()->SetBinLabel(2,"Multiple"); + fAnalysisSettings->Fill(1.5,fMultiple); + fAnalysisSettings->GetXaxis()->SetBinLabel(3,"r_{0}"); + fAnalysisSettings->Fill(2.5,fR0); + fAnalysisSettings->GetXaxis()->SetBinLabel(4,"Print RF results"); + fAnalysisSettings->Fill(3.5,fPrintFinalResults[0]); + fAnalysisSettings->GetXaxis()->SetBinLabel(5,"Print RP results"); + fAnalysisSettings->Fill(4.5,fPrintFinalResults[1]); + fAnalysisSettings->GetXaxis()->SetBinLabel(6,"Print POI results"); + fAnalysisSettings->Fill(5.5,fPrintFinalResults[2]); + fHistList->Add(fAnalysisSettings); + +} // end of void AliFlowAnalysisWithCumulants::BookProfileHoldingSettings() + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::BookCommonHistograms() +{ + // Book common control histograms and common histograms for final results. + + // Common control histogram: + TString commonHistsName = "AliFlowCommonHistGFC"; + fCommonHists = new AliFlowCommonHist(commonHistsName.Data()); + fHistList->Add(fCommonHists); + // Common histograms for final results from 2nd order GFC: + TString commonHistResults2ndOrderName = "AliFlowCommonHistResults2ndOrderGFC"; + fCommonHistsResults2nd = new AliFlowCommonHistResults(commonHistResults2ndOrderName.Data()); + fHistList->Add(fCommonHistsResults2nd); + // Common histograms for final results from 4th order GFC: + TString commonHistResults4thOrderName = "AliFlowCommonHistResults4thOrderGFC"; + fCommonHistsResults4th = new AliFlowCommonHistResults(commonHistResults4thOrderName.Data()); + fHistList->Add(fCommonHistsResults4th); + // Common histograms for final results from 6th order GFC: + TString commonHistResults6thOrderName = "AliFlowCommonHistResults6thOrderGFC"; + fCommonHistsResults6th = new AliFlowCommonHistResults(commonHistResults6thOrderName.Data()); + fHistList->Add(fCommonHistsResults6th); + // Common histograms for final results from 8th order GFC: + TString commonHistResults8thOrderName = "AliFlowCommonHistResults8thOrderGFC"; + fCommonHistsResults8th = new AliFlowCommonHistResults(commonHistResults8thOrderName.Data()); + fHistList->Add(fCommonHistsResults8th); + +} // end of void AliFlowAnalysisWithCumulants::BookCommonHistograms() + +//================================================================================================================ + +void AliFlowAnalysisWithCumulants::CheckPointersUsedInMake() +{ + // Check pointers used in method Make(). + + if(!fCommonHists) + { + cout<GetHistMultRP())) + { + cout<GetHistMultRP) is NULL in CPUIF() !!!!"<GetBinContent(1); + fMultiple = (Int_t)fAnalysisSettings->GetBinContent(2); + fR0 = (Double_t)fAnalysisSettings->GetBinContent(3); + fPrintFinalResults[0] = (Bool_t)fAnalysisSettings->GetBinContent(4); + fPrintFinalResults[1] = (Bool_t)fAnalysisSettings->GetBinContent(5); + fPrintFinalResults[2] = (Bool_t)fAnalysisSettings->GetBinContent(6); - AliCumulantsFunctions finalResults(fIntFlowGenFun,fIntFlowGenFun4,fIntFlowGenFun6,fIntFlowGenFun8,fIntFlowGenFun16,fAvMultIntFlow4GFC, fAvMultIntFlow6GFC,fAvMultIntFlow8GFC,fAvMultIntFlow16GFC,fDiffFlowPtRPGenFunRe,fDiffFlowPtRPGenFunIm,fPtBinRPNoOfParticles, fDiffFlowEtaRPGenFunRe,fDiffFlowEtaRPGenFunIm,fEtaBinRPNoOfParticles,fDiffFlowPtPOIGenFunRe,fDiffFlowPtPOIGenFunIm,fPtBinPOINoOfParticles, fDiffFlowEtaPOIGenFunRe,fDiffFlowEtaPOIGenFunIm,fEtaBinPOINoOfParticles,fIntFlowResultsGFC,fDiffFlowResults2ndOrderGFC,fDiffFlowResults4thOrderGFC,fDiffFlowResults6thOrderGFC,fDiffFlowResults8thOrderGFC, fAvMultIntFlowGFC,fQVectorComponentsGFC,fAverageOfSquaredWeight,fCommonHistsResults2nd, fCommonHistsResults4th,fCommonHistsResults6th,fCommonHistsResults8th,fCommonHists); - - finalResults.Calculate(); -} +} // end of AliFlowAnalysisWithCumulants::AccessSettings() //================================================================================================================ void AliFlowAnalysisWithCumulants::WriteHistograms(TString* outputFileName) { - //store the final results in output .root file + // Store the final results in output .root file. + TFile *output = new TFile(outputFileName->Data(),"RECREATE"); //output->WriteObject(fHistList, "cobjGFC","SingleKey"); fHistList->SetName("cobjGFC"); fHistList->SetOwner(kTRUE); fHistList->Write(fHistList->GetName(), TObject::kSingleKey); delete output; -} -//================================================================================================================ +} // end of void AliFlowAnalysisWithCumulants::WriteHistograms(TString* outputFileName) //================================================================================================================ void AliFlowAnalysisWithCumulants::WriteHistograms(TString outputFileName) { - //store the final results in output .root file + // Store the final results in output .root file. + TFile *output = new TFile(outputFileName.Data(),"RECREATE"); //output->WriteObject(fHistList, "cobjGFC","SingleKey"); fHistList->SetName("cobjGFC"); fHistList->SetOwner(kTRUE); fHistList->Write(fHistList->GetName(), TObject::kSingleKey); delete output; -} -//================================================================================================================ +} // end of void AliFlowAnalysisWithCumulants::WriteHistograms(TString outputFileName) //================================================================================================================ void AliFlowAnalysisWithCumulants::WriteHistograms(TDirectoryFile *outputFileName) { - //store the final results in output .root file + // Store the final results in output .root file. + fHistList->SetName("cobjGFC"); fHistList->SetOwner(kTRUE); outputFileName->Add(fHistList); outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey); -} + +} // end of void AliFlowAnalysisWithCumulants::WriteHistograms(TDirectoryFile *outputFileName) //================================================================================================================ + + diff --git a/PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithCumulants.h b/PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithCumulants.h index 4e86c7af78d..d46c42fed84 100644 --- a/PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithCumulants.h +++ b/PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithCumulants.h @@ -4,24 +4,27 @@ * $Id$ */ -/******************************** - * flow analysis with cumulants * - * * - * author: Ante Bilandzic * - * (anteb@nikhef.nl) * - *******************************/ +/************************************************* + * Flow analysis with cumulants. In this class * + * cumulants are calculated by making use of the * + * formalism of generating functions proposed by * + * Ollitrault et al. * + * * + * Author: Ante Bilandzic * + * (abilandzic@gmail.com) * + *************************************************/ #ifndef ALIFLOWANALYSISWITHCUMULANTS_H #define ALIFLOWANALYSISWITHCUMULANTS_H -#include "AliFlowCommonConstants.h" -#include "AliFlowCumuConstants.h" +#include "AliFlowCommonConstants.h" +#include "TMatrixD.h" -class TObjArray; class TList; class TFile; class TH1D; +class TH1F; class TProfile; class TProfile2D; class TProfile3D; @@ -39,238 +42,247 @@ class AliFlowAnalysisWithCumulants{ public: AliFlowAnalysisWithCumulants(); virtual ~AliFlowAnalysisWithCumulants(); - + // 0.) Methods called in the constructor: + virtual void InitializeArrays(); + // 1.) Method Init() and methods called within Init(): virtual void Init(); + virtual void CrossCheckSettings(); + virtual void AccessConstants(); + virtual void BookAndNestAllLists(); + virtual void BookProfileHoldingSettings(); + virtual void BookCommonHistograms(); + virtual void BookAndFillWeightsHistograms(); + virtual void BookEverythingForReferenceFlow(); + virtual void BookEverythingForDiffFlow(); + virtual void StoreReferenceFlowFlags(); + virtual void StoreDiffFlowFlags(); + virtual void BookEverythingForTuning(); + // 2.) Method Make() and methods called within Make(): virtual void Make(AliFlowEventSimple* anEvent); - virtual void GetOutputHistograms(TList *outputListHistos); //get pointers to all output histograms (called before Finish()) - virtual void Finish(); - virtual void WriteHistograms(TString* outputFileName); + virtual void CheckPointersUsedInMake(); + virtual void FillGeneratingFunctionForReferenceFlow(AliFlowEventSimple *anEvent); + virtual void FillQvectorComponents(AliFlowEventSimple *anEvent); + virtual void FillGeneratingFunctionForDiffFlow(AliFlowEventSimple *anEvent); + virtual void FillGeneratingFunctionsForDifferentTuningParameters(AliFlowEventSimple *anEvent); + // 3.) Method Finish() and methods called within Finish(): + virtual void Finish(); + virtual void CheckPointersUsedInFinish(); + virtual void AccessSettings(); + virtual void GetAvMultAndNoOfEvts(); + virtual void CalculateCumulantsForReferenceFlow(); + virtual void CalculateReferenceFlow(); + virtual void CalculateReferenceFlowError(); + virtual void FillCommonHistResultsForReferenceFlow(); + virtual void CalculateCumulantsForDiffFlow(TString rpPoi,TString ptEta); + virtual void CalculateDifferentialFlow(TString rpPoi,TString ptEta); + virtual void CalculateDifferentialFlowErrors(TString rpPoi,TString ptEta); + virtual void FillCommonHistResultsForDifferentialFlow(TString rpPoi); + virtual void CalculateIntegratedFlow(TString rpPoi); // to be improved (add also possibility to integrate over eta yield) + virtual void PrintFinalResults(TString rpPoi); + virtual void FinalizeTuning(); + // 4.) Method GetOutputHistograms() and method called within it: + virtual void GetOutputHistograms(TList *outputListHistos); + virtual void GetPointersForBaseHistograms(); + virtual void GetPointersForCommonControlHistograms(); + virtual void GetPointersForCommonResultsHistograms(); + virtual void GetPointersForParticleWeightsHistograms(); + virtual void GetPointersForReferenceFlowObjects(); + virtual void GetPointersForDiffFlowObjects(); + virtual void WriteHistograms(TString *outputFileName); virtual void WriteHistograms(TString outputFileName); virtual void WriteHistograms(TDirectoryFile *outputFileName); - -//---------------------------------------------------------------------------------------------------------------- -// setters and getters -//---------------------------------------------------------------------------------------------------------------- - TList* GetHistList() const {return this->fHistList;} - - void SetWeightsList(TList* wlist) {this->fWeightsList = wlist;} - TList* GetWeightsList() const {return this->fWeightsList;} - - void SetIntFlowResults(TH1D* const ifr) {this->fIntFlowResultsGFC = ifr;}; - TH1D* GetIntFlowResults() const {return this->fIntFlowResultsGFC;}; - - void SetDiffFlowResults2nd(TH1D* const diff2nd) {this->fDiffFlowResults2ndOrderGFC = diff2nd;}; - TH1D* GetDiffFlowResults2nd() const {return this->fDiffFlowResults2ndOrderGFC;}; - - void SetDiffFlowResults4th(TH1D* const diff4th) {this->fDiffFlowResults4thOrderGFC = diff4th;}; - TH1D* GetDiffFlowResults4th() const {return this->fDiffFlowResults4thOrderGFC;}; - - void SetDiffFlowResults6th(TH1D* const diff6th) {this->fDiffFlowResults6thOrderGFC = diff6th;}; - TH1D* GetDiffFlowResults6th() const {return this->fDiffFlowResults6thOrderGFC;}; - - void SetDiffFlowResults8th(TH1D* const diff8th) {this->fDiffFlowResults8thOrderGFC = diff8th;}; - TH1D* GetDiffFlowResults8th() const {return this->fDiffFlowResults8thOrderGFC;}; - - void SetCommonHistsResults2nd(AliFlowCommonHistResults* const chr2nd) {this->fCommonHistsResults2nd = chr2nd;}; - AliFlowCommonHistResults* GetCommonHistsResults2nd() const {return this->fCommonHistsResults2nd;}; - - void SetCommonHistsResults4th(AliFlowCommonHistResults* const chr4th) {this->fCommonHistsResults4th = chr4th;}; - AliFlowCommonHistResults* GetCommonHistsResults4th() const {return this->fCommonHistsResults4th;}; - - void SetCommonHistsResults6th(AliFlowCommonHistResults* const chr6th) {this->fCommonHistsResults6th = chr6th;}; - AliFlowCommonHistResults* GetCommonHistsResults6th() const {return this->fCommonHistsResults6th;}; - - void SetCommonHistsResults8th(AliFlowCommonHistResults* const chr8th) {this->fCommonHistsResults8th = chr8th;}; - AliFlowCommonHistResults* GetCommonHistsResults8th() const {return this->fCommonHistsResults8th;}; - - void SetIntFlowGenFun(TProfile2D* const ifgf) {this->fIntFlowGenFun = ifgf;}; - TProfile2D* GetIntFlowGenFun() const {return this->fIntFlowGenFun;}; - - void SetIntFlowGenFun4(TProfile2D* const ifgf4) {this->fIntFlowGenFun4 = ifgf4;}; //(only for other system of Eq.) - TProfile2D* GetIntFlowGenFun4() const {return this->fIntFlowGenFun4;}; //(only for other system of Eq.) - - void SetIntFlowGenFun6(TProfile2D* const ifgf6) {this->fIntFlowGenFun6 = ifgf6;}; //(only for other system of Eq.) - TProfile2D* GetIntFlowGenFun6() const {return this->fIntFlowGenFun6;}; //(only for other system of Eq.) - - void SetIntFlowGenFun8(TProfile2D* const ifgf8) {this->fIntFlowGenFun8 = ifgf8;}; //(only for other system of Eq.) - TProfile2D* GetIntFlowGenFun8() const {return this->fIntFlowGenFun8;}; //(only for other system of Eq.) - - void SetIntFlowGenFun16(TProfile2D* const ifgf16) {this->fIntFlowGenFun16 = ifgf16;}; //(only for other system of Eq.) - TProfile2D* GetIntFlowGenFun16() const {return this->fIntFlowGenFun16;}; //(only for other system of Eq.) - - void SetAverageMultiplicity4(TProfile* const am4) {this->fAvMultIntFlow4GFC = am4;}; //(only for other system of Eq.) - TProfile* GetAverageMultiplicity4() const {return this->fAvMultIntFlow4GFC;}; //(only for other system of Eq.) - - void SetAverageMultiplicity6(TProfile* const am6) {this->fAvMultIntFlow6GFC = am6;}; //(only for other system of Eq.) - TProfile* GetAverageMultiplicity6() const {return this->fAvMultIntFlow6GFC;}; //(only for other system of Eq.) - - void SetAverageMultiplicity8(TProfile* const am8) {this->fAvMultIntFlow8GFC = am8;}; //(only for other system of Eq.) - TProfile* GetAverageMultiplicity8() const {return this->fAvMultIntFlow8GFC;}; //(only for other system of Eq.) - - void SetAverageMultiplicity16(TProfile* const am16) {this->fAvMultIntFlow16GFC = am16;}; //(only for other system of Eq.) - TProfile* GetAverageMultiplicity16() const {return this->fAvMultIntFlow16GFC;}; //(only for other system of Eq.) - - //---------------------------------------------------------------------------------------------------------- - //RP = Reaction Plane (RP denotes particles used to determine the reaction plane) - void SetDiffFlowPtRPGenFunRe(TProfile3D* const dfRPPtgfRe) {this->fDiffFlowPtRPGenFunRe = dfRPPtgfRe;}; - TProfile3D* GetDiffFlowPtRPGenFunRe() const {return this->fDiffFlowPtRPGenFunRe;}; - - void SetDiffFlowPtRPGenFunIm(TProfile3D* const dfRPPtgfIm) {this->fDiffFlowPtRPGenFunIm = dfRPPtgfIm;}; - TProfile3D* GetDiffPtRPFlowGenFunIm() const {return this->fDiffFlowPtRPGenFunIm;}; - - void SetNumberOfParticlesPerPtBinRP(TProfile* const nopppbRP) {this->fPtBinRPNoOfParticles = nopppbRP;}; - TProfile* GetNumberOfParticlesPerPtBinRP() const {return this->fPtBinRPNoOfParticles;}; - - void SetDiffFlowEtaRPGenFunRe(TProfile3D* const dfEtaRPgfRe) {this->fDiffFlowEtaRPGenFunRe = dfEtaRPgfRe;}; - TProfile3D* GetDiffFlowEtaRPGenFunRe() const {return this->fDiffFlowEtaRPGenFunRe;}; - - void SetDiffFlowEtaRPGenFunIm(TProfile3D* const dfEtaRPgfIm) {this->fDiffFlowEtaRPGenFunIm = dfEtaRPgfIm;}; - TProfile3D* GetDiffFlowEtaRPGenFunIm() const {return this->fDiffFlowEtaRPGenFunIm;}; - - void SetNumberOfParticlesPerEtaBinRP(TProfile* const noppebRP) {this->fEtaBinRPNoOfParticles = noppebRP;}; - TProfile* GetNumberOfParticlesPerEtaBinRP() const {return this->fEtaBinRPNoOfParticles;}; - //---------------------------------------------------------------------------------------------------------- - - //---------------------------------------------------------------------------------------------------------- - //POI = Particles Of Interest (POI denotes particles of interest for the final physical results for int. and diff. flow) - void SetDiffFlowPtPOIGenFunRe(TProfile3D* const dfPOIPtgfRe) {this->fDiffFlowPtPOIGenFunRe = dfPOIPtgfRe;}; - TProfile3D* GetDiffFlowPtPOIGenFunRe() const {return this->fDiffFlowPtPOIGenFunRe;}; - - void SetDiffFlowPtPOIGenFunIm(TProfile3D* const dfPOIPtgfIm) {this->fDiffFlowPtPOIGenFunIm = dfPOIPtgfIm;}; - TProfile3D* GetDiffPtPOIFlowGenFunIm() const {return this->fDiffFlowPtPOIGenFunIm;}; - - void SetNumberOfParticlesPerPtBinPOI(TProfile* const nopppbPOI) {this->fPtBinPOINoOfParticles = nopppbPOI;}; - TProfile* GetNumberOfParticlesPerPtBinPOI() const {return this->fPtBinPOINoOfParticles;}; - - void SetDiffFlowEtaPOIGenFunRe(TProfile3D* const dfEtaPOIgfRe) {this->fDiffFlowEtaPOIGenFunRe = dfEtaPOIgfRe;}; - TProfile3D* GetDiffFlowEtaPOIGenFunRe() const {return this->fDiffFlowEtaPOIGenFunRe;}; - - void SetDiffFlowEtaPOIGenFunIm(TProfile3D* const dfEtaPOIgfIm) {this->fDiffFlowEtaPOIGenFunIm = dfEtaPOIgfIm;}; - TProfile3D* GetDiffFlowEtaPOIGenFunIm() const {return this->fDiffFlowEtaPOIGenFunIm;}; - - void SetNumberOfParticlesPerEtaBinPOI(TProfile* const noppebPOI) {this->fEtaBinPOINoOfParticles = noppebPOI;}; - TProfile* GetNumberOfParticlesPerEtaBinPOI() const {return this->fEtaBinPOINoOfParticles;}; - //---------------------------------------------------------------------------------------------------------- - - void SetAverageMultiplicity(TProfile* const am) {this->fAvMultIntFlowGFC = am;}; - TProfile* GetAverageMultiplicity() const {return this->fAvMultIntFlowGFC;}; - - void SetQVectorComponents(TProfile* const sqvc) {this->fQVectorComponentsGFC = sqvc;}; - TProfile* GetQVectorComponents() const {return this->fQVectorComponentsGFC;}; - - void SetCommonHists(AliFlowCommonHist* const ch) {this->fCommonHists = ch;}; - AliFlowCommonHist* GetCommonHists() const {return this->fCommonHists;}; - + // 5.) Other methods: + // ... + // 6.) Setters and getters: + // 6.0.) base: + void SetHistList(TList* const hl) {this->fHistList = hl;} + TList* GetHistList() const {return this->fHistList;} + void SetHistListName(const char *hln) {this->fHistListName->Append(*hln);}; + TString *GetHistListName() const {return this->fHistListName;}; + void SetAnalysisSettings(TProfile* const as) {this->fAnalysisSettings = as;}; + TProfile* GetAnalysisSettings() const {return this->fAnalysisSettings;}; + // 6.1.) common: + void SetCommonHists(AliFlowCommonHist* const ch) {this->fCommonHists = ch;}; + AliFlowCommonHist* GetCommonHists() const {return this->fCommonHists;}; + void SetCommonHistsResults2nd(AliFlowCommonHistResults* const chr2nd) {this->fCommonHistsResults2nd = chr2nd;}; + AliFlowCommonHistResults* GetCommonHistsResults2nd() const {return this->fCommonHistsResults2nd;}; + void SetCommonHistsResults4th(AliFlowCommonHistResults* const chr4th) {this->fCommonHistsResults4th = chr4th;}; + AliFlowCommonHistResults* GetCommonHistsResults4th() const {return this->fCommonHistsResults4th;}; + void SetCommonHistsResults6th(AliFlowCommonHistResults* const chr6th) {this->fCommonHistsResults6th = chr6th;}; + AliFlowCommonHistResults* GetCommonHistsResults6th() const {return this->fCommonHistsResults6th;}; + void SetCommonHistsResults8th(AliFlowCommonHistResults* const chr8th) {this->fCommonHistsResults8th = chr8th;}; + AliFlowCommonHistResults* GetCommonHistsResults8th() const {return this->fCommonHistsResults8th;}; + void SetHarmonic(Int_t const harmonic) {this->fHarmonic = harmonic;}; + Int_t GetHarmonic() const {return this->fHarmonic;}; + void SetMultiple(Int_t const multiple) {this->fMultiple = multiple;}; + Int_t GetMultiple() const {return this->fMultiple;}; + void SetR0(Double_t const r0) {this->fR0 = r0;}; + Double_t GetR0() const {return this->fR0;}; + void SetPrintFinalResults(Bool_t const printOrNot, Int_t const i) {this->fPrintFinalResults[i] = printOrNot;}; + Bool_t GetPrintFinalResults(Int_t i) const {return this->fPrintFinalResults[i];}; + // 6.2.0.) particle weights: + void SetWeightsList(TList* const wlist) {this->fWeightsList = (TList*)wlist->Clone();} + TList* GetWeightsList() const {return this->fWeightsList;} void SetUsePhiWeights(Bool_t const uPhiW) {this->fUsePhiWeights = uPhiW;}; Bool_t GetUsePhiWeights() const {return this->fUsePhiWeights;}; - void SetUsePtWeights(Bool_t const uPtW) {this->fUsePtWeights = uPtW;}; Bool_t GetUsePtWeights() const {return this->fUsePtWeights;}; - void SetUseEtaWeights(Bool_t const uEtaW) {this->fUseEtaWeights = uEtaW;}; Bool_t GetUseEtaWeights() const {return this->fUseEtaWeights;}; - + void SetUseParticleWeights(TProfile* const uPW) {this->fUseParticleWeights = uPW;}; + TProfile* GetUseParticleWeights() const {return this->fUseParticleWeights;}; + void SetPhiWeights(TH1F* const histPhiWeights) {this->fPhiWeights = histPhiWeights;}; + TH1F* GetPhiWeights() const {return this->fPhiWeights;}; + void SetPtWeights(TH1D* const histPtWeights) {this->fPtWeights = histPtWeights;}; + TH1D* GetPtWeights() const {return this->fPtWeights;}; + void SetEtaWeights(TH1D* const histEtaWeights) {this->fEtaWeights = histEtaWeights;}; + TH1D* GetEtaWeights() const {return this->fEtaWeights;}; + // 6.2.1.) event weights: + void SetMultiplicityWeight(const char *multiplicityWeight) {*this->fMultiplicityWeight = multiplicityWeight;}; + // 6.3.) reference flow: + void SetReferenceFlowFlags(TProfile* const rff) {this->fReferenceFlowFlags = rff;}; + TProfile* GetReferenceFlowFlags() const {return this->fReferenceFlowFlags;}; + void SetnBinsMult(Int_t const nbm) {this->fnBinsMult = nbm;}; + Int_t GetnBinsMult() const {return this->fnBinsMult;}; + void SetMinMult(Double_t const minm) {this->fMinMult = minm;}; + Double_t GetMinMult() const {return this->fMinMult;}; + void SetMaxMult(Double_t const maxm) {this->fMaxMult = maxm;}; + Double_t GetMaxMult() const {return this->fMaxMult;}; + // 6.3.0.) profiles: + void SetReferenceFlowGenFun(TProfile2D* const rfgf) {this->fReferenceFlowGenFun = rfgf;}; + TProfile2D* GetReferenceFlowGenFun() const {return this->fReferenceFlowGenFun;}; + void SetQvectorComponents(TProfile* const qvc) {this->fQvectorComponents = qvc;}; + TProfile* GetQvectorComponents() const {return this->fQvectorComponents;}; void SetAverageOfSquaredWeight(TProfile* const aosw) {this->fAverageOfSquaredWeight = aosw;}; TProfile* GetSumOfSquaredWeight() const {return this->fAverageOfSquaredWeight;}; -//---------------------------------------------------------------------------------------------------------------- - + // 6.3.1.) results: + void SetReferenceFlowCumulants(TH1D* const rfc) {this->fReferenceFlowCumulants = rfc;}; + TH1D* GetReferenceFlowCumulants() const {return this->fReferenceFlowCumulants;}; + void SetReferenceFlow(TH1D* const rf) {this->fReferenceFlow = rf;}; + TH1D* GetReferenceFlow() const {return this->fReferenceFlow;}; + void SetChi(TH1D* const c) {this->fChi = c;}; + TH1D* GetChi() const {return this->fChi;}; + // 6.4.) differential flow: + void SetDiffFlowFlags(TProfile* const dff) {this->fDiffFlowFlags = dff;}; + TProfile* GetDiffFlowFlags() const {return this->fDiffFlowFlags;}; + // 6.4.0.) profiles: + void SetDiffFlowGenFun(TProfile3D* const dfgf, Int_t const ri, Int_t const rp, Int_t const pe) {this->fDiffFlowGenFun[ri][rp][pe] = dfgf;}; + TProfile3D* GetDiffFlowGenFun(Int_t const ri, Int_t const rp, Int_t const pe) const {return this->fDiffFlowGenFun[ri][rp][pe];}; + void SetNoOfParticlesInBin(TProfile* const nopib, Int_t const rp, Int_t const pe) {this->fNoOfParticlesInBin[rp][pe] = nopib;}; + TProfile* GetNoOfParticlesInBin(Int_t const rp, Int_t const pe) const {return this->fNoOfParticlesInBin[rp][pe];}; + // 6.4.1.) results: + void SetDiffFlowCumulants(TH1D* const dfc, Int_t const rp, Int_t const pe, Int_t const co) {this->fDiffFlowCumulants[rp][pe][co] = dfc;}; + TH1D* GetDiffFlowCumulants(Int_t const rp, Int_t const pe, Int_t const co) const {return this->fDiffFlowCumulants[rp][pe][co];}; + void SetDiffFlow(TH1D* const df, Int_t const rp, Int_t const pe, Int_t const co) {this->fDiffFlow[rp][pe][co] = df;}; + TH1D* GetDiffFlow(Int_t const rp, Int_t const pe, Int_t const co) const {return this->fDiffFlow[rp][pe][co];}; + // 6.x.) Tuning the interpolating parameter r0 and using cutoff at different order in series: + void SetTuningFlags(TProfile* const tf) {this->fTuningFlags = tf;}; + TProfile* GetTuningFlags() const {return this->fTuningFlags;}; + void SetTuneParameters(Bool_t const tp) {this->fTuneParameters = tp;}; + Bool_t GetTuneParameters() const {return this->fTuneParameters;}; + void SetTuningR0(Double_t const tr0, Int_t const r) {this->fTuningR0[r] = tr0;}; + Double_t GetTuningR0(Int_t const r) const {return this->fTuningR0[r];}; + // 6.x.0.) profiles: + void SetTuningGenFun(TProfile2D* const tgf, Int_t const r, Int_t const pq) {this->fTuningGenFun[r][pq] = tgf;}; + TProfile2D* GetTuningGenFun(Int_t const r, Int_t const pq) const {return this->fTuningGenFun[r][pq];}; + // 6.x.1.) results: + void SetTuningCumulants(TH1D* const tc, Int_t const r, Int_t const pq) {this->fTuningCumulants[r][pq] = tc;}; + TH1D* GetTuningCumulants(Int_t const r, Int_t const pq) const {return this->fTuningCumulants[r][pq];}; + void SetTuningFlow(TH1D* const tf, Int_t const r, Int_t const pq) {this->fTuningFlow[r][pq] = tf;}; + TH1D* GetTuningFlow(Int_t const r, Int_t const pq) const {return this->fTuningFlow[r][pq];}; + private: AliFlowAnalysisWithCumulants(const AliFlowAnalysisWithCumulants& afawc); - AliFlowAnalysisWithCumulants& operator=(const AliFlowAnalysisWithCumulants& afawc); - AliFlowTrackSimple* fTrack; //track - static const Int_t fgkQmax = AliFlowCumuConstants::fgkQmax; //needed for numerics - static const Int_t fgkPmax = AliFlowCumuConstants::fgkPmax; //needed for numerics - static const Int_t fgkQmax4 = AliFlowCumuConstants::fgkQmax4; //needed for numerics (only for different system of Eq.) - static const Int_t fgkPmax4 = AliFlowCumuConstants::fgkPmax4; //needed for numerics (only for different system of Eq.) - static const Int_t fgkQmax6 = AliFlowCumuConstants::fgkQmax6; //needed for numerics (only for different system of Eq.) - static const Int_t fgkPmax6 = AliFlowCumuConstants::fgkPmax6; //needed for numerics (only for different system of Eq.) - static const Int_t fgkQmax8 = AliFlowCumuConstants::fgkQmax8; //needed for numerics (only for different system of Eq.) - static const Int_t fgkPmax8 = AliFlowCumuConstants::fgkPmax8; //needed for numerics (only for different system of Eq.) - static const Int_t fgkQmax16 = AliFlowCumuConstants::fgkQmax16; //needed for numerics (only for different system of Eq.) - static const Int_t fgkPmax16 = AliFlowCumuConstants::fgkPmax16; //needed for numerics (only for different system of Eq.) - static const Int_t fgkFlow = AliFlowCumuConstants::fgkFlow; //integrated flow coefficient to be calculated - static const Int_t fgkMltpl = AliFlowCumuConstants::fgkMltpl; //the multiple in p=m*n (diff. flow) - TList* fHistList; //list to hold all output histograms - TList* fWeightsList; //list to hold all histograms with weights - - Double_t fR0; //needed for numerics - Double_t fPtMax; //maximum pt - Double_t fPtMin; //minimum pt - Double_t fBinWidthPt; //width of pt bin (in GeV) - Int_t fgknBinsPt; //number of pt bins - - Double_t fEtaMax; //maximum pt - Double_t fEtaMin; //minimum pt - Double_t fBinWidthEta; //width of pt bin (in GeV) - Int_t fgknBinsEta; //number of pt bins - - Double_t fAvQx; // - Double_t fAvQy; // - Double_t fAvQ2x; //<(Q_x)^2> - Double_t fAvQ2y; //<(Q_y)^2> - - TProfile* fAvMultIntFlowGFC; //average selected multiplicity - - TProfile* fQVectorComponentsGFC; //averages of Q-vector components (1st bin: , 2nd bin: , 3rd bin: <(Q_x)^2>, 4th bin: <(Q_y)^2>) + AliFlowAnalysisWithCumulants& operator=(const AliFlowAnalysisWithCumulants& afawc); + // 0.) Base: + TList *fHistList; // base list to hold all output objects + TString *fHistListName; // name of base list + TProfile *fAnalysisSettings; // profile to hold analysis settings + // 1.) Common: + AliFlowCommonHist *fCommonHists; // common control histograms (filled only with events with 3 or more tracks for 3-p correlators) + AliFlowCommonHistResults *fCommonHistsResults2nd; // common result histograms for 2nd order cumulant + AliFlowCommonHistResults *fCommonHistsResults4th; // common result histograms for 4th order cumulant + AliFlowCommonHistResults *fCommonHistsResults6th; // common result histograms for 6th order cumulant + AliFlowCommonHistResults *fCommonHistsResults8th; // common result histograms for 8th order cumulant + Int_t fnBinsPhi; // number of phi bins + Double_t fPhiMin; // minimum phi + Double_t fPhiMax; // maximum phi + Double_t fPhiBinWidth; // bin width for phi histograms + Int_t fnBinsPt; // number of pt bins + Double_t fPtMin; // minimum pt + Double_t fPtMax; // maximum pt + Double_t fPtBinWidth; // bin width for pt histograms + Int_t fnBinsEta; // number of eta bins + Double_t fEtaMin; // minimum eta + Double_t fEtaMax; // maximum eta + Double_t fEtaBinWidth; // bin width for eta histograms + Int_t fHarmonic; // harmonic + Int_t fMultiple; // the multiple m in p=m*n, where n is harmonic (relevant for differential flow) + Double_t fR0; // r_{0} parameter + Bool_t fPrintFinalResults[3]; // print on the screen the final results [0=RF,1=RP,2=POI] + // 2a.) Particle weights: + TList *fWeightsList; // list to hold all histograms with particle weights: fUseParticleWeights, fPhiWeights, fPtWeights and fEtaWeights + Bool_t fUsePhiWeights; // use phi weights + Bool_t fUsePtWeights; // use pt weights + Bool_t fUseEtaWeights; // use eta weights + TProfile *fUseParticleWeights; // profile with three bins to hold values of fUsePhiWeights, fUsePtWeights and fUseEtaWeights + TH1F *fPhiWeights; // histogram holding phi weights + TH1D *fPtWeights; // histogram holding phi weights + TH1D *fEtaWeights; // histogram holding phi weights + // 2b.) Event weights: + TString *fMultiplicityWeight; // event-by-event weight for reference flow generating function (can be "unit" or "multiplicity") + // 3.) Reference flow: + // 3a.) lists: + TList *fReferenceFlowList; // list to hold all histograms and profiles relevant for reference flow + TList *fReferenceFlowProfiles; // list to hold all profiles relevant for reference flow + TList *fReferenceFlowResults; // list to hold all histograms with final results relevant for reference flow + // 3b.) flags: + TProfile *fReferenceFlowFlags; // profile to hold all flags for reference flow + Int_t fnBinsMult; // number of multiplicity bins for flow analysis versus multiplicity + Double_t fMinMult; // minimal multiplicity for flow analysis versus multiplicity + Double_t fMaxMult; // maximal multiplicity for flow analysis versus multiplicity + // 3c.) event-by-event quantities: + TMatrixD *fGEBE; // reference flow generating function only for current event + // 3d.) profiles: + TProfile2D *fReferenceFlowGenFun; // all-event average of the generating function used to calculate reference flow + TProfile *fQvectorComponents; // averages of Q-vector components (1st bin: , 2nd bin: , 3rd bin: <(Q_x)^2>, 4th bin: <(Q_y)^2>) + TProfile *fAverageOfSquaredWeight; // <>, where w = wPhi*wPt*wEta + // 3e.) results: + Double_t fAvM; // average multiplicity + Int_t fnEvts; // number of events + TH1D *fReferenceFlowCumulants; // final results for isotropic cumulants for reference flow + TH1D *fReferenceFlow; // final results for reference flow + TH1D *fChi; // final results for resolution + // 4.) Differential flow: + // 4a.) lists: + TList *fDiffFlowList; // list to hold all histograms and profiles relevant for differential flow + TList *fDiffFlowProfiles; // list to hold all profiles relevant for differential flow + TList *fDiffFlowResults; // list to hold all histograms with final results relevant for differential flow + // 4b.) flags: + TProfile *fDiffFlowFlags; // profile to hold all flags for reference flow + // 4c.) profiles: + TProfile3D *fDiffFlowGenFun[2][2][2]; // all-event avarage of generating function used for differential flow [0=Re,1=Im][0=RP,1=POI][0=pt,1=eta] + TProfile *fNoOfParticlesInBin[2][2]; // number of particles in pt/eta bin for RPs/POIs [0=RP,1=POI][0=pt,1=eta] + // 4d.) results: + TH1D *fDiffFlowCumulants[2][2][4]; // isotropic differential flow cumulants [0=RP,1=POI][0=pt,1=eta][cumulant order] + TH1D *fDiffFlow[2][2][4]; // differential flow [0=RP,1=POI][0=pt,1=eta][cumulant order] + // x.) Tuning the interpolating parameter r0 and using cutoff at different order in series: + // xa.) lists: + TList *fTuningList; // list to hold all histograms and profiles relevant for tuning + TList *fTuningProfiles; // list to hold all profiles relevant for tuning + TList *fTuningResults; // list to hold all histograms with final results relevant for tuning + // xb.) flags: + TProfile *fTuningFlags; // profile to hold all flags for tuning + Bool_t fTuneParameters; // tune r0 and cut series at different order + Double_t fTuningR0[10]; // different r0 values (at maximum 10 different values allowed) + // xc.) profiles: + TProfile2D *fTuningGenFun[10][5]; // generating function G evaluated for 10 different r0s and 5 different sets of (pmax,qmax) + // xd.) results: + TH1D *fTuningCumulants[10][5]; // isotropic cumulants for reference flow for 10 different r0s and 5 different sets of (pmax,qmax) + TH1D *fTuningFlow[10][5]; // reference flow for 10 different r0s and 5 different sets of (pmax,qmax) - TH1D* fIntFlowResultsGFC; //integrated flow final results - - TH1D* fDiffFlowResults2ndOrderGFC; //differential flow final results (2nd order estimate) - TH1D* fDiffFlowResults4thOrderGFC; //differential flow final results (4th order estimate) - TH1D* fDiffFlowResults6thOrderGFC; //differential flow final results (6th order estimate) - TH1D* fDiffFlowResults8thOrderGFC; //differential flow final results (8th order estimate) - - AliFlowCommonHistResults* fCommonHistsResults2nd; //final results for 2nd order int. and diff. flow stored in the common histograms - AliFlowCommonHistResults* fCommonHistsResults4th; //final results for 4th order int. and diff. flow stored in the common histograms - AliFlowCommonHistResults* fCommonHistsResults6th; //final results for 6th order int. and diff. flow stored in the common histograms - AliFlowCommonHistResults* fCommonHistsResults8th; //final results for 8th order int. and diff. flow stored in the common histograms - - TProfile2D* fIntFlowGenFun; //avarage of the generating function for integrated flow - - TProfile2D* fIntFlowGenFun4; //avarage of the generating function for integrated flow (only for different system of Eq.) - TProfile2D* fIntFlowGenFun6; //avarage of the generating function for integrated flow (only for different system of Eq.) - TProfile2D* fIntFlowGenFun8; //avarage of the generating function for integrated flow (only for different system of Eq.) - TProfile2D* fIntFlowGenFun16; //avarage of the generating function for integrated flow (only for different system of Eq.) - - TProfile* fAvMultIntFlow4GFC; //average selected multiplicity (only for different system of Eq.) - TProfile* fAvMultIntFlow6GFC; //average selected multiplicity (only for different system of Eq.) - TProfile* fAvMultIntFlow8GFC; //average selected multiplicity (only for different system of Eq.) - TProfile* fAvMultIntFlow16GFC; //average selected multiplicity (only for different system of Eq.) - - //RP = Reaction Plane (RP denotes particles used to determine the reaction plane) - TProfile3D* fDiffFlowPtRPGenFunRe; //avarage of the generating function for differential flow in Pt (real part) - TProfile3D* fDiffFlowPtRPGenFunIm; //avarage of the generating function for differential flow in Pt (imaginary part) - TProfile* fPtBinRPNoOfParticles; //number of particles per pt bin - TProfile3D* fDiffFlowEtaRPGenFunRe; //avarage of the generating function for differential flow in Eta (real part) - TProfile3D* fDiffFlowEtaRPGenFunIm; //avarage of the generating function for differential flow in Eta (imaginary part) - TProfile* fEtaBinRPNoOfParticles; //number of particles per eta bin - - //POI = Particles Of Interest (POI denotes particles of interest for the final physical results for int. and diff. flow) - TProfile3D* fDiffFlowPtPOIGenFunRe; //avarage of the generating function for differential flow in Pt (real part) - TProfile3D* fDiffFlowPtPOIGenFunIm; //avarage of the generating function for differential flow in Pt (imaginary part) - TProfile* fPtBinPOINoOfParticles; //number of particles per pt bin - TProfile3D* fDiffFlowEtaPOIGenFunRe; //avarage of the generating function for differential flow in Eta (real part) - TProfile3D* fDiffFlowEtaPOIGenFunIm; //avarage of the generating function for differential flow in Eta (imaginary part) - TProfile* fEtaBinPOINoOfParticles; //number of particles per eta bin - - /* - TProfile2D *fDiffFlowGenFunRe0,*fDiffFlowGenFunRe1,*fDiffFlowGenFunRe2,*fDiffFlowGenFunRe3;//differential flow - TProfile2D *fDiffFlowGenFunRe4,*fDiffFlowGenFunRe5,*fDiffFlowGenFunRe6,*fDiffFlowGenFunRe7;//differential flow - TProfile2D *fDiffFlowGenFunIm0,*fDiffFlowGenFunIm1,*fDiffFlowGenFunIm2,*fDiffFlowGenFunIm3;//differential flow - TProfile2D *fDiffFlowGenFunIm4,*fDiffFlowGenFunIm5,*fDiffFlowGenFunIm6,*fDiffFlowGenFunIm7;//differential flow - */ - - AliFlowCommonHist* fCommonHists; //common control histograms - - Bool_t fOtherEquations; //numerical equations for cumulants solved up to different highest order - - Bool_t fUsePhiWeights; //phi weights - Bool_t fUsePtWeights; //v_2(pt) weights - Bool_t fUseEtaWeights; //v_2(eta) weights - - TProfile* fAverageOfSquaredWeight; // - ClassDef(AliFlowAnalysisWithCumulants, 0); + }; //================================================================================================================ diff --git a/PWG2/FLOW/AliFlowCommon/AliFlowCumuConstants.cxx b/PWG2/FLOW/AliFlowCommon/AliFlowCumuConstants.cxx deleted file mode 100644 index c1d84587196..00000000000 --- a/PWG2/FLOW/AliFlowCommon/AliFlowCumuConstants.cxx +++ /dev/null @@ -1,61 +0,0 @@ -/************************************************************************** - * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * - * * - * Author: The ALICE Off-line Project. * - * Contributors are mentioned in the code where appropriate. * - * * - * Permission to use, copy, modify and distribute this software and its * - * documentation strictly for non-commercial purposes is hereby granted * - * without fee, provided that the above copyright notice appears in all * - * copies and that both the copyright notice and this permission notice * - * appear in the supporting documentation. The authors make no claims * - * about the suitability of this software for any purpose. It is * - * provided "as is" without express or implied warranty. * - **************************************************************************/ - -#include -#include "AliFlowCumuConstants.h" - -ClassImp(AliFlowCumuConstants) -// Description: constants for the Cumulant flow analysis - -AliFlowCumuConstants* AliFlowCumuConstants::fgPMasterConfig = NULL; - -//______________________________________________________________________________ -AliFlowCumuConstants::AliFlowCumuConstants(): - TNamed(), - fQmax(11), - fPmax(5), - fQmax4(5), - fPmax4(2), - fQmax6(7), - fPmax6(3), - fQmax8(9), - fPmax8(4), - fQmax16(17), - fPmax16(8), - fFlow(2), - fMltpl(1), - fBinWidth(0.1), - fR0(2.2), - fPtMax(3.1), - fPtMin(0.0), - fOtherEquations(kFALSE) -{ - //def ctor -} - -//______________________________________________________________________________ -AliFlowCumuConstants::~AliFlowCumuConstants() -{ - //dtor -} - -//______________________________________________________________________________ -AliFlowCumuConstants* AliFlowCumuConstants::GetMaster() -{ - //get global master config - if (!fgPMasterConfig) fgPMasterConfig = new AliFlowCumuConstants(); - return fgPMasterConfig; -} - diff --git a/PWG2/FLOW/AliFlowCommon/AliFlowCumuConstants.h b/PWG2/FLOW/AliFlowCommon/AliFlowCumuConstants.h deleted file mode 100644 index 63c110fa108..00000000000 --- a/PWG2/FLOW/AliFlowCommon/AliFlowCumuConstants.h +++ /dev/null @@ -1,111 +0,0 @@ -/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * - * See cxx source for full Copyright notice */ - -/* $Id$ */ - -#ifndef ALIFLOWCUMUCONSTANTS_H -#define ALIFLOWCUMUCONSTANTS_H - -#include - -// Description: constants for the LYZ flow makers -// Author: Ante Bilandzic, Nikhef -// modified: Mikolaj Krzewicki, Nikhef (mikolaj.krzewicki@cern.ch) - -#include - -class AliFlowCumuConstants : public TNamed { - -public: - AliFlowCumuConstants(); - virtual ~AliFlowCumuConstants(); - static AliFlowCumuConstants* GetMaster(); - - //static consts for compatibility only! - static const Int_t fgkQmax = 11; // Qmax - static const Int_t fgkPmax = 5; // Pmax - static const Int_t fgkQmax4 = 5; // Qmax4 - static const Int_t fgkPmax4 = 2; // Pmax4 - static const Int_t fgkQmax6 = 7; // Qmax6 - static const Int_t fgkPmax6 = 3; // Pmax6 - static const Int_t fgkQmax8 = 9; // Qmax8 - static const Int_t fgkPmax8 = 4; // Pmax8 - static const Int_t fgkQmax16 = 17; // Qmax16 - static const Int_t fgkPmax16 = 8; // Pmax16 - static const Int_t fgkFlow = 2; // flow estimate - static const Int_t fgkMltpl = 1; // multiple - //static const Double_t fgkBinWidth = 0.1; // - //static const Double_t fgkR0 = 2.2; // - //static const Double_t fgkPtMax = 3.1; // - //static const Double_t fgkPtMin = 0.0; // - //static const Bool_t fgkOtherEquations = kFALSE; // - - Int_t GetQmax() const {return fQmax;} - Int_t GetPmax() const {return fPmax;} - Int_t GetQmax4() const {return fQmax4;} - Int_t GetPmax4() const {return fPmax4;} - Int_t GetQmax6() const {return fQmax6;} - Int_t GetPmax6() const {return fPmax6;} - Int_t GetQmax8() const {return fQmax8;} - Int_t GetPmax8() const {return fPmax8;} - Int_t GetQmax16() const {return fQmax16;} - Int_t GetPmax16() const {return fPmax16;} - Int_t GetFlow() const {return fFlow;} - Int_t GetMltpl() const {return fMltpl;} - Double_t GetBinWidth() const {return fBinWidth;} - Double_t GetR0() const {return fR0;} - Double_t GetPtMax() const {return fPtMax;} - Double_t GetPtMin() const {return fPtMin;} - Bool_t GetOtherEquations() const {return fOtherEquations;} - - void SetQmax( Int_t i ) {fQmax = i;} - void SetPmax( Int_t i ) {fPmax = i;} - void SetQmax4( Int_t i ) {fQmax4 = i;} - void SetPmax4( Int_t i ) {fPmax4 = i;} - void SetQmax6( Int_t i ) {fQmax6 = i;} - void SetPmax6( Int_t i ) {fPmax6 = i;} - void SetQmax8( Int_t i ) {fQmax8 = i;} - void SetPmax8( Int_t i ) {fPmax8 = i;} - void SetQmax16( Int_t i ) {fQmax16 = i;} - void SetPmax16( Int_t i ) {fPmax16 = i;} - void SetFlow( Int_t i ) {fFlow = i;} - void SetMltpl( Int_t i ) {fMltpl = i;} - void SetBinWidth( Double_t i ) {fBinWidth = i;} - void SetR0( Double_t i ) {fR0 = i;} - void SetPtMax( Double_t i ) {fPtMax = i;} - void SetPtMin( Double_t i ) {fPtMin = i;} - void SetOtherEquations( Bool_t i ) {fOtherEquations = i;} - -private: - AliFlowCumuConstants& operator= (const AliFlowCumuConstants& c); - AliFlowCumuConstants(const AliFlowCumuConstants& a); - - Int_t fQmax; // Qmax - Int_t fPmax; // Pmax - Int_t fQmax4; // Qmax4 - Int_t fPmax4; // Pmax4 - Int_t fQmax6; // Qmax6 - Int_t fPmax6; // Pmax6 - Int_t fQmax8; // Qmax8 - Int_t fPmax8; // Pmax8 - Int_t fQmax16; // Qmax16 - Int_t fPmax16; // Pmax16 - Int_t fFlow; // flow estimate - Int_t fMltpl; // multiple - - // Histograms limits - Double_t fBinWidth; // bin width - Double_t fR0; // R0 - Double_t fPtMax; // Pt max - Double_t fPtMin; // Pt min - - // Other numerical equations for cumulants - Bool_t fOtherEquations; // othe equations - - static AliFlowCumuConstants* fgPMasterConfig; // master config - - ClassDef(AliFlowCumuConstants,1) // macro for rootcint -}; - -#endif - diff --git a/PWG2/PWG2flowCommonLinkDef.h b/PWG2/PWG2flowCommonLinkDef.h index 699d8a5cf90..cf84b45d345 100644 --- a/PWG2/PWG2flowCommonLinkDef.h +++ b/PWG2/PWG2flowCommonLinkDef.h @@ -6,7 +6,6 @@ #pragma link C++ namespace AliFlowCommonConstants; #pragma link C++ namespace AliFlowLYZConstants; -#pragma link C++ namespace AliFlowCumuConstants; #pragma link C++ class AliFlowVector+; #pragma link C++ class AliFlowTrackSimple+; @@ -21,7 +20,6 @@ #pragma link C++ class AliFlowLYZHist1+; #pragma link C++ class AliFlowLYZHist2+; -#pragma link C++ class AliCumulantsFunctions+; #pragma link C++ class AliFlowLYZEventPlane+; #pragma link C++ class AliFlowAnalysisWithMCEventPlane+; diff --git a/PWG2/libPWG2flowCommon.pkg b/PWG2/libPWG2flowCommon.pkg index 655be1f25a9..9fc49138432 100644 --- a/PWG2/libPWG2flowCommon.pkg +++ b/PWG2/libPWG2flowCommon.pkg @@ -6,7 +6,6 @@ SRCS= FLOW/AliFlowCommon/AliFlowEventSimple.cxx \ FLOW/AliFlowCommon/AliFlowVector.cxx \ FLOW/AliFlowCommon/AliFlowCommonConstants.cxx \ FLOW/AliFlowCommon/AliFlowLYZConstants.cxx \ - FLOW/AliFlowCommon/AliFlowCumuConstants.cxx \ FLOW/AliFlowCommon/AliFlowEventSimpleMakerOnTheFly.cxx \ FLOW/AliFlowCommon/AliFlowCommonHist.cxx \ FLOW/AliFlowCommon/AliFlowCommonHistResults.cxx \ @@ -19,7 +18,6 @@ SRCS= FLOW/AliFlowCommon/AliFlowEventSimple.cxx \ FLOW/AliFlowCommon/AliFlowAnalysisWithLeeYangZeros.cxx \ FLOW/AliFlowCommon/AliFlowAnalysisWithCumulants.cxx \ FLOW/AliFlowCommon/AliFlowAnalysisWithQCumulants.cxx \ - FLOW/AliFlowCommon/AliCumulantsFunctions.cxx \ FLOW/AliFlowCommon/AliFlowAnalysisWithFittingQDistribution.cxx \ FLOW/AliFlowCommon/AliFlowAnalysisWithMixedHarmonics.cxx \ FLOW/AliFlowCommon/AliFlowAnalysisWithNestedLoops.cxx -- 2.39.3