ATLAS Offline Software
Pythia8_i.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 #include "Pythia8_i/Pythia8_i.h"
8 
10 #include <boost/algorithm/string.hpp>
11 #include <boost/lexical_cast.hpp>
12 
13 // calls to fortran routines
16 
17 #include <sstream>
18 // For limits
19 #include <limits>
20 
21 
22 namespace {
23  // Name of random number stream
24  std::string s_pythia_stream{"PYTHIA8_INIT"};
25 }
26 
27 
28 // fix Pythia8 shower weights change in conventions
29 #define PYTHIA8_NWEIGHTS nWeights
30 #define PYTHIA8_WEIGHT weight
31 #define PYTHIA8_WLABEL weightLabel
32 #define PYTHIA8_CONVERSION 1.0
33 
34 #ifdef PYTHIA_VERSION_INTEGER
35  #undef PYTHIA8_NWEIGHTS
36  #undef PYTHIA8_WEIGHT
37  #undef PYTHIA8_WLABEL
38  #if PYTHIA_VERSION_INTEGER > 8303
39  #define PYTHIA8_NWEIGHTS nWeightGroups
40  #else
41  #define PYTHIA8_NWEIGHTS nVariationGroups
42  #endif
43  #define PYTHIA8_WEIGHT getGroupWeight
44  #define PYTHIA8_WLABEL getGroupName
45 #endif
46 
51 namespace {
52 
53 std::string py8version()
54 {
55  static const std::string incdir (PY8INCLUDE_DIR);
56  std::string::size_type pos = incdir.find ("/pythia8/");
57  if (pos == std::string::npos) return "";
58  pos += 9;
59  std::string::size_type pos2 = incdir.find ("/", pos);
60  if (pos2 == std::string::npos) pos2 = incdir.size();
61  return incdir.substr (pos, pos2-pos);
62 }
63 
64 
65 }
66 
68 Pythia8_i::Pythia8_i(const std::string &name, ISvcLocator *pSvcLocator)
69 : GenModule(name, pSvcLocator)
70 {
71  m_particleIDs["PROTON"] = PROTON;
72  m_particleIDs["ANTIPROTON"] = ANTIPROTON;
73  m_particleIDs["ELECTRON"] = ELECTRON;
74  m_particleIDs["POSITRON"] = POSITRON;
75  m_particleIDs["NEUTRON"] = NEUTRON;
76  m_particleIDs["ANTINEUTRON"] = ANTINEUTRON;
77  m_particleIDs["MUON"] = MUON;
78  m_particleIDs["ANTIMUON"] = ANTIMUON;
79  m_particleIDs["LEAD"] = LEAD;
80  m_particleIDs["OXYGEN"] = OXYGEN;
81  m_particleIDs["HELIUM"] = HELIUM;
82 
83 }
84 
87 
88  #ifndef PYTHIA8_3SERIES
89 // if(m_userHookPtr != 0) delete m_userHookPtr;
90 
91 // for(UserHooksPtrType ptr: m_userHooksPtrs){
92 // delete ptr;
93 // }
94  #endif
95 }
96 
99 
100  ATH_MSG_DEBUG("Pythia8_i from genInitialize()");
101 
102  ATH_MSG_INFO("XML Path is " + xmlpath());
103  m_pythia = std::make_unique<Pythia8::Pythia> (xmlpath());
104 #ifdef HEPMC3
105  m_runinfo = std::make_shared<HepMC3::GenRunInfo>();
107  struct HepMC3::GenRunInfo::ToolInfo generator={std::string("Pythia8"),py8version(),std::string("Used generator")};
108  m_runinfo->tools().push_back(generator);
109 #endif
110 
111  bool canInit = true;
112 
113  m_version = m_pythia->settings.parm("Pythia:versionNumber");
114 
115  s_pythia_stream = "PYTHIA8_INIT";
116 
117  // By default add "nominal" to the list of shower weight names
119  m_showerWeightNames.insert(m_showerWeightNames.begin(), "nominal");
120 
121  // We do explicitly set tune 4C, since it is the starting point for many other tunes
122  // Tune 4C for pp collisions
123  m_pythia->readString("Tune:pp = 5");
124 
125  // also use CTEQ6L1 from LHAPDF 6 by default
126  // can be over-written using JO
127  m_pythia->readString("PDF:pSet= LHAPDF6:cteq6l1");
128 
129  // switch off verbose event print out
130  m_pythia->readString("Next:numberShowEvent = 0");
131 
132  // Add flag to switch off from JO the Pythia8ToHepMC::print_inconsistency internal variable
133  m_pythia->settings.addFlag("AthenaPythia8ToHepMC:print_inconsistency",true);
134 
135  // Add UserHooks first because these potentially add new settings that must exist prior to parsing commands
136 
137  for(const auto &hook: m_userHooks){
138  ATH_MSG_INFO("Adding user hook " + hook + ".");
139  bool canSetHook = true;
140  if (hook == "SuppressSmallPT") {
142  canSetHook=m_pythia->setUserHooksPtr(PYTHIA8_PTRWRAP(m_SuppressSmallPT));
143 }
144  else {
146  canSetHook = m_pythia->addUserHooksPtr(m_userHooksPtrs.back());
147 }
148  if(!canSetHook){
149  ATH_MSG_ERROR("Unable to set requested user hook.");
150  ATH_MSG_ERROR("Pythia 8 initialisation will FAIL!");
151  canInit=false;
152  }
153  }
154 
155  for(const std::pair<const std::string, double> &param : Pythia8_UserHooks::UserHooksFactory::userSettings<double>()){
156  m_pythia->settings.addParm(param.first, param.second, false, false, 0., 0.);
157  }
158 
159  for(const std::pair<const std::string, int> &param : Pythia8_UserHooks::UserHooksFactory::userSettings<int>()){
160  m_pythia->settings.addMode(param.first, param.second, false, false, 0., 0.);
161  }
162 
163  for(const std::pair<const std::string, bool> &param : Pythia8_UserHooks::UserHooksFactory::userSettings<bool>()){
164  m_pythia->settings.addFlag(param.first, param.second);
165  }
166 
167  for(const std::pair<const std::string, std::string> &param : Pythia8_UserHooks::UserHooksFactory::userSettings<std::string>()){
168  m_pythia->settings.addWord(param.first, param.second);
169  }
170 
171  if( ! m_athenaTool.empty() ){
172  if(m_athenaTool.retrieve().isFailure()){
173  ATH_MSG_ERROR("Unable to retrieve Athena Tool for custom Pythia processing");
174  return StatusCode::FAILURE;
175  }
176  else {
177  StatusCode status = m_athenaTool->InitializePythiaInfo(*m_pythia);
178  if(status != StatusCode::SUCCESS) return status;
179  }
180  }
181 
182  // Now apply the settings from the JO
183  for(const std::string &cmd : m_commands){
184 
185  if(cmd.compare("")==0) continue;
186  else if (cmd.find("Beams:id") != std::string::npos ) {
187  ATH_MSG_ERROR("With command '" << cmd << "' you are trying to set a beam different from p: please use the Beam1/Beam2 properties instead:");
188  ATH_MSG_ERROR(" example: genSeq.Pythia8.Beam1 = 'LEAD'");
189  ATH_MSG_ERROR(" genSeq.Pythia8.Beam2 = 'ANTINEUTRON'");
190  return StatusCode::FAILURE;
191  }
192  else if (cmd.find("Beams:frameType") != std::string::npos ) {
194  ATH_MSG_WARNING(" Found an explicit 'Beams:frameType' command: this will override transform beams momenta/energy parameters, regardless of its requested value. ");
195  }
196 
197  bool read = m_pythia->readString(cmd);
198 
199  if(!read){
200  ATH_MSG_ERROR("Pythia could not understand the command '"<< cmd<<"'");
201  return StatusCode::FAILURE;
202  }
203  }
204 
207 
208  ATH_MSG_INFO("Beam1 = "<<beam1);
209  ATH_MSG_INFO("Beam2 = "<<beam2);
210 
211  if(beam1 == 0 || beam2 == 0){
212  ATH_MSG_ERROR("Invalid beam particle!");
213  return StatusCode::FAILURE;
214  }
215 
216 
217 
218  if(m_useRndmGenSvc){
219 
220  ATH_MSG_INFO(" !!!!!!!!!!!! WARNING ON PYTHIA RANDOM NUMBERS !!!!!!!!!! ");
221  ATH_MSG_INFO(" THE ATHENA SERVICE AthRNGenSvc IS USED.");
222  ATH_MSG_INFO(" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ");
223 
224  m_atlasRndmEngine = std::make_shared<customRndm>();
225  CLHEP::HepRandomEngine* rndmEngine = getRandomEngineDuringInitialize(s_pythia_stream, m_randomSeed, m_dsid); // NOT THREAD-SAFE
226  m_atlasRndmEngine->init(rndmEngine);
227 #if PYTHIA_VERSION_INTEGER >= 8310
228  m_pythia->setRndmEnginePtr(m_atlasRndmEngine);
229 #else
230  m_pythia->setRndmEnginePtr(m_atlasRndmEngine.get());
231 #endif
232  s_pythia_stream = "PYTHIA8";
233  }else{
234  ATH_MSG_INFO(" !!!!!!!!!!!! WARNING ON PYTHIA RANDOM NUMBERS !!!!!!!!!! ");
235  ATH_MSG_INFO(" THE STANDARD PYTHIA8 RANDOM NUMBER SERVICE IS USED.");
236  ATH_MSG_INFO(" THE ATHENA SERVICE AthRNGSvc IS ***NOT*** USED.");
237  ATH_MSG_INFO(" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ");
238  }
239 
240  if(m_userProcess != ""){
242  }
243 
244 
245 
246  if(m_userResonances.value() != ""){
247 
248  std::vector<std::string> resonanceArgs;
249 
250  boost::split(resonanceArgs, m_userResonances.value(), boost::is_any_of(":"));
251  if(resonanceArgs.size() != 2){
252  ATH_MSG_ERROR("Cannot Understand UserResonance job option!");
253  ATH_MSG_ERROR("You should specify it as a 'name:id1,id2,id3...'");
254  ATH_MSG_ERROR("Where name is the name of your UserResonance, and id1,id2,id3 are a comma separated list of PDG IDs to which it is applied");
255  canInit = false;
256  }
257  std::vector<std::string> resonanceIds;
258  boost::split(resonanceIds, resonanceArgs.back(), boost::is_any_of(","));
259  if(resonanceIds.size()==0){
260  ATH_MSG_ERROR("You did not specifiy any PDG ids to which your user resonance width should be applied!");
261  ATH_MSG_ERROR("You should specify a list as 'name:id1,id2,id3...'");
262  ATH_MSG_ERROR("Where name is the name of your UserResonance, and id1,id2,id3 are a comma separated list of PDG IDs to which it is applied");
263  canInit=false;
264  }
265 
266  for(std::vector<std::string>::const_iterator sId = resonanceIds.begin();
267  sId != resonanceIds.end(); ++sId){
268  int idResIn = boost::lexical_cast<int>(*sId);
269  m_userResonancePtrs.push_back(Pythia8_UserResonance::UserResonanceFactory::create(resonanceArgs.front(), idResIn));
270  }
271 
272  for(std::shared_ptr<Pythia8::ResonanceWidths>& resonance : m_userResonancePtrs) {
273 #if PYTHIA_VERSION_INTEGER >= 8310
274  m_pythia->setResonancePtr(resonance);
275 #else
276  m_pythia->setResonancePtr(resonance.get());
277 #endif
278  }
279 
280  }
281 
282  if(m_particleDataFile != "") {
283  if(!m_pythia->particleData.reInit(m_particleDataFile, true)){
284  ATH_MSG_ERROR("Unable to read requested particle data table: " + m_particleDataFile + " !!");
285  ATH_MSG_ERROR("Pythia 8 initialisation will FAIL!");
286  canInit = false;
287  }
288  }
289 
290  if(m_lheFile != ""){
291  if(m_procPtr){
292  ATH_MSG_ERROR("Both LHE file and user process have been specified");
293  ATH_MSG_ERROR("LHE input does not make sense with a user process!");
294  canInit = false;
295  }
296 
297  canInit = canInit && m_pythia->readString("Beams:frameType = 4");
298  canInit = canInit && m_pythia->readString("Beams:LHEF = " + m_lheFile);
299  if(!canInit){
300  ATH_MSG_ERROR("Unable to read requested LHE file: " + m_lheFile + " !");
301  ATH_MSG_ERROR("Pythia 8 initialisation will FAIL!");
302  }
303  }else{
305  canInit = canInit && m_pythia->readString("Beams:frameType = 1");
306  canInit = canInit && m_pythia->readString("Beams:eCM = " + std::to_string(m_collisionEnergy));
307  }
308  canInit = canInit && m_pythia->readString("Beams:idA = " + std::to_string(beam1));
309  canInit = canInit && m_pythia->readString("Beams:idB = " + std::to_string(beam2));
310  }
311 
312  if(m_procPtr){
313 #if PYTHIA_VERSION_INTEGER >= 8310
314  if(!m_pythia->setSigmaPtr(m_procPtr))
315 #else
316  if(!m_pythia->setSigmaPtr(m_procPtr.get()))
317 #endif
318  {
319  ATH_MSG_ERROR("Unable to set requested user process: " + m_userProcess + " !!");
320  ATH_MSG_ERROR("Pythia 8 initialisation will FAIL!");
321  canInit = false;
322  }
323  }
324 
325  StatusCode returnCode = StatusCode::SUCCESS;
326 
327  m_pythia->particleData.listXML(m_outputParticleDataFile.value().substr(0,m_outputParticleDataFile.value().find("xml"))+"orig.xml");
328  m_pythia->settings.writeFile("Settings_before.log",true);
329 
330  if(canInit){
331  canInit = m_pythia->init();
332  }
333 
334  if(!canInit){
335  returnCode = StatusCode::FAILURE;
336  ATH_MSG_ERROR(" *** Unable to initialise Pythia !! ***");
337  }
338 
339  m_pythia->particleData.listXML(m_outputParticleDataFile);
340  m_pythia->settings.writeFile("Settings_after.log",true);
341 
342  //counter for event failures;
343  m_failureCount = 0;
344 
346 
347  // Set set_print_inconsistency to Athena corresponding flag (allowing to change it from JO)
348  m_pythiaToHepMC.set_print_inconsistency( m_pythia->settings.flag("AthenaPythia8ToHepMC:print_inconsistency") );
349 
350  m_pythiaToHepMC.set_store_pdf(true);
351  m_pythiaToHepMC.set_store_weights(false);
352 
353  return returnCode;
354 }
355 
358 
359  ATH_MSG_DEBUG(">>> Pythia8_i from callGenerator");
361  //Re-seed the random number stream
362  long seeds[7];
363  const EventContext& ctx = Gaudi::Hive::currentContext();
364  ATHRNG::calculateSeedsMC21(seeds, s_pythia_stream, ctx.eventID().event_number(), m_dsid, m_randomSeed);
365  m_atlasRndmEngine->getEngine()->setSeeds(seeds, 0); // NOT THREAD-SAFE
366  m_seeds.clear();
367  m_seeds.push_back(seeds[0]);
368  m_seeds.push_back(seeds[1]);
369  }
370 
371  bool status = m_pythia->next();
372 
373  StatusCode returnCode = StatusCode::SUCCESS;
374 
375  if(!status){
376  ++m_failureCount;
378  ATH_MSG_ERROR("Exceeded the max number of consecutive event failures.");
379  returnCode = StatusCode::FAILURE;
380  }else{
381  ATH_MSG_INFO("Event generation failed - re-trying.");
382  returnCode = this->callGenerator();
383  }
384  }
385 
386  ATH_MSG_DEBUG("Now checking with Tool.empty() second time");
387  if( ! m_athenaTool.empty() ){
388  StatusCode stat = m_athenaTool->ModifyPythiaEvent(*m_pythia);
389  if(stat != StatusCode::SUCCESS) returnCode = stat;
390  }
391 
392  m_failureCount = 0;
393 
394  m_nAccepted += 1.;
395  // some CKKWL merged events have zero weight (or unfilled event).
396  // start again with such events
397 
398  double eventWeight = m_pythia->info.mergingWeight()*m_pythia->info.weight();
399 
400 
401  if(returnCode != StatusCode::FAILURE &&
402  // Here this is a double, but it may be converted into float at some point downstream
403  (std::abs(eventWeight) < std::numeric_limits<float>::min() ||
404  m_pythia->event.size() < 2)){
405 
406  returnCode = this->callGenerator();
407  } else if ( std::abs(eventWeight) < std::numeric_limits<float>::min() &&
408  std::abs(eventWeight) > std::numeric_limits<double>::min() ){
409  ATH_MSG_WARNING("Found event weight " << eventWeight << " between the float and double precision limits. Rejecting event.");
410  } else {
411  m_nMerged += eventWeight;
413 
414  // For FxFx cross section:
415  m_sigmaTotal+=m_pythia->info.weight();
416  }
417 
418  return returnCode;
419 }
420 
423 
424  ATH_MSG_DEBUG(">>> Pythia8_i from fillEvt");
425 
426  evt->set_event_number(m_internal_event_number);
427 
428  // if using "getGroupWeight" and | lhastrategy | = 4, then need to convert mb to pb ( done otherwise when calling info.weight(), [...] )
429  if( m_internal_event_number == 1 && std::abs(m_pythia->info.lhaStrategy()) == 4 ) {
431  ATH_MSG_DEBUG(" LHA strategy needs a conversion to fix Pythia8 shower weights bug(s) equal to " << m_conversion);
432  }
433 
434 
435  if(m_pythia->event.size() < 2){
436  ATH_MSG_ERROR("Something wrong with this event - it contains fewer than 2 particles!");
437  ATH_MSG_ERROR("internal event number is "<<m_internal_event_number);
438  return StatusCode::FAILURE;
439  }
440 
442 
443  // In case we're asked, save the LHE event as an attribute. This only works in HEPMC3.
444  if(m_saveLHE){
445 #ifdef HEPMC3
446  auto evtlhe = std::make_shared<HepMC::GenEvent>();
447  m_pythiaToHepMC.fill_next_event(m_pythia->process, evtlhe.get(), m_internal_event_number, &m_pythia->info, &m_pythia->settings);
448  auto extra = std::make_shared<HepMC::ShortEventAttribute>(evtlhe.get());
449  evt->add_attribute("LHERecord", extra);
450 #else
451  ATH_MSG_WARNING("LHE record saving requested, but not implemented with HEPMC2 functionality.");
452 #endif
453  }
454 
455  // in debug mode you can check whether the pdf information is stored
456  if(evt->pdf_info()){
457 #ifdef HEPMC3
458  ATH_MSG_DEBUG("PDFinfo id1:" << evt->pdf_info()->parton_id[0]);
459  ATH_MSG_DEBUG("PDFinfo id2:" << evt->pdf_info()->parton_id[1]);
460  ATH_MSG_DEBUG("PDFinfo x1:" << evt->pdf_info()->x[0]);
461  ATH_MSG_DEBUG("PDFinfo x2:" << evt->pdf_info()->x[1]);
462  ATH_MSG_DEBUG("PDFinfo scalePDF:" << evt->pdf_info()->scale);
463  ATH_MSG_DEBUG("PDFinfo pdf1:" << evt->pdf_info()->pdf_id[0]);
464  ATH_MSG_DEBUG("PDFinfo pdf2:" << evt->pdf_info()->pdf_id[1]);
465 #else
466  ATH_MSG_DEBUG("PDFinfo id1:" << evt->pdf_info()->id1());
467  ATH_MSG_DEBUG("PDFinfo id2:" << evt->pdf_info()->id2());
468  ATH_MSG_DEBUG("PDFinfo x1:" << evt->pdf_info()->x1());
469  ATH_MSG_DEBUG("PDFinfo x2:" << evt->pdf_info()->x2());
470  ATH_MSG_DEBUG("PDFinfo scalePDF:" << evt->pdf_info()->scalePDF());
471  ATH_MSG_DEBUG("PDFinfo pdf1:" << evt->pdf_info()->pdf1());
472  ATH_MSG_DEBUG("PDFinfo pdf2:" << evt->pdf_info()->pdf2());
473 #endif
474  }
475  else
476  ATH_MSG_DEBUG("No PDF information available in HepMC::GenEvent!");
477 
478  // set the randomseeds
480 
481  // fill weights
483 
484  return StatusCode::SUCCESS;
485 }
486 
489 
490  // reset
491  m_mergingWeight = m_pythia->info.mergingWeightNLO(); // first initialisation, good for the nominal weight in any case (see Merging:includeWeightInXsection).
492  m_enhanceWeight = 1.0; // better to keep the enhancement from UserHooks in a separate variable, to keep the code clear
493 
494 #ifndef PYTHIA8_304SERIES
495  // include Enhance userhook weight
496  for(const auto &hook: m_userHooksPtrs) {
497  if (hook->canEnhanceEmission()) {
498  m_enhanceWeight *= hook->getEnhancedEventWeight();
499  }
500  }
501 #endif // not PYTHIA8_304SERIES
502 
503  // DO NOT try to distinguish between phase space and merging contributions: only their product contains both, so always multiply them in order to be robust
504  double eventWeight = m_pythia->info.weight() * m_mergingWeight * m_enhanceWeight;
505 
506  std::map<std::string,double> fWeights; // temporary weight vector, needed while managing hepmc2 and 3 simultaneously
507 
508  // save as Default the usual pythia8 weight, including enhancements
509  // NB: info.weight() is built upon info.eventWeightLHEF (but includes shower "decorations"), not from the rwgt LHE vector.
510  size_t atlas_specific_weights = 1;
511  fWeights["Default"]=eventWeight;
512  if(m_internal_event_number == 1) {
513  ATH_MSG_INFO("found LHEF version "<< m_pythia->info.LHEFversion() );
514  if (m_pythia->info.LHEFversion() > 1) m_doLHE3Weights = true;
515  m_weightIDs.clear();
516  m_weightIDs.push_back("Default"); // notice that we are assigning a weightname = weightID
517  }
518 
519  // Now systematic weights: a) generator and b) shower variations
520  // merging implies reading from lhe files (while the opposite is not true, ok, but still works)
521  if (m_lheFile!="") {
522  // remember that the true merging component can be disentangled only when there are LHE input events
523 
524  // make sure to get the merging weight correctly: this is needed -only- for systematic weights (generator and parton shower variations)
525  // NB: most robust way, regardless of the specific value of Merging:includeWeightInXsection (if on, info.mergingWeight() is always 1.0 ! )
526  m_mergingWeight *= ( m_pythia->info.weight() / m_pythia->info.eventWeightLHEF ); // this is the actual merging component, regardless of conventions,
527  // needed by systematic variation weights
528  // this includes rare compensation or selection-biased weights
529  }
530 
531 
532  ATH_MSG_DEBUG("Event weights: enhancing weight, merging weight, total weight = "<<m_enhanceWeight<<", "<<m_mergingWeight<<", "<<eventWeight);
533 
534  // we already filled the "Default" weight
535  std::vector<std::string>::const_iterator id = m_weightIDs.begin()+atlas_specific_weights;
536 
537  // a) fill generator weights
538  if(m_pythia->info.getWeightsDetailedSize() != 0){
539  for(std::map<std::string, Pythia8::LHAwgt>::const_iterator wgt = m_pythia->info.rwgt->wgts.begin();
540  wgt != m_pythia->info.rwgt->wgts.end(); ++wgt){
541 
542  if(m_internal_event_number == 1){
543  m_weightIDs.push_back(wgt->first);
544  }else{
545  if(*id != wgt->first){
546  ATH_MSG_ERROR("Mismatch in LHE3 weight id. Found "<<wgt->first<<", expected "<<*id);
547  return StatusCode::FAILURE;
548  }
549  ++id;
550  }
551 
552  std::map<std::string, Pythia8::LHAweight>::const_iterator weightName = m_pythia->info.init_weights->find(wgt->first);
553  if(weightName != m_pythia->info.init_weights->end()){
554  fWeights[weightName->second.contents] = m_mergingWeight * m_enhanceWeight * wgt->second.contents;
555  }else{
556  // if no name was assigned/found, assign a weightname = weightid
557  fWeights[wgt->first] = m_mergingWeight * m_enhanceWeight * wgt->second.contents;
558  }
559 
560  }
561  }
562 
563  // b) shower weight variations
564  // always start from second shower weight: the first was saved already as "Default"
565 
566  for(int iw = 1; iw < m_pythia->info.PYTHIA8_NWEIGHTS(); ++iw){
567 
568  std::string wtName = ((int)m_showerWeightNames.size() == m_pythia->info.PYTHIA8_NWEIGHTS())? m_showerWeightNames[iw]: "ShowerWt_" +std::to_string(iw);
569 
570  fWeights[wtName] = m_mergingWeight * m_enhanceWeight * m_pythia->info.PYTHIA8_WEIGHT(iw)*m_conversion;
571 
572  if(m_internal_event_number == 1) {
573  m_weightIDs.push_back(wtName);
574  ATH_MSG_DEBUG("Shower weight name, value, conversion: "<<wtName<<", "<< fWeights[wtName] <<","<<m_conversion );
575  }
576  }
577 
578  // Now save, as a backup, also the LHEF weight without enhancements (useful for possible, rare, non-default cases)
579  // to make clear that it should not be used in analyses, save it always with a -10 factor (so that its goal is clear even if not checking its name)
580  if (m_lheFile!="")
581  {
582  fWeights["EXTRA_bare_LHE_weight"]=(-10.0)*m_pythia->info.eventWeightLHEF;
583  if(m_internal_event_number == 1) m_weightIDs.push_back("EXTRA_bare_LHE_weight");
584  }
585 
586  // Sad, but needed: create a string vector with acceptable order of weight names
587  if(m_internal_event_number == 1){
588  m_weightNames.clear();
589  for(const std::string &id : m_weightIDs){
590  if(m_doLHE3Weights) {
591  std::map<std::string, Pythia8::LHAweight>::const_iterator weight = m_pythia->info.init_weights->find(id);
592  if(weight != m_pythia->info.init_weights->end()) m_weightNames.push_back(weight->second.contents);
593  else m_weightNames.push_back(id);
594  }
595  else {
596  m_weightNames.push_back(id);
597  }
598  }
599  if (m_weightNames.size() != fWeights.size() ) {
600  ATH_MSG_ERROR("Something wrong when building list of weight names: " << m_weightNames.size() << " vs "<< fWeights.size() << ", exiting ...");
601  return StatusCode::FAILURE;
602  }
603  }
604 
605  // The following depends on the specific hepmc2/3 implementation
606 #ifdef HEPMC3
607  if (!evt->run_info()) {
608  evt->set_run_info(m_runinfo);
609  }
610  evt->run_info()->set_weight_names(m_weightNames);
611 
612  // for the first event, weight EXTRA_bare_LHE_weight is not present in evt->weights(), so we need to book a place for it by hand
613  if (m_internal_event_number == 1 && evt->run_info()->weight_names().size() == evt->weights().size()+1 ) {
614  evt->weights().push_back(1.0);
615  }
616 
617  // added conversion GeV -> MeV to ensure correct units
618  evt->set_units(HepMC3::Units::MEV, HepMC3::Units::MM);
619 
620  evt->weights().resize(fWeights.size(), 1.0);
621  for (auto w: fWeights) {
622  evt->weight(w.first)=w.second;
623  }
624 
625  auto beams=evt->beams();
626  ATH_MSG_DEBUG( " Energy of the beams " << beams[0]->momentum().e() );
627 
628 //uncomment to list HepMC events
629 // std::cout << " print::listing Pythia8 " << std::endl;
630 // HepMC3::Print::listing(std::cout, *evt);
631 
632 #else
633  evt->weights().clear();
634  for (const auto& w: m_weightNames ) {evt->weights()[w]=fWeights[w];}
635  auto beams=evt->beam_particles();
636  ATH_MSG_DEBUG( " Energy of the beams " << beams.first->momentum().e() );
637 
638 //uncomment to list HepMC events
639 // std::cout << " print::listing Pythia8 " << std::endl;
640 // evt->print();
641 #endif
642 
643  return StatusCode::SUCCESS;
644 }
645 
648 
649  ATH_MSG_INFO(">>> Pythia8_i from genFinalize");
650  m_pythia->stat();
651 
652  Pythia8::Info info = m_pythia->info;
653  double xs = info.sigmaGen(); // in mb
654 
656  const double accfactor = m_nMerged / info.nAccepted();
657  ATH_MSG_DEBUG("Multiplying cross-section by CKKWL merging acceptance of "<<m_nMerged <<"/" <<info.nAccepted() << " = " << accfactor
658  << ": " << xs << " -> " << xs*accfactor);
659  xs *= accfactor;
660  }
661 
662  if(m_doFxFxXS){
663  ATH_MSG_DEBUG("Using FxFx cross section recipe: xs = "<< m_sigmaTotal << " / " << 1e9*info.nTried());
664  xs = m_sigmaTotal / (1e9*info.nTried());
665  std::cout << "Using FxFx cross section recipe: xs = "<< m_sigmaTotal << " / " << 1e9*info.nTried() << std::endl;
666  }
667 
668  if( ! m_athenaTool.empty()){
669  double xsmod = m_athenaTool->CrossSectionScaleFactor();
670  ATH_MSG_DEBUG("Multiplying cross-section by Pythia Modifier tool factor " << xsmod );
671  xs *= xsmod;
672  }
673 
674  xs *= 1000. * 1000.;//convert to nb
675 
676  std::cout << "MetaData: cross-section (nb)= " << xs <<std::endl;
677  std::cout << "MetaData: generator= Pythia 8." << std::string(py8version()) <<std::endl;
678 
679  if(m_doLHE3Weights || m_weightIDs.size()>1 ){
680  std::cout<<"MetaData: weights = ";
681  for (const auto& w : m_weightNames ) std::cout<< w <<" | ";
682  std::cout<<std::endl;
683  }
684 
685  if (m_computeEfficiency){
686  if (info.nTried()>0) ATH_MSG_INFO("Pythia8 efficiency (nAccepted/nTried %) = " << info.nAccepted()*100./info.nTried());
687  else ATH_MSG_INFO("Pythia8 efficiency cannot be computed, nTried <=0");
688  }
689 
690 if(m_useRndmGenSvc){
691  ATH_MSG_INFO("Number of random numbers produced " << m_atlasRndmEngine->getRNCalls());
692  if (m_useReseed) ATH_MSG_INFO("Each event was reseeded ");
693  else ATH_MSG_INFO("Events were not reseeded ");
694 }
695 
696  return StatusCode::SUCCESS;
697 }
698 
701  return m_version;
702 }
703 
705 const std::string& Pythia8_i::pythia_stream() {
706  return s_pythia_stream;
707 }
708 
710 std::string Pythia8_i::xmlpath(){
711 
712 
713  std::string foundpath = "";
714 
715 // Try to find the xmldoc directory using PathResolver:
716  foundpath = PathResolverFindCalibDirectory( "Pythia8/xmldoc" );
717 
718  return foundpath;
719 }
720 
722 
723 // fix Pythia8 shower weights change in conventions
724 #ifdef PYTHIA8_NWEIGHTS
725  #undef PYTHIA8_NWEIGHTS
726  #undef PYTHIA8_WEIGHT
727  #undef PYTHIA8_WLABEL
728  #undef PYTHIA8_CONVERSION
729 #endif
Pythia8_i::NEUTRON
@ NEUTRON
Definition: Pythia8_i.h:120
read
IovVectorMap_t read(const Folder &theFolder, const SelectionCriterion &choice, const unsigned int limit=10)
Definition: openCoraCool.cxx:569
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
Pythia8_i::m_maxFailures
UnsignedIntegerProperty m_maxFailures
Definition: Pythia8_i.h:94
Pythia8_i::m_enhanceWeight
double m_enhanceWeight
Definition: Pythia8_i.h:162
Pythia8_i::m_nMerged
double m_nMerged
Definition: Pythia8_i.h:135
Pythia8_i::m_saveLHE
BooleanProperty m_saveLHE
Definition: Pythia8_i.h:172
ForwardTracker::beam1
@ beam1
Definition: ForwardTrackerConstants.h:13
Pythia8_i::fillWeights
virtual StatusCode fillWeights(HepMC::GenEvent *evt)
Definition: Pythia8_i.cxx:488
Pythia8_i::MUON
@ MUON
Definition: Pythia8_i.h:120
Pythia8_i::m_version
double m_version
Definition: Pythia8_i.h:114
Pythia8_i::m_procPtr
std::shared_ptr< Pythia8::Sigma2Process > m_procPtr
Definition: Pythia8_i.h:148
ForwardTracker::beam2
@ beam2
Definition: ForwardTrackerConstants.h:13
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
Pythia8_i::xmlpath
static std::string xmlpath()
Definition: Pythia8_i.cxx:710
Pythia8_i::m_particleDataFile
StringProperty m_particleDataFile
Definition: Pythia8_i.h:158
Pythia8_i::m_weightNames
std::vector< std::string > m_weightNames
Definition: Pythia8_i.h:164
Pythia8_i::m_useReseed
BooleanProperty m_useReseed
Definition: Pythia8_i.h:99
rerun_display.cmd
string cmd
Definition: rerun_display.py:67
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
Pythia8_i::m_failureCount
unsigned int m_failureCount
Definition: Pythia8_i.h:139
Pythia8_i::ELECTRON
@ ELECTRON
Definition: Pythia8_i.h:120
Pythia8_i::m_athenaTool
PublicToolHandle< IPythia8Custom > m_athenaTool
Definition: Pythia8_i.h:170
MM
@ MM
Definition: RegSelEnums.h:38
Pythia8_i::pythiaVersion
double pythiaVersion() const
Definition: Pythia8_i.cxx:700
Pythia8_i::m_userHooks
StringArrayProperty m_userHooks
Definition: Pythia8_i.h:102
LArG4FSStartPointFilter.evt
evt
Definition: LArG4FSStartPointFilter.py:42
Pythia8_i::m_pt0timesMPI
DoubleProperty m_pt0timesMPI
Definition: Pythia8_i.h:104
Pythia8_i.h
Pythia8_i::m_showerWeightNamesProp
StringArrayProperty m_showerWeightNamesProp
Definition: Pythia8_i.h:168
Pythia8_i::ANTIPROTON
@ ANTIPROTON
Definition: Pythia8_i.h:120
Pythia8_i::ANTINEUTRON
@ ANTINEUTRON
Definition: Pythia8_i.h:120
Pythia8_i::m_userHooksPtrs
std::vector< UserHooksPtrType > m_userHooksPtrs
Definition: Pythia8_i.h:150
PathResolverFindCalibDirectory
std::string PathResolverFindCalibDirectory(const std::string &logical_file_name)
Definition: PathResolver.cxx:288
Pythia8_i::m_doFxFxXS
BooleanProperty m_doFxFxXS
Definition: Pythia8_i.h:132
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:190
Pythia8_UserResonance::UserResonanceFactory::create
static std::shared_ptr< ResonanceWidths > create(const std::string &name, int pdgid)
Call this with the name of the ResonanceWidth and PDG ID to which it will be applied e....
Definition: UserResonanceFactory.cxx:13
Pythia8_i::m_doLHE3Weights
bool m_doLHE3Weights
Definition: Pythia8_i.h:165
GenModule
Base class for common behaviour of generator interfaces.
Definition: GenModule.h:39
Pythia8_i::m_userResonancePtrs
std::vector< std::shared_ptr< Pythia8::ResonanceWidths > > m_userResonancePtrs
Definition: Pythia8_i.h:154
Pythia8_i::fillEvt
virtual StatusCode fillEvt(HepMC::GenEvent *evt)
For filling the HepMC event object.
Definition: Pythia8_i.cxx:422
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
ParticleGun_EoverP_Config.momentum
momentum
Definition: ParticleGun_EoverP_Config.py:63
Pythia8_i::m_numberAlphaS
DoubleProperty m_numberAlphaS
Definition: Pythia8_i.h:105
McEventCollection.h
Pythia8_i::m_showerWeightNames
std::vector< std::string > m_showerWeightNames
Definition: Pythia8_i.h:167
Pythia8_UserProcess::UserProcessFactory::create
static std::shared_ptr< Sigma2Process > create(const std::string &procName)
Definition: UserProcessFactory.cxx:11
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
Pythia8_i::OXYGEN
@ OXYGEN
Definition: Pythia8_i.h:120
Pythia8_i::m_atlasRndmEngine
std::shared_ptr< customRndm > m_atlasRndmEngine
Definition: Pythia8_i.h:97
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
PYTHIA8_CONVERSION
#define PYTHIA8_CONVERSION
Definition: Pythia8_i.cxx:32
Pythia8_i::callGenerator
virtual StatusCode callGenerator()
For calling the generator on each iteration of the event loop.
Definition: Pythia8_i.cxx:357
UserResonanceFactory.h
UserProcessFactory.h
Pythia8_i::genFinalize
virtual StatusCode genFinalize()
For finalising the generator, if required.
Definition: Pythia8_i.cxx:647
Pythia8_i::m_nAccepted
double m_nAccepted
Definition: Pythia8_i.h:134
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
ppToJpsimu0mu0.SuppressSmallPT
SuppressSmallPT
Definition: ppToJpsimu0mu0.py:26
beamspotman.stat
stat
Definition: beamspotman.py:264
Pythia8_i::m_beam2
StringProperty m_beam2
Definition: Pythia8_i.h:126
Pythia8_i::m_lheFile
StringProperty m_lheFile
Definition: Pythia8_i.h:129
Pythia8_i::m_sigmaTotal
double m_sigmaTotal
Definition: Pythia8_i.h:136
Pythia8_i::m_SuppressSmallPT
Pythia8::SuppressSmallPT * m_SuppressSmallPT
Definition: Pythia8_i.h:176
python.handimod.extra
int extra
Definition: handimod.py:521
PathResolver.h
id
SG::auxid_t id
Definition: Control/AthContainers/Root/debug.cxx:239
Pythia8_i::m_weightIDs
std::vector< std::string > m_weightIDs
Definition: Pythia8_i.h:163
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
MagicNumbers.h
ActsTrk::to_string
std::string to_string(const DetectorType &type)
Definition: GeometryDefs.h:34
Pythia8_i::~Pythia8_i
~Pythia8_i()
Definition: Pythia8_i.cxx:86
Pythia8_i::HELIUM
@ HELIUM
Definition: Pythia8_i.h:120
Pythia8_i::LEAD
@ LEAD
Definition: Pythia8_i.h:120
PYTHIA8_PTRWRAP
#define PYTHIA8_PTRWRAP(A)
Definition: UserHooksFactory.h:14
RNGWrapper.h
Pythia8_i::PDGID
PDGID
Definition: Pythia8_i.h:120
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:16
Pythia8_i::Pythia8_i
Pythia8_i(const std::string &name, ISvcLocator *pSvcLocator)
Definition: Pythia8_i.cxx:68
Pythia8_i::m_userProcess
StringProperty m_userProcess
Definition: Pythia8_i.h:145
Pythia8_i::m_seeds
std::vector< long int > m_seeds
Definition: Pythia8_i.h:143
Pythia8_i::m_particleIDs
std::map< std::string, PDGID > m_particleIDs
Definition: Pythia8_i.h:141
Pythia8_i::POSITRON
@ POSITRON
Definition: Pythia8_i.h:120
Pythia8_i::m_mergingWeight
double m_mergingWeight
Definition: Pythia8_i.h:161
mc.generator
generator
Configure Herwig7 These are the commands corresponding to what would go into the regular Herwig infil...
Definition: mc.MGH7_FxFx_H71-DEFAULT_test.py:18
Pythia8_i::ANTIMUON
@ ANTIMUON
Definition: Pythia8_i.h:120
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
GenModule::m_randomSeed
IntegerProperty m_randomSeed
Seed for random number engine.
Definition: GenModule.h:84
Pythia8_i::m_outputParticleDataFile
StringProperty m_outputParticleDataFile
Definition: Pythia8_i.h:159
Pythia8_i::m_conversion
double m_conversion
Definition: Pythia8_i.h:137
Pythia8_i::m_dsid
IntegerProperty m_dsid
Definition: Pythia8_i.h:101
Pythia8_i::m_override_transform_beamenergy
bool m_override_transform_beamenergy
Definition: Pythia8_i.h:127
ATHRNG::calculateSeedsMC21
void calculateSeedsMC21(long *seeds, const std::string &algName, uint64_t ev, uint64_t run, uint64_t offset=0)
Set the random seed using a string (e.g.
Definition: RNGWrapper.cxx:37
Pythia8_i::m_sameAlphaSAsMPI
BooleanProperty m_sameAlphaSAsMPI
Definition: Pythia8_i.h:106
Pythia8_i::m_userResonances
StringProperty m_userResonances
Definition: Pythia8_i.h:152
Pythia8_i::m_doCKKWLAcceptance
BooleanProperty m_doCKKWLAcceptance
Definition: Pythia8_i.h:131
Pythia8_i::m_commands
StringArrayProperty m_commands
Definition: Pythia8_i.h:116
merge.status
status
Definition: merge.py:16
Pythia8_i::m_collisionEnergy
DoubleProperty m_collisionEnergy
Definition: Pythia8_i.h:122
Pythia8_i::m_useRndmGenSvc
BooleanProperty m_useRndmGenSvc
Definition: Pythia8_i.h:96
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:198
Pythia8_i::PROTON
@ PROTON
Definition: Pythia8_i.h:120
Pythia8_i::m_beam1
StringProperty m_beam1
Definition: Pythia8_i.h:125
Pythia8_UserHooks::UserHooksFactory::create
static UserHooks * create(const std::string &hookName)
Definition: UserHooksFactory.cxx:10
Pythia8_i::genInitialize
virtual StatusCode genInitialize()
For initializing the generator, if required.
Definition: Pythia8_i.cxx:98
python.ParticleTypeUtil.info
def info
Definition: ParticleTypeUtil.py:87
Trk::split
@ split
Definition: LayerMaterialProperties.h:38
Pythia8_i::m_internal_event_number
int m_internal_event_number
Definition: Pythia8_i.h:112
Pythia8_i::pythia_stream
static const std::string & pythia_stream()
Definition: Pythia8_i.cxx:705
Pythia8_i::m_pythia
std::unique_ptr< Pythia8::Pythia > m_pythia
Definition: Pythia8_i.h:92
Pythia8_i::m_computeEfficiency
BooleanProperty m_computeEfficiency
Definition: Pythia8_i.h:133
GenModule::getRandomEngineDuringInitialize
CLHEP::HepRandomEngine * getRandomEngineDuringInitialize(const std::string &streamName, unsigned long int randomSeedOffset, unsigned int conditionsRun=1, unsigned int lbn=1) const
Definition: GenModule.cxx:53
Pythia8_i::m_pythiaToHepMC
HepMC::Pythia8ToHepMC m_pythiaToHepMC
Definition: Pythia8_i.h:93
HepMC::set_random_states
void set_random_states(GenEvent *e, std::vector< T > a)
Definition: GenEvent.h:647