Plotting Pulses Without Cuts¶
Below are all of our import statements. You probably won't need to mess with these, but if you want to use a new python package put it here!
%matplotlib inline
import h5py
import numpy as np
import matplotlib.pyplot as plt
import helperfunctions as hf
Convert a raw data file to a set of pulse blocks¶
Be sure to change to the run number you want to study
%run ../Scripts/convert_hdf5.py -r Run_1235
%run ../Physics/pulseAnalysis.py -r Run_1235
Get the directories in which we store our processed pulses and pulse blocks¶
dataDir contains our processed, uninterlaced pulses
blocksDir contains our interlaced pulse blocks
run = 'Run_1235'
dataDir = hf.checkDir(hf.getRootDir()+"/Data/Processed/Pulses/"+run+"/")
blocksDir = hf.checkDir(hf.getRootDir()+"/Physics/Blocks_Final/"+run+"/")
Here we print out what amplitudes are available in this dataset, as well as which amplitude we are currently looking at
amps = hf.getAmps(dataDir)
print('The available amplitudes are:', ', '.join(amps))
amp = amps[0] # if your data has more than one amplitude, changing 0 will allow you to look at those too
print('The amplitude you are using is: ' + amp)
Loading in a dataset¶
Here we finally load in a dataset. inFile is our .hdf5 file containing our pulse block, and channel1_data is a numpy array of pulses from channel 1 of the first COLUTA chip
inFile = h5py.File(blocksDir+"Blocks_Amp"+amp+"_1x.hdf5", "r")
channel1_data = inFile["/coluta1/channel1/samples"][()][:,0:1200]
Let's see what the shape of this dataset is
channel1_data.shape
This means our dataset stored in channel1_data has 1000 rows, and each row has 1200 entries. So how does this relate to the data we got from the testboard? We took 1000 measurements with the testboard, and we used data from each measurement to create a pulse that was 1200 samples long.
Making Plots¶
Now let's look at what a single pulse looks like
plt.plot(channel1_data[0])
plt.title('COLUTA 1, Channel 1, Measurement 0, Amplitude ' + amp)
plt.show()
This looks good. We can clearly see the pulse rising, falling to some constant negative value, then rising to a the baseline at 0. There is a fair amount of noise though, so let's get rid of that by averaging our 1000 pulses
channel1_data_averaged = np.mean(channel1_data,axis=0)
plt.plot(channel1_data_averaged)
plt.title('COLUTA 1, Channel 1, Averaged Pulse, Amplitude ' + amp)
plt.show()
That seems to have gotten rid of a fair amount of the noise. Since we can, let's take the derivative.
channel1_data_derivative = np.gradient(channel1_data_averaged)
plt.plot(channel1_data_derivative)
plt.title('COLUTA 1, Channel 1, Pulse Derivative, Amplitude ' + amp)
plt.show()
It's important to see how noisy our data is, so let's plot the standard deviation of each point in channel1_data
plt.plot(np.std(channel1_data,axis=0))
plt.title('COLUTA 1, Channel 1, Std Dev, Amplitude ' + amp)
plt.show()
You might notice that there's a region from around 200 to 600 with noticeably higher noise. Let's spend some time looking at that.
First we'll plot the above noise plot (shifted and scaled to be visible) and the pulse on the same axes.
plt.plot((np.std(channel1_data,axis=0)-1.2)*500, label = 'Std. Dev. of Pulse')
plt.plot(channel1_data_averaged, label = 'Averaged Pulse')
plt.title('COLUTA 1, Channel 1, Std Dev Comparison, Amplitude ' + amp)
plt.legend()
plt.show()
Hmm, seems like the noise is higher when the pulse is nonzero. Let's also plot the noise alongside the derivative.
plt.plot((np.std(channel1_data,axis=0)-1.2)*20, label = 'Std. Dev. of Pulse')
plt.plot(channel1_data_derivative, label = 'Pulse Derivative')
plt.title('COLUTA 1, Channel 1, Std Dev Comparison, Amplitude ' + amp)
plt.legend()
plt.show()
Okay, here we can see that it actually looks like the noise increases when the derivative is nonzero. It's always interesting what you can learn by making comparison plots!
Now we'll look at various points along the pulse. These points are listed in points_dict
points_dict = {
"20":"Baseline Before Pulse",
"142":"Rising Edge",
"165":"Peak",
"195":"Falling Edge",
"500":"Negative Lobe",
"1000":"Baseline After Pulse"
}
for point, description in points_dict.items():
data = channel1_data[:,int(point)]
ax=plt.subplot(111)
plt.hist(data, bins = range(int(min(data)), int(max(data))+1))
plt.title("Sample "+point+": "+description+", Channel "+str(channel+1))
plt.text(0.95, 0.85, r'$\mu = $'+str(np.round(np.mean(data), 3)), horizontalalignment='right',
verticalalignment='center', transform=ax.transAxes)
plt.text(0.95, 0.75, r'$\sigma = $'+str(np.round(np.std(data), 3)), horizontalalignment='right',
verticalalignment='center', transform=ax.transAxes)
plt.show()
There's quite a bit happening to make these plots, so make sure you know what the code in this above block is doing. As always when coding, Google is your friend.
Now try making a few plots of your own. I suggest:
- The derivative of a single, unaveraged pulse. This will look like nonsense, and that's why we average out the noise!
- The std. dev. of the pulse alongside the pulse derivative squared. We expect the noise to approximately follow $$ \sigma_s^2 = \sigma_n^2 + \sigma_t^2 \Big( \frac{\partial V}{\partial t} \Big)^2 $$ where the various $\sigma$'s are noise values characteristic of the testboard. Does your plot agree with this?
- You may have noticed in these noise histograms that there are some weird gaps. This is what we've been calling the "even/odd" problem. Make a short script to look at every sample in our data and determine if it's even or odd, then print out the final tally.