Pythia8Interface/plugins/SuepDecay.cc

0001 #include "SuepDecay.h"
0002
0003 SuepDecay::SuepDecay(const edm::ParameterSet& iConfig)
0004     : idMediator_(iConfig.getParameter<int>("idMediator")),
0005       idDark_(iConfig.getParameter<int>("idDark")),
0006       temperature_(iConfig.getParameter<double>("temperature")) {}
0007
0008 bool SuepDecay::initAfterBeams() {
0009   mDark_ = particleDataPtr->m0(idDark_);
0010   bool medDecay = particleDataPtr->mayDecay(idMediator_);
0011   if (!medDecay) {
0012     loggerPtr->errorMsg("SuepDecay::initAfterBeams", "mediator decay should be enabled");
0013     return false;
0014   }
0015
0016   //construct the shower helper
0017   suep_shower_ = std::make_unique<SuepShower>(mDark_, temperature_, rndmPtr);
0018
0019   return true;
0020 }
0021
0022 //based on https://gitlab.com/simonknapen/suep_generator/-/blob/master/suep_main.cc:AttachSuepShower
0023 bool SuepDecay::doVetoProcessLevel(Pythia8::Event& event) {
0024   Pythia8::Vec4 pMediator, pDark;
0025
0026   // Find the mediator in the event
0027   for (int i = 0; i < event.size(); ++i) {
0028     //mediator w/ distinct daughters = last copy (decayed)
0029     if (event[i].id() == idMediator_ && event[i].daughter1() != event[i].daughter2() && event[i].daughter1() > 0 &&
0030         event[i].daughter2() > 0) {
0031       pMediator = event[i].p();
0032       // undo mediator decay
0033       event[i].undoDecay();
0034
0035       // Generate the shower, output are 4 vectors in the rest frame of the shower, adding energy here avoids issues if scalar is off-shell
0036       std::vector<Pythia8::Vec4> suep_shower_fourmomenta = suep_shower_->generateShower(pMediator.mCalc());
0037       // Loop over hidden sector mesons and append to the event
0038       int firstDaughter = event.size();
0039       for (auto& pDark : suep_shower_fourmomenta) {
0040         // Boost to the lab frame, i.e. apply the mediator boost
0041         pDark.bst(pMediator);
0042         // Append particle to the event w/ hidden meson pdg code. Magic number 91 means it is produced as a normal decay product
0043         event.append(idDark_, 91, i, 0, 0, 0, 0, 0, pDark.px(), pDark.py(), pDark.pz(), pDark.e(), mDark_);
0044       }
0045
0046       // Change the status code of the mediator to reflect that it has decayed.
0047       event[i].statusNeg();
0048
0049       //set daughters of the mediator: daughter1 < daughter2 > 0 -> the particle has a range of decay products from daughter1 to daughter2
0050       event[i].daughters(firstDaughter, event.size() - 1);
0051       break;
0052     }
0053   }
0054
0055   //allow event to continue
0056   return false;
0057 }