Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-05-02 05:10:04

0001 // -*- C++ -*-
0002 //
0003 // Package:    test/EmbeddingLHEProducer
0004 // Class:      EmbeddingLHEProducer
0005 //
0006 /**\class EmbeddingLHEProducer EmbeddingLHEProducer.cc test/EmbeddingLHEProducer/plugins/EmbeddingLHEProducer.cc
0007 
0008  Description: [one line class summary]
0009 
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Stefan Wayand
0015 //         Created:  Wed, 13 Jan 2016 08:15:01 GMT
0016 //
0017 //
0018 
0019 // system include files
0020 #include "TLorentzVector.h"
0021 #include <algorithm>
0022 #include <fstream>
0023 #include <iostream>
0024 #include <iterator>
0025 #include <memory>
0026 #include <string>
0027 
0028 // user include files
0029 #include "FWCore/Framework/interface/Frameworkfwd.h"
0030 #include "FWCore/Framework/interface/one/EDProducer.h"
0031 
0032 #include "FWCore/Framework/interface/Event.h"
0033 #include "FWCore/Framework/interface/MakerMacros.h"
0034 #include "FWCore/Framework/interface/Run.h"
0035 
0036 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0037 
0038 #include "DataFormats/Math/interface/LorentzVector.h"
0039 #include "DataFormats/PatCandidates/interface/Muon.h"
0040 #include "DataFormats/VertexReco/interface/Vertex.h"
0041 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0042 
0043 #include "SimDataFormats/GeneratorProducts/interface/LHECommonBlocks.h"
0044 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
0045 #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
0046 #include "SimDataFormats/GeneratorProducts/interface/LHEXMLStringProduct.h"
0047 #include "SimDataFormats/GeneratorProducts/interface/LesHouches.h"
0048 
0049 #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
0050 #include "GeneratorInterface/LHEInterface/interface/LHEReader.h"
0051 #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
0052 
0053 #include "CLHEP/Random/RandExponential.h"
0054 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0055 #include "FWCore/ServiceRegistry/interface/Service.h"
0056 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0057 #include "FWCore/Utilities/interface/StreamID.h"
0058 
0059 namespace CLHEP {
0060   class HepRandomEngine;
0061 }
0062 
0063 class EmbeddingLHEProducer : public edm::one::EDProducer<edm::BeginRunProducer, edm::EndRunProducer> {
0064 public:
0065   explicit EmbeddingLHEProducer(const edm::ParameterSet &);
0066   ~EmbeddingLHEProducer() override;
0067 
0068   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0069 
0070 private:
0071   void beginJob() override;
0072   void produce(edm::Event &, const edm::EventSetup &) override;
0073   void endJob() override;
0074 
0075   void beginRunProduce(edm::Run &run, edm::EventSetup const &es) override;
0076   void endRunProduce(edm::Run &, edm::EventSetup const &) override;
0077 
0078   void fill_lhe_from_mumu(TLorentzVector &positiveLepton,
0079                           TLorentzVector &negativeLepton,
0080                           lhef::HEPEUP &outlhe,
0081                           CLHEP::HepRandomEngine *engine);
0082   void fill_lhe_with_particle(lhef::HEPEUP &outlhe, TLorentzVector &particle, int pdgid, double spin, double ctau);
0083 
0084   void transform_mumu_to_tautau(TLorentzVector &positiveLepton, TLorentzVector &negativeLepton);
0085   const reco::Candidate *find_original_muon(const reco::Candidate *muon);
0086   void assign_4vector(TLorentzVector &Lepton, const pat::Muon *muon, std::string FSRmode);
0087   void mirror(TLorentzVector &positiveLepton, TLorentzVector &negativeLepton);
0088   void InitialRecoCorrection(TLorentzVector &positiveLepton, TLorentzVector &negativeLepton);
0089   void rotate180(TLorentzVector &positiveLepton, TLorentzVector &negativeLepton);
0090 
0091   LHERunInfoProduct::Header give_slha();
0092 
0093   edm::EDGetTokenT<edm::View<pat::Muon>> muonsCollection_;
0094   edm::EDGetTokenT<reco::VertexCollection> vertexCollection_;
0095   int particleToEmbed_;
0096   bool mirror_, rotate180_, InitialRecoCorrection_;
0097   static constexpr double tauMass_ = 1.77682;
0098   static constexpr double muonMass_ = 0.1057;
0099   static constexpr double elMass_ = 0.00051;
0100   const int embeddingParticles[3]{11, 13, 15};
0101 
0102   std::ofstream file;
0103   bool write_lheout;
0104 
0105   // instead of reconstruted 4vectors of muons generated are taken to study FSR. Possible modes:
0106   // afterFSR - muons without FSR photons taken into account
0107   // beforeFSR - muons with FSR photons taken into account
0108   std::string studyFSRmode_;
0109 };
0110 
0111 //
0112 // constructors and destructor
0113 //
0114 EmbeddingLHEProducer::EmbeddingLHEProducer(const edm::ParameterSet &iConfig) {
0115   // register your products
0116   produces<LHEEventProduct>();
0117   produces<LHERunInfoProduct, edm::Transition::BeginRun>();
0118   produces<math::XYZTLorentzVectorD>("vertexPosition");
0119 
0120   muonsCollection_ = consumes<edm::View<pat::Muon>>(iConfig.getParameter<edm::InputTag>("src"));
0121   vertexCollection_ = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"));
0122   particleToEmbed_ = iConfig.getParameter<int>("particleToEmbed");
0123   mirror_ = iConfig.getParameter<bool>("mirror");
0124   InitialRecoCorrection_ = iConfig.getParameter<bool>("InitialRecoCorrection");
0125   rotate180_ = iConfig.getParameter<bool>("rotate180");
0126   studyFSRmode_ = iConfig.getUntrackedParameter<std::string>("studyFSRmode", "");
0127 
0128   write_lheout = false;
0129   std::string lhe_ouputfile = iConfig.getUntrackedParameter<std::string>("lhe_outputfilename", "");
0130   if (!lhe_ouputfile.empty()) {
0131     write_lheout = true;
0132     file.open(lhe_ouputfile, std::fstream::out | std::fstream::trunc);
0133   }
0134 
0135   // check if particle can be embedded
0136   if (std::find(std::begin(embeddingParticles), std::end(embeddingParticles), particleToEmbed_) ==
0137       std::end(embeddingParticles)) {
0138     throw cms::Exception("Configuration") << "The given particle to embed is not in the list of allowed particles.";
0139   }
0140 
0141   edm::Service<edm::RandomNumberGenerator> rng;
0142   if (!rng.isAvailable()) {
0143     throw cms::Exception("Configuration") << "The EmbeddingLHEProducer requires the RandomNumberGeneratorService\n"
0144                                              "which is not present in the configuration file. \n"
0145                                              "You must add the service\n"
0146                                              "in the configuration file or remove the modules that require it.";
0147   }
0148 }
0149 
0150 EmbeddingLHEProducer::~EmbeddingLHEProducer() {}
0151 
0152 //
0153 // member functions
0154 //
0155 
0156 // ------------ method called to produce the data  ------------
0157 void EmbeddingLHEProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0158   using namespace edm;
0159 
0160   edm::Service<edm::RandomNumberGenerator> rng;
0161   CLHEP::HepRandomEngine *engine = &rng->getEngine(iEvent.streamID());
0162 
0163   edm::Handle<edm::View<pat::Muon>> muonHandle;
0164   iEvent.getByToken(muonsCollection_, muonHandle);
0165   edm::View<pat::Muon> coll_muons = *muonHandle;
0166 
0167   Handle<std::vector<reco::Vertex>> coll_vertices;
0168   iEvent.getByToken(vertexCollection_, coll_vertices);
0169 
0170   TLorentzVector positiveLepton, negativeLepton;
0171   bool mu_plus_found = false;
0172   bool mu_minus_found = false;
0173   lhef::HEPEUP hepeup;
0174   hepeup.IDPRUP = 0;
0175   hepeup.XWGTUP = 1;
0176   hepeup.SCALUP = -1;
0177   hepeup.AQEDUP = -1;
0178   hepeup.AQCDUP = -1;
0179   // Assuming Pt-Order
0180   for (edm::View<pat::Muon>::const_iterator muon = coll_muons.begin(); muon != coll_muons.end(); ++muon) {
0181     if (muon->charge() == 1 && !mu_plus_found) {
0182       assign_4vector(positiveLepton, &(*muon), studyFSRmode_);
0183       mu_plus_found = true;
0184     } else if (muon->charge() == -1 && !mu_minus_found) {
0185       assign_4vector(negativeLepton, &(*muon), studyFSRmode_);
0186       mu_minus_found = true;
0187     } else if (mu_minus_found && mu_plus_found)
0188       break;
0189   }
0190   InitialRecoCorrection(
0191       positiveLepton,
0192       negativeLepton);  // corrects Z mass peak to take into account smearing happening due to first muon reconstruction in the selection step
0193   mirror(positiveLepton, negativeLepton);                    // if no mirror, function does nothing.
0194   rotate180(positiveLepton, negativeLepton);                 // if no rotate180, function does nothing
0195   transform_mumu_to_tautau(positiveLepton, negativeLepton);  // if MuonEmbedding, function does nothing.
0196   fill_lhe_from_mumu(positiveLepton, negativeLepton, hepeup, engine);
0197 
0198   double originalXWGTUP_ = 1.;
0199   std::unique_ptr<LHEEventProduct> product(new LHEEventProduct(hepeup, originalXWGTUP_));
0200 
0201   if (write_lheout)
0202     std::copy(product->begin(), product->end(), std::ostream_iterator<std::string>(file));
0203 
0204   iEvent.put(std::move(product));
0205   // Saving vertex position
0206   std::unique_ptr<math::XYZTLorentzVectorD> vertex_position(
0207       new math::XYZTLorentzVectorD(coll_vertices->at(0).x(), coll_vertices->at(0).y(), coll_vertices->at(0).z(), 0.0));
0208   iEvent.put(std::move(vertex_position), "vertexPosition");
0209 }
0210 
0211 // ------------ method called once each job just before starting event loop  ------------
0212 void EmbeddingLHEProducer::beginJob() {}
0213 
0214 // ------------ method called once each job just after ending the event loop  ------------
0215 void EmbeddingLHEProducer::endJob() {}
0216 
0217 // ------------ method called when starting to processes a run  ------------
0218 
0219 void EmbeddingLHEProducer::beginRunProduce(edm::Run &run, edm::EventSetup const &) {
0220   // fill HEPRUP common block and store in edm::Run
0221   lhef::HEPRUP heprup;
0222 
0223   // set number of processes: 1 for Z to tau tau
0224   heprup.resize(1);
0225 
0226   // master switch for event weight iterpretation (same as in officially produced DYJets LHE files)
0227   heprup.IDWTUP = 3;
0228 
0229   // Information for first process (Z to tau tau), for now only placeholder:
0230   heprup.XSECUP[0] = 1.;
0231   heprup.XERRUP[0] = 0;
0232   heprup.XMAXUP[0] = 1;
0233   heprup.LPRUP[0] = 1;
0234 
0235   std::unique_ptr<LHERunInfoProduct> runInfo(new LHERunInfoProduct(heprup));
0236   runInfo->addHeader(give_slha());
0237 
0238   if (write_lheout)
0239     std::copy(runInfo->begin(), runInfo->end(), std::ostream_iterator<std::string>(file));
0240   run.put(std::move(runInfo));
0241 }
0242 
0243 void EmbeddingLHEProducer::endRunProduce(edm::Run &run, edm::EventSetup const &es) {
0244   if (write_lheout) {
0245     file << LHERunInfoProduct::endOfFile();
0246     file.close();
0247   }
0248 }
0249 
0250 void EmbeddingLHEProducer::fill_lhe_from_mumu(TLorentzVector &positiveLepton,
0251                                               TLorentzVector &negativeLepton,
0252                                               lhef::HEPEUP &outlhe,
0253                                               CLHEP::HepRandomEngine *engine) {
0254   TLorentzVector Z = positiveLepton + negativeLepton;
0255   int leptonPDGID = particleToEmbed_;
0256 
0257   static constexpr double tau_ctau0 = 8.71100e-02;  // mm (for Pythia)
0258   double tau_ctau_p = tau_ctau0 * CLHEP::RandExponential::shoot(engine);
0259   double tau_ctau_n = tau_ctau0 * CLHEP::RandExponential::shoot(engine);
0260 
0261   fill_lhe_with_particle(outlhe, Z, 23, 9.0, 0);
0262   fill_lhe_with_particle(outlhe, positiveLepton, -leptonPDGID, 1.0, tau_ctau_p);
0263   fill_lhe_with_particle(outlhe, negativeLepton, leptonPDGID, -1.0, tau_ctau_n);
0264 
0265   return;
0266 }
0267 
0268 void EmbeddingLHEProducer::fill_lhe_with_particle(
0269     lhef::HEPEUP &outlhe, TLorentzVector &particle, int pdgid, double spin, double ctau) {
0270   // Pay attention to different index conventions:
0271   // 'particleindex' follows usual C++ index conventions starting at 0 for a list.
0272   // 'motherindex' follows the LHE index conventions: 0 is for 'not defined', so the listing starts at 1.
0273   // That means: LHE index 1 == C++ index 0.
0274   int particleindex = outlhe.NUP;
0275   outlhe.resize(outlhe.NUP + 1);
0276 
0277   outlhe.PUP[particleindex][0] = particle.Px();
0278   outlhe.PUP[particleindex][1] = particle.Py();
0279   outlhe.PUP[particleindex][2] = particle.Pz();
0280   outlhe.PUP[particleindex][3] = particle.E();
0281   outlhe.PUP[particleindex][4] = particle.M();
0282   outlhe.IDUP[particleindex] = pdgid;
0283   outlhe.SPINUP[particleindex] = spin;
0284   outlhe.VTIMUP[particleindex] = ctau;
0285 
0286   outlhe.ICOLUP[particleindex].first = 0;
0287   outlhe.ICOLUP[particleindex].second = 0;
0288 
0289   if (std::abs(pdgid) == 23) {
0290     outlhe.MOTHUP[particleindex].first = 0;  // No Mother
0291     outlhe.MOTHUP[particleindex].second = 0;
0292     outlhe.ISTUP[particleindex] = 2;  // status
0293   }
0294 
0295   if (std::find(std::begin(embeddingParticles), std::end(embeddingParticles), std::abs(pdgid)) !=
0296       std::end(embeddingParticles)) {
0297     outlhe.MOTHUP[particleindex].first = 1;   // Mother is the Z (first partile)
0298     outlhe.MOTHUP[particleindex].second = 1;  // Mother is the Z (first partile)
0299 
0300     outlhe.ISTUP[particleindex] = 1;  // status
0301   }
0302 
0303   return;
0304 }
0305 
0306 void EmbeddingLHEProducer::transform_mumu_to_tautau(TLorentzVector &positiveLepton, TLorentzVector &negativeLepton) {
0307   // No corrections applied for muon embedding
0308   double lep_mass;
0309   if (particleToEmbed_ == 11) {
0310     lep_mass = elMass_;
0311   } else if (particleToEmbed_ == 15) {
0312     lep_mass = tauMass_;
0313   } else {
0314     return;
0315   }
0316 
0317   TLorentzVector Z = positiveLepton + negativeLepton;
0318 
0319   TVector3 boost_from_Z_to_LAB = Z.BoostVector();
0320   TVector3 boost_from_LAB_to_Z = -Z.BoostVector();
0321   // Boosting the two leptons to Z restframe, then both are back to back. This means, same 3-momentum squared
0322   positiveLepton.Boost(boost_from_LAB_to_Z);
0323   negativeLepton.Boost(boost_from_LAB_to_Z);
0324 
0325   // Energy of tau = 0.5*Z-mass
0326   double lep_mass_squared = lep_mass * lep_mass;
0327   double lep_energy_squared = 0.25 * Z.M2();
0328   double lep_3momentum_squared = lep_energy_squared - lep_mass_squared;
0329   if (lep_3momentum_squared < 0) {
0330     edm::LogWarning("TauEmbedding") << "3-Momentum squared is negative";
0331     return;
0332   }
0333 
0334   // Computing scale, applying it on the 3-momenta and building new 4 momenta of the taus
0335   double scale = std::sqrt(lep_3momentum_squared / positiveLepton.Vect().Mag2());
0336   positiveLepton.SetPxPyPzE(scale * positiveLepton.Px(),
0337                             scale * positiveLepton.Py(),
0338                             scale * positiveLepton.Pz(),
0339                             std::sqrt(lep_energy_squared));
0340   negativeLepton.SetPxPyPzE(scale * negativeLepton.Px(),
0341                             scale * negativeLepton.Py(),
0342                             scale * negativeLepton.Pz(),
0343                             std::sqrt(lep_energy_squared));
0344 
0345   // Boosting the new taus back to LAB frame
0346   positiveLepton.Boost(boost_from_Z_to_LAB);
0347   negativeLepton.Boost(boost_from_Z_to_LAB);
0348 
0349   return;
0350 }
0351 
0352 void EmbeddingLHEProducer::assign_4vector(TLorentzVector &Lepton, const pat::Muon *muon, std::string FSRmode) {
0353   if ("afterFSR" == FSRmode && muon->genParticle() != nullptr) {
0354     const reco::GenParticle *afterFSRMuon = muon->genParticle();
0355     Lepton.SetPxPyPzE(
0356         afterFSRMuon->p4().px(), afterFSRMuon->p4().py(), afterFSRMuon->p4().pz(), afterFSRMuon->p4().e());
0357   } else if ("beforeFSR" == FSRmode && muon->genParticle() != nullptr) {
0358     const reco::Candidate *beforeFSRMuon = find_original_muon(muon->genParticle());
0359     Lepton.SetPxPyPzE(
0360         beforeFSRMuon->p4().px(), beforeFSRMuon->p4().py(), beforeFSRMuon->p4().pz(), beforeFSRMuon->p4().e());
0361   } else {
0362     Lepton.SetPxPyPzE(muon->p4().px(), muon->p4().py(), muon->p4().pz(), muon->p4().e());
0363   }
0364   return;
0365 }
0366 
0367 const reco::Candidate *EmbeddingLHEProducer::find_original_muon(const reco::Candidate *muon) {
0368   if (muon->mother(0) == nullptr)
0369     return muon;
0370   if (muon->pdgId() == muon->mother(0)->pdgId())
0371     return find_original_muon(muon->mother(0));
0372   else
0373     return muon;
0374 }
0375 
0376 void EmbeddingLHEProducer::rotate180(TLorentzVector &positiveLepton, TLorentzVector &negativeLepton) {
0377   if (!rotate180_)
0378     return;
0379   edm::LogInfo("TauEmbedding") << "Applying 180<C2><B0> rotation";
0380   // By construction, the 3-momenta of mu-, mu+ and Z are in one plane.
0381   // That means, one vector for perpendicular projection can be used for both leptons.
0382   TLorentzVector Z = positiveLepton + negativeLepton;
0383 
0384   edm::LogInfo("TauEmbedding") << "MuMinus before. Pt: " << negativeLepton.Pt() << " Eta: " << negativeLepton.Eta()
0385                                << " Phi: " << negativeLepton.Phi() << " Mass: " << negativeLepton.M();
0386 
0387   TVector3 Z3 = Z.Vect();
0388   TVector3 positiveLepton3 = positiveLepton.Vect();
0389   TVector3 negativeLepton3 = negativeLepton.Vect();
0390 
0391   TVector3 p3_perp = positiveLepton3 - positiveLepton3.Dot(Z3) / Z3.Dot(Z3) * Z3;
0392   p3_perp = p3_perp.Unit();
0393 
0394   positiveLepton3 -= 2 * positiveLepton3.Dot(p3_perp) * p3_perp;
0395   negativeLepton3 -= 2 * negativeLepton3.Dot(p3_perp) * p3_perp;
0396 
0397   positiveLepton.SetVect(positiveLepton3);
0398   negativeLepton.SetVect(negativeLepton3);
0399 
0400   edm::LogInfo("TauEmbedding") << "MuMinus after. Pt: " << negativeLepton.Pt() << " Eta: " << negativeLepton.Eta()
0401                                << " Phi: " << negativeLepton.Phi() << " Mass: " << negativeLepton.M();
0402 
0403   return;
0404 }
0405 
0406 void EmbeddingLHEProducer::InitialRecoCorrection(TLorentzVector &positiveLepton, TLorentzVector &negativeLepton) {
0407   if (!InitialRecoCorrection_)
0408     return;
0409   edm::LogInfo("TauEmbedding") << "Applying initial reconstruction correction";
0410   TLorentzVector Z = positiveLepton + negativeLepton;
0411   edm::LogInfo("TauEmbedding") << "MuMinus before. Pt: " << negativeLepton.Pt() << " Mass: " << negativeLepton.M();
0412   float diLeptonMass = (positiveLepton + negativeLepton).M();
0413   if (diLeptonMass > 60. && diLeptonMass < 122.) {
0414     static constexpr float zmass = 91.1876;
0415     float correction_deviation =
0416         5.;  // to ensure only a correction that drops corresponding to a Gaussian with mean zmass and std. dev. 5 GeV
0417     double EmbeddingCorrection =
0418         1.138;  // value derived by function fitting to fold embedded mass spectrum back to original selection when using mu -> mu embedding
0419     EmbeddingCorrection /=
0420         (EmbeddingCorrection -
0421          (EmbeddingCorrection - 1.) * exp(-pow((diLeptonMass - zmass), 2.) / (2. * pow(correction_deviation, 2.))));
0422     EmbeddingCorrection = ((diLeptonMass + (EmbeddingCorrection - 1.) * zmass) / (diLeptonMass * EmbeddingCorrection));
0423     double correctedpositiveLeptonEnergy =
0424         std::sqrt((pow(muonMass_, 2) / pow(EmbeddingCorrection, 2)) + pow(positiveLepton.Px(), 2) +
0425                   pow(positiveLepton.Py(), 2) + pow(positiveLepton.Pz(), 2));
0426     double correctednegativeLeptonEnergy =
0427         std::sqrt((pow(muonMass_, 2) / pow(EmbeddingCorrection, 2)) + pow(negativeLepton.Px(), 2) +
0428                   pow(negativeLepton.Py(), 2) + pow(negativeLepton.Pz(), 2));
0429     positiveLepton.SetE(correctedpositiveLeptonEnergy);
0430     negativeLepton.SetE(correctednegativeLeptonEnergy);
0431     positiveLepton *= EmbeddingCorrection;
0432     negativeLepton *= EmbeddingCorrection;
0433     edm::LogInfo("TauEmbedding") << "MuMinus after. Pt: " << negativeLepton.Pt() << " Mass: " << negativeLepton.M();
0434   }
0435   return;
0436 }
0437 
0438 void EmbeddingLHEProducer::mirror(TLorentzVector &positiveLepton, TLorentzVector &negativeLepton) {
0439   if (!mirror_)
0440     return;
0441   edm::LogInfo("TauEmbedding") << "Applying mirroring";
0442   TLorentzVector Z = positiveLepton + negativeLepton;
0443 
0444   edm::LogInfo("TauEmbedding") << "MuMinus before. Pt: " << negativeLepton.Pt() << " Eta: " << negativeLepton.Eta()
0445                                << " Phi: " << negativeLepton.Phi() << " Mass: " << negativeLepton.M();
0446 
0447   TVector3 Z3 = Z.Vect();
0448   TVector3 positiveLepton3 = positiveLepton.Vect();
0449   TVector3 negativeLepton3 = negativeLepton.Vect();
0450 
0451   TVector3 beam(0., 0., 1.);
0452   TVector3 perpToZandBeam = Z3.Cross(beam).Unit();
0453 
0454   positiveLepton3 -= 2 * positiveLepton3.Dot(perpToZandBeam) * perpToZandBeam;
0455   negativeLepton3 -= 2 * negativeLepton3.Dot(perpToZandBeam) * perpToZandBeam;
0456 
0457   positiveLepton.SetVect(positiveLepton3);
0458   negativeLepton.SetVect(negativeLepton3);
0459 
0460   edm::LogInfo("TauEmbedding") << "MuMinus after. Pt: " << negativeLepton.Pt() << " Eta: " << negativeLepton.Eta()
0461                                << " Phi: " << negativeLepton.Phi() << " Mass: " << negativeLepton.M();
0462 
0463   return;
0464 }
0465 
0466 LHERunInfoProduct::Header EmbeddingLHEProducer::give_slha() {
0467   LHERunInfoProduct::Header slhah("slha");
0468 
0469   slhah.addLine("######################################################################\n");
0470   slhah.addLine("## PARAM_CARD AUTOMATICALY GENERATED BY MG5 FOLLOWING UFO MODEL   ####\n");
0471   slhah.addLine("######################################################################\n");
0472   slhah.addLine("##                                                                  ##\n");
0473   slhah.addLine("##  Width set on Auto will be computed following the information    ##\n");
0474   slhah.addLine("##        present in the decay.py files of the model.               ##\n");
0475   slhah.addLine("##        See  arXiv:1402.1178 for more details.                    ##\n");
0476   slhah.addLine("##                                                                  ##\n");
0477   slhah.addLine("######################################################################\n");
0478   slhah.addLine("\n");
0479   slhah.addLine("###################################\n");
0480   slhah.addLine("## INFORMATION FOR MASS\n");
0481   slhah.addLine("###################################\n");
0482   slhah.addLine("Block mass \n");
0483   slhah.addLine("    6 1.730000e+02 # MT \n");
0484   slhah.addLine("   15 1.777000e+00 # MTA \n");
0485   slhah.addLine("   23 9.118800e+01 # MZ \n");
0486   slhah.addLine("   25 1.250000e+02 # MH \n");
0487   slhah.addLine("## Dependent parameters, given by model restrictions.\n");
0488   slhah.addLine("## Those values should be edited following the \n");
0489   slhah.addLine("## analytical expression. MG5 ignores those values \n");
0490   slhah.addLine("## but they are important for interfacing the output of MG5\n");
0491   slhah.addLine("## to external program such as Pythia.\n");
0492   slhah.addLine("  1 0.000000 # d : 0.0 \n");
0493   slhah.addLine("  2 0.000000 # u : 0.0 \n");
0494   slhah.addLine("  3 0.000000 # s : 0.0 \n");
0495   slhah.addLine("  4 0.000000 # c : 0.0 \n");
0496   slhah.addLine("  5 0.000000 # b : 0.0 \n");
0497   slhah.addLine("  11 0.000000 # e- : 0.0 \n");
0498   slhah.addLine("  12 0.000000 # ve : 0.0 \n");
0499   slhah.addLine("  13 0.000000 # mu- : 0.0 \n");
0500   slhah.addLine("  14 0.000000 # vm : 0.0 \n");
0501   slhah.addLine("  16 0.000000 # vt : 0.0 \n");
0502   slhah.addLine("  21 0.000000 # g : 0.0 \n");
0503   slhah.addLine("  22 0.000000 # a : 0.0 \n");
0504   slhah.addLine(
0505       "  24 80.419002 # w+ : cmath.sqrt(MZ__exp__2/2. + cmath.sqrt(MZ__exp__4/4. - "
0506       "(aEW*cmath.pi*MZ__exp__2)/(Gf*sqrt__2))) \n");
0507   slhah.addLine("\n");
0508   slhah.addLine("###################################\n");
0509   slhah.addLine("## INFORMATION FOR SMINPUTS\n");
0510   slhah.addLine("###################################\n");
0511   slhah.addLine("Block sminputs \n");
0512   slhah.addLine("    1 1.325070e+02 # aEWM1 \n");
0513   slhah.addLine("    2 1.166390e-05 # Gf \n");
0514   slhah.addLine("    3 1.180000e-01 # aS \n");
0515   slhah.addLine("\n");
0516   slhah.addLine("###################################\n");
0517   slhah.addLine("## INFORMATION FOR WOLFENSTEIN\n");
0518   slhah.addLine("###################################\n");
0519   slhah.addLine("Block wolfenstein \n");
0520   slhah.addLine("    1 2.253000e-01 # lamWS \n");
0521   slhah.addLine("    2 8.080000e-01 # AWS \n");
0522   slhah.addLine("    3 1.320000e-01 # rhoWS \n");
0523   slhah.addLine("    4 3.410000e-01 # etaWS \n");
0524   slhah.addLine("\n");
0525   slhah.addLine("###################################\n");
0526   slhah.addLine("## INFORMATION FOR YUKAWA\n");
0527   slhah.addLine("###################################\n");
0528   slhah.addLine("Block yukawa \n");
0529   slhah.addLine("    6 1.730000e+02 # ymt \n");
0530   slhah.addLine("   15 1.777000e+00 # ymtau \n");
0531   slhah.addLine("\n");
0532   slhah.addLine("###################################\n");
0533   slhah.addLine("## INFORMATION FOR DECAY\n");
0534   slhah.addLine("###################################\n");
0535   slhah.addLine("DECAY   6 1.491500e+00 # WT \n");
0536   slhah.addLine("DECAY  15 2.270000e-12 # WTau \n");
0537   slhah.addLine("DECAY  23 2.441404e+00 # WZ \n");
0538   slhah.addLine("DECAY  24 2.047600e+00 # WW \n");
0539   slhah.addLine("DECAY  25 6.382339e-03 # WH \n");
0540   slhah.addLine("## Dependent parameters, given by model restrictions.\n");
0541   slhah.addLine("## Those values should be edited following the \n");
0542   slhah.addLine("## analytical expression. MG5 ignores those values \n");
0543   slhah.addLine("## but they are important for interfacing the output of MG5\n");
0544   slhah.addLine("## to external program such as Pythia.\n");
0545   slhah.addLine("DECAY  1 0.000000 # d : 0.0 \n");
0546   slhah.addLine("DECAY  2 0.000000 # u : 0.0 \n");
0547   slhah.addLine("DECAY  3 0.000000 # s : 0.0 \n");
0548   slhah.addLine("DECAY  4 0.000000 # c : 0.0 \n");
0549   slhah.addLine("DECAY  5 0.000000 # b : 0.0 \n");
0550   slhah.addLine("DECAY  11 0.000000 # e- : 0.0 \n");
0551   slhah.addLine("DECAY  12 0.000000 # ve : 0.0 \n");
0552   slhah.addLine("DECAY  13 0.000000 # mu- : 0.0 \n");
0553   slhah.addLine("DECAY  14 0.000000 # vm : 0.0 \n");
0554   slhah.addLine("DECAY  16 0.000000 # vt : 0.0 \n");
0555   slhah.addLine("DECAY  21 0.000000 # g : 0.0 \n");
0556   slhah.addLine("DECAY  22 0.000000 # a : 0.0\n");
0557 
0558   return slhah;
0559 }
0560 
0561 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0562 void EmbeddingLHEProducer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0563   // The following says we do not know what parameters are allowed so do no validation
0564   //  Please change this to state exactly what you do use, even if it is no parameters
0565   edm::ParameterSetDescription desc;
0566   desc.setUnknown();
0567   descriptions.addDefault(desc);
0568 }
0569 
0570 // define this as a plug-in
0571 DEFINE_FWK_MODULE(EmbeddingLHEProducer);