mbd/PLOTS/CalAnalyzer.C

0001 #include <iostream>
0002 #include <iomanip>
0003 #include <fstream>
0004 #include <vector>
0005 #include <map>
0006 #include <string>
0007 #include "TH1F.h"
0008 #include "TH1D.h"
0009 #include "TGraph.h"
0010 #include "TCanvas.h"
0011
0012 class CalAnalyzer {
0013 private:
0014
0015   const int NPMT{128};
0016   int NPAR{9};
0017
0018   std::string run_name;
0019
0020   struct CalibrationStats {
0021         double mean;
0022         double rms;
0023         int count;
0024   };
0025
0026   // Data storage: [RunID][ChannelIdx]
0027   std::map<int, std::vector<std::vector<double>>> dataStore;
0028
0029   // Stats: [ParamIdx][ChannelIdx]
0030   std::vector<std::vector<CalibrationStats>> cleanStats;
0031
0032   // Plots: [ParamIdx][ChannelIdx]
0033   std::vector<std::vector<TH1F*>> hists;
0034   std::vector<std::vector<TGraph*>> graphs;
0035
0036 public:
0037   explicit CalAnalyzer(std::string_view input_name = "none")
0038   : run_name{input_name}
0039   {
0040     if ( run_name.find("auau") != std::string::npos )
0041     {
0042       NPAR = 5;
0043     }
0044
0045     // Initialize 2D vectors for hists and graphs: NPAR parameters x NPMT channels
0046     hists.resize(NPAR, std::vector<TH1F*>(NPMT));
0047     graphs.resize(NPAR, std::vector<TGraph*>(NPMT));
0048     cleanStats.resize(NPAR, std::vector<CalibrationStats>(NPMT));
0049
0050     TString name;
0051     TString title;
0052     for (int ch = 0; ch < NPMT; ch++)
0053     {
0054       for (int p = 0; p < NPAR; p++)
0055       {
0056         name = "h_ch"; name += ch; name += "_p"; name += p;
0057         title = "Dist: Param "; title += p; title += ", Ch "; title += ch; title += ";Value;Entries";
0058         hists[p][ch] = new TH1F(name, title, 200, 0., -1.);
0059
0060         graphs[p][ch] = new TGraph();
0061         name = "g_ch"; name += ch; name += "_p"; name += p;
0062         graphs[p][ch]->SetName( name );
0063         title = "Param "; title += p; title += ", Ch "; title += ch; title += " vs Run;Run ID;Value";
0064         graphs[p][ch]->SetTitle( title );
0065         graphs[p][ch]->SetMarkerStyle(20);
0066         graphs[p][ch]->SetMarkerSize(0.2);
0067       }
0068     }
0069   }
0070
0071   ~CalAnalyzer()
0072   {
0073     for (int ch = 0; ch < NPMT; ch++)
0074     {
0075       for (int p = 0; p < NPAR; p++)
0076       {
0077         delete hists[p][ch];
0078         delete graphs[p][ch];
0079       }
0080     }
0081   }
0082
0083   CalAnalyzer(const CalAnalyzer&) = delete;
0084   CalAnalyzer& operator=(const CalAnalyzer&) = delete;
0085   CalAnalyzer(CalAnalyzer&&) = delete;
0086   CalAnalyzer& operator=(CalAnalyzer&&) = delete;
0087
0088
0089   void LoadRun(int runID, const std::string& filename)
0090   {
0091     std::ifstream file(filename);
0092     if (!file.is_open()) {
0093       std::cerr << "Error: Could not open " << filename << std::endl;
0094       return;
0095     }
0096
0097     std::vector<std::vector<double>> channels;
0098     channels.resize(NPAR, std::vector<double>(NPMT));
0099     int pmtch{-1};
0100     while (file >> pmtch)
0101     {
0102       //cout << runID << "\t" << pmtch << endl;
0103       if (pmtch >= 0 && pmtch < NPMT)
0104       {
0105         // other stuff
0106         double seed_threshold, seed_minrej, seed_natpeak, seed_maxrej, seed_qmax; // NOLINT
0107         file >> seed_threshold >> seed_minrej >> seed_natpeak >> seed_maxrej >> seed_qmax;
0108
0109         //cout << "xxx " << pmtch << "\t" << seed_qmax << endl;
0110         for (int p = 0; p < NPAR; ++p) {
0111           double val;
0112           file >> val;
0113
0114           /*
0115              if ( p==1 && pmtch==0 )
0116              {
0117              cout << "xxx " << runID << "\t" << p << "\t" << pmtch << "\t" << val << endl;
0118              }
0119              */
0120
0121           channels[p][pmtch] = val;
0122
0123           // Update Histogram
0124           hists[p][pmtch]->Fill( val );
0125
0126           // Update Graph
0127           int nextPt = graphs[p][pmtch]->GetN();
0128           // should add second axis with runnums
0129           //graphs[p][pmtch]->SetPoint(nextPt, (double)runID, val);
0130           graphs[p][pmtch]->SetPoint(nextPt, (double)nextPt, val);
0131         }
0132       }
0133     }
0134     dataStore[runID] = channels;
0135
0136
0137     file.close();
0138   }
0139
0140   void CalculateCleanStats(double sigmaEquivalent = 3.0)
0141   {
0142     // Set Min and Max of graphs
0143     for (int ipmt=0; ipmt<NPMT; ipmt++)
0144     {
0145       for (int ipar=0; ipar<NPAR; ipar++)
0146       {
0147         CalculateCleanStats(ipar,ipmt,sigmaEquivalent);
0148       }
0149     }
0150   }
0151
0152   /**
0153    * Calculates Mean and RMS after removing outliers
0154    * @param p parameter index (0-NPAR)
0155    * @param ch pmtch (0-NPMT)
0156    * @param sigmaEquivalent number of standard deviations to keep (usually 3.0)
0157    */
0158   void CalculateCleanStats(int p, int ch, double sigmaEquivalent = 3.0)
0159   {
0160     //cout << "In CalculateCleanStats " << ch << "\t" << p << endl;
0161     if (p < 0 || p >= NPAR || ch < 0 || ch >= NPMT) return;
0162
0163     std::vector<double> values;
0164     for (auto const& [run, runData] : dataStore)
0165     {
0166       if ( std::isnan(runData[p][ch]) )
0167       {
0168         continue;
0169       }
0170
0171       values.push_back(runData[p][ch]);
0172     }
0173
0174     if (values.empty()) return;
0175
0176     // 1. Find the Median
0177     std::sort(values.begin(), values.end());
0178     double middle = values[values.size() / 2];
0179
0180     // 2. Find the Median Absolute Deviation (MAD)
0181     std::vector<double> diffs;
0182     for (double v : values)
0183     {
0184       diffs.push_back(std::abs(v - middle));  // NOLINT
0185     }
0186     std::sort(diffs.begin(), diffs.end());
0187     double mad = diffs[diffs.size() / 2];
0188
0189     // 3. Convert MAD to a "Sigma Equivalent"
0190     // For normal distributions, 1 sigma approx 1.4826 * MAD
0191     double robustSigma = 1.4826 * mad;
0192     if (robustSigma == 0) robustSigma = 1e-9; // Avoid division by zero
0193
0194     // 4. Pass 2: Calculate Mean/RMS only for points within the robust window
0195     double sum = 0;
0196     double sumSq = 0;
0197     int n = 0;
0198     for (double v : values) {
0199       if (std::abs(v - middle) < (sigmaEquivalent * robustSigma)) {
0200         sum += v;
0201         sumSq += v * v;
0202         n++;
0203       }
0204     }
0205
0206     if (n > 0) {
0207       double mean = sum / n;
0208       double rms = std::sqrt(std::abs((sumSq / n) - (mean * mean)));
0209       cleanStats[p][ch] = {mean, rms, n};
0210
0211       //cout << "aaa " << ch << "\t" << p << "\t" << middle << "\t" << mad << "\t" << mean << "\t" << rms << "\t" << n << endl;
0212
0213       // Update the to focus on the "Clean" area
0214       graphs[p][ch]->SetMinimum(mean - 12*rms);
0215       graphs[p][ch]->SetMaximum(mean + 10*rms);
0216
0217       delete hists[p][ch];
0218       TString name = "h_ch"; name += ch; name += "_"; name += p;
0219       TString title = "Dist: Param"; title += p; title += ", Ch "; title += ch; title += ";Value;Entries";
0220       hists[p][ch] = new TH1F(name, title, 200, mean-12*rms, mean+10*rms);
0221       for (double v : values)
0222       {
0223         hists[p][ch]->Fill( v );
0224       }
0225     }
0226
0227   }
0228
0229   CalibrationStats GetCleanStats(int ch, int p) { return cleanStats[p][ch]; }
0230
0231   /**
0232    * Saves the calculated means to a text file.
0233    * Format: ch  par0  par1  par2  ...  NPAR
0234    */
0235   void SaveCleanMeans(const std::string& outFilename = "temp_mipseeds.txt")
0236   {
0237     std::ofstream outFile(outFilename);
0238     if (!outFile.is_open()) {
0239       std::cerr << "Error: Could not open " << outFilename << " for writing." << std::endl;
0240       return;
0241     }
0242
0243     /*
0244     // Header
0245     outFile << "# ch";
0246     for (int p = 0; p < NPAR; ++p)
0247     {
0248       outFile << "      par" << p;
0249     }
0250     outFile << std::endl;
0251     */
0252
0253     // Set fixed precision for clean columns
0254     outFile << std::fixed << std::setprecision(6);
0255
0256     for (int ch = 0; ch < NPMT; ++ch) {
0257       outFile << std::setw(3) << ch; // Channel number column
0258
0259       for (int p = 0; p < NPAR; ++p) {
0260         // Access the clean mean we calculated using the MAD algorithm
0261         double meanVal = cleanStats[p][ch].mean;
0262
0263         // If stats haven't been calculated for this channel yet,
0264         // fallback to the raw histogram mean
0265         if (cleanStats[p][ch].count == 0) {
0266           meanVal = hists[p][ch]->GetMean();
0267         }
0268
0269         outFile << "  " << std::setw(10) << meanVal;
0270       }
0271       outFile << std::endl;
0272     }
0273
0274     outFile.close();
0275     std::cout << "Successfully saved clean parameters to: " << outFilename << std::endl;
0276   }
0277
0278   // Accessors
0279   double GetVal(int runID, int ch, int p)
0280   {
0281     return (dataStore.count(runID)) ? dataStore[runID][p][ch] : -1e9;
0282   }
0283
0284   TH1* GetHist(int ch, int p) { return hists[p][ch]; }
0285
0286   TGraph* GetGraph(int ch, int p) { return graphs[p][ch]; }
0287
0288   void DrawSummary(int ch, int p)
0289   {
0290     TString name = "c_"; name += ch; name += "_"; name += p;
0291     TCanvas *c = new TCanvas(name, "Calibration Summary", 800, 400);
0292     c->Divide(2, 1);
0293     c->cd(1);
0294     hists[p][ch]->Draw();
0295     c->cd(2);
0296     graphs[p][ch]->Draw("AP");
0297   }
0298
0299   /**
0300    * Generates a multi-page PDF for each parameter.
0301    * Each page contains the TH1 and TGraph for one channel.
0302    */
0303   void SaveParameterPDFs()
0304   {
0305     TCanvas *c = new TCanvas("c_print", "Printing Canvas", 1000, 500);
0306     c->Divide(2, 1);
0307
0308     // num mip fits out of range
0309     TString mipfitname = "mbd_"; mipfitname += run_name; mipfitname += "_mip.result";
0310     std::ofstream mipfile( mipfitname.Data() );
0311
0312     int par_mip = 6;
0313     if ( run_name.find("auau") != std::string::npos )
0314     {
0315       par_mip = 1;
0316     }
0317
0318     for (int p = 0; p < NPAR; ++p)
0319     {
0320       TString pdfname = "mbd_"; pdfname += run_name; pdfname += "_par"; pdfname += p; pdfname += ".pdf";
0321
0322       // Start the PDF file
0323       c->Print(pdfname + "[");
0324
0325       for (int ch = 0; ch < NPMT; ++ch)
0326       {
0327         // Left side: Histogram
0328         c->cd(1);
0329         hists[p][ch]->SetLineColor(kBlue);
0330         hists[p][ch]->Draw();
0331
0332         // Right side: TGraph (Run-by-Run)
0333         c->cd(2);
0334         graphs[p][ch]->SetMarkerColor(kRed);
0335         //graphs[p][ch]->SetDrawOption("AP");
0336         graphs[p][ch]->Draw("AP");
0337
0338         // Print current page
0339         c->Print(pdfname);
0340
0341         if ( p == par_mip )
0342         {
0343           int oflowbin = hists[p][ch]->GetNbinsX() + 1;
0344           mipfile << ch << "\t" << hists[p][ch]->GetBinContent(0) << "\t" << hists[p][ch]->GetBinContent(oflowbin) << std::endl;
0345         }
0346
0347       }
0348
0349       // Close the PDF file
0350       c->Print(pdfname + "]");
0351       std::cout << "Created: " << pdfname << std::endl;
0352
0353       if ( p==6 )
0354       {
0355         mipfile.close();
0356       }
0357     }
0358     delete c;
0359   }
0360 };