/*
 * TpcClusterCalibDataset.cxx
 *
 *  Created on: Apr 24, 2015
 *      Author: mberger
 */

#include "TpcClusterCalibDataset.h"
#include "TH2.h"
#include "TH1.h"
#include <iostream>
#include <vector>
#include <math.h>

TpcClusterCalibDataset::TpcClusterCalibDataset(int nbinsZ, double startZ, double endZ)
:	fnZbins(nbinsZ),
 	fminZ(startZ),
 	fmaxZ(endZ),
 	fnCluster(0),
 	fuseTheta(false),
 	fnThetabins(0),
 	fminTheta(0),
 	fmaxTheta(0),
 	fnClusterTheta(0),
 	fbinWidthTheta(0)

{
	fDataA1 = new std::vector<double> [nbinsZ];
	fDataA2 = new std::vector<double> [nbinsZ];

	fbinWidth=(endZ-startZ)/nbinsZ;

	fDataThetaA1=NULL;
	fDataThetaA2=NULL;

	fDataZThetaA1=NULL;
	fDataZThetaA2=NULL;
}

TpcClusterCalibDataset::TpcClusterCalibDataset(int nbinsZ, double startZ, double endZ, int nbinsTheta, double startTheta, double endTheta)
:	fnZbins(nbinsZ),
 	fminZ(startZ),
 	fmaxZ(endZ),
 	fnCluster(0),
 	fuseTheta(true),
 	fnThetabins(nbinsTheta),
 	fminTheta(startTheta),
 	fmaxTheta(endTheta),
 	fnClusterTheta(0)

{
	fbinWidth=(endZ-startZ)/nbinsZ;
	fDataA1 = new std::vector<double> [nbinsZ];
	fDataA2 = new std::vector<double> [nbinsZ];

	fbinWidthTheta=(endTheta-startTheta)/nbinsTheta;
	fDataThetaA1 = new std::vector<double> [nbinsTheta];
	fDataThetaA2 = new std::vector<double> [nbinsTheta];

	fDataZThetaA1 = new std::vector<double> [nbinsZ*nbinsTheta];
	fDataZThetaA2 = new std::vector<double> [nbinsZ*nbinsTheta];
	/*
	for (int i=0;i++;i<nbinsZ)
	{
		fDataZThetaA1[i]=new std::vector<double> [nbinsTheta];
		fDataZThetaA2[i]=new std::vector<double> [nbinsTheta];
	}
	 */
}

TpcClusterCalibDataset::TpcClusterCalibDataset(const TpcClusterCalibDataset& other)
: fnZbins(other.fnZbins),
  fminZ(other.fminZ),
  fmaxZ(other.fmaxZ),
  fbinWidth(other.fbinWidth),
  fnCluster(other.fnCluster),
  fuseTheta(other.fuseTheta),
  fnThetabins(other.fnThetabins),
  fminTheta(other.fminTheta),
  fmaxTheta(other.fmaxTheta),
  fnClusterTheta(other.fnClusterTheta),
  fbinWidthTheta(other.fbinWidthTheta)
{
	//be aware that this copy constructor doesn't copy the data!
	fDataA1 = new std::vector<double> [other.fnZbins];
	fDataA2 = new std::vector<double> [other.fnZbins];

	fDataThetaA1 = new std::vector<double> [other.fnThetabins];
	fDataThetaA2 = new std::vector<double> [other.fnThetabins];

	fDataZThetaA1 = new std::vector<double> [other.fnZbins*other.fnThetabins];
	fDataZThetaA2 = new std::vector<double> [other.fnZbins*other.fnThetabins];
	/*
	for (int i=0;i++;i<other.fnZbins)
	{
		fDataZThetaA1[i]=new std::vector<double> [other.fnThetabins];
		fDataZThetaA2[i]=new std::vector<double> [other.fnThetabins];
	}
	 */
}

TpcClusterCalibDataset::~TpcClusterCalibDataset()
{

	delete[] fDataA1;
	delete[] fDataA2;

	if (fDataThetaA1!=NULL)
		delete[] fDataThetaA1;

	if (fDataThetaA2!=NULL)
		delete[] fDataThetaA2;

	if(fDataZThetaA1!=NULL)
	{
		delete[] fDataZThetaA1;
	}
	if(fDataZThetaA2!=NULL)
	{
		delete[] fDataZThetaA2;
	}
}

void TpcClusterCalibDataset::setTheta(int nbinsTheta, double startTheta, double endTheta)
{
	fuseTheta=true;
	fnThetabins=nbinsTheta;
	fminTheta=startTheta;
	fmaxTheta=endTheta;

	fbinWidthTheta=(fmaxTheta-fminTheta)/fnThetabins;

	fDataThetaA1 = new std::vector<double> [fnThetabins];
	fDataThetaA2 = new std::vector<double> [fnThetabins];
}

void TpcClusterCalibDataset::Reset(TString opt)
{
	for(int i=0;i<fnZbins;i++)
	{
		fDataA1[i].clear();
		fDataA2[i].clear();
	}

	if (fDataThetaA1!=NULL)
	{
		for(int i=0;i<fnThetabins;i++)
		{
			fDataThetaA1[i].clear();
			fDataThetaA2[i].clear();
		}
	}

	if (fDataThetaA1!=NULL)
	{
		for(int i=0;i<fnZbins;i++)
		{
			for(int j=0;j<fnThetabins;j++)
			{
				fDataZThetaA1[i*fnZbins+j].clear();
				fDataZThetaA2[i*fnZbins+j].clear();
			}
		}
	}
}

void TpcClusterCalibDataset::Add(const TpcClusterCalibDataset& other)
{
	fnCluster+=other.fnCluster;
	for(int i=0;i<other.fnZbins;i++)
	{
		for(int j=0;j<other.fDataA1[i].size();j++)
		{
			fDataA1[i].push_back(other.fDataA1[i][j]);
			fDataA2[i].push_back(other.fDataA2[i][j]);
		}
	}

	fnClusterTheta+=other.fnClusterTheta;
	if(fDataThetaA1!=NULL)
	{
		for(int i=0;i<other.fnThetabins;i++)
		{
			for(int j=0;j<other.fDataThetaA1[i].size();j++)
			{
				fDataThetaA1[i].push_back(other.fDataThetaA1[i][j]);
				fDataThetaA2[i].push_back(other.fDataThetaA2[i][j]);
			}
		}
	}
	if(fDataZThetaA1!=NULL)
	{
		for(int i=0;i<other.fnZbins;i++)
		{
			for(int j=0;j<other.fnThetabins;j++)
			{
				for (int k=0;k<other.fDataZThetaA1[i*fnZbins+j].size();k++)
				{
					fDataZThetaA1[i*fnZbins+j].push_back(other.fDataZThetaA1[i*fnZbins+j][k]);
					fDataZThetaA2[i*fnZbins+j].push_back(other.fDataZThetaA2[i*fnZbins+j][k]);
				}
			}
		}
	}

}

void TpcClusterCalibDataset::Fill(int iset, double zpos, double pull)
{
	int zbin = int((zpos-fminZ)/fbinWidth);

	if(zbin>fnZbins)
	{
		std::cout<<"TpcClusterCalibDataset::Fill: zbin "<<zbin<<" not available\n";
		return;
	}

	if (iset==1)
		fDataA1[zbin].push_back(pull);
	else if(iset==2)
		fDataA2[zbin].push_back(pull);
	else
		std::cout<<"TpcClusterCalibDataset::Fill Dataset "<<iset<<" is not defined\n";
}

void TpcClusterCalibDataset::FillTheta(int iset, double theta, double pull)
{
	int tbin= int((theta-fminTheta)/fbinWidthTheta);
	if (tbin>fnThetabins)
	{
		std::cout<<"TpcClusterCalibDataset::Fill: thetaBin "<<tbin<<" not available\n";
		return;
	}

	if (iset==1 and  fDataThetaA1!=NULL)
		fDataThetaA1[tbin].push_back(pull);
	else if(iset==2 and fDataThetaA2!=NULL)
		fDataThetaA2[tbin].push_back(pull);
	else
		std::cout<<"TpcClusterCalibDataset::Fill Dataset Theta "<<iset<<" is not defined\n";
}

void TpcClusterCalibDataset::FillZTheta(int iset, double zpos,double theta, double pull)
{
	int zbin = int((zpos-fminZ)/fbinWidth);
	int tbin= int((theta-fminTheta)/fbinWidthTheta);
	if(zbin>fnZbins)
	{
		std::cout<<"TpcClusterCalibDataset::FillZTheta: zbin "<<zbin<<" not available\n";
		return;
	}
	if (tbin>fnThetabins)
	{
		std::cout<<"TpcClusterCalibDataset::FillZTheta: thetaBin "<<tbin<<" not available\n";
		return;
	}
	if (iset==1)
		fDataZThetaA1[zbin*fnZbins+tbin].push_back(pull);
	else if(iset==2)
		fDataZThetaA2[zbin*fnZbins+tbin].push_back(pull);
	else
		std::cout<<"TpcClusterCalibDataset::Fill Dataset Z Theta "<<iset<<" is not defined\n";
}

TH1D TpcClusterCalibDataset::ProjectionY(int iset, TString name, int firstBin,int lastBin)
{
	TString tit="projection ";
	tit+=firstBin;
	tit+="_";
	tit+=lastBin;

	TH1D proj(name.Data(),tit.Data(),100,-3,3);

	//std::cout<<"projection from "<<firstBin<<" to "<<lastBin<<" from total "<<fnThetabins<<"\n";
	for(int i=firstBin;i<=lastBin;i++)
	{
		//std::cout<<"at bin "<<i<<"\n";

		if(iset==1 and  fDataThetaA1!=NULL)
		{
			//std::cout<<"A1: "<<fDataThetaA1[i].size()<<"\n";
			for(int j=0;j<fDataThetaA1[i].size();j++)
			{
				//std::cout<<" A1 at entry "<<j<<"\n";
				proj.Fill(fDataThetaA1[i][j]);
			}
		}
		else if(iset==2 and fDataThetaA2!=NULL)
		{
			//std::cout<<"A2: "<<fDataThetaA2[i].size()<<"\n";
			for(int j=0;j<fDataThetaA2[i].size();j++)
			{
				//std::cout<<" A2 at entry "<<j<<"\n";
				proj.Fill(fDataThetaA2[i][j]);
			}
		}
		else
			std::cout<<"TpcClusterCalibDataset::ProjectionY Dataset "<<iset<<" is not defined\n";
	}
	return proj;
}

std::vector<double>* TpcClusterCalibDataset::getZbinData(int iset, int ibin)
{
	if(iset==1)
		return (&fDataA1[ibin]);
	else if(iset==2)
		return (&fDataA2[ibin]);

	std::cout<<"TpcClusterCalibDataset::getZbinData Dataset "<<iset<<" is not defined\n";
	return NULL;
}

double TpcClusterCalibDataset::getRMS(int iset, int firstBin,int lastBin, double mean)
{
	double rms=0;
	if(mean == -99999)
		mean=getMean(iset,firstBin,lastBin);

	int normalize=0;
	for(int i=firstBin;i<=lastBin;i++)
		if(iset==1)
			for(int j=0;j<fDataA1[i].size();j++)
			{
				rms+=(fDataA1[i][j]-mean)*(fDataA1[i][j]-mean);
				normalize+=1;
			}
		else if(iset==2)
			for(int j=0;j<fDataA2[i].size();j++)
			{
				rms+=(fDataA2[i][j]-mean)*(fDataA2[i][j]-mean);
				normalize+=1;
			}
		else
			std::cout<<"TpcClusterCalibDataset::getRMS Dataset "<<iset<<" is not defined\n";
	if (normalize!=0)
		rms/=normalize;
	rms=sqrt(rms);
	return rms;
}

double TpcClusterCalibDataset::getRMS20N(int iset, int firstBin,int lastBin)
{
	double rms=0;
	for(int i=firstBin;i<=lastBin;i++)
		if(iset==1)
			for(int j=0;j<fDataA1[i].size();j++)
			{
				rms+=(fDataA1[i][j])*(fDataA1[i][j]);
			}
		else if(iset==2)
			for(int j=0;j<fDataA2[i].size();j++)
			{
				rms+=(fDataA2[i][j])*(fDataA2[i][j]);
			}
		else
			std::cout<<"TpcClusterCalibDataset::getRMS Dataset "<<iset<<" is not defined\n";
	return rms;
}

double TpcClusterCalibDataset::getMean(int iset, int firstBin,int lastBin)
{
	double mean=0;
	int normalize=0;
	for(int i=firstBin;i<=lastBin;i++)
		if(iset==1)
			for(int j=0;j<fDataA1[i].size();j++)
			{
				mean+=fDataA1[i][j];
				normalize+=1;
			}
		else if(iset==2)
			for(int j=0;j<fDataA2[i].size();j++)
			{
				mean+=fDataA2[i][j];
				normalize+=1;
			}
		else
			std::cout<<"TpcClusterCalibDataset::getRMS Dataset "<<iset<<" is not defined\n";
	if (normalize!=0)
		mean/=normalize;
	return mean;
}

double TpcClusterCalibDataset::getRMSError(int iset, int firstBin,int lastBin)
{
	// The right formula for RMS error depends on 4th momentum (see Kendall-Stuart Vol 1 page 243)
	// formula valid for only gaussian distribution ( 4-th momentum =  3 * sigma^4 )
	int nentries=0;
	for(int i=firstBin;i<=lastBin;i++)
		if(iset==1)
			for(int j=0;j<fDataA1[i].size();j++)
			{
				nentries++;
			}
		else if(iset==2)
			for(int j=0;j<fDataA2[i].size();j++)
			{
				nentries++;
			}
		else
			std::cout<<"TpcClusterCalibDataset::getRMSError Dataset "<<iset<<" is not defined\n";
	double rms=getRMS(iset,firstBin,lastBin);
	return ( nentries > 0 ? sqrt( (rms*rms) /(2*nentries) ) : 0. );
}

double TpcClusterCalibDataset::getMeanError(int iset, int firstBin,int lastBin)
{
	int nentries=0;
	for(int i=firstBin;i<=lastBin;i++)
		if(iset==1)
			nentries+=fDataA1[i].size();
	//for(int j=0;j<fDataA1[i].size();j++)
	//{
	//	nentries++;
	//	}
		else if(iset==2)
			nentries+=fDataA2[i].size();
	/*
			for(int j=0;j<fDataA2[i].size();j++)
			{
				nentries++;
			}
	 */
		else
			std::cout<<"TpcClusterCalibDataset::getMeanError Dataset "<<iset<<" is not defined\n";
	double rms=getRMS(iset,firstBin,lastBin);
	//std::cout<<"get mean err: "<<nentries<<" "<<rms<<"\n";
	return(nentries>0?rms/sqrt(nentries):0);
}

TH2D TpcClusterCalibDataset::getHist(int iset)
{
	TString name="Pulls Axis ";
	name += iset;
	TH2D hist(name,name,75,0,75,400,-20,20);

	for(int i=0;i<fnZbins;i++)
	{
		if(iset==1)
			for(int j=0;j<fDataA1[i].size();j++)
			{
				hist.Fill(i,fDataA1[i][j]);
			}
		else if(iset==2)
			for(int j=0;j<fDataA2[i].size();j++)
			{
				hist.Fill(i,fDataA2[i][j]);
			}
		else
			std::cout<<"TpcClusterCalibDataset::getHist Dataset "<<iset<<" is not defined\n";
	}
	return hist;
}

int TpcClusterCalibDataset::getEventsInBin(int firstBin, int lastBin)
{
	int n=0;
	for(int i=firstBin;i<=lastBin;i++)
	{
		n+=fDataA1[i].size();
	}
	return n;
}

void TpcClusterCalibDataset::getMeanRMS2N(double& mean1, double& mean2,double& rms2n1, double& rms2n2, double& rms02n1,double& rms02n2, double& resSum1, double& resSum2, int* nentries, int firstBin, int lastBin)
{
	mean1=0;
	mean2=0;
	rms2n1=0;
	rms2n2=0;
	rms02n1=0;
	rms02n2=0;
	resSum1=0;
	resSum2=0;
	nentries[0]=nentries[1]=0;

	for(int i=firstBin;i<=lastBin;i++)
	{
		for(int j=0;j<fDataA1[i].size();j++)
		{
			mean1+=fDataA1[i][j];
			mean2+=fDataA2[i][j];
			nentries[0]+=1;
			nentries[1]+=1;
		}
	}
	if(nentries[0]!=0)
		mean1/=nentries[0];
	else
		mean1=0;

	if (nentries[1]!=0)
		mean2/=nentries[1];
	else
		mean2=0;
	//get N*(rms)**2

	nentries[0]=nentries[1]=0;
	for(int i=firstBin;i<=lastBin;i++)
	{
		for(int j=0;j<fDataA1[i].size();j++)
		{
			if(fabs(fDataA1[i][j]-mean1)<3)
			{
				rms2n1+=pow((fDataA1[i][j]-mean1),2);
				rms02n1+=pow(fDataA1[i][j],2);
				resSum1+=(fDataA1[i][j]-mean1);
				nentries[0]+=1;
			}

			if(fabs(fDataA2[i][j]-mean2)<3)
			{
				rms2n2+=pow((fDataA2[i][j]-mean2),2);
				rms02n2+=pow(fDataA2[i][j],2);
				resSum2+=(fDataA2[i][j]-mean2);
				nentries[1]+=1;
			}
		}
	}

}

void TpcClusterCalibDataset::getMeanRMS2N_Theta(double& mean1, double& mean2,double& rms2n1, double& rms2n2, double& rms02n1,double& rms02n2, double& resSum1, double& resSum2, int* nentries, int firstBin, int lastBin)
{
	mean1=0;
	mean2=0;
	rms2n1=0;
	rms2n2=0;
	rms02n1=0;
	rms02n2=0;
	resSum1=0;
	resSum2=0;
	nentries[0]=nentries[1]=0;

	if(lastBin>=fnThetabins)
	{
		std::cout<<"requested bin out of range ("<<lastBin<<")\n";
		return;
	}
	if(firstBin<0)
	{
		std::cout<<"requested bin out of range ("<<firstBin<<")\n";
		return;
	}

	if(fDataThetaA1==NULL or fDataThetaA2==NULL)
	{
		std::cout<<"data not present\n";
		return;
	}


	//calculate means
	for(int i=firstBin;i<=lastBin;i++)
	{
		for(int j=0;j<fDataThetaA1[i].size();j++)
		{
			mean1+=fDataThetaA1[i][j];
			mean2+=fDataThetaA2[i][j];
			nentries[0]+=1;
			nentries[1]+=1;
		}
	}
	if(nentries[0]!=0)
		mean1/=nentries[0];
	else
		mean1=0;

	if (nentries[1]!=0)
		mean2/=nentries[1];
	else
		mean2=0;
	//get N*(rms)**2
	nentries[0]=nentries[1]=0;
	for(int i=firstBin;i<=lastBin;i++)
	{
		for(int j=0;j<fDataThetaA1[i].size();j++)
		{
			if(fabs(fDataThetaA1[i][j]-mean1)<3)
			{
				rms2n1+=pow((fDataThetaA1[i][j]-mean1),2);
				rms02n1+=pow(fDataThetaA1[i][j],2);
				resSum1+=(fDataThetaA1[i][j]-mean1);
				nentries[0]+=1;
			}

			if(fabs(fDataThetaA2[i][j]-mean2)<3)
			{
				rms2n2+=pow((fDataThetaA2[i][j]-mean2),2);
				rms02n2+=pow(fDataThetaA2[i][j],2);
				resSum2+=(fDataThetaA2[i][j]-mean2);
				nentries[1]+=1;
			}
		}
	}

}

void TpcClusterCalibDataset::getMeanRMS2N_ZTheta(double& mean1, double& mean2,double& rms2n1, double& rms2n2, double& rms02n1,double& rms02n2, double& resSum1, double& resSum2, int* nentries, int firstBinZ, int lastBinZ, int firstBinT, int lastBinT)
{
	mean1=0;
	mean2=0;
	rms2n1=0;
	rms2n2=0;
	rms02n1=0;
	rms02n2=0;
	resSum1=0;
	resSum2=0;
	nentries[0]=nentries[1]=0;

	if(lastBinT>=fnThetabins)
	{
		std::cout<<"requested Theta bin out of range ("<<lastBinT<<")\n";
		return;
	}
	if(firstBinT<0)
	{
		std::cout<<"requested Theta bin out of range ("<<firstBinT<<")\n";
		return;
	}
	if(lastBinZ>=fnZbins)
	{
		std::cout<<"requested Z bin out of range ("<<lastBinZ<<")\n";
		return;
	}
	if(firstBinZ<0)
	{
		std::cout<<"requested Z bin out of range ("<<firstBinZ<<")\n";
		return;
	}

	if(fDataZThetaA1==NULL or fDataZThetaA2==NULL)
	{
		std::cout<<"Z-Theta data not present\n";
		return;
	}

	//calculate means
	for(int i=firstBinZ;i<=lastBinZ;i++)
	{
		for (int j=firstBinT;j<=lastBinT;j++)
		{
			for(int k=0;k<fDataZThetaA1[i*fnZbins+j].size();k++)
			{
				mean1+=fDataZThetaA1[i*fnZbins+j][k];
				mean2+=fDataZThetaA2[i*fnZbins+j][k];
				nentries[0]+=1;
				nentries[1]+=1;
			}
		}
	}
	if(nentries[0]!=0)
		mean1/=nentries[0];
	else
		mean1=0;

	if (nentries[1]!=0)
		mean2/=nentries[1];
	else
		mean2=0;
	//get N*(rms)**2
	nentries[0]=nentries[1]=0;
	for(int i=firstBinZ;i<=lastBinZ;i++)
	{
		for(int j=firstBinT;j<=lastBinT;j++)
		{
			for(int k=0;k<fDataZThetaA1[i*fnZbins+j].size();k++)
			{
				if(fabs(fDataZThetaA1[i*fnZbins+j][k]-mean1)<3)
				{
					rms2n1+=pow((fDataZThetaA1[i*fnZbins+j][k]-mean1),2);
					rms02n1+=pow(fDataZThetaA1[i*fnZbins+j][k],2);
					resSum1+=(fDataZThetaA1[i*fnZbins+j][k]-mean1);
					nentries[0]+=1;
				}

				if(fabs(fDataZThetaA2[i*fnZbins+j][k]-mean2)<3)
				{
					rms2n2+=pow((fDataZThetaA2[i*fnZbins+j][k]-mean2),2);
					rms02n2+=pow(fDataZThetaA2[i*fnZbins+j][k],2);
					resSum2+=(fDataZThetaA2[i*fnZbins+j][k]-mean2);
					nentries[1]+=1;
				}
			}
		}
	}
}