6#ifndef ATLASHEPMC_GENEVENT_H
7#define ATLASHEPMC_GENEVENT_H
9#include "HepMC3/GenEvent.h"
10#include "HepMC3/GenHeavyIon.h"
11#include "HepMC3/GenPdfInfo.h"
12#include "HepMC3/PrintStreams.h"
28#include <unordered_map>
31inline std::vector<HepMC3::GenParticlePtr>::const_iterator
begin(HepMC3::GenEvent& e) {
return e.particles().begin(); }
32inline std::vector<HepMC3::GenParticlePtr>::const_iterator
end(HepMC3::GenEvent& e) {
return e.particles().end(); }
33inline std::vector<HepMC3::ConstGenParticlePtr>::const_iterator
begin(
const HepMC3::GenEvent& e) {
return e.particles().begin(); }
34inline std::vector<HepMC3::ConstGenParticlePtr>::const_iterator
end(
const HepMC3::GenEvent& e) {
return e.particles().end(); }
38using Print = HepMC3::Print;
39using GenHeavyIon = HepMC3::GenHeavyIon;
40using GenEvent = HepMC3::GenEvent;
42class ShortEventAttribute :
public HepMC3::Attribute {
44 ShortEventAttribute(): HepMC3::Attribute(){}
45 ShortEventAttribute(
const HepMC3::GenEvent* e): HepMC3::Attribute(){ from_event(e); }
47 bool from_event(
const HepMC3::GenEvent* e){
48 NUP=
e->particles().size();
50 XWGTUP =
e->weights().size() ?
e->weights()[0] : 1.0;
52 IDPRUP = A_signal_process_id ? A_signal_process_id->value() : 0;
54 SCALUP = A_event_scale ? A_event_scale->value() : 0;
56 AQCDUP = A_alphaQCD ? A_alphaQCD->value() : 0;
58 AQEDUP = A_alphaQED ? A_alphaQED->value() : 0;
59 for (
int i = 0;
i < NUP; ++
i ){
60 const auto & thisParticle =
e->particles().at(i);
61 PUP[
i][0] = thisParticle->momentum().px();
62 PUP[
i][1] = thisParticle->momentum().py();
63 PUP[
i][2] = thisParticle->momentum().pz();
64 PUP[
i][3] = thisParticle->momentum().e();
65 PUP[
i][4] = thisParticle->momentum().m();
66 IDUP[
i] = thisParticle->pdg_id();
67 const auto pv = thisParticle->production_vertex();
68 const auto ev = thisParticle->end_vertex();
69 if (pv &&
ev ) ISTUP[
i] = 2;
70 if (pv && !
ev ) ISTUP[
i] = 1;
71 if (thisParticle->status() == 4 || !pv ) ISTUP[
i] = -1;
78 if (pv && !
pv->particles_in().empty()) {
79 l =
pv->particles_in().front()->id();
81 for (
const auto& p :
pv->particles_in()) {
82 const int id =
p->id();
93 bool from_string(
const std::string &att)
override {
94 std::istringstream iss(att);
104 for (
int i = 0;
i < NUP; ++
i ){
107 iss >> MOTHUP[
i].first;
108 iss >> MOTHUP[
i].second;
109 iss >> ICOLUP[
i].first;
110 iss >> ICOLUP[
i].second;
123 bool to_string(std::string &fl)
const override {
124 std::ostringstream
file;
125 file <<
" " << std::setw(4) << NUP
126 <<
" " << std::setw(6) << IDPRUP
127 <<
" " << std::setw(14) << XWGTUP
128 <<
" " << std::setw(14) << SCALUP
129 <<
" " << std::setw(14) << AQEDUP
130 <<
" " << std::setw(14) << AQCDUP <<
"\n";
131 for (
int i = 0;
i < NUP; ++
i )
132 file <<
" " << std::setw(8) << IDUP[
i]
133 <<
" " << std::setw(2) << ISTUP[
i]
134 <<
" " << std::setw(4) << MOTHUP[
i].first
135 <<
" " << std::setw(4) << MOTHUP[
i].second
136 <<
" " << std::setw(4) << ICOLUP[
i].first
137 <<
" " << std::setw(4) << ICOLUP[
i].second
138 <<
" " << std::setw(14) << PUP[
i][0]
139 <<
" " << std::setw(14) << PUP[
i][1]
140 <<
" " << std::setw(14) << PUP[
i][2]
141 <<
" " << std::setw(14) << PUP[
i][3]
142 <<
" " << std::setw(14) << PUP[
i][4]
143 <<
" " << std::setw(1) << VTIMUP[
i]
144 <<
" " << std::setw(1) << SPINUP[
i] << std::endl;
154 PUP.resize(NUP, std::vector<double>(5));
165 std::vector<long> IDUP{};
166 std::vector<int> ISTUP{};
167 std::vector< std::pair<int,int> > MOTHUP{};
168 std::vector< std::pair<int,int> > ICOLUP{};
169 std::vector< std::vector<double> > PUP{};
170 std::vector<double> VTIMUP{};
171 std::vector<double> SPINUP{};
175class GenEventBarcodes :
public HepMC3::Attribute
178 virtual bool from_string(
const std::string & )
override {
return false; }
179 virtual bool to_string(std::string &att)
const override
181 att =
"GenEventBarcodes";
185 const auto it = m_vertexBC.find (
id);
186 if (it != m_vertexBC.end())
return it->second;
190 const auto it = m_vertexBC.find (
id);
191 if (it != m_vertexBC.end())
return it->second;
195 const auto it = m_particleBC.find (
id);
196 if (it != m_particleBC.end())
return it->second;
200 const auto it = m_particleBC.find (
id);
201 if (it != m_particleBC.end())
return it->second;
205 void add (GenVertexPtr p) {
209 m_vertexBC[
barcode->value()] = std::move(p);
213 void remove (GenVertexPtr p) {
217 const auto it = m_vertexBC.find (
barcode->value());
218 if (it != m_vertexBC.end()) {
219 m_vertexBC.erase (it);
224 void add (GenParticlePtr p) {
228 m_particleBC[
barcode->value()] = std::move(p);
232 void remove (GenParticlePtr p) {
236 const auto it = m_particleBC.find (
barcode->value());
237 if (it != m_particleBC.end()) {
238 m_particleBC.erase (it);
243 std::map<int, ConstGenVertexPtr> barcode_to_vertex_map()
const {
244 std::map<int, ConstGenVertexPtr> ret;
245 for (
const auto &bcvertpair: m_vertexBC)
246 ret.insert({bcvertpair.first, std::const_pointer_cast<const HepMC3::GenVertex>(bcvertpair.second)});
249 std::map<int, ConstGenParticlePtr> barcode_to_particle_map()
const {
250 std::map<int, ConstGenParticlePtr> ret;
251 for (
const auto &bcpartpair: m_particleBC)
252 ret.insert({bcpartpair.first, std::const_pointer_cast<const HepMC3::GenParticle>(bcpartpair.second)});
255 std::map<int,int> id_to_barcode_map()
const {
256 std::map<int, int> ret;
257 for (
const auto &bcvertpair: m_vertexBC) ret.insert({bcvertpair.second->id(), bcvertpair.first});
258 for (
const auto &bcpartpair: m_particleBC) ret.insert({bcpartpair.second->id(), bcpartpair.first});
263 void fillAttribute(GenEvent* e) {
264 const auto eventAttributes =
e->attributes();
266 const bool hasBarcodeAttribute = barcodeAttributeIt != eventAttributes.end();
270 if (hasBarcodeAttribute && barcodeAttributeIt->second.count(i) > 0) {
271 const auto &
ptr = barcodeAttributeIt->second.at(i);
272 if (
ptr->is_parsed()) {
273 m_particleBC[
static_cast<HepMC3::IntAttribute*
>(
ptr.get())->value()] =
ptr->particle();
276 m_particleBC[std::atoi(
ptr->unparsed_string().c_str())] =
ptr->particle();
282 const auto &vertices =
e->vertices();
283 for (
size_t i = 1;
i <= vertices.size();
i++) {
284 if (hasBarcodeAttribute && barcodeAttributeIt->second.count(-i) > 0) {
285 const auto &
ptr = barcodeAttributeIt->second.at(-i);
286 if (
ptr->is_parsed()) {
287 m_vertexBC[
static_cast<HepMC3::IntAttribute*
>(
ptr.get())->value()] =
ptr->vertex();
290 m_vertexBC[std::atoi(
ptr->unparsed_string().c_str())] =
ptr->vertex();
293 m_vertexBC[-
i] = vertices[
i-1];
300 std::unordered_map<int, GenVertexPtr> m_vertexBC;
301 std::unordered_map<int, GenParticlePtr> m_particleBC;
317 : Attribute(), m_shadow(
p) {}
320 virtual bool from_string(
const std::string &)
override {
return false; }
323 virtual bool to_string(std::string &att)
const override
325 att =
"ShadowParticle";
344 return at ?
at->value() :
e->event_number();
347inline std::map<std::string, std::size_t> weights_map(
const HepMC3::GenEvent* e) {
348 std::map<std::string, std::size_t> ret;
349 auto run =
e->run_info();
350 if (!
run)
return ret;
351 const std::vector<std::string>
names =
run->weight_names();
352 for (
const auto& name: names) ret[
name] =
run->weight_index(name);
356inline std::vector<HepMC3::GenParticlePtr>::const_iterator
begin(HepMC3::GenEvent& e) {
return e.particles().begin(); }
357inline std::vector<HepMC3::GenParticlePtr>::const_iterator
end(HepMC3::GenEvent& e) {
return e.particles().end(); }
358inline std::vector<HepMC3::ConstGenParticlePtr>::const_iterator
begin(
const HepMC3::GenEvent& e) {
return e.particles().begin(); }
359inline std::vector<HepMC3::ConstGenParticlePtr>::const_iterator
end(
const HepMC3::GenEvent& e) {
return e.particles().end(); }
362 GenEvent*
e =
new GenEvent();
363 std::shared_ptr<HepMC3::IntAttribute> signal_process_id_A = std::make_shared<HepMC3::IntAttribute>(
signal_process_id);
366 e->set_event_number(event_number);
371 GenEvent*
e =
new GenEvent();
372 e->set_event_number(inEvt->event_number());
373 e->weights() = inEvt->weights();
374 auto a_mpi = inEvt->attribute<HepMC3::IntAttribute>(
HepMCStr::mpi);
375 if (a_mpi)
e->add_attribute(
HepMCStr::mpi, std::make_shared<HepMC3::IntAttribute>(*a_mpi));
379 if (a_event_scale)
e->add_attribute(
HepMCStr::event_scale, std::make_shared<HepMC3::DoubleAttribute>(*a_event_scale));
381 if (a_alphaQCD)
e->add_attribute(
HepMCStr::alphaQCD, std::make_shared<HepMC3::DoubleAttribute>(*a_alphaQCD));
383 if (a_alphaQED)
e->add_attribute(
HepMCStr::alphaQED, std::make_shared<HepMC3::DoubleAttribute>(*a_alphaQED));
384 auto a_pi = inEvt->pdf_info();
385 if (a_pi)
e->set_pdf_info(std::make_shared<HepMC3::GenPdfInfo>(*a_pi));
386 auto a_hi = inEvt->heavy_ion();
387 if (a_hi)
e->set_heavy_ion(std::make_shared<HepMC3::GenHeavyIon>(*a_hi));
389 if (a_random_states)
e->add_attribute(
HepMCStr::random_states, std::make_shared<HepMC3::VectorLongIntAttribute>(*a_random_states));
397 barcodes = std::make_shared<GenEventBarcodes>();
412 const auto eventAttributes =
e->attributes();
414 const bool hasBarcodeAttribute = barcodeAttributeIt != eventAttributes.end();
416 const auto &vertices =
e->vertices();
417 if (hasBarcodeAttribute) {
418 for (
size_t i = 1;
i <= vertices.size();
i++) {
419 const auto &ptrIt = barcodeAttributeIt->second.find(-i);
420 if (ptrIt != barcodeAttributeIt->second.end()) {
421 const auto &
ptr = ptrIt->second;
422 if (
ptr->is_parsed()) {
423 if (
id ==
static_cast<HepMC3::IntAttribute*
>(
ptr.get())->value()) {
424 return ptr->vertex();
428 if (
id == std::atoi(
ptr->unparsed_string().c_str())) {
429 return ptr->vertex();
436 if (-
id > 0 && -
id <=
static_cast<int>(vertices.size())) {
438 return vertices[-
id-1];
441 return HepMC3::ConstGenVertexPtr();
452 const auto eventAttributes =
e->attributes();
454 const bool hasBarcodeAttribute = barcodeAttributeIt != eventAttributes.end();
457 if (hasBarcodeAttribute) {
459 const auto &ptrIt = barcodeAttributeIt->second.find(i);
460 if (ptrIt != barcodeAttributeIt->second.end()) {
461 const auto &
ptr = ptrIt->second;
462 if (
ptr->is_parsed()) {
463 if (
id ==
static_cast<HepMC3::IntAttribute*
>(
ptr.get())->value()) {
464 return ptr->particle();
468 if (
id == std::atoi(
ptr->unparsed_string().c_str())) {
469 return ptr->particle();
476 if (
id > 0 &&
id <=
static_cast<int>(
particles.size())) {
481 return HepMC3::ConstGenParticlePtr();
492 const auto eventAttributes =
e->attributes();
494 const bool hasBarcodeAttribute = barcodeAttributeIt != eventAttributes.end();
496 const auto &vertices =
e->vertices();
497 if (hasBarcodeAttribute) {
498 for (
size_t i = 1;
i <= vertices.size();
i++) {
499 const auto &ptrIt = barcodeAttributeIt->second.find(-i);
500 if (ptrIt != barcodeAttributeIt->second.end()) {
501 const auto &
ptr = ptrIt->second;
502 if (
ptr->is_parsed()) {
503 if (
id ==
static_cast<HepMC3::IntAttribute*
>(
ptr.get())->value()) {
504 return ptr->vertex();
508 if (
id == std::atoi(
ptr->unparsed_string().c_str())) {
509 return ptr->vertex();
516 if (-
id > 0 && -
id <=
static_cast<int>(vertices.size())) {
518 return vertices[-
id-1];
521 return HepMC3::GenVertexPtr();
532 const auto eventAttributes =
e->attributes();
534 const bool hasBarcodeAttribute = barcodeAttributeIt != eventAttributes.end();
537 if (hasBarcodeAttribute) {
539 const auto &ptrIt = barcodeAttributeIt->second.find(i);
540 if (ptrIt != barcodeAttributeIt->second.end()) {
541 const auto &
ptr = ptrIt->second;
542 if (
ptr->is_parsed()) {
543 if (
id ==
static_cast<HepMC3::IntAttribute*
>(
ptr.get())->value()) {
544 return ptr->particle();
548 if (
id == std::atoi(
ptr->unparsed_string().c_str())) {
549 return ptr->particle();
556 if (
id > 0 &&
id <=
static_cast<int>(
particles.size())) {
561 return HepMC3::GenParticlePtr();
564inline int mpi(
const GenEvent& evt) {
565 std::shared_ptr<HepMC3::IntAttribute> A_mpi =
evt.attribute<HepMC3::IntAttribute>(
HepMCStr::mpi);
566 return A_mpi ? (A_mpi->value()) : 0;
568inline int mpi(
const GenEvent* evt) {
569 std::shared_ptr<HepMC3::IntAttribute> A_mpi =
evt->attribute<HepMC3::IntAttribute>(
HepMCStr::mpi);
570 return A_mpi ? (A_mpi->value()) : 0;
575 return A_signal_process_id ? (A_signal_process_id->value()) : 0;
579 return A_signal_process_id ? (A_signal_process_id->value()) : 0;
582 std::shared_ptr<HepMC3::IntAttribute>
signal_process_id = std::make_shared<HepMC3::IntAttribute>(i);
585inline void set_mpi(GenEvent* e,
const int i = 0) {
586 std::shared_ptr<HepMC3::IntAttribute>
mpi = std::make_shared<HepMC3::IntAttribute>(i);
593 if (!v || !e)
return;
601 if (!e)
return false;
603 for (
const auto& p :
e->beams()) {
if (
p->status() == 4) ++nBeams; }
604 if (nBeams != 2)
return false;
609 if (!
p->parent_event())
return false;
612 barcodes = std::make_shared<GenEventBarcodes>();
616 const bool ret =
p->add_attribute(
HepMCStr::barcode, std::make_shared<HepMC3::IntAttribute>(i));
632 if (num > std::numeric_limits<int>::max())
return false;
633 e->set_event_number((
int)num);
637 return e->event_number();
639inline GenEvent::particle_iterator
begin(HepMC::GenEvent& e) {
return e.particles_begin(); }
640inline GenEvent::particle_iterator
end(HepMC::GenEvent& e) {
return e.particles_end(); }
641inline GenEvent::particle_const_iterator
begin(
const HepMC::GenEvent& e) {
return e.particles_begin(); }
642inline GenEvent::particle_const_iterator
end(
const HepMC::GenEvent& e) {
return e.particles_end(); }
643inline GenEvent*
newGenEvent(
const int a,
const int b) {
return new GenEvent(
a,b); }
646inline GenVertex*
barcode_to_vertex(
const GenEvent* e,
int id) {
return e->barcode_to_vertex(
id);}
648inline int mpi(
const GenEvent& e) {
651inline int mpi(
const GenEvent* e) {
655 return e.signal_process_id();
658 return e->signal_process_id();
661 e->set_signal_process_id(i);
663inline void set_mpi(GenEvent* e,
const int i) {
667 e->set_random_states(
a);
670 e->set_signal_process_vertex(v);
673 HepMC::GenEvent* outEvt =
new HepMC::GenEvent( inEvt->signal_process_id(), inEvt->event_number() );
674 outEvt->set_mpi ( inEvt->mpi() );
675 outEvt->set_event_scale ( inEvt->event_scale() );
676 outEvt->set_alphaQCD ( inEvt->alphaQCD() );
677 outEvt->set_alphaQED ( inEvt->alphaQED() );
678 outEvt->weights() = inEvt->weights();
679 outEvt->set_random_states( inEvt->random_states() );
680 if (
nullptr != inEvt->heavy_ion() ) {
681 outEvt->set_heavy_ion ( *inEvt->heavy_ion() );
683 if (
nullptr != inEvt->pdf_info() ) {
684 outEvt->set_pdf_info ( *inEvt->pdf_info() );
695inline void line(std::ostream& os,
const GenEvent& e) {e.print(os);}
696inline void line(std::ostream& os,
const GenEvent* e) {e->print(os);}
697inline void content(std::ostream& os,
const GenEvent& e) {e.print(os);}
698inline void content(std::ostream& os,
const GenEvent* e) {e->print(os);}
bool add(const std::string &hname, TKey *tobj)
mapped_type at(key_type key) const
Look up an element in the map.
const std::string event_scale
const std::string barcode
const std::string long_long_event_number
const std::string signal_process_vertex
const std::string signal_process_id
const std::string random_states
const std::string alphaQED
const std::string alphaQCD
const std::string barcodes
const std::string ShadowParticle
void line(std::ostream &os, const GenEvent &e)
void content(std::ostream &os, const GenEvent &e)
int signal_process_id(const GenEvent &e)
void set_mpi(GenEvent *e, const int i)
bool set_ll_event_number(HepMC::GenEvent *e, long long int num)
void set_signal_process_vertex(GenEvent *e, T v)
GenParticle * barcode_to_particle(const GenEvent *e, int id)
GenEvent::particle_iterator begin(HepMC::GenEvent &e)
GenEvent::particle_iterator end(HepMC::GenEvent &e)
void set_random_states(GenEvent *e, std::vector< T > a)
GenVertex * barcode_to_vertex(const GenEvent *e, int id)
long long int get_ll_event_number(const HepMC::GenEvent *e)
HepMC::GenVertex * GenVertexPtr
bool suggest_barcode(T &p, int i)
void set_signal_process_id(GenEvent *e, const int i)
GenEvent * copyemptyGenEvent(const GenEvent *inEvt)
void fillBarcodesAttribute(GenEvent *)
GenParticle * GenParticlePtr
const GenParticle * ConstGenParticlePtr
GenVertex * signal_process_vertex(const GenEvent *e)
int mpi(const GenEvent &e)
bool valid_beam_particles(const GenEvent *e)
const HepMC::GenVertex * ConstGenVertexPtr
GenEvent * newGenEvent(const int a, const int b)
l
Printing final latex table to .tex output file.
virtual bool resize(size_t sz) override
Change the size of all aux data vectors.
int run(int argc, char *argv[])