wavecatcher-analysis
Loading...
Searching...
No Matches
ReadRun.h
Go to the documentation of this file.
1
2#ifndef _ReadRun
3#define _ReadRun
4
5// contact: christian.scharf@physik.hu-berlin.de
6
7#include <Math/Minimizer.h>
8#include <Math/Factory.h>
9#include <Math/Functor.h>
10#include <TClonesArray.h>
11#include <TObjString.h>
12#include <TObject.h>
13#include <TLine.h>
14#include <TGraph.h>
15#include <TGraph2D.h>
16#include <TMultiGraph.h>
17#include <TString.h>
18#include <TCanvas.h>
19#include <TFile.h>
20#include <TTree.h>
21#include <TH1D.h>
22#include <TMath.h>
23#include <TF1.h>
24#include <TStyle.h>
25#include <TH1F.h>
26#include <TH2D.h>
27#include <TEfficiency.h>
28#include <TLegend.h>
29#include <THStack.h>
30#include <THistPainter.h>
31#include <TText.h>
32#include <TFitResultPtr.h>
33#include <TFitResult.h>
34#include <TSpline.h>
35#include <TPaveStats.h>
36#include <TError.h> // root verbosity level
37#include <TSystem.h>
38#include <TROOT.h>
39#include <TLatex.h>
40
41//C, C++
42#include <stdio.h>
43#include <future>
44#include <assert.h>
45#include <iostream>
46#include <iomanip>
47#include <fstream>
48#include <vector>
49#include <stdlib.h>
50#include <sstream>
51#include <stdexcept>
52//multithreading (avoids finding omp.h manually)
53#ifndef OMP_FUNC
54extern "C" {
57}
58#endif
59
60#include "utils/FitFunctions.h"
61#include "utils/Helpers.h"
62#include "utils/Filters.h"
63
64using namespace std;
65
66class ReadRun {
67private:
68#pragma pack(1) // byte padding suppression for WaveCatcher data format
69 // structs copied from
70 // WaveCatcher binary -> root converter
71 // by manu chauveau@cenbg.in2p3.fr
72 struct event_data
73 {
74 int EventNumber;
75 double EpochTime;
76 unsigned int Year;
77 unsigned int Month;
78 unsigned int Day;
79 unsigned int Hour;
80 unsigned int Minute;
81 unsigned int Second;
82 unsigned int Millisecond;
83 unsigned long long int TDCsamIndex;
84 int nchannelstored;
85 };
86
87 struct channel_data_without_measurement
88 {
89 int channel;
90 int EventIDsamIndex;
91 int FirstCellToPlotsamIndex;
92 };
93
94 struct channel_data_with_measurement
95 {
96 int channel;
97 int EventIDsamIndex;
98 int FirstCellToPlotsamIndex;
99 float MeasuredBaseline;
100 float AmplitudeValue;
101 float ComputedCharge;
102 float RiseTimeInstant;
103 float FallTimeInstant;
104 float RawTriggerRate;
105 };
106#pragma pack() // byte padding suppression for WaveCatcher data format
107
108protected:
117
119 virtual void checkData(bool isBaselineCorrection = false) const {
121 if (eventnr_storage.empty()) {
122 throw runtime_error(
123 "Error: No data has been loaded yet!\n \
124 Please call ReadFile() before calling any functions which manipulate data.\n \
125 Aborting execution."
126 );
127 }
128 }
129
130public:
133
136
139
140 // plots amplValuessum
141 void PlotChannelSums(bool = false, bool = false, double = 0., double = 0., int = 2);
142
143 void PlotChannelAverages(bool = false);
144
145 TH2F* WFHeatmapChannel(int, float = -20, float = 200, int = 880);
146 void PlotWFHeatmaps(float = -20, float = 200, int = 880, string = "", float = 0, EColorPalette = kRainBow);
147
148 // baseline correction (shifts all waveforms individually in y)
149 void CorrectBaseline(float, float = -999);
150 void CorrectBaseline_function(vector<float>&, float, float, int);
151
152 void CorrectBaselineMinSlopeRMS(vector<float>, double = 0, int = 2);
153 void CorrectBaselineMinSlopeRMS(vector<float>, double, int, int);
154 void CorrectBaselineMinSlopeRMS(int = 100, bool = false, double = 0.5, int = 0, int = 0, int = 2);
155
156 void CorrectBaselineMin(vector<float>, double = 0, int = 2);
157 void CorrectBaselineMin(vector<float>, double, int, int);
158 void CorrectBaselineMin(int = 100, double = 0.5, int = 0, int = 0, int = 2);
159
160 // functions to check baseline correction results
161 TH1F* WFProjectionChannel(int, int = 0, int = 1024, float = -50, float = 50, int = 200);
162 void PrintWFProjection(float = 0, float = 320, float = -50, float = 50, int = 200);
163
164 TH1F* BaselineCorrectionResults(int, int, float = -5, float = 5, int = 200);
165 void PrintBaselineCorrectionResults(float = -5, float = 5, int = 200);
166
167 // get timing of peaks
168 void GetTimingCFD(float = .3, float = 100, float = 140, double = 0., bool = true, int = 2, bool = false, bool = false);
169 void SkipEventsTimeDiffCut(int, int, double, double, bool = false);
170
171 // average all waveforms to simplify peak ID
172 void SmoothAll(double, int);
173 void SmoothAll(double, string = "Gaus");
174 void FilterAll(double = .3, double = .9, double = .2);
175 void ShiftAllToAverageCF();
176
177 // functions for charge spectrum
178 array<int, 3> GetIntWindow(TH1F*, float, float, float, float, int = 0);
179 array<int, 3> GetIntWindow(const vector<float>&, float, float, float, float, int = 0);
180 array<int, 3> GetIntWindow(const vector<float>&, int, int, int, int, int = 0);
181 float GetPeakIntegral(TH1F*, float, float, float, float, int = 0);
182 float GetPeakIntegral(const vector<float>&, float, float, float, float, int = 0);
183 void PrintChargeSpectrumWF(float, float, float = 0, float = 300, int = 1, float = 0., float = 0., float = 0., float = 0.);
184 TH1F* ChargeSpectrum(int, float, float, float = 0, float = 300, float = -50, float = 600, int = 750);
185 void PrintChargeSpectrum(float, float, float = 0, float = 300, float = -50, float = 600, int = 750, float = 0., float = 0., int = 99, int = 0, bool = false);
188
189
190 float* ChargeList(int, float = 20, float = 80, float = 0, float = 300, bool = 1);
191 void SaveChargeLists(float = 20, float = 80, float = 0, float = 300, bool = 1);
192 void ChargeCorrelation(float, float, float, float, float, float, int, int, int, bool = false);
193
194 // SiPM specific
195 void PrintDCR(float = 15, float = 85, float = 0, float = 300, double = 3); // based on PMT::PrintChargeSpectrumPMTthreshold()
196
197 // functions for time distribution
198 TH1F* TimeDist(int, float = 0, float = 300, float = 0, float = 300, int = 100, int = 0, float = .3);
199 void PrintTimeDist(float = 0, float = 300, float = 0, float = 300, int = 100, int = 0, float = .3);
200 TGraph2D* MaxDist(int, float = 0, float = 300);
201 void PrintMaxDist(float = 0, float = 300);
202
203 TH1F* His_GetTimingCFD(int, float, float, int = -999);
204 void Print_GetTimingCFD(float = 100, float = 140, int = 0, int = -999, string = "S", bool = true);
205 TH1F* His_GetTimingCFD_diff(vector<int>, vector<int>, float, float, int = -999);
206 void Print_GetTimingCFD_diff(vector<int>, vector<int>, float = 100, float = 140, int = 0, int = -999, float = -999, float = -999, string = "RS", bool = true);
207
208
209 // helper functions
210 TH1F* Getwf(int); // waveform number
211 virtual TH1F* Getwf(int, int, int = 1); // channel, eventnr, color
212 template<typename T>
213 T* getx(double shift = 0.); // x values
214 double* gety(int); // y values for waveform index
215 double* gety(int, int); // y values for waveform(channel, event)
216
217 static pair<float, bool> LinearInterpolation(float, float, float, float, float, bool = false); // linear interpolation
218
219 virtual int GetWaveformIndex(int, int); // get index of waveform from event and channel
220 int GetEventIndex(unsigned int); // get index of a triggered event (finds the correct event if files are not read sequentially)
221 int GetChannelIndex(int); // get index of a certain channel
222 virtual int GetCurrentChannel(int); // get index of channel for a certain waveform
223 virtual int GetCurrentEvent(int); // get index of event for a certain waveform
224
225 bool PlotChannel(int); // check if channel should be plotted
226
227 vector<bool> PolarityMap(bool, int); // check if channel should be inverted
228
229 int CheckBoundsX(int); // Check if index exists in time of waveforms
230 int TimeToIndex(float); // Convert time to the bin number of the waveform
231 float IndexToTime(int); // Convert the bin number of the waveform to the time
232
239 ReadRun(int last_bin_file = 0, int first_bin_file = 0);
240
241 virtual void ReadFile(string, bool = false, int = 9, string = "out.root", bool = false, long long = -1);
242
243 virtual ~ReadRun();
244
248 string data_path;
249
258
259
267
280
282 int nevents = 0;
286 int nwf;
290 int maxNWF = 1e7;
291
293 float SP = .3125;
295 float SP_inv = 1 / .3125;
299 float DAQ_factor = 250. / 4096.;
303 int nChannelsWC = 64;
304
309
317
322
325
328
333 virtual int Nevents_good(int = 0);
334
335 void SkipEventsPerChannel(vector<float>, float = 0, float = 0, bool = false); // in case you want to have indiviual thresholds in individual channels
336 void IntegralFilter(vector<float>, vector<bool>, float, float, float = 50, float = 250, bool = false, bool = false); // Same as SkipEventsPerChannel() but filtering all events with integrals <(>) threshold
337 void PrintSkippedEvents();
338 void UnskipAll();
339 virtual bool SkipEvent(int, int = -1);
340
343
358
361
362 //other controls
363
368 float tCutg;
370 float tCutEndg;
371
376};
377#endif
int omp_get_max_threads()
Main class containing the file reader and most analysis functions.
int omp_get_thread_num()
int FirstBinFileToRead
Number first of .bin file to be read in.
Definition ReadRun.h:257
static pair< float, bool > LinearInterpolation(float, float, float, float, float, bool=false)
Simple linear interpolation for x.
Definition ReadRun.cc:2817
TH1F * Getwf(int)
Helper that returns the waveform histogram for a certain waveform number number.
Definition ReadRun.cc:2683
int maxNWF
Maximum possible number of waveforms in data: number of active channels x number of events.
Definition ReadRun.h:290
void Print_GetTimingCFD(float=100, float=140, int=0, int=-999, string="S", bool=true)
Plot results of GetTimingCFD()
Definition ReadRun.cc:2429
int PlotChannelAverages_cnt
Index for multiple executions of the same plotting function.
Definition ReadRun.h:112
void ChargeCorrelation(float, float, float, float, float, float, int, int, int, bool=false)
Plot correlation of integrals/amplitudes between two channels.
Definition ReadRun.cc:1856
float tCutEndg
End of time window for baseline correction when using ReadRun::Using_BaselineCorrection_in_file_loop.
Definition ReadRun.h:370
TH1F * ChargeSpectrum(int, float, float, float=0, float=300, float=-50, float=600, int=750)
Histogram of the "charge" spectrum for one channel.
Definition ReadRun.cc:1900
float DAQ_factor
DAQ conversion factor for wavecatcher output to mV.
Definition ReadRun.h:299
void CorrectBaselineMin(vector< float >, double=0, int=2)
Baseline correction using minimum sum ( mean) in range for correction.
Definition ReadRun.cc:940
int nevents
Number of triggered events in data.
Definition ReadRun.h:282
float IndexToTime(int)
Convert the bin number of the waveform to the time of the left bin edge
Definition ReadRun.cc:2872
vector< bool > skip_event
Stores the event numbers which should be skipped in the analysis.
Definition ReadRun.h:332
double * gety(int)
Get array of y values for a certain waveform.
Definition ReadRun.cc:2727
vector< int > active_channels
Stores the numbers of the active channels.
Definition ReadRun.h:311
vector< vector< float > > rundata
Stores waveforms.
Definition ReadRun.h:132
void CorrectBaseline_function(vector< float > &, float, float, int)
Helper function called by CorrectBaseline()
Definition ReadRun.cc:714
vector< bool > PolarityMap(bool, int)
Channel map of polarity changes during reading. For parameters see ReadFile().
Definition ReadRun.cc:2831
int TimeToIndex(float)
Convert time to the bin number of the waveform.
Definition ReadRun.cc:2865
float GetPeakIntegral(TH1F *, float, float, float, float, int=0)
Calculate the integral around a peak with several options explained in GetIntWindow().
Definition ReadRun.cc:1648
int GetEventIndex(unsigned int)
Returns index of a certain event number (if data files are read in parallel threads)
Definition ReadRun.cc:2755
TH1F * His_GetTimingCFD(int, float, float, int=-999)
Plot results of GetTimingCFD()
Definition ReadRun.cc:2404
void FilterAll(double=.3, double=.9, double=.2)
Filter all waveforms.
Definition ReadRun.cc:621
void PlotWFHeatmaps(float=-20, float=200, int=880, string="", float=0, EColorPalette=kRainBow)
Plot stacks of all non-skipped waveforms for all active channels.
Definition ReadRun.cc:521
int GetChannelIndex(int)
Match channel number (wavecatcher input channel) to channel index.
Definition ReadRun.cc:2773
float tCutg
Start of time window for baseline correction when using ReadRun::Using_BaselineCorrection_in_file_loo...
Definition ReadRun.h:368
float * ChargeList(int, float=20, float=80, float=0, float=300, bool=1)
Returns array with the individual "charge"/amplitude for all events of one channel.
Definition ReadRun.cc:1781
TFile * root_out
Stores results of analysis.
Definition ReadRun.h:360
int nChannelsWC
Wavecatcher hardware max. number of channels (reduce/increase if not using the 64 channel crate)
Definition ReadRun.h:303
void PrintBaselineCorrectionResults(float=-5, float=5, int=200)
Print histogram of the baseline correction values for all channels.
Definition ReadRun.cc:1156
virtual void ReadFile(string, bool=false, int=9, string="out.root", bool=false, long long=-1)
Routine to read files created by the wavecatcher.
Definition ReadRun.cc:67
virtual int GetWaveformIndex(int, int)
Returns index of a certain event number (if data files are read in parallel threads)
Definition ReadRun.cc:2748
vector< float > PrintChargeSpectrum_pars
Starting values of the fit parameters for PrintChargeSpectrum()
Definition ReadRun.h:187
bool Using_BaselineCorrection_in_file_loop
Set true for baseline correction during data reading. Needs to be called before ReadFile().
Definition ReadRun.h:366
void PrintWFProjection(float=0, float=320, float=-50, float=50, int=200)
Plots waveform projection histograms of all channels.
Definition ReadRun.cc:1089
void PlotChannelAverages(bool=false)
Plot averages only of the good, corrected waveforms for each channel.
Definition ReadRun.cc:396
bool PlotChannel(int)
Check if a channel index should be plotted according to ReadRun::plot_active_channels.
Definition ReadRun.cc:2801
array< int, 3 > GetIntWindow(TH1F *, float, float, float, float, int=0)
Determine indices for integration window for peaks.
Definition ReadRun.cc:1517
int PrintWFProjection_cnt
Index for multiple executions of the same plotting function.
Definition ReadRun.h:114
int start_read_at_channel
Do analysis only for limited range of channels to reduce memory usage.
Definition ReadRun.h:277
vector< int > maxSumBin
Stores bin numbers where the sum of waveforms have their maximum.
Definition ReadRun.h:308
void PrintChargeSpectrumWF(float, float, float=0, float=300, int=1, float=0., float=0., float=0., float=0.)
Plot waveforms of all channels for a given event number and add the determined integration windows to...
Definition ReadRun.cc:1695
vector< float > mean_integral
Stores the mean integral/lightyield from PrintChargeSpectrum() for all channels.
Definition ReadRun.h:327
bool discard_original_eventnr
Can be used to discard the original event numbering of the data.
Definition ReadRun.h:266
virtual int Nevents_good(int=0)
Number of good events that are not skipped.
Definition ReadRun.cc:1493
vector< vector< float > > amplValuessum
Collects sums of all waveforms for each channel.
Definition ReadRun.h:135
vector< int > switch_polarity_for_channels
Stores the channel number where the polarity should be inverted. Example use to switch polarity for c...
Definition ReadRun.h:321
TH1F * WFProjectionChannel(int, int=0, int=1024, float=-50, float=50, int=200)
Waveform projections for one channel.
Definition ReadRun.cc:1061
TGraph2D * MaxDist(int, float=0, float=300)
Finds maximum amplitude for a given channel in time window ["from", "to"] and creates 3d map of wavef...
Definition ReadRun.cc:2344
int nwf
Total number of waveforms read from data: number of active channels x number of events.
Definition ReadRun.h:286
TH1F * TimeDist(int, float=0, float=300, float=0, float=300, int=100, int=0, float=.3)
Time distribution of maximum, CFD, or 10% - 90% rise time in a certain time window.
Definition ReadRun.cc:2245
void GetTimingCFD(float=.3, float=100, float=140, double=0., bool=true, int=2, bool=false, bool=false)
Determine the timing of the maximum peak with constant fraction discrimination.
Definition ReadRun.cc:1208
void IntegralFilter(vector< float >, vector< bool >, float, float, float=50, float=250, bool=false, bool=false)
Skip events with threshold on integral.
Definition ReadRun.cc:1411
TH1F * BaselineCorrectionResults(int, int, float=-5, float=5, int=200)
Histograms of the contents of baseline_correction_result.
Definition ReadRun.cc:1128
void PrintDCR(float=15, float=85, float=0, float=300, double=3)
Calculate (SiPM) dark count rate.
Definition ReadRun.cc:2221
int PrintChargeSpectrum_cnt
Index for multiple executions of the same plotting function.
Definition ReadRun.h:110
void SmoothAll(double, int)
Smoothing all waveforms which are not skipped (for testing, careful when using for analysis!...
Definition ReadRun.cc:577
void CorrectBaseline(float, float=-999)
Baseline correction constant window.
Definition ReadRun.cc:689
int binNumber
Number of bins (usually 1024 samples per waveform).
Definition ReadRun.h:301
vector< vector< float > > PrintChargeSpectrum_cal
Calibration values to normalize charge spectrum to number of photoelectrons Chennels must be ordered ...
Definition ReadRun.h:375
vector< vector< float > > timing_results
Matrix to store timing of peaks from GetTimingCFD()
Definition ReadRun.h:355
void SaveChargeLists(float=20, float=80, float=0, float=300, bool=1)
Saves TGraphs to root file with the individual "charge"/amplitude for all events and all channels.
Definition ReadRun.cc:1808
vector< vector< float > > baseline_correction_result
Stores baseline correction results for CorrectBaseline() and related functions.
Definition ReadRun.h:342
float SP_inv
1/SP
Definition ReadRun.h:295
TH1F * His_GetTimingCFD_diff(vector< int >, vector< int >, float, float, int=-999)
Plot timing difference between the mean timings of two channel ranges.
Definition ReadRun.cc:2479
void SkipEventsTimeDiffCut(int, int, double, double, bool=false)
Skip events where the time difference between two channels is outside of specified range.
Definition ReadRun.cc:1305
int PlotWFHeatmaps_cnt
Index for multiple executions of the same plotting function.
Definition ReadRun.h:116
void SkipEventsPerChannel(vector< float >, float=0, float=0, bool=false)
Skip events above/below individual thresholds per channel.
Definition ReadRun.cc:1351
string data_path
Path to data.
Definition ReadRun.h:248
void PrintMaxDist(float=0, float=300)
Finds maximum amplitude for a given channel in time window ["from", "to"] and creates 3d map of wavef...
Definition ReadRun.cc:2381
void PrintChargeSpectrum(float, float, float=0, float=300, float=-50, float=600, int=750, float=0., float=0., int=99, int=0, bool=false)
Plots the "charge" spectrums of all channels.
Definition ReadRun.cc:1971
virtual void checkData(bool isBaselineCorrection=false) const
Primitive check to see if data has been loaded.
Definition ReadRun.h:119
int end_read_at_channel
See ReadRun::start_read_at_channel.
Definition ReadRun.h:279
int nchannels
Number of active channels in data.
Definition ReadRun.h:284
void UnskipAll()
Sets skip_event flag to false for all events, removing any previous cuts.
Definition ReadRun.cc:1476
float SP
Sampling: ns per bin of data, sampling rate 3.2 GS/s -> 0.3125 ns.
Definition ReadRun.h:293
vector< TFitResultPtr > timing_fit_results
Stores the fit results of Print_GetTimingCFD() for all channels.
Definition ReadRun.h:357
T * getx(double shift=0.)
Get array of x axis (time of the bin centers) for standard wavecatcher settings.
Definition ReadRun.cc:2713
virtual int GetCurrentEvent(int)
Get the current event index for a certain waveform index.
Definition ReadRun.cc:2795
virtual ~ReadRun()
Destructor.
Definition ReadRun.cc:319
virtual bool SkipEvent(int, int=-1)
Check if event should be skipped.
Definition ReadRun.cc:1485
void PlotChannelSums(bool=false, bool=false, double=0., double=0., int=2)
Plot sums of all raw waveforms for each channel.
Definition ReadRun.cc:338
int CheckBoundsX(int)
Check if index exists in time of waveforms.
Definition ReadRun.cc:2858
void PrintTimeDist(float=0, float=300, float=0, float=300, int=100, int=0, float=.3)
Time distribution of maximum, CFD, or 10% - 90% rise time in a certain time window.
Definition ReadRun.cc:2308
vector< int > plot_active_channels
Stores the numbers of the active channels which should be plotted.
Definition ReadRun.h:316
void ShiftAllToAverageCF()
This function shifts all waveforms to the average signal starting times for each channel.
Definition ReadRun.cc:641
vector< TFitResultPtr > fit_results
Stores the fit results of PrintChargeSpectrum() for all channels and all function calls in ascending ...
Definition ReadRun.h:324
virtual int GetCurrentChannel(int)
Get the current channel index for a certain waveform index.
Definition ReadRun.cc:2788
void PrintSkippedEvents()
Prints a list of all skipped events into the terminal for diagnostics.
Definition ReadRun.cc:1461
int LastBinFileToRead
Number of last .bin file to be read in.
Definition ReadRun.h:253
void CorrectBaselineMinSlopeRMS(vector< float >, double=0, int=2)
Baseline correction method searching for non-monotonic, rather constant regions.
Definition ReadRun.cc:768
TH2F * WFHeatmapChannel(int, float=-20, float=200, int=880)
2D histogram of all non-skipped waveforms for one channel
Definition ReadRun.cc:464
void Print_GetTimingCFD_diff(vector< int >, vector< int >, float=100, float=140, int=0, int=-999, float=-999, float=-999, string="RS", bool=true)
Plot timing difference between the mean timings of two channel ranges.
Definition ReadRun.cc:2574
vector< unsigned int > eventnr_storage
Events will be stored here in the order they have been read.
Definition ReadRun.h:138