---

HepEvt.cxx

Go to the documentation of this file.

// -*- C++ -*-  $Id: HepEvt.cxx,v 1.1.1.1 1999/12/20 22:28:41 burnett Exp $
//
// This file is part of Gismo 2
// inmplement HepEvt class
#include "gismo/HepEvt.h"
#include "facilities/error.h"
#include "gismo/Units.h"
#include "gismo/PDT.h"

#include <vector>
#include <string.h>
#include <stdio.h>
#include <assert.h>

static char line[133];


HepEvt::HepEvt(std::istream & in)
: Generator() // make a default something
, inputFile(in)
{
        if( !in || in.bad() )
                 FATAL("HepEvt: input file not open!");
        _hasVertices=false;
        for(;;)
        {
                in.getline(line,sizeof(line)-1);
                if(in.eof()){ nevhep()=-1; break;}
                if( strstr(line,"Event listing (HEP format"))
                {
                        nevhep()=0;
                        if( strstr(line, " with vertices"))
                           _hasVertices = true;
                        break;
                }
        }
        // readEvent();
        if( in.eof() )
                WARNING("HepEvt: no event data found!");

}

/*   start of record to decode:
                                                Event listing (HEP format with vertices)            Event:        1
         I  particle/jet     ISTHEP  IDHEP     JMOHEP      JDAHEP    PHEP(1,I)   PHEP(2,I)   PHEP(3,I)   PHEP(4,I)   PHEP(5,I)
                                                                                                                                                                          VHEP(1,I)   VHEP(2,I)   VHEP(3,I)   VHEP(4,I)

         1  (upsilon(4S))         2  70553     0     0     2     3      .00000      .00000     5.90889    12.11822    10.58000
                                                                                                                                                                                .000        .000        .000        .000
*/
int
 HepEvt::readEvent()
{
        char line[140];
        nhep() = 0;
        int i;
        for(;;)  // read forward until find line starting with  1
        {
                inputFile.getline(line,sizeof(line)-1);
                if( inputFile.eof() )
                        return 1;
                if( strlen(line)==0 ) continue;
                if( (i =atoi(line)) ==1 ) break;
        }
        // found event: now fill the Fortran common block
        char name[20];
        for(;atoi(line)>0;i++)
        {  int index;

                sscanf(line,"%d %s %D %D %D %D %D %D %g %g %g %g %g",
                                &index, &name, &isthep(i), &idhep(i),
                                jmohep(i),jmohep(i)+1, jdahep(i), jdahep(i)+1,
                                phep(i), phep(i)+1, phep(i)+2, phep(i)+3,  phep(i)+4);
                if( index!=i )
                        FATAL("Bad input data for HepEvt ");

                if( _hasVertices){
                  inputFile.getline(line,sizeof(line)-1);
                  sscanf(line,"%g %g %g %g", vhep(i),vhep(i)+1, vhep(i)+2, vhep(i)+3);
                }else{
                  vhep(i)[0]=vhep(i)[1]=vhep(i)[2]=vhep(i)[3];
                }

                inputFile.getline(line,sizeof(line)-1);

                // add the particle to the PDT, if not there
                PData* data = thePDT()->lookUp(idhep(i));
                if( !data )
                {
                  float mass = phep(i)[4];
                  char* pname=name;
                  if( name[0]=='(' && strlen(name)>0)
                  {
                          name[strlen(name)-1]='\0';
                          pname++;
                  }
                  thePDT()->addParticle(pname, idhep(i), 99, 99, mass );
                }
        }
        nhep() = i;
        nevhep()++;
        return (i==0); //nonzero means failure
}

static MCParticle* primary;
static void setupParticle(MCParticle& q, float* p, float * v)
{
  q.setMomentum(Vector(p[0],p[1],p[2]));
  double conv = Units::from_mm(1);
  q.setPosition( Point(conv*v[0],conv*v[1],conv*v[2]) );
  q.setT(Units::from_mm(v[3]));
}
MCParticle *
HepEvt::getEvent()
{
  removeChildren();

  //HepAList<MCParticle>
  std::vector<MCParticle*>  plist(20);  // corresponding list of created particles
  if( readEvent() )
         return 0;
  long primaryid = idhep(1);
  pData = thePDT()->lookUp(primaryid);
  setMass( phep(1)[4]);
  //THB: why?? _parent = 0;
  setupParticle(*this, phep(1),vhep(1));
  setStatus(MCParticle::DECAYED);
  plist.push_back(this);

  unsigned numsecond = unsigned(nhep());
  for( unsigned i=2; i<=numsecond; i++)
  {
          int statusCode = isthep(i);
//        if( statusCode !=1)
//                      continue;
          if( statusCode==0 ) continue;
          PData* data = thePDT()->lookUp(idhep(i));
          float mass = phep(i)[4];
          int parentIndex = static_cast<int>(jmohep(i)[0]-1);
          MCParticle* child = dynamic_cast<MCParticle*>
                          (plist[parentIndex]->addChild( data ,mass) );
          setupParticle(*child, phep(i),vhep(i));
          if( statusCode==2)
                  child->setStatus(MCParticle::DECAYED);
          plist.push_back(child);
  }
  return this;
}

void
HepEvt::propagate(const Medium * med)
{
        MCParticle::propagate(med);
}

---