1 const Int_t numberOfCentralityBins = 12;
2 TString centralityArray[numberOfCentralityBins] = {"0-100","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-100","0-1","1-2","2-3"};
5 const Int_t gRebin = 1;
6 void drawBalanceFunction2DPsi(const char* filename = "AnalysisResultsPsi.root",
9 const char* gCentralityEstimator = 0x0,
10 Bool_t kShowShuffled = kFALSE,
11 Bool_t kShowMixed = kTRUE,
12 Double_t psiMin = -0.5, Double_t psiMax = 0.5,
13 Double_t vertexZMin = -10.,
14 Double_t vertexZMax = 10.,
15 Double_t ptTriggerMin = -1.,
16 Double_t ptTriggerMax = -1.,
17 Double_t ptAssociatedMin = -1.,
18 Double_t ptAssociatedMax = -1.,
19 Bool_t kUseVzBinning = kTRUE,
20 Bool_t k2pMethod = kTRUE,
21 TString eventClass = "EventPlane") //Can be "EventPlane", "Centrality", "Multiplicity"
23 //Macro that draws the BF distributions for each centrality bin
24 //for reaction plane dependent analysis
25 //Author: Panos.Christakoglou@nikhef.nl
26 //Load the PWG2ebye library
27 gSystem->Load("libANALYSIS.so");
28 gSystem->Load("libANALYSISalice.so");
29 gSystem->Load("libEventMixing.so");
30 gSystem->Load("libCORRFW.so");
31 gSystem->Load("libPWGTools.so");
32 gSystem->Load("libPWGCFebye.so");
34 //gROOT->LoadMacro("~/SetPlotStyle.C");
36 gStyle->SetPalette(1,0);
38 //Prepare the objects and return them
39 TList *listBF = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,0);
40 TList *listBFShuffled = NULL;
41 if(kShowShuffled) listBFShuffled = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,1);
42 TList *listBFMixed = NULL;
43 if(kShowMixed) listBFMixed = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,2);
45 Printf("The TList object was not created");
49 draw(listBF,listBFShuffled,listBFMixed,gCentrality,gCentralityEstimator,
50 psiMin,psiMax,vertexZMin,vertexZMax,
51 ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,
52 kUseVzBinning,k2pMethod,eventClass);
55 //______________________________________________________//
56 TList *GetListOfObjects(const char* filename,
59 const char *gCentralityEstimator,
61 //Get the TList objects (QA, bf, bf shuffled)
65 TFile *f = TFile::Open(filename,"UPDATE");
66 if((!f)||(!f->IsOpen())) {
67 Printf("The file %s is not found. Aborting...",filename);
72 TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
74 Printf("The TDirectoryFile is not found. Aborting...",filename);
81 //cout<<"no shuffling - no mixing"<<endl;
82 listBFName = "listBFPsi_";
85 //cout<<"shuffling - no mixing"<<endl;
86 listBFName = "listBFPsiShuffled_";
89 //cout<<"no shuffling - mixing"<<endl;
90 listBFName = "listBFPsiMixed_";
92 listBFName += centralityArray[gCentrality-1];
94 listBFName += "_Bit"; listBFName += gBit; }
95 if(gCentralityEstimator) {
96 listBFName += "_"; listBFName += gCentralityEstimator;}
98 // histograms were already retrieved (in first iteration)
99 if(dir->Get(Form("%s_histograms",listBFName.Data()))){
100 //cout<<"second iteration"<<endl;
101 listBF = dynamic_cast<TList *>(dir->Get(Form("%s_histograms",listBFName.Data())));
104 // histograms were not yet retrieved (this is the first iteration)
107 listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
108 cout<<"======================================================="<<endl;
109 cout<<"List name: "<<listBF->GetName()<<endl;
115 histoName = "fHistP";
117 histoName = "fHistP_shuffle";
119 histoName = "fHistP";
120 if(gCentralityEstimator) histoName += gCentralityEstimator;
121 AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
123 Printf("fHistP %s not found!!!",histoName.Data());
126 fHistP->FillParent(); fHistP->DeleteContainers();
129 histoName = "fHistN";
131 histoName = "fHistN_shuffle";
133 histoName = "fHistN";
134 if(gCentralityEstimator) histoName += gCentralityEstimator;
135 AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
137 Printf("fHistN %s not found!!!",histoName.Data());
140 fHistN->FillParent(); fHistN->DeleteContainers();
143 histoName = "fHistPN";
145 histoName = "fHistPN_shuffle";
147 histoName = "fHistPN";
148 if(gCentralityEstimator) histoName += gCentralityEstimator;
149 AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
151 Printf("fHistPN %s not found!!!",histoName.Data());
154 fHistPN->FillParent(); fHistPN->DeleteContainers();
157 histoName = "fHistNP";
159 histoName = "fHistNP_shuffle";
161 histoName = "fHistNP";
162 if(gCentralityEstimator) histoName += gCentralityEstimator;
163 AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
165 Printf("fHistNP %s not found!!!",histoName.Data());
168 fHistNP->FillParent(); fHistNP->DeleteContainers();
171 histoName = "fHistPP";
173 histoName = "fHistPP_shuffle";
175 histoName = "fHistPP";
176 if(gCentralityEstimator) histoName += gCentralityEstimator;
177 AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
179 Printf("fHistPP %s not found!!!",histoName.Data());
182 fHistPP->FillParent(); fHistPP->DeleteContainers();
185 histoName = "fHistNN";
187 histoName = "fHistNN_shuffle";
189 histoName = "fHistNN";
190 if(gCentralityEstimator) histoName += gCentralityEstimator;
191 AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
193 Printf("fHistNN %s not found!!!",histoName.Data());
196 fHistNN->FillParent(); fHistNN->DeleteContainers();
199 listBF->Write(Form("%s_histograms",listBFName.Data()), TObject::kSingleKey);
208 //______________________________________________________//
209 void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
210 Int_t gCentrality, const char* gCentralityEstimator,
211 Double_t psiMin, Double_t psiMax,
214 Double_t ptTriggerMin, Double_t ptTriggerMax,
215 Double_t ptAssociatedMin, Double_t ptAssociatedMax,
216 Bool_t kUseVzBinning=kFALSE,
217 Bool_t k2pMethod = kFALSE, TString eventClass) {
226 //Printf("=================");
227 TString histoName = "fHistP";
228 if(gCentralityEstimator) histoName += gCentralityEstimator;
229 hP = (AliTHn*) listBF->FindObject(histoName.Data());
230 hP->SetName("gHistP");
231 histoName = "fHistN";
232 if(gCentralityEstimator) histoName += gCentralityEstimator;
233 hN = (AliTHn*) listBF->FindObject(histoName.Data());
234 hN->SetName("gHistN");
235 histoName = "fHistPN";
236 if(gCentralityEstimator) histoName += gCentralityEstimator;
237 hPN = (AliTHn*) listBF->FindObject(histoName.Data());
238 hPN->SetName("gHistPN");
239 histoName = "fHistNP";
240 if(gCentralityEstimator) histoName += gCentralityEstimator;
241 hNP = (AliTHn*) listBF->FindObject(histoName.Data());
242 hNP->SetName("gHistNP");
243 histoName = "fHistPP";
244 if(gCentralityEstimator) histoName += gCentralityEstimator;
245 hPP = (AliTHn*) listBF->FindObject(histoName.Data());
246 hPP->SetName("gHistPP");
247 histoName = "fHistNN";
248 if(gCentralityEstimator) histoName += gCentralityEstimator;
249 hNN = (AliTHn*) listBF->FindObject(histoName.Data());
250 hNN->SetName("gHistNN");
252 AliBalancePsi *b = new AliBalancePsi();
253 b->SetEventClass(eventClass);
260 if(kUseVzBinning) b->SetVertexZBinning(kTRUE);
263 //balance function shuffling
264 AliTHn *hPShuffled = NULL;
265 AliTHn *hNShuffled = NULL;
266 AliTHn *hPNShuffled = NULL;
267 AliTHn *hNPShuffled = NULL;
268 AliTHn *hPPShuffled = NULL;
269 AliTHn *hNNShuffled = NULL;
271 //listBFShuffled->ls();
272 histoName = "fHistP_shuffle";
273 if(gCentralityEstimator) histoName += gCentralityEstimator;
274 hPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
275 hPShuffled->SetName("gHistPShuffled");
276 histoName = "fHistN_shuffle";
277 if(gCentralityEstimator) histoName += gCentralityEstimator;
278 hNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
279 hNShuffled->SetName("gHistNShuffled");
280 histoName = "fHistPN_shuffle";
281 if(gCentralityEstimator) histoName += gCentralityEstimator;
282 hPNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
283 hPNShuffled->SetName("gHistPNShuffled");
284 histoName = "fHistNP_shuffle";
285 if(gCentralityEstimator) histoName += gCentralityEstimator;
286 hNPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
287 hNPShuffled->SetName("gHistNPShuffled");
288 histoName = "fHistPP_shuffle";
289 if(gCentralityEstimator) histoName += gCentralityEstimator;
290 hPPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
291 hPPShuffled->SetName("gHistPPShuffled");
292 histoName = "fHistNN_shuffle";
293 if(gCentralityEstimator) histoName += gCentralityEstimator;
294 hNNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
295 hNNShuffled->SetName("gHistNNShuffled");
297 AliBalancePsi *bShuffled = new AliBalancePsi();
298 bShuffled->SetEventClass(eventClass);
299 bShuffled->SetHistNp(hPShuffled);
300 bShuffled->SetHistNn(hNShuffled);
301 bShuffled->SetHistNpn(hPNShuffled);
302 bShuffled->SetHistNnp(hNPShuffled);
303 bShuffled->SetHistNpp(hPPShuffled);
304 bShuffled->SetHistNnn(hNNShuffled);
305 if(kUseVzBinning) bShuffled->SetVertexZBinning(kTRUE);
309 //balance function mixing
310 AliTHn *hPMixed = NULL;
311 AliTHn *hNMixed = NULL;
312 AliTHn *hPNMixed = NULL;
313 AliTHn *hNPMixed = NULL;
314 AliTHn *hPPMixed = NULL;
315 AliTHn *hNNMixed = NULL;
319 histoName = "fHistP";
320 if(gCentralityEstimator) histoName += gCentralityEstimator;
321 hPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
322 hPMixed->SetName("gHistPMixed");
323 histoName = "fHistN";
324 if(gCentralityEstimator) histoName += gCentralityEstimator;
325 hNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
326 hNMixed->SetName("gHistNMixed");
327 histoName = "fHistPN";
328 if(gCentralityEstimator) histoName += gCentralityEstimator;
329 hPNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
330 hPNMixed->SetName("gHistPNMixed");
331 histoName = "fHistNP";
332 if(gCentralityEstimator) histoName += gCentralityEstimator;
333 hNPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
334 histoName = "fHistNP";
335 if(gCentralityEstimator) histoName += gCentralityEstimator;
336 hNPMixed->SetName("gHistNPMixed");
337 histoName = "fHistPP";
338 if(gCentralityEstimator) histoName += gCentralityEstimator;
339 hPPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
340 hPPMixed->SetName("gHistPPMixed");
341 histoName = "fHistNN";
342 if(gCentralityEstimator) histoName += gCentralityEstimator;
343 hNNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
344 hNNMixed->SetName("gHistNNMixed");
346 AliBalancePsi *bMixed = new AliBalancePsi();
347 bMixed->SetEventClass(eventClass);
348 bMixed->SetHistNp(hPMixed);
349 bMixed->SetHistNn(hNMixed);
350 bMixed->SetHistNpn(hPNMixed);
351 bMixed->SetHistNnp(hNPMixed);
352 bMixed->SetHistNpp(hPPMixed);
353 bMixed->SetHistNnn(hNNMixed);
354 if(kUseVzBinning) bMixed->SetVertexZBinning(kTRUE);
358 TH2D *gHistBalanceFunction;
359 TH2D *gHistBalanceFunctionSubtracted;
360 TH2D *gHistBalanceFunctionShuffled;
361 TH2D *gHistBalanceFunctionMixed;
362 TString histoTitle, pngName;
364 if(eventClass == "Centrality"){
365 histoTitle = "Centrality: ";
366 histoTitle += psiMin;
368 histoTitle += psiMax;
370 histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
372 else if(eventClass == "Multiplicity"){
373 histoTitle = "Multiplicity: ";
374 histoTitle += psiMin;
376 histoTitle += psiMax;
377 histoTitle += " tracks";
378 histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
380 else{ // "EventPlane" (default)
381 histoTitle = "Centrality: ";
382 histoTitle += centralityArray[gCentrality-1];
384 if((psiMin == -0.5)&&(psiMax == 0.5))
385 histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
386 else if((psiMin == 0.5)&&(psiMax == 1.5))
387 histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
388 else if((psiMin == 1.5)&&(psiMax == 2.5))
389 histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
391 histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
396 gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
398 cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
402 gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
403 gHistBalanceFunction->SetTitle(histoTitle.Data());
404 gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
405 gHistBalanceFunction->SetName("gHistBalanceFunction");
411 gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
413 cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
417 gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
418 gHistBalanceFunctionShuffled->SetTitle(histoTitle.Data());
419 gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
420 gHistBalanceFunctionShuffled->SetName("gHistBalanceFunctionShuffled");
426 gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
428 cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
432 gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
433 gHistBalanceFunctionMixed->SetTitle(histoTitle.Data());
434 gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
435 gHistBalanceFunctionMixed->SetName("gHistBalanceFunctionMixed");
437 gHistBalanceFunctionSubtracted = dynamic_cast<TH2D *>(gHistBalanceFunction->Clone());
438 gHistBalanceFunctionSubtracted->Add(gHistBalanceFunctionMixed,-1);
439 gHistBalanceFunctionSubtracted->SetTitle(histoTitle.Data());
440 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
441 gHistBalanceFunctionSubtracted->SetName("gHistBalanceFunctionSubtracted");
445 TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
446 c1->SetFillColor(10);
447 c1->SetHighLightColor(10);
448 c1->SetLeftMargin(0.15);
449 gHistBalanceFunction->DrawCopy("lego2");
450 gPad->SetTheta(30); // default is 30
451 gPad->SetPhi(-60); // default is 30
453 TCanvas *c1a = new TCanvas("c1a","",600,0,600,500);
454 c1a->SetFillColor(10);
455 c1a->SetHighLightColor(10);
456 c1a->SetLeftMargin(0.15);
457 gHistBalanceFunction->DrawCopy("colz");
460 TCanvas *c2 = new TCanvas("c2","",100,100,600,500);
461 c2->SetFillColor(10);
462 c2->SetHighLightColor(10);
463 c2->SetLeftMargin(0.15);
464 gHistBalanceFunctionShuffled->DrawCopy("lego2");
465 gPad->SetTheta(30); // default is 30
466 gPad->SetPhi(-60); // default is 30
468 TCanvas *c2a = new TCanvas("c2a","",700,100,600,500);
469 c2a->SetFillColor(10);
470 c2a->SetHighLightColor(10);
471 c2a->SetLeftMargin(0.15);
472 gHistBalanceFunctionShuffled->DrawCopy("colz");
476 TCanvas *c3 = new TCanvas("c3","",200,200,600,500);
477 c3->SetFillColor(10);
478 c3->SetHighLightColor(10);
479 c3->SetLeftMargin(0.15);
480 gHistBalanceFunctionMixed->DrawCopy("lego2");
481 gPad->SetTheta(30); // default is 30
482 gPad->SetPhi(-60); // default is 30
484 TCanvas *c3a = new TCanvas("c3a","",800,200,600,500);
485 c3a->SetFillColor(10);
486 c3a->SetHighLightColor(10);
487 c3a->SetLeftMargin(0.15);
488 gHistBalanceFunctionMixed->DrawCopy("colz");
490 TCanvas *c4 = new TCanvas("c4","",300,300,600,500);
491 c4->SetFillColor(10);
492 c4->SetHighLightColor(10);
493 c4->SetLeftMargin(0.15);
494 gHistBalanceFunctionSubtracted->DrawCopy("lego2");
495 gPad->SetTheta(30); // default is 30
496 gPad->SetPhi(-60); // default is 30
498 TCanvas *c4a = new TCanvas("c4a","",900,300,600,500);
499 c4a->SetFillColor(10);
500 c4a->SetHighLightColor(10);
501 c4a->SetLeftMargin(0.15);
502 gHistBalanceFunctionSubtracted->DrawCopy("colz");
504 fitbalanceFunction(gCentrality, psiMin , psiMax,
505 ptTriggerMin, ptTriggerMax,
506 ptAssociatedMin, ptAssociatedMax,
507 gHistBalanceFunctionSubtracted,k2pMethod, eventClass);
510 TString newFileName = "balanceFunction2D.";
511 if(eventClass == "Centrality"){
512 newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
513 newFileName += ".PsiAll.PttFrom";
515 else if(eventClass == "Multiplicity"){
516 newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
517 newFileName += ".PsiAll.PttFrom";
519 else{ // "EventPlane" (default)
520 newFileName += "Centrality";
521 newFileName += gCentrality; newFileName += ".Psi";
522 if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
523 else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
524 else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
525 else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
526 else newFileName += "All.PttFrom";
528 newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
529 newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
530 newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
531 newFileName += Form("%.1f",ptAssociatedMax);
532 if(k2pMethod) newFileName += "_2pMethod";
534 // newFileName += "_";
535 // newFileName += Form("%.1f",psiMin);
536 // newFileName += "-";
537 // newFileName += Form("%.1f",psiMax);
538 // newFileName += ".root";
540 TFile *fOutput = new TFile(newFileName.Data(),"recreate");
542 /*hP->Write(); hN->Write();
543 hPN->Write(); hNP->Write();
544 hPP->Write(); hNN->Write();
545 hPShuffled->Write(); hNShuffled->Write();
546 hPNShuffled->Write(); hNPShuffled->Write();
547 hPPShuffled->Write(); hNNShuffled->Write();
548 hPMixed->Write(); hNMixed->Write();
549 hPNMixed->Write(); hNPMixed->Write();
550 hPPMixed->Write(); hNNMixed->Write();*/
551 gHistBalanceFunction->Write();
552 if(listBFShuffled) gHistBalanceFunctionShuffled->Write();
554 gHistBalanceFunctionMixed->Write();
555 gHistBalanceFunctionSubtracted->Write();
560 //____________________________________________________________//
561 void fitbalanceFunction(Int_t gCentrality = 1,
562 Double_t psiMin = -0.5, Double_t psiMax = 3.5,
563 Double_t ptTriggerMin = -1.,
564 Double_t ptTriggerMax = -1.,
565 Double_t ptAssociatedMin = -1.,
566 Double_t ptAssociatedMax = -1.,
568 Bool_t k2pMethod = kFALSE,
569 TString eventClass="EventPlane") {
570 //balancing charges: [1]*TMath::Exp(-0.5*TMath::Power(((x - [3])/[2]),2)-0.5*TMath::Power(((y - [5])/[4]),2))
571 //short range correlations: [6]*TMath::Exp(-0.5*TMath::Power(((x - [8])/[7]),2)-0.5*TMath::Power(((y - [10])/[9]),2))
572 cout<<"FITTING FUNCTION"<<endl;
574 TF2 *gFitFunction = new TF2("gFitFunction","[0] + [1]*TMath::Exp(-0.5*TMath::Power(((x - [3])/[2]),2)-0.5*TMath::Power(((y - [5])/[4]),2)) + [6]*TMath::Exp(-0.5*TMath::Power(((x - [8])/[7]),2)-0.5*TMath::Power(((y - [10])/[9]),2))",-1.2,1.2,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
575 gFitFunction->SetName("gFitFunction");
578 gFitFunction->SetParName(0,"N1");
579 gFitFunction->SetParameter(0,1.0);
581 //2D balance function
582 gFitFunction->SetParName(1,"N_{BF}");
583 gFitFunction->SetParameter(1,1.0);
584 gFitFunction->SetParLimits(1, 0., 100.);
585 gFitFunction->SetParName(2,"Sigma_{BF}(delta eta)");
586 gFitFunction->SetParameter(2,0.6);
587 gFitFunction->SetParLimits(2, 0., 1.);
588 gFitFunction->SetParName(3,"Mean_{BF}(delta eta)");
589 gFitFunction->SetParameter(3,0.0);
590 gFitFunction->SetParLimits(3, -0.2, 0.2);
591 gFitFunction->SetParName(4,"Sigma_{BF}(delta phi)");
592 gFitFunction->SetParameter(4,0.6);
593 gFitFunction->SetParLimits(4, 0., 1.);
594 gFitFunction->SetParName(5,"Mean_{BF}(delta phi)");
595 gFitFunction->SetParameter(5,0.0);
596 gFitFunction->SetParLimits(5, -0.2, 0.2);
598 //short range structure
599 gFitFunction->SetParName(6,"N_{SR}");
600 gFitFunction->SetParameter(6,5.0);
601 gFitFunction->SetParLimits(6, 0., 100.);
602 gFitFunction->SetParName(7,"Sigma_{SR}(delta eta)");
603 gFitFunction->SetParameter(7,0.01);
604 gFitFunction->SetParLimits(7, 0.0, 0.1);
605 gFitFunction->SetParName(8,"Mean_{SR}(delta eta)");
606 gFitFunction->SetParameter(8,0.0);
607 gFitFunction->SetParLimits(8, -0.01, 0.01);
608 gFitFunction->SetParName(9,"Sigma_{SR}(delta phi)");
609 gFitFunction->SetParameter(9,0.01);
610 gFitFunction->SetParLimits(9, 0.0, 0.1);
611 gFitFunction->SetParName(10,"Mean_{SR}(delta phi)");
612 gFitFunction->SetParameter(10,0.0);
613 gFitFunction->SetParLimits(10, -0.01, 0.01);
616 //Cloning the histogram
617 TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
618 gHistResidual->SetName("gHistResidual");
619 gHistResidual->Sumw2();
622 for(Int_t iAttempt = 0; iAttempt < 10; iAttempt++) {
623 gHist->Fit("gFitFunction","nm");
624 for(Int_t iParam = 0; iParam < 11; iParam++)
625 gFitFunction->SetParameter(iParam,gFitFunction->GetParameter(iParam));
627 cout<<"======================================================"<<endl;
628 cout<<"Fit chi2/ndf: "<<gFitFunction->GetChisquare()/gFitFunction->GetNDF()<<" - chi2: "<<gFitFunction->GetChisquare()<<" - ndf: "<<gFitFunction->GetNDF()<<endl;
629 cout<<"======================================================"<<endl;
631 //Getting the residual
632 gHistResidual->Add(gFitFunction,-1);
634 //Write to output file
635 TString newFileName = "balanceFunctionFit2D.";
636 if(eventClass == "Centrality"){
637 newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
638 newFileName += ".PsiAll.PttFrom";
640 else if(eventClass == "Multiplicity"){
641 newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
642 newFileName += ".PsiAll.PttFrom";
644 else{ // "EventPlane" (default)
645 newFileName += "Centrality";
646 newFileName += gCentrality; newFileName += ".Psi";
647 if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
648 else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
649 else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
650 else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
651 else newFileName += "All.PttFrom";
653 newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
654 newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
655 newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
656 newFileName += Form("%.1f",ptAssociatedMax);
657 if(k2pMethod) newFileName += "_2pMethod";
660 newFileName += Form("%.1f",psiMin);
662 newFileName += Form("%.1f",psiMax);
663 newFileName += ".root";
665 TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
667 gHistResidual->Write();
668 gFitFunction->Write();
672 //____________________________________________________________//
673 void drawBFPsi2D(const char* lhcPeriod = "LHC11h",
674 const char* gCentralityEstimator = "V0M",
676 const char* gEventPlaneEstimator = "VZERO",
677 Int_t gCentrality = 1,
678 Bool_t kShowShuffled = kFALSE,
679 Bool_t kShowMixed = kFALSE,
680 Double_t psiMin = -0.5, Double_t psiMax = 0.5,
681 Double_t ptTriggerMin = -1.,
682 Double_t ptTriggerMax = -1.,
683 Double_t ptAssociatedMin = -1.,
684 Double_t ptAssociatedMax = -1.,
685 Bool_t k2pMethod = kTRUE) {
686 //Macro that draws the BF distributions for each centrality bin
687 //for reaction plane dependent analysis
688 //Author: Panos.Christakoglou@nikhef.nl
689 TGaxis::SetMaxDigits(3);
692 TString filename = lhcPeriod;
693 filename += "/Centrality"; filename += gCentralityEstimator;
694 filename += "_Bit"; filename += gBit;
695 filename += "_"; filename += gEventPlaneEstimator;
696 filename +="/PttFrom";
697 filename += Form("%.1f",ptTriggerMin); filename += "To";
698 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
699 filename += Form("%.1f",ptAssociatedMin); filename += "To";
700 filename += Form("%.1f",ptAssociatedMax);
701 filename += "/balanceFunction2D.Centrality";
702 filename += gCentrality; filename += ".Psi";
703 if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
704 else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
705 else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
706 else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
707 else filename += "All.PttFrom";
708 filename += Form("%.1f",ptTriggerMin); filename += "To";
709 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
710 filename += Form("%.1f",ptAssociatedMin); filename += "To";
711 filename += Form("%.1f",ptAssociatedMax);
712 if(k2pMethod) filename += "_2pMethod";
715 filename += Form("%.1f",psiMin);
717 filename += Form("%.1f",psiMax);
721 TFile *f = TFile::Open(filename.Data());
722 if((!f)||(!f->IsOpen())) {
723 Printf("The file %s is not found. Aborting...",filename);
728 //Raw balance function
729 TH1D *gHistBalanceFunction = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunction"));
730 gHistBalanceFunction->SetStats(kFALSE);
731 gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
732 gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
733 gHistBalanceFunction->GetZaxis()->SetNdivisions(10);
734 gHistBalanceFunction->GetXaxis()->SetTitleOffset(1.3);
735 gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
736 gHistBalanceFunction->GetZaxis()->SetTitleOffset(1.3);
737 gHistBalanceFunction->GetXaxis()->SetTitle("#Delta #eta");
738 gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (rad)");
739 gHistBalanceFunction->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
741 //Shuffled balance function
743 TH1D *gHistBalanceFunctionShuffled = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionShuffled"));
744 gHistBalanceFunctionShuffled->SetStats(kFALSE);
745 gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
746 gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
747 gHistBalanceFunctionShuffled->GetZaxis()->SetNdivisions(10);
748 gHistBalanceFunctionShuffled->GetXaxis()->SetTitleOffset(1.3);
749 gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
750 gHistBalanceFunctionShuffled->GetZaxis()->SetTitleOffset(1.3);
751 gHistBalanceFunctionShuffled->GetXaxis()->SetTitle("#Delta #eta");
752 gHistBalanceFunctionShuffled->GetYaxis()->SetTitle("#Delta #varphi (rad)");
753 gHistBalanceFunctionShuffled->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
756 //Mixed balance function
758 TH1D *gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionMixed"));
759 gHistBalanceFunctionMixed->SetStats(kFALSE);
760 gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
761 gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
762 gHistBalanceFunctionMixed->GetZaxis()->SetNdivisions(10);
763 gHistBalanceFunctionMixed->GetXaxis()->SetTitleOffset(1.3);
764 gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
765 gHistBalanceFunctionMixed->GetZaxis()->SetTitleOffset(1.3);
766 gHistBalanceFunctionMixed->GetXaxis()->SetTitle("#Delta #eta");
767 gHistBalanceFunctionMixed->GetYaxis()->SetTitle("#Delta #varphi (rad)");
768 gHistBalanceFunctionMixed->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
771 //Subtracted balance function
773 TH1D *gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionSubtracted"));
774 gHistBalanceFunctionSubtracted->SetStats(kFALSE);
775 gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
776 gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
777 gHistBalanceFunctionSubtracted->GetZaxis()->SetNdivisions(10);
778 gHistBalanceFunctionSubtracted->GetXaxis()->SetTitleOffset(1.3);
779 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
780 gHistBalanceFunctionSubtracted->GetZaxis()->SetTitleOffset(1.3);
781 gHistBalanceFunctionSubtracted->GetXaxis()->SetTitle("#Delta #eta");
782 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitle("#Delta #varphi (rad)");
783 gHistBalanceFunctionSubtracted->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
788 TString centralityLatex = "Centrality: ";
789 centralityLatex += centralityArray[gCentrality-1];
790 centralityLatex += "%";
793 if((psiMin == -0.5)&&(psiMax == 0.5))
794 psiLatex = " -7.5^{o} < #varphi - #Psi_{2} < 7.5^{o}";
795 else if((psiMin == 0.5)&&(psiMax == 1.5))
796 psiLatex = " 37.5^{o} < #varphi - #Psi_{2} < 52.5^{o}";
797 else if((psiMin == 1.5)&&(psiMax == 2.5))
798 psiLatex = " 82.5^{o} < #varphi - #Psi_{2} < 97.5^{o}";
800 psiLatex = " 0^{o} < #varphi - #Psi_{2} < 180^{o}";
802 TString pttLatex = Form("%.1f",ptTriggerMin);
803 pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
804 pttLatex += " GeV/c";
806 TString ptaLatex = Form("%.1f",ptAssociatedMin);
807 ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
808 ptaLatex += " GeV/c";
810 TLatex *latexInfo1 = new TLatex();
811 latexInfo1->SetNDC();
812 latexInfo1->SetTextSize(0.045);
813 latexInfo1->SetTextColor(1);
816 TCanvas *c1 = new TCanvas("c1","Raw balance function 2D",0,0,600,500);
817 c1->SetFillColor(10); c1->SetHighLightColor(10);
818 c1->SetLeftMargin(0.17); c1->SetTopMargin(0.05);
819 gHistBalanceFunction->SetTitle("");
820 gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.4);
821 gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
822 gHistBalanceFunction->GetXaxis()->SetRangeUser(-1.4,1.4);
823 gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
824 gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (rad)");
825 gHistBalanceFunction->DrawCopy("lego2");
826 gPad->SetTheta(30); // default is 30
827 gPad->SetPhi(-60); // default is 30
830 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
831 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
832 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
833 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
835 TString pngName = "BalanceFunction2D.";
836 pngName += "Centrality";
837 pngName += gCentrality;
838 if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
839 else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
840 else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
841 else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.PttFrom";
842 else pngName += "All.PttFrom";
843 pngName += Form("%.1f",ptTriggerMin); pngName += "To";
844 pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
845 pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
846 pngName += Form("%.1f",ptAssociatedMax);
847 if(k2pMethod) pngName += "_2pMethod";
850 c1->SaveAs(pngName.Data());
853 TCanvas *c2 = new TCanvas("c2","Shuffled balance function 2D",100,100,600,500);
854 c2->SetFillColor(10); c2->SetHighLightColor(10);
855 c2->SetLeftMargin(0.17); c2->SetTopMargin(0.05);
856 gHistBalanceFunctionShuffled->SetTitle("Shuffled events");
857 gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.4);
858 gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
859 gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
860 gHistBalanceFunctionShuffled->DrawCopy("lego2");
861 gPad->SetTheta(30); // default is 30
862 gPad->SetPhi(-60); // default is 30
865 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
866 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
867 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
868 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
872 TCanvas *c3 = new TCanvas("c3","Mixed balance function 2D",200,200,600,500);
873 c3->SetFillColor(10); c3->SetHighLightColor(10);
874 c3->SetLeftMargin(0.17); c3->SetTopMargin(0.05);
875 gHistBalanceFunctionMixed->SetTitle("Mixed events");
876 gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.4);
877 gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
878 gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
879 gHistBalanceFunctionMixed->DrawCopy("lego2");
880 gPad->SetTheta(30); // default is 30
881 gPad->SetPhi(-60); // default is 30
884 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
885 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
886 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
887 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
891 TCanvas *c4 = new TCanvas("c4","Subtracted balance function 2D",300,300,600,500);
892 c4->SetFillColor(10); c4->SetHighLightColor(10);
893 c4->SetLeftMargin(0.17); c4->SetTopMargin(0.05);
894 gHistBalanceFunctionSubtracted->SetTitle("Subtracted balance function");
895 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.4);
896 gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
897 gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
898 gHistBalanceFunctionSubtracted->DrawCopy("lego2");
899 gPad->SetTheta(30); // default is 30
900 gPad->SetPhi(-60); // default is 30
903 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
904 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
905 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
906 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
910 //____________________________________________________________//
911 void drawProjections(const char* lhcPeriod = "LHC10h",
912 const char* gCentralityEstimator = "V0M",
914 const char* gEventPlaneEstimator = "VZERO",
915 Bool_t kProjectInEta = kFALSE,
918 Int_t gCentrality = 1,
919 Double_t psiMin = -0.5,
920 Double_t psiMax = 3.5,
921 Double_t vertexZMin = -10.,
922 Double_t vertexZMax = 10.,
923 Double_t ptTriggerMin = -1.,
924 Double_t ptTriggerMax = -1.,
925 Double_t ptAssociatedMin = -1.,
926 Double_t ptAssociatedMax = -1.,
927 Bool_t kUseZYAM = kFALSE,
928 Bool_t k2pMethod = kTRUE,
929 TString eventClass = "Centrality",
930 Bool_t bRootMoments = kFALSE) {
931 //Macro that draws the charge dependent correlation functions PROJECTIONS
932 //for each centrality bin for the different pT of trigger and
933 //associated particles
934 TGaxis::SetMaxDigits(3);
936 //first we need some libraries
937 gSystem->Load("libTree");
938 gSystem->Load("libGeom");
939 gSystem->Load("libVMC");
940 gSystem->Load("libXMLIO");
941 gSystem->Load("libPhysics");
943 gSystem->Load("libSTEERBase");
944 gSystem->Load("libESD");
945 gSystem->Load("libAOD");
947 gSystem->Load("libANALYSIS.so");
948 gSystem->Load("libANALYSISalice.so");
949 gSystem->Load("libEventMixing.so");
950 gSystem->Load("libCORRFW.so");
951 gSystem->Load("libPWGTools.so");
952 gSystem->Load("libPWGCFebye.so");
955 TString filename = "balanceFunction2D.";
956 if(eventClass == "Centrality"){
957 filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
958 filename += ".PsiAll.PttFrom";
960 else if(eventClass == "Multiplicity"){
961 filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
962 filename += ".PsiAll.PttFrom";
964 else{ // "EventPlane" (default)
965 filename += "Centrality";
966 filename += gCentrality; filename += ".Psi";
967 if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
968 else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
969 else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
970 else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
971 else filename += "All.PttFrom";
973 filename += Form("%.1f",ptTriggerMin); filename += "To";
974 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
975 filename += Form("%.1f",ptAssociatedMin); filename += "To";
976 filename += Form("%.1f",ptAssociatedMax);
977 if(k2pMethod) filename += "_2pMethod";
980 // filename += Form("%.1f",psiMin);
982 // filename += Form("%.1f",psiMax);
986 TFile *f = TFile::Open(filename.Data());
987 if((!f)||(!f->IsOpen())) {
988 Printf("The file %s is not found. Aborting...",filename);
994 TString centralityLatex = "Centrality: ";
995 centralityLatex += centralityArray[gCentrality-1];
996 centralityLatex += "%";
999 if((psiMin == -0.5)&&(psiMax == 0.5))
1000 psiLatex = " -7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o}";
1001 else if((psiMin == 0.5)&&(psiMax == 1.5))
1002 psiLatex = " 37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o}";
1003 else if((psiMin == 1.5)&&(psiMax == 2.5))
1004 psiLatex = " 82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o}";
1006 psiLatex = " 0^{o} < #varphi^{t} - #Psi_{2} < 180^{o}";
1008 TString pttLatex = Form("%.1f",ptTriggerMin);
1009 pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
1010 pttLatex += " GeV/c";
1012 TString ptaLatex = Form("%.1f",ptAssociatedMin);
1013 ptaLatex += " < p_{T,assoc} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1014 ptaLatex += " GeV/c";
1016 TLatex *latexInfo1 = new TLatex();
1017 latexInfo1->SetNDC();
1018 latexInfo1->SetTextSize(0.045);
1019 latexInfo1->SetTextColor(1);
1023 //============================================================//
1024 //Get subtracted and mixed balance function
1026 TH2D *gHistBalanceFunctionSubtracted2D = (TH2D*)f->Get("gHistBalanceFunctionSubtracted");
1027 TH2D *gHistBalanceFunctionMixed2D = (TH2D*)f->Get("gHistBalanceFunctionMixed");
1029 TH1D *gHistBalanceFunctionSubtracted = NULL;
1030 TH1D *gHistBalanceFunctionMixed = NULL;
1032 TH1D *gHistBalanceFunctionSubtracted_scale = NULL;
1033 TH1D *gHistBalanceFunctionMixed_scale = NULL;
1036 gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionX());
1037 gHistBalanceFunctionSubtracted->Scale(gHistBalanceFunctionSubtracted2D->GetYaxis()->GetBinWidth(1)); // to remove normalization to phi bin width
1038 gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionX());
1039 gHistBalanceFunctionMixed->Scale(gHistBalanceFunctionMixed2D->GetYaxis()->GetBinWidth(1)); // to remove normalization to phi bin width
1040 gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#eta)");
1041 gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#eta)");
1044 gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionY());
1045 gHistBalanceFunctionSubtracted->Scale(gHistBalanceFunctionSubtracted2D->GetXaxis()->GetBinWidth(1)); // to remove normalization to eta bin width
1046 gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionY());
1047 gHistBalanceFunctionMixed->Scale(gHistBalanceFunctionMixed2D->GetXaxis()->GetBinWidth(1)); // to remove normalization to eta bin width
1048 gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#varphi)");
1049 gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#varphi)");
1052 gHistBalanceFunctionSubtracted->SetMarkerStyle(20);
1053 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
1054 gHistBalanceFunctionSubtracted->SetName("gHistBalanceFunctionSubtracted");
1056 gHistBalanceFunctionMixed->SetMarkerStyle(25);
1057 gHistBalanceFunctionMixed->SetName("gHistBalanceFunctionMixed");
1059 TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
1060 c1->SetFillColor(10);
1061 c1->SetHighLightColor(10);
1062 c1->SetLeftMargin(0.15);
1063 gHistBalanceFunctionSubtracted->DrawCopy("E");
1064 gHistBalanceFunctionMixed->DrawCopy("ESAME");
1066 legend = new TLegend(0.18,0.62,0.45,0.82,"","brNDC");
1067 legend->SetTextSize(0.045);
1068 legend->SetTextFont(42);
1069 legend->SetBorderSize(0);
1070 legend->SetFillStyle(0);
1071 legend->SetFillColor(10);
1072 legend->SetMargin(0.25);
1073 legend->SetShadowColor(10);
1074 legend->AddEntry(gHistBalanceFunctionSubtracted,"Data","lp");
1075 legend->AddEntry(gHistBalanceFunctionMixed,"Mixed data","lp");
1078 pngName = "BalanceFunction.";
1079 if(eventClass == "Centrality"){
1080 pngName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1081 if(kProjectInEta) pngName += ".InDeltaEta.PsiAll.PttFrom";
1082 else pngName += ".InDeltaPhi.PsiAll.PttFrom";
1084 else if(eventClass == "Multiplicity"){
1085 pngName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1086 if(kProjectInEta) pngName += ".InDeltaEta.PsiAll.PttFrom";
1087 else pngName += ".InDeltaPhi.PsiAll.PttFrom";
1089 else{ // "EventPlane" (default)
1090 pngName += "Centrality";
1091 pngName += gCentrality;
1092 if(kProjectInEta) pngName += ".InDeltaEta.Psi";
1093 else pngName += ".InDeltaPhi.Psi";
1094 if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
1095 else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
1096 else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
1097 else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.PttFrom";
1098 else pngName += "All.PttFrom";
1100 pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1101 pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1102 pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1103 pngName += Form("%.1f",ptAssociatedMax);
1104 if(k2pMethod) pngName += "_2pMethod";
1107 pngName += Form("%.1f",psiMin); pngName += "-";
1108 pngName += Form("%.1f",psiMax);
1111 c1->SaveAs(pngName.Data());
1113 TString meanLatex, rmsLatex, skewnessLatex, kurtosisLatex;
1116 meanLatex = "#mu = ";
1117 meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMean());
1118 meanLatex += " #pm ";
1119 meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMeanError());
1121 rmsLatex = "#sigma = ";
1122 rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMS());
1123 rmsLatex += " #pm ";
1124 rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMSError());
1126 skewnessLatex = "S = ";
1127 skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(1));
1128 skewnessLatex += " #pm ";
1129 skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(11));
1131 kurtosisLatex = "K = ";
1132 kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(1));
1133 kurtosisLatex += " #pm ";
1134 kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(11));
1135 Printf("Mean: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetMean(),gHistBalanceFunctionSubtracted->GetMeanError());
1136 Printf("RMS: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetRMS(),gHistBalanceFunctionSubtracted->GetRMSError());
1137 Printf("Skeweness: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetSkewness(1),gHistBalanceFunctionSubtracted->GetSkewness(11));
1138 Printf("Kurtosis: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetKurtosis(1),gHistBalanceFunctionSubtracted->GetKurtosis(11));
1140 // calculate the moments by hand
1143 Double_t meanAnalytical, meanAnalyticalError;
1144 Double_t sigmaAnalytical, sigmaAnalyticalError;
1145 Double_t skewnessAnalytical, skewnessAnalyticalError;
1146 Double_t kurtosisAnalytical, kurtosisAnalyticalError;
1148 Int_t gDeltaEtaPhi = 2;
1149 if(kProjectInEta) gDeltaEtaPhi = 1;
1151 AliBalancePsi *bHelper = new AliBalancePsi;
1152 bHelper->GetMomentsAnalytical(gDeltaEtaPhi,gHistBalanceFunctionSubtracted,meanAnalytical,meanAnalyticalError,sigmaAnalytical,sigmaAnalyticalError,skewnessAnalytical,skewnessAnalyticalError,kurtosisAnalytical,kurtosisAnalyticalError);
1154 meanLatex = "#mu = ";
1155 meanLatex += Form("%.3f",meanAnalytical);
1156 meanLatex += " #pm ";
1157 meanLatex += Form("%.3f",meanAnalyticalError);
1159 rmsLatex = "#sigma = ";
1160 rmsLatex += Form("%.3f",sigmaAnalytical);
1161 rmsLatex += " #pm ";
1162 rmsLatex += Form("%.3f",sigmaAnalyticalError);
1164 skewnessLatex = "S = ";
1165 skewnessLatex += Form("%.3f",skewnessAnalytical);
1166 skewnessLatex += " #pm ";
1167 skewnessLatex += Form("%.3f",skewnessAnalyticalError);
1169 kurtosisLatex = "K = ";
1170 kurtosisLatex += Form("%.3f",kurtosisAnalytical);
1171 kurtosisLatex += " #pm ";
1172 kurtosisLatex += Form("%.3f",kurtosisAnalyticalError);
1173 Printf("Mean: %lf - Error: %lf",meanAnalytical, meanAnalyticalError);
1174 Printf("Sigma: %lf - Error: %lf",sigmaAnalytical, sigmaAnalyticalError);
1175 Printf("Skeweness: %lf - Error: %lf",skewnessAnalytical, skewnessAnalyticalError);
1176 Printf("Kurtosis: %lf - Error: %lf",kurtosisAnalytical, kurtosisAnalyticalError);
1179 TCanvas *c2 = new TCanvas("c2","",600,0,600,500);
1180 c2->SetFillColor(10);
1181 c2->SetHighLightColor(10);
1182 c2->SetLeftMargin(0.15);
1183 gHistBalanceFunctionSubtracted->DrawCopy("E");
1185 TLatex *latex = new TLatex();
1187 latex->SetTextSize(0.035);
1188 latex->SetTextColor(1);
1189 latex->DrawLatex(0.64,0.85,meanLatex.Data());
1190 latex->DrawLatex(0.64,0.81,rmsLatex.Data());
1191 latex->DrawLatex(0.64,0.77,skewnessLatex.Data());
1192 latex->DrawLatex(0.64,0.73,kurtosisLatex.Data());
1194 TString outFileName = filename;
1195 if(kProjectInEta) outFileName.ReplaceAll(".root","_DeltaEtaProjection.root");
1196 else outFileName.ReplaceAll(".root","_DeltaPhiProjection.root");
1197 TFile *fProjection = TFile::Open(outFileName.Data(),"recreate");
1198 gHistBalanceFunctionSubtracted->Write();
1199 gHistBalanceFunctionMixed->Write();
1200 fProjection->Close();