adding some drawing macros for balance functions
[u/mrichter/AliRoot.git] / PWGCF / EBYE / macros / drawCorrelationFunctionPsiChargeIndependent.C
1 const Int_t numberOfCentralityBins = 9;
2 TString centralityArray[numberOfCentralityBins] = {"0-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-100"};
3
4 const Int_t gRebin = 1;
5
6 void 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 //______________________________________________________//
54 TList *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 //______________________________________________________//
212 void 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
429     gHist[2] = bMixed->GetCorrelationFunctionChargeIndependent(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,normToTrig);
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 //____________________________________________________________//
536 void 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 //____________________________________________________________//
911 void 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 // }