[u/mrichter/AliRoot.git] / PWGCF / EBYE / LRC / AliLRCAnalysis.cxx

/**************************************************************************
 * Author: Andrey Ivanov.                                           *
 * 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.                  *
 **************************************************************************/
//-------------------------------------------------------------------------
//    Description: 
//    This class include into LRC library for Long-Range Correlation analysis
//    it is base class for NN, PtN, PtPt
//    implements base methods for thees classes
//    Origin: Petr Naumenko, SPbSU-CERN, Petr.Naoumenko@cern.ch ,
//    Andrey Ivanov (SPbSU-CERN), Igor Altsebeev (SPbSU-CERN) 
//-------------------------------------------------------------------------

#include "AliLRCAnalysis.h"
#include "TFile.h"
#include "AliLRCFit.h"
#include "TProfile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TPaveText.h"
#include "TF1.h"
#include "math.h"
#include "TStyle.h"
class gStyle;


ClassImp(AliLRCAnalysis) 

/******************************************************
 * AliLRCAnalysis class
 ******************************************************/

AliLRCAnalysis::AliLRCAnalysis():TObject(), fPrAbs(new TH1D()), fPrRel(new TH1D()), fPrf(new TH1D()), fPrb(new TH1D()), fileHist(new TFile()), fdptb(.0), fEntries(0), fSx((char*)" "), fSy((char*)" "), fxFitMin(.0), fxFitMax(.0), fNsigma(2.), farel(.0), fbrel(.0), farelError(.0), fbrelError(.0), fXi2rel(.0), faabs(.0), fbabs(.0), faabsError(.0), fbabsError(.0), fXi2abs(.0),fFitMethod(0){
//Empty constructor
}

AliLRCAnalysis::AliLRCAnalysis(const AliLRCAnalysis& a):TObject(a),fPrAbs(a.fPrAbs), fPrRel(a.fPrRel), fPrf(a.fPrf), fPrb(a.fPrb), fileHist(a.fileHist), fdptb(a.fdptb), fEntries(a.fEntries), fSx(a.fSx), fSy(a.fSy), fxFitMin(a.fxFitMin), fxFitMax(a.fxFitMax), fNsigma(a.fNsigma), farel(a.farel), fbrel(a.fbrel), farelError(a.farelError), fbrelError(a.fbrelError), fXi2rel(a.fXi2rel), faabs(a.faabs), fbabs(a.fbabs),faabsError(a.faabsError), fbabsError(a.fbabsError), fXi2abs(a.fXi2abs),fFitMethod(a.fFitMethod){
//Constructor
}

AliLRCAnalysis& AliLRCAnalysis::operator= (const AliLRCAnalysis& a){
//Operator =
	if(this!=&a){
		TObject::operator= (a);
		fPrAbs = a.fPrAbs;
		fPrRel = a.fPrRel;
		fPrf = a.fPrf;
		fPrb = a.fPrb;
		fileHist = a.fileHist;
		fSx = a.fSx;
		fSy = a.fSy;
		fdptb = a.fdptb;
		fEntries = a.fEntries;
		fxFitMin = a.fxFitMin;
		fxFitMax = a.fxFitMax;
		fNsigma = a.fNsigma;
 		farel = a.farel;
		fbrel = a.fbrel;
 		farelError = a.farelError;
		fbrelError = a.fbrelError;
		fXi2rel = a.fXi2rel;
		faabs = a.faabs;
		fbabs = a.fbabs; 
		faabsError = a.faabsError;
		fbabsError = a.fbabsError; 
		fXi2abs = a.fXi2abs;
		fFitMethod = a.fFitMethod;
	}
	return *this;
}

AliLRCAnalysis::~AliLRCAnalysis() {
//Destructor
    delete fPrAbs;
    delete fPrRel;
    delete fileHist;
	delete fPrf;
	delete fPrb;
}


void AliLRCAnalysis::CreateHist(char *name, char *nameAbs, char *nameRel, char *atitleF, char *atitleB, char *rtitleF, char *rtitleB, TH2D* sourceHist) {
//Create absolute and relation var histogramm
    TProfile* profX = (TProfile*) sourceHist->ProfileX(name, 1, sourceHist->GetNbinsY());
    fEntries = (int) sourceHist->GetEntries();
    fPrAbs = new TH1D(nameAbs, profX->GetTitle(), profX->GetXaxis()->GetNbins(), profX->GetXaxis()->GetXmin(), profX->GetXaxis()->GetXmax());
    fPrAbs->SetOption("E");
    fPrRel = new TH1D(nameRel, profX->GetTitle(), profX->GetXaxis()->GetNbins(), profX->GetXaxis()->GetXmin()/sourceHist->ProjectionX()->GetMean(), profX->GetXaxis()->GetXmax()/sourceHist->ProjectionX()->GetMean());
    fPrRel->SetOption("E");
    fPrAbs->GetXaxis()->SetTitle(atitleF);
    fPrAbs->GetYaxis()->SetTitle(atitleB);
    fPrRel->GetXaxis()->SetTitle(rtitleF);
    fPrRel->GetYaxis()->SetTitle(rtitleB);
    fPrf = (TH1D*) sourceHist->ProjectionX();
    fPrb = (TH1D*) sourceHist->ProjectionY();
    fPrf->GetXaxis()->SetTitle(atitleF);
    fPrf->GetYaxis()->SetTitle("Tracks");
    fPrb->GetXaxis()->SetTitle(atitleB);
    fPrb->GetYaxis()->SetTitle("Tracks");
    fSx = atitleF;
    fSy = atitleB;
    double mnf = fPrf->GetMean();
    double rms = fPrf->GetRMS();
    fxFitMin = mnf-fNsigma*rms;//sqrt(rms);
    fxFitMax = mnf+fNsigma*rms;//sqrt(rms);//mnf+2*sqrt(rms);
    //delete profX;

}

void AliLRCAnalysis::SetBinsRange(int binMin, int binMax){
//Set the bin range
	TH1D* h=fPrf;
    int n=h->GetNbinsX();
    fxFitMin = h->GetXaxis()->GetXmin();
    fxFitMax = h->GetXaxis()->GetXmax();
    double df = (fxFitMax-fxFitMin)/n;
    for(int i=1; i<=n; i++)
	if(h->GetBinContent(i) != 0){
		fxFitMin +=  (i + binMin) * df;
		break;
	}
    for(int i=1; i<=n; i++)
	if(h->GetBinContent(n-i) != 0){
		fxFitMax -= (i + binMax) * df;
		break;
	}
}

bool AliLRCAnalysis::SetFitRange(double xMin, double xMax){
//Set the fit range
	if(xMax < xMin){
		return false;
	}
	this->fxFitMin = xMin;
	this->fxFitMax = xMax;
	return true;
}
bool AliLRCAnalysis::SetFitRangeMin(double xMin){
//Set the fit range min
	fxFitMin = xMin;
	return true;
}

bool AliLRCAnalysis::SetFitRangeMax(double xMax){
//Set the fit range min
	//this->fxFitMax = xMax;

	TH1D* h=fPrf;
	double maxBorder = 0.;

	for ( int binI = h->GetNbinsX(); binI > 0; binI-- )
	{
		if ( h->GetBinContent(binI) != 0 )
		{
			maxBorder = h->GetBinLowEdge(binI+1) ;
			break;
		}	
	}
	fxFitMax = TMath::Max( fxFitMin, TMath::Min( xMax, maxBorder ) );

	return true;
}

void AliLRCAnalysis::SetFullFitRange(){
//Set fitting on full range
	TH1D* h=fPrf;
    fxFitMin = h->GetXaxis()->GetXmin();
	
	for ( int binI = h->GetNbinsX(); binI > 0; binI-- )
	{
		if ( h->GetBinContent(binI) != 0 )
		{
			fxFitMax = h->GetBinLowEdge(binI+1) ;
			break;
		}	
	}
}


void AliLRCAnalysis::SetXmin(double xMin){
	fxFitMin = xMin;
}

void AliLRCAnalysis::SetXmax(double xMax){
	fxFitMax = xMax;
}
void AliLRCAnalysis::SetNsigma(double nSigma){
//Sets fiting range in forward value RMSs
	fNsigma = nSigma;
	double mnf = fPrf->GetMean();
	double rms = fPrf->GetRMS();
	fxFitMin = mnf-fNsigma*rms;
	fxFitMax = mnf+fNsigma*rms;

//	cout << "### Fit params: N sigma is " << fNsigma 
//		<< ", mean = " << mnf << ", rms = " << rms << ", minFitX = " << fxFitMin << ", maxFitX = " << fxFitMax << endl;
}

void AliLRCAnalysis::SetFitMethod(int id){
//Choose fit method
	fFitMethod = id;
}
double AliLRCAnalysis::GetFitXmin() const
{
	return fxFitMin;
}

double AliLRCAnalysis::GetFitXmax() const 
{
	return fxFitMax;
}
double AliLRCAnalysis::GetArel() const {
	return farel;
}

double AliLRCAnalysis::GetBrel() const {
	return fbrel;
}

double AliLRCAnalysis::GetArelError() const {
	return farelError;
}

double AliLRCAnalysis::GetBrelError() const {
	return fbrelError;
}

double AliLRCAnalysis::GetXi2rel() const {
	return fXi2rel;
}

double AliLRCAnalysis::GetAabs() const {
	return faabs;
}

double AliLRCAnalysis::GetBabs() const {
	return fbabs;
}

double AliLRCAnalysis::GetAabsError() const {
	return faabsError;
}

double AliLRCAnalysis::GetBabsError() const {
	return fbabsError;
}
double AliLRCAnalysis::GetXi2abs() const {
	return fXi2abs;
}

void AliLRCAnalysis::Calculate(){
//Calculate all
	double mnf = fPrf->GetMean(); 

	// new fit for abs
	AliLRCFit *fit1 = new AliLRCFit(fPrRel, fxFitMin/mnf, fxFitMax/mnf, 1.);
	AliLRCFit *fit2 = new AliLRCFit(fPrAbs, fxFitMin, fxFitMax, 0.);
	farel = fit1->Geta();
	fbrel = fit1->Getb();

	faabs = fit2->Geta();
	fbabs = fit2->Getb();

	if ( fFitMethod==0 )
	{
		faabsError = fit2->Getda();
		fbabsError = fit2->Getdb();

		farelError = fit1->Getda();
		fbrelError = fit1->Getdb();
	}
	if ( fFitMethod==1 )
	{
		faabsError = fit2->Getda1();
		fbabsError = fit2->Getdb1();

		farelError = fit1->Getda1();
		fbrelError = fit1->Getdb1();
	}

	fXi2rel = fit1->Gethi2();
	fXi2abs = fit2->Gethi2();
}

void AliLRCAnalysis::DrawAbs() {
//Draw abs var hist with ALL info
	int * mas = new int [fgkPlotFlags];
	for ( int i = 0; i < fgkPlotFlags; i++ )
		mas[i] = 1;
	DrawAbsPure( mas, 1 );
	delete []mas;
}

void AliLRCAnalysis::DrawAbs( const int * const mas ) {
//Draw abs var hist with REQUESTED BY ARRAY info
	DrawAbsPure( mas, 1 );
}

void AliLRCAnalysis::DrawAbsPure( const int * const mDrawArray, bool drawPaveLabel ) {
	Calculate();
// Draw abs var histrogram
	DrawHist( mDrawArray, drawPaveLabel, faabs, fbabs, faabsError, fbabsError, fPrAbs, 0 );
}

void AliLRCAnalysis::DrawRel() {
//Draw rel var hist with ALL info
	int * mas = new int [fgkPlotFlags];
	for ( int i = 0; i < fgkPlotFlags; i++ )
		mas[i] = 1;
	DrawRelPure( mas, 1 );
	delete []mas;
}

void AliLRCAnalysis::DrawRel( const int * const mas ) {
//Draw rel var hist with REQUESTED BY ARRAY info
	DrawRelPure( mas, 1 );
}

void AliLRCAnalysis::DrawRelPure( const int * const mDrawArray, bool drawPaveLabel ) {
	Calculate();
// Draw rel var histogram
	DrawHist( mDrawArray, drawPaveLabel, farel, fbrel, farelError, fbrelError, fPrRel, 1 );
}
void AliLRCAnalysis::DrawHist( const int * const mDrawArray, bool drawPaveLabel, double aCoef, double bCoef,
		double aCoefError, double bCoefError, TH1D* profToDraw, int histType ) {
// Method called by DrawRelPure or DrawAbsPure to draw corresponding LRC  histo 
	double mnf = fPrf->GetMean();
	double y1, y2, x1, x2;
	Int_t i, n;
	char str[50];
    //mnf = fPrf->GetMean();
	profToDraw->SetStats(0);
    //profToDraw->Fit("pol1", "", "", fxFitMin, fxFitMax);
    //profToDraw->Fit("pol1");
    //TF1 *fit1 = profToDraw->GetFunction("pol1");

	//draw fit line
	TF1 *f1 = 0;
	if ( histType == 0 )	{
		f1 = new TF1( "f1", "[0]+[1]*x", fxFitMin, fxFitMax);
		f1->SetLineColor(kRed);
		f1->SetParameters(aCoef,bCoef); 
	}
	else if ( histType == 1 )	{
		f1 = new TF1( "f1", "[0]+[1]*(x-1)", fxFitMin/mnf, fxFitMax/mnf);
		f1->SetLineColor(kGreen);
		f1->SetParameters(aCoef,bCoef); 
	}
	//cout << " set draw params: a=" << aCoef << " b=" << bCoef << endl;

	y1=profToDraw->GetBinContent(1)-profToDraw->GetBinError(1);
	y2=profToDraw->GetBinContent(1)+profToDraw->GetBinError(1);
	n=profToDraw->GetNbinsX();

	for(i=2; i<=n; i++)
	{
        	if(profToDraw->GetBinContent(i)-profToDraw->GetBinError(i)<y1)
            		y1=profToDraw->GetBinContent(i)-profToDraw->GetBinError(i);
        	if(profToDraw->GetBinContent(i)+profToDraw->GetBinError(i)>y2)
            		y2=profToDraw->GetBinContent(i)+profToDraw->GetBinError(i);
	}

	profToDraw->DrawCopy();
	if(f1) f1->DrawCopy("same");
	else return;
	
	x1 = profToDraw->GetXaxis()->GetXmin();
	x2 = profToDraw->GetXaxis()->GetXmax();

	if ( drawPaveLabel )
	{	
		int nDatas = 0;
		for ( int j = 0; j < 9; j++)
			if ( mDrawArray[j] ) nDatas++;
		double aXshift = (x2-x1)/7;
		double aYshift = (y2-y1)/20;

		TPaveText *pt1 = new TPaveText(x1+(x2-x1)/2 + aXshift, y1+aYshift, x2-(x2-x1)/6 + aXshift, y1+(y2-y1)/3*2/9*nDatas + aYshift);
		snprintf(str,50, "Events = %i", fEntries);
		//sprintf(str, "Events = %i", fEntries);
		if ( mDrawArray[0] ) pt1->AddText(str);
		snprintf(str, 50,"a = %g #pm %g", GetRoundWithError( aCoef, aCoefError ), GetRoundWithPrecision(aCoefError, 2)); //fit1->GetParameter(0), fit1->GetParError(0));
		//sprintf(str, "a = %g #pm %g", GetRoundWithError( aCoef, aCoefError ), GetRoundWithPrecision(aCoefError, 2)); //fit1->GetParameter(0), fit1->GetParError(0));
		if ( mDrawArray[1] ) pt1->AddText(str);
		snprintf(str, 50,"b = %g #pm %g", GetRoundWithError( bCoef, bCoefError ), GetRoundWithPrecision(bCoefError, 2)); //fit1->GetParameter(1), fit1->GetParError(1));
		//sprintf(str, "b = %g #pm %g", GetRoundWithError( bCoef, bCoefError ), GetRoundWithPrecision(bCoefError, 2)); //fit1->GetParameter(1), fit1->GetParError(1));
		if ( mDrawArray[2] ) pt1->AddText(str);
		snprintf(str, 50,"#hat{#chi}^{2} = #chi^{2}/(n-2) = %g", GetRoundWithPrecision(fXi2abs, 3));
		//sprintf(str, "#hat{#chi}^{2} = #chi^{2}/(n-2) = %g", GetRoundWithPrecision(fXi2abs, 3));
		if ( mDrawArray[3] ) pt1->AddText(str);
		snprintf(str,50, "<%s> = %g " , fSx, GetRoundWithPrecision(fPrf->GetMean(), 3));
		//sprintf(str, "<%s> = %g " , fSx, GetRoundWithPrecision(fPrf->GetMean(), 3));
		if ( mDrawArray[4] ) pt1->AddText(str);
		
		snprintf(str,50, "<%s> = %g", fSy,  GetRoundWithPrecision(fPrb->GetMean(),3));
		//sprintf(str, "<%s> = %g", fSy,  GetRoundWithPrecision(fPrb->GetMean(),3));
		if ( mDrawArray[5] ) pt1->AddText(str);
		
		snprintf(str, 50,"<<%s>> = %g " , fSx, GetRoundWithPrecision(fPrf->GetRMS(), 3));
		//sprintf(str, "<<%s>> = %g " , fSx, GetRoundWithPrecision(fPrf->GetRMS(), 3));
		if ( mDrawArray[6] ) pt1->AddText(str);
		snprintf(str, 50,"<<%s>> = %g", fSy,  GetRoundWithPrecision(fPrb->GetRMS(), 3));
		//sprintf(str, "<<%s>> = %g", fSy,  GetRoundWithPrecision(fPrb->GetRMS(), 3));
		if ( mDrawArray[7] ) pt1->AddText(str);
		
		if ( fdptb ) {
			snprintf(str,50, "d%s = %g", fSy,  GetRoundWithPrecision(fdptb, 3));
			//sprintf(str, "d%s = %g", fSy,  GetRoundWithPrecision(fdptb, 3));
			if ( mDrawArray[8] ) pt1->AddText(str);
		}
		
		pt1->SetTextAlign(12);
		pt1->SetTextFont(42);
		pt1->SetFillColor(4000);
		//pt1->SetFillStyle(4100);
		pt1->SetShadowColor(4000);
		pt1->SetBorderSize(0);
		
		pt1->DrawClone("same");
	}
}

void AliLRCAnalysis::SetGraphics() const {
// Set root graph style
	TStyle tempSt;
    gStyle->SetCanvasColor(10);
    gStyle->SetFrameFillColor(10);
    gStyle->SetStatColor(0);
    gStyle->SetPadColor(0);
}

double AliLRCAnalysis::Integral(TH2D* const source, Int_t nbin) const {
// calculate the integrall for x bin and y bins of 2d histogramm
    double sum = 0;
    for(Int_t i = 1; i<=source->GetNbinsY(); i++) {
        sum += source->GetYaxis()->GetBinCenter(i)*source->GetBinContent(nbin, i);
    }
    return sum;
}

void AliLRCAnalysis::SetErrors(TH2D* source, const char *name){
// Calculate errors for NN
	TProfile* profX = (TProfile*) source->ProfileX(name, 1, source->GetNbinsY());
	for(int i = 0; i < profX->GetNbinsX(); i++)		
	{
		fPrAbs->SetBinContent(i, profX->GetBinContent(i));
		if(fPrf->GetBinContent(i)!=0)
		{
		   fPrAbs->SetBinError(i,sqrt(profX->GetBinContent(i)/fPrf->GetBinContent(i)));
		}
		fPrRel->SetBinContent(i, fPrAbs->GetBinContent(i)/fPrb->GetMean());
		fPrRel->SetBinError(i,fPrAbs->GetBinError(i)/fPrb->GetMean());	
	}	
}

void AliLRCAnalysis::SetErrors(TH2D* source, const char *name, double ptd, TH2D* nb){
//Calculate arrors for ptn and ptpt
	TProfile* profX = (TProfile*) source->ProfileX(name, 1, source->GetNbinsY());
	fdptb = ptd;
	double pt;
	for(int i = 0; i < profX->GetNbinsX(); i++)		
	{
		fPrAbs->SetBinContent(i, profX->GetBinContent(i));
		if(fPrf->GetBinContent(i)!=0)
		{
			  pt = profX->GetBinContent(i);
			  fPrAbs->SetBinError(i,ptd*sqrt(Integral(nb,i))/fPrf->GetBinContent(i));
		}
		fPrRel->SetBinContent(i, fPrAbs->GetBinContent(i)/fPrb->GetMean());
		fPrRel->SetBinError(i,fPrAbs->GetBinError(i)/fPrb->GetMean());	
	}	

}

void AliLRCAnalysis::SetErrors(TH2D* source, const char *name, double ptd, const TProfile* nb){
//Calculate arrors for ptn and ptpt
	TProfile* profX = (TProfile*) source->ProfileX(name, 1, source->GetNbinsY());
	fdptb = ptd;
	double pt;
	for(int i = 0; i < profX->GetNbinsX(); i++)		
	{
		fPrAbs->SetBinContent(i, profX->GetBinContent(i));
		if(fPrf->GetBinContent(i)!=0)
		{
			  pt = profX->GetBinContent(i);
			  fPrAbs->SetBinError(i,ptd*sqrt(nb->GetBinContent(i)/fPrf->GetBinContent(i)));
		}
		fPrRel->SetBinContent(i, fPrAbs->GetBinContent(i)/fPrb->GetMean());
		fPrRel->SetBinError(i,fPrAbs->GetBinError(i)/fPrb->GetMean());	
	}	

}

double AliLRCAnalysis::GetRoundWithError( double value, double error ) const {
//Rounding error and value with DEFAULT precision
	return GetRoundValueErrorPrecision( value, error, 2 );
}
double AliLRCAnalysis::GetRoundWithError( double value, double error, int pres ) const {
//Rounding error and value with REQUESTED precision
	return GetRoundValueErrorPrecision( value, error, pres );
}
double AliLRCAnalysis::GetRoundWithPrecision( double value, int pres ) const {
//Rounding error and value with requested precision
	return GetRoundValueErrorPrecision( value, 0, pres );
}

double AliLRCAnalysis::GetRoundValueErrorPrecision( double value, double error, int pres ) const {
	//Rounding error and value with requested precision
	//value == value, error == error
	//if single argument(value=0) - calculate without errors:

	int i = 0;
	double order = 1;
	bool noError = false;
	pres -= 1;
	
	//cout << "Before rounding: " << value << " " << error << endl;

	if ( !error)
	{
		error = value;
		noError = true;
	}
	while ( ((int)error)%10 == 0 && i < 6 )
	{
		i++;
		error*=10;
		order*=10;
	}

	for ( int j = 0; j < pres; j++ )
	{
		error*=10;
		order*=10;
		i++;
	}
	int adding = 0;
	if ( ((int)(error*10))%10 > 4 && ((int)(error*10))%10 != 9 ) //trouble: if we round 19 to 20 - zero disappeares!
		adding = 1;
	error = (double)((int)error + adding)/order;
	
	if ( noError )
	{
		//cout << "After rounding: " << error << endl;
		return error;
	}
	else
	{
		for ( int j = 0; j < i; j++ )
			value*=10;
		
		adding = 0;
		if ( ((int)(value*10))%10 > 4 && ((int)(value*10))%10 != 9 )
			adding = 1;
		value = (double)((int)value + adding)/order;
		
		//cout << "After rounding: " << value << " " << error << endl;
		return value; //taking into account ERROR
	}
}

TH1D* AliLRCAnalysis::GetAbsHisto() const
{
//Returns final histo in absolute variables
return fPrAbs;
}
TH1D* AliLRCAnalysis::GetRelHisto() const
{
//Returns final histo in relative variables
return fPrRel;
}

TH1D* AliLRCAnalysis::GetForwardValueDist() const
{
//Returns destribution of value used in forward window
return fPrf;
}
TH1D* AliLRCAnalysis::GetBackwardValueDist() const
{
//Returns destribution of value used in backward window
return fPrb;
}
Commit	Line	Data
d03e5de4	1	/**************************************************************************
ec6ab54b	2	* Author: Andrey Ivanov. *
d03e5de4	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	// Description:
	15	// This class include into LRC library for Long-Range Correlation analysis
	16	// it is base class for NN, PtN, PtPt
	17	// implements base methods for thees classes
ec6ab54b	18	// Origin: Petr Naumenko, SPbSU-CERN, Petr.Naoumenko@cern.ch ,
ec6ab54b	19	// Andrey Ivanov (SPbSU-CERN), Igor Altsebeev (SPbSU-CERN)
d03e5de4	20	//-------------------------------------------------------------------------
	21
	22	#include "AliLRCAnalysis.h"
65e5de6a	23	#include "TFile.h"
	24	#include "AliLRCFit.h"
	25	#include "TProfile.h"
	26	#include "TH1D.h"
	27	#include "TH2D.h"
	28	#include "TPaveText.h"
	29	#include "TF1.h"
ec6ab54b	30	#include "math.h"
65e5de6a	31	#include "TStyle.h"
65e5de6a	32	class gStyle;
ec6ab54b	33
d03e5de4	34
	35	ClassImp(AliLRCAnalysis)
	36
	37	/******************************************************
	38	* AliLRCAnalysis class
	39	******************************************************/
	40
ec6ab54b	41	AliLRCAnalysis::AliLRCAnalysis():TObject(), fPrAbs(new TH1D()), fPrRel(new TH1D()), fPrf(new TH1D()), fPrb(new TH1D()), fileHist(new TFile()), fdptb(.0), fEntries(0), fSx((char)" "), fSy((char)" "), fxFitMin(.0), fxFitMax(.0), fNsigma(2.), farel(.0), fbrel(.0), farelError(.0), fbrelError(.0), fXi2rel(.0), faabs(.0), fbabs(.0), faabsError(.0), fbabsError(.0), fXi2abs(.0),fFitMethod(0){
d03e5de4	42	//Empty constructor
	43	}
	44
ec6ab54b	45	AliLRCAnalysis::AliLRCAnalysis(const AliLRCAnalysis& a):TObject(a),fPrAbs(a.fPrAbs), fPrRel(a.fPrRel), fPrf(a.fPrf), fPrb(a.fPrb), fileHist(a.fileHist), fdptb(a.fdptb), fEntries(a.fEntries), fSx(a.fSx), fSy(a.fSy), fxFitMin(a.fxFitMin), fxFitMax(a.fxFitMax), fNsigma(a.fNsigma), farel(a.farel), fbrel(a.fbrel), farelError(a.farelError), fbrelError(a.fbrelError), fXi2rel(a.fXi2rel), faabs(a.faabs), fbabs(a.fbabs),faabsError(a.faabsError), fbabsError(a.fbabsError), fXi2abs(a.fXi2abs),fFitMethod(a.fFitMethod){
d03e5de4	46	//Constructor
	47	}
	48
	49	AliLRCAnalysis& AliLRCAnalysis::operator= (const AliLRCAnalysis& a){
	50	//Operator =
	51	if(this!=&a){
ec6ab54b	52	TObject::operator= (a);
d03e5de4	53	fPrAbs = a.fPrAbs;
	54	fPrRel = a.fPrRel;
	55	fPrf = a.fPrf;
	56	fPrb = a.fPrb;
	57	fileHist = a.fileHist;
	58	fSx = a.fSx;
	59	fSy = a.fSy;
	60	fdptb = a.fdptb;
	61	fEntries = a.fEntries;
d96e5666	62	fxFitMin = a.fxFitMin;
d96e5666	63	fxFitMax = a.fxFitMax;
ec6ab54b	64	fNsigma = a.fNsigma;
65e5de6a	65	farel = a.farel;
65e5de6a	66	fbrel = a.fbrel;
adf00ba6	67	farelError = a.farelError;
ec6ab54b	68	fbrelError = a.fbrelError;
65e5de6a	69	fXi2rel = a.fXi2rel;
	70	faabs = a.faabs;
	71	fbabs = a.fbabs;
adf00ba6	72	faabsError = a.faabsError;
ec6ab54b	73	fbabsError = a.fbabsError;
65e5de6a	74	fXi2abs = a.fXi2abs;
ec6ab54b	75	fFitMethod = a.fFitMethod;
d03e5de4	76	}
	77	return *this;
	78	}
	79
d03e5de4	80	AliLRCAnalysis::~AliLRCAnalysis() {
	81	//Destructor
	82	delete fPrAbs;
	83	delete fPrRel;
	84	delete fileHist;
65e5de6a	85	delete fPrf;
65e5de6a	86	delete fPrb;
d03e5de4	87	}
d03e5de4	88
d03e5de4	89
	90
	91	void AliLRCAnalysis::CreateHist(char name, char nameAbs, char nameRel, char atitleF, char atitleB, char rtitleF, char rtitleB, TH2D sourceHist) {
	92	//Create absolute and relation var histogramm
	93	TProfile* profX = (TProfile*) sourceHist->ProfileX(name, 1, sourceHist->GetNbinsY());
	94	fEntries = (int) sourceHist->GetEntries();
	95	fPrAbs = new TH1D(nameAbs, profX->GetTitle(), profX->GetXaxis()->GetNbins(), profX->GetXaxis()->GetXmin(), profX->GetXaxis()->GetXmax());
	96	fPrAbs->SetOption("E");
	97	fPrRel = new TH1D(nameRel, profX->GetTitle(), profX->GetXaxis()->GetNbins(), profX->GetXaxis()->GetXmin()/sourceHist->ProjectionX()->GetMean(), profX->GetXaxis()->GetXmax()/sourceHist->ProjectionX()->GetMean());
	98	fPrRel->SetOption("E");
	99	fPrAbs->GetXaxis()->SetTitle(atitleF);
	100	fPrAbs->GetYaxis()->SetTitle(atitleB);
	101	fPrRel->GetXaxis()->SetTitle(rtitleF);
	102	fPrRel->GetYaxis()->SetTitle(rtitleB);
	103	fPrf = (TH1D*) sourceHist->ProjectionX();
	104	fPrb = (TH1D*) sourceHist->ProjectionY();
	105	fPrf->GetXaxis()->SetTitle(atitleF);
	106	fPrf->GetYaxis()->SetTitle("Tracks");
	107	fPrb->GetXaxis()->SetTitle(atitleB);
	108	fPrb->GetYaxis()->SetTitle("Tracks");
	109	fSx = atitleF;
	110	fSy = atitleB;
d96e5666	111	double mnf = fPrf->GetMean();
ec6ab54b	112	double rms = fPrf->GetRMS();
	113	fxFitMin = mnf-fNsigma*rms;//sqrt(rms);
	114	fxFitMax = mnf+fNsigmarms;//sqrt(rms);//mnf+2sqrt(rms);
65e5de6a	115	//delete profX;
d96e5666	116
	117	}
	118
	119	void AliLRCAnalysis::SetBinsRange(int binMin, int binMax){
65e5de6a	120	//Set the bin range
d96e5666	121	TH1D* h=fPrf;
ec6ab54b	122	int n=h->GetNbinsX();
d96e5666	123	fxFitMin = h->GetXaxis()->GetXmin();
d96e5666	124	fxFitMax = h->GetXaxis()->GetXmax();
65e5de6a	125	double df = (fxFitMax-fxFitMin)/n;
65e5de6a	126	for(int i=1; i<=n; i++)
d96e5666	127	if(h->GetBinContent(i) != 0){
	128	fxFitMin += (i + binMin) * df;
	129	break;
	130	}
65e5de6a	131	for(int i=1; i<=n; i++)
65e5de6a	132	if(h->GetBinContent(n-i) != 0){
d96e5666	133	fxFitMax -= (i + binMax) * df;
	134	break;
	135	}
	136	}
	137
	138	bool AliLRCAnalysis::SetFitRange(double xMin, double xMax){
65e5de6a	139	//Set the fit range
d96e5666	140	if(xMax < xMin){
	141	return false;
	142	}
	143	this->fxFitMin = xMin;
	144	this->fxFitMax = xMax;
	145	return true;
	146	}
ec6ab54b	147	bool AliLRCAnalysis::SetFitRangeMin(double xMin){
	148	//Set the fit range min
	149	fxFitMin = xMin;
	150	return true;
	151	}
	152
	153	bool AliLRCAnalysis::SetFitRangeMax(double xMax){
	154	//Set the fit range min
	155	//this->fxFitMax = xMax;
	156
	157	TH1D* h=fPrf;
	158	double maxBorder = 0.;
	159
	160	for ( int binI = h->GetNbinsX(); binI > 0; binI-- )
	161	{
	162	if ( h->GetBinContent(binI) != 0 )
	163	{
	164	maxBorder = h->GetBinLowEdge(binI+1) ;
	165	break;
	166	}
	167	}
	168	fxFitMax = TMath::Max( fxFitMin, TMath::Min( xMax, maxBorder ) );
	169
	170	return true;
	171	}
	172
65e5de6a	173	void AliLRCAnalysis::SetFullFitRange(){
	174	//Set fitting on full range
	175	TH1D* h=fPrf;
	176	fxFitMin = h->GetXaxis()->GetXmin();
65e5de6a	177
	178	for ( int binI = h->GetNbinsX(); binI > 0; binI-- )
	179	{
	180	if ( h->GetBinContent(binI) != 0 )
	181	{
	182	fxFitMax = h->GetBinLowEdge(binI+1) ;
	183	break;
	184	}
	185	}
	186	}
	187
d96e5666	188
	189	void AliLRCAnalysis::SetXmin(double xMin){
	190	fxFitMin = xMin;
	191	}
	192
	193	void AliLRCAnalysis::SetXmax(double xMax){
	194	fxFitMax = xMax;
	195	}
ec6ab54b	196	void AliLRCAnalysis::SetNsigma(double nSigma){
	197	//Sets fiting range in forward value RMSs
	198	fNsigma = nSigma;
	199	double mnf = fPrf->GetMean();
	200	double rms = fPrf->GetRMS();
	201	fxFitMin = mnf-fNsigma*rms;
	202	fxFitMax = mnf+fNsigma*rms;
d96e5666	203
ec6ab54b	204	// cout << "### Fit params: N sigma is " << fNsigma
	205	// << ", mean = " << mnf << ", rms = " << rms << ", minFitX = " << fxFitMin << ", maxFitX = " << fxFitMax << endl;
	206	}
	207
	208	void AliLRCAnalysis::SetFitMethod(int id){
	209	//Choose fit method
	210	fFitMethod = id;
	211	}
	212	double AliLRCAnalysis::GetFitXmin() const
	213	{
	214	return fxFitMin;
	215	}
	216
	217	double AliLRCAnalysis::GetFitXmax() const
	218	{
	219	return fxFitMax;
	220	}
65e5de6a	221	double AliLRCAnalysis::GetArel() const {
65e5de6a	222	return farel;
d96e5666	223	}
d96e5666	224
65e5de6a	225	double AliLRCAnalysis::GetBrel() const {
65e5de6a	226	return fbrel;
d96e5666	227	}
d96e5666	228
65e5de6a	229	double AliLRCAnalysis::GetArelError() const {
65e5de6a	230	return farelError;
d96e5666	231	}
d96e5666	232
65e5de6a	233	double AliLRCAnalysis::GetBrelError() const {
65e5de6a	234	return fbrelError;
d96e5666	235	}
d96e5666	236
65e5de6a	237	double AliLRCAnalysis::GetXi2rel() const {
65e5de6a	238	return fXi2rel;
d96e5666	239	}
d96e5666	240
65e5de6a	241	double AliLRCAnalysis::GetAabs() const {
	242	return faabs;
	243	}
	244
	245	double AliLRCAnalysis::GetBabs() const {
	246	return fbabs;
	247	}
	248
	249	double AliLRCAnalysis::GetAabsError() const {
	250	return faabsError;
	251	}
	252
	253	double AliLRCAnalysis::GetBabsError() const {
	254	return fbabsError;
	255	}
	256	double AliLRCAnalysis::GetXi2abs() const {
	257	return fXi2abs;
d96e5666	258	}
	259
	260	void AliLRCAnalysis::Calculate(){
65e5de6a	261	//Calculate all
ec6ab54b	262	double mnf = fPrf->GetMean();
	263
	264	// new fit for abs
	265	AliLRCFit *fit1 = new AliLRCFit(fPrRel, fxFitMin/mnf, fxFitMax/mnf, 1.);
	266	AliLRCFit *fit2 = new AliLRCFit(fPrAbs, fxFitMin, fxFitMax, 0.);
65e5de6a	267	farel = fit1->Geta();
65e5de6a	268	fbrel = fit1->Getb();
ec6ab54b	269
	270	faabs = fit2->Geta();
	271	fbabs = fit2->Getb();
	272
	273	if ( fFitMethod==0 )
	274	{
	275	faabsError = fit2->Getda();
	276	fbabsError = fit2->Getdb();
	277
	278	farelError = fit1->Getda();
	279	fbrelError = fit1->Getdb();
	280	}
	281	if ( fFitMethod==1 )
	282	{
	283	faabsError = fit2->Getda1();
	284	fbabsError = fit2->Getdb1();
	285
	286	farelError = fit1->Getda1();
	287	fbrelError = fit1->Getdb1();
	288	}
	289
65e5de6a	290	fXi2rel = fit1->Gethi2();
ec6ab54b	291	fXi2abs = fit2->Gethi2();
d03e5de4	292	}
	293
	294	void AliLRCAnalysis::DrawAbs() {
adf00ba6	295	//Draw abs var hist with ALL info
ec6ab54b	296	int * mas = new int [fgkPlotFlags];
ec6ab54b	297	for ( int i = 0; i < fgkPlotFlags; i++ )
adf00ba6	298	mas[i] = 1;
	299	DrawAbsPure( mas, 1 );
	300	delete []mas;
65e5de6a	301	}
65e5de6a	302
ec6ab54b	303	void AliLRCAnalysis::DrawAbs( const int * const mas ) {
65e5de6a	304	//Draw abs var hist with REQUESTED BY ARRAY info
	305	DrawAbsPure( mas, 1 );
	306	}
	307
	308	void AliLRCAnalysis::DrawAbsPure( const int * const mDrawArray, bool drawPaveLabel ) {
ec6ab54b	309	Calculate();
65e5de6a	310	// Draw abs var histrogram
ec6ab54b	311	DrawHist( mDrawArray, drawPaveLabel, faabs, fbabs, faabsError, fbabsError, fPrAbs, 0 );
d03e5de4	312	}
	313
	314	void AliLRCAnalysis::DrawRel() {
adf00ba6	315	//Draw rel var hist with ALL info
ec6ab54b	316	int * mas = new int [fgkPlotFlags];
ec6ab54b	317	for ( int i = 0; i < fgkPlotFlags; i++ )
adf00ba6	318	mas[i] = 1;
	319	DrawRelPure( mas, 1 );
	320	delete []mas;
65e5de6a	321	}
65e5de6a	322
ec6ab54b	323	void AliLRCAnalysis::DrawRel( const int * const mas ) {
65e5de6a	324	//Draw rel var hist with REQUESTED BY ARRAY info
	325	DrawRelPure( mas, 1 );
	326	}
	327
	328	void AliLRCAnalysis::DrawRelPure( const int * const mDrawArray, bool drawPaveLabel ) {
ec6ab54b	329	Calculate();
65e5de6a	330	// Draw rel var histogram
ec6ab54b	331	DrawHist( mDrawArray, drawPaveLabel, farel, fbrel, farelError, fbrelError, fPrRel, 1 );
	332	}
	333	void AliLRCAnalysis::DrawHist( const int * const mDrawArray, bool drawPaveLabel, double aCoef, double bCoef,
	334	double aCoefError, double bCoefError, TH1D* profToDraw, int histType ) {
	335	// Method called by DrawRelPure or DrawAbsPure to draw corresponding LRC histo
	336	double mnf = fPrf->GetMean();
65e5de6a	337	double y1, y2, x1, x2;
	338	Int_t i, n;
	339	char str[50];
ec6ab54b	340	//mnf = fPrf->GetMean();
	341	profToDraw->SetStats(0);
	342	//profToDraw->Fit("pol1", "", "", fxFitMin, fxFitMax);
	343	//profToDraw->Fit("pol1");
	344	//TF1 *fit1 = profToDraw->GetFunction("pol1");
	345
	346	//draw fit line
c511520f	347	TF1 *f1 = 0;
ec6ab54b	348	if ( histType == 0 ) {
	349	f1 = new TF1( "f1", "[0]+[1]*x", fxFitMin, fxFitMax);
	350	f1->SetLineColor(kRed);
99cdfc57	351	f1->SetParameters(aCoef,bCoef);
ec6ab54b	352	}
	353	else if ( histType == 1 ) {
	354	f1 = new TF1( "f1", "[0]+[1]*(x-1)", fxFitMin/mnf, fxFitMax/mnf);
	355	f1->SetLineColor(kGreen);
99cdfc57	356	f1->SetParameters(aCoef,bCoef);
ec6ab54b	357	}
ec6ab54b	358	//cout << " set draw params: a=" << aCoef << " b=" << bCoef << endl;
	359
	360	y1=profToDraw->GetBinContent(1)-profToDraw->GetBinError(1);
	361	y2=profToDraw->GetBinContent(1)+profToDraw->GetBinError(1);
	362	n=profToDraw->GetNbinsX();
	363
	364	for(i=2; i<=n; i++)
	365	{
	366	if(profToDraw->GetBinContent(i)-profToDraw->GetBinError(i)<y1)
	367	y1=profToDraw->GetBinContent(i)-profToDraw->GetBinError(i);
	368	if(profToDraw->GetBinContent(i)+profToDraw->GetBinError(i)>y2)
	369	y2=profToDraw->GetBinContent(i)+profToDraw->GetBinError(i);
	370	}
65e5de6a	371
ec6ab54b	372	profToDraw->DrawCopy();
c511520f	373	if(f1) f1->DrawCopy("same");
c511520f	374	else return;
65e5de6a	375
ec6ab54b	376	x1 = profToDraw->GetXaxis()->GetXmin();
	377	x2 = profToDraw->GetXaxis()->GetXmax();
	378
65e5de6a	379	if ( drawPaveLabel )
	380	{
	381	int nDatas = 0;
	382	for ( int j = 0; j < 9; j++)
	383	if ( mDrawArray[j] ) nDatas++;
	384	double aXshift = (x2-x1)/7;
	385	double aYshift = (y2-y1)/20;
ec6ab54b	386
65e5de6a	387	TPaveText pt1 = new TPaveText(x1+(x2-x1)/2 + aXshift, y1+aYshift, x2-(x2-x1)/6 + aXshift, y1+(y2-y1)/32/9*nDatas + aYshift);
99cdfc57	388	snprintf(str,50, "Events = %i", fEntries);
99cdfc57	389	//sprintf(str, "Events = %i", fEntries);
65e5de6a	390	if ( mDrawArray[0] ) pt1->AddText(str);
99cdfc57	391	snprintf(str, 50,"a = %g #pm %g", GetRoundWithError( aCoef, aCoefError ), GetRoundWithPrecision(aCoefError, 2)); //fit1->GetParameter(0), fit1->GetParError(0));
99cdfc57	392	//sprintf(str, "a = %g #pm %g", GetRoundWithError( aCoef, aCoefError ), GetRoundWithPrecision(aCoefError, 2)); //fit1->GetParameter(0), fit1->GetParError(0));
65e5de6a	393	if ( mDrawArray[1] ) pt1->AddText(str);
99cdfc57	394	snprintf(str, 50,"b = %g #pm %g", GetRoundWithError( bCoef, bCoefError ), GetRoundWithPrecision(bCoefError, 2)); //fit1->GetParameter(1), fit1->GetParError(1));
99cdfc57	395	//sprintf(str, "b = %g #pm %g", GetRoundWithError( bCoef, bCoefError ), GetRoundWithPrecision(bCoefError, 2)); //fit1->GetParameter(1), fit1->GetParError(1));
65e5de6a	396	if ( mDrawArray[2] ) pt1->AddText(str);
99cdfc57	397	snprintf(str, 50,"#hat{#chi}^{2} = #chi^{2}/(n-2) = %g", GetRoundWithPrecision(fXi2abs, 3));
99cdfc57	398	//sprintf(str, "#hat{#chi}^{2} = #chi^{2}/(n-2) = %g", GetRoundWithPrecision(fXi2abs, 3));
65e5de6a	399	if ( mDrawArray[3] ) pt1->AddText(str);
99cdfc57	400	snprintf(str,50, "<%s> = %g " , fSx, GetRoundWithPrecision(fPrf->GetMean(), 3));
99cdfc57	401	//sprintf(str, "<%s> = %g " , fSx, GetRoundWithPrecision(fPrf->GetMean(), 3));
65e5de6a	402	if ( mDrawArray[4] ) pt1->AddText(str);
65e5de6a	403
99cdfc57	404	snprintf(str,50, "<%s> = %g", fSy, GetRoundWithPrecision(fPrb->GetMean(),3));
99cdfc57	405	//sprintf(str, "<%s> = %g", fSy, GetRoundWithPrecision(fPrb->GetMean(),3));
65e5de6a	406	if ( mDrawArray[5] ) pt1->AddText(str);
65e5de6a	407
99cdfc57	408	snprintf(str, 50,"<<%s>> = %g " , fSx, GetRoundWithPrecision(fPrf->GetRMS(), 3));
99cdfc57	409	//sprintf(str, "<<%s>> = %g " , fSx, GetRoundWithPrecision(fPrf->GetRMS(), 3));
65e5de6a	410	if ( mDrawArray[6] ) pt1->AddText(str);
99cdfc57	411	snprintf(str, 50,"<<%s>> = %g", fSy, GetRoundWithPrecision(fPrb->GetRMS(), 3));
99cdfc57	412	//sprintf(str, "<<%s>> = %g", fSy, GetRoundWithPrecision(fPrb->GetRMS(), 3));
65e5de6a	413	if ( mDrawArray[7] ) pt1->AddText(str);
	414
	415	if ( fdptb ) {
99cdfc57	416	snprintf(str,50, "d%s = %g", fSy, GetRoundWithPrecision(fdptb, 3));
99cdfc57	417	//sprintf(str, "d%s = %g", fSy, GetRoundWithPrecision(fdptb, 3));
65e5de6a	418	if ( mDrawArray[8] ) pt1->AddText(str);
65e5de6a	419	}
ec6ab54b	420
65e5de6a	421	pt1->SetTextAlign(12);
	422	pt1->SetTextFont(42);
	423	pt1->SetFillColor(4000);
	424	//pt1->SetFillStyle(4100);
	425	pt1->SetShadowColor(4000);
	426	pt1->SetBorderSize(0);
	427
ec6ab54b	428	pt1->DrawClone("same");
65e5de6a	429	}
d03e5de4	430	}
	431
	432	void AliLRCAnalysis::SetGraphics() const {
	433	// Set root graph style
65e5de6a	434	TStyle tempSt;
d03e5de4	435	gStyle->SetCanvasColor(10);
	436	gStyle->SetFrameFillColor(10);
	437	gStyle->SetStatColor(0);
	438	gStyle->SetPadColor(0);
	439	}
	440
65e5de6a	441	double AliLRCAnalysis::Integral(TH2D* const source, Int_t nbin) const {
d03e5de4	442	// calculate the integrall for x bin and y bins of 2d histogramm
	443	double sum = 0;
	444	for(Int_t i = 1; i<=source->GetNbinsY(); i++) {
	445	sum += source->GetYaxis()->GetBinCenter(i)*source->GetBinContent(nbin, i);
	446	}
	447	return sum;
	448	}
	449
	450	void AliLRCAnalysis::SetErrors(TH2D* source, const char *name){
	451	// Calculate errors for NN
	452	TProfile* profX = (TProfile*) source->ProfileX(name, 1, source->GetNbinsY());
	453	for(int i = 0; i < profX->GetNbinsX(); i++)
	454	{
	455	fPrAbs->SetBinContent(i, profX->GetBinContent(i));
	456	if(fPrf->GetBinContent(i)!=0)
	457	{
ec6ab54b	458	fPrAbs->SetBinError(i,sqrt(profX->GetBinContent(i)/fPrf->GetBinContent(i)));
d03e5de4	459	}
	460	fPrRel->SetBinContent(i, fPrAbs->GetBinContent(i)/fPrb->GetMean());
	461	fPrRel->SetBinError(i,fPrAbs->GetBinError(i)/fPrb->GetMean());
	462	}
	463	}
	464
	465	void AliLRCAnalysis::SetErrors(TH2D* source, const char name, double ptd, TH2D nb){
	466	//Calculate arrors for ptn and ptpt
	467	TProfile* profX = (TProfile*) source->ProfileX(name, 1, source->GetNbinsY());
	468	fdptb = ptd;
	469	double pt;
	470	for(int i = 0; i < profX->GetNbinsX(); i++)
	471	{
	472	fPrAbs->SetBinContent(i, profX->GetBinContent(i));
	473	if(fPrf->GetBinContent(i)!=0)
	474	{
	475	pt = profX->GetBinContent(i);
ec6ab54b	476	fPrAbs->SetBinError(i,ptd*sqrt(Integral(nb,i))/fPrf->GetBinContent(i));
d03e5de4	477	}
	478	fPrRel->SetBinContent(i, fPrAbs->GetBinContent(i)/fPrb->GetMean());
	479	fPrRel->SetBinError(i,fPrAbs->GetBinError(i)/fPrb->GetMean());
	480	}
	481
	482	}
	483
ec6ab54b	484	void AliLRCAnalysis::SetErrors(TH2D* source, const char name, double ptd, const TProfile nb){
d96e5666	485	//Calculate arrors for ptn and ptpt
	486	TProfile* profX = (TProfile*) source->ProfileX(name, 1, source->GetNbinsY());
	487	fdptb = ptd;
	488	double pt;
	489	for(int i = 0; i < profX->GetNbinsX(); i++)
	490	{
	491	fPrAbs->SetBinContent(i, profX->GetBinContent(i));
	492	if(fPrf->GetBinContent(i)!=0)
	493	{
	494	pt = profX->GetBinContent(i);
ec6ab54b	495	fPrAbs->SetBinError(i,ptd*sqrt(nb->GetBinContent(i)/fPrf->GetBinContent(i)));
d96e5666	496	}
	497	fPrRel->SetBinContent(i, fPrAbs->GetBinContent(i)/fPrb->GetMean());
	498	fPrRel->SetBinError(i,fPrAbs->GetBinError(i)/fPrb->GetMean());
	499	}
	500
	501	}
	502
ec6ab54b	503	double AliLRCAnalysis::GetRoundWithError( double value, double error ) const {
65e5de6a	504	//Rounding error and value with DEFAULT precision
	505	return GetRoundValueErrorPrecision( value, error, 2 );
	506	}
ec6ab54b	507	double AliLRCAnalysis::GetRoundWithError( double value, double error, int pres ) const {
65e5de6a	508	//Rounding error and value with REQUESTED precision
	509	return GetRoundValueErrorPrecision( value, error, pres );
	510	}
ec6ab54b	511	double AliLRCAnalysis::GetRoundWithPrecision( double value, int pres ) const {
65e5de6a	512	//Rounding error and value with requested precision
	513	return GetRoundValueErrorPrecision( value, 0, pres );
	514	}
	515
	516	double AliLRCAnalysis::GetRoundValueErrorPrecision( double value, double error, int pres ) const {
	517	//Rounding error and value with requested precision
	518	//value == value, error == error
	519	//if single argument(value=0) - calculate without errors:
	520
	521	int i = 0;
	522	double order = 1;
	523	bool noError = false;
	524	pres -= 1;
	525
ec6ab54b	526	//cout << "Before rounding: " << value << " " << error << endl;
65e5de6a	527
	528	if ( !error)
	529	{
	530	error = value;
	531	noError = true;
	532	}
ec6ab54b	533	while ( ((int)error)%10 == 0 && i < 6 )
65e5de6a	534	{
	535	i++;
	536	error*=10;
	537	order*=10;
	538	}
	539
	540	for ( int j = 0; j < pres; j++ )
	541	{
	542	error*=10;
	543	order*=10;
	544	i++;
	545	}
	546	int adding = 0;
	547	if ( ((int)(error10))%10 > 4 && ((int)(error10))%10 != 9 ) //trouble: if we round 19 to 20 - zero disappeares!
	548	adding = 1;
	549	error = (double)((int)error + adding)/order;
	550
	551	if ( noError )
	552	{
ec6ab54b	553	//cout << "After rounding: " << error << endl;
65e5de6a	554	return error;
	555	}
	556	else
	557	{
	558	for ( int j = 0; j < i; j++ )
	559	value*=10;
	560
	561	adding = 0;
	562	if ( ((int)(value10))%10 > 4 && ((int)(value10))%10 != 9 )
	563	adding = 1;
	564	value = (double)((int)value + adding)/order;
	565
ec6ab54b	566	//cout << "After rounding: " << value << " " << error << endl;
65e5de6a	567	return value; //taking into account ERROR
	568	}
	569	}
d03e5de4	570
ec6ab54b	571	TH1D* AliLRCAnalysis::GetAbsHisto() const
	572	{
	573	//Returns final histo in absolute variables
	574	return fPrAbs;
	575	}
	576	TH1D* AliLRCAnalysis::GetRelHisto() const
	577	{
	578	//Returns final histo in relative variables
	579	return fPrRel;
	580	}
	581
	582	TH1D* AliLRCAnalysis::GetForwardValueDist() const
	583	{
	584	//Returns destribution of value used in forward window
	585	return fPrf;
	586	}
	587	TH1D* AliLRCAnalysis::GetBackwardValueDist() const
	588	{
	589	//Returns destribution of value used in backward window
	590	return fPrb;
	591	}
	592
	593