1 | """ |
---|
2 | Python module |
---|
3 | |
---|
4 | This module provides ROOT IO interface for MCScore data |
---|
5 | |
---|
6 | [C] MCScoreROOTIO |
---|
7 | [f] loop_tree(tfile, analyze_vertex): |
---|
8 | |
---|
9 | Q, 2006 |
---|
10 | """ |
---|
11 | from array import array |
---|
12 | import string |
---|
13 | from Geant4.hepunit import * |
---|
14 | import mcscore |
---|
15 | import ROOT |
---|
16 | |
---|
17 | # ================================================================== |
---|
18 | # public symbols |
---|
19 | # ================================================================== |
---|
20 | __all__ = [ 'MCScoreROOTIO', 'loop_tree' ] |
---|
21 | |
---|
22 | |
---|
23 | # ================================================================== |
---|
24 | # class definition |
---|
25 | # ================================================================== |
---|
26 | # ------------------------------------------------------------------ |
---|
27 | # MCScoreROOTIO |
---|
28 | # ------------------------------------------------------------------ |
---|
29 | class MCScoreROOTIO: |
---|
30 | "ROOT IO interface for MCScore" |
---|
31 | def __init__(self, bsize=250): |
---|
32 | self.maxparticle = bsize # buffer size for #particles/vertex |
---|
33 | |
---|
34 | # ------------------------------------------------------------------ |
---|
35 | def define_tree(self): |
---|
36 | "define ROOT tree" |
---|
37 | # defining tree header... |
---|
38 | self.vtree = ROOT.TTree('vertex', 'mc vertex') |
---|
39 | self.ptree = ROOT.TTree('particle', 'mc particle') |
---|
40 | |
---|
41 | # vertex... |
---|
42 | self.a_x = array('d', [0.]); self.vtree.Branch('x', self.a_x, 'x/d') |
---|
43 | self.a_y = array('d', [0.]); self.vtree.Branch('y', self.a_y, 'y/d') |
---|
44 | self.a_z = array('d', [0.]); self.vtree.Branch('z', self.a_z, 'z/d') |
---|
45 | |
---|
46 | #global a_np, a_namelist, a_Z, a_A, a_ke, a_px, a_py, a_pz |
---|
47 | self.a_np = array('i', [0]); self.ptree.Branch('np', self.a_np, 'np/i') |
---|
48 | |
---|
49 | self.a_namelist = array('c', self.maxparticle*10*['\0']) |
---|
50 | # 10 characters/particle |
---|
51 | self.ptree.Branch('namelist', self.a_namelist, 'namelist/C') |
---|
52 | |
---|
53 | self.a_Z = array('i', self.maxparticle*[0]) |
---|
54 | self.ptree.Branch('Z', self.a_Z, 'Z[np]/I') |
---|
55 | |
---|
56 | self.a_A = array('i', self.maxparticle*[0]) |
---|
57 | self.ptree.Branch('A', self.a_A, 'A[np]/i') |
---|
58 | |
---|
59 | self.a_ke = array('d', self.maxparticle*[0.]) |
---|
60 | self.ptree.Branch('kE', self.a_ke, 'kE[np]/d') |
---|
61 | |
---|
62 | self.a_px = array('d', self.maxparticle*[0.]) |
---|
63 | self.ptree.Branch('px', self.a_px, 'px[np]/d') |
---|
64 | |
---|
65 | self.a_py = array('d', self.maxparticle*[0.]) |
---|
66 | self.ptree.Branch('py', self.a_py, 'py[np]/d') |
---|
67 | |
---|
68 | self.a_pz = array('d', self.maxparticle*[0.]) |
---|
69 | self.ptree.Branch('pz', self.a_pz, 'pz[np]/d') |
---|
70 | |
---|
71 | # ------------------------------------------------------------------ |
---|
72 | def fill_tree(self, vertex): |
---|
73 | "fill vertex information to ROOT tree" |
---|
74 | # ------------------------------------------------------------------ |
---|
75 | def push_pname(i0, pname): # local function |
---|
76 | n = len(pname) |
---|
77 | for i in xrange(n): |
---|
78 | self.a_namelist[i0+i] = pname[i] |
---|
79 | self.a_namelist[i0+n] = ' ' |
---|
80 | |
---|
81 | self.a_x[0] = vertex.x |
---|
82 | self.a_y[0] = vertex.y |
---|
83 | self.a_z[0] = vertex.z |
---|
84 | self.vtree.Fill() |
---|
85 | |
---|
86 | if vertex.nparticle > self.maxparticle: |
---|
87 | raise """ |
---|
88 | *** buffer overflow in #particles/vertex. |
---|
89 | *** please increment buffersize in MCScoreROOTIO(bsize). |
---|
90 | """, self.maxparticle |
---|
91 | |
---|
92 | idx_namelist = 0 |
---|
93 | self.a_np[0] = vertex.nparticle |
---|
94 | for ip in range(vertex.nparticle): |
---|
95 | particle = vertex.particle_list[ip] |
---|
96 | push_pname(idx_namelist, particle.name) |
---|
97 | idx_namelist += (len(particle.name)+1) |
---|
98 | self.a_Z[ip] = particle.Z |
---|
99 | self.a_A[ip] = particle.A |
---|
100 | self.a_ke[ip] = particle.kineticE |
---|
101 | self.a_px[ip] = particle.px |
---|
102 | self.a_py[ip] = particle.py |
---|
103 | self.a_pz[ip] = particle.pz |
---|
104 | |
---|
105 | self.a_namelist[idx_namelist] = '\0' |
---|
106 | self.ptree.Fill() |
---|
107 | |
---|
108 | |
---|
109 | # ================================================================== |
---|
110 | # functions |
---|
111 | # ================================================================== |
---|
112 | def loop_tree(tfile, analyze_vertex): |
---|
113 | """ |
---|
114 | loop ROOT tree in a ROOT file. |
---|
115 | * analyze_vertex: user function : analyze_vertex(MCVertex) |
---|
116 | """ |
---|
117 | avtree = tfile.Get("vertex") |
---|
118 | aptree = tfile.Get("particle") |
---|
119 | |
---|
120 | # reading vertex... |
---|
121 | n_vertex = avtree.GetEntries() |
---|
122 | for ivtx in xrange(n_vertex): |
---|
123 | avtree.GetEntry(ivtx) |
---|
124 | aptree.GetEntry(ivtx) |
---|
125 | |
---|
126 | # vertex |
---|
127 | vertex = MCScore.MCVertex(avtree.x, avtree.y, avtree.z) |
---|
128 | |
---|
129 | # reading secondary particles... |
---|
130 | nsec = aptree.np |
---|
131 | namelist = aptree.namelist |
---|
132 | pname = namelist.split() |
---|
133 | for ip in xrange(nsec): |
---|
134 | particle = MCScore.MCParticle(pname[ip], aptree.Z[ip], aptree.A[ip], |
---|
135 | aptree.kE[ip], aptree.px[ip], |
---|
136 | aptree.py[ip], aptree.pz[ip]) |
---|
137 | vertex.append_particle(particle) |
---|
138 | |
---|
139 | analyze_vertex(vertex) |
---|
140 | |
---|
141 | return n_vertex |
---|
142 | |
---|
143 | |
---|
144 | # ================================================================== |
---|
145 | # test |
---|
146 | # ================================================================== |
---|
147 | def test_t2root(): |
---|
148 | g = ROOT.TFile("reaction.root", 'recreate') |
---|
149 | |
---|
150 | rootio = MCScoreROOTIO() |
---|
151 | rootio.define_tree() |
---|
152 | |
---|
153 | f = open("reaction.dat") |
---|
154 | f.seek(0) |
---|
155 | |
---|
156 | while(1): |
---|
157 | vertex = MCScore.read_next_vertex(f) |
---|
158 | if vertex == 0: |
---|
159 | break |
---|
160 | |
---|
161 | # filling ... |
---|
162 | rootio.fill_tree(vertex) |
---|
163 | |
---|
164 | del vertex |
---|
165 | |
---|
166 | print ">>> EOF" |
---|
167 | f.close() |
---|
168 | |
---|
169 | # closing ROOT file |
---|
170 | g.Write() |
---|
171 | g.Close() |
---|
172 | |
---|
173 | # ------------------------------------------------------------------ |
---|
174 | def test_loop(): |
---|
175 | def my_analysis(vertex): |
---|
176 | vertex.printout() |
---|
177 | |
---|
178 | f = ROOT.TFile("reaction.root", 'read') |
---|
179 | nv = loop_tree(f, my_analysis) |
---|
180 | print "*** # of vertex= ", nv |
---|
181 | |
---|
182 | |
---|
183 | # ================================================================== |
---|
184 | # main |
---|
185 | # ================================================================== |
---|
186 | if __name__ == "__main__": |
---|
187 | test_t2root() |
---|
188 | test_loop() |
---|
189 | |
---|