Index: /fact/tools/pyscripts/examples/CalFitsPerformance.py =================================================================== --- /fact/tools/pyscripts/examples/CalFitsPerformance.py (revision 13503) +++ /fact/tools/pyscripts/examples/CalFitsPerformance.py (revision 13503) @@ -0,0 +1,78 @@ +#!/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/examples/CalFitsPerformanceWeave.py =================================================================== --- /fact/tools/pyscripts/examples/CalFitsPerformanceWeave.py (revision 13503) +++ /fact/tools/pyscripts/examples/CalFitsPerformanceWeave.py (revision 13503) @@ -0,0 +1,114 @@ +#!/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/examples/CalFitsTest.py =================================================================== --- /fact/tools/pyscripts/examples/CalFitsTest.py (revision 13503) +++ /fact/tools/pyscripts/examples/CalFitsTest.py (revision 13503) @@ -0,0 +1,58 @@ +#!/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/pyfact/calfits.h =================================================================== --- /fact/tools/pyscripts/pyfact/calfits.h (revision 13503) +++ /fact/tools/pyscripts/pyfact/calfits.h (revision 13503) @@ -0,0 +1,188 @@ +//******************************** +// +// 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