CMSSW Reference Manual

00001 
00019 #include "GeneratorInterface/MadGraphInterface/interface/MadGraphProducer.h"
00020 #include "GeneratorInterface/MadGraphInterface/interface/PYR.h"
00021 #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
00022 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
00023 #include "SimDataFormats/GeneratorProducts/interface/LHECommonBlocks.h"
00024 #include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
00025 #include "FWCore/Framework/interface/Event.h"
00026 #include "FWCore/Framework/interface/Run.h"
00027 #include "FWCore/ServiceRegistry/interface/Service.h"
00028 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
00029 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00030 
00031 #include <algorithm>
00032 #include <iterator>
00033 #include <iostream>
00034 #include <fstream> 
00035 #include <cstring>
00036 #include <cstdio>
00037 #include <string>
00038 
00039 #include "time.h"
00040 
00041 // Generator modifications
00042 #include "HepMC/PythiaWrapper6_2.h"
00043 //#include "CLHEP/HepMC/ConvertHEPEVT.h"
00044 #include "HepMC/IO_HEPEVT.h"
00045 //#include "CLHEP/HepMC/CBhepevt.h"
00046 
00047 // MCDB Interface 
00048 #include "GeneratorInterface/MadGraphInterface/interface/MCDBInterface.h"
00049 
00050 #define PYGIVE pygive_
00051 extern "C" {
00052   void PYGIVE(const char*,int length);
00053 }
00054 
00055 extern "C"{
00056  void eventtree_();
00057 }
00058 
00059 #define TXGIVE txgive_
00060 extern "C" {
00061   void TXGIVE(const char*,int length);
00062 }
00063 
00064 #define TXGIVE_INIT txgive_init_
00065 extern "C" {
00066   void TXGIVE_INIT();
00067 }
00068 
00069 /*
00070 ME2Pythia.f
00071       double precision etcjet,rclmax,etaclmax,qcut,clfact
00072       integer maxjets,minjets,iexcfile,ktsche
00073       common/MEMAIN/etcjet,rclmax,etaclmax,qcut,clfact,
00074      $   maxjets,minjets,iexcfile,ktsche
00075       DATA ktsche/0/
00076       DATA qcut,clfact/0d0,0d0/
00077 */
00078 extern "C" {
00079  extern struct MEMAIN{
00080  double etcjet;
00081  double rclmax;
00082  double etaclmax;
00083  double qcut;
00084  double clfact;
00085  int maxjets;
00086  int minjets;
00087  int iexcfile;
00088  int ktsche;
00089  }memain_;
00090 }
00091 
00092 /*
00093       LOGICAL minimalLH,externalLH,validLH
00094       common /SOURCEPRS/minimalLH,externalLH,validLH
00095 */
00096 extern "C"{
00097   extern struct SOURCEPRS {
00098     int minimalLH;
00099     int externalLH;
00100     int validLH;
00101   } sourceprs_;
00102 }
00103 
00104 // init LHNIN Fortran unit with file "fileName" to read.
00105 extern "C" {
00106   void mgopen_(const char *fname, int len);
00107   void mgclos_();
00108 }
00109 
00110 // turn ME2pythia.f aborts into an exception, shutting down EDM cleanly
00111 extern "C"{
00112   void pystop_(int *mcod)
00113   {
00114     throw cms::Exception("Generator|MadGraphProducer")
00115           << "PYSTOP called with code: " << *mcod << std::endl;
00116   }
00117 }
00118 
00119 //used for defaults - change these to defines? TODO
00120   static const unsigned long kNanoSecPerSec = 1000000000;
00121   static const unsigned long kAveEventPerSec = 200;
00122 
00123 using namespace edm;
00124 using namespace std;
00125 
00126 MadGraphProducer::MadGraphProducer( const ParameterSet & pset) :
00127  EDFilter(), evt(0),
00128  pythiaPylistVerbosity_ (pset.getUntrackedParameter<int>("pythiaPylistVerbosity",0)),
00129  pythiaHepMCVerbosity_ (pset.getUntrackedParameter<bool>("pythiaHepMCVerbosity",false)),
00130  maxEventsToPrint_ (pset.getUntrackedParameter<int>("maxEventsToPrint",0)),
00131  firstEvent_(pset.getUntrackedParameter<unsigned int>("firstEvent", 0)),
00132  MEMAIN_etaclmax(pset.getUntrackedParameter<double>("MEMAIN_etaclmax",0.)),
00133  MEMAIN_qcut(pset.getUntrackedParameter<double>("MEMAIN_qcut",0.)),
00134  MEMAIN_iexcfile(pset.getUntrackedParameter<unsigned int>("MEMAIN_iexcfile",0)),
00135  produceEventTreeFile_ (pset.getUntrackedParameter<bool>("produceEventTreeFile",false)),
00136  minimalLH_(pset.getUntrackedParameter<bool>("minimalLH",false)),
00137  initialized_(false),
00138  eventNumber_(firstEvent_),
00139  pset_(pset),
00140  useExternalGenerators_(false),
00141  useTauola_(false),
00142  useTauolaPolarization_(false)
00143 {
00144   edm::LogInfo("Generator|MadGraph")<<" initializing MadGraphProducer";
00145   pdf_info = new HepMC::PdfInfo();
00146 
00147   // check whether using minimal LH
00148   sourceprs_.minimalLH = minimalLH_; // pass to ME2pythia through common block
00149   sourceprs_.externalLH = true;
00150   if(minimalLH_) edm::LogInfo("Generator|MadGraph")<<" ----- Using minimal Les Houches Accord functionality - ignoring MadGraph specifics.";
00151   
00152   // first set to default values, mostly zeros
00153   memain_.etcjet=0.;
00154   memain_.rclmax=0.;
00155   memain_.etaclmax=0.;
00156   memain_.qcut=0.;
00157   memain_.clfact=0.;
00158   memain_.maxjets=0;
00159   memain_.minjets=0;
00160   memain_.iexcfile=0;
00161   memain_.ktsche=0;
00162   // then set (some) values from cards
00163   memain_.iexcfile=MEMAIN_iexcfile;
00164   memain_.etaclmax=MEMAIN_etaclmax;
00165   memain_.qcut=MEMAIN_qcut;
00166   // print out
00167   edm::LogInfo("Generator|MadGraph")<<"MEMAIN before ME2pythia initialization - etcjet ="<<memain_.etcjet<<" rclmax ="<<memain_.rclmax<<" etaclmax ="<<memain_.etaclmax<<" qcut ="<<memain_.qcut<<" clfact ="<<memain_.clfact<<" maxjets ="<<memain_.maxjets<<" minjets ="<<memain_.minjets<<" iexcfile ="<<memain_.iexcfile<<" ktsche ="<<memain_.ktsche;
00168 
00169   edm::LogInfo("Generator|MadGraph")<<"MadGraphProducer: initializing Pythia.";
00170   // PYLIST Verbosity Level
00171   // Valid PYLIST arguments are: 1, 2, 3, 5, 7, 11, 12, 13
00172   edm::LogInfo("Generator|MadGraph")<<"Pythia PYLIST verbosity level = " << pythiaPylistVerbosity_;
00173   // HepMC event verbosity Level
00174   edm::LogInfo("Generator|MadGraph")<<"Pythia HepMC verbosity = " << pythiaHepMCVerbosity_;
00175   //Max number of events printed on verbosity level 
00176   edm::LogInfo("Generator|MadGraph")<<"Number of events to be printed = " << maxEventsToPrint_;
00177   // max nr of events and first event
00178 //edm::LogInfo("Generator|MadGraph")<<"firstEvent / maxEvents = "<<firstEvent_<<" / "<< maxEvents();
00179   edm::LogInfo("Generator|MadGraph")<<"firstEvent / maxEvents = "<<firstEvent_<<" / 999999999";
00180 
00181   // Set PYTHIA parameters in a single ParameterSet
00182   ParameterSet pythia_params = pset.getParameter<ParameterSet>("PythiaParameters") ;
00183   
00184   // The parameter sets to be read (default, min bias, user ...) in the proper order.
00185   std::vector<std::string> setNames = pythia_params.getParameter<std::vector<std::string> >("parameterSets");
00186   
00187   // Loop over the sets
00188   for ( unsigned i=0; i<setNames.size(); ++i ) {
00189     std::string mySet = setNames[i];
00190     
00191     // Read the PYTHIA parameters for each set of parameters
00192     std::vector<std::string> pars = 
00193       pythia_params.getParameter<std::vector<std::string> >(mySet);
00194     
00195     edm::LogInfo("Generator|MadGraph")<<"----------------------------------------------";
00196     edm::LogInfo("Generator|MadGraph")<<"Read PYTHIA parameter set " << mySet;
00197     edm::LogInfo("Generator|MadGraph")<<"----------------------------------------------";
00198     
00199     // Loop over all parameters and stop in case of mistake
00200     for( std::vector<std::string>::const_iterator  
00201            itPar = pars.begin(); itPar != pars.end(); ++itPar ) {
00202       static std::string sRandomValueSetting("MRPY(1)");
00203       if( 0 == itPar->compare(0,sRandomValueSetting.size(),sRandomValueSetting) ) {
00204         throw edm::Exception(edm::errors::Configuration,"PythiaError")
00205           <<" attempted to set random number using pythia command 'MRPY(1)' this is not allowed.\n  Please use the RandomNumberGeneratorService to set the random number seed.";
00206       }
00207       if( ! call_pygive(*itPar) ) {
00208         throw edm::Exception(edm::errors::Configuration,"PythiaError") 
00209           <<" pythia did not accept the following \""<<*itPar<<"\"";
00210       }
00211     }
00212   }
00213 
00214   // TAUOLA, etc.
00215   //
00216   useExternalGenerators_ = pset.getUntrackedParameter<bool>("UseExternalGenerators",false);
00217   
00218   //Setting random numer seed
00219   edm::LogInfo("Generator|MadGraph")<<"----------------------------------------------";
00220   edm::LogInfo("Generator|MadGraph")<<"Setting Pythia random number seed";
00221   edm::LogInfo("Generator|MadGraph")<<"----------------------------------------------";
00222   edm::Service<RandomNumberGenerator> rng;
00223   randomEngine = fRandomEngine = &(rng->getEngine());
00224   uint32_t seed = rng->mySeed();
00225   std::ostringstream sRandomSet;
00226   sRandomSet <<"MRPY(1)="<<seed;
00227   call_pygive(sRandomSet.str());
00228   produces<HepMCProduct>();
00229   edm::LogInfo("Generator|MadGraph")<<"starting event generation ...";
00230 }
00231 
00232 
00233 MadGraphProducer::~MadGraphProducer(){
00234   edm::LogInfo("Generator|MadGraph")<<"event generation done.";
00235   delete pdf_info;
00236   clear();
00237   //  rm -f MGinput.dat;
00238 }
00239 
00240 void MadGraphProducer::clear() {
00241 }
00242 
00243 void MadGraphProducer::init() {
00244   initialized_ = true;
00245 
00246   // Call pythia initialisation with user defined upinit subroutine
00247   call_pyinit( "USER", "", "", 0.);
00248 
00249   if ( useExternalGenerators_ ) {
00250     // read External Generator parameters
00251     ParameterSet ext_gen_params =
00252        pset_.getParameter<ParameterSet>("ExternalGenerators") ;
00253     vector<string> extGenNames =
00254        ext_gen_params.getParameter< vector<string> >("parameterSets");
00255     for (unsigned int ip=0; ip<extGenNames.size(); ++ip )
00256     {
00257       string curSet = extGenNames[ip];
00258       ParameterSet gen_par_set =
00259          ext_gen_params.getUntrackedParameter< ParameterSet >(curSet);
00260 /*
00261      cout << "----------------------------------------------" << endl;
00262      cout << "Read External Generator parameter set "  << endl;
00263      cout << "----------------------------------------------" << endl;
00264 */
00265      if ( curSet == "Tauola" )
00266      {
00267         useTauola_ = true;
00268         if ( useTauola_ ) {
00269            cout << "--> use TAUOLA" << endl;
00270         } 
00271         useTauolaPolarization_ = gen_par_set.getParameter<bool>("UseTauolaPolarization");
00272         if ( useTauolaPolarization_ ) 
00273         {
00274            cout << "(Polarization effects enabled)" << endl;
00275            tauola_.enablePolarizationEffects();
00276         } 
00277         else 
00278         {
00279            cout << "(Polarization effects disabled)" << endl;
00280            tauola_.disablePolarizationEffects();
00281         }
00282         vector<string> cards = gen_par_set.getParameter< vector<string> >("InputCards");
00283         cout << "----------------------------------------------" << endl;
00284         cout << "Initializing Tauola" << endl;
00285         for( vector<string>::const_iterator
00286                 itPar = cards.begin(); itPar != cards.end(); ++itPar )
00287         {
00288            // call_txgive(*itPar);
00289            TXGIVE( (*itPar).c_str(), (*itPar).length() );
00290            cout << "     " <<  (*itPar).c_str() << endl;
00291         }
00292         tauola_.initialize();
00293         //call_pretauola(-1); // initialize TAUOLA package for tau decays
00294      }
00295     }
00296     // cout << "----------------------------------------------" << endl;
00297   }
00298 
00299 
00300   edm::LogInfo("Generator|MadGraph")<<"MEMAIN after ME2pythia initialization - etcjet ="<<memain_.etcjet<<" rclmax ="<<memain_.rclmax<<" etaclmax ="<<memain_.etaclmax<<" qcut ="<<memain_.qcut<<" clfact ="<<memain_.clfact<<" maxjets ="<<memain_.maxjets<<" minjets ="<<memain_.minjets<<" iexcfile ="<<memain_.iexcfile<<" ktsche ="<<memain_.ktsche;
00301 }
00302 
00303 bool MadGraphProducer::beginRun(Run& run, const EventSetup& es)
00304 {
00305   edm::Handle<LHERunInfoProduct> product;
00306   run.getByLabel("source", product);
00307   lhef::CommonBlocks::fillHEPRUP(&product->heprup());
00308 
00309   if (!minimalLH_) {
00310     const char *fname = std::tmpnam(NULL);
00311     std::ofstream file(fname, std::fstream::out | std::fstream::trunc);
00312     std::copy(product->begin(), product->init(),
00313               std::ostream_iterator<std::string>(file));
00314     file.close();
00315     mgopen_(fname, std::strlen(fname));
00316     std::remove(fname);
00317   }
00318 
00319   return true;
00320 }
00321 
00322 bool MadGraphProducer::endRun(Run& run, const EventSetup& es)
00323 {
00324   call_pystat(1);
00325   if ( useTauola_ ) {
00326     tauola_.print();
00327     //call_pretauola(1); // print TAUOLA decay statistics output
00328   }
00329 
00330   mgclos_();
00331   initialized_ = false;
00332   return true;
00333 }
00334 
00335 bool MadGraphProducer::filter(Event & e, const EventSetup& es) 
00336 {
00337   edm::Handle<LHEEventProduct> product;
00338   e.getByLabel("source", product);
00339   lhef::CommonBlocks::fillHEPEUP(&product->hepeup());
00340   sourceprs_.validLH = true;
00341   if (!initialized_)
00342     init();
00343 
00344   std::auto_ptr<HepMCProduct> bare_product(new HepMCProduct());  
00345 
00346   call_pyevnt();      // generate one event with Pythia
00347 
00348   if (hepeup_.nup == 0 || pypars.msti[0] == 1) {
00349     // the event got rejected by the matching and we have to filter
00350     edm::LogInfo("Generator|MadGraph") << "Event was skipped.";
00351     e.put(bare_product); // put an empty product, CMSSW filter will take care
00352     return false;
00353   }
00354 
00355   if ( useTauola_ ) {
00356     tauola_.processEvent();
00357     //call_pretauola(0); // generate tau decays with TAUOLA
00358   }
00359 
00360   call_pyhepc( 1 );
00361 
00362   if(produceEventTreeFile_) eventtree_(); // write out an event tree file
00363 
00364   HepMC::IO_HEPEVT conv;
00365   HepMC::GenEvent* evt = conv.read_next_event();
00366   evt->set_signal_process_id(pypars.msti[0]);
00367   evt->set_event_number(eventNumber_);
00368   ++eventNumber_;
00369 
00370   // store PDF information
00371   int id_1 = pypars.msti[14];
00372   int id_2 = pypars.msti[15];
00373   if (id_1 == 21) id_1 = 0;
00374   if (id_2 == 21) id_2 = 0;
00375   pdf_info->set_id1(id_1);
00376   pdf_info->set_id2(id_2);
00377   pdf_info->set_x1(pypars.pari[32]);
00378   pdf_info->set_x2(pypars.pari[33]);
00379   pdf_info->set_scalePDF(pypars.pari[20]);
00380   pdf_info->set_pdf1(pypars.pari[28]);
00381   pdf_info->set_pdf2(pypars.pari[29]);
00382   evt->set_pdf_info( *pdf_info);
00383 
00384   if(e.id().event() <= maxEventsToPrint_ && (pythiaPylistVerbosity_ || pythiaHepMCVerbosity_)) {
00385     // Prints PYLIST info
00386     if(pythiaPylistVerbosity_) {
00387       call_pylist(pythiaPylistVerbosity_);
00388     }
00389     // Prints HepMC event
00390     if(pythiaHepMCVerbosity_) {
00391       edm::LogInfo("Generator|MadGraph")<<"Event process = " << pypars.msti[0];
00392       evt->print();
00393     }
00394   }
00395 
00396   bare_product->addHepMCData(evt );
00397   e.put(bare_product);
00398 
00399   return true;
00400 }
00401 
00402 bool MadGraphProducer::call_pygive(const std::string& iParm ) {
00403   int numWarn = pydat1.mstu[26]; //# warnings
00404   int numErr = pydat1.mstu[22];// # errors
00405   //call the fortran routine pygive with a fortran string
00406   PYGIVE( iParm.c_str(), iParm.length() );
00407   //  PYGIVE( iParm );
00408   //if an error or warning happens it is problem
00409   return pydat1.mstu[26] == numWarn && pydat1.mstu[22] == numErr;
00410 }
MadGraphProducer.cc
