PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C

   1 const Int_t numberOfCentralityBins = 12;
   2 TString centralityArray[numberOfCentralityBins] = {"0-100","10-20","20-30","30-40","40-50","50-60","60-70","70-80","0-100","0-1","1-2","2-3"};
   3
   4
   5 const Int_t gRebin = 1;
   6 void drawBalanceFunction2DPsi(const char* filename = "AnalysisResultsPsi.root",
   7                               Int_t gCentrality = 1,
   8                               Int_t gBit = -1,
   9                               const char* gCentralityEstimator = 0x0,
  10                               Bool_t kShowShuffled = kFALSE,
  11                               Bool_t kShowMixed = kTRUE,
  12                               Double_t psiMin = -0.5, Double_t psiMax = 0.5,
  13                               Double_t vertexZMin = -10.,
  14                               Double_t vertexZMax = 10.,
  15                               Double_t ptTriggerMin = -1.,
  16                               Double_t ptTriggerMax = -1.,
  17                               Double_t ptAssociatedMin = -1.,
  18                               Double_t ptAssociatedMax = -1.,
  19                               Bool_t kUseVzBinning = kTRUE,
  20                               Bool_t k2pMethod = kTRUE,
  21                               TString eventClass = "EventPlane", //Can be "EventPlane", "Centrality", "Multiplicity"
  22                               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   //Macro that draws the balance functions PROJECTIONS
 938   //for each centrality bin for the different pT of trigger and
 939   //associated particles
 940   TGaxis::SetMaxDigits(3);
 941
 942   //first we need some libraries
 943   gSystem->Load("libTree");
 944   gSystem->Load("libGeom");
 945   gSystem->Load("libVMC");
 946   gSystem->Load("libXMLIO");
 947   gSystem->Load("libPhysics");
 948
 949   gSystem->Load("libSTEERBase");
 950   gSystem->Load("libESD");
 951   gSystem->Load("libAOD");
 952
 953   gSystem->Load("libANALYSIS.so");
 954   gSystem->Load("libANALYSISalice.so");
 955   gSystem->Load("libEventMixing.so");
 956   gSystem->Load("libCORRFW.so");
 957   gSystem->Load("libPWGTools.so");
 958   gSystem->Load("libPWGCFebye.so");
 959
 960   gStyle->SetOptStat(0);
 961
 962   //Get the input file
 963   TString filename = "balanceFunction2D.";
 964   if(eventClass == "Centrality"){
 965     filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
 966     filename += ".PsiAll.PttFrom";
 967   }
 968   else if(eventClass == "Multiplicity"){
 969     filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
 970     filename += ".PsiAll.PttFrom";
 971   }
 972   else{ // "EventPlane" (default)
 973     filename += "Centrality";
 974     filename += gCentrality; filename += ".Psi";
 975     if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
 976     else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
 977     else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
 978     else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
 979     else filename += "All.PttFrom";
 980   }
 981   filename += Form("%.1f",ptTriggerMin); filename += "To";
 982   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
 983   filename += Form("%.1f",ptAssociatedMin); filename += "To";
 984   filename += Form("%.1f",ptAssociatedMax);
 985   if(k2pMethod) filename += "_2pMethod";
 986
 987   filename += "_";
 988   filename += Form("%.1f",psiMin);
 989   filename += "-";
 990   filename += Form("%.1f",psiMax);
 991   filename += ".root";
 992
 993   //Open the file
 994   TFile *f = 0x0;
 995   if(!gHistBalanceFunction2D) {
 996     TFile::Open(filename.Data());
 997     if((!f)||(!f->IsOpen())) {
 998       Printf("The file %s is not found. Aborting...",filename);
 999       return listBF;
1000     }
1001     //f->ls();
1002   }
1003
1004   //Latex
1005   TString centralityLatex = "Centrality: ";
1006   centralityLatex += centralityArray[gCentrality-1];
1007   centralityLatex += "%";
1008
1009   TString psiLatex;
1010   if((psiMin == -0.5)&&(psiMax == 0.5))
1011     psiLatex = " -7.5^{o} < #varphi^{t} - #Psi_{2} < 7.5^{o}";
1012   else if((psiMin == 0.5)&&(psiMax == 1.5))
1013     psiLatex = " 37.5^{o} < #varphi^{t} - #Psi_{2} < 52.5^{o}";
1014   else if((psiMin == 1.5)&&(psiMax == 2.5))
1015     psiLatex = " 82.5^{o} < #varphi^{t} - #Psi_{2} < 97.5^{o}";
1016   else
1017     psiLatex = " 0^{o} < #varphi^{t} - #Psi_{2} < 180^{o}";
1018
1019   TString pttLatex = Form("%.1f",ptTriggerMin);
1020   pttLatex += " < p_{T,trig} < "; pttLatex += Form("%.1f",ptTriggerMax);
1021   pttLatex += " GeV/c";
1022
1023   TString ptaLatex = Form("%.1f",ptAssociatedMin);
1024   ptaLatex += " < p_{T,assoc} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1025   ptaLatex += " GeV/c";
1026
1027   TLatex *latexInfo1 = new TLatex();
1028   latexInfo1->SetNDC();
1029   latexInfo1->SetTextSize(0.045);
1030   latexInfo1->SetTextColor(1);
1031
1032
1033   //============================================================//
1034   //Get subtracted and mixed balance function
1035   TH2D *gHistBalanceFunctionSubtracted2D = 0x0;
1036   TH2D *gHistBalanceFunctionMixed2D      = 0x0;
1037   if(!gHistBalanceFunction2D) {
1038     gHistBalanceFunctionSubtracted2D = (TH2D*)f->Get("gHistBalanceFunctionSubtracted");
1039     gHistBalanceFunctionMixed2D      = (TH2D*)f->Get("gHistBalanceFunctionMixed");
1040   }
1041   else {
1042     gHistBalanceFunctionSubtracted2D = dynamic_cast<TH2D*>(gHistBalanceFunction2D->Clone());
1043     gHistBalanceFunctionMixed2D = dynamic_cast<TH2D*>(gHistBalanceFunction2D->Clone());
1044   }
1045
1046   TH1D *gHistBalanceFunctionSubtracted = NULL;
1047   TH1D *gHistBalanceFunctionMixed      = NULL;
1048
1049   TH1D *gHistBalanceFunctionSubtracted_scale = NULL;
1050   TH1D *gHistBalanceFunctionMixed_scale      = NULL;
1051
1052   if(kProjectInEta){
1053     gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionX("gHistBalanceFunctionSubtractedEta",binMin,binMax));
1054     gHistBalanceFunctionSubtracted->Scale(gHistBalanceFunctionSubtracted2D->GetYaxis()->GetBinWidth(1));   // to remove normalization to phi bin width
1055     gHistBalanceFunctionMixed      = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionX("gHistBalanceFunctionMixedEta",binMin,binMax));
1056     gHistBalanceFunctionMixed->Scale(gHistBalanceFunctionMixed2D->GetYaxis()->GetBinWidth(1));   // to remove normalization to phi bin width
1057     gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#eta)");
1058     gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#eta)");
1059   }
1060   else{
1061     gHistBalanceFunctionSubtracted = dynamic_cast<TH1D *>(gHistBalanceFunctionSubtracted2D->ProjectionY("gHistBalanceFunctionSubtractedPhi",binMin,binMax));
1062     gHistBalanceFunctionSubtracted->Scale(gHistBalanceFunctionSubtracted2D->GetXaxis()->GetBinWidth(1));   // to remove normalization to eta bin width
1063     gHistBalanceFunctionMixed      = dynamic_cast<TH1D *>(gHistBalanceFunctionMixed2D->ProjectionY("gHistBalanceFunctionMixedPhi",binMin,binMax));
1064     gHistBalanceFunctionMixed->Scale(gHistBalanceFunctionMixed2D->GetXaxis()->GetBinWidth(1));   // to remove normalization to eta bin width
1065     gHistBalanceFunctionSubtracted->SetTitle("B(#Delta#varphi)");
1066     gHistBalanceFunctionMixed->SetTitle("B_{mix}(#Delta#varphi)");
1067   }
1068
1069   gHistBalanceFunctionSubtracted->SetMarkerStyle(20);
1070   gHistBalanceFunctionSubtracted->GetYaxis()->SetTitleOffset(1.3);
1071
1072   gHistBalanceFunctionMixed->SetMarkerStyle(25);
1073
1074   TCanvas *c1 = new TCanvas("c1","",0,0,600,500);
1075   c1->SetFillColor(10);
1076   c1->SetHighLightColor(10);
1077   c1->SetLeftMargin(0.15);
1078   gHistBalanceFunctionSubtracted->DrawCopy("E");
1079   gHistBalanceFunctionMixed->DrawCopy("ESAME");
1080
1081   legend = new TLegend(0.18,0.62,0.45,0.82,"","brNDC");
1082   legend->SetTextSize(0.045);
1083   legend->SetTextFont(42);
1084   legend->SetBorderSize(0);
1085   legend->SetFillStyle(0);
1086   legend->SetFillColor(10);
1087   legend->SetMargin(0.25);
1088   legend->SetShadowColor(10);
1089   legend->AddEntry(gHistBalanceFunctionSubtracted,"Data","lp");
1090   legend->AddEntry(gHistBalanceFunctionMixed,"Mixed data","lp");
1091   legend->Draw();
1092
1093
1094   TString meanLatex, rmsLatex, skewnessLatex, kurtosisLatex;
1095
1096   if(bRootMoments){
1097
1098     // need to restrict to near side in case of dphi
1099     if(!kProjectInEta) gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-TMath::Pi()/2,TMath::Pi()/2);
1100
1101     if(bReduceRangeForMoments){
1102
1103       // default (for 2<pT,assoc<3<pT,trig<4)
1104       Double_t rangeReduced = 1.0;
1105
1106       //for 3<pT,assoc<8<pT,trig<15
1107       if(ptTriggerMax>10){
1108         rangeReduced = 0.3;
1109       }
1110
1111       cout<<"Use reduced range = "<<rangeReduced<<endl;
1112       gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-rangeReduced,rangeReduced);
1113     }
1114
1115     meanLatex = "#mu = ";
1116     meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMean());
1117     meanLatex += " #pm ";
1118     meanLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetMeanError());
1119
1120     rmsLatex = "#sigma = ";
1121     rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMS());
1122     rmsLatex += " #pm ";
1123     rmsLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetRMSError());
1124
1125     skewnessLatex = "S = ";
1126     skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(1));
1127     skewnessLatex += " #pm ";
1128     skewnessLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetSkewness(11));
1129
1130     kurtosisLatex = "K = ";
1131     kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(1));
1132     kurtosisLatex += " #pm ";
1133     kurtosisLatex += Form("%.3f",gHistBalanceFunctionSubtracted->GetKurtosis(11));
1134
1135     Printf("ROOT Mean: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetMean(),gHistBalanceFunctionSubtracted->GetMeanError());
1136     Printf("ROOT RMS: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetRMS(),gHistBalanceFunctionSubtracted->GetRMSError());
1137     Printf("ROOT Skeweness: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetSkewness(1),gHistBalanceFunctionSubtracted->GetSkewness(11));
1138     Printf("ROOT Kurtosis: %lf - Error: %lf",gHistBalanceFunctionSubtracted->GetKurtosis(1),gHistBalanceFunctionSubtracted->GetKurtosis(11));
1139
1140
1141     // store in txt files
1142     TString meanFileName = filename;
1143     if(kProjectInEta)
1144       meanFileName= "deltaEtaProjection_Mean.txt";
1145     //meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
1146     else
1147       meanFileName = "deltaPhiProjection_Mean.txt";
1148       //meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
1149     ofstream fileMean(meanFileName.Data(),ios::app);
1150     fileMean << " " << gHistBalanceFunctionSubtracted->GetMean() << " " <<gHistBalanceFunctionSubtracted->GetMeanError()<<endl;
1151     fileMean.close();
1152
1153     TString rmsFileName = filename;
1154     if(kProjectInEta)
1155       rmsFileName = "deltaEtaProjection_Rms.txt";
1156       //rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
1157     else
1158       rmsFileName = "deltaPhiProjection_Rms.txt";
1159       //rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
1160     ofstream fileRms(rmsFileName.Data(),ios::app);
1161     fileRms << " " << gHistBalanceFunctionSubtracted->GetRMS() << " " <<gHistBalanceFunctionSubtracted->GetRMSError()<<endl;
1162     fileRms.close();
1163
1164     TString skewnessFileName = filename;
1165     if(kProjectInEta)
1166       skewnessFileName = "deltaEtaProjection_Skewness.txt";
1167       //skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
1168     else
1169       skewnessFileName = "deltaPhiProjection_Skewness.txt";
1170     //skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
1171     ofstream fileSkewness(skewnessFileName.Data(),ios::app);
1172     fileSkewness << " " << gHistBalanceFunctionSubtracted->GetSkewness(1) << " " <<gHistBalanceFunctionSubtracted->GetSkewness(11)<<endl;
1173     fileSkewness.close();
1174
1175     TString kurtosisFileName = filename;
1176     if(kProjectInEta)
1177       kurtosisFileName = "deltaEtaProjection_Kurtosis.txt";
1178       //kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
1179     else
1180       kurtosisFileName = "deltaPhiProjection_Kurtosis.txt";
1181       //kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
1182     ofstream fileKurtosis(kurtosisFileName.Data(),ios::app);
1183     fileKurtosis << " " << gHistBalanceFunctionSubtracted->GetKurtosis(1) << " " <<gHistBalanceFunctionSubtracted->GetKurtosis(11)<<endl;
1184     fileKurtosis.close();
1185
1186     if(!kProjectInEta) gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-TMath::Pi()/2,3.*TMath::Pi()/2);
1187     else               gHistBalanceFunctionSubtracted->GetXaxis()->SetRangeUser(-1.6,1.6);
1188   }
1189   // calculate the moments by hand
1190   else{
1191
1192     Double_t meanAnalytical, meanAnalyticalError;
1193     Double_t sigmaAnalytical, sigmaAnalyticalError;
1194     Double_t skewnessAnalytical, skewnessAnalyticalError;
1195     Double_t kurtosisAnalytical, kurtosisAnalyticalError;
1196
1197     Int_t gDeltaEtaPhi = 2;
1198     if(kProjectInEta) gDeltaEtaPhi = 1;
1199
1200     AliBalancePsi *bHelper = new AliBalancePsi;
1201     bHelper->GetMomentsAnalytical(gDeltaEtaPhi,gHistBalanceFunctionSubtracted,kUseZYAM,meanAnalytical,meanAnalyticalError,sigmaAnalytical,sigmaAnalyticalError,skewnessAnalytical,skewnessAnalyticalError,kurtosisAnalytical,kurtosisAnalyticalError);
1202
1203     meanLatex = "#mu = ";
1204     meanLatex += Form("%.3f",meanAnalytical);
1205     meanLatex += " #pm ";
1206     meanLatex += Form("%.3f",meanAnalyticalError);
1207
1208     rmsLatex = "#sigma = ";
1209     rmsLatex += Form("%.3f",sigmaAnalytical);
1210     rmsLatex += " #pm ";
1211     rmsLatex += Form("%.3f",sigmaAnalyticalError);
1212
1213     skewnessLatex = "S = ";
1214     skewnessLatex += Form("%.3f",skewnessAnalytical);
1215     skewnessLatex += " #pm ";
1216     skewnessLatex += Form("%.3f",skewnessAnalyticalError);
1217
1218     kurtosisLatex = "K = ";
1219     kurtosisLatex += Form("%.3f",kurtosisAnalytical);
1220     kurtosisLatex += " #pm ";
1221     kurtosisLatex += Form("%.3f",kurtosisAnalyticalError);
1222     Printf("Mean: %lf - Error: %lf",meanAnalytical, meanAnalyticalError);
1223     Printf("Sigma: %lf - Error: %lf",sigmaAnalytical, sigmaAnalyticalError);
1224     Printf("Skeweness: %lf - Error: %lf",skewnessAnalytical, skewnessAnalyticalError);
1225     Printf("Kurtosis: %lf - Error: %lf",kurtosisAnalytical, kurtosisAnalyticalError);
1226
1227     // store in txt files
1228     TString meanFileName = filename;
1229     if(kProjectInEta)
1230       meanFileName = "deltaEtaProjection_Mean.txt";
1231       //meanFileName.ReplaceAll(".root","_DeltaEtaProjection_Mean.txt");
1232     else
1233       meanFileName = "deltaPhiProjection_Mean.txt";
1234       //meanFileName.ReplaceAll(".root","_DeltaPhiProjection_Mean.txt");
1235     ofstream fileMean(meanFileName.Data(),ios::app);
1236     fileMean << " " << meanAnalytical << " " <<meanAnalyticalError<<endl;
1237     fileMean.close();
1238
1239     TString rmsFileName = filename;
1240     if(kProjectInEta)
1241       rmsFileName = "deltaEtaProjection_Rms.txt";
1242       //rmsFileName.ReplaceAll(".root","_DeltaEtaProjection_Rms.txt");
1243     else
1244       rmsFileName = "deltaPhiProjection_Rms.txt";
1245 //rmsFileName.ReplaceAll(".root","_DeltaPhiProjection_Rms.txt");
1246     ofstream fileRms(rmsFileName.Data(),ios::app);
1247     fileRms << " " << sigmaAnalytical << " " <<sigmaAnalyticalError<<endl;
1248     fileRms.close();
1249
1250     TString skewnessFileName = filename;
1251     if(kProjectInEta)
1252       skewnessFileName = "deltaEtaProjection_Skewness.txt";
1253       //skewnessFileName.ReplaceAll(".root","_DeltaEtaProjection_Skewness.txt");
1254     else
1255       skewnessFileName = "deltaPhiProjection_Skewness.txt";
1256       //skewnessFileName.ReplaceAll(".root","_DeltaPhiProjection_Skewness.txt");
1257     ofstream fileSkewness(skewnessFileName.Data(),ios::app);
1258     fileSkewness << " " << skewnessAnalytical << " " <<skewnessAnalyticalError<<endl;
1259     fileSkewness.close();
1260
1261     TString kurtosisFileName = filename;
1262     if(kProjectInEta)
1263       kurtosisFileName = "deltaEtaProjection_Kurtosis.txt";
1264       //kurtosisFileName.ReplaceAll(".root","_DeltaEtaProjection_Kurtosis.txt");
1265     else
1266       kurtosisFileName = "deltaPhiProjection_Kurtosis.txt";
1267       //kurtosisFileName.ReplaceAll(".root","_DeltaPhiProjection_Kurtosis.txt");
1268     ofstream fileKurtosis(kurtosisFileName.Data(),ios::app);
1269     fileKurtosis << " " << kurtosisAnalytical << " " <<kurtosisAnalyticalError<<endl;
1270     fileKurtosis.close();
1271   }
1272
1273   // Weighted mean as calculated for 1D analysis
1274   Double_t weightedMean, weightedMeanError;
1275   if(bReduceRangeForMoments){
1276     GetWeightedMean1D(gHistBalanceFunctionSubtracted,kProjectInEta,1,-1,weightedMean,weightedMeanError,rangeReduced);
1277   }
1278   else{
1279     GetWeightedMean1D(gHistBalanceFunctionSubtracted,kProjectInEta,1,-1,weightedMean,weightedMeanError,-1.);
1280   }
1281   Printf("Weighted Mean: %lf - Error: %lf",weightedMean, weightedMeanError);
1282
1283   // store in txt files
1284   TString weightedMeanFileName = filename;
1285   if(kProjectInEta)
1286     weightedMeanFileName = "deltaEtaProjection_WeightedMean.txt";
1287     //weightedMeanFileName.ReplaceAll(".root","_DeltaEtaProjection_WeightedMean.txt");
1288   else
1289     weightedMeanFileName = "deltaPhiProjection_WeightedMean.txt";
1290     //weightedMeanFileName.ReplaceAll(".root","_DeltaPhiProjection_WeightedMean.txt");
1291   ofstream fileWeightedMean(weightedMeanFileName.Data(),ios::app);
1292   fileWeightedMean << " " << weightedMean << " " <<weightedMeanError<<endl;
1293   fileWeightedMean.close();
1294
1295
1296   TCanvas *c2 = new TCanvas("c2","",600,0,600,500);
1297   c2->SetFillColor(10);
1298   c2->SetHighLightColor(10);
1299   c2->SetLeftMargin(0.15);
1300   gHistBalanceFunctionSubtracted->DrawCopy("E");
1301
1302   TLatex *latex = new TLatex();
1303   latex->SetNDC();
1304   latex->SetTextSize(0.035);
1305   latex->SetTextColor(1);
1306   latex->DrawLatex(0.64,0.85,meanLatex.Data());
1307   latex->DrawLatex(0.64,0.81,rmsLatex.Data());
1308   latex->DrawLatex(0.64,0.77,skewnessLatex.Data());
1309   latex->DrawLatex(0.64,0.73,kurtosisLatex.Data());
1310
1311   TString pngName = filename;
1312   if(kProjectInEta) pngName.ReplaceAll(".root","_DeltaEtaProjection.png");
1313   else              pngName.ReplaceAll(".root","_DeltaPhiProjection.png");
1314
1315   c2->SaveAs(pngName.Data());
1316
1317   TString outFileName = filename;
1318   if(kProjectInEta) outFileName.ReplaceAll(".root","_DeltaEtaProjection.root");
1319   else              outFileName.ReplaceAll(".root","_DeltaPhiProjection.root");
1320   TFile *fProjection = TFile::Open(outFileName.Data(),"recreate");
1321   gHistBalanceFunctionSubtracted->Write();
1322   gHistBalanceFunctionMixed->Write();
1323   fProjection->Close();
1324 }
1325
1326 //____________________________________________________________//
1327 void drawBFPsi2DFromCorrelationFunctions(const char* lhcPeriod = "LHC10h",
1328                                          Int_t gTrainNumber = 64,
1329                                          const char* gCentralityEstimator = "V0M",
1330                                          Int_t gBit = 128,
1331                                          const char* gEventPlaneEstimator = "VZERO",
1332                                          Int_t gCentrality = 1,
1333                                          Double_t psiMin = -0.5, Double_t psiMax = 3.5,
1334                                          Double_t vertexZMin = -10.,
1335                                          Double_t vertexZMax = 10.,
1336                                          Double_t ptTriggerMin = -1.,
1337                                          Double_t ptTriggerMax = -1.,
1338                                          Double_t ptAssociatedMin = -1.,
1339                                          Double_t ptAssociatedMax = -1.,
1340                                          TString eventClass = "Multiplicity",
1341                                          Bool_t bRootMoments = kFALSE,
1342                                          Bool_t bFullPhiForEtaProjection = kFALSE,
1343                                          Bool_t bReduceRangeForMoments = kFALSE
1344 ) {
1345   //Macro that draws the BF distributions for each centrality bin
1346   //for reaction plane dependent analysis
1347   //Author: Panos.Christakoglou@nikhef.nl
1348   TGaxis::SetMaxDigits(3);
1349   gStyle->SetPalette(55,0);
1350
1351   //Get the input file
1352   TString filename = lhcPeriod;
1353   if(lhcPeriod != ""){
1354     //filename += "/Train"; filename += gTrainNumber;
1355     filename +="/PttFrom";
1356     filename += Form("%.1f",ptTriggerMin); filename += "To";
1357     filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1358     filename += Form("%.1f",ptAssociatedMin); filename += "To";
1359     filename += Form("%.1f",ptAssociatedMax);
1360     filename += "/correlationFunction.";
1361   }
1362   else{
1363     filename += "correlationFunction.";
1364   }
1365   if(eventClass == "Centrality"){
1366     filename += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1367     filename += ".PsiAll.PttFrom";
1368   }
1369   else if(eventClass == "Multiplicity"){
1370     filename += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1371     filename += ".PsiAll.PttFrom";
1372   }
1373   else{ // "EventPlane" (default)
1374     filename += "Centrality";
1375     filename += gCentrality; filename += ".Psi";
1376     if((psiMin == -0.5)&&(psiMax == 0.5)) filename += "InPlane.Ptt";
1377     else if((psiMin == 0.5)&&(psiMax == 1.5)) filename += "Intermediate.Ptt";
1378     else if((psiMin == 1.5)&&(psiMax == 2.5)) filename += "OutOfPlane.Ptt";
1379     else if((psiMin == 2.5)&&(psiMax == 3.5)) filename += "Rest.PttFrom";
1380     else filename += "All.PttFrom";
1381   }
1382   filename += Form("%.1f",ptTriggerMin); filename += "To";
1383   filename += Form("%.1f",ptTriggerMax); filename += "PtaFrom";
1384   filename += Form("%.1f",ptAssociatedMin); filename += "To";
1385   filename += Form("%.1f",ptAssociatedMax);
1386   filename += "_";
1387   filename += Form("%.1f",psiMin);
1388   filename += "-";
1389   filename += Form("%.1f",psiMax);
1390   filename += ".root";
1391
1392   //Open the file
1393   TFile *f = TFile::Open(filename.Data());
1394   if((!f)||(!f->IsOpen())) {
1395     Printf("The file %s is not found. Aborting...",filename);
1396     return listBF;
1397   }
1398   //f->ls();
1399
1400   TH2D *gHistPN = dynamic_cast<TH2D *>(f->Get("gHistPNCorrelationFunctions"));
1401   if(!gHistPN) return;
1402   TH2D *gHistNP = dynamic_cast<TH2D *>(f->Get("gHistNPCorrelationFunctions"));
1403   if(!gHistNP) return;
1404   TH2D *gHistPP = dynamic_cast<TH2D *>(f->Get("gHistPPCorrelationFunctions"));
1405   if(!gHistPP) return;
1406   TH2D *gHistNN = dynamic_cast<TH2D *>(f->Get("gHistNNCorrelationFunctions"));
1407   if(!gHistNN) return;
1408
1409   gHistPN->Sumw2();
1410   gHistPP->Sumw2();
1411   gHistPN->Add(gHistPP,-1);
1412   gHistNP->Sumw2();
1413   gHistNN->Sumw2();
1414   gHistNP->Add(gHistNN,-1);
1415   gHistPN->Add(gHistNP);
1416   //gHistPN->Scale(0.5);//not needed anymore, since pT,assoc < pT,trig in all pT bins
1417   TH2D *gHistBalanceFunction2D = dynamic_cast<TH2D *>(gHistPN->Clone());
1418   gHistBalanceFunction2D->SetStats(kFALSE);
1419   gHistBalanceFunction2D->GetXaxis()->SetTitle("#Delta#eta");
1420   gHistBalanceFunction2D->GetYaxis()->SetTitle("#Delta#varphi (rad)");
1421   gHistBalanceFunction2D->GetZaxis()->SetTitle("B(#Delta#eta,#Delta#varphi)");
1422
1423   //Draw the results
1424   TCanvas *c0 = new TCanvas("c0","Balance function 2D",0,0,600,500);
1425   c0->SetFillColor(10); c0->SetHighLightColor(10);
1426   c0->SetLeftMargin(0.17); c0->SetTopMargin(0.05);
1427   gHistBalanceFunction2D->SetTitle("");
1428   gHistBalanceFunction2D->GetZaxis()->SetTitleOffset(1.4);
1429   gHistBalanceFunction2D->GetZaxis()->SetNdivisions(10);
1430   gHistBalanceFunction2D->GetYaxis()->SetTitleOffset(1.4);
1431   gHistBalanceFunction2D->GetYaxis()->SetNdivisions(10);
1432   gHistBalanceFunction2D->GetXaxis()->SetRangeUser(-1.4,1.4);
1433   gHistBalanceFunction2D->GetXaxis()->SetNdivisions(10);
1434   gHistBalanceFunction2D->DrawCopy("lego2");
1435   gPad->SetTheta(30); // default is 30
1436   gPad->SetPhi(-60); // default is 30
1437   gPad->Update();
1438
1439   TString multLatex = Form("Multiplicity: %.1f - %.1f",psiMin,psiMax);
1440
1441   TString pttLatex = Form("%.1f",ptTriggerMin);
1442   pttLatex += " < p_{T}^{t} < "; pttLatex += Form("%.1f",ptTriggerMax);
1443   pttLatex += " GeV/c";
1444
1445   TString ptaLatex = Form("%.1f",ptAssociatedMin);
1446   ptaLatex += " < p_{T}^{a} < "; ptaLatex += Form("%.1f",ptAssociatedMax);
1447   ptaLatex += " GeV/c";
1448
1449   TLatex *latexInfo1 = new TLatex();
1450   latexInfo1->SetNDC();
1451   latexInfo1->SetTextSize(0.045);
1452   latexInfo1->SetTextColor(1);
1453   latexInfo1->DrawLatex(0.54,0.88,multLatex.Data());
1454   latexInfo1->DrawLatex(0.54,0.82,pttLatex.Data());
1455   latexInfo1->DrawLatex(0.54,0.76,ptaLatex.Data());
1456
1457   TString pngName = "BalanceFunction2D.";
1458   pngName += Form("%s: %.1f - %.1f",eventClass.Data(),psiMin,psiMax);
1459   pngName += ".PttFrom";
1460   pngName += Form("%.1f",ptTriggerMin); pngName += "To";
1461   pngName += Form("%.1f",ptTriggerMax); pngName += "PtaFrom";
1462   pngName += Form("%.1f",ptAssociatedMin); pngName += "To";
1463   pngName += Form("%.1f",ptAssociatedMax);
1464   pngName += ".png";
1465
1466   c0->SaveAs(pngName.Data());
1467
1468  // do the full range for the projection on eta (for cross checking with published results)
1469  if(bFullPhiForEtaProjection){
1470
1471    drawProjections(gHistBalanceFunction2D,
1472                    kTRUE,
1473                    1,72,
1474                    gCentrality,
1475                    psiMin,psiMax,
1476                    ptTriggerMin,ptTriggerMax,
1477                    ptAssociatedMin,ptAssociatedMax,
1478                    kTRUE,
1479                    eventClass.Data(),
1480                    bRootMoments,
1481                    kFALSE,
1482                    bReduceRangeForMoments);
1483  }
1484  else{
1485   drawProjections(gHistBalanceFunction2D,
1486                    kTRUE,
1487                    1,36,
1488                    gCentrality,
1489                    psiMin,psiMax,
1490                    ptTriggerMin,ptTriggerMax,
1491                    ptAssociatedMin,ptAssociatedMax,
1492                    kTRUE,
1493                    eventClass.Data(),
1494                   bRootMoments,
1495                   kFALSE,
1496                   bReduceRangeForMoments);
1497  }
1498
1499   drawProjections(gHistBalanceFunction2D,
1500                   kFALSE,
1501                   1,80,
1502                   gCentrality,
1503                   psiMin,psiMax,
1504                   ptTriggerMin,ptTriggerMax,
1505                   ptAssociatedMin,ptAssociatedMax,
1506                   kTRUE,
1507                   eventClass.Data(),
1508                   bRootMoments,
1509                   kFALSE,
1510                   bReduceRangeForMoments);
1511
1512   TString outFileName = filename;
1513   outFileName.ReplaceAll("correlationFunction","balanceFunction2D");
1514   gHistBalanceFunction2D->SetName("gHistBalanceFunctionSubtracted");
1515   TFile *fOut = TFile::Open(outFileName.Data(),"recreate");
1516   gHistBalanceFunction2D->Write();
1517   fOut->Close();
1518
1519 }
1520
1521
1522 //____________________________________________________________//
1523 void mergeBFPsi2D(TString momDirectory = "./",
1524                   TString directory1 = "",
1525                   TString directory2 = "",
1526                   TString directory3 = "",
1527                   TString directory4 = "",
1528                   Int_t gCentrality = 1,
1529                   Double_t psiMin = -0.5, Double_t psiMax = 3.5,
1530                   TString  eventClass = "Centrality",
1531                   Double_t ptTriggerMin = -1.,
1532                   Double_t ptTriggerMax = -1.,
1533                   Double_t ptAssociatedMin = -1.,
1534                   Double_t ptAssociatedMax = -1.,
1535                   Bool_t bReduceRangeForMoments = kFALSE
1536 ) {
1537   //Macro that draws the BF distributions for each centrality bin
1538   //for reaction plane dependent analysis
1539   //Author: Panos.Christakoglou@nikhef.nl
1540   TGaxis::SetMaxDigits(3);
1541   gStyle->SetPalette(55,0);
1542
1543   cout<<"REDUCE"<<bReduceRangeForMoments<<"  "<<endl;
1544
1545   const Int_t nMaxDirectories = 4; // maximum number of directories to merge (set to 4 for now)
1546   TString sDirectory[nMaxDirectories];
1547   Int_t nDirectories = nMaxDirectories;
1548
1549   TString sFileName[nMaxDirectories];
1550   TFile *fBF[nMaxDirectories];
1551   TH2D  *hBF[nMaxDirectories];
1552   Double_t entries[nMaxDirectories];
1553
1554   TFile *fOut;
1555   TH2D  *hBFOut;
1556   Double_t entriesOut = 0.;
1557
1558   // find out how many directories we have to merge
1559  if(directory1==""){
1560     nDirectories=0;
1561   }
1562   else if(directory2==""){
1563     nDirectories=1;
1564     sDirectory[0]=directory1;
1565   }
1566   else if(directory3==""){
1567     nDirectories=2;
1568     sDirectory[0]=directory1;
1569     sDirectory[1]=directory2;
1570   }
1571   else if(directory4==""){
1572     nDirectories=3;
1573     sDirectory[0]=directory1;
1574     sDirectory[1]=directory2;
1575     sDirectory[2]=directory3;
1576   }
1577   else{
1578     nDirectories=nMaxDirectories;
1579     sDirectory[0]=directory1;
1580     sDirectory[1]=directory2;
1581     sDirectory[2]=directory3;
1582     sDirectory[3]=directory4;
1583   }
1584
1585  for(Int_t iDir = 0; iDir < nDirectories; iDir++){
1586
1587    //Get the input file
1588    sFileName[iDir] = sDirectory[iDir];
1589    sFileName[iDir] += "/Centrality"; sFileName[iDir] += gCentrality;
1590    sFileName[iDir] += "/";
1591 sFileName[iDir] += momDirectory;
1592    sFileName[iDir] += "/balanceFunction2D.";
1593
1594   if(eventClass == "Centrality"){
1595     sFileName[iDir] += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1596     sFileName[iDir] += ".PsiAll.PttFrom";
1597   }
1598   else if(eventClass == "Multiplicity"){
1599     sFileName[iDir] += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1600     sFileName[iDir] += ".PsiAll.PttFrom";
1601   }
1602   else{ // "EventPlane" (default)
1603     sFileName[iDir] += "Centrality";
1604     sFileName[iDir] += gCentrality; sFileName[iDir] += ".Psi";
1605     if((psiMin == -0.5)&&(psiMax == 0.5)) sFileName[iDir] += "InPlane.Ptt";
1606     else if((psiMin == 0.5)&&(psiMax == 1.5)) sFileName[iDir] += "Intermediate.Ptt";
1607     else if((psiMin == 1.5)&&(psiMax == 2.5)) sFileName[iDir] += "OutOfPlane.Ptt";
1608     else if((psiMin == 2.5)&&(psiMax == 3.5)) sFileName[iDir] += "Rest.PttFrom";
1609     else sFileName[iDir] += "All.PttFrom";
1610   }
1611   sFileName[iDir] += Form("%.1f",ptTriggerMin); sFileName[iDir] += "To";
1612   sFileName[iDir] += Form("%.1f",ptTriggerMax); sFileName[iDir] += "PtaFrom";
1613   sFileName[iDir] += Form("%.1f",ptAssociatedMin); sFileName[iDir] += "To";
1614   sFileName[iDir] += Form("%.1f",ptAssociatedMax);
1615   sFileName[iDir] += "_";
1616   sFileName[iDir] += Form("%.1f",psiMin);
1617   sFileName[iDir] += "-";
1618   sFileName[iDir] += Form("%.1f",psiMax);
1619   sFileName[iDir] += ".root";
1620
1621   //Open the file
1622   fBF[iDir] = TFile::Open(sFileName[iDir].Data());
1623   if((!fBF[iDir])||(!fBF[iDir]->IsOpen())) {
1624     Printf("The file %s is not found. Not used...",sFileName[iDir]);
1625     continue;
1626   }
1627   //fBF[iDir]->ls();
1628
1629   hBF[iDir]     = (TH2D*)fBF[iDir]->Get("gHistBalanceFunctionSubtracted");
1630   if(!hBF[iDir]) continue;
1631   entries[iDir] = (Double_t) hBF[iDir]->GetEntries();
1632   entriesOut += entries[iDir];
1633   cout<<" BF histogram "<<iDir<<" has "<<entries[iDir] <<" entries."<<endl;
1634
1635   // scaling and adding (for average)
1636   hBF[iDir]->Scale(entries[iDir]);
1637   if(iDir==0) hBFOut = (TH2D*)hBF[iDir]->Clone("gHistBalanceFunctionSubtractedOut");
1638   else hBFOut->Add(hBF[iDir]);
1639
1640  }
1641
1642  // scaling with all entries
1643  hBFOut->Scale(1./entriesOut);
1644
1645  drawProjections(hBFOut,
1646                   kTRUE,
1647                   1,36,
1648                   gCentrality,
1649                   psiMin,psiMax,
1650                   ptTriggerMin,ptTriggerMax,
1651                   ptAssociatedMin,ptAssociatedMax,
1652                   kTRUE,
1653                   eventClass,
1654                   kTRUE,
1655                   kFALSE,
1656                   bReduceRangeForMoments);
1657
1658   drawProjections(hBFOut,
1659                   kFALSE,
1660                   1,80,
1661                   gCentrality,
1662                   psiMin,psiMax,
1663                   ptTriggerMin,ptTriggerMax,
1664                   ptAssociatedMin,ptAssociatedMax,
1665                   kTRUE,
1666                   eventClass.Data(),
1667                   kTRUE,
1668                   kFALSE,
1669                   bReduceRangeForMoments);
1670
1671  // output
1672  TString outFileName = "balanceFunction2D.";
1673
1674  if(eventClass == "Centrality"){
1675    outFileName += Form("Centrality%.1fTo%.1f",psiMin,psiMax);
1676    outFileName += ".PsiAll.PttFrom";
1677  }
1678  else if(eventClass == "Multiplicity"){
1679    outFileName += Form("Multiplicity%.0fTo%.0f",psiMin,psiMax);
1680    outFileName += ".PsiAll.PttFrom";
1681   }
1682  else{ // "EventPlane" (default)
1683    outFileName += "Centrality";
1684    outFileName += gCentrality; outFileName += ".Psi";
1685    if((psiMin == -0.5)&&(psiMax == 0.5)) outFileName += "InPlane.Ptt";
1686    else if((psiMin == 0.5)&&(psiMax == 1.5)) outFileName += "Intermediate.Ptt";
1687    else if((psiMin == 1.5)&&(psiMax == 2.5)) outFileName += "OutOfPlane.Ptt";
1688    else if((psiMin == 2.5)&&(psiMax == 3.5)) outFileName += "Rest.PttFrom";
1689    else outFileName += "All.PttFrom";
1690  }
1691  outFileName += Form("%.1f",ptTriggerMin); outFileName += "To";
1692  outFileName += Form("%.1f",ptTriggerMax); outFileName += "PtaFrom";
1693  outFileName += Form("%.1f",ptAssociatedMin); outFileName += "To";
1694  outFileName += Form("%.1f",ptAssociatedMax);
1695  outFileName += "_";
1696  outFileName += Form("%.1f",psiMin);
1697  outFileName += "-";
1698  outFileName += Form("%.1f",psiMax);
1699  outFileName += ".root";
1700
1701  hBFOut->SetName("gHistBalanceFunctionSubtracted");
1702  fBFOut = TFile::Open(outFileName.Data(),"recreate");
1703  hBFOut->Write();
1704  fBFOut->Close();
1705
1706 }
1707
1708 //____________________________________________________________________//
1709 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.) {
1710
1711   //Prints the calculated width of the BF and its error
1712   Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
1713   Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
1714   Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
1715   Double_t deltaBalP2 = 0.0, integral = 0.0;
1716   Double_t deltaErrorNew = 0.0;
1717
1718   //Retrieve this variables from Histogram
1719   Int_t fNumberOfBins = gHistBalance->GetNbinsX();
1720   if(fStopBin > -1)   fNumberOfBins = fStopBin;
1721   Double_t fP2Step    = gHistBalance->GetBinWidth(1); // assume equal binning!
1722   Double_t currentBinCenter = 0.;
1723
1724   for(Int_t i = fStartBin; i <= fNumberOfBins; i++) {
1725
1726     // in order to recover the |abs| in the 1D analysis
1727     currentBinCenter = gHistBalance->GetBinCenter(i);
1728     if(kProjectInEta){
1729       if(currentBinCenter < 0) currentBinCenter = -currentBinCenter;
1730     }
1731     else{
1732       if(currentBinCenter < 0) currentBinCenter = -currentBinCenter;
1733       if(currentBinCenter > TMath::Pi()) currentBinCenter = 2 * TMath::Pi() - currentBinCenter;
1734     }
1735
1736     if(rangeReduced > 0 && currentBinCenter > rangeReduced )
1737       continue;
1738
1739     gSumXi += currentBinCenter;
1740     gSumBi += gHistBalance->GetBinContent(i);
1741     gSumBiXi += gHistBalance->GetBinContent(i)*currentBinCenter;
1742     gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(currentBinCenter,2);
1743     gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(currentBinCenter,2);
1744     gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
1745     gSumXi2DeltaBi2 += TMath::Power(currentBinCenter,2) * TMath::Power(gHistBalance->GetBinError(i),2);
1746
1747     deltaBalP2 += fP2Step*TMath::Power(gHistBalance->GetBinError(i),2);
1748     integral += fP2Step*gHistBalance->GetBinContent(i);
1749   }
1750   for(Int_t i = fStartBin; i < fNumberOfBins; i++)
1751     deltaErrorNew += gHistBalance->GetBinError(i)*(currentBinCenter*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
1752
1753   Double_t integralError = TMath::Sqrt(deltaBalP2);
1754
1755   Double_t delta = gSumBiXi / gSumBi;
1756   Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
1757
1758   WM  = delta;
1759   WME = deltaError;
1760 }
1761