From a3c3fe645d38879857163d91fd2254649af4aab8 Mon Sep 17 00:00:00 2001 From: AndreaPareti Date: Mon, 13 Apr 2026 16:14:46 +0200 Subject: [PATCH 1/5] First tests with sim output digitization, including FERS trigger logic --- analysis/HidraAna.C | 11 +- analysis/HidraTB25Ana.C | 837 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 845 insertions(+), 3 deletions(-) create mode 100644 analysis/HidraTB25Ana.C diff --git a/analysis/HidraAna.C b/analysis/HidraAna.C index 65298e4..94146f0 100644 --- a/analysis/HidraAna.C +++ b/analysis/HidraAna.C @@ -42,13 +42,15 @@ void GetSiPMcoordinate(int TowID, int rowID, int colID_original, double &SiPM_X, double colID = (static_cast(grouping)-1)/2 + channel*grouping; if(fiber == "S"){ + //SiPM_X = +moduleX/2 - tuberadius - (tuberadius*2)*colID; + //SiPM_Y = -moduleY/2 + tuberadius + (sq3*tuberadius)*rowID+tuberadius*(2.*sq3m1-1.); SiPM_X = +moduleX/2 - tuberadius - (tuberadius*2)*colID; - SiPM_Y = -moduleY/2 + tuberadius + (sq3*tuberadius)*rowID+tuberadius*(2.*sq3m1-1.); + SiPM_Y = -moduleY/2 + tuberadius + (sq3*tuberadius)*rowID+tuberadius*(2.*sq3m1-1.); //std::cout << "S X: " << SiPM_X << " Y: " << SiPM_Y << std::endl; } if(fiber == "C"){ - SiPM_X = +moduleX/2 - tuberadius - tuberadius - (tuberadius*2)*colID ; + SiPM_X = +moduleX/2 + tuberadius - (tuberadius*2)*colID ; SiPM_Y = -moduleY/2 + tuberadius + (sq3*tuberadius)*rowID+tuberadius*(2.*sq3m1-1.); //std::cout << "C X: " << SiPM_X << " Y: " << SiPM_Y << std::endl; } @@ -172,8 +174,10 @@ void HidraAna(double energy, const string intup){ //SipmMapS->Fill( modcol[towID]*NofFiberscolumn + colID, modrow[towID]*NofFibersrow+rowID, content); SipmMapS->Fill( colID, towID*NofFibersrow+rowID, content); SciSiPMCoordinates->Fill(SiPM_X, SiPM_Y, content); + std::cout << "S channel: " << "x: " << SiPM_X << "\t y: " << SiPM_Y << "\t content: " << content << std::endl; } + for(unsigned int N=0; Nsize(); N++){ // Loop over SiPMs - S Fibers double content = CSiPM->at(N)/cerPheGeV; totcer+=content; @@ -185,6 +189,7 @@ void HidraAna(double energy, const string intup){ double SiPM_X, SiPM_Y; GetSiPMcoordinate(towID, rowID, colID, SiPM_X, SiPM_Y, "C", grouping); CerSiPMCoordinates->Fill(SiPM_X, SiPM_Y, content); + std::cout << "C channel: " << "x: " << SiPM_X << "\t y: " << SiPM_Y << "\t content: " << content << std::endl; } @@ -196,7 +201,7 @@ void HidraAna(double energy, const string intup){ leakene->Fill(lenergy/1000/energy); chidist->Fill((totsci-ecalo)/(totcer-ecalo)); //std::cout << "totSci: " << totsci << "\t totCer: " << totcer << std::endl; - //break; + break; } diff --git a/analysis/HidraTB25Ana.C b/analysis/HidraTB25Ana.C new file mode 100644 index 0000000..3f80140 --- /dev/null +++ b/analysis/HidraTB25Ana.C @@ -0,0 +1,837 @@ +// file: HidraAna.C + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "HidraGeo.h" +#include + +using json = nlohmann::json; + +// Global Constants and parameters +// +// Fallback per-channel thresholds used only if a FERS id is missing in fers_to_thr_map +const double FersThresholdS = 0.07; // GeV +const double FersThresholdC = 0.07; // GeV + +// Minimum visible signal per channel (>=) +// Set to 0.0 or 1.0 to disable the cut. +const double MinVisiblePheS = 2.0; // photoelectrons +const double MinVisiblePheC = 2.0; // photoelectrons + +// Scintillator on even rows (start from 0), Cerenkov on odd rows +const unsigned int grouping = 8; +// Write event display every N entries (0 to disable) +const unsigned int EventDisplayEvery = 100; +// A FERS is activated if at least this many channels exceed threshold +const unsigned int FersMultiplicity = 2; + +const double chi = 0.38; +//const double sciPheGeV = 119.001; // tb24 +//const double cerPheGeV = 29.4; // tb24 +const double sciPheGeV = 178.501; +const double cerPheGeV = 43; + +const double NofSipmCells_sci = 7772; +const double NofSipmCells_cer = 3443; +const double sci_pde = 0.22; +const double cer_pde = 0.38; + +const double elcont = 1.005; +const double picont = 1.028; + +// In this parametrised simulation, rawPhe is stored as "impinging optical photons" +// To account for multiple photons impinging on the same SiPM cell, +// we apply saturation correction to get "fired cells", which is the actual observable. +const bool ApplySaturation = true; + +double apply_sipm_saturation_from_pe(double pe, + double phePerGeV, + double nCells) +{ + if (pe <= 0.0 || phePerGeV <= 0.0 || nCells <= 0.0) { + return 0.0; + } + + const double firedCells = -nCells * std::expm1(-pe / nCells); + return firedCells / phePerGeV; +} + + + +std::map fers_to_thr_map = { + { 1, 0.17 }, + { 2, 0.20 }, + { 3, 0.30 }, + { 4, 0.28 }, + { 5, 0.105}, + { 6, 0.095}, + { 7, 0.08 }, + { 8, 0.08 }, + { 9, 0.05 }, + {10, 0.07 }, + {11, 0.05 }, + {12, 0.025}, + {13, 0.06 }, + {14, 0.05 }, + {15, 0.12 }, + {16, 0.06 } +}; + +struct SipmMapEntry { + int boardID = -1; + std::string type; + int row = -1; + int column = -1; + double x = 0.0; + double y = 0.0; + int fersId = -1; + int ch = -1; + std::string module_name; + double x_local = 0.0; + double y_local = 0.0; +}; + +struct FersChannelHit { + uint64_t key = 0; + double rawPhe = 0.0; + double signal = 0.0; +}; + +struct EventDisplayHit { + uint64_t key = 0; + double x = 0.0; + double y = 0.0; + double rawPhe = 0.0; + double signal = 0.0; +}; + +struct GaussianFitSummary { + double mean = std::numeric_limits::quiet_NaN(); + double meanErr = std::numeric_limits::quiet_NaN(); + double rms = std::numeric_limits::quiet_NaN(); + double rmsErr = std::numeric_limits::quiet_NaN(); +}; + +struct HistMomentSummary { + double mean = std::numeric_limits::quiet_NaN(); + double meanErr = std::numeric_limits::quiet_NaN(); + double rms = std::numeric_limits::quiet_NaN(); + double rmsErr = std::numeric_limits::quiet_NaN(); +}; + +using SipmLookup = std::unordered_map; +using FersKey = uint64_t; + +std::string make_sipm_key(const std::string& tower, + const std::string& type, + int row, + int column) +{ + std::ostringstream os; + os << tower << '|' << type << '|' << row << '|' << column; + return os.str(); +} + +FersKey make_fers_key(int boardID, int fersId) +{ + return (static_cast(static_cast(boardID)) << 32) | + static_cast(fersId); +} + +int grouped_column(int colID_original, unsigned int groupingValue) +{ + return static_cast(colID_original / groupingValue); +} + +int grouped_to_json_column(int groupedCol, unsigned int groupingValue) +{ + return groupedCol * static_cast(groupingValue) + static_cast(groupingValue) / 2; +} + +double GetFersThreshold(int fersId, double fallbackThreshold) +{ + if (fersId < 0) { + return fallbackThreshold; + } + + const auto it = fers_to_thr_map.find(static_cast(fersId)); + return (it != fers_to_thr_map.end()) ? it->second : fallbackThreshold; +} + +SipmLookup LoadSipmMap(const std::string& jsonPath) +{ + std::ifstream in(jsonPath); + if (!in) { + throw std::runtime_error("Cannot open SiPM map JSON: " + jsonPath); + } + + json j; + in >> j; + + SipmLookup lookup; + for (auto it = j.begin(); it != j.end(); ++it) { + const json& node = it.value(); + + SipmMapEntry entry; + entry.boardID = node.value("boardID", -1); + entry.type = node.value("type", ""); + entry.row = node.value("row", -1); + entry.column = node.value("column", -1); + entry.x = node.value("x", 0.0); + entry.y = node.value("y", 0.0); + entry.fersId = node.value("fersId", -1); + entry.ch = node.value("ch", -1); + entry.module_name = node.value("module_name", ""); + entry.x_local = node.value("x_local", 0.0); + entry.y_local = node.value("y_local", 0.0); + + lookup[make_sipm_key(entry.module_name, entry.type, entry.row, entry.column)] = entry; + } + + return lookup; +} + +const SipmMapEntry* FindSipmInfo(const SipmLookup& lookup, + const std::string& tower, + const std::string& type, + int row, + int groupedCol, + unsigned int /*groupingValue*/) +{ + const int jsonColumn = groupedCol; + const auto key = make_sipm_key(tower, type, row, jsonColumn); + const auto it = lookup.find(key); + return (it == lookup.end()) ? nullptr : &it->second; +} + +bool IsActivatedFers(const std::unordered_map& fersChannelsOverThreshold, + FersKey key) +{ + const auto it = fersChannelsOverThreshold.find(key); + return it != fersChannelsOverThreshold.end() && + it->second >= FersMultiplicity; +} + +std::string make_event_display_name(const std::string& prefix, unsigned int entry) +{ + std::ostringstream os; + os << prefix << "_evt" << std::setw(6) << std::setfill('0') << entry; + return os.str(); +} + +TH2F* CreateEventDisplayHist(const std::string& name, + const std::string& title, + unsigned int groupingValue) +{ + const int xBins = static_cast(NofSiPMTowersX * NofFiberscolumn / groupingValue); + const int yBins = NofSiPMTowersY * NofFibersrow; + + return new TH2F(name.c_str(), title.c_str(), + xBins, + -NofSiPMTowersX * moduleX / 2.0, NofSiPMTowersX * moduleX / 2.0, + yBins, + -NofSiPMTowersY * moduleY / 2.0, NofSiPMTowersY * moduleY / 2.0); +} + +void FillActivatedEventDisplay(TH2F* hist, + const std::vector& hits, + const std::unordered_map& fersChannelsOverThreshold, + double minVisiblePhe) +{ + if (!hist) { + return; + } + + for (const auto& hit : hits) { + if (IsActivatedFers(fersChannelsOverThreshold, hit.key) && + hit.rawPhe >= minVisiblePhe) { + hist->Fill(hit.x, hit.y, hit.signal); + } + } +} + +void GetSiPMcoordinate(int TowID, + int rowID, + int colID_original, + double& SiPM_X, + double& SiPM_Y, + std::string fiber, + unsigned int groupingValue) +{ + double TowerOffsetY = -((NoModulesSiPM - 1) * moduleY) / 2 + TowID * moduleY; + unsigned int channel = static_cast(colID_original / groupingValue); + double colID = (static_cast(groupingValue) - 1) / 2 + channel * groupingValue; + + if (fiber == "S") { + SiPM_X = +moduleX / 2 - tuberadius - (tuberadius * 2) * colID; + SiPM_Y = -moduleY / 2 + tuberadius + (sq3 * tuberadius) * rowID + tuberadius * (2. * sq3m1 - 1.); + } + + if (fiber == "C") { + SiPM_X = +moduleX / 2 + tuberadius - (tuberadius * 2) * colID; + SiPM_Y = -moduleY / 2 + tuberadius + (sq3 * tuberadius) * rowID + tuberadius * (2. * sq3m1 - 1.); + } + + SiPM_Y = TowerOffsetY + SiPM_Y; +} + +std::string simTower_to_tbTower(std::string simTower) +{ + return std::to_string(306 + std::stoi(simTower)); +} + +GaussianFitSummary FitGaussianSummary(TH1* hist) +{ + GaussianFitSummary out; + if (!hist || hist->GetEntries() < 10) { + return out; + } + + hist->Fit("gaus", "Q0"); + TF1* fit = hist->GetFunction("gaus"); + if (!fit) { + return out; + } + + out.mean = fit->GetParameter(1); + out.meanErr = fit->GetParError(1); + out.rms = fit->GetParameter(2); + out.rmsErr = fit->GetParError(2); + return out; +} + +HistMomentSummary GetHistMomentSummary(const TH1* hist) +{ + HistMomentSummary out; + if (!hist || hist->GetEntries() < 1) { + return out; + } + + out.mean = hist->GetMean(); + out.meanErr = hist->GetMeanError(); + out.rms = hist->GetRMS(); + out.rmsErr = hist->GetRMSError(); + return out; +} + +std::string CsvEscape(const std::string& value) +{ + std::string escaped = "\""; + for (char c : value) { + if (c == '"') { + escaped += "\"\""; + } else { + escaped += c; + } + } + escaped += "\""; + return escaped; +} + +void AppendSummaryCsvLocked(const std::string& csvPath, + const std::string& inputFile, + double truthEnergy, + int nentries, + const GaussianFitSummary& sFit, + const GaussianFitSummary& cFit, + const GaussianFitSummary& combFit, + const GaussianFitSummary& combChiFit, + const HistMomentSummary& sciX, + const HistMomentSummary& sciY, + const HistMomentSummary& cerX, + const HistMomentSummary& cerY, + const GaussianFitSummary& sFitFersOn, + const GaussianFitSummary& cFitFersOn) +{ + const int fd = open(csvPath.c_str(), O_CREAT | O_WRONLY | O_APPEND, 0644); + if (fd == -1) { + throw std::runtime_error("Cannot open summary CSV: " + csvPath); + } + + if (flock(fd, LOCK_EX) != 0) { + close(fd); + throw std::runtime_error("Cannot lock summary CSV: " + csvPath); + } + + struct stat st; + if (fstat(fd, &st) != 0) { + flock(fd, LOCK_UN); + close(fd); + throw std::runtime_error("Cannot stat summary CSV: " + csvPath); + } + + std::ostringstream out; + out << std::setprecision(17); + + if (st.st_size == 0) { + out + << "input_file,truth_energy,nentries," + << "s_fit_mean,s_fit_mean_err,s_fit_rms,s_fit_rms_err," + << "c_fit_mean,c_fit_mean_err,c_fit_rms,c_fit_rms_err," + << "comb_fit_mean,comb_fit_mean_err,comb_fit_rms,comb_fit_rms_err," + << "combchi_fit_mean,combchi_fit_mean_err,combchi_fit_rms,combchi_fit_rms_err," + << "res_sci_x_mean,res_sci_x_mean_err,res_sci_x_rms,res_sci_x_rms_err," + << "res_sci_y_mean,res_sci_y_mean_err,res_sci_y_rms,res_sci_y_rms_err," + << "res_cer_x_mean,res_cer_x_mean_err,res_cer_x_rms,res_cer_x_rms_err," + << "res_cer_y_mean,res_cer_y_mean_err,res_cer_y_rms,res_cer_y_rms_err," + << "s_fit_fers_on_mean,s_fit_fers_on_mean_err,s_fit_fers_on_rms,s_fit_fers_on_rms_err," + << "c_fit_fers_on_mean,c_fit_fers_on_mean_err,c_fit_fers_on_rms,c_fit_fers_on_rms_err\n"; + } + + out + << CsvEscape(inputFile) << ',' + << truthEnergy << ',' + << nentries << ',' + << sFit.mean << ',' << sFit.meanErr << ',' << sFit.rms << ',' << sFit.rmsErr << ',' + << cFit.mean << ',' << cFit.meanErr << ',' << cFit.rms << ',' << cFit.rmsErr << ',' + << combFit.mean << ',' << combFit.meanErr << ',' << combFit.rms << ',' << combFit.rmsErr << ',' + << combChiFit.mean << ',' << combChiFit.meanErr << ',' << combChiFit.rms << ',' << combChiFit.rmsErr << ',' + << sciX.mean << ',' << sciX.meanErr << ',' << sciX.rms << ',' << sciX.rmsErr << ',' + << sciY.mean << ',' << sciY.meanErr << ',' << sciY.rms << ',' << sciY.rmsErr << ',' + << cerX.mean << ',' << cerX.meanErr << ',' << cerX.rms << ',' << cerX.rmsErr << ',' + << cerY.mean << ',' << cerY.meanErr << ',' << cerY.rms << ',' << cerY.rmsErr << ',' + << sFitFersOn.mean << ',' << sFitFersOn.meanErr << ',' << sFitFersOn.rms << ',' << sFitFersOn.rmsErr << ',' + << cFitFersOn.mean << ',' << cFitFersOn.meanErr << ',' << cFitFersOn.rms << ',' << cFitFersOn.rmsErr << '\n'; + + const std::string text = out.str(); + const ssize_t written = write(fd, text.c_str(), text.size()); + if (written != static_cast(text.size())) { + flock(fd, LOCK_UN); + close(fd); + throw std::runtime_error("Cannot write summary CSV: " + csvPath); + } + + flock(fd, LOCK_UN); + close(fd); +} + +void HidraTB25Ana(double energy, const std::string& input) +{ + const std::string sipmMapPath = + "/home/apareti/HidraSim2025/TBDataPreparation/2025_SPS/MapAndCalibration/sipm_map.json"; + + SipmLookup sipmLookup; + try { + sipmLookup = LoadSipmMap(sipmMapPath); + std::cout << "Loaded " << sipmLookup.size() << " SiPM map entries from " + << sipmMapPath << std::endl; + } catch (const std::exception& e) { + std::cerr << e.what() << std::endl; + return; + } + + const std::string infile = "../build/" + input; + std::cout << "Using file: " << infile << std::endl; + + TFile* simfile = TFile::Open(infile.c_str(), "READ"); + if (!simfile || simfile->IsZombie()) { + std::cerr << "Cannot open input file " << infile << std::endl; + return; + } + + TTree* simtree = static_cast(simfile->Get("DREMTubesout")); + if (!simtree) { + std::cerr << "Cannot find TTree DREMTubesout in " << infile << std::endl; + simfile->Close(); + delete simfile; + return; + } + + std::ostringstream os; + os << energy; + const std::string enstr = os.str(); + const std::string outfile = "hidra" + enstr + ".root"; + TFile f(outfile.c_str(), "RECREATE"); + + TDirectory* eventDisplayDir = f.mkdir("EventDisplays"); + f.cd(); + + int modcol[NofModulesX * NofModulesY]; + int modrow[NofModulesX * NofModulesY]; + for (int i = 0; i < NofModulesX * NofModulesY; i++) { + int row = i / NofModulesX; + int col = i % NofModulesX; + int imod = modflag[i]; + if (imod >= 0) { + modcol[imod] = col; + modrow[imod] = row; + } + } + + const double bmin = energy - 0.4 * std::sqrt(energy) * 10.; + const double bmax = energy + 0.4 * std::sqrt(energy) * 10.; + + auto sciene = new TH1F("sciene", "sciene", 100, bmin, bmax); + auto cerene = new TH1F("cerene", "cerene", 100, bmin, bmax); + auto totene = new TH1F("totene", "totene", 100, bmin, bmax); + auto totenec = new TH1F("totenec", "totenec", 100, bmin, bmax); + auto totdep = new TH1F("totdep", "totdep", 100, 0., bmax); + auto leakene = new TH1F("leakene", "leakene", 100, 0., 0.1); + auto chidist = new TH1F("chidist", "chidist", 100, 0., 1.); + auto mapcalo = new TH2F("mapcalo", "mapcalo", + NofModulesX, 0., NofModulesX, + NofModulesY, 0., NofModulesY); + auto SipmMapS = new TH2F("SipmMapS", "SipmS; Col; Row", + NofSiPMTowersX * NofFiberscolumn, 0, NofSiPMTowersX * NofFiberscolumn, + NofSiPMTowersY * NofFibersrow, 0, NofSiPMTowersY * NofFibersrow); + auto SipmMapC = new TH2F("SipmMapC", "SipmC; Col; Row", + NofSiPMTowersX * NofFiberscolumn, 0, NofSiPMTowersX * NofFiberscolumn, + NofSiPMTowersY * NofFibersrow, 0, NofSiPMTowersY * NofFibersrow); + + auto SciSiPMCoordinates = new TH2F("SciSiPMCoordinates", "Sci SiPM Coordinates; X [mm]; Y[mm]", + NofSiPMTowersX * NofFiberscolumn / grouping, + -NofSiPMTowersX * moduleX / 2, NofSiPMTowersX * moduleX / 2, + NofSiPMTowersY * NofFibersrow, + -NofSiPMTowersY * moduleY / 2, NofSiPMTowersY * moduleY / 2); + auto CerSiPMCoordinates = new TH2F("CerSiPMCoordinates", "Cer SiPM Coordinates; X [mm]; Y[mm]", + NofSiPMTowersX * NofFiberscolumn / grouping, + -NofSiPMTowersX * moduleX / 2, NofSiPMTowersX * moduleX / 2, + NofSiPMTowersY * NofFibersrow, + -NofSiPMTowersY * moduleY / 2, NofSiPMTowersY * moduleY / 2); + + auto ResidualHistSciX = new TH1F("ResidualHistSciX", "Residual Sci X; Residual [mm]; Entries", + 100, -10, 10); + auto ResidualHistSciY = new TH1F("ResidualHistSciY", "Residual Sci Y; Residual [mm]; Entries", + 100, -10, 10); + auto ResidualHistCerX = new TH1F("ResidualHistCerX", "Residual Cer X; Residual [mm]; Entries", + 100, -10, 10); + auto ResidualHistCerY = new TH1F("ResidualHistCerY", "Residual Cer Y; Residual [mm]; Entries", + 100, -10, 10); + + auto SignalSfersID = new TH1F("SignalSfersID", "Signal S fibers FERS ID; FERS; Signal", 16, 0, 16); + auto SignalCfersID = new TH1F("SignalCfersID", "Signal C fibers FERS ID; FERS; Signal", 16, 0, 16); + + auto scieneFersOn = new TH1F("scieneFersOn", "S energy for Activated FERS", 100, 0., bmax); + auto cereneFersOn = new TH1F("cereneFersOn", "C energy for Activated FERS", 100, 0., bmax); + + const int nentries = simtree->GetEntries(); + std::cout << "Entries " << nentries << std::endl; + + int pdg; + simtree->SetBranchAddress("PrimaryPDGID", &pdg); + + double venergy; + simtree->SetBranchAddress("PrimaryParticleEnergy", &venergy); + + double lenergy; + simtree->SetBranchAddress("EscapedEnergyl", &lenergy); + + double denergy; + simtree->SetBranchAddress("EscapedEnergyd", &denergy); + + double edep; + simtree->SetBranchAddress("EnergyTot", &edep); + + double Stot; + simtree->SetBranchAddress("NofPMTScinDet", &Stot); + + double Ctot; + simtree->SetBranchAddress("NofPMTCherDet", &Ctot); + + double PSdep; + simtree->SetBranchAddress("PSEnergy", &PSdep); + + double beamX; + simtree->SetBranchAddress("PrimaryX", &beamX); + + double beamY; + simtree->SetBranchAddress("PrimaryY", &beamY); + + std::vector* TowerE = nullptr; + simtree->SetBranchAddress("VecTowerE", &TowerE); + + std::vector* SPMT = nullptr; + simtree->SetBranchAddress("VecSPMT", &SPMT); + + std::vector* CPMT = nullptr; + simtree->SetBranchAddress("VecCPMT", &CPMT); + + std::vector* SSiPM = nullptr; + simtree->SetBranchAddress("VectorSignals", &SSiPM); + + std::vector* CSiPM = nullptr; + simtree->SetBranchAddress("VectorSignalsCher", &CSiPM); + + for (unsigned int i = 0; i < static_cast(nentries); i++) { + simtree->GetEntry(i); + + const bool writeEventDisplay = + (EventDisplayEvery > 0) && (((i + 1) % EventDisplayEvery) == 0); + + double ecalo = energy - lenergy / 1000.0; + double totsci = 0.; + double totcer = 0.; + double tottow = 0.; + + double barX_sci = 0.; + double barY_sci = 0.; + double barX_cer = 0.; + double barY_cer = 0.; + + std::unordered_map fersChannelsOverThreshold; + std::vector sciHitsPerFers; + std::vector cerHitsPerFers; + std::vector eventDisplaySciHits; + std::vector eventDisplayCerHits; + + sciHitsPerFers.reserve(SSiPM->size()); + cerHitsPerFers.reserve(CSiPM->size()); + + if (writeEventDisplay) { + eventDisplaySciHits.reserve(SSiPM->size()); + eventDisplayCerHits.reserve(CSiPM->size()); + } + + for (unsigned int j = 0; j < SPMT->size(); j++) { + totsci += SPMT->at(j) / sciPheGeV; + totcer += CPMT->at(j) / cerPheGeV; + tottow += TowerE->at(j); + mapcalo->Fill(modcol[j], modrow[j], TowerE->at(j) / 1000. / nentries); + } + + double sciPosWeight = 0.0; + double cerPosWeight = 0.0; + + for (unsigned int n = 0; n < SSiPM->size(); n++) { + const double rawPhe = SSiPM->at(n); + const double content = rawPhe / sciPheGeV; + + totsci += content; + + const unsigned int towID = static_cast(n / (NofFiberscolumn * NofFibersrow / 2)); + const unsigned int SiPMID = n % (NofFiberscolumn * NofFibersrow / 2); + const unsigned int colID = NofFiberscolumn - static_cast(SiPMID / (NofFibersrow / 2)); + const unsigned int rowID = 2 * static_cast(SiPMID % (NofFibersrow / 2)); + + double SiPM_X = 0.0; + double SiPM_Y = 0.0; + + GetSiPMcoordinate(towID, rowID, colID, SiPM_X, SiPM_Y, "S", grouping); + + const std::string tbTower = simTower_to_tbTower(std::to_string(towID)); + const int groupedCol = grouped_column(static_cast(colID), grouping); + + SipmMapS->Fill(colID, towID * NofFibersrow + rowID, content); + + const SipmMapEntry* info = + FindSipmInfo(sipmLookup, tbTower, "S", static_cast(rowID), groupedCol, grouping); + + if (info) { + SciSiPMCoordinates->Fill(info->x, info->y, content); + SignalSfersID->Fill(info->fersId, content / nentries); + barX_sci += info->x * content; + barY_sci += info->y * content; + sciPosWeight += content; + + if (info->boardID >= 0 && info->fersId >= 0) { + const FersKey key = make_fers_key(info->boardID, info->fersId); + const double threshold = GetFersThreshold(info->fersId, FersThresholdS); + + sciHitsPerFers.push_back({key, rawPhe, content}); + + if (content > threshold) { + ++fersChannelsOverThreshold[key]; + } + + if (writeEventDisplay) { + eventDisplaySciHits.push_back({key, info->x, info->y, rawPhe, content}); + } + } + } + } + + for (unsigned int n = 0; n < CSiPM->size(); n++) { + const double rawPhe = CSiPM->at(n); + const double content = rawPhe / cerPheGeV; + + totcer += content; + + const unsigned int towID = static_cast(n / (NofFiberscolumn * NofFibersrow / 2)); + const unsigned int SiPMID = n % (NofFiberscolumn * NofFibersrow / 2); + const unsigned int colID = NofFiberscolumn - static_cast(SiPMID / (NofFibersrow / 2)); + const unsigned int rowID = 2 * static_cast(SiPMID % (NofFibersrow / 2)) + 1; + + SipmMapC->Fill(colID, towID * NofFibersrow + rowID, content); + + const std::string tbTower = simTower_to_tbTower(std::to_string(towID)); + const int groupedCol = grouped_column(static_cast(colID), grouping); + + const SipmMapEntry* info = + FindSipmInfo(sipmLookup, tbTower, "C", static_cast(rowID), groupedCol, grouping); + + if (info) { + CerSiPMCoordinates->Fill(info->x, info->y, content); + SignalCfersID->Fill(info->fersId, content / nentries); + barX_cer += info->x * content; + barY_cer += info->y * content; + cerPosWeight += content; + + if (info->boardID >= 0 && info->fersId >= 0) { + const FersKey key = make_fers_key(info->boardID, info->fersId); + const double threshold = GetFersThreshold(info->fersId, FersThresholdC); + + cerHitsPerFers.push_back({key, rawPhe, content}); + + if (content > threshold) { + ++fersChannelsOverThreshold[key]; + } + + if (writeEventDisplay) { + eventDisplayCerHits.push_back({key, info->x, info->y, rawPhe, content}); + } + } + } + } + + double totsciFersOn = 0.; + for (const auto& hit : sciHitsPerFers) { + if (IsActivatedFers(fersChannelsOverThreshold, hit.key) && + hit.rawPhe >= MinVisiblePheS) { + if(ApplySaturation){ + totsciFersOn += apply_sipm_saturation_from_pe(hit.rawPhe, sciPheGeV, NofSipmCells_sci); + } + else{ + totsciFersOn += hit.signal; + } + } + } + + double totcerFersOn = 0.; + for (const auto& hit : cerHitsPerFers) { + if (IsActivatedFers(fersChannelsOverThreshold, hit.key) && + hit.rawPhe >= MinVisiblePheC ) { + if(ApplySaturation){ + totcerFersOn += apply_sipm_saturation_from_pe(hit.rawPhe, cerPheGeV, NofSipmCells_cer); + } + else{ + totcerFersOn += hit.signal; + } + } + } + + if (sciPosWeight > 0.) { + ResidualHistSciX->Fill((barX_sci / sciPosWeight) - beamX); + ResidualHistSciY->Fill((barY_sci / sciPosWeight) - beamY); + } + + if (cerPosWeight > 0.) { + ResidualHistCerX->Fill((barX_cer / cerPosWeight) - beamX); + ResidualHistCerY->Fill((barY_cer / cerPosWeight) - beamY); + } + + if (writeEventDisplay && eventDisplayDir) { + eventDisplayDir->cd(); + + std::ostringstream sciTitle; + sciTitle << "Sci event display, entry " << i + << " (grouped channels, activated FERS only); X [mm]; Y [mm]"; + TH2F* sciEventDisplay = CreateEventDisplayHist( + make_event_display_name("EventDisplaySci", i), + sciTitle.str(), + grouping + ); + + std::ostringstream cerTitle; + cerTitle << "Cer event display, entry " << i + << " (grouped channels, activated FERS only); X [mm]; Y [mm]"; + TH2F* cerEventDisplay = CreateEventDisplayHist( + make_event_display_name("EventDisplayCer", i), + cerTitle.str(), + grouping + ); + + FillActivatedEventDisplay( + sciEventDisplay, + eventDisplaySciHits, + fersChannelsOverThreshold, + MinVisiblePheS + ); + FillActivatedEventDisplay( + cerEventDisplay, + eventDisplayCerHits, + fersChannelsOverThreshold, + MinVisiblePheC + ); + + sciEventDisplay->Write(); + cerEventDisplay->Write(); + + delete sciEventDisplay; + delete cerEventDisplay; + + f.cd(); + } + + scieneFersOn->Fill(totsciFersOn); + cereneFersOn->Fill(totcerFersOn); + + sciene->Fill(totsci); + cerene->Fill(totcer); + //totene->Fill(elcont * 0.5 * (totsci + totcer)); + totene->Fill(0.5 * (totsci + totcer)); + totenec->Fill(picont * (totsci - chi * totcer) / (1 - chi)); + totdep->Fill(tottow / 1000.); + leakene->Fill(lenergy / 1000. / energy); + chidist->Fill((totsci - ecalo) / (totcer - ecalo)); + } + + const GaussianFitSummary sFit = FitGaussianSummary(sciene); + const GaussianFitSummary cFit = FitGaussianSummary(cerene); + const GaussianFitSummary combFit = FitGaussianSummary(totene); + const GaussianFitSummary combChiFit = FitGaussianSummary(totenec); + + const HistMomentSummary resSciX = GetHistMomentSummary(ResidualHistSciX); + const HistMomentSummary resSciY = GetHistMomentSummary(ResidualHistSciY); + const HistMomentSummary resCerX = GetHistMomentSummary(ResidualHistCerX); + const HistMomentSummary resCerY = GetHistMomentSummary(ResidualHistCerY); + + const GaussianFitSummary sFitFersOn = FitGaussianSummary(scieneFersOn); + const GaussianFitSummary cFitFersOn = FitGaussianSummary(cereneFersOn); + + try { + AppendSummaryCsvLocked("hidra_summary.csv", + input, + energy, + nentries, + sFit, + cFit, + combFit, + combChiFit, + resSciX, + resSciY, + resCerX, + resCerY, + sFitFersOn, + cFitFersOn); + } catch (const std::exception& e) { + std::cerr << e.what() << std::endl; + } + + f.Write(); + f.Close(); + + simfile->Close(); + delete simfile; +} \ No newline at end of file From 7b3b27e2d165b7f218bfe6aa7496ab612d5c7db3 Mon Sep 17 00:00:00 2001 From: AndreaPareti Date: Mon, 13 Apr 2026 16:20:34 +0200 Subject: [PATCH 2/5] Add plotting functions --- analysis/plot_hidra.py | 484 +++++++++++++++++++++++++++++++++++++++++ analysis/run_ana.sh | 41 ++++ 2 files changed, 525 insertions(+) create mode 100644 analysis/plot_hidra.py create mode 100755 analysis/run_ana.sh diff --git a/analysis/plot_hidra.py b/analysis/plot_hidra.py new file mode 100644 index 0000000..cb0f3ad --- /dev/null +++ b/analysis/plot_hidra.py @@ -0,0 +1,484 @@ +# file: plot_hidra_summary.py +import os +import sys +from array import array + +import numpy as np +import pandas as pd +import ROOT + + +ROOT.gROOT.SetBatch(True) +ROOT.gStyle.SetOptStat(0) + +#scaling_s = 1.16 +#scaling_c = 1.255 +scaling_s = 1.125 +scaling_c = 1.155 +myOutDir = "plots/Attenuation7m_2phe" +extraTex2 = "Fibre attenuation set to 10m" +#myOutDir = "plots/AttenuationTB24_2phe" +extraTex1 = "Min visible signal: >= 2 phe" +#extraTex2 = "Fibre attenuation tuned to TB24 data: 3.7 m (Sci), 3.9 m (Cer)" +extraTex3 = "Correction for multiple photons on the same SiPM cell: OFF" + +def load_summary(csv_path: str) -> pd.DataFrame: + data = pd.read_csv(csv_path) + data = data.sort_values("truth_energy").reset_index(drop=True) + return data + + +def safe_ratio(num: np.ndarray, den: np.ndarray) -> np.ndarray: + out = np.full_like(num, np.nan, dtype=float) + mask = np.isfinite(num) & np.isfinite(den) & (den != 0.0) + out[mask] = num[mask] / den[mask] + return out + + +def safe_ratio_error( + num: np.ndarray, + num_err: np.ndarray, + den: np.ndarray, + den_err: np.ndarray, +) -> np.ndarray: + out = np.full_like(num, np.nan, dtype=float) + mask = ( + np.isfinite(num) + & np.isfinite(num_err) + & np.isfinite(den) + & np.isfinite(den_err) + & (num != 0.0) + & (den != 0.0) + ) + ratio = np.full_like(num, np.nan, dtype=float) + ratio[mask] = num[mask] / den[mask] + rel = np.full_like(num, np.nan, dtype=float) + rel[mask] = np.sqrt((num_err[mask] / num[mask]) ** 2 + (den_err[mask] / den[mask]) ** 2) + out[mask] = ratio[mask] * rel[mask] + return out + + +def to_root_arrays(x, y, ex, ey): + return ( + array("d", [float(v) for v in x]), + array("d", [float(v) for v in y]), + array("d", [float(v) for v in ex]), + array("d", [float(v) for v in ey]), + ) + + +def make_graph( + x: np.ndarray, + y: np.ndarray, + ey: np.ndarray, + color: int, + marker: int, + name: str, + line_style: int = 1, +): + ex = np.zeros_like(x, dtype=float) + mask = np.isfinite(x) & np.isfinite(y) & np.isfinite(ey) + x_use = x[mask] + y_use = y[mask] + ex_use = ex[mask] + ey_use = ey[mask] + + if len(x_use) == 0: + return None + + x_arr, y_arr, ex_arr, ey_arr = to_root_arrays(x_use, y_use, ex_use, ey_use) + graph = ROOT.TGraphErrors(len(x_use), x_arr, y_arr, ex_arr, ey_arr) + graph.SetName(name) + graph.SetMarkerColor(color) + graph.SetLineColor(color) + graph.SetMarkerStyle(marker) + graph.SetMarkerSize(1.2) + graph.SetLineWidth(2) + graph.SetLineStyle(line_style) + return graph + + +def fit_resolution_graph(graph: ROOT.TGraphErrors, color: int, fit_name: str): + if graph is None or graph.GetN() < 2: + return None, None + + fit = ROOT.TF1(fit_name, "pol1", 0.0, 1.0) + fit.SetLineColor(color) + fit.SetLineWidth(2) + result = graph.Fit(fit, "QS") + if int(result) != 0: + return fit, None + return fit, result + + +def draw_header(x_left: float, y_top: float, subtitle: str): + tex0 = ROOT.TLatex(x_left, y_top, "HidraSim") + tex0.SetNDC() + tex0.SetTextFont(72) + tex0.SetTextSize(0.042) + tex0.Draw() + + tex1 = ROOT.TLatex(x_left + 0.16, y_top, "Preliminary") + tex1.SetNDC() + tex1.SetTextFont(42) + tex1.SetTextSize(0.042) + tex1.Draw() + + tex2 = ROOT.TLatex(x_left, y_top - 0.055, subtitle) + tex2.SetNDC() + tex2.SetTextFont(42) + tex2.SetTextSize(0.032) + tex2.Draw() + + +def setup_canvas(name: str, width: int = 800, height: int = 650): + canvas = ROOT.TCanvas(name, name, width, height) + canvas.SetFillColor(0) + canvas.SetBorderMode(0) + canvas.SetBorderSize(2) + canvas.SetTickx(1) + canvas.SetTicky(1) + canvas.SetLeftMargin(0.13) + canvas.SetRightMargin(0.05) + canvas.SetTopMargin(0.09) + canvas.SetBottomMargin(0.11) + canvas.SetFrameBorderMode(0) + return canvas + + +def plot_energy_resolution(data: pd.DataFrame, out_dir: str): + energy = data["truth_energy"].to_numpy(dtype=float) + inv_sqrt_e = 1.0 / np.sqrt(energy) + + s_res = safe_ratio( + data["s_fit_rms"].to_numpy(dtype=float), + data["s_fit_mean"].to_numpy(dtype=float), + ) + s_res_err = safe_ratio_error( + data["s_fit_rms"].to_numpy(dtype=float), + data["s_fit_rms_err"].to_numpy(dtype=float), + data["s_fit_mean"].to_numpy(dtype=float), + data["s_fit_mean_err"].to_numpy(dtype=float), + ) + + c_res = safe_ratio( + data["c_fit_rms"].to_numpy(dtype=float), + data["c_fit_mean"].to_numpy(dtype=float), + ) + c_res_err = safe_ratio_error( + data["c_fit_rms"].to_numpy(dtype=float), + data["c_fit_rms_err"].to_numpy(dtype=float), + data["c_fit_mean"].to_numpy(dtype=float), + data["c_fit_mean_err"].to_numpy(dtype=float), + ) + + comb_res = safe_ratio( + data["comb_fit_rms"].to_numpy(dtype=float), + data["comb_fit_mean"].to_numpy(dtype=float), + ) + comb_res_err = safe_ratio_error( + data["comb_fit_rms"].to_numpy(dtype=float), + data["comb_fit_rms_err"].to_numpy(dtype=float), + data["comb_fit_mean"].to_numpy(dtype=float), + data["comb_fit_mean_err"].to_numpy(dtype=float), + ) + + canvas = setup_canvas("c_resolution") + mg = ROOT.TMultiGraph() + mg.SetTitle("Energy resolution;1/#sqrt{E_{beam}} [GeV^{-1/2}];#sigma(E)/E") + + graph_comb = make_graph(inv_sqrt_e, comb_res, comb_res_err, ROOT.kGreen + 2, 20, "graph_comb") + graph_s = make_graph(inv_sqrt_e, s_res, s_res_err, ROOT.kRed + 1, 21, "graph_s") + graph_c = make_graph(inv_sqrt_e, c_res, c_res_err, ROOT.kBlue + 1, 22, "graph_c") + + for graph in (graph_comb, graph_s, graph_c): + if graph is not None: + mg.Add(graph, "PE") + + mg.Draw("A") + mg.GetXaxis().SetTitleOffset(1.15) + mg.GetYaxis().SetTitleOffset(1.35) + + fit_comb, _ = fit_resolution_graph(graph_comb, ROOT.kGreen + 2, "fit_comb") + fit_s, _ = fit_resolution_graph(graph_s, ROOT.kRed + 1, "fit_s") + fit_c, _ = fit_resolution_graph(graph_c, ROOT.kBlue + 1, "fit_c") + + + alignLeft = 0.18 + alignTop = 0.85 + alignRight = 0.925 + tex0 = ROOT.TLatex(alignLeft, alignTop, "HidraSim") + tex0.SetNDC() + tex0.SetTextFont(72); + tex0.SetTextSize(0.042); + tex0.SetLineWidth(2); + tex0.Draw("same"); + #tex1 = ROOT.TLatex(alignLeft+0.16,alignTop,"Preliminary "); + tex1 = ROOT.TLatex(alignLeft+0.16,alignTop,"Work in progress"); + tex1.SetNDC(); + tex1.SetTextFont(42); + tex1.SetTextSize(0.042); + tex1.SetLineWidth(2); + tex1.Draw("same"); + + legend = ROOT.TLegend(0.16, 0.57, 0.72, 0.82) + legend.SetBorderSize(0) + legend.SetFillStyle(0) + legend.SetTextFont(42) + legend.SetTextSize(0.03) + + if fit_comb is not None: + legend.AddEntry( + fit_comb, + "Combined: #sigma/E = %.2f%%/#sqrt{E} %+ .2f%%" + % (100.0 * fit_comb.GetParameter(1), 100.0 * fit_comb.GetParameter(0)), + "l", + ) + if fit_s is not None: + legend.AddEntry( + fit_s, + "S: #sigma/E = %.2f%%/#sqrt{E} %+ .2f%%" + % (100.0 * fit_s.GetParameter(1), 100.0 * fit_s.GetParameter(0)), + "l", + ) + if fit_c is not None: + legend.AddEntry( + fit_c, + "C: #sigma/E = %.2f%%/#sqrt{E} %+ .2f%%" + % (100.0 * fit_c.GetParameter(1), 100.0 * fit_c.GetParameter(0)), + "l", + ) + + legend.Draw() + draw_header(0.16, 0.84, "Energy resolution from Gaussian fits") + canvas.SaveAs(os.path.join(out_dir, "energy_resolution.pdf")) + + +def plot_energy_linearity(data: pd.DataFrame, out_dir: str): + energy = data["truth_energy"].to_numpy(dtype=float) + + s_mean = data["s_fit_mean"].to_numpy(dtype=float) + s_mean_err = data["s_fit_mean_err"].to_numpy(dtype=float) + c_mean = data["c_fit_mean"].to_numpy(dtype=float) + c_mean_err = data["c_fit_mean_err"].to_numpy(dtype=float) + comb_mean = data["comb_fit_mean"].to_numpy(dtype=float) + comb_mean_err = data["comb_fit_mean_err"].to_numpy(dtype=float) + + s_fers_on_mean = data["s_fit_fers_on_mean"].to_numpy(dtype=float) + s_fers_on_mean_err = data["s_fit_fers_on_mean_err"].to_numpy(dtype=float) + c_fers_on_mean = data["c_fit_fers_on_mean"].to_numpy(dtype=float) + c_fers_on_mean_err = data["c_fit_fers_on_mean_err"].to_numpy(dtype=float) + + #s_lin = safe_ratio(s_mean - energy, energy) + s_lin = safe_ratio(s_mean, energy) + s_lin_err = safe_ratio(s_mean_err, energy) + + #c_lin = safe_ratio(c_mean - energy, energy) + c_lin = safe_ratio(c_mean, energy) + c_lin_err = safe_ratio(c_mean_err, energy) + + #comb_lin = safe_ratio(comb_mean - energy, energy) + comb_lin = safe_ratio(comb_mean, energy) + comb_lin_err = safe_ratio(comb_mean_err, energy) + + #s_fers_on_lin = safe_ratio(s_fers_on_mean - energy, energy) + s_fers_on_lin = safe_ratio(s_fers_on_mean*scaling_s, energy) + s_fers_on_lin_err = safe_ratio(s_fers_on_mean_err, energy) + + #c_fers_on_lin = safe_ratio(c_fers_on_mean - energy, energy) + c_fers_on_lin = safe_ratio(c_fers_on_mean*scaling_c, energy) + c_fers_on_lin_err = safe_ratio(c_fers_on_mean_err, energy) + + canvas = setup_canvas("c_linearity") + mg = ROOT.TMultiGraph() + #mg.SetTitle("Energy linearity;E_{beam} [GeV];(E_{fit}-E_{beam})/E_{beam}") + mg.SetTitle("Energy linearity;E_{beam} [GeV];E_{fit} / E_{beam}") + + graph_comb = make_graph(energy, comb_lin, comb_lin_err, ROOT.kGreen + 2, 20, "lin_comb") + graph_s = make_graph(energy, s_lin, s_lin_err, ROOT.kRed + 1, 21, "lin_s") + graph_c = make_graph(energy, c_lin, c_lin_err, ROOT.kBlue + 1, 22, "lin_c") + graph_s_fers_on = make_graph( + energy, + s_fers_on_lin, + s_fers_on_lin_err, + ROOT.kRed + 1, + 25, + "lin_s_fers_on", + line_style=1, + ) + graph_c_fers_on = make_graph( + energy, + c_fers_on_lin, + c_fers_on_lin_err, + ROOT.kBlue + 1, + 26, + "lin_c_fers_on", + line_style=1, + ) + + #for graph in (graph_comb, graph_s, graph_c, graph_s_fers_on, graph_c_fers_on): + for graph in (graph_comb, graph_s, graph_c): + if graph is not None: + mg.Add(graph, "PE") + + mg.Add(graph_s_fers_on, "PEL") + mg.Add(graph_c_fers_on, "PEL") + + mg.Draw("A") + mg.GetXaxis().SetTitleOffset(1.07) + mg.GetYaxis().SetTitleOffset(1.2) + #mg.SetMinimum(-0.40) + #mg.SetMaximum(0.08) + + mg.SetMinimum(+0.9) + mg.SetMaximum(+1.14) + + xmin = float(np.nanmin(energy)) + xmax = float(np.nanmax(energy)) + #for y_val, color in ((0.0, ROOT.kGreen + 3), (0.01, ROOT.kGray + 2), (-0.01, ROOT.kGray + 2)): + # line = ROOT.TLine(xmin, y_val, xmax, y_val) + # line.SetLineStyle(2) + # line.SetLineColor(color) + # line.Draw() + + + line1=ROOT.TLine(mg.GetXaxis().GetXmin(), 1.01, mg.GetXaxis().GetXmax(), 1.01) + line2=ROOT.TLine(mg.GetXaxis().GetXmin(), 0.99, mg.GetXaxis().GetXmax(), 0.99) + line3=ROOT.TLine(mg.GetXaxis().GetXmin(), 1.0, mg.GetXaxis().GetXmax(), 1.0) + line1.SetLineStyle(2) + line2.SetLineStyle(2) + line3.SetLineStyle(2) + line3.SetLineColor(ROOT.kGreen+3) + + line1.Draw("same") + line2.Draw("same") + line3.Draw("same") + line2.Draw("same") + line3.Draw("same") + + alignLeft = 0.18 + alignTop = 0.85 + alignRight = 0.925 + tex0 = ROOT.TLatex(alignLeft, alignTop, "HidraSim") + tex0.SetNDC() + tex0.SetTextFont(72); + tex0.SetTextSize(0.042); + tex0.SetLineWidth(2); + tex0.Draw("same"); + #tex1 = ROOT.TLatex(alignLeft+0.16,alignTop,"Preliminary "); + tex1 = ROOT.TLatex(alignLeft+0.16,alignTop,"Work in progress"); + tex1.SetNDC(); + tex1.SetTextFont(42); + tex1.SetTextSize(0.042); + tex1.SetLineWidth(2); + tex1.Draw("same"); + tex2 = ROOT.TLatex(alignLeft,alignTop-0.04, extraTex2) + tex2.SetNDC(); + tex2.SetTextFont(42); + tex2.SetTextSize(0.03); + tex2.SetLineWidth(2); + tex2.Draw("same"); + tex3 = ROOT.TLatex(alignLeft,alignTop-0.08, extraTex1) + tex3.SetNDC(); + tex3.SetTextFont(42); + tex3.SetTextSize(0.03); + tex3.SetLineWidth(2); + tex3.Draw("same"); + + + legend = ROOT.TLegend(0.450, 0.17, 0.83, 0.37) + legend.SetBorderSize(0) + legend.SetFillStyle(0) + legend.SetTextFont(42) + legend.SetTextSize(0.03) + legend.SetHeader("Fers thresholds tuned to TB data") # Set the header for the legend + if graph_comb is not None: + legend.AddEntry(graph_comb, "Combined - all fibres", "pe") + if graph_s is not None: + legend.AddEntry(graph_s, "S - all fibres", "pe") + if graph_c is not None: + legend.AddEntry(graph_c, "C - all fibres", "pe") + if graph_s_fers_on is not None: + #legend.AddEntry(graph_s_fers_on, "S - FERS On only, Thr: >= 1 phe #times 1.1", "pe") + legend.AddEntry(graph_s_fers_on, "S - FERS On only #times " + str(scaling_s), "ple") + if graph_c_fers_on is not None: + #legend.AddEntry(graph_c_fers_on, "C - FERS On only, Thr: >= 1 phe #times 1.1", "pe") + legend.AddEntry(graph_c_fers_on, "C - FERS On only #times " + str(scaling_c), "ple") + legend.Draw() + + draw_header(0.16, 0.84, "Linearity from fitted energy means") + canvas.SaveAs(os.path.join(out_dir, "energy_linearity.pdf")) + + +def plot_spatial_resolution(data: pd.DataFrame, out_dir: str): + energy = data["truth_energy"].to_numpy(dtype=float) + + configs = [ + ( + "spatial_resolution_x.pdf", + "Spatial resolution X;E_{beam} [GeV];Residual RMS [mm]", + data["res_sci_x_rms"].to_numpy(dtype=float), + data["res_sci_x_rms_err"].to_numpy(dtype=float), + data["res_cer_x_rms"].to_numpy(dtype=float), + data["res_cer_x_rms_err"].to_numpy(dtype=float), + "Residual RMS along X", + ), + ( + "spatial_resolution_y.pdf", + "Spatial resolution Y;E_{beam} [GeV];Residual RMS [mm]", + data["res_sci_y_rms"].to_numpy(dtype=float), + data["res_sci_y_rms_err"].to_numpy(dtype=float), + data["res_cer_y_rms"].to_numpy(dtype=float), + data["res_cer_y_rms_err"].to_numpy(dtype=float), + "Residual RMS along Y", + ), + ] + + for file_name, title, sci_rms, sci_rms_err, cer_rms, cer_rms_err, subtitle in configs: + canvas = setup_canvas(file_name.replace(".pdf", "")) + mg = ROOT.TMultiGraph() + mg.SetTitle(title) + + graph_sci = make_graph(energy, sci_rms, sci_rms_err, ROOT.kRed + 1, 21, f"sci_{file_name}") + graph_cer = make_graph(energy, cer_rms, cer_rms_err, ROOT.kBlue + 1, 22, f"cer_{file_name}") + + if graph_sci is not None: + mg.Add(graph_sci, "PE") + if graph_cer is not None: + mg.Add(graph_cer, "PE") + + mg.Draw("A") + mg.GetXaxis().SetTitleOffset(1.15) + mg.GetYaxis().SetTitleOffset(1.35) + + legend = ROOT.TLegend(0.16, 0.70, 0.34, 0.82) + legend.SetBorderSize(0) + legend.SetTextFont(42) + legend.SetTextSize(0.03) + if graph_sci is not None: + legend.AddEntry(graph_sci, "Sci", "pe") + if graph_cer is not None: + legend.AddEntry(graph_cer, "Cer", "pe") + legend.Draw() + + draw_header(0.16, 0.84, subtitle) + canvas.SaveAs(os.path.join(out_dir, file_name)) + + +def main(): + csv_path = sys.argv[1] if len(sys.argv) > 1 else "hidra_summary.csv" + out_dir = sys.argv[2] if len(sys.argv) > 2 else myOutDir + + os.makedirs(out_dir, exist_ok=True) + + data = load_summary(csv_path) + plot_energy_resolution(data, out_dir) + plot_energy_linearity(data, out_dir) + plot_spatial_resolution(data, out_dir) + + print(f"Read: {csv_path}") + print(f"Wrote plots to: {out_dir}") + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/analysis/run_ana.sh b/analysis/run_ana.sh new file mode 100755 index 0000000..eacb4f3 --- /dev/null +++ b/analysis/run_ana.sh @@ -0,0 +1,41 @@ +# file: run_hidra_parallel.sh +#!/usr/bin/env bash +set -euo pipefail + +energies=(10 20 30 40 60 80 100 120) +files=( + "OutputElectrons_26_04_08/DREMTubesout_Run0.root" + "OutputElectrons_26_04_08/DREMTubesout_Run1.root" + "OutputElectrons_26_04_08/DREMTubesout_Run2.root" + "OutputElectrons_26_04_08/DREMTubesout_Run3.root" + "OutputElectrons_26_04_08/DREMTubesout_Run4.root" + "OutputElectrons_26_04_08/DREMTubesout_Run5.root" + "OutputElectrons_26_04_08/DREMTubesout_Run6.root" + "OutputElectrons_26_04_08/DREMTubesout_Run7.root" +) + +summary_csv="hidra_summary.csv" +pids=() + +rm -f "${summary_csv}" + +for i in "${!energies[@]}"; do + energy="${energies[$i]}" + input_file="${files[$i]}" + log_file="run_${i}_${energy}GeV.log" + + echo "Starting ${input_file} at ${energy} GeV" + root -l -b -q "HidraTB25Ana.C(${energy}, \"${input_file}\")" >"${log_file}" 2>&1 & + pids+=("$!") +done + +status=0 + +for pid in "${pids[@]}"; do + if ! wait "${pid}"; then + status=1 + fi +done + +echo "Done. Shared summary written to: ${summary_csv}" +exit "${status}" \ No newline at end of file From 629c3fef06cfa5be95bf634deec7f29d0c65defe Mon Sep 17 00:00:00 2001 From: AndreaPareti Date: Mon, 20 Apr 2026 12:36:06 +0200 Subject: [PATCH 3/5] Update on Sim comparison: include noise-like effect, phe losses ecc --- analysis/HidraTB25Ana.C | 573 +++++++++++++++++++++++++++++++--------- analysis/plot_hidra.py | 43 ++- 2 files changed, 484 insertions(+), 132 deletions(-) diff --git a/analysis/HidraTB25Ana.C b/analysis/HidraTB25Ana.C index 3f80140..8495dd1 100644 --- a/analysis/HidraTB25Ana.C +++ b/analysis/HidraTB25Ana.C @@ -1,11 +1,12 @@ -// file: HidraAna.C +// file: HidraTB25Ana.C #include #include #include +#include #include #include - +#include #include #include #include @@ -26,32 +27,30 @@ #include #include "HidraGeo.h" #include +#include +#include using json = nlohmann::json; // Global Constants and parameters // -// Fallback per-channel thresholds used only if a FERS id is missing in fers_to_thr_map -const double FersThresholdS = 0.07; // GeV -const double FersThresholdC = 0.07; // GeV - // Minimum visible signal per channel (>=) // Set to 0.0 or 1.0 to disable the cut. -const double MinVisiblePheS = 2.0; // photoelectrons -const double MinVisiblePheC = 2.0; // photoelectrons +const double MinVisiblePheS = 0.0; // photoelectrons +const double MinVisiblePheC = 0.0; // photoelectrons // Scintillator on even rows (start from 0), Cerenkov on odd rows const unsigned int grouping = 8; // Write event display every N entries (0 to disable) const unsigned int EventDisplayEvery = 100; -// A FERS is activated if at least this many channels exceed threshold +// A FERS is ON if at least 2 smeared channels across S or C (both are good) are above that FERS threshold. const unsigned int FersMultiplicity = 2; const double chi = 0.38; -//const double sciPheGeV = 119.001; // tb24 -//const double cerPheGeV = 29.4; // tb24 -const double sciPheGeV = 178.501; -const double cerPheGeV = 43; +const double sciPheGeV = 119.001; // tb24 +const double cerPheGeV = 29.4; // tb24 +//const double sciPheGeV = 178.501; +//const double cerPheGeV = 43; const double NofSipmCells_sci = 7772; const double NofSipmCells_cer = 3443; @@ -61,10 +60,27 @@ const double cer_pde = 0.38; const double elcont = 1.005; const double picont = 1.028; +// Separate fallback noise sigmas for S and C +const unsigned int NoiseRandomSeed = 12345; // fixed seed for reproducibility + // In this parametrised simulation, rawPhe is stored as "impinging optical photons" -// To account for multiple photons impinging on the same SiPM cell, +// To account for multiple photons impinging on the same SiPM cell, // we apply saturation correction to get "fired cells", which is the actual observable. const bool ApplySaturation = true; +// +const bool ApplyTbNoise = true; +// decide whether to print per-event FERS channel counts (for debugging) +bool printFersLog = false; + + + +// lambda = mean number of loss opportunities per channel. +// Probability to lose 1 pe = 1 - exp(-lambda). +const double SciOnePeLossLambda = 3; // set > 0 to enable +const double CerOnePeLossLambda = 3; // set > 0 to enable + + + double apply_sipm_saturation_from_pe(double pe, double phePerGeV, @@ -78,8 +94,6 @@ double apply_sipm_saturation_from_pe(double pe, return firedCells / phePerGeV; } - - std::map fers_to_thr_map = { { 1, 0.17 }, { 2, 0.20 }, @@ -99,6 +113,61 @@ std::map fers_to_thr_map = { {16, 0.06 } }; +using ChannelNoiseVector = std::vector; + +std::string Trim(const std::string& value) +{ + const std::size_t first = value.find_first_not_of(" \t\r\n"); + if (first == std::string::npos) { + return ""; + } + + const std::size_t last = value.find_last_not_of(" \t\r\n"); + return value.substr(first, last - first + 1); +} + +ChannelNoiseVector LoadNoiseValuesFromCsv(const std::string& csvPath) +{ + std::ifstream in(csvPath); + if (!in) { + throw std::runtime_error("Cannot open noise CSV: " + csvPath); + } + + ChannelNoiseVector values; + std::string line; + + while (std::getline(in, line)) { + const std::string trimmed = Trim(line); + if (trimmed.empty()) { + continue; + } + + const std::size_t commaPos = trimmed.find(','); + const std::string token = + Trim(trimmed.substr(0, commaPos == std::string::npos ? trimmed.size() : commaPos)); + + values.push_back(std::stod(token)); + } + + if (values.empty()) { + throw std::runtime_error("Noise CSV is empty: " + csvPath); + } + + return values; +} + +double GetChannelNoiseSigma(const ChannelNoiseVector& noiseValues, int ch) +{ + if (ch < 0 || static_cast(ch) >= noiseValues.size()) { + std::ostringstream os; + os << "Missing per-channel noise sigma for channel " << ch + << " (CSV size = " << noiseValues.size() << ")"; + throw std::runtime_error(os.str()); + } + + return noiseValues[static_cast(ch)]; +} + struct SipmMapEntry { int boardID = -1; std::string type; @@ -144,6 +213,121 @@ struct HistMomentSummary { using SipmLookup = std::unordered_map; using FersKey = uint64_t; +struct ChannelKey { + FersKey fersKey = 0; + int ch = -1; + + bool operator==(const ChannelKey& other) const + { + return fersKey == other.fersKey && ch == other.ch; + } +}; + +struct ChannelKeyHash { + std::size_t operator()(const ChannelKey& key) const noexcept + { + const std::size_t h1 = std::hash{}(key.fersKey); + const std::size_t h2 = std::hash{}(key.ch); + return h1 ^ (h2 + 0x9e3779b97f4a7c15ULL + (h1 << 6U) + (h1 >> 2U)); + } +}; + +struct ChannelAccumulatedSignal { + FersKey fersKey = 0; + int fersId = -1; + int ch = -1; + double threshold = 0.0; + double thresholdSignal = 0.0; + double outputSignal = 0.0; + double smearedThresholdSignal = 0.0; + double smearedOutputSignal = 0.0; +}; + +using ChannelSignalMap = + std::unordered_map; + +void SmearAndCountChannelsOverThreshold( + ChannelSignalMap& channelSignals, + std::unordered_map& fersChannelsOverThreshold, + TRandom3& rng, + const ChannelNoiseVector& channelNoiseValues) +{ + for (auto& [channelKey, channel] : channelSignals) { + const double sigma = GetChannelNoiseSigma(channelNoiseValues, channel.ch); + + double noise = 0.; + if(ApplyTbNoise){ + noise = (sigma > 0.0) ? rng.Gaus(0.0, sigma) : 0.0; + } + + + channel.smearedThresholdSignal = channel.thresholdSignal + noise; + channel.smearedOutputSignal = channel.outputSignal + noise; + + // FERS activation logic: if smeared threshold signal is above threshold, + // count this channel for FERS multiplicity + if (channel.smearedThresholdSignal > channel.threshold) { + ++fersChannelsOverThreshold[channel.fersKey]; + } + } +} + + + + + + + + +int board_id_from_fers_key(FersKey key) +{ + return static_cast(static_cast(key >> 32)); +} + +void PrintPerEventFersChannelCounts( + unsigned int eventIndex, + const ChannelSignalMap& sciChannelSignals, + const ChannelSignalMap& cerChannelSignals, + const std::unordered_map& fersChannelsOverThreshold) +{ + std::map, FersKey> fersSeen; + + auto collectFers = [&](const ChannelSignalMap& channelSignals) { + for (const auto& [channelKey, channel] : channelSignals) { + const int boardID = board_id_from_fers_key(channel.fersKey); + fersSeen[{boardID, channel.fersId}] = channel.fersKey; + } + }; + + collectFers(sciChannelSignals); + collectFers(cerChannelSignals); + + if (printFersLog) { + std::cout << "Event " << eventIndex << ":\n"; + for (const auto& [id, key] : fersSeen) { + const int boardID = id.first; + const int fersId = id.second; + + const auto it = fersChannelsOverThreshold.find(key); + const unsigned int nChannelsOverThreshold = + (it != fersChannelsOverThreshold.end()) ? it->second : 0U; + + std::cout + << " board " << boardID + << ", FERS " << fersId + << " -> channels above threshold = " << nChannelsOverThreshold; + + if (nChannelsOverThreshold >= FersMultiplicity) { + std::cout << " [ON]"; + } else { + std::cout << " [OFF]"; + } + + std::cout << '\n'; + } + } +} + std::string make_sipm_key(const std::string& tower, const std::string& type, int row, @@ -170,16 +354,128 @@ int grouped_to_json_column(int groupedCol, unsigned int groupingValue) return groupedCol * static_cast(groupingValue) + static_cast(groupingValue) / 2; } -double GetFersThreshold(int fersId, double fallbackThreshold) +double GetFersThreshold(int fersId) { if (fersId < 0) { - return fallbackThreshold; + throw std::runtime_error("Invalid FERS id: " + std::to_string(fersId)); } const auto it = fers_to_thr_map.find(static_cast(fersId)); - return (it != fers_to_thr_map.end()) ? it->second : fallbackThreshold; + if (it == fers_to_thr_map.end()) { + throw std::runtime_error("Missing threshold for FERS id: " + std::to_string(fersId)); + } + + return it->second; +} + +void AccumulateChannelSignal(ChannelSignalMap& channelSignals, + const SipmMapEntry& info, + double thresholdSignalContribution, + double outputSignalContribution) +{ + if (info.boardID < 0 || info.fersId < 0 || info.ch < 0) { + return; + } + + const FersKey fersKey = make_fers_key(info.boardID, info.fersId); + const ChannelKey channelKey{fersKey, info.ch}; + + auto [it, inserted] = channelSignals.emplace( + channelKey, + ChannelAccumulatedSignal{ + fersKey, + info.fersId, + info.ch, + GetFersThreshold(info.fersId), + 0.0, + 0.0, + 0.0, + 0.0 + } + ); + + it->second.thresholdSignal += thresholdSignalContribution; + it->second.outputSignal += outputSignalContribution; +} + +bool IsActivatedFers(const std::unordered_map& fersChannelsOverThreshold, + FersKey key) +{ + const auto it = fersChannelsOverThreshold.find(key); + return it != fersChannelsOverThreshold.end() && + it->second >= FersMultiplicity; } +double SumActivatedSmearedChannelOutput( + const ChannelSignalMap& channelSignals, + const std::unordered_map& fersChannelsOverThreshold) +{ + double total = 0.0; + + for (const auto& [channelKey, channel] : channelSignals) { + if (!IsActivatedFers(fersChannelsOverThreshold, channel.fersKey)) { + continue; + } + + total += channel.smearedOutputSignal; + } + + return total; +} + +unsigned int CountActivatedFers( + const std::unordered_map& fersChannelsOverThreshold) +{ + unsigned int activated = 0; + + for (const auto& [fersKey, nChannels] : fersChannelsOverThreshold) { + if (nChannels >= FersMultiplicity) { + ++activated; + } + } + + return activated; +} + + + +void ApplyPoissonOnePeLossToActivatedChannels( + ChannelSignalMap& channelSignals, + const std::unordered_map& fersChannelsOverThreshold, + TRandom3& rng, + double phePerGeV, + double lossLambda) +{ + if (lossLambda <= 0.0 || phePerGeV <= 0.0) { + return; + } + + const double onePeInSignalUnits = 1.0 / phePerGeV; + + for (auto& [channelKey, channel] : channelSignals) { + if (!IsActivatedFers(fersChannelsOverThreshold, channel.fersKey)) { + continue; + } + + const int nLosses = rng.Poisson(lossLambda); + if (nLosses <= 0) { + continue; + } + + //channel.smearedOutputSignal -= onePeInSignalUnits; + channel.smearedOutputSignal -= nLosses * onePeInSignalUnits; + + if (channel.smearedOutputSignal < 0.0) { + channel.smearedOutputSignal = 0.0; + } + } +} + + + + +/**************************************************************************/ +/************ Utility functions for SiPM map and event display ************/ SipmLookup LoadSipmMap(const std::string& jsonPath) { std::ifstream in(jsonPath); @@ -226,14 +522,6 @@ const SipmMapEntry* FindSipmInfo(const SipmLookup& lookup, return (it == lookup.end()) ? nullptr : &it->second; } -bool IsActivatedFers(const std::unordered_map& fersChannelsOverThreshold, - FersKey key) -{ - const auto it = fersChannelsOverThreshold.find(key); - return it != fersChannelsOverThreshold.end() && - it->second >= FersMultiplicity; -} - std::string make_event_display_name(const std::string& prefix, unsigned int entry) { std::ostringstream os; @@ -297,6 +585,11 @@ void GetSiPMcoordinate(int TowID, SiPM_Y = TowerOffsetY + SiPM_Y; } + + + +/***********************************************************************/ +/************ Utility functions for analysis and output ************/ std::string simTower_to_tbTower(std::string simTower) { return std::to_string(306 + std::stoi(simTower)); @@ -427,6 +720,17 @@ void AppendSummaryCsvLocked(const std::string& csvPath, close(fd); } + + + + +/*********************************/ +/*** + * Main analysis function + * @param energy: beam energy in GeV (used for histogram ranges and output naming) + * @param input: name of the input ROOT file (relative to ../build/) + */ +/*********************************/ void HidraTB25Ana(double energy, const std::string& input) { const std::string sipmMapPath = @@ -468,6 +772,27 @@ void HidraTB25Ana(double energy, const std::string& input) TDirectory* eventDisplayDir = f.mkdir("EventDisplays"); f.cd(); + TRandom3 rng(NoiseRandomSeed); + + const std::string sciNoiseCsvPath = "S_sipm_std_run868_HG_4fers.csv"; + const std::string cerNoiseCsvPath = "C_sipm_std_run868_HG_4fers.csv"; + + ChannelNoiseVector sciChannelNoiseValues; + ChannelNoiseVector cerChannelNoiseValues; + + try { + sciChannelNoiseValues = LoadNoiseValuesFromCsv(sciNoiseCsvPath); + cerChannelNoiseValues = LoadNoiseValuesFromCsv(cerNoiseCsvPath); + + std::cout << "Loaded " << sciChannelNoiseValues.size() + << " S-channel noise values from " << sciNoiseCsvPath << std::endl; + std::cout << "Loaded " << cerChannelNoiseValues.size() + << " C-channel noise values from " << cerNoiseCsvPath << std::endl; + } catch (const std::exception& e) { + std::cerr << e.what() << std::endl; + return; + } + int modcol[NofModulesX * NofModulesY]; int modrow[NofModulesX * NofModulesY]; for (int i = 0; i < NofModulesX * NofModulesY; i++) { @@ -520,59 +845,30 @@ void HidraTB25Ana(double energy, const std::string& input) auto ResidualHistCerY = new TH1F("ResidualHistCerY", "Residual Cer Y; Residual [mm]; Entries", 100, -10, 10); - auto SignalSfersID = new TH1F("SignalSfersID", "Signal S fibers FERS ID; FERS; Signal", 16, 0, 16); - auto SignalCfersID = new TH1F("SignalCfersID", "Signal C fibers FERS ID; FERS; Signal", 16, 0, 16); + auto SignalSfersID = new TH1F("SignalSfersID", "Signal S fibers FERS ID; FERS; Signal", 16, 0.5, 16.5); + auto SignalCfersID = new TH1F("SignalCfersID", "Signal C fibers FERS ID; FERS; Signal", 16, 0.5, 16.5); - auto scieneFersOn = new TH1F("scieneFersOn", "S energy for Activated FERS", 100, 0., bmax); - auto cereneFersOn = new TH1F("cereneFersOn", "C energy for Activated FERS", 100, 0., bmax); + auto scieneFersOn = new TH1F("scieneFersOn", "S energy for Activated FERS", 100, -50., bmax); + auto cereneFersOn = new TH1F("cereneFersOn", "C energy for Activated FERS", 100, -50., bmax); const int nentries = simtree->GetEntries(); std::cout << "Entries " << nentries << std::endl; - int pdg; - simtree->SetBranchAddress("PrimaryPDGID", &pdg); - - double venergy; - simtree->SetBranchAddress("PrimaryParticleEnergy", &venergy); - - double lenergy; - simtree->SetBranchAddress("EscapedEnergyl", &lenergy); - - double denergy; - simtree->SetBranchAddress("EscapedEnergyd", &denergy); - - double edep; - simtree->SetBranchAddress("EnergyTot", &edep); - - double Stot; - simtree->SetBranchAddress("NofPMTScinDet", &Stot); - - double Ctot; - simtree->SetBranchAddress("NofPMTCherDet", &Ctot); - - double PSdep; - simtree->SetBranchAddress("PSEnergy", &PSdep); - - double beamX; - simtree->SetBranchAddress("PrimaryX", &beamX); - - double beamY; - simtree->SetBranchAddress("PrimaryY", &beamY); - - std::vector* TowerE = nullptr; - simtree->SetBranchAddress("VecTowerE", &TowerE); - - std::vector* SPMT = nullptr; - simtree->SetBranchAddress("VecSPMT", &SPMT); - - std::vector* CPMT = nullptr; - simtree->SetBranchAddress("VecCPMT", &CPMT); - - std::vector* SSiPM = nullptr; - simtree->SetBranchAddress("VectorSignals", &SSiPM); - - std::vector* CSiPM = nullptr; - simtree->SetBranchAddress("VectorSignalsCher", &CSiPM); + int pdg; simtree->SetBranchAddress("PrimaryPDGID", &pdg); + double venergy; simtree->SetBranchAddress("PrimaryParticleEnergy", &venergy); + double lenergy; simtree->SetBranchAddress("EscapedEnergyl", &lenergy); + double denergy; simtree->SetBranchAddress("EscapedEnergyd", &denergy); + double edep; simtree->SetBranchAddress("EnergyTot", &edep); + double Stot; simtree->SetBranchAddress("NofPMTScinDet", &Stot); + double Ctot; simtree->SetBranchAddress("NofPMTCherDet", &Ctot); + double PSdep; simtree->SetBranchAddress("PSEnergy", &PSdep); + double beamX; simtree->SetBranchAddress("PrimaryX", &beamX); + double beamY; simtree->SetBranchAddress("PrimaryY", &beamY); + std::vector* TowerE = nullptr; simtree->SetBranchAddress("VecTowerE", &TowerE); + std::vector* SPMT = nullptr; simtree->SetBranchAddress("VecSPMT", &SPMT); + std::vector* CPMT = nullptr; simtree->SetBranchAddress("VecCPMT", &CPMT); + std::vector* SSiPM = nullptr; simtree->SetBranchAddress("VectorSignals", &SSiPM); + std::vector* CSiPM = nullptr; simtree->SetBranchAddress("VectorSignalsCher", &CSiPM); for (unsigned int i = 0; i < static_cast(nentries); i++) { simtree->GetEntry(i); @@ -591,13 +887,14 @@ void HidraTB25Ana(double energy, const std::string& input) double barY_cer = 0.; std::unordered_map fersChannelsOverThreshold; - std::vector sciHitsPerFers; - std::vector cerHitsPerFers; std::vector eventDisplaySciHits; std::vector eventDisplayCerHits; - sciHitsPerFers.reserve(SSiPM->size()); - cerHitsPerFers.reserve(CSiPM->size()); + ChannelSignalMap sciChannelSignals; + ChannelSignalMap cerChannelSignals; + + sciChannelSignals.reserve(SSiPM->size()); + cerChannelSignals.reserve(CSiPM->size()); if (writeEventDisplay) { eventDisplaySciHits.reserve(SSiPM->size()); @@ -614,6 +911,8 @@ void HidraTB25Ana(double energy, const std::string& input) double sciPosWeight = 0.0; double cerPosWeight = 0.0; + + // Start loop over SiPM signals - S channels for (unsigned int n = 0; n < SSiPM->size(); n++) { const double rawPhe = SSiPM->at(n); const double content = rawPhe / sciPheGeV; @@ -628,15 +927,14 @@ void HidraTB25Ana(double energy, const std::string& input) double SiPM_X = 0.0; double SiPM_Y = 0.0; - GetSiPMcoordinate(towID, rowID, colID, SiPM_X, SiPM_Y, "S", grouping); + GetSiPMcoordinate(towID, rowID, colID, SiPM_X, SiPM_Y, "S", grouping); // Not used after JSON SiPM map implementation, kept for reference const std::string tbTower = simTower_to_tbTower(std::to_string(towID)); const int groupedCol = grouped_column(static_cast(colID), grouping); SipmMapS->Fill(colID, towID * NofFibersrow + rowID, content); - const SipmMapEntry* info = - FindSipmInfo(sipmLookup, tbTower, "S", static_cast(rowID), groupedCol, grouping); + const SipmMapEntry* info = FindSipmInfo(sipmLookup, tbTower, "S", static_cast(rowID), groupedCol, grouping); if (info) { SciSiPMCoordinates->Fill(info->x, info->y, content); @@ -647,21 +945,25 @@ void HidraTB25Ana(double energy, const std::string& input) if (info->boardID >= 0 && info->fersId >= 0) { const FersKey key = make_fers_key(info->boardID, info->fersId); - const double threshold = GetFersThreshold(info->fersId, FersThresholdS); - - sciHitsPerFers.push_back({key, rawPhe, content}); + const double outputSignal = ApplySaturation ? apply_sipm_saturation_from_pe(rawPhe, sciPheGeV, NofSipmCells_sci) : content; - if (content > threshold) { - ++fersChannelsOverThreshold[key]; - } + AccumulateChannelSignal( + sciChannelSignals, + *info, + content, + outputSignal + ); if (writeEventDisplay) { eventDisplaySciHits.push_back({key, info->x, info->y, rawPhe, content}); } } } - } + } // end loop over S SiPM channels + + + // Start loop over SiPM signals - C channels for (unsigned int n = 0; n < CSiPM->size(); n++) { const double rawPhe = CSiPM->at(n); const double content = rawPhe / cerPheGeV; @@ -673,7 +975,7 @@ void HidraTB25Ana(double energy, const std::string& input) const unsigned int colID = NofFiberscolumn - static_cast(SiPMID / (NofFibersrow / 2)); const unsigned int rowID = 2 * static_cast(SiPMID % (NofFibersrow / 2)) + 1; - SipmMapC->Fill(colID, towID * NofFibersrow + rowID, content); + SipmMapC->Fill(colID, towID * NofFibersrow + rowID, content); // Not used after JSON SiPM map implementation, kept for reference const std::string tbTower = simTower_to_tbTower(std::to_string(towID)); const int groupedCol = grouped_column(static_cast(colID), grouping); @@ -690,46 +992,72 @@ void HidraTB25Ana(double energy, const std::string& input) if (info->boardID >= 0 && info->fersId >= 0) { const FersKey key = make_fers_key(info->boardID, info->fersId); - const double threshold = GetFersThreshold(info->fersId, FersThresholdC); - - cerHitsPerFers.push_back({key, rawPhe, content}); - - if (content > threshold) { - ++fersChannelsOverThreshold[key]; - } + const double outputSignal = + ApplySaturation + ? apply_sipm_saturation_from_pe(rawPhe, cerPheGeV, NofSipmCells_cer) + : content; + + AccumulateChannelSignal( + cerChannelSignals, + *info, + content, + outputSignal + ); if (writeEventDisplay) { eventDisplayCerHits.push_back({key, info->x, info->y, rawPhe, content}); } } } - } + } // end loop over C SiPM channels - double totsciFersOn = 0.; - for (const auto& hit : sciHitsPerFers) { - if (IsActivatedFers(fersChannelsOverThreshold, hit.key) && - hit.rawPhe >= MinVisiblePheS) { - if(ApplySaturation){ - totsciFersOn += apply_sipm_saturation_from_pe(hit.rawPhe, sciPheGeV, NofSipmCells_sci); - } - else{ - totsciFersOn += hit.signal; - } - } - } + SmearAndCountChannelsOverThreshold( + sciChannelSignals, + fersChannelsOverThreshold, + rng, + sciChannelNoiseValues + ); - double totcerFersOn = 0.; - for (const auto& hit : cerHitsPerFers) { - if (IsActivatedFers(fersChannelsOverThreshold, hit.key) && - hit.rawPhe >= MinVisiblePheC ) { - if(ApplySaturation){ - totcerFersOn += apply_sipm_saturation_from_pe(hit.rawPhe, cerPheGeV, NofSipmCells_cer); - } - else{ - totcerFersOn += hit.signal; - } - } - } + SmearAndCountChannelsOverThreshold( + cerChannelSignals, + fersChannelsOverThreshold, + rng, + cerChannelNoiseValues + ); + + + + // Apply 1-pe loss only after FERS activation is known. + ApplyPoissonOnePeLossToActivatedChannels( + sciChannelSignals, + fersChannelsOverThreshold, + rng, + sciPheGeV, + SciOnePeLossLambda + ); + + ApplyPoissonOnePeLossToActivatedChannels( + cerChannelSignals, + fersChannelsOverThreshold, + rng, + cerPheGeV, + CerOnePeLossLambda + ); + + + + PrintPerEventFersChannelCounts( + i, + sciChannelSignals, + cerChannelSignals, + fersChannelsOverThreshold + ); + + const double totsciFersOn = + SumActivatedSmearedChannelOutput(sciChannelSignals, fersChannelsOverThreshold); + + const double totcerFersOn = + SumActivatedSmearedChannelOutput(cerChannelSignals, fersChannelsOverThreshold); if (sciPosWeight > 0.) { ResidualHistSciX->Fill((barX_sci / sciPosWeight) - beamX); @@ -789,7 +1117,6 @@ void HidraTB25Ana(double energy, const std::string& input) sciene->Fill(totsci); cerene->Fill(totcer); - //totene->Fill(elcont * 0.5 * (totsci + totcer)); totene->Fill(0.5 * (totsci + totcer)); totenec->Fill(picont * (totsci - chi * totcer) / (1 - chi)); totdep->Fill(tottow / 1000.); diff --git a/analysis/plot_hidra.py b/analysis/plot_hidra.py index cb0f3ad..3ffb893 100644 --- a/analysis/plot_hidra.py +++ b/analysis/plot_hidra.py @@ -11,16 +11,29 @@ ROOT.gROOT.SetBatch(True) ROOT.gStyle.SetOptStat(0) +#scaling_s = 1.07 +#scaling_c = 1.05 + #scaling_s = 1.16 #scaling_c = 1.255 -scaling_s = 1.125 -scaling_c = 1.155 -myOutDir = "plots/Attenuation7m_2phe" -extraTex2 = "Fibre attenuation set to 10m" -#myOutDir = "plots/AttenuationTB24_2phe" -extraTex1 = "Min visible signal: >= 2 phe" -#extraTex2 = "Fibre attenuation tuned to TB24 data: 3.7 m (Sci), 3.9 m (Cer)" -extraTex3 = "Correction for multiple photons on the same SiPM cell: OFF" + +scaling_s = 1.14 +scaling_c = 1.27 + +#myOutDir = "plots/Attenuation7m_2phe" +myOutDir = "plots/AttenuationTB24" + +#extraTex2 = "Fibre attenuation set to 10m" +extraTex2 = "Fibre attenuation tuned to TB24 data: 3.7 m (Sci), 3.9 m (Cer)" + +extraTex1 = "TB25 HG data used for tuning channel noise" +#extraTex1 = "No noise" + +extraTex5 = "Correction for multiple photons on the same SiPM cell: ON" +#extraTex5 = " " + +extraTex4 = "Poissonian probability to lose phe per channel: #mu = 3 (Sci), 3 (Cer)" +#extraTex4 = " " def load_summary(csv_path: str) -> pd.DataFrame: data = pd.read_csv(csv_path) @@ -372,7 +385,7 @@ def plot_energy_linearity(data: pd.DataFrame, out_dir: str): tex1.SetTextSize(0.042); tex1.SetLineWidth(2); tex1.Draw("same"); - tex2 = ROOT.TLatex(alignLeft,alignTop-0.04, extraTex2) + tex2 = ROOT.TLatex(alignLeft,alignTop-0.12, extraTex2) tex2.SetNDC(); tex2.SetTextFont(42); tex2.SetTextSize(0.03); @@ -384,6 +397,18 @@ def plot_energy_linearity(data: pd.DataFrame, out_dir: str): tex3.SetTextSize(0.03); tex3.SetLineWidth(2); tex3.Draw("same"); + tex4 = ROOT.TLatex(alignLeft,alignTop-0.04, extraTex4) + tex4.SetNDC(); + tex4.SetTextFont(42); + tex4.SetTextSize(0.03); + tex4.SetLineWidth(2); + tex4.Draw("same"); + tex5 = ROOT.TLatex(alignLeft,alignTop-0.16, extraTex5) + tex5.SetNDC(); + tex5.SetTextFont(42); + tex5.SetTextSize(0.03); + tex5.SetLineWidth(2); + tex5.Draw("same"); legend = ROOT.TLegend(0.450, 0.17, 0.83, 0.37) From deab1ce5797818a3b87d7ff5513685d959706e54 Mon Sep 17 00:00:00 2001 From: AndreaPareti Date: Fri, 15 May 2026 19:23:08 +0200 Subject: [PATCH 4/5] Updates to analysis scripts and TB25 sim --- analysis/HidraAna.C | 244 ++++++- analysis/HidraGeo.h | 82 ++- analysis/HidraTB25Ana.C | 983 +++++++++++++++++++-------- include/DREMTubesGeoPar.hh | 44 +- src/DREMTubesDetectorConstruction.cc | 48 +- 5 files changed, 1072 insertions(+), 329 deletions(-) diff --git a/analysis/HidraAna.C b/analysis/HidraAna.C index 94146f0..9c7a26e 100644 --- a/analysis/HidraAna.C +++ b/analysis/HidraAna.C @@ -16,6 +16,8 @@ // #include #include +#include +#include #include #include #include @@ -23,6 +25,10 @@ #include #include #include +#include +#include +#include +#include // include file with geometry of module #include "HidraGeo.h" @@ -31,10 +37,97 @@ // row ID input ranges from 0 to 31 // need to know fiber type to distinguish global row const unsigned int grouping = 8; -void GetSiPMcoordinate(int TowID, int rowID, int colID_original, double &SiPM_X, double &SiPM_Y, std::string fiber, unsigned int grouping) + +bool IsSiPMModule(int moduleID) +{ + for(int i=0; iGetXaxis()->FindBin(towerXMin + 1.e-9)); + const int lastXBin = std::min(hist->GetNbinsX(), hist->GetXaxis()->FindBin(towerXMax - 1.e-9)); + const int firstYBin = std::max(1, hist->GetYaxis()->FindBin(towerYMin + 1.e-9)); + const int lastYBin = std::min(hist->GetNbinsY(), hist->GetYaxis()->FindBin(towerYMax - 1.e-9)); + + for(int xbin=firstXBin; xbin<=lastXBin; xbin++){ + for(int ybin=firstYBin; ybin<=lastYBin; ybin++){ + const double x = hist->GetXaxis()->GetBinCenter(xbin); + const double y = hist->GetYaxis()->GetBinCenter(ybin); + if(IsInsideAnySiPMModuleCell(x, y, modcol, modrow)) continue; + hist->SetBinContent(xbin, ybin, content); + } + } +} + +void GetSiPMcoordinate(int TowID, int rowID, int colID_original, double &SiPM_X, double &SiPM_Y, std::string fiber, unsigned int grouping, int modcol_sipm, int modrow_sipm) +{ + + // Calculate tower offsets based on the actual column and row positions of the SiPM module + // Module spacing in detector: X = dtubeX*NofFiberscolumn, Y = dtubeY*NofFibersrow + // Placement follows: m_x = dtubeX*NofFiberscolumn*((NofmodulesX-1)/2 - column) + // m_y = dtubeY*NofFibersrow*(row - (NofmodulesY-1)/2) + double TowerOffsetX = dtubeX*NofFiberscolumn*((NofmodulesX-1.0)/2.0 - modcol_sipm); + double TowerOffsetY = dtubeY*NofFibersrow*(modrow_sipm - (NofmodulesY-1.0)/2.0); + // after grouping, there are nfibercolumns/grouping channels unsigned int channel = static_cast(colID_original/grouping); // colID should be such that the corresponding coordinate is in the middle of the channel, @@ -50,11 +143,13 @@ void GetSiPMcoordinate(int TowID, int rowID, int colID_original, double &SiPM_X, } if(fiber == "C"){ - SiPM_X = +moduleX/2 + tuberadius - (tuberadius*2)*colID ; + // Cherenkov fibers are shifted by one additional tube radius in X compared to scintillating fibers + SiPM_X = +moduleX/2 - (tuberadius*2) - (tuberadius*2)*colID ; SiPM_Y = -moduleY/2 + tuberadius + (sq3*tuberadius)*rowID+tuberadius*(2.*sq3m1-1.); //std::cout << "C X: " << SiPM_X << " Y: " << SiPM_Y << std::endl; } + SiPM_X = TowerOffsetX + SiPM_X; SiPM_Y = TowerOffsetY + SiPM_Y; } @@ -62,7 +157,7 @@ void GetSiPMcoordinate(int TowID, int rowID, int colID_original, double &SiPM_X, -void HidraAna(double energy, const string intup){ +void HidraAna(double energy, const string intup, unsigned int EventDisplayEvery = 100){ //Open ntuples string infile = "../build/"+intup; std::cout<<"Using file: "<=0){ modcol[imod]=col; @@ -105,12 +201,19 @@ void HidraAna(double energy, const string intup){ auto totdep = new TH1F("totdep", "totdep",100,0.,bmax); auto leakene = new TH1F("leakene", "leakene",100,0.,0.1); auto chidist = new TH1F("chidist", "chidist",100,0.,1.); - auto mapcalo = new TH2F("mapcalo", "mapcalo",NofModulesX,0.,NofModulesX,NofModulesY,0.,NofModulesY); - auto SipmMapS = new TH2F("SipmMapS", "SipmS; Col; Row", NofSiPMTowersX*NofFiberscolumn, 0, NofSiPMTowersX*NofFiberscolumn, NofSiPMTowersY*NofFibersrow/2, 0, NofSiPMTowersY*NofFibersrow); - auto SipmMapC = new TH2F("SipmMapC", "SipmC; Col; Row", NofSiPMTowersX*NofFiberscolumn, 0, NofSiPMTowersX*NofFiberscolumn, NofSiPMTowersY*NofFibersrow/2, 0, NofSiPMTowersY*NofFibersrow); + auto mapcalo = new TH2F("mapcalo", "mapcalo",NofmodulesX,0.,NofmodulesX,NofmodulesY,0.,NofmodulesY); + auto CaloCoordinatesMap = new TH2F("CaloCoordinatesMap", "Calo Coordinates Map; X [mm]; Y [mm]", + NofmodulesX, -dtubeX*NofFiberscolumn*(NofmodulesX/2.0), dtubeX*NofFiberscolumn*(NofmodulesX/2.0), + NofmodulesY, -dtubeY*NofFibersrow*(NofmodulesY/2.0), dtubeY*NofFibersrow*(NofmodulesY/2.0)); + auto SipmMapS = new TH2F("SipmMapS", "SipmS; Col; Row", NofmodulesX*NofFiberscolumn, 0, NofmodulesX*NofFiberscolumn, NofmodulesY*NofFibersrow/2, 0, NofmodulesY*NofFibersrow); + auto SipmMapC = new TH2F("SipmMapC", "SipmC; Col; Row", NofmodulesX*NofFiberscolumn, 0, NofmodulesX*NofFiberscolumn, NofmodulesY*NofFibersrow/2, 0, NofmodulesY*NofFibersrow); + + // Front view histograms: X-Y coordinates of the fiber hits as seen from +z looking towards detector + double detectorExtentX = dtubeX * NofFiberscolumn * NofmodulesX / 2.0; + double detectorExtentY = dtubeY * NofFibersrow * NofmodulesY / 2.0; - auto SciSiPMCoordinates = new TH2F("SciSiPMCoordinates", "Sci SiPM Coordinates; X [mm]; Y[mm]", NofSiPMTowersX*NofFiberscolumn/grouping, -NofSiPMTowersX*moduleX/2, NofSiPMTowersX*moduleX/2, NofSiPMTowersY*NofFibersrow/2, -NofSiPMTowersY*moduleY/2, NofSiPMTowersY*moduleY/2); - auto CerSiPMCoordinates = new TH2F("CerSiPMCoordinates", "Cer SiPM Coordinates; X [mm]; Y[mm]", NofSiPMTowersX*NofFiberscolumn/grouping, -NofSiPMTowersX*moduleX/2, NofSiPMTowersX*moduleX/2, NofSiPMTowersY*NofFibersrow/2, -NofSiPMTowersY*moduleY/2, NofSiPMTowersY*moduleY/2); + auto SciSiPMCoordinates = new TH2F("SciSiPMCoordinates", "Sci SiPM Coordinates; X [mm]; Y[mm]", NofmodulesX*NofFiberscolumn/grouping, -detectorExtentX, detectorExtentX, NofmodulesY*NofFibersrow/2, -detectorExtentY, detectorExtentY); + auto CerSiPMCoordinates = new TH2F("CerSiPMCoordinates", "Cer SiPM Coordinates; X [mm]; Y[mm]", NofmodulesX*NofFiberscolumn/grouping, -detectorExtentX, detectorExtentX, NofmodulesY*NofFibersrow/2, -detectorExtentY, detectorExtentY); int nentries=simtree->GetEntries(); std::cout<<"Entries "<GetEntries(); i++){ double ecalo=energy-lenergy/1000; simtree->GetEntry(i); + const bool writeEventDisplay = (EventDisplayEvery > 0) && (((i + 1) % EventDisplayEvery) == 0); + TH2F* eventDisplaySci = nullptr; + TH2F* eventDisplayCer = nullptr; + + if(writeEventDisplay){ + std::ostringstream sciTitle; + sciTitle << "Scintillation event display, entry " << i + << "; X [mm]; Y [mm]"; + eventDisplaySci = CreateHybridEventDisplay( + MakeEventDisplayName("EventDisplaySci", i), + sciTitle.str(), + detectorExtentX, + detectorExtentY + ); + + std::ostringstream cerTitle; + cerTitle << "Cherenkov event display, entry " << i + << "; X [mm]; Y [mm]"; + eventDisplayCer = CreateHybridEventDisplay( + MakeEventDisplayName("EventDisplayCer", i), + cerTitle.str(), + detectorExtentX, + detectorExtentY + ); + } + double totsci=0.; double totcer=0.; double tottow=0.; + + + // Sum energy over all MiniModules for(unsigned int j=0; jsize(); j++){ - totsci+=SPMT->at(j)/sciPheGeV; - totcer+=CPMT->at(j)/cerPheGeV; + const double sciTowerContent = SPMT->at(j)/sciPheGeV; + const double cerTowerContent = CPMT->at(j)/cerPheGeV; + sciphe_raw += SPMT->at(j); + cerphe_raw += CPMT->at(j); + totsci+=sciTowerContent; + totcer+=cerTowerContent; tottow+=TowerE->at(j); mapcalo->Fill(modcol[j],modrow[j],TowerE->at(j)/1000/nentries); + double towerX = dtubeX * NofFiberscolumn * ((NofmodulesX-1.0)/2.0 - modcol[j]); + double towerY = dtubeY * NofFibersrow * (modrow[j] - (NofmodulesY-1.0)/2.0); + CaloCoordinatesMap->Fill(towerX, towerY, TowerE->at(j)/1000/nentries); + + if(writeEventDisplay && !IsSiPMModule(j)){ + FillTowerPatch(eventDisplaySci, j, modcol, modrow, sciTowerContent); + FillTowerPatch(eventDisplayCer, j, modcol, modrow, cerTowerContent); + } } for(unsigned int N=0; Nsize(); N++){ // Loop over SiPMs - S Fibers double content = SSiPM->at(N)/sciPheGeV; totsci+=content; + sciphe_raw += SSiPM->at(N); unsigned int towID = static_cast( N/(NofFiberscolumn*NofFibersrow/2) ); unsigned int SiPMID = N%(NofFiberscolumn*NofFibersrow/2); unsigned int colID = static_cast(SiPMID/(NofFibersrow/2)); unsigned int rowID = 2*static_cast(SiPMID%(NofFibersrow/2)); + // Get actual module ID from SiPM tower index + int actual_mod_id = SiPMMod[towID]; + int modcol_sipm = modcol[actual_mod_id]; + int modrow_sipm = modrow[actual_mod_id]; // Get coordinate double SiPM_X, SiPM_Y; - GetSiPMcoordinate(towID, rowID, colID, SiPM_X, SiPM_Y, "S", grouping); + GetSiPMcoordinate(towID, rowID, colID, SiPM_X, SiPM_Y, "S", grouping, modcol_sipm, modrow_sipm); //std::cout << towID << "\tSiPM ID: " << SiPMID << "\tcolID: " << colID << "\trowID: " << rowID << "\tX: " << SiPM_X << "\tY: " << SiPM_Y << std::endl; - //std::cout << "Tower: " << towID << "\tModule row: " << modrow[towID] << "\tSiPMID: " << SiPMID << "\tColumn: " << colID << "\tRow: " << rowID << "\tTotal row: " << towID*NofFibersrow+rowID << std::endl; - //SipmMapS->Fill( modcol[towID]*NofFiberscolumn + colID, modrow[towID]*NofFibersrow+rowID, content); - SipmMapS->Fill( colID, towID*NofFibersrow+rowID, content); + //std::cout << "Tower: " << towID << "\tModule row: " << modrow_sipm << "\tSiPMID: " << SiPMID << "\tColumn: " << colID << "\tRow: " << rowID << "\tTotal row: " << towID*NofFibersrow+rowID << std::endl; + //SipmMapS->Fill( modcol_sipm*NofFiberscolumn + colID, modrow_sipm*NofFibersrow+rowID, content); + SipmMapS->Fill( modcol_sipm*NofFiberscolumn + colID, modrow_sipm*NofFibersrow+rowID, content); SciSiPMCoordinates->Fill(SiPM_X, SiPM_Y, content); - std::cout << "S channel: " << "x: " << SiPM_X << "\t y: " << SiPM_Y << "\t content: " << content << std::endl; + if(writeEventDisplay){ + eventDisplaySci->Fill(SiPM_X, SiPM_Y, content); + } + //std::cout << "S channel: " << "x: " << SiPM_X << "\t y: " << SiPM_Y << "\t content: " << content << std::endl; } - for(unsigned int N=0; Nsize(); N++){ // Loop over SiPMs - S Fibers + for(unsigned int N=0; Nsize(); N++){ // Loop over SiPMs - C Fibers double content = CSiPM->at(N)/cerPheGeV; totcer+=content; + cerphe_raw += CSiPM->at(N); unsigned int towID = static_cast( N/(NofFiberscolumn*NofFibersrow/2) ); unsigned int SiPMID = N%(NofFiberscolumn*NofFibersrow/2); unsigned int colID = static_cast(SiPMID/(NofFibersrow/2)); unsigned int rowID = 2*static_cast(SiPMID%(NofFibersrow/2)) + 1; // Cerenkov fibres on odd rows - SipmMapC->Fill( colID, towID*NofFibersrow+rowID, content); + // Get actual module ID from SiPM tower index + int actual_mod_id = SiPMMod[towID]; + int modcol_sipm = modcol[actual_mod_id]; + int modrow_sipm = modrow[actual_mod_id]; + SipmMapC->Fill( modcol_sipm*NofFiberscolumn + colID, modrow_sipm*NofFibersrow+rowID, content); double SiPM_X, SiPM_Y; - GetSiPMcoordinate(towID, rowID, colID, SiPM_X, SiPM_Y, "C", grouping); - CerSiPMCoordinates->Fill(SiPM_X, SiPM_Y, content); - std::cout << "C channel: " << "x: " << SiPM_X << "\t y: " << SiPM_Y << "\t content: " << content << std::endl; + GetSiPMcoordinate(towID, rowID, colID, SiPM_X, SiPM_Y, "C", grouping, modcol_sipm, modrow_sipm); + CerSiPMCoordinates->Fill(SiPM_X, SiPM_Y, content); + if(writeEventDisplay){ + eventDisplayCer->Fill(SiPM_X, SiPM_Y, content); + } + //std::cout << "C channel: " << "x: " << SiPM_X << "\t y: " << SiPM_Y << "\t content: " << content << std::endl; + + } + if(writeEventDisplay){ + if(eventDisplayDir){ + eventDisplayDir->cd(); + eventDisplaySci->Write(); + eventDisplayCer->Write(); + f.cd(); + } + delete eventDisplaySci; + delete eventDisplayCer; } sciene->Fill(totsci); @@ -201,10 +384,11 @@ void HidraAna(double energy, const string intup){ leakene->Fill(lenergy/1000/energy); chidist->Fill((totsci-ecalo)/(totcer-ecalo)); //std::cout << "totSci: " << totsci << "\t totCer: " << totcer << std::endl; - break; + //break; + } + std::cout << "Phe/GeV Sci: " << sciphe_raw/(venergy/1000)/nentries << "\t Phe/GeV Cer: " << cerphe_raw/(venergy/1000)/nentries << std::endl; - /* totenec->Fit("gaus","Q",""); TF1 *fit1 = totenec->GetFunction("gaus"); diff --git a/analysis/HidraGeo.h b/analysis/HidraGeo.h index 98193f0..fbbe50d 100644 --- a/analysis/HidraGeo.h +++ b/analysis/HidraGeo.h @@ -14,7 +14,7 @@ const double tuberadius = 1.0; - +/* // 70 modules -> TB25 const int NofModulesX = 5; const int NofModulesY = 18; @@ -46,11 +46,91 @@ const int NofFibersrow = 16; const int NoModulesActive=80; const double moduleZ = (2500.); const bool irot=false; +const int NoFibersTower=NofFiberscolumn*NofFibersrow/2;*/ + + + + + + +// 70 modules -> TB25 +const int NofmodulesX = 5; +const int NofmodulesY = 18; +// 80 miniM -> But different geometry-> move one minimodules from row to corner, horizontal placement +const int modflag[90]={-1, -1, 0, -1, -1, + -1, 1, 2, 3, -1, + -1, 4, 5, 6, -1, + -1, 7, 8, 9, -1, + 10, 11, 12, 13, 14, + 15, 16, 17, 18, 19, + 20, 21, 22, 23, 24, + 25, 26, 27, 28, 29, + 30, 31, 32, 33, 34, + 35, 36, 37, 38, 39, + 40, 41, 42, 43, 44, + 45, 46, 47, 48, 49, + 50, 51, 52, 53, 54, + 55, 56, 57, 58, 59, + -1, 60, 61, 62, -1, + -1, 63, 64, 65, -1, + -1, 66, 67, 68, -1, + -1, -1, 69, -1, -1}; +const int NoModulesSiPM=8; +const int NofSiPMTowersX = 1; +const int NofSiPMTowersY = 8; +const int SiPMMod[8]={17,22,27,32,37,42,47,52}; +const int NofFiberscolumn = 64; +const int NofFibersrow = 16; +const int NoModulesActive=70; +const double moduleZ = (2500.); +const bool irot=false; const int NoFibersTower=NofFiberscolumn*NofFibersrow/2; + +/* +// 80 modules -> TB26 (HiDRa geometry, SiPMs placed in two towers) +const int NofModulesX = 5; +const int NofModulesY = 20; +const int modflag[100]={-1, -1, 0, -1, -1, + -1, 1, 2, 3, -1, + -1, 4, 5, 6, -1, + -1, 7, 8, 9, -1, + 10, 1, 12, 13, 14, + 15, 16, 17, 18, 19, + 20, 21, 22, 23, 24, + 25, 26, 27, 28, 29, + 30, 31, 32, 33, 34, + 35, 36, 37, 38, 39, + 40, 41, 42, 43, 44, + 45, 46, 47, 48, 49, + 50, 51, 52, 53, 54, + 55, 56, 57, 58, 59, + 60, 61, 62, 63, 64, + 65, 66, 67, 68, 69, + -1, 70, 71, 72, -1, + -1, 73, 74, 75, -1, + -1, 76, 77, 78, -1, + -1, -1, 79, -1, -1}; +const int NoModulesSiPM=10; +const int SiPMMod[10]={27,32,37,42,47,52, 26, 31, 36, 41}; + +const int NofFiberscolumn = 64; +const int NofFibersrow = 16; +const int NofSiPMTowersX = 2; +const int NofSiPMTowersY = 10; +const int NoModulesActive=80; +const double moduleZ = (2500.); +const bool irot=false; +const int NoFibersTower=NofFiberscolumn*NofFibersrow/2;*/ + + + + + + const int NofSciSiPM=(NofFibersrow*NofFiberscolumn*NoModulesSiPM)/2; const int NofCerSiPM=(NofFibersrow*NofFiberscolumn*NoModulesSiPM)/2; diff --git a/analysis/HidraTB25Ana.C b/analysis/HidraTB25Ana.C index 8495dd1..0c18608 100644 --- a/analysis/HidraTB25Ana.C +++ b/analysis/HidraTB25Ana.C @@ -1,5 +1,3 @@ -// file: HidraTB25Ana.C - #include #include #include @@ -20,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -32,23 +31,20 @@ using json = nlohmann::json; -// Global Constants and parameters -// -// Minimum visible signal per channel (>=) -// Set to 0.0 or 1.0 to disable the cut. -const double MinVisiblePheS = 0.0; // photoelectrons -const double MinVisiblePheC = 0.0; // photoelectrons +// ----------------------------------------------------------------------------- +// Global constants and parameters +// ----------------------------------------------------------------------------- + +const double MinVisiblePheS = 0.0; +const double MinVisiblePheC = 0.0; -// Scintillator on even rows (start from 0), Cerenkov on odd rows const unsigned int grouping = 8; -// Write event display every N entries (0 to disable) const unsigned int EventDisplayEvery = 100; -// A FERS is ON if at least 2 smeared channels across S or C (both are good) are above that FERS threshold. const unsigned int FersMultiplicity = 2; const double chi = 0.38; const double sciPheGeV = 119.001; // tb24 -const double cerPheGeV = 29.4; // tb24 +const double cerPheGeV = 29.4; // tb24 //const double sciPheGeV = 178.501; //const double cerPheGeV = 43; @@ -60,28 +56,80 @@ const double cer_pde = 0.38; const double elcont = 1.005; const double picont = 1.028; -// Separate fallback noise sigmas for S and C -const unsigned int NoiseRandomSeed = 12345; // fixed seed for reproducibility -// In this parametrised simulation, rawPhe is stored as "impinging optical photons" -// To account for multiple photons impinging on the same SiPM cell, -// we apply saturation correction to get "fired cells", which is the actual observable. +enum class NoiseCorrelationMode { + UncorrelatedByChannel, + CorrelatedWithinFers +}; + +enum class NoiseDistributionMode { + Gaussian, + LogNormal +}; + +///////////////////////////////////////////////////////// +// Paremeters for changing simulation behaviour const bool ApplySaturation = true; -// -const bool ApplyTbNoise = true; -// decide whether to print per-event FERS channel counts (for debugging) -bool printFersLog = false; +const bool ApplyTbNoise = false; +const bool AddPedestalToEnergyContribution = false; +const bool ApplyPedestalSubtraction = false; +const double ScalePedestalSubtractionFactorS = 1.0; // 1.09 for correlated +const double ScalePedestalSubtractionFactorC = 1.0; // 1.055 for correlated +//const double ScalePedestalSubtractionFactorS = 1.1; // 1.09 for correlated +//const double ScalePedestalSubtractionFactorC = 1.055; // 1.055 for correlated +const double FersThresholdScaleFactor = 1.; +const NoiseCorrelationMode NoiseMode = NoiseCorrelationMode::UncorrelatedByChannel; // use uncorrelated noise +//const NoiseCorrelationMode NoiseMode = NoiseCorrelationMode::CorrelatedWithinFers; // use correlated noise within FERS + +//const NoiseDistributionMode NoiseDistribution = NoiseDistributionMode::LogNormal; +const NoiseDistributionMode NoiseDistribution = NoiseDistributionMode::Gaussian; +const bool printSmearingLog = false; +const bool printFersLog = false; +/////////////////////////////////////////////////////////// +const unsigned int ChannelsPerFers = 64; +const unsigned int TotalFers = 16; +const unsigned int TotalCalibChannels = ChannelsPerFers * TotalFers; -// lambda = mean number of loss opportunities per channel. -// Probability to lose 1 pe = 1 - exp(-lambda). -const double SciOnePeLossLambda = 3; // set > 0 to enable -const double CerOnePeLossLambda = 3; // set > 0 to enable + +using NoiseGeneratorBank = std::vector; + +// ----------------------------------------------------------------------------- +// Utility helpers +// ----------------------------------------------------------------------------- + +// ----------------------------------------------------------------------------- +// 3) Add deterministic RNG builders +// Note: ROOT treats seed=0 specially, so channel seed uses idx+1. +// ----------------------------------------------------------------------------- + +unsigned int MakeChannelNoiseSeed(unsigned int calibIndex){return calibIndex + 1U;} +unsigned int MakeFersNoiseSeed(unsigned int fersId){return fersId;} + +NoiseGeneratorBank BuildChannelNoiseGenerators() +{ + NoiseGeneratorBank generators(TotalCalibChannels); + for (unsigned int calibIndex = 0; calibIndex < TotalCalibChannels; ++calibIndex) { + generators[calibIndex].SetSeed(MakeChannelNoiseSeed(calibIndex)); + } + return generators; +} + +NoiseGeneratorBank BuildFersNoiseGenerators() +{ + NoiseGeneratorBank generators(TotalFers + 1U); + for (unsigned int fersId = 1; fersId <= TotalFers; ++fersId) { + generators[fersId].SetSeed(MakeFersNoiseSeed(fersId)); + } + return generators; +} + + double apply_sipm_saturation_from_pe(double pe, double phePerGeV, double nCells) @@ -94,6 +142,8 @@ double apply_sipm_saturation_from_pe(double pe, return firedCells / phePerGeV; } + +// FERS activation values (pedestal subtracted) std::map fers_to_thr_map = { { 1, 0.17 }, { 2, 0.20 }, @@ -113,62 +163,37 @@ std::map fers_to_thr_map = { {16, 0.06 } }; -using ChannelNoiseVector = std::vector; - -std::string Trim(const std::string& value) -{ - const std::size_t first = value.find_first_not_of(" \t\r\n"); - if (first == std::string::npos) { - return ""; - } - - const std::size_t last = value.find_last_not_of(" \t\r\n"); - return value.substr(first, last - first + 1); -} - -ChannelNoiseVector LoadNoiseValuesFromCsv(const std::string& csvPath) -{ - std::ifstream in(csvPath); - if (!in) { - throw std::runtime_error("Cannot open noise CSV: " + csvPath); - } - - ChannelNoiseVector values; - std::string line; - - while (std::getline(in, line)) { - const std::string trimmed = Trim(line); - if (trimmed.empty()) { - continue; - } - - const std::size_t commaPos = trimmed.find(','); - const std::string token = - Trim(trimmed.substr(0, commaPos == std::string::npos ? trimmed.size() : commaPos)); - - values.push_back(std::stod(token)); - } - - if (values.empty()) { - throw std::runtime_error("Noise CSV is empty: " + csvPath); - } - - return values; -} - -double GetChannelNoiseSigma(const ChannelNoiseVector& noiseValues, int ch) -{ - if (ch < 0 || static_cast(ch) >= noiseValues.size()) { - std::ostringstream os; - os << "Missing per-channel noise sigma for channel " << ch - << " (CSV size = " << noiseValues.size() << ")"; - throw std::runtime_error(os.str()); - } +/* +// FERS activation values ( non-pedestal subtracted) +std::map fers_to_thr_map = { +{1, 0.2225}, +{2, 0.24}, +{3, 0.3508}, +{4, 0.3508}, +{5, 0.1641}, +{6, 0.1525}, +{7, 0.1405}, +{8, 0.0942}, +{9, 0.0592}, +{10, 0.0592}, +{11, 0.0358}, +{12, 0.0825}, +{13, 0.0825}, +{14, 0.0825}, +{15, 0.1413}, +{16, 0.1413} +};*/ + +struct PedestalHgEntry { + double medianAdc = 0.0; + double rmsAdc = 0.0; +}; - return noiseValues[static_cast(ch)]; -} +using PedestalHgVector = std::vector; +using AdcToGeVVector = std::vector; struct SipmMapEntry { + int calibIndex = -1; int boardID = -1; std::string type; int row = -1; @@ -180,12 +205,11 @@ struct SipmMapEntry { std::string module_name; double x_local = 0.0; double y_local = 0.0; -}; -struct FersChannelHit { - uint64_t key = 0; - double rawPhe = 0.0; - double signal = 0.0; + double pedestalMedianHgAdc = 0.0; + double pedestalRmsHgAdc = 0.0; + double pedestalMedianHgGeV = 0.0; + double pedestalRmsHgGeV = 0.0; }; struct EventDisplayHit { @@ -213,6 +237,7 @@ struct HistMomentSummary { using SipmLookup = std::unordered_map; using FersKey = uint64_t; +// Constrct a unique key for each FERS+channel combination for easier lookup and accumulation struct ChannelKey { FersKey fersKey = 0; int ch = -1; @@ -236,144 +261,218 @@ struct ChannelAccumulatedSignal { FersKey fersKey = 0; int fersId = -1; int ch = -1; + int calibIndex = -1; + bool isSci = false; double threshold = 0.0; double thresholdSignal = 0.0; double outputSignal = 0.0; + double pedestalMean = 0.0; + double noiseSigma = 0.0; double smearedThresholdSignal = 0.0; double smearedOutputSignal = 0.0; }; + using ChannelSignalMap = std::unordered_map; -void SmearAndCountChannelsOverThreshold( - ChannelSignalMap& channelSignals, - std::unordered_map& fersChannelsOverThreshold, - TRandom3& rng, - const ChannelNoiseVector& channelNoiseValues) +std::string make_sipm_key(const std::string& tower, + const std::string& type, + int row, + int column) { - for (auto& [channelKey, channel] : channelSignals) { - const double sigma = GetChannelNoiseSigma(channelNoiseValues, channel.ch); + std::ostringstream os; + os << tower << '|' << type << '|' << row << '|' << column; + return os.str(); +} - double noise = 0.; - if(ApplyTbNoise){ - noise = (sigma > 0.0) ? rng.Gaus(0.0, sigma) : 0.0; - } +FersKey make_fers_key(int boardID, int fersId) +{ + return (static_cast(static_cast(boardID)) << 32) | + static_cast(fersId); +} +int board_id_from_fers_key(FersKey key) +{ + return static_cast(static_cast(key >> 32)); +} - channel.smearedThresholdSignal = channel.thresholdSignal + noise; - channel.smearedOutputSignal = channel.outputSignal + noise; +int grouped_column(int colID_original, unsigned int groupingValue) +{ + return static_cast(colID_original / groupingValue); +} - // FERS activation logic: if smeared threshold signal is above threshold, - // count this channel for FERS multiplicity - if (channel.smearedThresholdSignal > channel.threshold) { - ++fersChannelsOverThreshold[channel.fersKey]; - } +double GetFersThreshold(int fersId) +{ + if (fersId < 0) { + throw std::runtime_error("Invalid FERS id: " + std::to_string(fersId)); + } + + const auto it = fers_to_thr_map.find(static_cast(fersId)); + if (it == fers_to_thr_map.end()) { + throw std::runtime_error("Missing threshold for FERS id: " + std::to_string(fersId)); } + + //return it->second; + return it->second * FersThresholdScaleFactor; // try scaling to consider pedestal non-subtracted values } +// ----------------------------------------------------------------------------- +// Calibration loaders +// ----------------------------------------------------------------------------- +PedestalHgVector LoadHgPedestalsFromJson(const std::string& jsonPath) +{ + std::ifstream in(jsonPath); + if (!in) { + throw std::runtime_error("Cannot open pedestal JSON: " + jsonPath); + } + json j; + in >> j; + if (!j.is_object()) { + throw std::runtime_error("Pedestal JSON must be an object: " + jsonPath); + } + std::size_t maxIndex = 0; + for (auto it = j.begin(); it != j.end(); ++it) { + const int idx = std::stoi(it.key()); + if (idx < 0) { + throw std::runtime_error("Negative pedestal index in: " + jsonPath); + } + maxIndex = std::max(maxIndex, static_cast(idx)); + } + PedestalHgVector values(maxIndex + 1); + std::vector found(maxIndex + 1, false); + for (auto it = j.begin(); it != j.end(); ++it) { + const int idx = std::stoi(it.key()); + const json& node = it.value(); -int board_id_from_fers_key(FersKey key) -{ - return static_cast(static_cast(key >> 32)); + values[static_cast(idx)].medianAdc = node.value("median_HG", 0.0); + values[static_cast(idx)].rmsAdc = node.value("iqr_eff_HG", 0.0); + found[static_cast(idx)] = true; + } + + for (std::size_t i = 0; i < found.size(); ++i) { + if (!found[i]) { + std::ostringstream os; + os << "Missing HG pedestal entry for channel " << i + << " in " << jsonPath; + throw std::runtime_error(os.str()); + } + } + + return values; } -void PrintPerEventFersChannelCounts( - unsigned int eventIndex, - const ChannelSignalMap& sciChannelSignals, - const ChannelSignalMap& cerChannelSignals, - const std::unordered_map& fersChannelsOverThreshold) +AdcToGeVVector LoadAdcToGeVFromJson(const std::string& jsonPath) { - std::map, FersKey> fersSeen; - - auto collectFers = [&](const ChannelSignalMap& channelSignals) { - for (const auto& [channelKey, channel] : channelSignals) { - const int boardID = board_id_from_fers_key(channel.fersKey); - fersSeen[{boardID, channel.fersId}] = channel.fersKey; - } - }; + std::ifstream in(jsonPath); + if (!in) { + throw std::runtime_error("Cannot open ADC->GeV JSON: " + jsonPath); + } - collectFers(sciChannelSignals); - collectFers(cerChannelSignals); + json j; + in >> j; - if (printFersLog) { - std::cout << "Event " << eventIndex << ":\n"; - for (const auto& [id, key] : fersSeen) { - const int boardID = id.first; - const int fersId = id.second; + if (!j.is_array()) { + throw std::runtime_error("ADC->GeV JSON must be an array: " + jsonPath); + } - const auto it = fersChannelsOverThreshold.find(key); - const unsigned int nChannelsOverThreshold = - (it != fersChannelsOverThreshold.end()) ? it->second : 0U; + AdcToGeVVector values; + values.reserve(j.size()); - std::cout - << " board " << boardID - << ", FERS " << fersId - << " -> channels above threshold = " << nChannelsOverThreshold; - - if (nChannelsOverThreshold >= FersMultiplicity) { - std::cout << " [ON]"; - } else { - std::cout << " [OFF]"; - } + for (const auto& item : j) { + values.push_back(item.get()); + } - std::cout << '\n'; - } + if (values.empty()) { + throw std::runtime_error("ADC->GeV JSON is empty: " + jsonPath); } -} -std::string make_sipm_key(const std::string& tower, - const std::string& type, - int row, - int column) -{ - std::ostringstream os; - os << tower << '|' << type << '|' << row << '|' << column; - return os.str(); + return values; } -FersKey make_fers_key(int boardID, int fersId) +// not used in current version +std::size_t BuildCalibrationIndex(const SipmMapEntry& entry) { - return (static_cast(static_cast(boardID)) << 32) | - static_cast(fersId); -} + if (entry.fersId <= 0 || entry.fersId > static_cast(TotalFers)) { + std::ostringstream os; + os << "Invalid FERS id " << entry.fersId; + throw std::runtime_error(os.str()); + } -int grouped_column(int colID_original, unsigned int groupingValue) -{ - return static_cast(colID_original / groupingValue); -} + if (entry.ch < 0 || entry.ch >= static_cast(ChannelsPerFers)) { + std::ostringstream os; + os << "Invalid channel " << entry.ch + << " for FERS " << entry.fersId; + throw std::runtime_error(os.str()); + } -int grouped_to_json_column(int groupedCol, unsigned int groupingValue) -{ - return groupedCol * static_cast(groupingValue) + static_cast(groupingValue) / 2; + const std::size_t idx = + (static_cast(entry.fersId) - 1U) * ChannelsPerFers + + static_cast(entry.ch); + + if (idx >= TotalCalibChannels) { + std::ostringstream os; + os << "Calibration index out of range: " << idx + << " from FERS " << entry.fersId + << ", ch " << entry.ch; + throw std::runtime_error(os.str()); + } + + return idx; } -double GetFersThreshold(int fersId) +void AttachHgPedestalNoiseToSipmMap(SipmLookup& lookup, + const PedestalHgVector& pedestalsHg, + const AdcToGeVVector& adcToGeV) { - if (fersId < 0) { - throw std::runtime_error("Invalid FERS id: " + std::to_string(fersId)); + if (pedestalsHg.size() != TotalCalibChannels) { + std::ostringstream os; + os << "Expected " << TotalCalibChannels + << " HG pedestal entries, got " << pedestalsHg.size(); + throw std::runtime_error(os.str()); } - const auto it = fers_to_thr_map.find(static_cast(fersId)); - if (it == fers_to_thr_map.end()) { - throw std::runtime_error("Missing threshold for FERS id: " + std::to_string(fersId)); + if (adcToGeV.size() != TotalCalibChannels) { + std::ostringstream os; + os << "Expected " << TotalCalibChannels + << " ADC->GeV entries, got " << adcToGeV.size(); + throw std::runtime_error(os.str()); } - return it->second; + for (auto& [key, entry] : lookup) { + if (entry.calibIndex < 0 || entry.calibIndex >= static_cast(TotalCalibChannels)) { + std::ostringstream os; + os << "Invalid calibIndex " << entry.calibIndex + << " for map entry " << key; + throw std::runtime_error(os.str()); + } + + const std::size_t idx = static_cast(entry.calibIndex); + + entry.pedestalMedianHgAdc = pedestalsHg[idx].medianAdc; + entry.pedestalRmsHgAdc = pedestalsHg[idx].rmsAdc; + entry.pedestalMedianHgGeV = pedestalsHg[idx].medianAdc * adcToGeV[idx]; + entry.pedestalRmsHgGeV = pedestalsHg[idx].rmsAdc * adcToGeV[idx]; + } } + + +// ----------------------------------------------------------------------------- +// Signal accumulation in each channel and smearing +// Apply saturation to a single signal path for both threshold and output consistency +// ----------------------------------------------------------------------------- void AccumulateChannelSignal(ChannelSignalMap& channelSignals, const SipmMapEntry& info, - double thresholdSignalContribution, - double outputSignalContribution) + double signalContribution) { - if (info.boardID < 0 || info.fersId < 0 || info.ch < 0) { + if (info.boardID < 0 || info.fersId < 0 || info.ch < 0 || info.calibIndex < 0) { return; } @@ -386,18 +485,192 @@ void AccumulateChannelSignal(ChannelSignalMap& channelSignals, fersKey, info.fersId, info.ch, + info.calibIndex, + info.type == "S", GetFersThreshold(info.fersId), 0.0, 0.0, + info.pedestalMedianHgGeV, + info.pedestalRmsHgGeV, 0.0, 0.0 } ); - it->second.thresholdSignal += thresholdSignalContribution; - it->second.outputSignal += outputSignalContribution; + if (!inserted && it->second.calibIndex != info.calibIndex) { + std::ostringstream os; + os << "Inconsistent calibIndex for accumulated channel: existing=" + << it->second.calibIndex << ", new=" << info.calibIndex + << " for FERS " << info.fersId << ", ch " << info.ch; + throw std::runtime_error(os.str()); + } + + // Use single signal for both threshold and output to ensure consistency + it->second.thresholdSignal += signalContribution; + it->second.outputSignal += signalContribution; +} + +double DrawUnitGaussianNoiseForChannel(const ChannelAccumulatedSignal& channel, + NoiseGeneratorBank& channelNoiseGenerators) +{ + if (channel.calibIndex < 0 || + channel.calibIndex >= static_cast(channelNoiseGenerators.size())) { + std::ostringstream os; + os << "Invalid calibIndex " << channel.calibIndex + << " while drawing channel noise"; + throw std::runtime_error(os.str()); + } + + return channelNoiseGenerators[static_cast(channel.calibIndex)].Gaus(0.0, 1.0); +} + +double DrawUnitGaussianNoiseForFers( + const ChannelAccumulatedSignal& channel, + NoiseGeneratorBank& fersNoiseGenerators, + std::unordered_map& fersEventNoiseCache) +{ + const auto it = fersEventNoiseCache.find(channel.fersKey); + if (it != fersEventNoiseCache.end()) { + return it->second; + } + + if (channel.fersId <= 0 || + channel.fersId >= static_cast(fersNoiseGenerators.size())) + { + throw std::runtime_error( + "Invalid FERS id " + std::to_string(channel.fersId) + + " while drawing correlated FERS noise" + ); + } + + + const double sharedUnitNoise = + fersNoiseGenerators[static_cast(channel.fersId)].Gaus(0.0, 1.0); + + fersEventNoiseCache[channel.fersKey] = sharedUnitNoise; + return sharedUnitNoise; +} + + +struct LogNormalParameters { + double mu = 0.0; + double sigma = 0.0; +}; + +LogNormalParameters GetLogNormalParametersFromMeanRms(double mean, double rms) +{ + if (mean <= 0.0 || rms <= 0.0) { + return {0.0, 0.0}; + } + + const double variance = rms * rms; + + LogNormalParameters pars; + pars.sigma = std::sqrt(std::log(1.0 + variance / (mean * mean))); + pars.mu = std::log(mean) - 0.5 * pars.sigma * pars.sigma; + + return pars; } +double DrawLogNormalWithMeanRms(TRandom3& rng, double mean, double rms) +{ + if (mean <= 0.0) { + return 0.0; + } + + if (rms <= 0.0) { + return mean; + } + + const LogNormalParameters pars = + GetLogNormalParametersFromMeanRms(mean, rms); + + return std::exp(rng.Gaus(pars.mu, pars.sigma)); +} + + + +// ----------------------------------------------------------------------------- +//Smearing: noise first, activation decision first, no pedestal subtraction here +void SmearAndCountChannelsOverThreshold( + ChannelSignalMap& channelSignals, + std::unordered_map& fersChannelsOverThreshold, + NoiseCorrelationMode noiseMode, + NoiseGeneratorBank& channelNoiseGenerators, + NoiseGeneratorBank& fersNoiseGenerators, + std::unordered_map& fersEventNoiseCache) +{ + + for (auto& [channelKey, channel] : channelSignals) { + + double pedestalContribution = 0.0; + + if (AddPedestalToEnergyContribution) { + + if (!ApplyTbNoise || channel.noiseSigma <= 0.0) { + pedestalContribution = channel.pedestalMean; + } + + else if (NoiseDistribution == NoiseDistributionMode::Gaussian) { + double unitNoise = 0.0; + + if (noiseMode == NoiseCorrelationMode::UncorrelatedByChannel) { + unitNoise = DrawUnitGaussianNoiseForChannel( + channel, + channelNoiseGenerators + ); + } else if (noiseMode == NoiseCorrelationMode::CorrelatedWithinFers) { + unitNoise = DrawUnitGaussianNoiseForFers( + channel, + fersNoiseGenerators, + fersEventNoiseCache + ); + } + + pedestalContribution = + channel.pedestalMean + unitNoise * channel.noiseSigma; + } + + else if (NoiseDistribution == NoiseDistributionMode::LogNormal) { + if (noiseMode == NoiseCorrelationMode::UncorrelatedByChannel) { + pedestalContribution = DrawLogNormalWithMeanRms( + channelNoiseGenerators[static_cast(channel.calibIndex)], + channel.pedestalMean, + channel.noiseSigma + ); + } else if (noiseMode == NoiseCorrelationMode::CorrelatedWithinFers) { + pedestalContribution = DrawLogNormalWithMeanRms( + fersNoiseGenerators[static_cast(channel.fersId)], + channel.pedestalMean, + channel.noiseSigma + ); + } + } + } + + channel.smearedThresholdSignal = + channel.thresholdSignal + pedestalContribution; + + channel.smearedOutputSignal = + channel.outputSignal + pedestalContribution; + + + + + + if (channel.smearedThresholdSignal > channel.threshold) { + ++fersChannelsOverThreshold[channel.fersKey]; + } + + + + + } +} + + + + bool IsActivatedFers(const std::unordered_map& fersChannelsOverThreshold, FersKey key) { @@ -438,35 +711,115 @@ unsigned int CountActivatedFers( } - -void ApplyPoissonOnePeLossToActivatedChannels( +// ----------------------------------------------------------------------------- +// 4) subtract pedestal only for activated FERS +void SubtractPedestalMedianFromActivatedChannels( ChannelSignalMap& channelSignals, + const std::unordered_map& fersChannelsOverThreshold) +{ + for (auto& [channelKey, channel] : channelSignals) { + if (!IsActivatedFers(fersChannelsOverThreshold, channel.fersKey)) { + continue; + } + + if(ApplyPedestalSubtraction) + { + const double scale = channel.isSci + ? ScalePedestalSubtractionFactorS + : ScalePedestalSubtractionFactorC; + + //channel.smearedOutputSignal -= channel.pedestalMean; + if(channel.outputSignal>0 and printSmearingLog){ + std::cout << "Subtracting pedestal median for FERS " << channel.fersId + << ", ch " << channel.ch + << ", true signal: " << channel.outputSignal + << "\t smeared signal: " << channel.smearedOutputSignal + << "\t pedestal median: " << channel.pedestalMean + << "\t scale: " << scale << std::endl;} + channel.smearedOutputSignal -= channel.pedestalMean * scale; + } + // Allow negative values: downward fluctuations naturally subtract from total energy + } +} + +void FillActivatedChannelOutputArrays( + const ChannelSignalMap& channelSignals, const std::unordered_map& fersChannelsOverThreshold, - TRandom3& rng, - double phePerGeV, - double lossLambda) + std::vector& activatedChannelOutput, + std::vector& activatedChannelOutputTrue, + std::vector& activatedChannelOutputSmeared) { - if (lossLambda <= 0.0 || phePerGeV <= 0.0) { - return; + if (activatedChannelOutput.size() != activatedChannelOutputTrue.size() || + activatedChannelOutput.size() != activatedChannelOutputSmeared.size()) { + throw std::runtime_error("Activated channel output arrays must have the same size"); } - const double onePeInSignalUnits = 1.0 / phePerGeV; + for (const auto& [channelKey, channel] : channelSignals) { + if (channel.calibIndex < 0 || + channel.calibIndex >= static_cast(activatedChannelOutput.size())) { + std::ostringstream os; + os << "Invalid calibIndex " << channel.calibIndex + << " while filling activated channel output arrays"; + throw std::runtime_error(os.str()); + } - for (auto& [channelKey, channel] : channelSignals) { if (!IsActivatedFers(fersChannelsOverThreshold, channel.fersKey)) { continue; } - const int nLosses = rng.Poisson(lossLambda); - if (nLosses <= 0) { - continue; + const std::size_t idx = static_cast(channel.calibIndex); + + activatedChannelOutputTrue[idx] = channel.outputSignal; + activatedChannelOutputSmeared[idx] = channel.smearedOutputSignal; + + // keep the legacy branch identical to the smeared output + activatedChannelOutput[idx] = channel.smearedOutputSignal; + } +} + + + + +void PrintPerEventFersChannelCounts( + unsigned int eventIndex, + const ChannelSignalMap& sciChannelSignals, + const ChannelSignalMap& cerChannelSignals, + const std::unordered_map& fersChannelsOverThreshold) +{ + std::map, FersKey> fersSeen; + + auto collectFers = [&](const ChannelSignalMap& channelSignals) { + for (const auto& [channelKey, channel] : channelSignals) { + const int boardID = board_id_from_fers_key(channel.fersKey); + fersSeen[{boardID, channel.fersId}] = channel.fersKey; } + }; + + collectFers(sciChannelSignals); + collectFers(cerChannelSignals); + + if (printFersLog) { + std::cout << "Event " << eventIndex << ":\n"; + for (const auto& [id, key] : fersSeen) { + const int boardID = id.first; + const int fersId = id.second; + + const auto it = fersChannelsOverThreshold.find(key); + const unsigned int nChannelsOverThreshold = + (it != fersChannelsOverThreshold.end()) ? it->second : 0U; - //channel.smearedOutputSignal -= onePeInSignalUnits; - channel.smearedOutputSignal -= nLosses * onePeInSignalUnits; + std::cout + << " board " << boardID + << ", FERS " << fersId + << " -> channels above threshold = " << nChannelsOverThreshold; + + if (nChannelsOverThreshold >= FersMultiplicity) { + std::cout << " [ON]"; + } else { + std::cout << " [OFF]"; + } - if (channel.smearedOutputSignal < 0.0) { - channel.smearedOutputSignal = 0.0; + std::cout << '\n'; } } } @@ -474,8 +827,9 @@ void ApplyPoissonOnePeLossToActivatedChannels( -/**************************************************************************/ -/************ Utility functions for SiPM map and event display ************/ +// ----------------------------------------------------------------------------- +// SiPM map and event display helpers +// ----------------------------------------------------------------------------- SipmLookup LoadSipmMap(const std::string& jsonPath) { std::ifstream in(jsonPath); @@ -491,6 +845,7 @@ SipmLookup LoadSipmMap(const std::string& jsonPath) const json& node = it.value(); SipmMapEntry entry; + entry.calibIndex = std::stoi(it.key()); entry.boardID = node.value("boardID", -1); entry.type = node.value("type", ""); entry.row = node.value("row", -1); @@ -585,11 +940,10 @@ void GetSiPMcoordinate(int TowID, SiPM_Y = TowerOffsetY + SiPM_Y; } +// ----------------------------------------------------------------------------- +// Analysis/output helpers +// ----------------------------------------------------------------------------- - - -/***********************************************************************/ -/************ Utility functions for analysis and output ************/ std::string simTower_to_tbTower(std::string simTower) { return std::to_string(306 + std::stoi(simTower)); @@ -720,32 +1074,45 @@ void AppendSummaryCsvLocked(const std::string& csvPath, close(fd); } +// ----------------------------------------------------------------------------- +// Main analysis +// ----------------------------------------------------------------------------- - - - -/*********************************/ -/*** - * Main analysis function - * @param energy: beam energy in GeV (used for histogram ranges and output naming) - * @param input: name of the input ROOT file (relative to ../build/) - */ -/*********************************/ void HidraTB25Ana(double energy, const std::string& input) { + // read SiPM map and calibration files const std::string sipmMapPath = "/home/apareti/HidraSim2025/TBDataPreparation/2025_SPS/MapAndCalibration/sipm_map.json"; + const std::string pedestalJsonPath = + "/home/apareti/HidraSim2025/TBDataPreparation/2025_SPS/MapAndCalibration/SiPM_pedestals_v1.json"; + + const std::string adcToGeVJsonPath = + "/home/apareti/HidraSim2025/TBDataPreparation/2025_SPS/MapAndCalibration/SiPM_ADCtoGeV_v1.json"; + + // Prepare lookup map for SiPMs, with attached pedestal and noise info + // key: "tower|type|row|column", value: map entry with FERS/channel info and pedestal/noise SipmLookup sipmLookup; try { sipmLookup = LoadSipmMap(sipmMapPath); - std::cout << "Loaded " << sipmLookup.size() << " SiPM map entries from " - << sipmMapPath << std::endl; + + const PedestalHgVector pedestalsHg = LoadHgPedestalsFromJson(pedestalJsonPath); + const AdcToGeVVector adcToGeV = LoadAdcToGeVFromJson(adcToGeVJsonPath); + + AttachHgPedestalNoiseToSipmMap(sipmLookup, pedestalsHg, adcToGeV); + + std::cout << "Loaded " << sipmLookup.size() + << " SiPM map entries from " << sipmMapPath << std::endl; + std::cout << "Loaded " << pedestalsHg.size() + << " HG pedestal entries from " << pedestalJsonPath << std::endl; + std::cout << "Loaded " << adcToGeV.size() + << " ADC->GeV entries from " << adcToGeVJsonPath << std::endl; } catch (const std::exception& e) { std::cerr << e.what() << std::endl; return; } + // read input file and prepare output file const std::string infile = "../build/" + input; std::cout << "Using file: " << infile << std::endl; @@ -772,32 +1139,35 @@ void HidraTB25Ana(double energy, const std::string& input) TDirectory* eventDisplayDir = f.mkdir("EventDisplays"); f.cd(); - TRandom3 rng(NoiseRandomSeed); + TTree* activatedChannelsTree = + new TTree("ActivatedChannels", "Per-event activated channel output"); - const std::string sciNoiseCsvPath = "S_sipm_std_run868_HG_4fers.csv"; - const std::string cerNoiseCsvPath = "C_sipm_std_run868_HG_4fers.csv"; + unsigned int outEvent = 0; + std::vector activatedChannelOutput(TotalCalibChannels, 0.0); + std::vector activatedChannelOutputTrue(TotalCalibChannels, 0.0); + std::vector activatedChannelOutputSmeared(TotalCalibChannels, 0.0); - ChannelNoiseVector sciChannelNoiseValues; - ChannelNoiseVector cerChannelNoiseValues; + activatedChannelsTree->Branch("event", &outEvent); + activatedChannelsTree->Branch("activatedChannelOutput", &activatedChannelOutput); + activatedChannelsTree->Branch("activatedChannelOutputTrue", &activatedChannelOutputTrue); + activatedChannelsTree->Branch("activatedChannelOutputSmeared", &activatedChannelOutputSmeared); - try { - sciChannelNoiseValues = LoadNoiseValuesFromCsv(sciNoiseCsvPath); - cerChannelNoiseValues = LoadNoiseValuesFromCsv(cerNoiseCsvPath); + // --------------------------------------- + // create the RNG banks once + // Prepare random generators for noise smearing: + // one per calibration channel for uncorrelated noise, one per FERS for correlated noise + // ---------------------------------------- + NoiseGeneratorBank channelNoiseGenerators = BuildChannelNoiseGenerators(); + NoiseGeneratorBank fersNoiseGenerators = BuildFersNoiseGenerators(); - std::cout << "Loaded " << sciChannelNoiseValues.size() - << " S-channel noise values from " << sciNoiseCsvPath << std::endl; - std::cout << "Loaded " << cerChannelNoiseValues.size() - << " C-channel noise values from " << cerNoiseCsvPath << std::endl; - } catch (const std::exception& e) { - std::cerr << e.what() << std::endl; - return; - } - int modcol[NofModulesX * NofModulesY]; - int modrow[NofModulesX * NofModulesY]; - for (int i = 0; i < NofModulesX * NofModulesY; i++) { - int row = i / NofModulesX; - int col = i % NofModulesX; + + // Load module mapping to fill the event display maps + int modcol[NofmodulesX * NofmodulesY]; + int modrow[NofmodulesX * NofmodulesY]; + for (int i = 0; i < NofmodulesX * NofmodulesY; i++) { + int row = i / NofmodulesX; + int col = i % NofmodulesX; int imod = modflag[i]; if (imod >= 0) { modcol[imod] = col; @@ -805,6 +1175,9 @@ void HidraTB25Ana(double energy, const std::string& input) } } + // ---------------------------------------- + // Prepare output histograms + // ---------------------------------------- const double bmin = energy - 0.4 * std::sqrt(energy) * 10.; const double bmax = energy + 0.4 * std::sqrt(energy) * 10.; @@ -816,8 +1189,8 @@ void HidraTB25Ana(double energy, const std::string& input) auto leakene = new TH1F("leakene", "leakene", 100, 0., 0.1); auto chidist = new TH1F("chidist", "chidist", 100, 0., 1.); auto mapcalo = new TH2F("mapcalo", "mapcalo", - NofModulesX, 0., NofModulesX, - NofModulesY, 0., NofModulesY); + NofmodulesX, 0., NofmodulesX, + NofmodulesY, 0., NofmodulesY); auto SipmMapS = new TH2F("SipmMapS", "SipmS; Col; Row", NofSiPMTowersX * NofFiberscolumn, 0, NofSiPMTowersX * NofFiberscolumn, NofSiPMTowersY * NofFibersrow, 0, NofSiPMTowersY * NofFibersrow); @@ -848,8 +1221,15 @@ void HidraTB25Ana(double energy, const std::string& input) auto SignalSfersID = new TH1F("SignalSfersID", "Signal S fibers FERS ID; FERS; Signal", 16, 0.5, 16.5); auto SignalCfersID = new TH1F("SignalCfersID", "Signal C fibers FERS ID; FERS; Signal", 16, 0.5, 16.5); - auto scieneFersOn = new TH1F("scieneFersOn", "S energy for Activated FERS", 100, -50., bmax); - auto cereneFersOn = new TH1F("cereneFersOn", "C energy for Activated FERS", 100, -50., bmax); + auto scieneFersOn = new TH1F("scieneFersOn", "S energy for Activated FERS", 100, 0., bmax); + auto cereneFersOn = new TH1F("cereneFersOn", "C energy for Activated FERS", 100, 0., bmax); + + auto h2_SvsC = new TH2F("SvsC", "S vs C energy (all fibres); S/E; C/E", 100, -0.5, 1.5, 100, 0., 1.5); + auto h2_SvsCFersOn = new TH2F("SvsCFersOn", "S vs C energy for Activated FERS; S/E; C/E", 100, 0., 1.5, 100, 0., 1.5); + + // Count number of activated channels per FERS, to check the multiplicity distribution and the effect of different thresholds + auto ActivatedChannelsPerFers = new TH2F("ActivatedChannelsPerFers", "Activated channels per FERS; FERS ID; Channels above threshold",16, 0.5, 16.5, 65, -0.5, 64.5); + const int nentries = simtree->GetEntries(); std::cout << "Entries " << nentries << std::endl; @@ -873,6 +1253,7 @@ void HidraTB25Ana(double energy, const std::string& input) for (unsigned int i = 0; i < static_cast(nentries); i++) { simtree->GetEntry(i); + // store one event display every N events const bool writeEventDisplay = (EventDisplayEvery > 0) && (((i + 1) % EventDisplayEvery) == 0); @@ -886,7 +1267,13 @@ void HidraTB25Ana(double energy, const std::string& input) double barX_cer = 0.; double barY_cer = 0.; + outEvent = i; + // prepare the per-event data structures for the smearing and activation steps + std::fill(activatedChannelOutput.begin(), activatedChannelOutput.end(), 0.0); + std::fill(activatedChannelOutputTrue.begin(), activatedChannelOutputTrue.end(), 0.0); + std::fill(activatedChannelOutputSmeared.begin(), activatedChannelOutputSmeared.end(), 0.0); std::unordered_map fersChannelsOverThreshold; + std::vector eventDisplaySciHits; std::vector eventDisplayCerHits; @@ -901,6 +1288,7 @@ void HidraTB25Ana(double energy, const std::string& input) eventDisplayCerHits.reserve(CSiPM->size()); } + // Fill energy in PMT towers for (unsigned int j = 0; j < SPMT->size(); j++) { totsci += SPMT->at(j) / sciPheGeV; totcer += CPMT->at(j) / cerPheGeV; @@ -911,47 +1299,62 @@ void HidraTB25Ana(double energy, const std::string& input) double sciPosWeight = 0.0; double cerPosWeight = 0.0; + // Prepare shared random map to have correlated noise between S/C channels + std::unordered_map fersEventNoiseCache; + - // Start loop over SiPM signals - S channels + // Start looping on SiPM signals for (unsigned int n = 0; n < SSiPM->size(); n++) { - const double rawPhe = SSiPM->at(n); - const double content = rawPhe / sciPheGeV; + const double rawPhe = SSiPM->at(n); // number of photoelectrons + const double content = rawPhe / sciPheGeV; // converted to GeV totsci += content; const unsigned int towID = static_cast(n / (NofFiberscolumn * NofFibersrow / 2)); const unsigned int SiPMID = n % (NofFiberscolumn * NofFibersrow / 2); - const unsigned int colID = NofFiberscolumn - static_cast(SiPMID / (NofFibersrow / 2)); + //const unsigned int colID = NofFiberscolumn - static_cast(SiPMID / (NofFibersrow / 2)); + const unsigned int colID = NofFiberscolumn -1U - static_cast(SiPMID / (NofFibersrow / 2)); const unsigned int rowID = 2 * static_cast(SiPMID % (NofFibersrow / 2)); double SiPM_X = 0.0; double SiPM_Y = 0.0; - GetSiPMcoordinate(towID, rowID, colID, SiPM_X, SiPM_Y, "S", grouping); // Not used after JSON SiPM map implementation, kept for reference + GetSiPMcoordinate(towID, rowID, colID, SiPM_X, SiPM_Y, "S", grouping); // Currently not used, JSON map is read instead - const std::string tbTower = simTower_to_tbTower(std::to_string(towID)); - const int groupedCol = grouped_column(static_cast(colID), grouping); + const std::string tbTower = simTower_to_tbTower(std::to_string(towID)); // convert from sim tower ID to test beam tower ID + const int groupedCol = grouped_column(static_cast(colID), grouping); // single fibre to channel mapping SipmMapS->Fill(colID, towID * NofFibersrow + rowID, content); - const SipmMapEntry* info = FindSipmInfo(sipmLookup, tbTower, "S", static_cast(rowID), groupedCol, grouping); + // now associate the SiPM signal to the corresponding FERS and channel + const SipmMapEntry* info = + FindSipmInfo(sipmLookup, tbTower, "S", static_cast(rowID), groupedCol, grouping); if (info) { - SciSiPMCoordinates->Fill(info->x, info->y, content); + SciSiPMCoordinates->Fill(info->x, info->y, content); // add fibre content to mapped channel signal SignalSfersID->Fill(info->fersId, content / nentries); barX_sci += info->x * content; barY_sci += info->y * content; sciPosWeight += content; + // accumulate the signal in the corresponding channel for smearing and activation steps if (info->boardID >= 0 && info->fersId >= 0) { - const FersKey key = make_fers_key(info->boardID, info->fersId); - const double outputSignal = ApplySaturation ? apply_sipm_saturation_from_pe(rawPhe, sciPheGeV, NofSipmCells_sci) : content; + // create the key to identify the FERS this SiPM belongs to + const FersKey key = make_fers_key(info->boardID, info->fersId); + + // Apply saturation consistently for both threshold and output energy + const double saturatedSignal = + ApplySaturation + ? apply_sipm_saturation_from_pe(rawPhe, sciPheGeV, NofSipmCells_sci) + : content; + // accumulate the signal for this channel + // reads threshold values for this channel from the map + // and updates the signal and FERS->channel count maps AccumulateChannelSignal( sciChannelSignals, *info, - content, - outputSignal + saturatedSignal ); if (writeEventDisplay) { @@ -959,11 +1362,8 @@ void HidraTB25Ana(double energy, const std::string& input) } } } - } // end loop over S SiPM channels - - + } - // Start loop over SiPM signals - C channels for (unsigned int n = 0; n < CSiPM->size(); n++) { const double rawPhe = CSiPM->at(n); const double content = rawPhe / cerPheGeV; @@ -972,10 +1372,11 @@ void HidraTB25Ana(double energy, const std::string& input) const unsigned int towID = static_cast(n / (NofFiberscolumn * NofFibersrow / 2)); const unsigned int SiPMID = n % (NofFiberscolumn * NofFibersrow / 2); - const unsigned int colID = NofFiberscolumn - static_cast(SiPMID / (NofFibersrow / 2)); + //const unsigned int colID = NofFiberscolumn - static_cast(SiPMID / (NofFibersrow / 2)); + const unsigned int colID = NofFiberscolumn -1U - static_cast(SiPMID / (NofFibersrow / 2)); const unsigned int rowID = 2 * static_cast(SiPMID % (NofFibersrow / 2)) + 1; - SipmMapC->Fill(colID, towID * NofFibersrow + rowID, content); // Not used after JSON SiPM map implementation, kept for reference + SipmMapC->Fill(colID, towID * NofFibersrow + rowID, content); const std::string tbTower = simTower_to_tbTower(std::to_string(towID)); const int groupedCol = grouped_column(static_cast(colID), grouping); @@ -990,9 +1391,12 @@ void HidraTB25Ana(double energy, const std::string& input) barY_cer += info->y * content; cerPosWeight += content; + // accumulate the signal for this channel if (info->boardID >= 0 && info->fersId >= 0) { const FersKey key = make_fers_key(info->boardID, info->fersId); - const double outputSignal = + + // Apply saturation consistently for both threshold and output energy + const double saturatedSignal = ApplySaturation ? apply_sipm_saturation_from_pe(rawPhe, cerPheGeV, NofSipmCells_cer) : content; @@ -1000,8 +1404,7 @@ void HidraTB25Ana(double energy, const std::string& input) AccumulateChannelSignal( cerChannelSignals, *info, - content, - outputSignal + saturatedSignal ); if (writeEventDisplay) { @@ -1009,41 +1412,74 @@ void HidraTB25Ana(double energy, const std::string& input) } } } - } // end loop over C SiPM channels + } // end C Sipm loop - SmearAndCountChannelsOverThreshold( + + //------------------------------ + // FERS ACTIVATION LOGIC + //------------------------------ + + // Now apply noise smearing and thresholding to the accumulated channel signals + // and fill the activated channel output array for this event + // Count number of channels over threshold per FERS + SmearAndCountChannelsOverThreshold( // S channel sciChannelSignals, fersChannelsOverThreshold, - rng, - sciChannelNoiseValues + NoiseMode, + channelNoiseGenerators, + fersNoiseGenerators, + fersEventNoiseCache ); - SmearAndCountChannelsOverThreshold( + SmearAndCountChannelsOverThreshold( // C channel cerChannelSignals, fersChannelsOverThreshold, - rng, - cerChannelNoiseValues + NoiseMode, + channelNoiseGenerators, + fersNoiseGenerators, + fersEventNoiseCache ); + for (unsigned int fersId = 1; fersId <= TotalFers; ++fersId) { + unsigned int nOverThreshold = 0; + + for (const auto& [fersKey, count] : fersChannelsOverThreshold) { + if (static_cast(static_cast(fersKey)) == fersId) { + nOverThreshold += count; + } + } + + ActivatedChannelsPerFers->Fill(fersId, nOverThreshold); + } + + // For activated FERS channels, subtract the pedestal median + SubtractPedestalMedianFromActivatedChannels( + sciChannelSignals, + fersChannelsOverThreshold + ); + SubtractPedestalMedianFromActivatedChannels( + cerChannelSignals, + fersChannelsOverThreshold + ); - // Apply 1-pe loss only after FERS activation is known. - ApplyPoissonOnePeLossToActivatedChannels( + FillActivatedChannelOutputArrays( sciChannelSignals, fersChannelsOverThreshold, - rng, - sciPheGeV, - SciOnePeLossLambda + activatedChannelOutput, + activatedChannelOutputTrue, + activatedChannelOutputSmeared ); - ApplyPoissonOnePeLossToActivatedChannels( + FillActivatedChannelOutputArrays( cerChannelSignals, fersChannelsOverThreshold, - rng, - cerPheGeV, - CerOnePeLossLambda + activatedChannelOutput, + activatedChannelOutputTrue, + activatedChannelOutputSmeared ); + activatedChannelsTree->Fill(); PrintPerEventFersChannelCounts( @@ -1112,6 +1548,8 @@ void HidraTB25Ana(double energy, const std::string& input) f.cd(); } + + scieneFersOn->Fill(totsciFersOn); cereneFersOn->Fill(totcerFersOn); @@ -1122,6 +1560,13 @@ void HidraTB25Ana(double energy, const std::string& input) totdep->Fill(tottow / 1000.); leakene->Fill(lenergy / 1000. / energy); chidist->Fill((totsci - ecalo) / (totcer - ecalo)); + + + h2_SvsC->Fill(totsci / energy, totcer / energy); + h2_SvsCFersOn->Fill(totsciFersOn / energy, totcerFersOn / energy); + + + //break; } const GaussianFitSummary sFit = FitGaussianSummary(sciene); @@ -1161,4 +1606,4 @@ void HidraTB25Ana(double energy, const std::string& input) simfile->Close(); delete simfile; -} \ No newline at end of file +} diff --git a/include/DREMTubesGeoPar.hh b/include/DREMTubesGeoPar.hh index 25a6b95..cb374d4 100644 --- a/include/DREMTubesGeoPar.hh +++ b/include/DREMTubesGeoPar.hh @@ -9,7 +9,7 @@ const G4int NofLeakCounterLayers = 4; const G4bool PreShowerIn = false; const G4bool LeakageCounterIn = true; -const G4bool TruthLeakageIn = false; +const G4bool TruthLeakageIn = true; @@ -78,7 +78,7 @@ const G4int modflag[90]={-1, -1, 0, -1, -1, -1, 1, 2, 3, -1, -1, 4, 5, 6, -1, -1, 7, 8, 9, -1, - 10, 1, 12, 13, 14, + 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, @@ -96,12 +96,50 @@ const G4int NoModulesSiPM=8; const G4int SiPMMod[8]={17,22,27,32,37,42,47,52}; const G4int NofFiberscolumn = 64; const G4int NofFibersrow = 16; -const G4int NoModulesActive=80; +const G4int NoModulesActive=70; const G4double moduleZ = (2500.)*mm; const G4bool irot=false; const G4int NoFibersTower=NofFiberscolumn*NofFibersrow/2; +/* +// 80 modules -> TB26 (HiDRa geometry, SiPMs placed in two towers) +const G4int NofmodulesX = 5; +const G4int NofmodulesY = 20; +const G4int modflag[100]={-1, -1, 0, -1, -1, + -1, 1, 2, 3, -1, + -1, 4, 5, 6, -1, + -1, 7, 8, 9, -1, + 10, 11, 12, 13, 14, + 15, 16, 17, 18, 19, + 20, 21, 22, 23, 24, + 25, 26, 27, 28, 29, + 30, 31, 32, 33, 34, + 35, 36, 37, 38, 39, + 40, 41, 42, 43, 44, + 45, 46, 47, 48, 49, + 50, 51, 52, 53, 54, + 55, 56, 57, 58, 59, + 60, 61, 62, 63, 64, + 65, 66, 67, 68, 69, + -1, 70, 71, 72, -1, + -1, 73, 74, 75, -1, + -1, 76, 77, 78, -1, + -1, -1, 79, -1, -1}; +const G4int NoModulesSiPM=10; +//const G4int SiPMMod[10]={27,32,37,42,47,52,53, 48, 43, 38}; +const G4int SiPMMod[10]={27,32,37,42,47,52, 26, 31, 36, 41}; +const G4int NofFiberscolumn = 64; +const G4int NofFibersrow = 16; +const G4int NoModulesActive=80; +const G4double moduleZ = (2500.)*mm; +const G4bool irot=false; +const G4int NoFibersTower=NofFiberscolumn*NofFibersrow/2;*/ + + + + + /* // 80 modules -> Hidra Final geometry diff --git a/src/DREMTubesDetectorConstruction.cc b/src/DREMTubesDetectorConstruction.cc index f2985c7..95d7225 100644 --- a/src/DREMTubesDetectorConstruction.cc +++ b/src/DREMTubesDetectorConstruction.cc @@ -689,12 +689,14 @@ G4VPhysicalVolume* DREMTubesDetectorConstruction::DefineVolumes() { // Calorimeter placement (with rotation wrt beam axis) // - G4RotationMatrix rotm = G4RotationMatrix(); // Negative sign to make angles same as TB platform G4double xrot = - fVerrot; - G4double yrot = + fOrzrot; - rotm.rotateX(xrot); - rotm.rotateY(yrot); + G4double yrot = - fOrzrot; + G4RotationMatrix rotX = G4RotationMatrix(); + G4RotationMatrix rotY = G4RotationMatrix(); + rotX.rotateX(xrot); + rotY.rotateY(yrot); + G4RotationMatrix rotm = rotY * rotX; G4ThreeVector position; @@ -708,8 +710,6 @@ G4VPhysicalVolume* DREMTubesDetectorConstruction::DefineVolumes() { // // Build closed tube for detailed leakage study // - // comment for later introducing boolean variable - // to include truth leakage counters (or not) if(TruthLeakageIn) { G4double leakradint=sqrt( (caloBoxX)*(caloBoxX)+(caloBoxY)*(caloBoxY)); // Added *1.2 wrt Giacomo's @@ -725,10 +725,10 @@ G4VPhysicalVolume* DREMTubesDetectorConstruction::DefineVolumes() { "leakageabsorberl"); G4VisAttributes* LkVisAttl = new G4VisAttributes(G4Colour(0.0,0.8,0.0)); //green LkVisAttl->SetVisibility(true); - LkVisAttl->SetForceWireframe(true); + //LkVisAttl->SetForceWireframe(true); LkVisAttl->SetForceSolid(true); - //leakageabsorberlLV->SetVisAttributes(LkVisAttl); - leakageabsorberlLV->SetVisAttributes(invisibleAttr); + leakageabsorberlLV->SetVisAttributes(LkVisAttl); + //leakageabsorberlLV->SetVisAttributes(invisibleAttr); //leakageabsorberlLV->SetVisAttributes(G4VisAttributes::Invisible); new G4PVPlacement( transform, @@ -747,11 +747,12 @@ G4VPhysicalVolume* DREMTubesDetectorConstruction::DefineVolumes() { "leakageabsorberd"); G4VisAttributes* LkVisAttd = new G4VisAttributes(G4Colour(0.0,0.4,0.0)); //green LkVisAttd->SetVisibility(true); - LkVisAttd->SetForceWireframe(true); + //LkVisAttd->SetForceWireframe(true); LkVisAttd->SetForceSolid(true); - //leakageabsorberdLV->SetVisAttributes(LkVisAttd); - leakageabsorberdLV->SetVisAttributes(invisibleAttr); - //leakageabsorberdLV->SetVisAttributes(G4VisAttributes::Invisible); + leakageabsorberdLV->SetVisAttributes(LkVisAttd); + + + //leakageabsorberdLV->SetVisAttributes(invisibleAttr); G4ThreeVector positiond; positiond.setX(caloBoxZ*sin(xrot)); positiond.setY(-caloBoxZ*sin(yrot)); @@ -1193,8 +1194,7 @@ G4LogicalVolume* DREMTubesDetectorConstruction::constructscinfiber(double tolera G4LogicalVolume* logic_S_fiber = new G4LogicalVolume(S_fiber, absorberMaterial, "S_fiber"); - //logic_S_fiber->SetVisAttributes(G4VisAttributes::Invisible); - logic_S_fiber->SetVisAttributes(invisibleAttr); + //logic_S_fiber->SetVisAttributes(invisibleAttr); G4Tubs* Abs_S_fiber = new G4Tubs("Abs_Scin_fiber", claddingradiusmax, tuberadius, fiberZ/2,0.,2.*pi); @@ -1222,9 +1222,8 @@ G4LogicalVolume* DREMTubesDetectorConstruction::constructscinfiber(double tolera ScincoreVisAtt->SetVisibility(true); ScincoreVisAtt->SetForceWireframe(true); ScincoreVisAtt->SetForceSolid(true); - //logic_Core_S_fiber->SetVisAttributes(ScincoreVisAtt); - //logic_Core_S_fiber->SetVisAttributes(G4VisAttributes::Invisible); - logic_Core_S_fiber->SetVisAttributes(invisibleAttr); + logic_Core_S_fiber->SetVisAttributes(ScincoreVisAtt); + //logic_Core_S_fiber->SetVisAttributes(invisibleAttr); G4ThreeVector vec_Core_S; vec_Core_S.setX(0.); @@ -1253,8 +1252,7 @@ G4LogicalVolume* DREMTubesDetectorConstruction::constructscinfiber(double tolera ScincladVisAtt->SetVisibility(true); ScincladVisAtt->SetForceWireframe(true); ScincladVisAtt->SetForceSolid(true); - logic_Clad_S_fiber->SetVisAttributes(ScincladVisAtt); - //logic_Clad_S_fiber->SetVisAttributes(G4VisAttributes::Invisible); + //logic_Clad_S_fiber->SetVisAttributes(ScincladVisAtt); logic_Clad_S_fiber->SetVisAttributes(invisibleAttr); @@ -1278,7 +1276,7 @@ G4LogicalVolume* DREMTubesDetectorConstruction::constructscinfiber(double tolera TubeVisAtt->SetForceWireframe(true); TubeVisAtt->SetForceSolid(true); logic_Abs_S_fiber->SetVisAttributes(TubeVisAtt); - //logic_Abs_S_fiber->SetVisAttributes(G4VisAttributes::Invisible); + //logic_Abs_S_fiber->SetVisAttributes(invisibleAttr); return logic_S_fiber; @@ -1308,8 +1306,7 @@ G4LogicalVolume* DREMTubesDetectorConstruction::constructcherfiber(double tolera absorberMaterial, "C_fiber"); - //logic_C_fiber->SetVisAttributes(G4VisAttributes::Invisible); - logic_C_fiber->SetVisAttributes(invisibleAttr); + //logic_C_fiber->SetVisAttributes(invisibleAttr); G4Tubs* Abs_C_fiber = new G4Tubs("Abs_Cher_fiber", claddingradiusmax, tuberadius, fiberZ/2,0.,2.*pi); @@ -1336,8 +1333,7 @@ G4LogicalVolume* DREMTubesDetectorConstruction::constructcherfiber(double tolera ChercoreVisAtt->SetForceWireframe(true); ChercoreVisAtt->SetForceSolid(true); logic_Core_C_fiber->SetVisAttributes(ChercoreVisAtt); - //logic_Core_C_fiber->SetVisAttributes(G4VisAttributes::Invisible); - logic_Core_C_fiber->SetVisAttributes(invisibleAttr); + //logic_Core_C_fiber->SetVisAttributes(invisibleAttr); G4ThreeVector vec_Core_C; vec_Core_C.setX(0.); @@ -1388,7 +1384,7 @@ G4LogicalVolume* DREMTubesDetectorConstruction::constructcherfiber(double tolera TubeVisAtt->SetForceWireframe(true); TubeVisAtt->SetForceSolid(true); logic_Abs_C_fiber->SetVisAttributes(TubeVisAtt); - //logic_Abs_C_fiber->SetVisAttributes(G4VisAttributes::Invisible); + //logic_Abs_C_fiber->SetVisAttributes(invisibleAttr); return logic_C_fiber; From 5d793d12fd3ba0d07eea8124ec2c3b9510e91393 Mon Sep 17 00:00:00 2001 From: Andrea Pareti Date: Sat, 16 May 2026 11:46:01 +0200 Subject: [PATCH 5/5] Add TB25 analysis files. --- analysis/HidraGeo.h | 37 - analysis/README.md | 53 + .../C_sipm_std_run868_HG_4fers.csv | 512 ++++++ analysis/tb25_analysis/HidraTB25Ana.C | 1623 +++++++++++++++++ analysis/tb25_analysis/README.md | 81 + .../S_sipm_std_run868_HG_4fers.csv | 512 ++++++ analysis/tb25_analysis/plot_hidra.py | 529 ++++++ analysis/tb25_analysis/run_ana.sh | 56 + 8 files changed, 3366 insertions(+), 37 deletions(-) create mode 100644 analysis/README.md create mode 100644 analysis/tb25_analysis/C_sipm_std_run868_HG_4fers.csv create mode 100644 analysis/tb25_analysis/HidraTB25Ana.C create mode 100644 analysis/tb25_analysis/README.md create mode 100644 analysis/tb25_analysis/S_sipm_std_run868_HG_4fers.csv create mode 100644 analysis/tb25_analysis/plot_hidra.py create mode 100755 analysis/tb25_analysis/run_ana.sh diff --git a/analysis/HidraGeo.h b/analysis/HidraGeo.h index fbbe50d..995622a 100644 --- a/analysis/HidraGeo.h +++ b/analysis/HidraGeo.h @@ -14,42 +14,6 @@ const double tuberadius = 1.0; -/* -// 70 modules -> TB25 -const int NofModulesX = 5; -const int NofModulesY = 18; -// 80 miniM -> But different geometry-> move one minimodules from row to corner, horizontal placement -const int modflag[90]={-1, -1, 0, -1, -1, - -1, 1, 2, 3, -1, - -1, 4, 5, 6, -1, - -1, 7, 8, 9, -1, - 10, 1, 12, 13, 14, - 15, 16, 17, 18, 19, - 20, 21, 22, 23, 24, - 25, 26, 27, 28, 29, - 30, 31, 32, 33, 34, - 35, 36, 37, 38, 39, - 40, 41, 42, 43, 44, - 45, 46, 47, 48, 49, - 50, 51, 52, 53, 54, - 55, 56, 57, 58, 59, - -1, 60, 61, 62, -1, - -1, 63, 64, 65, -1, - -1, 66, 67, 68, -1, - -1, -1, 69, -1, -1}; -const int NoModulesSiPM=8; -const int SiPMMod[8]={17,22,27,32,37,42,47,52}; -const int NofSiPMTowersX = 1; -const int NofSiPMTowersY = 8; -const int NofFiberscolumn = 64; -const int NofFibersrow = 16; -const int NoModulesActive=80; -const double moduleZ = (2500.); -const bool irot=false; -const int NoFibersTower=NofFiberscolumn*NofFibersrow/2;*/ - - - @@ -116,7 +80,6 @@ const int modflag[100]={-1, -1, 0, -1, -1, -1, -1, 79, -1, -1}; const int NoModulesSiPM=10; const int SiPMMod[10]={27,32,37,42,47,52, 26, 31, 36, 41}; - const int NofFiberscolumn = 64; const int NofFibersrow = 16; const int NofSiPMTowersX = 2; diff --git a/analysis/README.md b/analysis/README.md new file mode 100644 index 0000000..3adb5c6 --- /dev/null +++ b/analysis/README.md @@ -0,0 +1,53 @@ +# HiDRa Analysis + +This directory contains the shared geometry header and the standard ROOT macro +used to analyze HidraSim simulation ntuples. + +## Standard Simulation Analysis + +`HidraAna.C` is the general analysis macro for HidraSim simulation ntuples. +It reads the simulation output tree, reconstructs the calorimeter response, and +produces ROOT histograms for the requested beam energy. + +The workflow is: + +1. Open the input ROOT file and read the `DREMTubes` simulation tree. +2. Load the detector geometry from `HidraGeo.h`. +3. Loop over events and accumulate PMT tower signals and SiPM fiber signals. +4. Convert scintillation and Cherenkov SiPM photoelectron yields to energy. +5. Reconstruct SiPM channel positions from the simulation vector index. +6. Fill calorimeter maps, SiPM coordinate maps, energy histograms, and optional + event-display histograms. + +Inputs: + +- Beam energy, passed as the first macro argument. +- HidraSim ROOT ntuple, passed as the second macro argument. +- Geometry constants from `HidraGeo.h`. + +Outputs: + +- A ROOT histogram file named from the input energy, for example `hidra10.root`. +- Energy-response histograms for scintillation, Cherenkov, and combined signals. +- Calorimeter tower maps and SiPM coordinate maps. +- Optional event-display histograms written at the configured event interval. + +Typical usage from this directory: + +```bash +root -l -b -q 'HidraAna.C(energy, "input.root")' +``` + +The macro uses `HidraGeo.h` for the detector geometry constants. SiPM channels +are decoded from the simulation vector index using the same convention used when +the SiPM copy numbers are assigned in the Geant4 geometry. The resulting SiPM +coordinates are filled in the simulation coordinate system. + +Generated ROOT files, logs, summaries, compiled ROOT dictionaries, and plots are +analysis outputs and should normally not be committed. + +## TB25 Analysis + +TB25-specific analysis files live in `tb25_analysis/`. See +`tb25_analysis/README.md` for the TB25 workflow, mapping inputs, and configurable +analysis switches. diff --git a/analysis/tb25_analysis/C_sipm_std_run868_HG_4fers.csv b/analysis/tb25_analysis/C_sipm_std_run868_HG_4fers.csv new file mode 100644 index 0000000..e947196 --- /dev/null +++ b/analysis/tb25_analysis/C_sipm_std_run868_HG_4fers.csv @@ -0,0 +1,512 @@ +5.724510147096531498e-02 +6.626481881954096598e-02 +4.219266674907513892e-02 +5.506907649265719434e-02 +3.828036407840942207e-02 +3.325525206966294806e-02 +4.082193006097444010e-02 +6.550448163148787917e-02 +3.351868509432122389e-02 +3.777186892071342228e-02 +3.950095309251594466e-02 +4.750018269345222810e-02 +2.927544623466722673e-02 +3.918231467835065623e-02 +2.935932839417737181e-02 +4.812432938103271274e-02 +5.197307691505153454e-02 +5.503661834623573695e-02 +6.169490779247599477e-02 +6.527661140581030030e-02 +5.673161406795954737e-02 +4.275658237245028026e-02 +5.077181188265304790e-02 +7.262914259140776374e-02 +3.957163554787955506e-02 +2.192126562879497670e-02 +3.634946642524859717e-02 +3.722150116307548773e-02 +3.884236125127533967e-02 +3.850775848833330739e-02 +3.687327807498621407e-02 +1.967238078286064756e-02 +5.455986133496494583e-02 +5.453699567458629438e-02 +5.421876021515959626e-02 +5.269804681269302238e-02 +5.176939442094376231e-02 +5.248959903370935265e-02 +5.199133299278355103e-02 +4.955606533351627346e-02 +6.366525881301898782e-02 +4.735075862419670101e-02 +4.344455288088732575e-02 +6.765026749847224963e-02 +4.041414351815468542e-02 +6.612504710715799705e-02 +3.484061551300083465e-02 +4.567376507873280311e-02 +3.917559955857029957e-02 +3.084431485486385099e-02 +4.161250258399920576e-02 +4.140564825953479272e-02 +4.604652886085825597e-02 +3.657446085851994572e-02 +3.948420239267604048e-02 +3.392497594174766901e-02 +2.497390225862996735e-02 +2.252307155325400725e-02 +2.953926353172835889e-02 +2.957793492988061018e-02 +2.584566708553267048e-02 +2.135107003527773536e-02 +2.676080714936554528e-02 +2.816250759466798828e-02 +3.840676223511878262e-02 +2.825119247307884829e-02 +3.260046666038951546e-02 +3.361184401064538918e-02 +3.495451056560599373e-02 +3.449584222260310140e-02 +3.340828531977423771e-02 +4.819817769256029560e-02 +2.687605649414885611e-02 +2.914960105904445056e-02 +2.382188605593601322e-02 +2.196424650053465516e-02 +2.544843172959233152e-02 +2.489044023104545686e-02 +2.765693640125998706e-02 +3.893792139355987253e-02 +3.309090214499994104e-02 +4.570901987336374594e-02 +3.864267265323068773e-02 +4.553617016225698644e-02 +3.789906916018015437e-02 +3.376308207956683904e-02 +3.436962793932484689e-02 +4.403998564819255196e-02 +2.924144954536263968e-02 +2.953465318640857706e-02 +2.957507931862856168e-02 +2.866946419461026754e-02 +2.868352462993253860e-02 +2.775919779905623025e-02 +2.891570578029308350e-02 +2.310679051008434140e-02 +4.110266549299177591e-02 +3.994999639863443253e-02 +4.124992736214831701e-02 +1.587423288062141391e-01 +3.883536902696321097e-02 +3.130048592316463096e-02 +3.920810200080476066e-02 +3.998971800670938354e-02 +4.625241374628175778e-02 +4.445138901924772235e-02 +3.620208343244758070e-02 +3.632772345809844733e-02 +3.671747556434799831e-02 +3.845764747047757115e-02 +6.145077643050742694e-02 +6.275794857747662991e-02 +2.710450300058916615e-02 +3.065827940918849220e-02 +2.536134713180719299e-02 +2.609161710875338228e-02 +3.360501315729845506e-02 +2.513741550056875287e-02 +2.826820743903556143e-02 +3.906692318720579593e-02 +2.680673785641322221e-02 +2.580401352856335692e-02 +1.893414107117860970e-02 +2.221255482078032112e-02 +2.258584521393766295e-02 +2.334000482788323519e-02 +2.460353891712580704e-02 +3.182025817813655816e-02 +2.780335689706518384e-02 +2.850155868089976688e-02 +3.448930640387523061e-01 +3.872042441310070277e-02 +7.482972612259421175e-02 +5.953770467594093896e-02 +4.253641180749766215e-02 +9.154786050670177255e-02 +5.295876919698517399e-02 +2.909451822480694019e-02 +3.136401769353826252e-02 +3.030808474481942738e-02 +3.467710111865229816e-02 +4.655097125465134517e-02 +2.858261509313269125e-02 +3.478119031103214798e-02 +5.174373288531999848e-02 +4.619986204275049191e-02 +5.756246700871312655e-02 +3.090070256722287148e-02 +2.766248883398981714e-02 +3.470740033149637138e-02 +2.921072027466995596e-02 +3.671956361253017609e-02 +4.847421118111735344e-02 +5.060359576925324565e-02 +5.030162881708962613e-02 +5.097297653138799234e-02 +5.148656975574820993e-02 +5.392520998977159280e-02 +1.982846895695027437e-02 +4.850220646860814444e-02 +5.993692198636508400e-02 +6.236806902353336096e-02 +6.126754416677564208e-02 +6.109755252327509628e-02 +6.149938472607571760e-02 +6.442749975356647230e-02 +1.669416233165254768e-02 +5.900366206928353563e-02 +3.888738773793400799e-02 +7.833453394728602526e-02 +6.152518163335098794e-02 +5.552640810714965652e-02 +3.466281066131313843e-02 +6.471905471898907469e-02 +8.477333985380124115e-02 +4.027559093387316630e-02 +2.781580662698647510e-02 +2.716275949857578642e-02 +2.117938563068803848e-02 +2.546848095678794568e-02 +2.218224828664567411e-02 +2.408474158392464454e-02 +1.752127770161348783e-02 +2.900937462036274825e-02 +2.558424203980210615e-02 +2.667613875327943390e-02 +2.285644816781575231e-02 +3.444787877105025958e-02 +2.842969367144984952e-02 +2.781806309997611470e-02 +2.278747128677800424e-02 +2.273800364621981346e-02 +3.292241114871598956e-02 +4.093571200855365688e-02 +3.109280884970710979e-02 +3.056388062361755503e-02 +3.810492097436095349e-02 +5.530068398851914480e-02 +5.820703624452350233e-02 +6.747047243447470977e-02 +2.792753288729967268e-02 +3.541408227176504747e-02 +3.391639642090304996e-02 +2.375279936207745635e-02 +3.263363643401630959e-02 +3.798020168967643911e-02 +5.955960382315609453e-02 +3.424656448798864272e-02 +3.062339375136216568e-02 +2.595868238723082755e-02 +2.613576126660265009e-02 +1.968839166084699471e-02 +3.490278595207796375e-02 +4.272263345338596363e-02 +4.767457268524539354e-02 +4.549593408553304769e-02 +3.078137774123668527e-02 +3.257842107217198691e-02 +3.182607189668774927e-02 +3.232012287919943522e-02 +1.222717036686074045e-02 +3.433084125235238981e-02 +5.535088045983847793e-02 +5.566248950566251436e-02 +4.099238239669900780e-02 +4.027985464143709476e-02 +4.136480723571822615e-02 +4.063477032594655508e-02 +3.856713581820427755e-02 +3.997427441927237091e-02 +4.028391216756185833e-02 +4.045847654924513565e-02 +3.162392675681734849e-02 +4.676123999405545212e-02 +3.393514438022732921e-02 +2.793973124991238136e-02 +3.339246885513522112e-02 +4.045348921353596888e-02 +4.467152006517137453e-02 +4.661152077344749706e-02 +2.687799224442116719e-02 +2.237280095246063638e-02 +2.682206211049998484e-02 +2.458039734822748193e-02 +3.420416366071603320e-02 +3.084417065671561653e-02 +3.621234613946675007e-02 +2.541826400767688982e-02 +2.493502246859492441e-02 +2.557193102204225604e-02 +3.798420372928026922e-02 +6.787629505046859080e-02 +1.910916678702530291e-02 +3.008899343866349579e-02 +3.071585519915779813e-02 +2.995715399879163299e-02 +3.515818834380970287e-02 +3.101214685620806016e-02 +2.821956214113520781e-02 +3.143376370926833502e-02 +2.668317594594997108e-02 +2.417700882136975676e-02 +3.879502116499371211e-02 +3.691681493199174685e-02 +2.180097024707033232e-02 +2.026414859014154593e-02 +2.451970544616823011e-02 +2.503096109134547026e-02 +2.407318994384816113e-02 +2.500608560614900680e-02 +2.299322113522302841e-02 +1.892804334540312133e-02 +2.725836347263946302e-02 +2.685727256201346155e-02 +2.952309269055473223e-02 +2.593627948951328713e-02 +3.616302266905261492e-02 +3.134467149926647989e-02 +2.809790820998999544e-02 +2.421627309095179120e-02 +2.577070534151945994e-02 +2.763960765635891173e-02 +2.707892688516730026e-02 +2.758239805679607234e-02 +2.728592195636053708e-02 +2.679830856676626571e-02 +2.548074077734096068e-02 +2.261212822877539663e-02 +3.305548531147328717e-02 +3.307818544975543729e-02 +3.371539065836282045e-02 +3.281858345644109859e-02 +3.333752707418618372e-02 +3.353920405835733587e-02 +3.295948611821102048e-02 +3.078360009393782751e-02 +3.256214117851691930e-02 +2.987447180973953897e-02 +2.737944443558125787e-02 +2.830966316696533272e-02 +2.656029599119809348e-02 +2.621970821083659947e-02 +2.768836955786336401e-02 +2.683037691511970177e-02 +2.764467472827880287e-02 +2.269076984437770270e-02 +2.694729857899708439e-02 +2.012418259804852130e-02 +2.795227339971490324e-02 +2.093752932706900716e-02 +2.802158101751617628e-02 +3.278989953312570066e-02 +2.008038307797321148e-02 +1.744044472579010449e-02 +1.748139758668832214e-02 +2.249132534141337278e-02 +1.975647148884353624e-02 +1.549082445414413645e-01 +2.108635651297914235e-02 +1.981187057135514948e-02 +2.226239097489192756e-02 +2.301031811679525585e-02 +2.167962947027838455e-02 +2.132140713511673780e-02 +2.270707206846134599e-02 +2.236745312466420152e-02 +2.158779844282380025e-02 +4.918417591331687190e-03 +2.000148343947480076e-02 +1.840708980520614282e-02 +2.110507013553907310e-02 +2.357740836022443817e-02 +1.961735548133048540e-02 +2.093439001425866780e-02 +2.159572522720769674e-02 +2.360380493184898593e-02 +2.092547233022221814e-02 +1.871652173160147370e-02 +1.688097987517104409e-02 +1.861866982709256199e-02 +1.547999544062114360e-02 +1.486256425789687051e-02 +1.555089471850900663e-02 +1.817379751459282863e-02 +1.670471832983368848e-02 +1.616124092435900633e-02 +1.732943630776628136e-02 +1.623603154062895279e-02 +1.591139894430858986e-02 +1.800282118758941841e-02 +1.427327439601875150e-02 +1.408566792305422297e-02 +2.068681520542113225e-02 +2.244182665819811559e-02 +2.503532433376525965e-02 +2.593273429932693372e-02 +2.077254636508810842e-02 +2.104592821805566227e-02 +2.022741370846423137e-02 +2.273147960422196995e-02 +1.713856362243489778e-02 +2.131707246231656291e-02 +2.142076409668075554e-02 +2.106060801220243500e-02 +2.106417197436198999e-02 +1.705112785781429455e-02 +1.532317722387000966e-02 +2.039227698851745971e-02 +2.520659348475314671e-02 +1.913494670223685537e-02 +2.162006407736902572e-02 +2.142434670357279067e-02 +1.666780736122279102e-02 +2.034743811298880822e-02 +1.535548641682486087e-02 +1.890882409275085330e-02 +1.707840925629109524e-02 +1.689179761760981058e-02 +1.625243318093048256e-02 +1.620857637287674916e-02 +1.645092711557727752e-02 +1.726608768074250691e-02 +1.660050430096490087e-02 +1.548794694891195839e-02 +2.870611986956998762e-02 +2.945458783169282654e-02 +2.899514619838083201e-02 +2.661236024326397651e-02 +3.000279383192859181e-02 +2.945963217241364188e-02 +2.822770944415662134e-02 +2.641137280083557734e-02 +3.396396406969628967e-02 +3.687468567388973301e-02 +3.369747874328042625e-02 +3.087412716493228432e-02 +2.616534630859552665e-02 +3.609775970114016685e-02 +2.944279377745562254e-02 +2.978318430829905022e-02 +2.574997524344562264e-02 +1.718285540831799782e-02 +2.178924725667299653e-02 +2.210567529823684665e-02 +1.672196922779749587e-02 +1.555702654288827064e-02 +1.820258090511459731e-02 +1.985293577423892161e-02 +1.965293497206659426e-02 +2.110069272057699466e-02 +2.035995561034906545e-02 +1.408185600528237547e-02 +1.664650730245143331e-02 +1.578114609026634238e-02 +1.667394966117062907e-02 +1.913048679183957171e-02 +2.544338591734448748e-02 +3.346471255409810830e-02 +3.343706716330910322e-02 +2.989762218418577297e-02 +2.460628854561830794e-02 +2.641006296210516788e-02 +2.477084064900300642e-02 +2.588506997943394727e-02 +2.338868647343537999e-02 +2.626587235998870465e-02 +3.156085432942037222e-02 +1.171126968019255282e-01 +2.402704529490212737e-02 +1.666990286876634045e-02 +2.894243531960747837e-02 +1.998438641717136324e-02 +2.777944951319335012e-02 +2.402693800423626180e-02 +3.333226182127600745e-02 +2.501410924714922088e-02 +1.915412700192369727e-02 +1.670355116940346746e-02 +2.855173007172775307e-02 +2.234250923119979193e-02 +2.173065724112052677e-02 +2.076036134141011347e-02 +2.262979460028160736e-02 +2.327813380938380944e-02 +2.198490905867101289e-02 +2.200227418726311002e-02 +2.095901423846064510e-02 +1.842230088712391856e-02 +3.405018691350287369e-02 +3.327348811135773243e-02 +3.380149877991046881e-02 +3.288806274311180389e-02 +3.379411356633452851e-02 +3.439701296928251995e-02 +3.181861044551604484e-02 +3.012124974815062750e-02 +2.008662358649957597e-02 +1.878858808771618691e-02 +2.814384812841709549e-02 +1.779090132615699163e-02 +2.498236841731857338e-02 +2.196810696362214155e-02 +1.362424813904945277e-02 +2.186339995731857375e-02 +1.840118150366428298e-02 +1.566788909835893293e-02 +1.543085385235618348e-02 +1.616013804313997521e-02 +1.867212907277575676e-02 +2.493122275399653129e-02 +1.019894150957608749e-02 +1.841774311675976536e-02 +1.465829549110628759e-02 +1.527281075184000296e-02 +1.954381180327853865e-02 +2.273356830662013120e-02 +2.118592073232443071e-02 +1.880183963419137030e-02 +2.378655853426279651e-02 +2.126848810924386601e-02 +1.988176927090141294e-02 +2.275291508909026150e-02 +1.864752823085735756e-02 +2.372510964400102035e-02 +2.137962335699018779e-02 +2.032723715463959427e-02 +1.742686820164202066e-02 +2.407129576913031102e-02 +2.137689622731629366e-02 +1.984072180545004560e-02 +1.636113665372297687e-02 +2.097922949022315922e-02 +1.982514610857237514e-02 +1.985282555155657574e-02 +1.352882424445514022e-02 +1.517018940973106390e-02 +4.986117562357225003e-02 +1.652583115114606566e-02 +3.588334441257649138e-02 +2.142703043640323615e-02 +3.614748318853107378e-02 +2.005102421717045724e-01 +4.216499255039212607e-02 +2.370986665528723031e-02 +2.794529381084682670e-02 +2.767464642676755598e-02 +2.335316665571461961e-02 +2.854242559137647878e-02 +2.803681480509577495e-02 +2.775352328972690782e-02 +2.674819475604402197e-02 +2.351668746023789433e-02 diff --git a/analysis/tb25_analysis/HidraTB25Ana.C b/analysis/tb25_analysis/HidraTB25Ana.C new file mode 100644 index 0000000..901e40b --- /dev/null +++ b/analysis/tb25_analysis/HidraTB25Ana.C @@ -0,0 +1,1623 @@ +/*********************************************/ +/* HidraTB25Ana.C */ +/*Analysis file for TB25 HidraSim simulation */ +/* Contact: A. Pareti (andrea.pareti@cern.ch)*/ +/*********************************************/ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "../HidraGeo.h" +#include +#include +#include + +using json = nlohmann::json; + +// ----------------------------------------------------------------------------- +// Global constants and parameters +// ----------------------------------------------------------------------------- + +const unsigned int grouping = 8; +const unsigned int EventDisplayEvery = 100; +const unsigned int FersMultiplicity = 2; + +const double chi = 0.38; +const double sciPheGeV = 119.001; // tb24 attenuation +const double cerPheGeV = 29.4; // tb24 attenuation +//const double sciPheGeV = 178.501; // test with 10m attenuation -> remove att. effect +//const double cerPheGeV = 43; // test with 10m attenuation -> remove att. effect + +const double NofSipmCells_sci = 7772; +const double NofSipmCells_cer = 3443; +const double sci_pde = 0.22; // increased wrt to datasheet due to higher overvoltage in TB25 +const double cer_pde = 0.38; // increased wrt to datasheet due to higher overvoltage in TB25 + +// Containment correction factors +const double elcont = 1.005; +const double picont = 1.028; + + +enum class NoiseCorrelationMode { + UncorrelatedByChannel, + CorrelatedWithinFers +}; + +enum class NoiseDistributionMode { + Gaussian, + LogNormal +}; + +///////////////////////////////////////////////////////// +// Paremeters for changing simulation behaviour +// Include/exclude: - saturation, noise, pedestal subtraction +// Noise parameters: - noise correlation mode, noise distribution +// Can also scale pedestal subtraction to try matching observed data linearity +const bool ApplySaturation = true; +const bool ApplyTbNoise = false; +const bool AddPedestalToEnergyContribution = false; +const bool ApplyPedestalSubtraction = false; +const double ScalePedestalSubtractionFactorS = 1.0; +const double ScalePedestalSubtractionFactorC = 1.0; +//const double ScalePedestalSubtractionFactorS = 1.1; // 1.09 for correlated +//const double ScalePedestalSubtractionFactorC = 1.055; // 1.055 for correlated +const double FersThresholdScaleFactor = 1.; +const NoiseCorrelationMode NoiseMode = NoiseCorrelationMode::UncorrelatedByChannel; // use uncorrelated noise +//const NoiseCorrelationMode NoiseMode = NoiseCorrelationMode::CorrelatedWithinFers; // use correlated noise within FERS + +//const NoiseDistributionMode NoiseDistribution = NoiseDistributionMode::LogNormal; +const NoiseDistributionMode NoiseDistribution = NoiseDistributionMode::Gaussian; +const bool printSmearingLog = false; +const bool printFersLog = false; + +/////////////////////////////////////////////////////////// + +const unsigned int ChannelsPerFers = 64; +const unsigned int TotalFers = 16; +const unsigned int TotalCalibChannels = ChannelsPerFers * TotalFers; + + + + + + +using NoiseGeneratorBank = std::vector; + +// ----------------------------------------------------------------------------- +// Utility helpers +// ----------------------------------------------------------------------------- + +// ----------------------------------------------------------------------------- +// 3) Add deterministic RNG builders +// Note: ROOT treats seed=0 specially, so channel seed uses idx+1. +// ----------------------------------------------------------------------------- + +unsigned int MakeChannelNoiseSeed(unsigned int calibIndex){return calibIndex + 1U;} +unsigned int MakeFersNoiseSeed(unsigned int fersId){return fersId;} + +NoiseGeneratorBank BuildChannelNoiseGenerators() +{ + NoiseGeneratorBank generators(TotalCalibChannels); + for (unsigned int calibIndex = 0; calibIndex < TotalCalibChannels; ++calibIndex) { + generators[calibIndex].SetSeed(MakeChannelNoiseSeed(calibIndex)); + } + return generators; +} + +NoiseGeneratorBank BuildFersNoiseGenerators() +{ + NoiseGeneratorBank generators(TotalFers + 1U); + for (unsigned int fersId = 1; fersId <= TotalFers; ++fersId) { + generators[fersId].SetSeed(MakeFersNoiseSeed(fersId)); + } + return generators; +} + + +double apply_sipm_saturation_from_pe(double pe, + double phePerGeV, + double nCells) +{ + if (pe <= 0.0 || phePerGeV <= 0.0 || nCells <= 0.0) { + return 0.0; + } + + const double firedCells = -nCells * std::expm1(-pe / nCells); + return firedCells / phePerGeV; +} + + +// FERS activation values (pedestal subtracted) -> to be used +std::map fers_to_thr_map = { + { 1, 0.17 }, + { 2, 0.20 }, + { 3, 0.30 }, + { 4, 0.28 }, + { 5, 0.105}, + { 6, 0.095}, + { 7, 0.08 }, + { 8, 0.08 }, + { 9, 0.05 }, + {10, 0.07 }, + {11, 0.05 }, + {12, 0.025}, + {13, 0.06 }, + {14, 0.05 }, + {15, 0.12 }, + {16, 0.06 } +}; + +/* +// FERS activation values ( non-pedestal subtracted) +// -> preliminary, measured threhsold values after pedestal subtraction with muons +std::map fers_to_thr_map = { +{1, 0.2225}, +{2, 0.24}, +{3, 0.3508}, +{4, 0.3508}, +{5, 0.1641}, +{6, 0.1525}, +{7, 0.1405}, +{8, 0.0942}, +{9, 0.0592}, +{10, 0.0592}, +{11, 0.0358}, +{12, 0.0825}, +{13, 0.0825}, +{14, 0.0825}, +{15, 0.1413}, +{16, 0.1413} +};*/ + +struct PedestalHgEntry { + double medianAdc = 0.0; + double rmsAdc = 0.0; +}; + +using PedestalHgVector = std::vector; +using AdcToGeVVector = std::vector; + +struct SipmMapEntry { + int calibIndex = -1; + int boardID = -1; + std::string type; + int row = -1; + int column = -1; + double x = 0.0; + double y = 0.0; + int fersId = -1; + int ch = -1; + std::string module_name; + double x_local = 0.0; + double y_local = 0.0; + + double pedestalMedianHgAdc = 0.0; + double pedestalRmsHgAdc = 0.0; + double pedestalMedianHgGeV = 0.0; + double pedestalRmsHgGeV = 0.0; +}; + +struct EventDisplayHit { + uint64_t key = 0; + double x = 0.0; + double y = 0.0; + double rawPhe = 0.0; + double signal = 0.0; +}; + +struct GaussianFitSummary { + double mean = std::numeric_limits::quiet_NaN(); + double meanErr = std::numeric_limits::quiet_NaN(); + double rms = std::numeric_limits::quiet_NaN(); + double rmsErr = std::numeric_limits::quiet_NaN(); +}; + +struct HistMomentSummary { + double mean = std::numeric_limits::quiet_NaN(); + double meanErr = std::numeric_limits::quiet_NaN(); + double rms = std::numeric_limits::quiet_NaN(); + double rmsErr = std::numeric_limits::quiet_NaN(); +}; + +using SipmLookup = std::unordered_map; +using FersKey = uint64_t; + +// Constrct a unique key for each FERS+channel combination for easier lookup and accumulation +struct ChannelKey { + FersKey fersKey = 0; + int ch = -1; + + bool operator==(const ChannelKey& other) const + { + return fersKey == other.fersKey && ch == other.ch; + } +}; + +struct ChannelKeyHash { + std::size_t operator()(const ChannelKey& key) const noexcept + { + const std::size_t h1 = std::hash{}(key.fersKey); + const std::size_t h2 = std::hash{}(key.ch); + return h1 ^ (h2 + 0x9e3779b97f4a7c15ULL + (h1 << 6U) + (h1 >> 2U)); + } +}; + +struct ChannelAccumulatedSignal { + FersKey fersKey = 0; + int fersId = -1; + int ch = -1; + int calibIndex = -1; + bool isSci = false; + double threshold = 0.0; + double thresholdSignal = 0.0; + double outputSignal = 0.0; + double pedestalMean = 0.0; + double noiseSigma = 0.0; + double smearedThresholdSignal = 0.0; + double smearedOutputSignal = 0.0; +}; + + +using ChannelSignalMap = + std::unordered_map; + +std::string make_sipm_key(const std::string& tower, + const std::string& type, + int row, + int column) +{ + std::ostringstream os; + os << tower << '|' << type << '|' << row << '|' << column; + return os.str(); +} + +FersKey make_fers_key(int boardID, int fersId) +{ + return (static_cast(static_cast(boardID)) << 32) | + static_cast(fersId); +} + +int board_id_from_fers_key(FersKey key) +{ + return static_cast(static_cast(key >> 32)); +} + +int grouped_column(int colID_original, unsigned int groupingValue) +{ + return static_cast(colID_original / groupingValue); +} + +double GetFersThreshold(int fersId) +{ + if (fersId < 0) { + throw std::runtime_error("Invalid FERS id: " + std::to_string(fersId)); + } + + const auto it = fers_to_thr_map.find(static_cast(fersId)); + if (it == fers_to_thr_map.end()) { + throw std::runtime_error("Missing threshold for FERS id: " + std::to_string(fersId)); + } + + //return it->second; + return it->second * FersThresholdScaleFactor; // try scaling to consider pedestal non-subtracted values +} + +// ----------------------------------------------------------------------------- +// Calibration loaders +// ----------------------------------------------------------------------------- + +PedestalHgVector LoadHgPedestalsFromJson(const std::string& jsonPath) +{ + std::ifstream in(jsonPath); + if (!in) { + throw std::runtime_error("Cannot open pedestal JSON: " + jsonPath); + } + + json j; + in >> j; + + if (!j.is_object()) { + throw std::runtime_error("Pedestal JSON must be an object: " + jsonPath); + } + + std::size_t maxIndex = 0; + for (auto it = j.begin(); it != j.end(); ++it) { + const int idx = std::stoi(it.key()); + if (idx < 0) { + throw std::runtime_error("Negative pedestal index in: " + jsonPath); + } + maxIndex = std::max(maxIndex, static_cast(idx)); + } + + PedestalHgVector values(maxIndex + 1); + std::vector found(maxIndex + 1, false); + + for (auto it = j.begin(); it != j.end(); ++it) { + const int idx = std::stoi(it.key()); + const json& node = it.value(); + + values[static_cast(idx)].medianAdc = node.value("median_HG", 0.0); + values[static_cast(idx)].rmsAdc = node.value("iqr_eff_HG", 0.0); + found[static_cast(idx)] = true; + } + + for (std::size_t i = 0; i < found.size(); ++i) { + if (!found[i]) { + std::ostringstream os; + os << "Missing HG pedestal entry for channel " << i + << " in " << jsonPath; + throw std::runtime_error(os.str()); + } + } + + return values; +} + +AdcToGeVVector LoadAdcToGeVFromJson(const std::string& jsonPath) +{ + std::ifstream in(jsonPath); + if (!in) { + throw std::runtime_error("Cannot open ADC->GeV JSON: " + jsonPath); + } + + json j; + in >> j; + + if (!j.is_array()) { + throw std::runtime_error("ADC->GeV JSON must be an array: " + jsonPath); + } + + AdcToGeVVector values; + values.reserve(j.size()); + + for (const auto& item : j) { + values.push_back(item.get()); + } + + if (values.empty()) { + throw std::runtime_error("ADC->GeV JSON is empty: " + jsonPath); + } + + return values; +} + +// not used in current version +std::size_t BuildCalibrationIndex(const SipmMapEntry& entry) +{ + if (entry.fersId <= 0 || entry.fersId > static_cast(TotalFers)) { + std::ostringstream os; + os << "Invalid FERS id " << entry.fersId; + throw std::runtime_error(os.str()); + } + + if (entry.ch < 0 || entry.ch >= static_cast(ChannelsPerFers)) { + std::ostringstream os; + os << "Invalid channel " << entry.ch + << " for FERS " << entry.fersId; + throw std::runtime_error(os.str()); + } + + const std::size_t idx = + (static_cast(entry.fersId) - 1U) * ChannelsPerFers + + static_cast(entry.ch); + + if (idx >= TotalCalibChannels) { + std::ostringstream os; + os << "Calibration index out of range: " << idx + << " from FERS " << entry.fersId + << ", ch " << entry.ch; + throw std::runtime_error(os.str()); + } + + return idx; +} + +void AttachHgPedestalNoiseToSipmMap(SipmLookup& lookup, + const PedestalHgVector& pedestalsHg, + const AdcToGeVVector& adcToGeV) +{ + if (pedestalsHg.size() != TotalCalibChannels) { + std::ostringstream os; + os << "Expected " << TotalCalibChannels + << " HG pedestal entries, got " << pedestalsHg.size(); + throw std::runtime_error(os.str()); + } + + if (adcToGeV.size() != TotalCalibChannels) { + std::ostringstream os; + os << "Expected " << TotalCalibChannels + << " ADC->GeV entries, got " << adcToGeV.size(); + throw std::runtime_error(os.str()); + } + + for (auto& [key, entry] : lookup) { + if (entry.calibIndex < 0 || entry.calibIndex >= static_cast(TotalCalibChannels)) { + std::ostringstream os; + os << "Invalid calibIndex " << entry.calibIndex + << " for map entry " << key; + throw std::runtime_error(os.str()); + } + + const std::size_t idx = static_cast(entry.calibIndex); + + entry.pedestalMedianHgAdc = pedestalsHg[idx].medianAdc; + entry.pedestalRmsHgAdc = pedestalsHg[idx].rmsAdc; + entry.pedestalMedianHgGeV = pedestalsHg[idx].medianAdc * adcToGeV[idx]; + entry.pedestalRmsHgGeV = pedestalsHg[idx].rmsAdc * adcToGeV[idx]; + } +} + + + +// ----------------------------------------------------------------------------- +// Signal accumulation in each channel and smearing +// Apply saturation to a single signal path for both threshold and output consistency +// ----------------------------------------------------------------------------- +void AccumulateChannelSignal(ChannelSignalMap& channelSignals, + const SipmMapEntry& info, + double signalContribution) +{ + if (info.boardID < 0 || info.fersId < 0 || info.ch < 0 || info.calibIndex < 0) { + return; + } + + const FersKey fersKey = make_fers_key(info.boardID, info.fersId); + const ChannelKey channelKey{fersKey, info.ch}; + + auto [it, inserted] = channelSignals.emplace( + channelKey, + ChannelAccumulatedSignal{ + fersKey, + info.fersId, + info.ch, + info.calibIndex, + info.type == "S", + GetFersThreshold(info.fersId), + 0.0, + 0.0, + info.pedestalMedianHgGeV, + info.pedestalRmsHgGeV, + 0.0, + 0.0 + } + ); + + if (!inserted && it->second.calibIndex != info.calibIndex) { + std::ostringstream os; + os << "Inconsistent calibIndex for accumulated channel: existing=" + << it->second.calibIndex << ", new=" << info.calibIndex + << " for FERS " << info.fersId << ", ch " << info.ch; + throw std::runtime_error(os.str()); + } + + // Use single signal for both threshold and output to ensure consistency + it->second.thresholdSignal += signalContribution; + it->second.outputSignal += signalContribution; +} + +double DrawUnitGaussianNoiseForChannel(const ChannelAccumulatedSignal& channel, + NoiseGeneratorBank& channelNoiseGenerators) +{ + if (channel.calibIndex < 0 || + channel.calibIndex >= static_cast(channelNoiseGenerators.size())) { + std::ostringstream os; + os << "Invalid calibIndex " << channel.calibIndex + << " while drawing channel noise"; + throw std::runtime_error(os.str()); + } + + return channelNoiseGenerators[static_cast(channel.calibIndex)].Gaus(0.0, 1.0); +} + +double DrawUnitGaussianNoiseForFers( + const ChannelAccumulatedSignal& channel, + NoiseGeneratorBank& fersNoiseGenerators, + std::unordered_map& fersEventNoiseCache) +{ + const auto it = fersEventNoiseCache.find(channel.fersKey); + if (it != fersEventNoiseCache.end()) { + return it->second; + } + + if (channel.fersId <= 0 || + channel.fersId >= static_cast(fersNoiseGenerators.size())) + { + throw std::runtime_error( + "Invalid FERS id " + std::to_string(channel.fersId) + + " while drawing correlated FERS noise" + ); + } + + + const double sharedUnitNoise = + fersNoiseGenerators[static_cast(channel.fersId)].Gaus(0.0, 1.0); + + fersEventNoiseCache[channel.fersKey] = sharedUnitNoise; + return sharedUnitNoise; +} + + +struct LogNormalParameters { + double mu = 0.0; + double sigma = 0.0; +}; + +LogNormalParameters GetLogNormalParametersFromMeanRms(double mean, double rms) +{ + if (mean <= 0.0 || rms <= 0.0) { + return {0.0, 0.0}; + } + + const double variance = rms * rms; + + LogNormalParameters pars; + pars.sigma = std::sqrt(std::log(1.0 + variance / (mean * mean))); + pars.mu = std::log(mean) - 0.5 * pars.sigma * pars.sigma; + + return pars; +} + +double DrawLogNormalWithMeanRms(TRandom3& rng, double mean, double rms) +{ + if (mean <= 0.0) { + return 0.0; + } + + if (rms <= 0.0) { + return mean; + } + + const LogNormalParameters pars = + GetLogNormalParametersFromMeanRms(mean, rms); + + return std::exp(rng.Gaus(pars.mu, pars.sigma)); +} + + + +// ----------------------------------------------------------------------------- +//Smearing: noise first, activation decision first, no pedestal subtraction here +void SmearAndCountChannelsOverThreshold( + ChannelSignalMap& channelSignals, + std::unordered_map& fersChannelsOverThreshold, + NoiseCorrelationMode noiseMode, + NoiseGeneratorBank& channelNoiseGenerators, + NoiseGeneratorBank& fersNoiseGenerators, + std::unordered_map& fersEventNoiseCache) +{ + + for (auto& [channelKey, channel] : channelSignals) { + + double pedestalContribution = 0.0; + + if (AddPedestalToEnergyContribution) { + + if (!ApplyTbNoise || channel.noiseSigma <= 0.0) { + pedestalContribution = channel.pedestalMean; + } + + else if (NoiseDistribution == NoiseDistributionMode::Gaussian) { + double unitNoise = 0.0; + + if (noiseMode == NoiseCorrelationMode::UncorrelatedByChannel) { + unitNoise = DrawUnitGaussianNoiseForChannel( + channel, + channelNoiseGenerators + ); + } else if (noiseMode == NoiseCorrelationMode::CorrelatedWithinFers) { + unitNoise = DrawUnitGaussianNoiseForFers( + channel, + fersNoiseGenerators, + fersEventNoiseCache + ); + } + + pedestalContribution = + channel.pedestalMean + unitNoise * channel.noiseSigma; + } + + else if (NoiseDistribution == NoiseDistributionMode::LogNormal) { + if (noiseMode == NoiseCorrelationMode::UncorrelatedByChannel) { + pedestalContribution = DrawLogNormalWithMeanRms( + channelNoiseGenerators[static_cast(channel.calibIndex)], + channel.pedestalMean, + channel.noiseSigma + ); + } else if (noiseMode == NoiseCorrelationMode::CorrelatedWithinFers) { + pedestalContribution = DrawLogNormalWithMeanRms( + fersNoiseGenerators[static_cast(channel.fersId)], + channel.pedestalMean, + channel.noiseSigma + ); + } + } + } + + channel.smearedThresholdSignal = + channel.thresholdSignal + pedestalContribution; + + channel.smearedOutputSignal = + channel.outputSignal + pedestalContribution; + + + + + + if (channel.smearedThresholdSignal > channel.threshold) { + ++fersChannelsOverThreshold[channel.fersKey]; + } + + + + + } +} + + + + +bool IsActivatedFers(const std::unordered_map& fersChannelsOverThreshold, + FersKey key) +{ + const auto it = fersChannelsOverThreshold.find(key); + return it != fersChannelsOverThreshold.end() && + it->second >= FersMultiplicity; +} + +double SumActivatedSmearedChannelOutput( + const ChannelSignalMap& channelSignals, + const std::unordered_map& fersChannelsOverThreshold) +{ + double total = 0.0; + + for (const auto& [channelKey, channel] : channelSignals) { + if (!IsActivatedFers(fersChannelsOverThreshold, channel.fersKey)) { + continue; + } + + total += channel.smearedOutputSignal; + } + + return total; +} + +unsigned int CountActivatedFers( + const std::unordered_map& fersChannelsOverThreshold) +{ + unsigned int activated = 0; + + for (const auto& [fersKey, nChannels] : fersChannelsOverThreshold) { + if (nChannels >= FersMultiplicity) { + ++activated; + } + } + + return activated; +} + + +// ----------------------------------------------------------------------------- +// 4) subtract pedestal only for activated FERS +void SubtractPedestalMedianFromActivatedChannels( + ChannelSignalMap& channelSignals, + const std::unordered_map& fersChannelsOverThreshold) +{ + for (auto& [channelKey, channel] : channelSignals) { + if (!IsActivatedFers(fersChannelsOverThreshold, channel.fersKey)) { + continue; + } + + if(ApplyPedestalSubtraction) + { + const double scale = channel.isSci + ? ScalePedestalSubtractionFactorS + : ScalePedestalSubtractionFactorC; + + //channel.smearedOutputSignal -= channel.pedestalMean; + if(channel.outputSignal>0 and printSmearingLog){ + std::cout << "Subtracting pedestal median for FERS " << channel.fersId + << ", ch " << channel.ch + << ", true signal: " << channel.outputSignal + << "\t smeared signal: " << channel.smearedOutputSignal + << "\t pedestal median: " << channel.pedestalMean + << "\t scale: " << scale << std::endl;} + channel.smearedOutputSignal -= channel.pedestalMean * scale; + } + // Allow negative values: downward fluctuations naturally subtract from total energy + } +} + +void FillActivatedChannelOutputArrays( + const ChannelSignalMap& channelSignals, + const std::unordered_map& fersChannelsOverThreshold, + std::vector& activatedChannelOutput, + std::vector& activatedChannelOutputTrue, + std::vector& activatedChannelOutputSmeared) +{ + if (activatedChannelOutput.size() != activatedChannelOutputTrue.size() || + activatedChannelOutput.size() != activatedChannelOutputSmeared.size()) { + throw std::runtime_error("Activated channel output arrays must have the same size"); + } + + for (const auto& [channelKey, channel] : channelSignals) { + if (channel.calibIndex < 0 || + channel.calibIndex >= static_cast(activatedChannelOutput.size())) { + std::ostringstream os; + os << "Invalid calibIndex " << channel.calibIndex + << " while filling activated channel output arrays"; + throw std::runtime_error(os.str()); + } + + if (!IsActivatedFers(fersChannelsOverThreshold, channel.fersKey)) { + continue; + } + + const std::size_t idx = static_cast(channel.calibIndex); + + activatedChannelOutputTrue[idx] = channel.outputSignal; + activatedChannelOutputSmeared[idx] = channel.smearedOutputSignal; + + // keep the legacy branch identical to the smeared output + activatedChannelOutput[idx] = channel.smearedOutputSignal; + } +} + + + + +void PrintPerEventFersChannelCounts( + unsigned int eventIndex, + const ChannelSignalMap& sciChannelSignals, + const ChannelSignalMap& cerChannelSignals, + const std::unordered_map& fersChannelsOverThreshold) +{ + std::map, FersKey> fersSeen; + + auto collectFers = [&](const ChannelSignalMap& channelSignals) { + for (const auto& [channelKey, channel] : channelSignals) { + const int boardID = board_id_from_fers_key(channel.fersKey); + fersSeen[{boardID, channel.fersId}] = channel.fersKey; + } + }; + + collectFers(sciChannelSignals); + collectFers(cerChannelSignals); + + if (printFersLog) { + std::cout << "Event " << eventIndex << ":\n"; + for (const auto& [id, key] : fersSeen) { + const int boardID = id.first; + const int fersId = id.second; + + const auto it = fersChannelsOverThreshold.find(key); + const unsigned int nChannelsOverThreshold = + (it != fersChannelsOverThreshold.end()) ? it->second : 0U; + + std::cout + << " board " << boardID + << ", FERS " << fersId + << " -> channels above threshold = " << nChannelsOverThreshold; + + if (nChannelsOverThreshold >= FersMultiplicity) { + std::cout << " [ON]"; + } else { + std::cout << " [OFF]"; + } + + std::cout << '\n'; + } + } +} + + + + +// ----------------------------------------------------------------------------- +// SiPM map and event display helpers +// ----------------------------------------------------------------------------- +SipmLookup LoadSipmMap(const std::string& jsonPath) +{ + std::ifstream in(jsonPath); + if (!in) { + throw std::runtime_error("Cannot open SiPM map JSON: " + jsonPath); + } + + json j; + in >> j; + + SipmLookup lookup; + for (auto it = j.begin(); it != j.end(); ++it) { + const json& node = it.value(); + + SipmMapEntry entry; + entry.calibIndex = std::stoi(it.key()); + entry.boardID = node.value("boardID", -1); + entry.type = node.value("type", ""); + entry.row = node.value("row", -1); + entry.column = node.value("column", -1); + entry.x = node.value("x", 0.0); + entry.y = node.value("y", 0.0); + entry.fersId = node.value("fersId", -1); + entry.ch = node.value("ch", -1); + entry.module_name = node.value("module_name", ""); + entry.x_local = node.value("x_local", 0.0); + entry.y_local = node.value("y_local", 0.0); + + lookup[make_sipm_key(entry.module_name, entry.type, entry.row, entry.column)] = entry; + } + + return lookup; +} + +const SipmMapEntry* FindSipmInfo(const SipmLookup& lookup, + const std::string& tower, + const std::string& type, + int row, + int groupedCol, + unsigned int /*groupingValue*/) +{ + const int jsonColumn = groupedCol; + const auto key = make_sipm_key(tower, type, row, jsonColumn); + const auto it = lookup.find(key); + return (it == lookup.end()) ? nullptr : &it->second; +} + +std::string make_event_display_name(const std::string& prefix, unsigned int entry) +{ + std::ostringstream os; + os << prefix << "_evt" << std::setw(6) << std::setfill('0') << entry; + return os.str(); +} + +TH2F* CreateEventDisplayHist(const std::string& name, + const std::string& title, + unsigned int groupingValue) +{ + const int xBins = static_cast(NofSiPMTowersX * NofFiberscolumn / groupingValue); + const int yBins = NofSiPMTowersY * NofFibersrow; + + return new TH2F(name.c_str(), title.c_str(), + xBins, + -NofSiPMTowersX * moduleX / 2.0, NofSiPMTowersX * moduleX / 2.0, + yBins, + -NofSiPMTowersY * moduleY / 2.0, NofSiPMTowersY * moduleY / 2.0); +} + +void FillActivatedEventDisplay(TH2F* hist, + const std::vector& hits, + const std::unordered_map& fersChannelsOverThreshold) +{ + if (!hist) { + return; + } + + for (const auto& hit : hits) { + if (IsActivatedFers(fersChannelsOverThreshold, hit.key)) { + hist->Fill(hit.x, hit.y, hit.signal); + } + } +} + +void GetSiPMcoordinate(int TowID, + int rowID, + int colID_original, + double& SiPM_X, + double& SiPM_Y, + std::string fiber, + unsigned int groupingValue) +{ + const double towerPitchY = dtubeY * NofFibersrow; + double TowerOffsetY = -((NoModulesSiPM - 1) * towerPitchY) / 2 + TowID * towerPitchY; + unsigned int channel = static_cast(colID_original / groupingValue); + double colID = (static_cast(groupingValue) - 1) / 2 + channel * groupingValue; + + if (fiber == "S") { + SiPM_X = +moduleX / 2 - tuberadius - (tuberadius * 2) * colID; + SiPM_Y = -moduleY / 2 + tuberadius + (sq3 * tuberadius) * rowID + tuberadius * (2. * sq3m1 - 1.); + } + + if (fiber == "C") { + SiPM_X = +moduleX / 2 - (tuberadius * 2) - (tuberadius * 2) * colID; + SiPM_Y = -moduleY / 2 + tuberadius + (sq3 * tuberadius) * rowID + tuberadius * (2. * sq3m1 - 1.); + } + + SiPM_Y = TowerOffsetY + SiPM_Y; +} + +// ----------------------------------------------------------------------------- +// Analysis/output helpers +// ----------------------------------------------------------------------------- + +std::string simTower_to_tbTower(std::string simTower) +{ + return std::to_string(306 + std::stoi(simTower)); +} + +GaussianFitSummary FitGaussianSummary(TH1* hist) +{ + GaussianFitSummary out; + if (!hist || hist->GetEntries() < 10) { + return out; + } + + hist->Fit("gaus", "Q0"); + TF1* fit = hist->GetFunction("gaus"); + if (!fit) { + return out; + } + + out.mean = fit->GetParameter(1); + out.meanErr = fit->GetParError(1); + out.rms = fit->GetParameter(2); + out.rmsErr = fit->GetParError(2); + return out; +} + +HistMomentSummary GetHistMomentSummary(const TH1* hist) +{ + HistMomentSummary out; + if (!hist || hist->GetEntries() < 1) { + return out; + } + + out.mean = hist->GetMean(); + out.meanErr = hist->GetMeanError(); + out.rms = hist->GetRMS(); + out.rmsErr = hist->GetRMSError(); + return out; +} + +std::string CsvEscape(const std::string& value) +{ + std::string escaped = "\""; + for (char c : value) { + if (c == '"') { + escaped += "\"\""; + } else { + escaped += c; + } + } + escaped += "\""; + return escaped; +} + +void AppendSummaryCsvLocked(const std::string& csvPath, + const std::string& inputFile, + double truthEnergy, + int nentries, + const GaussianFitSummary& sFit, + const GaussianFitSummary& cFit, + const GaussianFitSummary& combFit, + const GaussianFitSummary& combChiFit, + const HistMomentSummary& sciX, + const HistMomentSummary& sciY, + const HistMomentSummary& cerX, + const HistMomentSummary& cerY, + const GaussianFitSummary& sFitFersOn, + const GaussianFitSummary& cFitFersOn) +{ + const int fd = open(csvPath.c_str(), O_CREAT | O_WRONLY | O_APPEND, 0644); + if (fd == -1) { + throw std::runtime_error("Cannot open summary CSV: " + csvPath); + } + + if (flock(fd, LOCK_EX) != 0) { + close(fd); + throw std::runtime_error("Cannot lock summary CSV: " + csvPath); + } + + struct stat st; + if (fstat(fd, &st) != 0) { + flock(fd, LOCK_UN); + close(fd); + throw std::runtime_error("Cannot stat summary CSV: " + csvPath); + } + + std::ostringstream out; + out << std::setprecision(17); + + if (st.st_size == 0) { + out + << "input_file,truth_energy,nentries," + << "s_fit_mean,s_fit_mean_err,s_fit_rms,s_fit_rms_err," + << "c_fit_mean,c_fit_mean_err,c_fit_rms,c_fit_rms_err," + << "comb_fit_mean,comb_fit_mean_err,comb_fit_rms,comb_fit_rms_err," + << "combchi_fit_mean,combchi_fit_mean_err,combchi_fit_rms,combchi_fit_rms_err," + << "res_sci_x_mean,res_sci_x_mean_err,res_sci_x_rms,res_sci_x_rms_err," + << "res_sci_y_mean,res_sci_y_mean_err,res_sci_y_rms,res_sci_y_rms_err," + << "res_cer_x_mean,res_cer_x_mean_err,res_cer_x_rms,res_cer_x_rms_err," + << "res_cer_y_mean,res_cer_y_mean_err,res_cer_y_rms,res_cer_y_rms_err," + << "s_fit_fers_on_mean,s_fit_fers_on_mean_err,s_fit_fers_on_rms,s_fit_fers_on_rms_err," + << "c_fit_fers_on_mean,c_fit_fers_on_mean_err,c_fit_fers_on_rms,c_fit_fers_on_rms_err\n"; + } + + out + << CsvEscape(inputFile) << ',' + << truthEnergy << ',' + << nentries << ',' + << sFit.mean << ',' << sFit.meanErr << ',' << sFit.rms << ',' << sFit.rmsErr << ',' + << cFit.mean << ',' << cFit.meanErr << ',' << cFit.rms << ',' << cFit.rmsErr << ',' + << combFit.mean << ',' << combFit.meanErr << ',' << combFit.rms << ',' << combFit.rmsErr << ',' + << combChiFit.mean << ',' << combChiFit.meanErr << ',' << combChiFit.rms << ',' << combChiFit.rmsErr << ',' + << sciX.mean << ',' << sciX.meanErr << ',' << sciX.rms << ',' << sciX.rmsErr << ',' + << sciY.mean << ',' << sciY.meanErr << ',' << sciY.rms << ',' << sciY.rmsErr << ',' + << cerX.mean << ',' << cerX.meanErr << ',' << cerX.rms << ',' << cerX.rmsErr << ',' + << cerY.mean << ',' << cerY.meanErr << ',' << cerY.rms << ',' << cerY.rmsErr << ',' + << sFitFersOn.mean << ',' << sFitFersOn.meanErr << ',' << sFitFersOn.rms << ',' << sFitFersOn.rmsErr << ',' + << cFitFersOn.mean << ',' << cFitFersOn.meanErr << ',' << cFitFersOn.rms << ',' << cFitFersOn.rmsErr << '\n'; + + const std::string text = out.str(); + const ssize_t written = write(fd, text.c_str(), text.size()); + if (written != static_cast(text.size())) { + flock(fd, LOCK_UN); + close(fd); + throw std::runtime_error("Cannot write summary CSV: " + csvPath); + } + + flock(fd, LOCK_UN); + close(fd); +} + + + + + + + + + + +// ----------------------------------------------------------------------------- +// Main analysis +// ----------------------------------------------------------------------------- + +void HidraTB25Ana(double energy, const std::string& input) +{ + // read SiPM map and calibration files + const std::string sipmMapPath = + "../../../TBDataPreparation/2025_SPS/MapAndCalibration/sipm_map.json"; + + const std::string pedestalJsonPath = + "../../../TBDataPreparation/2025_SPS/MapAndCalibration/SiPM_pedestals_v1.json"; + + const std::string adcToGeVJsonPath = + "../../../TBDataPreparation/2025_SPS/MapAndCalibration/SiPM_ADCtoGeV_v1.json"; + + // Prepare lookup map for SiPMs, with attached pedestal and noise info + // key: "tower|type|row|column", value: map entry with FERS/channel info and pedestal/noise + SipmLookup sipmLookup; + try { + sipmLookup = LoadSipmMap(sipmMapPath); + + const PedestalHgVector pedestalsHg = LoadHgPedestalsFromJson(pedestalJsonPath); + const AdcToGeVVector adcToGeV = LoadAdcToGeVFromJson(adcToGeVJsonPath); + + AttachHgPedestalNoiseToSipmMap(sipmLookup, pedestalsHg, adcToGeV); + + std::cout << "Loaded " << sipmLookup.size() + << " SiPM map entries from " << sipmMapPath << std::endl; + std::cout << "Loaded " << pedestalsHg.size() + << " HG pedestal entries from " << pedestalJsonPath << std::endl; + std::cout << "Loaded " << adcToGeV.size() + << " ADC->GeV entries from " << adcToGeVJsonPath << std::endl; + } catch (const std::exception& e) { + std::cerr << e.what() << std::endl; + return; + } + + // read input file and prepare output file + //const std::string infile = "../build/" + input; + const std::string infile = "" + input; + std::cout << "Using file: " << infile << std::endl; + + TFile* simfile = TFile::Open(infile.c_str(), "READ"); + if (!simfile || simfile->IsZombie()) { + std::cerr << "Cannot open input file " << infile << std::endl; + return; + } + + TTree* simtree = static_cast(simfile->Get("DREMTubesout")); + if (!simtree) { + std::cerr << "Cannot find TTree DREMTubesout in " << infile << std::endl; + simfile->Close(); + delete simfile; + return; + } + + std::ostringstream os; + os << energy; + const std::string enstr = os.str(); + const std::string outfile = "hidra" + enstr + ".root"; + TFile f(outfile.c_str(), "RECREATE"); + + TDirectory* eventDisplayDir = f.mkdir("EventDisplays"); + f.cd(); + + TTree* activatedChannelsTree = + new TTree("ActivatedChannels", "Per-event activated channel output"); + + unsigned int outEvent = 0; + std::vector activatedChannelOutput(TotalCalibChannels, 0.0); + std::vector activatedChannelOutputTrue(TotalCalibChannels, 0.0); + std::vector activatedChannelOutputSmeared(TotalCalibChannels, 0.0); + + activatedChannelsTree->Branch("event", &outEvent); + activatedChannelsTree->Branch("activatedChannelOutput", &activatedChannelOutput); + activatedChannelsTree->Branch("activatedChannelOutputTrue", &activatedChannelOutputTrue); + activatedChannelsTree->Branch("activatedChannelOutputSmeared", &activatedChannelOutputSmeared); + + // --------------------------------------- + // create the RNG banks once + // Prepare random generators for noise smearing: + // one per calibration channel for uncorrelated noise, one per FERS for correlated noise + // ---------------------------------------- + NoiseGeneratorBank channelNoiseGenerators = BuildChannelNoiseGenerators(); + NoiseGeneratorBank fersNoiseGenerators = BuildFersNoiseGenerators(); + + + + // Load module mapping to fill the event display maps + int modcol[NofmodulesX * NofmodulesY]; + int modrow[NofmodulesX * NofmodulesY]; + for (int i = 0; i < NofmodulesX * NofmodulesY; i++) { + int row = i / NofmodulesX; + int col = i % NofmodulesX; + int imod = modflag[i]; + if (imod >= 0) { + modcol[imod] = col; + modrow[imod] = row; + } + } + + // ---------------------------------------- + // Prepare output histograms + // ---------------------------------------- + const double bmin = energy - 0.4 * std::sqrt(energy) * 10.; + const double bmax = energy + 0.4 * std::sqrt(energy) * 10.; + + auto sciene = new TH1F("sciene", "sciene", 100, bmin, bmax); + auto cerene = new TH1F("cerene", "cerene", 100, bmin, bmax); + auto totene = new TH1F("totene", "totene", 100, bmin, bmax); + auto totenec = new TH1F("totenec", "totenec", 100, bmin, bmax); + auto totdep = new TH1F("totdep", "totdep", 100, 0., bmax); + auto leakene = new TH1F("leakene", "leakene", 100, 0., 0.1); + auto chidist = new TH1F("chidist", "chidist", 100, 0., 1.); + auto mapcalo = new TH2F("mapcalo", "mapcalo", + NofmodulesX, 0., NofmodulesX, + NofmodulesY, 0., NofmodulesY); + auto SipmMapS = new TH2F("SipmMapS", "SipmS; Col; Row", + NofSiPMTowersX * NofFiberscolumn, 0, NofSiPMTowersX * NofFiberscolumn, + NofSiPMTowersY * NofFibersrow, 0, NofSiPMTowersY * NofFibersrow); + auto SipmMapC = new TH2F("SipmMapC", "SipmC; Col; Row", + NofSiPMTowersX * NofFiberscolumn, 0, NofSiPMTowersX * NofFiberscolumn, + NofSiPMTowersY * NofFibersrow, 0, NofSiPMTowersY * NofFibersrow); + + auto SciSiPMCoordinates = new TH2F("SciSiPMCoordinates", "Sci SiPM Coordinates; X [mm]; Y[mm]", + NofSiPMTowersX * NofFiberscolumn / grouping, + -NofSiPMTowersX * moduleX / 2, NofSiPMTowersX * moduleX / 2, + NofSiPMTowersY * NofFibersrow, + -NofSiPMTowersY * moduleY / 2, NofSiPMTowersY * moduleY / 2); + auto CerSiPMCoordinates = new TH2F("CerSiPMCoordinates", "Cer SiPM Coordinates; X [mm]; Y[mm]", + NofSiPMTowersX * NofFiberscolumn / grouping, + -NofSiPMTowersX * moduleX / 2, NofSiPMTowersX * moduleX / 2, + NofSiPMTowersY * NofFibersrow, + -NofSiPMTowersY * moduleY / 2, NofSiPMTowersY * moduleY / 2); + + auto ResidualHistSciX = new TH1F("ResidualHistSciX", "Residual Sci X; Residual [mm]; Entries", + 100, -10, 10); + auto ResidualHistSciY = new TH1F("ResidualHistSciY", "Residual Sci Y; Residual [mm]; Entries", + 100, -10, 10); + auto ResidualHistCerX = new TH1F("ResidualHistCerX", "Residual Cer X; Residual [mm]; Entries", + 100, -10, 10); + auto ResidualHistCerY = new TH1F("ResidualHistCerY", "Residual Cer Y; Residual [mm]; Entries", + 100, -10, 10); + + auto SignalSfersID = new TH1F("SignalSfersID", "Signal S fibers FERS ID; FERS; Signal", 16, 0.5, 16.5); + auto SignalCfersID = new TH1F("SignalCfersID", "Signal C fibers FERS ID; FERS; Signal", 16, 0.5, 16.5); + + auto scieneFersOn = new TH1F("scieneFersOn", "S energy for Activated FERS", 100, 0., bmax); + auto cereneFersOn = new TH1F("cereneFersOn", "C energy for Activated FERS", 100, 0., bmax); + + auto h2_SvsC = new TH2F("SvsC", "S vs C energy (all fibres); S/E; C/E", 100, 0., 1.5, 100, 0., 1.5); + auto h2_SvsCFersOn = new TH2F("SvsCFersOn", "S vs C energy for Activated FERS; S/E; C/E", 100, 0., 1.5, 100, 0., 1.5); + + // Count number of activated channels per FERS, to check the multiplicity distribution and the effect of different thresholds + auto ActivatedChannelsPerFers = new TH2F("ActivatedChannelsPerFers", "Activated channels per FERS; FERS ID; Channels above threshold",16, 0.5, 16.5, 65, -0.5, 64.5); + + + const int nentries = simtree->GetEntries(); + std::cout << "Entries " << nentries << std::endl; + + int pdg; simtree->SetBranchAddress("PrimaryPDGID", &pdg); + double venergy; simtree->SetBranchAddress("PrimaryParticleEnergy", &venergy); + double lenergy; simtree->SetBranchAddress("EscapedEnergyl", &lenergy); + double denergy; simtree->SetBranchAddress("EscapedEnergyd", &denergy); + double edep; simtree->SetBranchAddress("EnergyTot", &edep); + double Stot; simtree->SetBranchAddress("NofPMTScinDet", &Stot); + double Ctot; simtree->SetBranchAddress("NofPMTCherDet", &Ctot); + double PSdep; simtree->SetBranchAddress("PSEnergy", &PSdep); + double beamX; simtree->SetBranchAddress("PrimaryX", &beamX); + double beamY; simtree->SetBranchAddress("PrimaryY", &beamY); + std::vector* TowerE = nullptr; simtree->SetBranchAddress("VecTowerE", &TowerE); + std::vector* SPMT = nullptr; simtree->SetBranchAddress("VecSPMT", &SPMT); + std::vector* CPMT = nullptr; simtree->SetBranchAddress("VecCPMT", &CPMT); + std::vector* SSiPM = nullptr; simtree->SetBranchAddress("VectorSignals", &SSiPM); + std::vector* CSiPM = nullptr; simtree->SetBranchAddress("VectorSignalsCher", &CSiPM); + + for (unsigned int i = 0; i < static_cast(nentries); i++) { + simtree->GetEntry(i); + + // store one event display every N events + const bool writeEventDisplay = + (EventDisplayEvery > 0) && (((i + 1) % EventDisplayEvery) == 0); + + double ecalo = energy - lenergy / 1000.0; + double totsci = 0.; + double totcer = 0.; + double tottow = 0.; + + double barX_sci = 0.; + double barY_sci = 0.; + double barX_cer = 0.; + double barY_cer = 0.; + + outEvent = i; + // prepare the per-event data structures for the smearing and activation steps + std::fill(activatedChannelOutput.begin(), activatedChannelOutput.end(), 0.0); + std::fill(activatedChannelOutputTrue.begin(), activatedChannelOutputTrue.end(), 0.0); + std::fill(activatedChannelOutputSmeared.begin(), activatedChannelOutputSmeared.end(), 0.0); + std::unordered_map fersChannelsOverThreshold; + + std::vector eventDisplaySciHits; + std::vector eventDisplayCerHits; + + ChannelSignalMap sciChannelSignals; + ChannelSignalMap cerChannelSignals; + + sciChannelSignals.reserve(SSiPM->size()); + cerChannelSignals.reserve(CSiPM->size()); + + if (writeEventDisplay) { + eventDisplaySciHits.reserve(SSiPM->size()); + eventDisplayCerHits.reserve(CSiPM->size()); + } + + // Fill energy in PMT towers + for (unsigned int j = 0; j < SPMT->size(); j++) { + totsci += SPMT->at(j) / sciPheGeV; + totcer += CPMT->at(j) / cerPheGeV; + tottow += TowerE->at(j); + mapcalo->Fill(modcol[j], modrow[j], TowerE->at(j) / 1000. / nentries); + } + + double sciPosWeight = 0.0; + double cerPosWeight = 0.0; + + // Prepare shared random map to have correlated noise between S/C channels + std::unordered_map fersEventNoiseCache; + + + // Start looping on SiPM signals + for (unsigned int n = 0; n < SSiPM->size(); n++) { + const double rawPhe = SSiPM->at(n); // number of photoelectrons + const double content = rawPhe / sciPheGeV; // converted to GeV + + totsci += content; + + const unsigned int towID = static_cast(n / (NofFiberscolumn * NofFibersrow / 2)); + const unsigned int SiPMID = n % (NofFiberscolumn * NofFibersrow / 2); + //const unsigned int colID = NofFiberscolumn - static_cast(SiPMID / (NofFibersrow / 2)); + const unsigned int colID = NofFiberscolumn -1U - static_cast(SiPMID / (NofFibersrow / 2)); + const unsigned int rowID = 2 * static_cast(SiPMID % (NofFibersrow / 2)); + + double SiPM_X = 0.0; + double SiPM_Y = 0.0; + + GetSiPMcoordinate(towID, rowID, colID, SiPM_X, SiPM_Y, "S", grouping); // Currently not used, JSON map is read instead + + const std::string tbTower = simTower_to_tbTower(std::to_string(towID)); // convert from sim tower ID to test beam tower ID + const int groupedCol = grouped_column(static_cast(colID), grouping); // single fibre to channel mapping + + SipmMapS->Fill(colID, towID * NofFibersrow + rowID, content); + + // now associate the SiPM signal to the corresponding FERS and channel + const SipmMapEntry* info = + FindSipmInfo(sipmLookup, tbTower, "S", static_cast(rowID), groupedCol, grouping); + + if (info) { + SciSiPMCoordinates->Fill(info->x, info->y, content); // add fibre content to mapped channel signal + SignalSfersID->Fill(info->fersId, content / nentries); + barX_sci += info->x * content; + barY_sci += info->y * content; + sciPosWeight += content; + + // accumulate the signal in the corresponding channel for smearing and activation steps + if (info->boardID >= 0 && info->fersId >= 0) { + // create the key to identify the FERS this SiPM belongs to + const FersKey key = make_fers_key(info->boardID, info->fersId); + + // Apply saturation consistently for both threshold and output energy + const double saturatedSignal = + ApplySaturation + ? apply_sipm_saturation_from_pe(rawPhe, sciPheGeV, NofSipmCells_sci) + : content; + + // accumulate the signal for this channel + // reads threshold values for this channel from the map + // and updates the signal and FERS->channel count maps + AccumulateChannelSignal( + sciChannelSignals, + *info, + saturatedSignal + ); + + if (writeEventDisplay) { + eventDisplaySciHits.push_back({key, info->x, info->y, rawPhe, content}); + } + } + } + } + + for (unsigned int n = 0; n < CSiPM->size(); n++) { + const double rawPhe = CSiPM->at(n); + const double content = rawPhe / cerPheGeV; + + totcer += content; + + const unsigned int towID = static_cast(n / (NofFiberscolumn * NofFibersrow / 2)); + const unsigned int SiPMID = n % (NofFiberscolumn * NofFibersrow / 2); + //const unsigned int colID = NofFiberscolumn - static_cast(SiPMID / (NofFibersrow / 2)); + const unsigned int colID = NofFiberscolumn -1U - static_cast(SiPMID / (NofFibersrow / 2)); + const unsigned int rowID = 2 * static_cast(SiPMID % (NofFibersrow / 2)) + 1; + + SipmMapC->Fill(colID, towID * NofFibersrow + rowID, content); + + const std::string tbTower = simTower_to_tbTower(std::to_string(towID)); + const int groupedCol = grouped_column(static_cast(colID), grouping); + + const SipmMapEntry* info = + FindSipmInfo(sipmLookup, tbTower, "C", static_cast(rowID), groupedCol, grouping); + + if (info) { + CerSiPMCoordinates->Fill(info->x, info->y, content); + SignalCfersID->Fill(info->fersId, content / nentries); + barX_cer += info->x * content; + barY_cer += info->y * content; + cerPosWeight += content; + + // accumulate the signal for this channel + if (info->boardID >= 0 && info->fersId >= 0) { + const FersKey key = make_fers_key(info->boardID, info->fersId); + + // Apply saturation consistently for both threshold and output energy + const double saturatedSignal = + ApplySaturation + ? apply_sipm_saturation_from_pe(rawPhe, cerPheGeV, NofSipmCells_cer) + : content; + + AccumulateChannelSignal( + cerChannelSignals, + *info, + saturatedSignal + ); + + if (writeEventDisplay) { + eventDisplayCerHits.push_back({key, info->x, info->y, rawPhe, content}); + } + } + } + } // end C Sipm loop + + + //------------------------------ + // FERS ACTIVATION LOGIC + //------------------------------ + + // Now apply noise smearing and thresholding to the accumulated channel signals + // and fill the activated channel output array for this event + // Count number of channels over threshold per FERS + SmearAndCountChannelsOverThreshold( // S channel + sciChannelSignals, + fersChannelsOverThreshold, + NoiseMode, + channelNoiseGenerators, + fersNoiseGenerators, + fersEventNoiseCache + ); + + SmearAndCountChannelsOverThreshold( // C channel + cerChannelSignals, + fersChannelsOverThreshold, + NoiseMode, + channelNoiseGenerators, + fersNoiseGenerators, + fersEventNoiseCache + ); + + for (unsigned int fersId = 1; fersId <= TotalFers; ++fersId) { + unsigned int nOverThreshold = 0; + + for (const auto& [fersKey, count] : fersChannelsOverThreshold) { + if (static_cast(static_cast(fersKey)) == fersId) { + nOverThreshold += count; + } + } + + ActivatedChannelsPerFers->Fill(fersId, nOverThreshold); + } + + + // For activated FERS channels, subtract the pedestal median + SubtractPedestalMedianFromActivatedChannels( + sciChannelSignals, + fersChannelsOverThreshold + ); + SubtractPedestalMedianFromActivatedChannels( + cerChannelSignals, + fersChannelsOverThreshold + ); + + FillActivatedChannelOutputArrays( + sciChannelSignals, + fersChannelsOverThreshold, + activatedChannelOutput, + activatedChannelOutputTrue, + activatedChannelOutputSmeared + ); + + FillActivatedChannelOutputArrays( + cerChannelSignals, + fersChannelsOverThreshold, + activatedChannelOutput, + activatedChannelOutputTrue, + activatedChannelOutputSmeared + ); + + activatedChannelsTree->Fill(); + + + PrintPerEventFersChannelCounts( + i, + sciChannelSignals, + cerChannelSignals, + fersChannelsOverThreshold + ); + + const double totsciFersOn = + SumActivatedSmearedChannelOutput(sciChannelSignals, fersChannelsOverThreshold); + + const double totcerFersOn = + SumActivatedSmearedChannelOutput(cerChannelSignals, fersChannelsOverThreshold); + + if (sciPosWeight > 0.) { + ResidualHistSciX->Fill((barX_sci / sciPosWeight) - beamX); + ResidualHistSciY->Fill((barY_sci / sciPosWeight) - beamY); + } + + if (cerPosWeight > 0.) { + ResidualHistCerX->Fill((barX_cer / cerPosWeight) - beamX); + ResidualHistCerY->Fill((barY_cer / cerPosWeight) - beamY); + } + + if (writeEventDisplay && eventDisplayDir) { + eventDisplayDir->cd(); + + std::ostringstream sciTitle; + sciTitle << "Sci event display, entry " << i + << " (grouped channels, activated FERS only); X [mm]; Y [mm]"; + TH2F* sciEventDisplay = CreateEventDisplayHist( + make_event_display_name("EventDisplaySci", i), + sciTitle.str(), + grouping + ); + + std::ostringstream cerTitle; + cerTitle << "Cer event display, entry " << i + << " (grouped channels, activated FERS only); X [mm]; Y [mm]"; + TH2F* cerEventDisplay = CreateEventDisplayHist( + make_event_display_name("EventDisplayCer", i), + cerTitle.str(), + grouping + ); + + FillActivatedEventDisplay( + sciEventDisplay, + eventDisplaySciHits, + fersChannelsOverThreshold + ); + FillActivatedEventDisplay( + cerEventDisplay, + eventDisplayCerHits, + fersChannelsOverThreshold + ); + + sciEventDisplay->Write(); + cerEventDisplay->Write(); + + delete sciEventDisplay; + delete cerEventDisplay; + + f.cd(); + } + + + + scieneFersOn->Fill(totsciFersOn); + cereneFersOn->Fill(totcerFersOn); + + sciene->Fill(totsci); + cerene->Fill(totcer); + totene->Fill(0.5 * (totsci + totcer)); + totenec->Fill(picont * (totsci - chi * totcer) / (1 - chi)); + totdep->Fill(tottow / 1000.); + leakene->Fill(lenergy / 1000. / energy); + chidist->Fill((totsci - ecalo) / (totcer - ecalo)); + + + h2_SvsC->Fill(totsci / energy, totcer / energy); + h2_SvsCFersOn->Fill(totsciFersOn / energy, totcerFersOn / energy); + + + //break; + } + + const GaussianFitSummary sFit = FitGaussianSummary(sciene); + const GaussianFitSummary cFit = FitGaussianSummary(cerene); + const GaussianFitSummary combFit = FitGaussianSummary(totene); + const GaussianFitSummary combChiFit = FitGaussianSummary(totenec); + + const HistMomentSummary resSciX = GetHistMomentSummary(ResidualHistSciX); + const HistMomentSummary resSciY = GetHistMomentSummary(ResidualHistSciY); + const HistMomentSummary resCerX = GetHistMomentSummary(ResidualHistCerX); + const HistMomentSummary resCerY = GetHistMomentSummary(ResidualHistCerY); + + const GaussianFitSummary sFitFersOn = FitGaussianSummary(scieneFersOn); + const GaussianFitSummary cFitFersOn = FitGaussianSummary(cereneFersOn); + + try { + AppendSummaryCsvLocked("hidra_summary.csv", + input, + energy, + nentries, + sFit, + cFit, + combFit, + combChiFit, + resSciX, + resSciY, + resCerX, + resCerY, + sFitFersOn, + cFitFersOn); + } catch (const std::exception& e) { + std::cerr << e.what() << std::endl; + } + + f.Write(); + f.Close(); + + simfile->Close(); + delete simfile; +} diff --git a/analysis/tb25_analysis/README.md b/analysis/tb25_analysis/README.md new file mode 100644 index 0000000..6bcbace --- /dev/null +++ b/analysis/tb25_analysis/README.md @@ -0,0 +1,81 @@ +# TB25 Analysis + +This directory contains the TB25-specific ROOT analysis macro, helper scripts, +and CSV mapping inputs. + +Main files: + +- `HidraTB25Ana.C`: TB25 analysis macro. +- `run_ana.sh`: helper script to run the macro over the configured + energy/input-file list. +- `plot_hidra.py`: plotting helper for the TB25 summary outputs. +- `*_sipm_std_run868_HG_4fers.csv`: TB25 SiPM CSV mapping files. + +## Workflow + +The TB25 workflow is: + +1. Open the input HidraSim ROOT file and read the simulated PMT and SiPM signal + vectors. +2. Load the SiPM hardware/channel map, high-gain pedestal/noise values, and + ADC-to-GeV calibration constants from JSON files. +3. Convert simulation SiPM indices to the TB25 tower, row, and grouped-column + convention used by the FERS map. +4. Associate each simulated SiPM signal to a mapped FERS channel. +5. Optionally apply SiPM saturation, pedestal terms, and electronics noise. +6. Count channels above threshold and mark a FERS as active when the configured + multiplicity requirement is satisfied. +7. Fill both all-signal histograms and FERS-activated histograms, then append + fit and resolution summaries to `hidra_summary.csv`. + +## Configuration + +The main behavior switches are defined near the top of `HidraTB25Ana.C`: + +- `ApplySaturation`: apply or skip SiPM finite-cell saturation. +- `ApplyTbNoise`: include or exclude the TB pedestal/noise smearing. +- `NoiseMode`: choose `UncorrelatedByChannel` or `CorrelatedWithinFers`. +- `NoiseDistribution`: choose `Gaussian` or `LogNormal` noise sampling. +- `AddPedestalToEnergyContribution`: add the pedestal contribution to the + channel signal before threshold/output evaluation. +- `ApplyPedestalSubtraction`: subtract the pedestal median from activated + channels. +- `ScalePedestalSubtractionFactorS` and `ScalePedestalSubtractionFactorC`: + scale the scintillation and Cherenkov pedestal subtraction. +- `FersThresholdScaleFactor`: scale the FERS activation thresholds. +- `FersMultiplicity`: number of channels above threshold required to activate a + FERS. +- `EventDisplayEvery`: event-display writing interval; set to `0` to disable. +- `printSmearingLog` and `printFersLog`: enable detailed debug prints. + +## Inputs + +The macro includes the shared detector geometry from the parent analysis +directory: + +```cpp +#include "../HidraGeo.h" +``` + +It reads JSON mapping and calibration inputs from: + +```text +../../../TBDataPreparation/2025_SPS/MapAndCalibration/ +``` + +when run from `HidraSim/analysis/tb25_analysis/`. These JSON files are kept with +the TB data-preparation inputs rather than duplicated here. + +The local CSV files in this directory are TB25-specific SiPM mapping inputs. + +## Usage + +Typical usage: + +```bash +cd HidraSim/analysis/tb25_analysis +./run_ana.sh +``` + +The script writes `hidra_summary.csv`, per-energy ROOT files, and run logs in +this directory. These are generated outputs and should normally not be committed. diff --git a/analysis/tb25_analysis/S_sipm_std_run868_HG_4fers.csv b/analysis/tb25_analysis/S_sipm_std_run868_HG_4fers.csv new file mode 100644 index 0000000..75da715 --- /dev/null +++ b/analysis/tb25_analysis/S_sipm_std_run868_HG_4fers.csv @@ -0,0 +1,512 @@ +2.305711819222606290e-02 +1.513391876707899895e-02 +1.322797414971654643e-02 +1.748574847610178534e-02 +2.473508375359059416e-02 +1.772337411965317633e-02 +2.685229938967816088e-02 +1.894335822301219940e-02 +2.646827908244578711e-02 +1.787763218381282701e-02 +1.765511576221477641e-02 +1.941687866630632808e-02 +2.132065844534831653e-02 +6.908711114226655237e-02 +2.397455499668809178e-02 +2.149104328600777783e-02 +2.807175448417143643e-02 +1.815141692784673702e-02 +1.502438916659084786e-02 +1.628697516341394771e-02 +2.279128628388118680e-02 +1.777556194621867855e-02 +2.500008995263142825e-02 +1.951038065719794615e-02 +1.334707227201344971e-02 +8.542301152815205342e-03 +8.077926783646978068e-03 +1.199056140858039202e-02 +9.913152847711980201e-03 +9.841874027004374886e-03 +9.788394072587879025e-03 +5.553031634572793557e-03 +4.310381431361559601e-02 +2.843450274843038750e-02 +1.554175003535850606e-02 +1.689996879994793214e-02 +1.571916021678879116e-02 +1.460279913532018699e-02 +1.601764863640543421e-02 +2.663592482500250561e-02 +1.989546894320525991e-02 +1.618380023731210565e-02 +2.605065710227762524e-02 +1.802501911265464232e-02 +1.585062338886101402e-02 +1.835630905533446738e-02 +1.628331667416600681e-02 +1.557137457250674928e-02 +1.394487818303505149e-02 +7.207567989508252039e-03 +7.213063611096495048e-03 +6.871005545480772872e-03 +6.533379766360193804e-03 +7.820471267818172867e-03 +6.528948933004081676e-03 +9.210675781694387804e-03 +1.293024013243798676e-02 +8.958105140251620274e-03 +8.651857544665630950e-03 +8.008420667889329381e-03 +8.414976108832803239e-03 +6.518481943235037400e-03 +6.637519791023680296e-03 +8.904482670210612724e-03 +2.117591332194612014e-02 +2.278314971128930083e-02 +2.074859796508380164e-02 +2.479078573245764930e-02 +1.887718399186507495e-02 +2.362071863464529006e-02 +1.845152149818079720e-02 +2.234567508403728989e-02 +1.522235494533204746e-02 +1.370621155651514760e-02 +1.359882856972681291e-02 +1.406114323102437488e-02 +1.340202500507947994e-02 +1.282744614739213945e-02 +1.296867788808374507e-02 +1.342790410117034386e-02 +1.662191803169347018e-02 +1.644534155074779680e-02 +1.967452557474560007e-02 +1.663856020587760332e-02 +1.264382552607514268e-02 +1.751217414914104167e-02 +1.045308516007653794e-02 +1.780369082117239793e-02 +9.995692024969641404e-03 +9.801574135469311513e-03 +9.226870815066854151e-03 +8.292994762413628471e-03 +7.707003426454046152e-03 +8.738835598011097591e-03 +7.271192199666008660e-03 +9.995663233697028860e-03 +1.983369513105603943e-02 +1.964245289136174605e-02 +1.351460851184235236e-02 +1.798723647318080376e-02 +1.900846002355956882e-02 +2.110259496272508251e-02 +1.755056493698323519e-02 +2.650439796531516012e-02 +1.827810536973928021e-02 +1.371160058529192691e-02 +1.817006888319001173e-02 +1.787794187567192869e-02 +2.038027302066183072e-02 +1.918368144521478968e-02 +1.647609314320732168e-02 +2.584813168879387252e-02 +1.378231225131894104e-02 +1.249301140951730379e-02 +1.209633864143224054e-02 +1.485955938826596066e-02 +1.081625277141968350e-02 +1.059905995974842405e-02 +8.833176474664543656e-03 +1.246949350409607901e-02 +1.766366072298494044e-02 +1.427840135036185024e-02 +8.221649985049228962e-03 +7.599065416730870631e-03 +6.912038162458175941e-03 +8.342663149717318566e-03 +6.989825588642209486e-03 +7.037982400303029358e-03 +2.412477138406342661e-02 +1.739111432728398382e-02 +1.581495789607646058e-02 +1.602875730825044018e-02 +1.912181127073478773e-02 +2.220194310585115252e-02 +2.258202309378515987e-02 +2.046572840021396861e-02 +1.492940714292235842e-02 +1.240909791738495707e-02 +1.101700913055672464e-02 +1.103376583429502289e-02 +1.124554521973124277e-02 +1.855674431465702615e-02 +1.263330854357613484e-02 +1.587274009157202403e-02 +1.335923271868529570e-02 +9.191038188588386559e-03 +1.115878861149216181e-02 +1.241219909331188996e-02 +1.387865277425220749e-02 +1.627357536395626636e-02 +1.392994436254543951e-02 +1.334619370985222660e-02 +9.909997779014448685e-03 +8.062106348138141018e-03 +7.749318672996121336e-03 +8.562968949429251159e-03 +1.147088660154488293e-02 +1.143108930740633668e-02 +1.113921647191978836e-02 +1.377319651134702345e-02 +2.475516539613116865e-02 +2.445156333010631636e-02 +2.336674636807934205e-02 +2.292646781866428465e-02 +2.517701805900934842e-02 +2.522008421061388422e-02 +2.292792404201430059e-02 +2.962390255188046204e-02 +2.700358099634503350e-02 +2.342058201670429302e-02 +1.918068878133140712e-02 +1.796453753128103223e-02 +2.137267082126124984e-02 +1.947867991954907799e-02 +2.139949906898428941e-02 +2.867160000334039935e-02 +1.086023284498679846e-02 +8.268878707467557099e-03 +8.902699695261405768e-03 +8.818291532764522406e-03 +9.185187025374774369e-03 +1.163105199025892683e-02 +1.079053128672334376e-02 +1.229092785844998828e-02 +7.056131993434780827e-03 +6.144392719391784247e-03 +6.581609389146724810e-03 +7.182904403282258654e-03 +7.732560230494253325e-03 +7.574433014469375558e-03 +8.753308008037335997e-03 +1.043332549227720014e-02 +9.895796328711846043e-03 +1.058573193643532612e-02 +6.906838296398618536e-03 +1.021787102435910838e-02 +1.002398349518047842e-02 +9.608535068042867860e-03 +1.220496738631933996e-02 +1.075756723489917910e-02 +1.238548922243091696e-02 +1.090578000018627841e-02 +1.035105430697457092e-02 +1.088182013690566080e-02 +1.117903535193046256e-02 +9.514469851133302084e-03 +1.239855057815660716e-02 +1.328643477543714480e-02 +8.277842858585689190e-03 +8.441245876613621968e-03 +6.861916256192850810e-03 +6.633276995921845443e-03 +6.576223564555145652e-03 +6.745101686861997708e-03 +1.005192023049534394e-02 +8.968547990593846858e-03 +6.661154235134736194e-03 +6.010436009689256419e-03 +5.748889139891848188e-03 +6.004298793773534676e-03 +5.584439924888111480e-03 +6.108141821150703948e-03 +7.346114779663830835e-03 +6.513897460875326881e-03 +1.853227255488703457e-02 +1.488299669090168906e-02 +1.692569867393687449e-02 +1.553095516023080339e-02 +1.721873821372871954e-02 +1.881267836869665935e-02 +1.909896521789889237e-02 +1.692722250442919157e-02 +1.137719604358441994e-02 +1.190873052123977982e-02 +1.446966871241184456e-02 +1.083470457121271137e-02 +1.150397091896960926e-02 +1.124537632248278317e-02 +1.166837607623026109e-02 +1.021950675111312455e-02 +8.168476327853228108e-03 +7.765799276999278862e-03 +7.843921777191786848e-03 +7.979609477286618568e-03 +8.090189575124093549e-03 +6.520339764218185258e-03 +6.994684470905278457e-03 +8.613730433064301431e-03 +7.537929578054690556e-03 +5.955814618447395105e-03 +6.074847374171108566e-03 +4.727209091953086216e-03 +4.710738375629699125e-03 +4.765086801794133665e-03 +5.751658816097853129e-03 +5.013253946602004747e-03 +6.380107708887659823e-03 +6.283077365514594913e-03 +6.411662995805026376e-03 +4.948500624741559474e-03 +6.482745164347274332e-03 +6.530409306650171442e-03 +6.528005871575447150e-03 +8.184673870362146886e-03 +5.543755283669403108e-03 +5.263487847102016857e-03 +5.602470835259222050e-03 +5.801557075547157660e-03 +5.681603317773607255e-03 +5.894059298559701371e-03 +6.176153705498186978e-03 +7.931601935006863530e-03 +5.660088815580100602e-03 +5.131665230867867022e-03 +5.110726205349474735e-03 +6.676013825940165826e-03 +6.444071805365590190e-03 +5.984282067307487331e-03 +5.993471381568860057e-03 +6.670018839776525009e-03 +3.782599880564759334e-03 +3.974237473048423605e-03 +3.662691509877761577e-03 +4.027108426257884621e-03 +4.303781082125520373e-03 +4.480871052183140607e-03 +3.965431974044360958e-03 +4.392286963650139306e-03 +6.184622764178267403e-03 +6.013513912279356069e-03 +6.021792350829145886e-03 +5.528743518603385088e-03 +6.699758456027955951e-03 +6.558577298769911051e-03 +5.487175701999091187e-03 +6.735153509998123393e-03 +6.391699633148381632e-03 +6.244296604960549411e-03 +7.006589459966516831e-03 +7.242363654527707133e-03 +7.953130394460852098e-03 +7.110639713351774548e-03 +7.588193530517970721e-03 +6.756346361173008915e-03 +5.023556395088599191e-03 +4.044902617698074507e-03 +4.677528392349453515e-03 +4.415027332990929289e-03 +4.840107596119327255e-03 +3.984419629017027162e-03 +3.911790997471015566e-03 +4.814894850238894165e-03 +4.104739731557826080e-03 +3.473353092460602747e-03 +3.435840228656407511e-03 +2.814826107128193776e-03 +3.372026298662513747e-03 +3.752924435012840784e-03 +3.424910641907093212e-03 +3.781358053167990264e-03 +1.049836435239131713e-02 +8.428537639687884095e-03 +7.453954381164024437e-03 +8.715443186412045323e-03 +1.300400913042182478e-02 +1.370761546869756799e-02 +1.144888661231404856e-02 +1.339809186027742308e-02 +7.931992892653498550e-03 +8.759571799394794134e-03 +6.552735388859336728e-03 +7.504874440129137049e-03 +1.056015135836597431e-02 +1.039188274675288855e-02 +1.029214812752084554e-02 +1.031189104172789903e-02 +6.832993277748782186e-03 +6.680629073738941597e-03 +5.174307447200923379e-03 +7.237762807931391942e-03 +9.079897798787402510e-03 +8.857693808180152825e-03 +9.209236159901004748e-03 +7.864511764586626028e-03 +4.790264486035778271e-03 +4.700186863059992691e-03 +3.252935802239306437e-03 +4.879036846275416721e-03 +5.270850498356439984e-03 +5.518182576382086972e-03 +5.680978346061065810e-03 +5.424155258726289473e-03 +9.415126518705100797e-03 +8.825515751642117393e-03 +7.984000949793573429e-03 +8.520956943976336004e-03 +1.411807389219540618e-02 +1.203657348394869368e-02 +1.637816789355909380e-02 +1.489666080017368333e-02 +6.274639237867831593e-03 +6.192008121624848795e-03 +6.758886231487167516e-03 +6.649544199006813920e-03 +1.311975886409485689e-02 +1.040616460643882871e-02 +1.120060211604187450e-02 +1.021946667001645027e-02 +7.305653972644995421e-03 +6.776017918032531569e-03 +6.488083403127588127e-03 +6.367969795907413993e-03 +1.055364555422462318e-02 +1.014940123207166804e-02 +8.302771089180559227e-03 +1.108019900149887806e-02 +5.157919413158772344e-03 +5.115367280213827840e-03 +5.517694954183883019e-03 +6.124937113565938340e-03 +7.248984045825345238e-03 +6.490781536198100832e-03 +6.862011607353920087e-03 +8.074574499057655605e-03 +1.284168217096736710e-02 +1.300156896997084130e-02 +1.203814371077377779e-02 +1.261696251048423367e-02 +1.260275800811621466e-02 +1.323940657153829109e-02 +1.422902108923636277e-02 +1.351971305379055625e-02 +9.404559215329774680e-03 +8.449114240700356104e-03 +8.290211801694856789e-03 +1.179531746489938249e-02 +9.749861358929667529e-03 +1.093203720664497579e-02 +9.322632623168467317e-03 +1.155455776087093399e-02 +7.928931948303501870e-03 +6.384816734770768290e-03 +6.655391804855696064e-03 +5.184129660965529698e-03 +5.155646370577120913e-03 +4.712463283277314868e-03 +5.184479305140322382e-03 +5.003169928967682312e-03 +4.856584900374296548e-03 +4.499122001822833915e-03 +3.755274734636015500e-03 +3.265351224668859575e-03 +3.260370710549323654e-03 +3.146172437359927385e-03 +3.272941988945381013e-03 +3.104389511305521293e-03 +1.672253041174403751e-02 +9.674226930046858283e-03 +8.647479381133556589e-03 +7.897325814738049773e-03 +7.758161384118163022e-03 +7.452523704086221577e-03 +6.353117153714962773e-03 +5.783044120893023228e-03 +1.027319720109367077e-02 +7.750207318899597977e-03 +2.396833115782897455e-02 +1.754324814529978696e-02 +6.410710113099798643e-03 +5.932994326603271799e-03 +5.569931647048666570e-03 +5.085573701463923756e-03 +7.325792867343913946e-03 +6.919951150376287517e-03 +2.779127494487095007e-02 +2.623649214513916819e-02 +6.438175890028154449e-03 +5.118435739067504417e-03 +4.294310852008550485e-03 +4.297854342555858069e-03 +5.745901178233893229e-03 +4.944807916624286738e-03 +4.758350133728087364e-03 +4.257194911016243082e-03 +3.949824724384824910e-03 +3.814810448360730129e-03 +3.495641865089077049e-03 +3.119232456837268529e-03 +4.525317918259464181e-03 +8.992474427447154361e-03 +8.731321535428141281e-03 +8.741141194956657604e-03 +8.327958075125905385e-03 +9.173955491847372856e-03 +8.682416278553792971e-03 +8.257545358401181812e-03 +6.548489272638183083e-03 +7.175652503177763590e-03 +5.543196594748405273e-03 +6.705449826937393172e-03 +6.063890376349194596e-03 +4.888938722169963377e-03 +5.590706459942880631e-03 +4.977937323488388840e-03 +5.159317183750596406e-03 +4.777580480459079661e-03 +3.754851699602590291e-03 +3.876975083104591303e-03 +3.808475529082814420e-03 +3.098876890795967780e-03 +3.337784460160160300e-03 +3.700097775799019058e-03 +2.897961981629366014e-03 +3.097121192917526560e-03 +2.865205478282426260e-03 +2.570180592701415759e-03 +2.805494301010433983e-03 +2.258558342117879085e-03 +2.673001486401912607e-03 +3.020266305424277199e-03 +9.641925193708532696e-03 +7.892657166720188841e-03 +6.472612719408980288e-03 +7.033539413702497717e-03 +5.702659558842836270e-03 +5.497026834873058182e-03 +4.937383503280632496e-03 +5.445138550447543857e-03 +5.276092712300537585e-03 +5.195695760332323318e-03 +4.468370981774457025e-03 +4.276701333630262246e-03 +4.840967736026589342e-03 +5.455757452292591715e-03 +4.302742683086949009e-03 +4.619187514228566271e-03 +6.569277213095098093e-03 +6.374616240988035566e-03 +5.730342855778345137e-03 +6.442361865297583333e-03 +5.317987871509457197e-03 +5.305939148522820746e-03 +3.812425052652455675e-03 +4.219485726344291386e-03 +5.537800871496355824e-03 +4.698651454016443234e-03 +4.616501262338533994e-03 +3.993114683439616397e-03 +3.851139679999344544e-03 +3.439144424124739984e-03 +3.622062453994871322e-03 +3.781768839251361571e-03 diff --git a/analysis/tb25_analysis/plot_hidra.py b/analysis/tb25_analysis/plot_hidra.py new file mode 100644 index 0000000..7644565 --- /dev/null +++ b/analysis/tb25_analysis/plot_hidra.py @@ -0,0 +1,529 @@ +# file: plot_hidra_summary.py +import os +import sys +from array import array + +import numpy as np +import pandas as pd +import ROOT + + +ROOT.gROOT.SetBatch(True) +ROOT.gStyle.SetOptStat(0) + +scaling_s = 1.0 +scaling_c = 1.0 + +#scaling_s = 1.1 +#scaling_c = 1.1 + +#scaling_s = 1.14 +#scaling_c = 1.27 + +#scaling_s = 1.16 +#scaling_c = 1.2 +#scaling_s = 1.12 +#scaling_c = 1.18 + +#scaling_s = 1.13 +#scaling_c = 1.14 + +myOutDir = "plots/Attenuation10m" +#myOutDir = "plots/AttenuationTB24" + +extraTex2 = "Fibre attenuation set to 10m" +#extraTex2 = "Fibre attenuation tuned to TB24 data: 3.7 m (Sci), 3.9 m (Cer)" + +#extraTex1 = "" +#extraTex1 = "Fully correlated gaussian noise between FERS" +extraTex1 = "Uncorrelated gaussian noise" +#extraTex1 = "Fully correlated log-normal noise between FERS" +#extraTex1 = "Uncorrelated log-normal noise" +#extraTex1 = "TB25 HG data used for tuning channel noise" +#extraTex1 = "FERS Thresholds without pedestal subtraction" +#extraTex1 = "FERS Thresholds pedestal-subtracted" + +#extraTex5 = " " +extraTex5 = "Correction for multiple photons on the same SiPM cell: ON" + +extraTex4 = " " +#extraTex4 = "FERS Thresholds without pedestal subtraction" +extraTex4 = "FERS Thresholds pedestal-subtracted" +#extraTex4 = "Pedestal subtraction: (S: #times1.045 , C: #times1.035 mean pedestal)" +#extraTex4 = "Pedestal subtraction: (S: #times1.10 , C: #times1.055 mean pedestal)" +#extraTex4 = "Pedestal subtraction: ON" +#extraTex4 = "Pedestal (ADC -> GeV) subtraction: ON (0.4x mean pedestal)" +#extraTex4 = "Poissonian probability to lose phe per channel: #mu = 3 (Sci), 3 (Cer)" + +def load_summary(csv_path: str) -> pd.DataFrame: + data = pd.read_csv(csv_path) + data = data.sort_values("truth_energy").reset_index(drop=True) + return data + + +def safe_ratio(num: np.ndarray, den: np.ndarray) -> np.ndarray: + out = np.full_like(num, np.nan, dtype=float) + mask = np.isfinite(num) & np.isfinite(den) & (den != 0.0) + out[mask] = num[mask] / den[mask] + return out + + +def safe_ratio_error( + num: np.ndarray, + num_err: np.ndarray, + den: np.ndarray, + den_err: np.ndarray, +) -> np.ndarray: + out = np.full_like(num, np.nan, dtype=float) + mask = ( + np.isfinite(num) + & np.isfinite(num_err) + & np.isfinite(den) + & np.isfinite(den_err) + & (num != 0.0) + & (den != 0.0) + ) + ratio = np.full_like(num, np.nan, dtype=float) + ratio[mask] = num[mask] / den[mask] + rel = np.full_like(num, np.nan, dtype=float) + rel[mask] = np.sqrt((num_err[mask] / num[mask]) ** 2 + (den_err[mask] / den[mask]) ** 2) + out[mask] = ratio[mask] * rel[mask] + return out + + +def to_root_arrays(x, y, ex, ey): + return ( + array("d", [float(v) for v in x]), + array("d", [float(v) for v in y]), + array("d", [float(v) for v in ex]), + array("d", [float(v) for v in ey]), + ) + + +def make_graph( + x: np.ndarray, + y: np.ndarray, + ey: np.ndarray, + color: int, + marker: int, + name: str, + line_style: int = 1, +): + ex = np.zeros_like(x, dtype=float) + mask = np.isfinite(x) & np.isfinite(y) & np.isfinite(ey) + x_use = x[mask] + y_use = y[mask] + ex_use = ex[mask] + ey_use = ey[mask] + + if len(x_use) == 0: + return None + + x_arr, y_arr, ex_arr, ey_arr = to_root_arrays(x_use, y_use, ex_use, ey_use) + graph = ROOT.TGraphErrors(len(x_use), x_arr, y_arr, ex_arr, ey_arr) + graph.SetName(name) + graph.SetMarkerColor(color) + graph.SetLineColor(color) + graph.SetMarkerStyle(marker) + graph.SetMarkerSize(1.2) + graph.SetLineWidth(2) + graph.SetLineStyle(line_style) + return graph + + +def fit_resolution_graph(graph: ROOT.TGraphErrors, color: int, fit_name: str): + if graph is None or graph.GetN() < 2: + return None, None + + fit = ROOT.TF1(fit_name, "pol1", 0.0, 1.0) + fit.SetLineColor(color) + fit.SetLineWidth(2) + result = graph.Fit(fit, "QS") + if int(result) != 0: + return fit, None + return fit, result + + +def draw_header(x_left: float, y_top: float, subtitle: str): + tex0 = ROOT.TLatex(x_left, y_top, "HidraSim") + tex0.SetNDC() + tex0.SetTextFont(72) + tex0.SetTextSize(0.042) + tex0.Draw() + + tex1 = ROOT.TLatex(x_left + 0.16, y_top, "Preliminary") + tex1.SetNDC() + tex1.SetTextFont(42) + tex1.SetTextSize(0.042) + tex1.Draw() + + tex2 = ROOT.TLatex(x_left, y_top - 0.055, subtitle) + tex2.SetNDC() + tex2.SetTextFont(42) + tex2.SetTextSize(0.032) + tex2.Draw() + + +def setup_canvas(name: str, width: int = 800, height: int = 650): + canvas = ROOT.TCanvas(name, name, width, height) + canvas.SetFillColor(0) + canvas.SetBorderMode(0) + canvas.SetBorderSize(2) + canvas.SetTickx(1) + canvas.SetTicky(1) + canvas.SetLeftMargin(0.13) + canvas.SetRightMargin(0.05) + canvas.SetTopMargin(0.09) + canvas.SetBottomMargin(0.11) + canvas.SetFrameBorderMode(0) + return canvas + + +def plot_energy_resolution(data: pd.DataFrame, out_dir: str): + energy = data["truth_energy"].to_numpy(dtype=float) + inv_sqrt_e = 1.0 / np.sqrt(energy) + + s_res = safe_ratio( + data["s_fit_rms"].to_numpy(dtype=float), + data["s_fit_mean"].to_numpy(dtype=float), + ) + s_res_err = safe_ratio_error( + data["s_fit_rms"].to_numpy(dtype=float), + data["s_fit_rms_err"].to_numpy(dtype=float), + data["s_fit_mean"].to_numpy(dtype=float), + data["s_fit_mean_err"].to_numpy(dtype=float), + ) + + c_res = safe_ratio( + data["c_fit_rms"].to_numpy(dtype=float), + data["c_fit_mean"].to_numpy(dtype=float), + ) + c_res_err = safe_ratio_error( + data["c_fit_rms"].to_numpy(dtype=float), + data["c_fit_rms_err"].to_numpy(dtype=float), + data["c_fit_mean"].to_numpy(dtype=float), + data["c_fit_mean_err"].to_numpy(dtype=float), + ) + + comb_res = safe_ratio( + data["comb_fit_rms"].to_numpy(dtype=float), + data["comb_fit_mean"].to_numpy(dtype=float), + ) + comb_res_err = safe_ratio_error( + data["comb_fit_rms"].to_numpy(dtype=float), + data["comb_fit_rms_err"].to_numpy(dtype=float), + data["comb_fit_mean"].to_numpy(dtype=float), + data["comb_fit_mean_err"].to_numpy(dtype=float), + ) + + canvas = setup_canvas("c_resolution") + mg = ROOT.TMultiGraph() + mg.SetTitle("Energy resolution;1/#sqrt{E_{beam}} [GeV^{-1/2}];#sigma(E)/E") + + graph_comb = make_graph(inv_sqrt_e, comb_res, comb_res_err, ROOT.kGreen + 2, 20, "graph_comb") + graph_s = make_graph(inv_sqrt_e, s_res, s_res_err, ROOT.kRed + 1, 21, "graph_s") + graph_c = make_graph(inv_sqrt_e, c_res, c_res_err, ROOT.kBlue + 1, 22, "graph_c") + + for graph in (graph_comb, graph_s, graph_c): + if graph is not None: + mg.Add(graph, "PE") + + mg.Draw("A") + mg.GetXaxis().SetTitleOffset(1.15) + mg.GetYaxis().SetTitleOffset(1.35) + + fit_comb, _ = fit_resolution_graph(graph_comb, ROOT.kGreen + 2, "fit_comb") + fit_s, _ = fit_resolution_graph(graph_s, ROOT.kRed + 1, "fit_s") + fit_c, _ = fit_resolution_graph(graph_c, ROOT.kBlue + 1, "fit_c") + + + alignLeft = 0.18 + alignTop = 0.85 + alignRight = 0.925 + tex0 = ROOT.TLatex(alignLeft, alignTop, "HidraSim") + tex0.SetNDC() + tex0.SetTextFont(72); + tex0.SetTextSize(0.042); + tex0.SetLineWidth(2); + tex0.Draw("same"); + #tex1 = ROOT.TLatex(alignLeft+0.16,alignTop,"Preliminary "); + tex1 = ROOT.TLatex(alignLeft+0.16,alignTop,"Work in progress"); + tex1.SetNDC(); + tex1.SetTextFont(42); + tex1.SetTextSize(0.042); + tex1.SetLineWidth(2); + tex1.Draw("same"); + + legend = ROOT.TLegend(0.16, 0.57, 0.72, 0.82) + legend.SetBorderSize(0) + legend.SetFillStyle(0) + legend.SetTextFont(42) + legend.SetTextSize(0.03) + + if fit_comb is not None: + legend.AddEntry( + fit_comb, + "Combined: #sigma/E = %.2f%%/#sqrt{E} %+ .2f%%" + % (100.0 * fit_comb.GetParameter(1), 100.0 * fit_comb.GetParameter(0)), + "l", + ) + if fit_s is not None: + legend.AddEntry( + fit_s, + "S: #sigma/E = %.2f%%/#sqrt{E} %+ .2f%%" + % (100.0 * fit_s.GetParameter(1), 100.0 * fit_s.GetParameter(0)), + "l", + ) + if fit_c is not None: + legend.AddEntry( + fit_c, + "C: #sigma/E = %.2f%%/#sqrt{E} %+ .2f%%" + % (100.0 * fit_c.GetParameter(1), 100.0 * fit_c.GetParameter(0)), + "l", + ) + + legend.Draw() + draw_header(0.16, 0.84, "Energy resolution from Gaussian fits") + canvas.SaveAs(os.path.join(out_dir, "energy_resolution.pdf")) + + +def plot_energy_linearity(data: pd.DataFrame, out_dir: str): + energy = data["truth_energy"].to_numpy(dtype=float) + + s_mean = data["s_fit_mean"].to_numpy(dtype=float) + s_mean_err = data["s_fit_mean_err"].to_numpy(dtype=float) + c_mean = data["c_fit_mean"].to_numpy(dtype=float) + c_mean_err = data["c_fit_mean_err"].to_numpy(dtype=float) + comb_mean = data["comb_fit_mean"].to_numpy(dtype=float) + comb_mean_err = data["comb_fit_mean_err"].to_numpy(dtype=float) + + s_fers_on_mean = data["s_fit_fers_on_mean"].to_numpy(dtype=float) + s_fers_on_mean_err = data["s_fit_fers_on_mean_err"].to_numpy(dtype=float) + c_fers_on_mean = data["c_fit_fers_on_mean"].to_numpy(dtype=float) + c_fers_on_mean_err = data["c_fit_fers_on_mean_err"].to_numpy(dtype=float) + + #s_lin = safe_ratio(s_mean - energy, energy) + s_lin = safe_ratio(s_mean, energy) + s_lin_err = safe_ratio(s_mean_err, energy) + + #c_lin = safe_ratio(c_mean - energy, energy) + c_lin = safe_ratio(c_mean, energy) + c_lin_err = safe_ratio(c_mean_err, energy) + + #comb_lin = safe_ratio(comb_mean - energy, energy) + comb_lin = safe_ratio(comb_mean, energy) + comb_lin_err = safe_ratio(comb_mean_err, energy) + + #s_fers_on_lin = safe_ratio(s_fers_on_mean - energy, energy) + s_fers_on_lin = safe_ratio(s_fers_on_mean*scaling_s, energy) + s_fers_on_lin_err = safe_ratio(s_fers_on_mean_err, energy) + + #c_fers_on_lin = safe_ratio(c_fers_on_mean - energy, energy) + c_fers_on_lin = safe_ratio(c_fers_on_mean*scaling_c, energy) + c_fers_on_lin_err = safe_ratio(c_fers_on_mean_err, energy) + + canvas = setup_canvas("c_linearity") + mg = ROOT.TMultiGraph() + #mg.SetTitle("Energy linearity;E_{beam} [GeV];(E_{fit}-E_{beam})/E_{beam}") + mg.SetTitle("Energy linearity;E_{beam} [GeV];E_{fit} / E_{beam}") + + graph_comb = make_graph(energy, comb_lin, comb_lin_err, ROOT.kGreen + 2, 20, "lin_comb") + graph_s = make_graph(energy, s_lin, s_lin_err, ROOT.kRed + 1, 21, "lin_s") + graph_c = make_graph(energy, c_lin, c_lin_err, ROOT.kBlue + 1, 22, "lin_c") + graph_s_fers_on = make_graph( + energy, + s_fers_on_lin, + s_fers_on_lin_err, + ROOT.kRed + 1, + 25, + "lin_s_fers_on", + line_style=1, + ) + graph_c_fers_on = make_graph( + energy, + c_fers_on_lin, + c_fers_on_lin_err, + ROOT.kBlue + 1, + 26, + "lin_c_fers_on", + line_style=1, + ) + + #for graph in (graph_comb, graph_s, graph_c, graph_s_fers_on, graph_c_fers_on): + for graph in (graph_comb, graph_s, graph_c): + if graph is not None: + mg.Add(graph, "PE") + + mg.Add(graph_s_fers_on, "PEL") + mg.Add(graph_c_fers_on, "PEL") + + mg.Draw("A") + mg.GetXaxis().SetTitleOffset(1.07) + mg.GetYaxis().SetTitleOffset(1.2) + #mg.SetMinimum(-0.40) + #mg.SetMaximum(0.08) + + mg.SetMinimum(+0.9) + mg.SetMaximum(+1.14) + + xmin = float(np.nanmin(energy)) + xmax = float(np.nanmax(energy)) + #for y_val, color in ((0.0, ROOT.kGreen + 3), (0.01, ROOT.kGray + 2), (-0.01, ROOT.kGray + 2)): + # line = ROOT.TLine(xmin, y_val, xmax, y_val) + # line.SetLineStyle(2) + # line.SetLineColor(color) + # line.Draw() + + + line1=ROOT.TLine(mg.GetXaxis().GetXmin(), 1.01, mg.GetXaxis().GetXmax(), 1.01) + line2=ROOT.TLine(mg.GetXaxis().GetXmin(), 0.99, mg.GetXaxis().GetXmax(), 0.99) + line3=ROOT.TLine(mg.GetXaxis().GetXmin(), 1.0, mg.GetXaxis().GetXmax(), 1.0) + line1.SetLineStyle(2) + line2.SetLineStyle(2) + line3.SetLineStyle(2) + line3.SetLineColor(ROOT.kGreen+3) + + line1.Draw("same") + line2.Draw("same") + line3.Draw("same") + line2.Draw("same") + line3.Draw("same") + + alignLeft = 0.18 + alignTop = 0.85 + alignRight = 0.925 + tex0 = ROOT.TLatex(alignLeft, alignTop, "HidraSim") + tex0.SetNDC() + tex0.SetTextFont(72); + tex0.SetTextSize(0.042); + tex0.SetLineWidth(2); + tex0.Draw("same"); + #tex1 = ROOT.TLatex(alignLeft+0.16,alignTop,"Preliminary "); + tex1 = ROOT.TLatex(alignLeft+0.16,alignTop,"Work in progress"); + tex1.SetNDC(); + tex1.SetTextFont(42); + tex1.SetTextSize(0.042); + tex1.SetLineWidth(2); + tex1.Draw("same"); + tex2 = ROOT.TLatex(alignLeft,alignTop-0.12, extraTex2) + tex2.SetNDC(); + tex2.SetTextFont(42); + tex2.SetTextSize(0.03); + tex2.SetLineWidth(2); + tex2.Draw("same"); + tex3 = ROOT.TLatex(alignLeft,alignTop-0.08, extraTex1) + tex3.SetNDC(); + tex3.SetTextFont(42); + tex3.SetTextSize(0.03); + tex3.SetLineWidth(2); + tex3.Draw("same"); + tex4 = ROOT.TLatex(alignLeft,alignTop-0.04, extraTex4) + tex4.SetNDC(); + tex4.SetTextFont(42); + tex4.SetTextSize(0.03); + tex4.SetLineWidth(2); + tex4.Draw("same"); + tex5 = ROOT.TLatex(alignLeft,alignTop-0.16, extraTex5) + tex5.SetNDC(); + tex5.SetTextFont(42); + tex5.SetTextSize(0.03); + tex5.SetLineWidth(2); + tex5.Draw("same"); + + + legend = ROOT.TLegend(0.450, 0.17, 0.83, 0.37) + legend.SetBorderSize(0) + legend.SetFillStyle(0) + legend.SetTextFont(42) + legend.SetTextSize(0.03) + #legend.SetHeader("Fers thresholds tuned to TB data") # Set the header for the legend + if graph_comb is not None: + legend.AddEntry(graph_comb, "Combined - all fibres", "pe") + if graph_s is not None: + legend.AddEntry(graph_s, "S - all fibres", "pe") + if graph_c is not None: + legend.AddEntry(graph_c, "C - all fibres", "pe") + if graph_s_fers_on is not None: + #legend.AddEntry(graph_s_fers_on, "S - FERS On only, Thr: >= 1 phe #times 1.1", "pe") + legend.AddEntry(graph_s_fers_on, "S - FERS On only #times " + str(scaling_s), "ple") + if graph_c_fers_on is not None: + #legend.AddEntry(graph_c_fers_on, "C - FERS On only, Thr: >= 1 phe #times 1.1", "pe") + legend.AddEntry(graph_c_fers_on, "C - FERS On only #times " + str(scaling_c), "ple") + legend.Draw() + + draw_header(0.16, 0.84, "Linearity from fitted energy means") + canvas.SaveAs(os.path.join(out_dir, "energy_linearity.pdf")) + + +def plot_spatial_resolution(data: pd.DataFrame, out_dir: str): + energy = data["truth_energy"].to_numpy(dtype=float) + + configs = [ + ( + "spatial_resolution_x.pdf", + "Spatial resolution X;E_{beam} [GeV];Residual RMS [mm]", + data["res_sci_x_rms"].to_numpy(dtype=float), + data["res_sci_x_rms_err"].to_numpy(dtype=float), + data["res_cer_x_rms"].to_numpy(dtype=float), + data["res_cer_x_rms_err"].to_numpy(dtype=float), + "Residual RMS along X", + ), + ( + "spatial_resolution_y.pdf", + "Spatial resolution Y;E_{beam} [GeV];Residual RMS [mm]", + data["res_sci_y_rms"].to_numpy(dtype=float), + data["res_sci_y_rms_err"].to_numpy(dtype=float), + data["res_cer_y_rms"].to_numpy(dtype=float), + data["res_cer_y_rms_err"].to_numpy(dtype=float), + "Residual RMS along Y", + ), + ] + + for file_name, title, sci_rms, sci_rms_err, cer_rms, cer_rms_err, subtitle in configs: + canvas = setup_canvas(file_name.replace(".pdf", "")) + mg = ROOT.TMultiGraph() + mg.SetTitle(title) + + graph_sci = make_graph(energy, sci_rms, sci_rms_err, ROOT.kRed + 1, 21, f"sci_{file_name}") + graph_cer = make_graph(energy, cer_rms, cer_rms_err, ROOT.kBlue + 1, 22, f"cer_{file_name}") + + if graph_sci is not None: + mg.Add(graph_sci, "PE") + if graph_cer is not None: + mg.Add(graph_cer, "PE") + + mg.Draw("A") + mg.GetXaxis().SetTitleOffset(1.15) + mg.GetYaxis().SetTitleOffset(1.35) + + legend = ROOT.TLegend(0.16, 0.70, 0.34, 0.82) + legend.SetBorderSize(0) + legend.SetTextFont(42) + legend.SetTextSize(0.03) + if graph_sci is not None: + legend.AddEntry(graph_sci, "Sci", "pe") + if graph_cer is not None: + legend.AddEntry(graph_cer, "Cer", "pe") + legend.Draw() + + draw_header(0.16, 0.84, subtitle) + canvas.SaveAs(os.path.join(out_dir, file_name)) + + +def main(): + csv_path = sys.argv[1] if len(sys.argv) > 1 else "hidra_summary.csv" + out_dir = sys.argv[2] if len(sys.argv) > 2 else myOutDir + + os.makedirs(out_dir, exist_ok=True) + + data = load_summary(csv_path) + plot_energy_resolution(data, out_dir) + plot_energy_linearity(data, out_dir) + plot_spatial_resolution(data, out_dir) + + print(f"Read: {csv_path}") + print(f"Wrote plots to: {out_dir}") + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/analysis/tb25_analysis/run_ana.sh b/analysis/tb25_analysis/run_ana.sh new file mode 100755 index 0000000..9f9b5b2 --- /dev/null +++ b/analysis/tb25_analysis/run_ana.sh @@ -0,0 +1,56 @@ +#!/usr/bin/env bash +# file: run_ana.sh +set -euo pipefail + +energies=(10 20 30 40 60 80 100 120) +files=( + "../../build/OutputElectrons_26_04_10/DREMTubesout_Run0.root" + "../../build/OutputElectrons_26_04_10/DREMTubesout_Run1.root" + "../../build/OutputElectrons_26_04_10/DREMTubesout_Run2.root" + "../../build/OutputElectrons_26_04_10/DREMTubesout_Run3.root" + "../../build/OutputElectrons_26_04_10/DREMTubesout_Run4.root" + "../../build/OutputElectrons_26_04_10/DREMTubesout_Run5.root" + "../../build/OutputElectrons_26_04_10/DREMTubesout_Run6.root" + "../../build/OutputElectrons_26_04_10/DREMTubesout_Run7.root" +) + +summary_csv="hidra_summary.csv" +pids=() + +rm -f "${summary_csv}" + +for i in "${!energies[@]}"; do + energy="${energies[$i]}" + input_file="${files[$i]}" + log_file="run_${i}_${energy}GeV.log" + + echo "Starting ${input_file} at ${energy} GeV" + root -l -b -q "HidraTB25Ana.C(${energy}, \"${input_file}\")" >"${log_file}" 2>&1 & + pids+=("$!") +done + +status=0 + +for pid in "${pids[@]}"; do + if ! wait "${pid}"; then + status=1 + fi +done + +echo "Done. Shared summary written to: ${summary_csv}" +exit "${status}" + + + + + + + +# "DREMTubesout_Run0.root" +# "DREMTubesout_Run1.root" +# "DREMTubesout_Run2.root" +# "DREMTubesout_Run3.root" +# "DREMTubesout_Run4.root" +# "DREMTubesout_Run5.root" +# "DREMTubesout_Run6.root" +# "DREMTubesout_Run7.root"