22 #include "GaudiKernel/ThreadLocalContext.h"
25 static const std::set<std::string> attributes_to_ignore {
27 "flows",
"flow1",
"flow2",
"flow3",
32 "filterWeight",
"filterHT",
"filterMET",
33 "event_scale",
"alphaQCD",
"alphaQED",
"random_states",
"weights",
34 "GenCrossSection",
"GenPdfInfo",
"GenHeavyIon"};
42 m_isPileup(false),m_hepMCWeightSvc(
"HepMCWeightSvc",
"McEventCollectionCnv_p6")
47 m_isPileup(false),m_hepMCWeightSvc(
"HepMCWeightSvc",
"McEventCollectionCnv_p6")
72 const EventContext& ctx = Gaudi::Hive::currentContext();
74 msg <<
MSG::DEBUG <<
"Loading McEventCollection from persistent state..."
85 datapools.
part.prepareToAdd(nParts);
86 const unsigned int nEvts = persObj->
m_genEvents.size();
87 datapools.
evt.prepareToAdd(nEvts);
90 for ( std::vector<GenEvent_p6>::const_iterator
96 HepMC::GenEvent * genEvt(
nullptr);
98 genEvt =
new HepMC::GenEvent();
103 genEvt->add_attribute (
"barcodes", std::make_shared<HepMC::GenEventBarcodes>());
110 genEvt->add_attribute(
"signal_process_id", std::make_shared<HepMC3::IntAttribute>(persEvt.
m_signalProcessId));
112 genEvt->add_attribute(
"mpi", std::make_shared<HepMC3::IntAttribute>(persEvt.
m_mpi));
113 genEvt->add_attribute(
"event_scale", std::make_shared<HepMC3::DoubleAttribute>(persEvt.
m_eventScale));
114 genEvt->add_attribute(
"alphaQCD", std::make_shared<HepMC3::DoubleAttribute>(persEvt.
m_alphaQCD));
115 genEvt->add_attribute(
"alphaQED", std::make_shared<HepMC3::DoubleAttribute>(persEvt.
m_alphaQED));
116 genEvt->add_attribute(
"filterWeight", std::make_shared<HepMC3::DoubleAttribute>(persEvt.
m_filterWeight));
117 genEvt->add_attribute(
"filterHT", std::make_shared<HepMC3::DoubleAttribute>(persEvt.m_filterHT));
118 genEvt->add_attribute(
"filterMET", std::make_shared<HepMC3::DoubleAttribute>(persEvt.m_filterMET));
120 genEvt->add_attribute(
"random_states", std::make_shared<HepMC3::VectorLongIntAttribute>(persEvt.
m_randomStates));
122 genEvt->set_units(
static_cast<HepMC3::Units::MomentumUnit
>(persEvt.
m_momentumUnit),
123 static_cast<HepMC3::Units::LengthUnit
>(persEvt.
m_lengthUnit));
126 if(!genEvt->run_info()) {
127 HepMC3::GenRunInfoData ri_read;
129 ri_read.tool_name = std::vector<std::string>();
130 ri_read.tool_version = std::vector<std::string>();
131 ri_read.tool_description = std::vector<std::string>();
134 auto ri = std::make_shared<HepMC3::GenRunInfo>();
135 ri->read_data(ri_read);
136 genEvt->set_run_info(std::move(ri));
141 auto cs = std::make_shared<HepMC3::GenCrossSection>();
143 genEvt->set_cross_section(cs);
144 if(
static_cast<bool>(xsection[0]) )
145 cs->set_cross_section(xsection[2],xsection[1]);
147 cs->set_cross_section(-1.0, -1.0);
152 auto hi = std::make_shared<HepMC3::GenHeavyIon>();
153 const std::vector<float>& hIon = persEvt.
m_heavyIon;
156 static_cast<int>(hIon[12]),
157 static_cast<int>(hIon[11]),
158 static_cast<int>(hIon[10]),
159 static_cast<int>(hIon[9]),
160 static_cast<int>(hIon[8]),
161 static_cast<int>(hIon[7]),
162 static_cast<int>(hIon[6]),
163 static_cast<int>(hIon[5]),
164 static_cast<int>(hIon[4]),
169 genEvt->set_heavy_ion(std::move(hi));
179 pi->set(
static_cast<int>(
pdf[8]),
180 static_cast<int>(
pdf[7]),
186 static_cast<int>(
pdf[6]),
187 static_cast<int>(
pdf[5]));
188 genEvt->set_pdf_info(std::move(
pi));
198 std::map<int, HepMC::GenVertexPtr> brc_to_vertex;
202 for (
unsigned int iVtx = persEvt.
m_verticesBegin; iVtx != endVtx; ++iVtx ) {
204 brc_to_vertex[persObj->
m_genVertices[iVtx].m_barcode] = std::move(vtx);
209 if ( sigProcVtx != 0 && brc_to_vertex.count(sigProcVtx) ) {
214 for (
auto &
p : partToEndVtx) {
215 if ( brc_to_vertex.count(
p.second) ) {
216 auto decayVtx = brc_to_vertex[
p.second];
217 decayVtx->add_particle_in(
p.first );
219 msg << MSG::ERROR <<
"GenParticle points to null end vertex !!" <<
endmsg;
225 if ( beamPart1 != 0 && beamPart2 != 0 ) {
233 genEvt->m_mpi = persEvt.
m_mpi;
237 genEvt->m_signal_process_vertex = 0;
238 genEvt->m_beam_particle_1 = 0;
239 genEvt->m_beam_particle_2 = 0;
242 genEvt->m_vertex_barcodes.clear();
243 genEvt->m_particle_barcodes.clear();
244 genEvt->m_momentum_unit =
static_cast<HepMC::Units::MomentumUnit
>(persEvt.
m_momentumUnit);
245 genEvt->m_position_unit =
static_cast<HepMC::Units::LengthUnit
>(persEvt.
m_lengthUnit);
251 if( genEvt->m_cross_section )
252 delete genEvt->m_cross_section;
253 genEvt->m_cross_section = 0;
256 genEvt->m_cross_section =
new HepMC::GenCrossSection();
258 if(
static_cast<bool>(xsection[0]) )
259 genEvt->m_cross_section->set_cross_section(xsection[2],xsection[1]);
263 if(genEvt->m_heavy_ion )
264 delete genEvt->m_heavy_ion;
265 genEvt->m_heavy_ion = 0;
267 const std::vector<float>& hIon = persEvt.
m_heavyIon;
268 genEvt->m_heavy_ion =
new HepMC::HeavyIon
270 static_cast<int>(hIon[12]),
271 static_cast<int>(hIon[11]),
272 static_cast<int>(hIon[10]),
273 static_cast<int>(hIon[9]),
274 static_cast<int>(hIon[8]),
275 static_cast<int>(hIon[7]),
276 static_cast<int>(hIon[6]),
277 static_cast<int>(hIon[5]),
278 static_cast<int>(hIon[4]),
288 if(genEvt->m_pdf_info)
289 delete genEvt->m_pdf_info;
290 genEvt->m_pdf_info = 0;
293 genEvt->m_pdf_info =
new HepMC::PdfInfo
295 static_cast<int>(
pdf[8]),
296 static_cast<int>(
pdf[7]),
302 static_cast<int>(
pdf[6]),
303 static_cast<int>(
pdf[5])
318 for (
unsigned int iVtx= persEvt.
m_verticesBegin; iVtx != endVtx; ++iVtx ) {
327 if ( sigProcVtx != 0 ) {
328 genEvt->set_signal_process_vertex( genEvt->barcode_to_vertex( sigProcVtx ) );
333 p = partToEndVtx.begin(),
334 endItr = partToEndVtx.end();
339 decayVtx->add_particle_in(
p->first );
342 <<
"GenParticle points to null end vertex !!"
350 if ( beamPart1 != 0 && beamPart2 !=0 ) {
351 genEvt->set_beam_particles(genEvt->barcode_to_particle(beamPart1),
352 genEvt->barcode_to_particle(beamPart2));
360 msg <<
MSG::DEBUG <<
"Loaded McEventCollection from persistent state [OK]"
368 const EventContext& ctx = Gaudi::Hive::currentContext();
370 msg <<
MSG::DEBUG <<
"Creating persistent state of McEventCollection..."
374 const std::pair<unsigned int,unsigned int>
stats = nbrParticlesAndVertices( transObj );
384 const HepMC::GenEvent* genEvt = *itr;
387 auto ri = genEvt->run_info();
388 HepMC3::GenRunInfoData ri_data;
390 ri->write_data(ri_data);
391 if (!ri_data.weight_names.empty()) {
392 m_hepMCWeightSvc->setWeightNames( names_to_name_index_map(ri_data.weight_names), ctx ).ignore();
400 auto A_mpi=genEvt->attribute<HepMC3::IntAttribute>(
"mpi");
401 auto A_signal_process_id=genEvt->attribute<HepMC3::IntAttribute>(
"signal_process_id");
402 auto A_event_scale=genEvt->attribute<HepMC3::DoubleAttribute>(
"event_scale");
403 auto A_alphaQCD=genEvt->attribute<HepMC3::DoubleAttribute>(
"alphaQCD");
404 auto A_alphaQED=genEvt->attribute<HepMC3::DoubleAttribute>(
"alphaQED");
405 auto A_filterWeight=genEvt->attribute<HepMC3::DoubleAttribute>(
"filterWeight");
406 auto A_filterHT=genEvt->attribute<HepMC3::DoubleAttribute>(
"filterHT");
407 auto A_filterMET=genEvt->attribute<HepMC3::DoubleAttribute>(
"filterMET");
412 auto event_spv = genEvt->attribute<HepMC3::IntAttribute>(
"signal_process_vertex");
418 auto A_random_states=genEvt->attribute<HepMC3::VectorLongIntAttribute>(
"random_states");
419 auto beams=genEvt->beams();
421 emplace_back(A_signal_process_id?(A_signal_process_id->value()):-1,
422 genEvt->event_number(),
423 A_mpi?(A_mpi->value()):-1,
424 A_event_scale?(A_event_scale->value()):0.0,
425 A_alphaQCD?(A_alphaQCD->value()):0.0,
426 A_alphaQED?(A_alphaQED->value()):0.0,
427 A_filterWeight?(A_filterWeight->value()):1.0,
428 A_filterHT?(A_filterHT->value()):-13.,
429 A_filterMET?(A_filterMET->value()):-13.0,
434 A_random_states?(A_random_states->value()):std::vector<long>(),
435 std::vector<double>(),
436 std::vector<float>(),
437 std::vector<double>(),
438 genEvt->momentum_unit(),
439 genEvt->length_unit(),
441 nPersVtx + genEvt->vertices().size(),
443 nPersParts + genEvt->particles().size() );
446 std::map< std::string, std::map<int, std::shared_ptr<HepMC3::Attribute> > > e_atts = genEvt->attributes();
450 for (
auto& attmap: e_atts) {
451 if (attributes_to_ignore.count(attmap.first))
continue;
452 if (attmap.first ==
"ShadowParticle")
continue;
453 if (attmap.first ==
"ShadowParticleId")
continue;
454 for (
auto& att: attmap.second) {
458 att.second->to_string(
st);
478 if (genEvt->cross_section()) {
479 auto cs=genEvt->cross_section();
499 if (genEvt->heavy_ion()) {
500 auto hi=genEvt->heavy_ion();
502 std::vector<float>& heavyIon = persEvt.
m_heavyIon;
504 heavyIon[12] =
static_cast<float>(hi->Ncoll_hard);
505 heavyIon[11] =
static_cast<float>(hi->Npart_proj);
506 heavyIon[10] =
static_cast<float>(hi->Npart_targ);
507 heavyIon[9] =
static_cast<float>(hi->Ncoll);
508 heavyIon[8] =
static_cast<float>(hi->spectator_neutrons);
509 heavyIon[7] =
static_cast<float>(hi->spectator_protons);
510 heavyIon[6] =
static_cast<float>(hi->N_Nwounded_collisions);
511 heavyIon[5] =
static_cast<float>(hi->Nwounded_N_collisions);
512 heavyIon[4] =
static_cast<float>(hi->Nwounded_Nwounded_collisions);
513 heavyIon[3] = hi->impact_parameter;
514 heavyIon[2] = hi->event_plane_angle;
515 heavyIon[1] = hi->eccentricity;
516 heavyIon[0] = hi->sigma_inel_NN;
520 if (genEvt->pdf_info()) {
521 auto pi=genEvt->pdf_info();
523 std::vector<double>& pdfinfo = persEvt.
m_pdfinfo;
525 pdfinfo[8] =
static_cast<double>(
pi->parton_id[0]);
526 pdfinfo[7] =
static_cast<double>(
pi->parton_id[1]);
527 pdfinfo[6] =
static_cast<double>(
pi->pdf_id[0]);
528 pdfinfo[5] =
static_cast<double>(
pi->pdf_id[1]);
529 pdfinfo[4] =
pi->x[0];
530 pdfinfo[3] =
pi->x[1];
531 pdfinfo[2] =
pi->scale;
532 pdfinfo[1] =
pi->xf[0];
533 pdfinfo[0] =
pi->xf[1];
537 for (
const auto&
v: genEvt->vertices()) {
541 const int signalProcessVtx = genEvt->m_signal_process_vertex
542 ? genEvt->m_signal_process_vertex->barcode()
544 const int beamParticle1Barcode = genEvt->m_beam_particle_1
545 ? genEvt->m_beam_particle_1->barcode()
547 const int beamParticle2Barcode = genEvt->m_beam_particle_2
548 ? genEvt->m_beam_particle_2->barcode()
552 m_hepMCWeightSvc->setWeightNames( genEvt->m_weights.m_names, ctx ).ignore();
556 push_back(
GenEvent_p6( genEvt->m_signal_process_id,
557 genEvt->m_event_number,
559 genEvt->m_event_scale,
564 beamParticle1Barcode,
565 beamParticle2Barcode,
566 genEvt->m_weights.m_weights,
567 genEvt->m_random_states,
568 std::vector<double>(),
569 std::vector<float>(),
570 std::vector<double>(),
571 genEvt->m_momentum_unit,
572 genEvt->m_position_unit,
574 nPersVtx + genEvt->vertices_size(),
576 nPersParts + genEvt->particles_size() ) );
578 if (genEvt->m_cross_section) {
582 crossSection[2] = genEvt->m_cross_section->m_cross_section;
583 crossSection[1] = genEvt->m_cross_section->m_cross_section_error;
584 crossSection[0] =
static_cast<double>(genEvt->m_cross_section->m_is_set);
588 if (genEvt->m_heavy_ion) {
590 std::vector<float>& heavyIon = persEvt.
m_heavyIon;
592 heavyIon[12] =
static_cast<float>(genEvt->m_heavy_ion->m_Ncoll_hard);
593 heavyIon[11] =
static_cast<float>(genEvt->m_heavy_ion->m_Npart_proj);
594 heavyIon[10] =
static_cast<float>(genEvt->m_heavy_ion->m_Npart_targ);
595 heavyIon[9] =
static_cast<float>(genEvt->m_heavy_ion->m_Ncoll);
596 heavyIon[8] =
static_cast<float>(genEvt->m_heavy_ion->m_spectator_neutrons);
597 heavyIon[7] =
static_cast<float>(genEvt->m_heavy_ion->m_spectator_protons);
598 heavyIon[6] =
static_cast<float>(genEvt->m_heavy_ion->m_N_Nwounded_collisions);
599 heavyIon[5] =
static_cast<float>(genEvt->m_heavy_ion->m_Nwounded_N_collisions);
600 heavyIon[4] =
static_cast<float>(genEvt->m_heavy_ion->m_Nwounded_Nwounded_collisions);
601 heavyIon[3] = genEvt->m_heavy_ion->m_impact_parameter;
602 heavyIon[2] = genEvt->m_heavy_ion->m_event_plane_angle;
603 heavyIon[1] = genEvt->m_heavy_ion->m_eccentricity;
604 heavyIon[0] = genEvt->m_heavy_ion->m_sigma_inel_NN;
608 if (genEvt->m_pdf_info) {
610 std::vector<double>& pdfinfo = persEvt.
m_pdfinfo;
612 pdfinfo[8] =
static_cast<double>(genEvt->m_pdf_info->m_id1);
613 pdfinfo[7] =
static_cast<double>(genEvt->m_pdf_info->m_id2);
614 pdfinfo[6] =
static_cast<double>(genEvt->m_pdf_info->m_pdf_id1);
615 pdfinfo[5] =
static_cast<double>(genEvt->m_pdf_info->m_pdf_id2);
616 pdfinfo[4] = genEvt->m_pdf_info->m_x1;
617 pdfinfo[3] = genEvt->m_pdf_info->m_x2;
618 pdfinfo[2] = genEvt->m_pdf_info->m_scalePDF;
619 pdfinfo[1] = genEvt->m_pdf_info->m_pdf1;
620 pdfinfo[0] = genEvt->m_pdf_info->m_pdf2;
624 const HepMC::GenEvent::vertex_const_iterator endVtx=genEvt->vertices_end();
625 for ( HepMC::GenEvent::vertex_const_iterator
i = genEvt->vertices_begin();
653 vtx->set_position(HepMC::FourVector( persVtx.
m_x , persVtx.
m_y , persVtx.
m_z ,persVtx.
m_t ));
658 vtx->add_attribute(
"weights",std::make_shared<HepMC3::VectorDoubleAttribute>(
weights));
662 for (
unsigned int i = 0;
i != nPartsIn; ++
i ) {
668 for (
unsigned int i = 0;
i != nPartsOut; ++
i ) {
672 vtx->m_position.setX( persVtx.
m_x );
673 vtx->m_position.setY( persVtx.
m_y );
674 vtx->m_position.setZ( persVtx.
m_z );
675 vtx->m_position.setT( persVtx.
m_t );
676 vtx->m_particles_in.clear();
677 vtx->m_particles_out.clear();
679 vtx->m_weights.m_weights.reserve( persVtx.
m_weights.size() );
680 vtx->m_weights.m_weights.assign ( persVtx.
m_weights.begin(),
687 for (
unsigned int i = 0;
i != nPartsIn; ++
i ) {
695 for (
unsigned int i = 0;
i != nPartsOut; ++
i ) {
719 p->add_attribute(
"phi",std::make_shared<HepMC3::DoubleAttribute>(persPart.
m_phiPolarization));
720 p->add_attribute(
"theta",std::make_shared<HepMC3::DoubleAttribute>(persPart.
m_thetaPolarization));
730 double temp_e = std::sqrt( (
long double)(persPart.
m_px)*persPart.
m_px +
731 (
long double)(persPart.
m_py)*persPart.
m_py +
732 (
long double)(persPart.
m_pz)*persPart.
m_pz +
733 (
long double)(persPart.
m_m) *persPart.
m_m );
734 p->set_momentum( HepMC::FourVector(persPart.
m_px,persPart.
m_py,persPart.
m_pz,temp_e));
736 const int signM2 = ( persPart.
m_m >= 0. ? 1 : -1 );
737 const double persPart_ene =
738 std::sqrt( std::abs((
long double)(persPart.
m_px)*persPart.
m_px +
739 (
long double)(persPart.
m_py)*persPart.
m_py +
740 (
long double)(persPart.
m_pz)*persPart.
m_pz +
741 signM2* (
long double)(persPart.
m_m)* persPart.
m_m));
742 const int signEne = ( persPart.
m_recoMethod == 1 ? 1 : -1 );
743 p->set_momentum(HepMC::FourVector( persPart.
m_px,
746 signEne * persPart_ene ));
750 std::vector<int> flows;
751 const unsigned int nFlow = persPart.
m_flow.size();
752 for (
unsigned int iFlow= 0; iFlow != nFlow; ++iFlow ) {
753 flows.push_back(persPart.
m_flow[iFlow].second );
756 p->add_attribute(
"flows", std::make_shared<HepMC3::VectorIntAttribute>(flows));
762 p->m_production_vertex = 0;
774 p->m_momentum.setPx( persPart.
m_px);
775 p->m_momentum.setPy( persPart.
m_py);
776 p->m_momentum.setPz( persPart.
m_pz);
777 double temp_e = std::sqrt( (
long double)(persPart.
m_px)*persPart.
m_px +
778 (
long double)(persPart.
m_py)*persPart.
m_py +
779 (
long double)(persPart.
m_pz)*persPart.
m_pz +
780 (
long double)(persPart.
m_m) *persPart.
m_m );
781 p->m_momentum.setE( temp_e);
783 const int signM2 = ( persPart.
m_m >= 0. ? 1 : -1 );
784 const double persPart_ene =
785 std::sqrt( std::abs((
long double)(persPart.
m_px)*persPart.
m_px +
786 (
long double)(persPart.
m_py)*persPart.
m_py +
787 (
long double)(persPart.
m_pz)*persPart.
m_pz +
788 signM2* (
long double)(persPart.
m_m)* persPart.
m_m));
789 const int signEne = ( persPart.
m_recoMethod == 1 ? 1 : -1 );
790 p->m_momentum.set( persPart.
m_px,
793 signEne * persPart_ene );
797 const unsigned int nFlow = persPart.
m_flow.size();
799 for (
unsigned int iFlow= 0; iFlow != nFlow; ++iFlow ) {
800 p->m_flow.set_icode( persPart.
m_flow[iFlow].first,
801 persPart.
m_flow[iFlow].second );
816 const HepMC::FourVector& position = vtx->position();
817 auto A_weights=vtx->attribute<HepMC3::VectorDoubleAttribute>(
"weights");
818 auto A_barcode=vtx->attribute<HepMC3::IntAttribute>(
"barcode");
821 auto weights_d = A_weights->value();
822 for (
auto&
w: weights_d)
weights.push_back(
w);
831 A_barcode?(A_barcode->value()):vtx->id() );
836 for (
const auto&
p: vtx->particles_in()) {
837 if ( !
p->production_vertex() ||
p->production_vertex()->id() == 0 ) {
843 for (
const auto&
p: vtx->particles_out()) {
852 const HepMC::FourVector& position = vtx.m_position;
859 vtx.m_weights.m_weights.begin(),
860 vtx.m_weights.m_weights.end(),
865 const std::vector<HepMC::GenParticlePtr>::const_iterator endInVtx = vtx.m_particles_in.end();
867 for ( std::vector<HepMC::GenParticlePtr>::const_iterator
p = vtx.m_particles_in.begin();
870 if ( 0 == (*p)->production_vertex() ) {
875 const std::vector<HepMC::GenParticlePtr>::const_iterator endOutVtx = vtx.m_particles_out.end();
877 for ( std::vector<HepMC::GenParticlePtr>::const_iterator
p = vtx.m_particles_out.begin();
891 const HepMC::FourVector
mom =
p->momentum();
892 const double ene =
mom.e();
893 const double m2 =
mom.m2();
896 const bool useP2M2 = !(
m2 > 0) &&
898 !(std::abs(
m2) < 2.0*DBL_EPSILON*ene*ene);
899 auto A_flows=
p->attribute<HepMC3::VectorIntAttribute>(
"flows");
900 auto A_phi=
p->attribute<HepMC3::DoubleAttribute>(
"phi");
901 auto A_theta=
p->attribute<HepMC3::DoubleAttribute>(
"theta");
903 const short recoMethod = ( !useP2M2 ? 0: ( ene >= 0.? 1: 2 ) );
905 emplace_back(
mom.px(),
911 A_flows?(A_flows->value().size()):0,
912 A_theta?(A_theta->value()):0.0,
913 A_phi?(A_phi->value()):0.0,
920 std::vector< std::pair<int,int> > flow_hepmc2;
921 if(A_flows) flow_hepmc2=vector_to_vector_int_int(A_flows->value());
922 persEvt.
m_genParticles.back().m_flow.assign( flow_hepmc2.begin(),flow_hepmc2.end() );
933 const HepMC::FourVector&
mom =
p.m_momentum;
934 const double ene =
mom.e();
935 const double m2 =
mom.m2();
938 const bool useP2M2 = !(
m2 > 0) &&
940 !(std::abs(
m2) < 2.0*DBL_EPSILON*ene*ene);
942 const short recoMethod = ( !useP2M2
956 p.m_polarization.theta(),
957 p.m_polarization.phi(),
958 p.m_production_vertex
959 ?
p.m_production_vertex->barcode()
962 ?
p.m_end_vertex->barcode()