nhcal_basic_distribution/scripts/NeutronThresholdUtil.h

0001 // NeutronThresholdsUtil.h
0002 #ifndef NEUTRON_THRESHOLDS_UTIL_H
0003 #define NEUTRON_THRESHOLDS_UTIL_H
0004
0005 #include <vector>
0006 #include <iostream>
0007 #include <TH2D.h>
0008
0009
0010 // Constants for the neutron threshold analysis
0011 namespace NeutronThresholds {
0012     const double E_MIP = 0.00075;
0013
0014     // Time thresholds in ns
0015     const std::vector<double> T_MAX = {25, 50, 100, 200, 500};
0016
0017     // Energy thresholds based on MIP
0018     const std::vector<double> E_TH = {0.5 * E_MIP, 0.25 * E_MIP, 0.1 * E_MIP, 0.05 * E_MIP, 0};
0019
0020     // Initialize the hits passed counters
0021     inline std::vector<std::vector<double>> createHitsPassedMatrix() {
0022         return std::vector<std::vector<double>>(5, std::vector<double>(5, 0));
0023     }
0024
0025     // Process the contributions for a hit
0026     template<typename ContribCollection>
0027     void processContributions(const ContribCollection& contrib, std::vector<std::vector<double>>& hits_passed, std::vector<std::vector<double>>& hits_passed_telap,
0028         TH1* h_contrib_time = nullptr, TH1* h_contrib_energy = nullptr, double E_TH_correction = 1.0, bool debug = false
0029     ) {
0030         std::vector<double> E_sum(5, 0);
0031         std::vector<double> E_sum_telap(5, 0);
0032         double t_min = 1100;
0033
0034         for (unsigned c = 0; c < contrib.size(); ++c) {
0035             if (contrib.at(c).getTime() < t_min) { t_min = contrib.at(c).getTime(); }
0036         }
0037
0038         if(debug) std::cout << contrib.at(0).getTime() << std::endl;
0039
0040         for (unsigned c = 0; c < contrib.size(); ++c) {
0041             if(debug) std::cout << "hit time = " << contrib.at(c).getTime() << std::endl;
0042
0043             // double t_min_per_contrib = std::min(1100, contrib.at(c).getTime());
0044
0045             if (h_contrib_time) h_contrib_time->Fill(contrib.at(c).getTime());
0046             if (h_contrib_energy) h_contrib_energy->Fill(contrib.at(c).getEnergy());
0047
0048             for (int itm = 0; itm < T_MAX.size(); itm++) {
0049                 if (contrib.at(c).getTime() < T_MAX[itm]) {
0050                     E_sum[itm] += contrib.at(c).getEnergy();
0051                 }
0052                 if ((contrib.at(c).getTime() - t_min) < T_MAX[itm]) {
0053                     E_sum_telap[itm] += contrib.at(c).getEnergy();
0054                 }
0055             }
0056         }
0057
0058         for (int iet = 0; iet < E_TH.size(); iet++) {
0059             for (int ies = 0; ies < E_sum.size(); ies++) {
0060                 if (E_sum[ies] > E_TH[iet]*E_TH_correction) {
0061                     hits_passed[iet][ies] += 1;
0062                 }
0063                 if (E_sum_telap[ies] > E_TH[iet]*E_TH_correction) {
0064                     hits_passed_telap[iet][ies] += 1;
0065                 }
0066             }
0067         }
0068     }
0069
0070     // Fill histograms with the hits_passed data
0071     void fillThresholdHistograms(
0072         const std::vector<std::vector<double>>& hits_passed,
0073         const std::vector<std::vector<double>>& hits_passed_telap,
0074         TH2* h_energy_vs_time,
0075         TH2* h_energy_vs_telap,
0076         TH2* h_energy_vs_time_total,
0077         double E_TH_correction = 1.0
0078     ) {
0079         for (int iet = 0; iet < E_TH.size(); iet++) {
0080             for (int itm = 0; itm < T_MAX.size(); itm++) {
0081                 if (hits_passed[iet][itm] > 0) {
0082                     h_energy_vs_time->Fill(E_TH[iet]*E_TH_correction, T_MAX[itm]);
0083                 }
0084                 if (hits_passed_telap[iet][itm] > 0) {
0085                     h_energy_vs_telap->Fill(E_TH[iet]*E_TH_correction, T_MAX[itm]);
0086                 }
0087                 h_energy_vs_time_total->Fill(E_TH[iet]*E_TH_correction, T_MAX[itm]);
0088             }
0089         }
0090     }
0091 }
0092
0093 #endif // NEUTRON_THRESHOLDS_UTIL_H