From beb2b1fb475063a96d3112c07c9ff14fb90479d3 Mon Sep 17 00:00:00 2001 From: sudipbhattarai Date: Sun, 31 May 2026 10:20:02 -0400 Subject: [PATCH 1/3] Fixed the merge conflict on position.xml file for stackplanes --- geometry/mollerParallel.gdml | 362 +++++++++++++++++++++++++-- geometry/positions.xml | 61 ++++- macros/showermax/run_stackplanes.mac | 54 ++++ 3 files changed, 457 insertions(+), 20 deletions(-) create mode 100644 macros/showermax/run_stackplanes.mac diff --git a/geometry/mollerParallel.gdml b/geometry/mollerParallel.gdml index 83348920e..14e1c9911 100644 --- a/geometry/mollerParallel.gdml +++ b/geometry/mollerParallel.gdml @@ -103,7 +103,7 @@ - + @@ -539,25 +539,200 @@ - + - + - + - + - + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -1408,22 +1583,175 @@ - - - + + + + + + + + + + + - - - + + + + - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/geometry/positions.xml b/geometry/positions.xml index 63e7267ed..1e9ef6ce6 100644 --- a/geometry/positions.xml +++ b/geometry/positions.xml @@ -74,6 +74,61 @@ - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/macros/showermax/run_stackplanes.mac b/macros/showermax/run_stackplanes.mac new file mode 100644 index 000000000..602297ec8 --- /dev/null +++ b/macros/showermax/run_stackplanes.mac @@ -0,0 +1,54 @@ +# Example file + +# store tracks +#/tracking/storeTrajectory 1 + +# This must be called before initialize +/remoll/geometry/setfile geometry/mollerMother.gdml +/remoll/parallel/setfile geometry/mollerParallel.gdml + +#/remoll/physlist/register QGSP_BERT_HP +/remoll/physlist/parallel/enable + +# This must be explicitly called +/run/initialize + +/remoll/printgeometry true + +/control/execute macros/load_magnetic_fieldmaps.mac + +# Raster and initial angle stuff +/remoll/oldras true +/remoll/rasx 5 mm +/remoll/rasy 5 mm + +/remoll/evgen/set moller +/remoll/evgen/beamPolarization +L +# /remoll/field/equationtype 2 +# /remoll/field/steppertype 2 +# /remoll/field/print + +/remoll/beamene 11 GeV +/remoll/beamcurr 65 microampere + +/remoll/SD/disable_all +# /remoll/SD/enable 28 +# /remoll/SD/detect lowenergyneutral 28 +# /remoll/SD/detect secondaries 28 +# /remoll/SD/detect boundaryhits 28 + +/remoll/SD/enable_range 73001 73028 +/remoll/SD/detect_range secondaries 73001 73028 +/remoll/SD/detect_range boundaryhits 73001 73028 + +# /process/list + +# Specify random number seed +/remoll/seed 123456 + +/remoll/filename remollout.root + +#/tracking/verbose 2 + +#/remoll/printtargetinfo +/run/beamOn 1000 From daadca5224b4b34df3887670f74e1bb01b8c8ef5 Mon Sep 17 00:00:00 2001 From: sudipbhattarai Date: Sun, 31 May 2026 10:58:05 -0400 Subject: [PATCH 2/3] Added and relocated the response paramater files Positive and negative charged particles have their own duplicate copy files. e.g e- and e+ have idential individual files --- .../data/fit_param_xy_e+_ifarm.csv} | 0 .../data/fit_param_xy_e-_ifarm.csv | 18 ++++++++++++++++++ .../data/fit_param_xy_gamma_ifarm.csv | 0 .../data/fit_param_xy_mu+_ifarm.csv} | 0 .../data/fit_param_xy_mu-_ifarm.csv | 14 ++++++++++++++ .../data/fit_param_xy_neutron_ifarm.csv | 0 .../data/fit_param_xy_pi+_ifarm.csv} | 0 .../data/fit_param_xy_pi-_ifarm.csv | 14 ++++++++++++++ .../scripts/plot_showermax_response.cpp | 0 .../scripts/remollToQsim.hh | 0 .../scripts/shower-max_resposne_lookup.hh | 0 11 files changed, 46 insertions(+) rename analysis/showermax/{data/fit_param_xy_e-_ifarm.csv => response_lookup/data/fit_param_xy_e+_ifarm.csv} (100%) create mode 100644 analysis/showermax/response_lookup/data/fit_param_xy_e-_ifarm.csv rename analysis/showermax/{ => response_lookup}/data/fit_param_xy_gamma_ifarm.csv (100%) rename analysis/showermax/{data/fit_param_xy_mu-_ifarm.csv => response_lookup/data/fit_param_xy_mu+_ifarm.csv} (100%) create mode 100644 analysis/showermax/response_lookup/data/fit_param_xy_mu-_ifarm.csv rename analysis/showermax/{ => response_lookup}/data/fit_param_xy_neutron_ifarm.csv (100%) rename analysis/showermax/{data/fit_param_xy_pi-_ifarm.csv => response_lookup/data/fit_param_xy_pi+_ifarm.csv} (100%) create mode 100644 analysis/showermax/response_lookup/data/fit_param_xy_pi-_ifarm.csv rename analysis/showermax/{ => response_lookup}/scripts/plot_showermax_response.cpp (100%) rename analysis/showermax/{ => response_lookup}/scripts/remollToQsim.hh (100%) rename analysis/showermax/{ => response_lookup}/scripts/shower-max_resposne_lookup.hh (100%) diff --git a/analysis/showermax/data/fit_param_xy_e-_ifarm.csv b/analysis/showermax/response_lookup/data/fit_param_xy_e+_ifarm.csv similarity index 100% rename from analysis/showermax/data/fit_param_xy_e-_ifarm.csv rename to analysis/showermax/response_lookup/data/fit_param_xy_e+_ifarm.csv diff --git a/analysis/showermax/response_lookup/data/fit_param_xy_e-_ifarm.csv b/analysis/showermax/response_lookup/data/fit_param_xy_e-_ifarm.csv new file mode 100644 index 000000000..37a91e120 --- /dev/null +++ b/analysis/showermax/response_lookup/data/fit_param_xy_e-_ifarm.csv @@ -0,0 +1,18 @@ +energy,xp0,xp1,yp0,yp1,yp2 +5,0.00186854,-6.27248e-07,0.000318276,3.09224e-07,1.11394e-07 +10,0.0102587,-1.50709e-05,0.0182153,2.54559e-05,-1.01119e-06 +50,0.617914,0.00109035,0.734966,-0.000662327,-2.13036e-05 +100,2.40497,0.00304654,2.52607,-0.000402899,-2.82908e-05 +500,17.2653,0.0150206,17.7501,-0.00199242,-0.000123669 +1000,35.8243,0.0354252,36.7149,-0.00443617,-0.000235146 +2000,71.2615,0.0664176,73.2177,-0.0101103,-0.000481341 +1000,35.8243,0.0354252,36.7149,-0.00443617,-0.000235146 +3000,105.66,0.102495,107.507,-0.00432239,-0.000539539 +4000,138.983,0.129233,141.393,-0.00765265,-0.000714903 +5000,170.902,0.156594,174.967,-0.0165793,-0.00102413 +6000,202.591,0.195464,206.781,-0.0185614,-0.00107666 +7000,233.861,0.214678,238.794,-0.0255901,-0.00132374 +8000,265.736,0.253171,270.522,-0.021348,-0.001311 +9000,295.394,0.2864,301.713,-0.0221357,-0.00161815 +10000,324.714,0.318387,331.45,-0.0336423,-0.00174544 +11000,354.866,0.339425,362.352,-0.00987083,-0.0018857 diff --git a/analysis/showermax/data/fit_param_xy_gamma_ifarm.csv b/analysis/showermax/response_lookup/data/fit_param_xy_gamma_ifarm.csv similarity index 100% rename from analysis/showermax/data/fit_param_xy_gamma_ifarm.csv rename to analysis/showermax/response_lookup/data/fit_param_xy_gamma_ifarm.csv diff --git a/analysis/showermax/data/fit_param_xy_mu-_ifarm.csv b/analysis/showermax/response_lookup/data/fit_param_xy_mu+_ifarm.csv similarity index 100% rename from analysis/showermax/data/fit_param_xy_mu-_ifarm.csv rename to analysis/showermax/response_lookup/data/fit_param_xy_mu+_ifarm.csv diff --git a/analysis/showermax/response_lookup/data/fit_param_xy_mu-_ifarm.csv b/analysis/showermax/response_lookup/data/fit_param_xy_mu-_ifarm.csv new file mode 100644 index 000000000..ef8a33b0f --- /dev/null +++ b/analysis/showermax/response_lookup/data/fit_param_xy_mu-_ifarm.csv @@ -0,0 +1,14 @@ +energy,xp0,xp1,yp0,yp1,yp2 +500,7.61904,0.00705542,7.68358,8.26442e-05,-2.35436e-05 +1000,9.84448,0.00476633,9.85678,-0.000641187,-7.83977e-06 +2000,10.6694,0.00424422,10.6791,0.000134362,-1.42913e-06 +1000,9.84448,0.00476633,9.85678,-0.000641187,-7.83977e-06 +3000,10.902,0.00580823,10.7998,-0.000183595,1.279e-05 +4000,10.9523,0.00511675,10.9356,-0.000334526,-9.12267e-07 +5000,11.1126,0.00453355,11.0124,8.11234e-05,1.29421e-05 +6000,11.0563,0.00719545,11.0066,0.000110291,4.91498e-06 +7000,11.3055,0.00669282,11.1802,0.000123959,1.57476e-05 +8000,11.1125,0.00717593,10.9525,-0.000299558,2.83419e-05 +9000,11.2697,0.00621689,11.1181,-0.000846437,2.12847e-05 +10000,11.2074,0.00588017,11.126,-0.0010563,1.76701e-05 +11000,11.2681,0.00494686,11.1262,0.000572583,2.20093e-05 diff --git a/analysis/showermax/data/fit_param_xy_neutron_ifarm.csv b/analysis/showermax/response_lookup/data/fit_param_xy_neutron_ifarm.csv similarity index 100% rename from analysis/showermax/data/fit_param_xy_neutron_ifarm.csv rename to analysis/showermax/response_lookup/data/fit_param_xy_neutron_ifarm.csv diff --git a/analysis/showermax/data/fit_param_xy_pi-_ifarm.csv b/analysis/showermax/response_lookup/data/fit_param_xy_pi+_ifarm.csv similarity index 100% rename from analysis/showermax/data/fit_param_xy_pi-_ifarm.csv rename to analysis/showermax/response_lookup/data/fit_param_xy_pi+_ifarm.csv diff --git a/analysis/showermax/response_lookup/data/fit_param_xy_pi-_ifarm.csv b/analysis/showermax/response_lookup/data/fit_param_xy_pi-_ifarm.csv new file mode 100644 index 000000000..41d6e7307 --- /dev/null +++ b/analysis/showermax/response_lookup/data/fit_param_xy_pi-_ifarm.csv @@ -0,0 +1,14 @@ +energy,xp0,xp1,yp0,yp1,yp2 +500,5.39623,0.00672054,5.5963,-0.000690162,-4.39051e-05 +1000,8.70975,0.00426233,8.69977,0.000156259,-7.83212e-06 +2000,10.7281,0.00291869,10.8914,0.000375828,-4.35732e-05 +1000,8.70975,0.00426233,8.69977,0.000156259,-7.83212e-06 +3000,12.4285,0.00493266,12.6518,0.000701055,-4.96781e-05 +4000,14.5609,-0.000216201,15.0292,-0.00106373,-0.000120579 +5000,17.187,0.00532813,17.9249,0.000601087,-0.000154178 +6000,21.536,0.00222366,21.9014,-0.00511615,-0.000135389 +7000,23.5107,0.015075,23.7392,-0.00262259,-6.06112e-05 +8000,25.8184,0.015339,26.8405,-0.00330766,-0.000222387 +9000,28.3697,0.0249847,28.6818,0.0040931,-0.000111408 +10000,29.5291,0.0173609,29.3689,-0.00523963,-4.497e-05 +11000,30.8697,0.0487093,31.7276,0.00668593,-0.000143964 diff --git a/analysis/showermax/scripts/plot_showermax_response.cpp b/analysis/showermax/response_lookup/scripts/plot_showermax_response.cpp similarity index 100% rename from analysis/showermax/scripts/plot_showermax_response.cpp rename to analysis/showermax/response_lookup/scripts/plot_showermax_response.cpp diff --git a/analysis/showermax/scripts/remollToQsim.hh b/analysis/showermax/response_lookup/scripts/remollToQsim.hh similarity index 100% rename from analysis/showermax/scripts/remollToQsim.hh rename to analysis/showermax/response_lookup/scripts/remollToQsim.hh diff --git a/analysis/showermax/scripts/shower-max_resposne_lookup.hh b/analysis/showermax/response_lookup/scripts/shower-max_resposne_lookup.hh similarity index 100% rename from analysis/showermax/scripts/shower-max_resposne_lookup.hh rename to analysis/showermax/response_lookup/scripts/shower-max_resposne_lookup.hh From cd04cd67e1ae0912a82518f7b025faf60cc8c319 Mon Sep 17 00:00:00 2001 From: sudipbhattarai Date: Sun, 31 May 2026 10:59:42 -0400 Subject: [PATCH 3/3] Updated lookup table scripts Now the implementation of the response lookup table is more straightforward; just add the header file and use the functions from there. --- .../scripts/plot_showermax_response.cpp | 64 ++- .../response_lookup/scripts/remollToQsim.hh | 3 +- .../scripts/shower-max_resposne_lookup.hh | 404 ++++++++++++++---- 3 files changed, 339 insertions(+), 132 deletions(-) diff --git a/analysis/showermax/response_lookup/scripts/plot_showermax_response.cpp b/analysis/showermax/response_lookup/scripts/plot_showermax_response.cpp index 74590c508..c86c6fea0 100644 --- a/analysis/showermax/response_lookup/scripts/plot_showermax_response.cpp +++ b/analysis/showermax/response_lookup/scripts/plot_showermax_response.cpp @@ -5,6 +5,7 @@ // Purpose: Read the remoll rootfile and use shower-max lookup table to plot the PE response. // ----------------------------------------------------------------------------------- +#include "RtypesCore.h" #include "TChain.h" #include "TH1.h" #include "TROOT.h" @@ -17,6 +18,11 @@ #include "./remollToQsim.hh" #include "./shower-max_resposne_lookup.hh" +#ifdef __REMOLLTYPES__ // for clangd LSP support + #include "../include/remolltypes.hh" +#endif // __REMOLLTYPES__ + + // Namespace list to use using std::cout, std::cerr, std::endl; using std::vector, std::string; @@ -24,37 +30,17 @@ using std::ifstream, std::ofstream; using std::pair, std::map; // Main function -void plot_showermax_response(){ - // Genetate random values for x, y, theta, phi - - // Get the fit parameters - vector> fit_data_electron = retrieve_fit_data("e-"); - vector> fit_data_gamma = retrieve_fit_data("gamma"); - vector> fit_data_mu = retrieve_fit_data("mu-"); - vector> fit_data_pi = retrieve_fit_data("pi-"); - vector> fit_data_neutron = retrieve_fit_data("neutron"); - - // Make a map of pid and fit_data - map>> map_fit_data; - map_fit_data[11] = fit_data_electron; - map_fit_data[-11] = fit_data_electron; - map_fit_data[22] = fit_data_gamma; - map_fit_data[13] = fit_data_mu; - map_fit_data[-13] = fit_data_mu; - map_fit_data[211] = fit_data_pi; - map_fit_data[-211] = fit_data_pi; - map_fit_data[2112] = fit_data_neutron; - - // Test pe response - // cout << "PE: " << get_PE_response(map_fit_data[11], 1234, 0, 0) << endl; +void plot_showermax_response(TString rootFile = "~/moller/softwares/remoll-zone/rootfiles/sm_stack/stack_v27/sm_tqStack_3sector_moller_50k_1001.root") { + + // Testing shower-max response + ShowermaxLookup::getPeResponse(73015, 11, 1000, 0.0, 1080.0); // Remoll file path goes here int goodFileCount = 1; // Number of good (non-corrupted) files - TString rootFile= "~/moller/softwares/remoll-zone/rootfiles/smStack_v23/sm_tqStack_3sector_moller_50k_1001.root"; // Define histograms TH1D* h_hitRate = new TH1D("hitRate", "Rate weighted hit; hit.r [mm]; rate[GHz/65uA]", 100, 1000, 1200); - TH1D* h_hitRateSmPE = new TH1D("hitRateSmPE", "Rate weighted hit; hit.r [mm]; rate [GHz/65uA/PE]", 100, 1000, 1200); + TH1D* h_hitRateSmPE = new TH1D("hitRateSmPE", "PE*Rate weighted hit; hit.r [mm]; rate [PE*GHz/65uA]", 100, 1000, 1200); // Declare TChain TChain* T = new TChain("T"); @@ -74,7 +60,7 @@ void plot_showermax_response(){ T->SetBranchAddress("rate", &fRate); // T->SetBranchAddress("ev", &fEv); - //This loop goes over all the events in the root files + // This loop goes over all the events in the root files for (Long64_t iEvent = 0; iEvent < nEntries; iEvent++){ if (iEvent % (nEntries/10) == 0) cout << "Analyzed " << iEvent << " events." << endl; @@ -92,15 +78,16 @@ void plot_showermax_response(){ rate = fRate/1.0e9/goodFileCount; // convert to GHz/uA //Fill histograms with proper cuts - bool all_cuts = (det == 30 && //det number 30 is SM plane + bool all_cuts = (det >= 73001 && det <= 73028 && //det number 30 is SM plane hitr>1020 && hitr<1180)&& (pid==11 || pid==-11 || pid==22 || pid==13 || pid==-13 || pid==211 || pid==-211 || pid==2112) && //particle selection - energy>2 && //energy cut + energy>10 && //energy cut hitpz>0; //particle coming from upstream (pz>0) if (all_cuts) { - std::pair qsimxy = ConvertRemollToQsim(hitx, hity); - pe = get_PE_response(map_fit_data[pid], energy, qsimxy.first, qsimxy.second); + // Double_t pe = ShowermaxLookup::getPeResponseWithoutLP(pid, energy, hitx, hity); + Double_t pe = ShowermaxLookup::getPeResponse(det, pid, energy, hitx, hity); + std::cout << "pid: " << pid << ", energy: " << energy << ", pe: " << pe << std::endl; h_hitRate->Fill(hitr, rate); h_hitRateSmPE->Fill(hitr,rate*pe); @@ -109,18 +96,19 @@ void plot_showermax_response(){ } // Print the rate in shower-max for script validation - double rateTotal = h_hitRate->Integral(); // in GHz + double rateTotal = h_hitRate->Integral(); // in GH cout << "Accepted rate: " << rateTotal << " GHz" << endl; double cathCurrent = h_hitRateSmPE->Integral()*1e9*1.6e-19*1e9; // in nA - cout << "Cathode current: " << cathCurrent << " nA" << endl; + cout << "Total cathode current: " << cathCurrent << " nA" << endl; // Plot the histograms - TCanvas* c1 = new TCanvas("c1", "c1", 1300, 500); - c1->Divide(2, 1); - c1->cd(1); - h_hitRate->Draw("hist E"); - c1->cd(2); - h_hitRateSmPE->Draw("hist E"); + // TCanvas* c1 = new TCanvas("c1", "c1", 1300, 500); + // c1->Divide(2, 1); + // c1->cd(1); + // h_hitRate->Draw("hist E"); + // c1->cd(2); + // h_hitRateSmPE->Draw("hist E"); } + diff --git a/analysis/showermax/response_lookup/scripts/remollToQsim.hh b/analysis/showermax/response_lookup/scripts/remollToQsim.hh index 08922907a..77da2d5d9 100644 --- a/analysis/showermax/response_lookup/scripts/remollToQsim.hh +++ b/analysis/showermax/response_lookup/scripts/remollToQsim.hh @@ -30,7 +30,7 @@ inline std::pair ConvertRemollToQsim(double x_remoll, double y_r // Adjust // Adjust phi_remoll to be in the range [0, 2pi) - if (phi_remoll < 0) { phi_remoll += 2 * TMath::Pi(); } + if (phi_remoll < 0) phi_remoll += 2 * TMath::Pi(); // Calculate qsim plane angles std::vector qsimAngles = CalculateQsimPlaneAngles(); @@ -67,3 +67,4 @@ inline std::pair ConvertRemollToQsim(double x_remoll, double y_r } #endif // REMOLL_TO_QSIM_HH + diff --git a/analysis/showermax/response_lookup/scripts/shower-max_resposne_lookup.hh b/analysis/showermax/response_lookup/scripts/shower-max_resposne_lookup.hh index e36ec3876..15bd92dd3 100644 --- a/analysis/showermax/response_lookup/scripts/shower-max_resposne_lookup.hh +++ b/analysis/showermax/response_lookup/scripts/shower-max_resposne_lookup.hh @@ -1,144 +1,362 @@ #ifndef SHOWER_MAX_RESPONSE_LOOKUP_H #define SHOWER_MAX_RESPONSE_LOOKUP_H -#include #include -#include +#include +#include #include #include #include +#include +#include +#include +#include + #include "RtypesCore.h" #include "TF2.h" +# include "remollToQsim.hh" -using std::cout, std::cerr, std::endl; -using std::vector, std::string; -using std::ifstream, std::ofstream; -using std::stringstream, std::getline; -using std::pair; +namespace ShowermaxLookup{ -const int nParticles = 5; //(e-, gamma, mu-, pi-, neutron) -inline string particleList[] = {"e-", "gamma", "mu-", "pi-", "neutron"}; -inline const Double_t energyList[] = {5, 10, 50, 100, 500, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000}; -inline std::pair get_energy_bounds(Double_t energy) { - Double_t lower = 0; - Double_t upper = 0; - int energyListSize = sizeof(energyList) / sizeof(energyList[0]); +/* ============================================================ + Particle definitions + ============================================================ */ - for (int i = 0; i < energyListSize; ++i) { - if (energyList[i] <= energy) { - lower = energyList[i]; - } - if (energyList[i] > energy) { - upper = energyList[i]; +constexpr int nParticles = 8; + +// PID convention (you may adjust to match Geant4 if needed) +inline const int particlePID[nParticles] = { + 11, // e- + -11, // e+ + 22, // gamma + 13, // mu- + -13, // mu+ + 211, // pi- + -211, // pi+ + 2112 // neutron +}; + +inline const char* particleName[nParticles] = { + "e-", + "e+", + "gamma", + "mu-", + "mu+", + "pi-", + "pi+", + "neutron" +}; + +/* ============================================================ + Energy grid + ============================================================ */ + +inline const Double_t energyList[] = { + 5, 10, 50, 100, 500, + 1000, 2000, 3000, 4000, + 5000, 6000, 7000, 8000, 9000 +}; + +constexpr int nEnergies = + sizeof(energyList) / sizeof(energyList[0]); + +/* ============================================================ + Utility: PID → particle index + ============================================================ */ + +inline int getParticleIndexFromPid(int pid) +{ + for (int i = 0; i < nParticles; ++i) { + if (particlePID[i] == pid) + return i; + } + throw std::runtime_error("Unsupported particle PID"); +} + +/* ============================================================ + Utility: energy bounds + ============================================================ */ + +inline std::pair +getEnergyBounds(Double_t energy) +{ + Double_t low = energyList[0]; + Double_t high = energyList[nEnergies - 1]; + + for (int i = 0; i < nEnergies - 1; ++i) { + if (energy >= energyList[i] && + energy < energyList[i + 1]) { + low = energyList[i]; + high = energyList[i + 1]; break; } } - - return {lower, upper}; + return {low, high}; } -inline vector> retrieve_fit_data(const string& particleName) { - string csv_file_name = Form("../data/fit_param_xy_%s_ifarm.csv", particleName.c_str()); - ifstream file(csv_file_name); - vector> fit_params; +/* ============================================================ + CSV loader (called once per particle) + ============================================================ */ + +inline std::vector> +load_fit_data(const std::string& particle) +{ + std::string filename = + Form("./data/fit_param_xy_%s_ifarm.csv", + particle.c_str()); + std::ifstream file(filename); if (!file.is_open()) { - cerr << "Error: Unable to open the file: " << csv_file_name << endl; - exit(1); + throw std::runtime_error("Cannot open file: " + filename); } - string line; - getline(file, line); - while (getline(file, line)) { - if (line.empty() || line[0] == '#') { continue; } + std::vector> data; + std::string line; - stringstream ss(line); - string token; - vector row; + // Skip header + std::getline(file, line); - while (getline(ss, token, ',')) { + while (std::getline(file, line)) { + if (line.empty() || line[0] == '#') continue; + + std::stringstream ss(line); + std::string token; + std::vector row; + + while (std::getline(ss, token, ',')) { row.push_back(std::stod(token)); } - fit_params.push_back(row); + data.push_back(row); } - - file.close(); - return fit_params; + return data; } -inline TF2* create_fit_function(Double_t energy, const vector>& fit_data) { - Double_t xWeight = 0.50; - Double_t yWeight = 0.50; - string fit_function = "[0]*([1] + x*[2]) + [3]*([4] + y*[5] + y*y*[6])"; - - TF2* fitFormula = new TF2("fitFormula", fit_function.c_str()); - - for (const auto& data : fit_data) { - if (data[0] == energy) { - Double_t xp0 = data[1]; - Double_t xp1 = data[2]; - Double_t yp0 = data[3]; - Double_t yp1 = data[4]; - Double_t yp2 = data[5]; - Double_t tp0 = data[6]; - - fitFormula->SetParameter(0, xWeight); - fitFormula->SetParameter(1, xp0); - fitFormula->SetParameter(2, xp1); - fitFormula->SetParameter(3, yWeight); - fitFormula->SetParameter(4, yp0); - fitFormula->SetParameter(5, yp1); - fitFormula->SetParameter(6, yp2); +/* ============================================================ + Per-particle TF2 cache + ============================================================ */ + +class ShowerMaxParticleCache +{ +public: + explicit ShowerMaxParticleCache(const std::string& particle) + : fit_data_(load_fit_data(particle)) + {} + + const TF2* getFit(Double_t energy) const + { + auto it = fit_cache_.find(energy); + if (it != fit_cache_.end()) + return it->second.get(); + + return build_fit(energy); + } + +private: + const TF2* build_fit(Double_t energy) const + { + static const char* formula = + "[0]*([1] + x*[2]) + [3]*([4] + y*[5] + y*y*[6])"; + + auto fit = std::make_unique( + Form("shmax_%g", energy), + formula + ); + + constexpr Double_t xWeight = 0.5; + constexpr Double_t yWeight = 0.5; + + for (const auto& row : fit_data_) { + if (row[0] == energy) { + fit->SetParameters( + xWeight, + row[1], row[2], + yWeight, + row[3], row[4], row[5] + ); + break; + } } + + auto* ptr = fit.get(); + fit_cache_[energy] = std::move(fit); + return ptr; } - return fitFormula; + +private: + std::vector> fit_data_; + mutable std::map> fit_cache_; +}; + +/* ============================================================ + Global cache manager (one per particle) + ============================================================ */ + +inline ShowerMaxParticleCache& +getParticleCache(int pid) +{ + static std::map caches; + + auto it = caches.find(pid); + if (it != caches.end()) + return it->second; + + int idx = getParticleIndexFromPid(pid); + caches.emplace( + pid, + ShowerMaxParticleCache(particleName[idx]) + ); + + return caches.at(pid); } -inline Double_t interpolate_fit_values(Double_t energy, Double_t x, Double_t y) { - vector> fit_data = retrieve_fit_data("e-"); +/* ============================================================ + Public API — safe for hit loops + ============================================================ */ + +inline Double_t +getPeResponseWithoutLP( + int pid, + Double_t energy, + Double_t x, + Double_t y +) +{ + std::pair qsimxy = ConvertRemollToQsim(x,y); - pair energy_bounds = get_energy_bounds(energy); - Double_t lower = energy_bounds.first; - Double_t upper = energy_bounds.second; + auto& cache = getParticleCache(pid); + auto [eLow, eHigh] = getEnergyBounds(energy); - TF2* fitLower = create_fit_function(lower, fit_data); - TF2* fitUpper = create_fit_function(upper, fit_data); + const TF2* fLow = cache.getFit(eLow); + const TF2* fHigh = cache.getFit(eHigh); - Double_t valueLower = fitLower->Eval(x, y); - Double_t valueUpper = fitUpper->Eval(x, y); + const Double_t vLow = fLow->Eval(qsimxy.first, qsimxy.second); + const Double_t vHigh = fHigh->Eval(qsimxy.first, qsimxy.second); - Double_t value = valueLower + (valueUpper - valueLower) * (energy - lower) / (upper - lower); + Double_t value = vLow + (vHigh - vLow) + * (energy - eLow) / (eHigh - eLow); + + if (value < 0.0) value = 0.0; return value; } -inline Double_t get_PE_response(vector> fit_data, Double_t energy, Double_t x, Double_t y) { - // Get the energy bounds - pair energy_bounds = get_energy_bounds(energy); - Double_t lowerBoundE = energy_bounds.first; - Double_t upperBoundE = energy_bounds.second; - // cout << "Energy bounds: " << lowerBoundE << " " << upperBoundE << endl; +/*============================================================ +Return long-pass preserve factor based on angular sector +Return values: + open sector : 0.22 with 450 nm LP filter + closed/transition : 0.33 with 425 nm LP filter +============================================================*/ +inline Double_t getLongPassPreserveFactor(Double_t x, Double_t y) { + constexpr double PI = M_PI; + std::cout << "Pi = " << PI << std::endl; - // Fit function - // cout << "Creating fit function for energy: " << lowerBoundE << " MeV" << endl; - TF2* fitFnLower_e = create_fit_function(lowerBoundE, fit_data); + Double_t phi = atan2(y, x); + std::cout << "phi = " << phi * 180 / PI << std::endl; - // cout << "Ceating fit function for energy: " << upperBoundE << " MeV" << endl; - TF2* fitFnUpper = create_fit_function(upperBoundE, fit_data); + // Convert to [0, 2*pi) + if (phi < 0.0) + phi += 2.0 * PI; - // From fit function, get the value of the function at given x and y - Double_t valueLower = fitFnLower_e->Eval(x, y); - Double_t valueUpper = fitFnUpper->Eval(x, y); + // --------------------------------------------------------- + // Sector definitions + // --------------------------------------------------------- + constexpr Int_t nSectors = 28; + constexpr Double_t sectorWidth = 2.0 * PI / nSectors; + constexpr Double_t startPhi = -sectorWidth / 2.0; - // cout << "Evaluated values (upper and lower): " << valueLower << ", " << valueUpper << endl; + // Shift phi relative to sector start + Double_t shiftedPhi = phi + startPhi; - // Interpolate the values - Double_t value = (valueUpper - valueLower) / (upperBoundE - lowerBoundE) * (energy - lowerBoundE) + valueLower; + std::cout << "shiftedPhi = " << shiftedPhi * 180 / PI << std::endl; - // cout << "Estimated PEs: " << value << endl; - return value; + // Wrap into [0, 2*pi) + while (shiftedPhi < 0.0) + shiftedPhi += 2.0 * PI; + + while (shiftedPhi >= 2.0 * PI) + shiftedPhi -= 2.0 * PI; + + // Determine sector index: 0 -> 27 + Int_t sector = static_cast(shiftedPhi / sectorWidth); + + std::cout << "sector = " << sector << std::endl; + + // --------------------------------------------------------- + // Pattern: 0 is closed, 1 is transition, 2 is open and 3 is transition + // --------------------------------------------------------- + Int_t pattern = sector % 4; + + std::cout << "pattern = " << pattern << std::endl; + + // Open sector + if (pattern == 2) + return 0.22; + + // Closed + transition sectors + return 0.33; } +/*============================================================ + * Return long-pass preserve factor based on angular sectors (stackplanes det Num) + * Return values: 0.22 for open sector, 0.33 for closed/transition sector +============================================================*/ +inline Double_t calc_LP_preserve(Int_t detNum) +{ + // Expected detector range: 73001 - 73028 (28 stackplanes) + constexpr Int_t detMin = 73001; + constexpr Int_t detMax = 73028; + + // Validate detector number + if (detNum < detMin || detNum > detMax) { + std::cerr << "Invalid detector number: " << detNum << std::endl; + return -1.0; + } + + // Convert to sector index: 1 -> 28 + Int_t sector = detNum - 73000; + std::cout << "sector = " << sector << std::endl; + + // Sector pattern: + // 1 : closed + // 2 : transition + // 3 : open + // 4 : transition + // repeated every 4 sectors + // + // Open sectors preserve 22% with 450 nm LP filter + // Closed + transition sectors preserve 33% with 425 nm LP filter + + Double_t LP_preserve = (sector % 4 == 3) ? 0.22 : 0.33; + std::cout << "LP_preserve = " << LP_preserve << std::endl; + + return LP_preserve; +} + +/* ============================================================ + Public API — safe for hit loops + ============================================================ */ + +inline Double_t getPeResponse(Int_t detNum, Int_t pid, Double_t energy, Double_t x, Double_t y) { + std::pair qsimxy = ConvertRemollToQsim(x,y); + // Double_t lpFactor = getLongPassPreserveFactor(x, y); + Double_t lpFactor = calc_LP_preserve(detNum); + + auto& cache = getParticleCache(pid); + auto [eLow, eHigh] = getEnergyBounds(energy); + + const TF2* fLow = cache.getFit(eLow); + const TF2* fHigh = cache.getFit(eHigh); + + const Double_t vLow = fLow->Eval(qsimxy.first, qsimxy.second); + const Double_t vHigh = fHigh->Eval(qsimxy.first, qsimxy.second); + + Double_t value = vLow + (vHigh - vLow) * (energy - eLow) / (eHigh - eLow); + Double_t valueLP = value * lpFactor; + + if (valueLP < 0.0) valueLP = 0.0; + return valueLP; +} +}; #endif // SHOWER_MAX_RESPONSE_LOOKUP_H