CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SherpaHadronizer.cc
Go to the documentation of this file.
1 #include <cassert>
2 #include <iostream>
3 #include <sstream>
4 #include <string>
5 #include <memory>
6 #include <cstdint>
7 #include <vector>
8 
9 
10 #include "SHERPA/Main/Sherpa.H"
11 #include "ATOOLS/Math/Random.H"
12 
13 #include "ATOOLS/Org/Run_Parameter.H"
14 #include "ATOOLS/Org/MyStrStream.H"
15 #include "ATOOLS/Org/CXXFLAGS.H"
16 #include "ATOOLS/Org/CXXFLAGS_PACKAGES.H"
17 #include "ATOOLS/Org/My_MPI.H"
18 
24 
25 #include "CLHEP/Random/RandomEngine.h"
26 
27 
28 //This unnamed namespace is used (instead of static variables) to pass the
29 //randomEngine passed to doSetRandomEngine to the External Random
30 //Number Generator CMS_SHERPA_RNG of sherpa
31 //The advantage of the unnamed namespace over static variables is
32 //that it is only accessible in this file
33 
34 namespace {
35  CLHEP::HepRandomEngine* ExternalEngine=nullptr;
36  CLHEP::HepRandomEngine* GetExternalEngine() { return ExternalEngine; }
37  void SetExternalEngine(CLHEP::HepRandomEngine* v) { ExternalEngine=v; }
38 }
39 
41 public:
42  SherpaHadronizer(const edm::ParameterSet &params);
43  ~SherpaHadronizer() override;
44 
45  bool readSettings( int ) { return true; }
47  bool declareStableParticles(const std::vector<int> &pdgIds);
48  bool declareSpecialSettings( const std::vector<std::string>& ) { return true; }
49  void statistics();
51  bool decay();
53  bool residualDecay();
54  void finalizeEvent();
55  GenLumiInfoHeader *getGenLumiInfoHeader() const override;
56  const char *classname() const { return "SherpaHadronizer"; }
57 
58 
59 private:
60 
61  void doSetRandomEngine(CLHEP::HepRandomEngine* v) override;
62 
70  unsigned int maxEventsToPrint;
71  std::vector<std::string> arguments;
72  SHERPA::Sherpa *Generator = new SHERPA::Sherpa();
75  std::vector<std::string> weightlist;
76  std::vector<std::string> variationweightlist;
77 };
78 
79 
80 
81 
82 class CMS_SHERPA_RNG: public ATOOLS::External_RNG {
83 public:
84 
86  edm::LogVerbatim("SherpaHadronizer") << "Use stored reference for the external RNG";
87  setRandomEngine(GetExternalEngine());
88  }
89  void setRandomEngine(CLHEP::HepRandomEngine* v) { randomEngine = v; }
90 
91 private:
92  double Get() override;
93  CLHEP::HepRandomEngine* randomEngine;
94 };
95 
96 
97 void SherpaHadronizer::doSetRandomEngine(CLHEP::HepRandomEngine* v) {
98  CMS_SHERPA_RNG* cmsSherpaRng = dynamic_cast<CMS_SHERPA_RNG*>(ATOOLS::ran->GetExternalRng());
99  if (cmsSherpaRng ==nullptr) {
100  //First time call to this function makes the interface store the reference in the unnamed namespace
101  if (!isRNGinitialized){
102  isRNGinitialized=true;
103  edm::LogVerbatim("SherpaHadronizer") << "Store assigned reference of the randomEngine";
104  SetExternalEngine(v);
105  // Throw exception if there is no reference to an external RNG and it is not the first call!
106  } else {
108  << "The Sherpa interface got a randomEngine reference but there is no reference to the external RNG to hand it over to\n";
109  }
110  } else {
111  cmsSherpaRng->setRandomEngine(v);
112  }
113 }
114 
116  BaseHadronizer(params),
117  SherpaParameterSet(params.getParameter<edm::ParameterSet>("SherpaParameters")),
119 {
120  if (!params.exists("SherpaProcess")) SherpaProcess="";
121  else SherpaProcess=params.getParameter<std::string>("SherpaProcess");
122  if (!params.exists("SherpaPath")) SherpaPath="";
123  else SherpaPath=params.getParameter<std::string>("SherpaPath");
124  if (!params.exists("SherpaPathPiece")) SherpaPathPiece="";
125  else SherpaPathPiece=params.getParameter<std::string>("SherpaPathPiece");
126  if (!params.exists("SherpaResultDir")) SherpaResultDir="Result";
127  else SherpaResultDir=params.getParameter<std::string>("SherpaResultDir");
128  if (!params.exists("SherpaDefaultWeight")) SherpaDefaultWeight=1.;
129  else SherpaDefaultWeight=params.getParameter<double>("SherpaDefaultWeight");
130  if (!params.exists("maxEventsToPrint")) maxEventsToPrint=0;
131  else maxEventsToPrint=params.getParameter<int>("maxEventsToPrint");
132 // if hepmcextendedweights is used the event weights have to be reordered ( unordered list can be accessed via event->weights().write() )
133 // two lists have to be provided:
134 // 1) SherpaWeights
135 // - containing nominal event weight, combined matrix element and phase space weight, event normalization, and possibly other sherpa weights
136 // 2) SherpaVariationsWeights
137 // - containing weights from scale and PDF variations ( have to be defined in the runcard )
138 // - in case of unweighted events these weights are also divided by the event normalization (see list 1 )
139 // Sherpa Documentation: http://sherpa.hepforge.org/doc/SHERPA-MC-2.2.0.html#Scale-and-PDF-variations
140  if (!params.exists("SherpaWeightsBlock")) {
141  rearrangeWeights=false;
142  } else {
143  rearrangeWeights=true;
144  edm::ParameterSet WeightsBlock = params.getParameter<edm::ParameterSet>("SherpaWeightsBlock");
145  if (WeightsBlock.exists("SherpaWeights"))
146  weightlist=WeightsBlock.getParameter< std::vector<std::string> >("SherpaWeights");
147  else
148  throw cms::Exception("SherpaInterface") <<"SherpaWeights does not exists in SherpaWeightsBlock" << std::endl;
149  if (WeightsBlock.exists("SherpaVariationWeights"))
150  variationweightlist=WeightsBlock.getParameter< std::vector<std::string> >("SherpaVariationWeights");
151  else
152  throw cms::Exception("SherpaInterface") <<"SherpaVariationWeights does not exists in SherpaWeightsBlock" << std::endl;
153  edm::LogVerbatim("SherpaHadronizer") << "SherpaHadronizer will try rearrange the event weights according to SherpaWeights and SherpaVariationWeights";
154  }
155 
156 
157  spf::SherpackFetcher Fetcher(params);
158  int retval=Fetcher.Fetch();
159  if (retval != 0) {
160  throw cms::Exception("SherpaInterface") <<"SherpaHadronizer: Preparation of Sherpack failed ... ";
161  }
162  // The ids (names) of parameter sets to be read (Analysis,Run) to create Analysis.dat, Run.dat
163  //They are given as a vstring.
164  std::vector<std::string> setNames = SherpaParameterSet.getParameter<std::vector<std::string> >("parameterSets");
165  //Loop all set names...
166  for ( unsigned i=0; i<setNames.size(); ++i ) {
167  // ...and read the parameters for each set given in vstrings
168  std::vector<std::string> pars = SherpaParameterSet.getParameter<std::vector<std::string> >(setNames[i]);
169  edm::LogVerbatim("SherpaHadronizer") << "Write Sherpa parameter set " << setNames[i] <<" to "<<setNames[i]<<".dat ";
170  std::string datfile = SherpaPath + "/" + setNames[i] +".dat";
171  std::ofstream os(datfile.c_str());
172  // Loop over all strings and write the according *.dat
173  for(std::vector<std::string>::const_iterator itPar = pars.begin(); itPar != pars.end(); ++itPar ) {
174  os<<(*itPar)<<std::endl;
175  }
176  }
177 
178  //To be conform to the default Sherpa usage create a command line:
179  //name of executable (only for demonstration, could also be empty)
180  std::string shRun = "./Sherpa";
181  //Path where the Sherpa libraries are stored
182  std::string shPath = "PATH=" + SherpaPath;
183  // new for Sherpa 1.3.0, suggested by authors
184  std::string shPathPiece = "PATH_PIECE=" + SherpaPathPiece;
185  //Path where results are stored
186  std::string shRes = "RESULT_DIRECTORY=" + SherpaResultDir; // from Sherpa 1.2.0 on
187  //Name of the external random number class
188  std::string shRng = "EXTERNAL_RNG=CMS_SHERPA_RNG";
189 
190  //create the command line
191  arguments.push_back(shRun);
192  arguments.push_back(shPath);
193  arguments.push_back(shPathPiece);
194  arguments.push_back(shRes);
195  arguments.push_back(shRng);
196  isInitialized=false;
197  //initialization of Sherpa moved to initializeForInternalPartons
198 }
199 
201 {
202  Generator->~Sherpa();
203  #ifdef USING__MPI
204  MPI::Finalize();
205  #endif
206 }
207 
209 {
210  //initialize Sherpa but only once
211  if (!isInitialized){
212  int argc=arguments.size();
213  char* argv[argc];
214  for (int l=0; l<argc; l++) argv[l]=(char*)arguments[l].c_str();
215  #ifdef USING__MPI
216  MPI::Init();
217  #endif
218  Generator->InitializeTheRun(argc,argv);
219  Generator->InitializeTheEventHandler();
220  isInitialized=true;
221  }
222  return true;
223 }
224 
225 
226 bool SherpaHadronizer::declareStableParticles(const std::vector<int> &pdgIds)
227 {
228  return false;
229 }
230 
231 
233 {
234  //calculate statistics
235  Generator->SummarizeRun();
236 
237  //get the xsec & err
238  double xsec_val = Generator->TotalXS();
239  double xsec_err = Generator->TotalErr();
240 
241  //set the internal cross section in pb in GenRunInfoProduct
242  runInfo().setInternalXSec(GenRunInfoProduct::XSec(xsec_val,xsec_err));
243 }
244 
245 
247 {
248  //get the next event and check if it produced
249  bool rc = false;
250  int itry = 0;
251  bool gen_event = true;
252  while((itry < 3) && gen_event){
253  try{
254  rc = Generator->GenerateOneEvent();
255  gen_event = false;
256  } catch(...){
257  ++itry;
258  std::cerr << "Exception from Generator->GenerateOneEvent() catch. Call # "
259  << itry << " for this event\n";
260  }
261  }
262  if (rc) {
263  //convert it to HepMC2
264  HepMC::GenEvent* evt = new HepMC::GenEvent();
265  Generator->FillHepMCEvent(*evt);
266 
267  // in case of unweighted events sherpa puts the max weight as event weight.
268  // this is not optimal, we want 1 for unweighted events, so we check
269  // whether we are producing unweighted events ("EVENT_GENERATION_MODE" == "1")
270  // the information about the weights to the HepMC weight vector:
271  // [0] event weight
272  // [1] combined matrix element and phase space weight (missing only PDF information, thus directly suitable for PDF reweighting)
273  // [2] event weight normalisation (in case of unweighted events event weights of ~ +/-1 can be obtained by (event weight)/(event weight normalisation))
274  // [3] number of trials.
275  // see also: https://sherpa.hepforge.org/doc/SHERPA-MC-2.1.0.html#Event-output-formats
276  bool unweighted = false;
277  double weight_normalization = -1;
278  if(ATOOLS::ToType<int>(ATOOLS::rpa->gen.Variable("EVENT_GENERATION_MODE")) == 1){
279  if (evt->weights().size()>2) {
280  unweighted = true;
281  weight_normalization = evt->weights()[2];
282  }else{
283  throw cms::Exception("SherpaInterface") <<"Requested unweighted production. Missing normalization weight." << std::endl;
284  }
285  }
286 
287  // vector to fill new weights in correct order
288  std::vector<double> newWeights;
289  if (rearrangeWeights){
290  for ( auto &i : weightlist ) {
291  if (evt->weights().has_key(i)) {
292  newWeights.push_back(evt->weights()[i]);
293  } else {
294  throw cms::Exception("SherpaInterface") <<"Missing weights! Key " << i << " not found, please check the weight definition!" << std::endl;
295  }
296  }
297  for ( auto &i : variationweightlist ) {
298  if (evt->weights().has_key(i)) {
299  if(unweighted){
300  newWeights.push_back(evt->weights()[i]/weight_normalization);
301  }else{
302  newWeights.push_back(evt->weights()[i]);
303  }
304  } else {
305  throw cms::Exception("SherpaInterface") <<"Missing weights! Key " << i << " not found, please check the weight definition!" << std::endl;
306  }
307 
308  }
309 
310  //Change original weights for reordered ones
311  evt->weights().clear();
312  for (auto& elem: newWeights) {
313  evt->weights().push_back(elem);
314  }
315  }
316 
317  if(unweighted){
318  evt->weights()[0]/=weight_normalization;
319  }
320  resetEvent(evt);
321  return true;
322  }
323  else {
324  return false;
325  }
326 }
327 
329 {
330  return true;
331 }
332 
334 {
335  return true;
336 }
337 
339 {
340  //******** Verbosity *******
341  if (maxEventsToPrint > 0) {
343  event()->print();
344  }
345 }
346 
347 
348 //GETTER for the external random numbers
349 DECLARE_GETTER(CMS_SHERPA_RNG,"CMS_SHERPA_RNG",ATOOLS::External_RNG,ATOOLS::RNG_Key);
350 
351 ATOOLS::External_RNG *ATOOLS::Getter<ATOOLS::External_RNG,ATOOLS::RNG_Key,CMS_SHERPA_RNG>::operator()(const ATOOLS::RNG_Key &) const
352 { return new CMS_SHERPA_RNG(); }
353 
354 void ATOOLS::Getter<ATOOLS::External_RNG,ATOOLS::RNG_Key,CMS_SHERPA_RNG>::PrintInfo(std::ostream &str,const size_t) const
355 { str<<"CMS_SHERPA_RNG interface"; }
356 
358  if(randomEngine == nullptr) {
360  << "The Sherpa code attempted to a generate random number while\n"
361  << "the engine pointer was null. This might mean that the code\n"
362  << "was modified to generate a random number outside the event and\n"
363  << "beginLuminosityBlock methods, which is not allowed.\n";
364  }
365  return randomEngine->flat();
366 
367 }
369  GenLumiInfoHeader *genLumiInfoHeader = BaseHadronizer::getGenLumiInfoHeader();
370 
371  if(rearrangeWeights){
372  edm::LogPrint("SherpaHadronizer") << "The order of event weights was changed!" ;
373  for(auto &i: weightlist){
374  genLumiInfoHeader->weightNames().push_back(i);
375  edm::LogVerbatim("SherpaHadronizer") << i;
376  }
377  for(auto &i: variationweightlist) {
378  genLumiInfoHeader->weightNames().push_back(i);
379  edm::LogVerbatim("SherpaHadronizer") << i;
380  }
381  }
382  return genLumiInfoHeader;
383 }
385 
CLHEP::HepRandomEngine * randomEngine
Definition: Dummies.cc:7
T getParameter(std::string const &) const
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
CLHEP::HepRandomEngine * randomEngine
BaseHadronizer(edm::ParameterSet const &ps)
const std::vector< std::string > & weightNames() const
bool exists(std::string const &parameterName) const
checks if a parameter exists
std::vector< std::string > weightlist
std::string SherpaPathPiece
unsigned int maxEventsToPrint
edm::GeneratorFilter< SherpaHadronizer, gen::ExternalDecayDriver > SherpaGeneratorFilter
#define nullptr
void setInternalXSec(const XSec &xsec)
bool declareStableParticles(const std::vector< int > &pdgIds)
std::auto_ptr< HepMC::GenEvent > & event()
GenRunInfoProduct & runInfo()
edm::ParameterSet SherpaParameterSet
bool initializeForInternalPartons()
void setRandomEngine(CLHEP::HepRandomEngine *v)
std::vector< std::string > variationweightlist
bool generatePartonsAndHadronize()
std::string SherpaPath
def elem(elemtype, innerHTML='', html_class='', kwargs)
Definition: HTMLExport.py:18
~SherpaHadronizer() override
bool declareSpecialSettings(const std::vector< std::string > &)
void setRandomEngine(CLHEP::HepRandomEngine *v)
void doSetRandomEngine(CLHEP::HepRandomEngine *v) override
std::string SherpaProcess
GenLumiInfoHeader * getGenLumiInfoHeader() const override
std::string SherpaResultDir
HLT enums.
std::vector< std::string > arguments
SherpaHadronizer(const edm::ParameterSet &params)
void resetEvent(HepMC::GenEvent *event)
DECLARE_GETTER(CMS_SHERPA_RNG,"CMS_SHERPA_RNG", ATOOLS::External_RNG, ATOOLS::RNG_Key)
const char * classname() const
std::string SherpaChecksum
double Get() override