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);
86 #ifndef PYTHIA8_3SERIES
103 m_runinfo = std::make_shared<HepMC3::GenRunInfo>();
105 struct HepMC3::GenRunInfo::ToolInfo
generator={std::string(
"Pythia8"),py8version(),std::string(
"Used generator")};
113 s_pythia_stream =
"PYTHIA8_INIT";
121 m_pythia->readString(
"Tune:pp = 5");
125 m_pythia->readString(
"PDF:pSet= LHAPDF6:cteq6l1");
128 m_pythia->readString(
"Next:numberShowEvent = 0");
131 m_pythia->settings.addFlag(
"AthenaPythia8ToHepMC:print_inconsistency",
true);
137 bool canSetHook =
true;
138 if (hook ==
"SuppressSmallPT") {
153 for(
const std::pair<const std::string, double> ¶m : Pythia8_UserHooks::UserHooksFactory::userSettings<double>()){
154 m_pythia->settings.addParm(param.first, param.second,
false,
false, 0., 0.);
157 for(
const std::pair<const std::string, int> ¶m : Pythia8_UserHooks::UserHooksFactory::userSettings<int>()){
158 m_pythia->settings.addMode(param.first, param.second,
false,
false, 0., 0.);
161 for(
const std::pair<const std::string, bool> ¶m : Pythia8_UserHooks::UserHooksFactory::userSettings<bool>()){
162 m_pythia->settings.addFlag(param.first, param.second);
165 for(
const std::pair<const std::string, std::string> ¶m : Pythia8_UserHooks::UserHooksFactory::userSettings<std::string>()){
166 m_pythia->settings.addWord(param.first, param.second);
171 ATH_MSG_ERROR(
"Unable to retrieve Athena Tool for custom Pythia processing");
172 return StatusCode::FAILURE;
183 if(
cmd.compare(
"")==0)
continue;
184 else if (
cmd.find(
"Beams:id") != std::string::npos ) {
185 ATH_MSG_ERROR(
"With command '" <<
cmd <<
"' you are trying to set a beam different from p: please use the Beam1/Beam2 properties instead:");
188 return StatusCode::FAILURE;
190 else if (
cmd.find(
"Beams:frameType") != std::string::npos ) {
192 ATH_MSG_WARNING(
" Found an explicit 'Beams:frameType' command: this will override transform beams momenta/energy parameters, regardless of its requested value. ");
199 return StatusCode::FAILURE;
211 return StatusCode::FAILURE;
218 ATH_MSG_INFO(
" !!!!!!!!!!!! WARNING ON PYTHIA RANDOM NUMBERS !!!!!!!!!! ");
219 ATH_MSG_INFO(
" THE ATHENA SERVICE AthRNGenSvc IS USED.");
220 ATH_MSG_INFO(
" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ");
226 #if PYTHIA_VERSION_INTEGER >= 8310
231 s_pythia_stream =
"PYTHIA8";
233 ATH_MSG_INFO(
" !!!!!!!!!!!! WARNING ON PYTHIA RANDOM NUMBERS !!!!!!!!!! ");
234 ATH_MSG_INFO(
" THE STANDARD PYTHIA8 RANDOM NUMBER SERVICE IS USED.");
235 ATH_MSG_INFO(
" THE ATHENA SERVICE AthRNGSvc IS ***NOT*** USED.");
236 ATH_MSG_INFO(
" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ");
247 std::vector<std::string> resonanceArgs;
250 if(resonanceArgs.size() != 2){
251 ATH_MSG_ERROR(
"Cannot Understand UserResonance job option!");
252 ATH_MSG_ERROR(
"You should specify it as a 'name:id1,id2,id3...'");
253 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");
256 std::vector<std::string> resonanceIds;
257 boost::split(resonanceIds, resonanceArgs.back(), boost::is_any_of(
","));
258 if(resonanceIds.size()==0){
259 ATH_MSG_ERROR(
"You did not specifiy any PDG ids to which your user resonance width should be applied!");
260 ATH_MSG_ERROR(
"You should specify a list as 'name:id1,id2,id3...'");
261 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");
265 for(std::vector<std::string>::const_iterator sId = resonanceIds.begin();
266 sId != resonanceIds.end(); ++sId){
267 int idResIn = boost::lexical_cast<int>(*sId);
272 #if PYTHIA_VERSION_INTEGER >= 8310
273 m_pythia->setResonancePtr(resonance);
275 m_pythia->setResonancePtr(resonance.get());
291 ATH_MSG_ERROR(
"Both LHE file and user process have been specified");
292 ATH_MSG_ERROR(
"LHE input does not make sense with a user process!");
296 canInit = canInit &&
m_pythia->readString(
"Beams:frameType = 4");
304 canInit = canInit &&
m_pythia->readString(
"Beams:frameType = 1");
312 #if PYTHIA_VERSION_INTEGER >= 8310
327 m_pythia->settings.writeFile(
"Settings_before.log",
true);
339 m_pythia->settings.writeFile(
"Settings_after.log",
true);
347 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.");
398 double eventWeight =
m_pythia->info.mergingWeight()*
m_pythia->info.weight();
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;
466 ATH_MSG_DEBUG(
"No PDF information available in HepMC::GenEvent!");
474 return StatusCode::SUCCESS;
484 #ifndef PYTHIA8_304SERIES
487 if (hook->canEnhanceEmission()) {
491 #endif // not PYTHIA8_304SERIES
496 std::map<std::string,double> fWeights;
500 size_t atlas_specific_weights = 1;
501 fWeights[
"Default"]=eventWeight;
525 std::vector<std::string>::const_iterator
id =
m_weightIDs.begin()+atlas_specific_weights;
528 if(
m_pythia->info.getWeightsDetailedSize() != 0){
529 for(std::map<std::string, Pythia8::LHAwgt>::const_iterator wgt =
m_pythia->info.rwgt->wgts.begin();
530 wgt !=
m_pythia->info.rwgt->wgts.end(); ++wgt){
535 if(*
id != wgt->first){
536 ATH_MSG_ERROR(
"Mismatch in LHE3 weight id. Found "<<wgt->first<<
", expected "<<*
id);
537 return StatusCode::FAILURE;
542 std::map<std::string, Pythia8::LHAweight>::const_iterator weightName =
m_pythia->info.init_weights->find(wgt->first);
543 if(weightName !=
m_pythia->info.init_weights->end()){
556 for(
int iw = 1; iw <
m_pythia->info.PYTHIA8_NWEIGHTS(); ++iw){
572 fWeights[
"AUX_bare_not_for_analyses"]=(-10.0)*
m_pythia->info.eventWeightLHEF;
581 std::map<std::string, Pythia8::LHAweight>::const_iterator
weight =
m_pythia->info.init_weights->find(
id);
590 ATH_MSG_ERROR(
"Something wrong when building list of weight names: " <<
m_weightNames.size() <<
" vs "<< fWeights.size() <<
", exiting ...");
591 return StatusCode::FAILURE;
597 if (!
evt->run_info()) {
598 evt->set_run_info(m_runinfo);
604 evt->weights().push_back(1.0);
610 evt->weights().resize(fWeights.size(), 1.0);
611 for (
auto w: fWeights) {
612 evt->weight(
w.first)=
w.second;
615 auto beams=
evt->beams();
623 evt->weights().clear();
625 auto beams=
evt->beam_particles();
626 ATH_MSG_DEBUG(
" Energy of the beams " << beams.first->momentum().e() );
633 return StatusCode::SUCCESS;
643 double xs =
info.sigmaGen();
647 ATH_MSG_DEBUG(
"Multiplying cross-section by CKKWL merging acceptance of "<<
m_nMerged <<
"/" <<
info.nAccepted() <<
" = " << accfactor
648 <<
": " << xs <<
" -> " << xs*accfactor);
655 std::cout <<
"Using FxFx cross section recipe: xs = "<<
m_sigmaTotal <<
" / " << 1e9*
info.nTried() << std::endl;
660 ATH_MSG_DEBUG(
"Multiplying cross-section by Pythia Modifier tool factor " << xsmod );
666 std::cout <<
"MetaData: cross-section (nb)= " << xs <<std::endl;
667 std::cout <<
"MetaData: generator= Pythia 8." << std::string(py8version()) <<std::endl;
670 std::cout<<
"MetaData: weights = ";
672 std::cout<<std::endl;
676 if (
info.nTried()>0)
ATH_MSG_INFO(
"Pythia8 efficiency (nAccepted/nTried %) = " <<
info.nAccepted()*100./
info.nTried());
677 else ATH_MSG_INFO(
"Pythia8 efficiency cannot be computed, nTried <=0");
681 return StatusCode::SUCCESS;
688 HepMC::GenEvent *procEvent =
new HepMC::GenEvent();
695 for(
auto p: *procEvent){
700 for(
auto&
v: procEvent->vertices())
v->set_status(1);
701 auto beams=
evt->beams();
702 auto procBeams=procEvent->beams();
704 for (
auto p: procBeams[0]->end_vertex()->particles_out()) beams[0]->end_vertex()->add_particle_out(
p);
705 for (
auto p: procBeams[1]->end_vertex()->particles_out()) beams[1]->end_vertex()->add_particle_out(
p);
709 HepMC::GenEvent *procEvent =
new HepMC::GenEvent(
evt->momentum_unit(),
evt->length_unit());
716 }
catch(HepMC::PartonEndVertexException &ignoreIt){}
718 for(HepMC::GenEvent::particle_iterator
p = procEvent->particles_begin();
719 p != procEvent->particles_end(); ++
p){
723 std::vector<HepMC::GenParticle*> beams;
724 std::vector<HepMC::GenParticle*> procBeams;
725 beams.push_back(
evt->beam_particles().first);
726 beams.push_back(
evt->beam_particles().second);
728 if(beams[0]->
momentum().
pz() * procEvent->beam_particles().first->momentum().pz() > 0.){
729 procBeams.push_back(procEvent->beam_particles().first);
730 procBeams.push_back(procEvent->beam_particles().second);
732 procBeams.push_back(procEvent->beam_particles().second);
733 procBeams.push_back(procEvent->beam_particles().first);
736 std::map<const HepMC::GenVertex*, HepMC::GenVertex*> vtxCopies;
738 for(HepMC::GenEvent::vertex_const_iterator
v = procEvent->vertices_begin();
739 v != procEvent->vertices_end(); ++
v ) {
740 if(*
v == procBeams[0]->end_vertex() || *
v == procBeams[1]->end_vertex())
continue;
741 HepMC::GenVertex* vCopy =
new HepMC::GenVertex((*v)->position(), (*v)->id(), (*v)->weights());
742 vCopy->suggest_barcode(-(
evt->vertices_size()));
744 vtxCopies[*
v] = vCopy;
745 evt->add_vertex(vCopy);
748 for(HepMC::GenEvent::particle_const_iterator
p = procEvent->particles_begin();
749 p != procEvent->particles_end(); ++
p ){
750 if((*p)->is_beam())
continue;
753 pCopy->suggest_barcode(
evt->particles_size());
756 for(
size_t ii =0; ii != 2; ++ii){
757 if((*p)->production_vertex() == procBeams[ii]->end_vertex()){
758 beams[ii]->end_vertex()->add_particle_out(pCopy);
763 vit = vtxCopies.find((*p)->production_vertex());
764 if(vit != vtxCopies.end()) vit->second->add_particle_out(pCopy);
768 vit = vtxCopies.find((*p)->end_vertex());
769 if(vit != vtxCopies.end()) vit->second->add_particle_in(pCopy);
783 return s_pythia_stream;
790 std::string foundpath =
"";
801 #ifdef PYTHIA8_NWEIGHTS
802 #undef PYTHIA8_NWEIGHTS
803 #undef PYTHIA8_WEIGHT
804 #undef PYTHIA8_WLABEL
805 #undef PYTHIA8_CONVERSION