Changeset 13120 for fact/tools/pyscripts/examples/introductory_course/produce_npz.py

Timestamp:

03/15/12 13:57:33 (13 years ago)

Author:

neise

Message:

tidied up a bit. moved examples to examples and deleted scripts, I never used anymore

Location:

fact/tools/pyscripts/examples/introductory_course

Files:

• . (added)
• produce_npz.py (moved) (moved from fact/tools/pyscripts/produce_npz.py ) (1 diff)

fact/tools/pyscripts/examples/introductory_course/produce_npz.py

@@ -4 +4 @@
 # TU Dortmund
 #
-# test for freshly implemented methods in pyfact.
-# subject to frequent change!
+# example for storing numpy arrays containing FACT analysis results
+# to file
 
-from pyfact import *
-#from savecam import *
+#import this at first, please
+from pyfact_rename import RawData
 
+import os.path
+import numpy as np
 
-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'
+from drs_spikes import DRSSpikes
+from fir_filter import SlidingAverage
+from extractor  import GlobalMaxFinder
+from cleaners   import AmplitudeCleaner
 
-# access the data
-rd = rawdata( dfname, calfname )
+##############################################################################
+confirm_next_step = False# this is for user interaction
 
-print 'dfname:  ', dfname
-print 'calfname:', calfname
-print 'NEvents: ', rd.NEvents
+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)
 
-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
+run = RawData(data_file_name, calib_file_name)
+despike = DRSSpikes()
+smooth = SlidingAverage(8)
+extract = GlobalMaxFinder(40,200)
+cleaner = AmplitudeCleaner(45, 18)
 
-#print (rd.filterSlidingAverage.__doc__)
-#print (rd.filterCFD.__doc__)
-#print (rd.findPeak.__doc__)
-#print (rd.sumAroundPeak.__doc__)
+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)
 
-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)
+        # 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 -------------------------
 
-np.savez(outfname, amplitude=maxAmp, time=maxPos, integral=integ)
+outfname = raw_input("please enter name for outputfile")
+np.savez(outfname, Area=areas, Size=sizes)