PWGCF/EBYE/macros/drawBalanceFunctionPsi.C

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