25 #include "GaudiKernel/ThreadLocalContext.h"
35 m_isPileup(false),m_hepMCWeightSvc(
"HepMCWeightSvc",
"McEventCollectionCnv_p4")
40 m_isPileup(false),m_hepMCWeightSvc(
"HepMCWeightSvc",
"McEventCollectionCnv_p4")
47 Base_t::operator=( rhs );
68 const EventContext& ctx = Gaudi::Hive::currentContext();
70 msg <<
MSG::DEBUG <<
"Loading McEventCollection from persistent state..."
82 datapools.
part.prepareToAdd(nParts);
83 const unsigned int nEvts = persObj->
m_genEvents.size();
84 datapools.
evt.prepareToAdd(nEvts);
87 for ( std::vector<GenEvent_p4>::const_iterator
94 HepMC::GenEvent * genEvt(
nullptr);
97 genEvt =
new HepMC::GenEvent();
104 genEvt->add_attribute (
"barcodes", std::make_shared<HepMC::GenEventBarcodes>());
105 genEvt->add_attribute(
"signal_process_id", std::make_shared<HepMC3::IntAttribute>(persEvt.
m_signalProcessId));
107 genEvt->add_attribute(
"event_scale", std::make_shared<HepMC3::DoubleAttribute>(persEvt.
m_eventScale));
108 genEvt->add_attribute(
"alphaQCD", std::make_shared<HepMC3::DoubleAttribute>(persEvt.
m_alphaQCD));
109 genEvt->add_attribute(
"alphaQED", std::make_shared<HepMC3::DoubleAttribute>(persEvt.
m_alphaQED));
111 genEvt->add_attribute(
"random_states", std::make_shared<HepMC3::VectorLongIntAttribute>(persEvt.
m_randomStates));
113 if(!genEvt->run_info()) genEvt->set_run_info(std::make_shared<HepMC3::GenRunInfo>());
114 if(genEvt->run_info()) genEvt->run_info()->set_weight_names(
m_hepMCWeightSvc->weightNameVec(ctx));
123 static_cast<int>(
pdf[6]),
124 static_cast<int>(
pdf[5]),
130 genEvt->set_pdf_info(
pi);
141 std::map<int, HepMC::GenVertexPtr> brc_to_vertex;
144 for (
unsigned int iVtx= persEvt.
m_verticesBegin; iVtx != endVtx; ++iVtx )
152 if ( sigProcVtx != 0 && brc_to_vertex.count(sigProcVtx) ) {
157 for (
auto &
p : partToEndVtx) {
158 if ( brc_to_vertex.count(
p.second) ) {
159 auto decayVtx = brc_to_vertex[
p.second];
160 decayVtx->add_particle_in(
p.first );
162 msg << MSG::ERROR <<
"GenParticle points to null end vertex !!" <<
endmsg;
171 genEvt->m_signal_process_vertex = 0;
174 genEvt->m_vertex_barcodes.clear();
175 genEvt->m_particle_barcodes.clear();
180 delete genEvt->m_pdf_info; genEvt->m_pdf_info = 0;
184 genEvt->m_pdf_info =
new HepMC::PdfInfo
185 (
static_cast<int>(
pdf[6]),
186 static_cast<int>(
pdf[5]),
206 for (
unsigned int iVtx= persEvt.
m_verticesBegin; iVtx != endVtx; ++iVtx )
216 if ( sigProcVtx != 0 )
224 p = partToEndVtx.begin(),
225 endItr = partToEndVtx.end();
232 decayVtx->add_particle_in(
p->first );
237 <<
"GenParticle points to null end vertex !!"
245 msg <<
MSG::DEBUG <<
"Loaded McEventCollection from persistent state [OK]"
253 const EventContext& ctx = Gaudi::Hive::currentContext();
255 msg <<
MSG::DEBUG <<
"Creating persistent state of McEventCollection..."
259 const std::pair<unsigned int,unsigned int>
stats = nbrParticlesAndVertices( transObj );
271 const HepMC::GenEvent* genEvt = *itr;
273 if (genEvt->run_info()) {
274 if (!genEvt->run_info()->weight_names().empty()) {
275 m_hepMCWeightSvc->setWeightNames( names_to_name_index_map(genEvt->weight_names()), ctx ).ignore();
282 auto A_signal_process_id=genEvt->attribute<HepMC3::IntAttribute>(
"signal_process_id");
283 auto A_event_scale=genEvt->attribute<HepMC3::DoubleAttribute>(
"event_scale");
284 auto A_alphaQCD=genEvt->attribute<HepMC3::DoubleAttribute>(
"alphaQCD");
285 auto A_alphaQED=genEvt->attribute<HepMC3::DoubleAttribute>(
"alphaQED");
287 auto A_random_states=genEvt->attribute<HepMC3::VectorLongIntAttribute>(
"random_states");
290 emplace_back( A_signal_process_id?(A_signal_process_id->value()):0,
291 genEvt->event_number(),
292 A_event_scale?(A_event_scale->value()):0.0,
293 A_alphaQCD?(A_alphaQCD->value()):0.0,
294 A_alphaQED?(A_alphaQED->value()):0.0,
297 std::vector<double>(),
298 A_random_states?(A_random_states->value()):std::vector<long>(),
300 nPersVtx + genEvt->vertices().size(),
302 nPersParts + genEvt->particles().size() );
305 if (genEvt->pdf_info())
307 auto pi=genEvt->pdf_info();
309 std::vector<double>& pdfinfo = persEvt.
m_pdfinfo;
311 pdfinfo[6] =
static_cast<double>(
pi->parton_id[0]);
312 pdfinfo[5] =
static_cast<double>(
pi->parton_id[1]);
313 pdfinfo[4] =
pi->x[0];
314 pdfinfo[3] =
pi->x[1];
315 pdfinfo[2] =
pi->scale;
316 pdfinfo[1] =
pi->xf[0];
317 pdfinfo[0] =
pi->xf[1];
320 for (
const auto&
v: genEvt->vertices())
328 const HepMC::GenEvent* genEvt = *itr;
329 const int signalProcessVtx = genEvt->m_signal_process_vertex
330 ? genEvt->m_signal_process_vertex->barcode()
333 m_hepMCWeightSvc->setWeightNames( genEvt->m_weights.m_names, ctx ).ignore();
335 push_back(
GenEvent_p4( genEvt->m_signal_process_id,
336 genEvt->m_event_number,
337 genEvt->m_event_scale,
341 genEvt->m_weights.m_weights,
342 std::vector<double>(),
343 genEvt->m_random_states,
345 nPersVtx + genEvt->vertices_size(),
347 nPersParts + genEvt->particles_size() ) );
349 if (genEvt->m_pdf_info)
352 std::vector<double>& pdfinfo = persEvt.
m_pdfinfo;
354 pdfinfo[6] =
static_cast<double>(genEvt->m_pdf_info->m_id1);
355 pdfinfo[5] =
static_cast<double>(genEvt->m_pdf_info->m_id2);
356 pdfinfo[4] = genEvt->m_pdf_info->m_x1;
357 pdfinfo[3] = genEvt->m_pdf_info->m_x2;
358 pdfinfo[2] = genEvt->m_pdf_info->m_scalePDF;
359 pdfinfo[1] = genEvt->m_pdf_info->m_pdf1;
360 pdfinfo[0] = genEvt->m_pdf_info->m_pdf2;
364 const HepMC::GenEvent::vertex_const_iterator endVtx=genEvt->vertices_end();
365 for ( HepMC::GenEvent::vertex_const_iterator
i = genEvt->vertices_begin();
375 msg <<
MSG::DEBUG <<
"Created persistent state of HepMC::GenEvent [OK]"
397 vtx->set_position(HepMC::FourVector( persVtx.
m_x , persVtx.
m_y , persVtx.
m_z ,persVtx.
m_t ));
401 vtx->add_attribute(
"weights",std::make_shared<HepMC3::VectorDoubleAttribute>(
weights));
407 for (
unsigned int i = 0;
i != nPartsIn; ++
i )
414 for (
unsigned int i = 0;
i != nPartsOut; ++
i )
419 vtx->m_position.setX( persVtx.
m_x );
420 vtx->m_position.setY( persVtx.
m_y );
421 vtx->m_position.setZ( persVtx.
m_z );
422 vtx->m_position.setT( persVtx.
m_t );
423 vtx->m_particles_in.clear();
424 vtx->m_particles_out.clear();
426 vtx->m_weights.m_weights.reserve( persVtx.
m_weights.size() );
427 vtx->m_weights.m_weights.assign ( persVtx.
m_weights.begin(),
434 for (
unsigned int i = 0;
i != nPartsIn; ++
i )
443 for (
unsigned int i = 0;
i != nPartsOut; ++
i )
472 p->add_attribute(
"phi",std::make_shared<HepMC3::DoubleAttribute>(persPart.
m_phiPolarization));
473 p->add_attribute(
"theta",std::make_shared<HepMC3::DoubleAttribute>(persPart.
m_thetaPolarization));
483 double temp_e = std::sqrt( (
long double)(persPart.
m_px)*persPart.
m_px +
484 (
long double)(persPart.
m_py)*persPart.
m_py +
485 (
long double)(persPart.
m_pz)*persPart.
m_pz +
486 (
long double)(persPart.
m_m) *persPart.
m_m );
487 p->set_momentum( HepMC::FourVector(persPart.
m_px,persPart.
m_py,persPart.
m_pz,temp_e));
491 const int signM2 = ( persPart.
m_m >= 0. ? 1 : -1 );
492 const double persPart_ene =
493 std::sqrt( std::abs((
long double)(persPart.
m_px)*persPart.
m_px +
494 (
long double)(persPart.
m_py)*persPart.
m_py +
495 (
long double)(persPart.
m_pz)*persPart.
m_pz +
496 signM2* (
long double)(persPart.
m_m)* persPart.
m_m));
497 const int signEne = ( persPart.
m_recoMethod == 1 ? 1 : -1 );
498 p->set_momentum( HepMC::FourVector( persPart.
m_px,
501 signEne * persPart_ene ));
505 std::vector<int> flows;
506 const unsigned int nFlow = persPart.
m_flow.size();
507 for (
unsigned int iFlow= 0; iFlow != nFlow; ++iFlow ) {
508 flows.push_back(persPart.
m_flow[iFlow].second );
511 p->add_attribute(
"flows", std::make_shared<HepMC3::VectorIntAttribute>(flows));
517 p->m_production_vertex = 0;
529 p->m_momentum.setPx( persPart.
m_px);
530 p->m_momentum.setPy( persPart.
m_py);
531 p->m_momentum.setPz( persPart.
m_pz);
532 double temp_e = std::sqrt( (
long double)(persPart.
m_px)*persPart.
m_px +
533 (
long double)(persPart.
m_py)*persPart.
m_py +
534 (
long double)(persPart.
m_pz)*persPart.
m_pz +
535 (
long double)(persPart.
m_m) *persPart.
m_m );
536 p->m_momentum.setE( temp_e);
540 const int signM2 = ( persPart.
m_m >= 0. ? 1 : -1 );
541 const double persPart_ene =
542 std::sqrt( std::abs((
long double)(persPart.
m_px)*persPart.
m_px +
543 (
long double)(persPart.
m_py)*persPart.
m_py +
544 (
long double)(persPart.
m_pz)*persPart.
m_pz +
545 signM2* (
long double)(persPart.
m_m)* persPart.
m_m));
546 const int signEne = ( persPart.
m_recoMethod == 1 ? 1 : -1 );
547 p->m_momentum.set( persPart.
m_px,
550 signEne * persPart_ene );
554 const unsigned int nFlow = persPart.
m_flow.size();
556 for (
unsigned int iFlow= 0; iFlow != nFlow; ++iFlow )
558 p->m_flow.set_icode( persPart.
m_flow[iFlow].first,
559 persPart.
m_flow[iFlow].second );
575 const HepMC::FourVector& position = vtx->position();
576 auto A_weights=vtx->attribute<HepMC3::VectorDoubleAttribute>(
"weights");
577 auto A_barcode=vtx->attribute<HepMC3::IntAttribute>(
"barcode");
580 auto weights_d = A_weights->value();
581 for (
auto&
w: weights_d)
weights.push_back(
w);
590 A_barcode?(A_barcode->value()):vtx->id()
595 for (
const auto&
p: vtx->particles_in())
597 if ( !
p->production_vertex() ||
p->production_vertex()->id() == 0 )
603 for (
const auto&
p: vtx->particles_out())
612 const HepMC::FourVector& position = vtx.m_position;
619 vtx.m_weights.m_weights.begin(),
620 vtx.m_weights.m_weights.end(),
625 const std::vector<HepMC::GenParticlePtr>::const_iterator endInVtx = vtx.m_particles_in.end();
627 for ( std::vector<HepMC::GenParticlePtr>::const_iterator
p = vtx.m_particles_in.begin();
631 if ( 0 == (*p)->production_vertex() )
637 const std::vector<HepMC::GenParticlePtr>::const_iterator endOutVtx = vtx.m_particles_out.end();
639 for ( std::vector<HepMC::GenParticlePtr>::const_iterator
p = vtx.m_particles_out.begin();
654 const HepMC::FourVector&
mom =
p->momentum();
655 const double ene =
mom.e();
656 const double m2 =
mom.m2();
659 const bool useP2M2 = !(
m2 > 0) &&
661 !(std::abs(
m2) < 2.0*DBL_EPSILON*ene*ene);
663 const short recoMethod = ( !useP2M2 ? 0: ( ene >= 0. ? 1 : 2 ) );
664 auto A_theta=
p->attribute<HepMC3::DoubleAttribute>(
"theta");
665 auto A_phi=
p->attribute<HepMC3::DoubleAttribute>(
"phi");
666 auto A_flows=
p->attribute<HepMC3::VectorIntAttribute>(
"flows");
675 A_flows?(A_flows->value().size()):0,
676 A_theta?(A_theta->value()):0.0,
677 A_phi?(A_phi->value()):0.0,
683 std::vector< std::pair<int,int> > flow_hepmc2;
684 if(A_flows) flow_hepmc2=vector_to_vector_int_int(A_flows->value());
685 persEvt.
m_genParticles.back().m_flow.assign( flow_hepmc2.begin(),flow_hepmc2.end() );
694 const HepMC::FourVector&
mom =
p.m_momentum;
695 const double ene =
mom.e();
696 const double m2 =
mom.m2();
699 const bool useP2M2 = !(
m2 > 0) &&
701 !(std::abs(
m2) < 2.0*DBL_EPSILON*ene*ene);
703 const short recoMethod = ( !useP2M2
717 p.m_polarization.theta(),
718 p.m_polarization.phi(),
719 p.m_production_vertex
720 ?
p.m_production_vertex->barcode()
723 ?
p.m_end_vertex->barcode()