If you use the package ngl_resum, please cite doi:10.1007/JHEP09(2020)029.
To have this example working as a jupyter notebook, one needs to have the packages numpy, physt and - obviously - ngl_resum installed. The easiest way to do this is to use pip install ngl_resum
. Details may be found here: https://packaging.python.org/tutorials/installing-packages/#use-pip-for-installing
We start by importing the package ngl_resum and numpy:
import ngl_resum as ngl
import numpy as np
We start with the FourVector
class. As its name suggests, this class is used to describe fourvectors and contains some information specifically used in collider physics. To instantiate a FourVector
, we have to feed all four components of it. Let us define fvA
with energy $e$ energyFvA
and momenta $p_i$ iMomFvA
:
energyFvA=6.9
xMomFvA=4.2
yMomFvA=3.5
zMomFvA=1.2
fvA=ngl.FourVector(energyFvA,xMomFvA,yMomFvA,zMomFvA)
print(fvA)
[6.9,4.2,3.5,1.2]
Of course, we can access the four individual coordinates:
print("energy: ",fvA.e,"\nx-mom.: ",fvA.px,\
"\ny-mom.: ",fvA.py,"\nz-mom.: ",fvA.pz)
energy: 6.9 x-mom.: 4.2 y-mom.: 3.5 z-mom.: 1.2
In some cases, it may become useful to have the momentum vector displayed as a numpy-array:
print("numpy-array: ",fvA.vec)
numpy-array: [6.9 4.2 3.5 1.2]
The two angles $\theta$ (angle between the z-axis or beam-axis)and $\phi$ (angle in the x-y-plane) of the three-vector are attributes:
print("theta: ",fvA.theta,"\nphi: ",fvA.phi)
theta: 1.3547308176908472 phi: 0.6947382761967031
We can also access the mass and the velocity of the particle:
print("mass: ",fvA.m,"\nvelocity: ",fvA.beta)
mass: 4.034848200366403 velocity: 0.811205911255451
The length of the spatial vector $\sqrt{p_x^2+p_y^2+p_z^2}$ can also be accessed:
print("np.sqrt(px*px+py*py+pz*pz): ",fvA.absSpace)
np.sqrt(px*px+py*py+pz*pz): 5.597320787662612
Collider-physics specific attributes are the transverse energy $E_T$, the transverse momentum $p_T$, the rapidity $y=\frac{1}{2}\ln\frac{e+p_z}{e-p_z}$ and the pseudorapidity $\eta=-\ln\left(\tan\frac{\theta}{2}\right)$:
print("transverse energy: ",fvA.eT,"\ntransverse momentum: ",fvA.pT,\
"\nrapidity: ",fvA.rap,"\npseudorapidity: ",fvA.pseudorap)
transverse energy: 6.7395647606579505 transverse momentum: 5.4671747731346585 rapidity: 0.1756989434189443 pseudorapidity: 0.2177665445807071
Let us instantiate another FourVector
:
fvB=ngl.FourVector(5,4,3,0)
We can now add and subtract the fourvectors with the usual operators:
print("fvA+fvB:",fvA+fvB,"\nfvA-fvB: ",fvA-fvB)
fvA+fvB: [11.9,8.2,6.5,1.2] fvA-fvB: [1.9000000000000004,0.20000000000000018,0.5,1.2]
The multiplication operator $*$ can be used to give a scalar product of two fourvectors or as the multiplication of a scalar. Note that we use te mostly-minus metric (+ - - -):
print("fvA*fvB:",fvA*fvB,"\n10*fvA: ",10*fvA)
fvA*fvB: 7.199999999999999 10*fvA: [69.0,42.0,35.0,12.0]
The division by a scalar does work, too:
print("fvA/10:",fvA/10)
fvA/10: [0.6900000000000001,0.42000000000000004,0.35,0.12]
To measure the squared angular distance $\Delta R^2=\Delta\phi^2+\Delta\eta^2=|\phi_A-\phi_B|^2+|\eta_A-\eta_B|^2$ between two fourvectors, we can use
print("deltaR^2:",fvA.R2(fvB)," or ",fvB.R2(fvA))
deltaR^2: 0.050047515262147006 or 0.050047515262147006
while the cosine of the spatial angle between two fourvectors can be accessed by
print("cosTheta:",fvA.costheta(fvB)," or ",fvB.costheta(fvA))
cosTheta: 0.9754666932856004 or 0.9754666932856004
We can check, whether the FourVector
is massive or massless. A FourVector
a
is treated as massive (or time-like), if a*a
is larger than $10^{-7}$, otherwise it is treated as massless (or light-like):
fvA.isMassive()
True
fvA.isMassless()
False
fvB.isMassive()
False
fvB.isMassless()
True
We can check whether two FourVectors are the same. We consider two FourVectors a
and b
to be the same, if (a-b).e^2+(a-b).px^2+(a-b).py^2+(a-b).pz^2
is smaller than $10^{-10}$ to account for rounding errors. To apply this check, we can use
fvA.isSame(fvB)
False
fvA.isSame(fvA+ngl.FourVector(0.000000001,0.000000001,0.000000001,0))
True
One last method of the FourVector
class is going to be the tensor product fvA$_\mu$fvB$_\nu$ (which is probably not going to be used that often), given by
print("fvA.metric.fvB:\n",fvA.tensorProd(fvB))
fvA.metric.fvB: [[ 34.5 -27.6 -20.7 -0. ] [ 21. -16.8 -12.6 -0. ] [ 17.5 -14. -10.5 -0. ] [ 6. -4.8 -3.6 -0. ]]
The Boost
class takes care of the boosting procedure as described in Section 3.1 doi:10.1007/JHEP09(2020)029. It takes two arbitrary momentum vectors from the lab frame and creates the boost from lab frame to the frame where these two fourvectors are back-to-back alongside the z-axis. For transparency we take the first two vectors of the event from Appendix A of above article, as given in (A.3):
p1=ngl.FourVector(504.7,125.6,82.44,-450.4)
u1=p1/p1.e
u2=ngl.FourVector(1,0,0,-1)
Now to get the boost accounting for the boost from the lab frame to the frame where u1
is alongside the positive z-axis and u2
alongside the negative z-axis, we just have to instantiate a Boost
with the two fourvectors as arguments:
bst=ngl.Boost(u1,u2)
We have access to each single transformation $X$, $B$ and $Z$ as defined in Section 3.1 and thoroughly explained with an example in the pages after (A.3) of doi:10.1007/JHEP09(2020)029. Note, that each of these matrices are 4x4 numpy arrays.
We start out with the boost $X$, which is a rotation that puts the added two initial fourvectors along the x-axis:
np.dot(bst.X,(u1+u2).vec)
array([2.00000000e+00, 1.91568102e+00, 0.00000000e+00, 2.18575158e-16])
Now we apply the boost $B$, which removes the spatial component of above fourvector:
np.dot(bst.B,np.dot(bst.X,(u1+u2).vec))
array([ 5.74600946e-01, -8.88178420e-16, 0.00000000e+00, 2.18575158e-16])
Finally, we apply a second rotation $Z$ that puts the two initial vectors alongside the z-axis, with u1
in the positive direction:
np.dot(bst.Z,np.dot(bst.B,np.dot(bst.X,u1.vec)))
array([ 3.87360275e-01, -3.53883589e-16, -8.32667268e-17, 1.87240671e-01])
np.dot(bst.Z,np.dot(bst.B,np.dot(bst.X,u2.vec)))
array([ 1.87240671e-01, -4.68375339e-16, 0.00000000e+00, -1.87240671e-01])
As we can see, the two fourvectors are now nicely aligned back-to-back. The full boost is also stored in an attribute:
u1prime=np.dot(bst.LABtoCMS,u1.vec)
u1prime
array([3.87360275e-01, 0.00000000e+00, 1.11022302e-16, 1.87240671e-01])
The same goes for the inverse boost:
u1BoostedBack=np.dot(bst.CMStoLAB,u1prime)
u1BoostedBack
array([ 1. , 0.24886071, 0.16334456, -0.89241133])
u1BoostedBack-u1.vec
array([ 4.44089210e-16, -1.11022302e-16, 2.77555756e-17, -2.22044605e-16])
Of course, the numpy arrays are not easy to handle as we have to keep in mind which variable is a FourVector
and which one a numpy array. We have a shortcut to erase this problem. We can take the boost and apply it on any FourVector
and get back the FourVector
in the new frame as follows:
fvu1prime=bst.boostLABtoCMS(u1)
fvu1prime
[0.38736027472048606,0.0,1.1102230246251565e-16,0.18724067086957508]
fvu1BoostedBack=bst.boostCMStoLAB(fvu1prime)
fvu1BoostedBack
[1.0000000000000004,0.24886070933227647,0.16334456112542106,-0.8924113334654252]
fvu1BoostedBack.isSame(u1)
True
The Hist
class acts as an adapter to the physt package and is immensly based on it. It accounts for the histograms $R(t)$ that is the result of the resummation (see (4.3) of doi:10.1007/JHEP09(2020)029).
When initializing a Hist
, we at least need to provide a number of bins nbins
and a maximal value for $t$, tmax
:
nbins=10
tmax=0.1
hst=ngl.Hist(nbins, tmax)
hst
t | entries --------|---------- 0.0050 | 0.000000 0.0150 | 0.000000 0.0250 | 0.000000 0.0350 | 0.000000 0.0450 | 0.000000 0.0550 | 0.000000 0.0650 | 0.000000 0.0750 | 0.000000 0.0850 | 0.000000 0.0950 | 0.000000
Another possibility is to also calculate an error estimate of each bin:
hstErr=ngl.Hist(nbins,tmax,errorHistCalc=True)
hstErr
t | entries | error --------|----------|--------- 0.0050 | 0.000000 | 0.000000 0.0150 | 0.000000 | 0.000000 0.0250 | 0.000000 | 0.000000 0.0350 | 0.000000 | 0.000000 0.0450 | 0.000000 | 0.000000 0.0550 | 0.000000 | 0.000000 0.0650 | 0.000000 | 0.000000 0.0750 | 0.000000 | 0.000000 0.0850 | 0.000000 | 0.000000 0.0950 | 0.000000 | 0.000000
Compared to the rest of the classes we switch the order and postpone the discussion of the attributes to after the methods and operations. This is due to the fact that the methods are mainly used to populate the histograms to avoid looking at a wall of zeroes in the discussion of the attributes.
For the sake of streamlinedness, we will stop discussing the more involved case of the error estimation at this point. The discussion thereof is quite involved and provides little insight. We will look at the extraction of the error from the histogram later.
We have two functions that can be used to set a whole histogram to zero or to one:
hst.setOne()
hst
t | entries --------|---------- 0.0050 | 1.000000 0.0150 | 1.000000 0.0250 | 1.000000 0.0350 | 1.000000 0.0450 | 1.000000 0.0550 | 1.000000 0.0650 | 1.000000 0.0750 | 1.000000 0.0850 | 1.000000 0.0950 | 1.000000
hst.setZero()
hst
t | entries --------|---------- 0.0050 | 0.000000 0.0150 | 0.000000 0.0250 | 0.000000 0.0350 | 0.000000 0.0450 | 0.000000 0.0550 | 0.000000 0.0650 | 0.000000 0.0750 | 0.000000 0.0850 | 0.000000 0.0950 | 0.000000
To populate the histogram, we can use addToBin
by specifying the value of t
at which we add a weight w
:
tVal=0.053
w=0.696969
hst.addToBin(tVal,w)
hst
t | entries --------|---------- 0.0050 | 0.000000 0.0150 | 0.000000 0.0250 | 0.000000 0.0350 | 0.000000 0.0450 | 0.000000 0.0550 | 0.696969 0.0650 | 0.000000 0.0750 | 0.000000 0.0850 | 0.000000 0.0950 | 0.000000
Let us start with a disclaimer - we assume the two histograms to be initialized with the same number of bin nbins
and the same maximal $t$ value tmax
. If this is not the case, anything might happen.
To show the use of some operators, we will populate two histograms with some random numbers:
hst1=ngl.Hist(nbins, tmax)
hst2=ngl.Hist(nbins, tmax)
for i in range(0,50):
hst1.addToBin(tmax*np.random.random_sample(),np.random.random_sample())
hst2.addToBin(tmax*np.random.random_sample(),np.random.random_sample())
Let us have a look at how these historgams are populated:
hst1
t | entries --------|---------- 0.0050 | 2.102598 0.0150 | 5.474357 0.0250 | 3.765161 0.0350 | 1.790246 0.0450 | 0.792018 0.0550 | 2.373341 0.0650 | 1.990888 0.0750 | 1.367947 0.0850 | 4.934581 0.0950 | 1.853173
hst2
t | entries --------|---------- 0.0050 | 2.918347 0.0150 | 3.740859 0.0250 | 3.559869 0.0350 | 1.760655 0.0450 | 2.113326 0.0550 | 4.627729 0.0650 | 2.237888 0.0750 | 0.591951 0.0850 | 4.630542 0.0950 | 1.203527
We can add and subtract the histograms. This sums (or subtracts) the entry of each bin:
hst1+hst2
t | entries --------|---------- 0.0050 | 5.020945 0.0150 | 9.215216 0.0250 | 7.325029 0.0350 | 3.550901 0.0450 | 2.905344 0.0550 | 7.001070 0.0650 | 4.228776 0.0750 | 1.959898 0.0850 | 9.565123 0.0950 | 3.056700
hst1-hst2
t | entries --------|---------- 0.0050 | -0.815749 0.0150 | 1.733499 0.0250 | 0.205292 0.0350 | 0.029590 0.0450 | -1.321308 0.0550 | -2.254388 0.0650 | -0.247000 0.0750 | 0.775997 0.0850 | 0.304039 0.0950 | 0.649647
The multiplication of two Hists multiplies the entry of each bin, and we can also multiply with a scalar (which multiplicates each entry by the scalar):
hst1*hst2
t | entries --------|---------- 0.0050 | 6.136110 0.0150 | 20.478796 0.0250 | 13.403477 0.0350 | 3.152005 0.0450 | 1.673791 0.0550 | 10.983180 0.0650 | 4.455385 0.0750 | 0.809758 0.0850 | 22.849784 0.0950 | 2.230344
10*hst1
t | entries --------|---------- 0.0050 | 21.025980 0.0150 | 54.743571 0.0250 | 37.651608 0.0350 | 17.902456 0.0450 | 7.920177 0.0550 | 23.733412 0.0650 | 19.908881 0.0750 | 13.679475 0.0850 | 49.345810 0.0950 | 18.531734
Division by a scalar is possible as well:
hst1/10
t | entries --------|---------- 0.0050 | 0.210260 0.0150 | 0.547436 0.0250 | 0.376516 0.0350 | 0.179025 0.0450 | 0.079202 0.0550 | 0.237334 0.0650 | 0.199089 0.0750 | 0.136795 0.0850 | 0.493458 0.0950 | 0.185317
Let us finally look at how to access the data in Hist
. The Histograms certainly knows about the number of bins and maximal $t$:
print("nbins: ",hst1.nbins,"\ntmax: ",hst1.tmax)
nbins: 10 tmax: 0.1
To access the entries of the histogram, we can do so:
hst1.entries
array([2.10259796, 5.47435712, 3.76516075, 1.79024562, 0.7920177 , 2.37334117, 1.99088811, 1.36794745, 4.93458098, 1.85317337])
The bin values can be read out as well. As it is sometimes useful to have the lower or upper bin boundary or the central value, we have created access to all of them:
print("central bin values: ",hst1.centerBinValue)
print("lower bin boundary: ",hst1.lowerBinBoundary)
print("upper bin boundary: ",hst1.upperBinBoundary)
central bin values: [0.005, 0.015, 0.025, 0.035, 0.045, 0.055, 0.065, 0.07500000000000001, 0.08499999999999999, 0.095] lower bin boundary: [0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09] upper bin boundary: [0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1]
Now let us have a quick glance at the error estimations of the bins. Note, that while the representation of the histogram comes with the error itself, due to intricacies of the error computation when multiplying histograms (which is used in the showering procedure), we are only able to access the squared of the error estimate. To illustrate this, let us first populize the hstErr
from above and show its representation:
for i in range(0,50):
hstErr.addToBin(tmax*np.random.random_sample(),np.random.random_sample())
hstErr
t | entries | error --------|----------|--------- 0.0050 | 1.295799 | 0.997647 0.0150 | 0.542311 | 0.527699 0.0250 | 1.156891 | 0.837240 0.0350 | 3.134448 | 1.461181 0.0450 | 1.987160 | 1.212088 0.0550 | 1.799335 | 1.103540 0.0650 | 3.308246 | 1.490670 0.0750 | 4.888158 | 1.768580 0.0850 | 5.766296 | 2.122759 0.0950 | 0.334684 | 0.296820
Now we will access the squared error of the bins:
hstErr.squaredError
array([0.99530047, 0.27846598, 0.70097064, 2.1350506 , 1.46915656, 1.21780011, 2.2220969 , 3.127874 , 4.50610368, 0.08810237])
This might seem a little bit odd. For more insight into the error handling we refer to the documentation of the example codes.
Now we get to the core of the package. An instance of Event
contains all the reevant information to start one showering. We will not go too deep into details here, but mainly refer to the two example codes. To instantiate an Event
, one can either feed it
feedDipole
parameter), orpylhe.LHEEvent
read-in via pylhe.readLHE
(using the eventFromFile
parameter)To each of those we have an example code.While we will not explain every attribute and method of this class in detail, we still want to give an overview of some intricacies.
First of all we want to discuss the feedDipole feature. Via Event(feedDipole=dipole)
we can set up the showering of one single dipole, consisting of two fourVectors in an array:
leg1=ngl.FourVector(1,0,0,0.5)
leg2=ngl.FourVector(1,0,0,-0.5)
dipole=[leg1,leg2]
evDip=ngl.Event(feedDipole=dipole)
This feature is straightforward enough. We have set up this dipole for showering.
The other feature of Event
is the more intricate showering of an event read in from a .lhe-file. After feeding the pylhe.LHEEvent
into the parameter eventFromFile
we have two additional options, namely
An example where we shower both the dipoles formed by the incoming-intermediate particles and the intermediate-outgoing particles is given for example in Section 5 of doi:10.1007/JHEP09(2020)029 .
To keep this documentation simple, we will not actually read in a .lhe-file, and therefore can not go hands-on here. If you want to play around with this feature, we suggest you move to the example code.
To instantiate an Event
by an event=pylhe.LHEEvent
, we have to use evLHE=ngl.Event(eventFromFile="pylhe.LHEEvent")
.
It sets up the Event
with the default case of color-sorting the incoming and outgoing particles. To set up the showering of the incoming-intermediate particle dipoles, we have to use evLHE=ngl.Event(eventFromFile=pylhe.LHEEvent, productionDipoles='intermediate')
,
and if we not only want to have the production dipoles showered, but also the dipoles associated to the decay, we have to use evLHE=ngl.Event(eventFromFile=pylhe.LHEEvent,productionDipoles='intermediate',decayDipoles=True)
.
Note that you will probably seldom use these additional features and most oftenly only have to use evLHE=ngl.Event(eventFromFile="pylhe.LHEEvent")
, except if you work with top quarks.
Note that if you instantiate an Event
using a pylhe.LHEEvent
, you have access to the weight as well as an array of the FourVector
of each kind of particles in the form of an attribute. You can access the weight of the event via
ev=ngl.Event(eventFromFile="pylhe.LHEEvent")
ev.weight
If you want the fourvectors of all incoming up-type quarks and antiquarks, you can access them via
ev=ngl.Event(eventFromFile="pylhe.LHEEvent")
ev.incomingUp
.
In the same way, you can access ev.statusType
, with status
being incoming
,intermediate
or outgoing
, and Type
being Down
, Up
, Strange
, Charm
, Bottom
, Top
, Electron
, ENeutrino
, Muon
, MNeutrino
, Tau
, TNeutrino
, Gluon
, Photon
, ZBoson
, WBoson
or Higgs
.
Just as the Event
class, OutsideRegion
is very specific to the observable you are considering. Its whole purpose is to tell the Shower
, whether a FourVector
is pointing into the region where it gets vetoed. The nomenclature of outside comes from the textbook example of the interjet energy flow, where radiation that is not inside the jets gets vetoed.
We can initiate the OutsideRegion
with or without an Event
. Whether you should feed it an Event
or not depends on whether the region where you want to veto radiation depends on the distribution of the outgoing particles. We have one example code each for the usage with and without an Event
.
An instance of OutsideRegion
doesn't do anything. It containes the stub of a method outside(self,v)
which needs to be implemented by you. To do so, you need to write a method
def _outside(self,v):
#
# Code that checks whether v is a FourVector
# landing outside.
#
retVal=True # if v outside
retVal=False # if v not outside
return (retVal)
and - after creating an instance of OutsideRegion
outsideRegion=ngl.OutsideRegion()
exchange the stub of outside(self,v)
of your instance outsideRegion
to your method by invoking
outsideRegion.outside = _outside.__get__(outsideRegion,ngl.OutsideRegion)
For more details we refer to the two example codes which show the handling of the OutsideRegion
class.
Let us now get to the core of our resummation precedure, the showering of an Event
. To instantiate a Shower
, we feed it the following parameters (most of which come with a default choice):
# event: Event,
# outsideRegion: OutsideRegion,
# nsh: int=50,
# nbins: int=100,
# tmax: float=0.1,
# cut: float=5.0,
# fixedOrderExpansion: bool=True,
# virtualSubtracted: bool=False
Of these parameters, event
and outsideRegion
unsurprisingly contain the Event
to shower with the respective OutsideRegion
under consideration. The number of showerings you want to apply on the event
is fed in via nsh
(50
by default). To create the Hist
which will eventually be the result of the resummation, we can change the nbins
(100
by default) and tmax
(0.1
by default). We can also change the collinear cutoff cut
(5.0
by default), which corresponds to $\eta_{max}$ as discussed in (A.14) of doi:10.1007/JHEP09(2020)029. Finally, we can decide whether or not we want to calculate the first two expansion parameters of the resummation as given in (4.3) of doi:10.1007/JHEP09(2020)029 in fixedOrderExpansion
(True
by default). The last option is virtualSubtracted
will most likely have to be turned off (as is the default), it is used to subtract the global one-loop part from the soft anomalous dimension as discussed in (3.5) of doi:10.1007/JHEP04(2019)020.
Instead of going through the details of the showering precedure we refer to Appendix A of doi:10.1007/JHEP09(2020)029, where this is explained in a very detailed fashion.