Index: ct/tools/pyscripts/camplotter.py =================================================================== --- /fact/tools/pyscripts/camplotter.py (revision 13119) +++ (revision ) @@ -1,134 +1,0 @@ -import matplotlib.pyplot as plt -import numpy as np -import numpy.ma as ma -import os, sys - -class camplotter( object ): - def __init__( self ): - chid, y,x,ye,xe,yh,xh,softid,hardid = np.loadtxt("map_dn.txt",unpack=True) - self.fig = plt.figure() - self.xe = xe - self.ye = ye - self.H = (6,0,30./180.*3.1415926) - - coor2chid = {} - chid2coor = [] - for i in range ( len(chid) ): - coor2chid[ (x[i],y[i]) ] = chid[i] - chid2coor.append( (x[i],y[i]) ) - -# print coor2chid -# print chid2coor - - - -######################################################################### -######################################################################### - -if __name__ == '__main__': - - inputfname = 'nofile.npz' - - if (len(sys.argv) > 1): - inputfname = sys.argv[1] - else: - print 'Usage ', sys.argv[0], 'input-file-path' - exit (0) - - cplt = camplotter() - - # np.savez ( filename, amplitude=maxAmp, time=maxPos, integral=integ) - npz = np.load( inputfname ) - print npz.files - - amp = npz['amplitude'] - int = npz['integral'] - tim = npz['time'] - - files = [] - fig = plt.figure(figsize=(24,16)) -# plt.title('left amplitude ... right integral') - ax1 = fig.add_subplot(231, aspect='equal') - ax2 = fig.add_subplot(232, aspect='equal') - ax3 = fig.add_subplot(233, aspect='equal') - ax4 = fig.add_subplot(234) - ax5 = fig.add_subplot(235) - ax6 = fig.add_subplot(236) - - ax1.grid(True) - ax2.grid(True) - ax3.grid(True) - ax4.grid(True) - ax5.grid(True) - ax6.grid(True) - - thr = 50 - hithr = 80 - lothr = 40 - print '------------ hiTHR=', hithr, 'loTHR=', lothr ,'-------------------' - print '------------', len(amp[:,0]), '-------------------' - -# for i in range(len(data[:,0])): - os.system("mkdir -p pngs") - for i in range(50): - - mamp = ma.masked_less(amp[i,:], thr) - mx = ma.masked_where(mamp.mask == True, cplt.xe) - my = ma.masked_where(mamp.mask == True, cplt.ye) - - core = ma.masked_less(amp[i,:], hithr) - edgecand = ma.masked_less(amp[i,:], lothr) - #edgecand = ma.masked_where(core.mask == True, edgecand) - - mxc = ma.masked_where(core.mask == True, cplt.xe).compressed() - myc = ma.masked_where(core.mask == True, cplt.ye).compressed() - mxe = ma.masked_where(edgecand.mask == True, cplt.xe).compressed() - mye = ma.masked_where(edgecand.mask == True, cplt.ye).compressed() - - - ax1.cla() - ax2.cla() - ax3.cla() - ax1.scatter(cplt.xe,cplt.ye,s=65,alpha=0.75,marker=cplt.H, c=amp[i,:], edgecolors='none', label='amp') - #ax1.scatter(mx.compressed() , my.compressed() ,s=65,alpha=0.75,marker=cplt.H, facecolors='none', linewidths=3, edgecolors='r') - if ( len(mxe) > 0 ): - ax1.scatter(mxe , mye ,s=65,alpha=0.75,marker=cplt.H, facecolors='none', linewidths=2, edgecolors='r') - if ( len(mxc) > 0 ): - ax1.scatter(mxc , myc ,s=65,alpha=0.75,marker=cplt.H, facecolors='none', linewidths=3, edgecolors='k') - ax1.axis([-22,22,-22,22]) - - ax2.scatter(cplt.xe,cplt.ye,s=70,alpha=0.75,marker=cplt.H, c=int[i,:], edgecolors='none', label='inte') - ax2.axis([-22,22,-22,22]) - ax3.scatter(cplt.xe,cplt.ye,s=70,alpha=0.75,marker=cplt.H, c=tim[i,:], edgecolors='none', label='time') - ax3.axis([-22,22,-22,22]) - - ax4.cla() - ax4.hist(amp[i,:], 100, facecolor='r', alpha=0.75) - ax4.set_yscale("log", nonposy='clip') - - ax5.cla() - ax5.hist(int[i,:], 100, facecolor='g', alpha=0.75) - ax5.set_yscale("log", nonposy='clip') - - ax6.cla() - ax6.hist(tim[i,:], 100,facecolor='b', alpha=0.75) - ax6.set_yscale("log", nonposy='clip') - - - - fname = 'pngs/_tmp%03d.png'%i - print 'Saving frame', fname - fig.savefig(fname) - - files.append(fname) - - # in order to make a movie using mencoder -- uncomment the following two lines - # - #print 'Making movie animation.mpg - this make take a while' - #os.system("mencoder 'mf://pngs/_tmp*.png' -mf type=png:fps=2 -ovc lavc -lavcopts vcodec=wmv2 -oac copy -o animation.mpg") - - # in order to make an animated gif using ImageMagick -- uncomment the following two lines - # - #print 'Making animated gif animation.gif - this make take a while' - #os.system("convert -delay 33 -loop 0 pngs/_tmp*.png animation.gif") - Index: /fact/tools/pyscripts/examples/how_to_make_movies.py =================================================================== --- /fact/tools/pyscripts/examples/how_to_make_movies.py (revision 13120) +++ /fact/tools/pyscripts/examples/how_to_make_movies.py (revision 13120) @@ -0,0 +1,134 @@ +import matplotlib.pyplot as plt +import numpy as np +import numpy.ma as ma +import os, sys + +class camplotter( object ): + def __init__( self ): + chid, y,x,ye,xe,yh,xh,softid,hardid = np.loadtxt("map_dn.txt",unpack=True) + self.fig = plt.figure() + self.xe = xe + self.ye = ye + self.H = (6,0,30./180.*3.1415926) + + coor2chid = {} + chid2coor = [] + for i in range ( len(chid) ): + coor2chid[ (x[i],y[i]) ] = chid[i] + chid2coor.append( (x[i],y[i]) ) + +# print coor2chid +# print chid2coor + + + +######################################################################### +######################################################################### + +if __name__ == '__main__': + + inputfname = 'nofile.npz' + + if (len(sys.argv) > 1): + inputfname = sys.argv[1] + else: + print 'Usage ', sys.argv[0], 'input-file-path' + exit (0) + + cplt = camplotter() + + # np.savez ( filename, amplitude=maxAmp, time=maxPos, integral=integ) + npz = np.load( inputfname ) + print npz.files + + amp = npz['amplitude'] + int = npz['integral'] + tim = npz['time'] + + files = [] + fig = plt.figure(figsize=(24,16)) +# plt.title('left amplitude ... right integral') + ax1 = fig.add_subplot(231, aspect='equal') + ax2 = fig.add_subplot(232, aspect='equal') + ax3 = fig.add_subplot(233, aspect='equal') + ax4 = fig.add_subplot(234) + ax5 = fig.add_subplot(235) + ax6 = fig.add_subplot(236) + + ax1.grid(True) + ax2.grid(True) + ax3.grid(True) + ax4.grid(True) + ax5.grid(True) + ax6.grid(True) + + thr = 50 + hithr = 80 + lothr = 40 + print '------------ hiTHR=', hithr, 'loTHR=', lothr ,'-------------------' + print '------------', len(amp[:,0]), '-------------------' + +# for i in range(len(data[:,0])): + os.system("mkdir -p pngs") + for i in range(50): + + mamp = ma.masked_less(amp[i,:], thr) + mx = ma.masked_where(mamp.mask == True, cplt.xe) + my = ma.masked_where(mamp.mask == True, cplt.ye) + + core = ma.masked_less(amp[i,:], hithr) + edgecand = ma.masked_less(amp[i,:], lothr) + #edgecand = ma.masked_where(core.mask == True, edgecand) + + mxc = ma.masked_where(core.mask == True, cplt.xe).compressed() + myc = ma.masked_where(core.mask == True, cplt.ye).compressed() + mxe = ma.masked_where(edgecand.mask == True, cplt.xe).compressed() + mye = ma.masked_where(edgecand.mask == True, cplt.ye).compressed() + + + ax1.cla() + ax2.cla() + ax3.cla() + ax1.scatter(cplt.xe,cplt.ye,s=65,alpha=0.75,marker=cplt.H, c=amp[i,:], edgecolors='none', label='amp') + #ax1.scatter(mx.compressed() , my.compressed() ,s=65,alpha=0.75,marker=cplt.H, facecolors='none', linewidths=3, edgecolors='r') + if ( len(mxe) > 0 ): + ax1.scatter(mxe , mye ,s=65,alpha=0.75,marker=cplt.H, facecolors='none', linewidths=2, edgecolors='r') + if ( len(mxc) > 0 ): + ax1.scatter(mxc , myc ,s=65,alpha=0.75,marker=cplt.H, facecolors='none', linewidths=3, edgecolors='k') + ax1.axis([-22,22,-22,22]) + + ax2.scatter(cplt.xe,cplt.ye,s=70,alpha=0.75,marker=cplt.H, c=int[i,:], edgecolors='none', label='inte') + ax2.axis([-22,22,-22,22]) + ax3.scatter(cplt.xe,cplt.ye,s=70,alpha=0.75,marker=cplt.H, c=tim[i,:], edgecolors='none', label='time') + ax3.axis([-22,22,-22,22]) + + ax4.cla() + ax4.hist(amp[i,:], 100, facecolor='r', alpha=0.75) + ax4.set_yscale("log", nonposy='clip') + + ax5.cla() + ax5.hist(int[i,:], 100, facecolor='g', alpha=0.75) + ax5.set_yscale("log", nonposy='clip') + + ax6.cla() + ax6.hist(tim[i,:], 100,facecolor='b', alpha=0.75) + ax6.set_yscale("log", nonposy='clip') + + + + fname = 'pngs/_tmp%03d.png'%i + print 'Saving frame', fname + fig.savefig(fname) + + files.append(fname) + + # in order to make a movie using mencoder -- uncomment the following two lines + # + #print 'Making movie animation.mpg - this make take a while' + #os.system("mencoder 'mf://pngs/_tmp*.png' -mf type=png:fps=2 -ovc lavc -lavcopts vcodec=wmv2 -oac copy -o animation.mpg") + + # in order to make an animated gif using ImageMagick -- uncomment the following two lines + # + #print 'Making animated gif animation.gif - this make take a while' + #os.system("convert -delay 33 -loop 0 pngs/_tmp*.png animation.gif") + Index: /fact/tools/pyscripts/examples/introductory_course/produce_npz.py =================================================================== --- /fact/tools/pyscripts/examples/introductory_course/produce_npz.py (revision 13120) +++ /fact/tools/pyscripts/examples/introductory_course/produce_npz.py (revision 13120) @@ -0,0 +1,77 @@ +#!/usr/bin/python +# +# Dominik Neise +# TU Dortmund +# +# example for storing numpy arrays containing FACT analysis results +# to file + +#import this at first, please +from pyfact_rename import RawData + +import os.path +import numpy as np + +from drs_spikes import DRSSpikes +from fir_filter import SlidingAverage +from extractor import GlobalMaxFinder +from cleaners import AmplitudeCleaner + +############################################################################## +confirm_next_step = False# this is for user interaction + +data_file_name = '/media/daten_platte/FACT/data/20120229_144.fits.gz' +if not os.path.isfile(data_file_name): + print 'not able to find file:', data_file_name + sys.exit(-1) +calib_file_name = '/media/daten_platte/FACT/data/20120229_132.drs.fits.gz' +if not os.path.isfile(calib_file_name ): + print 'not able to find file:', calib_file_name + sys.exit(-1) + +run = RawData(data_file_name, calib_file_name) +despike = DRSSpikes() +smooth = SlidingAverage(8) +extract = GlobalMaxFinder(40,200) +cleaner = AmplitudeCleaner(45, 18) + +areas = [] +sizes = [] +for data,startcell,tt in run: + # trigger type 4 means 'physics event' + if tt==4: + data = despike(data) + data = smooth(data) + amplitude, time_of_max = extract(data) + survivors = cleaner(amplitude, return_bool_mask=False) + + # this is not python like, should be done in a single line + # adding up the amplitude of all survivors + size = 0 + for pixel in survivors: + size += amplitude[pixel] + + area = len(survivors) + + if area > 0: + areas.append( area ) + sizes.append( size ) + + # This is for ---------- USER INTERACTION ------------------------- + if confirm_next_step: + user_input = raw_input("'q'-quit, 'r'-run, anything else goes one step") + number=None + try: + number=int(user_input) + except: + number=None + if user_input.find('q') != -1: + sys.exit(0) + elif user_input.find('r') != -1: + confirm_next_step = False + elif number!=None: + run += number + # ---------------END OF USER INTERACTION ------------------------- + +outfname = raw_input("please enter name for outputfile") +np.savez(outfname, Area=areas, Size=sizes) Index: ct/tools/pyscripts/examples/test_dneise.py =================================================================== --- /fact/tools/pyscripts/examples/test_dneise.py (revision 13119) +++ (revision ) @@ -1,52 +1,0 @@ -#!/usr/bin/python -# -# Dominik Neise -# TU Dortmund -# -# test for freshly implemented methods in pyfact. -# subject to frequent change! - -from pyfact import * - -dfname = '/data03/fact-construction/raw/2011/11/24/20111124_121.fits' -calfname = '/data03/fact-construction/raw/2011/11/24/20111124_111.drs.fits' -# access the data -rd = rawdata( dfname, calfname ) - -print 'dfname: ', dfname -print 'calfname:', calfname -print 'NEvents: ', rd.NEvents - -def loop_acal( Neve = 1000 ): - """ - bla - """ - if rd.NEvents < Neve: - print 'data file contains not enough events' - exit( 0 ) - else: - for i in range( Neve ): - #if np.mod(i,1000): - # print 'Event: ', i - rd.next() - - rd.filterSlidingAverage() - print rd.smoothData - rd.filterCFD() - print rd.cfdData - rd.findPeak() - print rd.maxPos - print rd.maxAmp - rd.sumAroundPeak() - print rd.integral - #print len(rd.maxPos) - #print len(rd.maxAmp) - print 'Trigger Type', rd.trigType - -print (rd.filterSlidingAverage.__doc__) -print (rd.filterCFD.__doc__) -print (rd.findPeak.__doc__) -print (rd.sumAroundPeak.__doc__) - -loop_acal(2) - Index: ct/tools/pyscripts/produce_npz.py =================================================================== --- /fact/tools/pyscripts/produce_npz.py (revision 13119) +++ (revision ) @@ -1,64 +1,0 @@ -#!/usr/bin/python -# -# Dominik Neise -# TU Dortmund -# -# test for freshly implemented methods in pyfact. -# subject to frequent change! - -from pyfact import * -#from savecam import * - - -dfname = '/data00/fact-construction/raw/2011/11/24/20111124_081.fits.gz' -outfname = '20111124_081.npz' -calfname = '/data00/fact-construction/raw/2011/11/24/20111124_071.drs.fits.gz' - -# access the data -rd = rawdata( dfname, calfname ) - -print 'dfname: ', dfname -print 'calfname:', calfname -print 'NEvents: ', rd.NEvents - -def loop_acal( Neve = 1000 ): - """ - bla - """ - if rd.NEvents < Neve: - print 'data file contains not enough events' - exit( 0 ) - else: - for i in range( Neve ): - if (np.mod(i,100)==0) : - print 'Event: ', i - rd.next() - - rd.filterSlidingAverage() -# print rd.smoothData - rd.filterCFD() -# print rd.cfdData - rd.findPeak() - list_maxPos.append(rd.maxPos) - list_maxAmp.append(rd.maxAmp) - rd.sumAroundPeak() - list_integral.append(rd.integral) -# plotincam(rd.integral, 'evt_'+str(i)+'.pdf') - #print len(rd.maxPos) - #print len(rd.maxAmp) -# print 'Trigger Type', rd.trigType - -#print (rd.filterSlidingAverage.__doc__) -#print (rd.filterCFD.__doc__) -#print (rd.findPeak.__doc__) -#print (rd.sumAroundPeak.__doc__) - -list_maxPos = [] -list_maxAmp = [] -list_integral = [] -loop_acal(rd.NEvents) -maxAmp = np.vstack(list_maxAmp) -maxPos = np.vstack(list_maxPos) -integ = np.vstack(list_integral) - -np.savez(outfname, amplitude=maxAmp, time=maxPos, integral=integ)