Index: /fact/tools/pyscripts/sandbox/dneise/gapd_template_fit/EdgeOverlayTemplate_AllPixel_1440_0.csv =================================================================== --- /fact/tools/pyscripts/sandbox/dneise/gapd_template_fit/EdgeOverlayTemplate_AllPixel_1440_0.csv (revision 14253) +++ /fact/tools/pyscripts/sandbox/dneise/gapd_template_fit/EdgeOverlayTemplate_AllPixel_1440_0.csv (revision 14253) @@ -0,0 +1,310 @@ +### point-set of a single photon pulse template +### template determined with pulse overlay at: Edge +### Slice's Amplitude determined by calculating the +### value of maximum propability of slice -> AmplitudeMax +### mean of slice -> AmplitudeMean +### median of slice -> AmplitudeMedian +### for each slice +### Pixel number (CHid): 1440 +## +##time [slices],AmplitudeMax [mV],AmplitudeMean [mV],AmplitudeMedian [mV] +1,-1.25,-0.652001,-0.75 +2,-1.25,-1.07648,-1.25 +3,-1.25,-1.04914,-1.25 +4,-1.25,-0.823347,-1.25 +5,-1.25,-0.864903,-1.25 +6,-1.25,-1.05582,-1.25 +7,-1.25,-0.90244,-1.25 +8,-1.25,-0.903465,-1.25 +9,-1.25,-1.13651,-1.25 +10,-1.25,-1.127,-1.25 +11,-1.25,-0.736942,-0.75 +12,-1.25,-1.17842,-1.25 +13,-1.25,-1.16816,-1.25 +14,-1.25,-0.953418,-1.25 +15,-1.25,-0.971187,-1.25 +16,-1.25,-1.14135,-1.25 +17,-1.25,-0.966177,-1.25 +18,-1.25,-0.964521,-1.25 +19,-1.25,-1.17778,-1.25 +20,-1.25,-1.16627,-1.25 +21,-1.25,-0.759321,-0.75 +22,-1.25,-1.17925,-1.25 +23,-1.25,-1.14588,-1.25 +24,-1.25,-0.911606,-1.25 +25,-1.25,-0.930985,-1.25 +26,-1.25,-1.09513,-1.25 +27,-1.25,-0.928466,-1.25 +28,-1.25,-0.916932,-1.25 +29,-1.25,-1.13936,-1.25 +30,-1.25,-1.12685,-1.25 +31,-1.25,-0.749854,-0.75 +32,-1.25,-1.22809,-1.25 +33,-1.25,-1.24384,-1.25 +34,-1.25,-1.05792,-1.25 +35,-1.25,-1.09758,-1.25 +36,-1.25,-1.26678,-1.25 +37,-1.25,-1.09035,-1.25 +38,-1.25,-1.06052,-1.25 +39,-1.25,-1.26049,-1.25 +40,-1.25,-1.22482,-1.25 +41,-1.25,-0.784976,-0.75 +42,-1.25,-1.17956,-1.25 +43,-1.25,-1.11673,-1.25 +44,-1.25,-0.880103,-1.25 +45,-1.25,-0.887426,-1.25 +46,-1.25,-1.07458,-1.25 +47,-1.25,-0.914927,-1.25 +48,-1.25,-0.921556,-1.25 +49,-1.25,-1.15693,-1.25 +50,-1.25,-1.1531,-1.25 +51,-1.25,-0.754377,-0.75 +52,-1.25,-1.18882,-1.25 +53,-1.25,-1.18121,-1.25 +54,-1.25,-0.944888,-1.25 +55,-1.25,-0.928015,-1.25 +56,-1.25,-0.938964,-1.25 +57,-0.75,-0.439044,-0.75 +58,-0.25,0.114483,-0.25 +59,0.25,0.679467,0.75 +60,1.75,1.72504,1.75 +61,3.25,3.36502,3.25 +62,4.25,4.23151,4.25 +63,5.75,5.50953,5.75 +64,7.25,7.42987,7.25 +65,8.25,8.3014,8.25 +66,8.75,8.91587,8.75 +67,9.25,9.58884,9.75 +68,9.75,9.88877,9.75 +69,9.75,9.77901,9.75 +70,9.75,9.54593,9.75 +71,10.75,10.8013,10.75 +72,9.25,9.08734,9.25 +73,8.75,8.98959,8.75 +74,8.75,8.85727,8.75 +75,8.25,8.49497,8.25 +76,7.75,7.94158,7.75 +77,7.75,7.72934,7.75 +78,7.25,7.47874,7.25 +79,6.75,6.45184,6.25 +80,6.25,6.24493,6.25 +81,6.25,6.41068,6.25 +82,5.75,5.75264,5.75 +83,5.75,5.64393,5.75 +84,5.75,5.77257,5.75 +85,5.75,5.68198,5.75 +86,5.25,5.39562,5.25 +87,5.25,5.47681,5.25 +88,5.25,5.36306,5.25 +89,4.75,4.98915,4.75 +90,4.75,4.86087,4.75 +91,4.75,5.15261,5.25 +92,4.25,4.56613,4.25 +93,4.25,4.46258,4.25 +94,4.25,4.58308,4.25 +95,4.25,4.45715,4.25 +96,3.75,4.14818,3.75 +97,3.75,4.22064,3.75 +98,3.75,4.09372,3.75 +99,3.25,3.73819,3.75 +100,3.25,3.59881,3.25 +101,3.25,3.90122,3.75 +102,2.75,3.32751,3.25 +103,2.75,3.24553,2.75 +104,2.75,3.39886,3.25 +105,2.75,3.30552,2.75 +106,2.75,3.05532,2.75 +107,2.75,3.172,2.75 +108,2.25,3.13819,2.75 +109,2.25,2.8701,2.25 +110,2.25,2.83843,2.25 +111,2.25,3.22495,2.75 +112,2.25,2.72449,2.25 +113,1.75,2.65978,2.25 +114,1.75,2.81165,2.25 +115,1.75,2.70799,2.25 +116,1.75,2.44228,1.75 +117,1.75,2.53925,1.75 +118,1.75,2.4836,1.75 +119,1.25,2.18157,1.75 +120,1.25,2.1165,1.75 +121,1.75,2.47537,1.75 +122,1.25,1.97604,1.25 +123,1.25,1.93857,1.25 +124,1.25,2.11392,1.25 +125,1.25,2.05541,1.25 +126,0.75,1.81714,1.25 +127,0.75,1.94734,1.25 +128,0.75,1.91158,1.25 +129,0.75,1.6408,0.75 +130,0.75,1.61944,0.75 +131,0.75,2.0136,1.25 +132,0.75,1.54993,0.75 +133,0.75,1.53822,0.75 +134,0.75,1.74684,1.25 +135,0.75,1.69413,0.75 +136,0.25,1.47862,0.75 +137,0.75,1.61765,0.75 +138,0.25,1.59497,0.75 +139,0.25,1.31473,0.75 +140,0.25,1.26668,0.75 +141,0.25,1.63274,0.75 +142,0.25,1.13263,0.25 +143,0.25,1.10589,0.25 +144,0.25,1.30305,0.75 +145,0.25,1.24738,0.25 +146,-0.25,1.0436,0.25 +147,0.25,1.19519,0.25 +148,-0.25,1.19467,0.25 +149,-0.25,0.955759,0.25 +150,-0.25,0.946793,0.25 +151,0.25,1.35701,0.75 +152,-0.25,0.888685,0.25 +153,-0.25,0.88642,0.25 +154,-0.25,1.08269,0.25 +155,-0.25,1.04041,0.25 +156,-0.25,0.845374,0.25 +157,-0.25,1.00098,0.25 +158,-0.25,1.00575,0.25 +159,-0.25,0.763485,0.25 +160,-0.25,0.752686,-0.25 +161,-0.25,1.12951,0.25 +162,-0.75,0.666418,-0.25 +163,-0.75,0.664558,-0.25 +164,-0.25,0.844339,0.25 +165,-0.25,0.812723,0.25 +166,-0.75,0.601164,-0.25 +167,-0.75,0.761883,-0.25 +168,-0.75,0.748405,-0.25 +169,-0.75,0.516972,-0.25 +170,-0.75,0.515768,-0.25 +171,-0.25,0.911085,0.25 +172,-0.75,0.480652,-0.25 +173,-0.75,0.482902,-0.25 +174,-0.75,0.694839,-0.25 +175,-0.75,0.660046,-0.25 +176,-0.75,0.473287,-0.25 +177,-0.75,0.638261,-0.25 +178,-0.75,0.6346,-0.25 +179,-0.75,0.420355,-0.25 +180,-0.75,0.415532,-0.25 +181,-0.75,0.802615,0.25 +182,-0.75,0.351846,-0.25 +183,-0.75,0.352349,-0.25 +184,-0.75,0.55064,-0.25 +185,-0.75,0.507186,-0.25 +186,-0.75,0.327336,-0.25 +187,-0.75,0.472693,-0.25 +188,-0.75,0.480292,-0.25 +189,-0.75,0.254696,-0.75 +190,-0.75,0.260976,-0.75 +191,-0.75,0.667115,-0.25 +192,-1.25,0.228633,-0.75 +193,-0.75,0.251781,-0.75 +194,-0.75,0.439191,-0.25 +195,-0.75,0.422729,-0.25 +196,-0.75,0.233759,-0.75 +197,-0.75,0.400405,-0.25 +198,-0.75,0.41445,-0.25 +199,-1.25,0.189499,-0.75 +200,-1.25,0.203606,-0.75 +201,-0.75,0.591403,-0.25 +202,-1.25,0.161781,-0.75 +203,-1.25,0.170651,-0.75 +204,-0.75,0.384269,-0.25 +205,-0.75,0.35237,-0.25 +206,-1.25,0.170947,-0.75 +207,-0.75,0.343599,-0.25 +208,-0.75,0.335233,-0.25 +209,-1.25,0.116996,-0.75 +210,-1.25,0.120795,-0.75 +211,-0.75,0.512014,-0.25 +212,-1.25,0.0777615,-0.75 +213,-1.25,0.0914797,-0.75 +214,-0.75,0.304652,-0.25 +215,-1.25,0.2757,-0.75 +216,-1.25,0.0988412,-0.75 +217,-0.75,0.269365,-0.75 +218,-1.25,0.278669,-0.75 +219,-1.25,0.0757899,-0.75 +220,-1.25,0.0745676,-0.75 +221,-0.75,0.481154,-0.25 +222,-1.25,0.0453272,-0.75 +223,-1.25,0.0607856,-0.75 +224,-1.25,0.267645,-0.75 +225,-1.25,0.231283,-0.75 +226,-1.25,0.0648088,-0.75 +227,-1.25,0.223246,-0.75 +228,-1.25,0.239797,-0.75 +229,-1.25,0.0247974,-0.75 +230,-1.25,0.0365836,-0.75 +231,-0.75,0.439639,-0.25 +232,-1.25,0.00652546,-0.75 +233,-1.25,0.0306397,-0.75 +234,-1.25,0.23178,-0.75 +235,-1.25,0.189928,-0.75 +236,-1.25,0.0203009,-0.75 +237,-1.25,0.184404,-0.75 +238,-1.25,0.196727,-0.75 +239,-1.25,-0.0161804,-0.75 +240,-1.25,-0.00531809,-0.75 +241,-0.75,0.389288,-0.25 +242,-1.25,-0.0384582,-0.75 +243,-1.25,-0.019923,-0.75 +244,-1.25,0.186787,-0.75 +245,-1.25,0.173486,-0.75 +246,-1.25,-0.00191226,-0.75 +247,-1.25,0.17546,-0.75 +248,-1.25,0.181321,-0.75 +249,-1.25,-0.0250881,-0.75 +250,-1.25,-0.0269648,-0.75 +251,-0.75,0.373565,-0.25 +252,-1.25,-0.0513396,-0.75 +253,-1.25,-0.0371323,-0.75 +254,-1.25,0.175768,-0.75 +255,-1.25,0.15032,-0.75 +256,-1.25,-0.025751,-0.75 +257,-1.25,0.138683,-0.75 +258,-1.25,0.147341,-0.75 +259,-1.25,-0.0613598,-0.75 +260,-1.25,-0.0523961,-0.75 +261,-1.25,0.34831,-0.25 +262,-1.25,-0.0804417,-0.75 +263,-1.25,-0.0546038,-0.75 +264,-1.25,0.142936,-0.75 +265,-1.25,0.120325,-0.75 +266,-1.25,-0.0509319,-0.75 +267,-1.25,0.109485,-0.75 +268,-1.25,0.128321,-0.75 +269,-1.25,-0.0867021,-0.75 +270,-1.25,-0.0672052,-0.75 +271,-0.75,0.333344,-0.25 +272,-1.25,-0.0939269,-0.75 +273,-1.25,-0.0733059,-0.75 +274,-1.25,0.134929,-0.75 +275,-1.25,0.1132,-0.75 +276,-1.25,-0.0736547,-0.75 +277,-1.25,0.0993889,-0.75 +278,-1.25,0.109251,-0.75 +279,-1.25,-0.1051,-0.75 +280,-1.25,-0.0906483,-0.75 +281,-1.25,0.301502,-0.25 +282,-1.25,-0.115394,-0.75 +283,-1.25,-0.110289,-0.75 +284,-1.25,0.105354,-0.75 +285,-1.25,0.0790554,-0.75 +286,-1.25,-0.101445,-0.75 +287,-1.25,0.0875081,-0.75 +288,-1.25,0.101048,-0.75 +289,-1.25,-0.106498,-0.75 +290,-1.25,-0.097941,-0.75 +291,-1.25,0.298563,-0.25 +292,-1.25,-0.118258,-0.75 +293,-1.25,-0.111428,-0.75 +294,-1.25,0.10625,-0.75 +295,-1.25,0.0758478,-0.75 +296,-1.25,-0.0868433,-0.75 +297,-1.25,0.0774326,-0.75 +298,-1.25,0.0968001,-0.75 +299,-1.25,-0.10879,-0.75 +300,-1.25,-0.0958457,-0.75 Index: /fact/tools/pyscripts/sandbox/dneise/gapd_template_fit/gapdpulse.py =================================================================== --- /fact/tools/pyscripts/sandbox/dneise/gapd_template_fit/gapdpulse.py (revision 14253) +++ /fact/tools/pyscripts/sandbox/dneise/gapd_template_fit/gapdpulse.py (revision 14253) @@ -0,0 +1,176 @@ +#!/usr/bin/python -itt + +import numpy as np + +from scipy import optimize +from numpy import * + +class Parameter: + def __init__(self, value): + self.value = value + + def set(self, value): + self.value = value + + def __call__(self): + return self.value + +def fit(function, parameters, y, x = None): + def f(params): + i = 0 + for p in parameters: + p.set(params[i]) + i += 1 + return y - function(x) + + if x is None: x = arange(y.shape[0]) + p = [param() for param in parameters] + optimize.leastsq(f, p) + + +class gapdpulse( object ): + + bsl = -1.0 + height = 18.0 + start = 50. + rising = 1.2 + tau1 = 7. + tau2 = 25. + noise = 0.0001 + sampling = False + + def __call__(self, x): + if type(x) != type(np.empty(1)): + x = np.array(x) + val = np.zeros(x.shape) + zeros = np.zeros(x.shape) + + #shortcuts + bsl = self.bsl() + height = self.height() + start = self.start() + rising = self.rising() + tau1 = self.tau1() + tau2 = self.tau2() + + val += bsl + + + + # this is not really the maximum, but about it.... so I call it max + max = start + rising * tau1 + + e1 = height * (1 - np.exp(-(x-start)/tau1 ) ) + e2 = -1 * height * (1 - np.exp(-(x-max)/tau2 ) ) + val += np.where( x>start , e1, zeros) + val += np.where( x>max , e2, zeros) + + val += np.random.normal(0, self.noise, x.shape) + + #sampling + self.presampling = val.copy() + if self.sampling == True: + val = np.round(val*2)/2.0 + self.postsampling = val.copy() + return val + + + + + + + +# Target function, 6 parameters +def fitfunc( p, x): + + # give names to the parameters + bsl = p[0] + height = p[1] + start = p[2] + rising = p[3] + tau1 = p[4] + tau2 = p[5] + + # these are just helper variables + max = start + rising * tau1 + e1 = height * (1 - np.exp(-(x-start)/tau1 ) ) # the rising edge + e2 = -1 * height * (1 - np.exp(-(x-max)/tau2 ) )# the falling edge + + y = bsl + if x > start: + y += e1 + if x > max: + y += e2 + return y + +def fitfunc2(p, x): + # give names to the parameters + bsl = p[0] # just the baseline + height = p[1] # this is not the pulse height, but the maximum of the discharging edge + # it *would* be the pulse height only in case the + # recharging process would set in later... + + start = p[2] # start of avalanche aka start of discharge of diode capacitance + stop = p[3] # stop of avalanche, or start of recharging of diode + tau1 = p[4] # time constant of discharging process + tau2 = p[5] # time constant of recharging process + + # these are just helper variables + + e1 = height * (1 - np.exp(-(x-start)/tau1 ) ) + e2 = -1 * height * (1 - np.exp(-(x-stop)/tau2 ) ) + zeros = np.zeros(x.shape) + + y = np.zeros(x.shape) + y += bsl + y += np.where( x>start , e1, zeros) + y += np.where( x>stop , e2, zeros) + return y + + + +# Distance to the Targetfunction +def errfunc( p, x , y): + return fitfunc2(p, x) - y + + + +if __name__ == '__main__': + import matplotlib.pyplot as plt + plt.ion() + from gen_TH1D import SignalGeneratorCSV + + g = SignalGeneratorCSV('EdgeOverlayTemplate_AllPixel_1440_0.csv') + + x = np.linspace(0,300,300,False) + + fig = plt.figure() + plt.grid(True) + + plt.plot(x,g.csv_maxprob, label='max prob') + + pnames = ['bsl', 'height', 'start','stop','tau1','tau2'] + + p0 = [-1., 10., 50., 70, 30., 30.] # Initial guess for the parameters + + # plot function before fitting + #plt.plot( x, fitfunc2( p0, x) , '.:') + + p1, success = optimize.leastsq(errfunc, p0[:], args=(x[0:300], g.csv_maxprob[0:300])) + + # plot function after fitting + plt.plot( x, fitfunc2( p1, x) , '.:') + print 'Start values:' + for par,name in zip(p0,pnames): + print name, ':', par + print + print + print '------------------' + print 'Fit parameters' + for par,name in zip(p1,pnames): + print name, ':', par + print '------------------' + + print 'Fit successfull?:',bool(success) + + Index: /fact/tools/pyscripts/sandbox/dneise/gapd_template_fit/gen_TH1D.py =================================================================== --- /fact/tools/pyscripts/sandbox/dneise/gapd_template_fit/gen_TH1D.py (revision 14253) +++ /fact/tools/pyscripts/sandbox/dneise/gapd_template_fit/gen_TH1D.py (revision 14253) @@ -0,0 +1,425 @@ +#!/usr/bin/python -tt +# +# Dominik Neise, Werner Lustermann +# TU Dortmund, ETH Zurich +# +import numpy as np +import sys +from ROOT import * +from scipy import interpolate + +from plotters import Plotter +import matplotlib.pyplot as plt + +class SignalGenerator(object): + """ Signal Generator + generates signals for testing several helper classes like: + * fir filters + * signal extractors + """ + + def __init__(self, option_str = 'len 100 noise 3', name = 'SignalGenerator'): + """ initialize the generator + sets default signal to generate + """ + self.__module__ = 'generator' + self.option_str = option_str.lower() + self.options = make_options_from_str(option_str) + self.parse_options() + self.name = name + + def parse_options(self): + o = self.options #shortcut + if 'len' in o: + self.npoints = int(o['len'][0]) + else: + self.npoints = 100 + if 'noise' in o: + self.sigma = float(o['noise'][0]) + else: + self.sigma = 1 + if 'bsl' in o: + self.bsl = float(o['bsl'][0]) + else: + self.bsl = -0.5 + + if 'step' in o: + self.step_height = float(o['step'][0]) + self.step_start = int(o['step'][1]) + self.step_stop = int(o['step'][2]) + + if 'triangle' in o: + self.pulses = [] + # append 1st pulse to list of pulses + self.pulses.append( ( float(o['triangle'][0]) , float(o['triangle'][1]), int(o['triangle'][2]), int(o['triangle'][3]) ) ) + number_of_pulses_after_1st = (len(o['triangle'])-4)/2 + for i in range(number_of_pulses_after_1st): + self.pulses.append( ( float(o['triangle'][2*i+4]) , float(o['triangle'][2*i+5]), int(o['triangle'][2]), int(o['triangle'][3]) ) ) + + if 'spike' in o: + self.spikes = [] + for i in range(len(o['spike'])/2): + self.spikes.append( ( int(o['spike'][2*i]), float(o['spike'][2*i+1]) ) ) + + def __call__(self, option_str = ''): + if option_str: + self.option_str = option_str.lower() + self.options = make_options_from_str(self.option_str) + self.parse_options() + + signal = np.zeros(self.npoints) + signal += self.bsl + signal += np.random.randn(self.npoints) * self.sigma + if 'step' in self.options: + signal[self.step_start:self.step_stop] += self.step_height + if 'triangle' in self.options: + for pulse in self.pulses: + pos = pulse[0] + height = pulse[1] + rise = pulse[2] + fall = pulse[3] + start = pos - rise + stop = pos + fall + signal[start:pos] += np.linspace(0., height, rise) + signal[pos:stop] += np.linspace(height, 0. , fall) + if 'spike' in self.options: + for spike in self.spikes: + signal[spike[0]] += spike[1] + return signal + + def __str__(self): + s = self.name + '\n' + s += 'possible options and parameters\n' + s += ' * len: number of samples (100)\n' + s += ' * noise: sigma (1)\n' + s += ' * bsl: level (-0.5)\n' + s += ' * step: height, start, end\n' + s += ' * triangle: pos height risingedge, fallingedge [pos height ...]\n' + s += ' * spike: pos height [pos height ...]\n' + + s += 'current options are:\n' + for key in self.options.keys(): + s += key + ':' + str(self.options[key]) + '\n' + return s + + +class SignalGeneratorCSV(object): + + def __init__(self, file_name, option_str = 'len 100 noise 3', name = 'SignalGenerator'): + time, maxprob, mean, median = np.loadtxt( file_name, delimiter=',', unpack=True) + csv_data = median + + self.csv_maxprob = maxprob + self.csv_mean = mean + self.csv_median = median + + # shift the input data such, that the left baseline is equal to zero + csv_data = csv_data - csv_data[:50].mean() + + # I want the data to be longer than they are, so I add 1000 slices left and right + self.csv_data = np.zeros(2000 + len(csv_data)) + + for i in range(len(csv_data)): + self.csv_data[1000+i] = csv_data[i] + for i in range(1000): + self.csv_data[i] = csv_data[:50].mean() + self.csv_data[-1*i] = csv_data[-50:].mean() + + + x = np.arange(0,len(self.csv_data),1) + x2 = np.arange(0,len(self.csv_data),0.1) + + # fermi dirac distribution for weighting the input data. + # we do not believe the right tail is really higher than zero, so we make it go to zero. + weight_function_right = lambda x : 1. / (1 + np.exp((x - 1180.)/100.)) + weight_function_left = lambda x : 1. / (1 + np.exp((-1*(x - 1045.))/5.)) + + weights = np.array(map(weight_function_right, x)) + weights *= np.array(map(weight_function_left, x)) + weights *= 8 + weights [1200:] += 3 + weights += 2 + #y_weighted = weights * y + + #print len(weights) == len(x) + #print len(weights) == len(self.csv_data) + #print weights > 0 + + tck= interpolate.splrep(x, self.csv_data, w=weights) + self.spline_int = tck + + + x2 = np.arange(1000,1500,1) + t,c,k = tck + y2 = interpolate.splev(x2,(t,c,k),der=0) + + + #tck = interpolate.splrep(x, self.csv_data ,w=weights, s=50) + #xnew = np.arange(0,len(self.csv_data), 0.1) + #y = interpolate.splev(x,tck,der=0) + #y2 = interpolate.splev(x2,tck,der=0) + + + #plt.figure() + #plt.plot(x,self.csv_data,'.:b',x,weights,'m',x2,y2,'.:r') + #plt.show() + #raw_input('stop') + + # csv data was downshifted, I shift it up here + self.csv_data = csv_data - csv_data.min() + + + self.__module__ = 'CSV generator' + self.option_str = option_str.lower() + self.options = make_options_from_str(option_str) + self.parse_options() + self.name = name + + def parse_options(self): + o = self.options #shortcut + if 'len' in o: + self.npoints = int(o['len'][0]) + else: + self.npoints = 1024 + if 'noise' in o: + self.sigma = float(o['noise'][0]) + else: + self.sigma = 1 + if 'bsl' in o: + self.bsl = float(o['bsl'][0]) + else: + self.bsl = -0.5 + + if 'step' in o: + self.step_height = float(o['step'][0]) + self.step_start = int(o['step'][1]) + self.step_stop = int(o['step'][2]) + + if 'triangle' in o: + self.pulses = [] + # append 1st pulse to list of pulses + self.pulses.append( ( float(o['triangle'][0]) , float(o['triangle'][1]), int(o['triangle'][2]), int(o['triangle'][3]) ) ) + number_of_pulses_after_1st = (len(o['triangle'])-4)/2 + for i in range(number_of_pulses_after_1st): + self.pulses.append( ( float(o['triangle'][2*i+4]) , float(o['triangle'][2*i+5]), int(o['triangle'][2]), int(o['triangle'][3]) ) ) + + if 'spike' in o: + self.spikes = [] + for i in range(len(o['spike'])/2): + self.spikes.append( ( int(o['spike'][2*i]), float(o['spike'][2*i+1]) ) ) + + if 'csv' in o: + self.csvs = [] + for i in range(len(o['csv'])/2): + time = int( o['csv'][2*i] ) + amplitude = float( o['csv'][2*i+1] ) + self.csvs.append( (time, amplitude) ) + if 'rate' in o: + self.rate = float(o['rate'][0]) + + if 'signal' in o: + self.amplitude = float(o['signal'][0]) + self.position = int(o['signal'][1]) + + def __call__(self, option_str = ''): + if option_str: + self.option_str = option_str.lower() + self.options = make_options_from_str(self.option_str) + self.parse_options() + + signal = np.zeros(self.npoints) + signal += self.bsl + + # shortcut + o = self.options + + if 'step' in o: + signal[self.step_start:self.step_stop] += self.step_height + if 'triangle' in o: + for pulse in o: + pos = pulse[0] + height = pulse[1] + rise = pulse[2] + fall = pulse[3] + start = pos - rise + stop = pos + fall + signal[start:pos] += np.linspace(0., height, rise) + signal[pos:stop] += np.linspace(height, 0. , fall) + if 'spike' in o: + for spike in self.spikes: + signal[spike[0]] += spike[1] + + if 'csv' in o: + for csv in self.csvs: + amplitude = csv[1] + time = csv[0] + + #csv_data = self.csv_data.copy() + #scale + #csv_data *= amplitude + + x = np.arange(0,300,1) + t,c,k = self.spline_int + d = c.copy() + d *= amplitude + csv_data = interpolate.splev(x,(t,d,k),der=0) + + # add shifted + signal[time:time+len(csv_data)] += csv_data + + if 'rate' in o: + self._add_template_pulses(signal, self.rate) + + if 'signal' in o: + self._add_signal(signal) + + # add noise + if 'noise' in o: + signal += + np.random.normal(0.0,self.sigma, signal.shape) + + return signal + + def _add_signal(self, signal): + a = self.amplitude + p = self.position + + #csv_data = self.csv_data + #csv = csv_data.copy() * a + + x = np.arange(1000,1500,1) + t,c,k = self.spline_int + d = c.copy() + d *= a + csv = interpolate.splev(x,(t,d,k),der=0) + + npoints = self.npoints + 2 * len(csv) + internal_signal = np.zeros(npoints) + + internal_signal[p-70+len(csv):p+2*len(csv)-70] += csv + + signal += internal_signal[len(csv):-len(csv)] + + def _add_template_pulses(self, signal, rate): + + #csv_data = self.csv_data + x = np.arange(1000,1500,1) + t,c,k = self.spline_int + d = c.copy() + csv_data = interpolate.splev(x,(t,d,k),der=0) + + + period = 2.*1e9/(float(rate)*1e6) # in slices + # in order to simulate pulses which are just at the beginning and + # at the end of the signal, we pre- and append some points + npoints = self.npoints + 2 * len(csv_data) + + internal_signal = np.zeros(npoints) + + d=np.random.exponential(period) + pulse_positions = [] + while d 1: + num_events = int(sys.argv[1]) + if len(sys.argv) > 2: + phe = float(sys.argv[2]) + if len(sys.argv) > 3: + rate = float(sys.argv[3]) + if len(sys.argv) > 4: + noise = float(sys.argv[4]) + + + #less interesting + length = 1000 + baseline = 0.0 + signal_pos = 65 + + generator_string = 'len ' + str(length) + ' ' + generator_string += 'bsl ' + str(baseline) + ' ' + generator_string += 'signal ' + str(phe) + ' ' + str(signal_pos) + ' ' + generator_string += 'noise ' + str(noise) + ' ' + generator_string += 'rate ' + str(rate) + ' ' + gen = SignalGeneratorCSV('PulseTemplate_PointSet_0.csv', generator_string) + print gen + + rootfilename = str(phe)+'phe' + rootfilename += '_'+str(rate)+'MHz' + rootfilename += '_'+str(noise)+'mV' + rootfilename += '_'+str(num_events) + rootfilename += '.root' + + + + file = TFile(rootfilename, 'recreate') + + hist = TH1D('hist','to be done', length, -0.5,length-0.5) + hist.SetStats(kFALSE) + hist.SetXTitle('time in slices') + hist.SetYTitle('amplitude in mV - no baseline adjustment done yet :-( ') + + for event in range(10): + signal = gen() + + histogram_title = str(phe)+'phe @' + str(signal_pos) + histogram_title += ' and ' + str(rate) +'MHz dark count rate' + histogram_title += ', electronics noise with sigma=' + str(noise) +'mV' + histogram_title += ' #' + str(event) + hist.SetTitle(histogram_title ) + + for i,s in enumerate(signal): + hist.SetBinContent(i+1,s) + hist.Write() + hist.Reset() + file.Close()