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       // safety check:
1105       // default (for 2<pT,assoc<3<pT,trig<4)
1106       Double_t rangeReduced = 1.2;
1107       //for 3<pT,assoc<8<pT,trig<15
1108       if(ptTriggerMax>10){
1109         rangeReduced = 0.35;
1110       }
1111
1112       // new define range by fit!
1113       gHistBalanceFunctionSubtracted->Fit("gaus","","");
1114       Double_t sigmaGaus = gHistBalanceFunctionSubtracted->GetFunction("gaus")->GetParameter(2);
1115
1116       // if safety check OK, set rangeReduced to 3 sigma of gauss fit
1117       if(3*sigmaGaus > rangeReduced){
1118         cout<<"RANGE REDUCE ERROR: "<<3*sigmaGaus<<" > "<<rangeReduced<<endl;
1119       }
1120       else{
1121         rangeReduced = 3*sigmaGaus;
1122       }
1123
1124       cout<<"Use reduced range = "<<rangeReduced<<endl;
1125       gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-rangeReduced,rangeReduced);
1126     }
1127
1128     meanLatex = "#mu = ";
1129     meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMean());
1130     meanLatex += " #pm ";
1131     meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMeanError());
1132
1133     rmsLatex = "#sigma = ";
1134     rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMS());
1135     rmsLatex += " #pm ";
1136     rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMSError());
1137
1138     skewnessLatex = "S = ";
1139     skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(1));
1140     skewnessLatex += " #pm ";
1141     skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(11));
1142
1143     kurtosisLatex = "K = ";
1144     kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(1));
1145     kurtosisLatex += " #pm ";
1146     kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(11));
1147
1148     Printf("ROOT Mean: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetMean(),gHistBalanceFunctionSubtracted->GetMeanError());
1149     Printf("ROOT RMS: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetRMS(),gHistBalanceFunctionSubtracted->GetRMSError());
1150     Printf("ROOT Skeweness: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetSkewness(1),gHistBalanceFunctionSubtracted->GetSkewness(11));
1151     Printf("ROOT Kurtosis: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetKurtosis(1),gHistBalanceFunctionSubtracted->GetKurtosis(11));
1152
1153
1154     // store in txt files
1155     TString meanFileName = filename;
1156     if(kProjectInEta)
1157       meanFileName= "deltaEtaProjection_Mean.txt";
1158     //meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
1159     else
1160       meanFileName = "deltaPhiProjection_Mean.txt";
1161       //meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
1162     ofstream fileMean(meanFileName.Data(),ios::app);
1163     fileMean << " " << gHistBalanceFunctionSubtracted->GetMean() << " " <<gHistBalanceFunctionSubtracted->GetMeanError()<<endl;
1164     fileMean.close();
1165
1166     TString rmsFileName = filename;
1167     if(kProjectInEta)
1168       rmsFileName = "deltaEtaProjection_Rms.txt";
1169       //rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
1170     else
1171       rmsFileName = "deltaPhiProjection_Rms.txt";
1172       //rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
1173     ofstream fileRms(rmsFileName.Data(),ios::app);
1174     fileRms << " " << gHistBalanceFunctionSubtracted->GetRMS() << " " <<gHistBalanceFunctionSubtracted->GetRMSError()<<endl;
1175     fileRms.close();
1176
1177     TString skewnessFileName = filename;
1178     if(kProjectInEta)
1179       skewnessFileName = "deltaEtaProjection_Skewness.txt";
1180       //skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
1181     else
1182       skewnessFileName = "deltaPhiProjection_Skewness.txt";
1183     //skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
1184     ofstream fileSkewness(skewnessFileName.Data(),ios::app);
1185     fileSkewness << " " << gHistBalanceFunctionSubtracted->GetSkewness(1) << " " <<gHistBalanceFunctionSubtracted->GetSkewness(11)<<endl;
1186     fileSkewness.close();
1187
1188     TString kurtosisFileName = filename;
1189     if(kProjectInEta)
1190       kurtosisFileName = "deltaEtaProjection_Kurtosis.txt";
1191       //kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
1192     else
1193       kurtosisFileName = "deltaPhiProjection_Kurtosis.txt";
1194       //kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
1195     ofstream fileKurtosis(kurtosisFileName.Data(),ios::app);
1196     fileKurtosis << " " << gHistBalanceFunctionSubtracted->GetKurtosis(1) << " " <<gHistBalanceFunctionSubtracted->GetKurtosis(11)<<endl;
1197     fileKurtosis.close();
1198
1199     if(!kProjectInEta) gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-TMath::Pi()/2,3.*TMath::Pi()/2);
1200     else               gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-1.6,1.6);
1201   }
1202   // calculate the moments by hand
1203   else{
1204
1205     Double_t meanAnalytical, meanAnalyticalError;
1206     Double_t sigmaAnalytical, sigmaAnalyticalError;
1207     Double_t skewnessAnalytical, skewnessAnalyticalError;
1208     Double_t kurtosisAnalytical, kurtosisAnalyticalError;
1209
1210     Int_t gDeltaEtaPhi = 2;
1211     if(kProjectInEta) gDeltaEtaPhi = 1;
1212
1213     AliBalancePsi *bHelper = new AliBalancePsi;
1214     bHelper->GetMomentsAnalytical(gDeltaEtaPhi,gHistBalanceFunctionSubtracted,kUseZYAM,meanAnalytical,meanAnalyticalError,sigmaAnalytical,sigmaAnalyticalError,skewnessAnalytical,skewnessAnalyticalError,kurtosisAnalytical,kurtosisAnalyticalError);
1215
1216     meanLatex = "#mu = ";
1217     meanLatex += Form("%.3f",meanAnalytical);
1218     meanLatex += " #pm ";
1219     meanLatex += Form("%.3f",meanAnalyticalError);
1220
1221     rmsLatex = "#sigma = ";
1222     rmsLatex += Form("%.3f",sigmaAnalytical);
1223     rmsLatex += " #pm ";
1224     rmsLatex += Form("%.3f",sigmaAnalyticalError);
1225
1226     skewnessLatex = "S = ";
1227     skewnessLatex += Form("%.3f",skewnessAnalytical);
1228     skewnessLatex += " #pm ";
1229     skewnessLatex += Form("%.3f",skewnessAnalyticalError);
1230
1231     kurtosisLatex = "K = ";
1232     kurtosisLatex += Form("%.3f",kurtosisAnalytical);
1233     kurtosisLatex += " #pm ";
1234     kurtosisLatex += Form("%.3f",kurtosisAnalyticalError);
1235     Printf("Mean: %lf - Error: %lf",meanAnalytical, meanAnalyticalError);
1236     Printf("Sigma: %lf - Error: %lf",sigmaAnalytical, sigmaAnalyticalError);
1237     Printf("Skeweness: %lf - Error: %lf",skewnessAnalytical, skewnessAnalyticalError);
1238     Printf("Kurtosis: %lf - Error: %lf",kurtosisAnalytical, kurtosisAnalyticalError);
1239
1240     // store in txt files
1241     TString meanFileName = filename;
1242     if(kProjectInEta)
1243       meanFileName = "deltaEtaProjection_Mean.txt";
1244       //meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
1245     else
1246       meanFileName = "deltaPhiProjection_Mean.txt";
1247       //meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
1248     ofstream fileMean(meanFileName.Data(),ios::app);
1249     fileMean << " " << meanAnalytical << " " <<meanAnalyticalError<<endl;
1250     fileMean.close();
1251
1252     TString rmsFileName = filename;
1253     if(kProjectInEta)
1254       rmsFileName = "deltaEtaProjection_Rms.txt";
1255       //rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
1256     else
1257       rmsFileName = "deltaPhiProjection_Rms.txt";
1258 //rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
1259     ofstream fileRms(rmsFileName.Data(),ios::app);
1260     fileRms << " " << sigmaAnalytical << " " <<sigmaAnalyticalError<<endl;
1261     fileRms.close();
1262
1263     TString skewnessFileName = filename;
1264     if(kProjectInEta)
1265       skewnessFileName = "deltaEtaProjection_Skewness.txt";
1266       //skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
1267     else
1268       skewnessFileName = "deltaPhiProjection_Skewness.txt";
1269       //skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
1270     ofstream fileSkewness(skewnessFileName.Data(),ios::app);
1271     fileSkewness << " " << skewnessAnalytical << " " <<skewnessAnalyticalError<<endl;
1272     fileSkewness.close();
1273
1274     TString kurtosisFileName = filename;
1275     if(kProjectInEta)
1276       kurtosisFileName = "deltaEtaProjection_Kurtosis.txt";
1277       //kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
1278     else
1279       kurtosisFileName = "deltaPhiProjection_Kurtosis.txt";
1280       //kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
1281     ofstream fileKurtosis(kurtosisFileName.Data(),ios::app);
1282     fileKurtosis << " " << kurtosisAnalytical << " " <<kurtosisAnalyticalError<<endl;
1283     fileKurtosis.close();
1284   }
1285
1286   // Weighted mean as calculated for 1D analysis
1287   Double_t weightedMean, weightedMeanError;
1288   if(bReduceRangeForMoments){
1289     GetWeightedMean1D(gHistBalanceFunctionSubtracted,kProjectInEta,1,-1,weightedMean,weightedMeanError,rangeReduced);
1290   }
1291   else{
1292     GetWeightedMean1D(gHistBalanceFunctionSubtracted,kProjectInEta,1,-1,weightedMean,weightedMeanError,-1.);
1293   }
1294   Printf("Weighted Mean: %lf - Error: %lf",weightedMean, weightedMeanError);
1295
1296   // store in txt files
1297   TString weightedMeanFileName = filename;
1298   if(kProjectInEta)
1299     weightedMeanFileName = "deltaEtaProjection_WeightedMean.txt";
1300     //weightedMeanFileName.ReplaceAll(".root","_DeltaEtaProjection_WeightedMean.txt");
1301   else
1302     weightedMeanFileName = "deltaPhiProjection_WeightedMean.txt";
1303     //weightedMeanFileName.ReplaceAll(".root","_DeltaPhiProjection_WeightedMean.txt");
1304   ofstream fileWeightedMean(weightedMeanFileName.Data(),ios::app);
1305   fileWeightedMean << " " << weightedMean << " " <<weightedMeanError<<endl;
1306   fileWeightedMean.close();
1307
1308   TCanvas *c2 = new TCanvas("c2","",600,0,600,500);
1309   c2->SetFillColor(10);
1310   c2->SetHighLightColor(10);
1311   c2->SetLeftMargin(0.15);
1312   gHistBalanceFunctionSubtracted->DrawCopy("E");
1313
1314   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++//
1315   Double_t fwhm_spline = 0.;
1316   Double_t fwhmError = 0.;
1317   if (bFWHM){
1318     AliBalancePsi *bHelper_1 = new AliBalancePsi;
1319     bHelper_1->GetFWHM(kProjectInEta,gHistBalanceFunctionSubtracted,fwhm_spline,fwhmError);
1320     Printf("FWHM: %lf - Error: %lf",fwhm_spline, fwhmError);
1321
1322     //store and in txt files FWHM
1323     TString sigmaFileName = filename;
1324     if(kProjectInEta) {
1325       sigmaFileName = "deltaEtaProjection_FWHM.txt";
1326     }
1327     else{
1328       sigmaFileName = "deltaPhiProjection_FWHM.txt";
1329     }
1330     ofstream fileSigma(sigmaFileName.Data(),ios::app);
1331     fileSigma << " " << fwhm_spline << " " <<fwhmError<<endl;
1332     fileSigma.close();
1333
1334     gHistBalanceFunctionSubtracted->SetLineColor(2);
1335     gHistBalanceFunctionSubtracted->SetMarkerColor(2);
1336     gHistBalanceFunctionSubtracted->Draw("SAME");
1337   }
1338   //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++//
1339
1340   TLatex *latex = new TLatex();
1341   latex->SetNDC();
1342   latex->SetTextSize(0.035);
1343   latex->SetTextColor(1);
1344   latex->DrawLatex(0.64,0.85,meanLatex.Data());
1345   latex->DrawLatex(0.64,0.81,rmsLatex.Data());
1346   latex->DrawLatex(0.64,0.77,skewnessLatex.Data());
1347   latex->DrawLatex(0.64,0.73,kurtosisLatex.Data());
1348
1349   TString pngName = filename;
1350   if(kProjectInEta) pngName.ReplaceAll(".root","_DeltaEtaProjection.png");
1351   else              pngName.ReplaceAll(".root","_DeltaPhiProjection.png");
1352
1353   c2->SaveAs(pngName.Data());
1354
1355   TString outFileName = filename;
1356   if(kProjectInEta) outFileName.ReplaceAll(".root","_DeltaEtaProjection.root");
1357   else              outFileName.ReplaceAll(".root","_DeltaPhiProjection.root");
1358   TFile *fProjection = TFile::Open(outFileName.Data(),"recreate");
1359   gHistBalanceFunctionSubtracted->Write();
1360   gHistBalanceFunctionMixed->Write();
1361   fProjection->Close();
1362 }
1363
1364 //____________________________________________________________//
1365 void drawBFPsi2DFromCorrelationFunctions(const char* lhcPeriod = "LHC10h",
1366                                          Int_t gTrainNumber = 64,
1367                                          const char* gCentralityEstimator = "V0M",
1368                                          Int_t gBit = 128,
1369                                          const char* gEventPlaneEstimator = "VZERO",
1370                                          Int_t gCentrality = 1,
1371                                          Double_t psiMin = -0.5, Double_t psiMax = 3.5,
1372                                          Double_t vertexZMin = -10.,
1373                                          Double_t vertexZMax = 10.,
1374                                          Double_t ptTriggerMin = -1.,
1375                                          Double_t ptTriggerMax = -1.,
1376                                          Double_t ptAssociatedMin = -1.,
1377                                          Double_t ptAssociatedMax = -1.,
1378                                          TString eventClass = "Multiplicity",
1379                                          Bool_t bRootMoments = kFALSE,
1380                                          Bool_t bFullPhiForEtaProjection = kFALSE,
1381                                          Bool_t bReduceRangeForMoments = kFALSE,
1382                                          Bool_t bFWHM = kFALSE) {
1383   //Macro that draws the BF distributions for each centrality bin
1384   //for reaction plane dependent analysis
1385   //Author: Panos.Christakoglou@nikhef.nl
1386   TGaxis::SetMaxDigits(3);
1387   gStyle->SetPalette(55,0);
1388
1389   //Get the input file
1390   TString filename = lhcPeriod;
1391   if(lhcPeriod != ""){
1392     //filename += "/Train"; filename += gTrainNumber;
1393     filename +="/PttFrom";
1394     filename += Form("%.1f",ptTriggerMin); filename += "To";
1395     filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1396     filename += Form("%.1f",ptAssociatedMin); filename += "To";
1397     filename += Form("%.1f",ptAssociatedMax);
1398     filename += "/correlationFunction.";
1399   }
1400   else{
1401     filename += "correlationFunction.";
1402   }
1403   if(eventClass == "Centrality"){
1404     filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1405     filename += ".PsiAll.PttFrom";
1406   }
1407   else if(eventClass == "Multiplicity"){
1408     filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1409     filename += ".PsiAll.PttFrom";
1410   }
1411   else{ // "EventPlane" (default)
1412     filename += "Centrality";
1413     filename += gCentrality; filename += ".Psi";
1414     if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
1415     else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
1416     else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
1417     else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
1418     else filename += "All.PttFrom";
1419   }
1420   filename += Form("%.1f",ptTriggerMin); filename += "To";
1421   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1422   filename += Form("%.1f",ptAssociatedMin); filename += "To";
1423   filename += Form("%.1f",ptAssociatedMax);
1424   filename += "_";
1425   filename += Form("%.1f",psiMin);
1426   filename += "-";
1427   filename += Form("%.1f",psiMax);
1428   filename += ".root";
1429
1430   //Open the file
1431   TFile *f = TFile::Open(filename.Data());
1432   if((!f)||(!f->IsOpen())) {
1433     Printf("The file %s is not found. Aborting...",filename);
1434     return listBF;
1435   }
1436   //f->ls();
1437
1438   TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
1439   if(!gHistPN) return;
1440   TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
1441   if(!gHistNP) return;
1442   TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
1443   if(!gHistPP) return;
1444   TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
1445   if(!gHistNN) return;
1446
1447   // in order to get unzoomed (in older versions used smaller user ranger)
1448   Int_t binMinX = gHistPN->GetXaxis()->FindBin(gHistPN->GetXaxis()->GetXmin()+0.00001);
1449   Int_t binMaxX = gHistPN->GetXaxis()->FindBin(gHistPN->GetXaxis()->GetXmax()-0.00001);
1450   Int_t binMinY = gHistPN->GetYaxis()->FindBin(gHistPN->GetYaxis()->GetXmin()+0.00001);
1451   Int_t binMaxY = gHistPN->GetYaxis()->FindBin(gHistPN->GetYaxis()->GetXmax()-0.00001);
1452   gHistPN->GetXaxis()->SetRange(binMinX,binMaxX); gHistPN->GetYaxis()->SetRange(binMinY,binMaxY);
1453   gHistNP->GetXaxis()->SetRange(binMinX,binMaxX); gHistNP->GetYaxis()->SetRange(binMinY,binMaxY);
1454   gHistPP->GetXaxis()->SetRange(binMinX,binMaxX); gHistPP->GetYaxis()->SetRange(binMinY,binMaxY);
1455   gHistNN->GetXaxis()->SetRange(binMinX,binMaxX); gHistNN->GetYaxis()->SetRange(binMinY,binMaxY);
1456
1457   gHistPN->Sumw2();
1458   gHistPP->Sumw2();
1459   gHistPN->Add(gHistPP,-1);
1460   gHistNP->Sumw2();
1461   gHistNN->Sumw2();
1462   gHistNP->Add(gHistNN,-1);
1463   gHistPN->Add(gHistNP);
1464   //gHistPN->Scale(0.5);//not needed anymore, since pT,assoc < pT,trig in all pT bins
1465   TH2D *gHistBalanceFunction2D = dynamic_cast<TH2D *>(gHistPN->Clone());
1466   gHistBalanceFunction2D->SetStats(kFALSE);
1467   gHistBalanceFunction2D->GetXaxis()->SetTitle("#Delta#eta");
1468   gHistBalanceFunction2D->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1469   gHistBalanceFunction2D->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
1470
1471   //Draw the results
1472   TCanvas *c0 = new TCanvas("c0","Balance function 2D",0,0,600,500);
1473   c0->SetFillColor(10); c0->SetHighLightColor(10);
1474   c0->SetLeftMargin(0.17); c0->SetTopMargin(0.05);
1475   gHistBalanceFunction2D->SetTitle("");
1476   gHistBalanceFunction2D->GetZaxis()->SetTitleOffset(1.4);
1477   gHistBalanceFunction2D->GetZaxis()->SetNdivisions(10);
1478   gHistBalanceFunction2D->GetYaxis()->SetTitleOffset(1.4);
1479   gHistBalanceFunction2D->GetYaxis()->SetNdivisions(10);
1480   //gHistBalanceFunction2D->GetXaxis()->SetRangeUser(-1.4,1.4);
1481   gHistBalanceFunction2D->GetXaxis()->SetNdivisions(10);
1482   gHistBalanceFunction2D->DrawCopy("lego2");
1483   gPad->SetTheta(30); // default is 30
1484   gPad->SetPhi(-60); // default is 30
1485   gPad->Update();
1486
1487   TString multLatex = Form("Multiplicity: %.1f - %.1f",psiMin,psiMax);
1488
1489   TString pttLatex = Form("%.1f",ptTriggerMin);
1490   pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
1491   pttLatex += " GeV/c";
1492
1493   TString ptaLatex = Form("%.1f",ptAssociatedMin);
1494   ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1495   ptaLatex += " GeV/c";
1496
1497   TLatex *latexInfo1 = new TLatex();
1498   latexInfo1->SetNDC();
1499   latexInfo1->SetTextSize(0.045);
1500   latexInfo1->SetTextColor(1);
1501   latexInfo1->DrawLatex(0.54,0.88,multLatex.Data());
1502   latexInfo1->DrawLatex(0.54,0.82,pttLatex.Data());
1503   latexInfo1->DrawLatex(0.54,0.76,ptaLatex.Data());
1504
1505   TString pngName = "BalanceFunction2D.";
1506   pngName += Form("%s: %.1f - %.1f",eventClass.Data(),psiMin,psiMax);
1507   pngName += ".PttFrom";
1508   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1509   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1510   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1511   pngName += Form("%.1f",ptAssociatedMax);
1512   pngName += ".png";
1513
1514   c0->SaveAs(pngName.Data());
1515
1516  // do the full range for the projection on eta (for cross checking with published results)
1517  if(bFullPhiForEtaProjection){
1518
1519    drawProjections(gHistBalanceFunction2D,
1520                    kTRUE,
1521                    1,72,
1522                    gCentrality,
1523                    psiMin,psiMax,
1524                    ptTriggerMin,ptTriggerMax,
1525                    ptAssociatedMin,ptAssociatedMax,
1526                    kTRUE,
1527                    eventClass.Data(),
1528                    bRootMoments,
1529                    kFALSE,
1530                    bReduceRangeForMoments,
1531                    bFWHM);
1532  }
1533  else{
1534    drawProjections(gHistBalanceFunction2D,
1535                    kTRUE,
1536                    1,36,
1537                    gCentrality,
1538                    psiMin,psiMax,
1539                    ptTriggerMin,ptTriggerMax,
1540                    ptAssociatedMin,ptAssociatedMax,
1541                    kTRUE,
1542                    eventClass.Data(),
1543                    bRootMoments,
1544                    kFALSE,
1545                    bReduceRangeForMoments,
1546                    bFWHM);
1547  }
1548
1549   drawProjections(gHistBalanceFunction2D,
1550                   kFALSE,
1551                   1,80,
1552                   gCentrality,
1553                   psiMin,psiMax,
1554                   ptTriggerMin,ptTriggerMax,
1555                   ptAssociatedMin,ptAssociatedMax,
1556                   kTRUE,
1557                   eventClass.Data(),
1558                   bRootMoments,
1559                   kFALSE,
1560                   bReduceRangeForMoments,
1561                   bFWHM);
1562
1563   TString outFileName = filename;
1564   outFileName.ReplaceAll("correlationFunction","balanceFunction2D");
1565   gHistBalanceFunction2D->SetName("gHistBalanceFunctionSubtracted");
1566   TFile *fOut = TFile::Open(outFileName.Data(),"recreate");
1567   gHistBalanceFunction2D->Write();
1568   fOut->Close();
1569
1570 }
1571
1572
1573 //____________________________________________________________//
1574 void mergeBFPsi2D(TString momDirectory = "./",
1575                   TString directory1 = "",
1576                   TString directory2 = "",
1577                   TString directory3 = "",
1578                   TString directory4 = "",
1579                   Int_t gCentrality = 1,
1580                   Double_t psiMin = -0.5, Double_t psiMax = 3.5,
1581                   TString  eventClass = "Centrality",
1582                   Double_t ptTriggerMin = -1.,
1583                   Double_t ptTriggerMax = -1.,
1584                   Double_t ptAssociatedMin = -1.,
1585                   Double_t ptAssociatedMax = -1.,
1586                   Bool_t bReduceRangeForMoments = kFALSE,
1587                   Bool_t bFWHM = kFALSE
1588 ) {
1589   //Macro that draws the BF distributions for each centrality bin
1590   //for reaction plane dependent analysis
1591   //Author: Panos.Christakoglou@nikhef.nl
1592   TGaxis::SetMaxDigits(3);
1593   gStyle->SetPalette(55,0);
1594
1595   const Int_t nMaxDirectories = 4; // maximum number of directories to merge (set to 4 for now)
1596   TString sDirectory[nMaxDirectories];
1597   Int_t nDirectories = nMaxDirectories;
1598
1599   TString sFileName[nMaxDirectories];
1600   TFile *fBF[nMaxDirectories];
1601   TH2D  *hBF[nMaxDirectories];
1602   Double_t entries[nMaxDirectories];
1603
1604   TFile *fOut;
1605   TH2D  *hBFOut;
1606   Double_t entriesOut = 0.;
1607
1608   // find out how many directories we have to merge
1609  if(directory1==""){
1610     nDirectories=0;
1611   }
1612   else if(directory2==""){
1613     nDirectories=1;
1614     sDirectory[0]=directory1;
1615   }
1616   else if(directory3==""){
1617     nDirectories=2;
1618     sDirectory[0]=directory1;
1619     sDirectory[1]=directory2;
1620   }
1621   else if(directory4==""){
1622     nDirectories=3;
1623     sDirectory[0]=directory1;
1624     sDirectory[1]=directory2;
1625     sDirectory[2]=directory3;
1626   }
1627   else{
1628     nDirectories=nMaxDirectories;
1629     sDirectory[0]=directory1;
1630     sDirectory[1]=directory2;
1631     sDirectory[2]=directory3;
1632     sDirectory[3]=directory4;
1633   }
1634
1635  for(Int_t iDir = 0; iDir < nDirectories; iDir++){
1636
1637    //Get the input file
1638    sFileName[iDir] = sDirectory[iDir];
1639    sFileName[iDir] += "/Centrality"; sFileName[iDir] += gCentrality;
1640    sFileName[iDir] += "/";
1641 sFileName[iDir] += momDirectory;
1642    sFileName[iDir] += "/balanceFunction2D.";
1643
1644   if(eventClass == "Centrality"){
1645     sFileName[iDir] += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1646     sFileName[iDir] += ".PsiAll.PttFrom";
1647   }
1648   else if(eventClass == "Multiplicity"){
1649     sFileName[iDir] += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1650     sFileName[iDir] += ".PsiAll.PttFrom";
1651   }
1652   else{ // "EventPlane" (default)
1653     sFileName[iDir] += "Centrality";
1654     sFileName[iDir] += gCentrality; sFileName[iDir] += ".Psi";
1655     if((psiMin == -0.5)&&(psiMax == 0.5)) sFileName[iDir] += "InPlane.Ptt";
1656     else if((psiMin == 0.5)&&(psiMax == 1.5)) sFileName[iDir] += "Intermediate.Ptt";
1657     else if((psiMin == 1.5)&&(psiMax == 2.5)) sFileName[iDir] += "OutOfPlane.Ptt";
1658     else if((psiMin == 2.5)&&(psiMax == 3.5)) sFileName[iDir] += "Rest.PttFrom";
1659     else sFileName[iDir] += "All.PttFrom";
1660   }
1661   sFileName[iDir] += Form("%.1f",ptTriggerMin); sFileName[iDir] += "To";
1662   sFileName[iDir] += Form("%.1f",ptTriggerMax); sFileName[iDir] += "PtaFrom";
1663   sFileName[iDir] += Form("%.1f",ptAssociatedMin); sFileName[iDir] += "To";
1664   sFileName[iDir] += Form("%.1f",ptAssociatedMax);
1665   sFileName[iDir] += "_";
1666   sFileName[iDir] += Form("%.1f",psiMin);
1667   sFileName[iDir] += "-";
1668   sFileName[iDir] += Form("%.1f",psiMax);
1669   sFileName[iDir] += ".root";
1670
1671   //Open the file
1672   fBF[iDir] = TFile::Open(sFileName[iDir].Data());
1673   if((!fBF[iDir])||(!fBF[iDir]->IsOpen())) {
1674     Printf("The file %s is not found. Not used...",sFileName[iDir]);
1675     continue;
1676   }
1677   //fBF[iDir]->ls();
1678
1679   hBF[iDir]     = (TH2D*)fBF[iDir]->Get("gHistBalanceFunctionSubtracted");
1680   if(!hBF[iDir]) continue;
1681   entries[iDir] = (Double_t) hBF[iDir]->GetEntries();
1682   entriesOut += entries[iDir];
1683   cout<<" BF histogram "<<iDir<<" has "<<entries[iDir] <<" entries."<<endl;
1684
1685   // scaling and adding (for average)
1686   hBF[iDir]->Scale(entries[iDir]);
1687   if(iDir==0) hBFOut = (TH2D*)hBF[iDir]->Clone("gHistBalanceFunctionSubtractedOut");
1688   else hBFOut->Add(hBF[iDir]);
1689
1690  }
1691
1692  // scaling with all entries
1693  hBFOut->Scale(1./entriesOut);
1694
1695  drawProjections(hBFOut,
1696                  kTRUE,
1697                  1,36,
1698                  gCentrality,
1699                  psiMin,psiMax,
1700                  ptTriggerMin,ptTriggerMax,
1701                  ptAssociatedMin,ptAssociatedMax,
1702                  kTRUE,
1703                  eventClass,
1704                  kTRUE,
1705                  kFALSE,
1706                  bReduceRangeForMoments,
1707                  bFWHM);
1708
1709   drawProjections(hBFOut,
1710                   kFALSE,
1711                   1,80,
1712                   gCentrality,
1713                   psiMin,psiMax,
1714                   ptTriggerMin,ptTriggerMax,
1715                   ptAssociatedMin,ptAssociatedMax,
1716                   kTRUE,
1717                   eventClass.Data(),
1718                   kTRUE,
1719                   kFALSE,
1720                   bReduceRangeForMoments,
1721                   bFWHM);
1722
1723  // output
1724  TString outFileName = "balanceFunction2D.";
1725
1726  if(eventClass == "Centrality"){
1727    outFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1728    outFileName += ".PsiAll.PttFrom";
1729  }
1730  else if(eventClass == "Multiplicity"){
1731    outFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1732    outFileName += ".PsiAll.PttFrom";
1733   }
1734  else{ // "EventPlane" (default)
1735    outFileName += "Centrality";
1736    outFileName += gCentrality; outFileName += ".Psi";
1737    if((psiMin == -0.5)&&(psiMax == 0.5)) outFileName += "InPlane.Ptt";
1738    else if((psiMin == 0.5)&&(psiMax == 1.5)) outFileName += "Intermediate.Ptt";
1739    else if((psiMin == 1.5)&&(psiMax == 2.5)) outFileName += "OutOfPlane.Ptt";
1740    else if((psiMin == 2.5)&&(psiMax == 3.5)) outFileName += "Rest.PttFrom";
1741    else outFileName += "All.PttFrom";
1742  }
1743  outFileName += Form("%.1f",ptTriggerMin); outFileName += "To";
1744  outFileName += Form("%.1f",ptTriggerMax); outFileName += "PtaFrom";
1745  outFileName += Form("%.1f",ptAssociatedMin); outFileName += "To";
1746  outFileName += Form("%.1f",ptAssociatedMax);
1747  outFileName += "_";
1748  outFileName += Form("%.1f",psiMin);
1749  outFileName += "-";
1750  outFileName += Form("%.1f",psiMax);
1751  outFileName += ".root";
1752
1753  hBFOut->SetName("gHistBalanceFunctionSubtracted");
1754  fBFOut = TFile::Open(outFileName.Data(),"recreate");
1755  hBFOut->Write();
1756  fBFOut->Close();
1757
1758 }
1759
1760 //____________________________________________________________________//
1761 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.) {
1762
1763   //Prints the calculated width of the BF and its error
1764   Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
1765   Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
1766   Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
1767   Double_t deltaBalP2 = 0.0, integral = 0.0;
1768   Double_t deltaErrorNew = 0.0;
1769
1770   //Retrieve this variables from Histogram
1771   Int_t fNumberOfBins = gHistBalance->GetNbinsX();
1772   if(fStopBin > -1)   fNumberOfBins = fStopBin;
1773   Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!
1774   Double_t currentBinCenter = 0.;
1775
1776   for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
1777
1778     // in order to recover the |abs| in the 1D analysis
1779     currentBinCenter = gHistBalance->GetBinCenter(i);
1780     if(kProjectInEta){
1781       if(currentBinCenter < 0) currentBinCenter = -currentBinCenter;
1782     }
1783     else{
1784       if(currentBinCenter < 0) currentBinCenter = -currentBinCenter;
1785       if(currentBinCenter > TMath::Pi()) currentBinCenter = 2 * TMath::Pi() - currentBinCenter;
1786     }
1787
1788     if(rangeReduced > 0 && currentBinCenter > rangeReduced )
1789       continue;
1790
1791     gSumXi += currentBinCenter;
1792     gSumBi += gHistBalance->GetBinContent(i);
1793     gSumBiXi += gHistBalance->GetBinContent(i)*currentBinCenter;
1794     gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(currentBinCenter,2);
1795     gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(currentBinCenter,2);
1796     gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
1797     gSumXi2DeltaBi2 += TMath::Power(currentBinCenter,2) * TMath::Power(gHistBalance->GetBinError(i),2);
1798
1799     deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
1800     integral += fP2Step*gHistBalance->GetBinContent(i);
1801   }
1802   for(Int_t i = fStartBin; i < fNumberOfBins; i++)
1803     deltaErrorNew += gHistBalance->GetBinError(i)*(currentBinCenter*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
1804
1805   Double_t integralError = TMath::Sqrt(deltaBalP2);
1806
1807   Double_t delta = gSumBiXi / gSumBi;
1808   Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
1809
1810   WM  = delta;
1811   WME = deltaError;
1812 }
1813