You are a guest. Restricted access. Read more.
SCaVis manual

Particle Physics

Showing data from HepData database

The best approach is to download “scavis” Jython script from Durham HepData database and run it in the SCaVis editor or as a batch script. This is the best supported approach which also allows to deal with statistical and systematical errors separately since data are imported to the P1D object. You can further modify P1D graphical attributes or replace HPlot canvas with other canvaces (like HPlotJa). You can also combine statistical and systematic errors using the methods of the P1D class.

Here is a step-by step instruction:

You will see an image with the HepData points.

Reading ROOT files

You can read ROOT files with histograms or TGraph as discussed in reading root files. SCaVis use FreeHep Java API to read such files and convert to SCaVis native representation.

Reading events with Monte Carlo simulation

SCaVis can be used to read and process Monte Carlo simulation events with particle collisions. The files can be download Monte Carlo events can be downloaded from HepSim Monte Carlo repository supported at the ANL laboratory. Read the HepSim manual that shows how to read Monte Carlo simulations (including NLO) data. The HepSim also includes scripts illustrating how to access 4-momentum information on particles.

You can create a list with particles from HepSim and ProMC files using hephysics.hepsim.PromcUtil class. You can create a lists for fast manipulation with hephysics.jet.FastParticle class, or list with the standard hephysics.particle.LParticle class.

Reading PDG table

Use SCaVis to read PDG table. Here is a handy script that reads the file “mass_width_2008.csv” located inside the SCaVis directory “macros/data”.

Code example

  Download for this example is disabled for non-members
 1: # Physics/HEP | P | 1.7 |  Matt Williams  | Reading information on partions from PDG data (2008)
 2:
 3: # -*- coding: utf-8 -*-
 4:
 5: """
 6: Particle data from the PDG
 7: """
 8:
 9: class Particle(object):
10:     """
11:     An elementary or composite particle
12:     """
13:     def __init__(self, mass, masspositiveerror, massnegativeerror, width, widthpositiveerror, widthnegativeerror, isospin, gparity, totalspin, spaceparity, chargeconjugation, antiparticleflag, number, charge, rank, particlestatus, name, quarkcontent):
14:         self.number = number
15:         self.name = name
16:         self.mass = mass
17:         self.masspositiveerror = masspositiveerror
18:         self.massnegativeerror = massnegativeerror
19:         self.width = width
20:         self.widthpositiveerror = widthpositiveerror
21:         self.widthnegativeerror = widthnegativeerror
22:         self.isospin = isospin
23:         self.gparity = gparity
24:         self.totalspin = totalspin
25:         self.spaceparity = spaceparity
26:         self.chargeconjugation = chargeconjugation
27:         self.antiparticleflag = antiparticleflag
28:         self.charge = charge
29:         self.rank = rank
30:         self.particlestatus = particlestatus
31:         self.quarkcontent = quarkcontent
32:     def add_decay(self, decay):
33:         pass
34:     def __str__(self):
35:         return self.name+' : '+self.quarkcontent
36:     def __repr__(self):
37:         return '<Particle '+self.name+'>'
38:
39: class Decay(object):
40:     def __init__(self):
41:         pass
42:
43: file=SystemDir+fSep+"macros"+fSep+"examples"+fSep+"data"+fSep+"mass_width_2008.csv"
44: pdginfofile = open("mass_width_2008.csv")
45: particles = []
46: for line in pdginfofile:
47:     if line[0] == '*':
48:         continue
49:     line = line.split(',')
50:     if len(line) != 18:
51:         continue
52:     line = [element.strip() for element in line]
53:     print line
54:     particle = Particle(line[0], line[1], line[2], line[3], line[4], line[5], line[6], line[7], line[8], line[9], line[10], line[11], line[12], line[13], line[14], line[15], line[16], line[17])
55:     particles.append(particle)
56:
57: #for particle in particles: print particle
58: print particles

This examples create “particle” objects with all needed informations (masses, spins, charges etc). The script assumes that you run it inside the SCaVis which defines the path to the mass_width_2008.csv file. Look at examples in the Scientific data analysis using Jython and Java

Lorentz representation of a particle

You can represent a relativistic particle using a Lorentz representation. Use the class hephysics.particle.LParticle which defines a particle in terms of the four-Lorentz vector. A particle is characterized by either (px,py,pz,E) or (x,y,z,time).

For example, below we define a “top” quark with the mass 170 GeV and with the momentum pZ=400 GeV

from hephysics.particle import *
Top=LParticle("top quark",170)
Top.setV3( 0,0,400 )
print Top.toString()

Now we are ready to go with a small Monte Carlo simulation for decay of the top particle. Below we calculate the maximum opening angle between its decay products

Code example

  Download for this example is disabled for non-members
 1: # Physics/HEP | C | 1.7 | S.Chekanov | Feynman diagram for double Higgs production
 2:
 3: from java.awt import Color
 4: from jhplot  import *
 5: from java.util import Random
 6: from hephysics.particle import *
 7: from math import *
 8:
 9: P=400
10: # define initial  Lorentz particle
11: Top=LParticle("Top quark",170)
12: Top.setV3( 0,0,P )
13: print Top.toString()
14:
15: # define Lorentz particles
16: W = LParticle("W",80)
17: b = LParticle("b",5)
18: q1 = LParticle("q1",0.004)
19: q2 = LParticle("q2",0.004)
20:
21: r = Random()
22: PI=3.1415926
23: PI2=2*PI
24: deg= 180 / PI
25:
26: h1 = H1D("maxAngle",100, 0,180)
27: for i in range(10000):
28:      if (i%1000==0): print "event=",i
29:
30:      Top.twoBodyDecay(W, b,1);
31:      theta = acos(2.0*r.nextDouble()-1.0 )
32:      phi   = PI2 *r.nextDouble()
33:      Top.setThetaPhiP(theta,phi,P)
34:      W.boost( Top )
35:      b.boost( Top )
36:      # W decay
37:      W.twoBodyDecay(q1,q2,1)
38:      # the boost
39:      q1.boost( W )
40:      q2.boost( W )
41:
42:      # print q1.toString()
43:      p=P0D("max angle")
44:      p.add(W.angle(q1)*deg)
45:      p.add(W.angle(q2)*deg)
46:      pmax=p.getMax()
47:      h1.fill(pmax)
48:
49: c1 = HPlot("Canvas",600,400,1, 1)
50: c1.setGTitle("Max angle t &rarr; W b &rarr; b q #bar{q}")
51: c1.visible(1)
52: c1.setRange(0,180,0,1000)
53: c1.setNameX("angle_{max} [deg]")
54: c1.setNameY("Events")
55:
56: h1.setFill(1)
57: h1.setFillColor(Color.blue)
58: c1.draw(h1)
59:
60: # export to EPS figure
61: c1.export("angle"+".eps")

The output of the code is below:

Feynman diagrams

One can draw Feynman diagrams programically using Java, Python (Jython) or any scripting language supported by SCaVis. Look at the section Drawing diagrams which shows how to draw diagrams and export to a number of image formats.

SCaVis includes Redberry computer algebra (CA) system designed for algebraic manipulations with tensors. It supports symbolic algebra, tensor-specific transformations, simplification, Feynman diagrams and one-loop counterterms Please read “Introduction to Redberry: the computer algebra system designed for tensor manipulation” by D. A. Bolotin and S. V. Poslavsky (arXiv:1302.1219 [cs.SC]).

Currently, the description of Redberry is beyond the scope of this manual. You can run the Groovy tutorials that are given in the Redberry link.

Jet algorithms

SCaVis includes a number of jet algorithms implemented in Java. This is a list of the algorithms:

  • hephysics.jet.KTjet - typical jet clustering algorithms (kT, anti-kT, Cambridge/Aachen) used for pp collision using pseudorapidity for distance measure (fast).
  • hephysics.jet.SCJet - jet clustering algorithms (kT, anti-kT, Cambridge/Aachen) used for pp collision using rapidity for distance measure.
  • e+e- jet algorithms including Jade, Geneva and Durham

The hephysics.jet.KTjet class is interfaces with the LParticle class. You can run the longitudinally invariant kT algorithm in inclusive mode. Jet can be clustered using Cambridge/Aachen or anti-kT approaches. It uses Eta-Phi space to make jets (not Rapidity-Phi). Look at the KT jet package summary.

Let us give an example of how to visualize anti-kT jets in Eta-Phi plane using Jython script:

Code example

  Download for this example is disabled for non-members
 1: # Physics/HEP | C | 2.2 | S.Chekanov | Running the anti-Kt jet algorithm for high energy particle collisions
 2:
 3: from java.awt import Color
 4: from java.util import ArrayList
 5: from jhplot  import *
 6: from hephysics.jet import *
 7: from hephysics.particle import *
 8: from jhplot.shapes import *
 9:
10: plist=ArrayList()
11: # add px,py,pz,e
12: p=ParticleD(-0.880,-1.233,-157.431,157.439); plist.add(p)
13: p=ParticleD(10.0, -2.3211,-15.0785, 18.250); plist.add(p)
14: p=ParticleD(552.60,-143.38,127.405, 584.95); plist.add(p)
15: p=ParticleD(-58.56, 13.672,-60.729, 85.479); plist.add(p)
16: p=ParticleD(249.62, 59.98,-252.05,  359.78); plist.add(p)
17: p=ParticleD(8.31, -1.26, -13.692,   16.072); plist.add(p)
18: p=ParticleD(-132.34, 31.19,-133.81, 190.7);  plist.add(p)
19: p=ParticleD(31.629, -7.989, 7.4822, 33.47);  plist.add(p)
20:
21: R=0.6 # jet radius. Use anti-kT jets
22: ktjet = SCJet(R, 1, -1, 200.0,True);
23: ktjet.buildJets(plist);
24: ktjet.printJets()
25: jets=ktjet.getJetsSorted()
26:
27: c1 = HPlot("Canvas",600,600)
28: c1.setGTitle("anti-kT jets")
29: c1.visible(1)
30: c1.setRange(-4,4,-7,7)
31: c1.setNameX("&eta;")
32: c1.setNameY("&phi; [rad]")
33:
34: pp=P1D("particles")
35:
36: for i in range(plist.size()):
37:       p=plist.get(i)
38:       eta=p.getRapidity()
39:       phi=p.getPhi()
40:       pp.add(eta,phi)
41:       # psize=int(p.getEt()/10.)
42:       # pp.setSymbolSize(2)
43:
44: # show jets
45: for k in range(jets.size()):
46:     jet=P1D("jet Nr "+str(k))
47:     jet.setColor(Color.red)
48:     jet.setSymbolSize(7)
49:     jj=jets.get(k)
50:     etaJ=jj.getRapidity()
51:     phiJ=jj.phi()
52:     jet.add(etaJ,phiJ)
53:     c1.draw(jet)
54:     cic= Circle(etaJ, phiJ, R)
55:     cic.setFill(0)
56:     cic.setColor(Color.red)
57:     cic.setTransparency(0.5)
58:     c1.add(cic)
59:
60: c1.draw(pp)

This example will show this canvas with particles and jets plotted on top of each other:

man/science/physics/hep.txt · Last modified: 2015/04/15 19:20 by admin
CC Attribution-Share Alike 3.0 Unported
Powered by PHP Driven by DokuWiki Recent changes RSS feed Valid CSS Valid XHTML 1.0 Valid HTML5