PWGCF/EBYE/macros/drawCorrelationFunctionPsi.C

   1 const Int_t numberOfCentralityBins = 12;
   2 //TString centralityArray[numberOfCentralityBins] = {"0-4","4-5","6-14","30-40","40-50","50-60","60-70","70-80","0-100","0-1","1-2","2-3"};
   3 TString centralityArray[numberOfCentralityBins] = {"0-100","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-100","0-1","1-2","2-3"};
   4
   5 const Int_t gRebin = 1;
   6
   7 void drawCorrelationFunctionPsi(const char* filename = "AnalysisResultsPsi_train_06feb.root",
   8                                 Int_t gCentrality = 1,
   9                                 Int_t gBit = -1,
  10                                 const char* gCentralityEstimator = 0x0,
  11                                 Bool_t kShowShuffled = kFALSE,
  12                                 Bool_t kShowMixed = kTRUE,
  13                                 Double_t psiMin = -0.5,
  14                                 Double_t psiMax = 3.5,
  15                                 Double_t vertexZMin = -10.,
  16                                 Double_t vertexZMax = 10.,
  17                                 Double_t ptTriggerMin = -1.,
  18                                 Double_t ptTriggerMax = -1.,
  19                                 Double_t ptAssociatedMin = -1.,
  20                                 Double_t ptAssociatedMax = -1.,
  21                                 Bool_t normToTrig = kTRUE,
  22                                 Bool_t kUseVzBinning = kTRUE,
  23                                 Int_t rebinEta = 1,
  24                                 Int_t rebinPhi = 1,
  25                                 TString eventClass = "EventPlane", //Can be "EventPlane", "Centrality", "Multiplicity"
  26                                 Bool_t bToy = kFALSE,
  27                                 TString listNameAdd = "")
  28 {
  29   //Macro that draws the correlation functions from the balance function
  30   //analysis vs the reaction plane
  31   //Author: Panos.Christakoglou@nikhef.nl
  32   //gROOT->LoadMacro("~/SetPlotStyle.C");
  33   //SetPlotStyle();
  34   gStyle->SetPalette(1,0);
  35
  36   //Load the PWG2ebye library
  37   gSystem->Load("libCore.so");
  38   gSystem->Load("libGeom.so");
  39   gSystem->Load("libVMC.so");
  40   gSystem->Load("libPhysics.so");
  41   gSystem->Load("libTree.so");
  42   gSystem->Load("libSTEERBase.so");
  43   gSystem->Load("libESD.so");
  44   gSystem->Load("libAOD.so");
  45
  46   gSystem->Load("libANALYSIS.so");
  47   gSystem->Load("libANALYSISalice.so");
  48   gSystem->Load("libEventMixing.so");
  49   gSystem->Load("libCORRFW.so");
  50   gSystem->Load("libPWGTools.so");
  51   gSystem->Load("libPWGCFebye.so");
  52
  53   //Prepare the objects and return them
  54   TList *listQA = NULL;
  55   TList *list = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,0,bToy,listNameAdd);
  56   TList *listShuffled = NULL;
  57   if(kShowShuffled) listShuffled = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,1,bToy,listNameAdd);
  58   TList *listMixed = NULL;
  59   if(kShowMixed) listMixed = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,2,bToy,listNameAdd);
  60
  61   // else  get the QA histograms (for convolution)
  62   else{
  63     //Open the file again
  64     TFile *f = TFile::Open(filename,"UPDATE");
  65     if((!f)||(!f->IsOpen())) {
  66       Printf("The file %s is not found.",filename);
  67     }
  68
  69
  70     TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
  71     if(!dir) {
  72       Printf("The TDirectoryFile is not found.",filename);
  73     }
  74
  75     TString listQAName = "listQAPsi_";
  76
  77     listQAName += centralityArray[gCentrality-1];
  78     if(gBit > -1) {
  79       listQAName += "_Bit"; listQAName += gBit; }
  80     if(gCentralityEstimator) {
  81       listQAName += "_"; listQAName += gCentralityEstimator;}
  82
  83     listQA = (TList*)dir->Get(Form("%s",listQAName.Data()));
  84     if(!listQA) {
  85       Printf("TList QA not found!!!");
  86     }
  87   }
  88
  89   if(!list) {
  90     Printf("The TList object was not created");
  91     return;
  92   }
  93   else
  94     draw(list,listShuffled,listMixed,listQA,
  95          gCentralityEstimator,gCentrality,psiMin,psiMax,vertexZMin,vertexZMax,
  96          ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,normToTrig,kUseVzBinning,rebinEta,rebinPhi,eventClass);
  97 }
  98
  99 //______________________________________________________//
 100 TList *GetListOfObjects(const char* filename,
 101                         Int_t gCentrality,
 102                         Int_t gBit,
 103                         const char *gCentralityEstimator,
 104                         Int_t kData = 1,
 105                         Bool_t bToy = kFALSE,
 106                         TString listNameAdd = "") {
 107   //Get the TList objects (QA, bf, bf shuffled)
 108   TList *listBF = 0x0;
 109
 110   //Open the file
 111   TFile *f = TFile::Open(filename,"UPDATE");
 112   if((!f)||(!f->IsOpen())) {
 113     Printf("The file %s is not found. Aborting...",filename);
 114     return listBF;
 115   }
 116   //f->ls();
 117
 118   TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
 119   if(!dir) {
 120     Printf("The TDirectoryFile is not found. Aborting...",filename);
 121     return listBF;
 122   }
 123   //dir->ls();
 124
 125   TString listBFName;
 126   if(kData == 0) {
 127     //cout<<"no shuffling - no mixing"<<endl;
 128     listBFName = "listBFPsi";
 129   }
 130   else if(kData == 1) {
 131     //cout<<"shuffling - no mixing"<<endl;
 132     listBFName = "listBFPsiShuffled";
 133   }
 134   else if(kData == 2) {
 135     //cout<<"no shuffling - mixing"<<endl;
 136     listBFName = "listBFPsiMixed";
 137   }
 138
 139   // different list names in case of toy model
 140   if(!bToy){
 141     listBFName += "_";
 142     listBFName += centralityArray[gCentrality-1];
 143     if(gBit > -1) {
 144       listBFName += "_Bit"; listBFName += gBit; }
 145     if(gCentralityEstimator) {
 146       listBFName += "_"; listBFName += gCentralityEstimator;}
 147   }
 148   else{
 149     listBFName.ReplaceAll("Psi","");
 150   }
 151   listBFName += listNameAdd;
 152
 153   // histograms were already retrieved (in first iteration)
 154   if(dir->Get(Form("%s_histograms",listBFName.Data()))){
 155     listBF = dynamic_cast<TList *>(dir->Get(Form("%s_histograms",listBFName.Data())));
 156   }
 157
 158   // histograms were not yet retrieved (this is the first iteration)
 159   else{
 160
 161     listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
 162     if(!listBF) dir->ls();
 163     cout<<"======================================================="<<endl;
 164     cout<<"List name (control): "<<listBFName.Data()<<endl;
 165     cout<<"List name: "<<listBF->GetName()<<endl;
 166     //listBF->ls();
 167
 168     //Get the histograms
 169     TString histoName;
 170     if(kData == 0)
 171       histoName = "fHistP";
 172     else if(kData == 1)
 173       histoName = "fHistP_shuffle";
 174     else if(kData == 2)
 175       histoName = "fHistP";
 176     if(gCentralityEstimator)
 177       histoName += gCentralityEstimator;
 178     AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 179     if(!fHistP) {
 180       Printf("fHistP %s not found!!!",histoName.Data());
 181       break;
 182     }
 183     fHistP->FillParent(); fHistP->DeleteContainers();
 184
 185     if(kData == 0)
 186       histoName = "fHistN";
 187     if(kData == 1)
 188       histoName = "fHistN_shuffle";
 189     if(kData == 2)
 190       histoName = "fHistN";
 191     if(gCentralityEstimator)
 192       histoName += gCentralityEstimator;
 193     AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 194     if(!fHistN) {
 195       Printf("fHistN %s not found!!!",histoName.Data());
 196       break;
 197     }
 198     fHistN->FillParent(); fHistN->DeleteContainers();
 199
 200     if(kData == 0)
 201       histoName = "fHistPN";
 202     if(kData == 1)
 203       histoName = "fHistPN_shuffle";
 204     if(kData == 2)
 205       histoName = "fHistPN";
 206     if(gCentralityEstimator)
 207       histoName += gCentralityEstimator;
 208     AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 209     if(!fHistPN) {
 210       Printf("fHistPN %s not found!!!",histoName.Data());
 211       break;
 212     }
 213     fHistPN->FillParent(); fHistPN->DeleteContainers();
 214
 215     if(kData == 0)
 216       histoName = "fHistNP";
 217     if(kData == 1)
 218       histoName = "fHistNP_shuffle";
 219     if(kData == 2)
 220       histoName = "fHistNP";
 221     if(gCentralityEstimator)
 222       histoName += gCentralityEstimator;
 223     AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 224     if(!fHistNP) {
 225       Printf("fHistNP %s not found!!!",histoName.Data());
 226       break;
 227     }
 228     fHistNP->FillParent(); fHistNP->DeleteContainers();
 229
 230     if(kData == 0)
 231       histoName = "fHistPP";
 232     if(kData == 1)
 233       histoName = "fHistPP_shuffle";
 234     if(kData == 2)
 235       histoName = "fHistPP";
 236     if(gCentralityEstimator)
 237       histoName += gCentralityEstimator;
 238     AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 239     if(!fHistPP) {
 240       Printf("fHistPP %s not found!!!",histoName.Data());
 241       break;
 242     }
 243     fHistPP->FillParent(); fHistPP->DeleteContainers();
 244
 245     if(kData == 0)
 246       histoName = "fHistNN";
 247     if(kData == 1)
 248       histoName = "fHistNN_shuffle";
 249     if(kData == 2)
 250       histoName = "fHistNN";
 251     if(gCentralityEstimator)
 252       histoName += gCentralityEstimator;
 253     AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 254     if(!fHistNN) {
 255       Printf("fHistNN %s not found!!!",histoName.Data());
 256       break;
 257     }
 258     fHistNN->FillParent(); fHistNN->DeleteContainers();
 259
 260     dir->cd();
 261     listBF->Write(Form("%s_histograms",listBFName.Data()), TObject::kSingleKey);
 262
 263   }// first iteration
 264
 265   f->Close();
 266
 267   return listBF;
 268 }
 269
 270 //______________________________________________________//
 271 void draw(TList *list, TList *listBFShuffled, TList *listBFMixed,
 272           TList *listQA,
 273           const char *gCentralityEstimator,
 274           Int_t gCentrality, Double_t psiMin, Double_t psiMax,
 275           Double_t vertexZMin,Double_t vertexZMax,
 276           Double_t ptTriggerMin, Double_t ptTriggerMax,
 277           Double_t ptAssociatedMin, Double_t ptAssociatedMax,
 278           Bool_t normToTrig,
 279           Bool_t kUseVzBinning,
 280           Int_t rebinEta, Int_t rebinPhi,TString eventClass) {
 281   //Draws the correlation functions for every centrality bin
 282   //(+-), (-+), (++), (--)
 283   AliTHn *hP = NULL;
 284   AliTHn *hN = NULL;
 285   AliTHn *hPN = NULL;
 286   AliTHn *hNP = NULL;
 287   AliTHn *hPP = NULL;
 288   AliTHn *hNN = NULL;
 289
 290   TString gHistPname = "fHistP";
 291   if(gCentralityEstimator) gHistPname += gCentralityEstimator;
 292   TString gHistNname = "fHistN";
 293   if(gCentralityEstimator) gHistNname += gCentralityEstimator;
 294   TString gHistPNname = "fHistPN";
 295   if(gCentralityEstimator) gHistPNname += gCentralityEstimator;
 296   TString gHistNPname = "fHistNP";
 297   if(gCentralityEstimator) gHistNPname += gCentralityEstimator;
 298   TString gHistPPname = "fHistPP";
 299   if(gCentralityEstimator) gHistPPname += gCentralityEstimator;
 300   TString gHistNNname = "fHistNN";
 301   if(gCentralityEstimator) gHistNNname += gCentralityEstimator;
 302
 303   hP = (AliTHn*) list->FindObject(gHistPname.Data());
 304   hN = (AliTHn*) list->FindObject(gHistNname.Data());
 305   hPN = (AliTHn*) list->FindObject(gHistPNname.Data());
 306   hNP = (AliTHn*) list->FindObject(gHistNPname.Data());
 307   hPP = (AliTHn*) list->FindObject(gHistPPname.Data());
 308   hNN = (AliTHn*) list->FindObject(gHistNNname.Data());
 309   hNN->Print();
 310
 311
 312   //Create the AliBalancePsi object and fill it with the AliTHn objects
 313   AliBalancePsi *b = new AliBalancePsi();
 314   b->SetEventClass(eventClass);
 315   b->SetHistNp(hP);
 316   b->SetHistNn(hN);
 317   b->SetHistNpn(hPN);
 318   b->SetHistNnp(hNP);
 319   b->SetHistNpp(hPP);
 320   b->SetHistNnn(hNN);
 321   if(kUseVzBinning) b->SetVertexZBinning(kTRUE);
 322
 323   //balance function shuffling
 324   AliTHn *hPShuffled = NULL;
 325   AliTHn *hNShuffled = NULL;
 326   AliTHn *hPNShuffled = NULL;
 327   AliTHn *hNPShuffled = NULL;
 328   AliTHn *hPPShuffled = NULL;
 329   AliTHn *hNNShuffled = NULL;
 330   if(listBFShuffled) {
 331     //listBFShuffled->ls();
 332
 333     gHistPname = "fHistP_shuffle";
 334     if(gCentralityEstimator) gHistPname += gCentralityEstimator;
 335     gHistNname = "fHistN_shuffle";
 336     if(gCentralityEstimator) gHistNname += gCentralityEstimator;
 337     gHistPNname = "fHistPN_shuffle";
 338     if(gCentralityEstimator) gHistPNname += gCentralityEstimator;
 339     gHistNPname = "fHistNP_shuffle";
 340     if(gCentralityEstimator) gHistNPname += gCentralityEstimator;
 341     gHistPPname = "fHistPP_shuffle";
 342     if(gCentralityEstimator) gHistPPname += gCentralityEstimator;
 343     gHistNNname = "fHistNN_shuffle";
 344     if(gCentralityEstimator) gHistNNname += gCentralityEstimator;
 345
 346     hPShuffled = (AliTHn*) listBFShuffled->FindObject(gHistPname.Data());
 347     hPShuffled->SetName("gHistPShuffled");
 348     hNShuffled = (AliTHn*) listBFShuffled->FindObject(gHistNname.Data());
 349     hNShuffled->SetName("gHistNShuffled");
 350     hPNShuffled = (AliTHn*) listBFShuffled->FindObject(gHistPNname.Data());
 351     hPNShuffled->SetName("gHistPNShuffled");
 352     hNPShuffled = (AliTHn*) listBFShuffled->FindObject(gHistNPname.Data());
 353     hNPShuffled->SetName("gHistNPShuffled");
 354     hPPShuffled = (AliTHn*) listBFShuffled->FindObject(gHistPPname.Data());
 355     hPPShuffled->SetName("gHistPPShuffled");
 356     hNNShuffled = (AliTHn*) listBFShuffled->FindObject(gHistNNname.Data());
 357     hNNShuffled->SetName("gHistNNShuffled");
 358
 359     AliBalancePsi *bShuffled = new AliBalancePsi();
 360     bShuffled->SetEventClass(eventClass);
 361     bShuffled->SetHistNp(hPShuffled);
 362     bShuffled->SetHistNn(hNShuffled);
 363     bShuffled->SetHistNpn(hPNShuffled);
 364     bShuffled->SetHistNnp(hNPShuffled);
 365     bShuffled->SetHistNpp(hPPShuffled);
 366     bShuffled->SetHistNnn(hNNShuffled);
 367     if(kUseVzBinning) bShuffled->SetVertexZBinning(kTRUE);
 368   }
 369
 370   //balance function mixing
 371   AliTHn *hPMixed = NULL;
 372   AliTHn *hNMixed = NULL;
 373   AliTHn *hPNMixed = NULL;
 374   AliTHn *hNPMixed = NULL;
 375   AliTHn *hPPMixed = NULL;
 376   AliTHn *hNNMixed = NULL;
 377
 378   if(listBFMixed) {
 379     //listBFMixed->ls();
 380
 381     gHistPname = "fHistP";
 382     if(gCentralityEstimator) gHistPname += gCentralityEstimator;
 383     gHistNname = "fHistN";
 384     if(gCentralityEstimator) gHistNname += gCentralityEstimator;
 385     gHistPNname = "fHistPN";
 386     if(gCentralityEstimator) gHistPNname += gCentralityEstimator;
 387     gHistNPname = "fHistNP";
 388     if(gCentralityEstimator) gHistNPname += gCentralityEstimator;
 389     gHistPPname = "fHistPP";
 390     if(gCentralityEstimator) gHistPPname += gCentralityEstimator;
 391     gHistNNname = "fHistNN";
 392     if(gCentralityEstimator) gHistNNname += gCentralityEstimator;
 393     hPMixed = (AliTHn*) listBFMixed->FindObject(gHistPname.Data());
 394     hPMixed->SetName("gHistPMixed");
 395     hNMixed = (AliTHn*) listBFMixed->FindObject(gHistNname.Data());
 396     hNMixed->SetName("gHistNMixed");
 397     hPNMixed = (AliTHn*) listBFMixed->FindObject(gHistPNname.Data());
 398     hPNMixed->SetName("gHistPNMixed");
 399     hNPMixed = (AliTHn*) listBFMixed->FindObject(gHistNPname.Data());
 400     hNPMixed->SetName("gHistNPMixed");
 401     hPPMixed = (AliTHn*) listBFMixed->FindObject(gHistPPname.Data());
 402     hPPMixed->SetName("gHistPPMixed");
 403     hNNMixed = (AliTHn*) listBFMixed->FindObject(gHistNNname.Data());
 404     hNNMixed->SetName("gHistNNMixed");
 405
 406     AliBalancePsi *bMixed = new AliBalancePsi();
 407     bMixed->SetEventClass(eventClass);
 408     bMixed->SetHistNp(hPMixed);
 409     bMixed->SetHistNn(hNMixed);
 410     bMixed->SetHistNpn(hPNMixed);
 411     bMixed->SetHistNnp(hNPMixed);
 412     bMixed->SetHistNpp(hPPMixed);
 413     bMixed->SetHistNnn(hNNMixed);
 414     if(kUseVzBinning) bMixed->SetVertexZBinning(kTRUE);
 415   }
 416
 417
 418   TH2D *gHistPN[4];
 419   TH2D *gHistNP[4];
 420   TH2D *gHistPP[4];
 421   TH2D *gHistNN[4];
 422
 423   TCanvas *cPN[4];
 424   TCanvas *cNP[4];
 425   TCanvas *cPP[4];
 426   TCanvas *cNN[4];
 427   TString histoTitle, pngName;
 428
 429   // if no mixing then divide by convoluted histograms
 430   if(!listBFMixed && listQA){
 431
 432     if(!listQA->FindObject("fHistEtaPhiPos") || !listQA->FindObject("fHistEtaPhiNeg")){
 433
 434       Printf("fHistEtaPhiPos or fHistEtaPhiNeg not found! --> no convolution");
 435       listQA = NULL;
 436
 437     }
 438     else{
 439
 440       TH2D* fHistPos = (TH2D*)((TH3D*)listQA->FindObject("fHistEtaPhiPos"))->Project3D("xy");
 441       fHistPos->GetYaxis()->SetRangeUser(-0.79,0.79);
 442
 443       TH2D* fHistNeg = (TH2D*)((TH3D*)listQA->FindObject("fHistEtaPhiNeg"))->Project3D("xy");
 444       fHistNeg->GetYaxis()->SetRangeUser(-0.79,0.79);
 445
 446       gHistPN[2] = convolute2D(fHistPos, fHistNeg, "hConvPN");
 447       gHistPN[2]->Scale(1./gHistPN[2]->GetBinContent(gHistPN[2]->FindBin(0,0)));
 448
 449       gHistNP[2] = convolute2D(fHistNeg, fHistPos, "hConvNP");
 450       gHistNP[2]->Scale(1./gHistNP[2]->GetBinContent(gHistNP[2]->FindBin(0,0)));
 451
 452       gHistNN[2] = convolute2D(fHistNeg, fHistNeg, "hConvNN");
 453       gHistNN[2]->Scale(1./gHistNN[2]->GetBinContent(gHistNN[2]->FindBin(0,0)));
 454
 455       gHistPP[2] = convolute2D(fHistPos, fHistPos, "hConvPP");
 456       gHistPP[2]->Scale(1./gHistPP[2]->GetBinContent(gHistPP[2]->FindBin(0,0)));
 457     }
 458   }
 459
 460   //(+-)
 461   if(eventClass == "Centrality"){
 462     histoTitle = "(+-) | Centrality: ";
 463     histoTitle += psiMin;
 464     histoTitle += " - ";
 465     histoTitle += psiMax;
 466     histoTitle += " % ";
 467     histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
 468   }
 469   else if(eventClass == "Multiplicity"){
 470     histoTitle = "(+-) | Multiplicity: ";
 471     histoTitle += psiMin;
 472     histoTitle += " - ";
 473     histoTitle += psiMax;
 474     histoTitle += " tracks";
 475     histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
 476   }
 477   else{ // "EventPlane" (default)
 478     histoTitle = "(+-) | Centrality: ";
 479     histoTitle += centralityArray[gCentrality-1];
 480     histoTitle += "%";
 481     if((psiMin == -0.5)&&(psiMax == 0.5))
 482       histoTitle += " (-7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o})";
 483     else if((psiMin == 0.5)&&(psiMax == 1.5))
 484       histoTitle += " (37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o})";
 485     else if((psiMin == 1.5)&&(psiMax == 2.5))
 486       histoTitle += " (82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o})";
 487     else
 488       histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
 489   }
 490   gHistPN[0] = b->GetCorrelationFunctionPN(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 491   if(rebinEta > 1 || rebinPhi > 1){
 492     gHistPN[0]->Rebin2D(rebinEta,rebinPhi);
 493     gHistPN[0]->Scale(1./(Double_t)(rebinEta*rebinPhi));
 494   }
 495   gHistPN[0]->GetYaxis()->SetTitleOffset(1.5);
 496   gHistPN[0]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 497   gHistPN[0]->SetTitle(histoTitle.Data());
 498   cPN[0] = new TCanvas("cPN0","",0,0,600,500);
 499   cPN[0]->SetFillColor(10); cPN[0]->SetHighLightColor(10);
 500   gHistPN[0]->DrawCopy("surf1fb");
 501   gPad->SetTheta(30); // default is 30
 502   //gPad->SetPhi(130); // default is 30
 503   gPad->SetPhi(-60); // default is 30
 504   gPad->Update();
 505   pngName = "DeltaPhiDeltaEta.Centrality";
 506   pngName += centralityArray[gCentrality-1];
 507   pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 508   pngName += ".PositiveNegative.png";
 509   //cPN[0]->SaveAs(pngName.Data());
 510
 511   if(listBFShuffled) {
 512
 513     histoTitle.ReplaceAll("(+-)","(+-) shuffled");
 514
 515     gHistPN[1] = bShuffled->GetCorrelationFunctionPN(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 516     if(rebinEta > 1 || rebinPhi > 1){
 517       gHistPN[1]->Rebin2D(rebinEta,rebinPhi);
 518       gHistPN[1]->Scale(1./(Double_t)(rebinEta*rebinPhi));
 519     }
 520     gHistPN[1]->GetYaxis()->SetTitleOffset(1.5);
 521     gHistPN[1]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 522     gHistPN[1]->SetTitle(histoTitle.Data());
 523     cPN[1] = new TCanvas("cPN1","",0,100,600,500);
 524     cPN[1]->SetFillColor(10);
 525     cPN[1]->SetHighLightColor(10);
 526     gHistPN[1]->DrawCopy("surf1fb");
 527     gPad->SetTheta(30); // default is 30
 528     //gPad->SetPhi(130); // default is 30
 529     gPad->SetPhi(-60); // default is 30
 530     gPad->Update();
 531     pngName = "DeltaPhiDeltaEtaShuffled.Centrality";
 532     pngName += centralityArray[gCentrality-1];
 533     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 534     pngName += ".PositiveNegative.png";
 535     //cPN[1]->SaveAs(pngName.Data());
 536   }
 537
 538   if(listBFMixed) {
 539
 540     if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(+-) shuffled","(+-) mixed");
 541     else                                histoTitle.ReplaceAll("(+-)","(+-) mixed");
 542
 543     // if normalization to trigger then do not divide Event mixing by number of trigger particles
 544     gHistPN[2] = bMixed->GetCorrelationFunctionPN(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 545     if(rebinEta > 1 || rebinPhi > 1){
 546       gHistPN[2]->Rebin2D(rebinEta,rebinPhi);
 547       gHistPN[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));
 548     }
 549
 550     // normalization to 1 at (0,0) --> Jan Fietes method
 551     if(normToTrig){
 552       Double_t mixedNorm = gHistPN[2]->Integral(gHistPN[2]->GetXaxis()->FindBin(0-10e-5),gHistPN[2]->GetXaxis()->FindBin(0+10e-5),1,gHistPN[2]->GetNbinsX());
 553       mixedNorm /= gHistPN[2]->GetNbinsY()*(gHistPN[2]->GetXaxis()->FindBin(0.01) - gHistPN[2]->GetXaxis()->FindBin(-0.01) + 1);
 554       gHistPN[2]->Scale(1./mixedNorm);
 555     }
 556
 557     gHistPN[2]->GetYaxis()->SetTitleOffset(1.5);
 558     gHistPN[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 559     gHistPN[2]->SetTitle(histoTitle.Data());
 560     cPN[2] = new TCanvas("cPN2","",0,200,600,500);
 561     cPN[2]->SetFillColor(10);
 562     cPN[2]->SetHighLightColor(10);
 563     gHistPN[2]->DrawCopy("surf1fb");
 564     gPad->SetTheta(30); // default is 30
 565     //gPad->SetPhi(130); // default is 30
 566     gPad->SetPhi(-60); // default is 30
 567     gPad->Update();
 568     pngName = "DeltaPhiDeltaEtaMixed.Centrality";
 569     pngName += centralityArray[gCentrality-1];
 570     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 571     pngName += ".PositiveNegative.png";
 572     //cPN[2]->SaveAs(pngName.Data());
 573
 574     //Correlation function (+-)
 575     gHistPN[3] = b->GetCorrelationFunction("PN",psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 576     gHistPN[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
 577     if(normToTrig)
 578       gHistPN[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
 579     else
 580       gHistPN[3]->GetZaxis()->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
 581     //gHistPN[3]->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
 582     cPN[3] = new TCanvas("cPN3","",0,300,600,500);
 583     cPN[3]->SetFillColor(10);
 584     cPN[3]->SetHighLightColor(10);
 585     gHistPN[3]->DrawCopy("surf1fb");
 586     gPad->SetTheta(30); // default is 30
 587     //gPad->SetPhi(130); // default is 30
 588     gPad->SetPhi(-60); // default is 30
 589     gPad->Update();
 590     pngName = "CorrelationFunction.Centrality";
 591     pngName += centralityArray[gCentrality-1];
 592     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 593     pngName += ".PositiveNegative.png";
 594     //cPN[3]->SaveAs(pngName.Data());
 595   }
 596   // if no mixing then divide by convoluted histograms
 597   else if(listQA){
 598
 599     if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(+-) shuffled","(+-) convoluted");
 600     else                                histoTitle.ReplaceAll("(+-)","(+-) convoluted");
 601
 602     if(rebinEta > 1 || rebinPhi > 1){
 603       gHistPN[2]->Rebin2D(rebinEta,rebinPhi);
 604       gHistPN[0]->Scale(1./(Double_t)(rebinEta*rebinPhi));
 605     }
 606
 607     // normalization to 1 at (0,0) --> Jan Fietes method
 608     if(normToTrig){
 609       Double_t mixedNorm = gHistPN[2]->Integral(gHistPN[2]->GetXaxis()->FindBin(0-10e-5),gHistPN[2]->GetXaxis()->FindBin(0+10e-5),1,gHistPN[2]->GetNbinsX());
 610       mixedNorm /= gHistPN[2]->GetNbinsY()*(gHistPN[2]->GetXaxis()->FindBin(0.01) - gHistPN[2]->GetXaxis()->FindBin(-0.01) + 1);
 611       gHistPN[2]->Scale(1./mixedNorm);
 612     }
 613
 614     gHistPN[2]->GetYaxis()->SetTitleOffset(1.5);
 615     gHistPN[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 616     gHistPN[2]->SetTitle(histoTitle.Data());
 617     cPN[2] = new TCanvas("cPN2","",0,200,600,500);
 618     cPN[2]->SetFillColor(10);
 619     cPN[2]->SetHighLightColor(10);
 620     gHistPN[2]->DrawCopy("surf1fb");
 621     gPad->SetTheta(30); // default is 30
 622     //gPad->SetPhi(130); // default is 30
 623     gPad->SetPhi(-60); // default is 30
 624     gPad->Update();
 625     pngName = "DeltaPhiDeltaEtaMixed.Centrality";
 626     pngName += centralityArray[gCentrality-1];
 627     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 628     pngName += ".PositiveNegative.png";
 629     //cPN[2]->SaveAs(pngName.Data());
 630
 631     //Correlation function (+-)
 632     gHistPN[3] = dynamic_cast<TH2D *>(gHistPN[0]->Clone());
 633     gHistPN[3]->Divide(gHistPN[2]);
 634     gHistPN[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
 635     if(normToTrig)
 636       gHistPN[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
 637     else
 638       gHistPN[3]->GetZaxis()->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
 639     //gHistPN[3]->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
 640     cPN[3] = new TCanvas("cPN3","",0,300,600,500);
 641     cPN[3]->SetFillColor(10);
 642     cPN[3]->SetHighLightColor(10);
 643     gHistPN[3]->DrawCopy("surf1fb");
 644     gPad->SetTheta(30); // default is 30
 645     //gPad->SetPhi(130); // default is 30
 646     gPad->SetPhi(-60); // default is 30
 647     gPad->Update();
 648     pngName = "CorrelationFunction.Centrality";
 649     pngName += centralityArray[gCentrality-1];
 650     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 651     pngName += ".PositiveNegative.png";
 652     //cPN[3]->SaveAs(pngName.Data());
 653   }
 654
 655   //(-+)
 656   if(eventClass == "Centrality"){
 657     histoTitle = "(-+) | Centrality: ";
 658     histoTitle += psiMin;
 659     histoTitle += " - ";
 660     histoTitle += psiMax;
 661     histoTitle += " % ";
 662     histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
 663   }
 664   else if(eventClass == "Multiplicity"){
 665     histoTitle = "(-+) | Multiplicity: ";
 666     histoTitle += psiMin;
 667     histoTitle += " - ";
 668     histoTitle += psiMax;
 669     histoTitle += " tracks";
 670     histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
 671   }
 672   else{ // "EventPlane" (default)
 673     histoTitle = "(-+) | Centrality: ";
 674     histoTitle += centralityArray[gCentrality-1];
 675     histoTitle += "%";
 676     if((psiMin == -0.5)&&(psiMax == 0.5))
 677       histoTitle += " (-7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o})";
 678     else if((psiMin == 0.5)&&(psiMax == 1.5))
 679       histoTitle += " (37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o})";
 680     else if((psiMin == 1.5)&&(psiMax == 2.5))
 681       histoTitle += " (82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o})";
 682     else
 683       histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
 684   }
 685
 686   gHistNP[0] = b->GetCorrelationFunctionNP(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 687   if(rebinEta > 1 || rebinPhi > 1){
 688     gHistNP[0]->Rebin2D(rebinEta,rebinPhi);
 689     gHistNP[0]->Scale(1./(Double_t)(rebinEta*rebinPhi));
 690   }
 691   gHistNP[0]->GetYaxis()->SetTitleOffset(1.5);
 692   gHistNP[0]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 693   gHistNP[0]->SetTitle(histoTitle.Data());
 694   cNP[0] = new TCanvas("cNP0","",100,0,600,500);
 695   cNP[0]->SetFillColor(10);
 696   cNP[0]->SetHighLightColor(10);
 697   gHistNP[0]->DrawCopy("surf1fb");
 698   gPad->SetTheta(30); // default is 30
 699   //gPad->SetPhi(130); // default is 30
 700   gPad->SetPhi(-60); // default is 30
 701   gPad->Update();
 702   pngName = "DeltaPhiDeltaEta.Centrality";
 703   pngName += centralityArray[gCentrality-1];
 704   pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 705   pngName += ".NegativePositive.png";
 706   //cNP[0]->SaveAs(pngName.Data());
 707
 708   if(listBFShuffled) {
 709
 710     histoTitle.ReplaceAll("(-+)","(-+) shuffled");
 711
 712     gHistNP[1] = bShuffled->GetCorrelationFunctionNP(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 713     if(rebinEta > 1 || rebinPhi > 1){
 714       gHistNP[1]->Rebin2D(rebinEta,rebinPhi);
 715       gHistNP[1]->Scale(1./(Double_t)(rebinEta*rebinPhi));
 716     }
 717     gHistNP[1]->GetYaxis()->SetTitleOffset(1.5);
 718     gHistNP[1]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 719     gHistNP[1]->SetTitle(histoTitle.Data());
 720     cNP[1] = new TCanvas("cNP1","",100,100,600,500);
 721     cNP[1]->SetFillColor(10);
 722     cNP[1]->SetHighLightColor(10);
 723     gHistNP[1]->DrawCopy("surf1fb");
 724     gPad->SetTheta(30); // default is 30
 725     //gPad->SetPhi(130); // default is 30
 726     gPad->SetPhi(-60); // default is 30
 727     gPad->Update();
 728     pngName = "DeltaPhiDeltaEtaShuffled.Centrality";
 729     pngName += centralityArray[gCentrality-1];
 730     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 731     pngName += ".NegativePositive.png";
 732     //cNP[1]->SaveAs(pngName.Data());
 733   }
 734
 735   if(listBFMixed) {
 736
 737     if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(-+) shuffled","(-+) mixed");
 738     else                                histoTitle.ReplaceAll("(-+)","(-+) mixed");
 739
 740     // if normalization to trigger then do not divide Event mixing by number of trigger particles
 741     gHistNP[2] = bMixed->GetCorrelationFunctionNP(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 742     if(rebinEta > 1 || rebinPhi > 1){
 743       gHistNP[2]->Rebin2D(rebinEta,rebinPhi);
 744       gHistNP[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));
 745     }
 746     // normalization to 1 at (0,0) --> Jan Fietes method
 747     if(normToTrig){
 748       Double_t mixedNorm = gHistNP[2]->Integral(gHistNP[2]->GetXaxis()->FindBin(0-10e-5),gHistNP[2]->GetXaxis()->FindBin(0+10e-5),1,gHistNP[2]->GetNbinsX());
 749       mixedNorm /= gHistNP[2]->GetNbinsY()*(gHistNP[2]->GetXaxis()->FindBin(0.01) - gHistNP[2]->GetXaxis()->FindBin(-0.01) + 1);
 750       gHistNP[2]->Scale(1./mixedNorm);
 751     }
 752
 753     gHistNP[2]->GetYaxis()->SetTitleOffset(1.5);
 754     gHistNP[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 755     gHistNP[2]->SetTitle(histoTitle.Data());
 756     cNP[2] = new TCanvas("cNP2","",100,200,600,500);
 757     cNP[2]->SetFillColor(10);
 758     cNP[2]->SetHighLightColor(10);
 759     gHistNP[2]->DrawCopy("surf1fb");
 760     gPad->SetTheta(30); // default is 30
 761     //gPad->SetPhi(130); // default is 30
 762     gPad->SetPhi(-60); // default is 30
 763     gPad->Update();
 764     pngName = "DeltaPhiDeltaEtaMixed.Centrality";
 765     pngName += centralityArray[gCentrality-1];
 766     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 767     pngName += ".NegativePositive.png";
 768     //cNP[2]->SaveAs(pngName.Data());
 769
 770     //Correlation function (-+)
 771     gHistNP[3] = b->GetCorrelationFunction("NP",psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 772     gHistNP[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
 773     if(normToTrig)
 774       gHistNP[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
 775     else
 776       gHistNP[3]->GetZaxis()->SetTitle("C_{-+}(#Delta#eta,#Delta#varphi)");
 777     //gHistNP[3]->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
 778     cNP[3] = new TCanvas("cNP3","",100,300,600,500);
 779     cNP[3]->SetFillColor(10);
 780     cNP[3]->SetHighLightColor(10);
 781     gHistNP[3]->DrawCopy("surf1fb");
 782     gPad->SetTheta(30); // default is 30
 783     //gPad->SetPhi(130); // default is 30
 784     gPad->SetPhi(-60); // default is 30
 785     gPad->Update();
 786     pngName = "CorrelationFunction.Centrality";
 787     pngName += centralityArray[gCentrality-1];
 788     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 789     pngName += ".NegativePositive.png";
 790     //cNP[3]->SaveAs(pngName.Data());
 791   }
 792   // if no mixing then divide by convoluted histograms
 793   else if(listQA){
 794
 795     if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(-+) shuffled","(-+) convoluted");
 796     else                                histoTitle.ReplaceAll("(-+)","(-+) convoluted");
 797
 798     if(rebinEta > 1 || rebinPhi > 1){
 799       gHistNP[2]->Rebin2D(rebinEta,rebinPhi);
 800       gHistNP[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));
 801     }
 802
 803     // normalization to 1 at (0,0) --> Jan Fietes method
 804     if(normToTrig){
 805       Double_t mixedNorm = gHistNP[2]->Integral(gHistNP[2]->GetXaxis()->FindBin(0-10e-5),gHistNP[2]->GetXaxis()->FindBin(0+10e-5),1,gHistNP[2]->GetNbinsX());
 806       mixedNorm /= gHistNP[2]->GetNbinsY()*(gHistNP[2]->GetXaxis()->FindBin(0.01) - gHistNP[2]->GetXaxis()->FindBin(-0.01) + 1);
 807       gHistNP[2]->Scale(1./mixedNorm);
 808     }
 809
 810     gHistNP[2]->GetYaxis()->SetTitleOffset(1.5);
 811     gHistNP[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 812     gHistNP[2]->SetTitle(histoTitle.Data());
 813     cNP[2] = new TCanvas("cNP2","",100,200,600,500);
 814     cNP[2]->SetFillColor(10);
 815     cNP[2]->SetHighLightColor(10);
 816     gHistNP[2]->DrawCopy("surf1fb");
 817     gPad->SetTheta(30); // default is 30
 818     //gPad->SetPhi(130); // default is 30
 819     gPad->SetPhi(-60); // default is 30
 820     gPad->Update();
 821     pngName = "DeltaPhiDeltaEtaMixed.Centrality";
 822     pngName += centralityArray[gCentrality-1];
 823     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 824     pngName += ".NegativePositive.png";
 825     //cNP[2]->SaveAs(pngName.Data());
 826
 827     //Correlation function (-+)
 828     gHistNP[3] = dynamic_cast<TH2D *>(gHistNP[0]->Clone());
 829     gHistNP[3]->Divide(gHistNP[2]);
 830     gHistNP[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
 831     if(normToTrig)
 832       gHistNP[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
 833     else
 834       gHistNP[3]->GetZaxis()->SetTitle("C_{-+}(#Delta#eta,#Delta#varphi)");
 835     //gHistNP[3]->SetTitle("C_{-+}(#Delta#eta,#Delta#varphi)");
 836     cNP[3] = new TCanvas("cNP3","",100,300,600,500);
 837     cNP[3]->SetFillColor(10);
 838     cNP[3]->SetHighLightColor(10);
 839     gHistNP[3]->DrawCopy("surf1fb");
 840     gPad->SetTheta(30); // default is 30
 841     //gPad->SetPhi(130); // default is 30
 842     gPad->SetPhi(-60); // default is 30
 843     gPad->Update();
 844     pngName = "CorrelationFunction.Centrality";
 845     pngName += centralityArray[gCentrality-1];
 846     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 847     pngName += ".NegativePositive.png";
 848     //cNP[3]->SaveAs(pngName.Data());
 849   }
 850
 851
 852   //(++)
 853   if(eventClass == "Centrality"){
 854     histoTitle = "(++) | Centrality: ";
 855     histoTitle += psiMin;
 856     histoTitle += " - ";
 857     histoTitle += psiMax;
 858     histoTitle += " % ";
 859     histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
 860   }
 861   else if(eventClass == "Multiplicity"){
 862     histoTitle = "(++) | Multiplicity: ";
 863     histoTitle += psiMin;
 864     histoTitle += " - ";
 865     histoTitle += psiMax;
 866     histoTitle += " tracks";
 867     histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
 868   }
 869   else{ // "EventPlane" (default)
 870     histoTitle = "(++) | Centrality: ";
 871     histoTitle += centralityArray[gCentrality-1];
 872     histoTitle += "%";
 873     if((psiMin == -0.5)&&(psiMax == 0.5))
 874       histoTitle += " (-7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o})";
 875     else if((psiMin == 0.5)&&(psiMax == 1.5))
 876       histoTitle += " (37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o})";
 877     else if((psiMin == 1.5)&&(psiMax == 2.5))
 878       histoTitle += " (82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o})";
 879     else
 880       histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
 881   }
 882
 883   gHistPP[0] = b->GetCorrelationFunctionPP(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 884   if(rebinEta > 1 || rebinPhi > 1){
 885     gHistPP[0]->Rebin2D(rebinEta,rebinPhi);
 886     gHistPP[0]->Scale(1./(Double_t)(rebinEta*rebinPhi));
 887   }
 888   gHistPP[0]->GetYaxis()->SetTitleOffset(1.5);
 889   gHistPP[0]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 890   gHistPP[0]->SetTitle(histoTitle.Data());
 891   cPP[0] = new TCanvas("cPP0","",200,0,600,500);
 892   cPP[0]->SetFillColor(10);
 893   cPP[0]->SetHighLightColor(10);
 894   gHistPP[0]->DrawCopy("surf1fb");
 895   gPad->SetTheta(30); // default is 30
 896   //gPad->SetPhi(130); // default is 30
 897   gPad->SetPhi(-60); // default is 30
 898   gPad->Update();
 899   pngName = "DeltaPhiDeltaEta.Centrality";
 900   pngName += centralityArray[gCentrality-1];
 901   pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 902   pngName += ".PositivePositive.png";
 903   //cPP[0]->SaveAs(pngName.Data());
 904
 905   if(listBFShuffled) {
 906
 907     histoTitle.ReplaceAll("(++)","(++) shuffled");
 908
 909     gHistPP[1] = bShuffled->GetCorrelationFunctionPP(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 910     if(rebinEta > 1 || rebinPhi > 1){
 911       gHistPP[1]->Rebin2D(rebinEta,rebinPhi);
 912       gHistPP[1]->Scale(1./(Double_t)(rebinEta*rebinPhi));
 913     }
 914     gHistPP[1]->GetYaxis()->SetTitleOffset(1.5);
 915     gHistPP[1]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 916     gHistPP[1]->SetTitle(histoTitle.Data());
 917     cPP[1] = new TCanvas("cPP1","",200,100,600,500);
 918     cPP[1]->SetFillColor(10);
 919     cPP[1]->SetHighLightColor(10);
 920     gHistPP[1]->DrawCopy("surf1fb");
 921     gPad->SetTheta(30); // default is 30
 922     //gPad->SetPhi(130); // default is 30
 923     gPad->SetPhi(-60); // default is 30
 924     gPad->Update();
 925     pngName = "DeltaPhiDeltaEtaShuffled.Centrality";
 926     pngName += centralityArray[gCentrality-1];
 927     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 928     pngName += ".PositivePositive.png";
 929     //cPP[1]->SaveAs(pngName.Data());
 930   }
 931
 932   if(listBFMixed) {
 933
 934     if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(++) shuffled","(++) mixed");
 935     else                                histoTitle.ReplaceAll("(++)","(++) mixed");
 936
 937     // if normalization to trigger then do not divide Event mixing by number of trigger particles
 938     gHistPP[2] = bMixed->GetCorrelationFunctionPP(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 939     if(rebinEta > 1 || rebinPhi > 1){
 940       gHistPP[2]->Rebin2D(rebinEta,rebinPhi);
 941       gHistPP[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));
 942     }
 943     // normalization to 1 at (0,0) --> Jan Fietes method
 944     if(normToTrig){
 945       Double_t mixedNorm = gHistPP[2]->Integral(gHistPP[2]->GetXaxis()->FindBin(0-10e-5),gHistPP[2]->GetXaxis()->FindBin(0+10e-5),1,gHistPP[2]->GetNbinsX());
 946       mixedNorm /= gHistPP[2]->GetNbinsY()*(gHistPP[2]->GetXaxis()->FindBin(0.01) - gHistPP[2]->GetXaxis()->FindBin(-0.01) + 1);
 947       gHistPP[2]->Scale(1./mixedNorm);
 948     }
 949
 950     gHistPP[2]->GetYaxis()->SetTitleOffset(1.5);
 951     gHistPP[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 952     gHistPP[2]->SetTitle(histoTitle.Data());
 953     cPP[2] = new TCanvas("cPP2","",200,200,600,500);
 954     cPP[2]->SetFillColor(10);
 955     cPP[2]->SetHighLightColor(10);
 956     gHistPP[2]->DrawCopy("surf1fb");
 957     gPad->SetTheta(30); // default is 30
 958     //gPad->SetPhi(130); // default is 30
 959     gPad->SetPhi(-60); // default is 30
 960     gPad->Update();
 961     pngName = "DeltaPhiDeltaEtaMixed.Centrality";
 962     pngName += centralityArray[gCentrality-1];
 963     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 964     pngName += ".PositivePositive.png";
 965     //cPP[2]->SaveAs(pngName.Data());
 966
 967     //Correlation function (++)
 968     gHistPP[3] = b->GetCorrelationFunction("PP",psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 969     gHistPP[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
 970     if(normToTrig)
 971       gHistPP[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
 972     else
 973       gHistPP[3]->GetZaxis()->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
 974     // gHistPP[3]->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
 975     cPP[3] = new TCanvas("cPP3","",200,300,600,500);
 976     cPP[3]->SetFillColor(10);
 977     cPP[3]->SetHighLightColor(10);
 978     gHistPP[3]->DrawCopy("surf1fb");
 979     gPad->SetTheta(30); // default is 30
 980     //gPad->SetPhi(130); // default is 30
 981     gPad->SetPhi(-60); // default is 30
 982     gPad->Update();
 983     pngName = "CorrelationFunction.Centrality";
 984     pngName += centralityArray[gCentrality-1];
 985     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 986     pngName += ".PositivePositive.png";
 987     //cPP[3]->SaveAs(pngName.Data());
 988   }
 989   // if no mixing then divide by convoluted histograms
 990   else if(listQA){
 991
 992     if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(++) shuffled","(++) convoluted");
 993     else                                histoTitle.ReplaceAll("(++)","(++) convoluted");
 994
 995     if(rebinEta > 1 || rebinPhi > 1){
 996       gHistPP[2]->Rebin2D(rebinEta,rebinPhi);
 997       gHistPP[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));
 998     }
 999     // normalization to 1 at (0,0) --> Jan Fietes method
1000     if(normToTrig){
1001       Double_t mixedNorm = gHistPP[2]->Integral(gHistPP[2]->GetXaxis()->FindBin(0-10e-5),gHistPP[2]->GetXaxis()->FindBin(0+10e-5),1,gHistPP[2]->GetNbinsX());
1002       mixedNorm /= gHistPP[2]->GetNbinsY()*(gHistPP[2]->GetXaxis()->FindBin(0.01) - gHistPP[2]->GetXaxis()->FindBin(-0.01) + 1);
1003       gHistPP[2]->Scale(1./mixedNorm);
1004     }
1005
1006     gHistPP[2]->GetYaxis()->SetTitleOffset(1.5);
1007     gHistPP[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
1008     gHistPP[2]->SetTitle(histoTitle.Data());
1009     cPP[2] = new TCanvas("cPP2","",200,200,600,500);
1010     cPP[2]->SetFillColor(10);
1011     cPP[2]->SetHighLightColor(10);
1012     gHistPP[2]->DrawCopy("surf1fb");
1013     gPad->SetTheta(30); // default is 30
1014     //gPad->SetPhi(130); // default is 30
1015     gPad->SetPhi(-60); // default is 30
1016     gPad->Update();
1017     pngName = "DeltaPhiDeltaEtaMixed.Centrality";
1018     pngName += centralityArray[gCentrality-1];
1019     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
1020     pngName += ".PositivePositive.png";
1021     //cPP[2]->SaveAs(pngName.Data());
1022
1023     //Correlation function (++)
1024     gHistPP[3] = dynamic_cast<TH2D *>(gHistPP[0]->Clone());
1025     gHistPP[3]->Divide(gHistPP[2]);
1026     gHistPP[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
1027     if(normToTrig)
1028       gHistPP[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
1029     else
1030       gHistPP[3]->GetZaxis()->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
1031     //gHistPP[3]->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
1032     cPP[3] = new TCanvas("cPP3","",200,300,600,500);
1033     cPP[3]->SetFillColor(10);
1034     cPP[3]->SetHighLightColor(10);
1035     gHistPP[3]->DrawCopy("surf1fb");
1036     gPad->SetTheta(30); // default is 30
1037     //gPad->SetPhi(130); // default is 30
1038     gPad->SetPhi(-60); // default is 30
1039     gPad->Update();
1040     pngName = "CorrelationFunction.Centrality";
1041     pngName += centralityArray[gCentrality-1];
1042     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
1043     pngName += ".PositivePositive.png";
1044     //cPP[3]->SaveAs(pngName.Data());
1045   }
1046
1047   //(--)
1048   if(eventClass == "Centrality"){
1049     histoTitle = "(--) | Centrality: ";
1050     histoTitle += psiMin;
1051     histoTitle += " - ";
1052     histoTitle += psiMax;
1053     histoTitle += " % ";
1054     histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
1055   }
1056   else if(eventClass == "Multiplicity"){
1057     histoTitle = "(--) | Multiplicity: ";
1058     histoTitle += psiMin;
1059     histoTitle += " - ";
1060     histoTitle += psiMax;
1061     histoTitle += " tracks";
1062     histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
1063   }
1064   else{ // "EventPlane" (default)
1065     histoTitle = "(--) | Centrality: ";
1066     histoTitle += centralityArray[gCentrality-1];
1067     histoTitle += "%";
1068     if((psiMin == -0.5)&&(psiMax == 0.5))
1069       histoTitle += " (-7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o})";
1070     else if((psiMin == 0.5)&&(psiMax == 1.5))
1071       histoTitle += " (37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o})";
1072     else if((psiMin == 1.5)&&(psiMax == 2.5))
1073       histoTitle += " (82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o})";
1074     else
1075       histoTitle += " (0^{o} < #varphi^{t} - #Psi_{2} < 180^{o})";
1076   }
1077
1078   gHistNN[0] = b->GetCorrelationFunctionNN(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1079   if(rebinEta > 1 || rebinPhi > 1){
1080     gHistNN[0]->Rebin2D(rebinEta,rebinPhi);
1081     gHistNN[0]->Scale(1./(Double_t)(rebinEta*rebinPhi));
1082   }
1083   gHistNN[0]->GetYaxis()->SetTitleOffset(1.5);
1084   gHistNN[0]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
1085   gHistNN[0]->SetTitle(histoTitle.Data());
1086   cNN[0] = new TCanvas("cNN0","",300,0,600,500);
1087   cNN[0]->SetFillColor(10);
1088   cNN[0]->SetHighLightColor(10);
1089   gHistNN[0]->DrawCopy("surf1fb");
1090   gPad->SetTheta(30); // default is 30
1091   gPad->SetPhi(-60); // default is 30
1092   //gPad->SetPhi(-60); // default is 30
1093   gPad->Update();
1094   pngName = "DeltaPhiDeltaEta.Centrality";
1095   pngName += centralityArray[gCentrality-1];
1096   pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
1097   pngName += ".NegativeNegative.png";
1098   //cNN[0]->SaveAs(pngName.Data());
1099
1100   if(listBFShuffled) {
1101
1102     histoTitle.ReplaceAll("(--)","(--) shuffled");
1103
1104     gHistNN[1] = bShuffled->GetCorrelationFunctionNN(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1105     if(rebinEta > 1 || rebinPhi > 1){
1106       gHistNN[1]->Rebin2D(rebinEta,rebinPhi);
1107       gHistNN[1]->Scale(1./(Double_t)(rebinEta*rebinPhi));
1108     }
1109     gHistNN[1]->GetYaxis()->SetTitleOffset(1.5);
1110     gHistNN[1]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
1111     gHistNN[1]->SetTitle(histoTitle.Data());
1112     cNN[1] = new TCanvas("cNN1","",300,100,600,500);
1113     cNN[1]->SetFillColor(10);
1114     cNN[1]->SetHighLightColor(10);
1115     gHistNN[1]->DrawCopy("surf1fb");
1116     gPad->SetTheta(30); // default is 30
1117     //gPad->SetPhi(130); // default is 30
1118     gPad->SetPhi(-60); // default is 30
1119     gPad->Update();
1120     pngName = "DeltaPhiDeltaEtaShuffled.Centrality";
1121     pngName += centralityArray[gCentrality-1];
1122     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
1123     pngName += ".NegativeNegative.png";
1124     //cNN[1]->SaveAs(pngName.Data());
1125   }
1126
1127   if(listBFMixed) {
1128
1129     if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(--) shuffled","(--) mixed");
1130     else                                histoTitle.ReplaceAll("(--)","(--) mixed");
1131
1132     // if normalization to trigger then do not divide Event mixing by number of trigger particles
1133     gHistNN[2] = bMixed->GetCorrelationFunctionNN(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
1134     if(rebinEta > 1 || rebinPhi > 1){
1135       gHistNN[2]->Rebin2D(rebinEta,rebinPhi);
1136       gHistNN[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));
1137     }
1138     // normalization to 1 at (0,0) --> Jan Fietes method
1139     if(normToTrig){
1140       Double_t mixedNorm = gHistNN[2]->Integral(gHistNN[2]->GetXaxis()->FindBin(0-10e-5),gHistNN[2]->GetXaxis()->FindBin(0+10e-5),1,gHistNN[2]->GetNbinsX());
1141       mixedNorm /= gHistNN[2]->GetNbinsY()*(gHistNN[2]->GetXaxis()->FindBin(0.01) - gHistNN[2]->GetXaxis()->FindBin(-0.01) + 1);
1142       gHistNN[2]->Scale(1./mixedNorm);
1143     }
1144
1145     gHistNN[2]->GetYaxis()->SetTitleOffset(1.5);
1146     gHistNN[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
1147     gHistNN[2]->SetTitle(histoTitle.Data());
1148     cNN[2] = new TCanvas("cNN2","",300,200,600,500);
1149     cNN[2]->SetFillColor(10);
1150     cNN[2]->SetHighLightColor(10);
1151     gHistNN[2]->DrawCopy("surf1fb");
1152     gPad->SetTheta(30); // default is 30
1153     //gPad->SetPhi(130); // default is 30
1154     gPad->SetPhi(-60); // default is 30
1155     gPad->Update();
1156     pngName = "DeltaPhiDeltaEtaMixed.Centrality";
1157     pngName += centralityArray[gCentrality-1];
1158     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
1159     pngName += ".NegativeNegative.png";
1160     //cNN[2]->SaveAs(pngName.Data());
1161
1162     //Correlation function (--)
1163     gHistNN[3] = b->GetCorrelationFunction("NN",psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
1164     gHistNN[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
1165     if(normToTrig)
1166       gHistNN[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
1167     else
1168       gHistNN[3]->GetZaxis()->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
1169     // gHistNN[3]->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
1170     cNN[3] = new TCanvas("cNN3","",300,300,600,500);
1171     cNN[3]->SetFillColor(10);
1172     cNN[3]->SetHighLightColor(10);
1173     gHistNN[3]->DrawCopy("surf1fb");
1174     gPad->SetTheta(30); // default is 30
1175     //gPad->SetPhi(130); // default is 30
1176     gPad->SetPhi(-60); // default is 30
1177     gPad->Update();
1178     pngName = "CorrelationFunction.Centrality";
1179     pngName += centralityArray[gCentrality-1];
1180     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
1181     pngName += ".NegativeNegative.png";
1182     //cNN[3]->SaveAs(pngName.Data());
1183   }
1184   // if no mixing then divide by convoluted histograms
1185   else if(listQA){
1186
1187     if(histoTitle.Contains("shuffled")) histoTitle.ReplaceAll("(--) shuffled","(--) convoluted");
1188     else                                histoTitle.ReplaceAll("(--)","(--) convoluted");
1189
1190     if(rebinEta > 1 || rebinPhi > 1){
1191       gHistNN[2]->Rebin2D(rebinEta,rebinPhi);
1192       gHistNN[2]->Scale(1./(Double_t)(rebinEta*rebinPhi));
1193     }
1194     // normalization to 1 at (0,0) --> Jan Fietes method
1195     if(normToTrig){
1196       Double_t mixedNorm = gHistNN[2]->Integral(gHistNN[2]->GetXaxis()->FindBin(0-10e-5),gHistNN[2]->GetXaxis()->FindBin(0+10e-5),1,gHistNN[2]->GetNbinsX());
1197       mixedNorm /= gHistNN[2]->GetNbinsY()*(gHistNN[2]->GetXaxis()->FindBin(0.01) - gHistNN[2]->GetXaxis()->FindBin(-0.01) + 1);
1198       gHistNN[2]->Scale(1./mixedNorm);
1199     }
1200
1201     gHistNN[2]->GetYaxis()->SetTitleOffset(1.5);
1202     gHistNN[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
1203     gHistNN[2]->SetTitle(histoTitle.Data());
1204     cNN[2] = new TCanvas("cNN2","",300,200,600,500);
1205     cNN[2]->SetFillColor(10);
1206     cNN[2]->SetHighLightColor(10);
1207     gHistNN[2]->DrawCopy("surf1fb");
1208     gPad->SetTheta(30); // default is 30
1209     //gPad->SetPhi(130); // default is 30
1210     gPad->SetPhi(-60); // default is 30
1211     gPad->Update();
1212     pngName = "DeltaPhiDeltaEtaMixed.Centrality";
1213     pngName += centralityArray[gCentrality-1];
1214     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
1215     pngName += ".NegativeNegative.png";
1216     //cNN[2]->SaveAs(pngName.Data());
1217
1218     //Correlation function (--)
1219     gHistNN[3] = dynamic_cast<TH2D *>(gHistNN[0]->Clone());
1220     gHistNN[3]->Divide(gHistNN[2]);
1221     gHistNN[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
1222     if(normToTrig)
1223       gHistNN[3]->GetZaxis()->SetTitle("#frac{1}{N_{trig}}#frac{d^{2}N_{assoc}}{d#Delta#eta#Delta#varphi} (rad^{-1})");
1224     else
1225       gHistNN[3]->GetZaxis()->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
1226     //gHistNN[3]->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
1227     cNN[3] = new TCanvas("cNN3","",300,300,600,500);
1228     cNN[3]->SetFillColor(10);
1229     cNN[3]->SetHighLightColor(10);
1230     gHistNN[3]->DrawCopy("surf1fb");
1231     gPad->SetTheta(30); // default is 30
1232     //gPad->SetPhi(130); // default is 30
1233     gPad->SetPhi(-60); // default is 30
1234     gPad->Update();
1235     pngName = "CorrelationFunction.Centrality";
1236     pngName += centralityArray[gCentrality-1];
1237     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
1238     pngName += ".NegativeNegative.png";
1239     //cNN[3]->SaveAs(pngName.Data());
1240   }
1241
1242   //Write to output file
1243   TString newFileName = "correlationFunction.";
1244   if(eventClass == "Centrality"){
1245     newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1246     newFileName += ".PsiAll.PttFrom";
1247   }
1248   else if(eventClass == "Multiplicity"){
1249     newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1250     newFileName += ".PsiAll.PttFrom";
1251   }
1252   else{ // "EventPlane" (default)
1253     newFileName += "Centrality";
1254     newFileName += gCentrality; newFileName += ".Psi";
1255     if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
1256     else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
1257     else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
1258     else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
1259     else newFileName += "All.PttFrom";
1260   }
1261   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
1262   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
1263   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
1264   newFileName += Form("%.1f",ptAssociatedMax);
1265
1266   newFileName += "_";
1267   newFileName += Form("%.1f",psiMin);
1268   newFileName += "-";
1269   newFileName += Form("%.1f",psiMax);
1270   newFileName += ".root";
1271
1272   TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
1273   gHistPN[0]->SetName("gHistPNRaw"); gHistPN[0]->Write();
1274   gHistNP[0]->SetName("gHistNPRaw"); gHistNP[0]->Write();
1275   gHistPP[0]->SetName("gHistPPRaw"); gHistPP[0]->Write();
1276   gHistNN[0]->SetName("gHistNNRaw"); gHistNN[0]->Write();
1277   if(listBFShuffled) {
1278     gHistPN[1]->SetName("gHistPNShuffled"); gHistPN[1]->Write();
1279     gHistNP[1]->SetName("gHistNPShuffled"); gHistNP[1]->Write();
1280     gHistPP[1]->SetName("gHistPPShuffled"); gHistPP[1]->Write();
1281     gHistNN[1]->SetName("gHistNNShuffled"); gHistNN[1]->Write();
1282   }
1283   if(listBFMixed || (!listBFMixed&&listQA)) {
1284     gHistPN[2]->SetName("gHistPNMixed"); gHistPN[2]->Write();
1285     gHistNP[2]->SetName("gHistNPMixed"); gHistNP[2]->Write();
1286     gHistPP[2]->SetName("gHistPPMixed"); gHistPP[2]->Write();
1287     gHistNN[2]->SetName("gHistNNMixed"); gHistNN[2]->Write();
1288
1289     gHistPN[3]->SetName("gHistPNCorrelationFunctions"); gHistPN[3]->Write();
1290     gHistNP[3]->SetName("gHistNPCorrelationFunctions"); gHistNP[3]->Write();
1291     gHistPP[3]->SetName("gHistPPCorrelationFunctions"); gHistPP[3]->Write();
1292     gHistNN[3]->SetName("gHistNNCorrelationFunctions"); gHistNN[3]->Write();
1293   }
1294   newFile->Close();
1295
1296   // some cleaning
1297   for(Int_t i = 0; i < 4; i++){
1298
1299     if(!listBFShuffled && i == 1) continue;
1300     if(!listBFMixed && (i == 2 || i == 3)) continue;
1301
1302     if(gHistPP[i]) delete gHistPP[i];
1303     if(gHistPN[i]) delete gHistPN[i];
1304     if(gHistNP[i]) delete gHistNP[i];
1305     if(gHistNN[i]) delete gHistNN[i];
1306
1307     if(cPN[i]) delete cPN[i];
1308     if(cNP[i]) delete cNP[i];
1309     if(cPP[i]) delete cPP[i];
1310     if(cNN[i]) delete cNN[i];
1311   }
1312
1313   delete hP;
1314   delete hN;
1315   delete hPP;
1316   delete hPN;
1317   delete hNP;
1318   delete hNN;
1319
1320   delete hPMixed;
1321   delete hNMixed;
1322   delete hPPMixed;
1323   delete hPNMixed;
1324   delete hNPMixed;
1325   delete hNNMixed;
1326
1327   delete hPShuffled;
1328   delete hNShuffled;
1329   delete hPPShuffled;
1330   delete hPNShuffled;
1331   delete hNPShuffled;
1332   delete hNNShuffled;
1333
1334 }
1335
1336 //____________________________________________________________//
1337 void drawCorrelationFunctions(const char* lhcPeriod = "LHC10h",
1338                               const char* gCentralityEstimator = "V0M",
1339                               Int_t gBit = 128,
1340                               const char* gEventPlaneEstimator = "VZERO",
1341                               Int_t gCentrality = 1,
1342                               Double_t psiMin = -0.5, Double_t psiMax = 3.5,
1343                               Double_t vertexZMin = -10.,
1344                               Double_t vertexZMax = 10.,
1345                               Double_t ptTriggerMin = -1.,
1346                               Double_t ptTriggerMax = -1.,
1347                               Double_t ptAssociatedMin = -1.,
1348                               Double_t ptAssociatedMax = -1.,
1349                               Bool_t kFit = kFALSE) {
1350   //Macro that draws the charge dependent correlation functions
1351   //for each centrality bin for the different pT of trigger and
1352   //associated particles
1353   //Author: Panos.Christakoglou@nikhef.nl
1354   TGaxis::SetMaxDigits(3);
1355
1356   //Get the input file
1357   /* TString filename = lhcPeriod;
1358   filename += "/Centrality"; filename += gCentralityEstimator;
1359   filename += "_Bit"; filename += gBit;
1360   filename += "_"; filename += gEventPlaneEstimator;
1361   filename +="/PttFrom";
1362   filename += Form("%.1f",ptTriggerMin); filename += "To";
1363   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1364   filename += Form("%.1f",ptAssociatedMin); filename += "To";
1365   filename += Form("%.1f",ptAssociatedMax); */
1366
1367   TString filename = "correlationFunction.Centrality";
1368   filename += gCentrality; filename += ".Psi";
1369   if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
1370   else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
1371   else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
1372   else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
1373   else filename += "All.PttFrom";
1374   filename += Form("%.1f",ptTriggerMin); filename += "To";
1375   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1376   filename += Form("%.1f",ptAssociatedMin); filename += "To";
1377   filename += Form("%.1f",ptAssociatedMax);
1378
1379   filename += "_";
1380   filename += Form("%.1f",psiMin);
1381   filename += "-";
1382   filename += Form("%.1f",psiMax);
1383   filename += ".root";
1384
1385   //Open the file
1386   TFile *f = TFile::Open(filename.Data());
1387   if((!f)||(!f->IsOpen())) {
1388     Printf("The file %s is not found. Aborting...",filename);
1389     return;
1390   }
1391   //f->ls();
1392
1393   //Latex
1394   TString centralityLatex = "Centrality: ";
1395   centralityLatex += centralityArray[gCentrality-1];
1396   centralityLatex += "%";
1397
1398   TString psiLatex;
1399   if((psiMin == -0.5)&&(psiMax == 0.5))
1400     psiLatex = " -7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o}";
1401   else if((psiMin == 0.5)&&(psiMax == 1.5))
1402     psiLatex = " 37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o}";
1403   else if((psiMin == 1.5)&&(psiMax == 2.5))
1404     psiLatex = " 82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o}";
1405   else
1406     psiLatex = " 0^{o} < #varphi^{t} - #Psi_{2} < 180^{o}";
1407
1408   TString pttLatex = Form("%.1f",ptTriggerMin);
1409   pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
1410   pttLatex += " GeV/c";
1411
1412   TString ptaLatex = Form("%.1f",ptAssociatedMin);
1413   ptaLatex += " < p_{T,assoc} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1414   ptaLatex += " GeV/c";
1415
1416   TLatex *latexInfo1 = new TLatex();
1417   latexInfo1->SetNDC();
1418   latexInfo1->SetTextSize(0.045);
1419   latexInfo1->SetTextColor(1);
1420
1421   TString pngName;
1422
1423   //============================================================//
1424   //Get the +- correlation function
1425   TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
1426   gHistPN->SetStats(kFALSE);
1427   gHistPN->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
1428   gHistPN->GetXaxis()->SetRangeUser(-1.4,1.4);
1429   gHistPN->GetXaxis()->CenterTitle();
1430   gHistPN->GetXaxis()->SetTitleOffset(1.2);
1431   gHistPN->GetYaxis()->CenterTitle();
1432   gHistPN->GetYaxis()->SetTitleOffset(1.2);
1433   gHistPN->GetZaxis()->SetTitleOffset(1.5);
1434   TCanvas *cPN = new TCanvas("cPN","",0,0,600,500);
1435   cPN->SetFillColor(10); cPN->SetHighLightColor(10);
1436   cPN->SetLeftMargin(0.15);
1437   gHistPN->DrawCopy("surf1fb");
1438
1439   gPad->SetTheta(30); // default is 30
1440   gPad->SetPhi(-60); // default is 30
1441   gPad->Update();
1442
1443   latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
1444   latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
1445   latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
1446   latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
1447
1448   pngName = "CorrelationFunction.Centrality";
1449   pngName += centralityArray[gCentrality-1];
1450   pngName += ".Psi";
1451   if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
1452   else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
1453   else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
1454   else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
1455   else pngName += "All.PttFrom";
1456   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1457   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1458   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1459   pngName += Form("%.1f",ptAssociatedMax);
1460   pngName += ".PositiveNegative.png";
1461   cPN->SaveAs(pngName.Data());
1462   if(kFit)
1463     fitCorrelationFunctions(gCentrality, psiMin, psiMax,vertexZMin, vertexZMax,
1464                             ptTriggerMin,ptTriggerMax,
1465                             ptAssociatedMin, ptAssociatedMax,gHistPN);
1466
1467   //============================================================//
1468   //Get the -+ correlation function
1469   TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
1470   gHistNP->SetStats(kFALSE);
1471   gHistNP->SetTitle("C_{-+}(#Delta#eta,#Delta#varphi)");
1472   gHistNP->GetXaxis()->SetRangeUser(-1.4,1);
1473   gHistNP->GetXaxis()->CenterTitle();
1474   gHistNP->GetXaxis()->SetTitleOffset(1.2);
1475   gHistNP->GetYaxis()->CenterTitle();
1476   gHistNP->GetYaxis()->SetTitleOffset(1.2);
1477   gHistNP->GetZaxis()->SetTitleOffset(1.5);
1478   TCanvas *cNP = new TCanvas("cNP","",50,50,600,500);
1479   cNP->SetFillColor(10); cNP->SetHighLightColor(10);
1480   cNP->SetLeftMargin(0.15);
1481   gHistNP->DrawCopy("surf1fb");
1482
1483   gPad->SetTheta(30); // default is 30
1484   gPad->SetPhi(-60); // default is 30
1485   gPad->Update();
1486
1487   latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
1488   latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
1489   latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
1490   latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
1491
1492   pngName = "CorrelationFunction.Centrality";
1493   pngName += centralityArray[gCentrality-1];
1494   pngName += ".Psi";
1495   if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
1496   else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
1497   else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
1498   else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
1499   else pngName += "All.PttFrom";
1500   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1501   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1502   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1503   pngName += Form("%.1f",ptAssociatedMax);
1504   pngName += ".NegativePositive.png";
1505   cNP->SaveAs(pngName.Data());
1506
1507   if(kFit)
1508     fitCorrelationFunctions(gCentrality, psiMin, psiMax,vertexZMin, vertexZMax,
1509                             ptTriggerMin,ptTriggerMax,
1510                             ptAssociatedMin, ptAssociatedMax,gHistNP);
1511
1512   //============================================================//
1513   //Get the ++ correlation function
1514   TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
1515   gHistPP->SetStats(kFALSE);
1516   gHistPP->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
1517   gHistPP->GetXaxis()->SetRangeUser(-1.4,1.4);
1518   gHistPP->GetXaxis()->CenterTitle();
1519   gHistPP->GetXaxis()->SetTitleOffset(1.2);
1520   gHistPP->GetYaxis()->CenterTitle();
1521   gHistPP->GetYaxis()->SetTitleOffset(1.2);
1522   gHistPP->GetZaxis()->SetTitleOffset(1.5);
1523   TCanvas *cPP = new TCanvas("cPP","",100,100,600,500);
1524   cPP->SetFillColor(10); cPP->SetHighLightColor(10);
1525   cPP->SetLeftMargin(0.15);
1526   gHistPP->DrawCopy("surf1fb");
1527
1528   gPad->SetTheta(30); // default is 30
1529   gPad->SetPhi(-60); // default is 30
1530   gPad->Update();
1531
1532   latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
1533   latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
1534   latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
1535   latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
1536
1537   pngName = "CorrelationFunction.Centrality";
1538   pngName += centralityArray[gCentrality-1];
1539   pngName += ".Psi";
1540   if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
1541   else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
1542   else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
1543   else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
1544   else pngName += "All.PttFrom";
1545   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1546   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1547   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1548   pngName += Form("%.1f",ptAssociatedMax);
1549   pngName += ".PositivePositive.png";
1550   cPP->SaveAs(pngName.Data());
1551
1552   if(kFit)
1553     fitCorrelationFunctions(gCentrality, psiMin, psiMax,vertexZMin, vertexZMax,
1554                             ptTriggerMin,ptTriggerMax,
1555                             ptAssociatedMin, ptAssociatedMax,gHistPP);
1556
1557   //============================================================//
1558   //Get the -- correlation function
1559   TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
1560   gHistNN->SetStats(kFALSE);
1561   gHistNN->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
1562   gHistNN->GetXaxis()->SetRangeUser(-1.4,1.4);
1563   gHistNN->GetXaxis()->CenterTitle();
1564   gHistNN->GetXaxis()->SetTitleOffset(1.2);
1565   gHistNN->GetYaxis()->CenterTitle();
1566   gHistNN->GetYaxis()->SetTitleOffset(1.2);
1567   gHistNN->GetZaxis()->SetTitleOffset(1.5);
1568   TCanvas *cNN = new TCanvas("cNN","",150,150,600,500);
1569   cNN->SetFillColor(10); cNN->SetHighLightColor(10);
1570   cNN->SetLeftMargin(0.15);
1571   gHistNN->DrawCopy("surf1fb");
1572
1573   gPad->SetTheta(30); // default is 30
1574   gPad->SetPhi(-60); // default is 30
1575   gPad->Update();
1576
1577   latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
1578   latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
1579   latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
1580   latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
1581
1582   pngName = "CorrelationFunction.Centrality";
1583   pngName += centralityArray[gCentrality-1];
1584   pngName += ".Psi";
1585   if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
1586   else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
1587   else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
1588   else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
1589   else pngName += "All.PttFrom";
1590   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1591   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1592   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1593   pngName += Form("%.1f",ptAssociatedMax);
1594   pngName += ".NegativeNegative.png";
1595   cNN->SaveAs(pngName.Data());
1596
1597   if(kFit)
1598     fitCorrelationFunctions(gCentrality, psiMin, psiMax,vertexZMin, vertexZMax,
1599                             ptTriggerMin,ptTriggerMax,
1600                             ptAssociatedMin, ptAssociatedMax,gHistNN);
1601 }
1602
1603 //____________________________________________________________//
1604 void drawProjections(Bool_t kProjectInEta = kFALSE,
1605                      Int_t binMin = 1,
1606                      Int_t binMax = 80,
1607                      Int_t gCentrality = 1,
1608                      Double_t psiMin = -0.5,
1609                      Double_t psiMax = 3.5,
1610                      Double_t ptTriggerMin = -1.,
1611                      Double_t ptTriggerMax = -1.,
1612                      Double_t ptAssociatedMin = -1.,
1613                      Double_t ptAssociatedMax = -1.,
1614                      Bool_t kUseZYAM = kFALSE,
1615                      TString eventClass = "Centrality") {
1616   //Macro that draws the charge dependent correlation functions PROJECTIONS
1617   //for each centrality bin for the different pT of trigger and
1618   //associated particles
1619   TGaxis::SetMaxDigits(3);
1620
1621   //Get the input file
1622   /*TString filename = lhcPeriod;
1623   filename += "/Centrality"; filename += gCentralityEstimator;
1624   filename += "_Bit"; filename += gBit;
1625   filename += "_"; filename += gEventPlaneEstimator;
1626   filename +="/PttFrom";
1627   filename += Form("%.1f",ptTriggerMin); filename += "To";
1628   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1629   filename += Form("%.1f",ptAssociatedMin); filename += "To";
1630   filename += Form("%.1f",ptAssociatedMax); */
1631
1632   TString filename = "correlationFunction.";
1633   if(eventClass == "Centrality"){
1634     filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1635     filename += ".PsiAll.PttFrom";
1636   }
1637   else if(eventClass == "Multiplicity"){
1638     filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1639     filename += ".PsiAll.PttFrom";
1640   }
1641   else{ // "EventPlane" (default)
1642     filename += "Centrality";
1643     filename += gCentrality; filename += ".Psi";
1644     if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
1645     else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
1646     else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
1647     else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
1648     else filename += "All.PttFrom";
1649   }
1650   filename += Form("%.1f",ptTriggerMin); filename += "To";
1651   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1652   filename += Form("%.1f",ptAssociatedMin); filename += "To";
1653   filename += Form("%.1f",ptAssociatedMax);
1654   //if(k2pMethod) filename += "_2pMethod";
1655
1656   filename += "_";
1657   filename += Form("%.1f",psiMin);
1658   filename += "-";
1659   filename += Form("%.1f",psiMax);
1660
1661   filename += ".root";
1662
1663   //Open the file
1664   TFile *f = TFile::Open(filename.Data());
1665   if((!f)||(!f->IsOpen())) {
1666     Printf("The file %s is not found. Aborting...",filename);
1667     return listBF;
1668   }
1669   //f->ls();
1670
1671   //Latex
1672   TString centralityLatex = "Centrality: ";
1673   if(eventClass == "Centrality")
1674     centralityLatex += Form("%.0f - %.0f",psiMin,psiMax);
1675   else
1676     centralityLatex += centralityArray[gCentrality-1];
1677   centralityLatex += " %";
1678
1679   TString psiLatex;
1680   if((psiMin == -0.5)&&(psiMax == 0.5))
1681     psiLatex = " -7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o}";
1682   else if((psiMin == 0.5)&&(psiMax == 1.5))
1683     psiLatex = " 37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o}";
1684   else if((psiMin == 1.5)&&(psiMax == 2.5))
1685     psiLatex = " 82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o}";
1686   else
1687     psiLatex = " 0^{o} < #varphi^{t} - #Psi_{2} < 180^{o}";
1688
1689   TString pttLatex = Form("%.1f",ptTriggerMin);
1690   pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
1691   pttLatex += " GeV/c";
1692
1693   TString ptaLatex = Form("%.1f",ptAssociatedMin);
1694   ptaLatex += " < p_{T,assoc} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1695   ptaLatex += " GeV/c";
1696
1697   TLatex *latexInfo1 = new TLatex();
1698   latexInfo1->SetNDC();
1699   latexInfo1->SetTextSize(0.045);
1700   latexInfo1->SetTextColor(1);
1701
1702   TString pngName;
1703
1704   //============================================================//
1705   //Get the +- correlation function
1706   TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
1707   gHistPN->SetStats(kFALSE);
1708   gHistPN->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
1709   gHistPN->GetXaxis()->SetRangeUser(-1.4,1.4);
1710   gHistPN->GetXaxis()->CenterTitle();
1711   gHistPN->GetXaxis()->SetTitleOffset(1.2);
1712   gHistPN->GetYaxis()->CenterTitle();
1713   gHistPN->GetYaxis()->SetTitleOffset(1.2);
1714   gHistPN->GetZaxis()->SetTitleOffset(1.5);
1715   TCanvas *cPN = new TCanvas("cPN","",0,0,600,500);
1716   cPN->SetFillColor(10);
1717   cPN->SetHighLightColor(10);
1718   cPN->SetLeftMargin(0.15);
1719
1720   //================//
1721   TH1D* gHistPNprojection = 0x0;
1722   Double_t sum = 0.0;
1723   Double_t gError = 0.0;
1724   Int_t nCounter = 0;
1725   //projection in delta eta
1726   if(kProjectInEta) {
1727     gHistPNprojection = new TH1D("gHistPNprojection","",gHistPN->GetNbinsX(),gHistPN->GetXaxis()->GetXmin(),gHistPN->GetXaxis()->GetXmax());
1728     for(Int_t iBinX = 1; iBinX <= gHistPN->GetNbinsX(); iBinX++) {
1729       sum = 0.; gError = 0.0; nCounter = 0;
1730       for(Int_t iBinY = binMin; iBinY <= binMax; iBinY++) {
1731         sum += gHistPN->GetBinContent(iBinX,iBinY);
1732         if(gHistPN->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
1733         Double_t exy = gHistPN->GetCellError(iBinX,iBinY);
1734         gError += exy*exy;
1735       }
1736       if(nCounter != 0) {
1737         sum /= nCounter;
1738         gError = TMath::Sqrt(gError)/nCounter;
1739       }
1740       gHistPNprojection->SetBinContent(iBinX,sum);
1741       gHistPNprojection->SetBinError(iBinX,gError);
1742     }
1743     gHistPNprojection->GetXaxis()->SetRangeUser(-1.4,1.4);
1744     if(kUseZYAM)
1745       gHistPNprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#eta} - b_{ZYAM}");
1746     else
1747       gHistPNprojection->GetYaxis()->SetTitle("C_{+-}(#Delta#eta)");
1748     gHistPNprojection->GetXaxis()->SetTitle("#Delta#eta");
1749   }//projection in delta eta
1750   //projection in delta phi
1751   else {
1752     gHistPNprojection = new TH1D("gHistPNprojection","",gHistPN->GetNbinsY(),gHistPN->GetYaxis()->GetXmin(),gHistPN->GetYaxis()->GetXmax());
1753     for(Int_t iBinY = 1; iBinY <= gHistPN->GetNbinsY(); iBinY++) {
1754       sum = 0.; gError = 0.0; nCounter = 0;
1755       for(Int_t iBinX = binMin; iBinX <= binMax; iBinX++) {
1756       //for(Int_t iBinX = 1; iBinX <= gHistPN->GetNbinsX(); iBinX++) {
1757         sum += gHistPN->GetBinContent(iBinX,iBinY);
1758         if(gHistPN->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
1759         Double_t exy = gHistPN->GetCellError(iBinX,iBinY);
1760         gError += exy*exy;
1761       }
1762       if(nCounter != 0) {
1763         sum /= nCounter;
1764         gError = TMath::Sqrt(gError)/nCounter;
1765       }
1766       gHistPNprojection->SetBinContent(iBinY,sum);
1767       gHistPNprojection->SetBinError(iBinY,gError);
1768     }
1769     if(kUseZYAM)
1770       gHistPNprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#varphi} - b_{ZYAM} (rad^{-1})");
1771     else
1772       gHistPNprojection->GetYaxis()->SetTitle("C_{+-}(#Delta#varphi)");
1773     gHistPNprojection->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1774   }
1775
1776   //ZYAM
1777   if(kUseZYAM) {
1778     Double_t reference = gHistPNprojection->GetBinContent(gHistPNprojection->GetMinimumBin());
1779     for(Int_t iBinX = 1; iBinX <= gHistPNprojection->GetNbinsX(); iBinX++)
1780       gHistPNprojection->SetBinContent(iBinX,gHistPNprojection->GetBinContent(iBinX) - reference);
1781   }
1782
1783   gHistPNprojection->GetYaxis()->SetTitleOffset(1.5);
1784   gHistPNprojection->SetMarkerStyle(20);
1785   gHistPNprojection->SetStats(kFALSE);
1786   gHistPNprojection->DrawCopy("E");
1787   //=================//
1788
1789   gPad->SetTheta(30); // default is 30
1790   gPad->SetPhi(-60); // default is 30
1791   gPad->Update();
1792
1793   latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
1794   latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
1795   latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
1796   latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
1797
1798   pngName = filename;
1799   if(kProjectInEta) pngName.ReplaceAll(".root","_DeltaEtaProjection");
1800   else              pngName.ReplaceAll(".root","_DeltaPhiProjection");
1801   pngName += ".PositiveNegative.png";
1802   cPN->SaveAs(pngName.Data());
1803
1804   //============================================================//
1805   //Get the -+ correlation function
1806   TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
1807   gHistNP->SetStats(kFALSE);
1808   gHistNP->SetTitle("C_{-+}(#Delta#eta,#Delta#varphi)");
1809   gHistNP->GetXaxis()->SetRangeUser(-1.4,1.4);
1810   gHistNP->GetXaxis()->CenterTitle();
1811   gHistNP->GetXaxis()->SetTitleOffset(1.2);
1812   gHistNP->GetYaxis()->CenterTitle();
1813   gHistNP->GetYaxis()->SetTitleOffset(1.2);
1814   gHistNP->GetZaxis()->SetTitleOffset(1.5);
1815   TCanvas *cNP = new TCanvas("cNP","",50,50,600,500);
1816   cNP->SetFillColor(10); cNP->SetHighLightColor(10);
1817   cNP->SetLeftMargin(0.15);
1818
1819   //================//
1820   TH1D* gHistNPprojection = 0x0;
1821   Double_t sum = 0.0;
1822   Double_t gError = 0.0;
1823   Int_t nCounter = 0;
1824   if(kProjectInEta) {
1825     gHistNPprojection = new TH1D("gHistNPprojection","",gHistNP->GetNbinsX(),gHistNP->GetXaxis()->GetXmin(),gHistNP->GetXaxis()->GetXmax());
1826     for(Int_t iBinX = 1; iBinX <= gHistNP->GetNbinsX(); iBinX++) {
1827       sum = 0.; gError = 0.0; nCounter = 0;
1828       for(Int_t iBinY = binMin; iBinY <= binMax; iBinY++) {
1829       //for(Int_t iBinY = 1; iBinY <= gHistNP->GetNbinsY(); iBinY++) {
1830         sum += gHistNP->GetBinContent(iBinX,iBinY);
1831         if(gHistNP->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
1832         Double_t exy = gHistNP->GetCellError(iBinX,iBinY);
1833         gError += exy*exy;
1834       }
1835       if(nCounter != 0) {
1836         sum /= nCounter;
1837         gError = TMath::Sqrt(gError)/nCounter;
1838       }
1839       gHistNPprojection->SetBinContent(iBinX,sum);
1840       gHistNPprojection->SetBinError(iBinX,gError);
1841     }
1842     gHistNPprojection->GetXaxis()->SetRangeUser(-1.4,1.4);
1843     if(kUseZYAM)
1844       gHistNPprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#eta} - b_{ZYAM}");
1845     else
1846       gHistNPprojection->GetYaxis()->SetTitle("C_{-+}(#Delta#eta)");
1847     gHistNPprojection->GetXaxis()->SetTitle("#Delta#eta");
1848   }
1849   else {
1850     gHistNPprojection = new TH1D("gHistNPprojection","",gHistNP->GetNbinsY(),gHistNP->GetYaxis()->GetXmin(),gHistNP->GetYaxis()->GetXmax());
1851     for(Int_t iBinY = 1; iBinY <= gHistNP->GetNbinsY(); iBinY++) {
1852       sum = 0.; gError = 0.0; nCounter = 0;
1853       for(Int_t iBinX = binMin; iBinX <= binMax; iBinX++) {
1854         //for(Int_t iBinX = 1; iBinX <= gHistNP->GetNbinsX(); iBinX++) {
1855         sum += gHistNP->GetBinContent(iBinX,iBinY);
1856         if(gHistNP->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
1857         Double_t exy = gHistNP->GetCellError(iBinX,iBinY);
1858         gError += exy*exy;
1859       }
1860       if(nCounter != 0) {
1861         sum /= nCounter;
1862         gError = TMath::Sqrt(gError)/nCounter;
1863       }
1864       gHistNPprojection->SetBinContent(iBinY,sum);
1865       gHistNPprojection->SetBinError(iBinY,gError);
1866     }
1867     if(kUseZYAM)
1868       gHistNPprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#varphi} - b_{ZYAM} (rad^{-1})");
1869     else
1870       gHistNPprojection->GetYaxis()->SetTitle("C_{-+}(#Delta#varphi)");
1871     gHistNPprojection->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1872   }
1873   //ZYAM
1874   if(kUseZYAM) {
1875     Double_t reference = gHistNPprojection->GetBinContent(gHistNPprojection->GetMinimumBin());
1876   for(Int_t iBinX = 1; iBinX <= gHistNPprojection->GetNbinsX(); iBinX++)
1877     gHistNPprojection->SetBinContent(iBinX,gHistNPprojection->GetBinContent(iBinX) - reference);
1878   }
1879
1880   gHistNPprojection->GetYaxis()->SetTitleOffset(1.5);
1881   gHistNPprojection->SetMarkerStyle(20);
1882   gHistNPprojection->SetStats(kFALSE);
1883   gHistNPprojection->DrawCopy("E");
1884   //================//
1885
1886   gPad->SetTheta(30); // default is 30
1887   gPad->SetPhi(-60); // default is 30
1888   gPad->Update();
1889
1890   latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
1891   latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
1892   latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
1893   latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
1894
1895   pngName = filename;
1896   if(kProjectInEta) pngName.ReplaceAll(".root","_DeltaEtaProjection");
1897   else              pngName.ReplaceAll(".root","_DeltaPhiProjection");
1898   pngName += ".NegativePositive.png";
1899   cNP->SaveAs(pngName.Data());
1900
1901   //============================================================//
1902   //Get the ++ correlation function
1903   TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
1904   gHistPP->SetStats(kFALSE);
1905   gHistPP->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
1906   gHistPP->GetXaxis()->SetRangeUser(-1.4,1.4);
1907   gHistPP->GetXaxis()->CenterTitle();
1908   gHistPP->GetXaxis()->SetTitleOffset(1.2);
1909   gHistPP->GetYaxis()->CenterTitle();
1910   gHistPP->GetYaxis()->SetTitleOffset(1.2);
1911   gHistPP->GetZaxis()->SetTitleOffset(1.5);
1912   TCanvas *cPP = new TCanvas("cPP","",100,100,600,500);
1913   cPP->SetFillColor(10); cPP->SetHighLightColor(10);
1914   cPP->SetLeftMargin(0.15);
1915
1916   //=================//
1917   TH1D* gHistPPprojection = 0x0;
1918   Double_t sum = 0.0;
1919   Double_t gError = 0.0;
1920   Int_t nCounter = 0;
1921   if(kProjectInEta) {
1922     gHistPPprojection = new TH1D("gHistPPprojection","",gHistPP->GetNbinsX(),gHistPP->GetXaxis()->GetXmin(),gHistPP->GetXaxis()->GetXmax());
1923     for(Int_t iBinX = 1; iBinX <= gHistPP->GetNbinsX(); iBinX++) {
1924       sum = 0.; gError = 0.0; nCounter = 0;
1925       for(Int_t iBinY = binMin; iBinY <= binMax; iBinY++) {
1926         //for(Int_t iBinY = 1; iBinY <= gHistPP->GetNbinsY(); iBinY++) {
1927         sum += gHistPP->GetBinContent(iBinX,iBinY);
1928         if(gHistPP->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
1929         Double_t exy = gHistPP->GetCellError(iBinX,iBinY);
1930         gError += exy*exy;
1931       }
1932       if(nCounter != 0) {
1933         sum /= nCounter;
1934         gError = TMath::Sqrt(gError)/nCounter;
1935       }
1936       gHistPPprojection->SetBinContent(iBinX,sum);
1937       gHistPPprojection->SetBinError(iBinX,gError);
1938     }
1939     gHistPPprojection->GetXaxis()->SetRangeUser(-1.4,1.4);
1940     if(kUseZYAM)
1941       gHistPPprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#eta} - b_{ZYAM}");
1942     else
1943       gHistPPprojection->GetYaxis()->SetTitle("C_{++}(#Delta#eta)");
1944     gHistPPprojection->GetXaxis()->SetTitle("#Delta#eta");
1945   }
1946   else {
1947     gHistPPprojection = new TH1D("gHistPPprojection","",gHistPP->GetNbinsY(),gHistPP->GetYaxis()->GetXmin(),gHistPP->GetYaxis()->GetXmax());
1948     for(Int_t iBinY = 1; iBinY <= gHistPP->GetNbinsY(); iBinY++) {
1949       sum = 0.; gError = 0.0; nCounter = 0;
1950       for(Int_t iBinX = binMin; iBinX <= binMax; iBinX++) {
1951         //for(Int_t iBinX = 1; iBinX <= gHistPP->GetNbinsX(); iBinX++) {
1952         sum += gHistPP->GetBinContent(iBinX,iBinY);
1953         if(gHistPP->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
1954         Double_t exy = gHistPP->GetCellError(iBinX,iBinY);
1955         gError += exy*exy;
1956       }
1957       if(nCounter != 0) {
1958         sum /= nCounter;
1959         gError = TMath::Sqrt(gError)/nCounter;
1960       }
1961       gHistPPprojection->SetBinContent(iBinY,sum);
1962       gHistPPprojection->SetBinError(iBinY,gError);
1963     }
1964     if(kUseZYAM)
1965       gHistPPprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#varphi} - b_{ZYAM} (rad^{-1})");
1966     else
1967       gHistPPprojection->GetYaxis()->SetTitle("C_{++}(#Delta#varphi)");
1968     gHistPPprojection->GetXaxis()->SetTitle("#Delta#varphi (rad)");
1969   }
1970   //ZYAM
1971   if(kUseZYAM) {
1972     Double_t reference = gHistPPprojection->GetBinContent(gHistPPprojection->GetMinimumBin());
1973     for(Int_t iBinX = 1; iBinX <= gHistPPprojection->GetNbinsX(); iBinX++)
1974       gHistPPprojection->SetBinContent(iBinX,gHistPPprojection->GetBinContent(iBinX) - reference);
1975   }
1976
1977   gHistPPprojection->GetYaxis()->SetTitleOffset(1.5);
1978   gHistPPprojection->SetMarkerStyle(20);
1979   gHistPPprojection->SetStats(kFALSE);
1980   gHistPPprojection->DrawCopy("E");
1981   //================//
1982
1983   gPad->SetTheta(30); // default is 30
1984   gPad->SetPhi(-60); // default is 30
1985   gPad->Update();
1986
1987   latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
1988   latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
1989   latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
1990   latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
1991
1992   pngName = filename;
1993   if(kProjectInEta) pngName.ReplaceAll(".root","_DeltaEtaProjection");
1994   else              pngName.ReplaceAll(".root","_DeltaPhiProjection");
1995   pngName += ".PositivePositive.png";
1996   cPP->SaveAs(pngName.Data());
1997
1998   //============================================================//
1999   //Get the -- correlation function
2000   TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
2001   gHistNN->SetStats(kFALSE);
2002   gHistNN->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
2003   gHistNN->GetXaxis()->SetRangeUser(-1.4,1.4);
2004   gHistNN->GetXaxis()->CenterTitle();
2005   gHistNN->GetXaxis()->SetTitleOffset(1.2);
2006   gHistNN->GetYaxis()->CenterTitle();
2007   gHistNN->GetYaxis()->SetTitleOffset(1.2);
2008   gHistNN->GetZaxis()->SetTitleOffset(1.5);
2009   TCanvas *cNN = new TCanvas("cNN","",150,150,600,500);
2010   cNN->SetFillColor(10); cNN->SetHighLightColor(10);
2011   cNN->SetLeftMargin(0.15);
2012
2013   //=================//
2014   TH1D* gHistNNprojection = 0x0;
2015   Double_t sum = 0.0;
2016   Double_t gError = 0.0;
2017   Int_t nCounter = 0;
2018   if(kProjectInEta) {
2019     gHistNNprojection = new TH1D("gHistNNprojection","",gHistNN->GetNbinsX(),gHistNN->GetXaxis()->GetXmin(),gHistNN->GetXaxis()->GetXmax());
2020     for(Int_t iBinX = 1; iBinX <= gHistNN->GetNbinsX(); iBinX++) {
2021       sum = 0.; gError = 0.0; nCounter = 0;
2022       for(Int_t iBinY = binMin; iBinY <= binMax; iBinY++) {
2023         //for(Int_t iBinY = 1; iBinY <= gHistNN->GetNbinsY(); iBinY++) {
2024         sum += gHistNN->GetBinContent(iBinX,iBinY);
2025         if(gHistNN->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
2026         Double_t exy = gHistNN->GetCellError(iBinX,iBinY);
2027         gError += exy*exy;
2028       }
2029       if(nCounter != 0) {
2030         sum /= nCounter;
2031         gError = TMath::Sqrt(gError)/nCounter;
2032       }
2033       gHistNNprojection->SetBinContent(iBinX,sum);
2034       gHistNNprojection->SetBinError(iBinX,gError);
2035     }
2036     gHistNNprojection->GetXaxis()->SetRangeUser(-1.4,1.4);
2037     if(kUseZYAM)
2038       gHistNNprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#eta} - b_{ZYAM}");
2039     else
2040       gHistNNprojection->GetYaxis()->SetTitle("C_{--}(#Delta#eta)");
2041     gHistNNprojection->GetXaxis()->SetTitle("#Delta#eta");
2042   }
2043   else {
2044     gHistNNprojection = new TH1D("gHistNNprojection","",gHistNN->GetNbinsY(),gHistNN->GetYaxis()->GetXmin(),gHistNN->GetYaxis()->GetXmax());
2045     for(Int_t iBinY = 1; iBinY <= gHistNN->GetNbinsY(); iBinY++) {
2046       sum = 0.; gError = 0.0; nCounter = 0;
2047       for(Int_t iBinX = binMin; iBinX <= binMax; iBinX++) {
2048         //for(Int_t iBinX = 1; iBinX <= gHistNN->GetNbinsX(); iBinX++) {
2049         sum += gHistNN->GetBinContent(iBinX,iBinY);
2050         if(gHistNN->GetBinContent(iBinX,iBinY) != 0.) nCounter += 1;
2051         Double_t exy = gHistNN->GetCellError(iBinX,iBinY);
2052         gError += exy*exy;
2053       }
2054       if(nCounter != 0) {
2055         sum /= nCounter;
2056         gError = TMath::Sqrt(gError)/nCounter;
2057       }
2058       gHistNNprojection->SetBinContent(iBinY,sum);
2059       gHistNNprojection->SetBinError(iBinY,gError);
2060     }
2061     if(kUseZYAM)
2062       gHistNNprojection->GetYaxis()->SetTitle("#frac{1}{N_{trig}}#frac{dN_{assoc}}{#Delta#varphi} - b_{ZYAM} (rad^{-1})");
2063     else
2064       gHistNNprojection->GetYaxis()->SetTitle("C_{--}(#Delta#varphi)");
2065     gHistNNprojection->GetXaxis()->SetTitle("#Delta#varphi (rad)");
2066   }
2067   //ZYAM
2068   if(kUseZYAM) {
2069     Double_t reference = gHistNNprojection->GetBinContent(gHistNNprojection->GetMinimumBin());
2070     for(Int_t iBinX = 1; iBinX <= gHistNNprojection->GetNbinsX(); iBinX++)
2071       gHistNNprojection->SetBinContent(iBinX,gHistNNprojection->GetBinContent(iBinX) - reference);
2072   }
2073
2074   gHistNNprojection->GetYaxis()->SetTitleOffset(1.5);
2075   gHistNNprojection->SetMarkerStyle(20);
2076   gHistNNprojection->SetStats(kFALSE);
2077   gHistNNprojection->DrawCopy("E");
2078   //=================//
2079
2080   gPad->SetTheta(30); // default is 30
2081   gPad->SetPhi(-60); // default is 30
2082   gPad->Update();
2083
2084   latexInfo1->DrawLatex(0.6,0.95,centralityLatex.Data());
2085   latexInfo1->DrawLatex(0.6,0.89,psiLatex.Data());
2086   latexInfo1->DrawLatex(0.6,0.83,pttLatex.Data());
2087   latexInfo1->DrawLatex(0.6,0.77,ptaLatex.Data());
2088
2089   pngName = filename;
2090   if(kProjectInEta) pngName.ReplaceAll(".root","_DeltaEtaProjection");
2091   else              pngName.ReplaceAll(".root","_DeltaPhiProjection");
2092   pngName += ".NegativeNegative.png";
2093   cNN->SaveAs(pngName.Data());
2094
2095   TString outFileName = filename;
2096   if(kProjectInEta) outFileName.ReplaceAll(".root","_DeltaEtaProjection.root");
2097   else              outFileName.ReplaceAll(".root","_DeltaPhiProjection.root");
2098   TFile *fProjection = TFile::Open(outFileName.Data(),"recreate");
2099   gHistNPprojection->Write();
2100   gHistPNprojection->Write();
2101   gHistNNprojection->Write();
2102   gHistPPprojection->Write();
2103   fProjection->Close();
2104 }
2105
2106 //____________________________________________________________//
2107 void fitCorrelationFunctions(Int_t gCentrality = 1,
2108                              Double_t psiMin = -0.5, Double_t psiMax = 3.5,
2109                              Double_t vertexZMin = -10.,Double_t vertexZMax = 10.,
2110                              Double_t ptTriggerMin = -1.,
2111                              Double_t ptTriggerMax = -1.,
2112                              Double_t ptAssociatedMin = -1.,
2113                              Double_t ptAssociatedMax = -1.,
2114                              TH2D *gHist) {
2115
2116   cout<<"FITTING FUNCTION"<<endl;
2117
2118   //near side peak: [1]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4]))
2119   //away side ridge: [5]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((y-TMath::Pi())/[6]),2)),[7]))
2120   //longitudinal ridge: [8]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[9]),2)),[10]))
2121   //wing structures: [11]*TMath::Power(x,2)
2122   //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))
2123   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.);
2124   gFitFunction->SetName("gFitFunction");
2125
2126
2127   //Normalization
2128   gFitFunction->SetParName(0,"N1");
2129   //near side peak
2130   gFitFunction->SetParName(1,"N_{near side}");gFitFunction->FixParameter(1,0.0);
2131   gFitFunction->SetParName(2,"Sigma_{near side}(delta eta)"); gFitFunction->FixParameter(2,0.0);
2132   gFitFunction->SetParName(3,"Sigma_{near side}(delta phi)"); gFitFunction->FixParameter(3,0.0);
2133   gFitFunction->SetParName(4,"Exponent_{near side}"); gFitFunction->FixParameter(4,1.0);
2134   //away side ridge
2135   gFitFunction->SetParName(5,"N_{away side}"); gFitFunction->FixParameter(5,0.0);
2136   gFitFunction->SetParName(6,"Sigma_{away side}(delta phi)"); gFitFunction->FixParameter(6,0.0);
2137   gFitFunction->SetParName(7,"Exponent_{away side}"); gFitFunction->FixParameter(7,1.0);
2138   //longitudinal ridge
2139   gFitFunction->SetParName(8,"N_{long. ridge}"); gFitFunction->SetParameter(8,0.05);//
2140   gFitFunction->FixParameter(8,0.0);
2141   gFitFunction->SetParName(9,"Sigma_{long. ridge}(delta eta)"); gFitFunction->FixParameter(9,0.0);
2142   gFitFunction->SetParName(10,"Exponent_{long. ridge}"); gFitFunction->FixParameter(10,1.0);
2143   //wing structures
2144   gFitFunction->SetParName(11,"N_{wing}"); gFitFunction->FixParameter(11,0.0);
2145   //flow contribution
2146   gFitFunction->SetParName(12,"N_{flow}"); gFitFunction->SetParameter(12,0.2);gFitFunction->SetParLimits(12,0.0,10);
2147   gFitFunction->SetParName(13,"V1"); gFitFunction->SetParameter(13,0.005);gFitFunction->SetParLimits(13,0.0,10);
2148   gFitFunction->SetParName(14,"V2"); gFitFunction->SetParameter(14,0.1);gFitFunction->SetParLimits(14,0.0,10);
2149   gFitFunction->SetParName(15,"V3"); gFitFunction->SetParameter(15,0.05);gFitFunction->SetParLimits(15,0.0,10);
2150   gFitFunction->SetParName(16,"V4"); gFitFunction->SetParameter(16,0.005);gFitFunction->SetParLimits(16,0.0,10);
2151   gFitFunction->SetParName(17,"Offset"); gFitFunction->FixParameter(17,0.0);
2152
2153   // flow parameters
2154   Double_t fNV = 0.;
2155   Double_t fV1 = 0.;
2156   Double_t fV2 = 0.;
2157   Double_t fV3 = 0.;
2158   Double_t fV4 = 0.;
2159
2160   //Fitting the correlation function (first the edges to extract flow)
2161   gHist->Fit("gFitFunction","nm","",1.0,1.6);
2162
2163   fNV += gFitFunction->GetParameter(12);
2164   fV1 += gFitFunction->GetParameter(13);
2165   fV2 += gFitFunction->GetParameter(14);
2166   fV3 += gFitFunction->GetParameter(15);
2167   fV4 += gFitFunction->GetParameter(16);
2168
2169   gHist->Fit("gFitFunction","nm","",-1.6,-1.0);
2170
2171   fNV += gFitFunction->GetParameter(12);
2172   fV1 += gFitFunction->GetParameter(13);
2173   fV2 += gFitFunction->GetParameter(14);
2174   fV3 += gFitFunction->GetParameter(15);
2175   fV4 += gFitFunction->GetParameter(16);
2176
2177   // Now fit the whole with fixed flow
2178   gFitFunction->FixParameter(12,fNV/2.);
2179   gFitFunction->FixParameter(13,fV1/2.);
2180   gFitFunction->FixParameter(14,fV2/2.);
2181   gFitFunction->FixParameter(15,fV3/2.);
2182   gFitFunction->FixParameter(16,fV4/2.);
2183
2184   gFitFunction->ReleaseParameter(0);gFitFunction->SetParameter(0,1.0);
2185   gFitFunction->ReleaseParameter(1);gFitFunction->SetParameter(1,0.3);
2186   gFitFunction->ReleaseParameter(2);gFitFunction->SetParameter(2,0.3);gFitFunction->SetParLimits(2,0.05,0.7);
2187   gFitFunction->ReleaseParameter(3);gFitFunction->SetParameter(3,0.3);gFitFunction->SetParLimits(3,0.05,1.7);
2188   //gFitFunction->ReleaseParameter(4);gFitFunction->SetParameter(4,1.0);gFitFunction->SetParLimits(4,0.0,2.0);
2189   gFitFunction->ReleaseParameter(5);gFitFunction->SetParameter(5,1.0);//gFitFunction->SetParLimits(5,0.0,10);
2190   gFitFunction->ReleaseParameter(6);gFitFunction->SetParameter(6,0.5);//gFitFunction->SetParLimits(6,0.0,10);
2191   //gFitFunction->ReleaseParameter(7);gFitFunction->SetParameter(7,1.0);gFitFunction->SetParLimits(7,0.0,2.0);
2192   //gFitFunction->ReleaseParameter(8);gFitFunction->SetParameter(8,0.05);
2193   //gFitFunction->ReleaseParameter(9);gFitFunction->SetParameter(9,0.6);gFitFunction->SetParLimits(9,0.1,10.0);
2194   //gFitFunction->ReleaseParameter(10);gFitFunction->SetParameter(10,1.0);gFitFunction->SetParLimits(10,0.0,2.0);
2195   gFitFunction->ReleaseParameter(17);gFitFunction->SetParameter(17,0.7);gFitFunction->SetParLimits(17,0.0,0.9);
2196
2197   gHist->Fit("gFitFunction","nm");
2198
2199
2200   //Cloning the histogram
2201   TString histoName = gHist->GetName(); histoName += "Fit";
2202   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());
2203   TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
2204   gHistResidual->SetName("gHistResidual");
2205   gHistResidual->Sumw2();
2206
2207   for(Int_t iBinDeltaEta = 1; iBinDeltaEta <= gHist->GetNbinsX(); iBinDeltaEta++) {
2208     for(Int_t iBinDeltaPhi = 1; iBinDeltaPhi <= gHist->GetNbinsY(); iBinDeltaPhi++) {
2209       gHistFit->SetBinContent(iBinDeltaEta,iBinDeltaPhi,gFitFunction->Eval(gHist->GetXaxis()->GetBinCenter(iBinDeltaEta),gHist->GetYaxis()->GetBinCenter(iBinDeltaPhi)));
2210     }
2211   }
2212   gHistResidual->Add(gHistFit,-1);
2213
2214   //Write to output file
2215   TString newFileName = "correlationFunctionFit";
2216   if(histoName.Contains("PN")) newFileName += "PN";
2217   else if(histoName.Contains("NP")) newFileName += "NP";
2218   else if(histoName.Contains("PP")) newFileName += "PP";
2219   else if(histoName.Contains("NN")) newFileName += "NN";
2220   newFileName += ".Centrality";
2221   newFileName += gCentrality; newFileName += ".Psi";
2222   if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
2223   else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
2224   else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
2225   else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
2226   else newFileName += "All.PttFrom";
2227   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
2228   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
2229   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
2230   newFileName += Form("%.1f",ptAssociatedMax);
2231
2232   newFileName += "_";
2233   newFileName += Form("%.1f",psiMin);
2234   newFileName += "-";
2235   newFileName += Form("%.1f",psiMax);
2236
2237   newFileName += ".root";
2238   TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
2239   gHist->Write();
2240   gHistFit->Write();
2241   gHistResidual->Write();
2242   gFitFunction->Write();
2243   newFile->Close();
2244
2245
2246 }
2247
2248 //____________________________________________________________//
2249 TH2D *convolute2D(TH2D* h1, TH2D* h2, TString hname) {
2250   //
2251   // this function does the convolution of 2 eta or phi "efficiencies" in a deta or dphi "efficiency"
2252   // and returns a new histogram which is normalized to the number of bin combinations
2253
2254   // histogram definition
2255   TH2D *hConv = NULL;
2256   hConv = new TH2D(hname.Data(),hname.Data(), h1->GetNbinsY(),-2,2,h1->GetNbinsX(),-TMath::Pi()/2.,3*TMath::Pi()/2.);
2257
2258   Double_t x1 = 0.;
2259   Double_t x2 = 0.;
2260   Double_t y1 = 0.;
2261   Double_t y2 = 0.;
2262   Double_t z1 = 0.;
2263   Double_t z2 = 0.;
2264   Double_t dx = 0.;
2265   Double_t dy = 0.;
2266   Double_t dz = 0.;
2267
2268
2269   // convolution
2270   for(Int_t i = 0; i < h1->GetNbinsX(); i ++){
2271     cout<<"CONVOLUTING BIN = "<<i<<" of "<<h1->GetNbinsX()<<endl;
2272     for(Int_t j = 0; j < h1->GetNbinsY(); j ++){
2273       for(Int_t k = 0; k < h2->GetNbinsX(); k ++){
2274         for(Int_t l = 0; l < h2->GetNbinsY(); l ++){
2275
2276           x1 = (Double_t)h1->GetXaxis()->GetBinCenter(i+1);
2277           y1 = (Double_t)h1->GetYaxis()->GetBinCenter(j+1);
2278           x2 = (Double_t)h2->GetXaxis()->GetBinCenter(k+1);
2279           y2 = (Double_t)h2->GetYaxis()->GetBinCenter(l+1);
2280
2281           z1 = (Double_t)h1->GetBinContent(i+1,j+1);
2282           z2 = (Double_t)h2->GetBinContent(k+1,l+1);
2283
2284           // need the gymnastics to keep the same binning
2285           dx = x1 - x2 - (h1->GetXaxis()->GetBinWidth(1)/2.);
2286           if(gRandom->Gaus() > 0) dy = y1 - y2 + (h1->GetYaxis()->GetBinWidth(1)/2.);
2287           else                    dy = y1 - y2 - (h1->GetYaxis()->GetBinWidth(1)/2.);
2288
2289           if(dx>3./2.*TMath::Pi())  dx = dx - 2.*TMath::Pi();
2290           if(dx<-1./2.*TMath::Pi()) dx = 2*TMath::Pi() + dx;
2291
2292           dz = z1*z2;
2293
2294           hConv->Fill(dy,dx,dz);
2295
2296         }
2297       }
2298     }
2299   }
2300
2301   return hConv;
2302
2303 }