@@ -11,40 +11,53 @@ ENABLE_LOGGING();
1111DataFitter::DataFitter () {};
1212
1313void DataFitter::sort_and_append (std::vector<int >& inv_vrefs,
14- std::vector<int >& pedestals,
14+ std::vector<int >& pedestals,
15+ std::vector<double >& stds,
1516 int & step) {
1617 static auto the_log_{::pflib::logging::get (" inv_vref_scan:sort" )};
1718 // We ignore first and last elements since they miss derivs
1819 struct DerivPoint {
1920 int i;
2021 double LH;
22+ double LH_std;
2123 double RH;
24+ double RH_std;
2225 };
2326 std::vector<DerivPoint> slope_points;
2427
2528 double flat_threshold = 0.05 * step;
2629 std::vector<double > LH_derivs;
30+ std::vector<double > LH_stds;
2731 std::vector<double > RH_derivs;
2832 for (int i = 1 ; i < inv_vrefs.size ()-1 ; i++) {
2933 double LH = static_cast <double >(pedestals[i] - pedestals[i-1 ]) / (inv_vrefs[i] - inv_vrefs[i-1 ]);
3034 double RH = static_cast <double >(pedestals[i] - pedestals[i+1 ]) / (inv_vrefs[i] - inv_vrefs[i+1 ]);
3135
32- // Threshold check. CMS uses 0.05. Need to consider stepsize
36+ // Threshold check. CMS uses 0.05. This value fits with my analysis as well.
3337 if (std::abs (LH) < flat_threshold || std::abs (RH) < flat_threshold) { // flat regime
3438 nonlinear_.push_back ({inv_vrefs[i], pedestals[i], LH, RH});
35- } else { // we're in a linear regime or there's an outlier
36- slope_points.push_back ({i, LH, RH});
39+ } else { // we're in a linear regime or there's outliers
40+ double LH_err = std::abs (stds[i] - stds[i-1 ]) / (inv_vrefs[i] - inv_vrefs[i-1 ]);
41+ double RH_err = std::abs (stds[i] - stds[i+1 ]) / (inv_vrefs[i] - inv_vrefs[i+1 ]);
42+ slope_points.push_back ({i, LH, LH_err, RH, RH_err});
3743 LH_derivs.push_back (LH);
44+ LH_stds.push_back (err);
3845 RH_derivs.push_back (RH);
3946 }
4047 }
41- // Now we get the slopey region, removing outliers by using the median
48+ // Now we get the slopey region. CMS removes outliers by using the ADC median.
49+ // From my analysis I chose to consider the std of each point, which is huge for outliers
50+ // compared to the points in the linear region. We set a selection of linear points when
51+ // they have a std < median(std).
52+ // We could use both LH and RH derivs to improve selections.
53+
54+
55+ LH_std_median_ = pflib::utility::median (LH_stds);
4256 LH_median_ = pflib::utility::median (LH_derivs);
43- RH_median_ = pflib::utility::median (RH_derivs);
4457
4558 for (const auto & p : slope_points) {
4659 if ((std::abs (p.LH - LH_median_) < 0.3 *std::abs (LH_median_)) &&
47- std::abs (p. RH - RH_median_) < 0.3 * std::abs (RH_median_) ) { // Linear regime.
60+ p. LH_err < LH_std_median_ ) { // Linear regime.
4861 pflib_log (info) << " inv_vref is : " << inv_vrefs[p.i ];
4962 linear_.push_back ({inv_vrefs[p.i ], pedestals[p.i ], p.LH , p.RH });
5063 }
@@ -101,7 +114,9 @@ static void inv_vref_scan_writer(Target* tgt, pflib::ROC& roc, size_t nevents,
101114 1 /* dummy */ );
102115
103116 std::vector<int > pedestals_l0;
117+ std::vector<double > stds_l0;
104118 std::vector<int > pedestals_l1;
119+ std::vector<double > stds_l1;
105120 std::vector<int > inv_vrefs;
106121
107122 int step = 1 ;
@@ -141,15 +156,17 @@ static void inv_vref_scan_writer(Target* tgt, pflib::ROC& roc, size_t nevents,
141156 " Unable to get adc for the cofigured format" );
142157 }
143158 }
144- pedestals_l0.push_back (pflib::utility::median (adcs_l0));
159+ pedestals_l0.push_back (pflib::utility::median (adcs_l0));
160+ stds_l0.push_back (pflib::utility::std (adcs_l0));
145161 pedestals_l1.push_back (pflib::utility::median (adcs_l1));
162+ stds_l1.push_back (pflib::utility::std (adcs_l1));
146163 inv_vrefs.push_back (inv_vref);
147164 }
148165 // sort data and fit
149166 DataFitter fitter_l0;
150167 DataFitter fitter_l1;
151- fitter_l0.sort_and_append (inv_vrefs, pedestals_l0, step);
152- fitter_l1.sort_and_append (inv_vrefs, pedestals_l1, step);
168+ fitter_l0.sort_and_append (inv_vrefs, pedestals_l0, stds_l0, step);
169+ fitter_l1.sort_and_append (inv_vrefs, pedestals_l1, stds_l1, step);
153170 int inv_vref_l0 = fitter_l0.fit (target_adc);
154171 int inv_vref_l1 = fitter_l1.fit (target_adc);
155172
0 commit comments