//
//  Sample test program for running EvtGen
//  
//  Created 1999-12-27 Ryd/Lange   
//  Modified nby N. Davidson to demonstrate MC-TESTER with EvtGenLHC

#include "CLHEP/Vector/LorentzVector.h"
#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/EvtRandom.hh"
#include "EvtGenBase/EvtReport.hh"
#include <string>

//StdHep header file
#include "stdhep.h"

//MC-TESTER header file
#include "Generate.h"

using std::endl;
using std::cout;
using namespace CLHEP;

//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[]){

  //Initialize MC-TESTER
  MC_Initialize();

  EvtStdHep evtstdhep;
  EvtParticle *parent;

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

  static EvtId B_plus=EvtPDL::getId(std::string("B+"));

  int nEvents=500000;

  // Loop to create nEvents, starting from an Upsilon(4S)
  for(int i=0;i<nEvents;i++){
    if(i%1000==0)
      cout<<i<<" ("<<i*100.0/nEvents<<"%)"<<endl;

    // Set up the parent particle
    EvtVector4R pInit(EvtVector4R(EvtPDL::getMass(B_plus),0.0,0.0,0.0));
    parent=EvtParticleFactory::particleFactory(B_plus,pInit);

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

    evtstdhep.init();
    
    // Write out the results
    parent->makeStdHep(evtstdhep);
    
    if(i==0) //print 1st event
      report(INFO,"EvtGen") << evtstdhep << endl;

    //Copy results to hepevt common block from stdhep
    hepevt_.nevhep=i;   //hepevt_ is defined in stdhep.h
    hepevt_.nhep=evtstdhep.getNPart();    
    for(int pcle = 0; pcle < evtstdhep.getNPart(); pcle++){
      hepevt_.isthep[pcle]=evtstdhep.getIStat(pcle);
      hepevt_.idhep[pcle]=evtstdhep.getStdHepID(pcle);

      hepevt_.jmohep[pcle][0]=evtstdhep.getFirstMother(pcle)+1; //first pcle at 1 not 0.
      hepevt_.jmohep[pcle][1]=evtstdhep.getLastMother(pcle)+1;
      hepevt_.jdahep[pcle][0]=evtstdhep.getFirstDaughter(pcle)+1;
      hepevt_.jdahep[pcle][1]=evtstdhep.getLastDaughter(pcle)+1;

      /* 4-Momentum, mass */
      hepevt_.phep[pcle][0]=evtstdhep.getP4(pcle).get(1); //px
      hepevt_.phep[pcle][1]=evtstdhep.getP4(pcle).get(2); //py
      hepevt_.phep[pcle][2]=evtstdhep.getP4(pcle).get(3); //pz
      hepevt_.phep[pcle][3]=evtstdhep.getP4(pcle).get(0); //e
      hepevt_.phep[pcle][4]=evtstdhep.getP4(pcle).mass();

      /* Vertex information */
      hepevt_.vhep[pcle][0]=evtstdhep.getX4(pcle).get(1); //x
      hepevt_.vhep[pcle][1]=evtstdhep.getX4(pcle).get(2); //y
      hepevt_.vhep[pcle][2]=evtstdhep.getX4(pcle).get(2); //z
      hepevt_.vhep[pcle][3]=evtstdhep.getX4(pcle).get(0); //t
    }

    //Analyse the B+ decay
    MC_Analyze(521);

    //clean up
    parent->deleteTree();  
  }
  
  //Finalize MC-TESTER
  MC_Finalize();

  return 1;
}