PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C

   1 const Int_t numberOfCentralityBins = 12;
   2 TString centralityArray[numberOfCentralityBins] = {"0-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-100","0-1","1-2","2-3"};
   3
   4
   5 const Int_t gRebin = 1;
   6 void drawBalanceFunction2DPsi(const char* filename = "AnalysisResultsPsi.root",
   7                               Int_t gCentrality = 1,
   8                               Int_t gBit = -1,
   9                               const char* gCentralityEstimator = 0x0,
  10                               Bool_t kShowShuffled = kFALSE,
  11                               Bool_t kShowMixed = kTRUE,
  12                               Double_t psiMin = -0.5, Double_t psiMax = 0.5,
  13                               Double_t ptTriggerMin = -1.,
  14                               Double_t ptTriggerMax = -1.,
  15                               Double_t ptAssociatedMin = -1.,
  16                               Double_t ptAssociatedMax = -1.,
  17                               Bool_t k2pMethod = kFALSE,
  18                               TString eventClass = "EventPlane") //Can be "EventPlane", "Centrality", "Multiplicity"
  19 {
  20   //Macro that draws the BF distributions for each centrality bin
  21   //for reaction plane dependent analysis
  22   //Author: Panos.Christakoglou@nikhef.nl
  23   //Load the PWG2ebye library
  24   gSystem->Load("libANALYSIS.so");
  25   gSystem->Load("libANALYSISalice.so");
  26   gSystem->Load("libEventMixing.so");
  27   gSystem->Load("libCORRFW.so");
  28   gSystem->Load("libPWGTools.so");
  29   gSystem->Load("libPWGCFebye.so");
  30
  31   gROOT->LoadMacro("~/SetPlotStyle.C");
  32   SetPlotStyle();
  33   gStyle->SetPalette(1,0);
  34
  35   //Prepare the objects and return them
  36   TList *listBF = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,0);
  37   TList *listBFShuffled = NULL;
  38   if(kShowShuffled) listBFShuffled = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,1);
  39   TList *listBFMixed = NULL;
  40   if(kShowMixed) listBFMixed = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,2);
  41   if(!listBF) {
  42     Printf("The TList object was not created");
  43     return;
  44   }
  45   else
  46     draw(listBF,listBFShuffled,listBFMixed,gCentrality,gCentralityEstimator,
  47          psiMin,psiMax,
  48          ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,
  49          k2pMethod,eventClass);
  50 }
  51
  52 //______________________________________________________//
  53 TList *GetListOfObjects(const char* filename,
  54                         Int_t gCentrality,
  55                         Int_t gBit,
  56                         const char *gCentralityEstimator,
  57                         Int_t kData = 1) {
  58   //Get the TList objects (QA, bf, bf shuffled)
  59   TList *listBF = 0x0;
  60
  61   //Open the file
  62   TFile *f = TFile::Open(filename,"UPDATE");
  63   if((!f)||(!f->IsOpen())) {
  64     Printf("The file %s is not found. Aborting...",filename);
  65     return listBF;
  66   }
  67   //f->ls();
  68
  69   TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
  70   if(!dir) {
  71     Printf("The TDirectoryFile is not found. Aborting...",filename);
  72     return listBF;
  73   }
  74   //dir->ls();
  75
  76   TString listBFName;
  77   if(kData == 0) {
  78     //cout<<"no shuffling - no mixing"<<endl;
  79     listBFName = "listBFPsi_";
  80   }
  81   else if(kData == 1) {
  82     //cout<<"shuffling - no mixing"<<endl;
  83     listBFName = "listBFPsiShuffled_";
  84   }
  85   else if(kData == 2) {
  86     //cout<<"no shuffling - mixing"<<endl;
  87     listBFName = "listBFPsiMixed_";
  88   }
  89   listBFName += centralityArray[gCentrality-1];
  90   if(gBit > -1) {
  91     listBFName += "_Bit"; listBFName += gBit; }
  92   if(gCentralityEstimator) {
  93     listBFName += "_"; listBFName += gCentralityEstimator;}
  94
  95   // histograms were already retrieved (in first iteration)
  96   if(dir->Get(Form("%s_histograms",listBFName.Data()))){
  97     //cout<<"second iteration"<<endl;
  98     listBF = dynamic_cast<TList *>(dir->Get(Form("%s_histograms",listBFName.Data())));
  99   }
 100
 101   // histograms were not yet retrieved (this is the first iteration)
 102   else{
 103
 104     listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
 105     cout<<"======================================================="<<endl;
 106     cout<<"List name: "<<listBF->GetName()<<endl;
 107     //listBF->ls();
 108
 109     //Get the histograms
 110     TString histoName;
 111     if(kData == 0)
 112       histoName = "fHistP";
 113     else if(kData == 1)
 114       histoName = "fHistP_shuffle";
 115     else if(kData == 2)
 116       histoName = "fHistP";
 117     if(gCentralityEstimator) histoName += gCentralityEstimator;
 118     AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 119     if(!fHistP) {
 120       Printf("fHistP %s not found!!!",histoName.Data());
 121       break;
 122     }
 123     fHistP->FillParent(); fHistP->DeleteContainers();
 124
 125     if(kData == 0)
 126       histoName = "fHistN";
 127     if(kData == 1)
 128       histoName = "fHistN_shuffle";
 129     if(kData == 2)
 130       histoName = "fHistN";
 131     if(gCentralityEstimator) histoName += gCentralityEstimator;
 132     AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 133     if(!fHistN) {
 134       Printf("fHistN %s not found!!!",histoName.Data());
 135       break;
 136     }
 137     fHistN->FillParent(); fHistN->DeleteContainers();
 138
 139     if(kData == 0)
 140       histoName = "fHistPN";
 141     if(kData == 1)
 142       histoName = "fHistPN_shuffle";
 143     if(kData == 2)
 144       histoName = "fHistPN";
 145     if(gCentralityEstimator) histoName += gCentralityEstimator;
 146     AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 147     if(!fHistPN) {
 148       Printf("fHistPN %s not found!!!",histoName.Data());
 149       break;
 150     }
 151     fHistPN->FillParent(); fHistPN->DeleteContainers();
 152
 153     if(kData == 0)
 154       histoName = "fHistNP";
 155     if(kData == 1)
 156       histoName = "fHistNP_shuffle";
 157     if(kData == 2)
 158       histoName = "fHistNP";
 159     if(gCentralityEstimator) histoName += gCentralityEstimator;
 160     AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 161     if(!fHistNP) {
 162       Printf("fHistNP %s not found!!!",histoName.Data());
 163       break;
 164     }
 165     fHistNP->FillParent(); fHistNP->DeleteContainers();
 166
 167     if(kData == 0)
 168       histoName = "fHistPP";
 169     if(kData == 1)
 170       histoName = "fHistPP_shuffle";
 171     if(kData == 2)
 172       histoName = "fHistPP";
 173     if(gCentralityEstimator) histoName += gCentralityEstimator;
 174     AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 175     if(!fHistPP) {
 176       Printf("fHistPP %s not found!!!",histoName.Data());
 177       break;
 178     }
 179     fHistPP->FillParent(); fHistPP->DeleteContainers();
 180
 181     if(kData == 0)
 182       histoName = "fHistNN";
 183     if(kData == 1)
 184       histoName = "fHistNN_shuffle";
 185     if(kData == 2)
 186       histoName = "fHistNN";
 187     if(gCentralityEstimator) histoName += gCentralityEstimator;
 188     AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 189     if(!fHistNN) {
 190       Printf("fHistNN %s not found!!!",histoName.Data());
 191       break;
 192     }
 193     fHistNN->FillParent(); fHistNN->DeleteContainers();
 194
 195     dir->cd();
 196     listBF->Write(Form("%s_histograms",listBFName.Data()), TObject::kSingleKey);
 197
 198   }// first iteration
 199
 200   f->Close();
 201
 202   return listBF;
 203 }
 204
 205 //______________________________________________________//
 206 void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
 207           Int_t gCentrality, const char* gCentralityEstimator,
 208           Double_t psiMin, Double_t psiMax,
 209           Double_t ptTriggerMin, Double_t ptTriggerMax,
 210           Double_t ptAssociatedMin, Double_t ptAssociatedMax,
 211           Bool_t k2pMethod = kFALSE, TString eventClass) {
 212   //balance function
 213   AliTHn *hP = NULL;
 214   AliTHn *hN = NULL;
 215   AliTHn *hPN = NULL;
 216   AliTHn *hNP = NULL;
 217   AliTHn *hPP = NULL;
 218   AliTHn *hNN = NULL;
 219   //listBF->ls();
 220   //Printf("=================");
 221   TString histoName = "fHistP";
 222   if(gCentralityEstimator) histoName += gCentralityEstimator;
 223   hP = (AliTHn*) listBF->FindObject(histoName.Data());
 224   hP->SetName("gHistP");
 225   histoName = "fHistN";
 226   if(gCentralityEstimator) histoName += gCentralityEstimator;
 227   hN = (AliTHn*) listBF->FindObject(histoName.Data());
 228   hN->SetName("gHistN");
 229   histoName = "fHistPN";
 230   if(gCentralityEstimator) histoName += gCentralityEstimator;
 231   hPN = (AliTHn*) listBF->FindObject(histoName.Data());
 232   hPN->SetName("gHistPN");
 233   histoName = "fHistNP";
 234   if(gCentralityEstimator) histoName += gCentralityEstimator;
 235   hNP = (AliTHn*) listBF->FindObject(histoName.Data());
 236   hNP->SetName("gHistNP");
 237   histoName = "fHistPP";
 238   if(gCentralityEstimator) histoName += gCentralityEstimator;
 239   hPP = (AliTHn*) listBF->FindObject(histoName.Data());
 240   hPP->SetName("gHistPP");
 241   histoName = "fHistNN";
 242   if(gCentralityEstimator) histoName += gCentralityEstimator;
 243   hNN = (AliTHn*) listBF->FindObject(histoName.Data());
 244   hNN->SetName("gHistNN");
 245
 246   AliBalancePsi *b = new AliBalancePsi();
 247   b->SetEventClass(eventClass);
 248   b->SetHistNp(hP);
 249   b->SetHistNn(hN);
 250   b->SetHistNpn(hPN);
 251   b->SetHistNnp(hNP);
 252   b->SetHistNpp(hPP);
 253   b->SetHistNnn(hNN);
 254
 255   //balance function shuffling
 256   AliTHn *hPShuffled = NULL;
 257   AliTHn *hNShuffled = NULL;
 258   AliTHn *hPNShuffled = NULL;
 259   AliTHn *hNPShuffled = NULL;
 260   AliTHn *hPPShuffled = NULL;
 261   AliTHn *hNNShuffled = NULL;
 262   if(listBFShuffled) {
 263     //listBFShuffled->ls();
 264     histoName = "fHistP_shuffle";
 265     if(gCentralityEstimator) histoName += gCentralityEstimator;
 266     hPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
 267     hPShuffled->SetName("gHistPShuffled");
 268     histoName = "fHistN_shuffle";
 269     if(gCentralityEstimator) histoName += gCentralityEstimator;
 270     hNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
 271     hNShuffled->SetName("gHistNShuffled");
 272     histoName = "fHistPN_shuffle";
 273     if(gCentralityEstimator) histoName += gCentralityEstimator;
 274     hPNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
 275     hPNShuffled->SetName("gHistPNShuffled");
 276     histoName = "fHistNP_shuffle";
 277     if(gCentralityEstimator) histoName += gCentralityEstimator;
 278     hNPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
 279     hNPShuffled->SetName("gHistNPShuffled");
 280     histoName = "fHistPP_shuffle";
 281     if(gCentralityEstimator) histoName += gCentralityEstimator;
 282     hPPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
 283     hPPShuffled->SetName("gHistPPShuffled");
 284     histoName = "fHistNN_shuffle";
 285     if(gCentralityEstimator) histoName += gCentralityEstimator;
 286     hNNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
 287     hNNShuffled->SetName("gHistNNShuffled");
 288
 289     AliBalancePsi *bShuffled = new AliBalancePsi();
 290     bShuffled->SetEventClass(eventClass);
 291     bShuffled->SetHistNp(hPShuffled);
 292     bShuffled->SetHistNn(hNShuffled);
 293     bShuffled->SetHistNpn(hPNShuffled);
 294     bShuffled->SetHistNnp(hNPShuffled);
 295     bShuffled->SetHistNpp(hPPShuffled);
 296     bShuffled->SetHistNnn(hNNShuffled);
 297   }
 298
 299   //balance function mixing
 300   AliTHn *hPMixed = NULL;
 301   AliTHn *hNMixed = NULL;
 302   AliTHn *hPNMixed = NULL;
 303   AliTHn *hNPMixed = NULL;
 304   AliTHn *hPPMixed = NULL;
 305   AliTHn *hNNMixed = NULL;
 306
 307   if(listBFMixed) {
 308     //listBFMixed->ls();
 309     histoName = "fHistP";
 310     if(gCentralityEstimator) histoName += gCentralityEstimator;
 311     hPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
 312     hPMixed->SetName("gHistPMixed");
 313     histoName = "fHistN";
 314     if(gCentralityEstimator) histoName += gCentralityEstimator;
 315     hNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
 316     hNMixed->SetName("gHistNMixed");
 317     histoName = "fHistPN";
 318     if(gCentralityEstimator) histoName += gCentralityEstimator;
 319     hPNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
 320     hPNMixed->SetName("gHistPNMixed");
 321     histoName = "fHistNP";
 322     if(gCentralityEstimator) histoName += gCentralityEstimator;
 323     hNPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
 324     histoName = "fHistNP";
 325     if(gCentralityEstimator) histoName += gCentralityEstimator;
 326     hNPMixed->SetName("gHistNPMixed");
 327     histoName = "fHistPP";
 328     if(gCentralityEstimator) histoName += gCentralityEstimator;
 329     hPPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
 330     hPPMixed->SetName("gHistPPMixed");
 331     histoName = "fHistNN";
 332     if(gCentralityEstimator) histoName += gCentralityEstimator;
 333     hNNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
 334     hNNMixed->SetName("gHistNNMixed");
 335
 336     AliBalancePsi *bMixed = new AliBalancePsi();
 337     bMixed->SetEventClass(eventClass);
 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 *gHistBalanceFunction;
 347   TH2D *gHistBalanceFunctionSubtracted;
 348   TH2D *gHistBalanceFunctionShuffled;
 349   TH2D *gHistBalanceFunctionMixed;
 350   TString histoTitle, pngName;
 351
 352   if(eventClass == "Centrality"){
 353     histoTitle = "Centrality: ";
 354     histoTitle += psiMin;
 355     histoTitle += " - ";
 356     histoTitle += psiMax;
 357     histoTitle += " % ";
 358     histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 359   }
 360   else if(eventClass == "Multiplicity"){
 361     histoTitle = "Multiplicity: ";
 362     histoTitle += psiMin;
 363     histoTitle += " - ";
 364     histoTitle += psiMax;
 365     histoTitle += " tracks";
 366     histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 367   }
 368   else{ // "EventPlane" (default)
 369     histoTitle = "Centrality: ";
 370     histoTitle += centralityArray[gCentrality-1];
 371     histoTitle += "%";
 372     if((psiMin == -0.5)&&(psiMax == 0.5))
 373       histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 374     else if((psiMin == 0.5)&&(psiMax == 1.5))
 375       histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 376     else if((psiMin == 1.5)&&(psiMax == 2.5))
 377       histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 378     else
 379       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 380   }
 381
 382   if(k2pMethod)
 383     if(bMixed)
 384       gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 385     else{
 386       cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
 387       return;
 388     }
 389   else
 390     gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 391   gHistBalanceFunction->SetTitle(histoTitle.Data());
 392   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
 393   gHistBalanceFunction->SetName("gHistBalanceFunction");
 394
 395   if(listBFShuffled) {
 396
 397     if(k2pMethod)
 398       if(bMixed)
 399         gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 400       else{
 401         cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
 402         return;
 403       }
 404     else
 405       gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 406     gHistBalanceFunctionShuffled->SetTitle(histoTitle.Data());
 407     gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
 408     gHistBalanceFunctionShuffled->SetName("gHistBalanceFunctionShuffled");
 409   }
 410
 411   if(listBFMixed) {
 412     if(k2pMethod)
 413       if(bMixed)
 414         gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 415       else{
 416         cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
 417         return;
 418       }
 419     else
 420       gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 421     gHistBalanceFunctionMixed->SetTitle(histoTitle.Data());
 422     gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
 423     gHistBalanceFunctionMixed->SetName("gHistBalanceFunctionMixed");
 424
 425     gHistBalanceFunctionSubtracted = dynamic_cast<TH2D *>(gHistBalanceFunction->Clone());
 426     gHistBalanceFunctionSubtracted->Add(gHistBalanceFunctionMixed,-1);
 427     gHistBalanceFunctionSubtracted->SetTitle(histoTitle.Data());
 428     gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
 429     gHistBalanceFunctionSubtracted->SetName("gHistBalanceFunctionSubtracted");
 430   }
 431
 432   //Draw the results
 433   TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
 434   c1->SetFillColor(10);
 435   c1->SetHighLightColor(10);
 436   c1->SetLeftMargin(0.15);
 437   gHistBalanceFunction->DrawCopy("lego2");
 438   gPad->SetTheta(30); // default is 30
 439   gPad->SetPhi(-60); // default is 30
 440   gPad->Update();
 441   TCanvas *c1a = new TCanvas("c1a","",600,0,600,500);
 442   c1a->SetFillColor(10);
 443   c1a->SetHighLightColor(10);
 444   c1a->SetLeftMargin(0.15);
 445   gHistBalanceFunction->DrawCopy("colz");
 446
 447   if(listBFShuffled) {
 448     TCanvas *c2 = new TCanvas("c2","",100,100,600,500);
 449     c2->SetFillColor(10);
 450     c2->SetHighLightColor(10);
 451     c2->SetLeftMargin(0.15);
 452     gHistBalanceFunctionShuffled->DrawCopy("lego2");
 453     gPad->SetTheta(30); // default is 30
 454     gPad->SetPhi(-60); // default is 30
 455     gPad->Update();
 456     TCanvas *c2a = new TCanvas("c2a","",700,100,600,500);
 457     c2a->SetFillColor(10);
 458     c2a->SetHighLightColor(10);
 459     c2a->SetLeftMargin(0.15);
 460     gHistBalanceFunctionShuffled->DrawCopy("colz");
 461   }
 462
 463   if(listBFMixed) {
 464     TCanvas *c3 = new TCanvas("c3","",200,200,600,500);
 465     c3->SetFillColor(10);
 466     c3->SetHighLightColor(10);
 467     c3->SetLeftMargin(0.15);
 468     gHistBalanceFunctionMixed->DrawCopy("lego2");
 469     gPad->SetTheta(30); // default is 30
 470     gPad->SetPhi(-60); // default is 30
 471     gPad->Update();
 472     TCanvas *c3a = new TCanvas("c3a","",800,200,600,500);
 473     c3a->SetFillColor(10);
 474     c3a->SetHighLightColor(10);
 475     c3a->SetLeftMargin(0.15);
 476     gHistBalanceFunctionMixed->DrawCopy("colz");
 477
 478     TCanvas *c4 = new TCanvas("c4","",300,300,600,500);
 479     c4->SetFillColor(10);
 480     c4->SetHighLightColor(10);
 481     c4->SetLeftMargin(0.15);
 482     gHistBalanceFunctionSubtracted->DrawCopy("lego2");
 483     gPad->SetTheta(30); // default is 30
 484     gPad->SetPhi(-60); // default is 30
 485     gPad->Update();
 486     TCanvas *c4a = new TCanvas("c4a","",900,300,600,500);
 487     c4a->SetFillColor(10);
 488     c4a->SetHighLightColor(10);
 489     c4a->SetLeftMargin(0.15);
 490     gHistBalanceFunctionSubtracted->DrawCopy("colz");
 491
 492     fitbalanceFunction(gCentrality, psiMin , psiMax,
 493                        ptTriggerMin, ptTriggerMax,
 494                        ptAssociatedMin, ptAssociatedMax,
 495                        gHistBalanceFunctionSubtracted,k2pMethod, eventClass);
 496   }
 497
 498   TString newFileName = "balanceFunction2D.";
 499   if(eventClass == "Centrality"){
 500     newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
 501     newFileName += ".PsiAll.PttFrom";
 502   }
 503   else if(eventClass == "Multiplicity"){
 504     newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
 505     newFileName += ".PsiAll.PttFrom";
 506   }
 507   else{ // "EventPlane" (default)
 508     newFileName += "Centrality";
 509     newFileName += gCentrality; newFileName += ".Psi";
 510     if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
 511     else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
 512     else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
 513     else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
 514     else newFileName += "All.PttFrom";
 515   }
 516   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
 517   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
 518   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
 519   newFileName += Form("%.1f",ptAssociatedMax);
 520   if(k2pMethod) newFileName += "_2pMethod";
 521   newFileName += ".root";
 522
 523   TFile *fOutput = new TFile(newFileName.Data(),"recreate");
 524   fOutput->cd();
 525   /*hP->Write(); hN->Write();
 526   hPN->Write(); hNP->Write();
 527   hPP->Write(); hNN->Write();
 528   hPShuffled->Write(); hNShuffled->Write();
 529   hPNShuffled->Write(); hNPShuffled->Write();
 530   hPPShuffled->Write(); hNNShuffled->Write();
 531   hPMixed->Write(); hNMixed->Write();
 532   hPNMixed->Write(); hNPMixed->Write();
 533   hPPMixed->Write(); hNNMixed->Write();*/
 534   gHistBalanceFunction->Write();
 535   if(listBFShuffled) gHistBalanceFunctionShuffled->Write();
 536   if(listBFMixed) {
 537     gHistBalanceFunctionMixed->Write();
 538     gHistBalanceFunctionSubtracted->Write();
 539   }
 540   fOutput->Close();
 541 }
 542
 543 //____________________________________________________________//
 544 void fitbalanceFunction(Int_t gCentrality = 1,
 545                         Double_t psiMin = -0.5, Double_t psiMax = 3.5,
 546                         Double_t ptTriggerMin = -1.,
 547                         Double_t ptTriggerMax = -1.,
 548                         Double_t ptAssociatedMin = -1.,
 549                         Double_t ptAssociatedMax = -1.,
 550                         TH2D *gHist,
 551                         Bool_t k2pMethod = kFALSE,
 552                         TString eventClass="EventPlane") {
 553
 554   cout<<"FITTING FUNCTION"<<endl;
 555
 556   //near side peak: [1]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[2]),2)+0.5*TMath::Power((y/[3]),2)),[4]))
 557   //away side ridge: [5]*TMath::Exp(-TMath::Power((0.5*TMath::Power(((y-TMath::Pi())/[6]),2)),[7]))
 558   //longitudinal ridge: [8]*TMath::Exp(-TMath::Power((0.5*TMath::Power((x/[9]),2)),[10]))
 559   //wing structures: [11]*TMath::Power(x,2)
 560   //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))
 561   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.);
 562   gFitFunction->SetName("gFitFunction");
 563
 564
 565   //Normalization
 566   gFitFunction->SetParName(0,"N1");
 567   //near side peak
 568   gFitFunction->SetParName(1,"N_{near side}");gFitFunction->FixParameter(1,0.0);
 569   gFitFunction->SetParName(2,"Sigma_{near side}(delta eta)"); gFitFunction->FixParameter(2,0.0);
 570   gFitFunction->SetParName(3,"Sigma_{near side}(delta phi)"); gFitFunction->FixParameter(3,0.0);
 571   gFitFunction->SetParName(4,"Exponent_{near side}"); gFitFunction->FixParameter(4,1.0);
 572   //away side ridge
 573   gFitFunction->SetParName(5,"N_{away side}"); gFitFunction->FixParameter(5,0.0);
 574   gFitFunction->SetParName(6,"Sigma_{away side}(delta phi)"); gFitFunction->FixParameter(6,0.0);
 575   gFitFunction->SetParName(7,"Exponent_{away side}"); gFitFunction->FixParameter(7,1.0);
 576   //longitudinal ridge
 577   gFitFunction->SetParName(8,"N_{long. ridge}"); gFitFunction->SetParameter(8,0.05);//
 578   gFitFunction->FixParameter(8,0.0);
 579   gFitFunction->SetParName(9,"Sigma_{long. ridge}(delta eta)"); gFitFunction->FixParameter(9,0.0);
 580   gFitFunction->SetParName(10,"Exponent_{long. ridge}"); gFitFunction->FixParameter(10,1.0);
 581   //wing structures
 582   gFitFunction->SetParName(11,"N_{wing}"); gFitFunction->FixParameter(11,0.0);
 583   //flow contribution
 584   gFitFunction->SetParName(12,"N_{flow}"); gFitFunction->FixParameter(12,0.0);
 585   gFitFunction->SetParName(13,"V1"); gFitFunction->FixParameter(13,0.0);
 586   gFitFunction->SetParName(14,"V2"); gFitFunction->FixParameter(14,0.0);
 587   gFitFunction->SetParName(15,"V3"); gFitFunction->FixParameter(15,0.0);
 588   gFitFunction->SetParName(16,"V4"); gFitFunction->FixParameter(16,0.0);
 589   gFitFunction->SetParName(17,"Offset"); gFitFunction->FixParameter(17,0.0);
 590
 591   // just release the near side peak
 592   gFitFunction->ReleaseParameter(0);gFitFunction->SetParameter(0,1.0);
 593   gFitFunction->ReleaseParameter(1);gFitFunction->SetParameter(1,0.3);gFitFunction->SetParLimits(1,0.0,100000);
 594   gFitFunction->ReleaseParameter(2);gFitFunction->SetParameter(2,0.3);gFitFunction->SetParLimits(2,0.05,0.7);
 595   gFitFunction->ReleaseParameter(3);gFitFunction->SetParameter(3,0.3);gFitFunction->SetParLimits(3,0.05,1.7);
 596
 597   gHist->Fit("gFitFunction","nm");
 598
 599
 600   //Cloning the histogram
 601   TString histoName = gHist->GetName(); histoName += "Fit";
 602   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());
 603   TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
 604   gHistResidual->SetName("gHistResidual");
 605   gHistResidual->Sumw2();
 606
 607   for(Int_t iBinDeltaEta = 1; iBinDeltaEta <= gHist->GetNbinsX(); iBinDeltaEta++) {
 608     for(Int_t iBinDeltaPhi = 1; iBinDeltaPhi <= gHist->GetNbinsY(); iBinDeltaPhi++) {
 609       gHistFit->SetBinContent(iBinDeltaEta,iBinDeltaPhi,gFitFunction->Eval(gHist->GetXaxis()->GetBinCenter(iBinDeltaEta),gHist->GetYaxis()->GetBinCenter(iBinDeltaPhi)));
 610     }
 611   }
 612   gHistResidual->Add(gHistFit,-1);
 613
 614   //Write to output file
 615   TString newFileName = "balanceFunctionFit2D.";
 616   if(eventClass == "Centrality"){
 617     newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
 618     newFileName += ".PsiAll.PttFrom";
 619   }
 620   else if(eventClass == "Multiplicity"){
 621     newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
 622     newFileName += ".PsiAll.PttFrom";
 623   }
 624   else{ // "EventPlane" (default)
 625     newFileName += "Centrality";
 626     newFileName += gCentrality; newFileName += ".Psi";
 627     if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
 628     else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
 629     else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
 630     else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
 631     else newFileName += "All.PttFrom";
 632   }
 633   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
 634   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
 635   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
 636   newFileName += Form("%.1f",ptAssociatedMax);
 637   if(k2pMethod) newFileName += "_2pMethod";
 638   newFileName += ".root";
 639   TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
 640   gHist->Write();
 641   gHistFit->Write();
 642   gHistResidual->Write();
 643   gFitFunction->Write();
 644   newFile->Close();
 645 }
 646
 647 //____________________________________________________________//
 648 void drawBFPsi2D(const char* lhcPeriod = "LHC11h",
 649                  const char* gCentralityEstimator = "V0M",
 650                  Int_t gBit = 128,
 651                  const char* gEventPlaneEstimator = "VZERO",
 652                  Int_t gCentrality = 1,
 653                  Bool_t kShowShuffled = kFALSE,
 654                  Bool_t kShowMixed = kFALSE,
 655                  Double_t psiMin = -0.5, Double_t psiMax = 0.5,
 656                  Double_t ptTriggerMin = -1.,
 657                  Double_t ptTriggerMax = -1.,
 658                  Double_t ptAssociatedMin = -1.,
 659                  Double_t ptAssociatedMax = -1.,
 660                  Bool_t k2pMethod = kTRUE) {
 661   //Macro that draws the BF distributions for each centrality bin
 662   //for reaction plane dependent analysis
 663   //Author: Panos.Christakoglou@nikhef.nl
 664   TGaxis::SetMaxDigits(3);
 665
 666   //Get the input file
 667   TString filename = lhcPeriod;
 668   filename += "/Centrality"; filename += gCentralityEstimator;
 669   filename += "_Bit"; filename += gBit;
 670   filename += "_"; filename += gEventPlaneEstimator;
 671   filename +="/PttFrom";
 672   filename += Form("%.1f",ptTriggerMin); filename += "To";
 673   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
 674   filename += Form("%.1f",ptAssociatedMin); filename += "To";
 675   filename += Form("%.1f",ptAssociatedMax);
 676   filename += "/balanceFunction2D.Centrality";
 677   filename += gCentrality; filename += ".Psi";
 678   if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
 679   else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
 680   else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
 681   else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
 682   else filename += "All.PttFrom";
 683   filename += Form("%.1f",ptTriggerMin); filename += "To";
 684   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
 685   filename += Form("%.1f",ptAssociatedMin); filename += "To";
 686   filename += Form("%.1f",ptAssociatedMax);
 687   if(k2pMethod) filename += "_2pMethod";
 688   filename += ".root";
 689
 690   //Open the file
 691   TFile *f = TFile::Open(filename.Data());
 692   if((!f)||(!f->IsOpen())) {
 693     Printf("The file %s is not found. Aborting...",filename);
 694     return listBF;
 695   }
 696   //f->ls();
 697
 698   //Raw balance function
 699   TH1D *gHistBalanceFunction = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunction"));
 700   gHistBalanceFunction->SetStats(kFALSE);
 701   gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
 702   gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
 703   gHistBalanceFunction->GetZaxis()->SetNdivisions(10);
 704   gHistBalanceFunction->GetXaxis()->SetTitleOffset(1.3);
 705   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
 706   gHistBalanceFunction->GetZaxis()->SetTitleOffset(1.3);
 707   gHistBalanceFunction->GetXaxis()->SetTitle("#Delta #eta");
 708   gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 709   gHistBalanceFunction->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
 710
 711   //Shuffled balance function
 712   if(kShowShuffled) {
 713     TH1D *gHistBalanceFunctionShuffled = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionShuffled"));
 714     gHistBalanceFunctionShuffled->SetStats(kFALSE);
 715     gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
 716     gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
 717     gHistBalanceFunctionShuffled->GetZaxis()->SetNdivisions(10);
 718     gHistBalanceFunctionShuffled->GetXaxis()->SetTitleOffset(1.3);
 719     gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
 720     gHistBalanceFunctionShuffled->GetZaxis()->SetTitleOffset(1.3);
 721     gHistBalanceFunctionShuffled->GetXaxis()->SetTitle("#Delta #eta");
 722     gHistBalanceFunctionShuffled->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 723     gHistBalanceFunctionShuffled->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
 724   }
 725
 726   //Mixed balance function
 727   if(kShowMixed) {
 728     TH1D *gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionMixed"));
 729     gHistBalanceFunctionMixed->SetStats(kFALSE);
 730     gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
 731     gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
 732     gHistBalanceFunctionMixed->GetZaxis()->SetNdivisions(10);
 733     gHistBalanceFunctionMixed->GetXaxis()->SetTitleOffset(1.3);
 734     gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
 735     gHistBalanceFunctionMixed->GetZaxis()->SetTitleOffset(1.3);
 736     gHistBalanceFunctionMixed->GetXaxis()->SetTitle("#Delta #eta");
 737     gHistBalanceFunctionMixed->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 738     gHistBalanceFunctionMixed->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
 739   }
 740
 741   //Subtracted balance function
 742   if(kShowMixed) {
 743     TH1D *gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionSubtracted"));
 744     gHistBalanceFunctionSubtracted->SetStats(kFALSE);
 745     gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
 746     gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
 747     gHistBalanceFunctionSubtracted->GetZaxis()->SetNdivisions(10);
 748     gHistBalanceFunctionSubtracted->GetXaxis()->SetTitleOffset(1.3);
 749     gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
 750     gHistBalanceFunctionSubtracted->GetZaxis()->SetTitleOffset(1.3);
 751     gHistBalanceFunctionSubtracted->GetXaxis()->SetTitle("#Delta #eta");
 752     gHistBalanceFunctionSubtracted->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 753     gHistBalanceFunctionSubtracted->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
 754   }
 755
 756   TString pngName;
 757
 758   TString centralityLatex = "Centrality: ";
 759   centralityLatex += centralityArray[gCentrality-1];
 760   centralityLatex += "%";
 761
 762   TString psiLatex;
 763   if((psiMin == -0.5)&&(psiMax == 0.5))
 764     psiLatex = " -7.5^{o} < #varphi - #Psi_{2} < 7.5^{o}";
 765   else if((psiMin == 0.5)&&(psiMax == 1.5))
 766     psiLatex = " 37.5^{o} < #varphi - #Psi_{2} < 52.5^{o}";
 767   else if((psiMin == 1.5)&&(psiMax == 2.5))
 768     psiLatex = " 82.5^{o} < #varphi - #Psi_{2} < 97.5^{o}";
 769   else
 770     psiLatex = " 0^{o} < #varphi - #Psi_{2} < 180^{o}";
 771
 772   TString pttLatex = Form("%.1f",ptTriggerMin);
 773   pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
 774   pttLatex += " GeV/c";
 775
 776   TString ptaLatex = Form("%.1f",ptAssociatedMin);
 777   ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
 778   ptaLatex += " GeV/c";
 779
 780   TLatex *latexInfo1 = new TLatex();
 781   latexInfo1->SetNDC();
 782   latexInfo1->SetTextSize(0.045);
 783   latexInfo1->SetTextColor(1);
 784
 785   //Draw the results
 786   TCanvas *c1 = new TCanvas("c1","Raw balance function 2D",0,0,600,500);
 787   c1->SetFillColor(10); c1->SetHighLightColor(10);
 788   c1->SetLeftMargin(0.17); c1->SetTopMargin(0.05);
 789   gHistBalanceFunction->SetTitle("");
 790   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.4);
 791   gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
 792   gHistBalanceFunction->GetXaxis()->SetRangeUser(-1.4,1.4);
 793   gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
 794   gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 795   gHistBalanceFunction->DrawCopy("lego2");
 796   gPad->SetTheta(30); // default is 30
 797   gPad->SetPhi(-60); // default is 30
 798   gPad->Update();
 799
 800   latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
 801   latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
 802   latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
 803   latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
 804
 805   TString pngName = "BalanceFunction2D.";
 806   pngName += "Centrality";
 807   pngName += gCentrality;
 808   if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
 809   else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
 810   else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
 811   else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.PttFrom";
 812   else pngName += "All.PttFrom";
 813   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
 814   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
 815   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
 816   pngName += Form("%.1f",ptAssociatedMax);
 817   if(k2pMethod) pngName += "_2pMethod";
 818   pngName += ".png";
 819
 820   c1->SaveAs(pngName.Data());
 821
 822   if(kShowShuffled) {
 823     TCanvas *c2 = new TCanvas("c2","Shuffled balance function 2D",100,100,600,500);
 824     c2->SetFillColor(10); c2->SetHighLightColor(10);
 825     c2->SetLeftMargin(0.17); c2->SetTopMargin(0.05);
 826     gHistBalanceFunctionShuffled->SetTitle("Shuffled events");
 827     gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.4);
 828     gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
 829     gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
 830     gHistBalanceFunctionShuffled->DrawCopy("lego2");
 831     gPad->SetTheta(30); // default is 30
 832     gPad->SetPhi(-60); // default is 30
 833     gPad->Update();
 834
 835     latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
 836     latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
 837     latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
 838     latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
 839   }
 840
 841   if(kShowMixed) {
 842     TCanvas *c3 = new TCanvas("c3","Mixed balance function 2D",200,200,600,500);
 843     c3->SetFillColor(10); c3->SetHighLightColor(10);
 844     c3->SetLeftMargin(0.17); c3->SetTopMargin(0.05);
 845     gHistBalanceFunctionMixed->SetTitle("Mixed events");
 846     gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.4);
 847     gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
 848     gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
 849     gHistBalanceFunctionMixed->DrawCopy("lego2");
 850     gPad->SetTheta(30); // default is 30
 851     gPad->SetPhi(-60); // default is 30
 852     gPad->Update();
 853
 854     latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
 855     latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
 856     latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
 857     latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
 858   }
 859
 860   if(kShowMixed) {
 861     TCanvas *c4 = new TCanvas("c4","Subtracted balance function 2D",300,300,600,500);
 862     c4->SetFillColor(10); c4->SetHighLightColor(10);
 863     c4->SetLeftMargin(0.17); c4->SetTopMargin(0.05);
 864     gHistBalanceFunctionSubtracted->SetTitle("Subtracted balance function");
 865     gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.4);
 866     gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
 867     gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
 868     gHistBalanceFunctionSubtracted->DrawCopy("lego2");
 869     gPad->SetTheta(30); // default is 30
 870     gPad->SetPhi(-60); // default is 30
 871     gPad->Update();
 872
 873     latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
 874     latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
 875     latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
 876     latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
 877   }
 878 }