PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C

   1 const Int_t numberOfCentralityBins = 12;
   2 TString centralityArray[numberOfCentralityBins] = {"0-80","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-100","0-1","1-2","2-3"};
   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                               Bool_t bToy = kFALSE)
  23 {
  24   //Macro that draws the BF distributions for each centrality bin
  25   //for reaction plane dependent analysis
  26   //Author: Panos.Christakoglou@nikhef.nl
  27   //Load the PWG2ebye library
  28   gSystem->Load("libANALYSIS.so");
  29   gSystem->Load("libANALYSISalice.so");
  30   gSystem->Load("libEventMixing.so");
  31   gSystem->Load("libCORRFW.so");
  32   gSystem->Load("libPWGTools.so");
  33   gSystem->Load("libPWGCFebye.so");
  34
  35   //gROOT->LoadMacro("~/SetPlotStyle.C");
  36   //SetPlotStyle();
  37   gStyle->SetPalette(1,0);
  38
  39   //Prepare the objects and return them
  40   TList *listBF = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,0,bToy);
  41   TList *listBFShuffled = NULL;
  42   if(kShowShuffled) listBFShuffled = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,1,bToy);
  43   TList *listBFMixed = NULL;
  44   if(kShowMixed) listBFMixed = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,2,bToy);
  45   if(!listBF) {
  46     Printf("The TList object was not created");
  47     return;
  48   }
  49   else
  50     draw(listBF,listBFShuffled,listBFMixed,gCentrality,gCentralityEstimator,
  51          psiMin,psiMax,vertexZMin,vertexZMax,
  52          ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,
  53          kUseVzBinning,k2pMethod,eventClass);
  54 }
  55
  56 //______________________________________________________//
  57 TList *GetListOfObjects(const char* filename,
  58                         Int_t gCentrality,
  59                         Int_t gBit,
  60                         const char *gCentralityEstimator,
  61                         Int_t kData = 1,
  62                         Bool_t bToy = kFALSE) {
  63   //Get the TList objects (QA, bf, bf shuffled)
  64   TList *listBF = 0x0;
  65
  66   //Open the file
  67   TFile *f = TFile::Open(filename,"UPDATE");
  68   if((!f)||(!f->IsOpen())) {
  69     Printf("The file %s is not found. Aborting...",filename);
  70     return listBF;
  71   }
  72   //f->ls();
  73
  74   TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
  75   if(!dir) {
  76     Printf("The TDirectoryFile is not found. Aborting...",filename);
  77     return listBF;
  78   }
  79   //dir->ls();
  80
  81   TString listBFName;
  82   if(kData == 0) {
  83     //cout<<"no shuffling - no mixing"<<endl;
  84     listBFName = "listBFPsi";
  85   }
  86   else if(kData == 1) {
  87     //cout<<"shuffling - no mixing"<<endl;
  88     listBFName = "listBFPsiShuffled";
  89   }
  90   else if(kData == 2) {
  91     //cout<<"no shuffling - mixing"<<endl;
  92     listBFName = "listBFPsiMixed";
  93   }
  94
  95   // different list names in case of toy model
  96   if(!bToy){
  97     listBFName += "_";
  98     listBFName += centralityArray[gCentrality-1];
  99     if(gBit > -1) {
 100       listBFName += "_Bit"; listBFName += gBit; }
 101     if(gCentralityEstimator) {
 102       listBFName += "_"; listBFName += gCentralityEstimator;}
 103   }
 104   else{
 105     listBFName.ReplaceAll("Psi","");
 106   }
 107
 108   // histograms were already retrieved (in first iteration)
 109   if(dir->Get(Form("%s_histograms",listBFName.Data()))){
 110     //cout<<"second iteration"<<endl;
 111     listBF = dynamic_cast<TList *>(dir->Get(Form("%s_histograms",listBFName.Data())));
 112   }
 113
 114   // histograms were not yet retrieved (this is the first iteration)
 115   else{
 116
 117     listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
 118     cout<<"======================================================="<<endl;
 119     cout<<"List name: "<<listBF->GetName()<<endl;
 120     //listBF->ls();
 121
 122     //Get the histograms
 123     TString histoName;
 124     if(kData == 0)
 125       histoName = "fHistP";
 126     else if(kData == 1)
 127       histoName = "fHistP_shuffle";
 128     else if(kData == 2)
 129       histoName = "fHistP";
 130     if(gCentralityEstimator) histoName += gCentralityEstimator;
 131     AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 132     if(!fHistP) {
 133       Printf("fHistP %s not found!!!",histoName.Data());
 134       break;
 135     }
 136     fHistP->FillParent(); fHistP->DeleteContainers();
 137
 138     if(kData == 0)
 139       histoName = "fHistN";
 140     if(kData == 1)
 141       histoName = "fHistN_shuffle";
 142     if(kData == 2)
 143       histoName = "fHistN";
 144     if(gCentralityEstimator) histoName += gCentralityEstimator;
 145     AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 146     if(!fHistN) {
 147       Printf("fHistN %s not found!!!",histoName.Data());
 148       break;
 149     }
 150     fHistN->FillParent(); fHistN->DeleteContainers();
 151
 152     if(kData == 0)
 153       histoName = "fHistPN";
 154     if(kData == 1)
 155       histoName = "fHistPN_shuffle";
 156     if(kData == 2)
 157       histoName = "fHistPN";
 158     if(gCentralityEstimator) histoName += gCentralityEstimator;
 159     AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 160     if(!fHistPN) {
 161       Printf("fHistPN %s not found!!!",histoName.Data());
 162       break;
 163     }
 164     fHistPN->FillParent(); fHistPN->DeleteContainers();
 165
 166     if(kData == 0)
 167       histoName = "fHistNP";
 168     if(kData == 1)
 169       histoName = "fHistNP_shuffle";
 170     if(kData == 2)
 171       histoName = "fHistNP";
 172     if(gCentralityEstimator) histoName += gCentralityEstimator;
 173     AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 174     if(!fHistNP) {
 175       Printf("fHistNP %s not found!!!",histoName.Data());
 176       break;
 177     }
 178     fHistNP->FillParent(); fHistNP->DeleteContainers();
 179
 180     if(kData == 0)
 181       histoName = "fHistPP";
 182     if(kData == 1)
 183       histoName = "fHistPP_shuffle";
 184     if(kData == 2)
 185       histoName = "fHistPP";
 186     if(gCentralityEstimator) histoName += gCentralityEstimator;
 187     AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 188     if(!fHistPP) {
 189       Printf("fHistPP %s not found!!!",histoName.Data());
 190       break;
 191     }
 192     fHistPP->FillParent(); fHistPP->DeleteContainers();
 193
 194     if(kData == 0)
 195       histoName = "fHistNN";
 196     if(kData == 1)
 197       histoName = "fHistNN_shuffle";
 198     if(kData == 2)
 199       histoName = "fHistNN";
 200     if(gCentralityEstimator) histoName += gCentralityEstimator;
 201     AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 202     if(!fHistNN) {
 203       Printf("fHistNN %s not found!!!",histoName.Data());
 204       break;
 205     }
 206     fHistNN->FillParent(); fHistNN->DeleteContainers();
 207
 208     dir->cd();
 209     listBF->Write(Form("%s_histograms",listBFName.Data()), TObject::kSingleKey);
 210
 211   }// first iteration
 212
 213   f->Close();
 214
 215   return listBF;
 216 }
 217
 218 //______________________________________________________//
 219 void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
 220           Int_t gCentrality, const char* gCentralityEstimator,
 221           Double_t psiMin, Double_t psiMax,
 222           Double_t vertexZMin,
 223           Double_t vertexZMax,
 224           Double_t ptTriggerMin, Double_t ptTriggerMax,
 225           Double_t ptAssociatedMin, Double_t ptAssociatedMax,
 226           Bool_t kUseVzBinning=kFALSE,
 227           Bool_t k2pMethod = kFALSE, TString eventClass) {
 228   //balance function
 229   AliTHn *hP = NULL;
 230   AliTHn *hN = NULL;
 231   AliTHn *hPN = NULL;
 232   AliTHn *hNP = NULL;
 233   AliTHn *hPP = NULL;
 234   AliTHn *hNN = NULL;
 235   //listBF->ls();
 236   //Printf("=================");
 237   TString histoName = "fHistP";
 238   if(gCentralityEstimator) histoName += gCentralityEstimator;
 239   hP = (AliTHn*) listBF->FindObject(histoName.Data());
 240   hP->SetName("gHistP");
 241   histoName = "fHistN";
 242   if(gCentralityEstimator) histoName += gCentralityEstimator;
 243   hN = (AliTHn*) listBF->FindObject(histoName.Data());
 244   hN->SetName("gHistN");
 245   histoName = "fHistPN";
 246   if(gCentralityEstimator) histoName += gCentralityEstimator;
 247   hPN = (AliTHn*) listBF->FindObject(histoName.Data());
 248   hPN->SetName("gHistPN");
 249   histoName = "fHistNP";
 250   if(gCentralityEstimator) histoName += gCentralityEstimator;
 251   hNP = (AliTHn*) listBF->FindObject(histoName.Data());
 252   hNP->SetName("gHistNP");
 253   histoName = "fHistPP";
 254   if(gCentralityEstimator) histoName += gCentralityEstimator;
 255   hPP = (AliTHn*) listBF->FindObject(histoName.Data());
 256   hPP->SetName("gHistPP");
 257   histoName = "fHistNN";
 258   if(gCentralityEstimator) histoName += gCentralityEstimator;
 259   hNN = (AliTHn*) listBF->FindObject(histoName.Data());
 260   hNN->SetName("gHistNN");
 261
 262   AliBalancePsi *b = new AliBalancePsi();
 263   b->SetEventClass(eventClass);
 264   b->SetHistNp(hP);
 265   b->SetHistNn(hN);
 266   b->SetHistNpn(hPN);
 267   b->SetHistNnp(hNP);
 268   b->SetHistNpp(hPP);
 269   b->SetHistNnn(hNN);
 270   if(kUseVzBinning) b->SetVertexZBinning(kTRUE);
 271
 272
 273   //balance function shuffling
 274   AliTHn *hPShuffled = NULL;
 275   AliTHn *hNShuffled = NULL;
 276   AliTHn *hPNShuffled = NULL;
 277   AliTHn *hNPShuffled = NULL;
 278   AliTHn *hPPShuffled = NULL;
 279   AliTHn *hNNShuffled = NULL;
 280   if(listBFShuffled) {
 281     //listBFShuffled->ls();
 282     histoName = "fHistP_shuffle";
 283     if(gCentralityEstimator) histoName += gCentralityEstimator;
 284     hPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
 285     hPShuffled->SetName("gHistPShuffled");
 286     histoName = "fHistN_shuffle";
 287     if(gCentralityEstimator) histoName += gCentralityEstimator;
 288     hNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
 289     hNShuffled->SetName("gHistNShuffled");
 290     histoName = "fHistPN_shuffle";
 291     if(gCentralityEstimator) histoName += gCentralityEstimator;
 292     hPNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
 293     hPNShuffled->SetName("gHistPNShuffled");
 294     histoName = "fHistNP_shuffle";
 295     if(gCentralityEstimator) histoName += gCentralityEstimator;
 296     hNPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
 297     hNPShuffled->SetName("gHistNPShuffled");
 298     histoName = "fHistPP_shuffle";
 299     if(gCentralityEstimator) histoName += gCentralityEstimator;
 300     hPPShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
 301     hPPShuffled->SetName("gHistPPShuffled");
 302     histoName = "fHistNN_shuffle";
 303     if(gCentralityEstimator) histoName += gCentralityEstimator;
 304     hNNShuffled = (AliTHn*) listBFShuffled->FindObject(histoName.Data());
 305     hNNShuffled->SetName("gHistNNShuffled");
 306
 307     AliBalancePsi *bShuffled = new AliBalancePsi();
 308     bShuffled->SetEventClass(eventClass);
 309     bShuffled->SetHistNp(hPShuffled);
 310     bShuffled->SetHistNn(hNShuffled);
 311     bShuffled->SetHistNpn(hPNShuffled);
 312     bShuffled->SetHistNnp(hNPShuffled);
 313     bShuffled->SetHistNpp(hPPShuffled);
 314     bShuffled->SetHistNnn(hNNShuffled);
 315   if(kUseVzBinning) bShuffled->SetVertexZBinning(kTRUE);
 316
 317   }
 318
 319   //balance function mixing
 320   AliTHn *hPMixed = NULL;
 321   AliTHn *hNMixed = NULL;
 322   AliTHn *hPNMixed = NULL;
 323   AliTHn *hNPMixed = NULL;
 324   AliTHn *hPPMixed = NULL;
 325   AliTHn *hNNMixed = NULL;
 326
 327   if(listBFMixed) {
 328     //listBFMixed->ls();
 329     histoName = "fHistP";
 330     if(gCentralityEstimator) histoName += gCentralityEstimator;
 331     hPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
 332     hPMixed->SetName("gHistPMixed");
 333     histoName = "fHistN";
 334     if(gCentralityEstimator) histoName += gCentralityEstimator;
 335     hNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
 336     hNMixed->SetName("gHistNMixed");
 337     histoName = "fHistPN";
 338     if(gCentralityEstimator) histoName += gCentralityEstimator;
 339     hPNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
 340     hPNMixed->SetName("gHistPNMixed");
 341     histoName = "fHistNP";
 342     if(gCentralityEstimator) histoName += gCentralityEstimator;
 343     hNPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
 344     histoName = "fHistNP";
 345     if(gCentralityEstimator) histoName += gCentralityEstimator;
 346     hNPMixed->SetName("gHistNPMixed");
 347     histoName = "fHistPP";
 348     if(gCentralityEstimator) histoName += gCentralityEstimator;
 349     hPPMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
 350     hPPMixed->SetName("gHistPPMixed");
 351     histoName = "fHistNN";
 352     if(gCentralityEstimator) histoName += gCentralityEstimator;
 353     hNNMixed = (AliTHn*) listBFMixed->FindObject(histoName.Data());
 354     hNNMixed->SetName("gHistNNMixed");
 355
 356     AliBalancePsi *bMixed = new AliBalancePsi();
 357     bMixed->SetEventClass(eventClass);
 358     bMixed->SetHistNp(hPMixed);
 359     bMixed->SetHistNn(hNMixed);
 360     bMixed->SetHistNpn(hPNMixed);
 361     bMixed->SetHistNnp(hNPMixed);
 362     bMixed->SetHistNpp(hPPMixed);
 363     bMixed->SetHistNnn(hNNMixed);
 364     if(kUseVzBinning) bMixed->SetVertexZBinning(kTRUE);
 365
 366   }
 367
 368   TH2D *gHistBalanceFunction;
 369   TH2D *gHistBalanceFunctionSubtracted;
 370   TH2D *gHistBalanceFunctionShuffled;
 371   TH2D *gHistBalanceFunctionMixed;
 372   TString histoTitle, pngName;
 373
 374   if(eventClass == "Centrality"){
 375     histoTitle = "Centrality: ";
 376     histoTitle += psiMin;
 377     histoTitle += " - ";
 378     histoTitle += psiMax;
 379     histoTitle += " % ";
 380     histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 381   }
 382   else if(eventClass == "Multiplicity"){
 383     histoTitle = "Multiplicity: ";
 384     histoTitle += psiMin;
 385     histoTitle += " - ";
 386     histoTitle += psiMax;
 387     histoTitle += " tracks";
 388     histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 389   }
 390   else{ // "EventPlane" (default)
 391     histoTitle = "Centrality: ";
 392     histoTitle += centralityArray[gCentrality-1];
 393     histoTitle += "%";
 394     if((psiMin == -0.5)&&(psiMax == 0.5))
 395       histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 396     else if((psiMin == 0.5)&&(psiMax == 1.5))
 397       histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 398     else if((psiMin == 1.5)&&(psiMax == 2.5))
 399       histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 400     else
 401       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 402   }
 403
 404   if(k2pMethod)
 405     if(bMixed)
 406       gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 407     else{
 408       cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
 409       return;
 410     }
 411   else
 412     gHistBalanceFunction = b->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 413   gHistBalanceFunction->SetTitle(histoTitle.Data());
 414   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
 415   gHistBalanceFunction->SetName("gHistBalanceFunction");
 416
 417   if(listBFShuffled) {
 418
 419     if(k2pMethod)
 420       if(bMixed)
 421         gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 422       else{
 423         cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
 424         return;
 425       }
 426     else
 427       gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 428     gHistBalanceFunctionShuffled->SetTitle(histoTitle.Data());
 429     gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
 430     gHistBalanceFunctionShuffled->SetName("gHistBalanceFunctionShuffled");
 431   }
 432
 433   if(listBFMixed) {
 434     if(k2pMethod)
 435       if(bMixed)
 436         gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi2pMethod(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 437       else{
 438         cerr<<"NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
 439         return;
 440       }
 441     else
 442       gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionDeltaEtaDeltaPhi(psiMin,psiMax,vertexZMin,vertexZMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 443     gHistBalanceFunctionMixed->SetTitle(histoTitle.Data());
 444     gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
 445     gHistBalanceFunctionMixed->SetName("gHistBalanceFunctionMixed");
 446
 447     gHistBalanceFunctionSubtracted = dynamic_cast<TH2D *>(gHistBalanceFunction->Clone());
 448     gHistBalanceFunctionSubtracted->Add(gHistBalanceFunctionMixed,-1);
 449     gHistBalanceFunctionSubtracted->SetTitle(histoTitle.Data());
 450     gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
 451     gHistBalanceFunctionSubtracted->SetName("gHistBalanceFunctionSubtracted");
 452   }
 453
 454   //Draw the results
 455   TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
 456   c1->SetFillColor(10);
 457   c1->SetHighLightColor(10);
 458   c1->SetLeftMargin(0.15);
 459   gHistBalanceFunction->DrawCopy("lego2");
 460   gPad->SetTheta(30); // default is 30
 461   gPad->SetPhi(-60); // default is 30
 462   gPad->Update();
 463   TCanvas *c1a = new TCanvas("c1a","",600,0,600,500);
 464   c1a->SetFillColor(10);
 465   c1a->SetHighLightColor(10);
 466   c1a->SetLeftMargin(0.15);
 467   gHistBalanceFunction->DrawCopy("colz");
 468
 469   if(listBFShuffled) {
 470     TCanvas *c2 = new TCanvas("c2","",100,100,600,500);
 471     c2->SetFillColor(10);
 472     c2->SetHighLightColor(10);
 473     c2->SetLeftMargin(0.15);
 474     gHistBalanceFunctionShuffled->DrawCopy("lego2");
 475     gPad->SetTheta(30); // default is 30
 476     gPad->SetPhi(-60); // default is 30
 477     gPad->Update();
 478     TCanvas *c2a = new TCanvas("c2a","",700,100,600,500);
 479     c2a->SetFillColor(10);
 480     c2a->SetHighLightColor(10);
 481     c2a->SetLeftMargin(0.15);
 482     gHistBalanceFunctionShuffled->DrawCopy("colz");
 483   }
 484
 485   if(listBFMixed) {
 486     TCanvas *c3 = new TCanvas("c3","",200,200,600,500);
 487     c3->SetFillColor(10);
 488     c3->SetHighLightColor(10);
 489     c3->SetLeftMargin(0.15);
 490     gHistBalanceFunctionMixed->DrawCopy("lego2");
 491     gPad->SetTheta(30); // default is 30
 492     gPad->SetPhi(-60); // default is 30
 493     gPad->Update();
 494     TCanvas *c3a = new TCanvas("c3a","",800,200,600,500);
 495     c3a->SetFillColor(10);
 496     c3a->SetHighLightColor(10);
 497     c3a->SetLeftMargin(0.15);
 498     gHistBalanceFunctionMixed->DrawCopy("colz");
 499
 500     TCanvas *c4 = new TCanvas("c4","",300,300,600,500);
 501     c4->SetFillColor(10);
 502     c4->SetHighLightColor(10);
 503     c4->SetLeftMargin(0.15);
 504     gHistBalanceFunctionSubtracted->DrawCopy("lego2");
 505     gPad->SetTheta(30); // default is 30
 506     gPad->SetPhi(-60); // default is 30
 507     gPad->Update();
 508     TCanvas *c4a = new TCanvas("c4a","",900,300,600,500);
 509     c4a->SetFillColor(10);
 510     c4a->SetHighLightColor(10);
 511     c4a->SetLeftMargin(0.15);
 512     gHistBalanceFunctionSubtracted->DrawCopy("colz");
 513
 514     fitbalanceFunction(gCentrality, psiMin , psiMax,
 515                        ptTriggerMin, ptTriggerMax,
 516                        ptAssociatedMin, ptAssociatedMax,
 517                        gHistBalanceFunctionSubtracted,k2pMethod, eventClass);
 518   }
 519
 520   TString newFileName = "balanceFunction2D.";
 521   if(eventClass == "Centrality"){
 522     newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
 523     newFileName += ".PsiAll.PttFrom";
 524   }
 525   else if(eventClass == "Multiplicity"){
 526     newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
 527     newFileName += ".PsiAll.PttFrom";
 528   }
 529   else{ // "EventPlane" (default)
 530     newFileName += "Centrality";
 531     newFileName += gCentrality; newFileName += ".Psi";
 532     if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
 533     else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
 534     else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
 535     else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
 536     else newFileName += "All.PttFrom";
 537   }
 538   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
 539   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
 540   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
 541   newFileName += Form("%.1f",ptAssociatedMax);
 542   if(k2pMethod) newFileName += "_2pMethod";
 543
 544   newFileName += "_";
 545   newFileName += Form("%.1f",psiMin);
 546   newFileName += "-";
 547   newFileName += Form("%.1f",psiMax);
 548   newFileName += ".root";
 549
 550   TFile *fOutput = new TFile(newFileName.Data(),"recreate");
 551   fOutput->cd();
 552   /*hP->Write(); hN->Write();
 553   hPN->Write(); hNP->Write();
 554   hPP->Write(); hNN->Write();
 555   hPShuffled->Write(); hNShuffled->Write();
 556   hPNShuffled->Write(); hNPShuffled->Write();
 557   hPPShuffled->Write(); hNNShuffled->Write();
 558   hPMixed->Write(); hNMixed->Write();
 559   hPNMixed->Write(); hNPMixed->Write();
 560   hPPMixed->Write(); hNNMixed->Write();*/
 561   gHistBalanceFunction->Write();
 562   if(listBFShuffled) gHistBalanceFunctionShuffled->Write();
 563   if(listBFMixed) {
 564     gHistBalanceFunctionMixed->Write();
 565     gHistBalanceFunctionSubtracted->Write();
 566   }
 567   fOutput->Close();
 568 }
 569
 570 //____________________________________________________________//
 571 void fitbalanceFunction(Int_t gCentrality = 1,
 572                         Double_t psiMin = -0.5, Double_t psiMax = 3.5,
 573                         Double_t ptTriggerMin = -1.,
 574                         Double_t ptTriggerMax = -1.,
 575                         Double_t ptAssociatedMin = -1.,
 576                         Double_t ptAssociatedMax = -1.,
 577                         TH2D *gHist,
 578                         Bool_t k2pMethod = kFALSE,
 579                         TString eventClass="EventPlane") {
 580   //balancing charges: [1]*TMath::Exp(-0.5*TMath::Power(((x - [3])/[2]),2)-0.5*TMath::Power(((y - [5])/[4]),2))
 581   //short range correlations: [6]*TMath::Exp(-0.5*TMath::Power(((x - [8])/[7]),2)-0.5*TMath::Power(((y - [10])/[9]),2))
 582   cout<<"FITTING FUNCTION"<<endl;
 583
 584   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.);
 585   gFitFunction->SetName("gFitFunction");
 586
 587   //Normalization
 588   gFitFunction->SetParName(0,"N1");
 589   gFitFunction->SetParameter(0,1.0);
 590
 591   //2D balance function
 592   gFitFunction->SetParName(1,"N_{BF}");
 593   gFitFunction->SetParameter(1,1.0);
 594   gFitFunction->SetParLimits(1, 0., 100.);
 595   gFitFunction->SetParName(2,"Sigma_{BF}(delta eta)");
 596   gFitFunction->SetParameter(2,0.6);
 597   gFitFunction->SetParLimits(2, 0., 1.);
 598   gFitFunction->SetParName(3,"Mean_{BF}(delta eta)");
 599   gFitFunction->SetParameter(3,0.0);
 600   gFitFunction->SetParLimits(3, -0.2, 0.2);
 601   gFitFunction->SetParName(4,"Sigma_{BF}(delta phi)");
 602   gFitFunction->SetParameter(4,0.6);
 603   gFitFunction->SetParLimits(4, 0., 1.);
 604   gFitFunction->SetParName(5,"Mean_{BF}(delta phi)");
 605   gFitFunction->SetParameter(5,0.0);
 606   gFitFunction->SetParLimits(5, -0.2, 0.2);
 607
 608   //short range structure
 609   gFitFunction->SetParName(6,"N_{SR}");
 610   gFitFunction->SetParameter(6,5.0);
 611   gFitFunction->SetParLimits(6, 0., 100.);
 612   gFitFunction->SetParName(7,"Sigma_{SR}(delta eta)");
 613   gFitFunction->SetParameter(7,0.01);
 614   gFitFunction->SetParLimits(7, 0.0, 0.1);
 615   gFitFunction->SetParName(8,"Mean_{SR}(delta eta)");
 616   gFitFunction->SetParameter(8,0.0);
 617   gFitFunction->SetParLimits(8, -0.01, 0.01);
 618   gFitFunction->SetParName(9,"Sigma_{SR}(delta phi)");
 619   gFitFunction->SetParameter(9,0.01);
 620   gFitFunction->SetParLimits(9, 0.0, 0.1);
 621   gFitFunction->SetParName(10,"Mean_{SR}(delta phi)");
 622   gFitFunction->SetParameter(10,0.0);
 623   gFitFunction->SetParLimits(10, -0.01, 0.01);
 624
 625
 626   //Cloning the histogram
 627   TH2D *gHistResidual = dynamic_cast<TH2D *>(gHist->Clone());
 628   gHistResidual->SetName("gHistResidual");
 629   gHistResidual->Sumw2();
 630
 631   //Fitting the 2D bf
 632   for(Int_t iAttempt = 0; iAttempt < 10; iAttempt++) {
 633     gHist->Fit("gFitFunction","nm");
 634     for(Int_t iParam = 0; iParam < 11; iParam++)
 635       gFitFunction->SetParameter(iParam,gFitFunction->GetParameter(iParam));
 636   }
 637   cout<<"======================================================"<<endl;
 638   cout<<"Fit chi2/ndf: "<<gFitFunction->GetChisquare()/gFitFunction->GetNDF()<<" - chi2: "<<gFitFunction->GetChisquare()<<" - ndf: "<<gFitFunction->GetNDF()<<endl;
 639   cout<<"======================================================"<<endl;
 640
 641   //Getting the residual
 642   gHistResidual->Add(gFitFunction,-1);
 643
 644   //Write to output file
 645   TString newFileName = "balanceFunctionFit2D.";
 646   if(eventClass == "Centrality"){
 647     newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
 648     newFileName += ".PsiAll.PttFrom";
 649   }
 650   else if(eventClass == "Multiplicity"){
 651     newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
 652     newFileName += ".PsiAll.PttFrom";
 653   }
 654   else{ // "EventPlane" (default)
 655     newFileName += "Centrality";
 656     newFileName += gCentrality; newFileName += ".Psi";
 657     if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
 658     else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
 659     else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
 660     else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
 661     else newFileName += "All.PttFrom";
 662   }
 663   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
 664   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
 665   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
 666   newFileName += Form("%.1f",ptAssociatedMax);
 667   if(k2pMethod) newFileName += "_2pMethod";
 668
 669   newFileName += "_";
 670   newFileName += Form("%.1f",psiMin);
 671   newFileName += "-";
 672   newFileName += Form("%.1f",psiMax);
 673   newFileName += ".root";
 674
 675   TFile *newFile = TFile::Open(newFileName.Data(),"recreate");
 676   gHist->Write();
 677   gHistResidual->Write();
 678   gFitFunction->Write();
 679   newFile->Close();
 680 }
 681
 682 //____________________________________________________________//
 683 void drawBFPsi2D(const char* lhcPeriod = "LHC11h",
 684                  const char* gCentralityEstimator = "V0M",
 685                  Int_t gBit = 128,
 686                  const char* gEventPlaneEstimator = "VZERO",
 687                  Int_t gCentrality = 1,
 688                  Bool_t kShowShuffled = kFALSE,
 689                  Bool_t kShowMixed = kFALSE,
 690                  Double_t psiMin = -0.5, Double_t psiMax = 0.5,
 691                  Double_t ptTriggerMin = -1.,
 692                  Double_t ptTriggerMax = -1.,
 693                  Double_t ptAssociatedMin = -1.,
 694                  Double_t ptAssociatedMax = -1.,
 695                  Bool_t k2pMethod = kTRUE) {
 696   //Macro that draws the BF distributions for each centrality bin
 697   //for reaction plane dependent analysis
 698   //Author: Panos.Christakoglou@nikhef.nl
 699   TGaxis::SetMaxDigits(3);
 700
 701   //Get the input file
 702   TString filename = lhcPeriod;
 703   filename += "/Centrality"; filename += gCentralityEstimator;
 704   filename += "_Bit"; filename += gBit;
 705   filename += "_"; filename += gEventPlaneEstimator;
 706   filename +="/PttFrom";
 707   filename += Form("%.1f",ptTriggerMin); filename += "To";
 708   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
 709   filename += Form("%.1f",ptAssociatedMin); filename += "To";
 710   filename += Form("%.1f",ptAssociatedMax);
 711   filename += "/balanceFunction2D.Centrality";
 712   filename += gCentrality; filename += ".Psi";
 713   if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
 714   else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
 715   else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
 716   else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
 717   else filename += "All.PttFrom";
 718   filename += Form("%.1f",ptTriggerMin); filename += "To";
 719   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
 720   filename += Form("%.1f",ptAssociatedMin); filename += "To";
 721   filename += Form("%.1f",ptAssociatedMax);
 722   if(k2pMethod) filename += "_2pMethod";
 723
 724   filename += "_";
 725   filename += Form("%.1f",psiMin);
 726   filename += "-";
 727   filename += Form("%.1f",psiMax);
 728   filename += ".root";
 729
 730   //Open the file
 731   TFile *f = TFile::Open(filename.Data());
 732   if((!f)||(!f->IsOpen())) {
 733     Printf("The file %s is not found. Aborting...",filename);
 734     return listBF;
 735   }
 736   //f->ls();
 737
 738   //Raw balance function
 739   TH1D *gHistBalanceFunction = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunction"));
 740   gHistBalanceFunction->SetStats(kFALSE);
 741   gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
 742   gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
 743   gHistBalanceFunction->GetZaxis()->SetNdivisions(10);
 744   gHistBalanceFunction->GetXaxis()->SetTitleOffset(1.3);
 745   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
 746   gHistBalanceFunction->GetZaxis()->SetTitleOffset(1.3);
 747   gHistBalanceFunction->GetXaxis()->SetTitle("#Delta #eta");
 748   gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 749   gHistBalanceFunction->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
 750
 751   //Shuffled balance function
 752   if(kShowShuffled) {
 753     TH1D *gHistBalanceFunctionShuffled = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionShuffled"));
 754     gHistBalanceFunctionShuffled->SetStats(kFALSE);
 755     gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
 756     gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
 757     gHistBalanceFunctionShuffled->GetZaxis()->SetNdivisions(10);
 758     gHistBalanceFunctionShuffled->GetXaxis()->SetTitleOffset(1.3);
 759     gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.3);
 760     gHistBalanceFunctionShuffled->GetZaxis()->SetTitleOffset(1.3);
 761     gHistBalanceFunctionShuffled->GetXaxis()->SetTitle("#Delta #eta");
 762     gHistBalanceFunctionShuffled->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 763     gHistBalanceFunctionShuffled->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
 764   }
 765
 766   //Mixed balance function
 767   if(kShowMixed) {
 768     TH1D *gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionMixed"));
 769     gHistBalanceFunctionMixed->SetStats(kFALSE);
 770     gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
 771     gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
 772     gHistBalanceFunctionMixed->GetZaxis()->SetNdivisions(10);
 773     gHistBalanceFunctionMixed->GetXaxis()->SetTitleOffset(1.3);
 774     gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.3);
 775     gHistBalanceFunctionMixed->GetZaxis()->SetTitleOffset(1.3);
 776     gHistBalanceFunctionMixed->GetXaxis()->SetTitle("#Delta #eta");
 777     gHistBalanceFunctionMixed->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 778     gHistBalanceFunctionMixed->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
 779   }
 780
 781   //Subtracted balance function
 782   if(kShowMixed) {
 783     TH1D *gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionSubtracted"));
 784     gHistBalanceFunctionSubtracted->SetStats(kFALSE);
 785     gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
 786     gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
 787     gHistBalanceFunctionSubtracted->GetZaxis()->SetNdivisions(10);
 788     gHistBalanceFunctionSubtracted->GetXaxis()->SetTitleOffset(1.3);
 789     gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
 790     gHistBalanceFunctionSubtracted->GetZaxis()->SetTitleOffset(1.3);
 791     gHistBalanceFunctionSubtracted->GetXaxis()->SetTitle("#Delta #eta");
 792     gHistBalanceFunctionSubtracted->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 793     gHistBalanceFunctionSubtracted->GetZaxis()->SetTitle("B(#Delta #eta, #Delta #varphi)");
 794   }
 795
 796   TString pngName;
 797
 798   TString centralityLatex = "Centrality: ";
 799   centralityLatex += centralityArray[gCentrality-1];
 800   centralityLatex += "%";
 801
 802   TString psiLatex;
 803   if((psiMin == -0.5)&&(psiMax == 0.5))
 804     psiLatex = " -7.5^{o} < #varphi - #Psi_{2} < 7.5^{o}";
 805   else if((psiMin == 0.5)&&(psiMax == 1.5))
 806     psiLatex = " 37.5^{o} < #varphi - #Psi_{2} < 52.5^{o}";
 807   else if((psiMin == 1.5)&&(psiMax == 2.5))
 808     psiLatex = " 82.5^{o} < #varphi - #Psi_{2} < 97.5^{o}";
 809   else
 810     psiLatex = " 0^{o} < #varphi - #Psi_{2} < 180^{o}";
 811
 812   TString pttLatex = Form("%.1f",ptTriggerMin);
 813   pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
 814   pttLatex += " GeV/c";
 815
 816   TString ptaLatex = Form("%.1f",ptAssociatedMin);
 817   ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
 818   ptaLatex += " GeV/c";
 819
 820   TLatex *latexInfo1 = new TLatex();
 821   latexInfo1->SetNDC();
 822   latexInfo1->SetTextSize(0.045);
 823   latexInfo1->SetTextColor(1);
 824
 825   //Draw the results
 826   TCanvas *c1 = new TCanvas("c1","Raw balance function 2D",0,0,600,500);
 827   c1->SetFillColor(10); c1->SetHighLightColor(10);
 828   c1->SetLeftMargin(0.17); c1->SetTopMargin(0.05);
 829   gHistBalanceFunction->SetTitle("");
 830   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.4);
 831   gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
 832   gHistBalanceFunction->GetXaxis()->SetRangeUser(-1.4,1.4);
 833   gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
 834   gHistBalanceFunction->GetYaxis()->SetTitle("#Delta #varphi (rad)");
 835   gHistBalanceFunction->DrawCopy("lego2");
 836   gPad->SetTheta(30); // default is 30
 837   gPad->SetPhi(-60); // default is 30
 838   gPad->Update();
 839
 840   latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
 841   latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
 842   latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
 843   latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
 844
 845   TString pngName = "BalanceFunction2D.";
 846   pngName += "Centrality";
 847   pngName += gCentrality;
 848   if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
 849   else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
 850   else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
 851   else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.PttFrom";
 852   else pngName += "All.PttFrom";
 853   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
 854   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
 855   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
 856   pngName += Form("%.1f",ptAssociatedMax);
 857   if(k2pMethod) pngName += "_2pMethod";
 858   pngName += ".png";
 859
 860   c1->SaveAs(pngName.Data());
 861
 862   if(kShowShuffled) {
 863     TCanvas *c2 = new TCanvas("c2","Shuffled balance function 2D",100,100,600,500);
 864     c2->SetFillColor(10); c2->SetHighLightColor(10);
 865     c2->SetLeftMargin(0.17); c2->SetTopMargin(0.05);
 866     gHistBalanceFunctionShuffled->SetTitle("Shuffled events");
 867     gHistBalanceFunctionShuffled->GetYaxis()->SetTitleOffset(1.4);
 868     gHistBalanceFunctionShuffled->GetYaxis()->SetNdivisions(10);
 869     gHistBalanceFunctionShuffled->GetXaxis()->SetNdivisions(10);
 870     gHistBalanceFunctionShuffled->DrawCopy("lego2");
 871     gPad->SetTheta(30); // default is 30
 872     gPad->SetPhi(-60); // default is 30
 873     gPad->Update();
 874
 875     latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
 876     latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
 877     latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
 878     latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
 879   }
 880
 881   if(kShowMixed) {
 882     TCanvas *c3 = new TCanvas("c3","Mixed balance function 2D",200,200,600,500);
 883     c3->SetFillColor(10); c3->SetHighLightColor(10);
 884     c3->SetLeftMargin(0.17); c3->SetTopMargin(0.05);
 885     gHistBalanceFunctionMixed->SetTitle("Mixed events");
 886     gHistBalanceFunctionMixed->GetYaxis()->SetTitleOffset(1.4);
 887     gHistBalanceFunctionMixed->GetYaxis()->SetNdivisions(10);
 888     gHistBalanceFunctionMixed->GetXaxis()->SetNdivisions(10);
 889     gHistBalanceFunctionMixed->DrawCopy("lego2");
 890     gPad->SetTheta(30); // default is 30
 891     gPad->SetPhi(-60); // default is 30
 892     gPad->Update();
 893
 894     latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
 895     latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
 896     latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
 897     latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
 898   }
 899
 900   if(kShowMixed) {
 901     TCanvas *c4 = new TCanvas("c4","Subtracted balance function 2D",300,300,600,500);
 902     c4->SetFillColor(10); c4->SetHighLightColor(10);
 903     c4->SetLeftMargin(0.17); c4->SetTopMargin(0.05);
 904     gHistBalanceFunctionSubtracted->SetTitle("Subtracted balance function");
 905     gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.4);
 906     gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
 907     gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
 908     gHistBalanceFunctionSubtracted->DrawCopy("lego2");
 909     gPad->SetTheta(30); // default is 30
 910     gPad->SetPhi(-60); // default is 30
 911     gPad->Update();
 912
 913     latexInfo1->DrawLatex(0.64,0.88,centralityLatex.Data());
 914     latexInfo1->DrawLatex(0.64,0.82,psiLatex.Data());
 915     latexInfo1->DrawLatex(0.64,0.76,pttLatex.Data());
 916     latexInfo1->DrawLatex(0.64,0.70,ptaLatex.Data());
 917   }
 918 }
 919
 920 //____________________________________________________________//
 921 void drawProjections(TH2D *gHistBalanceFunction2D = 0x0,
 922                      Bool_t kProjectInEta = kFALSE,
 923                      Int_t binMin = 1,
 924                      Int_t binMax = 80,
 925                      Int_t gCentrality = 1,
 926                      Double_t psiMin = -0.5,
 927                      Double_t psiMax = 3.5,
 928                      Double_t ptTriggerMin = -1.,
 929                      Double_t ptTriggerMax = -1.,
 930                      Double_t ptAssociatedMin = -1.,
 931                      Double_t ptAssociatedMax = -1.,
 932                      Bool_t k2pMethod = kTRUE,
 933                      TString eventClass = "Centrality",
 934                      Bool_t bRootMoments = kFALSE,
 935                      Bool_t kUseZYAM = kFALSE,
 936                      Bool_t bReduceRangeForMoments = kFALSE,
 937                      Bool_t bFWHM = kFALSE) {
 938   //Macro that draws the balance functions PROJECTIONS
 939   //for each centrality bin for the different pT of trigger and
 940   //associated particles
 941   TGaxis::SetMaxDigits(3);
 942
 943   //first we need some libraries
 944   gSystem->Load("libTree");
 945   gSystem->Load("libGeom");
 946   gSystem->Load("libVMC");
 947   gSystem->Load("libXMLIO");
 948   gSystem->Load("libPhysics");
 949
 950   gSystem->Load("libSTEERBase");
 951   gSystem->Load("libESD");
 952   gSystem->Load("libAOD");
 953
 954   gSystem->Load("libANALYSIS.so");
 955   gSystem->Load("libANALYSISalice.so");
 956   gSystem->Load("libEventMixing.so");
 957   gSystem->Load("libCORRFW.so");
 958   gSystem->Load("libPWGTools.so");
 959   gSystem->Load("libPWGCFebye.so");
 960
 961   gStyle->SetOptStat(0);
 962
 963   //Get the input file
 964   TString filename = "balanceFunction2D.";
 965   if(eventClass == "Centrality"){
 966     filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
 967     filename += ".PsiAll.PttFrom";
 968   }
 969   else if(eventClass == "Multiplicity"){
 970     filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
 971     filename += ".PsiAll.PttFrom";
 972   }
 973   else{ // "EventPlane" (default)
 974     filename += "Centrality";
 975     filename += gCentrality; filename += ".Psi";
 976     if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
 977     else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
 978     else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
 979     else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
 980     else filename += "All.PttFrom";
 981   }
 982   filename += Form("%.1f",ptTriggerMin); filename += "To";
 983   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
 984   filename += Form("%.1f",ptAssociatedMin); filename += "To";
 985   filename += Form("%.1f",ptAssociatedMax);
 986   if(k2pMethod) filename += "_2pMethod";
 987
 988   filename += "_";
 989   filename += Form("%.1f",psiMin);
 990   filename += "-";
 991   filename += Form("%.1f",psiMax);
 992   filename += ".root";
 993
 994   //Open the file
 995   TFile *f = 0x0;
 996   if(!gHistBalanceFunction2D) {
 997     TFile::Open(filename.Data());
 998     if((!f)||(!f->IsOpen())) {
 999       Printf("The file %s is not found. Aborting...",filename);
1000       return listBF;
1001     }
1002     //f->ls();
1003   }
1004
1005   //Latex
1006   TString centralityLatex = "Centrality: ";
1007   centralityLatex += centralityArray[gCentrality-1];
1008   centralityLatex += "%";
1009
1010   TString psiLatex;
1011   if((psiMin == -0.5)&&(psiMax == 0.5))
1012     psiLatex = " -7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o}";
1013   else if((psiMin == 0.5)&&(psiMax == 1.5))
1014     psiLatex = " 37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o}";
1015   else if((psiMin == 1.5)&&(psiMax == 2.5))
1016     psiLatex = " 82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o}";
1017   else
1018     psiLatex = " 0^{o} < #varphi^{t} - #Psi_{2} < 180^{o}";
1019
1020   TString pttLatex = Form("%.1f",ptTriggerMin);
1021   pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
1022   pttLatex += " GeV/c";
1023
1024   TString ptaLatex = Form("%.1f",ptAssociatedMin);
1025   ptaLatex += " < p_{T,assoc} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1026   ptaLatex += " GeV/c";
1027
1028   TLatex *latexInfo1 = new TLatex();
1029   latexInfo1->SetNDC();
1030   latexInfo1->SetTextSize(0.045);
1031   latexInfo1->SetTextColor(1);
1032
1033
1034   //============================================================//
1035   //Get subtracted and mixed balance function
1036   TH2D *gHistBalanceFunctionSubtracted2D = 0x0;
1037   TH2D *gHistBalanceFunctionMixed2D      = 0x0;
1038   if(!gHistBalanceFunction2D) {
1039     gHistBalanceFunctionSubtracted2D = (TH2D*)f->Get("gHistBalanceFunctionSubtracted");
1040     gHistBalanceFunctionMixed2D      = (TH2D*)f->Get("gHistBalanceFunctionMixed");
1041   }
1042   else {
1043     gHistBalanceFunctionSubtracted2D = dynamic_cast<TH2D*>(gHistBalanceFunction2D->Clone());
1044     gHistBalanceFunctionMixed2D = dynamic_cast<TH2D*>(gHistBalanceFunction2D->Clone());
1045   }
1046
1047   TH1D *gHistBalanceFunctionSubtracted = NULL;
1048   TH1D *gHistBalanceFunctionMixed      = NULL;
1049
1050   TH1D *gHistBalanceFunctionSubtracted_scale = NULL;
1051   TH1D *gHistBalanceFunctionMixed_scale      = NULL;
1052
1053   if(kProjectInEta){
1054     gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionX("gHistBalanceFunctionSubtractedEta",binMin,binMax));
1055     gHistBalanceFunctionSubtracted->Scale(gHistBalanceFunctionSubtracted2D->GetYaxis()->GetBinWidth(1));   // to remove normalization to phi bin width
1056     gHistBalanceFunctionMixed      = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionX("gHistBalanceFunctionMixedEta",binMin,binMax));
1057     gHistBalanceFunctionMixed->Scale(gHistBalanceFunctionMixed2D->GetYaxis()->GetBinWidth(1));   // to remove normalization to phi bin width
1058     gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#eta)");
1059     gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#eta)");
1060   }
1061   else{
1062     gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionY("gHistBalanceFunctionSubtractedPhi",binMin,binMax));
1063     gHistBalanceFunctionSubtracted->Scale(gHistBalanceFunctionSubtracted2D->GetXaxis()->GetBinWidth(1));   // to remove normalization to eta bin width
1064     gHistBalanceFunctionMixed      = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionY("gHistBalanceFunctionMixedPhi",binMin,binMax));
1065     gHistBalanceFunctionMixed->Scale(gHistBalanceFunctionMixed2D->GetXaxis()->GetBinWidth(1));   // to remove normalization to eta bin width
1066     gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#varphi)");
1067     gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#varphi)");
1068   }
1069
1070   gHistBalanceFunctionSubtracted->SetMarkerStyle(20);
1071   gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
1072
1073   gHistBalanceFunctionMixed->SetMarkerStyle(25);
1074
1075   TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
1076   c1->SetFillColor(10);
1077   c1->SetHighLightColor(10);
1078   c1->SetLeftMargin(0.15);
1079   gHistBalanceFunctionSubtracted->DrawCopy("E");
1080   gHistBalanceFunctionMixed->DrawCopy("ESAME");
1081
1082   legend = new TLegend(0.18,0.62,0.45,0.82,"","brNDC");
1083   legend->SetTextSize(0.045);
1084   legend->SetTextFont(42);
1085   legend->SetBorderSize(0);
1086   legend->SetFillStyle(0);
1087   legend->SetFillColor(10);
1088   legend->SetMargin(0.25);
1089   legend->SetShadowColor(10);
1090   legend->AddEntry(gHistBalanceFunctionSubtracted,"Data","lp");
1091   legend->AddEntry(gHistBalanceFunctionMixed,"Mixed data","lp");
1092   legend->Draw();
1093
1094
1095   TString meanLatex, rmsLatex, skewnessLatex, kurtosisLatex;
1096
1097   if(bRootMoments){
1098
1099     // need to restrict to near side in case of dphi
1100     if(!kProjectInEta) gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-TMath::Pi()/2,TMath::Pi()/2);
1101
1102     if(bReduceRangeForMoments){
1103
1104       // default (for 2<pT,assoc<3<pT,trig<4)
1105       Double_t rangeReduced = 1.0;
1106
1107       //for 3<pT,assoc<8<pT,trig<15
1108       if(ptTriggerMax>10){
1109         rangeReduced = 0.3;
1110       }
1111
1112       cout<<"Use reduced range = "<<rangeReduced<<endl;
1113       gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-rangeReduced,rangeReduced);
1114     }
1115
1116     meanLatex = "#mu = ";
1117     meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMean());
1118     meanLatex += " #pm ";
1119     meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMeanError());
1120
1121     rmsLatex = "#sigma = ";
1122     rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMS());
1123     rmsLatex += " #pm ";
1124     rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMSError());
1125
1126     skewnessLatex = "S = ";
1127     skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(1));
1128     skewnessLatex += " #pm ";
1129     skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(11));
1130
1131     kurtosisLatex = "K = ";
1132     kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(1));
1133     kurtosisLatex += " #pm ";
1134     kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(11));
1135
1136     Printf("ROOT Mean: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetMean(),gHistBalanceFunctionSubtracted->GetMeanError());
1137     Printf("ROOT RMS: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetRMS(),gHistBalanceFunctionSubtracted->GetRMSError());
1138     Printf("ROOT Skeweness: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetSkewness(1),gHistBalanceFunctionSubtracted->GetSkewness(11));
1139     Printf("ROOT Kurtosis: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetKurtosis(1),gHistBalanceFunctionSubtracted->GetKurtosis(11));
1140
1141
1142     // store in txt files
1143     TString meanFileName = filename;
1144     if(kProjectInEta)
1145       meanFileName= "deltaEtaProjection_Mean.txt";
1146     //meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
1147     else
1148       meanFileName = "deltaPhiProjection_Mean.txt";
1149       //meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
1150     ofstream fileMean(meanFileName.Data(),ios::app);
1151     fileMean << " " << gHistBalanceFunctionSubtracted->GetMean() << " " <<gHistBalanceFunctionSubtracted->GetMeanError()<<endl;
1152     fileMean.close();
1153
1154     TString rmsFileName = filename;
1155     if(kProjectInEta)
1156       rmsFileName = "deltaEtaProjection_Rms.txt";
1157       //rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
1158     else
1159       rmsFileName = "deltaPhiProjection_Rms.txt";
1160       //rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
1161     ofstream fileRms(rmsFileName.Data(),ios::app);
1162     fileRms << " " << gHistBalanceFunctionSubtracted->GetRMS() << " " <<gHistBalanceFunctionSubtracted->GetRMSError()<<endl;
1163     fileRms.close();
1164
1165     TString skewnessFileName = filename;
1166     if(kProjectInEta)
1167       skewnessFileName = "deltaEtaProjection_Skewness.txt";
1168       //skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
1169     else
1170       skewnessFileName = "deltaPhiProjection_Skewness.txt";
1171     //skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
1172     ofstream fileSkewness(skewnessFileName.Data(),ios::app);
1173     fileSkewness << " " << gHistBalanceFunctionSubtracted->GetSkewness(1) << " " <<gHistBalanceFunctionSubtracted->GetSkewness(11)<<endl;
1174     fileSkewness.close();
1175
1176     TString kurtosisFileName = filename;
1177     if(kProjectInEta)
1178       kurtosisFileName = "deltaEtaProjection_Kurtosis.txt";
1179       //kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
1180     else
1181       kurtosisFileName = "deltaPhiProjection_Kurtosis.txt";
1182       //kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
1183     ofstream fileKurtosis(kurtosisFileName.Data(),ios::app);
1184     fileKurtosis << " " << gHistBalanceFunctionSubtracted->GetKurtosis(1) << " " <<gHistBalanceFunctionSubtracted->GetKurtosis(11)<<endl;
1185     fileKurtosis.close();
1186
1187     if(!kProjectInEta) gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-TMath::Pi()/2,3.*TMath::Pi()/2);
1188     else               gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-1.6,1.6);
1189   }
1190   // calculate the moments by hand
1191   else{
1192
1193     Double_t meanAnalytical, meanAnalyticalError;
1194     Double_t sigmaAnalytical, sigmaAnalyticalError;
1195     Double_t skewnessAnalytical, skewnessAnalyticalError;
1196     Double_t kurtosisAnalytical, kurtosisAnalyticalError;
1197
1198     Int_t gDeltaEtaPhi = 2;
1199     if(kProjectInEta) gDeltaEtaPhi = 1;
1200
1201     AliBalancePsi *bHelper = new AliBalancePsi;
1202     bHelper->GetMomentsAnalytical(gDeltaEtaPhi,gHistBalanceFunctionSubtracted,kUseZYAM,meanAnalytical,meanAnalyticalError,sigmaAnalytical,sigmaAnalyticalError,skewnessAnalytical,skewnessAnalyticalError,kurtosisAnalytical,kurtosisAnalyticalError);
1203
1204     meanLatex = "#mu = ";
1205     meanLatex += Form("%.3f",meanAnalytical);
1206     meanLatex += " #pm ";
1207     meanLatex += Form("%.3f",meanAnalyticalError);
1208
1209     rmsLatex = "#sigma = ";
1210     rmsLatex += Form("%.3f",sigmaAnalytical);
1211     rmsLatex += " #pm ";
1212     rmsLatex += Form("%.3f",sigmaAnalyticalError);
1213
1214     skewnessLatex = "S = ";
1215     skewnessLatex += Form("%.3f",skewnessAnalytical);
1216     skewnessLatex += " #pm ";
1217     skewnessLatex += Form("%.3f",skewnessAnalyticalError);
1218
1219     kurtosisLatex = "K = ";
1220     kurtosisLatex += Form("%.3f",kurtosisAnalytical);
1221     kurtosisLatex += " #pm ";
1222     kurtosisLatex += Form("%.3f",kurtosisAnalyticalError);
1223     Printf("Mean: %lf - Error: %lf",meanAnalytical, meanAnalyticalError);
1224     Printf("Sigma: %lf - Error: %lf",sigmaAnalytical, sigmaAnalyticalError);
1225     Printf("Skeweness: %lf - Error: %lf",skewnessAnalytical, skewnessAnalyticalError);
1226     Printf("Kurtosis: %lf - Error: %lf",kurtosisAnalytical, kurtosisAnalyticalError);
1227
1228     // store in txt files
1229     TString meanFileName = filename;
1230     if(kProjectInEta)
1231       meanFileName = "deltaEtaProjection_Mean.txt";
1232       //meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
1233     else
1234       meanFileName = "deltaPhiProjection_Mean.txt";
1235       //meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
1236     ofstream fileMean(meanFileName.Data(),ios::app);
1237     fileMean << " " << meanAnalytical << " " <<meanAnalyticalError<<endl;
1238     fileMean.close();
1239
1240     TString rmsFileName = filename;
1241     if(kProjectInEta)
1242       rmsFileName = "deltaEtaProjection_Rms.txt";
1243       //rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
1244     else
1245       rmsFileName = "deltaPhiProjection_Rms.txt";
1246 //rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
1247     ofstream fileRms(rmsFileName.Data(),ios::app);
1248     fileRms << " " << sigmaAnalytical << " " <<sigmaAnalyticalError<<endl;
1249     fileRms.close();
1250
1251     TString skewnessFileName = filename;
1252     if(kProjectInEta)
1253       skewnessFileName = "deltaEtaProjection_Skewness.txt";
1254       //skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
1255     else
1256       skewnessFileName = "deltaPhiProjection_Skewness.txt";
1257       //skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
1258     ofstream fileSkewness(skewnessFileName.Data(),ios::app);
1259     fileSkewness << " " << skewnessAnalytical << " " <<skewnessAnalyticalError<<endl;
1260     fileSkewness.close();
1261
1262     TString kurtosisFileName = filename;
1263     if(kProjectInEta)
1264       kurtosisFileName = "deltaEtaProjection_Kurtosis.txt";
1265       //kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
1266     else
1267       kurtosisFileName = "deltaPhiProjection_Kurtosis.txt";
1268       //kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
1269     ofstream fileKurtosis(kurtosisFileName.Data(),ios::app);
1270     fileKurtosis << " " << kurtosisAnalytical << " " <<kurtosisAnalyticalError<<endl;
1271     fileKurtosis.close();
1272   }
1273
1274   // Weighted mean as calculated for 1D analysis
1275   Double_t weightedMean, weightedMeanError;
1276   if(bReduceRangeForMoments){
1277     GetWeightedMean1D(gHistBalanceFunctionSubtracted,kProjectInEta,1,-1,weightedMean,weightedMeanError,rangeReduced);
1278   }
1279   else{
1280     GetWeightedMean1D(gHistBalanceFunctionSubtracted,kProjectInEta,1,-1,weightedMean,weightedMeanError,-1.);
1281   }
1282   Printf("Weighted Mean: %lf - Error: %lf",weightedMean, weightedMeanError);
1283
1284   // store in txt files
1285   TString weightedMeanFileName = filename;
1286   if(kProjectInEta)
1287     weightedMeanFileName = "deltaEtaProjection_WeightedMean.txt";
1288     //weightedMeanFileName.ReplaceAll(".root","_DeltaEtaProjection_WeightedMean.txt");
1289   else
1290     weightedMeanFileName = "deltaPhiProjection_WeightedMean.txt";
1291     //weightedMeanFileName.ReplaceAll(".root","_DeltaPhiProjection_WeightedMean.txt");
1292   ofstream fileWeightedMean(weightedMeanFileName.Data(),ios::app);
1293   fileWeightedMean << " " << weightedMean << " " <<weightedMeanError<<endl;
1294   fileWeightedMean.close();
1295
1296   TCanvas *c2 = new TCanvas("c2","",600,0,600,500);
1297   c2->SetFillColor(10);
1298   c2->SetHighLightColor(10);
1299   c2->SetLeftMargin(0.15);
1300   gHistBalanceFunctionSubtracted->DrawCopy("E");
1301
1302   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++//
1303   Double_t fwhm_spline = 0.;
1304   Double_t fwhmError = 0.;
1305   if (bFWHM){
1306     AliBalancePsi *bHelper_1 = new AliBalancePsi;
1307     bHelper_1->GetFWHM(kProjectInEta,gHistBalanceFunctionSubtracted,fwhm_spline,fwhmError);
1308     Printf("FWHM: %lf - Error: %lf",fwhm_spline, fwhmError);
1309
1310     //store and in txt files FWHM
1311     TString sigmaFileName = filename;
1312     if(kProjectInEta) {
1313       sigmaFileName = "deltaEtaProjection_FWHM.txt";
1314     }
1315     else{
1316       sigmaFileName = "deltaPhiProjection_FWHM.txt";
1317     }
1318     ofstream fileSigma(sigmaFileName.Data(),ios::app);
1319     fileSigma << " " << fwhm_spline << " " <<fwhmError<<endl;
1320     fileSigma.close();
1321
1322     gHistBalanceFunctionSubtracted->SetLineColor(2);
1323     gHistBalanceFunctionSubtracted->SetMarkerColor(2);
1324     gHistBalanceFunctionSubtracted->Draw("SAME");
1325   }
1326   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++//
1327
1328   TLatex *latex = new TLatex();
1329   latex->SetNDC();
1330   latex->SetTextSize(0.035);
1331   latex->SetTextColor(1);
1332   latex->DrawLatex(0.64,0.85,meanLatex.Data());
1333   latex->DrawLatex(0.64,0.81,rmsLatex.Data());
1334   latex->DrawLatex(0.64,0.77,skewnessLatex.Data());
1335   latex->DrawLatex(0.64,0.73,kurtosisLatex.Data());
1336
1337   TString pngName = filename;
1338   if(kProjectInEta) pngName.ReplaceAll(".root","_DeltaEtaProjection.png");
1339   else              pngName.ReplaceAll(".root","_DeltaPhiProjection.png");
1340
1341   c2->SaveAs(pngName.Data());
1342
1343   TString outFileName = filename;
1344   if(kProjectInEta) outFileName.ReplaceAll(".root","_DeltaEtaProjection.root");
1345   else              outFileName.ReplaceAll(".root","_DeltaPhiProjection.root");
1346   TFile *fProjection = TFile::Open(outFileName.Data(),"recreate");
1347   gHistBalanceFunctionSubtracted->Write();
1348   gHistBalanceFunctionMixed->Write();
1349   fProjection->Close();
1350 }
1351
1352 //____________________________________________________________//
1353 void drawBFPsi2DFromCorrelationFunctions(const char* lhcPeriod = "LHC10h",
1354                                          Int_t gTrainNumber = 64,
1355                                          const char* gCentralityEstimator = "V0M",
1356                                          Int_t gBit = 128,
1357                                          const char* gEventPlaneEstimator = "VZERO",
1358                                          Int_t gCentrality = 1,
1359                                          Double_t psiMin = -0.5, Double_t psiMax = 3.5,
1360                                          Double_t vertexZMin = -10.,
1361                                          Double_t vertexZMax = 10.,
1362                                          Double_t ptTriggerMin = -1.,
1363                                          Double_t ptTriggerMax = -1.,
1364                                          Double_t ptAssociatedMin = -1.,
1365                                          Double_t ptAssociatedMax = -1.,
1366                                          TString eventClass = "Multiplicity",
1367                                          Bool_t bRootMoments = kFALSE,
1368                                          Bool_t bFullPhiForEtaProjection = kFALSE,
1369                                          Bool_t bReduceRangeForMoments = kFALSE,
1370                                          Bool_t bFWHM = kFALSE) {
1371   //Macro that draws the BF distributions for each centrality bin
1372   //for reaction plane dependent analysis
1373   //Author: Panos.Christakoglou@nikhef.nl
1374   TGaxis::SetMaxDigits(3);
1375   gStyle->SetPalette(55,0);
1376
1377   //Get the input file
1378   TString filename = lhcPeriod;
1379   if(lhcPeriod != ""){
1380     //filename += "/Train"; filename += gTrainNumber;
1381     filename +="/PttFrom";
1382     filename += Form("%.1f",ptTriggerMin); filename += "To";
1383     filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1384     filename += Form("%.1f",ptAssociatedMin); filename += "To";
1385     filename += Form("%.1f",ptAssociatedMax);
1386     filename += "/correlationFunction.";
1387   }
1388   else{
1389     filename += "correlationFunction.";
1390   }
1391   if(eventClass == "Centrality"){
1392     filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1393     filename += ".PsiAll.PttFrom";
1394   }
1395   else if(eventClass == "Multiplicity"){
1396     filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1397     filename += ".PsiAll.PttFrom";
1398   }
1399   else{ // "EventPlane" (default)
1400     filename += "Centrality";
1401     filename += gCentrality; filename += ".Psi";
1402     if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
1403     else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
1404     else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
1405     else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
1406     else filename += "All.PttFrom";
1407   }
1408   filename += Form("%.1f",ptTriggerMin); filename += "To";
1409   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1410   filename += Form("%.1f",ptAssociatedMin); filename += "To";
1411   filename += Form("%.1f",ptAssociatedMax);
1412   filename += "_";
1413   filename += Form("%.1f",psiMin);
1414   filename += "-";
1415   filename += Form("%.1f",psiMax);
1416   filename += ".root";
1417
1418   //Open the file
1419   TFile *f = TFile::Open(filename.Data());
1420   if((!f)||(!f->IsOpen())) {
1421     Printf("The file %s is not found. Aborting...",filename);
1422     return listBF;
1423   }
1424   //f->ls();
1425
1426   TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
1427   if(!gHistPN) return;
1428   TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
1429   if(!gHistNP) return;
1430   TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
1431   if(!gHistPP) return;
1432   TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
1433   if(!gHistNN) return;
1434
1435   gHistPN->Sumw2();
1436   gHistPP->Sumw2();
1437   gHistPN->Add(gHistPP,-1);
1438   gHistNP->Sumw2();
1439   gHistNN->Sumw2();
1440   gHistNP->Add(gHistNN,-1);
1441   gHistPN->Add(gHistNP);
1442   //gHistPN->Scale(0.5);//not needed anymore, since pT,assoc < pT,trig in all pT bins
1443   TH2D *gHistBalanceFunction2D = dynamic_cast<TH2D *>(gHistPN->Clone());
1444   gHistBalanceFunction2D->SetStats(kFALSE);
1445   gHistBalanceFunction2D->GetXaxis()->SetTitle("#Delta#eta");
1446   gHistBalanceFunction2D->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1447   gHistBalanceFunction2D->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
1448
1449   //Draw the results
1450   TCanvas *c0 = new TCanvas("c0","Balance function 2D",0,0,600,500);
1451   c0->SetFillColor(10); c0->SetHighLightColor(10);
1452   c0->SetLeftMargin(0.17); c0->SetTopMargin(0.05);
1453   gHistBalanceFunction2D->SetTitle("");
1454   gHistBalanceFunction2D->GetZaxis()->SetTitleOffset(1.4);
1455   gHistBalanceFunction2D->GetZaxis()->SetNdivisions(10);
1456   gHistBalanceFunction2D->GetYaxis()->SetTitleOffset(1.4);
1457   gHistBalanceFunction2D->GetYaxis()->SetNdivisions(10);
1458   gHistBalanceFunction2D->GetXaxis()->SetRangeUser(-1.4,1.4);
1459   gHistBalanceFunction2D->GetXaxis()->SetNdivisions(10);
1460   gHistBalanceFunction2D->DrawCopy("lego2");
1461   gPad->SetTheta(30); // default is 30
1462   gPad->SetPhi(-60); // default is 30
1463   gPad->Update();
1464
1465   TString multLatex = Form("Multiplicity: %.1f - %.1f",psiMin,psiMax);
1466
1467   TString pttLatex = Form("%.1f",ptTriggerMin);
1468   pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
1469   pttLatex += " GeV/c";
1470
1471   TString ptaLatex = Form("%.1f",ptAssociatedMin);
1472   ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1473   ptaLatex += " GeV/c";
1474
1475   TLatex *latexInfo1 = new TLatex();
1476   latexInfo1->SetNDC();
1477   latexInfo1->SetTextSize(0.045);
1478   latexInfo1->SetTextColor(1);
1479   latexInfo1->DrawLatex(0.54,0.88,multLatex.Data());
1480   latexInfo1->DrawLatex(0.54,0.82,pttLatex.Data());
1481   latexInfo1->DrawLatex(0.54,0.76,ptaLatex.Data());
1482
1483   TString pngName = "BalanceFunction2D.";
1484   pngName += Form("%s: %.1f - %.1f",eventClass.Data(),psiMin,psiMax);
1485   pngName += ".PttFrom";
1486   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1487   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1488   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1489   pngName += Form("%.1f",ptAssociatedMax);
1490   pngName += ".png";
1491
1492   c0->SaveAs(pngName.Data());
1493
1494  // do the full range for the projection on eta (for cross checking with published results)
1495  if(bFullPhiForEtaProjection){
1496
1497    drawProjections(gHistBalanceFunction2D,
1498                    kTRUE,
1499                    1,72,
1500                    gCentrality,
1501                    psiMin,psiMax,
1502                    ptTriggerMin,ptTriggerMax,
1503                    ptAssociatedMin,ptAssociatedMax,
1504                    kTRUE,
1505                    eventClass.Data(),
1506                    bRootMoments,
1507                    kFALSE,
1508                    bReduceRangeForMoments,
1509                    bFWHM);
1510  }
1511  else{
1512    drawProjections(gHistBalanceFunction2D,
1513                    kTRUE,
1514                    1,36,
1515                    gCentrality,
1516                    psiMin,psiMax,
1517                    ptTriggerMin,ptTriggerMax,
1518                    ptAssociatedMin,ptAssociatedMax,
1519                    kTRUE,
1520                    eventClass.Data(),
1521                    bRootMoments,
1522                    kFALSE,
1523                    bReduceRangeForMoments,
1524                    bFWHM);
1525  }
1526
1527   drawProjections(gHistBalanceFunction2D,
1528                   kFALSE,
1529                   1,80,
1530                   gCentrality,
1531                   psiMin,psiMax,
1532                   ptTriggerMin,ptTriggerMax,
1533                   ptAssociatedMin,ptAssociatedMax,
1534                   kTRUE,
1535                   eventClass.Data(),
1536                   bRootMoments,
1537                   kFALSE,
1538                   bReduceRangeForMoments,
1539                   bFWHM);
1540
1541   TString outFileName = filename;
1542   outFileName.ReplaceAll("correlationFunction","balanceFunction2D");
1543   gHistBalanceFunction2D->SetName("gHistBalanceFunctionSubtracted");
1544   TFile *fOut = TFile::Open(outFileName.Data(),"recreate");
1545   gHistBalanceFunction2D->Write();
1546   fOut->Close();
1547
1548 }
1549
1550
1551 //____________________________________________________________//
1552 void mergeBFPsi2D(TString momDirectory = "./",
1553                   TString directory1 = "",
1554                   TString directory2 = "",
1555                   TString directory3 = "",
1556                   TString directory4 = "",
1557                   Int_t gCentrality = 1,
1558                   Double_t psiMin = -0.5, Double_t psiMax = 3.5,
1559                   TString  eventClass = "Centrality",
1560                   Double_t ptTriggerMin = -1.,
1561                   Double_t ptTriggerMax = -1.,
1562                   Double_t ptAssociatedMin = -1.,
1563                   Double_t ptAssociatedMax = -1.,
1564                   Bool_t bReduceRangeForMoments = kFALSE
1565 ) {
1566   //Macro that draws the BF distributions for each centrality bin
1567   //for reaction plane dependent analysis
1568   //Author: Panos.Christakoglou@nikhef.nl
1569   TGaxis::SetMaxDigits(3);
1570   gStyle->SetPalette(55,0);
1571
1572   cout<<"REDUCE"<<bReduceRangeForMoments<<"  "<<endl;
1573
1574   const Int_t nMaxDirectories = 4; // maximum number of directories to merge (set to 4 for now)
1575   TString sDirectory[nMaxDirectories];
1576   Int_t nDirectories = nMaxDirectories;
1577
1578   TString sFileName[nMaxDirectories];
1579   TFile *fBF[nMaxDirectories];
1580   TH2D  *hBF[nMaxDirectories];
1581   Double_t entries[nMaxDirectories];
1582
1583   TFile *fOut;
1584   TH2D  *hBFOut;
1585   Double_t entriesOut = 0.;
1586
1587   // find out how many directories we have to merge
1588  if(directory1==""){
1589     nDirectories=0;
1590   }
1591   else if(directory2==""){
1592     nDirectories=1;
1593     sDirectory[0]=directory1;
1594   }
1595   else if(directory3==""){
1596     nDirectories=2;
1597     sDirectory[0]=directory1;
1598     sDirectory[1]=directory2;
1599   }
1600   else if(directory4==""){
1601     nDirectories=3;
1602     sDirectory[0]=directory1;
1603     sDirectory[1]=directory2;
1604     sDirectory[2]=directory3;
1605   }
1606   else{
1607     nDirectories=nMaxDirectories;
1608     sDirectory[0]=directory1;
1609     sDirectory[1]=directory2;
1610     sDirectory[2]=directory3;
1611     sDirectory[3]=directory4;
1612   }
1613
1614  for(Int_t iDir = 0; iDir < nDirectories; iDir++){
1615
1616    //Get the input file
1617    sFileName[iDir] = sDirectory[iDir];
1618    sFileName[iDir] += "/Centrality"; sFileName[iDir] += gCentrality;
1619    sFileName[iDir] += "/";
1620 sFileName[iDir] += momDirectory;
1621    sFileName[iDir] += "/balanceFunction2D.";
1622
1623   if(eventClass == "Centrality"){
1624     sFileName[iDir] += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1625     sFileName[iDir] += ".PsiAll.PttFrom";
1626   }
1627   else if(eventClass == "Multiplicity"){
1628     sFileName[iDir] += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1629     sFileName[iDir] += ".PsiAll.PttFrom";
1630   }
1631   else{ // "EventPlane" (default)
1632     sFileName[iDir] += "Centrality";
1633     sFileName[iDir] += gCentrality; sFileName[iDir] += ".Psi";
1634     if((psiMin == -0.5)&&(psiMax == 0.5)) sFileName[iDir] += "InPlane.Ptt";
1635     else if((psiMin == 0.5)&&(psiMax == 1.5)) sFileName[iDir] += "Intermediate.Ptt";
1636     else if((psiMin == 1.5)&&(psiMax == 2.5)) sFileName[iDir] += "OutOfPlane.Ptt";
1637     else if((psiMin == 2.5)&&(psiMax == 3.5)) sFileName[iDir] += "Rest.PttFrom";
1638     else sFileName[iDir] += "All.PttFrom";
1639   }
1640   sFileName[iDir] += Form("%.1f",ptTriggerMin); sFileName[iDir] += "To";
1641   sFileName[iDir] += Form("%.1f",ptTriggerMax); sFileName[iDir] += "PtaFrom";
1642   sFileName[iDir] += Form("%.1f",ptAssociatedMin); sFileName[iDir] += "To";
1643   sFileName[iDir] += Form("%.1f",ptAssociatedMax);
1644   sFileName[iDir] += "_";
1645   sFileName[iDir] += Form("%.1f",psiMin);
1646   sFileName[iDir] += "-";
1647   sFileName[iDir] += Form("%.1f",psiMax);
1648   sFileName[iDir] += ".root";
1649
1650   //Open the file
1651   fBF[iDir] = TFile::Open(sFileName[iDir].Data());
1652   if((!fBF[iDir])||(!fBF[iDir]->IsOpen())) {
1653     Printf("The file %s is not found. Not used...",sFileName[iDir]);
1654     continue;
1655   }
1656   //fBF[iDir]->ls();
1657
1658   hBF[iDir]     = (TH2D*)fBF[iDir]->Get("gHistBalanceFunctionSubtracted");
1659   if(!hBF[iDir]) continue;
1660   entries[iDir] = (Double_t) hBF[iDir]->GetEntries();
1661   entriesOut += entries[iDir];
1662   cout<<" BF histogram "<<iDir<<" has "<<entries[iDir] <<" entries."<<endl;
1663
1664   // scaling and adding (for average)
1665   hBF[iDir]->Scale(entries[iDir]);
1666   if(iDir==0) hBFOut = (TH2D*)hBF[iDir]->Clone("gHistBalanceFunctionSubtractedOut");
1667   else hBFOut->Add(hBF[iDir]);
1668
1669  }
1670
1671  // scaling with all entries
1672  hBFOut->Scale(1./entriesOut);
1673
1674  drawProjections(hBFOut,
1675                  kTRUE,
1676                  1,36,
1677                  gCentrality,
1678                  psiMin,psiMax,
1679                  ptTriggerMin,ptTriggerMax,
1680                  ptAssociatedMin,ptAssociatedMax,
1681                  kTRUE,
1682                  eventClass,
1683                  kTRUE,
1684                  kFALSE,
1685                  bReduceRangeForMoments,
1686                  bFWHM);
1687
1688   drawProjections(hBFOut,
1689                   kFALSE,
1690                   1,80,
1691                   gCentrality,
1692                   psiMin,psiMax,
1693                   ptTriggerMin,ptTriggerMax,
1694                   ptAssociatedMin,ptAssociatedMax,
1695                   kTRUE,
1696                   eventClass.Data(),
1697                   kTRUE,
1698                   kFALSE,
1699                   bReduceRangeForMoments,
1700                   bFWHM);
1701
1702  // output
1703  TString outFileName = "balanceFunction2D.";
1704
1705  if(eventClass == "Centrality"){
1706    outFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1707    outFileName += ".PsiAll.PttFrom";
1708  }
1709  else if(eventClass == "Multiplicity"){
1710    outFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1711    outFileName += ".PsiAll.PttFrom";
1712   }
1713  else{ // "EventPlane" (default)
1714    outFileName += "Centrality";
1715    outFileName += gCentrality; outFileName += ".Psi";
1716    if((psiMin == -0.5)&&(psiMax == 0.5)) outFileName += "InPlane.Ptt";
1717    else if((psiMin == 0.5)&&(psiMax == 1.5)) outFileName += "Intermediate.Ptt";
1718    else if((psiMin == 1.5)&&(psiMax == 2.5)) outFileName += "OutOfPlane.Ptt";
1719    else if((psiMin == 2.5)&&(psiMax == 3.5)) outFileName += "Rest.PttFrom";
1720    else outFileName += "All.PttFrom";
1721  }
1722  outFileName += Form("%.1f",ptTriggerMin); outFileName += "To";
1723  outFileName += Form("%.1f",ptTriggerMax); outFileName += "PtaFrom";
1724  outFileName += Form("%.1f",ptAssociatedMin); outFileName += "To";
1725  outFileName += Form("%.1f",ptAssociatedMax);
1726  outFileName += "_";
1727  outFileName += Form("%.1f",psiMin);
1728  outFileName += "-";
1729  outFileName += Form("%.1f",psiMax);
1730  outFileName += ".root";
1731
1732  hBFOut->SetName("gHistBalanceFunctionSubtracted");
1733  fBFOut = TFile::Open(outFileName.Data(),"recreate");
1734  hBFOut->Write();
1735  fBFOut->Close();
1736
1737 }
1738
1739 //____________________________________________________________________//
1740 void GetWeightedMean1D(TH1D *gHistBalance, Bool_t kProjectInEta = kTRUE, Int_t fStartBin = 1, Int_t fStopBin = -1, Double_t &WM, Double_t &WME,Double_t rangeReduced = -1.) {
1741
1742   //Prints the calculated width of the BF and its error
1743   Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
1744   Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
1745   Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
1746   Double_t deltaBalP2 = 0.0, integral = 0.0;
1747   Double_t deltaErrorNew = 0.0;
1748
1749   //Retrieve this variables from Histogram
1750   Int_t fNumberOfBins = gHistBalance->GetNbinsX();
1751   if(fStopBin > -1)   fNumberOfBins = fStopBin;
1752   Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!
1753   Double_t currentBinCenter = 0.;
1754
1755   for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
1756
1757     // in order to recover the |abs| in the 1D analysis
1758     currentBinCenter = gHistBalance->GetBinCenter(i);
1759     if(kProjectInEta){
1760       if(currentBinCenter < 0) currentBinCenter = -currentBinCenter;
1761     }
1762     else{
1763       if(currentBinCenter < 0) currentBinCenter = -currentBinCenter;
1764       if(currentBinCenter > TMath::Pi()) currentBinCenter = 2 * TMath::Pi() - currentBinCenter;
1765     }
1766
1767     if(rangeReduced > 0 && currentBinCenter > rangeReduced )
1768       continue;
1769
1770     gSumXi += currentBinCenter;
1771     gSumBi += gHistBalance->GetBinContent(i);
1772     gSumBiXi += gHistBalance->GetBinContent(i)*currentBinCenter;
1773     gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(currentBinCenter,2);
1774     gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(currentBinCenter,2);
1775     gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
1776     gSumXi2DeltaBi2 += TMath::Power(currentBinCenter,2) * TMath::Power(gHistBalance->GetBinError(i),2);
1777
1778     deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
1779     integral += fP2Step*gHistBalance->GetBinContent(i);
1780   }
1781   for(Int_t i = fStartBin; i < fNumberOfBins; i++)
1782     deltaErrorNew += gHistBalance->GetBinError(i)*(currentBinCenter*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
1783
1784   Double_t integralError = TMath::Sqrt(deltaBalP2);
1785
1786   Double_t delta = gSumBiXi / gSumBi;
1787   Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
1788
1789   WM  = delta;
1790   WME = deltaError;
1791 }
1792