PWGCF/EBYE/macros/drawBalanceFunction2DPsi.C

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