]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Changes in drawing macros:
[u/mrichter/AliRoot.git] / PWGCF / EBYE / macros / drawBalanceFunction2DPsi.C
CommitLineData
e4391c02 1const Int_t numberOfCentralityBins = 12;
7827f01e 2TString centralityArray[numberOfCentralityBins] = {"0-80","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-100","0-1","1-2","2-3"};
e4391c02 3
6acdbcb2 4
5const Int_t gRebin = 1;
6void drawBalanceFunction2DPsi(const char* filename = "AnalysisResultsPsi.root",
eb63b883 7 Int_t gCentrality = 1,
5de9ad1a 8 Int_t gBit = -1,
9 const char* gCentralityEstimator = 0x0,
52daf7b2 10 Bool_t kShowShuffled = kFALSE,
11 Bool_t kShowMixed = kTRUE,
eb63b883 12 Double_t psiMin = -0.5, Double_t psiMax = 0.5,
34616eda 13 Double_t vertexZMin = -10.,
14 Double_t vertexZMax = 10.,
6acdbcb2 15 Double_t ptTriggerMin = -1.,
16 Double_t ptTriggerMax = -1.,
17 Double_t ptAssociatedMin = -1.,
f0c5040c 18 Double_t ptAssociatedMax = -1.,
832b21d8 19 Bool_t kUseVzBinning = kTRUE,
20 Bool_t k2pMethod = kTRUE,
3de51a46 21 TString eventClass = "EventPlane", //Can be "EventPlane", "Centrality", "Multiplicity"
22 Bool_t bToy = kFALSE)
d67eae53 23{
6acdbcb2 24 //Macro that draws the BF distributions for each centrality bin
25 //for reaction plane dependent analysis
26 //Author: Panos.Christakoglou@nikhef.nl
27 //Load the PWG2ebye library
28 gSystem->Load("libANALYSIS.so");
29 gSystem->Load("libANALYSISalice.so");
30 gSystem->Load("libEventMixing.so");
31 gSystem->Load("libCORRFW.so");
32 gSystem->Load("libPWGTools.so");
33 gSystem->Load("libPWGCFebye.so");
34
9a0495b0 35 //gROOT->LoadMacro("~/SetPlotStyle.C");
36 //SetPlotStyle();
8795e993 37 gStyle->SetPalette(1,0);
15dd45a0 38
6acdbcb2 39 //Prepare the objects and return them
3de51a46 40 TList *listBF = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,0,bToy);
52daf7b2 41 TList *listBFShuffled = NULL;
3de51a46 42 if(kShowShuffled) listBFShuffled = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,1,bToy);
52daf7b2 43 TList *listBFMixed = NULL;
3de51a46 44 if(kShowMixed) listBFMixed = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,2,bToy);
6acdbcb2 45 if(!listBF) {
46 Printf("The TList object was not created");
47 return;
48 }
49 else
27634c13 50 draw(listBF,listBFShuffled,listBFMixed,gCentrality,gCentralityEstimator,
34616eda 51 psiMin,psiMax,vertexZMin,vertexZMax,
f0c5040c 52 ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,
832b21d8 53 kUseVzBinning,k2pMethod,eventClass);
6acdbcb2 54}
55
56//______________________________________________________//
5365d1d7 57TList *GetListOfObjects(const char* filename,
58 Int_t gCentrality,
59 Int_t gBit,
60 const char *gCentralityEstimator,
3de51a46 61 Int_t kData = 1,
62 Bool_t bToy = kFALSE) {
6acdbcb2 63 //Get the TList objects (QA, bf, bf shuffled)
6acdbcb2 64 TList *listBF = 0x0;
6acdbcb2 65
66 //Open the file
5365d1d7 67 TFile *f = TFile::Open(filename,"UPDATE");
6acdbcb2 68 if((!f)||(!f->IsOpen())) {
69 Printf("The file %s is not found. Aborting...",filename);
70 return listBF;
71 }
72 //f->ls();
73
74 TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
75 if(!dir) {
76 Printf("The TDirectoryFile is not found. Aborting...",filename);
77 return listBF;
78 }
79 //dir->ls();
80
81 TString listBFName;
eb63b883 82 if(kData == 0) {
83 //cout<<"no shuffling - no mixing"<<endl;
3de51a46 84 listBFName = "listBFPsi";
eb63b883 85 }
86 else if(kData == 1) {
87 //cout<<"shuffling - no mixing"<<endl;
3de51a46 88 listBFName = "listBFPsiShuffled";
eb63b883 89 }
90 else if(kData == 2) {
91 //cout<<"no shuffling - mixing"<<endl;
3de51a46 92 listBFName = "listBFPsiMixed";
93 }
94
95 // different list names in case of toy model
96 if(!bToy){
dd46a0c3 97 listBFName += "_";
3de51a46 98 listBFName += centralityArray[gCentrality-1];
99 if(gBit > -1) {
100 listBFName += "_Bit"; listBFName += gBit; }
101 if(gCentralityEstimator) {
102 listBFName += "_"; listBFName += gCentralityEstimator;}
103 }
104 else{
105 listBFName.ReplaceAll("Psi","");
eb63b883 106 }
6acdbcb2 107
5365d1d7 108 // histograms were already retrieved (in first iteration)
109 if(dir->Get(Form("%s_histograms",listBFName.Data()))){
110 //cout<<"second iteration"<<endl;
111 listBF = dynamic_cast<TList *>(dir->Get(Form("%s_histograms",listBFName.Data())));
6acdbcb2 112 }
6acdbcb2 113
5365d1d7 114 // histograms were not yet retrieved (this is the first iteration)
115 else{
116
117 listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
118 cout<<"======================================================="<<endl;
5365d1d7 119 cout<<"List name: "<<listBF->GetName()<<endl;
120 //listBF->ls();
6acdbcb2 121
5365d1d7 122 //Get the histograms
123 TString histoName;
124 if(kData == 0)
27634c13 125 histoName = "fHistP";
5365d1d7 126 else if(kData == 1)
27634c13 127 histoName = "fHistP_shuffle";
5365d1d7 128 else if(kData == 2)
27634c13 129 histoName = "fHistP";
130 if(gCentralityEstimator) histoName += gCentralityEstimator;
5365d1d7 131 AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
132 if(!fHistP) {
133 Printf("fHistP %s not found!!!",histoName.Data());
134 break;
135 }
136 fHistP->FillParent(); fHistP->DeleteContainers();
137
138 if(kData == 0)
27634c13 139 histoName = "fHistN";
5365d1d7 140 if(kData == 1)
27634c13 141 histoName = "fHistN_shuffle";
5365d1d7 142 if(kData == 2)
27634c13 143 histoName = "fHistN";
144 if(gCentralityEstimator) histoName += gCentralityEstimator;
5365d1d7 145 AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
146 if(!fHistN) {
147 Printf("fHistN %s not found!!!",histoName.Data());
148 break;
149 }
150 fHistN->FillParent(); fHistN->DeleteContainers();
151
152 if(kData == 0)
27634c13 153 histoName = "fHistPN";
5365d1d7 154 if(kData == 1)
27634c13 155 histoName = "fHistPN_shuffle";
5365d1d7 156 if(kData == 2)
27634c13 157 histoName = "fHistPN";
158 if(gCentralityEstimator) histoName += gCentralityEstimator;
5365d1d7 159 AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
160 if(!fHistPN) {
161 Printf("fHistPN %s not found!!!",histoName.Data());
162 break;
163 }
164 fHistPN->FillParent(); fHistPN->DeleteContainers();
165
166 if(kData == 0)
27634c13 167 histoName = "fHistNP";
5365d1d7 168 if(kData == 1)
27634c13 169 histoName = "fHistNP_shuffle";
5365d1d7 170 if(kData == 2)
27634c13 171 histoName = "fHistNP";
172 if(gCentralityEstimator) histoName += gCentralityEstimator;
5365d1d7 173 AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
174 if(!fHistNP) {
175 Printf("fHistNP %s not found!!!",histoName.Data());
176 break;
177 }
178 fHistNP->FillParent(); fHistNP->DeleteContainers();
179
180 if(kData == 0)
27634c13 181 histoName = "fHistPP";
5365d1d7 182 if(kData == 1)
27634c13 183 histoName = "fHistPP_shuffle";
5365d1d7 184 if(kData == 2)
27634c13 185 histoName = "fHistPP";
186 if(gCentralityEstimator) histoName += gCentralityEstimator;
5365d1d7 187 AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
188 if(!fHistPP) {
189 Printf("fHistPP %s not found!!!",histoName.Data());
190 break;
191 }
192 fHistPP->FillParent(); fHistPP->DeleteContainers();
193
194 if(kData == 0)
27634c13 195 histoName = "fHistNN";
5365d1d7 196 if(kData == 1)
27634c13 197 histoName = "fHistNN_shuffle";
5365d1d7 198 if(kData == 2)
27634c13 199 histoName = "fHistNN";
200 if(gCentralityEstimator) histoName += gCentralityEstimator;
5365d1d7 201 AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
202 if(!fHistNN) {
203 Printf("fHistNN %s not found!!!",histoName.Data());
204 break;
205 }
206 fHistNN->FillParent(); fHistNN->DeleteContainers();
207
208 dir->cd();
209 listBF->Write(Form("%s_histograms",listBFName.Data()), TObject::kSingleKey);
210
211 }// first iteration
6acdbcb2 212
5365d1d7 213 f->Close();
6acdbcb2 214
215 return listBF;
216}
217
218//______________________________________________________//
eb63b883 219void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
27634c13 220 Int_t gCentrality, const char* gCentralityEstimator,
221 Double_t psiMin, Double_t psiMax,
34616eda 222 Double_t vertexZMin,
223 Double_t vertexZMax,
6acdbcb2 224 Double_t ptTriggerMin, Double_t ptTriggerMax,
f0c5040c 225 Double_t ptAssociatedMin, Double_t ptAssociatedMax,
832b21d8 226 Bool_t kUseVzBinning=kFALSE,
d67eae53 227 Bool_t k2pMethod = kFALSE, TString eventClass) {
6acdbcb2 228 //balance function
229 AliTHn *hP = NULL;
230 AliTHn *hN = NULL;
231 AliTHn *hPN = NULL;
232 AliTHn *hNP = NULL;
233 AliTHn *hPP = NULL;
234 AliTHn *hNN = NULL;
235 //listBF->ls();
236 //Printf("=================");
27634c13 237 TString histoName = "fHistP";
238 if(gCentralityEstimator) histoName += gCentralityEstimator;
239 hP = (AliTHn*) listBF->FindObject(histoName.Data());
52daf7b2 240 hP->SetName("gHistP");
27634c13 241 histoName = "fHistN";
242 if(gCentralityEstimator) histoName += gCentralityEstimator;
243 hN = (AliTHn*) listBF->FindObject(histoName.Data());
52daf7b2 244 hN->SetName("gHistN");
27634c13 245 histoName = "fHistPN";
246 if(gCentralityEstimator) histoName += gCentralityEstimator;
247 hPN = (AliTHn*) listBF->FindObject(histoName.Data());
52daf7b2 248 hPN->SetName("gHistPN");
27634c13 249 histoName = "fHistNP";
250 if(gCentralityEstimator) histoName += gCentralityEstimator;
251 hNP = (AliTHn*) listBF->FindObject(histoName.Data());
52daf7b2 252 hNP->SetName("gHistNP");
27634c13 253 histoName = "fHistPP";
254 if(gCentralityEstimator) histoName += gCentralityEstimator;
255 hPP = (AliTHn*) listBF->FindObject(histoName.Data());
52daf7b2 256 hPP->SetName("gHistPP");
27634c13 257 histoName = "fHistNN";
258 if(gCentralityEstimator) histoName += gCentralityEstimator;
259 hNN = (AliTHn*) listBF->FindObject(histoName.Data());
52daf7b2 260 hNN->SetName("gHistNN");
6acdbcb2 261
262 AliBalancePsi *b = new AliBalancePsi();
d67eae53 263 b->SetEventClass(eventClass);
6acdbcb2 264 b->SetHistNp(hP);
265 b->SetHistNn(hN);
266 b->SetHistNpn(hPN);
267 b->SetHistNnp(hNP);
268 b->SetHistNpp(hPP);
269 b->SetHistNnn(hNN);
832b21d8 270 if(kUseVzBinning) b->SetVertexZBinning(kTRUE);
271
6acdbcb2 272
273 //balance function shuffling
274 AliTHn *hPShuffled = NULL;
275 AliTHn *hNShuffled = NULL;
276 AliTHn *hPNShuffled = NULL;
277 AliTHn *hNPShuffled = NULL;
278 AliTHn *hPPShuffled = NULL;
279 AliTHn *hNNShuffled = NULL;
52daf7b2 280 if(listBFShuffled) {
281 //listBFShuffled->ls();
abf0c200 282 histoName = "fHistP_shuffle";
27634c13 283 if(gCentralityEstimator) histoName += gCentralityEstimator;
284 hPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
52daf7b2 285 hPShuffled->SetName("gHistPShuffled");
abf0c200 286 histoName = "fHistN_shuffle";
27634c13 287 if(gCentralityEstimator) histoName += gCentralityEstimator;
288 hNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
52daf7b2 289 hNShuffled->SetName("gHistNShuffled");
abf0c200 290 histoName = "fHistPN_shuffle";
27634c13 291 if(gCentralityEstimator) histoName += gCentralityEstimator;
292 hPNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
52daf7b2 293 hPNShuffled->SetName("gHistPNShuffled");
abf0c200 294 histoName = "fHistNP_shuffle";
27634c13 295 if(gCentralityEstimator) histoName += gCentralityEstimator;
296 hNPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
52daf7b2 297 hNPShuffled->SetName("gHistNPShuffled");
abf0c200 298 histoName = "fHistPP_shuffle";
27634c13 299 if(gCentralityEstimator) histoName += gCentralityEstimator;
300 hPPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
52daf7b2 301 hPPShuffled->SetName("gHistPPShuffled");
abf0c200 302 histoName = "fHistNN_shuffle";
27634c13 303 if(gCentralityEstimator) histoName += gCentralityEstimator;
304 hNNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
52daf7b2 305 hNNShuffled->SetName("gHistNNShuffled");
306
307 AliBalancePsi *bShuffled = new AliBalancePsi();
d67eae53 308 bShuffled->SetEventClass(eventClass);
52daf7b2 309 bShuffled->SetHistNp(hPShuffled);
310 bShuffled->SetHistNn(hNShuffled);
311 bShuffled->SetHistNpn(hPNShuffled);
312 bShuffled->SetHistNnp(hNPShuffled);
313 bShuffled->SetHistNpp(hPPShuffled);
314 bShuffled->SetHistNnn(hNNShuffled);
832b21d8 315 if(kUseVzBinning) bShuffled->SetVertexZBinning(kTRUE);
316
52daf7b2 317 }
6acdbcb2 318
eb63b883 319 //balance function mixing
320 AliTHn *hPMixed = NULL;
321 AliTHn *hNMixed = NULL;
322 AliTHn *hPNMixed = NULL;
323 AliTHn *hNPMixed = NULL;
324 AliTHn *hPPMixed = NULL;
325 AliTHn *hNNMixed = NULL;
eb63b883 326
52daf7b2 327 if(listBFMixed) {
328 //listBFMixed->ls();
27634c13 329 histoName = "fHistP";
330 if(gCentralityEstimator) histoName += gCentralityEstimator;
331 hPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
52daf7b2 332 hPMixed->SetName("gHistPMixed");
27634c13 333 histoName = "fHistN";
334 if(gCentralityEstimator) histoName += gCentralityEstimator;
335 hNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
52daf7b2 336 hNMixed->SetName("gHistNMixed");
27634c13 337 histoName = "fHistPN";
338 if(gCentralityEstimator) histoName += gCentralityEstimator;
339 hPNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
52daf7b2 340 hPNMixed->SetName("gHistPNMixed");
27634c13 341 histoName = "fHistNP";
342 if(gCentralityEstimator) histoName += gCentralityEstimator;
343 hNPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
344 histoName = "fHistNP";
345 if(gCentralityEstimator) histoName += gCentralityEstimator;
52daf7b2 346 hNPMixed->SetName("gHistNPMixed");
27634c13 347 histoName = "fHistPP";
348 if(gCentralityEstimator) histoName += gCentralityEstimator;
349 hPPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
52daf7b2 350 hPPMixed->SetName("gHistPPMixed");
27634c13 351 histoName = "fHistNN";
352 if(gCentralityEstimator) histoName += gCentralityEstimator;
353 hNNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
52daf7b2 354 hNNMixed->SetName("gHistNNMixed");
355
356 AliBalancePsi *bMixed = new AliBalancePsi();
d67eae53 357 bMixed->SetEventClass(eventClass);
52daf7b2 358 bMixed->SetHistNp(hPMixed);
359 bMixed->SetHistNn(hNMixed);
360 bMixed->SetHistNpn(hPNMixed);
361 bMixed->SetHistNnp(hNPMixed);
362 bMixed->SetHistNpp(hPPMixed);
363 bMixed->SetHistNnn(hNNMixed);
832b21d8 364 if(kUseVzBinning) bMixed->SetVertexZBinning(kTRUE);
365
52daf7b2 366 }
eb63b883 367
6acdbcb2 368 TH2D *gHistBalanceFunction;
eb63b883 369 TH2D *gHistBalanceFunctionSubtracted;
6acdbcb2 370 TH2D *gHistBalanceFunctionShuffled;
eb63b883 371 TH2D *gHistBalanceFunctionMixed;
6acdbcb2 372 TString histoTitle, pngName;
6acdbcb2 373
d67eae53 374 if(eventClass == "Centrality"){
375 histoTitle = "Centrality: ";
376 histoTitle += psiMin;
377 histoTitle += " - ";
378 histoTitle += psiMax;
379 histoTitle += " % ";
52daf7b2 380 histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
d67eae53 381 }
382 else if(eventClass == "Multiplicity"){
383 histoTitle = "Multiplicity: ";
384 histoTitle += psiMin;
385 histoTitle += " - ";
386 histoTitle += psiMax;
387 histoTitle += " tracks";
388 histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
389 }
390 else{ // "EventPlane" (default)
391 histoTitle = "Centrality: ";
392 histoTitle += centralityArray[gCentrality-1];
393 histoTitle += "%";
394 if((psiMin == -0.5)&&(psiMax == 0.5))
395 histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
396 else if((psiMin == 0.5)&&(psiMax == 1.5))
397 histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
398 else if((psiMin == 1.5)&&(psiMax == 2.5))
399 histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
400 else
401 histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
402 }
403
f0c5040c 404 if(k2pMethod)
405 if(bMixed)
34616eda 406 gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
f0c5040c 407 else{
408 cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
409 return;
410 }
411 else
34616eda 412 gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
6acdbcb2 413 gHistBalanceFunction->SetTitle(histoTitle.Data());
414 gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
eb63b883 415 gHistBalanceFunction->SetName("gHistBalanceFunction");
416
52daf7b2 417 if(listBFShuffled) {
f0c5040c 418
419 if(k2pMethod)
420 if(bMixed)
34616eda 421 gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
f0c5040c 422 else{
423 cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
424 return;
425 }
426 else
34616eda 427 gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
52daf7b2 428 gHistBalanceFunctionShuffled->SetTitle(histoTitle.Data());
429 gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
430 gHistBalanceFunctionShuffled->SetName("gHistBalanceFunctionShuffled");
431 }
eb63b883 432
52daf7b2 433 if(listBFMixed) {
f0c5040c 434 if(k2pMethod)
435 if(bMixed)
34616eda 436 gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
f0c5040c 437 else{
438 cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
439 return;
440 }
441 else
34616eda 442 gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
52daf7b2 443 gHistBalanceFunctionMixed->SetTitle(histoTitle.Data());
444 gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
445 gHistBalanceFunctionMixed->SetName("gHistBalanceFunctionMixed");
446
447 gHistBalanceFunctionSubtracted = dynamic_cast<TH2D *>(gHistBalanceFunction->Clone());
448 gHistBalanceFunctionSubtracted->Add(gHistBalanceFunctionMixed,-1);
449 gHistBalanceFunctionSubtracted->SetTitle(histoTitle.Data());
450 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
451 gHistBalanceFunctionSubtracted->SetName("gHistBalanceFunctionSubtracted");
452 }
eb63b883 453
454 //Draw the results
455 TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
6acdbcb2 456 c1->SetFillColor(10);
457 c1->SetHighLightColor(10);
458 c1->SetLeftMargin(0.15);
eb63b883 459 gHistBalanceFunction->DrawCopy("lego2");
5de9ad1a 460 gPad->SetTheta(30); // default is 30
461 gPad->SetPhi(-60); // default is 30
462 gPad->Update();
eb63b883 463 TCanvas *c1a = new TCanvas("c1a","",600,0,600,500);
464 c1a->SetFillColor(10);
465 c1a->SetHighLightColor(10);
466 c1a->SetLeftMargin(0.15);
467 gHistBalanceFunction->DrawCopy("colz");
468
52daf7b2 469 if(listBFShuffled) {
470 TCanvas *c2 = new TCanvas("c2","",100,100,600,500);
471 c2->SetFillColor(10);
472 c2->SetHighLightColor(10);
473 c2->SetLeftMargin(0.15);
474 gHistBalanceFunctionShuffled->DrawCopy("lego2");
5de9ad1a 475 gPad->SetTheta(30); // default is 30
476 gPad->SetPhi(-60); // default is 30
477 gPad->Update();
52daf7b2 478 TCanvas *c2a = new TCanvas("c2a","",700,100,600,500);
479 c2a->SetFillColor(10);
480 c2a->SetHighLightColor(10);
481 c2a->SetLeftMargin(0.15);
482 gHistBalanceFunctionShuffled->DrawCopy("colz");
483 }
eb63b883 484
52daf7b2 485 if(listBFMixed) {
486 TCanvas *c3 = new TCanvas("c3","",200,200,600,500);
487 c3->SetFillColor(10);
488 c3->SetHighLightColor(10);
489 c3->SetLeftMargin(0.15);
490 gHistBalanceFunctionMixed->DrawCopy("lego2");
5de9ad1a 491 gPad->SetTheta(30); // default is 30
492 gPad->SetPhi(-60); // default is 30
493 gPad->Update();
52daf7b2 494 TCanvas *c3a = new TCanvas("c3a","",800,200,600,500);
495 c3a->SetFillColor(10);
496 c3a->SetHighLightColor(10);
497 c3a->SetLeftMargin(0.15);
498 gHistBalanceFunctionMixed->DrawCopy("colz");
eb63b883 499
52daf7b2 500 TCanvas *c4 = new TCanvas("c4","",300,300,600,500);
501 c4->SetFillColor(10);
502 c4->SetHighLightColor(10);
503 c4->SetLeftMargin(0.15);
504 gHistBalanceFunctionSubtracted->DrawCopy("lego2");
5de9ad1a 505 gPad->SetTheta(30); // default is 30
506 gPad->SetPhi(-60); // default is 30
507 gPad->Update();
52daf7b2 508 TCanvas *c4a = new TCanvas("c4a","",900,300,600,500);
509 c4a->SetFillColor(10);
510 c4a->SetHighLightColor(10);
511 c4a->SetLeftMargin(0.15);
512 gHistBalanceFunctionSubtracted->DrawCopy("colz");
742af4bd 513
bd36d661 514 fitbalanceFunction(gCentrality, psiMin , psiMax,
742af4bd 515 ptTriggerMin, ptTriggerMax,
516 ptAssociatedMin, ptAssociatedMax,
d67eae53 517 gHistBalanceFunctionSubtracted,k2pMethod, eventClass);
52daf7b2 518 }
eb63b883 519
d67eae53 520 TString newFileName = "balanceFunction2D.";
521 if(eventClass == "Centrality"){
522 newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
523 newFileName += ".PsiAll.PttFrom";
524 }
525 else if(eventClass == "Multiplicity"){
526 newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
527 newFileName += ".PsiAll.PttFrom";
528 }
529 else{ // "EventPlane" (default)
530 newFileName += "Centrality";
531 newFileName += gCentrality; newFileName += ".Psi";
532 if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
533 else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
534 else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
535 else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
536 else newFileName += "All.PttFrom";
537 }
648f1a5a 538 newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
539 newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
540 newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
541 newFileName += Form("%.1f",ptAssociatedMax);
f0c5040c 542 if(k2pMethod) newFileName += "_2pMethod";
a9f20288 543
a2da950e 544 newFileName += "_";
545 newFileName += Form("%.1f",psiMin);
546 newFileName += "-";
547 newFileName += Form("%.1f",psiMax);
548 newFileName += ".root";
eb63b883 549
550 TFile *fOutput = new TFile(newFileName.Data(),"recreate");
551 fOutput->cd();
52daf7b2 552 /*hP->Write(); hN->Write();
553 hPN->Write(); hNP->Write();
554 hPP->Write(); hNN->Write();
555 hPShuffled->Write(); hNShuffled->Write();
556 hPNShuffled->Write(); hNPShuffled->Write();
557 hPPShuffled->Write(); hNNShuffled->Write();
558 hPMixed->Write(); hNMixed->Write();
559 hPNMixed->Write(); hNPMixed->Write();
560 hPPMixed->Write(); hNNMixed->Write();*/
eb63b883 561 gHistBalanceFunction->Write();
52daf7b2 562 if(listBFShuffled) gHistBalanceFunctionShuffled->Write();
563 if(listBFMixed) {
564 gHistBalanceFunctionMixed->Write();
565 gHistBalanceFunctionSubtracted->Write();
566 }
eb63b883 567 fOutput->Close();
6acdbcb2 568}
569
742af4bd 570//____________________________________________________________//
571void fitbalanceFunction(Int_t gCentrality = 1,
572 Double_t psiMin = -0.5, Double_t psiMax = 3.5,
573 Double_t ptTriggerMin = -1.,
574 Double_t ptTriggerMax = -1.,
575 Double_t ptAssociatedMin = -1.,
576 Double_t ptAssociatedMax = -1.,
d67eae53 577 TH2D *gHist,
578 Bool_t k2pMethod = kFALSE,
579 TString eventClass="EventPlane") {
48f521eb 580 //balancing charges: [1]*TMath::Exp(-0.5*TMath::Power(((x - [3])/[2]),2)-0.5*TMath::Power(((y - [5])/[4]),2))
581 //short range correlations: [6]*TMath::Exp(-0.5*TMath::Power(((x - [8])/[7]),2)-0.5*TMath::Power(((y - [10])/[9]),2))
742af4bd 582 cout<<"FITTING FUNCTION"<<endl;
583
48f521eb 584 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.);
742af4bd 585 gFitFunction->SetName("gFitFunction");
586
742af4bd 587 //Normalization
588 gFitFunction->SetParName(0,"N1");
48f521eb 589 gFitFunction->SetParameter(0,1.0);
590
591 //2D balance function
592 gFitFunction->SetParName(1,"N_{BF}");
593 gFitFunction->SetParameter(1,1.0);
594 gFitFunction->SetParLimits(1, 0., 100.);
595 gFitFunction->SetParName(2,"Sigma_{BF}(delta eta)");
596 gFitFunction->SetParameter(2,0.6);
597 gFitFunction->SetParLimits(2, 0., 1.);
598 gFitFunction->SetParName(3,"Mean_{BF}(delta eta)");
599 gFitFunction->SetParameter(3,0.0);
600 gFitFunction->SetParLimits(3, -0.2, 0.2);
601 gFitFunction->SetParName(4,"Sigma_{BF}(delta phi)");
602 gFitFunction->SetParameter(4,0.6);
603 gFitFunction->SetParLimits(4, 0., 1.);
604 gFitFunction->SetParName(5,"Mean_{BF}(delta phi)");
605 gFitFunction->SetParameter(5,0.0);
606 gFitFunction->SetParLimits(5, -0.2, 0.2);
607
608 //short range structure
609 gFitFunction->SetParName(6,"N_{SR}");
610 gFitFunction->SetParameter(6,5.0);
611 gFitFunction->SetParLimits(6, 0., 100.);
612 gFitFunction->SetParName(7,"Sigma_{SR}(delta eta)");
613 gFitFunction->SetParameter(7,0.01);
614 gFitFunction->SetParLimits(7, 0.0, 0.1);
615 gFitFunction->SetParName(8,"Mean_{SR}(delta eta)");
616 gFitFunction->SetParameter(8,0.0);
617 gFitFunction->SetParLimits(8, -0.01, 0.01);
618 gFitFunction->SetParName(9,"Sigma_{SR}(delta phi)");
619 gFitFunction->SetParameter(9,0.01);
620 gFitFunction->SetParLimits(9, 0.0, 0.1);
621 gFitFunction->SetParName(10,"Mean_{SR}(delta phi)");
622 gFitFunction->SetParameter(10,0.0);
623 gFitFunction->SetParLimits(10, -0.01, 0.01);
742af4bd 624
625
626 //Cloning the histogram
742af4bd 627 TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
628 gHistResidual->SetName("gHistResidual");
629 gHistResidual->Sumw2();
630
48f521eb 631 //Fitting the 2D bf
632 for(Int_t iAttempt = 0; iAttempt < 10; iAttempt++) {
633 gHist->Fit("gFitFunction","nm");
634 for(Int_t iParam = 0; iParam < 11; iParam++)
635 gFitFunction->SetParameter(iParam,gFitFunction->GetParameter(iParam));
742af4bd 636 }
48f521eb 637 cout<<"======================================================"<<endl;
638 cout<<"Fit chi2/ndf: "<<gFitFunction->GetChisquare()/gFitFunction->GetNDF()<<" - chi2: "<<gFitFunction->GetChisquare()<<" - ndf: "<<gFitFunction->GetNDF()<<endl;
639 cout<<"======================================================"<<endl;
640
641 //Getting the residual
642 gHistResidual->Add(gFitFunction,-1);
742af4bd 643
644 //Write to output file
d67eae53 645 TString newFileName = "balanceFunctionFit2D.";
646 if(eventClass == "Centrality"){
647 newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
648 newFileName += ".PsiAll.PttFrom";
649 }
650 else if(eventClass == "Multiplicity"){
651 newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
652 newFileName += ".PsiAll.PttFrom";
653 }
654 else{ // "EventPlane" (default)
655 newFileName += "Centrality";
656 newFileName += gCentrality; newFileName += ".Psi";
657 if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
658 else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
659 else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
660 else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
661 else newFileName += "All.PttFrom";
662 }
742af4bd 663 newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
664 newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
665 newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
666 newFileName += Form("%.1f",ptAssociatedMax);
d67eae53 667 if(k2pMethod) newFileName += "_2pMethod";
a9f20288 668
669 newFileName += "_";
670 newFileName += Form("%.1f",psiMin);
671 newFileName += "-";
672 newFileName += Form("%.1f",psiMax);
742af4bd 673 newFileName += ".root";
a9f20288 674
742af4bd 675 TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
676 gHist->Write();
742af4bd 677 gHistResidual->Write();
678 gFitFunction->Write();
679 newFile->Close();
742af4bd 680}
681
eb63b883 682//____________________________________________________________//
648f1a5a 683void drawBFPsi2D(const char* lhcPeriod = "LHC11h",
bd36d661 684 const char* gCentralityEstimator = "V0M",
685 Int_t gBit = 128,
686 const char* gEventPlaneEstimator = "VZERO",
648f1a5a 687 Int_t gCentrality = 1,
688 Bool_t kShowShuffled = kFALSE,
689 Bool_t kShowMixed = kFALSE,
eb63b883 690 Double_t psiMin = -0.5, Double_t psiMax = 0.5,
691 Double_t ptTriggerMin = -1.,
692 Double_t ptTriggerMax = -1.,
693 Double_t ptAssociatedMin = -1.,
bd36d661 694 Double_t ptAssociatedMax = -1.,
695 Bool_t k2pMethod = kTRUE) {
eb63b883 696 //Macro that draws the BF distributions for each centrality bin
697 //for reaction plane dependent analysis
698 //Author: Panos.Christakoglou@nikhef.nl
db7174c0 699 TGaxis::SetMaxDigits(3);
eb63b883 700
701 //Get the input file
bd36d661 702 TString filename = lhcPeriod;
703 filename += "/Centrality"; filename += gCentralityEstimator;
704 filename += "_Bit"; filename += gBit;
705 filename += "_"; filename += gEventPlaneEstimator;
706 filename +="/PttFrom";
648f1a5a 707 filename += Form("%.1f",ptTriggerMin); filename += "To";
708 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
709 filename += Form("%.1f",ptAssociatedMin); filename += "To";
710 filename += Form("%.1f",ptAssociatedMax);
711 filename += "/balanceFunction2D.Centrality";
eb63b883 712 filename += gCentrality; filename += ".Psi";
713 if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
714 else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
715 else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
648f1a5a 716 else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
717 else filename += "All.PttFrom";
718 filename += Form("%.1f",ptTriggerMin); filename += "To";
db7174c0 719 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
648f1a5a 720 filename += Form("%.1f",ptAssociatedMin); filename += "To";
bd36d661 721 filename += Form("%.1f",ptAssociatedMax);
722 if(k2pMethod) filename += "_2pMethod";
a9f20288 723
724 filename += "_";
725 filename += Form("%.1f",psiMin);
726 filename += "-";
727 filename += Form("%.1f",psiMax);
bd36d661 728 filename += ".root";
eb63b883 729
730 //Open the file
731 TFile *f = TFile::Open(filename.Data());
732 if((!f)||(!f->IsOpen())) {
733 Printf("The file %s is not found. Aborting...",filename);
734 return listBF;
6acdbcb2 735 }
eb63b883 736 //f->ls();
737
738 //Raw balance function
739 TH1D *gHistBalanceFunction = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunction"));
740 gHistBalanceFunction->SetStats(kFALSE);
741 gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
742 gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
743 gHistBalanceFunction->GetZaxis()->SetNdivisions(10);
744 gHistBalanceFunction->GetXaxis()->SetTitleOffset(1.3);
745 gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
746 gHistBalanceFunction->GetZaxis()->SetTitleOffset(1.3);
747 gHistBalanceFunction->GetXaxis()->SetTitle("#Delta #eta");
27634c13 748 gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (rad)");
15dd45a0 749 gHistBalanceFunction->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
eb63b883 750
751 //Shuffled balance function
648f1a5a 752 if(kShowShuffled) {
753 TH1D *gHistBalanceFunctionShuffled = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionShuffled"));
754 gHistBalanceFunctionShuffled->SetStats(kFALSE);
755 gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
756 gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
757 gHistBalanceFunctionShuffled->GetZaxis()->SetNdivisions(10);
758 gHistBalanceFunctionShuffled->GetXaxis()->SetTitleOffset(1.3);
759 gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
760 gHistBalanceFunctionShuffled->GetZaxis()->SetTitleOffset(1.3);
761 gHistBalanceFunctionShuffled->GetXaxis()->SetTitle("#Delta #eta");
27634c13 762 gHistBalanceFunctionShuffled->GetYaxis()->SetTitle("#Delta #varphi (rad)");
648f1a5a 763 gHistBalanceFunctionShuffled->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
764 }
eb63b883 765
766 //Mixed balance function
648f1a5a 767 if(kShowMixed) {
768 TH1D *gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionMixed"));
769 gHistBalanceFunctionMixed->SetStats(kFALSE);
770 gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
771 gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
772 gHistBalanceFunctionMixed->GetZaxis()->SetNdivisions(10);
773 gHistBalanceFunctionMixed->GetXaxis()->SetTitleOffset(1.3);
774 gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
775 gHistBalanceFunctionMixed->GetZaxis()->SetTitleOffset(1.3);
776 gHistBalanceFunctionMixed->GetXaxis()->SetTitle("#Delta #eta");
27634c13 777 gHistBalanceFunctionMixed->GetYaxis()->SetTitle("#Delta #varphi (rad)");
648f1a5a 778 gHistBalanceFunctionMixed->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
779 }
eb63b883 780
781 //Subtracted balance function
648f1a5a 782 if(kShowMixed) {
783 TH1D *gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionSubtracted"));
784 gHistBalanceFunctionSubtracted->SetStats(kFALSE);
785 gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
786 gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
787 gHistBalanceFunctionSubtracted->GetZaxis()->SetNdivisions(10);
788 gHistBalanceFunctionSubtracted->GetXaxis()->SetTitleOffset(1.3);
789 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
790 gHistBalanceFunctionSubtracted->GetZaxis()->SetTitleOffset(1.3);
791 gHistBalanceFunctionSubtracted->GetXaxis()->SetTitle("#Delta #eta");
27634c13 792 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitle("#Delta #varphi (rad)");
648f1a5a 793 gHistBalanceFunctionSubtracted->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
794 }
795
eb63b883 796 TString pngName;
797
798 TString centralityLatex = "Centrality: ";
799 centralityLatex += centralityArray[gCentrality-1];
800 centralityLatex += "%";
801
802 TString psiLatex;
803 if((psiMin == -0.5)&&(psiMax == 0.5))
804 psiLatex = " -7.5^{o} < #varphi - #Psi_{2} < 7.5^{o}";
805 else if((psiMin == 0.5)&&(psiMax == 1.5))
806 psiLatex = " 37.5^{o} < #varphi - #Psi_{2} < 52.5^{o}";
807 else if((psiMin == 1.5)&&(psiMax == 2.5))
808 psiLatex = " 82.5^{o} < #varphi - #Psi_{2} < 97.5^{o}";
809 else
810 psiLatex = " 0^{o} < #varphi - #Psi_{2} < 180^{o}";
811
812 TString pttLatex = Form("%.1f",ptTriggerMin);
813 pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
814 pttLatex += " GeV/c";
815
816 TString ptaLatex = Form("%.1f",ptAssociatedMin);
817 ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
818 ptaLatex += " GeV/c";
819
820 TLatex *latexInfo1 = new TLatex();
821 latexInfo1->SetNDC();
15dd45a0 822 latexInfo1->SetTextSize(0.045);
eb63b883 823 latexInfo1->SetTextColor(1);
824
825 //Draw the results
826 TCanvas *c1 = new TCanvas("c1","Raw balance function 2D",0,0,600,500);
827 c1->SetFillColor(10); c1->SetHighLightColor(10);
828 c1->SetLeftMargin(0.17); c1->SetTopMargin(0.05);
bd36d661 829 gHistBalanceFunction->SetTitle("");
eb63b883 830 gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.4);
831 gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
53ff9638 832 //gHistBalanceFunction->GetXaxis()->SetRangeUser(-1.4,1.4);
eb63b883 833 gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
27634c13 834 gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (rad)");
eb63b883 835 gHistBalanceFunction->DrawCopy("lego2");
5de9ad1a 836 gPad->SetTheta(30); // default is 30
837 gPad->SetPhi(-60); // default is 30
838 gPad->Update();
eb63b883 839
15dd45a0 840 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
841 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
842 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
843 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
eb63b883 844
bd36d661 845 TString pngName = "BalanceFunction2D.";
846 pngName += "Centrality";
847 pngName += gCentrality;
848 if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
849 else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
850 else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
851 else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.PttFrom";
852 else pngName += "All.PttFrom";
853 pngName += Form("%.1f",ptTriggerMin); pngName += "To";
854 pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
855 pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
856 pngName += Form("%.1f",ptAssociatedMax);
857 if(k2pMethod) pngName += "_2pMethod";
858 pngName += ".png";
859
860 c1->SaveAs(pngName.Data());
861
648f1a5a 862 if(kShowShuffled) {
863 TCanvas *c2 = new TCanvas("c2","Shuffled balance function 2D",100,100,600,500);
864 c2->SetFillColor(10); c2->SetHighLightColor(10);
865 c2->SetLeftMargin(0.17); c2->SetTopMargin(0.05);
866 gHistBalanceFunctionShuffled->SetTitle("Shuffled events");
867 gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.4);
868 gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
869 gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
870 gHistBalanceFunctionShuffled->DrawCopy("lego2");
5de9ad1a 871 gPad->SetTheta(30); // default is 30
872 gPad->SetPhi(-60); // default is 30
873 gPad->Update();
874
648f1a5a 875 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
876 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
877 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
878 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
879 }
eb63b883 880
648f1a5a 881 if(kShowMixed) {
882 TCanvas *c3 = new TCanvas("c3","Mixed balance function 2D",200,200,600,500);
883 c3->SetFillColor(10); c3->SetHighLightColor(10);
884 c3->SetLeftMargin(0.17); c3->SetTopMargin(0.05);
885 gHistBalanceFunctionMixed->SetTitle("Mixed events");
886 gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.4);
887 gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
888 gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
889 gHistBalanceFunctionMixed->DrawCopy("lego2");
5de9ad1a 890 gPad->SetTheta(30); // default is 30
891 gPad->SetPhi(-60); // default is 30
892 gPad->Update();
893
648f1a5a 894 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
895 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
896 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
897 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
898 }
eb63b883 899
648f1a5a 900 if(kShowMixed) {
901 TCanvas *c4 = new TCanvas("c4","Subtracted balance function 2D",300,300,600,500);
902 c4->SetFillColor(10); c4->SetHighLightColor(10);
903 c4->SetLeftMargin(0.17); c4->SetTopMargin(0.05);
904 gHistBalanceFunctionSubtracted->SetTitle("Subtracted balance function");
905 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.4);
906 gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
907 gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
908 gHistBalanceFunctionSubtracted->DrawCopy("lego2");
5de9ad1a 909 gPad->SetTheta(30); // default is 30
910 gPad->SetPhi(-60); // default is 30
911 gPad->Update();
912
648f1a5a 913 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
914 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
915 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
916 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
917 }
6acdbcb2 918}
17c69b4e 919
920//____________________________________________________________//
79dc837f 921void drawProjections(TH2D *gHistBalanceFunction2D = 0x0,
922 Bool_t kProjectInEta = kFALSE,
17c69b4e 923 Int_t binMin = 1,
924 Int_t binMax = 80,
925 Int_t gCentrality = 1,
926 Double_t psiMin = -0.5,
927 Double_t psiMax = 3.5,
17c69b4e 928 Double_t ptTriggerMin = -1.,
929 Double_t ptTriggerMax = -1.,
930 Double_t ptAssociatedMin = -1.,
931 Double_t ptAssociatedMax = -1.,
17c69b4e 932 Bool_t k2pMethod = kTRUE,
933 TString eventClass = "Centrality",
1822002d 934 Bool_t bRootMoments = kFALSE,
40a8b898 935 Bool_t kUseZYAM = kFALSE,
7827f01e 936 Bool_t bReduceRangeForMoments = kFALSE,
937 Bool_t bFWHM = kFALSE) {
438dab2c 938 //Macro that draws the balance functions PROJECTIONS
17c69b4e 939 //for each centrality bin for the different pT of trigger and
940 //associated particles
941 TGaxis::SetMaxDigits(3);
942
943 //first we need some libraries
a9f20288 944 gSystem->Load("libTree");
945 gSystem->Load("libGeom");
946 gSystem->Load("libVMC");
947 gSystem->Load("libXMLIO");
948 gSystem->Load("libPhysics");
949
950 gSystem->Load("libSTEERBase");
951 gSystem->Load("libESD");
952 gSystem->Load("libAOD");
953
17c69b4e 954 gSystem->Load("libANALYSIS.so");
955 gSystem->Load("libANALYSISalice.so");
956 gSystem->Load("libEventMixing.so");
957 gSystem->Load("libCORRFW.so");
958 gSystem->Load("libPWGTools.so");
959 gSystem->Load("libPWGCFebye.so");
960
438dab2c 961 gStyle->SetOptStat(0);
962
17c69b4e 963 //Get the input file
964 TString filename = "balanceFunction2D.";
965 if(eventClass == "Centrality"){
966 filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
967 filename += ".PsiAll.PttFrom";
968 }
969 else if(eventClass == "Multiplicity"){
970 filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
971 filename += ".PsiAll.PttFrom";
972 }
973 else{ // "EventPlane" (default)
974 filename += "Centrality";
975 filename += gCentrality; filename += ".Psi";
976 if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
977 else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
978 else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
979 else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
980 else filename += "All.PttFrom";
981 }
982 filename += Form("%.1f",ptTriggerMin); filename += "To";
983 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
984 filename += Form("%.1f",ptAssociatedMin); filename += "To";
985 filename += Form("%.1f",ptAssociatedMax);
986 if(k2pMethod) filename += "_2pMethod";
a9f20288 987
a2da950e 988 filename += "_";
989 filename += Form("%.1f",psiMin);
990 filename += "-";
991 filename += Form("%.1f",psiMax);
17c69b4e 992 filename += ".root";
993
994 //Open the file
79dc837f 995 TFile *f = 0x0;
996 if(!gHistBalanceFunction2D) {
997 TFile::Open(filename.Data());
998 if((!f)||(!f->IsOpen())) {
999 Printf("The file %s is not found. Aborting...",filename);
1000 return listBF;
1001 }
1002 //f->ls();
17c69b4e 1003 }
79dc837f 1004
17c69b4e 1005 //Latex
1006 TString centralityLatex = "Centrality: ";
1007 centralityLatex += centralityArray[gCentrality-1];
1008 centralityLatex += "%";
1009
1010 TString psiLatex;
1011 if((psiMin == -0.5)&&(psiMax == 0.5))
1012 psiLatex = " -7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o}";
1013 else if((psiMin == 0.5)&&(psiMax == 1.5))
1014 psiLatex = " 37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o}";
1015 else if((psiMin == 1.5)&&(psiMax == 2.5))
1016 psiLatex = " 82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o}";
1017 else
1018 psiLatex = " 0^{o} < #varphi^{t} - #Psi_{2} < 180^{o}";
1019
1020 TString pttLatex = Form("%.1f",ptTriggerMin);
1021 pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
1022 pttLatex += " GeV/c";
1023
1024 TString ptaLatex = Form("%.1f",ptAssociatedMin);
1025 ptaLatex += " < p_{T,assoc} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1026 ptaLatex += " GeV/c";
1027
1028 TLatex *latexInfo1 = new TLatex();
1029 latexInfo1->SetNDC();
1030 latexInfo1->SetTextSize(0.045);
1031 latexInfo1->SetTextColor(1);
1032
17c69b4e 1033
1034 //============================================================//
1035 //Get subtracted and mixed balance function
79dc837f 1036 TH2D *gHistBalanceFunctionSubtracted2D = 0x0;
1037 TH2D *gHistBalanceFunctionMixed2D = 0x0;
1038 if(!gHistBalanceFunction2D) {
1039 gHistBalanceFunctionSubtracted2D = (TH2D*)f->Get("gHistBalanceFunctionSubtracted");
1040 gHistBalanceFunctionMixed2D = (TH2D*)f->Get("gHistBalanceFunctionMixed");
1041 }
1042 else {
1043 gHistBalanceFunctionSubtracted2D = dynamic_cast<TH2D*>(gHistBalanceFunction2D->Clone());
1044 gHistBalanceFunctionMixed2D = dynamic_cast<TH2D*>(gHistBalanceFunction2D->Clone());
1045 }
17c69b4e 1046
1047 TH1D *gHistBalanceFunctionSubtracted = NULL;
1048 TH1D *gHistBalanceFunctionMixed = NULL;
1049
a9f20288 1050 TH1D *gHistBalanceFunctionSubtracted_scale = NULL;
1051 TH1D *gHistBalanceFunctionMixed_scale = NULL;
1052
17c69b4e 1053 if(kProjectInEta){
86044267 1054 gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionX("gHistBalanceFunctionSubtractedEta",binMin,binMax));
a9f20288 1055 gHistBalanceFunctionSubtracted->Scale(gHistBalanceFunctionSubtracted2D->GetYaxis()->GetBinWidth(1)); // to remove normalization to phi bin width
86044267 1056 gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionX("gHistBalanceFunctionMixedEta",binMin,binMax));
a9f20288 1057 gHistBalanceFunctionMixed->Scale(gHistBalanceFunctionMixed2D->GetYaxis()->GetBinWidth(1)); // to remove normalization to phi bin width
17c69b4e 1058 gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#eta)");
1059 gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#eta)");
1060 }
1061 else{
86044267 1062 gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionY("gHistBalanceFunctionSubtractedPhi",binMin,binMax));
a9f20288 1063 gHistBalanceFunctionSubtracted->Scale(gHistBalanceFunctionSubtracted2D->GetXaxis()->GetBinWidth(1)); // to remove normalization to eta bin width
86044267 1064 gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionY("gHistBalanceFunctionMixedPhi",binMin,binMax));
a9f20288 1065 gHistBalanceFunctionMixed->Scale(gHistBalanceFunctionMixed2D->GetXaxis()->GetBinWidth(1)); // to remove normalization to eta bin width
17c69b4e 1066 gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#varphi)");
1067 gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#varphi)");
1068 }
1069
1070 gHistBalanceFunctionSubtracted->SetMarkerStyle(20);
1071 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
17c69b4e 1072
1073 gHistBalanceFunctionMixed->SetMarkerStyle(25);
17c69b4e 1074
1075 TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
1076 c1->SetFillColor(10);
1077 c1->SetHighLightColor(10);
1078 c1->SetLeftMargin(0.15);
1079 gHistBalanceFunctionSubtracted->DrawCopy("E");
a9f20288 1080 gHistBalanceFunctionMixed->DrawCopy("ESAME");
17c69b4e 1081
1082 legend = new TLegend(0.18,0.62,0.45,0.82,"","brNDC");
1083 legend->SetTextSize(0.045);
1084 legend->SetTextFont(42);
1085 legend->SetBorderSize(0);
1086 legend->SetFillStyle(0);
1087 legend->SetFillColor(10);
1088 legend->SetMargin(0.25);
1089 legend->SetShadowColor(10);
1090 legend->AddEntry(gHistBalanceFunctionSubtracted,"Data","lp");
1091 legend->AddEntry(gHistBalanceFunctionMixed,"Mixed data","lp");
1092 legend->Draw();
1093
17c69b4e 1094
17c69b4e 1095 TString meanLatex, rmsLatex, skewnessLatex, kurtosisLatex;
1096
1097 if(bRootMoments){
732acb5b 1098
1099 // need to restrict to near side in case of dphi
1100 if(!kProjectInEta) gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-TMath::Pi()/2,TMath::Pi()/2);
1101
53ff9638 1102 if(bReduceRangeForMoments){
40a8b898 1103
53ff9638 1104 // safety check:
40a8b898 1105 // default (for 2<pT,assoc<3<pT,trig<4)
53ff9638 1106 Double_t rangeReduced = 1.2;
40a8b898
MW		 1107 //for 3<pT,assoc<8<pT,trig<15
1108 if(ptTriggerMax>10){
53ff9638 1109 rangeReduced = 0.35;
1110 }
1111
1112 // new define range by fit!
1113 gHistBalanceFunctionSubtracted->Fit("gaus","","");
1114 Double_t sigmaGaus = gHistBalanceFunctionSubtracted->GetFunction("gaus")->GetParameter(2);
1115
1116 // if safety check OK, set rangeReduced to 3 sigma of gauss fit
1117 if(3*sigmaGaus > rangeReduced){
1118 cout<<"RANGE REDUCE ERROR: "<<3*sigmaGaus<<" > "<<rangeReduced<<endl;
1119 }
1120 else{
1121 rangeReduced = 3*sigmaGaus;
40a8b898
MW		 1122 }
1123
1124 cout<<"Use reduced range = "<<rangeReduced<<endl;
1125 gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-rangeReduced,rangeReduced);
1126 }
1127
17c69b4e 1128 meanLatex = "#mu = ";
1129 meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMean());
1130 meanLatex += " #pm ";
1131 meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMeanError());
1132
1133 rmsLatex = "#sigma = ";
1134 rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMS());
1135 rmsLatex += " #pm ";
1136 rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMSError());
1137
1138 skewnessLatex = "S = ";
1139 skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(1));
1140 skewnessLatex += " #pm ";
1141 skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(11));
1142
1143 kurtosisLatex = "K = ";
1144 kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(1));
1145 kurtosisLatex += " #pm ";
1146 kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(11));
a2da950e 1147
732acb5b 1148 Printf("ROOT Mean: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetMean(),gHistBalanceFunctionSubtracted->GetMeanError());
1149 Printf("ROOT RMS: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetRMS(),gHistBalanceFunctionSubtracted->GetRMSError());
1150 Printf("ROOT Skeweness: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetSkewness(1),gHistBalanceFunctionSubtracted->GetSkewness(11));
1151 Printf("ROOT Kurtosis: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetKurtosis(1),gHistBalanceFunctionSubtracted->GetKurtosis(11));
438dab2c 1152
732acb5b 1153
438dab2c 1154 // store in txt files
438dab2c 1155 TString meanFileName = filename;
79dc837f 1156 if(kProjectInEta)
23a0fc15 1157 meanFileName= "deltaEtaProjection_Mean.txt";
1158 //meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
79dc837f 1159 else
23a0fc15 1160 meanFileName = "deltaPhiProjection_Mean.txt";
1161 //meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
79dc837f 1162 ofstream fileMean(meanFileName.Data(),ios::app);
438dab2c 1163 fileMean << " " << gHistBalanceFunctionSubtracted->GetMean() << " " <<gHistBalanceFunctionSubtracted->GetMeanError()<<endl;
1164 fileMean.close();
1165
1166 TString rmsFileName = filename;
79dc837f 1167 if(kProjectInEta)
23a0fc15 1168 rmsFileName = "deltaEtaProjection_Rms.txt";
1169 //rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
79dc837f 1170 else
732acb5b 1171 rmsFileName = "deltaPhiProjection_Rms.txt";
23a0fc15 1172 //rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
79dc837f 1173 ofstream fileRms(rmsFileName.Data(),ios::app);
438dab2c 1174 fileRms << " " << gHistBalanceFunctionSubtracted->GetRMS() << " " <<gHistBalanceFunctionSubtracted->GetRMSError()<<endl;
1175 fileRms.close();
1176
1177 TString skewnessFileName = filename;
79dc837f 1178 if(kProjectInEta)
23a0fc15 1179 skewnessFileName = "deltaEtaProjection_Skewness.txt";
1180 //skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
79dc837f 1181 else
23a0fc15 1182 skewnessFileName = "deltaPhiProjection_Skewness.txt";
1183 //skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
79dc837f 1184 ofstream fileSkewness(skewnessFileName.Data(),ios::app);
438dab2c 1185 fileSkewness << " " << gHistBalanceFunctionSubtracted->GetSkewness(1) << " " <<gHistBalanceFunctionSubtracted->GetSkewness(11)<<endl;
1186 fileSkewness.close();
1187
1188 TString kurtosisFileName = filename;
79dc837f 1189 if(kProjectInEta)
23a0fc15 1190 kurtosisFileName = "deltaEtaProjection_Kurtosis.txt";
1191 //kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
79dc837f 1192 else
23a0fc15 1193 kurtosisFileName = "deltaPhiProjection_Kurtosis.txt";
1194 //kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
79dc837f 1195 ofstream fileKurtosis(kurtosisFileName.Data(),ios::app);
438dab2c 1196 fileKurtosis << " " << gHistBalanceFunctionSubtracted->GetKurtosis(1) << " " <<gHistBalanceFunctionSubtracted->GetKurtosis(11)<<endl;
1197 fileKurtosis.close();
732acb5b 1198
1199 if(!kProjectInEta) gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-TMath::Pi()/2,3.*TMath::Pi()/2);
40a8b898 1200 else gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-1.6,1.6);
17c69b4e 1201 }
1202 // calculate the moments by hand
1203 else{
1204
1205 Double_t meanAnalytical, meanAnalyticalError;
1206 Double_t sigmaAnalytical, sigmaAnalyticalError;
1207 Double_t skewnessAnalytical, skewnessAnalyticalError;
1208 Double_t kurtosisAnalytical, kurtosisAnalyticalError;
1822002d 1209
17c69b4e 1210 Int_t gDeltaEtaPhi = 2;
1211 if(kProjectInEta) gDeltaEtaPhi = 1;
1212
1213 AliBalancePsi *bHelper = new AliBalancePsi;
1822002d 1214 bHelper->GetMomentsAnalytical(gDeltaEtaPhi,gHistBalanceFunctionSubtracted,kUseZYAM,meanAnalytical,meanAnalyticalError,sigmaAnalytical,sigmaAnalyticalError,skewnessAnalytical,skewnessAnalyticalError,kurtosisAnalytical,kurtosisAnalyticalError);
17c69b4e 1215
1216 meanLatex = "#mu = ";
1217 meanLatex += Form("%.3f",meanAnalytical);
1218 meanLatex += " #pm ";
1219 meanLatex += Form("%.3f",meanAnalyticalError);
1220
1221 rmsLatex = "#sigma = ";
1222 rmsLatex += Form("%.3f",sigmaAnalytical);
1223 rmsLatex += " #pm ";
1224 rmsLatex += Form("%.3f",sigmaAnalyticalError);
1225
1226 skewnessLatex = "S = ";
1227 skewnessLatex += Form("%.3f",skewnessAnalytical);
1228 skewnessLatex += " #pm ";
1229 skewnessLatex += Form("%.3f",skewnessAnalyticalError);
1230
1231 kurtosisLatex = "K = ";
1232 kurtosisLatex += Form("%.3f",kurtosisAnalytical);
1233 kurtosisLatex += " #pm ";
1234 kurtosisLatex += Form("%.3f",kurtosisAnalyticalError);
1235 Printf("Mean: %lf - Error: %lf",meanAnalytical, meanAnalyticalError);
1236 Printf("Sigma: %lf - Error: %lf",sigmaAnalytical, sigmaAnalyticalError);
1237 Printf("Skeweness: %lf - Error: %lf",skewnessAnalytical, skewnessAnalyticalError);
1238 Printf("Kurtosis: %lf - Error: %lf",kurtosisAnalytical, kurtosisAnalyticalError);
438dab2c 1239
1240 // store in txt files
1241 TString meanFileName = filename;
79dc837f 1242 if(kProjectInEta)
23a0fc15 1243 meanFileName = "deltaEtaProjection_Mean.txt";
1244 //meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
79dc837f 1245 else
23a0fc15 1246 meanFileName = "deltaPhiProjection_Mean.txt";
1247 //meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
79dc837f 1248 ofstream fileMean(meanFileName.Data(),ios::app);
438dab2c 1249 fileMean << " " << meanAnalytical << " " <<meanAnalyticalError<<endl;
1250 fileMean.close();
1251
1252 TString rmsFileName = filename;
79dc837f 1253 if(kProjectInEta)
23a0fc15 1254 rmsFileName = "deltaEtaProjection_Rms.txt";
1255 //rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
79dc837f 1256 else
23a0fc15 1257 rmsFileName = "deltaPhiProjection_Rms.txt";
1258//rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
79dc837f 1259 ofstream fileRms(rmsFileName.Data(),ios::app);
438dab2c 1260 fileRms << " " << sigmaAnalytical << " " <<sigmaAnalyticalError<<endl;
1261 fileRms.close();
1262
1263 TString skewnessFileName = filename;
79dc837f 1264 if(kProjectInEta)
23a0fc15 1265 skewnessFileName = "deltaEtaProjection_Skewness.txt";
1266 //skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
79dc837f 1267 else
23a0fc15 1268 skewnessFileName = "deltaPhiProjection_Skewness.txt";
1269 //skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
79dc837f 1270 ofstream fileSkewness(skewnessFileName.Data(),ios::app);
438dab2c 1271 fileSkewness << " " << skewnessAnalytical << " " <<skewnessAnalyticalError<<endl;
1272 fileSkewness.close();
1273
1274 TString kurtosisFileName = filename;
79dc837f 1275 if(kProjectInEta)
23a0fc15 1276 kurtosisFileName = "deltaEtaProjection_Kurtosis.txt";
1277 //kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
79dc837f 1278 else
23a0fc15 1279 kurtosisFileName = "deltaPhiProjection_Kurtosis.txt";
1280 //kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
79dc837f 1281 ofstream fileKurtosis(kurtosisFileName.Data(),ios::app);
438dab2c 1282 fileKurtosis << " " << kurtosisAnalytical << " " <<kurtosisAnalyticalError<<endl;
1283 fileKurtosis.close();
17c69b4e 1284 }
1285
3de51a46 1286 // Weighted mean as calculated for 1D analysis
1287 Double_t weightedMean, weightedMeanError;
40a8b898
MW		 1288 if(bReduceRangeForMoments){
1289 GetWeightedMean1D(gHistBalanceFunctionSubtracted,kProjectInEta,1,-1,weightedMean,weightedMeanError,rangeReduced);
1290 }
1291 else{
1292 GetWeightedMean1D(gHistBalanceFunctionSubtracted,kProjectInEta,1,-1,weightedMean,weightedMeanError,-1.);
1293 }
3de51a46 1294 Printf("Weighted Mean: %lf - Error: %lf",weightedMean, weightedMeanError);
1295
1296 // store in txt files
1297 TString weightedMeanFileName = filename;
1298 if(kProjectInEta)
1299 weightedMeanFileName = "deltaEtaProjection_WeightedMean.txt";
1300 //weightedMeanFileName.ReplaceAll(".root","_DeltaEtaProjection_WeightedMean.txt");
1301 else
1302 weightedMeanFileName = "deltaPhiProjection_WeightedMean.txt";
1303 //weightedMeanFileName.ReplaceAll(".root","_DeltaPhiProjection_WeightedMean.txt");
1304 ofstream fileWeightedMean(weightedMeanFileName.Data(),ios::app);
1305 fileWeightedMean << " " << weightedMean << " " <<weightedMeanError<<endl;
1306 fileWeightedMean.close();
438dab2c 1307
17c69b4e 1308 TCanvas *c2 = new TCanvas("c2","",600,0,600,500);
1309 c2->SetFillColor(10);
1310 c2->SetHighLightColor(10);
1311 c2->SetLeftMargin(0.15);
1312 gHistBalanceFunctionSubtracted->DrawCopy("E");
7827f01e 1313
1314 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++//
1315 Double_t fwhm_spline = 0.;
1316 Double_t fwhmError = 0.;
1317 if (bFWHM){
1318 AliBalancePsi *bHelper_1 = new AliBalancePsi;
1319 bHelper_1->GetFWHM(kProjectInEta,gHistBalanceFunctionSubtracted,fwhm_spline,fwhmError);
1320 Printf("FWHM: %lf - Error: %lf",fwhm_spline, fwhmError);
1321
1322 //store and in txt files FWHM
1323 TString sigmaFileName = filename;
1324 if(kProjectInEta) {
1325 sigmaFileName = "deltaEtaProjection_FWHM.txt";
1326 }
1327 else{
1328 sigmaFileName = "deltaPhiProjection_FWHM.txt";
1329 }
1330 ofstream fileSigma(sigmaFileName.Data(),ios::app);
1331 fileSigma << " " << fwhm_spline << " " <<fwhmError<<endl;
1332 fileSigma.close();
1333
1334 gHistBalanceFunctionSubtracted->SetLineColor(2);
1335 gHistBalanceFunctionSubtracted->SetMarkerColor(2);
1336 gHistBalanceFunctionSubtracted->Draw("SAME");
1337 }
1338 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++//
438dab2c 1339
17c69b4e 1340 TLatex *latex = new TLatex();
1341 latex->SetNDC();
1342 latex->SetTextSize(0.035);
1343 latex->SetTextColor(1);
1344 latex->DrawLatex(0.64,0.85,meanLatex.Data());
1345 latex->DrawLatex(0.64,0.81,rmsLatex.Data());
1346 latex->DrawLatex(0.64,0.77,skewnessLatex.Data());
1347 latex->DrawLatex(0.64,0.73,kurtosisLatex.Data());
1348
438dab2c 1349 TString pngName = filename;
1350 if(kProjectInEta) pngName.ReplaceAll(".root","_DeltaEtaProjection.png");
1351 else pngName.ReplaceAll(".root","_DeltaPhiProjection.png");
1352
1353 c2->SaveAs(pngName.Data());
1354
17c69b4e 1355 TString outFileName = filename;
1356 if(kProjectInEta) outFileName.ReplaceAll(".root","_DeltaEtaProjection.root");
1357 else outFileName.ReplaceAll(".root","_DeltaPhiProjection.root");
1358 TFile *fProjection = TFile::Open(outFileName.Data(),"recreate");
1359 gHistBalanceFunctionSubtracted->Write();
1360 gHistBalanceFunctionMixed->Write();
1361 fProjection->Close();
1362}
79dc837f 1363
1364//____________________________________________________________//
1365void drawBFPsi2DFromCorrelationFunctions(const char* lhcPeriod = "LHC10h",
1366 Int_t gTrainNumber = 64,
1367 const char* gCentralityEstimator = "V0M",
1368 Int_t gBit = 128,
1369 const char* gEventPlaneEstimator = "VZERO",
1370 Int_t gCentrality = 1,
1371 Double_t psiMin = -0.5, Double_t psiMax = 3.5,
1372 Double_t vertexZMin = -10.,
1373 Double_t vertexZMax = 10.,
1374 Double_t ptTriggerMin = -1.,
1375 Double_t ptTriggerMax = -1.,
1376 Double_t ptAssociatedMin = -1.,
86044267 1377 Double_t ptAssociatedMax = -1.,
732acb5b 1378 TString eventClass = "Multiplicity",
a91351d9 1379 Bool_t bRootMoments = kFALSE,
40a8b898 1380 Bool_t bFullPhiForEtaProjection = kFALSE,
7827f01e 1381 Bool_t bReduceRangeForMoments = kFALSE,
1382 Bool_t bFWHM = kFALSE) {
79dc837f 1383 //Macro that draws the BF distributions for each centrality bin
1384 //for reaction plane dependent analysis
1385 //Author: Panos.Christakoglou@nikhef.nl
1386 TGaxis::SetMaxDigits(3);
1387 gStyle->SetPalette(55,0);
1388
1389 //Get the input file
1390 TString filename = lhcPeriod;
86044267 1391 if(lhcPeriod != ""){
1392 //filename += "/Train"; filename += gTrainNumber;
1393 filename +="/PttFrom";
1394 filename += Form("%.1f",ptTriggerMin); filename += "To";
1395 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1396 filename += Form("%.1f",ptAssociatedMin); filename += "To";
1397 filename += Form("%.1f",ptAssociatedMax);
1398 filename += "/correlationFunction.";
1399 }
1400 else{
1401 filename += "correlationFunction.";
1402 }
1403 if(eventClass == "Centrality"){
1404 filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1405 filename += ".PsiAll.PttFrom";
1406 }
1407 else if(eventClass == "Multiplicity"){
1408 filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1409 filename += ".PsiAll.PttFrom";
1410 }
1411 else{ // "EventPlane" (default)
1412 filename += "Centrality";
1413 filename += gCentrality; filename += ".Psi";
1414 if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
1415 else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
1416 else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
1417 else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
1418 else filename += "All.PttFrom";
1419 }
79dc837f 1420 filename += Form("%.1f",ptTriggerMin); filename += "To";
1421 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1422 filename += Form("%.1f",ptAssociatedMin); filename += "To";
1423 filename += Form("%.1f",ptAssociatedMax);
1424 filename += "_";
1425 filename += Form("%.1f",psiMin);
1426 filename += "-";
1427 filename += Form("%.1f",psiMax);
1428 filename += ".root";
1429
1430 //Open the file
1431 TFile *f = TFile::Open(filename.Data());
1432 if((!f)||(!f->IsOpen())) {
1433 Printf("The file %s is not found. Aborting...",filename);
1434 return listBF;
1435 }
1436 //f->ls();
1437
1438 TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
1439 if(!gHistPN) return;
1440 TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
1441 if(!gHistNP) return;
1442 TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
1443 if(!gHistPP) return;
1444 TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
1445 if(!gHistNN) return;
1446
53ff9638 1447 // in order to get unzoomed (in older versions used smaller user ranger)
1448 Int_t binMinX = gHistPN->GetXaxis()->FindBin(gHistPN->GetXaxis()->GetXmin()+0.00001);
1449 Int_t binMaxX = gHistPN->GetXaxis()->FindBin(gHistPN->GetXaxis()->GetXmax()-0.00001);
1450 Int_t binMinY = gHistPN->GetYaxis()->FindBin(gHistPN->GetYaxis()->GetXmin()+0.00001);
1451 Int_t binMaxY = gHistPN->GetYaxis()->FindBin(gHistPN->GetYaxis()->GetXmax()-0.00001);
1452 gHistPN->GetXaxis()->SetRange(binMinX,binMaxX); gHistPN->GetYaxis()->SetRange(binMinY,binMaxY);
1453 gHistNP->GetXaxis()->SetRange(binMinX,binMaxX); gHistNP->GetYaxis()->SetRange(binMinY,binMaxY);
1454 gHistPP->GetXaxis()->SetRange(binMinX,binMaxX); gHistPP->GetYaxis()->SetRange(binMinY,binMaxY);
1455 gHistNN->GetXaxis()->SetRange(binMinX,binMaxX); gHistNN->GetYaxis()->SetRange(binMinY,binMaxY);
1456
79dc837f 1457 gHistPN->Sumw2();
1458 gHistPP->Sumw2();
1459 gHistPN->Add(gHistPP,-1);
1460 gHistNP->Sumw2();
1461 gHistNN->Sumw2();
1462 gHistNP->Add(gHistNN,-1);
1463 gHistPN->Add(gHistNP);
3aa5487c 1464 //gHistPN->Scale(0.5);//not needed anymore, since pT,assoc < pT,trig in all pT bins
79dc837f 1465 TH2D *gHistBalanceFunction2D = dynamic_cast<TH2D *>(gHistPN->Clone());
1466 gHistBalanceFunction2D->SetStats(kFALSE);
1467 gHistBalanceFunction2D->GetXaxis()->SetTitle("#Delta#eta");
1468 gHistBalanceFunction2D->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1469 gHistBalanceFunction2D->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
1470
1471 //Draw the results
1472 TCanvas *c0 = new TCanvas("c0","Balance function 2D",0,0,600,500);
1473 c0->SetFillColor(10); c0->SetHighLightColor(10);
1474 c0->SetLeftMargin(0.17); c0->SetTopMargin(0.05);
1475 gHistBalanceFunction2D->SetTitle("");
1476 gHistBalanceFunction2D->GetZaxis()->SetTitleOffset(1.4);
1477 gHistBalanceFunction2D->GetZaxis()->SetNdivisions(10);
1478 gHistBalanceFunction2D->GetYaxis()->SetTitleOffset(1.4);
1479 gHistBalanceFunction2D->GetYaxis()->SetNdivisions(10);
53ff9638 1480 //gHistBalanceFunction2D->GetXaxis()->SetRangeUser(-1.4,1.4);
79dc837f 1481 gHistBalanceFunction2D->GetXaxis()->SetNdivisions(10);
1482 gHistBalanceFunction2D->DrawCopy("lego2");
1483 gPad->SetTheta(30); // default is 30
1484 gPad->SetPhi(-60); // default is 30
1485 gPad->Update();
1486
1487 TString multLatex = Form("Multiplicity: %.1f - %.1f",psiMin,psiMax);
1488
1489 TString pttLatex = Form("%.1f",ptTriggerMin);
1490 pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
1491 pttLatex += " GeV/c";
1492
1493 TString ptaLatex = Form("%.1f",ptAssociatedMin);
1494 ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1495 ptaLatex += " GeV/c";
1496
1497 TLatex *latexInfo1 = new TLatex();
1498 latexInfo1->SetNDC();
1499 latexInfo1->SetTextSize(0.045);
1500 latexInfo1->SetTextColor(1);
1501 latexInfo1->DrawLatex(0.54,0.88,multLatex.Data());
1502 latexInfo1->DrawLatex(0.54,0.82,pttLatex.Data());
1503 latexInfo1->DrawLatex(0.54,0.76,ptaLatex.Data());
1504
1505 TString pngName = "BalanceFunction2D.";
86044267 1506 pngName += Form("%s: %.1f - %.1f",eventClass.Data(),psiMin,psiMax);
79dc837f 1507 pngName += ".PttFrom";
1508 pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1509 pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1510 pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1511 pngName += Form("%.1f",ptAssociatedMax);
1512 pngName += ".png";
1513
1514 c0->SaveAs(pngName.Data());
1515
a91351d9
MW		 1516 // do the full range for the projection on eta (for cross checking with published results)
1517 if(bFullPhiForEtaProjection){
1518
1519 drawProjections(gHistBalanceFunction2D,
1520 kTRUE,
1521 1,72,
1522 gCentrality,
1523 psiMin,psiMax,
1524 ptTriggerMin,ptTriggerMax,
1525 ptAssociatedMin,ptAssociatedMax,
1526 kTRUE,
1527 eventClass.Data(),
40a8b898
MW		 1528 bRootMoments,
1529 kFALSE,
7827f01e 1530 bReduceRangeForMoments,
1531 bFWHM);
a91351d9
MW		 1532 }
1533 else{
7827f01e 1534 drawProjections(gHistBalanceFunction2D,
a91351d9
MW		 1535 kTRUE,
1536 1,36,
1537 gCentrality,
1538 psiMin,psiMax,
1539 ptTriggerMin,ptTriggerMax,
1540 ptAssociatedMin,ptAssociatedMax,
1541 kTRUE,
1542 eventClass.Data(),
7827f01e 1543 bRootMoments,
1544 kFALSE,
1545 bReduceRangeForMoments,
1546 bFWHM);
a91351d9 1547 }
79dc837f 1548
1549 drawProjections(gHistBalanceFunction2D,
86044267 1550 kFALSE,
1551 1,80,
1552 gCentrality,
1553 psiMin,psiMax,
1554 ptTriggerMin,ptTriggerMax,
1555 ptAssociatedMin,ptAssociatedMax,
1556 kTRUE,
1557 eventClass.Data(),
40a8b898
MW		 1558 bRootMoments,
1559 kFALSE,
7827f01e 1560 bReduceRangeForMoments,
1561 bFWHM);
3de51a46 1562
1563 TString outFileName = filename;
1564 outFileName.ReplaceAll("correlationFunction","balanceFunction2D");
1565 gHistBalanceFunction2D->SetName("gHistBalanceFunctionSubtracted");
1566 TFile *fOut = TFile::Open(outFileName.Data(),"recreate");
1567 gHistBalanceFunction2D->Write();
1568 fOut->Close();
1569
1570}
1571
6f14379e 1572
1573//____________________________________________________________//
1574void mergeBFPsi2D(TString momDirectory = "./",
1575 TString directory1 = "",
1576 TString directory2 = "",
1577 TString directory3 = "",
1578 TString directory4 = "",
1579 Int_t gCentrality = 1,
1580 Double_t psiMin = -0.5, Double_t psiMax = 3.5,
1581 TString eventClass = "Centrality",
1582 Double_t ptTriggerMin = -1.,
1583 Double_t ptTriggerMax = -1.,
1584 Double_t ptAssociatedMin = -1.,
40a8b898 1585 Double_t ptAssociatedMax = -1.,
864839b8 1586 Bool_t bReduceRangeForMoments = kFALSE,
1587 Bool_t bFWHM = kFALSE
6f14379e 1588) {
1589 //Macro that draws the BF distributions for each centrality bin
1590 //for reaction plane dependent analysis
1591 //Author: Panos.Christakoglou@nikhef.nl
1592 TGaxis::SetMaxDigits(3);
1593 gStyle->SetPalette(55,0);
1594
1595 const Int_t nMaxDirectories = 4; // maximum number of directories to merge (set to 4 for now)
1596 TString sDirectory[nMaxDirectories];
1597 Int_t nDirectories = nMaxDirectories;
1598
1599 TString sFileName[nMaxDirectories];
1600 TFile *fBF[nMaxDirectories];
1601 TH2D *hBF[nMaxDirectories];
1602 Double_t entries[nMaxDirectories];
1603
1604 TFile *fOut;
1605 TH2D *hBFOut;
1606 Double_t entriesOut = 0.;
1607
1608 // find out how many directories we have to merge
1609 if(directory1==""){
1610 nDirectories=0;
1611 }
1612 else if(directory2==""){
1613 nDirectories=1;
1614 sDirectory[0]=directory1;
1615 }
1616 else if(directory3==""){
1617 nDirectories=2;
1618 sDirectory[0]=directory1;
1619 sDirectory[1]=directory2;
1620 }
1621 else if(directory4==""){
1622 nDirectories=3;
1623 sDirectory[0]=directory1;
1624 sDirectory[1]=directory2;
1625 sDirectory[2]=directory3;
1626 }
1627 else{
1628 nDirectories=nMaxDirectories;
1629 sDirectory[0]=directory1;
1630 sDirectory[1]=directory2;
1631 sDirectory[2]=directory3;
1632 sDirectory[3]=directory4;
1633 }
1634
1635 for(Int_t iDir = 0; iDir < nDirectories; iDir++){
1636
1637 //Get the input file
1638 sFileName[iDir] = sDirectory[iDir];
1639 sFileName[iDir] += "/Centrality"; sFileName[iDir] += gCentrality;
1640 sFileName[iDir] += "/";
1641sFileName[iDir] += momDirectory;
1642 sFileName[iDir] += "/balanceFunction2D.";
1643
1644 if(eventClass == "Centrality"){
1645 sFileName[iDir] += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1646 sFileName[iDir] += ".PsiAll.PttFrom";
1647 }
1648 else if(eventClass == "Multiplicity"){
1649 sFileName[iDir] += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1650 sFileName[iDir] += ".PsiAll.PttFrom";
1651 }
1652 else{ // "EventPlane" (default)
1653 sFileName[iDir] += "Centrality";
1654 sFileName[iDir] += gCentrality; sFileName[iDir] += ".Psi";
1655 if((psiMin == -0.5)&&(psiMax == 0.5)) sFileName[iDir] += "InPlane.Ptt";
1656 else if((psiMin == 0.5)&&(psiMax == 1.5)) sFileName[iDir] += "Intermediate.Ptt";
1657 else if((psiMin == 1.5)&&(psiMax == 2.5)) sFileName[iDir] += "OutOfPlane.Ptt";
1658 else if((psiMin == 2.5)&&(psiMax == 3.5)) sFileName[iDir] += "Rest.PttFrom";
1659 else sFileName[iDir] += "All.PttFrom";
1660 }
1661 sFileName[iDir] += Form("%.1f",ptTriggerMin); sFileName[iDir] += "To";
1662 sFileName[iDir] += Form("%.1f",ptTriggerMax); sFileName[iDir] += "PtaFrom";
1663 sFileName[iDir] += Form("%.1f",ptAssociatedMin); sFileName[iDir] += "To";
1664 sFileName[iDir] += Form("%.1f",ptAssociatedMax);
1665 sFileName[iDir] += "_";
1666 sFileName[iDir] += Form("%.1f",psiMin);
1667 sFileName[iDir] += "-";
1668 sFileName[iDir] += Form("%.1f",psiMax);
1669 sFileName[iDir] += ".root";
1670
1671 //Open the file
1672 fBF[iDir] = TFile::Open(sFileName[iDir].Data());
1673 if((!fBF[iDir])||(!fBF[iDir]->IsOpen())) {
1674 Printf("The file %s is not found. Not used...",sFileName[iDir]);
1675 continue;
1676 }
1677 //fBF[iDir]->ls();
1678
1679 hBF[iDir] = (TH2D*)fBF[iDir]->Get("gHistBalanceFunctionSubtracted");
1680 if(!hBF[iDir]) continue;
1681 entries[iDir] = (Double_t) hBF[iDir]->GetEntries();
1682 entriesOut += entries[iDir];
1683 cout<<" BF histogram "<<iDir<<" has "<<entries[iDir] <<" entries."<<endl;
1684
1685 // scaling and adding (for average)
1686 hBF[iDir]->Scale(entries[iDir]);
1687 if(iDir==0) hBFOut = (TH2D*)hBF[iDir]->Clone("gHistBalanceFunctionSubtractedOut");
1688 else hBFOut->Add(hBF[iDir]);
1689
1690 }
1691
1692 // scaling with all entries
1693 hBFOut->Scale(1./entriesOut);
40a8b898
MW		 1694
1695 drawProjections(hBFOut,
7827f01e 1696 kTRUE,
1697 1,36,
1698 gCentrality,
1699 psiMin,psiMax,
1700 ptTriggerMin,ptTriggerMax,
1701 ptAssociatedMin,ptAssociatedMax,
1702 kTRUE,
1703 eventClass,
1704 kTRUE,
1705 kFALSE,
1706 bReduceRangeForMoments,
1707 bFWHM);
6f14379e 1708
1709 drawProjections(hBFOut,
1710 kFALSE,
1711 1,80,
1712 gCentrality,
1713 psiMin,psiMax,
1714 ptTriggerMin,ptTriggerMax,
1715 ptAssociatedMin,ptAssociatedMax,
1716 kTRUE,
1717 eventClass.Data(),
40a8b898
MW		 1718 kTRUE,
1719 kFALSE,
7827f01e 1720 bReduceRangeForMoments,
1721 bFWHM);
6f14379e 1722
1723 // output
1724 TString outFileName = "balanceFunction2D.";
1725
1726 if(eventClass == "Centrality"){
1727 outFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1728 outFileName += ".PsiAll.PttFrom";
1729 }
1730 else if(eventClass == "Multiplicity"){
1731 outFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1732 outFileName += ".PsiAll.PttFrom";
1733 }
1734 else{ // "EventPlane" (default)
1735 outFileName += "Centrality";
1736 outFileName += gCentrality; outFileName += ".Psi";
1737 if((psiMin == -0.5)&&(psiMax == 0.5)) outFileName += "InPlane.Ptt";
1738 else if((psiMin == 0.5)&&(psiMax == 1.5)) outFileName += "Intermediate.Ptt";
1739 else if((psiMin == 1.5)&&(psiMax == 2.5)) outFileName += "OutOfPlane.Ptt";
1740 else if((psiMin == 2.5)&&(psiMax == 3.5)) outFileName += "Rest.PttFrom";
1741 else outFileName += "All.PttFrom";
1742 }
1743 outFileName += Form("%.1f",ptTriggerMin); outFileName += "To";
1744 outFileName += Form("%.1f",ptTriggerMax); outFileName += "PtaFrom";
1745 outFileName += Form("%.1f",ptAssociatedMin); outFileName += "To";
1746 outFileName += Form("%.1f",ptAssociatedMax);
1747 outFileName += "_";
1748 outFileName += Form("%.1f",psiMin);
1749 outFileName += "-";
1750 outFileName += Form("%.1f",psiMax);
1751 outFileName += ".root";
1752
1753 hBFOut->SetName("gHistBalanceFunctionSubtracted");
1754 fBFOut = TFile::Open(outFileName.Data(),"recreate");
1755 hBFOut->Write();
1756 fBFOut->Close();
1757
1758}
1759
3de51a46 1760//____________________________________________________________________//
40a8b898 1761void GetWeightedMean1D(TH1D *gHistBalance, Bool_t kProjectInEta = kTRUE, Int_t fStartBin = 1, Int_t fStopBin = -1, Double_t &WM, Double_t &WME,Double_t rangeReduced = -1.) {
3de51a46 1762
1763 //Prints the calculated width of the BF and its error
1764 Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
1765 Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
1766 Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
1767 Double_t deltaBalP2 = 0.0, integral = 0.0;
1768 Double_t deltaErrorNew = 0.0;
1769
1770 //Retrieve this variables from Histogram
1771 Int_t fNumberOfBins = gHistBalance->GetNbinsX();
1772 if(fStopBin > -1) fNumberOfBins = fStopBin;
1773 Double_t fP2Step = gHistBalance->GetBinWidth(1); // assume equal binning!
1774 Double_t currentBinCenter = 0.;
1775
1776 for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
1777
1778 // in order to recover the |abs| in the 1D analysis
1779 currentBinCenter = gHistBalance->GetBinCenter(i);
1780 if(kProjectInEta){
1781 if(currentBinCenter < 0) currentBinCenter = -currentBinCenter;
1782 }
1783 else{
1784 if(currentBinCenter < 0) currentBinCenter = -currentBinCenter;
1785 if(currentBinCenter > TMath::Pi()) currentBinCenter = 2 * TMath::Pi() - currentBinCenter;
1786 }
1787
40a8b898
MW		 1788 if(rangeReduced > 0 && currentBinCenter > rangeReduced )
1789 continue;
1790
3de51a46 1791 gSumXi += currentBinCenter;
1792 gSumBi += gHistBalance->GetBinContent(i);
1793 gSumBiXi += gHistBalance->GetBinContent(i)*currentBinCenter;
1794 gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(currentBinCenter,2);
1795 gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(currentBinCenter,2);
1796 gSumDeltaBi2 += TMath::Power(gHistBalance->GetBinError(i),2);
1797 gSumXi2DeltaBi2 += TMath::Power(currentBinCenter,2) * TMath::Power(gHistBalance->GetBinError(i),2);
1798
1799 deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
1800 integral += fP2Step*gHistBalance->GetBinContent(i);
1801 }
1802 for(Int_t i = fStartBin; i < fNumberOfBins; i++)
1803 deltaErrorNew += gHistBalance->GetBinError(i)*(currentBinCenter*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
1804
1805 Double_t integralError = TMath::Sqrt(deltaBalP2);
1806
1807 Double_t delta = gSumBiXi / gSumBi;
1808 Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
1809
1810 WM = delta;
1811 WME = deltaError;
79dc837f 1812}
3de51a46 1813