#include "mypart.h"
#include <TLorentzVector.h>
#include <iostream>
#include <math.h>

using namespace std ;

mypart::~mypart() {}
mypart::mypart():
  _ch(0), _n(0) , _zv(0), _ez(0), _x2(0)
 { _fv.SetXYZM(0,0,0,0) ;}

mypart::mypart( int charge, const TLorentzVector FourVector, double ZV, double eZ ):
  _ch(charge), _n(1) , _zv(ZV), _ez(eZ), _x2(0)
 { _fv = FourVector ;}

mypart::mypart( int charge, const TLorentzVector FourVector, double ZV, double eZ, double x2 ):
  _ch(charge), _n(1), _zv(ZV), _ez(eZ), _x2(x2)
 { _fv = FourVector ;}

mypart::mypart( const mypart & P ) { copy( P ) ;}

mypart::mypart( int charge, double px, double py, double pz, double mass, double ZV, double eZ ):
  _ch(charge), _n(1), _zv(ZV), _ez(eZ), _x2(0)
{ _fv.SetXYZM( px, py, pz, mass )  ;}


mypart& mypart::operator=( const mypart& P ) 
{
  copy( P ) ;
  return *this;
}

mypart& mypart::operator+( const mypart& P ) 
{  
  mypart *result = new mypart ;
  *result = *this ;
  *result += P ;
  return *result;
}
mypart& mypart::operator +=( const mypart& P ) 
{
  _n += P._n ;
  _ch += P._ch ;
  _fv += P._fv ;
  ComputeCommonVertex( P ) ;
  return *this;
}

void mypart::ComputeCommonVertex( const mypart& P ) 
{
  if( P._ez == 0 ) return ;
 
  if( _ez == 0 ) 
    {
      _ez = P._ez ;
      _zv = P._zv ;
      _x2 = P._x2 ;
      return ;
    }

  double delta = _zv-P._zv ;
  double chi2 = delta*delta/(_ez*_ez + P._ez*P._ez ) ;
  double sum = _zv/_ez/_ez + P._zv/P._ez/P._ez ;
  double wgt = 1./_ez/_ez + 1./P._ez/P._ez ;

  _zv = sum/wgt ;
  _ez = sqrt(1./wgt) ;
  _x2 = chi2+_x2+P._x2 ;
} 
TLorentzVector mypart::GetFourVector() { return _fv ; }
int    mypart::GetCharge(){ return _ch ; }
int    mypart::GetMultiplicity() { return _n ; }
double mypart::GetVertexPosition(){ return _zv ; }
double mypart::GetVertexError(){ return _ez ; }
double mypart::GetChiSquare() { return _x2 ; }
void mypart::copy( const mypart& P ) 
{
  _n = P._n ;
  _ch = P._ch ;
  _zv = P._zv ;
  _ez = P._ez ;
  _fv = P._fv ;
  _x2 = P._x2 ;
  return ;
}