plot.py
Go to the documentation of this file.
1 
2 from wirecell import units
3 import matplotlib.pyplot as plt
4 import matplotlib.patches as mpatches
5 import numpy
6 from collections import defaultdict
7 import math
8 
9 def plot_polyline(pts):
10  cmap = plt.get_cmap('seismic')
11  npts = len(pts)
12  colors = [cmap(i) for i in numpy.linspace(0, 1, npts)]
13  for ind, (p1, p2) in enumerate(zip(pts[:-1], pts[1:])):
14  x = numpy.asarray((p1.x, p2.x))
15  y = numpy.asarray((p1.y, p2.y))
16  plt.plot(x, y, linewidth=ind+1)
17 
18 
19 
20 def oneplane(store, iplane, segments=None):
21  '''
22  Plot one plane of wires.
23 
24  This plot is in protodune-numbers document.
25  '''
26  fig,axes = plt.subplots(nrows=1, ncols=3)
27 
28  uvw = "UVW"
29 
30  widths = [1, 2, 3]
31  wire_stride=20;
32 
33  iface = 0
34  face = store.faces[iface]
35 
36  cmap = plt.get_cmap('rainbow')
37 
38  for iplane in range(3):
39  ax = axes[iplane]
40 
41  planeid = face.planes[iplane]
42  plane = store.planes[planeid]
43 
44  wires = [w for w in plane.wires[::wire_stride]]
45 
46  nwires = len(wires)
47  colors = [cmap(i) for i in numpy.linspace(0, 1, nwires)]
48 
49  for wcount, wind in enumerate(wires):
50 
51  wire = store.wires[wind]
52  if segments and not wire.segment in segments:
53  continue
54 
55  color = colors[wcount]
56  if not iplane:
57  color = colors[nwires-wcount-1]
58 
59  p1 = store.points[wire.tail]
60  p2 = store.points[wire.head]
61  width = widths[wire.segment]
62  ax.plot((p1.z/units.meter, p2.z/units.meter), (p1.y/units.meter, p2.y/units.meter),
63  linewidth = width, color=color)
64 
65  ax.locator_params(axis='x', nbins=5)
66  ax.set_aspect('equal', 'box')
67  ax.set_xlabel("Z [meter]")
68  if not iplane:
69  ax.set_ylabel("Y [meter]")
70  ax.set_title("plane %d/%s" % (iplane,uvw[iplane]))
71 
72  return fig,ax
73 
74 def select_channels(store, pdffile, channels, labels=True):
75  '''
76  Plot wires for select channels.
77  '''
78  channels = set(channels)
79  bychan = defaultdict(list)
80 
81  # find selected wires and their wire-in-plane index
82  # fixme: there should be a better way!
83  for anode in store.anodes:
84  for iface in anode.faces:
85  face = store.faces[iface]
86  for iplane in face.planes:
87  plane = store.planes[iplane]
88  for wip,wind in enumerate(plane.wires):
89  wire = store.wires[wind]
90  if wire.channel in channels:
91  bychan[wire.channel].append((wire, wip))
92 
93  fig, ax = plt.subplots(nrows=1, ncols=1)
94 
95  for ch,wws in sorted(bychan.items()):
96  wws.sort(key=lambda ww: ww[0].segment)
97  for wire, wip in wws:
98  p1 = store.points[wire.tail]
99  p2 = store.points[wire.head]
100  width = wire.segment + 1
101  ax.plot((p1.z/units.meter, p2.z/units.meter),
102  (p1.y/units.meter, p2.y/units.meter), linewidth = width)
103  x = p2.z/units.meter
104  y = p2.y/units.meter
105  if x > 0:
106  hal="left"
107  else:
108  hal="right"
109  if labels:
110  t='w:%d ch:%d\nident:%d seg:%d' % \
111  (wip, wire.channel, wire.ident, wire.segment)
112  ax.text(x, y, t,
113  horizontalalignment=hal,
114  bbox=dict(facecolor='yellow', alpha=0.5, pad=10))
115  ax.set_xlabel("Z [meter]")
116  ax.set_ylabel("Y [meter]")
117  fig.savefig(pdffile)
118 
119 
120 
121 def allplanes(store, pdffile):
122  '''
123  Plot each plane of wires on a page of a PDF file.
124  '''
125  wire_step = 10 # how many wires to skip
126 
127  from matplotlib.backends.backend_pdf import PdfPages
128 
129  # make some global pltos
130  all_wire_x1 = list()
131  all_wire_x2 = list()
132  all_wire_z1 = list()
133  all_wire_z2 = list()
134  all_wire_anode = list()
135 
136  plane_colors=["blue","red","black"]
137 
138  with PdfPages(pdffile) as pdf:
139  # make channel plots
140  for anode in store.anodes:
141  edge_z = list()
142  edge_x = list()
143  edge_n = list()
144  edge_s = list()
145 
146  for iface in anode.faces:
147  face = store.faces[iface]
148 
149  for iplane in face.planes:
150  plane = store.planes[iplane]
151  wires_in_plane = [store.wires[wind] for wind in plane.wires]
152  wires = [w for w in wires_in_plane if w.segment == 0]
153  def pt(w): return store.points[w.head]
154  wires.sort(key=lambda a: pt(a).z)
155  def add_edge(w):
156  p = pt(w)
157  print (p,w.channel)
158  edge_z.append(p.z/units.m)
159  edge_x.append(p.x/units.m)
160  edge_n.append(w.channel)
161  edge_s.append('f%d p%d c%d wid%d' % (face.ident, plane.ident, w.channel, w.ident))
162  add_edge(wires[0])
163  add_edge(wires[-1])
164 
165  fig, ax = plt.subplots(nrows=1, ncols=1)
166  ax.scatter(edge_z, edge_x, s=1, c='red', marker='.')
167  for i,(z,x,s) in enumerate(zip(edge_z, edge_x, edge_s)):
168  hal = ["left","right"][i%2]
169  ax.text(z, x, s, horizontalalignment=hal,
170  bbox=dict(facecolor='yellow', alpha=0.5, pad=1))
171  for i in range(len(edge_n)//2):
172  z = 0.5*(edge_z[2*i]+edge_z[2*i+1])
173  n = 1+abs(edge_n[2*i] - edge_n[2*i+1])
174  x = edge_x[2*i]
175  ax.text(z, x, str(n), horizontalalignment='center',
176  bbox=dict(facecolor='yellow', alpha=0.5, pad=1))
177 
178 
179  ax.set_title("Edge Channels AnodeID: %d" % (anode.ident))
180  ax.set_xlabel("Z [meter]")
181  ax.set_ylabel("X [meter]")
182  pdf.savefig(fig)
183  plt.close()
184 
185  for anode in store.anodes:
186  seg_x1 = [list(),list(),list()]
187  seg_x2 = [list(),list(),list()]
188  seg_z1 = [list(),list(),list()]
189  seg_z2 = [list(),list(),list()]
190  seg_col = [list(),list(),list()]
191  for iface in anode.faces:
192  face = store.faces[iface]
193  for iplane in face.planes:
194  plane = store.planes[iplane]
195  for wind in plane.wires:
196  wire = store.wires[wind]
197  p1 = store.points[wire.tail]
198  p2 = store.points[wire.head]
199  seg = wire.segment
200  seg_x1[seg].append(p1.x/units.meter)
201  seg_x2[seg].append(p2.x/units.meter)
202  seg_z1[seg].append(p1.z/units.meter)
203  seg_z2[seg].append(p2.z/units.meter)
204  seg_col[seg].append(plane_colors[iplane%(len(plane_colors))])
205  continue # wires
206  continue # planes
207  continue # faces
208 
209  fig, axes = plt.subplots(nrows=3, ncols=1)
210  for seg in range(3):
211  ax = axes[seg]
212  ax.scatter(seg_z2[seg], seg_x2[seg], c=seg_col[seg], s=1, marker='.')
213  ax.set_title("AnodeID %d wires, seg %d, head (%d wires)" %
214  (anode.ident, seg, len(seg_col[seg])))
215  plt.tight_layout()
216  pdf.savefig(fig)
217  plt.close()
218 
219  fig, axes = plt.subplots(nrows=3, ncols=1)
220  for seg in range(3):
221  ax = axes[seg]
222  ax.scatter(seg_z1[seg], seg_x1[seg], c=seg_col[seg], s=1, marker='.')
223  ax.set_title("AnodeID %d wires, seg %d, tail (%d wires)" %
224  (anode.ident, seg, len(seg_col[seg])))
225  plt.tight_layout()
226  pdf.savefig(fig)
227  plt.close()
228  continue # anodes
229 
230  for anode in store.anodes:
231  wire_x1 = list()
232  wire_x2 = list()
233  wire_z1 = list()
234  wire_z2 = list()
235  wire_anode = list()
236 
237 
238  for iface in anode.faces:
239  face = store.faces[iface]
240 
241  first_wires = list()
242 
243  for iplane in face.planes:
244  plane = store.planes[iplane]
245 
246  print ("anodeID:%d faceID:%d planeID:%d" %
247  (anode.ident, face.ident, plane.ident))
248 
249  first_wires.append(plane.wires[:2])
250 
251  fig, ax = plt.subplots(nrows=1, ncols=1)
252  ax.set_aspect('equal','box')
253  for wind in plane.wires[::wire_step]:
254  wire = store.wires[wind]
255  p1 = store.points[wire.tail]
256  p2 = store.points[wire.head]
257  width = wire.segment + .1
258  ax.plot((p1.z/units.meter, p2.z/units.meter),
259  (p1.y/units.meter, p2.y/units.meter), linewidth = width)
260  wire_x1.append(p1.x/units.meter)
261  wire_z1.append(p1.z/units.meter)
262  wire_x2.append(p2.x/units.meter)
263  wire_z2.append(p2.z/units.meter)
264  wire_anode.append(anode.ident)
265 
266  wirex = None
267  for wcount, wind in enumerate([plane.wires[0], plane.wires[-1]]):
268  wire = store.wires[wind]
269  print ("\twcount:%d wind:%d wident:%d chan:%d" % (wcount,wind,wire.ident,wire.channel))
270  p1 = store.points[wire.tail]
271  p2 = store.points[wire.head]
272  x = p2.z/units.meter
273  y = p2.y/units.meter
274  wirex = p2.x/units.meter
275  hal="center"
276  # if wcount == 1:
277  # hal = "right"
278 
279  t='%s wid:%d ch:%d' %(["beg","end"][wcount], wire.ident, wire.channel)
280  ax.text(x, y, t,
281  horizontalalignment=hal,
282  bbox=dict(facecolor='yellow', alpha=0.5, pad=10))
283 
284  # if anode.ident==1 and face.ident==1:
285  # for wcount, wind in enumerate(plane.wires):
286  # wire = store.wires[wind]
287  # if plane.ident==0 and wire.ident == 491 or \
288  # plane.ident==1 and wire.ident == 514 or \
289  # plane.ident==2 and wire.ident == 256:
290 
291  # p1 = store.points[wire.tail]
292  # p2 = store.points[wire.head]
293  # print("A1, F1, P%d, W: %s [%s -> %s]" %(plane.ident, wire, p1, p2))
294 
295 
296 
297  ax.set_xlabel("Z [meter]")
298  ax.set_ylabel("Y [meter]")
299  ax.set_title("AnodeID %d, FaceID %d, PlaneID %d every %dth wire, x=%.3fm" % \
300  (anode.ident, face.ident, plane.ident, wire_step, wirex))
301  pdf.savefig(fig)
302  plt.close()
303  continue # over planes
304 
305  # plot directions
306  fig, axes = plt.subplots(nrows=2, ncols=2)
307  for iplane,winds in enumerate(first_wires):
308  plane_color = "red green blue".split()[iplane]
309  w0 = store.wires[winds[0]]
310  h0 = numpy.asarray(store.points[w0.head])
311  t0 = numpy.asarray(store.points[w0.tail])
312  w1 = store.wires[winds[1]]
313  h1 = numpy.asarray(store.points[w1.head])
314  t1 = numpy.asarray(store.points[w1.tail])
315 
316  c0 = 0.5*(h0+t0)
317  c1 = 0.5*(h1+t1)
318 
319  # print ("A%d F%d c0=%s c1=%s" % (anode.ident, face.ident, c0, c1))
320  # print (winds[0], w0, winds[1], w1)
321  # print ("h0=%s t0=%s" % (h0, t0))
322  # print ("h1=%s t1=%s" % (h1, t1))
323 
324  w = h0-t0 # wire direction
325  w = w/math.sqrt(numpy.sum(w*w))
326 
327  r = c1-c0 # roughly in the pitch direction
328  r = r/math.sqrt(numpy.sum(r*r))
329 
330  x = numpy.cross(w,r) # drift direction
331  x = x/math.sqrt(numpy.sum(x*x))
332  p = numpy.cross(x,w) # really pitch direction
333 
334  for ipt, pt in enumerate([x,w,p]):
335  axes[0,0].arrow(0,0, pt[2], pt[1], color=plane_color, linewidth=ipt+1) # 0
336  axes[0,1].arrow(0,0, pt[2], pt[0], color=plane_color, linewidth=ipt+1) # 1
337  axes[1,0].arrow(0,0, pt[0], pt[1], color=plane_color, linewidth=ipt+1) # 2
338 
339  for a,t in zip(axes.flatten(),["Z vs Y","X vs Y","Z vs X","none"]):
340  a.set_aspect('equal')
341  a.set_title("%s anode: %d, face: %d" % (t, anode.ident, face.ident))
342  a.set_ylim(-1.1,1.1)
343  a.set_xlim(-1.1,1.1)
344 
345  plt.tight_layout()
346  pdf.savefig(fig)
347  plt.close()
348 
349  continue # over faces
350 
351 
352 
353  fig, ax = plt.subplots(nrows=1, ncols=1)
354  ax.scatter(wire_z1, wire_x1,s=1, c=wire_anode, marker='.')
355  ax.set_title("AnodeID %d wires, tail" % anode.ident)
356  pdf.savefig(fig)
357  plt.close()
358  fig, ax = plt.subplots(nrows=1, ncols=1)
359  ax.scatter(wire_z2, wire_x2,s=1, c=wire_anode, marker='.')
360  ax.set_title("AnodeID %d wires, head" % anode.ident)
361  pdf.savefig(fig)
362  plt.close()
363  all_wire_x1 += wire_x1
364  all_wire_z1 += wire_z1
365  all_wire_x2 += wire_x2
366  all_wire_z2 += wire_z2
367  all_wire_anode += wire_anode
368 
369  continue # over anodes
370  fig, ax = plt.subplots(nrows=1, ncols=1)
371  ax.scatter(all_wire_z1, all_wire_x1,s=1, c=all_wire_anode,marker='.')
372  ax.set_title("All wires, tail")
373  pdf.savefig(fig)
374  plt.close()
375  fig, ax = plt.subplots(nrows=1, ncols=1)
376  ax.scatter(all_wire_z2, all_wire_x2,s=1, c=all_wire_anode,marker='.')
377  ax.set_title("All wires, head")
378  pdf.savefig(fig)
379  plt.close()
380 
381 
382 def face_in_allplanes(store, iface=0, segments=None):
383  fig = plt.figure()
384  ax = fig.add_subplot(111)
385 
386  face = store.faces[iface]
387  for planeid in face.planes:
388  plane = store.planes[planeid]
389  for wind in plane.wires[::20]:
390  wire = store.wires[wind]
391  if segments and not wire.segment in segments:
392  continue
393 
394  p1 = store.points[wire.tail]
395  p2 = store.points[wire.head]
396 
397  width = wire.segment + 1
398  ax.plot((p1.z, p2.z), (p1.y, p2.y), linewidth = width)
399 
400  return fig,ax
401 
402 def allwires(store):
403  fig = plt.figure()
404  ax = fig.add_subplot(111)
405 
406  iplane=0
407  plane = store.planes[store.faces[0].planes[iplane]]
408  wires = [store.wires[w] for w in plane.wires]
409  nwires = len(wires)
410 
411  cmap = plt.get_cmap('seismic')
412  colors = [cmap(i) for i in numpy.linspace(0, 1, nwires)]
413 
414  for iwire, wire in enumerate(wires):
415  p1 = store.points[wire.tail]
416  p2 = store.points[wire.head]
417 
418  color = colors[iwire]
419  ax.plot((p1.z, p2.z), (p1.y, p2.y), color=color)
420 
421  return fig,ax
422 
423 
424 
425 
426 
427 
428 
429 
430 
431 
432 
433 
434 
435 
436 
437 
438 
439 
440 
441 
442 
443 def plot_rect(rect, color="black"):
444  ax = plt.axes()
445  ax.add_patch(mpatches.Rectangle(rect.ll, rect.width, rect.height,
446  color=color, fill=False))
447  ax.set_xlabel("APA-local Z")
448  ax.set_ylabel("APA-local Y")
449  ax.set_title("Looking in anti-drift direction")
450 
451 
452 def plot_polyline(pts):
453  cmap = plt.get_cmap('seismic')
454  npts = len(pts)
455  colors = [cmap(i) for i in numpy.linspace(0, 1, npts)]
456  for ind, (p1, p2) in enumerate(zip(pts[:-1], pts[1:])):
457  x = numpy.asarray((p1.x, p2.x))
458  y = numpy.asarray((p1.y, p2.y))
459  plt.plot(x, y, linewidth=ind+1)
460 
461 
462 def plotwires(wires):
463  cmap = plt.get_cmap('seismic')
464  nwires = len(wires)
465 
466  chans = [w[2] for w in wires]
467  minchan = min(chans)
468  maxchan = max(chans)
469  nchans = maxchan - minchan + 1
470 
471  colors = [cmap(i) for i in numpy.linspace(0, 1, nchans)]
472  for ind, one in enumerate(wires):
473  pitch, side, ch, seg, p1, p2 = one
474  linestyle = 'solid'
475  if side < 0:
476  linestyle = 'dashed'
477  color = colors[ch-minchan]
478 
479  x = numpy.asarray((p1.x, p2.x))
480  y = numpy.asarray((p1.y, p2.y))
481  plt.plot(x, y, color=color, linewidth = seg+1, linestyle = linestyle)
482 
483 def plot_wires_sparse(wires, indices, group_size=40):
484  for ind in indices:
485  plotwires([w for w in wires if w[2]%group_size == ind])
486 
487 
488 def plot_some():
489  rect = Rectangle(6.0, 10.0)
490  plt.clf()
491  direc = Point(1,-1);
492  for offset in numpy.linspace(.1, 6, 60):
493  start = Point(-3.0 + offset, 5.0)
494  ray = Ray(start, start+direc)
495  pts = wrap_one(ray, rect)
496  plot_polyline(pts)
497 
498 
499 
500 
501 def plot_wires(wobj, wire_filter=None):
502  bbmin, bbmax = wobj.bounding_box
503  xmin, xmax = bbmin[2],bbmax[2]
504  ymin, ymax = bbmin[1],bbmax[1]
505  dx = xmax-xmin
506  dy = ymax-ymin
507  wires = wobj.wires
508 
509  #print (xmin,ymin), (dx,dy)
510  #print bbmin, bbmax
511 
512  wirenums = [w.wire for w in wires]
513  minwire = min(wirenums)
514  maxwire = max(wirenums)
515  nwires = maxwire-minwire+1
516 
517  if wire_filter:
518  wires = [w for w in wires if wire_filter(w)]
519  print ("filter leaves %d wires" % len(wires))
520  ax = plt.axes()
521  ax.set_aspect('equal', 'box') #'datalim')
522  ax.add_patch(mpatches.Rectangle((xmin, ymin), dx, dy,
523  color="black", fill=False))
524 
525  cmap = plt.get_cmap('rainbow') # seismic is bluewhitered
526 
527  colors = [cmap(i) for i in numpy.linspace(0, 1, nwires)]
528  for ind, one in enumerate(wires):
529  color = colors[one.wire-minwire]
530  x = numpy.asarray((one.beg[2], one.end[2]))
531  y = numpy.asarray((one.beg[1], one.end[1]))
532  plt.plot(x, y, color=color)
533 
534  plt.plot([ xmin + 0.5*dx ], [ ymin + 0.5*dy ], "o")
535 
536  plt.axis([xmin,xmax,ymin,ymax])
537 
def plot_polyline(pts)
Definition: plot.py:9
def plot_wires(wobj, wire_filter=None)
Definition: plot.py:501
def select_channels(store, pdffile, channels, labels=True)
Definition: plot.py:74
def plot_wires_sparse(wires, indices, group_size=40)
Definition: plot.py:483
auto enumerate(Iterables &&...iterables)
Range-for loop helper tracking the number of iteration.
Definition: enumerate.h:69
def plotwires(wires)
Definition: plot.py:462
T abs(T value)
def face_in_allplanes(store, iface=0, segments=None)
Definition: plot.py:382
def wrap_one(start_ray, rect)
Definition: generator.py:166
def allwires(store)
Definition: plot.py:402
def oneplane(store, iplane, segments=None)
Definition: plot.py:20
def plot_rect(rect, color="black")
Definition: plot.py:443
static int max(int a, int b)
def allplanes(store, pdffile)
Definition: plot.py:121
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
auto zip(Iterables &&...iterables)
Range-for loop helper iterating across many collections at the same time.
Definition: zip.h:295
void split(std::string const &s, char c, OutIter dest)
Definition: split.h:35
static QCString str