7 #include "AthLinks/ElementLink.h"
9 #include "GaudiKernel/MsgStream.h"
10 #include "GaudiKernel/DataSvc.h"
11 #include "GaudiKernel/PhysicalConstants.h"
54 xAODTruthCnvAlg::xAODTruthCnvAlg(
const string&
name, ISvcLocator* svcLoc )
56 , m_metaStore(
"MetaDataStore",
name )
57 , m_firstBeginRun(true)
67 ATH_MSG_FATAL(
"Contradictory xAOD truth pile-up setting: all pile-up AND in-time alone requested simultaneously. Check settings." );
68 return StatusCode::FAILURE;
92 incSvc->addListener(
this,
"BeginRun", 10);
107 return StatusCode::SUCCESS;
114 ATH_CHECK(truthLinkVec.
record(std::make_unique<xAODTruthParticleLinkVector>()));
121 return StatusCode::FAILURE;
131 ATH_CHECK(xTruthEventContainer.
record(std::make_unique<xAOD::TruthEventContainer>(),
132 std::make_unique<xAOD::TruthEventAuxContainer>()));
139 ATH_CHECK(xTruthPileupEventContainer.
record(std::make_unique<xAOD::TruthPileupEventContainer>(),
140 std::make_unique<xAOD::TruthPileupEventAuxContainer>()));
146 ATH_CHECK(xTruthParticleContainer.
record(std::make_unique<xAOD::TruthParticleContainer>(),
147 std::make_unique<xAOD::TruthParticleAuxContainer>()));
151 ATH_CHECK(xTruthVertexContainer.
record(std::make_unique<xAOD::TruthVertexContainer>(),
152 std::make_unique<xAOD::TruthVertexAuxContainer>()));
157 bool hadLHERecord =
false;
193 bool newAttributesPresent(
false);
195 for (
unsigned int cntr = 0; cntr < mcColl->
size(); ++cntr) {
196 const HepMC::GenEvent* genEvt = (*mcColl)[cntr];
197 bool isSignalProcess(
false);
199 isSignalProcess=
true;
201 auto bunchCrossingTime = genEvt->attribute<HepMC3::IntAttribute>(
"BunchCrossingTime");
202 if (bunchCrossingTime) {
203 newAttributesPresent =
true;
216 isSignalProcess=
false;
218 auto bunchCrossingTime = genEvt->attribute<HepMC3::IntAttribute>(
"BunchCrossingTime");
219 if (bunchCrossingTime) {
233 if (newAttributesPresent) {
261 if (isSignalProcess) {
262 xTruthEventContainer->
push_back( xTruthEvent );
283 return StatusCode::FAILURE;
290 for (
const double&
w : genEvt->weights())
weights.push_back((
float)(
w));
296 auto const hiInfo = genEvt->heavy_ion();
334 auto const pdfInfo = genEvt->pdf_info();
363 auto lhe_record_attribute = genEvt->attribute<HepMC::ShortEventAttribute>(
"LHERecord");
369 ATH_CHECK(xTruthLHEParticleContainer.
record(std::make_unique<xAOD::TruthParticleContainer>(),
370 std::make_unique<xAOD::TruthParticleAuxContainer>()));
373 for (
int nPart=0;nPart<lhe_record_attribute->NUP;++nPart) {
377 xTruthLHEParticleContainer->
push_back( xTruthParticle );
379 xTruthParticle->
setPdgId( lhe_record_attribute->IDUP[nPart] );
381 xTruthParticle->
setStatus( lhe_record_attribute->ISTUP[nPart] );
382 xTruthParticle->
setPx( lhe_record_attribute->PUP[nPart][0] );
383 xTruthParticle->
setPy( lhe_record_attribute->PUP[nPart][1] );
384 xTruthParticle->
setPz( lhe_record_attribute->PUP[nPart][2] );
385 xTruthParticle->
setE( lhe_record_attribute->PUP[nPart][3] );
386 xTruthParticle->
setM( lhe_record_attribute->PUP[nPart][4] );
389 else if (hadLHERecord){
390 ATH_MSG_WARNING(
"Truth record appeared to have two LHE records; this should not be possible");
394 ATH_MSG_WARNING(
"HEPMC2 does not support LHE truth record storage. Skipping.");
398 xTruthPileupEventContainer->
push_back( xTruthPileupEvent );
407 vector<HepMC::ConstGenVertexPtr> vertices;
412 if (disconnectedSignalProcessVtx) {
413 if (disconnectedSignalProcessVtx->particles_in_size() == 0 && disconnectedSignalProcessVtx->particles_out_size() == 0 ) {
415 vertices.push_back (disconnectedSignalProcessVtx);
418 ATH_MSG_WARNING(
"Signal process vertex pointer not valid in HepMC Collection for GenEvent #" << cntr <<
" / " << mcColl->
size());
422 pair<HepMC::ConstGenParticlePtr,HepMC::ConstGenParticlePtr> beamParticles;
423 bool genEvt_valid_beam_particles=
false;
425 auto beamParticles_vec = genEvt->beams();
426 genEvt_valid_beam_particles=(beamParticles_vec.size()>1);
427 if (genEvt_valid_beam_particles){beamParticles.first=beamParticles_vec[0]; beamParticles.second=beamParticles_vec[1]; }
429 auto bcmapatt = genEvt->attribute<HepMC::GenEventBarcodes>(
"barcodes");
430 if (!bcmapatt)
ATH_MSG_ERROR(
"TruthParticleCnvTool.cxx: Event does not contain barcodes attribute");
431 std::map<int, HepMC3::ConstGenParticlePtr> bcmap = bcmapatt->barcode_to_particle_map();
432 xTruthParticleContainer->
reserve(bcmap.size());
433 for (
const auto &[genPartBarcode,
part]: bcmap) {
435 genEvt_valid_beam_particles=genEvt->valid_beam_particles();
436 if ( genEvt_valid_beam_particles ) beamParticles = genEvt->beam_particles();
437 xTruthParticleContainer->
reserve(genEvt->particles_size());
438 for (
auto part: *genEvt) {
444 xTruthParticleContainer->
push_back( xTruthParticle );
456 if (genEvt_valid_beam_particles) {
457 if (isSignalProcess) {
463 auto productionVertex =
part->production_vertex();
467 if (productionVertex && productionVertex->parent_event() !=
nullptr) {
469 if (
parts.incoming.empty() &&
parts.outgoing.empty())
470 vertices.push_back (productionVertex);
471 parts.outgoingEL.push_back(eltp);
472 parts.outgoing.push_back(xTruthParticle);
478 auto decayVertex =
part->end_vertex();
481 if (
parts.incoming.empty() &&
parts.outgoing.empty())
482 vertices.push_back (decayVertex);
483 parts.incomingEL.push_back(eltp);
484 parts.incoming.push_back(xTruthParticle);
491 xTruthVertexContainer->
reserve(vertices.size());
492 for (
const auto&
vertex : vertices) {
497 xTruthVertexContainer->
push_back( xTruthVertex );
516 if (isSignalProcess)
delete xTruthPileupEvent;
517 if (!isSignalProcess)
delete xTruthEvent;
525 return StatusCode::SUCCESS;
532 if(inputStore.retrieve().isFailure()) {
537 if(inputStore->retrieve(tagInfo,
"/TagInfo").isFailure()) {
538 ATH_MSG_WARNING(
"Failed to retrieve /TagInfo metadata from the input store");
541 if(tagInfo->payloadContainer()->size()>0) {
543 if(tagInfoPayload->
size()>0) {
545 if (al.exists(
"lhefGenerator")){
549 if (al.exists(
"generators")){
553 if (al.exists(
"evgenProcess")){
557 if (al.exists(
"evgenTune")){
561 if (al.exists(
"hardPDF")){
565 if (al.exists(
"softPDF")){
582 tv->
setX(gv->position().x());
583 tv->
setY(gv->position().y());
584 tv->
setZ(gv->position().z());
585 tv->
setT(gv->position().t());
591 tp->setPdgId(gp->pdg_id());
596 if (pol.is_defined()) {
601 tp->setM(gp->generated_mass());
602 tp->setPx(gp->momentum().px());
603 tp->setPy(gp->momentum().py());
604 tp->setPz(gp->momentum().pz());
605 tp->setE(gp->momentum().e());
611 const std::string& metaName)
615 auto md = std::make_unique<xAOD::TruthMetaDataContainer>();
618 auto aux = std::make_unique<xAOD::TruthMetaDataAuxContainer>();
619 md->setStore( aux.get() );
622 CHECK( metaStore->record( std::move (aux), metaName +
"Aux." ) );
623 CHECK( metaStore->record( std::move (md), metaName ) );
624 return StatusCode::SUCCESS;
630 const HepMC::GenEvent& genEvt,
639 m_tmd->push_back (std::make_unique <xAOD::TruthMetaData>());
645 std::vector<std::string> orderedWeightNameVec;
646 if (!genEvt.run_info()) {
649 if (!genEvt.run_info()->weight_names().empty()) {
650 orderedWeightNameVec = genEvt.weight_names();
661 const auto& weightNameMap = genEvt.weights().m_names;
662 std::vector<std::string> orderedWeightNameVec;
663 orderedWeightNameVec.reserve( weightNameMap.size() );
664 for (
const auto&
entry: weightNameMap) {
665 orderedWeightNameVec.push_back(
entry.first);
670 std::sort(orderedWeightNameVec.begin(), orderedWeightNameVec.end(),
671 [&](
const std::string&
i,
const std::string& j){return weightNameMap.at(i) < weightNameMap.at(j);});
689 if(!metaFields.
hardPDF.empty()) {
692 if(!metaFields.
softPDF.empty()) {
697 return StatusCode::SUCCESS;