R__LOAD_LIBRARY ($PLUTOLIBDIR/libPluto.so)

int gamma_pairproduction () 
{
  double mN = 0.938, A = 40, mass_nucl = mN * A; 

  PStaticData *sd = makeStaticData() ;
  sd->AddParticle (-1, "A", mass_nucl );
  sd->SetParticleBaryon ("A", 1);

  PReaction* R = new PReaction (
    "1.0", "g", "A", "A e+ e-", "gamma_pairprod_ttree");
  R->Print ();

  TFile* file_TNtuple = new TFile 
        ("gamma_pairprod_tntuple.root", "RECREATE");
  TNtuple* my_ntuple = new TNtuple 
        ("ntu",  "My Ntuple", "pxep:pyep:pzep:pxem:pyem:pzem:pxA:pyA:pzA" );

  R->Do ("RTD = TMath::RadToDeg()");
  R->Do ("pxep = [e+]->Px() * RTD");
  R->Do ("pyep = [e+]->Py() * RTD");
  R->Do ("pzep = [e+]->Pz() * RTD");
  R->Do ("pxem = [e-]->Px() * RTD");
  R->Do ("pyem = [e-]->Py() * RTD");
  R->Do ("pzem = [e-]->Pz() * RTD");
  R->Do ("pxA  = [A]->Px()  * RTD");
  R->Do ("pyA  = [A]->Py()  * RTD");
  R->Do ("pzA  = [A]->Pz()  * RTD");
  R->Output (my_ntuple);
  R->Loop (10000, 1);

  file_TNtuple->Write ();
  file_TNtuple->Close ();
  return 0;
}