// #-> charm.py -- simulate j/psi electroproduction at JLab
// c++ version i hope

// from Numeric import *  [ left over from python version ]
// from vec4 import *
// import RNG

// #-> simulation driver parameters

//main(int argc, char **argv)
int jeff1()
{
  float mj = 3096.88;      //   # J/psi mass in MeV
  float me =    0.510999;  //   # electron mass
  float mp =  938.272;

  float beamen = 11.0 * 1000.0; // # GeV * 1000 MeV/GeV

  float ll_the = (3.1415926 / 180) *  1.0; // last number is e- angle LL in deg
  float ul_the = (3.1415926 / 180) * 20.0; // last number is e- angle UL in deg

// above is stored in radians since we will need to take cosines
// below is not

  float ll_phe = -93.0; // LL on e- azimuthal angle (degrees now)
  float ul_phe = -87.0;

  float ll_pe = 2275; // LL on scat electron momentum (MeV/c)
  float ul_pe = 2780;

  int howmany = 10000; //  # desired number of events

// #-> set up random generators for electron.  strategy: create sample
// #-> arrays at start of program, just take values from them during
// #-> execution

  TF1 *thegen = new TF1("thegen","sin(x)",ll_the,ul_the);
  TF1 *phegen = new TF1("phegen","1.0",   ll_phe,ul_phe);
  TF1 *pegen  = new TF1("pegen", "1.0",   ll_pe, ul_pe );

  // #-> set up four-momenta for beam and target, which should not change

  TLorentzVector beam;  beam.SetVectM(TVector3(0,0,beamen),me);
  TLorentzVector targ(0.0,0.0,0.0,mp);

  TFile hfile("charm.root","RECREATE","Charm");

  TTree *tree = new TTree("T","Charm Electroproduction");

  int nev = 0;
  TLorentzVector scat(1,1,1,1);
  TLorentzVector X;
  
  TLorentzVector *pbeam = &beam;
  TLorentzVector *ptarg = &targ;
  TLorentzVector *pscat = &scat;
  TLorentzVector *pX    = &X;

  tree->Branch("beam","TLorentzVector",&pbeam,32000,0);
  tree->Branch("targ","TLorentzVector",&ptarg,32000,0);
  tree->Branch("scat","TLorentzVector",&pscat,32000,0);
  tree->Branch("ppsi","TLorentzVector",&pX,32000,0);

  for ( int nev = 1; nev <= howmany; nev++) {
    float thtmp = thegen.GetRandom();
    float phtmp = phegen.GetRandom();
    float ptmp  = pegen.GetRandom();

    scat.SetVectM(TVector3(ptmp*sin(thtmp)*cos(phtmp),
                           ptmp*sin(thtmp)*sin(phtmp),
                           ptmp*cos(thtmp)),
                  me);

    X = beam + targ - scat;

    tree->Fill();
  }

  hfile.Write();
  hfile.Close();

  return 0;
}