#include <stdlib.h>
#include <iostream>
#include <math.h>
void pimuDecay(double Epion, double* a, double* b);
double rand01();
double qcm(double ma, double mb, double mc);

int main()
{
    double a[4], b[4];
    for (int i=0; i<100000; i++) {
        double pzpi = 20.;
//      pzpi = 19.*rand01() + 1.;
        pimuDecay(pzpi, a, b);
           std::cout << a[0] << " " << b[0] << std::endl;
     }
    return 0;
} 

void pimuDecay(double Epion, double* a, double* b)
{
    double acm0, acmx, acmy, acmz, acm;
    double bcm0, bcmx, bcmy, bcmz;
    double costhetacm, sinthetacm, phicm;
    double const mpion=0.13957, ma=0.10566, mb=0.;
    double const pi = 3.141592;
    double gamma, beta;

// in the center of mass system (CMS) 
    acm = qcm(mpion,ma,mb);
// phase-space sampling d(costhetacm)d(phi)
    costhetacm = -1.0 + 2.0 * rand01();
    sinthetacm = sqrt (1.0- costhetacm*costhetacm);
    phicm = -pi *+ 2.*pi*rand01();
    acm0 = sqrt(ma*ma + acm*acm);
    acmx = acm * sinthetacm * cos(phicm);
    acmy = acm * sinthetacm * sin(phicm);
    acmz = acm * costhetacm; 
    bcm0 = sqrt(mb*mb + acm*acm);
    bcmx = -acmx;
    bcmy = -acmy; 
    bcmz = -acmz; 
     
// Lorentz transformation from CMS to laboratory system
    gamma = Epion/mpion;
    beta = sqrt(1-1/pow(gamma,2));
    a[0] = gamma * (acm0 + beta * acmz); 
    a[1] = acmx;
    a[2] = acmy;
    a[3] = gamma * (acmz + beta * acm0); 
    b[0] = gamma * (bcm0 + beta * bcmz); 
    b[1] = bcmx;
    b[2] = bcmy;
    b[3] = gamma * (bcmz + beta * bcm0); 
    
  return ;
}

double qcm(double ma, double mb, double mc)
{
    return  sqrt(pow((ma*ma-mb*mb-mc*mc),2)-4.*mb*mb*mc*mc)/(2*ma);
}

double rand01()
{
  return  ( (double)rand() / ((double)(RAND_MAX)+(double)(1)) );
}