tgraphs_to_ttree.C

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <filesystem>
 
std::map<int, float> thresholds;
std::map<int, float> mean_1pe;
std::map<int, float> sigma_1pe;
void set_thresholds()
{
    thresholds[51] = 0.25;
    thresholds[52] = 0.30;
    thresholds[53] = 0.35;
    thresholds[55] = 0.45;
    thresholds[57] = 0.55;
    mean_1pe[51] = 1.72219e-01;
    mean_1pe[52] = 2.24019e-01;
    mean_1pe[53] = 2.73532e-01;
    mean_1pe[55] = 3.74145e-01;
    mean_1pe[57] = 4.60333e-01;
    sigma_1pe[51] = 4.21654e-03;
    sigma_1pe[52] = 4.22809e-03;
    sigma_1pe[53] = 5.20792e-03;
    sigma_1pe[55] = 5.19898e-03;
    sigma_1pe[57] = 6.05599e-03;
}
 
std::pair<float, float> find_first_amplitude(TGraph *gTarget, int start_point = 5000, bool maximum = true)
{
    auto maximum_points = gTarget->GetN();
    auto amplitude = 0.;
    auto amplitude_x = 0.;
    for (auto iPnt = start_point; iPnt <= maximum_points; iPnt++)
    {
        auto current_x = gTarget->GetPointX(iPnt);
        auto current_y = gTarget->GetPointY(iPnt);
        if (current_y <= -0.99)
            continue;
        if ((maximum ? (amplitude < current_y) : (amplitude > current_y)))
        {
            amplitude = current_y;
            amplitude_x = iPnt;
        }
        if ((maximum ? ((current_y - amplitude) < -0.02) : ((current_y - amplitude) > 0.02)) && (fabs(amplitude) > 0.1))
            break;
    }
    return {amplitude_x, amplitude};
}
 
float measure_baseline_signal(TGraph *gTarget, int max_point = 3000)
{
    auto average_baseline = 0.;
    for (auto iPnt = 1; iPnt <= max_point; iPnt++)
    {
        average_baseline += gTarget->GetPointY(iPnt) / max_point;
    }
    return average_baseline;
}
 
float measure_baseline_timing(TGraph *gTarget)
{
    auto maximum_points = gTarget->GetN();
    auto average_baseline = 0.;
    for (auto iPnt = 1; iPnt <= maximum_points; iPnt++)
    {
        // average_baseline += gTarget->GetPointY(iPnt) / max_point;
    }
    return average_baseline;
}
 
float find_time_signal(TGraph *gTarget, float baseline_level, std::pair<float, float> amplitude_full)
{
    auto max_x = 0.;
    auto min_x = 0.;
    for (auto iPnt = amplitude_full.first; iPnt > 0; iPnt--)
    {
        auto current_x = gTarget->GetPointX(iPnt);
        auto current_y = gTarget->GetPointY(iPnt);
        if ((current_y < amplitude_full.second * 0.8) && (max_x == 0))
            max_x = current_x;
        if ((current_y < amplitude_full.second * 0.2) && (min_x == 0))
            min_x = current_x;
        if ((min_x != 0) && (max_x != 0))
            break;
    }
    gTarget->Fit("pol1", "Q", "", min_x, max_x);
    auto pol1 = gTarget->GetFunction("pol1");
    auto q = pol1->GetParameter(0);
    auto m = pol1->GetParameter(1);
    auto timing_x = (baseline_level - q) / m;
    return timing_x;
}
 
float find_time_timing(TGraph *gTarget, float baseline_level, std::pair<float, float> amplitude_full)
{
    auto max_x = 0.;
    auto min_x = 0.;
    for (auto iPnt = amplitude_full.first; iPnt > 0; iPnt--)
    {
        auto current_x = gTarget->GetPointX(iPnt);
        auto current_y = gTarget->GetPointY(iPnt);
        if ((current_y > amplitude_full.second * 0.8) && (max_x == 0))
            max_x = current_x;
        if ((current_y > amplitude_full.second * 0.2) && (min_x == 0))
            min_x = current_x;
        if ((min_x != 0) && (max_x != 0))
            break;
    }
    gTarget->Fit("pol1", "Q", "", min_x, max_x);
    auto pol1 = gTarget->GetFunction("pol1");
    auto q = pol1->GetParameter(0);
    auto m = pol1->GetParameter(1);
    auto timing_x = (baseline_level - q) / m;
    return timing_x;
}
 
float find_time_signal_1(TGraph *gTarget, float baseline_level, std::pair<float, float> amplitude_full)
{
    for (auto iPnt = amplitude_full.first; iPnt > 0; iPnt--)
    {
        auto current_x = gTarget->GetPointX(iPnt);
        auto current_y = gTarget->GetPointY(iPnt);
        if ((current_y < amplitude_full.second * 0.5))
            return current_x;
    }
    return -1.;
}
 
float find_time_timing_1(TGraph *gTarget, float baseline_level, std::pair<float, float> amplitude_full)
{
    for (auto iPnt = amplitude_full.first; iPnt > 0; iPnt--)
    {
        auto current_x = gTarget->GetPointX(iPnt);
        auto current_y = gTarget->GetPointY(iPnt);
        if ((current_y > amplitude_full.second * 0.5))
            return current_x;
    }
    return -1.;
}
 
void tgraphs_to_ttree(std::string input_file, int vbias = 0, std::string output_file = "ttree.root")
{
    set_thresholds();
    std::map<int, TGraph *> _TimingGraphs;
    std::map<int, TGraph *> _SignalGraphs;
    std::map<std::string, TH1F *> _TH1F;
    _TH1F["hDelta"] = new TH1F("hDelta", "hDelta", 2000, 50, 100);
    _TH1F["hDelta_1"] = new TH1F("hDelta_1", "hDelta_1", 2000, 50, 100);
    auto timing_channel_graphs = -1;
    auto signal_channel_graphs = -1;
    TFile *fin = new TFile(input_file.c_str());
    while (true)
    {
        timing_channel_graphs++;
        signal_channel_graphs++;
        std::string timing_label = "timing_" + std::to_string(timing_channel_graphs);
        std::string signal_label = "signal_" + std::to_string(signal_channel_graphs);
        _TimingGraphs[timing_channel_graphs] = (TGraph *)(fin->Get(timing_label.c_str()));
        _SignalGraphs[signal_channel_graphs] = (TGraph *)(fin->Get(signal_label.c_str()));
        if (!_TimingGraphs[timing_channel_graphs] || !_SignalGraphs[signal_channel_graphs])
            break;
    }
    //  Output files
    TFile *fout = new TFile(output_file.c_str(), "RECREATE");
    float signal_amplitude = 0.;
    float timing_amplitude = 0.;
    float signal_time = 0.;
    float timing_time = 0.;
    float signal_time_1 = 0.;
    float timing_time_1 = 0.;
    float delta = 0.;
    float delta_1 = 0.;
    TTree *tree = new TTree("tree", "tree");
    tree->Branch("signal_amplitude", &signal_amplitude, "signal_amplitude/F");
    tree->Branch("timing_amplitude", &timing_amplitude, "timing_amplitude/F");
    tree->Branch("signal_time", &signal_time, "signal_time/F");
    tree->Branch("timing_time", &timing_time, "timing_time/F");
    tree->Branch("signal_time_1", &signal_time_1, "signal_time_1/F");
    tree->Branch("timing_time_1", &timing_time_1, "timing_time_1/F");
    tree->Branch("delta", &delta, "delta/F");
    tree->Branch("delta_1", &delta_1, "delta_1/F");
    int iTer = -1;
    for (auto not_used : _TimingGraphs)
    {
        iTer++;
        auto timing_graph = _TimingGraphs[iTer];
        auto signal_graph = _SignalGraphs[iTer];
        if (!timing_graph || !signal_graph)
            continue;
        //  Amplitude of signal
        auto amplitude_full = find_first_amplitude(signal_graph);
        auto baseline_level = measure_baseline_signal(signal_graph);
        signal_amplitude = amplitude_full.second - baseline_level;
        //  Cut on 1-pe
         if (fabs(amplitude_full.second - mean_1pe[vbias]) > 10 * sigma_1pe[vbias])
            continue;
        //  Amplitude of timing
        auto baseline_level_timing = measure_baseline_timing(timing_graph);
        auto amplitude_full_timing = find_first_amplitude(timing_graph, 0, false);
        timing_amplitude = amplitude_full_timing.second - baseline_level_timing;
        //  Cut on good timing
        if (amplitude_full_timing.second > -0.6)
            continue;
        //  Timing of signal
        signal_time = find_time_signal(signal_graph, baseline_level, amplitude_full);
        signal_time_1 = find_time_signal_1(signal_graph, baseline_level, amplitude_full);
        //  Timing of timing_1
        timing_time = find_time_timing(timing_graph, baseline_level_timing, amplitude_full_timing);
        timing_time_1 = find_time_timing_1(timing_graph, baseline_level_timing, amplitude_full_timing);
        //  Delta
        delta = signal_time - timing_time;
        delta_1 = signal_time_1 - timing_time_1;
        //  Histogram
        _TH1F["hDelta"]->Fill(delta*1.e9);
        _TH1F["hDelta_1"]->Fill(delta_1*1.e9);
        //  Tree fill
        tree->Fill();
    }
    tree->Write();
    for (auto [label, histo] : _TH1F)
        histo->Write();
    fout->Close();
}