]> 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
-add random variable
[u/mrichter/AliRoot.git] / PWGCF / EBYE / macros / drawBalanceFunction2DPsi.C
CommitLineData
e4391c02 1const Int_t numberOfCentralityBins = 12;
a9f20288 2TString centralityArray[numberOfCentralityBins] = {"0-100","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-100","0-1","1-2","2-3"};
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);
27634c13 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,
935 Bool_t kUseZYAM = kFALSE) {
438dab2c 936 //Macro that draws the balance functions PROJECTIONS
17c69b4e 937 //for each centrality bin for the different pT of trigger and
938 //associated particles
939 TGaxis::SetMaxDigits(3);
940
941 //first we need some libraries
a9f20288 942 gSystem->Load("libTree");
943 gSystem->Load("libGeom");
944 gSystem->Load("libVMC");
945 gSystem->Load("libXMLIO");
946 gSystem->Load("libPhysics");
947
948 gSystem->Load("libSTEERBase");
949 gSystem->Load("libESD");
950 gSystem->Load("libAOD");
951
17c69b4e 952 gSystem->Load("libANALYSIS.so");
953 gSystem->Load("libANALYSISalice.so");
954 gSystem->Load("libEventMixing.so");
955 gSystem->Load("libCORRFW.so");
956 gSystem->Load("libPWGTools.so");
957 gSystem->Load("libPWGCFebye.so");
958
438dab2c 959 gStyle->SetOptStat(0);
960
17c69b4e 961 //Get the input file
962 TString filename = "balanceFunction2D.";
963 if(eventClass == "Centrality"){
964 filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
965 filename += ".PsiAll.PttFrom";
966 }
967 else if(eventClass == "Multiplicity"){
968 filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
969 filename += ".PsiAll.PttFrom";
970 }
971 else{ // "EventPlane" (default)
972 filename += "Centrality";
973 filename += gCentrality; filename += ".Psi";
974 if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
975 else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
976 else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
977 else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
978 else filename += "All.PttFrom";
979 }
980 filename += Form("%.1f",ptTriggerMin); filename += "To";
981 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
982 filename += Form("%.1f",ptAssociatedMin); filename += "To";
983 filename += Form("%.1f",ptAssociatedMax);
984 if(k2pMethod) filename += "_2pMethod";
a9f20288 985
a2da950e 986 filename += "_";
987 filename += Form("%.1f",psiMin);
988 filename += "-";
989 filename += Form("%.1f",psiMax);
17c69b4e 990 filename += ".root";
991
992 //Open the file
79dc837f 993 TFile *f = 0x0;
994 if(!gHistBalanceFunction2D) {
995 TFile::Open(filename.Data());
996 if((!f)||(!f->IsOpen())) {
997 Printf("The file %s is not found. Aborting...",filename);
998 return listBF;
999 }
1000 //f->ls();
17c69b4e 1001 }
79dc837f 1002
17c69b4e 1003 //Latex
1004 TString centralityLatex = "Centrality: ";
1005 centralityLatex += centralityArray[gCentrality-1];
1006 centralityLatex += "%";
1007
1008 TString psiLatex;
1009 if((psiMin == -0.5)&&(psiMax == 0.5))
1010 psiLatex = " -7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o}";
1011 else if((psiMin == 0.5)&&(psiMax == 1.5))
1012 psiLatex = " 37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o}";
1013 else if((psiMin == 1.5)&&(psiMax == 2.5))
1014 psiLatex = " 82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o}";
1015 else
1016 psiLatex = " 0^{o} < #varphi^{t} - #Psi_{2} < 180^{o}";
1017
1018 TString pttLatex = Form("%.1f",ptTriggerMin);
1019 pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
1020 pttLatex += " GeV/c";
1021
1022 TString ptaLatex = Form("%.1f",ptAssociatedMin);
1023 ptaLatex += " < p_{T,assoc} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1024 ptaLatex += " GeV/c";
1025
1026 TLatex *latexInfo1 = new TLatex();
1027 latexInfo1->SetNDC();
1028 latexInfo1->SetTextSize(0.045);
1029 latexInfo1->SetTextColor(1);
1030
17c69b4e 1031
1032 //============================================================//
1033 //Get subtracted and mixed balance function
79dc837f 1034 TH2D *gHistBalanceFunctionSubtracted2D = 0x0;
1035 TH2D *gHistBalanceFunctionMixed2D = 0x0;
1036 if(!gHistBalanceFunction2D) {
1037 gHistBalanceFunctionSubtracted2D = (TH2D*)f->Get("gHistBalanceFunctionSubtracted");
1038 gHistBalanceFunctionMixed2D = (TH2D*)f->Get("gHistBalanceFunctionMixed");
1039 }
1040 else {
1041 gHistBalanceFunctionSubtracted2D = dynamic_cast<TH2D*>(gHistBalanceFunction2D->Clone());
1042 gHistBalanceFunctionMixed2D = dynamic_cast<TH2D*>(gHistBalanceFunction2D->Clone());
1043 }
17c69b4e 1044
1045 TH1D *gHistBalanceFunctionSubtracted = NULL;
1046 TH1D *gHistBalanceFunctionMixed = NULL;
1047
a9f20288 1048 TH1D *gHistBalanceFunctionSubtracted_scale = NULL;
1049 TH1D *gHistBalanceFunctionMixed_scale = NULL;
1050
17c69b4e 1051 if(kProjectInEta){
86044267 1052 gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionX("gHistBalanceFunctionSubtractedEta",binMin,binMax));
a9f20288 1053 gHistBalanceFunctionSubtracted->Scale(gHistBalanceFunctionSubtracted2D->GetYaxis()->GetBinWidth(1)); // to remove normalization to phi bin width
86044267 1054 gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionX("gHistBalanceFunctionMixedEta",binMin,binMax));
a9f20288 1055 gHistBalanceFunctionMixed->Scale(gHistBalanceFunctionMixed2D->GetYaxis()->GetBinWidth(1)); // to remove normalization to phi bin width
17c69b4e 1056 gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#eta)");
1057 gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#eta)");
1058 }
1059 else{
86044267 1060 gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionY("gHistBalanceFunctionSubtractedPhi",binMin,binMax));
a9f20288 1061 gHistBalanceFunctionSubtracted->Scale(gHistBalanceFunctionSubtracted2D->GetXaxis()->GetBinWidth(1)); // to remove normalization to eta bin width
86044267 1062 gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionY("gHistBalanceFunctionMixedPhi",binMin,binMax));
a9f20288 1063 gHistBalanceFunctionMixed->Scale(gHistBalanceFunctionMixed2D->GetXaxis()->GetBinWidth(1)); // to remove normalization to eta bin width
17c69b4e 1064 gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#varphi)");
1065 gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#varphi)");
1066 }
1067
1068 gHistBalanceFunctionSubtracted->SetMarkerStyle(20);
1069 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
17c69b4e 1070
1071 gHistBalanceFunctionMixed->SetMarkerStyle(25);
17c69b4e 1072
1073 TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
1074 c1->SetFillColor(10);
1075 c1->SetHighLightColor(10);
1076 c1->SetLeftMargin(0.15);
1077 gHistBalanceFunctionSubtracted->DrawCopy("E");
a9f20288 1078 gHistBalanceFunctionMixed->DrawCopy("ESAME");
17c69b4e 1079
1080 legend = new TLegend(0.18,0.62,0.45,0.82,"","brNDC");
1081 legend->SetTextSize(0.045);
1082 legend->SetTextFont(42);
1083 legend->SetBorderSize(0);
1084 legend->SetFillStyle(0);
1085 legend->SetFillColor(10);
1086 legend->SetMargin(0.25);
1087 legend->SetShadowColor(10);
1088 legend->AddEntry(gHistBalanceFunctionSubtracted,"Data","lp");
1089 legend->AddEntry(gHistBalanceFunctionMixed,"Mixed data","lp");
1090 legend->Draw();
1091
17c69b4e 1092
17c69b4e 1093 TString meanLatex, rmsLatex, skewnessLatex, kurtosisLatex;
1094
1095 if(bRootMoments){
732acb5b 1096
1097 // need to restrict to near side in case of dphi
1098 if(!kProjectInEta) gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-TMath::Pi()/2,TMath::Pi()/2);
1099
17c69b4e 1100 meanLatex = "#mu = ";
1101 meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMean());
1102 meanLatex += " #pm ";
1103 meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMeanError());
1104
1105 rmsLatex = "#sigma = ";
1106 rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMS());
1107 rmsLatex += " #pm ";
1108 rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMSError());
1109
1110 skewnessLatex = "S = ";
1111 skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(1));
1112 skewnessLatex += " #pm ";
1113 skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(11));
1114
1115 kurtosisLatex = "K = ";
1116 kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(1));
1117 kurtosisLatex += " #pm ";
1118 kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(11));
a2da950e 1119
732acb5b 1120 Printf("ROOT Mean: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetMean(),gHistBalanceFunctionSubtracted->GetMeanError());
1121 Printf("ROOT RMS: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetRMS(),gHistBalanceFunctionSubtracted->GetRMSError());
1122 Printf("ROOT Skeweness: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetSkewness(1),gHistBalanceFunctionSubtracted->GetSkewness(11));
1123 Printf("ROOT Kurtosis: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetKurtosis(1),gHistBalanceFunctionSubtracted->GetKurtosis(11));
438dab2c 1124
732acb5b 1125
438dab2c 1126 // store in txt files
438dab2c 1127 TString meanFileName = filename;
79dc837f 1128 if(kProjectInEta)
23a0fc15 1129 meanFileName= "deltaEtaProjection_Mean.txt";
1130 //meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
79dc837f 1131 else
23a0fc15 1132 meanFileName = "deltaPhiProjection_Mean.txt";
1133 //meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
79dc837f 1134 ofstream fileMean(meanFileName.Data(),ios::app);
438dab2c 1135 fileMean << " " << gHistBalanceFunctionSubtracted->GetMean() << " " <<gHistBalanceFunctionSubtracted->GetMeanError()<<endl;
1136 fileMean.close();
1137
1138 TString rmsFileName = filename;
79dc837f 1139 if(kProjectInEta)
23a0fc15 1140 rmsFileName = "deltaEtaProjection_Rms.txt";
1141 //rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
79dc837f 1142 else
732acb5b 1143 rmsFileName = "deltaPhiProjection_Rms.txt";
23a0fc15 1144 //rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
79dc837f 1145 ofstream fileRms(rmsFileName.Data(),ios::app);
438dab2c 1146 fileRms << " " << gHistBalanceFunctionSubtracted->GetRMS() << " " <<gHistBalanceFunctionSubtracted->GetRMSError()<<endl;
1147 fileRms.close();
1148
1149 TString skewnessFileName = filename;
79dc837f 1150 if(kProjectInEta)
23a0fc15 1151 skewnessFileName = "deltaEtaProjection_Skewness.txt";
1152 //skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
79dc837f 1153 else
23a0fc15 1154 skewnessFileName = "deltaPhiProjection_Skewness.txt";
1155 //skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
79dc837f 1156 ofstream fileSkewness(skewnessFileName.Data(),ios::app);
438dab2c 1157 fileSkewness << " " << gHistBalanceFunctionSubtracted->GetSkewness(1) << " " <<gHistBalanceFunctionSubtracted->GetSkewness(11)<<endl;
1158 fileSkewness.close();
1159
1160 TString kurtosisFileName = filename;
79dc837f 1161 if(kProjectInEta)
23a0fc15 1162 kurtosisFileName = "deltaEtaProjection_Kurtosis.txt";
1163 //kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
79dc837f 1164 else
23a0fc15 1165 kurtosisFileName = "deltaPhiProjection_Kurtosis.txt";
1166 //kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
79dc837f 1167 ofstream fileKurtosis(kurtosisFileName.Data(),ios::app);
438dab2c 1168 fileKurtosis << " " << gHistBalanceFunctionSubtracted->GetKurtosis(1) << " " <<gHistBalanceFunctionSubtracted->GetKurtosis(11)<<endl;
1169 fileKurtosis.close();
732acb5b 1170
1171 if(!kProjectInEta) gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-TMath::Pi()/2,3.*TMath::Pi()/2);
1172
17c69b4e 1173 }
1174 // calculate the moments by hand
1175 else{
1176
1177 Double_t meanAnalytical, meanAnalyticalError;
1178 Double_t sigmaAnalytical, sigmaAnalyticalError;
1179 Double_t skewnessAnalytical, skewnessAnalyticalError;
1180 Double_t kurtosisAnalytical, kurtosisAnalyticalError;
1822002d 1181
17c69b4e 1182 Int_t gDeltaEtaPhi = 2;
1183 if(kProjectInEta) gDeltaEtaPhi = 1;
1184
1185 AliBalancePsi *bHelper = new AliBalancePsi;
1822002d 1186 bHelper->GetMomentsAnalytical(gDeltaEtaPhi,gHistBalanceFunctionSubtracted,kUseZYAM,meanAnalytical,meanAnalyticalError,sigmaAnalytical,sigmaAnalyticalError,skewnessAnalytical,skewnessAnalyticalError,kurtosisAnalytical,kurtosisAnalyticalError);
17c69b4e 1187
1188 meanLatex = "#mu = ";
1189 meanLatex += Form("%.3f",meanAnalytical);
1190 meanLatex += " #pm ";
1191 meanLatex += Form("%.3f",meanAnalyticalError);
1192
1193 rmsLatex = "#sigma = ";
1194 rmsLatex += Form("%.3f",sigmaAnalytical);
1195 rmsLatex += " #pm ";
1196 rmsLatex += Form("%.3f",sigmaAnalyticalError);
1197
1198 skewnessLatex = "S = ";
1199 skewnessLatex += Form("%.3f",skewnessAnalytical);
1200 skewnessLatex += " #pm ";
1201 skewnessLatex += Form("%.3f",skewnessAnalyticalError);
1202
1203 kurtosisLatex = "K = ";
1204 kurtosisLatex += Form("%.3f",kurtosisAnalytical);
1205 kurtosisLatex += " #pm ";
1206 kurtosisLatex += Form("%.3f",kurtosisAnalyticalError);
1207 Printf("Mean: %lf - Error: %lf",meanAnalytical, meanAnalyticalError);
1208 Printf("Sigma: %lf - Error: %lf",sigmaAnalytical, sigmaAnalyticalError);
1209 Printf("Skeweness: %lf - Error: %lf",skewnessAnalytical, skewnessAnalyticalError);
1210 Printf("Kurtosis: %lf - Error: %lf",kurtosisAnalytical, kurtosisAnalyticalError);
438dab2c 1211
1212 // store in txt files
1213 TString meanFileName = filename;
79dc837f 1214 if(kProjectInEta)
23a0fc15 1215 meanFileName = "deltaEtaProjection_Mean.txt";
1216 //meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
79dc837f 1217 else
23a0fc15 1218 meanFileName = "deltaPhiProjection_Mean.txt";
1219 //meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
79dc837f 1220 ofstream fileMean(meanFileName.Data(),ios::app);
438dab2c 1221 fileMean << " " << meanAnalytical << " " <<meanAnalyticalError<<endl;
1222 fileMean.close();
1223
1224 TString rmsFileName = filename;
79dc837f 1225 if(kProjectInEta)
23a0fc15 1226 rmsFileName = "deltaEtaProjection_Rms.txt";
1227 //rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
79dc837f 1228 else
23a0fc15 1229 rmsFileName = "deltaPhiProjection_Rms.txt";
1230//rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
79dc837f 1231 ofstream fileRms(rmsFileName.Data(),ios::app);
438dab2c 1232 fileRms << " " << sigmaAnalytical << " " <<sigmaAnalyticalError<<endl;
1233 fileRms.close();
1234
1235 TString skewnessFileName = filename;
79dc837f 1236 if(kProjectInEta)
23a0fc15 1237 skewnessFileName = "deltaEtaProjection_Skewness.txt";
1238 //skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
79dc837f 1239 else
23a0fc15 1240 skewnessFileName = "deltaPhiProjection_Skewness.txt";
1241 //skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
79dc837f 1242 ofstream fileSkewness(skewnessFileName.Data(),ios::app);
438dab2c 1243 fileSkewness << " " << skewnessAnalytical << " " <<skewnessAnalyticalError<<endl;
1244 fileSkewness.close();
1245
1246 TString kurtosisFileName = filename;
79dc837f 1247 if(kProjectInEta)
23a0fc15 1248 kurtosisFileName = "deltaEtaProjection_Kurtosis.txt";
1249 //kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
79dc837f 1250 else
23a0fc15 1251 kurtosisFileName = "deltaPhiProjection_Kurtosis.txt";
1252 //kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
79dc837f 1253 ofstream fileKurtosis(kurtosisFileName.Data(),ios::app);
438dab2c 1254 fileKurtosis << " " << kurtosisAnalytical << " " <<kurtosisAnalyticalError<<endl;
1255 fileKurtosis.close();
17c69b4e 1256 }
1257
3de51a46 1258 // Weighted mean as calculated for 1D analysis
1259 Double_t weightedMean, weightedMeanError;
1260 GetWeightedMean1D(gHistBalanceFunctionSubtracted,kProjectInEta,1,-1,weightedMean,weightedMeanError);
1261 Printf("Weighted Mean: %lf - Error: %lf",weightedMean, weightedMeanError);
1262
1263 // store in txt files
1264 TString weightedMeanFileName = filename;
1265 if(kProjectInEta)
1266 weightedMeanFileName = "deltaEtaProjection_WeightedMean.txt";
1267 //weightedMeanFileName.ReplaceAll(".root","_DeltaEtaProjection_WeightedMean.txt");
1268 else
1269 weightedMeanFileName = "deltaPhiProjection_WeightedMean.txt";
1270 //weightedMeanFileName.ReplaceAll(".root","_DeltaPhiProjection_WeightedMean.txt");
1271 ofstream fileWeightedMean(weightedMeanFileName.Data(),ios::app);
1272 fileWeightedMean << " " << weightedMean << " " <<weightedMeanError<<endl;
1273 fileWeightedMean.close();
438dab2c 1274
1275
17c69b4e 1276 TCanvas *c2 = new TCanvas("c2","",600,0,600,500);
1277 c2->SetFillColor(10);
1278 c2->SetHighLightColor(10);
1279 c2->SetLeftMargin(0.15);
1280 gHistBalanceFunctionSubtracted->DrawCopy("E");
438dab2c 1281
17c69b4e 1282 TLatex *latex = new TLatex();
1283 latex->SetNDC();
1284 latex->SetTextSize(0.035);
1285 latex->SetTextColor(1);
1286 latex->DrawLatex(0.64,0.85,meanLatex.Data());
1287 latex->DrawLatex(0.64,0.81,rmsLatex.Data());
1288 latex->DrawLatex(0.64,0.77,skewnessLatex.Data());
1289 latex->DrawLatex(0.64,0.73,kurtosisLatex.Data());
1290
438dab2c 1291 TString pngName = filename;
1292 if(kProjectInEta) pngName.ReplaceAll(".root","_DeltaEtaProjection.png");
1293 else pngName.ReplaceAll(".root","_DeltaPhiProjection.png");
1294
1295 c2->SaveAs(pngName.Data());
1296
17c69b4e 1297 TString outFileName = filename;
1298 if(kProjectInEta) outFileName.ReplaceAll(".root","_DeltaEtaProjection.root");
1299 else outFileName.ReplaceAll(".root","_DeltaPhiProjection.root");
1300 TFile *fProjection = TFile::Open(outFileName.Data(),"recreate");
1301 gHistBalanceFunctionSubtracted->Write();
1302 gHistBalanceFunctionMixed->Write();
1303 fProjection->Close();
1304}
79dc837f 1305
1306//____________________________________________________________//
1307void drawBFPsi2DFromCorrelationFunctions(const char* lhcPeriod = "LHC10h",
1308 Int_t gTrainNumber = 64,
1309 const char* gCentralityEstimator = "V0M",
1310 Int_t gBit = 128,
1311 const char* gEventPlaneEstimator = "VZERO",
1312 Int_t gCentrality = 1,
1313 Double_t psiMin = -0.5, Double_t psiMax = 3.5,
1314 Double_t vertexZMin = -10.,
1315 Double_t vertexZMax = 10.,
1316 Double_t ptTriggerMin = -1.,
1317 Double_t ptTriggerMax = -1.,
1318 Double_t ptAssociatedMin = -1.,
86044267 1319 Double_t ptAssociatedMax = -1.,
732acb5b 1320 TString eventClass = "Multiplicity",
1321 Bool_t bRootMoments = kFALSE
86044267 1322) {
79dc837f 1323 //Macro that draws the BF distributions for each centrality bin
1324 //for reaction plane dependent analysis
1325 //Author: Panos.Christakoglou@nikhef.nl
1326 TGaxis::SetMaxDigits(3);
1327 gStyle->SetPalette(55,0);
1328
1329 //Get the input file
1330 TString filename = lhcPeriod;
86044267 1331 if(lhcPeriod != ""){
1332 //filename += "/Train"; filename += gTrainNumber;
1333 filename +="/PttFrom";
1334 filename += Form("%.1f",ptTriggerMin); filename += "To";
1335 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1336 filename += Form("%.1f",ptAssociatedMin); filename += "To";
1337 filename += Form("%.1f",ptAssociatedMax);
1338 filename += "/correlationFunction.";
1339 }
1340 else{
1341 filename += "correlationFunction.";
1342 }
1343 if(eventClass == "Centrality"){
1344 filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1345 filename += ".PsiAll.PttFrom";
1346 }
1347 else if(eventClass == "Multiplicity"){
1348 filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1349 filename += ".PsiAll.PttFrom";
1350 }
1351 else{ // "EventPlane" (default)
1352 filename += "Centrality";
1353 filename += gCentrality; filename += ".Psi";
1354 if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
1355 else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
1356 else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
1357 else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
1358 else filename += "All.PttFrom";
1359 }
79dc837f 1360 filename += Form("%.1f",ptTriggerMin); filename += "To";
1361 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1362 filename += Form("%.1f",ptAssociatedMin); filename += "To";
1363 filename += Form("%.1f",ptAssociatedMax);
1364 filename += "_";
1365 filename += Form("%.1f",psiMin);
1366 filename += "-";
1367 filename += Form("%.1f",psiMax);
1368 filename += ".root";
1369
1370 //Open the file
1371 TFile *f = TFile::Open(filename.Data());
1372 if((!f)||(!f->IsOpen())) {
1373 Printf("The file %s is not found. Aborting...",filename);
1374 return listBF;
1375 }
1376 //f->ls();
1377
1378 TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
1379 if(!gHistPN) return;
1380 TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
1381 if(!gHistNP) return;
1382 TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
1383 if(!gHistPP) return;
1384 TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
1385 if(!gHistNN) return;
1386
1387 gHistPN->Sumw2();
1388 gHistPP->Sumw2();
1389 gHistPN->Add(gHistPP,-1);
1390 gHistNP->Sumw2();
1391 gHistNN->Sumw2();
1392 gHistNP->Add(gHistNN,-1);
1393 gHistPN->Add(gHistNP);
1394 gHistPN->Scale(0.5);
1395 TH2D *gHistBalanceFunction2D = dynamic_cast<TH2D *>(gHistPN->Clone());
1396 gHistBalanceFunction2D->SetStats(kFALSE);
1397 gHistBalanceFunction2D->GetXaxis()->SetTitle("#Delta#eta");
1398 gHistBalanceFunction2D->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1399 gHistBalanceFunction2D->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
1400
1401 //Draw the results
1402 TCanvas *c0 = new TCanvas("c0","Balance function 2D",0,0,600,500);
1403 c0->SetFillColor(10); c0->SetHighLightColor(10);
1404 c0->SetLeftMargin(0.17); c0->SetTopMargin(0.05);
1405 gHistBalanceFunction2D->SetTitle("");
1406 gHistBalanceFunction2D->GetZaxis()->SetTitleOffset(1.4);
1407 gHistBalanceFunction2D->GetZaxis()->SetNdivisions(10);
1408 gHistBalanceFunction2D->GetYaxis()->SetTitleOffset(1.4);
1409 gHistBalanceFunction2D->GetYaxis()->SetNdivisions(10);
1410 gHistBalanceFunction2D->GetXaxis()->SetRangeUser(-1.4,1.4);
1411 gHistBalanceFunction2D->GetXaxis()->SetNdivisions(10);
1412 gHistBalanceFunction2D->DrawCopy("lego2");
1413 gPad->SetTheta(30); // default is 30
1414 gPad->SetPhi(-60); // default is 30
1415 gPad->Update();
1416
1417 TString multLatex = Form("Multiplicity: %.1f - %.1f",psiMin,psiMax);
1418
1419 TString pttLatex = Form("%.1f",ptTriggerMin);
1420 pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
1421 pttLatex += " GeV/c";
1422
1423 TString ptaLatex = Form("%.1f",ptAssociatedMin);
1424 ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1425 ptaLatex += " GeV/c";
1426
1427 TLatex *latexInfo1 = new TLatex();
1428 latexInfo1->SetNDC();
1429 latexInfo1->SetTextSize(0.045);
1430 latexInfo1->SetTextColor(1);
1431 latexInfo1->DrawLatex(0.54,0.88,multLatex.Data());
1432 latexInfo1->DrawLatex(0.54,0.82,pttLatex.Data());
1433 latexInfo1->DrawLatex(0.54,0.76,ptaLatex.Data());
1434
1435 TString pngName = "BalanceFunction2D.";
86044267 1436 pngName += Form("%s: %.1f - %.1f",eventClass.Data(),psiMin,psiMax);
79dc837f 1437 pngName += ".PttFrom";
1438 pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1439 pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1440 pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1441 pngName += Form("%.1f",ptAssociatedMax);
1442 pngName += ".png";
1443
1444 c0->SaveAs(pngName.Data());
1445
1446 drawProjections(gHistBalanceFunction2D,
1447 kTRUE,
1448 1,36,
1449 gCentrality,
1450 psiMin,psiMax,
1451 ptTriggerMin,ptTriggerMax,
1452 ptAssociatedMin,ptAssociatedMax,
1453 kTRUE,
86044267 1454 eventClass,
732acb5b 1455 bRootMoments);
79dc837f 1456
1457 drawProjections(gHistBalanceFunction2D,
86044267 1458 kFALSE,
1459 1,80,
1460 gCentrality,
1461 psiMin,psiMax,
1462 ptTriggerMin,ptTriggerMax,
1463 ptAssociatedMin,ptAssociatedMax,
1464 kTRUE,
1465 eventClass.Data(),
732acb5b 1466 bRootMoments);
3de51a46 1467
1468 TString outFileName = filename;
1469 outFileName.ReplaceAll("correlationFunction","balanceFunction2D");
1470 gHistBalanceFunction2D->SetName("gHistBalanceFunctionSubtracted");
1471 TFile *fOut = TFile::Open(outFileName.Data(),"recreate");
1472 gHistBalanceFunction2D->Write();
1473 fOut->Close();
1474
1475}
1476
1477//____________________________________________________________________//
1478void GetWeightedMean1D(TH1D *gHistBalance, Bool_t kProjectInEta = kTRUE, Int_t fStartBin = 1, Int_t fStopBin = -1, Double_t &WM, Double_t &WME) {
1479
1480 //Prints the calculated width of the BF and its error
1481 Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
1482 Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
1483 Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
1484 Double_t deltaBalP2 = 0.0, integral = 0.0;
1485 Double_t deltaErrorNew = 0.0;
1486
1487 //Retrieve this variables from Histogram
1488 Int_t fNumberOfBins = gHistBalance->GetNbinsX();
1489 if(fStopBin > -1) fNumberOfBins = fStopBin;
1490 Double_t fP2Step = gHistBalance->GetBinWidth(1); // assume equal binning!
1491 Double_t currentBinCenter = 0.;
1492
1493 for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
1494
1495 // in order to recover the |abs| in the 1D analysis
1496 currentBinCenter = gHistBalance->GetBinCenter(i);
1497 if(kProjectInEta){
1498 if(currentBinCenter < 0) currentBinCenter = -currentBinCenter;
1499 }
1500 else{
1501 if(currentBinCenter < 0) currentBinCenter = -currentBinCenter;
1502 if(currentBinCenter > TMath::Pi()) currentBinCenter = 2 * TMath::Pi() - currentBinCenter;
1503 }
1504
1505 gSumXi += currentBinCenter;
1506 gSumBi += gHistBalance->GetBinContent(i);
1507 gSumBiXi += gHistBalance->GetBinContent(i)*currentBinCenter;
1508 gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(currentBinCenter,2);
1509 gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(currentBinCenter,2);
1510 gSumDeltaBi2 += TMath::Power(gHistBalance->GetBinError(i),2);
1511 gSumXi2DeltaBi2 += TMath::Power(currentBinCenter,2) * TMath::Power(gHistBalance->GetBinError(i),2);
1512
1513 deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
1514 integral += fP2Step*gHistBalance->GetBinContent(i);
1515 }
1516 for(Int_t i = fStartBin; i < fNumberOfBins; i++)
1517 deltaErrorNew += gHistBalance->GetBinError(i)*(currentBinCenter*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
1518
1519 Double_t integralError = TMath::Sqrt(deltaBalP2);
1520
1521 Double_t delta = gSumBiXi / gSumBi;
1522 Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
1523
1524 WM = delta;
1525 WME = deltaError;
79dc837f 1526}
3de51a46 1527