sEPDChannelAnalysis/macro/analyze_channels.C

0001 #include <TFile.h>
0002 #include <TH1F.h>
0003 #include <TTree.h>
0004 #include <TSystemDirectory.h>
0005 #include <TSystemFile.h>
0006 #include <TF1.h>
0007 #include <vector>
0008 #include <string>
0009 #include <regex>
0010
0011 const int NUM_CHANNELS = 744;
0012 const double FIT_MIN = 50;
0013 const double FIT_MAX = 1000;
0014 const bool verbose = false;
0015
0016 struct ChannelData {
0017     float max_mpv = -1;
0018     float min_mpv = 10000;
0019     TH1F* max_hist = nullptr;
0020     TH1F* min_hist = nullptr;
0021     TF1* minFitFunc = nullptr;
0022     TF1* maxFitFunc = nullptr;
0023 };
0024
0025 // Lorentzian fit function
0026 Double_t mipFitFunction(Double_t* x, Double_t* par) {
0027     return par[0] * (par[1]/TMath::Pi()) * (1/((x[0]-par[2])*(x[0]-par[2]) + par[1]*par[1]));
0028 }
0029
0030 // Extract run number from filename
0031 /*
0032 int extract_run_number(const std::string& filename) {
0033     std::regex pattern("run_(\\d+)\\.root_histos\\.root");
0034     std::smatch match;
0035     if (std::regex_search(filename, match, pattern)) {
0036         return std::stoi(match[1]);
0037     }
0038     return -1;
0039 }*/
0040
0041
0042 int extract_run_number(const std::string& filename) {
0043     // Enhanced pattern to match different file naming conventions
0044     std::regex pattern("run[_-]?(\\d+).*\\.root");  // Matches:
0045                                                      // run123.root
0046                                                      // run_123_histos.root
0047                                                      // RUN-123.root
0048     std::smatch match;
0049
0050     if (std::regex_search(filename, match, pattern)) {
0051         try {
0052             return std::stoi(match[1].str());
0053         } catch (...) {
0054             std::cerr << "Invalid run number in: " << filename << std::endl;
0055             return -1;
0056         }
0057     }
0058     std::cerr << "No run number found in: " << filename << std::endl;
0059     return -1;
0060 }
0061
0062
0063 int count_root_files(const std::string& path) {
0064     int count = 0;
0065     TSystemDirectory dir("", path.c_str());
0066     TList* entries = dir.GetListOfFiles();
0067     TIter next(entries);
0068     TSystemFile* entry;
0069
0070     while ((entry = static_cast<TSystemFile*>(next()))) {
0071         std::string name = entry->GetName();
0072         if (name == "." || name == "..") continue;
0073
0074         std::string full_path = path + "/" + name;
0075
0076         if (entry->IsDirectory()) {
0077             count += count_root_files(full_path);
0078         } else if (name.find("_histos.root") != std::string::npos) {
0079             count++;
0080         }
0081     }
0082     delete entries;
0083     return count;
0084 }
0085
0086 void update_progress(int processed, int total) {
0087     const int bar_width = 50;
0088     float progress = static_cast<float>(processed)/total;
0089     int pos = static_cast<int>(bar_width * progress);
0090
0091     std::cerr << "\r[";
0092     for(int i = 0; i < bar_width; ++i) {
0093         if(i < pos) std::cerr << "=";
0094         else if(i == pos) std::cerr << ">";
0095         else std::cerr << " ";
0096     }
0097     std::cerr << "] " << int(progress * 100.0) << "% "
0098               << "(" << processed << "/" << total << ")";
0099     std::cerr.flush();
0100 }
0101
0102 // Recursive directory scanner
0103 void process_directory(const std::string& path,
0104                       std::vector<ChannelData>& channels,
0105                       TTree* data_tree,
0106                       int& run_number,
0107                       std::vector<float>& current_mpv,
0108                       int& processed_files,
0109                       int total_files) {
0110     TSystemDirectory dir("", path.c_str());
0111     TList* entries = dir.GetListOfFiles();
0112     TIter next(entries);
0113     TSystemFile* entry;
0114
0115     while ((entry = static_cast<TSystemFile*>(next()))) {
0116         std::string name = entry->GetName();
0117         if (name == "." || name == "..") continue;
0118
0119         std::string full_path = path + "/" + name;
0120
0121         if (entry->IsDirectory()) {
0122             process_directory(full_path, channels, data_tree, run_number,
0123                              current_mpv, processed_files, total_files);
0124         }
0125         else if (name.find("_histos.root") != std::string::npos) {
0126             TFile* file = TFile::Open(full_path.c_str(), "READ");
0127             if (!file || file->IsZombie()) {
0128                 std::cerr << "\nError opening: " << full_path << std::endl;
0129                 continue;
0130             }
0131
0132             run_number = extract_run_number(name);
0133
0134             if(verbose) {
0135                 std::cout << "Run number = " << run_number  << std::endl;
0136             }
0137
0138             std::fill(current_mpv.begin(), current_mpv.end(), -1.0f);
0139
0140
0141             for (int ch = 0; ch < NUM_CHANNELS; ch++) {
0142                 TH1F* hist = dynamic_cast<TH1F*>(file->Get(Form("channel%d", ch)));
0143                 if (!hist || hist->GetEntries() < 100) continue;
0144
0145
0146                 //TF1 fitfunc("mipfit", mipFitFunction, FIT_MIN, FIT_MAX, 3);
0147                 //try a landau distribution
0148                 TF1 fitfunc("mipfit","landau",FIT_MIN,FIT_MAX);
0149                 fitfunc.SetParameters(2200, 100, 100);
0150                 hist->Fit(&fitfunc, "QRN0");
0151                 //current_mpv[ch] = fitfunc.GetParameter(2);
0152                 current_mpv[ch] = fitfunc.GetParameter(1);
0153
0154
0155                 //float mpv = fitfunc.GetParameter(2);
0156                 float mpv = fitfunc.GetParameter(1);
0157
0158                 if(verbose){
0159                     std::cout << "MPV for channel " << ch << " = " << mpv << std::endl;
0160                 }
0161
0162
0163                 // Clone and manage histograms properly
0164                 TH1F* hist_clone = static_cast<TH1F*>(hist->Clone());
0165                 hist_clone->SetDirectory(0);
0166
0167                 TF1* fit_clone = static_cast<TF1*>(fitfunc.Clone(
0168                     Form("ch%d_fit_clone", ch)));
0169
0170
0171                 if (mpv > channels[ch].max_mpv) {
0172                     delete channels[ch].max_hist;
0173                     delete channels[ch].maxFitFunc;
0174
0175                     channels[ch].max_mpv = mpv;
0176                     channels[ch].max_hist = hist_clone;
0177                     channels[ch].maxFitFunc = fit_clone;
0178                 } else if (mpv < channels[ch].min_mpv) {
0179                     delete channels[ch].min_hist;
0180                     delete channels[ch].minFitFunc;
0181
0182                     channels[ch].min_mpv = mpv;
0183                     channels[ch].min_hist = hist_clone;
0184                     channels[ch].minFitFunc = fit_clone;
0185                 } else {
0186                     delete hist_clone; // Only keep needed clones
0187                     delete fit_clone;
0188                 }
0189             }
0190
0191             data_tree->Fill();
0192
0193             file->Close();
0194             delete file;
0195
0196             processed_files++;
0197             update_progress(processed_files, total_files);
0198         }
0199     }
0200     // Do NOT delete entries - ROOT manages this memory
0201 }
0202
0203
0204 void analyze_channels() {
0205     std::vector<ChannelData> channels(NUM_CHANNELS);
0206     TFile output("/sphenix/user/ecroft/channel_analysis_landau.root", "RECREATE");
0207
0208     // Create main data tree
0209     TTree* data_tree = new TTree("rundata", "Per-run channel MPV values");
0210     int run_number;
0211     std::vector<float> current_mpv(NUM_CHANNELS, -1.0f);
0212
0213     data_tree->Branch("run", &run_number);
0214     data_tree->Branch("mpv_values", &current_mpv);
0215
0216     // Count files first for progress tracking
0217     const std::string input_dir = "/sphenix/tg/tg01/bulk/ecroft/sEPD/test/";
0218     int total_files = count_root_files(input_dir);
0219     int processed_files = 0;
0220
0221     if(total_files == 0) {
0222         std::cerr << "No files found!" << std::endl;
0223         return;
0224     }
0225
0226     std::cerr << "Processing " << total_files << " files..." << std::endl;
0227
0228     // Process files recursively from main directory
0229     process_directory(input_dir, channels, data_tree, run_number,
0230                      current_mpv, processed_files, total_files);
0231
0232     std::cerr << std::endl;
0233     // Save results
0234     output.cd();
0235     data_tree->Write();
0236
0237     // Create summary histograms
0238     TH1F h_maxmpv("max_mpv", "Maximum MPV per channel", NUM_CHANNELS, 0, NUM_CHANNELS);
0239     TH1F h_minmpv("min_mpv", "Minimum MPV per channel", NUM_CHANNELS, 0, NUM_CHANNELS);
0240
0241     // Save extremal histograms and fill summary
0242     for (int ch = 0; ch < NUM_CHANNELS; ch++) {
0243         h_maxmpv.SetBinContent(ch+1, channels[ch].max_mpv);
0244         h_minmpv.SetBinContent(ch+1, channels[ch].min_mpv);
0245
0246         if (channels[ch].max_hist) {
0247             channels[ch].max_hist->Write(Form("ch%d_max_hist", ch));
0248             delete channels[ch].max_hist;
0249         }
0250         if (channels[ch].min_hist) {
0251             channels[ch].min_hist->Write(Form("ch%d_min_hist", ch));
0252             delete channels[ch].min_hist;
0253         }
0254
0255         if (channels[ch].minFitFunc) {
0256             channels[ch].minFitFunc->Write(Form("ch%d_min_fit", ch));
0257             delete channels[ch].minFitFunc;
0258         }
0259
0260         if (channels[ch].maxFitFunc) {
0261             channels[ch].maxFitFunc->Write(Form("ch%d_max_fit", ch));
0262             delete channels[ch].maxFitFunc;
0263         }
0264     }
0265
0266     h_maxmpv.Write();
0267     h_minmpv.Write();
0268     output.Close();
0269 }