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