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