[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 <TAxis.h>
#include <TH1D.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"

ClassImp(AliBalance)

//____________________________________________________________________//
AliBalance::AliBalance() :
  TObject(), 
  fAnalysisLevel("ESD"),
  fAnalyzedEvents(0) {
  // Default constructor
 
  for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
    if(i == 6) {
      fNumberOfBins[i] = 180;
      fP1Start[i]      = 0.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];
    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;
      fB[i][j] = 0.0;
      ferror[i][j] = 0.0;
    } 
  }

  TString gAnalysisType[ANALYSIS_TYPES] = {"y","eta","qlong","qout","qside","qinv","phi"};
  TString histName;
  for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
    histName = "fHistP"; histName += gAnalysisType[iAnalysisType];
    fHistP[iAnalysisType] = new TH1D(histName.Data(),"",100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
    histName = "fHistN"; histName += gAnalysisType[iAnalysisType];
    fHistN[iAnalysisType] = new TH1D(histName.Data(),"",100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
  
    histName = "fHistPN"; histName += gAnalysisType[iAnalysisType];
    fHistPN[iAnalysisType] = new TH1D(histName.Data(),"",fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
    histName = "fHistNP"; histName += gAnalysisType[iAnalysisType];
    fHistNP[iAnalysisType] = new TH1D(histName.Data(),"",fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
    histName = "fHistPP"; histName += gAnalysisType[iAnalysisType];
    fHistPP[iAnalysisType] = new TH1D(histName.Data(),"",fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
    histName = "fHistNN"; histName += gAnalysisType[iAnalysisType];
    fHistNN[iAnalysisType] = new TH1D(histName.Data(),"",fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);

  }
}


//____________________________________________________________________//
AliBalance::AliBalance(const AliBalance& balance):
  TObject(balance), fAnalysisLevel(balance.fAnalysisLevel),
  fAnalyzedEvents(balance.fAnalyzedEvents) {
  //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];

    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;
      fB[i][j] = 0.0;
      ferror[i][j] = 0.0;
    } 
  }
 }
 

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

}


//____________________________________________________________________//
void AliBalance::SetNumberOfBins(Int_t ibin, Int_t ibins) {
  // Sets the number of bins for the analyzed interval
  // Set the same Information for all analyses
  if(ibin == -1){             
    for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
      fNumberOfBins[i] = ibins;
    }
  }
  // Set the Information for one analysis
  else if(ibin > -1 && ibin < ANALYSIS_TYPES){
    fNumberOfBins[ibin] = ibins;
  }
  else{
    AliError("Wrong ANALYSIS number!");
  }
}

//____________________________________________________________________//
void AliBalance::SetInterval(Double_t p1Start, Double_t p1Stop,
			     Int_t ibin, Double_t p2Start, Double_t p2Stop) {
  // Sets the analyzed interval. 
  // Set the same Information for all analyses
  if(ibin == -1){             
    for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
      fP1Start[i] = p1Start;
      fP1Stop[i] = p1Stop;
      fP2Start[i] = p2Start;
      fP2Stop[i] = p2Stop;
      fP2Step[i] = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins[i];
    }
  }
  // Set the Information for one analysis
  else if(ibin > -1 && ibin < ANALYSIS_TYPES){
      fP1Start[ibin] = p1Start;
      fP1Stop[ibin] = p1Stop;
      fP2Start[ibin] = p2Start;
      fP2Stop[ibin] = p2Stop;
      fP2Step[ibin] = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins[ibin];
  }
  else{
    AliError("Wrong ANALYSIS number!");
  }
}


//____________________________________________________________________//
void AliBalance::PrintAnalysisSettings() {
  
  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(TObjArray *gTrackArray) {
  // Calculates the balance function
  fAnalyzedEvents++;
  Int_t i = 0 , j = 0;
  Int_t iBin = 0;
  
  AliVParticle* track = 0;
  AliVParticle* track1 = 0;
  AliVParticle* track2 = 0;
    
  //Printf("(AliBalance) Number of tracks: %d",gTrackArray->GetEntries());
  Int_t gNtrack = gTrackArray->GetEntries();
  for(i = 0; i < gNtrack; i++) {
    if(fAnalysisLevel == "ESD")
      track = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
    else if(fAnalysisLevel == "AOD")
      track = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
    else if(fAnalysisLevel == "MC")
      track = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));

    if(track) {
      Short_t charge  = track->Charge();
      Double_t pseudorapidity = track->Eta();
      Double_t rapidity = track->Y();
      Double_t phi = TMath::ATan(track->Py()/track->Px())*180.0/TMath::Pi();
      
      //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 == 1) {
	  if((pseudorapidity >= fP1Start[iAnalysisType]) && (pseudorapidity <= fP1Stop[iAnalysisType])) {
	    if(charge > 0) {
	      fNp[iAnalysisType] += 1.;
	      fHistP[iAnalysisType]->Fill(pseudorapidity);
	    }//charge > 0
	    if(charge < 0) {
	      fNn[iAnalysisType] += 1.;
	      fHistN[iAnalysisType]->Fill(pseudorapidity);
	    }//charge < 0
	  }//p1 interval check
	}//analysis type: eta
	if(iAnalysisType == 6) {
	  if((phi >= fP1Start[iAnalysisType]) && (phi <= fP1Stop[iAnalysisType])) {
	    if(charge > 0) {
	      fNp[iAnalysisType] += 1.;
	      fHistP[iAnalysisType]->Fill(phi);
	    }//charge > 0
	    if(charge < 0) {
	      fNn[iAnalysisType] += 1.;
	      fHistN[iAnalysisType]->Fill(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(rapidity);
	    }//charge > 0
	    if(charge < 0) {
	      fNn[iAnalysisType] += 1.;
	      fHistN[iAnalysisType]->Fill(rapidity);
	    }//charge < 0
	  }//p1 interval check
	}//analysis type: y, qside, qout, qlong, qinv
      }//analysis type loop
    }//track object valid
    else continue;
  }
  //Printf("Np: %lf - Nn: %lf",fNp,fNn);

  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 p1      = 0.;
  Double_t pX1     = 0., pY1     = 0., pZ1     = 0.;
  Double_t eta1 = 0., rap1 = 0.;
  Double_t energy1 = 0.;
  Double_t phi1    = 0.;

  Short_t charge2  = 0;
  Double_t p2      = 0.;
  Double_t pX2     = 0., pY2     = 0., pZ2     = 0.;
  Double_t eta2 = 0., rap2 = 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++) {
    if(fAnalysisLevel == "ESD")
      track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
    else if(fAnalysisLevel == "AOD")
      track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
    else if(fAnalysisLevel == "MC")
      track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));


    if(track1) {
      charge1 = track1->Charge();
      p1      = track1->P();
      pX1     = track1->Px();
      pY1     = track1->Py();
      pZ1     = track1->Pz();
      energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
			    TMath::Power(track1->M(),2));
      phi1    = track1->Phi();

    }
    else continue;
    
    for(j = 0; j < i; j++) {
      if(fAnalysisLevel == "ESD")
	track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
      else if(fAnalysisLevel == "AOD")
	track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
      else if(fAnalysisLevel == "MC")
	track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
      
      if(track2) {
	charge2 = track2->Charge();
	p2      = track2->P();
	pX2     = track2->Px();
	pY2     = track2->Py();
	pZ2     = track2->Pz();
	energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
			      TMath::Power(track2->M(),2));
	phi2    = track2->Phi();

	// filling the arrays

	// RAPIDITY 
	rap1 = 0.5*log((energy1 + pZ1)/(energy1 - pZ1)); 
	rap2 = 0.5*log((energy2 + pZ2)/(energy2 - pZ2)); 
	dy = TMath::Abs(rap1 - rap2);

	// Eta
	eta1 = track1->Eta();
	eta2 = track2->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 snn = TMath::Power(eTot,2) - TMath::Power(pxTot,2) - TMath::Power(pyTot,2) - TMath::Power(pzTot,2);
	Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
	
	qLong = TMath::Abs(eTot*qzTot - pzTot*q0Tot)/TMath::Sqrt(snn + TMath::Power(ptTot,2));
	
	//qout
	qOut = TMath::Sqrt(snn/(snn + TMath::Power(ptTot,2))) * TMath::Abs(pxTot*qxTot + pyTot*qyTot)/ptTot;
	
	//qside
	qSide = TMath::Abs(pxTot*qyTot - pyTot*qxTot)/ptTot;
	
	//qinv
	qInv = TMath::Sqrt(TMath::Abs(-TMath::Power(q0Tot,2) +TMath::Power(qxTot,2) +TMath::Power(qyTot,2) +TMath::Power(qzTot,2)));
	
	//phi
	dphi = phi1 - phi2;

	//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 == kRapidity) {
	    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(dy);
		}
		if((charge1 < 0)&&(charge2 < 0)) {
		  fNnn[kRapidity][iBin] += 1.;
		  fHistNN[kRapidity]->Fill(dy);
		}
		if((charge1 > 0)&&(charge2 < 0)) {
		  fNpn[kRapidity][iBin] += 1.;
		  fHistPN[kRapidity]->Fill(dy);
		}
		if((charge1 < 0)&&(charge2 > 0)) {
		  fNpn[kRapidity][iBin] += 1.;
		  fHistPN[kRapidity]->Fill(dy);
		}
	      }
	    }
	  }//rapidity
	
	  if(iAnalysisType == 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(deta);
		}
		if((charge1 < 0)&&(charge2 < 0)) {
		  fNnn[kEta][iBin] += 1.;
		  fHistNN[kEta]->Fill(deta);
		}
		if((charge1 > 0)&&(charge2 < 0)) {
		  fNpn[kEta][iBin] += 1.;
		  fHistPN[kEta]->Fill(deta);
		}
		if((charge1 < 0)&&(charge2 > 0)) {
		  fNpn[kEta][iBin] += 1.;
		  fHistPN[kEta]->Fill(deta);
		}
	      }
	    }
	  }//pseudorapidity
	  
	  // Qlong, out, side, inv
	  // thresholds missing!
	  if(iAnalysisType == kQlong) {
	    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(qLong);
		}
		if((charge1 < 0)&&(charge2 < 0)) {
		  fNnn[kQlong][iBin] += 1.;
		  fHistNN[kQlong]->Fill(qLong);
		}
		if((charge1 > 0)&&(charge2 < 0)) {
		  fNpn[kQlong][iBin] += 1.;
		  fHistPN[kQlong]->Fill(qLong);
		}
		if((charge1 < 0)&&(charge2 > 0)) {
		  fNpn[kQlong][iBin] += 1.;
		  fHistPN[kQlong]->Fill(qLong);
		}
	      }
	    }
	  }//qLong

	  if(iAnalysisType == kQout) {
	    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(qOut);
		}
		if((charge1 < 0)&&(charge2 < 0)) {
		  fNnn[kQout][iBin] += 1.;
		  fHistNN[kQout]->Fill(qOut);
		}
		if((charge1 > 0)&&(charge2 < 0)) {
		  fNpn[kQout][iBin] += 1.;
		  fHistPN[kQout]->Fill(qOut);
		}
		if((charge1 < 0)&&(charge2 > 0)) {
		  fNpn[kQout][iBin] += 1.;
		  fHistPN[kQout]->Fill(qOut);
		}
	      }
	    }
	  }//qOut

	  if(iAnalysisType == kQside) {
	    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(qSide);
		}
		if((charge1 < 0)&&(charge2 < 0)) {
		  fNnn[kQside][iBin] += 1.;
		  fHistNN[kQside]->Fill(qSide);
		}
		if((charge1 > 0)&&(charge2 < 0)) {
		  fNpn[kQside][iBin] += 1.;
		  fHistPN[kQside]->Fill(qSide);
		}
		if((charge1 < 0)&&(charge2 > 0)) {
		  fNpn[kQside][iBin] += 1.;
		  fHistPN[kQside]->Fill(qSide);
		}
	      }
	    }
	  }//qSide
	
	  if(iAnalysisType == kQinv) {
	    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(qInv);
		}
		if((charge1 < 0)&&(charge2 < 0)) {
		  fNnn[kQinv][iBin] += 1.;
		  fHistNN[kQinv]->Fill(qInv);
		}
		if((charge1 > 0)&&(charge2 < 0)) {
		  fNpn[kQinv][iBin] += 1.;
		  fHistPN[kQinv]->Fill(qInv);
		}
		if((charge1 < 0)&&(charge2 > 0)) {
		  fNpn[kQinv][iBin] += 1.;
		  fHistPN[kQinv]->Fill(qInv);
		}
	      }
	    }
	  }//qInv

	  // Phi
	  if(iAnalysisType == 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(dphi);
		}
		if((charge1 < 0)&&(charge2 < 0)) {
		  fNnn[kPhi][iBin] += 1.;
		  fHistNN[kPhi]->Fill(dphi);
		}
		if((charge1 > 0)&&(charge2 < 0)) {
		  fNpn[kPhi][iBin] += 1.;
		  fHistPN[kPhi]->Fill(dphi);
		}
		if((charge1 < 0)&&(charge2 > 0)) {
		  fNpn[kPhi][iBin] += 1.;
		  fHistPN[kPhi]->Fill(dphi);
		}
	      }
	    }
	  }//phi
	}//analysis type loop
      }//track2 valid
      else continue;
    }//end of 2nd particle loop
  }//end of 1st particle loop
  //Printf("Number of analyzed events: %i",fAnalyzedEvents);
}  


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

  return ferror[a][p2];
}
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>
d560b581	26	#include <TAxis.h>
b966652f	27	#include <TH1D.h>
6c178944	28	#include <TLorentzVector.h>
5c33329d	29	#include <TObjArray.h>
6c178944	30	#include <TGraphErrors.h>
9d1f0df5	31	#include <TString.h>
6c178944	32
cd54a838	33	#include "AliVParticle.h"
cd54a838	34	#include "AliMCParticle.h"
5c33329d	35	#include "AliESDtrack.h"
cd54a838	36	#include "AliAODTrack.h"
cd54a838	37
6c178944	38	#include "AliBalance.h"
	39
	40	ClassImp(AliBalance)
	41
7f0257ea	42	//____________________________________________________________________//
5d534bd4	43	AliBalance::AliBalance() :
5d534bd4	44	TObject(),
3feee083	45	fAnalysisLevel("ESD"),
b966652f	46	fAnalyzedEvents(0) {
6c178944	47	// Default constructor
3feee083	48
b966652f	49	for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
	50	if(i == 6) {
	51	fNumberOfBins[i] = 180;
	52	fP1Start[i] = 0.0;
	53	fP1Stop[i] = 360.0;
	54	fP2Start[i] = -360.0;
	55	fP2Stop[i] = 360.0;
	56	fP2Step[i] = 0.1;
	57	}
	58	else {
	59	fNumberOfBins[i] = 20;
	60	fP1Start[i] = -1.0;
	61	fP1Stop[i] = 1.0;
	62	fP2Start[i] = 0.0;
	63	fP2Stop[i] = 2.0;
	64	}
	65	fP2Step[i] = TMath::Abs(fP2Start - fP2Stop) / (Double_t)fNumberOfBins[i];
	66	fNn[i] = 0.0;
	67	fNp[i] = 0.0;
3feee083	68
	69	for(Int_t j = 0; j < MAXIMUM_NUMBER_OF_STEPS; j++) {
	70	fNpp[i][j] = .0;
	71	fNnn[i][j] = .0;
	72	fNpn[i][j] = .0;
	73	fB[i][j] = 0.0;
	74	ferror[i][j] = 0.0;
	75	}
09bb7bf4	76	}
b966652f	77
	78	TString gAnalysisType[ANALYSIS_TYPES] = {"y","eta","qlong","qout","qside","qinv","phi"};
	79	TString histName;
	80	for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
	81	histName = "fHistP"; histName += gAnalysisType[iAnalysisType];
	82	fHistP[iAnalysisType] = new TH1D(histName.Data(),"",100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
	83	histName = "fHistN"; histName += gAnalysisType[iAnalysisType];
	84	fHistN[iAnalysisType] = new TH1D(histName.Data(),"",100,fP1Start[iAnalysisType],fP1Stop[iAnalysisType]);
	85
	86	histName = "fHistPN"; histName += gAnalysisType[iAnalysisType];
	87	fHistPN[iAnalysisType] = new TH1D(histName.Data(),"",fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
	88	histName = "fHistNP"; histName += gAnalysisType[iAnalysisType];
	89	fHistNP[iAnalysisType] = new TH1D(histName.Data(),"",fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
	90	histName = "fHistPP"; histName += gAnalysisType[iAnalysisType];
	91	fHistPP[iAnalysisType] = new TH1D(histName.Data(),"",fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
	92	histName = "fHistNN"; histName += gAnalysisType[iAnalysisType];
	93	fHistNN[iAnalysisType] = new TH1D(histName.Data(),"",fNumberOfBins[iAnalysisType],fP2Start[iAnalysisType],fP2Stop[iAnalysisType]);
	94
	95	}
6c178944	96	}
6c178944	97
6c178944	98
7f0257ea	99	//____________________________________________________________________//
7f0257ea	100	AliBalance::AliBalance(const AliBalance& balance):
9d1f0df5	101	TObject(balance), fAnalysisLevel(balance.fAnalysisLevel),
b966652f	102	fAnalyzedEvents(balance.fAnalyzedEvents) {
7f0257ea	103	//copy constructor
b966652f	104	for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
	105	fNn[i] = balance.fNn[i];
	106	fNp[i] = balance.fNp[i];
3feee083	107
b966652f	108	fP1Start[i] = balance.fP1Start[i];
b966652f	109	fP1Stop[i] = balance.fP1Stop[i];
3feee083	110	fNumberOfBins[i] = balance.fNumberOfBins[i];
	111	fP2Start[i] = balance.fP2Start[i];
	112	fP2Stop[i] = balance.fP2Stop[i];
	113	fP2Step[i] = balance.fP2Step[i];
	114
b966652f	115	fHistP[i] = balance.fHistP[i];
	116	fHistN[i] = balance.fHistN[i];
	117	fHistPN[i] = balance.fHistPN[i];
	118	fHistNP[i] = balance.fHistNP[i];
	119	fHistPP[i] = balance.fHistPP[i];
	120	fHistNN[i] = balance.fHistNN[i];
	121
3feee083	122	for(Int_t j = 0; j < MAXIMUM_NUMBER_OF_STEPS; j++) {
	123	fNpp[i][j] = .0;
	124	fNnn[i][j] = .0;
	125	fNpn[i][j] = .0;
	126	fB[i][j] = 0.0;
	127	ferror[i][j] = 0.0;
	128	}
	129	}
	130	}
	131
7f0257ea	132
	133	//____________________________________________________________________//
	134	AliBalance::~AliBalance() {
6c178944	135	// Destructor
6c178944	136
7f0257ea	137	}
7f0257ea	138
7f0257ea	139
7f0257ea	140	//____________________________________________________________________//
3feee083	141	void AliBalance::SetNumberOfBins(Int_t ibin, Int_t ibins) {
6c178944	142	// Sets the number of bins for the analyzed interval
3feee083	143	// Set the same Information for all analyses
3feee083	144	if(ibin == -1){
b966652f	145	for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
3feee083	146	fNumberOfBins[i] = ibins;
	147	}
	148	}
	149	// Set the Information for one analysis
b966652f	150	else if(ibin > -1 && ibin < ANALYSIS_TYPES){
3feee083	151	fNumberOfBins[ibin] = ibins;
	152	}
	153	else{
	154	AliError("Wrong ANALYSIS number!");
	155	}
6c178944	156	}
6c178944	157
7f0257ea	158	//____________________________________________________________________//
b966652f	159	void AliBalance::SetInterval(Double_t p1Start, Double_t p1Stop,
b966652f	160	Int_t ibin, Double_t p2Start, Double_t p2Stop) {
6c178944	161	// Sets the analyzed interval.
3feee083	162	// Set the same Information for all analyses
3feee083	163	if(ibin == -1){
b966652f	164	for(Int_t i = 0; i < ANALYSIS_TYPES; i++){
	165	fP1Start[i] = p1Start;
	166	fP1Stop[i] = p1Stop;
3feee083	167	fP2Start[i] = p2Start;
	168	fP2Stop[i] = p2Stop;
	169	fP2Step[i] = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins[i];
	170	}
7f0257ea	171	}
3feee083	172	// Set the Information for one analysis
b966652f	173	else if(ibin > -1 && ibin < ANALYSIS_TYPES){
	174	fP1Start[ibin] = p1Start;
	175	fP1Stop[ibin] = p1Stop;
3feee083	176	fP2Start[ibin] = p2Start;
	177	fP2Stop[ibin] = p2Stop;
	178	fP2Step[ibin] = TMath::Abs(p2Start - p2Stop) / (Double_t)fNumberOfBins[ibin];
7f0257ea	179	}
3feee083	180	else{
3feee083	181	AliError("Wrong ANALYSIS number!");
7f0257ea	182	}
6c178944	183	}
6c178944	184
3feee083	185
898d5b28	186	//____________________________________________________________________//
09bb7bf4	187	void AliBalance::PrintAnalysisSettings() {
09bb7bf4	188
	189	Printf("======================================");
	190	Printf("Analysis level: %s",fAnalysisLevel.Data());
3feee083	191	Printf("======================================");
b966652f	192	for(Int_t ibin = 0; ibin < ANALYSIS_TYPES; ibin++){
3feee083	193	Printf("Interval info for variable %d",ibin);
	194	Printf("Analyzed interval (min.): %lf",fP2Start[ibin]);
	195	Printf("Analyzed interval (max.): %lf",fP2Stop[ibin]);
	196	Printf("Number of bins: %d",fNumberOfBins[ibin]);
	197	Printf("Step: %lf",fP2Step[ibin]);
	198	Printf(" ");
	199	}
09bb7bf4	200	Printf("======================================");
898d5b28	201	}
898d5b28	202
7f0257ea	203	//____________________________________________________________________//
5c33329d	204	void AliBalance::CalculateBalance(TObjArray *gTrackArray) {
6c178944	205	// Calculates the balance function
	206	fAnalyzedEvents++;
	207	Int_t i = 0 , j = 0;
3feee083	208	Int_t iBin = 0;
5c33329d	209
cd54a838	210	AliVParticle* track = 0;
	211	AliVParticle* track1 = 0;
	212	AliVParticle* track2 = 0;
9d1f0df5	213
2cd42194	214	//Printf("(AliBalance) Number of tracks: %d",gTrackArray->GetEntries());
5c33329d	215	Int_t gNtrack = gTrackArray->GetEntries();
5c33329d	216	for(i = 0; i < gNtrack; i++) {
9d1f0df5	217	if(fAnalysisLevel == "ESD")
9d1f0df5	218	track = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
cd54a838	219	else if(fAnalysisLevel == "AOD")
	220	track = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
	221	else if(fAnalysisLevel == "MC")
	222	track = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
c511520f	223
c511520f	224	if(track) {
b966652f	225	Short_t charge = track->Charge();
	226	Double_t pseudorapidity = track->Eta();
	227	Double_t rapidity = track->Y();
	228	Double_t phi = TMath::ATan(track->Py()/track->Px())*180.0/TMath::Pi();
	229
	230	//0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
	231	for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
	232	if(iAnalysisType == 1) {
	233	if((pseudorapidity >= fP1Start[iAnalysisType]) && (pseudorapidity <= fP1Stop[iAnalysisType])) {
	234	if(charge > 0) {
	235	fNp[iAnalysisType] += 1.;
	236	fHistP[iAnalysisType]->Fill(pseudorapidity);
	237	}//charge > 0
	238	if(charge < 0) {
	239	fNn[iAnalysisType] += 1.;
	240	fHistN[iAnalysisType]->Fill(pseudorapidity);
	241	}//charge < 0
	242	}//p1 interval check
	243	}//analysis type: eta
	244	if(iAnalysisType == 6) {
	245	if((phi >= fP1Start[iAnalysisType]) && (phi <= fP1Stop[iAnalysisType])) {
	246	if(charge > 0) {
	247	fNp[iAnalysisType] += 1.;
	248	fHistP[iAnalysisType]->Fill(phi);
	249	}//charge > 0
	250	if(charge < 0) {
	251	fNn[iAnalysisType] += 1.;
	252	fHistN[iAnalysisType]->Fill(phi);
	253	}//charge < 0
	254	}//p1 interval check
	255	}//analysis type: phi
	256	else {
	257	if((rapidity >= fP1Start[iAnalysisType]) && (rapidity <= fP1Stop[iAnalysisType])) {
	258	if(charge > 0) {
	259	fNp[iAnalysisType] += 1.;
	260	fHistP[iAnalysisType]->Fill(rapidity);
	261	}//charge > 0
	262	if(charge < 0) {
	263	fNn[iAnalysisType] += 1.;
	264	fHistN[iAnalysisType]->Fill(rapidity);
	265	}//charge < 0
	266	}//p1 interval check
	267	}//analysis type: y, qside, qout, qlong, qinv
	268	}//analysis type loop
	269	}//track object valid
c511520f	270	else continue;
7f0257ea	271	}
5c33329d	272	//Printf("Np: %lf - Nn: %lf",fNp,fNn);
6c178944	273
b966652f	274	Double_t dy = 0., deta = 0.;
	275	Double_t qLong = 0., qOut = 0., qSide = 0., qInv = 0.;
	276	Double_t dphi = 0.;
	277
	278	Short_t charge1 = 0;
	279	Double_t p1 = 0.;
	280	Double_t pX1 = 0., pY1 = 0., pZ1 = 0.;
	281	Double_t eta1 = 0., rap1 = 0.;
	282	Double_t energy1 = 0.;
	283	Double_t phi1 = 0.;
	284
	285	Short_t charge2 = 0;
	286	Double_t p2 = 0.;
	287	Double_t pX2 = 0., pY2 = 0., pZ2 = 0.;
	288	Double_t eta2 = 0., rap2 = 0.;
	289	Double_t energy2 = 0.;
	290	Double_t phi2 = 0.;
6c178944	291	//0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
3feee083	292	for(i = 1; i < gNtrack; i++) {
	293	if(fAnalysisLevel == "ESD")
	294	track1 = dynamic_cast<AliESDtrack *>(gTrackArray->At(i));
	295	else if(fAnalysisLevel == "AOD")
	296	track1 = dynamic_cast<AliAODTrack *>(gTrackArray->At(i));
	297	else if(fAnalysisLevel == "MC")
	298	track1 = dynamic_cast<AliMCParticle *>(gTrackArray->At(i));
	299
3feee083	300
	301	if(track1) {
	302	charge1 = track1->Charge();
	303	p1 = track1->P();
	304	pX1 = track1->Px();
	305	pY1 = track1->Py();
	306	pZ1 = track1->Pz();
	307	energy1 = TMath::Sqrt(TMath::Power(track1->P(),2) +
	308	TMath::Power(track1->M(),2));
	309	phi1 = track1->Phi();
c511520f	310
7f0257ea	311	}
3feee083	312	else continue;
	313
	314	for(j = 0; j < i; j++) {
9d1f0df5	315	if(fAnalysisLevel == "ESD")
3feee083	316	track2 = dynamic_cast<AliESDtrack *>(gTrackArray->At(j));
cd54a838	317	else if(fAnalysisLevel == "AOD")
3feee083	318	track2 = dynamic_cast<AliAODTrack *>(gTrackArray->At(j));
cd54a838	319	else if(fAnalysisLevel == "MC")
3feee083	320	track2 = dynamic_cast<AliMCParticle *>(gTrackArray->At(j));
	321
	322	if(track2) {
	323	charge2 = track2->Charge();
	324	p2 = track2->P();
	325	pX2 = track2->Px();
	326	pY2 = track2->Py();
	327	pZ2 = track2->Pz();
	328	energy2 = TMath::Sqrt(TMath::Power(track2->P(),2) +
	329	TMath::Power(track2->M(),2));
	330	phi2 = track2->Phi();
	331
	332	// filling the arrays
	333
	334	// RAPIDITY
b966652f	335	rap1 = 0.5*log((energy1 + pZ1)/(energy1 - pZ1));
	336	rap2 = 0.5*log((energy2 + pZ2)/(energy2 - pZ2));
	337	dy = TMath::Abs(rap1 - rap2);
	338
3feee083	339	// Eta
b966652f	340	eta1 = track1->Eta();
	341	eta2 = track2->Eta();
	342	deta = TMath::Abs(eta1 - eta2);
3feee083	343
b966652f	344	//qlong
3feee083	345	Double_t eTot = energy1 + energy2;
	346	Double_t pxTot = pX1 + pX2;
	347	Double_t pyTot = pY1 + pY2;
	348	Double_t pzTot = pZ1 + pZ2;
	349	Double_t q0Tot = energy1 - energy2;
	350	Double_t qxTot = pX1 - pX2;
	351	Double_t qyTot = pY1 - pY2;
	352	Double_t qzTot = pZ1 - pZ2;
b966652f	353	Double_t snn = TMath::Power(eTot,2) - TMath::Power(pxTot,2) - TMath::Power(pyTot,2) - TMath::Power(pzTot,2);
3feee083	354	Double_t ptTot = TMath::Sqrt( TMath::Power(pxTot,2) + TMath::Power(pyTot,2));
3feee083	355
b966652f	356	qLong = TMath::Abs(eTotqzTot - pzTotq0Tot)/TMath::Sqrt(snn + TMath::Power(ptTot,2));
	357
	358	//qout
	359	qOut = TMath::Sqrt(snn/(snn + TMath::Power(ptTot,2))) * TMath::Abs(pxTotqxTot + pyTotqyTot)/ptTot;
	360
	361	//qside
	362	qSide = TMath::Abs(pxTotqyTot - pyTotqxTot)/ptTot;
3feee083	363
b966652f	364	//qinv
	365	qInv = TMath::Sqrt(TMath::Abs(-TMath::Power(q0Tot,2) +TMath::Power(qxTot,2) +TMath::Power(qyTot,2) +TMath::Power(qzTot,2)));
	366
	367	//phi
	368	dphi = phi1 - phi2;
	369
	370	//0:y - 1:eta - 2:Qlong - 3:Qout - 4:Qside - 5:Qinv - 6:phi
	371	for(Int_t iAnalysisType = 0; iAnalysisType < ANALYSIS_TYPES; iAnalysisType++) {
	372	if(iAnalysisType == kRapidity) {
	373	if( dy > fP2Start[kRapidity] && dy < fP2Stop[kRapidity]){
	374	iBin = Int_t((dy-fP2Start[kRapidity])/fP2Step[kRapidity]);
	375	if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
	376	if((charge1 > 0)&&(charge2 > 0)) {
	377	fNpp[kRapidity][iBin] += 1.;
	378	fHistPP[kRapidity]->Fill(dy);
	379	}
	380	if((charge1 < 0)&&(charge2 < 0)) {
	381	fNnn[kRapidity][iBin] += 1.;
	382	fHistNN[kRapidity]->Fill(dy);
	383	}
	384	if((charge1 > 0)&&(charge2 < 0)) {
	385	fNpn[kRapidity][iBin] += 1.;
	386	fHistPN[kRapidity]->Fill(dy);
	387	}
	388	if((charge1 < 0)&&(charge2 > 0)) {
	389	fNpn[kRapidity][iBin] += 1.;
	390	fHistPN[kRapidity]->Fill(dy);
	391	}
	392	}
	393	}
	394	}//rapidity
	395
	396	if(iAnalysisType == kEta) {
	397	if( deta > fP2Start[kEta] && deta < fP2Stop[kEta]){
	398	iBin = Int_t((deta-fP2Start[kEta])/fP2Step[kEta]);
	399	if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
	400	if((charge1 > 0)&&(charge2 > 0)) {
	401	fNpp[kEta][iBin] += 1.;
	402	fHistPP[kEta]->Fill(deta);
	403	}
	404	if((charge1 < 0)&&(charge2 < 0)) {
	405	fNnn[kEta][iBin] += 1.;
	406	fHistNN[kEta]->Fill(deta);
	407	}
	408	if((charge1 > 0)&&(charge2 < 0)) {
	409	fNpn[kEta][iBin] += 1.;
	410	fHistPN[kEta]->Fill(deta);
	411	}
	412	if((charge1 < 0)&&(charge2 > 0)) {
	413	fNpn[kEta][iBin] += 1.;
	414	fHistPN[kEta]->Fill(deta);
	415	}
	416	}
	417	}
	418	}//pseudorapidity
	419
	420	// Qlong, out, side, inv
	421	// thresholds missing!
	422	if(iAnalysisType == kQlong) {
	423	if( qLong > fP2Start[kQlong] && qLong < fP2Stop[kQlong]){
	424	iBin = Int_t((qLong-fP2Start[kQlong])/fP2Step[kQlong]);
	425	if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
	426	if((charge1 > 0)&&(charge2 > 0)) {
	427	fNpp[kQlong][iBin] += 1.;
428	fHistPP[kQlong]->Fill(qLong);
429	}
430	if((charge1 < 0)&&(charge2 < 0)) {
431	fNnn[kQlong][iBin] += 1.;
432	fHistNN[kQlong]->Fill(qLong);
433	}
434	if((charge1 > 0)&&(charge2 < 0)) {
435	fNpn[kQlong][iBin] += 1.;
436	fHistPN[kQlong]->Fill(qLong);
437	}
438	if((charge1 < 0)&&(charge2 > 0)) {
439	fNpn[kQlong][iBin] += 1.;
440	fHistPN[kQlong]->Fill(qLong);
441	}
442	}
443	}
444	}//qLong
445
446	if(iAnalysisType == kQout) {
447	if( qOut > fP2Start[kQout] && qOut < fP2Stop[kQout]){
448	iBin = Int_t((qOut-fP2Start[kQout])/fP2Step[kQout]);
449	if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
450	if((charge1 > 0)&&(charge2 > 0)) {
451	fNpp[kQout][iBin] += 1.;
452	fHistPP[kQout]->Fill(qOut);
453	}
454	if((charge1 < 0)&&(charge2 < 0)) {
455	fNnn[kQout][iBin] += 1.;
456	fHistNN[kQout]->Fill(qOut);
457	}
458	if((charge1 > 0)&&(charge2 < 0)) {
459	fNpn[kQout][iBin] += 1.;
460	fHistPN[kQout]->Fill(qOut);
461	}
462	if((charge1 < 0)&&(charge2 > 0)) {
463	fNpn[kQout][iBin] += 1.;
464	fHistPN[kQout]->Fill(qOut);
465	}
466	}
467	}
468	}//qOut
469
470	if(iAnalysisType == kQside) {
471	if( qSide > fP2Start[kQside] && qSide < fP2Stop[kQside]){
472	iBin = Int_t((qSide-fP2Start[kQside])/fP2Step[kQside]);
473	if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
474	if((charge1 > 0)&&(charge2 > 0)) {
475	fNpp[kQside][iBin] += 1.;
476	fHistPP[kQside]->Fill(qSide);
477	}
478	if((charge1 < 0)&&(charge2 < 0)) {
479	fNnn[kQside][iBin] += 1.;
480	fHistNN[kQside]->Fill(qSide);
481	}
482	if((charge1 > 0)&&(charge2 < 0)) {
483	fNpn[kQside][iBin] += 1.;
484	fHistPN[kQside]->Fill(qSide);
485	}
486	if((charge1 < 0)&&(charge2 > 0)) {
487	fNpn[kQside][iBin] += 1.;
488	fHistPN[kQside]->Fill(qSide);
489	}
490	}
491	}
492	}//qSide
493
494	if(iAnalysisType == kQinv) {
495	if( qInv > fP2Start[kQinv] && qInv < fP2Stop[kQinv]){
496	iBin = Int_t((qInv-fP2Start[kQinv])/fP2Step[kQinv]);
497	if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
498	if((charge1 > 0)&&(charge2 > 0)) {
499	fNpp[kQinv][iBin] += 1.;
500	fHistPP[kQinv]->Fill(qInv);
501	}
502	if((charge1 < 0)&&(charge2 < 0)) {
503	fNnn[kQinv][iBin] += 1.;
504	fHistNN[kQinv]->Fill(qInv);
505	}
506	if((charge1 > 0)&&(charge2 < 0)) {
507	fNpn[kQinv][iBin] += 1.;
508	fHistPN[kQinv]->Fill(qInv);
509	}
510	if((charge1 < 0)&&(charge2 > 0)) {
511	fNpn[kQinv][iBin] += 1.;
512	fHistPN[kQinv]->Fill(qInv);
513	}
514	}
515	}
516	}//qInv
517
518	// Phi
519	if(iAnalysisType == kPhi) {
520	if( dphi > fP2Start[kPhi] && dphi < fP2Stop[kPhi]){
521	iBin = Int_t((dphi-fP2Start[kPhi])/fP2Step[kPhi]);
522	if(iBin >=0 && iBin < MAXIMUM_NUMBER_OF_STEPS){
523	if((charge1 > 0)&&(charge2 > 0)) {
524	fNpp[kPhi][iBin] += 1.;
525	fHistPP[kPhi]->Fill(dphi);
526	}
527	if((charge1 < 0)&&(charge2 < 0)) {
528	fNnn[kPhi][iBin] += 1.;
529	fHistNN[kPhi]->Fill(dphi);
530	}
531	if((charge1 > 0)&&(charge2 < 0)) {
532	fNpn[kPhi][iBin] += 1.;
533	fHistPN[kPhi]->Fill(dphi);
534	}
535	if((charge1 < 0)&&(charge2 > 0)) {
536	fNpn[kPhi][iBin] += 1.;
537	fHistPN[kPhi]->Fill(dphi);
538	}
539	}
540	}
541	}//phi
542	}//analysis type loop
543	}//track2 valid
3feee083	544	else continue;
b966652f	545	}//end of 2nd particle loop
b966652f	546	}//end of 1st particle loop
3feee083	547	//Printf("Number of analyzed events: %i",fAnalyzedEvents);
3feee083	548	}
1d8792e7	549
1d8792e7	550
1d8792e7	551
7f0257ea	552	//____________________________________________________________________//
3feee083	553	Double_t AliBalance::GetBalance(Int_t a, Int_t p2) {
6c178944	554	// Returns the value of the balance function in bin p2
b966652f	555	fB[a][p2] = 0.5(((fNpn[a][p2] - 2.0fNnn[a][p2])/fNn[a]) + ((fNpn[a][p2] - 2.0*fNpp[a][p2])/fNp[a]))/fP2Step[a];
7f0257ea	556
3feee083	557	return fB[a][p2];
6c178944	558	}
6c178944	559
7f0257ea	560	//____________________________________________________________________//
3feee083	561	Double_t AliBalance::GetError(Int_t a, Int_t p2) {
6c178944	562	// Returns the error on the BF value for bin p2
b966652f	563	ferror[a][p2] = TMath::Sqrt( Double_t(fNpp[a][p2])/(Double_t(fNp[a])Double_t(fNp[a])) + Double_t(fNnn[a][p2])/(Double_t(fNn[a])Double_t(fNn[a])) + Double_t(fNpn[a][p2])*TMath::Power((0.5/Double_t(fNp[a]) + 0.5/Double_t(fNn[a])),2))/fP2Step[a];
6c178944	564
3feee083	565	return ferror[a][p2];
6c178944	566	}
6c178944	567
898d5b28	568