#include <iostream>
#include <fstream>
using namespace std;

#include <math.h>
#include "psifamily.h"

Psifamily::Psifamily(Inputparameters& input,Beambeamsystem& bbsystem):Gammaanarrowvm(input,bbsystem)
{
//Defining width and mass...

	switch(input.getpidtest()){
		case StarlightConstants::JPSI:
                        width=0.000091;
                        mass=3.09692;
                break;
                case StarlightConstants::JPSI2S:
                        width=0.000337;
                        mass=3.686093;
                break;
		default: cout<<"This PSI Family Member Has Not Been Defined, Psifamily::Psifamily()"<<endl;
	}

}

Psifamily::~Psifamily()
{
}

double Psifamily::gettheta(StarlightConstants::particle ipid)
{
//should probably merge the psi fmaily back to the vm stuff.

//This depends on the decay angular distribution
//Valid for J/Psi, Psi(2s)?
//cout<<"Psi family theta"<<endl;
double theta=0.;
double xtest=0.;
double dndtheta=0.;
        L200td:
          theta = StarlightConstants::pi*Randy.Rndom();//random()/(RAND_MAX+1.0);
          xtest = Randy.Rndom();//random()/(RAND_MAX+1.0);
          //  Follow distribution for helicity +/-1
          //  Eq. 19 of J. Breitweg et al., Eur. Phys. J. C2, 247 (1998)//Does Not Apply for J/psi?
          //  SRK 11/14/2000

          dndtheta = sin(theta)*(1.+((cos(theta))*(cos(theta))));
          if(xtest > dndtheta)
            goto L200td;
        return theta;
}

double Psifamily::getdaughtermass(StarlightConstants::particle &ipid)
{
	double ytest=0.,mdec=0.;
	//  decays 50% to e+/e-, 50% to mu+/mu-
        ytest = Randy.Rndom();//random()/(RAND_MAX+1.0);
        if(ytest >= 0.5)
        {
	        mdec = StarlightConstants::mel;
        	ipid = StarlightConstants::ELECTRON;
        }
        else
        {
	        mdec = StarlightConstants::mmu;
        	ipid = StarlightConstants::MUON;
        }
	return mdec;


}
