#!/usr/bin/python -itt
#
# Dominik Neise
#
# cleaning a small step towards the truth
from pyfact import RawData
import os.path
import matplotlib.pyplot as plt
import numpy as np
from fir_filter import *
from extractor import *
from drs_spikes import *
from plotters import *
import time as t
from coor import Coordinator
confirm_next_step = False# this is for user interaction

data_file_name = '/media/daten_platte/FACT/data/20120229_144.fits.gz'
calib_file_name = '/media/daten_platte/FACT/data/20120229_132.drs.fits.gz'
if not os.path.isfile(data_file_name):
    print 'not able to find file:', data_file_name
    sys.exit(-1)
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)

#plotA = CamPlotter('amplitudes')
#plotT = CamPlotter('times')
#plotCA = CamPlotter('cleaned amplitudes')

#plotArea = HistPlotter('area', 1440, (0,1440) )
#plotSize = HistPlotter('size', 1000, (0,10000) )

#funnyPlot = Plotter('area vs log(size')


coreTHR = 45 # copied from F. Temme
edgeTHR = 18 # copied from F. Temme

# get dictionary of next neighbors
nn = Coordinator().nn

print 'Core threshold=', coreTHR
print 'Edge threshold=', edgeTHR
print '*'*70

areas = []
sizes = []

for data,startcell,tt in run:
    if tt==4:
        data = despike(data)
        data = smooth(data)
        amplitude, time_of_max = extract(data)
        
        #plotA.name='amplitudes EvtID:' + str(run.event_id.value) + ' TT:' + str(tt)
        #plotA(amplitude)
        #plotT(time_of_max)
        
        # Here we start the cleaning ... just like that...
        #print 'cleaning away all pixel <' , coreTHR
        core_chid_candidates = np.where( amplitude > coreTHR)[0]
        core_chids = [] # coor chids, which survive Gauks step 1
        survivors = [] # coor & edge pixel 
        #print 'number of core candidates:', len(core_chid_candidates)
        
        #print 'throwing away all pixel w/o any neighbor' 
        # get rid of single coor pics
        for cand in core_chid_candidates:
            neighbor_found = False
            # loop over all neigbors of cand'idate
            for n in nn[cand]:
                if n in core_chid_candidates:
                    neighbor_found = True
            if neighbor_found:
                core_chids.append(cand)
                
        #print 'after deletion of single coor pixels'
        
        
        #add edge pixel to the edge of the coors
        #print 'resurrecting edge pixels ... i.e. all pixel >', edgeTHR
        survivors = core_chids[:]
        for coor in core_chids:
            for n in nn[coor]:
                # if neighbor is a core pixel, then do nothing
                if n in core_chids:
                    pass
                elif amplitude[n] > edgeTHR:
                    survivors.append(n)
        survivors = np.array(survivors, dtype=int)
        print '#survivors', len(survivors), 'evtID', run.event_id.value
        
        #plotCA(data=amplitude, mask=survivors)
        
        size = 0
        for pixel in survivors:
            size += amplitude[pixel]
        
        if len(survivors) > 0:
            areas.append( len(survivors) )
            sizes.append( size )
        
        
        #plotArea(areas, 'areas of ' + str(run.event_id.value) + 'events')
        #plotSize(sizes, 'sizes of ' + str(run.event_id.value) + 'events')
        
        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
                
        total_event_number = run.event_id.value

plt.ion()
myfig = plt.figure()
myn = myfig.number
logsize = np.log10(np.array(sizes))
areas = np.array(areas)

plt.figure(myn)
plt.title('area vs. log10(size) of '+ str(total_event_number) + 'events')
plt.xlabel('log10(size/1mV)')
plt.ylabel('area [#pixel]')
plt.plot( logsize,areas, '.')