#! /usr/bin/env python """ Monte Carlo decays @author: Jorgen Midtbo """ import numpy as np import matplotlib.pyplot as plt from math import pi np.random.seed(2) # set seed for reproducibility # ====== define useful functions ======== def proj(v,u): # projects v onto u if np.linalg.norm(u) > 0: return np.dot(u,np.transpose(v))/float(np.dot(u,u.T)) * u else: return u # Contraction of two four-vectors a and b (inner product in Minkowski space) def minkowskidot(a,b): return float(a[0,0]*b[0,0]-a[0,1]*b[0,1]-a[0,2]*b[0,2]-a[0,3]*b[0,3]) # Squared four vector def minkowskinorm(a): return minkowskidot(a,a) # Finds four-momenta of particles 2&3 going back-to-back from # decay of particle 1 in the frame where particle 1 has four-momentum P1 # # particle 1 = decaying particle # particle 2 & particle 3 = decay products # primed system is rest frame of particle 1, unprimed is lab frame # rotated system is at rest in lab system, # but rotated so particle one goes in +x direction def decayfun(m1,P1,m2,m3): p1 = P1[0,1:4] p1abs = np.sqrt( float( np.dot( p1 , np.transpose(p1) ) ) ) # 3-momentum # of particle 1 in # lab frame # == Kinematical decay in RF of particle 1 == p2absprime = np.sqrt( (m1**2-m2**2-m3**2)**2 - (2*m2*m3)**2 ) / (2*m1) # abs-val of 3-momentum of particle 2/3 in RF of particle 1 # Pick a point uniformly on a sphere U, V = np.random.uniform(0,1,2) phi = 2*pi*U theta = np.arccos(2*V-1) # Calculate cartesian 3- and 4-momentum of particle 2&3 E2prime = np.sqrt( p2absprime**2 + m2**2 ) E3prime = np.sqrt( p2absprime**2 + m3**2 ) p2prime = np.matrix([ p2absprime*np.sin(theta)*np.cos(phi) , p2absprime*np.sin(theta)*np.sin(phi) , p2absprime*np.cos(theta) ]) P2prime = np.matrix([ E2prime , p2prime[0,0] , p2prime[0,1] , p2prime[0,2] ]) P3prime = np.matrix([ E3prime , -p2prime[0,0] , -p2prime[0,1] , -p2prime[0,2] ]) # == Back-transform to lab frame == # First check whether it is necessary to boost if p1abs > 1e-10: # Lorentz boost along x-direction to get to rotated lab frame # (lab frame moves in negative x direction) vlab = -p1abs/np.sqrt(p1abs**2 + m1**2) # velocity of particle 1 in lab frame gamma = 1/np.sqrt(1-vlab**2) P2rot = np.matrix([ gamma*(P2prime[0,0] - vlab*P2prime[0,1]) , gamma*(P2prime[0,1] - vlab*P2prime[0,0]) , P2prime[0,2] , P2prime[0,3] ]) P3rot = np.matrix([ gamma*(P3prime[0,0] - vlab*P3prime[0,1]) , gamma*(P3prime[0,1] - vlab*P3prime[0,0]) , P3prime[0,2] , P3prime[0,3] ]) # == Rotate back to lab frame == # Calculate the unit vectors of the rotated system axes in terms of lab axes # The definition is that x axis is along p1. # For the other axes we must make a choice - y&z directions are undetermined, # only the yz plane is determined from x choice. But since we have drawn # random angles and the yz plane is not boosted, the choice does not matter # as long as we are consistent from event to event. # So we pick two vectors orthogonal to p1 and do Gram-Schmidt orthogonalization: v1 = p1 v2 = np.matrix([ p1[0,1] , -p1[0,0] , 0 ]) v3 = np.matrix([ p1[0,2] , 0 , -p1[0,0] ]) u1 = v1 u2 = v2 - proj(v2,u1) u3 = v3 - proj(v3,u1) - proj(v3,u2) xrot = u1/np.linalg.norm(u1) if np.linalg.norm(u1) > 0 else np.matrix([0,0,1]) yrot = u2/np.linalg.norm(u2) if np.linalg.norm(u2) > 0 else np.matrix([0,1,0]) zrot = u3/np.linalg.norm(u3) if np.linalg.norm(u3) > 0 else np.matrix([1,0,0]) # Form a matrix T which takes a vector in the lab basis to a vector # in the rotated basis by T = np.concatenate( (xrot , yrot , zrot) , axis=0 ) # What we need is to rotate from rotated basis to lab basis, so we need the inverse # - which is the transpose, since rotation matrices are orthogonal. # Also, to ease calculation, we let T be the 3x3 submatrix of T4, setting the [0,0] #component of T4 to 1 to leave time component invariant under this spatial rotation T4 = np.matrix([[1, 0, 0, 0], [0,T[0,0],T[0,1],T[0,2]], [0,T[1,0],T[1,1],T[1,2]], [0,T[2,0],T[2,1],T[2,2]] ]) P2 = T4.T*P2rot.T P3 = T4.T*P3rot.T P2 = P2.T P3 = P3.T # If it was unneccessary, i.e. decay happened in lab frame, then else: P2 = P2prime P3 = P3prime # Finished! return P2, P3 # ====== Main program ====== # Define masses mE = 600 # GeV/c^2 mD = 500 # GeV/c^2 mC = 200 # GeV/c^2 mB = 150 # GeV/c^2 mA = 100 # GeV/c^2 md = 0 mc = 0 mb = 1.8 # GeV/c^2 ma = 1.8 # GeV/c^2 # Allocate lists to store invariant masses mab_list = [] mac_list = [] mad_list = [] mbc_list = [] mbd_list = [] mcd_list = [] N = 10000 # How much loop? for i in range(N): # Decay E from rest through chain pE = np.matrix([mE,0,0,0]) pd, pD = decayfun(mE,pE,md,mD) pc, pC = decayfun(mD,pD,mc,mC) pb, pB = decayfun(mC,pC,mb,mB) pa, pA = decayfun(mB,pB,ma,mA) # Make all invariant mass combinations mab_list.append( np.sqrt(minkowskinorm(pa+pb)) ) mac_list.append( np.sqrt(minkowskinorm(pa+pc)) ) mad_list.append( np.sqrt(minkowskinorm(pa+pd)) ) mbc_list.append( np.sqrt(minkowskinorm(pb+pc)) ) mbd_list.append( np.sqrt(minkowskinorm(pb+pd)) ) mcd_list.append( np.sqrt(minkowskinorm(pc+pd)) ) # Plot results as normalized histograms bins = np.linspace(0, 350, 151) # Histogram bins mab_hist = np.histogram(mab_list, bins=bins, normed=True) mac_hist = np.histogram(mac_list, bins=bins, normed=True) mad_hist = np.histogram(mad_list, bins=bins, normed=True) mbc_hist = np.histogram(mbc_list, bins=bins, normed=True) mbd_hist = np.histogram(mbd_list, bins=bins, normed=True) mcd_hist = np.histogram(mcd_list, bins=bins, normed=True) f, ax = plt.subplots(1,1) ax.step(mab_hist[1][:-1], mab_hist[0], where="post", label="ab") ax.step(mac_hist[1][:-1], mac_hist[0], where="post", label="ac") ax.step(mad_hist[1][:-1], mad_hist[0], where="post", label="ad") ax.step(mbc_hist[1][:-1], mbc_hist[0], where="post", label="bc") ax.step(mbd_hist[1][:-1], mbd_hist[0], where="post", label="bd") ax.step(mcd_hist[1][:-1], mcd_hist[0], where="post", label="cd") ax.set_xlabel("$m_{ij}\,\mathrm{[GeV/c^2]}$") ax.set_ylabel("Normalized to one") ax.legend() plt.show()