244{
245
246
249
250
252
253
255
256
258
259
266
270
272
274 ion->Ncoll_hard=static_cast<int>(jatt);
275 ion->Npart_proj=
static_cast<int>(
np);
276 ion->Npart_targ=
static_cast<int>(
nt);
277 ion->Ncoll=static_cast<int>(n0+n10+n01+n11);
278 ion->N_Nwounded_collisions=static_cast<int>(n01);
279 ion->Nwounded_N_collisions=static_cast<int>(n10);
280 ion->Nwounded_Nwounded_collisions=static_cast<int>(n11);
281 ion->spectator_neutrons=-1;
282 ion->spectator_protons=-1;
283 ion->impact_parameter=
b;
284 ion->event_plane_angle=bphi;
285 ion->event_plane_angle=-1;
286 ion->sigma_inel_NN=sigmainel;
287 evt->set_heavy_ion(std::move(ion));
288
289
290
291 bool keptHistory = (
m_hiparnt.ihpr2(21) == 1);
292
293
294
296
297
298
299 std::vector<int> partOriginVertex_vec(numHijingPart, 0);
300 std::vector<int> partDecayVertex_vec(numHijingPart, -1);
301 std::vector<HepMC::GenParticlePtr> particleHepPartPtr_vec(numHijingPart, nullptr);
302
303
304
305 std::vector<HepMC::GenVertexPtr> vertexPtrVec;
306
307
308
309 CLHEP::HepLorentzVector newVertex;
310 newVertex = CLHEP::HepLorentzVector(0.,0.,0.,0.);
311
313 else if(
m_sel) newVertex = CLHEP::HepLorentzVector(
m_x,
m_y,
m_z, 0.);
314
316
318 vertexPtrVec.push_back(v1);
319
321 if (
m_frame ==
"CMS " ) eproj = eproj / 2.;
322
323 int proj_id = 2212;
324 if (
m_proj ==
"PBAR " ) {
325 proj_id = -2212;
326 }
else if (
m_proj ==
"N " ) {
327 proj_id = 2112;
328 }
else if (
m_proj ==
"NBAR " ) {
329 proj_id = -2112;
330 }
else if (
m_proj ==
"PI+ " ) {
331 proj_id = 211;
332 }
else if (
m_proj ==
"PI- " ) {
333 proj_id = -211;
334 }
else if (
m_proj ==
"A " ) {
335 proj_id = 3000000 +
m_iap;
336 }
338 v1->add_particle_in( part_p );
339
340 double etarg = 0.;
342
343 int targ_id = 2212;
344 if (
m_targ ==
"PBAR " ) {
345 targ_id = -2212;
346 }
else if (
m_targ ==
"N " ) {
347 targ_id = 2112;
348 }
else if (
m_targ ==
"NBAR " ) {
349 targ_id = -2112;
350 }
else if (
m_targ ==
"PI+ " ) {
351 targ_id = 211;
352 }
else if (
m_targ ==
"PI- " ) {
353 targ_id = -211;
354 }
else if (
m_targ ==
"A " ) {
355 targ_id = 3000000 +
m_iat;
356 }
358 v1->add_particle_in( part_t );
359
360 evt->set_beam_particles(std::move(part_p),std::move(part_t));
361
363 << " px, "
364 << " py, "
365 << " pz, "
366 << " Id, "
367 << " Source, "
368 << " Parent, "
369 << " Status, " );
370
371 bool inconsistency = false;
372
373
374
376 {
377
378
380
381
382
383 int parentIndex =
m_himain2.katt(i, 3) - 1;
384 int parentOriginIndex = 0;
385 int parentDecayIndex = -1;
386
387
388
389 if (parentIndex >= 0)
390 {
391 parentOriginIndex = partOriginVertex_vec[parentIndex];
392 parentDecayIndex = partDecayVertex_vec[parentIndex];
393 }
394
395
396
397
398
404 }
405
407 << std::setw(5) << i << ","
408 << std::setw(7) <<
m_himain2.patt(i, 1) <<
", "
409 << std::setw(7) <<
m_himain2.patt(i, 2) <<
", "
410 << std::setw(7) <<
m_himain2.patt(i, 3) <<
", "
411 << std::setw(7) <<
m_himain2.katt(i, 1) <<
", "
412 << std::setw(7) <<
m_himain2.katt(i, 2) <<
", "
413 << std::setw(7) <<
m_himain2.katt(i, 3) <<
", "
414 << std::setw(7) <<
m_himain2.katt(i, 4) <<
", "
415 << std::setw(7) <<
m_himain2.vatt(i, 1) <<
", "
416 << std::setw(7) <<
m_himain2.vatt(i, 2) <<
", "
417 << std::setw(7) <<
m_himain2.vatt(i, 3) );
418
421
422
423
424 int particleVertexIndex = 0;
425
426
427
428 if (keptHistory)
429 {
430
431
432
433 if (parentDecayIndex != -1)
434 {
435
436
437 HepGeom::Point3D<double> vertex_pos(vertexPtrVec[parentDecayIndex]->
position().
x(),
438 vertexPtrVec[parentDecayIndex]->
position().
y(),
439 vertexPtrVec[parentDecayIndex]->
position().
z());
440 double distance = vertex_pos.distance(particleStart.vect());
441
442
444 {
445
446
447 ATH_MSG_WARNING(
" Inconsistency in Hijing particle vertexing, particle # " << i
448 << " starting point (x,y,z) = ("
449 << particleStart.x() << ", "
450 << particleStart.y() << ", "
451 << particleStart.z() << ") "
452 << " a distance " << distance << " away from parent decay point " );
453
454
455
456
457 log << MSG::WARNING <<
" Parent decay vertex: (x,y,z) = " << vertexPtrVec[parentDecayIndex]->position().x()
458 << ", " << vertexPtrVec[parentDecayIndex]->position().y()
459 << ", " << vertexPtrVec[parentDecayIndex]->position().z()
460 << ", associated daughter IDs = ";
461
462 auto vertexPtrVec_particles_out_const_begin=vertexPtrVec[parentDecayIndex]->particles_out().begin();
463 auto vertexPtrVec_particles_out_const_end=vertexPtrVec[parentDecayIndex]->particles_out().end();
464 for (auto iter = vertexPtrVec_particles_out_const_begin;
465 iter != vertexPtrVec_particles_out_const_end;
467 {
468 log << (*iter) <<
", ";
469 }
470
472 inconsistency = true;
473 }
474
475
476
477 particleVertexIndex = parentDecayIndex;
478 }
479 else
480 {
481
482
483
484 HepGeom::Point3D<double> vertex_pos(vertexPtrVec[parentOriginIndex]->
position().
x(),
485 vertexPtrVec[parentOriginIndex]->
position().
y(),
486 vertexPtrVec[parentOriginIndex]->
position().
z());
487 double distance = vertex_pos.distance(particleStart.vect());
488
490 {
491
492
493
494
495 ATH_MSG_WARNING(
"HIJING BUG:: Particle found with displaced vertex but no parent " );
496 ATH_MSG_WARNING(
" Particle parameters: " << std::fixed << std::setprecision(2)
497 << std::setw(5) << i << ","
498 << std::setw(7) <<
m_himain2.patt(i, 1) <<
", "
499 << std::setw(7) <<
m_himain2.patt(i, 2) <<
", "
500 << std::setw(7) <<
m_himain2.patt(i, 3) <<
", "
501 << std::setw(7) <<
m_himain2.katt(i, 1) <<
", "
502 << std::setw(7) <<
m_himain2.katt(i, 2) <<
", "
503 << std::setw(7) <<
m_himain2.katt(i, 3) <<
", "
504 << std::setw(7) <<
m_himain2.katt(i, 4) <<
", "
505 << std::setw(7) <<
m_himain2.vatt(i, 1) <<
", "
506 << std::setw(7) <<
m_himain2.vatt(i, 2) <<
", "
507 << std::setw(7) <<
m_himain2.vatt(i, 3) );
508
509
510
511 particleVertexIndex = parentOriginIndex;
512
513 }
515 {
516
517
519 evt->add_vertex(newVertex_p);
520 vertexPtrVec.push_back(newVertex_p);
521 particleVertexIndex = vertexPtrVec.size() - 1;
522
523
524
525 partDecayVertex_vec[parentIndex] = particleVertexIndex;
526
527
528
529
530 newVertex_p->add_particle_in(particleHepPartPtr_vec[parentIndex]);
531 }
532 else {
533
534
535 particleVertexIndex = parentOriginIndex;
536 }
537 }
538 }
539 else
540 {
541
542
543 int foundVert = -1;
544 for (unsigned int ivert = 0; ivert < vertexPtrVec.size(); ivert++)
545 {
546
547 HepGeom::Point3D<double> vertex_pos(vertexPtrVec[ivert]->
position().
x(),
550 double distance = vertex_pos.distance(particleStart.vect());
552 {
553 foundVert = ivert;
554 break;
555 }
556 }
557
558 if (foundVert == -1)
559 {
560
561
563
564 evt->add_vertex(newVertex_p);
565 vertexPtrVec.push_back(std::move(newVertex_p));
566 particleVertexIndex = vertexPtrVec.size() - 1;
567 }
568 else
569 {
570 particleVertexIndex = foundVert;
571 }
572
573 }
574
575
576
578 int particleStatus = 1;
579 if (
m_himain2.katt(i, 4) == 11) particleStatus = 2;
580
581
582
587
589 particleStatus);
590
591
592
593
594
595 particleHepPartPtr_vec[
i-1] = newParticle_p;
596
597
598
599 vertexPtrVec[particleVertexIndex]->add_particle_out(std::move(newParticle_p));
600 partOriginVertex_vec[
i-1] = particleVertexIndex;
601
602 }
603
604 if (inconsistency)
605 {
607 {
608
609
611
612 ATH_MSG_WARNING(
" Hijing.cxx inconsistency: " << std::fixed << std::setprecision(2)
613
614 << std::setw(5) << i << ","
615 << std::setw(7) <<
m_himain2.patt(i, 1) <<
", "
616 << std::setw(7) <<
m_himain2.patt(i, 2) <<
", "
617 << std::setw(7) <<
m_himain2.patt(i, 3) <<
", "
618 << std::setw(7) <<
m_himain2.katt(i, 1) <<
", "
619 << std::setw(7) <<
m_himain2.katt(i, 2) <<
", "
620 << std::setw(7) <<
m_himain2.katt(i, 3) <<
", "
621 << std::setw(7) <<
m_himain2.katt(i, 4) <<
", "
622 << std::setw(7) <<
m_himain2.vatt(i, 1) <<
", "
623 << std::setw(7) <<
m_himain2.vatt(i, 2) <<
", "
624 << std::setw(7) <<
m_himain2.vatt(i, 3) );
625 }
626 }
627
628
629
630 for (
int isg = 1; isg <=
m_hijjet2.nsg(); isg++)
631 {
632 for (
int jparton = 1; jparton <=
m_hijjet2.njsg(isg); jparton++)
633 {
638
640 double e = std::sqrt(ptsq + pz*pz + mass*mass);
641 double mt = std::sqrt(ptsq + mass*mass);
642 double pt = std::sqrt(ptsq);
643
644 double pseud = 0.5*std::log((e + pz)/(e - pz));
645
646 int partonId =
m_hijjet2.k2sg(isg, jparton);
647
649 {
651 << ", pt = " << pt
652 << ", mt = " << mt
653 << ", eta = " << pseud );
654
655
656
658 }
659 }
660 }
661
662
663
666
667 for (int iproj = 1; iproj <= iap; iproj++)
668 {
669 for (
int jparton = 1; jparton <=
m_hijjet1.npj(iproj); jparton++)
670 {
675
677 double e = std::sqrt(ptsq + pz*pz + mass*mass);
678 double mt = std::sqrt(ptsq + mass*mass);
679 double pt = std::sqrt(ptsq);
680
681 double pseud = 0.5*std::log((e + pz)/(e - pz));
682
683 int partonId =
m_hijjet1.kfpj(iproj, jparton);
684
686 {
688 << ", pt = " << pt
689 << ", mt = " << mt
690 << ", eta = " << pseud );
691
692
693
695 }
696 }
697 }
698
699
700
701 for (int itarg = 1; itarg <= iat; itarg++)
702 {
703 for (
int jparton = 1; jparton <=
m_hijjet1.ntj(itarg); jparton++)
704 {
709
711 double e = std::sqrt(ptsq + pz*pz + mass*mass);
712 double mt = std::sqrt(ptsq + mass*mass);
713 double pt = std::sqrt(ptsq);
714
715 double pseud = 0.5*std::log((e + pz)/(e - pz));
716
717 int partonId =
m_hijjet1.kftj(itarg, jparton);
718
720 {
722 << ", pt = " << pt
723 << ", mt = " << mt
724 << ", eta = " << pseud );
725
726
727
729 }
730 }
731 }
732
733
734 evt->set_units(HepMC3::Units::MEV, HepMC3::Units::MM);
735
736
739 double ranz = CLHEP::RandFlat::shoot(p_Engine);
740
741 if (ranz < 0.5) {
742
743 for(auto pitr: *evt){
744 tmpmom= pitr->momentum();
745 tmpmom.setX(-tmpmom.x());
746 tmpmom.setY(-tmpmom.y());
747 tmpmom.setZ(-tmpmom.z());
748 tmpmom.setT(tmpmom.t());
749 pitr->set_momentum(tmpmom);
750 }
751 }
752 }
753
754
755 return StatusCode::SUCCESS;
756}
#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
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.