R__LOAD_LIBRARY ($PLUTOLIBDIR/libPluto.so)

int compton () 
{ 
  /*****   Setting up the Compton reaction:

         gamma (10 MeV) + e-  --->   gamma + e-
                                                  *****/

  PReaction *R = new PReaction ("0.010",       /* Ekin of beam particle */
				"g", "e-", 	         /* Substrates */
				"g e-"   ,               /* Products   */
				"compton_ttree");        /* TTree output file name */
  R->Print ();



  /***   TNtuple as a quick Event Database   ***/

  TFile *file_TNtuple = new TFile ("compton_tntuple.root", "RECREATE");

  TNtuple *my_ntuple = new TNtuple ("ntu",  "Monitor ntuple", 
				    "thgam:phgam:thel:phel:opang" );

  gDirectory->pwd();

  R->Do ("RTD = TMath::RadToDeg() ");
  R->Do ("thgam = [g]->Theta()     * RTD");
  R->Do ("phgam = [g]->Phi()       * RTD");
  R->Do ("thel  = [e-]->Theta()    * RTD");
  R->Do ("phel  = [e-]->Phi()      * RTD");
  R->Do ("opang = [g]->Angle([e-]) * RTD");

  /***   Adding TNtuple to the list of filled objects during loop  ***/

  R->Output (my_ntuple ,
             "if (thgam>20 && thgam < 50) " );              /* Optional filter via "if" */


  /*****   Additional filtering of the basic TTree      *****/
  /*       Each #variable is treated as filter. 0 = Fail    */

  R->Do ( "#thetaIsAccepted = 1; if (thgam<20 || thgam > 50) ; #thetaIsAccepted  = 0" );


  /***** Quick histogram. For example, Opening angle, with some filter.  ***/

  TH1F *my_histo = new TH1F ("myhisto", "Opening angle", 90, 0., 90.);

  R->Do ( my_histo, "if thgam > 20 && thgam < 50 ; _x = opang" ); 

					                 /*  _x is the histogrammed variable */
						         /*  Possible filtering via "if"     */

  /***   Run the reaction   ***/

  R->Loop (10000, 1);

  gDirectory->cd ("compton_tntuple.root:/");
  my_histo ->Write();
  my_ntuple->Write();
  file_TNtuple->Close ();

  return 0;
}