+//____________________________________________________________//
+void drawCorrelationFunctions(const char* lhcPeriod = "LHC11h",
+ Int_t gTrainID = 208,
+ Int_t gCentrality = 1,
+ Double_t psiMin = -0.5, Double_t psiMax = 3.5,
+ Double_t ptTriggerMin = -1.,
+ Double_t ptTriggerMax = -1.,
+ Double_t ptAssociatedMin = -1.,
+ Double_t ptAssociatedMax = -1.) {
+ //Macro that draws the charge dependent correlation functions
+ //for each centrality bin for the different pT of trigger and
+ //associated particles
+ //Author: Panos.Christakoglou@nikhef.nl
+ TGaxis::SetMaxDigits(3);
+
+ //Get the input file
+ TString filename = "PbPb/"; filename += lhcPeriod;
+ filename +="/Train"; filename += gTrainID;
+ filename +="/Centrality"; filename += gCentrality;
+ filename += "/correlationFunction.Centrality";
+ filename += gCentrality; filename += ".Psi";
+ if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
+ else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
+ else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
+ else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
+ else filename += "All.PttFrom";
+ filename += Form("%.1f",ptTriggerMin); filename += "To";
+ filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
+ filename += Form("%.1f",ptAssociatedMin); filename += "To";
+ filename += Form("%.1f",ptAssociatedMax);
+ filename += ".root";
+
+ //Open the file
+ TFile *f = TFile::Open(filename.Data());
+ if((!f)||(!f->IsOpen())) {
+ Printf("The file %s is not found. Aborting...",filename);
+ return listBF;
+ }
+ //f->ls();
+
+ //Latex
+ TString centralityLatex = "Centrality: ";
+ centralityLatex += centralityArray[gCentrality-1];
+ centralityLatex += "%";
+
+ TString psiLatex;
+ if((psiMin == -0.5)&&(psiMax == 0.5))
+ psiLatex = " -7.5^{o} < #varphi - #Psi_{2} < 7.5^{o}";
+ else if((psiMin == 0.5)&&(psiMax == 1.5))
+ psiLatex = " 37.5^{o} < #varphi - #Psi_{2} < 52.5^{o}";
+ else if((psiMin == 1.5)&&(psiMax == 2.5))
+ psiLatex = " 82.5^{o} < #varphi - #Psi_{2} < 97.5^{o}";
+ else
+ psiLatex = " 0^{o} < #varphi - #Psi_{2} < 180^{o}";
+
+ TString pttLatex = Form("%.1f",ptTriggerMin);
+ pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
+ pttLatex += " GeV/c";
+
+ TString ptaLatex = Form("%.1f",ptAssociatedMin);
+ ptaLatex += " < p_{T,assoc} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
+ ptaLatex += " GeV/c";
+
+ TLatex *latexInfo1 = new TLatex();
+ latexInfo1->SetNDC();
+ latexInfo1->SetTextSize(0.045);
+ latexInfo1->SetTextColor(1);
+
+ TString pngName;
+
+ //============================================================//
+ //Get the +- correlation function
+ TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
+ gHistPN->SetStats(kFALSE);
+ gHistPN->SetTitle("");
+ gHistPN->GetXaxis()->SetRangeUser(-1.45,1.45);
+ gHistPN->GetXaxis()->CenterTitle();
+ gHistPN->GetXaxis()->SetTitleOffset(1.2);
+ gHistPN->GetYaxis()->CenterTitle();
+ gHistPN->GetYaxis()->SetTitleOffset(1.2);
+ gHistPN->GetZaxis()->SetTitleOffset(1.2);
+ TCanvas *cPN = new TCanvas("cPN","",0,0,600,500);
+ cPN->SetFillColor(10); cPN->SetHighLightColor(10);
+ cPN->SetLeftMargin(0.15);
+ gHistPN->DrawCopy("surf1fb");
+ gPad->SetTheta(30); // default is 30
+ gPad->SetPhi(-60); // default is 30
+ gPad->Update();
+
+ latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
+ //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
+ latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
+ latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
+
+ pngName = "CorrelationFunction.Centrality";
+ pngName += centralityArray[gCentrality-1];
+ pngName += ".Psi";
+ if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
+ else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
+ else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
+ else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
+ else pngName += "All.PttFrom";
+ pngName += Form("%.1f",ptTriggerMin); pngName += "To";
+ pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
+ pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
+ pngName += Form("%.1f",ptAssociatedMax);
+ pngName += ".PositiveNegative.png";
+ cPN->SaveAs(pngName.Data());
+ fitCorrelationFunctions(gCentrality, psiMin, psiMax,
+ ptTriggerMin,ptTriggerMax,
+ ptAssociatedMin, ptAssociatedMax,gHistPN);
+ //============================================================//
+ //Get the -+ correlation function
+ TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
+ gHistNP->SetStats(kFALSE);
+ gHistNP->SetTitle("");
+ gHistNP->GetXaxis()->SetRangeUser(-1.45,1.45);
+ gHistNP->GetXaxis()->CenterTitle();
+ gHistNP->GetXaxis()->SetTitleOffset(1.2);
+ gHistNP->GetYaxis()->CenterTitle();
+ gHistNP->GetYaxis()->SetTitleOffset(1.2);
+ gHistNP->GetZaxis()->SetTitleOffset(1.2);
+ TCanvas *cNP = new TCanvas("cNP","",50,50,600,500);
+ cNP->SetFillColor(10); cNP->SetHighLightColor(10);
+ cNP->SetLeftMargin(0.15);
+ gHistNP->DrawCopy("surf1fb");
+ gPad->SetTheta(30); // default is 30
+ gPad->SetPhi(-60); // default is 30
+ gPad->Update();
+
+ latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
+ //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
+ latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
+ latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
+
+ pngName = "CorrelationFunction.Centrality";
+ pngName += centralityArray[gCentrality-1];
+ pngName += ".Psi";
+ if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
+ else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
+ else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
+ else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
+ else pngName += "All.PttFrom";
+ pngName += Form("%.1f",ptTriggerMin); pngName += "To";
+ pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
+ pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
+ pngName += Form("%.1f",ptAssociatedMax);
+ pngName += ".NegativePositive.png";
+ cNP->SaveAs(pngName.Data());
+
+ //============================================================//
+ //Get the ++ correlation function
+ TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
+ gHistPP->SetStats(kFALSE);
+ gHistPP->SetTitle("");
+ gHistPP->GetXaxis()->SetRangeUser(-1.45,1.45);
+ gHistPP->GetXaxis()->CenterTitle();
+ gHistPP->GetXaxis()->SetTitleOffset(1.2);
+ gHistPP->GetYaxis()->CenterTitle();
+ gHistPP->GetYaxis()->SetTitleOffset(1.2);
+ gHistPP->GetZaxis()->SetTitleOffset(1.2);
+ TCanvas *cPP = new TCanvas("cPP","",100,100,600,500);
+ cPP->SetFillColor(10); cPP->SetHighLightColor(10);
+ cPP->SetLeftMargin(0.15);
+ gHistPP->DrawCopy("surf1fb");
+ gPad->SetTheta(30); // default is 30
+ gPad->SetPhi(-60); // default is 30
+ gPad->Update();
+
+ latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
+ //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
+ latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
+ latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
+
+ pngName = "CorrelationFunction.Centrality";
+ pngName += centralityArray[gCentrality-1];
+ pngName += ".Psi";
+ if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
+ else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
+ else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
+ else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
+ else pngName += "All.PttFrom";
+ pngName += Form("%.1f",ptTriggerMin); pngName += "To";
+ pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
+ pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
+ pngName += Form("%.1f",ptAssociatedMax);
+ pngName += ".PositivePositive.png";
+ cPP->SaveAs(pngName.Data());
+
+ //============================================================//
+ //Get the -- correlation function
+ TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
+ gHistNN->SetStats(kFALSE);
+ gHistNN->SetTitle("");
+ gHistNN->GetXaxis()->SetRangeUser(-1.45,1.45);
+ gHistNN->GetXaxis()->CenterTitle();
+ gHistNN->GetXaxis()->SetTitleOffset(1.2);
+ gHistNN->GetYaxis()->CenterTitle();
+ gHistNN->GetYaxis()->SetTitleOffset(1.2);
+ gHistNN->GetZaxis()->SetTitleOffset(1.2);
+ TCanvas *cNN = new TCanvas("cNN","",150,150,600,500);
+ cNN->SetFillColor(10); cNN->SetHighLightColor(10);
+ cNN->SetLeftMargin(0.15);
+ gHistNN->DrawCopy("surf1fb");
+ gPad->SetTheta(30); // default is 30
+ gPad->SetPhi(-60); // default is 30
+ gPad->Update();
+
+ latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
+ //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
+ latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
+ latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
+
+ pngName = "CorrelationFunction.Centrality";
+ pngName += centralityArray[gCentrality-1];
+ pngName += ".Psi";
+ if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
+ else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
+ else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
+ else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
+ else pngName += "All.PttFrom";
+ pngName += Form("%.1f",ptTriggerMin); pngName += "To";
+ pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
+ pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
+ pngName += Form("%.1f",ptAssociatedMax);
+ pngName += ".NegativeNegative.png";
+ cNN->SaveAs(pngName.Data());
+}
+
+//____________________________________________________________//
+void fitCorrelationFunctions(Int_t gCentrality = 1,
+ Double_t psiMin = -0.5, Double_t psiMax = 3.5,
+ Double_t ptTriggerMin = -1.,
+ Double_t ptTriggerMax = -1.,
+ Double_t ptAssociatedMin = -1.,
+ Double_t ptAssociatedMax = -1.,
+ TH2D *gHist) {
+
+ cout<<"FITTING FUNCTION"<<endl;
+
+ //near side peak: [1]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4]))
+ //away side ridge: [5]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((y-TMath::Pi())/[6]),2)),[7]))
+ //longitudinal ridge: [8]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[9]),2)),[10]))
+ //wing structures: [11]*TMath::Power(x,2)
+ //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))
+ 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/[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.);
+ gFitFunction->SetName("gFitFunction");
+ //Normalization
+ gFitFunction->SetParName(0,"N1"); gFitFunction->SetParameter(0,1.0);
+ //near side peak
+ gFitFunction->SetParName(1,"N_{near side}"); gFitFunction->SetParameter(1,0.3);
+ gFitFunction->SetParName(2,"Sigma_{near side}(delta eta)"); gFitFunction->SetParameter(2,0.3);
+ gFitFunction->SetParName(3,"Sigma_{near side}(delta phi)"); gFitFunction->SetParameter(3,0.1);
+ gFitFunction->SetParName(4,"Exponent_{near side}"); gFitFunction->SetParameter(4,1.1);
+ //away side ridge
+ gFitFunction->SetParName(5,"N_{away side}"); gFitFunction->SetParameter(5,0.1);
+ gFitFunction->SetParName(6,"Sigma_{away side}(delta phi)"); gFitFunction->SetParameter(6,1.1);
+ gFitFunction->SetParName(7,"Exponent_{away side}"); gFitFunction->SetParameter(7,1.0);
+ //longitudianl ridge
+ gFitFunction->SetParName(8,"N_{long. ridge}"); gFitFunction->SetParameter(8,0.05);
+ gFitFunction->SetParName(9,"Sigma_{long. ridge}(delta eta)"); gFitFunction->SetParameter(9,0.6);
+ gFitFunction->SetParName(10,"Exponent_{long. ridge}"); gFitFunction->SetParameter(10,1.0);
+ //wing structures
+ gFitFunction->SetParName(11,"N_{wing}"); gFitFunction->SetParameter(11,0.01);
+ //flow contribution
+ gFitFunction->SetParName(12,"N_{flow}"); gFitFunction->SetParameter(12,0.2);
+ gFitFunction->SetParName(13,"V1"); gFitFunction->SetParameter(13,0.005);
+ gFitFunction->SetParName(14,"V2"); gFitFunction->SetParameter(14,0.1);
+ gFitFunction->SetParName(15,"V3"); gFitFunction->SetParameter(15,0.05);
+ gFitFunction->SetParName(16,"V4"); gFitFunction->SetParameter(16,0.005);
+
+ //Fitting the correlation function
+ gHist->Fit("gFitFunction","nm");
+
+ //Cloning the histogram
+ TString histoName = gHist->GetName(); histoName += "Fit";
+ 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());
+ TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
+ gHistResidual->SetName("gHistResidual");
+ gHistResidual->Sumw2();
+
+ for(Int_t iBinDeltaEta = 1; iBinDeltaEta <= gHist->GetNbinsX(); iBinDeltaEta++) {
+ for(Int_t iBinDeltaPhi = 1; iBinDeltaPhi <= gHist->GetNbinsY(); iBinDeltaPhi++) {
+ gHistFit->SetBinContent(iBinDeltaEta,iBinDeltaPhi,gFitFunction->Eval(gHist->GetXaxis()->GetBinCenter(iBinDeltaEta),gHist->GetYaxis()->GetBinCenter(iBinDeltaPhi)));
+ }
+ }
+ gHistResidual->Add(gHistFit,-1);
+
+ //Write to output file
+ TString newFileName = "correlationFunctionFit";
+ if(histoName.Contains("PN")) newFileName += "PN";
+ else if(histoName.Contains("NP")) newFileName += "NP";
+ else if(histoName.Contains("PP")) newFileName += "PP";
+ else if(histoName.Contains("NN")) newFileName += "NN";
+ newFileName += ".Centrality";
+ newFileName += gCentrality; newFileName += ".Psi";
+ if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
+ else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
+ else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
+ else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
+ else newFileName += "All.PttFrom";
+ newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
+ newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
+ newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
+ newFileName += Form("%.1f",ptAssociatedMax);
+ newFileName += ".root";
+ TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
+ gHist->Write();
+ gHistFit->Write();
+ gHistResidual->Write();
+ gFitFunction->Write();
+ newFile->Close();
+
+
+}