HepMC.tar

HepMC/000755 000765 000024 00000000000 13310701701 012013 5ustar00sa639staff000000 000000 HepMC/PdfInfo.h000644 000765 000024 00000015500 13310701701 013512 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HEPMC_PDF_INFO_H
#define HEPMC_PDF_INFO_H

//////////////////////////////////////////////////////////////////////////
// garren@fnal.gov, July 2006
//
// Additional PDF information
//////////////////////////////////////////////////////////////////////////

namespace HepMC {

//! The PdfInfo class stores PDF information

///
/// \class  PdfInfo
/// HepMC::PdfInfo stores additional PDF information for a GenEvent.
/// Creation and use of this information is optional. 
///
///   - int id1; 	// flavour code of first parton
///   - int id2; 	// flavour code of second parton 
///   - int pdf_id1; 	// LHAPDF set id of first parton (zero by default)
///   - int pdf_id2; 	// LHAPDF set id of second parton  (zero by default)
///   - double x1;  	// fraction of beam momentum carried by first parton ("beam side") 
///   - double x2;  	// fraction of beam momentum carried by second parton ("target side") 
///   - double scalePDF;  //  Q-scale used in evaluation of PDF's   (in GeV) 
///   - double pdf1;  	// PDF (id1, x1, Q) 
///   - double pdf2;  	// PDF (id2, x2, Q)   
/// 
/// Input parton flavour codes id1 & id2 are expected to obey the 
/// PDG code conventions, especially g = 21.
/// 
/// The contents of pdf1 and pdf2 are expected to be x*f(x).
/// The LHAPDF set ids are the entries in the first column of
/// http:///projects.hepforge.org/lhapdf/PDFsets.index
///
class PdfInfo {

public:
  // ---  birth/death:
  //
  /// default constructor
  PdfInfo()
    : m_id1(0), 
      m_id2(0),
      m_pdf_id1(0),
      m_pdf_id2(0),
      m_x1(0),
      m_x2(0),
      m_scalePDF(0),
      m_pdf1(0),
      m_pdf2(0)
    {}

  /// all values EXCEPT pdf_id1 and pdf_id2 must be provided
  PdfInfo( int i1, int i2, double x1, double x2,
           double q, double p1, double p2,
	   int pdf_id1 = 0, int pdf_id2 = 0 );

  ~PdfInfo() {}
    
  // ---  copying:
  //
  PdfInfo( PdfInfo const & orig );	//!< copy constructor
  PdfInfo &  operator = ( PdfInfo const & rhs ); //!< make a copy
  void swap( PdfInfo & other );	//!< swap two PdfInfo objects

  // ---  equivalence:
  //
  bool    operator==( const PdfInfo& ) const; //!< check for equality
  bool    operator!=( const PdfInfo& ) const; //!< check for inequality

  // ---  accessors:
    /// flavour code of first parton
    int   id1()                   const { return m_id1; }
    /// flavour code of second parton 
    int   id2()                   const { return m_id2; }
    /// LHAPDF set id of first parton
    int   pdf_id1()                   const { return m_pdf_id1; }
    /// LHAPDF set id of second parton
    int   pdf_id2()                   const { return m_pdf_id2; }
    /// fraction of beam momentum carried by first parton ("beam side") 
    double x1()                   const { return m_x1; } 
    /// fraction of beam momentum carried by second parton ("target side") 
    double x2()                   const { return m_x2; }
    ///  Q-scale used in evaluation of PDF's   (in GeV) 
    double scalePDF()             const { return m_scalePDF;  }
    /// PDF (id1, x1, Q) - x*f(x)
    double pdf1()                 const { return m_pdf1; }
    /// PDF (id2, x2, Q) - x*f(x)
    double pdf2()                 const { return m_pdf2; }

    /// verify that the instance contains non-zero information
    bool  is_valid()                     const;

  // ---  mutators:
    /// set flavour code of first parton
    void   set_id1(const int &i)              { m_id1=i; }
    /// set flavour code of second parton 
    void   set_id2(const int &i)              { m_id2=i; }
    /// set LHAPDF set id of first parton
    void   set_pdf_id1(const int &i)              { m_pdf_id1=i; }
    /// set LHAPDF set id of second parton 
    void   set_pdf_id2(const int &i)              { m_pdf_id2=i; }
    /// set fraction of beam momentum carried by first parton ("beam side") 
    void   set_x1(const double &f)            { m_x1=f; }
    /// set fraction of beam momentum carried by second parton ("target side") 
    void   set_x2(const double &f)            { m_x2=f; }
    /// set Q-scale used in evaluation of PDF's   (in GeV) 
    void   set_scalePDF(const double &f)      { m_scalePDF=f;  }
    /// set x*f(x) of first parton
    void   set_pdf1(const double &f)          { m_pdf1=f; }
    /// set x*f(x) of second parton
    void   set_pdf2(const double &f)          { m_pdf2=f; }

private: // data members
    int   m_id1; 
    int   m_id2; 
    int   m_pdf_id1; 
    int   m_pdf_id2; 
    double m_x1;
    double m_x2;
    double m_scalePDF; 
    double m_pdf1;
    double m_pdf2;

};

// Free Functions

// IO
std::ostream & operator << (std::ostream &, PdfInfo const *);
std::istream & operator >> (std::istream &, PdfInfo *);

// inline operators
inline PdfInfo::PdfInfo( int i1, int i2, double xin1, double xin2,
        		 double q, double p1, double p2,
			 int pid1, int pid2 )
    : m_id1(i1), 
      m_id2(i2),
      m_pdf_id1(pid1),
      m_pdf_id2(pid2),
      m_x1(xin1),
      m_x2(xin2),
      m_scalePDF(q),
      m_pdf1(p1),
      m_pdf2(p2)
   {}

inline PdfInfo::PdfInfo( PdfInfo const & orig )
    : m_id1(orig.m_id1), 
      m_id2(orig.m_id2),
      m_pdf_id1(orig.m_pdf_id1), 
      m_pdf_id2(orig.m_pdf_id2),
      m_x1(orig.m_x1),
      m_x2(orig.m_x2),
      m_scalePDF(orig.m_scalePDF),
      m_pdf1(orig.m_pdf1),
      m_pdf2(orig.m_pdf2)
   {}

inline PdfInfo &  PdfInfo::operator = ( PdfInfo const & rhs ) 
{
  PdfInfo temp( rhs );
  swap( temp );
  return *this;
}

inline void PdfInfo::swap( PdfInfo & other ) 
{
  std::swap(m_id1, other.m_id1); 
  std::swap(m_id2, other.m_id2);
  std::swap(m_pdf_id1, other.m_pdf_id1); 
  std::swap(m_pdf_id2, other.m_pdf_id2);
  std::swap(m_x1, other.m_x1);
  std::swap(m_x2, other.m_x2);
  std::swap(m_scalePDF, other.m_scalePDF);
  std::swap(m_pdf1, other.m_pdf1);
  std::swap(m_pdf2, other.m_pdf2);
}

inline bool    PdfInfo::operator==( const PdfInfo& a ) const
{
    /// equality requires that each member match
    return (    a.id1() == this->id1() 
             && a.id2() == this->id2() 
             && a.pdf_id1() == this->pdf_id1() 
             && a.pdf_id2() == this->pdf_id2() 
	     && a.x1() == this->x1() 
	     && a.x2() == this->x2()
	     && a.scalePDF() == this->scalePDF() 
	     && a.pdf1() == this->pdf1()
	     && a.pdf2() == this->pdf2() );
}

inline bool    PdfInfo::operator!=( const PdfInfo& a ) const
{
    /// any nonmatching member generates inequality
    return !( a == *this );
}

inline bool  PdfInfo::is_valid() const
{
    if( m_id1 != 0 ) return true;
    if( m_id2 != 0 ) return true;
    if( m_pdf_id1 != 0 ) return true;
    if( m_pdf_id2 != 0 ) return true;
    if( m_x1 != 0 ) return true;
    if( m_x2 != 0 ) return true;
    if( m_scalePDF != 0 ) return true;
    if( m_pdf1 != 0 ) return true;
    if( m_pdf2 != 0 ) return true;
    return false;
}

} // HepMC

#endif  // HEPMC_PDF_INFO_H
HepMC/SimpleVector.icc000644 000765 000024 00000007704 13310701701 015117 0ustar00sa639staff000000 000000 //////////////////////////////////////////////////////////////////////////
// SimpleVector.icc
//////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////
// garren@fnal.gov, July 2006
//
//
//////////////////////////////////////////////////////////////////////////

#include <cmath>
#include <algorithm>	// for swap

namespace HepMC {

//////////////////////////////////////////////////////////////////////////
//  FourVector inline methods
//////////////////////////////////////////////////////////////////////////

inline void FourVector::swap( FourVector & other ) {
    std::swap( m_x, other.m_x );
    std::swap( m_y, other.m_y );
    std::swap( m_z, other.m_z );
    std::swap( m_t, other.m_t );
}

inline FourVector & FourVector::operator=(const FourVector & v) {
  m_x = v.x();
  m_y = v.y();
  m_z = v.z();
  m_t = v.t();
  return *this;
}

inline void FourVector::set(double xin, double yin, double zin, double  tin) {
  m_x = xin;
  m_y = yin;
  m_z = zin;
  m_t = tin;
}

inline double FourVector::m2() const {
  return m_t*m_t - (m_x*m_x + m_y*m_y + m_z*m_z);
}

inline double FourVector::m() const {
  double mm = m2();
  return mm < 0.0 ? -std::sqrt(-mm) : std::sqrt(mm);
}

inline double FourVector::perp2() const { return m_x*m_x + m_y*m_y; }

inline double FourVector::perp() const { return std::sqrt(perp2()); }

inline double FourVector::theta() const	 {
  return m_x == 0.0 && m_y == 0.0 && m_z == 0.0 ? 0.0 : std::atan2(perp(),m_z);
}

inline double FourVector::phi() const {
  return m_x == 0.0 && m_y == 0.0 ? 0.0 : std::atan2(m_y,m_x);
}

inline double  FourVector::rho() const { 
return std::sqrt( m_x*m_x + m_y*m_y + m_z*m_z );
}

inline bool FourVector::operator == (const FourVector & v) const {
  return (v.x()==x() && v.y()==y() && v.z()==z() && v.t()==t()) ? true : false;
}

inline bool FourVector::operator != (const FourVector & v) const {
  return (v.x()!=x() || v.y()!=y() || v.z()!=z() || v.t()!=t()) ? true : false;
}

inline double FourVector::pseudoRapidity() const {
  double m1 = std::sqrt( m_x*m_x + m_y*m_y + m_z*m_z );
  if ( m1==  0   ) return  0.0;   
  if ( m1==  z() ) return  1.0E72;
  if ( m1== -z() ) return -1.0E72;
  return 0.5*log( (m1+z())/(m1-z()) );
}

inline double FourVector::eta()    const { return pseudoRapidity();}


//////////////////////////////////////////////////////////////////////////
//  ThreeVector inline methods
//////////////////////////////////////////////////////////////////////////

inline void ThreeVector::swap( ThreeVector & other ) {
    std::swap( m_x, other.m_x );
    std::swap( m_y, other.m_y );
    std::swap( m_z, other.m_z );
}

inline double ThreeVector::theta() const	 {
  return m_x == 0.0 && m_y == 0.0 && m_z == 0.0 ? 0.0 : std::atan2(perp(),m_z);
}

inline double ThreeVector::phi() const {
  return m_x == 0.0 && m_y == 0.0 ? 0.0 : std::atan2(m_y,m_x);
}

inline double ThreeVector::r()    const { 
return std::sqrt( m_x*m_x + m_y*m_y + m_z*m_z );
}

inline void ThreeVector::set(double xin, double yin, double zin) {
  m_x = xin;
  m_y = yin;
  m_z = zin;
}

inline void ThreeVector::setPhi(double ph) { 
  double xy   = perp();
  setX(xy*std::cos(ph));
  setY(xy*std::sin(ph));
}

inline void ThreeVector::setTheta(double th) { 
  double ma   = r();
  double ph   = phi();
  setX(ma*std::sin(th)*std::cos(ph));
  setY(ma*std::sin(th)*std::sin(ph));
  setZ(ma*std::cos(th));
}

inline double ThreeVector::perp2() const { return m_x*m_x + m_y*m_y; }

inline double ThreeVector::perp() const { return std::sqrt(perp2()); }

inline ThreeVector & ThreeVector::operator = (const ThreeVector & p) {
  m_x = p.x();
  m_y = p.y();
  m_z = p.z();
  return *this;
}


inline bool ThreeVector::operator == (const ThreeVector& v) const {
  return (v.x()==x() && v.y()==y() && v.z()==z()) ? true : false;
}

inline bool ThreeVector::operator != (const ThreeVector& v) const {
  return (v.x()!=x() || v.y()!=y() || v.z()!=z()) ? true : false;
}

} // HepMC
HepMC/enable_if.h000644 000765 000024 00000002047 13310701701 014073 0ustar00sa639staff000000 000000 #ifndef ENABLE_IF_H
#define ENABLE_IF_H
// author: Walter Brown

// ----------------------------------------------------------------------
// prolog

namespace HepMC  {
namespace detail {


// ----------------------------------------------------------------------
// enable_if<>

/// internal - used to decide if a class is arithmetic
template< bool, class >
struct enable_if
{ };

/// internal - use if class T is arithmetic
template< class T >
struct enable_if<true, T>
{
  typedef  T  type;	//!< check type of class T
};


// ----------------------------------------------------------------------
// disable_if<>

/// internal - used by SimpleVector to decide if a class is arithmetic
template< bool, class >
struct disable_if
{ };

/// internal - used by SimpleVector to decide if a class is arithmetic
template< class T >
struct disable_if<false, T>
{
  typedef  T  type;	//!< check type of class T
};


// ----------------------------------------------------------------------
// epilog

}  // namespace detail
}  // namespace HepMC

#endif  // ENABLE_IF_H
HepMC/Version.h000644 000765 000024 00000001661 13310701701 013615 0ustar00sa639staff000000 000000 #ifndef HEPMC_VERSION_H
#define HEPMC_VERSION_H
// ----------------------------------------------------------------------
//
// Version.h
// Author:  Lynn Garren
//
//  for now, these are free functions
//
// ----------------------------------------------------------------------

#include <string>
#include <iostream>
#include "HepMC/HepMCDefs.h"

namespace HepMC {

void version( std::ostream & os = std::cout );			//!< print HepMC version
void writeVersion( std::ostream & os );	//!< write HepMC version to os
std::string versionName( );	//!< return HepMC version

inline std::string versionName( )
{
    return HEPMC_VERSION;
}

inline void version( std::ostream & os )
{
    os << " --------------- HepMC Version " << versionName()
       << " --------------- " << std::endl;
}

inline void writeVersion( std::ostream & os )
{
    os << "             HepMC Version: " << versionName() << std::endl;
}

}	// HepMC

#endif // HEPMC_VERSION_H
HepMC/IO_Exception.h000644 000765 000024 00000002401 13310701701 014506 0ustar00sa639staff000000 000000 #ifndef HEPMC_IO_EXCEPTION_H
#define HEPMC_IO_EXCEPTION_H
// ----------------------------------------------------------------------
//
// IO_Exception.h
// Author:  Lynn Garren
//
// IO exception handling
// IO_GenEvent, etc. catch the throw and set data members with the error type and message 
// Some of the messages are constructed with transient information 
//      (e.g., contents of a bad GenParticle)
//
// ----------------------------------------------------------------------



#include <stdexcept>
 
namespace HepMC {

//! IO exception handling

///
/// \class  IO_Exception
/// IO_GenEvent, etc. catch the throw and set data members with the error type and message 
/// Some of the messages are constructed with transient information 
///      (e.g., contents of a bad GenParticle)
class IO_Exception : public std::runtime_error {
public:
  IO_Exception(const std::string & msg) 
  : std::runtime_error(msg) { }

  /// IO error types
  enum ErrorType{ OK,
                  NullEvent, 
                  WrongFileType, 
                  MissingStartKey, 
		  EndOfStream, 
		  EndKeyMismatch, 
		  MissingEndKey, 
		  InvalidData,
                  InputAndOutput,
		  BadOutputStream,
		  BadInputStream };

};

}	// namespace HepMC

#endif // HEPMC_IO_EXCEPTION_H
HepMC/IO_BaseClass.h000644 000765 000024 00000007300 13310701701 014413 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HEPMC_IO_BASECLASS_H
#define HEPMC_IO_BASECLASS_H

//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, November 1999, refer to:
// M. Dobbs and J.B. Hansen, "The HepMC C++ Monte Carlo Event Record for
// High Energy Physics", Computer Physics Communications (to be published).
//
// event input/output base class
//////////////////////////////////////////////////////////////////////////
//
// class from which all input/output classes shall inherit from.
// i.e.: if you want to write events to hbook ntuples,
//              then inherit from this class and re-define read_event()
//              and write_event()
//
// (Possible extension: Could make this an input iterator)
//

#include <iostream>
#include "HepMC/GenEvent.h"

namespace HepMC {

    //! all input/output classes inherit from IO_BaseClass

    ///
    /// \class  IO_BaseClass
    /// If you want to write a new IO class, 
    /// then inherit from this class and re-define read_event()
    /// and write_event()
    ///
    class IO_BaseClass {
    public:
        virtual ~IO_BaseClass() {}

	/// write this GenEvent
	virtual void write_event( const GenEvent* ) =0;
	/// fill this GenEvent
	virtual bool fill_next_event( GenEvent* ) =0;
	/// write output to ostr
	virtual void print( std::ostream& ostr = std::cout ) const;
	//
	// the read_next_event() differs from
	// the fill_***() methods in that it creates a new event
	// before calling the  corresponding fill_*** method
	// (they are not intended to be over-ridden)
	GenEvent*    read_next_event();  //!< do not over-ride
	//
	// The overloaded stream operators >>,<< are identical to
	//   read_next_event and write_event methods respectively.
	//   (or read_particle_data_table and write_particle_data_table)
	// the event argument for the overloaded stream operators is a pointer,
	// which is passed by reference.
	//  i.e.  GenEvent* evt; 
	//        io >> evt; 
	// will give the expected result.
	// (note: I don't see any reason to have separate const and non-const
	//  versions of operator<<, but the pedantic ansi standard insists 
	//  on it) 
	/// the same as read_next_event
	virtual       GenEvent*& operator>>( GenEvent*& );
	/// the same as write_event
	virtual const GenEvent*& operator<<( const GenEvent*& );
	/// the same as write_event
	virtual       GenEvent*& operator<<( GenEvent*& );
    };

    //////////////
    // Inlines  //
    //////////////

    inline GenEvent* IO_BaseClass::read_next_event() {
	/// creates a new event and fills it by calling 
	/// the sister method read_next_event( GenEvent* )
	// 
        // 1. create an empty event container
        GenEvent* evt = new GenEvent();
	// 2. fill the evt container - if the read is successful, return the
	//    pointer, otherwise return null and delete the evt
	if ( fill_next_event( evt ) ) return evt;
	// note: the below delete is only reached if read fails
	//       ... thus there is not much overhead in new then delete 
	//       since this statement is rarely reached
	delete evt;
	return 0;
    }

    inline void IO_BaseClass::print( std::ostream& ostr ) const { 
	ostr << "IO_BaseClass: abstract parent I/O class. " <<  std::endl;
    }

    inline GenEvent*& IO_BaseClass::operator>>( GenEvent*& evt ){
	evt = read_next_event();
	return evt;
    }

    inline const GenEvent*& IO_BaseClass::operator<<(
					      const GenEvent*& evt ) {
	write_event( evt );
	return evt;
    }

    inline GenEvent*& IO_BaseClass::operator<<( GenEvent*& evt ) {
	write_event( evt );
	return evt;
    }

} // HepMC

#endif  // HEPMC_IO_BASECLASS_H
//--------------------------------------------------------------------------



HepMC/PythiaWrapper6_4.h000644 000765 000024 00000014070 13310701701 015276 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef PYTHIA_WRAPPER_6_H
#define PYTHIA_WRAPPER_6_H

//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, November 2000
// Version 6.200 update October 2001
// Wrapper for FORTRAN version of Pythia 
// This wrapper is NOT intended as a part of HepMC - it is only supplied
// for your convenience.
//////////////////////////////////////////////////////////////////////////
// 
// A simple example of calling Pythia from C++ using this header file is
// given in test/test_PythiaWrapper.cxx
//
// Note the pyhepc routine is used by Pythia to fill
// the HEPEVT common block uses double precision and 4000 entries.
//

#include <ctype.h>
#include <cstring>

//--------------------------------------------------------------------------
// Initialization routine


extern "C" {
    void initpydata(void);
}
#define initpydata initpydata_

//--------------------------------------------------------------------------
// PYTHIA Common Block Declarations

const int pyjets_maxn =4000;
extern "C" {
    extern struct {
	int n, npad, k[5][pyjets_maxn];
	double p[5][pyjets_maxn], v[5][pyjets_maxn];
    } pyjets_;
}
#define pyjets pyjets_

extern "C" {
    extern struct {
        int mstu[200];
        double paru[200];
        int mstj[200];
        double parj[200];
    } pydat1_;
}
#define pydat1 pydat1_

extern "C" {
    extern struct {
	int kchg[4][500];
	double pmas[4][500], parf[2000], vckm[4][4];  
    } pydat2_;
}
#define pydat2 pydat2_

extern "C" {
    extern struct {
	int mdcy[3][500], mdme[2][8000];
	double brat[8000];
	int kfdp[5][8000];
    } pydat3_;
}
#define pydat3 pydat3_

extern "C" {
    extern struct {
	int mrpy[6];
	double rrpy[100];
    } pydatr_;
}
#define pydatr pydatr_

extern "C" {
    extern struct {
	int msel, mselpd, msub[500], kfin[81][2];
	double ckin[200];
    } pysubs_;
}
#define pysubs pysubs_

extern "C" {
    extern struct {
	int mstp[200];
	double parp[200];
	int msti[200];
	double pari[200];
    } pypars_;
}
#define pypars pypars_

extern "C" {
    extern struct {
	int mint[400];
	double vint[400];
    } pyint1_;
}
#define pyint1 pyint1_

extern "C" {
    extern struct {
	int iset[500], kfpr[2][500];
	double coef[20][500];
	int icol[2][4][40];
    } pyint2_;
}
#define pyint2 pyint2_

extern "C" {
    extern struct pin3 {
        double xsfx[81][2];	// Fortran is xsfx(2,-40:40)
	int isig[3][1000];
	double sigh[1000];
    } pyint3_;
}
#define pyint3 pyint3_

extern "C" {
    extern struct {
        int mwid[500];
	double wids[5][500];
    } pyint4_;
}
#define pyint4 pyint4_

extern "C" {
    extern struct pin5 {
	int ngenpd, ngen[3][501];	// Fortran is ngen(0:500,3)
	double xsec[3][501];		// Fortran is xsec(0:500,3)
    } pyint5_;
}
#define pyint5 pyint5_

extern "C" {
    extern struct pin7 {
        double sigt[6][7][7];	// Fortran is sigt(0:6,0:6,0:5)
    } pyint7_;
}
#define pyint7 pyint7_

extern "C" {
    extern struct pin8 {
	double xpvmd[13];	// Fortran is xpvmd(-6:6)
	double xpanl[13];	// Fortran is xpanl(-6:6)
	double xpanh[13];	// Fortran is xpanh(-6:6)
	double xpbeh[13];	// Fortran is xpbeh(-6:6)
	double xpdir[13];	// Fortran is xpdir(-6:6)
    } pyint8_;
}
#define pyint8 pyint8_

extern "C" {
    extern struct pin9 {
	double vxpvmd[13];	// Fortran is vxpvmd(-6:6)
	double vxpanl[13];	// Fortran is vxpanl(-6:6)
	double vxpanh[13];	// Fortran is vxpanh(-6:6)
	double vxpdgm[13];	// Fortran is vxpdgm(-6:6)
    } pyint9_;
}
#define pyint9 pyint9_

extern "C" {
    extern struct pssm {
	int imss[100];		// Fortran is imss(0:99)
	double rmss[100];	// Fortran is rmss(0:99)
    } pyssm_;
}
#define pyssm pyssm_

extern "C" {
    extern struct {
        double zmix[4][4];
	double umix[2][2];
	double vmix[2][2];
	double smz[4];
	double smw[2];
	double sfmix[4][16];
        double zmixi[4][4];
	double umixi[2][2];
	double vmixi[2][2];
    } pyssmt_;
}
#define pyssmt pyssmt_

extern "C" {
    extern struct {
        double rvlam[3][3][3];
        double rvlamp[3][3][3];
        double rvlamb[3][3][3];
    } pymsrv_;
}
#define pymsrv pymsrv_

extern "C" {
    extern struct prvnv {
        double ab[2][16][2];
	double rms[4];		// Fortran is rms(0:3)
	double res[5][6];
	int idr;
	int idr2;
	double dcmass;
	int kfr[3];
    } pyrvnv_;
}
#define pyrvnv pyrvnv_

extern "C" {
    extern struct prvpm {
        double rm[4];		// Fortran is rm(0:3)
	double a[2];
	double b[2];
	double resm[2];
	double resw[2];
	bool mflag;
    } pyrvpm_;
}
#define pyrvpm pyrvpm_

extern "C" {
    extern struct {
        double xxm[20];
    } pyints_;
}
#define pyints pyints_

extern "C" {
    extern struct {
        double x1;
    } pyg2dx_;
}
#define pyg2dx pyg2dx_

//--------------------------------------------------------------------------
// PYTHIA routines declaration

#define pyhepc pyhepc_ 
#define pyinit pyinit_
#define pylist pylist_
#define pystat pystat_
#define pyevnt pyevnt_
#define upinit upinit_
#define upevnt upevnt_
    extern "C" {
	void pyhepc(int*);
	void pyinit(const char*,const char*,const char*,double*,int,int,int);
	void pylist(int*);
	void pystat(int*);
	void pyevnt();
	void upinit();
	void upevnt();
    }

// define methods to hide the subtle syntax necessary to call fortran from C++
inline void call_pyhepc( int mode ){ pyhepc( &mode ); }
inline void call_pyinit( const char* frame, const char* beam, const char* target,
		  double win ) 
{ pyinit( frame,beam,target,&win,strlen(frame),strlen(beam),strlen(target) ); }
inline void call_pylist( int mode ){ pylist( &mode ); }
inline void call_pystat( int mode ){ pystat( &mode ); }
inline void call_pyevnt(){ pyevnt(); }


//--------------------------------------------------------------------------
// PYTHIA block data
// ( with gcc it works to initialize the block data by calling 
//   "pydata();" at beginning, but this fails for f77, so the fortran routine
//   initpydata.f is supplied ... call it instead for platform independent
//   behaviour )

#define pydata pydata_
extern "C" {
    void pydata(void);
}

#endif  // PYTHIA_WRAPPER_6_H
//--------------------------------------------------------------------------
HepMC/GenParticle.h000644 000765 000024 00000023665 13310701701 014375 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HEPMC_GEN_PARTICLE_H
#define HEPMC_GEN_PARTICLE_H

//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, September 1999, refer to:
// M. Dobbs and J.B. Hansen, "The HepMC C++ Monte Carlo Event Record for
// High Energy Physics", Computer Physics Communications (to be published).
//
// particle within an event coming in/out of a vertex
// particle is the basic building block or unit of the event record
//////////////////////////////////////////////////////////////////////////
//
// example:
//      GenParticle* p = new GenParticle( FourVector(1,1,1,3), 11, 1 );
// creates a particle with 4-vector (p,E)=1,1,1,3 - with pdg id 11 (electron)
// and give this particle status =1.
//
// the pointers to end/production vertices can only be set by the
//  vertices themselves - thus to set the production vertex for a particle,
//  you add the particle to that vertex with GenVertex::add_particle_out()
//
// We decide not to have a separate 4 vector for the momentum 
//  at decay time (which MC++ includes to allow dE/dX losses etc). 
//  If you want that, just add a decay vertex with the
//  same particle (modified momentum) going out
//

#include "HepMC/Flow.h"
#include "HepMC/Polarization.h"
#include "HepMC/SimpleVector.h"
#include "HepMC/IteratorRange.h"
#include <iostream>
#ifdef _WIN32
#define hepmc_uint64_t  __int64
#else
#include <stdint.h>	// for uint64_t
#define hepmc_uint64_t   uint64_t
#endif

namespace HepMC {

    class GenVertex;
    class GenEvent; 

    class GenParticleProductionRange;
    class ConstGenParticleProductionRange;
    class GenParticleEndRange;
    class ConstGenParticleEndRange;

    //! The GenParticle class contains information about generated particles

    ///
    /// \class GenParticle 
    /// HepMC::GenParticle 
    /// contains momentum, generated mass, particle ID, decay status, 
    /// flow, polarization, pointers to production and decay vertices
    /// and a unique barcode identfier.
    ///
    class GenParticle {

	friend class GenVertex; // so vertex can set decay/production vertexes
	friend class GenEvent;  // so event can set the barCodes
	/// print particle
	friend std::ostream& operator<<( std::ostream&, const GenParticle& );

    public:
        /// default constructor
        GenParticle(void);
        /// constructor requires momentum and particle ID
 	GenParticle( const FourVector& momentum, int pdg_id,
		     int status = 0, const Flow& itsflow = Flow(),
		     const Polarization& polar = Polarization(0,0) );
	GenParticle( const GenParticle& inparticle ); //!< shallow copy.
	virtual ~GenParticle();

        void swap( GenParticle & other); //!< swap
	GenParticle& operator=( const GenParticle& inparticle ); //!< shallow.
        /// check for equality
	bool         operator==( const GenParticle& ) const;
        /// check for inequality
	bool         operator!=( const GenParticle& ) const;

	/// dump this particle's full info to ostr
	void       print( std::ostream& ostr = std::cout ) const; 

	operator HepMC::FourVector() const; //!< conversion operator

	////////////////////
	// access methods //
	////////////////////

	/// standard 4 momentum
	const FourVector &          momentum() const;
	/// particle ID
	int                  pdg_id() const;
	/// HEPEVT decay status
	int                  status() const;
	/// particle flow
	const Flow &         flow() const;
	/// particle flow index
	int                  flow( int code_index ) const;
        /// polarization information
	const Polarization & polarization() const;
	/// pointer to the production vertex
	GenVertex*           production_vertex() const;
	/// pointer to the decay vertex
	GenVertex*           end_vertex() const;
	/// pointer to the event that owns this particle
	GenEvent*            parent_event() const;

        /// Because of precision issues, the generated mass is not always the 
	/// same as the mass calculated from the momentum 4 vector.
        /// If the generated mass has been set, then generated_mass() 
	/// returns that value.
        /// If the generated mass has not been set, then generated_mass() 
	/// returns the mass calculated from the momentum 4 vector.
        double               generated_mass() const; //!< mass as generated

	/// generatedMass() is included for backwards compatibility with CLHEP HepMC
        double               generatedMass() const { return generated_mass(); }


	///
	/// The barcode is the particle's reference number, every vertex in the
	/// event has a unique barcode. Particle barcodes are positive numbers,
	/// vertex barcodes are negative numbers.
	/// 
	/// Please note that the barcodes are intended for internal use within 
	/// HepMC as a unique identifier for the particles and vertices.
	/// Using the barcode to encode extra information is an abuse of 
	/// the barcode data member and causes confusion among users. 
	/// 
	int                  barcode() const; //!< particle barcode
	
	/// Convenience method.  Returns true if status==1
	bool                 is_undecayed() const;
	/// Convenience method.  Returns true if status==2
	bool                 has_decayed() const;
	/// Convenience method.  Returns true if status==4
	/// Note that using status 4 for beam particles is a new convention which
	/// may not have been implemented by the code originating this GenEvent.
	bool                 is_beam() const;

	/// incoming particle range
	GenParticleProductionRange particles_in( IteratorRange range = relatives );
	/// incoming particle range
	ConstGenParticleProductionRange particles_in( IteratorRange range = relatives ) const;
	/// outgoing particle range
	GenParticleEndRange particles_out( IteratorRange range = relatives );
	/// outgoing particle range
	ConstGenParticleEndRange particles_out( IteratorRange range = relatives ) const;

	/////////////////////
	// mutator methods //
	/////////////////////

	/// In general there is no reason to "suggest_barcode"
	bool                 suggest_barcode( int the_bar_code );

	void   set_momentum( const FourVector& vec4 ); //!< set standard 4 momentum
	void   set_pdg_id( int id ); //!< set particle ID
	void   set_status( int status = 0 ); //!< set decay status
	void   set_flow( const Flow& f ); //!< set particle flow
	void   set_flow( int code_index, int code = 0 ); //!< set particle flow index
	/// set polarization
	void   set_polarization( const Polarization& pol = Polarization(0,0) );
        ///  If you do not call set_generated_mass(), then 
        ///  generated_mass() will simply return the mass calculated from momentum()
        void   set_generated_mass( const double & m ); //!< define the actual generated mass

	///  setGeneratedMass() is included for backwards compatibility with CLHEP HepMC
        void   setGeneratedMass( const double & m )  
	                 { return set_generated_mass(m); }

    protected: // for internal use only by friend GenVertex class

	//static unsigned int counter(); //!< temporary for debugging

        /// set production vertex - for internal use only
	void   set_production_vertex_( GenVertex* productionvertex = 0);
        /// set decay vertex - for internal use only
	void   set_end_vertex_( GenVertex* decayvertex = 0 );
	void   set_barcode_( int the_bar_code ); //!< for use by GenEvent only

        /// scale the momentum vector and generated mass 
        /// this method is only for use by GenEvent
	void convert_momentum( const double& );

    private:
	FourVector       m_momentum;          // momentum vector
	int              m_pdg_id;            // id according to PDG convention
	int              m_status;            // As defined for HEPEVT
	Flow             m_flow;
	Polarization     m_polarization;
	GenVertex*       m_production_vertex; // null if vacuum or beam
	GenVertex*       m_end_vertex;        // null if not-decayed
	int              m_barcode;           // unique identifier in the event
        double           m_generated_mass;    // mass of this particle when it was generated

	//static unsigned int       s_counter;
    };  

    //////////////
    // INLINES  //
    //////////////

    inline GenParticle::operator HepMC::FourVector() const 
    { return m_momentum; }

    inline const FourVector & GenParticle::momentum() const 
    { return m_momentum; }

    inline int GenParticle::pdg_id() const { return m_pdg_id; }

    inline int GenParticle::status() const { return m_status; }

    inline GenVertex* GenParticle::production_vertex() const 
    { return m_production_vertex; }

    inline GenVertex* GenParticle::end_vertex() const { return m_end_vertex; }

    inline const Flow & GenParticle::flow() const { return m_flow; }

    inline int GenParticle::flow( int code_index ) const
    { return m_flow.icode( code_index ); }

    inline const Polarization & GenParticle::polarization() const 
    { return m_polarization; }

    inline void GenParticle::set_momentum( const FourVector& vec4 )
    { m_momentum = vec4; }

    inline void GenParticle::set_pdg_id( int id ) { m_pdg_id = id; }

    inline void GenParticle::set_status( int st ) { m_status = st; }

    inline void GenParticle::set_flow( const Flow& f ) { m_flow = f; }

    inline void GenParticle::set_flow( int code_index, int code ) 
    {
	if ( code == 0 ) { 
	    m_flow.set_unique_icode( code_index );
	} else { 
	    m_flow.set_icode( code_index, code );
	}
    }

    inline void GenParticle::set_polarization( const Polarization& polar )
    { m_polarization = polar; }

    inline int  GenParticle::barcode() const { return m_barcode; }

    inline void GenParticle::set_barcode_( int bc ) { m_barcode = bc; }

    inline bool GenParticle::is_undecayed() const {
        return ( m_status==1 ) ?  true : false;
    }
    inline bool GenParticle::has_decayed() const {
        return ( m_status==2 ) ?  true : false;
    }
    inline bool GenParticle::is_beam() const {
        return ( m_status==4 ) ?  true : false;
    }

} // HepMC

#endif  // HEPMC_GEN_PARTICLE_H
//--------------------------------------------------------------------------

HepMC/WeightContainer.h000644 000765 000024 00000015775 13310701701 015275 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HEPMC_WEIGHT_CONTAINER_H
#define HEPMC_WEIGHT_CONTAINER_H

//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, November 2000, refer to:
// M. Dobbs and J.B. Hansen, "The HepMC C++ Monte Carlo Event Record for
// High Energy Physics", Computer Physics Communications (to be published).
//
// Container for the Weights associated with an event or vertex.
//
// This implementation adds a map-like interface in addition to the 
// vector-like interface.
//////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <vector>
#include <string>
#include <map>

namespace HepMC {

    //! Container for the Weights associated with an event or vertex.

    ///
    /// \class  WeightContainer
    /// This class has both map-like and vector-like functionality.
    /// Named weights are now supported.
    class WeightContainer {
	friend class GenEvent;

    public:
        /// defining the size type used by vector and map
	typedef std::size_t size_type;
        /// iterator for the weight container
	typedef std::vector<double>::iterator iterator;
        /// const iterator for the weight container
	typedef std::vector<double>::const_iterator const_iterator;
	
        /// default constructor
	explicit WeightContainer( size_type n = 0, double value = 0. );
        /// construct from a vector of weights
	WeightContainer( const std::vector<double>& weights );
        /// copy
	WeightContainer( const WeightContainer& in );
	~WeightContainer();

        /// swap
        void swap( WeightContainer & other);
        /// copy assignment
	WeightContainer& operator=( const WeightContainer& );
        /// alternate assignment using a vector of doubles
	WeightContainer& operator=( const std::vector<double>& in );

        /// print weights
	void          print( std::ostream& ostr = std::cout ) const;
        /// write weights in a readable table
	void          write( std::ostream& ostr = std::cout ) const;

        /// size of weight container
	size_type     size() const;
	/// return true if weight container is empty
	bool          empty() const;
	/// push onto weight container
	void          push_back( const double& );
	/// pop from weight container
	void          pop_back();
	/// clear the weight container
	void          clear();

	/// check to see if a name exists in the map
	bool          has_key( const std::string& s ) const;

        /// access the weight container
	double&       operator[]( size_type n );  // unchecked access
        /// access the weight container
	const double& operator[]( size_type n ) const;
        /// access the weight container
	double&       operator[]( const std::string& s );  // unchecked access
        /// access the weight container
	const double& operator[]( const std::string& s ) const;

        /// equality
	bool operator==( const WeightContainer & ) const;
        /// inequality
	bool operator!=( const WeightContainer & ) const;
 	
	/// returns the first element
	double&       front();
	/// returns the first element
	const double& front() const;   
	/// returns the last element
	double&       back();
	/// returns the last element
	const double& back() const;

	/// begining of the weight container
	iterator            begin();
	/// end of the weight container
	iterator            end();
	/// begining of the weight container
	const_iterator      begin() const;
	/// end of the weight container
	const_iterator      end() const;

    private:
        // for internal use only

        /// maplike iterator for the weight container
	/// for internal use only
	typedef std::map<std::string,size_type>::iterator       map_iterator;
        /// const iterator for the weight container
	/// for internal use only
	typedef std::map<std::string,size_type>::const_iterator const_map_iterator;
	/// begining of the weight container
	/// for internal use only
	map_iterator            map_begin();
	/// end of the weight container
	/// for internal use only
	map_iterator            map_end();
	/// begining of the weight container
	/// for internal use only
	const_map_iterator      map_begin() const;
	/// end of the weight container
	/// for internal use only
	const_map_iterator      map_end() const;
	
	/// used by the constructors to set initial names
	/// for internal use only
	void set_default_names( size_type n );
	
    private:
	std::vector<double>          m_weights;
	std::map<std::string,size_type> m_names;
    };

    ///////////////////////////
    // INLINES               //
    ///////////////////////////

    inline WeightContainer::WeightContainer( const WeightContainer& in )
	: m_weights(in.m_weights), m_names(in.m_names)
    {}

    inline WeightContainer::~WeightContainer() {}

    inline void WeightContainer::swap( WeightContainer & other)
    { 
        m_weights.swap( other.m_weights ); 
        m_names.swap( other.m_names ); 
    }

    inline WeightContainer& WeightContainer::operator=
    ( const WeightContainer& in ) {
        /// best practices implementation
	WeightContainer tmp( in );
	swap( tmp );
	return *this;
    }

    inline WeightContainer& WeightContainer::operator=
    ( const std::vector<double>& in ) {
        /// best practices implementation
	WeightContainer tmp( in );
	swap( tmp );
	return *this;
    }

    inline WeightContainer::size_type WeightContainer::size() const { return m_weights.size(); }

    inline bool WeightContainer::empty() const { return m_weights.empty(); }

    inline void WeightContainer::clear() 
    { 
	m_weights.clear(); 
	m_names.clear(); 
    }

    inline double& WeightContainer::operator[]( size_type n ) 
    { return m_weights[n]; }

    inline const double& WeightContainer::operator[]( size_type n ) const
    { return m_weights[n]; }

    inline double& WeightContainer::front() { return m_weights.front(); }

    inline const double& WeightContainer::front() const 
    { return m_weights.front(); }

    inline double& WeightContainer::back() { return m_weights.back(); }

    inline const double& WeightContainer::back() const 
    { return m_weights.back(); }

    inline WeightContainer::iterator WeightContainer::begin() 
    { return m_weights.begin(); }

    inline WeightContainer::iterator WeightContainer::end() 
    { return m_weights.end(); }

    inline WeightContainer::const_iterator WeightContainer::begin() const 
    { return m_weights.begin(); }

    inline WeightContainer::const_iterator WeightContainer::end() const 
    { return m_weights.end(); }

    inline WeightContainer::map_iterator WeightContainer::map_begin() 
    { return m_names.begin(); }

    inline WeightContainer::map_iterator WeightContainer::map_end() 
    { return m_names.end(); }

    inline WeightContainer::const_map_iterator WeightContainer::map_begin() const 
    { return m_names.begin(); }

    inline WeightContainer::const_map_iterator WeightContainer::map_end() const 
    { return m_names.end(); }

} // HepMC

#endif  // HEPMC_WEIGHT_CONTAINER_H
//--------------------------------------------------------------------------



HepMC/PythiaWrapper6_4_WIN32.h000644 000765 000024 00000010065 13310701701 016160 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifdef  _WIN32  // This version is for Windows MS Visual C++ only.
#ifndef PYTHIA_WRAPPER_6_H
#define PYTHIA_WRAPPER_6_H

//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, November 2000
// Version 6.200 update October 2001
// Wrapper for FORTRAN version of Pythia 
// The _WIN32 version is provided by Witold Pokorski 
//       <Witold.Pokorski@Cern.CH>, 2002-01-22, and is the version 
//       which should be used when compiling on Windows MS Visual C++.
// This wrapper is NOT intended as a part of HepMC - it is only supplied
// for your convenience.
//////////////////////////////////////////////////////////////////////////
// 
// A simple example of calling Pythia from C++ using this header file is
// given in test/test_PythiaWrapper.cxx
//
// Note the pyhepc routine is used by Pythia to fill
// the HEPEVT common block uses double precision and 4000 entries.
//

#include <ctype.h>
#include <string.h>

//--------------------------------------------------------------------------
// Initialization routine


extern "C" {
    void INITPYDATA(void);
}
#define initpydata INITPYDATA


//--------------------------------------------------------------------------
// PYTHIA Common Block Declarations

const int pyjets_maxn =4000;
struct PYJETS_DEF {
	int n, npad, k[5][pyjets_maxn];
	double p[5][pyjets_maxn], v[5][pyjets_maxn];
};

struct PYDAT1_DEF{
        int mstu[200];
        double paru[200];
        int mstj[200];
        double parj[200];
};

struct PYDAT2_DEF{
	int kchg[4][500];
	double pmas[4][500], parf[2000], vckm[4][4];  
};

struct PYDAT3_DEF{
	int mdcy[3][500], mdme[2][8000];
	double brat[8000];
	int kfdp[5][8000];
};

struct PYDATR_DEF{
	int mrpy[6];
	double rrpy[100];
    };

struct PYSUBS_DEF{
	int msel, mselpd, msub[500], kfin[81][2];
	double ckin[200];
};

struct PYPARS_DEF{
	int mstp[200];
	double parp[200];
	int msti[200];
	double pari[200];
};

struct PYINT1_DEF{
	int mint[400];
	double vint[400];
};

struct PYINT2_DEF{
	int iset[500], kfpr[2][500];
	double coef[20][500];
	int icol[2][4][40];       // was [320] was [40][4][2]
};

struct PYINT5_DEF{
	int ngenpd, ngen[3][501];
	double xsec[3][501];
};


extern "C" PYSUBS_DEF PYSUBS;
extern "C" PYJETS_DEF PYJETS;
extern "C" PYDAT1_DEF PYDAT1;
extern "C" PYDAT2_DEF PYDAT2;
extern "C" PYPARS_DEF PYPARS;
extern "C" PYDATR_DEF PYDATR;
extern "C" PYDAT3_DEF PYDAT3;
extern "C" PYINT1_DEF PYINT1;
extern "C" PYINT2_DEF PYINT2;
extern "C" PYINT5_DEF PYINT5;


#define pysubs PYSUBS
#define pyjets PYJETS
#define pydat1 PYDAT1
#define pydat2 PYDAT2
#define pypars PYPARS
#define pydatr PYDATR
#define pydat3 PYDAT3
#define pyint1 PYINT1
#define pyint2 PYINT2
#define pyint5 PYINT5



//--------------------------------------------------------------------------
// PYTHIA routines declaration


extern "C" {
	void __stdcall PYHEPC(int*);
	void __stdcall PYINIT(const char*,int,const char*,int,const char*,int,double*);
	void __stdcall PYLIST(int*);
	void __stdcall PYSTAT(int*);
	void __stdcall PYEVNT();
	void __stdcall UPINIT();
	void __stdcall UPEVNT();
}

//--------------------------------------------------------------------------
// PYTHIA block data
// ( with gcc it works to initialize the block data by calling 
//   "pydata();" at beginning, but this fails for f77, so the fortran routine
//   initpydata.f is supplied ... call it instead for platform independent
//   behaviour )


extern "C" {
    void __stdcall PYDATA(void);
}

inline void call_pyhepc( int mode ){ PYHEPC( &mode ); }
inline void call_pyinit( const char* frame, const char* beam, const char* target,
		  double win ) 
{ PYINIT( frame,strlen(frame),beam,strlen(beam),target,strlen(target),&win); }
inline void call_pylist( int mode ){ PYLIST( &mode ); }
inline void call_pystat( int mode ){ PYSTAT( &mode ); }
inline void call_pyevnt(){ PYEVNT(); }
inline void call_upinit(){ UPINIT(); }
inline void call_upevnt(){ UPEVNT(); }




#endif  // PYTHIA_WRAPPER_6_H
#endif  // _WIN32
//--------------------------------------------------------------------------
HepMC/IO_AsciiParticles.h000644 000765 000024 00000007536 13310701701 015465 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HEPMC_IO_ASCIIPARTICLES_H
#define HEPMC_IO_ASCIIPARTICLES_H

//////////////////////////////////////////////////////////////////////////
// Mikhail.Kirsanov@Cern.CH, 2006
// event input/output in ascii format for eye and machine reading
//////////////////////////////////////////////////////////////////////////
//
// Strategy for reading or writing events as machine readable
//  ascii to a file. When instantiating, the mode of file to be created 
//  must be specified. Options are:
//      std::ios::in     open file for input 
//      std::ios::out    open file for output
//      std::ios::trunc  erase old file when opening (i.e. ios::out|ios::trunc
//                    removes oldfile, and creates a new one for output )
//      std::ios::app    append output to end of file
//  for the purposes of this class, simultaneous input and output mode 
//  ( std::ios::in | std::ios::out ) is not allowed.
// 
// Event listings are preceded by the key:
//  "HepMC::IO_AsciiParticles-START_EVENT_LISTING\n"
//  and terminated by the key:
//  "HepMC::IO_AsciiParticles-END_EVENT_LISTING\n"
// Comments are allowed. They need not be preceded by anything, though if
//  a comment is written using write_comment( const string ) then it will be
//  preceded by "HepMC::IO_AsciiParticles-COMMENT\n"
// Each event, vertex, particle, particle data is preceded by 
//  "E ","V ","P ","D "    respectively.
// Comments may appear anywhere in the file -- so long as they do not contain
//  any of the 4 start/stop keys.
//

#include <fstream>
#include <string>
#include <map>
#include <vector>
#include "HepMC/IO_BaseClass.h"

namespace HepMC {

    class GenEvent;
    class GenVertex;
    class GenParticle;

    //! event input/output in ascii format for eye and machine reading

    ///
    /// \class IO_AsciiParticles
    /// Strategy for reading or writing events as machine readable
    ///  ascii to a file. When instantiating, the mode of file to be created 
    ///  must be specified. 
    ///
    class IO_AsciiParticles : public IO_BaseClass {
    public:
        /// constructor requiring a file name and std::ios mode
	IO_AsciiParticles( const char* filename="IO_AsciiParticles.dat", 
		  std::ios::openmode mode=std::ios::out );
	virtual       ~IO_AsciiParticles();

        /// write this event
  	void          write_event( const GenEvent* evt );
        /// get the next event
	bool          fill_next_event( GenEvent* evt );
	/// insert a comment directly into the output file --- normally you
	///  only want to do this at the beginning or end of the file. All
	///  comments are preceded with "HepMC::IO_AsciiParticles-COMMENT\n"
  	void          write_comment( const std::string comment );

        /// set output precision
        void          setPrecision(int iprec);

	int           rdstate() const;  //!< check the state of the IO stream
	void          clear();  //!< clear the IO stream

        /// write to ostr
	void          print( std::ostream& ostr = std::cout ) const;

    protected: // for internal use only
	/// write end tag
	bool          write_end_listing();
    private: // use of copy constructor is not allowed
	IO_AsciiParticles( const IO_AsciiParticles& ) : IO_BaseClass() {}
    private: // data members
    int                 m_precision;
	std::ios::openmode  m_mode;
	std::fstream*       m_file;
    std::ostream*       m_outstream;
	bool                m_finished_first_event_io;
    };

    //////////////
    // Inlines  //
    //////////////

    inline int  IO_AsciiParticles::rdstate() const { return (int)m_file->rdstate(); }
    inline void IO_AsciiParticles::clear() { m_file->clear(); }
    inline void IO_AsciiParticles::setPrecision(int iprec) { m_precision=iprec; }

} // HepMC

#endif  // HEPMC_IO_ASCIIPARTICLES_H
//--------------------------------------------------------------------------
HepMC/GenCrossSection.h000644 000765 000024 00000006564 13310701701 015247 0ustar00sa639staff000000 000000 #ifndef HEPMC_GEN_CROSS_SECTION_H
#define HEPMC_GEN_CROSS_SECTION_H

//--------------------------------------------------------------------------
//////////////////////////////////////////////////////////////////////////
// garren@fnal.gov, May 2009
// 
//////////////////////////////////////////////////////////////////////////
//--------------------------------------------------------------------------

#include <iostream>

namespace HepMC {

//! The GenCrossSection class stores the generated cross section

///
/// \class  GenCrossSection
/// HepMC::GenCrossSection is used to store the generated cross section.
/// This class is meant to be used to pass, on an event by event basis,
/// the current best guess of the total cross section. 
/// It is expected that the final cross section will be stored elsewhere.
///
///   - double cross_section;  		// cross section in pb
///   - double cross_section_error;  	// error associated with this cross section   
/// 
/// The units of cross_section and cross_section_error are expected to be pb.
/// 
/// GenCrossSection information will be written if GenEvent contains a pointer 
/// to a valid GenCrossSection object.
///
class GenCrossSection {

public:
  GenCrossSection()
    : m_cross_section(0),
      m_cross_section_error(0),
      m_is_set(false)
    {}
  ~GenCrossSection() {}

  GenCrossSection( GenCrossSection const & orig ); //!< copy

  void swap( GenCrossSection & other); //!< swap
  GenCrossSection &  operator = ( GenCrossSection const & rhs ); //!< shallow
  /// check for equality
  bool         operator==( const GenCrossSection& ) const;
  /// check for inequality
  bool         operator!=( const GenCrossSection& ) const;


  // ---  accessors:

  /// cross section in pb
  double cross_section()          const { return m_cross_section; }
  /// error associated with this cross section in pb
  double cross_section_error()    const { return m_cross_section_error; }

  /// True if the cross section has been set.  False by default.
  bool   is_set()                 const { return m_is_set; }

  // ---  mutators:
  /// Set cross section and error in pb
  void   set_cross_section( double xs, double xs_err );
  /// set cross section in pb
  void   set_cross_section( double );
  /// set error associated with this cross section in pb
  void   set_cross_section_error( double );
  /// Clear all GenCrossSection info 
  /// (disables output of GenCrossSection until the cross section is set again)
  void   clear();
 
  // ---  I/O:
  /// write to an output stream
  std::ostream &  write( std::ostream & ) const;
  /// read from an input stream
  std::istream &  read( std::istream & );

private: // data members
    double m_cross_section;
    double m_cross_section_error;
    bool   m_is_set;

};

//
// streaming I/O

inline std::ostream & operator << ( std::ostream & os, GenCrossSection & xs )
{ return xs.write(os); }

inline std::istream & operator >> ( std::istream & is, GenCrossSection & xs )
{ return xs.read(is); }

//
// inline methods

inline void GenCrossSection::set_cross_section( double xs, double xserr ) { 
  set_cross_section(xs);
  set_cross_section_error(xserr); 
}

inline void GenCrossSection::set_cross_section( double xs )        
{
  m_cross_section = xs;
  m_is_set = true;
}

inline void GenCrossSection::set_cross_section_error( double xserr )  
{
  m_cross_section_error = xserr;
}

} // HepMC

#endif  // HEPMC_GEN_CROSS_SECTION_H
HepMC/HepMCDefs.h000644 000765 000024 00000003170 13310701701 013723 0ustar00sa639staff000000 000000 #ifndef HEPMC_DEFS_H
#define HEPMC_DEFS_H
// ----------------------------------------------------------------------
//
// HepMCDefs.h
// Author:  Lynn Garren
//
//  Define various useful macros
//  Most allow users to check for various code features
//
// ----------------------------------------------------------------------

// the HeavyIon class is available in HepMC  
#ifndef HEPMC_HAS_HEAVY_ION
#define HEPMC_HAS_HEAVY_ION
#endif

// the IO_Ascii class is NOT available in HepMC   
#ifndef HEPMC_IO_ASCII_REMOVED
#define HEPMC_IO_ASCII_REMOVED
#endif

// the ParticleData class is NOT available in HepMC   
#ifndef HEPMC_PARTICLE_DATA_REMOVED
#define HEPMC_PARTICLE_DATA_REMOVED
#endif

// the IO_GenEvent class is available in HepMC   
#ifndef HEPMC_HAS_IO_GENEVENT
#define HEPMC_HAS_IO_GENEVENT
#endif

// the PdfInfo class is available in HepMC  
#ifndef HEPMC_HAS_PDF_INFO
#define HEPMC_HAS_PDF_INFO
#endif

// HepMC uses SimpleVector (FourVector) to store momentum and position  
#ifndef HEPMC_HAS_SIMPLE_VECTOR
#define HEPMC_HAS_SIMPLE_VECTOR
#endif

// units are defined in HepMC  
#ifndef HEPMC_HAS_UNITS
#define HEPMC_HAS_UNITS
#endif

// the GenCrossSection class is available in HepMC  
#ifndef HEPMC_HAS_CROSS_SECTION
#define HEPMC_HAS_CROSS_SECTION
#endif

// the iterator range classes are available in HepMC  
#ifndef HEPMC_HAS_ITERATOR_RANGES
#define HEPMC_HAS_ITERATOR_RANGES
#endif

// the HepMC::WeightContainer class allows named weights
#ifndef HEPMC_HAS_NAMED_WEIGHTS
#define HEPMC_HAS_NAMED_WEIGHTS
#endif

// define the version of HepMC. 
#ifndef HEPMC_VERSION
#define HEPMC_VERSION "2.06.08"
#endif

#endif  // HEPMC_DEFS_H
HepMC/HEPEVT_Wrapper.h000644 000765 000024 00000051575 13310701701 014674 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------

#ifndef HEPEVT_EntriesAllocation
#define HEPEVT_EntriesAllocation 200000
#endif  // HEPEVT_EntriesAllocation

//--------------------------------------------------------------------------
#ifndef HEPMC_HEPEVT_COMMON_H
#define HEPMC_HEPEVT_COMMON_H
//////////////////////////////////////////////////////////////////////////
//
//      PARAMETER (NMXHEP=2000) 
//      COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP), 
//     &        JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
/**********************************************************/
/*           D E S C R I P T I O N :                      */
/*--------------------------------------------------------*/
/* NEVHEP          - event number (or some special meaning*/
/*                    (see documentation for details)     */
/* NHEP            - actual number of entries in current  */
/*                    event.                              */
/* ISTHEP[IHEP]    - status code for IHEP'th entry - see  */
/*                    documentation for details           */
/* IDHEP [IHEP]    - IHEP'th particle identifier according*/
/*                    to PDG.                             */
/* JMOHEP[IHEP][0] - pointer to position of 1st mother    */
/* JMOHEP[IHEP][1] - pointer to position of 2nd mother    */
/* JDAHEP[IHEP][0] - pointer to position of 1st daughter  */
/* JDAHEP[IHEP][1] - pointer to position of 2nd daughter  */
/* PHEP  [IHEP][0] - X momentum                           */
/* PHEP  [IHEP][1] - Y momentum                           */
/* PHEP  [IHEP][2] - Z momentum                           */
/* PHEP  [IHEP][3] - Energy                               */
/* PHEP  [IHEP][4] - Mass                                 */
/* VHEP  [IHEP][0] - X vertex                             */
/* VHEP  [IHEP][1] - Y vertex                             */
/* VHEP  [IHEP][2] - Z vertex                             */
/* VHEP  [IHEP][3] - production time                      */
/*========================================================*/
// Remember, array(1) is the first entry in a fortran array, array[0] is the
//           first entry in a C array.
//
// This interface to HEPEVT common block treats the block as
// an array of bytes --- the precision and number of entries 
// is determined "on the fly" by the wrapper and used to decode
// each entry.
//
// HEPEVT_EntriesAllocation is the maximum size of the HEPEVT common block 
//   that can be interfaced.
//   It is NOT the actual size of the HEPEVT common used in each
//   individual application. The actual size can be changed on
//   the fly using HEPEVT_Wrapper::set_max_number_entries().
// Thus HEPEVT_EntriesAllocation should typically be set
// to the maximum possible number of entries --- 10000 is a good choice
// (and is the number used by ATLAS versions of Pythia).
//
// Note: a statement like    *( (int*)&hepevt.data[0] )
//      takes the memory address of the first byte in HEPEVT,
//      interprets it as an integer pointer, 
//      and dereferences the pointer.
//      i.e. it returns an integer corresponding to nevhep
//

#include <ctype.h>

    const unsigned int hepevt_bytes_allocation = 
                sizeof(long int) * ( 2 + 6 * HEPEVT_EntriesAllocation )
                + sizeof(double) * ( 9 * HEPEVT_EntriesAllocation );


#ifdef _WIN32 // Platform: Windows MS Visual C++
struct HEPEVT_DEF{
        char data[hepevt_bytes_allocation];
    };
extern "C" HEPEVT_DEF HEPEVT;
#define hepevt HEPEVT

#else
extern "C" {
    extern struct {
	char data[hepevt_bytes_allocation];
    } hepevt_;
}
#define hepevt hepevt_

#endif // Platform

#endif  // HEPMC_HEPEVT_COMMON_H

//--------------------------------------------------------------------------
#ifndef HEPMC_HEPEVT_WRAPPER_H
#define HEPMC_HEPEVT_WRAPPER_H

//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, April 24, 2000, refer to:
// M. Dobbs and J.B. Hansen, "The HepMC C++ Monte Carlo Event Record for
// High Energy Physics", Computer Physics Communications (to be published).
//
// Generic Wrapper for the fortran HEPEVT common block
// This class is intended for static use only - it makes no sense to 
// instantiate it.
// Updated: June 30, 2000 (static initialization moved to separate .cxx file)
//////////////////////////////////////////////////////////////////////////
//
// The index refers to the fortran style index: 
// i.e. index=1 refers to the first entry in the HEPEVT common block.
// all indices must be >0
// number_entries --> integer between 0 and max_number_entries() giving total
//                    number of sequential particle indices
// first_parent/child --> index of first mother/child if there is one, 
//                        zero otherwise
// last_parent/child --> if number children is >1, address of last parent/child
//                       if number of children is 1, same as first_parent/child
//                       if there are no children, returns zero.
// is_double_precision --> T or F depending if floating point variables 
//                         are 8 or 4 bytes
//

#include <iostream>
#include <cstdio>       // needed for formatted output using sprintf 

namespace HepMC {

    //! Generic Wrapper for the fortran HEPEVT common block
    
    /// \class HEPEVT_Wrapper
    /// This class is intended for static use only - it makes no sense to 
    /// instantiate it.
    ///
    class HEPEVT_Wrapper {
    public:

        /// write information from HEPEVT common block
	static void print_hepevt( std::ostream& ostr = std::cout );
	/// write particle information to ostr
	static void print_hepevt_particle( int index, 
					   std::ostream& ostr = std::cout );
        static bool is_double_precision();  //!< True if common block uses double

        /// check for problems with HEPEVT common block
	static bool check_hepevt_consistency( std::ostream& ostr = std::cout );

        /// set all entries in HEPEVT to zero
	static void zero_everything();

	////////////////////
	// Access Methods //
	////////////////////
        static int    event_number();             //!< event number
        static int    number_entries();           //!< num entries in current evt
        static int    status( int index );        //!< status code
        static int    id( int index );            //!< PDG particle id
        static int    first_parent( int index );  //!< index of 1st mother
        static int    last_parent( int index );   //!< index of last mother
	static int    number_parents( int index ); //!< number of parents 
        static int    first_child( int index );   //!< index of 1st daughter
        static int    last_child( int index );    //!< index of last daughter
	static int    number_children( int index ); //!< number of children
        static double px( int index );            //!< X momentum       
        static double py( int index );            //!< Y momentum   
        static double pz( int index );            //!< Z momentum   
        static double e( int index );             //!< Energy
        static double m( int index );             //!< generated mass
        static double x( int index );             //!< X Production vertex
        static double y( int index );             //!< Y Production vertex
        static double z( int index );             //!< Z Production vertex
        static double t( int index );             //!< production time

	////////////////////
	// Set Methods    //
	////////////////////

	/// set event number
        static void set_event_number( int evtno );
	/// set number of entries in HEPEVT
        static void set_number_entries( int noentries );
	/// set particle status
        static void set_status( int index, int status );
	/// set particle ID
        static void set_id( int index, int id );
	/// define parents of a particle
        static void set_parents( int index, int firstparent, int lastparent );
	/// define children of a particle
        static void set_children( int index, int firstchild, int lastchild );
	/// set particle momentum
        static void set_momentum( int index, double px, double py,
				  double pz, double e );
	/// set particle mass
        static void set_mass( int index, double mass );
	/// set particle production vertex
        static void set_position( int index, double x, double y, double z, 
				  double t );
	//////////////////////
	// HEPEVT Floorplan //
	//////////////////////
	static unsigned int sizeof_int();  //!< size of integer in bytes
	static unsigned int sizeof_real(); //!< size of real in bytes
        static int  max_number_entries();  //!< size of common block
	static void set_sizeof_int(unsigned int);  //!< define size of integer
	static void set_sizeof_real(unsigned int); //!< define size of real
	static void set_max_number_entries(unsigned int); //!< define size of common block

    protected:
        /// navigate a byte array
	static double byte_num_to_double( unsigned int );
        /// navigate a byte array
	static int    byte_num_to_int( unsigned int );
        /// pretend common block is an array of bytes
	static void   write_byte_num( double, unsigned int );
        /// pretend common block is an array of bytes
	static void   write_byte_num( int, unsigned int );
	/// print output legend
	static void   print_legend( std::ostream& ostr = std::cout );

    private:
	static unsigned int s_sizeof_int;
	static unsigned int s_sizeof_real;
	static unsigned int s_max_number_entries;

    }; 

    //////////////////////////////
    // HEPEVT Floorplan Inlines //
    //////////////////////////////
    inline unsigned int HEPEVT_Wrapper::sizeof_int(){ return s_sizeof_int; }

    inline unsigned int HEPEVT_Wrapper::sizeof_real(){ return s_sizeof_real; }

    inline int HEPEVT_Wrapper::max_number_entries() 
    { return (int)s_max_number_entries; }

    inline void HEPEVT_Wrapper::set_sizeof_int( unsigned int size ) 
    {
	if ( size != sizeof(short int) && size != sizeof(long int) && size != sizeof(int) ) {
	    std::cerr << "HepMC is not able to handle integers "
		      << " of size other than 2 or 4."
		      << " You requested: " << size << std::endl;
	}
	s_sizeof_int = size;
    }

    inline void HEPEVT_Wrapper::set_sizeof_real( unsigned int size ) {
	if ( size != sizeof(float) && size != sizeof(double) ) {
	    std::cerr << "HepMC is not able to handle floating point numbers"
		      << " of size other than 4 or 8."
		      << " You requested: " << size << std::endl;
	}
	s_sizeof_real = size;
    }

    inline void HEPEVT_Wrapper::set_max_number_entries( unsigned int size ) {
	s_max_number_entries = size;
    }

    inline double HEPEVT_Wrapper::byte_num_to_double( unsigned int b ) {
	if ( b >= hepevt_bytes_allocation ) std::cerr 
		  << "HEPEVT_Wrapper: requested hepevt data exceeds allocation"
		  << std::endl;
	if ( s_sizeof_real == sizeof(float) ) {
	    float* myfloat = (float*)&hepevt.data[b];
	    return (double)(*myfloat);
	} else if ( s_sizeof_real == sizeof(double) ) {
	    double* mydouble = (double*)&hepevt.data[b];
	    return (*mydouble);
	} else {
	    std::cerr 
		<< "HEPEVT_Wrapper: illegal floating point number length." 
		<< s_sizeof_real << std::endl;
	}
	return 0;
    }

    inline int HEPEVT_Wrapper::byte_num_to_int( unsigned int b ) {
	if ( b >= hepevt_bytes_allocation ) std::cerr 
		  << "HEPEVT_Wrapper: requested hepevt data exceeds allocation"
		  << std::endl;
	if ( s_sizeof_int == sizeof(short int) ) {
	    short int* myshortint = (short int*)&hepevt.data[b];
	    return (int)(*myshortint);
	} else if ( s_sizeof_int == sizeof(long int) ) {
	    long int* mylongint = (long int*)&hepevt.data[b];
	    return (*mylongint);
       // on some 64 bit machines, int, short, and long are all different
	} else if ( s_sizeof_int == sizeof(int) ) {
	    int* myint = (int*)&hepevt.data[b];
	    return (*myint);
	} else {
	    std::cerr 
		<< "HEPEVT_Wrapper: illegal integer number length." 
		<< s_sizeof_int << std::endl;
	}
	return 0;
    }

    inline void HEPEVT_Wrapper::write_byte_num( double in, unsigned int b ) {
	if ( b >= hepevt_bytes_allocation ) std::cerr 
		  << "HEPEVT_Wrapper: requested hepevt data exceeds allocation"
		  << std::endl;
	if ( s_sizeof_real == sizeof(float) ) {
	    float* myfloat = (float*)&hepevt.data[b];
	    (*myfloat) = (float)in;
	} else if ( s_sizeof_real == sizeof(double) ) {
	    double* mydouble = (double*)&hepevt.data[b];
	    (*mydouble) = (double)in;
	} else {
	    std::cerr 
		<< "HEPEVT_Wrapper: illegal floating point number length." 
		<< s_sizeof_real << std::endl;
	}
    }

    inline void HEPEVT_Wrapper::write_byte_num( int in, unsigned int b ) {
	if ( b >= hepevt_bytes_allocation ) std::cerr 
		  << "HEPEVT_Wrapper: requested hepevt data exceeds allocation"
		  << std::endl;
	if ( s_sizeof_int == sizeof(short int) ) {
	    short int* myshortint = (short int*)&hepevt.data[b];
	    (*myshortint) = (short int)in;
	} else if ( s_sizeof_int == sizeof(long int) ) {
	    long int* mylongint = (long int*)&hepevt.data[b];
	    (*mylongint) = (int)in;
       // on some 64 bit machines, int, short, and long are all different
	} else if ( s_sizeof_int == sizeof(int) ) {
	    int* myint = (int*)&hepevt.data[b];
	    (*myint) = (int)in;
	} else {
	    std::cerr 
		<< "HEPEVT_Wrapper: illegal integer number length." 
		<< s_sizeof_int << std::endl;
	}
    }

    //////////////
    // INLINES  //
    //////////////

    inline bool HEPEVT_Wrapper::is_double_precision() 
    { 
	// true if 8byte floating point numbers are used in the HepEVT common.
	return ( sizeof(double) == sizeof_real() );
    }

    inline int HEPEVT_Wrapper::event_number()
    { return byte_num_to_int(0); }

    inline int HEPEVT_Wrapper::number_entries() 
    { 
	int nhep = byte_num_to_int( 1*sizeof_int() );
	return ( nhep <= max_number_entries() ?
		 nhep : max_number_entries() );
    }

    inline int HEPEVT_Wrapper::status( int index )   
    { return byte_num_to_int( (2+index-1) * sizeof_int() ); }

    inline int HEPEVT_Wrapper::id( int index )
    { 
	return byte_num_to_int( (2+max_number_entries()+index-1) 
				* sizeof_int() ); 
    }

    inline int HEPEVT_Wrapper::first_parent( int index )
    { 
	int parent = byte_num_to_int( (2+2*max_number_entries()+2*(index-1)) 
				      * sizeof_int() ); 
	return ( parent > 0 && parent <= number_entries() ) ?
					 parent : 0; 
    }

    inline int HEPEVT_Wrapper::last_parent( int index )
    { 
	// Returns the Index of the LAST parent in the HEPEVT record
	// for particle with Index index.
	// If there is only one parent, the last parent is forced to 
	// be the same as the first parent.
	// If there are no parents for this particle, both the first_parent
	// and the last_parent with return 0.
	// Error checking is done to ensure the parent is always
	// within range ( 0 <= parent <= nhep )
	//
	int firstparent = first_parent(index);
	int parent = byte_num_to_int( (2+2*max_number_entries()+2*(index-1)+1) 
				      * sizeof_int() ); 
	return ( parent > firstparent && parent <= number_entries() ) 
						   ? parent : firstparent; 
    }

    inline int HEPEVT_Wrapper::number_parents( int index ) {
	int firstparent = first_parent(index);
	return ( firstparent>0 ) ? 
	    ( 1+last_parent(index)-firstparent ) : 0;
    }

    inline int HEPEVT_Wrapper::first_child( int index )
    { 
	int child = byte_num_to_int( (2+4*max_number_entries()+2*(index-1)) 
				     * sizeof_int() ); 
	return ( child > 0 && child <= number_entries() ) ?
				       child : 0; 
    }

    inline int HEPEVT_Wrapper::last_child( int index )
    { 
	// Returns the Index of the LAST child in the HEPEVT record
	// for particle with Index index.
	// If there is only one child, the last child is forced to 
	// be the same as the first child.
	// If there are no children for this particle, both the first_child
	// and the last_child with return 0.
	// Error checking is done to ensure the child is always
	// within range ( 0 <= parent <= nhep )
	//
	int firstchild = first_child(index);
	int child = byte_num_to_int( (2+4*max_number_entries()+2*(index-1)+1) 
				     * sizeof_int() ); 
	return ( child > firstchild && child <= number_entries() ) 
						? child : firstchild;
    }

    inline int HEPEVT_Wrapper::number_children( int index ) 
    {
	int firstchild = first_child(index);
	return ( firstchild>0 ) ? 
	    ( 1+last_child(index)-firstchild ) : 0;
    }

    inline double HEPEVT_Wrapper::px( int index )
    { 
	return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
				 + (5*(index-1)+0) *sizeof_real() );
    }

    inline double HEPEVT_Wrapper::py( int index )
    { 
	return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
				 + (5*(index-1)+1) *sizeof_real() );
    }


    inline double HEPEVT_Wrapper::pz( int index )
    { 
	return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
				 + (5*(index-1)+2) *sizeof_real() );
    }

    inline double HEPEVT_Wrapper::e( int index )
    { 
	return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
				 + (5*(index-1)+3) *sizeof_real() );
    }

    inline double HEPEVT_Wrapper::m( int index )
    { 
	return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
				 + (5*(index-1)+4) *sizeof_real() );
    }

    inline double HEPEVT_Wrapper::x( int index )
    { 
	return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
				   + ( 5*max_number_entries()
				       + (4*(index-1)+0) ) *sizeof_real() );
    }

    inline double HEPEVT_Wrapper::y( int index )
    { 
	return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
				   + ( 5*max_number_entries()
				       + (4*(index-1)+1) ) *sizeof_real() );
    }

    inline double HEPEVT_Wrapper::z( int index )
    { 
	return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
				   + ( 5*max_number_entries()
				       + (4*(index-1)+2) ) *sizeof_real() );
    }

    inline double HEPEVT_Wrapper::t( int index )
    { 
	return byte_num_to_double( (2+6*max_number_entries())*sizeof_int()
				   + ( 5*max_number_entries()
				       + (4*(index-1)+3) ) *sizeof_real() );
    }

    inline void HEPEVT_Wrapper::set_event_number( int evtno ) 
    { write_byte_num( evtno, 0 ); }

    inline void HEPEVT_Wrapper::set_number_entries( int noentries ) 
    { write_byte_num( noentries, 1*sizeof_int() ); }

    inline void HEPEVT_Wrapper::set_status( int index, int status ) 
    {
        if ( index <= 0 || index > max_number_entries() ) return;
	write_byte_num( status, (2+index-1) * sizeof_int() );
    }

    inline void HEPEVT_Wrapper::set_id( int index, int id ) 
    {
        if ( index <= 0 || index > max_number_entries() ) return;
	write_byte_num( id, (2+max_number_entries()+index-1) *sizeof_int() );
    }

    inline void HEPEVT_Wrapper::set_parents( int index, int firstparent, 
					     int lastparent ) 
    {
        if ( index <= 0 || index > max_number_entries() ) return;
	write_byte_num( firstparent, (2+2*max_number_entries()+2*(index-1)) 
			             *sizeof_int() );
	write_byte_num( lastparent, (2+2*max_number_entries()+2*(index-1)+1) 
				    * sizeof_int() );
    }
    
    inline void HEPEVT_Wrapper::set_children( int index, int firstchild, 
					      int lastchild ) 
    {
        if ( index <= 0 || index > max_number_entries() ) return;
	write_byte_num( firstchild, (2+4*max_number_entries()+2*(index-1)) 
				     *sizeof_int() );
	write_byte_num( lastchild, (2+4*max_number_entries()+2*(index-1)+1) 
				    *sizeof_int() );
    }

    inline void HEPEVT_Wrapper::set_momentum( int index, double px, 
					      double py, double pz, double e ) 
    {
        if ( index <= 0 || index > max_number_entries() ) return;
	write_byte_num( px, (2+6*max_number_entries()) *sizeof_int()
			    + (5*(index-1)+0) *sizeof_real() );
	write_byte_num( py, (2+6*max_number_entries())*sizeof_int()
			    + (5*(index-1)+1) *sizeof_real() );
	write_byte_num( pz, (2+6*max_number_entries())*sizeof_int()
			    + (5*(index-1)+2) *sizeof_real() );
	write_byte_num( e,  (2+6*max_number_entries())*sizeof_int()
			    + (5*(index-1)+3) *sizeof_real() );
    }

    inline void HEPEVT_Wrapper::set_mass( int index, double mass ) 
    {
        if ( index <= 0 || index > max_number_entries() ) return;
	write_byte_num( mass, (2+6*max_number_entries())*sizeof_int()
			      + (5*(index-1)+4) *sizeof_real() );
    }

    inline void HEPEVT_Wrapper::set_position( int index, double x, double y,
					      double z, double t ) 
    {
        if ( index <= 0 || index > max_number_entries() ) return;
	write_byte_num( x, (2+6*max_number_entries())*sizeof_int()
			   + ( 5*max_number_entries()
			       + (4*(index-1)+0) ) *sizeof_real() );
	write_byte_num( y, (2+6*max_number_entries())*sizeof_int()
			   + ( 5*max_number_entries()
			       + (4*(index-1)+1) ) *sizeof_real() );
	write_byte_num( z, (2+6*max_number_entries())*sizeof_int()
			   + ( 5*max_number_entries()
			       + (4*(index-1)+2) ) *sizeof_real() );
	write_byte_num( t, (2+6*max_number_entries())*sizeof_int()
			   + ( 5*max_number_entries()
			       + (4*(index-1)+3) ) *sizeof_real() );
    }

} // HepMC

#endif  // HEPMC_HEPEVT_WRAPPER_H
//--------------------------------------------------------------------------

HepMC/TempParticleMap.h000644 000765 000024 00000004015 13310701701 015213 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HEPMC_TempParticleMap_H
#define HEPMC_TempParticleMap_H

//////////////////////////////////////////////////////////////////////////
// garren@fnal.gov, October 2007
//
// Used by IO classes
//////////////////////////////////////////////////////////////////////////

#include <map>

namespace HepMC {

    class GenParticle;

    //! TempParticleMap is a temporary GenParticle* container used during input.

    ///
    /// \class  TempParticleMap
    /// Used by IO classes for recoverable particle ordering.
    /// Map GenParticle* against both outgoing vertex and particle order.
    ///
    class TempParticleMap {
    public:
        typedef std::map<HepMC::GenParticle*,int> TempMap;
        typedef std::map<int,HepMC::GenParticle*> TempOrderMap;
	typedef TempMap::iterator     TempMapIterator;
	typedef TempOrderMap::iterator  orderIterator;
	
	TempParticleMap() 
	: m_particle_to_end_vertex(), m_particle_order() {}
	
	~TempParticleMap() {}
	
	TempMapIterator begin() { return m_particle_to_end_vertex.begin(); }
	TempMapIterator end() { return m_particle_to_end_vertex.end(); }
	orderIterator order_begin() { return m_particle_order.begin(); }
	orderIterator order_end() { return m_particle_order.end(); }
	
	int end_vertex( GenParticle* );
	
	void addEndParticle( GenParticle*, int& );
	
    private:
        TempMap       m_particle_to_end_vertex;
	TempOrderMap  m_particle_order;
    };
    
    inline int TempParticleMap::end_vertex( GenParticle* p )
    { 
        //return m_particle_to_end_vertex[p]->second; 
	TempMapIterator it = m_particle_to_end_vertex.find(p);
	if( it == end() ) return 0;
	return m_particle_to_end_vertex[p];
    }

    inline void TempParticleMap::addEndParticle( GenParticle* p, int& end_vtx_code )
    {
	m_particle_order[p->barcode()] = p;
        m_particle_to_end_vertex[p] = end_vtx_code;
    }        

} // HepMC

#endif  // HEPMC_TempParticleMap_H
//--------------------------------------------------------------------------
HepMC/IO_HEPEVT.h000644 000765 000024 00000013004 13310701701 013544 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HEPMC_IO_HEPEVT_H
#define HEPMC_IO_HEPEVT_H

//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, January 2000, refer to:
// M. Dobbs and J.B. Hansen, "The HepMC C++ Monte Carlo Event Record for
// High Energy Physics", Computer Physics Communications (to be published).
//
// HEPEVT IO class
//////////////////////////////////////////////////////////////////////////
//
// Important note: This class uses HepMC::HEPEVT_Wrapper which is an
//                 interface to the fortran77 HEPEVT common block.
//                 The precision and number of entries in the F77 common 
//                 block can be specified. See HepMC/HEPEVT_Wrapper.h.
//                 You will very likely have to specify these values for your
//                 application.
//
//

#include <map>
#include <vector>
#include "HepMC/IO_BaseClass.h"
#include "HepMC/HEPEVT_Wrapper.h"

namespace HepMC {

    class GenEvent;
    class GenVertex;
    class GenParticle;

    //! HEPEVT IO class

    ///
    /// \class  IO_HEPEVT
    /// IO class for reading the standard HEPEVT common block.
    ///
    class IO_HEPEVT : public IO_BaseClass {
    public:
	IO_HEPEVT();
	virtual           ~IO_HEPEVT();
	bool              fill_next_event( GenEvent* );
	void              write_event( const GenEvent* );
	void              print( std::ostream& ostr = std::cout ) const;
	
	// see comments below for these switches.
	/// default is false
	bool              trust_both_mothers_and_daughters() const;
	/// default is true
	bool              trust_mothers_before_daughters() const;
	/// default is true
	bool              print_inconsistency_errors() const;
	/// default is true
	bool              trust_beam_particles() const;
        /// define mother daughter trust rules
	void              set_trust_mothers_before_daughters( bool b = true );
        /// define mother daughter trust rules
	void              set_trust_both_mothers_and_daughters( bool b = false );
	/// Since HEPEVT has bi-directional pointers, it is possible that
	/// the mother/daughter pointers are inconsistent (though physically
	/// speaking this should never happen). In practise it happens often.
	/// When a conflict occurs (i.e. when mother/daughter pointers are in 
	/// disagreement, where an empty (0) pointer is not considered a 
	/// disagreement) an error is printed. These errors can be turned off 
	/// with:            myio_hepevt.set_print_inconsistency_errors(0);
	/// but it is STRONGLY recommended that you print the HEPEVT 
	/// common and understand the inconsistency BEFORE you turn off the
	/// errors. The messages are there for a reason [remember, there is
	/// no message printed when the information is missing, ... only when
	/// is it inconsistent. User beware.]
	/// You can inspect the HEPEVT common block for inconsistencies with
	///   HEPEVT_Wrapper::check_hepevt_consistency()
	///
	/// There is a switch controlling whether the mother pointers or
	/// the daughters are to be trusted.
	/// For example, in Pythia the mother information is always correctly
	/// included, but the daughter information is often left unfilled: in
	/// this case we want to trust the mother pointers and not necessarily
	/// the daughters. [THIS IS THE DEFAULT]. Unfortunately the reverse
	/// happens for the stdhep(2001) translation of Isajet, so we need
	/// an option to toggle the choices.
	void              set_print_inconsistency_errors( bool b = true );
        /// declare whether or not beam particles exist
	void              set_trust_beam_particles( bool b = true );

    protected: // for internal use only
        /// create a GenParticle
	GenParticle* build_particle( int index );
        /// create a production vertex
	void build_production_vertex( 
	    int i,std::vector<HepMC::GenParticle*>& hepevt_particle, GenEvent* evt );
        /// create an end vertex
	void build_end_vertex( 
	    int i, std::vector<HepMC::GenParticle*>& hepevt_particle, GenEvent* evt );
        /// find this particle in the particle map
	int  find_in_map( 
	    const std::map<HepMC::GenParticle*,int>& m, GenParticle* p) const;

    private: // use of copy constructor is not allowed
	IO_HEPEVT( const IO_HEPEVT& ) : IO_BaseClass() {}

    private: // data members

	bool m_trust_mothers_before_daughters;
	bool m_trust_both_mothers_and_daughters;
	bool m_print_inconsistency_errors; 
	bool m_trust_beam_particles;
    };

    ////////////////////////////
    // INLINES access methods //
    ////////////////////////////
    inline bool IO_HEPEVT::trust_both_mothers_and_daughters() const 
    { return m_trust_both_mothers_and_daughters; }
	
    inline bool IO_HEPEVT::trust_mothers_before_daughters() const 
    { return m_trust_mothers_before_daughters; }

    inline bool IO_HEPEVT::print_inconsistency_errors() const
    { return m_print_inconsistency_errors; }

    inline void IO_HEPEVT::set_trust_both_mothers_and_daughters( bool b )
    { m_trust_both_mothers_and_daughters = b; }

    inline void IO_HEPEVT::set_trust_mothers_before_daughters( bool b )
    { m_trust_mothers_before_daughters = b; }

    inline void IO_HEPEVT::set_print_inconsistency_errors( bool b  )
    { m_print_inconsistency_errors = b; }

    inline bool IO_HEPEVT::trust_beam_particles() const
    { return m_trust_beam_particles; }

    inline void IO_HEPEVT::set_trust_beam_particles( bool b )
    { m_trust_beam_particles = b; }

} // HepMC

#endif  // HEPMC_IO_HEPEVT_H
//--------------------------------------------------------------------------
HepMC/GenEvent.h000644 000765 000024 00000106660 13310701701 013710 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HEPMC_GEN_EVENT_H
#define HEPMC_GEN_EVENT_H

//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, September 1999, refer to:
// M. Dobbs and J.B. Hansen, "The HepMC C++ Monte Carlo Event Record for
// High Energy Physics", Computer Physics Communications (to be published).
//
// Event record for MC generators (for use at any stage of generation)
//////////////////////////////////////////////////////////////////////////
//
// This class is intended as both a "container class" ( to store a MC
//  event for interface between MC generators and detector simulation )
//  and also as a "work in progress class" ( that could be used inside
//  a generator and modified as the event is built ).
//
// Iterators are provided which allow the user to easily obtain a
//  list of particles or vertices in an event --- this list can be filled
//  subject to some sort of selection criteria. Examples are given below
//  ( see HepMC::copy_if and std::copy )

///
/// \namespace HepMC
/// All classes in the HepMC packages are in the HepMC namespace 
///
namespace HepMC {

    // To create a list from an iterator, use: (i.e. for a list of particles);
    // #include <algorithm>
    //     list<GenParticle*> thelist;
    //     copy( evt->particles_begin(), evt->particles_end(), 
    //           back_inserter(thelist) );
    // to create a list subject to a condition (predicate) use:
    //     list<GenParticle*> thelist;
    //     HepMC::copy_if( evt->particles_begin(), evt->particles_end(), 
    //                     back_inserter(thelist), is_photon() );
    // where is_photon() is a predicate like:
    //     class is_photon {
    //       public:
    //         bool operator() ( GenParticle const * p ) {
    //             if ( p && p->pdg_id() == 22 ) return true;
    //             return false;
    //         }
    //     };
    // which the user defines herself.

    /// define the type of iterator to use
    template <class InputIterator, class OutputIterator, class Predicate>
    void copy_if( InputIterator first, InputIterator last, OutputIterator out,
		  Predicate pred ) {
	for ( ; first != last; ++first ) { if ( pred(*first) ) out = *first; }
    }
} // HepMC

// Since a container of all vertices in the event is maintained, the time
//  required to loop over all vertices (or particles) is very fast -- and 
//  the user does not gain much by first making his own list.
//  (this is not true for the GenVertex:: versions of these iterators, which
//   allow you to specify the vertex starting point and range)

// Data Members:
// signal_process_id()   The integer ID that uniquely specifies this signal
//                       process, i.e. MSUB in Pythia. It is necessary to
//                       package this with each event rather than with the run
//                       because many processes may be generated within one
//                       run.
// event_number()        Strictly speaking we cannot think of any reason that
//                       an event would need to know its own event number, it
//                       is more likely something that would be assigned by
//                       a database. It is included anyway (tradition?) since
//                       we expect it may be useful for debugging. It can
//                       be reset later by a database.
// mpi()                 The number of multi parton interactions in the event.
//                       This is NOT beam pileup.  Set to -1 by default.
// beam_particles()      A pair of pointers to the incoming beam particles.
// signal_process_vertex() pointer to the vertex containing the signal process
// weights()             Vector of doubles which specify th weight of the evnt,
//                       the first entry will be the "event weight" used for
//                       hit and miss etc., but a general vector is used to
//                       allow for reweighting etc. We envision a list of
//                       WeightTags to be included with a run class which
//                       would specify the meaning of the Weights .
// random_states()       Vector of integers which specify the random number 
//                       generator's state for this event. It is left to the
//                       generator to make use of this. We envision a vector of
//                       RndmStatesTags to be included with a run class which
//                       would specify the meaning of the random_states.
//
///////////////////////
// Memory allocation //
///////////////////////
// -When a vertex (particle) is added to a event (vertex), it is "adopted" 
//  and becomes the responsibility of the event (vertex) to delete that 
//  particle. 
// -objects responsible for deleting memory:
//    -events delete included vertices
//    -each vertex deletes its outgoing particles which do not have decay
//     vertices
//    -each vertex deletes its incoming particles which do not
//     have creation vertices 
//
////////////////////////
// About the Barcodes //
////////////////////////
// - each vertex or particle has a barcode, which is just an integer which
//   uniquely identifies it inside the event (i.e. there is a one to one
//   mapping between particle memory addresses and particle barcodes... and 
//   the same applied for vertices)
// - The value of a barcode has NO MEANING and NO ORDER!
//   For the user's convenience, when an event is read in via an IO_method
//   from an indexed list (like the HEPEVT common block), then the index will
//   become the barcode for that particle.
// - particle barcodes are always positive integers
//   vertex barcodes are always negative integers
//   The barcodes are chosen and set automatically when a vertex or particle
//   comes under the ownership of an event (i.e. it is contained in an event).
// - You can tell when a particle or vertex is owned, because its 
//   parent_event() return value will return a pointer to the event which owns
//   it (or null if its an orphan).
// - Please note that the barcodes are intended for internal use within HepMC 
//   as a unique identifier for the particles and vertices.
//   Using the barcode to encode extra information is an abuse of 
//   the barcode data member and causes confusion among users. 
// 

#include "HepMC/GenVertex.h"
#include "HepMC/GenParticle.h"
#include "HepMC/WeightContainer.h"
#include "HepMC/GenCrossSection.h"
#include "HepMC/HeavyIon.h"
#include "HepMC/PdfInfo.h"
#include "HepMC/Units.h"
#include "HepMC/HepMCDefs.h"
#include <map>
#include <string>
#include <vector>
#include <algorithm>
#include <iostream>

namespace HepMC {

    class GenEventVertexRange;
    class ConstGenEventVertexRange;
    class GenEventParticleRange;
    class ConstGenEventParticleRange;

    //! The GenEvent class is the core of HepMC

    ///
    /// \class GenEvent 
    /// HepMC::GenEvent contains information about generated particles.
    /// GenEvent is structured as a set of vertices which contain the particles.
    ///
    class GenEvent {
	friend class GenParticle;
	friend class GenVertex;  
    public:
        /// default constructor creates null pointers to HeavyIon, PdfInfo, and GenCrossSection
	GenEvent( int signal_process_id = 0, int event_number = 0,
		  GenVertex* signal_vertex = 0,
		  const WeightContainer& weights = std::vector<double>(),
		  const std::vector<long>& randomstates = std::vector<long>(),
		  Units::MomentumUnit = Units::default_momentum_unit(), 
		  Units::LengthUnit = Units::default_length_unit() );
        /// explicit constructor that takes HeavyIon and PdfInfo
	GenEvent( int signal_process_id, int event_number,
		  GenVertex* signal_vertex, const WeightContainer& weights,
		  const std::vector<long>& randomstates,
		  const HeavyIon& ion, const PdfInfo& pdf,
		  Units::MomentumUnit = Units::default_momentum_unit(), 
		  Units::LengthUnit = Units::default_length_unit() );
        /// constructor requiring units - all else is default
	GenEvent( Units::MomentumUnit, Units::LengthUnit,
	          int signal_process_id = 0, int event_number = 0,
		  GenVertex* signal_vertex = 0,
		  const WeightContainer& weights = std::vector<double>(),
		  const std::vector<long>& randomstates = std::vector<long>() );
        /// explicit constructor with units first that takes HeavyIon and PdfInfo
	GenEvent( Units::MomentumUnit, Units::LengthUnit,
	          int signal_process_id, int event_number,
		  GenVertex* signal_vertex, const WeightContainer& weights,
		  const std::vector<long>& randomstates,
		  const HeavyIon& ion, const PdfInfo& pdf );
	GenEvent( const GenEvent& inevent );          //!< deep copy
	GenEvent& operator=( const GenEvent& inevent ); //!< make a deep copy
	virtual ~GenEvent(); //!<deletes all vertices/particles in this evt

        void swap( GenEvent & other );  //!< swap
    
	void print( std::ostream& ostr = std::cout ) const; //!< dumps to ostr
	void print_version( std::ostream& ostr = std::cout ) const; //!< dumps release version to ostr

        /// assign a barcode to a particle
	GenParticle* barcode_to_particle( int barCode ) const;
        /// assign a barcode to a vertex
	GenVertex*   barcode_to_vertex(   int barCode ) const;

	////////////////////
	// access methods //
	////////////////////

	int signal_process_id() const; //!<  unique signal process id
	int event_number() const; //!<  event number
	int mpi() const;          //!<  number of multi parton interactions
	double event_scale() const; //!< energy scale, see hep-ph/0109068
	double alphaQCD() const; //!<  QCD coupling, see hep-ph/0109068
	double alphaQED() const; //!<  QED coupling, see hep-ph/0109068
        /// pointer to the vertex containing the signal process
	GenVertex* signal_process_vertex() const;
	/// test to see if we have two valid beam particles
	bool valid_beam_particles() const;
	/// pair of pointers to the two incoming beam particles
	std::pair<HepMC::GenParticle*,HepMC::GenParticle*> beam_particles() const;
	/// check GenEvent for validity
	/// A GenEvent is presumed valid if it has particles and/or vertices.
	bool is_valid() const;

	/// direct access to the weights container is allowed. 
	/// Thus you can use myevt.weights()[2];
	/// to access element 2 of the weights.
	/// or use myevt.weights().push_back( mywgt ); to add an element.
	/// and you can set the weights with myevt.weights() = myvector;
	WeightContainer&        weights(); //!< direct access to WeightContainer
	const WeightContainer&  weights() const; //!< direct access to WeightContainer

	/// access the GenCrossSection container if it exists
	GenCrossSection const *     cross_section() const;
	GenCrossSection*            cross_section();
	/// access the HeavyIon container if it exists
	HeavyIon const *         heavy_ion() const;
	HeavyIon*                heavy_ion();
	/// access the PdfInfo container if it exists
	PdfInfo const *          pdf_info() const;
	PdfInfo*                 pdf_info();

	/// vector of integers containing information about the random state
	const std::vector<long>& random_states() const;

        /// how many particle barcodes exist?
	int     particles_size() const;
        /// return true if there are no particle barcodes
	bool    particles_empty() const;
        /// how many vertex barcodes exist?
	int     vertices_size() const;
        /// return true if there are no vertex barcodes
	bool    vertices_empty() const;

	/// Write the unit information to an output stream.  
	/// If the output stream is not defined, use std::cout.
        void write_units( std::ostream & os = std::cout ) const; 
	/// If the cross section is defined,
	/// write the cross section information to an output stream.  
	/// If the output stream is not defined, use std::cout.
	void write_cross_section( std::ostream& ostr = std::cout ) const;

	/// Units used by the GenParticle momentum FourVector.
	Units::MomentumUnit momentum_unit() const;
	/// Units used by the GenVertex position FourVector.
	Units::LengthUnit   length_unit()   const;
	
	std::ostream& write(std::ostream&);
	std::istream& read(std::istream&);

	/////////////////////
	// mutator methods //
	/////////////////////

	bool    add_vertex( GenVertex* vtx );    //!< adds to evt and adopts
	bool    remove_vertex( GenVertex* vtx ); //!< erases vtx from evt
	void    clear();                         //!< empties the entire event
	
	void set_signal_process_id( int id ); //!< set unique signal process id
	void set_event_number( int eventno ); //!< set event number
	void set_mpi( int  ); //!< set number of multi parton interactions
	void set_event_scale( double scale ); //!< set energy scale
	void set_alphaQCD( double a ); //!< set QCD coupling
	void set_alphaQED( double a ); //!< set QED coupling

        /// set pointer to the vertex containing the signal process
	void set_signal_process_vertex( GenVertex* );
	/// set incoming beam particles
	bool set_beam_particles(GenParticle*, GenParticle*);
        /// use a pair of GenParticle*'s to set incoming beam particles
	bool set_beam_particles(std::pair<HepMC::GenParticle*,HepMC::GenParticle*> const &);
	/// provide random state information
	void set_random_states( const std::vector<long>& randomstates );

	/// provide a pointer to the GenCrossSection container
	void set_cross_section( const GenCrossSection& );
	/// provide a pointer to the HeavyIon container
	void set_heavy_ion( const HeavyIon& ion );
	/// provide a pointer to the PdfInfo container
	void set_pdf_info( const PdfInfo& p );
	
	/// set the units using enums
	/// This method will convert momentum and position data if necessary
	void use_units( Units::MomentumUnit, Units::LengthUnit );
	/// set the units using strings
	/// the string must match the enum exactly
	/// This method will convert momentum and position data if necessary
        void use_units( std::string&, std::string& );
	
	/// set the units using enums
	/// This method will NOT convert momentum and position data
	void define_units( Units::MomentumUnit, Units::LengthUnit );
	/// set the units using strings
	/// the string must match the enum exactly
	/// This method will NOT convert momentum and position data
        void define_units( std::string&, std::string& );
	
	/// vertex range
	GenEventVertexRange vertex_range();
	/// vertex range
	ConstGenEventVertexRange vertex_range() const;
	/// particle range
	GenEventParticleRange particle_range();
	/// particle range
	ConstGenEventParticleRange particle_range() const;

    public:
	///////////////////////////////
	// vertex_iterators          //
	///////////////////////////////
	// Note:  the XXX_iterator is "resolvable" as XXX_const_iterator, but 
	//  not the reverse, which is consistent with STL, 
	//  see Musser, Derge, Saini 2ndEd. p. 69,70.

	//!  const vertex iterator

	/// \class  vertex_const_iterator
	/// HepMC::GenEvent::vertex_const_iterator
	/// is used to iterate over all vertices in the event.
	class vertex_const_iterator :
	  public std::iterator<std::forward_iterator_tag,HepMC::GenVertex*,ptrdiff_t>{
	    // Iterates over all vertices in this event
	public:
	    /// constructor requiring vertex information
	    vertex_const_iterator(
		const 
		std::map<int,HepMC::GenVertex*,std::greater<int> >::const_iterator& i)
		: m_map_iterator(i) {}
	    vertex_const_iterator() {}
	    /// copy constructor
	    vertex_const_iterator( const vertex_const_iterator& i )
		{ *this = i; }
	    virtual ~vertex_const_iterator() {}
	    /// make a copy
	    vertex_const_iterator&  operator=( const vertex_const_iterator& i )
		{ m_map_iterator = i.m_map_iterator; return *this; }
	    /// return a pointer to a GenVertex
	    GenVertex* operator*(void) const { return m_map_iterator->second; }
	    /// Pre-fix increment
	    vertex_const_iterator&  operator++(void)  //Pre-fix increment 
		{ ++m_map_iterator; return *this; }
	    /// Post-fix increment
	    vertex_const_iterator   operator++(int)   //Post-fix increment
		{ vertex_const_iterator out(*this); ++(*this); return out; }
	    /// equality
	    bool  operator==( const vertex_const_iterator& a ) const
		{ return m_map_iterator == a.m_map_iterator; }
	    /// inequality
	    bool  operator!=( const vertex_const_iterator& a ) const
		{ return !(m_map_iterator == a.m_map_iterator); }
	protected:
	    /// const iterator to a vertex map
	    std::map<int,HepMC::GenVertex*,std::greater<int> >::const_iterator 
	                                                        m_map_iterator;
	private:
	    /// Pre-fix increment -- is not allowed
	    vertex_const_iterator&  operator--(void);
	    /// Post-fix increment -- is not allowed
	    vertex_const_iterator   operator--(int);
	};
	friend class vertex_const_iterator;
	/// begin vertex iteration
	vertex_const_iterator      vertices_begin() const
	    { return GenEvent::vertex_const_iterator( 
		m_vertex_barcodes.begin() ); }
	/// end vertex iteration
	vertex_const_iterator      vertices_end() const
	    { return GenEvent::vertex_const_iterator(
		m_vertex_barcodes.end() ); }


	//!  non-const vertex iterator

	/// \class  vertex_iterator
	/// HepMC::GenEvent::vertex_iterator
	/// is used to iterate over all vertices in the event.
	class vertex_iterator :
	  public std::iterator<std::forward_iterator_tag,HepMC::GenVertex*,ptrdiff_t>{
	    // Iterates over all vertices in this event
	public:
	    /// constructor requiring vertex information
	    vertex_iterator( 
		const 
		std::map<int,HepMC::GenVertex*,std::greater<int> >::iterator& i )
		: m_map_iterator( i ) {}
	    vertex_iterator() {}
	    /// copy constructor
	    vertex_iterator( const vertex_iterator& i ) { *this = i; }
	    virtual ~vertex_iterator() {}
	    /// make a copy
	    vertex_iterator&  operator=( const vertex_iterator& i ) {
		m_map_iterator = i.m_map_iterator;
		return *this;
	    }
	    /// const vertex iterator
	    operator vertex_const_iterator() const
		{ return vertex_const_iterator(m_map_iterator); }
	    /// return a pointer to a GenVertex
	    GenVertex*        operator*(void) const
		{ return m_map_iterator->second; }
	    /// Pre-fix increment
	    vertex_iterator&  operator++(void)  //Pre-fix increment 
		{ ++m_map_iterator;     return *this; }
	    /// Post-fix increment
	    vertex_iterator   operator++(int)   //Post-fix increment
		{ vertex_iterator out(*this); ++(*this); return out; }
	    /// equality
	    bool              operator==( const vertex_iterator& a ) const
		{ return m_map_iterator == a.m_map_iterator; }
	    /// inequality
	    bool              operator!=( const vertex_iterator& a ) const
		{ return !(m_map_iterator == a.m_map_iterator); }
	protected:
	    /// iterator to the vertex map
	    std::map<int,HepMC::GenVertex*,std::greater<int> >::iterator 
	                                                       m_map_iterator;
	private:
	    /// Pre-fix increment
	    vertex_iterator&  operator--(void);
	    /// Post-fix increment
	    vertex_iterator   operator--(int);

	};
	friend class vertex_iterator;
	/// begin vertex iteration
	vertex_iterator            vertices_begin() 
	    { return GenEvent::vertex_iterator( 
		m_vertex_barcodes.begin() ); }
	/// end vertex iteration
        vertex_iterator            vertices_end()
	    { return GenEvent::vertex_iterator(
		m_vertex_barcodes.end() ); }

    public:
	///////////////////////////////
	// particle_iterator         //
	///////////////////////////////
	// Example of iterating over all particles in the event:
	//      for ( GenEvent::particle_const_iterator p = particles_begin();
	//            p != particles_end(); ++p ) {
	//         (*p)->print();
	//      }
	//

	//!  const particle iterator

	/// \class  particle_const_iterator
	/// HepMC::GenEvent::particle_const_iterator 
	/// is used to iterate over all particles in the event.
	class particle_const_iterator :
	  public std::iterator<std::forward_iterator_tag,HepMC::GenParticle*,ptrdiff_t>{
	    // Iterates over all vertices in this event
	public:
	    /// iterate over particles
	    particle_const_iterator(
		const std::map<int,HepMC::GenParticle*>::const_iterator& i )
		: m_map_iterator(i) {}
	    particle_const_iterator() {}
	    /// copy constructor
	    particle_const_iterator( const particle_const_iterator& i )
		{ *this = i; }
	    virtual ~particle_const_iterator() {}
	    /// make a copy
	    particle_const_iterator& operator=(
		const particle_const_iterator& i )
		{ m_map_iterator = i.m_map_iterator; return *this; }
	    /// return a pointer to GenParticle
	    GenParticle*        operator*(void) const
		{ return m_map_iterator->second; }
	    /// Pre-fix increment
	    particle_const_iterator&  operator++(void)  //Pre-fix increment 
		{ ++m_map_iterator; return *this; }
	    /// Post-fix increment
	    particle_const_iterator   operator++(int)   //Post-fix increment
		{ particle_const_iterator out(*this); ++(*this); return out; }
	    /// equality
	    bool  operator==( const particle_const_iterator& a ) const
		{ return m_map_iterator == a.m_map_iterator; }
	    /// inequality
	    bool  operator!=( const particle_const_iterator& a ) const
		{ return !(m_map_iterator == a.m_map_iterator); }
	protected:
	    /// const iterator to the GenParticle map
	    std::map<int,HepMC::GenParticle*>::const_iterator m_map_iterator;
	private:
	    /// Pre-fix increment
	    particle_const_iterator&  operator--(void);
	    /// Post-fix increment
	    particle_const_iterator   operator--(int);
	};	
	friend class particle_const_iterator;
	/// begin particle iteration
	particle_const_iterator      particles_begin() const
	    { return GenEvent::particle_const_iterator( 
		m_particle_barcodes.begin() ); }
	/// end particle iteration
	particle_const_iterator      particles_end() const
	    { return GenEvent::particle_const_iterator(
		m_particle_barcodes.end() ); }

	//!  non-const particle iterator

	/// \class  particle_iterator
	/// HepMC::GenEvent::particle_iterator 
	/// is used to iterate over all particles in the event.
 	class particle_iterator :
	  public std::iterator<std::forward_iterator_tag,HepMC::GenParticle*,ptrdiff_t>{
	    // Iterates over all vertices in this event
	public:
	    /// iterate over particles
	    particle_iterator( const std::map<int,HepMC::GenParticle*>::iterator& i )
		: m_map_iterator( i ) {}
	    particle_iterator() {}
	    /// copy constructor
	    particle_iterator( const particle_iterator& i ) { *this = i; }
	    virtual ~particle_iterator() {}
	    /// make a copy
	    particle_iterator&  operator=( const particle_iterator& i ) {
		m_map_iterator = i.m_map_iterator;
		return *this;
	    }
	    /// const particle iterator
	    operator particle_const_iterator() const
		{ return particle_const_iterator(m_map_iterator); }
	    /// return pointer to GenParticle
	    GenParticle*        operator*(void) const
		{ return m_map_iterator->second; }
            /// Pre-fix increment
	    particle_iterator&  operator++(void) 
		{ ++m_map_iterator;     return *this; }
            /// Post-fix increment
	    particle_iterator   operator++(int)   
		{ particle_iterator out(*this); ++(*this); return out; }
            /// equality
	    bool              operator==( const particle_iterator& a ) const
		{ return m_map_iterator == a.m_map_iterator; }
            /// inequality
	    bool              operator!=( const particle_iterator& a ) const
		{ return !(m_map_iterator == a.m_map_iterator); }
	protected:
	    /// iterator for GenParticle map
	    std::map<int,HepMC::GenParticle*>::iterator m_map_iterator;
	private:
            /// Pre-fix increment
	    particle_iterator&  operator--(void);
            /// Post-fix increment
	    particle_iterator   operator--(int);
	};
	friend class particle_iterator;
	/// begin particle iteration
	particle_iterator particles_begin() 
	    { return GenEvent::particle_iterator(
		m_particle_barcodes.begin() ); }
	/// end particle iteration
        particle_iterator particles_end()
	    { return GenEvent::particle_iterator(
		m_particle_barcodes.end() ); }

	////////////////////////////////////////////////
    protected:
	//
	// Following methods intended for use by GenParticle/Vertex classes:
	// In general there is no reason they should be used elsewhere.
	/// set the barcode - intended for use by GenParticle
	bool         set_barcode( GenParticle* p, int suggested_barcode =false );
	/// set the barcode - intended for use by GenVertex
	bool         set_barcode( GenVertex*   v, int suggested_barcode =false );
	///  intended for use by GenParticle
	void         remove_barcode( GenParticle* p );
	///  intended for use by GenVertex
	void         remove_barcode( GenVertex*   v );

   	void delete_all_vertices(); //!<delete all vertices owned by this event

     private: // methods
        /// internal method used when converting momentum units
        bool use_momentum_unit( Units::MomentumUnit );
        bool use_momentum_unit( std::string& );
        /// internal method used when converting length units
	bool use_length_unit( Units::LengthUnit );
	bool use_length_unit( std::string& );
	
	// the following internal methods are used by read() and write()

	/// send the beam particles to ASCII output
	std::ostream & write_beam_particles( std::ostream &, 
                	     std::pair<HepMC::GenParticle *,HepMC::GenParticle *> );
	/// send a GenVertex to ASCII output
	std::ostream & write_vertex( std::ostream &, GenVertex const * );
	/// send a GenParticle to ASCII output
	std::ostream & write_particle( std::ostream&, GenParticle const * );
	/// find the file type
	std::istream & find_file_type( std::istream & );
	/// find the key at the end of the block
	std::istream & find_end_key( std::istream &, int & );
        /// get unit information from ASCII input
        std::istream & read_units( std::istream & );
	/// get weight names from ASCII input
        std::istream & read_weight_names( std::istream & );
	/// read the event header line
        std::istream & process_event_line( std::istream &, int &, int &, int &, int & );

    private: // data members
	int                   m_signal_process_id;
	int                   m_event_number;  
	int                   m_mpi;        // number of multi paricle interactions
	double                m_event_scale;// energy scale, see hep-ph/0109068
	double                m_alphaQCD;   // QCD coupling, see hep-ph/0109068
	double                m_alphaQED;   // QED coupling, see hep-ph/0109068
	GenVertex*            m_signal_process_vertex;
	GenParticle*          m_beam_particle_1;
	GenParticle*          m_beam_particle_2;
	WeightContainer       m_weights; // weights for this event first weight
	                                 // is used by default for hit and miss
	std::vector<long> m_random_states; // container of rndm num 
	                                       // generator states

	std::map< int,HepMC::GenVertex*,std::greater<int> >   m_vertex_barcodes;
	std::map< int,HepMC::GenParticle*,std::less<int> >    m_particle_barcodes;
	GenCrossSection*         m_cross_section; 	      // undefined by default
	HeavyIon*             m_heavy_ion; 	      // undefined by default
	PdfInfo*              m_pdf_info; 	      // undefined by default
	Units::MomentumUnit   m_momentum_unit;    // default value set by configure switch
	Units::LengthUnit     m_position_unit;    // default value set by configure switch

    };


    ///////////////////////////
    // IO Free Functions     //
    ///////////////////////////
  
    /// standard streaming IO output operator
    std::ostream & operator << (std::ostream &, GenEvent &);
    /// standard streaming IO input operator
    std::istream & operator >> (std::istream &, GenEvent &);
    /// set the units for this input stream
    std::istream & set_input_units(std::istream &, 
                                   Units::MomentumUnit, Units::LengthUnit);
    /// Explicitly write the begin block lines that IO_GenEvent uses
    std::ostream & write_HepMC_IO_block_begin(std::ostream & );
    /// Explicitly write the end block line that IO_GenEvent uses
    std::ostream & write_HepMC_IO_block_end(std::ostream & );


    ///////////////////////////
    // INLINE Free Functions //
    ///////////////////////////

    // Implemented in terms of GenEvent::use_...
    inline GenEvent& convert_units(GenEvent & evt, Units::MomentumUnit m, Units::LengthUnit l)
    {
      evt.use_units(m, l);
      return evt;
    }

    ///////////////////////////
    // INLINE Access Methods //
    ///////////////////////////

    ///  The integer ID that uniquely specifies this signal
    ///  process, i.e. MSUB in Pythia. It is necessary to
    ///  package this with each event rather than with the run
    ///  because many processes may be generated within one run.
    inline int GenEvent::signal_process_id() const 
    { return m_signal_process_id; }

    inline int GenEvent::event_number() const { return m_event_number; }

    /// Returns the number of multi parton interactions in the event.
    /// This number is -1 if it is not set.
    inline int GenEvent::mpi() const { return m_mpi; }

    inline double GenEvent::event_scale() const { return m_event_scale; }

    inline double GenEvent::alphaQCD() const { return m_alphaQCD; }

    inline double GenEvent::alphaQED() const { return m_alphaQED; }
 
    inline GenVertex* GenEvent::signal_process_vertex() const {
	/// returns a (mutable) pointer to the signal process vertex
	return m_signal_process_vertex;
    }  

    inline WeightContainer& GenEvent::weights() { return m_weights; }

    inline const WeightContainer& GenEvent::weights() const 
    { return m_weights; }

    inline GenCrossSection const * GenEvent::cross_section() const 
    { return m_cross_section; }

    inline GenCrossSection*  GenEvent::cross_section()  
    { return m_cross_section; }

    inline HeavyIon const * GenEvent::heavy_ion() const 
    { return m_heavy_ion; }

    inline HeavyIon*  GenEvent::heavy_ion()  
    { return m_heavy_ion; }

    inline PdfInfo const * GenEvent::pdf_info() const 
    { return m_pdf_info; }

    inline PdfInfo*  GenEvent::pdf_info()  
    { return m_pdf_info; }

    ///  Vector of integers which specify the random number 
    ///  generator's state for this event. It is left to the
    ///  generator to make use of this. We envision a vector of
    ///  RndmStatesTags to be included with a run class which
    ///  would specify the meaning of the random_states.
    inline const std::vector<long>& GenEvent::random_states() const 
    { return m_random_states; }

    inline void GenEvent::set_signal_process_id( int id )
    { m_signal_process_id = id; }

    inline void GenEvent::set_event_number( int eventno )
    { m_event_number = eventno; }

    /// Use this to set the number of multi parton interactions in each event.
    inline void GenEvent::set_mpi( int nmpi )
    { m_mpi = nmpi; }


    inline void GenEvent::set_event_scale( double sc ) { m_event_scale = sc; }

    inline void GenEvent::set_alphaQCD( double a ) { m_alphaQCD = a; }

    inline void GenEvent::set_alphaQED( double a ) { m_alphaQED = a; }

    inline void GenEvent::set_signal_process_vertex( GenVertex* vtx ) {
	m_signal_process_vertex = vtx;
	if ( m_signal_process_vertex ) add_vertex( m_signal_process_vertex );
    }

    inline void GenEvent::set_cross_section( const GenCrossSection& xs )
    { 
        delete m_cross_section;
        m_cross_section = new GenCrossSection(xs); 
    }

    inline void GenEvent::set_heavy_ion( const HeavyIon& ion )
    { 
        delete m_heavy_ion;
        m_heavy_ion = new HeavyIon(ion); 
    }

    inline void GenEvent::set_pdf_info( const PdfInfo& p )
    { 
        delete m_pdf_info;
        m_pdf_info = new PdfInfo(p); 
    }

    inline void GenEvent::set_random_states( const std::vector<long>&
					     randomstates )
    { m_random_states = randomstates; }

    inline void GenEvent::remove_barcode( GenParticle* p )
    { m_particle_barcodes.erase( p->barcode() ); }

    inline void GenEvent::remove_barcode( GenVertex* v )
    { m_vertex_barcodes.erase( v->barcode() ); }

    /// Each vertex or particle has a barcode, which is just an integer which
    /// uniquely identifies it inside the event (i.e. there is a one to one
    /// mapping between particle memory addresses and particle barcodes... and 
    /// the same applied for vertices).
    ///
    /// The value of a barcode has NO MEANING and NO ORDER!
    /// For the user's convenience, when an event is read in via an IO_method
    /// from an indexed list (like the HEPEVT common block), then the index will
    /// become the barcode for that particle.
    ///
    /// Particle barcodes are always positive integers.
    /// The barcodes are chosen and set automatically when a vertex or particle
    /// comes under the ownership of an event (i.e. it is contained in an event).
    /// 
    /// Please note that the barcodes are intended for internal use within 
    /// HepMC as a unique identifier for the particles and vertices.
    /// Using the barcode to encode extra information is an abuse of 
    /// the barcode data member and causes confusion among users. 
    inline GenParticle* GenEvent::barcode_to_particle( int barCode ) const
    { 
	std::map<int,HepMC::GenParticle*>::const_iterator i 
	    = m_particle_barcodes.find(barCode);
	return ( i != m_particle_barcodes.end() ) ? (*i).second : 0;
    }

    /// Each vertex or particle has a barcode, which is just an integer which
    /// uniquely identifies it inside the event (i.e. there is a one to one
    /// mapping between particle memory addresses and particle barcodes... and 
    /// the same applied for vertices).
    ///
    /// The value of a barcode has NO MEANING and NO ORDER!
    /// For the user's convenience, when an event is read in via an IO_method
    /// from an indexed list (like the HEPEVT common block), then the index will
    /// become the barcode for that particle.
    ///
    /// Vertex barcodes are always negative integers.
    /// The barcodes are chosen and set automatically when a vertex or particle
    /// comes under the ownership of an event (i.e. it is contained in an event).
    /// 
    /// Please note that the barcodes are intended for internal use within 
    /// HepMC as a unique identifier for the particles and vertices.
    /// Using the barcode to encode extra information is an abuse of 
    /// the barcode data member and causes confusion among users. 
    inline GenVertex* GenEvent::barcode_to_vertex( int barCode ) const
    {
	std::map<int,GenVertex*,std::greater<int> >::const_iterator i 
	    = m_vertex_barcodes.find(barCode);
	return ( i != m_vertex_barcodes.end() ) ? (*i).second : 0;
    }

    inline int GenEvent::particles_size() const {
	return (int)m_particle_barcodes.size();
    }
    inline bool GenEvent::particles_empty() const {
	return (bool)m_particle_barcodes.empty();
    }
    inline int GenEvent::vertices_size() const {
	return (int)m_vertex_barcodes.size();
    }
    inline bool GenEvent::vertices_empty() const {
	return (bool)m_vertex_barcodes.empty();
    }
    
    // beam particles
    inline std::pair<HepMC::GenParticle *,HepMC::GenParticle *> GenEvent::beam_particles() const {
        return std::pair<GenParticle *,GenParticle *> (m_beam_particle_1, m_beam_particle_2);
    }

    // units
    inline Units::MomentumUnit GenEvent::momentum_unit() const {
        return m_momentum_unit; 
    }
    inline Units::LengthUnit   GenEvent::length_unit()   const {
        return m_position_unit; 
    }
    
    inline void GenEvent::use_units( Units::MomentumUnit new_m, Units::LengthUnit new_l ) { 
       use_momentum_unit( new_m );
       use_length_unit( new_l );
    }
    
    inline void GenEvent::use_units( std::string& new_m, std::string& new_l ) { 
       use_momentum_unit( new_m );
       use_length_unit( new_l );
    }
    
    inline void GenEvent::define_units( Units::MomentumUnit new_m, Units::LengthUnit new_l ) { 
	m_momentum_unit = new_m; 
	m_position_unit = new_l; 
    }

} // HepMC

#endif  // HEPMC_GEN_EVENT_H

//--------------------------------------------------------------------------


HepMC/StreamInfo.h000644 000765 000024 00000012005 13310701701 014231 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HEPMC_STREAM_INFO_H
#define HEPMC_STREAM_INFO_H

//////////////////////////////////////////////////////////////////////////
// garren@fnal.gov, March 2009
//
// This class contains the extra information needed when using streaming IO
//////////////////////////////////////////////////////////////////////////

#include <string>
#include "HepMC/Units.h"

namespace HepMC {

/// The known_io enum is used to track which type of input is being read
enum known_io { gen=1, ascii, extascii, ascii_pdt, extascii_pdt };

//! StreamInfo contains extra information needed when using streaming IO.

///
/// \class  StreamInfo
/// This class contains the extra information needed when using streaming IO
/// to process HepMC GenEvents
///
class StreamInfo {
public:
    /// default constructor
    StreamInfo( );
    /// destructor
    ~StreamInfo() {}

    /// IO_GenEvent begin event block key
    std::string IO_GenEvent_Key()          const { return m_io_genevent_start; }
    /// IO_GenEvent end event block key
    std::string IO_GenEvent_End()          const { return m_io_genevent_end; }

    /// IO_Ascii begin event block key
    /// IO_Ascii has been removed, but we want to be able to read 
    /// existing files written by IO_Ascii
    std::string IO_Ascii_Key()             const { return m_io_ascii_start; }
    /// IO_Ascii end event block key
    std::string IO_Ascii_End()             const { return m_io_ascii_end; }
    /// IO_Ascii begin particle data block key
    std::string IO_Ascii_PDT_Key()             const { return m_io_ascii_pdt_start; }
    /// IO_Ascii end particle data block key
    std::string IO_Ascii_PDT_End()             const { return m_io_ascii_pdt_end; }

    /// IO_ExtendedAscii begin event block key
    /// IO_ExtendedAscii has been removed, but we want to be able to read 
    /// existing files written by IO_ExtendedAscii
    std::string IO_ExtendedAscii_Key()     const { return m_io_extendedascii_start; }
    /// IO_ExtendedAscii end event block key
    std::string IO_ExtendedAscii_End()     const { return m_io_extendedascii_end; }
    /// IO_ExtendedAscii begin particle data block key
    std::string IO_ExtendedAscii_PDT_Key()             const { return m_io_extendedascii_pdt_start; }
    /// IO_ExtendedAscii end particle data block key
    std::string IO_ExtendedAscii_PDT_End()             const { return m_io_extendedascii_pdt_end; }

    /// get IO type
    int io_type() const { return m_io_type; }
    /// set IO type
    void set_io_type( int );

    /// true if the stream has a file type key
    /// has_key is true by default
    bool has_key() const { return m_has_key; }
    /// set to false if the stream does not have a file type key
    void set_has_key( bool );
    
    /// get the I/O momentum units
    Units::MomentumUnit io_momentum_unit() const { return m_io_momentum_unit; }
    /// get the I/O length units
    Units::LengthUnit io_position_unit() const { return m_io_position_unit; }

    /// get the I/O stream id
    /// This is used for sanity checking.
    int stream_id() const { return m_stream_id; }
    
    /// Special information is processed the first time we use the IO
    bool finished_first_event() const { return m_finished_first_event_io; }
    /// Special information is processed the first time we use the IO
    void set_finished_first_event( bool b ) { m_finished_first_event_io = b; }

    /// needed when reading a file without units if those units are 
    /// different than the declared default units 
    /// (e.g., the default units are MeV, but the file was written with GeV)
    /// This method is not necessary if the units are written in the file
    void use_input_units( Units::MomentumUnit, Units::LengthUnit );
    
    /// reading_event_header will return true when streaming input is 
    /// processing the GenEvent header information
    bool reading_event_header();
    /// set the reading_event_header flag
    void set_reading_event_header(bool);

private: // data members
    bool        m_finished_first_event_io;
    // GenEvent I/O method keys
    std::string m_io_genevent_start;
    std::string m_io_ascii_start;
    std::string m_io_extendedascii_start;
    std::string m_io_genevent_end;
    std::string m_io_ascii_end;
    std::string m_io_extendedascii_end;
    // particle data I/O method keys
    std::string m_io_ascii_pdt_start;
    std::string m_io_extendedascii_pdt_start;
    std::string m_io_ascii_pdt_end;
    std::string m_io_extendedascii_pdt_end;
    // io information
    int         m_io_type;
    bool        m_has_key;
    // default io units - used only when reading a file with no units
    Units::MomentumUnit m_io_momentum_unit;
    Units::LengthUnit   m_io_position_unit;
    // used to keep identify the I/O stream
    unsigned int m_stream_id;
    static unsigned int m_stream_counter;
    // used to keep track when reading event
    bool m_reading_event_header;

};

} // HepMC

#endif  // HEPMC_STREAM_INFO_H
//--------------------------------------------------------------------------
HepMC/IO_HERWIG.h000644 000765 000024 00000015051 13310701701 013542 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HEPMC_IO_HERWIG_H
#define HEPMC_IO_HERWIG_H

//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, October 2002, refer to:
// M. Dobbs and J.B. Hansen, "The HepMC C++ Monte Carlo Event Record for
// High Energy Physics", Computer Physics Communications (to be published).
//
// IO class for reading the (non-standard) HEPEVT common block from 
//  the Herwig monte carlo program.
// Notes:
//   - The HERWIG HEPEVT common block is non-standard, primarily because it 
//     contains some color flow information. When you call IO_HERWIG, the 
//     HEPEVT common block is transformed to the standard. THIS CHANGES THE
//     CONTENT of HEPEVT!.
//   - The HERWIG HEPEVT common block has some EXTRA non-physical ENTRIES 
//     (such as CMS frame, HARD subprocess, and CONE).
//     These are removed by IO_HERWIG. Thus the HepMC event will APPEAR
//     to have fewer particles in it that herwig did.
//     There is a switch m_no_gaps_in_barcodes. For
//       true  - then the extra particles are removed from HEPEVT, with 
//               the result that the HepMC barcodes will be sequential, with 
//               no gaps.
//       false - the barcodes will correspond directly to the HEPEVT index, but
//               there will be gaps ... ie some barcodes will be unassigned.
//       this switch requested by I Hinchliffe, October 31, 2002
//   - some of the Herwig GLUON SPLITTING products are not properly documented
//     in hepevt. I was unable to repair this in a simple and robust way. 
//     Therefore some of the gluon splitting products will be orphans 
//     in the HepMC output. 
//   - Herwig uses      HEPEVT_Wrapper::set_max_number_entries(4000);
//  	                HEPEVT_Wrapper::set_sizeof_real(8);
//     which are the defaults for HEPEVT_Wrapper.
//////////////////////////////////////////////////////////////////////////
//

#include <set>
#include <vector>
#include "HepMC/IO_BaseClass.h"
#include "HepMC/HEPEVT_Wrapper.h"

namespace HepMC {

    class GenEvent;
    class GenVertex;
    class GenParticle;

    //! IO_HERWIG is used to get Herwig information

    ///
    /// \class  IO_HERWIG
    /// IO class for reading the HEPEVT common block from 
    ///  the Herwig monte carlo program.
    ///
    class IO_HERWIG : public IO_BaseClass {
    public:
	IO_HERWIG();
	virtual           ~IO_HERWIG();
        /// get the next event
	bool              fill_next_event( GenEvent* );
        /// write to ostr
	void              print( std::ostream& ostr = std::cout ) const;
        /// this information is dubious
	double            interfaces_to_version_number() const {return 6.400;}
		
	// see comments below for these switches.
	/// default is true
	bool              print_inconsistency_errors() const;
	/// decide whether or not to print inconsistency errors
	void              set_print_inconsistency_errors( bool b = true );

        /// ask how to deal with extra non-physical pseudo particles
	bool              no_gaps_in_barcodes() const 
	                     { return m_no_gaps_in_barcodes; }
	/// The HERWIG HEPEVT common block has some EXTRA non-physical ENTRIES 
	/// (such as CMS frame, HARD subprocess, and CONE).
	/// These are removed by IO_HERWIG. Thus the HepMC event will APPEAR
	/// to have fewer particles in it that herwig did.
	/// There is a switch m_no_gaps_in_barcodes. For
	///   true  - then the extra particles are removed from HEPEVT, with 
	///             the result that the HepMC barcodes will be sequential, with 
	///             no gaps.
	///   false - the barcodes will correspond directly to the HEPEVT index, but
	///             there will be gaps ... ie some barcodes will be unassigned.
	///   this switch requested by I Hinchliffe, October 31, 2002
	void              set_no_gaps_in_barcodes( bool a ) 
	                     { m_no_gaps_in_barcodes=a; }

    protected: // for internal use only
        /// default is true
	bool              trust_both_mothers_and_daughters() const;
        /// default is false
	bool              trust_mothers_before_daughters() const;
        /// define mother daughter trust rules
	void              set_trust_mothers_before_daughters( bool b = true );
        /// define mother daughter trust rules
	void              set_trust_both_mothers_and_daughters( bool b = false );

        /// make a particle
	GenParticle* build_particle( int index );
	/// make a production vertex
	void         build_production_vertex( 
	    int i,std::vector<GenParticle*>& hepevt_particle, GenEvent* evt );
	/// make a decay vertex
	void         build_end_vertex( 
	    int i, std::vector<GenParticle*>& hepevt_particle, GenEvent* evt );
        /// find this particle in the map
	int          find_in_map( 
	    const std::map<GenParticle*,int>& m, GenParticle* p) const;

        /// make the HERWIG HEPEVT common block look like the standard
	void repair_hepevt() const;
	/// deal with artifacts of repairing HEPEVT
	void remove_gaps_in_hepevt() const;
	/// zero out a HEPEVT pseudo particle
	void zero_hepevt_entry( int i ) const;
	/// translate particle ID
	int  translate_herwig_to_pdg_id( int i ) const;

    private: // following are not implemented for Herwig
	virtual void write_event( const GenEvent* ){}

    private: // use of copy constructor is not allowed
	IO_HERWIG( const IO_HERWIG& ) : IO_BaseClass() {}

    private: // data members
	bool              m_trust_mothers_before_daughters;
	bool              m_trust_both_mothers_and_daughters;
	bool              m_print_inconsistency_errors; 
	bool              m_no_gaps_in_barcodes;
	std::vector<int>  m_herwig_to_pdg_id;
	std::set<int>     m_no_antiparticles;
    };

    ////////////////////////////
    // INLINES access methods //
    ////////////////////////////
    inline bool IO_HERWIG::trust_both_mothers_and_daughters() const 
    { return m_trust_both_mothers_and_daughters; }
	
    inline bool IO_HERWIG::trust_mothers_before_daughters() const 
    { return m_trust_mothers_before_daughters; }

    inline bool IO_HERWIG::print_inconsistency_errors() const
    { return m_print_inconsistency_errors; }

    inline void IO_HERWIG::set_trust_both_mothers_and_daughters( bool b )
    { m_trust_both_mothers_and_daughters = b; }

    inline void IO_HERWIG::set_trust_mothers_before_daughters( bool b )
    { m_trust_mothers_before_daughters = b; }

    inline void IO_HERWIG::set_print_inconsistency_errors( bool b  )
    { m_print_inconsistency_errors = b; }

} // HepMC

#endif  // HEPMC_IO_HERWIG_H
//--------------------------------------------------------------------------
HepMC/StreamHelpers.h000644 000765 000024 00000004773 13310701701 014755 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HEPMC_STREAM_HELPERS_H
#define HEPMC_STREAM_HELPERS_H

//////////////////////////////////////////////////////////////////////////
// garren@fnal.gov, March 2009
//
// This header contains helper functions used by streaming IO
//////////////////////////////////////////////////////////////////////////

#include <ostream>
#include <istream>

#include "HepMC/GenEvent.h"
#include "HepMC/TempParticleMap.h"

namespace HepMC {

namespace detail {

/// used by IO_GenEvent constructor
std::ostream & establish_output_stream_info( std::ostream & );
/// used by IO_GenEvent constructor
std::istream & establish_input_stream_info( std::istream & );

/// get a GenVertex from ASCII input
/// TempParticleMap is used to track the associations of particles with vertices
std::istream & read_vertex( std::istream &, TempParticleMap &, GenVertex * );

/// get a GenParticle from ASCII input
/// TempParticleMap is used to track the associations of particles with vertices
std::istream & read_particle( std::istream&, TempParticleMap &, GenParticle * );

/// write a double - for internal use by streaming IO
inline std::ostream & output( std::ostream & os, const double& d ) {
    if( os  ) {
	if ( d == 0. ) {
	    os << ' ' << (int)0;
	} else {
	    os << ' ' << d;
	}
    }
    return os;
}

/// write a float - for internal use by streaming IO
inline std::ostream & output( std::ostream & os, const float& d ) {
    if( os  ) {
	if ( d == 0. ) {
	    os << ' ' << (int)0;
	} else {
	    os << ' ' << d;
	}
    }
    return os;
}

/// write an int - for internal use by streaming IO
inline std::ostream & output( std::ostream & os, const int& i ) { 
    if( os  ) {
	if ( i == 0. ) {
	    os << ' ' << (int)0;
	} else {
            os << ' ' << i; 
	}
    }
    return os;
}

/// write a long - for internal use by streaming IO
inline std::ostream & output( std::ostream & os, const long& i ) {
    if( os  ) {
	if ( i == 0. ) {
	    os << ' ' << (int)0;
	} else {
            os << ' ' << i; 
	}
    }
    return os;
}

/// write a single char - for internal use by streaming IO
inline std::ostream & output( std::ostream & os, const char& c ) {
    if( os  ) {
	if ( c ) {
            os << c; 
	} else {
	    os << ' ' ;
	}
    }
    return os;
}

/// used to read to the end of a bad event
std::istream & find_event_end( std::istream & );

} // detail

} // HepMC

#endif  // HEPMC_STREAM_HELPERS_H
//--------------------------------------------------------------------------
HepMC/Flow.h000644 000765 000024 00000017574 13310701701 013111 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HEPMC_FLOW_H
#define HEPMC_FLOW_H

//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, January 2000, refer to:
// M. Dobbs and J.B. Hansen, "The HepMC C++ Monte Carlo Event Record for
// High Energy Physics", Computer Physics Communications (to be published).
//
// particle's flow object
// keeps track of an arbitrary number of flow patterns within a graph 
// (i.e. color flow, charge flow, lepton number flow, ...) 
// Flow patterns are coded with an integer, in the same manner as in Herwig.
// Note: 0 is NOT allowed as code index nor as flow code since it
//       is used to indicate null.
//////////////////////////////////////////////////////////////////////////

// This class can be used to keep track of flow patterns within 
//  a graph. An example is color flow. If we have two quarks going through
//  an s-channel gluon to form two more quarks:
//
//  \q1       /q3   then we can keep track of the color flow with the
//   \_______/      HepMC::Flow class as follows: 
//   /   g   \. 
//  /q2       \q4
//
//  lets say the color flows from q2-->g-->q3  and q1-->g-->q4
//  the individual colors are unimportant, but the flow pattern is.
//  We can capture this flow by assigning the first pattern (q2-->g-->q3)
//  a unique (arbitrary) flow code 678 and the second pattern (q1-->g-->q4)
//  flow code 269  ( you can ask HepMC::Flow to choose
//  a unique code for you using Flow::set_unique_icode() ).
//  The first two code indices are reserved for color codes, so we store 
//  these codes with the particles as follows:
//    q2->flow().set_icode(1,678);
//    g->flow().set_icode(1,678);
//    q3->flow().set_icode(1,678);
//    q1->flow().set_icode(1,269);
//    g->flow().set_icode(2,269);
//    q4->flow().set_icode(1,269);
//  later on if we wish to know the color partner of q1 we can ask for a list
//  of all particles connected via this code to q1 which do have less than 
//  2 color partners using:
//    vector<GenParticle*> result=q1->dangling_connected_partners(q1->icode(1),1,2);
//  this will return a list containing q1 and q4.
//    vector<GenParticle*> result=q1->connected_partners(q1->icode(1),1,2);
//  would return a list containing q1, g, and q4.
//

#include <iostream>
#include <map>
#include <vector>

namespace HepMC {

    class GenParticle;

    //! The flow object

    ///
    /// \class  Flow
    /// The particle's flow object
    /// keeps track of an arbitrary number of flow patterns within a graph 
    /// (i.e. color flow, charge flow, lepton number flow, ...) 
    /// Flow patterns are coded with an integer, in the same manner as in Herwig.
    class Flow {

        /// for printing
	friend std::ostream& operator<<( std::ostream& ostr, const Flow& f );
	
    public:
        /// default constructor
	Flow( GenParticle* particle_owner = 0 );
	/// copy
	Flow( const Flow& );
	virtual         ~Flow();
        /// swap
        void swap( Flow & other);
	/// make a copy
	Flow&           operator=( const Flow& );
	/// equality
	bool            operator==( const Flow& a ) const; //compares only flow
	/// inequality
	bool            operator!=( const Flow& a ) const; //patterns not owner

        /// print Flow information to ostr
	void            print( std::ostream& ostr = std::cout ) const;

	/// returns all connected particles which have "code" in any  of the 
	///  num_indices beginning with index code_index.
	std::vector<HepMC::GenParticle*> connected_partners( int code, int code_index =1,
						   int num_indices = 2 ) const;
	/// same as connected_partners, but returns only those particles which
	///  are connected to <=1 other particles (i.e. the flow line "dangles"
	///  at these particles)
	std::vector<HepMC::GenParticle*> dangling_connected_partners( int code, 
			       int code_index = 1, int num_indices = 2 ) const;

	////////////////////
	// access methods //
	////////////////////

	/// find particle owning this Flow
	const GenParticle* particle_owner() const;
	/// flow code
	int             icode( int code_index = 1 ) const;
	/// set flow code
	Flow            set_icode( int code_index, int code );
	/// set unique flow code
	Flow            set_unique_icode( int code_index = 1 );

	//////////////////////
	// container access //
	//////////////////////

        /// return true if there is no flow container
	bool            empty() const;
	/// size of flow pattern container
	int             size() const;
	/// clear flow patterns
        void            clear();
	/// empty flow pattern container
	bool            erase( int code_index );

        /// iterator for flow pattern container
        typedef std::map<int,int>::iterator       iterator;
        /// const iterator for flow pattern container
        typedef std::map<int,int>::const_iterator const_iterator;
	/// beginning of flow pattern container
        iterator            begin();
	/// end of flow pattern container
        iterator            end();
	/// beginning of flow pattern container
        const_iterator      begin() const;
	/// end of flow pattern container
        const_iterator      end() const;

    protected: // intended for internal use only
        /// for internal use only
	void            connected_partners( std::vector<HepMC::GenParticle*>* output, 
					    int code,
					    int code_index,
					    int num_indices ) const;
        /// for internal use only
	void            dangling_connected_partners( std::vector<HepMC::GenParticle*>* 
						     output, 
						     std::vector<HepMC::GenParticle*>*
						     visited_particles, 
						     int code, int code_index, 
						     int num_indices ) const; 
    private:
	GenParticle*         m_particle_owner;
	std::map<int,int> m_icode; // stores flow patterns as(code_index,icode)
    };  

    ///////////////////////////
    // INLINE Access Methods //
    ///////////////////////////

    inline const GenParticle* Flow::particle_owner() const {
	return m_particle_owner;
    }
    inline int Flow::icode( int code_index ) const {
	std::map<int,int>::const_iterator a = m_icode.find(code_index);
	return a==m_icode.end() ? 0 : (*a).second;
    }
    inline Flow Flow::set_icode( int code_index, int code ) {
	m_icode[code_index] = code;
	return *this;
    }
    inline Flow Flow::set_unique_icode( int flow_num ) {
	/// use this method if you want to assign a unique flow code, but
	/// do not want the burden of choosing it yourself
	m_icode[flow_num] = size_t(this);
	return *this;
    }
    inline bool Flow::empty() const { return (bool)m_icode.empty(); }
    inline int Flow::size() const { return (int)m_icode.size(); }
    inline void Flow::clear() { m_icode.clear(); }
    inline bool Flow::erase( int code_index ) {
	// this will return true if the number of elements removed is nonzero
	return m_icode.erase( code_index )==0 ? false : true ;
    }
    inline Flow::iterator Flow::begin() { return m_icode.begin(); }
    inline Flow::iterator Flow::end() { return m_icode.end(); }
    inline Flow::const_iterator Flow::begin() const { return m_icode.begin(); }
    inline Flow::const_iterator Flow::end() const { return m_icode.end(); }

    ///////////////////////////
    // INLINE Operators      //
    ///////////////////////////

    inline bool Flow::operator==( const Flow& a ) const {
	/// equivalent flows have the same flow codes for all flow_numbers 
	/// (i.e. their m_icode maps are identical), but they need not have the
	/// same m_particle owner
	return (m_icode == a.m_icode);
    }
    inline bool Flow::operator!=( const Flow& a ) const {
	return !( *this == a );
    }
    inline Flow& Flow::operator=( const Flow& inflow ) {
	/// copies only the m_icode ... not the particle_owner
	/// this is intuitive behaviour so you can do
	/// oneparticle->flow() = otherparticle->flow()
	//
	m_icode = inflow.m_icode;
	return *this;
    }

} // HepMC

#endif  // HEPMC_FLOW_H
//--------------------------------------------------------------------------

HepMC/HeavyIon.h000644 000765 000024 00000026570 13310701701 013720 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HEPMC_HEAVY_ION_H
#define HEPMC_HEAVY_ION_H

//////////////////////////////////////////////////////////////////////////
// garren@fnal.gov, February 2006
//
// Additional information storage for Heavy Ion generators
//////////////////////////////////////////////////////////////////////////
//
//      int   Ncoll_hard          // Number of hard scatterings
//      int   Npart_proj          // Number of projectile participants
//      int   Npart_targ          // Number of target participants
//      int   Ncoll               // Number of NN (nucleon-nucleon) collisions
//      int   N_Nwounded_collisions        // Number of N-Nwounded collisions
//      int   Nwounded_N_collisions        // Number of Nwounded-N collisons
//      int   Nwounded_Nwounded_collisions // Number of Nwounded-Nwounded collisions
//      int   spectator_neutrons           // Number of spectator neutrons
//      int   spectator_protons            // Number of spectator protons
//      float impact_parameter        // Impact Parameter(fm) of collision
//      float event_plane_angle       // Azimuthal angle of event plane
//      float eccentricity            // eccentricity of participating nucleons
//                                        in the transverse plane 
//                                        (as in phobos nucl-ex/0510031) 
//      float sigma_inel_NN           // nucleon-nucleon inelastic 
//                                        (including diffractive) cross-section
//
//////////////////////////////////////////////////////////////////////////
// Feb. 17, 2006: adjust names according to suggestions from Heavy Ion users
// Feb.  7, 2006: first pass at making method names consistent with existing
//                HepMC code
//////////////////////////////////////////////////////////////////////////

namespace HepMC {


//! The HeavyIon class stores information about heavy ions

///
/// \class  HeavyIon
/// HepMC::HeavyIon provides
/// additional information storage for Heavy Ion generators in GenEvent.
/// Creation and use of this information is optional. 
///
class HeavyIon {

public:
  // ---  birth/death:
  //
  /// default constructor
  HeavyIon()
    : m_Ncoll_hard(0), 
      m_Npart_proj(0),
      m_Npart_targ(0),
      m_Ncoll(0),
      m_spectator_neutrons(0),
      m_spectator_protons(0),
      m_N_Nwounded_collisions(0),
      m_Nwounded_N_collisions(0),
      m_Nwounded_Nwounded_collisions(0),
      m_impact_parameter(0),
      m_event_plane_angle(0),
      m_eccentricity(0),
      m_sigma_inel_NN(0)
    {}

  /// The first 6 values must be provided.
  HeavyIon( int nh, int np, int nt, int nc, int ns, int nsp,
                   int nnw=0, int nwn=0, int nwnw=0, 
		   float im=0., float pl=0., float ec=0., float s=0. );

  ~HeavyIon() {}
    
  // ---  copying:
  //
  HeavyIon( HeavyIon const & orig );	//!< copy constructor
  HeavyIon &  operator = ( HeavyIon const & rhs ); //!< make a copy
  void swap( HeavyIon & other );	//!< swap two HeavyIon objects

  // ---  equivalence:
  //
  bool    operator==( const HeavyIon& ) const; //!< check for equality
  bool    operator!=( const HeavyIon& ) const; //!< check for inequality

  // ---  accessors:
    /// Number of hard scatterings
    int   Ncoll_hard()                   const { return m_Ncoll_hard; }
    /// Number of projectile participants
    int   Npart_proj()                   const { return m_Npart_proj; }
    /// Number of target participants
    int   Npart_targ()                   const { return m_Npart_targ; }
    /// Number of NN (nucleon-nucleon) collisions
    int   Ncoll()                        const { return m_Ncoll; }
    /// Number of spectator neutrons
    int   spectator_neutrons()           const { return m_spectator_neutrons; }
    /// Number of spectator protons
    int   spectator_protons()            const { return m_spectator_protons; }
    /// Number of N-Nwounded collisions
    int   N_Nwounded_collisions()        const { return m_N_Nwounded_collisions; }
    /// Number of Nwounded-N collisons
    int   Nwounded_N_collisions()        const { return m_Nwounded_N_collisions; }
    /// Number of Nwounded-Nwounded collisions
    int   Nwounded_Nwounded_collisions() const { return m_Nwounded_Nwounded_collisions; }
    /// Impact Parameter(in fm) of collision
    float impact_parameter()             const { return m_impact_parameter; }
    /// Azimuthal angle of event plane
    float event_plane_angle()            const { return m_event_plane_angle; }
    /// eccentricity of participating nucleons in the transverse plane 
    /// (as in phobos nucl-ex/0510031)
    float eccentricity()                 const { return m_eccentricity;  }
    /// nucleon-nucleon inelastic (including diffractive) cross-section
    float sigma_inel_NN()                const { return m_sigma_inel_NN; }

    /// verify that the instance contains non-zero information
    bool  is_valid()                     const;

  // ---  mutators:
    /// set number of hard scatterings
    void   set_Ncoll_hard(const int &i)              { m_Ncoll_hard=i; }
    /// set number of projectile participants
    void   set_Npart_proj(const int &i)              { m_Npart_proj=i; }
    /// set number of target participants
    void   set_Npart_targ(const int &i)              { m_Npart_targ=i; }
    /// set number of NN (nucleon-nucleon) collisions
    void   set_Ncoll(const int &i)                   { m_Ncoll=i; }
    /// set number of spectator neutrons
    void   set_spectator_neutrons(const int &i)      { m_spectator_neutrons=i; }
    /// set number of spectator protons
    void   set_spectator_protons(const int &i)       { m_spectator_protons=i; }
    /// set number of N-Nwounded collisions
    void   set_N_Nwounded_collisions(const int &i)   { m_N_Nwounded_collisions=i; }
    /// set number of Nwounded-N collisons
    void   set_Nwounded_N_collisions(const int &i)   { m_Nwounded_N_collisions=i; }
    /// set number of Nwounded-Nwounded collisions
    void   set_Nwounded_Nwounded_collisions(const int &i) 
                                          { m_Nwounded_Nwounded_collisions=i; }
    /// set Impact Parameter in fm
    void   set_impact_parameter(const float &f)      { m_impact_parameter=f; }
    /// set azimuthal angle of event plane
    void   set_event_plane_angle(const float &f)     { m_event_plane_angle=f; }
    /// set eccentricity of participating nucleons in the transverse plane 
    void   set_eccentricity(const float &f)          { m_eccentricity=f;  }
    /// set nucleon-nucleon inelastic cross-section
    void   set_sigma_inel_NN(const float &f)         { m_sigma_inel_NN=f; }

private: // data members
    int   m_Ncoll_hard; 
    int   m_Npart_proj; 
    int   m_Npart_targ;
    int   m_Ncoll;
    int   m_spectator_neutrons;
    int   m_spectator_protons;
    int   m_N_Nwounded_collisions;
    int   m_Nwounded_N_collisions;
    int   m_Nwounded_Nwounded_collisions;
    float m_impact_parameter;
    float m_event_plane_angle;
    float m_eccentricity; 
    float m_sigma_inel_NN;

};

// Free Functions

/// Write the contents of HeavyIon to an output stream.
std::ostream & operator << (std::ostream &, HeavyIon const *);
/// Read the contents of HeavyIon from an input stream.
std::istream & operator >> (std::istream &, HeavyIon *);

// inline operators
  /// Required members are
  /// the number of hard scatterings,
  /// the number of projectile participants.
  /// the number of target participants.
  /// the number of nucleon-nucleon collisions,
  /// the number of spectator neutrons, and
  /// the number of spectator protons.
inline HeavyIon::HeavyIon( int nh, int np, int nt, int nc, int ns, int nsp,
                   int nnw, int nwn, int nwnw, 
		   float im, float pl, float ec, float s )
    : m_Ncoll_hard(nh), 
      m_Npart_proj(np),
      m_Npart_targ(nt),
      m_Ncoll(nc),
      m_spectator_neutrons(ns),
      m_spectator_protons(nsp),
      m_N_Nwounded_collisions(nnw),
      m_Nwounded_N_collisions(nwn),
      m_Nwounded_Nwounded_collisions(nwnw),
      m_impact_parameter(im),
      m_event_plane_angle(pl),
      m_eccentricity(ec),
      m_sigma_inel_NN(s)
   {}

inline HeavyIon::HeavyIon( HeavyIon const & orig )
    : m_Ncoll_hard(orig.m_Ncoll_hard), 
      m_Npart_proj(orig.m_Npart_proj),
      m_Npart_targ(orig.m_Npart_targ),
      m_Ncoll(orig.m_Ncoll),
      m_spectator_neutrons(orig.m_spectator_neutrons),
      m_spectator_protons(orig.m_spectator_protons),
      m_N_Nwounded_collisions(orig.m_N_Nwounded_collisions),
      m_Nwounded_N_collisions(orig.m_Nwounded_N_collisions),
      m_Nwounded_Nwounded_collisions(orig.m_Nwounded_Nwounded_collisions),
      m_impact_parameter(orig.m_impact_parameter),
      m_event_plane_angle(orig.m_event_plane_angle),
      m_eccentricity(orig.m_eccentricity),
      m_sigma_inel_NN(orig.m_sigma_inel_NN)
   {}

inline HeavyIon &  HeavyIon::operator = ( HeavyIon const & rhs ) 
{
  HeavyIon temp( rhs );
  swap( temp );
  return *this;
}

inline void HeavyIon::swap( HeavyIon & other ) 
{
  std::swap(m_Ncoll_hard, other.m_Ncoll_hard); 
  std::swap(m_Npart_proj, other.m_Npart_proj);
  std::swap(m_Npart_targ, other.m_Npart_targ);
  std::swap(m_Ncoll, other.m_Ncoll);
  std::swap(m_N_Nwounded_collisions, other.m_N_Nwounded_collisions);
  std::swap(m_Nwounded_N_collisions, other.m_Nwounded_N_collisions);
  std::swap(m_Nwounded_Nwounded_collisions, other.m_Nwounded_Nwounded_collisions);
  std::swap(m_spectator_neutrons, other.m_spectator_neutrons);
  std::swap(m_spectator_protons, other.m_spectator_protons);
  std::swap(m_impact_parameter, other.m_impact_parameter);
  std::swap(m_event_plane_angle, other.m_event_plane_angle);
  std::swap(m_eccentricity, other.m_eccentricity);
  std::swap(m_sigma_inel_NN, other.m_sigma_inel_NN);
}

inline bool    HeavyIon::operator==( const HeavyIon& a ) const
{
    /// equality requires that each member match
    return ( a.Ncoll_hard() == this->Ncoll_hard() 
             && a.Npart_proj() == this->Npart_proj() 
             && a.Npart_targ() == this->Npart_targ() 
	     && a.Ncoll() == this->Ncoll() 
             && a.N_Nwounded_collisions() == this->N_Nwounded_collisions() 
	     && a.Nwounded_N_collisions() == this->Nwounded_N_collisions() 
             && a.Nwounded_Nwounded_collisions() == this->Nwounded_Nwounded_collisions() 
	     && a.spectator_neutrons() == this->spectator_neutrons() 
             && a.spectator_protons() == this->spectator_protons() 
	     && a.impact_parameter() == this->impact_parameter() 
	     && a.event_plane_angle() == this->event_plane_angle()
	     && a.eccentricity() == this->eccentricity() 
	     && a.sigma_inel_NN() == this->sigma_inel_NN() );
}

inline bool    HeavyIon::operator!=( const HeavyIon& a ) const
{
    /// any nonmatching member generates inequality
    return !( a == *this );
}

inline bool  HeavyIon::is_valid() const
{
    if( m_Ncoll_hard != 0 ) return true;
    if( m_Npart_proj != 0 ) return true;
    if( m_Npart_targ != 0 ) return true;
    if( m_Ncoll != 0 ) return true;
    if( m_spectator_neutrons != 0 ) return true;
    if( m_spectator_protons != 0 ) return true;
    if( m_N_Nwounded_collisions != 0 ) return true;
    if( m_Nwounded_N_collisions != 0 ) return true;
    if( m_Nwounded_Nwounded_collisions != 0 ) return true;
    if( m_impact_parameter != 0 ) return true;
    if( m_event_plane_angle != 0 ) return true;
    if( m_eccentricity != 0 ) return true;
    if( m_sigma_inel_NN != 0 ) return true;
    return false;
}

} // HepMC

#endif  // HEPMC_HEAVY_ION_H
HepMC/Polarization.h000644 000765 000024 00000007762 13310701701 014653 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HEPMC_POLARIZATION_H
#define HEPMC_POLARIZATION_H

//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, September 1999, refer to:
// M. Dobbs and J.B. Hansen, "The HepMC C++ Monte Carlo Event Record for
// High Energy Physics", Computer Physics Communications (to be published).
//
// Polarization object for a particle. All angles are in radians.
//////////////////////////////////////////////////////////////////////////

#include "HepMC/SimpleVector.h"
#include <iostream>
#include <cmath>

namespace HepMC {

    static const double HepMC_pi = 3.14159265358979323846;  // copy of pi from CLHEP
    
    //! The Polarization class stores theta and phi for a GenParticle

    ///
    /// \class  Polarization
    /// HepMC::Polarization stores a particle's theta and phi in radians.
    /// Use of this information is optional. 
    /// By default, the polarization is set to zero.
    ///
    class Polarization {

        /// print polarization information
	friend std::ostream& operator<<( std::ostream&, const Polarization& );

    public:
        /// default constructor
	Polarization( );
	/// constructor requiring at least one value
	Polarization( double theta, double phi = 0 );
	/// construct from another polarization object
	Polarization( const Polarization& inpolar );
	/// construct using the polar and azimuthal angles from a ThreeVector
	Polarization( const ThreeVector& vec3in );
	virtual       ~Polarization() {}

        /// swap
        void swap( Polarization & other);
	/// make a copy
	Polarization& operator=( const Polarization& inpolar );
	/// equality requires that theta and phi are equal
	bool          operator==( const Polarization& ) const;
	/// inequality results if either theta or phi differ
	bool          operator!=( const Polarization& ) const;

        /// print theta and phi
	void          print( std::ostream& ostr = std::cout ) const;
    
	////////////////////
	// access methods //
	////////////////////
	double        theta() const;    //!< returns polar angle in radians
	double        phi() const;      //!< returns azimuthal angle in radians
	ThreeVector   normal3d() const; //!< unit 3 vector for easy manipulation
	bool          is_defined() const;   //!< returns true if the Polarization has been defined

        /// set polar angle in radians 
	double        set_theta( double theta );
        /// set azimuthal angle in radians 
	double        set_phi( double phi );
	/// set both polar and azimuthal angles in radians
	void          set_theta_phi( double theta, double phi );
	/// sets polarization according to direction of 3 vec
	ThreeVector   set_normal3d( const ThreeVector& vec3in ); 
	/// declares the Polarization as undefined and zeros the values
	void          set_undefined();

    private:
    	/// private method to return a polar angle in the correct range
	double valid_theta( double theta );
    	/// private method to return an azimuthal angle in the correct range
	double valid_phi( double phi );

    private:
	double m_theta; //polar angle of polarization in radians 0< theta <pi
	double m_phi;   //azimuthal angle of polarization in rad. 0< phi <2pi
	bool   m_defined; //used to flag if the Polarization has been defined
    };

    ///////////////////////////
    // INLINE Access Methods //
    ///////////////////////////

    inline double Polarization::theta() const { return m_theta; }
    inline double Polarization::phi() const { return m_phi; }

    ///////////////////////////
    // INLINE Operators      //
    ///////////////////////////

    inline bool Polarization::operator==( const Polarization& a ) const 
    {
	return ( a.theta() == this->theta() && a.phi() == this->phi() && a.is_defined() == this->is_defined() );
    }

    inline bool Polarization::operator!=(const Polarization& a ) const 
    {
	return !( a == *this );
    }

} // HepMC

#endif  // HEPMC_POLARIZATION_H
//--------------------------------------------------------------------------
HepMC/IteratorRange.h000644 000765 000024 00000001430 13310701701 014730 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HEPMC_ITERATOR_RANGE_H
#define HEPMC_ITERATOR_RANGE_H

//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, September 1999, refer to:
// M. Dobbs and J.B. Hansen, "The HepMC C++ Monte Carlo Event Record for
// High Energy Physics", Computer Physics Communications (to be published).
//
// Iterator ranges used by GenVertex
//////////////////////////////////////////////////////////////////////////


namespace HepMC {

    /// type of iteration
    enum IteratorRange { parents, children, family, 
			 ancestors, descendants, relatives };
} // HepMC

#endif  // HEPMC_ITERATOR_RANGE_H
//--------------------------------------------------------------------------
HepMC/is_arithmetic.h000644 000765 000024 00000003702 13310701701 015012 0ustar00sa639staff000000 000000 #ifndef IS_ARITHMETIC
#define IS_ARITHMETIC
// author: Walter Brown

// ----------------------------------------------------------------------
// prolog

namespace HepMC  {

///
/// \namespace detail
/// internal namespace
///
namespace detail  {


// ----------------------------------------------------------------------
// is_arithmetic<>

/// undefined and therefore non-arithmetic
template< class T >
  struct is_arithmetic
{
  static  bool const  value = false;
};

/// character is arithmetic
template<>
  struct is_arithmetic<char>
{ static  bool const  value = true; };

/// unsigned character is arithmetic
template<>
  struct is_arithmetic<unsigned char>
{ static  bool const  value = true; };

/// signed character is arithmetic
template<>
  struct is_arithmetic<signed char>
{ static  bool const  value = true; };

/// short is arithmetic
template<>
  struct is_arithmetic<short>
{ static  bool const  value = true; };

/// unsigned short is arithmetic
template<>
  struct is_arithmetic<unsigned short>
{ static  bool const  value = true; };

/// int is arithmetic
template<>
  struct is_arithmetic<int>
{ static  bool const  value = true; };

/// unsigned int is arithmetic
template<>
  struct is_arithmetic<unsigned int>
{ static  bool const  value = true; };

/// long is arithmetic
template<>
  struct is_arithmetic<long>
{ static  bool const  value = true; };

/// unsigned long is arithmetic
template<>
  struct is_arithmetic<unsigned long>
{ static  bool const  value = true; };

/// float is arithmetic
template<>
  struct is_arithmetic<float>
{ static  bool const  value = true; };

/// double is arithmetic
template<>
  struct is_arithmetic<double>
{ static  bool const  value = true; };

/// long double is arithmetic
template<>
  struct is_arithmetic<long double>
{ static  bool const  value = true; };


// ----------------------------------------------------------------------
// epilog

}  // namespace detail
}  // namespace HepMC

#endif  // IS_ARITHMETIC
HepMC/GenRanges.h000644 000765 000024 00000025303 13310701701 014040 0ustar00sa639staff000000 000000 #ifndef HEPMC_GEN_EVENT_ITERATORS_H
#define HEPMC_GEN_EVENT_ITERATORS_H

//--------------------------------------------------------------------------
//////////////////////////////////////////////////////////////////////////
// garren@fnal.gov, May 2009
// 
//////////////////////////////////////////////////////////////////////////
//--------------------------------------------------------------------------

#include <stdexcept>

#include "HepMC/GenEvent.h"
#include "HepMC/GenVertex.h"

namespace HepMC {

//! GenEventVertexRange acts like a collection of vertices

///
/// \class  GenEventVertexRange
/// HepMC::GenEventVertexRange is used to mimic a collection of
/// vertices for ease of use - especially with utilities such as 
/// the Boost foreach funtion
///
class GenEventVertexRange {

public:

  /// the constructor requires a GenEvent
  GenEventVertexRange( GenEvent & e ) : m_event(e) {}
  /// 
  GenEvent::vertex_iterator begin() { return m_event.vertices_begin(); }
  GenEvent::vertex_iterator end()   { return m_event.vertices_end(); }

private:
  /// Because the class contains a reference, assignments are not allowed.
  /// However, we need the copy constructor for GenEvent::vertex_range().
  GenEventVertexRange& operator=( GenEventVertexRange & );

private:
  GenEvent & m_event;

};

//! ConstGenEventVertexRange acts like a collection of vertices

///
/// \class  ConstGenEventVertexRange
/// HepMC::ConstGenEventVertexRange is used to mimic a collection of
/// vertices for ease of use - especially with utilities such as 
/// the Boost foreach funtion
/// This is the const partner of GenEventVertexRange
///
class ConstGenEventVertexRange {

public:

  /// the constructor requires a const GenEvent
  ConstGenEventVertexRange( GenEvent const & e ) : m_event(e) {}
  /// 
  GenEvent::vertex_const_iterator begin() const { return m_event.vertices_begin(); }
  GenEvent::vertex_const_iterator end()   const { return m_event.vertices_end(); }

private:
  /// Because the class contains a reference, assignments are not allowed.
  /// However, we need the copy constructor for GenEvent::vertex_range().
  ConstGenEventVertexRange& operator=( ConstGenEventVertexRange & );

private:
  GenEvent const & m_event;

};

//! GenEventParticleRange acts like a collection of particles

///
/// \class  GenEventParticleRange
/// HepMC::GenEventParticleRange is used to mimic a collection of
/// particles for ease of use - especially with utilities such as 
/// the Boost foreach funtion
///
class GenEventParticleRange {

public:

  /// the constructor requires a GenEvent
  GenEventParticleRange( GenEvent & e ) : m_event(e) {}
  /// 
  GenEvent::particle_iterator begin() { return m_event.particles_begin(); }
  GenEvent::particle_iterator end()   { return m_event.particles_end(); }

private:
  /// Because the class contains a reference, assignments are not allowed.
  /// However, we need the copy constructor for GenEvent::particle_range().
  GenEventParticleRange& operator=( GenEventParticleRange & );

private:
  GenEvent & m_event;

};

//! ConstGenEventParticleRange acts like a collection of particles

///
/// \class  ConstGenEventParticleRange
/// HepMC::ConstGenEventParticleRange is used to mimic a collection of
/// particles for ease of use - especially with utilities such as 
/// the Boost foreach funtion
/// This is the const partner of GenEventParticleRange
///
class ConstGenEventParticleRange {

public:

  /// the constructor requires a const GenEvent
  ConstGenEventParticleRange( GenEvent const & e ) : m_event(e) {}
  /// 
  GenEvent::particle_const_iterator begin() const { return m_event.particles_begin(); }
  GenEvent::particle_const_iterator end()   const { return m_event.particles_end(); }

private:
  /// Because the class contains a reference, assignments are not allowed.
  /// However, we need the copy constructor for GenEvent::particle_range().
  ConstGenEventParticleRange& operator=( ConstGenEventParticleRange & );

private:
  GenEvent const & m_event;

};

//! GenVertexParticleRange acts like a collection of particles

///
/// \class  GenVertexParticleRange
/// HepMC::GenVertexParticleRange is used to mimic a collection of
/// particles for ease of use - especially with utilities such as 
/// the Boost foreach funtion
///
class GenVertexParticleRange {

public:

  /// the constructor requires a GenVertex
  GenVertexParticleRange( GenVertex & v, IteratorRange range = relatives ) 
  : m_vertex(v),m_range(range) {}
  /// 
  GenVertex::particle_iterator begin() { return m_vertex.particles_begin(m_range); }
  GenVertex::particle_iterator end()   { return m_vertex.particles_end(m_range); }

private:
  /// Because the class contains a reference, assignments are not allowed.
  /// However, we need the copy constructor for GenVertex::particles().
  GenVertexParticleRange& operator=( GenVertexParticleRange & );

private:
  GenVertex     & m_vertex;
  IteratorRange   m_range;

};

//! GenParticleProductionRange acts like a collection of particles

///
/// \class  GenParticleProductionRange
/// HepMC::GenParticleProductionRange is used to mimic a collection of
/// particles associated with the particle's production vertex for ease of use 
/// Utilities such as the Boost foreach funtion will want to use this class.
///
class GenParticleProductionRange {

public:

  /// the constructor requires a GenParticle
  GenParticleProductionRange( GenParticle const & p, IteratorRange range = relatives ) 
  : m_particle(p),m_range(range) {}
  /// begin iterator throws an error if the particle production_vertex is undefined
  GenVertex::particle_iterator begin();
  /// end iterator throws an error if the particle production_vertex is undefined
  GenVertex::particle_iterator end(); 

private:
  /// Because the class contains a reference, assignments are not allowed.
  /// However, we need the copy constructor for GenVertex::particles_in().
  GenParticleProductionRange& operator=( GenParticleProductionRange & );

private:
  GenParticle const & m_particle;
  IteratorRange       m_range;

};

class ConstGenParticleProductionRange {

public:

  /// the constructor requires a GenParticle
  ConstGenParticleProductionRange( GenParticle const & p, IteratorRange range = relatives ) 
  : m_particle(p),m_range(range) {}
  /// begin iterator throws an error if the particle production_vertex is undefined
  GenVertex::particle_iterator begin();
  /// end iterator throws an error if the particle production_vertex is undefined
  GenVertex::particle_iterator end(); 

private:
  /// Because the class contains a reference, assignments are not allowed.
  /// However, we need the copy constructor for GenVertex::particles_in().
  ConstGenParticleProductionRange& operator=( ConstGenParticleProductionRange & );

private:
  GenParticle const & m_particle;
  IteratorRange       m_range;

};

//! GenParticleEndRange acts like a collection of particles

///
/// \class  GenParticleEndRange
/// HepMC::GenParticleEndRange is used to mimic a collection of
/// particles associated with the particle's end vertex for ease of use 
/// Utilities such as the Boost foreach funtion will want to use this class.
///
class GenParticleEndRange {

public:

  /// the constructor requires a GenParticle
  GenParticleEndRange( GenParticle const & p, IteratorRange range = relatives ) 
  : m_particle(p),m_range(range) {}
  /// begin iterator throws an error if the particle end_vertex is undefined
  GenVertex::particle_iterator begin();
  /// end iterator throws an error if the particle end_vertex is undefined
  GenVertex::particle_iterator end(); 

private:
  /// Because the class contains a reference, assignments are not allowed.
  /// However, we need the copy constructor for GenVertex::particles_out().
  GenParticleEndRange& operator=( GenParticleEndRange & );

private:
  GenParticle const & m_particle;
  IteratorRange       m_range;

};

class ConstGenParticleEndRange {

public:

  /// the constructor requires a GenParticle
  ConstGenParticleEndRange( GenParticle const & p, IteratorRange range = relatives ) 
  : m_particle(p),m_range(range) {}
  /// begin iterator throws an error if the particle end_vertex is undefined
  GenVertex::particle_iterator begin();
  /// end iterator throws an error if the particle end_vertex is undefined
  GenVertex::particle_iterator end(); 

private:
  /// Because the class contains a reference, assignments are not allowed.
  /// However, we need the copy constructor for GenVertex::particles_out().
  ConstGenParticleEndRange& operator=( ConstGenParticleEndRange & );

private:
  GenParticle const & m_particle;
  IteratorRange       m_range;

};


inline GenVertex::particle_iterator GenParticleProductionRange::begin() 
{ 
    if ( ! m_particle.production_vertex() ) 
        throw(std::range_error("GenParticleProductionRange: GenParticle has no production_vertex"));
    return m_particle.production_vertex()->particles_begin(m_range); 
}

inline GenVertex::particle_iterator GenParticleProductionRange::end()
{ 
    if ( ! m_particle.production_vertex() ) 
        throw(std::range_error("GenParticleProductionRange: GenParticle has no production_vertex"));
    return m_particle.production_vertex()->particles_end(m_range); 
}


inline GenVertex::particle_iterator ConstGenParticleProductionRange::begin() 
{ 
    if ( ! m_particle.production_vertex() ) 
        throw(std::range_error("ConstGenParticleProductionRange: GenParticle has no production_vertex"));
    return m_particle.production_vertex()->particles_begin(m_range); 
}

inline GenVertex::particle_iterator ConstGenParticleProductionRange::end()
{ 
    if ( ! m_particle.production_vertex() ) 
        throw(std::range_error("ConstGenParticleProductionRange: GenParticle has no production_vertex"));
    return m_particle.production_vertex()->particles_end(m_range); 
}

inline GenVertex::particle_iterator GenParticleEndRange::begin() 
{ 
    if ( ! m_particle.end_vertex() ) 
        throw(std::range_error("GenParticleEndRange: GenParticle has no end_vertex"));
    return m_particle.end_vertex()->particles_begin(m_range); 
}
inline GenVertex::particle_iterator GenParticleEndRange::end()
{ 
    if ( ! m_particle.end_vertex() ) 
        throw(std::range_error("GenParticleEndRange: GenParticle has no end_vertex"));
    return m_particle.end_vertex()->particles_end(m_range); 
}

inline GenVertex::particle_iterator ConstGenParticleEndRange::begin() 
{ 
    if ( ! m_particle.end_vertex() ) 
        throw(std::range_error("ConstGenParticleEndRange: GenParticle has no end_vertex"));
    return m_particle.end_vertex()->particles_begin(m_range); 
}
inline GenVertex::particle_iterator ConstGenParticleEndRange::end()
{ 
    if ( ! m_particle.end_vertex() ) 
        throw(std::range_error("ConstGenParticleEndRange: GenParticle has no end_vertex"));
    return m_particle.end_vertex()->particles_end(m_range); 
}

} // HepMC

#endif  // HEPMC_GEN_EVENT_ITERATORS_H
HepMC/PythiaWrapper.h000644 000765 000024 00000003005 13310701701 014761 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef PYTHIA_WRAPPER_H
#define PYTHIA_WRAPPER_H

//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, November 2000
// Just a link to whichever pythia version is current.
//////////////////////////////////////////////////////////////////////////

// This pre-compiler directive is included (2002-01-16) to allow compatibility
// with MS Visual C++, which interfaces to fortran in a different manner.
// For it to work you need to define the _WIN32 variable when compiling.
#ifdef _WIN32 // Platform: Windows MS Visual C++
#include "HepMC/PythiaWrapper6_4_WIN32.h"

#else // Generic version, tested on Linux ecgs/gcc
#include "HepMC/PythiaWrapper6_4.h"

#endif // Platform

#include <cmath>

#include "HepMC/GenCrossSection.h"

namespace HepMC {

/// calculate the Pythia cross section and statistical error
inline GenCrossSection getPythiaCrossSection() {

  GenCrossSection xsec;
  // xsec(0,2) contains the sum of differential cross sections in mb
  // ngen(0,2) contains the combined number of generated events
  // convert to pb (HepMC convention)
  double xsecval = pyint5.xsec[2][0] * 1.0e9;
  // statistical error
  double xsecerr = xsecval / std::sqrt( (double)pyint5.ngen[2][0] );
  // set and return cross section information
  xsec.set_cross_section(xsecval, xsecerr);
  return xsec;
}



} // HepMC

#endif  // PYTHIA_WRAPPER_H
//--------------------------------------------------------------------------
HepMC/HerwigWrapper.h000644 000765 000024 00000013433 13310701701 014756 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HERWIG_WRAPPER_H
#define HERWIG_WRAPPER_H

//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, November 200July 2002
// Wrapper for FORTRAN version of Herwig
// Wrapper written for Herwig 6.4
//////////////////////////////////////////////////////////////////////////

#include <ctype.h>

// This pre-compiler directive is included (2002-01-16) to allow compatibility
// with MS Visual C++, which interfaces to fortran in a different manner.
// For it to work you need to define the _WIN32 variable when compiling.
#ifdef _WIN32 // Platform: Windows MS Visual C++

//  Sorry, there is NO version currently available for Vusual C++.

#else // Generic version, tested on Linux ecgs/gcc

//--------------------------------------------------------------------------
// HERWIG Common Block Declarations

//        COMMON/HWPROC/EBEAM1,EBEAM2,PBEAM1,PBEAM2,IPROC,MAXEV
extern "C" {
    extern struct {
	double EBEAM1,EBEAM2,PBEAM1,PBEAM2;
        int IPROC,MAXEV;
    } hwproc_;
}
#define hwproc hwproc_

//      COMMON/HWBEAM/IPART1,IPART2
extern "C" {
    extern struct {
        int IPART1, IPART2;
    } hwbeam_;
}
#define hwbeam hwbeam_

//        CHARACTER*8 PART1,PART2
//        COMMON/HWBMCH/PART1,PART2
extern "C" {
    extern struct {
	char PART1[8],PART2[8];
    } hwbmch_;
}
#define hwbmch hwbmch_

//  COMMON/HWEVNT/AVWGT,EVWGT,GAMWT,TLOUT,WBIGST,WGTMAX,WGTSUM,WSQSUM,
//       & IDHW(NMXHEP),IERROR,ISTAT,LWEVT,MAXER,MAXPR,NOWGT,NRN(2),NUMER,
//       & NUMERU,NWGTS,GENSOF
const int herwig_hepevt_size = 4000;
extern "C" {
    extern struct hwgev {
	double AVWGT,EVWGT,GAMWT,TLOUT,WBIGST,WGTMAX,WGTSUM,WSQSUM;
	int IDHW[herwig_hepevt_size],IERROR,ISTAT,LWEVT,MAXER,MAXPR;
	int NOWGT,NRN[2],NUMER,NUMERU,NWGTS;
	int GENSOF; //Beware! in F77 this is logical
    } hwevnt_;
}
#define hwevnt hwevnt_

//  C Basic parameters (and quantities derived from them)
//        COMMON/HWPRAM/AFCH(16,2),ALPHEM,B1LIM,BETAF,BTCLM,CAFAC,CFFAC,
//       & CLMAX,CLPOW,CLSMR(2),CSPEED,ENSOF,ETAMIX,F0MIX,F1MIX,F2MIX,GAMH,
//       & GAMW,GAMZ,GAMZP,GEV2NB,H1MIX,PDIQK,PGSMX,PGSPL(4),PHIMIX,PIFAC,
//       & PRSOF,PSPLT(2),PTRMS,PXRMS,QCDL3,QCDL5,QCDLAM,QDIQK,QFCH(16),QG,
//       & QSPAC,QV,SCABI,SWEIN,TMTOP,VFCH(16,2),VCKM(3,3),VGCUT,VQCUT,
//       & VPCUT,ZBINM,EFFMIN,OMHMIX,ET2MIX,PH3MIX,GCUTME,
//       & IOPREM,IPRINT,ISPAC,LRSUD,LWSUD,MODPDF(2),NBTRY,NCOLO,NCTRY,
//       & NDTRY,NETRY,NFLAV,NGSPL,NSTRU,NSTRY,NZBIN,IOP4JT(2),NPRFMT,
//       & AZSOFT,AZSPIN,CLDIR(2),HARDME,NOSPAC,PRNDEC,PRVTX,SOFTME,ZPRIME,
//       & PRNDEF,PRNTEX,PRNWEB

extern "C" {
    extern struct {
	double AFCH[2][16],ALPHEM,B1LIM,BETAF,BTCLM,CAFAC,CFFAC,
	    CLMAX,CLPOW,CLSMR[2],CSPEED,ENSOF,ETAMIX,F0MIX,F1MIX,F2MIX,GAMH,
	    GAMW,GAMZ,GAMZP,GEV2NB,H1MIX,PDIQK,PGSMX,PGSPL[4],PHIMIX,PIFAC,
	    PRSOF,PSPLT[2],PTRMS,PXRMS,QCDL3,QCDL5,QCDLAM,QDIQK,QFCH[16],QG,
	    QSPAC,QV,SCABI,SWEIN,TMTOP,VFCH[2][16],VCKM[3][3],VGCUT,VQCUT,
	    VPCUT,ZBINM,EFFMIN,OMHMIX,ET2MIX,PH3MIX,GCUTME;
	int IOPREM,IPRINT,ISPAC,LRSUD,LWSUD,MODPDF[2],NBTRY,NCOLO,NCTRY,
	    NDTRY,NETRY,NFLAV,NGSPL,NSTRU,NSTRY,NZBIN,IOP4JT[2],NPRFMT;
	int AZSOFT,AZSPIN,CLDIR[2],HARDME,NOSPAC,PRNDEC,PRVTX,SOFTME,
	    ZPRIME,PRNDEF,PRNTEX,PRNWEB; //Beware! in F77 these are logical
    } hwpram_;
}
#define hwpram hwpram_

//--------------------------------------------------------------------------
// HERWIG routines declaration

#define hwigin hwigin_  // initialise other common blocks
#define hwigup hwigup_  // initialise HepUP run common block
#define hwuinc hwuinc_  // compute parameter-dependent constants
#define hwusta hwusta_  // call hwusta to make any particle stable
#define hweini hweini_  // initialise elementary process
#define hwuine hwuine_  // initialise event
#define hwepro hwepro_  // generate HERWIG hard subprocess
#define hwupro hwupro_  // read USER hard subprocess from HepUP event common
#define hwbgen hwbgen_  // generate parton cascades
#define hwdhob hwdhob_  // do heavy object decays
#define hwcfor hwcfor_  // do cluster hadronization
#define hwcdec hwcdec_  // do cluster decay
#define hwdhad hwdhad_  // do unstable particle decays
#define hwdhvy hwdhvy_  // do heavy flavour decays
#define hwmevt hwmevt_  // add soft underlying event if needed 
#define hwufne hwufne_  // event generation completed, wrap up event .... 
#define hwefin hwefin_  // terminate elementary process

#define hwudpr hwudpr_  // prints out particle/decay properties
#define hwuepr hwuepr_  // prints out event data
#define hwupup hwupup_  // prints out HepEUP user common block event data
#define hwegup hwegup_  // terminal calculations to replace HWEFIN for HepUP
    extern "C" {
        void hwigin(void);
        void hwigup(void);
        void hwuinc(void);
        void hwusta(const char*,int);
        void hweini(void);
        void hwuine(void);
        void hwepro(void);
        void hwupro(void);
        void hwbgen(void);
        void hwdhob(void);
        void hwcfor(void);
        void hwcdec(void);
        void hwdhad(void);
        void hwdhvy(void);
        void hwmevt(void);
        void hwufne(void);
        void hwefin(void);
        void hwudpr(void);
        void hwuepr(void);
	void hwupup(void);
	void hwegup(void);
    }

//--------------------------------------------------------------------------
// HERWIG block data
// ( with gcc it works to initialize the block data by calling 
//   "hwudat();" at beginning. )

#define hwudat hwudat_
extern "C" {
    void hwudat(void);
}


#include "HepMC/GenCrossSection.h"

namespace HepMC {

/// calculate the Herwig cross section and statistical error
GenCrossSection getHerwigCrossSection(int ngen);

} // HepMC

#endif // Platform

#endif  // HERWIG_WRAPPER_H

//--------------------------------------------------------------------------
HepMC/GenVertex.h000644 000765 000024 00000046662 13310701701 014111 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HEPMC_GEN_VERTEX_H
#define HEPMC_GEN_VERTEX_H

//////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, September 1999, refer to:
// M. Dobbs and J.B. Hansen, "The HepMC C++ Monte Carlo Event Record for
// High Energy Physics", Computer Physics Communications (to be published).
//
// GenVertex within an event
// A vertex is indirectly (via particle "edges") linked to other 
//   vertices ("nodes") to form a composite "graph"
//////////////////////////////////////////////////////////////////////////

// --> HANDLE COMPILER INCONSISTENCIES
// This pre-compiler directive is included (2002-01-16) to allow compatibility
// with several compilers.
// Mar 27, 2004: HepMC is now standard compliant only. 
//   I've removed forward_iterator, and it will no longer compile on gcc < 3. 
#ifdef __SUNPRO_CC    // Solaris CC 5.2
#define NEED_SOLARIS_FRIEND_FEATURE
#endif // Platform

#include "HepMC/WeightContainer.h"
#include "HepMC/SimpleVector.h"
#include "HepMC/IteratorRange.h"
#include <iostream>
#include <iterator>
#include <vector>
#include <set>
#include <algorithm>
#include <cstddef>

namespace HepMC {

    class GenVertexParticleRange;
    class GenParticleProductionRange;
    class ConstGenParticleProductionRange;
    class GenParticleEndRange;
    class ConstGenParticleEndRange;

    class GenParticle;
    class GenEvent;

    //! GenVertex contains information about decay vertices.

    ///
    /// \class GenVertex 
    /// HepMC::GenVertex contains the position in space and time of a decay.
    /// It also contains lists of incoming and outgoing particles.  
    ///
    class GenVertex {

        /// print vertex information
	friend std::ostream& operator<<( std::ostream&, const GenVertex& );
	friend class GenEvent;

#ifdef NEED_SOLARIS_FRIEND_FEATURE
	// This bit of ugly code is only for CC-5.2 compiler. 
	// M.Dobbs 2002/02/19 
	// It is not needed by linux gcc, nor Windows Visual C++.
    public: 
	class vertex_iterator;
	friend class vertex_iterator;
	class particle_iterator;
	friend class particle_iterator;
#endif // NEED_SOLARIS_FRIEND_FEATURE

    public:
        /// default constructor
	GenVertex( const FourVector& position =FourVector(0,0,0,0),
		   int id = 0, 
		   const WeightContainer& weights = std::vector<double>() );
	GenVertex( const GenVertex& invertex );            //!< shallow copy
	virtual    ~GenVertex();

        void swap( GenVertex & other); //!< swap
	GenVertex& operator= ( const GenVertex& invertex ); //!< shallow
	bool       operator==( const GenVertex& a ) const; //!< equality
	bool       operator!=( const GenVertex& a ) const; //!< inequality
	void       print( std::ostream& ostr = std::cout ) const; //!< print vertex information

	double     check_momentum_conservation() const;//!< |Sum (three_mom_in-three_mom_out)|

	/// add incoming particle
	void       add_particle_in( GenParticle* inparticle );
	/// add outgoing particle
	void       add_particle_out( GenParticle* outparticle );
	/// remove_particle finds *particle in the in and/or out list and
	///  removes it from these lists ... it DOES NOT DELETE THE PARTICLE 
	///  or its relations. You could delete the particle too as follows:
	///      delete vtx->remove_particle( particle );
	GenParticle* remove_particle( GenParticle* particle ); //!< remove a particle

	operator    HepMC::FourVector() const; //!< conversion operator
	operator   HepMC::ThreeVector() const; //!< conversion operator

	////////////////////
	// access methods //
	////////////////////

	/// pointer to the event that owns this vertex
	GenEvent*               parent_event() const;
	/// vertex position
	ThreeVector             point3d() const;
	/// vertex position and time
	const FourVector &      position() const;
	/// set vertex position and time
	void                    set_position( const FourVector& position = FourVector(0,0,0,0) );
	/// we don't define what you use the id for -- but we imagine,
	/// for example it might code the meaning of the weights()
	int                     id() const;  //!< vertex ID
	void                    set_id( int id );  //!< set vertex ID

	///
	/// The barcode is the vertex's reference number, every vertex in the
	/// event has a unique barcode. Vertex barcodes are negative numbers,
	/// particle barcodes are positive numbers.
	/// 
	/// Please note that the barcodes are intended for internal use within 
	/// HepMC as a unique identifier for the particles and vertices.
	/// Using the barcode to encode extra information is an abuse of 
	/// the barcode data member and causes confusion among users. 
	/// 
	int                     barcode() const; //!< unique identifier

	/// In general there is no reason to "suggest_barcode"
	bool                    suggest_barcode( int the_bar_code );
	
	/// direct access to the weights container is allowed. 
	WeightContainer&        weights();
	/// const direct access to the weights container
	const WeightContainer&  weights() const;

	/// particle range
	GenVertexParticleRange particles( IteratorRange range = relatives );
	/// incoming particle range
	GenParticleProductionRange particles_in( GenParticle&, IteratorRange range = relatives );
	/// incoming particle range
	ConstGenParticleProductionRange particles_in( GenParticle const &, IteratorRange range = relatives ) const;
	/// outgoing particle range
	GenParticleEndRange particles_out( GenParticle&, IteratorRange range = relatives );
	/// outgoing particle range
	ConstGenParticleEndRange particles_out( GenParticle const &, IteratorRange range = relatives ) const;

	////////////////////
	// Iterators      // users should use prefer to use particle_iterator
	////////////////////    

	/// const iterator for incoming particles
	typedef std::vector<HepMC::GenParticle*>::const_iterator 
	particles_in_const_iterator;
	/// const iterator for outgoing particles
	typedef std::vector<HepMC::GenParticle*>::const_iterator 
	particles_out_const_iterator;
	/// begin iteration of incoming particles
	particles_in_const_iterator         particles_in_const_begin() const;
	/// end iteration of incoming particles
	particles_in_const_iterator         particles_in_const_end() const;
	/// begin iteration of outgoing particles
	particles_out_const_iterator        particles_out_const_begin() const;
	/// end iteration of outgoing particles
	particles_out_const_iterator        particles_out_const_end() const;
	/// number of incoming particles
	int                                 particles_in_size() const;
	/// number of outgoing particles
	int                                 particles_out_size() const;

    protected:
	//static unsigned int     counter(); //!< temporary for debugging

	/// only the GenEvent (friend) is allowed to set the parent_event,
	///  and barcode. It is done automatically anytime you add a 
	///  vertex to an event
	void                    set_parent_event_( GenEvent* evt ); //!< set parent event
	void                    set_barcode_( int the_bar_code ); //!< set identifier
	void                    change_parent_event_( GenEvent* evt ); //!< for use with swap

	/////////////////////////////
	// edge_iterator           // (protected - for internal use only)
	/////////////////////////////
	// If the user wants the functionality of the edge_iterator, he should
	// use particle_iterator with IteratorRange = family, parents, children
	//

	//!  edge iterator

	/// \class  edge_iterator
	/// iterate over the family of edges connected to m_vertex begins 
	/// with parents (incoming particles) then children (outgoing)
	/// This is not a recursive iterator ... it is a building block
	/// for the public iterators and is intended for internal use only.
	/// The acceptable Iterator Ranges are: family, parents, children
	class edge_iterator :
	  public std::iterator<std::forward_iterator_tag,HepMC::GenParticle*,ptrdiff_t>{
	public:
	    edge_iterator();
	    /// used to set limits on the iteration
	    edge_iterator( const GenVertex& vtx, IteratorRange range =family );
	    /// copy
	    edge_iterator( const edge_iterator& p );
	    virtual        ~edge_iterator();
	    /// make a copy
	    edge_iterator& operator=( const edge_iterator& p );
	    /// return a pointer to a particle
	    GenParticle*      operator*(void) const;
	    /// Pre-fix increment 
	    edge_iterator& operator++(void); // Pre-fix increment 
	    /// Post-fix increment
	    edge_iterator  operator++(int);   // Post-fix increment
	    /// equality
	    bool           operator==( const edge_iterator& a ) const;
	    /// inequality
	    bool           operator!=( const edge_iterator& a ) const;
	    /// true if parent of root vtx
	    bool           is_parent() const;
	    /// true if child of root vtx
	    bool           is_child() const;
	    /// root vertex of this iteration
	    const GenVertex*  vertex_root() const;
	private:
	    /// Pre-fix increment -- is not allowed
	    edge_iterator& operator--(void);
	    /// Post-fix increment -- is not allowed
	    edge_iterator  operator--(int);
	private:
	    const GenVertex*  m_vertex;
	    IteratorRange  m_range;
	    std::vector<HepMC::GenParticle*>::const_iterator m_set_iter;
	    bool           m_is_inparticle_iter;
	    bool           m_is_past_end;
	};
	friend class edge_iterator;
	/// size
	int              edges_size( IteratorRange range = family ) const;
	/// begin range
	edge_iterator    edges_begin( IteratorRange range = family) const;
	/// end range
	edge_iterator    edges_end( IteratorRange /* dummy_range */ ) const;

    public:
	///////////////////////////////
	// vertex_iterator           //
	///////////////////////////////

	//!  vertex iterator

	/// \class  vertex_iterator
	/// Iterates over all vertices connected via a graph to this vertex.
	/// this is made friend to that it can access protected edge
	/// iterator the range can be IteratorRange= ( parents, children, 
	/// family, ancestors, descendants, relatives )
	/// example for range=descendants the iterator 
	/// will return all vertices
	/// which are children (connected by an outgoing particle edge),
	/// grandchildren, great-grandchildren, etc. of this vertex
	/// In all cases the iterator always returns this vertex
	/// (returned last).
	/// The algorithm is accomplished by converting the graph to a tree
	/// (by "chopping" the edges connecting to an already visited
	/// vertex) and returning the vertices in POST ORDER traversal.
	///
	class vertex_iterator :
	  public std::iterator<std::forward_iterator_tag,HepMC::GenVertex*,ptrdiff_t>{
	public:
	    vertex_iterator();
	    /// used to set limits on the iteration
	    vertex_iterator( GenVertex& vtx_root, IteratorRange range );
	    /// next constructor is intended for internal use only
	    vertex_iterator( GenVertex& vtx_root, IteratorRange range,
			     std::set<const HepMC::GenVertex*>& visited_vertices );
            /// copy
	    vertex_iterator( const vertex_iterator& v_iter );
	    virtual             ~vertex_iterator();
	    /// make a copy
	    vertex_iterator&    operator=( const vertex_iterator& );
	    /// return a pointer to a vertex
	    GenVertex*          operator*(void) const;
	    /// Pre-fix increment 
	    vertex_iterator&    operator++(void);  //Pre-fix increment 
	    /// Post-fix increment
	    vertex_iterator     operator++(int);   //Post-fix increment
	    /// equality
	    bool                operator==( const vertex_iterator& ) const;
	    /// inequality
	    bool                operator!=( const vertex_iterator& ) const;
	    /// vertex that this iterator begins from
	    GenVertex*          vertex_root() const;
	    /// iterator range
	    IteratorRange       range() const;
	    /// intended for internal use only.
	    void                copy_with_own_set( const vertex_iterator& 
						   v_iter,
						   std::set<const HepMC::GenVertex*>& 
						   visited_vertices );

	protected:                  // intended for internal use only
	    /// non-null if recursive iter. created
	    GenVertex* follow_edge_(); 
	    /// copy recursive iterator
	    void    copy_recursive_iterator_( const vertex_iterator* 
					      recursive_v_iter );
	private:
	    /// Pre-fix increment -- is not allowed 
	    vertex_iterator&    operator--(void);
	    /// Post-fix increment -- is not allowed 
	    vertex_iterator     operator--(int);

	private:
	    GenVertex*       m_vertex;   // the vertex associated to this iter
	    IteratorRange    m_range;
	    std::set<const HepMC::GenVertex*>* m_visited_vertices;
	    bool             m_it_owns_set;  // true if it is responsible for 
	                                     // deleting the visited vertex set
	    edge_iterator    m_edge; // particle edge pointing to return vtx
	    vertex_iterator* m_recursive_iterator;
	};	
	friend class vertex_iterator;
	/// begin vertex range
	vertex_iterator     vertices_begin( IteratorRange range = relatives );
	/// end vertex range
	vertex_iterator     vertices_end( IteratorRange /* dummy_range */ );
 
    public:
	///////////////////////////////
	// particle_iterator         //
	///////////////////////////////

	//!  particle iterator

	/// \class  particle_iterator
	/// Iterates over all particles connected via a graph.
	/// by iterating through all vertices in the m_range. For each
	/// vertex it returns orphaned parent particles 
	/// (i.e. parents without production vertices) 
	/// then children ... in this way each particle is associated
	/// to exactly one vertex and so it is returned exactly once.
	/// Is made friend so that it can access protected edge iterator
 	class particle_iterator :
	  public std::iterator<std::forward_iterator_tag,GenParticle*,ptrdiff_t>{
	public:
	    particle_iterator();
	    /// used to set limits on the iteration
	    particle_iterator( GenVertex& vertex_root, IteratorRange range );
	    /// copy
	    particle_iterator( const particle_iterator& );
	    virtual             ~particle_iterator();
	    /// make a copy
	    particle_iterator&  operator=( const particle_iterator& );
	    /// return a pointer to a particle
	    GenParticle*        operator*(void) const;
	    /// Pre-fix increment 
	    particle_iterator&  operator++(void); 
	    /// Post-fix increment
	    particle_iterator   operator++(int); 
	    /// equality
	    bool                operator==( const particle_iterator& ) const;
	    /// inequality
	    bool                operator!=( const particle_iterator& ) const;
	protected:
	    GenParticle*        advance_to_first_(); //!< "first" particle
	private:
	    vertex_iterator     m_vertex_iterator;
	    edge_iterator       m_edge;     // points to the return
	};
	friend class particle_iterator;
	/// begin particle range
	particle_iterator       particles_begin( IteratorRange range 
						 = relatives );
	/// end particle range
	particle_iterator       particles_end( IteratorRange 
					       /* dummy_range */ );

	////////////////////////////////////////////////
    protected: 
        /// for internal use only
	void delete_adopted_particles();
	/// for internal use only - remove particle from incoming list
	void remove_particle_in( GenParticle* );
	/// for internal use only - remove particle from outgoing list
	void remove_particle_out( GenParticle* );
	/// scale the position vector
        /// this method is only for use by GenEvent
	void convert_position( const double& );

    private: // GenVertex data members
	FourVector              m_position;      //4-vec of vertex [mm]
	std::vector<HepMC::GenParticle*>  m_particles_in;  //all incoming particles
	std::vector<HepMC::GenParticle*>  m_particles_out; //all outgoing particles
	int                  m_id;
	WeightContainer      m_weights;       // weights for this vtx
	GenEvent*            m_event;
	int                  m_barcode;   // unique identifier in the event

	//static unsigned int  s_counter;
    };  

    ////////////////////////////
    // INLINES access methods //
    ////////////////////////////

    inline GenVertex::operator HepMC::FourVector() const { return position(); }

    inline GenVertex::operator HepMC::ThreeVector() const { return point3d(); }

    inline const FourVector & GenVertex::position() const { return m_position; }

    inline GenEvent* GenVertex::parent_event() const { return m_event; }

    inline ThreeVector GenVertex::point3d() const { 
	return ThreeVector(m_position.x(),m_position.y(),m_position.z()); 
    }

    inline int GenVertex::id() const { return m_id; }

    inline int  GenVertex::barcode() const { return m_barcode; }
    inline void GenVertex::set_barcode_( int bc ) { m_barcode = bc; }

    inline WeightContainer& GenVertex::weights() { return m_weights; }

    inline const WeightContainer& GenVertex::weights() const 
    { return m_weights; }

    inline void GenVertex::set_position( const FourVector& pos ) {
	m_position = pos;
    }

    inline void GenVertex::set_id( int pid ) { m_id = pid; }

    //////////////
    // INLINES  //
    //////////////

    inline GenVertex::particles_in_const_iterator 
    GenVertex::particles_in_const_begin() const { 
	return m_particles_in.begin(); 
    }

    inline GenVertex::particles_in_const_iterator 
    GenVertex::particles_in_const_end() const { 
	return m_particles_in.end(); 
    }

    inline GenVertex::particles_out_const_iterator 
    GenVertex::particles_out_const_begin() const { 
	return m_particles_out.begin();
    }

    inline GenVertex::particles_out_const_iterator 
    GenVertex::particles_out_const_end() const {	
	return m_particles_out.end(); 
    }

    inline int GenVertex::particles_in_size() const {
	return m_particles_in.size(); 
    }

    inline int GenVertex::particles_out_size() const {
	return m_particles_out.size(); 
    }	

    inline bool GenVertex::edge_iterator::operator==( 
	const edge_iterator& a ) const { 
	return **this == *a; 
    }

    inline bool GenVertex::edge_iterator::operator!=(
	const edge_iterator& a ) const { 
	return !(**this == *a); 
    }

    inline const GenVertex* GenVertex::edge_iterator::vertex_root() const {
	return m_vertex;
    }

    inline GenVertex::edge_iterator GenVertex::edges_begin( IteratorRange 
						      range ) const {
	return GenVertex::edge_iterator(*this, range);
    }

    inline GenVertex::edge_iterator GenVertex::edges_end( IteratorRange 
						    /* dummy_range */ ) const {
	return GenVertex::edge_iterator();
    }

    inline bool GenVertex::vertex_iterator::operator==( 
	const vertex_iterator& a ) const {
	return **this == *a; 
    }

    inline bool GenVertex::vertex_iterator::operator!=( 
	const vertex_iterator& a ) const {
	return !(**this == *a); 
    }

    inline GenVertex* GenVertex::vertex_iterator::vertex_root() const {
	return m_vertex; 
    }

    inline IteratorRange GenVertex::vertex_iterator::range() const {
	return m_range; 
    }

    inline GenVertex::vertex_iterator GenVertex::vertices_begin( 
	IteratorRange range ){
	// this is not const because the it could return itself
	return vertex_iterator( *this, range );
    }

    inline GenVertex::vertex_iterator GenVertex::vertices_end( 
	IteratorRange /* dummy_range */ ) {
	return vertex_iterator();
    }

    inline bool GenVertex::particle_iterator::operator==( 
	const particle_iterator& a ) const {
	return **this == *a; 
    }

    inline bool GenVertex::particle_iterator::operator!=( 
	const particle_iterator& a ) const {
	return !(**this == *a); 
    }

    inline GenVertex::particle_iterator GenVertex::particles_begin( 
	IteratorRange range ) {
	return particle_iterator( *this, range );
    }

    inline GenVertex::particle_iterator GenVertex::particles_end(
	IteratorRange /* dummy_range */ ){
	return particle_iterator();
    }

} // HepMC

#endif  // HEPMC_GEN_VERTEX_H
//--------------------------------------------------------------------------




HepMC/dummy.F000644 000765 000024 00000002601 13310701701 013254 0ustar00sa639staff000000 000000       BLOCK DATA DUMMY
C     Initialize the HEPEVT COMMON block
C     This helps to resolve the hepevt_ symbol when we link the library

#include "stdhep.inc"

      DATA nevhep/0/                   ! serial number
      DATA nhep/0/                     ! number of particles
      DATA (isthep(I),I=1,nmxhep)/nmxhep*0/           ! status code
      DATA (idhep(I),I=1,nmxhep)/nmxhep*0/            ! particle ident KF
      DATA (jmohep(1,I),I=1,nmxhep)/nmxhep*0/         ! parent particles
      DATA (jmohep(2,I),I=1,nmxhep)/nmxhep*0/         ! parent particles
      DATA (jdahep(1,I),I=1,nmxhep)/nmxhep*0/         ! childreen particles
      DATA (jdahep(2,I),I=1,nmxhep)/nmxhep*0/         ! childreen particles
      DATA (phep(1,I),I=1,nmxhep)/nmxhep*0/           ! four-momentum, mass [GeV]
      DATA (phep(2,I),I=1,nmxhep)/nmxhep*0/           ! four-momentum, mass [GeV]
      DATA (phep(3,I),I=1,nmxhep)/nmxhep*0/           ! four-momentum, mass [GeV]
      DATA (phep(4,I),I=1,nmxhep)/nmxhep*0/           ! four-momentum, mass [GeV]
      DATA (phep(5,I),I=1,nmxhep)/nmxhep*0/           ! four-momentum, mass [GeV]
      DATA (vhep(1,I),I=1,nmxhep)/nmxhep*0/           ! vertex [mm]
      DATA (vhep(2,I),I=1,nmxhep)/nmxhep*0/           ! vertex [mm]
      DATA (vhep(3,I),I=1,nmxhep)/nmxhep*0/           ! vertex [mm]
      DATA (vhep(4,I),I=1,nmxhep)/nmxhep*0/           ! vertex [mm]

      END
HepMC/CompareGenEvent.h000644 000765 000024 00000001666 13310701701 015217 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HEPMC_COMPARE_GENEVENT_H
#define HEPMC_COMPARE_GENEVENT_H

/////////////////////////////////////////////////////////////////////////
// CompareGenEvent.h
//
// garren@fnal.gov, January 2008
// Free functions used to compare two copies of GenEvent
//////////////////////////////////////////////////////////////////////////
//

#include <iostream>

#include "HepMC/GenEvent.h"

namespace HepMC {

bool compareGenEvent( GenEvent*, GenEvent* );
bool compareSignalProcessVertex( GenEvent*, GenEvent* );
bool compareBeamParticles( GenEvent*, GenEvent* );
bool compareWeights( GenEvent*, GenEvent* );
bool compareVertices( GenEvent*, GenEvent* );
bool compareParticles( GenEvent*, GenEvent* );
bool compareVertex( GenVertex* v1, GenVertex* v2 );

} // HepMC

#endif  // HEPMC_COMPARE_GENEVENT_H
//--------------------------------------------------------------------------
HepMC/Units.h000644 000765 000024 00000002436 13310701701 013273 0ustar00sa639staff000000 000000 #ifndef HEPMC_UNITS_H
#define HEPMC_UNITS_H

//--------------------------------------------------------------------------
// Units.h
// Author:  A. Buckley, D. Grellscheid
//
// units used by a GenEvent
// The default units are set by a configure switch at compile time in Units.cc.
//--------------------------------------------------------------------------

#include <iostream>
#include <string>

namespace HepMC {

  ///
  /// \namespace Units
  /// Allow units to be specified within HepMC.
  /// The default units are set at compile time. 
  ///
  namespace Units {

    // Convention: if both types are passed, MomentumUnit always goes first.
    enum MomentumUnit { MEV, GEV };	//!< momentum units
    enum LengthUnit   { MM, CM };	//!< position units
    
    LengthUnit   default_length_unit();		//!< default unit is defined by configure
    MomentumUnit default_momentum_unit();	//!< default unit is defined by configure

    // helper functions
    std::string name( MomentumUnit );	//!< convert enum to string
    std::string name( LengthUnit );	//!< convert enum to string

    /// scaling factor relative to MeV
    double conversion_factor( MomentumUnit from, MomentumUnit to ); 
    double conversion_factor( LengthUnit from, LengthUnit to );

  }	// Units
}	// HepMC

#endif // HEPMC_UNITS_H
HepMC/IO_GenEvent.h000644 000765 000024 00000011765 13310701701 014300 0ustar00sa639staff000000 000000 //--------------------------------------------------------------------------
#ifndef HEPMC_IO_GENEVENT_H
#define HEPMC_IO_GENEVENT_H

//////////////////////////////////////////////////////////////////////////
// garren@fnal.gov, July 2007
// with input from Gavin Salam, salam@lpthe.jussieu.fr
//
// event input/output in ascii format for machine reading
// This class persists all information found in a GenEvent
//////////////////////////////////////////////////////////////////////////

#include <fstream>
#include <string>
#include <map>
#include <vector>
#include "HepMC/IO_BaseClass.h"
#include "HepMC/IO_Exception.h"
#include "HepMC/Units.h"

namespace HepMC {

class GenEvent;
class GenVertex;
class GenParticle;
class HeavyIon;
class PdfInfo;

//! IO_GenEvent also deals with HeavyIon and PdfInfo 

///
/// \class  IO_GenEvent
/// event input/output in ascii format for machine reading
/// extended format contains HeavyIon and PdfInfo classes
///
/// Strategy for reading or writing events using iostreams
/// When instantiating with a file name, the mode of file to be created 
///  must be specified. Options are:
///      std::ios::in     open file for input 
///      std::ios::out    open file for output
///      std::ios::trunc  erase old file when opening (i.e. ios::out|ios::trunc
///                    removes oldfile, and creates a new one for output )
///      std::ios::app    append output to end of file
///  for the purposes of this class, simultaneous input and output mode 
///  ( std::ios::in | std::ios::out ) is not allowed.
/// 
/// Event listings are preceded by the key:
///  "HepMC::IO_GenEvent-START_EVENT_LISTING\n"
///  and terminated by the key:
///  "HepMC::IO_GenEvent-END_EVENT_LISTING\n"
/// GenParticle Data tables are preceded by the key:
///  "HepMC::IO_GenEvent-START_PARTICLE_DATA\n"
///  and terminated by the key:
///  "HepMC::IO_GenEvent-END_PARTICLE_DATA\n"
/// Comments are allowed. They need not be preceded by anything, though if
///  a comment is written using write_comment( const string ) then it will be
///  preceded by "HepMC::IO_GenEvent-COMMENT\n"
/// Each event, vertex, particle, particle data, heavy ion, or pdf info line
///  is preceded by   "E ","V ","P ","D ","H ","F "    respectively.
/// Comments may appear anywhere in the file -- so long as they do not contain
///  any of the start/stop keys.
///
class IO_GenEvent : public IO_BaseClass {
public:
    /// constructor requiring a file name and std::ios mode
    IO_GenEvent( const std::string& filename="IO_GenEvent.dat", 
	      std::ios::openmode mode=std::ios::out );
    /// constructor requiring an input stream
    IO_GenEvent( std::istream & );
    /// constructor requiring an output stream
    IO_GenEvent( std::ostream & );
    virtual       ~IO_GenEvent();

    /// write this event
    void          write_event( const GenEvent* evt );
    /// get the next event
    bool          fill_next_event( GenEvent* evt );
    /// insert a comment directly into the output file --- normally you
    ///  only want to do this at the beginning or end of the file. All
    ///  comments are preceded with "HepMC::IO_GenEvent-COMMENT\n"
    void          write_comment( const std::string comment );

    int           rdstate() const;  //!< check the state of the IO stream
    void          clear();  //!< clear the IO stream

    /// write to ostr
    void          print( std::ostream& ostr = std::cout ) const;

    /// needed when reading a file without units if those units are 
    /// different than the declared default units 
    /// (e.g., the default units are MeV, but the file was written with GeV)
    /// This method is not necessary if the units are written in the file
    void use_input_units( Units::MomentumUnit, Units::LengthUnit );

    /// set output precision
    /// The default precision is 16.
    void precision( int );
	
    /// integer (enum) associated with read error
    int           error_type()    const;
    /// the read error message string
    const std::string & error_message() const;

private: // use of copy constructor is not allowed
    IO_GenEvent( const IO_GenEvent& ) : IO_BaseClass() {}

private: // data members
    std::ios::openmode  m_mode;
    std::fstream        m_file;
    std::ostream *      m_ostr;
    std::istream *      m_istr;
    std::ios *          m_iostr;
    bool                m_have_file;
    IO_Exception::ErrorType m_error_type;
    std::string         m_error_message;

};

//////////////
// Inlines  //
//////////////

inline int  IO_GenEvent::rdstate() const { 
    int state;
    if( m_istr ) {
	state =  (int)m_istr->rdstate();
    } else {
	state =  (int)m_ostr->rdstate();
    }
    return state; 
}

inline void IO_GenEvent::clear() { 
    if( m_istr ) {
	m_istr->clear();
    } else {
	m_ostr->clear();
    }
}

inline int IO_GenEvent::error_type() const {
    return m_error_type;
}

inline const std::string & IO_GenEvent::error_message() const {
    return m_error_message;
}

} // HepMC

#endif  // HEPMC_IO_GENEVENT_H
//--------------------------------------------------------------------------
HepMC/SimpleVector.h000644 000765 000024 00000014334 13310701701 014605 0ustar00sa639staff000000 000000 //////////////////////////////////////////////////////////////////////////
// SimpleVector.h
//////////////////////////////////////////////////////////////////////////
#ifndef  HEPMC_SIMPLEVECTOR_H
#define  HEPMC_SIMPLEVECTOR_H

//////////////////////////////////////////////////////////////////////////
// garren@fnal.gov, July 2006
//
// This header provides a place to hold the doubles which are part of one of 
// three types of physics vectors:
//    momentum 4 vector 
//    position or displacement 4 vector
//    position or displacement 3 vector
//
// For compatibility with existing code, 
// the basic expected geometrical access methods are povided 
// Also, both FourVector and ThreeVector have a templated constructor that will 
// take another vector (HepLorentzVector, GenVector, ...)
//    --> this vector must have the following methods: x(), y(), z()
//    -->  FourVector also requires the t() method
//
//////////////////////////////////////////////////////////////////////////


#include "HepMC/enable_if.h"
#include "HepMC/is_arithmetic.h"


namespace HepMC {

//! FourVector is a simple representation of a physics 4 vector

///
/// \class  FourVector
/// For compatibility with existing code, 
/// the basic expected geometrical access methods are povided.
/// Also, there is a templated constructor that will 
/// take another vector (HepLorentzVector, GenVector, ...)
/// which must have the following methods: x(), y(), z(), t().
///
class FourVector {

public:

  /// constructor requiring at least x, y, and z
  FourVector( double xin, double yin, double zin, double tin=0) 
  : m_x(xin), m_y(yin), m_z(zin), m_t(tin) {}

  /// constructor requiring only t 
  FourVector(double tin)
  : m_x(0), m_y(0), m_z(0), m_t(tin) {}

  FourVector() 
  : m_x(0), m_y(0), m_z(0), m_t(0) {}

  /// templated constructor
  /// this is used ONLY if T is not arithmetic
  template <class T >
  FourVector( const T& v,
         typename detail::disable_if< detail::is_arithmetic<T>::value, void >::type * = 0 )
  : m_x(v.x()), m_y(v.y()), m_z(v.z()), m_t(v.t()) {}

  /// copy constructor
  FourVector(const FourVector & v)
  : m_x(v.x()), m_y(v.y()), m_z(v.z()), m_t(v.t()) {}

  void swap( FourVector & other );  //!< swap

  double px() const { return m_x; }  //!< return px
  double py() const { return m_y; }  //!< return py
  double pz() const { return m_z; }  //!< return pz
  double e()  const { return m_t; }  //!< return E

  double x() const { return m_x; }  //!< return x
  double y() const { return m_y; }  //!< return y
  double z() const { return m_z; }  //!< return z
  double t() const { return m_t; }  //!< return t

  double m2() const;  //!< Invariant mass squared.
  double m() const;   //!< Invariant mass. If m2() is negative then -sqrt(-m2()) is returned.

  double perp2() const;  //!< Transverse component of the spatial vector squared.
  double perp() const;   //!< Transverse component of the spatial vector (R in cylindrical system).

  // Get spatial vector components in spherical coordinate system.
  double theta() const;  //!< The polar angle.
  double phi() const;  //!< The azimuth angle.
  double rho() const;  //!< spatial vector component magnitude

  FourVector & operator = (const FourVector &); //!< make a copy

  bool operator == (const FourVector &) const; //!< equality
  bool operator != (const FourVector &) const; //!< inequality

  double pseudoRapidity() const;  //!< Returns the pseudo-rapidity, i.e. -ln(tan(theta/2))
  double eta() const;             //!< Pseudorapidity (of the space part)

  /// set x, y, z, and t
  void set        (double x, double y, double z, double  t);

  void setX(double xin) { m_x=xin; }  //!< set x
  void setY(double yin) { m_y=yin; }  //!< set y
  void setZ(double zin) { m_z=zin; }  //!< set z
  void setT(double tin) { m_t=tin; }  //!< set t

  void setPx(double xin) { m_x=xin; }  //!< set px
  void setPy(double yin) { m_y=yin; }  //!< set py
  void setPz(double zin) { m_z=zin; }  //!< set pz
  void setE(double tin)  { m_t=tin; }  //!< set E

private:

  double m_x;
  double m_y;
  double m_z;
  double m_t;
  
};

//! ThreeVector is a simple representation of a position or displacement 3 vector

///
/// \class  ThreeVector
/// For compatibility with existing code, 
/// the basic expected geometrical access methods are povided.
/// Also, there is a templated constructor that will 
/// take another vector (HepLorentzVector, GenVector, ...)
/// which must have the following methods: x(), y(), z().
///
class ThreeVector {

public:

  /// construct using x, y, and z (only x is required)
  ThreeVector( double xin, double yin =0, double zin =0 ) 
  : m_x(xin), m_y(yin), m_z(zin) {}

  ThreeVector( ) 
  : m_x(0), m_y(0), m_z(0) {}
  
  /// templated constructor
  /// this is used ONLY if T is not arithmetic
  template <class T >
  ThreeVector( const T& v,
         typename detail::disable_if< detail::is_arithmetic<T>::value, void >::type * = 0 )
  : m_x(v.x()), m_y(v.y()), m_z(v.z()) {}

  /// copy constructor
  ThreeVector(const ThreeVector & v)
  : m_x(v.x()), m_y(v.y()), m_z(v.z()) {}

  void swap( ThreeVector & other );  //!< swap

  double x() const { return m_x; }  //!< return x
  double y() const { return m_y; }  //!< return y
  double z() const { return m_z; }  //!< return z

  void setX(double xin) { m_x=xin; }  //!< set x
  void setY(double yin) { m_y=yin; }  //!< set y
  void setZ(double zin) { m_z=zin; }  //!< set z
  void set( double x, double y, double z);   //!< set x, y, and z

  double phi()   const;  //!< The azimuth angle.
  double theta() const;  //!< The polar angle.
  double r()     const;  //!< The magnitude

  void setPhi(double);  //!< Set phi keeping magnitude and theta constant (BaBar).
  void setTheta(double);  //!< Set theta keeping magnitude and phi constant (BaBar).

  double perp2() const;  //!< The transverse component squared (rho^2 in cylindrical coordinate system).
  double perp() const;  //!< The transverse component (rho in cylindrical coordinate system).

  ThreeVector & operator = (const ThreeVector &); //!< make a copy

  bool operator == (const ThreeVector &) const; //!< equality
  bool operator != (const ThreeVector &) const; //!< inequality

private:

  double m_x;
  double m_y;
  double m_z;

};  


} // HepMC

#include "HepMC/SimpleVector.icc"

#endif  // HEPMC_SIMPLEVECTOR_H

HepMC/SearchVector.h000644 000765 000024 00000001427 13310701701 014560 0ustar00sa639staff000000 000000 #ifndef SearchVector_h
#define SearchVector_h
// ----------------------------------------------------------------------
//
// SearchVector.h
// Author: Lynn Garren
//
// Utilities to search std::vector<GenParticle*> a GenParticle instance
// ----------------------------------------------------------------------

#include "HepMC/GenVertex.h"
#include "HepMC/GenParticle.h"

namespace HepMC {
  
/// returns true if it cannot find GenParticle* in the vector
bool not_in_vector( std::vector<HepMC::GenParticle*>*, GenParticle* );

/// Returns the index of a GenParticle* within a vector.
/// Returns -1 if GenParticle* is not in the vector.
std::vector<HepMC::GenParticle*>::iterator already_in_vector( std::vector<HepMC::GenParticle*>*, GenParticle* );

}	// HepMC

#endif // SearchVector_h