R__LOAD_LIBRARY ($PLUTOLIBDIR/libPluto.so)

int elastic_AA_collision () 
{
  double A1 = 40, A2 = 40;
  double mN = 0.9315, m_nucl1 = mN * A1, m_nucl2 = mN * A2;
  char Tbeam[10];  sprintf (Tbeam, "%f", 1.00 * A1 );

  /* Simulate Ca+Ca -> Ca+Ca elastic collision */

  PStaticData *sd = makeStaticData() ;
  sd->AddParticle (-1, "A1", mass_nucl1);     // -1 : PID next in list
  sd->AddParticle (-1, "A2", mass_nucl2);    

  PReaction *R = new PReaction (Tbeam,        /* Ekin of beam particle */
		"A1", "A2", 	              /* Substrates */
		"A1 A2"   ,                   /* Products   */
		"elastic_AA");                /* TTree output file name */
  R->Print ();

  TFile *file_TNtuple = new TFile ("elastic_AA.root", "RECREATE");
  TNtuple *my_ntuple = new TNtuple ("ntu", 
          "Elastic", "th1:ph1:th2:ph2:opang" );
  R->Do ("PI = TMath::RadToDeg()");
  R->Do ("th1 = [A1]->Theta()       * PI");
  R->Do ("ph1 = [A1]->Phi()         * PI");
  R->Do ("th2 = [A2]->Theta()       * PI");
  R->Do ("ph2 = [A2]->Phi()         * PI");
  R->Do ("opang = [A1]->Angle([A2]) * PI");
  R->Output (my_ntuple);

  R->Loop (10000, 1);
  file_TNtuple->Write ();
  file_TNtuple->Close ();
  return 0;
}