mbd/MBDTRKZ/plot_mbdtrkzdiff.C

0001 // plot_mbdtrkzdiff.C
0002 //
0003 // Reads a list of mbdtrk_vertex.root files, fits the highest peak in
0004 // h_mbdtrkz with a Gaussian, and plots the mean vs run number.
0005 // Run number is extracted from the parent directory name
0006 //
0007 // Usage:
0008 //   root -b -q 'plot_mbdtrkzdiff.C("MBDTRKZ/f.list")'
0009
0010 #include <fstream>
0011 #include <iostream>
0012 #include <string>
0013 #include <vector>
0014 #include <format>
0015
0016 #include <TCanvas.h>
0017 #include <TF1.h>
0018 #include <TFile.h>
0019 #include <TGraphErrors.h>
0020 #include <TH1F.h>
0021 #include <TLatex.h>
0022 #include <TStyle.h>
0023
0024 // Extract the run number from a path of the form .../MBDTRKZ/<run>/<seg>/filename
0025 // Returns -1 on failure.
0026 int run_from_path(const std::string &path)
0027 {
0028   // find last two '/' to get the directory name
0029   size_t slash2 = path.rfind('/');
0030   if (slash2 == std::string::npos) return -1;
0031   size_t slash1 = path.rfind('/', slash2 - 1);
0032   std::string rundir = (slash1 != std::string::npos)
0033                        ? path.substr(slash1 + 1, slash2 - slash1 - 1)
0034                        : path.substr(0, slash2);
0035   size_t us = rundir.find('_');
0036   if (us == std::string::npos) return -1;
0037   try { return std::stoi(rundir.substr(0, us)); }
0038   catch (...) { return -1; }
0039 }
0040
0041 void plot_mbdtrkzdiff(const std::string &filelist = "MBDTRKZ/f.list")
0042 {
0043   gStyle->SetOptStat(0);
0044   gStyle->SetOptFit(111111);
0045
0046   // --- read file list ---
0047   std::vector<std::string> files;
0048   {
0049     std::ifstream ifs(filelist);
0050     if (!ifs)
0051     {
0052       std::cerr << "Cannot open " << filelist << std::endl;
0053       return;
0054     }
0055     std::string line;
0056     while (std::getline(ifs, line))
0057     {
0058       if (!line.empty()) files.push_back(line);
0059     }
0060   }
0061
0062   if (files.empty())
0063   {
0064     std::cerr << "No files in " << filelist << std::endl;
0065     return;
0066   }
0067
0068   auto *g_mbdtrkzdiff = new TGraphErrors();
0069   g_mbdtrkzdiff->SetName("g_mbdtrkzdiff");
0070   g_mbdtrkzdiff->SetTitle("MBD - Tracker z-vertex;run;#Deltaz [cm]");
0071
0072   int ipt = 0;
0073   for (const auto &fname : files)
0074   {
0075     int runnumber = run_from_path(fname);
0076     if (runnumber < 0)
0077     {
0078       std::cerr << "WARNING: cannot parse run number from: " << fname << std::endl;
0079       continue;
0080     }
0081
0082     TFile *f = TFile::Open(fname.c_str(), "READ");
0083     if (!f || f->IsZombie())
0084     {
0085       std::cerr << "WARNING: cannot open " << fname << std::endl;
0086       continue;
0087     }
0088
0089     auto *h = dynamic_cast<TH1F *>(f->Get("h_mbdtrkz"));
0090     if (!h)
0091     {
0092       std::cerr << "WARNING: h_mbdtrkz not found in " << fname << std::endl;
0093       f->Close();
0094       continue;
0095     }
0096     h->SetDirectory(nullptr);
0097     f->Close();
0098
0099     if (h->GetEntries() == 0)
0100     {
0101       std::cerr << "WARNING: h_mbdtrkz empty in " << fname << std::endl;
0102       delete h;
0103       continue;
0104     }
0105
0106     // locate highest peak
0107     int    peak_bin = h->GetMaximumBin();
0108     double peak_x   = h->GetBinCenter(peak_bin);
0109
0110     // first pass: rough Gaussian fit ±3 cm around peak
0111     TF1 gfit("gfit", "gaus", peak_x - 3.0, peak_x + 3.0);
0112     int status = h->Fit(&gfit, "RQ0");
0113     if (status != 0)
0114     {
0115       std::cerr << "WARNING: initial fit failed for run " << runnumber << std::endl;
0116       delete h;
0117       continue;
0118     }
0119
0120     double sigma = gfit.GetParameter(2);
0121     if (sigma <= 0)
0122     {
0123       std::cerr << "WARNING: bad sigma from initial fit for run " << runnumber << std::endl;
0124       delete h;
0125       continue;
0126     }
0127
0128     // second pass: refit within ±2σ of first-pass mean
0129     double mean1 = gfit.GetParameter(1);
0130     gfit.SetRange(mean1 - 2.0 * sigma, mean1 + 2.0 * sigma);
0131     status = h->Fit(&gfit, "RQ0");
0132     if (status != 0)
0133     {
0134       std::cerr << "WARNING: second fit failed for run " << runnumber << std::endl;
0135       delete h;
0136       continue;
0137     }
0138
0139     double mean     = gfit.GetParameter(1);
0140     double mean_err = gfit.GetParError(1);
0141
0142     g_mbdtrkzdiff->SetPoint(ipt, runnumber, mean);
0143     g_mbdtrkzdiff->SetPointError(ipt, 0., mean_err);
0144     ++ipt;
0145
0146     delete h;
0147   }
0148
0149   if (g_mbdtrkzdiff->GetN() == 0)
0150   {
0151     std::cerr << "No points in graph — nothing to plot." << std::endl;
0152     return;
0153   }
0154
0155   g_mbdtrkzdiff->Sort();
0156
0157   // pol0 fit
0158   TF1 *fpol0 = new TF1("fpol0", "pol0", g_mbdtrkzdiff->GetX()[0]-10, g_mbdtrkzdiff->GetX()[g_mbdtrkzdiff->GetN() - 1]+10);
0159   fpol0->SetLineColor(kRed);
0160   fpol0->SetLineWidth(2);
0161   g_mbdtrkzdiff->Fit(fpol0, "QR");
0162
0163   // draw
0164   auto *c = new TCanvas("c_mbdtrkzdiff", "MBD-Tracker dz", 900, 600);
0165   c->SetLeftMargin(0.12);
0166   c->SetBottomMargin(0.12);
0167
0168   g_mbdtrkzdiff->SetMarkerStyle(20);
0169   g_mbdtrkzdiff->SetMarkerSize(0.8);
0170   g_mbdtrkzdiff->Draw("AP");
0171
0172   fpol0->Draw("same");
0173
0174   TLatex tex;
0175   tex.SetNDC();
0176   tex.SetTextSize(0.035);
0177   tex.DrawLatex(0.15, 0.85, std::format("pol0 fit: {:.3f} #pm {:.3f} cm", fpol0->GetParameter(0), fpol0->GetParError(0)).c_str() );
0178
0179   c->SaveAs("mbdtrkzdiff.pdf");
0180   std::cout << "Saved mbdtrkzdiff.pdf" << std::endl;
0181 }