ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Hijing.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Hijing.cxx
1 // Allows the user to generate hijing events and store the results in
2 // a HepMC file.
3 //
4 // Inspired by code from ATLAS. Thanks!
5 //
6 #include "HiMain1.h"
7 #include "HiMain2.h"
8 #include "HiParnt.h"
9 #include "HiStrng.h"
10 #include "HijCrdn.h"
11 #include "HijJet1.h"
12 #include "HijJet2.h"
13 #include "HijJet4.h"
14 #include "RanSeed.h"
15 
16 #include <HepMC/GenEvent.h>
17 #include <HepMC/GenParticle.h>
18 #include <HepMC/GenVertex.h>
19 #include <HepMC/IO_AsciiParticles.h>
20 #include <HepMC/IO_GenEvent.h>
21 
22 #include <CLHEP/Geometry/Point3D.h>
24 #include <CLHEP/Random/RandFlat.h>
27 
28 #include <boost/algorithm/string.hpp>
29 #include <boost/foreach.hpp>
30 #include <boost/lexical_cast.hpp>
31 #include <boost/property_tree/ptree.hpp>
32 #include <boost/property_tree/xml_parser.hpp>
33 
34 #include <cmath>
35 #include <cstdlib>
36 #include <iostream>
37 #include <memory>
38 #include <numeric>
39 #include <random>
40 #include <string>
41 #include <vector>
42 
43 #define f2cFortran
44 #define gFortran
45 // cppcheck-suppress *
46 #include <cfortran.h>
47 
48 PROTOCCALLSFSUB3(HIJING, hijing, STRING, FLOAT, FLOAT)
49 #define HIJING(FRAME, BMIN0, BMAX0) \
50  CCALLSFSUB3(HIJING, hijing, STRING, FLOAT, FLOAT, FRAME, BMIN0, BMAX0)
51 
52 PROTOCCALLSFSUB8(HIJSET, hijset, FLOAT, STRING, STRING, STRING, INT, INT, INT, INT)
53 #define HIJSET(EFRM, FRAME, PROJ, TARG, IAP, IZP, IAT, IZT) \
54  CCALLSFSUB8(HIJSET, hijset, FLOAT, STRING, STRING, STRING, INT, INT, INT, INT, \
55  EFRM, FRAME, PROJ, TARG, IAP, IZP, IAT, IZT)
56 
57 using namespace std;
58 using namespace boost;
59 using namespace boost::property_tree;
60 
61 void hijfst_control(int, vector<string>, vector<float>, vector<int>, vector<float>, vector<float>, vector<float>);
62 
64 
65 float atl_ran(int *)
66 {
67  return (float) CLHEP::RandFlat::shoot(engine);
68 }
69 // This prevents cppcheck to flag the next line as error
70 // cppcheck-suppress *
71 FCALLSCFUN1(FLOAT, atl_ran, ATL_RAN, atl_ran, PINT)
72 
73 typedef HepGeom::Point3D<double> HepPoint3D;
74 typedef HepMC::GenEvent::particle_iterator piter;
75 
76 int fillEvent(HepMC::GenEvent *evt);
77 
78 // Accessor to HIJING Options and parameters COMMON
80 
81 // Accessor to HIJING Random number seed COMMON
83 
84 // Accessors to HIJING Event information COMMONS
92 
93 float efrm;
94 std::string m_frame;
95 std::string m_proj;
96 std::string m_targ;
97 int iap;
98 int iat;
99 int izp;
100 int izt;
101 int spec;
102 double m_vertexOffsetCut = 1.0E-7;
104 
105 int main(int argc, char **argv)
106 {
107  string config_filename = "sHijing.xml";
108  string output;
109  long randomSeed = 0;
110  bool randomseed_set = false;
111  unsigned int N = 1;
112  bool NEvents_set = false;
113  for (int i = 1; i < argc; ++i)
114  {
115  std::string optionstring = argv[i];
116  if (optionstring == "-h")
117  {
118  cout << endl
119  << "Usage: sHijing <config xmlfile [sHijing.xml]>" << endl;
120  cout << endl;
121  cout << "Parameters:" << endl;
122  cout << "-n <number of events [1]>" << endl;
123  cout << "-o <outputfile [sHijing.dat]>" << endl;
124  cout << "-s <random seet [std::random_device]>" << endl;
125  exit(0);
126  }
127  else if (optionstring == "-o")
128  {
129  if (i + 1 < argc)
130  { // Make sure we aren't at the end of argv!
131  output = argv[++i]; // Increment 'i' so we get the argument
132  }
133  else
134  { // Uh-oh, there was no argument to the destination option.
135  std::cerr << "-o option requires one argument." << std::endl;
136  exit(1);
137  }
138  continue;
139  }
140  else if (optionstring == "-s")
141  {
142  if (i + 1 < argc)
143  { // Make sure we aren't at the end of argv!
144  randomSeed = std::stol(argv[++i]); // Increment 'i' so get the argument.
145  randomseed_set = true;
146  }
147  else
148  { // Uh-oh, there was no argument to the destination option.
149  std::cerr << "-s option requires one argument." << std::endl;
150  exit(1);
151  }
152  continue;
153  }
154  else if (optionstring == "-n")
155  {
156  if (i + 1 < argc)
157  { // Make sure we aren't at the end of argv!
158  N = std::stoul(argv[++i]); // Increment 'i' so get the argument.
159  NEvents_set = true;
160  }
161  else
162  { // Uh-oh, there was no argument to the destination option.
163  std::cerr << "-s option requires one argument." << std::endl;
164  exit(1);
165  }
166  continue;
167  }
168  else
169  {
170  config_filename = argv[i];
171  }
172  }
173  char frame[] = " ";
174  char proj[] = " ";
175  char targ[] = " ";
176 
177  using boost::property_tree::ptree;
178  iptree pt, null;
179  std::ifstream config_file(config_filename);
180 
181  if (config_file)
182  {
183  // Read XML configuration file.
184  read_xml(config_file, pt);
185  cout << "using config file: " << config_filename << endl;
186  }
187  else
188  {
189  cout << "no xml config file - using internal values" << endl;
190  }
191  efrm = pt.get("HIJING.EFRM", 200.0);
192  m_frame = pt.get("HIJING.FRAME", "CMS");
193  m_proj = pt.get("HIJING.PROJ", "A");
194  m_targ = pt.get("HIJING.TARG", "A");
195  iap = pt.get("HIJING.IAP", 197);
196  izp = pt.get("HIJING.IZP", 79);
197  iat = pt.get("HIJING.IAT", 197);
198  izt = pt.get("HIJING.IZT", 79);
199  float bmin = pt.get("HIJING.BMIN", 0.0);
200  float bmax = pt.get("HIJING.BMAX", 0.0);
201  if (!NEvents_set)
202  {
203  N = pt.get("HIJING.N", 1);
204  }
205  keepSpectators = pt.get("HIJING.KEEP_SPECTATORS", 1);
206  if (output.empty())
207  {
208  output = pt.get("HIJING.OUTPUT", "sHijing.dat");
209  }
210  if (!randomseed_set)
211  {
212  std::random_device rdev;
213  randomSeed = pt.get("HIJING.RANDOM.SEED", rdev());
214  }
215  engine = new CLHEP::MTwistEngine(randomSeed);
216 
217  // See if there are any sections for HIPR1, IHPR2
218  iptree &pr1 = pt.get_child("HIJING.HIPR1", null);
219  BOOST_FOREACH (iptree::value_type &v, pr1)
220  {
221  int key = boost::lexical_cast<int>(v.first.data());
222  float value = boost::lexical_cast<float>(v.second.data());
223  m_hiparnt.hipr1(key) = value;
224  }
225  iptree &pr2 = pt.get_child("HIJING.IHPR2", null);
226  BOOST_FOREACH (iptree::value_type &v, pr2)
227  {
228  int key = boost::lexical_cast<int>(v.first.data());
229  int value = boost::lexical_cast<int>(v.second.data());
230  m_hiparnt.ihpr2(key) = value;
231  }
232 
233  int fastjet_enable_p = 0;
234  std::vector<string> algorithm_v;
235  std::vector<float> R_v;
236  std::vector<int> PID_v;
237  std::vector<float> EtaMin_v;
238  std::vector<float> EtaMax_v;
239  std::vector<float> EtMin_v;
240 
241  iptree &it = pt.get_child("HIJING.FASTJET", null);
242  BOOST_FOREACH (iptree::value_type &v, it)
243  {
244  if (to_upper_copy(v.first) != "ALGORITHM") continue;
245  algorithm_v.push_back(to_upper_copy(v.second.get("NAME", "ANTIKT")));
246  R_v.push_back(v.second.get("R", 0.2));
247  PID_v.push_back(v.second.get("PID", 2000000));
248 
249  EtaMin_v.push_back(v.second.get("EtaMin", -2));
250  EtaMax_v.push_back(v.second.get("EtaMax", 2));
251  EtMin_v.push_back(v.second.get("EtMin", 5));
252 
253  fastjet_enable_p = 1;
254  }
255  cout << "seed: " << randomSeed << endl;
256  cout << "output: " << output << endl;
257  cout << "Number of Events: " << N << endl;
258 
259  hijfst_control(fastjet_enable_p, algorithm_v, R_v, PID_v, EtaMin_v, EtaMax_v, EtMin_v);
260 
261  // The call to Hijing needs simple C-style strings.
262  m_frame.copy(frame, m_frame.size());
263  m_proj.copy(proj, m_proj.size());
264  m_targ.copy(targ, m_targ.size());
265 
266  HIJSET(efrm, frame, proj, targ, iap, izp, iat, izt);
267 
268  // int status;
269  HepMC::IO_GenEvent ascii_io(output.c_str(), std::ios::out);
270  unsigned int events = 0;
271  do
272  {
273  HIJING(frame, bmin, bmax);
274 
275  if (m_himain1.ierrstat() != 0)
276  {
277  continue;
278  }
279  events++;
280 
281  HepMC::GenEvent *evt = new HepMC::GenEvent();
282  evt->use_units(HepMC::Units::GEV, HepMC::Units::MM);
283  evt->set_event_number(events);
284 
285  fillEvent(evt);
286 
287  ascii_io << evt;
288 
289  delete evt;
290  } while (events < N);
291 
292  return 0;
293 }
294 
295 int fillEvent(HepMC::GenEvent *evt)
296 {
297  int np = m_himain1.np();
298  int nt = m_himain1.nt();
299  int n0 = m_himain1.n0();
300  int n01 = m_himain1.n01();
301  int n10 = m_himain1.n10();
302  int n11 = m_himain1.n11();
303 
304  // The Hijing documentation is a little unclear about this, but the
305  // four values, n0, n01, n10 and n11, count collisions between
306  // nucleons that have already (or have not already had) collisions
307  // with other nucleons. N_coll, the way we use it, is the sum of
308  // all these values.
309  int ncoll = n0 + n01 + n10 + n11;
310 
311  int natt = m_himain1.natt();
312  float b = m_hiparnt.hint1(19);
313  float bphi = m_hiparnt.hint1(20);
314 
315  // Determine the participant eccentricity.
316  std::vector<float> x;
317  std::vector<float> y;
318  float tx, ty;
319  float bbx = b * cos(bphi);
320  float bby = b * sin(bphi);
321  // std::cout << "x,y,c" << std::endl;
322  for (int i = 0; i < m_hiparnt.ihnt2(1); i++)
323  {
324  if (m_histrng.nfp(i, 5) > 2)
325  {
326  tx = m_hijcrdn.yp(1, i) + bbx;
327  ty = m_hijcrdn.yp(2, i) + bby;
328  x.push_back(tx);
329  y.push_back(ty);
330 
331  // std::cout << tx << "," << ty << "," << "0" << std::endl;
332  }
333  }
334  for (int i = 0; i < m_hiparnt.ihnt2(3); i++)
335  {
336  if (m_histrng.nft(i, 5) > 2)
337  {
338  tx = m_hijcrdn.yt(1, i);
339  ty = m_hijcrdn.yt(2, i);
340  x.push_back(tx);
341  y.push_back(ty);
342 
343  // std::cout << tx << "," << ty << "," << "1" << std::endl;
344  }
345  }
346 
347  float e_part = 0.0;
348  if (x.size() != 0)
349  {
350  float N = x.size();
351  float xa = std::accumulate(x.begin(), x.end(), 0) / N;
352  float ya = std::accumulate(y.begin(), y.end(), 0) / N;
353 
354  float sxx = 0;
355  float syy = 0;
356  float sxy = 0;
357  for (int i = 0; i < N; i++)
358  {
359  sxx += (x[i] - xa) * (x[i] - xa);
360  syy += (y[i] - ya) * (y[i] - ya);
361  sxy += (x[i] - xa) * (y[i] - ya);
362  }
363  sxx /= N;
364  syy /= N;
365  sxy /= N;
366  e_part = std::sqrt((sxx - syy) * (sxx - syy) + 4 * sxy * sxy) / (sxx + syy);
367  }
368 
369  // Need to calculate a few things missing from this list
370  HepMC::HeavyIon hi(1, np, nt, ncoll, 0, 0, n01, n10, n11, b, bphi, e_part, 42.0);
371 
372  evt->set_heavy_ion(hi);
373 
374  // Set the random generator seeds. How do people handle the fact
375  // that CLHEP produced unsigned longs and HepMC wants signed ones?
376  // evt->set_random_states(engine->put());
377 
378  // Did we keep decay history?
379  bool keptHistory = (m_hiparnt.ihpr2(21) == 1);
380 
381  // The number of particles in the hijing output
382  int numHijingPart = m_himain1.natt();
383 
384  // Vectors that will keep track of where particles originate from and die
385  std::vector<int> partOriginVertex_vec;
386  std::vector<int> partDecayVertex_vec;
387  std::vector<HepMC::GenParticle *> particleHepPartPtr_vec;
388 
389  partOriginVertex_vec.assign(numHijingPart, 0);
390  partDecayVertex_vec.assign(numHijingPart, -1);
391  particleHepPartPtr_vec.assign(numHijingPart, (HepMC::GenParticle *) 0);
392 
393  // Vector that will keep pointers to generated vertices
394  std::vector<HepMC::GenVertex *> vertexPtrVec;
395 
396  // Create the event vertex
397  CLHEP::HepLorentzVector newVertex;
398  newVertex = CLHEP::HepLorentzVector(0., 0., 0., 0.);
399  HepMC::GenVertex *v1 = new HepMC::GenVertex(newVertex);
400 
401  evt->add_vertex(v1);
402  vertexPtrVec.push_back(v1);
403 
404  double eproj = (m_frame == "CMS") ? efrm / 2.0 : efrm;
405  int proj_id = (m_proj == "A") ? 3000000 + iap : 2212;
406  v1->add_particle_in(new HepMC::GenParticle(CLHEP::HepLorentzVector(0., 0., eproj, eproj),
407  proj_id, 101));
408 
409  double etarg = (m_frame == "CMS") ? efrm / 2.0 : efrm;
410  int targ_id = (m_targ == "A") ? 3000000 + iap : 2212;
411  v1->add_particle_in(new HepMC::GenParticle(CLHEP::HepLorentzVector(0., 0., -etarg, etarg),
412  targ_id, 102));
413 
414  // Loop on all Hijing particles and put them all as outgoing from the event vertex
415  for (int i = 1; i <= natt; ++i)
416  {
417  if (m_himain2.katt(i, 4) == 103)
418  {
419  // We handle jets a little differently.
421  m_himain2.patt(i, 2),
422  m_himain2.patt(i, 3),
423  m_himain2.patt(i, 4));
424  v1->add_particle_out(new HepMC::GenParticle(jetP4,
425  m_himain2.katt(i, 1),
426  m_himain2.katt(i, 4)));
427 
428  continue;
429  }
430 
431  // Skip non-interacting projectile and target nucleons
432  if (!keepSpectators &&
433  ((m_himain2.katt(i, 2) == 0) || (m_himain2.katt(i, 2)) == 10))
434  continue;
435 
436  // Find the vertex of the parent particle
437  int parentIndex = m_himain2.katt(i, 3) - 1;
438  int parentOriginIndex = 0;
439  int parentDecayIndex = -1;
440 
441  // If the particle has a true parent, record its vertex info
442  if (parentIndex >= 0)
443  {
444  parentOriginIndex = partOriginVertex_vec[parentIndex];
445  parentDecayIndex = partDecayVertex_vec[parentIndex];
446  }
447 
448  // A CLHEP vector containing the particle start point
449  CLHEP::HepLorentzVector particleStart(m_himain2.vatt(i, 1),
450  m_himain2.vatt(i, 2),
451  m_himain2.vatt(i, 3),
452  m_himain2.vatt(i, 4));
453 
454  // By default, the particle originates from the primary vertex
455  int particleVertexIndex = 0;
456 
457  // Have we kept the history?
458  if (keptHistory)
459  {
460  // Check to see if we've already generated a decay vertex
461  // for this parent
462  if (parentDecayIndex != -1)
463  {
464  // Make sure it is consistent with this track origin
465  HepPoint3D vertex_pos(vertexPtrVec[parentDecayIndex]->point3d().x(),
466  vertexPtrVec[parentDecayIndex]->point3d().y(),
467  vertexPtrVec[parentDecayIndex]->point3d().z());
468  double distance = vertex_pos.distance(particleStart.vect());
469 
470  if (distance > m_vertexOffsetCut)
471  {
472  HepMC::GenVertex::particles_out_const_iterator iter;
473  for (iter = vertexPtrVec[parentDecayIndex]->particles_out_const_begin();
474  iter != vertexPtrVec[parentDecayIndex]->particles_out_const_end();
475  iter++)
476  {
477  std::cout << (*iter)->barcode() << ", ";
478  }
479 
480  std::cout << std::endl;
481  }
482 
483  // Nonetheless set the parent decay vertex to be this particle's vertex
484  particleVertexIndex = parentDecayIndex;
485  }
486  else
487  {
488  // Now compare the distance between the vertex FROM
489  // which the parent originates and the start of this
490  // particle
491  HepPoint3D vertex_pos(vertexPtrVec[parentOriginIndex]->point3d().x(),
492  vertexPtrVec[parentOriginIndex]->point3d().y(),
493  vertexPtrVec[parentOriginIndex]->point3d().z());
494  double distance = vertex_pos.distance(particleStart.vect());
495 
496  if (distance > m_vertexOffsetCut && parentIndex == -1)
497  {
498  // *** Explicitly handle Hijing bug which generates
499  // *** particle with displaced vertex but no parent
500  // *** ????
501 
502  std::cout
503  << "HIJING BUG:: Particle found with displaced vertex but no parent "
504  << std::endl;
505  }
506  if (distance > m_vertexOffsetCut || parentOriginIndex != 0)
507  {
508  // We need to create a new vertex
509  HepMC::GenVertex *newVertex_p = new HepMC::GenVertex(particleStart);
510 
511  evt->add_vertex(newVertex_p);
512  vertexPtrVec.push_back(newVertex_p);
513  particleVertexIndex = vertexPtrVec.size() - 1;
514 
515  // Now we indicate that the parent has a decay vertex
516  partDecayVertex_vec[parentIndex] = particleVertexIndex;
517 
518  // Now tell the vertex about the particle that created it
519  newVertex_p->add_particle_in(particleHepPartPtr_vec[parentIndex]);
520  }
521  else
522  {
523  // Assign the particle to its parent's vertex
524  particleVertexIndex = parentOriginIndex;
525  }
526  }
527  }
528  else
529  {
530  // We have to brute-force search for a vertex that might match this particle
531  int foundVert = -1;
532  for (unsigned int ivert = 0; ivert < vertexPtrVec.size(); ivert++)
533  {
534  HepPoint3D vertex_pos(vertexPtrVec[ivert]->point3d().x(),
535  vertexPtrVec[ivert]->point3d().y(),
536  vertexPtrVec[ivert]->point3d().z());
537  double distance = vertex_pos.distance(particleStart.vect());
538  if (distance < m_vertexOffsetCut)
539  {
540  foundVert = ivert;
541  break;
542  }
543  }
544 
545  if (foundVert == -1)
546  {
547  // We need to create a new vertex
548  HepMC::GenVertex *newVertex_p = new HepMC::GenVertex(particleStart);
549 
550  evt->add_vertex(newVertex_p);
551  vertexPtrVec.push_back(newVertex_p);
552  particleVertexIndex = vertexPtrVec.size() - 1;
553  }
554  else
555  {
556  particleVertexIndex = foundVert;
557  }
558  }
559 
560  // If the Hijing particle has decayed, set the status appropriately
561  int particleId = m_himain2.katt(i, 1);
562  int particleStatus = 1;
563  if (m_himain2.katt(i, 4) == 11)
564  {
565  particleStatus = 2;
566  }
567 
568  // Create the new particle
569  CLHEP::HepLorentzVector particleP4(m_himain2.patt(i, 1),
570  m_himain2.patt(i, 2),
571  m_himain2.patt(i, 3),
572  m_himain2.patt(i, 4));
573 
574  HepMC::GenParticle *newParticle_p = new HepMC::GenParticle(particleP4,
575  particleId,
576  particleStatus);
577 
578  // Record the particle in the vector of pointers (ostensibly we
579  // only need this when we have the history but for simplicity
580  // always do it)
581  particleHepPartPtr_vec[i - 1] = newParticle_p;
582 
583  // Now add the particle to its vertex
584  vertexPtrVec[particleVertexIndex]->add_particle_out(newParticle_p);
585  partOriginVertex_vec[i - 1] = particleVertexIndex;
586  }
587 
588  return 0;
589 }