PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C

   1 const Int_t numberOfCentralityBins = 12;
   2 TString centralityArray[numberOfCentralityBins] = {"0-100","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 vertexZMin = -10.,
  14                               Double_t vertexZMax = 10.,
  15                               Double_t ptTriggerMin = -1.,
  16                               Double_t ptTriggerMax = -1.,
  17                               Double_t ptAssociatedMin = -1.,
  18                               Double_t ptAssociatedMax = -1.,
  19                               Bool_t kUseVzBinning = kTRUE,
  20                               Bool_t k2pMethod = kTRUE,
  21                               TString eventClass = "EventPlane") //Can be "EventPlane", "Centrality", "Multiplicity"
  22 {
  23   //Macro that draws the BF distributions for each centrality bin
  24   //for reaction plane dependent analysis
  25   //Author: Panos.Christakoglou@nikhef.nl
  26   //Load the PWG2ebye library
  27   gSystem->Load("libANALYSIS.so");
  28   gSystem->Load("libANALYSISalice.so");
  29   gSystem->Load("libEventMixing.so");
  30   gSystem->Load("libCORRFW.so");
  31   gSystem->Load("libPWGTools.so");
  32   gSystem->Load("libPWGCFebye.so");
  33
  34   //gROOT->LoadMacro("~/SetPlotStyle.C");
  35   //SetPlotStyle();
  36   gStyle->SetPalette(1,0);
  37
  38   //Prepare the objects and return them
  39   TList *listBF = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,0);
  40   TList *listBFShuffled = NULL;
  41   if(kShowShuffled) listBFShuffled = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,1);
  42   TList *listBFMixed = NULL;
  43   if(kShowMixed) listBFMixed = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,2);
  44   if(!listBF) {
  45     Printf("The TList object was not created");
  46     return;
  47   }
  48   else
  49     draw(listBF,listBFShuffled,listBFMixed,gCentrality,gCentralityEstimator,
  50          psiMin,psiMax,vertexZMin,vertexZMax,
  51          ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,
  52          kUseVzBinning,k2pMethod,eventClass);
  53 }
  54
  55 //______________________________________________________//
  56 TList *GetListOfObjects(const char* filename,
  57                         Int_t gCentrality,
  58                         Int_t gBit,
  59                         const char *gCentralityEstimator,
  60                         Int_t kData = 1) {
  61   //Get the TList objects (QA, bf, bf shuffled)
  62   TList *listBF = 0x0;
  63
  64   //Open the file
  65   TFile *f = TFile::Open(filename,"UPDATE");
  66   if((!f)||(!f->IsOpen())) {
  67     Printf("The file %s is not found. Aborting...",filename);
  68     return listBF;
  69   }
  70   //f->ls();
  71
  72   TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
  73   if(!dir) {
  74     Printf("The TDirectoryFile is not found. Aborting...",filename);
  75     return listBF;
  76   }
  77   //dir->ls();
  78
  79   TString listBFName;
  80   if(kData == 0) {
  81     //cout<<"no shuffling - no mixing"<<endl;
  82     listBFName = "listBFPsi_";
  83   }
  84   else if(kData == 1) {
  85     //cout<<"shuffling - no mixing"<<endl;
  86     listBFName = "listBFPsiShuffled_";
  87   }
  88   else if(kData == 2) {
  89     //cout<<"no shuffling - mixing"<<endl;
  90     listBFName = "listBFPsiMixed_";
  91   }
  92   listBFName += centralityArray[gCentrality-1];
  93   if(gBit > -1) {
  94     listBFName += "_Bit"; listBFName += gBit; }
  95   if(gCentralityEstimator) {
  96     listBFName += "_"; listBFName += gCentralityEstimator;}
  97
  98   // histograms were already retrieved (in first iteration)
  99   if(dir->Get(Form("%s_histograms",listBFName.Data()))){
 100     //cout<<"second iteration"<<endl;
 101     listBF = dynamic_cast<TList *>(dir->Get(Form("%s_histograms",listBFName.Data())));
 102   }
 103
 104   // histograms were not yet retrieved (this is the first iteration)
 105   else{
 106
 107     listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
 108     cout<<"======================================================="<<endl;
 109     cout<<"List name: "<<listBF->GetName()<<endl;
 110     //listBF->ls();
 111
 112     //Get the histograms
 113     TString histoName;
 114     if(kData == 0)
 115       histoName = "fHistP";
 116     else if(kData == 1)
 117       histoName = "fHistP_shuffle";
 118     else if(kData == 2)
 119       histoName = "fHistP";
 120     if(gCentralityEstimator) histoName += gCentralityEstimator;
 121     AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 122     if(!fHistP) {
 123       Printf("fHistP %s not found!!!",histoName.Data());
 124       break;
 125     }
 126     fHistP->FillParent(); fHistP->DeleteContainers();
 127
 128     if(kData == 0)
 129       histoName = "fHistN";
 130     if(kData == 1)
 131       histoName = "fHistN_shuffle";
 132     if(kData == 2)
 133       histoName = "fHistN";
 134     if(gCentralityEstimator) histoName += gCentralityEstimator;
 135     AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 136     if(!fHistN) {
 137       Printf("fHistN %s not found!!!",histoName.Data());
 138       break;
 139     }
 140     fHistN->FillParent(); fHistN->DeleteContainers();
 141
 142     if(kData == 0)
 143       histoName = "fHistPN";
 144     if(kData == 1)
 145       histoName = "fHistPN_shuffle";
 146     if(kData == 2)
 147       histoName = "fHistPN";
 148     if(gCentralityEstimator) histoName += gCentralityEstimator;
 149     AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 150     if(!fHistPN) {
 151       Printf("fHistPN %s not found!!!",histoName.Data());
 152       break;
 153     }
 154     fHistPN->FillParent(); fHistPN->DeleteContainers();
 155
 156     if(kData == 0)
 157       histoName = "fHistNP";
 158     if(kData == 1)
 159       histoName = "fHistNP_shuffle";
 160     if(kData == 2)
 161       histoName = "fHistNP";
 162     if(gCentralityEstimator) histoName += gCentralityEstimator;
 163     AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 164     if(!fHistNP) {
 165       Printf("fHistNP %s not found!!!",histoName.Data());
 166       break;
 167     }
 168     fHistNP->FillParent(); fHistNP->DeleteContainers();
 169
 170     if(kData == 0)
 171       histoName = "fHistPP";
 172     if(kData == 1)
 173       histoName = "fHistPP_shuffle";
 174     if(kData == 2)
 175       histoName = "fHistPP";
 176     if(gCentralityEstimator) histoName += gCentralityEstimator;
 177     AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 178     if(!fHistPP) {
 179       Printf("fHistPP %s not found!!!",histoName.Data());
 180       break;
 181     }
 182     fHistPP->FillParent(); fHistPP->DeleteContainers();
 183
 184     if(kData == 0)
 185       histoName = "fHistNN";
 186     if(kData == 1)
 187       histoName = "fHistNN_shuffle";
 188     if(kData == 2)
 189       histoName = "fHistNN";
 190     if(gCentralityEstimator) histoName += gCentralityEstimator;
 191     AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 192     if(!fHistNN) {
 193       Printf("fHistNN %s not found!!!",histoName.Data());
 194       break;
 195     }
 196     fHistNN->FillParent(); fHistNN->DeleteContainers();
 197
 198     dir->cd();
 199     listBF->Write(Form("%s_histograms",listBFName.Data()), TObject::kSingleKey);
 200
 201   }// first iteration
 202
 203   f->Close();
 204
 205   return listBF;
 206 }
 207
 208 //______________________________________________________//
 209 void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
 210           Int_t gCentrality, const char* gCentralityEstimator,
 211           Double_t psiMin, Double_t psiMax,
 212           Double_t vertexZMin,
 213           Double_t vertexZMax,
 214           Double_t ptTriggerMin, Double_t ptTriggerMax,
 215           Double_t ptAssociatedMin, Double_t ptAssociatedMax,
 216           Bool_t kUseVzBinning=kFALSE,
 217           Bool_t k2pMethod = kFALSE, TString eventClass) {
 218   //balance function
 219   AliTHn *hP = NULL;
 220   AliTHn *hN = NULL;
 221   AliTHn *hPN = NULL;
 222   AliTHn *hNP = NULL;
 223   AliTHn *hPP = NULL;
 224   AliTHn *hNN = NULL;
 225   //listBF->ls();
 226   //Printf("=================");
 227   TString histoName = "fHistP";
 228   if(gCentralityEstimator) histoName += gCentralityEstimator;
 229   hP = (AliTHn*) listBF->FindObject(histoName.Data());
 230   hP->SetName("gHistP");
 231   histoName = "fHistN";
 232   if(gCentralityEstimator) histoName += gCentralityEstimator;
 233   hN = (AliTHn*) listBF->FindObject(histoName.Data());
 234   hN->SetName("gHistN");
 235   histoName = "fHistPN";
 236   if(gCentralityEstimator) histoName += gCentralityEstimator;
 237   hPN = (AliTHn*) listBF->FindObject(histoName.Data());
 238   hPN->SetName("gHistPN");
 239   histoName = "fHistNP";
 240   if(gCentralityEstimator) histoName += gCentralityEstimator;
 241   hNP = (AliTHn*) listBF->FindObject(histoName.Data());
 242   hNP->SetName("gHistNP");
 243   histoName = "fHistPP";
 244   if(gCentralityEstimator) histoName += gCentralityEstimator;
 245   hPP = (AliTHn*) listBF->FindObject(histoName.Data());
 246   hPP->SetName("gHistPP");
 247   histoName = "fHistNN";
 248   if(gCentralityEstimator) histoName += gCentralityEstimator;
 249   hNN = (AliTHn*) listBF->FindObject(histoName.Data());
 250   hNN->SetName("gHistNN");
 251
 252   AliBalancePsi *b = new AliBalancePsi();
 253   b->SetEventClass(eventClass);
 254   b->SetHistNp(hP);
 255   b->SetHistNn(hN);
 256   b->SetHistNpn(hPN);
 257   b->SetHistNnp(hNP);
 258   b->SetHistNpp(hPP);
 259   b->SetHistNnn(hNN);
 260   if(kUseVzBinning) b->SetVertexZBinning(kTRUE);
 261
 262
 263   //balance function shuffling
 264   AliTHn *hPShuffled = NULL;
 265   AliTHn *hNShuffled = NULL;
 266   AliTHn *hPNShuffled = NULL;
 267   AliTHn *hNPShuffled = NULL;
 268   AliTHn *hPPShuffled = NULL;
 269   AliTHn *hNNShuffled = NULL;
 270   if(listBFShuffled) {
 271     //listBFShuffled->ls();
 272     histoName = "fHistP_shuffle";
 273     if(gCentralityEstimator) histoName += gCentralityEstimator;
 274     hPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
 275     hPShuffled->SetName("gHistPShuffled");
 276     histoName = "fHistN_shuffle";
 277     if(gCentralityEstimator) histoName += gCentralityEstimator;
 278     hNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
 279     hNShuffled->SetName("gHistNShuffled");
 280     histoName = "fHistPN_shuffle";
 281     if(gCentralityEstimator) histoName += gCentralityEstimator;
 282     hPNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
 283     hPNShuffled->SetName("gHistPNShuffled");
 284     histoName = "fHistNP_shuffle";
 285     if(gCentralityEstimator) histoName += gCentralityEstimator;
 286     hNPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
 287     hNPShuffled->SetName("gHistNPShuffled");
 288     histoName = "fHistPP_shuffle";
 289     if(gCentralityEstimator) histoName += gCentralityEstimator;
 290     hPPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
 291     hPPShuffled->SetName("gHistPPShuffled");
 292     histoName = "fHistNN_shuffle";
 293     if(gCentralityEstimator) histoName += gCentralityEstimator;
 294     hNNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
 295     hNNShuffled->SetName("gHistNNShuffled");
 296
 297     AliBalancePsi *bShuffled = new AliBalancePsi();
 298     bShuffled->SetEventClass(eventClass);
 299     bShuffled->SetHistNp(hPShuffled);
 300     bShuffled->SetHistNn(hNShuffled);
 301     bShuffled->SetHistNpn(hPNShuffled);
 302     bShuffled->SetHistNnp(hNPShuffled);
 303     bShuffled->SetHistNpp(hPPShuffled);
 304     bShuffled->SetHistNnn(hNNShuffled);
 305   if(kUseVzBinning) bShuffled->SetVertexZBinning(kTRUE);
 306
 307   }
 308
 309   //balance function mixing
 310   AliTHn *hPMixed = NULL;
 311   AliTHn *hNMixed = NULL;
 312   AliTHn *hPNMixed = NULL;
 313   AliTHn *hNPMixed = NULL;
 314   AliTHn *hPPMixed = NULL;
 315   AliTHn *hNNMixed = NULL;
 316
 317   if(listBFMixed) {
 318     //listBFMixed->ls();
 319     histoName = "fHistP";
 320     if(gCentralityEstimator) histoName += gCentralityEstimator;
 321     hPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
 322     hPMixed->SetName("gHistPMixed");
 323     histoName = "fHistN";
 324     if(gCentralityEstimator) histoName += gCentralityEstimator;
 325     hNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
 326     hNMixed->SetName("gHistNMixed");
 327     histoName = "fHistPN";
 328     if(gCentralityEstimator) histoName += gCentralityEstimator;
 329     hPNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
 330     hPNMixed->SetName("gHistPNMixed");
 331     histoName = "fHistNP";
 332     if(gCentralityEstimator) histoName += gCentralityEstimator;
 333     hNPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
 334     histoName = "fHistNP";
 335     if(gCentralityEstimator) histoName += gCentralityEstimator;
 336     hNPMixed->SetName("gHistNPMixed");
 337     histoName = "fHistPP";
 338     if(gCentralityEstimator) histoName += gCentralityEstimator;
 339     hPPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
 340     hPPMixed->SetName("gHistPPMixed");
 341     histoName = "fHistNN";
 342     if(gCentralityEstimator) histoName += gCentralityEstimator;
 343     hNNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
 344     hNNMixed->SetName("gHistNNMixed");
 345
 346     AliBalancePsi *bMixed = new AliBalancePsi();
 347     bMixed->SetEventClass(eventClass);
 348     bMixed->SetHistNp(hPMixed);
 349     bMixed->SetHistNn(hNMixed);
 350     bMixed->SetHistNpn(hPNMixed);
 351     bMixed->SetHistNnp(hNPMixed);
 352     bMixed->SetHistNpp(hPPMixed);
 353     bMixed->SetHistNnn(hNNMixed);
 354     if(kUseVzBinning) bMixed->SetVertexZBinning(kTRUE);
 355
 356   }
 357
 358   TH2D *gHistBalanceFunction;
 359   TH2D *gHistBalanceFunctionSubtracted;
 360   TH2D *gHistBalanceFunctionShuffled;
 361   TH2D *gHistBalanceFunctionMixed;
 362   TString histoTitle, pngName;
 363
 364   if(eventClass == "Centrality"){
 365     histoTitle = "Centrality: ";
 366     histoTitle += psiMin;
 367     histoTitle += " - ";
 368     histoTitle += psiMax;
 369     histoTitle += " % ";
 370     histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 371   }
 372   else if(eventClass == "Multiplicity"){
 373     histoTitle = "Multiplicity: ";
 374     histoTitle += psiMin;
 375     histoTitle += " - ";
 376     histoTitle += psiMax;
 377     histoTitle += " tracks";
 378     histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 379   }
 380   else{ // "EventPlane" (default)
 381     histoTitle = "Centrality: ";
 382     histoTitle += centralityArray[gCentrality-1];
 383     histoTitle += "%";
 384     if((psiMin == -0.5)&&(psiMax == 0.5))
 385       histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 386     else if((psiMin == 0.5)&&(psiMax == 1.5))
 387       histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 388     else if((psiMin == 1.5)&&(psiMax == 2.5))
 389       histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 390     else
 391       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 392   }
 393
 394   if(k2pMethod)
 395     if(bMixed)
 396       gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 397     else{
 398       cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
 399       return;
 400     }
 401   else
 402     gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 403   gHistBalanceFunction->SetTitle(histoTitle.Data());
 404   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
 405   gHistBalanceFunction->SetName("gHistBalanceFunction");
 406
 407   if(listBFShuffled) {
 408
 409     if(k2pMethod)
 410       if(bMixed)
 411         gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 412       else{
 413         cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
 414         return;
 415       }
 416     else
 417       gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 418     gHistBalanceFunctionShuffled->SetTitle(histoTitle.Data());
 419     gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
 420     gHistBalanceFunctionShuffled->SetName("gHistBalanceFunctionShuffled");
 421   }
 422
 423   if(listBFMixed) {
 424     if(k2pMethod)
 425       if(bMixed)
 426         gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 427       else{
 428         cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
 429         return;
 430       }
 431     else
 432       gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 433     gHistBalanceFunctionMixed->SetTitle(histoTitle.Data());
 434     gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
 435     gHistBalanceFunctionMixed->SetName("gHistBalanceFunctionMixed");
 436
 437     gHistBalanceFunctionSubtracted = dynamic_cast<TH2D *>(gHistBalanceFunction->Clone());
 438     gHistBalanceFunctionSubtracted->Add(gHistBalanceFunctionMixed,-1);
 439     gHistBalanceFunctionSubtracted->SetTitle(histoTitle.Data());
 440     gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
 441     gHistBalanceFunctionSubtracted->SetName("gHistBalanceFunctionSubtracted");
 442   }
 443
 444   //Draw the results
 445   TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
 446   c1->SetFillColor(10);
 447   c1->SetHighLightColor(10);
 448   c1->SetLeftMargin(0.15);
 449   gHistBalanceFunction->DrawCopy("lego2");
 450   gPad->SetTheta(30); // default is 30
 451   gPad->SetPhi(-60); // default is 30
 452   gPad->Update();
 453   TCanvas *c1a = new TCanvas("c1a","",600,0,600,500);
 454   c1a->SetFillColor(10);
 455   c1a->SetHighLightColor(10);
 456   c1a->SetLeftMargin(0.15);
 457   gHistBalanceFunction->DrawCopy("colz");
 458
 459   if(listBFShuffled) {
 460     TCanvas *c2 = new TCanvas("c2","",100,100,600,500);
 461     c2->SetFillColor(10);
 462     c2->SetHighLightColor(10);
 463     c2->SetLeftMargin(0.15);
 464     gHistBalanceFunctionShuffled->DrawCopy("lego2");
 465     gPad->SetTheta(30); // default is 30
 466     gPad->SetPhi(-60); // default is 30
 467     gPad->Update();
 468     TCanvas *c2a = new TCanvas("c2a","",700,100,600,500);
 469     c2a->SetFillColor(10);
 470     c2a->SetHighLightColor(10);
 471     c2a->SetLeftMargin(0.15);
 472     gHistBalanceFunctionShuffled->DrawCopy("colz");
 473   }
 474
 475   if(listBFMixed) {
 476     TCanvas *c3 = new TCanvas("c3","",200,200,600,500);
 477     c3->SetFillColor(10);
 478     c3->SetHighLightColor(10);
 479     c3->SetLeftMargin(0.15);
 480     gHistBalanceFunctionMixed->DrawCopy("lego2");
 481     gPad->SetTheta(30); // default is 30
 482     gPad->SetPhi(-60); // default is 30
 483     gPad->Update();
 484     TCanvas *c3a = new TCanvas("c3a","",800,200,600,500);
 485     c3a->SetFillColor(10);
 486     c3a->SetHighLightColor(10);
 487     c3a->SetLeftMargin(0.15);
 488     gHistBalanceFunctionMixed->DrawCopy("colz");
 489
 490     TCanvas *c4 = new TCanvas("c4","",300,300,600,500);
 491     c4->SetFillColor(10);
 492     c4->SetHighLightColor(10);
 493     c4->SetLeftMargin(0.15);
 494     gHistBalanceFunctionSubtracted->DrawCopy("lego2");
 495     gPad->SetTheta(30); // default is 30
 496     gPad->SetPhi(-60); // default is 30
 497     gPad->Update();
 498     TCanvas *c4a = new TCanvas("c4a","",900,300,600,500);
 499     c4a->SetFillColor(10);
 500     c4a->SetHighLightColor(10);
 501     c4a->SetLeftMargin(0.15);
 502     gHistBalanceFunctionSubtracted->DrawCopy("colz");
 503
 504     fitbalanceFunction(gCentrality, psiMin , psiMax,
 505                        ptTriggerMin, ptTriggerMax,
 506                        ptAssociatedMin, ptAssociatedMax,
 507                        gHistBalanceFunctionSubtracted,k2pMethod, eventClass);
 508   }
 509
 510   TString newFileName = "balanceFunction2D.";
 511   if(eventClass == "Centrality"){
 512     newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
 513     newFileName += ".PsiAll.PttFrom";
 514   }
 515   else if(eventClass == "Multiplicity"){
 516     newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
 517     newFileName += ".PsiAll.PttFrom";
 518   }
 519   else{ // "EventPlane" (default)
 520     newFileName += "Centrality";
 521     newFileName += gCentrality; newFileName += ".Psi";
 522     if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
 523     else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
 524     else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
 525     else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
 526     else newFileName += "All.PttFrom";
 527   }
 528   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
 529   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
 530   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
 531   newFileName += Form("%.1f",ptAssociatedMax);
 532   if(k2pMethod) newFileName += "_2pMethod";
 533
 534   newFileName += "_";
 535   newFileName += Form("%.1f",psiMin);
 536   newFileName += "-";
 537   newFileName += Form("%.1f",psiMax);
 538   newFileName += ".root";
 539
 540   TFile *fOutput = new TFile(newFileName.Data(),"recreate");
 541   fOutput->cd();
 542   /*hP->Write(); hN->Write();
 543   hPN->Write(); hNP->Write();
 544   hPP->Write(); hNN->Write();
 545   hPShuffled->Write(); hNShuffled->Write();
 546   hPNShuffled->Write(); hNPShuffled->Write();
 547   hPPShuffled->Write(); hNNShuffled->Write();
 548   hPMixed->Write(); hNMixed->Write();
 549   hPNMixed->Write(); hNPMixed->Write();
 550   hPPMixed->Write(); hNNMixed->Write();*/
 551   gHistBalanceFunction->Write();
 552   if(listBFShuffled) gHistBalanceFunctionShuffled->Write();
 553   if(listBFMixed) {
 554     gHistBalanceFunctionMixed->Write();
 555     gHistBalanceFunctionSubtracted->Write();
 556   }
 557   fOutput->Close();
 558 }
 559
 560 //____________________________________________________________//
 561 void fitbalanceFunction(Int_t gCentrality = 1,
 562                         Double_t psiMin = -0.5, Double_t psiMax = 3.5,
 563                         Double_t ptTriggerMin = -1.,
 564                         Double_t ptTriggerMax = -1.,
 565                         Double_t ptAssociatedMin = -1.,
 566                         Double_t ptAssociatedMax = -1.,
 567                         TH2D *gHist,
 568                         Bool_t k2pMethod = kFALSE,
 569                         TString eventClass="EventPlane") {
 570   //balancing charges: [1]*TMath::Exp(-0.5*TMath::Power(((x - [3])/[2]),2)-0.5*TMath::Power(((y - [5])/[4]),2))
 571   //short range correlations: [6]*TMath::Exp(-0.5*TMath::Power(((x - [8])/[7]),2)-0.5*TMath::Power(((y - [10])/[9]),2))
 572   cout<<"FITTING FUNCTION"<<endl;
 573
 574   TF2 *gFitFunction = new TF2("gFitFunction","[0] + [1]*TMath::Exp(-0.5*TMath::Power(((x - [3])/[2]),2)-0.5*TMath::Power(((y - [5])/[4]),2)) + [6]*TMath::Exp(-0.5*TMath::Power(((x - [8])/[7]),2)-0.5*TMath::Power(((y - [10])/[9]),2))",-1.2,1.2,-TMath::Pi()/2.,3.*TMath::Pi()/2.);
 575   gFitFunction->SetName("gFitFunction");
 576
 577   //Normalization
 578   gFitFunction->SetParName(0,"N1");
 579   gFitFunction->SetParameter(0,1.0);
 580
 581   //2D balance function
 582   gFitFunction->SetParName(1,"N_{BF}");
 583   gFitFunction->SetParameter(1,1.0);
 584   gFitFunction->SetParLimits(1, 0., 100.);
 585   gFitFunction->SetParName(2,"Sigma_{BF}(delta eta)");
 586   gFitFunction->SetParameter(2,0.6);
 587   gFitFunction->SetParLimits(2, 0., 1.);
 588   gFitFunction->SetParName(3,"Mean_{BF}(delta eta)");
 589   gFitFunction->SetParameter(3,0.0);
 590   gFitFunction->SetParLimits(3, -0.2, 0.2);
 591   gFitFunction->SetParName(4,"Sigma_{BF}(delta phi)");
 592   gFitFunction->SetParameter(4,0.6);
 593   gFitFunction->SetParLimits(4, 0., 1.);
 594   gFitFunction->SetParName(5,"Mean_{BF}(delta phi)");
 595   gFitFunction->SetParameter(5,0.0);
 596   gFitFunction->SetParLimits(5, -0.2, 0.2);
 597
 598   //short range structure
 599   gFitFunction->SetParName(6,"N_{SR}");
 600   gFitFunction->SetParameter(6,5.0);
 601   gFitFunction->SetParLimits(6, 0., 100.);
 602   gFitFunction->SetParName(7,"Sigma_{SR}(delta eta)");
 603   gFitFunction->SetParameter(7,0.01);
 604   gFitFunction->SetParLimits(7, 0.0, 0.1);
 605   gFitFunction->SetParName(8,"Mean_{SR}(delta eta)");
 606   gFitFunction->SetParameter(8,0.0);
 607   gFitFunction->SetParLimits(8, -0.01, 0.01);
 608   gFitFunction->SetParName(9,"Sigma_{SR}(delta phi)");
 609   gFitFunction->SetParameter(9,0.01);
 610   gFitFunction->SetParLimits(9, 0.0, 0.1);
 611   gFitFunction->SetParName(10,"Mean_{SR}(delta phi)");
 612   gFitFunction->SetParameter(10,0.0);
 613   gFitFunction->SetParLimits(10, -0.01, 0.01);
 614
 615
 616   //Cloning the histogram
 617   TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
 618   gHistResidual->SetName("gHistResidual");
 619   gHistResidual->Sumw2();
 620
 621   //Fitting the 2D bf
 622   for(Int_t iAttempt = 0; iAttempt < 10; iAttempt++) {
 623     gHist->Fit("gFitFunction","nm");
 624     for(Int_t iParam = 0; iParam < 11; iParam++)
 625       gFitFunction->SetParameter(iParam,gFitFunction->GetParameter(iParam));
 626   }
 627   cout<<"======================================================"<<endl;
 628   cout<<"Fit chi2/ndf: "<<gFitFunction->GetChisquare()/gFitFunction->GetNDF()<<" - chi2: "<<gFitFunction->GetChisquare()<<" - ndf: "<<gFitFunction->GetNDF()<<endl;
 629   cout<<"======================================================"<<endl;
 630
 631   //Getting the residual
 632   gHistResidual->Add(gFitFunction,-1);
 633
 634   //Write to output file
 635   TString newFileName = "balanceFunctionFit2D.";
 636   if(eventClass == "Centrality"){
 637     newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
 638     newFileName += ".PsiAll.PttFrom";
 639   }
 640   else if(eventClass == "Multiplicity"){
 641     newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
 642     newFileName += ".PsiAll.PttFrom";
 643   }
 644   else{ // "EventPlane" (default)
 645     newFileName += "Centrality";
 646     newFileName += gCentrality; newFileName += ".Psi";
 647     if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
 648     else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
 649     else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
 650     else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
 651     else newFileName += "All.PttFrom";
 652   }
 653   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
 654   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
 655   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
 656   newFileName += Form("%.1f",ptAssociatedMax);
 657   if(k2pMethod) newFileName += "_2pMethod";
 658
 659   newFileName += "_";
 660   newFileName += Form("%.1f",psiMin);
 661   newFileName += "-";
 662   newFileName += Form("%.1f",psiMax);
 663   newFileName += ".root";
 664
 665   TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
 666   gHist->Write();
 667   gHistResidual->Write();
 668   gFitFunction->Write();
 669   newFile->Close();
 670 }
 671
 672 //____________________________________________________________//
 673 void drawBFPsi2D(const char* lhcPeriod = "LHC11h",
 674                  const char* gCentralityEstimator = "V0M",
 675                  Int_t gBit = 128,
 676                  const char* gEventPlaneEstimator = "VZERO",
 677                  Int_t gCentrality = 1,
 678                  Bool_t kShowShuffled = kFALSE,
 679                  Bool_t kShowMixed = kFALSE,
 680                  Double_t psiMin = -0.5, Double_t psiMax = 0.5,
 681                  Double_t ptTriggerMin = -1.,
 682                  Double_t ptTriggerMax = -1.,
 683                  Double_t ptAssociatedMin = -1.,
 684                  Double_t ptAssociatedMax = -1.,
 685                  Bool_t k2pMethod = kTRUE) {
 686   //Macro that draws the BF distributions for each centrality bin
 687   //for reaction plane dependent analysis
 688   //Author: Panos.Christakoglou@nikhef.nl
 689   TGaxis::SetMaxDigits(3);
 690
 691   //Get the input file
 692   TString filename = lhcPeriod;
 693   filename += "/Centrality"; filename += gCentralityEstimator;
 694   filename += "_Bit"; filename += gBit;
 695   filename += "_"; filename += gEventPlaneEstimator;
 696   filename +="/PttFrom";
 697   filename += Form("%.1f",ptTriggerMin); filename += "To";
 698   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
 699   filename += Form("%.1f",ptAssociatedMin); filename += "To";
 700   filename += Form("%.1f",ptAssociatedMax);
 701   filename += "/balanceFunction2D.Centrality";
 702   filename += gCentrality; filename += ".Psi";
 703   if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
 704   else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
 705   else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
 706   else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
 707   else filename += "All.PttFrom";
 708   filename += Form("%.1f",ptTriggerMin); filename += "To";
 709   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
 710   filename += Form("%.1f",ptAssociatedMin); filename += "To";
 711   filename += Form("%.1f",ptAssociatedMax);
 712   if(k2pMethod) filename += "_2pMethod";
 713
 714   filename += "_";
 715   filename += Form("%.1f",psiMin);
 716   filename += "-";
 717   filename += Form("%.1f",psiMax);
 718   filename += ".root";
 719
 720   //Open the file
 721   TFile *f = TFile::Open(filename.Data());
 722   if((!f)||(!f->IsOpen())) {
 723     Printf("The file %s is not found. Aborting...",filename);
 724     return listBF;
 725   }
 726   //f->ls();
 727
 728   //Raw balance function
 729   TH1D *gHistBalanceFunction = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunction"));
 730   gHistBalanceFunction->SetStats(kFALSE);
 731   gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
 732   gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
 733   gHistBalanceFunction->GetZaxis()->SetNdivisions(10);
 734   gHistBalanceFunction->GetXaxis()->SetTitleOffset(1.3);
 735   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
 736   gHistBalanceFunction->GetZaxis()->SetTitleOffset(1.3);
 737   gHistBalanceFunction->GetXaxis()->SetTitle("#Delta #eta");
 738   gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 739   gHistBalanceFunction->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
 740
 741   //Shuffled balance function
 742   if(kShowShuffled) {
 743     TH1D *gHistBalanceFunctionShuffled = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionShuffled"));
 744     gHistBalanceFunctionShuffled->SetStats(kFALSE);
 745     gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
 746     gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
 747     gHistBalanceFunctionShuffled->GetZaxis()->SetNdivisions(10);
 748     gHistBalanceFunctionShuffled->GetXaxis()->SetTitleOffset(1.3);
 749     gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
 750     gHistBalanceFunctionShuffled->GetZaxis()->SetTitleOffset(1.3);
 751     gHistBalanceFunctionShuffled->GetXaxis()->SetTitle("#Delta #eta");
 752     gHistBalanceFunctionShuffled->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 753     gHistBalanceFunctionShuffled->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
 754   }
 755
 756   //Mixed balance function
 757   if(kShowMixed) {
 758     TH1D *gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionMixed"));
 759     gHistBalanceFunctionMixed->SetStats(kFALSE);
 760     gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
 761     gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
 762     gHistBalanceFunctionMixed->GetZaxis()->SetNdivisions(10);
 763     gHistBalanceFunctionMixed->GetXaxis()->SetTitleOffset(1.3);
 764     gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
 765     gHistBalanceFunctionMixed->GetZaxis()->SetTitleOffset(1.3);
 766     gHistBalanceFunctionMixed->GetXaxis()->SetTitle("#Delta #eta");
 767     gHistBalanceFunctionMixed->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 768     gHistBalanceFunctionMixed->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
 769   }
 770
 771   //Subtracted balance function
 772   if(kShowMixed) {
 773     TH1D *gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionSubtracted"));
 774     gHistBalanceFunctionSubtracted->SetStats(kFALSE);
 775     gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
 776     gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
 777     gHistBalanceFunctionSubtracted->GetZaxis()->SetNdivisions(10);
 778     gHistBalanceFunctionSubtracted->GetXaxis()->SetTitleOffset(1.3);
 779     gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
 780     gHistBalanceFunctionSubtracted->GetZaxis()->SetTitleOffset(1.3);
 781     gHistBalanceFunctionSubtracted->GetXaxis()->SetTitle("#Delta #eta");
 782     gHistBalanceFunctionSubtracted->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 783     gHistBalanceFunctionSubtracted->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
 784   }
 785
 786   TString pngName;
 787
 788   TString centralityLatex = "Centrality: ";
 789   centralityLatex += centralityArray[gCentrality-1];
 790   centralityLatex += "%";
 791
 792   TString psiLatex;
 793   if((psiMin == -0.5)&&(psiMax == 0.5))
 794     psiLatex = " -7.5^{o} < #varphi - #Psi_{2} < 7.5^{o}";
 795   else if((psiMin == 0.5)&&(psiMax == 1.5))
 796     psiLatex = " 37.5^{o} < #varphi - #Psi_{2} < 52.5^{o}";
 797   else if((psiMin == 1.5)&&(psiMax == 2.5))
 798     psiLatex = " 82.5^{o} < #varphi - #Psi_{2} < 97.5^{o}";
 799   else
 800     psiLatex = " 0^{o} < #varphi - #Psi_{2} < 180^{o}";
 801
 802   TString pttLatex = Form("%.1f",ptTriggerMin);
 803   pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
 804   pttLatex += " GeV/c";
 805
 806   TString ptaLatex = Form("%.1f",ptAssociatedMin);
 807   ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
 808   ptaLatex += " GeV/c";
 809
 810   TLatex *latexInfo1 = new TLatex();
 811   latexInfo1->SetNDC();
 812   latexInfo1->SetTextSize(0.045);
 813   latexInfo1->SetTextColor(1);
 814
 815   //Draw the results
 816   TCanvas *c1 = new TCanvas("c1","Raw balance function 2D",0,0,600,500);
 817   c1->SetFillColor(10); c1->SetHighLightColor(10);
 818   c1->SetLeftMargin(0.17); c1->SetTopMargin(0.05);
 819   gHistBalanceFunction->SetTitle("");
 820   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.4);
 821   gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
 822   gHistBalanceFunction->GetXaxis()->SetRangeUser(-1.4,1.4);
 823   gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
 824   gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 825   gHistBalanceFunction->DrawCopy("lego2");
 826   gPad->SetTheta(30); // default is 30
 827   gPad->SetPhi(-60); // default is 30
 828   gPad->Update();
 829
 830   latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
 831   latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
 832   latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
 833   latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
 834
 835   TString pngName = "BalanceFunction2D.";
 836   pngName += "Centrality";
 837   pngName += gCentrality;
 838   if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
 839   else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
 840   else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
 841   else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.PttFrom";
 842   else pngName += "All.PttFrom";
 843   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
 844   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
 845   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
 846   pngName += Form("%.1f",ptAssociatedMax);
 847   if(k2pMethod) pngName += "_2pMethod";
 848   pngName += ".png";
 849
 850   c1->SaveAs(pngName.Data());
 851
 852   if(kShowShuffled) {
 853     TCanvas *c2 = new TCanvas("c2","Shuffled balance function 2D",100,100,600,500);
 854     c2->SetFillColor(10); c2->SetHighLightColor(10);
 855     c2->SetLeftMargin(0.17); c2->SetTopMargin(0.05);
 856     gHistBalanceFunctionShuffled->SetTitle("Shuffled events");
 857     gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.4);
 858     gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
 859     gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
 860     gHistBalanceFunctionShuffled->DrawCopy("lego2");
 861     gPad->SetTheta(30); // default is 30
 862     gPad->SetPhi(-60); // default is 30
 863     gPad->Update();
 864
 865     latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
 866     latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
 867     latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
 868     latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
 869   }
 870
 871   if(kShowMixed) {
 872     TCanvas *c3 = new TCanvas("c3","Mixed balance function 2D",200,200,600,500);
 873     c3->SetFillColor(10); c3->SetHighLightColor(10);
 874     c3->SetLeftMargin(0.17); c3->SetTopMargin(0.05);
 875     gHistBalanceFunctionMixed->SetTitle("Mixed events");
 876     gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.4);
 877     gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
 878     gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
 879     gHistBalanceFunctionMixed->DrawCopy("lego2");
 880     gPad->SetTheta(30); // default is 30
 881     gPad->SetPhi(-60); // default is 30
 882     gPad->Update();
 883
 884     latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
 885     latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
 886     latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
 887     latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
 888   }
 889
 890   if(kShowMixed) {
 891     TCanvas *c4 = new TCanvas("c4","Subtracted balance function 2D",300,300,600,500);
 892     c4->SetFillColor(10); c4->SetHighLightColor(10);
 893     c4->SetLeftMargin(0.17); c4->SetTopMargin(0.05);
 894     gHistBalanceFunctionSubtracted->SetTitle("Subtracted balance function");
 895     gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.4);
 896     gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
 897     gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
 898     gHistBalanceFunctionSubtracted->DrawCopy("lego2");
 899     gPad->SetTheta(30); // default is 30
 900     gPad->SetPhi(-60); // default is 30
 901     gPad->Update();
 902
 903     latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
 904     latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
 905     latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
 906     latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
 907   }
 908 }
 909
 910 //____________________________________________________________//
 911 void drawProjections(TH2D *gHistBalanceFunction2D = 0x0,
 912                      Bool_t kProjectInEta = kFALSE,
 913                      Int_t binMin = 1,
 914                      Int_t binMax = 80,
 915                      Int_t gCentrality = 1,
 916                      Double_t psiMin = -0.5,
 917                      Double_t psiMax = 3.5,
 918                      Double_t ptTriggerMin = -1.,
 919                      Double_t ptTriggerMax = -1.,
 920                      Double_t ptAssociatedMin = -1.,
 921                      Double_t ptAssociatedMax = -1.,
 922                      Bool_t k2pMethod = kTRUE,
 923                      TString eventClass = "Centrality",
 924                      Bool_t bRootMoments = kFALSE) {
 925   //Macro that draws the balance functions PROJECTIONS
 926   //for each centrality bin for the different pT of trigger and
 927   //associated particles
 928   TGaxis::SetMaxDigits(3);
 929
 930   //first we need some libraries
 931   gSystem->Load("libTree");
 932   gSystem->Load("libGeom");
 933   gSystem->Load("libVMC");
 934   gSystem->Load("libXMLIO");
 935   gSystem->Load("libPhysics");
 936
 937   gSystem->Load("libSTEERBase");
 938   gSystem->Load("libESD");
 939   gSystem->Load("libAOD");
 940
 941   gSystem->Load("libANALYSIS.so");
 942   gSystem->Load("libANALYSISalice.so");
 943   gSystem->Load("libEventMixing.so");
 944   gSystem->Load("libCORRFW.so");
 945   gSystem->Load("libPWGTools.so");
 946   gSystem->Load("libPWGCFebye.so");
 947
 948   gStyle->SetOptStat(0);
 949
 950   //Get the input file
 951   TString filename = "balanceFunction2D.";
 952   if(eventClass == "Centrality"){
 953     filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
 954     filename += ".PsiAll.PttFrom";
 955   }
 956   else if(eventClass == "Multiplicity"){
 957     filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
 958     filename += ".PsiAll.PttFrom";
 959   }
 960   else{ // "EventPlane" (default)
 961     filename += "Centrality";
 962     filename += gCentrality; filename += ".Psi";
 963     if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
 964     else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
 965     else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
 966     else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
 967     else filename += "All.PttFrom";
 968   }
 969   filename += Form("%.1f",ptTriggerMin); filename += "To";
 970   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
 971   filename += Form("%.1f",ptAssociatedMin); filename += "To";
 972   filename += Form("%.1f",ptAssociatedMax);
 973   if(k2pMethod) filename += "_2pMethod";
 974
 975   filename += "_";
 976   filename += Form("%.1f",psiMin);
 977   filename += "-";
 978   filename += Form("%.1f",psiMax);
 979   filename += ".root";
 980
 981   //Open the file
 982   TFile *f = 0x0;
 983   if(!gHistBalanceFunction2D) {
 984     TFile::Open(filename.Data());
 985     if((!f)||(!f->IsOpen())) {
 986       Printf("The file %s is not found. Aborting...",filename);
 987       return listBF;
 988     }
 989     //f->ls();
 990   }
 991
 992   //Latex
 993   TString centralityLatex = "Centrality: ";
 994   centralityLatex += centralityArray[gCentrality-1];
 995   centralityLatex += "%";
 996
 997   TString psiLatex;
 998   if((psiMin == -0.5)&&(psiMax == 0.5))
 999     psiLatex = " -7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o}";
1000   else if((psiMin == 0.5)&&(psiMax == 1.5))
1001     psiLatex = " 37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o}";
1002   else if((psiMin == 1.5)&&(psiMax == 2.5))
1003     psiLatex = " 82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o}";
1004   else
1005     psiLatex = " 0^{o} < #varphi^{t} - #Psi_{2} < 180^{o}";
1006
1007   TString pttLatex = Form("%.1f",ptTriggerMin);
1008   pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
1009   pttLatex += " GeV/c";
1010
1011   TString ptaLatex = Form("%.1f",ptAssociatedMin);
1012   ptaLatex += " < p_{T,assoc} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1013   ptaLatex += " GeV/c";
1014
1015   TLatex *latexInfo1 = new TLatex();
1016   latexInfo1->SetNDC();
1017   latexInfo1->SetTextSize(0.045);
1018   latexInfo1->SetTextColor(1);
1019
1020
1021   //============================================================//
1022   //Get subtracted and mixed balance function
1023   TH2D *gHistBalanceFunctionSubtracted2D = 0x0;
1024   TH2D *gHistBalanceFunctionMixed2D      = 0x0;
1025   if(!gHistBalanceFunction2D) {
1026     gHistBalanceFunctionSubtracted2D = (TH2D*)f->Get("gHistBalanceFunctionSubtracted");
1027     gHistBalanceFunctionMixed2D      = (TH2D*)f->Get("gHistBalanceFunctionMixed");
1028   }
1029   else {
1030     gHistBalanceFunctionSubtracted2D = dynamic_cast<TH2D*>(gHistBalanceFunction2D->Clone());
1031     gHistBalanceFunctionMixed2D = dynamic_cast<TH2D*>(gHistBalanceFunction2D->Clone());
1032   }
1033
1034   TH1D *gHistBalanceFunctionSubtracted = NULL;
1035   TH1D *gHistBalanceFunctionMixed      = NULL;
1036
1037   TH1D *gHistBalanceFunctionSubtracted_scale = NULL;
1038   TH1D *gHistBalanceFunctionMixed_scale      = NULL;
1039
1040   if(kProjectInEta){
1041     gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionX("gHistBalanceFunctionSubtractedEta",binMin,binMax));
1042     gHistBalanceFunctionSubtracted->Scale(gHistBalanceFunctionSubtracted2D->GetYaxis()->GetBinWidth(1));   // to remove normalization to phi bin width
1043     gHistBalanceFunctionMixed      = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionX("gHistBalanceFunctionMixedEta",binMin,binMax));
1044     gHistBalanceFunctionMixed->Scale(gHistBalanceFunctionMixed2D->GetYaxis()->GetBinWidth(1));   // to remove normalization to phi bin width
1045     gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#eta)");
1046     gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#eta)");
1047   }
1048   else{
1049     gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionY("gHistBalanceFunctionSubtractedPhi",binMin,binMax));
1050     gHistBalanceFunctionSubtracted->Scale(gHistBalanceFunctionSubtracted2D->GetXaxis()->GetBinWidth(1));   // to remove normalization to eta bin width
1051     gHistBalanceFunctionMixed      = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionY("gHistBalanceFunctionMixedPhi",binMin,binMax));
1052     gHistBalanceFunctionMixed->Scale(gHistBalanceFunctionMixed2D->GetXaxis()->GetBinWidth(1));   // to remove normalization to eta bin width
1053     gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#varphi)");
1054     gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#varphi)");
1055   }
1056
1057   gHistBalanceFunctionSubtracted->SetMarkerStyle(20);
1058   gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
1059
1060   gHistBalanceFunctionMixed->SetMarkerStyle(25);
1061
1062   TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
1063   c1->SetFillColor(10);
1064   c1->SetHighLightColor(10);
1065   c1->SetLeftMargin(0.15);
1066   gHistBalanceFunctionSubtracted->DrawCopy("E");
1067   gHistBalanceFunctionMixed->DrawCopy("ESAME");
1068
1069   legend = new TLegend(0.18,0.62,0.45,0.82,"","brNDC");
1070   legend->SetTextSize(0.045);
1071   legend->SetTextFont(42);
1072   legend->SetBorderSize(0);
1073   legend->SetFillStyle(0);
1074   legend->SetFillColor(10);
1075   legend->SetMargin(0.25);
1076   legend->SetShadowColor(10);
1077   legend->AddEntry(gHistBalanceFunctionSubtracted,"Data","lp");
1078   legend->AddEntry(gHistBalanceFunctionMixed,"Mixed data","lp");
1079   legend->Draw();
1080
1081
1082   TString meanLatex, rmsLatex, skewnessLatex, kurtosisLatex;
1083
1084   if(bRootMoments){
1085     meanLatex = "#mu = ";
1086     meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMean());
1087     meanLatex += " #pm ";
1088     meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMeanError());
1089
1090     rmsLatex = "#sigma = ";
1091     rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMS());
1092     rmsLatex += " #pm ";
1093     rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMSError());
1094
1095     skewnessLatex = "S = ";
1096     skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(1));
1097     skewnessLatex += " #pm ";
1098     skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(11));
1099
1100     kurtosisLatex = "K = ";
1101     kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(1));
1102     kurtosisLatex += " #pm ";
1103     kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(11));
1104
1105     Printf("Mean: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetMean(),gHistBalanceFunctionSubtracted->GetMeanError());
1106     Printf("RMS: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetRMS(),gHistBalanceFunctionSubtracted->GetRMSError());
1107     Printf("Skeweness: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetSkewness(1),gHistBalanceFunctionSubtracted->GetSkewness(11));
1108     Printf("Kurtosis: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetKurtosis(1),gHistBalanceFunctionSubtracted->GetKurtosis(11));
1109
1110
1111     // store in txt files
1112     TString meanFileName = filename;
1113     if(kProjectInEta)
1114       meanFileName= "deltaEtaProjection_Mean.txt";
1115     //meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
1116     else
1117       meanFileName = "deltaPhiProjection_Mean.txt";
1118       //meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
1119     ofstream fileMean(meanFileName.Data(),ios::app);
1120     fileMean << " " << gHistBalanceFunctionSubtracted->GetMean() << " " <<gHistBalanceFunctionSubtracted->GetMeanError()<<endl;
1121     fileMean.close();
1122
1123     TString rmsFileName = filename;
1124     if(kProjectInEta)
1125       rmsFileName = "deltaEtaProjection_Rms.txt";
1126       //rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
1127     else
1128       rmsFileName = "deltaPhiProjection_Rms.txt");
1129       //rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
1130     ofstream fileRms(rmsFileName.Data(),ios::app);
1131     fileRms << " " << gHistBalanceFunctionSubtracted->GetRMS() << " " <<gHistBalanceFunctionSubtracted->GetRMSError()<<endl;
1132     fileRms.close();
1133
1134     TString skewnessFileName = filename;
1135     if(kProjectInEta)
1136       skewnessFileName = "deltaEtaProjection_Skewness.txt";
1137       //skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
1138     else
1139       skewnessFileName = "deltaPhiProjection_Skewness.txt";
1140     //skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
1141     ofstream fileSkewness(skewnessFileName.Data(),ios::app);
1142     fileSkewness << " " << gHistBalanceFunctionSubtracted->GetSkewness(1) << " " <<gHistBalanceFunctionSubtracted->GetSkewness(11)<<endl;
1143     fileSkewness.close();
1144
1145     TString kurtosisFileName = filename;
1146     if(kProjectInEta)
1147       kurtosisFileName = "deltaEtaProjection_Kurtosis.txt";
1148       //kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
1149     else
1150       kurtosisFileName = "deltaPhiProjection_Kurtosis.txt";
1151       //kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
1152     ofstream fileKurtosis(kurtosisFileName.Data(),ios::app);
1153     fileKurtosis << " " << gHistBalanceFunctionSubtracted->GetKurtosis(1) << " " <<gHistBalanceFunctionSubtracted->GetKurtosis(11)<<endl;
1154     fileKurtosis.close();
1155   }
1156   // calculate the moments by hand
1157   else{
1158
1159     Double_t meanAnalytical, meanAnalyticalError;
1160     Double_t sigmaAnalytical, sigmaAnalyticalError;
1161     Double_t skewnessAnalytical, skewnessAnalyticalError;
1162     Double_t kurtosisAnalytical, kurtosisAnalyticalError;
1163
1164     Int_t gDeltaEtaPhi = 2;
1165     if(kProjectInEta) gDeltaEtaPhi = 1;
1166
1167     AliBalancePsi *bHelper = new AliBalancePsi;
1168     bHelper->GetMomentsAnalytical(gDeltaEtaPhi,gHistBalanceFunctionSubtracted,meanAnalytical,meanAnalyticalError,sigmaAnalytical,sigmaAnalyticalError,skewnessAnalytical,skewnessAnalyticalError,kurtosisAnalytical,kurtosisAnalyticalError);
1169
1170     meanLatex = "#mu = ";
1171     meanLatex += Form("%.3f",meanAnalytical);
1172     meanLatex += " #pm ";
1173     meanLatex += Form("%.3f",meanAnalyticalError);
1174
1175     rmsLatex = "#sigma = ";
1176     rmsLatex += Form("%.3f",sigmaAnalytical);
1177     rmsLatex += " #pm ";
1178     rmsLatex += Form("%.3f",sigmaAnalyticalError);
1179
1180     skewnessLatex = "S = ";
1181     skewnessLatex += Form("%.3f",skewnessAnalytical);
1182     skewnessLatex += " #pm ";
1183     skewnessLatex += Form("%.3f",skewnessAnalyticalError);
1184
1185     kurtosisLatex = "K = ";
1186     kurtosisLatex += Form("%.3f",kurtosisAnalytical);
1187     kurtosisLatex += " #pm ";
1188     kurtosisLatex += Form("%.3f",kurtosisAnalyticalError);
1189     Printf("Mean: %lf - Error: %lf",meanAnalytical, meanAnalyticalError);
1190     Printf("Sigma: %lf - Error: %lf",sigmaAnalytical, sigmaAnalyticalError);
1191     Printf("Skeweness: %lf - Error: %lf",skewnessAnalytical, skewnessAnalyticalError);
1192     Printf("Kurtosis: %lf - Error: %lf",kurtosisAnalytical, kurtosisAnalyticalError);
1193
1194     // store in txt files
1195     TString meanFileName = filename;
1196     if(kProjectInEta)
1197       meanFileName = "deltaEtaProjection_Mean.txt";
1198       //meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
1199     else
1200       meanFileName = "deltaPhiProjection_Mean.txt";
1201       //meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
1202     ofstream fileMean(meanFileName.Data(),ios::app);
1203     fileMean << " " << meanAnalytical << " " <<meanAnalyticalError<<endl;
1204     fileMean.close();
1205
1206     TString rmsFileName = filename;
1207     if(kProjectInEta)
1208       rmsFileName = "deltaEtaProjection_Rms.txt";
1209       //rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
1210     else
1211       rmsFileName = "deltaPhiProjection_Rms.txt";
1212 //rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
1213     ofstream fileRms(rmsFileName.Data(),ios::app);
1214     fileRms << " " << sigmaAnalytical << " " <<sigmaAnalyticalError<<endl;
1215     fileRms.close();
1216
1217     TString skewnessFileName = filename;
1218     if(kProjectInEta)
1219       skewnessFileName = "deltaEtaProjection_Skewness.txt";
1220       //skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
1221     else
1222       skewnessFileName = "deltaPhiProjection_Skewness.txt";
1223       //skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
1224     ofstream fileSkewness(skewnessFileName.Data(),ios::app);
1225     fileSkewness << " " << skewnessAnalytical << " " <<skewnessAnalyticalError<<endl;
1226     fileSkewness.close();
1227
1228     TString kurtosisFileName = filename;
1229     if(kProjectInEta)
1230       kurtosisFileName = "deltaEtaProjection_Kurtosis.txt";
1231       //kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
1232     else
1233       kurtosisFileName = "deltaPhiProjection_Kurtosis.txt";
1234       //kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
1235     ofstream fileKurtosis(kurtosisFileName.Data(),ios::app);
1236     fileKurtosis << " " << kurtosisAnalytical << " " <<kurtosisAnalyticalError<<endl;
1237     fileKurtosis.close();
1238   }
1239
1240
1241
1242   TCanvas *c2 = new TCanvas("c2","",600,0,600,500);
1243   c2->SetFillColor(10);
1244   c2->SetHighLightColor(10);
1245   c2->SetLeftMargin(0.15);
1246   gHistBalanceFunctionSubtracted->DrawCopy("E");
1247
1248   TLatex *latex = new TLatex();
1249   latex->SetNDC();
1250   latex->SetTextSize(0.035);
1251   latex->SetTextColor(1);
1252   latex->DrawLatex(0.64,0.85,meanLatex.Data());
1253   latex->DrawLatex(0.64,0.81,rmsLatex.Data());
1254   latex->DrawLatex(0.64,0.77,skewnessLatex.Data());
1255   latex->DrawLatex(0.64,0.73,kurtosisLatex.Data());
1256
1257   TString pngName = filename;
1258   if(kProjectInEta) pngName.ReplaceAll(".root","_DeltaEtaProjection.png");
1259   else              pngName.ReplaceAll(".root","_DeltaPhiProjection.png");
1260
1261   c2->SaveAs(pngName.Data());
1262
1263   TString outFileName = filename;
1264   if(kProjectInEta) outFileName.ReplaceAll(".root","_DeltaEtaProjection.root");
1265   else              outFileName.ReplaceAll(".root","_DeltaPhiProjection.root");
1266   TFile *fProjection = TFile::Open(outFileName.Data(),"recreate");
1267   gHistBalanceFunctionSubtracted->Write();
1268   gHistBalanceFunctionMixed->Write();
1269   fProjection->Close();
1270 }
1271
1272 //____________________________________________________________//
1273 void drawBFPsi2DFromCorrelationFunctions(const char* lhcPeriod = "LHC10h",
1274                                          Int_t gTrainNumber = 64,
1275                                          const char* gCentralityEstimator = "V0M",
1276                                          Int_t gBit = 128,
1277                                          const char* gEventPlaneEstimator = "VZERO",
1278                                          Int_t gCentrality = 1,
1279                                          Double_t psiMin = -0.5, Double_t psiMax = 3.5,
1280                                          Double_t vertexZMin = -10.,
1281                                          Double_t vertexZMax = 10.,
1282                                          Double_t ptTriggerMin = -1.,
1283                                          Double_t ptTriggerMax = -1.,
1284                                          Double_t ptAssociatedMin = -1.,
1285                                          Double_t ptAssociatedMax = -1.,
1286                                          TString eventClass = "Multiplicity"
1287 ) {
1288   //Macro that draws the BF distributions for each centrality bin
1289   //for reaction plane dependent analysis
1290   //Author: Panos.Christakoglou@nikhef.nl
1291   TGaxis::SetMaxDigits(3);
1292   gStyle->SetPalette(55,0);
1293
1294   //Get the input file
1295   TString filename = lhcPeriod;
1296   if(lhcPeriod != ""){
1297     //filename += "/Train"; filename += gTrainNumber;
1298     filename +="/PttFrom";
1299     filename += Form("%.1f",ptTriggerMin); filename += "To";
1300     filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1301     filename += Form("%.1f",ptAssociatedMin); filename += "To";
1302     filename += Form("%.1f",ptAssociatedMax);
1303     filename += "/correlationFunction.";
1304   }
1305   else{
1306     filename += "correlationFunction.";
1307   }
1308   if(eventClass == "Centrality"){
1309     filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1310     filename += ".PsiAll.PttFrom";
1311   }
1312   else if(eventClass == "Multiplicity"){
1313     filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1314     filename += ".PsiAll.PttFrom";
1315   }
1316   else{ // "EventPlane" (default)
1317     filename += "Centrality";
1318     filename += gCentrality; filename += ".Psi";
1319     if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
1320     else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
1321     else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
1322     else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
1323     else filename += "All.PttFrom";
1324   }
1325   filename += Form("%.1f",ptTriggerMin); filename += "To";
1326   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1327   filename += Form("%.1f",ptAssociatedMin); filename += "To";
1328   filename += Form("%.1f",ptAssociatedMax);
1329   filename += "_";
1330   filename += Form("%.1f",psiMin);
1331   filename += "-";
1332   filename += Form("%.1f",psiMax);
1333   filename += ".root";
1334
1335   //Open the file
1336   TFile *f = TFile::Open(filename.Data());
1337   if((!f)||(!f->IsOpen())) {
1338     Printf("The file %s is not found. Aborting...",filename);
1339     return listBF;
1340   }
1341   //f->ls();
1342
1343   TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
1344   if(!gHistPN) return;
1345   TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
1346   if(!gHistNP) return;
1347   TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
1348   if(!gHistPP) return;
1349   TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
1350   if(!gHistNN) return;
1351
1352   gHistPN->Sumw2();
1353   gHistPP->Sumw2();
1354   gHistPN->Add(gHistPP,-1);
1355   gHistNP->Sumw2();
1356   gHistNN->Sumw2();
1357   gHistNP->Add(gHistNN,-1);
1358   gHistPN->Add(gHistNP);
1359   gHistPN->Scale(0.5);
1360   TH2D *gHistBalanceFunction2D = dynamic_cast<TH2D *>(gHistPN->Clone());
1361   gHistBalanceFunction2D->SetStats(kFALSE);
1362   gHistBalanceFunction2D->GetXaxis()->SetTitle("#Delta#eta");
1363   gHistBalanceFunction2D->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1364   gHistBalanceFunction2D->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
1365
1366   //Draw the results
1367   TCanvas *c0 = new TCanvas("c0","Balance function 2D",0,0,600,500);
1368   c0->SetFillColor(10); c0->SetHighLightColor(10);
1369   c0->SetLeftMargin(0.17); c0->SetTopMargin(0.05);
1370   gHistBalanceFunction2D->SetTitle("");
1371   gHistBalanceFunction2D->GetZaxis()->SetTitleOffset(1.4);
1372   gHistBalanceFunction2D->GetZaxis()->SetNdivisions(10);
1373   gHistBalanceFunction2D->GetYaxis()->SetTitleOffset(1.4);
1374   gHistBalanceFunction2D->GetYaxis()->SetNdivisions(10);
1375   gHistBalanceFunction2D->GetXaxis()->SetRangeUser(-1.4,1.4);
1376   gHistBalanceFunction2D->GetXaxis()->SetNdivisions(10);
1377   gHistBalanceFunction2D->DrawCopy("lego2");
1378   gPad->SetTheta(30); // default is 30
1379   gPad->SetPhi(-60); // default is 30
1380   gPad->Update();
1381
1382   TString multLatex = Form("Multiplicity: %.1f - %.1f",psiMin,psiMax);
1383
1384   TString pttLatex = Form("%.1f",ptTriggerMin);
1385   pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
1386   pttLatex += " GeV/c";
1387
1388   TString ptaLatex = Form("%.1f",ptAssociatedMin);
1389   ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1390   ptaLatex += " GeV/c";
1391
1392   TLatex *latexInfo1 = new TLatex();
1393   latexInfo1->SetNDC();
1394   latexInfo1->SetTextSize(0.045);
1395   latexInfo1->SetTextColor(1);
1396   latexInfo1->DrawLatex(0.54,0.88,multLatex.Data());
1397   latexInfo1->DrawLatex(0.54,0.82,pttLatex.Data());
1398   latexInfo1->DrawLatex(0.54,0.76,ptaLatex.Data());
1399
1400   TString pngName = "BalanceFunction2D.";
1401   pngName += Form("%s: %.1f - %.1f",eventClass.Data(),psiMin,psiMax);
1402   pngName += ".PttFrom";
1403   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1404   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1405   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1406   pngName += Form("%.1f",ptAssociatedMax);
1407   pngName += ".png";
1408
1409   c0->SaveAs(pngName.Data());
1410
1411   drawProjections(gHistBalanceFunction2D,
1412                   kTRUE,
1413                   1,36,
1414                   gCentrality,
1415                   psiMin,psiMax,
1416                   ptTriggerMin,ptTriggerMax,
1417                   ptAssociatedMin,ptAssociatedMax,
1418                   kTRUE,
1419                   eventClass,
1420                   kFALSE);
1421
1422   drawProjections(gHistBalanceFunction2D,
1423                   kFALSE,
1424                   1,80,
1425                   gCentrality,
1426                   psiMin,psiMax,
1427                   ptTriggerMin,ptTriggerMax,
1428                   ptAssociatedMin,ptAssociatedMax,
1429                   kTRUE,
1430                   eventClass.Data(),
1431                   kFALSE);
1432 }