PWG2/FLOW/AliFlowAnalysisWithCumulants.cxx

   1 #define AliFlowAnalysisWithCumulants_cxx
   2
   3 #include "Riostream.h"
   4 #include "AliFlowCommonConstants.h"
   5 #include "AliFlowCommonHist.h"
   6 #include "AliFlowCommonHistResults.h"
   7 #include "TChain.h"
   8 #include "TFile.h"
   9 #include "TParticle.h"
  10
  11 #include "AliFlowEventSimple.h"
  12 #include "AliFlowTrackSimple.h"
  13 #include "AliFlowAnalysisWithCumulants.h"
  14 #include "AliFlowCumuConstants.h"
  15
  16 class TH1;
  17 class TGraph;
  18 class TPave;
  19 class TLatex;
  20 class TMarker;
  21 class TRandom3;
  22 class TObjArray;
  23 class TList;
  24 class TCanvas;
  25 class TSystem;
  26 class TROOT;
  27
  28 class AliFlowVector;
  29 class TVector;
  30
  31 //*******************************
  32 // flow analysis with cumulants *
  33 // author: Ante Bilandzic       *
  34 // email: anteb@nikhef.nl       *
  35 //*******************************
  36
  37 ClassImp(AliFlowAnalysisWithCumulants)
  38
  39 //________________________________________________________________________
  40
  41 AliFlowAnalysisWithCumulants::AliFlowAnalysisWithCumulants():
  42   fEvent(0),
  43   fTrack(0),
  44   fnEvts(0),
  45   fnPrim(0),
  46   fAvM(0),
  47   fR0(0),
  48   fPtMax(0),
  49   fPtMin(0),
  50   fBinWidth(0),
  51   fAvQx(0),
  52   fAvQy(0),
  53   fAvQ2x(0),
  54   fAvQ2y(0),
  55   fHistFileName(0),
  56   fHistFile(0),
  57   fCommonHists(0),
  58   fCommonHistsRes2(0),
  59   fCommonHistsRes4(0),
  60   fCommonHistsRes6(0),
  61   fCommonHistsRes8(0)
  62   {
  63    //constructor
  64    fR0=AliFlowCumuConstants::fgR0;
  65    fPtMax=AliFlowCommonConstants::GetPtMax();
  66    fPtMin=AliFlowCommonConstants::GetPtMin();
  67    fBinWidth=(fPtMax-fPtMin)/fgknBins;
  68
  69    for(Int_t n=0;n<fgknBins;n++){
  70     fBinEventEntries[n]=0;
  71     fBinNoOfParticles[n]=0;
  72     fBinMeanPt[n]=0;
  73     for(Int_t p=0;p<fgkPmax;p++){
  74      for(Int_t q=0;q<fgkQmax;q++){
  75       fAvG[p][q]=0;
  76       fBinEventDRe[n][p][q]=0;
  77       fBinEventDIm[n][p][q]=0;
  78      }
  79     }
  80    }
  81   }
  82
  83 //________________________________________________________________________
  84
  85 AliFlowAnalysisWithCumulants::AliFlowAnalysisWithCumulants(const AliFlowAnalysisWithCumulants&):
  86   fEvent(0),
  87   fTrack(0),
  88   fnEvts(0),
  89   fnPrim(0),
  90   fAvM(0),
  91   fR0(0),
  92   fPtMax(0),
  93   fPtMin(0),
  94   fBinWidth(0),
  95   fAvQx(0),
  96   fAvQy(0),
  97   fAvQ2x(0),
  98   fAvQ2y(0),
  99   fHistFileName(0),
 100   fHistFile(0),
 101   fCommonHists(0),
 102   fCommonHistsRes2(0),
 103   fCommonHistsRes4(0),
 104   fCommonHistsRes6(0),
 105   fCommonHistsRes8(0)
 106   {
 107   //copy constructor
 108    fR0=AliFlowCumuConstants::fgR0;
 109    fPtMax=AliFlowCommonConstants::GetPtMax();
 110    fPtMin=AliFlowCommonConstants::GetPtMin();
 111    fBinWidth=(fPtMax-fPtMin)/fgknBins;
 112
 113    for(Int_t n=0;n<fgknBins;n++){
 114     fBinEventEntries[n]=0;
 115     fBinNoOfParticles[n]=0;
 116     fBinMeanPt[n]=0;
 117     for(Int_t p=0;p<fgkPmax;p++){
 118      for(Int_t q=0;q<fgkQmax;q++){
 119       fAvG[p][q]=0;
 120       fBinEventDRe[n][p][q]=0;
 121       fBinEventDIm[n][p][q]=0;
 122      }
 123     }
 124    }
 125   }
 126
 127 AliFlowAnalysisWithCumulants::~AliFlowAnalysisWithCumulants(){
 128 //desctructor
 129 }
 130
 131 AliFlowAnalysisWithCumulants& AliFlowAnalysisWithCumulants::operator=(const AliFlowAnalysisWithCumulants&)
 132 {
 133  return *this;
 134 }
 135
 136 //___________________________________________________________________________
 137 void AliFlowAnalysisWithCumulants::CreateOutputObjects(){
 138  //output histograms
 139  TString fHistFileName = "cumulants.root";
 140  fHistFile = new TFile(fHistFileName.Data(),"RECREATE");
 141  fCommonHists = new AliFlowCommonHist("Cumulants");//control histograms
 142  fCommonHistsRes2 = new AliFlowCommonHistResults("Cumulants2");
 143  fCommonHistsRes4 = new AliFlowCommonHistResults("Cumulants4");
 144  fCommonHistsRes6 = new AliFlowCommonHistResults("Cumulants6");
 145  fCommonHistsRes8 = new AliFlowCommonHistResults("Cumulants8");
 146 }
 147
 148 //________________________________________________________________________
 149 void AliFlowAnalysisWithCumulants::Exec(AliFlowEventSimple* fEvent) {
 150   //running over data
 151
 152   fCommonHists->FillControlHistograms(fEvent);
 153
 154   Double_t fG[fgkPmax][fgkQmax];//generating function for integrated flow
 155   for(Int_t p=0;p<fgkPmax;p++){
 156    for(Int_t q=0;q<fgkQmax;q++){
 157     fG[p][q]=1.;
 158    }
 159   }
 160
 161   //---------------------------------------------------------
 162   //Q vector stuff
 163   AliFlowVector fQVector;
 164   fQVector.Set(0.,0.);
 165   fQVector.SetMult(0);
 166
 167   fQVector=fEvent->GetQ();//get the Q vector for this event
 168
 169   fAvQx+=fQVector.X();
 170   fAvQy+=fQVector.Y();
 171   fAvQ2x+=pow(fQVector.X(),2.);
 172   fAvQ2y+=pow(fQVector.Y(),2.);
 173   //----------------------------------------------------------
 174
 175   Int_t fnPrim = fEvent->NumberOfTracks();
 176   Int_t fEventNSelTracksIntFlow = fEvent->GetEventNSelTracksIntFlow();
 177   Int_t fSelTracksIntFlow = 0;
 178
 179   cout<<"Number of input tracks for cumulant analysis: "<<fnPrim<<endl;
 180   cout<<"Number of selected tracks for cumulant analysis: "<<fEventNSelTracksIntFlow<<endl;
 181
 182     //------------------------------------------------------------------------------------
 183     //STARTING THE FIRST LOOP (CALCULATING THE GENERATING FUNCTION FOR INTEGRATED FLOW)
 184     for(Int_t i=0;i<fnPrim;i++){
 185     fTrack=fEvent->GetTrack(i);
 186      if(fTrack&&fTrack->UseForIntegratedFlow()){
 187       fSelTracksIntFlow++;
 188       for(Int_t p=0;p<fgkPmax;p++){
 189        for(Int_t q=0;q<fgkQmax;q++){
 190         fG[p][q]*=(1.+(2.*fR0*sqrt(p+1.)/fEventNSelTracksIntFlow)*cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax));
 191        }
 192       }
 193      }
 194     }
 195     // ENDING THE FIRST LOOP OVER TRACKS
 196     //------------------------------------------------------------------------------------
 197
 198
 199   //------------------------------------------------------------------------------------
 200   //avarage multiplicity
 201   fAvM+=fSelTracksIntFlow;
 202   //avarage of the generating function for integrated flow
 203   for(Int_t p=0;p<fgkPmax;p++){
 204    for(Int_t q=0;q<fgkQmax;q++){
 205     fAvG[p][q]+=1.*fG[p][q];
 206    }
 207   }
 208   //------------------------------------------------------------------------------------
 209
 210
 211   //STARTING WITH DIFFERENTIAL FLOW...
 212   Int_t fBinEntries[fgknBins]={0};//stores the number of particles per bin for the current event
 213   Double_t fNumDRe[fgknBins][fgkPmax][fgkQmax]={0.};//real part of the numerator of D (see relation (10) in PG)
 214   Double_t fNumDIm[fgknBins][fgkPmax][fgkQmax]={0.};//imaginary part of the numerator D
 215
 216   //------------------------------------------------------------------------------------------------
 217   //STARTING THE SECOND LOOP OVER TRACKS (CALCULATING THE GENERATING FUNCTION FOR DIFFERENTIAL FLOW)
 218   for(Int_t i=0;i<fnPrim;i++){
 219     fTrack=fEvent->GetTrack(i);
 220     if (fTrack && fTrack->UseForDifferentialFlow()) {
 221       Int_t fBin=TMath::Nint(floor(fTrack->Pt()/fBinWidth));
 222       if(fBin>=fgknBins)continue;//ignoring the particles with pt>ptMax
 223       fBinNoOfParticles[fBin]++;
 224       fBinEntries[fBin]++;
 225       fBinMeanPt[fBin]+=fTrack->Pt();
 226       for(Int_t p=0;p<fgkPmax;p++){
 227         for(Int_t q=0;q<fgkQmax;q++){
 228           fNumDRe[fBin][p][q]+=fG[p][q]*cos(fgkMltpl*fgkFlow*fTrack->Phi())/
 229             (1.+(2.*fR0*sqrt(p+1.)/fSelTracksIntFlow) *
 230              cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax));
 231           fNumDIm[fBin][p][q]+=fG[p][q]*sin(fgkMltpl*fgkFlow*fTrack->Phi())/
 232             (1.+(2.*fR0*sqrt(p+1.)/fSelTracksIntFlow) *
 233              cos(fgkFlow*fTrack->Phi()-2.*q*TMath::Pi()/fgkQmax));
 234         }
 235       }
 236     }
 237   }
 238   //ENDING THE SECOND LOOP OVER TRACKS
 239   //-----------------------------------------------------------------------------------------------
 240
 241
 242
 243   //----------------------------------------------------------
 244   //AVARAGING OVER ALL pt BINS WITHIN ONE EVENT
 245   for(Int_t b=0;b<fgknBins;b++){
 246     if(fBinEntries[b]==0)continue;
 247     fBinEventEntries[b]++;
 248     for(Int_t p=0;p<fgkPmax;p++){
 249       for(Int_t q=0;q<fgkQmax;q++){
 250         fBinEventDRe[b][p][q]+=fNumDRe[b][p][q]/fBinEntries[b];
 251         fBinEventDIm[b][p][q]+=fNumDIm[b][p][q]/fBinEntries[b];
 252       }
 253     }
 254   }
 255   //----------------------------------------------------------
 256 }
 257
 258 //________________________________________________________________________
 259 void AliFlowAnalysisWithCumulants::Terminate(Int_t fCount){
 260   //final results
 261   cout<<""<<endl;
 262   cout<<"***************************************"<<endl;
 263   cout<<"**** results of cumulant analysis: ****"<<endl;
 264   cout<<"***************************************"<<endl;
 265   cout<<""<<endl;
 266   cout<<"nEvts = "<<fCount<<endl;
 267
 268   Int_t fnEvts=fCount;
 269
 270   //final avarage multiplicity
 271   fAvM/=fnEvts;
 272
 273   //final avarage of generating function for the integrated flow
 274   for(Int_t p=0;p<fgkPmax;p++){
 275    for(Int_t q=0;q<fgkQmax;q++){
 276     fAvG[p][q]/=fnEvts;
 277    }
 278   }
 279
 280   //final avarage of the Q vector stuff
 281   fAvQx/=fnEvts;
 282   fAvQy/=fnEvts;
 283   fAvQ2x/=fnEvts;
 284   fAvQ2y/=fnEvts;
 285
 286   /////////////////////////////////////////////////////////////////////////////
 287   //////////////////gen. function for the cumulants////////////////////////////
 288   /////////////////////////////////////////////////////////////////////////////
 289
 290   Double_t fC[fgkPmax][fgkQmax]={0.};
 291   for (Int_t p=0;p<fgkPmax;p++){
 292     for (Int_t q=0;q<fgkQmax;q++){
 293       fC[p][q]=1.*fAvM*(pow(fAvG[p][q],(1./fAvM))-1.);
 294     }
 295   }
 296
 297   /////////////////////////////////////////////////////////////////////////////
 298   ///////avaraging the gen. function for the cumulants over azimuth////////////
 299   /////////////////////////////////////////////////////////////////////////////
 300
 301   Double_t fCAv[fgkPmax]={0.};
 302   for (Int_t p=0;p<fgkPmax;p++){
 303     Double_t fTempHere=0.;
 304     for (Int_t q=0;q<fgkQmax;q++){
 305       fTempHere+=1.*fC[p][q];
 306     }
 307     fCAv[p]=1.*fTempHere/fgkQmax;
 308   }
 309
 310   /////////////////////////////////////////////////////////////////////////////
 311   //////////////////////////////////final results//////////////////////////////
 312   /////////////////////////////////////////////////////////////////////////////
 313
 314   Double_t fCumulant[fgkPmax];//c_{iFlow}\{2(p+1)\}
 315
 316   fCumulant[0]=(1./(fR0*fR0)) * (8.*fCAv[0] - 14.*fCAv[1] + (56./3.)*fCAv[2] - (35./2.)*fCAv[3] +
 317                               (56./5.)*fCAv[4] - (14./3.)*fCAv[5] + (8./7.)*fCAv[6] - (1./8.)*fCAv[7]);
 318
 319   fCumulant[1]=(1./pow(fR0,4.)) * ((-1924./35.)*fCAv[0] + (621./5.)*fCAv[1] - (8012./45.)*fCAv[2] +
 320                                  (691./4.)*fCAv[3] - (564./5.)*fCAv[4] + (2143./45.)*fCAv[5] -
 321                                  (412./35.)*fCAv[6] + (363./280.)*fCAv[7]);
 322
 323   fCumulant[2]=(1./pow(fR0,6.)) * (349.*fCAv[0] - (18353./20.)*fCAv[1] + (7173./5.)*fCAv[2] -
 324                                  1457.*fCAv[3] + (4891./5.)*fCAv[4] - (1683./4.)*fCAv[5] +
 325                                  (527./5.)*fCAv[6] - (469./40.)*fCAv[7]);
 326
 327   fCumulant[3]=(1./pow(fR0,8.)) * ((-10528./5.)*fCAv[0] + (30578./5.)*fCAv[1] - (51456./5.)*fCAv[2] +
 328                                  10993.*fCAv[3] - (38176./5.)*fCAv[4] + (16818./5.)*fCAv[5] -
 329                                  (4288./5.)*fCAv[6] + (967./10.)*fCAv[7]);
 330
 331   fCumulant[4]=(1./pow(fR0,10.)) * (11500.*fCAv[0] - 35800.*fCAv[1] + 63900.*fCAv[2] - 71600.*fCAv[3] +
 332                                   51620.*fCAv[4] - 23400.*fCAv[5] + 6100.*fCAv[6] - 700.*fCAv[7]);
 333
 334   fCumulant[5]=(1./pow(fR0,12.)) * (-52560.*fCAv[0] + 172080.*fCAv[1] - 321840.*fCAv[2] + 376200.*fCAv[3] -
 335                                   281520.*fCAv[4] + 131760.*fCAv[5] - 35280.*fCAv[6] + 4140.*fCAv[7]);
 336
 337   fCumulant[6]=(1./pow(fR0,14.)) * (176400.*fCAv[0] - 599760.*fCAv[1] + 1164240.*fCAv[2] - 1411200.*fCAv[3] +
 338                                   1093680.*fCAv[4] - 529200.*fCAv[5] + 146160.*fCAv[6] - 17640.*fCAv[7]);
 339
 340   fCumulant[7]=(1./pow(fR0,16.)) * (-322560*fCAv[0] + 1128960.*fCAv[1] - 2257920.*fCAv[2] + 2822400.*fCAv[3] -
 341                                   2257920.*fCAv[4] + 1128960.*fCAv[5] - 322560.*fCAv[6] + 40320.*fCAv[7]);
 342
 343   cout<<""<<endl;
 344   cout<<"***************************"<<endl;
 345   cout<<"cumulants:"<<endl;
 346
 347   cout<<" c_"<<fgkFlow<<"{2} = "<<fCumulant[0]<<endl;
 348   cout<<" c_"<<fgkFlow<<"{4} = "<<fCumulant[1]<<endl;
 349   cout<<" c_"<<fgkFlow<<"{6} = "<<fCumulant[2]<<endl;
 350   cout<<" c_"<<fgkFlow<<"{8} = "<<fCumulant[3]<<endl;
 351   cout<<"c_"<<fgkFlow<<"{10} = "<<fCumulant[4]<<endl;
 352   cout<<"c_"<<fgkFlow<<"{12} = "<<fCumulant[5]<<endl;
 353   cout<<"c_"<<fgkFlow<<"{14} = "<<fCumulant[6]<<endl;
 354   cout<<"c_"<<fgkFlow<<"{16} = "<<fCumulant[7]<<endl;
 355   cout<<""<<endl;
 356   cout<<"integrated flow: "<<endl;
 357
 358   Double_t fV2=0.,fV4=0.,fV6=0.,fV8=0.;
 359   Double_t fSdQ[4]={0.};
 360   Double_t fChiQ[4]={0.};
 361    if (fCumulant[0]>=0.){
 362     fV2=sqrt(fCumulant[0]);
 363     fChiQ[0]=fAvM*fV2/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV2*fAvM,2.),0.5);
 364     fSdQ[0]=pow(((1./(2.*fAvM*fnEvts))*((1.+1.*pow(fChiQ[0],2))/(1.*pow(fChiQ[0],2)))),0.5);
 365     cout<<" v_"<<fgkFlow<<"{2} = "<<100.*fV2<<"%, chi{2} = "<<fChiQ[0]<<", sd{2} = "<<100.*fSdQ[0]<<"%"<<endl;
 366     fCommonHistsRes2->FillIntegratedFlow(100.*fV2,100.*fSdQ[0]);
 367     fCommonHistsRes2->FillChi(fChiQ[0]);
 368    } else {
 369     cout<<" v_"<<fgkFlow<<"{2} = Im"<<endl;
 370   }
 371   if (fCumulant[1]<=0.){
 372     fV4=pow(-fCumulant[1],(1./4.));
 373     fChiQ[1]=fAvM*fV4/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV4*fAvM,2.),0.5);
 374     fSdQ[1]=(1./(pow(2.*fAvM*fnEvts,.5)))*pow((1.+2.*pow(fChiQ[1],2)+(1./4.)*pow(fChiQ[1],4.)+(1./4.)*pow(fChiQ[1],6.))/((1./4.)*pow(fChiQ[1],6.)),.5);
 375     cout<<" v_"<<fgkFlow<<"{4} = "<<100.*fV4<<"%, chi{4} = "<<fChiQ[1]<<", sd{4} = "<<100.*fSdQ[1]<<"%"<<endl;
 376     fCommonHistsRes4->FillIntegratedFlow(100.*fV4,100.*fSdQ[1]);
 377     fCommonHistsRes4->FillChi(fChiQ[1]);
 378   } else {
 379     cout<<" v_"<<fgkFlow<<"{4} = Im"<<endl;
 380   }
 381   if (fCumulant[2]>=0.){
 382     fV6=pow((1./4.)*fCumulant[2],(1./6.));
 383     fChiQ[2]=fAvM*fV6/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV6*fAvM,2.),0.5);
 384     fSdQ[2]=(1./(pow(2.*fAvM*fnEvts,.5)))*pow((3.+18.*pow(fChiQ[2],2)+9.*pow(fChiQ[2],4.)+28.*pow(fChiQ[2],6.)+12.*pow(fChiQ[2],8.)+24.*pow(fChiQ[2],10.))/(24.*pow(fChiQ[2],10.)),.5);
 385     cout<<" v_"<<fgkFlow<<"{6} = "<<100.*fV6<<"%, chi{6} = "<<fChiQ[2]<<", sd{6} = "<<100.*fSdQ[2]<<"%"<<endl;
 386     fCommonHistsRes6->FillIntegratedFlow(100.*fV6,100.*fSdQ[2]);
 387     fCommonHistsRes6->FillChi(fChiQ[2]);
 388   } else {
 389     cout<<" v_"<<fgkFlow<<"{6} = Im"<<endl;
 390   }
 391   if (fCumulant[3]<=0.){
 392     fV8=pow(-(1./33.)*fCumulant[3],(1./8.));
 393     fChiQ[3]=fAvM*fV8/pow(fAvQ2x+fAvQ2y-pow(fAvQx,2.)-pow(fAvQy,2.)-pow(fV8*fAvM,2.),0.5);
 394     fSdQ[3]=(1./(pow(2.*fAvM*fnEvts,.5)))*pow((12.+96.*pow(fChiQ[3],2)+72.*pow(fChiQ[3],4.)+304.*pow(fChiQ[3],6.)+257.*pow(fChiQ[3],8.)+804.*pow(fChiQ[3],10.)+363.*pow(fChiQ[3],12.)+726.*pow(fChiQ[3],14.))/(726.*pow(fChiQ[3],14.)),.5);
 395     cout<<" v_"<<fgkFlow<<"{8} = "<<100.*fV8<<"%, chi{8} = "<<fChiQ[3]<<", sd{8} = "<<100.*fSdQ[3]<<"%"<<endl;
 396      fCommonHistsRes8->FillIntegratedFlow(100.*fV8,100.*fSdQ[3]);
 397      fCommonHistsRes8->FillChi(fChiQ[3]);
 398   } else {
 399     cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl;
 400
 401   }
 402   if (fCumulant[4]>=0.){
 403     cout<<"v_"<<fgkFlow<<"{10} = "<<100.*pow((1./456.)*fCumulant[4],(1./10.))<<"%"<<endl;
 404   } else {
 405     cout<<"v_"<<fgkFlow<<"{10} = Im"<<endl;
 406   }
 407   if (fCumulant[5]<=0.){
 408     cout<<"v_"<<fgkFlow<<"{12} = "<<100.*pow(-(1./9460.)*fCumulant[5],(1./12.))<<"%"<<endl;
 409   } else {
 410     cout<<"v_"<<fgkFlow<<"{12} = Im"<<endl;
 411   }
 412   if (fCumulant[6]>=0.){
 413     cout<<"v_"<<fgkFlow<<"{14} = "<<100.*pow((1./274800.)*fCumulant[6],(1./14.))<<"%"<<endl;
 414   } else {
 415     cout<<"v_"<<fgkFlow<<"{14} = Im"<<endl;
 416   }
 417   if (fCumulant[7]<=0.){
 418     cout<<"v_"<<fgkFlow<<"{16} = "<<100.*pow(-(1./10643745.)*fCumulant[7],(1./16.))<<"%"<<endl;
 419   } else {
 420     cout<<"v_"<<fgkFlow<<"{16} = Im"<<endl;
 421   }
 422   cout<<"***************************"<<endl;
 423
 424   cout<<""<<endl;
 425   cout<<"continuing with calculations for differential flow..."<<endl;
 426   cout<<""<<endl;
 427
 428   Double_t fBinEventDReAv[fgknBins][fgkPmax][fgkQmax]={0.};
 429   Double_t fBinEventDImAv[fgknBins][fgkPmax][fgkQmax]={0.};
 430
 431   for(Int_t b=0;b<fgknBins;b++){
 432     if(fBinEventEntries[b]==0) continue;
 433     for(Int_t p=0;p<fgkPmax;p++){
 434       for(Int_t q=0;q<fgkQmax;q++){
 435         fBinEventDReAv[b][p][q]=fBinEventDRe[b][p][q]/fBinEventEntries[b];//avarage of the real part of numerator in relation (10) in PG
 436         fBinEventDImAv[b][p][q]=fBinEventDIm[b][p][q]/fBinEventEntries[b];//avarage of the imaginary part of numerator in relation (10) in PG
 437       }
 438     }
 439   }
 440
 441   Double_t fX[fgknBins][fgkPmax][fgkQmax]={0.};//see the text bellow relation (11) in PG
 442   Double_t fY[fgknBins][fgkPmax][fgkQmax]={0.};
 443
 444   for(Int_t b=0;b<fgknBins;b++){
 445     for(Int_t p=0;p<fgkPmax;p++){
 446       for(Int_t q=0;q<fgkQmax;q++){
 447         fX[b][p][q]=fBinEventDReAv[b][p][q]/fAvG[p][q];
 448         fY[b][p][q]=fBinEventDImAv[b][p][q]/fAvG[p][q];
 449       }
 450     }
 451   }
 452   cout<<""<<endl;
 453   cout<<"I have calculated X and Y."<<endl;
 454   cout<<""<<endl;
 455
 456   Double_t fD[fgknBins][fgkPmax]={0.};//implementing relation (11) from PG
 457
 458   for (Int_t b=0;b<fgknBins;b++){
 459     for (Int_t p=0;p<fgkPmax;p++){
 460       Double_t fTempHere3=0.;
 461       for (Int_t q=0;q<fgkQmax;q++){
 462         fTempHere3+=cos(fgkMltpl*2.*q*TMath::Pi()/fgkQmax)*fX[b][p][q] + sin(fgkMltpl*2.*q*TMath::Pi()/fgkQmax)*fY[b][p][q];
 463       }
 464       fD[b][p]=1.*(pow(fR0*pow(p+1.,.5),fgkMltpl)/fgkQmax)*fTempHere3;
 465     }
 466   }
 467
 468   cout<<""<<endl;
 469   cout<<"calculating differential cumulants now..."<<endl;
 470   cout<<""<<endl;
 471
 472   Double_t fDiffCumulant2[fgknBins]={0.};//implementing relation (12) from PG
 473   Double_t fDiffCumulant4[fgknBins]={0.};
 474   Double_t fDiffCumulant6[fgknBins]={0.};
 475   Double_t fDiffCumulant8[fgknBins]={0.};
 476
 477   for (Int_t b=0;b<fgknBins;b++){
 478     fDiffCumulant2[b]=(1./(fR0*fR0))*(4.*fD[b][0]-3.*fD[b][1]+(4./3.)*fD[b][2]-(1./4.)*fD[b][3]);
 479     fDiffCumulant4[b]=(1./pow(fR0,4.))*((-26./3.)*fD[b][0]+(19./2.)*fD[b][1]-(14./3.)*fD[b][2]+(11./12.)*fD[b][3]);
 480     fDiffCumulant6[b]=(1./pow(fR0,6.))*(18.*fD[b][0]-24.*fD[b][1]+14.*fD[b][2]-3.*fD[b][3]);
 481     fDiffCumulant8[b]=(1./pow(fR0,8.))*(-24.*fD[b][0]+36.*fD[b][1]-24.*fD[b][2]+6.*fD[b][3]);
 482   }
 483
 484   Double_t fv2[fgknBins],fv4[fgknBins],fv6[fgknBins],fv8[fgknBins];
 485   Double_t fAvPt[fgknBins];
 486   Double_t fSddiff2[fgknBins],fSddiff4[fgknBins];
 487
 488   cout<<"number of pt bins: "<<fgknBins<<endl;
 489   cout<<"****************************************"<<endl;
 490   for (Int_t b=0;b<fgknBins;b++){
 491     if(fBinNoOfParticles[b]==0)continue;
 492     fAvPt[b]=fBinMeanPt[b]/fBinNoOfParticles[b];
 493     cout<<"pt bin: "<<b*fBinWidth<<"-"<<(b+1)*fBinWidth<<" GeV"<<endl;
 494     cout<<"number of particles in this pt bin: "<<fBinNoOfParticles[b]<<endl;
 495     cout<<"mean pt in this bin: "<<fAvPt[b]<<" GeV"<<endl;
 496     if(fCumulant[0]>=0){
 497       fv2[b]=100.*fDiffCumulant2[b]/pow(fCumulant[0],.5);
 498       fSddiff2[b]=pow((1./(2.*fBinNoOfParticles[b]))*((1.+pow(fChiQ[0],2.))/pow(fChiQ[0],2.)),0.5);
 499       cout<<"v'_2/2{2} = "<<fv2[b]<<"%, "<<" "<<"sd{2} = "<<100.*fSddiff2[b]<<"%"<<endl;
 500       fCommonHistsRes2->FillDifferentialFlow(b+1,fv2[b],100.*fSddiff2[b]);
 501     }else{
 502       cout<<"v'_2/2{2} = Im"<<endl;
 503     }
 504     if(fCumulant[1]<=0){
 505       fv4[b]=-100.*fDiffCumulant4[b]/pow(-fCumulant[1],.75);
 506       fSddiff4[b]=pow((1./(2.*fBinNoOfParticles[b]))*((2.+6.*pow(fChiQ[1],2.)+pow(fChiQ[1],4.)+pow(fChiQ[1],6.))/pow(fChiQ[1],6.)),0.5);
 507       cout<<"v'_2/2{4} = "<<fv4[b]<<"%, "<<" "<<"sd{4} = "<<100.*fSddiff4[b]<<"%"<<endl;
 508       fCommonHistsRes4->FillDifferentialFlow(b+1,fv4[b],100.*fSddiff4[b]);
 509     }else{
 510       cout<<"v'_2/2{4} = Im"<<endl;
 511     }
 512     if(fCumulant[2]>=0){
 513       cout<<"v'_2/2{6} = "<<100.*fDiffCumulant6[b]/(4.*pow((1./4.)*fCumulant[2],(5./6.)))<<"%"<<endl;
 514       fv6[b]=100.*fDiffCumulant6[b]/(4.*pow((1./4.)*fCumulant[2],(5./6.)));
 515       fCommonHistsRes6->FillDifferentialFlow(b+1,fv6[b],0.);
 516     }else{
 517       cout<<"v'_2/2{6} = Im"<<endl;
 518     }
 519     if(fCumulant[3]<=0){
 520       cout<<"v'_2/2{8} = "<<-100.*fDiffCumulant8[b]/(33.*pow(-(1./33.)*fCumulant[3],(7./8.)))<<"%"<<endl;
 521       fv8[b]=-100.*fDiffCumulant8[b]/(33.*pow(-(1./33.)*fCumulant[3],(7./8.)));
 522       fCommonHistsRes8->FillDifferentialFlow(b+1,fv8[b],0.);
 523     }else{
 524       cout<<"v'_2/2{8} = Im"<<endl;
 525     }
 526     cout<<"****************************************"<<endl;
 527   }
 528 }
 529
 530
 531