#include <iostream>
#include <fstream>
using namespace std;
                                                                                
#include <math.h>
//#include <cmath>
//#include <vector>
//#include "starlightconstants.h"
#include "gammaavm.h"
#include "gammaacrosssection.h"

Gammaavectormeson::Gammaavectormeson(Inputparameters& input,Beambeamsystem& bbsystem):Eventchannel(input,bbsystem)//:Readinluminosity(input),bbs(bbsystem)
{
	VMNPT=input.getNPT();
	VMWmax=input.getWmax();
	VMWmin=input.getWmin();
	VMYmax=input.getYmax();
	VMYmin=-1.*VMYmax;
	VMnumw=input.getnumw();
	VMnumy=input.getnumy();
	VMgamma_em=input.getgamma_em();
	VMinterferencemode=input.getinterferencemode();
	VMbslope=0.;//Will define in wide/narrow constructor
	VMpidtest=input.getpidtest();
	VMptmax=input.getmaximuminterpt();
	VMdpt=input.getdpt();
	Randy.SetSeed(input.getseed());
	switch(VMpidtest){
                case StarlightConstants::RHO:
                case StarlightConstants::RHOZEUS:
			width=0.1507;
			mass=0.7685;
		break;
                case StarlightConstants::OMEGA:
			width=0.00843;
			mass=0.78194;
                break;
                case StarlightConstants::PHI:
			width=0.00443;
			mass=1.019413;
                break;
               /* case StarlightConstants::JPSI:
			width=0.000091;
			mass=3.09692;
		break;
                case StarlightConstants::JPSI2S:
			width=0.000337;
			mass=3.686093;
		break;
		*/
                case StarlightConstants::UPSILON:
			width=0.00005402;
			mass=9.46030;
		break;
                case StarlightConstants::UPSILON2S:
			width=0.00003198;
			mass=10.02326;
		break;
                case StarlightConstants::UPSILON3S:
			width=0.00002032;
			mass=10.3552;
                break;
                default: cout<<"No proper vector meson defined, gammaavectormeson::gammaavectormeson"<<endl;
         }


}

Gammaavectormeson::~Gammaavectormeson()
{

}

void Gammaavectormeson::pickwy(double &W, double &Y)
{

        //cout <<"momentapickwya"<<endl;
  double dW, dY, xw,xy,xtest;
  int  IW,IY;

  dW = (VMWmax-VMWmin)/double(VMnumw);
  dY = (VMYmax-VMYmin)/double(VMnumy);
  //       >> DRAW xw,xy
 L201pwy:

  xw = Randy.Rndom();// random()/(RAND_MAX+1.0);

  W = VMWmin + xw*(VMWmax-VMWmin);
//      cout <<"W= "<<W<<endl;
  //  C  protect against sub-threshold events

  if (W < 2*StarlightConstants::mpi)
    goto L201pwy;

  IW = int((W-VMWmin)/dW); //+ 1;
  xy = Randy.Rndom();//random()/(RAND_MAX+1.0);
  Y = VMYmin + xy*(VMYmax-VMYmin);
  IY = int((Y-VMYmin)/dY); //+ 1;
//  cout <<"Y IW IY= "<<Y<<" "<<IW<<" "<<IY<<endl;
  //       >> Check if this matches the map
    //if(Y < 0.) goto L201pwy;  commented out in sergey's

  xtest = Randy.Rndom();//random()/(RAND_MAX+1.0);

  if( xtest > Farray[IW][IY] )
    goto L201pwy;

 //       >> W,Y should now have the right distribution
//      cout <<"pick xw,xy,xtest,w,y= "<<xw<<" "<<" "<<xy<<" "<<xtest<<" "<<W<<" "<<Y<<endl;

//      Egam = 0.5*W*exp(Y);
//      Epom = 0.5*W*exp(-Y);

//      if(Z2==1&&Epom>30.) goto L201pwy;

}         


                                                                                                       
void Gammaavectormeson::twodecay(StarlightConstants::particle &ipid,double E,double W,double px0,double py0,double pz0,double &px1,double &py1,double&pz1,double &px2,double &py2,double &pz2,int &iFbadevent)
{
                                                                                                                
  //     This routine decays a particle into two particles of mass mdec,
  //     taking spin into account

  //  subroutine twodecay (ipid,E,W,px0,py0,pz0,mdec,px1,py1,pz1,E1,px2,py2,pz2,E2,iFbadevent)
  double pmag, anglelep[20001],ytest=0.;
  double phi,theta,xtest,dndtheta,Ecm;
  double  betax,betay,betaz;
  double mdec=0.0;
  double E1=0.0,E2=0.0;
  //    set the mass of the daughter particles
	mdec=getdaughtermass(ipid);
        //     calculate the magnitude of the momenta
              if(W < 2*mdec)
                {
                  cout<<" ERROR: W="<<W<<endl;
                  iFbadevent = 1;
                  return;
                }
              pmag = sqrt(W*W/4. - mdec*mdec);
                                                                                                                
        //     pick an orientation, based on the spin
        //      phi has a flat distribution in 2*pi
        phi = Randy.Rndom()*2.*StarlightConstants::pi;//(random()/(RAND_MAX+1.0))* 2.*pi;
                                                                                                                
        //     find theta, the angle between one of the outgoing particles and
        //    the beamline, in the frame of the two photons

//Make a get theta function with dowhile...this would allow it to be overloaded?

      //if(spin == 1.) theta=gettheta();
      //if(spin != 1) cout<<" This model cannot yet handle this spin value: "<<spin<<endl;
	//We know it is spin==1...
//	cout<<"Grabbing theta"<<endl;
	theta=gettheta(ipid);
//	cout<<"Theta grabbed"<<endl;
      
      //     compute unboosted momenta
      px1 = sin(theta)*cos(phi)*pmag;
      py1 = sin(theta)*sin(phi)*pmag;
      pz1 = cos(theta)*pmag;
      px2 = -px1;
      py2 = -py1;
      pz2 = -pz1;
//      cout <<"2decay1 px1 py1 pz1= "<<px1<<" "<<py1<<" "<<pz1<<endl;
      //        compute energies
      //      if(ip.eq.33) then
                                                                                                                
      //Changed mass to W 11/9/2000 SRK
//      cout<<"spin1, theta,xtest, px1,py1,pz1: "<<theta<<" "<<xtest<<" "<<px1<<" "<<py1<<" "<<pz1<<endl;
      Ecm = sqrt(W*W+px0*px0+py0*py0+pz0*pz0);
                                                                                                                
          //      else
          //      Ecm = sqrt(E**2 + px0**2 + py0**2 + pz0**2)
          //      endif
      E1 = sqrt(mdec*mdec+px1*px1+py1*py1+pz1*pz1);
      E2 = sqrt(mdec*mdec+px2*px2+py2*py2+pz2*pz2);
//      cout<<"Ecm E1 E2= "<<Ecm<<" "<<E1<<" "<<E2<<endl;
                                                                                                                
      //     Lorentz transform into the lab frame
      // betax,betay,betaz are the boost of the complete system
      betax = -(px0/Ecm);
      betay = -(py0/Ecm);
      betaz = -(pz0/Ecm);
//      cout <<"2decaybetax Ecm px1 px0= "<<betax<<" "<<Ecm<<" "<<px1<<" "<<px0<<endl;
      transform (betax,betay,betaz,E1,px1,py1,pz1,iFbadevent);
      transform (betax,betay,betaz,E2,px2,py2,pz2,iFbadevent);
//      cout <<"2decay2 px1= "<<px1<<endl;
//      cout<<"iFbadevent: "<<iFbadevent<<endl;
      if(iFbadevent == 1)
        return;

//      cout<<"twodecayend px1 py1 pz1 E1 px2 py2 pz2 E2= "<<px1<<" "<<py1<<" "<<pz1<<" "<<E1<<" "<<px2<<" "<<py2<<" "<<pz2<<" "<<E2<<endl;

}
double Gammaavectormeson::getdaughtermass(StarlightConstants::particle &ipid)
{

	//This will return the daughter particles mass, and the final particles outputed id...
	double ytest=0.,mdec=0.;

	switch(VMpidtest){
		case StarlightConstants::RHO:
		case StarlightConstants::RHOZEUS:
		case StarlightConstants::OMEGA:
			 mdec = StarlightConstants::mpi;
			 ipid = StarlightConstants::PION;
		break;
		case StarlightConstants::PHI:
			 mdec = StarlightConstants::mK;
			 ipid = StarlightConstants::KAONCHARGE;
		break;
		//case StarlightConstants::JPSI:
		//case StarlightConstants::JPSI2S:
		case StarlightConstants::UPSILON:
		case StarlightConstants::UPSILON2S:
		case StarlightConstants::UPSILON3S:
	            //  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;
              	//     }
		break;
		default: cout<<"No daughtermass defined, gammaavectormeson::getdaughtermass"<<endl;
         }

	return mdec;
}
double Gammaavectormeson::gettheta(StarlightConstants::particle ipid)
{
//This depends on the decay angular distribution
//Valid for rho, phi, omega.
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)
          //  SRK 11/14/2000
//        cout<<"spin1"<<endl;
	switch(ipid){

	case StarlightConstants::MUON:
	case StarlightConstants::ELECTRON:
	//primarily for upsilon/j/psi.  VM->ee/mumu
		dndtheta = sin(theta)*(1.+((cos(theta))*(cos(theta))));
	break;

	case StarlightConstants::PION:
	case StarlightConstants::KAONCHARGE:
	//rhos etc
          dndtheta= sin(theta)*(1.-((cos(theta))*(cos(theta))));
	break;
	
	default: cout<<"No proper theta dependence defined, check gammaavectormeson::gettheta"<<endl;
	}//end of switch

          if(xtest > dndtheta)
            goto L200td;

	return theta;

}

/*void Gammaavectormeson::transform(double betax,double betay,double betaz,double &E,double &px,double &py,double &pz,int &iFbadevent)
{

  //     carries out a lorentz transform of the frame.  (Not a
  //     boost!)
  double beta,gamma,gob;
  double E0,px0,py0,pz0;

  E0 = E;
  px0 = px;
  py0 = py;
  pz0 = pz;
//cout <<"transform betax, y, z="<<betax<<" "<<betay<<" "<<betaz<<endl;
  beta = sqrt(betax*betax + betay*betay + betaz*betaz);
//cout <<"transform beta= "<<beta<<endl;
  if(beta >= 1.0)  iFbadevent = 1;
  gamma = 1./sqrt(1. - beta*beta);
  gob = (gamma-1)/(beta*beta);
  
  E = gamma*(E0 - betax*px0 - betay*py0 - betaz*pz0);

  px = -gamma*betax*E0 + (1. + gob*betax*betax)*px0  +  gob*betax*betay*py0 + gob*betax*betaz*pz0;
  
  py = -gamma*betay*E0 + gob*betay*betax*px0+  (1. + gob*betay*betay)*py0 + gob*betay*betaz*pz0;
  
  pz = -gamma*betaz*E0 +  gob*betaz*betax*px0 +  gob*betaz*betay*py0 + (1. + gob*betaz*betaz)*pz0;
 
}
*/
double Gammaavectormeson::getwidth()
{
return width;
}
double Gammaavectormeson::getmass()
{
return mass;
}

double Gammaavectormeson::getspin()
{
return 1.0; //VM spins are the same
}

void Gammaavectormeson::momenta(double W,double Y,double &E,double &px,double &py,double &pz,int &tcheck)//may want to remove double w, y//better way for bslope...maybe store it as a vm parameter instead of cs parameter...not sure.
{
// This subroutine calculates momentum and energy of vector meson
  //     given W and Y,   without interference.  Subroutine vmpt.f handles
  //     production with interference

  //    SUBROUTINE VMOMENTA(W,Y,E,px,py,pz,tcheck)
//cout <<"vmomentaa W= "<<W<<" Y= "<<Y<<endl;

  double dW,dY;

  double Egam,Epom,tmin,pt1,pt2,phi1,phi2;
  double px1,py1,px2,py2;
  double pt,xt,xtest;
  double photon_spectrum;
  double t1,t2;

  dW = (VMWmax-VMWmin)/double(VMnumw);
  dY  = (VMYmax-VMYmin)/double(VMnumy);

  //       >> Find Egam,Epom in CM frame
  Egam = 0.5*W*exp(Y);
  Epom = 0.5*W*exp(-Y);
  
  //       >> Draw pt1 (for the photon) according to photon pt spectrum (Ref ?)
 L202vm:
  xt = Randy.Rndom();//random()/(RAND_MAX+1.0);
  pt1  = 0.5*xt;
  //       >> Check tmin
  tmin = ((Egam/VMgamma_em)*(Egam/VMgamma_em));
  if(tmin > 0.5)
    {
      cout<< " WARNING: tmin= "<<tmin<<endl;
      cout<< " Will pick a new W,Y "<<endl;
      tcheck = 1;
      return;
    }
 
  xtest = Randy.Rndom();//random()/(RAND_MAX+1.0);
  t1 = tmin + pt1*pt1;
  photon_spectrum = (bbs.getBeam1().formfactor(t1)*bbs.getBeam1().formfactor(t1)*pt1*pt1*pt1)/(t1*t1);
                                                                                                                                  
  //cout << "photospeci,form: "<<photon_spectrum<<" "<<formf_nuc(t1,Z)<<endl;
  //       >> Normalize so that photon_spectrum always is < 1
  photon_spectrum = 16.*sqrt(tmin)*photon_spectrum/(3.*sqrt(3.));
                                                                                                                                  
        if( photon_spectrum >  1.0 )
          {
            cout<< "WARNING: photon pt spectrum error "<<"  photon_spectrum="<<photon_spectrum<<endl;
          }
        if( xtest > photon_spectrum )
          goto L202vm;
//      cout <<"PASS1"<<endl;
        phi1 = 2.*StarlightConstants::pi*Randy.Rndom();//random()/(RAND_MAX+1.0);
                                                                                                                                  
        //       >> Draw pt2 (for the Pomeron) according to form factor
        // Separate out pp
        if( bbs.getBeam1().getZin()==1 && bbs.getBeam1().getAin()==1) {
        //dsig/dt= exp(-VMbslope*t)
        xtest = Randy.Rndom();//random()/(RAND_MAX+1.0);
        t2 = (-1./VMbslope)*log(xtest);
        pt2 = sqrt(1.*t2);
        }
        else{
 L203vm:
        xt = Randy.Rndom(); //random()/(RAND_MAX+1.0);
        //dAu--Sergey
   if(bbs.getBeam2().getZin()==1&&bbs.getBeam2().getAin()==2){
        pt2  = 0.8*xt; //it was 0.5,0.8, 1.0  (Sergey)
   }
   else{
        pt2  = 0.5*xt;
    }
        //       >> Check tmin
        tmin = ((Epom/VMgamma_em)*(Epom/VMgamma_em));

        if(tmin > 0.5)
          {
            cout<<" WARNING: tmin= "<<tmin<<endl;
            cout<<" Will pick a new W,Y "<<endl;
            tcheck = 1;
            return;
          }

	xtest = Randy.Rndom();
	
//      xtest = random()/(RAND_MAX+1.0);
        t2 = tmin + pt2*pt2;
        //cout<<"vmom formfsqrdpt2= "<<formf_nuc(t2)*formf_nuc(t2)*pt2<<endl;
   if(bbs.getBeam2().getZin()==1&&bbs.getBeam2().getAin()==2){
        if(1.0 < bbs.getBeam2().formfactor(t2)*pt2)  cout <<"POMERON"<<endl;
        if( xtest > bbs.getBeam2().formfactor(t2)*pt2) goto L203vm;
    }
    else{
        if(1.0 < bbs.getBeam2().formfactor(t2)*bbs.getBeam2().formfactor(t2)*pt2) cout <<"POMERON:Sergey"<<endl;
        if( xtest > bbs.getBeam2().formfactor(t2)*bbs.getBeam2().formfactor(t2)*pt2 )
          goto L203vm;
    }//dAu else end
//      cout<<"PASS2"<<endl;
        }//else end from pp
        phi2 = 2.*StarlightConstants::pi*Randy.Rndom();//random()/(RAND_MAX+1.0);
       
        //       >> Compute px1,py1,px2,py2
        px1 = pt1*cos(phi1);
        py1 = pt1*sin(phi1);
        px2 = pt2*cos(phi2);
        py2 = pt2*sin(phi2);
        
        px = px1 + px2;
        py = py1 + py2;
        
        //       Compute vector sum Pt = Pt1 + Pt2 to find pt for the vector meson
        pt = sqrt( px*px + py*py );
       
        //      I guess W is the mass of the vector meson (not necessarily
        //      on-mass-shell), and E is the energy
        E  = sqrt(W*W+pt*pt)*cosh(Y);
        pz = sqrt(W*W+pt*pt)*sinh(Y);
        //      randomly choose to make pz negative 50% of the time
        if(bbs.getBeam2().getZin()==1&&bbs.getBeam2().getAin()==2){
        pz = -pz;
        }
        else{
        if (Randy.Rndom() >= 0.5) pz = -pz;
        }

//cout<<"vmomentab px= "<<px<<" py= "<<py<<" pz= "<<pz<<" E= "<<E<<endl;

}
void Gammaavectormeson::vmpt(double W,double Y,double &E,double &px,double &py, double &pz,int &tcheck)//may want to remove double w, y
{
//    This function calculates momentum and energy of vector meson
//     given W and Y, including interference.
//     It gets the pt distribution from a lookup table.
        double dW=0.,dY=0.,yleft=0.,yfract=0.,xpt=0.,pt1=0.,ptfract=0.,pt=0.,pt2=0.,theta=0.;
        int IY=0,j=0;
                                                                                                                                  
      dW = (VMWmax-VMWmin)/double(VMnumw);
                                                                                                                                  
      dY  = (VMYmax-VMYmin)/double(VMnumy);
                                                                                                                                  
//  Y is already fixed; choose a pt
//  Follow the approavh in pickwy.f
                                                                                                                                  
// in   fptarray(IY,pt) IY=1 corresponds to Y=0, IY=numy/2 corresponds to +y
                                                                                                                                  
      IY=int(fabs(Y)/dY);//+1;
      if (IY > (VMnumy/2)-1){
            IY=(VMnumy/2)-1;
            }
                                                                                                                                  
//      yfract=ABS(Y)-(2.*IY-numy)*dY
      yleft=fabs(Y)-(IY)*dY;
      yfract=yleft*dY;
                                                                                                                                  
//      cout << "vmpt IY,yleft,yfract,dY= "<<IY<<" "<<yleft<<" "<<yfract<<" "<<dY<<endl;
//      if (fract .lt. 0. .or. fract .gt. 1) write(6,23)fract,Y,IY
//   format(' error in vmpt.  fract=',E11.3,'  Y,IY=',F7.4,2X,I5)
                                                                                                                                  
//       write(6,111)Y,IY,ptmax,dpt
//  111  format(' VMPT. Y,IY,PTMAX,dpt=',F7.4,2X,I5,2X,F7.4,2X,F7.4)
                                                                                                                                  
      xpt=Randy.Rndom(); //random()/(RAND_MAX+1.0);
                                                                                                                                  
      for(j=0;j<VMNPT+1;j++)
      {
              if (xpt < fptarray[IY][j]) goto L60;
      }
L60:
                                                                                                                                  
//  now do linear interpolation - start with extremes
                                                                                                                                  
      if (j == 0){//1) {
               pt1=xpt/fptarray[IY][j]*VMdpt/2.;
               goto L80;
      }
      if (j == VMNPT){//+1){
         pt1=(VMptmax-VMdpt/2.) + VMdpt/2.*(xpt-fptarray[IY][j])/(1.-fptarray[IY][j]);
         goto L80;
      }
                                                                                                                                  
//  we're in the middle
                                                                                                                                  
      ptfract=(xpt-fptarray[IY][j])/(fptarray[IY][j+1]-fptarray[IY][j]);
      pt1=(j+1)*VMdpt+ptfract*VMdpt;
                                                                                                                                  
//  at an extreme in y?
      if (IY == (VMnumy/2)-1){
         pt=pt1;
         goto L120;
      }
                                                                                                                                  
                                                                                                                                  
L80:
//  interpolate in y repeat for next fractional y bin
                                                                                                                                  
      for(j=0;j<VMNPT+1;j++)
      {
              if (xpt < fptarray[IY+1][j]) goto L90;
      }
L90:
                                                                                                                                  
//  now do linear interpolation - start with extremes
                                                                                                                                  
      if (j == 0){//1){
         pt2=xpt/fptarray[IY+1][j]*VMdpt/2.;
         goto L100;
         }
      if (j == VMNPT){//+1){
         pt2=(VMptmax-VMdpt/2.) + VMdpt/2.*(xpt-fptarray[IY+1][j])/(1.-fptarray[IY+1][j]);
         goto L100;
      }
                                                                                                                                  
//  we're in the middle
                                                                                                                                  
      ptfract=(xpt-fptarray[IY+1][j])/(fptarray[IY+1][j+1]-fptarray[IY+1][j]);
      pt2=(j+1)*VMdpt+ptfract*VMdpt;
                                                                                                                                  
L100:
                                                                                                                                  
//  now interpolate in y
                                                                                                                                  
      pt=yfract*pt2+(1-yfract)*pt1;
                                                                                                                                  
L120:
                                                                                                                                  
//  we have a pt
                                                                                                                                  
      theta=2.*StarlightConstants::pi*Randy.Rndom();//(random()/(RAND_MAX+1.0))*2.*pi;
      px=pt*cos(theta);
      py=pt*sin(theta);
                                                                                                                                  
//      I guess W is the mass of the vector meson (not necessarily
//      on-mass-shell), and E is the energy
                                                                                                                                  
        E  = sqrt(W*W+pt*pt)*cosh(Y);
        pz = sqrt(W*W+pt*pt)*sinh(Y);
//      randomly choose to make pz negative 50% of the time
        //if (random()/(RAND_MAX+1.0) >= 0.5) pz = -pz;
	if(Randy.Rndom()>=0.5) pz = -pz;

//cout <<"vmpt2 px,py,pz,E,pt= "<<px<<" "<<py<<" "<<pz<<" "<<E<<" "<<pt<<endl;
//cout <<"vmpt2b ptfract,dpt,j,px,py= "<<ptfract<<" "<<dpt<<" "<<j<<" "<<px<<" "<<py<<endl;
//cout <<"vmpt2c pt1,pt2,yfract,pz,E= "<<pt1<<" "<<pt2<<" "<<yfract<<" "<<pz<<" "<<E<<endl;

}



StarlightConstants::event Gammaavectormeson::produceevent(int &ievent)
{
//returns the vector with the decay particles inside.
        StarlightConstants::event vmeson; //This object will store all the tracks for a single event
        double comenergy = 0.;
        double rapidity = 0.;
        double E = 0.;
        double momx=0.,momy=0.,momz=0.;
        int iFbadevent=0;
	int tcheck=0;
        StarlightConstants::particle ipid = StarlightConstants::UNKNOWN;
                                                                                                                                                           
        double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.;
//this function decays particles and writes events to a file
        //zero out the event structure
                vmeson.numberoftracks=0;
        for(int i=0;i<4;i++)
        {
                vmeson.px[i]=0.;
                vmeson.py[i]=0.;
                vmeson.pz[i]=0.;
                vmeson.fsparticle[i]=StarlightConstants::UNKNOWN;
                vmeson.charge[i]=0;
        }
                 
                pickwy(comenergy,rapidity);
                if(VMinterferencemode==0)momenta(comenergy,rapidity,E,momx,momy,momz,tcheck);
		//bslope...maybe this should be stored in vm instead of sigma...
		if(VMinterferencemode==1)vmpt(comenergy,rapidity,E,momx,momy,momz,tcheck); 
                twodecay(ipid,E,comenergy,momx,momy,momz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
        if (iFbadevent==0&&tcheck==0){
                int q1=0,q2=0;
                                                                                                                                                           
                double xtest = Randy.Rndom();//random()/(RAND_MAX+1.0);
                if (xtest<0.5)
                {
                        q1=1;
                        q2=-1;
                }
                else{
                        q1=-1;
                        q2=1;
                }
                vmeson.numberoftracks=2;//daughter particles are two tracks...
                vmeson.px[0]=px1;
                vmeson.py[0]=py1;
                vmeson.pz[0]=pz1;
                vmeson.fsparticle[0]=ipid;
                vmeson.charge[0]=q1;
                                                                                                                                                           
                vmeson.px[1]=px2;
                vmeson.py[1]=py2;
                vmeson.pz[1]=pz2;
                vmeson.fsparticle[1]=ipid;
                vmeson.charge[1]=q2;
                                                                                                                                                           
                ievent=ievent+1;
                }
        return vmeson;

}

Gammaanarrowvm::Gammaanarrowvm(Inputparameters& input,Beambeamsystem& bbsystem):Gammaavectormeson(input,bbsystem)
{
//Need to make sigma object/run it and read in luminosity tables.
//will just do that outside...of it?
	cout<<"Reading in luminosity tables. Gammaanarrowvm()"<<endl;
	read();
	cout<<"Creating and calculating crosssection. Gammaanarrowvm()"<<endl;
	Narrowresonancesigma sigma(input,bbsystem);
	sigma.crosssectioncalculation(bwnormsave);
	VMbslope=sigma.getbslope();
	
}
Gammaanarrowvm::~Gammaanarrowvm()
{
}

Gammaawidevm::Gammaawidevm(Inputparameters& input,Beambeamsystem& bbsystem):Gammaavectormeson(input,bbsystem)
{
	cout<<"Reading in luminosity tables. Gammaawidevm()"<<endl;
        read();
        cout<<"Creating and calculating crosssection. Gammaawidevm()"<<endl;
        Wideresonancesigma sigma(input,bbsystem);
        sigma.crosssectioncalculation(bwnormsave);
        VMbslope=sigma.getbslope();
}

/*double Gammaawidevm::gettheta()
{
//This depends on the decay angular distribution
//Valid for rho, phi, omega.
cout<<"WIDE THETA TEST LOOK AT ME PLZZZZZ"<<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)
          //  SRK 11/14/2000
//        cout<<"spin1"<<endl;
          dndtheta= sin(theta)*(1.-((cos(theta))*(cos(theta))));
          //    dndtheta = sin(theta)*cos(theta)*cos(theta)
          //if(pid==443||pid==444) dndtheta = sin(theta)*(1.+((cos(theta))*(cos(theta))));will add to psipsiprime code.
          if(xtest > dndtheta)
            goto L200td;
        return theta;
                                                                                                                                                      
}
*/


Gammaawidevm::~Gammaawidevm()
{

}


/*
StarlightConstants::event Gammaanarrowvm::produceevent(int &ievent)
{
//Same as for general wide.  We are just picking W and Y then decaying it...
//Psis will have their own class that inherits from this
}
*/

