diff --git a/Gen/LbPowheg/src/component/PowhegProduction.cpp b/Gen/LbPowheg/src/component/PowhegProduction.cpp index 52298cbd073cfc8d47553eb40bbf70c47b01c71f..10c9e6c4a7a25dc13ebbd953774967ae994c1730 100644 --- a/Gen/LbPowheg/src/component/PowhegProduction.cpp +++ b/Gen/LbPowheg/src/component/PowhegProduction.cpp @@ -90,7 +90,7 @@ StatusCode PowhegProduction::hardInitialize() { if (!m_beamTool) return Error("Beam tool not initialized."); Gaudi::XYZVector pBeam1, pBeam2; m_beamTool->getMeanBeams(pBeam1, pBeam2); - stringstream eBeam1("ebeam1 "), eBeam2("ebeam2 "); + std::stringstream eBeam1("ebeam1 "), eBeam2("ebeam2 "); eBeam1 << sqrt(pBeam1.Mag2()) / Gaudi::Units::GeV; eBeam2 << sqrt(pBeam2.Mag2()) / Gaudi::Units::GeV; m_defaultSettings.push_back(eBeam1.str()); diff --git a/Gen/LbPowheg/src/component/PowhegProduction.h b/Gen/LbPowheg/src/component/PowhegProduction.h index aba65928ab83a3242fb4a76796726b68f5d96afb..20ef72b1847d3c2e7465ac177db5f91e5087381b 100644 --- a/Gen/LbPowheg/src/component/PowhegProduction.h +++ b/Gen/LbPowheg/src/component/PowhegProduction.h @@ -60,10 +60,10 @@ public: private: // Members. - CommandVector m_defaultSettings; ///< The default settings. - string m_proc; ///< The POWHEG process. - map<string, int> m_procs; ///< The POWHEG process map (int is nFinal). - Pythia8::PowhegProcs *m_powheg; ///< The actual POWHEG process container. + CommandVector m_defaultSettings; ///< The default settings. + std::string m_proc; ///< The POWHEG process. + std::map<std::string, int> m_procs; /// The POWHEG process map for nFinal. + Pythia8::PowhegProcs *m_powheg; ///< The actual POWHEG process container. }; #endif // LBPOWHEG_POWHEGPRODUCTION_H diff --git a/Gen/LbPythia8/LbPythia8/Pythia8Production.h b/Gen/LbPythia8/LbPythia8/Pythia8Production.h index 770fb3358c91987c6a6edf31723ce1399fd86380..04c83dfd7646afa4e59dc23626acca2cbbb4b4c1 100755 --- a/Gen/LbPythia8/LbPythia8/Pythia8Production.h +++ b/Gen/LbPythia8/LbPythia8/Pythia8Production.h @@ -16,8 +16,6 @@ #include "Pythia8Plugins/LHAFortran.h" #include "Pythia8Plugins/HepMC2.h" -using namespace std; - /** * Production tool to generate events with Pythia 8. * @@ -37,10 +35,10 @@ using namespace std; */ class Pythia8Production : public GaudiTool, virtual public IProductionTool { public: - typedef vector<string> CommandVector ; + typedef std::vector<std::string> CommandVector ; /// Default constructor. - Pythia8Production(const string& type, const string& name, + Pythia8Production(const std::string& type, const std::string& name, const IInterface* parent); /// Default destructor. @@ -136,7 +134,7 @@ public: Pythia8::Event m_event; ///< The Pythia 8 event record. // Members needed externally. - string m_beamToolName; ///< The name of the beam tool. + std::string m_beamToolName; ///< The name of the beam tool. protected: @@ -159,14 +157,15 @@ protected: GaudiRandomForPythia8* m_randomEngine; ///< Random number generator. int m_nEvents; ///< Number of generated events. CommandVector m_userSettings; ///< The user settings vector. - string m_tuningFile; ///< The global tuning file. - string m_tuningUserFile; ///< The user tuning file. + std::string m_tuningFile; ///< The global tuning file. + std::string m_tuningUserFile; ///< The user tuning file. bool m_validate_HEPEVT; ///< Flag to validate the event. bool m_listAllParticles; ///< Flag to list all the particles. bool m_checkParticleProperties ; ///< Flag to check particle properties. bool m_showBanner; ///< Flag to print the Pythia 8 banner. ICounterLogFile* m_xmlLogTool; ///< The XML log file. - set<unsigned int> m_special; ///< The set of special particles. + std::set<unsigned int> m_special; ///< The set of special particles. + std::set<int> m_bws; ///< Set of particles with a valid BW. }; #endif // LBPYTHIA8_PYTHIA8PRODUCTION_H diff --git a/Gen/LbPythia8/src/Lib/Pythia8Production.cpp b/Gen/LbPythia8/src/Lib/Pythia8Production.cpp index 80ce91e5daa1a0e7b412a838091f626b115693a9..d5004ef2b0e9ec2d7d27128086054b0c476ab650 100755 --- a/Gen/LbPythia8/src/Lib/Pythia8Production.cpp +++ b/Gen/LbPythia8/src/Lib/Pythia8Production.cpp @@ -30,7 +30,8 @@ //============================================================================= // Default constructor. //============================================================================= -Pythia8Production::Pythia8Production(const string& type, const string& name, +Pythia8Production::Pythia8Production(const std::string& type, + const std::string& name, const IInterface* parent) : GaudiTool(type, name, parent), m_pythia(0), m_hooks(0), m_lhaup(0), m_beamTool(0), m_pythiaBeamTool(0), m_randomEngine(0), m_nEvents(0), @@ -123,14 +124,25 @@ StatusCode Pythia8Production::initialize() { m_xmlLogTool = tool<ICounterLogFile >("XmlCounterLogFile"); // Create the Pythia 8 generator. - string xmlpath("UNKNOWN" != System::getEnv("PYTHIA8XML") ? + std::string xmlpath("UNKNOWN" != System::getEnv("PYTHIA8XML") ? System::getEnv("PYTHIA8XML") : ""); m_pythia = new Pythia8::Pythia(xmlpath, m_showBanner); if (!m_pythia) return StatusCode::FAILURE; + // Create the Breit-Wigner map for checking particle widths later. + m_bws.clear(); + int id = m_pythia->particleData.nextId(1); + while (id != 0) { + if (!m_pythia->particleData.isResonance(id)) { + m_pythia->particleData.particleDataEntryPtr(id)->initBWmass(); + if (m_pythia->particleData.useBreitWigner(id)) m_bws.insert(id); + } + id = m_pythia->particleData.nextId(id); + } + // Add LhcbHooks parameters. Pythia8::Settings &set = m_pythia->settings; - string sm("StandardModel:"), mpi("MultiPartonInteractions:"), + std::string sm("StandardModel:"), mpi("MultiPartonInteractions:"), pre("LhcbHooks:"), parm("pT0Ref"); set.addParm(pre + parm, set.parm(mpi + parm), false, false, 0, 0); parm = "ecmRef"; @@ -185,19 +197,19 @@ StatusCode Pythia8Production::initializeGenerator() { // Turn off minimum bias if using LHAup. if (m_lhaup) { - vector<string> procs; procs.push_back("SoftQCD:"); + std::vector<std::string> procs; procs.push_back("SoftQCD:"); procs.push_back("HardQCD:"); procs.push_back("Onia:"); procs.push_back("Charmonium:"); procs.push_back("Bottomonium:"); for (unsigned int proc = 0; proc < procs.size(); ++proc) { - map<string, Pythia8::FVec> fvecs = + std::map<std::string, Pythia8::FVec> fvecs = m_pythia->settings.getFVecMap(procs[proc]); - map<string, Pythia8::Flag> flags = + std::map<std::string, Pythia8::Flag> flags = m_pythia->settings.getFlagMap(procs[proc]); - for(map<string, Pythia8::FVec>::iterator itr = fvecs.begin(); + for(std::map<std::string, Pythia8::FVec>::iterator itr = fvecs.begin(); itr != fvecs.end(); ++itr) - m_pythia->settings.fvec(itr->first, vector<bool> + m_pythia->settings.fvec(itr->first, std::vector<bool> (itr->second.valNow.size(), false)); - for(map<string, Pythia8::Flag>::iterator itr = flags.begin(); + for(std::map<std::string, Pythia8::Flag>::iterator itr = flags.begin(); itr != flags.end(); ++itr) m_pythia->settings.flag(itr->first, false); } @@ -212,7 +224,7 @@ StatusCode Pythia8Production::initializeGenerator() { // Check particle properties if requested. if (m_checkParticleProperties) { - int id = m_pythia->particleData.nextId(0); + int id = m_pythia->particleData.nextId(1); while (id != 0) { if (!m_pythia->particleData.hasChanged(id)) warning() << "Data for particle with ID " << id @@ -222,6 +234,31 @@ StatusCode Pythia8Production::initializeGenerator() { } } + // Check the Breit-Wigner mass thresholds. + for (std::set<int>::iterator id = m_bws.begin(); id != m_bws.end(); ++id) { + Pythia8::ParticleDataEntry *pde = + m_pythia->particleData.particleDataEntryPtr(*id); + if (pde->isResonance()) continue; + pde->initBWmass(); + if (pde->useBreitWigner()) continue; + double mThr(0), bRatSum(0), mThrSum(0); + for (int i = 0; i < int(pde->sizeChannels()); ++i) + if (pde->channel(i).onMode() > 0) { + bRatSum += pde->channel(i).bRatio(); + double mChannelSum = 0.; + for (int j = 0; j < pde->channel(i).multiplicity(); ++j) + mChannelSum += m_pythia->particleData.m0(pde->channel(i).product(j)); + mThrSum += pde->channel(i).bRatio() * mChannelSum; + } + mThr = (bRatSum == 0.) ? 0. : mThrSum / bRatSum; + if (mThr > pde->m0()) { + warning() << "The threshold mass for particle with ID " << *id << " is " + << mThr << " GeV but its nominal mass is " << pde->m0() + << " GeV; clearing its decay channels." << endmsg; + pde->clearChannels(); + } + } + // Initialize. if (m_pythia->init()) return StatusCode::SUCCESS; else return Error("Failed to initialize Pythia 8."); @@ -236,7 +273,7 @@ StatusCode Pythia8Production::finalize() { m_pythia->stat(); // Write the cross-sections to the XML log. - vector<int> codes = m_pythia->info.codesHard(); + std::vector<int> codes = m_pythia->info.codesHard(); for (unsigned int code = 0; code < codes.size(); ++code) m_xmlLogTool->addCrossSection(m_pythia->info.nameProc(codes[code]), codes[code], @@ -265,44 +302,39 @@ StatusCode Pythia8Production::generateEvent(HepMC::GenEvent* theEvent, if (!m_pythia->flag("HadronLevel:all")) m_event = m_pythia->event; ++m_nEvents; - IDataProviderSvc* fileRecordSvc = svc<IDataProviderSvc>("FileRecordDataSvc", true); - std::string FSRName = LHCb::GenFSRLocation::Default; - LHCb::GenFSR* genFSR = getIfExists<LHCb::GenFSR>(fileRecordSvc, FSRName); - int key = 0; - - vector<int> codes = m_pythia->info.codesHard(); - - // Store the minimum bias cross-section in the GenFSR - key = LHCb::CrossSectionsFSR::CrossSectionKeyToType("MBCrossSection"); - - if(genFSR->hasGenCounter(key+100)) - { - longlong count = genFSR->getGenCounterInfo(100+key).second; + // Grab the file service. + IDataProviderSvc* fileRecordSvc = svc<IDataProviderSvc> + ("FileRecordDataSvc", true); + LHCb::GenFSR* genFSR = getIfExists<LHCb::GenFSR> + (fileRecordSvc, LHCb::GenFSRLocation::Default); + + // Store the minimum bias cross-section in the GenFSR. + std::vector<int> codes = m_pythia->info.codesHard(); + int key = LHCb::CrossSectionsFSR::CrossSectionKeyToType("MBCrossSection"); + if (genFSR->hasGenCounter(key+100)) { + longlong count = genFSR->getGenCounterInfo(100 + key).second; count = m_pythia->info.nAccepted(key) - count; - if(count > 0) genFSR->incrementGenCounter(key+100, count); - } - else if (m_pythia->info.nAccepted(key) != 0) - genFSR->addGenCounter(100+key, m_pythia->info.nAccepted(key)); - - if(genFSR->hasCrossSection(key)) genFSR->eraseCrossSection(key); - genFSR->addCrossSection(key,LHCb::GenFSR::CrossValues("Total cross-section", m_pythia->info.sigmaGen(key))); - - // Store the others cross-sections in the GenFSR - for (unsigned int code = 0; code < codes.size(); ++code) - { + if (count > 0) genFSR->incrementGenCounter(key + 100, count); + } else if (m_pythia->info.nAccepted(key) != 0) + genFSR->addGenCounter(100 + key, m_pythia->info.nAccepted(key)); + if (genFSR->hasCrossSection(key)) genFSR->eraseCrossSection(key); + genFSR->addCrossSection + (key, LHCb::GenFSR::CrossValues("Total cross-section", + m_pythia->info.sigmaGen(key))); + + // Store the others cross-sections in the GenFSR. + for (unsigned int code = 0; code < codes.size(); ++code) { key = codes[code]; - - if(genFSR->hasGenCounter(key+100)) - { - longlong count = genFSR->getGenCounterInfo(100+key).second; + if (genFSR->hasGenCounter(key + 100)) { + longlong count = genFSR->getGenCounterInfo(100 + key).second; count = m_pythia->info.nAccepted(key) - count; - if(count > 0) genFSR->incrementGenCounter(key+100, count); - } - else if (m_pythia->info.nAccepted(key) != 0) - genFSR->addGenCounter(100+key, m_pythia->info.nAccepted(key)); - - if(genFSR->hasCrossSection(key)) genFSR->eraseCrossSection(key); - genFSR->addCrossSection(key,LHCb::GenFSR::CrossValues(m_pythia->info.nameProc(key),m_pythia->info.sigmaGen(key))); + if (count > 0) genFSR->incrementGenCounter(key + 100, count); + } else if (m_pythia->info.nAccepted(key) != 0) + genFSR->addGenCounter(100 + key, m_pythia->info.nAccepted(key)); + if (genFSR->hasCrossSection(key)) genFSR->eraseCrossSection(key); + genFSR->addCrossSection + (key, LHCb::GenFSR::CrossValues(m_pythia->info.nameProc(key), + m_pythia->info.sigmaGen(key))); } // Convert the event to HepMC and return. @@ -379,7 +411,7 @@ void Pythia8Production::updateParticleProperties(const LHCb::ParticleProperty* // Create the particle if needed. int id = pythia8Id(thePP); - string name = thePP->name(); + std::string name = thePP->name(); Pythia8::ParticleData &pd = m_pythia->particleData; if (id == 0) { const LHCb::ParticleID pid = thePP->pid(); @@ -397,13 +429,12 @@ void Pythia8Production::updateParticleProperties(const LHCb::ParticleProperty* pd.m0(id, thePP->mass() / Gaudi::Units::GeV); if (id == 6 || (id >= 23 && id <= 37)) return; double lifetime = thePP->lifetime()*Gaudi::Units::c_light; - if (lifetime <= 1.e-4 * Gaudi::Units::mm || - lifetime >= 1.e16 * Gaudi::Units::mm) lifetime = 0; + if (lifetime >= 1.e16 * Gaudi::Units::mm) lifetime = 0; double width = lifetime == 0 ? 0 : Gaudi::Units::hbarc / lifetime; - if (width < 1.5e-6*Gaudi::Units::GeV) {width = 0; pd.mMin(id, 0);} - else pd.mMin(id, (thePP->mass() - thePP->maxWidth()) / Gaudi::Units::GeV); + pd.mMin(id, (thePP->mass() - (thePP->maxWidth() ? thePP->maxWidth() : + 15*width)) / Gaudi::Units::GeV); + pd.mMax(id, (thePP->mass() + 15*width) / Gaudi::Units::GeV); pd.mWidth(id, width / Gaudi::Units::GeV); - pd.mMax(id, 0); pd.tau0(id, lifetime / Gaudi::Units::mm); }