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