ATLAS Offline Software
Public Member Functions | Private Member Functions | Private Attributes | List of all members
MuonCalib::AdaptiveResidualSmoothing Class Reference

#include <AdaptiveResidualSmoothing.h>

Collaboration diagram for MuonCalib::AdaptiveResidualSmoothing:

Public Member Functions

 AdaptiveResidualSmoothing ()
 Default constructor. More...
 
void clear ()
 clear the memory of the class More...
 
void addResidual (const double radius, const double residual)
 add the residual at the given radius More...
 
bool addResidualsFromSegment (MuonCalibSegment &seg, bool curved, double road_width)
 reconstruct the given segment and store the residuals; if curved is true a curved segment fit is performed, otherwise a straight segment is fitted to the drift radii; the user must set the road width used in the pattern recognition; returns true in case of success More...
 
RtRelationLookUp performSmoothing (const IRtRelation &rt_rel, unsigned int nb_entries_per_bin, bool fix_t_min, bool fix_t_max)
 use the stored residuals to improve the given r-t relationship to give smoother residuals; the user has to set the number of entries per radial bin; the user can request that the radii for the minimum and maximum drift time are untouched More...
 
RtRelationLookUp performSmoothing (const IRtRelation &rt_rel, const bool &fix_t_min, const bool &fix_t_max)
 use the stored residuals to improve the given r-t relationship to give smoother residuals; the user can request that the radii for the minimum and maximum drift time are untouched More...
 

Private Member Functions

double t_from_r (const IRtRelation &rt_rel, const double r)
 

Private Attributes

std::vector< DataPoint > m_residual_point
 
StraightPatRec m_sfitter
 
CurvedPatRec m_cfitter
 

Detailed Description

Definition at line 48 of file AdaptiveResidualSmoothing.h.

Constructor & Destructor Documentation

◆ AdaptiveResidualSmoothing()

AdaptiveResidualSmoothing::AdaptiveResidualSmoothing ( )

Default constructor.

Definition at line 25 of file AdaptiveResidualSmoothing.cxx.

25 {}

Member Function Documentation

◆ addResidual()

void AdaptiveResidualSmoothing::addResidual ( const double  radius,
const double  residual 
)

add the residual at the given radius

Definition at line 27 of file AdaptiveResidualSmoothing.cxx.

27  {
28  // make the data point //
29  CLHEP::HepVector point(2);
30  point[0] = radius;
31  point[1] = residual;
32  DataPoint residual_point(point, 0);
33 
34  // store the data point //
35  m_residual_point.push_back(residual_point);
36 
37  return;
38 }

◆ addResidualsFromSegment()

bool AdaptiveResidualSmoothing::addResidualsFromSegment ( MuonCalibSegment &  seg,
bool  curved,
double  road_width 
)

reconstruct the given segment and store the residuals; if curved is true a curved segment fit is performed, otherwise a straight segment is fitted to the drift radii; the user must set the road width used in the pattern recognition; returns true in case of success

Definition at line 39 of file AdaptiveResidualSmoothing.cxx.

39  {
41  // VARIABLES //
43 
44  CLHEP::HepVector point(2); // auxiliary point
45  IMdtPatRecFitter *fitter = nullptr;
46 
48  // RESIDUAL DETERMINATION AND STORAGE //
50 
51  // perform fit //
52  if (curved) {
53  // set road width /
54  m_cfitter.setRoadWidth(road_width);
55 
56  // set fitter pointer //
57  fitter = &m_cfitter;
58 
59  } else {
60  // set road width /
61  m_sfitter.setRoadWidth(road_width);
62 
63  // set fitter pointer //
64  fitter = &m_sfitter;
65  }
66 
67  fitter->SetRefineSegmentFlag(false);
68 
69  CurvedLine curved_track;
70  MTStraightLine straight_track;
71  IMdtSegmentFitter::HitSelection selection(seg.mdtHitsOnTrack(), 0);
72 
73  // perform the track fit; return in case of failure //
74 
75  // reject bad fits //
76  if (curved) {
77  if (!m_cfitter.fit(seg, selection, curved_track)) return false;
78  if (curved_track.chi2PerDegreesOfFreedom() > 50) { return false; }
79  } else {
80  if (!m_sfitter.fit(seg, selection, straight_track)) return false;
81  if (straight_track.chi2PerDegreesOfFreedom() > 50) { return false; }
82  }
83 
84  // store the residuals //
85  const MuonCalibSegment::MdtHitVec &trackHits = curved ? curved_track.trackHits() : straight_track.trackHits();
86  for (const auto & trackHit : trackHits) {
87  point[0] = trackHit->driftRadius();
88  point[1] = trackHit->driftRadius() - std::abs(trackHit->signedDistanceToTrack());
89  m_residual_point.emplace_back(point, 0);
90  }
91 
92  return true;
93 }

◆ clear()

void AdaptiveResidualSmoothing::clear ( )

clear the memory of the class

Definition at line 26 of file AdaptiveResidualSmoothing.cxx.

26 { m_residual_point.clear(); }

◆ performSmoothing() [1/2]

RtRelationLookUp AdaptiveResidualSmoothing::performSmoothing ( const IRtRelation &  rt_rel,
const bool &  fix_t_min,
const bool &  fix_t_max 
)

use the stored residuals to improve the given r-t relationship to give smoother residuals; the user can request that the radii for the minimum and maximum drift time are untouched

Definition at line 178 of file AdaptiveResidualSmoothing.cxx.

178  {
180  // VARIABLES //
182 
183  unsigned int nb_bins; // number of bins in r
184  unsigned int min_nb_entries_per_bin(20); // minimum number of
185  // entries per bin
186  double step_size; // real step size [mm]
187  double r_max(16.0); // maximum drift radius
188  double aux_rad; // auxiliary radius
189  double aux_corr; // auxiliary correction value
190  unsigned int bin_index;
191 
193  // DETERMINE THE NUMBER OF BINS //
195 
196  // take the sqrt from the number of residuals to determine the binsize.
197  // for low stat, a larger statistical error on the bins is allowed
198  // for example:
199  // 1000 segments -> 7.6% and 35 bins
200  // 5000 segments -> 5.1% and 64 bins
201  // 10000 segments-> 4.1% and 91 bins
202  nb_bins = static_cast<unsigned int>(std::sqrt(m_residual_point.size() / 6.));
203 
204  // calculate min_nb_per_bin
205  if (nb_bins == 0) { nb_bins = 1; }
206  min_nb_entries_per_bin = m_residual_point.size() / nb_bins;
207 
209  // RETURN IF THERE ARE NO RESIDUALS //
211 
212  if (m_residual_point.size() == 0) {
213  MsgStream log(Athena::getMessageSvc(), "AdaptiveResidualSmoothing");
214  log << MSG::WARNING << "performSmoothing() - No residuals are stored. no correction applied to r(t)!" << endmsg;
215  }
216 
217  step_size = r_max / static_cast<double>(nb_bins);
218 
220  // DETERMINE THE AVERAGE RESIDUAL VALUES IN THE RADIAL BINS //
222 
223  // vector for the correction function //
224  std::vector<SamplePoint> corr(nb_bins);
225  std::vector<double> radii(nb_bins);
226 
227  // sort the residuals point in the residual value for a simple outlyer //
228  // rejection performed later //
229  for (auto & k : m_residual_point) { k.setReferenceComponent(1); }
230  sort(m_residual_point.begin(), m_residual_point.end());
231 
232  // auxiliary data point arrays //
233  std::vector<std::vector<const DataPoint *> > sample_points_per_r_bin(nb_bins);
234 
235  // group data points in r bins //
236  for (auto & k : m_residual_point) {
237  if (std::abs(k.dataVector()[0]) >= r_max) { continue; }
238  bin_index = static_cast<unsigned int>(std::abs(k.dataVector()[0]) / step_size);
239  sample_points_per_r_bin[bin_index].push_back(&k);
240  }
241 
242  // loop over the r bins and determine the r-t corrections //
243  // static ofstream outfile("tmp.txt");
244  // outfile << endl;
245  for (unsigned int bin = 0; bin < nb_bins; bin++) {
246  aux_rad = 0.0;
247  aux_corr = 0.0;
248 
249  // remove the first and last 1% of the residuals as simple outlyer removal //
250  unsigned int k_min(static_cast<unsigned int>(0.01 * sample_points_per_r_bin[bin].size()));
251  unsigned int k_max(static_cast<unsigned int>(0.99 * sample_points_per_r_bin[bin].size()));
252  // outfile << bin << "\t" << (0.5+bin)*step_size << "\t"
253  // << k_min << "\t" << k_max << endl;
254 
255  // loop over residuals in a bin
256  for (unsigned int k = k_min; k < k_max; k++) {
257  aux_rad = aux_rad + sample_points_per_r_bin[bin][k]->dataVector()[0];
258  aux_corr = aux_corr + sample_points_per_r_bin[bin][k]->dataVector()[1];
259  }
260  if (k_max - k_min < 0.25 * min_nb_entries_per_bin) {
261  radii[bin] = (0.5 + bin) * step_size;
262  corr[bin] = SamplePoint(radii[bin], 0.0, 1.0);
263  } else {
264  radii[bin] = aux_rad / static_cast<double>(k_max - k_min);
265  corr[bin] = SamplePoint(radii[bin], aux_corr / static_cast<double>(k_max - k_min), 1.0);
266  }
267  }
268 
269  // correction polygon //
270  PolygonBase polygon(radii);
271  BaseFunctionFitter fitter(nb_bins);
272  fitter.fit_parameters(corr, 1, corr.size(), polygon);
273 
274  // create output r-t relationship //
275  std::vector<double> rt_params;
276  rt_params.push_back(rt_rel.tLower());
277  rt_params.push_back(0.01 * (rt_rel.tUpper() - rt_rel.tLower()));
278  // ofstream outfile("out2.txt");
279  for (double t = rt_rel.tLower(); t <= rt_rel.tUpper(); t = t + rt_params[1]) {
280  double delta_r(0.0);
281  for (unsigned int l = 0; l < radii.size(); l++) {
282  delta_r = delta_r + fitter.coefficients()[l] * polygon.value(l, rt_rel.radius(t));
283  }
284  if (rt_rel.radius(t) > r_max) { delta_r = 0.0; }
285  if (fix_t_min && (rt_rel.radius(t) < 0.5)) { delta_r = 0.0; }
286  if (fix_t_max && (rt_rel.radius(t) > 14.1)) { delta_r = 0.0; }
287  rt_params.push_back(rt_rel.radius(t) - delta_r);
288  }
289 
290  RtRelationLookUp improved_rt(rt_params);
291  // for (double t=rt_rel.tLower(); t<=rt_rel.tUpper(); t=t+rt_params[1]) {
292  // outfile << rt_rel.radius(t) << "\t" << improved_rt.radius(t)
293  // << "\t" << endl;
294  // }
295 
296  return improved_rt;
297 }

◆ performSmoothing() [2/2]

RtRelationLookUp AdaptiveResidualSmoothing::performSmoothing ( const IRtRelation &  rt_rel,
unsigned int  nb_entries_per_bin,
bool  fix_t_min,
bool  fix_t_max 
)

use the stored residuals to improve the given r-t relationship to give smoother residuals; the user has to set the number of entries per radial bin; the user can request that the radii for the minimum and maximum drift time are untouched

Definition at line 95 of file AdaptiveResidualSmoothing.cxx.

96  {
98  // VARIABLES //
100 
101  ConstantContentBinMaker bin_maker(m_residual_point, 0.001);
102  std::vector<unsigned int> ref_coord(1);
103 
105  // CHECKS //
107 
108  if (m_residual_point.size() == 0) {
109  throw std::runtime_error(
110  Form("File: %s, Line: %d\nAdaptiveResidualSmoothing::performSmoothing - No residuals stored.", __FILE__, __LINE__));
111  }
112 
113  if (m_residual_point.size() < nb_entries_per_bin) {
114  throw std::runtime_error(
115  Form("File: %s, Line: %d\nAdaptiveResidualSmoothing::performSmoothing - Not enough residuals stored.", __FILE__, __LINE__));
116  }
117 
119  // PERFORM BINNING //
121 
122  ref_coord[0] = 0;
123  bin_maker.binDataPoints(nb_entries_per_bin, ref_coord);
124 
126  // DETERMINE r-t COORECTION //
128 
129  // get a polygon through the bin centres //
130  std::vector<double> rad(bin_maker.getBins().size());
131  std::vector<SamplePoint> corr(bin_maker.getBins().size());
132  for (unsigned int k = 0; k < bin_maker.getBins().size(); k++) {
133  rad[k] = (bin_maker.getBins()[k])->centreOfGravity()[0];
134  corr[k].set_x1(rad[k]);
135  corr[k].set_x2((bin_maker.getBins()[k])->centreOfGravity()[1]);
136  corr[k].set_error(1.0);
137  }
138  std::sort(rad.begin(), rad.end());
139 
140  std::vector<double> radi;
141  radi.push_back(rad[0]);
142  unsigned int counter(0);
143  for (unsigned int k = 1; k < rad.size(); k++) {
144  if (rad[counter] < rad[k]) {
145  radi.push_back(rad[k]);
146  counter++;
147  }
148  }
149 
150  PolygonBase polygon(radi);
151  BaseFunctionFitter fitter(bin_maker.getBins().size());
152  fitter.fit_parameters(corr, 1, corr.size(), polygon);
153 
154  // create an improved r-t relationship //
155  std::vector<double> rt_params;
156  rt_params.push_back(rt_rel.tLower());
157  rt_params.push_back(0.01 * (rt_rel.tUpper() - rt_rel.tLower()));
158  // ofstream outfile("out2.txt");
159  for (double t = rt_rel.tLower(); t <= rt_rel.tUpper(); t = t + rt_params[1]) {
160  double delta_r(0.0);
161  for (unsigned int l = 0; l < radi.size(); l++) {
162  delta_r = delta_r + fitter.coefficients()[l] * polygon.value(l, rt_rel.radius(t));
163  }
164  if (fix_t_min && (rt_rel.radius(t) < 0.5)) { delta_r = 0.0; }
165  if (fix_t_max && (rt_rel.radius(t) > 14.1)) { delta_r = 0.0; }
166  rt_params.push_back(rt_rel.radius(t) - delta_r);
167  }
168 
169  RtRelationLookUp improved_rt(rt_params);
170  // for (double t=rt_rel.tLower(); t<=rt_rel.tUpper(); t=t+rt_params[1]) {
171  // outfile << rt_rel.radius(t) << "\t" << improved_rt.radius(t)
172  // << "\t" << endl;
173  // }
174 
175  return improved_rt;
176 }

◆ t_from_r()

double AdaptiveResidualSmoothing::t_from_r ( const IRtRelation &  rt_rel,
const double  r 
)
private

Definition at line 298 of file AdaptiveResidualSmoothing.cxx.

298  {
300  // VARIABLES //
302 
303  double precision(0.001); // spatial precision of the inversion
304  double t_max(rt_rel.tUpper()); // upper time search limit
305  double t_min(rt_rel.tLower()); // lower time search limit
306 
308  // SEARCH FOR THE CORRESPONDING DRIFT TIME //
310 
311  while (t_max - t_min > 0.1 && std::abs(rt_rel.radius(0.5 * (t_min + t_max)) - r) > precision) {
312  if (rt_rel.radius(0.5 * (t_min + t_max)) > r) {
313  t_max = 0.5 * (t_min + t_max);
314  } else {
315  t_min = 0.5 * (t_min + t_max);
316  }
317  }
318 
319  return 0.5 * (t_min + t_max);
320 }

Member Data Documentation

◆ m_cfitter

CurvedPatRec MuonCalib::AdaptiveResidualSmoothing::m_cfitter
private

Definition at line 86 of file AdaptiveResidualSmoothing.h.

◆ m_residual_point

std::vector<DataPoint> MuonCalib::AdaptiveResidualSmoothing::m_residual_point
private

Definition at line 84 of file AdaptiveResidualSmoothing.h.

◆ m_sfitter

StraightPatRec MuonCalib::AdaptiveResidualSmoothing::m_sfitter
private

Definition at line 85 of file AdaptiveResidualSmoothing.h.

The documentation for this class was generated from the following files:
beamspotman.r
def r
Definition: beamspotman.py:676
LArSamples::FitterData::fitter
const ShapeFitter * fitter
Definition: ShapeFitter.cxx:23
MuonCalib::StraightPatRec::setRoadWidth
void setRoadWidth(const double r_road_width)
set the road width for the pattern recognition = r_road_width
Definition: StraightPatRec.cxx:279
ClusterSeg::residual
@ residual
Definition: ClusterNtuple.h:20
MuonCalib::ConstantContentBinMaker
Definition: ConstantContentBinMaker.h:37
MuonCalib::CurvedLine
Definition: CurvedLine.h:30
MuonCalib::MuonCalibSegment::MdtHitVec
std::vector< MdtHitPtr > MdtHitVec
Definition: MuonCalibSegment.h:45
bin
Definition: BinsDiffFromStripMedian.h:43
MuonCalib::CurvedPatRec::fit
bool fit(MuonCalibSegment &r_segment) const
reconstruction of the track using all hits in the segment "r_segment", returns true in case of succes...
Definition: CurvedPatRec.cxx:23
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:158
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
MuonCalib::RtRelationLookUp
Equidistant look up table for rt-relations with the time as key.
Definition: RtRelationLookUp.h:23
Athena::getMessageSvc
IMessageSvc * getMessageSvc(bool quiet=false)
Definition: getMessageSvc.cxx:20
std::sort
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
Definition: DVL_algorithms.h:554
MuonCalib::AdaptiveResidualSmoothing::m_cfitter
CurvedPatRec m_cfitter
Definition: AdaptiveResidualSmoothing.h:86
MuonCalib::MTStraightLine::chi2PerDegreesOfFreedom
double chi2PerDegreesOfFreedom() const
Return chi2 / number of TrackHits - 2.
Definition: MTStraightLine.cxx:139
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
MuonCalib::BaseFunctionFitter
Definition: BaseFunctionFitter.h:39
MuonCalib::CurvedPatRec::setRoadWidth
void setRoadWidth(const double r_road_width)
set the road width [mm] for the pattern recognition = r_road_width
Definition: CurvedPatRec.cxx:21
MuonCalib::AdaptiveResidualSmoothing::m_residual_point
std::vector< DataPoint > m_residual_point
Definition: AdaptiveResidualSmoothing.h:84
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
MuonCalib::CurvedLine::trackHits
const MdtHitVec & trackHits() const
Definition: CurvedLine.cxx:189
MuonCalib::DataPoint
Definition: DataPoint.h:32
MuonCalib::IRtRelation::tUpper
virtual double tUpper() const =0
Returns the upper time covered by the r-t.
MuonCalib::MuonCalibSegment::mdtHitsOnTrack
unsigned int mdtHitsOnTrack() const
retrieve the number of MdtCalibHitBase s assigned to this segment
Definition: MuonCalibSegment.cxx:146
selection
const std::string selection
Definition: fbtTestBasics.cxx:74
MuonCalib::IRtRelation::radius
virtual double radius(double t) const =0
returns drift radius for a given time
MuonCalib::StraightPatRec::fit
bool fit(MuonCalibSegment &r_segment) const
reconstruction of the track using all hits in the segment "r_segment", returns true in case of succes...
Definition: StraightPatRec.cxx:289
MuonCalib::SamplePoint
Definition: SamplePoint.h:16
plotBeamSpotVxVal.bin
int bin
Definition: plotBeamSpotVxVal.py:83
ParticleGun_SamplingFraction.radius
radius
Definition: ParticleGun_SamplingFraction.py:96
MuonCalib::CurvedLine::chi2PerDegreesOfFreedom
double chi2PerDegreesOfFreedom() const
Return chi2 / number of TrackHits - 3.
Definition: CurvedLine.cxx:186
MuonCalib::IRtRelation::tLower
virtual double tLower() const =0
Returns the lower time covered by the r-t.
std::sort
void sort(typename std::reverse_iterator< DataModel_detail::iterator< DVL > > beg, typename std::reverse_iterator< DataModel_detail::iterator< DVL > > end, const Compare &comp)
Specialization of sort for DataVector/List.
Definition: DVL_algorithms.h:623
MuonCalib::IMdtPatRecFitter
Definition: IMdtPatRecFitter.h:18
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
MuonCalib::MTStraightLine
Definition: MTStraightLine.h:16
MuonCalib::MTStraightLine::trackHits
const MdtHitVec & trackHits() const
Definition: MTStraightLine.cxx:142
MuonCalib::PolygonBase
Definition: PolygonBase.h:42
test_pyathena.counter
counter
Definition: test_pyathena.py:15
python.SystemOfUnits.rad
int rad
Definition: SystemOfUnits.py:111
MuonCalib::AdaptiveResidualSmoothing::m_sfitter
StraightPatRec m_sfitter
Definition: AdaptiveResidualSmoothing.h:85
fitman.k
k
Definition: fitman.py:528
MuonCalib::IMdtSegmentFitter::HitSelection
std::vector< unsigned int > HitSelection
Definition: IMdtSegmentFitter.h:32
Generated on Fri Jan 10 2025 21:58:00 for ATLAS Offline Software by   doxygen 1.8.18