Skip to content
Snippets Groups Projects
Commit c1537c8d authored by Nicolas Koehler's avatar Nicolas Koehler Committed by Walter Lampl
Browse files

mdt calibration - add setMaxRadius method and control histograms

parent 04263fe9
5 merge requests!58791DataQualityConfigurations: Modify L1Calo config for web display,!46784MuonCondInterface: Enable thread-safety checking.,!46776Updated LArMonitoring config file for WD to match new files produced using MT,!45405updated ART test cron job,!42417Draft: DIRE and VINCIA Base Fragments for Pythia 8.3
......@@ -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
......
......@@ -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),
......
/*
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,
......
......@@ -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;
......
/*
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);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment