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