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