#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cstdio>
#include "particleType.h"
#include "particle.h"
#include "resonanceType.h"
using namespace std;

particleType *particle::fParticleType[particle::fMaxNumParticleType];

int particle::fNparticleType = 0;

int particle::FindParticle(const char *name){
  for(int i=0;i<fNparticleType;i++){
    const char *currentType = fParticleType[i]->getName();
    int k=0;
    while((currentType[k] == name[k]) && currentType[k] != '\0' && name[k] != '\0'){
      k++;
    } 
    if(currentType[k] == name[k]) return i;    
  }
  return -1; // if there's no match
}
//once a function returns a value, any other return in the same function is automatically thrown into hell
//that's why this FindParticle doesn't always return -1 as it would look like

particle::particle(){
 fIparticle = 0;
 fPx = 0.;
 fPy = 0.;
 fPz = 0.;
}

particle::particle(const particle& p){
 fIparticle = p.fIparticle;
 fPx = p.fPx;
 fPy = p.fPy;
 fPz = p.fPz;
}

particle::particle(int iParticle,double px,double py,double pz): fPx(px), fPy(py), fPz(pz){
 if (iParticle < fNparticleType && iParticle>=0) fIparticle = iParticle;
 else { cout << "particle type doesn't exist in the array" << endl; 
  fIparticle = -1;}
}

particle::particle(const char* name, double px,double py, double pz): fPx(px), fPy(py), fPz(pz){
 int flag=FindParticle(name);
 if (flag != -1) fIparticle = flag;
 else {
  cout << "particle type doesn't exist in the array" << endl;
  fIparticle = -1;}
}

int particle::addParticleType(const char* name, double mass, int charge, double width){
 if (fNparticleType < fMaxNumParticleType){
  int flag = FindParticle(name);
  if (flag != -1) {
   cout << "a particle of this type already exists in the array (nothing added)" << endl;
   return 2; }
  if (flag == -1 && width > 0){
   fParticleType[fNparticleType] = new resonanceType(name, mass, charge, width);
   fNparticleType++;}
  else if (flag == -1 && width < 0){
   fParticleType[fNparticleType] = new particleType(name, mass, charge);
   fNparticleType++;}}
 else {
  cout << "array is full (nothing added)" << endl;
  return 1;}
}

void particle::printParticleType() {
 if (fNparticleType){
  for (int i=0; i<fNparticleType; i++){
   if(!(fParticleType[i]->isRes())) fParticleType[i]->print();
   else ((resonanceType *) fParticleType[i]) -> print(); }
  cout << endl; }
 else cout << "no particle type added yet" << endl; 
}

void particle::print() const{
 if(fIparticle != -1){
  cout << "id : " << fIparticle << endl;
  cout << "name : " << fParticleType[fIparticle] -> getName() << endl;
  cout << "px, py, pz: " << fPx << ", " << fPy << ", " << fPz << endl; }
 else {
  cout << "particle type is not valid" << endl;
  cout << "px, py, pz: " << fPx << ", " << fPy << ", " << fPz << endl; }
}

bool particle::operator==(const particle &value){
 bool flag;
 if (fPx == value.fPx && fPy == value.fPy && fPy && fPz == value.fPz && fIparticle == value.fIparticle) flag = true;
 else flag = false;
 return flag; 
}

particle& particle::operator=(const particle &value){
  if (&value==this)return *this; 
  else {
  fPx = value.fPx;
  fPy = value.fPy;
  fPz = value.fPz;
  fIparticle = value.fIparticle;
  return *this; }
}
 
particle& particle::operator+=(const particle & value){
  fPx += value.fPx;
  fPy += value.fPy;
  fPz += value.fPz;
  return *this;
}

particle  particle::operator+(const particle &value) const{
  particle result(*this);  
  result += value;  
  return result;
}

double particle::GetMass() const {
 if (fIparticle > -1) {
  return fParticleType[fIparticle] -> getMass(); }
 else return 0.0;
}

double particle::GetEnergy() const {
 double mass = GetMass();
 return sqrt(mass*mass + fPx*fPx + fPy*fPy + fPz*fPz);
}

double particle::invMass(particle & p) const{
 double energy = GetEnergy() + p.GetEnergy();
 double p2 = pow(fPx+p.fPx, 2) + pow(fPy+p.fPy, 2) + pow(fPz+p.fPz, 2);
 return sqrt(energy*energy - p2); 
}

void particle::changeParticleType(int iParticle){
  if(iParticle >=0 && iParticle < fNparticleType) fIparticle = iParticle;
  else printf("Particle %d doesn't exist\n",iParticle);
}
 
void particle::changeParticleType(const char *name){
  int ip=FindParticle(name);
  if(ip != -1){
    fIparticle = ip;
  }
  else{
    printf("Particle \"%s\" doesn't exist in the stack\n",name);
  }
}

//implemented extra methods

int particle::Decay2body(particle &dau1,particle &dau2) const {
  if(GetMass() == 0.0){
    printf("Decayment cannot be preformed if mass is zero\n");
    return 1;
  }
  
  double massMot = GetMass();
  double massDau1 = dau1.GetMass();
  double massDau2 = dau2.GetMass();

  if(fIparticle > -1){ // add width effect

    // gaussian random numbers

    float x1, x2, w, y1, y2;
    
    double invnum = 1./RAND_MAX;
    do {
      x1 = 2.0 * rand()*invnum - 1.0;
      x2 = 2.0 * rand()*invnum - 1.0;
      w = x1 * x1 + x2 * x2;
    } while ( w >= 1.0 );
    
    w = sqrt( (-2.0 * log( w ) ) / w );
    y1 = x1 * w;
    y2 = x2 * w;

    massMot += fParticleType[fIparticle]->getWidth() * y1;

  }

  if(massMot < massDau1 + massDau2){
    printf("Decayment cannot be preformed because mass is too low in this channel\n");
    return 2;
  }
  
  double pout = sqrt((massMot*massMot - (massDau1+massDau2)*(massDau1+massDau2))*(massMot*massMot - (massDau1-massDau2)*(massDau1-massDau2)))/massMot*0.5;

  double norm = 2*M_PI/RAND_MAX;

  double phi = rand()*norm;
  double theta = rand()*norm*0.5 - M_PI/2.;
  dau1.setP(pout*sin(theta)*cos(phi),pout*sin(theta)*sin(phi),pout*cos(theta));
  dau2.setP(-pout*sin(theta)*cos(phi),-pout*sin(theta)*sin(phi),-pout*cos(theta));

  double energy = sqrt(fPx*fPx + fPy*fPy + fPz*fPz + massMot*massMot);

  double bx = fPx/energy;
  double by = fPy/energy;
  double bz = fPz/energy;

  dau1.Boost(bx,by,bz);
  dau2.Boost(bx,by,bz);

  return 0;
}
void particle::Boost(double bx, double by, double bz)
{

  double energy = GetEnergy();

  //Boost this Lorentz vector
  double b2 = bx*bx + by*by + bz*bz;
  double gamma = 1.0 / sqrt(1.0 - b2);
  double bp = bx*fPx + by*fPy + bz*fPz;
  double gamma2 = b2 > 0 ? (gamma - 1.0)/b2 : 0.0;

  fPx += gamma2*bp*bx + gamma*bx*energy;
  fPy += gamma2*bp*by + gamma*by*energy;
  fPz += gamma2*bp*bz + gamma*bz*energy;
}