Commit 532c1ad5 authored by James William Monk's avatar James William Monk Committed by Graeme Stewart
Browse files

CKKWLAcceptance with multiple LH processes (Pythia8_i-00-11-05)

  * Allow CKKWLAcceptance=True even with multiple sub-processes (previously failed the job)

2015-05-06  James Monk	<jmonk@cern.ch>
  * New JO switch CKKWLAcceptance to turn on the CKKWL matching acceptance in the cross section (off by default)

2015-04-15  James Monk <jmonk@cern.ch>
  * Track events with zero Sudakov merging weight and take them into account for cross section (affects MG5)
  * Not clear how to handle all cases with > 1 (sub) process per input LHE file, so cause those jobs to fail.

2015-03-07  James Monk <jmonk@cern.ch>
  * Update func name on version check
  * Tested WprimeFlat
  * Tag as Pythia8_i-00-11-03

2015-01-04  James Monk	<jmonk@cern.ch>
  * Added list map from Pythia version to max allowed Tune:pp
  * Added a check that the Tune:pp value does not exceed the allowed max (c.f. Monash with pre 8.185)
  * Updated the WprimeFlat UserHook to have a WprimeFlat:EnergyMode setting.  Default is 8, but 13 is a new option
parent 4361020f
......@@ -14,9 +14,11 @@
#include "CLHEP/Random/RandFlat.h"
#include "AthenaKernel/IAtRndmGenSvc.h"
#include <stdexcept>
/**
* Author: James Monk (jmonk@hep.ucl.ac.uk)
* Author: James Monk (jmonk@cern.ch)
*/
using std::string;
......@@ -58,6 +60,13 @@ public:
Pythia8_i(const string &name, ISvcLocator *pSvcLocator);
~Pythia8_i();
class CommandException : public std::runtime_error{
public:
CommandException(const string &cmd): std::runtime_error("Cannot interpret command: " + cmd){
}
};
virtual StatusCode genInitialize();
virtual StatusCode callGenerator();
......@@ -80,6 +89,8 @@ private:
// StoreGateSvc* m_sGateService;
// string m_mcEventKey;
static std::string findValue(const std::string &command, const std::string &key);
int m_internal_event_number;
......@@ -101,6 +112,10 @@ private:
std::string m_lheFile;
bool m_doCKKWLAcceptance;
double m_nAccepted;
double m_nMerged;
unsigned int m_maxFailures;
unsigned int m_failureCount;
......@@ -126,6 +141,8 @@ private:
std::string m_particleDataFile;
std::string m_outputParticleDataFile;
static int s_allowedTunes(double version);
};
#endif
......@@ -11,6 +11,8 @@
#include <boost/algorithm/string.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/assign/std/vector.hpp>
#include <boost/foreach.hpp>
#define foreach BOOST_FOREACH
// calls to fortran routines
#include "CLHEP/Random/RandFlat.h"
......@@ -31,6 +33,8 @@ m_pythia(xmlpath()),
m_internal_event_number(0),
m_version(-1.),
m_atlasRndmEngine(0),
m_nAccepted(0.),
m_nMerged(0.),
m_failureCount(0),
m_procPtr(0),
m_userHookPtr(0)
......@@ -43,6 +47,7 @@ m_userHookPtr(0)
declareProperty("Beam1", m_beam1 = "PROTON");
declareProperty("Beam2", m_beam2 = "PROTON");
declareProperty("LHEFile", m_lheFile = "");
declareProperty("CKKWLAcceptance", m_doCKKWLAcceptance = false);
declareProperty("MaxFailures", m_maxFailures = 10);//the max number of consecutive failures
declareProperty("UserProcess", m_userProcess="");
declareProperty("UserHook", m_userHook="");
......@@ -127,12 +132,27 @@ StatusCode Pythia8_i::genInitialize() {
}
// Now apply the settings from the JO
for(vector<string>::const_iterator cmd = m_commands.begin();
cmd != m_commands.end(); ++cmd){
bool read = m_pythia.readString(*cmd);
foreach(const string &cmd, m_commands){
if(cmd.compare("")==0) continue;
try{
string val = findValue(cmd, "Tune:pp");
if(val != ""){
int tune = boost::lexical_cast<int>(val);
if(tune > s_allowedTunes(m_version)){
ATH_MSG_ERROR("Tune setting is not known to this version of Pythia:" + cmd);
return StatusCode::FAILURE;
}
}
}catch(CommandException err){
ATH_MSG_ERROR(err.what());
return StatusCode::FAILURE;
}
bool read = m_pythia.readString(cmd);
if(!read){
ATH_MSG_ERROR("Pythia could not understand the command '"<< *cmd<<"'");
ATH_MSG_ERROR("Pythia could not understand the command '"<< cmd<<"'");
return StatusCode::FAILURE;
}
}
......@@ -311,12 +331,20 @@ StatusCode Pythia8_i::callGenerator(){
m_failureCount = 0;
m_nAccepted += 1.;
// some CKKWL merged events have zero weight (or unfilled event).
// start again with such events
double eventWeight = m_pythia.info.mergingWeight()*m_pythia.info.weight();
if(returnCode != StatusCode::FAILURE &&
(fabs(m_pythia.info.mergingWeight()*m_pythia.info.weight()) < 1.e-18 ||
m_pythia.event.size() < 2))
returnCode = this->callGenerator();
(fabs(eventWeight) < 1.e-18 ||
m_pythia.event.size() < 2)){
returnCode = this->callGenerator();
}else{
m_nMerged += eventWeight;
}
++m_internal_event_number;
......@@ -395,7 +423,12 @@ StatusCode Pythia8_i::genFinalize(){
m_pythia.stat();
Pythia8::Info info = m_pythia.info;
double xs = info.sigmaGen();// in mb
double xs = info.sigmaGen(); // in mb
if(m_doCKKWLAcceptance){
ATH_MSG_DEBUG("Multiplying cross-section by CKKWL merging acceptance of "<<m_nMerged <<"/" <<info.nAccepted());
xs *= m_nMerged / info.nAccepted();
}
xs *= 1000. * 1000.;//convert to nb
std::cout << "MetaData: cross-section (nb)= " << xs <<std::endl;
......@@ -409,6 +442,21 @@ double Pythia8_i::pythiaVersion()const{
return m_version;
}
////////////////////////////////////////////////////////////////////////////////
string Pythia8_i::findValue(const string &command, const string &key){
if(command.find(key) == std::string::npos) return "";
vector<string> splits;
boost::split(splits, command, boost::is_any_of("="));
if(splits.size() != 2){
throw Pythia8_i::CommandException(command);
}
boost::erase_all(splits[1], " ");
return splits[1];
}
////////////////////////////////////////////////////////////////////////////////
string Pythia8_i::xmlpath(){
......@@ -436,3 +484,29 @@ string Pythia8_i::xmlpath(){
return foundpath;
}
////////////////////////////////////////////////////////////////////////////////
int Pythia8_i::s_allowedTunes(double version){
static map<double, int> allowedTunes;
if(allowedTunes.size()==0){
allowedTunes[0.] = 0;
allowedTunes[8.126] = 1;
allowedTunes[8.139] = 2;
allowedTunes[8.140] = 4;
allowedTunes[8.145] = 5;
allowedTunes[8.153] = 6;
allowedTunes[8.160] = 7;
allowedTunes[8.165] = 11;
allowedTunes[8.183] = 14;
allowedTunes[8.2] = 17;
allowedTunes[8.204] = 18;
allowedTunes[8.205] = 32;
}
map<double, int>::const_iterator maxTune = allowedTunes.upper_bound(version + 0.0000001);
if(maxTune != allowedTunes.begin()) --maxTune;
return maxTune->second;
}
......@@ -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 WprimeFlat;
}
......@@ -18,7 +21,7 @@ namespace Pythia8 {
public:
// Constructor.
WprimeFlat(){}
WprimeFlat(): m_energyMode("WprimeFlat:EnergyMode", 8){}
// Destructor.
~WprimeFlat(){}
......@@ -44,15 +47,41 @@ namespace Pythia8 {
double GamMRat = wRes/mRes;
double sHat = phaseSpacePtr->sHat();
double weightBW = pow2(sHat - m2Res) + pow2(sHat * GamMRat);
double m = sqrt(sHat)/8000.0;
if(m < 0.0375)
weightBW *= 121.88e-12*exp(13.0*m);
else
weightBW *= 1.0e-12*exp(18.5*m-1.4*log(m));
double m = sqrt(sHat)/((double)m_energyMode(settingsPtr)*1000.0);
switch(m_energyMode(settingsPtr)){
case 8:
if(m < 0.0375){
weightBW *= 121.88e-12*exp(13.0*m);
}else{
weightBW *= 1.0e-12*exp(18.5*m-1.4*log(m));
}
break;
case 13:
if(m < 0.023){
weightBW *= 102.77e-12*exp(11.5*m);
}else if(m < 0.231){
weightBW *= 1.0e-12*exp(16.1*m-1.2*log(m));
}else{
weightBW *= 1.8675e-16*exp(31.7*m-4.6*log(m));
}
break;
default:
throw std::runtime_error("Unknown WprimeFlat:EnergyMode - should be either 8 or 13!");
}
return weightBW;
}
private:
/// User-settable mode to set the collision energy (default is 8)
Pythia8_UserHooks::UserSetting<int> m_energyMode;
};
} // end namespace Pythia8
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment