R__LOAD_LIBRARY ($PLUTOLIBDIR/libPluto.so)

int thermal_pi0_gg ()
{
   //
   // Thermal source with temperature(s) T1 (T2), frac*f(T1) + (1-frac)*f(T2)
   //
   // optional blast, optional polar anisotropies (A2,A4), optional flow (v1,v2)
   //

   PUtils* putils = new PUtils ();                   // Randomization of seed
   putils->SetSeed (0);
   gRandom->SetSeed();

   /***    Properties of particle's production from Fireball  ***/

   PFireball* source_pi0 = new PFireball (
	"pi0",  /* type of source particle            */
	0.040,  /* Beam Ekin (per nucleon) of AA collision             [GeV/A]     */
	0.015,  /* Temperature of thermal source of particle           [GeV]       */
	0.,     /* Possible addition of 2. source -> its temperature   [GeV]       */
	1.,     /* Weight of 1st source               */ 
	0.,     /* Possible addition of collective radial expansion,  Beta (v/c)   */
	0.,     /* Possible addition of polar angle anisotropy,  a2 coefficient    */
	0.,     /*                                               a4 coefficient    */
	0.,     /* Possible addition of anisotropy in phi angle, v1 coefficient    */
	0.      /*                                               v2 coefficient    */
       );
   source_pi0 ->Print ();

   /***        Fireball ---> pi0 production         ***/

   PParticle* pi0 = new PParticle ("pi0");
   PParticle* parttable_fireball_pi0[] = {source_pi0, pi0};
   PChannel*  channel_fireball_pi0   = new PChannel (parttable_fireball_pi0, 1, 1);

   /***    pi0 ---> gamma gamma decay   (BR = 99%)  ***/

   PParticle* g1  = new PParticle ("g");
   PParticle* g2  = new PParticle ("g");
   PParticle* parttable_pi0_gammagamma [] = {pi0, g1, g2};
   PChannel*  channel_pi0_gammagamma = new PChannel (parttable_pi0_gammagamma, 2, 1);

   /***    Building up a complete reaction          ***/

   PChannel*  table_of_channels[]   = { channel_fireball_pi0 , channel_pi0_gammagamma };

   PReaction* reaction_thermal_pi0 = new PReaction (
	table_of_channels, 		    /* Table of channels           */
	"thermal_pi0_gg_Eb40MeV_T15MeV",    /* Output file name            */
        2,                                  /* Number of channels involved */
        1, 0, 0, 0                          /* Flags                       */
       );

   reaction_thermal_pi0 ->Print ();
   reaction_thermal_pi0 ->loop (10000, 1);

   return 0;
}