PWGCF/EBYE/macros/drawBalanceFunctionPsi.C

   1 const Int_t numberOfCentralityBins = 11;
   2 TString centralityArray[numberOfCentralityBins] = {"0-10","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-1","1-2","0-100"};
   3
   4 void drawBalanceFunctionPsi(const char* filename = "AnalysisResultsPsi.root",
   5                             Int_t gCentrality = 1,
   6                             Int_t gDeltaEtaDeltaPhi = 2,
   7                             Int_t gBit = -1,
   8                             const char* gCentralityEstimator = 0x0,
   9                             Double_t psiMin = -0.5, Double_t psiMax = 0.5,
  10                             Double_t ptTriggerMin = -1.,
  11                             Double_t ptTriggerMax = -1.,
  12                             Double_t ptAssociatedMin = -1.,
  13                             Double_t ptAssociatedMax = -1.,
  14                             Bool_t k2pMethod = kFALSE,
  15                             Bool_t k2pMethod2D = kFALSE,
  16                             TString eventClass = "EventPlane") //Can be "EventPlane", "Centrality", "Multiplicity"
  17 {
  18   //Macro that draws the BF distributions for each centrality bin
  19   //for reaction plane dependent analysis
  20   //Author: Panos.Christakoglou@nikhef.nl
  21   //Load the PWG2ebye library
  22   gSystem->Load("libANALYSIS.so");
  23   gSystem->Load("libANALYSISalice.so");
  24   gSystem->Load("libEventMixing.so");
  25   gSystem->Load("libCORRFW.so");
  26   gSystem->Load("libPWGTools.so");
  27   gSystem->Load("libPWGCFebye.so");
  28
  29   //correction method check
  30   if(k2pMethod2D&&!k2pMethod){
  31     Printf("Chosen 2D 2particle correction method w/o 2particle correction --> not possible");
  32     return;
  33   }
  34
  35   //Prepare the objects and return them
  36   TList *listBF = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,0);
  37   TList *listBFShuffled = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,1);
  38   TList *listBFMixed = GetListOfObjects(filename,gCentrality,gBit,gCentralityEstimator,2);
  39   if(!listBF) {
  40     Printf("The TList object was not created");
  41     return;
  42   }
  43   else
  44     draw(listBF,listBFShuffled,listBFMixed,
  45          gCentrality,gDeltaEtaDeltaPhi,
  46          psiMin,psiMax,
  47          ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,
  48          k2pMethod,k2pMethod2D,eventClass);
  49 }
  50
  51 //______________________________________________________//
  52 TList *GetListOfObjects(const char* filename,
  53                         Int_t gCentrality,
  54                         Int_t gBit,
  55                         const char *gCentralityEstimator,
  56                         Int_t kData = 1) {
  57   //Get the TList objects (QA, bf, bf shuffled)
  58   TList *listBF = 0x0;
  59
  60   //Open the file
  61   TFile *f = TFile::Open(filename,"UPDATE");
  62   if((!f)||(!f->IsOpen())) {
  63     Printf("The file %s is not found. Aborting...",filename);
  64     return listBF;
  65   }
  66   //f->ls();
  67
  68   TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(f->Get("PWGCFEbyE.outputBalanceFunctionPsiAnalysis"));
  69   if(!dir) {
  70     Printf("The TDirectoryFile is not found. Aborting...",filename);
  71     return listBF;
  72   }
  73   //dir->ls();
  74
  75   TString listBFName;
  76   if(kData == 0) {
  77     //cout<<"no shuffling - no mixing"<<endl;
  78     listBFName = "listBFPsi_";
  79   }
  80   else if(kData == 1) {
  81     //cout<<"shuffling - no mixing"<<endl;
  82     listBFName = "listBFPsiShuffled_";
  83   }
  84   else if(kData == 2) {
  85     //cout<<"no shuffling - mixing"<<endl;
  86     listBFName = "listBFPsiMixed_";
  87   }
  88   listBFName += centralityArray[gCentrality-1];
  89   if(gBit > -1) {
  90     listBFName += "_Bit"; listBFName += gBit; }
  91   if(gCentralityEstimator) {
  92     listBFName += "_"; listBFName += gCentralityEstimator;}
  93
  94   // histograms were already retrieved (in first iteration)
  95   if(dir->Get(Form("%s_histograms",listBFName.Data()))){
  96     listBF = dynamic_cast<TList *>(dir->Get(Form("%s_histograms",listBFName.Data())));
  97   }
  98
  99   // histograms were not yet retrieved (this is the first iteration)
 100   else{
 101
 102     listBF = dynamic_cast<TList *>(dir->Get(listBFName.Data()));
 103     cout<<"======================================================="<<endl;
 104     cout<<"List name (check): "<<listBFName.Data()<<endl;
 105     cout<<"List name: "<<listBF->GetName()<<endl;
 106     //listBF->ls();
 107
 108     //Get the histograms
 109     TString histoName;
 110     if(kData == 0)
 111       histoName = "fHistPV0M";
 112     else if(kData == 1)
 113       histoName = "fHistP_shuffleV0M";
 114     else if(kData == 2)
 115       histoName = "fHistPV0M";
 116     AliTHn *fHistP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 117     if(!fHistP) {
 118       Printf("fHistP %s not found!!!",histoName.Data());
 119       break;
 120     }
 121     fHistP->FillParent(); fHistP->DeleteContainers();
 122
 123     if(kData == 0)
 124       histoName = "fHistNV0M";
 125     if(kData == 1)
 126       histoName = "fHistN_shuffleV0M";
 127     if(kData == 2)
 128       histoName = "fHistNV0M";
 129     AliTHn *fHistN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 130     if(!fHistN) {
 131       Printf("fHistN %s not found!!!",histoName.Data());
 132       break;
 133     }
 134     fHistN->FillParent(); fHistN->DeleteContainers();
 135
 136     if(kData == 0)
 137       histoName = "fHistPNV0M";
 138     if(kData == 1)
 139       histoName = "fHistPN_shuffleV0M";
 140     if(kData == 2)
 141       histoName = "fHistPNV0M";
 142     AliTHn *fHistPN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 143     if(!fHistPN) {
 144       Printf("fHistPN %s not found!!!",histoName.Data());
 145       break;
 146     }
 147     fHistPN->FillParent(); fHistPN->DeleteContainers();
 148
 149     if(kData == 0)
 150       histoName = "fHistNPV0M";
 151     if(kData == 1)
 152       histoName = "fHistNP_shuffleV0M";
 153     if(kData == 2)
 154       histoName = "fHistNPV0M";
 155     AliTHn *fHistNP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 156     if(!fHistNP) {
 157       Printf("fHistNP %s not found!!!",histoName.Data());
 158       break;
 159     }
 160     fHistNP->FillParent(); fHistNP->DeleteContainers();
 161
 162     if(kData == 0)
 163       histoName = "fHistPPV0M";
 164     if(kData == 1)
 165       histoName = "fHistPP_shuffleV0M";
 166     if(kData == 2)
 167       histoName = "fHistPPV0M";
 168     AliTHn *fHistPP = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 169     if(!fHistPP) {
 170       Printf("fHistPP %s not found!!!",histoName.Data());
 171       break;
 172     }
 173     fHistPP->FillParent(); fHistPP->DeleteContainers();
 174
 175     if(kData == 0)
 176       histoName = "fHistNNV0M";
 177     if(kData == 1)
 178       histoName = "fHistNN_shuffleV0M";
 179     if(kData == 2)
 180       histoName = "fHistNNV0M";
 181     AliTHn *fHistNN = dynamic_cast<AliTHn *>(listBF->FindObject(histoName.Data()));
 182     if(!fHistNN) {
 183       Printf("fHistNN %s not found!!!",histoName.Data());
 184       break;
 185     }
 186     fHistNN->FillParent(); fHistNN->DeleteContainers();
 187
 188     dir->cd();
 189     listBF->Write(Form("%s_histograms",listBFName.Data()), TObject::kSingleKey);
 190
 191   }// first iteration
 192
 193   f->Close();
 194
 195   return listBF;
 196 }
 197
 198 //______________________________________________________//
 199 void draw(TList *listBF, TList *listBFShuffled, TList *listBFMixed,
 200           Int_t gCentrality, Int_t gDeltaEtaDeltaPhi,
 201           Double_t psiMin, Double_t psiMax,
 202           Double_t ptTriggerMin, Double_t ptTriggerMax,
 203           Double_t ptAssociatedMin, Double_t ptAssociatedMax,
 204           Bool_t k2pMethod = kFALSE,Bool_t k2pMethod2D = kFALSE, TString eventClass="EventPlane") {
 205   gROOT->LoadMacro("~/SetPlotStyle.C");
 206   SetPlotStyle();
 207   gStyle->SetPalette(1,0);
 208
 209   const Int_t gRebin = gDeltaEtaDeltaPhi; //rebin by 2 the Delta phi projection
 210
 211   //balance function
 212   AliTHn *hP = NULL;
 213   AliTHn *hN = NULL;
 214   AliTHn *hPN = NULL;
 215   AliTHn *hNP = NULL;
 216   AliTHn *hPP = NULL;
 217   AliTHn *hNN = NULL;
 218   //listBF->ls();
 219   //Printf("=================");
 220   hP = (AliTHn*) listBF->FindObject("fHistPV0M");
 221   hN = (AliTHn*) listBF->FindObject("fHistNV0M");
 222   hPN = (AliTHn*) listBF->FindObject("fHistPNV0M");
 223   hNP = (AliTHn*) listBF->FindObject("fHistNPV0M");
 224   hPP = (AliTHn*) listBF->FindObject("fHistPPV0M");
 225   hNN = (AliTHn*) listBF->FindObject("fHistNNV0M");
 226
 227   AliBalancePsi *b = new AliBalancePsi();
 228   b->SetEventClass(eventClass);
 229   b->SetHistNp(hP);
 230   b->SetHistNn(hN);
 231   b->SetHistNpn(hPN);
 232   b->SetHistNnp(hNP);
 233   b->SetHistNpp(hPP);
 234   b->SetHistNnn(hNN);
 235
 236   //balance function shuffling
 237   AliTHn *hPShuffled = NULL;
 238   AliTHn *hNShuffled = NULL;
 239   AliTHn *hPNShuffled = NULL;
 240   AliTHn *hNPShuffled = NULL;
 241   AliTHn *hPPShuffled = NULL;
 242   AliTHn *hNNShuffled = NULL;
 243   //listBFShuffled->ls();
 244   hPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistP_shuffleV0M");
 245   hNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistN_shuffleV0M");
 246   hPNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPN_shuffleV0M");
 247   hNPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNP_shuffleV0M");
 248   hPPShuffled = (AliTHn*) listBFShuffled->FindObject("fHistPP_shuffleV0M");
 249   hNNShuffled = (AliTHn*) listBFShuffled->FindObject("fHistNN_shuffleV0M");
 250
 251   AliBalancePsi *bShuffled = new AliBalancePsi();
 252   bShuffled->SetEventClass(eventClass);
 253   bShuffled->SetHistNp(hPShuffled);
 254   bShuffled->SetHistNn(hNShuffled);
 255   bShuffled->SetHistNpn(hPNShuffled);
 256   bShuffled->SetHistNnp(hNPShuffled);
 257   bShuffled->SetHistNpp(hPPShuffled);
 258   bShuffled->SetHistNnn(hNNShuffled);
 259
 260   //balance function mixing
 261   AliTHn *hPMixed = NULL;
 262   AliTHn *hNMixed = NULL;
 263   AliTHn *hPNMixed = NULL;
 264   AliTHn *hNPMixed = NULL;
 265   AliTHn *hPPMixed = NULL;
 266   AliTHn *hNNMixed = NULL;
 267   //listBFMixed->ls();
 268   hPMixed = (AliTHn*) listBFMixed->FindObject("fHistPV0M");
 269   hNMixed = (AliTHn*) listBFMixed->FindObject("fHistNV0M");
 270   hPNMixed = (AliTHn*) listBFMixed->FindObject("fHistPNV0M");
 271   hNPMixed = (AliTHn*) listBFMixed->FindObject("fHistNPV0M");
 272   hPPMixed = (AliTHn*) listBFMixed->FindObject("fHistPPV0M");
 273   hNNMixed = (AliTHn*) listBFMixed->FindObject("fHistNNV0M");
 274
 275   AliBalancePsi *bMixed = new AliBalancePsi();
 276   bMixed->SetEventClass(eventClass);
 277   bMixed->SetHistNp(hPMixed);
 278   bMixed->SetHistNn(hNMixed);
 279   bMixed->SetHistNpn(hPNMixed);
 280   bMixed->SetHistNnp(hNPMixed);
 281   bMixed->SetHistNpp(hPPMixed);
 282   bMixed->SetHistNnn(hNNMixed);
 283
 284   TH1D *gHistBalanceFunction;
 285   TH1D *gHistBalanceFunctionShuffled;
 286   TH1D *gHistBalanceFunctionMixed;
 287   TH1D *gHistBalanceFunctionSubtracted;
 288   TString histoTitle, pngName;
 289   TLegend *legend;
 290
 291   if(eventClass == "Centrality"){
 292     histoTitle = "Centrality: ";
 293     histoTitle += psiMin;
 294     histoTitle += " - ";
 295     histoTitle += psiMax;
 296     histoTitle += " % ";
 297     histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 298   }
 299   else if(eventClass == "Multiplicity"){
 300     histoTitle = "Multiplicity: ";
 301     histoTitle += psiMin;
 302     histoTitle += " - ";
 303     histoTitle += psiMax;
 304     histoTitle += " tracks";
 305     histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 306   }
 307   else{ // "EventPlane" (default)
 308     histoTitle = "Centrality: ";
 309     histoTitle += centralityArray[gCentrality-1];
 310     histoTitle += "%";
 311     if((psiMin == -0.5)&&(psiMax == 0.5))
 312       histoTitle += " (-7.5^{o} < #varphi - #Psi_{2} < 7.5^{o})";
 313     else if((psiMin == 0.5)&&(psiMax == 1.5))
 314       histoTitle += " (37.5^{o} < #varphi - #Psi_{2} < 52.5^{o})";
 315     else if((psiMin == 1.5)&&(psiMax == 2.5))
 316       histoTitle += " (82.5^{o} < #varphi - #Psi_{2} < 97.5^{o})";
 317     else
 318       histoTitle += " (0^{o} < #varphi - #Psi_{2} < 180^{o})";
 319   }
 320
 321   //Raw balance function
 322   if(k2pMethod){
 323     if(bMixed){
 324       gHistBalanceFunction = b->GetBalanceFunctionHistogram2pMethod(0,gDeltaEtaDeltaPhi,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 325     }
 326     else{
 327       cerr<<"RAW: NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
 328       return;
 329     }
 330   }
 331   else if(k2pMethod2D){
 332     if(bMixed){
 333       if(gDeltaEtaDeltaPhi==1) //Delta eta
 334         gHistBalanceFunction = b->GetBalanceFunction1DFrom2D2pMethod(0,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 335       else //Delta phi
 336         gHistBalanceFunction = b->GetBalanceFunction1DFrom2D2pMethod(1,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 337     }
 338     else{
 339       cerr<<"RAW: NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
 340       return;
 341     }
 342   }
 343   else
 344     gHistBalanceFunction = b->GetBalanceFunctionHistogram(0,gDeltaEtaDeltaPhi,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 345   gHistBalanceFunction->SetMarkerStyle(20);
 346   gHistBalanceFunction->SetTitle(histoTitle.Data());
 347   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
 348   gHistBalanceFunction->SetName("gHistBalanceFunction");
 349
 350   //Shuffled balance function
 351   if(k2pMethod){
 352     if(bMixed)
 353       gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionHistogram2pMethod(0,gDeltaEtaDeltaPhi,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 354     else{
 355       cerr<<"SHUFFLE: NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
 356       return;
 357     }
 358   }
 359   else if(k2pMethod2D){
 360     if(bMixed){
 361       if(gDeltaEtaDeltaPhi==1) //Delta eta
 362         gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunction1DFrom2D2pMethod(0,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 363       else //Delta phi
 364         gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunction1DFrom2D2pMethod(1,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 365     }
 366     else{
 367       cerr<<"SHUFFLE: NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
 368       return;
 369     }
 370   }
 371   else
 372     gHistBalanceFunctionShuffled = bShuffled->GetBalanceFunctionHistogram(0,gDeltaEtaDeltaPhi,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 373   gHistBalanceFunctionShuffled->SetMarkerStyle(24);
 374   gHistBalanceFunctionShuffled->SetName("gHistBalanceFunctionShuffled");
 375
 376   //Mixed balance function
 377   if(k2pMethod){
 378     if(bMixed)
 379       gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionHistogram2pMethod(0,gDeltaEtaDeltaPhi,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 380     else{
 381       cerr<<"MIXED: NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
 382       return;
 383     }
 384   }
 385   else if(k2pMethod2D){
 386     if(bMixed){
 387       if(gDeltaEtaDeltaPhi==1) //Delta eta
 388         gHistBalanceFunctionMixed = bMixed->GetBalanceFunction1DFrom2D2pMethod(0,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 389       else //Delta phi
 390         gHistBalanceFunctionMixed = bMixed->GetBalanceFunction1DFrom2D2pMethod(1,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax,bMixed);
 391     }
 392     else{
 393       cerr<<"MIXED: NO MIXED BF BUT REQUESTED CORRECTING WITH IT! --> FAIL"<<endl;
 394       return;
 395     }
 396   }
 397   else
 398     gHistBalanceFunctionMixed = bMixed->GetBalanceFunctionHistogram(0,gDeltaEtaDeltaPhi,psiMin,psiMax,ptTriggerMin,ptTriggerMax,ptAssociatedMin,ptAssociatedMax);
 399   gHistBalanceFunctionMixed->SetMarkerStyle(25);
 400   gHistBalanceFunctionMixed->SetName("gHistBalanceFunctionMixed");
 401
 402   //Subtracted balance function
 403   gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunction->Clone());
 404   gHistBalanceFunctionSubtracted->Add(gHistBalanceFunctionMixed,-1);
 405   gHistBalanceFunctionSubtracted->Rebin(gRebin);
 406   gHistBalanceFunctionSubtracted->Scale(1./(Double_t)(gRebin));
 407   gHistBalanceFunctionSubtracted->SetMarkerStyle(20);
 408   gHistBalanceFunctionSubtracted->SetTitle(histoTitle.Data());
 409   gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
 410   gHistBalanceFunctionSubtracted->SetName("gHistBalanceFunctionSubtracted");
 411
 412   TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
 413   c1->SetFillColor(10);
 414   c1->SetHighLightColor(10);
 415   c1->SetLeftMargin(0.15);
 416   gHistBalanceFunction->DrawCopy("E");
 417   gHistBalanceFunctionShuffled->DrawCopy("ESAME");
 418   gHistBalanceFunctionMixed->DrawCopy("ESAME");
 419
 420   legend = new TLegend(0.18,0.62,0.45,0.82,"","brNDC");
 421   legend->SetTextSize(0.045);
 422   legend->SetTextFont(42);
 423   legend->SetBorderSize(0);
 424   legend->SetFillStyle(0);
 425   legend->SetFillColor(10);
 426   legend->SetMargin(0.25);
 427   legend->SetShadowColor(10);
 428   legend->AddEntry(gHistBalanceFunction,"Data","lp");
 429   legend->AddEntry(gHistBalanceFunctionShuffled,"Shuffled data","lp");
 430   legend->AddEntry(gHistBalanceFunctionMixed,"Mixed data","lp");
 431   legend->Draw();
 432
 433   pngName = "BalanceFunction.";
 434   if(eventClass == "Centrality"){
 435     pngName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
 436     if(gDeltaEtaDeltaPhi == 1) pngName += ".InDeltaEta.PsiAll.PttFrom";
 437     else if(gDeltaEtaDeltaPhi == 2) pngName += ".InDeltaPhi.PsiAll.PttFrom";
 438   }
 439   else if(eventClass == "Multiplicity"){
 440     pngName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
 441     if(gDeltaEtaDeltaPhi == 1) pngName += ".InDeltaEta.PsiAll.PttFrom";
 442     else if(gDeltaEtaDeltaPhi == 2) pngName += ".InDeltaPhi.PsiAll.PttFrom";
 443   }
 444   else{ // "EventPlane" (default)
 445     pngName += "Centrality";
 446     pngName += gCentrality;
 447     if(gDeltaEtaDeltaPhi == 1) pngName += ".InDeltaEta.Psi";
 448     else if(gDeltaEtaDeltaPhi == 2) pngName += ".InDeltaPhi.Psi";
 449     if((psiMin == -0.5)&&(psiMax == 0.5)) pngName += "InPlane.Ptt";
 450     else if((psiMin == 0.5)&&(psiMax == 1.5)) pngName += "Intermediate.Ptt";
 451     else if((psiMin == 1.5)&&(psiMax == 2.5)) pngName += "OutOfPlane.Ptt";
 452     else if((psiMin == 2.5)&&(psiMax == 3.5)) pngName += "Rest.PttFrom";
 453     else pngName += "All.PttFrom";
 454   }
 455   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
 456   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
 457   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
 458   pngName += Form("%.1f",ptAssociatedMax);
 459   if(k2pMethod2D) pngName += "_2pMethod2D";
 460   else if(k2pMethod) pngName += "_2pMethod";
 461   pngName += ".png";
 462
 463   c1->SaveAs(pngName.Data());
 464
 465   GetWeightedMean(gHistBalanceFunction);
 466   GetWeightedMean(gHistBalanceFunctionShuffled);
 467
 468   TString meanLatex, rmsLatex, skewnessLatex, kurtosisLatex;
 469   meanLatex = "#mu = ";
 470   meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMean());
 471   meanLatex += " #pm ";
 472   meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMeanError());
 473
 474   rmsLatex = "#sigma = ";
 475   rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMS());
 476   rmsLatex += " #pm ";
 477   rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMSError());
 478
 479   skewnessLatex = "S = ";
 480   skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(1));
 481   skewnessLatex += " #pm ";
 482   skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(11));
 483
 484   kurtosisLatex = "K = ";
 485   kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(1));
 486   kurtosisLatex += " #pm ";
 487   kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(11));
 488   Printf("Mean: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetMean(),gHistBalanceFunctionSubtracted->GetMeanError());
 489   Printf("RMS: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetRMS(),gHistBalanceFunctionSubtracted->GetRMSError());
 490   Printf("Skeweness: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetSkewness(1),gHistBalanceFunctionSubtracted->GetSkewness(11));
 491   Printf("Kurtosis: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetKurtosis(1),gHistBalanceFunctionSubtracted->GetKurtosis(11));
 492
 493   TCanvas *c2 = new TCanvas("c2","",600,0,600,500);
 494   c2->SetFillColor(10);
 495   c2->SetHighLightColor(10);
 496   c2->SetLeftMargin(0.15);
 497   gHistBalanceFunctionSubtracted->DrawCopy("E");
 498
 499   TLatex *latex = new TLatex();
 500   latex->SetNDC();
 501   latex->SetTextSize(0.035);
 502   latex->SetTextColor(1);
 503   latex->DrawLatex(0.64,0.85,meanLatex.Data());
 504   latex->DrawLatex(0.64,0.81,rmsLatex.Data());
 505   latex->DrawLatex(0.64,0.77,skewnessLatex.Data());
 506   latex->DrawLatex(0.64,0.73,kurtosisLatex.Data());
 507
 508   TString newFileName = "balanceFunction.";
 509   if(eventClass == "Centrality"){
 510     newFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
 511     if(gDeltaEtaDeltaPhi == 1) newFileName += ".InDeltaEta.PsiAll.PttFrom";
 512     else if(gDeltaEtaDeltaPhi == 2) newFileName += ".InDeltaPhi.PsiAll.PttFrom";
 513   }
 514   else if(eventClass == "Multiplicity"){
 515     newFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
 516     if(gDeltaEtaDeltaPhi == 1) newFileName += ".InDeltaEta.PsiAll.PttFrom";
 517     else if(gDeltaEtaDeltaPhi == 2) newFileName += ".InDeltaPhi.PsiAll.PttFrom";
 518   }
 519   else{ // "EventPlane" (default)
 520     newFileName += "Centrality";
 521     newFileName += gCentrality;
 522     if(gDeltaEtaDeltaPhi == 1) newFileName += ".InDeltaEta.Psi";
 523     else if(gDeltaEtaDeltaPhi == 2) newFileName += ".InDeltaPhi.Psi";
 524     if((psiMin == -0.5)&&(psiMax == 0.5)) newFileName += "InPlane.Ptt";
 525     else if((psiMin == 0.5)&&(psiMax == 1.5)) newFileName += "Intermediate.Ptt";
 526     else if((psiMin == 1.5)&&(psiMax == 2.5)) newFileName += "OutOfPlane.Ptt";
 527     else if((psiMin == 2.5)&&(psiMax == 3.5)) newFileName += "Rest.PttFrom";
 528     else newFileName += "All.PttFrom";
 529   }
 530   newFileName += Form("%.1f",ptTriggerMin); newFileName += "To";
 531   newFileName += Form("%.1f",ptTriggerMax); newFileName += "PtaFrom";
 532   newFileName += Form("%.1f",ptAssociatedMin); newFileName += "To";
 533   newFileName += Form("%.1f",ptAssociatedMax);
 534   if(k2pMethod2D) newFileName += "_2pMethod2D";
 535   else if(k2pMethod) newFileName += "_2pMethod";
 536   newFileName += ".root";
 537
 538   TFile *fOutput = new TFile(newFileName.Data(),"recreate");
 539   fOutput->cd();
 540   gHistBalanceFunction->Write();
 541   gHistBalanceFunctionShuffled->Write();
 542   gHistBalanceFunctionMixed->Write();
 543   gHistBalanceFunctionSubtracted->Write();
 544   fOutput->Close();
 545 }
 546
 547 //____________________________________________________________________//
 548 void GetWeightedMean(TH1D *gHistBalance, Int_t fStartBin = 1) {
 549   //Prints the calculated width of the BF and its error
 550   Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
 551   Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
 552   Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
 553   Double_t deltaBalP2 = 0.0, integral = 0.0;
 554   Double_t deltaErrorNew = 0.0;
 555
 556   //Retrieve this variables from Histogram
 557   Int_t fNumberOfBins = gHistBalance->GetNbinsX();
 558   Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!
 559
 560   //cout<<"=================================================="<<endl;
 561   //cout<<"RECALCULATION OF BF WIDTH (StartBin = "<<fStartBin<<")"<<endl;
 562   //cout<<"HISTOGRAM has "<<fNumberOfBins<<" bins with bin size of "<<fP2Step<<endl;
 563   //cout<<"=================================================="<<endl;
 564   for(Int_t i = 1; i <= fNumberOfBins; i++) {
 565     // this is to simulate |Delta eta| or |Delta phi|
 566     if(fNumberOfBins/2 - fStartBin + 1 < i && i < fNumberOfBins/2 + fStartBin ) continue;
 567
 568     //cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<TMath::Abs(gHistBalance->GetBinCenter(i))<<endl;
 569
 570     gSumXi += TMath::Abs(gHistBalance->GetBinCenter(i)); // this is to simulate |Delta eta| or |Delta phi|
 571     gSumBi += gHistBalance->GetBinContent(i);
 572     gSumBiXi += gHistBalance->GetBinContent(i)*TMath::Abs(gHistBalance->GetBinCenter(i));
 573     gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2);
 574     gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2);
 575     gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
 576     gSumXi2DeltaBi2 += TMath::Power(TMath::Abs(gHistBalance->GetBinCenter(i)),2) * TMath::Power(gHistBalance->GetBinError(i),2);
 577
 578     deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
 579     integral += fP2Step*gHistBalance->GetBinContent(i);
 580   }
 581   for(Int_t i = fStartBin; i < fNumberOfBins; i++)
 582     deltaErrorNew += gHistBalance->GetBinError(i)*(TMath::Abs(gHistBalance->GetBinCenter(i))*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
 583
 584   Double_t integralError = TMath::Sqrt(deltaBalP2);
 585
 586   Double_t delta = gSumBiXi / gSumBi;
 587   Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
 588   cout<<"=================================================="<<endl;
 589   cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
 590   cout<<"New error: "<<deltaErrorNew<<endl;
 591   cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
 592   cout<<"=================================================="<<endl;
 593   cout<<endl;
 594 }
 595
 596 //______________________________________________________//
 597 void drawBFPsi(const char* lhcPeriod = "LHC10h",
 598                Int_t gCentrality = 1,
 599                Int_t gDeltaEtaDeltaPhi = 2,
 600                Double_t psiMin = -0.5, Double_t psiMax = 0.5,
 601                Double_t ptTriggerMin = -1.,
 602                Double_t ptTriggerMax = -1.,
 603                Double_t ptAssociatedMin = -1.,
 604                Double_t ptAssociatedMax = -1.) {
 605   //Macro that draws the BF distributions for each centrality bin
 606   //for reaction plane dependent analysis
 607   //Author: Panos.Christakoglou@nikhef.nl
 608   gROOT->LoadMacro("~/SetPlotStyle.C");
 609   SetPlotStyle();
 610
 611   //Get the input file
 612   TString filename = lhcPeriod; filename +="/PttFrom";
 613   filename += Form("%.1f",ptTriggerMin); filename += "To";
 614   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
 615   filename += Form("%.1f",ptAssociatedMin); filename += "To";
 616   filename += Form("%.1f",ptAssociatedMax);
 617   filename += "/balanceFunction.Centrality";
 618   filename += gCentrality; filename += ".In";
 619   if(gDeltaEtaDeltaPhi == 1) filename += "DeltaEta.Psi";
 620   else if(gDeltaEtaDeltaPhi == 2) filename += "DeltaPhi.Psi";
 621   if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
 622   else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
 623   else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
 624   else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.Ptt";
 625   else filename += "All.PttFrom";
 626   filename += Form("%.1f",ptTriggerMin); filename += "To";
 627   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
 628   filename += Form("%.1f",ptAssociatedMin); filename += "To";
 629   filename += Form("%.1f",ptAssociatedMax);  filename += ".root";
 630
 631   //Open the file
 632   TFile *f = TFile::Open(filename.Data());
 633   if((!f)||(!f->IsOpen())) {
 634     Printf("The file %s is not found. Aborting...",filename);
 635     return listBF;
 636   }
 637   //f->ls();
 638
 639   //Raw balance function
 640   TH1D *gHistBalanceFunction = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunction"));
 641   gHistBalanceFunction->SetStats(kFALSE);
 642   gHistBalanceFunction->SetMarkerStyle(20);
 643   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.3);
 644   if(gDeltaEtaDeltaPhi == 2) {
 645     gHistBalanceFunction->GetYaxis()->SetTitle("B(#Delta #varphi)");
 646     gHistBalanceFunction->GetXaxis()->SetTitle("#Delta #varphi (deg.)");
 647   }
 648
 649   //Shuffled balance function
 650   TH1D *gHistBalanceFunctionShuffled = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionShuffled"));
 651   gHistBalanceFunction->SetStats(kFALSE);
 652   gHistBalanceFunctionShuffled->SetMarkerStyle(24);
 653
 654   //Mixed balance function
 655   TH1D *gHistBalanceFunctionMixed = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionMixed"));
 656   gHistBalanceFunction->SetStats(kFALSE);
 657   gHistBalanceFunctionMixed->SetMarkerStyle(25);
 658
 659   //Subtracted balance function
 660   TH1D *gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(f->Get("gHistBalanceFunctionSubtracted"));
 661   gHistBalanceFunctionSubtracted->SetStats(kFALSE);
 662   gHistBalanceFunctionSubtracted->SetMarkerStyle(20);
 663   gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
 664   if(gDeltaEtaDeltaPhi == 2) {
 665     gHistBalanceFunctionSubtracted->GetYaxis()->SetTitle("B(#Delta #varphi)");
 666     gHistBalanceFunctionSubtracted->GetXaxis()->SetTitle("#Delta #varphi (deg.)");
 667   }
 668
 669   TString pngName;
 670   TLegend *legend;
 671
 672   TString centralityLatex = "Centrality: ";
 673   centralityLatex += centralityArray[gCentrality-1];
 674   centralityLatex += "%";
 675
 676   TString psiLatex;
 677   if((psiMin == -0.5)&&(psiMax == 0.5))
 678     psiLatex = " -7.5^{o} < #varphi - #Psi_{2} < 7.5^{o}";
 679   else if((psiMin == 0.5)&&(psiMax == 1.5))
 680     psiLatex = " 37.5^{o} < #varphi - #Psi_{2} < 52.5^{o}";
 681   else if((psiMin == 1.5)&&(psiMax == 2.5))
 682     psiLatex = " 82.5^{o} < #varphi - #Psi_{2} < 97.5^{o}";
 683   else
 684     psiLatex = " 0^{o} < #varphi - #Psi_{2} < 180^{o}";
 685
 686   TString pttLatex = Form("%.1f",ptTriggerMin);
 687   pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
 688   pttLatex += " GeV/c";
 689
 690   TString ptaLatex = Form("%.1f",ptAssociatedMin);
 691   ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
 692   ptaLatex += " GeV/c";
 693
 694   //Draw the results
 695   TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
 696   c1->SetFillColor(10); c1->SetHighLightColor(10);
 697   c1->SetLeftMargin(0.17); c1->SetTopMargin(0.05);
 698   gHistBalanceFunction->SetTitle("");
 699   gHistBalanceFunction->GetYaxis()->SetTitleOffset(1.4);
 700   gHistBalanceFunction->GetYaxis()->SetNdivisions(10);
 701   gHistBalanceFunction->GetXaxis()->SetNdivisions(10);
 702   gHistBalanceFunction->DrawCopy("E");
 703   gHistBalanceFunctionShuffled->DrawCopy("ESAME");
 704   gHistBalanceFunctionMixed->DrawCopy("ESAME");
 705
 706   legend = new TLegend(0.2,0.72,0.45,0.92,"","brNDC");
 707   legend->SetTextSize(0.05);
 708   legend->SetTextFont(42);
 709   legend->SetBorderSize(0);
 710   legend->SetFillStyle(0);
 711   legend->SetFillColor(10);
 712   legend->SetMargin(0.25);
 713   legend->SetShadowColor(10);
 714   legend->AddEntry(gHistBalanceFunction,"Data","lp");
 715   legend->AddEntry(gHistBalanceFunctionShuffled,"Shuffled data","lp");
 716   legend->AddEntry(gHistBalanceFunctionMixed,"Mixed data","lp");
 717   legend->Draw();
 718
 719   TLatex *latexInfo1 = new TLatex();
 720   latexInfo1->SetNDC();
 721   latexInfo1->SetTextSize(0.04);
 722   latexInfo1->SetTextColor(1);
 723   latexInfo1->DrawLatex(0.58,0.88,centralityLatex.Data());
 724   latexInfo1->DrawLatex(0.58,0.82,psiLatex.Data());
 725   latexInfo1->DrawLatex(0.58,0.76,pttLatex.Data());
 726   latexInfo1->DrawLatex(0.58,0.70,ptaLatex.Data());
 727
 728   //pngName = "BalanceFunctionDeltaEta.Centrality";
 729   //pngName += centralityArray[gCentrality-1];
 730   //pngName += ".Psi"; //pngName += psiMin; pngName += "To"; pngName += psiMax;
 731   //pngName += ".png";
 732   //c1->SaveAs(pngName.Data());
 733
 734   TCanvas *c2 = new TCanvas("c2","",600,0,600,500);
 735   c2->SetFillColor(10); c2->SetHighLightColor(10);
 736   c2->SetLeftMargin(0.17); c2->SetTopMargin(0.05);
 737   gHistBalanceFunctionSubtracted->SetTitle("");
 738   gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.4);
 739   gHistBalanceFunctionSubtracted->GetYaxis()->SetNdivisions(10);
 740   gHistBalanceFunctionSubtracted->GetXaxis()->SetNdivisions(10);
 741   gHistBalanceFunctionSubtracted->DrawCopy("E");
 742
 743   //Opening the output ascii files
 744   TString filenameMean = "meanIn";
 745   TString filenameSigma = "sigmaIn";
 746   TString filenameSkewness = "skewnessIn";
 747   TString filenameKurtosis = "kurtosisIn";
 748   if(gDeltaEtaDeltaPhi == 1) {
 749     filenameMean += "DeltaEta.Psi"; filenameSigma += "DeltaEta.Psi";
 750     filenameSkewness += "DeltaEta.Psi"; filenameKurtosis += "DeltaEta.Psi";
 751   }
 752   else if(gDeltaEtaDeltaPhi == 2) {
 753     filenameMean += "DeltaPhi.Psi"; filenameSigma += "DeltaPhi.Psi";
 754     filenameSkewness += "DeltaPhi.Psi"; filenameKurtosis += "DeltaPhi.Psi";
 755   }
 756   if((psiMin == -0.5)&&(psiMax == 0.5)) {
 757     filenameMean += "InPlane.Ptt"; filenameSigma += "InPlane.Ptt";
 758     filenameSkewness += "InPlane.Ptt"; filenameKurtosis += "InPlane.Ptt";
 759   }
 760   else if((psiMin == 0.5)&&(psiMax == 1.5)) {
 761     filenameMean += "Intermediate.Ptt"; filenameSigma += "Intermediate.Ptt";
 762     filenameSkewness += "Intermediate.Ptt";
 763     filenameKurtosis += "Intermediate.Ptt";
 764   }
 765   else if((psiMin == 1.5)&&(psiMax == 2.5)) {
 766     filenameMean += "OutOfPlane.Ptt"; filenameSigma += "OutOfPlane.Ptt";
 767     filenameSkewness += "OutOfPlane.Ptt";
 768     filenameKurtosis += "OutOfPlane.Ptt";
 769   }
 770   else if((psiMin == 2.5)&&(psiMax == 3.5)) {
 771     filenameMean += "Rest.Ptt"; filenameSigma += "Rest.Ptt";
 772     filenameSkewness += "Rest.Ptt"; filenameKurtosis += "Rest.Ptt";
 773   }
 774   else {
 775     filenameMean += "All.Ptt"; filenameSigma += "All.Ptt";
 776     filenameSkewness += "All.Ptt"; filenameKurtosis += "All.Ptt";
 777   }
 778   filenameMean += Form("%.1f",ptTriggerMin); filenameMean += "To";
 779   filenameMean += Form("%.1f",ptTriggerMax); filenameMean += "PtaFrom";
 780   filenameMean += Form("%.1f",ptAssociatedMin); filenameMean += "To";
 781   filenameMean += Form("%.1f",ptAssociatedMax);  filenameMean += ".txt";
 782   filenameSigma += Form("%.1f",ptTriggerMin); filenameSigma += "To";
 783   filenameSigma += Form("%.1f",ptTriggerMax); filenameSigma += "PtaFrom";
 784   filenameSigma += Form("%.1f",ptAssociatedMin); filenameSigma += "To";
 785   filenameSigma += Form("%.1f",ptAssociatedMax);  filenameSigma += ".txt";
 786   filenameSkewness += Form("%.1f",ptTriggerMin); filenameSkewness += "To";
 787   filenameSkewness += Form("%.1f",ptTriggerMax); filenameSkewness += "PtaFrom";
 788   filenameSkewness += Form("%.1f",ptAssociatedMin); filenameSkewness += "To";
 789   filenameSkewness += Form("%.1f",ptAssociatedMax);
 790   filenameSkewness += ".txt";
 791   filenameKurtosis += Form("%.1f",ptTriggerMin); filenameKurtosis += "To";
 792   filenameKurtosis += Form("%.1f",ptTriggerMax); filenameKurtosis += "PtaFrom";
 793   filenameKurtosis += Form("%.1f",ptAssociatedMin); filenameKurtosis += "To";
 794   filenameKurtosis += Form("%.1f",ptAssociatedMax);
 795   filenameKurtosis += ".txt";
 796
 797   TString meanLatex, rmsLatex, skewnessLatex, kurtosisLatex;
 798   meanLatex = "#mu = ";
 799   meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMean());
 800   meanLatex += " #pm ";
 801   meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMeanError());
 802   ofstream fileMean(filenameMean.Data(),ios::app);
 803   fileMean << gCentrality << " " << gHistBalanceFunctionSubtracted->GetMean() << " " <<gHistBalanceFunctionSubtracted->GetMeanError()<<endl;
 804   fileMean.close();
 805
 806   rmsLatex = "#sigma = ";
 807   rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMS());
 808   rmsLatex += " #pm ";
 809   rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMSError());
 810   ofstream fileSigma(filenameSigma.Data(),ios::app);
 811   fileSigma << gCentrality << " " << gHistBalanceFunctionSubtracted->GetRMS() << " " <<gHistBalanceFunctionSubtracted->GetRMSError()<<endl;
 812   fileSigma.close();
 813
 814   skewnessLatex = "S = ";
 815   skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(1));
 816   skewnessLatex += " #pm ";
 817   skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(11));
 818   ofstream fileSkewness(filenameSkewness.Data(),ios::app);
 819   fileSkewness << gCentrality << " " << gHistBalanceFunctionSubtracted->GetSkewness(1) << " " <<gHistBalanceFunctionSubtracted->GetSkewness(11)<<endl;
 820   fileSkewness.close();
 821
 822   kurtosisLatex = "K = ";
 823   kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(1));
 824   kurtosisLatex += " #pm ";
 825   kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(11));
 826   ofstream fileKurtosis(filenameKurtosis.Data(),ios::app);
 827   fileKurtosis << gCentrality << " " << gHistBalanceFunctionSubtracted->GetKurtosis(1) << " " <<gHistBalanceFunctionSubtracted->GetKurtosis(11)<<endl;
 828   fileKurtosis.close();
 829
 830   Printf("Mean: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetMean(),gHistBalanceFunctionSubtracted->GetMeanError());
 831   Printf("RMS: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetRMS(),gHistBalanceFunctionSubtracted->GetRMSError());
 832   Printf("Skeweness: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetSkewness(1),gHistBalanceFunctionSubtracted->GetSkewness(11));
 833   Printf("Kurtosis: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetKurtosis(1),gHistBalanceFunctionSubtracted->GetKurtosis(11));
 834
 835   latexInfo1->DrawLatex(0.18,0.88,centralityLatex.Data());
 836   latexInfo1->DrawLatex(0.18,0.82,psiLatex.Data());
 837   latexInfo1->DrawLatex(0.18,0.76,pttLatex.Data());
 838   latexInfo1->DrawLatex(0.18,0.70,ptaLatex.Data());
 839
 840   TLatex *latexResults = new TLatex();
 841   latexResults->SetNDC();
 842   latexResults->SetTextSize(0.04);
 843   latexResults->SetTextColor(1);
 844   latexResults->DrawLatex(0.6,0.88,meanLatex.Data());
 845   latexResults->DrawLatex(0.6,0.82,rmsLatex.Data());
 846   latexResults->DrawLatex(0.6,0.76,skewnessLatex.Data());
 847   latexResults->DrawLatex(0.6,0.70,kurtosisLatex.Data());
 848
 849 }