source: fact/tools/pyscripts/sandbox/vogler/LP_template_1_time.py@ 14173

#!/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
# ********************************
#
# modified and adapted py Patrick Vogler
#
# ################################
from ROOT import gSystem
gSystem.Load("calfactfits_h.so")   # according to new naming scheme
from ROOT import *

import time

# c1 = TCanvas( 'c1', 'Example', 200, 10, 700, 500 )
# hpx = TH1F( 'hpx', 'px', 500,-50, 450 )



outfile = open(r'LP_template.txt','w')



#############
#
    # write to the file
#    outfile.write('The first number is: ' + str(num0) + '\n')
#    outfile.write('The second number is: ' + str(num1) + '\n')
#    outfile.write('The third number is: ' + str(num2) + '\n')
#    outfile.write('The last number is: ' + str(num3) + '\n')
#
#    outfile.write('The sum is: ' + str(num0 + num1 + num2 + num3) + '\n')
#    outfile.write('The product is: ' + str(num0 * num1 * num2 * num3) + '\n')

    # close the file

################







#define filenames

# 2012 04 17
#calibfilename = '/fact/raw/2012/04/17/20120417_003.drs.fits.gz'    # NROI 300, pedestal
#calibfilename = '/fact/raw/2012/04/17/20120417_033.drs.fits.gz'    # NROI 300, pedestal
calibfilename = '/fact/raw/2012/06/01/20120601_013.drs.fits.gz'    # NROI 300, pedestal

#datafilename = '/fact/raw/2012/04/17/20120417_015.fits.gz'         # NROI 300, LP_ext
#datafilename = '/fact/raw/2012/04/17/20120417_021.fits.gz'         # NROI 300, LP_ext
#datafilename = '/fact/raw/2012/04/17/20120417_036.fits.gz'         # NROI 300, LP_ext
#datafilename = '/fact/raw/2012/04/17/20120417_042.fits.gz'         # NROI 300, LP_ext

datafilename = '/fact/raw/2012/06/01/20120601_017.fits.gz'          # NROI 300, 5 minutes physics data with
                                                                    # interleaved LP_ext

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

from plotters import Plotter         # ADC display
#from plotters import CamPlotter      # event display
#from drs_spikes import DRSSpikes


print "Testing object creation: "
caltest = CalFactFits(datafilename, calibfilename)
npcalevent  = np.empty( caltest.npix * caltest.nroi, np.float64) #.reshape(caltest.npix ,caltest.nroi)
caltest.SetNpcaldataPtr(npcalevent)


numroi = np.int64(caltest.nroi)
numpix = np.int64(caltest.npix)
numevents = np.int64(caltest.nevents)

num_LP_ev = 0   # number of ext Lightpulser events
num_ped = 0     # number of pedestal events


print "Common variables run information: "
print "ROI: ", numroi
print "#Pix: ", numpix
print "Number of events: ", numevents
print


event = np.zeros((numpix, numroi))  # create an array to store an event in the format numpix * numroi (2-dim array)

run_average = np.zeros((numpix, numroi))   # create an array to store the "average event" of a run

data_average_run = np.zeros(numroi)



time_start = time.clock()
print "start loop", time_start


while caltest.GetCalEvent():    # Loop  ueber alle events
    # print npcalevent  ## Daten, Array   1 dim Laenge caltest.npix * caltest.nroi   1 Event


    if ((caltest.event_triggertype > 256) and (caltest.event_triggertype < 512)):  #ext LP event
        print 'event id:', caltest.event_id, '  Trigger Type:', caltest.event_triggertype
        num_LP_ev = num_LP_ev + 1
        event = np.reshape(npcalevent, (numpix, -1))     # bring the event the shape numpix * numroi (2-dim array)
        run_average += event

   # if(num_LP_ev==50):    # just for testing, not to read the whole file
    #    break

time_stop = time.clock()
print "loop finished", time_stop

print "Looped ... "

for i in range (0, (numpix - 1) ):

    # if (()and()and()and()and()and())

    data_average_run += run_average[i]


print "Looped second loop "


del caltest
print "caltest deleted "



print "Common variables run information: "
print "ROI: ", numroi
print "#Pix: ", numpix
print "Number of events: ", numevents
print "Number of external Lightpulser events: ", num_LP_ev
print


#####################################################################################################################
#make a Plotter class ... this is an easy way for plotting ... but there are
# many was to plot data ... without this class ... it was written for convenience, but there is no strong reason to use it...
#myplotter = Plotter('titel of the plot', xlabel='time in slices', ylabel='amplitude calibrated data ... in mV')

myplotter2 = Plotter('titel of the plot', xlabel='time in slices', ylabel='amplitude calibrated data ... in mV')

###################################################################################################################

myplotter2( data_average_run / (numpix * num_LP_ev) , 'average signal (all pixel) over the whole run' )



outfile.write('This file contains a template of the Lightpulser pulse, i.e. an average of one file. \n')
outfile.write('data file: ' + str(datafilename) + '\n')
outfile.write('DRS calib file: ' + str(calibfilename) + '\n')
outfile.write('Number of Lightpulser events in this run: ' + str(num_LP_ev) + '\n')
outfile.write('ROI: ' + str(numroi) + '\n')
outfile.write('\n')
outfile.write('slice   signal [mV] \n')
outfile.write('\n')


for k in range(0,(numroi)):
    outfile.write( str(k) + '   ' + str((data_average_run[k])/(10 * num_LP_ev)) +  '\n')

outfile.close

print
print "The loop over the whole file took",(time_stop-time_start),"seconds to complete."
print "processing rate: ",(numevents/(time_stop-time_start)),"events per second"


answer = raw_input('type "quit" to quit ')
while not 'quit' in answer:
     answer = raw_input('type "quit" to quit ')