From c1537c8d7a29094c38556ec742990f561f020e76 Mon Sep 17 00:00:00 2001
From: Nicolas Koehler <nicolas.koehler@cern.ch>
Date: Thu, 3 Sep 2020 14:02:32 +0000
Subject: [PATCH] mdt calibration - add setMaxRadius method and control
histograms
---
.../QuasianalyticLineReconstruction.h | 3 +
.../src/QuasianalyticLineReconstruction.cxx | 13 ++-
.../MdtCalibRt/RtCalibrationCurved.h | 41 +++----
.../MdtCalibRt/src/RtCalibrationCurved.cxx | 108 ++++++++----------
.../MdtStationT0Container.h | 19 +--
5 files changed, 78 insertions(+), 106 deletions(-)
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/QuasianalyticLineReconstruction.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/QuasianalyticLineReconstruction.h
index 30562bcca4df..5b79249f1f5d 100644
--- a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/QuasianalyticLineReconstruction.h
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/QuasianalyticLineReconstruction.h
@@ -119,6 +119,9 @@ public:
void setTimeOut(const double & time_out);
///< set the time-out for the track
///< finding to time_out (in seconds)
+ void setMaxRadius(const double& maxR);
+ ///< set the max innerRadius, default for MDT
+ ///< sMDT 7.1mm, MDT 14.6mm
// methods required by the base class "IMdtSegmentFitter" //
bool fit(MuonCalibSegment & r_segment) const;
///< reconstruction of the track using
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/QuasianalyticLineReconstruction.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/QuasianalyticLineReconstruction.cxx
index ab6263e05212..0828dcb44f58 100644
--- a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/QuasianalyticLineReconstruction.cxx
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/QuasianalyticLineReconstruction.cxx
@@ -80,7 +80,7 @@ void QuasianalyticLineReconstruction::init(const double & r_road_width) {
//:: SET THE MAXIMUM RADIUS ::
//::::::::::::::::::::::::::::
- m_r_max = 15.0;
+ m_r_max = 14.6;
//::::::::::::::::::::::::
//:: SET THE ROAD WIDTH ::
@@ -547,6 +547,11 @@ void QuasianalyticLineReconstruction::setTimeOut(const double & time_out) {
}
+void QuasianalyticLineReconstruction::setMaxRadius(const double& maxR) {
+ if (maxR<7 || maxR>15) throw std::runtime_error(Form("File: %s, Line: %d\nQuasianalyticLineReconstruction::setMaxRadius() - given radius %.3f not supported, neither MDT nor sMDT", __FILE__, __LINE__, maxR));
+ m_r_max = maxR;
+}
+
//*****************************************************************************
//:::::::::::::::::::::::
@@ -865,10 +870,8 @@ bool QuasianalyticLineReconstruction::fit(MuonCalibSegment & r_segment,
m_track_hits[k]->localPosition().z()),
xhat, null, null);
double d(std::abs(m_track.signDistFrom(wire)));
- m_chi2 = m_chi2+
- std::pow(d-m_track_hits[k]->driftRadius(),
- 2)/
- m_track_hits[k]->sigma2DriftRadius();
+ double r(std::abs(m_track_hits[k]->driftRadius()));
+ m_chi2 = m_chi2+std::pow(d-r,2)/m_track_hits[k]->sigma2DriftRadius();
}
r_segment.set(
m_chi2/static_cast<double>(m_nb_track_hits-2),
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationCurved.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationCurved.h
index adc89c73aa37..630df6a9c8c3 100644
--- a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationCurved.h
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationCurved.h
@@ -1,22 +1,10 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
-//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-// 10.05.2008, AUTHOR: OLIVER KORTNER
-// Modified: 01.03.2009 by O. Kortner, parabolic extrapolation added.
-// 15.03.2009 by O. Kortner, smoothing added.
-// 18.03.2009 by O. Kortner, method performParabolicExtrapolation
-// added.
-//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
#ifndef MuonCalib_RtCalibrationCurvedH
#define MuonCalib_RtCalibrationCurvedH
-//:::::::::::::::::::::::::::::::
-//:: CLASS RtCalibrationCurved ::
-//:::::::::::::::::::::::::::::::
-
/// \class RtCalibrationCurved
///
/// This class performs the autocalibration of MDT chambers from the residuals
@@ -28,14 +16,10 @@
///
/// \date 10.05.2008
-//::::::::::::::::::
-//:: HEADER FILES ::
-//::::::::::::::::::
-
-// STL //
#include <vector>
#include <string>
#include <list>
+#include <memory>
// CLHEP //
#include "CLHEP/Matrix/SymMatrix.h"
@@ -252,7 +236,7 @@ private:
IRtRelation * m_rt_new; // r-t as determined by the autocalibration
RtCalibrationOutput * m_output; // class holding the results of the
// autocalibration
- MultilayerRtDifference * m_multilayer_rt_difference;
+ std::unique_ptr<MultilayerRtDifference> m_multilayer_rt_difference;
// curved-segment fitting //
double m_r_max; // maximum value for accepted drift radii
CurvedPatRec *m_tracker; // curved segment finder (used for track fitting)
@@ -290,14 +274,17 @@ private:
// describing the curved line
// control histograms //
- TFile *m_tfile; // ROOT file
- TH1F *m_cut_evolution; // cut evolution histogram
- TH1F *m_nb_segment_hits; // number of hits on the segments
- TH1F *m_pull_initial; // initial pull distribution
- TH1F *m_pull_final; // final pull distribution after convergence
- TH2F *m_residuals_initial; // initial residual distribution
- TH2F *m_residuals_final; // final residual distribution after convergence
-
+ std::unique_ptr<TFile> m_tfile; // ROOT file
+ std::unique_ptr<TH1F> m_cut_evolution; // cut evolution histogram
+ std::unique_ptr<TH1F> m_nb_segment_hits; // number of hits on the segments
+ std::unique_ptr<TH1F> m_pull_initial; // initial pull distribution
+ std::unique_ptr<TH1F> m_pull_final; // final pull distribution after convergence
+ std::unique_ptr<TH2F> m_residuals_initial; // initial residual distribution
+ std::unique_ptr<TH2F> m_residuals_initial_all; // initial residual distribution before convergence
+ std::unique_ptr<TH2F> m_residuals_final; // final residual distribution after convergence
+ std::unique_ptr<TH2F> m_driftTime_initial; // final residual distribution after convergence
+ std::unique_ptr<TH2F> m_driftTime_final; // final residual distribution after convergence
+ std::unique_ptr<TH2F> m_adc_vs_residual_final; // final residual distribution after convergence
// private methods //
void init(const double & rt_accuracy, const unsigned int & func_type,
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtCalibrationCurved.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtCalibrationCurved.cxx
index 5db8b5eeb1a0..009f4ca9c933 100644
--- a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtCalibrationCurved.cxx
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtCalibrationCurved.cxx
@@ -109,16 +109,8 @@ RtCalibrationCurved::~RtCalibrationCurved(void) {
if (m_base_function!=0) {
delete m_base_function;
}
- delete m_multilayer_rt_difference;
- if (m_tfile!=0) {
+ if (m_tfile) {
m_tfile->Write();
- delete m_cut_evolution;
- delete m_nb_segment_hits;
- delete m_pull_initial;
- delete m_pull_final;
- delete m_residuals_initial;
- delete m_residuals_final;
- delete m_tfile;
}
}
@@ -207,32 +199,19 @@ void RtCalibrationCurved::switch_on_control_histograms(const std::string & file_
/////////////////////////////////////////////
m_control_histograms = true;
- m_tfile = new TFile(file_name.c_str(), "RECREATE");
+ m_tfile = std::make_unique<TFile>(file_name.c_str(), "RECREATE");
- m_cut_evolution = new TH1F("m_cut_evolution",
- "CUT EVOLUTION",
- 11, -0.5, 10.5);
+ m_cut_evolution = std::make_unique<TH1F>("m_cut_evolution","CUT EVOLUTION",11, -0.5, 10.5);
+ m_nb_segment_hits = std::make_unique<TH1F>("m_nb_segment_hits","NUMBER OF HITS ON THE REFITTED SEGMENTS",11, -0.5, 10.5);
+ m_pull_initial = std::make_unique<TH1F>("m_pull_initial","INITIAL PULL DISTRIBUTION",200, -5.05, 5.05);
+ m_residuals_initial = std::make_unique<TH2F>("m_residuals_initial","INITIAL OF THE REFITTED SEGMENTS",100, -0.5, 15.0, 300, -1.5, 1.5);
+ m_residuals_initial_all = std::make_unique<TH2F>("m_residuals_initial_all","INITIAL OF THE REFITTED SEGMENTS BEFORE CONVERGENCE",300, -15.0, 15.0, 1000, -5, 5);
+ m_residuals_final = std::make_unique<TH2F>("m_residuals_final","FINAL OF THE REFITTED SEGMENTS",100, -0.5, 15.0, 300, -1.5, 1.5);
+ m_driftTime_final = std::make_unique<TH2F>("m_driftTime_final","FINAL DRIFTTIME OF THE REFITTED SEGMENTS",300, -15.0, 15.0, 300,-15.0, 15.0);
+ m_driftTime_initial = std::make_unique<TH2F>("m_driftTime_initial","FINAL DRIFTTIME OF THE REFITTED SEGMENTS",300, -15.0, 15.0, 1300, -100, 1200);
+ m_adc_vs_residual_final = std::make_unique<TH2F>("m_adc_resi_initial","FINAL ADC VS RESIDUAL OF THE REFITTED SEGMENTS",1350, 0, 1350, 300,-15,15);
- m_nb_segment_hits = new TH1F("m_nb_segment_hits",
- "NUMBER OF HITS ON THE REFITTED SEGMENTS",
- 11, -0.5, 10.5);
-
- m_pull_initial = new TH1F("m_pull_initial",
- "INITIAL PULL DISTRIBUTION",
- 200, -5.05, 5.05);
-// m_pull_final = new TH1F("m_pull_final",
-// "FINAL PULL DISTRIBUTION",
-// 200, -5.05, 5.05);
-
- m_residuals_initial = new TH2F("m_residuals_initial",
- "INITIAL OF THE REFITTED SEGMENTS",
- 100, -0.5, 15.0, 300, -1.5, 1.5);
- m_residuals_final = new TH2F("m_residuals_final",
- "FINAL OF THE REFITTED SEGMENTS",
- 100, -0.5, 15.0, 300, -1.5, 1.5);
-
- delete m_multilayer_rt_difference;
- m_multilayer_rt_difference=new MultilayerRtDifference(10000, m_tfile);
+ m_multilayer_rt_difference=std::make_unique<MultilayerRtDifference>(10000, m_tfile.get());
return;
}
@@ -243,10 +222,9 @@ void RtCalibrationCurved::switch_on_control_histograms(const std::string & file_
//::::::::::::::::::::::::::::::::::::::::::
void RtCalibrationCurved::switch_off_control_histograms(void) {
m_control_histograms = false;
- if (m_tfile!=0) {
+ if (m_tfile) {
m_tfile->Write();
- delete m_multilayer_rt_difference;
- m_multilayer_rt_difference=new MultilayerRtDifference(10000);
+ m_multilayer_rt_difference=std::make_unique<MultilayerRtDifference>(10000);
}
return;
@@ -309,6 +287,18 @@ void RtCalibrationCurved::noSmoothing(void) {
//::::::::::::::::::::::::::::
const IMdtCalibrationOutput * RtCalibrationCurved::analyseSegments(const std::vector<MuonCalibSegment*> & seg) {
const IRtRelation *tmp_rt(0);
+
+/////////////////////
+// Initial RESIDUALS //
+/////////////////////
+ for (unsigned int k=0; k<seg.size(); k++) {
+ for (unsigned int l=0; l<seg[k]->hitsOnTrack(); l++) {
+ double rr = (seg[k]->mdtHOT())[l]->driftRadius();
+ double dd = (seg[k]->mdtHOT())[l]->signedDistanceToTrack();
+ if(m_residuals_final) m_residuals_initial_all->Fill(std::abs(dd), std::abs(rr)-std::abs(dd));
+ m_driftTime_initial->Fill(rr, dd);
+ }
+ }
//////////////////////////
// AUTOCALIBRATION LOOP //
@@ -352,13 +342,17 @@ const IMdtCalibrationOutput * RtCalibrationCurved::analyseSegments(const std::ve
//////////////////////////////////////////////
if (!m_do_smoothing) {
// final residuals //
- double r, d;
+ double r(0);
+ double d(0);
+ double adc(0);
for (unsigned int k=0; k<seg.size(); k++) {
for (unsigned int l=0; l<seg[k]->hitsOnTrack(); l++) {
- r = std::abs((seg[k]->mdtHOT())[l]->driftRadius());
- d = std::abs((seg[k]->mdtHOT())[l]->signedDistanceToTrack());
- if(m_residuals_final != NULL)
- m_residuals_final->Fill(d, r-d, 1.0);
+ adc = (seg[k]->mdtHOT())[l]->adcCount();
+ r = (seg[k]->mdtHOT())[l]->driftRadius();
+ d = (seg[k]->mdtHOT())[l]->signedDistanceToTrack();
+ if(m_residuals_final) m_residuals_final->Fill(std::abs(d), std::abs(r)-std::abs(d));
+ m_driftTime_final->Fill(r, d);
+ m_adc_vs_residual_final->Fill(adc,std::abs(r)-std::abs(d));
}
}
return getResults();
@@ -419,10 +413,9 @@ const IMdtCalibrationOutput * RtCalibrationCurved::analyseSegments(const std::ve
// final residuals //
for (unsigned int k=0; k<seg.size(); k++) {
for (unsigned int l=0; l<seg[k]->hitsOnTrack(); l++) {
- double r = std::abs((seg[k]->mdtHOT())[l]->driftRadius());
- double d = std::abs((seg[k]->mdtHOT())[l]->signedDistanceToTrack());
- if(m_residuals_final != NULL)
- m_residuals_final->Fill(d, r-d, 1.0);
+ double r = (seg[k]->mdtHOT())[l]->driftRadius();
+ double d = (seg[k]->mdtHOT())[l]->signedDistanceToTrack();
+ if(m_residuals_final) m_residuals_final->Fill(d, std::abs(r)-std::abs(d), 1.0);
}
}
return getResults();
@@ -461,10 +454,12 @@ const IMdtCalibrationOutput * RtCalibrationCurved::analyseSegments(const std::ve
/////////////////////
for (unsigned int k=0; k<seg.size(); k++) {
for (unsigned int l=0; l<seg[k]->hitsOnTrack(); l++) {
- double r = std::abs((seg[k]->mdtHOT())[l]->driftRadius());
- double d = std::abs((seg[k]->mdtHOT())[l]->signedDistanceToTrack());
- if(m_residuals_final != NULL)
- m_residuals_final->Fill(d, r-d, 1.0);
+ double adc = (seg[k]->mdtHOT())[l]->adcCount();
+ double r = (seg[k]->mdtHOT())[l]->driftRadius();
+ double d = (seg[k]->mdtHOT())[l]->signedDistanceToTrack();
+ if(m_residuals_final) m_residuals_final->Fill(d, std::abs(r)-std::abs(d));
+ m_driftTime_final->Fill(r, d);
+ m_adc_vs_residual_final->Fill(adc,std::abs(r)-std::abs(d));
}
}
@@ -890,7 +885,7 @@ bool RtCalibrationCurved::analyse(const std::vector<MuonCalibSegment*> & seg) {
// do not change the r-t and the endpoints //
if (((k==0 || x_r[k].x2()<0.5) && m_fix_min) ||
- ((k==nb_points || x_r[k].x2()>14.1) && m_fix_max)) {
+ ((k==nb_points || x_r[k].x2()>m_r_max) && m_fix_max)) {
r_corr = 0.0;
x_r[k].set_error(0.01);
}
@@ -935,7 +930,7 @@ bool RtCalibrationCurved::analyse(const std::vector<MuonCalibSegment*> & seg) {
max_k = rt_param.size()-1;
}
for (unsigned int k=min_k; k<max_k; k++) {
- x = (rt_param[k] - 7.6)/7.6;
+ x = (rt_param[k] - 0.5*m_r_max)/(0.5*m_r_max);
r_corr = m_alpha[0];
for (unsigned int l=1; l<m_order; l++) {
r_corr = r_corr+m_alpha[l]*
@@ -1057,16 +1052,9 @@ void RtCalibrationCurved::init(const double & rt_accuracy,
/////////////////////////////
// RESET PRIVATE VARIABLES //
/////////////////////////////
- m_rt=NULL;
+ m_rt=nullptr;
m_r_max = 15.0*CLHEP::mm;
m_control_histograms = false;
- m_tfile = 0;
- m_cut_evolution = 0;
- m_nb_segment_hits = 0;
- m_pull_initial = 0;
- m_pull_final = 0;
- m_residuals_initial = 0;
- m_residuals_final = 0;
m_nb_segments = 0;
m_nb_segments_used = 0;
m_iteration = 0;
@@ -1081,7 +1069,7 @@ void RtCalibrationCurved::init(const double & rt_accuracy,
m_fix_max = fix_max;
m_max_it = std::abs(max_it);
m_do_multilayer_rt_scale=do_multilayer_rt_scale;
- m_multilayer_rt_difference = new MultilayerRtDifference(10000);
+ m_multilayer_rt_difference = std::make_unique<MultilayerRtDifference>(10000);
m_tracker = new CurvedPatRec;
diff --git a/MuonSpectrometer/MuonCalib/MuonCalibStandAlone/MuonCalibStandAloneBase/MuonCalibStandAloneBase/MdtStationT0Container.h b/MuonSpectrometer/MuonCalib/MuonCalibStandAlone/MuonCalibStandAloneBase/MuonCalibStandAloneBase/MdtStationT0Container.h
index ef5752bcdefa..3bc432c61db0 100644
--- a/MuonSpectrometer/MuonCalib/MuonCalibStandAlone/MuonCalibStandAloneBase/MuonCalibStandAloneBase/MdtStationT0Container.h
+++ b/MuonSpectrometer/MuonCalib/MuonCalibStandAlone/MuonCalibStandAloneBase/MuonCalibStandAloneBase/MdtStationT0Container.h
@@ -1,11 +1,7 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
-//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-// 01.02.2007, AUTHOR: OLIVER KORTNER
-//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
#ifndef MuonCalib_MdtStationT0ContainerH
#define MuonCalib_MdtStationT0ContainerH
@@ -24,11 +20,6 @@
///
/// \date 01.02.2007
-//::::::::::::::::::
-//:: HEADER FILES ::
-//::::::::::::::::::
-
-// STL //
#include <vector>
#include <string>
@@ -39,16 +30,16 @@ class MdtStationT0Container {
public:
// Constructors //
MdtStationT0Container(void) :
-m_t0(2, std::vector< std::vector<double> >(4, std::vector<double>(72, 9e9))),
-m_adc(2, std::vector< std::vector<double> >(4, std::vector<double>(72, 9e9))),
+m_t0(2, std::vector< std::vector<double> >(4, std::vector<double>(78, 9e9))),
+m_adc(2, std::vector< std::vector<double> >(4, std::vector<double>(78, 9e9))),
m_t0_loaded(false)
{
}
///< Default constructor.
MdtStationT0Container(const std::string & file_name):
-m_t0(2, std::vector< std::vector<double> >(4, std::vector<double>(72, 9e9))),
-m_adc(2, std::vector< std::vector<double> >(4, std::vector<double>(72, 9e9))),
+m_t0(2, std::vector< std::vector<double> >(4, std::vector<double>(78, 9e9))),
+m_adc(2, std::vector< std::vector<double> >(4, std::vector<double>(78, 9e9))),
m_t0_loaded(false)
{
readT0File(file_name);
--
GitLab