From 85b32f809835ed87d7e16f070c8e7f44608c49ea Mon Sep 17 00:00:00 2001 From: Antonio Vilela Pereira Date: Wed, 25 Apr 2007 12:50:31 +0000 Subject: [PATCH] --- yaml --- r: 23565 b: "refs/heads/CMSSW_7_0_X" c: fa1fe5ed3b28110aa934eb767f89bdff8a52615c h: "refs/heads/CMSSW_7_0_X" i: 23563: e9ca5640f5b79a9473f4e60ac5b39ed6efd23711 --- [refs] | 2 +- .../PomwigInterface/BuildFile | 7 +- .../PomwigInterface/interface/PomwigFilter.h | 43 + .../PomwigInterface/interface/PomwigSource.h | 23 +- .../PomwigInterface/src/HEPEVT_Wrapper.cc | 228 ----- .../PomwigInterface/src/HEPEVT_Wrapper.h | 542 ------------ .../PomwigInterface/src/HepEvt.F | 10 - .../PomwigInterface/src/HerwigWrapper6_4.h | 148 ---- .../PomwigInterface/src/IO_BaseClass.h | 142 ---- .../PomwigInterface/src/IO_HERWIG.cc | 779 ------------------ .../PomwigInterface/src/IO_HERWIG.h | 128 --- .../PomwigInterface/src/ParticleData.h | 216 ----- .../PomwigInterface/src/ParticleDataTable.h | 233 ------ .../PomwigInterface/src/PomwigFilter.cc | 35 + .../PomwigInterface/src/PomwigSource.cc | 135 ++- .../PomwigInterface/src/SealModule.cc | 3 +- .../PomwigInterface/src/h1init.f | 14 + .../PomwigInterface/src/herwig.h | 10 +- .../PomwigInterface/src/hwabeg.f | 2 +- .../PomwigInterface/test/PomwigAnalyzer.cc | 23 +- .../PomwigInterface/test/PomwigAnalyzer.h | 6 +- .../PomwigInterface/test/Z2muAnalyzer.cc | 15 +- .../PomwigInterface/test/Z2muAnalyzer.h | 6 +- 23 files changed, 240 insertions(+), 2510 deletions(-) create mode 100644 trunk/GeneratorInterface/PomwigInterface/interface/PomwigFilter.h delete mode 100644 trunk/GeneratorInterface/PomwigInterface/src/HEPEVT_Wrapper.cc delete mode 100644 trunk/GeneratorInterface/PomwigInterface/src/HEPEVT_Wrapper.h delete mode 100644 trunk/GeneratorInterface/PomwigInterface/src/HepEvt.F delete mode 100644 trunk/GeneratorInterface/PomwigInterface/src/HerwigWrapper6_4.h delete mode 100644 trunk/GeneratorInterface/PomwigInterface/src/IO_BaseClass.h delete mode 100644 trunk/GeneratorInterface/PomwigInterface/src/IO_HERWIG.cc delete mode 100644 trunk/GeneratorInterface/PomwigInterface/src/IO_HERWIG.h delete mode 100644 trunk/GeneratorInterface/PomwigInterface/src/ParticleData.h delete mode 100644 trunk/GeneratorInterface/PomwigInterface/src/ParticleDataTable.h create mode 100644 trunk/GeneratorInterface/PomwigInterface/src/PomwigFilter.cc create mode 100644 trunk/GeneratorInterface/PomwigInterface/src/h1init.f diff --git a/[refs] b/[refs] index 1348ddce8e926..0bdc48bc4287f 100644 --- a/[refs] +++ b/[refs] @@ -1,2 +1,2 @@ --- -"refs/heads/CMSSW_7_0_X": f257566ccdcb74d942875dd7ee406acdb506742b +"refs/heads/CMSSW_7_0_X": fa1fe5ed3b28110aa934eb767f89bdff8a52615c diff --git a/trunk/GeneratorInterface/PomwigInterface/BuildFile b/trunk/GeneratorInterface/PomwigInterface/BuildFile index 80307b86a8a94..adebef8ad08f7 100644 --- a/trunk/GeneratorInterface/PomwigInterface/BuildFile +++ b/trunk/GeneratorInterface/PomwigInterface/BuildFile @@ -1,9 +1,10 @@ - + + - + @@ -11,7 +12,9 @@ + + diff --git a/trunk/GeneratorInterface/PomwigInterface/interface/PomwigFilter.h b/trunk/GeneratorInterface/PomwigInterface/interface/PomwigFilter.h new file mode 100644 index 0000000000000..421a94a07b9b6 --- /dev/null +++ b/trunk/GeneratorInterface/PomwigInterface/interface/PomwigFilter.h @@ -0,0 +1,43 @@ +#ifndef POMWIG_SOURCE_EMPTY_FILTER_H +#define POMWIG_SOURCE_EMPTY_FILTER_H + +// +// Original Author: Fabian Stoeckli +// Created: Mon Mar 12 17:34:14 CET 2007 +// $Id: Herwig6Filter.h,v 1.1 2007/03/15 10:19:11 fabstoec Exp $ +// +// + +// Filter to remove empty events produced with MC@NLO/HERWIG + +// Modified for POMWIG + +// system include files +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDFilter.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +// +// class declaration +// + +class PomwigFilter : public edm::EDFilter { + public: + explicit PomwigFilter(const edm::ParameterSet&); + ~PomwigFilter(); + + private: + virtual void beginJob(const edm::EventSetup&) ; + virtual bool filter(edm::Event&, const edm::EventSetup&); + virtual void endJob() ; + +}; + +#endif diff --git a/trunk/GeneratorInterface/PomwigInterface/interface/PomwigSource.h b/trunk/GeneratorInterface/PomwigInterface/interface/PomwigSource.h index d05a961c99041..b7528329a3894 100644 --- a/trunk/GeneratorInterface/PomwigInterface/interface/PomwigSource.h +++ b/trunk/GeneratorInterface/PomwigInterface/interface/PomwigSource.h @@ -1,11 +1,10 @@ #ifndef PomwigSource_h #define PomwigSource_h -/** \class PomwigSource +/** \class PomwigSource * * Generates Pomwig HepMC events * - * Based on Herwig6Source ***************************************/ @@ -14,16 +13,15 @@ #include "FWCore/ServiceRegistry/interface/Service.h" #include #include -#include "CLHEP/HepMC/GenEvent.h" +#include "HepMC/GenEvent.h" namespace edm { class PomwigSource : public GeneratedInputSource { + public: - /// Constructor PomwigSource(const ParameterSet &, const InputSourceDescription &); - /// Destructor virtual ~PomwigSource(); @@ -31,22 +29,23 @@ namespace edm virtual bool produce(Event & e); void clear(); - - bool hwgive(const std::string& iParm ); + + bool hwgive(const std::string& iParm ); bool setRngSeeds(int); - + HepMC::GenEvent *evt; - - /// Verbosity flag int herwigVerbosity_; bool herwigHepMCVerbosity_; int herwigLhapdfVerbosity_; int maxEventsToPrint_; double comenergy; - int diffTopology; - std::string lhapdfSetPath_; + bool useJimmy_; + bool doMPInteraction_; bool printCards_; + int numTrials_; + + int diffTopology; }; } diff --git a/trunk/GeneratorInterface/PomwigInterface/src/HEPEVT_Wrapper.cc b/trunk/GeneratorInterface/PomwigInterface/src/HEPEVT_Wrapper.cc deleted file mode 100644 index 6514f5192e296..0000000000000 --- a/trunk/GeneratorInterface/PomwigInterface/src/HEPEVT_Wrapper.cc +++ /dev/null @@ -1,228 +0,0 @@ -////////////////////////////////////////////////////////////////////////// -// Matt.Dobbs@Cern.CH June 30, 2000 -// Generic Wrapper for the fortran HEPEVT common block -// -// The static data member's initializations must be separate from .h file. -// -////////////////////////////////////////////////////////////////////////// - -#include "HEPEVT_Wrapper.h" - -namespace HepMC { - - //////////////////////////////////////// - // static data member initializations // - //////////////////////////////////////// - - unsigned int HEPEVT_Wrapper::s_sizeof_int = 4; - - unsigned int HEPEVT_Wrapper::s_sizeof_real = sizeof(double); - - unsigned int HEPEVT_Wrapper::s_max_number_entries = 4000; - - /////////////////// - // Print Methods // - /////////////////// - - void HEPEVT_Wrapper::print_hepevt( std::ostream& ostr ) - { - // dumps the content of this HEPEVT event to ostr (Width is 80) - ostr << "________________________________________" - << "________________________________________" << std::endl; - ostr << "***** HEPEVT Common Event#: " - << event_number() - << ", " << number_entries() << " particles (max " - << max_number_entries() << ") *****"; - if ( is_double_precision() ) { - ostr << " Double Precision" << std::endl; - } else { - ostr << " Single Precision" << std::endl; - } - 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_hepevt_particle( i, ostr ); - } - ostr << "________________________________________" - << "________________________________________" << std::endl; - } - - void HEPEVT_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; - } - - void HEPEVT_Wrapper::print_hepevt_particle( int i, std::ostream& ostr ) - { - // dumps the content HEPEVT particle entry i (Width is 120) - // here i is the C array index (i.e. it starts at 0 ... whereas the - // fortran array index starts at 1) So if there's 100 particles, the - // last valid index is 100-1=99 - 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; - } - - - bool HEPEVT_Wrapper::check_hepevt_consistency( std::ostream& os ) - { - // This method inspects the HEPEVT common block and looks for - // inconsistencies in the mother/daughter pointers - bool isConsistent=true; - char header[81]; - sprintf( header, - "\n\n\t**** WARNINGInconsistent HEPEVT input, Event %10d ****" - , HEPEVT_Wrapper::event_number() ); - - for ( int i = 1; i <= HEPEVT_Wrapper::number_entries(); ++i ) { - // 1. check its mothers - int moth1 = HEPEVT_Wrapper::first_parent( i ); - int moth2 = HEPEVT_Wrapper::last_parent( i ); - if ( moth2 last parent " << std::endl; - HEPEVT_Wrapper::print_hepevt_particle( i, os ); - } - for ( int m = moth1; m<=moth2 && m!=0; ++m ) { - if ( m>HEPEVT_Wrapper::number_entries() || m < 0 ) { - if ( isConsistent ) { - os << header << std::endl; - isConsistent = false; - print_legend(os); - } - os << "Inconsistent entry " << i - << " mother points out of range " << std::endl; - HEPEVT_Wrapper::print_hepevt_particle( i, os ); - } - int mChild1 = HEPEVT_Wrapper::first_child(m); - int mChild2 = HEPEVT_Wrapper::last_child(m); - // we don't consider null pointers as inconsistent - if ( mChild1==0 && mChild2==0 ) continue; - if ( imChild2 ) { - if ( isConsistent ) { - os << header << std::endl; - isConsistent = false; - print_legend(os); - } - os << "Inconsistent mother-daughter relationship between " - << i << " & " << m - << " (try !trust_mother)" << std::endl; - HEPEVT_Wrapper::print_hepevt_particle( i, os ); - HEPEVT_Wrapper::print_hepevt_particle( m, os ); - } - } - // 2. check its daughters - int dau1 = HEPEVT_Wrapper::first_child( i ); - int dau2 = HEPEVT_Wrapper::last_child( i ); - if ( dau2 last child " << std::endl; - HEPEVT_Wrapper::print_hepevt_particle( i, os ); - } - for ( int d = dau1; d<=dau2 && d!=0; ++d ) { - if ( d>HEPEVT_Wrapper::number_entries() || d < 0 ) { - if ( isConsistent ) { - os << header << std::endl; - isConsistent = false; - print_legend(os); - } - os << "Inconsistent entry " << i - << " child points out of range " << std::endl; - HEPEVT_Wrapper::print_hepevt_particle( i, os ); - } - int d_moth1 = HEPEVT_Wrapper::first_parent(d); - int d_moth2 = HEPEVT_Wrapper::last_parent(d); - // we don't consider null pointers as inconsistent - if ( d_moth1==0 && d_moth2==0 ) continue; - if ( id_moth2 ) { - if ( isConsistent ) { - os << header << std::endl; - isConsistent = false; - print_legend(os); - } - os << "Inconsistent mother-daughter relationship between " - << i << " & " << d - << " (try trust_mothers)"<< std::endl; - HEPEVT_Wrapper::print_hepevt_particle( i, os ); - HEPEVT_Wrapper::print_hepevt_particle( d, os ); - } - } - } - if (!isConsistent) { - os << "Above lists all the inconsistencies in the HEPEVT common " - << "\n block which has been provided as input to HepMC. " - << "\n HepMC WILL have trouble interpreting the mother-daughter" - << "\n relationships ... but all other information " - << "\n (4-vectors etc) will be correctly transferred." - << "\n In order for HepMC to be able to interpret the mother/" - << "\n daughter hierachy, it MUST be given consistent input." - << "\n This is one of the design criteria of HepMC: " - << "\n consistency is enforced by the code."; - os << "\nThere is a switch in IO_HEPEVT, set-able using " - << "\n IO_HEPEVT::set_trust_mothers_before_daughters( bool )" - << "\n which you may want to try."; - os << "\nNote: if HEPEVT common block has been filled by pythia" - << "\n pyhepc, then the switch MSTP(128)=2 should be used in" - << "\n pythia, which instructs pythia not to put multiple " - << "\n copies of resonances in the event record.\n"; - os << "To obtain a file summarizing the inconsistency, you should:" - << "\n\t ofstream myFile(\"myInconsistentEvent.txt\"); " - << "\n\t HEPEVT_Wrapper::check_hepevt_consistency(myFile); " - << "\n\t HEPEVT_Wrapper::print_hepevt(myFile); " - << "\n[now write the event to HepMC using something like" - << "\n\t\t myIO_HEPEVT->write_event(myEvent); ]" - << "\n\t myEvent->print( myFile ); " - << " // print event as HepMC sees it" - << "\n ------------------------- Thank-you. \n\n" << std::endl; - } - return isConsistent; - } - - void HEPEVT_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 ); - } - } - -} // HepMC - diff --git a/trunk/GeneratorInterface/PomwigInterface/src/HEPEVT_Wrapper.h b/trunk/GeneratorInterface/PomwigInterface/src/HEPEVT_Wrapper.h deleted file mode 100644 index ac34a91f4c16b..0000000000000 --- a/trunk/GeneratorInterface/PomwigInterface/src/HEPEVT_Wrapper.h +++ /dev/null @@ -1,542 +0,0 @@ -//-------------------------------------------------------------------------- - -#ifndef HEPEVT_EntriesAllocation -#define HEPEVT_EntriesAllocation 10000 -#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 - - const unsigned int hepevt_bytes_allocation = - sizeof(long int) * ( 2 + 4 * 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 -#include // needed for formatted output using sprintf - -namespace HepMC { - - class HEPEVT_Wrapper { - public: - - static void print_hepevt( std::ostream& ostr = std::cout ); - static void print_hepevt_particle( int index, - std::ostream& ostr = std::cout ); - static bool is_double_precision(); // True if common block uses double - - static bool check_hepevt_consistency( std::ostream& ostr = std::cout ); - - 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 ); - 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 ); - static double px( int index ); // X momentum - static double py( int index ); - static double pz( int index ); - 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 ); - static double z( int index ); - static double t( int index ); // production time - - //////////////////// - // 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 ); - ////////////////////// - // HEPEVT 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; - - }; - - ////////////////////////////// - // 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 -//-------------------------------------------------------------------------- - diff --git a/trunk/GeneratorInterface/PomwigInterface/src/HepEvt.F b/trunk/GeneratorInterface/PomwigInterface/src/HepEvt.F deleted file mode 100644 index c23d20cfb6ac3..0000000000000 --- a/trunk/GeneratorInterface/PomwigInterface/src/HepEvt.F +++ /dev/null @@ -1,10 +0,0 @@ - - - subroutine ihepevt - - -#include "CLHEP/HepMC/include/stdhep.inc" -#include "CLHEP/HepMC/include/hepev4.inc" - - return - end diff --git a/trunk/GeneratorInterface/PomwigInterface/src/HerwigWrapper6_4.h b/trunk/GeneratorInterface/PomwigInterface/src/HerwigWrapper6_4.h deleted file mode 100644 index 6c3c119c8ee1b..0000000000000 --- a/trunk/GeneratorInterface/PomwigInterface/src/HerwigWrapper6_4.h +++ /dev/null @@ -1,148 +0,0 @@ -//-------------------------------------------------------------------------- -#ifndef HERWIG_WRAPPER_H -#define HERWIG_WRAPPER_H - -////////////////////////////////////////////////////////////////////////// -// Matt.Dobbs@Cern.CH, April 2002 -// Wrapper for FORTRAN version of Herwig -////////////////////////////////////////////////////////////////////////// -// - -#include - -//-------------------------------------------------------------------------- -// 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_ - -// 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 { - 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); -} -/* -// Subroutine inside H1QCD -#define qcd_1994 qcd_1994_ -extern "C" { - void qcd_1994(double&,double&,double*,int); -} -// HWABEG to initialize H1 pomeron structure function -#define hwabeg hwabeg_ -extern "C" { - void hwabeg(); -} -*/ - -#endif // HERWIG_WRAPPER_H -//-------------------------------------------------------------------------- diff --git a/trunk/GeneratorInterface/PomwigInterface/src/IO_BaseClass.h b/trunk/GeneratorInterface/PomwigInterface/src/IO_BaseClass.h deleted file mode 100644 index 95b2523b156a8..0000000000000 --- a/trunk/GeneratorInterface/PomwigInterface/src/IO_BaseClass.h +++ /dev/null @@ -1,142 +0,0 @@ -//-------------------------------------------------------------------------- -#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 -#include "ParticleDataTable.h" -#include "CLHEP/HepMC/GenEvent.h" - -namespace HepMC { - - class IO_BaseClass { - public: - virtual ~IO_BaseClass() {} - - virtual void write_event( const GenEvent* ) =0; - virtual bool fill_next_event( GenEvent* ) =0; - virtual void write_particle_data_table( const ParticleDataTable* ) =0; - virtual bool fill_particle_data_table( ParticleDataTable* ) =0; - virtual void print( std::ostream& ostr = std::cout ) const; - // - // the read_next_event() and read_particle_data_table() differ from - // the fill_***() methods in that they create a new event or pdt - // before calling the correspondingfill_*** method - // (they are not intended to be over-ridden) - GenEvent* read_next_event(); - ParticleDataTable* read_particle_data_table(); - // - // 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) - virtual GenEvent*& operator>>( GenEvent*& ); - virtual const GenEvent*& operator<<( const GenEvent*& ); - virtual GenEvent*& operator<<( GenEvent*& ); - virtual ParticleDataTable*& operator>>( ParticleDataTable*& ); - virtual const ParticleDataTable*& operator<<( const - ParticleDataTable*& ); - virtual ParticleDataTable*& operator<<( ParticleDataTable*& ); - }; - - ////////////// - // 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 ParticleDataTable* IO_BaseClass::read_particle_data_table() { - // creates a new particle data table and fills it by calling - // the sister method read_particle_data_table( ParticleDataTable* ) - // - // 1. create an empty pdt - ParticleDataTable* pdt = new ParticleDataTable(); - // 2. fill the pdt container - if the read is successful, return the - // pointer, otherwise return null and delete the evt - if ( fill_particle_data_table( pdt ) ) return pdt; - // next statement is only reached if read fails - delete pdt; - 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; - } - - inline ParticleDataTable*& IO_BaseClass::operator>>( - ParticleDataTable*& pdt ){ - pdt = read_particle_data_table(); - return pdt; - } - - inline const ParticleDataTable*& IO_BaseClass::operator<<( - const ParticleDataTable*& pdt ) { - write_particle_data_table( pdt ); - return pdt; - } - - inline ParticleDataTable*& IO_BaseClass::operator<<( - ParticleDataTable*& pdt ) { - write_particle_data_table( pdt ); - return pdt; - } - -} // HepMC - -#endif // HEPMC_IO_BASECLASS_H -//-------------------------------------------------------------------------- - - - diff --git a/trunk/GeneratorInterface/PomwigInterface/src/IO_HERWIG.cc b/trunk/GeneratorInterface/PomwigInterface/src/IO_HERWIG.cc deleted file mode 100644 index 4fd4339e3aeec..0000000000000 --- a/trunk/GeneratorInterface/PomwigInterface/src/IO_HERWIG.cc +++ /dev/null @@ -1,779 +0,0 @@ -////////////////////////////////////////////////////////////////////////// -// Matt.Dobbs@Cern.CH, October 2002 -// Herwig 6.400 IO class -////////////////////////////////////////////////////////////////////////// - -#include "IO_HERWIG.h" -#include "CLHEP/HepMC/GenEvent.h" -#include "CLHEP/Vector/LorentzVector.h" -#include // needed for formatted output using sprintf - - -#define FourVector HepLorentzVector - -namespace HepMC { - - IO_HERWIG::IO_HERWIG() : m_trust_mothers_before_daughters(false), - m_trust_both_mothers_and_daughters(true), - m_print_inconsistency_errors(true), - m_no_gaps_in_barcodes(true), - m_herwig_to_pdg_id(100,0) - { - // These arrays are copied from Lynn Garren's stdhep 5.01. - // see http://www-pat.fnal.gov/stdhep.html - // Translation from HERWIG particle ID's to PDG particle ID's. - m_herwig_to_pdg_id[1] =1; - m_herwig_to_pdg_id[2] =2; - m_herwig_to_pdg_id[3] =3; - m_herwig_to_pdg_id[4] =4; - m_herwig_to_pdg_id[5] =5; - m_herwig_to_pdg_id[6] =6; - m_herwig_to_pdg_id[7] =7; - m_herwig_to_pdg_id[8] =8; - - m_herwig_to_pdg_id[11] =11; - m_herwig_to_pdg_id[12] =12; - m_herwig_to_pdg_id[13] =13; - m_herwig_to_pdg_id[14] =14; - m_herwig_to_pdg_id[15] =15; - m_herwig_to_pdg_id[16] =16; - - m_herwig_to_pdg_id[21] =21; - m_herwig_to_pdg_id[22] =22; - m_herwig_to_pdg_id[23] =23; - m_herwig_to_pdg_id[24] =24; - m_herwig_to_pdg_id[25] =25; - m_herwig_to_pdg_id[26] =51; // <-- - - m_herwig_to_pdg_id[32] =32; - m_herwig_to_pdg_id[35] =35; - m_herwig_to_pdg_id[36] =36; - m_herwig_to_pdg_id[37] =37; - m_herwig_to_pdg_id[39] =39; - - m_herwig_to_pdg_id[81] =81; - m_herwig_to_pdg_id[82] =82; - m_herwig_to_pdg_id[83] =83; - m_herwig_to_pdg_id[84] =84; - m_herwig_to_pdg_id[85] =85; - m_herwig_to_pdg_id[86] =86; - m_herwig_to_pdg_id[87] =87; - m_herwig_to_pdg_id[88] =88; - m_herwig_to_pdg_id[89] =89; - m_herwig_to_pdg_id[90] =90; - - m_herwig_to_pdg_id[91] =91; - m_herwig_to_pdg_id[92] =92; - m_herwig_to_pdg_id[93] =93; - m_herwig_to_pdg_id[94] =94; - m_herwig_to_pdg_id[95] =95; - m_herwig_to_pdg_id[96] =96; - m_herwig_to_pdg_id[97] =97; - m_herwig_to_pdg_id[98] =9920022; // <-- - m_herwig_to_pdg_id[99] =9922212; // <-- - - // These particle ID's have no antiparticle, so aren't allowed. - m_no_antiparticles.insert(-21); - m_no_antiparticles.insert(-22); - m_no_antiparticles.insert(-23); - m_no_antiparticles.insert(-25); - m_no_antiparticles.insert(-51); - m_no_antiparticles.insert(-35); - m_no_antiparticles.insert(-36); - } - - IO_HERWIG::~IO_HERWIG(){} - - void IO_HERWIG::print( std::ostream& ostr ) const { - ostr << "IO_HERWIG: reads an event from the FORTRAN Herwig HEPEVT " - << "common block. \n" - << " trust_mothers_before_daughters = " - << m_trust_mothers_before_daughters - << " trust_both_mothers_and_daughters = " - << m_trust_both_mothers_and_daughters - << " print_inconsistency_errors = " - << m_print_inconsistency_errors << std::endl; - } - - bool IO_HERWIG::fill_next_event( GenEvent* evt ) { - // read one event from the Herwig HEPEVT common block and fill GenEvent - // return T/F =success/failure - // - // 0. Test that evt pointer is not null and set event number - if ( !evt ) { - std::cerr - << "IO_HERWIG::fill_next_event error - passed null event." - << std::endl; - return 0; - } - - // 1. First we have to fix the HEPEVT input, which is all mucked up for - // herwig. - repair_hepevt(); - - evt->set_event_number( HEPEVT_Wrapper::event_number() ); - // - // 2. create a particle instance for each HEPEVT entry and fill a map - // create a vector which maps from the HEPEVT particle index to the - // GenParticle address - // (+1 in size accounts for hepevt_particle[0] which is unfilled) - std::vector hepevt_particle( - HEPEVT_Wrapper::number_entries()+1 ); - hepevt_particle[0] = 0; - for ( int i1 = 1; i1 <= HEPEVT_Wrapper::number_entries(); ++i1 ) { - hepevt_particle[i1] = build_particle(i1); - } - std::set new_vertices; - // - // 3. We need to take special care with the hard process - // vertex. The problem we are trying to avoid is when the - // partons entering the hard process also have daughters from - // the parton shower. When this happens, each one can get its - // own decay vertex, making it difficult to join them - // later. We handle it by joining them together first, then - // the other daughters get added on later. - // Find the partons entering the hard vertex (status codes 121, 122). - int index_121 = 0; - int index_122 = 0; - for ( int i = 1; i <=HEPEVT_Wrapper::number_entries(); i++ ) { - if ( HEPEVT_Wrapper::status(i)==121 ) index_121=i; - if ( HEPEVT_Wrapper::status(i)==122 ) index_122=i; - if ( index_121!=0 && index_122!=0 ) break; - } - if ( index_121 && index_122 ) { - GenVertex* hard_vtx = new GenVertex(); - hard_vtx->add_particle_in( hepevt_particle[index_121] ); - hard_vtx->add_particle_in( hepevt_particle[index_122] ); - // evt->add_vertex( hard_vtx ); // not necessary, its done in - // set_signal_process_vertex - evt->set_signal_process_vertex( hard_vtx ); - } - // - // 4. loop over HEPEVT particles AGAIN, this time creating vertices - for ( int i = 1; i <= HEPEVT_Wrapper::number_entries(); ++i ) { - // We go through and build EITHER the production or decay - // vertex for each entry in hepevt, depending on the switch - // m_trust_mothers_before_daughters (new 2001-02-28) - // Note: since the HEPEVT pointers are bi-directional, it is - /// sufficient to do one or the other. - // - // 3. Build the production_vertex (if necessary) - if ( m_trust_mothers_before_daughters || - m_trust_both_mothers_and_daughters ) { - build_production_vertex( i, hepevt_particle, evt ); - } - // - // 4. Build the end_vertex (if necessary) - // Identical steps as for production vertex - if ( !m_trust_mothers_before_daughters || - m_trust_both_mothers_and_daughters ) { - build_end_vertex( i, hepevt_particle, evt ); - } - } - // 5. 01.02.2000 - // handle the case of particles in HEPEVT which come from nowhere - - // i.e. particles without mothers or daughters. - // These particles need to be attached to a vertex, or else they - // will never become part of the event. check for this situation. - for ( int i3 = 1; i3 <= HEPEVT_Wrapper::number_entries(); ++i3 ) { - // Herwig also has some non-physical entries in HEPEVT - // like CMS, HARD, and CONE. These are flagged by - // repair_hepevt by making their status and id zero. We - // delete those particles here. - if ( hepevt_particle[i3] && !hepevt_particle[i3]->parent_event() - && !hepevt_particle[i3]->pdg_id() - && !hepevt_particle[i3]->status() ) { - //std::cout << "IO_HERWIG::fill_next_event is deleting null " - // << "particle" << std::endl; - //hepevt_particle[i3]->print(); - delete hepevt_particle[i3]; - } else if ( hepevt_particle[i3] && - !hepevt_particle[i3]->end_vertex() && - !hepevt_particle[i3]->production_vertex() ) { - GenVertex* prod_vtx = new GenVertex(); - prod_vtx->add_particle_out( hepevt_particle[i3] ); - evt->add_vertex( prod_vtx ); - } - } - return true; - } - - void IO_HERWIG::build_production_vertex(int i, - std::vector& - hepevt_particle, - GenEvent* evt ) { - // - // for particle in HEPEVT with index i, build a production vertex - // if appropriate, and add that vertex to the event - GenParticle* p = hepevt_particle[i]; - // a. search to see if a production vertex already exists - int mother = HEPEVT_Wrapper::first_parent(i); - GenVertex* prod_vtx = p->production_vertex(); - while ( !prod_vtx && mother > 0 ) { - prod_vtx = hepevt_particle[mother]->end_vertex(); - if ( prod_vtx ) prod_vtx->add_particle_out( p ); - // increment mother for next iteration - if ( ++mother > HEPEVT_Wrapper::last_parent(i) ) mother = 0; - } - // b. if no suitable production vertex exists - and the particle - // has atleast one mother or position information to store - - // make one - FourVector prod_pos( HEPEVT_Wrapper::x(i), HEPEVT_Wrapper::y(i), - HEPEVT_Wrapper::z(i), HEPEVT_Wrapper::t(i) - ); - if ( !prod_vtx && (HEPEVT_Wrapper::number_parents(i)>0 - || prod_pos!=FourVector(0,0,0,0)) ) - { - prod_vtx = new GenVertex(); - prod_vtx->add_particle_out( p ); - evt->add_vertex( prod_vtx ); - } - // c. if prod_vtx doesn't already have position specified, fill it - if ( prod_vtx && prod_vtx->position()==FourVector(0,0,0,0) ) { - prod_vtx->set_position( prod_pos ); - } - // d. loop over mothers to make sure their end_vertices are - // consistent - mother = HEPEVT_Wrapper::first_parent(i); - while ( prod_vtx && mother > 0 ) { - if ( !hepevt_particle[mother]->end_vertex() ) { - // if end vertex of the mother isn't specified, do it now - prod_vtx->add_particle_in( hepevt_particle[mother] ); - } else if (hepevt_particle[mother]->end_vertex() != prod_vtx ) { - // problem scenario --- the mother already has a decay - // vertex which differs from the daughter's produciton - // vertex. This means there is internal - // inconsistency in the HEPEVT event record. Print an - // error - // Note: we could provide a fix by joining the two - // vertices with a dummy particle if the problem - // arrises often with any particular generator. - if ( m_print_inconsistency_errors ) { - std::cerr - << "HepMC::IO_HERWIG: inconsistent mother/daugher " - << "information in HEPEVT event " - << HEPEVT_Wrapper::event_number() - << ". \n I recommend you try " - << "inspecting the event first with " - << "\n\tHEPEVT_Wrapper::check_hepevt_consistency()" - << "\n This warning can be turned off with the " - << "IO_HERWIG::print_inconsistency_errors switch." - << std::endl; - hepevt_particle[mother]->print(std::cerr); - std::cerr - << "problem vertices are: (prod_vtx, mother)" << std::endl; - if ( prod_vtx ) prod_vtx->print(std::cerr); - hepevt_particle[mother]->end_vertex()->print(std::cerr); - } - } - if ( ++mother > HEPEVT_Wrapper::last_parent(i) ) mother = 0; - } - } - - void IO_HERWIG::build_end_vertex - ( int i, std::vector& hepevt_particle, GenEvent* evt ) - { - // - // for particle in HEPEVT with index i, build an end vertex - // if appropriate, and add that vertex to the event - // Identical steps as for build_production_vertex - GenParticle* p = hepevt_particle[i]; - // a. - int daughter = HEPEVT_Wrapper::first_child(i); - GenVertex* end_vtx = p->end_vertex(); - while ( !end_vtx && daughter > 0 ) { - end_vtx = hepevt_particle[daughter]->production_vertex(); - if ( end_vtx ) end_vtx->add_particle_in( p ); - if ( ++daughter > HEPEVT_Wrapper::last_child(i) ) daughter = 0; - } - // b. (different from 3c. because HEPEVT particle can not know its - // decay position ) - if ( !end_vtx && HEPEVT_Wrapper::number_children(i)>0 ) { - end_vtx = new GenVertex(); - end_vtx->add_particle_in( p ); - evt->add_vertex( end_vtx ); - } - // c+d. loop over daughters to make sure their production vertices - // point back to the current vertex. - // We get the vertex position from the daughter as well. - daughter = HEPEVT_Wrapper::first_child(i); - while ( end_vtx && daughter > 0 ) { - if ( !hepevt_particle[daughter]->production_vertex() ) { - // if end vertex of the mother isn't specified, do it now - end_vtx->add_particle_out( hepevt_particle[daughter] ); - // - // 2001-03-29 M.Dobbs, fill vertex the position. - if ( end_vtx->position()==FourVector(0,0,0,0) ) { - FourVector prod_pos( HEPEVT_Wrapper::x(daughter), - HEPEVT_Wrapper::y(daughter), - HEPEVT_Wrapper::z(daughter), - HEPEVT_Wrapper::t(daughter) - ); - if ( prod_pos != FourVector(0,0,0,0) ) { - end_vtx->set_position( prod_pos ); - } - } - } else if (hepevt_particle[daughter]->production_vertex() - != end_vtx){ - // problem scenario --- the daughter already has a prod - // vertex which differs from the mother's end - // vertex. This means there is internal - // inconsistency in the HEPEVT event record. Print an - // error - if ( m_print_inconsistency_errors ) std::cerr - << "HepMC::IO_HERWIG: inconsistent mother/daugher " - << "information in HEPEVT event " - << HEPEVT_Wrapper::event_number() - << ". \n I recommend you try " - << "inspecting the event first with " - << "\n\tHEPEVT_Wrapper::check_hepevt_consistency()" - << "\n This warning can be turned off with the " - << "IO_HERWIG::print_inconsistency_errors switch." - << std::endl; - } - if ( ++daughter > HEPEVT_Wrapper::last_child(i) ) daughter = 0; - } - if ( !p->end_vertex() && !p->production_vertex() ) { - // Added 2001-11-04, to try and handle Isajet problems. - build_production_vertex( i, hepevt_particle, evt ); - } - } - - GenParticle* IO_HERWIG::build_particle( int index ) { - // Builds a particle object corresponding to index in HEPEVT - // - GenParticle* p - = new GenParticle( FourVector( HEPEVT_Wrapper::px(index), - HEPEVT_Wrapper::py(index), - HEPEVT_Wrapper::pz(index), - HEPEVT_Wrapper::e(index) ), - HEPEVT_Wrapper::id(index), - HEPEVT_Wrapper::status(index) ); - p->setGeneratedMass( HEPEVT_Wrapper::m(index) ); - p->suggest_barcode( index ); - return p; - } - - int IO_HERWIG::find_in_map( const std::map& m, - GenParticle* p) const { - std::map::const_iterator iter = m.find(p); - if ( iter == m.end() ) return 0; - return iter->second; - } - - void IO_HERWIG::repair_hepevt() const { - // This routine takes the HEPEVT common block as used in HERWIG, - // and converts it into the HEPEVT common block in the standard format - // - // This means it: - // - removes the color structure, which herwig overloads - // into the mother/daughter fields - // - zeros extra entries for hard subprocess, etc. - // - // - // Special HERWIG status codes - // 101,102 colliding beam particles - // 103 beam-beam collision CMS vector - // 120 hard subprocess CMS vector - // 121,122 hard subprocess colliding partons - // 123-129 hard subprocess outgoing particles - // 141-149 (ID=94) mirror image of hard subrpocess particles - // 100 (ID=0 cone) - // - // Special HERWIG particle id's - // 91 clusters - // 94 jets - // 0 others with no pdg code - - // Make sure hepvt isn't empty. - if ( HEPEVT_Wrapper::number_entries() <= 0 ) return; - - // Find the index of the beam-beam collision and of the hard subprocess - // Later we will assume that - // 101 ---> 121 \. - // X Hard subprocess - // 102 ---> 122 / - // - int index_collision = 0; - int index_hard = 0; - int index_101 = 0; - int index_102 = 0; - int index_121 = 0; - int index_122 = 0; - - for ( int i = 1; i <=HEPEVT_Wrapper::number_entries(); i++ ) { - if ( HEPEVT_Wrapper::status(i)==101 ) index_101=i; - if ( HEPEVT_Wrapper::status(i)==102 ) index_102=i; - if ( HEPEVT_Wrapper::status(i)==103 ) index_collision=i; - if ( HEPEVT_Wrapper::status(i)==120 ) index_hard=i; - if ( HEPEVT_Wrapper::status(i)==121 ) index_121=i; - if ( HEPEVT_Wrapper::status(i)==122 ) index_122=i; - if ( index_collision!=0 && index_hard!=0 && index_101!=0 && - index_102!=0 && index_121!=0 && index_122!=0 ) break; - } - - // The mother daughter information for the hard subprocess entry (120) - // IS correct, whereas the information for the particles participating - // in the hard subprocess contains instead the color flow relationships - // Transfer the hard subprocess info onto the other particles - // in the hard subprocess. - // - // We cannot specify daughters of the incoming hard process particles - // because they have some daughters (their showered versions) which - // are not adjacent in the particle record, so we cannot properly - // set the daughter indices in hepevt. - // - if (index_121) HEPEVT_Wrapper::set_parents(index_121, index_101, 0 ); - if (index_121) HEPEVT_Wrapper::set_children( index_121, 0, 0 ); - if (index_122) HEPEVT_Wrapper::set_parents(index_122, index_102, 0 ); - if (index_122) HEPEVT_Wrapper::set_children( index_122, 0, 0 ); - - for ( int i = HEPEVT_Wrapper::first_child(index_hard); - i <= HEPEVT_Wrapper::last_child(index_hard); i++ ) { - HEPEVT_Wrapper::set_parents( - i, HEPEVT_Wrapper::first_parent(index_hard), - HEPEVT_Wrapper::last_parent(index_hard) ); - - // When the direct descendants of the hard process are hadrons, - // then the 2nd child contains color flow information, and so - // we zero it. - // However, if the direct descendant is status=195, then it is - // a non-hadron, and so the 2nd child does contain real mother - // daughter relationships. ( particularly relevant for H->WW, - // April 18, 2003 ) - if ( HEPEVT_Wrapper::status(i) != 195 ) { - HEPEVT_Wrapper::set_children(i,HEPEVT_Wrapper::first_child(i),0); - } - } - - // now zero the collision and hard entries. - if (index_hard) zero_hepevt_entry(index_hard); - if (index_hard) zero_hepevt_entry(index_collision); - - // Loop over the particles individually and handle oddities - for ( int i=1; i <=HEPEVT_Wrapper::number_entries(); i++ ) { - - // ----------- Fix ID codes ---------- - // particles with ID=94 are mirror images of their mothers: - if ( HEPEVT_Wrapper::id(i)==94 ) { - HEPEVT_Wrapper::set_id( - i, HEPEVT_Wrapper::id( HEPEVT_Wrapper::first_parent(i) ) ); - } - - // ----------- fix STATUS codes ------ - // status=100 particles are "cones" which carry only color info - // throw them away - if ( HEPEVT_Wrapper::status(i)==100 ) zero_hepevt_entry(i); - - - // NOTE: status 101,102 particles are the beam particles. - // status 121,129 particles are the hard subprocess particles - // we choose to allow the herwig particles to have herwig - // specific codes, and so we don't bother to change these - // to status =3. - - - - - // ----------- fix some MOTHER/DAUGHTER relationships - // Whenever the mother points to the hard process, it is referring - // to a color flow, so we zero it. - if ( HEPEVT_Wrapper::last_parent(i)==index_hard ) { - HEPEVT_Wrapper::set_parents( - i, HEPEVT_Wrapper::first_parent(i), 0 ); - } - - // It makes no sense to have a mother that is younger than you are! - - if ( HEPEVT_Wrapper::first_parent(i) >= i ) { - HEPEVT_Wrapper::set_parents( i, 0, 0 ); - } - if ( HEPEVT_Wrapper::last_parent(i) >= i ) { - HEPEVT_Wrapper::set_parents( - i, HEPEVT_Wrapper::first_parent(i), 0 ); - } - - // Whenever the second mother/daughter has a lower index than the - // first, it means the second mother/daughter contains color - // info. Purge it. - if ( HEPEVT_Wrapper::last_parent(i) <= - HEPEVT_Wrapper::first_parent(i) ) { - HEPEVT_Wrapper::set_parents( - i, HEPEVT_Wrapper::first_parent(i), 0 ); - } - - if ( HEPEVT_Wrapper::last_child(i) <= - HEPEVT_Wrapper::first_child(i) ) { - HEPEVT_Wrapper::set_children( - i, HEPEVT_Wrapper::first_child(i), 0 ); - } - - // The mothers & daughters of a soft centre of mass (stat=170) seem - // to be correct, but they are out of sequence. The information is - // elsewhere in the event record, so zero it. - // - if ( HEPEVT_Wrapper::status(i) == 170 ) { - HEPEVT_Wrapper::set_parents( i, 0, 0 ); - HEPEVT_Wrapper::set_children( i, 0, 0 ); - } - - // Recognise clusters. - // Case 1: cluster has particle parents. - // Clusters normally DO point to its two - // correct mothers, but those 2 mothers are rarely adjacent in the - // event record ... so the mother information might say something - // like 123,48 where index123 and index48 really are the correct - // mothers... however the hepevt standard states that the mother - // pointers should give the index range. So we would have to - // reorder the event record and add entries if we wanted to use - // it. Instead we just zero the mothers, since all of that - // information is contained in the daughter information of the - // mothers. - // Case 2: cluster has a soft process centre of mass (stat=170) - // as parent. This is ok, keep it. - // - // Note if we were going directly to HepMC, then we could - // use this information properly! - - if ( HEPEVT_Wrapper::id(i)==91 ) { - // if the cluster comes from a SOFT (id=0,stat=170) - if ( HEPEVT_Wrapper::status(HEPEVT_Wrapper::first_parent(i)) - == 170 ) { - ; // In this case the mothers are ok - } else { - HEPEVT_Wrapper::set_parents( i, 0, 0 ); - } - } - } - - // ---------- Loop over the particles individually and look - // for mother/daughter inconsistencies. - // We consider a mother daughter relationship to be valid - // ONLy when the mother points to the daughter AND the - // daughter points back (true valid bidirectional - // pointers) OR when a one thing points to the other, but - // the other points to zero. If this isn't true, we zero - // the offending relationship. - - for ( int i=1; i <=HEPEVT_Wrapper::number_entries(); i++ ) { - // loop over parents - int ifirst = HEPEVT_Wrapper::first_parent(i); - int ilast = HEPEVT_Wrapper::last_parent(i); - if ( ilast == 0 ) ilast = HEPEVT_Wrapper::first_parent(i); - bool first_is_acceptable = true; - bool last_is_acceptable = true; - // check for out of range. - if ( ifirst>=i || ifirst<0 ) first_is_acceptable = false; - if ( ilast>=i || ilast=i ) {;} - else if ( HEPEVT_Wrapper::first_child(j) ==0 && - HEPEVT_Wrapper::last_child(j) ==0 ) {;} - - // Error Condition: - // modified by MADobbs@lbl.gov April 21, 2003 - // we distinguish between the first parent and all parents - // being incorrect - else if (j==ifirst) { first_is_acceptable = false; break; } - else { last_is_acceptable = false; break; } - } - } - // if any one of the mothers gave a bad outcome, zero all mothers - if ( !first_is_acceptable ) { - HEPEVT_Wrapper::set_parents( i, 0, 0 ); - } else if ( !last_is_acceptable ) { - HEPEVT_Wrapper::set_parents(i,HEPEVT_Wrapper::first_parent(i),0); - } - } - // Note: it's important to finish the mother loop, before - // starting the daughter loop ... since many mother relations - // will be zero'd which will validate the daughters.... i.e., - // we want relationships like: - // IHEP ID IDPDG IST MO1 MO2 DA1 DA2 - // 27 TQRK 6 3 26 26 30 30 - // 30 TQRK 6 155 26 11 31 32 - // to come out right. - - for ( int i=1; i <=HEPEVT_Wrapper::number_entries(); i++ ) { - // loop over daughters - int ifirst = HEPEVT_Wrapper::first_child(i); - int ilast = HEPEVT_Wrapper::last_child(i); - if ( ilast==0 ) ilast = HEPEVT_Wrapper::first_child(i); - bool is_acceptable = true; - // check for out of range. - if ( ifirst<=i || ifirst<0 ) is_acceptable = false; - if ( ilast<=i || ilast=i ) {;} - else if ( HEPEVT_Wrapper::first_parent(j) ==0 && - HEPEVT_Wrapper::last_parent(j) ==0 ) {;} - else { is_acceptable = false; } // error condition - } - } - // if any one of the children gave a bad outcome, zero all children - if ( !is_acceptable ) HEPEVT_Wrapper::set_children( i, 0, 0 ); - } - - // fixme - - for ( int i=1; i <=HEPEVT_Wrapper::number_entries(); i++ ) { - HEPEVT_Wrapper::set_id( - i, translate_herwig_to_pdg_id(HEPEVT_Wrapper::id(i)) ); - } - - - if ( m_no_gaps_in_barcodes ) remove_gaps_in_hepevt(); - } - - void IO_HERWIG::remove_gaps_in_hepevt() const { - // in this scenario, we do not allow there to be zero-ed - // entries in the HEPEVT common block, and so be reshuffle - // the common block, removing the zeero-ed entries as we - // go and making sure we keep the mother/daughter - // relationships appropriate - - std::vector mymap(HEPEVT_Wrapper::number_entries()+1,0); - int ilast = 0; - for ( int i=1; i <=HEPEVT_Wrapper::number_entries(); i++ ) { - if (HEPEVT_Wrapper::status(i)==0 && HEPEVT_Wrapper::id(i)==0) { - // we remove all entries for which stat=0, id=0 - mymap[i]=0; - } else { - ilast += 1; - if ( ilast != i ) { - HEPEVT_Wrapper::set_status(ilast, - HEPEVT_Wrapper::status(i) ); - HEPEVT_Wrapper::set_id(ilast, HEPEVT_Wrapper::id(i) ); - HEPEVT_Wrapper::set_parents( - ilast, - HEPEVT_Wrapper::first_parent(i), - HEPEVT_Wrapper::last_parent(i) ); - HEPEVT_Wrapper::set_children( - ilast, - HEPEVT_Wrapper::first_child(i), - HEPEVT_Wrapper::last_child(i) ); - HEPEVT_Wrapper::set_momentum( - ilast, - HEPEVT_Wrapper::px(i), HEPEVT_Wrapper::py(i), - HEPEVT_Wrapper::pz(i), HEPEVT_Wrapper::e(i) ); - HEPEVT_Wrapper::set_mass(ilast, HEPEVT_Wrapper::m(i) ); - HEPEVT_Wrapper::set_position( - ilast, HEPEVT_Wrapper::x(i),HEPEVT_Wrapper::y(i), - HEPEVT_Wrapper::z(i),HEPEVT_Wrapper::t(i) ); - } - mymap[i]=ilast; - } - } - - // M. Dobbs (from Borut) - April 26, to fix tauolo/herwig past - // the end problem with daughter pointers: - // HEPEVT_Wrapper::set_number_entries( ilast ); - - // Finally we need to re-map the mother/daughter pointers. - for ( int i=1; i <=ilast; i++ ) { - - HEPEVT_Wrapper::set_parents( - i, - mymap[HEPEVT_Wrapper::first_parent(i)], - mymap[HEPEVT_Wrapper::last_parent(i)] ); - HEPEVT_Wrapper::set_children( - i, - mymap[HEPEVT_Wrapper::first_child(i)], - mymap[HEPEVT_Wrapper::last_child(i)] ); - } - // M. Dobbs (from Borut, part B) - April 26, to fix tauolo/herwig past - // the end problem with daughter pointers: - HEPEVT_Wrapper::set_number_entries( ilast ); - } - - void IO_HERWIG::zero_hepevt_entry( int i ) const { - if ( i <=0 || i > HepMC::HEPEVT_Wrapper::max_number_entries() ) return; - HEPEVT_Wrapper::set_status( i, 0 ); - HEPEVT_Wrapper::set_id( i, 0 ); - HEPEVT_Wrapper::set_parents( i, 0, 0 ); - HEPEVT_Wrapper::set_children( i, 0, 0 ); - HEPEVT_Wrapper::set_momentum( i, 0, 0, 0, 0 ); - HEPEVT_Wrapper::set_mass( i, 0 ); - HEPEVT_Wrapper::set_position( i, 0, 0, 0, 0 ); - } - - int IO_HERWIG::translate_herwig_to_pdg_id( int id ) const { - // This routine is copied from Lynn Garren's stdhep 5.01. - // see http://www-pat.fnal.gov/stdhep.html - - // example -9922212 - int hwtran = id; // -9922212 - int ida = abs(id); // 9922212 - int j1 = ida%10; // 2 - int i1 = (ida/10)%10; // 1 - int i2 = (ida/100)%10; // 2 - int i3 = (ida/1000)%10; // 2 - //int i4 =(ida/10000)%10; // 2 - //int i5 =(ida/100000)%10; // 9 - //int k99 = (ida/100000)%100; // 9 - int ksusy = (ida/1000000)%10; // 0 - //int ku = (ida/10000000)%10; // 0 - int kqn = (ida/1000000000)%10; // 0 - - if ( kqn==1 ) { - // ions not recognized - hwtran=0; - if ( m_print_inconsistency_errors ) { - std::cerr << "IO_HERWIG::translate_herwig_to_pdg_id " << id - << "nonallowed ion" << std::endl; - } - } - else if (ida < 100) { - // Higgs, etc. - hwtran = m_herwig_to_pdg_id[ida]; - if ( id < 0 ) hwtran *= -1; - // check for illegal antiparticles - if ( id < 0 ) { - if ( hwtran>=-99 && hwtran<=-81) hwtran=0; - if ( m_no_antiparticles.count(hwtran) ) hwtran=0; - } - } - else if ( ksusy==1 || ksusy==2 ) { ; } - // SUSY - else if ( i1!=0 && i3!=0 && j1==2 ) {;} - // spin 1/2 baryons - else if ( i1!=0 && i3!=0 && j1==4 ) {;} - // spin 3/2 baryons - else if ( i1!=0 && i2!=0 && i3==0 ) { - // mesons - // check for illegal antiparticles - if ( i1==i2 && id<0) hwtran=0; - } - else if ( i2!=0 && i3!=0 && i1==0 ) {;} - // diquarks - else { - // undefined - hwtran=0; - } - - // check for illegal anti KS, KL - if ( id==-130 || id==-310 ) hwtran=0; - - if ( hwtran==0 && ida!=0 && m_print_inconsistency_errors ) { - std::cerr - << "IO_HERWIG::translate_herwig_to_pdg_id HERWIG particle " - << id << " translates to zero." << std::endl; - } - - return hwtran; - } - -} // HepMC - - - - diff --git a/trunk/GeneratorInterface/PomwigInterface/src/IO_HERWIG.h b/trunk/GeneratorInterface/PomwigInterface/src/IO_HERWIG.h deleted file mode 100644 index 61667cf5adc61..0000000000000 --- a/trunk/GeneratorInterface/PomwigInterface/src/IO_HERWIG.h +++ /dev/null @@ -1,128 +0,0 @@ -//-------------------------------------------------------------------------- -#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 -#include -#include "IO_BaseClass.h" -#include "HEPEVT_Wrapper.h" - -namespace HepMC { - - class GenEvent; - class GenVertex; - class GenParticle; - class ParticleDataTable; - - class IO_HERWIG : public IO_BaseClass { - public: - IO_HERWIG(); - virtual ~IO_HERWIG(); - bool fill_next_event( GenEvent* ); - void print( std::ostream& ostr = std::cout ) const; - double interfaces_to_version_number() const {return 6.400;} - - // see comments below for these switches. - bool print_inconsistency_errors() const; - void set_print_inconsistency_errors( bool b = 1 ); - - bool no_gaps_in_barcodes() const - { return m_no_gaps_in_barcodes; } - void set_no_gaps_in_barcodes( bool a ) - { m_no_gaps_in_barcodes=a; } - - protected: // for internal use only - bool trust_both_mothers_and_daughters() const; - bool trust_mothers_before_daughters() const; - void set_trust_mothers_before_daughters( bool b = 1 ); - void set_trust_both_mothers_and_daughters( bool b = 0 ); - - GenParticle* build_particle( int index ); - void build_production_vertex( - int i,std::vector& hepevt_particle, GenEvent* evt ); - void build_end_vertex( - int i, std::vector& hepevt_particle, GenEvent* evt ); - int find_in_map( - const std::map& m, GenParticle* p) const; - - void repair_hepevt() const; - void remove_gaps_in_hepevt() const; - void zero_hepevt_entry( int i ) const; - int translate_herwig_to_pdg_id( int i ) const; - - private: // following are not implemented for Herwig - virtual void write_event( const GenEvent* ){} - virtual void write_particle_data_table( const ParticleDataTable* ){} - virtual bool fill_particle_data_table( ParticleDataTable* ) - { return 0; } - - 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 m_herwig_to_pdg_id; - std::set 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 -//-------------------------------------------------------------------------- diff --git a/trunk/GeneratorInterface/PomwigInterface/src/ParticleData.h b/trunk/GeneratorInterface/PomwigInterface/src/ParticleData.h deleted file mode 100644 index cb01c46b2ed9f..0000000000000 --- a/trunk/GeneratorInterface/PomwigInterface/src/ParticleData.h +++ /dev/null @@ -1,216 +0,0 @@ -//-------------------------------------------------------------------------- -#ifndef HEPMC_PARTICLE_DATA_H -#define HEPMC_PARTICLE_DATA_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). -// -// GenParticle Data common to all particles of a given PDG id -////////////////////////////////////////////////////////////////////////// -// -// Units ID: defined by PDG group (particles are +ve, antiparticles are -ve) -// also consistent with the Pythia definitions -// See: http://d0lblt.lbl.gov/wwwpdg/mc_numbers.htm -// charge: fraction of proton charge -// mass in user defined energy units -// width: ( stored as cLifetime = hbar / Width ) -// cLifetime: c*time -// spin: fraction of photon spin (always a positive number) -// -// Default mass is 0. -// Default cLifetime is -1 which means stable (setting width = 0 gives this) -// (we define cLifetime = -1 --> width = 0 (i.e stable), -// width = -1 --> cLifetime = 0 (i.e. prompt) ) -// These defaults exist because many very basic MC generators -// may produce only massless stable particles in the event record. -// -// It is intended that a different ParticleData object is created for each -// particle and its anti-particle - useful for CP violation studies. -// -// There are few set methods for this class, there should be no reason -// to change anything after instantiating. If you need to, then -// create a new object and kill the old one. -// -// Example: -// HepMC::ParticleData* pd_electron = -// new HepMC::ParticleData("electron",11,-1,0.000511,-1,.5); -// A method is provided to allow you to set the lifetime from the -// width in the constructor -// Example: new HepMC::ParticleData("W+",24,+1,80.396, -// HepMC::clifetime_from_width(2.06),1); -// Example of finding a ParticleData object from its PDG ID -// in ParticleDataTable pdt: -// HepMC::ParticleData* electron = pdt.find(11); -// or if you just wanted two know the electron mass, you could do: -// pdt.find(11)->mass(); - -#include -#include -#include - -namespace HepMC { - - // hbar * c --> calculated with units of [mm*GeV] - static const double HepMC_hbarc = (6.6260755e-34 * (1.e-6/1.60217733e-19) / (2*3.14159265358979323846)) - * (2.99792458e+8 * 1000.) * 1.e+3; - - // if you want to instantiate the particle lifetime from its width, - // use this static method inside the constructor: - double clifetime_from_width( double width ); - - class ParticleData { - - friend std::ostream& operator<<( std::ostream&, const ParticleData& ); - - public: - ParticleData( std::string name, int id, double charge, double mass = 0, - double cLifetime = -1, double spin = 0 ); - ParticleData( const char* name, int id, double charge, double mass = 0, - double cLifetime = -1, double spin = 0 ); - virtual ~ParticleData(); - - bool operator==( const ParticleData& ) const; - bool operator!=( const ParticleData& ) const; - - void print( std::ostream& ostr = std::cout ) const; - - bool is_lepton() const; // true if charged lepton /neutrino - bool is_charged_lepton() const;// true if a charged lepton - bool is_em() const; // true if an electron or photon - bool is_neutrino() const; // true if a neutrino - bool is_hadron() const; // true if a hadron - bool is_boson() const; // true if a gauge or higgs boson - - //////////////////// - // access methods // - //////////////////// - std::string name() const; - int pdg_id() const; - double charge() const; - double mass() const; - double width() const; // width as calculated from clifetime - double clifetime() const; - double spin() const; - - void set_charge( double ); - void set_mass( double ); - void set_width( double ); - void set_clifetime( double ); - void set_spin( double ); - - protected: - static unsigned int counter(); // num ParticleData objects in memory - - // omits susy/excited/technicolor digit from returned ID - int model_independent_pdg_id_() const; - - private: - std::string m_name; // description of the particle according to PDG - // i.e. "Delta(1900) S_31" - int m_pdg_id; // PDG ID number (note we allow -ve) - int m_3charge;// 3*electric charge in units of proton charge - double m_mass; // nominal mass in user defined energy units - double m_clifetime; // [mm] - unsigned char m_2spin; // 2*spin (J) of particle - - static unsigned int s_counter; - }; - - /////////////////////////// - // INLINES // - /////////////////////////// - - inline bool ParticleData::is_lepton() const { - // true if a charged lepton or neutrino --> | 11,13,15,12,14,16,17,18 | - return ( abs(pdg_id()) >=11 && abs(pdg_id()) <= 18 ); - } - inline bool ParticleData::is_charged_lepton() const { - // true if a charged lepton --> | 11,13,15 | - return ( is_lepton() && abs(pdg_id())%2==1 ); - } - inline bool ParticleData::is_neutrino() const { - // true if a neutrino --> | 12,14,16 | - return ( is_lepton() && abs(pdg_id())%2==0 ); - } - inline bool ParticleData::is_em() const { - // true if an electron or photon --> | 11, 22 | - return ( abs(pdg_id()) == 11 || abs(pdg_id()) == 22 ); - } - inline bool ParticleData::is_hadron() const { - // true if a hadron --> q,g,meson,baryon - return ( abs(pdg_id()) <= 9 || abs(pdg_id()) == 21 - || abs(pdg_id()) >100 ); - } - inline bool ParticleData::is_boson() const { - // true if a gauge or higgs boson --> | 9, 21-39 | - return ( ( abs(pdg_id()) >20 && abs(pdg_id()) <=40 ) - || abs(pdg_id()) == 9 ); - } - - /////////////////////////// - // INLINE Access Methods // - /////////////////////////// - - inline std::string ParticleData::name() const { return m_name; } - inline int ParticleData::pdg_id() const { return m_pdg_id; } - inline double ParticleData::charge() const { - return ( (double)m_3charge )/3.; - } - inline double ParticleData::mass() const { return m_mass; } - inline double ParticleData::clifetime() const { return m_clifetime; } - inline double ParticleData::spin() const { return m_2spin/2.; } - inline void ParticleData::set_charge( double charge ) { - if ( charge > 0 ) { - m_3charge = (int)(3.*charge+.1); - } else if ( charge < 0. ) { - m_3charge = (int)(3.*charge-.1); - } else { - m_3charge = 0; - } - } - inline void ParticleData::set_mass( double its_mass ) { - m_mass = its_mass; - } - inline void ParticleData::set_width( double width ) { - if ( width > 0 ) { - m_clifetime = HepMC_hbarc/width; - } else if ( width == 0. ) { - m_clifetime = -1.; - } else { - m_clifetime = 0.; - } - } - inline void ParticleData::set_clifetime( double its_clifetime ) { - m_clifetime = its_clifetime; - } - inline void ParticleData::set_spin( double spin ) { - m_2spin = (unsigned char)(2.*spin+.1); - } - - /////////////////////////// - // INLINE Operators // - /////////////////////////// - - inline bool ParticleData::operator==( const ParticleData& a ) const { - // compares everything except the particle's name - return ( a.m_pdg_id != m_pdg_id || - a.m_mass != m_mass || - a.m_clifetime != m_clifetime || - a.m_3charge != m_3charge || - a.m_2spin != m_2spin ) ? 0 : 1; - } - - inline bool ParticleData::operator!=( const ParticleData& a ) const { - // compares everything except the particle's name - return ( a.pdg_id() != this->pdg_id() ); - } - -} // HepMC - -#endif // HEPMC_PARTICLE_DATA_H -//-------------------------------------------------------------------------- - - - diff --git a/trunk/GeneratorInterface/PomwigInterface/src/ParticleDataTable.h b/trunk/GeneratorInterface/PomwigInterface/src/ParticleDataTable.h deleted file mode 100644 index f97e3d4991c2e..0000000000000 --- a/trunk/GeneratorInterface/PomwigInterface/src/ParticleDataTable.h +++ /dev/null @@ -1,233 +0,0 @@ -//-------------------------------------------------------------------------- -#ifndef HEPMC_PARTICLE_DATA_TABLE_H -#define HEPMC_PARTICLE_DATA_TABLE_H - -////////////////////////////////////////////////////////////////////////// -// Matt.Dobbs@Cern.CH, Jan 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 GenParticle Data Instances --- basically just an interface -// to STL map -- the same naming conventions are used -// A GenParticle may belong to any number of ParticleDataTables. -// The ParticleDataTable does not own the ParticleData objects and will NOT -// delete them unless explicity told to do so with the delete_all method. -// Each ParticleData entry in the table MUST have a unique pdg_id (otherwise -// an attempt to insert it as a new entry will fail). -// Updated 2000.02.08 M.Dobbs added merge_table and -// make_antiparticles_from_particles -////////////////////////////////////////////////////////////////////////// - -#include -#include -#include // needed for formatted output using sprintf -#include "ParticleData.h" - -namespace HepMC { - - class ParticleDataTable { - - public: - ParticleDataTable( std::string description = std::string() ); - ParticleDataTable( const char description ); - ParticleDataTable( const ParticleDataTable& ); - virtual ~ParticleDataTable(); // Shallow: does not delete - // ParticleData entries - // shallow: does not copy the entries, only makes new pointers - ParticleDataTable& operator=( const ParticleDataTable& ); - - void make_antiparticles_from_particles(); - int merge_table( const ParticleDataTable& ); - - void print( std::ostream& ostr = std::cout ) const; - - void delete_all(); //delete all ParticleData instances in this table - void clear(); //clears table without deleting - - ParticleData* operator[]( int id ) const; - ParticleData* find( int id ) const; - int size() const; - bool empty() const; - bool insert( ParticleData* ); // true if successful - bool erase( ParticleData* ); // removes from table - // does not delete - bool erase( int id ); // removes from table - // does not delete - typedef std::map::iterator iterator; - typedef std::map::const_iterator const_iterator; - iterator begin(); - iterator end(); - const_iterator begin() const; - const_iterator end() const; - - //////////////////// - // access methods // - //////////////////// - std::string description() const; - void set_description( std::string ); - void set_description( const char ); - - private: - std::string m_description; - std::map m_data_table; - }; - - /////////////////////////// - // INLINES // - /////////////////////////// - - inline ParticleDataTable::ParticleDataTable( std::string description ) - : m_description(description) {} - - inline ParticleDataTable::ParticleDataTable( const char description ) { - m_description = description; - } - - inline ParticleDataTable::ParticleDataTable( const ParticleDataTable& pdt){ - *this = pdt; - } - - inline ParticleDataTable::~ParticleDataTable(){} - - inline ParticleDataTable& ParticleDataTable::operator=( const - ParticleDataTable& - pdt) { - m_description = pdt.m_description; - m_data_table = pdt.m_data_table; - return *this; - } - - inline void ParticleDataTable::make_antiparticles_from_particles() { - ParticleDataTable new_data; - for ( ParticleDataTable::iterator p = begin(); p != end(); ++p ) { - ParticleData* pdata = p->second; - if ( pdata->charge() ) { - new_data.insert( new ParticleData( pdata->name()+"~", - -1*pdata->pdg_id(), - -1.*pdata->charge(), - pdata->mass(), - pdata->clifetime(), - pdata->spin() )); - } - } - merge_table( new_data ); - } - - inline void ParticleDataTable::print( std::ostream& ostr ) const { - // prints a summary of all particle Data currently in memory - // - ostr << "________________________________________" - << "________________________________________\n"; - ostr << "ParticleData: ***** ParticleDataTable" - << " ***** ( " << size() - << " entries )\n"; - ostr << " Description: " << m_description << "\n"; - ostr << " PDG ID " << " PARTICLE NAME " - << "CHARGE" << " MASS " - << " C*LIFETIME (CM) " << " SPIN\n"; - for ( std::map< int,ParticleData* >::const_iterator pd - = m_data_table.begin(); pd != m_data_table.end(); pd++ ) { - ostr << *(pd->second) << "\n"; - } - ostr << "________________________________________" - << "________________________________________" << std::endl; - } - - inline ParticleData* ParticleDataTable::find( int id ) const { - // finds a ParticleData pointer corresponding to id IF it exists in - // the table. If not returns NULL - std::map::const_iterator iter - = m_data_table.find(id); - return ( iter == m_data_table.end() ) ? 0 : iter->second; - } - - inline ParticleData* ParticleDataTable::operator[]( int id ) const { - return find(id); - } - - inline int ParticleDataTable::size() const { - return (int)m_data_table.size(); - } - - inline bool ParticleDataTable::empty() const { - return (bool)m_data_table.empty(); - } - - inline bool ParticleDataTable::insert( ParticleData* pdata ) { - // inserts pdata in the table IFF pdata's id has not already been used. - // It does NOT replace entries with the same id. True if successful. - // If you wish to overwrite another entry, first use erase() - if ( m_data_table.count(pdata->pdg_id()) ) return 0; - return ( m_data_table[pdata->pdg_id()] = pdata ); // true is success - } - - inline bool ParticleDataTable::erase( ParticleData* pdata ) { - // removes from table does not delete - // returns True is an entry pdata existed in the table and was erased - return (bool)m_data_table.erase( pdata->pdg_id() ); - } - - - inline bool ParticleDataTable::erase( int id ) { - // removes from table does not delete - // returns True is an entry pdata existed in the table and was erased - return (bool)m_data_table.erase( id ); - } - - inline ParticleDataTable::iterator ParticleDataTable::begin() { - return m_data_table.begin(); - } - - inline ParticleDataTable::iterator ParticleDataTable::end() { - return m_data_table.end(); - } - - inline ParticleDataTable::const_iterator ParticleDataTable::begin() const { - return m_data_table.begin(); - } - - inline ParticleDataTable::const_iterator ParticleDataTable::end() const { - return m_data_table.end(); - } - - inline std::string ParticleDataTable::description() const { - return m_description; - } - - inline void ParticleDataTable::set_description( std::string description ) { - m_description = description; - } - - inline void ParticleDataTable::set_description( const char description ) { - m_description = description; - } - - inline void ParticleDataTable::delete_all() { - // deletes all ParticleData instances in this table - for ( std::map::iterator pd = m_data_table.begin(); - pd != m_data_table.end(); pd++) delete pd->second; - clear(); - } - - inline void ParticleDataTable::clear() { m_data_table.clear(); } - - inline int ParticleDataTable::merge_table( const ParticleDataTable& pdt ) { - // merges pdt into this table - // each entry from pdt is inserted only if this table does not - // already have an entry matching the ParticleData's id - // returns the number of new entries inserted into this table. - int count_number_insertions =0; - for ( ParticleDataTable::const_iterator p = pdt.begin(); - p != pdt.end(); ++p ) { - if ( insert(p->second) ) ++count_number_insertions; - } - return count_number_insertions; - } - -} // HepMC - -#endif // HEPMC_PARTICLE_DATA_TABLE_H -//-------------------------------------------------------------------------- - - - diff --git a/trunk/GeneratorInterface/PomwigInterface/src/PomwigFilter.cc b/trunk/GeneratorInterface/PomwigInterface/src/PomwigFilter.cc new file mode 100644 index 0000000000000..196f80bd2e4ec --- /dev/null +++ b/trunk/GeneratorInterface/PomwigInterface/src/PomwigFilter.cc @@ -0,0 +1,35 @@ +#include "GeneratorInterface/PomwigInterface/interface/PomwigFilter.h" +#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h" + +PomwigFilter::PomwigFilter(const edm::ParameterSet& ppp) +{} + +PomwigFilter::~PomwigFilter() +{} + +bool +PomwigFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + using namespace edm; + + std::vector< Handle > AllProds; + iEvent.getManyByType(AllProds); + + if(AllProds.size()==0) { + std::cout<<" Event is skipped and removed." << std::endl; + return false; + } + else return true; +} + + +void +PomwigFilter::beginJob(const edm::EventSetup&) +{ +} + +void +PomwigFilter::endJob() { +} + +//define this as a plug-in diff --git a/trunk/GeneratorInterface/PomwigInterface/src/PomwigSource.cc b/trunk/GeneratorInterface/PomwigInterface/src/PomwigSource.cc index 4e8517468ef4a..b0fc6466d62b5 100644 --- a/trunk/GeneratorInterface/PomwigInterface/src/PomwigSource.cc +++ b/trunk/GeneratorInterface/PomwigInterface/src/PomwigSource.cc @@ -1,12 +1,14 @@ /* * Original Author: Fabian Stoeckli * 26/09/06 - * Modified for Pomwig interface - * 02/2007 + * Modified for Pomwig + * 03/2007 Antonio.Vilela.Pereira@cern.ch */ #include "GeneratorInterface/PomwigInterface/interface/PomwigSource.h" + #include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h" + #include "FWCore/Framework/interface/Event.h" #include "FWCore/ServiceRegistry/interface/Service.h" #include "FWCore/Utilities/interface/RandomNumberGenerator.h" @@ -20,11 +22,11 @@ using namespace edm; using namespace std; -#include "HerwigWrapper6_4.h" -#include "IO_HERWIG.h" -#include "HEPEVT_Wrapper.h" +#include "HepMC/HerwigWrapper6_4.h" +#include "HepMC/IO_HERWIG.h" +#include "HepMC/HEPEVT_Wrapper.h" -// INCLUDE JIMMY,HERWIG,LHAPDF COMMON BLOCKS AND FUNTIONS +// INCLUDE JIMMY,HERWIG,LHAPDF,POMWIG COMMON BLOCKS AND FUNTIONS #include "herwig.h" @@ -35,8 +37,15 @@ extern"C" { void setherwpdf_(void); // function to chatch 'STOP' in original HWWARN void cmsending_(int*); + + // struct to check wheter HERWIG killed an event + extern struct { + double eventisok; + } eventstat_; } +#define eventstat eventstat_ + #define setpdfpath setpdfpath_ #define mysetpdfpath mysetpdfpath_ #define setlhaparm setlhaparm_ @@ -60,10 +69,15 @@ PomwigSource::PomwigSource( const ParameterSet & pset, herwigLhapdfVerbosity_ (pset.getUntrackedParameter("herwigLhapdfVerbosity",0)), maxEventsToPrint_ (pset.getUntrackedParameter("maxEventsToPrint",0)), comenergy(pset.getUntrackedParameter("comEnergy",14000.)), - diffTopology(pset.getUntrackedParameter("diffTopology",0)), lhapdfSetPath_(pset.getUntrackedParameter("lhapdfSetPath","")), - printCards_(pset.getUntrackedParameter("printCards",true)) + printCards_(pset.getUntrackedParameter("printCards",true)), + diffTopology(pset.getParameter("diffTopology")) { + //Don't want this for Pomwig + useJimmy_ = false; + doMPInteraction_ = false; + numTrials_ = 0; + cout << "----------------------------------------------" << endl; cout << "Initializing PomwigSource" << endl; cout << "----------------------------------------------" << endl; @@ -82,8 +96,6 @@ PomwigSource::PomwigSource( const ParameterSet & pset, cout << " LHAPDF verbosity level = " << herwigLhapdfVerbosity_ << endl; cout << " HepMC verbosity = " << herwigHepMCVerbosity_ << endl; cout << " Number of events to be printed = " << maxEventsToPrint_ << endl; - - //Doesn't make sense for POMWIG /*if(useJimmy_) { cout << " HERWIG will be using JIMMY for UE/MI." << endl; if(doMPInteraction_) @@ -112,14 +124,14 @@ PomwigSource::PomwigSource( const ParameterSet & pset, hwbmch.PART2[0] = 'P'; hwbmch.PART2[1] = ' '; //hwpram.MODPDF[0] = -1; - //hwpram.MODPDF[1] = 20060; + //hwpram.MODPDF[1] = 20060; break; case 2: //SD survive PART2 hwbmch.PART1[0] = 'P'; hwbmch.PART1[1] = ' '; hwbmch.PART2[0] = 'E'; hwbmch.PART2[1] = '-'; - //hwpram.MODPDF[0] = 20060; + //hwpram.MODPDF[0] = 20060; //hwpram.MODPDF[1] = -1; break; case 3: //Non diffractive @@ -127,28 +139,26 @@ PomwigSource::PomwigSource( const ParameterSet & pset, hwbmch.PART1[1] = ' '; hwbmch.PART2[0] = 'P'; hwbmch.PART2[1] = ' '; - //hwpram.MODPDF[0] = 20060; - //hwpram.MODPDF[1] = 20060; + //hwpram.MODPDF[0] = 20060; + //hwpram.MODPDF[1] = 20060; break; default: throw edm::Exception(edm::errors::Configuration,"HerwigError") - <<" Invalid Diff. Topology. Must be DPE(diffTopology = 0), SD particle 1 (diffTopology = 1), SD particle 2 (diffTop -ology = 2) and Non diffractive (diffTopology = 3)"; + <<" Invalid Diff. Topology. Must be DPE(diffTopology = 0), SD particle 1 (diffTopology = 1), SD particle 2 (diffTopology = 2) and Non diffractive (diffTopology = 3)"; break; } for(int i=2;i<8;++i){ hwbmch.PART1[i] = ' '; hwbmch.PART2[i] = ' ';} - hwproc.MAXEV = pset.getUntrackedParameter("maxEvents",10); - hwevnt.MAXER = hwproc.MAXEV/10; + hwevnt.MAXER = pset.getUntrackedParameter("maxErrors",hwproc.MAXEV/10); + //hwproc.MAXEV = pset.getUntrackedParameter("maxEvents",10); + //hwevnt.MAXER = hwproc.MAXEV/10; if(hwevnt.MAXER<10) hwevnt.MAXER = 10; - //Not using JIMMY //if(useJimmy_) jmparm.MSFLAG = 1; // initialize other common blocks ... hwigin(); - //Not using JIMMY //if(useJimmy_) jimmin(); // set some 'non-herwig' defaults @@ -219,17 +229,22 @@ ology = 2) and Non diffractive (diffTopology = 3)"; // HERWIG preparations ... hwuinc(); hwusta("PI0 ",1); - // Initialize H1 pomeron structure function - /*int ifit = 5; - double xp = 0.1; - double Q2 = 75.0; - double xpq[13]; - qcd_1994(xp,Q2,xpq,ifit);*/ - hwabeg(); - + if(diffTopology != 3){ + int ifit = pset.getParameter("h1fit"); + if((ifit <= 0)||(ifit >= 7)){ + throw edm::Exception(edm::errors::Configuration,"HerwigError") + <<" Attempted to set non existant H1 fit index. Has to be 1...6"; + } + cout << " IFIT = "<< ifit < bare_product(new HepMCProduct()); - + + /*int counter = 0; + double mpiok = 1.0; + + while(mpiok > 0.5 && counter < numTrials_) { + // so far event is ok :) + eventstat.eventisok = 0.0; + + // call herwig routines to create HEPEVT + hwuine(); + hwepro(); + hwbgen(); + + // call jimmy ... only if event is not killed yet by HERWIG + if(useJimmy_ && doMPInteraction_ && (eventstat.eventisok < 0.5)) { + mpiok = hwmsct_dummy(1.1); + } + else mpiok = 0.0; + counter++; + } + + // event after numTrials MP is not ok -> skip event + if(mpiok > 0.5) { + cout<<" JIMMY could not produce MI in "< 0.5) return true; - }*/ + hwbgen(); hwdhob(); hwcfor(); @@ -277,11 +313,18 @@ bool PomwigSource::produce(Event & e) { hwdhvy(); hwmevt(); hwufne(); + + // if event was killed by HERWIG; skip + if(eventstat.eventisok > 0.5) return true; - HepMC::GenEvent* evt = new HepMC::GenEvent(); + // HEPEVT is ok, create new HepMC event + evt = new HepMC::GenEvent(); bool ok = conv.fill_next_event( evt ); + // if conversion failed; throw excpetion and stop processing if(!ok) throw cms::Exception("HerwigError") <<" Conversion problems in event nr."<set_signal_process_id(hwproc.IPROC); evt->set_event_number(numberEventsInRun() - remainingEvents() - 1); @@ -290,8 +333,14 @@ bool PomwigSource::produce(Event & e) { << "----------------------" << endl; evt->print(); } - if(evt) bare_product->addHepMCData(evt ); - e.put(bare_product); + + // dummy if: event MUST be there + if(evt) { + auto_ptr bare_product(new HepMCProduct()); + bare_product->addHepMCData(evt ); + e.put(bare_product); + } + return true; } diff --git a/trunk/GeneratorInterface/PomwigInterface/src/SealModule.cc b/trunk/GeneratorInterface/PomwigInterface/src/SealModule.cc index e5a3b20e26d82..c65c30f7c3848 100644 --- a/trunk/GeneratorInterface/PomwigInterface/src/SealModule.cc +++ b/trunk/GeneratorInterface/PomwigInterface/src/SealModule.cc @@ -2,9 +2,10 @@ #include "FWCore/Framework/interface/InputSourceMacros.h" #include "FWCore/Framework/interface/MakerMacros.h" #include "GeneratorInterface/PomwigInterface/interface/PomwigSource.h" +#include "GeneratorInterface/PomwigInterface/interface/PomwigFilter.h" using edm::PomwigSource; DEFINE_SEAL_MODULE(); DEFINE_ANOTHER_FWK_INPUT_SOURCE(PomwigSource); - +DEFINE_ANOTHER_FWK_MODULE(PomwigFilter); diff --git a/trunk/GeneratorInterface/PomwigInterface/src/h1init.f b/trunk/GeneratorInterface/PomwigInterface/src/h1init.f new file mode 100644 index 0000000000000..223ff8a6497db --- /dev/null +++ b/trunk/GeneratorInterface/PomwigInterface/src/h1init.f @@ -0,0 +1,14 @@ + SUBROUTINE H1INIT(IFIT) + INCLUDE 'HERWIG65.INC' + INTEGER IFIT + DOUBLE PRECISION XPQ(-6:6),XP,Q2 + +c---Initialization of the QCD_1994 subroutine + Q2=75D0 +C IFIT=5 + XP=1D-1 + CALL QCD_1994(XP,Q2,XPQ,IFIT) + + + END + diff --git a/trunk/GeneratorInterface/PomwigInterface/src/herwig.h b/trunk/GeneratorInterface/PomwigInterface/src/herwig.h index 4ff28ce7bbba4..5f8b0b4e43148 100644 --- a/trunk/GeneratorInterface/PomwigInterface/src/herwig.h +++ b/trunk/GeneratorInterface/PomwigInterface/src/herwig.h @@ -324,17 +324,23 @@ extern"C" { //------------------------------ LHAPDF functions ------------------------------------------------- - //------------------------------ POMWIG functions ------------------------------------------------- // Subroutine inside H1QCD #define qcd_1994 qcd_1994_ extern "C" { - void qcd_1994(double&,double&,double*,int); + void qcd_1994(double&,double&,double*,int&); } // HWABEG to initialize H1 pomeron structure function #define hwabeg hwabeg_ extern "C" { void hwabeg(); } +// H1INIT initializes H1 structure function w/ IFIT as argument +#define h1init h1init_ +extern "C" { + void h1init(int&); +} + + #endif diff --git a/trunk/GeneratorInterface/PomwigInterface/src/hwabeg.f b/trunk/GeneratorInterface/PomwigInterface/src/hwabeg.f index e29c1dccb6d70..d8cf81417defe 100644 --- a/trunk/GeneratorInterface/PomwigInterface/src/hwabeg.f +++ b/trunk/GeneratorInterface/PomwigInterface/src/hwabeg.f @@ -1,7 +1,7 @@ C * 07/03/2003, Tibor Kucs C * User's routine for initialization C----------------------------------------------------------------------- - SUBROUTINE HWABEG + SUBROUTINE HWABEG() INCLUDE 'HERWIG65.INC' INTEGER IFIT DOUBLE PRECISION XPQ(-6:6),XP,Q2 diff --git a/trunk/GeneratorInterface/PomwigInterface/test/PomwigAnalyzer.cc b/trunk/GeneratorInterface/PomwigInterface/test/PomwigAnalyzer.cc index b28d4d0cf2c62..a7f1c85afd802 100644 --- a/trunk/GeneratorInterface/PomwigInterface/test/PomwigAnalyzer.cc +++ b/trunk/GeneratorInterface/PomwigInterface/test/PomwigAnalyzer.cc @@ -14,8 +14,10 @@ #include "FWCore/ParameterSet/interface/ParameterSet.h" -#include "CLHEP/HepMC/GenEvent.h" -#include "CLHEP/HepMC/GenParticle.h" +#include "HepMC/GenEvent.h" +#include "HepMC/GenParticle.h" + +//#include "CLHEP/Vector/LorentzVector.h" #include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h" @@ -51,17 +53,19 @@ PomwigAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) //std::vector protons; HepMC::GenParticle* proton1 = 0; HepMC::GenParticle* proton2 = 0; + double pz1max = 0.; + double pz2min = 0.; for(HepMC::GenEvent::particle_iterator it = evt->particles_begin(); it != evt->particles_end(); ++it) { - if(((*it)->pdg_id() == 2212)&&((*it)->status() == 1)&&((*it)->momentum().pz() > 5200.)){ - proton1 = *it; - } else if(((*it)->pdg_id() == 2212)&&((*it)->status() == 1)&&((*it)->momentum().pz() < -5200.)){ - proton2 = *it; + double pz = (*it)->momentum().pz(); + if(((*it)->pdg_id() == 2212)&&((*it)->status() == 1)&&(pz > 5200.)){ + if(pz > pz1max){proton1 = *it;pz1max=pz;} + } else if(((*it)->pdg_id() == 2212)&&((*it)->status() == 1)&&(pz < -5200.)){ + if(pz < pz2min){proton2 = *it;pz2min=pz;} } } if(proton1){ std::cout << "Proton 1: " << proton1->momentum().perp() << " " << proton1->momentum().pseudoRapidity() << " " << proton1->momentum().phi() << std::endl; - HepLorentzVector proton1in(0.,0.,7000.,7000.); - double t1 = (proton1->momentum() - proton1in).m2(); + double t1 = -proton1->momentum().perp2(); double xigen1 = 1 - proton1->momentum().pz()/7000.; hist_t->Fill(t1); hist_xigen->Fill(xigen1); @@ -69,8 +73,7 @@ PomwigAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) if(proton2){ std::cout << "Proton 2: " << proton2->momentum().perp() << " " << proton2->momentum().pseudoRapidity() << " " << proton2->momentum().phi() << std::endl; - HepLorentzVector proton2in(0.,0.,-7000.,7000.); - double t2 = (proton2->momentum() - proton2in).m2(); + double t2 = -proton2->momentum().perp2(); double xigen2 = 1 + proton2->momentum().pz()/7000.; hist_t->Fill(t2); hist_xigen->Fill(xigen2); diff --git a/trunk/GeneratorInterface/PomwigInterface/test/PomwigAnalyzer.h b/trunk/GeneratorInterface/PomwigInterface/test/PomwigAnalyzer.h index a1ba6228b1f5b..6842803728b98 100644 --- a/trunk/GeneratorInterface/PomwigInterface/test/PomwigAnalyzer.h +++ b/trunk/GeneratorInterface/PomwigInterface/test/PomwigAnalyzer.h @@ -9,9 +9,9 @@ #include "FWCore/ParameterSet/interface/ParameterSet.h" -#include "CLHEP/HepMC/WeightContainer.h" -#include "CLHEP/HepMC/GenEvent.h" -#include "CLHEP/HepMC/GenParticle.h" +#include "HepMC/WeightContainer.h" +#include "HepMC/GenEvent.h" +#include "HepMC/GenParticle.h" #include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h" diff --git a/trunk/GeneratorInterface/PomwigInterface/test/Z2muAnalyzer.cc b/trunk/GeneratorInterface/PomwigInterface/test/Z2muAnalyzer.cc index 2d5fbc940e1a7..d01859433562b 100644 --- a/trunk/GeneratorInterface/PomwigInterface/test/Z2muAnalyzer.cc +++ b/trunk/GeneratorInterface/PomwigInterface/test/Z2muAnalyzer.cc @@ -1,7 +1,7 @@ // // Original Author: Fabian Stoeckli // Created: Tue Nov 14 13:43:02 CET 2006 -// $Id: H4muAnalyzer.cc,v 1.2 2007/02/14 15:51:35 fabstoec Exp $ +// $Id: Z2muAnalyzer.cc,v 1.1 2007/03/07 17:05:31 antoniov Exp $ // // Modified for PomwigInterface test for Z/gamma* -> 2mu // 02/2007 @@ -24,8 +24,10 @@ #include "FWCore/ParameterSet/interface/ParameterSet.h" -#include "CLHEP/HepMC/GenEvent.h" -#include "CLHEP/HepMC/GenParticle.h" +#include "HepMC/GenEvent.h" +#include "HepMC/GenParticle.h" + +#include "CLHEP/Vector/LorentzVector.h" #include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h" @@ -65,12 +67,13 @@ Z2muAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) // if there are at least two muons // calculate invarant mass of first two and fill it into histogram - HepLorentzVector tot_momentum; double inv_mass = 0.0; std::cout<=2) { - tot_momentum = muons[0]->momentum(); - tot_momentum += muons[1]->momentum(); + HepLorentzVector tot_momentum(muons[0]->momentum().px() + muons[1]->momentum().px(), + muons[0]->momentum().py() + muons[1]->momentum().py(), + muons[0]->momentum().pz() + muons[1]->momentum().pz(), + muons[0]->momentum().e() + muons[1]->momentum().e()); inv_mass = sqrt(tot_momentum.m2()); } diff --git a/trunk/GeneratorInterface/PomwigInterface/test/Z2muAnalyzer.h b/trunk/GeneratorInterface/PomwigInterface/test/Z2muAnalyzer.h index 27170d7782aa4..a0d4515f7e9a7 100644 --- a/trunk/GeneratorInterface/PomwigInterface/test/Z2muAnalyzer.h +++ b/trunk/GeneratorInterface/PomwigInterface/test/Z2muAnalyzer.h @@ -9,9 +9,9 @@ #include "FWCore/ParameterSet/interface/ParameterSet.h" -#include "CLHEP/HepMC/WeightContainer.h" -#include "CLHEP/HepMC/GenEvent.h" -#include "CLHEP/HepMC/GenParticle.h" +#include "HepMC/WeightContainer.h" +#include "HepMC/GenEvent.h" +#include "HepMC/GenParticle.h" #include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"