CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/GeneratorInterface/SherpaInterface/src/SherpaHadronizer.cc

Go to the documentation of this file.
00001 #include <iostream>
00002 #include <sstream>
00003 #include <string>
00004 #include <memory>
00005 #include <stdint.h>
00006 
00007 
00008 #include "HepMC/GenEvent.h"
00009 
00010 #include "SHERPA/Main/Sherpa.H"
00011 #include "ATOOLS/Org/Message.H"
00012 #include "ATOOLS/Math/Random.H"
00013 #include "ATOOLS/Org/Exception.H"
00014 #include "ATOOLS/Org/Run_Parameter.H"
00015 #include "ATOOLS/Org/MyStrStream.H"
00016 #include "SHERPA/Tools/Input_Output_Handler.H"
00017 #include "SHERPA/Tools/HepMC2_Interface.H"
00018 
00019 #include "GeneratorInterface/Core/interface/ParameterCollector.h"
00020 #include "GeneratorInterface/Core/interface/BaseHadronizer.h"
00021 #include "GeneratorInterface/Core/interface/GeneratorFilter.h"
00022 #include "GeneratorInterface/Core/interface/HadronizerFilter.h"
00023 #include "GeneratorInterface/Core/interface/RNDMEngineAccess.h"
00024 
00025 
00026 class SherpaHadronizer : public gen::BaseHadronizer {
00027 public:
00028   SherpaHadronizer(const edm::ParameterSet &params);
00029   ~SherpaHadronizer();
00030   
00031   bool readSettings( int ) { return true; }
00032   bool initializeForInternalPartons();
00033   bool declareStableParticles(const std::vector<int> &pdgIds);
00034   bool declareSpecialSettings( const std::vector<std::string> ) { return true; }
00035   void statistics();
00036   bool generatePartonsAndHadronize();
00037   bool decay();
00038   bool residualDecay();
00039   void finalizeEvent();
00040   const char *classname() const { return "SherpaHadronizer"; }
00041   
00042 private:
00043   
00044   std::string SherpaPath;
00045   std::string SherpaPathPiece;
00046   std::string SherpaResultDir;
00047   double default_weight;
00048   edm::ParameterSet     SherpaParameter;
00049   unsigned int  maxEventsToPrint;
00050   
00051   SHERPA::Sherpa Generator;
00052   CLHEP::HepRandomEngine* randomEngine;
00053   
00054 };
00055 
00056 class CMS_SHERPA_RNG: public ATOOLS::External_RNG {
00057 public:
00058   CMS_SHERPA_RNG(){randomEngine = &gen::getEngineReference();};
00059 private: 
00060   CLHEP::HepRandomEngine* randomEngine;
00061   double Get();
00062 };
00063 
00064 
00065 
00066 SherpaHadronizer::SherpaHadronizer(const edm::ParameterSet &params) :
00067   BaseHadronizer(params),
00068   SherpaPath(params.getUntrackedParameter<std::string>("Path","SherpaRun")),
00069   SherpaPathPiece(params.getUntrackedParameter<std::string>("PathPiece","SherpaRun")),
00070   SherpaResultDir(params.getUntrackedParameter<std::string>("ResultDir","Result")),
00071   default_weight(params.getUntrackedParameter<double>("default_weight",1.)),
00072   SherpaParameter(params.getParameter<edm::ParameterSet>("SherpaParameters")),
00073   maxEventsToPrint(params.getUntrackedParameter<int>("maxEventsToPrint", 0))
00074 {
00075 
00076   // The ids (names) of parameter sets to be read (Analysis,Run) to create Analysis.dat, Run.dat
00077   //They are given as a vstring.  
00078   std::vector<std::string> setNames = SherpaParameter.getParameter<std::vector<std::string> >("parameterSets");
00079   //Loop all set names...
00080   for ( unsigned i=0; i<setNames.size(); ++i ) {
00081     // ...and read the parameters for each set given in vstrings
00082     std::vector<std::string> pars = SherpaParameter.getParameter<std::vector<std::string> >(setNames[i]);
00083     std::cout << "Write Sherpa parameter set " << setNames[i] <<" to "<<setNames[i]<<".dat "<<std::endl;
00084     std::string datfile =  SherpaPath + "/" + setNames[i] +".dat";
00085     std::ofstream os(datfile.c_str());  
00086     // Loop over all strings and write the according *.dat
00087     for(std::vector<std::string>::const_iterator itPar = pars.begin(); itPar != pars.end(); ++itPar ) {
00088       os<<(*itPar)<<std::endl;
00089     } 
00090   }
00091 
00092   //To be conform to the default Sherpa usage create a command line:
00093   //name of executable  (only for demonstration, could also be empty)
00094   std::string shRun  = "./Sherpa";
00095   //Path where the Sherpa libraries are stored
00096   std::string shPath = "PATH=" + SherpaPath;
00097   // new for Sherpa 1.3.0, suggested by authors
00098   std::string shPathPiece = "PATH_PIECE=" + SherpaPathPiece; 
00099   //Path where results are stored 
00100   std::string shRes  = "RESULT_DIRECTORY=" + SherpaResultDir; // from Sherpa 1.2.0 on
00101   //Name of the external random number class
00102   std::string shRng  = "EXTERNAL_RNG=CMS_SHERPA_RNG";
00103   
00104   //create the command line
00105   char* argv[5];
00106   argv[0]=(char*)shRun.c_str();
00107   argv[1]=(char*)shPath.c_str();
00108   argv[2]=(char*)shPathPiece.c_str();
00109   argv[3]=(char*)shRes.c_str();
00110   argv[4]=(char*)shRng.c_str();
00111   
00112   //initialize Sherpa with the command line
00113   Generator.InitializeTheRun(5,argv);
00114 }
00115 
00116 SherpaHadronizer::~SherpaHadronizer()
00117 {
00118 }
00119 
00120 bool SherpaHadronizer::initializeForInternalPartons()
00121 {
00122   
00123   //initialize Sherpa
00124   Generator.InitializeTheEventHandler();
00125   
00126   return true;
00127 }
00128 
00129 #if 0
00130 // naive Sherpa HepMC status fixup //FIXME 
00131 static int getStatus(const HepMC::GenParticle *p)
00132 {
00133   return status;
00134 }
00135 #endif
00136 
00137 //FIXME
00138 bool SherpaHadronizer::declareStableParticles(const std::vector<int> &pdgIds)
00139 {
00140 #if 0
00141   for(std::vector<int>::const_iterator iter = pdgIds.begin();
00142       iter != pdgIds.end(); ++iter)
00143     if (!markStable(*iter))
00144       return false;
00145   
00146   return true;
00147 #else
00148   return false;
00149 #endif
00150 }
00151 
00152 
00153 void SherpaHadronizer::statistics()
00154 {
00155   //calculate statistics
00156   Generator.SummarizeRun(); 
00157 
00158   //get the xsec from EventHandler
00159   SHERPA::Event_Handler* theEventHandler = Generator.GetEventHandler();
00160   double xsec_val = theEventHandler->TotalXS();
00161   double xsec_err = theEventHandler->TotalErr();
00162   
00163   //set the internal cross section in pb in GenRunInfoProduct 
00164   runInfo().setInternalXSec(GenRunInfoProduct::XSec(xsec_val,xsec_err));
00165 
00166 }
00167 
00168 
00169 bool SherpaHadronizer::generatePartonsAndHadronize()
00170 {
00171   //get the next event and check if it produced
00172   if (Generator.GenerateOneEvent()) { 
00173     //convert it to HepMC2
00174     SHERPA::Input_Output_Handler* ioh = Generator.GetIOHandler();
00175     SHERPA::HepMC2_Interface* hm2i = ioh->GetHepMC2Interface();
00176     //get the event weight from blobs
00177     ATOOLS::Blob_List* blobs = Generator.GetEventHandler()-> GetBlobs();
00178     ATOOLS::Blob* sp(blobs->FindFirst(ATOOLS::btp::Signal_Process));
00179     double weight((*sp)["Weight"]->Get<double>());
00180     double ef((*sp)["Enhance"]->Get<double>());
00181     // in case of unweighted events sherpa puts the max weight as event weight. 
00182     // this is not optimal, we want 1 for unweighted events, so we check 
00183     // whether we are producing unweighted events ("EVENT_GENERATION_MODE" == "1")
00184     if ( ATOOLS::ToType<int>( ATOOLS::rpa.gen.Variable("EVENT_GENERATION_MODE") ) == 1 ) {
00185       if (ef > 0.) {
00186         weight = default_weight/ef;
00187       } else {
00188         weight = -1234.;
00189       }
00190     }
00191     //create and empty event and then hand it to SherpaIOHandler to fill it
00192     HepMC::GenEvent* evt = new HepMC::GenEvent();
00193     hm2i->Sherpa2HepMC(blobs, *evt, weight);
00194     resetEvent(evt);         
00195     return true;
00196   }
00197   else {
00198     return false;
00199   }  
00200 }
00201 
00202 bool SherpaHadronizer::decay()
00203 {
00204         return true;
00205 }
00206 
00207 bool SherpaHadronizer::residualDecay()
00208 {
00209         return true;
00210 }
00211 
00212 void SherpaHadronizer::finalizeEvent()
00213 {
00214 #if 0
00215         for(HepMC::GenEvent::particle_iterator iter = event->particles_begin();
00216             iter != event->particles_end(); iter++)
00217                 (*iter)->set_status(getStatus(*iter));
00218 #endif
00219         //******** Verbosity *******
00220         if (maxEventsToPrint > 0) {
00221                 maxEventsToPrint--;
00222                 event()->print();               
00223         }
00224 }
00225 
00226 //GETTER for the external random numbers
00227 DECLARE_GETTER(CMS_SHERPA_RNG_Getter,"CMS_SHERPA_RNG",ATOOLS::External_RNG,ATOOLS::RNG_Key);
00228 
00229 ATOOLS::External_RNG *CMS_SHERPA_RNG_Getter::operator()(const ATOOLS::RNG_Key &) const
00230 { return new CMS_SHERPA_RNG(); }
00231 
00232 void CMS_SHERPA_RNG_Getter::PrintInfo(std::ostream &str,const size_t) const
00233 { str<<"CMS_SHERPA_RNG interface"; }
00234 
00235 double CMS_SHERPA_RNG::Get(){
00236    return randomEngine->flat();
00237    }
00238    
00239 #include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h"
00240 
00241 typedef edm::GeneratorFilter<SherpaHadronizer, gen::ExternalDecayDriver> SherpaGeneratorFilter;
00242 DEFINE_FWK_MODULE(SherpaGeneratorFilter);