ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4AllHepMCPileupInputManager.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4AllHepMCPileupInputManager.cc
2 
3 #include "PHHepMCGenEvent.h"
4 #include "PHHepMCGenEventMap.h"
5 #include "PHHepMCGenHelper.h" // for PHHepMCGenHelper, PHHepMCGen...
6 
7 #include <fun4all/Fun4AllBase.h> // for Fun4AllBase::VERBOSITY_SOME
9 
10 #include <phool/PHRandomSeed.h>
11 
12 #include <HepMC/GenEvent.h>
13 #include <HepMC/IO_GenEvent.h>
14 
15 #include <gsl/gsl_randist.h>
16 #include <gsl/gsl_rng.h>
17 
18 #include <cassert> // for assert
19 #include <fstream>
20 #include <iostream>
21 
23  const std::string &name, const std::string &nodename, const std::string &topnodename)
24  : Fun4AllHepMCInputManager(name, nodename, topnodename)
25 {
27  Repeat(1);
28 
31 
32  RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
33  unsigned int seed = PHRandomSeed(); // fixed seed is handled in this funtcion
34  gsl_rng_set(RandomGenerator, seed);
35 
36  return;
37 }
38 
40 {
41  gsl_rng_free(RandomGenerator);
42 }
43 
45 {
46  for (int i = 0; i < nevents; ++i)
47  {
49  {
51  }
52  int iret = run(1, true);
53  if (iret)
54  {
55  return iret;
56  }
57  }
58  return 0;
59 }
60 
61 int Fun4AllHepMCPileupInputManager::run(const int /*nevents*/, const bool skip)
62 {
63  if (_first_run)
64  {
65  _first_run = false;
66 
70 
71  if (Verbosity() >= VERBOSITY_SOME)
72  {
73  std::cout << "Fun4AllHepMCPileupInputManager::run:first event: - ";
74  std::cout << " _ave_coll_per_crossing = " << _ave_coll_per_crossing;
75  std::cout << " _min_crossing = " << _min_crossing;
76  std::cout << " _max_crossing = " << _max_crossing;
77  std::cout << ". Start first event." << std::endl;
78  }
79  }
80  if (m_SignalInputManager && !skip)
81  {
83  }
84  // if we have pushed back events we read them here and return
86  {
87  // if m_EventPushedBackFlag = -1 we have no events but still created the temporary HepMC file
88  // so we still need to remove it
89  // if m_EventPushedBackFlag > 0 we have events pushed back which need to be recovered
90  if (m_EventPushedBackFlag > 0)
91  {
92  HepMC::IO_GenEvent ascii_tmp_in(m_HepMCTmpFile, std::ios::in);
93  HepMC::GenEvent *evttmp = ascii_tmp_in.read_next_event();
94  while (ascii_tmp_in.rdstate() == 0)
95  {
96  if (evttmp->event_number() != m_SignalEventNumber)
97  {
98  double crossing_time = m_EventNumberMap.find(evttmp->event_number())->second;
99  InsertEvent(evttmp, crossing_time);
100  }
101  evttmp = ascii_tmp_in.read_next_event();
102  }
103  }
104  m_EventNumberMap.clear();
106  remove(m_HepMCTmpFile.c_str());
107  return 0;
108  }
109  // toss multiple crossings all the way back
110  for (int icrossing = _min_crossing; icrossing <= _max_crossing; ++icrossing)
111  {
112  double crossing_time = _time_between_crossings * icrossing;
113 
114  int ncollisions = gsl_ran_poisson(RandomGenerator, _ave_coll_per_crossing);
115  // if (icrossing == 0) --ncollisions; // Jin: this is not necessary.
116  // Triggering an interesting crossing does not change the distribution of the number of background collisions in that crossing
117 
118  for (int icollision = 0; icollision < ncollisions; ++icollision)
119  {
120  // loop until retrieve a valid event
121  while (true)
122  {
123  if (!IsOpen())
124  {
125  if (FileListEmpty())
126  {
127  if (Verbosity() > 0)
128  {
129  std::cout << Name() << ": No Input file open" << std::endl;
130  }
131  return -1;
132  }
133  else
134  {
135  if (OpenNextFile())
136  {
137  std::cout << Name() << ": No Input file from filelist opened" << std::endl;
138  return -1;
139  }
140  }
141  }
142  {
143  if (ReadOscar())
144  {
145  evt = ConvertFromOscar();
146  if (evt && m_SignalEventNumber == evt->event_number())
147  {
148  delete evt;
149  evt = ConvertFromOscar();
150  }
151  }
152  else
153  {
154  evt = ascii_in->read_next_event();
155  if (evt && m_SignalEventNumber == evt->event_number())
156  {
157  delete evt;
158  evt = ascii_in->read_next_event();
159  }
160  }
161  }
162 
163  if (!evt)
164  {
165  if (Verbosity() > 1)
166  {
167  std::cout << "error type: " << ascii_in->error_type()
168  << ", rdstate: " << ascii_in->rdstate() << std::endl;
169  }
170  fileclose();
171  }
172  else
173  {
174  if (Verbosity() > 0)
175  {
176  std::cout << "Fun4AllHepMCPileupInputManager::run::" << Name();
177  if (skip) std::cout << " skip";
178  std::cout << " hepmc evt no: " << evt->event_number() << std::endl;
179  }
180  events_total++;
181  events_thisfile++;
182  // check if the local SubsysReco discards this event
184  {
185  delete evt;
186  evt = nullptr;
187  ResetEvent();
188  }
189  else
190  {
191  m_EventNumberMap.insert(std::make_pair(evt->event_number(), crossing_time));
192  break; // got the evt, move on
193  }
194  }
195  } // loop until retrieve a valid event
196  if (!skip)
197  {
198  InsertEvent(evt, crossing_time);
199  }
200  else
201  {
202  delete evt;
203  evt = nullptr;
204  }
205  } // for (int icollision = 0; icollision < ncollisions; ++icollision)
206 
207  } // for (int icrossing = _min_crossing; icrossing <= _max_crossing; ++icrossing)
208  return 0;
209 }
210 
212 {
213  m_EventNumberMap.clear();
215  return 0;
216 }
217 
219 {
220  if (i == 1)
221  {
222  if (m_HepMCTmpFile.empty())
223  {
224  // we need to create this filename just once, we reuse it. Do it only if we need it
225  m_HepMCTmpFile = "/tmp/HepMCTmpPileEvent-" + Name() + "-" + std::to_string(getpid()) + ".hepmc";
226  }
228  HepMC::IO_GenEvent ascii_io(m_HepMCTmpFile, std::ios::out);
230  for (auto iter = geneventmap->begin(); iter != geneventmap->end(); ++iter)
231  {
232  if (m_EventNumberMap.find((iter->second)->getEvent()->event_number()) != m_EventNumberMap.end())
233  {
235  HepMC::GenEvent *evttmp = (iter->second)->getEvent();
236  ascii_io << evttmp;
237  }
238  }
239  return 0;
240  }
241  return -1;
242 }
243 int Fun4AllHepMCPileupInputManager::InsertEvent(HepMC::GenEvent *evt, const double crossing_time)
244 {
246  PHHepMCGenEvent *genevent = nullptr;
248  {
250 
251  genevent = geneventmap->insert_active_event(get_PHHepMCGenEvent_template() );
252  }
253  else
254  {
256 
257  genevent = geneventmap->insert_background_event(get_PHHepMCGenEvent_template() );
258  }
259  assert(genevent);
260  assert(evt);
261  genevent->addEvent(evt);
263  // place to the crossing center in time
264  genevent->moveVertex(0, 0, 0, crossing_time);
265  return 0;
266 }