232{
233
234
237
238
240
241
243
244
246
247
254
258
260
262 ion->Ncoll_hard=static_cast<int>(jatt);
263 ion->Npart_proj=
static_cast<int>(
np);
264 ion->Npart_targ=
static_cast<int>(
nt);
265 ion->Ncoll=static_cast<int>(n0+n10+n01+n11);
266 ion->N_Nwounded_collisions=static_cast<int>(n01);
267 ion->Nwounded_N_collisions=static_cast<int>(n10);
268 ion->Nwounded_Nwounded_collisions=static_cast<int>(n11);
269 ion->spectator_neutrons=-1;
270 ion->spectator_protons=-1;
271 ion->impact_parameter=
b;
272 ion->event_plane_angle=bphi;
273 ion->event_plane_angle=-1;
274 ion->sigma_inel_NN=sigmainel;
275 evt->set_heavy_ion(std::move(ion));
276
277
278
279 bool keptHistory = (
m_hiparnt.ihpr2(21) == 1);
280
281
282
284
285
286
287 std::vector<int> partOriginVertex_vec(numHijingPart, 0);
288 std::vector<int> partDecayVertex_vec(numHijingPart, -1);
289 std::vector<HepMC::GenParticlePtr> particleHepPartPtr_vec(numHijingPart, nullptr);
290
291
292
293 std::vector<HepMC::GenVertexPtr> vertexPtrVec;
294
295
296
297 CLHEP::HepLorentzVector newVertex;
298 newVertex = CLHEP::HepLorentzVector(0.,0.,0.,0.);
299
301 else if(
m_sel) newVertex = CLHEP::HepLorentzVector(
m_x,
m_y,
m_z, 0.);
302
304
306 vertexPtrVec.push_back(v1);
307
309 if (
m_frame ==
"CMS " ) eproj = eproj / 2.;
310
311 int proj_id = 2212;
312 if (
m_proj ==
"PBAR " ) {
313 proj_id = -2212;
314 }
else if (
m_proj ==
"N " ) {
315 proj_id = 2112;
316 }
else if (
m_proj ==
"NBAR " ) {
317 proj_id = -2112;
318 }
else if (
m_proj ==
"PI+ " ) {
319 proj_id = 211;
320 }
else if (
m_proj ==
"PI- " ) {
321 proj_id = -211;
322 }
else if (
m_proj ==
"A " ) {
323 proj_id = 3000000 +
m_iap;
324 }
326 v1->add_particle_in( part_p );
327
328 double etarg = 0.;
330
331 int targ_id = 2212;
332 if (
m_targ ==
"PBAR " ) {
333 targ_id = -2212;
334 }
else if (
m_targ ==
"N " ) {
335 targ_id = 2112;
336 }
else if (
m_targ ==
"NBAR " ) {
337 targ_id = -2112;
338 }
else if (
m_targ ==
"PI+ " ) {
339 targ_id = 211;
340 }
else if (
m_targ ==
"PI- " ) {
341 targ_id = -211;
342 }
else if (
m_targ ==
"A " ) {
343 targ_id = 3000000 +
m_iat;
344 }
346 v1->add_particle_in( part_t );
347
348 evt->set_beam_particles(std::move(part_p),std::move(part_t));
349
351 << " px, "
352 << " py, "
353 << " pz, "
354 << " Id, "
355 << " Source, "
356 << " Parent, "
357 << " Status, " );
358
359 bool inconsistency = false;
360
361
362
364 {
365
366
368
369
370
371 int parentIndex =
m_himain2.katt(i, 3) - 1;
372 int parentOriginIndex = 0;
373 int parentDecayIndex = -1;
374
375
376
377 if (parentIndex >= 0)
378 {
379 parentOriginIndex = partOriginVertex_vec[parentIndex];
380 parentDecayIndex = partDecayVertex_vec[parentIndex];
381 }
382
383
384
385
386
392 }
393
395 << std::setw(5) << i << ","
396 << std::setw(7) <<
m_himain2.patt(i, 1) <<
", "
397 << std::setw(7) <<
m_himain2.patt(i, 2) <<
", "
398 << std::setw(7) <<
m_himain2.patt(i, 3) <<
", "
399 << std::setw(7) <<
m_himain2.katt(i, 1) <<
", "
400 << std::setw(7) <<
m_himain2.katt(i, 2) <<
", "
401 << std::setw(7) <<
m_himain2.katt(i, 3) <<
", "
402 << std::setw(7) <<
m_himain2.katt(i, 4) <<
", "
403 << std::setw(7) <<
m_himain2.vatt(i, 1) <<
", "
404 << std::setw(7) <<
m_himain2.vatt(i, 2) <<
", "
405 << std::setw(7) <<
m_himain2.vatt(i, 3) );
406
409
410
411
412 int particleVertexIndex = 0;
413
414
415
416 if (keptHistory)
417 {
418
419
420
421 if (parentDecayIndex != -1)
422 {
423
424
425 HepGeom::Point3D<double> vertex_pos(vertexPtrVec[parentDecayIndex]->
position().
x(),
426 vertexPtrVec[parentDecayIndex]->
position().
y(),
427 vertexPtrVec[parentDecayIndex]->
position().
z());
428 double distance = vertex_pos.distance(particleStart.vect());
429
430
432 {
433
434
435 ATH_MSG_WARNING(
" Inconsistency in Hijing particle vertexing, particle # " << i
436 << " starting point (x,y,z) = ("
437 << particleStart.x() << ", "
438 << particleStart.y() << ", "
439 << particleStart.z() << ") "
440 << " a distance " << distance << " away from parent decay point " );
441
442
443
444
445 log << MSG::WARNING <<
" Parent decay vertex: (x,y,z) = " << vertexPtrVec[parentDecayIndex]->position().x()
446 << ", " << vertexPtrVec[parentDecayIndex]->position().y()
447 << ", " << vertexPtrVec[parentDecayIndex]->position().z()
448 << ", associated daughter IDs = ";
449
450 auto vertexPtrVec_particles_out_const_begin=vertexPtrVec[parentDecayIndex]->particles_out().begin();
451 auto vertexPtrVec_particles_out_const_end=vertexPtrVec[parentDecayIndex]->particles_out().end();
452 for (auto iter = vertexPtrVec_particles_out_const_begin;
453 iter != vertexPtrVec_particles_out_const_end;
455 {
456 log << (*iter) <<
", ";
457 }
458
460 inconsistency = true;
461 }
462
463
464
465 particleVertexIndex = parentDecayIndex;
466 }
467 else
468 {
469
470
471
472 HepGeom::Point3D<double> vertex_pos(vertexPtrVec[parentOriginIndex]->
position().
x(),
473 vertexPtrVec[parentOriginIndex]->
position().
y(),
474 vertexPtrVec[parentOriginIndex]->
position().
z());
475 double distance = vertex_pos.distance(particleStart.vect());
476
478 {
479
480
481
482
483 ATH_MSG_WARNING(
"HIJING BUG:: Particle found with displaced vertex but no parent " );
484 ATH_MSG_WARNING(
" Particle parameters: " << std::fixed << std::setprecision(2)
485 << std::setw(5) << i << ","
486 << std::setw(7) <<
m_himain2.patt(i, 1) <<
", "
487 << std::setw(7) <<
m_himain2.patt(i, 2) <<
", "
488 << std::setw(7) <<
m_himain2.patt(i, 3) <<
", "
489 << std::setw(7) <<
m_himain2.katt(i, 1) <<
", "
490 << std::setw(7) <<
m_himain2.katt(i, 2) <<
", "
491 << std::setw(7) <<
m_himain2.katt(i, 3) <<
", "
492 << std::setw(7) <<
m_himain2.katt(i, 4) <<
", "
493 << std::setw(7) <<
m_himain2.vatt(i, 1) <<
", "
494 << std::setw(7) <<
m_himain2.vatt(i, 2) <<
", "
495 << std::setw(7) <<
m_himain2.vatt(i, 3) );
496
497
498
499 particleVertexIndex = parentOriginIndex;
500
501 }
503 {
504
505
507 evt->add_vertex(newVertex_p);
508 vertexPtrVec.push_back(newVertex_p);
509 particleVertexIndex = vertexPtrVec.size() - 1;
510
511
512
513 partDecayVertex_vec[parentIndex] = particleVertexIndex;
514
515
516
517
518 newVertex_p->add_particle_in(particleHepPartPtr_vec[parentIndex]);
519 }
520 else {
521
522
523 particleVertexIndex = parentOriginIndex;
524 }
525 }
526 }
527 else
528 {
529
530
531 int foundVert = -1;
532 for (unsigned int ivert = 0; ivert < vertexPtrVec.size(); ivert++)
533 {
534
535 HepGeom::Point3D<double> vertex_pos(vertexPtrVec[ivert]->
position().
x(),
538 double distance = vertex_pos.distance(particleStart.vect());
540 {
541 foundVert = ivert;
542 break;
543 }
544 }
545
546 if (foundVert == -1)
547 {
548
549
551
552 evt->add_vertex(newVertex_p);
553 vertexPtrVec.push_back(std::move(newVertex_p));
554 particleVertexIndex = vertexPtrVec.size() - 1;
555 }
556 else
557 {
558 particleVertexIndex = foundVert;
559 }
560
561 }
562
563
564
566 int particleStatus = 1;
567 if (
m_himain2.katt(i, 4) == 11) particleStatus = 2;
568
569
570
575
577 particleStatus);
578
579
580
581
582
583 particleHepPartPtr_vec[
i-1] = newParticle_p;
584
585
586
587 vertexPtrVec[particleVertexIndex]->add_particle_out(std::move(newParticle_p));
588 partOriginVertex_vec[
i-1] = particleVertexIndex;
589
590 }
591
592 if (inconsistency)
593 {
595 {
596
597
599
600 ATH_MSG_WARNING(
" Hijing.cxx inconsistency: " << std::fixed << std::setprecision(2)
601
602 << std::setw(5) << i << ","
603 << std::setw(7) <<
m_himain2.patt(i, 1) <<
", "
604 << std::setw(7) <<
m_himain2.patt(i, 2) <<
", "
605 << std::setw(7) <<
m_himain2.patt(i, 3) <<
", "
606 << std::setw(7) <<
m_himain2.katt(i, 1) <<
", "
607 << std::setw(7) <<
m_himain2.katt(i, 2) <<
", "
608 << std::setw(7) <<
m_himain2.katt(i, 3) <<
", "
609 << std::setw(7) <<
m_himain2.katt(i, 4) <<
", "
610 << std::setw(7) <<
m_himain2.vatt(i, 1) <<
", "
611 << std::setw(7) <<
m_himain2.vatt(i, 2) <<
", "
612 << std::setw(7) <<
m_himain2.vatt(i, 3) );
613 }
614 }
615
616
617
618 for (
int isg = 1; isg <=
m_hijjet2.nsg(); isg++)
619 {
620 for (
int jparton = 1; jparton <=
m_hijjet2.njsg(isg); jparton++)
621 {
626
628 double e = std::sqrt(ptsq + pz*pz + mass*mass);
629 double mt = std::sqrt(ptsq + mass*mass);
630 double pt = std::sqrt(ptsq);
631
632 double pseud = 0.5*std::log((e + pz)/(e - pz));
633
634 int partonId =
m_hijjet2.k2sg(isg, jparton);
635
637 {
639 << ", pt = " << pt
640 << ", mt = " << mt
641 << ", eta = " << pseud );
642
643
644
646 }
647 }
648 }
649
650
651
654
655 for (int iproj = 1; iproj <= iap; iproj++)
656 {
657 for (
int jparton = 1; jparton <=
m_hijjet1.npj(iproj); jparton++)
658 {
663
665 double e = std::sqrt(ptsq + pz*pz + mass*mass);
666 double mt = std::sqrt(ptsq + mass*mass);
667 double pt = std::sqrt(ptsq);
668
669 double pseud = 0.5*std::log((e + pz)/(e - pz));
670
671 int partonId =
m_hijjet1.kfpj(iproj, jparton);
672
674 {
676 << ", pt = " << pt
677 << ", mt = " << mt
678 << ", eta = " << pseud );
679
680
681
683 }
684 }
685 }
686
687
688
689 for (int itarg = 1; itarg <= iat; itarg++)
690 {
691 for (
int jparton = 1; jparton <=
m_hijjet1.ntj(itarg); jparton++)
692 {
697
699 double e = std::sqrt(ptsq + pz*pz + mass*mass);
700 double mt = std::sqrt(ptsq + mass*mass);
701 double pt = std::sqrt(ptsq);
702
703 double pseud = 0.5*std::log((e + pz)/(e - pz));
704
705 int partonId =
m_hijjet1.kftj(itarg, jparton);
706
708 {
710 << ", pt = " << pt
711 << ", mt = " << mt
712 << ", eta = " << pseud );
713
714
715
717 }
718 }
719 }
720
721
722 evt->set_units(HepMC3::Units::MEV, HepMC3::Units::MM);
723
724
727 double ranz = CLHEP::RandFlat::shoot(p_Engine);
728
729 if (ranz < 0.5) {
730
731 for(auto pitr: *evt){
732 tmpmom= pitr->momentum();
733 tmpmom.setX(-tmpmom.x());
734 tmpmom.setY(-tmpmom.y());
735 tmpmom.setZ(-tmpmom.z());
736 tmpmom.setT(tmpmom.t());
737 pitr->set_momentum(tmpmom);
738 }
739 }
740 }
741
742
743 return StatusCode::SUCCESS;
744}
#define ATH_MSG_WARNING(x,...)
DoubleProperty m_vertexOffsetCut
DoubleProperty m_partonStoreMinPt
BooleanProperty m_keepAllDecayVertices
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
HepMC3::FourVector FourVector
HepMC3::GenHeavyIonPtr GenHeavyIonPtr
void set_signal_process_vertex(GenEvent *e, T &v)
void set_signal_process_id(GenEvent *e, const int i=0)
void set_random_states(GenEvent *e, std::vector< long int > &a)
GenParticlePtr newGenParticlePtr(const HepMC3::FourVector &mom=HepMC3::FourVector::ZERO_VECTOR(), int pid=0, int status=0)
HepMC3::GenParticlePtr GenParticlePtr
GenVertexPtr newGenVertexPtr(const HepMC3::FourVector &pos=HepMC3::FourVector::ZERO_VECTOR(), const int i=0)
HepMC3::GenVertexPtr GenVertexPtr
const Amg::Vector3D & position() const
Method to retrieve the position of the Intersection.