]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FLOW/AliFlowCommon/AliCumulantsFunctions.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
small mods for separate task
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliCumulantsFunctions.cxx
CommitLineData
2188af53 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 *********************************/
aaebd73d 24
25#define AliCumulantsFunctions_cxx
26
27#include "Riostream.h"
aaebd73d 28#include "TChain.h"
29#include "TFile.h"
2188af53 30#include "TList.h"
aaebd73d 31#include "TParticle.h"
32
33#include "TProfile.h"
34#include "TProfile2D.h"
35#include "TProfile3D.h"
2188af53 36#include "TH1.h"
37#include "TH1D.h"
38
aaebd73d 39#include "AliFlowEventSimple.h"
40#include "AliFlowTrackSimple.h"
41#include "AliFlowAnalysisWithCumulants.h"
42#include "AliFlowCumuConstants.h"
43#include "AliFlowCommonConstants.h"
813a4157 44#include "AliFlowCommonHist.h"
52021ae2 45#include "AliFlowCommonHistResults.h"
aaebd73d 46#include "AliCumulantsFunctions.h"
47
aaebd73d 48ClassImp(AliCumulantsFunctions)
49
2188af53 50//================================================================================================================_
aaebd73d 51
52AliCumulantsFunctions::AliCumulantsFunctions():
2188af53 53 fIntGenFun(NULL),
cb308e83 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.
813a4157 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),
2188af53 74 fifr(NULL),
75 fdfr2(NULL),
76 fdfr4(NULL),
77 fdfr6(NULL),
78 fdfr8(NULL),
79 fAvMult(NULL),
80 fQVector(NULL),
9dd53ff2 81 fAverageOfSquaredWeight(NULL),
52021ae2 82 fchr2nd(NULL),
83 fchr4th(NULL),
84 fchr6th(NULL),
813a4157 85 fchr8th(NULL),
86 fch(NULL)//common control histograms
52021ae2 87 /*
2188af53 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)
52021ae2 104 */
2188af53 105{
106 //default constructor
107}
aaebd73d 108
2188af53 109AliCumulantsFunctions::~AliCumulantsFunctions()
110{
111 //destructor
aaebd73d 112}
113
9dd53ff2 114AliCumulantsFunctions::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):
813a4157 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),
2188af53 136 fifr(ifr),
137 fdfr2(dfr2),
138 fdfr4(dfr4),
139 fdfr6(dfr6),
140 fdfr8(dfr8),
813a4157 141 fAvMult(avMult),
142 fQVector(qVector),
9dd53ff2 143 fAverageOfSquaredWeight(averageOfSquaredWeight),
52021ae2 144 fchr2nd(chr2nd),
145 fchr4th(chr4th),
146 fchr6th(chr6th),
813a4157 147 fchr8th(chr8th),
148 fch(ch)//common control histograms
52021ae2 149 /*
2188af53 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)
52021ae2 166 */
2188af53 167{
168 //custom constructor
169}
aaebd73d 170
2188af53 171//================================================================================================================
172
173void AliCumulantsFunctions::Calculate()
174{
175 //calculate cumulants and final integrated and differential flow estimates and store the results into output histograms
b6cd16a9 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
2188af53 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)
813a4157 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
52021ae2 191
813a4157 192 Double_t fR0=AliFlowCumuConstants::fgR0; //needed for numerics
cb308e83 193 //Double_t fPtMax=AliFlowCommonConstants::GetPtMax(); //maximum pt
194 //Double_t fPtMin=AliFlowCommonConstants::GetPtMin(); //minimum pt
813a4157 195 //Double_t fBinWidthPt=(fPtMax-fPtMin)/fgknBinsPt; //width of pt bin (in GeV)
b6cd16a9 196
197 Bool_t fOtherEquations=AliFlowCumuConstants::fgOtherEquations; //numerical equations for cumulants solved up to different highest order
2188af53 198
b6cd16a9 199 //avarage selected multiplicity
813a4157 200 Double_t dAvM=0.;
b6cd16a9 201 if(fAvMult)
202 {
813a4157 203 dAvM=fAvMult->GetBinContent(1);
b6cd16a9 204 }
205
52021ae2 206 //number of events
b6cd16a9 207 Int_t nEvents=0;
208 if(fAvMult)
209 {
210 nEvents=(Int_t)(fAvMult->GetBinEntries(1));
211 }
aaebd73d 212
2188af53 213 //<Q-vector stuff>
813a4157 214 Double_t dAvQx=0.,dAvQy=0.,dAvQ2x=0.,dAvQ2y=0.;
b6cd16a9 215 if(fQVector)
216 {
813a4157 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>
b6cd16a9 221 }
222
9dd53ff2 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
b6cd16a9 231 //<G[p][q]>
813a4157 232 Double_t dAvG[fgkPmax][fgkQmax]={{0.}};
cb308e83 233 Bool_t someAvGEntryIsNegative=kFALSE;
b6cd16a9 234 if(fIntGenFun)
235 {
236 for(Int_t p=0;p<fgkPmax;p++)
237 {
238 for(Int_t q=0;q<fgkQmax;q++)
239 {
813a4157 240 dAvG[p][q]=fIntGenFun->GetBinContent(p+1,q+1);
241 if(dAvG[p][q]<0.)
cb308e83 242 {
243 someAvGEntryIsNegative=kTRUE;
244 }
b6cd16a9 245 }
246 }
cb308e83 247 }
248
2188af53 249 /////////////////////////////////////////////////////////////////////////////
250 //////////////////gen. function for the cumulants////////////////////////////
251 /////////////////////////////////////////////////////////////////////////////
cb308e83 252
813a4157 253 Double_t dC[fgkPmax][fgkQmax]={{0.}};//C[p][q]
254 if(dAvM>0 && someAvGEntryIsNegative==kFALSE)
b6cd16a9 255 {
256 for(Int_t p=0;p<fgkPmax;p++)
257 {
258 for(Int_t q=0;q<fgkQmax;q++)
259 {
813a4157 260 dC[p][q]=1.*dAvM*(pow(dAvG[p][q],(1./dAvM))-1.);
6d4fa5d3 261 }
2188af53 262 }
263 }
cb308e83 264
2188af53 265 /////////////////////////////////////////////////////////////////////////////
266 ///////avaraging the gen. function for the cumulants over azimuth////////////
267 /////////////////////////////////////////////////////////////////////////////
aaebd73d 268
813a4157 269 Double_t dAvC[fgkPmax]={0.};//<C[p][q]>
270 Double_t tempHere=0.;
b6cd16a9 271 for(Int_t p=0;p<fgkPmax;p++)
272 {
813a4157 273 tempHere=0.;
b6cd16a9 274 for(Int_t q=0;q<fgkQmax;q++)
275 {
813a4157 276 tempHere+=1.*dC[p][q];
2188af53 277 }
813a4157 278 dAvC[p]=1.*tempHere/fgkQmax;
2188af53 279 }
280
281 /////////////////////////////////////////////////////////////////////////////
282 //////////////////////////////////final results//////////////////////////////
283 /////////////////////////////////////////////////////////////////////////////
284
285 Double_t cumulant[fgkPmax];//array to store various order cumulants
b6cd16a9 286
287 //system of eq. for the cumulants
813a4157 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]);
2188af53 293
b6cd16a9 294 /*
295 cout<<endl;
5b04fe11 296 cout<<"*********************************"<<endl;
b6cd16a9 297 cout<<"cumulants:"<<endl;
2188af53 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;
b6cd16a9 302 cout<<"c_"<<fgkFlow<<"{10} = "<<cumulant[4]<<endl;
303 cout<<endl;
304 */
2188af53 305
813a4157 306 Double_t dV2=0.,dV4=0.,dV6=0.,dV8=0.,dV10=0.;//integrated flow estimates
2188af53 307
b6cd16a9 308 if(cumulant[0]>=0.)
309 {
813a4157 310 dV2=pow(cumulant[0],(1./2.));
2188af53 311 }
b6cd16a9 312 if(cumulant[1]<=0.)
313 {
813a4157 314 dV4=pow(-cumulant[1],(1./4.));
2188af53 315 }
b6cd16a9 316 if(cumulant[2]>=0.)
317 {
813a4157 318 dV6=pow((1./4.)*cumulant[2],(1./6.));
2188af53 319 }
b6cd16a9 320 if(cumulant[3]<=0.)
321 {
813a4157 322 dV8=pow(-(1./33.)*cumulant[3],(1./8.));
2188af53 323 }
b6cd16a9 324 if(cumulant[4]>=0.)
325 {
813a4157 326 dV10=pow((1./456.)*cumulant[4],(1./10.));
2188af53 327 }
b6cd16a9 328
329 cout<<endl;
813a4157 330 cout<<"***************************************"<<endl;
331 cout<<"***************************************"<<endl;
9dd53ff2 332 cout<<"flow estimates from GF-cumulants: "<<endl;
b6cd16a9 333 cout<<endl;
334
813a4157 335 Double_t sdQ[4]={0.};
336 Double_t chiQ[4]={0.};
b6cd16a9 337
338 //v_2{2}
9dd53ff2 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.))
b6cd16a9 340 {
9dd53ff2 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
813a4157 342 if(chiQ[0])
cb308e83 343 {
813a4157 344 sdQ[0]=pow(((1./(2.*dAvM*nEvents))*((1.+2.*pow(chiQ[0],2))/(2.*pow(chiQ[0],2)))),0.5);
cb308e83 345 }
813a4157 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]);
b6cd16a9 350 //filling common histograms:
813a4157 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
b6cd16a9 360 }
361 else
362 {
813a4157 363 cout<<" v_"<<fgkFlow<<"{2} = Im"<<endl;
b6cd16a9 364 }
2188af53 365
b6cd16a9 366 //v_2{4}
9dd53ff2 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.))
b6cd16a9 368 {
9dd53ff2 369 chiQ[1]=(dAvM/dAvw2)*dV4/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4*(dAvM/dAvw2),2.),0.5);
813a4157 370 if(chiQ[1])
cb308e83 371 {
813a4157 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);
cb308e83 373 }
813a4157 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]);
b6cd16a9 378 //filling common histograms:
813a4157 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
b6cd16a9 388 }
389 else
390 {
813a4157 391 cout<<" v_"<<fgkFlow<<"{4} = Im"<<endl;
b6cd16a9 392 }
2188af53 393
b6cd16a9 394 //v_2{6}
9dd53ff2 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.))
b6cd16a9 396 {
9dd53ff2 397 chiQ[2]=(dAvM/dAvw2)*dV6/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV6*(dAvM/dAvw2),2.),0.5);
813a4157 398 if(chiQ[2])
cb308e83 399 {
813a4157 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);
cb308e83 401 }
813a4157 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]);
b6cd16a9 406 //filling common histograms:
813a4157 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
b6cd16a9 417 }
418 else
419 {
813a4157 420 cout<<" v_"<<fgkFlow<<"{6} = Im"<<endl;
b6cd16a9 421 }
2188af53 422
b6cd16a9 423 //v_2{8}
9dd53ff2 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.))
b6cd16a9 425 {
9dd53ff2 426 chiQ[3]=(dAvM/dAvw2)*dV8/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV8*(dAvM/dAvw2),2.),0.5);
813a4157 427 if(chiQ[3])
cb308e83 428 {
813a4157 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);
cb308e83 430 }
813a4157 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]);
b6cd16a9 435 //filling common histograms:
813a4157 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
b6cd16a9 444 }
445 else
446 {
813a4157 447 cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl;
b6cd16a9 448 }
2188af53 449
b6cd16a9 450 /*
451 //v_2{10}
813a4157 452 if (dAvM && cumulant[4]>=0.)
b6cd16a9 453 {
813a4157 454 cout<<"v_"<<fgkFlow<<"{10} = "<<dV10<<endl;
b6cd16a9 455 }
456 else
457 {
458 cout<<"v_"<<fgkFlow<<"{10} = Im"<<endl;
459 }
460 */
461
462 cout<<endl;
813a4157 463 cout<<" nEvts = "<<nEvents<<", AvM = "<<dAvM<<endl;
464 cout<<"***************************************"<<endl;
465 cout<<"***************************************"<<endl;
aaebd73d 466
813a4157 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);
9dd53ff2 584 if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2*(dAvM/dAvw2),2.)>0. && ptBinRPNoOfParticles[b]>0.)
813a4157 585 {
9dd53ff2 586 sdRPDiff2pt[b]=pow((1./(2.*ptBinRPNoOfParticles[b]))*((1.+pow(chiQ[0],2.))/pow(chiQ[0],2.)),0.5);
813a4157 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);
9dd53ff2 610 if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4*(dAvM/dAvw2),2.)>0.&&ptBinRPNoOfParticles[b]>0.) // to be improved
2f5e6dea 611 if (ptBinRPNoOfParticles[b]>0.)
813a4157 612 {
9dd53ff2 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);
813a4157 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
9dd53ff2 681 ///////////////////////////////////////////////////////////////////////////////////
682 //////////////////////// INTEGRATED FLOW CALCULATIONS (RP) /////////////////////////
683 ///////////////////////////////////////////////////////////////////////////////////
813a4157 684
9dd53ff2 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 }
813a4157 715
9dd53ff2 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 }
813a4157 731
9dd53ff2 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
813a4157 750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
b6cd16a9 773
774
775
813a4157 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);
9dd53ff2 830 if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2*(dAvM/dAvw2),2.)>0. && etaBinRPNoOfParticles[b]>0.) // to be improved
813a4157 831 {
9dd53ff2 832 sdRPDiff2eta[b]=pow((1./(2.*etaBinRPNoOfParticles[b]))*((1.+pow(chiQ[0],2.))/pow(chiQ[0],2.)),0.5);
813a4157 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);
9dd53ff2 856 if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4*(dAvM/dAvw2),2.)>0.&&etaBinRPNoOfParticles[b]>0.) // to be improved
813a4157 857 {
9dd53ff2 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);
813a4157 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 //===========================================================================
b6cd16a9 970
2188af53 971 /////////////////////////////////////////////////////////////////////////////
972 ///////////////////////DIFFERENTIAL FLOW CALCULATIONS////////////////////////
973 /////////////////////////////////////////////////////////////////////////////
aaebd73d 974
813a4157 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++)
b6cd16a9 981 {
4e1fb6e8 982 ptBinPOINoOfParticles[b]=fPtBinPOINoOfParticles->GetBinEntries(b+1);
b6cd16a9 983 for(Int_t p=0;p<fgkPmax;p++)
984 {
985 for(Int_t q=0;q<fgkQmax;q++)
986 {
813a4157 987 if(dAvG[p][q])
b6cd16a9 988 {
813a4157 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 }
b6cd16a9 997 }
998 }
999 }
813a4157 1000 }
52021ae2 1001
813a4157 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 {
4e1fb6e8 1009 etaBinPOINoOfParticles[b]=fEtaBinPOINoOfParticles->GetBinEntries(b+1);
813a4157 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
52021ae2 1029 /*
813a4157 1030 if(dAvM){
1031 for(Int_t b=0;b<fgknBinsPt;b++){cout<<"***************************************"<<endl;
1032 cout<<"***************************************"<<endl;
2188af53 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 }
6d4fa5d3 1061 }
52021ae2 1062 */
aaebd73d 1063
813a4157 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.};
aaebd73d 1101
813a4157 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);
9dd53ff2 1117 if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2*(dAvM/dAvw2),2.)>0.&&ptBinPOINoOfParticles[b]>0.)
813a4157 1118 {
9dd53ff2 1119 sdDiff2pt[b]=pow((1./(2.*ptBinPOINoOfParticles[b]))*((1.+pow(chiQ[0],2.))/pow(chiQ[0],2.)),0.5);
813a4157 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);
9dd53ff2 1143 if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4*(dAvM/dAvw2),2.)>0.&&ptBinPOINoOfParticles[b]>0.)
813a4157 1144 {
9dd53ff2 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);
813a4157 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;
aaebd73d 1162 }
813a4157 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
9dd53ff2 1215 ///////////////////////////////////////////////////////////////////////////////////
1216 //////////////////////// INTEGRATED FLOW CALCULATIONS (POI) ///////////////////////
1217 ///////////////////////////////////////////////////////////////////////////////////
813a4157 1218
1219 Double_t dV2POI=0., dV4POI=0., dV6POI=0., dV8POI=0.;
1220 Double_t dV2POIError=0., dV4POIError=0., dV6POIError=0., dV8POIError=0.;
9dd53ff2 1221 Double_t dSumOfYieldPOI=0.;
4e1fb6e8 1222 for (Int_t b=1;b<fgknBinsPt+1;b++)
813a4157 1223 {
1224 if(fch->GetHistPtDiff())
1225 {
9dd53ff2 1226 dSumOfYieldPOI+=(fch->GetHistPtDiff())->GetBinContent(b);
813a4157 1227 if(fchr2nd->GetHistDiffFlowPtPOI())
1228 {
1229 dV2POI+=((fchr2nd->GetHistDiffFlowPtPOI())->GetBinContent(b))*(fch->GetHistPtDiff())->GetBinContent(b);
4e1fb6e8 1230 dV2POIError+=pow((fch->GetHistPtDiff())->GetBinContent(b),2.)*pow((fchr2nd->GetHistDiffFlowPtPOI())->GetBinError(b),2.);
813a4157 1231 }
1232 if(fchr4th->GetHistDiffFlowPtPOI())
1233 {
1234 dV4POI+=((fchr4th->GetHistDiffFlowPtPOI())->GetBinContent(b))*(fch->GetHistPtDiff())->GetBinContent(b);
4e1fb6e8 1235 dV4POIError+=pow((fch->GetHistPtDiff())->GetBinContent(b),2.)*pow((fchr4th->GetHistDiffFlowPtPOI())->GetBinError(b),2.);
813a4157 1236 }
1237 if(fchr6th->GetHistDiffFlowPtPOI())
1238 {
1239 dV6POI+=((fchr6th->GetHistDiffFlowPtPOI())->GetBinContent(b))*(fch->GetHistPtDiff())->GetBinContent(b);
4e1fb6e8 1240 dV6POIError+=pow((fch->GetHistPtDiff())->GetBinContent(b),2.)*pow((fchr6th->GetHistDiffFlowPtPOI())->GetBinError(b),2.);
813a4157 1241 }
1242 if(fchr8th->GetHistDiffFlowPtPOI())
1243 {
1244 dV8POI+=((fchr8th->GetHistDiffFlowPtPOI())->GetBinContent(b))*(fch->GetHistPtDiff())->GetBinContent(b);
4e1fb6e8 1245 dV8POIError+=pow((fch->GetHistPtDiff())->GetBinContent(b),2.)*pow((fchr8th->GetHistDiffFlowPtPOI())->GetBinError(b),2.);
813a4157 1246 }
1247 }
1248 }
1249
9dd53ff2 1250 if(dSumOfYieldPOI)
813a4157 1251 {
9dd53ff2 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);
4e1fb6e8 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)
813a4157 1264 }
1265
1266 cout<<endl;
1267 cout<<"***************************************"<<endl;
1268 cout<<"***************************************"<<endl;
1269 cout<<"flow estimates from GF-cumulants (POI):"<<endl;
1270 cout<<endl;
1271
4e1fb6e8 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;
813a4157 1276
9dd53ff2 1277 cout<<endl;
1278 cout<<" nEvts = "<<(fch->GetHistMultDiff())->GetEntries()<<", AvM = "<<(fch->GetHistMultDiff())->GetMean()<<endl;
813a4157 1279 cout<<"***************************************"<<endl;
1280 cout<<"***************************************"<<endl;
9dd53ff2 1281 cout<<endl;
813a4157 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 }
aaebd73d 1298 }
1299
813a4157 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.};
aaebd73d 1305
813a4157 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]);
aaebd73d 1313 }
1314
813a4157 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)
aaebd73d 1326
813a4157 1327 //cout<<"number of eta bins: "<<fgknBinsEta<<endl;
2188af53 1328 //cout<<"****************************************"<<endl;
813a4157 1329 for (Int_t b=0;b<fgknBinsEta;b++){
1330 //cout<<"eta bin: "<<b*fBinWidthPt<<"-"<<(b+1)*fBinWidthPt<<" GeV"<<endl;
2188af53 1331
1332 //v'_{2/2}{2}
6d4fa5d3 1333 if(cumulant[0]>0)
2188af53 1334 {
813a4157 1335 v2eta[b]=etaDiffCumulant2[b]/pow(cumulant[0],0.5);
9dd53ff2 1336 if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2*(dAvM/dAvw2),2.)>0.&&etaBinPOINoOfParticles[b]>0.) // to be improved
2188af53 1337 {
9dd53ff2 1338 sdDiff2eta[b]=pow((1./(2.*etaBinPOINoOfParticles[b]))*((1.+pow(chiQ[0],2.))/pow(chiQ[0],2.)),0.5);
813a4157 1339 //cout<<"v'_2/2{2} = "<<v2eta[b]<<"%, "<<" "<<"sd{2} = "<<100.*sdDiff2eta[b]<<"%"<<endl;
9dd53ff2 1340 fdfr2->SetBinContent(b+1,v2eta[b]);
1341 fdfr2->SetBinError(b+1,sdDiff2eta[b]);
52021ae2 1342 //common histogram:
813a4157 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
2188af53 1351 } else {
aaebd73d 1352 //cout<<"v'_2/2{2} = Im"<<endl;
2188af53 1353 }
aaebd73d 1354 }else{
2188af53 1355 //cout<<"v'_2/2{2} = Im"<<endl;
aaebd73d 1356 }
2188af53 1357
1358 //v'_{2/2}{4}
6d4fa5d3 1359 if(cumulant[1]<0)
2188af53 1360 {
813a4157 1361 v4eta[b]=-etaDiffCumulant4[b]/pow(-cumulant[1],0.75);
9dd53ff2 1362 if (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4*(dAvM/dAvw2),2.)>0.&&etaBinPOINoOfParticles[b]>0.) // to be improved
2188af53 1363 {
9dd53ff2 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);
813a4157 1365 //cout<<"v'_2/2{4} = "<<v4eta[b]<<"%, "<<" "<<"sd{4} = "<<100.*sdDiff4eta[b]<<"%"<<endl;
9dd53ff2 1366 fdfr4->SetBinContent(b+1,v4eta[b]);
1367 fdfr4->SetBinError(b+1,sdDiff4eta[b]);
52021ae2 1368 //common histogram:
813a4157 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
2188af53 1379 } else {
1380 //cout<<"v'_2/2{4} = Im"<<endl;
1381 }
aaebd73d 1382 }else{
2188af53 1383 //cout<<"v'_2/2{4} = Im"<<endl;
aaebd73d 1384 }
2188af53 1385
1386 //v'_{2/2}{6}
6d4fa5d3 1387 if(cumulant[2]>0){
813a4157 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]);
52021ae2 1391 //common histogram:
813a4157 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
aaebd73d 1403 }else{
2188af53 1404 //cout<<"v'_2/2{6} = Im"<<endl;
aaebd73d 1405 }
2188af53 1406
1407 //v'_{2/2}{8}
6d4fa5d3 1408 if(cumulant[3]<0){
813a4157 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]);
52021ae2 1412 //common histogram:
813a4157 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
aaebd73d 1423 }else{
2188af53 1424 //cout<<"v'_2/2{8} = Im"<<endl;
aaebd73d 1425 }
2188af53 1426 //cout<<"****************************************"<<endl;
aaebd73d 1427 }
813a4157 1428 //-------------------------------------------------------------------------------------------------------------------------------
1429
1430
b6cd16a9 1431
cb308e83 1432
1433
1434
1435
1436
1437 //off the record: numerical equations for cumulants solved up to different highest order
b6cd16a9 1438 if(fOtherEquations)
1439 {
1440
1441 //==============================================================================================================================================
cb308e83 1442
b6cd16a9 1443 //up to 4th order
cb308e83 1444 //avarage selected multiplicity
813a4157 1445 Double_t dAvM4=0.;
cb308e83 1446 if(fAvMult4)
1447 {
813a4157 1448 dAvM4=fAvMult4->GetBinContent(1);
cb308e83 1449 }
1450
1451 //number of events
1452 Int_t nEvents4=0;
1453 if(fAvMult4)
1454 {
1455 nEvents4=(Int_t)(fAvMult4->GetBinEntries(1));
1456 }
1457
813a4157 1458 Double_t dAvG4[fgkPmax4][fgkQmax4]={{0.}};
cb308e83 1459 Bool_t someAvGEntryIsNegative4=kFALSE;
b6cd16a9 1460 for(Int_t p=0;p<fgkPmax4;p++)
1461 {
1462 for(Int_t q=0;q<fgkQmax4;q++)
1463 {
813a4157 1464 dAvG4[p][q]=fIntGenFun4->GetBinContent(p+1,q+1);
1465 if(dAvG4[p][q]<0.)
cb308e83 1466 {
1467 someAvGEntryIsNegative4=kTRUE;
1468 }
b6cd16a9 1469 }
1470 }
1471 /////////////////////////////////////////////////////////////////////////////
1472 //////////////////gen. function for the cumulants////////////////////////////
1473 /////////////////////////////////////////////////////////////////////////////
1474
813a4157 1475 Double_t dC4[fgkPmax4][fgkQmax4]={{0.}};//C[p][q]
1476 if(dAvM4>0 && someAvGEntryIsNegative4==kFALSE)
b6cd16a9 1477 {
1478 for (Int_t p=0;p<fgkPmax4;p++)
1479 {
1480 for (Int_t q=0;q<fgkQmax4;q++)
1481 {
813a4157 1482 dC4[p][q]=1.*dAvM4*(pow(dAvG4[p][q],(1./dAvM4))-1.);
b6cd16a9 1483 }
1484 }
1485 }
1486 /////////////////////////////////////////////////////////////////////////////
1487 ///////avaraging the gen. function for the cumulants over azimuth////////////
1488 /////////////////////////////////////////////////////////////////////////////
1489
813a4157 1490 Double_t dAvC4[fgkPmax4]={0.};//<C[p][q]>
b6cd16a9 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 {
813a4157 1496 tempHere4+=1.*dC4[p][q];
b6cd16a9 1497 }
813a4157 1498 dAvC4[p]=1.*tempHere4/fgkQmax4;
b6cd16a9 1499 }
1500
1501 /////////////////////////////////////////////////////////////////////////////
1502 //////////////////////////////////final results//////////////////////////////
1503 /////////////////////////////////////////////////////////////////////////////
1504
1505 Double_t cumulant4[fgkPmax4];//array to store various order cumulants
1506
813a4157 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]);
b6cd16a9 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
813a4157 1519 Double_t dV2o4=0.,dV4o4=0.;
b6cd16a9 1520
1521 if(cumulant4[0]>=0.)
1522 {
813a4157 1523 dV2o4 = pow(cumulant4[0],(1./2.));
b6cd16a9 1524 }
1525 if(cumulant4[1]<=0.)
1526 {
813a4157 1527 dV4o4 = pow(-cumulant4[1],(1./4.));
b6cd16a9 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
813a4157 1537 Double_t sdQo4[2]={0.};
1538 Double_t chiQo4[2]={0.};
b6cd16a9 1539
1540 //v_2{2}
813a4157 1541 if(nEvents4 && dAvM4 && cumulant4[0]>=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(cumulant4[0],(1./2.))*dAvM4,2.)>0.))
b6cd16a9 1542 {
813a4157 1543 chiQo4[0]=dAvM4*dV2o4/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2o4*dAvM4,2.),0.5);
1544 if(chiQo4[0])
cb308e83 1545 {
813a4157 1546 sdQo4[0]=pow(((1./(2.*dAvM4*nEvents4))*((1.+2.*pow(chiQo4[0],2))/(2.*pow(chiQo4[0],2)))),0.5);
cb308e83 1547 }
813a4157 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
b6cd16a9 1550 }
1551 else
1552 {
1553 cout<<" v_"<<fgkFlow<<"{2} = Im"<<endl;
1554 }
1555
1556 //v_2{4}
813a4157 1557 if(nEvents4 && dAvM4 && cumulant4[1]<=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(-cumulant4[1],(1./4.))*dAvM4,2.)>0.))
b6cd16a9 1558 {
813a4157 1559 chiQo4[1]=dAvM4*dV4o4/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4o4*dAvM4,2.),0.5);
1560 if(chiQo4[1])
cb308e83 1561 {
813a4157 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);
cb308e83 1563 }
813a4157 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
b6cd16a9 1566 }
1567 else
1568 {
1569 cout<<" v_"<<fgkFlow<<"{4} = Im"<<endl;
1570 }
1571
1572 cout<<endl;
813a4157 1573 cout<<" nEvts = "<<nEvents4<<", AvM = "<<dAvM4<<endl;
b6cd16a9 1574 cout<<"***********************************"<<endl;
1575 cout<<"***********************************"<<endl;
1576
1577 //==============================================================================================================================================
1578
1579 //up to 6th order
cb308e83 1580 //avarage selected multiplicity
813a4157 1581 Double_t dAvM6=0.;
cb308e83 1582 if(fAvMult6)
1583 {
813a4157 1584 dAvM6=fAvMult6->GetBinContent(1);
cb308e83 1585 }
1586
1587 //number of events
1588 Int_t nEvents6=0;
1589 if(fAvMult6)
1590 {
1591 nEvents6=(Int_t)(fAvMult6->GetBinEntries(1));
1592 }
1593
813a4157 1594 Double_t dAvG6[fgkPmax6][fgkQmax6]={{0.}};
cb308e83 1595 Bool_t someAvGEntryIsNegative6=kFALSE;
b6cd16a9 1596 for(Int_t p=0;p<fgkPmax6;p++)
1597 {
1598 for(Int_t q=0;q<fgkQmax6;q++)
1599 {
813a4157 1600 dAvG6[p][q]=fIntGenFun6->GetBinContent(p+1,q+1);
1601 if(dAvG6[p][q]<0.)
cb308e83 1602 {
1603 someAvGEntryIsNegative6=kTRUE;
1604 }
b6cd16a9 1605 }
1606 }
1607 /////////////////////////////////////////////////////////////////////////////
1608 //////////////////gen. function for the cumulants////////////////////////////
1609 /////////////////////////////////////////////////////////////////////////////
1610
813a4157 1611 Double_t dC6[fgkPmax6][fgkQmax6]={{0.}};//C[p][q]
1612 if(dAvM6>0 && someAvGEntryIsNegative6==kFALSE)
cb308e83 1613 {
1614 for (Int_t p=0;p<fgkPmax6;p++)
1615 {
1616 for (Int_t q=0;q<fgkQmax6;q++)
1617 {
813a4157 1618 dC6[p][q]=1.*dAvM6*(pow(dAvG6[p][q],(1./dAvM6))-1.);
b6cd16a9 1619 }
1620 }
1621 }
1622
1623 /////////////////////////////////////////////////////////////////////////////
1624 ///////avaraging the gen. function for the cumulants over azimuth////////////
1625 /////////////////////////////////////////////////////////////////////////////
1626
9dd53ff2 1627 Double_t dAvC6[fgkPmax6]={0.};//<etBinContent(1)C[p][q]>
813a4157 1628 Double_t tempHere6=0.;
b6cd16a9 1629 for (Int_t p=0;p<fgkPmax6;p++){
813a4157 1630 tempHere6=0.;
b6cd16a9 1631 for (Int_t q=0;q<fgkQmax6;q++){
813a4157 1632 tempHere6+=1.*dC6[p][q];
b6cd16a9 1633 }
813a4157 1634 dAvC6[p]=1.*tempHere6/fgkQmax6;
b6cd16a9 1635 }
1636
1637 /////////////////////////////////////////////////////////////////////////////
1638 //////////////////////////////////final results//////////////////////////////
1639 /////////////////////////////////////////////////////////////////////////////
1640
1641 Double_t cumulant6[fgkPmax6];//array to store various order cumulants
813a4157 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]);
b6cd16a9 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
813a4157 1656 Double_t dV2o6=0.,dV4o6=0.,dV6o6=0.;
b6cd16a9 1657
1658 if(cumulant6[0]>=0.)
1659 {
813a4157 1660 dV2o6 = pow(cumulant6[0],(1./2.));
b6cd16a9 1661 }
1662 if(cumulant6[1]<=0.)
1663 {
813a4157 1664 dV4o6 = pow(-cumulant6[1],(1./4.));
b6cd16a9 1665 }
1666 if(cumulant6[2]>=0.)
1667 {
813a4157 1668 dV6o6 = pow((1./4.)*cumulant6[2],(1./6.));
b6cd16a9 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
813a4157 1678 Double_t sdQo6[3]={0.};
1679 Double_t chiQo6[3]={0.};
b6cd16a9 1680
1681 //v_2{2}
813a4157 1682 if(nEvents6 && dAvM6 && cumulant6[0]>=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(cumulant6[0],(1./2.))*dAvM6,2.)>0.))
b6cd16a9 1683 {
813a4157 1684 chiQo6[0]=dAvM6*dV2o6/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2o6*dAvM6,2.),0.5);
1685 if(chiQo6[0])
cb308e83 1686 {
813a4157 1687 sdQo6[0]=pow(((1./(2.*dAvM6*nEvents6))*((1.+2.*pow(chiQo6[0],2))/(2.*pow(chiQo6[0],2)))),0.5);
cb308e83 1688 }
813a4157 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
b6cd16a9 1691 }
1692 else
1693 {
1694 cout<<" v_"<<fgkFlow<<"{2} = Im"<<endl;
1695 }
1696
1697 //v_2{4}
813a4157 1698 if(nEvents6 && dAvM6 && cumulant6[1]<=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(-cumulant6[1],(1./4.))*dAvM6,2.)>0.))
b6cd16a9 1699 {
813a4157 1700 chiQo6[1]=dAvM6*dV4o6/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4o6*dAvM6,2.),0.5);
1701 if(chiQo6[1])
cb308e83 1702 {
813a4157 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);
cb308e83 1704 }
813a4157 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
b6cd16a9 1707 }
1708 else
1709 {
1710 cout<<" v_"<<fgkFlow<<"{4} = Im"<<endl;
1711 }
1712
1713 //v_2{6}
813a4157 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.))
b6cd16a9 1715 {
813a4157 1716 chiQo6[2]=dAvM6*dV6/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV6o6*dAvM6,2.),0.5);
1717 if(chiQo6[2])
cb308e83 1718 {
813a4157 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);
cb308e83 1720 }
813a4157 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
b6cd16a9 1723 }
1724 else
1725 {
1726 cout<<" v_"<<fgkFlow<<"{6} = Im"<<endl;
1727 }
1728
1729 cout<<endl;
813a4157 1730 cout<<" nEvts = "<<nEvents6<<", AvM = "<<dAvM6<<endl;
b6cd16a9 1731 cout<<"***********************************"<<endl;
1732 cout<<"***********************************"<<endl;
1733
1734 //==============================================================================================================================================
1735
1736 //up to 8th order
cb308e83 1737 //avarage selected multiplicity
813a4157 1738 Double_t dAvM8=0.;
cb308e83 1739 if(fAvMult8)
1740 {
813a4157 1741 dAvM8=fAvMult8->GetBinContent(1);
cb308e83 1742 }
1743
1744 //number of events
1745 Int_t nEvents8=0;
1746 if(fAvMult8)
1747 {
1748 nEvents8=(Int_t)(fAvMult8->GetBinEntries(1));
1749 }
1750
813a4157 1751 Double_t dAvG8[fgkPmax8][fgkQmax8]={{0.}};
cb308e83 1752 Bool_t someAvGEntryIsNegative8=kFALSE;
b6cd16a9 1753 for(Int_t p=0;p<fgkPmax8;p++)
1754 {
1755 for(Int_t q=0;q<fgkQmax8;q++)
1756 {
813a4157 1757 dAvG8[p][q]=fIntGenFun8->GetBinContent(p+1,q+1);
1758 if(dAvG8[p][q]<0.)
cb308e83 1759 {
1760 someAvGEntryIsNegative8=kTRUE;
1761 }
b6cd16a9 1762 }
1763 }
1764 /////////////////////////////////////////////////////////////////////////////
1765 //////////////////gen. function for the cumulants////////////////////////////
1766 /////////////////////////////////////////////////////////////////////////////
1767
813a4157 1768 Double_t dC8[fgkPmax8][fgkQmax8]={{0.}};//C[p][q]
1769 if(dAvM8>0 && someAvGEntryIsNegative8==kFALSE)
cb308e83 1770 {
1771 for (Int_t p=0;p<fgkPmax8;p++)
1772 {
1773 for (Int_t q=0;q<fgkQmax8;q++)
1774 {
813a4157 1775 dC8[p][q]=1.*dAvM8*(pow(dAvG8[p][q],(1./dAvM8))-1.);
b6cd16a9 1776 }
1777 }
1778 }
1779
1780 /////////////////////////////////////////////////////////////////////////////
1781 ///////avaraging the gen. function for the cumulants over azimuth////////////
1782 /////////////////////////////////////////////////////////////////////////////
1783
813a4157 1784 Double_t dAvC8[fgkPmax8]={0.};//<C[p][q]>
1785 Double_t tempHere8=0.;
b6cd16a9 1786 for (Int_t p=0;p<fgkPmax8;p++)
1787 {
813a4157 1788 tempHere8=0.;
b6cd16a9 1789 for (Int_t q=0;q<fgkQmax8;q++)
1790 {
813a4157 1791 tempHere8+=1.*dC8[p][q];
b6cd16a9 1792 }
813a4157 1793 dAvC8[p]=1.*tempHere8/fgkQmax8;
b6cd16a9 1794 }
1795
1796 /////////////////////////////////////////////////////////////////////////////
1797 //////////////////////////////////final results//////////////////////////////
1798 /////////////////////////////////////////////////////////////////////////////
1799
1800 Double_t cumulant8[fgkPmax8];//array to store various order cumulants
813a4157 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]);
b6cd16a9 1805
1806 /* x0,y0 ~ 1/sqrt(p) (setting another complex mesh, doesn't work correctly)
813a4157 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]);
b6cd16a9 1811 */
1812
1813 /* x0,y0 ~ p (setting another complex mesh, doesn't work correctly)
813a4157 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]);
b6cd16a9 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
813a4157 1831 Double_t dV2o8=0.,dV4o8=0.,dV6o8=0.,dV8o8=0.;
b6cd16a9 1832
1833 if(cumulant8[0]>=0.)
1834 {
813a4157 1835 dV2o8 = pow(cumulant8[0],(1./2.));
b6cd16a9 1836 }
1837 if(cumulant8[1]<=0.)
1838 {
813a4157 1839 dV4o8 = pow(-cumulant8[1],(1./4.));
b6cd16a9 1840 }
1841 if(cumulant8[2]>=0.)
1842 {
813a4157 1843 dV6o8 = pow((1./4.)*cumulant8[2],(1./6.));
b6cd16a9 1844 }
1845 if(cumulant8[3]<=0.)
1846 {
813a4157 1847 dV8o8 = pow(-(1./33.)*cumulant8[3],(1./8.));
b6cd16a9 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
813a4157 1857 Double_t sdQo8[4]={0.};
1858 Double_t chiQo8[4]={0.};
b6cd16a9 1859
1860 //v_2{2}
813a4157 1861 if(nEvents8 && dAvM8 && cumulant8[0]>=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(cumulant8[0],(1./2.))*dAvM8,2.)>0.))
b6cd16a9 1862 {
813a4157 1863 chiQo8[0]=dAvM8*dV2o8/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2o8*dAvM8,2.),0.5);
1864 if(chiQo8[0])
cb308e83 1865 {
813a4157 1866 sdQo8[0]=pow(((1./(2.*dAvM8*nEvents8))*((1.+2.*pow(chiQo8[0],2.))/(2.*pow(chiQo8[0],2)))),0.5);
cb308e83 1867 }
813a4157 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
b6cd16a9 1870 }
1871 else
1872 {
1873 cout<<" v_"<<fgkFlow<<"{2} = Im"<<endl;
1874 }
1875
1876 //v_2{4}
813a4157 1877 if(nEvents8 && dAvM8 && cumulant8[1]<=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(-cumulant8[1],(1./4.))*dAvM8,2.)>0.))
b6cd16a9 1878 {
813a4157 1879 chiQo8[1]=dAvM8*dV4o8/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4o8*dAvM8,2.),0.5);
1880 if(chiQo8[1])
cb308e83 1881 {
813a4157 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);
cb308e83 1883 }
813a4157 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
b6cd16a9 1886 }
1887 else
1888 {
1889 cout<<" v_"<<fgkFlow<<"{4} = Im"<<endl;
1890 }
1891
1892 //v_2{6}
813a4157 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.))
b6cd16a9 1894 {
813a4157 1895 chiQo8[2]=dAvM8*dV6o8/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV6o8*dAvM8,2.),0.5);
1896 if(chiQo8[2])
cb308e83 1897 {
813a4157 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);
cb308e83 1899 }
813a4157 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
b6cd16a9 1902 }
1903 else
1904 {
1905 cout<<" v_"<<fgkFlow<<"{6} = Im"<<endl;
1906 }
1907
1908 //v_2{8}
813a4157 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.))
b6cd16a9 1910 {
813a4157 1911 chiQo8[3]=dAvM8*dV8o8/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV8o8*dAvM8,2.),0.5);
1912 if(chiQo8[3])
cb308e83 1913 {
813a4157 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);
cb308e83 1915 }
813a4157 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
b6cd16a9 1918 }
1919 else
1920 {
1921 cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl;
1922 }
1923
1924 cout<<endl;
813a4157 1925 cout<<" nEvts = "<<nEvents8<<", AvM = "<<dAvM8<<endl;
b6cd16a9 1926 cout<<"*********************************"<<endl;
1927
1928 //==============================================================================================================================================
1929
1930 //up to 16-th order:
cb308e83 1931 //avarage selected multiplicity
813a4157 1932 Double_t dAvM16=0.;
cb308e83 1933 if(fAvMult16)
1934 {
813a4157 1935 dAvM16=fAvMult16->GetBinContent(1);
cb308e83 1936 }
1937
1938 //number of events
1939 Int_t nEvents16=0;
1940 if(fAvMult16)
1941 {
1942 nEvents16=(Int_t)(fAvMult16->GetBinEntries(1));
1943 }
1944
813a4157 1945 Double_t dAvG16[fgkPmax16][fgkQmax16]={{0.}};
cb308e83 1946 Bool_t someAvGEntryIsNegative16=kFALSE;
b6cd16a9 1947 for(Int_t p=0;p<fgkPmax16;p++)
1948 {
1949 for(Int_t q=0;q<fgkQmax16;q++)
1950 {
813a4157 1951 dAvG16[p][q]=fIntGenFun16->GetBinContent(p+1,q+1);
1952 if(dAvG16[p][q]<0.)
cb308e83 1953 {
1954 someAvGEntryIsNegative16=kTRUE;
1955 }
b6cd16a9 1956 }
1957 }
1958
1959 /////////////////////////////////////////////////////////////////////////////
1960 //////////////////gen. function for the cumulants////////////////////////////
1961 /////////////////////////////////////////////////////////////////////////////
1962
813a4157 1963 Double_t dC16[fgkPmax16][fgkQmax16]={{0.}};//C16[p][q]
1964 if(dAvM16>0 && someAvGEntryIsNegative16==kFALSE)
b6cd16a9 1965 {
1966 for(Int_t p=0;p<fgkPmax16;p++)
1967 {
1968 for(Int_t q=0;q<fgkQmax16;q++)
1969 {
813a4157 1970 dC16[p][q]=1.*dAvM16*(pow(dAvG16[p][q],(1./dAvM16))-1.);
b6cd16a9 1971 }
1972 }
1973 }
1974
1975 /////////////////////////////////////////////////////////////////////////////
1976 ///////avaraging the gen. function for the cumulants over azimuth////////////
1977 /////////////////////////////////////////////////////////////////////////////
1978
813a4157 1979 Double_t dAvC16[fgkPmax16]={0.};//<C16[p][q]>
1980 Double_t tempHere16=0.;
b6cd16a9 1981 for (Int_t p=0;p<fgkPmax16;p++)
1982 {
813a4157 1983 tempHere16=0.;
b6cd16a9 1984 for (Int_t q=0;q<fgkQmax16;q++)
1985 {
813a4157 1986 tempHere16+=1.*dC16[p][q];
b6cd16a9 1987 }
813a4157 1988 dAvC16[p]=1.*tempHere16/fgkQmax16;
b6cd16a9 1989 }
1990
1991 /////////////////////////////////////////////////////////////////////////////
1992 //////////////////////////////////final results//////////////////////////////
1993 /////////////////////////////////////////////////////////////////////////////
1994
1995 Double_t cumulant16[fgkPmax16];//array to store various order cumulants
1996
813a4157 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]);
b6cd16a9 1999
813a4157 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]);
b6cd16a9 2003
813a4157 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]);
b6cd16a9 2007
813a4157 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]);
b6cd16a9 2011
813a4157 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]);
b6cd16a9 2014
813a4157 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]);
b6cd16a9 2017
813a4157 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]);
b6cd16a9 2020
813a4157 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]);
b6cd16a9 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
813a4157 2039 Double_t dV2o16=0.,dV4o16=0.,dV6o16=0.,dV8o16=0.,V10o16=0.,V12o16=0.,V14o16=0.,V16o16=0.;
b6cd16a9 2040
2041 if(cumulant16[0]>=0.)
2042 {
813a4157 2043 dV2o16 = pow(cumulant16[0],(1./2.));
b6cd16a9 2044 }
2045 if(cumulant16[1]<=0.)
2046 {
813a4157 2047 dV4o16 = pow(-cumulant16[1],(1./4.));
b6cd16a9 2048 }
2049 if(cumulant16[2]>=0.)
2050 {
813a4157 2051 dV6o16 = pow((1./4.)*cumulant16[2],(1./6.));
b6cd16a9 2052 }
2053 if(cumulant16[3]<=0.)
2054 {
813a4157 2055 dV8o16 = pow(-(1./33.)*cumulant16[3],(1./8.));
b6cd16a9 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
813a4157 2081 Double_t sdQo16[8]={0.};
2082 Double_t chiQo16[8]={0.};
b6cd16a9 2083
2084 //v_2{2}
813a4157 2085 if(nEvents16 && dAvM16 && cumulant16[0]>=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(cumulant16[0],(1./2.))*dAvM16,2.)>0.))
b6cd16a9 2086 {
813a4157 2087 chiQo16[0]=dAvM16*dV2o16/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV2o16*dAvM16,2.),0.5);
2088 if(chiQo16[0])
cb308e83 2089 {
813a4157 2090 sdQo16[0]=pow(((1./(2.*dAvM16*nEvents16))*((1.+2.*pow(chiQo16[0],2))/(2.*pow(chiQo16[0],2)))),0.5);
cb308e83 2091 }
813a4157 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
b6cd16a9 2094 }
2095 else
2096 {
2097 cout<<" v_"<<fgkFlow<<"{2} = Im"<<endl;
2098 }
2099
2100 //v_2{4}
813a4157 2101 if(nEvents16 && dAvM16 && cumulant16[1]<=0. && (dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(pow(-cumulant16[1],(1./4.))*dAvM16,2.)>0.))
b6cd16a9 2102 {
813a4157 2103 chiQo16[1]=dAvM16*dV4o16/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV4o16*dAvM16,2.),0.5);
2104 if(chiQo16[1])
cb308e83 2105 {
813a4157 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);
cb308e83 2107 }
813a4157 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
b6cd16a9 2110 }
2111 else
2112 {
2113 cout<<" v_"<<fgkFlow<<"{4} = Im"<<endl;
2114 }
2115
2116 //v_2{6}
813a4157 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.))
b6cd16a9 2118 {
813a4157 2119 chiQo16[2]=dAvM16*dV6o16/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV6o16*dAvM16,2.),0.5);
2120 if(chiQo16[2])
cb308e83 2121 {
813a4157 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);
cb308e83 2123 }
813a4157 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
b6cd16a9 2126 }
2127 else
2128 {
2129 cout<<" v_"<<fgkFlow<<"{6} = Im"<<endl;
2130 }
2131
2132 //v_2{8}
813a4157 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.))
b6cd16a9 2134 {
813a4157 2135 chiQo16[3]=dAvM16*dV8o16/pow(dAvQ2x+dAvQ2y-pow(dAvQx,2.)-pow(dAvQy,2.)-pow(dV8o16*dAvM16,2.),0.5);
2136 if(chiQo16[3])
cb308e83 2137 {
813a4157 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);
cb308e83 2139 }
813a4157 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
b6cd16a9 2142 }
2143 else
2144 {
2145 cout<<" v_"<<fgkFlow<<"{8} = Im"<<endl;
2146 }
2147
2148 //v_2{10}
813a4157 2149 if(nEvents16 && dAvM16 && cumulant16[4]>=0.)
b6cd16a9 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}
813a4157 2159 if(nEvents16 && dAvM16 && cumulant16[5]<=0.)
b6cd16a9 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}
813a4157 2169 if(nEvents16 && dAvM16 && cumulant16[6]>=0.)
b6cd16a9 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}
813a4157 2179 if(nEvents16 && dAvM16 && cumulant16[7]<=0.)
b6cd16a9 2180 {
2181 cout<<"v_"<<fgkFlow<<"{16} = "<<V16o16<<endl;
2182 }
2183 else
2184 {
2185 cout<<"v_"<<fgkFlow<<"{16} = Im"<<endl;
2186 }
2187
2188 cout<<endl;
813a4157 2189 cout<<" nEvts = "<<nEvents16<<", AvM = "<<dAvM16<<endl;
b6cd16a9 2190 cout<<"*********************************"<<endl;
2191
2192 //==============================================================================================================================================
2193
2194 }//end of if(otherEquations)
aaebd73d 2195}
2196
2197
2198