]> 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
Implemented projections for per trigger yield (Alis)
[u/mrichter/AliRoot.git] / PWGCF / EBYE / macros / drawBalanceFunction2DPsi.C
CommitLineData
d67eae53 1const Int_t numberOfCentralityBins = 11;
2TString centralityArray[numberOfCentralityBins] = {"0-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-1","1-2","0-100"};
6acdbcb2 3
4const Int_t gRebin = 1;
5void drawBalanceFunction2DPsi(const char* filename = "AnalysisResultsPsi.root",
eb63b883 6 Int_t gCentrality = 1,
5de9ad1a 7 Int_t gBit = -1,
8 const char* gCentralityEstimator = 0x0,
52daf7b2 9 Bool_t kShowShuffled = kFALSE,
10 Bool_t kShowMixed = kTRUE,
eb63b883 11 Double_t psiMin = -0.5, Double_t psiMax = 0.5,
6acdbcb2 12 Double_t ptTriggerMin = -1.,
13 Double_t ptTriggerMax = -1.,
14 Double_t ptAssociatedMin = -1.,
f0c5040c 15 Double_t ptAssociatedMax = -1.,
d67eae53 16 Bool_t k2pMethod = kFALSE,
17 TString eventClass = "EventPlane") //Can be "EventPlane", "Centrality", "Multiplicity"
18{
6acdbcb2 19 //Macro that draws the BF distributions for each centrality bin
20 //for reaction plane dependent analysis
21 //Author: Panos.Christakoglou@nikhef.nl
22 //Load the PWG2ebye library
23 gSystem->Load("libANALYSIS.so");
24 gSystem->Load("libANALYSISalice.so");
25 gSystem->Load("libEventMixing.so");
26 gSystem->Load("libCORRFW.so");
27 gSystem->Load("libPWGTools.so");
28 gSystem->Load("libPWGCFebye.so");
29
15dd45a0 30 gROOT->LoadMacro("~/SetPlotStyle.C");
31 SetPlotStyle();
8795e993 32 gStyle->SetPalette(1,0);
15dd45a0 33
6acdbcb2 34 //Prepare the objects and return them
5de9ad1a 35 TList *listBF = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,0);
52daf7b2 36 TList *listBFShuffled = NULL;
5de9ad1a 37 if(kShowShuffled) listBFShuffled = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,1);
52daf7b2 38 TList *listBFMixed = NULL;
5de9ad1a 39 if(kShowMixed) listBFMixed = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,2);
6acdbcb2 40 if(!listBF) {
41 Printf("The TList object was not created");
42 return;
43 }
44 else
eb63b883 45 draw(listBF,listBFShuffled,listBFMixed,gCentrality,psiMin,psiMax,
f0c5040c 46 ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,
d67eae53 47 k2pMethod,eventClass);
6acdbcb2 48}
49
50//______________________________________________________//
5365d1d7 51TList *GetListOfObjects(const char* filename,
52 Int_t gCentrality,
53 Int_t gBit,
54 const char *gCentralityEstimator,
5de9ad1a 55 Int_t kData = 1) {
6acdbcb2 56 //Get the TList objects (QA, bf, bf shuffled)
6acdbcb2 57 TList *listBF = 0x0;
6acdbcb2 58
59 //Open the file
5365d1d7 60 TFile *f = TFile::Open(filename,"UPDATE");
6acdbcb2 61 if((!f)||(!f->IsOpen())) {
62 Printf("The file %s is not found. Aborting...",filename);
63 return listBF;
64 }
65 //f->ls();
66
67 TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
68 if(!dir) {
69 Printf("The TDirectoryFile is not found. Aborting...",filename);
70 return listBF;
71 }
72 //dir->ls();
73
74 TString listBFName;
eb63b883 75 if(kData == 0) {
76 //cout<<"no shuffling - no mixing"<<endl;
77 listBFName = "listBFPsi_";
78 }
79 else if(kData == 1) {
80 //cout<<"shuffling - no mixing"<<endl;
81 listBFName = "listBFPsiShuffled_";
82 }
83 else if(kData == 2) {
84 //cout<<"no shuffling - mixing"<<endl;
85 listBFName = "listBFPsiMixed_";
86 }
87 listBFName += centralityArray[gCentrality-1];
5de9ad1a 88 if(gBit > -1) {
89 listBFName += "_Bit"; listBFName += gBit; }
90 if(gCentralityEstimator) {
91 listBFName += "_"; listBFName += gCentralityEstimator;}
6acdbcb2 92
5365d1d7 93 // histograms were already retrieved (in first iteration)
94 if(dir->Get(Form("%s_histograms",listBFName.Data()))){
95 //cout<<"second iteration"<<endl;
96 listBF = dynamic_cast<TList *>(dir->Get(Form("%s_histograms",listBFName.Data())));
6acdbcb2 97 }
6acdbcb2 98
5365d1d7 99 // histograms were not yet retrieved (this is the first iteration)
100 else{
101
102 listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
103 cout<<"======================================================="<<endl;
104 cout<<"List name (check): "<<listBFName.Data()<<endl;
105 cout<<"List name: "<<listBF->GetName()<<endl;
106 //listBF->ls();
6acdbcb2 107
5365d1d7 108 //Get the histograms
109 TString histoName;
110 if(kData == 0)
111 histoName = "fHistPV0M";
112 else if(kData == 1)
113 histoName = "fHistP_shuffleV0M";
114 else if(kData == 2)
115 histoName = "fHistPV0M";
116 AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
117 if(!fHistP) {
118 Printf("fHistP %s not found!!!",histoName.Data());
119 break;
120 }
121 fHistP->FillParent(); fHistP->DeleteContainers();
122
123 if(kData == 0)
124 histoName = "fHistNV0M";
125 if(kData == 1)
126 histoName = "fHistN_shuffleV0M";
127 if(kData == 2)
128 histoName = "fHistNV0M";
129 AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
130 if(!fHistN) {
131 Printf("fHistN %s not found!!!",histoName.Data());
132 break;
133 }
134 fHistN->FillParent(); fHistN->DeleteContainers();
135
136 if(kData == 0)
137 histoName = "fHistPNV0M";
138 if(kData == 1)
139 histoName = "fHistPN_shuffleV0M";
140 if(kData == 2)
141 histoName = "fHistPNV0M";
142 AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
143 if(!fHistPN) {
144 Printf("fHistPN %s not found!!!",histoName.Data());
145 break;
146 }
147 fHistPN->FillParent(); fHistPN->DeleteContainers();
148
149 if(kData == 0)
150 histoName = "fHistNPV0M";
151 if(kData == 1)
152 histoName = "fHistNP_shuffleV0M";
153 if(kData == 2)
154 histoName = "fHistNPV0M";
155 AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
156 if(!fHistNP) {
157 Printf("fHistNP %s not found!!!",histoName.Data());
158 break;
159 }
160 fHistNP->FillParent(); fHistNP->DeleteContainers();
161
162 if(kData == 0)
163 histoName = "fHistPPV0M";
164 if(kData == 1)
165 histoName = "fHistPP_shuffleV0M";
166 if(kData == 2)
167 histoName = "fHistPPV0M";
168 AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
169 if(!fHistPP) {
170 Printf("fHistPP %s not found!!!",histoName.Data());
171 break;
172 }
173 fHistPP->FillParent(); fHistPP->DeleteContainers();
174
175 if(kData == 0)
176 histoName = "fHistNNV0M";
177 if(kData == 1)
178 histoName = "fHistNN_shuffleV0M";
179 if(kData == 2)
180 histoName = "fHistNNV0M";
181 AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
182 if(!fHistNN) {
183 Printf("fHistNN %s not found!!!",histoName.Data());
184 break;
185 }
186 fHistNN->FillParent(); fHistNN->DeleteContainers();
187
188 dir->cd();
189 listBF->Write(Form("%s_histograms",listBFName.Data()), TObject::kSingleKey);
190
191 }// first iteration
6acdbcb2 192
5365d1d7 193 f->Close();
6acdbcb2 194
195 return listBF;
196}
197
198//______________________________________________________//
eb63b883 199void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
200 Int_t gCentrality, Double_t psiMin, Double_t psiMax,
6acdbcb2 201 Double_t ptTriggerMin, Double_t ptTriggerMax,
f0c5040c 202 Double_t ptAssociatedMin, Double_t ptAssociatedMax,
d67eae53 203 Bool_t k2pMethod = kFALSE, TString eventClass) {
6acdbcb2 204 //balance function
205 AliTHn *hP = NULL;
206 AliTHn *hN = NULL;
207 AliTHn *hPN = NULL;
208 AliTHn *hNP = NULL;
209 AliTHn *hPP = NULL;
210 AliTHn *hNN = NULL;
211 //listBF->ls();
212 //Printf("=================");
213 hP = (AliTHn*) listBF->FindObject("fHistPV0M");
52daf7b2 214 hP->SetName("gHistP");
6acdbcb2 215 hN = (AliTHn*) listBF->FindObject("fHistNV0M");
52daf7b2 216 hN->SetName("gHistN");
6acdbcb2 217 hPN = (AliTHn*) listBF->FindObject("fHistPNV0M");
52daf7b2 218 hPN->SetName("gHistPN");
6acdbcb2 219 hNP = (AliTHn*) listBF->FindObject("fHistNPV0M");
52daf7b2 220 hNP->SetName("gHistNP");
6acdbcb2 221 hPP = (AliTHn*) listBF->FindObject("fHistPPV0M");
52daf7b2 222 hPP->SetName("gHistPP");
6acdbcb2 223 hNN = (AliTHn*) listBF->FindObject("fHistNNV0M");
52daf7b2 224 hNN->SetName("gHistNN");
6acdbcb2 225
226 AliBalancePsi *b = new AliBalancePsi();
d67eae53 227 b->SetEventClass(eventClass);
6acdbcb2 228 b->SetHistNp(hP);
229 b->SetHistNn(hN);
230 b->SetHistNpn(hPN);
231 b->SetHistNnp(hNP);
232 b->SetHistNpp(hPP);
233 b->SetHistNnn(hNN);
234
235 //balance function shuffling
236 AliTHn *hPShuffled = NULL;
237 AliTHn *hNShuffled = NULL;
238 AliTHn *hPNShuffled = NULL;
239 AliTHn *hNPShuffled = NULL;
240 AliTHn *hPPShuffled = NULL;
241 AliTHn *hNNShuffled = NULL;
52daf7b2 242 if(listBFShuffled) {
243 //listBFShuffled->ls();
244
245 hPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistP_shuffleV0M");
246 hPShuffled->SetName("gHistPShuffled");
247 hNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistN_shuffleV0M");
248 hNShuffled->SetName("gHistNShuffled");
249 hPNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPN_shuffleV0M");
250 hPNShuffled->SetName("gHistPNShuffled");
251 hNPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNP_shuffleV0M");
252 hNPShuffled->SetName("gHistNPShuffled");
253 hPPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPP_shuffleV0M");
254 hPPShuffled->SetName("gHistPPShuffled");
255 hNNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNN_shuffleV0M");
256 hNNShuffled->SetName("gHistNNShuffled");
257
258 AliBalancePsi *bShuffled = new AliBalancePsi();
d67eae53 259 bShuffled->SetEventClass(eventClass);
52daf7b2 260 bShuffled->SetHistNp(hPShuffled);
261 bShuffled->SetHistNn(hNShuffled);
262 bShuffled->SetHistNpn(hPNShuffled);
263 bShuffled->SetHistNnp(hNPShuffled);
264 bShuffled->SetHistNpp(hPPShuffled);
265 bShuffled->SetHistNnn(hNNShuffled);
266 }
6acdbcb2 267
eb63b883 268 //balance function mixing
269 AliTHn *hPMixed = NULL;
270 AliTHn *hNMixed = NULL;
271 AliTHn *hPNMixed = NULL;
272 AliTHn *hNPMixed = NULL;
273 AliTHn *hPPMixed = NULL;
274 AliTHn *hNNMixed = NULL;
eb63b883 275
52daf7b2 276 if(listBFMixed) {
277 //listBFMixed->ls();
278
279 hPMixed = (AliTHn*) listBFMixed->FindObject("fHistPV0M");
280 hPMixed->SetName("gHistPMixed");
281 hNMixed = (AliTHn*) listBFMixed->FindObject("fHistNV0M");
282 hNMixed->SetName("gHistNMixed");
283 hPNMixed = (AliTHn*) listBFMixed->FindObject("fHistPNV0M");
284 hPNMixed->SetName("gHistPNMixed");
285 hNPMixed = (AliTHn*) listBFMixed->FindObject("fHistNPV0M");
286 hNPMixed->SetName("gHistNPMixed");
287 hPPMixed = (AliTHn*) listBFMixed->FindObject("fHistPPV0M");
288 hPPMixed->SetName("gHistPPMixed");
289 hNNMixed = (AliTHn*) listBFMixed->FindObject("fHistNNV0M");
290 hNNMixed->SetName("gHistNNMixed");
291
292 AliBalancePsi *bMixed = new AliBalancePsi();
d67eae53 293 bMixed->SetEventClass(eventClass);
52daf7b2 294 bMixed->SetHistNp(hPMixed);
295 bMixed->SetHistNn(hNMixed);
296 bMixed->SetHistNpn(hPNMixed);
297 bMixed->SetHistNnp(hNPMixed);
298 bMixed->SetHistNpp(hPPMixed);
299 bMixed->SetHistNnn(hNNMixed);
300 }
eb63b883 301
6acdbcb2 302 TH2D *gHistBalanceFunction;
eb63b883 303 TH2D *gHistBalanceFunctionSubtracted;
6acdbcb2 304 TH2D *gHistBalanceFunctionShuffled;
eb63b883 305 TH2D *gHistBalanceFunctionMixed;
6acdbcb2 306 TString histoTitle, pngName;
6acdbcb2 307
d67eae53 308 if(eventClass == "Centrality"){
309 histoTitle = "Centrality: ";
310 histoTitle += psiMin;
311 histoTitle += " - ";
312 histoTitle += psiMax;
313 histoTitle += " % ";
52daf7b2 314 histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
d67eae53 315 }
316 else if(eventClass == "Multiplicity"){
317 histoTitle = "Multiplicity: ";
318 histoTitle += psiMin;
319 histoTitle += " - ";
320 histoTitle += psiMax;
321 histoTitle += " tracks";
322 histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
323 }
324 else{ // "EventPlane" (default)
325 histoTitle = "Centrality: ";
326 histoTitle += centralityArray[gCentrality-1];
327 histoTitle += "%";
328 if((psiMin == -0.5)&&(psiMax == 0.5))
329 histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
330 else if((psiMin == 0.5)&&(psiMax == 1.5))
331 histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
332 else if((psiMin == 1.5)&&(psiMax == 2.5))
333 histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
334 else
335 histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
336 }
337
f0c5040c 338 if(k2pMethod)
339 if(bMixed)
340 gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
341 else{
342 cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
343 return;
344 }
345 else
346 gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
6acdbcb2 347 gHistBalanceFunction->SetTitle(histoTitle.Data());
348 gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
eb63b883 349 gHistBalanceFunction->SetName("gHistBalanceFunction");
350
52daf7b2 351 if(listBFShuffled) {
f0c5040c 352
353 if(k2pMethod)
354 if(bMixed)
355 gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
356 else{
357 cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
358 return;
359 }
360 else
361 gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
52daf7b2 362 gHistBalanceFunctionShuffled->SetTitle(histoTitle.Data());
363 gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
364 gHistBalanceFunctionShuffled->SetName("gHistBalanceFunctionShuffled");
365 }
eb63b883 366
52daf7b2 367 if(listBFMixed) {
f0c5040c 368 if(k2pMethod)
369 if(bMixed)
370 gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
371 else{
372 cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
373 return;
374 }
375 else
376 gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
52daf7b2 377 gHistBalanceFunctionMixed->SetTitle(histoTitle.Data());
378 gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
379 gHistBalanceFunctionMixed->SetName("gHistBalanceFunctionMixed");
380
381 gHistBalanceFunctionSubtracted = dynamic_cast<TH2D *>(gHistBalanceFunction->Clone());
382 gHistBalanceFunctionSubtracted->Add(gHistBalanceFunctionMixed,-1);
383 gHistBalanceFunctionSubtracted->SetTitle(histoTitle.Data());
384 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
385 gHistBalanceFunctionSubtracted->SetName("gHistBalanceFunctionSubtracted");
386 }
eb63b883 387
388 //Draw the results
389 TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
6acdbcb2 390 c1->SetFillColor(10);
391 c1->SetHighLightColor(10);
392 c1->SetLeftMargin(0.15);
eb63b883 393 gHistBalanceFunction->DrawCopy("lego2");
5de9ad1a 394 gPad->SetTheta(30); // default is 30
395 gPad->SetPhi(-60); // default is 30
396 gPad->Update();
eb63b883 397 TCanvas *c1a = new TCanvas("c1a","",600,0,600,500);
398 c1a->SetFillColor(10);
399 c1a->SetHighLightColor(10);
400 c1a->SetLeftMargin(0.15);
401 gHistBalanceFunction->DrawCopy("colz");
402
52daf7b2 403 if(listBFShuffled) {
404 TCanvas *c2 = new TCanvas("c2","",100,100,600,500);
405 c2->SetFillColor(10);
406 c2->SetHighLightColor(10);
407 c2->SetLeftMargin(0.15);
408 gHistBalanceFunctionShuffled->DrawCopy("lego2");
5de9ad1a 409 gPad->SetTheta(30); // default is 30
410 gPad->SetPhi(-60); // default is 30
411 gPad->Update();
52daf7b2 412 TCanvas *c2a = new TCanvas("c2a","",700,100,600,500);
413 c2a->SetFillColor(10);
414 c2a->SetHighLightColor(10);
415 c2a->SetLeftMargin(0.15);
416 gHistBalanceFunctionShuffled->DrawCopy("colz");
417 }
eb63b883 418
52daf7b2 419 if(listBFMixed) {
420 TCanvas *c3 = new TCanvas("c3","",200,200,600,500);
421 c3->SetFillColor(10);
422 c3->SetHighLightColor(10);
423 c3->SetLeftMargin(0.15);
424 gHistBalanceFunctionMixed->DrawCopy("lego2");
5de9ad1a 425 gPad->SetTheta(30); // default is 30
426 gPad->SetPhi(-60); // default is 30
427 gPad->Update();
52daf7b2 428 TCanvas *c3a = new TCanvas("c3a","",800,200,600,500);
429 c3a->SetFillColor(10);
430 c3a->SetHighLightColor(10);
431 c3a->SetLeftMargin(0.15);
432 gHistBalanceFunctionMixed->DrawCopy("colz");
eb63b883 433
52daf7b2 434 TCanvas *c4 = new TCanvas("c4","",300,300,600,500);
435 c4->SetFillColor(10);
436 c4->SetHighLightColor(10);
437 c4->SetLeftMargin(0.15);
438 gHistBalanceFunctionSubtracted->DrawCopy("lego2");
5de9ad1a 439 gPad->SetTheta(30); // default is 30
440 gPad->SetPhi(-60); // default is 30
441 gPad->Update();
52daf7b2 442 TCanvas *c4a = new TCanvas("c4a","",900,300,600,500);
443 c4a->SetFillColor(10);
444 c4a->SetHighLightColor(10);
445 c4a->SetLeftMargin(0.15);
446 gHistBalanceFunctionSubtracted->DrawCopy("colz");
742af4bd 447
448 fitbalanceFunction(gCentrality, psiMin = -0.5, psiMax,
449 ptTriggerMin, ptTriggerMax,
450 ptAssociatedMin, ptAssociatedMax,
d67eae53 451 gHistBalanceFunctionSubtracted,k2pMethod, eventClass);
52daf7b2 452 }
eb63b883 453
d67eae53 454 TString newFileName = "balanceFunction2D.";
455 if(eventClass == "Centrality"){
456 newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
457 newFileName += ".PsiAll.PttFrom";
458 }
459 else if(eventClass == "Multiplicity"){
460 newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
461 newFileName += ".PsiAll.PttFrom";
462 }
463 else{ // "EventPlane" (default)
464 newFileName += "Centrality";
465 newFileName += gCentrality; newFileName += ".Psi";
466 if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
467 else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
468 else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
469 else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
470 else newFileName += "All.PttFrom";
471 }
648f1a5a 472 newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
473 newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
474 newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
475 newFileName += Form("%.1f",ptAssociatedMax);
f0c5040c 476 if(k2pMethod) newFileName += "_2pMethod";
648f1a5a 477 newFileName += ".root";
eb63b883 478
479 TFile *fOutput = new TFile(newFileName.Data(),"recreate");
480 fOutput->cd();
52daf7b2 481 /*hP->Write(); hN->Write();
482 hPN->Write(); hNP->Write();
483 hPP->Write(); hNN->Write();
484 hPShuffled->Write(); hNShuffled->Write();
485 hPNShuffled->Write(); hNPShuffled->Write();
486 hPPShuffled->Write(); hNNShuffled->Write();
487 hPMixed->Write(); hNMixed->Write();
488 hPNMixed->Write(); hNPMixed->Write();
489 hPPMixed->Write(); hNNMixed->Write();*/
eb63b883 490 gHistBalanceFunction->Write();
52daf7b2 491 if(listBFShuffled) gHistBalanceFunctionShuffled->Write();
492 if(listBFMixed) {
493 gHistBalanceFunctionMixed->Write();
494 gHistBalanceFunctionSubtracted->Write();
495 }
eb63b883 496 fOutput->Close();
6acdbcb2 497}
498
742af4bd 499//____________________________________________________________//
500void fitbalanceFunction(Int_t gCentrality = 1,
501 Double_t psiMin = -0.5, Double_t psiMax = 3.5,
502 Double_t ptTriggerMin = -1.,
503 Double_t ptTriggerMax = -1.,
504 Double_t ptAssociatedMin = -1.,
505 Double_t ptAssociatedMax = -1.,
d67eae53 506 TH2D *gHist,
507 Bool_t k2pMethod = kFALSE,
508 TString eventClass="EventPlane") {
742af4bd 509
510 cout<<"FITTING FUNCTION"<<endl;
511
512 //near side peak: [1]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4]))
513 //away side ridge: [5]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((y-TMath::Pi())/[6]),2)),[7]))
514 //longitudinal ridge: [8]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[9]),2)),[10]))
515 //wing structures: [11]*TMath::Power(x,2)
516 //flow contribution (v1 up to v4): 2.*([12]*TMath::Cos(y) + [13]*TMath::Cos(2.*y) + [14]*TMath::Cos(3.*y) + [15]*TMath::Cos(4.*y))
517 TF2 *gFitFunction = new TF2("gFitFunction","[0]+[1]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4]))+[5]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((y-TMath::Pi())/[6]),2)),[7]))+[8]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((x+[17])/[9]),2)),[10]))+[8]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((x-[17])/[9]),2)),[10]))+[11]*TMath::Power(x,2)+2.*[12]*([13]*TMath::Cos(y) + [14]*TMath::Cos(2.*y) + [15]*TMath::Cos(3.*y) + [16]*TMath::Cos(4.*y))",-2.0,2.0,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
518 gFitFunction->SetName("gFitFunction");
519
520
521 //Normalization
522 gFitFunction->SetParName(0,"N1");
523 //near side peak
524 gFitFunction->SetParName(1,"N_{near side}");gFitFunction->FixParameter(1,0.0);
525 gFitFunction->SetParName(2,"Sigma_{near side}(delta eta)"); gFitFunction->FixParameter(2,0.0);
526 gFitFunction->SetParName(3,"Sigma_{near side}(delta phi)"); gFitFunction->FixParameter(3,0.0);
527 gFitFunction->SetParName(4,"Exponent_{near side}"); gFitFunction->FixParameter(4,1.0);
528 //away side ridge
529 gFitFunction->SetParName(5,"N_{away side}"); gFitFunction->FixParameter(5,0.0);
530 gFitFunction->SetParName(6,"Sigma_{away side}(delta phi)"); gFitFunction->FixParameter(6,0.0);
531 gFitFunction->SetParName(7,"Exponent_{away side}"); gFitFunction->FixParameter(7,1.0);
532 //longitudinal ridge
533 gFitFunction->SetParName(8,"N_{long. ridge}"); gFitFunction->SetParameter(8,0.05);//
534 gFitFunction->FixParameter(8,0.0);
535 gFitFunction->SetParName(9,"Sigma_{long. ridge}(delta eta)"); gFitFunction->FixParameter(9,0.0);
536 gFitFunction->SetParName(10,"Exponent_{long. ridge}"); gFitFunction->FixParameter(10,1.0);
537 //wing structures
538 gFitFunction->SetParName(11,"N_{wing}"); gFitFunction->FixParameter(11,0.0);
539 //flow contribution
540 gFitFunction->SetParName(12,"N_{flow}"); gFitFunction->FixParameter(12,0.0);
541 gFitFunction->SetParName(13,"V1"); gFitFunction->FixParameter(13,0.0);
542 gFitFunction->SetParName(14,"V2"); gFitFunction->FixParameter(14,0.0);
543 gFitFunction->SetParName(15,"V3"); gFitFunction->FixParameter(15,0.0);
544 gFitFunction->SetParName(16,"V4"); gFitFunction->FixParameter(16,0.0);
545 gFitFunction->SetParName(17,"Offset"); gFitFunction->FixParameter(17,0.0);
546
547 // just release the near side peak
548 gFitFunction->ReleaseParameter(0);gFitFunction->SetParameter(0,1.0);
549 gFitFunction->ReleaseParameter(1);gFitFunction->SetParameter(1,0.3);gFitFunction->SetParLimits(1,0.0,100);
550 gFitFunction->ReleaseParameter(2);gFitFunction->SetParameter(2,0.3);gFitFunction->SetParLimits(2,0.05,0.7);
551 gFitFunction->ReleaseParameter(3);gFitFunction->SetParameter(3,0.3);gFitFunction->SetParLimits(3,0.05,1.7);
552
553 gHist->Fit("gFitFunction","nm");
554
555
556 //Cloning the histogram
557 TString histoName = gHist->GetName(); histoName += "Fit";
558 TH2D *gHistFit = new TH2D(histoName.Data(),";#Delta#eta;#Delta#varphi (rad);C(#Delta#eta,#Delta#varphi)",gHist->GetNbinsX(),gHist->GetXaxis()->GetXmin(),gHist->GetXaxis()->GetXmax(),gHist->GetNbinsY(),gHist->GetYaxis()->GetXmin(),gHist->GetYaxis()->GetXmax());
559 TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
560 gHistResidual->SetName("gHistResidual");
561 gHistResidual->Sumw2();
562
563 for(Int_t iBinDeltaEta = 1; iBinDeltaEta <= gHist->GetNbinsX(); iBinDeltaEta++) {
564 for(Int_t iBinDeltaPhi = 1; iBinDeltaPhi <= gHist->GetNbinsY(); iBinDeltaPhi++) {
565 gHistFit->SetBinContent(iBinDeltaEta,iBinDeltaPhi,gFitFunction->Eval(gHist->GetXaxis()->GetBinCenter(iBinDeltaEta),gHist->GetYaxis()->GetBinCenter(iBinDeltaPhi)));
566 }
567 }
568 gHistResidual->Add(gHistFit,-1);
569
570 //Write to output file
d67eae53 571 TString newFileName = "balanceFunctionFit2D.";
572 if(eventClass == "Centrality"){
573 newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
574 newFileName += ".PsiAll.PttFrom";
575 }
576 else if(eventClass == "Multiplicity"){
577 newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
578 newFileName += ".PsiAll.PttFrom";
579 }
580 else{ // "EventPlane" (default)
581 newFileName += "Centrality";
582 newFileName += gCentrality; newFileName += ".Psi";
583 if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
584 else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
585 else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
586 else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
587 else newFileName += "All.PttFrom";
588 }
742af4bd 589 newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
590 newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
591 newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
592 newFileName += Form("%.1f",ptAssociatedMax);
d67eae53 593 if(k2pMethod) newFileName += "_2pMethod";
742af4bd 594 newFileName += ".root";
595 TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
596 gHist->Write();
597 gHistFit->Write();
598 gHistResidual->Write();
599 gFitFunction->Write();
600 newFile->Close();
601
602
603}
604
605
606
eb63b883 607//____________________________________________________________//
648f1a5a 608void drawBFPsi2D(const char* lhcPeriod = "LHC11h",
609 Int_t gCentrality = 1,
610 Bool_t kShowShuffled = kFALSE,
611 Bool_t kShowMixed = kFALSE,
eb63b883 612 Double_t psiMin = -0.5, Double_t psiMax = 0.5,
613 Double_t ptTriggerMin = -1.,
614 Double_t ptTriggerMax = -1.,
615 Double_t ptAssociatedMin = -1.,
616 Double_t ptAssociatedMax = -1.) {
617 //Macro that draws the BF distributions for each centrality bin
618 //for reaction plane dependent analysis
619 //Author: Panos.Christakoglou@nikhef.nl
db7174c0 620 TGaxis::SetMaxDigits(3);
eb63b883 621
622 //Get the input file
648f1a5a 623 TString filename = lhcPeriod; filename +="/PttFrom";
624 filename += Form("%.1f",ptTriggerMin); filename += "To";
625 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
626 filename += Form("%.1f",ptAssociatedMin); filename += "To";
627 filename += Form("%.1f",ptAssociatedMax);
628 filename += "/balanceFunction2D.Centrality";
eb63b883 629 filename += gCentrality; filename += ".Psi";
630 if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
631 else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
632 else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
648f1a5a 633 else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
634 else filename += "All.PttFrom";
635 filename += Form("%.1f",ptTriggerMin); filename += "To";
db7174c0 636 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
648f1a5a 637 filename += Form("%.1f",ptAssociatedMin); filename += "To";
638 filename += Form("%.1f",ptAssociatedMax); filename += ".root";
eb63b883 639
640 //Open the file
641 TFile *f = TFile::Open(filename.Data());
642 if((!f)||(!f->IsOpen())) {
643 Printf("The file %s is not found. Aborting...",filename);
644 return listBF;
6acdbcb2 645 }
eb63b883 646 //f->ls();
647
648 //Raw balance function
649 TH1D *gHistBalanceFunction = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunction"));
650 gHistBalanceFunction->SetStats(kFALSE);
651 gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
652 gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
653 gHistBalanceFunction->GetZaxis()->SetNdivisions(10);
654 gHistBalanceFunction->GetXaxis()->SetTitleOffset(1.3);
655 gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
656 gHistBalanceFunction->GetZaxis()->SetTitleOffset(1.3);
657 gHistBalanceFunction->GetXaxis()->SetTitle("#Delta #eta");
658 gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (deg.)");
15dd45a0 659 gHistBalanceFunction->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
eb63b883 660
661 //Shuffled balance function
648f1a5a 662 if(kShowShuffled) {
663 TH1D *gHistBalanceFunctionShuffled = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionShuffled"));
664 gHistBalanceFunctionShuffled->SetStats(kFALSE);
665 gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
666 gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
667 gHistBalanceFunctionShuffled->GetZaxis()->SetNdivisions(10);
668 gHistBalanceFunctionShuffled->GetXaxis()->SetTitleOffset(1.3);
669 gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
670 gHistBalanceFunctionShuffled->GetZaxis()->SetTitleOffset(1.3);
671 gHistBalanceFunctionShuffled->GetXaxis()->SetTitle("#Delta #eta");
672 gHistBalanceFunctionShuffled->GetYaxis()->SetTitle("#Delta #varphi (deg.)");
673 gHistBalanceFunctionShuffled->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
674 }
eb63b883 675
676 //Mixed balance function
648f1a5a 677 if(kShowMixed) {
678 TH1D *gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionMixed"));
679 gHistBalanceFunctionMixed->SetStats(kFALSE);
680 gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
681 gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
682 gHistBalanceFunctionMixed->GetZaxis()->SetNdivisions(10);
683 gHistBalanceFunctionMixed->GetXaxis()->SetTitleOffset(1.3);
684 gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
685 gHistBalanceFunctionMixed->GetZaxis()->SetTitleOffset(1.3);
686 gHistBalanceFunctionMixed->GetXaxis()->SetTitle("#Delta #eta");
687 gHistBalanceFunctionMixed->GetYaxis()->SetTitle("#Delta #varphi (deg.)");
688 gHistBalanceFunctionMixed->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
689 }
eb63b883 690
691 //Subtracted balance function
648f1a5a 692 if(kShowMixed) {
693 TH1D *gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionSubtracted"));
694 gHistBalanceFunctionSubtracted->SetStats(kFALSE);
695 gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
696 gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
697 gHistBalanceFunctionSubtracted->GetZaxis()->SetNdivisions(10);
698 gHistBalanceFunctionSubtracted->GetXaxis()->SetTitleOffset(1.3);
699 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
700 gHistBalanceFunctionSubtracted->GetZaxis()->SetTitleOffset(1.3);
701 gHistBalanceFunctionSubtracted->GetXaxis()->SetTitle("#Delta #eta");
702 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitle("#Delta #varphi (deg.)");
703 gHistBalanceFunctionSubtracted->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
704 }
705
eb63b883 706 TString pngName;
707
708 TString centralityLatex = "Centrality: ";
709 centralityLatex += centralityArray[gCentrality-1];
710 centralityLatex += "%";
711
712 TString psiLatex;
713 if((psiMin == -0.5)&&(psiMax == 0.5))
714 psiLatex = " -7.5^{o} < #varphi - #Psi_{2} < 7.5^{o}";
715 else if((psiMin == 0.5)&&(psiMax == 1.5))
716 psiLatex = " 37.5^{o} < #varphi - #Psi_{2} < 52.5^{o}";
717 else if((psiMin == 1.5)&&(psiMax == 2.5))
718 psiLatex = " 82.5^{o} < #varphi - #Psi_{2} < 97.5^{o}";
719 else
720 psiLatex = " 0^{o} < #varphi - #Psi_{2} < 180^{o}";
721
722 TString pttLatex = Form("%.1f",ptTriggerMin);
723 pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
724 pttLatex += " GeV/c";
725
726 TString ptaLatex = Form("%.1f",ptAssociatedMin);
727 ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
728 ptaLatex += " GeV/c";
729
730 TLatex *latexInfo1 = new TLatex();
731 latexInfo1->SetNDC();
15dd45a0 732 latexInfo1->SetTextSize(0.045);
eb63b883 733 latexInfo1->SetTextColor(1);
734
735 //Draw the results
736 TCanvas *c1 = new TCanvas("c1","Raw balance function 2D",0,0,600,500);
737 c1->SetFillColor(10); c1->SetHighLightColor(10);
738 c1->SetLeftMargin(0.17); c1->SetTopMargin(0.05);
739 gHistBalanceFunction->SetTitle("Raw balance function");
740 gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.4);
741 gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
742 gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
743 gHistBalanceFunction->DrawCopy("lego2");
5de9ad1a 744 gPad->SetTheta(30); // default is 30
745 gPad->SetPhi(-60); // default is 30
746 gPad->Update();
eb63b883 747
15dd45a0 748 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
749 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
750 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
751 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
eb63b883 752
648f1a5a 753 if(kShowShuffled) {
754 TCanvas *c2 = new TCanvas("c2","Shuffled balance function 2D",100,100,600,500);
755 c2->SetFillColor(10); c2->SetHighLightColor(10);
756 c2->SetLeftMargin(0.17); c2->SetTopMargin(0.05);
757 gHistBalanceFunctionShuffled->SetTitle("Shuffled events");
758 gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.4);
759 gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
760 gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
761 gHistBalanceFunctionShuffled->DrawCopy("lego2");
5de9ad1a 762 gPad->SetTheta(30); // default is 30
763 gPad->SetPhi(-60); // default is 30
764 gPad->Update();
765
648f1a5a 766 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
767 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
768 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
769 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
770 }
eb63b883 771
648f1a5a 772 if(kShowMixed) {
773 TCanvas *c3 = new TCanvas("c3","Mixed balance function 2D",200,200,600,500);
774 c3->SetFillColor(10); c3->SetHighLightColor(10);
775 c3->SetLeftMargin(0.17); c3->SetTopMargin(0.05);
776 gHistBalanceFunctionMixed->SetTitle("Mixed events");
777 gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.4);
778 gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
779 gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
780 gHistBalanceFunctionMixed->DrawCopy("lego2");
5de9ad1a 781 gPad->SetTheta(30); // default is 30
782 gPad->SetPhi(-60); // default is 30
783 gPad->Update();
784
648f1a5a 785 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
786 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
787 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
788 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
789 }
eb63b883 790
648f1a5a 791 if(kShowMixed) {
792 TCanvas *c4 = new TCanvas("c4","Subtracted balance function 2D",300,300,600,500);
793 c4->SetFillColor(10); c4->SetHighLightColor(10);
794 c4->SetLeftMargin(0.17); c4->SetTopMargin(0.05);
795 gHistBalanceFunctionSubtracted->SetTitle("Subtracted balance function");
796 gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.4);
797 gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
798 gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
799 gHistBalanceFunctionSubtracted->DrawCopy("lego2");
5de9ad1a 800 gPad->SetTheta(30); // default is 30
801 gPad->SetPhi(-60); // default is 30
802 gPad->Update();
803
648f1a5a 804 latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
805 latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
806 latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
807 latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
808 }
6acdbcb2 809}