//
//  Sample test program for running EvtGen
//  
//  Created 17/10/2006 by Stefano Spataro
//
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtPatches.hh"
#include <iostream>
#include "EvtGenBase/EvtParticleFactory.hh"
#include "EvtGenBase/EvtStdHep.hh"
#include "EvtGen/EvtGen.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtVector4R.hh"
#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtReport.hh"
#include <string>
#include <fstream>
#include <stdio.h>
#include <stdlib.h>
using std::endl;
using std::ofstream;
using std::cout;


//Define random number fcn used by Jetset
extern "C" {
  extern float rlu_();
  extern float begran_(int *);
}

float rlu_(){
  return EvtRandom::Flat();
}

float begran_(int *){
  return EvtRandom::Flat();
}


int main(int argc, char* argv[]){

  EvtStdHep evtstdhep;
  EvtParticle *parent;

  if (argc<3) {
    cout << "\nUSAGE: myEvtGen <particle> <dec-file> <# events> <pbar-mom/cms-energy>\n" << endl;
    cout << "  <particle> = particle type to decay, e.g. 'eta_c', 'pbarpSystem' etc."<<endl;               //argv[1]
    cout << "  <dec-file> = EvtGen decay file (.DEC) to use; see directory 'test' for examples"<<endl;     //argv[2]
    cout << "  <# events> = number of events to produce; default value = 10"<<endl;                    //argv[3]
    cout << "  <pbar-mom> = (>0) momentum of the pbar beam; (<0) negativ cms energy; default value = mass of <particle>"<<endl;  //argv[4]
    cout << "               mandatory, when <particle> = pbarpSystem\n"<<endl;
    cout << "There will be stored just events with output-particles below 50mrad theta!\n"<<endl;
    cout << "Output is stored in file 'output_cut.evt'.\n\n"<<endl;
    return 0;
  }

  //Initialize the generator - read in the decay table and particle properties
  EvtGen myGenerator("DECAY.DEC","evt.pdl");

  //If I wanted a user decay file, I would read it in now.
  myGenerator.readUDecay(argv[2]);

  static EvtId PART=EvtPDL::getId(std::string(argv[1]));

  int number=10;
  if (argc>=4) number=atoi(argv[3]);

  if (std::string(argv[1])=="pbarpSystem" && argc<5)
  {
    cout <<"\n******  FATAL ERROR: <particle> is 'pbarpSystem'; MUST give pbar momentum or cms energy!\n\n"<<endl;
    return 0;
  }
  //else if (std::string(argv[1])!="pbarpSystem" && argc>=5)
  // {
  //  cout <<"\n****** WARNING: overriding given momentum, setting cms energy to mass of "<<argv[1]<<".\n"<<endl;
  // }

  double val=-3.0969;
  double P = 0.0;
  double E = 0.0;
  double mp=0.93827;

  if (argc>=5) 
    val=atof(argv[4]);
  else
    val=-EvtPDL::getMass(PART);
  
  // val is the momentum of the pbar beam
  if (val>0){  
    P = val;
    E = mp+sqrt(P*P+mp*mp);
  }
  else  //val is -E_cm
  {
    val=-val;
    E = val*val/(2*mp);
    P = sqrt(E*E-val*val);
  }
  
  cout <<"\n\n############# Generating with following conditions:\n\n";
  cout <<"particle       : '"<<argv[1]<<"'"<<endl;
  cout <<"decay file     : "<<argv[2]<<endl;
  cout <<"incident 4-mom : ("<<E<<", 0, 0, "<<P<<"), m = "<<sqrt(E*E-P*P)<<endl;
  cout <<"# events       : "<<number<<"\n\n######################\n\n"<<endl;


  // Open the output file  in the format requested by PandaROOT
  ofstream out;
  out.open("output_cut.evt");
 
  // Loop to create nEvents, starting from an Upsilon(4S)
  int i;
  unsigned long j=0;
  for(i=0;i<number;i++){
    j++;
    // Set up the parent particle

    EvtVector4R pInit(E,  0.0001, -0.0000,  P);
    parent=EvtParticleFactory::particleFactory(PART,pInit);
    parent->setDiagonalSpinDensity();  

    // Generate the event
    myGenerator.generateDecay(parent);

    // Write out the results
    evtstdhep.init();
    parent->makeStdHep(evtstdhep);

    //CUT to theta-out < 3° ---------------------------------------------------
    bool isNice=false;
    EvtParticle* iter = parent->nextIter();
    while(iter){
      EvtVector4R imp = iter->getP4Lab();
      imp.get(0);
      double theta=1000*fabs(acos( sqrt(imp.get(3)*imp.get(3)/(imp.get(1)*imp.get(1)+imp.get(2)*imp.get(2)+imp.get(3)*imp.get(3))) ));
      if(theta<50)
        isNice=true;
      iter = iter->nextIter();
    }
    //CUT end -----------------------------------------------------------------
    
    // Write the output file
    if(isNice){
      out << i << "\t" << evtstdhep.getNPart();
      out <<evtstdhep<<endl;
    }else{
      i--;
    }

    //print out some status info
    if (i<10) report(INFO,"EvtGen") << "event Number\t"<< i << evtstdhep << endl;
    if (!((i+1)%100))  report(INFO,"EvtGen") << "event Number\t"<<i+1<<endl;

    parent->deleteTree();  
  }	
  out.close();
  cout <<"\n\n############# Generating done! ############\n\n";
  cout <<"# simulated decay's : "<<j<<endl;
  return 1;
}