[u/mrichter/AliRoot.git] / PWGCF / EBYE / macros / drawBalanceFunction2DPsi.C

const Int_t numberOfCentralityBins = 9;
TString centralityArray[numberOfCentralityBins] = {"0-5","5-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80"};
Double_t gMinCentrality[numberOfCentralityBins] = {0.,5.,10.,20.,30.,40.,50.,60.,70.};
Double_t gMaxCentrality[numberOfCentralityBins] = {5.,10.,20.,30.,40.,50.,60.,70.,80.};
TString gAnalysisType[7] = {"y","eta","qlong","qout","qside","qinv","phi"};

const Int_t gRebin = 1;
void drawBalanceFunction2DPsi(const char* filename = "AnalysisResultsPsi.root", 
			      Double_t psiMin = 0., Double_t psiMax = 7.5,
			      Double_t ptTriggerMin = -1.,
			      Double_t ptTriggerMax = -1.,
			      Double_t ptAssociatedMin = -1.,
			      Double_t ptAssociatedMax = -1.) {
  //Macro that draws the BF distributions for each centrality bin
  //for reaction plane dependent analysis
  //Author: Panos.Christakoglou@nikhef.nl
  //Load the PWG2ebye library
  gSystem->Load("libANALYSIS.so");
  gSystem->Load("libANALYSISalice.so");
  gSystem->Load("libEventMixing.so");
  gSystem->Load("libCORRFW.so");
  gSystem->Load("libPWGTools.so");
  gSystem->Load("libPWGCFebye.so");

  //Prepare the objects and return them
  TList *listBF = GetListOfObjects(filename);
  TList *listBFShuffled = GetListOfObjects(filename,kTRUE);
  if(!listBF) {
    Printf("The TList object was not created");
    return;
  }
  else 
    draw(listBF,listBFShuffled,psiMin,psiMax,
	 ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);  
}

//______________________________________________________//
TList *GetListOfObjects(const char* filename, Bool_t kShuffling = kFALSE) {
  //Get the TList objects (QA, bf, bf shuffled)
  TList *listQA = 0x0;
  TList *listBF = 0x0;
  TList *listBFShuffling = 0x0;
  
  //Open the file
  TFile *f = TFile::Open(filename);
  if((!f)||(!f->IsOpen())) {
    Printf("The file %s is not found. Aborting...",filename);
    return listBF;
  }
  //f->ls();
  
  TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
  if(!dir) {   
    Printf("The TDirectoryFile is not found. Aborting...",filename);
    return listBF;
  }
  //dir->ls();
  
  TString listBFName;
  if(!kShuffling) 
    listBFName = "listBFPsi";
  else if(kShuffling)
    listBFName = "listBFPsiShuffled";
  listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
  //listBF->ls();

  //Get the histograms
  TString histoName;
  if(!kShuffling)
    histoName = "fHistPV0M";
  else if(kShuffling)
    histoName = "fHistP_shuffleV0M";
  AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));  
  if(!fHistP) {
    Printf("fHistP %s not found!!!",histoName.Data());
    break;
  }
  fHistP->FillParent(); fHistP->DeleteContainers();

  if(!kShuffling)
    histoName = "fHistNV0M";
  else if(kShuffling)
    histoName = "fHistN_shuffleV0M";
  AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
  if(!fHistN) {
    Printf("fHistN %s not found!!!",histoName.Data());
    break;
  }
  fHistN->FillParent(); fHistN->DeleteContainers();
    
  if(!kShuffling)
    histoName = "fHistPNV0M";
  else if(kShuffling)
    histoName = "fHistPN_shuffleV0M";
  AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
  if(!fHistPN) {
    Printf("fHistPN %s not found!!!",histoName.Data());
    break;
  }
  fHistPN->FillParent(); fHistPN->DeleteContainers();
  
  if(!kShuffling)
    histoName = "fHistNPV0M";
  else if(kShuffling)
    histoName = "fHistNP_shuffleV0M";
  AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
  if(!fHistNP) {
    Printf("fHistNP %s not found!!!",histoName.Data());
    break;
  }
  fHistNP->FillParent(); fHistNP->DeleteContainers();

  if(!kShuffling)
    histoName = "fHistPPV0M";
  else if(kShuffling)
    histoName = "fHistPP_shuffleV0M";
  AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
  if(!fHistPP) {
    Printf("fHistPP %s not found!!!",histoName.Data());
    break;
  }
  fHistPP->FillParent(); fHistPP->DeleteContainers();

  if(!kShuffling)
    histoName = "fHistNNV0M";
  else if(kShuffling)
    histoName = "fHistNN_shuffleV0M";
  AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
  if(!fHistNN) {
    Printf("fHistNN %s not found!!!",histoName.Data());
    break;
  }
  fHistNN->FillParent(); fHistNN->DeleteContainers();
  
  return listBF;
}

//______________________________________________________//
void draw(TList *listBF, TList *listBFShuffled,
	  Double_t psiMin, Double_t psiMax,
	  Double_t ptTriggerMin, Double_t ptTriggerMax,
	  Double_t ptAssociatedMin, Double_t ptAssociatedMax) {
  gROOT->LoadMacro("~/SetPlotStyle.C");
  SetPlotStyle();
  gStyle->SetPalette(1,0);
  
  //balance function
  AliTHn *hP = NULL;
  AliTHn *hN = NULL;
  AliTHn *hPN = NULL;
  AliTHn *hNP = NULL;
  AliTHn *hPP = NULL;
  AliTHn *hNN = NULL;
  //listBF->ls();
  //Printf("=================");
  hP = (AliTHn*) listBF->FindObject("fHistPV0M");
  hN = (AliTHn*) listBF->FindObject("fHistNV0M");
  hPN = (AliTHn*) listBF->FindObject("fHistPNV0M");
  hNP = (AliTHn*) listBF->FindObject("fHistNPV0M");
  hPP = (AliTHn*) listBF->FindObject("fHistPPV0M");
  hNN = (AliTHn*) listBF->FindObject("fHistNNV0M");

  AliBalancePsi *b = new AliBalancePsi();
  b->SetHistNp(hP);
  b->SetHistNn(hN);
  b->SetHistNpn(hPN);
  b->SetHistNnp(hNP);
  b->SetHistNpp(hPP);
  b->SetHistNnn(hNN);

  //balance function shuffling
  AliTHn *hPShuffled = NULL;
  AliTHn *hNShuffled = NULL;
  AliTHn *hPNShuffled = NULL;
  AliTHn *hNPShuffled = NULL;
  AliTHn *hPPShuffled = NULL;
  AliTHn *hNNShuffled = NULL;
  //listBFShuffled->ls();
  hPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistP_shuffleV0M");
  hNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistN_shuffleV0M");
  hPNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPN_shuffleV0M");
  hNPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNP_shuffleV0M");
  hPPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPP_shuffleV0M");
  hNNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNN_shuffleV0M");

  AliBalancePsi *bShuffled = new AliBalancePsi();
  bShuffled->SetHistNp(hPShuffled);
  bShuffled->SetHistNn(hNShuffled);
  bShuffled->SetHistNpn(hPNShuffled);
  bShuffled->SetHistNnp(hNPShuffled);
  bShuffled->SetHistNpp(hPPShuffled);
  bShuffled->SetHistNnn(hNNShuffled);

  TH2D *gHistBalanceFunction;
  TH2D *gHistBalanceFunctionShuffled;
  TCanvas *c1;
  TString histoTitle, pngName;
  TLegend *legend;
  
  //loop over the centrality bins
  //for(Int_t iCentralityBin = 0; iCentralityBin < numberOfCentralityBins; iCentralityBin++) {  
  histoTitle = "Centrality: "; 
  histoTitle += centralityArray[6]; 
  histoTitle += "%";
  histoTitle += " | "; histoTitle += psiMin; 
  histoTitle += " < #phi - #Psi_{2} < "; histoTitle += psiMax;
  
  gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
  gHistBalanceFunction->SetTitle(histoTitle.Data());
  gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
  
  gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
  
  c1 = new TCanvas(histoTitle.Data(),"",0,0,600,500);
  c1->SetFillColor(10); 
  c1->SetHighLightColor(10);
  c1->SetLeftMargin(0.15);
  gHistBalanceFunction->Draw("lego2");
  //gHistBalanceFunctionShuffled->Draw("ESAME");
  
  legend = new TLegend(0.18,0.6,0.45,0.82,"","brNDC");
  legend->SetTextSize(0.045); 
  legend->SetTextFont(42); 
  legend->SetBorderSize(0);
  legend->SetFillStyle(0); 
  legend->SetFillColor(10);
  legend->SetMargin(0.25); 
  legend->SetShadowColor(10);
  legend->AddEntry(gHistBalanceFunction,"Data","lp");
  legend->AddEntry(gHistBalanceFunctionShuffled,"Shuffled data","lp");
  legend->Draw();
  
  pngName = "BalanceFunctionDeltaEta.Centrality"; 
  pngName += centralityArray[6]; 
  pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
  pngName += ".png";
  c1->SaveAs(pngName.Data());
}

//____________________________________________________________________//
void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1) {
  //Prints the calculated width of the BF and its error
  Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
  Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
  Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
  Double_t deltaBalP2 = 0.0, integral = 0.0;
  Double_t deltaErrorNew = 0.0;

  //Retrieve this variables from Histogram
  Int_t fNumberOfBins = gHistBalance->GetNbinsX();
  Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!
  
  cout<<"=================================================="<<endl;
  cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
  cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
  cout<<"=================================================="<<endl;
  for(Int_t i = 1; i <= fNumberOfBins; i++) {

    // this is to simulate |Delta eta| or |Delta phi|
    if(fNumberOfBins/2 - fStartBin + 1 < i && i < fNumberOfBins/2 + fStartBin ) continue;

    cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<TMath::Abs(gHistBalance->GetBinCenter(i))<<endl;

    gSumXi += TMath::Abs(gHistBalance->GetBinCenter(i)); // this is to simulate |Delta eta| or |Delta phi|
    gSumBi += gHistBalance->GetBinContent(i); 
    gSumBiXi += gHistBalance->GetBinContent(i)*TMath::Abs(gHistBalance->GetBinCenter(i));
    gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2);
    gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2);
    gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
    gSumXi2DeltaBi2 += TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2) * TMath::Power(gHistBalance->GetBinError(i),2);
    
    deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
    integral += fP2Step*gHistBalance->GetBinContent(i);
  }
  for(Int_t i = fStartBin; i < fNumberOfBins; i++)
    deltaErrorNew += gHistBalance->GetBinError(i)*(TMath::Abs(gHistBalance->GetBinCenter(i))*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
  
  Double_t integralError = TMath::Sqrt(deltaBalP2);
  
  Double_t delta = gSumBiXi / gSumBi;
  Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
  cout<<"=================================================="<<endl;
  cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
  cout<<"New error: "<<deltaErrorNew<<endl;
  cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
  cout<<"=================================================="<<endl;
  cout<<endl;
}
Commit	Line	Data
6acdbcb2	1	const Int_t numberOfCentralityBins = 9;
	2	TString centralityArray[numberOfCentralityBins] = {"0-5","5-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80"};
	3	Double_t gMinCentrality[numberOfCentralityBins] = {0.,5.,10.,20.,30.,40.,50.,60.,70.};
	4	Double_t gMaxCentrality[numberOfCentralityBins] = {5.,10.,20.,30.,40.,50.,60.,70.,80.};
	5	TString gAnalysisType[7] = {"y","eta","qlong","qout","qside","qinv","phi"};
	6
	7	const Int_t gRebin = 1;
	8	void drawBalanceFunction2DPsi(const char* filename = "AnalysisResultsPsi.root",
	9	Double_t psiMin = 0., Double_t psiMax = 7.5,
	10	Double_t ptTriggerMin = -1.,
	11	Double_t ptTriggerMax = -1.,
	12	Double_t ptAssociatedMin = -1.,
	13	Double_t ptAssociatedMax = -1.) {
	14	//Macro that draws the BF distributions for each centrality bin
	15	//for reaction plane dependent analysis
	16	//Author: Panos.Christakoglou@nikhef.nl
	17	//Load the PWG2ebye library
	18	gSystem->Load("libANALYSIS.so");
	19	gSystem->Load("libANALYSISalice.so");
	20	gSystem->Load("libEventMixing.so");
	21	gSystem->Load("libCORRFW.so");
	22	gSystem->Load("libPWGTools.so");
	23	gSystem->Load("libPWGCFebye.so");
	24
	25	//Prepare the objects and return them
	26	TList *listBF = GetListOfObjects(filename);
	27	TList *listBFShuffled = GetListOfObjects(filename,kTRUE);
	28	if(!listBF) {
	29	Printf("The TList object was not created");
	30	return;
	31	}
	32	else
	33	draw(listBF,listBFShuffled,psiMin,psiMax,
	34	ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
	35	}
	36
	37	//______________________________________________________//
	38	TList GetListOfObjects(const char filename, Bool_t kShuffling = kFALSE) {
	39	//Get the TList objects (QA, bf, bf shuffled)
	40	TList *listQA = 0x0;
	41	TList *listBF = 0x0;
	42	TList *listBFShuffling = 0x0;
	43
	44	//Open the file
	45	TFile *f = TFile::Open(filename);
	46	if((!f)\|\|(!f->IsOpen())) {
	47	Printf("The file %s is not found. Aborting...",filename);
	48	return listBF;
	49	}
	50	//f->ls();
	51
	52	TDirectoryFile dir = dynamic_cast<TDirectoryFile >(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
	53	if(!dir) {
	54	Printf("The TDirectoryFile is not found. Aborting...",filename);
	55	return listBF;
	56	}
	57	//dir->ls();
	58
	59	TString listBFName;
	60	if(!kShuffling)
	61	listBFName = "listBFPsi";
	62	else if(kShuffling)
	63	listBFName = "listBFPsiShuffled";
	64	listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
65	//listBF->ls();
66
67	//Get the histograms
68	TString histoName;
69	if(!kShuffling)
70	histoName = "fHistPV0M";
71	else if(kShuffling)
72	histoName = "fHistP_shuffleV0M";
73	AliTHn fHistP = dynamic_cast<AliTHn >(listBF->FindObject(histoName.Data()));
74	if(!fHistP) {
75	Printf("fHistP %s not found!!!",histoName.Data());
76	break;
77	}
78	fHistP->FillParent(); fHistP->DeleteContainers();
79
80	if(!kShuffling)
81	histoName = "fHistNV0M";
82	else if(kShuffling)
83	histoName = "fHistN_shuffleV0M";
84	AliTHn fHistN = dynamic_cast<AliTHn >(listBF->FindObject(histoName.Data()));
85	if(!fHistN) {
86	Printf("fHistN %s not found!!!",histoName.Data());
87	break;
88	}
89	fHistN->FillParent(); fHistN->DeleteContainers();
90
91	if(!kShuffling)
92	histoName = "fHistPNV0M";
93	else if(kShuffling)
94	histoName = "fHistPN_shuffleV0M";
95	AliTHn fHistPN = dynamic_cast<AliTHn >(listBF->FindObject(histoName.Data()));
96	if(!fHistPN) {
97	Printf("fHistPN %s not found!!!",histoName.Data());
98	break;
99	}
100	fHistPN->FillParent(); fHistPN->DeleteContainers();
101
102	if(!kShuffling)
103	histoName = "fHistNPV0M";
104	else if(kShuffling)
105	histoName = "fHistNP_shuffleV0M";
106	AliTHn fHistNP = dynamic_cast<AliTHn >(listBF->FindObject(histoName.Data()));
107	if(!fHistNP) {
108	Printf("fHistNP %s not found!!!",histoName.Data());
109	break;
110	}
111	fHistNP->FillParent(); fHistNP->DeleteContainers();
112
113	if(!kShuffling)
114	histoName = "fHistPPV0M";
115	else if(kShuffling)
116	histoName = "fHistPP_shuffleV0M";
117	AliTHn fHistPP = dynamic_cast<AliTHn >(listBF->FindObject(histoName.Data()));
118	if(!fHistPP) {
119	Printf("fHistPP %s not found!!!",histoName.Data());
120	break;
121	}
122	fHistPP->FillParent(); fHistPP->DeleteContainers();
123
124	if(!kShuffling)
125	histoName = "fHistNNV0M";
126	else if(kShuffling)
127	histoName = "fHistNN_shuffleV0M";
128	AliTHn fHistNN = dynamic_cast<AliTHn >(listBF->FindObject(histoName.Data()));
129	if(!fHistNN) {
130	Printf("fHistNN %s not found!!!",histoName.Data());
131	break;
132	}
133	fHistNN->FillParent(); fHistNN->DeleteContainers();
134
135	return listBF;
136	}
137
138	//______________________________________________________//
139	void draw(TList listBF, TList listBFShuffled,
140	Double_t psiMin, Double_t psiMax,
141	Double_t ptTriggerMin, Double_t ptTriggerMax,
142	Double_t ptAssociatedMin, Double_t ptAssociatedMax) {
143	gROOT->LoadMacro("~/SetPlotStyle.C");
144	SetPlotStyle();
145	gStyle->SetPalette(1,0);
146
147	//balance function
148	AliTHn *hP = NULL;
149	AliTHn *hN = NULL;
150	AliTHn *hPN = NULL;
151	AliTHn *hNP = NULL;
152	AliTHn *hPP = NULL;
153	AliTHn *hNN = NULL;
154	//listBF->ls();
155	//Printf("=================");
156	hP = (AliTHn*) listBF->FindObject("fHistPV0M");
157	hN = (AliTHn*) listBF->FindObject("fHistNV0M");
158	hPN = (AliTHn*) listBF->FindObject("fHistPNV0M");
159	hNP = (AliTHn*) listBF->FindObject("fHistNPV0M");
160	hPP = (AliTHn*) listBF->FindObject("fHistPPV0M");
161	hNN = (AliTHn*) listBF->FindObject("fHistNNV0M");
162
163	AliBalancePsi *b = new AliBalancePsi();
164	b->SetHistNp(hP);
165	b->SetHistNn(hN);
166	b->SetHistNpn(hPN);
167	b->SetHistNnp(hNP);
168	b->SetHistNpp(hPP);
169	b->SetHistNnn(hNN);
170
171	//balance function shuffling
172	AliTHn *hPShuffled = NULL;
173	AliTHn *hNShuffled = NULL;
174	AliTHn *hPNShuffled = NULL;
175	AliTHn *hNPShuffled = NULL;
176	AliTHn *hPPShuffled = NULL;
177	AliTHn *hNNShuffled = NULL;
178	//listBFShuffled->ls();
179	hPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistP_shuffleV0M");
180	hNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistN_shuffleV0M");
181	hPNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPN_shuffleV0M");
182	hNPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNP_shuffleV0M");
183	hPPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPP_shuffleV0M");
184	hNNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNN_shuffleV0M");
185
186	AliBalancePsi *bShuffled = new AliBalancePsi();
187	bShuffled->SetHistNp(hPShuffled);
188	bShuffled->SetHistNn(hNShuffled);
189	bShuffled->SetHistNpn(hPNShuffled);
190	bShuffled->SetHistNnp(hNPShuffled);
191	bShuffled->SetHistNpp(hPPShuffled);
192	bShuffled->SetHistNnn(hNNShuffled);
193
194	TH2D *gHistBalanceFunction;
195	TH2D *gHistBalanceFunctionShuffled;
196	TCanvas *c1;
197	TString histoTitle, pngName;
198	TLegend *legend;
199
200	//loop over the centrality bins
201	//for(Int_t iCentralityBin = 0; iCentralityBin < numberOfCentralityBins; iCentralityBin++) {
202	histoTitle = "Centrality: ";
203	histoTitle += centralityArray[6];
204	histoTitle += "%";
205	histoTitle += " \| "; histoTitle += psiMin;
206	histoTitle += " < #phi - #Psi_{2} < "; histoTitle += psiMax;
207
208	gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
209	gHistBalanceFunction->SetTitle(histoTitle.Data());
210	gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
211
212	gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
213
214	c1 = new TCanvas(histoTitle.Data(),"",0,0,600,500);
215	c1->SetFillColor(10);
216	c1->SetHighLightColor(10);
217	c1->SetLeftMargin(0.15);
218	gHistBalanceFunction->Draw("lego2");
219	//gHistBalanceFunctionShuffled->Draw("ESAME");
220
221	legend = new TLegend(0.18,0.6,0.45,0.82,"","brNDC");
222	legend->SetTextSize(0.045);
223	legend->SetTextFont(42);
224	legend->SetBorderSize(0);
225	legend->SetFillStyle(0);
226	legend->SetFillColor(10);
227	legend->SetMargin(0.25);
228	legend->SetShadowColor(10);
229	legend->AddEntry(gHistBalanceFunction,"Data","lp");
230	legend->AddEntry(gHistBalanceFunctionShuffled,"Shuffled data","lp");
231	legend->Draw();
232
233	pngName = "BalanceFunctionDeltaEta.Centrality";
234	pngName += centralityArray[6];
235	pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
236	pngName += ".png";
237	c1->SaveAs(pngName.Data());
238	}
239
240	//____________________________________________________________________//
241	void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1) {
242	//Prints the calculated width of the BF and its error
243	Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
244	Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
245	Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
246	Double_t deltaBalP2 = 0.0, integral = 0.0;
247	Double_t deltaErrorNew = 0.0;
248
249	//Retrieve this variables from Histogram
250	Int_t fNumberOfBins = gHistBalance->GetNbinsX();
251	Double_t fP2Step = gHistBalance->GetBinWidth(1); // assume equal binning!
252
253	cout<<"=================================================="<<endl;
254	cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
255	cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
256	cout<<"=================================================="<<endl;
257	for(Int_t i = 1; i <= fNumberOfBins; i++) {
258
259	// this is to simulate \|Delta eta\| or \|Delta phi\|
260	if(fNumberOfBins/2 - fStartBin + 1 < i && i < fNumberOfBins/2 + fStartBin ) continue;
261
262	cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<TMath::Abs(gHistBalance->GetBinCenter(i))<<endl;
263
264	gSumXi += TMath::Abs(gHistBalance->GetBinCenter(i)); // this is to simulate \|Delta eta\| or \|Delta phi\|
265	gSumBi += gHistBalance->GetBinContent(i);
266	gSumBiXi += gHistBalance->GetBinContent(i)*TMath::Abs(gHistBalance->GetBinCenter(i));
267	gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2);
268	gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2);
269	gSumDeltaBi2 += TMath::Power(gHistBalance->GetBinError(i),2);
270	gSumXi2DeltaBi2 += TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2) * TMath::Power(gHistBalance->GetBinError(i),2);
271
272	deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
273	integral += fP2Step*gHistBalance->GetBinContent(i);
274	}
275	for(Int_t i = fStartBin; i < fNumberOfBins; i++)
276	deltaErrorNew += gHistBalance->GetBinError(i)(TMath::Abs(gHistBalance->GetBinCenter(i))gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
277
278	Double_t integralError = TMath::Sqrt(deltaBalP2);
279
280	Double_t delta = gSumBiXi / gSumBi;
281	Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
282	cout<<"=================================================="<<endl;
283	cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
284	cout<<"New error: "<<deltaErrorNew<<endl;
285	cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
286	cout<<"=================================================="<<endl;
287	cout<<endl;
288	}