source: fact/tools/pyscripts/sandbox/vogler/CalFitsPerformanceWeave.py@ 16858

#!/usr/bin/python -tt
# ********************************
# Test script for the CalFits class
# 
# written by Thomas Kraehenbuehl, ETH Zurich
# tpk@phys.ethz.ch, +41 44 633 3973
# April 2012
# ********************************

datafilename = '/fact/raw/2012/04/17/20120417_004.fits.gz'
calibfilename = '/fact/raw/2012/04/17/20120417_003.drs.fits.gz'

import numpy as np
from scipy import weave
from scipy.weave import converters

from ROOT import gSystem
gSystem.Load("calfits_h.so")
from ROOT import *
print "Testing object creation: "
caltest = CalFits(datafilename,calibfilename)
npcalevent  = np.empty( caltest.npix * caltest.nroi, np.float64) #.reshape(caltest.npix ,caltest.nroi)
caltest.SetNpcaldataPtr(npcalevent)

print "Common variables:"
print "ROI: ", caltest.nroi
print "#Pix: ", caltest.npix
print "Number of events: ", caltest.nevents
print

print "Information per Event:"
caltest.GetCalEvent()
print "Event ID: ", caltest.event_id
print "Trigger type: ", caltest.event_triggertype
print "Uncalibrated data: ", caltest.event_data
print "Calibrated data: ", caltest.npcaldata
print "Board times: ", caltest.event_boardtimes
print "Trigger offsets: ", caltest.event_offset
print

print "Examples of other information"
print "Calibfile ROI: ", caltest.calib_nroi
print "Column size BaselineMean: ", caltest.calibfile.GetN("BaselineMean")
print "Datafile ROI: ", caltest.data_nroi
print "Data: ", caltest.datafile.GetN("Data")
print "StartCellData: ", caltest.datafile.GetN("StartCellData")
print "Direct datafile access: ", caltest.datafile.GetN("StartCellData")
print
print "Columns of the datafile: "
caltest.datafile.PrintColumns()

#AND WE HAVE A WINNER: 43 Hz with scipy.weave!
while caltest.GetCalEvent():
    c_integrals = np.zeros(40000, np.float64) #Allocate the memory for about 40000 singles
    if caltest.event_id>300:
        break
    
    cppcode = """
    //ToDo: get nroi and npix, check singles<40000
    double slice, lastslice, integral;
    int singles = 0;
    double threshold = 2.5;
    for(int px=0; px<npix; px++) {
        lastslice = npcalevent(px*nroi);
        for(int sl=1; sl<nroi-14; sl++) {
            slice = npcalevent(px*nroi+sl);
            if((lastslice<threshold)&&(slice>threshold)) {
                integral = 0;
                for(int l=5; l<15; l++) {
                    integral += npcalevent(px*nroi+sl+l);
                }
                c_integrals(singles) = integral;
                ++singles;
//                std::cout << sl << " ";
            }
            lastslice = slice;
        }
    }
//    std::cout << std::endl << singles << std::endl;
    """
    
    #Passing the caltest object to weave is not possible
    #The int() statement is necessary to add additional data type information to the variable (ie. otherwise some things would work in the C-code, other would not)
    nroi = int(caltest.nroi)
    npix = int(caltest.npix)
    
    weave.inline(cppcode, ['c_integrals', 'npcalevent','nroi','npix'], type_converters=converters.blitz)
    
#    print "CalFitsPerformanceWeave"
#    for i in range(10):
#        print i, npcalevent[i], c_integrals[i]
#    for x in c_integrals:
#        if x:
#            print x,

#    #Old attempts using numpy functions
#    npcalevent = npcalevent.reshape((caltest.npix, caltest.nroi))
#    npthr = npcalevent>2.5 #true if above threshold
#    npthr = npthr[:,1:-14] * np.logical_not(npthr[:,:-15]) #only true for the zero crossings, shape (1400,285) [smaller due to necessary integration range]
##    print [(x,y) for x,y in zip(npthr.nonzero()[0],npthr.nonzero()[1])] #print the coordinates, range 0-1399,0-284
#    
#    #Various versions to get the integral, all with approximately the same miserable speed (3 Hz), except for the last one (vectorized function) with ~4.3 Hz
#    #Missing: add deadtime after an integration, remove start & end of the array
##    integrals = [np.sum(npcalevent[x,y+5:y+15]) for x,y in np.transpose(npthr.nonzero())]
##    integrals = [np.sum(npcalevent[x,y+5:y+15]) for x,y in zip(npthr.nonzero()[0],npthr.nonzero()[1])]
##    integrals = [np.sum(npcalevent[x,y+5:y+15]) for x,y in itertools.izip(npthr.nonzero()[0],npthr.nonzero()[1])]
##    integrals = map((lambda index_a,index_b: np.sum(npcalevent[index_a,index_b+5:index_b+15])),npthr.nonzero()[0],npthr.nonzero()[1])
#    integrals = vecoffsum(npthr.nonzero()[0],npthr.nonzero()[1])
#    print len(integrals)

#    print caltest.event_id, caltest.event_triggertype, caltest.event_caldata[10]
#    pass

del caltest