Simple¶
Create a works space folder anywhere in your system orinside fitpack2. We will call such folder as
workspacein this tutorial.We need to create two files in the workspace
input.py: this file will define the scope of what is going to be fitted driver.py: this file wil serve as our main script
input.py¶
First we add some general QCD flags
#--fitting setups
conf['bootstrap'] = False
conf['flat par'] = False
conf['ftol'] = 1e-8
#--setups for DGLAP
conf['alphaSmode'] = 'backward'
#conf['alphaSmode'] = 'forward'
conf['order'] = 'NLO'
conf['Q20'] = 1.27**2
We now add addtional flags specific to DIS. In python we call DIS as IDIS to avoid clash with python’s internal library that is also called dis.
conf['tmc'] = False
conf['ht'] = False
conf['nuc'] = False
conf['offshell'] = False
conf['hq'] = False
Next we create a space for the datasets
#--datasets
conf['datasets']={}
and fill the entries with world DIS datasets
#--IDIS
conf['datasets'] = {}
conf['datasets']['idis'] = {}
conf['datasets']['idis']['xlsx'] = {}
conf['datasets']['idis']['xlsx'][10010] = 'idis/expdata/10010.xlsx' # proton | F2 | SLAC
conf['datasets']['idis']['xlsx'][10016] = 'idis/expdata/10016.xlsx' # proton | F2 | BCDMS
conf['datasets']['idis']['xlsx'][10020] = 'idis/expdata/10020.xlsx' # proton | F2 | NMC
conf['datasets']['idis']['xlsx'][10011] = 'idis/expdata/10011.xlsx' # deuteron | F2 | SLAC
conf['datasets']['idis']['xlsx'][10017] = 'idis/expdata/10017.xlsx' # deuteron | F2 | BCDMS
conf['datasets']['idis']['xlsx'][10021] = 'idis/expdata/10021.xlsx' # d/p | F2d/F2p | NMC
The datasets are located at fitpack/database and can be viewed with a spreadsheed reader. Next we specify DIS cuts
Q2cut=1.3**2
W2cut=10.0
conf['datasets']['idis']['filters']=[]
conf['datasets']['idis']['filters'].append("Q2>%f"%Q2cut)
conf['datasets']['idis']['filters'].append("W2>%f"%W2cut)
In addition we add normalization parameters to be fitted.
conf['datasets']['idis']['norm'] = {}
conf['datasets']['idis']['norm'][10010] ={'value': 1.03591e+00, 'min': 8.00000e-01, 'max': 1.20000e+00, 'fixed': False}
conf['datasets']['idis']['norm'][10016] ={'value': 9.88788e-01, 'min': 8.00000e-01, 'max': 1.20000e+00, 'fixed': False}
conf['datasets']['idis']['norm'][10020] ={'value': 1.02603e+00, 'min': 8.00000e-01, 'max': 1.20000e+00, 'fixed': False}
conf['datasets']['idis']['norm'][10011] ={'value': 1.03121e+00, 'min': 8.00000e-01, 'max': 1.20000e+00, 'fixed': False}
conf['datasets']['idis']['norm'][10017] ={'value': 1.01588e+00, 'min': 8.00000e-01, 'max': 1.20000e+00, 'fixed': False}
The min and max specify the allowed ranges for these parameters while the actual normalization uncertaninty is specified within each xlsx table. We proceed next to specify the paramters for the pdfs.
#--parameters
conf['params'] = {}
#--pdf parameters
conf['params']['pdf'] = {}
conf['params']['pdf']['g1 N'] ={'value': 3.09994e-01, 'min': None, 'max': None, 'fixed': True }
conf['params']['pdf']['g1 a'] ={'value': -5.20900e-01, 'min': -1.9, 'max': 1, 'fixed': False}
conf['params']['pdf']['g1 b'] ={'value': 4.29360e+00, 'min': 0, 'max': 10, 'fixed': False}
conf['params']['pdf']['uv1 N'] ={'value': 3.25322e-01, 'min': None, 'max': None, 'fixed': True }
conf['params']['pdf']['uv1 a'] ={'value': -2.14402e-01, 'min': -0.6, 'max': 1, 'fixed': False}
conf['params']['pdf']['uv1 b'] ={'value': 3.04406e+00, 'min': 0, 'max': 10, 'fixed': False}
conf['params']['pdf']['dv1 N'] ={'value': 1.06672e-01, 'min': None, 'max': None, 'fixed': True }
conf['params']['pdf']['dv1 a'] ={'value': -3.45404e-01, 'min': -0.6, 'max': 1, 'fixed': False}
conf['params']['pdf']['dv1 b'] ={'value': 4.48193e+00, 'min': 0, 'max': 10, 'fixed': False}
conf['params']['pdf']['db1 N'] ={'value': 3.65346e-02, 'min': 0, 'max': 1, 'fixed': False}
conf['params']['pdf']['db1 a'] ={'value': -9.35028e-01, 'min': -1, 'max': 1, 'fixed': False}
conf['params']['pdf']['db1 b'] ={'value': 4.48545e+00, 'min': 0, 'max': 10, 'fixed': False}
conf['params']['pdf']['ub1 N'] ={'value': 1.70043e-02, 'min': 0, 'max': 1, 'fixed': False}
conf['params']['pdf']['ub1 a'] ={'value': -1.00000e+00, 'min': -1, 'max': 1, 'fixed': False}
conf['params']['pdf']['ub1 b'] ={'value': 1.00000e+01, 'min': 0, 'max': 10, 'fixed': False}
conf['params']['pdf']['s1 N'] ={'value': 9.91077e-02, 'min': 0, 'max': 1, 'fixed': True}
conf['params']['pdf']['s1 a'] ={'value': 1.00000e+00, 'min': -0.6, 'max': 1, 'fixed': False}
conf['params']['pdf']['s1 b'] ={'value': 4.43290e+00, 'min': 0, 'max': 10, 'fixed': False}
conf['params']['pdf']['sb1 N'] ={'value': 2.96987e-02, 'min': 0, 'max': 1, 'fixed': False}
conf['params']['pdf']['sb1 a'] ={'value': -6.00000e-01, 'min': -0.6, 'max': 1, 'fixed': False}
conf['params']['pdf']['sb1 b'] ={'value': 3.56087e+00, 'min': 0, 'max': 10, 'fixed': False}
conf['params']['pdf']['sea1 N'] ={'value': 3.68792e-03, 'min': 0, 'max': 1, 'fixed': False}
conf['params']['pdf']['sea1 a'] ={'value': -1.87906e+00, 'min': -1.9, 'max': -1, 'fixed': False}
conf['params']['pdf']['sea1 b'] ={'value': 8.07746e+00, 'min': 0, 'max': 10, 'fixed': False}
conf['params']['pdf']['sea2 N'] ={'value': 3.68792e-03, 'min': 0, 'max': 1, 'fixed': 'sea1 N'}
conf['params']['pdf']['sea2 a'] ={'value': -1.87906e+00, 'min': -1.9, 'max': -1, 'fixed': 'sea1 a'}
conf['params']['pdf']['sea2 b'] ={'value': 8.07746e+00, 'min': 0, 'max': 10, 'fixed': 'sea1 b'}
The entries for fixed are set as follows:
False: the parameter is free to vary
True: the parameter is fixed to the numerical value in entry value
other: the parameter is fixed to be the same as another parmeter that has that entry name
driver.py¶
First make this a python executable by adding the following line at the begining
#!/usr/bin/env python
import sys,os
import numpy as np
import pylab as py
import numpy as np
import matplotlib
matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{amsmath}"]
matplotlib.rc('text',usetex=True)
import pylab as py
#--local
from tools.config import load_config,conf
from fitlib.resman import RESMAN
from fitlib import simple
Make the driver.py an excecutable
chmod +x driver.py
Then you should be able execute it as
./driver.py
If the imports does not work, it measn you have not soruces the setup.sh script.
Lets proceed to add some main routines
def main00():
np.random.seed(seed=100)
load_config('input.py')
nworkers=2
resman=RESMAN(nworkers,parallel=True)
resman.test(10)
resman.shutdown()
if __name__== "__main__":
main00()
This code will load the input.py and run the test method inside fitlib/resman.py
RESMAN is the master class controlling everything. The test evalute 10 times the residuals.
This gives you an idea for the speed at which the chi2 is evalauted.
Lets proceed to carryout some fits. For this add the following main in the code
def main01():
simple.MAXLIKE('input.py').run()
and replace the main part with
if __name__== "__main__":
main01()
The code will run showing some metrics for the chi2s, the values of the parameters etc.
Once completed it will generate an output.py which is essentially the same as the input.py
but with the parameters in the dictionaries updated.
For MC procedures, we will run the code slighly differently.
The simple example shows is very instructive when testing new parametrization, addition of new datasets or new observables.
Using the output.py it is simple to craft scripts to produce plots for the PDFs, observables etc.
For instance to examine and/or make plots for the observables that have been fitted we can craft the following
def main02():
load_config('output.py')
nworkers=2
resman=RESMAN(nworkers,parallel=True)
par=resman.parman.par
res,rres,nres=resman.get_residuals(par)
for idx in resman.idis_tabs:
tab=resman.idis_tabs[idx]
tab=pd.DataFrame(tab)
print(tab)
resman.shutdown()
Notice that we load the output.py instead of the input.py as we want to use the most update parameters from a fit.
The resman object will have the member idis_tabs provided IDIS datasets are present in the data entries of the
input file. This member is a dictionary containing each data set as a dictionary labeled by the idx variable.
Typically the tab will contain the following columns
obs: the lable for the experimental observables e.g. \(F_2, \sigma_{\rm red}\) etc
value: the quoted experimetal value of the observable
alpha: uncorrelated uncertainties added in quadrature
predictions: the theory prediction using the fitted pdfs
For each observable there are additional columns specifiying the kinematics of the observable.
For insrance, DIS tables will have the columns X and Q2.
In order to evaluate the PDFs that we have fitted we proceed as follow
def main03():
load_config('output.py')
resman=RESMAN(datasets=False)
pdf=conf['pdf']
xi=0.5
mu2=10.0
flav='i'
print(pdf.get_xF(xi,mu2,flav))
The PDF class (see qcdlib/pdf.py) is loaded from RESMAN into the global dictionary conf.
The parameters of for the PDFs are internally updated by the values in the input file.
The dictionarry conf is defined at tools/config.py. This conf dictioary is used as a global
object similar to common blocks in Fortran.