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