]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/macros/drawCorrelationFunctionPsiChargeIndependent.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
new task MK
[u/mrichter/AliRoot.git] / PWGCF / EBYE / macros / drawCorrelationFunctionPsiChargeIndependent.C
CommitLineData
5c43679b 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"};
a67d6cee 3
4const Int_t gRebin = 1;
5
6void drawCorrelationFunctionPsiChargeIndependent(const char* filename = "AnalysisResults.root",
7 Int_t gCentrality = 1,
8 Int_t gBit = -1,
9 const char* gCentralityEstimator = 0x0,
10 Bool_t kShowShuffled = kFALSE,
11 Bool_t kShowMixed = kTRUE,
12 Double_t psiMin = -0.5,
13 Double_t psiMax = 3.5,
34616eda 14 Double_t vertexZMin = -10.,
15 Double_t vertexZMax = 10.,
a67d6cee 16 Double_t ptTriggerMin = -1.,
17 Double_t ptTriggerMax = -1.,
18 Double_t ptAssociatedMin = -1.,
19 Double_t ptAssociatedMax = -1.,
20 Bool_t normToTrig = kFALSE,
21 Int_t rebinEta = 1,
22 Int_t rebinPhi = 1) {
23 //Macro that draws the correlation functions from the balance function
24 //analysis vs the reaction plane
25 //Author: Panos.Christakoglou@nikhef.nl
26 gROOT->LoadMacro("~/SetPlotStyle.C");
27 SetPlotStyle();
28 gStyle->SetPalette(1,0);
29
30 //Load the PWG2ebye library
31 gSystem->Load("libANALYSIS.so");
32 gSystem->Load("libANALYSISalice.so");
33 gSystem->Load("libEventMixing.so");
34 gSystem->Load("libCORRFW.so");
35 gSystem->Load("libPWGTools.so");
36 gSystem->Load("libPWGCFebye.so");
37
38 //Prepare the objects and return them
39 TList *list = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,0);
40 TList *listShuffled = NULL;
41 if(kShowShuffled) listShuffled = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,1);
42 TList *listMixed = NULL;
43 if(kShowMixed) listMixed = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,2);
44
45 if(!list) {
46 Printf("The TList object was not created");
47 return;
48 }
49 else
50 draw(list,listShuffled,listMixed,
34616eda 51 gCentralityEstimator,gCentrality,psiMin,psiMax,vertexZMin,vertexZMax,
a67d6cee 52 ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,normToTrig,rebinEta,rebinPhi);
53}
54
55//______________________________________________________//
56TList *GetListOfObjects(const char* filename,
57 Int_t gCentrality,
58 Int_t gBit,
59 const char *gCentralityEstimator,
60 Int_t kData = 1) {
61 //Get the TList objects (QA, bf, bf shuffled)
62 TList *listBF = 0x0;
63
64 //Open the file
65 TFile *f = TFile::Open(filename,"UPDATE");
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;
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];
93 if(gBit > -1) {
94 listBFName += "_Bit"; listBFName += gBit; }
95 if(gCentralityEstimator) {
96 listBFName += "_"; listBFName += gCentralityEstimator;}
97
98 // histograms were already retrieved (in first iteration)
99 if(dir->Get(Form("%s_histograms",listBFName.Data()))){
100 listBF = dynamic_cast<TList *>(dir->Get(Form("%s_histograms",listBFName.Data())));
101 }
102
103 // histograms were not yet retrieved (this is the first iteration)
104 else{
105
106 listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
107 cout<<"======================================================="<<endl;
108 cout<<"List name: "<<listBF->GetName()<<endl;
109 //listBF->ls();
110
111 //Get the histograms
112 TString histoName;
113 if(kData == 0)
114 histoName = "fHistP";
115 else if(kData == 1)
116 histoName = "fHistP_shuffle";
117 else if(kData == 2)
118 histoName = "fHistP";
119 if(gCentralityEstimator)
120 histoName += gCentralityEstimator;
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)
129 histoName = "fHistN";
130 if(kData == 1)
131 histoName = "fHistN_shuffle";
132 if(kData == 2)
133 histoName = "fHistN";
134 if(gCentralityEstimator)
135 histoName += gCentralityEstimator;
136 AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
137 if(!fHistN) {
138 Printf("fHistN %s not found!!!",histoName.Data());
139 break;
140 }
141 fHistN->FillParent(); fHistN->DeleteContainers();
142
143 if(kData == 0)
144 histoName = "fHistPN";
145 if(kData == 1)
146 histoName = "fHistPN_shuffle";
147 if(kData == 2)
148 histoName = "fHistPN";
149 if(gCentralityEstimator)
150 histoName += gCentralityEstimator;
151 AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
152 if(!fHistPN) {
153 Printf("fHistPN %s not found!!!",histoName.Data());
154 break;
155 }
156 fHistPN->FillParent(); fHistPN->DeleteContainers();
157
158 if(kData == 0)
159 histoName = "fHistNP";
160 if(kData == 1)
161 histoName = "fHistNP_shuffle";
162 if(kData == 2)
163 histoName = "fHistNP";
164 if(gCentralityEstimator)
165 histoName += gCentralityEstimator;
166 AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
167 if(!fHistNP) {
168 Printf("fHistNP %s not found!!!",histoName.Data());
169 break;
170 }
171 fHistNP->FillParent(); fHistNP->DeleteContainers();
172
173 if(kData == 0)
174 histoName = "fHistPP";
175 if(kData == 1)
176 histoName = "fHistPP_shuffle";
177 if(kData == 2)
178 histoName = "fHistPP";
179 if(gCentralityEstimator)
180 histoName += gCentralityEstimator;
181 AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
182 if(!fHistPP) {
183 Printf("fHistPP %s not found!!!",histoName.Data());
184 break;
185 }
186 fHistPP->FillParent(); fHistPP->DeleteContainers();
187
188 if(kData == 0)
189 histoName = "fHistNN";
190 if(kData == 1)
191 histoName = "fHistNN_shuffle";
192 if(kData == 2)
193 histoName = "fHistNN";
194 if(gCentralityEstimator)
195 histoName += gCentralityEstimator;
196 AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
197 if(!fHistNN) {
198 Printf("fHistNN %s not found!!!",histoName.Data());
199 break;
200 }
201 fHistNN->FillParent(); fHistNN->DeleteContainers();
202
203 dir->cd();
204 listBF->Write(Form("%s_histograms",listBFName.Data()), TObject::kSingleKey);
205
206 }// first iteration
207
208 f->Close();
209
210 return listBF;
211}
212
213//______________________________________________________//
214void draw(TList *list, TList *listBFShuffled, TList *listBFMixed,
215 const char *gCentralityEstimator,
216 Int_t gCentrality, Double_t psiMin, Double_t psiMax,
34616eda 217 Double_t vertexZMin,
218 Double_t vertexZMax,
a67d6cee 219 Double_t ptTriggerMin, Double_t ptTriggerMax,
220 Double_t ptAssociatedMin, Double_t ptAssociatedMax,
221 Bool_t normToTrig, Int_t rebinEta, Int_t rebinPhi) {
222 //Draws the correlation functions for every centrality bin
223 //(+-), (-+), (++), (--)
224 AliTHn *hP = NULL;
225 AliTHn *hN = NULL;
226 AliTHn *hPN = NULL;
227 AliTHn *hNP = NULL;
228 AliTHn *hPP = NULL;
229 AliTHn *hNN = NULL;
230
231 TString gHistPname = "fHistP";
232 if(gCentralityEstimator) gHistPname += gCentralityEstimator;
233 TString gHistNname = "fHistN";
234 if(gCentralityEstimator) gHistNname += gCentralityEstimator;
235 TString gHistPNname = "fHistPN";
236 if(gCentralityEstimator) gHistPNname += gCentralityEstimator;
237 TString gHistNPname = "fHistNP";
238 if(gCentralityEstimator) gHistNPname += gCentralityEstimator;
239 TString gHistPPname = "fHistPP";
240 if(gCentralityEstimator) gHistPPname += gCentralityEstimator;
241 TString gHistNNname = "fHistNN";
242 if(gCentralityEstimator) gHistNNname += gCentralityEstimator;
243
244 hP = (AliTHn*) list->FindObject(gHistPname.Data());
245 hN = (AliTHn*) list->FindObject(gHistNname.Data());
246 hPN = (AliTHn*) list->FindObject(gHistPNname.Data());
247 hNP = (AliTHn*) list->FindObject(gHistNPname.Data());
248 hPP = (AliTHn*) list->FindObject(gHistPPname.Data());
249 hNN = (AliTHn*) list->FindObject(gHistNNname.Data());
250
251 //Create the AliBalancePsi object and fill it with the AliTHn objects
252 AliBalancePsi *b = new AliBalancePsi();
253 b->SetHistNp(hP);
254 b->SetHistNn(hN);
255 b->SetHistNpn(hPN);
256 b->SetHistNnp(hNP);
257 b->SetHistNpp(hPP);
258 b->SetHistNnn(hNN);
259
260 //balance function shuffling
261 AliTHn *hPShuffled = NULL;
262 AliTHn *hNShuffled = NULL;
263 AliTHn *hPNShuffled = NULL;
264 AliTHn *hNPShuffled = NULL;
265 AliTHn *hPPShuffled = NULL;
266 AliTHn *hNNShuffled = NULL;
267 if(listBFShuffled) {
268 //listBFShuffled->ls();
269
270 gHistPname = "fHistP_shuffle";
271 if(gCentralityEstimator) gHistPname += gCentralityEstimator;
272 gHistNname = "fHistN_shuffle";
273 if(gCentralityEstimator) gHistNname += gCentralityEstimator;
274 gHistPNname = "fHistPN_shuffle";
275 if(gCentralityEstimator) gHistPNname += gCentralityEstimator;
276 gHistNPname = "fHistNP_shuffle";
277 if(gCentralityEstimator) gHistNPname += gCentralityEstimator;
278 gHistPPname = "fHistPP_shuffle";
279 if(gCentralityEstimator) gHistPPname += gCentralityEstimator;
280 gHistNNname = "fHistNN_shuffle";
281 if(gCentralityEstimator) gHistNNname += gCentralityEstimator;
282
283 hPShuffled = (AliTHn*) listBFShuffled->FindObject(gHistPname.Data());
284 hPShuffled->SetName("gHistPShuffled");
285 hNShuffled = (AliTHn*) listBFShuffled->FindObject(gHistNname.Data());
286 hNShuffled->SetName("gHistNShuffled");
287 hPNShuffled = (AliTHn*) listBFShuffled->FindObject(gHistPNname.Data());
288 hPNShuffled->SetName("gHistPNShuffled");
289 hNPShuffled = (AliTHn*) listBFShuffled->FindObject(gHistNPname.Data());
290 hNPShuffled->SetName("gHistNPShuffled");
291 hPPShuffled = (AliTHn*) listBFShuffled->FindObject(gHistPPname.Data());
292 hPPShuffled->SetName("gHistPPShuffled");
293 hNNShuffled = (AliTHn*) listBFShuffled->FindObject(gHistNNname.Data());
294 hNNShuffled->SetName("gHistNNShuffled");
295
296 AliBalancePsi *bShuffled = new AliBalancePsi();
297 bShuffled->SetHistNp(hPShuffled);
298 bShuffled->SetHistNn(hNShuffled);
299 bShuffled->SetHistNpn(hPNShuffled);
300 bShuffled->SetHistNnp(hNPShuffled);
301 bShuffled->SetHistNpp(hPPShuffled);
302 bShuffled->SetHistNnn(hNNShuffled);
303 }
304
305 //balance function mixing
306 AliTHn *hPMixed = NULL;
307 AliTHn *hNMixed = NULL;
308 AliTHn *hPNMixed = NULL;
309 AliTHn *hNPMixed = NULL;
310 AliTHn *hPPMixed = NULL;
311 AliTHn *hNNMixed = NULL;
312
313 if(listBFMixed) {
314 //listBFMixed->ls();
315
316 gHistPname = "fHistP";
317 if(gCentralityEstimator) gHistPname += gCentralityEstimator;
318 gHistNname = "fHistN";
319 if(gCentralityEstimator) gHistNname += gCentralityEstimator;
320 gHistPNname = "fHistPN";
321 if(gCentralityEstimator) gHistPNname += gCentralityEstimator;
322 gHistNPname = "fHistNP";
323 if(gCentralityEstimator) gHistNPname += gCentralityEstimator;
324 gHistPPname = "fHistPP";
325 if(gCentralityEstimator) gHistPPname += gCentralityEstimator;
326 gHistNNname = "fHistNN";
327 if(gCentralityEstimator) gHistNNname += gCentralityEstimator;
328 hPMixed = (AliTHn*) listBFMixed->FindObject(gHistPname.Data());
329 hPMixed->SetName("gHistPMixed");
330 hNMixed = (AliTHn*) listBFMixed->FindObject(gHistNname.Data());
331 hNMixed->SetName("gHistNMixed");
332 hPNMixed = (AliTHn*) listBFMixed->FindObject(gHistPNname.Data());
333 hPNMixed->SetName("gHistPNMixed");
334 hNPMixed = (AliTHn*) listBFMixed->FindObject(gHistNPname.Data());
335 hNPMixed->SetName("gHistNPMixed");
336 hPPMixed = (AliTHn*) listBFMixed->FindObject(gHistPPname.Data());
337 hPPMixed->SetName("gHistPPMixed");
338 hNNMixed = (AliTHn*) listBFMixed->FindObject(gHistNNname.Data());
339 hNNMixed->SetName("gHistNNMixed");
340
341 AliBalancePsi *bMixed = new AliBalancePsi();
342 bMixed->SetHistNp(hPMixed);
343 bMixed->SetHistNn(hNMixed);
344 bMixed->SetHistNpn(hPNMixed);
345 bMixed->SetHistNnp(hNPMixed);
346 bMixed->SetHistNpp(hPPMixed);
347 bMixed->SetHistNnn(hNNMixed);
348 }
349
5c43679b 350 TH2D *gHist[6];
a67d6cee 351
5c43679b 352 TCanvas *c[6];
a67d6cee 353 TString histoTitle, pngName;
354
355 // all charges together
356 histoTitle = "(+-) | Centrality: ";
357 histoTitle += centralityArray[gCentrality-1];
358 histoTitle += "%";
359 if((psiMin == -0.5)&&(psiMax == 0.5))
360 histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
361 else if((psiMin == 0.5)&&(psiMax == 1.5))
362 histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
363 else if((psiMin == 1.5)&&(psiMax == 2.5))
364 histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
365 else
366 histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
367
34616eda 368 gHist[0] = b->GetCorrelationFunctionChargeIndependent(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
faa03677 369 if(rebinEta > 1 || rebinPhi > 1){
370 gHist[0]->Rebin2D(rebinEta,rebinPhi);
371 gHist[0]->Scale(1./(Double_t)(rebinEta*rebinPhi));
372 }
a67d6cee 373 gHist[0]->GetYaxis()->SetTitleOffset(1.5);
374 gHist[0]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
375 gHist[0]->SetTitle(histoTitle.Data());
376 c[0] = new TCanvas("c0","",0,0,600,500);
377 c[0]->SetFillColor(10); c[0]->SetHighLightColor(10);
378 gHist[0]->DrawCopy("surf1fb");
379 gPad->SetTheta(30); // default is 30
380 //gPad->SetPhi(130); // default is 30
381 gPad->SetPhi(-60); // default is 30
382 gPad->Update();
383 pngName = "DeltaPhiDeltaEta.Centrality";
384 pngName += centralityArray[gCentrality-1];
385 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
386 pngName += ".PositiveNegative.png";
387 //c[0]->SaveAs(pngName.Data());
388
389 if(listBFShuffled) {
390 histoTitle = "(+-) shuffled | Centrality: ";
391 histoTitle += centralityArray[gCentrality-1];
392 histoTitle += "%";
393 if((psiMin == -0.5)&&(psiMax == 0.5))
394 histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
395 else if((psiMin == 0.5)&&(psiMax == 1.5))
396 histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
397 else if((psiMin == 1.5)&&(psiMax == 2.5))
398 histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
399 else
400 histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
401
34616eda 402 gHist[1] = bShuffled->GetCorrelationFunctionChargeIndependent(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
faa03677 403 if(rebinEta > 1 || rebinPhi > 1){
404 gHist[1]->Rebin2D(rebinEta,rebinPhi);
405 gHist[1]->Scale(1./(Double_t)(rebinEta*rebinPhi));
406 }
a67d6cee 407 gHist[1]->GetYaxis()->SetTitleOffset(1.5);
408 gHist[1]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
409 gHist[1]->SetTitle(histoTitle.Data());
410 c[1] = new TCanvas("c1","",0,100,600,500);
411 c[1]->SetFillColor(10);
412 c[1]->SetHighLightColor(10);
413 gHist[1]->DrawCopy("surf1fb");
414 gPad->SetTheta(30); // default is 30
415 //gPad->SetPhi(130); // default is 30
416 gPad->SetPhi(-60); // default is 30
417 gPad->Update();
418 pngName = "DeltaPhiDeltaEtaShuffled.Centrality";
419 pngName += centralityArray[gCentrality-1];
420 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
421 pngName += ".PositiveNegative.png";
422 //c[1]->SaveAs(pngName.Data());
423 }
424
425 if(listBFMixed) {
426 histoTitle = "(+-) mixed | Centrality: ";
427 histoTitle += centralityArray[gCentrality-1];
428 histoTitle += "%";
429 if((psiMin == -0.5)&&(psiMax == 0.5))
430 histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
431 else if((psiMin == 0.5)&&(psiMax == 1.5))
432 histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
433 else if((psiMin == 1.5)&&(psiMax == 2.5))
434 histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
435 else
436 histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
437
438 // if normalization to trigger then do not divide Event mixing by number of trigger particles
34616eda 439 gHist[2] = bMixed->GetCorrelationFunctionChargeIndependent(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
faa03677 440 if(rebinEta > 1 || rebinPhi > 1){
441 gHist[2]->Rebin2D(rebinEta,rebinPhi);
442 gHist[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));
443 }
a67d6cee 444 // normalization to 1 at (0,0) --> Jan Fietes method
445 if(normToTrig){
446 Double_t mixedNorm = gHist[2]->Integral(gHist[2]->GetXaxis()->FindBin(0-10e-5),gHist[2]->GetXaxis()->FindBin(0+10e-5),1,gHist[2]->GetNbinsX());
447 mixedNorm /= gHist[2]->GetNbinsY()*(gHist[2]->GetXaxis()->FindBin(0.01) - gHist[2]->GetXaxis()->FindBin(-0.01) + 1);
448 gHist[2]->Scale(1./mixedNorm);
449 }
450
451 gHist[2]->GetYaxis()->SetTitleOffset(1.5);
452 gHist[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
453 gHist[2]->SetTitle(histoTitle.Data());
454 c[2] = new TCanvas("c2","",0,200,600,500);
455 c[2]->SetFillColor(10);
456 c[2]->SetHighLightColor(10);
457 gHist[2]->DrawCopy("surf1fb");
458 gPad->SetTheta(30); // default is 30
459 //gPad->SetPhi(130); // default is 30
460 gPad->SetPhi(-60); // default is 30
461 gPad->Update();
462 pngName = "DeltaPhiDeltaEtaMixed.Centrality";
463 pngName += centralityArray[gCentrality-1];
464 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
465 pngName += ".PositiveNegative.png";
466 //c[2]->SaveAs(pngName.Data());
467
468 //Correlation function (+-)
34616eda 469 gHist[3] = b->GetCorrelationFunction("ALL",psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
a67d6cee 470 gHist[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
471 gHist[3]->GetZaxis()->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
472 c[3] = new TCanvas("c3","",0,300,600,500);
473 c[3]->SetFillColor(10);
474 c[3]->SetHighLightColor(10);
475 gHist[3]->DrawCopy("surf1fb");
476 gPad->SetTheta(30); // default is 30
477 //gPad->SetPhi(130); // default is 30
478 gPad->SetPhi(-60); // default is 30
479 gPad->Update();
480 pngName = "CorrelationFunction.Centrality";
481 pngName += centralityArray[gCentrality-1];
482 pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
483 pngName += ".PositiveNegative.png";
484 //c[3]->SaveAs(pngName.Data());
5c43679b 485
486 //Correlation function subtracted awayside (+-)
487 gHist[4] = dynamic_cast<TH2D *>(gHist[3]->Clone()); // this will be used to imitate twice the away-side
488 gHist[4]->GetXaxis()->SetRangeUser(-1.5,1.5);
489 gHist[4]->GetZaxis()->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
490 gHist[5] = dynamic_cast<TH2D *>(gHist[3]->Clone()); // this will be the subtracted one
491 gHist[5]->GetXaxis()->SetRangeUser(-1.5,1.5);
492 gHist[5]->GetZaxis()->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
493
494 //prepare the double away side histo
495 for(Int_t ix = 0; ix < gHist[4]->GetNbinsX(); ix++ ){
496 for(Int_t iy = 0; iy < gHist[4]->GetNbinsX(); iy++ ){
497 if(iy<gHist[4]->GetNbinsY()/2) gHist[4]->SetBinContent(ix+1,iy+1,gHist[4]->GetBinContent(ix+1,iy+1+gHist[4]->GetNbinsY()/2));
498 }
499 }
500 gHist[5]->Add(gHist[4],-1);
501
502 c[4] = new TCanvas("c3","",0,300,600,500);
503 c[4]->SetFillColor(10);
504 c[4]->SetHighLightColor(10);
505 gHist[5]->DrawCopy("surf1fb");
506 gPad->SetTheta(30); // default is 30
507 //gPad->SetPhi(130); // default is 30
508 gPad->SetPhi(-60); // default is 30
509 gPad->Update();
a67d6cee 510 }
511
512 //Write to output file
513 TString newFileName = "correlationFunctionChargeIndependent.Centrality";
514 newFileName += gCentrality; newFileName += ".Psi";
515 if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
516 else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
517 else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
518 else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
519 else newFileName += "All.PttFrom";
520 newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
521 newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
522 newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
523 newFileName += Form("%.1f",ptAssociatedMax);
524 newFileName += ".root";
525 TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
526 gHist[0]->SetName("gHistRaw"); gHist[0]->Write();
527 if(listBFShuffled) {
528 gHist[1]->SetName("gHistShuffled"); gHist[1]->Write();
529 }
530 if(listBFMixed) {
531 gHist[2]->SetName("gHistMixed"); gHist[2]->Write();
532
5c43679b 533 gHist[3]->SetName("gHistCorrelationFunctions");gHist[3]->Write();
534 gHist[4]->SetName("gHistCorrelationFunctionsAwaySide"); gHist[4]->Write();
535 gHist[5]->SetName("gHistCorrelationFunctionsSubtracted"); gHist[5]->Write();
a67d6cee 536 }
537 newFile->Close();
538
5c43679b 539 // // some cleaning
540 // for(Int_t i = 0; i < 6; i++){
a67d6cee 541
5c43679b 542 // if(!listBFShuffled && i == 1) continue;
543 // if(!listBFMixed && (i == 2 || i == 3 || i == 4 || i == 5)) continue;
a67d6cee 544
5c43679b 545 // if(gHist[i]) delete gHist[i];
a67d6cee 546
5c43679b 547 // if(c[i]) delete c[i];
548 // }
a67d6cee 549
550 delete hP;
551 delete hN;
552 delete hPP;
553 delete hPN;
554 delete hNP;
555 delete hNN;
556
557 delete hPMixed;
558 delete hNMixed;
559 delete hPPMixed;
560 delete hPNMixed;
561 delete hNPMixed;
562 delete hNNMixed;
563
564 delete hPShuffled;
565 delete hNShuffled;
566 delete hPPShuffled;
567 delete hPNShuffled;
568 delete hNPShuffled;
569 delete hNNShuffled;
570
571}
572
573//____________________________________________________________//
574void drawCorrelationFunctions(const char* lhcPeriod = "LHC11h",
575 Int_t gTrainID = 208,
576 Int_t gCentrality = 1,
577 Double_t psiMin = -0.5, Double_t psiMax = 3.5,
578 Double_t ptTriggerMin = -1.,
579 Double_t ptTriggerMax = -1.,
580 Double_t ptAssociatedMin = -1.,
581 Double_t ptAssociatedMax = -1.) {
582 //Macro that draws the charge dependent correlation functions
583 //for each centrality bin for the different pT of trigger and
584 //associated particles
585 //Author: Panos.Christakoglou@nikhef.nl
586 TGaxis::SetMaxDigits(3);
587
588 //Get the input file
589 TString filename = "PbPb/"; filename += lhcPeriod;
590 filename +="/Train"; filename += gTrainID;
591 filename +="/Centrality"; filename += gCentrality;
592 filename += "/correlationFunction.Centrality";
593 filename += gCentrality; filename += ".Psi";
594 if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
595 else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
596 else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
597 else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
598 else filename += "All.PttFrom";
599 filename += Form("%.1f",ptTriggerMin); filename += "To";
600 filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
601 filename += Form("%.1f",ptAssociatedMin); filename += "To";
602 filename += Form("%.1f",ptAssociatedMax);
603 filename += ".root";
604
605 //Open the file
606 TFile *f = TFile::Open(filename.Data());
607 if((!f)||(!f->IsOpen())) {
608 Printf("The file %s is not found. Aborting...",filename);
609 return listBF;
610 }
611 //f->ls();
612
613 //Latex
614 TString centralityLatex = "Centrality: ";
615 centralityLatex += centralityArray[gCentrality-1];
616 centralityLatex += "%";
617
618 TString psiLatex;
619 if((psiMin == -0.5)&&(psiMax == 0.5))
620 psiLatex = " -7.5^{o} < #varphi - #Psi_{2} < 7.5^{o}";
621 else if((psiMin == 0.5)&&(psiMax == 1.5))
622 psiLatex = " 37.5^{o} < #varphi - #Psi_{2} < 52.5^{o}";
623 else if((psiMin == 1.5)&&(psiMax == 2.5))
624 psiLatex = " 82.5^{o} < #varphi - #Psi_{2} < 97.5^{o}";
625 else
626 psiLatex = " 0^{o} < #varphi - #Psi_{2} < 180^{o}";
627
628 TString pttLatex = Form("%.1f",ptTriggerMin);
629 pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
630 pttLatex += " GeV/c";
631
632 TString ptaLatex = Form("%.1f",ptAssociatedMin);
633 ptaLatex += " < p_{T,assoc} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
634 ptaLatex += " GeV/c";
635
636 TLatex *latexInfo1 = new TLatex();
637 latexInfo1->SetNDC();
638 latexInfo1->SetTextSize(0.045);
639 latexInfo1->SetTextColor(1);
640
641 TString pngName;
642
643 //============================================================//
644 //Get the +- correlation function
645 TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
646 gHistPN->SetStats(kFALSE);
647 gHistPN->SetTitle("");
648 gHistPN->GetXaxis()->SetRangeUser(-1.45,1.45);
649 gHistPN->GetXaxis()->CenterTitle();
650 gHistPN->GetXaxis()->SetTitleOffset(1.2);
651 gHistPN->GetYaxis()->CenterTitle();
652 gHistPN->GetYaxis()->SetTitleOffset(1.2);
653 gHistPN->GetZaxis()->SetTitleOffset(1.2);
654 TCanvas *cPN = new TCanvas("cPN","",0,0,600,500);
655 cPN->SetFillColor(10); cPN->SetHighLightColor(10);
656 cPN->SetLeftMargin(0.15);
657 gHistPN->DrawCopy("surf1fb");
658 gPad->SetTheta(30); // default is 30
659 gPad->SetPhi(-60); // default is 30
660 gPad->Update();
661
662 latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
663 //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
664 latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
665 latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
666
667 pngName = "CorrelationFunctionChargeIndependent.Centrality";
668 pngName += centralityArray[gCentrality-1];
669 pngName += ".Psi";
670 if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
671 else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
672 else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
673 else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
674 else pngName += "All.PttFrom";
675 pngName += Form("%.1f",ptTriggerMin); pngName += "To";
676 pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
677 pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
678 pngName += Form("%.1f",ptAssociatedMax);
679 pngName += ".PositiveNegative.png";
680 cPN->SaveAs(pngName.Data());
681 fitCorrelationFunctions(gCentrality, psiMin, psiMax,
682 ptTriggerMin,ptTriggerMax,
683 ptAssociatedMin, ptAssociatedMax,gHistPN);
684 //============================================================//
685 //Get the -+ correlation function
686 TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
687 gHistNP->SetStats(kFALSE);
688 gHistNP->SetTitle("");
689 gHistNP->GetXaxis()->SetRangeUser(-1.45,1.45);
690 gHistNP->GetXaxis()->CenterTitle();
691 gHistNP->GetXaxis()->SetTitleOffset(1.2);
692 gHistNP->GetYaxis()->CenterTitle();
693 gHistNP->GetYaxis()->SetTitleOffset(1.2);
694 gHistNP->GetZaxis()->SetTitleOffset(1.2);
695 TCanvas *cNP = new TCanvas("cNP","",50,50,600,500);
696 cNP->SetFillColor(10); cNP->SetHighLightColor(10);
697 cNP->SetLeftMargin(0.15);
698 gHistNP->DrawCopy("surf1fb");
699 gPad->SetTheta(30); // default is 30
700 gPad->SetPhi(-60); // default is 30
701 gPad->Update();
702
703 latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
704 //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
705 latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
706 latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
707
708 pngName = "CorrelationFunction.Centrality";
709 pngName += centralityArray[gCentrality-1];
710 pngName += ".Psi";
711 if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
712 else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
713 else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
714 else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
715 else pngName += "All.PttFrom";
716 pngName += Form("%.1f",ptTriggerMin); pngName += "To";
717 pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
718 pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
719 pngName += Form("%.1f",ptAssociatedMax);
720 pngName += ".NegativePositive.png";
721 cNP->SaveAs(pngName.Data());
722
723 fitCorrelationFunctions(gCentrality, psiMin, psiMax,
724 ptTriggerMin,ptTriggerMax,
725 ptAssociatedMin, ptAssociatedMax,gHistNP);
726 //============================================================//
727 //Get the ++ correlation function
728 TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
729 gHistPP->SetStats(kFALSE);
730 gHistPP->SetTitle("");
731 gHistPP->GetXaxis()->SetRangeUser(-1.45,1.45);
732 gHistPP->GetXaxis()->CenterTitle();
733 gHistPP->GetXaxis()->SetTitleOffset(1.2);
734 gHistPP->GetYaxis()->CenterTitle();
735 gHistPP->GetYaxis()->SetTitleOffset(1.2);
736 gHistPP->GetZaxis()->SetTitleOffset(1.2);
737 TCanvas *cPP = new TCanvas("cPP","",100,100,600,500);
738 cPP->SetFillColor(10); cPP->SetHighLightColor(10);
739 cPP->SetLeftMargin(0.15);
740 gHistPP->DrawCopy("surf1fb");
741 gPad->SetTheta(30); // default is 30
742 gPad->SetPhi(-60); // default is 30
743 gPad->Update();
744
745 latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
746 //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
747 latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
748 latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
749
750 pngName = "CorrelationFunction.Centrality";
751 pngName += centralityArray[gCentrality-1];
752 pngName += ".Psi";
753 if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
754 else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
755 else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
756 else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
757 else pngName += "All.PttFrom";
758 pngName += Form("%.1f",ptTriggerMin); pngName += "To";
759 pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
760 pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
761 pngName += Form("%.1f",ptAssociatedMax);
762 pngName += ".PositivePositive.png";
763 cPP->SaveAs(pngName.Data());
764
765 fitCorrelationFunctions(gCentrality, psiMin, psiMax,
766 ptTriggerMin,ptTriggerMax,
767 ptAssociatedMin, ptAssociatedMax,gHistPP);
768 //============================================================//
769 //Get the -- correlation function
770 TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
771 gHistNN->SetStats(kFALSE);
772 gHistNN->SetTitle("");
773 gHistNN->GetXaxis()->SetRangeUser(-1.45,1.45);
774 gHistNN->GetXaxis()->CenterTitle();
775 gHistNN->GetXaxis()->SetTitleOffset(1.2);
776 gHistNN->GetYaxis()->CenterTitle();
777 gHistNN->GetYaxis()->SetTitleOffset(1.2);
778 gHistNN->GetZaxis()->SetTitleOffset(1.2);
779 TCanvas *cNN = new TCanvas("cNN","",150,150,600,500);
780 cNN->SetFillColor(10); cNN->SetHighLightColor(10);
781 cNN->SetLeftMargin(0.15);
782 gHistNN->DrawCopy("surf1fb");
783 gPad->SetTheta(30); // default is 30
784 gPad->SetPhi(-60); // default is 30
785 gPad->Update();
786
787 latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
788 //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
789 latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
790 latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
791
792 pngName = "CorrelationFunction.Centrality";
793 pngName += centralityArray[gCentrality-1];
794 pngName += ".Psi";
795 if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
796 else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
797 else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
798 else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
799 else pngName += "All.PttFrom";
800 pngName += Form("%.1f",ptTriggerMin); pngName += "To";
801 pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
802 pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
803 pngName += Form("%.1f",ptAssociatedMax);
804 pngName += ".NegativeNegative.png";
805 cNN->SaveAs(pngName.Data());
806
807 fitCorrelationFunctions(gCentrality, psiMin, psiMax,
808 ptTriggerMin,ptTriggerMax,
809 ptAssociatedMin, ptAssociatedMax,gHistNN);
810}
811
812// //____________________________________________________________//
813// void fitCorrelationFunctions(Int_t gCentrality = 1,
814// Double_t psiMin = -0.5, Double_t psiMax = 3.5,
815// Double_t ptTriggerMin = -1.,
816// Double_t ptTriggerMax = -1.,
817// Double_t ptAssociatedMin = -1.,
818// Double_t ptAssociatedMax = -1.,
819// TH2D *gHist) {
820
821// cout<<"FITTING FUNCTION (MW style)"<<endl;
822
823// //near side peak(s): [1]*(TMath::Exp(-TMath::Power((0.5*TMath::Power(((x-[5])/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4])) +
824// // TMath::Exp(-TMath::Power((0.5*TMath::Power(((x+[5])/[6]),2)+0.5*TMath::Power((y/[3]),2)),[4])))
825// //away side peak(s): [7]*(TMath::Exp(-TMath::Power((0.5*TMath::Power(((x-[11])/[8]),2)+0.5*TMath::Power((y/[9]),2)),[10])) +
826// // TMath::Exp(-TMath::Power((0.5*TMath::Power(((x+[11])/[12]),2)+0.5*TMath::Power((y/[9]),2)),[10])))
827// //flow contribution (v1 up to v4): 2.*[13]*([14]*TMath::Cos(y) + [15]*TMath::Cos(2.*y) + [16]*TMath::Cos(3.*y) + [17]*TMath::Cos(4.*y))
828
829
830
831// TF2 *gFitFunction = new TF2("gFitFunction","[0]+[1]*(TMath::Exp(-TMath::Power((0.5*TMath::Power(((x-[5])/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4])) + TMath::Exp(-TMath::Power((0.5*TMath::Power(((x+[5])/[6]),2)+0.5*TMath::Power((y/[3]),2)),[4]))) + [7]*(TMath::Exp(-TMath::Power((0.5*TMath::Power(((x-[11])/[8]),2)+0.5*TMath::Power(((y-TMath::Pi())/[9]),2)),[10])) + TMath::Exp(-TMath::Power((0.5*TMath::Power(((x+[11])/[12]),2)+0.5*TMath::Power(((y-TMath::Pi())/[9]),2)),[10]))) + 2.*[13]*([14]*TMath::Cos(y) + [15]*TMath::Cos(2.*y) + [16]*TMath::Cos(3.*y) + [17]*TMath::Cos(4.*y))",-2.0,2.0,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
832// gFitFunction->SetName("gFitFunction");
833
834
835// //Normalization
836// gFitFunction->SetParName(0,"N1");
837// //near side peak(s)
838// gFitFunction->SetParName(1,"N_{near side}");gFitFunction->FixParameter(1,0.0);
839// gFitFunction->SetParName(2,"Sigma_{near side}(delta eta 1)"); gFitFunction->FixParameter(2,0.0);
840// gFitFunction->SetParName(3,"Sigma_{near side}(delta phi)"); gFitFunction->FixParameter(3,0.0);
841// gFitFunction->SetParName(4,"Exponent_{near side}"); gFitFunction->FixParameter(4,1.0);
842// gFitFunction->SetParName(5,"Offset_{near side}"); gFitFunction->FixParameter(5,0.0);
843// gFitFunction->SetParName(6,"Sigma_{near side}(delta eta 2)"); gFitFunction->FixParameter(6,0.0);
844
845// //away side peak(s)
846// gFitFunction->SetParName(7,"N_{near side}");gFitFunction->FixParameter(7,0.0);
847// gFitFunction->SetParName(8,"Sigma_{near side}(delta eta 1)"); gFitFunction->FixParameter(8,0.0);
848// gFitFunction->SetParName(9,"Sigma_{near side}(delta phi)"); gFitFunction->FixParameter(9,0.0);
849// gFitFunction->SetParName(10,"Exponent_{near side}"); gFitFunction->FixParameter(10,1.0);
850// gFitFunction->SetParName(11,"Offset_{near side}"); gFitFunction->FixParameter(11,0.0);
851// gFitFunction->SetParName(12,"Sigma_{near side}(delta eta 2)"); gFitFunction->FixParameter(12,0.0);
852
853// //flow contribution
854// gFitFunction->SetParName(13,"N_{flow}"); gFitFunction->SetParameter(13,0.2);
855// gFitFunction->SetParName(14,"V1"); gFitFunction->SetParameter(14,0.005);
856// gFitFunction->SetParName(15,"V2"); gFitFunction->SetParameter(15,0.1);
857// gFitFunction->SetParName(16,"V3"); gFitFunction->SetParameter(16,0.05);
858// gFitFunction->SetParName(17,"V4"); gFitFunction->SetParameter(17,0.005);
859
860// // flow parameters
861// Double_t fNV = 0.;
862// Double_t fV1 = 0.;
863// Double_t fV2 = 0.;
864// Double_t fV3 = 0.;
865// Double_t fV4 = 0.;
866
867// //Fitting the correlation function (first the edges to extract flow)
868// gHist->Fit("gFitFunction","nm","",1.0,1.6);
869
870// fNV += gFitFunction->GetParameter(13);
871// fV1 += gFitFunction->GetParameter(14);
872// fV2 += gFitFunction->GetParameter(15);
873// fV3 += gFitFunction->GetParameter(16);
874// fV4 += gFitFunction->GetParameter(17);
875
876// gHist->Fit("gFitFunction","nm","",-1.6,-1.0);
877
878// fNV += gFitFunction->GetParameter(13);
879// fV1 += gFitFunction->GetParameter(14);
880// fV2 += gFitFunction->GetParameter(15);
881// fV3 += gFitFunction->GetParameter(16);
882// fV4 += gFitFunction->GetParameter(17);
883
884// // Now fit the whole with fixed flow
885// gFitFunction->FixParameter(13,fNV/2.);
886// gFitFunction->FixParameter(14,fV1/2.);
887// gFitFunction->FixParameter(15,fV2/2.);
888// gFitFunction->FixParameter(16,fV3/2.);
889// gFitFunction->FixParameter(17,fV4/2.);
890
891// gFitFunction->ReleaseParameter(1);gFitFunction->SetParameter(1,0.3);
892// gFitFunction->ReleaseParameter(2);gFitFunction->SetParameter(2,0.3);gFitFunction->SetParLimits(2,0.05,0.7);
893// gFitFunction->ReleaseParameter(3);gFitFunction->SetParameter(3,0.3);gFitFunction->SetParLimits(3,0.05,1.7);
894// gFitFunction->ReleaseParameter(5);gFitFunction->SetParameter(5,0.7);gFitFunction->SetParLimits(5,0.0,0.9);
895// gFitFunction->ReleaseParameter(6);gFitFunction->SetParameter(6,0.3);gFitFunction->SetParLimits(6,0.01,1.7);
896
897// gFitFunction->ReleaseParameter(7);gFitFunction->SetParameter(1,0.3);
898// gFitFunction->ReleaseParameter(8);gFitFunction->SetParameter(2,0.3);gFitFunction->SetParLimits(2,0.05,0.7);
899// gFitFunction->ReleaseParameter(9);gFitFunction->SetParameter(3,0.3);gFitFunction->SetParLimits(3,0.05,1.7);
900// gFitFunction->ReleaseParameter(11);gFitFunction->SetParameter(5,0.7);gFitFunction->SetParLimits(5,0.0,0.9);
901// gFitFunction->ReleaseParameter(12);gFitFunction->SetParameter(6,0.3);gFitFunction->SetParLimits(6,0.01,1.7);
902
903// gHist->Fit("gFitFunction","nm");
904
905
906// //Cloning the histogram
907// TString histoName = gHist->GetName(); histoName += "Fit";
908// 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());
909// TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
910// gHistResidual->SetName("gHistResidual");
911// gHistResidual->Sumw2();
912
913// for(Int_t iBinDeltaEta = 1; iBinDeltaEta <= gHist->GetNbinsX(); iBinDeltaEta++) {
914// for(Int_t iBinDeltaPhi = 1; iBinDeltaPhi <= gHist->GetNbinsY(); iBinDeltaPhi++) {
915// gHistFit->SetBinContent(iBinDeltaEta,iBinDeltaPhi,gFitFunction->Eval(gHist->GetXaxis()->GetBinCenter(iBinDeltaEta),gHist->GetYaxis()->GetBinCenter(iBinDeltaPhi)));
916// }
917// }
918// gHistResidual->Add(gHistFit,-1);
919
920// //Write to output file
921// TString newFileName = "correlationFunctionFit";
922// if(histoName.Contains("PN")) newFileName += "PN";
923// else if(histoName.Contains("NP")) newFileName += "NP";
924// else if(histoName.Contains("PP")) newFileName += "PP";
925// else if(histoName.Contains("NN")) newFileName += "NN";
926// newFileName += ".Centrality";
927// newFileName += gCentrality; newFileName += ".Psi";
928// if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
929// else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
930// else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
931// else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
932// else newFileName += "All.PttFrom";
933// newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
934// newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
935// newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
936// newFileName += Form("%.1f",ptAssociatedMax);
937// newFileName += ".root";
938// TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
939// gHist->Write();
940// gHistFit->Write();
941// gHistResidual->Write();
942// gFitFunction->Write();
943// newFile->Close();
944
945
946// }
947
948//____________________________________________________________//
949void fitCorrelationFunctions(Int_t gCentrality = 1,
950 Double_t psiMin = -0.5, Double_t psiMax = 3.5,
951 Double_t ptTriggerMin = -1.,
952 Double_t ptTriggerMax = -1.,
953 Double_t ptAssociatedMin = -1.,
954 Double_t ptAssociatedMax = -1.,
955 TH2D *gHist) {
956
957 cout<<"FITTING FUNCTION"<<endl;
958
959 //near side peak: [1]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4]))
960 //away side ridge: [5]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((y-TMath::Pi())/[6]),2)),[7]))
961 //longitudinal ridge: [8]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[9]),2)),[10]))
962 //wing structures: [11]*TMath::Power(x,2)
963 //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))
964 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.);
965 gFitFunction->SetName("gFitFunction");
966
967
968 //Normalization
969 gFitFunction->SetParName(0,"N1");
970 //near side peak
971 gFitFunction->SetParName(1,"N_{near side}");gFitFunction->FixParameter(1,0.0);
972 gFitFunction->SetParName(2,"Sigma_{near side}(delta eta)"); gFitFunction->FixParameter(2,0.0);
973 gFitFunction->SetParName(3,"Sigma_{near side}(delta phi)"); gFitFunction->FixParameter(3,0.0);
974 gFitFunction->SetParName(4,"Exponent_{near side}"); gFitFunction->FixParameter(4,1.0);
975 //away side ridge
976 gFitFunction->SetParName(5,"N_{away side}"); gFitFunction->FixParameter(5,0.0);
977 gFitFunction->SetParName(6,"Sigma_{away side}(delta phi)"); gFitFunction->FixParameter(6,0.0);
978 gFitFunction->SetParName(7,"Exponent_{away side}"); gFitFunction->FixParameter(7,1.0);
979 //longitudinal ridge
980 gFitFunction->SetParName(8,"N_{long. ridge}"); gFitFunction->SetParameter(8,0.05);//
981 gFitFunction->FixParameter(8,0.0);
982 gFitFunction->SetParName(9,"Sigma_{long. ridge}(delta eta)"); gFitFunction->FixParameter(9,0.0);
983 gFitFunction->SetParName(10,"Exponent_{long. ridge}"); gFitFunction->FixParameter(10,1.0);
984 //wing structures
985 gFitFunction->SetParName(11,"N_{wing}"); gFitFunction->FixParameter(11,0.0);
986 //flow contribution
987 gFitFunction->SetParName(12,"N_{flow}"); gFitFunction->SetParameter(12,0.2);gFitFunction->SetParLimits(12,0.0,10);
988 gFitFunction->SetParName(13,"V1"); gFitFunction->SetParameter(13,0.005);gFitFunction->SetParLimits(13,0.0,10);
989 gFitFunction->SetParName(14,"V2"); gFitFunction->SetParameter(14,0.1);gFitFunction->SetParLimits(14,0.0,10);
990 gFitFunction->SetParName(15,"V3"); gFitFunction->SetParameter(15,0.05);gFitFunction->SetParLimits(15,0.0,10);
991 gFitFunction->SetParName(16,"V4"); gFitFunction->SetParameter(16,0.005);gFitFunction->SetParLimits(16,0.0,10);
992 gFitFunction->SetParName(17,"Offset"); gFitFunction->FixParameter(17,0.0);
993
994 // flow parameters
995 Double_t fNV = 0.;
996 Double_t fV1 = 0.;
997 Double_t fV2 = 0.;
998 Double_t fV3 = 0.;
999 Double_t fV4 = 0.;
1000
1001 //Fitting the correlation function (first the edges to extract flow)
1002 gHist->Fit("gFitFunction","nm","",1.0,1.6);
1003
1004 fNV += gFitFunction->GetParameter(12);
1005 fV1 += gFitFunction->GetParameter(13);
1006 fV2 += gFitFunction->GetParameter(14);
1007 fV3 += gFitFunction->GetParameter(15);
1008 fV4 += gFitFunction->GetParameter(16);
1009
1010 gHist->Fit("gFitFunction","nm","",-1.6,-1.0);
1011
1012 fNV += gFitFunction->GetParameter(12);
1013 fV1 += gFitFunction->GetParameter(13);
1014 fV2 += gFitFunction->GetParameter(14);
1015 fV3 += gFitFunction->GetParameter(15);
1016 fV4 += gFitFunction->GetParameter(16);
1017
1018 // Now fit the whole with fixed flow
1019 gFitFunction->FixParameter(12,fNV/2.);
1020 gFitFunction->FixParameter(13,fV1/2.);
1021 gFitFunction->FixParameter(14,fV2/2.);
1022 gFitFunction->FixParameter(15,fV3/2.);
1023 gFitFunction->FixParameter(16,fV4/2.);
1024
1025 gFitFunction->ReleaseParameter(0);gFitFunction->SetParameter(0,1.0);
1026 gFitFunction->ReleaseParameter(1);gFitFunction->SetParameter(1,0.3);
1027 gFitFunction->ReleaseParameter(2);gFitFunction->SetParameter(2,0.3);gFitFunction->SetParLimits(2,0.05,0.7);
1028 gFitFunction->ReleaseParameter(3);gFitFunction->SetParameter(3,0.3);gFitFunction->SetParLimits(3,0.05,1.7);
1029 //gFitFunction->ReleaseParameter(4);gFitFunction->SetParameter(4,1.0);gFitFunction->SetParLimits(4,0.0,2.0);
1030 gFitFunction->ReleaseParameter(5);gFitFunction->SetParameter(5,1.0);//gFitFunction->SetParLimits(5,0.0,10);
1031 gFitFunction->ReleaseParameter(6);gFitFunction->SetParameter(6,0.5);//gFitFunction->SetParLimits(6,0.0,10);
1032 //gFitFunction->ReleaseParameter(7);gFitFunction->SetParameter(7,1.0);gFitFunction->SetParLimits(7,0.0,2.0);
1033 gFitFunction->ReleaseParameter(8);gFitFunction->SetParameter(8,0.05);
1034 gFitFunction->ReleaseParameter(9);gFitFunction->SetParameter(9,0.6);gFitFunction->SetParLimits(9,0.1,10.0);
1035 //gFitFunction->ReleaseParameter(10);gFitFunction->SetParameter(10,1.0);gFitFunction->SetParLimits(10,0.0,2.0);
1036 gFitFunction->ReleaseParameter(17);gFitFunction->SetParameter(17,0.7);gFitFunction->SetParLimits(17,0.0,0.9);
1037
1038 gHist->Fit("gFitFunction","nm");
1039
1040
1041 //Cloning the histogram
1042 TString histoName = gHist->GetName(); histoName += "Fit";
1043 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());
1044 TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
1045 gHistResidual->SetName("gHistResidual");
1046 gHistResidual->Sumw2();
1047
1048 for(Int_t iBinDeltaEta = 1; iBinDeltaEta <= gHist->GetNbinsX(); iBinDeltaEta++) {
1049 for(Int_t iBinDeltaPhi = 1; iBinDeltaPhi <= gHist->GetNbinsY(); iBinDeltaPhi++) {
1050 gHistFit->SetBinContent(iBinDeltaEta,iBinDeltaPhi,gFitFunction->Eval(gHist->GetXaxis()->GetBinCenter(iBinDeltaEta),gHist->GetYaxis()->GetBinCenter(iBinDeltaPhi)));
1051 }
1052 }
1053 gHistResidual->Add(gHistFit,-1);
1054
1055 //Write to output file
1056 TString newFileName = "correlationFunctionFit";
1057 if(histoName.Contains("PN")) newFileName += "PN";
1058 else if(histoName.Contains("NP")) newFileName += "NP";
1059 else if(histoName.Contains("PP")) newFileName += "PP";
1060 else if(histoName.Contains("NN")) newFileName += "NN";
1061 newFileName += ".Centrality";
1062 newFileName += gCentrality; newFileName += ".Psi";
1063 if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
1064 else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
1065 else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
1066 else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
1067 else newFileName += "All.PttFrom";
1068 newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
1069 newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
1070 newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
1071 newFileName += Form("%.1f",ptAssociatedMax);
1072 newFileName += ".root";
1073 TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
1074 gHist->Write();
1075 gHistFit->Write();
1076 gHistResidual->Write();
1077 gFitFunction->Write();
1078 newFile->Close();
1079
1080
1081}
1082
1083
1084
1085// //____________________________________________________________//
1086// void fitCorrelationFunctions(Int_t gCentrality = 1,
1087// Double_t psiMin = -0.5, Double_t psiMax = 3.5,
1088// Double_t ptTriggerMin = -1.,
1089// Double_t ptTriggerMax = -1.,
1090// Double_t ptAssociatedMin = -1.,
1091// Double_t ptAssociatedMax = -1.,
1092// TH2D *gHist) {
1093
1094// cout<<"FITTING FUNCTION (HOUSTON)"<<endl;
1095
1096// // Fit Function
1097// //x axis = delta_eta
1098// //y axis = delta_phi
1099// // [9]*exp(-1*pow(((x/[10])^2 + (y/[10])^2),0.5)) Hyper expo
1100// TF2 *fit1 = new TF2("fit1","[0] + [1]*cos(y) + [2]*cos(2*y) + [3]*cos(3*y) + [4]*cos(4*y) +[5]*cos(5*y)+ [6]*exp(-0.5*pow(((x/[7])^2 + (y/[8])^2),[11])) + [6]*exp(-0.5*pow(((x/[7])^2 + ((y-6.283)/[8])^2),[11]))+ [9]*exp(-1*((x/[10])^2 + (y/[10])^2)) ",-2.0,2.0,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
1101
1102
1103
1104// Double_t Parameters[] = {0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.3,0.3,1.0,0.1,0.1};
1105
1106// fit1->SetParameters(Parameters); // input pars from macro arguments
1107
1108// fit1->SetParName(0,"offset");
1109// fit1->SetParName(1,"v1");
1110// fit1->SetParName(2,"v2");
1111// fit1->SetParName(3,"v3");
1112// fit1->SetParName(4,"v4");
1113// fit1->SetParName(5,"v5");
1114// fit1->SetParName(6,"Ridgeamp");
1115// fit1->SetParName(7,"Ridgesigx");
1116// fit1->SetParName(8,"Ridgesigy");
1117// fit1->SetParName(9,"Expoamp");
1118// fit1->SetParName(10,"Exposig");
1119// fit1->SetParName(11,"Gausspara");
1120
1121
1122// //Fit Parameter Ranges
1123// fit1->SetParLimits(0,-2.0,2.0); //offset
1124// fit1->SetParLimits(1,-1.0,0.1); //v1
1125// fit1->SetParLimits(2,-1.6,0.9); //v2
1126// fit1->SetParLimits(3,0.0,0.5); //v3
1127// fit1->SetParLimits(4,0.0,0.9); //v4
1128// fit1->SetParLimits(5,0.0,0.9); //v5
1129// fit1->SetParLimits(6,0.0,1.5); //Ridgeamp
1130// fit1->SetParLimits(7,0.1,3.0); //Ridgesigx
1131// fit1->SetParLimits(8,0.1,2.0); //Ridgesigy
1132// fit1->SetParLimits(9,0.0,6.0); //Expoamp
1133// fit1->SetParLimits(10,0.05,0.5); //Exposig
1134// fit1->SetParLimits(11,0.0,2.0); //Gausspara
1135
1136// //Fitting the correlation function
1137// gHist->Fit("fit1","nm");
1138
1139// //Cloning the histogram
1140// TString histoName = gHist->GetName(); histoName += "Fit";
1141// 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());
1142// TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
1143// gHistResidual->SetName("gHistResidual");
1144// gHistResidual->Sumw2();
1145
1146// for(Int_t iBinDeltaEta = 1; iBinDeltaEta <= gHist->GetNbinsX(); iBinDeltaEta++) {
1147// for(Int_t iBinDeltaPhi = 1; iBinDeltaPhi <= gHist->GetNbinsY(); iBinDeltaPhi++) {
1148// gHistFit->SetBinContent(iBinDeltaEta,iBinDeltaPhi,fit1->Eval(gHist->GetXaxis()->GetBinCenter(iBinDeltaEta),gHist->GetYaxis()->GetBinCenter(iBinDeltaPhi)));
1149// }
1150// }
1151// gHistResidual->Add(gHistFit,-1);
1152
1153// //Write to output file
1154// TString newFileName = "correlationFunctionFit";
1155// if(histoName.Contains("PN")) newFileName += "PN";
1156// else if(histoName.Contains("NP")) newFileName += "NP";
1157// else if(histoName.Contains("PP")) newFileName += "PP";
1158// else if(histoName.Contains("NN")) newFileName += "NN";
1159// newFileName += ".Centrality";
1160// newFileName += gCentrality; newFileName += ".Psi";
1161// if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
1162// else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
1163// else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
1164// else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
1165// else newFileName += "All.PttFrom";
1166// newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
1167// newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
1168// newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
1169// newFileName += Form("%.1f",ptAssociatedMax);
1170// newFileName += ".root";
1171// TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
1172// gHist->Write();
1173// gHistFit->Write();
1174// gHistResidual->Write();
1175// fit1->Write();
1176// newFile->Close();
1177
1178
1179// }