LHEInterface/src/LHEEvent.cc

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