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  // Revert the recoil strategy to the old default of 1, ie. 'recoil to color'
136  // In 8.314, the default option was changed to 0, 'recoil to top' and we
137  // revert it back unless 'recoil to top'
138  // is explicitly needed for the samples
139  if (m_version > 8.313) m_pythia->readString("TimeShower:recoilStrategyRF = 1");
140 
141  // Add UserHooks first because these potentially add new settings that must exist prior to parsing commands
142 
143  for(const auto &hook: m_userHooks){
144  ATH_MSG_INFO("Adding user hook " + hook + ".");
145  bool canSetHook = true;
146  if (hook == "SuppressSmallPT") {
148  canSetHook=m_pythia->setUserHooksPtr(PYTHIA8_PTRWRAP(m_SuppressSmallPT));
149 }
150  else {
152  canSetHook = m_pythia->addUserHooksPtr(m_userHooksPtrs.back());
153 }
154  if(!canSetHook){
155  ATH_MSG_ERROR("Unable to set requested user hook.");
156  ATH_MSG_ERROR("Pythia 8 initialisation will FAIL!");
157  canInit=false;
158  }
159  }
160 
161  for(const std::pair<const std::string, double> &param : Pythia8_UserHooks::UserHooksFactory::userSettings<double>()){
162  m_pythia->settings.addParm(param.first, param.second, false, false, 0., 0.);
163  }
164 
165  for(const std::pair<const std::string, int> &param : Pythia8_UserHooks::UserHooksFactory::userSettings<int>()){
166  m_pythia->settings.addMode(param.first, param.second, false, false, 0., 0.);
167  }
168 
169  for(const std::pair<const std::string, bool> &param : Pythia8_UserHooks::UserHooksFactory::userSettings<bool>()){
170  m_pythia->settings.addFlag(param.first, param.second);
171  }
172 
173  for(const std::pair<const std::string, std::string> &param : Pythia8_UserHooks::UserHooksFactory::userSettings<std::string>()){
174  m_pythia->settings.addWord(param.first, param.second);
175  }
176 
177  if( ! m_athenaTool.empty() ){
178  if(m_athenaTool.retrieve().isFailure()){
179  ATH_MSG_ERROR("Unable to retrieve Athena Tool for custom Pythia processing");
180  return StatusCode::FAILURE;
181  }
182  else {
183  StatusCode status = m_athenaTool->InitializePythiaInfo(*m_pythia);
184  if(status != StatusCode::SUCCESS) return status;
185  }
186  }
187 
188  // Now apply the settings from the JO
189  for(const std::string &cmd : m_commands){
190 
191  if(cmd.compare("")==0) continue;
192  else if (cmd.find("Beams:id") != std::string::npos ) {
193  ATH_MSG_ERROR("With command '" << cmd << "' you are trying to set a beam different from p: please use the Beam1/Beam2 properties instead:");
194  ATH_MSG_ERROR(" example: genSeq.Pythia8.Beam1 = 'LEAD'");
195  ATH_MSG_ERROR(" genSeq.Pythia8.Beam2 = 'ANTINEUTRON'");
196  return StatusCode::FAILURE;
197  }
198  else if (cmd.find("Beams:frameType") != std::string::npos ) {
200  ATH_MSG_WARNING(" Found an explicit 'Beams:frameType' command: this will override transform beams momenta/energy parameters, regardless of its requested value. ");
201  }
202 
203  bool read = m_pythia->readString(cmd);
204 
205  if(!read){
206  ATH_MSG_ERROR("Pythia could not understand the command '"<< cmd<<"'");
207  return StatusCode::FAILURE;
208  }
209  }
210 
213 
214  ATH_MSG_INFO("Beam1 = "<<beam1);
215  ATH_MSG_INFO("Beam2 = "<<beam2);
216 
217  if(beam1 == 0 || beam2 == 0){
218  ATH_MSG_ERROR("Invalid beam particle!");
219  return StatusCode::FAILURE;
220  }
221 
222 
223 
224  if(m_useRndmGenSvc){
225 
226  ATH_MSG_INFO(" !!!!!!!!!!!! WARNING ON PYTHIA RANDOM NUMBERS !!!!!!!!!! ");
227  ATH_MSG_INFO(" THE ATHENA SERVICE AthRNGenSvc IS USED.");
228  ATH_MSG_INFO(" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ");
229 
230  m_atlasRndmEngine = std::make_shared<customRndm>();
231  CLHEP::HepRandomEngine* rndmEngine = getRandomEngineDuringInitialize(s_pythia_stream, m_randomSeed, m_dsid); // NOT THREAD-SAFE
232  m_atlasRndmEngine->init(rndmEngine);
233 #if PYTHIA_VERSION_INTEGER >= 8310
234  m_pythia->setRndmEnginePtr(m_atlasRndmEngine);
235 #else
236  m_pythia->setRndmEnginePtr(m_atlasRndmEngine.get());
237 #endif
238  s_pythia_stream = "PYTHIA8";
239  }else{
240  ATH_MSG_INFO(" !!!!!!!!!!!! WARNING ON PYTHIA RANDOM NUMBERS !!!!!!!!!! ");
241  ATH_MSG_INFO(" THE STANDARD PYTHIA8 RANDOM NUMBER SERVICE IS USED.");
242  ATH_MSG_INFO(" THE ATHENA SERVICE AthRNGSvc IS ***NOT*** USED.");
243  ATH_MSG_INFO(" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ");
244  }
245 
246  if(m_userProcess != ""){
248  }
249 
250 
251 
252  if(m_userResonances.value() != ""){
253 
254  std::vector<std::string> resonanceArgs;
255 
256  boost::split(resonanceArgs, m_userResonances.value(), boost::is_any_of(":"));
257  if(resonanceArgs.size() != 2){
258  ATH_MSG_ERROR("Cannot Understand UserResonance job option!");
259  ATH_MSG_ERROR("You should specify it as a 'name:id1,id2,id3...'");
260  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");
261  canInit = false;
262  }
263  std::vector<std::string> resonanceIds;
264  boost::split(resonanceIds, resonanceArgs.back(), boost::is_any_of(","));
265  if(resonanceIds.size()==0){
266  ATH_MSG_ERROR("You did not specifiy any PDG ids to which your user resonance width should be applied!");
267  ATH_MSG_ERROR("You should specify a list as 'name:id1,id2,id3...'");
268  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");
269  canInit=false;
270  }
271 
272  for(std::vector<std::string>::const_iterator sId = resonanceIds.begin();
273  sId != resonanceIds.end(); ++sId){
274  int idResIn = boost::lexical_cast<int>(*sId);
275  m_userResonancePtrs.push_back(Pythia8_UserResonance::UserResonanceFactory::create(resonanceArgs.front(), idResIn));
276  }
277 
278  for(std::shared_ptr<Pythia8::ResonanceWidths>& resonance : m_userResonancePtrs) {
279 #if PYTHIA_VERSION_INTEGER >= 8310
280  m_pythia->setResonancePtr(resonance);
281 #else
282  m_pythia->setResonancePtr(resonance.get());
283 #endif
284  }
285 
286  }
287 
288  if(m_particleDataFile != "") {
289  if(!m_pythia->particleData.reInit(m_particleDataFile, true)){
290  ATH_MSG_ERROR("Unable to read requested particle data table: " + m_particleDataFile + " !!");
291  ATH_MSG_ERROR("Pythia 8 initialisation will FAIL!");
292  canInit = false;
293  }
294  }
295 
296  if(m_lheFile != ""){
297  if(m_procPtr){
298  ATH_MSG_ERROR("Both LHE file and user process have been specified");
299  ATH_MSG_ERROR("LHE input does not make sense with a user process!");
300  canInit = false;
301  }
302 
303  canInit = canInit && m_pythia->readString("Beams:frameType = 4");
304  canInit = canInit && m_pythia->readString("Beams:LHEF = " + m_lheFile);
305  if(!canInit){
306  ATH_MSG_ERROR("Unable to read requested LHE file: " + m_lheFile + " !");
307  ATH_MSG_ERROR("Pythia 8 initialisation will FAIL!");
308  }
309  }else{
311  canInit = canInit && m_pythia->readString("Beams:frameType = 1");
312  canInit = canInit && m_pythia->readString("Beams:eCM = " + std::to_string(m_collisionEnergy));
313  }
314  canInit = canInit && m_pythia->readString("Beams:idA = " + std::to_string(beam1));
315  canInit = canInit && m_pythia->readString("Beams:idB = " + std::to_string(beam2));
316  }
317 
318  if(m_procPtr){
319 #if PYTHIA_VERSION_INTEGER >= 8310
320  if(!m_pythia->setSigmaPtr(m_procPtr))
321 #else
322  if(!m_pythia->setSigmaPtr(m_procPtr.get()))
323 #endif
324  {
325  ATH_MSG_ERROR("Unable to set requested user process: " + m_userProcess + " !!");
326  ATH_MSG_ERROR("Pythia 8 initialisation will FAIL!");
327  canInit = false;
328  }
329  }
330 
331  StatusCode returnCode = StatusCode::SUCCESS;
332 
333  m_pythia->particleData.listXML(m_outputParticleDataFile.value().substr(0,m_outputParticleDataFile.value().find("xml"))+"orig.xml");
334  m_pythia->settings.writeFile("Settings_before.log",true);
335 
336  if(canInit){
337  canInit = m_pythia->init();
338  }
339 
340  if(!canInit){
341  returnCode = StatusCode::FAILURE;
342  ATH_MSG_ERROR(" *** Unable to initialise Pythia !! ***");
343  }
344 
345  m_pythia->particleData.listXML(m_outputParticleDataFile);
346  m_pythia->settings.writeFile("Settings_after.log",true);
347 
348  //counter for event failures;
349  m_failureCount = 0;
350 
352 
353  // Set set_print_inconsistency to Athena corresponding flag (allowing to change it from JO)
354  m_pythiaToHepMC.set_print_inconsistency( m_pythia->settings.flag("AthenaPythia8ToHepMC:print_inconsistency") );
355 
356  m_pythiaToHepMC.set_store_pdf(true);
357  m_pythiaToHepMC.set_store_weights(false);
358 
359  return returnCode;
360 }
361 
364 
365  ATH_MSG_DEBUG(">>> Pythia8_i from callGenerator");
367  //Re-seed the random number stream
368  long seeds[7];
369  const EventContext& ctx = Gaudi::Hive::currentContext();
370  ATHRNG::calculateSeedsMC21(seeds, s_pythia_stream, ctx.eventID().event_number(), m_dsid, m_randomSeed);
371  m_atlasRndmEngine->getEngine()->setSeeds(seeds, 0); // NOT THREAD-SAFE
372  m_seeds.clear();
373  m_seeds.push_back(seeds[0]);
374  m_seeds.push_back(seeds[1]);
375  }
376 
377  bool status = m_pythia->next();
378 
379  StatusCode returnCode = StatusCode::SUCCESS;
380 
381  if(!status){
382  ++m_failureCount;
384  ATH_MSG_ERROR("Exceeded the max number of consecutive event failures.");
385  returnCode = StatusCode::FAILURE;
386  }else{
387  ATH_MSG_INFO("Event generation failed - re-trying.");
388  returnCode = this->callGenerator();
389  }
390  }
391 
392  ATH_MSG_DEBUG("Now checking with Tool.empty() second time");
393  if( ! m_athenaTool.empty() ){
394  StatusCode stat = m_athenaTool->ModifyPythiaEvent(*m_pythia);
395  if(stat != StatusCode::SUCCESS) returnCode = stat;
396  }
397 
398  m_failureCount = 0;
399 
400  m_nAccepted += 1.;
401  // some CKKWL merged events have zero weight (or unfilled event).
402  // start again with such events
403 
404  double eventWeight = m_pythia->info.mergingWeight()*m_pythia->info.weight();
405 
406 
407  if(returnCode != StatusCode::FAILURE &&
408  // Here this is a double, but it may be converted into float at some point downstream
409  (std::abs(eventWeight) < std::numeric_limits<float>::min() ||
410  m_pythia->event.size() < 2)){
411 
412  returnCode = this->callGenerator();
413  } else if ( std::abs(eventWeight) < std::numeric_limits<float>::min() &&
414  std::abs(eventWeight) > std::numeric_limits<double>::min() ){
415  ATH_MSG_WARNING("Found event weight " << eventWeight << " between the float and double precision limits. Rejecting event.");
416  } else {
417  m_nMerged += eventWeight;
419 
420  // For FxFx cross section:
421  m_sigmaTotal+=m_pythia->info.weight();
422  }
423 
424  return returnCode;
425 }
426 
429 
430  ATH_MSG_DEBUG(">>> Pythia8_i from fillEvt");
431 
432  evt->set_event_number(m_internal_event_number);
433 
434  // if using "getGroupWeight" and | lhastrategy | = 4, then need to convert mb to pb ( done otherwise when calling info.weight(), [...] )
435  if( m_internal_event_number == 1 && std::abs(m_pythia->info.lhaStrategy()) == 4 ) {
437  ATH_MSG_DEBUG(" LHA strategy needs a conversion to fix Pythia8 shower weights bug(s) equal to " << m_conversion);
438  }
439 
440 
441  if(m_pythia->event.size() < 2){
442  ATH_MSG_ERROR("Something wrong with this event - it contains fewer than 2 particles!");
443  ATH_MSG_ERROR("internal event number is "<<m_internal_event_number);
444  return StatusCode::FAILURE;
445  }
446 
448 
449  // In case we're asked, save the LHE event as an attribute. This only works in HEPMC3.
450  if(m_saveLHE){
451 #ifdef HEPMC3
452  auto evtlhe = std::make_shared<HepMC::GenEvent>();
453  m_pythiaToHepMC.fill_next_event(m_pythia->process, evtlhe.get(), m_internal_event_number, &m_pythia->info, &m_pythia->settings);
454  auto extra = std::make_shared<HepMC::ShortEventAttribute>(evtlhe.get());
455  evt->add_attribute("LHERecord", extra);
456 #else
457  ATH_MSG_WARNING("LHE record saving requested, but not implemented with HEPMC2 functionality.");
458 #endif
459  }
460 
461  // in debug mode you can check whether the pdf information is stored
462  if(evt->pdf_info()){
463 #ifdef HEPMC3
464  ATH_MSG_DEBUG("PDFinfo id1:" << evt->pdf_info()->parton_id[0]);
465  ATH_MSG_DEBUG("PDFinfo id2:" << evt->pdf_info()->parton_id[1]);
466  ATH_MSG_DEBUG("PDFinfo x1:" << evt->pdf_info()->x[0]);
467  ATH_MSG_DEBUG("PDFinfo x2:" << evt->pdf_info()->x[1]);
468  ATH_MSG_DEBUG("PDFinfo scalePDF:" << evt->pdf_info()->scale);
469  ATH_MSG_DEBUG("PDFinfo pdf1:" << evt->pdf_info()->pdf_id[0]);
470  ATH_MSG_DEBUG("PDFinfo pdf2:" << evt->pdf_info()->pdf_id[1]);
471 #else
472  ATH_MSG_DEBUG("PDFinfo id1:" << evt->pdf_info()->id1());
473  ATH_MSG_DEBUG("PDFinfo id2:" << evt->pdf_info()->id2());
474  ATH_MSG_DEBUG("PDFinfo x1:" << evt->pdf_info()->x1());
475  ATH_MSG_DEBUG("PDFinfo x2:" << evt->pdf_info()->x2());
476  ATH_MSG_DEBUG("PDFinfo scalePDF:" << evt->pdf_info()->scalePDF());
477  ATH_MSG_DEBUG("PDFinfo pdf1:" << evt->pdf_info()->pdf1());
478  ATH_MSG_DEBUG("PDFinfo pdf2:" << evt->pdf_info()->pdf2());
479 #endif
480  }
481  else
482  ATH_MSG_DEBUG("No PDF information available in HepMC::GenEvent!");
483 
484  // set the randomseeds
486 
487  // fill weights
489 
490  return StatusCode::SUCCESS;
491 }
492 
495 
496  // reset
497  m_mergingWeight = m_pythia->info.mergingWeightNLO(); // first initialisation, good for the nominal weight in any case (see Merging:includeWeightInXsection).
498  m_enhanceWeight = 1.0; // better to keep the enhancement from UserHooks in a separate variable, to keep the code clear
499 
500 #ifndef PYTHIA8_304SERIES
501  // include Enhance userhook weight
502  for(const auto &hook: m_userHooksPtrs) {
503  if (hook->canEnhanceEmission()) {
504  m_enhanceWeight *= hook->getEnhancedEventWeight();
505  }
506  }
507 #endif // not PYTHIA8_304SERIES
508 
509  // 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
510  double eventWeight = m_pythia->info.weight() * m_mergingWeight * m_enhanceWeight;
511 
512  std::map<std::string,double> fWeights; // temporary weight vector, needed while managing hepmc2 and 3 simultaneously
513 
514  // save as Default the usual pythia8 weight, including enhancements
515  // NB: info.weight() is built upon info.eventWeightLHEF (but includes shower "decorations"), not from the rwgt LHE vector.
516  size_t atlas_specific_weights = 1;
517  fWeights["Default"]=eventWeight;
518  if(m_internal_event_number == 1) {
519  ATH_MSG_INFO("found LHEF version "<< m_pythia->info.LHEFversion() );
520  if (m_pythia->info.LHEFversion() > 1) m_doLHE3Weights = true;
521  m_weightIDs.clear();
522  m_weightIDs.push_back("Default"); // notice that we are assigning a weightname = weightID
523  }
524 
525  // Now systematic weights: a) generator and b) shower variations
526  // merging implies reading from lhe files (while the opposite is not true, ok, but still works)
527  if (m_lheFile!="") {
528  // remember that the true merging component can be disentangled only when there are LHE input events
529 
530  // make sure to get the merging weight correctly: this is needed -only- for systematic weights (generator and parton shower variations)
531  // NB: most robust way, regardless of the specific value of Merging:includeWeightInXsection (if on, info.mergingWeight() is always 1.0 ! )
532  m_mergingWeight *= ( m_pythia->info.weight() / m_pythia->info.eventWeightLHEF ); // this is the actual merging component, regardless of conventions,
533  // needed by systematic variation weights
534  // this includes rare compensation or selection-biased weights
535  }
536 
537 
538  ATH_MSG_DEBUG("Event weights: enhancing weight, merging weight, total weight = "<<m_enhanceWeight<<", "<<m_mergingWeight<<", "<<eventWeight);
539 
540  // we already filled the "Default" weight
541  std::vector<std::string>::const_iterator id = m_weightIDs.begin()+atlas_specific_weights;
542 
543  // a) fill generator weights
544  if(m_pythia->info.getWeightsDetailedSize() != 0){
545  for(std::map<std::string, Pythia8::LHAwgt>::const_iterator wgt = m_pythia->info.rwgt->wgts.begin();
546  wgt != m_pythia->info.rwgt->wgts.end(); ++wgt){
547 
548  if(m_internal_event_number == 1){
549  m_weightIDs.push_back(wgt->first);
550  }else{
551  if(*id != wgt->first){
552  ATH_MSG_ERROR("Mismatch in LHE3 weight id. Found "<<wgt->first<<", expected "<<*id);
553  return StatusCode::FAILURE;
554  }
555  ++id;
556  }
557 
558  std::map<std::string, Pythia8::LHAweight>::const_iterator weightName = m_pythia->info.init_weights->find(wgt->first);
559  if(weightName != m_pythia->info.init_weights->end()){
560  fWeights[weightName->second.contents] = m_mergingWeight * m_enhanceWeight * wgt->second.contents;
561  }else{
562  // if no name was assigned/found, assign a weightname = weightid
563  fWeights[wgt->first] = m_mergingWeight * m_enhanceWeight * wgt->second.contents;
564  }
565 
566  }
567  }
568 
569  // b) shower weight variations
570  // always start from second shower weight: the first was saved already as "Default"
571 
572  for(int iw = 1; iw < m_pythia->info.PYTHIA8_NWEIGHTS(); ++iw){
573 
574  std::string wtName = ((int)m_showerWeightNames.size() == m_pythia->info.PYTHIA8_NWEIGHTS())? m_showerWeightNames[iw]: "ShowerWt_" +std::to_string(iw);
575 
576  fWeights[wtName] = m_mergingWeight * m_enhanceWeight * m_pythia->info.PYTHIA8_WEIGHT(iw)*m_conversion;
577 
578  if(m_internal_event_number == 1) {
579  m_weightIDs.push_back(wtName);
580  ATH_MSG_DEBUG("Shower weight name, value, conversion: "<<wtName<<", "<< fWeights[wtName] <<","<<m_conversion );
581  }
582  }
583 
584  // Now save, as a backup, also the LHEF weight without enhancements (useful for possible, rare, non-default cases)
585  // 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)
586  if (m_lheFile!="")
587  {
588  fWeights["EXTRA_bare_LHE_weight"]=(-10.0)*m_pythia->info.eventWeightLHEF;
589  if(m_internal_event_number == 1) m_weightIDs.push_back("EXTRA_bare_LHE_weight");
590  }
591 
592  // Sad, but needed: create a string vector with acceptable order of weight names
593  if(m_internal_event_number == 1){
594  m_weightNames.clear();
595  for(const std::string &id : m_weightIDs){
596  if(m_doLHE3Weights) {
597  std::map<std::string, Pythia8::LHAweight>::const_iterator weight = m_pythia->info.init_weights->find(id);
598  if(weight != m_pythia->info.init_weights->end()) m_weightNames.push_back(weight->second.contents);
599  else m_weightNames.push_back(id);
600  }
601  else {
602  m_weightNames.push_back(id);
603  }
604  }
605  if (m_weightNames.size() != fWeights.size() ) {
606  ATH_MSG_ERROR("Something wrong when building list of weight names: " << m_weightNames.size() << " vs "<< fWeights.size() << ", exiting ...");
607  return StatusCode::FAILURE;
608  }
609  }
610 
611  // The following depends on the specific hepmc2/3 implementation
612 #ifdef HEPMC3
613  if (!evt->run_info()) {
614  evt->set_run_info(m_runinfo);
615  }
616  evt->run_info()->set_weight_names(m_weightNames);
617 
618  // 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
619  if (m_internal_event_number == 1 && evt->run_info()->weight_names().size() == evt->weights().size()+1 ) {
620  evt->weights().push_back(1.0);
621  }
622 
623  // added conversion GeV -> MeV to ensure correct units
624  evt->set_units(HepMC3::Units::MEV, HepMC3::Units::MM);
625 
626  evt->weights().resize(fWeights.size(), 1.0);
627  for (auto w: fWeights) {
628  evt->weight(w.first)=w.second;
629  }
630 
631  auto beams=evt->beams();
632  ATH_MSG_DEBUG( " Energy of the beams " << beams[0]->momentum().e() );
633 
634 //uncomment to list HepMC events
635 // std::cout << " print::listing Pythia8 " << std::endl;
636 // HepMC3::Print::listing(std::cout, *evt);
637 
638 #else
639  evt->weights().clear();
640  for (const auto& w: m_weightNames ) {evt->weights()[w]=fWeights[w];}
641  auto beams=evt->beam_particles();
642  ATH_MSG_DEBUG( " Energy of the beams " << beams.first->momentum().e() );
643 
644 //uncomment to list HepMC events
645 // std::cout << " print::listing Pythia8 " << std::endl;
646 // evt->print();
647 #endif
648 
649  return StatusCode::SUCCESS;
650 }
651 
654 
655  ATH_MSG_INFO(">>> Pythia8_i from genFinalize");
656  m_pythia->stat();
657 
658  Pythia8::Info info = m_pythia->info;
659  double xs = info.sigmaGen(); // in mb
660 
662  const double accfactor = m_nMerged / info.nAccepted();
663  ATH_MSG_DEBUG("Multiplying cross-section by CKKWL merging acceptance of "<<m_nMerged <<"/" <<info.nAccepted() << " = " << accfactor
664  << ": " << xs << " -> " << xs*accfactor);
665  xs *= accfactor;
666  }
667 
668  if(m_doFxFxXS){
669  ATH_MSG_DEBUG("Using FxFx cross section recipe: xs = "<< m_sigmaTotal << " / " << 1e9*info.nTried());
670  xs = m_sigmaTotal / (1e9*info.nTried());
671  std::cout << "Using FxFx cross section recipe: xs = "<< m_sigmaTotal << " / " << 1e9*info.nTried() << std::endl;
672  }
673 
674  if( ! m_athenaTool.empty()){
675  double xsmod = m_athenaTool->CrossSectionScaleFactor();
676  ATH_MSG_DEBUG("Multiplying cross-section by Pythia Modifier tool factor " << xsmod );
677  xs *= xsmod;
678  }
679 
680  xs *= 1000. * 1000.;//convert to nb
681 
682  std::cout << "MetaData: cross-section (nb)= " << xs <<std::endl;
683  std::cout << "MetaData: generator= Pythia 8." << std::string(py8version()) <<std::endl;
684 
685  if(m_doLHE3Weights || m_weightIDs.size()>1 ){
686  std::cout<<"MetaData: weights = ";
687  for (const auto& w : m_weightNames ) std::cout<< w <<" | ";
688  std::cout<<std::endl;
689  }
690 
691  if (m_computeEfficiency){
692  if (info.nTried()>0) ATH_MSG_INFO("Pythia8 efficiency (nAccepted/nTried %) = " << info.nAccepted()*100./info.nTried());
693  else ATH_MSG_INFO("Pythia8 efficiency cannot be computed, nTried <=0");
694  }
695 
696 if(m_useRndmGenSvc){
697  ATH_MSG_INFO("Number of random numbers produced " << m_atlasRndmEngine->getRNCalls());
698  if (m_useReseed) ATH_MSG_INFO("Each event was reseeded ");
699  else ATH_MSG_INFO("Events were not reseeded ");
700 }
701 
702  return StatusCode::SUCCESS;
703 }
704 
707  return m_version;
708 }
709 
711 const std::string& Pythia8_i::pythia_stream() {
712  return s_pythia_stream;
713 }
714 
716 std::string Pythia8_i::xmlpath(){
717 
718 
719  std::string foundpath = "";
720 
721 // Try to find the xmldoc directory using PathResolver:
722  foundpath = PathResolverFindCalibDirectory( "Pythia8/xmldoc" );
723 
724  return foundpath;
725 }
726 
728 
729 // fix Pythia8 shower weights change in conventions
730 #ifdef PYTHIA8_NWEIGHTS
731  #undef PYTHIA8_NWEIGHTS
732  #undef PYTHIA8_WEIGHT
733  #undef PYTHIA8_WLABEL
734  #undef PYTHIA8_CONVERSION
735 #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:494
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:716
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:706
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:428
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:363
UserResonanceFactory.h
UserProcessFactory.h
Pythia8_i::genFinalize
virtual StatusCode genFinalize()
For finalising the generator, if required.
Definition: Pythia8_i.cxx:653
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:262
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:711
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