22 #include "GaudiKernel/ThreadLocalContext.h"
23 static const std::set<std::string> attributes_to_ignore {
25 "flows",
"flow1",
"flow2",
"flow3",
30 "filterWeight",
"filterHT",
"filterMET",
31 "event_scale",
"alphaQCD",
"alphaQED",
"random_states",
"weights",
32 "GenCrossSection",
"GenPdfInfo",
"GenHeavyIon"};
40 m_isPileup(false),m_hepMCWeightSvc(
"HepMCWeightSvc",
"McEventCollectionCnv_p7")
45 m_isPileup(false),m_hepMCWeightSvc(
"HepMCWeightSvc",
"McEventCollectionCnv_p7")
52 Base_t::operator=( rhs );
70 const EventContext& ctx = Gaudi::Hive::currentContext();
72 msg <<
MSG::DEBUG <<
"Loading McEventCollection from persistent state..."
83 datapools.
part.prepareToAdd(nParts);
84 const unsigned int nEvts = persObj->
m_genEvents.size();
85 datapools.
evt.prepareToAdd(nEvts);
88 for ( std::vector<GenEvent_p7>::const_iterator
94 HepMC::GenEvent * genEvt(
nullptr);
96 genEvt =
new HepMC::GenEvent();
101 genEvt->add_attribute (
"barcodes", std::make_shared<HepMC::GenEventBarcodes>());
108 genEvt->add_attribute(
"signal_process_id", std::make_shared<HepMC3::IntAttribute>(persEvt.
m_signalProcessId));
110 genEvt->add_attribute(
"mpi", std::make_shared<HepMC3::IntAttribute>(persEvt.
m_mpi));
111 genEvt->add_attribute(
"event_scale", std::make_shared<HepMC3::DoubleAttribute>(persEvt.
m_eventScale));
112 genEvt->add_attribute(
"alphaQCD", std::make_shared<HepMC3::DoubleAttribute>(persEvt.
m_alphaQCD));
113 genEvt->add_attribute(
"alphaQED", std::make_shared<HepMC3::DoubleAttribute>(persEvt.
m_alphaQED));
114 genEvt->add_attribute(
"filterWeight", std::make_shared<HepMC3::DoubleAttribute>(persEvt.
m_filterWeight));
115 genEvt->add_attribute(
"filterHT", std::make_shared<HepMC3::DoubleAttribute>(persEvt.m_filterHT));
116 genEvt->add_attribute(
"filterMET", std::make_shared<HepMC3::DoubleAttribute>(persEvt.m_filterMET));
118 genEvt->add_attribute(
"random_states", std::make_shared<HepMC3::VectorLongIntAttribute>(persEvt.
m_randomStates));
120 genEvt->set_units(
static_cast<HepMC3::Units::MomentumUnit
>(persEvt.
m_momentumUnit),
121 static_cast<HepMC3::Units::LengthUnit
>(persEvt.
m_lengthUnit));
124 if(!genEvt->run_info()) {
125 HepMC3::GenRunInfoData ri_read;
132 auto ri = std::make_shared<HepMC3::GenRunInfo>();
133 ri->read_data(ri_read);
134 genEvt->set_run_info(ri);
139 auto cs = std::make_shared<HepMC3::GenCrossSection>();
141 if(
static_cast<bool>(xsection[0]) )
142 cs->set_cross_section(xsection[2],xsection[1]);
144 cs->set_cross_section(-1.0, -1.0);
145 genEvt->set_cross_section(cs);
150 auto hi = std::make_shared<HepMC3::GenHeavyIon>();
151 const std::vector<float>& hIon = persEvt.
m_heavyIon;
154 static_cast<int>(hIon[12]),
155 static_cast<int>(hIon[11]),
156 static_cast<int>(hIon[10]),
157 static_cast<int>(hIon[9]),
158 static_cast<int>(hIon[8]),
159 static_cast<int>(hIon[7]),
160 static_cast<int>(hIon[6]),
161 static_cast<int>(hIon[5]),
162 static_cast<int>(hIon[4]),
167 genEvt->set_heavy_ion(hi);
177 pi->set(
static_cast<int>(
pdf[8]),
178 static_cast<int>(
pdf[7]),
184 static_cast<int>(
pdf[6]),
185 static_cast<int>(
pdf[5]));
186 genEvt->set_pdf_info(
pi);
196 std::map<int, HepMC::GenVertexPtr> brc_to_vertex;
200 for (
unsigned int iVtx = persEvt.
m_verticesBegin; iVtx != endVtx; ++iVtx ) {
207 if ( sigProcVtx != 0 && brc_to_vertex.count(sigProcVtx) ) {
212 for (
auto &
p : partToEndVtx) {
213 if ( brc_to_vertex.count(
p.second) ) {
214 auto decayVtx = brc_to_vertex[
p.second];
215 decayVtx->add_particle_in(
p.first );
217 msg << MSG::ERROR <<
"GenParticle points to null end vertex !!" <<
endmsg;
223 if ( beamPart1 != 0 && beamPart2 != 0 ) {
231 genEvt->m_mpi = persEvt.
m_mpi;
235 genEvt->m_signal_process_vertex = 0;
236 genEvt->m_beam_particle_1 = 0;
237 genEvt->m_beam_particle_2 = 0;
240 genEvt->m_vertex_barcodes.clear();
241 genEvt->m_particle_barcodes.clear();
242 genEvt->m_momentum_unit =
static_cast<HepMC::Units::MomentumUnit
>(persEvt.
m_momentumUnit);
243 genEvt->m_position_unit =
static_cast<HepMC::Units::LengthUnit
>(persEvt.
m_lengthUnit);
249 if( genEvt->m_cross_section )
250 delete genEvt->m_cross_section;
251 genEvt->m_cross_section = 0;
254 genEvt->m_cross_section =
new HepMC::GenCrossSection();
256 if(
static_cast<bool>(xsection[0]) )
257 genEvt->m_cross_section->set_cross_section(xsection[2],xsection[1]);
261 if(genEvt->m_heavy_ion )
262 delete genEvt->m_heavy_ion;
263 genEvt->m_heavy_ion = 0;
265 const std::vector<float>& hIon = persEvt.
m_heavyIon;
266 genEvt->m_heavy_ion =
new HepMC::HeavyIon
268 static_cast<int>(hIon[12]),
269 static_cast<int>(hIon[11]),
270 static_cast<int>(hIon[10]),
271 static_cast<int>(hIon[9]),
272 static_cast<int>(hIon[8]),
273 static_cast<int>(hIon[7]),
274 static_cast<int>(hIon[6]),
275 static_cast<int>(hIon[5]),
276 static_cast<int>(hIon[4]),
286 if(genEvt->m_pdf_info)
287 delete genEvt->m_pdf_info;
288 genEvt->m_pdf_info = 0;
291 genEvt->m_pdf_info =
new HepMC::PdfInfo
293 static_cast<int>(
pdf[8]),
294 static_cast<int>(
pdf[7]),
300 static_cast<int>(
pdf[6]),
301 static_cast<int>(
pdf[5])
316 for (
unsigned int iVtx= persEvt.
m_verticesBegin; iVtx != endVtx; ++iVtx ) {
325 if ( sigProcVtx != 0 ) {
326 genEvt->set_signal_process_vertex( genEvt->barcode_to_vertex( sigProcVtx ) );
331 p = partToEndVtx.begin(),
332 endItr = partToEndVtx.end();
337 decayVtx->add_particle_in(
p->first );
340 <<
"GenParticle points to null end vertex !!"
348 if ( beamPart1 != 0 && beamPart2 !=0 ) {
349 genEvt->set_beam_particles(genEvt->barcode_to_particle(beamPart1),
350 genEvt->barcode_to_particle(beamPart2));
358 msg <<
MSG::DEBUG <<
"Loaded McEventCollection from persistent state [OK]"
366 const EventContext& ctx = Gaudi::Hive::currentContext();
368 msg <<
MSG::DEBUG <<
"Creating persistent state of McEventCollection..."
372 const std::pair<unsigned int,unsigned int>
stats = nbrParticlesAndVertices( transObj );
382 const HepMC::GenEvent* genEvt = *itr;
384 auto ri = genEvt->run_info();
386 HepMC3::GenRunInfoData ri_data;
388 ri->write_data(ri_data);
389 if (!ri_data.weight_names.empty()) {
390 m_hepMCWeightSvc->setWeightNames( names_to_name_index_map(ri_data.weight_names), ctx ).ignore();
398 auto A_mpi=genEvt->attribute<HepMC3::IntAttribute>(
"mpi");
399 auto A_signal_process_id=genEvt->attribute<HepMC3::IntAttribute>(
"signal_process_id");
400 auto A_event_scale=genEvt->attribute<HepMC3::DoubleAttribute>(
"event_scale");
401 auto A_alphaQCD=genEvt->attribute<HepMC3::DoubleAttribute>(
"alphaQCD");
402 auto A_alphaQED=genEvt->attribute<HepMC3::DoubleAttribute>(
"alphaQED");
403 auto A_filterWeight=genEvt->attribute<HepMC3::DoubleAttribute>(
"filterWeight");
404 auto A_filterHT=genEvt->attribute<HepMC3::DoubleAttribute>(
"filterHT");
405 auto A_filterMET=genEvt->attribute<HepMC3::DoubleAttribute>(
"filterMET");
407 auto A_random_states=genEvt->attribute<HepMC3::VectorLongIntAttribute>(
"random_states");
408 auto beams=genEvt->beams();
410 emplace_back(A_signal_process_id?(A_signal_process_id->value()):-1,
411 genEvt->event_number(),
412 A_mpi?(A_mpi->value()):-1,
413 A_event_scale?(A_event_scale->value()):0.0,
414 A_alphaQCD?(A_alphaQCD->value()):0.0,
415 A_alphaQED?(A_alphaQED->value()):0.0,
416 A_filterWeight?(A_filterWeight->value()):1.0,
417 A_filterHT?(A_filterHT->value()):-13.,
418 A_filterMET?(A_filterMET->value()):-13.0,
423 A_random_states?(A_random_states->value()):std::vector<long>(),
424 std::vector<double>(),
425 std::vector<float>(),
426 std::vector<double>(),
427 genEvt->momentum_unit(),
428 genEvt->length_unit(),
430 nPersVtx + genEvt->vertices().size(),
432 nPersParts + genEvt->particles().size() );
435 std::map< std::string, std::map<int, std::shared_ptr<HepMC3::Attribute> > > e_atts = genEvt->attributes();
439 for (
auto& attmap: e_atts) {
440 if (attributes_to_ignore.count(attmap.first))
continue;
441 if (attmap.first ==
"ShadowParticle")
continue;
442 if (attmap.first ==
"ShadowParticleId")
continue;
443 for (
auto& att: attmap.second) {
447 att.second->to_string(st);
470 if (genEvt->cross_section()) {
471 auto cs=genEvt->cross_section();
491 if (genEvt->heavy_ion()) {
492 auto hi=genEvt->heavy_ion();
494 std::vector<float>& heavyIon = persEvt.
m_heavyIon;
496 heavyIon[12] =
static_cast<float>(hi->Ncoll_hard);
497 heavyIon[11] =
static_cast<float>(hi->Npart_proj);
498 heavyIon[10] =
static_cast<float>(hi->Npart_targ);
499 heavyIon[9] =
static_cast<float>(hi->Ncoll);
500 heavyIon[8] =
static_cast<float>(hi->spectator_neutrons);
501 heavyIon[7] =
static_cast<float>(hi->spectator_protons);
502 heavyIon[6] =
static_cast<float>(hi->N_Nwounded_collisions);
503 heavyIon[5] =
static_cast<float>(hi->Nwounded_N_collisions);
504 heavyIon[4] =
static_cast<float>(hi->Nwounded_Nwounded_collisions);
505 heavyIon[3] = hi->impact_parameter;
506 heavyIon[2] = hi->event_plane_angle;
507 heavyIon[1] = hi->eccentricity;
508 heavyIon[0] = hi->sigma_inel_NN;
512 if (genEvt->pdf_info()) {
513 auto pi=genEvt->pdf_info();
515 std::vector<double>& pdfinfo = persEvt.
m_pdfinfo;
517 pdfinfo[8] =
static_cast<double>(
pi->parton_id[0]);
518 pdfinfo[7] =
static_cast<double>(
pi->parton_id[1]);
519 pdfinfo[6] =
static_cast<double>(
pi->pdf_id[0]);
520 pdfinfo[5] =
static_cast<double>(
pi->pdf_id[1]);
521 pdfinfo[4] =
pi->x[0];
522 pdfinfo[3] =
pi->x[1];
523 pdfinfo[2] =
pi->scale;
524 pdfinfo[1] =
pi->xf[0];
525 pdfinfo[0] =
pi->xf[1];
529 for (
const auto&
v: genEvt->vertices()) {
533 const int signalProcessVtx = genEvt->m_signal_process_vertex
534 ? genEvt->m_signal_process_vertex->barcode()
536 const int beamParticle1Barcode = genEvt->m_beam_particle_1
537 ? genEvt->m_beam_particle_1->barcode()
539 const int beamParticle2Barcode = genEvt->m_beam_particle_2
540 ? genEvt->m_beam_particle_2->barcode()
544 m_hepMCWeightSvc->setWeightNames( genEvt->m_weights.m_names, ctx ).ignore();
548 push_back(
GenEvent_p7( genEvt->m_signal_process_id,
549 genEvt->m_event_number,
551 genEvt->m_event_scale,
556 beamParticle1Barcode,
557 beamParticle2Barcode,
558 genEvt->m_weights.m_weights,
559 genEvt->m_random_states,
560 std::vector<double>(),
561 std::vector<float>(),
562 std::vector<double>(),
563 genEvt->m_momentum_unit,
564 genEvt->m_position_unit,
566 nPersVtx + genEvt->vertices_size(),
568 nPersParts + genEvt->particles_size() ) );
570 if (genEvt->m_cross_section) {
574 crossSection[2] = genEvt->m_cross_section->m_cross_section;
575 crossSection[1] = genEvt->m_cross_section->m_cross_section_error;
576 crossSection[0] =
static_cast<double>(genEvt->m_cross_section->m_is_set);
580 if (genEvt->m_heavy_ion) {
582 std::vector<float>& heavyIon = persEvt.
m_heavyIon;
584 heavyIon[12] =
static_cast<float>(genEvt->m_heavy_ion->m_Ncoll_hard);
585 heavyIon[11] =
static_cast<float>(genEvt->m_heavy_ion->m_Npart_proj);
586 heavyIon[10] =
static_cast<float>(genEvt->m_heavy_ion->m_Npart_targ);
587 heavyIon[9] =
static_cast<float>(genEvt->m_heavy_ion->m_Ncoll);
588 heavyIon[8] =
static_cast<float>(genEvt->m_heavy_ion->m_spectator_neutrons);
589 heavyIon[7] =
static_cast<float>(genEvt->m_heavy_ion->m_spectator_protons);
590 heavyIon[6] =
static_cast<float>(genEvt->m_heavy_ion->m_N_Nwounded_collisions);
591 heavyIon[5] =
static_cast<float>(genEvt->m_heavy_ion->m_Nwounded_N_collisions);
592 heavyIon[4] =
static_cast<float>(genEvt->m_heavy_ion->m_Nwounded_Nwounded_collisions);
593 heavyIon[3] = genEvt->m_heavy_ion->m_impact_parameter;
594 heavyIon[2] = genEvt->m_heavy_ion->m_event_plane_angle;
595 heavyIon[1] = genEvt->m_heavy_ion->m_eccentricity;
596 heavyIon[0] = genEvt->m_heavy_ion->m_sigma_inel_NN;
600 if (genEvt->m_pdf_info) {
602 std::vector<double>& pdfinfo = persEvt.
m_pdfinfo;
604 pdfinfo[8] =
static_cast<double>(genEvt->m_pdf_info->m_id1);
605 pdfinfo[7] =
static_cast<double>(genEvt->m_pdf_info->m_id2);
606 pdfinfo[6] =
static_cast<double>(genEvt->m_pdf_info->m_pdf_id1);
607 pdfinfo[5] =
static_cast<double>(genEvt->m_pdf_info->m_pdf_id2);
608 pdfinfo[4] = genEvt->m_pdf_info->m_x1;
609 pdfinfo[3] = genEvt->m_pdf_info->m_x2;
610 pdfinfo[2] = genEvt->m_pdf_info->m_scalePDF;
611 pdfinfo[1] = genEvt->m_pdf_info->m_pdf1;
612 pdfinfo[0] = genEvt->m_pdf_info->m_pdf2;
616 const HepMC::GenEvent::vertex_const_iterator endVtx=genEvt->vertices_end();
617 for ( HepMC::GenEvent::vertex_const_iterator
i = genEvt->vertices_begin();
645 vtx->set_position(HepMC::FourVector( persVtx.
m_x , persVtx.
m_y , persVtx.
m_z ,persVtx.
m_t ));
647 vtx->set_status(persVtx.
m_id);
650 vtx->add_attribute(
"weights",std::make_shared<HepMC3::VectorDoubleAttribute>(
weights));
654 for (
unsigned int i = 0;
i != nPartsIn; ++
i ) {
660 for (
unsigned int i = 0;
i != nPartsOut; ++
i ) {
664 vtx->m_position.setX( persVtx.
m_x );
665 vtx->m_position.setY( persVtx.
m_y );
666 vtx->m_position.setZ( persVtx.
m_z );
667 vtx->m_position.setT( persVtx.
m_t );
668 vtx->m_particles_in.clear();
669 vtx->m_particles_out.clear();
670 vtx->m_id = persVtx.
m_id;
671 vtx->m_weights.m_weights.reserve( persVtx.
m_weights.size() );
672 vtx->m_weights.m_weights.assign ( persVtx.
m_weights.begin(),
679 for (
unsigned int i = 0;
i != nPartsIn; ++
i ) {
687 for (
unsigned int i = 0;
i != nPartsOut; ++
i ) {
711 p->add_attribute(
"phi",std::make_shared<HepMC3::DoubleAttribute>(persPart.
m_phiPolarization));
712 p->add_attribute(
"theta",std::make_shared<HepMC3::DoubleAttribute>(persPart.
m_thetaPolarization));
722 double temp_e = std::sqrt( (
long double)(persPart.
m_px)*persPart.
m_px +
723 (
long double)(persPart.
m_py)*persPart.
m_py +
724 (
long double)(persPart.
m_pz)*persPart.
m_pz +
725 (
long double)(persPart.
m_m) *persPart.
m_m );
726 p->set_momentum( HepMC::FourVector(persPart.
m_px,persPart.
m_py,persPart.
m_pz,temp_e));
728 const int signM2 = ( persPart.
m_m >= 0. ? 1 : -1 );
729 const double persPart_ene =
730 std::sqrt( std::abs((
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 signM2* (
long double)(persPart.
m_m)* persPart.
m_m));
734 const int signEne = ( persPart.
m_recoMethod == 1 ? 1 : -1 );
735 p->set_momentum(HepMC::FourVector( persPart.
m_px,
738 signEne * persPart_ene ));
742 std::vector<int> flows;
743 const unsigned int nFlow = persPart.
m_flow.size();
744 for (
unsigned int iFlow= 0; iFlow != nFlow; ++iFlow ) {
745 flows.push_back(persPart.
m_flow[iFlow].second );
748 p->add_attribute(
"flows", std::make_shared<HepMC3::VectorIntAttribute>(flows));
754 p->m_production_vertex = 0;
766 p->m_momentum.setPx( persPart.
m_px);
767 p->m_momentum.setPy( persPart.
m_py);
768 p->m_momentum.setPz( persPart.
m_pz);
769 double temp_e = std::sqrt( (
long double)(persPart.
m_px)*persPart.
m_px +
770 (
long double)(persPart.
m_py)*persPart.
m_py +
771 (
long double)(persPart.
m_pz)*persPart.
m_pz +
772 (
long double)(persPart.
m_m) *persPart.
m_m );
773 p->m_momentum.setE( temp_e);
775 const int signM2 = ( persPart.
m_m >= 0. ? 1 : -1 );
776 const double persPart_ene =
777 std::sqrt( std::abs((
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 signM2* (
long double)(persPart.
m_m)* persPart.
m_m));
781 const int signEne = ( persPart.
m_recoMethod == 1 ? 1 : -1 );
782 p->m_momentum.set( persPart.
m_px,
785 signEne * persPart_ene );
789 const unsigned int nFlow = persPart.
m_flow.size();
791 for (
unsigned int iFlow= 0; iFlow != nFlow; ++iFlow ) {
792 p->m_flow.set_icode( persPart.
m_flow[iFlow].first,
793 persPart.
m_flow[iFlow].second );
808 const HepMC::FourVector& position = vtx->position();
809 auto A_weights=vtx->attribute<HepMC3::VectorDoubleAttribute>(
"weights");
810 auto A_barcode=vtx->attribute<HepMC3::IntAttribute>(
"barcode");
813 auto weights_d = A_weights->value();
814 for (
auto&
w: weights_d)
weights.push_back(
w);
823 A_barcode?(A_barcode->value()):vtx->id() );
828 for (
const auto&
p: vtx->particles_in()) {
829 if ( !
p->production_vertex() ||
p->production_vertex()->id() == 0 ) {
835 for (
const auto&
p: vtx->particles_out()) {
844 const HepMC::FourVector& position = vtx.m_position;
851 vtx.m_weights.m_weights.begin(),
852 vtx.m_weights.m_weights.end(),
857 const std::vector<HepMC::GenParticlePtr>::const_iterator endInVtx = vtx.m_particles_in.end();
859 for ( std::vector<HepMC::GenParticlePtr>::const_iterator
p = vtx.m_particles_in.begin();
862 if ( 0 == (*p)->production_vertex() ) {
867 const std::vector<HepMC::GenParticlePtr>::const_iterator endOutVtx = vtx.m_particles_out.end();
869 for ( std::vector<HepMC::GenParticlePtr>::const_iterator
p = vtx.m_particles_out.begin();
883 const HepMC::FourVector
mom =
p->momentum();
884 const double ene =
mom.e();
885 const double m2 =
mom.m2();
888 const bool useP2M2 = !(
m2 > 0) &&
890 !(std::abs(
m2) < 2.0*DBL_EPSILON*ene*ene);
891 auto A_flows=
p->attribute<HepMC3::VectorIntAttribute>(
"flows");
892 auto A_phi=
p->attribute<HepMC3::DoubleAttribute>(
"phi");
893 auto A_theta=
p->attribute<HepMC3::DoubleAttribute>(
"theta");
895 const short recoMethod = ( !useP2M2 ? 0: ( ene >= 0.? 1: 2 ) );
897 emplace_back(
mom.px(),
903 A_flows?(A_flows->value().size()):0,
904 A_theta?(A_theta->value()):0.0,
905 A_phi?(A_phi->value()):0.0,
912 std::vector< std::pair<int,int> > flow_hepmc2;
913 if(A_flows) flow_hepmc2=vector_to_vector_int_int(A_flows->value());
914 persEvt.
m_genParticles.back().m_flow.assign( flow_hepmc2.begin(),flow_hepmc2.end() );
925 const HepMC::FourVector&
mom =
p.m_momentum;
926 const double ene =
mom.e();
927 const double m2 =
mom.m2();
930 const bool useP2M2 = !(
m2 > 0) &&
932 !(std::abs(
m2) < 2.0*DBL_EPSILON*ene*ene);
934 const short recoMethod = ( !useP2M2
948 p.m_polarization.theta(),
949 p.m_polarization.phi(),
950 p.m_production_vertex
951 ?
p.m_production_vertex->barcode()
954 ?
p.m_end_vertex->barcode()