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