Changeset 13096 for fact

Timestamp:

03/14/12 09:08:24 (13 years ago)

Author:

lusterma

Message:

modified working version - not finished

File:

• fact/tools/pyscripts/tools/spike_ana.py (modified) (6 diffs)

fact/tools/pyscripts/tools/spike_ana.py

@@ -4 +4 @@
 # ETH Zurich
 #
+from ROOT import TFile
 import numpy as np
 from pyfact_rename import RawData
@@ -9 +10 @@
 from drs_spikes import DRSSpikes
 
-from hist import histogramList, SaveHistograms
+from hist import hist, histogramList, SaveHistograms
 
 import matplotlib.pyplot as plt
 
-
-spike_histo = histogramList('spike_ana')
-spike_histo.h1d('N_candidates', 'N candidates', 
+spike_histo = {}
+spike_histo['cand'] = hist('N_candidates', 'N candidates',  
                 10000, 0., 10000., 'number of spike candidates', 'N / 1')
-spike_histo.h1d('N_singles', 'N doubles', 
+spike_histo['sing']= hist('N_singles', 'N singles', 
                 10000, 0., 10000., 'number of signel spikes', 'N / 1')
-spike_histo.h1d('N_doubles', 'N singles', 
+spike_histo['doub'] = hist('N_doubles', 'N doubles', 
                 10000, 0., 10000., 'number of double spikes', 'N / 1')
 
-spike_histo.h1d('ind_ampl_single', 'indicator amplitude singles', 
-                500, 0., 1000., 'indicator amplitude in mV', 'N / 2 mV')
-spike_histo.h1d('ind_ampl_double', 'indicator amplitude singles', 
-                500, 0., 1000., 'indicator amplitude in mV', 'N / 2 mV')
+
+spike_histo['all_ampl'] = hist('all_ampl', 'indicator all amplitude',                 1000, -1000., 1000., 'indicator amplitude in mV', 'N / 2 mV')
+spike_histo['cand_ampl'] = hist('cand_ampl', 'candidate amplitudes',                 500, 0., 1000., 'indicator amplitude in mV', 'N / 2 mV')
+spike_histo['sing_ampl'] = hist('cand_ampl_single', 'amplitude singles',                500, 0., 1000., 'indicator amplitude in mV', 'N / 2 mV')
+spike_histo['doub_ampl'] = hist('cand_ampl_double', 'amplitude doubles',                 500, 0., 1000., 'indicator amplitude in mV', 'N / 2 mV')
 
 spikes = {}
@@ -31 +32 @@
 spikes['sing'] = []
 spikes['doub'] = []
+spikes['cand_ampl'] = []
+spikes['sing_ampl'] = []
+spikes['doub_ampl'] = []
+#spikes['ind'] = []
 
 def spike_ana(candidates, singles, doubles, data, ind):
     """ call back function for DRSSpikes class """
-    
+
     #number of spikes
     #print len(candidates), len(singles), len(doubles)
@@ -40 +45 @@
     spikes['sing'].append(len(singles))
     spikes['doub'].append(len(doubles))
+    spike_histo['cand'](len(candidates))
+    spike_histo['sing'](len(singles))
+    spike_histo['doub'](len(doubles))
 
-    #print 'spikes[cand]', spikes['cand']
-    #spike_histo('N_candidates', len(candidates))
-    #spike_histo('N_singles', len(singles))
-    #spike_histo('N_doubles', len(doubles))
+    spike_histo['all_ampl'](ind)
+    spike_histo['cand_ampl'](ind[ np.array(candidates)+1])
+    spike_histo['sing_ampl'](ind[ np.array(singles)+1])
+    spike_histo['doub_ampl']( (ind[ np.array(doubles)+1] + ind[ np.array(doubles)+2]) / 2.)
 
-    #indicator
-    print [ind[singles[0]+i] for i in range(4)], '\n'
-    spikes['sing_ampl'].append(ind[singles+1])
+    spikes['cand_ampl'].extend( list(ind[ np.array(candidates)+1]) )
+    spikes['sing_ampl'].extend( list(ind[ np.array(singles)+1]) )
+    spikes['doub_ampl'].extend( list(
+            (ind[ np.array(doubles)+1] + ind[ np.array(doubles)+2]) / 2.) )
+    print 'len() ', len( spikes['sing_ampl'] )
+    #print spikes['sing_ampl']
+
     #spike_histo('ind_ampl_single', ind[single]
     #amplitudes
@@ -57 +69 @@
     #print singles[0], a, b, data[a,b-1], data[a,b], data[1,b+1]
 
-data_file_name = '/data00/fact-construction/raw/2011/11/24/20111124_111.fits.gz'
+data_file_name = '/data00/fact-construction/raw/2011/11/24/20111124_117.fits.gz'
 calib_file_name = '/data00/fact-construction/raw/2011/11/24/20111124_111.drs.fits.gz'
 run = RawData( data_file_name, calib_file_name )
 
-remove_spikes = DRSSpikes(7., user_action=spike_ana, debug=False)
+remove_spikes = DRSSpikes(5., user_action=spike_ana, debug=False)
 
-for event in range(10):
-    data, scell, tt = run.next()
+#for event in range(1000):
+ev = 0
+for data, scell, tt in run:
+    print 'event ', ev
+    ev += 1
+    if ev == 10:
+        break
     remove_spikes(data)
- 
+    
+
+root_file = TFile( 'spike_histo.root', 'RECREATE')
+    
+root_file.mkdir( 'spike_histo' )
+root_file.cd( 'spike_histo' )
+for key in spike_histo:
+    spike_histo[key].mkroot()
+    spike_histo[key].rh.Write()
+
+root_file.Close()
+
+plt.figure()
+plt.semilogy(spike_histo['all_ampl'].y, 'b', label = 'all')
+plt.semilogy(spike_histo['cand_ampl'].y, 'r', label='candidates')
+plt.title('Estimator values')
+plt.xlabel('estimator value')
+plt.ylabel('N / 1')
+plt.legend()
+plt.savefig('spike_histo.jpg')
+
 x = range(len(spikes['cand']))
-#print 'x: ', x
+
 plt.figure()
 plt.plot(x, np.array(spikes['cand']), 'b.', label='candidates')
 plt.plot(x, np.array(spikes['sing']), 'g.', label='singles')
 plt.plot(x, np.array(spikes['doub']), 'r.', label='doubles')
+plt.savefig('t.jpg')
 
 plt.title('Number of spikes')
@@ -94 +132 @@
 
 plt.figure()
+plt.hist(np.array(spikes['cand_ampl']), 100, color='blue',
+         label='candidates', histtype='step', log=True)
+plt.hist(np.array(spikes['sing_ampl']), 100, color='green',
+         label='singles', histtype='step', log=True)
+plt.hist(np.array(spikes['doub_ampl']), 100, color='red',
+         label='doubles', histtype='step', log=True)
+plt.grid(True)
+plt.title('Amplitude of spikes spike indicator')
+plt.legend()
+plt.savefig('spikes_ampl_hist.jpg')
 
 plt.show()
-
-
-#SaveHistograms( [spike_histo] ,'spike_histo.root')
-