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 <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 <assert.h>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include <stdlib.h>
#include <sstream>
#include <stdexcept>
 
#include "utils/FitFunctions.h"
#include "utils/Helpers.h"
#include "utils/Filters.h"
 
using namespace std;
 
class ReadRun {
private:
    int PrintChargeSpectrum_cnt;
    int PlotChannelAverages_cnt;
    int PrintWFProjection_cnt;
    int PlotWFHeatmaps_cnt;
 
 
#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;
        short waveform[1024];
    };
 
    struct channel_data_with_measurement
    {
        int channel;
        int EventIDsamIndex;
        int FirstCellToPlotsamIndex;
        float   MeasuredBaseline;
        float   AmplitudeValue;
        float   ComputedCharge;
        float   RiseTimeInstant;
        float   FallTimeInstant;
        float   RawTriggerRate;
        short waveform[1024];
    };
#pragma pack() // byte padding suppression for WaveCatcher data format
 
protected:
    void checkData() const {
        if (eventnr_storage.empty()) {
            throw std::runtime_error(
                "Error: No data has been loaded yet!\n \
                Please call ReadFile() before calling any functions which manipulate data.\n \
                Aborting execution.");
        }
    }
 
public:
    TClonesArray* rundata;
 
    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 = kGistEarth);
    
    // baseline correction (shifts all waveforms individually in y)
    void CorrectBaseline(float, float = -999);
    void CorrectBaseline_function(TH1F*, float, float, int);
 
    void CorrectBaselineMinSlopeRMS(vector<float>, double = 0, int = 2, int = 3);
    void CorrectBaselineMinSlopeRMS(int = 100, bool = false, double = 0.5, int = 0, int = 0, int = 2);
 
    void CorrectBaselineMin(vector<float>, double = 0, int = 2, int = 2);
    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);
    void SkipEventsTimeDiffCut(int, int, double, double, bool = false);
 
    void FractionEventsAboveThreshold(float = 4, bool = true, bool = true, double = 0., double = 0., 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
    int* GetIntWindow(TH1F*, float, float, float, float, int = 0);
    float GetPeakIntegral(TH1F*, 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
    TH1F* Getwf(int, int, int = 1);             // channel, eventnr, color
    double* getx(double = 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); // linear interpolation
    
    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
    int GetCurrentChannel(int);     // get index of channel for a certain waveform
    int GetCurrentEvent(int);       // get index of event for a certain waveform
    
    bool PlotChannel(int);          // check if channel should be plotted
    
 
    ReadRun(int max_no_of_bin_files_to_read = 0, int min_no_of_bin_files_to_read = 0);
    
    void ReadFile(string, bool = false, int = 9, string = "out.root", bool = false);
 
    virtual ~ReadRun();
 
    string data_path;
 
    int MaxNoOfBinFilesToRead;
    int MinNoOfBinFilesToRead;
 
 
    bool discard_original_eventnr = false;
 
    int start_read_at_channel = -1;
    int end_read_at_channel = -1;
 
    int nevents;
    int nchannels;
    int nwf;
    int maxNWF = 1e7;
 
    float SP = .3125;
    float DAQ_factor = 250. / 4096.;
    int binNumber = 1024;
    int nChannelsWC = 64;
 
    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;
    int Nevents_good();
 
    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();
 
    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;
 
    bool Shift_WFs_in_file_loop = false;
    float tWF_CF = 0.3;
    int tWF_CF_bin = 375;
    int tWF_CF_lo = 320;
    int tWF_CF_hi = 500;
 
    vector<vector<float>> PrintChargeSpectrum_cal = vector(nChannelsWC, vector<float>(2, 1));
};
#endif