CMS 3D CMS Logo

Py8InterfaceBase.cc
Go to the documentation of this file.
1 #include <memory>
2 #include <string>
3 
6 
9 #include <filesystem>
10 
11 // EvtGen plugin
12 //
13 //#include "Pythia8Plugins/EvtGen.h"
14 
15 using namespace Pythia8;
16 
17 namespace gen {
18 
19  Py8InterfaceBase::Py8InterfaceBase(edm::ParameterSet const& ps)
20  : BaseHadronizer(ps), useEvtGen(false), evtgenDecays(nullptr) {
21  fParameters = ps;
22 
23  pythiaPylistVerbosity = ps.getUntrackedParameter<int>("pythiaPylistVerbosity", 0);
24  pythiaHepMCVerbosity = ps.getUntrackedParameter<bool>("pythiaHepMCVerbosity", false);
25  pythiaHepMCVerbosityParticles = ps.getUntrackedParameter<bool>("pythiaHepMCVerbosityParticles", false);
26  maxEventsToPrint = ps.getUntrackedParameter<int>("maxEventsToPrint", 0);
27  p8RndmEngine_ = std::make_shared<P8RndmEngine>();
28 
30  ascii_io = new HepMC::IO_AsciiParticles("cout", std::ios::out);
31 
32  if (ps.exists("useEvtGenPlugin")) {
33  useEvtGen = true;
34  string evtgenpath(std::getenv("EVTGENDATA"));
35  evtgenDecFile = evtgenpath + string("/DECAY_2010.DEC");
36  evtgenPdlFile = evtgenpath + string("/evt.pdl");
37 
38  if (ps.exists("evtgenDecFile")) {
39  edm::FileInPath decay_table(ps.getParameter<std::string>("evtgenDecFile"));
40  evtgenDecFile = decay_table.fullPath();
41  }
42 
43  if (ps.exists("evtgenPdlFile")) {
44  edm::FileInPath pdt(ps.getParameter<std::string>("evtgenPdlFile"));
45  evtgenPdlFile = pdt.fullPath();
46  }
47 
48  if (ps.exists("evtgenUserFile")) {
49  std::vector<std::string> user_decays = ps.getParameter<std::vector<std::string> >("evtgenUserFile");
50  for (unsigned int i = 0; i < user_decays.size(); i++) {
51  edm::FileInPath user_decay(user_decays.at(i));
52  evtgenUserFiles.push_back(user_decay.fullPath());
53  }
54  //evtgenUserFiles = ps.getParameter< std::vector<std::string> >("evtgenUserFile");
55  }
56 
57  if (ps.exists("evtgenUserFileEmbedded")) {
58  std::vector<std::string> user_decay_lines =
59  ps.getParameter<std::vector<std::string> >("evtgenUserFileEmbedded");
60  char tempslhaname[] = "pythia8evtgenXXXXXX";
61  int fd = mkstemp(tempslhaname);
62 
63  for (unsigned int i = 0; i < user_decay_lines.size(); i++) {
64  user_decay_lines.at(i) += "\n";
65  write(fd, user_decay_lines.at(i).c_str(), user_decay_lines.at(i).size());
66  }
67  close(fd);
68  evtgenUserFiles.push_back(std::string(tempslhaname));
69  }
70  }
71  }
72 
74  //Pythia 8's default value for first argument to constructor
75  const string xmlDir = "../share/Pythia8/xmldoc";
76  bool printBanner = true;
77  if (fParameters.exists("printBanner")) {
78  printBanner = fParameters.getUntrackedParameter<bool>("printBanner");
79  }
80  if (!fMasterGen.get())
81  fMasterGen = std::make_unique<Pythia>(xmlDir, printBanner);
82  fDecayer = std::make_unique<Pythia>(xmlDir, printBanner);
83 
84  //add settings for resonance decay filter
85  fMasterGen->settings.addFlag("BiasedTauDecayer:filter", false);
86  fMasterGen->settings.addFlag("BiasedTauDecayer:eDecays", true);
87  fMasterGen->settings.addFlag("BiasedTauDecayer:muDecays", true);
88 
89  //add settings for resonance decay filter
90  fMasterGen->settings.addFlag("ResonanceDecayFilter:filter", false);
91  fMasterGen->settings.addFlag("ResonanceDecayFilter:exclusive", false);
92  fMasterGen->settings.addFlag("ResonanceDecayFilter:eMuAsEquivalent", false);
93  fMasterGen->settings.addFlag("ResonanceDecayFilter:eMuTauAsEquivalent", false);
94  fMasterGen->settings.addFlag("ResonanceDecayFilter:allNuAsEquivalent", false);
95  fMasterGen->settings.addFlag("ResonanceDecayFilter:udscAsEquivalent", false);
96  fMasterGen->settings.addFlag("ResonanceDecayFilter:udscbAsEquivalent", false);
97  fMasterGen->settings.addFlag("ResonanceDecayFilter:wzAsEquivalent", false);
98  fMasterGen->settings.addMVec("ResonanceDecayFilter:mothers", std::vector<int>(), false, false, 0, 0);
99  fMasterGen->settings.addMVec("ResonanceDecayFilter:daughters", std::vector<int>(), false, false, 0, 0);
100 
101  //add settings for PT filter
102  fMasterGen->settings.addFlag("PTFilter:filter", false);
103  fMasterGen->settings.addMode("PTFilter:quarkToFilter", 5, true, true, 3, 6);
104  fMasterGen->settings.addParm("PTFilter:scaleToFilter", 0.4, true, true, 0.0, 10.);
105  fMasterGen->settings.addParm("PTFilter:quarkRapidity", 10.0, true, true, 0.0, 10.);
106  fMasterGen->settings.addParm("PTFilter:quarkPt", -.1, true, true, -.1, 100.);
107 
108  //add settings for RecoilToTop tool
109  fMasterGen->settings.addFlag("TopRecoilHook:doTopRecoilIn", false);
110  fMasterGen->settings.addFlag("TopRecoilHook:useOldDipoleIn", false);
111  fMasterGen->settings.addFlag("TopRecoilHook:doListIn", false);
112 
113  //add settings for powheg resonance scale calculation
114  fMasterGen->settings.addFlag("POWHEGres:calcScales", false);
115  fMasterGen->settings.addFlag("POWHEG:bb4l", false);
116  fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:veto", false);
117  fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:vetoQED", false);
118  fMasterGen->settings.addFlag("POWHEG:bb4l:DEBUG", false);
119  fMasterGen->settings.addFlag("POWHEG:bb4l:ScaleResonance:veto", false);
120  fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:vetoDipoleFrame", false);
121  fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:pTpythiaVeto", false);
122  fMasterGen->settings.addParm("POWHEG:bb4l:pTminVeto", 10.0, true, true, 0.0, 10.);
123  fMasterGen->settings.addFlag("POWHEG:bb4l:vetoAllRadtypes", false);
124 
125  fMasterGen->setRndmEnginePtr(p8RndmEngine_);
126  fDecayer->setRndmEnginePtr(p8RndmEngine_);
127 
128  fMasterGen->readString("Next:numberShowEvent = 0");
129  fDecayer->readString("Next:numberShowEvent = 0");
130 
131  edm::ParameterSet currentParameters;
132  if (randomIndex() >= 0) {
133  std::vector<edm::ParameterSet> randomizedParameters =
134  fParameters.getParameter<std::vector<edm::ParameterSet> >("RandomizedParameters");
135  currentParameters = randomizedParameters[randomIndex()];
136  } else {
137  currentParameters = fParameters;
138  }
139 
140  ParameterCollector pCollector = currentParameters.getParameter<edm::ParameterSet>("PythiaParameters");
141 
142  for (ParameterCollector::const_iterator line = pCollector.begin(); line != pCollector.end(); ++line) {
143  if (line->find("Random:") != std::string::npos)
144  throw cms::Exception("PythiaError") << "Attempted to set random number "
145  "using Pythia commands. Please use "
146  "the RandomNumberGeneratorService."
147  << std::endl;
148 
149  if (!fMasterGen->readString(*line))
150  throw cms::Exception("PythiaError") << "Pythia 8 did not accept \"" << *line << "\"." << std::endl;
151 
152  if (line->find("ParticleDecays:") != std::string::npos) {
153  if (!fDecayer->readString(*line))
154  throw cms::Exception("PythiaError") << "Pythia 8 Decayer did not accept \"" << *line << "\"." << std::endl;
155  }
156  }
157 
158  slhafile_.clear();
159 
160  if (currentParameters.exists("SLHAFileForPythia8")) {
161  std::string slhafilenameshort = currentParameters.getParameter<std::string>("SLHAFileForPythia8");
162  edm::FileInPath f1(slhafilenameshort);
163 
164  fMasterGen->settings.mode("SLHA:readFrom", 2);
165  fMasterGen->settings.word("SLHA:file", f1.fullPath());
166  } else if (currentParameters.exists("SLHATableForPythia8")) {
167  std::string slhatable = currentParameters.getParameter<std::string>("SLHATableForPythia8");
168 
169  makeTmpSLHA(slhatable);
170  } else if (currentParameters.exists("SLHATreeForPythia8")) {
171  auto slhaReaderParams = currentParameters.getParameter<edm::ParameterSet>("SLHATreeForPythia8");
172  std::unique_ptr<SLHAReaderBase> reader =
173  SLHAReaderFactory::get()->create(slhaReaderParams.getParameter<std::string>("name"), slhaReaderParams);
174 
175  makeTmpSLHA(reader->getSLHA(currentParameters.getParameter<std::string>("ConfigDescription")));
176  }
177 
178  return true;
179  }
180 
182  char tempslhaname[] = "pythia8SLHAtableXXXXXX";
183  int fd = mkstemp(tempslhaname);
184  write(fd, slhatable.c_str(), slhatable.size());
185  close(fd);
186 
187  slhafile_ = tempslhaname;
188 
189  fMasterGen->settings.mode("SLHA:readFrom", 2);
190  fMasterGen->settings.word("SLHA:file", slhafile_);
191  }
192 
193  bool Py8InterfaceBase::declareStableParticles(const std::vector<int>& pdgIds) {
194  for (size_t i = 0; i < pdgIds.size(); i++) {
195  // FIXME: need to double check if PID's are the same in Py6 & Py8,
196  // because the HepPDT translation tool is actually for **Py6**
197  //
198  // well, actually it looks like Py8 operates in PDT id's rather than Py6's
199  //
200  // int PyID = HepPID::translatePDTtoPythia( pdgIds[i] );
201  int PyID = pdgIds[i];
202  std::ostringstream pyCard;
203  pyCard << PyID << ":mayDecay=false";
204 
205  if (fMasterGen->particleData.isParticle(PyID)) {
206  fMasterGen->readString(pyCard.str());
207  } else {
208  edm::LogWarning("DataNotUnderstood") << "Pythia8 does not "
209  << "recognize particle id = " << PyID << std::endl;
210  }
211  // alternative:
212  // set the 2nd input argument warn=false
213  // - this way Py8 will NOT print warnings about unknown particle code(s)
214  // fMasterPtr->readString( pyCard.str(), false )
215  }
216 
217  return true;
218  }
219 
220  bool Py8InterfaceBase::declareSpecialSettings(const std::vector<std::string>& settings) {
221  for (unsigned int iss = 0; iss < settings.size(); iss++) {
222  if (settings[iss].find("QED-brem-off") != std::string::npos) {
223  fMasterGen->readString("TimeShower:QEDshowerByL=off");
224  } else {
225  size_t fnd1 = settings[iss].find("Pythia8:");
226  if (fnd1 != std::string::npos) {
227  std::string value = settings[iss].substr(fnd1 + 8);
228  fDecayer->readString(value);
229  }
230  }
231  }
232  return true;
233  }
234 
236  fMasterGen->stat();
237  return;
238  }
239 
240 } // namespace gen
bool declareSpecialSettings(const std::vector< std::string > &)
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
std::shared_ptr< P8RndmEngine > p8RndmEngine_
std::unique_ptr< Pythia8::Pythia > fDecayer
void makeTmpSLHA(const std::string &)
std::string fullPath() const
Definition: FileInPath.cc:161
HepMC::IO_AsciiParticles * ascii_io
const_iterator begin() const
bool exists(std::string const &parameterName) const
checks if a parameter exists
const_iterator end() const
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
std::vector< std::string > evtgenUserFiles
T getUntrackedParameter(std::string const &, T const &) const
int randomIndex() const
unsigned int pythiaPylistVerbosity
Definition: value.py:1
bool declareStableParticles(const std::vector< int > &)
unsigned int maxEventsToPrint
std::unique_ptr< Pythia8::Pythia > fMasterGen
#define get
Log< level::Warning, false > LogWarning
edm::ParameterSet fParameters
fd
Definition: ztee.py:136