Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHPythia8.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHPythia8.cc
1 #include "PHPythia8.h"
2 
3 #include "PHPy8GenTrigger.h"
4 
5 #include <phhepmc/PHGenIntegral.h> // for PHGenIntegral
7 #include <phhepmc/PHHepMCGenHelper.h> // for PHHepMCGenHelper
8 
9 #include <fun4all/Fun4AllBase.h> // for Fun4AllBase::VERBO...
11 #include <fun4all/SubsysReco.h> // for SubsysReco
12 
13 #include <phool/PHCompositeNode.h>
14 #include <phool/PHIODataNode.h>
15 #include <phool/PHNode.h> // for PHNode
16 #include <phool/PHNodeIterator.h> // for PHNodeIterator
17 #include <phool/PHObject.h> // for PHObject
18 #include <phool/PHRandomSeed.h>
19 #include <phool/getClass.h>
20 #include <phool/phool.h> // for PHWHERE
21 
22 #pragma GCC diagnostic push
23 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
24 #include <HepMC/GenEvent.h>
25 #pragma GCC diagnostic pop
26 
27 #include <HepMC/Units.h> // for GEV, MM
28 #include <HepMC/WeightContainer.h> // for WeightContainer
29 
30 #include <Pythia8/Event.h> // for Event
31 #include <Pythia8/Info.h> // for Info
32 #include <Pythia8/Pythia.h>
33 #include <Pythia8Plugins/HepMC2.h>
34 
35 #include <boost/format.hpp>
36 
37 #include <cassert>
38 #include <cstdlib>
39 #include <iostream> // for operator<<, endl
40 
41 class PHHepMCGenEvent;
42 
43 using namespace std;
44 
46  : SubsysReco(name)
47  , m_EventCount(0)
48  , m_TriggersOR(true)
49  , m_TriggersAND(false)
50  , m_Pythia8(nullptr)
51  , m_ConfigFileName("phpythia8.cfg")
52  , m_Pythia8ToHepMC(nullptr)
53  , m_SaveEventWeightFlag(true)
54  , m_SaveIntegratedLuminosityFlag(true)
55  , m_IntegralNode(nullptr)
56 {
57  char *charPath = getenv("PYTHIA8");
58  if (!charPath)
59  {
60  cout << "PHPythia8::Could not find $PYTHIA8 path!" << endl;
61  return;
62  }
63 
64  std::string thePath(charPath);
65  thePath += "/xmldoc/";
66  m_Pythia8 = new Pythia8::Pythia(thePath.c_str());
67 
68  m_Pythia8ToHepMC = new HepMC::Pythia8ToHepMC();
69  m_Pythia8ToHepMC->set_store_proc(true);
70  m_Pythia8ToHepMC->set_store_pdf(true);
71  m_Pythia8ToHepMC->set_store_xsec(true);
72 
73  PHHepMCGenHelper::set_embedding_id(1); // default embedding ID to 1
74 }
75 
77 {
78  delete m_Pythia8;
79  delete m_Pythia8ToHepMC;
80 }
81 
83 {
84  if (!m_ConfigFileName.empty())
85  {
87  }
88  for (unsigned int j = 0; j < m_Commands.size(); j++)
89  {
90  m_Pythia8->readString(m_Commands[j]);
91  }
92 
93  create_node_tree(topNode);
94 
95  // PYTHIA8 has very specific requires for its random number range
96  // I map the designated unique seed from recoconst into something
97  // acceptable for PYTHIA8
98 
99  unsigned int seed = PHRandomSeed();
100 
101  if (seed > 900000000)
102  {
103  seed = seed % 900000000;
104  }
105 
106  if ((seed > 0) && (seed <= 900000000))
107  {
108  m_Pythia8->readString("Random:setSeed = on");
109  m_Pythia8->readString(str(boost::format("Random:seed = %1%") % seed));
110  }
111  else
112  {
113  cout << PHWHERE << " ERROR: seed " << seed << " is not valid" << endl;
114  exit(1);
115  }
116  // print out seed so we can make this is reproducible
117  cout << "PHPythia8 random seed: " << seed << endl;
118 
119  m_Pythia8->init();
120 
122 }
123 
124 int PHPythia8::End(PHCompositeNode * /*topNode*/)
125 {
126  if (Verbosity() >= VERBOSITY_MORE)
127  {
128  cout << "PHPythia8::End - I'm here!" << endl;
129  }
130 
131  if (Verbosity() >= VERBOSITY_SOME)
132  {
133  //-* dump out closing info (cross-sections, etc)
134  m_Pythia8->stat();
135 
136  //match pythia printout
137  cout << " | "
138  << " | " << endl;
139  cout << " PHPythia8::End - " << m_EventCount
140  << " events passed trigger" << endl;
141  cout << " Fraction passed: " << m_EventCount
142  << "/" << m_Pythia8->info.nAccepted()
143  << " = " << m_EventCount / float(m_Pythia8->info.nAccepted()) << endl;
144  cout << " *------- End PYTHIA Trigger Statistics ------------------------"
145  << "-------------------------------------------------* " << endl;
146 
147  if (m_IntegralNode)
148  {
149  cout << "Integral information on stored on node RUN/PHGenIntegral:" << endl;
151  cout << " *------- End PYTHIA Integral Node Print ------------------------"
152  << "-------------------------------------------------* " << endl;
153  }
154  }
155 
157 }
158 
159 //__________________________________________________________
160 int PHPythia8::read_config(const string &cfg_file)
161 {
162  m_ConfigFileName = cfg_file;
163 
164  if (Verbosity() >= VERBOSITY_SOME)
165  {
166  cout << "PHPythia8::read_config - Reading " << m_ConfigFileName << endl;
167  }
168 
169  ifstream infile(m_ConfigFileName);
170  if (infile.fail())
171  {
172  cout << "PHPythia8::read_config - Failed to open file " << m_ConfigFileName << endl;
173  exit(2);
174  }
175 
176  m_Pythia8->readFile(m_ConfigFileName);
177 
179 }
180 
181 //-* print pythia config info
183 {
184  m_Pythia8->info.list();
185 }
186 
188 {
189  if (Verbosity() >= VERBOSITY_MORE)
190  {
191  cout << "PHPythia8::process_event - event: " << m_EventCount << endl;
192  }
193 
194  bool passedGen = false;
195  bool passedTrigger = false;
196  int genCounter = 0;
197 
198  while (!passedTrigger)
199  {
200  ++genCounter;
201 
202  // generate another pythia event
203  while (!passedGen)
204  {
205  passedGen = m_Pythia8->next();
206  }
207 
208  // test trigger logic
209 
210  bool andScoreKeeper = true;
212  {
213  cout << "PHPythia8::process_event - triggersize: " << m_RegisteredTriggers.size() << endl;
214  }
215 
216  for (unsigned int tr = 0; tr < m_RegisteredTriggers.size(); tr++)
217  {
218  bool trigResult = m_RegisteredTriggers[tr]->Apply(m_Pythia8);
219 
221  {
222  cout << "PHPythia8::process_event trigger: "
223  << m_RegisteredTriggers[tr]->GetName() << " " << trigResult << endl;
224  }
225 
226  if (m_TriggersOR && trigResult)
227  {
228  passedTrigger = true;
229  break;
230  }
231  else if (m_TriggersAND)
232  {
233  andScoreKeeper &= trigResult;
234  }
235 
236  if (Verbosity() >= VERBOSITY_EVEN_MORE && !passedTrigger)
237  {
238  cout << "PHPythia8::process_event - failed trigger: "
239  << m_RegisteredTriggers[tr]->GetName() << endl;
240  }
241  }
242 
243  if ((andScoreKeeper && m_TriggersAND) || (m_RegisteredTriggers.size() == 0))
244  {
245  passedTrigger = true;
246  genCounter = 0;
247  }
248 
249  passedGen = false;
250  }
251 
252  // fill HepMC object with event & pass to
253 
254  HepMC::GenEvent *genevent = new HepMC::GenEvent(HepMC::Units::GEV, HepMC::Units::MM);
255  m_Pythia8ToHepMC->fill_next_event(*m_Pythia8, genevent, m_EventCount);
256  // Enable continuous reweighting by storing additional reweighting factor
258  {
259  genevent->weights().push_back(m_Pythia8->info.weight());
260  }
261 
262  /* pass HepMC to PHNode*/
263  PHHepMCGenEvent *success = PHHepMCGenHelper::insert_event(genevent);
264  if (!success)
265  {
266  cout << "PHPythia8::process_event - Failed to add event to HepMC record!" << endl;
268  }
269 
270  // print outs
271 
272  if (Verbosity() >= VERBOSITY_MORE)
273  {
274  cout << "PHPythia8::process_event - FINISHED WHOLE EVENT" << endl;
275  }
276  if (m_EventCount < 2 && Verbosity() >= VERBOSITY_SOME)
277  {
278  m_Pythia8->event.list();
279  }
280  if (m_EventCount >= 2 && Verbosity() >= VERBOSITY_A_LOT)
281  {
282  m_Pythia8->event.list();
283  }
284 
285  ++m_EventCount;
286 
287  // save statistics
288  if (m_IntegralNode)
289  {
292  m_IntegralNode->set_Sum_Of_Weight(m_Pythia8->info.weightSum());
293  m_IntegralNode->set_Integrated_Lumi(m_Pythia8->info.nAccepted() / (m_Pythia8->info.sigmaGen() * 1e9));
294  }
295 
297 }
298 
300 {
301  // HepMC IO
303 
304  PHNodeIterator iter(topNode);
305  PHCompositeNode *sumNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
306  if (!sumNode)
307  {
308  cout << PHWHERE << "RUN Node missing doing nothing" << endl;
310  }
312  {
313  m_IntegralNode = findNode::getClass<PHGenIntegral>(sumNode, "PHGenIntegral");
314  if (!m_IntegralNode)
315  {
316  m_IntegralNode = new PHGenIntegralv1("PHPythia8 with embedding ID of " + std::to_string(PHHepMCGenHelper::get_embedding_id()));
317  PHIODataNode<PHObject> *newmapnode = new PHIODataNode<PHObject>(m_IntegralNode, "PHGenIntegral", "PHObject");
318  sumNode->addNode(newmapnode);
319  }
320  else
321  {
322  cout << "PHPythia8::create_node_tree - Fatal Error - "
323  << "RUN/PHGenIntegral node already exist. "
324  << "It is messy to overwrite integrated luminosities. Please turn off this function in the macro with " << endl;
325  cout << " PHPythia8::save_integrated_luminosity(false);" << endl;
326  cout << "The current RUN/PHGenIntegral node is ";
327  m_IntegralNode->identify(cout);
328 
329  exit(EXIT_FAILURE);
330  }
332  }
334 }
335 
337 {
338  if (Verbosity() >= VERBOSITY_MORE)
339  {
340  cout << "PHPythia8::registerTrigger - trigger " << theTrigger->GetName() << " registered" << endl;
341  }
342  m_RegisteredTriggers.push_back(theTrigger);
343 }