[u/mrichter/AliRoot.git] / PWG2 / EBYE / AliBalance.cxx

/**************************************************************************
 * 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 <TLorentzVector.h>
#include <TGraphErrors.h>

#include "AliBalance.h"

ClassImp(AliBalance)

//____________________________________________________________________//
AliBalance::AliBalance() :
  TObject(), 
  fCharge(0), fNtrack(0), fV(0), fNumberOfBins(0),
  fAnalysisType(0), fAnalyzedEvents(0), fP2Start(0),
  fP2Stop(0), fP2Step(0), fNn(0), fNp(0) {
  // Default constructor
  for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++) {
    fNpp[i] = .0;
    fNnn[i] = .0;
    fNpn[i] = .0;
    fB[i] = 0.0;
    ferror[i] = 0.0;
  } 
}

//____________________________________________________________________//
AliBalance::AliBalance(Double_t p2Start, Double_t p2Stop, Int_t p2Bins) :
  TObject(), 
  fCharge(0), fNtrack(0), fV(0), 
  fNumberOfBins(p2Bins), fAnalysisType(0), 
  fAnalyzedEvents(0), fP2Start(p2Start), fP2Stop(p2Stop), 
  fP2Step(TMath::Abs(fP2Start - fP2Stop) / (Double_t)fNumberOfBins), 
  fNn(0), fNp(0) {
  // Constructor
  for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++) {
    fNpp[i] = .0;
    fNnn[i] = .0;
    fNpn[i] = .0;
    fB[i] = 0.0;
    ferror[i] = 0.0;
  } 
}

//____________________________________________________________________//
AliBalance::AliBalance(const AliBalance& balance):
  TObject(balance),
  fNtrack(balance.fNtrack),
  fNumberOfBins(balance.fNumberOfBins),
  fAnalysisType(balance.fAnalysisType),
  fAnalyzedEvents(balance.fAnalyzedEvents),
  fP2Start(balance.fP2Start),
  fP2Stop(balance.fP2Stop),
  fP2Step(balance.fP2Step),
  fNn(balance.fNn),
  fNp(balance.fNp),
  fCharge(0),
  fV(0)
{
  //copy constructor

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

  if (balance.fV) fV = new TLorentzVector(*(balance.fV));
  if (balance.fCharge) fCharge = new Double_t(*(balance.fCharge));
}

//____________________________________________________________________//
AliBalance::~AliBalance() {
  // Destructor
  delete fV;
  delete fCharge;
}

//____________________________________________________________________//
void AliBalance::SetNnn(Double_t *nn) {
  // Setter of the Nnn term
  for(Int_t i = 0; i < fNumberOfBins; i++) fNnn[i] = nn[i];
}

//____________________________________________________________________//
void AliBalance::SetNpp(Double_t *pp) {
  // Setter of the Npp term
  for(Int_t i = 0; i < fNumberOfBins; i++) fNpp[i] = pp[i];
}

//____________________________________________________________________//
void AliBalance::SetNpn(Double_t *pn) {
  // Setter of the Npn term
  for(Int_t i = 0; i < fNumberOfBins; i++) fNpn[i] = pn[i];
}

//____________________________________________________________________//
void AliBalance::SetNumberOfBins(Int_t ibins) {
  // Sets the number of bins for the analyzed interval
  fNumberOfBins = ibins;
}

//____________________________________________________________________//
void AliBalance::SetInterval(Double_t p2Start, Double_t p2Stop) {
  // Sets the analyzed interval. 
  // The analysis variable is set by SetAnalysisType
  fP2Start = p2Start;
  fP2Stop = p2Stop;
  fP2Step = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins;
}

//____________________________________________________________________//
void AliBalance::SetAnalysisType(Int_t iType) {
  //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
  this->fAnalysisType = iType; 
  if(fAnalysisType==0) {
    cout<<" ====================== "<<endl;
    cout<<"||Analysis selected: y||"<<endl;
    cout<<" ====================== "<<endl;
  } 
  else if(fAnalysisType==1) {
    cout<<" ======================== "<<endl;
    cout<<"||Analysis selected: eta||"<<endl;
    cout<<" ======================== "<<endl;
  }
  else if(fAnalysisType==2) {
    cout<<" ========================== "<<endl;
    cout<<"||Analysis selected: Qlong||"<<endl;
    cout<<" ========================== "<<endl;
  }
  else if(fAnalysisType==3) {
    cout<<" ========================= "<<endl;
    cout<<"||Analysis selected: Qout||"<<endl;
    cout<<" ========================= "<<endl;
  }
  else if(fAnalysisType==4) {
    cout<<" ========================== "<<endl;
    cout<<"||Analysis selected: Qside||"<<endl;
    cout<<" ========================== "<<endl;
  }
  else if(fAnalysisType==5) {
    cout<<" ========================= "<<endl;
    cout<<"||Analysis selected: Qinv||"<<endl;
    cout<<" ========================= "<<endl;
  }
  else if(fAnalysisType==6) {
    cout<<" ======================== "<<endl;
    cout<<"||Analysis selected: phi||"<<endl;
    cout<<" ======================== "<<endl;
  }
  else {
    cout<<"Selection of analysis mode failed!!!"<<endl;
    cout<<"Choices are: 0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi"<<endl;
    abort();
  }
}

//____________________________________________________________________//
void AliBalance::SetParticles(TLorentzVector *P, Double_t *charge, Int_t dim) {
  // Sets a new particle with given 4-momentum and charge.
  // dim is the size of the array of charges and corresponds
  // to the number of selected tracks.
  this->fV = P;
  this->fCharge = charge;
  fNtrack = dim;
}


//____________________________________________________________________//
void AliBalance::CalculateBalance() {
  // Calculates the balance function
  fAnalyzedEvents++;
  Int_t i = 0 , j = 0;
  Int_t ibin = 0;

  for(i = 0; i < fNtrack; i++) {
    if(fCharge[i] > 0) fNp += 1.;
    if(fCharge[i] < 0) fNn += 1.;
  }

  //0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
  if(fAnalysisType==0) {
    for(i = 1; i < fNtrack; i++) {
      for(j = 0; j < i; j++) {
	Double_t rap1 = 0.5*log((fV[i].E() + fV[i].Pz())/(fV[i].E() - fV[i].Pz())); 
	Double_t rap2 = 0.5*log((fV[j].E() + fV[j].Pz())/(fV[j].E() - fV[j].Pz())); 
	Double_t dy = TMath::Abs(rap1 - rap2);
	ibin = Int_t(dy/fP2Step);
	if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
	if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
	if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
	if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
      }
    }
  }//case 0
  if(fAnalysisType==1) {
    for(i = 1; i < fNtrack; i++) {
      for(j = 0; j < i; j++) {
	Double_t p1 = sqrt(pow(fV[i].Px(),2) + pow(fV[i].Py(),2) + pow(fV[i].Pz(),2)); 
	Double_t p2 = sqrt(pow(fV[j].Px(),2) + pow(fV[j].Py(),2) + pow(fV[j].Pz(),2));
	Double_t eta1 = 0.5*log((p1 + fV[i].Pz())/(p1 - fV[i].Pz())); 
	Double_t eta2 = 0.5*log((p2 + fV[j].Pz())/(p2 - fV[j].Pz())); 
	Double_t deta = TMath::Abs(eta1 - eta2);
	ibin = Int_t(deta/fP2Step);
	if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
	if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
	if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
	if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
      }
    }
  }//case 1
  if(fAnalysisType==2) {
    for(i = 1; i < fNtrack; i++) {
      for(j = 0; j < i; j++) {
	Double_t eTot = fV[i].E() + fV[j].E();
	Double_t pxTot = fV[i].Px() + fV[j].Px();
	Double_t pyTot = fV[i].Py() + fV[j].Py();
	Double_t pzTot = fV[i].Pz() + fV[j].Pz();
	Double_t q0Tot = fV[i].E() - fV[j].E();
	Double_t qzTot = fV[i].Pz() - fV[j].Pz();
	Double_t snn = pow(eTot,2) - pow(pxTot,2) - pow(pyTot,2) - pow(pzTot,2);
	Double_t ptTot = sqrt( pow(pxTot,2) + pow(pyTot,2));
	Double_t qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/sqrt(snn + pow(ptTot,2));
	ibin = Int_t(qLong/fP2Step);
	//cout<<ibin<<endl;
	if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
	if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
	if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
	if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
      }
    }
  }//case 2
  if(fAnalysisType==3) {
    for(i = 1; i < fNtrack; i++) {
      for(j = 0; j < i; j++) {
	Double_t eTot = fV[i].E() + fV[j].E();
	Double_t pxTot = fV[i].Px() + fV[j].Px();
	Double_t pyTot = fV[i].Py() + fV[j].Py();
	Double_t pzTot = fV[i].Pz() + fV[j].Pz();
	Double_t qxTot = fV[i].Px() - fV[j].Px();
	Double_t qyTot = fV[i].Py() - fV[j].Py();
	Double_t snn = pow(eTot,2) - pow(pxTot,2) - pow(pyTot,2) - pow(pzTot,2);
	Double_t ptTot = sqrt( pow(pxTot,2) + pow(pyTot,2));
	Double_t qOut = sqrt(snn/(snn + pow(ptTot,2))) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
	ibin = Int_t(qOut/fP2Step);
	//cout<<ibin<<endl;
	if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
	if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
	if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
	if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
      }
    }
  }//case 3
  if(fAnalysisType==4) {
    for(i = 1; i < fNtrack; i++) {
      for(j = 0; j < i; j++) {
	Double_t pxTot = fV[i].Px() + fV[j].Px();
	Double_t pyTot = fV[i].Py() + fV[j].Py();
	Double_t qxTot = fV[i].Px() - fV[j].Px();
	Double_t qyTot = fV[i].Py() - fV[j].Py();
	Double_t ptTot = sqrt( pow(pxTot,2) + pow(pyTot,2));
	Double_t qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
	ibin = Int_t(qSide/fP2Step);
	//cout<<ibin<<endl;
	if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
	if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
	if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
	if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
      }
    }
  }//case 4
  if(fAnalysisType==5) {
    for(i = 1; i < fNtrack; i++) {
      for(j = 0; j < i; j++) {
	Double_t q0Tot = fV[i].E() - fV[j].E();
	Double_t qxTot = fV[i].Px() - fV[j].Px();
	Double_t qyTot = fV[i].Py() - fV[j].Py();
	Double_t qzTot = fV[i].Pz() - fV[j].Pz();
	Double_t qInv = sqrt(TMath::Abs(-pow(q0Tot,2) +pow(qxTot,2) +pow(qyTot,2) +pow(qzTot,2)));
	ibin = Int_t(qInv/fP2Step);
	//cout<<ibin<<endl;
	if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
	if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
	if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
	if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
      }
    }
  }//case 5	
  if(fAnalysisType==6) {
    for(i = 1; i < fNtrack; i++) {
      for(j = 0; j < i; j++) {
	Double_t phi1 = TMath::ATan(fV[i].Py()/fV[i].Px())*180.0/TMath::Pi();
	Double_t phi2 = TMath::ATan(fV[j].Py()/fV[j].Px())*180.0/TMath::Pi();
	Double_t dphi = TMath::Abs(phi1 - phi2);
	ibin = Int_t(dphi/fP2Step);
	if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
	if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
	if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
	if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
      }
    }
  }//case 6
}

//____________________________________________________________________//
Double_t AliBalance::GetBalance(Int_t p2) {
  // Returns the value of the balance function in bin p2
  fB[p2] = 0.5*(((fNpn[p2] - 2.0*fNnn[p2])/fNn) + ((fNpn[p2] - 2.0*fNpp[p2])/fNp))/fP2Step;
  
  return fB[p2];
}
    
//____________________________________________________________________//
Double_t AliBalance::GetError(Int_t p2) {
  // Returns the error on the BF value for bin p2
  ferror[p2] = sqrt( Double_t(fNpp[p2])/(Double_t(fNp)*Double_t(fNp)) + Double_t(fNnn[p2])/(Double_t(fNn)*Double_t(fNn)) + Double_t(fNpn[p2])*pow((0.5/Double_t(fNp) + 0.5/Double_t(fNn)),2))/fP2Step;

  return ferror[p2];
}

//____________________________________________________________________//
void AliBalance::PrintResults() {
  // Prints the results
  Double_t x[MAXIMUM_NUMBER_OF_STEPS];
  Double_t fSumXi = 0.0, fSumBi = 0.0, fSumBiXi = 0.0;
  Double_t fSumBiXi2 = 0.0, fSumBi2Xi2 = 0.0;
  Double_t fSumDeltaBi2 = 0.0, fSumXi2DeltaBi2 = 0.0;
  Double_t deltaBalP2 = 0.0, integral = 0.0;
  Double_t deltaErrorNew = 0.0;

  cout<<"=================================================="<<endl;
  for(Int_t i = 0; i < fNumberOfBins; i++) { 
    x[i] = fP2Start + fP2Step*i + fP2Step/2;
    cout<<"B: "<<fB[i]<<"\t Error: "<<ferror[i]<<"\t bin: "<<x[i]<<endl;
  } 
  cout<<"=================================================="<<endl;
  for(Int_t i = 1; i < fNumberOfBins; i++) {
    fSumXi += x[i];
    fSumBi += fB[i];
    fSumBiXi += fB[i]*x[i];
    fSumBiXi2 += fB[i]*pow(x[i],2);
    fSumBi2Xi2 += pow(fB[i],2)*pow(x[i],2);
    fSumDeltaBi2 +=  pow(ferror[i],2);
    fSumXi2DeltaBi2 += pow(x[i],2) * pow(ferror[i],2);
    
    deltaBalP2 += fP2Step*pow(ferror[i],2);
    integral += fP2Step*fB[i];
  }
  for(Int_t i = 1; i < fNumberOfBins; i++) deltaErrorNew += ferror[i]*(x[i]*fSumBi - fSumBiXi)/pow(fSumBi,2);
   
  Double_t integralError = sqrt(deltaBalP2);
  
  Double_t delta = fSumBiXi / fSumBi;
  Double_t deltaError = (fSumBiXi / fSumBi) * sqrt(pow((sqrt(fSumXi2DeltaBi2)/fSumBiXi),2) + pow((fSumDeltaBi2/fSumBi),2) );
 
  cout<<"Analyzed events: "<<fAnalyzedEvents<<endl;
  cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
  cout<<"New error: "<<deltaErrorNew<<endl;
  cout<<"Interval: "<<integral<<"\t Error: "<<integralError<<endl;
  cout<<"=================================================="<<endl;
}
  
//____________________________________________________________________//
TGraphErrors *AliBalance::DrawBalance() {
  // 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 == 0)||(fNn == 0)) {
    cout<<"Couldn't find any particles in the analyzed interval!!!"<<endl;
    cout<<"Aborting....."<<endl;
    abort();
  }
  
  for(Int_t i = 0; i < fNumberOfBins; i++) {
    b[i] = GetBalance(i);
    ber[i] = GetError(i);
    x[i] = fP2Start + fP2Step*i + fP2Step/2;
    xer[i] = 0.0;
  }
  
  TGraphErrors *gr = new TGraphErrors(fNumberOfBins,x,b,xer,ber);
  gr->SetMarkerStyle(25);
  gr->GetXaxis()->SetTitleColor(1);
  if(fAnalysisType==0) {
    gr->GetXaxis()->SetTitle("#Delta y");
    gr->GetYaxis()->SetTitle("B(#Delta y)");
  }
  if(fAnalysisType==1) {
    gr->GetXaxis()->SetTitle("#Delta #eta");
    gr->GetYaxis()->SetTitle("B(#Delta #eta)");
  }
  if(fAnalysisType==2) {
    gr->GetXaxis()->SetTitle("Q_{long} [GeV]");
    gr->GetYaxis()->SetTitle("B(Q_{long})");
  }
  if(fAnalysisType==3) {
    gr->GetXaxis()->SetTitle("Q_{out} [GeV]");
    gr->GetYaxis()->SetTitle("B(Q_{out})");
  }
  if(fAnalysisType==4) {
    gr->GetXaxis()->SetTitle("Q_{side} [GeV]");
    gr->GetYaxis()->SetTitle("B(Q_{side})");
  }
  if(fAnalysisType==5) {
    gr->GetXaxis()->SetTitle("Q_{inv} [GeV]");
    gr->GetYaxis()->SetTitle("B(Q_{inv})");
  }
  if(fAnalysisType==6) {
    gr->GetXaxis()->SetTitle("#Delta #phi");
    gr->GetYaxis()->SetTitle("B(#Delta #phi)");
  }

  return gr;
}
Commit	Line	Data
6c178944	1	/**************************************************************************
	2	* Author: Panos Christakoglou. *
	3	* Contributors are mentioned in the code where appropriate. *
	4	* *
	5	* Permission to use, copy, modify and distribute this software and its *
	6	* documentation strictly for non-commercial purposes is hereby granted *
	7	* without fee, provided that the above copyright notice appears in all *
	8	* copies and that both the copyright notice and this permission notice *
	9	* appear in the supporting documentation. The authors make no claims *
	10	* about the suitability of this software for any purpose. It is *
	11	* provided "as is" without express or implied warranty. *
	12	**************************************************************************/
	13
	14	/* $Id$ */
	15
	16	//-----------------------------------------------------------------
	17	// Balance Function class
	18	// This is the class to deal with the Balance Function analysis
	19	// Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
	20	//-----------------------------------------------------------------
	21
	22
	23	//ROOT
	24	#include <Riostream.h>
	25	#include <TMath.h>
	26	#include <TLorentzVector.h>
	27	#include <TGraphErrors.h>
	28
	29	#include "AliBalance.h"
	30
	31	ClassImp(AliBalance)
	32
7f0257ea	33	//____________________________________________________________________//
5d534bd4	34	AliBalance::AliBalance() :
	35	TObject(),
	36	fCharge(0), fNtrack(0), fV(0), fNumberOfBins(0),
	37	fAnalysisType(0), fAnalyzedEvents(0), fP2Start(0),
	38	fP2Stop(0), fP2Step(0), fNn(0), fNp(0) {
6c178944	39	// Default constructor
7f0257ea	40	for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++) {
	41	fNpp[i] = .0;
	42	fNnn[i] = .0;
	43	fNpn[i] = .0;
	44	fB[i] = 0.0;
	45	ferror[i] = 0.0;
	46	}
6c178944	47	}
6c178944	48
7f0257ea	49	//____________________________________________________________________//
5d534bd4	50	AliBalance::AliBalance(Double_t p2Start, Double_t p2Stop, Int_t p2Bins) :
	51	TObject(),
	52	fCharge(0), fNtrack(0), fV(0),
	53	fNumberOfBins(p2Bins), fAnalysisType(0),
	54	fAnalyzedEvents(0), fP2Start(p2Start), fP2Stop(p2Stop),
	55	fP2Step(TMath::Abs(fP2Start - fP2Stop) / (Double_t)fNumberOfBins),
	56	fNn(0), fNp(0) {
6c178944	57	// Constructor
7f0257ea	58	for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++) {
	59	fNpp[i] = .0;
	60	fNnn[i] = .0;
	61	fNpn[i] = .0;
	62	fB[i] = 0.0;
	63	ferror[i] = 0.0;
	64	}
6c178944	65	}
6c178944	66
7f0257ea	67	//____________________________________________________________________//
7f0257ea	68	AliBalance::AliBalance(const AliBalance& balance):
e559055a	69	TObject(balance),
7f0257ea	70	fNtrack(balance.fNtrack),
	71	fNumberOfBins(balance.fNumberOfBins),
	72	fAnalysisType(balance.fAnalysisType),
	73	fAnalyzedEvents(balance.fAnalyzedEvents),
	74	fP2Start(balance.fP2Start),
	75	fP2Stop(balance.fP2Stop),
	76	fP2Step(balance.fP2Step),
	77	fNn(balance.fNn),
	78	fNp(balance.fNp),
	79	fCharge(0),
	80	fV(0)
6c178944	81	{
7f0257ea	82	//copy constructor
	83
	84	for(Int_t i = 0; i < MAXIMUM_NUMBER_OF_STEPS; i++) {
	85	fNpp[i] = .0;
	86	fNnn[i] = .0;
	87	fNpn[i] = .0;
	88	fB[i] = 0.0;
	89	ferror[i] = 0.0;
	90	}
	91
	92	if (balance.fV) fV = new TLorentzVector(*(balance.fV));
	93	if (balance.fCharge) fCharge = new Double_t(*(balance.fCharge));
	94	}
	95
	96	//____________________________________________________________________//
	97	AliBalance::~AliBalance() {
6c178944	98	// Destructor
7f0257ea	99	delete fV;
7f0257ea	100	delete fCharge;
6c178944	101	}
6c178944	102
7f0257ea	103	//____________________________________________________________________//
	104	void AliBalance::SetNnn(Double_t *nn) {
	105	// Setter of the Nnn term
	106	for(Int_t i = 0; i < fNumberOfBins; i++) fNnn[i] = nn[i];
	107	}
	108
	109	//____________________________________________________________________//
	110	void AliBalance::SetNpp(Double_t *pp) {
	111	// Setter of the Npp term
	112	for(Int_t i = 0; i < fNumberOfBins; i++) fNpp[i] = pp[i];
	113	}
	114
	115	//____________________________________________________________________//
	116	void AliBalance::SetNpn(Double_t *pn) {
	117	// Setter of the Npn term
	118	for(Int_t i = 0; i < fNumberOfBins; i++) fNpn[i] = pn[i];
	119	}
	120
	121	//____________________________________________________________________//
	122	void AliBalance::SetNumberOfBins(Int_t ibins) {
6c178944	123	// Sets the number of bins for the analyzed interval
7f0257ea	124	fNumberOfBins = ibins;
6c178944	125	}
6c178944	126
7f0257ea	127	//____________________________________________________________________//
7f0257ea	128	void AliBalance::SetInterval(Double_t p2Start, Double_t p2Stop) {
6c178944	129	// Sets the analyzed interval.
	130	// The analysis variable is set by SetAnalysisType
	131	fP2Start = p2Start;
	132	fP2Stop = p2Stop;
	133	fP2Step = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins;
	134	}
	135
7f0257ea	136	//____________________________________________________________________//
7f0257ea	137	void AliBalance::SetAnalysisType(Int_t iType) {
6c178944	138	//0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
6c178944	139	this->fAnalysisType = iType;
7f0257ea	140	if(fAnalysisType==0) {
	141	cout<<" ====================== "<<endl;
	142	cout<<"\|\|Analysis selected: y\|\|"<<endl;
	143	cout<<" ====================== "<<endl;
	144	}
	145	else if(fAnalysisType==1) {
	146	cout<<" ======================== "<<endl;
	147	cout<<"\|\|Analysis selected: eta\|\|"<<endl;
	148	cout<<" ======================== "<<endl;
	149	}
	150	else if(fAnalysisType==2) {
	151	cout<<" ========================== "<<endl;
	152	cout<<"\|\|Analysis selected: Qlong\|\|"<<endl;
	153	cout<<" ========================== "<<endl;
	154	}
	155	else if(fAnalysisType==3) {
	156	cout<<" ========================= "<<endl;
	157	cout<<"\|\|Analysis selected: Qout\|\|"<<endl;
	158	cout<<" ========================= "<<endl;
	159	}
	160	else if(fAnalysisType==4) {
	161	cout<<" ========================== "<<endl;
	162	cout<<"\|\|Analysis selected: Qside\|\|"<<endl;
	163	cout<<" ========================== "<<endl;
	164	}
	165	else if(fAnalysisType==5) {
	166	cout<<" ========================= "<<endl;
	167	cout<<"\|\|Analysis selected: Qinv\|\|"<<endl;
	168	cout<<" ========================= "<<endl;
	169	}
	170	else if(fAnalysisType==6) {
	171	cout<<" ======================== "<<endl;
	172	cout<<"\|\|Analysis selected: phi\|\|"<<endl;
	173	cout<<" ======================== "<<endl;
	174	}
	175	else {
	176	cout<<"Selection of analysis mode failed!!!"<<endl;
	177	cout<<"Choices are: 0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi"<<endl;
	178	abort();
	179	}
6c178944	180	}
6c178944	181
7f0257ea	182	//____________________________________________________________________//
7f0257ea	183	void AliBalance::SetParticles(TLorentzVector P, Double_t charge, Int_t dim) {
6c178944	184	// Sets a new particle with given 4-momentum and charge.
	185	// dim is the size of the array of charges and corresponds
	186	// to the number of selected tracks.
	187	this->fV = P;
	188	this->fCharge = charge;
	189	fNtrack = dim;
	190	}
	191
	192
7f0257ea	193	//____________________________________________________________________//
7f0257ea	194	void AliBalance::CalculateBalance() {
6c178944	195	// Calculates the balance function
	196	fAnalyzedEvents++;
	197	Int_t i = 0 , j = 0;
	198	Int_t ibin = 0;
	199
7f0257ea	200	for(i = 0; i < fNtrack; i++) {
	201	if(fCharge[i] > 0) fNp += 1.;
	202	if(fCharge[i] < 0) fNn += 1.;
	203	}
6c178944	204
6c178944	205	//0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
7f0257ea	206	if(fAnalysisType==0) {
	207	for(i = 1; i < fNtrack; i++) {
	208	for(j = 0; j < i; j++) {
	209	Double_t rap1 = 0.5*log((fV[i].E() + fV[i].Pz())/(fV[i].E() - fV[i].Pz()));
	210	Double_t rap2 = 0.5*log((fV[j].E() + fV[j].Pz())/(fV[j].E() - fV[j].Pz()));
	211	Double_t dy = TMath::Abs(rap1 - rap2);
	212	ibin = Int_t(dy/fP2Step);
	213	if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
	214	if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
	215	if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
	216	if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
	217	}
	218	}
	219	}//case 0
	220	if(fAnalysisType==1) {
	221	for(i = 1; i < fNtrack; i++) {
	222	for(j = 0; j < i; j++) {
	223	Double_t p1 = sqrt(pow(fV[i].Px(),2) + pow(fV[i].Py(),2) + pow(fV[i].Pz(),2));
	224	Double_t p2 = sqrt(pow(fV[j].Px(),2) + pow(fV[j].Py(),2) + pow(fV[j].Pz(),2));
	225	Double_t eta1 = 0.5*log((p1 + fV[i].Pz())/(p1 - fV[i].Pz()));
	226	Double_t eta2 = 0.5*log((p2 + fV[j].Pz())/(p2 - fV[j].Pz()));
	227	Double_t deta = TMath::Abs(eta1 - eta2);
	228	ibin = Int_t(deta/fP2Step);
	229	if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
	230	if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
	231	if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
	232	if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
	233	}
	234	}
	235	}//case 1
	236	if(fAnalysisType==2) {
	237	for(i = 1; i < fNtrack; i++) {
	238	for(j = 0; j < i; j++) {
	239	Double_t eTot = fV[i].E() + fV[j].E();
	240	Double_t pxTot = fV[i].Px() + fV[j].Px();
	241	Double_t pyTot = fV[i].Py() + fV[j].Py();
	242	Double_t pzTot = fV[i].Pz() + fV[j].Pz();
	243	Double_t q0Tot = fV[i].E() - fV[j].E();
	244	Double_t qzTot = fV[i].Pz() - fV[j].Pz();
	245	Double_t snn = pow(eTot,2) - pow(pxTot,2) - pow(pyTot,2) - pow(pzTot,2);
	246	Double_t ptTot = sqrt( pow(pxTot,2) + pow(pyTot,2));
	247	Double_t qLong = TMath::Abs(eTotqzTot - pzTotq0Tot)/sqrt(snn + pow(ptTot,2));
	248	ibin = Int_t(qLong/fP2Step);
	249	//cout<<ibin<<endl;
	250	if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
	251	if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
	252	if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
	253	if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
	254	}
	255	}
	256	}//case 2
	257	if(fAnalysisType==3) {
	258	for(i = 1; i < fNtrack; i++) {
	259	for(j = 0; j < i; j++) {
	260	Double_t eTot = fV[i].E() + fV[j].E();
	261	Double_t pxTot = fV[i].Px() + fV[j].Px();
	262	Double_t pyTot = fV[i].Py() + fV[j].Py();
	263	Double_t pzTot = fV[i].Pz() + fV[j].Pz();
	264	Double_t qxTot = fV[i].Px() - fV[j].Px();
	265	Double_t qyTot = fV[i].Py() - fV[j].Py();
	266	Double_t snn = pow(eTot,2) - pow(pxTot,2) - pow(pyTot,2) - pow(pzTot,2);
	267	Double_t ptTot = sqrt( pow(pxTot,2) + pow(pyTot,2));
	268	Double_t qOut = sqrt(snn/(snn + pow(ptTot,2))) * TMath::Abs(pxTotqxTot + pyTotqyTot)/ptTot;
	269	ibin = Int_t(qOut/fP2Step);
270	//cout<<ibin<<endl;
271	if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
272	if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
273	if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
274	if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
275	}
276	}
277	}//case 3
278	if(fAnalysisType==4) {
279	for(i = 1; i < fNtrack; i++) {
280	for(j = 0; j < i; j++) {
281	Double_t pxTot = fV[i].Px() + fV[j].Px();
282	Double_t pyTot = fV[i].Py() + fV[j].Py();
283	Double_t qxTot = fV[i].Px() - fV[j].Px();
284	Double_t qyTot = fV[i].Py() - fV[j].Py();
285	Double_t ptTot = sqrt( pow(pxTot,2) + pow(pyTot,2));
286	Double_t qSide = TMath::Abs(pxTotqyTot - pyTotqxTot)/ptTot;
287	ibin = Int_t(qSide/fP2Step);
288	//cout<<ibin<<endl;
289	if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
290	if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
291	if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
292	if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
293	}
294	}
295	}//case 4
296	if(fAnalysisType==5) {
297	for(i = 1; i < fNtrack; i++) {
298	for(j = 0; j < i; j++) {
299	Double_t q0Tot = fV[i].E() - fV[j].E();
300	Double_t qxTot = fV[i].Px() - fV[j].Px();
301	Double_t qyTot = fV[i].Py() - fV[j].Py();
302	Double_t qzTot = fV[i].Pz() - fV[j].Pz();
303	Double_t qInv = sqrt(TMath::Abs(-pow(q0Tot,2) +pow(qxTot,2) +pow(qyTot,2) +pow(qzTot,2)));
304	ibin = Int_t(qInv/fP2Step);
305	//cout<<ibin<<endl;
306	if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
307	if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
308	if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
309	if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
310	}
311	}
312	}//case 5
313	if(fAnalysisType==6) {
314	for(i = 1; i < fNtrack; i++) {
315	for(j = 0; j < i; j++) {
316	Double_t phi1 = TMath::ATan(fV[i].Py()/fV[i].Px())*180.0/TMath::Pi();
317	Double_t phi2 = TMath::ATan(fV[j].Py()/fV[j].Px())*180.0/TMath::Pi();
318	Double_t dphi = TMath::Abs(phi1 - phi2);
319	ibin = Int_t(dphi/fP2Step);
320	if((fCharge[i] > 0)&&(fCharge[j] > 0)) fNpp[ibin] += 1.;
321	if((fCharge[i] < 0)&&(fCharge[j] < 0)) fNnn[ibin] += 1.;
322	if((fCharge[i] > 0)&&(fCharge[j] < 0)) fNpn[ibin] += 1.;
323	if((fCharge[i] < 0)&&(fCharge[j] > 0)) fNpn[ibin] += 1.;
324	}
325	}
326	}//case 6
6c178944	327	}
6c178944	328
7f0257ea	329	//____________________________________________________________________//
7f0257ea	330	Double_t AliBalance::GetBalance(Int_t p2) {
6c178944	331	// Returns the value of the balance function in bin p2
6c178944	332	fB[p2] = 0.5(((fNpn[p2] - 2.0fNnn[p2])/fNn) + ((fNpn[p2] - 2.0*fNpp[p2])/fNp))/fP2Step;
7f0257ea	333
6c178944	334	return fB[p2];
	335	}
	336
7f0257ea	337	//____________________________________________________________________//
7f0257ea	338	Double_t AliBalance::GetError(Int_t p2) {
6c178944	339	// Returns the error on the BF value for bin p2
	340	ferror[p2] = sqrt( Double_t(fNpp[p2])/(Double_t(fNp)Double_t(fNp)) + Double_t(fNnn[p2])/(Double_t(fNn)Double_t(fNn)) + Double_t(fNpn[p2])*pow((0.5/Double_t(fNp) + 0.5/Double_t(fNn)),2))/fP2Step;
	341
	342	return ferror[p2];
	343	}
	344
7f0257ea	345	//____________________________________________________________________//
7f0257ea	346	void AliBalance::PrintResults() {
6c178944	347	// Prints the results
	348	Double_t x[MAXIMUM_NUMBER_OF_STEPS];
	349	Double_t fSumXi = 0.0, fSumBi = 0.0, fSumBiXi = 0.0;
	350	Double_t fSumBiXi2 = 0.0, fSumBi2Xi2 = 0.0;
	351	Double_t fSumDeltaBi2 = 0.0, fSumXi2DeltaBi2 = 0.0;
	352	Double_t deltaBalP2 = 0.0, integral = 0.0;
	353	Double_t deltaErrorNew = 0.0;
	354
	355	cout<<"=================================================="<<endl;
7f0257ea	356	for(Int_t i = 0; i < fNumberOfBins; i++) {
	357	x[i] = fP2Start + fP2Step*i + fP2Step/2;
	358	cout<<"B: "<<fB[i]<<"\t Error: "<<ferror[i]<<"\t bin: "<<x[i]<<endl;
	359	}
6c178944	360	cout<<"=================================================="<<endl;
7f0257ea	361	for(Int_t i = 1; i < fNumberOfBins; i++) {
	362	fSumXi += x[i];
	363	fSumBi += fB[i];
	364	fSumBiXi += fB[i]*x[i];
	365	fSumBiXi2 += fB[i]*pow(x[i],2);
	366	fSumBi2Xi2 += pow(fB[i],2)*pow(x[i],2);
	367	fSumDeltaBi2 += pow(ferror[i],2);
	368	fSumXi2DeltaBi2 += pow(x[i],2) * pow(ferror[i],2);
	369
	370	deltaBalP2 += fP2Step*pow(ferror[i],2);
	371	integral += fP2Step*fB[i];
	372	}
	373	for(Int_t i = 1; i < fNumberOfBins; i++) deltaErrorNew += ferror[i](x[i]fSumBi - fSumBiXi)/pow(fSumBi,2);
	374
	375	Double_t integralError = sqrt(deltaBalP2);
6c178944	376
7f0257ea	377	Double_t delta = fSumBiXi / fSumBi;
7f0257ea	378	Double_t deltaError = (fSumBiXi / fSumBi) * sqrt(pow((sqrt(fSumXi2DeltaBi2)/fSumBiXi),2) + pow((fSumDeltaBi2/fSumBi),2) );
6c178944	379
	380	cout<<"Analyzed events: "<<fAnalyzedEvents<<endl;
	381	cout<<"Width: "<<delta<<"\t Error: "<<deltaError<<endl;
	382	cout<<"New error: "<<deltaErrorNew<<endl;
	383	cout<<"Interval: "<<integral<<"\t Error: "<<integralError<<endl;
	384	cout<<"=================================================="<<endl;
	385	}
	386
7f0257ea	387	//____________________________________________________________________//
7f0257ea	388	TGraphErrors *AliBalance::DrawBalance() {
6c178944	389	// Draws the BF
	390	Double_t x[MAXIMUM_NUMBER_OF_STEPS];
	391	Double_t xer[MAXIMUM_NUMBER_OF_STEPS];
	392	Double_t b[MAXIMUM_NUMBER_OF_STEPS];
	393	Double_t ber[MAXIMUM_NUMBER_OF_STEPS];
	394
7f0257ea	395	if((fNp == 0)\|\|(fNn == 0)) {
	396	cout<<"Couldn't find any particles in the analyzed interval!!!"<<endl;
	397	cout<<"Aborting....."<<endl;
	398	abort();
	399	}
	400
	401	for(Int_t i = 0; i < fNumberOfBins; i++) {
	402	b[i] = GetBalance(i);
	403	ber[i] = GetError(i);
	404	x[i] = fP2Start + fP2Step*i + fP2Step/2;
	405	xer[i] = 0.0;
	406	}
	407
	408	TGraphErrors *gr = new TGraphErrors(fNumberOfBins,x,b,xer,ber);
	409	gr->SetMarkerStyle(25);
6c178944	410	gr->GetXaxis()->SetTitleColor(1);
7f0257ea	411	if(fAnalysisType==0) {
	412	gr->GetXaxis()->SetTitle("#Delta y");
	413	gr->GetYaxis()->SetTitle("B(#Delta y)");
	414	}
	415	if(fAnalysisType==1) {
	416	gr->GetXaxis()->SetTitle("#Delta #eta");
	417	gr->GetYaxis()->SetTitle("B(#Delta #eta)");
	418	}
	419	if(fAnalysisType==2) {
	420	gr->GetXaxis()->SetTitle("Q_{long} [GeV]");
	421	gr->GetYaxis()->SetTitle("B(Q_{long})");
	422	}
	423	if(fAnalysisType==3) {
	424	gr->GetXaxis()->SetTitle("Q_{out} [GeV]");
	425	gr->GetYaxis()->SetTitle("B(Q_{out})");
	426	}
	427	if(fAnalysisType==4) {
	428	gr->GetXaxis()->SetTitle("Q_{side} [GeV]");
	429	gr->GetYaxis()->SetTitle("B(Q_{side})");
	430	}
	431	if(fAnalysisType==5) {
	432	gr->GetXaxis()->SetTitle("Q_{inv} [GeV]");
	433	gr->GetYaxis()->SetTitle("B(Q_{inv})");
	434	}
	435	if(fAnalysisType==6) {
	436	gr->GetXaxis()->SetTitle("#Delta #phi");
	437	gr->GetYaxis()->SetTitle("B(#Delta #phi)");
	438	}
6c178944	439
	440	return gr;
	441	}