#include "LCDreadStdFile.h"

#include "LCDMcPart.h"

#include <iostream.h>
#include <assert.h>

//______________________________________________________________________
// 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(LCDreadStdFile)
 LCDreadStdFile::LCDreadStdFile(Char_t* filein, LCDEvent* event) {
  inputfile = filein;
  istr = 0;

  if (event == 0) {
    m_event = new LCDEvent();
    m_ownevent=kTRUE;
  } else {
    m_event = event;
    m_ownevent=kFALSE;
  }

  assert(! StdHepXdrReadInit(inputfile, 1, istr));
}
//
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 LCDreadStdFile::~LCDreadStdFile() {
  StdHepXdrEnd(istr);
  printf("Closing input file %s n",inputfile);
  if (m_ownevent) {
    m_event->Delete();
    delete m_event;
  }
}
//
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 Int_t LCDreadStdFile::readEvent() {
  // reading the binary file to fill the common block

  hepevt_.nhep = 0; Int_t 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(" unexpected end of file after %d eventsn",hepevt_.nhep);
    return 0;
  }
  return (hepevt_.nhep == 0);
}


 Int_t LCDreadStdFile::GetEvent(LCDEvent* event) {
  // From stdHEP file representation to MCpart instances
  if( readEvent() ) return 0;
  if (event) {
    MakeMcPart(event);
  } else if (m_event) {
    MakeMcPart(m_event);
  } else {
    fprintf(stderr,
	    "LCDreadStdFile::GetEvent Pointer of LCDEvent is not given!n");
    return 0;
  }
  return 1;

}
//
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 Int_t LCDreadStdFile::MakeMcPart(LCDEvent* event){
  // Makes McPart objects from the stdHEP common

  Int_t          npart = hepevt_.nhep;
  Int_t          i        ;
  LCDMcPart*     part     ;
  Int_t          status   ;
  Int_t          id       ;
  Double_t       charge   ;
  TVector3       position ; 
  TLorentzVector momentum4;

  TClonesArray* mclist=event->MCparticles();
  mclist->Delete();

  for(i = 0; i < npart; i++) {
    if (i > 3999) continue;
    //if (npart > 3999) continue;
    // create and fill the McPart
    new((*mclist)[i]) LCDMcPart();
    part = (LCDMcPart*)mclist->At(i);

    status = hepevt_.isthep[i];
    id     = hepevt_.idhep[i];
    charge = partchg(id);

    // 09/22/2000 T.Abe mm to cm conversion.
    position.SetXYZ(hepevt_.vhep[i][0]/10.0,
		    hepevt_.vhep[i][1]/10.0,
		    hepevt_.vhep[i][2]/10.0);
    momentum4.SetPxPyPzE(hepevt_.phep[i][0],
			 hepevt_.phep[i][1],
			 hepevt_.phep[i][2],
			 hepevt_.phep[i][3]);

    part->SetStatus(status);
    part->SetParticleID(id);
    part->SetCharge(charge);
    part->SetPosition(position);
    part->Set4Momentum(momentum4);

  }

  // Fill parent and daughter information
  Int_t i_parent;
  LCDMcPart* parent;
  for(i = 0; i < npart; i++){
    part = (LCDMcPart*)(event->MCparticles()->At(i));

    i_parent   =  hepevt_.jmohep[i][0] - 1;
    part->SetParentIdx(i_parent);      

    if (i_parent >= 0) {
      parent=(LCDMcPart*)(event->MCparticles()->At(i_parent));
      parent->SetTermPosition(*part->GetPositionPtr());
    }
  }

  return 1;
}