CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Pythia6Service.cc
Go to the documentation of this file.
1 #include <algorithm>
2 #include <functional>
3 #include <iostream>
4 #include <sstream>
5 #include <fstream>
6 #include <cmath>
7 #include <string>
8 #include <set>
9 
10 #include <boost/bind.hpp>
11 #include <boost/algorithm/string/classification.hpp>
12 #include <boost/algorithm/string/split.hpp>
13 
14 #include <CLHEP/Random/RandomEngine.h>
15 
17 
20 
22 
25 // #include "GeneratorInterface/Core/interface/ParameterCollector.h"
26 
27 extern "C"
28 {
29  void fioopn_( int* unit, const char* line, int length );
30  void fioopnw_( int* unit, const char* line, int length );
31  void fiocls_( int* unit );
32  void pyslha_( int*, int*, int* );
33  static int call_pyslha(int mupda, int kforig = 0)
34  {
35  int iretrn = 0;
36  pyslha_(&mupda, &kforig, &iretrn);
37  return iretrn;
38  }
39 
40  void pyupda_(int*, int*);
41  static void call_pyupda( int opt, int iunit ){
42  pyupda_( &opt, &iunit );
43  }
44 
45  double gen::pyr_(int *idummy)
46  {
47  // getInstance will throw if no one used enter/leave
48  // or this is the wrong caller class, like e.g. Herwig6Instance
49  Pythia6Service* service = FortranInstance::getInstance<Pythia6Service>();
50  return service->fRandomEngine->flat();
51  }
52 }
53 
54 using namespace gen;
55 using namespace edm;
56 
58 
60  : fRandomEngine(&getEngineReference()), fUnitSLHA(24), fUnitPYUPDA(25)
61 {
62 }
63 
65  : fRandomEngine(&getEngineReference()), fUnitSLHA(24), fUnitPYUPDA(25)
66 {
67  if (fPythia6Owner)
68  throw cms::Exception("PythiaError") <<
69  "Two Pythia6Service instances claiming Pythia6 ownership." <<
70  std::endl;
71 
72  fPythia6Owner = this;
73 
74 /*
75  ParameterCollector collector(ps.getParameter<edm::ParameterSet>("PythiaParameters"));
76 
77  fParamGeneral.clear();
78  fParamCSA.clear();
79  fParamSLHA.clear();
80 
81  fParamGeneral = std::vector<std::string>(collector.begin(), collector.end());
82  fParamCSA = std::vector<std::string>(collector.begin("CSAParameters"), collector.end());
83  fParamSLHA = std::vector<std::string>(collector.begin("SLHAParameters"), collector.end());
84 */
85 
86 
87  // Set PYTHIA parameters in a single ParameterSet
88  //
89  edm::ParameterSet pythia_params =
90  ps.getParameter<edm::ParameterSet>("PythiaParameters") ;
91 
92  // read and sort Pythia6 cards
93  //
94  std::vector<std::string> setNames =
95  pythia_params.getParameter<std::vector<std::string> >("parameterSets");
96 
97  // std::vector<std::string> paramLines;
98  fParamGeneral.clear();
99  fParamCSA.clear();
100  fParamSLHA.clear();
101  fParamPYUPDA.clear();
102 
103 
104  for(std::vector<std::string>::const_iterator iter = setNames.begin();
105  iter != setNames.end(); ++iter)
106  {
107  std::vector<std::string> lines =
108  pythia_params.getParameter< std::vector<std::string> >(*iter);
109 
110  for(std::vector<std::string>::const_iterator line = lines.begin();
111  line != lines.end(); ++line )
112  {
113  if (line->substr(0, 7) == "MRPY(1)")
114  throw cms::Exception("PythiaError") <<
115  "Attempted to set random number"
116  " using Pythia command 'MRPY(1)'."
117  " Please use the"
118  " RandomNumberGeneratorService." <<
119  std::endl;
120 
121  if ( *iter == "CSAParameters" )
122  {
123  fParamCSA.push_back(*line);
124  }
125  else if ( *iter == "SLHAParameters" )
126  {
127  fParamSLHA.push_back(*line);
128  }
129  else if ( *iter == "PYUPDAParameters" )
130  {
131  fParamPYUPDA.push_back(*line);
132  }
133  else
134  {
135  fParamGeneral.push_back(*line);
136  }
137  }
138  }
139 }
140 
142 {
143  if (fPythia6Owner == this)
144  fPythia6Owner = 0;
145 
146  fParamGeneral.clear();
147  fParamCSA.clear();
148  fParamSLHA.clear();
149  fParamPYUPDA.clear();
150 }
151 
153 {
155 
156  if (!fPythia6Owner) {
157  edm::LogInfo("Generator|Pythia6Interface") <<
158  "gen::Pythia6Service is going to initialise Pythia, as no other "
159  "instace has done so yet, and Pythia service routines have been "
160  "requested by a dummy instance." << std::endl;
161 
162  call_pygive("MSTU(12)=12345");
163  call_pyinit("NONE", "", "", 0.0);
164 
165  fPythia6Owner = this;
166  }
167 }
168 
170 {
171  // now pass general config cards
172  //
173  for(std::vector<std::string>::const_iterator iter = fParamGeneral.begin();
174  iter != fParamGeneral.end(); ++iter)
175  {
176  if (!call_pygive(*iter))
177  throw cms::Exception("PythiaError")
178  << "Pythia did not accept \""
179  << *iter << "\"." << std::endl;
180  }
181 
182  return ;
183 }
184 
186 {
187 
188  txgive_init_();
189 
190  for(std::vector<std::string>::const_iterator iter = fParamCSA.begin();
191  iter != fParamCSA.end(); ++iter)
192  {
193  txgive_( iter->c_str(), iter->length() );
194  }
195 
196  return ;
197 }
198 
199 void Pythia6Service::openSLHA( const char* file )
200 {
201 
202  std::ostringstream pyCard1 ;
203  pyCard1 << "IMSS(21)=" << fUnitSLHA;
204  call_pygive( pyCard1.str() );
205  std::ostringstream pyCard2 ;
206  pyCard2 << "IMSS(22)=" << fUnitSLHA;
207  call_pygive( pyCard2.str() );
208 
209  fioopn_( &fUnitSLHA, file, strlen(file) );
210 
211  return;
212 
213 }
214 
215 void Pythia6Service::openPYUPDA( const char* file, bool write_file )
216 {
217 
218  if (write_file) {
219  std::cout<<"=== WRITING PYUPDA FILE "<<file<<" ==="<<std::endl;
220  fioopnw_( &fUnitPYUPDA, file, strlen(file) );
221  // Write Pythia particle table to this card file.
223  } else {
224  std::cout<<"=== READING PYUPDA FILE "<<file<<" ==="<<std::endl;
225  fioopn_( &fUnitPYUPDA, file, strlen(file) );
226  // Update Pythia particle table with this card file.
228  }
229 
230  return;
231 
232 }
233 
235 {
236 
237  fiocls_(&fUnitSLHA);
238 
239  return;
240 }
241 
243 {
244 
246 
247  return;
248 
249 }
250 
252 {
253  for (std::vector<std::string>::const_iterator iter = fParamSLHA.begin();
254  iter != fParamSLHA.end(); iter++ )
255  {
256 
257  if( iter->find( "SLHAFILE", 0 ) == std::string::npos ) continue;
258  std::string::size_type start = iter->find_first_of( "=" ) + 1;
259  std::string::size_type end = iter->length() - 1;
260  std::string::size_type temp = iter->find_first_of( "'", start );
261  if( temp != std::string::npos ) {
262  start = temp + 1;
263  end = iter->find_last_of( "'" ) - 1;
264  }
265  start = iter->find_first_not_of( " ", start );
266  end = iter->find_last_not_of( " ", end );
267  //std::cout << " start, end = " << start << " " << end << std::endl;
268  std::string shortfile = iter->substr( start, end - start + 1 );
269  FileInPath f1( shortfile );
270  std::string file = f1.fullPath();
271 
272 /*
273  //
274  // override what might have be given via the external config
275  //
276  std::ostringstream pyCard ;
277  pyCard << "IMSS(21)=" << fUnitSLHA;
278  call_pygive( pyCard.str() );
279  pyCard << "IMSS(22)=" << fUnitSLHA;
280  call_pygive( pyCard.str() );
281 
282  fioopn_( &fUnitSLHA, file.c_str(), file.length() );
283 */
284 
285  openSLHA( file.c_str() );
286 
287  }
288 
289  return;
290 }
291 
292 void Pythia6Service::setPYUPDAParams(bool afterPyinit)
293 {
294  std::string shortfile;
295  bool write_file = false;
296  bool usePostPyinit = false;
297 
298  // std::cout<<"=== CALLING setPYUPDAParams === "<<afterPyinit<<" "<<fParamPYUPDA.size()<<std::endl;
299 
300  // This assumes that PYUPDAFILE only appears once ...
301 
302  for (std::vector<std::string>::const_iterator iter = fParamPYUPDA.begin();
303  iter != fParamPYUPDA.end(); iter++ )
304  {
305  // std::cout<<"PYUPDA check "<<*iter<<std::endl;
306  if( iter->find( "PYUPDAFILE", 0 ) != std::string::npos ) {
307  std::string::size_type start = iter->find_first_of( "=" ) + 1;
308  std::string::size_type end = iter->length() - 1;
309  std::string::size_type temp = iter->find_first_of( "'", start );
310  if( temp != std::string::npos ) {
311  start = temp + 1;
312  end = iter->find_last_of( "'" ) - 1;
313  }
314  start = iter->find_first_not_of( " ", start );
315  end = iter->find_last_not_of( " ", end );
316  //std::cout << " start, end = " << start << " " << end << std::endl;
317  shortfile = iter->substr( start, end - start + 1 );
318  } else if ( iter->find( "PYUPDAWRITE", 0 ) != std::string::npos ) {
319  write_file = true;
320  } else if ( iter->find( "PYUPDApostPYINIT", 0 ) != std::string::npos ) {
321  usePostPyinit = true;
322  }
323  }
324 
325  if (!shortfile.empty()) {
326  std::string file;
327  if (write_file) {
328  file = shortfile;
329  } else {
330  // If reading, get full path to file and require it to exist.
331  FileInPath f1( shortfile );
332  file = f1.fullPath();
333  }
334 
335  if (afterPyinit == usePostPyinit || (write_file && afterPyinit)) {
336  openPYUPDA( file.c_str(), write_file );
337  }
338  }
339 
340  return;
341 }
342 
343 void Pythia6Service::setSLHAFromHeader( const std::vector<std::string> &lines )
344 {
345 
346  std::set<std::string> blocks;
347  unsigned int model = 0, subModel = 0;
348 
349  const char *fname = std::tmpnam(NULL);
350  std::ofstream file(fname, std::fstream::out | std::fstream::trunc);
351  std::string block;
352  for(std::vector<std::string>::const_iterator iter = lines.begin();
353  iter != lines.end(); ++iter) {
354  file << *iter;
355 
356  std::string line = *iter;
357  std::transform(line.begin(), line.end(),
358  line.begin(), (int(*)(int))std::toupper);
359  std::string::size_type pos = line.find('#');
360  if (pos != std::string::npos)
361  line.resize(pos);
362 
363  if (line.empty())
364  continue;
365 
366  if (!boost::algorithm::is_space()(line[0])) {
367  std::vector<std::string> tokens;
368  boost::split(tokens, line,
369  boost::algorithm::is_space(),
370  boost::token_compress_on);
371  if (!tokens.size())
372  continue;
373  block.clear();
374  if (tokens.size() < 2)
375  continue;
376  if (tokens[0] == "BLOCK") {
377  block = tokens[1];
378  blocks.insert(block);
379  continue;
380  }
381 
382  if (tokens[0] == "DECAY") {
383  block = "DECAY";
384  blocks.insert(block);
385  }
386  } else if (block == "MODSEL") {
387  std::istringstream ss(line);
388  ss >> model >> subModel;
389  } else if (block == "SMINPUTS") {
390  std::istringstream ss(line);
391  int index;
392  double value;
393  ss >> index >> value;
394  switch(index) {
395  case 1:
396  pydat1_.paru[103 - 1] = 1.0 / value;
397  break;
398  case 2:
399  pydat1_.paru[105 - 1] = value;
400  break;
401  case 4:
402  pydat2_.pmas[0][23 - 1] = value;
403  break;
404  case 6:
405  pydat2_.pmas[0][6 - 1] = value;
406  break;
407  case 7:
408  pydat2_.pmas[0][15 - 1] = value;
409  break;
410  }
411  }
412  }
413  file.close();
414 
415  if (blocks.count("SMINPUTS"))
416  pydat1_.paru[102 - 1] = 0.5 - std::sqrt(0.25 -
417  pydat1_.paru[0] * M_SQRT1_2 *
418  pydat1_.paru[103 - 1] / pydat1_.paru[105 - 1] /
419  (pydat2_.pmas[0][23 - 1] * pydat2_.pmas[0][23 - 1]));
420 
421 /*
422  int unit = 24;
423  fioopn_(&unit, fname, std::strlen(fname));
424  std::remove(fname);
425 
426  call_pygive("IMSS(21)=24");
427  call_pygive("IMSS(22)=24");
428 */
429 
430  openSLHA( fname ) ;
431  std::remove( fname );
432 
433  if (model ||
434  blocks.count("HIGMIX") ||
435  blocks.count("SBOTMIX") ||
436  blocks.count("STOPMIX") ||
437  blocks.count("STAUMIX") ||
438  blocks.count("AMIX") ||
439  blocks.count("NMIX") ||
440  blocks.count("UMIX") ||
441  blocks.count("VMIX"))
442  call_pyslha(1);
443  if (model ||
444  blocks.count("QNUMBERS") ||
445  blocks.count("PARTICLE") ||
446  blocks.count("MINPAR") ||
447  blocks.count("EXTPAR") ||
448  blocks.count("SMINPUTS") ||
449  blocks.count("SMINPUTS"))
450  call_pyslha(0);
451  if (blocks.count("MASS"))
452  call_pyslha(5, 0);
453  if (blocks.count("DECAY"))
454  call_pyslha(2);
455 
456  return ;
457 
458 }
T getParameter(std::string const &) const
static void call_pyupda(int opt, int iunit)
void pyupda_(int *, int *)
void fioopn_(int *unit, const char *line, int length)
list file
Definition: dbtoweb.py:253
void fiocls_(int *unit)
bool call_pygive(const std::string &line)
#define NULL
Definition: scimark2.h:8
uint16_t size_type
return((rh^lh)&mask)
void fioopnw_(int *unit, const char *line, int length)
CLHEP::HepRandomEngine & getEngineReference()
void pyslha_(int *, int *, int *)
string unit
Definition: csvLumiCalc.py:46
T sqrt(T t)
Definition: SSEVec.h:28
void setSLHAFromHeader(const std::vector< std::string > &lines)
std::vector< std::string > fParamPYUPDA
static int call_pyslha(int mupda, int kforig=0)
virtual void enter()
virtual void enter()
void setPYUPDAParams(bool afterPyinit)
#define end
Definition: vmac.h:38
std::vector< std::string > fParamGeneral
block
Formating index page&#39;s pieces.
Definition: Association.py:187
static Pythia6Service * fPythia6Owner
tuple out
Definition: dbtoconf.py:99
struct gen::@242 pydat1_
void txgive_init_(void)
CLHEP::HepRandomEngine * fRandomEngine
std::vector< std::string > fParamCSA
void openSLHA(const char *)
string fname
main script
void openPYUPDA(const char *, bool write_file)
tuple cout
Definition: gather_cfg.py:41
std::vector< std::string > fParamSLHA
std::string fullPath() const
Definition: FileInPath.cc:170
double pyr_(int *idummy)
void txgive_(const char *, int)
double split
Definition: MVATrainer.cc:139