sEPDChannelAnalysis/macro/plot_channel_mpvs_RooFit.C

0001 #include <TFile.h>
0002 #include <TDirectory.h>
0003 #include <TGraph.h>
0004 #include <TCanvas.h>
0005 #include <TF1.h>
0006 #include <TFitResult.h>
0007 #include <TPaveText.h>
0008 #include <vector>
0009 #include <string>
0010 #include <regex>
0011 #include <TSystemDirectory.h>
0012 #include <TSystemFile.h>
0013
0014 void plot_channel_mpvs_fit(const char* input_dir = "/sphenix/tg/tg01/bulk/ecroft/sEPD/analysis_output",
0015                            const char* output_filename = "channel_mpvs.pdf") {
0016     // Data storage: vector per channel containing (run, mpv) pairs
0017     const int num_channels = 744;
0018     std::vector<std::vector<std::pair<int, float>>> channel_data(num_channels);
0019
0020     // Get list of run files
0021     TSystemDirectory dir("input_dir", input_dir);
0022     TList* files = dir.GetListOfFiles();
0023     TIter next(files);
0024     TSystemFile* file;
0025     int processed_runs = 0;
0026
0027     std::regex run_pattern("run_(\\d+)_analysis.root");
0028
0029     while ((file = (TSystemFile*)next())) {
0030         TString fname = file->GetName();
0031         std::string filename = fname.Data();
0032
0033         // Match run files
0034         std::smatch match;
0035         if (!std::regex_search(filename, match, run_pattern)) continue;
0036
0037         int run_number = std::stoi(match[1].str());
0038         TString full_path = Form("%s/%s", input_dir, fname.Data());
0039
0040         // Open run file
0041         TFile* run_file = TFile::Open(full_path, "READ");
0042         if (!run_file || run_file->IsZombie()) continue;
0043
0044         // Get run directory
0045         TDirectory* run_dir = (TDirectory*)run_file->Get(Form("run_%d", run_number));
0046         if (!run_dir) {
0047             run_file->Close();
0048             continue;
0049         }
0050
0051         // Process channels
0052         TIter ch_next(run_dir->GetListOfKeys());
0053         TKey* ch_key;
0054         while ((ch_key = (TKey*)ch_next())) {
0055             if (!TString(ch_key->GetName()).BeginsWith("ch_")) continue;
0056
0057             // Extract channel number
0058             TString ch_name = ch_key->GetName();
0059             int channel = TString(ch_name(3, 4)).Atoi();
0060
0061             // Get fit parameters
0062             TDirectory* ch_dir = (TDirectory*)ch_key->ReadObj();
0063             TVectorD* params = (TVectorD*)ch_dir->Get("fit_parameters");
0064             if (!params || params->GetNoElements() < 1) continue;
0065
0066             double mpv = (*params)[0];
0067             if (mpv > 0) {
0068                 channel_data[channel].emplace_back(run_number, mpv);
0069             }
0070         }
0071
0072         run_file->Close();
0073         delete run_file;
0074
0075         // Progress
0076         if (++processed_runs % 10 == 0) {
0077             std::cout << "Processed " << processed_runs << " runs ("
0078                       << channel_data[0].size() << " entries for channel 0)" << std::endl;
0079         }
0080     }
0081
0082     // Create plots
0083     TCanvas canvas("canvas", "Channel MPV Plots", 1200, 800);
0084     canvas.Print(Form("%s[", output_filename));
0085
0086     for (int ch = 0; ch < num_channels; ch++) {
0087         const auto& data = channel_data[ch];
0088         if (data.empty()) continue;
0089
0090         std::vector<double> runs, mpvs;
0091         for (const auto& point : data) {
0092             runs.push_back(point.first);
0093             mpvs.push_back(point.second);
0094         }
0095
0096         TGraph graph(runs.size(), runs.data(), mpvs.data());
0097         graph.SetTitle(Form("Channel %d MPV History;Run Number;MPV [ADC]", ch));
0098         graph.SetMarkerStyle(20);
0099         graph.SetMarkerSize(0.8);
0100         graph.SetMarkerColor(kBlue);
0101
0102         bool has_fit = false;
0103         double intercept = 0.0, slope = 0.0;
0104         double intercept_err = 0.0, slope_err = 0.0;
0105
0106         if (graph.GetN() >= 2) {
0107             // Perform linear fit
0108             TFitResultPtr fit_result = graph.Fit("pol1", "QS");
0109             if (fit_result->Status() == 0) {
0110                 has_fit = true;
0111                 intercept = fit_result->Parameter(0);
0112                 intercept_err = fit_result->ParError(0);
0113                 slope = fit_result->Parameter(1);
0114                 slope_err = fit_result->ParError(1);
0115             }
0116         }
0117
0118
0119         canvas.Clear();
0120         graph.Draw("AP");
0121
0122         if (has_fit) {
0123             // Draw the fit function
0124             TF1 *fit_func = graph.GetFunction("pol1");
0125             if (fit_func) {
0126                 fit_func->SetLineColor(kRed);
0127                 fit_func->Draw("same");
0128             }
0129
0130             // Add text box with fit parameters
0131             TPaveText *pt = new TPaveText(0.15, 0.7, 0.45, 0.9, "NDC");
0132             pt->SetFillColor(0);
0133             pt->SetBorderSize(1);
0134             pt->SetTextAlign(12);
0135             pt->SetTextSize(0.03);
0136             pt->AddText(Form("Intercept: %.2f #pm %.2f", intercept, intercept_err));
0137             pt->AddText(Form("Slope: %.4f #pm %.4f", slope, slope_err));
0138             pt->Draw();
0139         }
0140
0141         canvas.Update();
0142
0143
0144         canvas.Print(output_filename);
0145
0146         if (ch % 50 == 0) {
0147             std::cout << "Plotted " << ch << "/" << num_channels
0148                       << " channels (" << (ch*100/num_channels) << "%)" << std::endl;
0149         }
0150     }
0151
0152     canvas.Print(Form("%s]", output_filename));
0153
0154     // Cleanup
0155     //input_file->Close();
0156     //delete input_file;
0157     std::cout << "Saved plots to: " << output_filename << std::endl;
0158 }