//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
//      Implementation of class TpcCluster
//      see TpcCluster.h for details
//
// Environment:
//      Software developed for the PANDA Detector at FAIR.
//
// Author List:
//      Sebastian Neubert    TUM            (original author)
//
//
//-----------------------------------------------------------

// Panda Headers ----------------------

// This Class' Header ------------------
#include "TpcCluster.h"


// C/C++ Headers ----------------------
#include <iostream>
#include <set>
using namespace std;

// Collaborating Class Headers --------
#include "TVectorD.h"
#include "TMatrixDEigen.h"
#include "TMath.h"

// Class Member definitions -----------


//TODO: get rid of these horrible const_casts!!!!!! 

TpcCluster::TpcCluster()
  : fpos(0,0,0),
    fcogpos(0,0,0),
    fsig(0.1,0.1,0.01),fsigT(0,0,0), fsigA(0,0,0), fFieldCorr(0,0,0),fcov(3,3), fcovOnPlane(2,2),faxis(0,0,0),fhasaxis(false),
    famp(0), fnPad(0), fnPadX(0), fnPadY(0), findex(0), fsector(0), fsaturated(false),
    findexInTrack(-1),fphi(0),fhasSigA(false),fhasSigT(false),fhasSigC(false),feigenVect(3,3),fhasCov(false),fsigWrong(-1,-1,-1),
    ffitGood(false),fisGeomErr(0,0,0),feigenValOnPlane(0,0,0),fminSigma(0,0,0),fminNll(1000),fshapeCov(3,3),funwShapeCov(3,3),
    fshapeCovPlane(2,2),funwShapeCovPlane(2,2)
{}

TpcCluster::TpcCluster(const TpcCluster& clust) 
  : fpos(clust.fpos), fcogpos(clust.fcogpos),fsig(clust.fsig), fsigT(clust.fsigT),fcov(clust.fcov), fsigA(clust.fsigA),
    fsigC(clust.fsigC),faxis(clust.faxis),fFieldCorr(clust.fFieldCorr),
    famp(clust.famp), fhasaxis(clust.fhasaxis), fnPad(clust.fnPad), 
    fnPadX(clust.fnPadX),fnPadY(clust.fnPadY),
    findex(clust.findex),fsector(clust.fsector),findexInTrack(clust.findexInTrack),
    fmcid(clust.fmcid), digis(clust.digis), padids(clust.padids),fsaturated(clust.fsaturated),
    fphi(clust.fphi),fhasSigC(clust.fhasSigC),fhasSigA(clust.fhasSigA),fhasSigT(clust.fhasSigT),
    feigenVect(clust.feigenVect),fhasCov(clust.fhasCov),fsigWrong(clust.fsigWrong),
    fcovOnPlane(clust.fcovOnPlane),fmidPlane(clust.fmidPlane),ffitGood(clust.ffitGood),fisGeomErr(clust.fisGeomErr),
    feigenValOnPlane(clust.feigenValOnPlane),fminSigma(clust.fminSigma),fminNll(clust.fminNll),fshapeCov(clust.fshapeCov),
    funwShapeCov(clust.funwShapeCov),fshapeCovPlane(clust.fshapeCovPlane),funwShapeCovPlane(clust.funwShapeCovPlane),
    fprfWeights(clust.fprfWeights)
{}

TpcCluster::TpcCluster(const TVector3& Pos, double Amp, 
			     unsigned int Index)
  : fpos(Pos), fcogpos(0,0,0),fsig(0.1,0.1,0.01),fsigT(0,0,0), fsigA(0,0,0), fFieldCorr(0,0,0),fcov(3,3), fcovOnPlane(2,2),faxis(0,0,0),fhasaxis(false),famp(Amp),
    fnPad(0), fnPadX(0), fnPadY(0),findex(Index),fsector(0),fsaturated(false),findexInTrack(-1),fphi(0),
    fhasSigC(false),fhasSigA(false),fhasSigT(false),feigenVect(3,3),fhasCov(false),fsigWrong(-1,-1,-1),
    ffitGood(false),fisGeomErr(0,0,0),feigenValOnPlane(0,0,0),fminSigma(0,0,0),fminNll(1000),fshapeCov(3,3),funwShapeCov(3,3),
    fshapeCovPlane(2,2),funwShapeCovPlane(2,2)
{}

TpcCluster::TpcCluster(const TVector3& Pos, const TVector3& Sig, double Amp, 
			     unsigned int Index)
  : fpos(Pos), fcogpos(0,0,0),fsig(Sig) ,fsigT(0,0,0), fsigA(0,0,0),fFieldCorr(0,0,0), fcov(3,3), fcovOnPlane(2,2),faxis(0,0,0),fhasaxis(false),famp(Amp),
    fnPad(0), fnPadX(0), fnPadY(0),findex(Index),fsector(0),fsaturated(false),findexInTrack(-1),fphi(0),fhasSigC(false),
    fhasSigA(false),fhasSigT(false),feigenVect(3,3),fhasCov(false),fsigWrong(-1,-1,-1),ffitGood(false),fisGeomErr(0,0,0),
    feigenValOnPlane(0,0,0),fminSigma(0,0,0),fminNll(1000),fshapeCov(3,3),funwShapeCov(3,3),fshapeCovPlane(2,2),funwShapeCovPlane(2,2)
{}

void
TpcCluster::addDigi(const TpcDigi* d,double weight){
  if(!digis.count(d->index()))
  {
    digis[d->index()]=weight;
    padids.insert(d->padId());
  }else{
    std::cerr<<"WARNING: TpcCluster::addDigi(unsigned int digiIndex="
             <<d->index()<<",double weight="<<weight
             <<") called where sampleIndex is already existing, overwriting"<<std::endl;
    digis[d->index()]=weight;
    padids.insert(d->padId());
  }
  
}

TVector3
TpcCluster::calcAxis(unsigned int i,Bool_t norm) const {
  //assert(i<3);//proper error message here 
  //TVectorD eigenValues(3);
  //  std::cout<<"here:cluster";
  //fcov.Print();
  //TMatrixD eigenVec=fcov.EigenVectors(eigenValues);
  TMatrixD eigenVec(3,3);
  if (fhasaxis)
    eigenVec=feigenVect;
  else
    {
      TMatrixDEigen eigen(fcov);
      eigenVec=eigen.GetEigenVectors();
    }
  //  (eigen.GetEigenValues()).Print();
  //eigen.GetEigenValuesRe().Print();
  //  eigen.GetEigenValuesIm().Print();
  //  std::cout<<"here2:cluster\n";
  // eigenvalues are sorted according to their value
  // in descendign order -> first one is largest
  //  std::cout<<"cluster eigenvect\n";
  // eigenVec.Print();
  TVectorD a=TMatrixDColumn(eigenVec,i);
  TVector3 Axis(a[0],a[1],a[2]);
  if (norm)
    Axis.SetMag(1);
  return Axis;
}

unsigned int 
TpcCluster::get2DSize() { //return the 2-dimensional size of the cluster
  return padids.size();
}

double
TpcCluster::getUnsharedSize()
{
  double usize=0;
  std::map<unsigned int, double>::const_iterator it;
  for(it=digis.begin();it!=digis.end();it++)
    {
      usize+=it->second;
      //std::cout<<"clusterdingens: "<<it->second<<"\n";
    //    if (it->second >= 0.5)
    //usize++;
    }

  return usize;
}

TpcCluster::~TpcCluster(){}

std::ostream& operator<< (std::ostream& s, const TpcCluster& me){
  s << "TpcCluster\n"
    << "pos=("<<me.fpos.X()<<","<<me.fpos.Y()<<","<<me.fpos.Z()<<")\n"
    << "amp=" << me.famp << "\n"
    << "size=" << me.digis.size() <<"\n"
    << "index="<< me.findex;
  return s;
} 

std::vector<unsigned int>
TpcCluster::getDigiIDs() const {
  std::vector<unsigned int> ret;
  std::map<unsigned int, double>::const_iterator it;
  for(it=digis.begin(); it!=digis.end(); it++)
    ret.push_back(it->first);
  return ret;
}

const TVector3& TpcCluster::sig(int mode) const
{
  switch(mode)
    {
    case 0:
      return  fsig;
    case 1:
      if(fhasSigT)
    	  return fsigT;
      else
    	  return fsigWrong;
    case 2:
      if (fhasSigC)
    	  return fsigC;
      else
    	  return fsigWrong;
    case 3:
      if (fhasSigA)
    	  return fsigA;
      else
    	  return fsigWrong;
    default:
      std::cout<<"TpcCluster::be careful! returning default sig!\n";
      return fsig; 

    }

}

void TpcCluster::calcSigCA()
{
  TMatrixDSym cov_i=TMatrixDSym(3,fcov.GetMatrixArray());
  cov_i.InvertFast();
    //caclulate intersection of ellipse with x y and z axis:
  for(int i=0;i<3;++i)
    {
      TVectorD C(3);
      for(int j=0;j<3;++j)
	if (i==j)
	  C[j]=1;
	else
	  C[j]=0;
      
      Double_t lambda=cov_i.Similarity(C);
      
      fsigC[i]=1./TMath::Sqrt(lambda);
    }
  fhasSigC=true;
}

bool TpcCluster::hasGeomErr() const
{
	for (int i=0;i<3;i++)
	{
		if(fisGeomErr[i])
			return true;
	}
	return false;
}


ClassImp(TpcCluster)