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