]> 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
Extended pt range; corrected pi0/gamma weights
[u/mrichter/AliRoot.git] / PWGCF / EBYE / macros / drawBalanceFunction2DPsi.C
CommitLineData
e4391c02 1const Int_t numberOfCentralityBins = 12;
2TString centralityArray[numberOfCentralityBins] = {"0-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-100","0-1","1-2","2-3"};
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,
d67eae53 21 TString eventClass = "EventPlane") //Can be "EventPlane", "Centrality", "Multiplicity"
22{
6acdbcb2 23 //Macro that draws the BF distributions for each centrality bin
24 //for reaction plane dependent analysis
25 //Author: Panos.Christakoglou@nikhef.nl
26 //Load the PWG2ebye library
27 gSystem->Load("libANALYSIS.so");
28 gSystem->Load("libANALYSISalice.so");
29 gSystem->Load("libEventMixing.so");
30 gSystem->Load("libCORRFW.so");
31 gSystem->Load("libPWGTools.so");
32 gSystem->Load("libPWGCFebye.so");
33
9a0495b0 34 //gROOT->LoadMacro("~/SetPlotStyle.C");
35 //SetPlotStyle();
8795e993 36 gStyle->SetPalette(1,0);
15dd45a0 37
6acdbcb2 38 //Prepare the objects and return them
5de9ad1a 39 TList *listBF = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,0);
52daf7b2 40 TList *listBFShuffled = NULL;
5de9ad1a 41 if(kShowShuffled) listBFShuffled = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,1);
52daf7b2 42 TList *listBFMixed = NULL;
5de9ad1a 43 if(kShowMixed) listBFMixed = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,2);
6acdbcb2 44 if(!listBF) {
45 Printf("The TList object was not created");
46 return;
47 }
48 else
27634c13 49 draw(listBF,listBFShuffled,listBFMixed,gCentrality,gCentralityEstimator,
34616eda 50 psiMin,psiMax,vertexZMin,vertexZMax,
f0c5040c 51 ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,
832b21d8 52 kUseVzBinning,k2pMethod,eventClass);
6acdbcb2 53}
54
55//______________________________________________________//
5365d1d7 56TList *GetListOfObjects(const char* filename,
57 Int_t gCentrality,
58 Int_t gBit,
59 const char *gCentralityEstimator,
5de9ad1a 60 Int_t kData = 1) {
6acdbcb2 61 //Get the TList objects (QA, bf, bf shuffled)
6acdbcb2 62 TList *listBF = 0x0;
6acdbcb2 63
64 //Open the file
5365d1d7 65 TFile *f = TFile::Open(filename,"UPDATE");
6acdbcb2 66 if((!f)||(!f->IsOpen())) {
67 Printf("The file %s is not found. Aborting...",filename);
68 return listBF;
69 }
70 //f->ls();
71
72 TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
73 if(!dir) {
74 Printf("The TDirectoryFile is not found. Aborting...",filename);
75 return listBF;
76 }
77 //dir->ls();
78
79 TString listBFName;
eb63b883 80 if(kData == 0) {
81 //cout<<"no shuffling - no mixing"<<endl;
82 listBFName = "listBFPsi_";
83 }
84 else if(kData == 1) {
85 //cout<<"shuffling - no mixing"<<endl;
86 listBFName = "listBFPsiShuffled_";
87 }
88 else if(kData == 2) {
89 //cout<<"no shuffling - mixing"<<endl;
90 listBFName = "listBFPsiMixed_";
91 }
92 listBFName += centralityArray[gCentrality-1];
5de9ad1a 93 if(gBit > -1) {
94 listBFName += "_Bit"; listBFName += gBit; }
95 if(gCentralityEstimator) {
96 listBFName += "_"; listBFName += gCentralityEstimator;}
6acdbcb2 97
5365d1d7 98 // histograms were already retrieved (in first iteration)
99 if(dir->Get(Form("%s_histograms",listBFName.Data()))){
100 //cout<<"second iteration"<<endl;
101 listBF = dynamic_cast<TList *>(dir->Get(Form("%s_histograms",listBFName.Data())));
6acdbcb2 102 }
6acdbcb2 103
5365d1d7 104 // histograms were not yet retrieved (this is the first iteration)
105 else{
106
107 listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
108 cout<<"======================================================="<<endl;
5365d1d7 109 cout<<"List name: "<<listBF->GetName()<<endl;
110 //listBF->ls();
6acdbcb2 111
5365d1d7 112 //Get the histograms
113 TString histoName;
114 if(kData == 0)
27634c13 115 histoName = "fHistP";
5365d1d7 116 else if(kData == 1)
27634c13 117 histoName = "fHistP_shuffle";
5365d1d7 118 else if(kData == 2)
27634c13 119 histoName = "fHistP";
120 if(gCentralityEstimator) histoName += gCentralityEstimator;
5365d1d7 121 AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
122 if(!fHistP) {
123 Printf("fHistP %s not found!!!",histoName.Data());
124 break;
125 }
126 fHistP->FillParent(); fHistP->DeleteContainers();
127
128 if(kData == 0)
27634c13 129 histoName = "fHistN";
5365d1d7 130 if(kData == 1)
27634c13 131 histoName = "fHistN_shuffle";
5365d1d7 132 if(kData == 2)
27634c13 133 histoName = "fHistN";
134 if(gCentralityEstimator) histoName += gCentralityEstimator;
5365d1d7 135 AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
136 if(!fHistN) {
137 Printf("fHistN %s not found!!!",histoName.Data());
138 break;
139 }
140 fHistN->FillParent(); fHistN->DeleteContainers();
141
142 if(kData == 0)
27634c13 143 histoName = "fHistPN";
5365d1d7 144 if(kData == 1)
27634c13 145 histoName = "fHistPN_shuffle";
5365d1d7 146 if(kData == 2)
27634c13 147 histoName = "fHistPN";
148 if(gCentralityEstimator) histoName += gCentralityEstimator;
5365d1d7 149 AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
150 if(!fHistPN) {
151 Printf("fHistPN %s not found!!!",histoName.Data());
152 break;
153 }
154 fHistPN->FillParent(); fHistPN->DeleteContainers();
155
156 if(kData == 0)
27634c13 157 histoName = "fHistNP";
5365d1d7 158 if(kData == 1)
27634c13 159 histoName = "fHistNP_shuffle";
5365d1d7 160 if(kData == 2)
27634c13 161 histoName = "fHistNP";
162 if(gCentralityEstimator) histoName += gCentralityEstimator;
5365d1d7 163 AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
164 if(!fHistNP) {
165 Printf("fHistNP %s not found!!!",histoName.Data());
166 break;
167 }
168 fHistNP->FillParent(); fHistNP->DeleteContainers();
169
170 if(kData == 0)
27634c13 171 histoName = "fHistPP";
5365d1d7 172 if(kData == 1)
27634c13 173 histoName = "fHistPP_shuffle";
5365d1d7 174 if(kData == 2)
27634c13 175 histoName = "fHistPP";
176 if(gCentralityEstimator) histoName += gCentralityEstimator;
5365d1d7 177 AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
178 if(!fHistPP) {
179 Printf("fHistPP %s not found!!!",histoName.Data());
180 break;
181 }
182 fHistPP->FillParent(); fHistPP->DeleteContainers();
183
184 if(kData == 0)
27634c13 185 histoName = "fHistNN";
5365d1d7 186 if(kData == 1)
27634c13 187 histoName = "fHistNN_shuffle";
5365d1d7 188 if(kData == 2)
27634c13 189 histoName = "fHistNN";
190 if(gCentralityEstimator) histoName += gCentralityEstimator;
5365d1d7 191 AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
192 if(!fHistNN) {
193 Printf("fHistNN %s not found!!!",histoName.Data());
194 break;
195 }
196 fHistNN->FillParent(); fHistNN->DeleteContainers();
197
198 dir->cd();
199 listBF->Write(Form("%s_histograms",listBFName.Data()), TObject::kSingleKey);
200
201 }// first iteration
6acdbcb2 202
5365d1d7 203 f->Close();
6acdbcb2 204
205 return listBF;
206}
207
208//______________________________________________________//
eb63b883 209void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
27634c13 210 Int_t gCentrality, const char* gCentralityEstimator,
211 Double_t psiMin, Double_t psiMax,
34616eda 212 Double_t vertexZMin,
213 Double_t vertexZMax,
6acdbcb2 214 Double_t ptTriggerMin, Double_t ptTriggerMax,
f0c5040c 215 Double_t ptAssociatedMin, Double_t ptAssociatedMax,
832b21d8 216 Bool_t kUseVzBinning=kFALSE,
d67eae53 217 Bool_t k2pMethod = kFALSE, TString eventClass) {
6acdbcb2 218 //balance function
219 AliTHn *hP = NULL;
220 AliTHn *hN = NULL;
221 AliTHn *hPN = NULL;
222 AliTHn *hNP = NULL;
223 AliTHn *hPP = NULL;
224 AliTHn *hNN = NULL;
225 //listBF->ls();
226 //Printf("=================");
27634c13 227 TString histoName = "fHistP";
228 if(gCentralityEstimator) histoName += gCentralityEstimator;
229 hP = (AliTHn*) listBF->FindObject(histoName.Data());
52daf7b2 230 hP->SetName("gHistP");
27634c13 231 histoName = "fHistN";
232 if(gCentralityEstimator) histoName += gCentralityEstimator;
233 hN = (AliTHn*) listBF->FindObject(histoName.Data());
52daf7b2 234 hN->SetName("gHistN");
27634c13 235 histoName = "fHistPN";
236 if(gCentralityEstimator) histoName += gCentralityEstimator;
237 hPN = (AliTHn*) listBF->FindObject(histoName.Data());
52daf7b2 238 hPN->SetName("gHistPN");
27634c13 239 histoName = "fHistNP";
240 if(gCentralityEstimator) histoName += gCentralityEstimator;
241 hNP = (AliTHn*) listBF->FindObject(histoName.Data());
52daf7b2 242 hNP->SetName("gHistNP");
27634c13 243 histoName = "fHistPP";
244 if(gCentralityEstimator) histoName += gCentralityEstimator;
245 hPP = (AliTHn*) listBF->FindObject(histoName.Data());
52daf7b2 246 hPP->SetName("gHistPP");
27634c13 247 histoName = "fHistNN";
248 if(gCentralityEstimator) histoName += gCentralityEstimator;
249 hNN = (AliTHn*) listBF->FindObject(histoName.Data());
52daf7b2 250 hNN->SetName("gHistNN");
6acdbcb2 251
252 AliBalancePsi *b = new AliBalancePsi();
d67eae53 253 b->SetEventClass(eventClass);
6acdbcb2 254 b->SetHistNp(hP);
255 b->SetHistNn(hN);
256 b->SetHistNpn(hPN);
257 b->SetHistNnp(hNP);
258 b->SetHistNpp(hPP);
259 b->SetHistNnn(hNN);
832b21d8 260 if(kUseVzBinning) b->SetVertexZBinning(kTRUE);
261
6acdbcb2 262
263 //balance function shuffling
264 AliTHn *hPShuffled = NULL;
265 AliTHn *hNShuffled = NULL;
266 AliTHn *hPNShuffled = NULL;
267 AliTHn *hNPShuffled = NULL;
268 AliTHn *hPPShuffled = NULL;
269 AliTHn *hNNShuffled = NULL;
52daf7b2 270 if(listBFShuffled) {
271 //listBFShuffled->ls();
abf0c200 272 histoName = "fHistP_shuffle";
27634c13 273 if(gCentralityEstimator) histoName += gCentralityEstimator;
274 hPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
52daf7b2 275 hPShuffled->SetName("gHistPShuffled");
abf0c200 276 histoName = "fHistN_shuffle";
27634c13 277 if(gCentralityEstimator) histoName += gCentralityEstimator;
278 hNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
52daf7b2 279 hNShuffled->SetName("gHistNShuffled");
abf0c200 280 histoName = "fHistPN_shuffle";
27634c13 281 if(gCentralityEstimator) histoName += gCentralityEstimator;
282 hPNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
52daf7b2 283 hPNShuffled->SetName("gHistPNShuffled");
abf0c200 284 histoName = "fHistNP_shuffle";
27634c13 285 if(gCentralityEstimator) histoName += gCentralityEstimator;
286 hNPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
52daf7b2 287 hNPShuffled->SetName("gHistNPShuffled");
abf0c200 288 histoName = "fHistPP_shuffle";
27634c13 289 if(gCentralityEstimator) histoName += gCentralityEstimator;
290 hPPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
52daf7b2 291 hPPShuffled->SetName("gHistPPShuffled");
abf0c200 292 histoName = "fHistNN_shuffle";
27634c13 293 if(gCentralityEstimator) histoName += gCentralityEstimator;
294 hNNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
52daf7b2 295 hNNShuffled->SetName("gHistNNShuffled");
296
297 AliBalancePsi *bShuffled = new AliBalancePsi();
d67eae53 298 bShuffled->SetEventClass(eventClass);
52daf7b2 299 bShuffled->SetHistNp(hPShuffled);
300 bShuffled->SetHistNn(hNShuffled);
301 bShuffled->SetHistNpn(hPNShuffled);
302 bShuffled->SetHistNnp(hNPShuffled);
303 bShuffled->SetHistNpp(hPPShuffled);
304 bShuffled->SetHistNnn(hNNShuffled);
832b21d8 305 if(kUseVzBinning) bShuffled->SetVertexZBinning(kTRUE);
306
52daf7b2 307 }
6acdbcb2 308
eb63b883 309 //balance function mixing
310 AliTHn *hPMixed = NULL;
311 AliTHn *hNMixed = NULL;
312 AliTHn *hPNMixed = NULL;
313 AliTHn *hNPMixed = NULL;
314 AliTHn *hPPMixed = NULL;
315 AliTHn *hNNMixed = NULL;
eb63b883 316
52daf7b2 317 if(listBFMixed) {
318 //listBFMixed->ls();
27634c13 319 histoName = "fHistP";
320 if(gCentralityEstimator) histoName += gCentralityEstimator;
321 hPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
52daf7b2 322 hPMixed->SetName("gHistPMixed");
27634c13 323 histoName = "fHistN";
324 if(gCentralityEstimator) histoName += gCentralityEstimator;
325 hNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
52daf7b2 326 hNMixed->SetName("gHistNMixed");
27634c13 327 histoName = "fHistPN";
328 if(gCentralityEstimator) histoName += gCentralityEstimator;
329 hPNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
52daf7b2 330 hPNMixed->SetName("gHistPNMixed");
27634c13 331 histoName = "fHistNP";
332 if(gCentralityEstimator) histoName += gCentralityEstimator;
333 hNPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
334 histoName = "fHistNP";
335 if(gCentralityEstimator) histoName += gCentralityEstimator;
52daf7b2 336 hNPMixed->SetName("gHistNPMixed");
27634c13 337 histoName = "fHistPP";
338 if(gCentralityEstimator) histoName += gCentralityEstimator;
339 hPPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
52daf7b2 340 hPPMixed->SetName("gHistPPMixed");
27634c13 341 histoName = "fHistNN";
342 if(gCentralityEstimator) histoName += gCentralityEstimator;
343 hNNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
52daf7b2 344 hNNMixed->SetName("gHistNNMixed");
345
346 AliBalancePsi *bMixed = new AliBalancePsi();
d67eae53 347 bMixed->SetEventClass(eventClass);
52daf7b2 348 bMixed->SetHistNp(hPMixed);
349 bMixed->SetHistNn(hNMixed);
350 bMixed->SetHistNpn(hPNMixed);
351 bMixed->SetHistNnp(hNPMixed);
352 bMixed->SetHistNpp(hPPMixed);
353 bMixed->SetHistNnn(hNNMixed);
832b21d8 354 if(kUseVzBinning) bMixed->SetVertexZBinning(kTRUE);
355
52daf7b2 356 }
eb63b883 357
6acdbcb2 358 TH2D *gHistBalanceFunction;
eb63b883 359 TH2D *gHistBalanceFunctionSubtracted;
6acdbcb2 360 TH2D *gHistBalanceFunctionShuffled;
eb63b883 361 TH2D *gHistBalanceFunctionMixed;
6acdbcb2 362 TString histoTitle, pngName;
6acdbcb2 363
d67eae53 364 if(eventClass == "Centrality"){
365 histoTitle = "Centrality: ";
366 histoTitle += psiMin;
367 histoTitle += " - ";
368 histoTitle += psiMax;
369 histoTitle += " % ";
52daf7b2 370 histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
d67eae53 371 }
372 else if(eventClass == "Multiplicity"){
373 histoTitle = "Multiplicity: ";
374 histoTitle += psiMin;
375 histoTitle += " - ";
376 histoTitle += psiMax;
377 histoTitle += " tracks";
378 histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
379 }
380 else{ // "EventPlane" (default)
381 histoTitle = "Centrality: ";
382 histoTitle += centralityArray[gCentrality-1];
383 histoTitle += "%";
384 if((psiMin == -0.5)&&(psiMax == 0.5))
385 histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
386 else if((psiMin == 0.5)&&(psiMax == 1.5))
387 histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
388 else if((psiMin == 1.5)&&(psiMax == 2.5))
389 histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
390 else
391 histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
392 }
393
f0c5040c 394 if(k2pMethod)
395 if(bMixed)
34616eda 396 gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
f0c5040c 397 else{
398 cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
399 return;
400 }
401 else
34616eda 402 gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
6acdbcb2 403 gHistBalanceFunction->SetTitle(histoTitle.Data());
404 gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
eb63b883 405 gHistBalanceFunction->SetName("gHistBalanceFunction");
406
52daf7b2 407 if(listBFShuffled) {
f0c5040c 408
409 if(k2pMethod)
410 if(bMixed)
34616eda 411 gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
f0c5040c 412 else{
413 cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
414 return;
415 }
416 else
34616eda 417 gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
52daf7b2 418 gHistBalanceFunctionShuffled->SetTitle(histoTitle.Data());
419 gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
420 gHistBalanceFunctionShuffled->SetName("gHistBalanceFunctionShuffled");
421 }
eb63b883 422
52daf7b2 423 if(listBFMixed) {
f0c5040c 424 if(k2pMethod)
425 if(bMixed)
34616eda 426 gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
f0c5040c 427 else{
428 cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
429 return;
430 }
431 else
34616eda 432 gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
52daf7b2 433 gHistBalanceFunctionMixed->SetTitle(histoTitle.Data());
434 gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
435 gHistBalanceFunctionMixed->SetName("gHistBalanceFunctionMixed");
436
437 gHistBalanceFunctionSubtracted = dynamic_cast<TH2D *>(gHistBalanceFunction->Clone());
438 gHistBalanceFunctionSubtracted->Add(gHistBalanceFunctionMixed,-1);
439 gHistBalanceFunctionSubtracted->SetTitle(histoTitle.Data());
440 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
441 gHistBalanceFunctionSubtracted->SetName("gHistBalanceFunctionSubtracted");
442 }
eb63b883 443
444 //Draw the results
445 TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
6acdbcb2 446 c1->SetFillColor(10);
447 c1->SetHighLightColor(10);
448 c1->SetLeftMargin(0.15);
eb63b883 449 gHistBalanceFunction->DrawCopy("lego2");
5de9ad1a 450 gPad->SetTheta(30); // default is 30
451 gPad->SetPhi(-60); // default is 30
452 gPad->Update();
eb63b883 453 TCanvas *c1a = new TCanvas("c1a","",600,0,600,500);
454 c1a->SetFillColor(10);
455 c1a->SetHighLightColor(10);
456 c1a->SetLeftMargin(0.15);
457 gHistBalanceFunction->DrawCopy("colz");
458
52daf7b2 459 if(listBFShuffled) {
460 TCanvas *c2 = new TCanvas("c2","",100,100,600,500);
461 c2->SetFillColor(10);
462 c2->SetHighLightColor(10);
463 c2->SetLeftMargin(0.15);
464 gHistBalanceFunctionShuffled->DrawCopy("lego2");
5de9ad1a 465 gPad->SetTheta(30); // default is 30
466 gPad->SetPhi(-60); // default is 30
467 gPad->Update();
52daf7b2 468 TCanvas *c2a = new TCanvas("c2a","",700,100,600,500);
469 c2a->SetFillColor(10);
470 c2a->SetHighLightColor(10);
471 c2a->SetLeftMargin(0.15);
472 gHistBalanceFunctionShuffled->DrawCopy("colz");
473 }
eb63b883 474
52daf7b2 475 if(listBFMixed) {
476 TCanvas *c3 = new TCanvas("c3","",200,200,600,500);
477 c3->SetFillColor(10);
478 c3->SetHighLightColor(10);
479 c3->SetLeftMargin(0.15);
480 gHistBalanceFunctionMixed->DrawCopy("lego2");
5de9ad1a 481 gPad->SetTheta(30); // default is 30
482 gPad->SetPhi(-60); // default is 30
483 gPad->Update();
52daf7b2 484 TCanvas *c3a = new TCanvas("c3a","",800,200,600,500);
485 c3a->SetFillColor(10);
486 c3a->SetHighLightColor(10);
487 c3a->SetLeftMargin(0.15);
488 gHistBalanceFunctionMixed->DrawCopy("colz");
eb63b883 489
52daf7b2 490 TCanvas *c4 = new TCanvas("c4","",300,300,600,500);
491 c4->SetFillColor(10);
492 c4->SetHighLightColor(10);
493 c4->SetLeftMargin(0.15);
494 gHistBalanceFunctionSubtracted->DrawCopy("lego2");
5de9ad1a 495 gPad->SetTheta(30); // default is 30
496 gPad->SetPhi(-60); // default is 30
497 gPad->Update();
52daf7b2 498 TCanvas *c4a = new TCanvas("c4a","",900,300,600,500);
499 c4a->SetFillColor(10);
500 c4a->SetHighLightColor(10);
501 c4a->SetLeftMargin(0.15);
502 gHistBalanceFunctionSubtracted->DrawCopy("colz");
742af4bd 503
bd36d661 504 fitbalanceFunction(gCentrality, psiMin , psiMax,
742af4bd 505 ptTriggerMin, ptTriggerMax,
506 ptAssociatedMin, ptAssociatedMax,
d67eae53 507 gHistBalanceFunctionSubtracted,k2pMethod, eventClass);
52daf7b2 508 }
eb63b883 509
d67eae53 510 TString newFileName = "balanceFunction2D.";
511 if(eventClass == "Centrality"){
512 newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
513 newFileName += ".PsiAll.PttFrom";
514 }
515 else if(eventClass == "Multiplicity"){
516 newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
517 newFileName += ".PsiAll.PttFrom";
518 }
519 else{ // "EventPlane" (default)
520 newFileName += "Centrality";
521 newFileName += gCentrality; newFileName += ".Psi";
522 if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
523 else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
524 else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
525 else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
526 else newFileName += "All.PttFrom";
527 }
648f1a5a 528 newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
529 newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
530 newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
531 newFileName += Form("%.1f",ptAssociatedMax);
f0c5040c 532 if(k2pMethod) newFileName += "_2pMethod";
648f1a5a 533 newFileName += ".root";
eb63b883 534
535 TFile *fOutput = new TFile(newFileName.Data(),"recreate");
536 fOutput->cd();
52daf7b2 537 /*hP->Write(); hN->Write();
538 hPN->Write(); hNP->Write();
539 hPP->Write(); hNN->Write();
540 hPShuffled->Write(); hNShuffled->Write();
541 hPNShuffled->Write(); hNPShuffled->Write();
542 hPPShuffled->Write(); hNNShuffled->Write();
543 hPMixed->Write(); hNMixed->Write();
544 hPNMixed->Write(); hNPMixed->Write();
545 hPPMixed->Write(); hNNMixed->Write();*/
eb63b883 546 gHistBalanceFunction->Write();
52daf7b2 547 if(listBFShuffled) gHistBalanceFunctionShuffled->Write();
548 if(listBFMixed) {
549 gHistBalanceFunctionMixed->Write();
550 gHistBalanceFunctionSubtracted->Write();
551 }
eb63b883 552 fOutput->Close();
6acdbcb2 553}
554
742af4bd 555//____________________________________________________________//
556void fitbalanceFunction(Int_t gCentrality = 1,
557 Double_t psiMin = -0.5, Double_t psiMax = 3.5,
558 Double_t ptTriggerMin = -1.,
559 Double_t ptTriggerMax = -1.,
560 Double_t ptAssociatedMin = -1.,
561 Double_t ptAssociatedMax = -1.,
d67eae53 562 TH2D *gHist,
563 Bool_t k2pMethod = kFALSE,
564 TString eventClass="EventPlane") {
48f521eb 565 //balancing charges: [1]*TMath::Exp(-0.5*TMath::Power(((x - [3])/[2]),2)-0.5*TMath::Power(((y - [5])/[4]),2))
566 //short range correlations: [6]*TMath::Exp(-0.5*TMath::Power(((x - [8])/[7]),2)-0.5*TMath::Power(((y - [10])/[9]),2))
742af4bd 567 cout<<"FITTING FUNCTION"<<endl;
568
48f521eb 569 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 570 gFitFunction->SetName("gFitFunction");
571
742af4bd 572 //Normalization
573 gFitFunction->SetParName(0,"N1");
48f521eb 574 gFitFunction->SetParameter(0,1.0);
575
576 //2D balance function
577 gFitFunction->SetParName(1,"N_{BF}");
578 gFitFunction->SetParameter(1,1.0);
579 gFitFunction->SetParLimits(1, 0., 100.);
580 gFitFunction->SetParName(2,"Sigma_{BF}(delta eta)");
581 gFitFunction->SetParameter(2,0.6);
582 gFitFunction->SetParLimits(2, 0., 1.);
583 gFitFunction->SetParName(3,"Mean_{BF}(delta eta)");
584 gFitFunction->SetParameter(3,0.0);
585 gFitFunction->SetParLimits(3, -0.2, 0.2);
586 gFitFunction->SetParName(4,"Sigma_{BF}(delta phi)");
587 gFitFunction->SetParameter(4,0.6);
588 gFitFunction->SetParLimits(4, 0., 1.);
589 gFitFunction->SetParName(5,"Mean_{BF}(delta phi)");
590 gFitFunction->SetParameter(5,0.0);
591 gFitFunction->SetParLimits(5, -0.2, 0.2);
592
593 //short range structure
594 gFitFunction->SetParName(6,"N_{SR}");
595 gFitFunction->SetParameter(6,5.0);
596 gFitFunction->SetParLimits(6, 0., 100.);
597 gFitFunction->SetParName(7,"Sigma_{SR}(delta eta)");
598 gFitFunction->SetParameter(7,0.01);
599 gFitFunction->SetParLimits(7, 0.0, 0.1);
600 gFitFunction->SetParName(8,"Mean_{SR}(delta eta)");
601 gFitFunction->SetParameter(8,0.0);
602 gFitFunction->SetParLimits(8, -0.01, 0.01);
603 gFitFunction->SetParName(9,"Sigma_{SR}(delta phi)");
604 gFitFunction->SetParameter(9,0.01);
605 gFitFunction->SetParLimits(9, 0.0, 0.1);
606 gFitFunction->SetParName(10,"Mean_{SR}(delta phi)");
607 gFitFunction->SetParameter(10,0.0);
608 gFitFunction->SetParLimits(10, -0.01, 0.01);
742af4bd 609
610
611 //Cloning the histogram
742af4bd 612 TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
613 gHistResidual->SetName("gHistResidual");
614 gHistResidual->Sumw2();
615
48f521eb 616 //Fitting the 2D bf
617 for(Int_t iAttempt = 0; iAttempt < 10; iAttempt++) {
618 gHist->Fit("gFitFunction","nm");
619 for(Int_t iParam = 0; iParam < 11; iParam++)
620 gFitFunction->SetParameter(iParam,gFitFunction->GetParameter(iParam));
742af4bd 621 }
48f521eb 622 cout<<"======================================================"<<endl;
623 cout<<"Fit chi2/ndf: "<<gFitFunction->GetChisquare()/gFitFunction->GetNDF()<<" - chi2: "<<gFitFunction->GetChisquare()<<" - ndf: "<<gFitFunction->GetNDF()<<endl;
624 cout<<"======================================================"<<endl;
625
626 //Getting the residual
627 gHistResidual->Add(gFitFunction,-1);
742af4bd 628
629 //Write to output file
d67eae53 630 TString newFileName = "balanceFunctionFit2D.";
631 if(eventClass == "Centrality"){
632 newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
633 newFileName += ".PsiAll.PttFrom";
634 }
635 else if(eventClass == "Multiplicity"){
636 newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
637 newFileName += ".PsiAll.PttFrom";
638 }
639 else{ // "EventPlane" (default)
640 newFileName += "Centrality";
641 newFileName += gCentrality; newFileName += ".Psi";
642 if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
643 else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
644 else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
645 else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
646 else newFileName += "All.PttFrom";
647 }
742af4bd 648 newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
649 newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
650 newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
651 newFileName += Form("%.1f",ptAssociatedMax);
d67eae53 652 if(k2pMethod) newFileName += "_2pMethod";
742af4bd 653 newFileName += ".root";
654 TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
655 gHist->Write();
742af4bd 656 gHistResidual->Write();
657 gFitFunction->Write();
658 newFile->Close();
742af4bd 659}
660
eb63b883 661//____________________________________________________________//
648f1a5a 662void drawBFPsi2D(const char* lhcPeriod = "LHC11h",
bd36d661 663 const char* gCentralityEstimator = "V0M",
664 Int_t gBit = 128,
665 const char* gEventPlaneEstimator = "VZERO",
648f1a5a 666 Int_t gCentrality = 1,
667 Bool_t kShowShuffled = kFALSE,
668 Bool_t kShowMixed = kFALSE,
eb63b883 669 Double_t psiMin = -0.5, Double_t psiMax = 0.5,
670 Double_t ptTriggerMin = -1.,
671 Double_t ptTriggerMax = -1.,
672 Double_t ptAssociatedMin = -1.,
bd36d661 673 Double_t ptAssociatedMax = -1.,
674 Bool_t k2pMethod = kTRUE) {
eb63b883 675 //Macro that draws the BF distributions for each centrality bin
676 //for reaction plane dependent analysis
677 //Author: Panos.Christakoglou@nikhef.nl
db7174c0 678 TGaxis::SetMaxDigits(3);
eb63b883 679
680 //Get the input file
bd36d661 681 TString filename = lhcPeriod;
682 filename += "/Centrality"; filename += gCentralityEstimator;
683 filename += "_Bit"; filename += gBit;
684 filename += "_"; filename += gEventPlaneEstimator;
685 filename +="/PttFrom";
648f1a5a 686 filename += Form("%.1f",ptTriggerMin); filename += "To";
687 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
688 filename += Form("%.1f",ptAssociatedMin); filename += "To";
689 filename += Form("%.1f",ptAssociatedMax);
690 filename += "/balanceFunction2D.Centrality";
eb63b883 691 filename += gCentrality; filename += ".Psi";
692 if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
693 else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
694 else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
648f1a5a 695 else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
696 else filename += "All.PttFrom";
697 filename += Form("%.1f",ptTriggerMin); filename += "To";
db7174c0 698 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
648f1a5a 699 filename += Form("%.1f",ptAssociatedMin); filename += "To";
bd36d661 700 filename += Form("%.1f",ptAssociatedMax);
701 if(k2pMethod) filename += "_2pMethod";
702 filename += ".root";
eb63b883 703
704 //Open the file
705 TFile *f = TFile::Open(filename.Data());
706 if((!f)||(!f->IsOpen())) {
707 Printf("The file %s is not found. Aborting...",filename);
708 return listBF;
6acdbcb2 709 }
eb63b883 710 //f->ls();
711
712 //Raw balance function
713 TH1D *gHistBalanceFunction = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunction"));
714 gHistBalanceFunction->SetStats(kFALSE);
715 gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
716 gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
717 gHistBalanceFunction->GetZaxis()->SetNdivisions(10);
718 gHistBalanceFunction->GetXaxis()->SetTitleOffset(1.3);
719 gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
720 gHistBalanceFunction->GetZaxis()->SetTitleOffset(1.3);
721 gHistBalanceFunction->GetXaxis()->SetTitle("#Delta #eta");
27634c13 722 gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (rad)");
15dd45a0 723 gHistBalanceFunction->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
eb63b883 724
725 //Shuffled balance function
648f1a5a 726 if(kShowShuffled) {
727 TH1D *gHistBalanceFunctionShuffled = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionShuffled"));
728 gHistBalanceFunctionShuffled->SetStats(kFALSE);
729 gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
730 gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
731 gHistBalanceFunctionShuffled->GetZaxis()->SetNdivisions(10);
732 gHistBalanceFunctionShuffled->GetXaxis()->SetTitleOffset(1.3);
733 gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
734 gHistBalanceFunctionShuffled->GetZaxis()->SetTitleOffset(1.3);
735 gHistBalanceFunctionShuffled->GetXaxis()->SetTitle("#Delta #eta");
27634c13 736 gHistBalanceFunctionShuffled->GetYaxis()->SetTitle("#Delta #varphi (rad)");
648f1a5a 737 gHistBalanceFunctionShuffled->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
738 }
eb63b883 739
740 //Mixed balance function
648f1a5a 741 if(kShowMixed) {
742 TH1D *gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionMixed"));
743 gHistBalanceFunctionMixed->SetStats(kFALSE);
744 gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
745 gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
746 gHistBalanceFunctionMixed->GetZaxis()->SetNdivisions(10);
747 gHistBalanceFunctionMixed->GetXaxis()->SetTitleOffset(1.3);
748 gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
749 gHistBalanceFunctionMixed->GetZaxis()->SetTitleOffset(1.3);
750 gHistBalanceFunctionMixed->GetXaxis()->SetTitle("#Delta #eta");
27634c13 751 gHistBalanceFunctionMixed->GetYaxis()->SetTitle("#Delta #varphi (rad)");
648f1a5a 752 gHistBalanceFunctionMixed->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
753 }
eb63b883 754
755 //Subtracted balance function
648f1a5a 756 if(kShowMixed) {
757 TH1D *gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionSubtracted"));
758 gHistBalanceFunctionSubtracted->SetStats(kFALSE);
759 gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
760 gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
761 gHistBalanceFunctionSubtracted->GetZaxis()->SetNdivisions(10);
762 gHistBalanceFunctionSubtracted->GetXaxis()->SetTitleOffset(1.3);
763 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
764 gHistBalanceFunctionSubtracted->GetZaxis()->SetTitleOffset(1.3);
765 gHistBalanceFunctionSubtracted->GetXaxis()->SetTitle("#Delta #eta");
27634c13 766 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitle("#Delta #varphi (rad)");
648f1a5a 767 gHistBalanceFunctionSubtracted->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
768 }
769
eb63b883 770 TString pngName;
771
772 TString centralityLatex = "Centrality: ";
773 centralityLatex += centralityArray[gCentrality-1];
774 centralityLatex += "%";
775
776 TString psiLatex;
777 if((psiMin == -0.5)&&(psiMax == 0.5))
778 psiLatex = " -7.5^{o} < #varphi - #Psi_{2} < 7.5^{o}";
779 else if((psiMin == 0.5)&&(psiMax == 1.5))
780 psiLatex = " 37.5^{o} < #varphi - #Psi_{2} < 52.5^{o}";
781 else if((psiMin == 1.5)&&(psiMax == 2.5))
782 psiLatex = " 82.5^{o} < #varphi - #Psi_{2} < 97.5^{o}";
783 else
784 psiLatex = " 0^{o} < #varphi - #Psi_{2} < 180^{o}";
785
786 TString pttLatex = Form("%.1f",ptTriggerMin);
787 pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
788 pttLatex += " GeV/c";
789
790 TString ptaLatex = Form("%.1f",ptAssociatedMin);
791 ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
792 ptaLatex += " GeV/c";
793
794 TLatex *latexInfo1 = new TLatex();
795 latexInfo1->SetNDC();
15dd45a0 796 latexInfo1->SetTextSize(0.045);
eb63b883 797 latexInfo1->SetTextColor(1);
798
799 //Draw the results
800 TCanvas *c1 = new TCanvas("c1","Raw balance function 2D",0,0,600,500);
801 c1->SetFillColor(10); c1->SetHighLightColor(10);
802 c1->SetLeftMargin(0.17); c1->SetTopMargin(0.05);
bd36d661 803 gHistBalanceFunction->SetTitle("");
eb63b883 804 gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.4);
805 gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
27634c13 806 gHistBalanceFunction->GetXaxis()->SetRangeUser(-1.4,1.4);
eb63b883 807 gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
27634c13 808 gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (rad)");
eb63b883 809 gHistBalanceFunction->DrawCopy("lego2");
5de9ad1a 810 gPad->SetTheta(30); // default is 30
811 gPad->SetPhi(-60); // default is 30
812 gPad->Update();
eb63b883 813
15dd45a0 814 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
815 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
816 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
817 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
eb63b883 818
bd36d661 819 TString pngName = "BalanceFunction2D.";
820 pngName += "Centrality";
821 pngName += gCentrality;
822 if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
823 else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
824 else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
825 else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.PttFrom";
826 else pngName += "All.PttFrom";
827 pngName += Form("%.1f",ptTriggerMin); pngName += "To";
828 pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
829 pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
830 pngName += Form("%.1f",ptAssociatedMax);
831 if(k2pMethod) pngName += "_2pMethod";
832 pngName += ".png";
833
834 c1->SaveAs(pngName.Data());
835
648f1a5a 836 if(kShowShuffled) {
837 TCanvas *c2 = new TCanvas("c2","Shuffled balance function 2D",100,100,600,500);
838 c2->SetFillColor(10); c2->SetHighLightColor(10);
839 c2->SetLeftMargin(0.17); c2->SetTopMargin(0.05);
840 gHistBalanceFunctionShuffled->SetTitle("Shuffled events");
841 gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.4);
842 gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
843 gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
844 gHistBalanceFunctionShuffled->DrawCopy("lego2");
5de9ad1a 845 gPad->SetTheta(30); // default is 30
846 gPad->SetPhi(-60); // default is 30
847 gPad->Update();
848
648f1a5a 849 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
850 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
851 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
852 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
853 }
eb63b883 854
648f1a5a 855 if(kShowMixed) {
856 TCanvas *c3 = new TCanvas("c3","Mixed balance function 2D",200,200,600,500);
857 c3->SetFillColor(10); c3->SetHighLightColor(10);
858 c3->SetLeftMargin(0.17); c3->SetTopMargin(0.05);
859 gHistBalanceFunctionMixed->SetTitle("Mixed events");
860 gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.4);
861 gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
862 gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
863 gHistBalanceFunctionMixed->DrawCopy("lego2");
5de9ad1a 864 gPad->SetTheta(30); // default is 30
865 gPad->SetPhi(-60); // default is 30
866 gPad->Update();
867
648f1a5a 868 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
869 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
870 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
871 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
872 }
eb63b883 873
648f1a5a 874 if(kShowMixed) {
875 TCanvas *c4 = new TCanvas("c4","Subtracted balance function 2D",300,300,600,500);
876 c4->SetFillColor(10); c4->SetHighLightColor(10);
877 c4->SetLeftMargin(0.17); c4->SetTopMargin(0.05);
878 gHistBalanceFunctionSubtracted->SetTitle("Subtracted balance function");
879 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.4);
880 gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
881 gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
882 gHistBalanceFunctionSubtracted->DrawCopy("lego2");
5de9ad1a 883 gPad->SetTheta(30); // default is 30
884 gPad->SetPhi(-60); // default is 30
885 gPad->Update();
886
648f1a5a 887 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
888 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
889 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
890 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
891 }
6acdbcb2 892}
17c69b4e 893
894//____________________________________________________________//
895void drawProjections(const char* lhcPeriod = "LHC10h",
896 const char* gCentralityEstimator = "V0M",
897 Int_t gBit = 128,
898 const char* gEventPlaneEstimator = "VZERO",
899 Bool_t kProjectInEta = kFALSE,
900 Int_t binMin = 1,
901 Int_t binMax = 80,
902 Int_t gCentrality = 1,
903 Double_t psiMin = -0.5,
904 Double_t psiMax = 3.5,
905 Double_t vertexZMin = -10.,
906 Double_t vertexZMax = 10.,
907 Double_t ptTriggerMin = -1.,
908 Double_t ptTriggerMax = -1.,
909 Double_t ptAssociatedMin = -1.,
910 Double_t ptAssociatedMax = -1.,
911 Bool_t kUseZYAM = kFALSE,
912 Bool_t k2pMethod = kTRUE,
913 TString eventClass = "Centrality",
914 Bool_t bRootMoments = kFALSE) {
915 //Macro that draws the charge dependent correlation functions PROJECTIONS
916 //for each centrality bin for the different pT of trigger and
917 //associated particles
918 TGaxis::SetMaxDigits(3);
919
920 //first we need some libraries
921 gSystem->Load("libANALYSIS.so");
922 gSystem->Load("libANALYSISalice.so");
923 gSystem->Load("libEventMixing.so");
924 gSystem->Load("libCORRFW.so");
925 gSystem->Load("libPWGTools.so");
926 gSystem->Load("libPWGCFebye.so");
927
928 //Get the input file
929 TString filename = "balanceFunction2D.";
930 if(eventClass == "Centrality"){
931 filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
932 filename += ".PsiAll.PttFrom";
933 }
934 else if(eventClass == "Multiplicity"){
935 filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
936 filename += ".PsiAll.PttFrom";
937 }
938 else{ // "EventPlane" (default)
939 filename += "Centrality";
940 filename += gCentrality; filename += ".Psi";
941 if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
942 else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
943 else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
944 else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
945 else filename += "All.PttFrom";
946 }
947 filename += Form("%.1f",ptTriggerMin); filename += "To";
948 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
949 filename += Form("%.1f",ptAssociatedMin); filename += "To";
950 filename += Form("%.1f",ptAssociatedMax);
951 if(k2pMethod) filename += "_2pMethod";
952 filename += ".root";
953
954 //Open the file
955 TFile *f = TFile::Open(filename.Data());
956 if((!f)||(!f->IsOpen())) {
957 Printf("The file %s is not found. Aborting...",filename);
958 return listBF;
959 }
960 //f->ls();
961
962 //Latex
963 TString centralityLatex = "Centrality: ";
964 centralityLatex += centralityArray[gCentrality-1];
965 centralityLatex += "%";
966
967 TString psiLatex;
968 if((psiMin == -0.5)&&(psiMax == 0.5))
969 psiLatex = " -7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o}";
970 else if((psiMin == 0.5)&&(psiMax == 1.5))
971 psiLatex = " 37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o}";
972 else if((psiMin == 1.5)&&(psiMax == 2.5))
973 psiLatex = " 82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o}";
974 else
975 psiLatex = " 0^{o} < #varphi^{t} - #Psi_{2} < 180^{o}";
976
977 TString pttLatex = Form("%.1f",ptTriggerMin);
978 pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
979 pttLatex += " GeV/c";
980
981 TString ptaLatex = Form("%.1f",ptAssociatedMin);
982 ptaLatex += " < p_{T,assoc} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
983 ptaLatex += " GeV/c";
984
985 TLatex *latexInfo1 = new TLatex();
986 latexInfo1->SetNDC();
987 latexInfo1->SetTextSize(0.045);
988 latexInfo1->SetTextColor(1);
989
990 TString pngName;
991
992 //============================================================//
993 //Get subtracted and mixed balance function
994
995 TH2D *gHistBalanceFunctionSubtracted2D = (TH2D*)f->Get("gHistBalanceFunctionSubtracted");
996 TH2D *gHistBalanceFunctionMixed2D = (TH2D*)f->Get("gHistBalanceFunctionMixed");
997
998 TH1D *gHistBalanceFunctionSubtracted = NULL;
999 TH1D *gHistBalanceFunctionMixed = NULL;
1000
1001 if(kProjectInEta){
1002 gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionX());
1003 gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionX());
1004 gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#eta)");
1005 gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#eta)");
1006 }
1007 else{
1008 gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionY());
1009 gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionY());
1010 gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#varphi)");
1011 gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#varphi)");
1012 }
1013
1014 gHistBalanceFunctionSubtracted->SetMarkerStyle(20);
1015 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
1016 gHistBalanceFunctionSubtracted->SetName("gHistBalanceFunctionSubtracted");
1017
1018 gHistBalanceFunctionMixed->SetMarkerStyle(25);
1019 gHistBalanceFunctionMixed->SetName("gHistBalanceFunctionMixed");
1020
1021 TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
1022 c1->SetFillColor(10);
1023 c1->SetHighLightColor(10);
1024 c1->SetLeftMargin(0.15);
1025 gHistBalanceFunctionSubtracted->DrawCopy("E");
1026 gHistBalanceFunctionMixed->DrawCopy("E, SAME");
1027
1028 legend = new TLegend(0.18,0.62,0.45,0.82,"","brNDC");
1029 legend->SetTextSize(0.045);
1030 legend->SetTextFont(42);
1031 legend->SetBorderSize(0);
1032 legend->SetFillStyle(0);
1033 legend->SetFillColor(10);
1034 legend->SetMargin(0.25);
1035 legend->SetShadowColor(10);
1036 legend->AddEntry(gHistBalanceFunctionSubtracted,"Data","lp");
1037 legend->AddEntry(gHistBalanceFunctionMixed,"Mixed data","lp");
1038 legend->Draw();
1039
1040 pngName = "BalanceFunction.";
1041 if(eventClass == "Centrality"){
1042 pngName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1043 if(kProjectInEta) pngName += ".InDeltaEta.PsiAll.PttFrom";
1044 else pngName += ".InDeltaPhi.PsiAll.PttFrom";
1045 }
1046 else if(eventClass == "Multiplicity"){
1047 pngName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1048 if(kProjectInEta) pngName += ".InDeltaEta.PsiAll.PttFrom";
1049 else pngName += ".InDeltaPhi.PsiAll.PttFrom";
1050 }
1051 else{ // "EventPlane" (default)
1052 pngName += "Centrality";
1053 pngName += gCentrality;
1054 if(kProjectInEta) pngName += ".InDeltaEta.Psi";
1055 else pngName += ".InDeltaPhi.Psi";
1056 if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
1057 else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
1058 else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
1059 else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.PttFrom";
1060 else pngName += "All.PttFrom";
1061 }
1062 pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1063 pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1064 pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1065 pngName += Form("%.1f",ptAssociatedMax);
1066 if(k2pMethod) pngName += "_2pMethod";
1067 pngName += ".png";
1068
1069 c1->SaveAs(pngName.Data());
1070
1071 TString meanLatex, rmsLatex, skewnessLatex, kurtosisLatex;
1072
1073 if(bRootMoments){
1074 meanLatex = "#mu = ";
1075 meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMean());
1076 meanLatex += " #pm ";
1077 meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMeanError());
1078
1079 rmsLatex = "#sigma = ";
1080 rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMS());
1081 rmsLatex += " #pm ";
1082 rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMSError());
1083
1084 skewnessLatex = "S = ";
1085 skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(1));
1086 skewnessLatex += " #pm ";
1087 skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(11));
1088
1089 kurtosisLatex = "K = ";
1090 kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(1));
1091 kurtosisLatex += " #pm ";
1092 kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(11));
1093 Printf("Mean: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetMean(),gHistBalanceFunctionSubtracted->GetMeanError());
1094 Printf("RMS: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetRMS(),gHistBalanceFunctionSubtracted->GetRMSError());
1095 Printf("Skeweness: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetSkewness(1),gHistBalanceFunctionSubtracted->GetSkewness(11));
1096 Printf("Kurtosis: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetKurtosis(1),gHistBalanceFunctionSubtracted->GetKurtosis(11));
1097 }
1098 // calculate the moments by hand
1099 else{
1100
1101 Double_t meanAnalytical, meanAnalyticalError;
1102 Double_t sigmaAnalytical, sigmaAnalyticalError;
1103 Double_t skewnessAnalytical, skewnessAnalyticalError;
1104 Double_t kurtosisAnalytical, kurtosisAnalyticalError;
1105
1106 Int_t gDeltaEtaPhi = 2;
1107 if(kProjectInEta) gDeltaEtaPhi = 1;
1108
1109 AliBalancePsi *bHelper = new AliBalancePsi;
1110 bHelper->GetMomentsAnalytical(gDeltaEtaPhi,gHistBalanceFunctionSubtracted,meanAnalytical,meanAnalyticalError,sigmaAnalytical,sigmaAnalyticalError,skewnessAnalytical,skewnessAnalyticalError,kurtosisAnalytical,kurtosisAnalyticalError);
1111
1112 meanLatex = "#mu = ";
1113 meanLatex += Form("%.3f",meanAnalytical);
1114 meanLatex += " #pm ";
1115 meanLatex += Form("%.3f",meanAnalyticalError);
1116
1117 rmsLatex = "#sigma = ";
1118 rmsLatex += Form("%.3f",sigmaAnalytical);
1119 rmsLatex += " #pm ";
1120 rmsLatex += Form("%.3f",sigmaAnalyticalError);
1121
1122 skewnessLatex = "S = ";
1123 skewnessLatex += Form("%.3f",skewnessAnalytical);
1124 skewnessLatex += " #pm ";
1125 skewnessLatex += Form("%.3f",skewnessAnalyticalError);
1126
1127 kurtosisLatex = "K = ";
1128 kurtosisLatex += Form("%.3f",kurtosisAnalytical);
1129 kurtosisLatex += " #pm ";
1130 kurtosisLatex += Form("%.3f",kurtosisAnalyticalError);
1131 Printf("Mean: %lf - Error: %lf",meanAnalytical, meanAnalyticalError);
1132 Printf("Sigma: %lf - Error: %lf",sigmaAnalytical, sigmaAnalyticalError);
1133 Printf("Skeweness: %lf - Error: %lf",skewnessAnalytical, skewnessAnalyticalError);
1134 Printf("Kurtosis: %lf - Error: %lf",kurtosisAnalytical, kurtosisAnalyticalError);
1135 }
1136
1137 TCanvas *c2 = new TCanvas("c2","",600,0,600,500);
1138 c2->SetFillColor(10);
1139 c2->SetHighLightColor(10);
1140 c2->SetLeftMargin(0.15);
1141 gHistBalanceFunctionSubtracted->DrawCopy("E");
1142
1143 TLatex *latex = new TLatex();
1144 latex->SetNDC();
1145 latex->SetTextSize(0.035);
1146 latex->SetTextColor(1);
1147 latex->DrawLatex(0.64,0.85,meanLatex.Data());
1148 latex->DrawLatex(0.64,0.81,rmsLatex.Data());
1149 latex->DrawLatex(0.64,0.77,skewnessLatex.Data());
1150 latex->DrawLatex(0.64,0.73,kurtosisLatex.Data());
1151
1152 TString outFileName = filename;
1153 if(kProjectInEta) outFileName.ReplaceAll(".root","_DeltaEtaProjection.root");
1154 else outFileName.ReplaceAll(".root","_DeltaPhiProjection.root");
1155 TFile *fProjection = TFile::Open(outFileName.Data(),"recreate");
1156 gHistBalanceFunctionSubtracted->Write();
1157 gHistBalanceFunctionMixed->Write();
1158 fProjection->Close();
1159}