forked from key4hep/k4geo
/
lcio_particle_gun.py
147 lines (102 loc) · 3.45 KB
/
lcio_particle_gun.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
#####################################
#
# simple script to create lcio files with single particle
# events - modify as needed
# @author F.Gaede, DESY
# @date 1/07/2014
#
# initialize environment:
# export PYTHONPATH=${LCIO}/src/python:${ROOTSYS}/lib
#
#####################################
import math
import random
from array import array
# --- LCIO dependencies ---
from pyLCIO import UTIL, EVENT, IMPL, IO, IOIMPL
#---- number of events per momentum bin -----
nevt = 1000
outfile = "mcparticles.slcio"
#--------------------------------------------
wrt = IOIMPL.LCFactory.getInstance().createLCWriter( )
wrt.open( outfile , EVENT.LCIO.WRITE_NEW )
print( " opened outfile: " , outfile )
random.seed()
#========== particle properties ===================
#momenta = [ 1. , 3., 5., 10., 15., 25., 50., 100. ]
momenta = [ 5. ]
genstat = 1
pdg = -13
charge = +1.
#pdg = 211
mass = 0.105658
theta = 85./180. * math.pi
#theta = 20./180. * math.pi
decayLen = 1.e32
#=================================================
# write a RunHeader
run = IMPL.LCRunHeaderImpl()
run.setRunNumber( 0 )
run.parameters().setValue("Generator","${lcgeo}_DIR/examples/lcio_particle_gun.py")
run.parameters().setValue("PDG", pdg )
run.parameters().setValue("Charge", charge )
run.parameters().setValue("Mass", mass )
wrt.writeRunHeader( run )
#================================================
for p in momenta:
for j in range( 0, nevt ):
col = IMPL.LCCollectionVec( EVENT.LCIO.MCPARTICLE )
evt = IMPL.LCEventImpl()
evt.setEventNumber( j )
evt.addCollection( col , "MCParticle" )
phi = random.random() * math.pi * 2.
energy = math.sqrt( mass*mass + p * p )
px = p * math.cos( phi ) * math.sin( theta )
py = p * math.sin( phi ) * math.sin( theta )
pz = p * math.cos( theta )
momentum = array('f',[ px, py, pz ] )
epx = decayLen * math.cos( phi ) * math.sin( theta )
epy = decayLen * math.sin( phi ) * math.sin( theta )
epz = decayLen * math.cos( theta )
endpoint = array('d',[ epx, epy, epz ] )
#--------------- create MCParticle -------------------
mcp = IMPL.MCParticleImpl()
mcp.setGeneratorStatus( genstat )
mcp.setMass( mass )
mcp.setPDG( pdg )
mcp.setMomentum( momentum )
mcp.setCharge( charge )
if( decayLen < 1.e9 ) : # arbitrary ...
mcp.setEndpoint( endpoint )
#-------------------------------------------------------
#-------------------------------------------------------
col.addElement( mcp )
wrt.writeEvent( evt )
wrt.close()
#
# longer format: - use ".hepevt"
#
#
# int ISTHEP; // status code
# int IDHEP; // PDG code
# int JMOHEP1; // first mother
# int JMOHEP2; // last mother
# int JDAHEP1; // first daughter
# int JDAHEP2; // last daughter
# double PHEP1; // px in GeV/c
# double PHEP2; // py in GeV/c
# double PHEP3; // pz in GeV/c
# double PHEP4; // energy in GeV
# double PHEP5; // mass in GeV/c**2
# double VHEP1; // x vertex position in mm
# double VHEP2; // y vertex position in mm
# double VHEP3; // z vertex position in mm
# double VHEP4; // production time in mm/c
#
# inputFile >> ISTHEP >> IDHEP
# >> JMOHEP1 >> JMOHEP2
# >> JDAHEP1 >> JDAHEP2
# >> PHEP1 >> PHEP2 >> PHEP3
# >> PHEP4 >> PHEP5
# >> VHEP1 >> VHEP2 >> VHEP3
# >> VHEP4;