ReadRun.h

Go to the documentation of this file.

 
#ifndef _ReadRun
#define _ReadRun
 
// contact: christian.scharf@physik.hu-berlin.de
 
#include <Math/Minimizer.h>
#include <Math/Factory.h>
#include <Math/Functor.h>
#include <TClonesArray.h>
#include <TObjString.h>
#include <TObject.h>
#include <TLine.h>
#include <TGraph.h>
#include <TGraph2D.h>
#include <TMultiGraph.h>
#include <TString.h>
#include <TCanvas.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1D.h>
#include <TMath.h>
#include <TF1.h>
#include <TStyle.h>
#include <TH1F.h>
#include <TH2D.h>
#include <TEfficiency.h>
#include <TLegend.h>
#include <THStack.h>
#include <THistPainter.h>
#include <TText.h>
#include <TFitResultPtr.h>
#include <TFitResult.h>
#include <TSpline.h>
#include <TPaveStats.h> 
#include <TError.h>      // root verbosity level
#include <TSystem.h>
#include <TROOT.h>
#include <TLatex.h>
 
//C, C++
#include <stdio.h>
#include <future>
#include <assert.h>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include <stdlib.h>
#include <sstream>
#include <stdexcept>
//multithreading (avoids finding omp.h manually)
#ifndef OMP_FUNC
extern "C" {
    int omp_get_max_threads();
    int omp_get_thread_num();
}
#endif
 
#include "utils/FitFunctions.h"
#include "utils/Helpers.h"
#include "utils/Filters.h"
 
using namespace std;
 
class ReadRun {
private:
#pragma pack(1) // byte padding suppression for WaveCatcher data format
    // structs copied from
    // WaveCatcher binary -> root converter
    // by manu chauveau@cenbg.in2p3.fr
    struct event_data
    {
        int EventNumber;
        double EpochTime;
        unsigned int Year;
        unsigned int Month;
        unsigned int Day;
        unsigned int Hour;
        unsigned int Minute;
        unsigned int Second;
        unsigned int Millisecond;
        unsigned long long int TDCsamIndex;
        int nchannelstored;
    };
 
    struct channel_data_without_measurement
    {
        int channel;
        int EventIDsamIndex;
        int FirstCellToPlotsamIndex;
    };
 
    struct channel_data_with_measurement
    {
        int channel;
        int EventIDsamIndex;
        int FirstCellToPlotsamIndex;
        float   MeasuredBaseline;
        float   AmplitudeValue;
        float   ComputedCharge;
        float   RiseTimeInstant;
        float   FallTimeInstant;
        float   RawTriggerRate;
    };
#pragma pack() // byte padding suppression for WaveCatcher data format
 
protected:
    int PrintChargeSpectrum_cnt;
    int PlotChannelAverages_cnt;
    int PrintWFProjection_cnt;
    int PlotWFHeatmaps_cnt;
    
    virtual void checkData(bool isBaselineCorrection = false) const {
        (void)isBaselineCorrection;
        if (eventnr_storage.empty()) {
            throw runtime_error(
                "Error: No data has been loaded yet!\n \
                Please call ReadFile() before calling any functions which manipulate data.\n \
                Aborting execution."
            );
        }
    }
 
public:
    vector<vector<float>> rundata;
 
    vector<vector<float>> amplValuessum;
 
    vector<unsigned int> eventnr_storage;
 
    // plots amplValuessum
    void PlotChannelSums(bool = false, bool = false, double = 0., double = 0., int = 2);
 
    void PlotChannelAverages(bool = false);
 
    TH2F* WFHeatmapChannel(int, float = -20, float = 200, int = 880);
    void PlotWFHeatmaps(float = -20, float = 200, int = 880, string = "", float = 0, EColorPalette = kRainBow);
    
    // baseline correction (shifts all waveforms individually in y)
    void CorrectBaseline(float, float = -999);
    void CorrectBaseline_function(vector<float>&, float, float, int);
 
    void CorrectBaselineMinSlopeRMS(vector<float>, double = 0, int = 2);
    void CorrectBaselineMinSlopeRMS(vector<float>, double, int, int);
    void CorrectBaselineMinSlopeRMS(int = 100, bool = false, double = 0.5, int = 0, int = 0, int = 2);
 
    void CorrectBaselineMin(vector<float>, double = 0, int = 2);
    void CorrectBaselineMin(vector<float>, double, int, int);
    void CorrectBaselineMin(int = 100, double = 0.5, int = 0, int = 0, int = 2);
 
    // functions to check baseline correction results
    TH1F* WFProjectionChannel(int, int = 0, int = 1024, float = -50, float = 50, int = 200);
    void PrintWFProjection(float = 0, float = 320, float = -50, float = 50, int = 200);
 
    TH1F* BaselineCorrectionResults(int, int, float = -5, float = 5, int = 200);
    void PrintBaselineCorrectionResults(float = -5, float = 5, int = 200);
 
    // get timing of peaks
    void GetTimingCFD(float = .3, float = 100, float = 140, double = 0., bool = true, int = 2, bool = false, bool = false);
    void SkipEventsTimeDiffCut(int, int, double, double, bool = false);
 
    // average all waveforms to simplify peak ID
    void SmoothAll(double, int);
    void SmoothAll(double, string = "Gaus");
    void FilterAll(double = .3, double = .9, double = .2);
    void ShiftAllToAverageCF();
 
    // functions for charge spectrum
    array<int, 3> GetIntWindow(TH1F*, float, float, float, float, int = 0);
    array<int, 3> GetIntWindow(const vector<float>&, float, float, float, float, int = 0);
    array<int, 3> GetIntWindow(const vector<float>&, int, int, int, int, int = 0);
    float GetPeakIntegral(TH1F*, float, float, float, float, int = 0);
    float GetPeakIntegral(const vector<float>&, float, float, float, float, int = 0);
    void PrintChargeSpectrumWF(float, float, float = 0, float = 300, int = 1, float = 0., float = 0., float = 0., float = 0.);
    TH1F* ChargeSpectrum(int, float, float, float = 0, float = 300, float = -50, float = 600, int = 750);
    void PrintChargeSpectrum(float, float, float = 0, float = 300, float = -50, float = 600, int = 750, float = 0., float = 0., int = 99, int = 0, bool = false);
    vector<float> PrintChargeSpectrum_pars;
 
 
    float* ChargeList(int, float = 20, float = 80, float = 0, float = 300, bool = 1);
    void SaveChargeLists(float = 20, float = 80, float = 0, float = 300, bool = 1);
    void ChargeCorrelation(float, float, float, float, float, float, int, int, int, bool = false);
 
    // SiPM specific
    void PrintDCR(float = 15, float = 85, float = 0, float = 300, double = 3); // based on PMT::PrintChargeSpectrumPMTthreshold()
 
    // functions for time distribution
    TH1F* TimeDist(int, float = 0, float = 300, float = 0, float = 300, int = 100, int = 0, float = .3);
    void PrintTimeDist(float = 0, float = 300, float = 0, float = 300, int = 100, int = 0, float = .3);
    TGraph2D* MaxDist(int, float = 0, float = 300);
    void PrintMaxDist(float = 0, float = 300);
 
    TH1F* His_GetTimingCFD(int, float, float, int = -999);
    void Print_GetTimingCFD(float = 100, float = 140, int = 0, int = -999, string = "S", bool = true);
    TH1F* His_GetTimingCFD_diff(vector<int>, vector<int>, float, float, int = -999);
    void Print_GetTimingCFD_diff(vector<int>, vector<int>, float = 100, float = 140, int = 0, int = -999, float = -999, float = -999, string = "RS", bool = true);
 
    
    // helper functions
    TH1F* Getwf(int);                           // waveform number
    virtual TH1F* Getwf(int, int, int = 1);     // channel, eventnr, color
    template<typename T>
    T* getx(double shift = 0.);                 // x values
    double* gety(int);                          // y values for waveform index
    double* gety(int, int);                     // y values for waveform(channel, event)
 
    static pair<float, bool> LinearInterpolation(float, float, float, float, float, bool = false); // linear interpolation
    
    virtual int GetWaveformIndex(int, int); // get index of waveform from event and channel
    int GetEventIndex(unsigned int);        // get index of a triggered event (finds the correct event if files are not read sequentially)
    int GetChannelIndex(int);               // get index of a certain channel
    virtual int GetCurrentChannel(int);     // get index of channel for a certain waveform
    virtual int GetCurrentEvent(int);       // get index of event for a certain waveform
    
    bool PlotChannel(int);                  // check if channel should be plotted
    
    vector<bool> PolarityMap(bool, int);    // check if channel should be inverted
 
    int CheckBoundsX(int);                  // Check if index exists in time of waveforms
    int TimeToIndex(float);                 // Convert time to the bin number of the waveform
    float IndexToTime(int);                 // Convert the bin number of the waveform to the time
 
    ReadRun(int last_bin_file = 0, int first_bin_file = 0);
    
    virtual void ReadFile(string, bool = false, int = 9, string = "out.root", bool = false, long long = -1);
 
    virtual ~ReadRun();
 
    string data_path;
 
    int LastBinFileToRead;
    int FirstBinFileToRead;
 
 
    bool discard_original_eventnr = false;
 
    int start_read_at_channel = -1;
    int end_read_at_channel = -1;
 
    int nevents = 0;
    int nchannels;
    int nwf;
    int maxNWF = 1e7;
 
    float SP = .3125;
    float SP_inv = 1 / .3125;
    float DAQ_factor = 250. / 4096.;
    int binNumber;
    int nChannelsWC = 64;
 
    vector<int> maxSumBin;
 
    vector<int> active_channels;
    vector<int> plot_active_channels;
    
    vector<int> switch_polarity_for_channels;
 
    vector<TFitResultPtr> fit_results;
 
    vector<float> mean_integral;
 
    vector<bool> skip_event;
    virtual int Nevents_good(int = 0);
 
    void SkipEventsPerChannel(vector<float>, float = 0, float = 0, bool = false);  // in case you want to have indiviual thresholds in individual channels
    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
    void PrintSkippedEvents();
    void UnskipAll();
    virtual bool SkipEvent(int, int = -1);
 
    vector<vector<float>> baseline_correction_result;
 
    vector<vector<float>> timing_results;
    vector<TFitResultPtr> timing_fit_results;
 
    TFile* root_out;
 
    //other controls 
 
    bool Using_BaselineCorrection_in_file_loop = false;
    float tCutg;
    float tCutEndg;
 
    vector<vector<float>> PrintChargeSpectrum_cal = vector(nChannelsWC, vector<float>(2, 1));
};
#endif