UtilHepevt.cxx

Go to the documentation of this file.

#include "Util/UtilHepevt.h"
#include "Util/UtilPDG.h"

// for ::open, ::write, ::close
#include <unistd.h>
// for ssize_t
#include <sys/types.h>
#include <sys/stat.h>
// for scanf
#include <stdio.h>

#include <iostream>
#include <iomanip>
#include <vector>

#include "TMath.h"   // for units
#include "TClonesArray.h"
#include "TParticle.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TApplication.h"

// Hepevt.h is in $ROOTSYS/eg
#include "Hepevt.h"
const int mxnhep = 4000;  // this should have been part of Hepevt.h

int  UtilHepevt::getNHep() { return HEPEVT.nhep; }
int  UtilHepevt::getNevHep() { return HEPEVT.nevhep; }
void UtilHepevt::setNevHep(int evtnum) { HEPEVT.nevhep = evtnum; }

void UtilHepevt::clearCommon(bool fast) 
{
  HEPEVT.nevhep = 0;
  HEPEVT.nhep   = 0;
  if (fast) return;
  for (int i=0; i<mxnhep; ++i) {
    HEPEVT.isthep[i]     = 0;
    HEPEVT.idhep[i]      = 0;
    HEPEVT.jmohep[i][0]  = 0;
    HEPEVT.jmohep[i][1]  = 0;
    HEPEVT.jdahep[i][0]  = 0;
    HEPEVT.jdahep[i][1]  = 0;
    HEPEVT.phep[i][0]    = 0;
    HEPEVT.phep[i][1]    = 0;
    HEPEVT.phep[i][2]    = 0;
    HEPEVT.phep[i][3]    = 0;
    HEPEVT.phep[i][4]    = 0;
    HEPEVT.vhep[i][0]    = 0;
    HEPEVT.vhep[i][1]    = 0;
    HEPEVT.vhep[i][2]    = 0;
    HEPEVT.vhep[i][3]    = 0;
  }
}

TClonesArray* UtilHepevt::getEmptyStdHep(int dfltsiz)
{
  if (dfltsiz<0) dfltsiz = mxnhep;
  TClonesArray* tc = new TClonesArray("TParticle",dfltsiz);
  tc->SetName("StdHep");
  return tc;
}

int UtilHepevt::unpackCommon(TClonesArray* stdhep_ptr,
                               cvu_t vtxCvu, cvu_t tCvu)
{ 
  // Copy info from /HEPEVT/ common block into the supplied
  // TClonesArray of TParticles.

  if (stdhep_ptr == 0) return -1;  // user failed to supply an array
  TClonesArray &stdhep = *stdhep_ptr;

  double vtxscale = unitScale(vtxCvu);
  double tscale   = unitScale(tCvu);

  stdhep.Clear();  // clear what's already there

  int numpart = HEPEVT.nhep;
  for (int i = 0; i<numpart; ++i) {
    // force PDG # for ions into chosen canonical form
    int ipdg = UtilPDG::stdIonNumber(HEPEVT.idhep[i]);
    new(stdhep[i]) 
      TParticle(
                ipdg, // HEPEVT.idhep[i],   // pdg
                HEPEVT.isthep[i],  // status
                HEPEVT.jmohep[i][0]-1,   // FORTRAN to C index shift
                HEPEVT.jmohep[i][1]-1,   // FORTRAN to C index shift
                HEPEVT.jdahep[i][0]-1,   // FORTRAN to C index shift
                HEPEVT.jdahep[i][1]-1,   // FORTRAN to C index shift
                HEPEVT.phep[i][0],
                HEPEVT.phep[i][1],
                HEPEVT.phep[i][2],
                HEPEVT.phep[i][3],
                HEPEVT.vhep[i][0]*vtxscale,
                HEPEVT.vhep[i][1]*vtxscale,
                HEPEVT.vhep[i][2]*vtxscale,
                HEPEVT.vhep[i][3]*tscale
                );
    TParticle* apart = dynamic_cast<TParticle*>(stdhep[i]);
    apart->SetCalcMass(HEPEVT.phep[i][4]); // set the mass we found
  }

  return numpart;
}

int UtilHepevt::fillCommon(const TClonesArray* stdhep,
                           cvu_t vtxCvu, cvu_t tCvu)
{ 
  // Copy info from TClonesArray of TParticles into /HEPEVT/ common block.

  int numpart = stdhep->GetLast() + 1;  // count of entries in the array
  if (numpart > mxnhep) {
    std::cerr << "UtilHepevt::fillCommon() mxnhep=" << mxnhep 
              << " but StdHep has " << numpart << " entries," << std::endl
              << "  /HEPEVT/ copy will be incomplete." << std::endl;
    numpart = mxnhep;
  }

  double vtxscale = unitScale(vtxCvu);
  double tscale   = unitScale(tCvu);

  HEPEVT.nhep = numpart;

  for (int i=0; i<numpart; ++i) {
    TParticle* apart = dynamic_cast<TParticle*>(stdhep->At(i));
    if (!apart) {
      std::cerr << "UtilHepevt::fillCommon() "
                << "index " << i << " in StdHep is not a TParticle!" << std::endl;
      continue;  // blank entry? 
    }
    HEPEVT.isthep[i]    = apart->GetStatusCode();
    // force PDG # for ions into chosen canonical form
    int ipdg = UtilPDG::stdIonNumber(apart->GetPdgCode());
    HEPEVT.idhep[i]     = ipdg;  // apart->GetPdgCode();
    HEPEVT.jmohep[i][0] = apart->GetMother(0)   + 1; // C to FORTRAN index shift
    HEPEVT.jmohep[i][1] = apart->GetMother(1)   + 1; // C to FORTRAN index shift
    HEPEVT.jdahep[i][0] = apart->GetDaughter(0) + 1; // C to FORTRAN index shift
    HEPEVT.jdahep[i][1] = apart->GetDaughter(1) + 1; // C to FORTRAN index shift
    HEPEVT.phep[i][0]   = apart->Px();
    HEPEVT.phep[i][1]   = apart->Py();
    HEPEVT.phep[i][2]   = apart->Pz();
    HEPEVT.phep[i][3]   = apart->Energy();
    HEPEVT.phep[i][4]   = apart->GetCalcMass(); // is this the right one?
    HEPEVT.vhep[i][0]   = apart->Vx() * vtxscale;
    HEPEVT.vhep[i][1]   = apart->Vy() * vtxscale;
    HEPEVT.vhep[i][2]   = apart->Vz() * vtxscale;
    HEPEVT.vhep[i][3]   = apart->T()  * tscale;
  }
  return numpart;
}

int  UtilHepevt::openBinary(std::string fname, bool write)
{
  if (fname == "") {
    std::cout << "UtilHepevt::openBinary can't open nameless file " << std::endl;
    return -1;
  }
  std::cout << "UtilHepevt::openBinary file: " << fname << " for "
            << (write?"write":"read") << "." << std::endl;

  gSystem->ResetErrno();

  int flags =  O_RDWR | O_CREAT | O_TRUNC ;  // | O_BINARY;
  if (!write) flags = O_RDONLY;

  int mode  = S_IREAD | S_IWRITE;   // chmod bits

  int fd = ::open(fname.c_str(),flags,mode);

  if (fd < 0) {
    std::cerr << "Failed to open file: " << fname 
              << " Errno " << gSystem->GetErrno() << std::endl;
  }
  return fd;

}
 
void UtilHepevt::closeBinary(int& fd)
{
  if ( fd < 0 ) return;  // nothing open
  ::close(fd);
  fd = -1;
}

int UtilHepevt::readBinary(int& fd)
{ 
  // Read FORTRAN compatible binary file of /HEPEVT/
  HEPEVT.nhep = 0;
  if ( fd < 0 ) return -1;

  const ssize_t szi4 = sizeof(Int_t);
  const ssize_t szd8 = sizeof(Double_t);

  ssize_t szrecmrk   = szi4;  // initial guess: record marker size is i4

  Int_t    ibufhead[6];     // recsiz [0] nevhep nhep recsiz [0]
  Int_t    ibufrowsz1[2];   // recsiz [0] 
  Int_t    ibufrowdata[6];  // ist id jmo1 jmo2 jda1 jda2
  Double_t dbufrowdata[9];  // px py pz e m vx vy vz t
  Int_t    ibufrowsz2[2];   // recsiz [0] 

  ssize_t rhead = ::read(fd,ibufhead,2*szi4+2*szrecmrk); // assume 4byte record markers
  // std::cout << "readBinary rhead=" << rhead << " errno " << gSystem->GetErrno() << std::endl;
  if (rhead == 0) {
    // saw EOF
    closeBinary(fd);
    return 0;
  }
  if (ibufhead[0] != 2*sizeof(Int_t)) {  // first record should just be two ints nevhep,nhep
    std::cerr << "UtilHepevt::readBinary on fd=" << fd
              << " failed first recsize test " << ibufhead[0] << std::endl;
    return -2;
  }
  if (ibufhead[0] == ibufhead[3]) {
    HEPEVT.nevhep = ibufhead[1];
    HEPEVT.nhep   = ibufhead[2];
  } else {
    // oops .. 8bytes record markers
    szrecmrk = 2*szi4;
    rhead += ::read(fd,ibufhead+4*szi4,2*szi4);  // read two more ints to flesh out the szrecmrk=8
    if (ibufhead[0] != ibufhead[5]) {
      std::cerr << "UtilHepevt::readBinary on fd=" << fd << " tried recmrksiz=8, "
                << " failed 2nd recsize test " << ibufhead[0] 
                << " " << ibufhead[5] << std::endl;
      return -3;
    }
    HEPEVT.nevhep = ibufhead[2];
    HEPEVT.nhep   = ibufhead[3];
  }

  // now read data rows
  int numpart   = HEPEVT.nhep;
  for (int i = 0; i < numpart; ++i) {
    ssize_t rsz1  = ::read(fd,ibufrowsz1,szrecmrk);
    ssize_t rirow = ::read(fd,ibufrowdata,6*szi4);
    ssize_t rdrow = ::read(fd,dbufrowdata,9*szd8);
    ssize_t rsz2  = ::read(fd,ibufrowsz2,szrecmrk);
    if (ibufrowsz1[0] != 6*szi4+9*szd8 || ibufrowsz1[0] != ibufrowsz2[0] ||
        rsz1 != szrecmrk || rirow != 6*szi4 || rdrow != 9*szd8 || rsz2 != szrecmrk) {
      std::cerr << "UtilHepevt::readBinary on fd=" << fd
                << ", szrecmrk=" << szrecmrk << "failed.";
      std::cerr << " rsz1 " << rsz1 << " data: " << ibufrowsz1[0] << " " << ibufrowsz1[1] << std::endl;
      std::cerr << " rsz2 " << rsz2 << " data: " << ibufrowsz2[0] << " " << ibufrowsz2[1] << std::endl;
      std::cerr << " rirow " << rirow << " rdrow " << rdrow << std::endl;
      return -4;
    } // bad read of row due to rec marker problems or size problems
    HEPEVT.isthep[i]     = ibufrowdata[0];
    // force PDG # for ions into chosen canonical form
    int ipdg = UtilPDG::stdIonNumber(ibufrowdata[1]);
    HEPEVT.idhep[i]      = ipdg; // ibufrowdata[1];
    HEPEVT.jmohep[i][0]  = ibufrowdata[2];
    HEPEVT.jmohep[i][1]  = ibufrowdata[3];
    HEPEVT.jdahep[i][0]  = ibufrowdata[4];
    HEPEVT.jdahep[i][1]  = ibufrowdata[5];
    HEPEVT.phep[i][0]    = dbufrowdata[0];
    HEPEVT.phep[i][1]    = dbufrowdata[1];
    HEPEVT.phep[i][2]    = dbufrowdata[2];
    HEPEVT.phep[i][3]    = dbufrowdata[3];
    HEPEVT.phep[i][4]    = dbufrowdata[4];
    HEPEVT.vhep[i][0]    = dbufrowdata[5];
    HEPEVT.vhep[i][1]    = dbufrowdata[6];
    HEPEVT.vhep[i][2]    = dbufrowdata[7];
    HEPEVT.vhep[i][3]    = dbufrowdata[8];
  } // end of rows

  return numpart;  // success
}

bool UtilHepevt::writeBinary(int fd, int szrecmrk)
{ 
  // Write FORTRAN compatible binary file of /HEPEVT/

  const ssize_t szi4 = sizeof(Int_t);
  const ssize_t szd8 = sizeof(Double_t);

  if (szi4 != szrecmrk && 2*szi4 != szrecmrk) {
    std::cerr << "UtilHepevt::writeBinary of fd=" << fd
              << " was passed illegal szrecmrk=" << szrecmrk
              << " must be " << szi4 << " or "
              << 2*szi4 << "."  << std::endl;
    return false;
  }

  Int_t    ibufhead[6]    = { 0, 0, 0, 0, 0, 0 };  // recsiz [0] nevhep nhep recsiz [0]
  Int_t    ibufrowsz[2]   = { 0, 0 };              // recsiz [0] 
  Int_t    ibufrowdata[6] = { 0, 0, 0, 0, 0, 0 };  // ist id jmo1 jmo2 jda1 jda2
  Double_t dbufrowdata[9] = { 0 };                 // px py pz e m vx vy vz t

  int i = 0, sztowrite=2*szi4+2*szrecmrk;
  ssize_t whead, wsz1, wirow, wdrow, wsz2;

  // write the header record
  ibufhead[i] = 2*szi4;  // first record has only two ints
  if (szi4 != szrecmrk) i++;  // skip a word for 8 byte record markers
  ibufhead[i+1] = HEPEVT.nevhep;
  ibufhead[i+2] = HEPEVT.nhep;
  ibufhead[i+3] = 2*szi4;
  whead = ::write(fd,ibufhead,sztowrite);
  if (whead != sztowrite) {
    std::cerr << "UtilHepevt::writeBinary fd=" << fd << " failed first write of " 
              << sztowrite << ", just " << whead << "." << std::endl;
    return false;
  }

  // write individual rows
  const int recsiz = 6*szi4 + 9*szd8;
  int numpart = HEPEVT.nhep;
  for (int i=0; i<numpart; ++i) {
    ibufrowsz[0] = recsiz;

    ibufrowdata[0] = HEPEVT.isthep[i];
    // force PDG # for ions into chosen canonical form
    int ipdg = UtilPDG::stdIonNumber(HEPEVT.idhep[i]);
    ibufrowdata[1] = ipdg; // HEPEVT.idhep[i];
    ibufrowdata[2] = HEPEVT.jmohep[i][0];
    ibufrowdata[3] = HEPEVT.jmohep[i][1];
    ibufrowdata[4] = HEPEVT.jdahep[i][0];
    ibufrowdata[5] = HEPEVT.jdahep[i][1];

    dbufrowdata[0] = HEPEVT.phep[i][0];
    dbufrowdata[1] = HEPEVT.phep[i][1];
    dbufrowdata[2] = HEPEVT.phep[i][2];
    dbufrowdata[3] = HEPEVT.phep[i][3];
    dbufrowdata[4] = HEPEVT.phep[i][4];
    dbufrowdata[5] = HEPEVT.vhep[i][0];
    dbufrowdata[6] = HEPEVT.vhep[i][1];
    dbufrowdata[7] = HEPEVT.vhep[i][2];
    dbufrowdata[8] = HEPEVT.vhep[i][3];

    wsz1  = ::write(fd,ibufrowsz,szrecmrk);
    wirow = ::write(fd,ibufrowdata,6*szi4);
    wdrow = ::write(fd,dbufrowdata,9*szd8);
    wsz2  = ::write(fd,ibufrowsz,szrecmrk);
    //std::cout << " wsz " << wsz1 << " " << wsz2 << " " << szrecmrk << std::endl;
    if (wsz1 != szrecmrk || wsz2 != szrecmrk || wirow != 6*szi4 || wdrow != 9*szd8) {
      std::cerr << "UtilHepevt::writeBinary fd=" << fd << " failed sizerec write 1" << std::endl;
      std::cerr << " wsz1 " << wsz1 << " wsz2 " << wsz2 
                << " wirow " << wirow << " wdrow " << wdrow << "." << std::endl;
      return false;
    }

  } // end of rows
  return true;
}

void UtilHepevt::dump(const TClonesArray* stdhep, int mode)
{ 
  // Formatted dump of StdHep array.

   // mode < 0 --> nothing
   // mode = 0 --> suppress vtx info
   // mode = 1 --> include vtx info

   if (mode < 0 ) return;

   bool dovtx = (mode > 0);
   bool doheader = true;

   // Loop over array ... print each element
   int n = stdhep->GetLast()+1;
   printf("StdHep:  %d entries\n",n);
   for (Int_t i=0; i<n; i++) {
     TParticle* apart = dynamic_cast<TParticle*>(stdhep->At(i));
     if (!apart) {
       std::cout << "dumpStdhep entry " << i << " is empty" << std::endl;
       continue;
     }

     //rwh: apart->Print() -- does a cruddy job, do it ourselves
     
     // header before the first entry
     if (doheader) {
       doheader = kFALSE;
       const char* head1 =
         "ihep stat  type   [parent][chldrn]         px          py          pz     energy\n";
       // 12345678901234567890123456789012345678901234567890123456789012345678901234567890
       //                          1234123412341234 123456789A 123456789A  123456789A 12345678
       const char* head2 =
         "                                           vx          vy          vz      tprod\n";
       printf(head1);
       if (dovtx) printf(head2);
     }
     
     static Char_t alt[64];
     
     Int_t status = apart->GetStatusCode();
     // if there isn't a name construct one
     const Char_t* name = apart->GetName();
     if (name[0]=='X' && name[1]=='X' && name[2] == 'X') {
       sprintf(alt,"%d",apart->GetPdgCode());
       name = alt;
     }
     Int_t first_parent = apart->GetMother(0);
     //Int_t last_parent  = apart->GetMother(1);
     Int_t first_child  = apart->GetDaughter(0);
     Int_t last_child   = apart->GetDaughter(1);
     Double_t px = apart->Px();
     Double_t py = apart->Py();
     Double_t pz = apart->Pz();
     Double_t energy = apart->Energy();
     //Double_t calcmass = apart->GetCalcMass();
     //Double_t mass = (apart->GetPDG())?apart->GetMass():-99999.;
     Double_t vx = apart->Vx();
     Double_t vy = apart->Vy();
     Double_t vz = apart->Vz();
     Double_t tprod = apart->T();
     
     printf("%4d(%4d)%-10.10s %4d%4d%4d",
            i,status,name,first_parent,first_child,last_child);
     printf(" %11.4e %11.4e %11.4e%11.4e\n",px,py,pz,energy);
     if (dovtx) {
       printf("                                 ");
       printf(" %11.4e %11.4e %11.4e%11.4e\n",vx,vy,vz,tprod);
     }
   }  // loop over entries

}

double UtilHepevt::unitScale(UtilHepevt::cvu_t cvu)
{
  switch (cvu) {
  case k_none:     return    1.0;
  case k_mmtosec:  return    0.001 / TMath::C();
  case k_sectomm:  return 1000.0   * TMath::C();
  case k_cmtom:    return    0.01;
  case k_mtocm:    return  100.0;
  case k_mmtom:    return    0.001;
  case k_mtomm:    return 1000.0;
  default:
      std::cerr << "UtilHepevt::unitScale() passed invalid enum "
                << (int)cvu << "." << std::endl;
      return 1.0;
  }
}

int UtilHepevt::modStatusStdHep(TClonesArray* stdhep, std::string conditional,
                                 int istflag)
{

  TApplication* app = gROOT->GetApplication();
  std::string tmp;
  int nflagged = 0;
  int nhep = stdhep->GetLast()+1;
  for (int i=0; i<nhep; ++i) {
    TParticle* apart = dynamic_cast<TParticle*>(stdhep->At(i));
    int ipdg = apart->GetPdgCode();
    bool ischarm = UtilPDG::isCharm(ipdg);
    bool istau   = UtilPDG::isTau(ipdg);
    int  mom0    = apart->GetMother(0);
    int  mom1    = apart->GetMother(1);
    bool fromcharm = false;
    bool fromtau   = false;
    if (mom0 >= 0 && mom1 >= mom0 ) {
      for (int k=mom0; k<=mom1; ++k) {
        TParticle* mompart = dynamic_cast<TParticle*>(stdhep->At(k));
        int ipdgmom = mompart->GetPdgCode();
        if (UtilPDG::isCharm(ipdgmom)) fromcharm = true;
        if (UtilPDG::isTau(ipdgmom)  ) fromtau   = true;
      }
    }
    std::string fmt = "";
    fmt += " const TParticle* p = (const TParticle*)0x%lx;";
    fmt += " int indx=%d, ist=%d, ipdg=%d, mo0=%d, mo1=%d, da0=%d, da1=%d;";
    fmt += " bool ischarm=%s, istau=%s, fromcharm=%s, fromtau=%s;";
    fmt += " return %s; ";
    tmp = Form(fmt.c_str(),
               reinterpret_cast<unsigned long>(apart),
               i,
               apart->GetStatusCode(),
               ipdg,
               mom0,
               mom1,
               apart->GetDaughter(0),
               apart->GetDaughter(1),
               (ischarm?"true":"false"),
               (istau?"true":"false"),
               (fromcharm?"true":"false"),
               (fromtau?"true":"false"),
               conditional.c_str());
    int pval = app->ProcessLine(tmp.c_str());
    //std::cout << "ProcessLine: " << tmp << " returned " << pval << std::endl;
    if (pval) { apart->SetStatusCode(istflag); ++nflagged; }
  }
  return nflagged;
}

int UtilHepevt::squeezeStdHep(TClonesArray* stdhep, int istrm)
{

  int nhepold = stdhep->GetLast() + 1;  // count of entries
  int nhepnew = 0;
  std::vector<int> newindx(nhepold), ist(nhepold);
  std::vector<int> mo0(nhepold), mo1(nhepold), da0(nhepold), da1(nhepold);
  int i;
  int nremoved = 0;
  // first pass through decides who to keep and puts them in temporary
  for (i=0; i<nhepold; ++i) {
    TParticle* apart = dynamic_cast<TParticle*>(stdhep->At(i));
    ist[i]     = apart->GetStatusCode();
    mo0[i]     = apart->GetMother(0);
    mo1[i]     = apart->GetMother(1);
    da0[i]     = apart->GetDaughter(0);
    da1[i]     = apart->GetDaughter(1);
    newindx[i] = -1;
    if ( ist[i] != istrm ) { newindx[i] = nhepnew; ++nhepnew; }
  }
  
  // determine new mother/daughter indices
  for (i=nhepold-1; i>=0; --i) {
    if (newindx[i] < 0) {
      mo0[i] = -1;
      mo1[i] = -1;
      da0[i] = -1;
      da1[i] = -1;
      continue;  // not keeping it ...
    } else {
      if (mo0[i]>=0) mo0[i] = newindx[mo0[i]];
      if (mo1[i]>=0) mo1[i] = newindx[mo1[i]];
      if (da0[i]>=0) da0[i] = newindx[da0[i]];
      if (da1[i]>=0) da1[i] = newindx[da1[i]];
    }
    
  }

  // second pass through asks that the entry be removed, 
  // or fixes the mother/dauther indices
  for (i=nhepold-1; i>=0; --i) { 
    if ( newindx[i] < 0 ) {
      stdhep->RemoveAt(i); 
      ++nremoved;
    } else {
      TParticle* apart = dynamic_cast<TParticle*>(stdhep->At(i));
      apart->SetFirstMother(mo0[i]);
      apart->SetLastMother(mo1[i]);
      apart->SetFirstDaughter(da0[i]);
      apart->SetLastDaughter(da1[i]);
    }
  }
  if (nremoved>0) stdhep->Compress();

  return nremoved;
}

TClonesArray* UtilHepevt::readAscii(std::istream& istrm, 
                                    const TLorentzVector vtx_offset)
{
  // read event record from ascii stream 
  // return null pointer if EOF or error

  const int LINEMX = 1024;
  char buffer[ LINEMX ];
  std::string line;

  enum Efilefmt { kUnknown, kHeplst4, kExodus, kdump, kArgoneut };
  const char* fmtname[] = 
    { "unknown", "Heplst4", "Exodus", "dump", "Argoneut" };

  bool isbad = false;
  Efilefmt thisfilefmt = kUnknown;

  if ( ! istrm.good() ) {
    // bad input return null pointer
    std::cerr << "UtilHepevt::readAscii input stream !good()" << std::endl;
    return 0;
  }

  int nhep = 0;
  while ( thisfilefmt == kUnknown ) {
    // for now only handle file layouts that explicitly give row counts
    istrm.getline(buffer,LINEMX);
    //std::cout << "HEY: read line \"" << buffer << "\"" << std::endl;
    if ( ! istrm.good() ) return 0; 
    line = buffer;
    if ( line.find("StdHep:") != std::string::npos ) {
      std::cout << buffer << std::endl;
      thisfilefmt = kdump;
      if ( line.find("entries") == std::string::npos ) 
        thisfilefmt = kArgoneut;
      sscanf(buffer,"StdHep:%d",&nhep);
      istrm.getline(buffer,LINEMX);  // header line
      istrm.getline(buffer,LINEMX);  // header line
    } else if ( line.find("nhep =") != std::string::npos ) {
        std::cout << "** Saw nhep=" << std::endl;
        std::cout << buffer << std::endl;
        thisfilefmt = kExodus;
        sscanf(buffer,"nhep =%d",&nhep);
        istrm.getline(buffer,LINEMX);  // header line
        istrm.getline(buffer,LINEMX);  // header line
    }
  }
  //std::cout << "HEY: saw file type " << (int)thisfilefmt 
  //          << " " << fmtname[thisfilefmt] << " nhep = " 
  //          << nhep << std::endl;

  // get an empty TClonesArray of TParticles, we know how many entries
  TClonesArray* stdhep_ptr =  getEmptyStdHep(nhep);

  // no-entries return an empty list
  if ( nhep <= 0 ) return stdhep_ptr;

  TClonesArray &stdhep = *stdhep_ptr;

  double vtxscale = 1.0, tscale = 1.0;
  int indx_offset = 0;  // mother+daughter index ajustment FORTRAN vs. C
  switch ( thisfilefmt ) {
  case kHeplst4:  indx_offset = 1; vtxscale = 1.000; tscale = 1.000;  break; 
  default:        indx_offset = 0; 
  }

  // read the entries for this record
  for (int irow = 0; irow < nhep; ++irow) {
    // read first of pair
    char pname[ LINEMX ];
    int irowx, nitems;
    int isthep = -999, ipdg = -999;
    int jmo1 = -999, jmo2 = -999, jda1 = -999, jda2 = -999;
    double phep[5] = { 0, 0, 0, 0, 0 };
    double vhep[4] = { 0, 0, 0, 0 };


    istrm.getline(buffer,LINEMX);  // status, particle, energy line
    if ( ! istrm.good() ) {
      std::cerr << " irow " << irow << " != good()" << std::endl;
      isbad = true;
      break;
    }
    // parse the line
    if ( thisfilefmt == kdump ) {
      nitems = sscanf(buffer,"%d(%d)%s %d %d %d %lf %lf %lf %lf",
                      &irowx,&isthep,pname,&jmo1,&jda1,&jda2,
                      &phep[0],&phep[1],&phep[2],&phep[3]);
      jmo2 = jmo1;
      // translate pname to ipdg here ...
      std::cerr << "pname: " << pname << std::endl;
    } else if ( thisfilefmt == kExodus ) {
      nitems = sscanf(buffer,"%d(%d)%s<%d %d> %lf %lf %lf %lf",
                      &irowx,&isthep,pname,&jmo1,&jmo2,
                      &phep[0],&phep[1],&phep[2],&phep[3]);
      if (nitems != 9) {
        std::cerr << "irow " << irow << " read only " << nitems 
                  << " of 9 elements" << std::endl;
        isbad = true;
        break;
      }
    } else if ( thisfilefmt == kArgoneut ) {
      nitems = sscanf(buffer,"%d(%d)%d %d %d %d %lf %lf %lf %lf",
                      &irowx,&isthep,&ipdg,&jmo1,&jda1,&jda2,
                      &phep[0],&phep[1],&phep[2],&phep[3]);
      jmo2 = jmo1;
    }
    // read second of pair
    istrm.getline(buffer,LINEMX);  // vertex info
    if ( thisfilefmt == kdump ) {
    } else if ( thisfilefmt == kExodus ) {
    } else if ( thisfilefmt == kArgoneut ) {
      nitems = sscanf(buffer,"%lf %lf %lf %lf",
                      &vhep[0],&vhep[1],&vhep[2],&vhep[3]);
    }

    if (irowx != irow ) {
      std::cerr << " irow " << irow << " != irowx " << irowx << std::endl;
      isbad = true;
      break;
    }

    if ( ipdg > UtilPDG::kIonBase )
      ipdg = UtilPDG::stdIonNumber(ipdg);  // cannonical form
    new(stdhep[irow]) 
      TParticle(ipdg, // HEPEVT.idhep[i],   // pdg
                isthep,  // status
                jmo1 - indx_offset, // FORTRAN to C index shift
                jmo2 - indx_offset, // FORTRAN to C index shift
                jda1 - indx_offset, // FORTRAN to C index shift
                jda2 - indx_offset, // FORTRAN to C index shift
                phep[0],
                phep[1],
                phep[2],
                phep[3],
                vhep[0]*vtxscale + vtx_offset[0],
                vhep[1]*vtxscale + vtx_offset[1],
                vhep[2]*vtxscale + vtx_offset[2],
                vhep[3]*tscale   + vtx_offset[3]
                );

  }

  return stdhep_ptr;
}

void UtilHepevt::AppendStdHep(TClonesArray* stdhep_ptr, const TClonesArray* addon,
                              const TLorentzVector vtx_offset)
{

  if (stdhep_ptr == 0) return;  // user failed to supply an array
  TClonesArray &stdhep = *stdhep_ptr;

  int numpart = stdhep_ptr->GetLast() + 1; // count of entries in original list
  int nadd    = addon->GetLast() + 1;  // count of entries to add

  for (int i = 0; i<nadd; ++i) {
    TParticle* apart = dynamic_cast<TParticle*>(addon->At(i));
    if (!apart) {
      std::cerr << "UtilHepevt::AppendStdHep() "
                << "index " << i << " in StdHep is not a TParticle!" << std::endl;
      continue;  // blank entry? 
    }
    Int_t indx = numpart + i;
    // extract the information from the entry in the old list
    Int_t isthep    = apart->GetStatusCode();
    // force PDG # for ions into chosen canonical form
    Int_t ipdg = UtilPDG::stdIonNumber(apart->GetPdgCode());
    Int_t jmohep[2], jdahep[2];
    jmohep[0] = apart->GetMother(0);
    jmohep[1] = apart->GetMother(1);
    jdahep[0] = apart->GetDaughter(0);
    jdahep[1] = apart->GetDaughter(1);
    if ( jmohep[0] >= 0 ) jmohep[0] += numpart; // shift index
    if ( jmohep[1] >= 0 ) jmohep[1] += numpart; // shift index
    if ( jdahep[0] >= 0 ) jdahep[0] += numpart; // shift index
    if ( jdahep[1] >= 0 ) jdahep[1] += numpart; // shift index
    Double_t phep[5], vhep[4];
    phep[0]   = apart->Px();
    phep[1]   = apart->Py();
    phep[2]   = apart->Pz();
    phep[3]   = apart->Energy();
    phep[4]   = apart->GetCalcMass(); // is this the right one?
    vhep[0]   = apart->Vx() + vtx_offset.X();
    vhep[1]   = apart->Vy() + vtx_offset.Y();
    vhep[2]   = apart->Vz() + vtx_offset.Z();
    vhep[3]   = apart->T()  + vtx_offset.T();
    // create the new entry
    new(stdhep[indx]) TParticle(ipdg,    // pdg
                                isthep,  // status
                                jmohep[0], jmohep[1],
                                jdahep[0], jdahep[1],
                                phep[0], phep[1], phep[2], phep[3],
                                vhep[0], vhep[1], vhep[2], vhep[3]
                                );
    TParticle* newpart = dynamic_cast<TParticle*>(stdhep[indx]);
    newpart->SetCalcMass(phep[4]); // set the mass we found
  }

}

---