Source code for tools.nllFastWrapper

#!/usr/bin/env python3

"""
.. module:: nllFastWrapper
   :synopsis: This module provides methods to access the nllfast grid and
              compute k-factors (when available) to SUSY pair
              production cross sections.

.. moduleauthor:: Suchita Kulkarni <suchita.kulkarni@gmail.com>
.. moduleauthor:: Andre Lessa <lessa.a.p@gmail.com>
.. moduleauthor:: Wolfgang Waltenberger <wolfgang.waltenberger@gmail.com>

"""
from __future__ import print_function
import operator
import pyslha

try:
    import commands as executor
except ImportError:
    import subprocess as executor

squarks = [1000001,
           2000001,
           1000002,
           2000002,
           1000003,
           2000003,
           1000004,
           2000004]
antisquarks = list(map(operator.neg, squarks))
third = [1000005,
         2000005,
         1000006,
         2000006]
gluinos = [1000021]

import os
from smodels.tools.wrapperBase import WrapperBase
from smodels.base.smodelsLogging import logger


[docs]class NllFastWrapper(WrapperBase): """ An instance of this class represents the installation of nllfast. """ def __init__(self, sqrts, nllfastVersion, testParams, testCondition): """ :param sqrts: sqrt of s, in TeV, as an integer, :param nllfastVersion: version of the nllfast tool :param testParams: what are the test params we need to run things with? :param testCondition: the line that should be the last output line when running executable :srcPath: the path of the source code, for compilation """ WrapperBase.__init__(self) self.sqrts = int(sqrts) self.name = "nllfast%d" % sqrts self.nllfastVersion = nllfastVersion path = "<install>/smodels/lib/nllfast/nllfast-" location = path + self.nllfastVersion + "/" self.cdPath = self.absPath(location) self.executablePath = self.cdPath + "/nllfast_%dTeV" % self.sqrts self.testParams = testParams self.testCondition = testCondition self.srcPath = self.cdPath self.compiler = "gfortran" self.executable = "" def _interpolateKfactors( self, kFacsVector, xval): """ Interpolate a list of k-factor values from kFacsVector = [[x0,[k1,k2,..]], [x1,[k1,k2,..],...]. :returns: list of interpolated k-factor values at x-value xval """ import numpy kFacs = [] xpts = [x[0] for x in kFacsVector] ypts = [x[1] for x in kFacsVector] coeffs = numpy.matrix.transpose(numpy.polyfit(xpts, ypts, len(xpts) - 1)) for ik in range(len(ypts[0])): kfac = 0. for ip, coeff in enumerate(coeffs[ik]): kfac += coeff * xval ** (len(xpts) - 1 - ip) if kfac <= 0.: kfac = 1. kFacs.append(kfac) return kFacs def _getKfactorsFrom( self, output ): """ Read NLLfast output and return the k-factors. """ if not output: return False else: lines = output.split('\n') il = 0 line = lines[il] process = False while not "K_NLO" in line and il < len(lines) - 2: if "process" in line: process = line[line.find("process:") + 8:].replace(" ", "") il += 1 line = lines[il] if not process: return False # Line with header line = lines[il] # Count number of mass entries nmass = line.count('GeV') # Line with values line = lines[il + 2] data = [eval(x) for x in line.split()] if len(data) != nmass + 11: return False else: kFacs = data[9 + nmass:] return kFacs def _run ( self, this ): """ Run. Code taken from nllFast.runNLLfast Return the process name (in NLLfast notation) for the pair production of pIDs. :returns: None, if the particle ID pair is not contained in NLLfast """ current_dir = os.getcwd() try: os.chdir ( self.cdPath ) nll_output = executor.getoutput( this ) os.chdir(current_dir) return nll_output except Exception as e: ## make sure we always cd back! os.chdir(current_dir) raise e def _compute ( self, energy, pIDs, pdf, squarkmass, gluinomass ): process = self._getProcessName(pIDs) if process == "st": nll_run = "./nllfast_" + energy + " %s %s %s" % \ (process, pdf, squarkmass) else: nll_run = "./nllfast_" + energy + " %s %s %s %s" % \ (process, pdf, squarkmass, gluinomass) return self._run( nll_run ) def _getProcessName(self, pIDs): """ Return the process name (in NLLfast notation) for the pair production of pIDs. :returns: None, if the particle ID pair is not contained in NLLfast """ pid1, pid2 = sorted(pIDs) process = None # Obtain the type of process: # - gluino-gluino production = gg # - squark-antisquark = sb # - squark-squark = ss # - squark-gluino = sg # - antistop-stop # - antisbottom-sbottom = st if pid1 in antisquarks and pid2 in squarks: process = 'sb' elif abs(pid1) in squarks and abs(pid2) in squarks: process = 'ss' elif pid1 == pid2 and pid1 in gluinos: process = 'gg' elif abs(pid1) in squarks and pid2 == 1000021 or \ abs(pid2) in squarks and pid1 == 1000021: process = 'sg' elif abs(pid1) == pid2 and pid2 in third: process = 'st' return process def _getDecoupledKfactors( self, process, energy, pdf, mass ): """ Compute k-factors in the decoupled (squark or gluino) regime for the process. If a decoupled grid does not exist for the process, return None """ if process != 'sb' and process != 'gg': return None elif process == 'sb': process_dcpl = 'sdcpl' elif process == 'gg': process_dcpl = 'gdcpl' nll_run = "./nllfast_" + energy + " %s %s %s" % \ (process_dcpl, pdf, mass) e = energy.replace ( "TeV", "" ).replace ( "*", "" ) # tool = toolBox.ToolBox().get ( "nllfast%d" % int ( e ) ) nll_output = self._run(nll_run ) if "K_NLO" in nll_output: return self._getKfactorsFrom(nll_output) else: return None def _runForDecoupled ( self, energy, nllinput ): nll_run = "./nllfast_" + energy + " %s %s %s %s" % nllinput return self._run ( nll_run )
[docs] def getKfactorsFor( self, pIDs, slhafile, pdf='cteq' ): """ Read the NLLfast grid and returns a pair of k-factors (NLO and NLL) for the PIDs pair. :returns: k-factors = None, if NLLfast does not contain the process; uses the slhafile to obtain the SUSY spectrum. """ if not os.path.isfile(slhafile): logger.error("SLHA file %s not found", slhafile) return False energy = str(int(self.sqrts)) + 'TeV' # Get process name (in NLLfast notation) process = self._getProcessName(pIDs) if not process: # Return k-factors = None, if NLLfast does not have the process return (None, None) # Obtain relevant masses readfile = pyslha.readSLHAFile(slhafile) masses=readfile.blocks['MASS'] check_pids=squarks+gluinos+third for check in check_pids: if not check in masses.entries: logger.error ( "cannot compute k factor for pdgid %d: " \ " no particle mass given. will set mass to inf." % check ) masses.entries[check]=1.e10 gluinomass = abs(masses.entries[1000021]) squarkmass = sum([abs(masses.entries[pid]) for pid in squarks]) / 8. pid1, pid2 = sorted(pIDs) if pid1 in antisquarks and pid2 in squarks: squarkmass = (abs(masses.entries[abs(pid1)]) + abs(masses.entries[pid2])) / 2. elif pid1 in squarks and pid2 in squarks: squarkmass = (abs(masses.entries[pid1]) + abs(masses.entries[pid2])) / 2. elif abs(pid1) == pid2 and pid2 in third: squarkmass = abs(masses.entries[abs(pid1)]) #if tool == None: # logger.warning("No NLLfast data for sqrts = " + str(sqrts)) # return (None, None) nllpath = self.installDirectory() # self.pathOfExecutable() self.checkInstallation() nll_output = self._compute ( energy, pIDs, pdf, squarkmass, gluinomass ) # If run was successful, return k-factors: if "K_NLO" in nll_output: # NLLfast ran ok, try to get the k-factors kFacs = self._getKfactorsFrom(nll_output) if not kFacs or min(kFacs) <= 0.: logger.warning("Error obtaining k-factors") return (None, None) else: return kFacs # If run was not successful, check for decoupling error messages: elif not "too low/high" in nll_output.lower(): logger.warning("Error running NLLfast") return (None, None) # Check for decoupling cases with a decoupling grid (only for sb and gg) doDecoupling = False if "too low/high gluino" in nll_output.lower(): if gluinomass > 500. and process == 'sb': doDecoupling = True dcpl_mass = gluinomass elif "too low/high squark" in nll_output.lower(): if squarkmass > 500. and process == 'gg': doDecoupling = True dcpl_mass = squarkmass # If process do not have decoupled grids, return None: if not doDecoupling: if gluinomass < 99000. and squarkmass < 99000.: ## dont warn about obviously decoupled cases logger.warning("Masses of (q,g)=(%s,%s) out of NLLfast grid for %s, %s" % ( squarkmass, gluinomass, process, energy )) return (None, None) # Obtain k-factors from the NLLfast decoupled grid kfacs = self._getDecoupledKfactors(process,energy,pdf,min(gluinomass,squarkmass)) # Decoupling limit is satisfied, do not interpolate if not kfacs: logger.warning("Error obtaining k-factors from the NLLfast decoupled grid for " + process) return (None, None) elif dcpl_mass/min(gluinomass,squarkmass) > 10.: return kfacs # Interpolate between the non-decoupled and decoupled grids else: kFacsVector = [[10.*min(gluinomass,squarkmass),kfacs]] #First point for interpolation (decoupled grid) kfacs = None while not kfacs and dcpl_mass > 500.: dcpl_mass -= 100. # Reduce decoupled mass, until NLLfast produces results if process == 'sb': nllinput = (process, pdf, squarkmass, dcpl_mass) else: nllinput = (process, pdf, dcpl_mass, gluinomass) nll_output = self._runForDecoupled ( energy, nllinput ) kfacs = self._getKfactorsFrom(nll_output) kFacsVector.append([dcpl_mass, kfacs]) #Second point for interpolation (non-decoupled grid) if len(kFacsVector) < 2: logger.warning("Not enough points for interpolation in the decoupling " "limit") return (None, None) else: # Interpolate k-factors kFacs = self._interpolateKfactors(kFacsVector, max(squarkmass, gluinomass)) return kFacs
[docs]class NllFastWrapper7(NllFastWrapper): """ An instance of this class represents the installation of nllfast 7. """ def __init__(self): NllFastWrapper.__init__(self, 7, "1.2", testParams="gg cteq 500 600", testCondition="500. 600. 0.193E+00 " "0.450E+00 0.497E+00")
[docs]class NllFastWrapper8(NllFastWrapper): """ An instance of this class represents the installation of nllfast 8. """ def __init__(self): NllFastWrapper.__init__(self, 8, "2.1", testParams="gg cteq 500 600", testCondition="500. 600. 0.406E+00 " "0.873E+00 0.953E+00")
[docs]class NllFastWrapper13(NllFastWrapper): """ An instance of this class represents the installation of nllfast 8. """ def __init__(self): NllFastWrapper.__init__(self, 13, "3.1", testParams="gg cteq 500 600", testCondition="600. 0.394E+01 0.690E+01 " "0.731E+01 0.394E+00" )
nllFastTools = { 7 : NllFastWrapper7(), 8 : NllFastWrapper8(), 13 : NllFastWrapper13() } if __name__ == "__main__": for (sqrts, tool) in nllFastTools.items(): print("%s: installed in %s" % (tool.name, tool.installDirectory()))