]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/macros/drawBalanceFunctionPsi.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
correct efficiency/acceptance histograms (centrality dependence)
[u/mrichter/AliRoot.git] / PWGCF / EBYE / macros / drawBalanceFunctionPsi.C
CommitLineData
db7174c0 1const Int_t numberOfCentralityBins = 10;
2TString centralityArray[numberOfCentralityBins] = {"0-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-1","1-2"};
a38fd7f3 3
a38fd7f3 4void drawBalanceFunctionPsi(const char* filename = "AnalysisResultsPsi.root",
eb63b883 5 Int_t gCentrality = 1,
6 Int_t gDeltaEtaDeltaPhi = 2,
6acdbcb2 7 Double_t psiMin = -0.5, Double_t psiMax = 0.5,
8 Double_t ptTriggerMin = -1.,
9 Double_t ptTriggerMax = -1.,
10 Double_t ptAssociatedMin = -1.,
11 Double_t ptAssociatedMax = -1.) {
a38fd7f3 12 //Macro that draws the BF distributions for each centrality bin
13 //for reaction plane dependent analysis
14 //Author: Panos.Christakoglou@nikhef.nl
15 //Load the PWG2ebye library
16 gSystem->Load("libANALYSIS.so");
17 gSystem->Load("libANALYSISalice.so");
18 gSystem->Load("libEventMixing.so");
19 gSystem->Load("libCORRFW.so");
20 gSystem->Load("libPWGTools.so");
21 gSystem->Load("libPWGCFebye.so");
22
23 //Prepare the objects and return them
eb63b883 24 TList *listBF = GetListOfObjects(filename,gCentrality,0);
25 TList *listBFShuffled = GetListOfObjects(filename,gCentrality,1);
26 TList *listBFMixed = GetListOfObjects(filename,gCentrality,2);
a38fd7f3 27 if(!listBF) {
28 Printf("The TList object was not created");
29 return;
30 }
31 else
eb63b883 32 draw(listBF,listBFShuffled,listBFMixed,
33 gCentrality,gDeltaEtaDeltaPhi,
34 psiMin,psiMax,
6acdbcb2 35 ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
a38fd7f3 36}
37
38//______________________________________________________//
eb63b883 39TList *GetListOfObjects(const char* filename,
40 Int_t gCentrality, Int_t kData = 1) {
a38fd7f3 41 //Get the TList objects (QA, bf, bf shuffled)
eb63b883 42 //kData == 0: data
43 //kData == 1: shuffling
44 //kData == 2: mixing
a38fd7f3 45 TList *listQA = 0x0;
46 TList *listBF = 0x0;
a38fd7f3 47
48 //Open the file
49 TFile *f = TFile::Open(filename);
50 if((!f)||(!f->IsOpen())) {
51 Printf("The file %s is not found. Aborting...",filename);
52 return listBF;
53 }
f18dd679 54 f->ls();
a38fd7f3 55
6acdbcb2 56 TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
a38fd7f3 57 if(!dir) {
58 Printf("The TDirectoryFile is not found. Aborting...",filename);
59 return listBF;
60 }
f18dd679 61 dir->ls();
a38fd7f3 62
63 TString listBFName;
eb63b883 64 if(kData == 0) {
65 //cout<<"no shuffling - no mixing"<<endl;
66 listBFName = "listBFPsi_";
67 }
68 else if(kData == 1) {
69 //cout<<"shuffling - no mixing"<<endl;
70 listBFName = "listBFPsiShuffled_";
71 }
72 else if(kData == 2) {
73 //cout<<"no shuffling - mixing"<<endl;
74 listBFName = "listBFPsiMixed_";
75 }
76 listBFName += centralityArray[gCentrality-1];
f18dd679 77 listBFName += "_Bit128_V0M";
a38fd7f3 78 listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
eb63b883 79 cout<<"======================================================="<<endl;
80 cout<<"List name: "<<listBF->GetName()<<endl;
f18dd679 81 listBF->ls();
a38fd7f3 82
83 //Get the histograms
84 TString histoName;
eb63b883 85 if(kData == 0)
a38fd7f3 86 histoName = "fHistPV0M";
eb63b883 87 else if(kData == 1)
a38fd7f3 88 histoName = "fHistP_shuffleV0M";
eb63b883 89 else if(kData == 2)
90 histoName = "fHistPV0M";
a38fd7f3 91 AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
92 if(!fHistP) {
93 Printf("fHistP %s not found!!!",histoName.Data());
94 break;
95 }
96 fHistP->FillParent(); fHistP->DeleteContainers();
97
eb63b883 98 if(kData == 0)
a38fd7f3 99 histoName = "fHistNV0M";
eb63b883 100 if(kData == 1)
a38fd7f3 101 histoName = "fHistN_shuffleV0M";
eb63b883 102 if(kData == 2)
103 histoName = "fHistNV0M";
a38fd7f3 104 AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
105 if(!fHistN) {
106 Printf("fHistN %s not found!!!",histoName.Data());
107 break;
108 }
109 fHistN->FillParent(); fHistN->DeleteContainers();
110
eb63b883 111 if(kData == 0)
a38fd7f3 112 histoName = "fHistPNV0M";
eb63b883 113 if(kData == 1)
a38fd7f3 114 histoName = "fHistPN_shuffleV0M";
eb63b883 115 if(kData == 2)
116 histoName = "fHistPNV0M";
a38fd7f3 117 AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
118 if(!fHistPN) {
119 Printf("fHistPN %s not found!!!",histoName.Data());
120 break;
121 }
122 fHistPN->FillParent(); fHistPN->DeleteContainers();
123
eb63b883 124 if(kData == 0)
a38fd7f3 125 histoName = "fHistNPV0M";
eb63b883 126 if(kData == 1)
a38fd7f3 127 histoName = "fHistNP_shuffleV0M";
eb63b883 128 if(kData == 2)
129 histoName = "fHistNPV0M";
a38fd7f3 130 AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
131 if(!fHistNP) {
132 Printf("fHistNP %s not found!!!",histoName.Data());
133 break;
134 }
135 fHistNP->FillParent(); fHistNP->DeleteContainers();
136
eb63b883 137 if(kData == 0)
a38fd7f3 138 histoName = "fHistPPV0M";
eb63b883 139 if(kData == 1)
a38fd7f3 140 histoName = "fHistPP_shuffleV0M";
eb63b883 141 if(kData == 2)
142 histoName = "fHistPPV0M";
a38fd7f3 143 AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
144 if(!fHistPP) {
145 Printf("fHistPP %s not found!!!",histoName.Data());
146 break;
147 }
148 fHistPP->FillParent(); fHistPP->DeleteContainers();
149
eb63b883 150 if(kData == 0)
a38fd7f3 151 histoName = "fHistNNV0M";
eb63b883 152 if(kData == 1)
a38fd7f3 153 histoName = "fHistNN_shuffleV0M";
eb63b883 154 if(kData == 2)
155 histoName = "fHistNNV0M";
a38fd7f3 156 AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
157 if(!fHistNN) {
158 Printf("fHistNN %s not found!!!",histoName.Data());
159 break;
160 }
161 fHistNN->FillParent(); fHistNN->DeleteContainers();
162
163 return listBF;
164}
165
166//______________________________________________________//
eb63b883 167void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
168 Int_t gCentrality, Int_t gDeltaEtaDeltaPhi,
6acdbcb2 169 Double_t psiMin, Double_t psiMax,
170 Double_t ptTriggerMin, Double_t ptTriggerMax,
171 Double_t ptAssociatedMin, Double_t ptAssociatedMax) {
a38fd7f3 172 gROOT->LoadMacro("~/SetPlotStyle.C");
173 SetPlotStyle();
174 gStyle->SetPalette(1,0);
175
eb63b883 176 const Int_t gRebin = gDeltaEtaDeltaPhi; //rebin by 2 the Delta phi projection
177
a38fd7f3 178 //balance function
179 AliTHn *hP = NULL;
180 AliTHn *hN = NULL;
181 AliTHn *hPN = NULL;
182 AliTHn *hNP = NULL;
183 AliTHn *hPP = NULL;
184 AliTHn *hNN = NULL;
185 //listBF->ls();
186 //Printf("=================");
187 hP = (AliTHn*) listBF->FindObject("fHistPV0M");
188 hN = (AliTHn*) listBF->FindObject("fHistNV0M");
189 hPN = (AliTHn*) listBF->FindObject("fHistPNV0M");
190 hNP = (AliTHn*) listBF->FindObject("fHistNPV0M");
191 hPP = (AliTHn*) listBF->FindObject("fHistPPV0M");
192 hNN = (AliTHn*) listBF->FindObject("fHistNNV0M");
193
194 AliBalancePsi *b = new AliBalancePsi();
195 b->SetHistNp(hP);
196 b->SetHistNn(hN);
197 b->SetHistNpn(hPN);
198 b->SetHistNnp(hNP);
199 b->SetHistNpp(hPP);
200 b->SetHistNnn(hNN);
201
202 //balance function shuffling
203 AliTHn *hPShuffled = NULL;
204 AliTHn *hNShuffled = NULL;
205 AliTHn *hPNShuffled = NULL;
206 AliTHn *hNPShuffled = NULL;
207 AliTHn *hPPShuffled = NULL;
208 AliTHn *hNNShuffled = NULL;
209 //listBFShuffled->ls();
210 hPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistP_shuffleV0M");
211 hNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistN_shuffleV0M");
212 hPNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPN_shuffleV0M");
213 hNPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNP_shuffleV0M");
214 hPPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPP_shuffleV0M");
215 hNNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNN_shuffleV0M");
216
217 AliBalancePsi *bShuffled = new AliBalancePsi();
218 bShuffled->SetHistNp(hPShuffled);
219 bShuffled->SetHistNn(hNShuffled);
220 bShuffled->SetHistNpn(hPNShuffled);
221 bShuffled->SetHistNnp(hNPShuffled);
222 bShuffled->SetHistNpp(hPPShuffled);
223 bShuffled->SetHistNnn(hNNShuffled);
224
eb63b883 225 //balance function mixing
226 AliTHn *hPMixed = NULL;
227 AliTHn *hNMixed = NULL;
228 AliTHn *hPNMixed = NULL;
229 AliTHn *hNPMixed = NULL;
230 AliTHn *hPPMixed = NULL;
231 AliTHn *hNNMixed = NULL;
232 //listBFMixed->ls();
233 hPMixed = (AliTHn*) listBFMixed->FindObject("fHistPV0M");
234 hNMixed = (AliTHn*) listBFMixed->FindObject("fHistNV0M");
235 hPNMixed = (AliTHn*) listBFMixed->FindObject("fHistPNV0M");
236 hNPMixed = (AliTHn*) listBFMixed->FindObject("fHistNPV0M");
237 hPPMixed = (AliTHn*) listBFMixed->FindObject("fHistPPV0M");
238 hNNMixed = (AliTHn*) listBFMixed->FindObject("fHistNNV0M");
239
240 AliBalancePsi *bMixed = new AliBalancePsi();
241 bMixed->SetHistNp(hPMixed);
242 bMixed->SetHistNn(hNMixed);
243 bMixed->SetHistNpn(hPNMixed);
244 bMixed->SetHistNnp(hNPMixed);
245 bMixed->SetHistNpp(hPPMixed);
246 bMixed->SetHistNnn(hNNMixed);
247
6acdbcb2 248 TH1D *gHistBalanceFunction;
249 TH1D *gHistBalanceFunctionShuffled;
eb63b883 250 TH1D *gHistBalanceFunctionMixed;
251 TH1D *gHistBalanceFunctionSubtracted;
a38fd7f3 252 TString histoTitle, pngName;
6acdbcb2 253 TLegend *legend;
a38fd7f3 254
6acdbcb2 255 histoTitle = "Centrality: ";
eb63b883 256 histoTitle += centralityArray[gCentrality-1];
6acdbcb2 257 histoTitle += "%";
eb63b883 258 if((psiMin == -0.5)&&(psiMax == 0.5))
259 histoTitle += " (-7.5^{o} < #phi - #Psi_{2} < 7.5^{o})";
260 else if((psiMin == 0.5)&&(psiMax == 1.5))
261 histoTitle += " (37.5^{o} < #phi - #Psi_{2} < 52.5^{o})";
262 else if((psiMin == 1.5)&&(psiMax == 2.5))
263 histoTitle += " (82.5^{o} < #phi - #Psi_{2} < 97.5^{o})";
264 else
265 histoTitle += " (0^{o} < #phi - #Psi_{2} < 180^{o})";
6acdbcb2 266
eb63b883 267 //Raw balance function
268 gHistBalanceFunction = b->GetBalanceFunctionHistogram(0,gDeltaEtaDeltaPhi,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
6acdbcb2 269 gHistBalanceFunction->SetMarkerStyle(20);
270 gHistBalanceFunction->SetTitle(histoTitle.Data());
271 gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
eb63b883 272 gHistBalanceFunction->SetName("gHistBalanceFunction");
6acdbcb2 273
eb63b883 274 //Shuffled balance function
275 gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionHistogram(0,gDeltaEtaDeltaPhi,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
6acdbcb2 276 gHistBalanceFunctionShuffled->SetMarkerStyle(24);
eb63b883 277 gHistBalanceFunctionShuffled->SetName("gHistBalanceFunctionShuffled");
278
279 //Mixed balance function
280 gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionHistogram(0,gDeltaEtaDeltaPhi,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
281 gHistBalanceFunctionMixed->SetMarkerStyle(25);
282 gHistBalanceFunctionMixed->SetName("gHistBalanceFunctionMixed");
6acdbcb2 283
eb63b883 284 //Subtracted balance function
285 gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunction->Clone());
286 gHistBalanceFunctionSubtracted->Add(gHistBalanceFunctionMixed,-1);
287 gHistBalanceFunctionSubtracted->Rebin(gRebin);
288 gHistBalanceFunctionSubtracted->SetMarkerStyle(20);
289 gHistBalanceFunctionSubtracted->SetTitle(histoTitle.Data());
290 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
291 gHistBalanceFunctionSubtracted->SetName("gHistBalanceFunctionSubtracted");
292
293 TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
6acdbcb2 294 c1->SetFillColor(10);
295 c1->SetHighLightColor(10);
296 c1->SetLeftMargin(0.15);
297 gHistBalanceFunction->DrawCopy("E");
298 gHistBalanceFunctionShuffled->DrawCopy("ESAME");
eb63b883 299 gHistBalanceFunctionMixed->DrawCopy("ESAME");
6acdbcb2 300
eb63b883 301 legend = new TLegend(0.18,0.62,0.45,0.82,"","brNDC");
6acdbcb2 302 legend->SetTextSize(0.045);
303 legend->SetTextFont(42);
304 legend->SetBorderSize(0);
305 legend->SetFillStyle(0);
306 legend->SetFillColor(10);
307 legend->SetMargin(0.25);
308 legend->SetShadowColor(10);
309 legend->AddEntry(gHistBalanceFunction,"Data","lp");
310 legend->AddEntry(gHistBalanceFunctionShuffled,"Shuffled data","lp");
eb63b883 311 legend->AddEntry(gHistBalanceFunctionMixed,"Mixed data","lp");
6acdbcb2 312 legend->Draw();
313
314 pngName = "BalanceFunctionDeltaEta.Centrality";
eb63b883 315 pngName += centralityArray[gCentrality-1];
316 pngName += ".Psi"; //pngName += psiMin; pngName += "To"; pngName += psiMax;
6acdbcb2 317 pngName += ".png";
318 c1->SaveAs(pngName.Data());
319
320 GetWeightedMean(gHistBalanceFunction);
321 GetWeightedMean(gHistBalanceFunctionShuffled);
322
323 TString meanLatex, rmsLatex, skewnessLatex, kurtosisLatex;
324 meanLatex = "#mu = ";
eb63b883 325 meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMean());
6acdbcb2 326 meanLatex += " #pm ";
eb63b883 327 meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMeanError());
6acdbcb2 328
329 rmsLatex = "#sigma = ";
eb63b883 330 rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMS());
6acdbcb2 331 rmsLatex += " #pm ";
eb63b883 332 rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMSError());
6acdbcb2 333
334 skewnessLatex = "S = ";
eb63b883 335 skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(1));
6acdbcb2 336 skewnessLatex += " #pm ";
eb63b883 337 skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(11));
6acdbcb2 338
339 kurtosisLatex = "K = ";
eb63b883 340 kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(1));
6acdbcb2 341 kurtosisLatex += " #pm ";
eb63b883 342 kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(11));
343 Printf("Mean: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetMean(),gHistBalanceFunctionSubtracted->GetMeanError());
344 Printf("RMS: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetRMS(),gHistBalanceFunctionSubtracted->GetRMSError());
345 Printf("Skeweness: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetSkewness(1),gHistBalanceFunctionSubtracted->GetSkewness(11));
346 Printf("Kurtosis: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetKurtosis(1),gHistBalanceFunctionSubtracted->GetKurtosis(11));
347
348 TCanvas *c2 = new TCanvas("c2","",600,0,600,500);
349 c2->SetFillColor(10);
350 c2->SetHighLightColor(10);
351 c2->SetLeftMargin(0.15);
352 gHistBalanceFunctionSubtracted->DrawCopy("E");
6acdbcb2 353
354 TLatex *latex = new TLatex();
355 latex->SetNDC();
356 latex->SetTextSize(0.035);
357 latex->SetTextColor(1);
358 latex->DrawLatex(0.64,0.85,meanLatex.Data());
359 latex->DrawLatex(0.64,0.81,rmsLatex.Data());
360 latex->DrawLatex(0.64,0.77,skewnessLatex.Data());
361 latex->DrawLatex(0.64,0.73,kurtosisLatex.Data());
eb63b883 362
363 TString newFileName = "balanceFunction.Centrality";
364 newFileName += gCentrality; newFileName += ".In";
365 if(gDeltaEtaDeltaPhi == 1) newFileName += "DeltaEta.Psi";
366 else if(gDeltaEtaDeltaPhi == 2) newFileName += "DeltaPhi.Psi";
367 if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
368 else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
369 else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
648f1a5a 370 else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
371 else newFileName += "All.PttFrom";
372 newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
373 newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
374 newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
375 newFileName += Form("%.1f",ptAssociatedMax);
376 newFileName += ".root";
eb63b883 377
378 TFile *fOutput = new TFile(newFileName.Data(),"recreate");
379 fOutput->cd();
380 gHistBalanceFunction->Write();
381 gHistBalanceFunctionShuffled->Write();
382 gHistBalanceFunctionMixed->Write();
383 gHistBalanceFunctionSubtracted->Write();
384 fOutput->Close();
a38fd7f3 385}
386
387//____________________________________________________________________//
388void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1) {
389 //Prints the calculated width of the BF and its error
390 Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
391 Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
392 Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
393 Double_t deltaBalP2 = 0.0, integral = 0.0;
394 Double_t deltaErrorNew = 0.0;
395
396 //Retrieve this variables from Histogram
397 Int_t fNumberOfBins = gHistBalance->GetNbinsX();
398 Double_t fP2Step = gHistBalance->GetBinWidth(1); // assume equal binning!
399
400 cout<<"=================================================="<<endl;
401 cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
402 cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
403 cout<<"=================================================="<<endl;
404 for(Int_t i = 1; i <= fNumberOfBins; i++) {
a38fd7f3 405 // this is to simulate |Delta eta| or |Delta phi|
406 if(fNumberOfBins/2 - fStartBin + 1 < i && i < fNumberOfBins/2 + fStartBin ) continue;
407
408 cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<TMath::Abs(gHistBalance->GetBinCenter(i))<<endl;
409
410 gSumXi += TMath::Abs(gHistBalance->GetBinCenter(i)); // this is to simulate |Delta eta| or |Delta phi|
411 gSumBi += gHistBalance->GetBinContent(i);
412 gSumBiXi += gHistBalance->GetBinContent(i)*TMath::Abs(gHistBalance->GetBinCenter(i));
413 gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2);
414 gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2);
415 gSumDeltaBi2 += TMath::Power(gHistBalance->GetBinError(i),2);
416 gSumXi2DeltaBi2 += TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2) * TMath::Power(gHistBalance->GetBinError(i),2);
417
418 deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
419 integral += fP2Step*gHistBalance->GetBinContent(i);
420 }
421 for(Int_t i = fStartBin; i < fNumberOfBins; i++)
422 deltaErrorNew += gHistBalance->GetBinError(i)*(TMath::Abs(gHistBalance->GetBinCenter(i))*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
423
424 Double_t integralError = TMath::Sqrt(deltaBalP2);
425
426 Double_t delta = gSumBiXi / gSumBi;
427 Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
428 cout<<"=================================================="<<endl;
429 cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
430 cout<<"New error: "<<deltaErrorNew<<endl;
431 cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
432 cout<<"=================================================="<<endl;
433 cout<<endl;
434}
eb63b883 435
648f1a5a 436//______________________________________________________//
437void drawBFPsi(const char* lhcPeriod = "LHC10h",
438 Int_t gCentrality = 1,
eb63b883 439 Int_t gDeltaEtaDeltaPhi = 2,
440 Double_t psiMin = -0.5, Double_t psiMax = 0.5,
441 Double_t ptTriggerMin = -1.,
442 Double_t ptTriggerMax = -1.,
443 Double_t ptAssociatedMin = -1.,
444 Double_t ptAssociatedMax = -1.) {
445 //Macro that draws the BF distributions for each centrality bin
446 //for reaction plane dependent analysis
447 //Author: Panos.Christakoglou@nikhef.nl
448 gROOT->LoadMacro("~/SetPlotStyle.C");
449 SetPlotStyle();
450
451 //Get the input file
648f1a5a 452 TString filename = lhcPeriod; filename +="/PttFrom";
db7174c0 453 filename += Form("%.1f",ptTriggerMin); filename += "To";
454 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
455 filename += Form("%.1f",ptAssociatedMin); filename += "To";
456 filename += Form("%.1f",ptAssociatedMax);
457 filename += "/balanceFunction.Centrality";
eb63b883 458 filename += gCentrality; filename += ".In";
459 if(gDeltaEtaDeltaPhi == 1) filename += "DeltaEta.Psi";
460 else if(gDeltaEtaDeltaPhi == 2) filename += "DeltaPhi.Psi";
461 if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
462 else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
463 else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
648f1a5a 464 else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
465 else filename += "All.PttFrom";
db7174c0 466 filename += Form("%.1f",ptTriggerMin); filename += "To";
467 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
468 filename += Form("%.1f",ptAssociatedMin); filename += "To";
469 filename += Form("%.1f",ptAssociatedMax); filename += ".root";
eb63b883 470
471 //Open the file
472 TFile *f = TFile::Open(filename.Data());
473 if((!f)||(!f->IsOpen())) {
474 Printf("The file %s is not found. Aborting...",filename);
475 return listBF;
476 }
477 //f->ls();
478
479 //Raw balance function
480 TH1D *gHistBalanceFunction = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunction"));
481 gHistBalanceFunction->SetStats(kFALSE);
482 gHistBalanceFunction->SetMarkerStyle(20);
483 gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
484 if(gDeltaEtaDeltaPhi == 2) {
485 gHistBalanceFunction->GetYaxis()->SetTitle("B(#Delta #varphi)");
486 gHistBalanceFunction->GetXaxis()->SetTitle("#Delta #varphi (deg.)");
487 }
488
489 //Shuffled balance function
490 TH1D *gHistBalanceFunctionShuffled = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionShuffled"));
491 gHistBalanceFunction->SetStats(kFALSE);
492 gHistBalanceFunctionShuffled->SetMarkerStyle(24);
493
494 //Mixed balance function
495 TH1D *gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionMixed"));
496 gHistBalanceFunction->SetStats(kFALSE);
497 gHistBalanceFunctionMixed->SetMarkerStyle(25);
498
499 //Subtracted balance function
500 TH1D *gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionSubtracted"));
501 gHistBalanceFunctionSubtracted->SetStats(kFALSE);
502 gHistBalanceFunctionSubtracted->SetMarkerStyle(20);
503 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
504 if(gDeltaEtaDeltaPhi == 2) {
505 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitle("B(#Delta #varphi)");
506 gHistBalanceFunctionSubtracted->GetXaxis()->SetTitle("#Delta #varphi (deg.)");
507 }
508
509 TString pngName;
510 TLegend *legend;
511
512 TString centralityLatex = "Centrality: ";
513 centralityLatex += centralityArray[gCentrality-1];
514 centralityLatex += "%";
515
516 TString psiLatex;
517 if((psiMin == -0.5)&&(psiMax == 0.5))
518 psiLatex = " -7.5^{o} < #varphi - #Psi_{2} < 7.5^{o}";
519 else if((psiMin == 0.5)&&(psiMax == 1.5))
520 psiLatex = " 37.5^{o} < #varphi - #Psi_{2} < 52.5^{o}";
521 else if((psiMin == 1.5)&&(psiMax == 2.5))
522 psiLatex = " 82.5^{o} < #varphi - #Psi_{2} < 97.5^{o}";
523 else
524 psiLatex = " 0^{o} < #varphi - #Psi_{2} < 180^{o}";
525
526 TString pttLatex = Form("%.1f",ptTriggerMin);
527 pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
528 pttLatex += " GeV/c";
529
530 TString ptaLatex = Form("%.1f",ptAssociatedMin);
531 ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
532 ptaLatex += " GeV/c";
533
534 //Draw the results
535 TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
536 c1->SetFillColor(10); c1->SetHighLightColor(10);
537 c1->SetLeftMargin(0.17); c1->SetTopMargin(0.05);
538 gHistBalanceFunction->SetTitle("");
539 gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.4);
540 gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
541 gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
542 gHistBalanceFunction->DrawCopy("E");
543 gHistBalanceFunctionShuffled->DrawCopy("ESAME");
544 gHistBalanceFunctionMixed->DrawCopy("ESAME");
545
546 legend = new TLegend(0.2,0.72,0.45,0.92,"","brNDC");
547 legend->SetTextSize(0.05);
548 legend->SetTextFont(42);
549 legend->SetBorderSize(0);
550 legend->SetFillStyle(0);
551 legend->SetFillColor(10);
552 legend->SetMargin(0.25);
553 legend->SetShadowColor(10);
554 legend->AddEntry(gHistBalanceFunction,"Data","lp");
555 legend->AddEntry(gHistBalanceFunctionShuffled,"Shuffled data","lp");
556 legend->AddEntry(gHistBalanceFunctionMixed,"Mixed data","lp");
557 legend->Draw();
558
559 TLatex *latexInfo1 = new TLatex();
560 latexInfo1->SetNDC();
561 latexInfo1->SetTextSize(0.04);
562 latexInfo1->SetTextColor(1);
563 latexInfo1->DrawLatex(0.58,0.88,centralityLatex.Data());
564 latexInfo1->DrawLatex(0.58,0.82,psiLatex.Data());
565 latexInfo1->DrawLatex(0.58,0.76,pttLatex.Data());
566 latexInfo1->DrawLatex(0.58,0.70,ptaLatex.Data());
567
568 //pngName = "BalanceFunctionDeltaEta.Centrality";
569 //pngName += centralityArray[gCentrality-1];
570 //pngName += ".Psi"; //pngName += psiMin; pngName += "To"; pngName += psiMax;
571 //pngName += ".png";
572 //c1->SaveAs(pngName.Data());
573
574 TCanvas *c2 = new TCanvas("c2","",600,0,600,500);
575 c2->SetFillColor(10); c2->SetHighLightColor(10);
576 c2->SetLeftMargin(0.17); c2->SetTopMargin(0.05);
577 gHistBalanceFunctionSubtracted->SetTitle("");
578 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.4);
579 gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
580 gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
581 gHistBalanceFunctionSubtracted->DrawCopy("E");
582
db7174c0 583 //Opening the output ascii files
584 TString filenameMean = "meanIn";
585 TString filenameSigma = "sigmaIn";
586 TString filenameSkewness = "skewnessIn";
587 TString filenameKurtosis = "kurtosisIn";
588 if(gDeltaEtaDeltaPhi == 1) {
589 filenameMean += "DeltaEta.Psi"; filenameSigma += "DeltaEta.Psi";
590 filenameSkewness += "DeltaEta.Psi"; filenameKurtosis += "DeltaEta.Psi";
591 }
592 else if(gDeltaEtaDeltaPhi == 2) {
593 filenameMean += "DeltaPhi.Psi"; filenameSigma += "DeltaPhi.Psi";
594 filenameSkewness += "DeltaPhi.Psi"; filenameKurtosis += "DeltaPhi.Psi";
595 }
596 if((psiMin == -0.5)&&(psiMax == 0.5)) {
597 filenameMean += "InPlane.Ptt"; filenameSigma += "InPlane.Ptt";
598 filenameSkewness += "InPlane.Ptt"; filenameKurtosis += "InPlane.Ptt";
599 }
600 else if((psiMin == 0.5)&&(psiMax == 1.5)) {
601 filenameMean += "Intermediate.Ptt"; filenameSigma += "Intermediate.Ptt";
602 filenameSkewness += "Intermediate.Ptt";
603 filenameKurtosis += "Intermediate.Ptt";
604 }
605 else if((psiMin == 1.5)&&(psiMax == 2.5)) {
606 filenameMean += "OutOfPlane.Ptt"; filenameSigma += "OutOfPlane.Ptt";
607 filenameSkewness += "OutOfPlane.Ptt";
608 filenameKurtosis += "OutOfPlane.Ptt";
609 }
610 else if((psiMin == 2.5)&&(psiMax == 3.5)) {
611 filenameMean += "Rest.Ptt"; filenameSigma += "Rest.Ptt";
612 filenameSkewness += "Rest.Ptt"; filenameKurtosis += "Rest.Ptt";
613 }
614 else {
615 filenameMean += "All.Ptt"; filenameSigma += "All.Ptt";
616 filenameSkewness += "All.Ptt"; filenameKurtosis += "All.Ptt";
617 }
618 filenameMean += Form("%.1f",ptTriggerMin); filenameMean += "To";
619 filenameMean += Form("%.1f",ptTriggerMax); filenameMean += "PtaFrom";
620 filenameMean += Form("%.1f",ptAssociatedMin); filenameMean += "To";
621 filenameMean += Form("%.1f",ptAssociatedMax); filenameMean += ".txt";
622 filenameSigma += Form("%.1f",ptTriggerMin); filenameSigma += "To";
623 filenameSigma += Form("%.1f",ptTriggerMax); filenameSigma += "PtaFrom";
624 filenameSigma += Form("%.1f",ptAssociatedMin); filenameSigma += "To";
625 filenameSigma += Form("%.1f",ptAssociatedMax); filenameSigma += ".txt";
626 filenameSkewness += Form("%.1f",ptTriggerMin); filenameSkewness += "To";
627 filenameSkewness += Form("%.1f",ptTriggerMax); filenameSkewness += "PtaFrom";
628 filenameSkewness += Form("%.1f",ptAssociatedMin); filenameSkewness += "To";
629 filenameSkewness += Form("%.1f",ptAssociatedMax);
630 filenameSkewness += ".txt";
631 filenameKurtosis += Form("%.1f",ptTriggerMin); filenameKurtosis += "To";
632 filenameKurtosis += Form("%.1f",ptTriggerMax); filenameKurtosis += "PtaFrom";
633 filenameKurtosis += Form("%.1f",ptAssociatedMin); filenameKurtosis += "To";
634 filenameKurtosis += Form("%.1f",ptAssociatedMax);
635 filenameKurtosis += ".txt";
636
eb63b883 637 TString meanLatex, rmsLatex, skewnessLatex, kurtosisLatex;
638 meanLatex = "#mu = ";
639 meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMean());
640 meanLatex += " #pm ";
641 meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMeanError());
db7174c0 642 ofstream fileMean(filenameMean.Data(),ios::app);
643 fileMean << gCentrality << " " << gHistBalanceFunctionSubtracted->GetMean() << " " <<gHistBalanceFunctionSubtracted->GetMeanError()<<endl;
644 fileMean.close();
645
eb63b883 646 rmsLatex = "#sigma = ";
647 rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMS());
648 rmsLatex += " #pm ";
649 rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMSError());
db7174c0 650 ofstream fileSigma(filenameSigma.Data(),ios::app);
651 fileSigma << gCentrality << " " << gHistBalanceFunctionSubtracted->GetRMS() << " " <<gHistBalanceFunctionSubtracted->GetRMSError()<<endl;
652 fileSigma.close();
eb63b883 653
654 skewnessLatex = "S = ";
655 skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(1));
656 skewnessLatex += " #pm ";
657 skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(11));
db7174c0 658 ofstream fileSkewness(filenameSkewness.Data(),ios::app);
659 fileSkewness << gCentrality << " " << gHistBalanceFunctionSubtracted->GetSkewness(1) << " " <<gHistBalanceFunctionSubtracted->GetSkewness(11)<<endl;
660 fileSkewness.close();
661
eb63b883 662 kurtosisLatex = "K = ";
663 kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(1));
664 kurtosisLatex += " #pm ";
665 kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(11));
db7174c0 666 ofstream fileKurtosis(filenameKurtosis.Data(),ios::app);
667 fileKurtosis << gCentrality << " " << gHistBalanceFunctionSubtracted->GetKurtosis(1) << " " <<gHistBalanceFunctionSubtracted->GetKurtosis(11)<<endl;
668 fileKurtosis.close();
669
eb63b883 670 Printf("Mean: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetMean(),gHistBalanceFunctionSubtracted->GetMeanError());
671 Printf("RMS: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetRMS(),gHistBalanceFunctionSubtracted->GetRMSError());
672 Printf("Skeweness: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetSkewness(1),gHistBalanceFunctionSubtracted->GetSkewness(11));
673 Printf("Kurtosis: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetKurtosis(1),gHistBalanceFunctionSubtracted->GetKurtosis(11));
674
675 latexInfo1->DrawLatex(0.18,0.88,centralityLatex.Data());
676 latexInfo1->DrawLatex(0.18,0.82,psiLatex.Data());
677 latexInfo1->DrawLatex(0.18,0.76,pttLatex.Data());
678 latexInfo1->DrawLatex(0.18,0.70,ptaLatex.Data());
679
680 TLatex *latexResults = new TLatex();
681 latexResults->SetNDC();
682 latexResults->SetTextSize(0.04);
683 latexResults->SetTextColor(1);
684 latexResults->DrawLatex(0.6,0.88,meanLatex.Data());
685 latexResults->DrawLatex(0.6,0.82,rmsLatex.Data());
686 latexResults->DrawLatex(0.6,0.76,skewnessLatex.Data());
687 latexResults->DrawLatex(0.6,0.70,kurtosisLatex.Data());
688
689}