// 1999/02/26   Adapted for use in lcd simulations from Gismo 2,
// /afs/cern.ch/atlas/software/cvs/offline/arve/gismo/readStdFile.cxx,v 1.1.1.1
// 1997/08/14 14:04:13 lat 
//


#include "Event.h"
#include "TMap.h"
#include "readStdFile.h"

#ifdef WIN32
#define hepevt_ HEPEVT
#endif

struct hepevt HEPEVT;

//______________________________________________________________________
// readStdFile
//
// This class is one type of event source.  It can get an event from a stdHEP
// file and feed the MC particle information into instances of MC

ClassImp(readStdFile)
 readStdFile::readStdFile(char* filein) {
  inputfile = filein;
  istr = 0;

  assert(! StdHepXdrReadInit(inputfile, 1, istr));
}
//
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
 readStdFile::~readStdFile()
{
   StdHepXdrEnd(istr);
   printf("Closing input file %c n",inputfile);
}
//
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
 int readStdFile::readEvent() {
  // reading the binary file to fill the common block

  hepevt_.nhep = 0;
  int lbl = 1;	

  StdHepZero();
	
  ierr = StdHepXdrRead(&lbl,istr);

  if(lbl == 100){
    printf(" Reading StdHep begin run record.n");
    ierr = StdHepXdrRead(&lbl,istr);
  } 
  if (ierr == 0)
    printf(" at event %d with %d particlesn", hepevt_.nevhep, hepevt_.nhep);
  else {
    printf(" unexpected end of file after %d eventsn",hepevt_.nhep);
    return 0;
  }
  printf(" 1 event readn");
  return (hepevt_.nhep == 0);
}


 int readStdFile::getEvent(Event* event) {
  // From stdHEP file representation to MCpart instances

  if( readEvent() )
         return 0;
  makeMCs(event);
return 1;
}
//
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
 int readStdFile::makeMCs(Event* event){
  // This method makes McPart objects from the stdHEP common
  // Someday if we have an MCdiag which also wants to output to the
  // stdHEP common this method should be moved to a utility base
  // class from which readStdFile and MCdiag could both inherit.

  m_event = event;

  float pzero[4] = {0.};
  float vzero[4] = {0.};
  McPart* parents[4000] = {0};


  // create a dummy first particle at IP

  McPart* firstMC = new McPart();
  int fakeID = 999999, fakeStat = 0;
  float fakeChg = 0.;

  firstMC->SetUpParticle(fakeID,     // fake particle type
			 fakeChg,     // charge
			 fakeStat,    // status
			 0,           // no parent particle
			 pzero,        // set momentum to zero
			 vzero);        // set position to zero


  m_event->MCparticles()->Add(firstMC);

  parents[0] = firstMC;

  unsigned numsecond = unsigned(hepevt_.nhep);
  for( unsigned i = 0; i < numsecond; i++)
  {										  										  //for different generator output.
    int statusCode = hepevt_.isthep[i];
    if( statusCode==0 ) continue;

    float mass = hepevt_.phep[i][4];
    int parentIndex = static_cast<int>(hepevt_.jmohep[i][0]);
    if (parentIndex < 0) parentIndex = 0;
			  
    float pvect[4];
    float vvect[4];
    pvect[0] = static_cast<float>(hepevt_.phep[i][0]);
    pvect[1] = static_cast<float>(hepevt_.phep[i][1]);
    pvect[2] = static_cast<float>(hepevt_.phep[i][2]);
    pvect[3] = static_cast<float>(hepevt_.phep[i][3]);
    vvect[0] = static_cast<float>(hepevt_.vhep[i][0]);
    vvect[1] = static_cast<float>(hepevt_.vhep[i][1]);
    vvect[2] = static_cast<float>(hepevt_.vhep[i][2]);
    vvect[3] = static_cast<float>(hepevt_.vhep[i][3]); // not really used

    int id = hepevt_.idhep[i];

    // create and fill the McPart

    McPart* part = new McPart();

    McPart* p = parents[parentIndex];

    float chg = stdchg_(&id);

    part->SetUpParticle(id,         // particle type
			chg,         // charge
			statusCode,  // status from generator
			p,            // parent particle
			pvect,        // momentum
			vzero);        // set position to zero

    // set the parent's termination point to be the initial position of
    // its daughter.

    if ( p ) p->SetPosition(vvect);

    // add this guy to parentMap

    parents[i+1] = part;

    // and to the list of McPart's
    m_event->MCparticles()->Add(part);
  }

  return 1;
}