Index: fact/tools/pyscripts/sandbox/kraehenb/CalFitsPerformance.py =================================================================== --- fact/tools/pyscripts/sandbox/kraehenb/CalFitsPerformance.py (revision 13469) +++ (revision ) @@ -1,78 +1,0 @@ -#!/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 -import itertools - -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() - -def offsum(x,y): - return np.sum(npcalevent[x,y+5:y+15]) -vecoffsum = np.vectorize(offsum) #specifying otypes=[np.float/double] does not improve speed - -while caltest.GetCalEvent(): - if caltest.event_id>30: - break - 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 - - integrals = [] - #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) -# for i in range(10): -# print i, npcalevent[0,i], integrals[i] -# print caltest.event_id, caltest.event_triggertype, caltest.event_caldata[10] -# pass - -del caltest Index: fact/tools/pyscripts/sandbox/kraehenb/CalFitsPerformanceWeave.py =================================================================== --- fact/tools/pyscripts/sandbox/kraehenb/CalFitsPerformanceWeave.py (revision 13469) +++ (revision ) @@ -1,114 +1,0 @@ -#!/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; pxthreshold)) { - 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 Index: fact/tools/pyscripts/sandbox/kraehenb/CalFitsTest.py =================================================================== --- fact/tools/pyscripts/sandbox/kraehenb/CalFitsTest.py (revision 13469) +++ (revision ) @@ -1,58 +1,0 @@ -#!/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 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 "Calibrated data in numpy array: ", npcalevent -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() - -#while caltest.GetCalEvent(): -# if caltest.event_id>10: -# break -# print caltest.event_id, caltest.event_triggertype, caltest.event_caldata[10] -# pass -#print - -del caltest Index: fact/tools/pyscripts/sandbox/kraehenb/calfits.h =================================================================== --- fact/tools/pyscripts/sandbox/kraehenb/calfits.h (revision 13469) +++ (revision ) @@ -1,188 +1,0 @@ -//******************************** -// -// Class Calfits -// Wrapper class for fits.h or pyfits.h -// Provides fast access to calibrated events of FACT raw files -// -// written by Thomas Kraehenbuehl, ETH Zurich -// tpk@phys.ethz.ch, +41 44 633 3973 -// April 2012 -// -//******************************** -// -// Compilation (root in pyfact directory) -// root [0] gSystem->Load("/usr/lib64/libz.so"); -// root [1] .L pyfits.h++ -// root [2] .L ../sandbox/kraehenb/calfits.h++ -// -// Usage in Python: -// See CalFitsTest.py -// -//******************************** - -//ToDo: shared library creation debuggen - -#ifndef CALFITS_H -#define CALFITS_H - -#include -#include - -#ifndef __MARS__ -#include -#include -#include -#define gLog cerr -#define ___err___ "" -#define ___all___ "" -#else -#include "MLog.h" -#include "MLogManip.h" -#define ___err___ err -#define ___all___ all -#endif - -#ifdef __EXCEPTIONS -#include -#endif - -#include "../../pyfact/pyfits.h" - -class CalFits -{ -public: - //No standard constructor CalFits()! - - //Direct handlers of the files - fits datafile, calibfile; //Class name should be PyFits or better FactPyFits... - - //Basic file parameters - UInt_t calib_nroi, calib_npix; - UInt_t calib_blm_size, calib_gm_size, calib_tom_size; - UInt_t data_nroi, data_npix, data_ndata; - - //Common variables - UInt_t nroi, npix, nevents; - - //Calibration variables - float* calib_baselinemean; - float* calib_gainmean; - float* calib_triggeroffsetmean; - //Using instead of arrays makes no visible difference - //ToDo: use arrays of size 1440x1024 (x2 for wrap-arounds) and read all variables into those - - //Event variables - UInt_t event_id; - UShort_t event_triggertype; - int16_t* event_data; - int16_t* event_offset; - int32_t* event_boardtimes; - double* npcaldata; - - CalFits(const string &datafilename, const string &calibfilename) //Constructor with two filenames - : datafile(datafilename), - calibfile(calibfilename), - npcaldata(NULL) - { - //Read basic parameters of the two files -// std::cout << "...Reading basic file parameters..." << std::endl; - calib_nroi = calibfile.GetUInt("NROI"); - calib_npix = calibfile.GetUInt("NPIX"); - data_nroi = datafile.GetUInt("NROI"); - data_npix = datafile.GetUInt("NPIX"); - data_ndata = datafile.GetN("Data"); - - calib_blm_size = calibfile.GetN("BaselineMean")/calib_npix; - calib_gm_size = calibfile.GetN("GainMean")/calib_npix; - calib_tom_size = calibfile.GetN("TriggerOffsetMean")/calib_npix; - -// std::cout << "Column sizes: " << calib_blm_size << " " << calib_gm_size << " " << calib_tom_size << std::endl; - - //Define the common variables - if((calib_nroi==data_nroi)&&(calib_npix==data_npix)&&(data_nroi*data_npix==data_ndata)&&(calib_blm_size==calib_gm_size)&&(calib_tom_size==calib_nroi)) { - nroi = data_nroi; - npix = data_npix; - } - else { - ostringstream str; - str << "Data/calib file error: NROI mismatch, NPIX mismatch, data column size wrong or calib columns mismatch."; -#ifdef __EXCEPTIONS - throw runtime_error(str.str()); -#else - gLog << ___err___ << "ERROR - " << str.str() << endl; - return; -#endif - } - nevents = datafile.GetNumRows(); - - //Read the calibration data -// std::cout << "...Reading calibration data..." << std::endl; - calib_baselinemean = new float[calibfile.GetN("BaselineMean")]; - calibfile.SetPtrAddress("BaselineMean", calib_baselinemean, calibfile.GetN("BaselineMean")); - calib_gainmean = new float[calibfile.GetN("GainMean")]; - calibfile.SetPtrAddress("GainMean", calib_gainmean, calibfile.GetN("GainMean")); - calib_triggeroffsetmean = new float[calibfile.GetN("TriggerOffsetMean")]; - calibfile.SetPtrAddress("TriggerOffsetMean", calib_triggeroffsetmean, calibfile.GetN("TriggerOffsetMean")); - calibfile.GetRow(0); - - //Set the event pointers -// std::cout << "...Setting event pointers..." << std::endl; - datafile.SetRefAddress("EventNum", event_id); - datafile.SetRefAddress("TriggerType", event_triggertype); - - event_data = new int16_t[data_ndata]; - datafile.SetPtrAddress("Data", event_data, data_ndata); - - event_offset = new int16_t[datafile.GetN("StartCellData")]; - datafile.SetPtrAddress("StartCellData", event_offset, datafile.GetN("StartCellData")); - - event_boardtimes = new int32_t[datafile.GetN("BoardTime")]; - datafile.SetPtrAddress("BoardTime", event_boardtimes, datafile.GetN("BoardTime")); - } - - ~CalFits() //Standard destructor - { - delete[] calib_baselinemean; - delete[] calib_gainmean; - delete[] calib_triggeroffsetmean; - delete[] event_data; - delete[] event_offset; - delete[] event_boardtimes; - } - - bool GetCalEvent() //Read calibrated event into the event variables - { - if(!npcaldata) { - ostringstream str; - str << "Pointer to the calibrated data not initialized!"; -#ifdef __EXCEPTIONS - throw runtime_error(str.str()); -#else - gLog << ___err___ << "ERROR - " << str.str() << endl; - return false; -#endif - } - if(datafile.GetNextRow() == false) { -// std::cout << "Last event reached..." << std::endl; - return false; - } - else { - UInt_t drs_calib_offset; - for(UInt_t pixel=0;pixel