ReadRun.cc

Go to the documentation of this file.

 
 
#include "ReadRun.h"
 
ReadRun::ReadRun(int max_no_of_bin_files_to_read, int min_no_of_bin_files_to_read) {
 
    cout << "\ninitializing ..." << endl;
 
    if (max_no_of_bin_files_to_read > 0) {
        cout << "will read " << max_no_of_bin_files_to_read - min_no_of_bin_files_to_read << " .bin files from file number " 
            << min_no_of_bin_files_to_read << " to file number " << max_no_of_bin_files_to_read << endl;
    }
    if (min_no_of_bin_files_to_read > 0) discard_original_eventnr = true;
 
    ROOT::EnableImplicitMT();
    TH1::AddDirectory(kFALSE);
    // init counters
    nwf = 0;
    PrintChargeSpectrum_cnt = 0;
    PlotChannelAverages_cnt = 0;
    PrintWFProjection_cnt = 0;
    PlotWFHeatmaps_cnt = 0;
    MaxNoOfBinFilesToRead = max_no_of_bin_files_to_read;
    MinNoOfBinFilesToRead = min_no_of_bin_files_to_read;
 
    root_out = new TFile(); // init results file
}
 
void ReadRun::ReadFile(string path, bool change_polarity, int change_sign_from_to_ch_num, string out_file_name, bool debug) {
 
    if (path.back() != '/') path += '/'; // add '/' to path if not present
    // save output path
    data_path = path;
 
    if (out_file_name == "") out_file_name = "out.root";
    printf("+++ saving analysis results in '%s' ...\n\n", out_file_name.c_str());
    root_out = TFile::Open(out_file_name.c_str(), "recreate");
 
    rundata = new TClonesArray("TH1F", maxNWF); //raw data will be stored here as TH1F
    rundata->BypassStreamer(kFALSE);            //kFALSE: potentially faster read & write
    TClonesArray& testrundata = *rundata;
 
    // verbosity
    bool debug_header = debug;
    bool debug_data = debug;
 
    unsigned short output_channel;
    unsigned int output_event;
    unsigned short output_nbchannels;
    unsigned short read_channels = 0;
 
    size_t length_of_waveform = static_cast<size_t>(binNumber) * sizeof(short);
 
    amplValuessum = new float* [nChannelsWC]; //sum of all wf for each channel
    for (int i = 0; i < nChannelsWC; i++) {//init
        amplValuessum[i] = new float[binNumber]();
    }
 
    maxSumBin = new int[nChannelsWC];
 
    //Start reading the raw data from .bin files.
    stringstream inFileList;
    inFileList << Helpers::ListFiles(path.c_str(), ".bin"); //all *.bin* files in folder path
    if (debug) cout << inFileList.str() << endl;
    string fileName;
    int file_counter = -1;
    int wfcounter = 0;
    int event_counter = 0;
 
    while (inFileList >> fileName) {
        // file loop
        file_counter++;
        // read only fraction/batch of the .bin files for testing or to reduce memory usage
        if (MinNoOfBinFilesToRead > 0 && MinNoOfBinFilesToRead < MaxNoOfBinFilesToRead && file_counter < MinNoOfBinFilesToRead) continue;
        if (MaxNoOfBinFilesToRead > 0 && file_counter >= MaxNoOfBinFilesToRead) break;
 
        fileName = path + fileName;
        ifstream input_file(fileName.c_str(), ios::binary | ios::in);
 
        bool has_measurement = false;
 
        if (!input_file.is_open()) {
            printf("*** failed to open '%s'\n", fileName.c_str());
            continue;
        }
 
        if (file_counter < 10 || file_counter % 10 == 0 || debug) printf("+++ reading '%s' ...\n", fileName.c_str());
 
        // Header
        string header_line;
        // HEADER 1 //
        //
        // "=== DATA FILE SAVED WITH SOFTWARE VERSION: V?.??.? ==="
        //
        getline(input_file, header_line, '\n');
 
        if (debug_header) printf("%s\n", header_line.data());
        assert(header_line[0] == '=');
 
        size_t header_version_first = header_line.find_last_of('V');
        size_t header_version_last = header_line.find_first_of(' ', header_version_first);
        string software_version = header_line.substr(header_version_first, header_version_last - header_version_first);
        if (debug_header) printf("    |- data version = '%s'\n", software_version.data());
        // convert software version
        software_version.erase(0, 1);
        int v_major, v_minor, v_patch;
        istringstream software_version_iss(software_version);
        char dot_;
        software_version_iss >> v_major >> dot_ >> v_minor >> dot_ >> v_patch;
 
        // HEADER 2 //
        // "=== WAVECATCHER SYSTEM OF TYPE ?? WITH ?? CHANNELS AND GAIN: ??? ==="
        getline(input_file, header_line, '\n');
 
        if (debug_header) printf("%s\n", header_line.data());
        assert(header_line[0] == '=');
 
        // HEADER 3 //
        // === Rate coincidence masks ... === Posttrig in ns for SamBlock ... ===
        getline(input_file, header_line, '\n');
 
        if (debug_header) printf("%s\n", header_line.data());
        assert(header_line[0] == '=');
 
        // HEADER 4 //
        // V2.9.13: === DATA SAMPLES [1024] in Volts == NB OF CHANNELS ACQUIRED: 64 == Sampling Period: 312.5 ps == INL Correction: 1
        // V2.9.15: === DATA SAMPLES [1024] in Volts == NB OF CHANNELS ACQUIRED: 64 == Sampling Period: 312.5 ps == INL Correction: 1 == MEASUREMENTS: 0 ===
        getline(input_file, header_line, '\n');
 
        if (debug_header) printf("%s\n", header_line.data());
        assert(header_line[0] == '=');
 
        size_t nsamples_first = 1 + header_line.find_last_of('[');
        size_t nsamples_last = header_line.find_first_of(']', nsamples_first);
        string nsamples_str = header_line.substr(nsamples_first, nsamples_last - nsamples_first);
 
        binNumber = atoi(nsamples_str.data());
        if (debug_header) printf("    |- data sample  = %d\n", binNumber);
        if (binNumber != 1024) cout << "\nERROR: Measurement must have 1024 samples.\nCheck WaveCatcher setting!" << endl;
 
        size_t nchannels_first = 10 + header_line.find("ACQUIRED: ", nsamples_first);
        size_t nchannels_last = header_line.find_first_of(' ', nchannels_first);
        string nchannels_str = header_line.substr(nchannels_first, nchannels_last - nchannels_first);
 
        nchannels = atoi(nchannels_str.data());
        if (debug_header) printf("    |- nchannels    = %d\n", nchannels);
 
        size_t sp_first = 8 + header_line.find("Period:");
        size_t sp_last = header_line.find(" ps");
        float sampling_period = atof(header_line.substr(sp_first, sp_last - sp_first).data());
        if (debug_header) printf("sampling period  = %f ps\n", sampling_period);
        SP = sampling_period * 1e-3;
 
        // compatibility with older WC software versions
        if (v_major ==2 && v_minor == 9 && v_patch <= 13) {
            // V2.9.13 has always measurement stored (everything is set to 0 when disabled!)
            has_measurement = true; 
        }
        else {
            size_t has_measurement_first = 14 + header_line.find("MEASUREMENTS: ", nsamples_first);
            size_t has_measurement_last = header_line.find_first_of(' ', has_measurement_first);
            string has_measurement_str = header_line.substr(has_measurement_first, has_measurement_last - has_measurement_first);
            has_measurement = atoi(has_measurement_str.data());
        }
 
        if (debug_header) printf("    `- measurement  = %d\n", has_measurement);
 
        // end of header reader
 
        event_data an_event;
 
        while (input_file.read((char*)(&an_event), sizeof(an_event))) {
            //event loop
            if (debug_data) printf("%03d has %d channels\n", an_event.EventNumber, an_event.nchannelstored);
 
            output_event = an_event.EventNumber;
            output_nbchannels = an_event.nchannelstored;
 
            if (debug_data && output_event % 200 == 0) printf("EventNr: %d, nCh: %d\n", output_event, output_nbchannels);
            if (output_nbchannels > nChannelsWC) {
                cout << "ERROR:\nThe number of channels in the data is " << output_nbchannels 
                    << ", which is larger than the maximum allowed number of channels which is set to " << nChannelsWC 
                    << "\nPlease set the parameter nChannelsWC=" << output_nbchannels << endl;
            }
 
            // do analysis only for limited range of channels to reduce memory usage for large datasets with many channels and many events
            int start_at_ch = 0;
            if (start_read_at_channel < output_nbchannels && start_read_at_channel >= 0) start_at_ch = start_read_at_channel;
            int end_at_ch = output_nbchannels - 1;
            if (end_read_at_channel == -1 && start_read_at_channel != -1) end_read_at_channel = start_read_at_channel;
            else if (end_read_at_channel < output_nbchannels && end_read_at_channel >= 0) end_at_ch = end_read_at_channel;
            read_channels = end_at_ch - start_at_ch + 1;
 
            if (event_counter == 0) cout << "\nstart at ch " << start_at_ch << " end at ch " << end_at_ch << endl;
 
            for (int ch = 0; ch < output_nbchannels; ++ch) { // channel loop
                channel_data_with_measurement a_channel_data;
                channel_data_without_measurement a_channel_data_without_measurement;
 
                if (has_measurement) { // read with 'channel_data_with_measurement' struct
                    input_file.read((char*)(&a_channel_data), sizeof(channel_data_with_measurement));
                }
                else { // read with 'channel_data_without_measurement' struct
                    input_file.read((char*)(&a_channel_data_without_measurement), sizeof(channel_data_without_measurement));
 
                    // copy the content
                    a_channel_data.channel = a_channel_data_without_measurement.channel;
                    a_channel_data.EventIDsamIndex = a_channel_data_without_measurement.EventIDsamIndex;
                    a_channel_data.FirstCellToPlotsamIndex = a_channel_data_without_measurement.FirstCellToPlotsamIndex;
                    memcpy(a_channel_data.waveform, a_channel_data_without_measurement.waveform, length_of_waveform);
                }
 
 
                output_channel = a_channel_data.channel;
                if (debug_data) printf("- reading channel %d\n", output_channel);
 
                //---------------------------------------------------------------------------------------------------------------
                if (ch >= start_at_ch && ch <= end_at_ch) {
 
                    if (event_counter == 0) active_channels.push_back(static_cast<int>(output_channel));
 
                    TString name(Form("ch%02d_%05d", output_channel, an_event.EventNumber));
                    TString title(Form("ch%d, event %d;t [ns];U [mV]", output_channel, an_event.EventNumber));
                    auto hCh = (TH1F*)testrundata.ConstructedAt(wfcounter);
                    hCh->SetName(name.Data());
                    hCh->SetTitle(title.Data());
                    hCh->SetBins(binNumber, -0.5 * SP, 1023.5 * SP);
 
                    int nshift = 0;
                    if (Shift_WFs_in_file_loop) { // shift all waveforms to tWF_CF_bin
                        float max = 0.;
                        int nmax = 0;
                        float cf = tWF_CF;
                        int count_fall = 0;
 
                        // do a mini-baseline correction (needed in case a voltage offset is set in the wavecatcher)
                        short bsln = (a_channel_data.waveform[0] + a_channel_data.waveform[1] + a_channel_data.waveform[3]) / 3;
                        for (int lll = 0; lll < binNumber; lll++) a_channel_data.waveform[lll] -= bsln;
 
                        int g_max = TMath::MaxElement(binNumber, a_channel_data.waveform);
                        int g_min = TMath::Abs(TMath::MinElement(binNumber, a_channel_data.waveform));
                        int global_max = g_max < g_min ? g_min : g_max;
 
                        for (int s = tWF_CF_lo; s < tWF_CF_hi; ++s) {
                            if (max < TMath::Abs(a_channel_data.waveform[s])) {
                                max = TMath::Abs(a_channel_data.waveform[s]);
                                nmax = s;
                            }
 
                            // stop search if the current maximum is at least 0.5 * global maximum and if there are at 
                            // least three consecutive bins where the waveform amplitude is decreasing
                            if (max > .5 * global_max && s - nmax < 4 && TMath::Abs(a_channel_data.waveform[s + 2]) + TMath::Abs(a_channel_data.waveform[s + 3]) < TMath::Abs(a_channel_data.waveform[s]) + TMath::Abs(a_channel_data.waveform[s + 1])) count_fall++;
                            else count_fall = 0;
                            if (count_fall == 3) break;
                        }
                        cf *= max;
 
                        for (int s = nmax; s > tWF_CF_lo; --s) {
                            if (cf >= TMath::Abs(a_channel_data.waveform[s])) {
                                if (cf - TMath::Abs(a_channel_data.waveform[s]) < TMath::Abs(a_channel_data.waveform[s - 1]) - cf) nshift = tWF_CF_bin - s;
                                else nshift = tWF_CF_bin - s - 1;
                                break;
                            }
                        }
                    }
 
                    // loop to fill waveform histograms
                    bool change_pol_ch = false;
                    float val = 0;
                    for (int s = 0; s < binNumber; ++s) {
                        // shift all waveforms
                        int shiftind = s - nshift;
                        if (shiftind < 0) shiftind += (binNumber - 1);
                        else if (shiftind > binNumber - 1) shiftind -= (binNumber - 1);
                        val = static_cast<float>(a_channel_data.waveform[shiftind]) * DAQ_factor;
 
                        // change the polarity for certain channels since our PMTs have negative signals while our SiPMs have positive signals
                        if (s == 0 && change_polarity) {
                            if ((output_channel >= change_sign_from_to_ch_num)
                                || (change_sign_from_to_ch_num < 0 && output_channel <= abs(change_sign_from_to_ch_num))) {
                                change_pol_ch = true;
                            }
                            else if (!switch_polarity_for_channels.empty()
                                && find(switch_polarity_for_channels.begin(), switch_polarity_for_channels.end(), output_channel) 
                                != switch_polarity_for_channels.end()) {
                                change_pol_ch = true;
                            }
                        }
                        if (change_pol_ch == true) val *= -1;
 
                        // fill waveforms
                        hCh->SetBinContent(s + 1, val);
 
                        // channel sums
                        amplValuessum[ch][s] += val;
                    }
 
                    // baseline correction
                    if (Using_BaselineCorrection_in_file_loop) {
                        CorrectBaseline_function(hCh, tCutg, tCutEndg, wfcounter);
                    }
 
                    wfcounter++;
                }//--------------------------------------------------------------------------------------------------------------
 
            } // for ch
 
            skip_event.push_back(false);
            if (!discard_original_eventnr) eventnr_storage.push_back(output_event); // Stores the current WaveCatcher event number
            else eventnr_storage.push_back(event_counter);
            event_counter++;
        } // while an_event
 
        input_file.close();
    } // for file_id
 
    // in case there are empty channels, nchannels is the number of channels which contain data
    nchannels = read_channels;
 
    // get bins where the sum spectrum has its maximum for runs with fixed trigger delay and fixed 
    // integration window relative to the max of the sum spectrum (not working for DC measurement)
    for (int ch = 0; ch < nchannels; ch++) {
        float max = 0.;
        for (int i = 0; i < binNumber; i++) {
            if (amplValuessum[ch][i] > max) {
                max = amplValuessum[ch][i];
                maxSumBin[ch] = i;
            }
        }
    }
 
    nevents = event_counter;
    nwf = wfcounter;
 
    printf("Finished reading %d files containing %d events with %d channels.\n\n", file_counter, nevents, nchannels);
}
 
ReadRun::~ReadRun() {
    //delete rundata;
    //plot_active_channels.clear();
    if (root_out->IsOpen()) root_out->Close();
    cout << "\ndeleting nothing currently..." << endl;
}
 
void ReadRun::PlotChannelSums(bool smooth, bool normalize, double shift, double sigma, int smooth_method) {
 
    double* xv = getx(shift);
    auto mgsums = new TMultiGraph();
    mgsums->SetTitle("channel sums; t [ns]; amplitude [mV]");
    if (normalize) mgsums->SetTitle("channel sums; t [ns]; amplitude [arb.]");
 
    double max = 0., min = 0.;
 
    for (int i = 0; i < nchannels; i++) {
        if (PlotChannel(i)) {
            double* yv = new double[binNumber];
            for (int k = 0; k < binNumber; k++) yv[k] = static_cast<double>(amplValuessum[i][k]);
 
            if (smooth) Filters::SmoothArray(yv, binNumber, sigma, smooth_method);
 
            TGraph* gr = new TGraph(binNumber, xv, yv);
            delete[] yv;
 
            double tmp_min = TMath::MinElement(gr->GetN(), gr->GetY());
            if (tmp_min < min) min = tmp_min;
            double tmp_max = TMath::MaxElement(gr->GetN(), gr->GetY());
            if (tmp_max > max) max = tmp_max;
            if (normalize) {
                for (int j = 0; j < gr->GetN(); j++) gr->SetPointY(j, gr->GetPointY(j) / tmp_max);
            }
 
            TString name(Form("channel_%02d", active_channels[i]));
            TString title(Form("Channel %d", active_channels[i]));
            gr->SetName(name.Data());
            gr->SetTitle(title.Data());
            gr->SetLineColor(Helpers::rcolor(i));
            gr->SetMarkerColor(Helpers::rcolor(i));
            mgsums->Add(gr);
        }
    }
    delete[] xv;
 
    auto sumc = new TCanvas("Sums", "", 600, 400);
    mgsums->Draw("AL");
    if (normalize) mgsums->GetYaxis()->SetRangeUser(-0.2, 1);
    else mgsums->GetYaxis()->SetRangeUser(min, max);
    sumc->BuildLegend(0.85, 0.70, .99, .95);
    root_out->WriteObject(mgsums, "channelsums");
    root_out->WriteObject(sumc, "channelsums_c");
}
 
void ReadRun::PlotChannelAverages(bool normalize) {
 
    double* xv = getx(0);
    auto mgav = new TMultiGraph();
    mgav->SetTitle("channel averages; t [ns]; amplitude [mV]");
    if (normalize) mgav->SetTitle("channel averages; t[ns]; amplitude[arb.]");
 
    double max = 0., min = 0.;
 
    for (int i = 0; i < nchannels; i++) {
        if (PlotChannel(i)) {
 
            double* yv = new double[binNumber]();
 
            for (int j = 0; j < nevents; j++) {
                if (!skip_event[j]) {
                    double* y_tmp = gety(i, j);
                    for (int k = 0; k < binNumber; k++) yv[k] += y_tmp[k];
                    delete[] y_tmp;
                }
            }
 
            double norm = static_cast<double>(Nevents_good());
            for (int k = 0; k < binNumber; k++) yv[k] /= norm;
 
            auto gr = new TGraph(binNumber, xv, yv);
            delete[] yv;
 
            double tmp_min = TMath::MinElement(gr->GetN(), gr->GetY());
            if (tmp_min < min) min = tmp_min;
            double tmp_max = TMath::MaxElement(gr->GetN(), gr->GetY());
            if (tmp_max > max) max = tmp_max;
            if (normalize) {
                for (int j = 0; j < gr->GetN(); j++) gr->SetPointY(j, gr->GetPointY(j) / tmp_max);
            }
 
            TString name(Form("channel_%02d", active_channels[i]));
            TString title(Form("Channel %d", active_channels[i]));
            gr->SetName(name.Data());
            gr->SetTitle(title.Data());
            gr->SetLineColor(Helpers::rcolor(i));
            gr->SetMarkerColor(Helpers::rcolor(i));
            mgav->Add(gr);
        }
    }
    delete[] xv;
 
    string cname("Averages_" + to_string(PlotChannelAverages_cnt++));
    auto avc = new TCanvas(cname.c_str(), cname.c_str(), 600, 400);
    mgav->Draw("AL");
    if (normalize) mgav->GetYaxis()->SetRangeUser(-0.2, 1);
    else mgav->GetYaxis()->SetRangeUser(min, max);
    avc->BuildLegend(0.85, 0.70, .99, .95);
    root_out->WriteObject(mgav, ("channelaverages" + to_string(PlotChannelAverages_cnt)).c_str());
    root_out->WriteObject(avc, ("channelaverages_c" + to_string(PlotChannelAverages_cnt)).c_str());
}
 
 
TH2F* ReadRun::WFHeatmapChannel(int channel_index, float ymin, float ymax, int n_bins_y) {
 
    TString name(Form("channel__%02d", active_channels[channel_index]));
    TH2F* h2 = new TH2F(name.Data(), name.Data(), binNumber, 0, SP * static_cast<float>(binNumber), n_bins_y, ymin, ymax);
    h2->GetXaxis()->SetTitle("t [ns]");
    h2->GetYaxis()->SetTitle("I [arb.]");
    h2->GetZaxis()->SetTitle("entries");
 
    for (int j = 0; j < nevents; j++) {
        if (!skip_event[j]) {
            auto his = Getwf(channel_index, j);
            for (int i = 0; i < binNumber; i++) h2->Fill(his->GetBinCenter(i), his->GetBinContent(i));
        }
    }
    return h2;
}
 
void ReadRun::PlotWFHeatmaps(float ymin, float ymax, int n_bins_y, string z_opt, float z_max, EColorPalette palette) {
 
    gStyle->SetOptStat(0);
    string name("waveforms_heatmap_" + to_string(PlotWFHeatmaps_cnt++));
    auto wfhm_c = new TCanvas(name.c_str(), name.c_str(), 600, 400);
    Helpers::SplitCanvas(wfhm_c, active_channels, plot_active_channels);
 
    int current_canvas = 0;
    for (int i = 0; i < nchannels; i++) {
        if (PlotChannel(i)) {
            wfhm_c->cd(++current_canvas);
            // formatting
            gPad->SetTopMargin(.1);
            gPad->SetBottomMargin(.1);
            gPad->SetLeftMargin(.15);
            gPad->SetRightMargin(.15);
            gStyle->SetPalette(palette);
 
            auto h2 = WFHeatmapChannel(i, ymin, ymax, n_bins_y);
            h2->SetContour(99);
            h2->SetStats(0);
            if (z_opt == "COLZ0") h2->Draw("COLZ0");
            else h2->Draw("CONT4Z");
            if (z_opt == "log") {
                gPad->SetLogz();
                if (z_max > 1) h2->GetZaxis()->SetRangeUser(1, z_max);
            }
            else if (z_max > 0) h2->GetZaxis()->SetRangeUser(0, z_max);
        }
    }
    wfhm_c->Update();
 
    root_out->WriteObject(wfhm_c, name.c_str());
}
 
void ReadRun::SmoothAll(double sigma, int method) {
    cout << "\nSmoothing all non-skipped waveforms:" << endl;
    for (int j = 0; j < nwf; j++) {
        if (!skip_event[GetCurrentEvent(j)]) {
            auto his = Getwf(j);
            double* yvals = Helpers::gety(his);
            Filters::SmoothArray(yvals, binNumber, sigma, method, SP);
            for (int i = 1; i <= his->GetNbinsX(); i++) his->SetBinContent(i, yvals[i - 1]);
            delete[] yvals;
        }
        Helpers::PrintProgressBar(j, nwf);
    }
}
 
void ReadRun::SmoothAll(double sigma, string method) {
    cout << "\nSmoothing all non-skipped waveforms:" << endl;
    for (int j = 0; j < nwf; j++) {
        if (!skip_event[GetCurrentEvent(j)]) {
            auto his = Getwf(j);
            double* yvals = Helpers::gety(his);
            Filters::SmoothArray(yvals, binNumber, sigma, method, SP);
            for (int i = 1; i <= binNumber; i++) his->SetBinContent(i, yvals[i - 1]);
            delete[] yvals;
        }
        Helpers::PrintProgressBar(j, nwf);
    }
}
 
void ReadRun::FilterAll(double sigma1, double sigma2, double factor) {
    cout << "\nFiltering all waveforms" << endl;
    for (int j = 0; j < nwf; j++) {
        auto his = Getwf(j);
        double* yvals = Helpers::gety(his);
        if (factor > 0 && factor <= 1) Filters::ResponseFilter(yvals, binNumber, sigma1, sigma2, factor, SP);
        else Filters::SecondOrderUnderdampedFilter(yvals, binNumber, sigma1, sigma2, abs(factor), SP);
        for (int i = 1; i <= binNumber; i++) his->SetBinContent(i, yvals[i - 1]);
        delete[] yvals;
        Helpers::PrintProgressBar(j, nwf);
    }
}
 
void ReadRun::ShiftAllToAverageCF() {
    cout << "\nShifting all waveforms to the average constant fraction time for each channel:" << endl;
 
    //call GetTimingCFD() in case it was not initialized
    if (static_cast<int>(timing_results.size()) == 0) GetTimingCFD();
 
    double* timing_mean = new double[nchannels]();
 
    for (int j = 0; j < nwf; j++) {
        if (!skip_event[GetCurrentEvent(j)]) timing_mean[GetCurrentChannel(j)] += timing_results[j][0];
    }
 
    double norm = static_cast<double>(Nevents_good());
 
    int* timing_mean_n = new int[nchannels];
    for (int i = 0; i < nchannels; i++) {
        timing_mean_n[i] = static_cast<int>(round(timing_mean[i] / norm));
    }
    delete[] timing_mean;
 
    for (int j = 0; j < nwf; j++) {
        if (!skip_event[GetCurrentEvent(j)]) {
            int shift = static_cast<int>(timing_results[j][0]) - timing_mean_n[GetCurrentChannel(j)];
            auto his = Getwf(j);
            Helpers::ShiftTH1F(his, shift);
        }
        Helpers::PrintProgressBar(j, nwf);
    }
    delete[] timing_mean_n;
}
 
void ReadRun::CorrectBaseline(float tCut, float tCutEnd) {
    checkData();
 
    cout << "\nPerforming simple baseline correction in fixed time window."
        << "This method is only suitable for measurements without dark counts!" << endl;
    tCutg = tCut;
    tCutEndg = tCutEnd;
    if (nwf == 0) {
        Using_BaselineCorrection_in_file_loop = true;
    }
    else {
        cout << "Baseline correction (" << nwf << " waveforms):" << endl;
        for (int j = 0; j < nwf; j++) {
            CorrectBaseline_function(Getwf(j), tCut, tCutEnd, j);
            Helpers::PrintProgressBar(j, nwf);
        }
    }
}
 
void ReadRun::CorrectBaseline_function(TH1F* his, float tCut, float tCutEnd, int nwaveform) {
    int iCut, iCutEnd;
    float corr = 0;
 
    iCut = his->GetXaxis()->FindBin(tCut);
 
    if (tCutEnd <= 0) { //
        corr = his->Integral(1, iCut) / static_cast<float>(iCut);
    }
    else {
        iCutEnd = his->GetXaxis()->FindBin(tCutEnd);
        corr = his->Integral(iCut, iCutEnd) / static_cast<float>(iCutEnd - iCut + 1);
    }
 
    // write corrected values to histograms
    if (tCut >= 0) {
        for (int i = 1; i <= his->GetNbinsX(); i++) his->SetBinContent(i, his->GetBinContent(i) - corr);
    }
 
    if (!Using_BaselineCorrection_in_file_loop) {
        baseline_correction_result.push_back(vector<float>());
        baseline_correction_result[nwaveform].push_back(corr);
        baseline_correction_result[nwaveform].push_back(0);
        baseline_correction_result[nwaveform].push_back(tCut);
        baseline_correction_result[nwaveform].push_back(tCutEnd);
    }
}
 
void ReadRun::CorrectBaselineMinSlopeRMS(vector<float> window, double sigma, int smooth_method, int increment) {
    checkData();
 
    cout << "\nBaseline correction (minimum slope variation method, " << nwf << " waveforms):" << endl;
    if (window.empty()) cout << "\nWarning: Window not set in CorrectBaselineMinSlopeRMS. Will use default values." << endl;
    if (sigma != 0.) cout << "\nNotification: Using smoothing in CorrectBaselineMinSlopeRMS." << endl;
    if (increment < 1) increment = 1;
 
    int nbins_average = !window.empty() ? static_cast<int>(round(abs(window[0]) / SP)) : static_cast<int>(50. / SP);
    int start_search_at = static_cast<int>(window.size()) > 1 ? static_cast<int>(round(abs(window[1]) / SP)) : 0;
    int end_search_at = static_cast<int>(window.size()) > 2 ? 
        min(static_cast<int>(round(abs(window[2]) / SP)), binNumber - 1) : binNumber - 1;
 
    int end_search_loop_at = end_search_at - start_search_at - nbins_average;
 
    // if no valid static search window is specified, it will be dynamic from 0 ns up to 8 ns before the global maximum
    bool search_relative_to_local_max = false;
    int min_distance_from_max = static_cast<int>(8. / SP);
    if (start_search_at < 0 || end_search_loop_at < increment) {
        search_relative_to_local_max = true;
        start_search_at = 0;
        cout << "\nNotification: Using dynamic search window in CorrectBaselineMinSlopeRMS." << endl;
    }
 
    int nbins_search = end_search_at - start_search_at;
 
    float minchange = 1.e99;
    float sum = 0, change = 0, minsumsq = 0, sqsum = 0, minsqsum = 0, corr = 0;
    int iintwindowstart = 0, imax = 0;
 
    for (int j = 0; j < nwf; j++) {
        minchange = 1.e99;
        minsumsq = 0, sqsum = 0, minsqsum = 0, corr = 0;
        iintwindowstart = 0;
 
        auto his = Getwf(j);
        if (search_relative_to_local_max) {
            imax = GetIntWindow(his, 0, 0, static_cast<float>(nbins_search + min_distance_from_max) * SP,
                static_cast<float>(end_search_at) * SP)[1];
            end_search_at = imax - min_distance_from_max;
            nbins_search = end_search_at;
            end_search_loop_at = nbins_search - nbins_average;
        }
 
        double* yvals = Helpers::gety(his, start_search_at, end_search_at);
        // smoothing suppresses variations in slope due to noise, so the method is potentially more sensitive to excluding peaks
        if (sigma > 0) Filters::SmoothArray(yvals, nbins_search, sigma, smooth_method);
        //calculate slope
        double slope[nbins_search - 1];
        for (int i = 0; i < nbins_search - 1; i++) slope[i] = yvals[i + 1] - yvals[i];
        delete[] yvals;
 
        //find window for correction
        for (int i = 0; i < end_search_loop_at; i += increment) { // recommendend max. 3 bins (~1 ns) 
            sum = 0., sqsum = 0., change = 0.;
 
            for (int k = i; k < nbins_average + i; k += increment) { // recommendend max. 3 bins (~1 ns)
                sum += slope[k];
                sqsum += (slope[k] * slope[k]);
            }
            change = sqsum + sum * sum;
 
            if (change < minchange) {
                minchange = change;
                iintwindowstart = i + start_search_at;
                minsumsq = sum * sum;
                minsqsum = sqsum;
            }
        }
        // do correction
        corr = his->Integral(iintwindowstart + 1, iintwindowstart + nbins_average + 1) / static_cast<float>(nbins_average + 1);
        for (int i = 1; i <= binNumber; i++) his->SetBinContent(i, his->GetBinContent(i) - corr);
 
        baseline_correction_result.push_back(vector<float>());
        baseline_correction_result[j].push_back(corr);
        baseline_correction_result[j].push_back(minchange);
        baseline_correction_result[j].push_back(static_cast<float>(iintwindowstart) * SP);
        baseline_correction_result[j].push_back(static_cast<float>(iintwindowstart + nbins_average) * SP);
        baseline_correction_result[j].push_back(minsumsq);
        baseline_correction_result[j].push_back(minsqsum);
        Helpers::PrintProgressBar(j, nwf);
    }
}
 
void ReadRun::CorrectBaselineMinSlopeRMS(int nIntegrationWindow, bool smooth, double sigma, int max_bin_for_baseline, int start_at, int smooth_method) {
    checkData();
 
    cout << "Notification: This is a deprecated version of CorrectBaselineMinSlopeRMS. "
        << "It will be removed in future releases. Parameter bool smooth=" << smooth << " will be ignored." << endl;
    vector<float> window;
    window.push_back(static_cast<float>(nIntegrationWindow) * SP);
    window.push_back(static_cast<float>(start_at) * SP);
    window.push_back(static_cast<float>(max_bin_for_baseline) * SP);
    CorrectBaselineMinSlopeRMS(window, sigma, smooth_method, 3);
}
 
void ReadRun::CorrectBaselineMin(vector<float> window, double sigma, int smooth_method, int increment) {
    checkData();
 
    cout << "\nBaseline correction (minimal sum method, " << nwf << " waveforms):" << endl;
    if (window.empty()) cout << "\nWarning: Window not set in CorrectBaselineMin. Will use default values." << endl;
    if (sigma != 0.) cout << "\nNotification: Using smoothing in CorrectBaselineMin." << endl;
    if (increment < 1) increment = 1;
 
    int nbins_average = !window.empty() ? static_cast<int>(round(abs(window[0]) / SP)) : static_cast<int>(50. / SP);
    int start_search_at = static_cast<int>(window.size()) > 1 ? static_cast<int>(round(abs(window[1]) / SP)) : 0;
    int end_search_at = static_cast<int>(window.size()) > 2 ? 
        min(static_cast<int>(round(abs(window[2]) / SP)), binNumber - 1) : binNumber - 1;
 
    int end_search_loop_at = end_search_at - start_search_at - nbins_average;
 
    // if no valid static search window is specified, it will be dynamic from 0 ns up to 8 ns before the global maximum
    bool search_relative_to_local_max = false;
    int min_distance_from_max = static_cast<int>(8. / SP);
    if (start_search_at < 0 || end_search_loop_at < increment) {
        search_relative_to_local_max = true;
        start_search_at = 0;
        cout << "\nNotification: Using dynamic search window in CorrectBaselineMin." << endl;
    }
 
    int nbins_search = end_search_at - start_search_at;
 
    float minchange = 1e9;
    float sum = 0, corr = 0;
    int iintwindowstart = 0, imax = 0;
 
    for (int j = 0; j < nwf; j++) {
        minchange = 1.e9;
        corr = 0;
        iintwindowstart = 0;
 
        auto his = Getwf(j);
        if (search_relative_to_local_max) {
            imax = GetIntWindow(his, 0, 0, (float)(nbins_search + min_distance_from_max) * SP, (float)(end_search_at)*SP)[1];
            end_search_at = imax - min_distance_from_max;
            nbins_search = end_search_at;
            end_search_loop_at = nbins_search - nbins_average;
        }
 
        double* yvals = Helpers::gety(his, start_search_at, end_search_at);
        // smoothing suppresses variations in slope due to noise, so the method is potentially more sensitive to excluding peaks
        if (sigma > 0) Filters::SmoothArray(yvals, nbins_search, sigma, smooth_method);
        //find window for correction
        for (int i = 0; i < end_search_loop_at; i += increment) { // recommendend max. 2 bins or min. method
            sum = 0;
            for (int k = i; k < nbins_average + i; k += increment) { // recommendend max. 2 bins
                sum += yvals[k];
            }
 
            if (sum < minchange) {
                minchange = sum;
                iintwindowstart = i + start_search_at;
            }
        }
        // do correction
        corr = his->Integral(iintwindowstart + 1, iintwindowstart + nbins_average + 1) / static_cast<float>(nbins_average + 1);
        for (int i = 1; i <= binNumber; i++) his->SetBinContent(i, his->GetBinContent(i) - corr);
        delete[] yvals;
 
        baseline_correction_result.push_back(vector<float>());
        baseline_correction_result[j].push_back(corr);
        baseline_correction_result[j].push_back(minchange);
        baseline_correction_result[j].push_back(static_cast<float>(iintwindowstart) * SP);
        baseline_correction_result[j].push_back(static_cast<float>(iintwindowstart + nbins_average) * SP);
        Helpers::PrintProgressBar(j, nwf);
    }
}
 
void ReadRun::CorrectBaselineMin(int nIntegrationWindow, double sigma, int max_bin_for_baseline, int start_at, int smooth_method) {
    checkData();
 
    cout << "Notification: This is a deprecated version of CorrectBaselineMin. It will be removed in future releases." << endl;
    vector<float> window;
    window.push_back(static_cast<float>(nIntegrationWindow) * SP);
    window.push_back(static_cast<float>(start_at) * SP);
    window.push_back(static_cast<float>(max_bin_for_baseline) * SP);
    CorrectBaselineMin(window, sigma, smooth_method, 2);
}
 
TH1F* ReadRun::WFProjectionChannel(int channel_index, int from_n, int to_n, float rangestart, float rangeend, int nbins) {
 
    TString name(Form("channel__%02d", active_channels[channel_index]));
    auto h1 = new TH1F(name.Data(), name.Data(), nbins, rangestart, rangeend);
 
    for (int j = 0; j < nevents; j++) {
        if (!skip_event[j]) {
            auto his = Getwf(channel_index, j);
            for (int i = from_n; i <= to_n; i++) h1->Fill(his->GetBinContent(i));
        }
    }
    return h1;
}
 
void ReadRun::PrintWFProjection(float from, float to, float rangestart, float rangeend, int nbins) {
    gStyle->SetOptFit(111);
    float default_rangestart = -10;
    float default_rangeend = 20;
    if (default_rangestart > rangestart) default_rangestart = rangestart;
    if (default_rangeend < rangeend) default_rangeend = rangeend;
    int default_nbins = static_cast<int>((default_rangeend - default_rangestart) * nbins / (rangeend - rangestart));
 
    string ctitle("WFProjection" + to_string(PrintWFProjection_cnt++));
    auto wf_projection_c = new TCanvas(ctitle.c_str(), ctitle.c_str(), 600, 400);
    Helpers::SplitCanvas(wf_projection_c, active_channels, plot_active_channels);
    int current_canvas = 0;
 
    int from_n = (from > 0 && from < static_cast<float>(binNumber) * SP) ? static_cast<int>(round(abs(from) / SP)) + 1 : 0;
    int to_n = (to > 0 && to < static_cast<float>(binNumber) * SP) ? static_cast<int>(round(abs(to) / SP)) + 1 : binNumber;
 
    for (int i = 0; i < nchannels; i++) {
        if (PlotChannel(i)) {
            wf_projection_c->cd(++current_canvas);
 
            auto his = WFProjectionChannel(i, from_n, to_n, default_rangestart, default_rangeend, default_nbins);
            his->GetYaxis()->SetTitle("#Entries");
            his->GetXaxis()->SetTitle("amplitude in mV");
            TString name(Form("WFProjection channel_%02d_%d", active_channels[i], PrintWFProjection_cnt));
            his->Draw();
            his->Fit("gaus", "WWM", "same");
            root_out->WriteObject(his, name.Data());
        }
    }
 
    Helpers::SetRangeCanvas(wf_projection_c, rangestart, rangeend);
    root_out->WriteObject(wf_projection_c, ("WFProjections" + to_string(PrintWFProjection_cnt)).c_str());
}
 
TH1F* ReadRun::BaselineCorrectionResults(int channel_index, int which, float rangestart, float rangeend, int nbins) {
    if (baseline_correction_result.empty() || static_cast<int>(baseline_correction_result[0].size()) < which - 1) {
        cout << "\nError: baseline_correction_result empty. Call baseline correction first." << endl;
    }
    TString name(Form("channel__%02d", active_channels[channel_index]));
    auto h1 = new TH1F(name.Data(), name.Data(), nbins, rangestart, rangeend);
 
    for (int j = 0; j < nevents; j++) {
        if (!skip_event[j]) h1->Fill(baseline_correction_result[j * nchannels + channel_index][which]);
    }
    return h1;
}
 
void ReadRun::PrintBaselineCorrectionResults(float rangestart, float rangeend, int nbins) {
    gStyle->SetOptFit(111);
    float default_rangestart = -10;
    float default_rangeend = 20;
    if (default_rangestart > rangestart) default_rangestart = rangestart;
    if (default_rangeend < rangeend) default_rangeend = rangeend;
    int default_nbins = static_cast<int>((default_rangeend - default_rangestart) * nbins / (rangeend - rangestart));
 
    string ctitle;
    ctitle = "Correction values in mV";
    auto blc_res_c = new TCanvas(ctitle.c_str(), ctitle.c_str(), 600, 400);
    Helpers::SplitCanvas(blc_res_c, active_channels, plot_active_channels);
    int current_canvas = 0;
 
    for (int i = 0; i < nchannels; i++) {
        if (PlotChannel(i)) {
            blc_res_c->cd(++current_canvas);
 
            auto his = BaselineCorrectionResults(i, 0, default_rangestart, default_rangeend, default_nbins);
            his->GetYaxis()->SetTitle("#Entries");
            his->GetXaxis()->SetTitle(ctitle.c_str());
            his->Draw();
            his->Fit("gaus", "WWM", "same");
        }
    }
    Helpers::SetRangeCanvas(blc_res_c, rangestart, rangeend);
    root_out->WriteObject(blc_res_c, ctitle.c_str());
}
 
void ReadRun::GetTimingCFD(float cf_r, float start_at_t, float end_at_t, double sigma, bool find_CF_from_start, int smooth_method, bool use_spline) {
 
    int start_at = max(static_cast<int>(floor(start_at_t / SP)), 0);
    int end_at = min(static_cast<int>(ceil(end_at_t / SP)), binNumber - 1);
    int n_range = end_at - start_at;
 
    cout << "\nGet timing at " << (cf_r > 0 && cf_r <= 1 ? "CF=" : "threshold=");
    printf("%.2f between %.2f ns and %.2f ns (%d waveforms):\n", cf_r, start_at_t, end_at_t, nwf);
 
    for (int j = 0; j < nwf; j++) {
        double* yvals = Helpers::gety(Getwf(j), start_at, end_at); // get range where to search for CFD for timing
 
        if (sigma > 0.) Filters::SmoothArray(yvals, n_range, sigma, smooth_method); // smoothing to suppress noise, will also change timing so use with care!
 
        float max = 0.;
        int n_max = 0;
        int i = 0;
        for (i = 0; i < n_range; i++) {
            if (yvals[i] > max) {
                max = yvals[i];
                n_max = i;
            }
        }
 
        float cf = cf_r;
        if (cf_r > 0 && cf_r <= 1) cf *= max;
 
        if (!find_CF_from_start) {
            i = n_max;
            while (yvals[i] > cf && i >= 0) i--;
            i++;
        }
        else {
            i = 0;
            while (yvals[i] < cf && i < n_max) i++;
            i--;
        }
 
        // do interpolation for cf
        pair<float, bool> lin_interpol_res = {0, true};
        lin_interpol_res = LinearInterpolation(cf, static_cast<float>(i), static_cast<float>(i + 1), yvals[i], yvals[i + 1]);
        // go to center of bin
        float interpol_bin = lin_interpol_res.first + .5;
 
        if (use_spline) { // use spline interpolation with tolerance epsilon*bin_size
            double epsilon = 1e-4;
            double x_low = interpol_bin - .5;
            double x_high = interpol_bin + .5;
 
            double* xvals = new double[n_range];
            for (int k = 0; k < n_range; k++) xvals[k] = static_cast<double>(k) + .5;
 
            TSpline5* wfspl = 0;
            wfspl = new TSpline5("wf_spline", xvals, yvals, n_range, "b1e1b2e2", 0., 0., 0., 0.);
            delete[] xvals;
 
            // using bisection method: halving search window until cf is less than epsilon bins from spline value
            while (x_high - x_low > epsilon) {
                double x_mid = (x_low + x_high) / 2;
                double f_mid = wfspl->Eval(x_mid);
                if (f_mid == cf) break;
 
                if (f_mid > cf) x_high = x_mid;
                else x_low = x_mid;
            }
            interpol_bin = (x_low + x_high) / 2;
            delete wfspl;
        }
 
        timing_results.push_back(vector<float>());
        timing_results[j].push_back(interpol_bin);                                          // the bin we looked for
        timing_results[j].push_back((interpol_bin + static_cast<float>(start_at)) * SP);    // cfd-time we looked for
        timing_results[j].push_back(max);                                                   // maximum value
        timing_results[j].push_back(n_max);                                                 // bin of maximum
        timing_results[j].push_back(cf);                                                    // constant fraction
        timing_results[j].push_back(static_cast<float>(start_at) * SP);                     // starting time
        timing_results[j].push_back(static_cast<float>(end_at) * SP);                       // end time
        timing_results[j].push_back(static_cast<float>(lin_interpol_res.second));           // flag will be 1 if linear interpolation worked
        delete[] yvals;
 
        Helpers::PrintProgressBar(j, nwf);
    }
}
 
void ReadRun::SkipEventsTimeDiffCut(int first_channel_abs, int second_channel_abs, double time_diff_min, double time_diff_max, bool verbose) {
 
    cout << "\n Removing events if the event-wise time difference between the main peaks in ch"
        << first_channel_abs << " and ch" << second_channel_abs << " is <" << setprecision(2)
        << time_diff_min << " ns or >" << time_diff_max << " ns" << endl;
 
    int counter = 0;
    int first_channel = GetChannelIndex(first_channel_abs);
    int second_channel = GetChannelIndex(second_channel_abs);
    int timing_results_size = static_cast<int>(timing_results.size());
 
    // call GetTimingCFD() in case it was not initialized
    if (static_cast<int>(timing_results.size()) == 0) GetTimingCFD();
 
    // loop through events, calculate timing difference between channels and compare with cuts
    for (int j = 0; j < nwf; j += nchannels) {
        if (!skip_event[GetCurrentEvent(j)]) {
            float time_diff = timing_results[j + second_channel][1] - timing_results[j + first_channel][1];
 
            if (j <= timing_results_size && (time_diff < time_diff_min || time_diff > time_diff_max)) {
                int currevent = eventnr_storage[GetCurrentEvent(j)];
                if (verbose) cout << "\nevent:\t" << currevent << "\tchannels:\t" << first_channel_abs << " & " << second_channel_abs << "\ttime diff:\t" << time_diff;
                skip_event[GetCurrentEvent(j)] = true;
                counter++;
            }
        }
        Helpers::PrintProgressBar(j, nwf);
    }
    cout << "\t" << counter << " events will be cut out of " << nevents << endl;
}
 
void ReadRun::FractionEventsAboveThreshold(float threshold, bool max, bool greater, double from, double to, bool verbose) {
 
    int occurences = 0;
    int occurences2ch = 0;
    int o2ch = 0;
    int currchannel = 0;
    int currevent = 0;
    int lastevent = 0;
    if (plot_active_channels.empty()) plot_active_channels = active_channels;
    vector<int> counter_abovethr(static_cast<int>(plot_active_channels.size()));
    // DORAMAS: It stores a counter of events above threshold for each channel that will be plotted
 
    cout << "\nSearching for fraction of events with " << (max ? "max" : "min") << (greater ? " > " : " < ") << threshold << " mV:" << endl;
 
    for (int j = 0; j < nwf; j++) {
        auto his = (TH1F*)(Getwf(j))->Clone(); // use Clone() to not change ranges of original histogram
 
        // set range where to search for amplitudes above threshold
        if (from >= 0 && to > 0) {
            his->GetXaxis()->SetRange(his->GetXaxis()->FindBin(from), his->GetXaxis()->FindBin(to));
        }
 
        currchannel = j - nchannels * GetCurrentEvent(j);
        if (find(plot_active_channels.begin(), plot_active_channels.end(), active_channels[currchannel]) != plot_active_channels.end()) {
            if ((max && greater && his->GetMaximum() > threshold) || (max && !greater && his->GetMaximum() < threshold) || (!max && greater && his->GetMinimum() > threshold) || (!max && !greater && his->GetMinimum() < threshold)) {
                currevent = eventnr_storage[GetCurrentEvent(j)];
 
                if (verbose) cout << "\nevent:\t" << currevent << "\tchannel:\t" << active_channels[currchannel];
 
                // We must use 'distance' to make sure the position in 'counter_above' matches with the corresponding channel's position at 'plot_active_channels'
                counter_abovethr[distance(plot_active_channels.begin(), find(plot_active_channels.begin(), plot_active_channels.end(), active_channels[currchannel]))] += 1;
                // This is to detect events w/ at least two channels above threshold
                if (lastevent != currevent) occurences += 1;
                if (lastevent == currevent && o2ch != occurences) {
                    occurences2ch += 1;
                    o2ch = occurences;
                }
                lastevent = currevent;
            }
        }
        delete his;
        Helpers::PrintProgressBar(j, nwf);
    }
 
    //  Loop to show the fraction of events above threshold for each channel that will be plotted
    for (int i = 0; i < static_cast<int>(plot_active_channels.size()); i++) {
        cout << "Fraction of events in channel " << plot_active_channels[i] << " above threshold: "
            << 100. * static_cast<float>(counter_abovethr[i]) / static_cast<float>(nevents) << "%\n";
    }
 
    cout << "Fraction of events w/ at least 2 channels above threshold: "
        << 100. * static_cast<float>(occurences2ch) / static_cast<float>(nevents) << "%\n"
        << "For a total of " << nevents << " events" << endl;
}
 
void ReadRun::SkipEventsPerChannel(vector<float> thresholds, float rangestart, float rangeend, bool verbose) { // merge with IntegralFilter()?
 
    if (thresholds.empty()) cout << "\nError: thresholds is empty";
    while (thresholds.size() <= active_channels.size()) thresholds.push_back(thresholds[0]);
 
    cout << "\n Removing events with individual amplitude threshold per channel:" << endl;
    int counter = 0;
    int n_thrshld = static_cast<int>(thresholds.size());
    int currchannel = 0;
 
    for (int j = 0; j < nwf; j++) {
        if (!skip_event[GetCurrentEvent(j)]) {
            currchannel = j - nchannels * GetCurrentEvent(j);
            if (currchannel < n_thrshld) {
                auto his = (TH1F*)(Getwf(j))->Clone(); // use Clone() to not change ranges of original histogram
                if (rangestart != 0 && rangeend != 0) his->GetXaxis()->SetRange(his->GetXaxis()->FindBin(rangestart), his->GetXaxis()->FindBin(rangeend));
 
                if (thresholds[currchannel] != 0. && !skip_event[GetCurrentEvent(j)] && ((thresholds[currchannel] > 0 && his->GetMaximum() > thresholds[currchannel]) || (thresholds[currchannel] < 0 && his->GetMinimum() < thresholds[currchannel]))) {
 
                    int currevent = eventnr_storage[GetCurrentEvent(j)];
                    if (verbose) cout << "\nevent:\t" << currevent << "\tchannel:\t" << active_channels[currchannel] << "\tthreshold\t" << thresholds[currchannel];
                    skip_event[GetCurrentEvent(j)] = true;
                    counter++;
                }
                delete his;
            }
        }
        Helpers::PrintProgressBar(j, nwf);
    }
 
    cout << "\t" << counter << " events will be cut out of " << nevents << endl;
}
 
void ReadRun::IntegralFilter(vector<float> thresholds, vector<bool> highlow, float windowlow, float windowhi, float start, float end, bool use_AND_condition, bool verbose) {
 
    if (thresholds.empty() || highlow.empty()) cout << "\nError: thresholds or highlow are empty";
    while (thresholds.size() <= active_channels.size()) { thresholds.push_back(thresholds[0]); }
    while (highlow.size() <= active_channels.size()) { highlow.push_back(highlow[0]); }
 
    cout << "\n\nRemoving events with individual integral threshold per channel:" << endl;
    int counter = 0;
    float integral = 0;
    int currevent_counter = 0;
    int currchannel = 0;
    int currevent = 0;
 
    for (int j = 0; j < nwf; j++) {
        currevent_counter = GetCurrentEvent(j);
 
        if (use_AND_condition || !skip_event[currevent_counter]) {
            currchannel = GetCurrentChannel(j);
 
            if (currchannel < static_cast<int>(thresholds.size())) {
                auto his = (TH1F*)(Getwf(j))->Clone(); // use Clone() to not change ranges of original histogram
                integral = GetPeakIntegral(his, windowlow, windowhi, start, end, currchannel);
 
                if (thresholds[currchannel] != 0 && !skip_event[currevent_counter] &&
                    ((highlow[currchannel] && integral > thresholds[currchannel]) || (!highlow[currchannel] && integral < thresholds[currchannel]))) {
                    currevent = eventnr_storage[currevent_counter];
                    if (verbose) cout << "\nevent:\t" << currevent << "\tchannel:\t" << active_channels[currchannel] << "\tthreshold\t" << thresholds[currchannel] << "\tintegral:\t" << integral;
                    skip_event[currevent_counter] = true;
                    while (floor((j + 1) / nchannels) == currevent_counter) j++;
                    counter++;
                }
                else if (use_AND_condition && thresholds[currchannel] != 0 && skip_event[currevent_counter] && !highlow[currchannel]) {
                    currevent = eventnr_storage[currevent_counter];
                    skip_event[currevent_counter] = false;
                    if (verbose) cout << "\nevent:\t" << currevent << "\tchannel:\t" << active_channels[currchannel] << "\thas been flagged good by integral:\t" << integral;
                }
                delete his;
            }
            Helpers::PrintProgressBar(j, nwf);
        }
    }
    cout << "\t" << counter << " events will be cut out of " << nevents << endl;
}
 
void ReadRun::PrintSkippedEvents() {
    int counter = 0;
    for (int j = 0; j < static_cast<int>(skip_event.size()); j++) {
        if (skip_event[j]) {
            int currevent = eventnr_storage[j];
            cout << "\nevent:\t" << currevent;
            counter++;
        }
    }
    cout << "\n\ntotal number of skipped events:\t" << counter << "\tout of:\t" << nevents << endl;
}
 
void ReadRun::UnskipAll() {
    for (int j = 0; j < static_cast<int>(skip_event.size()); j++) skip_event[j] = false;
    cout << "\n\nAll event cuts were removed" << endl;
}
 
int ReadRun::Nevents_good() {
    int nevents_good = 0;
    for (int i = 0; i < nevents; i++) if (!skip_event[i]) nevents_good++;
    return nevents_good;
}
 
// functions for charge spectrum
 
int* ReadRun::GetIntWindow(TH1F* his, float windowlow, float windowhi, float start, float end, int channel) {
 
    int istart, iend;
    int* foundindices = new int[3]();
 
    if (start < 0 || end < 0) {                         // fixed integration window relative to maximum of sum spectrum for each channel
        foundindices[1] = his->GetXaxis()->FindBin(his->GetXaxis()->GetBinCenter(maxSumBin[channel]) - windowlow);
        foundindices[2] = his->GetXaxis()->FindBin(his->GetXaxis()->GetBinCenter(maxSumBin[channel]) + windowhi);
    }
    else if (windowlow == start && windowhi == end) {   // fixed integration window for all channels
        foundindices[1] = his->GetXaxis()->FindBin(windowlow);
        foundindices[2] = his->GetXaxis()->FindBin(windowhi);
    }
    else {                                              // fixed integration window relative to maximum of each individual waveform
        istart = his->GetXaxis()->FindBin(start);
        iend = his->GetXaxis()->FindBin(end);
        foundindices[0] = istart;
 
        if (istart < 1 || iend > his->GetNbinsX()) {
            cout << "\nError: start or end out of range" << endl;
            return 0;
        }
 
        float max = -9.e99;
        float val = 0;
        for (int i = istart; i < iend; i++) {
            val = his->GetBinContent(i);
            if (val > max) {
                max = val;
                foundindices[0] = i;
            }
        }
 
        foundindices[1] = his->GetXaxis()->FindBin(his->GetXaxis()->GetBinCenter(foundindices[0]) - windowlow);
        foundindices[2] = his->GetXaxis()->FindBin(his->GetXaxis()->GetBinCenter(foundindices[0]) + windowhi);
    }
    return foundindices;
}
 
float ReadRun::GetPeakIntegral(TH1F* his, float windowlow, float windowhi, float start, float end, int channel_index) {
    int* windowind = GetIntWindow(his, windowlow, windowhi, start, end, channel_index); // find integration window
    string integral_option(""); // For amplitude -> unit[mV].
    if (windowind[1] != windowind[2]) integral_option = "width"; // 'width' (bin width) for integral -> unit[mV x ns].
    float integral = his->Integral(windowind[1], windowind[2], integral_option.c_str());
    delete[] windowind;
    return integral;
}
 
void ReadRun::PrintChargeSpectrumWF(float windowlow, float windowhi, float start, float end, int eventnr, float ymin, float ymax, float xmin, float xmax) {
 
    gStyle->SetOptStat(0);
 
    TString name(Form("waveforms_event__%05d", eventnr));
    auto intwinc = new TCanvas(name.Data(), name.Data(), 600, 400);
    Helpers::SplitCanvas(intwinc, active_channels, plot_active_channels);
    int event_index = GetEventIndex(eventnr);
 
    int current_canvas = 0;
    for (int i = 0; i < nchannels; i++) {
        if (PlotChannel(i)) {
            intwinc->cd(++current_canvas);
 
            auto his = Getwf(i, event_index);
            int* windowind = GetIntWindow(his, windowlow, windowhi, start, end, i);
            // create lines to indicate the integration window
            TLine* low = new TLine(his->GetXaxis()->GetBinCenter(windowind[1]), -5, his->GetXaxis()->GetBinCenter(windowind[1]), 10);
            low->SetLineColor(2);
            TLine* hi = new TLine(his->GetXaxis()->GetBinCenter(windowind[2]), -2, his->GetXaxis()->GetBinCenter(windowind[2]), 3);
            hi->SetLineColor(2);
            TLine* zero = new TLine(0, 0, 320, 0); // draw line at x=0 to check if baseline correction worked
            zero->SetLineColor(1);
            delete[] windowind;
 
            // drawing and formatting
            gPad->SetTopMargin(.01);
            int last_canvas = nchannels;
            if (!plot_active_channels.empty()) last_canvas = static_cast<int>(plot_active_channels.size());
            if (current_canvas == 1 && last_canvas < 4) gPad->SetLeftMargin(.15);
            if (current_canvas % 4 == 0 || current_canvas == last_canvas) gPad->SetRightMargin(.01);
            his->Draw("HIST");
            his->SetStats(0);
 
            low->Draw("same");
            hi->Draw("same");
            zero->Draw("same");
 
            // draw baseline parameters
            if (static_cast<int>(baseline_correction_result.size()) > event_index * nchannels + i) {
                TLine* baselinel = new TLine(baseline_correction_result[event_index * nchannels + i][2], -1, baseline_correction_result[event_index * nchannels + i][2], 1);
                baselinel->SetLineColor(6);
                baselinel->SetLineWidth(2);
                TLine* baselineh = new TLine(baseline_correction_result[event_index * nchannels + i][3], -1, baseline_correction_result[event_index * nchannels + i][3], 1);
                baselineh->SetLineColor(6);
                baselineh->SetLineWidth(2);
                TLine* baseline = new TLine(baseline_correction_result[event_index * nchannels + i][2], 0, baseline_correction_result[event_index * nchannels + i][3], 0);
                baseline->SetLineColor(6);
                TLine* correction_value = new TLine(baseline_correction_result[event_index * nchannels + i][2], baseline_correction_result[event_index * nchannels + i][0], baseline_correction_result[event_index * nchannels + i][3], baseline_correction_result[event_index * nchannels + i][0]);
                correction_value->SetLineColor(2);
 
                baselinel->Draw("same");
                baselineh->Draw("same");
                baseline->Draw("same");
                correction_value->Draw("same");
            }
 
            if (static_cast<int>(timing_results.size()) > event_index * nchannels + i) {
                TLine* timing = new TLine(timing_results[event_index * nchannels + i][1], -10, timing_results[event_index * nchannels + i][1], 100);
                timing->SetLineColor(9);
                timing->SetLineWidth(2);
                timing->Draw("same");
            }
 
            if (ymin != 0. && ymax != 0.) his->GetYaxis()->SetRangeUser(ymin, ymax); // fix y range for better comparison 
            if (xmin != 0. && xmax != 0.) his->GetXaxis()->SetRangeUser(xmin, xmax);
        }
    }
    intwinc->Update();
 
    root_out->WriteObject(intwinc, name.Data());
}
 
float* ReadRun::ChargeList(int channel_index, float windowlow, float windowhi, float start, float end, bool negative_vals) {
    float* charge_list = new float[nevents]();
    for (int j = 0; j < nevents; j++) {
        auto his = Getwf(j * nchannels + channel_index);
        charge_list[j] = GetPeakIntegral(his, windowlow, windowhi, start, end, channel_index);
        if (!negative_vals && charge_list[j] < 0.) charge_list[j] = 0.;
    }
    return charge_list;
}
 
void ReadRun::SaveChargeLists(float windowlow, float windowhi, float start, float end, bool negative_vals) {
    float* event_list = new float[nevents];
    for (int i = 0; i < nevents; i++) event_list[i] = static_cast<float>(i);
 
    auto charge_list_mg = new TMultiGraph();
    if (windowlow + windowhi > 0.) charge_list_mg->SetTitle("event-wise integrals; Event number; integral [mV#timesns]");
    else charge_list_mg->SetTitle("event-wise amplitudes; Event number; amplitude [mV]");
 
    for (int i = 0; i < nchannels; i++) {
        if (PlotChannel(i)) {
            TString name(Form("charge_list_ch_%02d", active_channels[i]));
            float* charge_list = ChargeList(i, windowlow, windowhi, start, end, negative_vals);
            TGraph* charge_list_graph = new TGraph(nevents, event_list, charge_list);
            charge_list_graph->SetLineWidth(0);
            charge_list_graph->SetMarkerStyle(2);
            charge_list_graph->SetMarkerColor(Helpers::rcolor(i));
            charge_list_graph->SetTitle(name.Data());
 
            //remove skipped events
            for (int j = 0; j < nevents; j++) {
                if (skip_event[j]) charge_list_graph->RemovePoint(j);
            }
 
            charge_list_mg->Add(charge_list_graph);
            root_out->WriteObject(charge_list_graph, name.Data());
            delete[] charge_list;
        }
    }
    root_out->WriteObject(charge_list_mg, "all_charge_lists");
    delete[] event_list;
}
 
void ReadRun::ChargeCorrelation(float windowlow, float windowhi, float start, float end, float rangestart, float rangeend, int nbins, int channel1, int channel2, bool ignore_skipped_events) {
    gStyle->SetOptStat(1111);
    stringstream name;
    name << "charge_correlation_ch" << channel1 << "_ch" << channel2;
    stringstream title;
    if (windowlow + windowhi > 0.) title << ";integral ch" << channel1 << " in mV#timesns;integral ch" << channel2 << " in mV#timesns;Entries";
    else title << ";amplitude ch" << channel1 << " in mV;amplitude ch" << channel2 << " in mV;Entries";
 
    auto charge_corr_canvas = new TCanvas(name.str().c_str(), "canvas", 600, 600);
    charge_corr_canvas->SetRightMargin(0.15);
 
    float* charge1 = ChargeList(GetChannelIndex(channel1), windowlow, windowhi, start, end);
    float* charge2 = ChargeList(GetChannelIndex(channel2), windowlow, windowhi, start, end);
 
    auto charge_corr = new TH2F(name.str().c_str(), title.str().c_str(), nbins, rangestart, rangeend, nbins, rangestart, rangeend);
    for (int i = 0; i < nevents; i++) {
        if (!ignore_skipped_events || !skip_event[i]) charge_corr->Fill(charge1[i], charge2[i]);
    }
 
    charge_corr->Draw("colz");
    root_out->WriteObject(charge_corr, name.str().c_str());
 
    charge_corr_canvas->Update();
    charge_corr_canvas->SetGrid();
    // move stat box out of the way (causing problems since May 23?)
    //TPaveStats* stat_box = (TPaveStats*)charge_corr->FindObject("stats"); 
    //stat_box->SetX1NDC(0.6); 
    //stat_box->SetX2NDC(0.85);
    charge_corr->SetStats(0);
    charge_corr_canvas->Modified();
    name << "_c";
    root_out->WriteObject(charge_corr_canvas, name.str().c_str());
    delete[] charge1;
    delete[] charge2;
}
 
TH1F* ReadRun::ChargeSpectrum(int channel_index, float windowlow, float windowhi, float start, float end, float rangestart, float rangeend, int nbins) {
 
    TString name(Form("channel__%02d", active_channels[channel_index]));
    TH1F* h1 = new TH1F(name.Data(), name.Data(), nbins, rangestart, rangeend);
 
    float pedestal = 0;
    float gain = 1;
    if (channel_index < static_cast<int>(PrintChargeSpectrum_cal.size())) {
        if (PrintChargeSpectrum_cal[channel_index][0] != 0) gain = PrintChargeSpectrum_cal[channel_index][0];
        if (PrintChargeSpectrum_cal[channel_index][1] != 1) pedestal = PrintChargeSpectrum_cal[channel_index][1];
    }
 
    for (int j = 0; j < nevents; j++) {
        if (!skip_event[j]) {
            auto his = Getwf(channel_index, j);
            h1->Fill((GetPeakIntegral(his, windowlow, windowhi, start, end, channel_index) - pedestal) / gain); 
        }
    }
    return h1;
}
 
void ReadRun::PrintChargeSpectrum(float windowlow, float windowhi, float start, float end, float rangestart, float rangeend, int nbins, float fitrangestart, float fitrangeend, int max_channel_nr_to_fit, int which_fitf, bool use_log_y) {
    // print ReadRun::ChargeSpectrum for all channels optimized for SiPM signals
    checkData();
    
    gStyle->SetOptStat("ne");
    gStyle->SetOptFit(1111);
 
    if (fitrangestart == 0.)    fitrangestart = rangestart;
    if (fitrangeend == 0.)      fitrangeend = rangeend;
 
    string ctitle("\"charge\" spectra" + to_string(PrintChargeSpectrum_cnt++));
    auto chargec = new TCanvas(ctitle.c_str(), ctitle.c_str(), 600, 400);
    Helpers::SplitCanvas(chargec, active_channels, plot_active_channels);
 
    cout << "\n\nThere is data recorded in " << active_channels.size() << " channels \n\n\n";
    int current_canvas = 0;
 
    float default_rangestart = -2000;
    float default_rangeend = 30000;
    if (default_rangestart > rangestart) default_rangestart = rangestart;
    if (default_rangeend < rangeend) default_rangeend = rangeend;
 
    int default_nbins = static_cast<int>((default_rangeend - default_rangestart) * nbins / (rangeend - rangestart));
 
    for (int i = 0; i < nchannels; i++) {
        if (PlotChannel(i)) {
            chargec->cd(++current_canvas);
            if (use_log_y) gPad->SetLogy();
 
            auto his = ChargeSpectrum(i, windowlow, windowhi, start, end, default_rangestart, default_rangeend, default_nbins);
            his->GetYaxis()->SetTitle("#Entries");
            if (windowlow + windowhi > 0.) his->GetXaxis()->SetTitle("Integral in mV#timesns");
            else his->GetXaxis()->SetTitle("Amplitude in mV");
            
            if (i < static_cast<int>(PrintChargeSpectrum_cal.size()) && PrintChargeSpectrum_cal[i][0] != 1) {
                cout << "Charge spectrum for channel index " << i << " will be normalized using a gain of "
                    << PrintChargeSpectrum_cal[i][0] << " and a pedestal value of " << PrintChargeSpectrum_cal[i][1] << endl;
                his->GetXaxis()->SetTitle("Number of photoelectrons");
            }
 
            //store the mean integral of each channel --> used for correction factors of phi_ew analysis
            mean_integral.push_back(his->GetMean());
 
            //fitting
            if (i < max_channel_nr_to_fit) {
                if (which_fitf == 0) {}
                else if (which_fitf == 1) { // landau gauss convolution for large number of photons
                    Fitf_langaus fitf;
                    int n_par = 4;
                    TF1* f = new TF1("fitf_langaus", fitf, fitrangestart, fitrangeend, n_par); f->SetLineColor(3); f->SetNpx(1000);
 
                    f->SetParName(0, "Width");              f->SetParameter(0, 35);
                    f->SetParName(1, "MPV");                f->SetParameter(1, 1000);
                    f->SetParName(2, "Area");               f->SetParameter(2, 10000);
                    f->SetParName(3, "#sigma_{Gauss}");     f->SetParameter(3, 100);
 
                    if (!PrintChargeSpectrum_pars.empty()) for (int j = 0; j < min(4, static_cast<int>(PrintChargeSpectrum_pars.size())); j++) f->SetParameter(j, PrintChargeSpectrum_pars[j]);
 
                    if (i < max_channel_nr_to_fit) {
                        cout << "\n\n---------------------- Fit for channel " << active_channels[i] << " ----------------------\n";
                        TFitResultPtr fresults = his->Fit(f, "LRS");
                        fit_results.push_back(fresults);
                    }
                }
                else if (which_fitf == 2) { // if pedestal is biased because of peak finder algorithm
                    Fitf_biased fitf;
                    int n_par = 9;
                    TF1* f = new TF1("fitf_biased", fitf, fitrangestart, fitrangeend, n_par); f->SetLineColor(3); f->SetNpx(1000);
 
                    f->SetParName(0, "N_{0}");              f->SetParameter(0, his->Integral());
                    f->SetParName(1, "#mu");                f->SetParameter(1, 2.);
                    f->SetParName(2, "#lambda");            f->SetParameter(2, .04);
                    f->SetParName(3, "#sigma_{0}");         f->SetParameter(3, 2.1);
                    f->SetParName(4, "#sigma_{1}");         f->SetParameter(4, 3.4); //f->SetParLimits(4, 1.e-9, 1.e3);
                    f->SetParName(5, "Gain");               f->SetParameter(5, 30.); //f->FixParameter(5, 10.);
                    f->SetParName(6, "Pedestal");           f->SetParameter(6, 2.);
                    f->SetParName(7, "norm_{0}");           f->SetParameter(7, 0.7);
                    f->SetParName(8, "x_{0}");              f->SetParameter(8, 5.);
 
                    if (!PrintChargeSpectrum_pars.empty()) for (int j = 0; j < min(n_par, static_cast<int>(PrintChargeSpectrum_pars.size())); j++) f->SetParameter(j, PrintChargeSpectrum_pars[j]);
 
                    if (i < max_channel_nr_to_fit) {
                        cout << "\n\n---------------------- Fit for channel " << active_channels[i] << " ----------------------\n";
                        TFitResultPtr fresults = his->Fit(f, "LRS");
                        fit_results.push_back(fresults);
                    }
 
                    // get number of excess events in the pedestal in the fit region. To get the absolute number of excess events the full pedestal needs to be inside of the fit range (fitrangestart, fitrangeend)
                    //double excessEventsInPedestal = f->Integral(fitrangestart, fitrangeend)/.3125;
                    //f->SetParameter(7, 1.);
                    //excessEventsInPedestal -= f->Integral(fitrangestart, fitrangeend)/.3125;
                    //cout << "\nNumber of excess events in the pedestal within the fit range:\t" << excessEventsInPedestal << "\n\n";
                }
                else if (which_fitf == 3) { // SiPM fit function with exponential delayed after pulsing
                    Fitf_full fitf;
                    int n_par = 9;
                    TF1* f = new TF1("fitf", fitf, fitrangestart, fitrangeend, n_par); f->SetLineColor(3); f->SetNpx(1000);
 
                    f->SetParName(0, "N_{0}");              f->SetParameter(0, his->Integral());
                    f->SetParName(1, "#mu");                f->SetParameter(1, 2.);
                    f->SetParName(2, "#lambda");            f->SetParameter(2, .04);
                    f->SetParName(3, "#sigma_{0}");         f->SetParameter(3, 2.1);
                    f->SetParName(4, "#sigma_{1}");         f->SetParameter(4, 3.4);
                    f->SetParName(5, "Gain");               f->SetParameter(5, 30.); //f->FixParameter(5, 40.);
                    f->SetParName(6, "Pedestal");           f->SetParameter(6, 2.);
                    f->SetParName(7, "#alpha");             f->SetParameter(7, .1); //f->FixParameter(7, .2);
                    f->SetParName(8, "#beta");              f->SetParameter(8, 80.); //f->FixParameter(8, 80);
 
                    if (!PrintChargeSpectrum_pars.empty()) for (int j = 0; j < min(n_par, static_cast<int>(PrintChargeSpectrum_pars.size())); j++) f->SetParameter(j, PrintChargeSpectrum_pars[j]);
 
                    if (i < max_channel_nr_to_fit) {
                        cout << "\n\n---------------------- Fit for channel " << active_channels[i] << " ----------------------\n";
                        TFitResultPtr fresults = his->Fit(f, "LRS");
                        fit_results.push_back(fresults);
                    }
                }
                else if (which_fitf == 4) { // ideal PMT fit function
                    Fitf_PMT_ideal fitf;
                    int n_par = 4;
                    TF1* f = new TF1("fitf", fitf, fitrangestart, fitrangeend, n_par); f->SetLineColor(3); f->SetNpx(1000);
 
                    f->SetParName(0, "N_{0}");              f->SetParameter(0, his->Integral());
                    f->SetParName(1, "#mu");                f->SetParameter(1, 1.);
                    f->SetParName(2, "#sigma");             f->SetParameter(2, 5.);
                    f->SetParName(3, "gain");               f->SetParameter(3, 10.);
 
                    if (!PrintChargeSpectrum_pars.empty()) for (int j = 0; j < min(n_par, static_cast<int>(PrintChargeSpectrum_pars.size())); j++) f->SetParameter(j, PrintChargeSpectrum_pars[j]);
 
                    if (i < max_channel_nr_to_fit) {
                        cout << "\n\n---------------------- Fit for channel " << active_channels[i] << " ----------------------\n";
                        TFitResultPtr fresults = his->Fit(f, "LRS");
                        fit_results.push_back(fresults);
                    }
                }
                else if (which_fitf == 5) { // PMT fit function
                    Fitf_PMT fitf;
                    int n_par = 8;
                    TF1* f = new TF1("fitf", fitf, fitrangestart, fitrangeend, n_par); f->SetLineColor(3); f->SetNpx(1000);
 
                    f->SetParName(0, "N_{0}");              f->SetParameter(0, his->Integral());
                    f->SetParName(1, "w");                  f->SetParameter(1, .05);    f->SetParLimits(1, 1.e-99, 4.e-1); //probability for type II BG
                    f->SetParName(2, "#alpha");             f->SetParameter(2, .05);    f->SetParLimits(2, 1.e-99, 5.e-2); //coefficient of exponential decrease of typ II BG
                    f->SetParName(3, "#sigma_{0}");         f->SetParameter(3, 5.);     f->SetParLimits(3, 1.e-9, 1.e3);
                    f->SetParName(4, "Q_{0}");              f->SetParameter(4, 0.);
                    f->SetParName(5, "#mu");                f->SetParameter(5, 1.);
                    f->SetParName(6, "#sigma_{1}");         f->SetParameter(6, 5.);     f->SetParLimits(6, 1.e-9, 1.e3);
                    f->SetParName(7, "Q_{1}");              f->SetParameter(7, 10.);
 
                    if (!PrintChargeSpectrum_pars.empty()) for (int j = 0; j < min(n_par, static_cast<int>(PrintChargeSpectrum_pars.size())); j++) f->SetParameter(j, PrintChargeSpectrum_pars[j]);
 
                    if (i < max_channel_nr_to_fit) {
                        cout << "\n\n---------------------- Fit for channel " << active_channels[i] << " ----------------------\n";
                        TFitResultPtr fresults = his->Fit(f, "LRS");
                        fit_results.push_back(fresults);
                    }
                }
                else if (which_fitf == 6) { // PMT fit function with biased pedestal
                    Fitf_PMT_pedestal fitf;
                    int n_par = 9;
                    TF1* f = new TF1("fitf", fitf, fitrangestart, fitrangeend, n_par); f->SetLineColor(3); f->SetNpx(1000);
 
                    f->SetParName(0, "A");                  f->SetParameter(0, his->Integral());
                    f->SetParName(1, "w");                  f->SetParameter(1, .05);    f->SetParLimits(1, 1.e-9, 4.e-1);   //probability for type II BG
                    f->SetParName(2, "#alpha");             f->SetParameter(2, .05);    f->SetParLimits(2, 1.e-9, 5.e-2);   //coefficient of exponential decrease of typ II BG
                    f->SetParName(3, "#sigma_{0}");         f->SetParameter(3, 5.);     f->SetParLimits(3, 1.e-9, 1.e3);
                    f->SetParName(4, "Q_{0}");              f->SetParameter(4, 0.);     f->SetParLimits(4, -1.e2, 1.e2);
                    f->SetParName(5, "#mu");                f->SetParameter(5, 1.);     f->SetParLimits(5, 1.e-9, 1.e2);
                    f->SetParName(6, "#sigma_{1}");         f->SetParameter(6, 5.);     f->SetParLimits(6, 1.e-9, 1.e3);
                    f->SetParName(7, "Q_{1}");              f->SetParameter(7, 10.);    f->SetParLimits(7, 1.e-9, 1.e9);
                    f->SetParName(8, "A_{0}");              f->SetParameter(8, 1.);     f->SetParLimits(8, 1.e-9, 1.e1);
 
                    if (!PrintChargeSpectrum_pars.empty()) for (int j = 0; j < min(n_par, static_cast<int>(PrintChargeSpectrum_pars.size())); j++) f->SetParameter(j, PrintChargeSpectrum_pars[j]);
 
                    if (i < max_channel_nr_to_fit) {
                        cout << "\n\n---------------------- Fit for channel " << active_channels[i] << " ----------------------\n";
                        TFitResultPtr fresults = his->Fit(f, "LRS");
                        fit_results.push_back(fresults);
                    }
                }
                else if (which_fitf == 7) { // default SiPM fit function + dark count spectrum (for lots of false triggers)
                    Fitf_plus_DC fitf;
                    int n_par = 9;
                    TF1* f = new TF1("fitf", fitf, fitrangestart, fitrangeend, n_par); f->SetLineColor(3); f->SetNpx(1000);
 
                    f->SetParName(0, "A");                  f->SetParameter(0, his->Integral());
                    f->SetParName(1, "w");                  f->SetParameter(1, .05);    f->SetParLimits(1, 1.e-9, 4.e-1);   //probability for type II BG
                    f->SetParName(2, "#alpha");             f->SetParameter(2, .05);    f->SetParLimits(2, 1.e-9, 5.e-2);   //coefficient of exponential decrease of typ II BG
                    f->SetParName(3, "#sigma_{0}");         f->SetParameter(3, 5.);     f->SetParLimits(3, 1.e-9, 1.e3);
                    f->SetParName(4, "Q_{0}");              f->SetParameter(4, 0.);     f->SetParLimits(4, -1.e2, 1.e2);
                    f->SetParName(5, "#mu");                f->SetParameter(5, 1.);     f->SetParLimits(5, 1.e-9, 1.e2);
                    f->SetParName(6, "#sigma_{1}");         f->SetParameter(6, 5.);     f->SetParLimits(6, 1.e-9, 1.e3);
                    f->SetParName(7, "#mu_darkcount");      f->SetParameter(7, .1);     f->SetParLimits(7, 1.e-9, 1.);
                    f->SetParName(8, "N_{0}_darkcount");    f->SetParameter(8, .05);    f->SetParLimits(8, 1.e-9, .3);
 
                    if (!PrintChargeSpectrum_pars.empty()) for (int j = 0; j < min(n_par, static_cast<int>(PrintChargeSpectrum_pars.size())); j++) f->SetParameter(j, PrintChargeSpectrum_pars[j]);
 
                    if (i < max_channel_nr_to_fit) {
                        cout << "\n\n---------------------- Fit for channel " << active_channels[i] << " ----------------------\n";
                        TFitResultPtr fresults = his->Fit(f, "LRS");
                        fit_results.push_back(fresults);
                    }
                }
                else { // default SiPM fit function
                    Fitf fitf;
                    int n_par = 7;
                    TF1* f = new TF1("fitf", fitf, fitrangestart, fitrangeend, n_par); f->SetLineColor(3); f->SetNpx(1000);
 
                    f->SetParName(0, "N_{0}");              f->SetParameter(0, his->Integral());
                    f->SetParName(1, "#mu");                f->SetParameter(1, 2.);
                    f->SetParName(2, "#lambda");            f->SetParameter(2, .04);
                    f->SetParName(3, "#sigma_{0}");         f->SetParameter(3, 2.1);
                    f->SetParName(4, "#sigma_{1}");         f->SetParameter(4, 3.4);
                    f->SetParName(5, "Gain");               f->SetParameter(5, 30.); //f->FixParameter(5, 40.);
                    f->SetParName(6, "Pedestal");           f->SetParameter(6, 2.);
 
                    if (!PrintChargeSpectrum_pars.empty()) for (int j = 0; j < min(n_par, static_cast<int>(PrintChargeSpectrum_pars.size())); j++) f->SetParameter(j, PrintChargeSpectrum_pars[j]);
 
                    if (i < max_channel_nr_to_fit) {
                        cout << "\n\n---------------------- Fit for channel " << active_channels[i] << " ----------------------\n";
                        TFitResultPtr fresults = his->Fit(f, "LRS");
                        fit_results.push_back(fresults);
                    }
                }
            }
            TString name(Form("ChargeSpectrum channel_%02d_%d", active_channels[i], PrintChargeSpectrum_cnt));
            root_out->WriteObject(his, name.Data());
            his->Draw();
        }
    }
 
    Helpers::SetRangeCanvas(chargec, rangestart, rangeend);
    root_out->WriteObject(chargec, ("ChargeSpectra" + to_string(PrintChargeSpectrum_cnt)).c_str());
}
 
void ReadRun::PrintDCR(float windowlow, float windowhi, float rangestart, float rangeend, double threshold) {
 
    string unit(" mV");
    if (windowlow + windowhi > 0.) unit = " mV*ns";
 
    for (int i = 0; i < nchannels; i++) {
        if (PlotChannel(i)) {
            auto his = ChargeSpectrum(i, windowlow, windowhi, rangestart, rangeend, rangestart, rangeend, 500);
 
            stringstream lonamerate;
            lonamerate << "<0.5 pe=" << threshold << unit << " -> " << his->Integral(his->GetXaxis()->FindBin(rangestart), his->GetXaxis()->FindBin(threshold)) / his->GetEntries() / (1.e-3 * (rangeend - rangestart)) << " MHz";
            cout << "\n" << lonamerate.str().c_str() << endl;
 
            stringstream hinamerate;
            hinamerate << ">0.5 pe=" << threshold << unit << " -> " << his->Integral(his->GetXaxis()->FindBin(threshold) + 1, his->GetXaxis()->FindBin(rangeend)) / his->GetEntries() / (1.e-3 * (rangeend - rangestart)) << " MHz";
            cout << "\n" << hinamerate.str().c_str() << endl;
        }
    }
}
 
TH1F* ReadRun::TimeDist(int channel_index, float from, float to, float rangestart, float rangeend, int nbins, int which, float cf_r) {
 
    TString name(Form("timedist_ch%02d", active_channels[channel_index]));
    auto h1 = new TH1F(name.Data(), name.Data(), nbins, rangestart, rangeend);
 
    for (int j = 0; j < nevents; j++) {
        if (!skip_event[j]) {
            auto his = (TH1F*)(Getwf(j * nchannels + channel_index));
 
            int from_n = his->GetXaxis()->FindBin(from);
 
            int max_n = GetIntWindow(his, 0, 0, from, to)[0];
            float max = his->GetBinContent(max_n);
 
            if (which == 0) { // time of maximum 
                h1->Fill(max);
            }
            else if (which == 1) { // time of 50% CFD
                do {
                    max_n--;
                } while (his->GetBinContent(max_n) >= cf_r * max && max_n > from_n);
                max_n++;
 
                h1->Fill(LinearInterpolation(cf_r * max, his->GetXaxis()->GetBinCenter(max_n - 1), his->GetXaxis()->GetBinCenter(max_n), his->GetBinContent(max_n - 1), his->GetBinContent(max_n)).first);
            }
            else { // 10%-90% rise time
                // search backwards from maximum
                int n10 = -1;
                int n90 = -1;
                do {
                    max_n--;
                    if (n10 == -1 && his->GetBinContent(max_n) >= .1 * max && his->GetBinContent(max_n - 1) <= .1 * max) n10 = max_n;
                    if (n90 == -1 && his->GetBinContent(max_n) >= .9 * max && his->GetBinContent(max_n - 1) <= .9 * max) n90 = max_n;
                } while (his->GetBinContent(max_n) <= max && max_n > from_n);
 
                float t10 = LinearInterpolation(.1 * max, his->GetXaxis()->GetBinCenter(n10 - 1), his->GetXaxis()->GetBinCenter(n10), his->GetBinContent(n10 - 1), his->GetBinContent(n10)).first;
                float t90 = LinearInterpolation(.9 * max, his->GetXaxis()->GetBinCenter(n90 - 1), his->GetXaxis()->GetBinCenter(n90), his->GetBinContent(n90 - 1), his->GetBinContent(n90)).first;
 
                h1->Fill(t90 - t10);
            }
        }
    }
    if (which == 1) h1->Fit("gaus", "WWM", "same");
    return h1;
}
 
void ReadRun::PrintTimeDist(float from, float to, float rangestart, float rangeend, int nbins, int which, float cf_r) {
    gStyle->SetOptStat(1111); // 11 is title + entries
 
    auto time_dist_c = new TCanvas("timing of maximum", "timing of maximum", 600, 400);
    Helpers::SplitCanvas(time_dist_c, active_channels, plot_active_channels);
 
    int current_canvas = 0;
 
    for (int i = 0; i < nchannels; i++) {
        if (PlotChannel(i)) {
            time_dist_c->cd(++current_canvas);
 
            auto his = TimeDist(i, from, to, rangestart, rangeend, nbins, which, cf_r);
            his->GetYaxis()->SetTitle("#Entries");
            his->GetXaxis()->SetTitle("time [ns]");
            his->Draw();
            stringstream name; name << "t_{max} for " << from << "<t<" << to << " ns";
            his->SetTitle(name.str().c_str());
 
            TString name_save(Form("TimeDist channel_%02d", active_channels[i]));
            root_out->WriteObject(his, name_save.Data());
        }
    }
 
    time_dist_c->Update();
    root_out->WriteObject(time_dist_c, "TimeDist");
}
 
TGraph2D* ReadRun::MaxDist(int channel_index, float from, float to) {
    // find maximum amplitude for a given channel in time window [from, to] and return 3d histogram with the number of bins nbinsy,z
 
    TString name(Form("maxdist_ch%02d", active_channels[channel_index]));
    TGraph2D* g3d = new TGraph2D((binNumber + 2) * nevents);
    g3d->SetTitle("waveforms; t [ns]; max. amplitude [mv]; amplitude [mV]");
 
    for (int j = 0; j < nevents; j++) {
        if (!skip_event[j]) {
            auto his = (TH1F*)(Getwf(j * nchannels + channel_index))->Clone();
            if (from >= 0 && to > 0) his->GetXaxis()->SetRange(his->GetXaxis()->FindBin(from), his->GetXaxis()->FindBin(to));
            double max = his->GetMaximum();
            for (int i = 0; i < binNumber; i++) g3d->SetPoint(j * binNumber + i, his->GetXaxis()->GetBinCenter(i), max, his->GetBinContent(i));
            delete his;
        }
    }
    root_out->WriteObject(g3d, name.Data());
    return g3d;
}
 
void ReadRun::PrintMaxDist(float from, float to) {
 
    auto max_dist_c = new TCanvas("wf grouped by maximum", "wf grouped by maximum", 600, 400);
    Helpers::SplitCanvas(max_dist_c, active_channels, plot_active_channels);
 
    int current_canvas = 0;
 
    for (int i = 0; i < nchannels; i++) {
        if (PlotChannel(i)) {
            max_dist_c->cd(++current_canvas);
            auto g3d = MaxDist(i, from, to);
            g3d->Draw();
        }
    }
    max_dist_c->Update();
    root_out->WriteObject(max_dist_c, "MaxDist");
}
 
TH1F* ReadRun::His_GetTimingCFD(int channel_index, float rangestart, float rangeend, int nbins) {
 
    if (nbins == -999) nbins = static_cast<int>((rangeend - rangestart) / SP);
 
    TString name(Form("GetTimingCFD_ch%02d", active_channels[channel_index]));
    auto h1 = new TH1F(name.Data(), name.Data(), nbins, rangestart, rangeend);
    for (int j = 0; j < nevents; j++) if (!skip_event[j]) h1->Fill(timing_results[j * nchannels + channel_index][1]);
    return h1;
}
 
void ReadRun::Print_GetTimingCFD(float rangestart, float rangeend, int do_fit, int nbins, string fitoption, bool set_errors) {
 
    // call GetTimingCFD() in case it was not initialized
    if (static_cast<int>(timing_results.size()) == 0) GetTimingCFD();
 
    gStyle->SetOptStat(1111);
    gStyle->SetOptFit(111);
 
    auto timing_cfd_c = new TCanvas("timing of cfd", "timing of cfd", 600, 400);
    Helpers::SplitCanvas(timing_cfd_c, active_channels, plot_active_channels);
    int current_canvas = 0;
 
    for (int i = 0; i < nchannels; i++) {
        if (PlotChannel(i)) {
            timing_cfd_c->cd(++current_canvas);
 
            auto his = His_GetTimingCFD(i, rangestart, rangeend, nbins);
            his->GetYaxis()->SetTitle("#Entries");
            his->GetXaxis()->SetTitle("time [ns]");
 
            if (set_errors) {
                int end = his->GetNbinsX();
                for (int i = 1; i <= end; i++) {
                    if (his->GetBinContent(i) < 2) his->SetBinError(i, 1);
                    else his->SetBinError(i, sqrt(his->GetBinContent(i)));
                }
            }
 
            his->Draw();
 
            if (do_fit == 1) {
                TFitResultPtr fresults = his->Fit("gaus", fitoption.c_str(), "same");
                timing_fit_results.push_back(fresults);
            }
 
            TString name_save(Form("Timing_cfd_channel_%02d", active_channels[i]));
            root_out->WriteObject(his, name_save.Data());
        }
    }
 
    timing_cfd_c->Update();
    root_out->WriteObject(timing_cfd_c, "TimingCFD");
}
 
TH1F* ReadRun::His_GetTimingCFD_diff(vector<int> channels1, vector<int> channels2, float rangestart, float rangeend, int nbins) {
 
    if (nbins == -999) nbins = static_cast<int>((rangeend - rangestart) / SP);
 
    stringstream name;
    name << "GetTimingCFD_diff <";
 
    // find channel indices and assemble title
    int counter = 0;
    for (int& entry : channels2) {
        if (counter > 0) name << "&";
        auto chin2 = find(active_channels.begin(), active_channels.end(), entry);
        if (chin2 != active_channels.end()) {
            name << "ch" << entry;
            entry = chin2 - active_channels.begin();
        }
        else cout << "\n\n ERROR: channels2 = " << entry << " does not exist in data. Check parameters for Print_GetTimingCFD_diff()\n\n";
        counter++;
    }
    name << ">-<";
    counter = 0;
 
    for (int& entry : channels1) {
        if (counter > 0) name << "&";
        auto chin1 = find(active_channels.begin(), active_channels.end(), entry);
        if (chin1 != active_channels.end()) {
            name << "ch" << entry;
            entry = chin1 - active_channels.begin();
        }
        else cout << "\n\n ERROR: channels1 = " << entry << " does not exist in data. Check parameters for Print_GetTimingCFD_diff()\n\n";
        counter++;
    }
    name << ">";
 
    // fill histogram
    auto h1 = new TH1F(name.str().c_str(), name.str().c_str(), nbins, rangestart, rangeend);
    for (int j = 0; j < nevents; j++) {
        if (!skip_event[j]) {
            float mean1 = 0., mean2 = 0., cnt1 = 0., cnt2 = 0.;
            for (int i : channels1) {
                mean1 += timing_results[j * nchannels + i][1];
                cnt1 += 1.;
            }
            for (int i : channels2) {
                mean2 += timing_results[j * nchannels + i][1];
                cnt2 += 1.;
            }
            h1->Fill(mean2 / cnt2 - mean1 / cnt1);
        }
    }
 
    return h1;
}
 
void ReadRun::Print_GetTimingCFD_diff(vector<int> channels1, vector<int> channels2, float rangestart, float rangeend, int do_fit, int nbins, float fitrangestart, float fitrangeend, string fitoption, bool set_errors) {
 
    // call GetTimingCFD() in case it was not initialized
    if (static_cast<int>(timing_results.size()) == 0) GetTimingCFD();
 
    if (fitrangestart == -999) {
        fitrangestart = rangestart;
        fitrangeend = rangeend;
    }
 
    //gStyle->SetOptStat(1111);
    gStyle->SetOptFit(111);
 
    auto timing_cfd_d_c = new TCanvas("timing of cfd diff", "timing of cfd diff", 600, 400);
 
    auto his = His_GetTimingCFD_diff(channels1, channels2, rangestart, rangeend, nbins);
    his->GetYaxis()->SetTitle("#Entries");
    his->GetXaxis()->SetTitle("time [ns]");
 
    if (set_errors) {
        int end = his->GetNbinsX();
        for (int i = 1; i <= end; i++) {
            if (his->GetBinContent(i) < 2) his->SetBinError(i, 1);
            else his->SetBinError(i, sqrt(his->GetBinContent(i)));
        }
    }
 
    his->Draw();
 
    double skewness = his->GetSkewness();
    
    if (do_fit == 1 || (do_fit == 2 && abs(skewness) < .15)) {
        // gauss (default)
        TFitResultPtr fresults = his->Fit("gaus", fitoption.c_str(), "same", fitrangestart, fitrangeend);
        timing_fit_results.push_back(fresults);
        if (do_fit == 2) cout << "\nWARNING: Print_GetTimingCFD_diff\nFITTING GAUSS INSTEAD OF GAUSS x EXP CONVOLUTION BC SYMMETRY" << endl;
    }
    else if (do_fit == 2) {
        // gauss x exp convolution (effective delay from random light path and/or self-absorption and reemission)
        Fitf_exp_gauss fitf_exp_gauss;
        auto expgconv = new TF1("exp x gauss convolution", fitf_exp_gauss, fitrangestart, fitrangeend, 4);
        expgconv->SetNpx(5000);
 
        // this parameter describes the sigma from different light paths 
        // and/or the effective decay time constant for self-absorption and reemission
        expgconv->SetParName(0, "#tau_{eff}");      expgconv->SetParameter(0, skewness);
        if (skewness > 0) expgconv->SetParLimits(0, .15, 5.);
        else expgconv->SetParLimits(0, -5., -.15);
        //expgconv->FixParameter(0, 1.55);
 
        expgconv->SetParName(1, "#sigma_{gaus}");       expgconv->SetParameter(1, his->GetStdDev());
        expgconv->SetParLimits(1, 1e-1, 7.);    //expgconv->FixParameter(1, .7);
 
        expgconv->SetParName(2, "t_{0}");       expgconv->SetParameter(2, his->GetMean());
        expgconv->SetParLimits(2, fitrangestart, fitrangeend);  //expgconv->FixParameter(2, 6.6);
 
        expgconv->SetParName(3, "norm");        expgconv->SetParameter(3, his->Integral("width"));
        expgconv->SetParLimits(3, 1., 1e8);     //expgconv->FixParameter(3, 105.5);
 
        TFitResultPtr fresults = his->Fit(expgconv, "SR", "same");
        timing_fit_results.push_back(fresults);
 
        // for the phi_ew-analysis: print out the time value of the maximum of the best fit --> used to determine timing cuts
        float t_of_maximum = expgconv->GetMaximumX(-5, 5);
        cout << "Maximum of the fit is at t=" << t_of_maximum << " ns and the ";
 
        double max_val = expgconv->GetMaximum();
        double fwhm_x1 = expgconv->GetX(max_val / 2, fitrangestart, fitrangeend);
        double fwhm_x2 = expgconv->GetX(max_val / 2, fwhm_x1 + 1e-3, fitrangeend);
        double fwhm = fwhm_x2 - fwhm_x1;
        auto fwhm_line = new TLine(fwhm_x1, max_val/2, fwhm_x2, max_val/2);
        fwhm_line->SetLineColor(2); fwhm_line->SetLineWidth(2);
        fwhm_line->Draw("same");
        cout << "FWHM=" << fwhm << " ns" << endl;
 
        // TLatex l;
        // l.SetTextSize(0.025);
        // l.DrawLatex(t_of_maximum/(rangeend - rangestart), 0.4, Form("FWHM = %.2f ns", fwhm));
 
        auto mean = new TLine(expgconv->GetParameter(2), 1e-2, expgconv->GetParameter(2), his->GetMaximum());
        mean->SetLineColor(1); mean->SetLineWidth(2);
        mean->Draw("same");
    }
    else if (do_fit == 3) {
        // sum of two gaussians (one as background estimate)
        auto two_gauss = new TF1("two gaussians", "gaus(0)+gaus(3)", rangestart, rangeend);
        two_gauss->SetTitle("Sum of two gauss");
        float posmax = his->GetXaxis()->GetBinCenter(his->GetMaximumBin());
        two_gauss->SetParameters(his->Integral("width"), posmax, 0.35, his->Integral("width") / 30, posmax, 2);
        two_gauss->SetParName(0, "norm_{peak}");        two_gauss->SetParName(1, "#mu_{peak}");         two_gauss->SetParName(2, "#sigma_{peak}");          two_gauss->SetParLimits(2, 1e-9, 1e2);
        two_gauss->SetParName(3, "norm_{background}");  two_gauss->SetParName(4, "#mu_{background}");   two_gauss->SetParName(5, "#sigma_{background}");    two_gauss->SetParLimits(5, 1e-9, 1e2);
        TFitResultPtr fresults = his->Fit(two_gauss, fitoption.c_str(), "same", fitrangestart, fitrangeend);
        timing_fit_results.push_back(fresults);
    }
 
    root_out->WriteObject(his, his->GetTitle());
    timing_cfd_d_c->Update();
    root_out->WriteObject(timing_cfd_d_c, "TimingCFD_diff");
}
 
 
 
 
 
 
TH1F* ReadRun::Getwf(int wf_nr) {
    checkData();
    return (TH1F*)rundata->At(wf_nr);
}
 
TH1F* ReadRun::Getwf(int channelnr, int eventnr, int color) {
    checkData();
    TH1F* his;
    his = (TH1F*)rundata->At(eventnr * nchannels + channelnr);
    his->SetLineColor(color);
    his->SetMarkerColor(color);
    return his;
}
 
double* ReadRun::getx(double shift) {
    double* xvals = new double[binNumber];
    for (int i = 0; i < binNumber; i++) {
        xvals[i] = static_cast<double>(SP) * static_cast<double>(i) + shift;
    }
    return xvals;
}
 
double* ReadRun::gety(int waveform_index) {
    auto his = Getwf(waveform_index);
    return Helpers::gety(his);
}
 
double* ReadRun::gety(int channelnr, int eventnr) {
    auto his = Getwf(channelnr, eventnr);
    return Helpers::gety(his);
}
 
int ReadRun::GetEventIndex(unsigned int eventnr) {
    auto event_pos = find(eventnr_storage.begin(), eventnr_storage.end(), eventnr);
    if (event_pos == eventnr_storage.end()) {
        cout << "WARNING: Event number " << eventnr << " for GetEventIndex() does not exist in data.\n"
            << "Please check the events in the data or set discard_original_eventnr = true before calling ReadFile()." << endl;
        
        if (static_cast<int>(eventnr) < nevents) {
            eventnr = eventnr_storage[eventnr];
            event_pos = find(eventnr_storage.begin(), eventnr_storage.end(), eventnr);
            cout << "Found event number " << eventnr << " in data and will use it." << endl;
        }
    }
    return static_cast<int>(distance(eventnr_storage.begin(), event_pos));
}
 
int ReadRun::GetChannelIndex(int channel_number) {
    int channel_index = -1;
    for (int i = 0; i < static_cast<int>(active_channels.size()); i++) {
        if (active_channels[i] == channel_number) channel_index = i;
    }
    if (channel_index == -1) {
        cout << "\n\n\tERROR: channel " << channel_number << " does not exist in data. Will continue with first channel\n\n";
        channel_index = 0;
    }
    return channel_index;
}
 
int ReadRun::GetCurrentChannel(int waveform_index) {
    return (waveform_index - nchannels * floor(waveform_index / nchannels));
}
 
int ReadRun::GetCurrentEvent(int waveform_index) {
    return floor(waveform_index / nchannels);
}
 
bool ReadRun::PlotChannel(int i) {
    if (plot_active_channels.empty() || find(plot_active_channels.begin(), plot_active_channels.end(), active_channels[i]) != plot_active_channels.end()) return true;
    else return false;
}
 
pair<float, bool> ReadRun::LinearInterpolation(float ym, float x1, float x2, float y1, float y2) {
    if (y1 == y2) return {(x1 + x2) / 2., false};
    else if (y1 > ym) {
        cout << "\nError in LinearInterpolation: Value ym=" << ym << " out of range (" << y1 << "|" << y2 << ").\n"
            << "Will return x1. Increase window for search.\n";
        return {x1, false};
    }
    else return {x1 + (ym - y1) * (x2 - x1) / (y2 - y1), true};
}