for (int i=0; i<100000; i++) {
        double pzpi = 19.*rand01() + 1.;
        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.) continue;
        pimuDecay(pzpi, a, b);
        std::cout<<pzpi<<"\t"<<decay<<"\t"<<a[0]<<"\t"<<b[0]<< std::endl;
     }