From 2c5ad5d1969a83ac6895daac9311ae4f789034ee Mon Sep 17 00:00:00 2001 From: James William Monk <james.william.monk@cern.ch> Date: Mon, 14 Sep 2015 12:21:06 +0200 Subject: [PATCH] improved weight names, JetMatching Hook, GravFlat Hook (Pythia8_i-00-12-02) * Updated weight storage so that names instead of ids are used in the HepMC record * tagging Pythia8_i-00-12-02 2015-09-08 James Monk <jmonk@cern.ch> * Added JetMatchingMadGraphHook.cxx for FxFx matching 2015-08-02 James Monk <jmonk@cern.ch> * Updated version of GravFlat UserHook from Jan Stark for 13 TeV running 2015-06-16 James Monk <jmonk@cern.ch> * reinstate use of map syntax so that weight names are stored in HepMC 2015-06-08 James Monk <jmonk@cern.ch> * Merging branch for Pythia 8.205 with multiple LHE3 weights onto trunk 2015-05-18 James Monk <jmonk@cern.ch> * Check order of LHE3 weights on each event is the same * write MetaData in finalize giving the weight descriptions 2015-03-18 James Monk <jmonk@cern.ch> ... (Long ChangeLog diff - truncated) --- Generators/Pythia8_i/Pythia8_i/Pythia8_i.h | 13 +- Generators/Pythia8_i/cmt/requirements | 2 +- Generators/Pythia8_i/src/Pythia8_i.cxx | 143 ++++++++--- .../Pythia8_i/src/UserHooks/GravFlat.cxx | 226 ++++++++++++------ .../src/UserHooks/JetMatchingMadGraphHook.cxx | 10 + 5 files changed, 281 insertions(+), 113 deletions(-) create mode 100644 Generators/Pythia8_i/src/UserHooks/JetMatchingMadGraphHook.cxx diff --git a/Generators/Pythia8_i/Pythia8_i/Pythia8_i.h b/Generators/Pythia8_i/Pythia8_i/Pythia8_i.h index 4d2c597f2e9..bbc7bade10c 100644 --- a/Generators/Pythia8_i/Pythia8_i/Pythia8_i.h +++ b/Generators/Pythia8_i/Pythia8_i/Pythia8_i.h @@ -8,7 +8,8 @@ #include "GeneratorModules/GenModule.h" #include "Pythia8/Pythia.h" -#include "Pythia8/Pythia8ToHepMC.h" +//#include "Pythia8/../Pythia8Plugins/HepMC2.h" +#include "Pythia8Plugins/HepMC2.h" // calls to fortran routines #include "CLHEP/Random/RandFlat.h" @@ -67,7 +68,7 @@ public: CommandException(const string &cmd): std::runtime_error("Cannot interpret command: " + cmd){ } }; - + virtual StatusCode genInitialize(); virtual StatusCode callGenerator(); virtual StatusCode fillEvt(GenEvent *evt); @@ -87,9 +88,6 @@ private: static std::string xmlpath(); -// StoreGateSvc* m_sGateService; -// string m_mcEventKey; - static std::string findValue(const std::string &command, const std::string &key); int m_internal_event_number; @@ -141,8 +139,11 @@ private: std::string m_particleDataFile; std::string m_outputParticleDataFile; - static int s_allowedTunes(double version); + vector<string> m_weightIDs; + bool m_doLHE3Weights; + static int s_allowedTunes(double version); + }; #endif diff --git a/Generators/Pythia8_i/cmt/requirements b/Generators/Pythia8_i/cmt/requirements index 055de758b42..0e0b237d439 100644 --- a/Generators/Pythia8_i/cmt/requirements +++ b/Generators/Pythia8_i/cmt/requirements @@ -31,7 +31,7 @@ end_private macro_append pythia8_ifiles 'Pythia8_i.cxx UserProcessFactory.cxx UserHooksFactory.cxx UserResonanceFactory.cxx' macro_append UserProcessFiles 'UserProcesses/Sigma2qqbar2emu.cxx UserProcesses/Sigma2qqbar2lStarlBar.cxx UserProcesses/Sigma2qqbar2lStarlStarBar.cxx' -macro_append UserHookFiles 'UserHooks/WZVetoedShower.cxx UserHooks/QCDVetoedShower.cxx UserHooks/PoWHEGVetoedShower.cxx UserHooks/GravFlat.cxx UserHooks/SuppressMPI.cxx UserHooks/EnhanceMPI.cxx UserHooks/ISRVetoedShower.cxx UserHooks/PTRelVetoedShower.cxx UserHooks/WprimeFlat.cxx UserHooks/WprimeWZFlat.cxx UserHooks/main31.cxx' +macro_append UserHookFiles 'UserHooks/WZVetoedShower.cxx UserHooks/QCDVetoedShower.cxx UserHooks/PoWHEGVetoedShower.cxx UserHooks/GravFlat.cxx UserHooks/JetMatchingMadGraphHook.cxx UserHooks/SuppressMPI.cxx UserHooks/EnhanceMPI.cxx UserHooks/ISRVetoedShower.cxx UserHooks/PTRelVetoedShower.cxx UserHooks/WprimeFlat.cxx UserHooks/WprimeWZFlat.cxx UserHooks/main31.cxx' macro_append UserResonanceFiles 'UserResonances/ResonanceExcitedCI.cxx UserResonances/ResonanceLQ.cxx' apply_pattern named_dual_use_library library="Pythia8_i" files="$(pythia8_ifiles) $(UserProcessFiles) $(UserHookFiles) $(UserResonanceFiles)" diff --git a/Generators/Pythia8_i/src/Pythia8_i.cxx b/Generators/Pythia8_i/src/Pythia8_i.cxx index 06fa8146a8e..bee5ba18ef3 100644 --- a/Generators/Pythia8_i/src/Pythia8_i.cxx +++ b/Generators/Pythia8_i/src/Pythia8_i.cxx @@ -37,7 +37,8 @@ m_nAccepted(0.), m_nMerged(0.), m_failureCount(0), m_procPtr(0), -m_userHookPtr(0) +m_userHookPtr(0), +m_doLHE3Weights(false) { declareProperty("Commands", m_commands); declareProperty("UserParams", m_userParams); @@ -92,43 +93,62 @@ StatusCode Pythia8_i::genInitialize() { Pythia8_i::pythia_stream = "PYTHIA8_INIT"; // We do explicitly set tune 4C, since it is the starting point for many other tunes - // Tune 4C for pp collisions (default from 8.150 anyway) + // Tune 4C for pp collisions m_pythia.readString("Tune:pp = 5"); - // but also use LHAPDF and CTEQ6L1 unless otherwise instructed in JO - if(m_useLHAPDF){ - m_pythia.readString("PDF:useLHAPDF = on"); - m_pythia.readString("PDF:LHAPDFset = cteq6ll.LHpdf"); - }else{ - m_pythia.readString("PDF:useLHAPDF = off"); + // also use CTEQ6L1 from LHAPDF 6 by default + // can be over-written using JO + m_pythia.readString("PDF:pSet= LHAPDF6:cteq6ll.LHpdf"); + + // have to find any old-style Pythia 8.18x PDF commands and convert them + foreach(string &cmd, m_commands){ + try{ + string val = findValue(cmd, "PDF:LHAPDFset"); + if(val != ""){ + cmd = "PDF:pSet = LHAPDF6:" + val; + }else{ + val = findValue(cmd, "BeamRemnants:reconnectRange"); + if(val != ""){ + cmd = "ColourReconnection:range = " + val; + m_commands += "ColourReconnection:mode = 0"; + } + } + + }catch(Pythia8_i::CommandException err){ + ATH_MSG_ERROR(err.what()); + return StatusCode::FAILURE; + } + + if(cmd.find("PDF:useLHAPDF") != std::string::npos){ + ATH_MSG_WARNING("PDF:useLHAPDF is deprecated and ignored."); + cmd=""; + } } - - for(vector<string>::const_iterator param = m_userParams.begin(); - param != m_userParams.end(); ++param){ + + foreach(const string ¶m, m_userParams){ vector<string> splits; - boost::split(splits, *param, boost::is_any_of("=")); + boost::split(splits, param, boost::is_any_of("=")); if(splits.size() != 2){ - ATH_MSG_ERROR("Cannot interpret user param command: " + *param); + ATH_MSG_ERROR("Cannot interpret user param command: " + param); return StatusCode::FAILURE; } boost::erase_all(splits[0], " "); m_pythia.settings.addParm(splits[0], 0., false, false, 0., 0.); - m_commands+=*param; + m_commands+=param; } - for(vector<string>::const_iterator mode = m_userModes.begin(); - mode != m_userModes.end(); ++mode){ + foreach(const string &mode, m_userModes){ vector<string> splits; - boost::split(splits, *mode, boost::is_any_of("=")); + boost::split(splits, mode, boost::is_any_of("=")); if(splits.size() != 2){ - ATH_MSG_ERROR("Cannot interpret user mode command: " + *mode); + ATH_MSG_ERROR("Cannot interpret user mode command: " + mode); return StatusCode::FAILURE; } boost::erase_all(splits[0], " "); m_pythia.settings.addMode(splits[0], 0, false, false, 0, 0); - m_commands+=*mode; + m_commands+=mode; } // Now apply the settings from the JO @@ -257,30 +277,41 @@ StatusCode Pythia8_i::genInitialize() { canInit = false; } } - + if(m_lheFile != ""){ - if(!m_pythia.init(m_lheFile)){ - ATH_MSG_ERROR("Unable to read requested LHE file: " + m_lheFile + " !!"); - ATH_MSG_ERROR("Pythia 8 initialisation will FAIL!"); + if(m_procPtr != 0){ + ATH_MSG_ERROR("Both LHE file and user process have been specified"); + ATH_MSG_ERROR("LHE input does not make sense with a user process!"); canInit = false; } - }else{ - if(m_procPtr != 0){ - if(!m_pythia.setSigmaPtr(m_procPtr)){ - ATH_MSG_ERROR("Unable to set requested user process: " + m_userProcess + " !!"); - ATH_MSG_ERROR("Pythia 8 initialisation will FAIL!"); - canInit = false; - } + + canInit = canInit && m_pythia.readString("Beams:frameType = 4"); + canInit = canInit && m_pythia.readString("Beams:LHEF = " + m_lheFile); + if(!canInit){ + ATH_MSG_ERROR("Unable to read requested LHE file: " + m_lheFile + " !"); + ATH_MSG_ERROR("Pythia 8 initialisation will FAIL!"); } - if(!m_pythia.init(beam1, beam2, m_collisionEnergy)){ - ATH_MSG_ERROR("Unable to initialise Pythia beams !!"); + }else{ + canInit = canInit && m_pythia.readString("Beams:frameType = 1"); + canInit = canInit && m_pythia.readString("Beams:idA = " + boost::lexical_cast<string>(beam1)); + canInit = canInit && m_pythia.readString("Beams:idB = " + boost::lexical_cast<string>(beam2)); + canInit = canInit && m_pythia.readString("Beams:eCM = " + boost::lexical_cast<string>(m_collisionEnergy)); + } + + if(m_procPtr != 0){ + if(!m_pythia.setSigmaPtr(m_procPtr)){ + ATH_MSG_ERROR("Unable to set requested user process: " + m_userProcess + " !!"); ATH_MSG_ERROR("Pythia 8 initialisation will FAIL!"); canInit = false; } } - + StatusCode returnCode = SUCCESS; + if(canInit){ + canInit = m_pythia.init(); + } + if(!canInit){ returnCode = StatusCode::FAILURE; ATH_MSG_ERROR(" *** Unable to initialise Pythia !! ***"); @@ -390,7 +421,35 @@ StatusCode Pythia8_i::fillEvt(GenEvent *evt){ // set the event weight evt->weights().clear(); - evt->weights().push_back(eventWeight); + + vector<string>::const_iterator id = m_weightIDs.begin(); + + if(m_pythia.info.getWeightsDetailedSize() != 0){ + for(map<string, Pythia8::LHAwgt>::const_iterator wgt = m_pythia.info.rwgt->wgts.begin(); + wgt != m_pythia.info.rwgt->wgts.end(); ++wgt){ + + if(m_internal_event_number == 1){ + m_doLHE3Weights = true; + m_weightIDs.push_back(wgt->first); + }else{ + if(*id != wgt->first){/// mismatch in weight name! + ATH_MSG_ERROR("Mismatch in LHE3 weight name. Found "<<wgt->first<<", expected "<<*id); + return StatusCode::FAILURE; + } + ++id; + } + + map<string, Pythia8::LHAweight>::const_iterator weightName = m_pythia.info.init_weights->find(wgt->first); + if(weightName != m_pythia.info.init_weights->end()){ + evt->weights()[weightName->second.contents] = mergingWeight * wgt->second.contents; + }else{ + evt->weights()[wgt->first] = mergingWeight * wgt->second.contents; + } + + } + }else{ + evt->weights().push_back(eventWeight); + } // Units correction #ifdef HEPMC_HAS_UNITS @@ -434,6 +493,22 @@ StatusCode Pythia8_i::genFinalize(){ std::cout << "MetaData: cross-section (nb)= " << xs <<std::endl; std::cout << "MetaData: generator= Pythia 8." << PY8VERSION <<std::endl; + if(m_doLHE3Weights){ + + std::cout<<"MetaData: weights = "; + foreach(const string &id, m_weightIDs){ + + map<string, Pythia8::LHAweight>::const_iterator weight = m_pythia.info.init_weights->find(id); + + if(weight != m_pythia.info.init_weights->end()){ + std::cout<<weight->second.contents<<" | "; + }else{ + std::cout<<"Unknown | "; + } + } + std::cout<<std::endl; + } + return SUCCESS; } diff --git a/Generators/Pythia8_i/src/UserHooks/GravFlat.cxx b/Generators/Pythia8_i/src/UserHooks/GravFlat.cxx index f83217d1a39..c7acad93b99 100644 --- a/Generators/Pythia8_i/src/UserHooks/GravFlat.cxx +++ b/Generators/Pythia8_i/src/UserHooks/GravFlat.cxx @@ -2,9 +2,12 @@ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration */ +#include "UserSetting.h" #include "Pythia8_i/UserHooksFactory.h" #include "Pythia8/PhaseSpace.h" +#include <stdexcept> + namespace Pythia8{ class GravFlat; } @@ -12,30 +15,30 @@ namespace Pythia8{ Pythia8_UserHooks::UserHooksFactory::Creator<Pythia8::GravFlat> gravFlatCreator("GravFlat"); namespace Pythia8 { - + class GravFlat : public UserHooks { public: // Constructor. - GravFlat(){} + GravFlat(): m_energyMode("GravFlat:EnergyMode", 8){} // Destructor. ~GravFlat(){} // Allow process cross section to be modified... virtual bool canModifySigma() {return true;} - + // ...which gives access to the event at the trial level, before selection. - virtual double multiplySigmaBy(const SigmaProcess* sigmaProcessPtr, - const PhaseSpace* phaseSpacePtr, - bool /* inEvent */) { + virtual double multiplySigmaBy(const SigmaProcess* sigmaProcessPtr, + const PhaseSpace* phaseSpacePtr, + bool /* inEvent */) { // All events should be 2 -> 1, but kill them if not. - if (sigmaProcessPtr->nFinal() != 1) return 0.; + if (sigmaProcessPtr->nFinal() != 1) return 0.; // Weight cross section with BW propagator, i.e. to remove it. // (inEvent = false for initialization). - // No inEvent criteria, want weight both for cross section + // No inEvent criteria, want weight both for cross section // and MC generation. int idRes = sigmaProcessPtr->resonanceA(); double mRes = particleDataPtr->m0(idRes); @@ -43,75 +46,154 @@ namespace Pythia8 { double m2Res = mRes*mRes; double GamMRat = wRes/mRes; double sHat = phaseSpacePtr->sHat(); - + double weightBW = pow2(sHat - m2Res) + pow2(sHat * GamMRat); - + double weightPL = 1; double rH = sqrt(sHat); - - if (rH >= 50 && rH < 2000){ - - // Completely Flat M=2 TeV, Above 0.05 TeV, Below 2 TeV - // 1e5 - double p0 = 12.4176; - double p1 = -0.364611; - double p2 = 0.00293827; - double p3 = 4.3106e-07; - double p4 = -2.61662e-09; - double p5 = 1.33037e-12; - double p6 = -2.07848e-16; - - weightPL = 1/(p0+(p1*rH)+(p2*pow(rH,2))+(p3*pow(rH,3))+(p4*pow(rH,4))+(p5*pow(rH,5))+(p6*pow(rH,6))); - - } - else if (rH >= 2000 && rH < 4000){ - - // Completely Flat M=2 TeV, Above 2 TeV, Below 4 TeV - // 1e5 - double p0 = 41237.9; - double p1 = -95.9745; - double p2 = 0.0979667; - double p3 = -5.14644e-05; - double p4 = 1.45857e-08; - double p5 = -2.13169e-12; - double p6 = 1.26473e-16; - - weightPL = 1/(p0+(p1*rH)+(p2*pow(rH,2))+(p3*pow(rH,3))+(p4*pow(rH,4))+(p5*pow(rH,5))+(p6*pow(rH,6))); - - } - else if (rH >= 4000 && rH < 5500){ - - // Completely Flat M=2 TeV, Above 4 TeV, Below 5.5 TeV - // 1e5 - double p0 = 1.11676e+06; - double p1 = -925.647; - double p2 = 0.311595; - double p3 = -5.38465e-05; - double p4 = 4.92286e-09; - double p5 = -2.1452e-13; - double p6 = 2.97112e-18; - - weightPL = 1/(p0+(p1*rH)+(p2*pow(rH,2))+(p3*pow(rH,3))+(p4*pow(rH,4))+(p5*pow(rH,5))+(p6*pow(rH,6))); - + + double weight=1; + switch(m_energyMode(settingsPtr)){ + case 8: + + if (rH >= 50 && rH < 2000){ + + // Completely Flat M=2 TeV, Above 0.05 TeV, Below 2 TeV + // 1e5 + double p0 = 12.4176; + double p1 = -0.364611; + double p2 = 0.00293827; + double p3 = 4.3106e-07; + double p4 = -2.61662e-09; + double p5 = 1.33037e-12; + double p6 = -2.07848e-16; + + weightPL = 1/(p0+(p1*rH)+(p2*pow(rH,2))+(p3*pow(rH,3))+(p4*pow(rH,4))+(p5*pow(rH,5))+(p6*pow(rH,6))); + + } + else if (rH >= 2000 && rH < 4000){ + + // Completely Flat M=2 TeV, Above 2 TeV, Below 4 TeV + // 1e5 + double p0 = 41237.9; + double p1 = -95.9745; + double p2 = 0.0979667; + double p3 = -5.14644e-05; + double p4 = 1.45857e-08; + double p5 = -2.13169e-12; + double p6 = 1.26473e-16; + + weightPL = 1/(p0+(p1*rH)+(p2*pow(rH,2))+(p3*pow(rH,3))+(p4*pow(rH,4))+(p5*pow(rH,5))+(p6*pow(rH,6))); + + } + else if (rH >= 4000 && rH < 5500){ + + // Completely Flat M=2 TeV, Above 4 TeV, Below 5.5 TeV + // 1e5 + double p0 = 1.11676e+06; + double p1 = -925.647; + double p2 = 0.311595; + double p3 = -5.38465e-05; + double p4 = 4.92286e-09; + double p5 = -2.1452e-13; + double p6 = 2.97112e-18; + + weightPL = 1/(p0+(p1*rH)+(p2*pow(rH,2))+(p3*pow(rH,3))+(p4*pow(rH,4))+(p5*pow(rH,5))+(p6*pow(rH,6))); + + } + else if (rH >= 5500 && rH <= 6500){ + + // Completely Flat M=2 TeV, Above 5.5 TeV, Below 6.5 TeV + // 1e5 + double p0 = 9.70964e+13; + double p1 = -4.17142e-03; + double p2 = -3.06415e-03; + + weightPL = 1/((p0*exp(p1*rH))+(p2*rH)); + + } + + weight = weightBW * weightPL; + return weight; + + case 13: + + if (rH >= 50 && rH < 2500){ + + // Completely flat M=5 TeV, k/M=0.1; above 0.05 TeV, below 2.5 TeV + // 8e5 + double p0 = 1.06963e+01; + double p1 = -1.86534e-01; + double p2 = 9.00278e-04; + double p3 = 1.01576e-06; + double p4 = -1.29332e-09; + double p5 = 4.64752e-13; + double p6 = -5.68297e-17; + + weightPL = 1/(p0+(p1*rH)+(p2*pow(rH,2))+(p3*pow(rH,3))+(p4*pow(rH,4))+(p5*pow(rH,5))+(p6*pow(rH,6))); + + } + else if (rH >= 2500 && rH < 6000){ + + // Completely flat M=5 TeV, k/M=0.1; above 2.5 TeV, below 6.0 TeV + // 8e5 + double p0 = -1.77843e+03; + double p1 = 3.70114e+00; + double p2 = -1.06644e-03; + double p3 = 3.77245e-08; + double p4 = 2.49922e-11; + double p5 = -3.96985e-15; + double p6 = 1.88504e-19; + + weightPL = 1/(p0+(p1*rH)+(p2*pow(rH,2))+(p3*pow(rH,3))+(p4*pow(rH,4))+(p5*pow(rH,5))+(p6*pow(rH,6))); + + } + else if (rH >= 6000 && rH < 8000){ + + // Completely flat M=5 TeV, k/M=0.1; above 6.0 TeV, below 8.0 TeV + // 8e5 + double p0 = 1.43663e+03; + double p1 = 3.40467e-01; + double p2 = -7.44453e-05; + double p3 = -1.08578e-08; + double p4 = 6.74486e-13; + double p5 = 2.82952e-16; + double p6 = -2.20149e-20; + + weightPL = 1/(p0+(p1*rH)+(p2*pow(rH,2))+(p3*pow(rH,3))+(p4*pow(rH,4))+(p5*pow(rH,5))+(p6*pow(rH,6))); + + } + else if (rH >= 8000 && rH <= 10500){ + + // Completely flat M=5 TeV, k/M=0.1; above 8.0 TeV, below 10.5 TeV + // 8e5 + double p0 = 1.21245e+03; + double p1 = -1.08389e-01; + double p2 = -1.02834e-05; + double p3 = 4.11376e-12; + double p4 = 8.55312e-14; + double p5 = 6.98307e-18; + double p6 = -6.52683e-22; + + weightPL = 1/(p0+(p1*rH)+(p2*pow(rH,2))+(p3*pow(rH,3))+(p4*pow(rH,4))+(p5*pow(rH,5))+(p6*pow(rH,6))); + + } + + weight = weightBW * weightPL; + return weight; + + default: + + throw std::runtime_error("Unknown GravFlat:EnergyMode - should be either 8 or 13!"); + } - else if (rH >= 5500 && rH <= 6500){ - - // Completely Flat M=2 TeV, Above 5.5 TeV, Below 6.5 TeV - // 1e5 - double p0 = 9.70964e+13; - double p1 = -4.17142e-03; - double p2 = -3.06415e-03; - - weightPL = 1/((p0*exp(p1*rH))+(p2*rH)); - - } - - double weight = weightBW * weightPL; - - return weight; - } + private: + + /// User-settable mode to set the collision energy (default is 8) + Pythia8_UserHooks::UserSetting<int> m_energyMode; + }; - + } // end namespace Pythia8 diff --git a/Generators/Pythia8_i/src/UserHooks/JetMatchingMadGraphHook.cxx b/Generators/Pythia8_i/src/UserHooks/JetMatchingMadGraphHook.cxx new file mode 100644 index 00000000000..3bf5d1cb101 --- /dev/null +++ b/Generators/Pythia8_i/src/UserHooks/JetMatchingMadGraphHook.cxx @@ -0,0 +1,10 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "Pythia8_i/UserHooksFactory.h" +#include "Pythia8Plugins/JetMatching.h" + +namespace Pythia8{ + Pythia8_UserHooks::UserHooksFactory::Creator<Pythia8::JetMatchingMadgraph> JetMatchingMadGraphCreator("JetMatchingMadgraph"); +} -- GitLab