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>
22 #include <CLHEP/Geometry/Point3D.h>
23 #include <CLHEP/Random/MTwistEngine.h>
24 #include <CLHEP/Random/RandFlat.h>
25 #include <CLHEP/Random/RandomEngine.h>
26 #include <CLHEP/Vector/LorentzVector.h>
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>
49 #define HIJING(FRAME, BMIN0, BMAX0) \
50 CCALLSFSUB3(HIJING, hijing, STRING, FLOAT, FLOAT, FRAME, BMIN0, BMAX0)
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)
58 using namespace boost;
59 using namespace boost::property_tree;
61 void hijfst_control(
int, vector<string>, vector<float>, vector<int>, vector<float>, vector<float>, vector<float>);
67 return (
float) CLHEP::RandFlat::shoot(
engine);
74 typedef HepMC::GenEvent::particle_iterator
piter;
105 int main(
int argc,
char **argv)
107 string config_filename =
"sHijing.xml";
110 bool randomseed_set =
false;
112 bool NEvents_set =
false;
113 for (
int i = 1;
i < argc; ++
i)
116 if (optionstring ==
"-h")
119 <<
"Usage: sHijing <config xmlfile [sHijing.xml]>" << 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;
127 else if (optionstring ==
"-o")
135 std::cerr <<
"-o option requires one argument." << std::endl;
140 else if (optionstring ==
"-s")
144 randomSeed = std::stol(argv[++
i]);
145 randomseed_set =
true;
149 std::cerr <<
"-s option requires one argument." << std::endl;
154 else if (optionstring ==
"-n")
158 N = std::stoul(argv[++
i]);
163 std::cerr <<
"-s option requires one argument." << std::endl;
170 config_filename = argv[
i];
177 using boost::property_tree::ptree;
184 read_xml(config_file, pt);
185 cout <<
"using config file: " << config_filename << endl;
189 cout <<
"no xml config file - using internal values" << endl;
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);
203 N = pt.get(
"HIJING.N", 1);
205 keepSpectators = pt.get(
"HIJING.KEEP_SPECTATORS", 1);
208 output = pt.get(
"HIJING.OUTPUT",
"sHijing.dat");
212 std::random_device rdev;
213 randomSeed = pt.get(
"HIJING.RANDOM.SEED", rdev());
215 engine =
new CLHEP::MTwistEngine(randomSeed);
218 iptree &pr1 = pt.get_child(
"HIJING.HIPR1", null);
221 int key = boost::lexical_cast<
int>(v.first.data());
222 float value = boost::lexical_cast<
float>(v.second.data());
225 iptree &pr2 = pt.get_child(
"HIJING.IHPR2", null);
228 int key = boost::lexical_cast<
int>(v.first.data());
229 int value = boost::lexical_cast<
int>(v.second.data());
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;
241 iptree &
it = pt.get_child(
"HIJING.FASTJET", null);
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));
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));
253 fastjet_enable_p = 1;
255 cout <<
"seed: " << randomSeed << endl;
256 cout <<
"output: " << output << endl;
257 cout <<
"Number of Events: " << N << endl;
259 hijfst_control(fastjet_enable_p, algorithm_v, R_v, PID_v, EtaMin_v, EtaMax_v, EtMin_v);
262 m_frame.copy(frame, m_frame.size());
263 m_proj.copy(proj, m_proj.size());
264 m_targ.copy(targ, m_targ.size());
266 HIJSET(efrm, frame, proj, targ, iap, izp, iat, izt);
273 HIJING(frame, bmin, bmax);
281 HepMC::GenEvent *evt =
new HepMC::GenEvent();
282 evt->use_units(HepMC::Units::GEV, HepMC::Units::MM);
283 evt->set_event_number(events);
290 }
while (events < N);
309 int ncoll = n0 + n01 + n10 + n11;
316 std::vector<float>
x;
317 std::vector<float>
y;
319 float bbx = b * cos(bphi);
320 float bby = b * sin(bphi);
351 float xa = std::accumulate(x.begin(), x.end(), 0) / N;
352 float ya = std::accumulate(y.begin(), y.end(), 0) / N;
357 for (
int i = 0;
i <
N;
i++)
359 sxx += (x[
i] - xa) * (x[
i] - xa);
360 syy += (y[
i] - ya) * (y[
i] - ya);
361 sxy += (x[
i] - xa) * (y[
i] - ya);
366 e_part = std::sqrt((sxx - syy) * (sxx - syy) + 4 * sxy * sxy) / (sxx + syy);
370 HepMC::HeavyIon hi(1, np, nt, ncoll, 0, 0, n01, n10, n11, b, bphi, e_part, 42.0);
372 evt->set_heavy_ion(hi);
385 std::vector<int> partOriginVertex_vec;
386 std::vector<int> partDecayVertex_vec;
387 std::vector<HepMC::GenParticle *> particleHepPartPtr_vec;
389 partOriginVertex_vec.assign(numHijingPart, 0);
390 partDecayVertex_vec.assign(numHijingPart, -1);
391 particleHepPartPtr_vec.assign(numHijingPart, (HepMC::GenParticle *) 0);
394 std::vector<HepMC::GenVertex *> vertexPtrVec;
397 CLHEP::HepLorentzVector newVertex;
398 newVertex = CLHEP::HepLorentzVector(0., 0., 0., 0.);
399 HepMC::GenVertex *
v1 =
new HepMC::GenVertex(newVertex);
402 vertexPtrVec.push_back(v1);
405 int proj_id = (
m_proj ==
"A") ? 3000000 +
iap : 2212;
406 v1->add_particle_in(
new HepMC::GenParticle(CLHEP::HepLorentzVector(0., 0., eproj, eproj),
410 int targ_id = (
m_targ ==
"A") ? 3000000 +
iap : 2212;
411 v1->add_particle_in(
new HepMC::GenParticle(CLHEP::HepLorentzVector(0., 0., -etarg, etarg),
415 for (
int i = 1;
i <= natt; ++
i)
424 v1->add_particle_out(
new HepMC::GenParticle(jetP4,
438 int parentOriginIndex = 0;
439 int parentDecayIndex = -1;
442 if (parentIndex >= 0)
444 parentOriginIndex = partOriginVertex_vec[parentIndex];
445 parentDecayIndex = partDecayVertex_vec[parentIndex];
455 int particleVertexIndex = 0;
462 if (parentDecayIndex != -1)
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());
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();
477 std::cout << (*iter)->barcode() <<
", ";
480 std::cout << std::endl;
484 particleVertexIndex = parentDecayIndex;
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());
503 <<
"HIJING BUG:: Particle found with displaced vertex but no parent "
509 HepMC::GenVertex *newVertex_p =
new HepMC::GenVertex(particleStart);
511 evt->add_vertex(newVertex_p);
512 vertexPtrVec.push_back(newVertex_p);
513 particleVertexIndex = vertexPtrVec.size() - 1;
516 partDecayVertex_vec[parentIndex] = particleVertexIndex;
519 newVertex_p->add_particle_in(particleHepPartPtr_vec[parentIndex]);
524 particleVertexIndex = parentOriginIndex;
532 for (
unsigned int ivert = 0; ivert < vertexPtrVec.size(); ivert++)
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());
548 HepMC::GenVertex *newVertex_p =
new HepMC::GenVertex(particleStart);
550 evt->add_vertex(newVertex_p);
551 vertexPtrVec.push_back(newVertex_p);
552 particleVertexIndex = vertexPtrVec.size() - 1;
556 particleVertexIndex = foundVert;
562 int particleStatus = 1;
574 HepMC::GenParticle *newParticle_p =
new HepMC::GenParticle(particleP4,
581 particleHepPartPtr_vec[
i - 1] = newParticle_p;
584 vertexPtrVec[particleVertexIndex]->add_particle_out(newParticle_p);
585 partOriginVertex_vec[
i - 1] = particleVertexIndex;