EPOS_Wrapper.h

Go to the documentation of this file.

//--------------------------------------------------------------------------
//THIS IS A BRUTAL COPY OF HEPEVT_Wrapper from HEPMC
//We need it because the EPOS generator needs a largeer version of HEPEVT to store the event
#ifndef EPOS_EntriesAllocation
#define EPOS_EntriesAllocation 99900
#endif  // EPOS_EntriesAllocation

//--------------------------------------------------------------------------
#ifndef HEPMC_EPOS_COMMON_H
#define HEPMC_EPOS_COMMON_H
//
//      PARAMETER (NMXHEP=2000)
//      COMMON/HEPCOM/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 EPOS 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.
//
// EPOS_EntriesAllocation is the maximum size of the EPOS common block
//   that can be interfaced.
//   It is NOT the actual size of the EPOS common used in each
//   individual application. The actual size can be changed on
//   the fly using EPOS_Wrapper::set_max_number_entries().
// Thus EPOS_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*)&hepcom.data[0] )
//      takes the memory address of the first byte in EPOS,
//      interprets it as an integer pointer,
//      and dereferences the pointer.
//      i.e. it returns an integer corresponding to nevhep
//

#include <cctype>

const unsigned int epos_bytes_allocation =
    sizeof(long int) * (2 + 6 * EPOS_EntriesAllocation) + sizeof(double) * (9 * EPOS_EntriesAllocation);

#ifdef _WIN32  // Platform: Windows MS Visual C++
struct HEPCOM_DEF {
  char data[epos_bytes_allocation];
};
extern "C" HEPCOM_DEF HEPCOM;
#define hepcom HEPCOM

#else
extern "C" {
extern struct { char data[epos_bytes_allocation]; } hepcom_;
}
#define hepcom hepcom_

#endif  // Platform

#endif  // HEPMC_EPOS_COMMON_H

//--------------------------------------------------------------------------
#ifndef HEPMC_EPOS_WRAPPER_H
#define HEPMC_EPOS_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 EPOS 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 EPOS 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 EPOS {


  class EPOS_Wrapper {
  public:
    static void print_hepcom(std::ostream& ostr = std::cout);
    static void print_hepcom_particle(int index, std::ostream& ostr = std::cout);
    static void zero_everything();

    // Access Methods //
    static int event_number();              
    static int number_entries();            
    static int status(int index);           
    static int id(int index);               
    static int first_parent(int index);     
    static int last_parent(int index);      
    static int number_parents(int index);   
    static int first_child(int index);      
    static int last_child(int index);       
    static int number_children(int index);  
    static double px(int index);            
    static double py(int index);            
    static double pz(int index);            
    static double e(int index);             
    static double m(int index);             
    static double x(int index);             
    static double y(int index);             
    static double z(int index);             
    static double t(int index);             

    // Set Methods    //

    static void set_event_number(int evtno);
    static void set_number_entries(int noentries);
    static void set_status(int index, int status);
    static void set_id(int index, int id);
    static void set_parents(int index, int firstparent, int lastparent);
    static void set_children(int index, int firstchild, int lastchild);
    static void set_momentum(int index, double px, double py, double pz, double e);
    static void set_mass(int index, double mass);
    static void set_position(int index, double x, double y, double z, double t);
    // EPOS Floorplan //
    static unsigned int sizeof_int();                  
    static unsigned int sizeof_real();                 
    static int max_number_entries();                   
    static void set_sizeof_int(unsigned int);          
    static void set_sizeof_real(unsigned int);         
    static void set_max_number_entries(unsigned int);  

  protected:
    static double byte_num_to_double(unsigned int);
    static int byte_num_to_int(unsigned int);
    static void write_byte_num(double, unsigned int);
    static void write_byte_num(int, unsigned int);
    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;
  };

  // EPOS Floorplan Inlines //
  inline unsigned int EPOS_Wrapper::sizeof_int() { return s_sizeof_int; }

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

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

  inline void EPOS_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 EPOS_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 EPOS_Wrapper::set_max_number_entries(unsigned int size) { s_max_number_entries = size; }

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

  inline int EPOS_Wrapper::byte_num_to_int(unsigned int b) {
    if (b >= epos_bytes_allocation)
      std::cerr << "EPOS_Wrapper: requested hepcom data exceeds allocation" << std::endl;
    if (s_sizeof_int == sizeof(short int)) {
      short int* myshortint = (short int*)&hepcom.data[b];
      return (int)(*myshortint);
    } else if (s_sizeof_int == sizeof(long int)) {
      long int* mylongint = (long int*)&hepcom.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*)&hepcom.data[b];
      return (*myint);
    } else {
      std::cerr << "EPOS_Wrapper: illegal integer number length." << s_sizeof_int << std::endl;
    }
    return 0;
  }

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

  inline void EPOS_Wrapper::write_byte_num(int in, unsigned int b) {
    if (b >= epos_bytes_allocation)
      std::cerr << "EPOS_Wrapper: requested hepcom data exceeds allocation" << std::endl;
    if (s_sizeof_int == sizeof(short int)) {
      short int* myshortint = (short int*)&hepcom.data[b];
      (*myshortint) = (short int)in;
    } else if (s_sizeof_int == sizeof(long int)) {
      long int* mylongint = (long int*)&hepcom.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*)&hepcom.data[b];
      (*myint) = (int)in;
    } else {
      std::cerr << "EPOS_Wrapper: illegal integer number length." << s_sizeof_int << std::endl;
    }
  }

  // INLINES  //

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

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

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

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

  inline int EPOS_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 EPOS_Wrapper::last_parent(int index) {
    // Returns the Index of the LAST parent in the EPOS 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 EPOS_Wrapper::number_parents(int index) {
    int firstparent = first_parent(index);
    return (firstparent > 0) ? (1 + last_parent(index) - firstparent) : 0;
  }

  inline int EPOS_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 EPOS_Wrapper::last_child(int index) {
    // Returns the Index of the LAST child in the EPOS 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 EPOS_Wrapper::number_children(int index) {
    int firstchild = first_child(index);
    return (firstchild > 0) ? (1 + last_child(index) - firstchild) : 0;
  }

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

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

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

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

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

  inline double EPOS_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 EPOS_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 EPOS_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 EPOS_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 EPOS_Wrapper::set_event_number(int evtno) { write_byte_num(evtno, 0); }

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

  inline void EPOS_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 EPOS_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 EPOS_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 EPOS_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 EPOS_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 EPOS_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 EPOS_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());
  }

  inline void EPOS_Wrapper::zero_everything() {
    set_event_number(0);
    set_number_entries(0);
    for (int i = 1; i <= max_number_entries(); ++i) {
      set_status(i, 0);
      set_id(i, 0);
      set_parents(i, 0, 0);
      set_children(i, 0, 0);
      set_momentum(i, 0, 0, 0, 0);
      set_mass(i, 0);
      set_position(i, 0, 0, 0, 0);
    }
  }

  inline void EPOS_Wrapper::print_hepcom(std::ostream& ostr) {
    ostr << "________________________________________"
         << "________________________________________" << std::endl;
    ostr << "***** HEPEVT Common Event#: " << event_number() << ", " << number_entries() << " particles (max "
         << max_number_entries() << ") *****";
    ostr << sizeof_int() << "-byte integers, " << sizeof_real() << "-byte floating point numbers, "
         << max_number_entries() << "-allocated entries." << std::endl;
    print_legend(ostr);
    ostr << "________________________________________"
         << "________________________________________" << std::endl;
    for (int i = 1; i <= number_entries(); ++i) {
      print_hepcom_particle(i, ostr);
    }
    ostr << "________________________________________"
         << "________________________________________" << std::endl;
  }

  inline void EPOS_Wrapper::print_hepcom_particle(int i, std::ostream& ostr) {
    char outline[81];
    sprintf(outline,
            "%4d %+4d %4d %4d    (%9.3g, %9.3g, %9.3g, %9.3g, %9.3g)",
            i,
            status(i),
            first_parent(i),
            first_child(i),
            px(i),
            py(i),
            pz(i),
            e(i),
            m(i));
    ostr << outline << "\n";
    sprintf(outline,
            "%+9d %4d %4d    (%9.3g, %9.3g, %9.3g, %9.3g)",
            // old version was:" (%+9.2e, %+9.2e, %+9.2e, %+9.2e)"
            id(i),
            last_parent(i),
            last_child(i),
            x(i),
            y(i),
            z(i),
            t(i));
    ostr << outline << std::endl;
  }

  inline void EPOS_Wrapper::print_legend(std::ostream& ostr) {
    char outline[81];
    sprintf(outline,
            "%4s %4s %4s %5s   %10s, %9s, %9s, %9s, %10s",
            "Indx",
            "Stat",
            "Par-",
            "chil-",
            "(  P_x",
            "P_y",
            "P_z",
            "Energy",
            "M ) ");
    ostr << outline << std::endl;
    sprintf(
        outline, "%9s %4s %4s    %10s, %9s, %9s, %9s) %9s", "ID ", "ents", "dren", "Prod (   X", "Y", "Z", "cT", "[mm]");
    ostr << outline << std::endl;
  }

}  // namespace EPOS

#endif  // HEPMC_EPOS_WRAPPER_H
//--------------------------------------------------------------------------