]>
Commit | Line | Data |
---|---|---|
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 | ||
16 | class TH1; | |
17 | class TGraph; | |
18 | class TPave; | |
19 | class TLatex; | |
20 | class TMarker; | |
21 | class TRandom3; | |
22 | class TObjArray; | |
23 | class TList; | |
24 | class TCanvas; | |
25 | class TSystem; | |
26 | class TROOT; | |
27 | ||
28 | class AliFlowVector; | |
29 | class TVector; | |
30 | ||
31 | //******************************* | |
32 | // flow analysis with cumulants * | |
33 | // author: Ante Bilandzic * | |
34 | // email: anteb@nikhef.nl * | |
35 | //******************************* | |
36 | ||
37 | ClassImp(AliFlowAnalysisWithCumulants) | |
38 | ||
39 | //________________________________________________________________________ | |
40 | ||
41 | AliFlowAnalysisWithCumulants::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 | ||
83 | AliFlowAnalysisWithCumulants::~AliFlowAnalysisWithCumulants(){ | |
84 | //desctructor | |
85 | } | |
86 | ||
87 | //___________________________________________________________________________ | |
88 | void 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 | //________________________________________________________________________ | |
100 | void 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 | //________________________________________________________________________ | |
210 | void 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 |