ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHPythia6.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHPythia6.cc
1 #include "PHPythia6.h"
2 #include "PHPy6GenTrigger.h"
3 
4 #include <phhepmc/PHHepMCGenHelper.h> // for PHHepMCGenHelper
5 
7 #include <fun4all/SubsysReco.h> // for SubsysReco
8 
9 #include <phool/PHRandomSeed.h>
10 #include <phool/phool.h> // for PHWHERE
11 
12 #include <Rtypes.h> // for Int_t Float_t
13 
14 #include <HepMC/GenEvent.h>
15 #include <HepMC/HEPEVT_Wrapper.h> // for HEPEVT_Wrappe
16 #include <HepMC/IO_BaseClass.h> // for IO_BaseClass
17 #include <HepMC/IO_GenEvent.h>
18 #include <HepMC/IO_HEPEVT.h>
19 #include <HepMC/PdfInfo.h> // for PdfInfo
20 #include <HepMC/PythiaWrapper.h>
21 #include <HepMC/PythiaWrapper6_4.h> // for (anonymous), pypars, pydat1
22 #include <HepMC/Units.h> // for GEV, MM
23 
24 #include <algorithm> // for transform
25 #include <cctype> // for tolower
26 #include <cmath> // for fmod
27 #include <cstdlib> // for exit, abs
28 #include <iostream> // for operator<<, endl, basic_ostream
29 #include <sstream>
30 
31 class PHHepMCGenEvent;
32 
33 #define pytune pytune_
34 extern "C" int pytune(int *itune);
35 
36 using namespace std;
37 
38 PHPythia6::PHPythia6(const std::string &name)
39  : SubsysReco(name)
40  , _eventcount(0)
41  , _geneventcount(0)
42  , _configFile("phpythia6.cfg")
43  , _save_ascii(false)
44  , _filename_ascii("pythia_hepmc.dat")
45  , _registeredTriggers()
46  , _triggersOR(true)
47  , _triggersAND(false)
48 {
49  PHHepMCGenHelper::set_embedding_id(1); // default embedding ID to 1
50 }
51 
53 {
54  /* Create node tree */
55  create_node_tree(topNode);
56 
57  /* event numbering will start from 1 */
58  _eventcount = 0;
59 
60  /* HepMC/example_MyPythia.cc:
61  *
62  * Pythia 6.1 uses HEPEVT with 4000 entries and 8-byte floating point
63  * numbers. We need to explicitly pass this information to the
64  * HEPEVT_Wrapper.
65  */
66  HepMC::HEPEVT_Wrapper::set_max_number_entries(4000);
67  HepMC::HEPEVT_Wrapper::set_sizeof_real(8);
68 
69  /* set pythia random number seed (mandatory!) */
70  int fSeed = PHRandomSeed();
71  fSeed = abs(fSeed) % 900000000;
72 
73  if ((fSeed >= 0) && (fSeed <= 900000000))
74  {
75  pydatr.mrpy[0] = fSeed; // set seed
76  pydatr.mrpy[1] = 0; // use new seed
77  }
78  else
79  {
80  cout << PHWHERE << " ERROR: seed " << fSeed << " is not valid" << endl;
81  exit(2);
82  }
83  // print out seed so we can make this is reproducible
84  if (Verbosity()) cout << "PHPythia6 random seed: " << fSeed << endl;
85 
86  /* read pythia configuration and initialize */
87  if (!_configFile.empty()) ReadConfig(_configFile);
88 
89  /* Initialization done */
91 }
92 
94 {
95  //........................................TERMINATION
96  // write out some information from Pythia to the screen
97  call_pystat(1);
98 
99  //match pythia printout
100  cout << " | "
101  << " | " << endl;
102  cout << " PHPythia6::End - " << _eventcount
103  << " events passed trigger" << endl;
104  cout << " Fraction passed: " << _eventcount
105  << "/" << _geneventcount
106  << " = " << _eventcount / float(_geneventcount) << endl;
107  cout << " *------- End PYTHIA Trigger Statistics ------------------------"
108  << "-------------------------------------------------* " << endl;
109 
111 }
112 
113 //__________________________________________________________
114 int PHPythia6::ReadConfig(const string &cfg_file)
115 {
116  if (cfg_file != "") _configFile = cfg_file;
117  cout << "PHPythia6::read_config - Reading " << _configFile << endl;
118 
119  ifstream infile(_configFile);
120  if (infile.fail())
121  {
122  cout << "PHPythia6::read_config - Failed to open file " << _configFile << endl;
123  exit(2);
124  }
125 
126  // initialize variables
127  Int_t _nevents(0);
128  Float_t _roots(0);
129  string _proj;
130  string _targ;
131  string _frame;
132 
133  string FullLine; // a complete line in the config file
134  string label; // the label
135 
136  int index = 999999;
137  double value = 1e9;
138 
139  // get one line first
140  getline(infile, FullLine);
141  while (!infile.eof())
142  {
143  // skip lines that begin with #, or "//"
144  if (FullLine[0] == '#' || FullLine.substr(0, 2) == "//")
145  {
146  getline(infile, FullLine);
147  continue;
148  }
149 
150  // make FullLine an istringstream
151  istringstream line(FullLine.c_str());
152 
153  // get label
154  line >> label;
155 
156  // to lower case
157  std::transform(label.begin(), label.end(), label.begin(), (int (*)(int)) std::tolower);
158 
159  // based on label, fill correct item
160  if (label == "nevents")
161  {
162  line >> _nevents;
163  cout << "nevents\t" << _nevents << endl;
164  }
165  else if (label == "roots")
166  {
167  line >> _roots;
168  cout << "roots\t" << _roots << endl;
169  }
170  else if (label == "proj")
171  {
172  line >> _proj;
173  cout << "proj\t" << _proj << endl;
174  }
175  else if (label == "targ")
176  {
177  line >> _targ;
178  cout << "targ\t" << _targ << endl;
179  }
180  else if (label == "frame")
181  {
182  line >> _frame;
183  cout << "frame\t" << _frame << endl;
184  }
185  else if (label == "p1" || label == "p2")
186  {
187  if (label == "p1") //Momentum of Projectile Beam (e- for e-p)
188  {
189  line >> index >> value;
190  //Index Options(3MOM): 1 = x-momentum; 2 = y-momentum; 3 = z-momentum
191  pyjets.p[index - 1][0] = value;
192  cout << "p1\t" << index << " " << value << endl;
193  }
194  if (label == "p2") //Momentum of Target Beam (p for e-p)
195  {
196  line >> index >> value;
197  //Index Options(3MOM): 1 = x-momentum; 2 = y-momentum; 3 = z-momentum
198  pyjets.p[index - 1][1] = value;
199  cout << "p2\t" << index << " " << value << endl;
200  }
201  /*
202  int entry = 0;
203  if ( label=="p1") entry = 1;
204  else if ( label=="p2") entry = 2;
205 
206  char temp_line[10000];
207  strcpy(temp_line,FullLine.c_str());
208  char *sptr = strtok(temp_line," \t"); // skip first word
209  sptr = strtok(NULL," \t");
210  cout << label;
211  int index = 1;
212  while ( sptr != NULL )
213  {
214  double val = atof(sptr);
215  cout << "Entry: " << entry << endl;
216  index++;
217  cout << "\t" << val;
218  sptr = strtok(NULL," \t");
219  }
220  cout << endl;*/
221  }
222  else if (label == "msel")
223  {
224  line >> value;
225  pysubs.msel = (int) value;
226  cout << "msel\t" << value << endl;
227  IntegerTest(value);
228  }
229  else if (label == "msub")
230  {
231  line >> index >> value;
232  // careful with C/F77 differences: arrays in C start at 0, F77 at 1,
233  // so we need to subtract 1 from the process #)
234  pysubs.msub[index - 1] = (int) value;
235  cout << "msub\t" << index << " " << value << endl;
236  IntegerTest(value);
237  }
238  else if (label == "mstp")
239  {
240  line >> index >> value;
241  pypars.mstp[index - 1] = (int) value;
242  cout << "mstp\t" << index << " " << value << endl;
243  IntegerTest(value);
244  }
245  else if (label == "mstj")
246  {
247  line >> index >> value;
248  pydat1.mstj[index - 1] = (int) value;
249  cout << "mstj\t" << index << " " << value << endl;
250  IntegerTest(value);
251  }
252  else if (label == "mstu")
253  {
254  line >> index >> value;
255  pydat1.mstu[index - 1] = (int) value;
256  cout << "mstu\t" << index << " " << value << endl;
257  IntegerTest(value);
258  }
259  else if (label == "ckin")
260  {
261  line >> index >> value;
262  pysubs.ckin[index - 1] = value;
263  cout << "ckin\t" << index << " " << value << endl;
264  }
265  else if (label == "parp")
266  {
267  line >> index >> value;
268  pypars.parp[index - 1] = value;
269  cout << "parp\t" << index << " " << value << endl;
270  }
271  else if (label == "parj")
272  {
273  line >> index >> value;
274  pydat1.parj[index - 1] = value;
275  cout << "parj\t" << index << " " << value << endl;
276  }
277  else if (label == "paru")
278  {
279  line >> index >> value;
280  pydat1.paru[index - 1] = value;
281  cout << "paru\t" << index << " " << value << endl;
282  }
283  else if (label == "parf")
284  {
285  line >> index >> value;
286  pydat2.parf[index - 1] = value;
287  cout << "parf\t" << index << " " << value << endl;
288  }
289  else if (label == "mdme")
290  {
291  int idc = 0; // decay channel
292  line >> idc >> index >> value;
293 
294  // if (ivalue==1/0) turn on/off decay channel idc
295  pydat3.mdme[index - 1][idc - 1] = value;
296  cout << "mdme\t" << idc << " " << index << " " << value << endl;
297  }
298  else if (label == "pmas")
299  {
300  int idc = 0;
301  line >> idc >> index >> value;
302 
303  pydat2.pmas[index - 1][idc - 1] = value;
304  cout << "pmas\t" << idc << " " << index << " " << value << endl;
305  }
306  else if (label == "pytune")
307  {
308  int ivalue;
309  line >> ivalue;
310  pytune(&ivalue);
311  cout << "pytune\t" << ivalue << endl;
312  }
313  else
314  {
315  // label was not understood
316  cout << "************************************************************" << endl;
317  cout << "PHPythia6::ReadConfig(), ERROR this option is not supported: " << FullLine << endl;
318  cout << "************************************************************" << endl;
319  }
320 
321  // get next line in file
322  getline(infile, FullLine);
323  }
324 
325  // Call pythia initialization
326  call_pyinit(_frame.c_str(), _proj.c_str(), _targ.c_str(), _roots);
327 
328  //call_pylist(12);
329 
330  infile.close();
331 
332  return _nevents;
333 }
334 
335 //-* print pythia config info
337 {
338 }
339 
341 {
342  if (Verbosity() > 1) cout << "PHPythia6::process_event - event: " << _eventcount << endl;
343 
344  bool passedTrigger = false;
345  std::vector<bool> theTriggerResults;
346  int genCounter = 0;
347 
348  /* based on HepMC/example_MyPythia.cc
349  *........................................HepMC INITIALIZATIONS
350  *
351  * Instantiate an IO strategy for reading from HEPEVT. */
352  HepMC::IO_HEPEVT hepevtio;
353  HepMC::GenEvent *evt;
354 
355  while (!passedTrigger)
356  {
357  ++genCounter;
358 
359  call_pyevnt(); // generate one event with Pythia
360  _geneventcount++;
361  // pythia pyhepc routine converts common PYJETS in common HEPEVT
362  call_pyhepc(1);
363  evt = hepevtio.read_next_event();
364 
365  // define the units (Pythia uses GeV and mm)
366  evt->use_units(HepMC::Units::GEV, HepMC::Units::MM);
367 
368  // add some information to the event
369  evt->set_event_number(_eventcount);
370 
371  /* process ID from pythia */
372  evt->set_signal_process_id(pypars.msti[1 - 1]);
373 
374  // set number of multi parton interactions
375  evt->set_mpi(pypars.msti[31 - 1]);
376 
377  // set cross section information
378  evt->set_cross_section(HepMC::getPythiaCrossSection());
379 
380  // Set the PDF information
381  HepMC::PdfInfo pdfinfo;
382  pdfinfo.set_x1(pypars.pari[33 - 1]);
383  pdfinfo.set_x2(pypars.pari[34 - 1]);
384  pdfinfo.set_scalePDF(pypars.pari[22 - 1]);
385  pdfinfo.set_id1(pypars.msti[15 - 1]);
386  pdfinfo.set_id2(pypars.msti[16 - 1]);
387  evt->set_pdf_info(pdfinfo);
388 
389  // test trigger logic
390 
391  bool andScoreKeeper = true;
392  if (Verbosity() > 2)
393  {
394  cout << "PHPythia6::process_event - triggersize: " << _registeredTriggers.size() << endl;
395  }
396 
397  for (unsigned int tr = 0; tr < _registeredTriggers.size(); tr++)
398  {
399  bool trigResult = _registeredTriggers[tr]->Apply(evt);
400 
401  if (Verbosity() > 2)
402  {
403  cout << "PHPythia6::process_event trigger: "
404  << _registeredTriggers[tr]->GetName() << " " << trigResult << endl;
405  }
406 
407  if (_triggersOR && trigResult)
408  {
409  passedTrigger = true;
410  break;
411  }
412  else if (_triggersAND)
413  {
414  andScoreKeeper &= trigResult;
415  }
416 
417  if (Verbosity() > 2 && !passedTrigger)
418  {
419  cout << "PHPythia8::process_event - failed trigger: "
420  << _registeredTriggers[tr]->GetName() << endl;
421  }
422  }
423 
424  if ((andScoreKeeper && _triggersAND) || (_registeredTriggers.size() == 0))
425  {
426  passedTrigger = true;
427  genCounter = 0;
428  }
429 
430  // delete failed events
431  if (!passedTrigger) delete evt;
432  }
433 
434  /* write the event out to the ascii files */
435  if (_save_ascii)
436  {
437  HepMC::IO_GenEvent ascii_io(_filename_ascii.c_str(), std::ios::out);
438  ascii_io << evt;
439  }
440 
441  /* pass HepMC to PHNode*/
442 
444  if (!success)
445  {
446  cout << "PHPythia6::process_event - Failed to add event to HepMC record!" << endl;
448  }
449  /* print outs*/
450  if (Verbosity() > 2) cout << "PHPythia6::process_event - FINISHED WHOLE EVENT" << endl;
451 
452  ++_eventcount;
454 }
455 
456 void PHPythia6::IntegerTest(double number)
457 {
458  if (fmod(number, 1.0) != 0)
459  {
460  cout << "Warning: Value " << number << " is not an integer." << endl;
461  cout << "This parameter requires an integer value." << endl;
462  cout << "Value of parameter truncated to " << (int) number << endl;
463 
464  //...End simulation if a double value is input for an integer parameter
465  // throw Fun4AllReturnCodes::ABORTRUN;
466  }
467  return;
468 }
469 
471 {
473 }
474 
476 {
477  if (Verbosity() > 1) cout << "PHPythia6::registerTrigger - trigger " << theTrigger->GetName() << " registered" << endl;
478  _registeredTriggers.push_back(theTrigger);
479 }