Input/src/HepMCFileReader.cc

0001
0002 /**
0003 *  See header file for a description of this class.
0004 *
0005 *
0006 *  \author Jo. Weng  - CERN, Ph Division & Uni Karlsruhe
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   // Implement a HepMCFileReader singleton.
0030
0031   if (instance_ == nullptr) {
0032     instance_ = new HepMCFileReader();
0033   }
0034   return instance_;
0035 }
0036
0037 //-------------------------------------------------------------------------
0038 HepMCFileReader::HepMCFileReader() : evt_(nullptr), input_(nullptr) {
0039   // Default constructor.
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   // work around a HepMC IO_ inheritence shortfall
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     //  printHepMcEvent();
0089     //          printEvent();
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 // Print out in old CMKIN style for comparisons
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       // tau=L/(gamma*beta*c)
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         // convert lab time to proper time
0149         double properTime = labTime / g->momentum().rho() * (g->generatedMass());
0150         // set the proper time in nanoseconds
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     // making a list of incoming particles of the vertices
0172     // so that the mother indices in HEPEVT can be filled properly
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       //particle_counter can be very large for heavy ions
0178       if (particle_counter >= index_to_particle.size()) {
0179         //make it large enough to hold up to this index
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     // daughters are entered only if they aren't a mother of
0186     // another vertex
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         //particle_counter can be very large for heavy ions
0193         if (particle_counter >= index_to_particle.size()) {
0194           //make it large enough to hold up to this index
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       //HepLorentzVector p = index_to_particle[j]->
0213       //production_vertex()->position();
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       //find # of 1. daughter
0229       int first_daughter =
0230           find_in_map(particle_to_index, *(index_to_particle[j]->end_vertex()->particles_begin(HepMC::children)));
0231       //cout <<"first_daughter "<< first_daughter <<  "num_daughters " << num_daughters << endl;
0232       HepMC::GenVertex::particle_iterator ic;
0233       int last_daughter = 0;
0234       //find # of last daughter
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 }