10 #include <boost/algorithm/string.hpp>
11 #include <boost/lexical_cast.hpp>
24 std::string s_pythia_stream{
"PYTHIA8_INIT"};
29 #define PYTHIA8_NWEIGHTS nWeights
30 #define PYTHIA8_WEIGHT weight
31 #define PYTHIA8_WLABEL weightLabel
32 #define PYTHIA8_CONVERSION 1.0
34 #ifdef PYTHIA_VERSION_INTEGER
35 #undef PYTHIA8_NWEIGHTS
38 #if PYTHIA_VERSION_INTEGER > 8303
39 #define PYTHIA8_NWEIGHTS nWeightGroups
41 #define PYTHIA8_NWEIGHTS nVariationGroups
43 #define PYTHIA8_WEIGHT getGroupWeight
44 #define PYTHIA8_WLABEL getGroupName
53 std::string py8version()
55 static const std::string incdir (PY8INCLUDE_DIR);
56 std::string::size_type
pos = incdir.find (
"/pythia8/");
57 if (
pos == std::string::npos)
return "";
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);
88 #ifndef PYTHIA8_3SERIES
105 m_runinfo = std::make_shared<HepMC3::GenRunInfo>();
107 struct HepMC3::GenRunInfo::ToolInfo
generator={std::string(
"Pythia8"),py8version(),std::string(
"Used generator")};
115 s_pythia_stream =
"PYTHIA8_INIT";
123 m_pythia->readString(
"Tune:pp = 5");
127 m_pythia->readString(
"PDF:pSet= LHAPDF6:cteq6l1");
130 m_pythia->readString(
"Next:numberShowEvent = 0");
133 m_pythia->settings.addFlag(
"AthenaPythia8ToHepMC:print_inconsistency",
true);
139 bool canSetHook =
true;
140 if (hook ==
"SuppressSmallPT") {
155 for(
const std::pair<const std::string, double> ¶m : Pythia8_UserHooks::UserHooksFactory::userSettings<double>()){
156 m_pythia->settings.addParm(param.first, param.second,
false,
false, 0., 0.);
159 for(
const std::pair<const std::string, int> ¶m : Pythia8_UserHooks::UserHooksFactory::userSettings<int>()){
160 m_pythia->settings.addMode(param.first, param.second,
false,
false, 0., 0.);
163 for(
const std::pair<const std::string, bool> ¶m : Pythia8_UserHooks::UserHooksFactory::userSettings<bool>()){
164 m_pythia->settings.addFlag(param.first, param.second);
167 for(
const std::pair<const std::string, std::string> ¶m : Pythia8_UserHooks::UserHooksFactory::userSettings<std::string>()){
168 m_pythia->settings.addWord(param.first, param.second);
173 ATH_MSG_ERROR(
"Unable to retrieve Athena Tool for custom Pythia processing");
174 return StatusCode::FAILURE;
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:");
190 return StatusCode::FAILURE;
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. ");
201 return StatusCode::FAILURE;
213 return StatusCode::FAILURE;
220 ATH_MSG_INFO(
" !!!!!!!!!!!! WARNING ON PYTHIA RANDOM NUMBERS !!!!!!!!!! ");
221 ATH_MSG_INFO(
" THE ATHENA SERVICE AthRNGenSvc IS USED.");
222 ATH_MSG_INFO(
" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ");
227 #if PYTHIA_VERSION_INTEGER >= 8310
232 s_pythia_stream =
"PYTHIA8";
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(
" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ");
248 std::vector<std::string> resonanceArgs;
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");
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");
266 for(std::vector<std::string>::const_iterator sId = resonanceIds.begin();
267 sId != resonanceIds.end(); ++sId){
268 int idResIn = boost::lexical_cast<int>(*sId);
273 #if PYTHIA_VERSION_INTEGER >= 8310
274 m_pythia->setResonancePtr(resonance);
276 m_pythia->setResonancePtr(resonance.get());
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!");
297 canInit = canInit &&
m_pythia->readString(
"Beams:frameType = 4");
305 canInit = canInit &&
m_pythia->readString(
"Beams:frameType = 1");
313 #if PYTHIA_VERSION_INTEGER >= 8310
328 m_pythia->settings.writeFile(
"Settings_before.log",
true);
335 returnCode = StatusCode::FAILURE;
340 m_pythia->settings.writeFile(
"Settings_after.log",
true);
348 m_pythiaToHepMC.set_print_inconsistency(
m_pythia->settings.flag(
"AthenaPythia8ToHepMC:print_inconsistency") );
363 const EventContext& ctx = Gaudi::Hive::currentContext();
378 ATH_MSG_ERROR(
"Exceeded the max number of consecutive event failures.");
379 returnCode = StatusCode::FAILURE;
389 if(stat != StatusCode::SUCCESS) returnCode =
stat;
398 double eventWeight =
m_pythia->info.mergingWeight()*
m_pythia->info.weight();
401 if(returnCode != StatusCode::FAILURE &&
409 ATH_MSG_WARNING(
"Found event weight " << eventWeight <<
" between the float and double precision limits. Rejecting event.");
436 ATH_MSG_ERROR(
"Something wrong with this event - it contains fewer than 2 particles!");
438 return StatusCode::FAILURE;
446 auto evtlhe = std::make_shared<HepMC::GenEvent>();
448 auto extra = std::make_shared<HepMC::ShortEventAttribute>(evtlhe.get());
449 evt->add_attribute(
"LHERecord",
extra);
451 ATH_MSG_WARNING(
"LHE record saving requested, but not implemented with HEPMC2 functionality.");
476 ATH_MSG_DEBUG(
"No PDF information available in HepMC::GenEvent!");
484 return StatusCode::SUCCESS;
494 #ifndef PYTHIA8_304SERIES
497 if (hook->canEnhanceEmission()) {
501 #endif // not PYTHIA8_304SERIES
506 std::map<std::string,double> fWeights;
510 size_t atlas_specific_weights = 1;
511 fWeights[
"Default"]=eventWeight;
535 std::vector<std::string>::const_iterator
id =
m_weightIDs.begin()+atlas_specific_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){
545 if(*
id != wgt->first){
546 ATH_MSG_ERROR(
"Mismatch in LHE3 weight id. Found "<<wgt->first<<
", expected "<<*
id);
547 return StatusCode::FAILURE;
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()){
566 for(
int iw = 1; iw <
m_pythia->info.PYTHIA8_NWEIGHTS(); ++iw){
582 fWeights[
"EXTRA_bare_LHE_weight"]=(-10.0)*
m_pythia->info.eventWeightLHEF;
591 std::map<std::string, Pythia8::LHAweight>::const_iterator
weight =
m_pythia->info.init_weights->find(
id);
600 ATH_MSG_ERROR(
"Something wrong when building list of weight names: " <<
m_weightNames.size() <<
" vs "<< fWeights.size() <<
", exiting ...");
601 return StatusCode::FAILURE;
607 if (!
evt->run_info()) {
608 evt->set_run_info(m_runinfo);
614 evt->weights().push_back(1.0);
620 evt->weights().resize(fWeights.size(), 1.0);
621 for (
auto w: fWeights) {
622 evt->weight(
w.first)=
w.second;
625 auto beams=
evt->beams();
633 evt->weights().clear();
635 auto beams=
evt->beam_particles();
636 ATH_MSG_DEBUG(
" Energy of the beams " << beams.first->momentum().e() );
643 return StatusCode::SUCCESS;
653 double xs =
info.sigmaGen();
657 ATH_MSG_DEBUG(
"Multiplying cross-section by CKKWL merging acceptance of "<<
m_nMerged <<
"/" <<
info.nAccepted() <<
" = " << accfactor
658 <<
": " << xs <<
" -> " << xs*accfactor);
665 std::cout <<
"Using FxFx cross section recipe: xs = "<<
m_sigmaTotal <<
" / " << 1e9*
info.nTried() << std::endl;
670 ATH_MSG_DEBUG(
"Multiplying cross-section by Pythia Modifier tool factor " << xsmod );
676 std::cout <<
"MetaData: cross-section (nb)= " << xs <<std::endl;
677 std::cout <<
"MetaData: generator= Pythia 8." << std::string(py8version()) <<std::endl;
680 std::cout<<
"MetaData: weights = ";
682 std::cout<<std::endl;
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");
696 return StatusCode::SUCCESS;
706 return s_pythia_stream;
713 std::string foundpath =
"";
724 #ifdef PYTHIA8_NWEIGHTS
725 #undef PYTHIA8_NWEIGHTS
726 #undef PYTHIA8_WEIGHT
727 #undef PYTHIA8_WLABEL
728 #undef PYTHIA8_CONVERSION