PWGCF/EBYE/macros/drawCorrelationFunctionPsi.C

   1 const Int_t numberOfCentralityBins = 8;
   2 TString centralityArray[numberOfCentralityBins] = {"0-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80"};
   3
   4 const Int_t gRebin = 1;
   5
   6 //____________________________________________________________//
   7 void drawCorrelationFunctionPsiAllPtCombinations(const char* filename = "AnalysisResults.root",
   8                                                  Int_t gCentrality = 1,
   9                                                  Int_t gBit = -1,
  10                                                  const char* gCentralityEstimator = 0x0,
  11                                                  Bool_t kShowShuffled = kFALSE,
  12                                                  Bool_t kShowMixed = kTRUE,
  13                                                  Double_t psiMin = -0.5,
  14                                                  Double_t psiMax = 3.5){
  15
  16   // this could also be retrieved directly from AliBalancePsi
  17   const Int_t kNPtBins = 16;
  18   Double_t ptBins[kNPtBins+1] = {0.2,0.6,1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,10.,12.,15.,20.};
  19
  20   cout<<"You have chosen to do all pT combinations --> this could take some time."<<endl;
  21
  22   for(Int_t iTrig = 0; iTrig < 2/*kNPtBins*/; iTrig++){
  23     for(Int_t iAssoc = 0; iAssoc < 2/*kNPtBins*/; iAssoc++){
  24       cout<<"================================================================="<<endl;
  25       cout<<"PROCESS NOW: "<<endl;
  26       cout<<" -> "<< ptBins[iTrig]<<" < pTtrig < "<<ptBins[iTrig+1]<<"   "<<ptBins[iAssoc]<<" < pTassoc < "<<ptBins[iAssoc+1]<<endl;
  27       drawCorrelationFunctionPsi(filename,gCentrality,gBit,gCentralityEstimator,kShowShuffled,kShowMixed,psiMin,psiMax,ptBins[iTrig],ptBins[iTrig+1],ptBins[iAssoc],ptBins[iAssoc+1]);
  28       cout<<"================================================================="<<endl;
  29       cout<<endl;
  30     }
  31   }
  32
  33 }
  34
  35 //____________________________________________________________//
  36 void drawCorrelationFunctionsAllPtCombinations(const char* lhcPeriod = "LHC11h",
  37                                                Int_t gTrainID = 208,
  38                                                Int_t gCentrality = 1,
  39                                                Double_t psiMin = -0.5, Double_t psiMax = 3.5)
  40 {
  41
  42  // this could also be retrieved directly from AliBalancePsi
  43   const Int_t kNPtBins = 16;
  44   Double_t ptBins[kNPtBins+1] = {0.2,0.6,1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,10.,12.,15.,20.};
  45
  46   cout<<"You have chosen to do all pT combinations --> this could take some time."<<endl;
  47
  48   for(Int_t iTrig = 0; iTrig < kNPtBins; iTrig++){
  49     for(Int_t iAssoc = 0; iAssoc < kNPtBins; iAssoc++){
  50       cout<<"================================================================="<<endl;
  51       cout<<"FIT NOW: "<<endl;
  52       cout<<" -> "<< ptBins[iTrig]<<" < pTtrig < "<<ptBins[iTrig+1]<<"   "<<ptBins[iAssoc]<<" < pTassoc < "<<ptBins[iAssoc+1]<<endl;
  53       drawCorrelationFunctions(lhcPeriod,gTrainID,gCentrality,psiMin,psiMax,ptBins[iTrig],ptBins[iTrig+1],ptBins[iAssoc],ptBins[iAssoc+1]);
  54       cout<<"================================================================="<<endl;
  55       cout<<endl;
  56     }
  57   }
  58
  59 }
  60
  61
  62 void drawCorrelationFunctionPsi(const char* filename = "AnalysisResults.root",
  63                                 Int_t gCentrality = 1,
  64                                 Int_t gBit = -1,
  65                                 const char* gCentralityEstimator = 0x0,
  66                                 Bool_t kShowShuffled = kFALSE,
  67                                 Bool_t kShowMixed = kTRUE,
  68                                 Double_t psiMin = -0.5,
  69                                 Double_t psiMax = 3.5,
  70                                 Double_t ptTriggerMin = -1.,
  71                                 Double_t ptTriggerMax = -1.,
  72                                 Double_t ptAssociatedMin = -1.,
  73                                 Double_t ptAssociatedMax = -1.) {
  74   //Macro that draws the correlation functions from the balance function
  75   //analysis vs the reaction plane
  76   //Author: Panos.Christakoglou@nikhef.nl
  77   gROOT->LoadMacro("~/SetPlotStyle.C");
  78   SetPlotStyle();
  79   gStyle->SetPalette(1,0);
  80
  81   //Load the PWG2ebye library
  82   gSystem->Load("libANALYSIS.so");
  83   gSystem->Load("libANALYSISalice.so");
  84   gSystem->Load("libEventMixing.so");
  85   gSystem->Load("libCORRFW.so");
  86   gSystem->Load("libPWGTools.so");
  87   gSystem->Load("libPWGCFebye.so");
  88
  89   //Prepare the objects and return them
  90   TList *list = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,0);
  91   TList *listShuffled = NULL;
  92   if(kShowShuffled) listShuffled = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,1);
  93   TList *listMixed = NULL;
  94   if(kShowMixed) listMixed = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,2);
  95
  96   if(!list) {
  97     Printf("The TList object was not created");
  98     return;
  99   }
 100   else
 101     draw(list,listShuffled,listMixed,
 102          gCentralityEstimator,gCentrality,psiMin,psiMax,
 103          ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 104 }
 105
 106 //______________________________________________________//
 107 TList *GetListOfObjects(const char* filename,
 108                         Int_t gCentrality,
 109                         Int_t gBit,
 110                         const char *gCentralityEstimator,
 111                         Int_t kData = 1) {
 112   //Get the TList objects (QA, bf, bf shuffled)
 113   TList *listBF = 0x0;
 114
 115   //Open the file
 116   TFile *f = TFile::Open(filename,"UPDATE");
 117   if((!f)||(!f->IsOpen())) {
 118     Printf("The file %s is not found. Aborting...",filename);
 119     return listBF;
 120   }
 121   //f->ls();
 122
 123   TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
 124   if(!dir) {
 125     Printf("The TDirectoryFile is not found. Aborting...",filename);
 126     return listBF;
 127   }
 128   //dir->ls();
 129
 130   TString listBFName;
 131   if(kData == 0) {
 132     //cout<<"no shuffling - no mixing"<<endl;
 133     listBFName = "listBFPsi_";
 134   }
 135   else if(kData == 1) {
 136     //cout<<"shuffling - no mixing"<<endl;
 137     listBFName = "listBFPsiShuffled_";
 138   }
 139   else if(kData == 2) {
 140     //cout<<"no shuffling - mixing"<<endl;
 141     listBFName = "listBFPsiMixed_";
 142   }
 143   listBFName += centralityArray[gCentrality-1];
 144   if(gBit > -1) {
 145     listBFName += "_Bit"; listBFName += gBit; }
 146   if(gCentralityEstimator) {
 147     listBFName += "_"; listBFName += gCentralityEstimator;}
 148
 149   // histograms were already retrieved (in first iteration)
 150   if(dir->Get(Form("%s_histograms",listBFName.Data()))){
 151     listBF = dynamic_cast<TList *>(dir->Get(Form("%s_histograms",listBFName.Data())));
 152   }
 153
 154   // histograms were not yet retrieved (this is the first iteration)
 155   else{
 156
 157     listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
 158     cout<<"======================================================="<<endl;
 159     cout<<"List name: "<<listBF->GetName()<<endl;
 160     //listBF->ls();
 161
 162     //Get the histograms
 163     TString histoName;
 164     if(kData == 0)
 165       histoName = "fHistP";
 166     else if(kData == 1)
 167       histoName = "fHistP_shuffle";
 168     else if(kData == 2)
 169       histoName = "fHistP";
 170     if(gCentralityEstimator)
 171       histoName += gCentralityEstimator;
 172     AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 173     if(!fHistP) {
 174       Printf("fHistP %s not found!!!",histoName.Data());
 175       break;
 176     }
 177     fHistP->FillParent(); fHistP->DeleteContainers();
 178
 179     if(kData == 0)
 180       histoName = "fHistN";
 181     if(kData == 1)
 182       histoName = "fHistN_shuffle";
 183     if(kData == 2)
 184       histoName = "fHistN";
 185     if(gCentralityEstimator)
 186       histoName += gCentralityEstimator;
 187     AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 188     if(!fHistN) {
 189       Printf("fHistN %s not found!!!",histoName.Data());
 190       break;
 191     }
 192     fHistN->FillParent(); fHistN->DeleteContainers();
 193
 194     if(kData == 0)
 195       histoName = "fHistPN";
 196     if(kData == 1)
 197       histoName = "fHistPN_shuffle";
 198     if(kData == 2)
 199       histoName = "fHistPN";
 200     if(gCentralityEstimator)
 201       histoName += gCentralityEstimator;
 202     AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 203     if(!fHistPN) {
 204       Printf("fHistPN %s not found!!!",histoName.Data());
 205       break;
 206     }
 207     fHistPN->FillParent(); fHistPN->DeleteContainers();
 208
 209     if(kData == 0)
 210       histoName = "fHistNP";
 211     if(kData == 1)
 212       histoName = "fHistNP_shuffle";
 213     if(kData == 2)
 214       histoName = "fHistNP";
 215     if(gCentralityEstimator)
 216       histoName += gCentralityEstimator;
 217     AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 218     if(!fHistNP) {
 219       Printf("fHistNP %s not found!!!",histoName.Data());
 220       break;
 221     }
 222     fHistNP->FillParent(); fHistNP->DeleteContainers();
 223
 224     if(kData == 0)
 225       histoName = "fHistPP";
 226     if(kData == 1)
 227       histoName = "fHistPP_shuffle";
 228     if(kData == 2)
 229       histoName = "fHistPP";
 230     if(gCentralityEstimator)
 231       histoName += gCentralityEstimator;
 232     AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 233     if(!fHistPP) {
 234       Printf("fHistPP %s not found!!!",histoName.Data());
 235       break;
 236     }
 237     fHistPP->FillParent(); fHistPP->DeleteContainers();
 238
 239     if(kData == 0)
 240       histoName = "fHistNN";
 241     if(kData == 1)
 242       histoName = "fHistNN_shuffle";
 243     if(kData == 2)
 244       histoName = "fHistNN";
 245     if(gCentralityEstimator)
 246       histoName += gCentralityEstimator;
 247     AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 248     if(!fHistNN) {
 249       Printf("fHistNN %s not found!!!",histoName.Data());
 250       break;
 251     }
 252     fHistNN->FillParent(); fHistNN->DeleteContainers();
 253
 254     dir->cd();
 255     listBF->Write(Form("%s_histograms",listBFName.Data()), TObject::kSingleKey);
 256
 257   }// first iteration
 258
 259   f->Close();
 260
 261   return listBF;
 262 }
 263
 264 //______________________________________________________//
 265 void draw(TList *list, TList *listBFShuffled, TList *listBFMixed,
 266           const char *gCentralityEstimator,
 267           Int_t gCentrality, Double_t psiMin, Double_t psiMax,
 268           Double_t ptTriggerMin, Double_t ptTriggerMax,
 269           Double_t ptAssociatedMin, Double_t ptAssociatedMax) {
 270   //Draws the correlation functions for every centrality bin
 271   //(+-), (-+), (++), (--)
 272   AliTHn *hP = NULL;
 273   AliTHn *hN = NULL;
 274   AliTHn *hPN = NULL;
 275   AliTHn *hNP = NULL;
 276   AliTHn *hPP = NULL;
 277   AliTHn *hNN = NULL;
 278
 279   TString gHistPname = "fHistP";
 280   if(gCentralityEstimator) gHistPname += gCentralityEstimator;
 281   TString gHistNname = "fHistN";
 282   if(gCentralityEstimator) gHistNname += gCentralityEstimator;
 283   TString gHistPNname = "fHistPN";
 284   if(gCentralityEstimator) gHistPNname += gCentralityEstimator;
 285   TString gHistNPname = "fHistNP";
 286   if(gCentralityEstimator) gHistNPname += gCentralityEstimator;
 287   TString gHistPPname = "fHistPP";
 288   if(gCentralityEstimator) gHistPPname += gCentralityEstimator;
 289   TString gHistNNname = "fHistNN";
 290   if(gCentralityEstimator) gHistNNname += gCentralityEstimator;
 291
 292   hP = (AliTHn*) list->FindObject(gHistPname.Data());
 293   hN = (AliTHn*) list->FindObject(gHistNname.Data());
 294   hPN = (AliTHn*) list->FindObject(gHistPNname.Data());
 295   hNP = (AliTHn*) list->FindObject(gHistNPname.Data());
 296   hPP = (AliTHn*) list->FindObject(gHistPPname.Data());
 297   hNN = (AliTHn*) list->FindObject(gHistNNname.Data());
 298
 299   //Create the AliBalancePsi object and fill it with the AliTHn objects
 300   AliBalancePsi *b = new AliBalancePsi();
 301   b->SetHistNp(hP);
 302   b->SetHistNn(hN);
 303   b->SetHistNpn(hPN);
 304   b->SetHistNnp(hNP);
 305   b->SetHistNpp(hPP);
 306   b->SetHistNnn(hNN);
 307
 308   //balance function shuffling
 309   AliTHn *hPShuffled = NULL;
 310   AliTHn *hNShuffled = NULL;
 311   AliTHn *hPNShuffled = NULL;
 312   AliTHn *hNPShuffled = NULL;
 313   AliTHn *hPPShuffled = NULL;
 314   AliTHn *hNNShuffled = NULL;
 315   if(listBFShuffled) {
 316     //listBFShuffled->ls();
 317
 318     gHistPname = "fHistP_shuffle";
 319     if(gCentralityEstimator) gHistPname += gCentralityEstimator;
 320     gHistNname = "fHistN_shuffle";
 321     if(gCentralityEstimator) gHistNname += gCentralityEstimator;
 322     gHistPNname = "fHistPN_shuffle";
 323     if(gCentralityEstimator) gHistPNname += gCentralityEstimator;
 324     gHistNPname = "fHistNP_shuffle";
 325     if(gCentralityEstimator) gHistNPname += gCentralityEstimator;
 326     gHistPPname = "fHistPP_shuffle";
 327     if(gCentralityEstimator) gHistPPname += gCentralityEstimator;
 328     gHistNNname = "fHistNN_shuffle";
 329     if(gCentralityEstimator) gHistNNname += gCentralityEstimator;
 330
 331     hPShuffled = (AliTHn*) listBFShuffled->FindObject(gHistPname.Data());
 332     hPShuffled->SetName("gHistPShuffled");
 333     hNShuffled = (AliTHn*) listBFShuffled->FindObject(gHistNname.Data());
 334     hNShuffled->SetName("gHistNShuffled");
 335     hPNShuffled = (AliTHn*) listBFShuffled->FindObject(gHistPNname.Data());
 336     hPNShuffled->SetName("gHistPNShuffled");
 337     hNPShuffled = (AliTHn*) listBFShuffled->FindObject(gHistNPname.Data());
 338     hNPShuffled->SetName("gHistNPShuffled");
 339     hPPShuffled = (AliTHn*) listBFShuffled->FindObject(gHistPPname.Data());
 340     hPPShuffled->SetName("gHistPPShuffled");
 341     hNNShuffled = (AliTHn*) listBFShuffled->FindObject(gHistNNname.Data());
 342     hNNShuffled->SetName("gHistNNShuffled");
 343
 344     AliBalancePsi *bShuffled = new AliBalancePsi();
 345     bShuffled->SetHistNp(hPShuffled);
 346     bShuffled->SetHistNn(hNShuffled);
 347     bShuffled->SetHistNpn(hPNShuffled);
 348     bShuffled->SetHistNnp(hNPShuffled);
 349     bShuffled->SetHistNpp(hPPShuffled);
 350     bShuffled->SetHistNnn(hNNShuffled);
 351   }
 352
 353   //balance function mixing
 354   AliTHn *hPMixed = NULL;
 355   AliTHn *hNMixed = NULL;
 356   AliTHn *hPNMixed = NULL;
 357   AliTHn *hNPMixed = NULL;
 358   AliTHn *hPPMixed = NULL;
 359   AliTHn *hNNMixed = NULL;
 360
 361   if(listBFMixed) {
 362     //listBFMixed->ls();
 363
 364     gHistPname = "fHistP";
 365     if(gCentralityEstimator) gHistPname += gCentralityEstimator;
 366     gHistNname = "fHistN";
 367     if(gCentralityEstimator) gHistNname += gCentralityEstimator;
 368     gHistPNname = "fHistPN";
 369     if(gCentralityEstimator) gHistPNname += gCentralityEstimator;
 370     gHistNPname = "fHistNP";
 371     if(gCentralityEstimator) gHistNPname += gCentralityEstimator;
 372     gHistPPname = "fHistPP";
 373     if(gCentralityEstimator) gHistPPname += gCentralityEstimator;
 374     gHistNNname = "fHistNN";
 375     if(gCentralityEstimator) gHistNNname += gCentralityEstimator;
 376     hPMixed = (AliTHn*) listBFMixed->FindObject(gHistPname.Data());
 377     hPMixed->SetName("gHistPMixed");
 378     hNMixed = (AliTHn*) listBFMixed->FindObject(gHistNname.Data());
 379     hNMixed->SetName("gHistNMixed");
 380     hPNMixed = (AliTHn*) listBFMixed->FindObject(gHistPNname.Data());
 381     hPNMixed->SetName("gHistPNMixed");
 382     hNPMixed = (AliTHn*) listBFMixed->FindObject(gHistNPname.Data());
 383     hNPMixed->SetName("gHistNPMixed");
 384     hPPMixed = (AliTHn*) listBFMixed->FindObject(gHistPPname.Data());
 385     hPPMixed->SetName("gHistPPMixed");
 386     hNNMixed = (AliTHn*) listBFMixed->FindObject(gHistNNname.Data());
 387     hNNMixed->SetName("gHistNNMixed");
 388
 389     AliBalancePsi *bMixed = new AliBalancePsi();
 390     bMixed->SetHistNp(hPMixed);
 391     bMixed->SetHistNn(hNMixed);
 392     bMixed->SetHistNpn(hPNMixed);
 393     bMixed->SetHistNnp(hNPMixed);
 394     bMixed->SetHistNpp(hPPMixed);
 395     bMixed->SetHistNnn(hNNMixed);
 396   }
 397
 398   TH2D *gHistPN[4];
 399   TH2D *gHistNP[4];
 400   TH2D *gHistPP[4];
 401   TH2D *gHistNN[4];
 402
 403   TCanvas *cPN[4];
 404   TCanvas *cNP[4];
 405   TCanvas *cPP[4];
 406   TCanvas *cNN[4];
 407   TString histoTitle, pngName;
 408
 409   //(+-)
 410   histoTitle = "(+-) | Centrality: ";
 411   histoTitle += centralityArray[gCentrality-1];
 412   histoTitle += "%";
 413   if((psiMin == -0.5)&&(psiMax == 0.5))
 414     histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 415   else if((psiMin == 0.5)&&(psiMax == 1.5))
 416     histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 417   else if((psiMin == 1.5)&&(psiMax == 2.5))
 418     histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 419   else
 420     histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 421
 422   gHistPN[0] = b->GetCorrelationFunctionPN(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 423   gHistPN[0]->GetYaxis()->SetTitleOffset(1.5);
 424   gHistPN[0]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 425   gHistPN[0]->SetTitle(histoTitle.Data());
 426   cPN[0] = new TCanvas("cPN0","",0,0,600,500);
 427   cPN[0]->SetFillColor(10); cPN[0]->SetHighLightColor(10);
 428   gHistPN[0]->DrawCopy("surf1fb");
 429   gPad->SetTheta(30); // default is 30
 430   //gPad->SetPhi(130); // default is 30
 431   gPad->SetPhi(-60); // default is 30
 432   gPad->Update();
 433   pngName = "DeltaPhiDeltaEta.Centrality";
 434   pngName += centralityArray[gCentrality-1];
 435   pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 436   pngName += ".PositiveNegative.png";
 437   //cPN[0]->SaveAs(pngName.Data());
 438
 439   if(listBFShuffled) {
 440     histoTitle = "(+-) shuffled | Centrality: ";
 441     histoTitle += centralityArray[gCentrality-1];
 442     histoTitle += "%";
 443     if((psiMin == -0.5)&&(psiMax == 0.5))
 444       histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 445     else if((psiMin == 0.5)&&(psiMax == 1.5))
 446       histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 447     else if((psiMin == 1.5)&&(psiMax == 2.5))
 448       histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 449     else
 450       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 451
 452     gHistPN[1] = bShuffled->GetCorrelationFunctionPN(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 453     gHistPN[1]->GetYaxis()->SetTitleOffset(1.5);
 454     gHistPN[1]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 455     gHistPN[1]->SetTitle(histoTitle.Data());
 456     cPN[1] = new TCanvas("cPN1","",0,100,600,500);
 457     cPN[1]->SetFillColor(10);
 458     cPN[1]->SetHighLightColor(10);
 459     gHistPN[1]->DrawCopy("surf1fb");
 460     gPad->SetTheta(30); // default is 30
 461     //gPad->SetPhi(130); // default is 30
 462     gPad->SetPhi(-60); // default is 30
 463     gPad->Update();
 464     pngName = "DeltaPhiDeltaEtaShuffled.Centrality";
 465     pngName += centralityArray[gCentrality-1];
 466     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 467     pngName += ".PositiveNegative.png";
 468     //cPN[1]->SaveAs(pngName.Data());
 469   }
 470
 471   if(listBFMixed) {
 472     histoTitle = "(+-) mixed | Centrality: ";
 473     histoTitle += centralityArray[gCentrality-1];
 474     histoTitle += "%";
 475     if((psiMin == -0.5)&&(psiMax == 0.5))
 476       histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 477     else if((psiMin == 0.5)&&(psiMax == 1.5))
 478       histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 479     else if((psiMin == 1.5)&&(psiMax == 2.5))
 480       histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 481     else
 482       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 483
 484     gHistPN[2] = bMixed->GetCorrelationFunctionPN(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 485     gHistPN[2]->GetYaxis()->SetTitleOffset(1.5);
 486     gHistPN[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 487     gHistPN[2]->SetTitle(histoTitle.Data());
 488     cPN[2] = new TCanvas("cPN2","",0,200,600,500);
 489     cPN[2]->SetFillColor(10);
 490     cPN[2]->SetHighLightColor(10);
 491     gHistPN[2]->DrawCopy("surf1fb");
 492     gPad->SetTheta(30); // default is 30
 493     //gPad->SetPhi(130); // default is 30
 494     gPad->SetPhi(-60); // default is 30
 495     gPad->Update();
 496     pngName = "DeltaPhiDeltaEtaMixed.Centrality";
 497     pngName += centralityArray[gCentrality-1];
 498     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 499     pngName += ".PositiveNegative.png";
 500     //cPN[2]->SaveAs(pngName.Data());
 501
 502     //Correlation function (+-)
 503     gHistPN[3] = dynamic_cast<TH2D *>(gHistPN[0]->Clone());
 504     gHistPN[3]->Divide(gHistPN[2]);
 505     gHistPN[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
 506     gHistPN[3]->GetZaxis()->SetTitle("C_{+-}(#Delta#eta,#Delta#varphi)");
 507     cPN[3] = new TCanvas("cPN3","",0,300,600,500);
 508     cPN[3]->SetFillColor(10);
 509     cPN[3]->SetHighLightColor(10);
 510     gHistPN[3]->DrawCopy("surf1fb");
 511     gPad->SetTheta(30); // default is 30
 512     //gPad->SetPhi(130); // default is 30
 513     gPad->SetPhi(-60); // default is 30
 514     gPad->Update();
 515     pngName = "CorrelationFunction.Centrality";
 516     pngName += centralityArray[gCentrality-1];
 517     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 518     pngName += ".PositiveNegative.png";
 519     //cPN[3]->SaveAs(pngName.Data());
 520   }
 521
 522   //(-+)
 523   histoTitle = "(-+) | Centrality: ";
 524   histoTitle += centralityArray[gCentrality-1];
 525   histoTitle += "%";
 526   if((psiMin == -0.5)&&(psiMax == 0.5))
 527     histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 528   else if((psiMin == 0.5)&&(psiMax == 1.5))
 529     histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 530   else if((psiMin == 1.5)&&(psiMax == 2.5))
 531     histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 532   else
 533     histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 534
 535   gHistNP[0] = b->GetCorrelationFunctionNP(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 536   gHistNP[0]->GetYaxis()->SetTitleOffset(1.5);
 537   gHistNP[0]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 538   gHistNP[0]->SetTitle(histoTitle.Data());
 539   cNP[0] = new TCanvas("cNP0","",100,0,600,500);
 540   cNP[0]->SetFillColor(10);
 541   cNP[0]->SetHighLightColor(10);
 542   gHistNP[0]->DrawCopy("surf1fb");
 543   gPad->SetTheta(30); // default is 30
 544   //gPad->SetPhi(130); // default is 30
 545   gPad->SetPhi(-60); // default is 30
 546   gPad->Update();
 547   pngName = "DeltaPhiDeltaEta.Centrality";
 548   pngName += centralityArray[gCentrality-1];
 549   pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 550   pngName += ".NegativePositive.png";
 551   //cNP[0]->SaveAs(pngName.Data());
 552
 553   if(listBFShuffled) {
 554     histoTitle = "(-+) shuffled | Centrality: ";
 555     histoTitle += centralityArray[gCentrality-1];
 556     histoTitle += "%";
 557     if((psiMin == -0.5)&&(psiMax == 0.5))
 558       histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 559     else if((psiMin == 0.5)&&(psiMax == 1.5))
 560       histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 561     else if((psiMin == 1.5)&&(psiMax == 2.5))
 562       histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 563     else
 564       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 565
 566     gHistNP[1] = bShuffled->GetCorrelationFunctionNP(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 567     gHistNP[1]->GetYaxis()->SetTitleOffset(1.5);
 568     gHistNP[1]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 569     gHistNP[1]->SetTitle(histoTitle.Data());
 570     cNP[1] = new TCanvas("cNP1","",100,100,600,500);
 571     cNP[1]->SetFillColor(10);
 572     cNP[1]->SetHighLightColor(10);
 573     gHistNP[1]->DrawCopy("surf1fb");
 574     gPad->SetTheta(30); // default is 30
 575     //gPad->SetPhi(130); // default is 30
 576     gPad->SetPhi(-60); // default is 30
 577     gPad->Update();
 578     pngName = "DeltaPhiDeltaEtaShuffled.Centrality";
 579     pngName += centralityArray[gCentrality-1];
 580     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 581     pngName += ".NegativePositive.png";
 582     //cNP[1]->SaveAs(pngName.Data());
 583   }
 584
 585   if(listBFMixed) {
 586     histoTitle = "(-+) mixed | Centrality: ";
 587     histoTitle += centralityArray[gCentrality-1];
 588     histoTitle += "%";
 589     if((psiMin == -0.5)&&(psiMax == 0.5))
 590       histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 591     else if((psiMin == 0.5)&&(psiMax == 1.5))
 592       histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 593     else if((psiMin == 1.5)&&(psiMax == 2.5))
 594       histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 595     else
 596       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 597
 598     gHistNP[2] = bMixed->GetCorrelationFunctionNP(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 599     gHistNP[2]->GetYaxis()->SetTitleOffset(1.5);
 600     gHistNP[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 601     gHistNP[2]->SetTitle(histoTitle.Data());
 602     cNP[2] = new TCanvas("cNP2","",100,200,600,500);
 603     cNP[2]->SetFillColor(10);
 604     cNP[2]->SetHighLightColor(10);
 605     gHistNP[2]->DrawCopy("surf1fb");
 606     gPad->SetTheta(30); // default is 30
 607     //gPad->SetPhi(130); // default is 30
 608     gPad->SetPhi(-60); // default is 30
 609     gPad->Update();
 610     pngName = "DeltaPhiDeltaEtaMixed.Centrality";
 611     pngName += centralityArray[gCentrality-1];
 612     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 613     pngName += ".NegativePositive.png";
 614     //cNP[2]->SaveAs(pngName.Data());
 615
 616     //Correlation function (-+)
 617     gHistNP[3] = dynamic_cast<TH2D *>(gHistNP[0]->Clone());
 618     gHistNP[3]->Divide(gHistNP[2]);
 619     gHistNP[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
 620     gHistNP[3]->GetZaxis()->SetTitle("C_{-+}(#Delta#eta,#Delta#varphi)");
 621     cNP[3] = new TCanvas("cNP3","",100,300,600,500);
 622     cNP[3]->SetFillColor(10);
 623     cNP[3]->SetHighLightColor(10);
 624     gHistNP[3]->DrawCopy("surf1fb");
 625     gPad->SetTheta(30); // default is 30
 626     //gPad->SetPhi(130); // default is 30
 627     gPad->SetPhi(-60); // default is 30
 628     gPad->Update();
 629     pngName = "CorrelationFunction.Centrality";
 630     pngName += centralityArray[gCentrality-1];
 631     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 632     pngName += ".NegativePositive.png";
 633     //cNP[3]->SaveAs(pngName.Data());
 634   }
 635
 636   //(++)
 637   histoTitle = "(++) | Centrality: ";
 638   histoTitle += centralityArray[gCentrality-1];
 639   histoTitle += "%";
 640   if((psiMin == -0.5)&&(psiMax == 0.5))
 641     histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 642   else if((psiMin == 0.5)&&(psiMax == 1.5))
 643     histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 644   else if((psiMin == 1.5)&&(psiMax == 2.5))
 645     histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 646   else
 647     histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 648
 649   gHistPP[0] = b->GetCorrelationFunctionPP(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 650   gHistPP[0]->GetYaxis()->SetTitleOffset(1.5);
 651   gHistPP[0]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 652   gHistPP[0]->SetTitle(histoTitle.Data());
 653   cPP[0] = new TCanvas("cPP0","",200,0,600,500);
 654   cPP[0]->SetFillColor(10);
 655   cPP[0]->SetHighLightColor(10);
 656   gHistPP[0]->DrawCopy("surf1fb");
 657   gPad->SetTheta(30); // default is 30
 658   //gPad->SetPhi(130); // default is 30
 659   gPad->SetPhi(-60); // default is 30
 660   gPad->Update();
 661   pngName = "DeltaPhiDeltaEta.Centrality";
 662   pngName += centralityArray[gCentrality-1];
 663   pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 664   pngName += ".PositivePositive.png";
 665   //cPP[0]->SaveAs(pngName.Data());
 666
 667   if(listBFShuffled) {
 668     histoTitle = "(++) shuffled | Centrality: ";
 669     histoTitle += centralityArray[gCentrality-1];
 670     histoTitle += "%";
 671     if((psiMin == -0.5)&&(psiMax == 0.5))
 672       histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 673     else if((psiMin == 0.5)&&(psiMax == 1.5))
 674       histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 675     else if((psiMin == 1.5)&&(psiMax == 2.5))
 676       histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 677     else
 678       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 679
 680     gHistPP[1] = bShuffled->GetCorrelationFunctionPP(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 681     gHistPP[1]->GetYaxis()->SetTitleOffset(1.5);
 682     gHistPP[1]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 683     gHistPP[1]->SetTitle(histoTitle.Data());
 684     cPP[1] = new TCanvas("cPP1","",200,100,600,500);
 685     cPP[1]->SetFillColor(10);
 686     cPP[1]->SetHighLightColor(10);
 687     gHistPP[1]->DrawCopy("surf1fb");
 688     gPad->SetTheta(30); // default is 30
 689     //gPad->SetPhi(130); // default is 30
 690     gPad->SetPhi(-60); // default is 30
 691     gPad->Update();
 692     pngName = "DeltaPhiDeltaEtaShuffled.Centrality";
 693     pngName += centralityArray[gCentrality-1];
 694     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 695     pngName += ".PositivePositive.png";
 696     //cPP[1]->SaveAs(pngName.Data());
 697   }
 698
 699   if(listBFMixed) {
 700     histoTitle = "(++) mixed | Centrality: ";
 701     histoTitle += centralityArray[gCentrality-1];
 702     histoTitle += "%";
 703     if((psiMin == -0.5)&&(psiMax == 0.5))
 704       histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 705     else if((psiMin == 0.5)&&(psiMax == 1.5))
 706       histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 707     else if((psiMin == 1.5)&&(psiMax == 2.5))
 708       histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 709     else
 710       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 711
 712     gHistPP[2] = bMixed->GetCorrelationFunctionPP(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 713     gHistPP[2]->GetYaxis()->SetTitleOffset(1.5);
 714     gHistPP[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 715     gHistPP[2]->SetTitle(histoTitle.Data());
 716     cPP[2] = new TCanvas("cPP2","",200,200,600,500);
 717     cPP[2]->SetFillColor(10);
 718     cPP[2]->SetHighLightColor(10);
 719     gHistPP[2]->DrawCopy("surf1fb");
 720     gPad->SetTheta(30); // default is 30
 721     //gPad->SetPhi(130); // default is 30
 722     gPad->SetPhi(-60); // default is 30
 723     gPad->Update();
 724     pngName = "DeltaPhiDeltaEtaMixed.Centrality";
 725     pngName += centralityArray[gCentrality-1];
 726     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 727     pngName += ".PositivePositive.png";
 728     //cPP[2]->SaveAs(pngName.Data());
 729
 730     //Correlation function (++)
 731     gHistPP[3] = dynamic_cast<TH2D *>(gHistPP[0]->Clone());
 732     gHistPP[3]->Divide(gHistPP[2]);
 733     gHistPP[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
 734     gHistPP[3]->GetZaxis()->SetTitle("C_{++}(#Delta#eta,#Delta#varphi)");
 735     cPP[3] = new TCanvas("cPP3","",200,300,600,500);
 736     cPP[3]->SetFillColor(10);
 737     cPP[3]->SetHighLightColor(10);
 738     gHistPP[3]->DrawCopy("surf1fb");
 739     gPad->SetTheta(30); // default is 30
 740     //gPad->SetPhi(130); // default is 30
 741     gPad->SetPhi(-60); // default is 30
 742     gPad->Update();
 743     pngName = "CorrelationFunction.Centrality";
 744     pngName += centralityArray[gCentrality-1];
 745     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 746     pngName += ".PositivePositive.png";
 747     //cPP[3]->SaveAs(pngName.Data());
 748   }
 749
 750   //(--)
 751   histoTitle = "(--) | Centrality: ";
 752   histoTitle += centralityArray[gCentrality-1];
 753   histoTitle += "%";
 754   if((psiMin == -0.5)&&(psiMax == 0.5))
 755     histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 756   else if((psiMin == 0.5)&&(psiMax == 1.5))
 757     histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 758   else if((psiMin == 1.5)&&(psiMax == 2.5))
 759     histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 760   else
 761     histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 762
 763   gHistNN[0] = b->GetCorrelationFunctionNN(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 764   gHistNN[0]->GetYaxis()->SetTitleOffset(1.5);
 765   gHistNN[0]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 766   gHistNN[0]->SetTitle(histoTitle.Data());
 767   cNN[0] = new TCanvas("cNN0","",300,0,600,500);
 768   cNN[0]->SetFillColor(10);
 769   cNN[0]->SetHighLightColor(10);
 770   gHistNN[0]->DrawCopy("surf1fb");
 771   gPad->SetTheta(30); // default is 30
 772   gPad->SetPhi(-60); // default is 30
 773   //gPad->SetPhi(-60); // default is 30
 774   gPad->Update();
 775   pngName = "DeltaPhiDeltaEta.Centrality";
 776   pngName += centralityArray[gCentrality-1];
 777   pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 778   pngName += ".NegativeNegative.png";
 779   //cNN[0]->SaveAs(pngName.Data());
 780
 781   if(listBFShuffled) {
 782     histoTitle = "(--) shuffled | Centrality: ";
 783     histoTitle += centralityArray[gCentrality-1];
 784     histoTitle += "%";
 785     if((psiMin == -0.5)&&(psiMax == 0.5))
 786       histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 787     else if((psiMin == 0.5)&&(psiMax == 1.5))
 788       histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 789     else if((psiMin == 1.5)&&(psiMax == 2.5))
 790       histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 791     else
 792       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 793
 794     gHistNN[1] = bShuffled->GetCorrelationFunctionNN(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 795     gHistNN[1]->GetYaxis()->SetTitleOffset(1.5);
 796     gHistNN[1]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 797     gHistNN[1]->SetTitle(histoTitle.Data());
 798     cNN[1] = new TCanvas("cNN1","",300,100,600,500);
 799     cNN[1]->SetFillColor(10);
 800     cNN[1]->SetHighLightColor(10);
 801     gHistNN[1]->DrawCopy("surf1fb");
 802     gPad->SetTheta(30); // default is 30
 803     //gPad->SetPhi(130); // default is 30
 804     gPad->SetPhi(-60); // default is 30
 805     gPad->Update();
 806     pngName = "DeltaPhiDeltaEtaShuffled.Centrality";
 807     pngName += centralityArray[gCentrality-1];
 808     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 809     pngName += ".NegativeNegative.png";
 810     //cNN[1]->SaveAs(pngName.Data());
 811   }
 812
 813   if(listBFMixed) {
 814     histoTitle = "(--) mixed | Centrality: ";
 815     histoTitle += centralityArray[gCentrality-1];
 816     histoTitle += "%";
 817     if((psiMin == -0.5)&&(psiMax == 0.5))
 818       histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 819     else if((psiMin == 0.5)&&(psiMax == 1.5))
 820       histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 821     else if((psiMin == 1.5)&&(psiMax == 2.5))
 822       histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 823     else
 824       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 825
 826     gHistNN[2] = bMixed->GetCorrelationFunctionNN(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 827     gHistNN[2]->GetYaxis()->SetTitleOffset(1.5);
 828     gHistNN[2]->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 829     gHistNN[2]->SetTitle(histoTitle.Data());
 830     cNN[2] = new TCanvas("cNN2","",300,200,600,500);
 831     cNN[2]->SetFillColor(10);
 832     cNN[2]->SetHighLightColor(10);
 833     gHistNN[2]->DrawCopy("surf1fb");
 834     gPad->SetTheta(30); // default is 30
 835     //gPad->SetPhi(130); // default is 30
 836     gPad->SetPhi(-60); // default is 30
 837     gPad->Update();
 838     pngName = "DeltaPhiDeltaEtaMixed.Centrality";
 839     pngName += centralityArray[gCentrality-1];
 840     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 841     pngName += ".NegativeNegative.png";
 842     //cNN[2]->SaveAs(pngName.Data());
 843
 844     //Correlation function (--)
 845     gHistNN[3] = dynamic_cast<TH2D *>(gHistNN[0]->Clone());
 846     gHistNN[3]->Divide(gHistNN[2]);
 847     gHistNN[3]->GetXaxis()->SetRangeUser(-1.5,1.5);
 848     gHistNN[3]->GetZaxis()->SetTitle("C_{--}(#Delta#eta,#Delta#varphi)");
 849     cNN[3] = new TCanvas("cNN3","",300,300,600,500);
 850     cNN[3]->SetFillColor(10);
 851     cNN[3]->SetHighLightColor(10);
 852     gHistNN[3]->DrawCopy("surf1fb");
 853     gPad->SetTheta(30); // default is 30
 854     //gPad->SetPhi(130); // default is 30
 855     gPad->SetPhi(-60); // default is 30
 856     gPad->Update();
 857     pngName = "CorrelationFunction.Centrality";
 858     pngName += centralityArray[gCentrality-1];
 859     pngName += ".Psi"; pngName += psiMin; pngName += "To"; pngName += psiMax;
 860     pngName += ".NegativeNegative.png";
 861     //cNN[3]->SaveAs(pngName.Data());
 862   }
 863
 864   //Write to output file
 865   TString newFileName = "correlationFunction.Centrality";
 866   newFileName += gCentrality; newFileName += ".Psi";
 867   if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
 868   else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
 869   else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
 870   else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
 871   else newFileName += "All.PttFrom";
 872   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
 873   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
 874   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
 875   newFileName += Form("%.1f",ptAssociatedMax);
 876   newFileName += ".root";
 877   TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
 878   gHistPN[0]->SetName("gHistPNRaw"); gHistPN[0]->Write();
 879   gHistNP[0]->SetName("gHistNPRaw"); gHistNP[0]->Write();
 880   gHistPP[0]->SetName("gHistPPRaw"); gHistPP[0]->Write();
 881   gHistNN[0]->SetName("gHistNNRaw"); gHistNN[0]->Write();
 882   if(listBFShuffled) {
 883     gHistPN[1]->SetName("gHistPNShuffled"); gHistPN[1]->Write();
 884     gHistNP[1]->SetName("gHistNPShuffled"); gHistNP[1]->Write();
 885     gHistPP[1]->SetName("gHistPPShuffled"); gHistPP[1]->Write();
 886     gHistNN[1]->SetName("gHistNNShuffled"); gHistNN[1]->Write();
 887   }
 888   if(listBFMixed) {
 889     gHistPN[2]->SetName("gHistPNMixed"); gHistPN[2]->Write();
 890     gHistNP[2]->SetName("gHistNPMixed"); gHistNP[2]->Write();
 891     gHistPP[2]->SetName("gHistPPMixed"); gHistPP[2]->Write();
 892     gHistNN[2]->SetName("gHistNNMixed"); gHistNN[2]->Write();
 893
 894     gHistPN[3]->SetName("gHistPNCorrelationFunctions"); gHistPN[3]->Write();
 895     gHistNP[3]->SetName("gHistNPCorrelationFunctions"); gHistNP[3]->Write();
 896     gHistPP[3]->SetName("gHistPPCorrelationFunctions"); gHistPP[3]->Write();
 897     gHistNN[3]->SetName("gHistNNCorrelationFunctions"); gHistNN[3]->Write();
 898   }
 899   newFile->Close();
 900
 901   // some cleaning
 902   for(Int_t i = 0; i < 4; i++){
 903
 904     if(!listBFShuffled && i == 1) continue;
 905     if(!listBFMixed && (i == 2 || i == 3)) continue;
 906
 907     if(gHistPP[i]) delete gHistPP[i];
 908     if(gHistPN[i]) delete gHistPN[i];
 909     if(gHistNP[i]) delete gHistNP[i];
 910     if(gHistNN[i]) delete gHistNN[i];
 911
 912     if(cPN[i]) delete cPN[i];
 913     if(cNP[i]) delete cNP[i];
 914     if(cPP[i]) delete cPP[i];
 915     if(cNN[i]) delete cNN[i];
 916   }
 917
 918   delete hP;
 919   delete hN;
 920   delete hPP;
 921   delete hPN;
 922   delete hNP;
 923   delete hNN;
 924
 925   delete hPMixed;
 926   delete hNMixed;
 927   delete hPPMixed;
 928   delete hPNMixed;
 929   delete hNPMixed;
 930   delete hNNMixed;
 931
 932   delete hPShuffled;
 933   delete hNShuffled;
 934   delete hPPShuffled;
 935   delete hPNShuffled;
 936   delete hNPShuffled;
 937   delete hNNShuffled;
 938
 939 }
 940
 941 //____________________________________________________________//
 942 void drawCorrelationFunctions(const char* lhcPeriod = "LHC11h",
 943                               Int_t gTrainID = 208,
 944                               Int_t gCentrality = 1,
 945                               Double_t psiMin = -0.5, Double_t psiMax = 3.5,
 946                               Double_t ptTriggerMin = -1.,
 947                               Double_t ptTriggerMax = -1.,
 948                               Double_t ptAssociatedMin = -1.,
 949                               Double_t ptAssociatedMax = -1.) {
 950   //Macro that draws the charge dependent correlation functions
 951   //for each centrality bin for the different pT of trigger and
 952   //associated particles
 953   //Author: Panos.Christakoglou@nikhef.nl
 954   TGaxis::SetMaxDigits(3);
 955
 956   //Get the input file
 957   TString filename = "PbPb/"; filename += lhcPeriod;
 958   filename +="/Train"; filename += gTrainID;
 959   filename +="/Centrality"; filename += gCentrality;
 960   filename += "/correlationFunction.Centrality";
 961   filename += gCentrality; filename += ".Psi";
 962   if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
 963   else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
 964   else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
 965   else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
 966   else filename += "All.PttFrom";
 967   filename += Form("%.1f",ptTriggerMin); filename += "To";
 968   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
 969   filename += Form("%.1f",ptAssociatedMin); filename += "To";
 970   filename += Form("%.1f",ptAssociatedMax);
 971   filename += ".root";
 972
 973   //Open the file
 974   TFile *f = TFile::Open(filename.Data());
 975   if((!f)||(!f->IsOpen())) {
 976     Printf("The file %s is not found. Aborting...",filename);
 977     return listBF;
 978   }
 979   //f->ls();
 980
 981   //Latex
 982   TString centralityLatex = "Centrality: ";
 983   centralityLatex += centralityArray[gCentrality-1];
 984   centralityLatex += "%";
 985
 986   TString psiLatex;
 987   if((psiMin == -0.5)&&(psiMax == 0.5))
 988     psiLatex = " -7.5^{o} < #varphi - #Psi_{2} < 7.5^{o}";
 989   else if((psiMin == 0.5)&&(psiMax == 1.5))
 990     psiLatex = " 37.5^{o} < #varphi - #Psi_{2} < 52.5^{o}";
 991   else if((psiMin == 1.5)&&(psiMax == 2.5))
 992     psiLatex = " 82.5^{o} < #varphi - #Psi_{2} < 97.5^{o}";
 993   else
 994     psiLatex = " 0^{o} < #varphi - #Psi_{2} < 180^{o}";
 995
 996   TString pttLatex = Form("%.1f",ptTriggerMin);
 997   pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
 998   pttLatex += " GeV/c";
 999
1000   TString ptaLatex = Form("%.1f",ptAssociatedMin);
1001   ptaLatex += " < p_{T,assoc} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1002   ptaLatex += " GeV/c";
1003
1004   TLatex *latexInfo1 = new TLatex();
1005   latexInfo1->SetNDC();
1006   latexInfo1->SetTextSize(0.045);
1007   latexInfo1->SetTextColor(1);
1008
1009   TString pngName;
1010
1011   //============================================================//
1012   //Get the +- correlation function
1013   TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
1014   gHistPN->SetStats(kFALSE);
1015   gHistPN->SetTitle("");
1016   gHistPN->GetXaxis()->SetRangeUser(-1.45,1.45);
1017   gHistPN->GetXaxis()->CenterTitle();
1018   gHistPN->GetXaxis()->SetTitleOffset(1.2);
1019   gHistPN->GetYaxis()->CenterTitle();
1020   gHistPN->GetYaxis()->SetTitleOffset(1.2);
1021   gHistPN->GetZaxis()->SetTitleOffset(1.2);
1022   TCanvas *cPN = new TCanvas("cPN","",0,0,600,500);
1023   cPN->SetFillColor(10); cPN->SetHighLightColor(10);
1024   cPN->SetLeftMargin(0.15);
1025   gHistPN->DrawCopy("surf1fb");
1026   gPad->SetTheta(30); // default is 30
1027   gPad->SetPhi(-60); // default is 30
1028   gPad->Update();
1029
1030   latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
1031   //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
1032   latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
1033   latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
1034
1035   pngName = "CorrelationFunction.Centrality";
1036   pngName += centralityArray[gCentrality-1];
1037   pngName += ".Psi";
1038   if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
1039   else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
1040   else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
1041   else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
1042   else pngName += "All.PttFrom";
1043   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1044   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1045   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1046   pngName += Form("%.1f",ptAssociatedMax);
1047   pngName += ".PositiveNegative.png";
1048   cPN->SaveAs(pngName.Data());
1049   fitCorrelationFunctions(gCentrality, psiMin, psiMax,
1050                           ptTriggerMin,ptTriggerMax,
1051                           ptAssociatedMin, ptAssociatedMax,gHistPN);
1052   //============================================================//
1053   //Get the -+ correlation function
1054   TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
1055   gHistNP->SetStats(kFALSE);
1056   gHistNP->SetTitle("");
1057   gHistNP->GetXaxis()->SetRangeUser(-1.45,1.45);
1058   gHistNP->GetXaxis()->CenterTitle();
1059   gHistNP->GetXaxis()->SetTitleOffset(1.2);
1060   gHistNP->GetYaxis()->CenterTitle();
1061   gHistNP->GetYaxis()->SetTitleOffset(1.2);
1062   gHistNP->GetZaxis()->SetTitleOffset(1.2);
1063   TCanvas *cNP = new TCanvas("cNP","",50,50,600,500);
1064   cNP->SetFillColor(10); cNP->SetHighLightColor(10);
1065   cNP->SetLeftMargin(0.15);
1066   gHistNP->DrawCopy("surf1fb");
1067   gPad->SetTheta(30); // default is 30
1068   gPad->SetPhi(-60); // default is 30
1069   gPad->Update();
1070
1071   latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
1072   //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
1073   latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
1074   latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
1075
1076   pngName = "CorrelationFunction.Centrality";
1077   pngName += centralityArray[gCentrality-1];
1078   pngName += ".Psi";
1079   if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
1080   else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
1081   else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
1082   else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
1083   else pngName += "All.PttFrom";
1084   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1085   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1086   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1087   pngName += Form("%.1f",ptAssociatedMax);
1088   pngName += ".NegativePositive.png";
1089   cNP->SaveAs(pngName.Data());
1090
1091   fitCorrelationFunctions(gCentrality, psiMin, psiMax,
1092                           ptTriggerMin,ptTriggerMax,
1093                           ptAssociatedMin, ptAssociatedMax,gHistNP);
1094   //============================================================//
1095   //Get the ++ correlation function
1096   TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
1097   gHistPP->SetStats(kFALSE);
1098   gHistPP->SetTitle("");
1099   gHistPP->GetXaxis()->SetRangeUser(-1.45,1.45);
1100   gHistPP->GetXaxis()->CenterTitle();
1101   gHistPP->GetXaxis()->SetTitleOffset(1.2);
1102   gHistPP->GetYaxis()->CenterTitle();
1103   gHistPP->GetYaxis()->SetTitleOffset(1.2);
1104   gHistPP->GetZaxis()->SetTitleOffset(1.2);
1105   TCanvas *cPP = new TCanvas("cPP","",100,100,600,500);
1106   cPP->SetFillColor(10); cPP->SetHighLightColor(10);
1107   cPP->SetLeftMargin(0.15);
1108   gHistPP->DrawCopy("surf1fb");
1109   gPad->SetTheta(30); // default is 30
1110   gPad->SetPhi(-60); // default is 30
1111   gPad->Update();
1112
1113   latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
1114   //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
1115   latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
1116   latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
1117
1118   pngName = "CorrelationFunction.Centrality";
1119   pngName += centralityArray[gCentrality-1];
1120   pngName += ".Psi";
1121   if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
1122   else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
1123   else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
1124   else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
1125   else pngName += "All.PttFrom";
1126   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1127   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1128   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1129   pngName += Form("%.1f",ptAssociatedMax);
1130   pngName += ".PositivePositive.png";
1131   cPP->SaveAs(pngName.Data());
1132
1133   fitCorrelationFunctions(gCentrality, psiMin, psiMax,
1134                           ptTriggerMin,ptTriggerMax,
1135                           ptAssociatedMin, ptAssociatedMax,gHistPP);
1136   //============================================================//
1137   //Get the -- correlation function
1138   TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
1139   gHistNN->SetStats(kFALSE);
1140   gHistNN->SetTitle("");
1141   gHistNN->GetXaxis()->SetRangeUser(-1.45,1.45);
1142   gHistNN->GetXaxis()->CenterTitle();
1143   gHistNN->GetXaxis()->SetTitleOffset(1.2);
1144   gHistNN->GetYaxis()->CenterTitle();
1145   gHistNN->GetYaxis()->SetTitleOffset(1.2);
1146   gHistNN->GetZaxis()->SetTitleOffset(1.2);
1147   TCanvas *cNN = new TCanvas("cNN","",150,150,600,500);
1148   cNN->SetFillColor(10); cNN->SetHighLightColor(10);
1149   cNN->SetLeftMargin(0.15);
1150   gHistNN->DrawCopy("surf1fb");
1151   gPad->SetTheta(30); // default is 30
1152   gPad->SetPhi(-60); // default is 30
1153   gPad->Update();
1154
1155   latexInfo1->DrawLatex(0.44,0.88,centralityLatex.Data());
1156   //latexInfo1->DrawLatex(0.44,0.82,psiLatex.Data());
1157   latexInfo1->DrawLatex(0.44,0.82,pttLatex.Data());
1158   latexInfo1->DrawLatex(0.44,0.76,ptaLatex.Data());
1159
1160   pngName = "CorrelationFunction.Centrality";
1161   pngName += centralityArray[gCentrality-1];
1162   pngName += ".Psi";
1163   if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
1164   else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
1165   else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
1166   else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.Ptt";
1167   else pngName += "All.PttFrom";
1168   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1169   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1170   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1171   pngName += Form("%.1f",ptAssociatedMax);
1172   pngName += ".NegativeNegative.png";
1173   cNN->SaveAs(pngName.Data());
1174
1175   fitCorrelationFunctions(gCentrality, psiMin, psiMax,
1176                           ptTriggerMin,ptTriggerMax,
1177                           ptAssociatedMin, ptAssociatedMax,gHistNN);
1178 }
1179
1180 //____________________________________________________________//
1181 void fitCorrelationFunctions(Int_t gCentrality = 1,
1182                              Double_t psiMin = -0.5, Double_t psiMax = 3.5,
1183                              Double_t ptTriggerMin = -1.,
1184                              Double_t ptTriggerMax = -1.,
1185                              Double_t ptAssociatedMin = -1.,
1186                              Double_t ptAssociatedMax = -1.,
1187                              TH2D *gHist) {
1188
1189   cout<<"FITTING FUNCTION"<<endl;
1190
1191   //near side peak: [1]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4]))
1192   //away side ridge: [5]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((y-TMath::Pi())/[6]),2)),[7]))
1193   //longitudinal ridge: [8]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[9]),2)),[10]))
1194   //wing structures: [11]*TMath::Power(x,2)
1195   //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))
1196   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.);
1197   gFitFunction->SetName("gFitFunction");
1198   //Normalization
1199   gFitFunction->SetParName(0,"N1"); gFitFunction->SetParameter(0,1.0);//gFitFunction->SetParLimits(0,0.0,100);
1200   //near side peak
1201   gFitFunction->SetParName(1,"N_{near side}"); gFitFunction->SetParameter(1,0.3);//gFitFunction->SetParLimits(1,0.0,100);//gFitFunction->FixParameter(1,0.0);
1202   gFitFunction->SetParName(2,"Sigma_{near side}(delta eta)"); gFitFunction->SetParameter(2,0.3);gFitFunction->SetParLimits(2,0.01,0.7);//gFitFunction->FixParameter(2,0.0);
1203   gFitFunction->SetParName(3,"Sigma_{near side}(delta phi)"); gFitFunction->SetParameter(3,0.3);gFitFunction->SetParLimits(3,0.01,1.7);//gFitFunction->FixParameter(3,0.0);
1204   gFitFunction->SetParName(4,"Exponent_{near side}"); gFitFunction->SetParameter(4,1.1);//gFitFunction->SetParLimits(4,0.1,10.0);
1205   gFitFunction->FixParameter(4,1.0);
1206   //away side ridge
1207   gFitFunction->SetParName(5,"N_{away side}"); gFitFunction->SetParameter(5,1.0);//gFitFunction->SetParLimits(5,0.0,100);gFitFunction->FixParameter(5,0.0);
1208   gFitFunction->SetParName(6,"Sigma_{away side}(delta phi)"); gFitFunction->SetParameter(6,0.5);//gFitFunction->SetParLimits(6,0.01,100.0);gFitFunction->FixParameter(6,0.0);
1209   gFitFunction->SetParName(7,"Exponent_{away side}"); gFitFunction->SetParameter(7,1.0);//gFitFunction->SetParLimits(7,0.1,10.0);
1210   gFitFunction->FixParameter(7,1.0);
1211   //longitudinal ridge
1212   gFitFunction->SetParName(8,"N_{long. ridge}"); gFitFunction->SetParameter(8,0.05);//gFitFunction->FixParameter(8,1.0);
1213   gFitFunction->SetParName(9,"Sigma_{long. ridge}(delta eta)"); gFitFunction->SetParameter(9,0.6);gFitFunction->SetParLimits(9,0.1,10.0);//gFitFunction->FixParameter(9,0.0);
1214   gFitFunction->SetParName(10,"Exponent_{long. ridge}"); gFitFunction->SetParameter(10,1.0);gFitFunction->FixParameter(10,1.0);
1215   //wing structures
1216   gFitFunction->SetParName(11,"N_{wing}"); gFitFunction->SetParameter(11,0.01);gFitFunction->FixParameter(11,0.0);
1217   //flow contribution
1218   gFitFunction->SetParName(12,"N_{flow}"); gFitFunction->SetParameter(12,0.2);//gFitFunction->SetParLimits(12,0.0,10.0);//gFitFunction->FixParameter(12,0.0);
1219   gFitFunction->SetParName(13,"V1"); gFitFunction->SetParameter(13,0.005);//gFitFunction->SetParLimits(13,0.0,10.0);//gFitFunction->FixParameter(13,0.0);
1220   gFitFunction->SetParName(14,"V2"); gFitFunction->SetParameter(14,0.1);//gFitFunction->SetParLimits(14,0.0,10.0);//gFitFunction->FixParameter(14,0.0);
1221   gFitFunction->SetParName(15,"V3"); gFitFunction->SetParameter(15,0.05);//gFitFunction->SetParLimits(15,0.0,10.0);//gFitFunction->FixParameter(15,0.0);
1222   gFitFunction->SetParName(16,"V4"); gFitFunction->SetParameter(16,0.005);//gFitFunction->SetParLimits(16,0.0,10.0);//gFitFunction->FixParameter(16,0.0);
1223   gFitFunction->SetParName(17,"Offset"); gFitFunction->SetParameter(17,0.005);gFitFunction->SetParLimits(16,0.0,0.7);//gFitFunction->FixParameter(17,0.0);
1224
1225   //Fitting the correlation function
1226   gHist->Fit("gFitFunction","nm");
1227
1228   //Cloning the histogram
1229   TString histoName = gHist->GetName(); histoName += "Fit";
1230   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());
1231   TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
1232   gHistResidual->SetName("gHistResidual");
1233   gHistResidual->Sumw2();
1234
1235   for(Int_t iBinDeltaEta = 1; iBinDeltaEta <= gHist->GetNbinsX(); iBinDeltaEta++) {
1236     for(Int_t iBinDeltaPhi = 1; iBinDeltaPhi <= gHist->GetNbinsY(); iBinDeltaPhi++) {
1237       gHistFit->SetBinContent(iBinDeltaEta,iBinDeltaPhi,gFitFunction->Eval(gHist->GetXaxis()->GetBinCenter(iBinDeltaEta),gHist->GetYaxis()->GetBinCenter(iBinDeltaPhi)));
1238     }
1239   }
1240   gHistResidual->Add(gHistFit,-1);
1241
1242   //Write to output file
1243   TString newFileName = "correlationFunctionFit";
1244   if(histoName.Contains("PN")) newFileName += "PN";
1245   else if(histoName.Contains("NP")) newFileName += "NP";
1246   else if(histoName.Contains("PP")) newFileName += "PP";
1247   else if(histoName.Contains("NN")) newFileName += "NN";
1248   newFileName += ".Centrality";
1249   newFileName += gCentrality; newFileName += ".Psi";
1250   if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
1251   else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
1252   else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
1253   else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
1254   else newFileName += "All.PttFrom";
1255   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
1256   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
1257   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
1258   newFileName += Form("%.1f",ptAssociatedMax);
1259   newFileName += ".root";
1260   TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
1261   gHist->Write();
1262   gHistFit->Write();
1263   gHistResidual->Write();
1264   gFitFunction->Write();
1265   newFile->Close();
1266
1267
1268 }
1269
1270
1271
1272 // //____________________________________________________________//
1273 // void fitCorrelationFunctions(Int_t gCentrality = 1,
1274 //                           Double_t psiMin = -0.5, Double_t psiMax = 3.5,
1275 //                           Double_t ptTriggerMin = -1.,
1276 //                           Double_t ptTriggerMax = -1.,
1277 //                           Double_t ptAssociatedMin = -1.,
1278 //                           Double_t ptAssociatedMax = -1.,
1279 //                           TH2D *gHist) {
1280
1281 //   cout<<"FITTING FUNCTION (HOUSTON)"<<endl;
1282
1283 //   // Fit Function
1284 //    //x axis = delta_eta
1285 //    //y axis = delta_phi
1286 //                                                                                                                                                                                         //  [9]*exp(-1*pow(((x/[10])^2 + (y/[10])^2),0.5)) Hyper expo
1287 //   TF2 *fit1 = new TF2("fit1","[0] + [1]*cos(y) + [2]*cos(2*y) + [3]*cos(3*y) + [4]*cos(4*y) +[5]*cos(5*y)+ [6]*exp(-0.5*pow(((x/[7])^2 + (y/[8])^2),[11])) + [6]*exp(-0.5*pow(((x/[7])^2 + ((y-6.283)/[8])^2),[11]))+ [9]*exp(-1*((x/[10])^2 + (y/[10])^2)) ",-2.0,2.0,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
1288
1289
1290
1291 //   Double_t Parameters[] = {0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.3,0.3,1.0,0.1,0.1};
1292
1293 //   fit1->SetParameters(Parameters);  // input pars from macro arguments
1294
1295 //   fit1->SetParName(0,"offset");
1296 //   fit1->SetParName(1,"v1");
1297 //   fit1->SetParName(2,"v2");
1298 //   fit1->SetParName(3,"v3");
1299 //   fit1->SetParName(4,"v4");
1300 //   fit1->SetParName(5,"v5");
1301 //   fit1->SetParName(6,"Ridgeamp");
1302 //   fit1->SetParName(7,"Ridgesigx");
1303 //   fit1->SetParName(8,"Ridgesigy");
1304 //   fit1->SetParName(9,"Expoamp");
1305 //   fit1->SetParName(10,"Exposig");
1306 //   fit1->SetParName(11,"Gausspara");
1307
1308
1309 //   //Fit Parameter Ranges
1310 //   fit1->SetParLimits(0,-2.0,2.0);   //offset
1311 //   fit1->SetParLimits(1,-1.0,0.1);   //v1
1312 //   fit1->SetParLimits(2,-1.6,0.9);   //v2
1313 //   fit1->SetParLimits(3,0.0,0.5);    //v3
1314 //   fit1->SetParLimits(4,0.0,0.9);    //v4
1315 //   fit1->SetParLimits(5,0.0,0.9);    //v5
1316 //   fit1->SetParLimits(6,0.0,1.5);    //Ridgeamp
1317 //   fit1->SetParLimits(7,0.1,3.0);    //Ridgesigx
1318 //   fit1->SetParLimits(8,0.1,2.0);    //Ridgesigy
1319 //   fit1->SetParLimits(9,0.0,6.0);      //Expoamp
1320 //   fit1->SetParLimits(10,0.05,0.5); //Exposig
1321 //   fit1->SetParLimits(11,0.0,2.0);    //Gausspara
1322
1323 //   //Fitting the correlation function
1324 //   gHist->Fit("fit1","nm");
1325
1326 //   //Cloning the histogram
1327 //   TString histoName = gHist->GetName(); histoName += "Fit";
1328 //   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());
1329 //   TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
1330 //   gHistResidual->SetName("gHistResidual");
1331 //   gHistResidual->Sumw2();
1332
1333 //   for(Int_t iBinDeltaEta = 1; iBinDeltaEta <= gHist->GetNbinsX(); iBinDeltaEta++) {
1334 //     for(Int_t iBinDeltaPhi = 1; iBinDeltaPhi <= gHist->GetNbinsY(); iBinDeltaPhi++) {
1335 //       gHistFit->SetBinContent(iBinDeltaEta,iBinDeltaPhi,fit1->Eval(gHist->GetXaxis()->GetBinCenter(iBinDeltaEta),gHist->GetYaxis()->GetBinCenter(iBinDeltaPhi)));
1336 //     }
1337 //   }
1338 //   gHistResidual->Add(gHistFit,-1);
1339
1340 //   //Write to output file
1341 //   TString newFileName = "correlationFunctionFit";
1342 //   if(histoName.Contains("PN")) newFileName += "PN";
1343 //   else if(histoName.Contains("NP")) newFileName += "NP";
1344 //   else if(histoName.Contains("PP")) newFileName += "PP";
1345 //   else if(histoName.Contains("NN")) newFileName += "NN";
1346 //   newFileName += ".Centrality";
1347 //   newFileName += gCentrality; newFileName += ".Psi";
1348 //   if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
1349 //   else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
1350 //   else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
1351 //   else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
1352 //   else newFileName += "All.PttFrom";
1353 //   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
1354 //   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
1355 //   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
1356 //   newFileName += Form("%.1f",ptAssociatedMax);
1357 //   newFileName += ".root";
1358 //   TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
1359 //   gHist->Write();
1360 //   gHistFit->Write();
1361 //   gHistResidual->Write();
1362 //   fit1->Write();
1363 //   newFile->Close();
1364
1365
1366 // }