202{
203
204
207
208
210
211
213
214
216
217
224
228
230
231 #ifdef HEPMC3
232 HepMC::GenHeavyIonPtr ion= std::make_shared<HepMC::GenHeavyIon>();
233 ion->Ncoll_hard=static_cast<int>(jatt);
234 ion->Npart_proj=
static_cast<int>(
np);
235 ion->Npart_targ=
static_cast<int>(
nt);
236 ion->Ncoll=static_cast<int>(n0+n10+n01+n11);
237 ion->N_Nwounded_collisions=static_cast<int>(n01);
238 ion->Nwounded_N_collisions=static_cast<int>(n10);
239 ion->Nwounded_Nwounded_collisions=static_cast<int>(n11);
240 ion->spectator_neutrons=-1;
241 ion->spectator_protons=-1;
242 ion->impact_parameter=
b;
243 ion->event_plane_angle=bphi;
244 ion->event_plane_angle=-1;
245 ion->sigma_inel_NN=sigmainel;
246 evt->set_heavy_ion(std::move(ion));
247#else
248 HepMC::HeavyIon ion
249 (
250 static_cast<int>(jatt),
251 static_cast<int>(np),
252 static_cast<int>(nt),
253 static_cast<int>(n0+n10+n01+n11),
254 static_cast<int>(-1),
255 static_cast<int>(-1),
256 static_cast<int>(n01),
257 static_cast<int>(n10),
258 static_cast<int>(n11),
259 b,
260 bphi,
261 -1,
262 sigmainel );
263
264 evt->set_heavy_ion(std::move(ion));
265 std::cout <<
" heavy ion " <<
evt->heavy_ion() << std::endl;
266#endif
267
268
269
270 bool keptHistory = (
m_hiparnt.ihpr2(21) == 1);
271
272
273
275
276
277
278 std::vector<int> partOriginVertex_vec(numHijingPart, 0);
279 std::vector<int> partDecayVertex_vec(numHijingPart, -1);
280 std::vector<HepMC::GenParticlePtr> particleHepPartPtr_vec(numHijingPart, nullptr);
281
282
283
284 std::vector<HepMC::GenVertexPtr> vertexPtrVec;
285
286
287
288 CLHEP::HepLorentzVector newVertex;
289 newVertex = CLHEP::HepLorentzVector(0.,0.,0.,0.);
290
292 else if(
m_sel) newVertex = CLHEP::HepLorentzVector(
m_x,
m_y,
m_z, 0.);
293
295
297 vertexPtrVec.push_back(v1);
298
300 if (
m_frame ==
"CMS " ) eproj = eproj / 2.;
301
302 int proj_id = 2212;
303 if (
m_proj ==
"PBAR " ) {
304 proj_id = -2212;
305 }
else if (
m_proj ==
"N " ) {
306 proj_id = 2112;
307 }
else if (
m_proj ==
"NBAR " ) {
308 proj_id = -2112;
309 }
else if (
m_proj ==
"PI+ " ) {
310 proj_id = 211;
311 }
else if (
m_proj ==
"PI- " ) {
312 proj_id = -211;
313 }
else if (
m_proj ==
"A " ) {
314 proj_id = 3000000 +
m_iap;
315 }
317 v1->add_particle_in( part_p );
318
319 double etarg = 0.;
321
322 int targ_id = 2212;
323 if (
m_targ ==
"PBAR " ) {
324 targ_id = -2212;
325 }
else if (
m_targ ==
"N " ) {
326 targ_id = 2112;
327 }
else if (
m_targ ==
"NBAR " ) {
328 targ_id = -2112;
329 }
else if (
m_targ ==
"PI+ " ) {
330 targ_id = 211;
331 }
else if (
m_targ ==
"PI- " ) {
332 targ_id = -211;
333 }
else if (
m_targ ==
"A " ) {
334 targ_id = 3000000 +
m_iat;
335 }
337 v1->add_particle_in( part_t );
338
339 evt->set_beam_particles(std::move(part_p),std::move(part_t));
340
342 << " px, "
343 << " py, "
344 << " pz, "
345 << " Id, "
346 << " Source, "
347 << " Parent, "
348 << " Status, " );
349
350 bool inconsistency = false;
351
352
353
355 {
356
357
359
360
361
362 int parentIndex =
m_himain2.katt(i, 3) - 1;
363 int parentOriginIndex = 0;
364 int parentDecayIndex = -1;
365
366
367
368 if (parentIndex >= 0)
369 {
370 parentOriginIndex = partOriginVertex_vec[parentIndex];
371 parentDecayIndex = partDecayVertex_vec[parentIndex];
372 }
373
374
375
376
377
383 }
384
386 << std::setw(5) << i << ","
387 << std::setw(7) <<
m_himain2.patt(i, 1) <<
", "
388 << std::setw(7) <<
m_himain2.patt(i, 2) <<
", "
389 << std::setw(7) <<
m_himain2.patt(i, 3) <<
", "
390 << std::setw(7) <<
m_himain2.katt(i, 1) <<
", "
391 << std::setw(7) <<
m_himain2.katt(i, 2) <<
", "
392 << std::setw(7) <<
m_himain2.katt(i, 3) <<
", "
393 << std::setw(7) <<
m_himain2.katt(i, 4) <<
", "
394 << std::setw(7) <<
m_himain2.vatt(i, 1) <<
", "
395 << std::setw(7) <<
m_himain2.vatt(i, 2) <<
", "
396 << std::setw(7) <<
m_himain2.vatt(i, 3) );
397
400
401
402
403 int particleVertexIndex = 0;
404
405
406
407 if (keptHistory)
408 {
409
410
411
412 if (parentDecayIndex != -1)
413 {
414
415
416 HepGeom::Point3D<double> vertex_pos(vertexPtrVec[parentDecayIndex]->position().
x(),
417 vertexPtrVec[parentDecayIndex]->position().
y(),
418 vertexPtrVec[parentDecayIndex]->position().
z());
419 double distance = vertex_pos.distance(particleStart.vect());
420
421
423 {
424
425
426 ATH_MSG_WARNING(
" Inconsistency in Hijing particle vertexing, particle # " << i
427 << " starting point (x,y,z) = ("
428 << particleStart.x() << ", "
429 << particleStart.y() << ", "
430 << particleStart.z() << ") "
431 << " a distance " << distance << " away from parent decay point " );
432
433
434
435
436 log << MSG::WARNING <<
" Parent decay vertex: (x,y,z) = " << vertexPtrVec[parentDecayIndex]->position().x()
437 << ", " << vertexPtrVec[parentDecayIndex]->position().y()
438 << ", " << vertexPtrVec[parentDecayIndex]->position().z()
439 << ", associated daughter IDs = ";
440
441#ifdef HEPMC3
442 auto vertexPtrVec_particles_out_const_begin=vertexPtrVec[parentDecayIndex]->particles_out().begin();
443 auto vertexPtrVec_particles_out_const_end=vertexPtrVec[parentDecayIndex]->particles_out().end();
444#else
445 auto vertexPtrVec_particles_out_const_begin=vertexPtrVec[parentDecayIndex]->particles_out_const_begin();
446 auto vertexPtrVec_particles_out_const_end=vertexPtrVec[parentDecayIndex]->particles_out_const_end();
447#endif
448 for (auto iter = vertexPtrVec_particles_out_const_begin;
449 iter != vertexPtrVec_particles_out_const_end;
451 {
452 log << (*iter) <<
", ";
453 }
454
456 inconsistency = true;
457 }
458
459
460
461 particleVertexIndex = parentDecayIndex;
462 }
463 else
464 {
465
466
467
468 HepGeom::Point3D<double> vertex_pos(vertexPtrVec[parentOriginIndex]->position().
x(),
469 vertexPtrVec[parentOriginIndex]->position().
y(),
470 vertexPtrVec[parentOriginIndex]->position().
z());
471 double distance = vertex_pos.distance(particleStart.vect());
472
474 {
475
476
477
478
479 ATH_MSG_WARNING(
"HIJING BUG:: Particle found with displaced vertex but no parent " );
480 ATH_MSG_WARNING(
" Particle parameters: " << std::fixed << std::setprecision(2)
481 << std::setw(5) << i << ","
482 << std::setw(7) <<
m_himain2.patt(i, 1) <<
", "
483 << std::setw(7) <<
m_himain2.patt(i, 2) <<
", "
484 << std::setw(7) <<
m_himain2.patt(i, 3) <<
", "
485 << std::setw(7) <<
m_himain2.katt(i, 1) <<
", "
486 << std::setw(7) <<
m_himain2.katt(i, 2) <<
", "
487 << std::setw(7) <<
m_himain2.katt(i, 3) <<
", "
488 << std::setw(7) <<
m_himain2.katt(i, 4) <<
", "
489 << std::setw(7) <<
m_himain2.vatt(i, 1) <<
", "
490 << std::setw(7) <<
m_himain2.vatt(i, 2) <<
", "
491 << std::setw(7) <<
m_himain2.vatt(i, 3) );
492
493
494
495 particleVertexIndex = parentOriginIndex;
496
497 }
499 {
500
501
503 evt->add_vertex(newVertex_p);
504 vertexPtrVec.push_back(newVertex_p);
505 particleVertexIndex = vertexPtrVec.size() - 1;
506
507
508
509 partDecayVertex_vec[parentIndex] = particleVertexIndex;
510
511
512
513
514 newVertex_p->add_particle_in(particleHepPartPtr_vec[parentIndex]);
515 }
516 else {
517
518
519 particleVertexIndex = parentOriginIndex;
520 }
521 }
522 }
523 else
524 {
525
526
527 int foundVert = -1;
528 for (unsigned int ivert = 0; ivert < vertexPtrVec.size(); ivert++)
529 {
530
531 HepGeom::Point3D<double> vertex_pos(vertexPtrVec[ivert]->position().
x(),
532 vertexPtrVec[ivert]->position().
y(),
533 vertexPtrVec[ivert]->position().
z());
534 double distance = vertex_pos.distance(particleStart.vect());
536 {
537 foundVert = ivert;
538 break;
539 }
540 }
541
542 if (foundVert == -1)
543 {
544
545
547
548 evt->add_vertex(newVertex_p);
549 vertexPtrVec.push_back(std::move(newVertex_p));
550 particleVertexIndex = vertexPtrVec.size() - 1;
551 }
552 else
553 {
554 particleVertexIndex = foundVert;
555 }
556
557 }
558
559
560
562 int particleStatus = 1;
563 if (
m_himain2.katt(i, 4) == 11) particleStatus = 2;
564
565
566
567 HepMC::FourVector particleP4(
m_himain2.patt(i, 1),
571
573 particleStatus);
574
575
576
577
578
579 particleHepPartPtr_vec[
i-1] = newParticle_p;
580
581
582
583 vertexPtrVec[particleVertexIndex]->add_particle_out(std::move(newParticle_p));
584 partOriginVertex_vec[
i-1] = particleVertexIndex;
585
586 }
587
588 if (inconsistency)
589 {
591 {
592
593
595
596 ATH_MSG_WARNING(
" Hijing.cxx inconsistency: " << std::fixed << std::setprecision(2)
597
598 << std::setw(5) << i << ","
599 << std::setw(7) <<
m_himain2.patt(i, 1) <<
", "
600 << std::setw(7) <<
m_himain2.patt(i, 2) <<
", "
601 << std::setw(7) <<
m_himain2.patt(i, 3) <<
", "
602 << std::setw(7) <<
m_himain2.katt(i, 1) <<
", "
603 << std::setw(7) <<
m_himain2.katt(i, 2) <<
", "
604 << std::setw(7) <<
m_himain2.katt(i, 3) <<
", "
605 << std::setw(7) <<
m_himain2.katt(i, 4) <<
", "
606 << std::setw(7) <<
m_himain2.vatt(i, 1) <<
", "
607 << std::setw(7) <<
m_himain2.vatt(i, 2) <<
", "
608 << std::setw(7) <<
m_himain2.vatt(i, 3) );
609 }
610 }
611
612
613
614 for (
int isg = 1; isg <=
m_hijjet2.nsg(); isg++)
615 {
616 for (
int jparton = 1; jparton <=
m_hijjet2.njsg(isg); jparton++)
617 {
622
624 double e = std::sqrt(ptsq + pz*pz + mass*mass);
625 double mt = std::sqrt(ptsq + mass*mass);
626 double pt = std::sqrt(ptsq);
627
628 double pseud = 0.5*std::log((e + pz)/(e - pz));
629
630 int partonId =
m_hijjet2.k2sg(isg, jparton);
631
633 {
635 << ", pt = " << pt
636 << ", mt = " << mt
637 << ", eta = " << pseud );
638
639
640
642 }
643 }
644 }
645
646
647
650
651 for (int iproj = 1; iproj <= iap; iproj++)
652 {
653 for (
int jparton = 1; jparton <=
m_hijjet1.npj(iproj); jparton++)
654 {
659
661 double e = std::sqrt(ptsq + pz*pz + mass*mass);
662 double mt = std::sqrt(ptsq + mass*mass);
663 double pt = std::sqrt(ptsq);
664
665 double pseud = 0.5*std::log((e + pz)/(e - pz));
666
667 int partonId =
m_hijjet1.kfpj(iproj, jparton);
668
670 {
672 << ", pt = " << pt
673 << ", mt = " << mt
674 << ", eta = " << pseud );
675
676
677
679 }
680 }
681 }
682
683
684
685 for (int itarg = 1; itarg <= iat; itarg++)
686 {
687 for (
int jparton = 1; jparton <=
m_hijjet1.ntj(itarg); jparton++)
688 {
693
695 double e = std::sqrt(ptsq + pz*pz + mass*mass);
696 double mt = std::sqrt(ptsq + mass*mass);
697 double pt = std::sqrt(ptsq);
698
699 double pseud = 0.5*std::log((e + pz)/(e - pz));
700
701 int partonId =
m_hijjet1.kftj(itarg, jparton);
702
704 {
706 << ", pt = " << pt
707 << ", mt = " << mt
708 << ", eta = " << pseud );
709
710
711
713 }
714 }
715 }
716
717#ifdef HEPMC3
718
719 evt->set_units(HepMC3::Units::MEV, HepMC3::Units::MM);
720#endif
721
722
724 HepMC::FourVector tmpmom(0.,0.,0.,0.);
725 double ranz = CLHEP::RandFlat::shoot(p_Engine);
726
727 if (ranz < 0.5) {
728
729 for(auto pitr: *evt){
730 tmpmom= pitr->momentum();
731 tmpmom.setX(-tmpmom.x());
732 tmpmom.setY(-tmpmom.y());
733 tmpmom.setZ(-tmpmom.z());
734 tmpmom.setT(tmpmom.t());
735 pitr->set_momentum(tmpmom);
736 }
737 }
738 }
739
740
741 return StatusCode::SUCCESS;
742}
#define ATH_MSG_WARNING(x)
virtual CLHEP::HepLorentzVector randomizeVertex(CLHEP::HepRandomEngine *engine)
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
void set_signal_process_vertex(GenEvent *e, T v)
void set_random_states(GenEvent *e, std::vector< T > a)
HepMC::GenVertex * GenVertexPtr
GenVertexPtr newGenVertexPtr(const HepMC::FourVector &pos=HepMC::FourVector(0.0, 0.0, 0.0, 0.0), const int i=0)
GenParticlePtr newGenParticlePtr(const HepMC::FourVector &mom=HepMC::FourVector(0.0, 0.0, 0.0, 0.0), int pid=0, int status=0)
void set_signal_process_id(GenEvent *e, const int i)
GenParticle * GenParticlePtr
msgSvc
Provide convenience handles for various services.