#include <iostream>
#include <iomanip>

#include "hepevt.hh"
#include "TParticle.h"
#include "TPrimaryInteraction.hh"


ClassImp(TPrimaryInteraction)

//_____________________________________________________________________________
TPrimaryInteraction::TPrimaryInteraction() {
  fParticleList = new TClonesArray("TParticle",100);
}

//_____________________________________________________________________________
TPrimaryInteraction::TPrimaryInteraction(int n) {
  
				// sometimes you know how many particles are 
				// there in the event

  fParticleList = new TClonesArray("TParticle",n);
}

//_____________________________________________________________________________
TPrimaryInteraction::~TPrimaryInteraction() {
  if (fParticleList) {
    fParticleList->Delete();
    delete fParticleList;
  }
}

//_____________________________________________________________________________
void TPrimaryInteraction::InitFromHepevt() {
  fParticleList->Clear();
  InitHepevt();
  Int_t numpart = gHepevt->NHEP;
  for (Int_t i = 0; i<numpart; i++) {
    new ((*fParticleList)[i]) TParticle(gHepevt->IDHEP[i],
					gHepevt->ISTHEP[i],
					gHepevt->JMOHEP[i][0]-1,
					gHepevt->JMOHEP[i][1]-1,
					gHepevt->JDAHEP[i][0]-1,
					gHepevt->JDAHEP[i][1]-1,
					gHepevt->PHEP[i][0],
					gHepevt->PHEP[i][1],
					gHepevt->PHEP[i][2],
					gHepevt->PHEP[i][3],
					gHepevt->VHEP[i][0],
					gHepevt->VHEP[i][1],
					gHepevt->VHEP[i][2],
					gHepevt->VHEP[i][3]);
  }
}

//_____________________________________________________________________________
void TPrimaryInteraction::FillHepevt() {
  TLorentzVector v;
  TLorentzVector mom;
  gHepevt->NHEP = NParticles();
  for (Int_t i = 0; i<gHepevt->NHEP; i++) {
    TParticle* p = Particle(i);
    gHepevt->IDHEP [i] = p->GetPdgCode();
    gHepevt->ISTHEP[i] = p->GetStatusCode();
    gHepevt->JMOHEP[i][0] = p->GetFirstMother  ()+1;
    gHepevt->JMOHEP[i][1] = p->GetSecondMother ()+1;
    gHepevt->JDAHEP[i][0] = p->GetFirstDaughter()+1;
    gHepevt->JDAHEP[i][1] = p->GetLastDaughter ()+1;

    p->Momentum(mom);
    gHepevt->PHEP[i][0]   = mom.Px();
    gHepevt->PHEP[i][1]   = mom.Py();
    gHepevt->PHEP[i][2]   = mom.Pz();
    gHepevt->PHEP[i][3]   = mom.Energy();
    gHepevt->PHEP[i][4]   = mom.M();
    p->ProductionVertex(v);
    gHepevt->VHEP[i][0]   = v.X();
    gHepevt->VHEP[i][1]   = v.Y();
    gHepevt->VHEP[i][2]   = v.Z();
    gHepevt->VHEP[i][3]   = v.T();
  }
}

//_____________________________________________________________________________
void TPrimaryInteraction::Print(Option_t* opt) const {
  // `a' = "all": print all the particles

  printf("TPrimaryInteraction %20s NParticles = %5i\n",GetName(),NParticles());

  if (strchr(opt,'a')) {
    for (int i=0; i<NParticles(); i++) {
      printf("%5d ",i);
      Particle(i)->Print();
    }
  }
}