PWGCF/EBYE/macros/drawCorrelationFunctionPsi.C

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