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