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 (void)
 Default constructor. More...
 
void clear (void)
 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 ( void  )

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 ( void  )

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::IRtRelation::tUpper
virtual double tUpper(void) const =0
ClusterSeg::residual
@ residual
Definition: ClusterNtuple.h:20
MuonCalib::ConstantContentBinMaker
Definition: ConstantContentBinMaker.h:37
MuonCalib::CurvedLine
Definition: CurvedLine.h:31
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:24
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::StraightPatRec::setRoadWidth
void setRoadWidth(const double &r_road_width)
set the road width for the pattern recognition = r_road_width
Definition: StraightPatRec.cxx:279
MuonCalib::BaseFunctionFitter
Definition: BaseFunctionFitter.h:47
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:198
MuonCalib::IRtRelation::tLower
virtual double tLower(void) const =0
MuonCalib::DataPoint
Definition: DataPoint.h:32
selection
std::string selection
Definition: fbtTestBasics.cxx:73
MuonCalib::MuonCalibSegment::mdtHitsOnTrack
unsigned int mdtHitsOnTrack() const
retrieve the number of MdtCalibHitBase s assigned to this segment
Definition: MuonCalibSegment.cxx:147
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::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:17
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:195
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:19
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 Tue Jul 2 2024 22:22:56 for ATLAS Offline Software by   doxygen 1.8.18