#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 const m_pion=0.13957;
    double const c=2.998E8;         // light velocity (m/sec)
    double const tau_pion=2.6E-8;   // pion lifetime (sec)
    double a[4], b[4];
    double Pmax = 0.33329;
    for (int i=0; i<100000; i++) {
//        double pzpi = 20.;
//        double pzpi = 19.*rand01() + 1.;
        double pzpi = pow((1.-3.*rand01()*Pmax),-1./3.);
        double e_pion = sqrt(pzpi*pzpi + m_pion*m_pion);
        double gamma = e_pion / m_pion;
        double decay = -log(rand01()) * gamma * c * tau_pion;
        if( decay < 130.) { 
//        if( decay > 130.) { continue;
          pimuDecay(pzpi, a, b);
//           std::cout << a[0] << " " << b[0] << std::endl;
          double pnu=sqrt(b[1]*b[1]+b[2]*b[2]+b[3]*b[3]);  // p neutrino
          double theta=acos(b[3]/pnu);
          if(theta > 35.e-3 && theta < 53.e-3){
//          if(theta < 0.035 || theta > 0.053) continue;
            std::cout<<pzpi<<"\t"<<decay<<"\t"<<a[0]<<"\t"<<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)) );
}