{
   const Int_t nevt = 1000;
   const Float_t kPi0Mass = 0.1349766; // GeV/c^2

   // pt distribution for pi0 particles
   TF1 *pt_func = new TF1("pt_func","x*pow(1+x,[0])",0,5);
   pt_func->SetParameter(0,-5);

   TFile *fout = new TFile("decay_info.root","RECREATE");
   
   // Make variables and set up tree
   Float_t E, px, py, pz;
   TTree *dec_tree = new TTree("dec_tree","Decay info");
   dec_tree->Branch("E", &E, "E/F");
   dec_tree->Branch("px", &px, "px/F");
   dec_tree->Branch("py", &py, "py/F");
   dec_tree->Branch("pz", &pz, "pz/F");

   // Event loop
   for (Int_t ievt = 0; ievt < nevt; ievt++) {

     // Generate pi0 momentum
     Float_t y = gRandom->Rndm()*2.0 - 1.0;  // rapidity; uniform in [-1,1]
     Float_t phi = gRandom->Rndm()*TMath::Pi()*2.0;  // uniform [0, 2pi]
     Float_t pt = pt_func->GetRandom();  // according to pt_func
     Float_t mt = TMath::Sqrt(kPi0Mass*kPi0Mass + pt*pt);

     // Calculate px, py, pz from y, phi, pt
     px = pt * TMath::Cos(phi);
     py = pt * TMath::Sin(phi);
     pz = mt * TMath::SinH(y);
     E = mt * TMath::CosH(y);

     // Add decay generation code here
 
     dec_tree->Fill();
   }
   dec_tree->Write();

}