generator.cc

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////
// Matt.Dobbs@Cern.CH, October 2002
// example of generating events with Herwig using HepMC/HerwigWrapper.h
// Events are read into the HepMC event record from the FORTRAN HEPEVT
// common block using the IO_HERWIG strategy.
//////////////////////////////////////////////////////////////////////////

#include <stdlib.h>
#include <time.h>

#include <stdio.h>
#include <iostream>
#include <fstream>
#include <sstream>

#include "MyHerwigWrapper.h"
#include "HepMC/IO_HERWIG.h"
//#include "HepMC/IO_Ascii.h"   
#include "HepMC/GenEvent.h"
#include "HepMC/HEPEVT_Wrapper.h"
#include "HepMC/SimpleVector.h"
#include "HepMC/IO_GenEvent.h"  //ap add


// include our own classes
#include "Configuration.h"
#include "readherwigConfigFile.h"

using namespace std;

// a dummy routine called on the termination of the run
extern "C" void hwaend_ (void); 
extern "C" {
  int hwgethepevtsize_(void);
  void hwsetmodbos_(int& i, int& ivalue);
  double hwgetcs(void);
  double hwgetcserr(void);
}

// ---------------------------------------------------------------------- 
int main(int argc, const char* argv[]) 
{
  if(argc <= 2 ) { //ap
    std::cout << "Usage: " << argv[0] << " <Configuration File>" << " <name of file to generate to>"<< std::endl; //ap
    return 0;
  }

  Configuration* currentconfig = new Configuration();
  if ( currentconfig->read( argv[ 1 ] ) ) {
    cout << " ... skip !" << endl;
    return 0;
  }

 

  //
  //........................................HEPEVT
  // herwig-6510-2 uses HEPEVT with 10000 entries and 8-byte floating point
  //  numbers. We need to explicitly pass this information to the 
  //  HEPEVT_Wrapper.
  //
  HepMC::HEPEVT_Wrapper::set_max_number_entries(hwgethepevtsize_());
  HepMC::HEPEVT_Wrapper::set_sizeof_real(8);

  //random seed
  hwevnt.nrn[0] = currentconfig->rseed;
  
  // number of events
  hwproc.maxev = int(currentconfig->nevents); 
  
  // tell it what the beam particles are:
  for ( unsigned int i = 0; i < 8; ++i ) {
    hwbmch.part1[i] = (i < 1) ? 'P' : ' ';
    hwbmch.part2[i] = (i < 1) ? 'P' : ' ';
  }

  // INITIALISE OTHER COMMON BLOCKS
  hwigin();
  jimmin();

  cout << "Configuration Parameter List:" << endl;
  for (unsigned int files=0; files<(currentconfig->ConfigFileNames).size(); files++)
    readConfigFile((currentconfig->ConfigFileNames)[files]);

  // master switches for ISR
  // true switches ISR off!
  //  hwpram.nospac=int(currentconfig->ISR); //switch on/off ISR
    
  // number of events to print
  hwevnt.maxpr = 1; 
  // compute parameter-dependent constants
  hwuinc(); 
  // initialise elementary process
  hweini(); 
  // initialise jimmy event
  jminit();
  //........................................HepMC INITIALIZATIONS
  //
  // Instantiate an IO strategy for reading from HEPEVT.
  HepMC::IO_HERWIG hepevtio;
  //
  //........................................EVENT LOOP
  HepMC::IO_GenEvent ascii_out(argv[2],std::ios::out);  //add ap
  for ( int iEvent = 1; iEvent <= hwproc.maxev; iEvent++ ) {
    if ( iEvent%50==1 ) cout << "Processing Event Number " 
                                  << iEvent << endl;
//     // initialise event
//     hwuine();
//     // generate hard subproces
//     hwepro();
//     // generate parton cascades
//     hwbgen();
//     // do heavy object decays
//     hwdhob();
//     // do cluster formation
//     hwcfor();
//     // do cluster decays
//     hwcdec();
//     // do unstable particle decays
//     hwdhad();
//     // do heavy flavour hadron decays
//     hwdhvy();
//     // add soft underlying event if needed
//     hwmevt();
//     // finish event
//     hwufne();

    // initialise event
    hwuine();
    // generate hard subproces
    hwepro();
    // generate parton cascades
    hwbgen();

    // Do Jimmy underlying event if required.
    int abort = 0;
    hwmsct(&abort);
    if (abort != 0) {
      cout << "Event aborted by hwmsct." << endl;
      hwufne();  // finish event
      continue;
    }

    //    if (abort == 0) {
      hwdhob();  // Do heavy quark decays
      hwcfor();  // Do cluster formation
      hwcdec();  // Do cluster decays
      hwdhad();  // Do unstable particle decays
      hwdhvy();  // Do heavy flavor decays
      hwmevt();  // Add soft underlying event
      hwufne();  // finish event
    //}
      if (hwevnt.ierror==-1 || hwevnt.ierror>99) {
        cout << "Event Error, Number: " << " Code: " << hwevnt.ierror << endl;
        continue;
      }
    
    // Print HEPEVT format event here 
    // HepMC::HEPEVT_Wrapper::print_hepevt();

    HepMC::GenEvent* evt = hepevtio.read_next_event();
    // add some information to the event
    evt->set_event_number(iEvent);
    evt->set_signal_process_id(20);
   
    // the status codes should not be changed, otherwise you change
    // the output of the MC generators (status code variabl) and this
    // should not be done
    //     //change the status code in HepMC format event.  
    //     for ( HepMC::GenEvent::particle_const_iterator p =
    //      evt->particles_begin(); p != evt->particles_end(); ++p ){
    //       if( abs((*p)->status())>1 )
    //  (*p)->set_status(2);
    //     }
       
    ascii_out << evt;    //add ap
      
    // we also need to delete the created event from memory
    delete evt;

  }
  //........................................TERMINATION
  hwefin();

 
  // clean up config object
  delete currentconfig;

 
  cout<<"Run successfully!"<<endl;
  return 0;
}

// a dummy user-supplied subroutine called when the run terminates:
void hwaend_ (void) { return; }

---