ROOT logo
/**************************************************************************
 * Author: Panos Christakoglou.                                           *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

//-----------------------------------------------------------------
//           Balance Function class
//   This is the class to deal with the Balance Function analysis
//   Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
//-----------------------------------------------------------------


//ROOT
#include <Riostream.h>
#include <TMath.h>
#include <TAxis.h>
#include <TFile.h>
#include <TF1.h>
#include <TH2D.h>
#include <TLorentzVector.h>
#include <TObjArray.h>
#include <TGraphErrors.h>
#include <TString.h>

#include "AliVParticle.h"
#include "AliMCParticle.h"
#include "AliESDtrack.h"
#include "AliAODTrack.h"

#include "AliBalance.h"

using std::cout;
using std::cerr;
using std::endl;

ClassImp(AliBalance)

//____________________________________________________________________//
AliBalance::AliBalance() :
  TObject(), 
  fShuffle(kFALSE),
  fHBTcut(kFALSE),
  fConversionCut(kFALSE),
  fAnalysisLevel("ESD"),
  fAnalyzedEvents(0) ,
  fCentralityId(0) ,
  fCentStart(0.),
  fCentStop(0.),
  fHistHBTbefore(NULL),
  fHistHBTafter(NULL),
  fHistConversionbefore(NULL),
  fHistConversionafter(NULL)
{
  // Default constructor
 
  for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
    if(i == 6) {
      fNumberOfBins[i] = 180;
      fP1Start[i]      = -360.0;
      fP1Stop[i]       = 360.0;
      fP2Start[i]      = -360.0;
      fP2Stop[i]       = 360.0;
      fP2Step[i]       = 0.1;
    }
    else {
      fNumberOfBins[i] = 20;
      fP1Start[i]      = -1.0;
      fP1Stop[i]       = 1.0;
      fP2Start[i]      = 0.0;
      fP2Stop[i]       = 2.0;
    }
    fP2Step[i] = TMath::Abs(fP2Start - fP2Stop) / (Double_t)fNumberOfBins[i];
    fCentStart = 0.;
    fCentStop  = 0.;

    fNn[i] = 0.0;
    fNp[i] = 0.0;

    for(Int_t j = 0; j < MAXIMUM_NUMBER_OF_STEPS; j++) {
      fNpp[i][j] = .0;
      fNnn[i][j] = .0;
      fNpn[i][j] = .0;
      fNnp[i][j] = .0;
      fB[i][j] = 0.0;
      ferror[i][j] = 0.0;
    }

    fHistP[i]  = NULL;
    fHistN[i]  = NULL;
    fHistPP[i] = NULL;
    fHistPN[i] = NULL;
    fHistNP[i] = NULL;
    fHistNN[i] = NULL;

  }
}


//____________________________________________________________________//
AliBalance::AliBalance(const AliBalance& balance):
  TObject(balance), 
  fShuffle(balance.fShuffle),
  fHBTcut(balance.fHBTcut), 
  fConversionCut(balance.fConversionCut), 
  fAnalysisLevel(balance.fAnalysisLevel),
  fAnalyzedEvents(balance.fAnalyzedEvents), 
  fCentralityId(balance.fCentralityId),
  fCentStart(balance.fCentStart),
  fCentStop(balance.fCentStop),
  fHistHBTbefore(balance.fHistHBTbefore),
  fHistHBTafter(balance.fHistHBTafter),
  fHistConversionbefore(balance.fHistConversionbefore),
  fHistConversionafter(balance.fHistConversionafter) {
  //copy constructor
  for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
    fNn[i] = balance.fNn[i];
    fNp[i] = balance.fNp[i];

    fP1Start[i]      = balance.fP1Start[i];
    fP1Stop[i]       = balance.fP1Stop[i];
    fNumberOfBins[i] = balance.fNumberOfBins[i];
    fP2Start[i]      = balance.fP2Start[i];
    fP2Stop[i]       = balance.fP2Stop[i];
    fP2Step[i]       = balance.fP2Step[i];
    fCentStart       = balance.fCentStart;
    fCentStop        = balance.fCentStop; 

    fHistP[i]        = balance.fHistP[i];
    fHistN[i]        = balance.fHistN[i];
    fHistPN[i]        = balance.fHistPN[i];
    fHistNP[i]        = balance.fHistNP[i];
    fHistPP[i]        = balance.fHistPP[i];
    fHistNN[i]        = balance.fHistNN[i];

    for(Int_t j = 0; j < MAXIMUM_NUMBER_OF_STEPS; j++) {
      fNpp[i][j] = .0;
      fNnn[i][j] = .0;
      fNpn[i][j] = .0;
      fNnp[i][j] = .0;
      fB[i][j] = 0.0;
      ferror[i][j] = 0.0;
    } 
  }
 }
 

//____________________________________________________________________//
AliBalance::~AliBalance() {
  // Destructor

  for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
 
    delete fHistP[i];
    delete fHistN[i];
    delete fHistPN[i];
    delete fHistNP[i];
    delete fHistPP[i];
    delete fHistNN[i];
  
  }
}

//____________________________________________________________________//
void AliBalance::SetInterval(Int_t iAnalysisType,
			     Double_t p1Start, Double_t p1Stop,
			     Int_t ibins, Double_t p2Start, Double_t p2Stop) {
  // Sets the analyzed interval. 
  // Set the same Information for all analyses

  if(iAnalysisType == -1){             
    for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
      fP1Start[i] = p1Start;
      fP1Stop[i] = p1Stop;
      fNumberOfBins[i] = ibins;
      fP2Start[i] = p2Start;
      fP2Stop[i] = p2Stop;
      fP2Step[i] = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins[i];
    }
  }
  // Set the Information for one analysis
  else if((iAnalysisType > -1) && (iAnalysisType < ANALYSIS_TYPES)) {
    fP1Start[iAnalysisType] = p1Start;
    fP1Stop[iAnalysisType] = p1Stop;
    fNumberOfBins[iAnalysisType] = ibins;
    fP2Start[iAnalysisType] = p2Start;
    fP2Stop[iAnalysisType] = p2Stop;
    fP2Step[iAnalysisType] = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins[iAnalysisType];
  }
  else {
    AliError("Wrong ANALYSIS number!");
  }
}

//____________________________________________________________________//
void AliBalance::InitHistograms() {
  //Initialize the histograms

  // global switch disabling the reference 
  // (to avoid "Replacing existing TH1" if several wagons are created in train)
  Bool_t oldStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE);

  TString histName;
  for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
    histName = "fHistP"; histName += kBFAnalysisType[iAnalysisType]; 
    if(fShuffle) histName.Append("_shuffle");
    if(fCentralityId) histName += fCentralityId.Data();
    fHistP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);

    histName = "fHistN"; histName += kBFAnalysisType[iAnalysisType]; 
    if(fShuffle) histName.Append("_shuffle");
    if(fCentralityId) histName += fCentralityId.Data();
    fHistN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
  
    histName = "fHistPN"; histName += kBFAnalysisType[iAnalysisType]; 
    if(fShuffle) histName.Append("_shuffle");
    if(fCentralityId) histName += fCentralityId.Data();
    fHistPN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
    
    histName = "fHistNP"; histName += kBFAnalysisType[iAnalysisType]; 
    if(fShuffle) histName.Append("_shuffle");
    if(fCentralityId) histName += fCentralityId.Data();
    fHistNP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
    
    histName = "fHistPP"; histName += kBFAnalysisType[iAnalysisType]; 
    if(fShuffle) histName.Append("_shuffle");
    if(fCentralityId) histName += fCentralityId.Data();
    fHistPP[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
    
    histName = "fHistNN"; histName += kBFAnalysisType[iAnalysisType]; 
    if(fShuffle) histName.Append("_shuffle");
    if(fCentralityId) histName += fCentralityId.Data();
    fHistNN[iAnalysisType] = new TH2D(histName.Data(),"",fCentStop-fCentStart,fCentStart,fCentStop,fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
  }

  // QA histograms
  fHistHBTbefore        = new TH2D("fHistHBTbefore","before HBT cut",200,0,2,200,0,200);
  fHistHBTafter         = new TH2D("fHistHBTafter","after HBT cut",200,0,2,200,0,200);
  fHistConversionbefore = new TH2D("fHistConversionbefore","before Conversion cut",200,0,2,200,0,200);
  fHistConversionafter  = new TH2D("fHistConversionafter","after Conversion cut",200,0,2,200,0,200);

  TH1::AddDirectory(oldStatus);

}

//____________________________________________________________________//
void AliBalance::PrintAnalysisSettings() {
  //prints the analysis settings
  
  Printf("======================================");
  Printf("Analysis level: %s",fAnalysisLevel.Data());
  Printf("======================================");
  for(Int_t ibin = 0; ibin < ANALYSIS_TYPES; ibin++){
    Printf("Interval info for variable %d",ibin);
    Printf("Analyzed interval (min.): %lf",fP2Start[ibin]);
    Printf("Analyzed interval (max.): %lf",fP2Stop[ibin]);
    Printf("Number of bins: %d",fNumberOfBins[ibin]);
    Printf("Step: %lf",fP2Step[ibin]);
    Printf("          ");
  }
  Printf("======================================");
}

//____________________________________________________________________//
void AliBalance::CalculateBalance(Float_t fCentrality,vector<Double_t> **chargeVector,Float_t bSign) {
  // Calculates the balance function
  fAnalyzedEvents++;
  Int_t i = 0 , j = 0;
  Int_t iBin = 0;

  // Initialize histograms if not done yet
  if(!fHistPN[0]){
    AliWarning("Histograms not yet initialized! --> Will be done now");
    AliWarning("This works only in local mode --> Add 'gBalance->InitHistograms()' in your configBalanceFunction");
    InitHistograms();
  }

  Int_t gNtrack = chargeVector[0]->size();
  //Printf("(AliBalance) Number of tracks: %d",gNtrack);

  for(i = 0; i < gNtrack;i++){
      Short_t charge          = chargeVector[0]->at(i);
      Double_t rapidity       = chargeVector[1]->at(i);
      Double_t pseudorapidity = chargeVector[2]->at(i);
      Double_t phi            = chargeVector[3]->at(i);
      
      //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
      for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
	if(iAnalysisType == kEta) {
	  if((pseudorapidity >= fP1Start[iAnalysisType]) && (pseudorapidity <= fP1Stop[iAnalysisType])) {
	    if(charge > 0) {
	      fNp[iAnalysisType] += 1.;
	      fHistP[iAnalysisType]->Fill(fCentrality,pseudorapidity);
	    }//charge > 0
	    if(charge < 0) {
	      fNn[iAnalysisType] += 1.;
	      fHistN[iAnalysisType]->Fill(fCentrality,pseudorapidity);
	    }//charge < 0
	  }//p1 interval check
	}//analysis type: eta
	else if(iAnalysisType == kPhi) {
	  if((phi >= fP1Start[iAnalysisType]) && (phi <= fP1Stop[iAnalysisType])) {
	    if(charge > 0) {
	      fNp[iAnalysisType] += 1.;
	      fHistP[iAnalysisType]->Fill(fCentrality,phi);
	    }//charge > 0
	    if(charge < 0) {
	      fNn[iAnalysisType] += 1.;
	      fHistN[iAnalysisType]->Fill(fCentrality,phi);
	    }//charge < 0
	  }//p1 interval check
	}//analysis type: phi
	else {
	  if((rapidity >= fP1Start[iAnalysisType]) && (rapidity <= fP1Stop[iAnalysisType])) {
	    if(charge > 0) {
	      fNp[iAnalysisType] += 1.;
	      fHistP[iAnalysisType]->Fill(fCentrality,rapidity);
	    }//charge > 0
	    if(charge < 0) {
	      fNn[iAnalysisType] += 1.;
	      fHistN[iAnalysisType]->Fill(fCentrality,rapidity);
	    }//charge < 0
	  }//p1 interval check
	}//analysis type: y, qside, qout, qlong, qinv
      }//analysis type loop
  }

  //Printf("Np: %lf - Nn: %lf",fNp[0],fNn[0]);

  Double_t dy = 0., deta = 0.;
  Double_t qLong = 0., qOut = 0., qSide = 0., qInv = 0.;
  Double_t dphi = 0.;

  Short_t charge1  = 0;
  Double_t eta1 = 0., rap1 = 0.;
  Double_t px1 = 0., py1 = 0., pz1 = 0.;
  Double_t pt1 = 0.;
  Double_t energy1 = 0.;
  Double_t phi1    = 0.;

  Short_t charge2  = 0;
  Double_t eta2 = 0., rap2 = 0.;
  Double_t px2 = 0., py2 = 0., pz2 = 0.;
  Double_t pt2 = 0.;
  Double_t energy2 = 0.;
  Double_t phi2    = 0.;
  //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
  for(i = 1; i < gNtrack; i++) {

      charge1 = chargeVector[0]->at(i);
      rap1    = chargeVector[1]->at(i);
      eta1    = chargeVector[2]->at(i);
      phi1    = chargeVector[3]->at(i);
      px1     = chargeVector[4]->at(i);
      py1     = chargeVector[5]->at(i);
      pz1     = chargeVector[6]->at(i);
      pt1     = chargeVector[7]->at(i);
      energy1 = chargeVector[8]->at(i);
    
    for(j = 0; j < i; j++) {
   
      charge2 = chargeVector[0]->at(j);
      rap2    = chargeVector[1]->at(j);
      eta2    = chargeVector[2]->at(j);
      phi2    = chargeVector[3]->at(j);
      px2     = chargeVector[4]->at(j);
      py2     = chargeVector[5]->at(j);
      pz2     = chargeVector[6]->at(j);
      pt2     = chargeVector[7]->at(j);
      energy2 = chargeVector[8]->at(j);

	// filling the arrays

	// RAPIDITY 
        dy = TMath::Abs(rap1 - rap2);

	// Eta
	deta = TMath::Abs(eta1 - eta2);

	//qlong
	Double_t eTot = energy1 + energy2;
	Double_t pxTot = px1 + px2;
	Double_t pyTot = py1 + py2;
	Double_t pzTot = pz1 + pz2;
	Double_t q0Tot = energy1 - energy2;
	Double_t qxTot = px1 - px2;
	Double_t qyTot = py1 - py2;
	Double_t qzTot = pz1 - pz2;

	Double_t eTot2 = eTot*eTot;
	Double_t pTot2 = pxTot*pxTot + pyTot*pyTot + pzTot*pzTot;
	Double_t pzTot2 = pzTot*pzTot;

	Double_t q0Tot2 = q0Tot*q0Tot;
	Double_t qTot2  = qxTot*qxTot + qyTot*qyTot + qzTot*qzTot;

	Double_t snn    = eTot2 - pTot2;
	Double_t ptTot2 = pTot2 - pzTot2 ;
	Double_t ptTot  = TMath::Sqrt( ptTot2 );
	
	qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/TMath::Sqrt(snn + ptTot2);
	
	//qout
	qOut = TMath::Sqrt(snn/(snn + ptTot2)) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
	
	//qside
	qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
	
	//qinv
	qInv = TMath::Sqrt(TMath::Abs(-q0Tot2 + qTot2 ));
	
	//phi
	dphi = TMath::Abs(phi1 - phi2);
	if(dphi>180) dphi = 360 - dphi;  //dphi should be between 0 and 180!

	// HBT like cut
	if(fHBTcut && charge1 * charge2 > 0){
	  //if( dphi < 3 || deta < 0.01 ){   // VERSION 1
	  //  continue;
	  
	  // VERSION 2 (Taken from DPhiCorrelations)
	  // the variables & cuthave been developed by the HBT group 
	  // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700

	  fHistHBTbefore->Fill(deta,dphi);
	  
	  // optimization
	  if (TMath::Abs(deta) < 0.02 * 2.5 * 3) //twoTrackEfficiencyCutValue = 0.02 [default for dphicorrelations]
	    {

	      // phi in rad
	      Float_t phi1rad = phi1*TMath::DegToRad();
	      Float_t phi2rad = phi2*TMath::DegToRad();

	      // check first boundaries to see if is worth to loop and find the minimum
	      Float_t dphistar1 = GetDPhiStar(phi1rad, pt1, charge1, phi2rad, pt2, charge2, 0.8, bSign);
	      Float_t dphistar2 = GetDPhiStar(phi1rad, pt1, charge1, phi2rad, pt2, charge2, 2.5, bSign);
	      
	      const Float_t kLimit = 0.02 * 3;
	      
	      Float_t dphistarminabs = 1e5;

	      if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0 )
		{
		  for (Double_t rad=0.8; rad<2.51; rad+=0.01) 
		    {
		      Float_t dphistar = GetDPhiStar(phi1rad, pt1, charge1, phi2rad, pt2, charge2, rad, bSign);
		      Float_t dphistarabs = TMath::Abs(dphistar);
		      
		      if (dphistarabs < dphistarminabs)
			{
			  dphistarminabs = dphistarabs;
			}
		    }
		
		  if (dphistarminabs < 0.02 && TMath::Abs(deta) < 0.02)
		    {
		      //AliInfo(Form("HBT: Removed track pair %d %d with [[%f %f]] %f %f %f | %f %f %d %f %f %d %f", i, j, deta, dphi, dphistarminabs, dphistar1, dphistar2, phi1rad, pt1, charge1, phi2rad, pt2, charge2, bSign));
		      continue;
		    }
		}
	    }
	  fHistHBTafter->Fill(deta,dphi);
	}
	
	// conversions
	if(fConversionCut){
	  if (charge1 * charge2 < 0)
	    {

	      fHistConversionbefore->Fill(deta,dphi);

	      Float_t m0 = 0.510e-3;
	      Float_t tantheta1 = 1e10;

	      // phi in rad
	      Float_t phi1rad = phi1*TMath::DegToRad();
	      Float_t phi2rad = phi2*TMath::DegToRad();
	      
	      if (eta1 < -1e-10 || eta1 > 1e-10)
		tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
	      
	      Float_t tantheta2 = 1e10;
	      if (eta2 < -1e-10 || eta2 > 1e-10)
		tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
	      
	      Float_t e1squ = m0 * m0 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
	      Float_t e2squ = m0 * m0 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);

	      Float_t masssqu = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( TMath::Cos(phi1rad - phi2rad) + 1.0 / tantheta1 / tantheta2 ) ) );

	      if (masssqu < 0.04*0.04){
		//AliInfo(Form("Conversion: Removed track pair %d %d with [[%f %f] %f %f] %d %d <- %f %f  %f %f   %f %f ", i, j, deta, dphi, masssqu, charge1, charge2,eta1,eta2,phi1,phi2,pt1,pt2));
		continue;
	      }
	      fHistConversionafter->Fill(deta,dphi);
	    }
	}
	
	
	//0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
	if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {

	  // rapidity
	  if( dy > fP2Start[kRapidity] && dy < fP2Stop[kRapidity]){
	    iBin = Int_t((dy-fP2Start[kRapidity])/fP2Step[kRapidity]);
	    if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
	      
	      if((charge1 > 0)&&(charge2 > 0)) {
		fNpp[kRapidity][iBin] += 1.;
		fHistPP[kRapidity]->Fill(fCentrality,dy);
	      }
	      else if((charge1 < 0)&&(charge2 < 0)) {
		fNnn[kRapidity][iBin] += 1.;
		fHistNN[kRapidity]->Fill(fCentrality,dy);
	      }
	      else if((charge1 > 0)&&(charge2 < 0)) {
		fNpn[kRapidity][iBin] += 1.;
		fHistPN[kRapidity]->Fill(fCentrality,dy);
	      }
	      else if((charge1 < 0)&&(charge2 > 0)) {
		fNpn[kRapidity][iBin] += 1.;
		    fHistPN[kRapidity]->Fill(fCentrality,dy);
	      }
	    }//BF binning check
	  }//p2 interval check
	}//p1 interval check
	
	// pseudorapidity
	if((eta1 >= fP1Start[kEta]) && (eta1 <= fP1Stop[kEta]) && (eta2 >= fP1Start[kEta]) && (eta2 <= fP1Stop[kEta])) {
	  if( deta > fP2Start[kEta] && deta < fP2Stop[kEta]){
	    iBin = Int_t((deta-fP2Start[kEta])/fP2Step[kEta]);	
	    if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
	      if((charge1 > 0)&&(charge2 > 0)) {
		fNpp[kEta][iBin] += 1.;
		fHistPP[kEta]->Fill(fCentrality,deta);
	      }
	      if((charge1 < 0)&&(charge2 < 0)) {
		fNnn[kEta][iBin] += 1.;
	   	    fHistNN[kEta]->Fill(fCentrality,deta);
	      }
	      if((charge1 > 0)&&(charge2 < 0)) {
		fNpn[kEta][iBin] += 1.;
		fHistPN[kEta]->Fill(fCentrality,deta);
	      }
	      if((charge1 < 0)&&(charge2 > 0)) {
		fNpn[kEta][iBin] += 1.;
	   	    fHistPN[kEta]->Fill(fCentrality,deta);
	      }
	    }//BF binning check
	  }//p2 interval check
	}//p1 interval check
	
	// Qlong, out, side, inv
	// Check the p1 intervall for rapidity here (like for single tracks above)
	if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
	  if( qLong > fP2Start[kQlong] && qLong < fP2Stop[kQlong]){
	    iBin = Int_t((qLong-fP2Start[kQlong])/fP2Step[kQlong]);	
	    if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
	      if((charge1 > 0)&&(charge2 > 0)) {
		fNpp[kQlong][iBin] += 1.;
		fHistPP[kQlong]->Fill(fCentrality,qLong);
	      }
	      if((charge1 < 0)&&(charge2 < 0)) {
		fNnn[kQlong][iBin] += 1.;
		fHistNN[kQlong]->Fill(fCentrality,qLong);
	      }
	      if((charge1 > 0)&&(charge2 < 0)) {
		fNpn[kQlong][iBin] += 1.;
		fHistPN[kQlong]->Fill(fCentrality,qLong);
	      }
	      if((charge1 < 0)&&(charge2 > 0)) {
		fNpn[kQlong][iBin] += 1.;
		fHistPN[kQlong]->Fill(fCentrality,qLong);
	      }
	    }//BF binning check
	  }//p2 interval check
	}//p1 interval check
	  
	if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
	  if( qOut > fP2Start[kQout] && qOut < fP2Stop[kQout]){
	    iBin = Int_t((qOut-fP2Start[kQout])/fP2Step[kQout]);	
	    if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
	      if((charge1 > 0)&&(charge2 > 0)) {
		fNpp[kQout][iBin] += 1.;
		fHistPP[kQout]->Fill(fCentrality,qOut);
	  	  }
	      if((charge1 < 0)&&(charge2 < 0)) {
		fNnn[kQout][iBin] += 1.;
		fHistNN[kQout]->Fill(fCentrality,qOut);
	      }
	      if((charge1 > 0)&&(charge2 < 0)) {
		fNpn[kQout][iBin] += 1.;
		fHistPN[kQout]->Fill(fCentrality,qOut);
	      }
	      if((charge1 < 0)&&(charge2 > 0)) {
		fNpn[kQout][iBin] += 1.;
		fHistPN[kQout]->Fill(fCentrality,qOut);
	      }
	    }//BF binning check
	  }//p2 interval check
	}//p1 interval check	
	
	if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
	  if( qSide > fP2Start[kQside] && qSide < fP2Stop[kQside]){
	    iBin = Int_t((qSide-fP2Start[kQside])/fP2Step[kQside]);	
	    if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
	      if((charge1 > 0)&&(charge2 > 0)) {
		fNpp[kQside][iBin] += 1.;
		fHistPP[kQside]->Fill(fCentrality,qSide);
	      }
	      if((charge1 < 0)&&(charge2 < 0)) {
		fNnn[kQside][iBin] += 1.;
		fHistNN[kQside]->Fill(fCentrality,qSide);
	      }
	      if((charge1 > 0)&&(charge2 < 0)) {
		fNpn[kQside][iBin] += 1.;
		fHistPN[kQside]->Fill(fCentrality,qSide);
	      }
	      if((charge1 < 0)&&(charge2 > 0)) {
		fNpn[kQside][iBin] += 1.;
		fHistPN[kQside]->Fill(fCentrality,qSide);
		  	  }
	    }//BF binning check
	  }//p2 interval check
	}//p1 interval check
	
	if((rap1 >= fP1Start[kRapidity]) && (rap1 <= fP1Stop[kRapidity]) && (rap2 >= fP1Start[kRapidity]) && (rap2 <= fP1Stop[kRapidity])) {
	  if( qInv > fP2Start[kQinv] && qInv < fP2Stop[kQinv]){
	    iBin = Int_t((qInv-fP2Start[kQinv])/fP2Step[kQinv]);	
	    if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
	      if((charge1 > 0)&&(charge2 > 0)) {
		fNpp[kQinv][iBin] += 1.;
		fHistPP[kQinv]->Fill(fCentrality,qInv);
	      }
	      if((charge1 < 0)&&(charge2 < 0)) {
		fNnn[kQinv][iBin] += 1.;
		fHistNN[kQinv]->Fill(fCentrality,qInv);
	      }
	      if((charge1 > 0)&&(charge2 < 0)) {
		fNpn[kQinv][iBin] += 1.;
		fHistPN[kQinv]->Fill(fCentrality,qInv);
	      }
	      if((charge1 < 0)&&(charge2 > 0)) {
		fNpn[kQinv][iBin] += 1.;
		fHistPN[kQinv]->Fill(fCentrality,qInv);
	      }
	    }//BF binning check
	  }//p2 interval check
	}//p1 interval check
	
	// Phi
	if((phi1 >= fP1Start[kPhi]) && (phi1 <= fP1Stop[kPhi]) && (phi2 >= fP1Start[kPhi]) && (phi2 <= fP1Stop[kPhi])) {
	  if( dphi > fP2Start[kPhi] && dphi < fP2Stop[kPhi]){
	    iBin = Int_t((dphi-fP2Start[kPhi])/fP2Step[kPhi]);	
	    if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
	      if((charge1 > 0)&&(charge2 > 0)) {
		fNpp[kPhi][iBin] += 1.;
		fHistPP[kPhi]->Fill(fCentrality,dphi);
	      }
	      if((charge1 < 0)&&(charge2 < 0)) {
		fNnn[kPhi][iBin] += 1.;
		fHistNN[kPhi]->Fill(fCentrality,dphi);
	      }
	      if((charge1 > 0)&&(charge2 < 0)) {
		fNpn[kPhi][iBin] += 1.;
		fHistPN[kPhi]->Fill(fCentrality,dphi);
	      }
	      if((charge1 < 0)&&(charge2 > 0)) {
		fNpn[kPhi][iBin] += 1.;
		fHistPN[kPhi]->Fill(fCentrality,dphi);
	      }
	    }//BF binning check
	  }//p2 interval check
	}//p1 interval check
    }//end of 2nd particle loop
  }//end of 1st particle loop
  //Printf("Number of analyzed events: %i",fAnalyzedEvents);
  //Printf("DeltaEta NN[0] = %.0f, PP[0] = %.0f, NP[0] = %.0f, PN[0] = %.0f",fNnn[kEta][0],fNpp[kEta][0],fNnp[kEta][0],fNpn[kEta][0]);
}  


//____________________________________________________________________//
Double_t AliBalance::GetBalance(Int_t iAnalysisType, Int_t p2) {
  // Returns the value of the balance function in bin p2
  fB[iAnalysisType][p2] = 0.5*(((fNpn[iAnalysisType][p2] - 2.*fNnn[iAnalysisType][p2])/fNn[iAnalysisType]) + ((fNpn[iAnalysisType][p2] - 2.*fNpp[iAnalysisType][p2])/fNp[iAnalysisType]))/fP2Step[iAnalysisType];
  
  return fB[iAnalysisType][p2];
}
    
//____________________________________________________________________//
Double_t AliBalance::GetError(Int_t iAnalysisType, Int_t p2) {        
  // Returns the error on the BF value for bin p2
  // The errors for fNn and fNp are neglected here (0.1 % of total error)
  /*ferror[iAnalysisType][p2] = TMath::Sqrt(Double_t(fNpp[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType]))
			      + Double_t(fNnn[iAnalysisType][p2])/(Double_t(fNn[iAnalysisType])*Double_t(fNn[iAnalysisType]))
			      + Double_t(fNpn[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType])) 
			      + Double_t(fNnp[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType]))
			      //+ TMath::Power(fNpn[iAnalysisType][p2]-fNpp[iAnalysisType][p2],2)/TMath::Power(Double_t(fNp[iAnalysisType]),3)
			      //+ TMath::Power(fNnp[iAnalysisType][p2]-fNnn[iAnalysisType][p2],2)/TMath::Power(Double_t(fNn[iAnalysisType]),3) 
			       ) /fP2Step[iAnalysisType];*/

  ferror[iAnalysisType][p2] = TMath::Sqrt( Double_t(fNpp[iAnalysisType][p2])/(Double_t(fNp[iAnalysisType])*Double_t(fNp[iAnalysisType])) + 
					   Double_t(fNnn[iAnalysisType][p2])/(Double_t(fNn[iAnalysisType])*Double_t(fNn[iAnalysisType])) + 
					   Double_t(fNpn[iAnalysisType][p2])*TMath::Power((0.5/Double_t(fNp[iAnalysisType]) + 0.5/Double_t(fNn[iAnalysisType])),2))/fP2Step[iAnalysisType];
  
  return ferror[iAnalysisType][p2];
}
//____________________________________________________________________//
TGraphErrors *AliBalance::DrawBalance(Int_t iAnalysisType) {

  // Draws the BF
  Double_t x[MAXIMUM_NUMBER_OF_STEPS];
  Double_t xer[MAXIMUM_NUMBER_OF_STEPS];
  Double_t b[MAXIMUM_NUMBER_OF_STEPS];
  Double_t ber[MAXIMUM_NUMBER_OF_STEPS];

  if((fNp[iAnalysisType] == 0)||(fNn[iAnalysisType] == 0)) {
    cerr<<"Couldn't find any particles in the analyzed interval!!!"<<endl;
    return NULL;
  }
  
  for(Int_t i = 0; i < fNumberOfBins[iAnalysisType]; i++) {
    b[i] = GetBalance(iAnalysisType,i);
    ber[i] = GetError(iAnalysisType,i);
    x[i] = fP2Start[iAnalysisType] + fP2Step[iAnalysisType]*i + fP2Step[iAnalysisType]/2;
    xer[i] = 0.0;
  }
  
  TGraphErrors *gr = new TGraphErrors(fNumberOfBins[iAnalysisType],x,b,xer,ber);
  gr->GetXaxis()->SetTitleColor(1);
  if(iAnalysisType==0) {
    gr->SetTitle("Balance function B(#Delta y)");
    gr->GetXaxis()->SetTitle("#Delta y");
    gr->GetYaxis()->SetTitle("B(#Delta y)");
  }
  if(iAnalysisType==1) {
    gr->SetTitle("Balance function B(#Delta #eta)");
    gr->GetXaxis()->SetTitle("#Delta #eta");
    gr->GetYaxis()->SetTitle("B(#Delta #eta)");
  }
  if(iAnalysisType==2) {
    gr->SetTitle("Balance function B(q_{long})");
    gr->GetXaxis()->SetTitle("q_{long} (GeV/c)");
    gr->GetYaxis()->SetTitle("B(q_{long}) ((GeV/c)^{-1})");
  }
  if(iAnalysisType==3) {
    gr->SetTitle("Balance function B(q_{out})");
    gr->GetXaxis()->SetTitle("q_{out} (GeV/c)");
    gr->GetYaxis()->SetTitle("B(q_{out}) ((GeV/c)^{-1})");
  }
  if(iAnalysisType==4) {
    gr->SetTitle("Balance function B(q_{side})");
    gr->GetXaxis()->SetTitle("q_{side} (GeV/c)");
    gr->GetYaxis()->SetTitle("B(q_{side}) ((GeV/c)^{-1})");
  }
  if(iAnalysisType==5) {
    gr->SetTitle("Balance function B(q_{inv})");
    gr->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
    gr->GetYaxis()->SetTitle("B(q_{inv}) ((GeV/c)^{-1})");
  }
  if(iAnalysisType==6) {
    gr->SetTitle("Balance function B(#Delta #phi)");
    gr->GetXaxis()->SetTitle("#Delta #phi");
    gr->GetYaxis()->SetTitle("B(#Delta #phi)");
  }

  return gr;
}

//____________________________________________________________________//
void AliBalance::PrintResults(Int_t iAnalysisType, TH1D *gHistBalance) {
  //Prints the calculated width of the BF and its error
  Double_t gSumXi = 0.0, gSumBi = 0.0, gSumBiXi = 0.0;
  Double_t gSumBiXi2 = 0.0, gSumBi2Xi2 = 0.0;
  Double_t gSumDeltaBi2 = 0.0, gSumXi2DeltaBi2 = 0.0;
  Double_t deltaBalP2 = 0.0, integral = 0.0;
  Double_t deltaErrorNew = 0.0;
  
  // cout<<"=================================================="<<endl;
  // for(Int_t i = 1; i <= fNumberOfBins[iAnalysisType]; i++) { 
  //   x[i-1] = fP2Start[iAnalysisType] + fP2Step[iAnalysisType]*i + fP2Step[iAnalysisType]/2;
  //   cout<<"B: "<<gHistBalance->GetBinContent(i)<<"\t Error: "<<gHistBalance->GetBinError(i)<<"\t bin: "<<gHistBalance->GetBinCenter(i)<<endl;
  // } 
  // cout<<"=================================================="<<endl;
  for(Int_t i = 2; i <= fNumberOfBins[iAnalysisType]; i++) {
    gSumXi += gHistBalance->GetBinCenter(i);
    gSumBi += gHistBalance->GetBinContent(i);
    gSumBiXi += gHistBalance->GetBinContent(i)*gHistBalance->GetBinCenter(i);
    gSumBiXi2 += gHistBalance->GetBinContent(i)*TMath::Power(gHistBalance->GetBinCenter(i),2);
    gSumBi2Xi2 += TMath::Power(gHistBalance->GetBinContent(i),2)*TMath::Power(gHistBalance->GetBinCenter(i),2);
    gSumDeltaBi2 +=  TMath::Power(gHistBalance->GetBinError(i),2);
    gSumXi2DeltaBi2 += TMath::Power(gHistBalance->GetBinCenter(i),2) * TMath::Power(gHistBalance->GetBinError(i),2);
    
    deltaBalP2 += fP2Step[iAnalysisType]*TMath::Power(gHistBalance->GetBinError(i),2);
    integral += fP2Step[iAnalysisType]*gHistBalance->GetBinContent(i);
  }
  for(Int_t i = 1; i < fNumberOfBins[iAnalysisType]; i++)
    deltaErrorNew += gHistBalance->GetBinError(i)*(gHistBalance->GetBinCenter(i)*gSumBi - gSumBiXi)/TMath::Power(gSumBi,2);
  
  Double_t integralError = TMath::Sqrt(deltaBalP2);
  integralError *= 1.0;

  Double_t delta = gSumBiXi / gSumBi; delta *= 1.0;
  Double_t deltaError = (gSumBiXi / gSumBi) * TMath::Sqrt(TMath::Power((TMath::Sqrt(gSumXi2DeltaBi2)/gSumBiXi),2) + TMath::Power((gSumDeltaBi2/gSumBi),2) );
  deltaError *= 1.0;
  // cout<<"Analysis type: "<<kBFAnalysisType[iAnalysisType].Data()<<endl;
  // cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
  // cout<<"New error: "<<deltaErrorNew<<endl;
  // cout<<"Integral: "<<integral<<"\t Error: "<<integralError<<endl;
  // cout<<"=================================================="<<endl;
}
 
//____________________________________________________________________//
TH1D *AliBalance::GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t etaWindow,Bool_t correctWithEfficiency, Bool_t correctWithAcceptanceOnly, Bool_t correctWithMixed, TH1D *hMixed[4]) {
  //Returns the BF histogram, extracted from the 6 TH2D objects 
  //(private members) of the AliBalance class.
  //
  // Acceptance correction: 
  // - only for analysis type = kEta
  // - only if etaWindow > 0 (default = -1.)
  // - calculated as proposed by STAR 
  //
  TString gAnalysisType[ANALYSIS_TYPES] = {"y","eta","qlong","qout","qside","qinv","phi"};
  TString histName = "gHistBalanceFunctionHistogram";
  histName += gAnalysisType[iAnalysisType];

  SetInterval(iAnalysisType, fHistP[iAnalysisType]->GetYaxis()->GetXmin(),
	      fHistP[iAnalysisType]->GetYaxis()->GetXmin(),
	      fHistPP[iAnalysisType]->GetNbinsY(),
	      fHistPP[iAnalysisType]->GetYaxis()->GetXmin(),
	      fHistPP[iAnalysisType]->GetYaxis()->GetXmax());

  // determine the projection thresholds
  Int_t binMinX, binMinY, binMinZ;
  Int_t binMaxX, binMaxY, binMaxZ;

  fHistPP[iAnalysisType]->GetBinXYZ(fHistPP[iAnalysisType]->FindBin(centrMin),binMinX,binMinY,binMinZ);
  fHistPP[iAnalysisType]->GetBinXYZ(fHistPP[iAnalysisType]->FindBin(centrMax),binMaxX,binMaxY,binMaxZ);

  TH1D *gHistBalanceFunctionHistogram = new TH1D(histName.Data(),"",fHistPP[iAnalysisType]->GetNbinsY(),fHistPP[iAnalysisType]->GetYaxis()->GetXmin(),fHistPP[iAnalysisType]->GetYaxis()->GetXmax());
  switch(iAnalysisType) {
  case kRapidity:
    gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta y");
    gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta y)");
    break;
  case kEta:
    gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #eta");
    gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #eta)");
    break;
  case kQlong:
    gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{long} (GeV/c)");
    gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{long})");
    break;
  case kQout:
    gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{out} (GeV/c)");
    gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{out})");
    break;
  case kQside:
    gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{side} (GeV/c)");
    gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{side})");
    break;
  case kQinv:
    gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
    gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(q_{inv})");
    break;
  case kPhi:
    gHistBalanceFunctionHistogram->GetXaxis()->SetTitle("#Delta #phi (deg.)");
    gHistBalanceFunctionHistogram->GetYaxis()->SetTitle("B(#Delta #phi)");
    break;
  default:
    break;
  }

  TH1D *hTemp1 = dynamic_cast<TH1D *>(fHistPN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistPN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
  TH1D *hTemp2 = dynamic_cast<TH1D *>(fHistPN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f_copy",fHistPN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
  TH1D *hTemp3 = dynamic_cast<TH1D *>(fHistNN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistNN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
  TH1D *hTemp4 = dynamic_cast<TH1D *>(fHistPP[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistPP[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
  TH1D *hTemp5 = dynamic_cast<TH1D *>(fHistN[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistN[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));
  TH1D *hTemp6 = dynamic_cast<TH1D *>(fHistP[iAnalysisType]->ProjectionY(Form("%s_Cent_%.0f_%.0f",fHistP[iAnalysisType]->GetName(),centrMin,centrMax),binMinX,binMaxX));

  // get the file with the efficiency matrices
  // withAcceptanceOnly: Data single distributions are normalized to 1 (efficiency not taken into account)
  // else : Data single distributions are normalized to give single particle efficiency of MC
  TFile *fEfficiencyMatrix = NULL;
  if(correctWithEfficiency || correctWithMixed){
    if(correctWithAcceptanceOnly) fEfficiencyMatrix = TFile::Open("$ALICE_ROOT/PWGCF/EBYE/macros/accOnlyFromConvolutionAllCent.root");
    else  fEfficiencyMatrix = TFile::Open("$ALICE_ROOT/PWGCF/EBYE/macros/effFromConvolutionAllCent.root");
    if(!fEfficiencyMatrix){
      AliError("Efficiency histogram file not found");
      return NULL;
    }
  }

  // do correction with the efficiency calculated from MC + Data (for single particles and two particle correlations)
  // - single particle efficiencies from MC (AliAnalysiTaskEfficiency)
  // - two particle efficiencies from convolution of data single particle distributions 
  //   (normalized to single particle efficiency)
  if(iAnalysisType == kEta && etaWindow > 0 && correctWithEfficiency && !correctWithMixed){

    TH1F* hEffP  = NULL;
    TH1F* hEffN  = NULL;
    TH1F* hEffPP = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffPP_Cent%.0f-%.0f_Data",centrMin,centrMax));
    TH1F* hEffNN = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffNN_Cent%.0f-%.0f_Data",centrMin,centrMax));
    TH1F* hEffPN = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffPN_Cent%.0f-%.0f_Data",centrMin,centrMax));

    // take the data distributions
    if(correctWithAcceptanceOnly){
      hEffP = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffP_Cent%.0f-%.0f_Data",centrMin,centrMax));
      hEffN = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffN_Cent%.0f-%.0f_Data",centrMin,centrMax));
    }
    // take the MC distributions
    else{
      hEffP = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffP_Cent%.0f-%.0f_MC",centrMin,centrMax));
      hEffN = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffN_Cent%.0f-%.0f_MC",centrMin,centrMax));
    }

    if( !hEffP || !hEffN || !hEffPP || !hEffNN || !hEffPN){
      AliError(Form("Efficiency (eta) histograms not found: etaEffPP_Cent%.0f-%.0f_Data",centrMin,centrMax));
      return NULL;
    }

    for(Int_t iBin = 0; iBin < hEffP->GetNbinsX(); iBin++){
      hTemp5->SetBinError(iBin+1,hTemp5->GetBinError(iBin+1)/hEffN->GetBinContent(hEffN->FindBin(hTemp5->GetBinCenter(iBin+1))));
      hTemp5->SetBinContent(iBin+1,hTemp5->GetBinContent(iBin+1)/hEffN->GetBinContent(hEffN->FindBin(hTemp5->GetBinCenter(iBin+1))));

      hTemp6->SetBinError(iBin+1,hTemp6->GetBinError(iBin+1)/hEffP->GetBinContent(hEffP->FindBin(hTemp6->GetBinCenter(iBin+1))));
      hTemp6->SetBinContent(iBin+1,hTemp6->GetBinContent(iBin+1)/hEffP->GetBinContent(hEffP->FindBin(hTemp6->GetBinCenter(iBin+1))));
    }
  
    for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){

      hTemp1->SetBinError(iBin+1,hTemp1->GetBinError(iBin+1)/hEffPN->GetBinContent(hEffPN->FindBin(hTemp1->GetBinCenter(iBin+1))));
      hTemp1->SetBinContent(iBin+1,hTemp1->GetBinContent(iBin+1)/hEffPN->GetBinContent(hEffPN->FindBin(hTemp1->GetBinCenter(iBin+1))));
      hTemp2->SetBinError(iBin+1,hTemp2->GetBinError(iBin+1)/hEffPN->GetBinContent(hEffPN->FindBin(hTemp2->GetBinCenter(iBin+1))));
      hTemp2->SetBinContent(iBin+1,hTemp2->GetBinContent(iBin+1)/hEffPN->GetBinContent(hEffPN->FindBin(hTemp2->GetBinCenter(iBin+1))));
      hTemp3->SetBinError(iBin+1,hTemp3->GetBinError(iBin+1)/hEffNN->GetBinContent(hEffNN->FindBin(hTemp3->GetBinCenter(iBin+1))));
      hTemp3->SetBinContent(iBin+1,hTemp3->GetBinContent(iBin+1)/hEffNN->GetBinContent(hEffNN->FindBin(hTemp3->GetBinCenter(iBin+1))));
      hTemp4->SetBinError(iBin+1,hTemp4->GetBinError(iBin+1)/hEffPP->GetBinContent(hEffPP->FindBin(hTemp4->GetBinCenter(iBin+1))));      
      hTemp4->SetBinContent(iBin+1,hTemp4->GetBinContent(iBin+1)/hEffPP->GetBinContent(hEffPP->FindBin(hTemp4->GetBinCenter(iBin+1))));
      
    }

    // TF1 *fPP = new TF1("fPP","pol1",0,1.6);  // phase space factor + efficiency for ++
    // fPP->SetParameters(0.736466,-0.461529);
    // TF1 *fNN = new TF1("fNN","pol1",0,1.6);  // phase space factor + efficiency for --
    // fNN->SetParameters(0.718616,-0.450473);
    // TF1 *fPN = new TF1("fPN","pol1",0,1.6);  // phase space factor + efficiency for +-
    // fPN->SetParameters(0.727507,-0.455981);
    
    // for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){
    //   hTemp1->SetBinContent(iBin+1,hTemp1->GetBinContent(iBin+1)/fPN->Eval(hTemp1->GetBinCenter(iBin+1)));
    //   hTemp1->SetBinError(iBin+1,hTemp1->GetBinError(iBin+1)/fPN->Eval(hTemp1->GetBinCenter(iBin+1)));
    //   hTemp2->SetBinContent(iBin+1,hTemp2->GetBinContent(iBin+1)/fPN->Eval(hTemp1->GetBinCenter(iBin+1)));
    //   hTemp2->SetBinError(iBin+1,hTemp2->GetBinError(iBin+1)/fPN->Eval(hTemp1->GetBinCenter(iBin+1)));
    //   hTemp3->SetBinContent(iBin+1,hTemp3->GetBinContent(iBin+1)/fNN->Eval(hTemp1->GetBinCenter(iBin+1)));
    //   hTemp3->SetBinError(iBin+1,hTemp3->GetBinError(iBin+1)/fNN->Eval(hTemp1->GetBinCenter(iBin+1)));
    //   hTemp4->SetBinContent(iBin+1,hTemp4->GetBinContent(iBin+1)/fPP->Eval(hTemp1->GetBinCenter(iBin+1)));
    //   hTemp4->SetBinError(iBin+1,hTemp4->GetBinError(iBin+1)/fPP->Eval(hTemp1->GetBinCenter(iBin+1)));
    // }      
  }

  // do correction with the efficiency calculated from MC + Data (for single particles and two particle correlations)
  // - single particle efficiencies from MC (AliAnalysiTaskEfficiency)
  // - two particle efficiencies from convolution of data single particle distributions 
  //   (normalized to single particle efficiency)  
  if(iAnalysisType == kPhi && correctWithEfficiency && !correctWithMixed){

    TH1F* hEffPhiP  = NULL;
    TH1F* hEffPhiN  = NULL;
    TH1F* hEffPhiPP = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffPP_Cent%.0f-%.0f_Data",centrMin,centrMax));
    TH1F* hEffPhiNN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffNN_Cent%.0f-%.0f_Data",centrMin,centrMax));
    TH1F* hEffPhiPN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffPN_Cent%.0f-%.0f_Data",centrMin,centrMax));

    // take the data distributions
    if(correctWithAcceptanceOnly){
      hEffPhiP = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffP_Cent%.0f-%.0f_Data",centrMin,centrMax));
      hEffPhiN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffN_Cent%.0f-%.0f_Data",centrMin,centrMax));
    }
    // take the MC distributions
    else{
      hEffPhiP = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffP_Cent%.0f-%.0f_MC",centrMin,centrMax));
      hEffPhiN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffN_Cent%.0f-%.0f_MC",centrMin,centrMax));
    }

    if( !hEffPhiP || !hEffPhiN || !hEffPhiPP || !hEffPhiNN || !hEffPhiPN){
      AliError("Efficiency (phi) histograms not found");
      return NULL;
    }

    for(Int_t iBin = 0; iBin < hEffPhiP->GetNbinsX(); iBin++){
      hTemp5->SetBinError(iBin+1,hTemp5->GetBinError(iBin+1)/hEffPhiN->GetBinContent(hEffPhiN->FindBin(hTemp5->GetBinCenter(iBin+1))));
      hTemp5->SetBinContent(iBin+1,hTemp5->GetBinContent(iBin+1)/hEffPhiN->GetBinContent(hEffPhiN->FindBin(hTemp5->GetBinCenter(iBin+1))));

      hTemp6->SetBinError(iBin+1,hTemp6->GetBinError(iBin+1)/hEffPhiP->GetBinContent(hEffPhiP->FindBin(hTemp6->GetBinCenter(iBin+1))));
      hTemp6->SetBinContent(iBin+1,hTemp6->GetBinContent(iBin+1)/hEffPhiP->GetBinContent(hEffPhiP->FindBin(hTemp6->GetBinCenter(iBin+1))));
    }

    for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){

      hTemp1->SetBinError(iBin+1,hTemp1->GetBinError(iBin+1)/hEffPhiPN->GetBinContent(hEffPhiPN->FindBin(hTemp1->GetBinCenter(iBin+1))));
      hTemp1->SetBinContent(iBin+1,hTemp1->GetBinContent(iBin+1)/hEffPhiPN->GetBinContent(hEffPhiPN->FindBin(hTemp1->GetBinCenter(iBin+1))));
      hTemp2->SetBinError(iBin+1,hTemp2->GetBinError(iBin+1)/hEffPhiPN->GetBinContent(hEffPhiPN->FindBin(hTemp2->GetBinCenter(iBin+1))));
      hTemp2->SetBinContent(iBin+1,hTemp2->GetBinContent(iBin+1)/hEffPhiPN->GetBinContent(hEffPhiPN->FindBin(hTemp2->GetBinCenter(iBin+1))));
      hTemp3->SetBinError(iBin+1,hTemp3->GetBinError(iBin+1)/hEffPhiNN->GetBinContent(hEffPhiNN->FindBin(hTemp3->GetBinCenter(iBin+1))));
      hTemp3->SetBinContent(iBin+1,hTemp3->GetBinContent(iBin+1)/hEffPhiNN->GetBinContent(hEffPhiNN->FindBin(hTemp3->GetBinCenter(iBin+1))));
      hTemp4->SetBinError(iBin+1,hTemp4->GetBinError(iBin+1)/hEffPhiPP->GetBinContent(hEffPhiPP->FindBin(hTemp4->GetBinCenter(iBin+1))));      
      hTemp4->SetBinContent(iBin+1,hTemp4->GetBinContent(iBin+1)/hEffPhiPP->GetBinContent(hEffPhiPP->FindBin(hTemp4->GetBinCenter(iBin+1))));
      
    }  
  }

  // do the correction with the event mixing directly!
  if(correctWithMixed){
  
    // take the MC distributions (for average efficiency)
    TH1F* hEffP = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffP_Cent%.0f-%.0f_MC",centrMin,centrMax));
    TH1F* hEffN = (TH1F*)fEfficiencyMatrix->Get(Form("etaEffN_Cent%.0f-%.0f_MC",centrMin,centrMax));  

    TH1F* hEffPP = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffPP_Cent%.0f-%.0f_Data",centrMin,centrMax));
    TH1F* hEffNN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffNN_Cent%.0f-%.0f_Data",centrMin,centrMax));
    TH1F* hEffPN = (TH1F*)fEfficiencyMatrix->Get(Form("phiEffPN_Cent%.0f-%.0f_Data",centrMin,centrMax));

    if( !hEffP || !hEffN){
      AliError(Form("Efficiency (eta) histograms not found: etaEffPP_Cent%.0f-%.0f_Data",centrMin,centrMax));
      return NULL;
    }

    if(hMixed[0] && hMixed[1] && hMixed[2] && hMixed[3]){
    
      // scale to average efficiency in the pt region (0.3-1.5) and |eta| < 0.8
      // by multiplying the average single particle efficiencies from HIJING
      // here we assume that the distributions are 1:
      // - in the integral for dphi (for averaging over sector structure)
      // - in the maximum for deta
      Double_t normPMC = (Double_t)hEffP->Integral()/(Double_t)hEffP->GetNbinsX();
      Double_t normNMC = (Double_t)hEffN->Integral()/(Double_t)hEffN->GetNbinsX();
      Double_t normPPMC = (Double_t)hEffPP->Integral()/(Double_t)hEffPP->GetNbinsX();
      Double_t normNNMC = (Double_t)hEffNN->Integral()/(Double_t)hEffNN->GetNbinsX();
      Double_t normPNMC = (Double_t)hEffPN->Integral()/(Double_t)hEffPN->GetNbinsX();

      hMixed[0]->Scale(normPNMC);
      hMixed[1]->Scale(normPNMC);
      hMixed[2]->Scale(normNNMC);
      hMixed[3]->Scale(normPPMC);

      // divide by event mixing
      hTemp1->Divide(hMixed[0]);
      hTemp2->Divide(hMixed[1]);
      hTemp3->Divide(hMixed[2]);
      hTemp4->Divide(hMixed[3]);

      // scale also single histograms with average efficiency
      hTemp5->Scale(1./normNMC);
      hTemp6->Scale(1./normPMC);

    }
    else{
      AliError("Correction with EventMixing requested, but not all Histograms there!");
      return NULL;
    }
  }


  if((hTemp1)&&(hTemp2)&&(hTemp3)&&(hTemp4)) {
    hTemp1->Sumw2();
    hTemp2->Sumw2();
    hTemp3->Sumw2();
    hTemp4->Sumw2();
    hTemp1->Add(hTemp3,-2.);
    hTemp1->Scale(1./hTemp5->Integral());
    hTemp2->Add(hTemp4,-2.);
    hTemp2->Scale(1./hTemp6->Integral());
    gHistBalanceFunctionHistogram->Add(hTemp1,hTemp2,1.,1.);
    gHistBalanceFunctionHistogram->Scale(0.5/fP2Step[iAnalysisType]);
  }

  // do the acceptance correction (only for Eta and etaWindow > 0)
  if(iAnalysisType == kEta && etaWindow > 0 && !correctWithEfficiency && !correctWithMixed){
    for(Int_t iBin = 0; iBin < gHistBalanceFunctionHistogram->GetNbinsX(); iBin++){
      
      Double_t notCorrected = gHistBalanceFunctionHistogram->GetBinContent(iBin+1);
      Double_t corrected    = notCorrected / (1 - (gHistBalanceFunctionHistogram->GetBinCenter(iBin+1))/ etaWindow );
      gHistBalanceFunctionHistogram->SetBinContent(iBin+1, corrected);
      gHistBalanceFunctionHistogram->SetBinError(iBin+1,corrected/notCorrected*gHistBalanceFunctionHistogram->GetBinError(iBin+1));
      
    }
  }
  
  if(fEfficiencyMatrix)   fEfficiencyMatrix->Close();

  PrintResults(iAnalysisType,gHistBalanceFunctionHistogram);

  return gHistBalanceFunctionHistogram;
}
 AliBalance.cxx:1
 AliBalance.cxx:2
 AliBalance.cxx:3
 AliBalance.cxx:4
 AliBalance.cxx:5
 AliBalance.cxx:6
 AliBalance.cxx:7
 AliBalance.cxx:8
 AliBalance.cxx:9
 AliBalance.cxx:10
 AliBalance.cxx:11
 AliBalance.cxx:12
 AliBalance.cxx:13
 AliBalance.cxx:14
 AliBalance.cxx:15
 AliBalance.cxx:16
 AliBalance.cxx:17
 AliBalance.cxx:18
 AliBalance.cxx:19
 AliBalance.cxx:20
 AliBalance.cxx:21
 AliBalance.cxx:22
 AliBalance.cxx:23
 AliBalance.cxx:24
 AliBalance.cxx:25
 AliBalance.cxx:26
 AliBalance.cxx:27
 AliBalance.cxx:28
 AliBalance.cxx:29
 AliBalance.cxx:30
 AliBalance.cxx:31
 AliBalance.cxx:32
 AliBalance.cxx:33
 AliBalance.cxx:34
 AliBalance.cxx:35
 AliBalance.cxx:36
 AliBalance.cxx:37
 AliBalance.cxx:38
 AliBalance.cxx:39
 AliBalance.cxx:40
 AliBalance.cxx:41
 AliBalance.cxx:42
 AliBalance.cxx:43
 AliBalance.cxx:44
 AliBalance.cxx:45
 AliBalance.cxx:46
 AliBalance.cxx:47
 AliBalance.cxx:48
 AliBalance.cxx:49
 AliBalance.cxx:50
 AliBalance.cxx:51
 AliBalance.cxx:52
 AliBalance.cxx:53
 AliBalance.cxx:54
 AliBalance.cxx:55
 AliBalance.cxx:56
 AliBalance.cxx:57
 AliBalance.cxx:58
 AliBalance.cxx:59
 AliBalance.cxx:60
 AliBalance.cxx:61
 AliBalance.cxx:62
 AliBalance.cxx:63
 AliBalance.cxx:64
 AliBalance.cxx:65
 AliBalance.cxx:66
 AliBalance.cxx:67
 AliBalance.cxx:68
 AliBalance.cxx:69
 AliBalance.cxx:70
 AliBalance.cxx:71
 AliBalance.cxx:72
 AliBalance.cxx:73
 AliBalance.cxx:74
 AliBalance.cxx:75
 AliBalance.cxx:76
 AliBalance.cxx:77
 AliBalance.cxx:78
 AliBalance.cxx:79
 AliBalance.cxx:80
 AliBalance.cxx:81
 AliBalance.cxx:82
 AliBalance.cxx:83
 AliBalance.cxx:84
 AliBalance.cxx:85
 AliBalance.cxx:86
 AliBalance.cxx:87
 AliBalance.cxx:88
 AliBalance.cxx:89
 AliBalance.cxx:90
 AliBalance.cxx:91
 AliBalance.cxx:92
 AliBalance.cxx:93
 AliBalance.cxx:94
 AliBalance.cxx:95
 AliBalance.cxx:96
 AliBalance.cxx:97
 AliBalance.cxx:98
 AliBalance.cxx:99
 AliBalance.cxx:100
 AliBalance.cxx:101
 AliBalance.cxx:102
 AliBalance.cxx:103
 AliBalance.cxx:104
 AliBalance.cxx:105
 AliBalance.cxx:106
 AliBalance.cxx:107
 AliBalance.cxx:108
 AliBalance.cxx:109
 AliBalance.cxx:110
 AliBalance.cxx:111
 AliBalance.cxx:112
 AliBalance.cxx:113
 AliBalance.cxx:114
 AliBalance.cxx:115
 AliBalance.cxx:116
 AliBalance.cxx:117
 AliBalance.cxx:118
 AliBalance.cxx:119
 AliBalance.cxx:120
 AliBalance.cxx:121
 AliBalance.cxx:122
 AliBalance.cxx:123
 AliBalance.cxx:124
 AliBalance.cxx:125
 AliBalance.cxx:126
 AliBalance.cxx:127
 AliBalance.cxx:128
 AliBalance.cxx:129
 AliBalance.cxx:130
 AliBalance.cxx:131
 AliBalance.cxx:132
 AliBalance.cxx:133
 AliBalance.cxx:134
 AliBalance.cxx:135
 AliBalance.cxx:136
 AliBalance.cxx:137
 AliBalance.cxx:138
 AliBalance.cxx:139
 AliBalance.cxx:140
 AliBalance.cxx:141
 AliBalance.cxx:142
 AliBalance.cxx:143
 AliBalance.cxx:144
 AliBalance.cxx:145
 AliBalance.cxx:146
 AliBalance.cxx:147
 AliBalance.cxx:148
 AliBalance.cxx:149
 AliBalance.cxx:150
 AliBalance.cxx:151
 AliBalance.cxx:152
 AliBalance.cxx:153
 AliBalance.cxx:154
 AliBalance.cxx:155
 AliBalance.cxx:156
 AliBalance.cxx:157
 AliBalance.cxx:158
 AliBalance.cxx:159
 AliBalance.cxx:160
 AliBalance.cxx:161
 AliBalance.cxx:162
 AliBalance.cxx:163
 AliBalance.cxx:164
 AliBalance.cxx:165
 AliBalance.cxx:166
 AliBalance.cxx:167
 AliBalance.cxx:168
 AliBalance.cxx:169
 AliBalance.cxx:170
 AliBalance.cxx:171
 AliBalance.cxx:172
 AliBalance.cxx:173
 AliBalance.cxx:174
 AliBalance.cxx:175
 AliBalance.cxx:176
 AliBalance.cxx:177
 AliBalance.cxx:178
 AliBalance.cxx:179
 AliBalance.cxx:180
 AliBalance.cxx:181
 AliBalance.cxx:182
 AliBalance.cxx:183
 AliBalance.cxx:184
 AliBalance.cxx:185
 AliBalance.cxx:186
 AliBalance.cxx:187
 AliBalance.cxx:188
 AliBalance.cxx:189
 AliBalance.cxx:190
 AliBalance.cxx:191
 AliBalance.cxx:192
 AliBalance.cxx:193
 AliBalance.cxx:194
 AliBalance.cxx:195
 AliBalance.cxx:196
 AliBalance.cxx:197
 AliBalance.cxx:198
 AliBalance.cxx:199
 AliBalance.cxx:200
 AliBalance.cxx:201
 AliBalance.cxx:202
 AliBalance.cxx:203
 AliBalance.cxx:204
 AliBalance.cxx:205
 AliBalance.cxx:206
 AliBalance.cxx:207
 AliBalance.cxx:208
 AliBalance.cxx:209
 AliBalance.cxx:210
 AliBalance.cxx:211
 AliBalance.cxx:212
 AliBalance.cxx:213
 AliBalance.cxx:214
 AliBalance.cxx:215
 AliBalance.cxx:216
 AliBalance.cxx:217
 AliBalance.cxx:218
 AliBalance.cxx:219
 AliBalance.cxx:220
 AliBalance.cxx:221
 AliBalance.cxx:222
 AliBalance.cxx:223
 AliBalance.cxx:224
 AliBalance.cxx:225
 AliBalance.cxx:226
 AliBalance.cxx:227
 AliBalance.cxx:228
 AliBalance.cxx:229
 AliBalance.cxx:230
 AliBalance.cxx:231
 AliBalance.cxx:232
 AliBalance.cxx:233
 AliBalance.cxx:234
 AliBalance.cxx:235
 AliBalance.cxx:236
 AliBalance.cxx:237
 AliBalance.cxx:238
 AliBalance.cxx:239
 AliBalance.cxx:240
 AliBalance.cxx:241
 AliBalance.cxx:242
 AliBalance.cxx:243
 AliBalance.cxx:244
 AliBalance.cxx:245
 AliBalance.cxx:246
 AliBalance.cxx:247
 AliBalance.cxx:248
 AliBalance.cxx:249
 AliBalance.cxx:250
 AliBalance.cxx:251
 AliBalance.cxx:252
 AliBalance.cxx:253
 AliBalance.cxx:254
 AliBalance.cxx:255
 AliBalance.cxx:256
 AliBalance.cxx:257
 AliBalance.cxx:258
 AliBalance.cxx:259
 AliBalance.cxx:260
 AliBalance.cxx:261
 AliBalance.cxx:262
 AliBalance.cxx:263
 AliBalance.cxx:264
 AliBalance.cxx:265
 AliBalance.cxx:266
 AliBalance.cxx:267
 AliBalance.cxx:268
 AliBalance.cxx:269
 AliBalance.cxx:270
 AliBalance.cxx:271
 AliBalance.cxx:272
 AliBalance.cxx:273
 AliBalance.cxx:274
 AliBalance.cxx:275
 AliBalance.cxx:276
 AliBalance.cxx:277
 AliBalance.cxx:278
 AliBalance.cxx:279
 AliBalance.cxx:280
 AliBalance.cxx:281
 AliBalance.cxx:282
 AliBalance.cxx:283
 AliBalance.cxx:284
 AliBalance.cxx:285
 AliBalance.cxx:286
 AliBalance.cxx:287
 AliBalance.cxx:288
 AliBalance.cxx:289
 AliBalance.cxx:290
 AliBalance.cxx:291
 AliBalance.cxx:292
 AliBalance.cxx:293
 AliBalance.cxx:294
 AliBalance.cxx:295
 AliBalance.cxx:296
 AliBalance.cxx:297
 AliBalance.cxx:298
 AliBalance.cxx:299
 AliBalance.cxx:300
 AliBalance.cxx:301
 AliBalance.cxx:302
 AliBalance.cxx:303
 AliBalance.cxx:304
 AliBalance.cxx:305
 AliBalance.cxx:306
 AliBalance.cxx:307
 AliBalance.cxx:308
 AliBalance.cxx:309
 AliBalance.cxx:310
 AliBalance.cxx:311
 AliBalance.cxx:312
 AliBalance.cxx:313
 AliBalance.cxx:314
 AliBalance.cxx:315
 AliBalance.cxx:316
 AliBalance.cxx:317
 AliBalance.cxx:318
 AliBalance.cxx:319
 AliBalance.cxx:320
 AliBalance.cxx:321
 AliBalance.cxx:322
 AliBalance.cxx:323
 AliBalance.cxx:324
 AliBalance.cxx:325
 AliBalance.cxx:326
 AliBalance.cxx:327
 AliBalance.cxx:328
 AliBalance.cxx:329
 AliBalance.cxx:330
 AliBalance.cxx:331
 AliBalance.cxx:332
 AliBalance.cxx:333
 AliBalance.cxx:334
 AliBalance.cxx:335
 AliBalance.cxx:336
 AliBalance.cxx:337
 AliBalance.cxx:338
 AliBalance.cxx:339
 AliBalance.cxx:340
 AliBalance.cxx:341
 AliBalance.cxx:342
 AliBalance.cxx:343
 AliBalance.cxx:344
 AliBalance.cxx:345
 AliBalance.cxx:346
 AliBalance.cxx:347
 AliBalance.cxx:348
 AliBalance.cxx:349
 AliBalance.cxx:350
 AliBalance.cxx:351
 AliBalance.cxx:352
 AliBalance.cxx:353
 AliBalance.cxx:354
 AliBalance.cxx:355
 AliBalance.cxx:356
 AliBalance.cxx:357
 AliBalance.cxx:358
 AliBalance.cxx:359
 AliBalance.cxx:360
 AliBalance.cxx:361
 AliBalance.cxx:362
 AliBalance.cxx:363
 AliBalance.cxx:364
 AliBalance.cxx:365
 AliBalance.cxx:366
 AliBalance.cxx:367
 AliBalance.cxx:368
 AliBalance.cxx:369
 AliBalance.cxx:370
 AliBalance.cxx:371
 AliBalance.cxx:372
 AliBalance.cxx:373
 AliBalance.cxx:374
 AliBalance.cxx:375
 AliBalance.cxx:376
 AliBalance.cxx:377
 AliBalance.cxx:378
 AliBalance.cxx:379
 AliBalance.cxx:380
 AliBalance.cxx:381
 AliBalance.cxx:382
 AliBalance.cxx:383
 AliBalance.cxx:384
 AliBalance.cxx:385
 AliBalance.cxx:386
 AliBalance.cxx:387
 AliBalance.cxx:388
 AliBalance.cxx:389
 AliBalance.cxx:390
 AliBalance.cxx:391
 AliBalance.cxx:392
 AliBalance.cxx:393
 AliBalance.cxx:394
 AliBalance.cxx:395
 AliBalance.cxx:396
 AliBalance.cxx:397
 AliBalance.cxx:398
 AliBalance.cxx:399
 AliBalance.cxx:400
 AliBalance.cxx:401
 AliBalance.cxx:402
 AliBalance.cxx:403
 AliBalance.cxx:404
 AliBalance.cxx:405
 AliBalance.cxx:406
 AliBalance.cxx:407
 AliBalance.cxx:408
 AliBalance.cxx:409
 AliBalance.cxx:410
 AliBalance.cxx:411
 AliBalance.cxx:412
 AliBalance.cxx:413
 AliBalance.cxx:414
 AliBalance.cxx:415
 AliBalance.cxx:416
 AliBalance.cxx:417
 AliBalance.cxx:418
 AliBalance.cxx:419
 AliBalance.cxx:420
 AliBalance.cxx:421
 AliBalance.cxx:422
 AliBalance.cxx:423
 AliBalance.cxx:424
 AliBalance.cxx:425
 AliBalance.cxx:426
 AliBalance.cxx:427
 AliBalance.cxx:428
 AliBalance.cxx:429
 AliBalance.cxx:430
 AliBalance.cxx:431
 AliBalance.cxx:432
 AliBalance.cxx:433
 AliBalance.cxx:434
 AliBalance.cxx:435
 AliBalance.cxx:436
 AliBalance.cxx:437
 AliBalance.cxx:438
 AliBalance.cxx:439
 AliBalance.cxx:440
 AliBalance.cxx:441
 AliBalance.cxx:442
 AliBalance.cxx:443
 AliBalance.cxx:444
 AliBalance.cxx:445
 AliBalance.cxx:446
 AliBalance.cxx:447
 AliBalance.cxx:448
 AliBalance.cxx:449
 AliBalance.cxx:450
 AliBalance.cxx:451
 AliBalance.cxx:452
 AliBalance.cxx:453
 AliBalance.cxx:454
 AliBalance.cxx:455
 AliBalance.cxx:456
 AliBalance.cxx:457
 AliBalance.cxx:458
 AliBalance.cxx:459
 AliBalance.cxx:460
 AliBalance.cxx:461
 AliBalance.cxx:462
 AliBalance.cxx:463
 AliBalance.cxx:464
 AliBalance.cxx:465
 AliBalance.cxx:466
 AliBalance.cxx:467
 AliBalance.cxx:468
 AliBalance.cxx:469
 AliBalance.cxx:470
 AliBalance.cxx:471
 AliBalance.cxx:472
 AliBalance.cxx:473
 AliBalance.cxx:474
 AliBalance.cxx:475
 AliBalance.cxx:476
 AliBalance.cxx:477
 AliBalance.cxx:478
 AliBalance.cxx:479
 AliBalance.cxx:480
 AliBalance.cxx:481
 AliBalance.cxx:482
 AliBalance.cxx:483
 AliBalance.cxx:484
 AliBalance.cxx:485
 AliBalance.cxx:486
 AliBalance.cxx:487
 AliBalance.cxx:488
 AliBalance.cxx:489
 AliBalance.cxx:490
 AliBalance.cxx:491
 AliBalance.cxx:492
 AliBalance.cxx:493
 AliBalance.cxx:494
 AliBalance.cxx:495
 AliBalance.cxx:496
 AliBalance.cxx:497
 AliBalance.cxx:498
 AliBalance.cxx:499
 AliBalance.cxx:500
 AliBalance.cxx:501
 AliBalance.cxx:502
 AliBalance.cxx:503
 AliBalance.cxx:504
 AliBalance.cxx:505
 AliBalance.cxx:506
 AliBalance.cxx:507
 AliBalance.cxx:508
 AliBalance.cxx:509
 AliBalance.cxx:510
 AliBalance.cxx:511
 AliBalance.cxx:512
 AliBalance.cxx:513
 AliBalance.cxx:514
 AliBalance.cxx:515
 AliBalance.cxx:516
 AliBalance.cxx:517
 AliBalance.cxx:518
 AliBalance.cxx:519
 AliBalance.cxx:520
 AliBalance.cxx:521
 AliBalance.cxx:522
 AliBalance.cxx:523
 AliBalance.cxx:524
 AliBalance.cxx:525
 AliBalance.cxx:526
 AliBalance.cxx:527
 AliBalance.cxx:528
 AliBalance.cxx:529
 AliBalance.cxx:530
 AliBalance.cxx:531
 AliBalance.cxx:532
 AliBalance.cxx:533
 AliBalance.cxx:534
 AliBalance.cxx:535
 AliBalance.cxx:536
 AliBalance.cxx:537
 AliBalance.cxx:538
 AliBalance.cxx:539
 AliBalance.cxx:540
 AliBalance.cxx:541
 AliBalance.cxx:542
 AliBalance.cxx:543
 AliBalance.cxx:544
 AliBalance.cxx:545
 AliBalance.cxx:546
 AliBalance.cxx:547
 AliBalance.cxx:548
 AliBalance.cxx:549
 AliBalance.cxx:550
 AliBalance.cxx:551
 AliBalance.cxx:552
 AliBalance.cxx:553
 AliBalance.cxx:554
 AliBalance.cxx:555
 AliBalance.cxx:556
 AliBalance.cxx:557
 AliBalance.cxx:558
 AliBalance.cxx:559
 AliBalance.cxx:560
 AliBalance.cxx:561
 AliBalance.cxx:562
 AliBalance.cxx:563
 AliBalance.cxx:564
 AliBalance.cxx:565
 AliBalance.cxx:566
 AliBalance.cxx:567
 AliBalance.cxx:568
 AliBalance.cxx:569
 AliBalance.cxx:570
 AliBalance.cxx:571
 AliBalance.cxx:572
 AliBalance.cxx:573
 AliBalance.cxx:574
 AliBalance.cxx:575
 AliBalance.cxx:576
 AliBalance.cxx:577
 AliBalance.cxx:578
 AliBalance.cxx:579
 AliBalance.cxx:580
 AliBalance.cxx:581
 AliBalance.cxx:582
 AliBalance.cxx:583
 AliBalance.cxx:584
 AliBalance.cxx:585
 AliBalance.cxx:586
 AliBalance.cxx:587
 AliBalance.cxx:588
 AliBalance.cxx:589
 AliBalance.cxx:590
 AliBalance.cxx:591
 AliBalance.cxx:592
 AliBalance.cxx:593
 AliBalance.cxx:594
 AliBalance.cxx:595
 AliBalance.cxx:596
 AliBalance.cxx:597
 AliBalance.cxx:598
 AliBalance.cxx:599
 AliBalance.cxx:600
 AliBalance.cxx:601
 AliBalance.cxx:602
 AliBalance.cxx:603
 AliBalance.cxx:604
 AliBalance.cxx:605
 AliBalance.cxx:606
 AliBalance.cxx:607
 AliBalance.cxx:608
 AliBalance.cxx:609
 AliBalance.cxx:610
 AliBalance.cxx:611
 AliBalance.cxx:612
 AliBalance.cxx:613
 AliBalance.cxx:614
 AliBalance.cxx:615
 AliBalance.cxx:616
 AliBalance.cxx:617
 AliBalance.cxx:618
 AliBalance.cxx:619
 AliBalance.cxx:620
 AliBalance.cxx:621
 AliBalance.cxx:622
 AliBalance.cxx:623
 AliBalance.cxx:624
 AliBalance.cxx:625
 AliBalance.cxx:626
 AliBalance.cxx:627
 AliBalance.cxx:628
 AliBalance.cxx:629
 AliBalance.cxx:630
 AliBalance.cxx:631
 AliBalance.cxx:632
 AliBalance.cxx:633
 AliBalance.cxx:634
 AliBalance.cxx:635
 AliBalance.cxx:636
 AliBalance.cxx:637
 AliBalance.cxx:638
 AliBalance.cxx:639
 AliBalance.cxx:640
 AliBalance.cxx:641
 AliBalance.cxx:642
 AliBalance.cxx:643
 AliBalance.cxx:644
 AliBalance.cxx:645
 AliBalance.cxx:646
 AliBalance.cxx:647
 AliBalance.cxx:648
 AliBalance.cxx:649
 AliBalance.cxx:650
 AliBalance.cxx:651
 AliBalance.cxx:652
 AliBalance.cxx:653
 AliBalance.cxx:654
 AliBalance.cxx:655
 AliBalance.cxx:656
 AliBalance.cxx:657
 AliBalance.cxx:658
 AliBalance.cxx:659
 AliBalance.cxx:660
 AliBalance.cxx:661
 AliBalance.cxx:662
 AliBalance.cxx:663
 AliBalance.cxx:664
 AliBalance.cxx:665
 AliBalance.cxx:666
 AliBalance.cxx:667
 AliBalance.cxx:668
 AliBalance.cxx:669
 AliBalance.cxx:670
 AliBalance.cxx:671
 AliBalance.cxx:672
 AliBalance.cxx:673
 AliBalance.cxx:674
 AliBalance.cxx:675
 AliBalance.cxx:676
 AliBalance.cxx:677
 AliBalance.cxx:678
 AliBalance.cxx:679
 AliBalance.cxx:680
 AliBalance.cxx:681
 AliBalance.cxx:682
 AliBalance.cxx:683
 AliBalance.cxx:684
 AliBalance.cxx:685
 AliBalance.cxx:686
 AliBalance.cxx:687
 AliBalance.cxx:688
 AliBalance.cxx:689
 AliBalance.cxx:690
 AliBalance.cxx:691
 AliBalance.cxx:692
 AliBalance.cxx:693
 AliBalance.cxx:694
 AliBalance.cxx:695
 AliBalance.cxx:696
 AliBalance.cxx:697
 AliBalance.cxx:698
 AliBalance.cxx:699
 AliBalance.cxx:700
 AliBalance.cxx:701
 AliBalance.cxx:702
 AliBalance.cxx:703
 AliBalance.cxx:704
 AliBalance.cxx:705
 AliBalance.cxx:706
 AliBalance.cxx:707
 AliBalance.cxx:708
 AliBalance.cxx:709
 AliBalance.cxx:710
 AliBalance.cxx:711
 AliBalance.cxx:712
 AliBalance.cxx:713
 AliBalance.cxx:714
 AliBalance.cxx:715
 AliBalance.cxx:716
 AliBalance.cxx:717
 AliBalance.cxx:718
 AliBalance.cxx:719
 AliBalance.cxx:720
 AliBalance.cxx:721
 AliBalance.cxx:722
 AliBalance.cxx:723
 AliBalance.cxx:724
 AliBalance.cxx:725
 AliBalance.cxx:726
 AliBalance.cxx:727
 AliBalance.cxx:728
 AliBalance.cxx:729
 AliBalance.cxx:730
 AliBalance.cxx:731
 AliBalance.cxx:732
 AliBalance.cxx:733
 AliBalance.cxx:734
 AliBalance.cxx:735
 AliBalance.cxx:736
 AliBalance.cxx:737
 AliBalance.cxx:738
 AliBalance.cxx:739
 AliBalance.cxx:740
 AliBalance.cxx:741
 AliBalance.cxx:742
 AliBalance.cxx:743
 AliBalance.cxx:744
 AliBalance.cxx:745
 AliBalance.cxx:746
 AliBalance.cxx:747
 AliBalance.cxx:748
 AliBalance.cxx:749
 AliBalance.cxx:750
 AliBalance.cxx:751
 AliBalance.cxx:752
 AliBalance.cxx:753
 AliBalance.cxx:754
 AliBalance.cxx:755
 AliBalance.cxx:756
 AliBalance.cxx:757
 AliBalance.cxx:758
 AliBalance.cxx:759
 AliBalance.cxx:760
 AliBalance.cxx:761
 AliBalance.cxx:762
 AliBalance.cxx:763
 AliBalance.cxx:764
 AliBalance.cxx:765
 AliBalance.cxx:766
 AliBalance.cxx:767
 AliBalance.cxx:768
 AliBalance.cxx:769
 AliBalance.cxx:770
 AliBalance.cxx:771
 AliBalance.cxx:772
 AliBalance.cxx:773
 AliBalance.cxx:774
 AliBalance.cxx:775
 AliBalance.cxx:776
 AliBalance.cxx:777
 AliBalance.cxx:778
 AliBalance.cxx:779
 AliBalance.cxx:780
 AliBalance.cxx:781
 AliBalance.cxx:782
 AliBalance.cxx:783
 AliBalance.cxx:784
 AliBalance.cxx:785
 AliBalance.cxx:786
 AliBalance.cxx:787
 AliBalance.cxx:788
 AliBalance.cxx:789
 AliBalance.cxx:790
 AliBalance.cxx:791
 AliBalance.cxx:792
 AliBalance.cxx:793
 AliBalance.cxx:794
 AliBalance.cxx:795
 AliBalance.cxx:796
 AliBalance.cxx:797
 AliBalance.cxx:798
 AliBalance.cxx:799
 AliBalance.cxx:800
 AliBalance.cxx:801
 AliBalance.cxx:802
 AliBalance.cxx:803
 AliBalance.cxx:804
 AliBalance.cxx:805
 AliBalance.cxx:806
 AliBalance.cxx:807
 AliBalance.cxx:808
 AliBalance.cxx:809
 AliBalance.cxx:810
 AliBalance.cxx:811
 AliBalance.cxx:812
 AliBalance.cxx:813
 AliBalance.cxx:814
 AliBalance.cxx:815
 AliBalance.cxx:816
 AliBalance.cxx:817
 AliBalance.cxx:818
 AliBalance.cxx:819
 AliBalance.cxx:820
 AliBalance.cxx:821
 AliBalance.cxx:822
 AliBalance.cxx:823
 AliBalance.cxx:824
 AliBalance.cxx:825
 AliBalance.cxx:826
 AliBalance.cxx:827
 AliBalance.cxx:828
 AliBalance.cxx:829
 AliBalance.cxx:830
 AliBalance.cxx:831
 AliBalance.cxx:832
 AliBalance.cxx:833
 AliBalance.cxx:834
 AliBalance.cxx:835
 AliBalance.cxx:836
 AliBalance.cxx:837
 AliBalance.cxx:838
 AliBalance.cxx:839
 AliBalance.cxx:840
 AliBalance.cxx:841
 AliBalance.cxx:842
 AliBalance.cxx:843
 AliBalance.cxx:844
 AliBalance.cxx:845
 AliBalance.cxx:846
 AliBalance.cxx:847
 AliBalance.cxx:848
 AliBalance.cxx:849
 AliBalance.cxx:850
 AliBalance.cxx:851
 AliBalance.cxx:852
 AliBalance.cxx:853
 AliBalance.cxx:854
 AliBalance.cxx:855
 AliBalance.cxx:856
 AliBalance.cxx:857
 AliBalance.cxx:858
 AliBalance.cxx:859
 AliBalance.cxx:860
 AliBalance.cxx:861
 AliBalance.cxx:862
 AliBalance.cxx:863
 AliBalance.cxx:864
 AliBalance.cxx:865
 AliBalance.cxx:866
 AliBalance.cxx:867
 AliBalance.cxx:868
 AliBalance.cxx:869
 AliBalance.cxx:870
 AliBalance.cxx:871
 AliBalance.cxx:872
 AliBalance.cxx:873
 AliBalance.cxx:874
 AliBalance.cxx:875
 AliBalance.cxx:876
 AliBalance.cxx:877
 AliBalance.cxx:878
 AliBalance.cxx:879
 AliBalance.cxx:880
 AliBalance.cxx:881
 AliBalance.cxx:882
 AliBalance.cxx:883
 AliBalance.cxx:884
 AliBalance.cxx:885
 AliBalance.cxx:886
 AliBalance.cxx:887
 AliBalance.cxx:888
 AliBalance.cxx:889
 AliBalance.cxx:890
 AliBalance.cxx:891
 AliBalance.cxx:892
 AliBalance.cxx:893
 AliBalance.cxx:894
 AliBalance.cxx:895
 AliBalance.cxx:896
 AliBalance.cxx:897
 AliBalance.cxx:898
 AliBalance.cxx:899
 AliBalance.cxx:900
 AliBalance.cxx:901
 AliBalance.cxx:902
 AliBalance.cxx:903
 AliBalance.cxx:904
 AliBalance.cxx:905
 AliBalance.cxx:906
 AliBalance.cxx:907
 AliBalance.cxx:908
 AliBalance.cxx:909
 AliBalance.cxx:910
 AliBalance.cxx:911
 AliBalance.cxx:912
 AliBalance.cxx:913
 AliBalance.cxx:914
 AliBalance.cxx:915
 AliBalance.cxx:916
 AliBalance.cxx:917
 AliBalance.cxx:918
 AliBalance.cxx:919
 AliBalance.cxx:920
 AliBalance.cxx:921
 AliBalance.cxx:922
 AliBalance.cxx:923
 AliBalance.cxx:924
 AliBalance.cxx:925
 AliBalance.cxx:926
 AliBalance.cxx:927
 AliBalance.cxx:928
 AliBalance.cxx:929
 AliBalance.cxx:930
 AliBalance.cxx:931
 AliBalance.cxx:932
 AliBalance.cxx:933
 AliBalance.cxx:934
 AliBalance.cxx:935
 AliBalance.cxx:936
 AliBalance.cxx:937
 AliBalance.cxx:938
 AliBalance.cxx:939
 AliBalance.cxx:940
 AliBalance.cxx:941
 AliBalance.cxx:942
 AliBalance.cxx:943
 AliBalance.cxx:944
 AliBalance.cxx:945
 AliBalance.cxx:946
 AliBalance.cxx:947
 AliBalance.cxx:948
 AliBalance.cxx:949
 AliBalance.cxx:950
 AliBalance.cxx:951
 AliBalance.cxx:952
 AliBalance.cxx:953
 AliBalance.cxx:954
 AliBalance.cxx:955
 AliBalance.cxx:956
 AliBalance.cxx:957
 AliBalance.cxx:958
 AliBalance.cxx:959
 AliBalance.cxx:960
 AliBalance.cxx:961
 AliBalance.cxx:962
 AliBalance.cxx:963
 AliBalance.cxx:964
 AliBalance.cxx:965
 AliBalance.cxx:966
 AliBalance.cxx:967
 AliBalance.cxx:968
 AliBalance.cxx:969
 AliBalance.cxx:970
 AliBalance.cxx:971
 AliBalance.cxx:972
 AliBalance.cxx:973
 AliBalance.cxx:974
 AliBalance.cxx:975
 AliBalance.cxx:976
 AliBalance.cxx:977
 AliBalance.cxx:978
 AliBalance.cxx:979
 AliBalance.cxx:980
 AliBalance.cxx:981
 AliBalance.cxx:982
 AliBalance.cxx:983
 AliBalance.cxx:984
 AliBalance.cxx:985
 AliBalance.cxx:986
 AliBalance.cxx:987
 AliBalance.cxx:988
 AliBalance.cxx:989
 AliBalance.cxx:990
 AliBalance.cxx:991
 AliBalance.cxx:992
 AliBalance.cxx:993
 AliBalance.cxx:994
 AliBalance.cxx:995
 AliBalance.cxx:996
 AliBalance.cxx:997
 AliBalance.cxx:998
 AliBalance.cxx:999
 AliBalance.cxx:1000
 AliBalance.cxx:1001
 AliBalance.cxx:1002
 AliBalance.cxx:1003
 AliBalance.cxx:1004
 AliBalance.cxx:1005
 AliBalance.cxx:1006
 AliBalance.cxx:1007
 AliBalance.cxx:1008
 AliBalance.cxx:1009
 AliBalance.cxx:1010
 AliBalance.cxx:1011
 AliBalance.cxx:1012
 AliBalance.cxx:1013
 AliBalance.cxx:1014
 AliBalance.cxx:1015
 AliBalance.cxx:1016
 AliBalance.cxx:1017
 AliBalance.cxx:1018
 AliBalance.cxx:1019
 AliBalance.cxx:1020
 AliBalance.cxx:1021
 AliBalance.cxx:1022
 AliBalance.cxx:1023
 AliBalance.cxx:1024
 AliBalance.cxx:1025
 AliBalance.cxx:1026
 AliBalance.cxx:1027
 AliBalance.cxx:1028
 AliBalance.cxx:1029
 AliBalance.cxx:1030
 AliBalance.cxx:1031
 AliBalance.cxx:1032
 AliBalance.cxx:1033
 AliBalance.cxx:1034
 AliBalance.cxx:1035
 AliBalance.cxx:1036
 AliBalance.cxx:1037
 AliBalance.cxx:1038
 AliBalance.cxx:1039
 AliBalance.cxx:1040
 AliBalance.cxx:1041
 AliBalance.cxx:1042
 AliBalance.cxx:1043
 AliBalance.cxx:1044
 AliBalance.cxx:1045
 AliBalance.cxx:1046
 AliBalance.cxx:1047
 AliBalance.cxx:1048
 AliBalance.cxx:1049
 AliBalance.cxx:1050
 AliBalance.cxx:1051
 AliBalance.cxx:1052
 AliBalance.cxx:1053
 AliBalance.cxx:1054
 AliBalance.cxx:1055
 AliBalance.cxx:1056
 AliBalance.cxx:1057
 AliBalance.cxx:1058
 AliBalance.cxx:1059
 AliBalance.cxx:1060
 AliBalance.cxx:1061
 AliBalance.cxx:1062
 AliBalance.cxx:1063
 AliBalance.cxx:1064
 AliBalance.cxx:1065
 AliBalance.cxx:1066
 AliBalance.cxx:1067
 AliBalance.cxx:1068
 AliBalance.cxx:1069
 AliBalance.cxx:1070
 AliBalance.cxx:1071
 AliBalance.cxx:1072
 AliBalance.cxx:1073
 AliBalance.cxx:1074
 AliBalance.cxx:1075
 AliBalance.cxx:1076
 AliBalance.cxx:1077
 AliBalance.cxx:1078
 AliBalance.cxx:1079
 AliBalance.cxx:1080
 AliBalance.cxx:1081
 AliBalance.cxx:1082
 AliBalance.cxx:1083
 AliBalance.cxx:1084
 AliBalance.cxx:1085
 AliBalance.cxx:1086
 AliBalance.cxx:1087
 AliBalance.cxx:1088
 AliBalance.cxx:1089
 AliBalance.cxx:1090
 AliBalance.cxx:1091
 AliBalance.cxx:1092
 AliBalance.cxx:1093
 AliBalance.cxx:1094
 AliBalance.cxx:1095
 AliBalance.cxx:1096
 AliBalance.cxx:1097
 AliBalance.cxx:1098
 AliBalance.cxx:1099
 AliBalance.cxx:1100
 AliBalance.cxx:1101
 AliBalance.cxx:1102
 AliBalance.cxx:1103
 AliBalance.cxx:1104
 AliBalance.cxx:1105