]>
Commit | Line | Data |
---|---|---|
a663ab5a | 1 | // Macro plotCumulants.C is used to represent true correlations, a.k.a. cumulants, |
2 | // as a function of reference multiplicity. Different order cumulants provide | |
3 | // independent estimates of the flow harmonics. In the ideal case when only flow | |
4 | // correlations are present we have: | |
5 | // c{2} = v^2, c{4} = -v^4, c{6} = 4v^6 and c{8} = -33v^8. | |
6 | // In practice cumulants can are calculated via Q-vectors (QC) or via formalism | |
7 | // of generating functions (GFC). | |
8 | ||
9dcf9913 | 9 | // Set how many output analysis files in total you want to access: |
34f4fef7 | 10 | const Int_t nFiles = 2; |
9dcf9913 | 11 | |
12 | // Set how many of those output analysis files you want to represent with a mesh (usually used to represent results of simulations): | |
a9a2547e | 13 | const Int_t nSim = 0; |
dfb33a36 | 14 | |
9dcf9913 | 15 | // Set paths of all output analysis files (first the ones to be represented with mesh (simulations), then the ones to be represented with markers (real data)) |
34f4fef7 | 16 | TString files[nFiles] = {"trackletsCorrectedNUA","default"}; |
63cdf10a | 17 | |
9dcf9913 | 18 | // Set analysis types for all output analysis files (can be "ESD","AOD","MC",""): |
34f4fef7 | 19 | TString type[nFiles] = {"ESD","ESD"}; |
9dcf9913 | 20 | |
21 | // Set mesh color: | |
a9a2547e | 22 | Int_t meshColor[nSim] = {}; |
9dcf9913 | 23 | |
24 | // Set marker styles: | |
34f4fef7 | 25 | Int_t markerStyle[nFiles-nSim] = {kFullSquare,kOpenSquare}; |
9dcf9913 | 26 | |
c109b8c7 | 27 | // Set marker colors: |
34f4fef7 | 28 | Int_t markerColor[nFiles-nSim] = {kBlack,kRed}; |
c109b8c7 | 29 | |
9dcf9913 | 30 | // Set legend entries: |
34f4fef7 | 31 | TString legendEntry[nFiles] = {"",""}; |
63cdf10a | 32 | |
83cac266 | 33 | // Set if you want to rebin the histograms into wider multiplicity bins (set for each cumulant order separately): |
34f4fef7 | 34 | Bool_t rebin = kFALSE; |
35 | Int_t nMergedBins[4] = {10,10,10,10}; // set how many original multiplicity bins will be merged into 1 new one | |
9dcf9913 | 36 | |
c109b8c7 | 37 | // Set if you whish to plot cumulants versus <reference multiplicity> (by default they are plotted versus # of RPs): |
34f4fef7 | 38 | Bool_t plotCumulantsVsReferenceMultiplicity = kFALSE; |
a9a2547e | 39 | Bool_t showReferenceMultiplicityVsNoOfRPs = kFALSE; |
c109b8c7 | 40 | |
9dcf9913 | 41 | // Set flow values whose theoretical contribution to cumulants will be shown on the plots with the straight coloured lines: |
a9a2547e | 42 | Bool_t showTheoreticalLines = kFALSE; |
63cdf10a | 43 | const Int_t nFlowValues = 1; |
a663ab5a | 44 | Double_t v[nFlowValues] = {0.05}; |
63cdf10a | 45 | Int_t lineColor[nFlowValues] = {kRed}; |
9dcf9913 | 46 | |
47 | // If the statistical error of 6th and 8th order cumulant is huge you may prefer not to show them: | |
16eed38e | 48 | Bool_t plotOnly2ndAnd4thOrderCumulant = kFALSE; |
49 | ||
50 | // Set if you want independent canvas with results for reference flow vs multiplicity: | |
51 | Bool_t showRefFlowVsM = kTRUE; | |
34f4fef7 | 52 | Int_t refFlowVsMMarkerStyle[nFiles] = {kFullSquare,kOpenSquare}; // marker style is different for different file |
53 | Int_t refFlowVsMMarkerColor[4] = {kBlack,kRed,kBlue,kGreen+2}; // marker color is different for different cumulant order | |
54 | Int_t refFlowVsMeshOrder = -1; // set results of which order will be plotted as mesh in all pads [0=2nd, 1=4th, 2=6th, 3=8th, -1=do not show] | |
55 | Int_t refFlowVsMMeshColor = kRed-10; // mesh color for above specified order | |
16eed38e | 56 | Double_t refFlowVsMxRange[2] = {1.,11000.}; // x range on the plots for reference multiplicity vs M |
57 | Double_t refFlowVsMyRange[2] = {0.0,0.194}; // x range on the plots for reference multiplicity vs M | |
9dcf9913 | 58 | |
a9a2547e | 59 | // Set if you want to show theoretical curves for the toy model: |
60 | Bool_t showToyModel = kFALSE; | |
61 | const Int_t nToyModels = 2; // number of toy models with different parameters k, vn, v2n | |
62 | Double_t k[nToyModels] = {2.,4.}; | |
63 | Double_t vn[nToyModels] = {0.25,0.25}; | |
64 | Double_t v2n[nToyModels] = {0.0,0.0}; | |
65 | Int_t lineColorToyModel[nToyModels] = {kBlack,kRed}; | |
66 | ||
9dcf9913 | 67 | // For comparison sake show also GFC results with dotted line: |
68 | Bool_t showAlsoGFCResults = kFALSE; | |
69 | Int_t gfcLineStyle = 3; | |
70 | ||
dfb33a36 | 71 | // For comparison sake show also Monte Carlo QC results with coloured mesh: |
72 | Bool_t showAlsoMCEPResults = kFALSE; | |
73 | Bool_t showOnlyMCEPResults = kFALSE; | |
a9a2547e | 74 | Int_t mcepMeshColor[nSim] = {}; |
dfb33a36 | 75 | |
76 | // Set method names which calculate cumulants vs multiplicity (do not touch these settings unless you are looking for a trouble): | |
77 | const Int_t nMethods = 3; | |
78 | TString method[nMethods] = {"QC","GFC","MCEP"}; | |
9dcf9913 | 79 | |
80 | TFile *commonOutputFiles[nFiles] = {NULL}; // common output files "AnalysisResults.root" | |
81 | TList *lists[nFiles][nMethods] = {{NULL}}; // lists cobj<method> holding objects with results for each method | |
82 | TH1D *cumulantsVsM[nFiles][nMethods][4] = {{{NULL}}}; // histograms with results for cumulants vs multiplicity (4 stands for 4 cumulant orders) | |
16eed38e | 83 | TH1D *refFlowVsM[nFiles][nMethods][4] = {{{NULL}}}; // histograms with results for reference flow vs multiplicity (4 stands for 4 cumulant orders) |
63cdf10a | 84 | TGraph *lines[nFlowValues][4] = {{NULL}}; // lines denoting theoretical flow contribution to cumulants |
85 | TProfile *refMultVsNoOfRPs[nFiles] = {NULL}; // <reference multipicity> versus # of RPs | |
a9a2547e | 86 | TF1 *toyModel[2][nToyModels] = {{NULL}}; // [cumulant order][number of toy models] |
9dcf9913 | 87 | |
88 | // Ranges for plots: | |
89 | Double_t xMin[4]={0.}; | |
90 | Double_t xMax[4]={0.}; | |
91 | Double_t yMin[4]={0.}; | |
92 | Double_t yMax[4]={0.}; | |
93 | ||
94 | enum libModes {mLocal,mLocalSource}; | |
95 | ||
96 | void plotCumulants(Int_t analysisMode=mLocal) | |
97 | { | |
98 | // analysisMode: if analysisMode = mLocal -> analyze data on your computer using aliroot | |
99 | // if analysisMode = mLocalSource -> analyze data on your computer using root + source files | |
a663ab5a | 100 | |
101 | // Cross-check user settings: | |
102 | CrossCheckSettings(); | |
63cdf10a | 103 | |
9dcf9913 | 104 | // Load needed libraries: |
105 | LoadLibrariesPC(analysisMode); | |
106 | ||
dfb33a36 | 107 | // Global settings which will affect all plots: |
108 | GlobalSettings(); | |
109 | ||
9dcf9913 | 110 | // Access all common output files: |
111 | TString commonOutputFileName = "AnalysisResults.root"; | |
112 | AccessCommonOutputFiles(commonOutputFileName); | |
113 | ||
114 | // Get from common output files the lists holding histograms for each method: | |
115 | GetLists(); | |
116 | ||
117 | // Get histograms with results for cumulants vs multiplicity: | |
118 | GetHistograms(); | |
16eed38e | 119 | |
9dcf9913 | 120 | // Determine ranges for plots: |
121 | DetermineMinMax(); | |
dfb33a36 | 122 | |
9dcf9913 | 123 | // Make lines which indicate theoretical contributions of flow to cumulants: |
124 | Lines(); | |
125 | ||
126 | // Print number of events and average multiplicities for each common output file: | |
127 | Print(); | |
128 | ||
16eed38e | 129 | // Make plots for cumulants vs multiplicity: // to be improved |
130 | PlotCumulantsVsM(); | |
131 | ||
132 | // Make plots for reference flow vs multiplicity: | |
133 | if(showRefFlowVsM){PlotRefFlowVsM();} | |
a663ab5a | 134 | |
34f4fef7 | 135 | } // end of void plotCumulants(Int_t analysisMode=mLocal) |
136 | ||
9dcf9913 | 137 | // ===================================================================================== |
138 | ||
16eed38e | 139 | void PlotRefFlowVsM() |
9dcf9913 | 140 | { |
34f4fef7 | 141 | // Make plots for reference flow vs multiplicity. |
142 | ||
143 | // Calculate reference flow from cumulants vs M: | |
144 | CalculateReferenceFlowFromCumulantsVsM(); | |
145 | ||
16eed38e | 146 | TCanvas *cRefFlowVsM = new TCanvas("cRefFlowVsM","Reference Flow"); |
147 | cRefFlowVsM->Divide(2,2); | |
148 | ||
149 | TLegend *lRefFlowVsM = new TLegend(0.1,0.7,0.33,0.9); | |
150 | lRefFlowVsM->SetFillStyle(0); | |
34f4fef7 | 151 | lRefFlowVsM->SetHeader(" Therminator (20-30%)"); |
16eed38e | 152 | |
153 | TString refFlowVsMFlag[4] = {"v_{2}{2,QC}","v_{2}{4,QC}","v_{2}{6,QC}","v_{2}{8,QC}"}; | |
154 | ||
34f4fef7 | 155 | for(Int_t f=0;f<nFiles;f++) |
16eed38e | 156 | { |
34f4fef7 | 157 | for(Int_t m=0;m<1+(Int_t)(showAlsoGFCResults)+(Int_t)(showAlsoMCEPResults);m++) |
16eed38e | 158 | { |
34f4fef7 | 159 | for(Int_t co=0;co<4;co++) // cumulant order |
160 | { | |
161 | cRefFlowVsM->cd(co+1); | |
162 | // Style histogram: | |
163 | if(f==0) // superimpose histograms from other files on top of this one | |
164 | { | |
165 | TH1D *styleHist = (TH1D*)StyleHist(refFlowVsMFlag[co].Data(),co)->Clone(); | |
166 | if(styleHist) | |
167 | { | |
168 | styleHist->GetXaxis()->SetRangeUser(refFlowVsMxRange[0],refFlowVsMxRange[1]); | |
169 | styleHist->GetYaxis()->SetRangeUser(refFlowVsMyRange[0],refFlowVsMyRange[1]); | |
170 | styleHist->Draw(); | |
171 | } | |
172 | } // end of if(f==0) | |
173 | // Plot first the mesh for reference flow of specified order in all pads: | |
174 | if(refFlowVsMeshOrder != -1) | |
175 | { | |
176 | TGraph *rfMeshOutOf = GetErrorMesh(refFlowVsM[f][m][refFlowVsMeshOrder]); | |
177 | if(rfMeshOutOf) | |
178 | { | |
179 | rfMeshOutOf->SetFillColor(refFlowVsMMeshColor); | |
180 | rfMeshOutOf->Draw("lfsame"); | |
181 | if(refFlowVsMeshOrder==co){lRefFlowVsM->AddEntry(rfMeshOutOf,Form("v_{2}{%i,QC} stat. error",2*(refFlowVsMeshOrder+1)),"f");} | |
182 | } | |
183 | } // end of if(refFlowVsMeshOrder != -1) | |
184 | // Plot reference flow vs M: | |
185 | refFlowVsM[f][m][co]->SetMarkerStyle(refFlowVsMMarkerStyle[f]); | |
186 | refFlowVsM[f][m][co]->SetMarkerColor(refFlowVsMMarkerColor[co]); | |
187 | refFlowVsM[f][m][co]->Draw("e1same"); | |
188 | lRefFlowVsM->AddEntry(refFlowVsM[f][m][co],Form("v_{2}{%i,QC}, %s",2*(co+1),files[f].Data()),"p"); | |
189 | } // end of for(Int_t co=0;co<4;co++) // cumulant order | |
190 | } // end of for(Int_t m=0;m<nMethods;m++) | |
191 | } // end of for(Int_t f=0;f<nFiles;f++) | |
16eed38e | 192 | |
34f4fef7 | 193 | // Draw a common legend in the 1st pad: |
16eed38e | 194 | cRefFlowVsM->cd(1); |
195 | lRefFlowVsM->Draw("same"); | |
196 | ||
197 | } // end of void PlotRefFlowVsM() | |
198 | ||
199 | // ===================================================================================== | |
200 | ||
201 | void PlotCumulantsVsM() | |
202 | { | |
203 | // Make plots for cumulants vs multiplicity: // to be improved | |
9dcf9913 | 204 | |
205 | TCanvas *c = NULL; | |
206 | Int_t coMax = 0; | |
207 | if(!plotOnly2ndAnd4thOrderCumulant) | |
208 | { | |
209 | c = new TCanvas("c","cumulants"); | |
210 | c->Divide(2,2); | |
211 | coMax = 4; | |
212 | } else | |
213 | { | |
214 | c = new TCanvas("c","cumulants",1200,500); | |
215 | c->Divide(2,1); | |
216 | coMax = 2; | |
a9a2547e | 217 | } |
9dcf9913 | 218 | |
219 | TLegend *legend = new TLegend(0.1,0.7,0.33,0.9); | |
220 | legend->SetFillStyle(0); | |
63cdf10a | 221 | //legend->SetHeader(" minClustersTpcRP"); |
9dcf9913 | 222 | |
223 | TString qcFlag[4] = {"QC{2}","QC{4}","QC{6}","QC{8}"}; | |
224 | ||
225 | for(Int_t co=0;co<coMax;co++) // cumulant order | |
226 | { | |
227 | c->cd(co+1); | |
16eed38e | 228 | TH1D *styleHist = (TH1D*) StyleHist(qcFlag[co].Data(),co)->Clone(); |
229 | if(styleHist){styleHist->Draw();} | |
a9a2547e | 230 | if(co==0) |
231 | { | |
232 | if(showToyModel) | |
233 | { | |
234 | for(Int_t ntm=0;ntm<nToyModels;ntm++) | |
235 | { | |
236 | TF1 *tm = ToyModel(co,k[ntm],vn[ntm],v2n[ntm]); | |
237 | if(tm) | |
238 | { | |
239 | tm->SetLineColor(lineColorToyModel[ntm]); | |
240 | tm->Draw("same"); | |
241 | } | |
242 | } // for(Int_t ntm=0;ntm<nToyModels;ntm++) | |
243 | } // end of if(showToyModel) | |
244 | } // end of if(co==0) | |
245 | if(co==1) | |
246 | { | |
247 | if(showToyModel) | |
248 | { | |
249 | for(Int_t ntm=0;ntm<nToyModels;ntm++) | |
250 | { | |
251 | TF1 *tm = ToyModel(co,k[ntm],vn[ntm],v2n[ntm]); | |
252 | if(tm) | |
253 | { | |
254 | tm->SetLineColor(lineColorToyModel[ntm]); | |
255 | tm->Draw("same"); | |
256 | } | |
257 | } // for(Int_t ntm=0;ntm<nToyModels;ntm++) | |
258 | } // end of if(showToyModel) | |
259 | } // end of if(co==1) | |
260 | ||
9dcf9913 | 261 | // simulations: |
262 | for(Int_t s=0;s<nSim;s++) | |
263 | { | |
dfb33a36 | 264 | // Monte Carlo QC: |
265 | if(showAlsoMCEPResults) | |
266 | { | |
267 | TGraph *mcepQC = GetErrorMesh(cumulantsVsM[s][2][co]); | |
268 | if(mcepQC) | |
269 | { | |
270 | mcepQC->SetFillColor(mcepMeshColor[s]); | |
271 | mcepQC->Draw("lfsame"); | |
272 | } | |
273 | if(co==0){legend->AddEntry(mcepQC,Form("%s (MC)",legendEntry[s].Data()),"f");} | |
274 | } // end of if(showAlsoMCEPResults) | |
275 | // QC: | |
9dcf9913 | 276 | TGraph *errorMesh = GetErrorMesh(cumulantsVsM[s][0][co]); |
277 | if(errorMesh) | |
278 | { | |
279 | errorMesh->SetFillColor(meshColor[s]); | |
dfb33a36 | 280 | if(!(showOnlyMCEPResults && showAlsoMCEPResults)){errorMesh->Draw("lfsame");} |
9dcf9913 | 281 | } |
dfb33a36 | 282 | if(co==0 && !(showOnlyMCEPResults && showAlsoMCEPResults)){legend->AddEntry(errorMesh,legendEntry[s].Data(),"f");} |
9dcf9913 | 283 | } // end of if(Int_t s=0;s<nSim;s++) |
dfb33a36 | 284 | // Theoretical lines: |
285 | if(showTheoreticalLines) | |
286 | { | |
287 | for(Int_t fv=0;fv<nFlowValues;fv++) | |
288 | { | |
289 | lines[fv][co]->Draw("lsame"); | |
290 | if(co==0){legend->AddEntry(lines[fv][co],Form("v_{2} = %g",v[fv]),"l");} | |
291 | } | |
292 | } | |
9dcf9913 | 293 | // data: |
294 | for(Int_t f=nSim;f<nFiles;f++) | |
295 | { | |
9dcf9913 | 296 | // QC results: |
297 | if(cumulantsVsM[f][0][co]) | |
298 | { | |
299 | cumulantsVsM[f][0][co]->Draw("e1same"); | |
c109b8c7 | 300 | cumulantsVsM[f][0][co]->SetMarkerStyle(markerStyle[f-nSim]); |
301 | cumulantsVsM[f][0][co]->SetMarkerColor(markerColor[f-nSim]); | |
34f4fef7 | 302 | if(co==0) |
9dcf9913 | 303 | { |
304 | if(showAlsoGFCResults) | |
305 | { | |
306 | legend->AddEntry(cumulantsVsM[f][0][co],Form("%s (QC)",legendEntry[f].Data()),"p"); | |
307 | } else | |
308 | { | |
309 | legend->AddEntry(cumulantsVsM[f][0][co],legendEntry[f].Data(),"p"); | |
310 | } | |
311 | } | |
312 | } | |
313 | // GFC results: | |
314 | if(showAlsoGFCResults && cumulantsVsM[f][1][co]) | |
315 | { | |
316 | cumulantsVsM[f][1][co]->Draw("lsame"); | |
317 | cumulantsVsM[f][1][co]->SetLineStyle(gfcLineStyle); | |
318 | if(co==0){legend->AddEntry(cumulantsVsM[f][1][co],Form("%s (GFC)",legendEntry[f].Data()),"l");} | |
319 | } | |
320 | } | |
321 | // Draw legend: | |
16eed38e | 322 | if(co==0){legend->Draw("same");} |
9dcf9913 | 323 | } // end of for(Int_t co=0;co<4;co++) // cumulant order |
a663ab5a | 324 | |
c109b8c7 | 325 | // Plot also <reference multiplicity> vs # of RPs: |
a663ab5a | 326 | if(plotCumulantsVsReferenceMultiplicity && showReferenceMultiplicityVsNoOfRPs) |
c109b8c7 | 327 | { |
63cdf10a | 328 | TCanvas *cRefMultVsNoOfRPs = new TCanvas("cRefMultVsNoOfRPs","#LTreference multiplicity#GT vs # of RPs",1200,600); |
329 | cRefMultVsNoOfRPs->Divide(nFiles,1); | |
330 | for(Int_t f=0;f<nFiles;f++) | |
c109b8c7 | 331 | { |
63cdf10a | 332 | cRefMultVsNoOfRPs->cd(f+1); |
333 | if(refMultVsNoOfRPs[f]) | |
334 | { | |
335 | refMultVsNoOfRPs[f]->SetTitle(legendEntry[f].Data()); | |
336 | refMultVsNoOfRPs[f]->GetXaxis()->SetRangeUser(0,refMultVsNoOfRPs[f]->FindLastBinAbove()); | |
337 | refMultVsNoOfRPs[f]->Draw(); | |
338 | } else | |
339 | { | |
340 | cout<<endl; | |
341 | cout<<"WARNING: refMultVsNoOfRPs[f] is NULL in Plot(), f = "<<f<<" !!!!"<<endl; | |
342 | cout<<endl; | |
343 | } | |
344 | } | |
c109b8c7 | 345 | } |
346 | ||
16eed38e | 347 | } // end of void PlotCumulantsVsM() |
9dcf9913 | 348 | |
349 | // ===================================================================================== | |
350 | ||
a663ab5a | 351 | void CrossCheckSettings() |
352 | { | |
353 | // Cross-check user settings in this method. | |
354 | ||
355 | if(showAlsoGFCResults && rebin) | |
356 | { | |
357 | cout<<endl; | |
358 | cout<<" WARNING: Rebinning in M not supported for GFC yet !!!!"<<endl; | |
359 | cout<<endl; | |
360 | exit(0); | |
361 | } | |
362 | ||
363 | if(showAlsoGFCResults && plotCumulantsVsReferenceMultiplicity) | |
364 | { | |
365 | cout<<endl; | |
366 | cout<<" WARNING: Showing GFC versus <reference multiplicity> not supported yet !!!!"<<endl; | |
367 | cout<<endl; | |
368 | exit(0); | |
369 | } | |
370 | ||
371 | } // end of void CrossCheckSettings() | |
372 | ||
373 | // ===================================================================================== | |
374 | ||
9dcf9913 | 375 | void Lines() |
376 | { | |
377 | // Make lines denoting theoretical contribution of flow to cumulants. | |
9dcf9913 | 378 | |
379 | for(Int_t fv=0;fv<nFlowValues;fv++) | |
380 | { | |
63cdf10a | 381 | lines[fv][0] = new TGraph(2); |
382 | lines[fv][0]->SetPoint(0,xMin[0],pow(v[fv],2)); | |
383 | lines[fv][0]->SetPoint(1,xMax[0]+0.5,pow(v[fv],2)); | |
9dcf9913 | 384 | lines[fv][0]->SetLineColor(lineColor[fv]); |
63cdf10a | 385 | lines[fv][1] = new TGraph(2); |
386 | lines[fv][1]->SetPoint(0,xMin[1],-pow(v[fv],4)); | |
387 | lines[fv][1]->SetPoint(1,xMax[1]+0.5,-pow(v[fv],4)); | |
388 | lines[fv][1]->SetLineColor(lineColor[fv]); | |
389 | lines[fv][2] = new TGraph(2); | |
390 | lines[fv][2]->SetPoint(0,xMin[2],4.*pow(v[fv],6)); | |
391 | lines[fv][2]->SetPoint(1,xMax[2]+0.5,4.*pow(v[fv],6)); | |
392 | lines[fv][2]->SetLineColor(lineColor[fv]); | |
393 | lines[fv][3] = new TGraph(2); | |
394 | lines[fv][3]->SetPoint(0,xMin[3],-33.*pow(v[fv],8)); | |
395 | lines[fv][3]->SetPoint(1,xMax[3]+0.5,-33.*pow(v[fv],8)); | |
9dcf9913 | 396 | lines[fv][3]->SetLineColor(lineColor[fv]); |
397 | } | |
398 | ||
399 | } // end of void Lines() | |
400 | ||
401 | // ===================================================================================== | |
402 | ||
a9a2547e | 403 | TF1* ToyModel(Int_t co, Double_t k, Double_t vn, Double_t v2n) |
404 | { | |
405 | // Make theoretical curves for the toy model. | |
406 | ||
407 | TF1 *function = NULL; | |
408 | ||
409 | if(co==0) | |
410 | { | |
411 | function = new TF1("function","([1]*[1]*(x-0.5-[0])+[0]-1)/((x-0.5)-1)",1.5,1000); | |
412 | function->SetParameter(0,k); | |
413 | function->SetParameter(1,vn); | |
414 | } | |
415 | else if(co==1) | |
416 | { | |
417 | function = new TF1("function","-1.*([1]*[1]*[1]*[1]*([0]-(x-0.5))*(12.*[0]-6.*[0]*[0]-12.*(x-0.5)-5.*[0]*(x-0.5)+6.*[0]*[0]*(x-0.5)+9.*(x-0.5)*(x-0.5)-3.*[0]*(x-0.5)*(x-0.5)-(x-0.5)*(x-0.5)*(x-0.5))+[1]*[1]*4.*([0]-1.)*(x-0.5-[0])*([0]*(x-0.5-1.)-2.*(x-0.5-2.))-[1]*[1]*[2]*2.*(x-0.5-1.)*([0]-1.)*(x-0.5-[0])*(x-0.5-2.*[0])-[2]*[2]*([0]-1.)*([0]-1.)*((x-0.5)-[0])*(x-0.5-1.)+([0]-1.)*(-6.+9*[0]-[0]*[0]+2.*(x-0.5)-5.*[0]*(x-0.5)+[0]*[0]*(x-0.5)))/((x-0.5-1.)*(x-0.5-1.)*(x-0.5-2.)*(x-0.5-3.))",3.5,100); | |
418 | function->SetParameter(0,k); | |
419 | function->SetParameter(1,vn); | |
420 | function->SetParameter(2,v2n); | |
421 | } | |
422 | ||
423 | return function; | |
424 | ||
425 | } // end of TF1* ToyModel(Int_t co) | |
426 | ||
427 | // ===================================================================================== | |
428 | ||
34f4fef7 | 429 | void CalculateReferenceFlowFromCumulantsVsM() |
430 | { | |
431 | // Calculate reference flow from cumulants vs M: | |
432 | ||
433 | for(Int_t f=0;f<nFiles;f++) | |
434 | { | |
435 | for(Int_t m=0;m<1+(Int_t)(showAlsoGFCResults)+(Int_t)(showAlsoMCEPResults);m++) | |
436 | { | |
437 | for(Int_t co=0;co<4;co++) // cumulant order | |
438 | { | |
439 | if(cumulantsVsM[f][m][co]) | |
440 | { | |
441 | refFlowVsM[f][m][co] = (TH1D*)cumulantsVsM[f][m][co]->Clone(Form("%i,%i,%i",f,m,co)); | |
442 | if(!refFlowVsM[f][m][co]){cout<<" WARNING: "<<Form("refFlowVsM[%i][%i][%i]",f,m,co)<<" is NULL in PlotRefFlowVsM() !!!!"<<exit(0)<<endl;} | |
443 | refFlowVsM[f][m][co]->Reset(); // to have empty histogram with the same binning as the one with cumulants | |
444 | Int_t nBins = refFlowVsM[f][m][co]->GetNbinsX(); | |
445 | for(Int_t b=1;b<=nBins;b++) | |
446 | { | |
447 | Double_t qcVsM = cumulantsVsM[f][m][co]->GetBinContent(b); // QC vs M | |
448 | Double_t qcErrorVsM = cumulantsVsM[f][m][co]->GetBinError(b); // QC vs M stat. error | |
449 | Double_t vVsM = 0.; // reference flow vs M | |
450 | Double_t vErrorVsM = 0.; // reference flow vs M stat. error | |
451 | if(co==0) // 2nd order | |
452 | { | |
453 | if(qcVsM>0.) | |
454 | { | |
455 | vVsM = pow(qcVsM,1./2.); | |
456 | vErrorVsM = (1./2.)*pow(qcVsM,-1./2.)*qcErrorVsM; | |
457 | } | |
458 | } // end of if(co==0) 2nd order | |
459 | else if(co==1) // 4th order | |
460 | { | |
461 | if(qcVsM<0.) | |
462 | { | |
463 | vVsM = pow(-1.*qcVsM,1./4.); | |
464 | vErrorVsM = (1./4.)*pow(-qcVsM,-3./4.)*qcErrorVsM; | |
465 | } | |
466 | } // end of if(co==1) 4th order | |
467 | else if(co==2) // 6th order | |
468 | { | |
469 | if(qcVsM>0.) | |
470 | { | |
471 | vVsM = pow((1./4.)*qcVsM,1./6.); | |
472 | vErrorVsM = (1./6.)*pow(2.,-1./3.)*pow(qcVsM,-5./6.)*qcErrorVsM; | |
473 | } | |
474 | } // end of if(co==2) 6th order | |
475 | else if(co==3) // 8th order | |
476 | { | |
477 | if(qcVsM<0.) | |
478 | { | |
479 | vVsM = pow((-1./33.)*qcVsM,1./8.); | |
480 | vErrorVsM = (1./8.)*pow(33.,-1./8.)*pow(-qcVsM,-7./8.)*qcErrorVsM; | |
481 | } | |
482 | } // end of if(co==3) 8th order | |
483 | // Store final results and statisticial errror for reference flow: | |
484 | refFlowVsM[f][m][co]->SetBinContent(b,vVsM); | |
485 | refFlowVsM[f][m][co]->SetBinError(b,vErrorVsM); | |
486 | } // end of for(Int_t b=1;b<=nBins;b++) | |
487 | } else | |
488 | { | |
489 | cout<<endl; | |
490 | cout<<" WARNING: "<<Form("cumulantsVsM[%i][%i][%i]",f,m,co)<<" is NULL in CalculateReferenceFlowFromCumulantsVsM() !!!!"<<endl; | |
491 | cout<<endl; | |
492 | } | |
493 | } // end of for(Int_t co=0;co<4;co++) // cumulant order | |
494 | } // end of for(Int_t m=0;m<nMethods;m++) | |
495 | } // end of for(Int_t f=0;f<nFiles;f++) | |
496 | ||
497 | } // end of void CalculateReferenceFlowFromCumulantsVsM() | |
498 | ||
499 | // ===================================================================================== | |
500 | ||
9dcf9913 | 501 | void Print() |
502 | { | |
503 | // Print number of events and average multiplicities for each common output file. | |
504 | ||
505 | cout<<endl; | |
506 | cout<<"Accessed files:"<<endl; | |
507 | cout<<endl; | |
508 | for(Int_t f=0;f<nFiles;f++) | |
509 | { | |
510 | cout<<commonOutputFiles[f]->GetName()<<endl; | |
511 | for(Int_t m=0;m<nMethods;m++) | |
512 | { | |
63cdf10a | 513 | AliFlowCommonHist *commonHist = NULL; |
514 | if(lists[f][m]) | |
515 | { | |
516 | commonHist = dynamic_cast<AliFlowCommonHist*> (lists[f][m]->FindObject(Form("AliFlowCommonHist%s",method[m].Data()))); | |
517 | } | |
518 | Double_t nEvts = 0.; | |
519 | Double_t AvM = 0.; | |
9dcf9913 | 520 | if(commonHist && commonHist->GetHistMultRP()) |
521 | { | |
522 | nEvts = commonHist->GetHistMultRP()->GetEntries(); | |
523 | AvM = commonHist->GetHistMultRP()->GetMean(); | |
524 | } | |
525 | if(!(strcmp(method[m].Data(),"QC"))) | |
526 | { | |
527 | cout<<Form("%s:",method[m].Data())<<" <M> = "<<AvM<<", N = "<<nEvts<<endl; | |
528 | } | |
529 | if(!(strcmp(method[m].Data(),"GFC")) && showAlsoGFCResults) | |
530 | { | |
531 | cout<<Form("%s:",method[m].Data())<<" <M> = "<<AvM<<", N = "<<nEvts<<endl; | |
532 | } | |
63cdf10a | 533 | } // end of for(Int_t m=0;m<nMethods;m++) |
9dcf9913 | 534 | cout<<endl; |
535 | } // end of for(Int_t f=0;f<nFiles;f++) | |
536 | ||
537 | } // end of void Print() | |
538 | ||
539 | // ===================================================================================== | |
540 | ||
541 | void DetermineMinMax() | |
542 | { | |
543 | // Determine ranges for plots. | |
544 | ||
545 | for(Int_t co=0;co<4;co++) | |
546 | { | |
547 | xMin[co] = 0.; yMin[co] = 44.; | |
548 | xMax[co] = -440000.; yMax[co] = -44.; | |
549 | } | |
550 | ||
551 | Double_t tfc[nFlowValues][4] = {{0.}}; // theoretical flow contribution | |
552 | for(Int_t fv=0;fv<nFlowValues;fv++) | |
553 | { | |
554 | tfc[fv][0] = pow(v[fv],2); | |
555 | tfc[fv][1] = -pow(v[fv],4); | |
556 | tfc[fv][2] = 4.*pow(v[fv],6); | |
557 | tfc[fv][3] = -33.*pow(v[fv],8); | |
558 | } | |
559 | ||
560 | for(Int_t f=0;f<nFiles;f++) | |
561 | { | |
a663ab5a | 562 | for(Int_t m=0;m<1+(Int_t)(showAlsoGFCResults)+(Int_t)(showAlsoMCEPResults);m++) |
9dcf9913 | 563 | { |
564 | for(Int_t co=0;co<4;co++) | |
565 | { | |
566 | if(cumulantsVsM[f][m][co]) | |
567 | { | |
63cdf10a | 568 | Int_t nBins = cumulantsVsM[f][m][co]->GetXaxis()->GetNbins(); |
569 | for(Int_t b=1;b<=nBins;b++) | |
9dcf9913 | 570 | { |
571 | Double_t result = cumulantsVsM[f][m][co]->GetBinContent(b); | |
572 | Double_t error = cumulantsVsM[f][m][co]->GetBinError(b); | |
573 | if(TMath::Abs(result)>1.e-44 && TMath::Abs(error)>1.e-44) | |
574 | { | |
575 | // y-axis: | |
576 | if(yMin[co] > result-error){yMin[co] = result-error;} // min value | |
a663ab5a | 577 | if(yMax[co] < result+error){yMax[co] = result+error;} // max value |
9dcf9913 | 578 | // x-axis: |
a663ab5a | 579 | if(xMax[co] < cumulantsVsM[f][m][co]->GetBinLowEdge(b+1)){xMax[co] = cumulantsVsM[f][m][co]->GetBinLowEdge(b+1);} |
9dcf9913 | 580 | } |
581 | } // end of for(Int_t b=1;b<=cumulantsVsM[f][m][co]->GetXaxis()->GetNbins();b++) | |
582 | // theoretical contributions: | |
63cdf10a | 583 | if(showTheoreticalLines) |
9dcf9913 | 584 | { |
63cdf10a | 585 | for(Int_t fv=0;fv<nFlowValues;fv++) |
586 | { | |
587 | if(yMin[co] > tfc[fv][co]) {yMin[co] = tfc[fv][co];} // min value | |
588 | if(yMax[co] < tfc[fv][co]) {yMax[co] = tfc[fv][co];} // max value | |
589 | } // end of for(Int_t fv=0;fv<nFlowValues;fv++) | |
590 | } // end of if(showTheoreticalLines) | |
9dcf9913 | 591 | } // end of if(cumulantsVsM[f][m][co]) |
592 | } // end of for(Int_t co=0;co<4;co++) | |
593 | } // end of for(Int_t m=0;m<nMethods;m++) | |
594 | } // end of for(Int_t f=0;f<nFiles;f++) | |
595 | ||
596 | } // end of void DetermineMinMax() | |
597 | ||
598 | // ===================================================================================== | |
599 | ||
600 | void GetHistograms() | |
601 | { | |
a663ab5a | 602 | // Get all histograms and profiles. |
603 | ||
604 | // a) Get profiles holding results for <refMult> vs number of Reference Particles (RPs); | |
16eed38e | 605 | // b) Get histograms holding results for cumulants and reference flow vs multiplicity. |
c109b8c7 | 606 | |
a663ab5a | 607 | // a) Get profiles holding results for <refMult> vs number of Reference Particles (RPs): |
608 | if(plotCumulantsVsReferenceMultiplicity){GetProfileRefMultVsNoOfRPs();} | |
9dcf9913 | 609 | |
16eed38e | 610 | // b) Get histograms holding results for cumulants and reference flow vs multiplicity: |
9dcf9913 | 611 | TString qcFlag[4] = {"QC{2}","QC{4}","QC{6}","QC{8}"}; |
612 | TString gfcFlag[4] = {"GFC{2}","GFC{4}","GFC{6}","GFC{8}"}; | |
613 | for(Int_t f=0;f<nFiles;f++) | |
614 | { | |
615 | for(Int_t m=0;m<nMethods;m++) | |
616 | { | |
617 | TList *temp = NULL; | |
63cdf10a | 618 | if(!(strcmp(method[m].Data(),"QC")) && lists[f][m]) |
9dcf9913 | 619 | { |
620 | temp = dynamic_cast<TList*> (lists[f][m]->FindObject("Integrated Flow")); | |
621 | if(temp) {temp = dynamic_cast<TList*> (temp->FindObject("Results"));} | |
622 | if(temp) | |
623 | { | |
624 | for(Int_t co=0;co<4;co++) | |
625 | { | |
16eed38e | 626 | // Cumulants vs multiplicity: |
9dcf9913 | 627 | cumulantsVsM[f][m][co] = dynamic_cast<TH1D*> (temp->FindObject(Form("fIntFlowQcumulantsVsM, %s",qcFlag[co].Data()))); |
dfb33a36 | 628 | if(plotCumulantsVsReferenceMultiplicity && cumulantsVsM[f][m][co]) |
629 | { | |
630 | cumulantsVsM[f][m][co] = Map(cumulantsVsM[f][m][co],f); | |
631 | } | |
6324dcfd | 632 | if(rebin && cumulantsVsM[f][m][co]) |
63cdf10a | 633 | { |
83cac266 | 634 | cumulantsVsM[f][m][co] = Rebin(cumulantsVsM[f][m][co],co); |
63cdf10a | 635 | } |
a663ab5a | 636 | } // end of for(Int_t co=0;co<4;co++) |
637 | } // end of if(temp) | |
9dcf9913 | 638 | } // end of if(!(strcmp(method[m].Data(),"QC"))) |
63cdf10a | 639 | else if(!(strcmp(method[m].Data(),"GFC")) && lists[f][m]) |
9dcf9913 | 640 | { |
641 | temp = dynamic_cast<TList*> (lists[f][m]->FindObject("Reference Flow")); | |
642 | if(temp) {temp = dynamic_cast<TList*> (temp->FindObject("Results"));} | |
643 | if(temp) | |
644 | { | |
645 | for(Int_t co=0;co<4;co++) | |
646 | { | |
647 | cumulantsVsM[f][m][co] = dynamic_cast<TH1D*> (temp->FindObject(Form("fReferenceFlowCumulantsVsM, %s",gfcFlag[co].Data()))); | |
a663ab5a | 648 | if(showAlsoGFCResults && !cumulantsVsM[f][m][co]) |
c109b8c7 | 649 | { |
a663ab5a | 650 | cout<<endl; |
651 | cout<<Form(" WARNING: Couldn't access histogram fReferenceFlowCumulantsVsM, %s in the ",gfcFlag[co].Data())<<endl; | |
652 | cout<<" file "<<commonOutputFiles[f]->GetName()<<" !!!!"<<endl; | |
653 | cout<<" Did you enable calculation of GFC vs M in the analysis which produced this file?"<<endl; | |
654 | cout<<endl; | |
c109b8c7 | 655 | } |
dfb33a36 | 656 | if(plotCumulantsVsReferenceMultiplicity && cumulantsVsM[f][m][co]) |
657 | { | |
658 | cumulantsVsM[f][m][co] = Map(cumulantsVsM[f][m][co],f); | |
a663ab5a | 659 | } |
660 | if(rebin && cumulantsVsM[f][m][co]) | |
661 | { | |
662 | cumulantsVsM[f][m][co] = Rebin(cumulantsVsM[f][m][co],co); | |
663 | } | |
63cdf10a | 664 | } // end of for(Int_t co=0;co<4;co++) |
665 | } // end of if(temp) | |
dfb33a36 | 666 | } // end of else if(!(strcmp(method[m].Data(),"GFC"))) |
667 | else if(!(strcmp(method[m].Data(),"MCEP")) && lists[f][m]) | |
668 | { | |
669 | TProfile *mcepVsM = dynamic_cast<TProfile*> (lists[f][m]->FindObject("FlowPro_VsM_MCEP")); | |
670 | if(!mcepVsM) | |
671 | { | |
672 | cout<<endl; | |
673 | cout<<" WARNING: Couldn't access TProfile FlowPro_VsM_MCEP in the file "<<commonOutputFiles[f]->GetName()<<" !!!!"<<endl; | |
674 | cout<<endl; | |
675 | return; | |
676 | } | |
677 | for(Int_t co=0;co<4;co++) | |
678 | { | |
679 | cumulantsVsM[f][m][co] = new TH1D("",Form("Monte Carlo QC{%i} #font[72]{vs} M",2*(co+1)), | |
680 | mcepVsM->GetNbinsX(),mcepVsM->GetBinLowEdge(1),mcepVsM->GetBinLowEdge(1+mcepVsM->GetNbinsX())); | |
681 | } // end of for(Int_t co=0;co<4;co++) | |
682 | for(Int_t b=1;b<=mcepVsM->GetNbinsX();b++) | |
683 | { | |
684 | Double_t v = mcepVsM->GetBinContent(b); | |
685 | Double_t vError = mcepVsM->GetBinError(b); | |
686 | if(TMath::Abs(v)<1.e-44 || TMath::Abs(vError)<1.e-44){continue;} | |
687 | // Monte Carlo QC{2}: | |
688 | Double_t qc2 = pow(v,2.); | |
689 | Double_t qc2Error = 2.*TMath::Abs(v)*vError; // error propagation for f(x) = x^2 | |
690 | cumulantsVsM[f][m][0]->SetBinContent(b,qc2); | |
691 | cumulantsVsM[f][m][0]->SetBinError(b,qc2Error); | |
692 | // Monte Carlo QC{4}: | |
693 | Double_t qc4 = -pow(v,4.); | |
694 | Double_t qc4Error = 4.*TMath::Abs(pow(v,3.))*vError; // error propagation for f(x) = -x^4 | |
695 | cumulantsVsM[f][m][1]->SetBinContent(b,qc4); | |
696 | cumulantsVsM[f][m][1]->SetBinError(b,qc4Error); | |
697 | // Monte Carlo QC{6}: | |
698 | Double_t qc6 = 4.*pow(v,6.); | |
699 | Double_t qc6Error = 24.*TMath::Abs(pow(v,5.))*vError; // error propagation for f(x) = 4x^6 | |
700 | cumulantsVsM[f][m][2]->SetBinContent(b,qc6); | |
701 | cumulantsVsM[f][m][2]->SetBinError(b,qc6Error); | |
702 | // Monte Carlo QC{8}: | |
703 | Double_t qc8 = -33.*pow(v,8.); | |
704 | Double_t qc8Error = 264.*TMath::Abs(pow(v,7.))*vError; // error propagation for f(x) = -33x^8 | |
705 | cumulantsVsM[f][m][3]->SetBinContent(b,qc8); | |
706 | cumulantsVsM[f][m][3]->SetBinError(b,qc8Error); | |
707 | } // end of for(Int_t b=1;b<=mcepVsM->GetNbinsX();b++) | |
708 | for(Int_t co=0;co<4;co++) | |
709 | { | |
710 | if(plotCumulantsVsReferenceMultiplicity && cumulantsVsM[f][m][co]) | |
711 | { | |
712 | cumulantsVsM[f][m][co] = Map(cumulantsVsM[f][m][co],f); | |
713 | } | |
714 | if(rebin && cumulantsVsM[f][m][co]) | |
715 | { | |
716 | cumulantsVsM[f][m][co] = Rebin(cumulantsVsM[f][m][co],co); | |
717 | } | |
718 | } // end of for(Int_t co=0;co<4;co++) | |
719 | } // end of else if(!(strcmp(method[m].Data(),"MCEP"))) | |
9dcf9913 | 720 | } // end of for(Int_t m=0;m<nMethods;m++) |
721 | } // end of for(Int_t f=0;f<nFiles;f++) | |
722 | ||
6324dcfd | 723 | Int_t counter = 0; |
724 | for(Int_t f=0;f<nFiles;f++) | |
725 | { | |
726 | for(Int_t m=0;m<nMethods;m++) | |
727 | { | |
728 | for(Int_t co=0;co<4;co++) | |
729 | { | |
730 | if(cumulantsVsM[f][m][co]){counter++;} | |
731 | } | |
732 | } | |
733 | } | |
734 | ||
735 | if(counter == 0) | |
736 | { | |
737 | cout<<endl; | |
738 | cout<<" WARNING: Couldn't access a single histogram with results vs multiplicity !!!!"<<endl; | |
739 | cout<<" Did you enable this calculation before running over data?"<<endl; | |
740 | cout<<endl; | |
741 | exit(0); | |
742 | } | |
743 | ||
9dcf9913 | 744 | } // end of void GetHistograms() |
745 | ||
746 | // ===================================================================================== | |
747 | ||
a663ab5a | 748 | void GetProfileRefMultVsNoOfRPs() |
c109b8c7 | 749 | { |
a663ab5a | 750 | // Get profiles holding results for <reference multiplicity> vs number of Reference Particles (RPs). |
c109b8c7 | 751 | |
a663ab5a | 752 | // Set here from which method's output file this profile will be accessed: |
753 | TString methodName = "QC"; // Alternatives are GFC and MCEP | |
754 | Int_t i = -1; | |
755 | for(Int_t m=0;m<nMethods;m++) | |
756 | { | |
757 | if(method[m] == methodName){i=m;} | |
758 | } | |
759 | if(i==-1) | |
760 | { | |
761 | cout<<endl; | |
762 | cout<<" WARNING: Unknown method name in GetProfileRefMultVsNoOfRPs() !!!!"<<endl; | |
763 | cout<<" Try something else for TString methodName in GetProfileRefMultVsNoOfRPs()."<<endl; | |
764 | cout<<endl;exit(0); | |
765 | } | |
766 | ||
63cdf10a | 767 | for(Int_t f=0;f<nFiles;f++) |
c109b8c7 | 768 | { |
63cdf10a | 769 | AliFlowCommonHist *commonHist = NULL; |
a663ab5a | 770 | if(lists[f][i]) |
63cdf10a | 771 | { |
a663ab5a | 772 | commonHist = dynamic_cast<AliFlowCommonHist*> (lists[f][i]->FindObject(Form("AliFlowCommonHist%s",method[i].Data()))); |
63cdf10a | 773 | } else |
774 | { | |
775 | cout<<endl; | |
a663ab5a | 776 | cout<<Form(" WARNING: lists[%i][%i] is NULL in GetProfileRefMultVsNoOfRPs() !!!!",f,i)<<endl; |
777 | cout<<Form(" Did you use method %s in the analysis which has produced output",method[i].Data())<<endl; | |
778 | cout<<Form(" file %s?",commonOutputFiles[f]->GetName())<<endl; | |
63cdf10a | 779 | } |
780 | if(commonHist && commonHist->GetRefMultVsNoOfRPs()) | |
781 | { | |
782 | refMultVsNoOfRPs[f] = commonHist->GetRefMultVsNoOfRPs(); | |
783 | } else | |
784 | { | |
785 | cout<<endl; | |
a663ab5a | 786 | cout<<" WARNING: commonHist && commonHist->GetRefMultVsNoOfRPs() is NULL in GetProfileRefMultVsNoOfRPs() !!!!"<<endl; |
63cdf10a | 787 | cout<<endl; |
788 | } | |
a663ab5a | 789 | } // end of for(Int_t f=0;f<nFiles;f++) |
dfb33a36 | 790 | |
791 | return; | |
c109b8c7 | 792 | |
a663ab5a | 793 | } // end of void GetProfileRefMultVsNoOfRPs() |
c109b8c7 | 794 | |
795 | // ===================================================================================== | |
796 | ||
dfb33a36 | 797 | TH1D* Map(TH1D *hist, Int_t f) |
c109b8c7 | 798 | { |
799 | // Map cumulant versus # of RPs into cumulants versus <reference multiplicity>. | |
800 | ||
63cdf10a | 801 | if(!refMultVsNoOfRPs[f]) |
802 | { | |
803 | cout<<endl; | |
a663ab5a | 804 | cout<<Form(" WARNING: refMultVsNoOfRPs[%i] is NULL in Map(...) !!!!",f)<<endl; |
63cdf10a | 805 | cout<<endl; |
806 | } | |
807 | ||
a663ab5a | 808 | Int_t rpMinBin = refMultVsNoOfRPs[f]->FindFirstBinAbove(); // FindFirstBinAbove(Double_t threshold = 0, Int_t axis = 1) |
dfb33a36 | 809 | Int_t rpMaxBin = refMultVsNoOfRPs[f]->FindLastBinAbove(); // FindLastBinAbove(Double_t threshold = 0, Int_t axis = 1) |
a663ab5a | 810 | Int_t rmMinBin = 440000;; |
dfb33a36 | 811 | Int_t rmMaxBin = (Int_t)TMath::Floor(refMultVsNoOfRPs[f]->GetMaximum()); |
a663ab5a | 812 | for(Int_t rpBin=rpMinBin;rpBin<=rpMaxBin;rpBin++) // non-empty # of RPs bins |
813 | { | |
814 | if(refMultVsNoOfRPs[f]->GetBinContent(rpBin)>0. && refMultVsNoOfRPs[f]->GetBinContent(rpBin)<rmMinBin) | |
815 | { | |
816 | rmMinBin = (Int_t)TMath::Floor(refMultVsNoOfRPs[f]->GetBinContent(rpBin)); | |
817 | } | |
818 | } // end of for(Int_t rpBin=rpMinBin;rpBin<=rpMaxBin;rpBin++) // non-empty # of RPs bins | |
dfb33a36 | 819 | |
c109b8c7 | 820 | if(hist) |
821 | { | |
822 | temp = (TH1D*) hist->Clone(); | |
dfb33a36 | 823 | temp->Reset(); |
a663ab5a | 824 | for(Int_t rmBin=rmMinBin;rmBin<=rmMaxBin;rmBin++) // reference multiplicity bins |
dfb33a36 | 825 | { |
826 | Double_t value = 0.; | |
827 | Double_t error = 0.; | |
828 | Double_t dSum1 = 0.; // sum value_i/(error_i)^2 | |
829 | Double_t dSum2 = 0.; // sum 1/(error_i)^2 | |
a663ab5a | 830 | for(Int_t rpBin=rpMinBin;rpBin<=rpMaxBin;rpBin++) // # of RPs bins |
dfb33a36 | 831 | { |
832 | if((Int_t)TMath::Floor(refMultVsNoOfRPs[f]->GetBinContent(rpBin)) >= temp->GetBinLowEdge(rmBin) && | |
833 | (Int_t)TMath::Floor(refMultVsNoOfRPs[f]->GetBinContent(rpBin)) < temp->GetBinLowEdge(rmBin+1)) | |
834 | { | |
835 | value = hist->GetBinContent(rpBin); | |
836 | error = hist->GetBinError(rpBin); | |
837 | if(error>0.) | |
838 | { | |
839 | dSum1+=value/(error*error); | |
840 | dSum2+=1./(error*error); | |
841 | } | |
842 | } | |
843 | } // end of for(Int_t rpBin=1;rpBin<=nBins;rpBin++) // # of RPs bins | |
844 | if(dSum2>0.) | |
845 | { | |
846 | temp->SetBinContent(rmBin,dSum1/dSum2); | |
847 | temp->SetBinError(rmBin,pow(1./dSum2,0.5)); | |
848 | } | |
849 | } // end of for(Int_t rmBin=1;rmBin<=nBins;rmBin++) // reference multiplicity bins | |
850 | } // end of if(hist) | |
851 | ||
852 | return temp; | |
63cdf10a | 853 | |
854 | } // end of Map() | |
855 | ||
856 | // ===================================================================================== | |
857 | ||
83cac266 | 858 | TH1D* Rebin(TH1D *hist, Int_t co) |
63cdf10a | 859 | { |
860 | // Rebin original histograms. | |
861 | ||
83cac266 | 862 | if(nMergedBins[co] == 0) |
63cdf10a | 863 | { |
864 | cout<<endl; | |
83cac266 | 865 | cout<<Form(" WARNING: nMergedBins[%i] == 0 !!!!",co)<<endl; |
63cdf10a | 866 | cout<<endl; |
867 | exit(0); | |
868 | } | |
869 | if(!hist) | |
870 | { | |
871 | cout<<endl; | |
872 | cout<<" WARNING: hist is NULL in Rebin() !!!!"<<endl; | |
873 | cout<<endl; | |
874 | exit(0); | |
875 | } | |
876 | ||
877 | Int_t nBinsOld = hist->GetXaxis()->GetNbins(); | |
878 | if(nBinsOld == 0){cout<<" WARNING: nBinsOld == 0 !!!!"<<endl;exit(0);} | |
879 | Double_t xMinOld = hist->GetXaxis()->GetXmin(); | |
880 | Double_t xMaxOld = hist->GetXaxis()->GetXmax(); | |
881 | Double_t binWidthOld = (xMaxOld-xMinOld)/nBinsOld; | |
83cac266 | 882 | Int_t nBinsNew = TMath::Floor(nBinsOld/nMergedBins[co]); |
63cdf10a | 883 | Double_t xMinNew = xMinOld; |
83cac266 | 884 | Double_t xMaxNew = xMinOld + nBinsNew*nMergedBins[co]*binWidthOld; |
c109b8c7 | 885 | |
63cdf10a | 886 | TH1D *temp = new TH1D("","",nBinsNew,xMinNew,xMaxNew); // rebinned histogram |
887 | Int_t binNew = 1; | |
888 | Double_t value = 0.; | |
889 | Double_t error = 0.; | |
890 | Double_t dSum1 = 0.; // sum value_i/(error_i)^2 | |
891 | Double_t dSum2 = 0.; // sum 1/(error_i)^2 | |
892 | for(Int_t b=1;b<=nBinsOld;b++) | |
893 | { | |
894 | value = hist->GetBinContent(b); | |
895 | error = hist->GetBinError(b); | |
896 | if(error>0.) | |
897 | { | |
898 | dSum1+=value/(error*error); | |
899 | dSum2+=1./(error*error); | |
900 | } | |
83cac266 | 901 | if(b%nMergedBins[co] == 0) |
63cdf10a | 902 | { |
903 | if(dSum2>0.) | |
904 | { | |
905 | temp->SetBinContent(binNew,dSum1/dSum2); | |
906 | temp->SetBinError(binNew,pow(1./dSum2,0.5)); | |
907 | } | |
908 | binNew++; | |
909 | dSum1 = 0.; | |
910 | dSum2 = 0.; | |
83cac266 | 911 | } // end of if(b%nMergedBins[co] == 0) |
63cdf10a | 912 | } // end of for(Int_t b=1;b<=nBinsOld;b++) |
c109b8c7 | 913 | |
63cdf10a | 914 | return temp; |
915 | ||
83cac266 | 916 | } // end of TH1D* Rebin(TH1D *hist, Int_t co) |
c109b8c7 | 917 | |
918 | // ===================================================================================== | |
919 | ||
9dcf9913 | 920 | TGraphErrors* GetGraphErrors(Int_t bin, Int_t nFiles, TH1D** qc) |
921 | { | |
922 | TGraphErrors *ge = new TGraphErrors(nFiles); | |
923 | for(Int_t f=0;f<nFiles;f++) | |
924 | { | |
925 | ge->SetPoint(f,f+0.5,qc[f]->GetBinContent(bin+1)); | |
926 | ge->SetPointError(f,0,qc[f]->GetBinError(bin+1)); | |
927 | } | |
928 | ||
929 | return ge; | |
930 | ||
931 | } // end of TGraphErrors* GetGraphErrors(Int_t bin, Int_t nFiles, TH1D** qc) | |
932 | ||
933 | // ===================================================================================== | |
934 | ||
935 | TH1D* StyleHist(TString yAxisTitle, Int_t co) | |
936 | { | |
937 | // Style histogram. | |
938 | ||
16eed38e | 939 | TH1D *styleHist = new TH1D(yAxisTitle.Data(),"",(Int_t)xMax[co],0,xMax[co]); |
9dcf9913 | 940 | // y-axis: |
941 | styleHist->GetYaxis()->SetRangeUser(yMin[co],yMax[co]); | |
63cdf10a | 942 | |
c109b8c7 | 943 | if(plotCumulantsVsReferenceMultiplicity) |
944 | { | |
945 | styleHist->GetXaxis()->SetTitle("#LTreference multiplicity#GT"); | |
946 | } else | |
947 | { | |
948 | styleHist->GetXaxis()->SetTitle("# of RPs"); | |
949 | } | |
950 | ||
9dcf9913 | 951 | styleHist->GetYaxis()->SetTitle(yAxisTitle.Data()); |
952 | ||
953 | return styleHist; | |
954 | ||
955 | } // end of TH1D* StyleHist(TString yAxisTitle, Int_t co) | |
956 | ||
957 | // =========================================================================================== | |
958 | ||
959 | TGraph* GetErrorMesh(TH1D *hist) | |
960 | { | |
961 | // Error mesh. | |
962 | ||
9dcf9913 | 963 | if(hist) |
964 | { | |
965 | Int_t nBins = hist->GetNbinsX(); | |
dfb33a36 | 966 | Double_t value = 0.; |
967 | Double_t error = 0.; | |
9dcf9913 | 968 | // Counting the non-empty bins: |
dfb33a36 | 969 | Int_t nNonEmptyBins = 0; |
9dcf9913 | 970 | for(Int_t i=1;i<nBins+1;i++) |
971 | { | |
dfb33a36 | 972 | value = hist->GetBinContent(i); |
973 | error = hist->GetBinError(i); | |
974 | if(TMath::Abs(value)>0.0 && error>0.0) | |
9dcf9913 | 975 | { |
976 | nNonEmptyBins++; | |
977 | } | |
dfb33a36 | 978 | } // end of for(Int_t i=1;i<nBins+1;i++) |
979 | // Error mesh: | |
980 | TGraph *errorMesh = new TGraph(2*nNonEmptyBins+1); | |
981 | Int_t count = 0; | |
982 | Double_t binCenter = 0.; | |
9dcf9913 | 983 | for(Int_t i=1;i<nBins+1;i++) |
984 | { | |
9dcf9913 | 985 | value = hist->GetBinContent(i); |
dfb33a36 | 986 | error = hist->GetBinError(i); |
987 | // Setting up the the mesh: | |
988 | if(TMath::Abs(value)>0.0 && error>0.0) | |
9dcf9913 | 989 | { |
dfb33a36 | 990 | binCenter = hist->GetBinCenter(i); |
991 | errorMesh->SetPoint(count,binCenter,value+error); | |
992 | errorMesh->SetPoint(2*nNonEmptyBins-(count+1),binCenter,value-error); | |
993 | count++; | |
994 | } | |
995 | } // end of for(Int_t i=1;i<nBins+1;i++) | |
9dcf9913 | 996 | // Closing the mesh area: |
dfb33a36 | 997 | Double_t xStart = 0.; |
998 | Double_t yStart = 0.; | |
999 | errorMesh->GetPoint(0,xStart,yStart); | |
1000 | errorMesh->SetPoint(2*nNonEmptyBins,xStart,yStart); | |
9dcf9913 | 1001 | } // end if(hist) |
1002 | ||
9dcf9913 | 1003 | return errorMesh; |
1004 | ||
1005 | } // end of TGraph* GetErrorMesh(TH1D *hist) | |
1006 | ||
1007 | // =========================================================================================== | |
1008 | ||
1009 | void GlobalSettings() | |
1010 | { | |
1011 | // Settings which will affect all plots. | |
1012 | ||
1013 | gROOT->SetStyle("Plain"); // default color is white instead of gray | |
1014 | gStyle->SetOptStat(0); // remove stat. box from all histos | |
63cdf10a | 1015 | TGaxis::SetMaxDigits(4); |
9dcf9913 | 1016 | |
1017 | } // end of void GlobalSettings() | |
1018 | ||
1019 | // =========================================================================================== | |
1020 | ||
1021 | void GetLists() | |
1022 | { | |
1023 | // Get from common output files the lists holding histograms for each method. | |
1024 | ||
1025 | TString fileName[nFiles][nMethods]; | |
1026 | TDirectoryFile *dirFile[nFiles][nMethods] = {{NULL}}; | |
1027 | TString listName[nFiles][nMethods]; | |
1028 | for(Int_t f=0;f<nFiles;f++) | |
1029 | { | |
1030 | for(Int_t i=0;i<nMethods;i++) | |
1031 | { | |
1032 | // Form a file name for each method: | |
1033 | fileName[f][i]+="output"; | |
1034 | fileName[f][i]+=method[i].Data(); | |
1035 | fileName[f][i]+="analysis"; | |
1036 | fileName[f][i]+=type[f].Data(); | |
1037 | // Access this file: | |
1038 | if(commonOutputFiles[f]){dirFile[f][i] = (TDirectoryFile*)commonOutputFiles[f]->FindObjectAny(fileName[f][i].Data());} | |
1039 | // Form a list name for each method: | |
1040 | listName[f][i]+="cobj"; | |
1041 | listName[f][i]+=method[i].Data(); | |
1042 | // Access this lists: | |
1043 | if(dirFile[f][i]) | |
1044 | { | |
1045 | dirFile[f][i]->GetObject(listName[f][i].Data(),lists[f][i]); | |
63cdf10a | 1046 | if(!lists[f][i]) |
1047 | { | |
1048 | cout<<endl; | |
1049 | cout<<" WARNING: Couldn't find a list "<<listName[f][i].Data()<<" in "<<commonOutputFiles[f]->GetName()<<" !!!!"<<endl; | |
1050 | cout<<" Did you use method "<<method[i].Data()<<" in the analysis which produced this file?"<<endl; | |
63cdf10a | 1051 | } |
9dcf9913 | 1052 | } else |
1053 | { | |
63cdf10a | 1054 | cout<<endl; |
dfb33a36 | 1055 | cout<<" WARNING: Couldn't find a file "<<fileName[f][i].Data()<<".root in "<<commonOutputFiles[f]->GetName()<<" !!!!"<<endl; |
63cdf10a | 1056 | cout<<endl; |
9dcf9913 | 1057 | } |
1058 | } // end of for(Int_t i=0;i<nMethods;i++) | |
1059 | } // end of for(Int_t f=0;f<nFiles;f++) | |
1060 | ||
1061 | } // end of void GetLists() | |
1062 | ||
1063 | // =========================================================================================== | |
1064 | ||
1065 | void AccessCommonOutputFiles(TString commonOutputFileName) | |
1066 | { | |
1067 | // Access all output files. | |
a663ab5a | 1068 | |
9dcf9913 | 1069 | for(Int_t f=0;f<nFiles;f++) |
1070 | { | |
1071 | if(!(gSystem->AccessPathName(Form("%s/%s/%s",gSystem->pwd(),files[f].Data(),commonOutputFileName.Data()),kFileExists))) | |
1072 | { | |
1073 | commonOutputFiles[f] = TFile::Open(Form("%s/%s/%s",gSystem->pwd(),files[f].Data(),commonOutputFileName.Data()),"READ"); | |
1074 | } else | |
1075 | { | |
1076 | cout<<endl; | |
a663ab5a | 1077 | cout<<" WARNING: Couldn't find the file "<<Form("%s/%s/%s",gSystem->pwd(),files[f].Data(),commonOutputFileName.Data())<<" !!!!"<<endl; |
1078 | cout<<" Did you specify correctly all paths in TString files[nFiles]?"<<endl; | |
9dcf9913 | 1079 | cout<<endl; |
1080 | exit(0); | |
1081 | } | |
1082 | } // end of for(Int_t f=0;f<nFiles;f++) | |
1083 | ||
1084 | } // void AccessCommonOutputFiles(TString commonOutputFileName); | |
1085 | ||
1086 | // =========================================================================================== | |
1087 | ||
1088 | void LoadLibrariesPC(const libModes analysisMode) { | |
1089 | ||
1090 | //-------------------------------------- | |
1091 | // Load the needed libraries most of them already loaded by aliroot | |
1092 | //-------------------------------------- | |
1093 | //gSystem->Load("libTree"); | |
1094 | gSystem->Load("libGeom"); | |
1095 | gSystem->Load("libVMC"); | |
1096 | gSystem->Load("libXMLIO"); | |
1097 | gSystem->Load("libPhysics"); | |
1098 | ||
1099 | //---------------------------------------------------------- | |
1100 | // >>>>>>>>>>> Local mode <<<<<<<<<<<<<< | |
1101 | //---------------------------------------------------------- | |
1102 | if (analysisMode==mLocal) { | |
1103 | //-------------------------------------------------------- | |
1104 | // If you want to use already compiled libraries | |
1105 | // in the aliroot distribution | |
1106 | //-------------------------------------------------------- | |
1107 | ||
1108 | //================================================================================== | |
1109 | //load needed libraries: | |
1110 | gSystem->AddIncludePath("-I$ROOTSYS/include"); | |
1111 | //gSystem->Load("libTree"); | |
1112 | ||
1113 | // for AliRoot | |
1114 | gSystem->AddIncludePath("-I$ALICE_ROOT/include"); | |
1115 | //gSystem->Load("libANALYSIS"); | |
1116 | gSystem->Load("libPWG2flowCommon"); | |
1117 | //cerr<<"libPWG2flowCommon loaded ..."<<endl; | |
1118 | ||
1119 | } | |
1120 | ||
1121 | else if (analysisMode==mLocalSource) { | |
1122 | ||
1123 | // In root inline compile | |
1124 | ||
3b67e1f4 | 1125 | // Constants |
9dcf9913 | 1126 | gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+"); |
1127 | gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+"); | |
9dcf9913 | 1128 | |
1129 | // Flow event | |
1130 | gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+"); | |
1131 | gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+"); | |
9dcf9913 | 1132 | gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+"); |
3b67e1f4 | 1133 | gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+"); |
1134 | ||
9dcf9913 | 1135 | // Output histosgrams |
1136 | gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+"); | |
1137 | gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+"); | |
1138 | gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+"); | |
1139 | gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+"); | |
1140 | ||
1141 | cout << "finished loading macros!" << endl; | |
1142 | ||
1143 | } | |
1144 | ||
1145 | } // end of void LoadLibraries(const libModes analysisMode) |