6#ifndef ATLASHEPMC_GENEVENT_H
7#define ATLASHEPMC_GENEVENT_H
11#include "HepMC3/GenEvent.h"
12#include "HepMC3/GenHeavyIon.h"
13#include "HepMC3/GenPdfInfo.h"
14#include "HepMC3/PrintStreams.h"
20#include <unordered_map>
23inline std::vector<HepMC3::GenParticlePtr>::const_iterator
begin(HepMC3::GenEvent& e) {
return e.particles().begin(); }
24inline std::vector<HepMC3::GenParticlePtr>::const_iterator
end(HepMC3::GenEvent& e) {
return e.particles().end(); }
25inline std::vector<HepMC3::ConstGenParticlePtr>::const_iterator
begin(
const HepMC3::GenEvent& e) {
return e.particles().begin(); }
26inline std::vector<HepMC3::ConstGenParticlePtr>::const_iterator
end(
const HepMC3::GenEvent& e) {
return e.particles().end(); }
30using Print=HepMC3::Print;
31using GenHeavyIon=HepMC3::GenHeavyIon;
32using GenEvent=HepMC3::GenEvent;
34class ShortEventAttribute :
public HepMC3::Attribute {
36 ShortEventAttribute():HepMC3::Attribute(){}
37 ShortEventAttribute(
const HepMC3::GenEvent* e):HepMC3::Attribute(){ from_event(e); }
39 bool from_event(
const HepMC3::GenEvent* e){
40 NUP=
e->particles().size();
42 XWGTUP =
e->weights().size() ?
e->weights()[0] : 1.0;
43 auto A_signal_process_id=
e->attribute<HepMC3::IntAttribute>(
"signal_process_id");
44 IDPRUP = A_signal_process_id?A_signal_process_id->value() : 0;
45 auto A_event_scale=
e->attribute<HepMC3::DoubleAttribute>(
"event_scale");
46 SCALUP = A_event_scale? A_event_scale->value():0;
47 auto A_alphaQCD=
e->attribute<HepMC3::DoubleAttribute>(
"alphaQCD");
48 AQCDUP = A_alphaQCD? A_alphaQCD->value():0;
49 auto A_alphaQED=
e->attribute<HepMC3::DoubleAttribute>(
"alphaQED");
50 AQEDUP = A_alphaQED? A_alphaQED->value():0;
52 for (
int i = 0;
i < NUP; ++
i ){
53 PUP[
i][0] =
e->particles().at(i)->momentum().px();
54 PUP[
i][1] =
e->particles().at(i)->momentum().py();
55 PUP[
i][2] =
e->particles().at(i)->momentum().pz();
56 PUP[
i][3] =
e->particles().at(i)->momentum().e();
57 PUP[
i][4] =
e->particles().at(i)->momentum().m();
58 IDUP[
i] =
e->particles().at(i)->pdg_id();
59 auto pv =
e->particles().at(i)->production_vertex();
60 auto ev =
e->particles().at(i)->end_vertex();
61 if (pv &&
ev ) ISTUP[
i] = 2;
62 if (pv && !
ev ) ISTUP[
i] = 1;
63 if (
e->particles().at(i)->status() == 4 || !pv ) ISTUP[
i] = -1;
64 auto flow1 =
e->particles().at(i)->attribute<HepMC3::IntAttribute>(
"flow1");
65 auto flow2 =
e->particles().at(i)->attribute<HepMC3::IntAttribute>(
"flow2");
66 ICOLUP[
i].first = flow1 ? flow1->value() : 0;
67 ICOLUP[
i].second = flow2 ? flow2->value() : 0;
69 if (pv)
for (
auto p:
pv->particles_in()) {
l= std::min(l,
p->id());
h=std::max(h,
p->id());}
76 bool from_string(
const std::string &att)
override {
77 std::istringstream iss(att);
87 for (
int i = 0;
i < NUP; ++
i ){
90 iss >> MOTHUP[
i].first;
91 iss >> MOTHUP[
i].second;
92 iss >> ICOLUP[
i].first;
93 iss >> ICOLUP[
i].second;
106 bool to_string(std::string &fl)
const override {
107 std::ostringstream
file;
108 file <<
" " << std::setw(4) << NUP
109 <<
" " << std::setw(6) << IDPRUP
110 <<
" " << std::setw(14) << XWGTUP
111 <<
" " << std::setw(14) << SCALUP
112 <<
" " << std::setw(14) << AQEDUP
113 <<
" " << std::setw(14) << AQCDUP <<
"\n";
114 for (
int i = 0;
i < NUP; ++
i )
115 file <<
" " << std::setw(8) << IDUP[
i]
116 <<
" " << std::setw(2) << ISTUP[
i]
117 <<
" " << std::setw(4) << MOTHUP[
i].first
118 <<
" " << std::setw(4) << MOTHUP[
i].second
119 <<
" " << std::setw(4) << ICOLUP[
i].first
120 <<
" " << std::setw(4) << ICOLUP[
i].second
121 <<
" " << std::setw(14) << PUP[
i][0]
122 <<
" " << std::setw(14) << PUP[
i][1]
123 <<
" " << std::setw(14) << PUP[
i][2]
124 <<
" " << std::setw(14) << PUP[
i][3]
125 <<
" " << std::setw(14) << PUP[
i][4]
126 <<
" " << std::setw(1) << VTIMUP[
i]
127 <<
" " << std::setw(1) << SPINUP[
i] << std::endl;
137 PUP.resize(NUP, std::vector<double>(5));
148 std::vector<long> IDUP{};
149 std::vector<int> ISTUP{};
150 std::vector< std::pair<int,int> > MOTHUP{};
151 std::vector< std::pair<int,int> > ICOLUP{};
152 std::vector< std::vector<double> > PUP{};
153 std::vector<double> VTIMUP{};
154 std::vector<double> SPINUP{};
158class GenEventBarcodes :
public HepMC3::Attribute
161 virtual bool from_string(
const std::string & )
override {
return false; }
162 virtual bool to_string(std::string &att)
const override
164 att =
"GenEventBarcodes";
168 auto it = m_vertexBC.find (
id);
169 if (it != m_vertexBC.end())
return it->second;
173 auto it = m_vertexBC.find (
id);
174 if (it != m_vertexBC.end())
return it->second;
178 auto it = m_particleBC.find (
id);
179 if (it != m_particleBC.end())
return it->second;
183 auto it = m_particleBC.find (
id);
184 if (it != m_particleBC.end())
return it->second;
188 void add (GenVertexPtr p) {
190 auto barcode =
p->attribute<HepMC3::IntAttribute>(
"barcode");
192 m_vertexBC[
barcode->value()] = std::move(p);
196 void remove (GenVertexPtr p) {
198 auto barcode =
p->attribute<HepMC3::IntAttribute>(
"barcode");
200 auto it = m_vertexBC.find (
barcode->value());
201 if (it != m_vertexBC.end()) {
202 m_vertexBC.erase (it);
207 void add (GenParticlePtr p) {
209 auto barcode =
p->attribute<HepMC3::IntAttribute>(
"barcode");
211 m_particleBC[
barcode->value()] = std::move(p);
215 void remove (GenParticlePtr p) {
217 auto barcode =
p->attribute<HepMC3::IntAttribute>(
"barcode");
219 auto it = m_particleBC.find (
barcode->value());
220 if (it != m_particleBC.end()) {
221 m_particleBC.erase (it);
226 std::map<int, ConstGenVertexPtr> barcode_to_vertex_map()
const {
227 std::map<int, ConstGenVertexPtr> ret;
228 for (
const auto &bcvertpair: m_vertexBC)
229 ret.insert({bcvertpair.first,std::const_pointer_cast<const HepMC3::GenVertex>(bcvertpair.second)});
232 std::map<int, ConstGenParticlePtr> barcode_to_particle_map()
const {
233 std::map<int, ConstGenParticlePtr> ret;
234 for (
const auto &bcpartpair: m_particleBC)
235 ret.insert({bcpartpair.first,std::const_pointer_cast<const HepMC3::GenParticle>(bcpartpair.second)});
238 std::map<int,int> id_to_barcode_map()
const {
239 std::map<int, int> ret;
240 for (
const auto &bcvertpair: m_vertexBC) ret.insert({bcvertpair.second->id(),bcvertpair.first});
241 for (
const auto &bcpartpair: m_particleBC) ret.insert({bcpartpair.second->id(),bcpartpair.first});
246 void fillAttribute(GenEvent* e) {
247 const auto eventAttributes =
e->attributes();
248 const auto barcodeAttributeIt = eventAttributes.find(
"barcode");
249 const bool hasBarcodeAttribute = barcodeAttributeIt != eventAttributes.end();
253 if (hasBarcodeAttribute && barcodeAttributeIt->second.count(i) > 0) {
254 const auto &
ptr = barcodeAttributeIt->second.at(i);
255 if (
ptr->is_parsed()) {
256 m_particleBC[
static_cast<HepMC3::IntAttribute*
>(
ptr.get())->value()] =
ptr->particle();
259 m_particleBC[std::atoi(
ptr->unparsed_string().c_str())] =
ptr->particle();
265 const auto &vertices =
e->vertices();
266 for (
size_t i = 1;
i <= vertices.size();
i++) {
267 if (hasBarcodeAttribute && barcodeAttributeIt->second.count(-i) > 0) {
268 const auto &
ptr = barcodeAttributeIt->second.at(-i);
269 if (
ptr->is_parsed()) {
270 m_vertexBC[
static_cast<HepMC3::IntAttribute*
>(
ptr.get())->value()] =
ptr->vertex();
273 m_vertexBC[std::atoi(
ptr->unparsed_string().c_str())] =
ptr->vertex();
276 m_vertexBC[-
i] = vertices[
i-1];
283 std::unordered_map<int, GenVertexPtr> m_vertexBC;
284 std::unordered_map<int, GenParticlePtr> m_particleBC;
292class ShadowParticle :
public HepMC3::Attribute {
299 ShadowParticle(ConstGenParticlePtr p)
300 : Attribute(), m_shadow(
p) {}
303 virtual bool from_string(
const std::string &)
override {
return false; }
306 virtual bool to_string(std::string &att)
const override
308 att =
"ShadowParticle";
322 e->add_attribute(
"long_long_event_number", std::make_shared<HepMC3::LongLongAttribute>(num));
326 auto at =
e->attribute<HepMC3::LongLongAttribute>(
"long_long_event_number");
327 return at?at->value():
e->event_number();
330inline std::map<std::string, std::size_t> weights_map(
const HepMC3::GenEvent* e) {
331 std::map<std::string, std::size_t> ret;
332 auto run =
e->run_info();
333 if (!run)
return ret;
334 std::vector<std::string>
names =
run->weight_names();
335 for (
const auto& name: names) ret[
name] =
run->weight_index(name);
339inline std::vector<HepMC3::GenParticlePtr>::const_iterator
begin(HepMC3::GenEvent& e) {
return e.particles().begin(); }
340inline std::vector<HepMC3::GenParticlePtr>::const_iterator
end(HepMC3::GenEvent& e) {
return e.particles().end(); }
341inline std::vector<HepMC3::ConstGenParticlePtr>::const_iterator
begin(
const HepMC3::GenEvent& e) {
return e.particles().begin(); }
342inline std::vector<HepMC3::ConstGenParticlePtr>::const_iterator
end(
const HepMC3::GenEvent& e) {
return e.particles().end(); }
345 GenEvent*
e=
new GenEvent();
346 std::shared_ptr<HepMC3::IntAttribute> signal_process_id_A = std::make_shared<HepMC3::IntAttribute>(
signal_process_id);
347 e->add_attribute(
"signal_process_id",signal_process_id_A);
348 e->add_attribute(
"barcodes", std::make_shared<GenEventBarcodes>());
349 e->set_event_number(event_number);
354 GenEvent*
e=
new GenEvent();
355 e->set_event_number(inEvt->event_number());
356 e->weights()=inEvt->weights();
357 auto a_mpi = inEvt->attribute<HepMC3::IntAttribute>(
"mpi");
358 if (a_mpi)
e->add_attribute(
"mpi",std::make_shared<HepMC3::IntAttribute>(*a_mpi));
359 auto a_signal_process_id = inEvt->attribute<HepMC3::IntAttribute>(
"signal_process_id");
360 if (a_signal_process_id)
e->add_attribute(
"signal_process_id",std::make_shared<HepMC3::IntAttribute>(*a_signal_process_id));
361 auto a_event_scale = inEvt->attribute<HepMC3::DoubleAttribute>(
"event_scale");
362 if (a_event_scale)
e->add_attribute(
"event_scale",std::make_shared<HepMC3::DoubleAttribute>(*a_event_scale));
363 auto a_alphaQCD = inEvt->attribute<HepMC3::DoubleAttribute>(
"alphaQCD");
364 if (a_alphaQCD)
e->add_attribute(
"alphaQCD",std::make_shared<HepMC3::DoubleAttribute>(*a_alphaQCD));
365 auto a_alphaQED = inEvt->attribute<HepMC3::DoubleAttribute>(
"alphaQED");
366 if (a_alphaQED)
e->add_attribute(
"alphaQED",std::make_shared<HepMC3::DoubleAttribute>(*a_alphaQED));
367 auto a_pi = inEvt->pdf_info();
368 if (a_pi)
e->set_pdf_info(std::make_shared<HepMC3::GenPdfInfo>(*a_pi));
369 auto a_hi = inEvt->heavy_ion();
370 if (a_hi)
e->set_heavy_ion(std::make_shared<HepMC3::GenHeavyIon>(*a_hi));
371 auto a_random_states = inEvt->attribute<HepMC3::VectorLongIntAttribute>(
"random_states");
372 if (a_random_states)
e->add_attribute(
"random_states",std::make_shared<HepMC3::VectorLongIntAttribute>(*a_random_states));
373 e->add_attribute(
"barcodes", std::make_shared<GenEventBarcodes>());
378 auto barcodes =
e->attribute<GenEventBarcodes> (
"barcodes");
380 barcodes = std::make_shared<GenEventBarcodes>();
381 e->add_attribute(
"barcodes", barcodes);
389 const auto &
barcodes =
e->attribute<GenEventBarcodes> (
"barcodes");
395 const auto eventAttributes =
e->attributes();
396 const auto barcodeAttributeIt = eventAttributes.find(
"barcode");
397 const bool hasBarcodeAttribute = barcodeAttributeIt != eventAttributes.end();
399 const auto &vertices =
e->vertices();
400 if (hasBarcodeAttribute) {
401 for (
size_t i = 1;
i <= vertices.size();
i++) {
402 const auto &ptrIt = barcodeAttributeIt->second.find(-i);
403 if (ptrIt != barcodeAttributeIt->second.end()) {
404 const auto &
ptr = ptrIt->second;
405 if (
ptr->is_parsed()) {
406 if (
id ==
static_cast<HepMC3::IntAttribute*
>(
ptr.get())->value()) {
407 return ptr->vertex();
411 if (
id == std::atoi(
ptr->unparsed_string().c_str())) {
412 return ptr->vertex();
419 if (-
id > 0 && -
id <=
static_cast<int>(vertices.size())) {
420 if (!vertices[-
id-1]->attribute<HepMC3::IntAttribute>(
"barcode")) {
421 return vertices[-
id-1];
424 return HepMC3::ConstGenVertexPtr();
429 const auto &
barcodes =
e->attribute<GenEventBarcodes> (
"barcodes");
435 const auto eventAttributes =
e->attributes();
436 const auto barcodeAttributeIt = eventAttributes.find(
"barcode");
437 const bool hasBarcodeAttribute = barcodeAttributeIt != eventAttributes.end();
440 if (hasBarcodeAttribute) {
442 const auto &ptrIt = barcodeAttributeIt->second.find(i);
443 if (ptrIt != barcodeAttributeIt->second.end()) {
444 const auto &
ptr = ptrIt->second;
445 if (
ptr->is_parsed()) {
446 if (
id ==
static_cast<HepMC3::IntAttribute*
>(
ptr.get())->value()) {
447 return ptr->particle();
451 if (
id == std::atoi(
ptr->unparsed_string().c_str())) {
452 return ptr->particle();
459 if (
id > 0 &&
id <=
static_cast<int>(
particles.size())) {
460 if (!particles[
id-1]->attribute<HepMC3::IntAttribute>(
"barcode")) {
464 return HepMC3::ConstGenParticlePtr();
469 const auto &
barcodes =
e->attribute<GenEventBarcodes> (
"barcodes");
475 const auto eventAttributes =
e->attributes();
476 const auto barcodeAttributeIt = eventAttributes.find(
"barcode");
477 const bool hasBarcodeAttribute = barcodeAttributeIt != eventAttributes.end();
479 const auto &vertices =
e->vertices();
480 if (hasBarcodeAttribute) {
481 for (
size_t i = 1;
i <= vertices.size();
i++) {
482 const auto &ptrIt = barcodeAttributeIt->second.find(-i);
483 if (ptrIt != barcodeAttributeIt->second.end()) {
484 const auto &
ptr = ptrIt->second;
485 if (
ptr->is_parsed()) {
486 if (
id ==
static_cast<HepMC3::IntAttribute*
>(
ptr.get())->value()) {
487 return ptr->vertex();
491 if (
id == std::atoi(
ptr->unparsed_string().c_str())) {
492 return ptr->vertex();
499 if (-
id > 0 && -
id <=
static_cast<int>(vertices.size())) {
500 if (!vertices[-
id-1]->attribute<HepMC3::IntAttribute>(
"barcode")) {
501 return vertices[-
id-1];
504 return HepMC3::GenVertexPtr();
509 const auto &
barcodes =
e->attribute<GenEventBarcodes> (
"barcodes");
515 const auto eventAttributes =
e->attributes();
516 const auto barcodeAttributeIt = eventAttributes.find(
"barcode");
517 const bool hasBarcodeAttribute = barcodeAttributeIt != eventAttributes.end();
520 if (hasBarcodeAttribute) {
522 const auto &ptrIt = barcodeAttributeIt->second.find(i);
523 if (ptrIt != barcodeAttributeIt->second.end()) {
524 const auto &
ptr = ptrIt->second;
525 if (
ptr->is_parsed()) {
526 if (
id ==
static_cast<HepMC3::IntAttribute*
>(
ptr.get())->value()) {
527 return ptr->particle();
531 if (
id == std::atoi(
ptr->unparsed_string().c_str())) {
532 return ptr->particle();
539 if (
id > 0 &&
id <=
static_cast<int>(
particles.size())) {
540 if (!particles[
id-1]->attribute<HepMC3::IntAttribute>(
"barcode")) {
544 return HepMC3::GenParticlePtr();
547inline int mpi(
const GenEvent& evt) {
548 std::shared_ptr<HepMC3::IntAttribute> A_mpi=
evt.attribute<HepMC3::IntAttribute>(
"mpi");
549 return A_mpi?(A_mpi->value()):0;
551inline int mpi(
const GenEvent* evt) {
552 std::shared_ptr<HepMC3::IntAttribute> A_mpi=
evt->attribute<HepMC3::IntAttribute>(
"mpi");
553 return A_mpi?(A_mpi->value()):0;
557 std::shared_ptr<HepMC3::IntAttribute> A_signal_process_id=
evt.attribute<HepMC3::IntAttribute>(
"signal_process_id");
558 return A_signal_process_id?(A_signal_process_id->value()):0;
561 std::shared_ptr<HepMC3::IntAttribute> A_signal_process_id=
evt->attribute<HepMC3::IntAttribute>(
"signal_process_id");
562 return A_signal_process_id?(A_signal_process_id->value()):0;
565 std::shared_ptr<HepMC3::IntAttribute>
signal_process_id = std::make_shared<HepMC3::IntAttribute>(i);
568inline void set_mpi(GenEvent* e,
const int i=0) {
569 std::shared_ptr<HepMC3::IntAttribute>
mpi = std::make_shared<HepMC3::IntAttribute>(i);
570 e->add_attribute(
"mpi",std::move(
mpi));
573 e->add_attribute(
"random_states",std::make_shared<HepMC3::VectorLongIntAttribute>(
a));
576 if (!v || !e)
return;
579 v->add_attribute(
"signal_process_vertex",std::make_shared<HepMC3::IntAttribute>(1));
582inline GenVertexPtr signal_process_vertex(GenEvent* e) {
for (
auto v:
e->vertices())
if (
v->attribute<HepMC3::IntAttribute>(
"signal_process_vertex"))
return v;
return nullptr; }
584 if (!e)
return false;
586 for (
const auto& p :
e->beams()) {
if (
p->status() == 4) ++nBeams; }
587 if (nBeams != 2)
return false;
592 if (!
p->parent_event())
return false;
593 auto barcodes =
p->parent_event()->template attribute<GenEventBarcodes> (
"barcodes");
595 barcodes = std::make_shared<GenEventBarcodes>();
596 p->parent_event()->add_attribute(
"barcodes", barcodes);
599 bool ret =
p->add_attribute(
"barcode",std::make_shared<HepMC3::IntAttribute>(i));
615 if (num > std::numeric_limits<int>::max())
return false;
616 e->set_event_number((
int)num);
620 return e->event_number();
622inline GenEvent::particle_iterator
begin(HepMC::GenEvent& e) {
return e.particles_begin(); }
623inline GenEvent::particle_iterator
end(HepMC::GenEvent& e) {
return e.particles_end(); }
624inline GenEvent::particle_const_iterator
begin(
const HepMC::GenEvent& e) {
return e.particles_begin(); }
625inline GenEvent::particle_const_iterator
end(
const HepMC::GenEvent& e) {
return e.particles_end(); }
626inline GenEvent*
newGenEvent(
const int a,
const int b ) {
return new GenEvent(
a,b); }
629inline GenVertex*
barcode_to_vertex(
const GenEvent* e,
int id ) {
return e->barcode_to_vertex(
id);}
631inline int mpi(
const GenEvent& e) {
634inline int mpi(
const GenEvent* e) {
638 return e.signal_process_id();
641 return e->signal_process_id();
644 e->set_signal_process_id(i);
646inline void set_mpi(GenEvent* e,
const int i) {
650 e->set_random_states(
a);
653 e->set_signal_process_vertex(v);
656 HepMC::GenEvent* outEvt =
new HepMC::GenEvent( inEvt->signal_process_id(), inEvt->event_number() );
657 outEvt->set_mpi ( inEvt->mpi() );
658 outEvt->set_event_scale ( inEvt->event_scale() );
659 outEvt->set_alphaQCD ( inEvt->alphaQCD() );
660 outEvt->set_alphaQED ( inEvt->alphaQED() );
661 outEvt->weights() = inEvt->weights();
662 outEvt->set_random_states( inEvt->random_states() );
663 if (
nullptr != inEvt->heavy_ion() ) {
664 outEvt->set_heavy_ion ( *inEvt->heavy_ion() );
666 if (
nullptr != inEvt->pdf_info() ) {
667 outEvt->set_pdf_info ( *inEvt->pdf_info() );
678inline void line(std::ostream& os,
const GenEvent& e) {e.print(os);}
679inline void line(std::ostream& os,
const GenEvent* e) {e->print(os);}
680inline void content(std::ostream& os,
const GenEvent& e) {e.print(os);}
681inline void content(std::ostream& os,
const GenEvent* e) {e->print(os);}
bool add(const std::string &hname, TKey *tobj)
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.