<file_sep>#! /usr/bin/env python3

from __future__ import division
from __future__ import print_function
ls

import os, sys, inspect

cmd_subfolder = os.path.realpath(os.path.abspath(os.path.join(os.path.split(inspect.getfile(inspect.currentframe()))[0],"src")))
if cmd_subfolder not in sys.path:
  sys.path.insert(0, cmd_subfolder)

from TLeCroy import *

# --------------------------------
# constants to be set by the user
# --------------------------------
channelNr = 4

traceNr = 18


 
#inName = "--testpulse--"
#inName = "--testpulse--autosave--"
inName = "--testpulse--autosave--segments--"


#inDir = "/Users/Lennard/Bachelor_Thesis_LF/BA_Lennard_Franz/cs137-longzylinder/"

inDir = "cs137-longzylinder/"


#inDir = "/Users/Lennard/desktop/osci/lennard/segment-flat-CeBr-CS137/"
#inDir = "/Users/Lennard/desktop/osci/lennard/segment-CeBr-cs137/"
#inDir = "/Users/Lennard/desktop/osci/lennard/cebr-big-4ch-cs137/"
#inDir = "/volumes/lennad/Lennard Franz/cs137-longzylinder-shaper/"
#inDir = "/Users/Lennard/desktop/osci/lennard/cs137-longzylinder/"
#inDir = "/Users/Lennard/desktop/osci/lennard/cs137-filter/"
#inDir = "/Users/Lennard/desktop/autosave-segments/"
#inDir = "/Users/Lennard/Desktop/osci/lennard/cebr-cs137/"

#inDir = "/home/lfranz/work/osci/lennard/segment-longcable-CeBr-CS137/"
#inDir = "/home/lfranz/work/osci/lennard/segment-flat-plexi-background/"
#inDir = "/home/lfranz/work/osci/lennard/segment-CeBr-ra226/"
#inDir = "/home/lfranz/work/osci/lennard/segment/"
#inDir = "/home/lfranz/work/autosave-segments/"
#inDir = "/ZIH.fast/projects/d-lab_data/MAPMT/2019_gsi/segment/"
#inDir = "/ZIH.fast/projects/d-lab_data/MAPMT/2019_gsi/segments-new/"
#inDir = "/ZIH.fast/projects/d-lab_data/MAPMT/2019_gsi/segments-new/Autosave-new/"
# inDir = "/ZIH.fast/projects/d-lab_data/MAPMT/2019_gsi/Autosave/"
# --------------------------------

fileName = inDir + "C" + str(channelNr).zfill(1) + inName + str(traceNr).zfill(5) + ".trc"
TLC = TLeCroy(fileName, debug=True)
TLC.PrintPrivate()
#TLC.PlotTrace(raw=False)

header = TLC.GetHeader()
x, y = TLC.GetTrace()
y = y - header['VERTICAL_OFFSET']
#y = -y
nano = 1e-9

import numpy as np
import scipy.integrate as integrate
 
#print(header['VERTICAL_OFFSET'])
#print('   vertical gain:', header['VERTICAL_GAIN'])
#print('----------------')
seqLength = header['WAVE_ARRAY_COUNT'] // header['SUBARRAY_COUNT']
seqFirst = header['FIRST_VALID_PNT']
print('  LENGTH', seqLength)
print(len(x)/1000,len(y)/1000)


import matplotlib.pyplot as plt
#plt.plot(x / nano, y)


idxHigh = np.zeros(header['SUBARRAY_COUNT'],dtype=np.int)
idxLow = np.zeros(header['SUBARRAY_COUNT'],dtype=np.int)
mean = np.zeros(header['SUBARRAY_COUNT'])
std = np.zeros(header['SUBARRAY_COUNT'])
area =  np.zeros(header['SUBARRAY_COUNT'])
area_all = np.zeros(header['SUBARRAY_COUNT'])
A = np.zeros(header['SUBARRAY_COUNT'])
area_bool = np.zeros(header['SUBARRAY_COUNT'])
area_mean = np.zeros(header['SUBARRAY_COUNT'])
area_all = np.zeros(header['SUBARRAY_COUNT'])

for seqNr in range(1, header['SUBARRAY_COUNT'] + 1):
  #bondarys of sequences
  idxHigh[seqNr-1] = seqFirst +( seqLength * seqNr - 1)
  idxLow[seqNr-1] = idxHigh[seqNr-1] - seqLength + 1
  #finding basline
  #y[...] has to be adjusted depending vertical position of the puls                     
  from scipy.stats import norm
  mean[seqNr-1], std[seqNr-1] = norm.fit(y[ idxLow[seqNr-1] : idxLow[seqNr-1] + int(0.35*seqLength)])

  #print(seqNr, idxLow, idxHigh)
mean_all = sum(mean)/header['SUBARRAY_COUNT']
std_all = sum(std)/header['SUBARRAY_COUNT'] 
#print(mean_all,std_all)
#plt.hist(mean, bins=100)


#plt.hist(y[idxLow[seqNr-1] : idxLow[seqNr-1]+int(0.35*seqLength)])


for seqNr in range(1, header['SUBARRAY_COUNT'] + 1):
 # if mean[seqNr -1] > 0.0285 and mean[seqNr - 1] < 0.0315 and np.min( y[idxLow[seqNr-1] : idxHigh[seqNr-1]]-mean[seqNr-1]) > -0.0035:
    
    #area of the pulse
  area[seqNr-1] = integrate.trapz((y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]- mean[seqNr-1]) > 2*std[seqNr-1] ,
                  x[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]/nano)  
 
    # area_all[seqNr-1] = integrate.trapz((y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]] - mean_all) > 2*std_all ,
                     # x[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]/nano)
    
  #boolArr = np.where((y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]- mean[seqNr-1]) > (3*std[seqNr-1]))
  boolArr = (y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]- mean[seqNr-1]) > (3*std[seqNr-1])
  #print((boolArr))
  #area_bool[seqNr-1] = np.sum(y[boolArr[0]])*(header['HORIZ_INTERVAL']/nano) - (mean[seqNr-1])*(len(boolArr[0])*header['HORIZ_INTERVAL']/nano)  
  #area_bool[seqNr-1] = np.sum(y[boolArr[0]])
  area_bool[seqNr-1] = np.sum((y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]*(header['HORIZ_INTERVAL']/nano)),where = boolArr) - mean[seqNr-1]*np.sum(boolArr)*header['HORIZ_INTERVAL']/nano


  #print(np.trapz(y[boolArr[0]],x[boolArr[0]])/nano,area_bool[seqNr-1])
  #print(np.sum(y[idxLow[seqNr-1]:idxHigh[seqNr-1]]>(mean[seqNr-1]+3*std[seqNr-1])),len(boolArr[0])) 
  area_bool_positive = area_bool > 0.0 
  
  
                    
  #print(seqNr, idxLow, idxHigh)
  #area = integrate.trapz(y[idxLow:idxHigh],x[idxLow:idxHigh]/nano)
  #print(area)
  if seqNr > 18  and seqNr < 20 :
    #plt.plot(x[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]  / nano, y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]-mean[seqNr-1], linewidth=0.5, color='k')
    plt.plot(x[idxLow[0]:idxHigh[0]]  / nano, (y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]-mean[seqNr-1])*1000, linewidth=1, color='k')
  
  #area_mean[seqNr-1]= mean[seqNr-1]*len(boolArr[0])*header['HORIZ_INTERVAL']/nano
  area_all[seqNr-1] = np.sum((y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]*(header['HORIZ_INTERVAL']/nano)),where = boolArr)
  #print(area_all[seqNr-1])


#plt.hist(area_mean)
#plt.hist(area_all)
#print(np.sum(y[idxLow[seqNr-1]:idxHigh[seqNr-1]]>(mean[seqNr-1]+3*std[seqNr-1])))
#print(area_bool)


plt.xlabel('time (ns)')
plt.ylabel('voltage (mV)')
plt.grid()
#plt.plot(x[idxLow[428]:idxHigh[428]]/nano , y[idxLow[428]:idxHigh[428]]-mean[428])

#x=np.linspace(0,1000.1e-9,num=10001,endpoint=False)
#print(x)
#plt.plot(x , y[idxLow[3]:idxHigh[3]]-mean[3])
#area_1 = integrate.trapz((y[idxLow[1]:idxHigh[1]] - mean_all) - 1*std_all,  x[idxLow[1]:idxHigh[1]]/nano)
    
            
#print(area_1)

#print("mean all:",mean_all,"std all:", std_all)

#print(mean,std) 

#print("-----------------------")
#print("area corr" , area)
#print("std" , std)
#print("mean" , mean)
#print("area uncorr" , #area_1)
#print("-----------------------")

plt.xlim(-50,150)



#ax1.hist(area,bins=50)
mean_E,std_E = norm.fit(area) 
#print("---------------------")
#print("mean:", mean_E, "std:", std_E)
#print("---------------------")

#plt.hist(area_bool ,bins=50)
#mean_E_all, std_E_all = norm.fit(area_bool)
#print("---------------------")
#print("mean:", mean_E_all, "std:", std_E_all)
#print("---------------------")


#import tikzplotlib
#tikzplotlib.save("pulse_mit_shaper.tex")



plt.show()

<file_sep>import numpy as np
import matplotlib.pyplot as plt

#inDir = "GSI/bi207-1660mV.txt"
#inDir = "GSI/cs137-1660mV-neu.txt"
inDir = "GSI/co60-1660mV.txt"

#inDir = "test-frequenzgenerator-pos-pulse.txt"
#inDir = "GSI/cs137-noshaper-1400mV.txt"


data = np.loadtxt(inDir)   
x = data[:,0]
y = data[:,3]
    
  
#a = 32.63 
#b = 0.169

#x = a*np.exp(b*x)
  
print(x)
  

N_bins = 500
bin_width = int(10000/N_bins)



idxLow = np.linspace(0,9999 - (bin_width-1), N_bins, dtype = np.int)
idxHigh = np.linspace((bin_width - 1),9999, N_bins, dtype = np.int)
idxMid = np.linspace((bin_width/2 - 1),9999 - (bin_width/2) , N_bins, dtype = np.int)



bins = np.zeros(N_bins)
n = np.zeros(N_bins) 

#print(idxLow)
#print(idxMid,len(idxMid))
#print(idxHigh)

for i in range(0,N_bins,1):
    #print(i)
    bins[i]= x[idxMid[i]]
    n[i] = np.sum(y[idxLow[i]:idxHigh[i]])
    

#print(n)
#print(bins)
    
bins_new = np.zeros(10000)
n_new = np.zeros(10000)
    


for i in range(0,N_bins,1):
    bins_new[idxLow[i] : idxHigh[i]] = x[idxMid[i]]
    n_new[idxLow[i]:idxHigh[i]] = np.sum(y[idxLow[i]:idxHigh[i]])

#print(bins_new)
#print(n_new)

a = 32.63 
b = 0.169

calib = a*np.exp(b*bins)

calib_bins = calib[calib<1600]

N = n[calib<1600]

#plt.plot(x,y)

n,bins,patches=plt.hist(bins, bins, weights=n, histtype = 'step', color='k')

#n,bins,patches=plt.hist(calib_bins, bins, weights=N, histtype = 'step', color='k')

from scipy.signal import find_peaks
peaks,properties = find_peaks(n,width=3,height=600, rel_height=0.7)
print("index of peaks:", peaks)
print("X Values of peaks:", x[peaks] )
print(properties)

from scipy.optimize import curve_fit

def gauss (x,a,sigma,mu):
    return a*np.exp(-(x-mu)**2/(2*sigma**2)) 



for peakNr in range(1, len(peaks) ):
    
    #peakNr = 1 
    fit_width = 0.5
    
    x_values = bins[np.int(peaks[peakNr] - (fit_width) * properties["widths"][peakNr]) :np.int( peaks[peakNr]+ (fit_width+0.15) * (properties["widths"][peakNr])) ]
    y_values = n[np.int(peaks[peakNr] - (fit_width) * properties["widths"][peakNr]) :np.int( peaks[peakNr]+ (fit_width+0.15) * (properties["widths"][peakNr]))  ]
  
    a_bound = [0, 1.1*np.max(y_values)]
    sigma_bound = [0, 1]
    mu_bound = [np.min(x_values), np.max(x_values)]


    popt, pcov = curve_fit(gauss, x_values , y_values, bounds= ([a_bound[0],sigma_bound[0],mu_bound[0]], [a_bound[1],sigma_bound[1],mu_bound[1]] ) ) 


    #plt.plot(x, gauss(x, *popt), label='gauss Fit', color='g')
    #plt.plot((popt[2], popt[2]), (0, 1500), 'r-')
    
    #statistic error of curve fit 
    perr = np.sqrt(np.diag(pcov))
  
    #FWHM
    FWHM = 2* np.sqrt(2*np.log(2)) * popt[1]
    delta_FWHM = 2* np.sqrt(2*np.log(2)) * perr[1]

  
    print("Peak Nr:" , peakNr) 
    print("Fit Bereich:" , mu_bound)
    print("Energy:" , popt[2], "statistic error of energy:" , perr[2], "FWHM:", FWHM, "statistic error of FWHM", delta_FWHM )
  

for peakNr in range(1, len(peaks) ):
    
    #peakNr = 0 
    fit_width = 0.1
    
    x_values = bins[np.int(peaks[peakNr] - (fit_width) * properties["widths"][peakNr]) :np.int( peaks[peakNr]+ (fit_width+0.2) * (properties["widths"][peakNr])) ]
    y_values = n[np.int(peaks[peakNr] - (fit_width) * properties["widths"][peakNr]) :np.int( peaks[peakNr]+ (fit_width+0.2) * (properties["widths"][peakNr]))  ]
  
    a_bound = [0, 1.1*np.max(y_values)]
    sigma_bound = [0, 1]
    mu_bound = [np.min(x_values), np.max(x_values)]


    popt, pcov = curve_fit(gauss, x_values , y_values, bounds= ([a_bound[0],sigma_bound[0],mu_bound[0]], [a_bound[1],sigma_bound[1],mu_bound[1]] ) ) 


    plt.plot(x, gauss(x, *popt), label='gauss Fit', color='g')
    plt.plot((popt[2], popt[2]), (0, 1500), 'r-')
    
    #statistic error of curve fit 
    perr = np.sqrt(np.diag(pcov))
  
    #FWHM
    FWHM = 2* np.sqrt(2*np.log(2)) * popt[1]
    delta_FWHM = 2* np.sqrt(2*np.log(2)) * perr[1]

  
    print("Peak Nr:" , peakNr) 
    print("Fit Bereich:" , mu_bound)
    print("Energy:" , popt[2], "statistic error of energy:" , perr[2], "FWHM:", FWHM, "statistic error of FWHM", delta_FWHM )
  





xticks = np.arange(0,110,10)
plt.xticks(xticks)

plt.xlim(0,100)
#plt.legend()
plt.grid()
plt.title("Histogramm der ToT")
plt.xlabel("ToT (ns)")
plt.ylabel("Anzahl N")

plt.xlim(0,30)
#plt.ylim(0,1250)

#import tikzplotlib
#tikzplotlib.save("GSI_generator.tex")




plt.show()




# FÜr den fi der mit dem Frequenzgernator erstellten Histogramme ergibt sich eine
# FWHM von 0.25 bis 0.35 bei 1 ns ungenauigkeit des genrator 
# Trailing und falling von je  2,5ns beachten 
#
#
#
#
#
#
#
#
#








<file_sep>
import iotool as io

class TLeCroy:

  def __init__(self, fname, debug=None, raw=None):
    self._FileName = fname
    if raw:
      self._Raw = raw
    else:
      self._Raw = False
    self._Info = dict() 
    self._Header = dict() 
    self._x, self._y, self._Header, self._Info = io.readTrc(self._FileName, raw = self._Raw)
    if debug:
      self._Debug = debug
    else:
      self._Debug = False
  

  def PrintPrivate(self):
    print(" *** INFO::lecroy:PrintPrivate ***")
    print("  self._FileName        " + str(self._FileName))
    print("  self._Debug           " + str(self._Debug))
    print("  self._Raw             " + str(self._Raw))
    print("  ---")
    print("   self._Header")
    for item in self._Header.items():
      print(item)
    print("  ---")
    print("  self._Info")
    for item in self._Info.items():
      print(item)
    print(" ------------------------")
    if self._Debug:
      print(type(self._y[0]))


  def GetTrace(self):
    return self._x, self._y


  def GetY(self):
    return self._y


  def GetHeader(self):
    return self._Header


  def GetInfo(self):
    return self._Info


  def PlotTrace(self, raw=True):
    import matplotlib.pyplot as plt
    if raw:
      plt.plot(self._x, self._y, 'k', label=self._FileName)
      plt.ylabel('amplitude (a.u.)')
      plt.xlabel('time (a.u.)')
    else:
      milli = 1/1000.0
      plt.plot(self._x, self._y / milli, 'k', label=self._FileName)
      plt.ylabel('volts (mV)')
      plt.xlabel('time (s)')

    plt.title(self._FileName)
    plt.show()


<file_sep>#! /usr/bin/env python3

from __future__ import division
from __future__ import print_function

import os, sys, inspect

cmd_subfolder = os.path.realpath(os.path.abspath(os.path.join(os.path.split(inspect.getfile(inspect.currentframe()))[0],"src")))
if cmd_subfolder not in sys.path:
  sys.path.insert(0, cmd_subfolder)

from TLeCroy import *   # lecroy.py 
import matplotlib.pyplot as plt
from scipy.stats import norm
import numpy as np
import scipy.integrate as integrate
import matplotlib.mlab as mlab
from scipy.signal import find_peaks, peak_widths, peak_prominences 

# --------------------------------
# constants to be set by the user
# --------------------------------
channelNr = 1

 
StartTrace = 0
StopTrace =20



#inName = "--testpulse--"
#inName = "--testpulse--autosave--"
inName = "--testpulse--autosave--segments--"

#inDir = "/volumes/lennad/segment/segment-CeBr-cs137/"
#inDir = "/volumes/lennad/<NAME>/th-zylinder/"
inDir = "/volumes/lennad/<NAME>/cs137-longzylinder-shaper/"
#inDir = "/Users/Lennard/desktop/osci/lennard/segment-flat-CeBr-CS137/"
#inDir = "/Users/Lennard/desktop/osci/lennard/cs137-longzylinder-shaper/"
#inDir = "/Users/Lennard/destop/osci/lennard/cebr-big-4ch-cs137/"
#inDir = "/Users/Lennard/desktop/osci/lennard/segment-CeBr-ra226/"
#inDir = "/Users/Lennard/desktop/osci/lennard/ra226-filter/"
#inDir = "/Users/Lennard/desktop/autosave-segments/"
#inDir = "/Users/Lennard/desktop/osci/lennard/gag-am241/"


#inDir = "/home/lfranz/work/osci/lennard/segment-longcable-CeBr-CS137/"
#inDir = "/home/lfranz/work/osci/lennard/segment-flat-CeBr-CS137/"
#inDir = "/home/lfranz/work/osci/lennard/segment-flat-plexi-background/"
#inDir = "/home/lfranz/work/osci/lennard/segment-CeBr-am241/"
#inDir = "/home/lfranz/work/osci/lennard/segment-CeBr-cs137/"
#inDir = "/home/lfranz/work/osci/lennard/segment-CeBr-ra226/"

#inDir = "/home/lfranz/work/osci/lennard/segment/"
#inDir = "/home/lfranz/work/autosave-segments/"
#inDir = "/ZIH.fast/projects/d-lab_data/MAPMT/2019_gsi/segment/"
#inDir = "/ZIH.fast/projects/d-lab_data/MAPMT/2019_gsi/segments-new/"
#inDir = "/ZIH.fast/projects/d-lab_data/MAPMT/2019_gsi/segments-new/Autosave-new/"
# inDir = "/ZIH.fast/projects/d-lab_data/MAPMT/2019_gsi/Autosave/"
# --------------------------------

import numpy as np
import scipy.integrate as integrate





#empty array for area of sequnces for every trace
A = np.empty([0])
print(inDir)
for traceNr in range(StartTrace, StopTrace + 1, 1):
  print(" INFO: working on trace: ", "C"+ str(channelNr).zfill(1) + inName + str(traceNr).zfill(5) + ".trc")
 
  fileName = inDir + "C" + str(channelNr).zfill(1) + inName + str(traceNr).zfill(5) + ".trc"
  TLC = TLeCroy(fileName, debug=True)
  #TLC.PrintPrivate()
  # TLC.PlotTrace(raw=False)

  header = TLC.GetHeader()
  x, y = TLC.GetTrace()
  y = y - header['VERTICAL_OFFSET']
  y = -y 
  nano = 1e-9

  #print('----------------')
  #print('   vertical gain:', header['VERTICAL_GAIN'])
  #print('----------------')
  seqLength = header['WAVE_ARRAY_COUNT'] // header['SUBARRAY_COUNT']
  seqFirst = header['FIRST_VALID_PNT']
  #print('  LENGTH', seqLength)

  #plt.plot(x / nano, y)
  
  #arrays  
  idxHigh = np.zeros(header['SUBARRAY_COUNT'],dtype=np.int)
  idxLow = np.zeros(header['SUBARRAY_COUNT'],dtype=np.int)
  mean = np.zeros(header['SUBARRAY_COUNT'])
  std = np.zeros(header['SUBARRAY_COUNT'])
  area =  np.zeros(header['SUBARRAY_COUNT'])
  area_all = np.zeros(header['SUBARRAY_COUNT'])
  area_all_sum = np.zeros(header['SUBARRAY_COUNT'])
  area_bool =np.zeros(header['SUBARRAY_COUNT'])
  tot = np.zeros(header['SUBARRAY_COUNT'])
  A = np.zeros(header['SUBARRAY_COUNT'])

# determine background
  for seqNr in range(1, header['SUBARRAY_COUNT'] + 1):
    #bondarys of sequences
    idxHigh[seqNr-1] = seqFirst +( seqLength * seqNr - 1)
    idxLow[seqNr-1] = idxHigh[seqNr-1] - seqLength + 1
    #finding basline
    #y[...] has to be adjusted depending vertical position of the puls                     
    mean[seqNr-1], std[seqNr-1] = norm.fit(y[ idxLow[seqNr-1] : idxLow[seqNr-1] + int(0.35*seqLength)])

  mean_all = sum(mean)/header['SUBARRAY_COUNT']
  std_all = sum(std)/header['SUBARRAY_COUNT'] 
  #Histogram of values for baseline 
  #plt.hist(y[ idxLow[seqNr-1] : idxLow[seqNr-1] + int(0.35*seqLength)])
  

#determination of area of pulses
  for seqNr in range(1, header['SUBARRAY_COUNT'] + 1):
   # if mean[seqNr -1] > 0.0285 and mean[seqNr - 1] < 0.0315 and np.min( y[idxLow[seqNr-1] : idxHigh[seqNr-1]]-mean[seqNr-1]) > -0.0035:

#area of the pulse
    #area[seqNr-1] = np.trapz((y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]- mean[seqNr-1]) > 3*std[seqNr-1] ,
                 #x[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]/nano)  
 
    #area_all[seqNr-1] = np.sum(y[idxLow[seqNr - 1] : idxHigh[seqNr - 1]])/(nano/header['HORIZ_INTERVAL']) - ((mean_all)*seqLength*header['HORIZ_INTERVAL']/nano)
   
    #print(area_all_sum)
    #print(area_all)
    #condition to differentiate values of pulse from background 
    #boolArr = np.where((y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]- mean[seqNr - 1]) > (3*std[seqNr -1]))
    boolArr = (y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]- mean[seqNr-1]) > (3*std[seqNr-1])

    #sum of all values which meet condition and scaling 
    #substracting of area belonging to background
    
    #print(y[boolArr[0]] - mean[seqNr-1])
    #area_bool[seqNr-1] = np.sum(y[boolArr[0]])*(header['HORIZ_INTERVAL']/nano) - (mean[seqNr-1])*(len(boolArr[0])*header['HORIZ_INTERVAL']/nano)
    
    A[seqNr-1]= np.max(y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]])
    
    area_bool[seqNr-1] = np.sum((y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]*(header['HORIZ_INTERVAL']/nano)),where = boolArr) - mean[seqNr-1]*np.sum(boolArr)*header['HORIZ_INTERVAL']/nano
    tot[seqNr-1] = (np.sum(boolArr)-1)*header['HORIZ_INTERVAL']/nano

    #print(area_bool[seqNr - 1],tot[seqNr - 1] )
#plt.plot(tot,area_bool, 'bo', ms=1) 


#print(A)

#E_real = np.array([569, 662,1173,1332,1063])
#E_fit = np.array([6.485, 7.2341,10.7222,11.7675,10.2165])


E_real = np.array ([662,1173,1332,569,1063])
E_fit = np.array ([0.7535521180502402,1.1657715187599953,1.273437535904615,0.701803960984309,1.1086378999695787])


from scipy import stats 
slope, intercept, r_value, p_value, std_err = stats.linregress(E_fit, E_real)

print(slope)
print(intercept)

E = intercept+slope*area_bool

print(E)


plt.plot(tot,area_bool,'bo',ms=1.5,color='k', label='Daten' )



#plt.hist2d(tot,area_bool, bins= 50)
#plt.yscale('log')
#plt.xscale('log')

from scipy.optimize import curve_fit

def fit(x,a,b):
    return a*np.exp(b*x)

popt, pcov = curve_fit(fit, tot ,area_bool, bounds= ([0,0],[1000,1]),)

perr = np.sqrt(np.diag(pcov))


print(popt)

print(perr)


x_fit = np.linspace(50,106,1000)


plt.plot(x_fit, fit(x_fit, *popt), label='exp. Fit', color='r')


plt.text(52,0.75, '$A(T)=a \cdot e^{(b \cdot x)}$')
plt.text(52,0.70, '$a= 0.053 \ {nVs}$')
plt.text(52,0.65, '$b= 0.027 \ ns^{-1} $')


plt.ylim(0.15,1)

plt.legend()
plt.grid()
plt.xlabel('ToT (ns)')
plt.ylabel('Pulsfläche A (nVs) ')


#import tikzplotlib
#tikzplotlib.save("area_tot.tex")



plt.show()
<file_sep>""" 
Little helper function to load data from a .trc binary file.
This is the file format used by LeCroy oscilloscopes.
<NAME> 09/2015
"""
import datetime
import numpy as np
import struct

def readTrc(fName, raw=None):
    """
        Reads .trc binary files from LeCroy Oscilloscopes.
        Decoding is based on LECROY_2_3 template.
        [More info](http://forums.ni.com/attachments/ni/60/4652/2/LeCroyWaveformTemplate_2_3.pdf)
        
        Parameters
        -----------       
        fName = filename of the .trc file
        
        Returns
        -----------       
        x: array with sample times [s],
        
        y: array with sample  values [V],
        
        d: dictionary with metadata
        
        
        <NAME> 09/2015
    """
    
    with open(fName, "rb") as fid:
        data = fid.read(50).decode()
        wdOffset = data.find('WAVEDESC')
        
        #-------------------------------
        # Get binary format / endianess
        #-------------------------------
        if readX( fid, '?', wdOffset + 32 ):  #16 or 8 bit sample format?
            smplFmt = "int16"
        else:
            smplFmt = "int8"
        if readX( fid, '?', wdOffset + 34 ):  #Big or little endian?
            endi = "<"
        else:
            endi = ">"
            
        #---------------------------------
        # Get length of blocks and arrays
        #---------------------------------
        lWAVE_DESCRIPTOR = readX( fid, endi+"l", wdOffset + 36 )
        lUSER_TEXT       = readX( fid, endi+"l", wdOffset + 40 )
        lTRIGTIME_ARRAY  = readX( fid, endi+"l", wdOffset + 48 )
        lRIS_TIME_ARRAY  = readX( fid, endi+"l", wdOffset + 52 )
        lWAVE_ARRAY_1    = readX( fid, endi+"l", wdOffset + 60 )
        lWAVE_ARRAY_2    = readX( fid, endi+"l", wdOffset + 64 )


        #--------------------------
        # Save formats and lengths
        #--------------------------
        info = dict()

        info["wdOffset"] = wdOffset
        info["smplFmt"] = smplFmt
        info["endi"] = endi
        info["lWAVE_DESCRIPTOR"] = lWAVE_DESCRIPTOR
        info["lUSER_TEXT"] = lUSER_TEXT
        info["lTRIGTIME_ARRAY"] = lTRIGTIME_ARRAY
        info["lRIS_TIME_ARRAY"] = lRIS_TIME_ARRAY
        info["lWAVE_ARRAY_1"] = lWAVE_ARRAY_1
        info["lWAVE_ARRAY_2"] = lWAVE_ARRAY_2

        
        #--------------------------------------
        # Get Instrument info from file header
        #--------------------------------------
        h = dict()  #Will store all the extracted Metadata

        h["INSTRUMENT_NAME"]  = readX( fid, "16s",    wdOffset + 76 ).decode().split('\x00')[0]
        h["INSTRUMENT_NUMBER"]= readX( fid, endi+"l", wdOffset + 92 )
        h["TRACE_LABEL"]      = readX( fid, "16s",    wdOffset + 96 ).decode().split('\x00')[0]
        
        #------------------------------------
        # Get Waveform info from file header
        #------------------------------------
        h["WAVE_ARRAY_COUNT"] = readX( fid, endi+"l", wdOffset + 116 )
        h["PNTS_PER_SCREEN"]  = readX( fid, endi+"l", wdOffset + 120 )
        h["FIRST_VALID_PNT"]  = readX( fid, endi+"l", wdOffset + 124 )
        h["LAST_VALID_PNT"]   = readX( fid, endi+"l", wdOffset + 128 )
        h["FIRST_POINT"]      = readX( fid, endi+"l", wdOffset + 132 )
        h["SPARSING_FACTOR"]  = readX( fid, endi+"l", wdOffset + 136 )
        h["SEGMENT_INDEX"]    = readX( fid, endi+"l", wdOffset + 140 )
        h["SUBARRAY_COUNT"]   = readX( fid, endi+"l", wdOffset + 144 )
        h["SWEEPS_PER_ACQ"]   = readX( fid, endi+"l", wdOffset + 148 )
        h["POINTS_PER_PAIR"]  = readX( fid, endi+"h", wdOffset + 152 )
        h["PAIR_OFFSET"]      = readX( fid, endi+"h", wdOffset + 154 )
        h["VERTICAL_GAIN"]    = readX( fid, endi+"f", wdOffset + 156 ) #to get floating values from raw data :
        h["VERTICAL_OFFSET"]  = readX( fid, endi+"f", wdOffset + 160 ) #VERTICAL_GAIN * data - VERTICAL_OFFSET 
        h["MAX_VALUE"]        = readX( fid, endi+"f", wdOffset + 164 )
        h["MIN_VALUE"]        = readX( fid, endi+"f", wdOffset + 168 )
        h["NOMINAL_BITS"]     = readX( fid, endi+"h", wdOffset + 172 )
        h["NOM_SUBARRAY_COUNT"]= readX( fid, endi+"h",wdOffset + 174 )
        h["HORIZ_INTERVAL"]   = readX( fid, endi+"f", wdOffset + 176 ) #sampling interval for time domain waveforms 
        h["HORIZ_OFFSET"]     = readX( fid, endi+"d", wdOffset + 180 ) #trigger offset for the first sweep of the trigger, seconds between the trigger and the first data point 
        h["PIXEL_OFFSET"]     = readX( fid, endi+"d", wdOffset + 188 )
        h["VERTUNIT"]         = readX( fid, "48s", wdOffset + 196 ).decode().split('\x00')[0]
        h["HORUNIT"]          = readX( fid, "48s", wdOffset + 244 ).decode().split('\x00')[0]
        h["HORIZ_UNCERTAINTY"]= readX( fid, endi+"f", wdOffset + 292 )
        h["TRIGGER_TIME"]     = getTimeStamp( fid, endi, wdOffset + 296 )
        Year, Month, Day, hour, minute, second = getTimeStampRaw( fid, endi, wdOffset + 296 )
        h["TS_Y"]             = Year
        h["TS_M"]             = Month
        h["TS_D"]             = Day
        h["TS_h"]             = hour
        h["TS_m"]             = minute
        h["TS_s"]             = second
        h["ACQ_DURATION"]     = readX( fid, endi+"f", wdOffset + 312 )
        h["RECORD_TYPE"]      = ["single_sweep","interleaved","histogram","graph","filter_coefficient","complex","extrema","sequence_obsolete","centered_RIS","peak_detect"][ readX( fid, endi+"H", wdOffset + 316 ) ]
        h["PROCESSING_DONE"]  = ["no_processing","fir_filter","interpolated","sparsed","autoscaled","no_result","rolling","cumulative"][ readX( fid, endi+"H", wdOffset + 318 ) ]
        h["RIS_SWEEPS"]       = readX( fid, endi+"h", wdOffset + 322 )
        h["TIMEBASE"]         = ['1_ps/div', '2_ps/div', '5_ps/div', '10_ps/div', '20_ps/div', '50_ps/div', '100_ps/div', '200_ps/div', '500_ps/div', '1_ns/div', '2_ns/div', '5_ns/div', '10_ns/div', '20_ns/div', '50_ns/div', '100_ns/div', '200_ns/div', '500_ns/div', '1_us/div', '2_us/div', '5_us/div', '10_us/div', '20_us/div', '50_us/div', '100_us/div', '200_us/div', '500_us/div', '1_ms/div', '2_ms/div', '5_ms/div', '10_ms/div', '20_ms/div', '50_ms/div', '100_ms/div', '200_ms/div', '500_ms/div', '1_s/div', '2_s/div', '5_s/div', '10_s/div', '20_s/div', '50_s/div', '100_s/div', '200_s/div', '500_s/div', '1_ks/div', '2_ks/div', '5_ks/div', 'EXTERNAL'][ readX( fid, endi+"H", wdOffset + 324 ) ]
        h["VERT_COUPLING"]    = ['DC_50_Ohms', 'ground', 'DC_1MOhm', 'ground', 'AC,_1MOhm'][ readX( fid, endi+"H", wdOffset + 326 ) ]
        h["PROBE_ATT"]        = readX( fid, endi+"f", wdOffset + 328 )
        h["FIXED_VERT_GAIN"]  = ['1_uV/div','2_uV/div','5_uV/div','10_uV/div','20_uV/div','50_uV/div','100_uV/div','200_uV/div','500_uV/div','1_mV/div','2_mV/div','5_mV/div','10_mV/div','20_mV/div','50_mV/div','100_mV/div','200_mV/div','500_mV/div','1_V/div','2_V/div','5_V/div','10_V/div','20_V/div','50_V/div','100_V/div','200_V/div','500_V/div','1_kV/div'][ readX( fid, endi+"H", wdOffset + 332 ) ]
        h["BANDWIDTH_LIMIT"]  = ['off', 'on'][ readX( fid, endi+"H", wdOffset + 334 ) ]
        h["VERTICAL_VERNIER"] = readX( fid, endi+"f", wdOffset + 336 )
        h["ACQ_VERT_OFFSET"]  = readX( fid, endi+"f", wdOffset + 340 )
        h["WAVE_SOURCE"]      = readX( fid, endi+"H", wdOffset + 344 )
        h["USER_TEXT"]        = readX( fid, "{0}s".format(lUSER_TEXT), wdOffset + lWAVE_DESCRIPTOR ).decode().split('\x00')[0]

        #----------------------------------------------------------
        # Get main sample data with the help of numpys .fromfile()
        #----------------------------------------------------------
        fid.seek( wdOffset + lWAVE_DESCRIPTOR + lUSER_TEXT + lTRIGTIME_ARRAY + lRIS_TIME_ARRAY ) #Seek to WAVE_ARRAY_1
        y = np.fromfile( fid, smplFmt, lWAVE_ARRAY_1 )
        if endi == ">":
            y.byteswap( True )
        if raw:
          x = np.arange(1,len(y)+1)
        else:
          y = h["VERTICAL_GAIN"] * y - h["VERTICAL_OFFSET"]
          x = np.arange(1,len(y)+1)*h["HORIZ_INTERVAL"] + h["HORIZ_OFFSET"]
    return x, y, h, info


def readX( fid, fmt, adr=None ):
    """ extract a byte / word / float / double from the binary file """

    #  struct format definitions
    #  Format	C Type	            Python type	        Standard size
    #  x	pad byte            no value	 	 
    #  c	char	            string of length 1  1
    #  b	signed char	    integer	        1
    #  B	unsigned char	    integer    	        1
    #  ?	_Bool	            bool                1
    #  h	short	            integer	        2
    #  H	unsigned short	    integer	        2
    #  i	int	            integer	        4
    #  I	unsigned int	    integer	        4
    #  l	long	            integer	        4
    #  L	unsigned long	    integer	        4
    #  q	long long	    integer	        8
    #  Q	unsigned long long  integer	        8
    #  f	float	            float	        4
    #  d	double	            float	        8
    #  s	char[]	            string	 	 
    #  p	char[]	            string	 	 
    #  P	void *	            integer

    nBytes = struct.calcsize( fmt )
    if adr is not None:
        fid.seek( adr )
    s = struct.unpack( fmt, fid.read( nBytes ) )
    if(type(s) == tuple):
        return s[0]
    else:
        return s


def getTimeStamp( fid, endi, adr, debug = None):
    """ extract a timestamp from the binary file """
    s = readX( fid, endi+"d", adr )
    m = readX( fid, endi+"b" )
    h = readX( fid, endi+"b" )
    D = readX( fid, endi+"b" )
    M = readX( fid, endi+"b" )
    Y = readX( fid, endi+"h" )
    trigTs = datetime.datetime(Y, M, D, h, m, int(s), int((s-int(s))*1e6) )
    if debug:
      print("  INFO::time", Y, M, D, h, m, s, int(s), int((s-int(s))*1e6))
    return trigTs

def getTimeStampRaw( fid, endi, adr, debug = None):
    """ extract a timestamp from the binary file """
    s = readX( fid, endi+"d", adr )
    m = readX( fid, endi+"b" )
    h = readX( fid, endi+"b" )
    D = readX( fid, endi+"b" )
    M = readX( fid, endi+"b" )
    Y = readX( fid, endi+"h" )
    if debug:
      print("  INFO time raw", Y, M, D, h, m, s)
    return Y, M, D, h, m, s
<file_sep>import numpy as np 
import matplotlib.pyplot as plt 
from scipy import stats 

#zylinder ohne  shaper
#E_real = np.array([80,186,242,295,352,32,662,1173,1332,1063,569])
#E_fit = np.array([0.97844,2.4133,3.137,3.7597,4.3779,0.3457,7.2341,10.7222,11.7675,10.2165,6.485])

#E_real = np.array([569, 662,1173,1332,1063])
#E_fit = np.array([6.472, 7.19,10.703,11.733,10.19])
#FWHM = np.array([0.51, 0.44, 0.87, 0.64, 0.48]) 

# zylinder mit shaper
E_real = np.array([662,1173,1332,569,1063])
E_fit = np.array([0.753,1.163,1.271,0.7018,1.1080])
 
FWHM_E = np.array([0.052, 0.069,0.0645, 0.065, 0.055])


#mit shaper ToT
#ToT_real = np.array([662, 1063,1252])
#ToT_fit = np.array([91.00116606316232,95.21278835637344,107.05266288124014])

#FWHM_T = np.array([9.798621093190247,15.77351053281185,22.63164384355587])
#x_err = np.array([0,0,79.5])



#tot mit shaper final

ToT_real = np.array([662,727,1069,1252])
ToT_fit = np.array([77.29,80.00,88.00,91.91]) 

FWHM_T = np.array([4.26,5,4,11]) 

#print(ToT_real)
#print(ToT_fit)

#GSI ToT

ToT_GSI = np.array([18.53,18.91,20.73,21.19,22.76])
E_real_GSI = np.array([569,662,1063,1173,1332])
#E_fit_GSI = np.array([0.701803960984309,0.7535521180502402,1.1657715187599953,1.273437535904615])
#ToT_GSI = np.array([18.53,18.95,21.2,22.75])

FWHM_GSI = np.array([1.88,1.68,2.81,3.31,3.46])



#bi207 1069keV 20,2

#plt.plot(E_fit,E_real,ls='None',marker='.')
#plt.errorbar(E_fit,E_real, xerr= 0.5*FWHM_E , fmt='.',capsize = 3,color='k', label='zugeordnete Energie') 

#plt.plot(ToT_fit,ToT_real, ls='None', marker='.')
#plt.errorbar(ToT_fit,ToT_real, xerr=0.5*FWHM_T,fmt='.',capsize = 3, label='zugeordnete Energie', color ='k') 

#plt.plot(test_tot,test_E, ls='None', marker='.')

plt.errorbar(ToT_GSI,E_real_GSI,xerr=0.5*FWHM_GSI,fmt='.',capsize = 3,color='k', label='zugeordnete Energie') 


#Energieauflösung

#E_fit_input = E_fit
#E_real_input = E_real
#FWHM_input = FWHM

#delta_E_fit = np.zeros([len(E_fit_input)])
#delta_E_real = np.zeros(len(E_fit_input))
#for i in range(0, len(E_fit_input), 1):
 #   print(i)
  #  delta_E_fit[i] = FWHM_input[i]/E_fit_input[i] 
   # delta_E_real[i] = delta_E_fit[i]*E_real_input[i]

#plt.plot(E_real, delta_E_fit, ls='None', marker = '.')


#print('%:',delta_E_fit)
#print('absolut:', delta_E_real)
#print(delta_E_fit*100)

from scipy.optimize import curve_fit

def fit_exp(x,a,b):
    return a*np.exp(b*x)

def fit_log(x,a,b):
    return a*np.log(x)+b

popt, pcov = curve_fit(fit_exp, ToT_GSI , E_real_GSI, bounds=([15,0],[35,1]) ,method ='trf')

#popt,pcov = curve_fit(fit, E_fit, E_real , bounds=([0,0,0],[1000,10,100]), method='trf'  )

perr = np.sqrt(np.diag(pcov))

print(popt)
print(perr)



#slope, intercept, r_value, p_value, std_err = stats.linregress(E_fit, E_real)

slope, intercept, r_value, p_value, std_err = stats.linregress(ToT_GSI, E_real_GSI)


fit_data = np.array([slope, intercept])
#print(slope)
#print(intercept)
#print(r_value)
#print(p_value)
#print(std_err)


x_fit_E = np.linspace(17,24,10)

plt.plot(x_fit_E, intercept + slope*x_fit_E, 'r',label='linearer Fit')

x_fit = np.linspace(17,24,1000)

#plt.plot(x_fit, fit_exp(x_fit, *popt), label='exp. Fit',color='r')

plt.legend()
plt.grid()
#plt.text(470,1.2, 'Q(E) = %5.5fE + %5.5f' %tuple(fit_data))
plt.text(17,1200, '$E(T) = m \cdot T + n$')
#plt.text(17,1200, '$E(T)=a \cdot e^{(b \cdot x)}$')
plt.text(17,1125, '$m= 187.69{keV/ns} $')
plt.text(17,1050, '$n= -2873.54 keV $')



plt.title("")
plt.ylabel("Energie E (keV)")
plt.xlabel("ToT T (ns)")

#import tikzplotlib

#tikzplotlib.save("calib_ToT_GSI_linear.tex")



plt.show()
<file_sep>#! /usr/bin/env python3

from __future__ import division
from __future__ import print_function

import os, sys, inspect

cmd_subfolder = os.path.realpath(os.path.abspath(os.path.join(os.path.split(inspect.getfile(inspect.currentframe()))[0],"src")))
if cmd_subfolder not in sys.path:
  sys.path.insert(0, cmd_subfolder)

from TLeCroy import *   # lecroy.py 
import matplotlib.pyplot as plt
from scipy.stats import norm
import numpy as np
import scipy.integrate as integrate
import matplotlib.mlab as mlab
from scipy.signal import find_peaks, peak_widths, peak_prominences 

# --------------------------------
# constants to be set by the user
# --------------------------------
channelNr = 4

 
StartTrace = 0
StopTrace =20



#inName = "--testpulse--"
#inName = "--testpulse--autosave--"
inName = "--testpulse--autosave--segments--"

inDir = "/Users/Lennard/Bachelor_Thesis_LF/BA_Lennard_Franz/cs137-longzylinder/"
#inDir = "/Users/Lennard/Bachelor_Thesis_LF/BA_Lennard_Franz/ra226-longzylinder/"
#inDir = "/volumes/lennad/segment/segment-CeBr-cs137/"
#inDir = "/volumes/lennad/<NAME>/th-zylinder/"
#inDir = "/volumes/lennad/<NAME>/cs137-longzylinder-shaper/"
#inDir = "/Users/Lennard/desktop/osci/lennard/segment-flat-CeBr-CS137/"
#inDir = "/Users/Lennard/desktop/osci/lennard/cs137-longzylinder-shaper/"
#inDir = "/Users/Lennard/destop/osci/lennard/cebr-big-4ch-cs137/"
#inDir = "/Users/Lennard/desktop/osci/lennard/segment-CeBr-ra226/"
#inDir = "/Users/Lennard/desktop/osci/lennard/ra226-filter/"
#inDir = "/Users/Lennard/desktop/autosave-segments/"
#inDir = "/Users/Lennard/desktop/osci/lennard/gag-am241/"


#inDir = "/home/lfranz/work/osci/lennard/segment-longcable-CeBr-CS137/"
#inDir = "/home/lfranz/work/osci/lennard/segment-flat-CeBr-CS137/"
#inDir = "/home/lfranz/work/osci/lennard/segment-flat-plexi-background/"
#inDir = "/home/lfranz/work/osci/lennard/segment-CeBr-am241/"
#inDir = "/home/lfranz/work/osci/lennard/segment-CeBr-cs137/"
#inDir = "/home/lfranz/work/osci/lennard/segment-CeBr-ra226/"

#inDir = "/home/lfranz/work/osci/lennard/segment/"
#inDir = "/home/lfranz/work/autosave-segments/"
#inDir = "/ZIH.fast/projects/d-lab_data/MAPMT/2019_gsi/segment/"
#inDir = "/ZIH.fast/projects/d-lab_data/MAPMT/2019_gsi/segments-new/"
#inDir = "/ZIH.fast/projects/d-lab_data/MAPMT/2019_gsi/segments-new/Autosave-new/"
# inDir = "/ZIH.fast/projects/d-lab_data/MAPMT/2019_gsi/Autosave/"
# --------------------------------

import numpy as np
import scipy.integrate as integrate


#empty array for area of sequnces for every trace
A = np.empty([0])
print(inDir)
for traceNr in range(StartTrace, StopTrace + 1, 1):
  print(" INFO: working on trace: ", "C"+ str(channelNr).zfill(1) + inName + str(traceNr).zfill(5) + ".trc")
 
  fileName = inDir + "C" + str(channelNr).zfill(1) + inName + str(traceNr).zfill(5) + ".trc"
  TLC = TLeCroy(fileName, debug=True)
  #TLC.PrintPrivate()
  # TLC.PlotTrace(raw=False)

  header = TLC.GetHeader()
  x, y = TLC.GetTrace()
  y = y - header['VERTICAL_OFFSET']
  y = -y 
  nano = 1e-9

  #print('----------------')
  #print('   vertical gain:', header['VERTICAL_GAIN'])
  #print('----------------')
  seqLength = header['WAVE_ARRAY_COUNT'] // header['SUBARRAY_COUNT']
  seqFirst = header['FIRST_VALID_PNT']
  #print('  LENGTH', seqLength)

  #plt.plot(x / nano, y)
  
  #arrays  
  idxHigh = np.zeros(header['SUBARRAY_COUNT'],dtype=np.int)
  idxLow = np.zeros(header['SUBARRAY_COUNT'],dtype=np.int)
  mean = np.zeros(header['SUBARRAY_COUNT'])
  std = np.zeros(header['SUBARRAY_COUNT'])
  area =  np.zeros(header['SUBARRAY_COUNT'])
  area_all = np.zeros(header['SUBARRAY_COUNT'])
  area_all_sum = np.zeros(header['SUBARRAY_COUNT'])
  area_bool =np.zeros(header['SUBARRAY_COUNT'])
  tot = np.zeros(header['SUBARRAY_COUNT'])

# determine background
  for seqNr in range(1, header['SUBARRAY_COUNT'] + 1):
    #bondarys of sequences
    idxHigh[seqNr-1] = seqFirst +( seqLength * seqNr - 1)
    idxLow[seqNr-1] = idxHigh[seqNr-1] - seqLength + 1
    #finding basline
    #y[...] has to be adjusted depending vertical position of the puls                     
    mean[seqNr-1], std[seqNr-1] = norm.fit(y[ idxLow[seqNr-1] : idxLow[seqNr-1] + int(0.35*seqLength)])

  mean_all = sum(mean)/header['SUBARRAY_COUNT']
  std_all = sum(std)/header['SUBARRAY_COUNT'] 
  #Histogram of values for baseline 
  #plt.hist(y[ idxLow[seqNr-1] : idxLow[seqNr-1] + int(0.35*seqLength)])
  

#determination of area of pulses
  for seqNr in range(1, header['SUBARRAY_COUNT'] + 1):
   # if mean[seqNr -1] > 0.0285 and mean[seqNr - 1] < 0.0315 and np.min( y[idxLow[seqNr-1] : idxHigh[seqNr-1]]-mean[seqNr-1]) > -0.0035:

#area of the pulse
    #area[seqNr-1] = np.trapz((y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]- mean[seqNr-1]) > 3*std[seqNr-1] ,
                 #x[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]/nano)  
 
    #area_all[seqNr-1] = np.sum(y[idxLow[seqNr - 1] : idxHigh[seqNr - 1]])/(nano/header['HORIZ_INTERVAL']) - ((mean_all)*seqLength*header['HORIZ_INTERVAL']/nano)
   
    #print(area_all_sum)
    #print(area_all)
    #condition to differentiate values of pulse from background 
    #boolArr = np.where((y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]- mean[seqNr - 1]) > (3*std[seqNr -1]))
    boolArr = (y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]- mean[seqNr-1]) > (3*std[seqNr-1])

    #sum of all values which meet condition and scaling 
    #substracting of area belonging to background
    
    #print(y[boolArr[0]] - mean[seqNr-1])
    #area_bool[seqNr-1] = np.sum(y[boolArr[0]])*(header['HORIZ_INTERVAL']/nano) - (mean[seqNr-1])*(len(boolArr[0])*header['HORIZ_INTERVAL']/nano)
    



    area_bool[seqNr-1] = np.sum((y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]*(header['HORIZ_INTERVAL']/nano)),where = boolArr) - mean[seqNr-1]*np.sum(boolArr)*header['HORIZ_INTERVAL']/nano
    tot[seqNr-1] = (np.sum(boolArr)-1)*header['HORIZ_INTERVAL']/nano

    #print(area_bool[seqNr - 1],tot[seqNr - 1] )

    


    #area_bool[seqNr-1] = np.sum(y[boolArr[0]]) 
     #- (mean[seqNr-1])*(len(boolArr[0]))


    #print(len(boolArr[0])*header['HORIZ_INTERVAL']/nano)    
    #plt.hist(len(boolArr[0])*header['HORIZ_INTERVAL']/nano)

    #plt.plot(x[idxLow[seqNr - 1]:idxHigh[seqNr - 1]] / nano, y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]- mean_all)
 # plt.plot(x[idxLow[2]:idxHigh[2]] / nano, y[idxLow[2]:idxHigh[2]]- mean_all)
  


  #print("-----------------------")
  #print("mean all:",mean_all,"sigma all:", std_all)
  #print("array of mean of segments:",mean,"array of sigma of segments",std) 
  #print("-----------------------")
  #appending area of segments to A for every trace
  #length of A is number of sequnces * number of traces 
  #print(y[boolArr[0]] - mean[seqNr-1])
  A = np.append(A,area_bool)

#<NAME>
E_real = np.array([569, 662,1173,1332,1063])
E_fit = np.array([6.485, 7.2341,10.7222,11.7675,10.2165])

#E_real = np.array([80,186,242,295,352,32,662,1173,1332,1063,569])
#E_fit = np.array([0.97844,2.4133,3.137,3.7597,4.3779,0.3457,7.2341,10.7222,11.7675,10.2165,6.485])


#mit shaper
#E_real = np.array([662,1173,1332,569,1063])
#E_fit = np.array([0.7535521180502402,1.1657715187599953,1.273437535904615,0.701803960984309,1.1086378999695787])



from scipy import stats 
slope, intercept, r_value, p_value, std_err = stats.linregress(E_fit, E_real)

#shaper wolfram alpha
#slope = 1298.14
#intercept = -339.482

#shaper python
#slope = 1298.93696929096
#intercept= -340.41065484554497

#no shaper wolfram
#slope = 143.049
#intercept = -363.694


print(intercept)
print(slope)

E = intercept+slope*A
print(E)




#histogram of all sequences of all traces       
n,bins,patches = plt.hist(E,bins=500,histtype='step', color='k')
#print(n,bins,patches)
plt.grid(True)

bin_width= (bins[1]-bins[0])

#plt.xlim(0.,15.0)

peaks,properties = find_peaks(n,width=3, distance=1,height=300, rel_height=0.7)
print("bins of peaks:", peaks)
print("X Values of peaks:", bins[peaks] )
print(properties)

#print(bins[peaks])
#print(bins[np.int(peaks[6] - 1.5* properties["widths"][6]) :np.int( peaks[6]+ 1.5* (properties["widths"][6])) ]) 
#print(peaks[2] - 1.5* properties["widths"][2])
#print(np.int(peaks[2] - 1.5* properties["widths"][2]))

#print(n,bins)

#peakNr = 6
#print (bins[np.int(peaks[peakNr] - 1.5 * properties["widths"][peakNr]) :np.int( peaks[peakNr]+ 1.5 * (properties["widths"][peakNr])) ])
#print (n[np.int(peaks[peakNr] - 1.5* properties["widths"][peakNr]) :np.int( peaks[peakNr]+ 1.5* (properties["widths"][peakNr]))  ] )












from scipy.optimize import curve_fit
def gauss (x,a,sigma,mu):
  #return (1/np.sqrt(2*np.pi*sigma**2))* np.exp(- (x-mu)**2 / 2*sigma**2)
  return a*np.exp(-(x-mu)**2/(2*sigma**2)) 


bins = np.delete(bins,[len(bins)-1])


for peakNr in range(1, len(peaks) ):
    
  fit_width = 0.7
  x_values = bins[np.int(peaks[peakNr] - fit_width * properties["widths"][peakNr]) :np.int( peaks[peakNr]+ fit_width * (properties["widths"][peakNr])) ]
  y_values = n[np.int(peaks[peakNr] - fit_width * properties["widths"][peakNr]) :np.int( peaks[peakNr]+ fit_width * (properties["widths"][peakNr]))  ]


  #print(x_values)
  #print(y_values)

  a_bound = [0, 1.1*np.max(y_values)]
  sigma_bound = [0, 50]
  mu_bound = [np.min(x_values), np.max(x_values)]

  #print(a_bound)
  #print(sigma_bound)
  #print(mu_bound)

  #popt, pcov = curve_fit(gauss, bins[np.int(peaks[peakNr] - 1.5 * properties["widths"][peakNr]) :np.int( peaks[peakNr]+ 1.5 * (properties["widths"][peakNr])) ] , n[np.int(peaks[peakNr] - 1.5* properties["widths"][peakNr]) :np.int( peaks[peakNr]+ 1.5* (properties["widths"][peakNr]))  ] )

  popt, pcov = curve_fit(gauss, x_values , y_values, bounds= ([a_bound[0],sigma_bound[0],mu_bound[0]], [a_bound[1],sigma_bound[1],mu_bound[1]] ) ) 

  plt.plot(bins, gauss(bins, *popt),color='g')
  plt.plot((popt[2], popt[2]), (0, 1600), 'r-')
  
  #statistic error of curve fit 
  perr = np.sqrt(np.diag(pcov))
  
  #FWHM
  FWHM = 2* np.sqrt(2*np.log(2)) * popt[1]
  delta_FWHM = 2* np.sqrt(2*np.log(2)) * perr[1]

  
  print("Peak Nr:" , peakNr) 
  print("Energy:" , popt[2], "statistic error of energy:" , perr[2], "FWHM:", FWHM, "statistic error of FWHM", delta_FWHM )
  




plt.title("Histogramm der deponierten Ladung im Detektor")
plt.xlabel("deponierte Ladung (nVs)")
plt.ylabel("Anzahl N ")
plt.xlim(-500,1000)
plt.ylim(0,1300)

#determination of mean and sigma of data in histogram
#MEAN, STD = norm.fit(A)
#print("--------------------------")
#print("Entries:" , len(A))
#print("--------------------------")
#print("Mean of pulse energy :", MEAN ,"Sigma of pulse energy:", STD) 
#print("--------------------------")

plt.show()
<file_sep>#! /usr/bin/env python3

from __future__ import division
from __future__ import print_function

import os, sys, inspect

cmd_subfolder = os.path.realpath(os.path.abspath(os.path.join(os.path.split(inspect.getfile(inspect.currentframe()))[0],"src")))
if cmd_subfolder not in sys.path:
  sys.path.insert(0, cmd_subfolder)

from TLeCroy import *   # lecroy.py 
import matplotlib.pyplot as plt
from scipy.stats import norm
import numpy as np
import scipy.integrate as integrate
import matplotlib.mlab as mlab
from scipy.signal import find_peaks

# --------------------------------
# constants to be set by the user
# --------------------------------
channelNr = 4

 
StartTrace = 0
StopTrace = 25



#inName = "--testpulse--"
#inName = "--testpulse--autosave--"
inName = "--testpulse--autosave--segments--"

inDir = "/Users/Lennard/Bachelor_Thesis_LF/BA_Lennard_Franz/cs137-longzylinder/"

#inDir = "/volumes/lennad/segment/segment-CeBr-cs137/"
#inDir = "/volumes/lennad/<NAME>anz/cs137-longzylinder/"
#inDir = "/volumes/lennad/<NAME>anz/th232-cebr-4x4/"
#inDir = "/Users/Lennard/desktop/osci/lennard/segment-flat-CeBr-CS137/"
#inDir = "/Users/Lennard/desktop/osci/lennard/cs137-longzylinder-shaper/"
#inDir = "/Users/Lennard/destop/osci/lennard/cebr-big-4ch-cs137/"
#inDir = "/Users/Lennard/desktop/osci/lennard/segment-CeBr-am241/"
#inDir = "/Users/Lennard/desktop/osci/lennard/ra226-filter/"
#inDir = "/Users/Lennard/desktop/autosave-segments/"
#inDir = "/Users/Lennard/desktop/osci/lennard/gag-am241/"


#inDir = "/home/lfranz/work/osci/lennard/segment-longcable-CeBr-CS137/"
#inDir = "/home/lfranz/work/osci/lennard/segment-flat-CeBr-CS137/"
#inDir = "/home/lfranz/work/osci/lennard/segment-flat-plexi-background/"
#inDir = "/home/lfranz/work/osci/lennard/segment-CeBr-am241/"
#inDir = "/home/lfranz/work/osci/lennard/segment-CeBr-cs137/"
#inDir = "/home/lfranz/work/osci/lennard/segment-CeBr-ra226/"

#inDir = "/home/lfranz/work/osci/lennard/segment/"
#inDir = "/home/lfranz/work/autosave-segments/"
#inDir = "/ZIH.fast/projects/d-lab_data/MAPMT/2019_gsi/segment/"
#inDir = "/ZIH.fast/projects/d-lab_data/MAPMT/2019_gsi/segments-new/"
#inDir = "/ZIH.fast/projects/d-lab_data/MAPMT/2019_gsi/segments-new/Autosave-new/"
# inDir = "/ZIH.fast/projects/d-lab_data/MAPMT/2019_gsi/Autosave/"
# --------------------------------

import numpy as np
import scipy.integrate as integrate

#empty array for area of sequnces for every trace
A = np.empty([0])
print(inDir)
for traceNr in range(StartTrace, StopTrace + 1, 1):
  print(" INFO: working on trace: ", "C"+ str(channelNr).zfill(1) + inName + str(traceNr).zfill(5) + ".trc")
 
  fileName = inDir + "C" + str(channelNr).zfill(1) + inName + str(traceNr).zfill(5) + ".trc"
  TLC = TLeCroy(fileName, debug=True)
  #TLC.PrintPrivate()
  # TLC.PlotTrace(raw=False)

  header = TLC.GetHeader()
  x, y = TLC.GetTrace()
  y = y - header['VERTICAL_OFFSET']
  y = -y 
  nano = 1e-9

  #print('----------------')
  #print('   vertical gain:', header['VERTICAL_GAIN'])
  #print('----------------')
  seqLength = header['WAVE_ARRAY_COUNT'] // header['SUBARRAY_COUNT']
  seqFirst = header['FIRST_VALID_PNT']
  #print('  LENGTH', seqLength)

  #plt.plot(x / nano, y)
  
  #arrays  
  idxHigh = np.zeros(header['SUBARRAY_COUNT'],dtype=np.int)
  idxLow = np.zeros(header['SUBARRAY_COUNT'],dtype=np.int)
  mean = np.zeros(header['SUBARRAY_COUNT'])
  std = np.zeros(header['SUBARRAY_COUNT'])
  area =  np.zeros(header['SUBARRAY_COUNT'])
  area_all = np.zeros(header['SUBARRAY_COUNT'])
  area_all_sum = np.zeros(header['SUBARRAY_COUNT'])
  area_bool =np.zeros(header['SUBARRAY_COUNT'])
  
# determine background
  for seqNr in range(1, header['SUBARRAY_COUNT'] + 1):
    #bondarys of sequences
    idxHigh[seqNr-1] = seqFirst +( seqLength * seqNr - 1)
    idxLow[seqNr-1] = idxHigh[seqNr-1] - seqLength + 1
    #finding basline
    #y[...] has to be adjusted depending vertical position of the puls                     
    mean[seqNr-1], std[seqNr-1] = norm.fit(y[ idxLow[seqNr-1] : idxLow[seqNr-1] + int(0.35*seqLength)])

  mean_all = sum(mean)/header['SUBARRAY_COUNT']
  std_all = sum(std)/header['SUBARRAY_COUNT'] 
  #Histogram of values for baseline 
  #plt.hist(y[ idxLow[seqNr-1] : idxLow[seqNr-1] + int(0.35*seqLength)])
  

#determination of area of pulses
  for seqNr in range(1, header['SUBARRAY_COUNT'] + 1):
   # if mean[seqNr -1] > 0.0285 and mean[seqNr - 1] < 0.0315 and np.min( y[idxLow[seqNr-1] : idxHigh[seqNr-1]]-mean[seqNr-1]) > -0.0035:

#area of the pulse
    #area[seqNr-1] = np.trapz((y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]- mean[seqNr-1]) > 3*std[seqNr-1] ,
                 #x[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]/nano)  
 
    #area_all[seqNr-1] = np.sum(y[idxLow[seqNr - 1] : idxHigh[seqNr - 1]])/(nano/header['HORIZ_INTERVAL']) - ((mean_all)*seqLength*header['HORIZ_INTERVAL']/nano)
   
    #print(area_all_sum)
    #print(area_all)
    #condition to differentiate values of pulse from background 
    #boolArr = np.where((y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]- mean[seqNr - 1]) > (3*std[seqNr -1]))
    boolArr = (y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]- mean[seqNr-1]) > (3*std[seqNr-1])

    #sum of all values which meet condition and scaling 
    #substracting of area belonging to background
    
    #print(y[boolArr[0]] - mean[seqNr-1])
    #area_bool[seqNr-1] = np.sum(y[boolArr[0]])*(header['HORIZ_INTERVAL']/nano) - (mean[seqNr-1])*(len(boolArr[0])*header['HORIZ_INTERVAL']/nano)
    area_bool[seqNr-1] = np.sum((y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]*(header['HORIZ_INTERVAL']/nano)),where = boolArr) - mean[seqNr-1]*np.sum(boolArr)*header['HORIZ_INTERVAL']/nano

    #area_bool[seqNr-1] = np.sum(y[boolArr[0]]) 
     #- (mean[seqNr-1])*(len(boolArr[0]))


    #print(len(boolArr[0])*header['HORIZ_INTERVAL']/nano)    
    #plt.hist(len(boolArr[0])*header['HORIZ_INTERVAL']/nano)

    #plt.plot(x[idxLow[seqNr - 1]:idxHigh[seqNr - 1]] / nano, y[idxLow[seqNr - 1]:idxHigh[seqNr - 1]]- mean_all)
 # plt.plot(x[idxLow[2]:idxHigh[2]] / nano, y[idxLow[2]:idxHigh[2]]- mean_all)
  


  #print("-----------------------")
  #print("mean all:",mean_all,"sigma all:", std_all)
  #print("array of mean of segments:",mean,"array of sigma of segments",std) 
  #print("-----------------------")
  #appending area of segments to A for every trace
  #length of A is number of sequnces * number of traces 
  #print(y[boolArr[0]] - mean[seqNr-1])
  A = np.append(A,area_bool)




#peaks = find_peaks_cwt(A)
#print(peaks)
#histogram of all sequences of all traces       
n,bins,patches=plt.hist(A,bins=1000,histtype='step', color = 'k')
#print(n,bins)
plt.grid(True)

plt.xlim(0.,10.0)

peaks = find_peaks(n,height=200, width=2, distance=5)
print( bins[peaks[0]])


plt.title("histogram of pulse energy")
plt.xlabel("energy [nVs]")
plt.ylabel("counts N")
#determination of mean and sigma of data in histogram
MEAN, STD = norm.fit(A)
#print("--------------------------")
#print("Entries:" , len(A))
#print("--------------------------")
#print("Mean of pulse energy :", MEAN ,"Sigma of pulse energy:", STD)
#print("--------------------------")

plt.show()
