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