File indexing completed on 2024-04-06 12:13:50
0001 #include <algorithm>
0002 #include <iomanip>
0003 #include <iostream>
0004 #include <memory>
0005
0006 #include <cmath>
0007 #include <cstring>
0008 #include <memory>
0009 #include <sstream>
0010 #include <string>
0011 #include <vector>
0012
0013 #include "HepMC/GenEvent.h"
0014 #include "HepMC/GenVertex.h"
0015 #include "HepMC/GenParticle.h"
0016 #include "HepMC/SimpleVector.h"
0017
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019
0020 #include "SimDataFormats/GeneratorProducts/interface/LesHouches.h"
0021
0022 #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
0023 #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
0024
0025 static int skipWhitespace(std::istream &in) {
0026 int ch;
0027 do {
0028 ch = in.get();
0029 } while (std::isspace(ch));
0030 if (ch != std::istream::traits_type::eof())
0031 in.putback(ch);
0032 return ch;
0033 }
0034
0035 namespace lhef {
0036
0037 LHEEvent::LHEEvent(const std::shared_ptr<LHERunInfo> &runInfo, std::istream &in)
0038 : runInfo(runInfo),
0039 weights_(0),
0040 counted(false),
0041 readAttemptCounter(0),
0042 npLO_(-99),
0043 npNLO_(-99),
0044 evtnum_(-1)
0045
0046 {
0047 hepeup.NUP = 0;
0048 hepeup.XPDWUP.first = hepeup.XPDWUP.second = 0.0;
0049
0050 in >> hepeup.NUP >> hepeup.IDPRUP >> hepeup.XWGTUP >> hepeup.SCALUP >> hepeup.AQEDUP >> hepeup.AQCDUP;
0051 if (!in.good())
0052 throw cms::Exception("InvalidFormat") << "Les Houches file contained invalid"
0053 " event header."
0054 << std::endl;
0055
0056
0057 originalXWGTUP_ = hepeup.XWGTUP;
0058
0059 int idwtup = runInfo->getHEPRUP()->IDWTUP;
0060 if (idwtup >= 0 && hepeup.XWGTUP < 0) {
0061 edm::LogWarning("Generator|LHEInterface") << "Non-allowed negative event weight encountered." << std::endl;
0062 hepeup.XWGTUP = std::abs(hepeup.XWGTUP);
0063 }
0064
0065 if (std::abs(idwtup) == 3 && std::abs(hepeup.XWGTUP) != 1.) {
0066 edm::LogInfo("Generator|LHEInterface") << "Event weight not set to one for abs(IDWTUP) == 3" << std::endl;
0067 hepeup.XWGTUP = hepeup.XWGTUP > 0. ? 1.0 : -1.0;
0068 }
0069
0070 hepeup.resize();
0071
0072 for (int i = 0; i < hepeup.NUP; i++) {
0073 in >> hepeup.IDUP[i] >> hepeup.ISTUP[i] >> hepeup.MOTHUP[i].first >> hepeup.MOTHUP[i].second >>
0074 hepeup.ICOLUP[i].first >> hepeup.ICOLUP[i].second >> hepeup.PUP[i][0] >> hepeup.PUP[i][1] >>
0075 hepeup.PUP[i][2] >> hepeup.PUP[i][3] >> hepeup.PUP[i][4] >> hepeup.VTIMUP[i] >> hepeup.SPINUP[i];
0076 if (!in.good())
0077 throw cms::Exception("InvalidFormat") << "Les Houches file contained invalid event"
0078 " in particle line "
0079 << (i + 1) << "." << std::endl;
0080 }
0081
0082 while (skipWhitespace(in) == '#') {
0083 std::string line;
0084 std::getline(in, line);
0085 std::istringstream ss(line);
0086 std::string tag;
0087 ss >> tag;
0088 if (tag == "#pdf") {
0089 pdf = std::make_unique<PDF>();
0090 ss >> pdf->id.first >> pdf->id.second >> pdf->x.first >> pdf->x.second >> pdf->scalePDF >> pdf->xPDF.first >>
0091 pdf->xPDF.second;
0092 if (ss.bad()) {
0093 edm::LogWarning("Generator|LHEInterface") << "Les Houches event contained"
0094 " unparseable PDF information."
0095 << std::endl;
0096 pdf.reset();
0097 } else
0098 continue;
0099 }
0100 comments.push_back(line + "\n");
0101 }
0102
0103 if (!in.eof())
0104 edm::LogWarning("Generator|LHEInterface") << "Les Houches file contained spurious"
0105 " content after event data."
0106 << std::endl;
0107 }
0108
0109 LHEEvent::LHEEvent(const std::shared_ptr<LHERunInfo> &runInfo, const HEPEUP &hepeup)
0110 : runInfo(runInfo), hepeup(hepeup), counted(false), readAttemptCounter(0), npLO_(-99), npNLO_(-99), evtnum_(-1) {}
0111
0112 LHEEvent::LHEEvent(const std::shared_ptr<LHERunInfo> &runInfo,
0113 const HEPEUP &hepeup,
0114 const LHEEventProduct::PDF *pdf,
0115 const std::vector<std::string> &comments)
0116 : runInfo(runInfo),
0117 hepeup(hepeup),
0118 pdf(pdf ? new PDF(*pdf) : nullptr),
0119 comments(comments),
0120 counted(false),
0121 readAttemptCounter(0),
0122 npLO_(-99),
0123 npNLO_(-99),
0124 evtnum_(-1) {}
0125
0126 LHEEvent::LHEEvent(const std::shared_ptr<LHERunInfo> &runInfo, const LHEEventProduct &product)
0127 : runInfo(runInfo),
0128 hepeup(product.hepeup()),
0129 pdf(product.pdf() ? new PDF(*product.pdf()) : nullptr),
0130 weights_(product.weights()),
0131 comments(product.comments_begin(), product.comments_end()),
0132 counted(false),
0133 readAttemptCounter(0),
0134 originalXWGTUP_(product.originalXWGTUP()),
0135 scales_(product.scales()),
0136 npLO_(product.npLO()),
0137 npNLO_(product.npNLO()),
0138 evtnum_(product.evtnum()) {}
0139
0140 LHEEvent::~LHEEvent() {}
0141
0142 void LHEEvent::removeParticle(HEPEUP &hepeup, int index) {
0143 index--;
0144 if (index < 0 || index >= hepeup.NUP) {
0145 edm::LogError("Generator|LHEInterface")
0146 << "removeParticle: Index " << (index + 1) << " out of bounds." << std::endl;
0147 return;
0148 }
0149
0150 std::pair<int, int> mo = hepeup.MOTHUP[index];
0151
0152 hepeup.IDUP.erase(hepeup.IDUP.begin() + index);
0153 hepeup.ISTUP.erase(hepeup.ISTUP.begin() + index);
0154 hepeup.MOTHUP.erase(hepeup.MOTHUP.begin() + index);
0155 hepeup.ICOLUP.erase(hepeup.ICOLUP.begin() + index);
0156 hepeup.PUP.erase(hepeup.PUP.begin() + index);
0157 hepeup.VTIMUP.erase(hepeup.VTIMUP.begin() + index);
0158 hepeup.SPINUP.erase(hepeup.SPINUP.begin() + index);
0159
0160 hepeup.NUP--;
0161
0162 index++;
0163 for (int i = 0; i < hepeup.NUP; i++) {
0164 if (hepeup.MOTHUP[i].first == index) {
0165 if (hepeup.MOTHUP[i].second > 0)
0166 edm::LogError("Generator|LHEInterface")
0167 << "removeParticle: Particle " << (i + 2) << " has two mothers." << std::endl;
0168 hepeup.MOTHUP[i] = mo;
0169 }
0170
0171 if (hepeup.MOTHUP[i].first > index)
0172 hepeup.MOTHUP[i].first--;
0173 if (hepeup.MOTHUP[i].second > index)
0174 hepeup.MOTHUP[i].second--;
0175 }
0176 }
0177
0178 void LHEEvent::removeResonances(const std::vector<int> &ids) {
0179 for (int i = 1; i <= hepeup.NUP; i++) {
0180 int id = std::abs(hepeup.IDUP[i - 1]);
0181 if (hepeup.MOTHUP[i - 1].first > 0 && hepeup.MOTHUP[i - 1].second > 0 &&
0182 std::find(ids.begin(), ids.end(), id) != ids.end())
0183 removeParticle(hepeup, i--);
0184 }
0185 }
0186
0187 void LHEEvent::count(LHERunInfo::CountMode mode, double weight, double matchWeight) {
0188 if (counted)
0189 edm::LogWarning("Generator|LHEInterface") << "LHEEvent::count called twice on same event!" << std::endl;
0190
0191 runInfo->count(hepeup.IDPRUP, mode, hepeup.XWGTUP, weight, matchWeight);
0192
0193 counted = true;
0194 }
0195
0196 void LHEEvent::fillPdfInfo(HepMC::PdfInfo *info) const {
0197 if (pdf.get()) {
0198 info->set_id1(pdf->id.first);
0199 info->set_id2(pdf->id.second);
0200 info->set_x1(pdf->x.first);
0201 info->set_x2(pdf->x.second);
0202 info->set_pdf1(pdf->xPDF.first);
0203 info->set_pdf2(pdf->xPDF.second);
0204 info->set_scalePDF(pdf->scalePDF);
0205 } else if (hepeup.NUP >= 2) {
0206 const HEPRUP *heprup = getHEPRUP();
0207 info->set_id1(hepeup.IDUP[0] == 21 ? 0 : hepeup.IDUP[0]);
0208 info->set_id2(hepeup.IDUP[1] == 21 ? 0 : hepeup.IDUP[1]);
0209 info->set_x1(std::abs(hepeup.PUP[0][2] / heprup->EBMUP.first));
0210 info->set_x2(std::abs(hepeup.PUP[1][2] / heprup->EBMUP.second));
0211 info->set_pdf1(-1.0);
0212 info->set_pdf2(-1.0);
0213 info->set_scalePDF(hepeup.SCALUP);
0214 } else {
0215 info->set_x1(-1.0);
0216 info->set_x2(-1.0);
0217 info->set_pdf1(-1.0);
0218 info->set_pdf2(-1.0);
0219 info->set_scalePDF(hepeup.SCALUP);
0220 }
0221 }
0222
0223 void LHEEvent::fillEventInfo(HepMC::GenEvent *event) const {
0224 event->set_signal_process_id(hepeup.IDPRUP);
0225 event->set_event_scale(hepeup.SCALUP);
0226 event->set_alphaQED(hepeup.AQEDUP);
0227 event->set_alphaQCD(hepeup.AQCDUP);
0228 }
0229
0230 std::unique_ptr<HepMC::GenEvent> LHEEvent::asHepMCEvent() const {
0231 std::unique_ptr<HepMC::GenEvent> hepmc(new HepMC::GenEvent);
0232
0233 hepmc->set_signal_process_id(hepeup.IDPRUP);
0234 hepmc->set_event_scale(hepeup.SCALUP);
0235 hepmc->set_alphaQED(hepeup.AQEDUP);
0236 hepmc->set_alphaQCD(hepeup.AQCDUP);
0237
0238 unsigned int nup = hepeup.NUP;
0239
0240
0241 if (!nup) {
0242 edm::LogWarning("Generator|LHEInterface") << "Les Houches Event does not contain any partons. "
0243 << "Not much to convert.";
0244 return hepmc;
0245 }
0246
0247
0248 std::vector<HepMC::GenParticle *> genParticles;
0249 std::vector<HepMC::GenVertex *> genVertices;
0250
0251
0252 for (unsigned int i = 0; i < nup; i++)
0253 genParticles.push_back(makeHepMCParticle(i));
0254
0255
0256 for (unsigned int i = 0; i < nup; i++) {
0257 unsigned int mother1 = hepeup.MOTHUP.at(i).first;
0258 unsigned int mother2 = hepeup.MOTHUP.at(i).second;
0259 double cTau = hepeup.VTIMUP.at(i);
0260
0261
0262 if (mother1) {
0263 mother1--;
0264 if (mother2)
0265 mother2--;
0266 else
0267 mother2 = mother1;
0268
0269 HepMC::GenParticle *in_par = genParticles.at(mother1);
0270 HepMC::GenVertex *current_vtx = in_par->end_vertex();
0271
0272 if (!current_vtx) {
0273 current_vtx = new HepMC::GenVertex(HepMC::FourVector(0, 0, 0, cTau));
0274
0275
0276 genVertices.push_back(current_vtx);
0277 }
0278
0279 for (unsigned int j = mother1; j <= mother2; j++)
0280 if (!genParticles.at(j)->end_vertex())
0281 current_vtx->add_particle_in(genParticles.at(j));
0282
0283
0284 current_vtx->add_particle_out(genParticles.at(i));
0285 }
0286 }
0287
0288 checkHepMCTree(hepmc.get());
0289
0290
0291
0292
0293
0294 const HEPRUP *heprup = getHEPRUP();
0295
0296
0297 HepMC::GenParticle *b1 = new HepMC::GenParticle(
0298 HepMC::FourVector(0.0, 0.0, +heprup->EBMUP.first, heprup->EBMUP.first), heprup->IDBMUP.first);
0299 HepMC::GenParticle *b2 = new HepMC::GenParticle(
0300 HepMC::FourVector(0.0, 0.0, -heprup->EBMUP.second, heprup->EBMUP.second), heprup->IDBMUP.second);
0301 b1->set_status(3);
0302 b2->set_status(3);
0303
0304 HepMC::GenVertex *v1 = new HepMC::GenVertex();
0305 HepMC::GenVertex *v2 = new HepMC::GenVertex();
0306 v1->add_particle_in(b1);
0307 v2->add_particle_in(b2);
0308
0309 hepmc->add_vertex(v1);
0310 hepmc->add_vertex(v2);
0311 hepmc->set_beam_particles(b1, b2);
0312
0313
0314 if (genParticles.size() >= 2) {
0315 if (!genParticles.at(0)->production_vertex() && !genParticles.at(1)->production_vertex()) {
0316 v1->add_particle_out(genParticles.at(0));
0317 v2->add_particle_out(genParticles.at(1));
0318 } else
0319 edm::LogWarning("Generator|LHEInterface") << "Initial partons do already have a"
0320 " production vertex. "
0321 << std::endl
0322 << "Beam particles not connected.";
0323 } else
0324 edm::LogWarning("Generator|LHEInterface") << "Can't find any initial partons to be"
0325 " connected to the beam particles.";
0326
0327 for (std::vector<HepMC::GenVertex *>::const_iterator iter = genVertices.begin(); iter != genVertices.end(); ++iter)
0328 hepmc->add_vertex(*iter);
0329
0330
0331 for (unsigned int i = 0; i < nup; i++) {
0332 if (!genParticles.at(i)->parent_event()) {
0333 edm::LogWarning("Generator|LHEInterface") << "Not all LHE particles could be stored"
0334 " stored in the HepMC event. "
0335 << std::endl
0336 << "Check the mother-daughter relations"
0337 " in the given LHE input file.";
0338 break;
0339 }
0340 }
0341
0342 hepmc->set_signal_process_vertex(const_cast<HepMC::GenVertex *>(findSignalVertex(hepmc.get(), false)));
0343
0344 return hepmc;
0345 }
0346
0347 HepMC::GenParticle *LHEEvent::makeHepMCParticle(unsigned int i) const {
0348 HepMC::GenParticle *particle = new HepMC::GenParticle(
0349 HepMC::FourVector(hepeup.PUP.at(i)[0], hepeup.PUP.at(i)[1], hepeup.PUP.at(i)[2], hepeup.PUP.at(i)[3]),
0350 hepeup.IDUP.at(i));
0351
0352 int status = hepeup.ISTUP.at(i);
0353
0354 particle->set_generated_mass(hepeup.PUP.at(i)[4]);
0355 particle->set_status(status > 0 ? (status == 2 ? 3 : status) : 3);
0356
0357 return particle;
0358 }
0359
0360 bool LHEEvent::checkHepMCTree(const HepMC::GenEvent *event) {
0361 double px = 0, py = 0, pz = 0, E = 0;
0362
0363 for (HepMC::GenEvent::particle_const_iterator iter = event->particles_begin(); iter != event->particles_end();
0364 iter++) {
0365 int status = (*iter)->status();
0366 HepMC::FourVector fv = (*iter)->momentum();
0367
0368
0369 if (status == 3 && *iter != event->beam_particles().first && *iter != event->beam_particles().second) {
0370 px -= fv.px();
0371 py -= fv.py();
0372 pz -= fv.pz();
0373 E -= fv.e();
0374 }
0375
0376
0377 if (status == 1) {
0378 px += fv.px();
0379 py += fv.py();
0380 pz += fv.pz();
0381 E += fv.e();
0382 }
0383 }
0384
0385 if (px * px + py * py + pz * pz + E * E > 0.1) {
0386 edm::LogWarning("Generator|LHEInterface")
0387 << "Energy-momentum badly conserved. " << std::setprecision(3) << "sum p_i = [" << std::setw(7) << E << ", "
0388 << std::setw(7) << px << ", " << std::setw(7) << py << ", " << std::setw(7) << pz << "]";
0389
0390 return false;
0391 }
0392
0393 return true;
0394 }
0395
0396 const HepMC::GenVertex *LHEEvent::findSignalVertex(const HepMC::GenEvent *event, bool status3) {
0397 double largestMass2 = -9.0e-30;
0398 const HepMC::GenVertex *vertex = nullptr;
0399 for (HepMC::GenEvent::vertex_const_iterator iter = event->vertices_begin(); iter != event->vertices_end(); ++iter) {
0400 if ((*iter)->particles_in_size() < 2)
0401 continue;
0402 if ((*iter)->particles_out_size() < 1 ||
0403 ((*iter)->particles_out_size() == 1 &&
0404 (!(*(*iter)->particles_out_const_begin())->end_vertex() ||
0405 !(*(*iter)->particles_out_const_begin())->end_vertex()->particles_out_size())))
0406 continue;
0407
0408 double px = 0.0, py = 0.0, pz = 0.0, E = 0.0;
0409 bool hadStatus3 = false;
0410 for (HepMC::GenVertex::particles_out_const_iterator iter2 = (*iter)->particles_out_const_begin();
0411 iter2 != (*iter)->particles_out_const_end();
0412 ++iter2) {
0413 hadStatus3 = hadStatus3 || (*iter2)->status() == 3;
0414 px += (*iter2)->momentum().px();
0415 py += (*iter2)->momentum().py();
0416 pz += (*iter2)->momentum().pz();
0417 E += (*iter2)->momentum().e();
0418 }
0419 if (status3 && !hadStatus3)
0420 continue;
0421
0422 double mass2 = E * E - (px * px + py * py + pz * pz);
0423 if (mass2 > largestMass2) {
0424 vertex = *iter;
0425 largestMass2 = mass2;
0426 }
0427 }
0428
0429 return vertex;
0430 }
0431
0432 static void fixSubTree(HepMC::GenVertex *vertex,
0433 HepMC::FourVector &_time,
0434 std::set<const HepMC::GenVertex *> &visited) {
0435 HepMC::FourVector time = _time;
0436 HepMC::FourVector curTime = vertex->position();
0437 bool needsFixup = curTime.t() < time.t();
0438
0439 if (!visited.insert(vertex).second && !needsFixup)
0440 return;
0441
0442 if (needsFixup)
0443 vertex->set_position(time);
0444 else
0445 time = curTime;
0446
0447 for (HepMC::GenVertex::particles_out_const_iterator iter = vertex->particles_out_const_begin();
0448 iter != vertex->particles_out_const_end();
0449 ++iter) {
0450 HepMC::GenVertex *endVertex = (*iter)->end_vertex();
0451 if (endVertex)
0452 fixSubTree(endVertex, time, visited);
0453 }
0454 }
0455
0456 void LHEEvent::fixHepMCEventTimeOrdering(HepMC::GenEvent *event) {
0457 std::set<const HepMC::GenVertex *> visited;
0458 HepMC::FourVector zeroTime;
0459 for (HepMC::GenEvent::vertex_iterator iter = event->vertices_begin(); iter != event->vertices_end(); ++iter)
0460 fixSubTree(*iter, zeroTime, visited);
0461 }
0462
0463 }