R__LOAD_LIBRARY ($PLUTOLIBDIR/libPluto.so)

int pp_pppi0_gg_long () 
{
    //   p (Tb = 0.2797) + p  --->  composite_particle
    PParticle p1 ("p", 0.2797) , p2 ("p");
    PParticle composite = p1 + p2;

    //    composite_particle  --->  p + p + pi0

    PParticle p3 ("p"), p4 ("p"), pi0 ("pi0");

    PParticle* parttable_pi0[] = { &composite , &pi0 , &p3 , &p4 };

    PChannel* pi0_prod = new PChannel ( parttable_pi0 , 3 );

                                                      // "3" = no of decay particles
                                                      // other flags: obsolete
    //    Decay   pi0 ---> gamma+gamma

    PParticle g1 ("g") , g2 ("g");

    PParticle* parttable_gg[] = { &pi0, &g1 , &g2 };

    PChannel* pi0_ggdecay = new PChannel ( parttable_gg , 2 );

    //    Constructing the full reaction

    PChannel* table_of_channels[] = { pi0_prod , pi0_ggdecay };

    PReaction* my_reaction = new PReaction 
               ( table_of_channels , "pp_pppi0_gg_long", 2 , 1 , 0 , 0 , 0 );
/*
   Table of channels
   Output file name
   Number of channels in table
   Options:
     f0: 0=Tracked, 1=All particles in ROOT output file_name.root
     f1: unused
     f2: 0=Do not calculate, 1=Calculate production vertices for product particles
     f3: 0=Do not produce, 1=Produce GEANT ascii output file_name.evt (tracked particles always)
*/

    my_reaction->Print ();

    //    Generate events!

    my_reaction->Loop (10000, 1);			/* No of events, status @ end */

    return 0;
}