File indexing completed on 2024-04-06 12:19:01
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include <iostream>
0010 #include <iomanip>
0011 #include <string>
0012
0013 #include "HepMC/GenEvent.h"
0014 #include "HepMC/GenParticle.h"
0015 #include "HepMC/IO_GenEvent.h"
0016
0017 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0018
0019 #include "IOMC/Input/interface/HepMCFileReader.h"
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021 #include "FWCore/Utilities/interface/Exception.h"
0022
0023 using namespace std;
0024
0025 HepMCFileReader *HepMCFileReader::instance_ = nullptr;
0026
0027
0028 HepMCFileReader *HepMCFileReader::instance() {
0029
0030
0031 if (instance_ == nullptr) {
0032 instance_ = new HepMCFileReader();
0033 }
0034 return instance_;
0035 }
0036
0037
0038 HepMCFileReader::HepMCFileReader() : evt_(nullptr), input_(nullptr) {
0039
0040 if (instance_ == nullptr) {
0041 instance_ = this;
0042 } else {
0043 edm::LogError("HepMCFileReader") << "Constructing a new instance";
0044 }
0045 }
0046
0047
0048 HepMCFileReader::~HepMCFileReader() {
0049 edm::LogInfo("HepMCFileReader") << "Destructing HepMCFileReader";
0050
0051 instance_ = nullptr;
0052 delete input_.get();
0053 }
0054
0055
0056 void HepMCFileReader::initialize(const string &filename) {
0057 if (isInitialized()) {
0058 edm::LogError("HepMCFileReader") << "Was already initialized... reinitializing";
0059 delete input_.get();
0060 }
0061
0062 edm::LogInfo("HepMCFileReader") << "Opening file" << filename << "using HepMC::IO_GenEvent";
0063 input_ = new HepMC::IO_GenEvent(filename, std::ios::in);
0064
0065 if (rdstate() == std::ios::failbit) {
0066 throw cms::Exception("FileNotFound", "HepMCFileReader::initialize()") << "File " << filename << " was not found.\n";
0067 }
0068 }
0069
0070
0071 int HepMCFileReader::rdstate() const {
0072
0073
0074 HepMC::IO_GenEvent const *p = dynamic_cast<HepMC::IO_GenEvent const *>(input());
0075 if (p)
0076 return p->rdstate();
0077
0078 return std::ios::failbit;
0079 }
0080
0081
0082 bool HepMCFileReader::readCurrentEvent() {
0083 evt_ = input_->read_next_event();
0084 if (evt_) {
0085 edm::LogInfo("HepMCFileReader") << "| --- Event Nr. " << evt_->event_number() << " with " << evt_->particles_size()
0086 << " particles --- !";
0087 ReadStats();
0088
0089
0090 } else {
0091 edm::LogInfo("HepMCFileReader") << "Got no event" << endl;
0092 }
0093
0094 return evt_ != nullptr;
0095 }
0096
0097
0098 bool HepMCFileReader::setEvent(int event) { return true; }
0099
0100
0101 bool HepMCFileReader::printHepMcEvent() const {
0102 if (evt_ != nullptr)
0103 evt_->print();
0104 return true;
0105 }
0106
0107
0108 HepMC::GenEvent *HepMCFileReader::fillCurrentEventData() {
0109 readCurrentEvent();
0110 return evt_;
0111 }
0112
0113
0114
0115 void HepMCFileReader::printEvent() const {
0116 int mo1 = 0, mo2 = 0, da1 = 0, da2 = 0, status = 0, pid = 0;
0117 if (evt_ != nullptr) {
0118 cout << "---#-------pid--st---Mo1---Mo2---Da1---Da2------px------py------pz-------E-";
0119 cout << "------m---------x---------y---------z---------t-";
0120 cout << endl;
0121 cout.setf(ios::right, ios::adjustfield);
0122 for (int n = 1; n <= evt_->particles_size(); n++) {
0123 HepMC::GenParticle *g = index_to_particle[n];
0124 getStatsFromTuple(mo1, mo2, da1, da2, status, pid, n);
0125 cout << setw(4) << n << setw(8) << pid << setw(5) << status << setw(6) << mo1 << setw(6) << mo2 << setw(6) << da1
0126 << setw(6) << da2;
0127 cout.setf(ios::fixed, ios::floatfield);
0128 cout.setf(ios::right, ios::adjustfield);
0129 cout << setw(10) << setprecision(2) << g->momentum().x();
0130 cout << setw(8) << setprecision(2) << g->momentum().y();
0131 cout << setw(10) << setprecision(2) << g->momentum().z();
0132 cout << setw(8) << setprecision(2) << g->momentum().t();
0133 cout << setw(8) << setprecision(2) << g->generatedMass();
0134
0135 if (g->production_vertex() != nullptr && g->end_vertex() != nullptr && status == 2) {
0136 cout << setw(10) << setprecision(2) << g->production_vertex()->position().x();
0137 cout << setw(10) << setprecision(2) << g->production_vertex()->position().y();
0138 cout << setw(10) << setprecision(2) << g->production_vertex()->position().z();
0139
0140 double xm = g->production_vertex()->position().x();
0141 double ym = g->production_vertex()->position().y();
0142 double zm = g->production_vertex()->position().z();
0143 double xd = g->end_vertex()->position().x();
0144 double yd = g->end_vertex()->position().y();
0145 double zd = g->end_vertex()->position().z();
0146 double decl = sqrt((xd - xm) * (xd - xm) + (yd - ym) * (yd - ym) + (zd - zm) * (zd - zm));
0147 double labTime = decl / c_light;
0148
0149 double properTime = labTime / g->momentum().rho() * (g->generatedMass());
0150
0151 cout << setw(8) << setprecision(2) << properTime;
0152 } else {
0153 cout << setw(10) << setprecision(2) << 0.0;
0154 cout << setw(10) << setprecision(2) << 0.0;
0155 cout << setw(10) << setprecision(2) << 0.0;
0156 cout << setw(8) << setprecision(2) << 0.0;
0157 }
0158 cout << endl;
0159 }
0160 } else {
0161 cout << " HepMCFileReader: No event available !" << endl;
0162 }
0163 }
0164
0165
0166 void HepMCFileReader::ReadStats() {
0167 unsigned int particle_counter = 0;
0168 index_to_particle.reserve(evt_->particles_size() + 1);
0169 index_to_particle[0] = nullptr;
0170 for (HepMC::GenEvent::vertex_const_iterator v = evt_->vertices_begin(); v != evt_->vertices_end(); ++v) {
0171
0172
0173 for (HepMC::GenVertex::particles_in_const_iterator p1 = (*v)->particles_in_const_begin();
0174 p1 != (*v)->particles_in_const_end();
0175 ++p1) {
0176 ++particle_counter;
0177
0178 if (particle_counter >= index_to_particle.size()) {
0179
0180 index_to_particle.resize(particle_counter + 1);
0181 }
0182 index_to_particle[particle_counter] = *p1;
0183 particle_to_index[*p1] = particle_counter;
0184 }
0185
0186
0187 for (HepMC::GenVertex::particles_out_const_iterator p2 = (*v)->particles_out_const_begin();
0188 p2 != (*v)->particles_out_const_end();
0189 ++p2) {
0190 if (!(*p2)->end_vertex()) {
0191 ++particle_counter;
0192
0193 if (particle_counter >= index_to_particle.size()) {
0194
0195 index_to_particle.resize(particle_counter + 1);
0196 }
0197 index_to_particle[particle_counter] = *p2;
0198 particle_to_index[*p2] = particle_counter;
0199 }
0200 }
0201 }
0202 }
0203
0204
0205 void HepMCFileReader::getStatsFromTuple(int &mo1, int &mo2, int &da1, int &da2, int &status, int &pid, int j) const {
0206 if (!evt_) {
0207 cout << "HepMCFileReader: Got no event :-( Game over already ?" << endl;
0208 } else {
0209 status = index_to_particle[j]->status();
0210 pid = index_to_particle[j]->pdg_id();
0211 if (index_to_particle[j]->production_vertex()) {
0212
0213
0214
0215 int num_mothers = index_to_particle[j]->production_vertex()->particles_in_size();
0216 int first_mother =
0217 find_in_map(particle_to_index, *(index_to_particle[j]->production_vertex()->particles_in_const_begin()));
0218 int last_mother = first_mother + num_mothers - 1;
0219 if (first_mother == 0)
0220 last_mother = 0;
0221 mo1 = first_mother;
0222 mo2 = last_mother;
0223 } else {
0224 mo1 = 0;
0225 mo2 = 0;
0226 }
0227 if (!index_to_particle[j]->end_vertex()) {
0228
0229 int first_daughter =
0230 find_in_map(particle_to_index, *(index_to_particle[j]->end_vertex()->particles_begin(HepMC::children)));
0231
0232 HepMC::GenVertex::particle_iterator ic;
0233 int last_daughter = 0;
0234
0235 for (ic = index_to_particle[j]->end_vertex()->particles_begin(HepMC::children);
0236 ic != index_to_particle[j]->end_vertex()->particles_end(HepMC::children);
0237 ++ic)
0238 last_daughter = find_in_map(particle_to_index, *ic);
0239
0240 if (first_daughter == 0)
0241 last_daughter = 0;
0242 da1 = first_daughter;
0243 da2 = last_daughter;
0244 } else {
0245 da1 = 0;
0246 da2 = 0;
0247 }
0248 }
0249 }
0250
0251
0252 int HepMCFileReader::find_in_map(const std::map<HepMC::GenParticle *, int> &m, HepMC::GenParticle *p) const {
0253 std::map<HepMC::GenParticle *, int>::const_iterator iter = m.find(p);
0254 return (iter == m.end()) ? 0 : iter->second;
0255 }