PWG2/FLOW/AliFlowCommon/AliCumulantsFunctions.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 /**********************************
  17  * functions and equations needed *
  18  * for calculation of cumulants   *
  19  * and final flow estimates       *
  20  *                                *
  21  * author: Ante Bilandzic         *
  22  *          (anteb@nikhef.nl)     *
  23  *********************************/
  24
  25 #define AliCumulantsFunctions_cxx
  26
  27 #include "Riostream.h"
  28 #include "TChain.h"
  29 #include "TFile.h"
  30 #include "TList.h"
  31 #include "TParticle.h"
  32 #include <TMatrixD.h>
  33 #include <TVectorD.h>
  34
  35 #include "TProfile.h"
  36 #include "TProfile2D.h"
  37 #include "TProfile3D.h"
  38 #include "TH1.h"
  39 #include "TH1D.h"
  40
  41 #include "AliFlowEventSimple.h"
  42 #include "AliFlowTrackSimple.h"
  43 #include "AliFlowAnalysisWithCumulants.h"
  44 #include "AliFlowCumuConstants.h"
  45 #include "AliFlowCommonConstants.h"
  46 #include "AliFlowCommonHist.h"
  47 #include "AliFlowCommonHistResults.h"
  48 #include "AliCumulantsFunctions.h"
  49
  50 ClassImp(AliCumulantsFunctions)
  51
  52 //================================================================================================================_
  53
  54 AliCumulantsFunctions::AliCumulantsFunctions():
  55  fIntGenFun(NULL),
  56  fIntGenFun4(NULL),//only for other system of Eq.
  57  fIntGenFun6(NULL),//only for other system of Eq.
  58  fIntGenFun8(NULL),//only for other system of Eq.
  59  fIntGenFun16(NULL),//only for other system of Eq.
  60  fAvMult4(NULL),//only for other system of Eq.
  61  fAvMult6(NULL),//only for other system of Eq.
  62  fAvMult8(NULL),//only for other system of Eq.
  63  fAvMult16(NULL),//only for other system of Eq.
  64  fDiffPtRPGenFunRe(NULL),
  65  fDiffPtRPGenFunIm(NULL),
  66  fPtBinRPNoOfParticles(NULL),
  67  fDiffEtaRPGenFunRe(NULL),
  68  fDiffEtaRPGenFunIm(NULL),
  69  fEtaBinRPNoOfParticles(NULL),
  70  fDiffPtPOIGenFunRe(NULL),
  71  fDiffPtPOIGenFunIm(NULL),
  72  fPtBinPOINoOfParticles(NULL),
  73  fDiffEtaPOIGenFunRe(NULL),
  74  fDiffEtaPOIGenFunIm(NULL),
  75  fEtaBinPOINoOfParticles(NULL),
  76  fifr(NULL),
  77  fdfr2(NULL),
  78  fdfr4(NULL),
  79  fdfr6(NULL),
  80  fdfr8(NULL),
  81  fAvMult(NULL),
  82  fQVector(NULL),
  83  fAverageOfSquaredWeight(NULL),
  84  fchr2nd(NULL),
  85  fchr4th(NULL),
  86  fchr6th(NULL),
  87  fchr8th(NULL),
  88  fch(NULL)//common control histograms
  89  /*
  90  fdRe0(NULL),
  91  fdRe1(NULL),
  92  fdRe2(NULL),
  93  fdRe3(NULL),
  94  fdRe4(NULL),
  95  fdRe5(NULL),
  96  fdRe6(NULL),
  97  fdRe7(NULL),
  98  fdIm0(NULL),
  99  fdIm1(NULL),
 100  fdIm2(NULL),
 101  fdIm3(NULL),
 102  fdIm4(NULL),
 103  fdIm5(NULL),
 104  fdIm6(NULL),
 105  fdIm7(NULL)
 106  */
 107 {
 108  //default constructor
 109 }
 110
 111 AliCumulantsFunctions::~AliCumulantsFunctions()
 112 {
 113  //destructor
 114 }
 115
 116 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):
 117  fIntGenFun(intGenFun),
 118  fIntGenFun4(intGenFun4),//only for other system of Eq.
 119  fIntGenFun6(intGenFun6),//only for other system of Eq.
 120  fIntGenFun8(intGenFun8),//only for other system of Eq.
 121  fIntGenFun16(intGenFun16),//only for other system of Eq.
 122  fAvMult4(avMult4),//only for other system of Eq.
 123  fAvMult6(avMult6),//only for other system of Eq.
 124  fAvMult8(avMult8),//only for other system of Eq.
 125  fAvMult16(avMult16),//only for other system of Eq.
 126  fDiffPtRPGenFunRe(diffPtRPGenFunRe),
 127  fDiffPtRPGenFunIm(diffPtRPGenFunIm),
 128  fPtBinRPNoOfParticles(ptBinRPNoOfParticles),
 129  fDiffEtaRPGenFunRe(diffEtaRPGenFunRe),
 130  fDiffEtaRPGenFunIm(diffEtaRPGenFunIm),
 131  fEtaBinRPNoOfParticles(etaBinRPNoOfParticles),
 132  fDiffPtPOIGenFunRe(diffPtPOIGenFunRe),
 133  fDiffPtPOIGenFunIm(diffPtPOIGenFunIm),
 134  fPtBinPOINoOfParticles(ptBinPOINoOfParticles),
 135  fDiffEtaPOIGenFunRe(diffEtaPOIGenFunRe),
 136  fDiffEtaPOIGenFunIm(diffEtaPOIGenFunIm),
 137  fEtaBinPOINoOfParticles(etaBinPOINoOfParticles),
 138  fifr(ifr),
 139  fdfr2(dfr2),
 140  fdfr4(dfr4),
 141  fdfr6(dfr6),
 142  fdfr8(dfr8),
 143  fAvMult(avMult),
 144  fQVector(qVector),
 145  fAverageOfSquaredWeight(averageOfSquaredWeight),
 146  fchr2nd(chr2nd),
 147  fchr4th(chr4th),
 148  fchr6th(chr6th),
 149  fchr8th(chr8th),
 150  fch(ch)//common control histograms
 151  /*
 152  fdRe0(dRe0),
 153  fdRe1(dRe1),
 154  fdRe2(dRe2),
 155  fdRe3(dRe3),
 156  fdRe4(dRe4),
 157  fdRe5(dRe5),
 158  fdRe6(dRe6),
 159  fdRe7(dRe7),
 160  fdIm0(dIm0),
 161  fdIm1(dIm1),
 162  fdIm2(dIm2),
 163  fdIm3(dIm3),
 164  fdIm4(dIm4),
 165  fdIm5(dIm5),
 166  fdIm6(dIm6),
 167  fdIm7(dIm7)
 168  */
 169 {
 170  //custom constructor
 171 }
 172
 173 //================================================================================================================
 174
 175 void AliCumulantsFunctions::Calculate()
 176 {
 177  //calculate cumulants and final integrated and differential flow estimates and store the results into output histograms
 178  const Int_t cQmax=AliFlowCumuConstants::GetMaster()->GetQmax();     //needed for numerics
 179  const Int_t cPmax=AliFlowCumuConstants::GetMaster()->GetPmax();     //needed for numerics
 180  const Int_t cQmax4=AliFlowCumuConstants::GetMaster()->GetQmax4();   //needed for numerics
 181  const Int_t cPmax4=AliFlowCumuConstants::GetMaster()->GetPmax4();   //needed for numerics
 182  const Int_t cQmax6=AliFlowCumuConstants::GetMaster()->GetQmax6();   //needed for numerics
 183  const Int_t cPmax6=AliFlowCumuConstants::GetMaster()->GetPmax6();   //needed for numerics
 184  const Int_t cQmax8=AliFlowCumuConstants::GetMaster()->GetQmax8();   //needed for numerics
 185  const Int_t cPmax8=AliFlowCumuConstants::GetMaster()->GetPmax8();   //needed for numerics
 186  const Int_t cQmax16=AliFlowCumuConstants::GetMaster()->GetQmax16(); //needed for numerics
 187  const Int_t cPmax16=AliFlowCumuConstants::GetMaster()->GetPmax16(); //needed for numerics
 188
 189  const Int_t cFlow=AliFlowCumuConstants::GetMaster()->GetFlow();   //integrated flow coefficient to be calculated
 190  const Int_t cMltpl=AliFlowCumuConstants::GetMaster()->GetMltpl(); //the multiple in p=m*n (diff. flow)
 191  const Int_t cnBinsPt=100;                        //number of pt bins //to be improved
 192  const Int_t cnBinsEta=80;                        //number of eta bins //to be improved
 193
 194  Double_t fR0=AliFlowCumuConstants::GetMaster()->GetR0();              //needed for numerics
 195  //Double_t fPtMax=AliFlowCommonConstants::GetMaster()->GetPtMax(); //maximum pt
 196  //Double_t fPtMin=AliFlowCommonConstants::GetMaster()->GetPtMin(); //minimum pt
 197  //Double_t fBinWidthPt=(fPtMax-fPtMin)/cnBinsPt;    //width of pt bin (in GeV)
 198
 199  Bool_t fOtherEquations=AliFlowCumuConstants::GetMaster()->GetOtherEquations();     //numerical equations for cumulants solved up to different highest order
 200
 201  //avarage selected multiplicity
 202  Double_t dAvM=0.;
 203  if(fAvMult)
 204  {
 205   dAvM=fAvMult->GetBinContent(1);
 206  }
 207
 208  //number of events
 209  Int_t nEvents=0;
 210  if(fAvMult)
 211  {
 212   nEvents=(Int_t)(fAvMult->GetBinEntries(1));
 213  }
 214
 215  //<Q-vector stuff>
 216  Double_t dAvQx=0.,dAvQy=0.,dAvQ2x=0.,dAvQ2y=0.;
 217  if(fQVector)
 218  {
 219   dAvQx  = fQVector->GetBinContent(1); //<Q_x>
 220   dAvQy  = fQVector->GetBinContent(2); //<Q_y>
 221   dAvQ2x = fQVector->GetBinContent(3); //<(Q_x)^2>
 222   dAvQ2y = fQVector->GetBinContent(4); //<(Q_y)^2>
 223  }
 224
 225  //<w^2>
 226  Double_t dAvw2 = 1.;
 227  if(nEvents>0)
 228  {
 229   dAvw2 = fAverageOfSquaredWeight->GetBinContent(1);
 230   if(dAvw2 == 0)
 231   {
 232    cout<<endl;
 233    cout<<"WARNING: Average squared weight is 0 in GFC. Most probably one of the histograms in <weights.root> is empty. Nothing will be calculated !!!!"<<endl;
 234    cout<<endl;
 235   }
 236  }
 237
 238  //<G[p][q]>
 239  TMatrixD dAvG(cPmax, cQmax);
 240  Bool_t someAvGEntryIsNegative=kFALSE;
 241  if(fIntGenFun)
 242  {
 243   for(Int_t p=0;p<cPmax;p++)
 244   {
 245    for(Int_t q=0;q<cQmax;q++)
 246    {
 247     dAvG(p,q)=fIntGenFun->GetBinContent(p+1,q+1);
 248     if(dAvG(p,q)<0.)
 249     {
 250      someAvGEntryIsNegative=kTRUE;
 251     }
 252    }
 253   }
 254  }
 255
 256  /////////////////////////////////////////////////////////////////////////////
 257  //////////////////gen. function for the cumulants////////////////////////////
 258  /////////////////////////////////////////////////////////////////////////////
 259
 260  TMatrixD dC(cPmax, cQmax);//C[p][q]
 261  if(dAvM>0 && someAvGEntryIsNegative==kFALSE)
 262  {
 263   for(Int_t p=0;p<cPmax;p++)
 264   {
 265    for(Int_t q=0;q<cQmax;q++)
 266    {
 267     dC(p,q)=1.*dAvM*(pow(dAvG(p,q),(1./dAvM))-1.);
 268    }
 269   }
 270  }
 271
 272  /////////////////////////////////////////////////////////////////////////////
 273  ///////avaraging the gen. function for the cumulants over azimuth////////////
 274  /////////////////////////////////////////////////////////////////////////////
 275
 276  TVectorD dAvC(cPmax);//<C[p][q]>
 277  Double_t tempHere=0.;
 278  for(Int_t p=0;p<cPmax;p++)
 279  {
 280   tempHere=0.;
 281   for(Int_t q=0;q<cQmax;q++)
 282   {
 283    tempHere+=1.*dC(p,q);
 284   }
 285   dAvC[p]=1.*tempHere/cQmax;
 286  }
 287
 288  /////////////////////////////////////////////////////////////////////////////
 289  //////////////////////////////////final results//////////////////////////////
 290  /////////////////////////////////////////////////////////////////////////////
 291
 292  TVectorD cumulant(cPmax);//array to store various order cumulants
 293
 294  //system of eq. for the cumulants
 295  cumulant[0] = (-1./(60*fR0*fR0))*((-300.)*dAvC[0]+300.*dAvC[1]-200.*dAvC[2]+75.*dAvC[3]-12.*dAvC[4]);
 296  cumulant[1] = (-1./(6.*pow(fR0,4.)))*(154.*dAvC[0]-214.*dAvC[1]+156.*dAvC[2]-61.*dAvC[3]+10.*dAvC[4]);
 297  cumulant[2] = (3./(2.*pow(fR0,6.)))*(71.*dAvC[0]-118.*dAvC[1]+98.*dAvC[2]-41.*dAvC[3]+7.*dAvC[4]);
 298  cumulant[3] = (-24./pow(fR0,8.))*(14.*dAvC[0]-26.*dAvC[1]+24.*dAvC[2]-11.*dAvC[3]+2.*dAvC[4]);
 299  cumulant[4] = (120./pow(fR0,10.))*(5.*dAvC[0]-10.*dAvC[1]+10.*dAvC[2]-5.*dAvC[3]+1.*dAvC[4]);
 300
 301  /*
 302  cout<<endl;
 303  cout<<"*********************************"<<endl;
 304  cout<<"cumulants:"<<endl;
 305  cout<<" c_"<<cFlow<<"{2} = "<<cumulant[0]<<endl;
 306  cout<<" c_"<<cFlow<<"{4} = "<<cumulant[1]<<endl;
 307  cout<<" c_"<<cFlow<<"{6} = "<<cumulant[2]<<endl;
 308  cout<<" c_"<<cFlow<<"{8} = "<<cumulant[3]<<endl;
 309  cout<<"c_"<<cFlow<<"{10} = "<<cumulant[4]<<endl;
 310  cout<<endl;
 311  */
 312
 313  Double_t dV2=0.,dV4=0.,dV6=0.,dV8=0.,dV10=0.;//integrated flow estimates
 314
 315  if(cumulant[0]>=0.)
 316  {
 317   dV2=pow(cumulant[0],(1./2.));
 318  }
 319  if(cumulant[1]<=0.)
 320  {
 321   dV4=pow(-cumulant[1],(1./4.));
 322  }
 323  if(cumulant[2]>=0.)
 324  {
 325   dV6=pow((1./4.)*cumulant[2],(1./6.));
 326  }
 327  if(cumulant[3]<=0.)
 328  {
 329   dV8=pow(-(1./33.)*cumulant[3],(1./8.));
 330  }
 331  if(cumulant[4]>=0.)
 332  {
 333   dV10=pow((1./456.)*cumulant[4],(1./10.));
 334  }
 335
 336  cout<<endl;
 337  cout<<"***************************************"<<endl;
 338  cout<<"***************************************"<<endl;
 339  cout<<"flow estimates from GF-cumulants:      "<<endl;
 340  cout<<endl;
 341
 342  Double_t sdQ[4]={0.};
 343  Double_t chiQ[4]={0.};
 344
 345  //v_2{2}
 346  if(nEvents>0 && dAvM>0 && dAvw2>0 && cumulant[0]>=0.)
 347  {
 348   if((dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(cumulant[0],(1./2.))*(dAvM/dAvw2),2.)>0.))
 349   {
 350    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
 351   }
 352   if(chiQ[0])
 353   {
 354    sdQ[0]=pow(((1./(2.*dAvM*nEvents))*((1.+2.*pow(chiQ[0],2))/(2.*pow(chiQ[0],2)))),0.5);
 355   }
 356   cout<<"   v_"<<cFlow<<"{2} = "<<dV2<<" +/- "<<sdQ[0]<<endl;
 357   //cout<<" v_"<<cFlow<<"{2} = "<<dV2<<" +/- "<<sdQ[0]<<", chi{2} = "<<chiQ[0]<<endl;//printing also the chi
 358   fifr->SetBinContent(1,dV2);
 359   fifr->SetBinError(1,sdQ[0]);
 360   //filling common histograms:
 361   fchr2nd->FillIntegratedFlow(dV2,sdQ[0]);
 362   fchr2nd->FillChi(chiQ[0]);
 363
 364
 365   //abTempDeleteMe
 366   fchr2nd->FillIntegratedFlowRP(dV2,sdQ[0]);
 367   fchr2nd->FillChiRP(chiQ[0]);
 368   //abTempDeleteMe
 369
 370  }
 371  else
 372  {
 373   cout<<"   v_"<<cFlow<<"{2} = Im"<<endl;
 374  }
 375
 376  //v_2{4}
 377  if(nEvents>0 && dAvM>0 && dAvw2>0 && cumulant[1]<=0.)
 378  {
 379   if((dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(-cumulant[1],(1./4.))*(dAvM/dAvw2),2.)>0.))
 380   {
 381    chiQ[1]=(dAvM*dV4)/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4*dAvM/dAvw2,2.),0.5);
 382   }
 383   if(chiQ[1])
 384   {
 385    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);
 386   }
 387   cout<<"   v_"<<cFlow<<"{4} = "<<dV4<<" +/- "<<sdQ[1]<<endl;
 388   //cout<<" v_"<<cFlow<<"{4} = "<<dV4<<" +/- "<<sdQ[1]<<", chi{4} = "<<chiQ[1]<<endl;//printing also the chi
 389   fifr->SetBinContent(2,dV4);
 390   fifr->SetBinError(2,sdQ[1]);
 391   //filling common histograms:
 392   fchr4th->FillIntegratedFlow(dV4,sdQ[1]);
 393   fchr4th->FillChi(chiQ[1]);
 394
 395
 396   //abTempDeleteMe
 397   fchr4th->FillIntegratedFlowRP(dV4,sdQ[1]);
 398   fchr4th->FillChiRP(chiQ[1]);
 399   //abTempDeleteMe
 400
 401  }
 402  else
 403  {
 404   cout<<"   v_"<<cFlow<<"{4} = Im"<<endl;
 405  }
 406
 407  //v_2{6}
 408  if(nEvents>0 && dAvM>0 && dAvw2>0 && cumulant[2]>=0.)
 409  {
 410   if((dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow((1./4.)*cumulant[2],(1./6.))*(dAvM/dAvw2),2.)>0.))
 411   {
 412    chiQ[2]=(dAvM*dV6)/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV6*dAvM/dAvw2,2.),0.5);
 413   }
 414   if(chiQ[2])
 415   {
 416    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);
 417   }
 418   cout<<"   v_"<<cFlow<<"{6} = "<<dV6<<" +/- "<<sdQ[2]<<endl;
 419   //cout<<" v_"<<cFlow<<"{6} = "<<dV6<<" +/- "<<sdQ[2]<<", chi{6} = "<<chiQ[2]<<endl;//printing also the chi
 420   fifr->SetBinContent(3,dV6);
 421   fifr->SetBinError(3,sdQ[2]);
 422   //filling common histograms:
 423   fchr6th->FillIntegratedFlow(dV6,sdQ[2]);
 424   fchr6th->FillChi(chiQ[2]);
 425
 426
 427   //abTempDeleteMe
 428   fchr6th->FillIntegratedFlowRP(dV6,sdQ[2]);
 429   fchr6th->FillChiRP(chiQ[2]);
 430   //abTempDeleteMe
 431
 432
 433  }
 434  else
 435  {
 436   cout<<"   v_"<<cFlow<<"{6} = Im"<<endl;
 437  }
 438
 439  //v_2{8}
 440  if(nEvents>0 && dAvM>0 && dAvw2>0 && cumulant[3]<=0.)
 441  {
 442   if((dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(-(1./33.)*cumulant[3],(1./8.))*(dAvM/dAvw2),2.)>0.))
 443   {
 444    chiQ[3]=(dAvM*dV8)/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV8*dAvM/dAvw2,2.),0.5);
 445   }
 446   if(chiQ[3])
 447   {
 448    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);
 449   }
 450   cout<<"   v_"<<cFlow<<"{8} = "<<dV8<<" +/- "<<sdQ[3]<<endl;
 451   //cout<<" v_"<<cFlow<<"{8} = "<<dV8<<" +/- "<<sdQ[3]<<", chi{8} = "<<chiQ[3]<<endl;//printing also the chi
 452   fifr->SetBinContent(4,dV8);
 453   fifr->SetBinError(4,sdQ[3]);
 454   //filling common histograms:
 455   fchr8th->FillIntegratedFlow(dV8,sdQ[3]);
 456   fchr8th->FillChi(chiQ[3]);
 457
 458   //abTempDeleteMe
 459   fchr8th->FillIntegratedFlowRP(dV8,sdQ[3]);
 460   fchr8th->FillChiRP(chiQ[3]);
 461   //abTempDeleteMe
 462
 463  }
 464  else
 465  {
 466   cout<<"   v_"<<cFlow<<"{8} = Im"<<endl;
 467  }
 468
 469  /*
 470  //v_2{10}
 471  if (dAvM && cumulant[4]>=0.)
 472  {
 473   cout<<"v_"<<cFlow<<"{10} = "<<dV10<<endl;
 474  }
 475  else
 476  {
 477   cout<<"v_"<<cFlow<<"{10} = Im"<<endl;
 478  }
 479  */
 480
 481  cout<<endl;
 482  cout<<"      nEvts = "<<nEvents<<", AvM = "<<dAvM<<endl;
 483  cout<<"***************************************"<<endl;
 484  cout<<"***************************************"<<endl;
 485
 486
 487  //===========================================================================
 488  //                      RP = Reaction Plane
 489  //===========================================================================
 490
 491  /////////////////////////////////////////////////////////////////////////////
 492  ///////////////////////DIFFERENTIAL FLOW CALCULATIONS////////////////////////
 493  /////////////////////////////////////////////////////////////////////////////
 494
 495   TVectorD ptXRP(cnBinsPt * cPmax * cQmax);
 496   TVectorD ptYRP(cnBinsPt * cPmax * cQmax);
 497   TVectorD ptBinRPNoOfParticles(cnBinsPt);
 498
 499   //3D profiles (for pt)
 500   for(Int_t b=0;b<cnBinsPt;b++)
 501   {
 502    ptBinRPNoOfParticles[b]=fPtBinRPNoOfParticles->GetBinEntries(b+1);
 503
 504    for(Int_t p=0;p<cPmax;p++)
 505    {
 506     for(Int_t q=0;q<cQmax;q++)
 507     {
 508      if(dAvG(p,q))
 509      {
 510       if(fDiffPtRPGenFunRe)
 511       {
 512        ptXRP[index3d(b,p,q,cnBinsPt,cPmax)]=fDiffPtRPGenFunRe->GetBinContent(b+1,p+1,q+1)/dAvG(p,q);
 513       }
 514       if(fDiffPtRPGenFunIm)
 515       {
 516        ptYRP[index3d(b,p,q,cnBinsPt,cPmax)]=fDiffPtRPGenFunIm->GetBinContent(b+1,p+1,q+1)/dAvG(p,q);
 517       }
 518      }
 519     }
 520    }
 521   }
 522
 523   TVectorD etaXRP(cnBinsEta*cPmax*cQmax);
 524   TVectorD etaYRP(cnBinsEta*cPmax*cQmax);
 525   TVectorD etaBinRPNoOfParticles(cnBinsEta);
 526
 527   //3D profiles (for eta)
 528   for(Int_t b=0;b<cnBinsEta;b++)
 529   {
 530    etaBinRPNoOfParticles[b]=fEtaBinRPNoOfParticles->GetBinEntries(b);
 531    for(Int_t p=0;p<cPmax;p++)
 532    {
 533     for(Int_t q=0;q<cQmax;q++)
 534     {
 535      if(dAvG(p,q))
 536      {
 537       if(fDiffEtaRPGenFunRe)
 538       {
 539        etaXRP[index3d(b,p,q,cnBinsEta,cPmax)]=fDiffEtaRPGenFunRe->GetBinContent(b+1,p+1,q+1)/dAvG(p,q);
 540       }
 541       if(fDiffEtaRPGenFunIm)
 542       {
 543        etaYRP[index3d(b,p,q,cnBinsEta,cPmax)]=fDiffEtaRPGenFunIm->GetBinContent(b+1,p+1,q+1)/dAvG(p,q);
 544       }
 545      }
 546     }
 547    }
 548   }
 549
 550
 551   //-------------------------------------------------------------------------------------------------------------------------------
 552   //final results for differential flow (in pt):
 553   TMatrixD ptDRP(cnBinsPt, cPmax);
 554   Double_t tempSumForptDRP=0.;
 555   for (Int_t b=0;b<cnBinsPt;b++)
 556   {
 557    for (Int_t p=0;p<cPmax;p++)
 558    {
 559     tempSumForptDRP=0.;
 560     for (Int_t q=0;q<cQmax;q++)
 561     {
 562      tempSumForptDRP+=cos(cMltpl*2.*q*TMath::Pi()/cQmax)*ptXRP[index3d(b,p,q,cnBinsPt,cPmax)] + sin(cMltpl*2.*q*TMath::Pi()/cQmax)*ptYRP[index3d(b,p,q,cnBinsPt,cPmax)];
 563     }
 564      ptDRP(b,p)=1.*(pow(fR0*pow(p+1.0,0.5),cMltpl)/cQmax)*tempSumForptDRP;
 565    }
 566   }
 567
 568   Double_t ptRPDiffCumulant2[cnBinsPt]={0.};
 569   Double_t ptRPDiffCumulant4[cnBinsPt]={0.};
 570   Double_t ptRPDiffCumulant6[cnBinsPt]={0.};
 571   Double_t ptRPDiffCumulant8[cnBinsPt]={0.};
 572   Double_t ptRPDiffCumulant10[cnBinsPt]={0.};
 573
 574   for (Int_t b=0;b<cnBinsPt;b++)
 575   {
 576    ptRPDiffCumulant2[b]=(1./(fR0*fR0))*(5.*ptDRP(b,0)-5.*ptDRP(b,1)+(10./3.)*ptDRP(b,2)-(5./4.)*ptDRP(b,3)+(1./5.)*ptDRP(b,4));
 577    ptRPDiffCumulant4[b]=(1./pow(fR0,4.))*((-77./6.)*ptDRP(b,0)+(107./6.)*ptDRP(b,1)-(13./1.)*ptDRP(b,2)+(61./12.)*ptDRP(b,3)-(5./6.)*ptDRP(b,4));
 578    ptRPDiffCumulant6[b]=(1./pow(fR0,6.))*((71./2.)*ptDRP(b,0)-59.*ptDRP(b,1)+49.*ptDRP(b,2)-(41./2.)*ptDRP(b,3)+(7./2.)*ptDRP(b,4));
 579    ptRPDiffCumulant8[b]=(1./pow(fR0,8.))*(-84.*ptDRP(b,0)+156.*ptDRP(b,1)-144.*ptDRP(b,2)+66.*ptDRP(b,3)-12.*ptDRP(b,4));
 580    ptRPDiffCumulant10[b]=(1./pow(fR0,10.))*(120.*ptDRP(b,0)-240.*ptDRP(b,1)+240.*ptDRP(b,2)-120.*ptDRP(b,3)+24.*ptDRP(b,4));
 581   }
 582
 583   //diff. flow values per pt bin:
 584   Double_t v2ptRP[cnBinsPt]={0.};
 585   Double_t v4ptRP[cnBinsPt]={0.};
 586   Double_t v6ptRP[cnBinsPt]={0.};
 587   Double_t v8ptRP[cnBinsPt]={0.};
 588
 589   //errrors:
 590   Double_t sdRPDiff2pt[cnBinsPt]={0.};
 591   Double_t sdRPDiff4pt[cnBinsPt]={0.};
 592   //Double_t sdDiff6pt[cnBinsPt]={0.};//to be improved (calculation needed)
 593   //Double_t sdDiff8pt[cnBinsPt]={0.};//to be improved (calculation needed)
 594
 595   //cout<<"number of pt bins: "<<cnBinsPt<<endl;
 596   //cout<<"****************************************"<<endl;
 597   for (Int_t b=0;b<cnBinsPt;b++){
 598     //cout<<"pt bin: "<<b*fBinWidthPt<<"-"<<(b+1)*fBinWidthPt<<" GeV"<<endl;
 599
 600     //v'_{2/2}{2}
 601     if(cumulant[0]>0)
 602     {
 603       v2ptRP[b]=ptRPDiffCumulant2[b]/pow(cumulant[0],0.5);
 604       if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2*dAvM,2.)>0. && ptBinRPNoOfParticles[b]>0.)
 605       {
 606        if(chiQ[0]>0)
 607        {
 608         sdRPDiff2pt[b]=pow((1./(2.*ptBinRPNoOfParticles[b]))*((1.+pow(chiQ[0],2.))/pow(chiQ[0],2.)),0.5);
 609        }
 610        //cout<<"v'_2/2{2} = "<<v2ptRP[b]<<"%, "<<" "<<"sd{2} = "<<100.*sdRPDiff2pt[b]<<"%"<<endl;
 611        fdfr2->SetBinContent(b+1,v2ptRP[b]);
 612        fdfr2->SetBinError(b+1,sdRPDiff2pt[b]);
 613        //common histogram:
 614        fchr2nd->FillDifferentialFlow(b+1,v2ptRP[b],sdRPDiff2pt[b]);
 615
 616
 617        //abTempDeleteMe
 618        fchr2nd->FillDifferentialFlowPtRP(b+1,v2ptRP[b],sdRPDiff2pt[b]);
 619        //abTempDeleteMe
 620
 621       } else {
 622          //cout<<"v'_2/2{2} = Im"<<endl;
 623       }
 624     }else{
 625       //cout<<"v'_2/2{2} = Im"<<endl;
 626     }
 627
 628     //v'_{2/2}{4}
 629     if(cumulant[1]<0)
 630     {
 631       v4ptRP[b]=-ptRPDiffCumulant4[b]/pow(-cumulant[1],.75);
 632       if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4*dAvM,2.)>0.&&ptBinRPNoOfParticles[b]>0.) // to be improved
 633       {
 634        if(chiQ[1]>0)
 635        {
 636         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);
 637        }
 638        //cout<<"v'_2/2{4} = "<<v4ptRP[b]<<"%, "<<" "<<"sd{4} = "<<100.*sdRPDiff4pt[b]<<"%"<<endl;
 639        fdfr4->SetBinContent(b+1,v4ptRP[b]);
 640        fdfr4->SetBinError(b+1,sdRPDiff4pt[b]);
 641        //common histogram:
 642        fchr4th->FillDifferentialFlow(b+1,v4ptRP[b],sdRPDiff4pt[b]);
 643
 644
 645
 646        //abTempDeleteMe
 647        fchr4th->FillDifferentialFlowPtRP(b+1,v4ptRP[b],sdRPDiff4pt[b]);
 648        //abTempDeleteMe
 649
 650
 651
 652       } else {
 653          //cout<<"v'_2/2{4} = Im"<<endl;
 654       }
 655     }else{
 656       //cout<<"v'_2/2{4} = Im"<<endl;
 657     }
 658
 659     //v'_{2/2}{6}
 660     if(cumulant[2]>0){
 661       //cout<<"v'_2/2{6} = "<<100.*ptRPDiffCumulant6[b]/(4.*pow((1./4.)*cumulant[2],(5./6.)))<<"%"<<endl;
 662       v6ptRP[b]=ptRPDiffCumulant6[b]/(4.*pow((1./4.)*cumulant[2],(5./6.)));
 663       fdfr6->SetBinContent(b+1,v6ptRP[b]);
 664       //common histogram:
 665       fchr6th->FillDifferentialFlow(b+1,v6ptRP[b],0.);
 666
 667
 668
 669       //abTempDeleteMe
 670       fchr6th->FillDifferentialFlowPtRP(b+1,v6ptRP[b],0.);
 671       //abTempDeleteMe
 672
 673
 674
 675
 676     }else{
 677       //cout<<"v'_2/2{6} = Im"<<endl;
 678     }
 679
 680     //v'_{2/2}{8}
 681     if(cumulant[3]<0){
 682       //cout<<"v'_2/2{8} = "<<-100.*ptRPDiffCumulant8[b]/(33.*pow(-(1./33.)*cumulant[3],(7./8.)))<<"%"<<endl;
 683       v8ptRP[b]=-ptRPDiffCumulant8[b]/(33.*pow(-(1./33.)*cumulant[3],(7./8.)));
 684       fdfr8->SetBinContent(b+1,v8ptRP[b]);
 685       //common histogram:
 686       fchr8th->FillDifferentialFlow(b+1,v8ptRP[b],0.);
 687
 688
 689
 690       //abTempDeleteMe
 691       fchr8th->FillDifferentialFlowPtRP(b+1,v8ptRP[b],0.);
 692       //abTempDeleteMe
 693
 694
 695
 696     }else{
 697       //cout<<"v'_2/2{8} = Im"<<endl;
 698     }
 699     //cout<<"****************************************"<<endl;
 700   }
 701  //-------------------------------------------------------
 702
 703
 704
 705   ///////////////////////////////////////////////////////////////////////////////////
 706  //////////////////////// INTEGRATED FLOW CALCULATIONS (RP) /////////////////////////
 707  ///////////////////////////////////////////////////////////////////////////////////
 708
 709  Double_t dV2RP=0., dV4RP=0., dV6RP=0., dV8RP=0.;
 710  Double_t dV2RPError=0., dV4RPError=0., dV6RPError=0., dV8RPError=0.;
 711  Double_t dSumOfYieldRP=0.;
 712  for (Int_t b=1;b<cnBinsPt+1;b++)
 713  {
 714   if(fch->GetHistPtRP())
 715   {
 716    dSumOfYieldRP+=(fch->GetHistPtRP())->GetBinContent(b);
 717    if(fchr2nd->GetHistDiffFlowPtRP())
 718    {
 719     dV2RP+=((fchr2nd->GetHistDiffFlowPtRP())->GetBinContent(b))*(fch->GetHistPtRP())->GetBinContent(b);
 720     dV2RPError+=pow((fch->GetHistPtRP())->GetBinContent(b),2.)*pow((fchr2nd->GetHistDiffFlowPtRP())->GetBinError(b),2.);
 721    }
 722    if(fchr4th->GetHistDiffFlowPtRP())
 723    {
 724     dV4RP+=((fchr4th->GetHistDiffFlowPtRP())->GetBinContent(b))*(fch->GetHistPtRP())->GetBinContent(b);
 725     dV4RPError+=pow((fch->GetHistPtRP())->GetBinContent(b),2.)*pow((fchr4th->GetHistDiffFlowPtRP())->GetBinError(b),2.);
 726    }
 727    if(fchr6th->GetHistDiffFlowPtRP())
 728    {
 729     dV6RP+=((fchr6th->GetHistDiffFlowPtRP())->GetBinContent(b))*(fch->GetHistPtRP())->GetBinContent(b);
 730     dV6RPError+=pow((fch->GetHistPtRP())->GetBinContent(b),2.)*pow((fchr6th->GetHistDiffFlowPtRP())->GetBinError(b),2.);
 731    }
 732    if(fchr8th->GetHistDiffFlowPtRP())
 733    {
 734     dV8RP+=((fchr8th->GetHistDiffFlowPtRP())->GetBinContent(b))*(fch->GetHistPtRP())->GetBinContent(b);
 735     dV8RPError+=pow((fch->GetHistPtRP())->GetBinContent(b),2.)*pow((fchr8th->GetHistDiffFlowPtRP())->GetBinError(b),2.);
 736    }
 737   }
 738  }
 739
 740  if(dSumOfYieldRP)
 741  {
 742   dV2RP/=dSumOfYieldRP;
 743   dV2RPError/=(dSumOfYieldRP*dSumOfYieldRP);
 744   dV4RP/=dSumOfYieldRP;
 745   dV4RPError/=(dSumOfYieldRP*dSumOfYieldRP);
 746   dV6RP/=dSumOfYieldRP;
 747   dV6RPError/=(dSumOfYieldRP*dSumOfYieldRP);
 748   dV8RP/=dSumOfYieldRP;
 749   dV8RPError/=(dSumOfYieldRP*dSumOfYieldRP);
 750   fchr2nd->FillIntegratedFlowRP(dV2RP,pow(dV2RPError,0.5));
 751   fchr4th->FillIntegratedFlowRP(dV4RP,pow(dV4RPError,0.5));
 752   fchr6th->FillIntegratedFlowRP(dV6RP,pow(dV6RPError,0.5));//to be improved (errors needed)
 753   fchr8th->FillIntegratedFlowRP(dV8RP,pow(dV8RPError,0.5));//to be improved (errors needed)
 754  }
 755
 756  cout<<endl;
 757  cout<<"***************************************"<<endl;
 758  cout<<"***************************************"<<endl;
 759  cout<<"flow estimates from GF-cumulants (RP):"<<endl;
 760  cout<<endl;
 761
 762  cout<<"   v_"<<cFlow<<"{2} = "<<dV2RP<<" +/- "<<pow(dV2RPError,0.5)<<endl;
 763  cout<<"   v_"<<cFlow<<"{4} = "<<dV4RP<<" +/- "<<pow(dV4RPError,0.5)<<endl;
 764  cout<<"   v_"<<cFlow<<"{6} = "<<dV6RP<<" +/- "<<pow(dV6RPError,0.5)<<endl;
 765  cout<<"   v_"<<cFlow<<"{8} = "<<dV8RP<<" +/- "<<pow(dV8RPError,0.5)<<endl;
 766
 767  cout<<endl;
 768  cout<<"      nEvts = "<<(fch->GetHistMultRP())->GetEntries()<<", AvM = "<<(fch->GetHistMultRP())->GetMean()<<endl;
 769  cout<<"***************************************"<<endl;
 770  cout<<"***************************************"<<endl;
 771
 772
 773
 774
 775
 776
 777
 778
 779
 780
 781
 782
 783
 784
 785
 786
 787
 788
 789
 790
 791
 792
 793
 794
 795
 796
 797
 798
 799
 800
 801  //-------------------------------------------------------------------------------------------------------------------------------
 802   //final results for differential flow (in eta):
 803   TMatrixD etaDRP(cnBinsEta,cPmax);
 804   Double_t tempSumForEtaDRP=0.;
 805   for (Int_t b=0;b<cnBinsEta;b++)
 806   {
 807    for (Int_t p=0;p<cPmax;p++)
 808    {
 809     tempSumForEtaDRP=0.;
 810     for (Int_t q=0;q<cQmax;q++)
 811     {
 812      tempSumForEtaDRP+=cos(cMltpl*2.*q*TMath::Pi()/cQmax)*etaXRP[index3d(b,p,q,cnBinsEta,cPmax)] + sin(cMltpl*2.*q*TMath::Pi()/cQmax)*etaYRP[index3d(b,p,q,cnBinsEta,cPmax)];
 813     }
 814      etaDRP(b,p)=1.*(pow(fR0*pow(p+1.,.5),cMltpl)/cQmax)*tempSumForEtaDRP;
 815    }
 816   }
 817
 818   Double_t etaRPDiffCumulant2[cnBinsEta]={0.};
 819   Double_t etaRPDiffCumulant4[cnBinsEta]={0.};
 820   Double_t etaRPDiffCumulant6[cnBinsEta]={0.};
 821   Double_t etaRPDiffCumulant8[cnBinsEta]={0.};
 822   Double_t etaRPDiffCumulant10[cnBinsEta]={0.};
 823
 824   for (Int_t b=0;b<cnBinsEta;b++)
 825   {
 826    etaRPDiffCumulant2[b]=(1./(fR0*fR0))*(5.*etaDRP(b,0)-5.*etaDRP(b,1)+(10./3.)*etaDRP(b,2)-(5./4.)*etaDRP(b,3)+(1./5.)*etaDRP(b,4));
 827    etaRPDiffCumulant4[b]=(1./pow(fR0,4.))*((-77./6.)*etaDRP(b,0)+(107./6.)*etaDRP(b,1)-(13./1.)*etaDRP(b,2)+(61./12.)*etaDRP(b,3)-(5./6.)*etaDRP(b,4));
 828    etaRPDiffCumulant6[b]=(1./pow(fR0,6.))*((71./2.)*etaDRP(b,0)-59.*etaDRP(b,1)+49.*etaDRP(b,2)-(41./2.)*etaDRP(b,3)+(7./2.)*etaDRP(b,4));
 829    etaRPDiffCumulant8[b]=(1./pow(fR0,8.))*(-84.*etaDRP(b,0)+156.*etaDRP(b,1)-144.*etaDRP(b,2)+66.*etaDRP(b,3)-12.*etaDRP(b,4));
 830    etaRPDiffCumulant10[b]=(1./pow(fR0,10.))*(120.*etaDRP(b,0)-240.*etaDRP(b,1)+240.*etaDRP(b,2)-120.*etaDRP(b,3)+24.*etaDRP(b,4));
 831   }
 832
 833   //diff. flow values per eta bin:
 834   Double_t v2etaRP[cnBinsEta]={0.};
 835   Double_t v4etaRP[cnBinsEta]={0.};
 836   Double_t v6etaRP[cnBinsEta]={0.};
 837   Double_t v8etaRP[cnBinsEta]={0.};
 838
 839   //errrors:
 840   Double_t sdRPDiff2eta[cnBinsEta]={0.};
 841   Double_t sdRPDiff4eta[cnBinsEta]={0.};
 842   //Double_t sdDiff6eta[cnBinsEta]={0.};//to be improved (calculation needed)
 843   //Double_t sdDiff8eta[cnBinsEta]={0.};//to be improved (calculation needed)
 844
 845   //cout<<"number of eta bins: "<<cnBinsEta<<endl;
 846   //cout<<"****************************************"<<endl;
 847   for (Int_t b=0;b<cnBinsEta;b++){
 848     //cout<<"eta bin: "<<b*fBinWidthPt<<"-"<<(b+1)*fBinWidthPt<<" GeV"<<endl;
 849
 850     //v'_{2/2}{2}
 851     if(cumulant[0]>0)
 852     {
 853       v2etaRP[b]=etaRPDiffCumulant2[b]/pow(cumulant[0],0.5);
 854       if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2*dAvM,2.)>0. && etaBinRPNoOfParticles[b]>0.) // to be improved
 855       {
 856        if(chiQ[0]>0)
 857        {
 858         sdRPDiff2eta[b]=pow((1./(2.*etaBinRPNoOfParticles[b]))*((1.+pow(chiQ[0],2.))/pow(chiQ[0],2.)),0.5);
 859        }
 860        //cout<<"v'_2/2{2} = "<<v2etaRP[b]<<"%, "<<" "<<"sd{2} = "<<100.*sdDiff2eta[b]<<"%"<<endl;
 861        //fdfr2->SetBinContent(b+1,v2etaRP[b]);
 862        //fdfr2->SetBinError(b+1,sdDiff2eta[b]);
 863        //common histogram:
 864        //fchr2nd->FillDifferentialFlow(b+1,v2etaRP[b],sdDiff2eta[b]);
 865
 866
 867        //abTempDeleteMe
 868        fchr2nd->FillDifferentialFlowEtaRP(b+1,v2etaRP[b],sdRPDiff2eta[b]);
 869        //abTempDeleteMe
 870
 871
 872       } else {
 873          //cout<<"v'_2/2{2} = Im"<<endl;
 874       }
 875     }else{
 876       //cout<<"v'_2/2{2} = Im"<<endl;
 877     }
 878
 879     //v'_{2/2}{4}
 880     if(cumulant[1]<0)
 881     {
 882       v4etaRP[b]=-etaRPDiffCumulant4[b]/pow(-cumulant[1],0.75);
 883       if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4*dAvM,2.)>0.&&etaBinRPNoOfParticles[b]>0.) // to be improved
 884       {
 885        if(chiQ[1])
 886        {
 887         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);
 888        }
 889        //cout<<"v'_2/2{4} = "<<v4eta[b]<<"%, "<<" "<<"sd{4} = "<<100.*sdDiff4eta[b]<<"%"<<endl;
 890        //fdfr4->SetBinContent(b+1,v4eta[b]);
 891        //fdfr4->SetBinError(b+1,sdDiff4eta[b]);
 892        //common histogram:
 893        //fchr4th->FillDifferentialFlow(b+1,v4eta[b],sdDiff4eta[b]);
 894
 895
 896
 897        //abTempDeleteMe
 898        fchr4th->FillDifferentialFlowEtaRP(b+1,v4etaRP[b],sdRPDiff4eta[b]);
 899        //abTempDeleteMe
 900
 901
 902
 903       } else {
 904          //cout<<"v'_2/2{4} = Im"<<endl;
 905       }
 906     }else{
 907       //cout<<"v'_2/2{4} = Im"<<endl;
 908     }
 909
 910     //v'_{2/2}{6}
 911     if(cumulant[2]>0){
 912       //cout<<"v'_2/2{6} = "<<100.*etaRPDiffCumulant6[b]/(4.*pow((1./4.)*cumulant[2],(5./6.)))<<"%"<<endl;
 913       v6etaRP[b]=etaRPDiffCumulant6[b]/(4.*pow((1./4.)*cumulant[2],(5./6.)));
 914       //fdfr6->SetBinContent(b+1,v6eta[b]);
 915       //common histogram:
 916       //fchr6th->FillDifferentialFlow(b+1,v6eta[b],0.);
 917
 918
 919
 920       //abTempDeleteMe
 921       fchr6th->FillDifferentialFlowEtaRP(b+1,v6etaRP[b],0.);
 922       //abTempDeleteMe
 923
 924
 925
 926
 927     }else{
 928       //cout<<"v'_2/2{6} = Im"<<endl;
 929     }
 930
 931     //v'_{2/2}{8}
 932     if(cumulant[3]<0){
 933       //cout<<"v'_2/2{8} = "<<-100.*etaRPDiffCumulant8[b]/(33.*pow(-(1./33.)*cumulant[3],(7./8.)))<<"%"<<endl;
 934       v8etaRP[b]=-etaRPDiffCumulant8[b]/(33.*pow(-(1./33.)*cumulant[3],(7./8.)));
 935       //fdfr8->SetBinContent(b+1,v8eta[b]);
 936       //common histogram:
 937       //fchr8th->FillDifferentialFlow(b+1,v8eta[b],0.);
 938
 939
 940
 941       //abTempDeleteMe
 942       fchr8th->FillDifferentialFlowEtaRP(b+1,v8etaRP[b],0.);
 943       //abTempDeleteMe
 944
 945
 946
 947     }else{
 948       //cout<<"v'_2/2{8} = Im"<<endl;
 949     }
 950     //cout<<"****************************************"<<endl;
 951   }
 952  //-------------------------------------------------------------------------------------------------------------------------------
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997  //===========================================================================
 998  //                      POI = Particle of Interest
 999  //===========================================================================
1000
1001  /////////////////////////////////////////////////////////////////////////////
1002  ///////////////////////DIFFERENTIAL FLOW CALCULATIONS////////////////////////
1003  /////////////////////////////////////////////////////////////////////////////
1004
1005   TVectorD ptX(cnBinsPt*cPmax*cQmax);
1006   TVectorD ptY(cnBinsPt*cPmax*cQmax);
1007   TVectorD ptBinPOINoOfParticles(cnBinsPt);
1008
1009   //3D profiles (for pt)
1010   for(Int_t b=0;b<cnBinsPt;b++)
1011   {
1012    ptBinPOINoOfParticles[b]=fPtBinPOINoOfParticles->GetBinEntries(b+1);
1013    for(Int_t p=0;p<cPmax;p++)
1014    {
1015     for(Int_t q=0;q<cQmax;q++)
1016     {
1017      if(dAvG(p,q))
1018      {
1019       if(fDiffPtPOIGenFunRe)
1020       {
1021        ptX[index3d(b,p,q,cnBinsPt,cPmax)]=fDiffPtPOIGenFunRe->GetBinContent(b+1,p+1,q+1)/dAvG(p,q);
1022       }
1023       if(fDiffPtPOIGenFunIm)
1024       {
1025        ptY[index3d(b,p,q,cnBinsPt,cPmax)]=fDiffPtPOIGenFunIm->GetBinContent(b+1,p+1,q+1)/dAvG(p,q);
1026       }
1027      }
1028     }
1029    }
1030   }
1031
1032   TVectorD etaX(cnBinsEta*cPmax*cQmax);
1033   TVectorD etaY(cnBinsEta*cPmax*cQmax);
1034   TVectorD etaBinPOINoOfParticles(cnBinsEta);
1035
1036   //3D profiles (for eta)
1037   for(Int_t b=0;b<cnBinsEta;b++)
1038   {
1039    etaBinPOINoOfParticles[b]=fEtaBinPOINoOfParticles->GetBinEntries(b+1);
1040    for(Int_t p=0;p<cPmax;p++)
1041    {
1042     for(Int_t q=0;q<cQmax;q++)
1043     {
1044      if(dAvG(p,q))
1045      {
1046       if(fDiffEtaPOIGenFunRe)
1047       {
1048        etaX[index3d(b,p,q,cnBinsEta,cPmax)]=fDiffEtaPOIGenFunRe->GetBinContent(b+1,p+1,q+1)/dAvG(p,q);
1049       }
1050       if(fDiffEtaPOIGenFunIm)
1051       {
1052        etaY[index3d(b,p,q,cnBinsEta,cPmax)]=fDiffEtaPOIGenFunIm->GetBinContent(b+1,p+1,q+1)/dAvG(p,q);
1053       }
1054      }
1055     }
1056    }
1057   }
1058
1059   /*
1060   if(dAvM){
1061   for(Int_t b=0;b<cnBinsPt;b++){cout<<"***************************************"<<endl;
1062  cout<<"***************************************"<<endl;
1063     //for(Int_t p=0;p<cPmax;p++){
1064       for(Int_t q=0;q<cQmax;q++){
1065         X[b][0][q]=fdRe0->GetBinContent(b+1,q+1)/AvG[0][q];
1066         Y[b][0][q]=fdIm0->GetBinContent(b+1,q+1)/AvG[0][q];
1067         //--------------------------------------------------
1068         X[b][1][q]=fdRe1->GetBinContent(b+1,q+1)/AvG[1][q];
1069         Y[b][1][q]=fdIm1->GetBinContent(b+1,q+1)/AvG[1][q];
1070         //--------------------------------------------------
1071         X[b][2][q]=fdRe2->GetBinContent(b+1,q+1)/AvG[2][q];
1072         Y[b][2][q]=fdIm2->GetBinContent(b+1,q+1)/AvG[2][q];
1073         //--------------------------------------------------
1074         X[b][3][q]=fdRe3->GetBinContent(b+1,q+1)/AvG[3][q];
1075         Y[b][3][q]=fdIm3->GetBinContent(b+1,q+1)/AvG[3][q];
1076         //--------------------------------------------------
1077         X[b][4][q]=fdRe4->GetBinContent(b+1,q+1)/AvG[4][q];
1078         Y[b][4][q]=fdIm4->GetBinContent(b+1,q+1)/AvG[4][q];
1079         //--------------------------------------------------
1080         X[b][5][q]=fdRe5->GetBinContent(b+1,q+1)/AvG[5][q];
1081         Y[b][5][q]=fdIm5->GetBinContent(b+1,q+1)/AvG[5][q];
1082         //--------------------------------------------------
1083         X[b][6][q]=fdRe6->GetBinContent(b+1,q+1)/AvG[6][q];
1084         Y[b][6][q]=fdIm6->GetBinContent(b+1,q+1)/AvG[6][q];
1085         //--------------------------------------------------
1086         X[b][7][q]=fdRe7->GetBinContent(b+1,q+1)/AvG[7][q];
1087         Y[b][7][q]=fdIm7->GetBinContent(b+1,q+1)/AvG[7][q];
1088       }
1089     //}
1090   }
1091   }
1092   */
1093
1094   //-------------------------------------------------------------------------------------------------------------------------------
1095   //final results for differential flow (in pt):
1096   TMatrixD ptD(cnBinsPt,cPmax);
1097   Double_t tempSumForPtD=0.;
1098   for (Int_t b=0;b<cnBinsPt;b++)
1099   {
1100    for (Int_t p=0;p<cPmax;p++)
1101    {
1102     tempSumForPtD=0.;
1103     for (Int_t q=0;q<cQmax;q++)
1104     {
1105      tempSumForPtD+=cos(cMltpl*2.*q*TMath::Pi()/cQmax)*ptX[index3d(b,p,q,cnBinsPt,cPmax)] + sin(cMltpl*2.*q*TMath::Pi()/cQmax)*ptY[index3d(b,p,q,cnBinsPt,cPmax)];
1106     }
1107      ptD(b,p)=1.*(pow(fR0*pow(p+1.0,0.5),cMltpl)/cQmax)*tempSumForPtD;
1108    }
1109   }
1110
1111   Double_t ptDiffCumulant2[cnBinsPt]={0.};
1112   Double_t ptDiffCumulant4[cnBinsPt]={0.};
1113   Double_t ptDiffCumulant6[cnBinsPt]={0.};
1114   Double_t ptDiffCumulant8[cnBinsPt]={0.};
1115   Double_t ptDiffCumulant10[cnBinsPt]={0.};
1116
1117   for (Int_t b=0;b<cnBinsPt;b++)
1118   {
1119    ptDiffCumulant2[b]=(1./(fR0*fR0))*(5.*ptD(b,0)-5.*ptD(b,1)+(10./3.)*ptD(b,2)-(5./4.)*ptD(b,3)+(1./5.)*ptD(b,4));
1120    ptDiffCumulant4[b]=(1./pow(fR0,4.))*((-77./6.)*ptD(b,0)+(107./6.)*ptD(b,1)-(13./1.)*ptD(b,2)+(61./12.)*ptD(b,3)-(5./6.)*ptD(b,4));
1121    ptDiffCumulant6[b]=(1./pow(fR0,6.))*((71./2.)*ptD(b,0)-59.*ptD(b,1)+49.*ptD(b,2)-(41./2.)*ptD(b,3)+(7./2.)*ptD(b,4));
1122    ptDiffCumulant8[b]=(1./pow(fR0,8.))*(-84.*ptD(b,0)+156.*ptD(b,1)-144.*ptD(b,2)+66.*ptD(b,3)-12.*ptD(b,4));
1123    ptDiffCumulant10[b]=(1./pow(fR0,10.))*(120.*ptD(b,0)-240.*ptD(b,1)+240.*ptD(b,2)-120.*ptD(b,3)+24.*ptD(b,4));
1124   }
1125
1126   //diff. flow values per pt bin:
1127   Double_t v2pt[cnBinsPt]={0.};
1128   Double_t v4pt[cnBinsPt]={0.};
1129   Double_t v6pt[cnBinsPt]={0.};
1130   Double_t v8pt[cnBinsPt]={0.};
1131
1132   //errrors:
1133   Double_t sdDiff2pt[cnBinsPt]={0.};
1134   Double_t sdDiff4pt[cnBinsPt]={0.};
1135   //Double_t sdDiff6pt[cnBinsPt]={0.};//to be improved (calculation needed)
1136   //Double_t sdDiff8pt[cnBinsPt]={0.};//to be improved (calculation needed)
1137
1138   //cout<<"number of pt bins: "<<cnBinsPt<<endl;
1139   //cout<<"****************************************"<<endl;
1140
1141
1142     for (Int_t b=0;b<cnBinsPt;b++){
1143     //cout<<"pt bin: "<<b*fBinWidthPt<<"-"<<(b+1)*fBinWidthPt<<" GeV"<<endl;
1144
1145
1146     //v'_{2/2}{2}
1147     if(cumulant[0]>0)
1148     {
1149       v2pt[b]=ptDiffCumulant2[b]/pow(cumulant[0],0.5);
1150       if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2*dAvM,2.)>0.&&ptBinPOINoOfParticles[b]>0.)
1151       {
1152        if(chiQ[0]>0)
1153        {
1154         sdDiff2pt[b]=pow((1./(2.*ptBinPOINoOfParticles[b]))*((1.+pow(chiQ[0],2.))/pow(chiQ[0],2.)),0.5);
1155        }
1156        //cout<<"v'_2/2{2} = "<<v2pt[b]<<"%, "<<" "<<"sd{2} = "<<100.*sdDiff2pt[b]<<"%"<<endl;
1157        fdfr2->SetBinContent(b+1,v2pt[b]);
1158        fdfr2->SetBinError(b+1,sdDiff2pt[b]);
1159        //common histogram:
1160        fchr2nd->FillDifferentialFlow(b+1,v2pt[b],sdDiff2pt[b]);
1161
1162
1163        //abTempDeleteMe
1164        fchr2nd->FillDifferentialFlowPtPOI(b+1,v2pt[b],sdDiff2pt[b]);
1165        //abTempDeleteMe
1166
1167
1168       } else {
1169         //cout<<"v'_2/2{2} = Im"<<endl;
1170       }
1171     }else{
1172       //cout<<"v'_2/2{2} = Im"<<endl;
1173     }
1174
1175     //v'_{2/2}{4}
1176     if(cumulant[1]<0)
1177     {
1178       v4pt[b]=-ptDiffCumulant4[b]/pow(-cumulant[1],.75);
1179       if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4*dAvM,2.)>0.&&ptBinPOINoOfParticles[b]>0.)
1180       {
1181        if(chiQ[1]>0)
1182        {
1183         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);
1184        }
1185        //cout<<"v'_2/2{4} = "<<v4pt[b]<<"%, "<<" "<<"sd{4} = "<<100.*sdDiff4pt[b]<<"%"<<endl;
1186        fdfr4->SetBinContent(b+1,v4pt[b]);
1187        fdfr4->SetBinError(b+1,sdDiff4pt[b]);
1188        //common histogram:
1189        fchr4th->FillDifferentialFlow(b+1,v4pt[b],sdDiff4pt[b]);
1190
1191
1192
1193        //abTempDeleteMe
1194        fchr4th->FillDifferentialFlowPtPOI(b+1,v4pt[b],sdDiff4pt[b]);
1195        //abTempDeleteMe
1196
1197
1198
1199       } else {
1200          //cout<<"v'_2/2{4} = Im"<<endl;
1201       }
1202     }else{
1203       //cout<<"v'_2/2{4} = Im"<<endl;
1204     }
1205
1206     //v'_{2/2}{6}
1207     if(cumulant[2]>0){
1208       //cout<<"v'_2/2{6} = "<<100.*ptDiffCumulant6[b]/(4.*pow((1./4.)*cumulant[2],(5./6.)))<<"%"<<endl;
1209       v6pt[b]=ptDiffCumulant6[b]/(4.*pow((1./4.)*cumulant[2],(5./6.)));
1210       fdfr6->SetBinContent(b+1,v6pt[b]);
1211       //common histogram:
1212       fchr6th->FillDifferentialFlow(b+1,v6pt[b],0.);
1213
1214
1215
1216       //abTempDeleteMe
1217       fchr6th->FillDifferentialFlowPtPOI(b+1,v6pt[b],0.);
1218       //abTempDeleteMe
1219
1220
1221
1222
1223     }else{
1224       //cout<<"v'_2/2{6} = Im"<<endl;
1225     }
1226
1227     //v'_{2/2}{8}
1228     if(cumulant[3]<0){
1229       //cout<<"v'_2/2{8} = "<<-100.*ptDiffCumulant8[b]/(33.*pow(-(1./33.)*cumulant[3],(7./8.)))<<"%"<<endl;
1230       v8pt[b]=-ptDiffCumulant8[b]/(33.*pow(-(1./33.)*cumulant[3],(7./8.)));
1231       fdfr8->SetBinContent(b+1,v8pt[b]);
1232       //common histogram:
1233       fchr8th->FillDifferentialFlow(b+1,v8pt[b],0.);
1234
1235
1236
1237       //abTempDeleteMe
1238       fchr8th->FillDifferentialFlowPtPOI(b+1,v8pt[b],0.);
1239       //abTempDeleteMe
1240
1241
1242
1243     }else{
1244       //cout<<"v'_2/2{8} = Im"<<endl;
1245     }
1246     //cout<<"****************************************"<<endl;
1247   }
1248  //-------------------------------------------------------------------------------------------------------------------------------
1249
1250
1251
1252
1253
1254  ///////////////////////////////////////////////////////////////////////////////////
1255  //////////////////////// INTEGRATED FLOW CALCULATIONS (POI) ///////////////////////
1256  ///////////////////////////////////////////////////////////////////////////////////
1257
1258  Double_t dV2POI=0., dV4POI=0., dV6POI=0., dV8POI=0.;
1259  Double_t dV2POIError=0., dV4POIError=0., dV6POIError=0., dV8POIError=0.;
1260  Double_t dSumOfYieldPOI=0.;
1261  for (Int_t b=1;b<cnBinsPt+1;b++)
1262  {
1263   if(fch->GetHistPtPOI())
1264   {
1265    dSumOfYieldPOI+=(fch->GetHistPtPOI())->GetBinContent(b);
1266    if(fchr2nd->GetHistDiffFlowPtPOI())
1267    {
1268     dV2POI+=((fchr2nd->GetHistDiffFlowPtPOI())->GetBinContent(b))*(fch->GetHistPtPOI())->GetBinContent(b);
1269     dV2POIError+=pow((fch->GetHistPtPOI())->GetBinContent(b),2.)*pow((fchr2nd->GetHistDiffFlowPtPOI())->GetBinError(b),2.);
1270    }
1271    if(fchr4th->GetHistDiffFlowPtPOI())
1272    {
1273     dV4POI+=((fchr4th->GetHistDiffFlowPtPOI())->GetBinContent(b))*(fch->GetHistPtPOI())->GetBinContent(b);
1274     dV4POIError+=pow((fch->GetHistPtPOI())->GetBinContent(b),2.)*pow((fchr4th->GetHistDiffFlowPtPOI())->GetBinError(b),2.);
1275    }
1276    if(fchr6th->GetHistDiffFlowPtPOI())
1277    {
1278     dV6POI+=((fchr6th->GetHistDiffFlowPtPOI())->GetBinContent(b))*(fch->GetHistPtPOI())->GetBinContent(b);
1279     dV6POIError+=pow((fch->GetHistPtPOI())->GetBinContent(b),2.)*pow((fchr6th->GetHistDiffFlowPtPOI())->GetBinError(b),2.);
1280    }
1281    if(fchr8th->GetHistDiffFlowPtPOI())
1282    {
1283     dV8POI+=((fchr8th->GetHistDiffFlowPtPOI())->GetBinContent(b))*(fch->GetHistPtPOI())->GetBinContent(b);
1284     dV8POIError+=pow((fch->GetHistPtPOI())->GetBinContent(b),2.)*pow((fchr8th->GetHistDiffFlowPtPOI())->GetBinError(b),2.);
1285    }
1286   }
1287  }
1288
1289  if(dSumOfYieldPOI)
1290  {
1291   dV2POI/=dSumOfYieldPOI;
1292   dV2POIError/=(dSumOfYieldPOI*dSumOfYieldPOI);
1293   dV4POI/=dSumOfYieldPOI;
1294   dV4POIError/=(dSumOfYieldPOI*dSumOfYieldPOI);
1295   dV6POI/=dSumOfYieldPOI;
1296   dV6POIError/=(dSumOfYieldPOI*dSumOfYieldPOI);
1297   dV8POI/=dSumOfYieldPOI;
1298   dV8POIError/=(dSumOfYieldPOI*dSumOfYieldPOI);
1299   fchr2nd->FillIntegratedFlowPOI(dV2POI,pow(dV2POIError,0.5));
1300   fchr4th->FillIntegratedFlowPOI(dV4POI,pow(dV4POIError,0.5));
1301   fchr6th->FillIntegratedFlowPOI(dV6POI,pow(dV6POIError,0.5));//to be improved (errors needed)
1302   fchr8th->FillIntegratedFlowPOI(dV8POI,pow(dV8POIError,0.5));//to be improved (errors needed)
1303  }
1304
1305  cout<<endl;
1306  cout<<"***************************************"<<endl;
1307  cout<<"***************************************"<<endl;
1308  cout<<"flow estimates from GF-cumulants (POI):"<<endl;
1309  cout<<endl;
1310
1311  cout<<"   v_"<<cFlow<<"{2} = "<<dV2POI<<" +/- "<<pow(dV2POIError,0.5)<<endl;
1312  cout<<"   v_"<<cFlow<<"{4} = "<<dV4POI<<" +/- "<<pow(dV4POIError,0.5)<<endl;
1313  cout<<"   v_"<<cFlow<<"{6} = "<<dV6POI<<" +/- "<<pow(dV6POIError,0.5)<<endl;
1314  cout<<"   v_"<<cFlow<<"{8} = "<<dV8POI<<" +/- "<<pow(dV8POIError,0.5)<<endl;
1315
1316  cout<<endl;
1317  cout<<"      nEvts = "<<(fch->GetHistMultPOI())->GetEntries()<<", AvM = "<<(fch->GetHistMultPOI())->GetMean()<<endl;
1318  cout<<"***************************************"<<endl;
1319  cout<<"***************************************"<<endl;
1320  cout<<endl;
1321
1322   //-------------------------------------------------------------------------------------------------------------------------------
1323   //final results for differential flow (in eta):
1324   TMatrixD etaD(cnBinsEta, cPmax);
1325   Double_t tempSumForEtaD=0.;
1326   for (Int_t b=0;b<cnBinsEta;b++)
1327   {
1328    for (Int_t p=0;p<cPmax;p++)
1329    {
1330     tempSumForEtaD=0.;
1331     for (Int_t q=0;q<cQmax;q++)
1332     {
1333      tempSumForEtaD+=cos(cMltpl*2.*q*TMath::Pi()/cQmax)*etaX[index3d(b,p,q,cnBinsEta,cPmax)] + sin(cMltpl*2.*q*TMath::Pi()/cQmax)*etaY[index3d(b,p,q,cnBinsEta,cPmax)];
1334     }
1335      etaD(b,p)=1.*(pow(fR0*pow(p+1.,.5),cMltpl)/cQmax)*tempSumForEtaD;
1336    }
1337   }
1338
1339   Double_t etaDiffCumulant2[cnBinsEta]={0.};
1340   Double_t etaDiffCumulant4[cnBinsEta]={0.};
1341   Double_t etaDiffCumulant6[cnBinsEta]={0.};
1342   Double_t etaDiffCumulant8[cnBinsEta]={0.};
1343   Double_t etaDiffCumulant10[cnBinsEta]={0.};
1344
1345   for (Int_t b=0;b<cnBinsEta;b++)
1346   {
1347    etaDiffCumulant2[b]=(1./(fR0*fR0))*(5.*etaD(b,0)-5.*etaD(b,1)+(10./3.)*etaD(b,2)-(5./4.)*etaD(b,3)+(1./5.)*etaD(b,4));
1348    etaDiffCumulant4[b]=(1./pow(fR0,4.))*((-77./6.)*etaD(b,0)+(107./6.)*etaD(b,1)-(13./1.)*etaD(b,2)+(61./12.)*etaD(b,3)-(5./6.)*etaD(b,4));
1349    etaDiffCumulant6[b]=(1./pow(fR0,6.))*((71./2.)*etaD(b,0)-59.*etaD(b,1)+49.*etaD(b,2)-(41./2.)*etaD(b,3)+(7./2.)*etaD(b,4));
1350    etaDiffCumulant8[b]=(1./pow(fR0,8.))*(-84.*etaD(b,0)+156.*etaD(b,1)-144.*etaD(b,2)+66.*etaD(b,3)-12.*etaD(b,4));
1351    etaDiffCumulant10[b]=(1./pow(fR0,10.))*(120.*etaD(b,0)-240.*etaD(b,1)+240.*etaD(b,2)-120.*etaD(b,3)+24.*etaD(b,4));
1352   }
1353
1354   //diff. flow values per eta bin:
1355   Double_t v2eta[cnBinsEta]={0.};
1356   Double_t v4eta[cnBinsEta]={0.};
1357   Double_t v6eta[cnBinsEta]={0.};
1358   Double_t v8eta[cnBinsEta]={0.};
1359
1360   //errrors:
1361   Double_t sdDiff2eta[cnBinsEta]={0.};
1362   Double_t sdDiff4eta[cnBinsEta]={0.};
1363   //Double_t sdDiff6eta[cnBinsEta]={0.};//to be improved (calculation needed)
1364   //Double_t sdDiff8eta[cnBinsEta]={0.};//to be improved (calculation needed)
1365
1366   //cout<<"number of eta bins: "<<cnBinsEta<<endl;
1367   //cout<<"****************************************"<<endl;
1368   for (Int_t b=0;b<cnBinsEta;b++){
1369     //cout<<"eta bin: "<<b*fBinWidthPt<<"-"<<(b+1)*fBinWidthPt<<" GeV"<<endl;
1370
1371     //v'_{2/2}{2}
1372     if(cumulant[0]>0)
1373     {
1374       v2eta[b]=etaDiffCumulant2[b]/pow(cumulant[0],0.5);
1375       if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2*dAvM,2.)>0.&&etaBinPOINoOfParticles[b]>0.) // to be improved
1376       {
1377        if(chiQ[0]>0)
1378        {
1379         sdDiff2eta[b]=pow((1./(2.*etaBinPOINoOfParticles[b]))*((1.+pow(chiQ[0],2.))/pow(chiQ[0],2.)),0.5);
1380        }
1381        //cout<<"v'_2/2{2} = "<<v2eta[b]<<"%, "<<" "<<"sd{2} = "<<100.*sdDiff2eta[b]<<"%"<<endl;
1382        fdfr2->SetBinContent(b+1,v2eta[b]);
1383        fdfr2->SetBinError(b+1,sdDiff2eta[b]);
1384        //common histogram:
1385        //fchr2nd->FillDifferentialFlow(b+1,v2eta[b],sdDiff2eta[b]);
1386
1387
1388        //abTempDeleteMe
1389        fchr2nd->FillDifferentialFlowEtaPOI(b+1,v2eta[b],sdDiff2eta[b]);
1390        //abTempDeleteMe
1391
1392
1393       } else {
1394          //cout<<"v'_2/2{2} = Im"<<endl;
1395       }
1396     }else{
1397       //cout<<"v'_2/2{2} = Im"<<endl;
1398     }
1399
1400     //v'_{2/2}{4}
1401     if(cumulant[1]<0)
1402     {
1403       v4eta[b]=-etaDiffCumulant4[b]/pow(-cumulant[1],0.75);
1404       if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4*dAvM,2.)>0.&&etaBinPOINoOfParticles[b]>0.) // to be improved
1405       {
1406        if(chiQ[1]>0)
1407        {
1408         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);
1409        }
1410        //cout<<"v'_2/2{4} = "<<v4eta[b]<<"%, "<<" "<<"sd{4} = "<<100.*sdDiff4eta[b]<<"%"<<endl;
1411        fdfr4->SetBinContent(b+1,v4eta[b]);
1412        fdfr4->SetBinError(b+1,sdDiff4eta[b]);
1413        //common histogram:
1414        //fchr4th->FillDifferentialFlow(b+1,v4eta[b],sdDiff4eta[b]);
1415
1416
1417
1418        //abTempDeleteMe
1419        fchr4th->FillDifferentialFlowEtaPOI(b+1,v4eta[b],sdDiff4eta[b]);
1420        //abTempDeleteMe
1421
1422
1423
1424       } else {
1425          //cout<<"v'_2/2{4} = Im"<<endl;
1426       }
1427     }else{
1428       //cout<<"v'_2/2{4} = Im"<<endl;
1429     }
1430
1431     //v'_{2/2}{6}
1432     if(cumulant[2]>0){
1433       //cout<<"v'_2/2{6} = "<<100.*etaDiffCumulant6[b]/(4.*pow((1./4.)*cumulant[2],(5./6.)))<<"%"<<endl;
1434       v6eta[b]=etaDiffCumulant6[b]/(4.*pow((1./4.)*cumulant[2],(5./6.)));
1435       //fdfr6->SetBinContent(b+1,v6eta[b]);
1436       //common histogram:
1437       //fchr6th->FillDifferentialFlow(b+1,v6eta[b],0.);
1438
1439
1440
1441       //abTempDeleteMe
1442       fchr6th->FillDifferentialFlowEtaPOI(b+1,v6eta[b],0.);
1443       //abTempDeleteMe
1444
1445
1446
1447
1448     }else{
1449       //cout<<"v'_2/2{6} = Im"<<endl;
1450     }
1451
1452     //v'_{2/2}{8}
1453     if(cumulant[3]<0){
1454       //cout<<"v'_2/2{8} = "<<-100.*etaDiffCumulant8[b]/(33.*pow(-(1./33.)*cumulant[3],(7./8.)))<<"%"<<endl;
1455       v8eta[b]=-etaDiffCumulant8[b]/(33.*pow(-(1./33.)*cumulant[3],(7./8.)));
1456       //fdfr8->SetBinContent(b+1,v8eta[b]);
1457       //common histogram:
1458       //fchr8th->FillDifferentialFlow(b+1,v8eta[b],0.);
1459
1460
1461
1462       //abTempDeleteMe
1463       fchr8th->FillDifferentialFlowEtaPOI(b+1,v8eta[b],0.);
1464       //abTempDeleteMe
1465
1466
1467
1468     }else{
1469       //cout<<"v'_2/2{8} = Im"<<endl;
1470     }
1471     //cout<<"****************************************"<<endl;
1472   }
1473  //-------------------------------------------------------------------------------------------------------------------------------
1474
1475
1476
1477
1478
1479
1480
1481
1482  //off the record: numerical equations for cumulants solved up to different highest order
1483  if(fOtherEquations)
1484  {
1485
1486  //==============================================================================================================================================
1487
1488  //up to 4th order
1489  //avarage selected multiplicity
1490  Double_t dAvM4=0.;
1491  if(fAvMult4)
1492  {
1493   dAvM4=fAvMult4->GetBinContent(1);
1494  }
1495
1496  //number of events
1497  Int_t nEvents4=0;
1498  if(fAvMult4)
1499  {
1500   nEvents4=(Int_t)(fAvMult4->GetBinEntries(1));
1501  }
1502
1503  TMatrixD dAvG4(cPmax4,cQmax4);
1504  Bool_t someAvGEntryIsNegative4=kFALSE;
1505  for(Int_t p=0;p<cPmax4;p++)
1506  {
1507   for(Int_t q=0;q<cQmax4;q++)
1508   {
1509    dAvG4(p,q)=fIntGenFun4->GetBinContent(p+1,q+1);
1510    if(dAvG4(p,q)<0.)
1511    {
1512     someAvGEntryIsNegative4=kTRUE;
1513    }
1514   }
1515  }
1516  /////////////////////////////////////////////////////////////////////////////
1517  //////////////////gen. function for the cumulants////////////////////////////
1518  /////////////////////////////////////////////////////////////////////////////
1519
1520  TMatrixD dC4(cPmax4,cQmax4);//C[p][q]
1521  if(dAvM4>0 && someAvGEntryIsNegative4==kFALSE)
1522  {
1523   for (Int_t p=0;p<cPmax4;p++)
1524   {
1525    for (Int_t q=0;q<cQmax4;q++)
1526    {
1527     dC4(p,q)=1.*dAvM4*(pow(dAvG4(p,q),(1./dAvM4))-1.);
1528    }
1529   }
1530  }
1531  /////////////////////////////////////////////////////////////////////////////
1532  ///////avaraging the gen. function for the cumulants over azimuth////////////
1533  /////////////////////////////////////////////////////////////////////////////
1534
1535  TVectorD dAvC4(cPmax4);//<C[p][q]>
1536  for (Int_t p=0;p<cPmax4;p++)
1537  {
1538   Double_t tempHere4=0.;
1539   for (Int_t q=0;q<cQmax4;q++)
1540   {
1541    tempHere4+=1.*dC4(p,q);
1542   }
1543   dAvC4[p]=1.*tempHere4/cQmax4;
1544  }
1545
1546  /////////////////////////////////////////////////////////////////////////////
1547  //////////////////////////////////final results//////////////////////////////
1548  /////////////////////////////////////////////////////////////////////////////
1549
1550  TVectorD cumulant4(cPmax4);//array to store various order cumulants
1551
1552  cumulant4[0]=(1./(fR0*fR0))*(2.*dAvC[0]-(1./2.)*dAvC[1]);
1553  cumulant4[1]=(2./pow(fR0,4.))*((-2.)*dAvC[0]+1.*dAvC[1]);
1554
1555  /*
1556  cout<<endl;
1557  cout<<"*********************************"<<endl;
1558  cout<<"cumulants:"<<endl;
1559  cout<<" c_"<<cFlow<<"{2} = "<<cumulant4[0]<<endl;
1560  cout<<" c_"<<cFlow<<"{4} = "<<cumulant4[1]<<endl;
1561  cout<<endl;
1562  */
1563
1564  Double_t dV2o4=0.,dV4o4=0.;
1565
1566  if(cumulant4[0]>=0.)
1567  {
1568   dV2o4 = pow(cumulant4[0],(1./2.));
1569  }
1570  if(cumulant4[1]<=0.)
1571  {
1572   dV4o4 = pow(-cumulant4[1],(1./4.));
1573  }
1574
1575  cout<<endl;
1576  cout<<"***********************************"<<endl;
1577  cout<<"***********************************"<<endl;
1578  cout<<"flow estimates from GF-cumulants:"<<endl;
1579  cout<<"  (calculated up to 4th order)   "<<endl;
1580  cout<<endl;
1581
1582  Double_t sdQo4[2]={0.};
1583  Double_t chiQo4[2]={0.};
1584
1585  //v_2{2}
1586  if(nEvents4 && dAvM4 && cumulant4[0]>=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(cumulant4[0],(1./2.))*dAvM4,2.)>0.))
1587  {
1588   chiQo4[0]=dAvM4*dV2o4/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2o4*dAvM4,2.),0.5);
1589   if(chiQo4[0])
1590   {
1591    sdQo4[0]=pow(((1./(2.*dAvM4*nEvents4))*((1.+2.*pow(chiQo4[0],2))/(2.*pow(chiQo4[0],2)))),0.5);
1592   }
1593   cout<<" v_"<<cFlow<<"{2} = "<<dV2o4<<" +/- "<<sdQo4[0]<<endl;
1594   //cout<<" v_"<<cFlow<<"{2} = "<<dV2o4<<" +/- "<<sdQo4[0]<<", chi{2} = "<<chiQo4[0]<<endl;//printing also the chi
1595  }
1596  else
1597  {
1598   cout<<" v_"<<cFlow<<"{2} = Im"<<endl;
1599  }
1600
1601  //v_2{4}
1602  if(nEvents4 && dAvM4 && cumulant4[1]<=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(-cumulant4[1],(1./4.))*dAvM4,2.)>0.))
1603  {
1604   chiQo4[1]=dAvM4*dV4o4/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4o4*dAvM4,2.),0.5);
1605   if(chiQo4[1])
1606   {
1607    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);
1608   }
1609   cout<<" v_"<<cFlow<<"{4} = "<<dV4o4<<" +/- "<<sdQo4[1]<<endl;
1610   //cout<<" v_"<<cFlow<<"{4} = "<<dV4o4<<" +/- "<<sdQo4[1]<<", chi{4} = "<<chiQo4[1]<<endl;//printing also the chi
1611  }
1612  else
1613  {
1614   cout<<" v_"<<cFlow<<"{4} = Im"<<endl;
1615  }
1616
1617  cout<<endl;
1618  cout<<"   nEvts = "<<nEvents4<<", AvM = "<<dAvM4<<endl;
1619  cout<<"***********************************"<<endl;
1620  cout<<"***********************************"<<endl;
1621
1622  //==============================================================================================================================================
1623
1624  //up to 6th order
1625  //avarage selected multiplicity
1626  Double_t dAvM6=0.;
1627  if(fAvMult6)
1628  {
1629   dAvM6=fAvMult6->GetBinContent(1);
1630  }
1631
1632  //number of events
1633  Int_t nEvents6=0;
1634  if(fAvMult6)
1635  {
1636   nEvents6=(Int_t)(fAvMult6->GetBinEntries(1));
1637  }
1638
1639  TMatrixD dAvG6(cPmax6,cQmax6);
1640  Bool_t someAvGEntryIsNegative6=kFALSE;
1641  for(Int_t p=0;p<cPmax6;p++)
1642  {
1643   for(Int_t q=0;q<cQmax6;q++)
1644   {
1645    dAvG6(p,q)=fIntGenFun6->GetBinContent(p+1,q+1);
1646    if(dAvG6(p,q)<0.)
1647    {
1648     someAvGEntryIsNegative6=kTRUE;
1649    }
1650   }
1651  }
1652  /////////////////////////////////////////////////////////////////////////////
1653  //////////////////gen. function for the cumulants////////////////////////////
1654  /////////////////////////////////////////////////////////////////////////////
1655
1656  TMatrixD dC6(cPmax6,cQmax6);//C[p][q]
1657  if(dAvM6>0 && someAvGEntryIsNegative6==kFALSE)
1658  {
1659   for (Int_t p=0;p<cPmax6;p++)
1660   {
1661    for (Int_t q=0;q<cQmax6;q++)
1662    {
1663     dC6(p,q)=1.*dAvM6*(pow(dAvG6(p,q),(1./dAvM6))-1.);
1664    }
1665   }
1666  }
1667
1668  /////////////////////////////////////////////////////////////////////////////
1669  ///////avaraging the gen. function for the cumulants over azimuth////////////
1670  /////////////////////////////////////////////////////////////////////////////
1671
1672  TVectorD dAvC6(cPmax6);//<etBinContent(1)C[p][q]>
1673  Double_t tempHere6=0.;
1674  for (Int_t p=0;p<cPmax6;p++){
1675   tempHere6=0.;
1676   for (Int_t q=0;q<cQmax6;q++){
1677    tempHere6+=1.*dC6(p,q);
1678   }
1679   dAvC6[p]=1.*tempHere6/cQmax6;
1680  }
1681
1682  /////////////////////////////////////////////////////////////////////////////
1683  //////////////////////////////////final results//////////////////////////////
1684  /////////////////////////////////////////////////////////////////////////////
1685
1686  TVectorD cumulant6(cPmax6);//array to store various order cumulants
1687  cumulant6[0] = (1./(fR0*fR0))*(3.*dAvC[0]-(3./2.)*dAvC[1]+(1./3.)*dAvC[2]);
1688  cumulant6[1] = (2./pow(fR0,4.))*((-5.)*dAvC[0]+4.*dAvC[1]-1.*dAvC[2]);
1689  cumulant6[2] = (6./pow(fR0,6.))*(3.*dAvC[0]-3.*dAvC[1]+1.*dAvC[2]);
1690
1691  /*
1692  cout<<endl;
1693  cout<<"*********************************"<<endl;
1694  cout<<"cumulants:"<<endl;
1695  cout<<" c_"<<cFlow<<"{2} = "<<cumulant6[0]<<endl;
1696  cout<<" c_"<<cFlow<<"{4} = "<<cumulant6[1]<<endl;
1697  cout<<" c_"<<cFlow<<"{6} = "<<cumulant6[2]<<endl;
1698  cout<<endl;
1699  */
1700
1701  Double_t dV2o6=0.,dV4o6=0.,dV6o6=0.;
1702
1703  if(cumulant6[0]>=0.)
1704  {
1705   dV2o6 = pow(cumulant6[0],(1./2.));
1706  }
1707  if(cumulant6[1]<=0.)
1708  {
1709   dV4o6 = pow(-cumulant6[1],(1./4.));
1710  }
1711  if(cumulant6[2]>=0.)
1712  {
1713   dV6o6 = pow((1./4.)*cumulant6[2],(1./6.));
1714  }
1715
1716  cout<<endl;
1717  cout<<"***********************************"<<endl;
1718  cout<<"***********************************"<<endl;
1719  cout<<"flow estimates from GF-cumulants:"<<endl;
1720  cout<<"  (calculated up to 6th order)   "<<endl;
1721  cout<<endl;
1722
1723  Double_t sdQo6[3]={0.};
1724  Double_t chiQo6[3]={0.};
1725
1726  //v_2{2}
1727  if(nEvents6 && dAvM6 && cumulant6[0]>=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(cumulant6[0],(1./2.))*dAvM6,2.)>0.))
1728  {
1729   chiQo6[0]=dAvM6*dV2o6/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2o6*dAvM6,2.),0.5);
1730   if(chiQo6[0])
1731   {
1732    sdQo6[0]=pow(((1./(2.*dAvM6*nEvents6))*((1.+2.*pow(chiQo6[0],2))/(2.*pow(chiQo6[0],2)))),0.5);
1733   }
1734   cout<<" v_"<<cFlow<<"{2} = "<<dV2o6<<" +/- "<<sdQo6[0]<<endl;
1735   //cout<<" v_"<<cFlow<<"{2} = "<<dV2o6<<" +/- "<<sdQo6[0]<<", chi{2} = "<<chiQo6[0]<<endl;//printing also the chi
1736  }
1737  else
1738  {
1739   cout<<" v_"<<cFlow<<"{2} = Im"<<endl;
1740  }
1741
1742  //v_2{4}
1743  if(nEvents6 && dAvM6 && cumulant6[1]<=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(-cumulant6[1],(1./4.))*dAvM6,2.)>0.))
1744  {
1745   chiQo6[1]=dAvM6*dV4o6/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4o6*dAvM6,2.),0.5);
1746   if(chiQo6[1])
1747   {
1748    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);
1749   }
1750   cout<<" v_"<<cFlow<<"{4} = "<<dV4o6<<" +/- "<<sdQo6[1]<<endl;
1751   //cout<<" v_"<<cFlow<<"{4} = "<<dV4o6<<" +/- "<<sdQo6[1]<<", chi{4} = "<<chiQo6[1]<<endl;//printing also the chi
1752  }
1753  else
1754  {
1755   cout<<" v_"<<cFlow<<"{4} = Im"<<endl;
1756  }
1757
1758  //v_2{6}
1759  if(nEvents6 && dAvM6 && cumulant6[2]>=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow((1./4.)*cumulant6[2],(1./6.))*dAvM6,2.)>0.))
1760  {
1761   chiQo6[2]=dAvM6*dV6/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV6o6*dAvM6,2.),0.5);
1762   if(chiQo6[2])
1763   {
1764    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);
1765   }
1766    cout<<" v_"<<cFlow<<"{6} = "<<dV6o6<<" +/- "<<sdQo6[2]<<endl;
1767    //cout<<" v_"<<cFlow<<"{6} = "<<dV6o6<<" +/- "<<sdQo6[2]<<", chi{6} = "<<chiQo6[2]<<endl;//printing also the chi
1768  }
1769  else
1770  {
1771   cout<<" v_"<<cFlow<<"{6} = Im"<<endl;
1772  }
1773
1774  cout<<endl;
1775  cout<<"   nEvts = "<<nEvents6<<", AvM = "<<dAvM6<<endl;
1776  cout<<"***********************************"<<endl;
1777  cout<<"***********************************"<<endl;
1778
1779  //==============================================================================================================================================
1780
1781  //up to 8th order
1782  //avarage selected multiplicity
1783  Double_t dAvM8=0.;
1784  if(fAvMult8)
1785  {
1786   dAvM8=fAvMult8->GetBinContent(1);
1787  }
1788
1789  //number of events
1790  Int_t nEvents8=0;
1791  if(fAvMult8)
1792  {
1793   nEvents8=(Int_t)(fAvMult8->GetBinEntries(1));
1794  }
1795
1796  TMatrixD dAvG8(cPmax8,cQmax8);
1797  Bool_t someAvGEntryIsNegative8=kFALSE;
1798  for(Int_t p=0;p<cPmax8;p++)
1799  {
1800   for(Int_t q=0;q<cQmax8;q++)
1801   {
1802    dAvG8(p,q)=fIntGenFun8->GetBinContent(p+1,q+1);
1803    if(dAvG8(p,q)<0.)
1804    {
1805     someAvGEntryIsNegative8=kTRUE;
1806    }
1807   }
1808  }
1809  /////////////////////////////////////////////////////////////////////////////
1810  //////////////////gen. function for the cumulants////////////////////////////
1811  /////////////////////////////////////////////////////////////////////////////
1812
1813  TMatrixD dC8(cPmax8,cQmax8);//C[p][q]
1814  if(dAvM8>0 && someAvGEntryIsNegative8==kFALSE)
1815  {
1816   for (Int_t p=0;p<cPmax8;p++)
1817   {
1818    for (Int_t q=0;q<cQmax8;q++)
1819    {
1820     dC8(p,q)=1.*dAvM8*(pow(dAvG8(p,q),(1./dAvM8))-1.);
1821    }
1822   }
1823  }
1824
1825  /////////////////////////////////////////////////////////////////////////////
1826  ///////avaraging the gen. function for the cumulants over azimuth////////////
1827  /////////////////////////////////////////////////////////////////////////////
1828
1829  TVectorD dAvC8(cPmax8);//<C[p][q]>
1830  Double_t tempHere8=0.;
1831  for (Int_t p=0;p<cPmax8;p++)
1832  {
1833   tempHere8=0.;
1834   for (Int_t q=0;q<cQmax8;q++)
1835   {
1836    tempHere8+=1.*dC8(p,q);
1837   }
1838   dAvC8[p]=1.*tempHere8/cQmax8;
1839  }
1840
1841  /////////////////////////////////////////////////////////////////////////////
1842  //////////////////////////////////final results//////////////////////////////
1843  /////////////////////////////////////////////////////////////////////////////
1844
1845  TVectorD cumulant8(cPmax8);//array to store various order cumulants
1846  cumulant8[0] = (1./(fR0*fR0))*(4.*dAvC[0]-3.*dAvC[1]+(4./3.)*dAvC[2]-(1./4.)*dAvC[3]);
1847  cumulant8[1] = (1./pow(fR0,4.))*((-52./3.)*dAvC[0]+19.*dAvC[1]-(28./3.)*dAvC[2]+(11./6.)*dAvC[3]);
1848  cumulant8[2] = (3./pow(fR0,6.))*(18.*dAvC[0]-24.*dAvC[1]+14.*dAvC[2]-3.*dAvC[3]);
1849  cumulant8[3] = (24./pow(fR0,8.))*((-4.)*dAvC[0]+6.*dAvC[1]-4.*dAvC[2]+1.*dAvC[3]);
1850
1851  /* x0,y0 ~ 1/sqrt(p) (setting another complex mesh, doesn't work correctly)
1852  cumulant8[0] = (-1./(6.*fR0*fR0)) * (1.*dAvC[0] - 48.*dAvC[1] + 243.*dAvC[2] - 256.*dAvC[3]);
1853  cumulant8[1] = (2./pow(fR0,4.)) * (3.*dAvC[0] - 128.*dAvC[1] + 567.*dAvC[2] - 512.*dAvC[3]);
1854  cumulant8[2] = (-12./pow(fR0,6.)) * (13.*dAvC[0] - 456.*dAvC[1] + 1701.*dAvC[2] - 1408.*dAvC[3]);
1855  cumulant8[3] = (2304./pow(fR0,8.)) * (1.*dAvC[0] - 24.*dAvC[1] + 81.*dAvC[2] - 64.*dAvC[3]);
1856  */
1857
1858  /* x0,y0 ~ p (setting another complex mesh, doesn't work correctly)
1859  cumulant8[0] = (-1./(5040.*fR0*fR0)) * ((-8064.)*dAvC[0] + 1008.*dAvC[1] - 128.*dAvC[2] + 9.*dAvC[3]);
1860  cumulant8[1] = (1./(720.*pow(fR0,4.))) * (1952.*dAvC[0] - 676.*dAvC[1] + 96.*dAvC[2] - 7.*dAvC[3]);
1861  cumulant8[2] = (-1./(40.*pow(fR0,6.))) * ((-116.)*dAvC[0] + 52.*dAvC[1] - 12.*dAvC[2] + 1.*dAvC[3]);
1862  cumulant8[3] = (1./(35.*pow(fR0,8.))) * (56.*dAvC[0] - 28.*dAvC[1] + 8.*dAvC[2] - 1.*dAvC[3]);
1863  */
1864
1865  /*
1866  cout<<endl;
1867  cout<<"*********************************"<<endl;
1868  cout<<"cumulants8:"<<endl;
1869  cout<<" c_"<<cFlow<<"{2} = "<<cumulant8[0]<<endl;
1870  cout<<" c_"<<cFlow<<"{4} = "<<cumulant8[1]<<endl;
1871  cout<<" c_"<<cFlow<<"{6} = "<<cumulant8[2]<<endl;
1872  cout<<" c_"<<cFlow<<"{8} = "<<cumulant8[3]<<endl;
1873  cout<<endl;
1874  */
1875
1876  Double_t dV2o8=0.,dV4o8=0.,dV6o8=0.,dV8o8=0.;
1877
1878  if(cumulant8[0]>=0.)
1879  {
1880   dV2o8 = pow(cumulant8[0],(1./2.));
1881  }
1882  if(cumulant8[1]<=0.)
1883  {
1884   dV4o8 = pow(-cumulant8[1],(1./4.));
1885  }
1886  if(cumulant8[2]>=0.)
1887  {
1888   dV6o8 = pow((1./4.)*cumulant8[2],(1./6.));
1889  }
1890  if(cumulant8[3]<=0.)
1891  {
1892   dV8o8 = pow(-(1./33.)*cumulant8[3],(1./8.));
1893  }
1894
1895  cout<<endl;
1896  cout<<"***********************************"<<endl;
1897  cout<<"***********************************"<<endl;
1898  cout<<"flow estimates from GF-cumulants:"<<endl;
1899  cout<<"  (calculated up to 8th order)   "<<endl;
1900  cout<<endl;
1901
1902  Double_t sdQo8[4]={0.};
1903  Double_t chiQo8[4]={0.};
1904
1905  //v_2{2}
1906  if(nEvents8 && dAvM8 && cumulant8[0]>=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(cumulant8[0],(1./2.))*dAvM8,2.)>0.))
1907  {
1908   chiQo8[0]=dAvM8*dV2o8/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2o8*dAvM8,2.),0.5);
1909   if(chiQo8[0])
1910   {
1911    sdQo8[0]=pow(((1./(2.*dAvM8*nEvents8))*((1.+2.*pow(chiQo8[0],2.))/(2.*pow(chiQo8[0],2)))),0.5);
1912   }
1913   cout<<" v_"<<cFlow<<"{2} = "<<dV2o8<<" +/- "<<sdQo8[0]<<endl;
1914   //cout<<" v_"<<cFlow<<"{2} = "<<dV2o8<<" +/- "<<sdQo8[0]<<", chi{2} = "<<chiQo8[0]<<endl;//printing also the chi
1915  }
1916  else
1917  {
1918   cout<<" v_"<<cFlow<<"{2} = Im"<<endl;
1919  }
1920
1921  //v_2{4}
1922  if(nEvents8 && dAvM8 && cumulant8[1]<=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(-cumulant8[1],(1./4.))*dAvM8,2.)>0.))
1923  {
1924   chiQo8[1]=dAvM8*dV4o8/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4o8*dAvM8,2.),0.5);
1925   if(chiQo8[1])
1926   {
1927    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);
1928   }
1929   cout<<" v_"<<cFlow<<"{4} = "<<dV4o8<<" +/- "<<sdQo8[1]<<endl;
1930   //cout<<" v_"<<cFlow<<"{4} = "<<dV4o8<<" +/- "<<sdQo8[1]<<", chi{4} = "<<chiQo8[1]<<endl;//printing also the chi
1931  }
1932  else
1933  {
1934   cout<<" v_"<<cFlow<<"{4} = Im"<<endl;
1935  }
1936
1937  //v_2{6}
1938  if(nEvents8 && dAvM8 && cumulant8[2]>=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow((1./4.)*cumulant8[2],(1./6.))*dAvM8,2.)>0.))
1939  {
1940   chiQo8[2]=dAvM8*dV6o8/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV6o8*dAvM8,2.),0.5);
1941   if(chiQo8[2])
1942   {
1943    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);
1944   }
1945   cout<<" v_"<<cFlow<<"{6} = "<<dV6o8<<" +/- "<<sdQo8[2]<<endl;
1946   //cout<<" v_"<<cFlow<<"{6} = "<<dV6o8<<" +/- "<<sdQo8[2]<<", chi{6} = "<<chiQo8[2]<<endl;//printing also the chi
1947  }
1948  else
1949  {
1950   cout<<" v_"<<cFlow<<"{6} = Im"<<endl;
1951  }
1952
1953  //v_2{8}
1954  if(nEvents8 && dAvM8 && cumulant8[3]<=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(-(1./33.)*cumulant8[3],(1./8.))*dAvM8,2.)>0.))
1955  {
1956   chiQo8[3]=dAvM8*dV8o8/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV8o8*dAvM8,2.),0.5);
1957   if(chiQo8[3])
1958   {
1959    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);
1960   }
1961   cout<<" v_"<<cFlow<<"{8} = "<<dV8o8<<" +/- "<<sdQo8[3]<<endl;
1962   //cout<<" v_"<<cFlow<<"{8} = "<<dV8o8<<" +/- "<<sdQo8[3]<<", chi{8} = "<<chiQo8[3]<<endl;//printing also the chi
1963  }
1964  else
1965  {
1966   cout<<" v_"<<cFlow<<"{8} = Im"<<endl;
1967  }
1968
1969  cout<<endl;
1970  cout<<"   nEvts = "<<nEvents8<<", AvM = "<<dAvM8<<endl;
1971  cout<<"*********************************"<<endl;
1972
1973  //==============================================================================================================================================
1974
1975  //up to 16-th order:
1976  //avarage selected multiplicity
1977  Double_t dAvM16=0.;
1978  if(fAvMult16)
1979  {
1980   dAvM16=fAvMult16->GetBinContent(1);
1981  }
1982
1983  //number of events
1984  Int_t nEvents16=0;
1985  if(fAvMult16)
1986  {
1987   nEvents16=(Int_t)(fAvMult16->GetBinEntries(1));
1988  }
1989
1990  TMatrixD dAvG16(cPmax16,cQmax16);
1991  Bool_t someAvGEntryIsNegative16=kFALSE;
1992  for(Int_t p=0;p<cPmax16;p++)
1993  {
1994   for(Int_t q=0;q<cQmax16;q++)
1995   {
1996    dAvG16(p,q)=fIntGenFun16->GetBinContent(p+1,q+1);
1997    if(dAvG16(p,q)<0.)
1998    {
1999     someAvGEntryIsNegative16=kTRUE;
2000    }
2001   }
2002  }
2003
2004  /////////////////////////////////////////////////////////////////////////////
2005  //////////////////gen. function for the cumulants////////////////////////////
2006  /////////////////////////////////////////////////////////////////////////////
2007
2008  TMatrixD dC16(cPmax16,cQmax16);//C16[p][q]
2009  if(dAvM16>0 && someAvGEntryIsNegative16==kFALSE)
2010  {
2011   for(Int_t p=0;p<cPmax16;p++)
2012   {
2013    for(Int_t q=0;q<cQmax16;q++)
2014    {
2015     dC16(p,q)=1.*dAvM16*(pow(dAvG16(p,q),(1./dAvM16))-1.);
2016    }
2017   }
2018  }
2019
2020  /////////////////////////////////////////////////////////////////////////////
2021  ///////avaraging the gen. function for the cumulants over azimuth////////////
2022  /////////////////////////////////////////////////////////////////////////////
2023
2024  TVectorD dAvC16(cPmax16);//<C16[p][q]>
2025  Double_t tempHere16=0.;
2026  for (Int_t p=0;p<cPmax16;p++)
2027  {
2028   tempHere16=0.;
2029   for (Int_t q=0;q<cQmax16;q++)
2030   {
2031    tempHere16+=1.*dC16(p,q);
2032   }
2033   dAvC16[p]=1.*tempHere16/cQmax16;
2034  }
2035
2036  /////////////////////////////////////////////////////////////////////////////
2037  //////////////////////////////////final results//////////////////////////////
2038  /////////////////////////////////////////////////////////////////////////////
2039
2040  TVectorD cumulant16(cPmax16);//array to store various order cumulants
2041
2042  cumulant16[0] = (1./(fR0*fR0)) * (8.*dAvC16[0] - 14.*dAvC16[1] + (56./3.)*dAvC16[2] - (35./2.)*dAvC16[3] +
2043                               (56./5.)*dAvC16[4] - (14./3.)*dAvC16[5] + (8./7.)*dAvC16[6] - (1./8.)*dAvC16[7]);
2044
2045  cumulant16[1] = (1./pow(fR0,4.)) * ((-1924./35.)*dAvC16[0] + (621./5.)*dAvC16[1] - (8012./45.)*dAvC16[2] +
2046                                  (691./4.)*dAvC16[3] - (564./5.)*dAvC16[4] + (2143./45.)*dAvC16[5] -
2047                                  (412./35.)*dAvC16[6] + (363./280.)*dAvC16[7]);
2048
2049  cumulant16[2] = (1./pow(fR0,6.)) * (349.*dAvC16[0] - (18353./20.)*dAvC16[1] + (7173./5.)*dAvC16[2] -
2050                                  1457.*dAvC16[3] + (4891./5.)*dAvC16[4] - (1683./4.)*dAvC16[5] +
2051                                  (527./5.)*dAvC16[6] - (469./40.)*dAvC16[7]);
2052
2053  cumulant16[3] = (1./pow(fR0,8.)) * ((-10528./5.)*dAvC16[0] + (30578./5.)*dAvC16[1] - (51456./5.)*dAvC16[2] +
2054                                  10993.*dAvC16[3] - (38176./5.)*dAvC16[4] + (16818./5.)*dAvC16[5] -
2055                                  (4288./5.)*dAvC16[6] + (967./10.)*dAvC16[7]);
2056
2057  cumulant16[4] = (1./pow(fR0,10.)) * (11500.*dAvC16[0] - 35800.*dAvC16[1] + 63900.*dAvC16[2] - 71600.*dAvC16[3] +
2058                                   51620.*dAvC16[4] - 23400.*dAvC16[5] + 6100.*dAvC16[6] - 700.*dAvC16[7]);
2059
2060  cumulant16[5] = (1./pow(fR0,12.)) * (-52560.*dAvC16[0] + 172080.*dAvC16[1] - 321840.*dAvC16[2] + 376200.*dAvC16[3] -
2061                                   281520.*dAvC16[4] + 131760.*dAvC16[5] - 35280.*dAvC16[6] + 4140.*dAvC16[7]);
2062
2063  cumulant16[6] = (1./pow(fR0,14.)) * (176400.*dAvC16[0] - 599760.*dAvC16[1] + 1164240.*dAvC16[2] - 1411200.*dAvC16[3] +
2064                                   1093680.*dAvC16[4] - 529200.*dAvC16[5] + 146160.*dAvC16[6] - 17640.*dAvC16[7]);
2065
2066  cumulant16[7] = (1./pow(fR0,16.)) * (-322560*dAvC16[0] + 1128960.*dAvC16[1] - 2257920.*dAvC16[2] + 2822400.*dAvC16[3] -
2067                                   2257920.*dAvC16[4] + 1128960.*dAvC16[5] - 322560.*dAvC16[6] + 40320.*dAvC16[7]);
2068
2069  /*
2070  cout<<endl;
2071  cout<<"*********************************"<<endl;
2072  cout<<"cumulants:"<<endl;
2073  cout<<" c_"<<cFlow<<"{2} = "<<cumulant16[0]<<endl;
2074  cout<<" c_"<<cFlow<<"{4} = "<<cumulant16[1]<<endl;
2075  cout<<" c_"<<cFlow<<"{6} = "<<cumulant16[2]<<endl;
2076  cout<<" c_"<<cFlow<<"{8} = "<<cumulant16[3]<<endl;
2077  cout<<"c_"<<cFlow<<"{10} = "<<cumulant16[4]<<endl;
2078  cout<<"c_"<<cFlow<<"{12} = "<<cumulant16[5]<<endl;
2079  cout<<"c_"<<cFlow<<"{14} = "<<cumulant16[6]<<endl;
2080  cout<<"c_"<<cFlow<<"{16} = "<<cumulant16[7]<<endl;
2081  cout<<endl;
2082  */
2083
2084  Double_t dV2o16=0.,dV4o16=0.,dV6o16=0.,dV8o16=0.,V10o16=0.,V12o16=0.,V14o16=0.,V16o16=0.;
2085
2086  if(cumulant16[0]>=0.)
2087  {
2088   dV2o16 = pow(cumulant16[0],(1./2.));
2089  }
2090  if(cumulant16[1]<=0.)
2091  {
2092   dV4o16 = pow(-cumulant16[1],(1./4.));
2093  }
2094  if(cumulant16[2]>=0.)
2095  {
2096   dV6o16 = pow((1./4.)*cumulant16[2],(1./6.));
2097  }
2098  if(cumulant16[3]<=0.)
2099  {
2100   dV8o16 = pow(-(1./33.)*cumulant16[3],(1./8.));
2101  }
2102  if(cumulant16[4]>=0.)
2103  {
2104   V10o16 = pow((1./456.)*cumulant16[4],(1./10.));
2105  }
2106  if(cumulant16[5]<=0.)
2107  {
2108   V12o16 = pow(-(1./9460.)*cumulant16[5],(1./12.));
2109  }
2110  if(cumulant16[6]>=0.)
2111  {
2112   V14o16 = pow((1./274800.)*cumulant16[6],(1./14.));
2113  }
2114  if(cumulant16[7]<=0.)
2115  {
2116   V16o16 = pow(-(1./10643745.)*cumulant16[7],(1./16.));
2117  }
2118
2119  cout<<endl;
2120  cout<<"***********************************"<<endl;
2121  cout<<"***********************************"<<endl;
2122  cout<<"flow estimates from GF-cumulants:"<<endl;
2123  cout<<"  (calculated up to 16th order)   "<<endl;
2124  cout<<endl;
2125
2126  Double_t sdQo16[8]={0.};
2127  Double_t chiQo16[8]={0.};
2128
2129  //v_2{2}
2130  if(nEvents16 && dAvM16 && cumulant16[0]>=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(cumulant16[0],(1./2.))*dAvM16,2.)>0.))
2131  {
2132   chiQo16[0]=dAvM16*dV2o16/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2o16*dAvM16,2.),0.5);
2133   if(chiQo16[0])
2134   {
2135    sdQo16[0]=pow(((1./(2.*dAvM16*nEvents16))*((1.+2.*pow(chiQo16[0],2))/(2.*pow(chiQo16[0],2)))),0.5);
2136   }
2137   cout<<" v_"<<cFlow<<"{2} = "<<dV2o16<<" +/- "<<sdQo16[0]<<endl;
2138   //cout<<" v_"<<cFlow<<"{2} = "<<dV2o16<<" +/- "<<sdQo16[0]<<", chi{2} = "<<chiQo16[0]<<endl;//printing also the chi
2139  }
2140  else
2141  {
2142   cout<<" v_"<<cFlow<<"{2} = Im"<<endl;
2143  }
2144
2145  //v_2{4}
2146  if(nEvents16 && dAvM16 && cumulant16[1]<=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(-cumulant16[1],(1./4.))*dAvM16,2.)>0.))
2147  {
2148   chiQo16[1]=dAvM16*dV4o16/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4o16*dAvM16,2.),0.5);
2149   if(chiQo16[1])
2150   {
2151    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);
2152   }
2153   cout<<" v_"<<cFlow<<"{4} = "<<dV4o16<<" +/- "<<sdQo16[1]<<endl;
2154   //cout<<" v_"<<cFlow<<"{4} = "<<dV4o16<<" +/- "<<sdQo16[1]<<", chi{4} = "<<chiQo16[1]<<endl;//printing also the chi
2155  }
2156  else
2157  {
2158   cout<<" v_"<<cFlow<<"{4} = Im"<<endl;
2159  }
2160
2161  //v_2{6}
2162  if(nEvents16 && dAvM16 && cumulant16[2]>=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow((1./4.)*cumulant16[2],(1./6.))*dAvM16,2.)>0.))
2163  {
2164   chiQo16[2]=dAvM16*dV6o16/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV6o16*dAvM16,2.),0.5);
2165   if(chiQo16[2])
2166   {
2167    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);
2168   }
2169   cout<<" v_"<<cFlow<<"{6} = "<<dV6o16<<" +/- "<<sdQo16[2]<<endl;
2170   //cout<<" v_"<<cFlow<<"{6} = "<<dV6o16<<" +/- "<<sdQo16[2]<<", chi{6} = "<<chiQo16[2]<<endl;//printing also the chi
2171  }
2172  else
2173  {
2174   cout<<" v_"<<cFlow<<"{6} = Im"<<endl;
2175  }
2176
2177  //v_2{8}
2178  if(nEvents16 && dAvM16 && cumulant16[3]<=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(-(1./33.)*cumulant16[3],(1./8.))*dAvM16,2.)>0.))
2179  {
2180   chiQo16[3]=dAvM16*dV8o16/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV8o16*dAvM16,2.),0.5);
2181   if(chiQo16[3])
2182   {
2183    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);
2184   }
2185   cout<<" v_"<<cFlow<<"{8} = "<<dV8o16<<" +/- "<<sdQo16[3]<<endl;
2186   //cout<<" v_"<<cFlow<<"{8} = "<<dV8o16<<" +/- "<<sdQo16[3]<<", chi{8} = "<<chiQo16[3]<<endl;//printing also the chi
2187  }
2188  else
2189  {
2190   cout<<" v_"<<cFlow<<"{8} = Im"<<endl;
2191  }
2192
2193  //v_2{10}
2194  if(nEvents16 && dAvM16 && cumulant16[4]>=0.)
2195  {
2196   cout<<"v_"<<cFlow<<"{10} = "<<V10o16<<endl;
2197  }
2198  else
2199  {
2200   cout<<"v_"<<cFlow<<"{10} = Im"<<endl;
2201  }
2202
2203  //v_2{12}
2204  if(nEvents16 && dAvM16 && cumulant16[5]<=0.)
2205  {
2206   cout<<"v_"<<cFlow<<"{12} = "<<V12o16<<endl;
2207  }
2208  else
2209  {
2210   cout<<"v_"<<cFlow<<"{12} = Im"<<endl;
2211  }
2212
2213  //v_2{14}
2214  if(nEvents16 && dAvM16 && cumulant16[6]>=0.)
2215  {
2216   cout<<"v_"<<cFlow<<"{14} = "<<V14o16<<endl;
2217  }
2218  else
2219  {
2220   cout<<"v_"<<cFlow<<"{14} = Im"<<endl;
2221  }
2222
2223  //v_2{16}
2224  if(nEvents16 && dAvM16 && cumulant16[7]<=0.)
2225  {
2226   cout<<"v_"<<cFlow<<"{16} = "<<V16o16<<endl;
2227  }
2228  else
2229  {
2230   cout<<"v_"<<cFlow<<"{16} = Im"<<endl;
2231  }
2232
2233  cout<<endl;
2234  cout<<"   nEvts = "<<nEvents16<<", AvM = "<<dAvM16<<endl;
2235  cout<<"*********************************"<<endl;
2236
2237  //==============================================================================================================================================
2238
2239  }//end of if(otherEquations)
2240 }
2241
2242
2243