#! /usr/bin/python3 ############################################################# # # Einfuehrung in die Theoretische Teilchenphysik (ETTP/TTP0) # Karlsruher Institut fuer Technologie, WS 2017/18 # # Sheet 12, Exercise 1, Higgs-Produktion via W Fusion # # Adapted by: Cody B Duncan # Original author: Michael Rauch # ################################################################### # We begin with the not-so interesting stuff: setting up the machinery # to allow us to visualize the output, and defining constants, etc. # import the libraries needed for the problem import sys import random import numpy as np import matplotlib.pyplot as plt from matplotlib.backends.backend_pdf import PdfPages # Input values # TODO: set the number of integration points, centre of mass energy # and the number of bins in the plots # number of integration points # In Monte Carlo algorithms (of which this is one), the error in your # distribution decreases with the number of points sampled. The tradeoff # being that you have to sample a large number of points numpoints = ... # Choose a centre of mass energy (e.g. ILC energy = 500 GeV) # Note, we'll be working in units of GeV throughout the exercise sqrt_s = ... # Number of bins for the plot nbins = ... # Higgs mass mh = 125. # Fermi constant for the weak coupling GF = 1.16638e-5 # W boson mass mw = 80.385 # Plot range information # Note that these numbers are set for ILC sqrt_s # so you may need to change some of them for larger or smaller values of # sqrt_s xbinrange = { 'x3': (0,1), # x3 and x4 are fractions, so do NOT change range 'x4': (0,1), 'pt_nu': (0,500), # pt's can increase with sqrt_s so make sure to 'pt_nubar': (0,250), # change if needed. 'pt_H': (0,500), 'y_nu': (-5,5), # Rapidity limits (usually determined by detector coverage) 'y_nubar': (-5,5), # So no real point in changing 'y_H': (-2.5,2.5), 'm_nunubar': (0,1000), # Invariant mass of the nu, nubar pair, so will scale with centre of mass energy 'deltay_nunubar': (-10,10), # Likewise, this is the difference in rapidities 'deltaphi_nunubar': (-np.pi,np.pi) # Differences in azimuthal angle } # Arrays for storing future data bins = {} step = {} distributions = {} val = {} # We're going to create a plot for each of the axes defined above for key in xbinrange.keys(): # Get the max and min values of the range low, up = xbinrange[key] # Generate the bins and step sizes bins[key],step[key] = np.linspace(low,up,nbins+1,retstep=True) # Initialize the distributions to 0 distributions[key] = np.zeros(nbins+2) #################################################################### # Now we start setting up some physics input and defining the boundaries # of the problem we are considering. Here we build all the tools needed to # sample phase space points randomly # Derived quantities # TODO: fill in the derived quantities # Vacuum expectation value (vev) of the Higgs field # Instead of inserting numbers, use the ones we already have above # How is it related to the Fermi constant? v = ... # Weak coupling constant # How is it related to Fermi's constant and the W boson mass? g = ... # Kinematic upper bound on x_3 # Remember that it is the fraction of energy available (i.e. sqrt_s) after # producing the Higgs boson. As it's an upper bound, what will the theoretical # momenta of the other two (p4 and pH) be? onemsqs = ... # There is one degree of freedom that the matrix element is independent of: # the global azimuthal angle phi_3 around the beam (i.e. z) axis ## Thus the integral has been ## taken and the corresponding factor 2 pi put into the prefactor of the ## phase-space integral. This only randomizes the momenta, so we do not have a ## preferred direction in the x-y-plane. # TODO: define a function that does an azimuthal rotation of a momentum p, # remembering that the 0th component is time, and 1,2,3 # are x, y, z respectively. def axrotate(p, phi): # Do something and return the new momentum pRot return ... # We'll find we need a Lorentzian dot product in places, so let's # define it in a function # TODO: define the 4D dot product, remembering that the 0th component is # time, and 1,2,3 are x, y, z respectively. def dot(p1, p2): # Do something then return the result return ... # Define an array to store the initial momenta p1, p2, # and the final momenta p3, p4, p5. The initial are fixed, since they # are the colliding particles, and the final will be filled by # sampling phase-space points # TODO: fill in the initial, incoming momenta. I have put placeholders # but you should put in values derived from above # the float keyword argument is just to tell numpy.array that the # values being stored are floats, ie non-integer numbers p1 = np.array([E1, px1, py1, pz1], float) p2 = np.array([E2, px2, py2, pz2], float) # Final state: p3 = np.array([0, 0, 0, 0], float) p4 = np.array([0, 0, 0, 0], float) p5 = np.array([0, 0, 0, 0], float) pArray = [p1, p2, p3, p4, p5] # Initialize the important quantities: crosssection = 0 error = 0 ########################################################### # Now here is the core of this exercise. We will sample phase-space # randomly and add it to our distributions. # Loop over the phase-space points for point in range(numpoints): # A very simple, though useful progress bar if (point + 1)%10000 == 0: sys.stderr.write(f"\rEvent {point+1} / {numpoints}") # Every phase-space point we sample will have a weight associated to it # i.e. how probable is such an event # Global factor from the phase-space integral # TODO: What is the global factor? # Hint: it is related only to pi, s and some numbers weight = ... # Convert the randomly sample numbers into phase-space integration variables # We'll need to keep track of Jacobians from transforming the integrals # to integrals on the unit hypercube # TODO: uniformly sample x3. What's its upper bound? x3 = ... # Modify the weight by this upper bound weight *= onemsqs # Now sample x4, but note that it will be constrained by the energy left # TODO: Define x4. It will have some flat sampling of the energy left by # not using up the full upper bound of x3 # Hint: If you are stuck, look below at the weight factor and see # if you can justify why this is the case x4 = ... # Modify the weight by the factor that was randomly sampled weight *= 1-(mh/sqrts)**2/(1.-x3) - (onemsqs-x3) # Generate the costheta of p3 (again uniformly sampled.) # TODO: Sample costheta3 costheta3 = ... # Updated weight weight *= 2 # Generate the azimuthal angle of p4 (uniformly sampled) # TODO: Sample phi4 phi4 = ... # Updated weight weight *= 2.*np.pi # Shorthands for derived quantities # TODO: Fill in the missing derived quantities # Simple sintheta relation sintheta3 = ... # How is costheta4 related the quantities onemsqs, x3, and n4? costheta4 = ... # Simple sintheta relation sintheta4 = ... # Now we will generate the momenta in explicit detail # TODO: Fill in the missing values below #### # Momentum 3 (remember arrays are indexed 0, so the first index # will be one less, while the second index refers to the momentum's # own entries (E, px, py, pz)) # Energy of p3: pArray[2][0] = ... # x, y, z components in order pArray[2][1] = ... pArray[2][2] = ... pArray[2][3] = ... #### # Momentum 4 # E4: pArray[3][0] = ... # 3-Momentum pArray[3][1] = ... pArray[3][2] = ... pArray[3][3] = ... # TODO: Perform the global azimuthal angle rotation by first sampling # the global angle phi3 phi3 = ... pArray[2] = azrotate(pArray[2], phi3) pArray[3] = azrotate(pArray[3], phi3) # TODO: Calculate p5 by energy and momentum conservation pArray[4] = ... # TODO: Calculate the q vectors, i.e. the momenta of the exchanged W bosons # What momenta are the related to? Make sure to get the sign of them # right q1 = ... q2 = ... # Matrix Element time, as calculated in the exercise # TODO: insert the final matrix element from sheet 12, part 4 me = ... # Modify the weight of the sampled point: weight *= me # Include the flux factor of the incoming beams weight /= (4.*dot(pArray[1], pArray[1])) # Include the number of points we are considering weight /= numpoints # TODO: multiply the weight by the conversion factor from GeV^2 # (natural units) to femtobarns weight *= ... # Update the cross section with the weight and update the error crosssection += weight error += weight**2 ######### End of the physics ################# # Store values for the histograms # TODO: Fill in missing values val['x3'] = ... val['x4'] = ... # How is pT calculated? val['pt_nu'] = ... val['pt_nubar'] = ... val['pt_H'] = ... # How is rapidity defined? Hint, you may need np.log val['y_nu'] = ... val['y_nubar'] = ... val['y_H'] = ... # How is the invariant mass of a pair of particles defined val['m_nunubar'] = ... # Difference in the rapidities val['deltay_nunubar'] = val['y_nu']-val['y_nubar'] # Difference in the azimuthal angles. # How can you calculate the azimuthal only from the momenta # Hint, you may find the function np.arctan2 useful. val['deltaphi_nunubar'] = ... # Here we just do some modulus arithmetic to ensure we have values between # pi and -pi if val['deltaphi_nunubar'] > np.pi: val['deltaphi_nunubar'] = val['deltaphi_nunubar'] - 2*np.pi elif val['deltaphi_nunubar'] < -np.pi: val['deltaphi_nunubar'] = val['deltaphi_nunubar'] + 2*np.pi # Update the distributions with the corresponding weight for this # sampled phase-space point. for key in distributions.keys(): distributions[key][np.digitize(val[key],bins[key])] += weight ############################################################################ # The final part of our program, which calculates errors and finally # produces the plots we have been attempting to make. # error = sqrt( (-^2)/N ) error = np.sqrt(error-crosssection**2/numpoints) sys.stderr.write("\n") print("Cross section: {} +/- {} fb".format(crosssection,error)) pdfout = PdfPages('WFusionHiggs.pdf') # Loop over all the distributions we want to plot for key in distributions.keys(): plt.hist(bins[key][0:nbins]+0.5*step[key], bins=nbins, range=xbinrange[key], weights=distributions[key][1:nbins+1]) plt.xlabel(key) # save to pdf pdfout.savefig() # show on screen # plt.show() plt.close() pdfout.close() ########################################################################### # If you made it to the end, a hearty congratulations! I hope you enjoyed this # exercise, and that it worked. If either of those were untrue, please feel # free to drop me an email: cody.b.duncan@kit.edu. # It has been a pleasure to teach you all ETTP this semester. Thank you for # all your hard work, and good luck with any future physics courses. #