source: fact/tools/pyscripts/sandbox/vogler/DRS_End_Noise_RMS_4.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
# ********************************
#
# modified and adapted py Patrick Vogler
#
# ################################
from ROOT import gSystem
gSystem.Load("calfactfits_h.so")   # according to new naming scheme
from ROOT import *


#define filenames

# 2011 12 05
calibfilename = '/fact/raw/2012/06/01/20120601_013.drs.fits.gz'    # NROI 300, pedestal

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


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


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


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)

data_average_run = np.zeros(numroi)
over_thresh = np.zeros(numpix)  # count for every pixel how many times first slice over 50mV

run_rms = np.zeros((numpix, numroi))   # rms



print "... looping..."

while caltest.GetCalEvent():    # Loop  ueber alle events


    if (caltest.event_triggertype == 4):  #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)

        for j in range (0, (numpix - 1)):
            if ((event[j,0] >50)and((j!=863)and(j!=297)and(j!=868))): #exclude crazy pixel (hard ID)
                over_thresh[j]=over_thresh[j]+1
                run_rms[j, : ]=run_rms[j, : ] + (event[j, : ])**2


for k in range (0, (numpix - 1)):
    run_rms[k, : ]=run_rms[k, : ]/over_thresh[k]
run_rms = np.sqrt(run_rms)

for i in range (0, (numpix - 1) ):
    data_average_run = data_average_run + run_rms[i]**2

data_average_run = data_average_run/numpix
data_average_run = np.sqrt(data_average_run)

maxindex = over_thresh.argmax()

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



#####################################################################################################################
print "creating plotters ..."

#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')

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')
myplotter3 = Plotter('titel of the plot', xlabel='time in slices', ylabel='amplitude calibrated data ... in mV')
myplotter4 = Plotter('titel of the plot', xlabel='time in slices', ylabel='amplitude calibrated data ... in mV')

over_th_plotter = Plotter('first slice ofer 50mv', xlabel='pixel number', ylabel='number of times first slice over 50mV')

print "Plotters created "
###################################################################################################################

myplotter(run_rms[20], 'pix 20 rms' )
myplotter2(run_rms[21], 'pix 21 rms' )
myplotter3(run_rms[maxindex], 'rms: pixel most times over 50mV')

print "Pix number (HW ID) most times over 50mV ", maxindex

over_th_plotter(over_thresh, 'run 20111205_022')


myplotter4(data_average_run , 'RMS (all pixel) over the whole run' )


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